基于流固耦合的重載飛艇副氣囊形態(tài)分析與試驗研究

艾科熱木江·塞米， 袁行飛， 陳務(wù)軍， 王雪明， 張大旭

(1. 浙江大學 空間結(jié)構(gòu)研究中心，杭州 310058；2. 上海交通大學 空間結(jié)構(gòu)研究中心，上海 200240； 3. 中航工業(yè)特種飛行器研究所，湖北 荊門 448035)

飛艇是一種依靠氣囊中浮力工質(zhì)氣體提供靜浮力，依靠發(fā)動機和螺旋槳提供推力的飛行器。軟式飛艇由直接作為外蒙皮的主氣囊，以及產(chǎn)生內(nèi)外壓差來保持艇體外形的副氣囊組成。重載飛艇和平流層飛艇內(nèi)置氣囊采用的浮升氣為氦氣[1]。氦氣囊在不同高度、巡航階段呈現(xiàn)不同的充盈度和形態(tài)，氦氣囊充盈度將決定浮力的大小、浮心位置變化，進而影響飛艇姿態(tài)控制[2-3]，目前主要是基于經(jīng)驗考慮，尚難定量分析，但是針對新型重載飛艇或大型平流層飛艇，建立有效數(shù)值模擬方法尤為必要。針對氦氣囊或者囊體非穩(wěn)定形態(tài)的研究較少，文獻[4]結(jié)合Abaqus軟件和FLUENT軟件，對飛艇外氣囊與周圍空氣流場進行流固耦合分析。文獻[5]研究飛艇主氣囊的模態(tài)特性及試驗驗證模態(tài)分析方法。文獻[6]針對飛艇非保形升空時囊體非保形的形態(tài)進行模擬分析。文獻[7]以正高斯曲率氦氣囊為對象，采用CFX軟件分析不同仰角下的非飽和形態(tài)。重載飛艇內(nèi)置氦氣囊構(gòu)型復雜，升空和循環(huán)過程形態(tài)變化和力學行為復雜，包括大變形、褶皺、接觸干涉等非線性，以及內(nèi)部氦氣和外部空氣及與囊體耦合作用，然而相關(guān)理論與數(shù)值分析鮮有報道。

流固耦合[8]問題最初在航空航天領(lǐng)域被提出之后，多種分析理論和計算方法在交叉學科領(lǐng)域中取得了快速發(fā)展。文獻[9]提出了基于交錯積分的弱耦合方法，文獻[10]進一步研究了基于同步迭代的弱耦合方法，這些方法將固體域和流體域控制方程分開求解后進行CSD(computational structural dynamics)/CFD(computational fluid dynamics)之間數(shù)據(jù)傳遞。文獻[11]研究了在同一控制方程中求解所有變量的強耦合方法。然而此方法對復雜邊界的適應性低，計算工作量大等原因，大部分主流的CFD分析軟件采用基于弱耦合理論的雙向流固耦合分析方法。

向量式有限元(vector form of intrinsic finite element，VFIFE)[12-14]是基于向量式固體力學和有限元的數(shù)值計算方法。在研究幾何、材料、邊界非線性以及結(jié)構(gòu)靜動力問題比傳統(tǒng)有限元法具有優(yōu)越性?；谟邢拊ǖ挠嬎懔黧w力學計算方法對氣囊形態(tài)變化中可能出現(xiàn)的任意幾何形狀的適應性很好，對復雜邊界流域可以實現(xiàn)高階離散精度，因此非常適用氣囊的形態(tài)分析問題。

本文將向量式有限元擅長解決強非線性問題以及流體力學有限元法適用于任意幾何邊界條件的特點結(jié)合起來，采用雙向流固耦合理論分析飛艇內(nèi)置氦氣囊的不穩(wěn)定形態(tài)變化。同時開展縮尺比例模型試驗，與數(shù)值分析結(jié)果加以對比，驗證了本文提出的分析方法的適用性和準確性。

1 向量式有限元

氦氣囊是由聚酰胺纖維薄膜制作而成，本文采用常應變?nèi)切文卧M行氣囊建模。以下簡單介紹向量式有限元膜單元基本理論[15]。

1.1 質(zhì)點運動方程

