包文龍，賈賀，2，*，薛曉鵬，黃雪姣，高樹義，榮偉，王奇，吳壯志

1．北京空間機(jī)電研究所，北京 100094

2．南京航空航天大學(xué) 航空學(xué)院，南京 210016

3．中南大學(xué) 航空航天學(xué)院，長沙 410083

4．北京航空航天大學(xué) 計(jì)算機(jī)學(xué)院，北京 100191

環(huán)帆傘因開傘性能優(yōu)越、抗破壞能力強(qiáng)、工作可靠、阻力系數(shù)大和穩(wěn)定性好，被作為核心減速裝置廣泛應(yīng)用于各種航天器的回收著陸系統(tǒng)［1-3］。作為一種依靠其結(jié)構(gòu)形狀和性能特征來改變載荷體運(yùn)動(dòng)狀態(tài)的氣動(dòng)減速裝置，其氣動(dòng)特性的好壞尤為重要。已有試驗(yàn)表明，適當(dāng)改變傘衣的透氣性可以使阻力效率和開傘可靠性，以及穩(wěn)定性和開傘動(dòng)載之間達(dá)到一個(gè)微妙的平衡［4］。

近年來，隨著計(jì)算機(jī)技術(shù)的發(fā)展，國內(nèi)外已有諸多的研究團(tuán)隊(duì)利用數(shù)值模擬技術(shù)對環(huán)帆傘透氣性進(jìn)行研究。Tezduyar 等［5］將環(huán)帆傘傘衣劃分為12 個(gè)不同位置的同心圓區(qū)域，利用局部變化的織物透氣性來模擬整個(gè)環(huán)帆傘的透氣性，從而實(shí)現(xiàn)織物透氣性與結(jié)構(gòu)透氣性統(tǒng)一，該方法有助于減少環(huán)帆傘這一復(fù)雜傘型數(shù)值計(jì)算所需的網(wǎng) 格 數(shù) 量。 Takizawa 等［6］應(yīng) 用FSI（Fluid-Structure Interation）方法分析了3 種環(huán)帆傘設(shè)計(jì)構(gòu)型的性能，結(jié)果表明移除靠近傘衣底部的帆能有 效 提 高 靜 穩(wěn) 定 性。Greathouse 和Schwing［7］基于CFD（Computational Fluid Dynamics）方法對不同結(jié)構(gòu)的環(huán)帆傘開展性能分析，以研究結(jié)構(gòu)透氣性對環(huán)帆傘阻力特性和靜穩(wěn)定性的影響，并提出以最小阻力損失最大化靜穩(wěn)定性的設(shè)計(jì)方案。甘小嬌等［8-11］探究了部分結(jié)構(gòu)參數(shù)及織物透氣性對環(huán)帆傘氣動(dòng)性能的影響，但均選擇相對簡單的穩(wěn)降階段。部分學(xué)者針對環(huán)帆傘充氣過程進(jìn)行了數(shù)值模擬，但未能研究透氣性對開傘性能的影響［12-13］。

環(huán)帆傘的設(shè)計(jì)改進(jìn)不可避免地涉及透氣性，但當(dāng)前透氣性對其開傘過程的影響研究還不夠深入。另外，隨著中國航天技術(shù)的發(fā)展及空間站建設(shè)，返回艙的重量及體積將有所增加。而目前采用的單傘系統(tǒng)受工藝、成本、可靠性等因素影響已無法滿足工程需求，群傘系統(tǒng)獲得關(guān)注并已在新一代載人飛船試驗(yàn)船成功應(yīng)用，但也帶來開傘不同步和傘間干擾等問題。

為解決上述問題，國外經(jīng)多次空投及風(fēng)洞試驗(yàn)設(shè)計(jì)修改，采用了開“窗”結(jié)構(gòu)設(shè)計(jì)，發(fā)現(xiàn)該結(jié)構(gòu)有利于提高群傘系統(tǒng)的同步性和穩(wěn)定性［14-15］，圖1 為NASA 獵戶座飛船主傘采用的開“窗”結(jié)構(gòu)。中國目前對此結(jié)構(gòu)尚未應(yīng)用也無相關(guān)公開仿真研究成果。

