十字形傘開傘充氣過程數(shù)值仿真研究

謝淮 劉宇 王臻 昌飛

十字形傘開傘充氣過程數(shù)值仿真研究

謝淮 劉宇 王臻 昌飛

（北京空間機電研究所，北京 100094）

十字形降落傘是航空航天減速領域常用的傘型，開傘充氣過程是其工作過程中較為復雜的一個階段。為了研究十字形傘開傘充氣過程中的動力學特性，文章基于任意拉格朗日-歐拉方法對十字形傘開傘過程進行流固耦合仿真計算，并將計算結果與風洞試驗結果進行對比。對比發(fā)現(xiàn)：仿真計算得到的十字形傘穩(wěn)態(tài)外形與試驗結果一致，且仿真獲得的降落傘氣動載荷歷程曲線與試驗結果基本相符，載荷呈現(xiàn)出隨著開傘過程逐漸增大，在開傘結束達到最大值后略有降低并逐漸穩(wěn)定的趨勢。仿真結果還表明，在十字形傘開傘過程中，傘衣的最大應力點位于傘臂的中心區(qū)域，充氣順序的先后會影響傘臂的應力分布，降落傘穩(wěn)定后傘衣應力呈對稱分布。文章采用的仿真方法能較好地模擬十字形傘開傘充氣的動力學過程，得到的傘衣應力分布特點及影響因素可為十字形傘的設計與優(yōu)化提供參考。

十字形降落傘 計算流固耦合動力學 開傘 充氣

0 引言

降落傘是目前航空航天領域應用最廣泛的減速手段[1-2]，根據(jù)傘衣結構形式的不同，可以分為平面圓形傘、方形傘、帶條傘、環(huán)帆傘以及十字形降落傘等。相較于其他傘形，十字形降落傘具有工藝簡單、穩(wěn)定性好、開傘動載小等優(yōu)點，其獨特的結構使得它在高空高速降落和復雜的氣象條件下表現(xiàn)出色，成為眾多型號任務的理想選擇[3-4]，被廣泛用于軍事空投、緊急救援、探測器著陸以及科學實驗等領域。

由于十字形傘具有獨特的設計特點和廣泛的應用背景，對其力學特性進行深入研究具有重要意義[5]。在降落傘方案設計中，充氣過程中的力學特性[6]，包括充氣時間、氣動力變化以及傘衣載荷變化是需要考慮的關鍵因素[7]。相較于高成本的試驗手段[8-10]，仿真計算方法是目前研究降落傘開傘過程動力學的有效手段[11]。

流固耦合是降落傘開傘過程仿真的核心問題，因為降落傘傘衣材料為柔性結構，同時具有透氣性，這使得降落傘開傘過程流固耦合處理及其困難[12]。目前應用于降落傘開傘計算的流固耦合仿真方法有CFD-MSD（Computational Fluid Dynamics-Mass Spring Damper）方法[13]、DSD/SST（Deforming-Spatial- Domain/Stabilized Space-Time）方法[14]及浸入邊界法（Immersed Boundary Method，IBM）[15]等，但是這些方法均具有局限性：CFD-MSD方法側重于傘衣外形的計算，無法得到結構應力信息；DSD/SST方法側重于外部流場的數(shù)值計算，難以對充氣階段進行研究；IBM側重于固體的響應，對流場的細節(jié)計算不夠精確。近些年來，商業(yè)軟件LS-DYNA的流固耦合方法可以克服以上方法的缺點，滿足降落傘仿真的需要，因此被廣泛應用于亞聲速降落傘、超聲速降落傘等降落傘仿真領域[16-19]。在LS-DYNA的流固耦合（Fluid-Structure Interaction，F(xiàn)SI）方法中，流場單元通過任意拉格朗日-歐拉（Arbitrary Lagrange-Euler，ALE）方法進行求解，降落傘單元通過拉格朗日法進行求解，最后通過罰函數(shù)方法對流場及傘衣之間的相互作用力進行耦合計算。

本文基于ALE/FSI方法對十字形傘開傘過程進行仿真計算，并將計算結果與風洞試驗結果進行對比，驗證仿真模型的正確性，在此基礎上分析十字形傘開傘過程中載荷變化以及傘衣內(nèi)的應力分布變化的特點，從而能夠更深入地理解十字形傘開傘過程動力學特性，為其設計和應用提供更準確的指導。