向量式有限元首先將結(jié)構(gòu)分解為一組離散的空間點組成的質(zhì)點群，相鄰質(zhì)點之間的參數(shù)用一組標準化的內(nèi)插函數(shù)表示，從而避免了有限單元法用微分形式的連續(xù)函數(shù)描述結(jié)構(gòu)。其次，用一組時間點組成的控制方程描述空間點的運動過程。控制方程中的質(zhì)點間相互作用力僅與相對位置中的純變形有關(guān)，因此使結(jié)構(gòu)作虛擬逆向運動，從總位移中扣除剛體運動分量得到純變形，進而求解單元內(nèi)力，得到結(jié)構(gòu)的內(nèi)力分布。質(zhì)點運動的控制方程為牛頓第二定律方程式，即

(1)

豎向壓力按氣囊內(nèi)外氣體的密度差隨高度呈線性分布，同一高度氣體環(huán)向壓力近似認為處處相等，如圖1所示。氣囊在不穩(wěn)定形態(tài)變化中，某一高度下的內(nèi)外壓強相等，稱之為零壓面。由于氦氣密度比空氣小，氣囊在零壓面以上處于膨脹(超壓)狀態(tài)，而零壓面以下則處于收縮(負壓)狀態(tài)。考慮到幾何形狀、流場變化以及外界因素等將對實際壓力梯度產(chǎn)生較大的影響，因此模擬分析中先按假設(shè)的壓力梯度求出節(jié)點位移后接著通過流體力學計算獲得真實壓力分布，并修正節(jié)點外力矩陣。

圖1 壓力梯度示意圖Fig.1 Schematic diagram of pressure gradient

1.2 邊界非線性

氣囊在泄氣過程中將產(chǎn)生褶皺和碰撞接觸等復雜力學行為，因此單位時間內(nèi)需檢測接觸單元并對其進行碰撞響應處理。碰撞檢測是整個邊界非線性問題中最耗時的工作，本文首先對整體空間域進行較粗略的全局檢測，然后對各個子空間域進行精確的局部檢測。全局檢測是把整個囊體所在空間域劃分成若干大小相等的子域之后，確定所有節(jié)點和單元所屬子域，單元可按形心點坐標劃分所屬子域。然后存儲每個子域所包含的節(jié)點和單元編號。由于向量式有限元時間步長非常小，假設(shè)單元只可能與所在子域及周圍3×3×3子域內(nèi)的節(jié)點發(fā)生碰撞，從而大大減小了搜索計算量。局部檢測是前后兩個時刻的節(jié)點位于單元兩側(cè)，以及穿透點位于單元內(nèi)部時可認為發(fā)生碰撞，此時兩個三角形的接觸可轉(zhuǎn)化為點-三角形接觸問題。

確定發(fā)生碰撞的節(jié)點和單元之后，計算虛擬接觸力進行碰撞響應處理。本文采用罰函數(shù)法[16]計算法向排斥力和切向摩擦力，然后將該力反向施加到發(fā)生碰撞的節(jié)點和單元，修正節(jié)點力矩陣并重新求解質(zhì)點運動方程。

2 流體力學有限元法

流體力學有限元法的基本思想是將流體域劃分成有限個離散的子區(qū)域，對于各子域按照一定規(guī)則構(gòu)造基函數(shù)，對整個單元區(qū)域進行積分得出單元有限元方程，從而使用各子域差分方程的近似解逼近微分方程的精確解[17]。

2.1 離散方程

流體力學基本方程主要包括連續(xù)方程、動量方程、能量方程等，按研究對象的復雜程度可以考慮理想氣體狀態(tài)方程、組分方程和湍流方程使總控制方程閉合。使用φ表示通用變量，則離散化的單元矩陣[18]可表示為

[Te+Ae+De]{φe}={Se}

(2)

式中：T為瞬態(tài)項；A為對流項；D為擴散項；S為源項。將單元權(quán)函數(shù)We代入式(2)中的各項，在單元內(nèi)積分表達式為

(3)

(4)

(5)

(6)

式中：ρ,Ω,Γ,v分別為單元密度、體積、擴散系數(shù)和速度；矢量符號div為散度； grad為梯度。

2.2 CFD求解

ANSYS Flotran是Mechanical APDL內(nèi)部基于有限元法的流體力學計算模塊。與Fluent，CFX等主流流體計算軟件相比，F(xiàn)lotran具有以下特點：

(1) Flotran采用有限元法，而Fluent，CFX等軟件采用有限體積法。在本文流固耦合計算中，需要跟蹤單元網(wǎng)格上各節(jié)點的速度、壓強等變量信息。有限元法中的離散方程定義明確，與向量式有限元結(jié)合運用，使流固耦合計算更簡便。