圖1 獵戶座飛船主傘的開“窗”結(jié)構(gòu)Fig.1 ‘Windows’ structure of NASA Orion main parachute

同時(shí)，由于群傘間復(fù)雜氣動(dòng)干擾的存在，對其開傘過程數(shù)值仿真難度較大。而群傘中各個(gè)組成傘結(jié)構(gòu)相同，且很大程度上決定了群傘系統(tǒng)整體性能的優(yōu)劣。故本文采用流固耦合方法，對不同開“窗”結(jié)構(gòu)的環(huán)帆傘單傘在無限質(zhì)量情況下的開傘過程進(jìn)行數(shù)值模擬，進(jìn)而分析開“窗”結(jié)構(gòu)對其氣動(dòng)特性的影響，為后續(xù)從開“窗”角度設(shè)計(jì)優(yōu)化環(huán)帆傘，從而提高群傘系統(tǒng)性能提供參考。

1 研究對象

本文以某環(huán)帆傘為研究對象，開展開傘過程數(shù)值計(jì)算。其平面結(jié)構(gòu)及具體尺寸如圖2 和表1所示。同時(shí)，為方便后續(xù)分析，圖2（a）對環(huán)1～4和帆1～6 進(jìn)行了序號標(biāo)注。

圖2 環(huán)帆傘結(jié)構(gòu)Fig.2 Flat structure of ringsail canopy

表1 環(huán)帆傘結(jié)構(gòu)尺寸Table 1 Structural size of ringsail parachute

由于加強(qiáng)帶等細(xì)化結(jié)構(gòu)會(huì)影響流固耦合的仿真結(jié)果［16］，故本文的環(huán)帆傘模型由環(huán)、帆、傘繩、頂孔繩、徑向帶、緯向帶等結(jié)構(gòu)構(gòu)成，確保仿真模型與真實(shí)物理更為貼近。傘衣幅結(jié)構(gòu)和整體繩帶結(jié)構(gòu)如圖3、圖4 所示。

圖3 環(huán)帆傘傘衣幅結(jié)構(gòu)Fig.3 Structure of ringsail canopy gore

圖4 整體繩帶結(jié)構(gòu)Fig.4 Integral rope structure

環(huán)帆傘的環(huán)片和帆片為2 種不同的材料，經(jīng)測得厚度分別為0.15，0.10 mm。采用Ergun 公式對2 種材料的透氣性進(jìn)行描述［17-18］，可得

式中：a、b分別為描述織物透氣性的黏性系數(shù)和慣性系數(shù)，可經(jīng)試驗(yàn)獲得；e為傘衣厚度；vq為透過傘衣織物的氣流速度；Δp為傘衣織物的試驗(yàn)壓差。經(jīng)透氣量與壓差曲線雙點(diǎn)估計(jì)反算可得開傘過程2 種織物材料的黏性系數(shù)a和慣性系數(shù)b。

環(huán)片：

帆片：

通常認(rèn)為在最大投影直徑位置以下進(jìn)行開“窗”對阻力特性影響較小，該傘空投試驗(yàn)中的投影直徑Dt與結(jié)構(gòu)直徑Dj比為0.82，接近于帆5 位置。數(shù)量采用占傘衣幅數(shù)的百分比定義，根據(jù)圖1可知國外已應(yīng)用的開“窗”數(shù)量為20%?？紤]到本文基準(zhǔn)算例中的環(huán)帆傘傘衣幅數(shù)為24，且開“窗”數(shù)量需為整數(shù)，最終在不改變傘衣構(gòu)型設(shè)計(jì)的前提下，以25%和帆5 位置為基礎(chǔ)選取5 種不同的結(jié)構(gòu)分析開“窗”位置和數(shù)量對環(huán)帆傘氣動(dòng)特性的影響，各研究對象結(jié)構(gòu)參數(shù)如表2 所示。

表2 研究算例Table 2 Research examples

2 控制方程

采用任意拉格朗日-歐拉（ALE）方法對環(huán)帆傘開傘過程進(jìn)行流固耦合仿真［19-24］。

1） 流場域控制方程