1 研究方法

1.1 ALE方法

式中L為拉格朗日描述下的質點坐標；A為ALE描述下的質點坐標；和分別為質點運動速度和相應的ALE坐標點運動速度。由式（1）可以看出，當=時，有L=A，公式轉變?yōu)槔窭嗜招问?；?0時，右邊第一項為0，公式轉變?yōu)闅W拉形式。

基于ALE方法的N-S方程可表示為

通過算子分裂的方法[23]，式（2）可以分為兩步求解：拉格朗日步和輸運步。首先通過一個顯式有限元方法實施拉格朗日步，在該過程中網(wǎng)格點隨流體質點運動，即在方程中令=。拉格朗日步的方程為

質量方程也可以通過拉格朗日法求解，但是在此可使用更簡便和精確的方法

式中A,i表示矢量的第個分量，L,j表示矢量的第個分量，由于A和L均為三維矢量，故=1,2,3,=1,2,3。

接下來實施輸運步，網(wǎng)格進行運動，質量、動量和能量均因網(wǎng)格運動而產(chǎn)生通量。在本文中，考慮網(wǎng)格運動回初始構型位置，即完整的ALE過程完成后網(wǎng)格保持不變，因此輸運步的方程可表示為

1.2 結構控制方程

降落傘結構為柔性連續(xù)介質，其控制方程為

1.3 罰函數(shù)方法

在流固耦合問題中，流體與固體之間的相互作用是通過模擬兩者之間的相互影響來實現(xiàn)的。其中一種常用的方法是罰函數(shù)方法。在罰函數(shù)法中，為了模擬接觸力的相互作用，在接觸節(jié)點處使用一個類似于彈簧的力模型。在該力模型中，假設存在兩個接觸節(jié)點，引入一個假想的“彈簧”，它對兩個節(jié)點施加的力大小相等，方向相反。這個“彈簧”的作用類似于一個力的傳遞媒介，使得節(jié)點之間的接觸力可以相互傳遞。當一個節(jié)點施加力到彈簧上時，彈簧會將等大、相反方向的力傳遞給另一節(jié)點。這種模擬接觸力的處理方法實施簡單，可以使得節(jié)點之間的接觸行為在流固耦合問題的模擬中得到有效的考慮，同時保證了計算過程中動量及能量的守恒性。

在本文中用罰函數(shù)方法來解決兩個問題：傘衣表面的相互接觸（自接觸）和氣體與傘衣表面的接觸（流固耦合）。

1.3.1 傘衣自接觸

當未充滿狀態(tài)下傘衣受來流與傘繩共同作用發(fā)生變形時，傘衣單元之間會發(fā)生相互接觸，本文使用單面接觸算法模擬這一現(xiàn)象。以二維情況為例，定義穿透距離1為一個單元節(jié)點（Slave）與面（Master）之間的垂直距離，則

式中為面上一點到節(jié)點的矢徑；為面法向量。由1的定義可知，當1= 0時，節(jié)點與面發(fā)生接觸，當1<0時，節(jié)點與面發(fā)生穿透。

當穿透發(fā)生時，定義約束力1為

式中1為剛度系數(shù)。

1.3.2 流固耦合

在流固耦合模型中，流體單元為三維ALE實體單元，傘衣單元為二維拉格朗日殼單元，流體單元與傘衣單元之間同樣通過罰函數(shù)法耦合來傳遞作用力。當拉格朗日單元與ALE單元之間相對運動時，將拉格朗日單元節(jié)點視為從節(jié)點，不考慮網(wǎng)格運動的ALE描述與拉格朗日描述的單元節(jié)點坐標相等，考慮網(wǎng)格運動后可進行計算得到主節(jié)點，定義穿透位移2為

式中m，s分別為主節(jié)點及從節(jié)點運動速度。耦合的作用力定義為

式中2為耦合剛度系數(shù)。

1.4 傘衣透氣性方程

織物透氣量對傘衣充氣過程影響很大，而在充氣過程中隨著傘衣的內(nèi)外壓力差變化，織物透氣量也產(chǎn)生變化。傘衣的透氣速率與傘衣兩側壓差的關系可由經(jīng)驗公式表示為[24]

2 數(shù)值計算模型

