賈 賀，榮 偉，陳國良

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

0 引言

降落傘系統(tǒng)作為一種高效、可靠的回收裝置，在我國載人飛船、返回式衛(wèi)星的回收著陸以及美國火星探測器的著陸等過程中，都得到了廣泛的應用。充氣是降落傘整個工作過程中最為復雜的一個階段，也是研究人員迫切地想通過理論分析來解釋的[1-2]。目前，對于降落傘充氣性能（除最大開傘動載外）的了解或掌握主要是通過一些試驗結(jié)果的統(tǒng)計分析來獲取的，充氣性能與開傘條件之間的關系也是通過試驗結(jié)果統(tǒng)計、擬合出來的[3]。這使得充氣時間、傘衣投影面積變化的一些經(jīng)驗公式存在著一定的局限性，只能在某些情況下適用。

隨著“嫦娥一號”的成功發(fā)射，我國對于火星探測的研究即將進入實質(zhì)階段?；鹦谴髿鈱拥某煞?、物理性質(zhì)與地球的大氣層存在著較大的差別，因此，對于降落傘系統(tǒng)在火星大氣環(huán)境下的工作性能評定，不能像地球上回收航天器所用的降落傘系統(tǒng)那樣可以通過空投試驗來直接驗證，而需要通過理論分析或者數(shù)值仿真來進行。目前回收系統(tǒng)工作過程的數(shù)值仿真系統(tǒng)所用到降落傘開傘過程的一些參數(shù)，如充氣時間、傘衣投影面積的變化一般是采用經(jīng)驗公式來處理的。有關研究表明，對于火星這樣稀薄的大氣環(huán)境，這些常用的經(jīng)驗關系式將不再適用[4]。因此，如何通過數(shù)值仿真來獲取研究降落傘充氣過程性能的主要參數(shù)，已經(jīng)被提上了降落傘理論研究和工程實踐的日程。

進入21世紀以來，應用LS-DYNA有限元分析軟件來分析降落傘的工作過程（充氣階段、穩(wěn)定階段等）成為一個熱點。國外在這方面做了大量的工作，但是對降落傘充氣過程的仿真還不能完全地解決[5-8]。國內(nèi)使用 LS-DYNA軟件對降落傘充氣過程仿真研究尚屬空白。本文基于LS-DYNA軟件對數(shù)值模擬降落傘充氣過程的方法進行了初步研究。

1 降落傘的充氣過程

降落傘的充氣過程是指從傘系統(tǒng)全長地拉直到傘衣第一次充滿（傘衣投影直徑第一次能達到穩(wěn)定下降所需的傘衣投影直徑，不考慮過度充氣）為止的整個工作過程[1-2]，如圖1所示。

圖1 降落傘的充氣過程Fig. 1 Parachute inflation process

降落傘的充氣性能主要由臨界開傘速度、傘衣充滿條件、充氣時間、充氣距離、傘衣投影面積變化和最大開傘動載等參數(shù)來描述和表征[1,3]。本文基于LS-DYNA軟件對降落傘的充氣時間和傘衣投影面積變化的數(shù)值模擬分析，來研究降落傘的充氣過程。

2 LS-DYNA的流固耦合算法

LS-DYNA軟件是著名的顯式動力分析程序，具有強大精確的有限元仿真性能，擁有大量不同種類的單元模型、材料模型和算法選擇，能夠很方便地處理各種高度非線性問題。到目前為止，它是對降落傘充氣過程進行仿真的最優(yōu)工具之一。LS-DYNA 971是一個主要的程序版本。它增加了透氣性材料的流固耦合算法，這樣就可以通過控制流固耦合過程的各個參數(shù)，真實地模擬降落傘的充氣過程。

LS-DYNA與一般的CAE輔助分析程序操作過程相似，LS-DYNA的一個完整的顯式動力分析過程包括前處理、求解以及后處理3個基本操作環(huán)節(jié)，如圖2所示。

圖2 LS-DYNA的一般計算流程Fig. 2 Computation process of LS-DYNA

