水下航行體充氣上浮仿真方法研究

張曉光， 李斌， 黨會學(xué)， 溫金鵬， 孫潘

(1.西北工業(yè)大學(xué) 航空學(xué)院， 陜西 西安 710072； 2.長安大學(xué) 建筑工程學(xué)院， 陜西 西安 710061；3.中國工程物理研究院 總體工程研究所， 四川 綿陽 621900)

0 引言

水下無人航行體是指能在水下自主航行的智能化裝置。航行體的相關(guān)裝置成本很高，在科研試驗(yàn)或日常操演中，為了使航行體在航行結(jié)束后可靠上浮，以便安全打撈、回收，必須設(shè)計可靠的助浮裝置。充氣式上浮裝置可以提供較大正浮力，且氣囊折疊后所占空間小，質(zhì)量輕，調(diào)試、使用、維修方便，可重復(fù)使用，是一種理想的助浮裝置。

助浮裝置對航行體運(yùn)動軌跡有較大影響，為給予助浮裝置的研制提供參考，需進(jìn)行航行體6自由度(6DOF)運(yùn)動的仿真模擬。對水下航行體充氣上浮的研究，目前多采用試驗(yàn)方法。呂汝信[1]曾對氣囊式充氣上浮裝置進(jìn)行過系統(tǒng)研究，包括裝置的系統(tǒng)構(gòu)成、導(dǎo)流罩解脫技術(shù)等，并進(jìn)行了水中拖曳試驗(yàn)、導(dǎo)流罩干擾試驗(yàn)以及實(shí)航試驗(yàn)。甄文強(qiáng)等[2]利用試驗(yàn)設(shè)計(DOE)方法研究了無人航行器的充氣上浮過程，重點(diǎn)研究了氣囊充氣過程對潛航器運(yùn)動的影響。李春雨等[3]研究了浮囊整流罩的安全解脫設(shè)計及試驗(yàn)驗(yàn)證方法。但試驗(yàn)研究受試驗(yàn)投入、場地和測試設(shè)備等條件的影響，成本相對較高，并且水下試驗(yàn)過程的可測數(shù)據(jù)往往非常有限，難以滿足精細(xì)分析的需要。隨著計算流體動力學(xué)以及動網(wǎng)格技術(shù)的發(fā)展，水下航行體充氣上浮動力學(xué)仿真的相關(guān)文獻(xiàn)陸續(xù)發(fā)表。王棟等[4]通過有限元方法，建立三維造波水池模型，研究了浮囊式操雷在不同強(qiáng)度波浪作用下的漂浮過程，但不能模擬展開上浮過程。程文鑫等[5]、葉慧娟等[6]建立了魚雷用助浮環(huán)型氣囊的充氣展開動力學(xué)模型，研究了不同水深條件下氣囊單純充氣展開行為，但未研究氣囊體積膨脹對魚雷航行姿態(tài)的影響。

要同時進(jìn)行水下航行體氣囊展開和上浮過程的緊耦合行為分析，需要將氣囊展開動力學(xué)模型和航行體的水動力學(xué)模型同步進(jìn)行緊耦合求解。實(shí)現(xiàn)這一緊耦合求解，一方面要模擬復(fù)雜的折疊氣囊在充入氣流作用下的非規(guī)則展開過程，還需應(yīng)用動網(wǎng)格技術(shù)模擬水下航行體的6DOF運(yùn)動，以實(shí)現(xiàn)多學(xué)科協(xié)同仿真分析。現(xiàn)有公開發(fā)表的航行體水動力模型的研究有入水問題[7-8]、潛射導(dǎo)彈出水問題[9-10]、水下航行體的水下航行問題[11-12]等。如潘光等[13]利用有限體積法與非結(jié)構(gòu)動網(wǎng)格技術(shù)，運(yùn)用協(xié)調(diào)一致的壓力耦合方程組半隱式 (SIMPLEC) 方法和k-ε湍流模型(k為湍動能，ε為湍流耗散率)，對無人潛航器進(jìn)行水下對接過程的數(shù)值模擬，得到了在不同攻角下無人水下航行體流體動力隨著距離的變化規(guī)律。Nair等[14]則研究了軸對稱體進(jìn)出水的計算流體力學(xué)(CFD)方法。