2.1 研究對象

根據(jù)試驗任務，本文仿真計算所使用十字形傘外形尺寸如圖1所示。降落傘的材料選用織物透氣性參數(shù)為1=35.4和2=5.29×106/m的錦絲綢。

圖1 十字形傘外形示意

為降低模型復雜度，作以下基本假設：不考慮前體的影響，連接繩一端為傘繩集束點，另一端固定于流場中；傘衣部分只考慮加強帶，其余部件忽略；流場為亞聲速不可壓縮流場，且為理想氣體的穩(wěn)態(tài)定常流動；不考慮傘繩與流場之間的相互耦合；傘衣為各向同性材料。降落傘部件結構參數(shù)見表1。

表1 降落傘部件結構參數(shù)

Tab.1 Structural parameters of parachute components

2.2 網(wǎng)格劃分

2.2.1 降落傘網(wǎng)格

本文通過動力學方法建立降落傘初始模型，如圖2所示。首先建立十字形傘平面模型，傘衣為四邊形單元，傘繩為線單元（見圖2（a））；然后將傘衣中點固定，對連接繩一端施加遠離傘衣的位移載荷，在傘繩作用下，降落傘傘衣將以十字形折疊（見圖2（b））。為了保證展開后傘繩長度一致，初始時刻部分傘繩呈折線狀態(tài)。

圖2 十字形傘網(wǎng)格模型

2.2.2 流場網(wǎng)格

為了與試驗風洞尺寸一致，從而方便與試驗結果進行對比，將流場區(qū)域設置為六面體，尺寸為18 m×8 m×6 m。流場網(wǎng)格為六面體單元，在降落傘處對網(wǎng)格進行了十字形加密，總單元數(shù)為1.7×106個，如圖3所示。

圖3 流場網(wǎng)格示意

3 開傘過程仿真

3.1 工況設置

在建立數(shù)值計算模型的基礎上，為了使仿真計算結果更加接近十字形傘實際工作過程中的狀態(tài)，同時方便與已有的風洞試驗結果進行對比，本文設置以下計算工況：來流速度40 m/s，氣體密度1.35 kg/m3，溫度258 K；流場除入口外，其余壁面均為無反射邊界。此工況下獲得的仿真結果具有一定的代表性，可以反映十字形傘開傘充氣過程的動力學規(guī)律。

3.2 仿真結果分析

3.2.1 穩(wěn)態(tài)外形

十字形傘充氣過程結束后進入穩(wěn)態(tài)階段，截取穩(wěn)態(tài)階段傘衣狀態(tài)并與風洞試驗照片進行對比，如圖4所示，可見穩(wěn)態(tài)階段仿真所得十字形傘外形與試驗所得基本一致，說明本文采用的數(shù)學模型及數(shù)值方法能夠較好地模擬十字形傘實際工作狀況，結果具有可靠性。

圖4 仿真與風洞試驗的穩(wěn)態(tài)傘衣外形對比

圖5 仿真與風洞試驗的氣動載荷曲線對比

3.2.2 載荷歷程

提取連接繩末端所受拉力，以此表征降落傘所受氣動載荷。圖5為十字形傘開傘過程載荷變化仿真結果與試驗結果的對比，由圖可知，氣動載荷在開傘過程中逐漸增大，達到峰值后略有降低，隨后趨于穩(wěn)定，載荷值在一定范圍內(nèi)波動，仿真結果與試驗結果的趨勢一致。

定義開傘載荷達到最大值的時刻為開傘過程結束時刻，由圖5可得開傘時間仿真結果約為0.66 s，試驗結果約為0.72 s，誤差為8.3%；降落傘最大載荷仿真結果為5 070 N，試驗結果為4 466 N，誤差為13.5%；降落傘充氣穩(wěn)定后阻力隨時間的平均值仿真結果約為4 868 N，試驗結果約為4 278 N，誤差為13.8%。仿真獲得的3項結果均與試驗值存在一定誤差，原因可能是仿真中沒有考慮湍流的影響，誤差均在工程可接受的合理范圍內(nèi)。

3.2.3 傘衣應力