式中：ρf為流體密度；t為時(shí)間；f為體積力；E為能量；vf和vm分別是流體和網(wǎng)格節(jié)點(diǎn)的速度；σf為流 體 的 應(yīng) 力 張 量，σf=-pδ+μ(grad(vf)+grad(vf)T)；p為壓強(qiáng)；δ為Kroneckerδ函數(shù)；μ為動(dòng)力黏度。

2） 結(jié)構(gòu)域控制方程

3） 流固耦合實(shí)現(xiàn)

對上述控制方程采用罰函數(shù)法進(jìn)行結(jié)構(gòu)與流場的全耦合計(jì)算，主要策略是追蹤結(jié)構(gòu)節(jié)點(diǎn)與流體節(jié)點(diǎn)之間的相對位移d。若以dn表示時(shí)間t=tn時(shí)的相對位移，并采用中心差分法按時(shí)間遞增進(jìn)行求解，則按以下等式對其進(jìn)行更新：

式中：vr表示結(jié)構(gòu)節(jié)點(diǎn)與流體節(jié)點(diǎn)之間相對速度，；Δt為時(shí)間間隔。

當(dāng)發(fā)生穿透，即在ns·dn＜0 時(shí)進(jìn)行耦合，罰函數(shù)耦合力Fp分布到結(jié)構(gòu)節(jié)點(diǎn)與流體節(jié)點(diǎn)上，ns表示與結(jié)構(gòu)節(jié)點(diǎn)相連的結(jié)構(gòu)單元的平均法線方向；若未發(fā)生穿透，就不進(jìn)行任何操作。

罰函數(shù)耦合力的大小與穿透深度呈正比，即

式中：i為笛卡爾坐標(biāo)系中的3 個(gè)坐標(biāo)軸方向；ki為i方向上基于主從節(jié)點(diǎn)模型特性的剛度系數(shù)；di為i方向上的穿透深度；界面接觸耦合力Fi的大小通常與罰函數(shù)耦合力的大小相等，但由于本文考慮織物透氣性，需基于Ergun 公式對罰函數(shù)耦合力進(jìn)行修正，同時(shí)傘衣采用薄殼單元，故i=xn，即僅考慮法向的耦合：

式中：an和bn為法向的黏性系數(shù)和慣性系數(shù)；vn為法向的流體速度；ε為孔隙率；S為單元表面積。結(jié)構(gòu)節(jié)點(diǎn)的界面接觸耦合力Fs=-Fxn；而對于流體，界面接觸耦合力需基于流體節(jié)點(diǎn)j的形函數(shù)［25］Nj分配給各個(gè)節(jié)點(diǎn)，即Fjf=Nj Fxn。

3 仿真模型及驗(yàn)證

3.1 半折疊狀態(tài)下環(huán)帆傘模型

開傘過程流固耦合仿真中，通常需要有一定的進(jìn)氣口面積才能保證傘衣的順利充氣。本文選擇具有一定初始展開的傘衣外廓作為初始狀態(tài)（初始展開直徑與名義直徑比為0.047），即建立半折疊狀態(tài)下的傘衣仿真模型，繼而進(jìn)行后續(xù)模擬計(jì)算。本研究中半折疊狀態(tài)下的傘衣模型基于數(shù)值計(jì)算方法獲得，主要利用了頂孔繩結(jié)點(diǎn)和傘繩結(jié)點(diǎn)的同步位移。傘衣采用三角形單元，減小折疊過程中出現(xiàn)網(wǎng)格穿透的概率。數(shù)值計(jì)算方法相對于直接建模方法計(jì)算快，人工干預(yù)因素小。最終，經(jīng)數(shù)值計(jì)算后生成的傘衣半折疊模型如圖5 所示，網(wǎng)格數(shù)為26 040。

圖5 傘衣半折疊模型Fig.5 Half-folding model of canopy

3.2 流場域模型

