翼傘充氣過程的流固耦合方法數(shù)值仿真

張思宇，余莉，，*，劉鑫

（1.南京航空航天大學(xué) 飛行器環(huán)境控制與生命保障工業(yè)和信息化部重點(diǎn)實(shí)驗(yàn)室，南京210016；2.南京航空航天大學(xué) 航空學(xué)院，南京210016）

沖壓式翼傘相比于傳統(tǒng)降落傘，升阻比高、滑翔穩(wěn)定性能和操作性良好，廣泛應(yīng)用于航空航天等領(lǐng)域［1-3］，近年來成為各國(guó)研究重點(diǎn)。翼傘充氣成功與否決定著翼傘系統(tǒng)工作的成敗。

翼傘和常規(guī)降落傘一樣，在充氣展開過程中，由于本身的柔性特性，在氣流作用下柔性傘衣極易變形，工作特性發(fā)生改變，涉及到的流固耦合問題具有高度非線性、非定常特點(diǎn)，理論研究難度很大。隨著計(jì)算機(jī)的發(fā)展和數(shù)值精度的提高，數(shù)值模擬逐漸成為研究降落傘流固耦合性能的重要手段。Purvis［4］、潘星和曹義華［5］建立了降落傘結(jié)構(gòu)的彈簧阻尼多節(jié)點(diǎn)模型，與內(nèi)部簡(jiǎn)化流場(chǎng)進(jìn)行耦合計(jì)算?；跁r(shí)空離散的MSD-CFD（Mass Spring Damper-Computational Fluid Dynamics）耦合方法［6-8］側(cè)重于傘衣外形的計(jì)算，無法得到結(jié)構(gòu)應(yīng)力信息。Takizawa［9］、Stein［10］和Tezduyar［11］等將DSD/SST（Deform ing-Spatial-Domain/Stabilized Space-Time）有限元計(jì)算方法用于三維降落傘流固耦合仿真計(jì)算研究中，該方法側(cè)重于外部流場(chǎng)的數(shù)值計(jì)算，對(duì)初始充氣段缺少研究。Kim和Peskin［12］采用浸入邊界法進(jìn)行降落傘充氣過程的流固耦合模擬。Tutt［13］和Yu［14］等采用任意拉格朗日-歐拉（Arbitrary Lagrange-Euler，ALE）方法模擬了三維降落傘開傘過程。

但是，目前降落傘的流固耦合研究大多集中于平面圓形傘，針對(duì)翼傘充氣過程的三維流固耦合數(shù)值研究一直沒有突破性進(jìn)展，常常做出一定簡(jiǎn)化，主要面向剛性翼傘［15-16］或者基于松散耦合方法進(jìn)行穩(wěn)態(tài)定常計(jì)算。例如，F(xiàn)ogell等［17］采用弱耦合方法分析了翼傘在穩(wěn)定滑翔階段的氣動(dòng)外形。Kalro等［18］采用并行有限元計(jì)算對(duì)翼傘初始充氣階段進(jìn)行了數(shù)值模擬。汪龍芳和賀衛(wèi)亮［19］對(duì)定常狀況下翼傘的流固耦合變形問題進(jìn)行了三維數(shù)值模擬。張春和曹義華［20］基于弱耦合方法分析了翼傘的氣動(dòng)變形。

翼傘充氣過程大變形、大位移、多尺度流動(dòng)的非定常流固耦合研究基本處于空白階段。這主要是由于：翼傘為非展平曲面構(gòu)成的多氣室-雙翼面復(fù)雜結(jié)構(gòu)，之間由開孔翼肋隔成若干連通氣室，前緣充氣，后緣封閉，結(jié)構(gòu)的復(fù)雜性導(dǎo)致其無法采用常規(guī)方法建立折疊模型。此外，多氣室翼傘在充氣過程中多方向同時(shí)展開，是一個(gè)典型的非穩(wěn)態(tài)流場(chǎng)-柔性傘衣結(jié)構(gòu)的劇烈耦合問題，規(guī)律復(fù)雜。上述問題作為翼傘研究領(lǐng)域的難點(diǎn)，一直備受學(xué)者關(guān)注。