通過以上總結(jié)發(fā)現(xiàn)，目前對水下航行體充氣上浮的研究多集中于試驗(yàn)以及氣囊的研究，相關(guān)仿真研究主要涉及剛體的水動力學(xué)研究，同時考慮氣囊展開和航行體上浮姿態(tài)變化的仿真研究尚未有人開展。但是這一問題對于確保航行體的安全可靠回收不可避免，若氣囊體積膨脹帶來的增浮效應(yīng)不夠，導(dǎo)致航行體持續(xù)下沉，并進(jìn)而增加氣囊外部水壓，阻礙氣囊的后續(xù)展開，則可能導(dǎo)致回收失效。而本文研究的航行體上浮過程，不僅需要仿真航行體的剛體運(yùn)動，還需要仿真氣囊的局部大變形。制約這一問題進(jìn)行緊耦合分析的難點(diǎn)在于：現(xiàn)有算法無法同時考慮氣囊展開內(nèi)流場和外部水流場的相互作用，以及氣囊非規(guī)則大幅變形時動態(tài)外流場網(wǎng)格的同步生成和多重流體與固體(簡稱流固)耦合過程的迭代求解。根據(jù)對現(xiàn)有仿真軟件的調(diào)研分析表明，現(xiàn)有軟件無法直接實(shí)現(xiàn)上述過程的仿真。

本文綜合考慮到仿真的可實(shí)現(xiàn)性和計算結(jié)果的有效性，提出一種兩階段的松耦合求解方法：第1階段為折疊氣囊展開過程的模擬。此時不考慮航行體的運(yùn)動姿態(tài)，僅關(guān)注定深水壓作用下氣囊的折疊與展開，獲得助浮氣囊的折疊展開過程中的體積膨脹率時程曲線。第2階段為航行體上浮過程的姿態(tài)運(yùn)動模擬。航行體運(yùn)動過程主要受氣囊展開體積影響，第2階段仿真以第1階段所獲得每個時間步的體積膨脹為輸入，忽略氣囊內(nèi)凹折疊等非規(guī)則變形的影響，采用本文提出的動網(wǎng)格技術(shù)同步表征第1階段計算的氣囊展開過程，并實(shí)時更新水動力學(xué)模型的流場網(wǎng)格，采用單向流固耦合與航行力學(xué)多學(xué)科協(xié)同分析方法來實(shí)現(xiàn)航行體上浮姿態(tài)的分析。

1 充氣展開過程的建模分析

1.1 有限元模型的建立

整個水下航行體幾何尺寸如圖1(a)所示。幾何模型如圖1(b)所示，圖中紅色區(qū)域表示氣囊所在位置。航行體長度為l，直徑為d，氣囊所在區(qū)域?qū)挾葹?.072l，氣囊所在區(qū)域中部距航行體頭部0.146l.

圖1 水下航行體幾何尺寸與幾何模型Fig.1 Geometric dimensions and model of underwater vehicle

整個航行體的剛度相對較大，在本文分析航行體充氣上浮過程的動力學(xué)行為時忽略航行體自身的彈性。實(shí)際充氣上浮的過程包括：發(fā)動機(jī)停車—接收上浮指令—?dú)饽页錃狻獟佌髡帧A(yù)定膨脹體積達(dá)到—排氣穩(wěn)壓維持浮力—航行體上浮至水面。引言部分已經(jīng)指出，完全真實(shí)地考慮航行體水下環(huán)境，以緊耦合方式來仿真模擬航行體的氣囊充氣展開及上浮過程是不現(xiàn)實(shí)的。本文建模分析的第1階段是忽略航行體姿態(tài)變化對折疊氣囊展開過程的影響，僅考慮定水深恒壓環(huán)境，將氣囊展開子系統(tǒng)單獨(dú)分離出來，并應(yīng)用控制體積方法來建模分析折疊氣囊的充氣展開過程。有關(guān)折疊氣囊的充氣展開建模方法在相關(guān)文獻(xiàn)已經(jīng)有較為系統(tǒng)的描述[5-6]，在此簡要交代本文模型的建立，并給出計算結(jié)果，作為第2階段建模分析的輸入。

將展開子系統(tǒng)的各組分模型在三維建模軟件CATIA中建立幾何模型，劃分網(wǎng)格后在有限元前后處理軟件LS-PrePost中合并得到展開子系統(tǒng)的有限元模型。對于展開子系統(tǒng)，其有限元模型主要由折疊環(huán)形氣囊、氣瓶、充氣管以及水下航行體部分剛殼段4部分組成，如圖2所示。由于水下環(huán)形氣囊為三維氣囊，其一般采用較為復(fù)雜的環(huán)向折疊方式安裝在水下航行體上，幾何屬于不可直接展平氣囊，需采用初始矩陣法進(jìn)行有限元網(wǎng)格模型的建立[15-17]。對于該環(huán)形氣囊，為得到單元節(jié)點(diǎn)編號、連接方式一一對應(yīng)的參考網(wǎng)格和映射網(wǎng)格，考慮到該環(huán)形氣囊的參考網(wǎng)格和映射網(wǎng)格均為規(guī)則周期結(jié)構(gòu)，本文采用數(shù)學(xué)分析軟件MATLAB自編程序?qū)崿F(xiàn)兩套網(wǎng)格的建立，如圖3所示。對于深水壓力的影響，通過設(shè)定不同的環(huán)境壓強(qiáng)進(jìn)行模擬，其計算公式為p=pe+ρgh，pe為一個標(biāo)準(zhǔn)大氣壓，取101 300 Pa，ρ為水密度，h為氣囊所在水下深度，g為重力加速度，g=9.8 m/s2.