為了分析開傘充氣全過程傘衣內(nèi)的應力變化，在圖6中分別展示了充氣過程中（=0.47 s）、載荷達到最大值時刻（=0.66 s）及降落傘穩(wěn)定階段（=2.0 s）三個不同時刻傘衣的有效應力分布情況。圖中紅色區(qū)域即為應力集中區(qū)域，可見無論是充氣階段還是穩(wěn)定階段，應力最大點都位于傘衣臂中心附近區(qū)域。

圖6 不同時刻傘衣應力分布

圖7為十字形傘4個傘臂中心點（1～4）有效應力隨時間變化的曲線。由圖可知，在0.66 s之前即充氣過程中，由于4個傘衣臂張開先后順序不同，4個傘衣臂中心點有效應力雖然均為上升趨勢，但并不是同步上升，結合圖6（a）可知各傘衣臂的應力分布也不同；在0.66 s時降落傘充氣完成，載荷達到最大值，4個傘衣臂中心點的有效應力值幾乎同時上升到最大值；0.66 s以后即降落傘穩(wěn)定后，各傘衣臂中心點的有效應力值在達到最大值后均略有降低并趨于穩(wěn)定，結合圖6（c）可知傘衣應力分布趨于穩(wěn)定，且呈現(xiàn)對稱分布狀態(tài)。

圖7 傘臂中心點有效應力變化歷程

4 結束語

本文基于流固耦合方法對十字形降落傘開傘充氣過程進行了仿真。在本文的計算模型中，通過ALE方法求解流場N-S方程，通過拉格朗日方法求解降落傘結構變形，通過罰函數(shù)方法處理流場與降落傘之間的流固耦合及傘衣之間接觸。

仿真結果與風洞試驗結果進行對比分析，發(fā)現(xiàn)十字形傘充氣后的穩(wěn)態(tài)外形仿真結果與試驗基本一致，且開傘過程中載荷變化曲線與試驗結果基本吻合，證明了仿真方法的適用性與準確性。

仿真結果表明，十字形傘傘衣的最大應力點總是位于傘臂的中心位置，但不同傘臂的應力分布受充氣先后順序的影響；在降落傘達到穩(wěn)定狀態(tài)后，傘衣的應力分布呈現(xiàn)對稱分布。

本文的仿真計算方法可作為試驗手段的有效補充。通過對十字形傘傘衣應力分布情況的分析可以認識到開傘過程中存在的應力分布不均現(xiàn)象，有助于避免降落傘因應力集中而發(fā)生失穩(wěn)或破壞，可為十字形傘優(yōu)化改進提供參考。但是本文的模型未考慮湍流等因素的影響，仿真結果仍然存在一定誤差，未來可在本文工作的基礎上優(yōu)化計算模型，以提高仿真的精度。

[1] 王利榮. 降落傘理論及應用[M]. 北京:宇航出版社, 1997: 102-104. WANG Lirong. Parachute Theory and Application[M]. Beijing: China Astronautic Publishing House, 1997:102-104. (in Chinese)

[2] 賈華明, 楊霞, 李少騰, 等. 環(huán)帆傘技術與發(fā)展綜述[J]. 航天返回與遙感, 2021, 42(3): 41-51. JIA Huaming, YANG Xia, LI Shaoteng, et al. Overview of the Technology and Development of Ringsail Parachute[J]. Spacecraft Recovery & Remote Sensing, 2021, 42(3): 41-51. (in Chinese)

[3] 韓雅慧, 楊春信, 肖華軍, 等. 十字形傘的實驗研究進展及展望[J]. 兵工自動化, 2013, 32(3): 3-8, 20. HAN Yahui, YANG Chunxin, XIAO Huajun, et al. Development History and Expectation of Cross Parachute[J].Ordnance Industry Automation, 2013, 32(3): 3-8, 20. (in Chinese)

[4] 喻東明, 馬嘯民, 楊品. 亞聲速條件下十字形傘充滿時間系數(shù)的解算方法及仿真驗證[J]. 南京航空航天大學學報, 2021, 53(2): 177-181. YU Dongming, MA Xiaomin, YANG Pin. Solution of Inflation Coeffcient of Cruciform Parachute Under Subsonic Conditions and Simulation Verification[J]. Journal of Nanjing University of Aeronautics & Astronautics, 2021, 53(2): 177-181. (in Chinese)