在解決流固耦合問題時，需要一種將Lagrange和 Euler的優(yōu)點結(jié)合起來的方法，即Arbitrary Lagrange-Euler，簡稱 ALE方法[9,10]。LS-DYNA程序中采用關鍵字*CONSTRAINED_ LAGRANGE_IN_SOLID實現(xiàn)固體與流體間的耦合。采用流固耦合算法模擬問題時，往往要對Lagrange算法中的固體結(jié)構(gòu)進行約束，將固體結(jié)構(gòu)的相關參數(shù)傳遞給流體單元。在LS-DYNA中數(shù)值模擬降落傘的充氣過程時，固體結(jié)構(gòu)的主要約束方法是罰函數(shù)耦合約束。對于罰函數(shù)約束而言，其原理是耦合系數(shù)追蹤Lagrange節(jié)點（結(jié)構(gòu)，即從物質(zhì)）和Euler流體物質(zhì)（主物質(zhì)）位置間的相對位移d，如圖3所示[11]。

圖3 罰函數(shù)耦合算法Fig. 3 Penalty function coupled method

3 LS-DYNA仿真模型的建立

3.1 降落傘仿真模型的建立

3.1.1 基本假設

對降落傘仿真模型作如下假設：

1）仿真充氣過程前模型為軸對稱的結(jié)構(gòu)；

2）模型的初始形狀為有一定的進氣孔尺寸，截面為梅花形狀的模型；

3）模型只考慮傘繩的連接，對于傘衣徑向帶、頂孔繩等不作具體處理。

3.1.2 基本參數(shù)設置

1）降落傘幾何模型

降落傘仿真模型選用美國的經(jīng)典傘型—C-9平面圓形傘，如圖4(a)所示，具體結(jié)構(gòu)參數(shù)如表1所示。

表1 C-9傘結(jié)構(gòu)參數(shù)表Table 1 Structure parameters of C-9 parachute

完整的傘衣和傘繩幾何模型如圖4(b)所示。

2）材料參數(shù)的選擇

降落傘仿真模型的傘衣材料選用美軍標的MIL-C-7020 III材料，假設該材料為各向同性材料，密度為533.77 kg/m3，彈性模量為0.4309 GPa，泊松比為0.14，傘衣厚度為0.000 1 m。傘繩材料的密度為462.00 kg/m3，彈性模量為97.0 GPa。假設降落傘傘繩的橫截面為圓形，面積為4.91×10-6m2。

3）網(wǎng)格的劃分

LS-DYNA的前處理軟件用于對降落傘仿真模型進行有限單元的網(wǎng)格劃分。傘衣單元選用薄殼單元，傘繩單元選用離散的梁/索單元，得到降落傘的有限元模型如圖4(c)所示。

圖4 降落傘模型Fig. 4 Parachute model

4）K文件的參數(shù)設置

降落傘仿真模型的傘衣單元為薄殼單元（*ELEMENT_SHELL_THICKNESS），單元算法的控制關鍵字為*SECTION_SHELL。傘衣材料屬性由*MAT_FABRIC關鍵字來定義。傘繩為離散的梁/索單元（*ELEMENT_BEAM），單元算法的控制關鍵字為*SECTION_BEAM。傘繩材料屬性由*MAT_CABLE_DISCRETE_BEAM來定義。

傘繩匯交點采用關鍵字*BOUNDARY_ SPC_NODE來進行約束，即對匯交點的3個方向平動自由度進行約束。

5）求解文件的生成

降落傘仿真模型建立完成，生成LS-DYNA軟件可以求解的文件格式，記為K文件(parachute.k)。

3.2 流場仿真模型的建立

3.2.1 基本假設

對流場仿真模型作如下假設：

1）流場模型模擬的為不可壓縮流場，符合N-S方程；

2）流場模型中流體的速度方向恒定，且速度大小不變；

3）流場模型外表面的邊界約束條件與風洞壁性質(zhì)一致。

3.2.2 基本參數(shù)設置

1）流場幾何模型

流場仿真模型形狀選用圓柱形流場，如圖5(a)所示。綜合文獻[12?14]，其基本的幾何尺寸確定如下：流場半徑 Df=12 m；流場長度 Lf=31 m；流場前端面距傘繩匯交點的長度Lq=6 m；未充氣流場后端面距傘衣頂端的長度Lh=12 m。

運用LS-DYNA的前處理軟件建立流場幾何模型如圖5(b)所示，流體速度為80 m/s。

2）材料參數(shù)的選擇

流場仿真模型的流體定義為理想氣體，滿足理想氣體狀態(tài)方程。氣體的密度為1.18 kg/m3，動態(tài)粘度系數(shù)為1.745 6×10-5，氣體的比熱容比為γ=1.4。