(2) Flotran和MATLAB基于Fortran語言開發(fā)，而Fluent，CFX等軟件基于C語言開發(fā)。由于本文中向量式有限元采用自編MATLAB程序，因此與Flotran配合可使編程語言統(tǒng)一化，計算更簡便。

另外，F(xiàn)lotran不需要進行人機交互式界面操作，容易實現(xiàn)程序化模塊并被其他程序調(diào)用，因此本文采用ANSYS Flotran進行CFD求解。

Flotran的CFD求解過程是首先假定一個初始流速分布以求解壓力方程，得到速度和壓力后代入組分傳輸方程和湍流方程，然后用求出的湍流動能k和湍流動能耗散率ε更新有效黏度，進入下一步迭代直至各自由度收斂。

3 流固耦合分析

3.1 映射網(wǎng)格

在固體域中單元變形和節(jié)點位移遠小于流體域中的變化量，因此固體網(wǎng)格比流體網(wǎng)格更稀疏。此時需要插值函數(shù)將固體網(wǎng)格和流體網(wǎng)格上的節(jié)點力、位移建立映射關(guān)系，為CFD/CSD之間的數(shù)據(jù)交換起橋梁作用。本文采用一種基于薄板樣條(thin plate spline，TPS)函數(shù)[19]的插值方法。us,uf分別為固體網(wǎng)格和流體網(wǎng)格上的節(jié)點位移，則

uf=Hus

(7)

式中，H為插值矩陣。同理，Ps,Pf分別為固體網(wǎng)格和流體網(wǎng)格節(jié)點上的荷載，由虛功原理得出

δW=δufPf=δusPs

(8)

將式(7)代入式(8)可得

Ps=HTPf

(9)

3.2 雙向流固耦合

基于流固耦合的氦氣囊形態(tài)分析具體步驟如下：

步驟1將初始模型導入到固體模塊，設(shè)置材料屬性、幾何參數(shù)、約束條件等；同時將初始模型導入到流體模塊，進行前處理、設(shè)置邊界條件。

步驟2固體模塊按照氦氣充盈度和囊體體積確定零壓面高度，計算壓力梯度和初始節(jié)點外力矩陣。

步驟3中央差分法計算節(jié)點位移，更新節(jié)點坐標。若發(fā)生接觸碰撞，計算接觸力并修正節(jié)點力矩陣，重新計算節(jié)點位移。

步驟4從節(jié)點位移依次計算各單元的純變形、應變、應力和內(nèi)力，最后反作用于節(jié)點上得到節(jié)點內(nèi)力。

步驟5由重力和壓力梯度計算節(jié)點外力，集成節(jié)點力矩陣。

步驟6重復步驟3～步驟5，直至模型趨于平衡，然后將步驟3得到的節(jié)點位移通過映射網(wǎng)格傳送到流體模塊。

步驟7更新流體模塊中的節(jié)點坐標，依次求解N-S方程、組分傳輸方程和湍流方程直至各自由度收斂。

步驟8得到的節(jié)點壓力分布通過映射網(wǎng)格傳送到固體模塊修正壓力梯度，回到步驟5。

步驟9重復步驟6～步驟8，直至該模型趨于穩(wěn)定。

步驟10改變充盈度，回到步驟2，重新計算氣囊形態(tài)；如果出現(xiàn)發(fā)散即結(jié)束程序并調(diào)整初始參數(shù)。

具體計算流程如圖2所示，圖中左邊區(qū)域為固體模塊，在自編MATLAB程序中進行；右邊區(qū)域為流體模塊，在ANSYS軟件的Flotran模塊中進行。

圖2 雙向流固耦合流程圖Fig.2 Flow chart of two-way fluid-solid coupling

4 試驗設(shè)計與數(shù)值分析

為了研究氣囊隨氦氣充盈度變化導致的形態(tài)變化規(guī)律，本文采用飛艇內(nèi)置氣囊專用聚酰胺纖維材料制作了一20∶1比例仿真模型進行泄氣試驗，并采用本文建立的雙向流固耦合方法進行數(shù)值模擬，最后將模型試驗和數(shù)值模擬結(jié)果進行對比。

4.1 試驗設(shè)計

