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

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

（1.南京航空航天大學 飛行器環境控制與生命保障工業和信息化部重點實驗室，南京210016；2.南京航空航天大學 航空學院，南京210016）

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

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

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

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

為了解決上述問題，本文提出了基于自由曲面變形（Free Form Deformation，FFD）理論的多氣室柔性沖壓翼傘展向折疊建模方法，流體域通過時步更新實現隨傘載系統運動，基于ALE方法開展折疊翼傘充氣過程的非定常流固耦合數值計算，分析翼傘的非線性動力學行為，系統地探究非穩態流場-柔性傘衣結構之間的耦合機理，為翼傘安全設計提供理論依據。

1 計算模型

1.1 自由曲面變形理論

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

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

FFD控制晶格各節點的全局坐標為

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

式中：（i，j，k）為控制晶格各節點的局部標識，i＝1，2，…，l；j＝1，2，…，m；k＝1，2，…，n。

在FFD控制晶格中，翼傘任意一點的局部坐標是（w，s，h），0≤w，s，h≤1，有

該點的全局坐標為Q＝Q0+w W+s S+h H。本文采用Bernstein基函數建立控制晶格各節點位移ΔPi，j，k與ΔQ之間的映射函數關系：

1.2 流固耦合控制方程

流場域控制方程為

流場域采用無反射邊界條件，以消除流場邊界處的單元速度反射波干擾，防止邊界對計算產生影響。

翼傘傘衣、傘繩和流場發生耦合，采用拉格朗日方法計算，控制方程為

ALE網格運動方程為

式中：Li為拉格朗日坐標。

對控制方程式（4）～式（6）進行全耦合計算，采用中心差分時間顯式方法求解。對于每一個節點，流場和結構的位移x、速度u按式（7）進行更新：

式中：Fint、Fext分別為內、外力矢量；M 為質量對角矩陣。

1.3 流場域時步更新

翼傘充氣過程是有限質量充氣，為減小計算消耗，本文在每一時間步都對流場域進行更新，使其跟隨載荷對象一起進行空間運動。首先在結構網格上任意選擇不共線的3個節點A、B、C，坐標分別為KA、KB和KC，建立一個以A為原點的局部坐標系（見圖2），其3個軸矢量分別為

圖2 流場域時間步更新Fig.2 Time-step update of flow field

每個時間步后，局部坐標將跟隨結構網格發生位移，根據局部坐標在每個時間步前后的位移獲得變換矩陣T，則流場網格新的坐標系為

式中：［e*1e*2e*31］為更新后的齊次坐標；［e1e2e31］為更新前的齊次坐標。流體域通過時間步更新技術實現運動與重構。

2 研究對象與模型驗證

2.1 研究對象與計算工況

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

本文采用ALE流固耦合數值方法，將結構網格單元穿插于流場網格中，其流固耦合數值計算模型如圖5所示（規模為6WS×8WS×10WS，其中WS為傘衣弦長），其中結構附近流場進行加密處理，參數設置如表1所示。翼傘以0°迎角充氣，初始充氣速度為17m/s。本文基于ANSYS平臺計算，采用8核處理器，仿真0.6 s耗時約172 h。

圖3 沖壓式翼傘結構示意圖Fig.3 Schematic diagram of rammed parafoil structure

圖5 流固耦合數值計算模型Fig.5 Numerical calculation model of fluid-structure interaction

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

2.2 模型驗證

為驗證本文模型的合理性和準確性，采用文獻［22］的空投試驗數據進行驗證計算，同時進行流場網格的無關性驗證。本文采用3種不同網格密度的流場進行數值計算，對比結果如表2所示。

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

表2 網格無關性對比結果Table 2 Comparison results of grid independence

圖6 數值計算結果與空投試驗充滿外形對比Fig.6 Comparison of canopy shape between numerical calculation result and airdrop test

圖7 開傘動載隨時間的變化曲線Fig.7 Curves of opening load over time

3 研究結果

3.1 傘衣外形及流場變化

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

為進一步探究翼傘充氣過程中的結構及流場的非定常變化規律，選取a、b、c 3個截面（見圖9）進行研究，各個截面的結構及流場變化如圖10所示（各個時刻從上到下依次為a、b、c截面），為刻畫細節情況，各個時刻的圖形縮比尺寸不同。

圖8 充氣過程翼傘外形變化（范式等效應力云圖）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截面）：在翼傘初始折疊狀態時，a、b截面的速度場分布相似，前緣產生渦。隨翼傘逐漸充氣，上翼面的渦逐漸后移至脫體，前緣處不斷產生新的渦。在充氣過程中，a截面隨氣室充氣逐漸形成飽滿翼型，內腔壓力增加，傘衣剛度增加，類似剛性翼，外部流暢穩定。而b截面處于最外側，傘衣由于約束較少變形明顯，在翼傘充滿后翼型飽滿度不足，外部流場不穩定，傘衣變形狀態具有一定的隨機性。

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

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

3.2 氣室充氣規律

為進一步研究翼傘的各氣室充氣規律，將翼傘沿弦向分為前（前緣）、中（中部）、后（尾緣）三部分，沿展向將7個氣室逐次編號，兩側為氣室1和氣室7，中央氣室為4號。圖12為各氣室充滿時間和充滿寬度的變化規律。可以看出：

1）各氣室的充氣規律關于中央氣室具有對稱性，中央氣室充氣快且充滿外形飽滿，越靠近兩側充滿速率越慢，充氣效果越差。

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

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

3）由于尾緣上下翼面縫合，翼尖繞流導致傘衣變形，內外氣壓不穩定，兩側氣室的后部未完全充滿。

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

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

3.3 氣動特性分析

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

圖14 氣動特性參數的變化Fig.14 Changes in aerodynamic characteristic parameter

4 結 論

本文針對非定常充氣展開過程，數值模擬了翼傘的三維流固耦合動力學特性，得出：

1）翼傘充氣展開過程中穩定幅和前緣切口處產生應力集中現象，因此，在翼傘設計中應加強穩定幅與前緣切口處的材料強度。

2）翼傘充氣展開過程存在“翼尖上翹，中部凹陷”的翼傘尾流再附現象，造成傘衣內塌，導致系統穩定性受影響，因此在翼傘開傘過程中應加入引導傘以改善充氣狀況，設計合理的穩定幅可以減小翼尖渦的影響。

3）各氣室的充氣規律關于中央氣室具有對稱性，中央氣室充氣快且充滿外形飽滿，兩側充滿速率越慢，氣室飽滿度下降；翼傘尾部最先膨脹展開，前緣展開最慢，充滿后存在一定的回彈和氣室“鼓包”現象，充滿后展長小于設計值。

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