開“窗”結構對環帆傘開傘過程影響

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

1．北京空間機電研究所，北京 100094

2．南京航空航天大學 航空學院，南京 210016

3．中南大學 航空航天學院，長沙 410083

4．北京航空航天大學 計算機學院，北京 100191

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

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

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

為解決上述問題，國外經多次空投及風洞試驗設計修改，采用了開“窗”結構設計，發現該結構有利于提高群傘系統的同步性和穩定性［14-15］，圖1 為NASA 獵戶座飛船主傘采用的開“窗”結構。中國目前對此結構尚未應用也無相關公開仿真研究成果。

圖1 獵戶座飛船主傘的開“窗”結構Fig.1 ‘Windows’ structure of NASA Orion main parachute

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

1 研究對象

本文以某環帆傘為研究對象，開展開傘過程數值計算。其平面結構及具體尺寸如圖2 和表1所示。同時，為方便后續分析，圖2（a）對環1～4和帆1～6 進行了序號標注。

圖2 環帆傘結構Fig.2 Flat structure of ringsail canopy

表1 環帆傘結構尺寸Table 1 Structural size of ringsail parachute

由于加強帶等細化結構會影響流固耦合的仿真結果［16］，故本文的環帆傘模型由環、帆、傘繩、頂孔繩、徑向帶、緯向帶等結構構成，確保仿真模型與真實物理更為貼近。傘衣幅結構和整體繩帶結構如圖3、圖4 所示。

圖3 環帆傘傘衣幅結構Fig.3 Structure of ringsail canopy gore

圖4 整體繩帶結構Fig.4 Integral rope structure

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

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

環片：

帆片：

通常認為在最大投影直徑位置以下進行開“窗”對阻力特性影響較小，該傘空投試驗中的投影直徑Dt與結構直徑Dj比為0.82，接近于帆5 位置。數量采用占傘衣幅數的百分比定義，根據圖1可知國外已應用的開“窗”數量為20%。考慮到本文基準算例中的環帆傘傘衣幅數為24，且開“窗”數量需為整數，最終在不改變傘衣構型設計的前提下，以25%和帆5 位置為基礎選取5 種不同的結構分析開“窗”位置和數量對環帆傘氣動特性的影響，各研究對象結構參數如表2 所示。

表2 研究算例Table 2 Research examples

2 控制方程

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

1） 流場域控制方程

式中：ρf為流體密度；t為時間；f為體積力；E為能量；vf和vm分別是流體和網格節點的速度；σf為流 體 的 應 力 張 量，σf=-pδ+μ(grad(vf)+grad(vf)T)；p為壓強；δ為Kroneckerδ函數；μ為動力黏度。

2） 結構域控制方程

3） 流固耦合實現

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

式中：vr表示結構節點與流體節點之間相對速度，；Δt為時間間隔。

當發生穿透，即在ns·dn＜0 時進行耦合，罰函數耦合力Fp分布到結構節點與流體節點上，ns表示與結構節點相連的結構單元的平均法線方向；若未發生穿透，就不進行任何操作。

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

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

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

3 仿真模型及驗證

3.1 半折疊狀態下環帆傘模型

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

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

3.2 流場域模型

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

圖6 仿真模型網格結構Fig.6 Mesh structure of simulation model

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

3.3 模型驗證

計算中開傘速度為70 m/s，開傘高度為4 000 m，初始攻角為0°，仿真時間為3 s，研究對象為基準算例。以傘衣外形、阻力系數及開傘動載系數作為考察對象，比較仿真結果與理論及實際試驗之間的異同。

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

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

圖8 仿真計算結果和空投試驗張滿外形對比Fig.8 Comparison of canopy shape between numerical simulation and airdrop test

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

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

4 開“窗”結構影響

4.1 傘衣外形

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

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

4.2 氣動特性參數

圖12為各算例的氣動特性參數。根據圖12，可以發現開“窗”位置及數量均對該環帆傘的氣動特性產生一定影響。其中，圖12（d）阻力系數標準差表示各算例的平均阻力系數在1.5～3.0 s 內的波動程度。

圖12 各算例的氣動特性參數Fig.12 Aerodynamic characteristic parameters of research examples

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

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

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

4.3 擺動角度

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

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

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

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

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

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

5 結 論

開“窗”結構有助于提高群傘系統的同步性和穩定性。本文基于流固耦合方法針對單傘開展不同開“窗”結構的環帆傘開傘過程數值計算，探究了開“窗”位置和數量對環帆傘傘衣外形、氣動特性參數、擺動角度的影響，所獲結論如下：

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

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

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

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