圖2 展開子系統(tǒng)有限元模型Fig.2 Expanded system finite element model

圖3 環(huán)形氣囊參考網(wǎng)格與映射網(wǎng)格及折疊方式示意圖Fig.3 Reference grid, mapping grid and folding mode of annular airbag

1.2 氣囊展開過程仿真結(jié)果

對氣囊在水下50 m深度的充氣展開過程進(jìn)行仿真分析，得到其充氣展開過程變形圖，如圖4所示。隨著氣瓶充氣，囊內(nèi)壓強(qiáng)不斷升高，整流罩打開，氣囊不斷向外展開，直至2.325 s時氣囊充氣完成，氣囊的體積不再發(fā)生變化。充氣過程中氣囊體積膨脹率與時間的關(guān)系曲線如圖5所示。

圖4 氣囊充氣展開過程變形圖Fig.4 Deformation diagram of airbag during deploying

圖5 氣囊體積膨脹率隨時間變化曲線Fig.5 Changing curve of volume expansion rate of airbag

2 上浮過程的仿真分析方法

對于航行體充氣上浮運(yùn)動的仿真問題，主要是要解決氣囊局部膨脹運(yùn)動與航行體剛體運(yùn)動姿態(tài)變化的耦合運(yùn)動問題，本文提出采用流固耦合和航行力學(xué)的多學(xué)科協(xié)同分析方法來解決這一問題。為避免氣囊初始展開階段的非規(guī)則內(nèi)凹幾何帶來的網(wǎng)格負(fù)體積問題，本文提出采用規(guī)則弧線的旋成體來等效模擬每個時間步的氣囊體積膨脹行為。具體計算時，先根據(jù)圖5曲線提供的每個時間步氣囊體積，確定當(dāng)前時間步等效旋成體的子午線，計算氣囊處網(wǎng)格點(diǎn)的運(yùn)動；然后求解Navier-Stokes方程，得到膨脹過程每一時刻的流場解；再沿著航行體表面對壓力進(jìn)行積分求得航行體所受的力(矩)，將已求得的力(矩)代入6DOF剛體運(yùn)動方程，求解出航行體下一時刻的運(yùn)動學(xué)物理量；利用動網(wǎng)格技術(shù)中的彈簧光順法和局部網(wǎng)格重構(gòu)法，根據(jù)計算出的邊界值自動計算出航行體移動后的網(wǎng)格節(jié)點(diǎn)位置。重復(fù)上述過程，即可對氣囊展開、航行體上浮的整個過程進(jìn)行數(shù)值模擬。本文在CFD軟件Fluent的基礎(chǔ)上，利用用戶自定義函數(shù)(UDF)進(jìn)行編程，實(shí)現(xiàn)了這一耦合求解。求解時每一時間步計算分4個步驟：1)計算氣囊膨脹的等效運(yùn)動；2)求解非定常流場控制方程；3)求解航行體6DOF剛體運(yùn)動；4)更新流場動網(wǎng)格。

2.1 氣囊體積等效膨脹的動網(wǎng)格計算方法

考慮到動網(wǎng)格實(shí)現(xiàn)的可行性，本文將氣囊處網(wǎng)格的等效膨脹簡化處理為截面為圓弧的形狀，氣囊膨脹過程的體積時間曲線由1.2節(jié)氣囊的膨脹過程仿真結(jié)果可得。計算氣囊上網(wǎng)格點(diǎn)的運(yùn)動時需要輸入圓弧對應(yīng)的半徑，所以需要推導(dǎo)出氣囊體積Va(mm3)與截面圓弧半徑r(mm)的關(guān)系式：

當(dāng)圓弧為劣弧時，

(1)

當(dāng)圓弧為優(yōu)弧時，

(2)

具體計算氣囊處網(wǎng)格點(diǎn)的移動時，需要先定義一個結(jié)構(gòu)體數(shù)組，里面分別存儲編號、網(wǎng)格點(diǎn)的三坐標(biāo)(未膨脹前氣囊處網(wǎng)格點(diǎn)坐標(biāo)已知)以及網(wǎng)格點(diǎn)在圓弧上的角度比例，再進(jìn)行旋轉(zhuǎn)角度的計算，確定它在真實(shí)坐標(biāo)系下的坐標(biāo)。

2.1.1 角度比例的計算