通常來講，環(huán)帆傘開傘過程數(shù)值仿真中，流場邊界至降落傘距離在5D0以上為宜，但考慮到大變形流固耦合仿真所需的計(jì)算資源龐大，本文流場域取的相對較小，沿降落傘軸向、法向、側(cè)向的流場計(jì)算域尺寸為5D0×3D0×3D0。同時(shí)為更好地捕捉流場特征，在綜合考慮計(jì)算精度與時(shí)間成本的基礎(chǔ)上，對流場的網(wǎng)格劃分采用十字形加密方式，以在傘衣附近區(qū)域建立更為細(xì)密的網(wǎng)格，確保降落傘傘衣附近的流體網(wǎng)格尺寸接近降落傘傘衣單元尺寸。另外在傘衣加密部分沿流場兩邊軸向方向的網(wǎng)格采用過渡方式，即網(wǎng)格逐漸稀疏，以避免出現(xiàn)質(zhì)量較差的網(wǎng)格。最終流場域的六面體網(wǎng)格數(shù)為276 萬，建立的流場計(jì)算域及仿真模型整體網(wǎng)格結(jié)構(gòu)如圖6 所示。

圖6 仿真模型網(wǎng)格結(jié)構(gòu)Fig.6 Mesh structure of simulation model

仿真中，由于該環(huán)帆傘開傘時(shí)的來流馬赫數(shù)為0.216，因此整個(gè)流場仿真的速度區(qū)間為低速段（可以不考慮空氣壓縮性影響）。同時(shí)，由于傘衣面積較小，開傘過程通常在零點(diǎn)幾秒內(nèi)即已完成，空氣組分、密度、壓強(qiáng)等參數(shù)的時(shí)間變化相對平緩，因此可采用恒定大氣密度和壓強(qiáng)進(jìn)行仿真。流場頂端和壁面均添加無反射約束，用以模擬遠(yuǎn)場狀態(tài)，使仿真結(jié)果與環(huán)帆傘置于無限大流場的情形更為接近。

3.3 模型驗(yàn)證

計(jì)算中開傘速度為70 m/s，開傘高度為4 000 m，初始攻角為0°，仿真時(shí)間為3 s，研究對象為基準(zhǔn)算例。以傘衣外形、阻力系數(shù)及開傘動(dòng)載系數(shù)作為考察對象，比較仿真結(jié)果與理論及實(shí)際試驗(yàn)之間的異同。

數(shù)值模擬得到的傘衣外形如圖7 所示。由圖7 可以發(fā)現(xiàn)：來流經(jīng)傘衣底邊進(jìn)氣口進(jìn)入后，形成1 個(gè)由傘衣包裹的管狀區(qū)域（0.1 s），之后傘衣頂部的空氣開始聚集，張開程度逐漸超過了傘衣底部（0.2 s），但隨后受到結(jié)構(gòu)慣性以及張力阻礙減緩了擴(kuò)張（0.4 s），同時(shí)張力使得進(jìn)氣口進(jìn)一步擴(kuò)大，引起傘衣的迅速張滿（0.6 s）。傘衣外形的變化與降落傘充氣理論描述基本吻合。圖8 為仿真計(jì)算與空投試驗(yàn)獲得的環(huán)帆傘張滿外形對比，由圖8 可知，二者基本一致。

圖7 開傘過程中的傘衣外形變化Fig.7 Canopy shape changes during inflation

圖8 仿真計(jì)算結(jié)果和空投試驗(yàn)張滿外形對比Fig.8 Comparison of canopy shape between numerical simulation and airdrop test

圖9為傘衣阻力系數(shù)仿真時(shí)間歷程曲線。以1.5 s 后阻力系數(shù)的平均值作為平均阻力系數(shù)。由圖9 可知：開傘過程中傘衣最大阻力面積為49.371 m2，開傘結(jié)束后該環(huán)帆傘的平均阻力面積為36.075 m2，對應(yīng)的平均阻力系數(shù)為0.717 2，與空投試驗(yàn)值0.82 較為接近，偏差為12.5%?？紤]到影響傘衣阻力系數(shù)的因素眾多，包括傘衣折疊模型、網(wǎng)格劃分質(zhì)量、空投試驗(yàn)誤差、仿真算法適用性以及環(huán)境參數(shù)真實(shí)程度等，認(rèn)為以上偏差可以接受。另外，根據(jù)圖9 可得，開傘動(dòng)載系數(shù)kd為1.368，與環(huán)帆傘動(dòng)載系數(shù)≥1.1 相符［26］。綜上認(rèn)為，本文建立的仿真模型可以有效模擬環(huán)帆傘的開傘過程。

圖9 傘衣阻力系數(shù)仿真時(shí)間歷程曲線Fig.9 Time history curve of drag coefficient