為了解決上述問題，本文提出了基于自由曲面變形（Free Form Deformation，F(xiàn)FD）理論的多氣室柔性沖壓翼傘展向折疊建模方法，流體域通過時(shí)步更新實(shí)現(xiàn)隨傘載系統(tǒng)運(yùn)動(dòng)，基于ALE方法開展折疊翼傘充氣過程的非定常流固耦合數(shù)值計(jì)算，分析翼傘的非線性動(dòng)力學(xué)行為，系統(tǒng)地探究非穩(wěn)態(tài)流場(chǎng)-柔性傘衣結(jié)構(gòu)之間的耦合機(jī)理，為翼傘安全設(shè)計(jì)提供理論依據(jù)。

1 計(jì)算模型

1.1 自由曲面變形理論

翼傘為非展平曲面所構(gòu)成的多氣室-雙翼面結(jié)構(gòu)，為解決其展向折疊問題，本文根據(jù)FFD理論［21］，提出了多氣室柔性沖壓翼傘的展向折疊建模方法。將待變形的翼傘模型嵌入FFD控制晶格中，通過移動(dòng)晶格節(jié)點(diǎn)，將晶格變形傳遞給嵌入其中的翼傘模型，達(dá)到對(duì)其折疊建模的目的。

本文采用六面體控制晶格，如圖1所示，建立局部坐標(biāo)系WSH，W、S、H是局部坐標(biāo)系軸矢量，沿坐標(biāo)軸分別均勻布置l、m、n個(gè)控制節(jié)點(diǎn)，局部坐標(biāo)系原點(diǎn)在全局坐標(biāo)系下的坐標(biāo)記為Q0。

FFD控制晶格各節(jié)點(diǎn)的全局坐標(biāo)為

圖1 FFD六面體控制晶格Fig.1 FFD hexahedral control lattice

式中：（i，j，k）為控制晶格各節(jié)點(diǎn)的局部標(biāo)識(shí)，i＝1，2，…，l；j＝1，2，…，m；k＝1，2，…，n。

在FFD控制晶格中，翼傘任意一點(diǎn)的局部坐標(biāo)是（w，s，h），0≤w，s，h≤1，有

該點(diǎn)的全局坐標(biāo)為Q＝Q0+w W+s S+h H。本文采用Bernstein基函數(shù)建立控制晶格各節(jié)點(diǎn)位移ΔPi，j，k與ΔQ之間的映射函數(shù)關(guān)系：

1.2 流固耦合控制方程

流場(chǎng)域控制方程為

流場(chǎng)域采用無反射邊界條件，以消除流場(chǎng)邊界處的單元速度反射波干擾，防止邊界對(duì)計(jì)算產(chǎn)生影響。

翼傘傘衣、傘繩和流場(chǎng)發(fā)生耦合，采用拉格朗日方法計(jì)算，控制方程為

ALE網(wǎng)格運(yùn)動(dòng)方程為

式中：Li為拉格朗日坐標(biāo)。

對(duì)控制方程式（4）～式（6）進(jìn)行全耦合計(jì)算，采用中心差分時(shí)間顯式方法求解。對(duì)于每一個(gè)節(jié)點(diǎn)，流場(chǎng)和結(jié)構(gòu)的位移x、速度u按式（7）進(jìn)行更新：

式中：Fint、Fext分別為內(nèi)、外力矢量；M 為質(zhì)量對(duì)角矩陣。

1.3 流場(chǎng)域時(shí)步更新

翼傘充氣過程是有限質(zhì)量充氣，為減小計(jì)算消耗，本文在每一時(shí)間步都對(duì)流場(chǎng)域進(jìn)行更新，使其跟隨載荷對(duì)象一起進(jìn)行空間運(yùn)動(dòng)。首先在結(jié)構(gòu)網(wǎng)格上任意選擇不共線的3個(gè)節(jié)點(diǎn)A、B、C，坐標(biāo)分別為KA、KB和KC，建立一個(gè)以A為原點(diǎn)的局部坐標(biāo)系（見圖2），其3個(gè)軸矢量分別為

圖2 流場(chǎng)域時(shí)間步更新Fig.2 Time-step update of flow field

