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

張曉光， 李斌， 黨會學， 溫金鵬， 孫潘

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

0 引言

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

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

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

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

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

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

1.1 有限元模型的建立

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

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

整個航行體的剛度相對較大，在本文分析航行體充氣上浮過程的動力學行為時忽略航行體自身的彈性。實際充氣上浮的過程包括：發動機停車—接收上浮指令—氣囊充氣—拋整流罩—預定膨脹體積達到—排氣穩壓維持浮力—航行體上浮至水面。引言部分已經指出，完全真實地考慮航行體水下環境，以緊耦合方式來仿真模擬航行體的氣囊充氣展開及上浮過程是不現實的。本文建模分析的第1階段是忽略航行體姿態變化對折疊氣囊展開過程的影響，僅考慮定水深恒壓環境，將氣囊展開子系統單獨分離出來，并應用控制體積方法來建模分析折疊氣囊的充氣展開過程。有關折疊氣囊的充氣展開建模方法在相關文獻已經有較為系統的描述[5-6]，在此簡要交代本文模型的建立，并給出計算結果，作為第2階段建模分析的輸入。

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

圖2 展開子系統有限元模型Fig.2 Expanded system finite element model

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

1.2 氣囊展開過程仿真結果

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

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

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

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

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

2.1 氣囊體積等效膨脹的動網格計算方法

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

當圓弧為劣弧時，

(1)

當圓弧為優弧時，

(2)

具體計算氣囊處網格點的移動時，需要先定義一個結構體數組，里面分別存儲編號、網格點的三坐標(未膨脹前氣囊處網格點坐標已知)以及網格點在圓弧上的角度比例，再進行旋轉角度的計算，確定它在真實坐標系下的坐標。

2.1.1 角度比例的計算

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

圖6 氣囊處點坐標運動示意圖Fig.6 Schematic diagram of motion of airbag’s point coordinates

2.1.2 氣囊處網格點坐標的更新計算

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

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

(3)

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

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

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

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

2.2 流場求解方法

展開系統上浮運動的模擬是考慮剛體運動與氣囊體積變化的非定常流體動力學問題。在模擬計算過程中，需要對空間離散格式、壓力速度耦合格式、時間推進格式、湍流模型、考慮氣囊體積變化及剛體運動的動網格等進行相應的設置和調整，以實現大范圍、大變形動態運動過程中，能準確模擬氣囊展開、水下航行體的上浮過程。

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

(4)

式中：u為流速矢量；ug為動網格的網格速度；Γ為擴散系數；Sφ為φ的源項；?V表示控制體積V的邊界；A是面向量。

在方程中出現的時間導數用2階時間精度向后差分表示為

(5)

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

(6)

為滿足守恒方程(4)式，體積對時間的導數由(7)式進行計算：

(7)

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

2.3 水下航行體受力分析

水下航行體在完成任務后，調整航行的速度、姿態角，發動機停車，氣囊展開產生浮力帶動航行體上浮。上浮過程中其受力情況如圖8所示。圖8中R為航行體和氣囊受到的阻力，Ff為航行體與氣囊產生的浮力，mg為航行體自身重力。計算公式分別為

(8)

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

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

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

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

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

3.1 幾何模型和網格生成

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

在網格劃分方面，航行體采用結構化網格劃分，而在航行體與外流場之間又劃分了一塊橢球形區域，在航行體與橢球形之間的區域采用結構化六面體網格劃分，來進行氣囊展開的計算，且橢球域與航行體是綁定在一起的，所以航行體上浮時，橢球域與航行體做相同運動。對于航行體運動所需的空間區域，則采用四面體非結構化網格進行剖分，以充分發揮四面體網格在動網格方面的優勢；在六面體和四面體網格之間，采用五面體(金字塔)網格來實現過渡，如圖10所示。

圖10 計算域網格示意圖Fig.10 Schematic diagram of computational domain grid

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

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

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

3.2 計算結果分析

3.2.1 航行體運動軌跡與姿態分析

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

仿真獲得的氣囊膨脹和航行體上浮過程如圖13所示，圖14、圖15分別為航行體質心運動的平面軌跡投影以及三維軌跡。計算結果表明，在當前工況條件下，隨著氣囊的不斷膨脹，航行體所受浮力逐漸增大，氣囊帶動航行體由水平位置逐漸轉動到豎直位置，氣囊充氣展開提供的增浮作用可有效實現航行體的上浮回收。

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

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

圖15 航行體質心在三維空間中的運動軌跡Fig.15 Trajectory of center of mass of vehicle

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

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

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

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

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

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

3.2.3 不同時刻流場分析

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

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

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

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

3.2.4 試驗驗證

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

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

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

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

4 結論

本文主要進行了水下航行體利用氣囊充氣上浮的仿真方法研究，先由LS-DYNA軟件利用控制體積算法，模擬氣囊在水下的展開過程，得到氣囊展開過程的體積時間曲線。考慮到CFD計算模型的需要，保持氣囊的體積膨脹率同步等效，用規則的環型旋成體等效模擬氣囊的展開過程。利用動網格技術和自編UDF程序實現氣囊膨脹過程的表征，然后基于流體力學控制方程與剛體6DOF運動方程求解航行體的運動過程。得出以下主要結論：

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

2) 氣囊膨脹時，航行體所受x軸方向阻力會急劇增大，而后逐漸減小。在上浮過程中，航行體會因為漩渦結構的不對稱性而受到一定側向力的作用，上浮時并不是沿Oxy平面的對稱上浮運動，所以若用二維方法模擬航行體充氣上浮運動會不準確。

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

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