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

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

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

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

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

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

1 向量式有限元

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

1.1 質(zhì)點運動方程

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

(1)

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

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

1.2 邊界非線性

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

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

2 流體力學(xué)有限元法

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

2.1 離散方程

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

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

(2)

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

(3)

(4)

(5)

(6)

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

2.2 CFD求解

ANSYS Flotran是Mechanical APDL內(nèi)部基于有限元法的流體力學(xué)計算模塊。與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不需要進(jìn)行人機(jī)交互式界面操作，容易實現(xiàn)程序化模塊并被其他程序調(diào)用，因此本文采用ANSYS Flotran進(jìn)行CFD求解。

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

3 流固耦合分析

3.1 映射網(wǎng)格

在固體域中單元變形和節(jié)點位移遠(yuǎn)小于流體域中的變化量，因此固體網(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將初始模型導(dǎo)入到固體模塊，設(shè)置材料屬性、幾何參數(shù)、約束條件等；同時將初始模型導(dǎo)入到流體模塊，進(jìn)行前處理、設(shè)置邊界條件。

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

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

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

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

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

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

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

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

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

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

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

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

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

4.1 試驗設(shè)計

試驗設(shè)備主要包括囊體、氦氣瓶、導(dǎo)氣管、真空泵和固定裝置。囊體初始幾何模型如圖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。試驗?zāi)Ｐ图翱刂泣c布局如圖4所示。

圖4 試驗?zāi)Ｐ图翱刂泣cFig.4 Test model and control points

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

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

囊體的泄氣采用真空泵抽出氦氣，由于本試驗設(shè)計的氦氣囊囊體氣嘴處于囊體頂部，氦氣密度遠(yuǎn)小于空氣，抽氣過程中囊體內(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；應(yīng)力應(yīng)變輸出步數(shù)100；碰撞檢測步數(shù)10；褶皺處理步數(shù)1；虛擬阻尼100。氣囊的材料參數(shù)如表1所示。

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

流體模塊采用流體力學(xué)有限元法的ANSYS Flotran求解器中計算。泄氣時間和時間步長與固體模塊一致。氦氣域和空氣域均采用Fluid 142單元進(jìn)行計算。將初始模型劃分網(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模型，相應(yīng)的計算參數(shù)取FLOTRAN給定的默認(rèn)值。為了加快收斂速度，采用快速近似求解器求解壓力和速度方程。為了保證求解穩(wěn)定性，人工黏度取0.2以及壓力松弛系數(shù)取0.8。

4.3 對比分析

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

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

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

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

5 結(jié) 論

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