每個(gè)時(shí)間步后，局部坐標(biāo)將跟隨結(jié)構(gòu)網(wǎng)格發(fā)生位移，根據(jù)局部坐標(biāo)在每個(gè)時(shí)間步前后的位移獲得變換矩陣T，則流場(chǎng)網(wǎng)格新的坐標(biāo)系為

式中：［e*1e*2e*31］為更新后的齊次坐標(biāo)；［e1e2e31］為更新前的齊次坐標(biāo)。流體域通過時(shí)間步更新技術(shù)實(shí)現(xiàn)運(yùn)動(dòng)與重構(gòu)。

2 研究對(duì)象與模型驗(yàn)證

2.1 研究對(duì)象與計(jì)算工況

本文采用文獻(xiàn)［22］的沖壓式翼傘模型，如圖3所示，翼梢兩側(cè)有穩(wěn)定幅，32根傘繩，平鋪展長(zhǎng)5.48m，弦長(zhǎng)2.74m，繩長(zhǎng)3.36m，有7個(gè)氣室，每個(gè)氣室由一個(gè)非承載肋片分割成2個(gè)半氣室。翼傘的折疊模型如圖4所示。

本文采用ALE流固耦合數(shù)值方法，將結(jié)構(gòu)網(wǎng)格單元穿插于流場(chǎng)網(wǎng)格中，其流固耦合數(shù)值計(jì)算模型如圖5所示（規(guī)模為6WS×8WS×10WS，其中WS為傘衣弦長(zhǎng)），其中結(jié)構(gòu)附近流場(chǎng)進(jìn)行加密處理，參數(shù)設(shè)置如表1所示。翼傘以0°迎角充氣，初始充氣速度為17m/s。本文基于ANSYS平臺(tái)計(jì)算，采用8核處理器，仿真0.6 s耗時(shí)約172 h。

圖3 沖壓式翼傘結(jié)構(gòu)示意圖Fig.3 Schematic diagram of rammed parafoil structure

圖5 流固耦合數(shù)值計(jì)算模型Fig.5 Numerical calculation model of fluid-structure interaction

表1 流固耦合過程的數(shù)值模型信息Table 1 Numerical model information of fluid-structure interaction process

2.2 模型驗(yàn)證

為驗(yàn)證本文模型的合理性和準(zhǔn)確性，采用文獻(xiàn)［22］的空投試驗(yàn)數(shù)據(jù)進(jìn)行驗(yàn)證計(jì)算，同時(shí)進(jìn)行流場(chǎng)網(wǎng)格的無關(guān)性驗(yàn)證。本文采用3種不同網(wǎng)格密度的流場(chǎng)進(jìn)行數(shù)值計(jì)算，對(duì)比結(jié)果如表2所示。

隨網(wǎng)格數(shù)量增加，計(jì)算誤差減小，但消耗時(shí)長(zhǎng)增加，因此綜合考慮，本文選擇網(wǎng)格數(shù)量為849 000的模型進(jìn)行數(shù)值計(jì)算，得到的翼傘充滿外形和開傘動(dòng)載F（充氣過程中作用在傘衣上的縱向載荷）與空投試驗(yàn)的對(duì)比情況如圖6、圖7所示?？梢钥闯?，翼傘充滿外形與試驗(yàn)外形基本一致。本文開傘動(dòng)載計(jì)算結(jié)果的變化趨勢(shì)與空投試驗(yàn)變化規(guī)律相近，數(shù)值計(jì)算得到的開傘動(dòng)載峰值出現(xiàn)時(shí)刻略早于文獻(xiàn)［22］，數(shù)值為4 760.1 N，和文獻(xiàn)［22］試驗(yàn)數(shù)據(jù)（4 963.2 N）比較，誤差為4.09%，該誤差滿足柔性翼傘開傘計(jì)算的誤差要求。綜上認(rèn)為，本文基于多氣室雙翼面沖壓翼傘建立的數(shù)值計(jì)算模型可以有效模擬翼傘充氣階段工作過程。

表2 網(wǎng)格無關(guān)性對(duì)比結(jié)果Table 2 Comparison results of grid independence

圖6 數(shù)值計(jì)算結(jié)果與空投試驗(yàn)充滿外形對(duì)比Fig.6 Comparison of canopy shape between numerical calculation result and airdrop test