圖6為氣囊處點(diǎn)坐標(biāo)運(yùn)動示意圖，其中：C′D′段為氣囊所在位置，A、B分別為航行體的前后端點(diǎn)，計算時以C和D點(diǎn)的軸向位置為依據(jù)，根據(jù)每一個點(diǎn)在CD段上的距離，確定其在角度區(qū)域內(nèi)的比例因子；N′點(diǎn)在CD段上的對應(yīng)位置為N，其在氣囊膨脹過程中，位置不斷更新，分別到了N″、N?、N?′；E′點(diǎn)在CD段上的對應(yīng)位置為E，其在氣囊膨脹過程中，分別到了E″、E?、E?′.

圖6 氣囊處點(diǎn)坐標(biāo)運(yùn)動示意圖Fig.6 Schematic diagram of motion of airbag’s point coordinates

2.1.2 氣囊處網(wǎng)格點(diǎn)坐標(biāo)的更新計算

在計算網(wǎng)格點(diǎn)坐標(biāo)時，需要用到羅德里格旋轉(zhuǎn)公式：設(shè)v是一個三維空間向量，k是旋轉(zhuǎn)軸的單位向量，則v在右手螺旋定則意義下繞旋轉(zhuǎn)軸k旋轉(zhuǎn)角度θ得到的，向量可以由3個不共面的向量v、k和k×v構(gòu)成的標(biāo)架表示：

vrot=cosθv+(1-cosθ)(v·k)k+sinθk×v.

(3)

使用注意事項(xiàng)：所有的向量(轉(zhuǎn)軸、被旋轉(zhuǎn)向量)，均默認(rèn)向量的起始點(diǎn)為原點(diǎn)(0 m, 0 m, 0 m)，并且要求旋轉(zhuǎn)向量和被旋轉(zhuǎn)向量交于原點(diǎn)。如果要求真正的坐標(biāo)，還要加上原來的點(diǎn)坐標(biāo)。如果被旋轉(zhuǎn)向量和轉(zhuǎn)軸向量不相交，則要把被旋轉(zhuǎn)向量分解成兩個向量(被旋轉(zhuǎn)向量起點(diǎn)和轉(zhuǎn)軸向量起點(diǎn)構(gòu)成一個向量；被旋轉(zhuǎn)向量終點(diǎn)和轉(zhuǎn)軸向量起點(diǎn)構(gòu)成一個向量)，然后分別計算向量兩個點(diǎn)旋轉(zhuǎn)后的終點(diǎn)坐標(biāo)，再把這兩個終點(diǎn)坐標(biāo)重新構(gòu)成向量。把二維轉(zhuǎn)角映射到三維上去。在圖6所示的二維圖中，首先確定的向量為AN′、BN′，然后計算它們矢量積的向量AN′×BN′，再確定移動向量N′N″. 把BN′的B點(diǎn)平移到N′點(diǎn)，可構(gòu)成圖7所示的三角形。

圖7 平移BN′得到的三角形Fig.7 Triangle by translating BN′

以此三角形為依據(jù)，向量N′N″減去向量BN′就是剩下一邊的向量。利用余弦定理，可以確定相應(yīng)的角度θ，這樣就在二維平面條件下確定了夾角θ.

在三維條件下，對于每個網(wǎng)格結(jié)點(diǎn)位置，都有對應(yīng)的AN′、BN′，則利用矢量積AN′×BN′，然后繞著這個單位矢量，把BN′旋轉(zhuǎn)角度θ(在二維中已經(jīng)計算得到)，則可以得到沿著N′N″的向量方向；確定了矢量長度|BN′|/|N′N″|的值，則可以直接確定N′N″量各個分量值；然后利用ON″=ON′+N′N″，則可以確定N″的點(diǎn)坐標(biāo)。同樣的方法可求出其余點(diǎn)坐標(biāo)。

2.2 流場求解方法

展開系統(tǒng)上浮運(yùn)動的模擬是考慮剛體運(yùn)動與氣囊體積變化的非定常流體動力學(xué)問題。在模擬計算過程中，需要對空間離散格式、壓力速度耦合格式、時間推進(jìn)格式、湍流模型、考慮氣囊體積變化及剛體運(yùn)動的動網(wǎng)格等進(jìn)行相應(yīng)的設(shè)置和調(diào)整，以實(shí)現(xiàn)大范圍、大變形動態(tài)運(yùn)動過程中，能準(zhǔn)確模擬氣囊展開、水下航行體的上浮過程。

對于存在邊界運(yùn)動的控制體，一般標(biāo)量φ守恒方程的積分形式，在任意控制體積V上，其邊界移動都可以寫成：

(4)

式中：u為流速矢量；ug為動網(wǎng)格的網(wǎng)格速度；Γ為擴(kuò)散系數(shù)；Sφ為φ的源項(xiàng)；?V表示控制體積V的邊界；A是面向量。