4 開“窗”結(jié)構(gòu)影響

4.1 傘衣外形

圖10 為各算例開傘過程中的投影面積與名義面積之比隨時(shí)間變化曲線，同時(shí)仿真計(jì)算中t= 0.1，0.2，0.4，0.6 s時(shí)傘衣外形如圖11 所示。就仿真結(jié)果整體而言：開“窗”后的環(huán)帆傘傘衣投影面積與名義面積比與基準(zhǔn)算例相差最大不超過15%；采用不同開“窗”位置及數(shù)量的環(huán)帆傘開傘狀態(tài)良好，未出現(xiàn)非對稱充氣、無法充滿等現(xiàn)象。這是由于上述傘衣模型均采用對稱開“窗”方式，同時(shí)傘衣主要承力部件為加強(qiáng)帶，故未對傘衣外形在開傘過程中的變化造成較大影響。

圖11 不同開“窗”方式下開傘過程中的傘衣外形變化Fig.11 Canopy shape changes during inflation process of research examples

4.2 氣動(dòng)特性參數(shù)

圖12為各算例的氣動(dòng)特性參數(shù)。根據(jù)圖12，可以發(fā)現(xiàn)開“窗”位置及數(shù)量均對該環(huán)帆傘的氣動(dòng)特性產(chǎn)生一定影響。其中，圖12（d）阻力系數(shù)標(biāo)準(zhǔn)差表示各算例的平均阻力系數(shù)在1.5～3.0 s 內(nèi)的波動(dòng)程度。

圖12 各算例的氣動(dòng)特性參數(shù)Fig.12 Aerodynamic characteristic parameters of research examples

相比于基準(zhǔn)算例，數(shù)量不變，不同開“窗”位置變化如下：帆4 和帆5 位置的平均阻力系數(shù)分別減小1.56%、2.60%，帆3 位置平均阻力系數(shù)增加1.40%；帆 3、帆4、帆5 位置的開傘動(dòng)載系數(shù)依次減小2.5%、6.8%、4.5%；帆4 位置的阻力系數(shù)標(biāo)準(zhǔn)差增大27.9%，而帆3 和帆5 位置阻力系數(shù)標(biāo)準(zhǔn)差分別減小14.5%、21.1%。

相比于基準(zhǔn)算例，位置不變，不同開“窗”數(shù)量變化如下：25%、50%、75%數(shù)量的平均阻力系數(shù)依次減小2.6%、11.6%、12.8%；開傘動(dòng)載系數(shù)依次減小4.5%、4.8%、9.1%；阻力系數(shù)標(biāo)準(zhǔn)差依次減小21.1%、41.4%、36.2%。

綜上，選擇開“窗”位置時(shí)，相比于帆3 和帆4，在帆5 位置進(jìn)行開“窗”可在平均阻力系數(shù)略有下降，開傘動(dòng)載系數(shù)降低的同時(shí)，最大程度減小阻力系數(shù)的波動(dòng)；選擇開“窗”數(shù)量時(shí)，相比于50%和75%，以25%的數(shù)量進(jìn)行開“窗”可在開傘動(dòng)載系數(shù)降低、阻力系數(shù)波動(dòng)減小的同時(shí)，最大程度減小對平均阻力系數(shù)的影響。

4.3 擺動(dòng)角度

圖13 為擺動(dòng)角度仿真時(shí)間歷程曲線。根據(jù)圖13 可以看出：總體而言，開“窗”結(jié)構(gòu)會(huì)顯著減小 最大擺角，減小幅度在5°～10°，其中帆5 開“窗”50%和帆5 開“窗”75%時(shí)的減小幅度最大；開“窗”數(shù)量不變時(shí)，帆3 和帆5 位置的擺動(dòng)角度最小；開“窗”位置不變時(shí)，50%和75%數(shù)量的擺動(dòng)角度最小且基本相同。

受計(jì)算資源以及求解能力限制，流固耦合計(jì)算對流場中旋渦、壓力的捕捉精度有限。為分析上述曲線最大擺動(dòng)角度的變化原因，基于計(jì)算流體力學(xué)方法獲得如圖14、圖15 所示擺角20°時(shí)的流場速度分布云圖及壓力云圖。