圖7 開傘動(dòng)載隨時(shí)間的變化曲線Fig.7 Curves of opening load over time

3 研究結(jié)果

3.1 傘衣外形及流場(chǎng)變化

翼傘充氣展開過程為有限質(zhì)量充氣情況，許多隨機(jī)參數(shù)都會(huì)對(duì)充氣性能產(chǎn)生影響，本文未考慮風(fēng)場(chǎng)等隨機(jī)環(huán)境因素對(duì)充氣過程的影響。圖8為翼傘充氣過程的外形變化。在充氣初期，翼傘前緣切口進(jìn)氣，尾部最先膨脹展開。之后下翼面高壓氣流作用于穩(wěn)定幅，拉動(dòng)傘衣展向伸展，穩(wěn)定幅產(chǎn)生應(yīng)力集中現(xiàn)象。充氣過程中，前緣切口約束較少，且在氣室充滿后由于切口角度的存在，始終承受氣流吹襲，震顫頻率大，產(chǎn)生持續(xù)的應(yīng)力集中現(xiàn)象。因此，在翼傘設(shè)計(jì)中應(yīng)加強(qiáng)穩(wěn)定幅與前緣切口處的材料強(qiáng)度，避免傘衣破損。

為進(jìn)一步探究翼傘充氣過程中的結(jié)構(gòu)及流場(chǎng)的非定常變化規(guī)律，選取a、b、c 3個(gè)截面（見圖9）進(jìn)行研究，各個(gè)截面的結(jié)構(gòu)及流場(chǎng)變化如圖10所示（各個(gè)時(shí)刻從上到下依次為a、b、c截面），為刻畫細(xì)節(jié)情況，各個(gè)時(shí)刻的圖形縮比尺寸不同。

圖8 充氣過程翼傘外形變化（范式等效應(yīng)力云圖）Fig.8 Parafoil shape changes during inflation process（von M ises equivalent stress contours）

圖9 選取截面示意圖Fig.9 Schematic diagram of selected sections

首先，通過觀察圖10所示翼傘充氣過程的弦向變化（a、b截面）：在翼傘初始折疊狀態(tài)時(shí)，a、b截面的速度場(chǎng)分布相似，前緣產(chǎn)生渦。隨翼傘逐漸充氣，上翼面的渦逐漸后移至脫體，前緣處不斷產(chǎn)生新的渦。在充氣過程中，a截面隨氣室充氣逐漸形成飽滿翼型，內(nèi)腔壓力增加，傘衣剛度增加，類似剛性翼，外部流暢穩(wěn)定。而b截面處于最外側(cè)，傘衣由于約束較少變形明顯，在翼傘充滿后翼型飽滿度不足，外部流場(chǎng)不穩(wěn)定，傘衣變形狀態(tài)具有一定的隨機(jī)性。

圖10中翼傘充氣過程的展向變化（c截面）表明：在翼傘初始折疊狀態(tài)時(shí)，c截面周圍流場(chǎng)產(chǎn)生對(duì)稱繞流現(xiàn)象，此時(shí)傘衣折疊壓縮率大，下翼面形成高壓區(qū)（壓力云圖見圖11）。之后翼傘底部對(duì)稱渦開始破裂，氣流繞翼尖流出，翼傘展開充氣。該過程上翼面流速快，下翼面流速慢，上下翼面形成的壓力造成“翼尖上翹，中部凹陷”的翼傘尾流再附現(xiàn)象。因此在翼傘開傘過程中應(yīng)加入引導(dǎo)傘改善充氣狀況，引導(dǎo)傘拉出傘衣中部減輕氣室塌陷。此外，設(shè)計(jì)合理的穩(wěn)定幅，阻擋下翼面到上翼面的繞流，消弱翼尖渦強(qiáng)度，減輕干擾影響。此后，傘衣隨氣室充滿保持穩(wěn)定充滿翼型。

圖11 充氣過程壓力云圖Fig.11 Pressure contours of inflation process

3.2 氣室充氣規(guī)律