在方程中出現(xiàn)的時間導(dǎo)數(shù)用2階時間精度向后差分表示為

(5)

式中：n+1、n、n-1分別表示第n+1個時刻、第n個時刻、第n-1個時刻。控制體積由(6)式進(jìn)行迭代計算：

(6)

為滿足守恒方程(4)式，體積對時間的導(dǎo)數(shù)由(7)式進(jìn)行計算：

(7)

以上對積分形式的Navier-Stokes方程，基于有限體積方法用2階迎風(fēng)格式進(jìn)行了空間離散；然后采用2階歐拉向后隱式格式進(jìn)行時間推進(jìn)；最后采用Transition SST湍流模型來封閉Navier-Stokes方程。

2.3 水下航行體受力分析

水下航行體在完成任務(wù)后，調(diào)整航行的速度、姿態(tài)角，發(fā)動機(jī)停車，氣囊展開產(chǎn)生浮力帶動航行體上浮。上浮過程中其受力情況如圖8所示。圖8中R為航行體和氣囊受到的阻力，F(xiàn)f為航行體與氣囊產(chǎn)生的浮力，mg為航行體自身重力。計算公式分別為

(8)

式中：D、Df分別為航行體與氣囊的阻力系數(shù)；S、Sf分別為航行體與氣囊的浸濕面積；v為水下航行體的速度；Cx為航行體的x軸方向力系數(shù)；Cy為航行體y軸方向的力系數(shù)；Z為航行體運(yùn)動過程受到的側(cè)向力；Cz為航行體z軸方向的力系數(shù)；Vtot為航行體與氣囊的體積。

圖8 水下航行體受力示意圖Fig.8 Force diagram of underwater vehicle in the process of floating movement

氣囊展開過程中會產(chǎn)生逐漸增大的正浮力，使得水下航行體所受到的合力不斷變化，其運(yùn)動姿態(tài)也相應(yīng)地隨時間改變。現(xiàn)有CFD軟件中的6DOF模型主要用于模擬計算域中的剛體在流體力作用下的軌跡及姿態(tài)，相關(guān)計算均建立在物體邊界為剛性邊界的前提條件下。因此，為了模擬水下航行體由于氣囊展開帶動的上浮過程，需用自編UDF對航行體壁面積分，計算流體作用及重力作用共同產(chǎn)生的合力(對質(zhì)心的3個分力)和合力矩(繞體軸的3個分力矩)，再代入到6DOF剛體動力學(xué)方程中，計算出航行體的運(yùn)動軌跡及姿態(tài)。本文動網(wǎng)格設(shè)置與剛體6DOF運(yùn)動程序經(jīng)美國阿諾德工程發(fā)展中心(AEDC)[18]標(biāo)準(zhǔn)彈翼分離算例驗(yàn)證，準(zhǔn)確有效。整個仿真計算的流程圖如圖9所示。

圖9 航行體充氣上浮仿真分析流程圖Fig.9 Flow chart of simulation of underwater vehicle’s floating

3 水下航行體充氣上浮仿真計算

3.1 幾何模型和網(wǎng)格生成

水下航行體幾何模型如圖1所示，航行體模型為殼體，航行體質(zhì)心距頭部距離為0.466l，質(zhì)量為245 kg，航行體浮心距頭部0.464l. 氣囊充氣完成后，航行體與氣囊的總體積為0.294 m3，浮心距航行體頭部距離為0.407l.

在網(wǎng)格劃分方面，航行體采用結(jié)構(gòu)化網(wǎng)格劃分，而在航行體與外流場之間又劃分了一塊橢球形區(qū)域，在航行體與橢球形之間的區(qū)域采用結(jié)構(gòu)化六面體網(wǎng)格劃分，來進(jìn)行氣囊展開的計算，且橢球域與航行體是綁定在一起的，所以航行體上浮時，橢球域與航行體做相同運(yùn)動。對于航行體運(yùn)動所需的空間區(qū)域，則采用四面體非結(jié)構(gòu)化網(wǎng)格進(jìn)行剖分，以充分發(fā)揮四面體網(wǎng)格在動網(wǎng)格方面的優(yōu)勢；在六面體和四面體網(wǎng)格之間，采用五面體(金字塔)網(wǎng)格來實(shí)現(xiàn)過渡，如圖10所示。

圖10 計算域網(wǎng)格示意圖Fig.10 Schematic diagram of computational domain grid