試驗設(shè)備主要包括囊體、氦氣瓶、導氣管、真空泵和固定裝置。囊體初始幾何模型如圖3(a)所示。囊體前、后、側(cè)面均布置用于位移測量的靶點。囊體在泄氣過程中，如約束數(shù)量不夠，會出現(xiàn)一直晃動而影響靶點的位移測量。因此本試驗布置了囊體上的18個xyz方向全約束的固定點：頂部區(qū)域共三排固定點，每排三個；三角網(wǎng)處共計 9 個固定點沿三角線分布，如圖3(b)所示。

圖3 初始幾何模型與固定點示意圖Fig.3 Schematic diagram of initial geometric model and fixed points

本試驗中分別在氣囊正面、背面和側(cè)面上選取若干靶點作為控制點：正面-7，8，9；背面-21，23，30；側(cè)面-27。試驗模型及控制點布局如圖4所示。

圖4 試驗模型及控制點Fig.4 Test model and control points

采用三臺照相機對囊體同一個面上的靶點分別從三個角度進行拍攝，最后通過計算處理得到靶點相對標定面的實際空間位置。依照該方法可以得到囊體表面在不同泄氣狀態(tài)下相對同一標定面(即同一空間直角坐標系)的空間位置。對比初始泄氣狀態(tài)下的空間位置，可以得到不同泄氣狀態(tài)下的靶點的空間位移，從而得到囊體在不同氦氣充盈度下的形態(tài)位移變化。

本試驗中，首先將氦氣囊充氣至初始幾何參考模型的1.25倍，即125%充盈度，然后通過預設(shè)的固定點固定住氦氣囊進行泄氣試驗，測量氣囊上的控制點在不同氦氣充盈度下的位移，得到囊體隨著氦氣充盈度減小的形態(tài)變化規(guī)律。

囊體的泄氣采用真空泵抽出氦氣，由于本試驗設(shè)計的氦氣囊囊體氣嘴處于囊體頂部，氦氣密度遠小于空氣，抽氣過程中囊體內(nèi)剩余氦氣始終聚集于囊體頂部區(qū)域，因此最后囊體內(nèi)的氦氣可以被真空泵抽取完。試驗中采用分級泄氣控制方案，將氦氣囊從125%充盈度狀態(tài)抽氣直到抽干為止，一共分19級泄氣，第一次抽氣時間為5 min，其余各級抽取時間為3.6 min。泄氣試驗過程如圖5所示。

圖5 泄氣試驗過程Fig.5 Process of deflation experiment

4.2 數(shù)值模擬

固體模塊采用向量式有限元理論在自編MATLAB程序中計算。每次泄氣計算時間0.5 s；時間步長1×10-4s；各充盈度下的循環(huán)步數(shù)5 000；應力應變輸出步數(shù)100；碰撞檢測步數(shù)10；褶皺處理步數(shù)1；虛擬阻尼100。氣囊的材料參數(shù)如表1所示。

表1 氣囊膜材料參數(shù)Tab.1 Material properties of membrane

流體模塊采用流體力學有限元法的ANSYS Flotran求解器中計算。泄氣時間和時間步長與固體模塊一致。氦氣域和空氣域均采用Fluid 142單元進行計算。將初始模型劃分網(wǎng)格后總單元數(shù)量為165 988個，其中氣囊為3 928個，氦氣域為32 069個，空氣域為129 991個。時間步長為1×10-4s，計算時間為0.5 s。氣囊的泄氣過程屬于瞬態(tài)問題，打開大變形和動網(wǎng)格ALE開關(guān)。假設(shè)不可壓縮氣體，系統(tǒng)與外界沒有熱量交換，因此不考慮能量守恒方程。由于氣囊內(nèi)部氦氣與外界空氣之間存在質(zhì)的交換，考慮組分守恒方程。本文采用湍流模型中的RNG模型，相應的計算參數(shù)取FLOTRAN給定的默認值。為了加快收斂速度，采用快速近似求解器求解壓力和速度方程。為了保證求解穩(wěn)定性，人工黏度取0.2以及壓力松弛系數(shù)取0.8。

4.3 對比分析

基于本文建立的雙向流固耦合方法計算得到的氦氣囊形態(tài)變化趨勢與本文試驗結(jié)果基本保持一致。氦氣從頂部區(qū)域的氣嘴抽出時，囊體的收縮凹陷首先發(fā)生在下部區(qū)域，繼而發(fā)生褶皺，相反，氦氣囊頂部區(qū)域附近仍處于充盈狀態(tài)；隨著氦氣的持續(xù)抽出，囊體的整體充盈度不斷降低，囊體下部收縮凹陷程度不斷增加，褶皺增多，形成復雜不規(guī)則的收縮凹陷，同時囊體下部區(qū)域的收縮凹陷狀態(tài)逐漸擴展到上部區(qū)域，使得部分中、上部區(qū)域也產(chǎn)生收縮、凹陷和褶皺。在囊體的泄氣過程中，其固定點處由于受到近似固端約束始終處于繃緊狀態(tài)，其區(qū)域附近位移相對很小。