3）網(wǎng)格的劃分

同樣，在LS-DYNA前處理軟件中對流場的幾何模型進行有限單元的網(wǎng)格劃分。根據(jù)有限元網(wǎng)格劃分的基本原則，流場的有限元模型如圖 5(c)所示。

圖5 流場模型Fig. 5 Fluid model

4）K文件的參數(shù)設置

流場仿真模型的單元選擇體單元（*ELEMENT_ SOLID），流場仿真模型單元算法的關鍵字為*SECTION_SOLID。LS-DYNA要想真實地模擬降落傘的充氣過程，必須在流場前端增加一層流體壓力入口單元，主要作用是向流場中源源不斷地輸入恒壓、恒速的流體。流場單元算法參數(shù)的選擇對于解決流固耦合問題很關鍵，本文為流體單元選擇的是單點ALE方法的多物質(zhì)單元方程。對于已定義關鍵字*MAT_NULL的材料而言，必須包含狀態(tài)方程的定義，用來提供壓力行為應力組件，表示材料的體積變形與壓力之間的關系。這樣兩者合起來提供材料的整個應力張量。流體仿真模型的狀態(tài)方程用*EOS_LINEAR_POLYNOMIAL來定義。

流場仿真模型的流場入口單元區(qū)域用*BOUNDARY_PRESCRIBED_MOTION_SET來施加對單元的強制運動，使得流體壓力入口單元獲得恒定的速度。

5）求解文件的生成

流場仿真模型建立完成，最后生成LS-DYNA軟件可以求解的文件格式，記為K文件（fluid.k）。

4 LS-DYNA流固耦合控制參數(shù)

LS-DYNA軟件是一種比較理想的顯式求解工具，通過它的前處理器，生成一個可以遞交給LS-DYNA的文件，稱之為關鍵字文件或K文件，該文件是由一系列的關鍵字組成的[15]。LS-DYNA求解流固耦合問題，實質(zhì)上就是對流固耦合控制關鍵字進行參數(shù)設置的過程。

但是，降落傘充氣過程與一般的流固耦合問題不完全相同，如果按照一般的流固耦合問題處理，僅僅設置一般問題的基本控制關鍵字來仿充氣過程，可以發(fā)現(xiàn)傘衣沒有完全充開，最終的充氣形狀也很不理想，如圖 6所示?？梢?，使用LS-DYNA對充氣過程進行仿真時，不能僅僅采用一般問題的解決方法和思路，必須選擇適合充氣過程的控制關鍵字，并對關鍵字進行參數(shù)值的優(yōu)化處理。

圖6 基本控制參數(shù)下的傘衣最終充氣形狀Fig. 6 Final canopy shape with basic control parameters

通過大量的研究分析，得到在仿真降落傘充氣過程中，主要的流固耦合過程控制關鍵字及最優(yōu)化的控制參數(shù)值，其中主要的控制關鍵字如表2所示[16]?？梢姡x擇適合降落傘充氣過程的控制關鍵字及設置最優(yōu)化的參數(shù)值是仿真充氣過程的重點和難點，對仿真過程的實現(xiàn)起著決定性作用。

表2 流固耦合過程的主要控制關鍵字Table 2 Major keywords of Fluid-Structure Interaction(FSI) process

5 降落傘充氣過程的仿真結(jié)果及分析

5.1 降落傘充氣過程的數(shù)值模擬

圖7中(a) (i)表示的是在LS-DYNA中數(shù)值模擬降落傘充氣的全過程。

可以很清楚地看到，在充氣開始的時候，降落傘的充氣形狀變化比較慢。在0.02 s時，氣流使得降落傘的傘衣頂部開始發(fā)生形狀的變化；從0.02 0.10 s，傘衣頂部變化顯著，氣流將整個降落傘的頂部完全頂開。從0.10 0.16 s，氣流從傘衣頂部向傘衣周圍流動，使得傘衣的充氣形狀不斷發(fā)生變化；在0.16 s時，降落傘充氣達到最大外形形狀，即達到充滿狀態(tài)。從0.16 0.18 s降落傘的充氣形狀發(fā)生一次較為明顯的收縮。從0.18 0.26 s降落傘的充氣形狀會不斷發(fā)生微小的變化，傘衣產(chǎn)生較為明顯的“呼吸現(xiàn)象”。通過對傘衣充氣形狀的測量，發(fā)現(xiàn)在0.26 s時降落傘的充氣過程基本完成。