具體進(jìn)行網(wǎng)格劃分時，計算域尺寸的選取相當(dāng)重要，橢球域過大會導(dǎo)致外流場區(qū)域過大，橢球域過小會使氣囊處網(wǎng)格膨脹時擠壓六面體網(wǎng)格導(dǎo)致網(wǎng)格出現(xiàn)負(fù)體積。外流場區(qū)域過大會導(dǎo)致計算量成倍增加，外流場區(qū)域太小會導(dǎo)致計算出現(xiàn)“回流”現(xiàn)象，從而導(dǎo)致計算不收斂。經(jīng)過參考文獻(xiàn)和計算驗(yàn)證，橢球長度方向的直徑定為10l，寬度和高度方向的直徑定為5l. 外流場的尺寸定為100l×29l×57l，如圖11所示。整個模型的網(wǎng)格單元數(shù)大約為600萬。

圖11 計算域外流場示意圖Fig.11 External flow field of computational domain

本文算例航行體航行水深為50 m，發(fā)動機(jī)停車時，航行體前行速度為5 m/s. 進(jìn)行瞬態(tài)計算時，在保證計算收斂和計算精度的前提下，為提高計算效率，需加大時間步長以減少計算時間，本文通過多次收斂性檢驗(yàn)計算后，設(shè)定計算時間步長為1×10-3s. 在計算氣囊膨脹的過程中，即使對氣囊進(jìn)行了規(guī)則化處理，氣囊膨脹時對周圍的結(jié)構(gòu)化六面體網(wǎng)格擠壓還是很嚴(yán)重，在氣囊體積膨脹至52%時，六面體網(wǎng)格由于扭曲率過大而出現(xiàn)負(fù)體積。所以本文在出現(xiàn)負(fù)體積的時間步重新劃分圖10(c)橢球域內(nèi)的六面體網(wǎng)格，進(jìn)行重啟動計算。因?yàn)闄E球域內(nèi)的六面體網(wǎng)格主要進(jìn)行氣囊膨脹的計算，所以重新劃分此區(qū)域內(nèi)的網(wǎng)格并不會影響此后的計算結(jié)果。

3.2 計算結(jié)果分析

3.2.1 航行體運(yùn)動軌跡與姿態(tài)分析

圖12為航行體上浮過程的網(wǎng)格截面示意圖，其直觀展示了航行體充氣上浮過程中網(wǎng)格的變化過程，橢球域與航行體綁定在一起，它們具有相同的運(yùn)動軌跡與姿態(tài)。橢球域內(nèi)的結(jié)構(gòu)化六面體網(wǎng)格用于氣囊膨脹的計算，隨著航行體的轉(zhuǎn)動與上升，橢球域周圍的五面體網(wǎng)格及外流場的四面體網(wǎng)格也逐步變化。當(dāng)網(wǎng)格變形嚴(yán)重扭曲率過大時，采用網(wǎng)格光順與局部網(wǎng)格重構(gòu)的方式重新提高網(wǎng)格質(zhì)量。

仿真獲得的氣囊膨脹和航行體上浮過程如圖13所示，圖14、圖15分別為航行體質(zhì)心運(yùn)動的平面軌跡投影以及三維軌跡。計算結(jié)果表明，在當(dāng)前工況條件下，隨著氣囊的不斷膨脹，航行體所受浮力逐漸增大，氣囊?guī)雍叫畜w由水平位置逐漸轉(zhuǎn)動到豎直位置，氣囊充氣展開提供的增浮作用可有效實(shí)現(xiàn)航行體的上浮回收。

圖13 航行體充氣上浮過程Fig.13 Aeration and floating process of vehicle

圖14 航行體質(zhì)心位移隨時間的變化曲線Fig.14 Displacement of center of mass of vehicle over time

圖15 航行體質(zhì)心在三維空間中的運(yùn)動軌跡Fig.15 Trajectory of center of mass of vehicle

圖16 航行體質(zhì)心線速度隨時間的變化曲線Fig.16 Linear velocity of center of mass of vehicle over time

由圖14x軸方向質(zhì)心位移曲線可知，在大約2.3 s以前，航行體一直沿x軸正向方向運(yùn)動，而2.3 s后航行體又向x軸負(fù)向方向運(yùn)動。航行體質(zhì)心線速度如圖16所示，由圖可知，這是因?yàn)楹叫畜wx軸方向的速度逐漸減小，并且在2.5 s后，x軸方向的速度由0 m/s逐漸向負(fù)方向增加。航行體x軸方向的速度之所以會向負(fù)方向增加，是由于航行體上浮時受到水動力的作用，此時航行體又具有較大俯仰角，導(dǎo)致航行體出現(xiàn)“回旋”，航行體向后運(yùn)動。由航行體質(zhì)心在水深方向(y軸方向)的位移時間曲線可知，在前1.1 s由于航行體重力大于所受浮力，航行體一直處于下降趨勢，而后隨著氣囊膨脹使航行體浮力增大，航行體逐漸上升。航行體開始上升時的氣囊體積膨脹率為62%，由于氣囊膨脹速度相對較快，所以航行體浮力增加很快，航行體在水中并未下降很深便開始上浮。航行體上浮過程中已向z軸負(fù)方向偏移了0.4 m，這是由于航行體受到沿z軸負(fù)方向力的作用。