為進(jìn)一步研究翼傘的各氣室充氣規(guī)律，將翼傘沿弦向分為前（前緣）、中（中部）、后（尾緣）三部分，沿展向?qū)?個(gè)氣室逐次編號(hào)，兩側(cè)為氣室1和氣室7，中央氣室為4號(hào)。圖12為各氣室充滿時(shí)間和充滿寬度的變化規(guī)律?？梢钥闯觯?/p>

1）各氣室的充氣規(guī)律關(guān)于中央氣室具有對(duì)稱性，中央氣室充氣快且充滿外形飽滿，越靠近兩側(cè)充滿速率越慢，充氣效果越差。

2）在高壓氣流推動(dòng)作用下，進(jìn)入氣室的氣體涌入尾部，翼傘尾緣最先展向打開，充氣速率明顯快于其他部分，前緣展開最慢。充滿后中部充氣效果最好，氣室最為飽滿。

圖12 各氣室充滿時(shí)間和充滿寬度Fig.12 Inflation time and width of each air chamber

3）由于尾緣上下翼面縫合，翼尖繞流導(dǎo)致傘衣變形，內(nèi)外氣壓不穩(wěn)定，兩側(cè)氣室的后部未完全充滿。

充氣過程的翼展變化見圖13。翼傘的前緣、中部、尾緣前、中、后三部分于0.39、0.37、0.25 s分別達(dá)到翼展峰值5.05、5.40、5.20m。之后由于傘繩和傘衣的彈性作用，各部位充滿后產(chǎn)生一定的回彈和氣室“鼓包”現(xiàn)象，中部充滿翼展穩(wěn)定在5.0m左右，略小于設(shè)計(jì)值。

圖13 充氣過程的翼展變化Fig.13 Changes in wingspan during inflation process

3.3 氣動(dòng)特性分析

翼傘的氣動(dòng)特性參數(shù)變化情況如圖14所示。0.1 s之前，計(jì)算域尚未形成穩(wěn)定繞流流場(chǎng)，氣動(dòng)參數(shù)不具有參考意義。隨繞流流場(chǎng)形成，翼傘的升阻力特性參數(shù)變化趨勢(shì)一致，在0.28 s時(shí)達(dá)到峰值，此時(shí)傘衣投影面積最大。之后由于呼吸回彈現(xiàn)象，傘衣面積減小，氣動(dòng)系數(shù)減小，翼傘充滿穩(wěn)定后氣動(dòng)系數(shù)略有回升。充滿后傘載系統(tǒng)存在一定迎角，滑翔比CL／CD穩(wěn)定在2.24，升力系數(shù)CL為1.56，阻力系數(shù)CD為0.70。

圖14 氣動(dòng)特性參數(shù)的變化Fig.14 Changes in aerodynamic characteristic parameter

4 結(jié) 論

本文針對(duì)非定常充氣展開過程，數(shù)值模擬了翼傘的三維流固耦合動(dòng)力學(xué)特性，得出：

1）翼傘充氣展開過程中穩(wěn)定幅和前緣切口處產(chǎn)生應(yīng)力集中現(xiàn)象，因此，在翼傘設(shè)計(jì)中應(yīng)加強(qiáng)穩(wěn)定幅與前緣切口處的材料強(qiáng)度。

2）翼傘充氣展開過程存在“翼尖上翹，中部凹陷”的翼傘尾流再附現(xiàn)象，造成傘衣內(nèi)塌，導(dǎo)致系統(tǒng)穩(wěn)定性受影響，因此在翼傘開傘過程中應(yīng)加入引導(dǎo)傘以改善充氣狀況，設(shè)計(jì)合理的穩(wěn)定幅可以減小翼尖渦的影響。

3）各氣室的充氣規(guī)律關(guān)于中央氣室具有對(duì)稱性，中央氣室充氣快且充滿外形飽滿，兩側(cè)充滿速率越慢，氣室飽滿度下降；翼傘尾部最先膨脹展開，前緣展開最慢，充滿后存在一定的回彈和氣室“鼓包”現(xiàn)象，充滿后展長(zhǎng)小于設(shè)計(jì)值。

4）升阻力特性參數(shù)變化趨勢(shì)一致，傘衣充滿后翼傘滑翔比穩(wěn)定在2.24，升力系數(shù)為1.56，阻力系數(shù)為0.70。