圖7 LS-DYNA數(shù)值模擬的降落傘充氣過程Fig. 7 Parachute inflation simulation with LS-DYNA

5.2 降落傘的充氣時間

5.2.1 充氣時間的理論求解

基于 LS-DYNA對降落傘充氣過程的數(shù)值模擬，其實就是對降落傘做類似于風洞試驗的研究。因此，數(shù)值模擬降落傘的充氣過程是一個無限質(zhì)量條件下的充氣過程?？赏ㄟ^理論方法求解無限質(zhì)量條件下C-9平面圓形傘的充氣時間tf。

[17，18]可知，對于無限質(zhì)量條件下的致密織物的降落傘充氣過程而言，降落傘的充氣時間指的是傘衣第一次達到最大外形狀態(tài)的時間。

設降落傘的名義直徑為D0，有效透氣量為C，開始充氣時的速度為 v0，無因次充氣時間為 T0，對于基本方程

而言，當考慮無限質(zhì)量條件充氣時，可以把速度 v看作一個常數(shù)，并且等于充氣開始時的速度v0，則基本方程可以直接積分。把上述方程寫成積分形式為

如果只考慮名義直徑為D0的平面圓形傘，并假設降落傘傘衣在完全充滿時呈半球狀，則有

式中Vmax為降落傘傘衣充滿時所包容的空氣體積。

對式(2)進行積分，并把式(3)代入，得到充氣時間的表達式

已知降落傘的名義直徑D0、有效透氣量C和降落傘開始充氣時的速度v0，就可以直接計算出在無限質(zhì)量條件下降落傘的充氣時間 tf。數(shù)值模擬充氣過程的C-9傘的名義直徑D0為8.5 m，開始充氣時的速度v0為80 m/s，故求解降落傘的充氣時間轉(zhuǎn)化為求解其有效透氣量的問題。由于C-9傘傘衣材料選用的MIL-C-7020 III型材料，由圖8所示的曲線，可以得到其壓差和有效透氣量的關系。

圖8 MIL-C-7020 III材料壓差與有效透氣量曲線Fig. 8 MIL-C-7020 III material pressure vs effective permeability

由于模擬無限質(zhì)量條件下的充氣過程可看作為風洞試驗，大氣密度ρ為1.18 kg/m3，則壓差為

由圖8可知，C-9傘傘衣的有效透氣量C約為0.09。由式(4)計算出其在無限質(zhì)量條件下充氣時間為tf=0.21 s。

5.2.2 充氣時間的數(shù)值模擬值與理論求解值的比較分析

數(shù)值模擬 C-9傘傘衣第一次達到最大外形狀態(tài)的充氣時間為0.16 s，這與理論求解的充氣時間0.21 s的偏差為23%，是在可接受的范圍之內(nèi)。并且數(shù)值模擬充氣過程時，可以求解得到傘衣最終達到穩(wěn)定狀態(tài)時的充氣時間為0.26 s。

5.3 降落傘的投影面積變化

5.3.1 投影面積變化的數(shù)值模擬

數(shù)值模擬的 C-9傘傘衣底邊投影面積的形狀變化如圖9所示。

圖9 傘衣底邊投影形狀的變化Fig. 9 Canopy project shape variations

5.3.2 數(shù)值模擬與試驗結(jié)果的比較分析

從充氣過程開始到傘衣最終達到穩(wěn)定狀態(tài)，數(shù)值模擬 C-9傘充氣過程的傘衣投影面積比隨充氣時間的變化情況。如圖10(a)所示，C-9傘充氣過程的傘衣投影面積比隨無因次充氣時間的變化如圖10(b)所示。

圖10 傘衣投影面積比的變化Fig. 10 Canopy project area ratio variations

由圖 10(b)中可以看出，傘衣投影面積隨著充氣時間的延長而不斷增加，在T=1時達到最大值，投影面積與傘衣名義面積之比為 0.42。而在 T=1之后傘衣有一個明顯的收縮過程。而從圖 10(a)可以更加清楚地看到，在tf=0.21 s時，傘衣的投影面積變化趨勢基本穩(wěn)定，雖然由于傘衣的“呼吸現(xiàn)象”會使得傘衣的投影面積產(chǎn)生一些微小的變化，但是在tf=0.21 0.26 s之間，可以認為傘衣的投影面積不再發(fā)生變化，即降落傘達到充滿穩(wěn)定狀態(tài)。