圖17為航行體運(yùn)動過程中歐拉角隨時間的變化曲線，從中可以看出偏航角在前1.5 s幾乎沒有什么變化。偏航角在1.5 s后逐漸向負(fù)向增加，而后又逐漸沿正向增加。說明航行體先向左偏轉(zhuǎn)，向左偏轉(zhuǎn)了約19°后又向右偏轉(zhuǎn)了20°. 橫滾角開始時變化不大，在3 s后上浮過程才出現(xiàn)較大滾轉(zhuǎn)。雖然開始時浮力還很小，航行體還處于下降趨勢，但由于氣囊膨脹，浮心在質(zhì)心前方不斷向航行體頭部方向移動，浮力產(chǎn)生了抬頭力矩，航行體從開始時俯仰角就一直在增大，直至增大到90°. 而后由于水動力的作用，航行體還會接著向后俯仰，最終達(dá)到106°后又逐漸減小。結(jié)合偏航角的變化，可以看出航行體上浮時處于螺旋上升過程。

圖17 航行體歐拉角隨時間的變化曲線Fig.17 Euler angle of vehicle over time

3.2.2 航行體合力分量的分析與討論

圖18 航行體在3方向上的力系數(shù)Fig.18 Force coefficients of vehicle in three directions

圖19 初始傾角-5°時航行體上浮過程流線圖Fig.19 Streamline diagrams of floating process of vehicle at an initial inclination of -5°

如圖18所示，分別為航行體上浮時合力在地軸系上的分力的無量綱值。圖18(a)為航行體合力在x軸方向分量的變化曲線，開始時航行體所受x軸方向力急劇增大，大約0.5 s時，力增大到最大，隨后逐漸減小。因?yàn)殚_始時航行體具有“停車”速度，并且開始時氣囊膨脹較快，導(dǎo)致航行體所受阻力急劇增大。而后航行體逐漸減速，并且氣囊膨脹速度減慢，這使得航行體所受x軸方向力逐漸減小。圖18(b)為航行體的y軸方向力變化曲線，航行體在充氣上浮過程中，所受的y軸方向力先增大后減小。在前2.3 s，y軸方向力一直增大，并且在1.1 s以前，y軸方向力一直為負(fù)。這是因?yàn)殚_始時，航行體所受重力大于浮力，航行體在豎直方向的合力為負(fù)。隨著氣囊不斷膨脹，航行體俯仰角逐漸增大，1.1 s后航行體y軸方向力沿正向逐漸增大。結(jié)合圖11，航行體在1.1 s后便逐漸上浮，1.1 s時總浮力小于重力，而此時攻角為正，說明航行體上浮時受到水動力提供的豎向力作用。大約2.3 s后，y軸方向力逐漸減小。結(jié)合氣囊體積膨脹率曲線，在2.3 s時，氣囊充氣完成，航行體所受浮力不變。而此時航行體在y軸方向的阻力隨著航行體速度的增大而增大，所以航行體在y軸方向的合力逐漸減小。圖18(c)為航行體z軸方向力的變化曲線，航行體開始運(yùn)動時z軸方向幾乎不受力，而后在上浮過程中逐漸受到一定側(cè)向力作用，這也是為什么航行體上浮過程并不是沿Oxy面完全對稱的原因。但相較于x軸方向、y軸方向分力來說，z軸方向分力的大小相對較小，所以航行體上浮過程中在z軸方向的位移較小，如圖15所示。

3.2.3 不同時刻流場分析

航行體充氣上浮過程所受總力可分為體力與流體力，體力包括航行體自身時不變重力和由氣囊體積確定的時變浮力，流體力是由航行體的運(yùn)動與氣囊的充氣引起的。為了深入了解流體力變化的機(jī)理，繪制了航行體在不同初始俯仰角下充氣上浮過程中的空間流線變化圖，如圖19、圖20所示。

圖20 初始傾角5°時航行體上浮過程流線圖Fig.20 Streamline diagrams of floating process of vehicle at an initial inclination of 5°

