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

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

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

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

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

向量式有限元(vector form of intrinsic finite element，VFIFE)[12-14]是基于向量式固體力學和有限元的數值計算方法。在研究幾何、材料、邊界非線性以及結構靜動力問題比傳統有限元法具有優越性。基于有限元法的計算流體力學計算方法對氣囊形態變化中可能出現的任意幾何形狀的適應性很好，對復雜邊界流域可以實現高階離散精度，因此非常適用氣囊的形態分析問題。

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

1 向量式有限元

氦氣囊是由聚酰胺纖維薄膜制作而成，本文采用常應變三角形膜單元進行氣囊建模。以下簡單介紹向量式有限元膜單元基本理論[15]。

1.1 質點運動方程

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

(1)

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

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

1.2 邊界非線性

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

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

2 流體力學有限元法

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

2.1 離散方程

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

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

(2)

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

(3)

(4)

(5)

(6)

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

2.2 CFD求解

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

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

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

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

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

3 流固耦合分析

3.1 映射網格

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

uf=Hus

(7)

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

δW=δufPf=δusPs

(8)

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

Ps=HTPf

(9)

3.2 雙向流固耦合

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

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

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

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

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

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

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

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

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

步驟9重復步驟6～步驟8，直至該模型趨于穩定。

步驟10改變充盈度，回到步驟2，重新計算氣囊形態；如果出現發散即結束程序并調整初始參數。

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

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

4 試驗設計與數值分析

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

4.1 試驗設計

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

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

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

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

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

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

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

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

4.2 數值模擬

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

表1 氣囊膜材料參數Tab.1 Material properties of membrane

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

4.3 對比分析

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

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

圖6 充盈度為125%時流固界面傳遞信息Fig.6 Flow-solid transferred data at the filling ratio of 125%

圖7 充盈度為80%時流固界面傳遞信息Fig.7 Flow-solid transferred data at the filling ratio of 80%

圖8 充盈度為40%時流固界面傳遞信息Fig.8 Flow-solid transferred data at the filling ratio of 40%

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

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

5 結 論

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