與 D0=8.5 m的平面圓形傘試驗結(jié)果[1]對比可以發(fā)現(xiàn)，數(shù)值模擬C-9平面圓形傘充氣過程的傘衣投影面積變化規(guī)律與試驗得到的變化規(guī)律基本一致。在達到充氣時間時，試驗的傘衣投影面積與名義面積比為0.4，數(shù)值模擬的傘衣底邊投影面積與名義面積比為0.42，面積比基本相等。可見，數(shù)值模擬降落傘充氣過程，得到的傘衣投影面積的變化情況與試驗得到的結(jié)果基本一致。

6 結(jié)束語

本文基于 LS-DYNA軟件對降落傘的充氣過程進行了數(shù)值模擬，得到了降落傘的充氣時間和投影面積變化的規(guī)律，并對這些主要參數(shù)進行了分析，驗證了使用 LS-DYNA軟件模擬充氣過程的可行性，為降落傘的設計及其性能分析提供了一種新的分析手段。由于流固耦合過程控制關鍵字的參數(shù)值還沒有完全達到最優(yōu)化，導致降落傘充氣時間的誤差較大，這也是下一步工作需解決的主要問題。

參考文獻（References）

[1] 王利榮. 降落傘理論與應用[M]. 北京宇航出版社, 1997

[2] 《降落傘技術(shù)導論》編寫組. 降落傘技術(shù)導論[M]. 北京: 國防工業(yè)出版社, 1977

[3] Ewing E G, Bixby H W, Knacke T W. 回收系統(tǒng)設計指南[M]. 北京: 航空工業(yè)出版社, 1988

[4] 榮偉. 火星探測器減速著陸技術(shù)研究[D]. 中國空間技術(shù)研究院博士學位論文. 北京: 中國空間技術(shù)研究院, 2008

[5] Taylor A P. An investigation of the apparent mass of parachutes under post inflation dynamic loading through the use of fluid structure interaction simulations, AIAA 2003-2104[R], 2003

[6] Benjamin A T, Anthony P T. The use of LS-DYNA to simulate the inflation of a parachute canopy, AIAA 2005-1608[R], 2005

[7] Benjamin A T, Anthony P T, Jean C B, et al. The use of LS-DYNA to assess the performance of airborne system North America candidate ATPS main parachutes, AIAA 2005-1609[R], 2005

[8] Benjamin A T. The application of a new material porosity algorithm for parachute analysis[C]//9thInternational LS-DYNA Users Conference, 2006

[9] 李裕春, 等. ANSYS 11.0/LS-DYNA 基礎理論與工程實踐[M]. 北京: 中國水利水電出版社, 2008

[10] 趙海鷗. LS-DYNA 動力分析指南[M]. 北京: 兵器工業(yè)出版社, 2003

[11] John O H. Theoretical Manual[M]. Livermore Software Technology Corporation(LSTC), 1998

[12] Lingard J S, Darley M G. Simulation of parachute fluid structure interaction in supersonic flow, AIAA 2005-1607[R], 2005

[13] Lingard J S, Darley M G, Underwood J C. Simulation of Mars supersonic parachute performance and dynamics, AIAA 2007-2507[R], 2007

[14] Christine Espinosa, Yves de Lassat de Pressigny, Pascal Bordenave, et al. Fluid-structure interaction simulation of parachute dynamic behaviour, AIAA 2007-2510[R], 2007

[15] Livermore Software Technology Corporation(LSTC). LS-DYNA keyword user’s manual[G]. Version 971, 2007

[16] 賈賀. 基于LS-DYNA的降落傘充氣過程研究[D]. 中國空間技術(shù)研究院碩士學位論文. 北京: 中國空間技術(shù)研究院, 2009

[17] 彭勇, 宋旭民. 降落傘充氣時間的計算方法[J]. 航天返回與遙感, 2004, 25(1): 17-20

[18] 返回技術(shù)翻譯組. 氣動力減速器原理及設計手冊(上冊)[M]. 北京: 國防工業(yè)出版社, 1976