圖14 擺角20°時(shí)流場速度分布云圖Fig.14 Nephogram of flow field velocity distribution at 20°

圖15 擺角20°時(shí)，壓力分布云圖Fig.15 Nephogram of pressure distribution at 20°

根據(jù)圖14、圖15 可以看出：各算例的傘衣內(nèi)外部流場相似，傘衣外部左右兩側(cè)均存在一個(gè)范圍較大的高速區(qū)域（低壓區(qū)）。開“窗”以后，傘衣內(nèi)部的部分氣流從“窗口”流出，一方面，在傘衣左側(cè)與外部氣流匯合，阻擾外部高速區(qū)域生成并使其遠(yuǎn)離傘衣，致使從環(huán)縫、月牙縫流出的射流偏折角度減小，導(dǎo)致“窗口”上側(cè)傘衣幅外部壓力降低；另一方面，在傘衣右側(cè)，破壞了沿頂孔至傘衣底邊內(nèi)部流場的連續(xù)性，造成傘衣右側(cè)內(nèi)部整體壓力降低；二者共同作用下，降落傘回正作用增強(qiáng)，最大擺角減?。环? 較帆4 最大擺動(dòng)角度較小的原因可能是越靠近傘衣底部，傘衣左側(cè)“窗口”上側(cè)傘衣幅外部壓力降低的區(qū)域面積越大；而帆3 較帆5 最大擺動(dòng)角度較小的原因可能是越靠近傘衣頂部，傘衣右側(cè)沿頂孔至傘衣底邊內(nèi)部流場的連續(xù)性越差，致使右側(cè)內(nèi)部整體壓力下降越多。

開“窗”數(shù)量越多，傘衣內(nèi)外流場受上述作用影響越大，這可能是造成最大擺角大幅度減小的原因。同時(shí)50%和75%開“窗”數(shù)量時(shí)的傘衣內(nèi)外部流場基本相同，表明50%以上開“窗”數(shù)量已對傘衣隨機(jī)擺動(dòng)過程中的流場變化已無明顯影響，最大擺角因開“窗”數(shù)量增加而減小的程度已達(dá)極限。

5 結(jié) 論

開“窗”結(jié)構(gòu)有助于提高群傘系統(tǒng)的同步性和穩(wěn)定性。本文基于流固耦合方法針對單傘開展不同開“窗”結(jié)構(gòu)的環(huán)帆傘開傘過程數(shù)值計(jì)算，探究了開“窗”位置和數(shù)量對環(huán)帆傘傘衣外形、氣動(dòng)特性參數(shù)、擺動(dòng)角度的影響，所獲結(jié)論如下：

1） 由于為對稱開“窗”，且傘衣主要承力部件為加強(qiáng)帶，故不同的開“窗”位置和數(shù)量對開傘過程中傘衣外形變化影響較小，未出現(xiàn)非對稱充氣、無法充滿等現(xiàn)象。

2） 開“窗”數(shù)量相同時(shí)，在帆5 位置開“窗”可在平均阻力系數(shù)基本不變的情況下，降低開傘動(dòng)載系數(shù)，最大程度地減小阻力系數(shù)的波動(dòng)。

3） 開“窗”位置相同時(shí)，以25%數(shù)量開“窗”可在降低開傘動(dòng)載系數(shù)、減小阻力系數(shù)波動(dòng)的同時(shí)，最大程度地減小對平均阻力系數(shù)的影響。

4） 開“窗”以后，傘衣內(nèi)部氣流從“窗口”流出，一方面在傘衣左側(cè)與外部氣流匯合后，阻擾外部高速區(qū)域生成并使其遠(yuǎn)離傘衣，致使從環(huán)縫、月牙縫流出的射流偏折角度減小，導(dǎo)致“窗口”上側(cè)傘衣幅外部壓力降低；另一方面在傘衣右側(cè)，破壞了沿頂孔至傘衣底邊內(nèi)部流場的連續(xù)性，造成傘衣右側(cè)內(nèi)部整體壓力降低。二者共同作用下，環(huán)帆傘的最大擺角發(fā)生顯著減小，且減小幅度隨著開“窗”數(shù)量的增加而增大，但與開“窗”位置之間并無明顯規(guī)律。