從圖19可知：t=0 s時，航行體以-5°入射角移動，由于總力向下，航行體將向下運(yùn)動；在t=0.075 s時，俯仰角速度很小，因此流線在z=0 m平面上是對稱的，在t=0.20 s時，充氣氣囊“源”產(chǎn)生較大的俯仰力矩，但入射角相對較小，虛擬氣囊結(jié)構(gòu)的外傾角不足以產(chǎn)生流體分離，因此流體形態(tài)在z=0 m平面上也是對稱的，在t=0.40 s時，外傾角增大，航行體背風(fēng)面產(chǎn)生對稱分離渦對，在這3個時間點(diǎn)上，俯仰角速度沒有產(chǎn)生尾翼渦脫落，因此流型由平動決定。在t=0.65 s時，由于足夠大的俯仰角速度，非對稱分離渦從尾翼脫落，這種加速的俯仰角速度抑制了氣囊背風(fēng)面的分離渦；在t=1.00 s時，流型相似，但流動的非對稱性引入了橫滾與偏航角速度，從而在航行體頭部和氣囊上觀察到流動的非對稱性。氣囊背風(fēng)面渦對的抑制作用可以維持到俯仰角速度明顯減速為止；在t=1.40 s時，由于俯仰運(yùn)動引起的水阻尼，俯仰角速度明顯減小，因此附著的渦對與氣囊表面分離，形成馬蹄渦[19]。這種馬蹄渦對與尾翼脫落的非對稱渦相互作用，當(dāng)它們順流而下時，馬蹄渦對消失，在氣囊背風(fēng)側(cè)只能看到兩個分離的渦。由于航行體尾翼渦的非對稱性，背風(fēng)渦在運(yùn)動過程中(t=2.00 s、t=2.35 s和t=3.00 s)保持非對稱，并持續(xù)使航行體產(chǎn)生橫滾與偏航運(yùn)動。

從圖20初始傾角5°時流線的變化過程可以看出，流型的演化過程與-5°的情況類似，不同之處在于：與-5°情況相比，5°時航行體向上運(yùn)動導(dǎo)致俯仰運(yùn)動的角加速度增大，從而導(dǎo)致水流分離不對稱提前。

3.2.4 試驗(yàn)驗(yàn)證

根據(jù)本文仿真計算的工況，開展了一系列試驗(yàn)以驗(yàn)證本文的仿真結(jié)果。氣囊在陸地上的展開過程如圖21所示，展開過程與仿真結(jié)果基本相同，氣囊均勻向周邊展開。進(jìn)行水下試驗(yàn)時利用電磁計時器對氣囊展開過程進(jìn)行計時，研究了氣囊在不同深度下的展開過程，仿真與試驗(yàn)展開時間對比如表1所示，從中可以看出誤差均在6 %以內(nèi)，說明仿真結(jié)果的正確性。

圖21 氣囊陸地試驗(yàn)展開過程Fig.21 Deployment process of airbag on the land

進(jìn)行航行體上浮試驗(yàn)的測試時，由于測試設(shè)備的限制，只能利用電磁計時器測量航行體從氣囊開始充氣至上浮完成的時間，結(jié)果如表1所示，從中可以看出誤差均在7 %以內(nèi)，可驗(yàn)證仿真的正確性。

表1 氣囊展開與航行體上浮時間仿真與試驗(yàn)對比

4 結(jié)論

本文主要進(jìn)行了水下航行體利用氣囊充氣上浮的仿真方法研究，先由LS-DYNA軟件利用控制體積算法，模擬氣囊在水下的展開過程，得到氣囊展開過程的體積時間曲線。考慮到CFD計算模型的需要，保持氣囊的體積膨脹率同步等效，用規(guī)則的環(huán)型旋成體等效模擬氣囊的展開過程。利用動網(wǎng)格技術(shù)和自編UDF程序?qū)崿F(xiàn)氣囊膨脹過程的表征，然后基于流體力學(xué)控制方程與剛體6DOF運(yùn)動方程求解航行體的運(yùn)動過程。得出以下主要結(jié)論：

1) 在本文的計算初始條件下，氣囊2.3 s左右就充氣完成，航行體在1.1 s后便開始上浮，在2 s時航行體俯仰角就達(dá)到90°，此后航行體便處于螺旋上升過程。

2) 氣囊膨脹時，航行體所受x軸方向阻力會急劇增大，而后逐漸減小。在上浮過程中，航行體會因?yàn)殇鰷u結(jié)構(gòu)的不對稱性而受到一定側(cè)向力的作用，上浮時并不是沿Oxy平面的對稱上浮運(yùn)動，所以若用二維方法模擬航行體充氣上浮運(yùn)動會不準(zhǔn)確。

3) 航行體上浮時會受到水流提供的豎向力作用，所以上浮前應(yīng)盡量調(diào)整航行體的攻角為正，以便充分利用航行體的殘余速度，從而加快上浮。

該方法能較好地模擬航行體充氣上浮的三維運(yùn)動過程，為水下航行體上浮的工程設(shè)計提供參考，同時本文提出的剛體運(yùn)動與局部大變形耦合的仿真分析方法也為相似問題的解決提供了參考。