由于篇幅限制，本文僅給出充盈度分別為125%，80%，40%時，向量式有限元得到的氦氣囊位移云圖、ANSYS Flotran模塊得到的壓力云圖以及模型試驗拍攝的氦氣囊形態(tài)，如圖6～圖8所示。將不同氦氣充盈度下的氦氣囊模型試驗與數(shù)值分析結(jié)果相比較，得出泄氣過程中的氣囊形態(tài)變化規(guī)律是一致的，從而驗證了本文提出的雙向流固耦合方法的適用性和準確性。

圖6 充盈度為125%時流固界面?zhèn)鬟f信息Fig.6 Flow-solid transferred data at the filling ratio of 125%

圖7 充盈度為80%時流固界面?zhèn)鬟f信息Fig.7 Flow-solid transferred data at the filling ratio of 80%

圖8 充盈度為40%時流固界面?zhèn)鬟f信息Fig.8 Flow-solid transferred data at the filling ratio of 40%

圖9(a)～圖9(c)為三個靶點面(正、背、側(cè)面)上選取的控制點相對初始狀態(tài)的距離在泄氣過程中的變化曲線。下橫坐標是數(shù)值模擬結(jié)果中的氦氣充盈度，而上橫坐標是模型試驗中的泄氣序號。由于向量式有限元中的充盈度是根據(jù)每一步泄氣時初始零壓面的高度確定的，而迭代循環(huán)結(jié)束后的零壓面高度往往發(fā)生微小變化，因此得到的氦氣充盈度有±3%以內(nèi)的誤差。在模型試驗中，理論上根據(jù)真空泵的流量(5 m3/h)可以算出抽取一定量的氦氣需要的時間，從而控制每級泄氣抽出的氦氣體積。但是隨泄氣程度增大，抽氣所需時間變長，利用真空泵的流量計算抽取出來的氦氣體積的誤差越來越大。因此模型試驗和數(shù)值模擬的計算時刻不能一一對應。但是在同一張時間-位移圖中可以看出，所有控制點隨著氦氣囊泄氣程度的增大，其相對初始狀態(tài)的距離變化趨勢基本上保持一致，從而證明了本文提出數(shù)值模擬方法在研究氣囊隨氦氣充盈度變化導致的非穩(wěn)定形態(tài)變化規(guī)律的適用性。

圖9 控制點相對初始狀態(tài)的距離在泄氣過程中的變化曲線Fig.9 Displacement of control points from initial state at different filling ratio

5 結(jié) 論

本文為揭示重載飛艇氦氣囊力學行為，采用向量式有限元CSD和ANSYS Flotran的CFD，通過映射實現(xiàn)流固界面數(shù)據(jù)傳遞，考慮囊體大變形、褶皺、接觸非線性和氦氣浮力效應，建立了副氣囊流固耦合數(shù)值模擬方法。對氦氣囊泄氣過程進行了數(shù)值模擬，得到了氦氣囊在不同氦氣充盈度下的形態(tài)變化規(guī)律，囊體底部區(qū)域隨泄氣首先發(fā)生收縮凹陷，頂部區(qū)域附近變形較小；隨著泄氣量增加，囊體沿高度方向從下至上依次收縮凹陷，凹陷之間產(chǎn)生擠壓繼而形成褶皺；囊體表面控制點隨泄氣程度的增大，其相對初始狀態(tài)的距離呈現(xiàn)增大趨勢，部分控制點在泄氣過程中相對初始狀態(tài)的距離會有小幅度的降低。同時開展了氦氣囊的縮尺比例模型試驗，對囊體在泄氣過程中得到的位移測量值與數(shù)值模擬結(jié)果進行對比。結(jié)果表明，氦氣囊在不同充盈度下的形態(tài)變化與雙向流固耦合數(shù)值分析結(jié)果一致，證明了本文提出的數(shù)值模擬方法可有效用于氣囊隨氦氣充盈度變化導致的非穩(wěn)定形態(tài)變化規(guī)律的研究。