[5] 賈華明, 李健. 收口十字形降落傘充氣過程動力學建模與仿真[J]. 航天返回與遙感, 2012, 33(5): 16-23. JIA Huaming, LI Jian.Modeling and Simulation of Inflation Process Dynamics of Reefed Cruciform Parachute [J]. Spacecraft Recovery & Remote Sensing, 2012, 33(5): 16-23. (in Chinese)

[6] 賈賀, 包進進, 榮偉. 設計參數(shù)及大氣參數(shù)對降落傘充氣性能的影響[J]. 航天返回與遙感, 2020, 41(3): 28-36. JIA He, BAO Jinjin, RONG Wei. The Design and Atmospheric Parameters Influences on Parachute Inflation Performance[J]. Spacecraft Recovery & Remote Sensing, 2020, 41(3): 28-36. (in Chinese)

[7] 榮偉. 航天器進入下降與著陸技術[M]. 北京: 北京理工大學出版社, 2018 RONG Wei. Spacecraft Entry, Decent and Landing Technology[M]. Beijing: Science Press, 2018. (in Chinese)

[8] 王立武, 房冠輝, 李健, 等. 降落傘超聲速低動壓高空開傘試驗[J]. 航天返回與遙感, 2020, 41(3): 1-9. WANG Liwu, FANG Guanhui, LI Jian, et al. The Experiments of Parachute Deployment at Supersonic and Low Dynamic Pressure Condition[J]. Spacecraft Recovery & Remote Sensing, 2020, 41(3): 1-9. (in Chinese)

[9] 余莉, 明曉, 胡斌. 降落傘開傘過程的試驗研究[J]. 南京航空航天大學學報, 2006, 38(2): 176-180. YU Li, MING Xiao, HU Bin.Experimental Investigation in Parachute Opening Process[J]. Journal of Nanjing University of Aeronautics & Astronautics, 2006, 38(2):176-180. (in Chinese)

[10] 余莉, 騰海山, 明曉. 利用七孔探針對降落傘流場的試驗測量研究[J]. 南京航空航天大學學報, 2007, 39(5): 65-70. YU Li, TENG Haishan, MING Xiao. Seven-hole Probe Measurement on Parachute Flow Field[J]. Journal of Nanjing University of Aeronautics & Astronautics,2007, 39(5): 65-70. (in Chinese)

[11] CAO Y, NIE S, WU Z. Numerical Simulation of Parachute Inflation: A Methodological Review[J]. Proceedings of the Institution of Mechanical Engineers, Part G: Journal of Aerospace Engineering, 2019, 233(2): 736-766.

[12] 包文龍, 賈賀, 薛曉鵬, 等. 開“窗”結構對環(huán)帆傘穩(wěn)態(tài)氣動性能的影響[J]. 航天返回與遙感, 2022, 43(1): 1-11. BAO Wenlong, JIA He, XUE Xiaopeng, et al. The Influence of the‘Windows’Structure on the Steady-state Aerodynamic Performance of Ringsail Parachute[J]. Spacecraft Recovery & Remote Sensing, 2022, 43(1): 1-11. (in Chinese)

[13] STEEVES E C, BENNEY R J, STEIN K R. A Computational Model that Couples Aerodynamic and Structural Dynamic Behavior of Parachutes during the Opening Process: NASA-ADA264115[R]. Washingtong DC: NASA, 1993.

[14] STEIN K, BENNEY R, TEZDUYAR T, et al.Fluid-structure Interactions of a Cross Parachute: Numerical Simulation[J]. Computer Methods in Applied Mechanics & Engineering, 2001, 191(6/7):673-687.

[15] KIM Y, PESKIN C S. 3-D Parachute Simulation by the Immersed Boundary Method[J]. Computers & Fluids, 2009, 38(6): 1080-1090.

[16] TUTT B, TAYLOR A. The Use of LS-DYNA to Simulate the Inflation of a Parachute Canopy[C]// 18th AIAA Aerodynamic Decelerator Systems Technology Conference and Seminar, May 23-26, 2005, Munich, Germany.

[17] GAO X, ZHANG Q, TANG Q. Numerical Modelling of Mars Supersonic Disk-gap-band Parachute Inflation[J]. Advances in Space Research, 2016, 57(11): 2259-2272.

[18] COQUET Y, BORDENAVE P, CAPMAS G, et al. Improvements in Fluid Structure Interaction Simulations of Parachutes Using LS-Dyna?[C]// 21st AIAA Aerodynamic Decelerator Systems Technology Conference and Seminar, May 23-26, 2011, Dublin, Ireland.

[19] 賈賀, 榮偉, 陳國良. 基于LS-DYNA軟件的降落傘充氣過程仿真研究[J]. 航天器環(huán)境工程, 2010, 27(3): 367-373. JIA He, RONG Wei, CHEN Guoliang. The Simulations of Parachute Inflation Process by Using LS-DYNA[J]. Spacecraft Environment Engineering, 2010, 27(3): 367-373. (in Chinese)

[20] 高興龍. 基于ALE方法的開縫降落傘充氣過程研究[J]. 航天返回與遙感, 2013, 34(1): 13-20. GAO Xinglong. A Study on the Slots-parachute Inflation Based on ALE Method[J].Spacecraft Recovery & Remote Sensing, 2013, 34(1), 13-20. (in Chinese)

[21] ZHAO M, ZHANG S, ZHANG Z, et al. Modeling and Simulation of Supersonic Parachute Inflation in Mars Environment[C]// Proceedings of the 21st International Conference on MMESE, October 23-25, 2021, Beijing, China. Singapore: Springer, 2022: 445-451.

[22] HIRT C W, AMSDEN A A, COOK J L. An Arbitrary Lagrangian–Eulerian Computing Method for All Flow Speeds[J]. Journal of Computational Physics, 1997, 135(2): 203-216.

[23] BENSON D J. Computational Methods in Lagrangian and Eulerian Hydrocodes[J]. Computer Methods in Applied Mechanics & Engineering, 1992, 99(2/3): 235-394.

[24] LINGARD J, UNDERWOOD J. The Effects of Low Density Atmospheres on the Aerodynamic Coefficients of Parachutes[C]// 13th Aerodynamic Decelerator Systems Technology Conference, May15-18, 1995, Clearwater Beach, FL, USA.

Numerical Simulation Study of Cruciform Parachute Deployment and Inflation Process

XIE Huai LIU Yu WANG Zhen CHANG Fei

（Beijing Institute of Space Mechanics & Electricity, Beijing 100094, China）

The cruciform parachute is a commonly used parachute design in the field of aerospace deceleration，and the deployment and inflation stage is a complex stage of its operation. In order to study the dynamics characteristics of cruciform parachute deployment and inflation process, this paper employs the arbitrary Lagrangian-Eulerian method to carry out the fluid-structure interactionsimulation calculations on the deployment process of a cruciform parachute. The computed results are compared with wind tunnel test results. The comparison shows that the simulated steady-state shape of the parachute is consistent with the experimental results, and the simulated aerodynamic load history curve of the parachute is largely in agreement with the test results. The load exhibits a gradual increase during the deployment process, reaching its maximum value after the completion of the deployment, followed by a slight decrease and gradual stabilization. The simulation results also show that, during the cruciform parachute deployment process, the maximum stress point of the canopy is located in the central region of the canopy arms, and the inflation sequence affects the stress distribution along the arms. After the parachute stabilizes, the canopy stress exhibits a symmetrical distribution. The simulation method employed in this study effectively simulates the dynamic process of cruciform parachute deployment and inflation, and the characteristics of canopy stress distribution and influencing factors can provide valuable insights for the design and optimization of cruciform parachutes.

cruciform parachute; fluid-structure interaction (FSI); deployment; inflation

V445

A

1009-8518(2023)03-0032-09

10.3969/j.issn.1009-8518.2023.03.004

謝淮，男，1976年生，2006年獲國防科技大學管理科學與工程專業(yè)碩士學位，高級工程師。主要研究方向為航天器進入、減速與著陸技術。E-mail：xh7603@hotmail.com。

2023-01-06

工業(yè)和信息化部重點實驗室開放基金（KLAECLS-E-202004）

謝淮, 劉宇, 王臻, 等. 十字形傘開傘充氣過程數(shù)值仿真研究[J]. 航天返回與遙感, 2023, 44(3): 32-40.

XIE Huai, LIU Yu, WANG Zhen, et al. Numerical Simulation Study of Cruciform Parachute Deployment and Inflation Process[J]. Spacecraft Recovery & Remote Sensing, 2023, 44(3): 32-40. (in Chinese)

（編輯：陳艷霞）