胡俊,于菲,于勇
(北京理工大學(xué) 宇航學(xué)院,北京 100081)
降落傘作為一種可展開(kāi)式氣動(dòng)減速器,具有質(zhì)量小、折疊體積小、價(jià)格便宜、阻力特性好等優(yōu)點(diǎn),被廣泛應(yīng)用于空投、航天器回收、新型武器等領(lǐng)域.
降落傘工作過(guò)程包括拉直、充氣、減速和穩(wěn)降,其中的拉直和充氣階段歷時(shí)短、變形大,涉及幾何非線性和材料非線性并存的流固耦合問(wèn)題[1].降落傘的實(shí)驗(yàn)研究包括風(fēng)洞實(shí)驗(yàn)、飛行實(shí)驗(yàn)等方法.如 LEVIN 等[2]通過(guò)風(fēng)洞試驗(yàn)測(cè)量了5 種不同傘面形狀的十字型傘的氣動(dòng)力以及俯仰力矩,并與空投實(shí)驗(yàn)數(shù)據(jù)進(jìn)行了結(jié)果對(duì)比,得到不同形狀、不同臂長(zhǎng)比十字型傘的氣動(dòng)力規(guī)律.RENDER 等[3]、LUDTKE[4]、BRAUN等[5]主要對(duì)小型十字傘傘群進(jìn)行了一系列風(fēng)洞試驗(yàn),研究了物傘連接繩長(zhǎng)度和傘面之間系留繩長(zhǎng)度對(duì)傘群系統(tǒng)氣動(dòng)力的影響,對(duì)2 種不同尺寸的傘群系統(tǒng)進(jìn)行不同的雷諾數(shù)范圍內(nèi)的阻力測(cè)量,得出了傘群阻力系數(shù)與單傘阻力系數(shù)之間的關(guān)系式.NASA 蘭利研究中心[6]多次采用飛行實(shí)驗(yàn),對(duì)“行星進(jìn)入降落傘計(jì)劃”(planetary entry parachute program,PEPP)中不同構(gòu)型的降落傘,進(jìn)行有關(guān)穩(wěn)定性、阻力性能和降落傘尾跡的研究.
隨著計(jì)算機(jī)技術(shù)的發(fā)展,數(shù)值模擬也成為降落傘研究的一種手段.通過(guò)將穩(wěn)定的降落傘視為剛體,模擬研究降落傘的氣動(dòng)特性和周圍流場(chǎng).如XUE等[7-8]采用CFD 方法針對(duì)探測(cè)器超聲速降落傘系統(tǒng)進(jìn)行了詳細(xì)研究,觀察了探測(cè)器尾流與傘前激波、探測(cè)器激波與傘前激波的氣動(dòng)干擾以及傘繩激波、風(fēng)洞連桿激波對(duì)整體流場(chǎng)的影響情況.周彤等[9]、李堯堯等[10]也采用了CFD 的方法對(duì)傘-彈系統(tǒng)中的相關(guān)氣動(dòng)問(wèn)題進(jìn)行了研究.BARNHARDT 等[11]則采用分離渦模擬方法研究了NASA“火星科學(xué)實(shí)驗(yàn)室任務(wù)”相關(guān)的超聲速降落傘系統(tǒng),模擬了超聲速來(lái)流下的探測(cè)器-降落傘系統(tǒng),研究了前體探測(cè)器尾流對(duì)降落傘氣動(dòng)性能的影響,分析了尾流對(duì)傘衣壓力、阻力的影響以及傘衣帶底部有可能會(huì)出現(xiàn)部分坍塌的原因.
降落傘在穩(wěn)定減速和穩(wěn)定下降階段時(shí),將傘衣視為剛體是可行的,但是在降落傘的拉直和充氣階段,傘衣會(huì)出現(xiàn)結(jié)構(gòu)大變形,不能再將傘衣視為剛體,此時(shí)數(shù)值模擬主要采用任意拉格朗日-歐拉(arbitrary Lagrangian-Eulerian,ALE)法.2011 年TUTT 等[12]首次采用ALE 方法模擬降落傘有限質(zhì)量充氣過(guò)程,結(jié)果驗(yàn)證了ALE 方法對(duì)開(kāi)傘力、下降速度及降落傘呼吸頻率有較好的預(yù)測(cè),驗(yàn)證了ALE 流固耦合方法模擬有限質(zhì)量開(kāi)傘過(guò)程的可行性.GAO 等[13]采用ALE算法模擬了開(kāi)縫傘的有限質(zhì)量充氣過(guò)程,用顯式中心差分法求解了不可壓縮流動(dòng)的Navier-Stokes (N-S)方程,得到了傘衣變形過(guò)程、減速規(guī)律及開(kāi)傘過(guò)載.
降落傘作為氣動(dòng)減速器,其主要作用是對(duì)載荷進(jìn)行減速,避免載荷在到達(dá)落點(diǎn)時(shí)產(chǎn)生過(guò)大的沖擊,從而達(dá)到保護(hù)載荷的作用.但傳統(tǒng)的降落傘在運(yùn)動(dòng)過(guò)程中無(wú)法對(duì)物傘系統(tǒng)的運(yùn)動(dòng)狀態(tài)進(jìn)行主動(dòng)調(diào)整.為了提高物傘系統(tǒng)落點(diǎn)的精準(zhǔn)度,近年來(lái)開(kāi)始出現(xiàn)主動(dòng)控制降落傘的應(yīng)用研究.通過(guò)傘繩的收放,改變傘面構(gòu)型,從而影響改變傘面所受氣動(dòng)力,達(dá)到控制物傘系統(tǒng)運(yùn)動(dòng)狀態(tài)的目的.
降落傘的主動(dòng)控制最早由美國(guó)提出,運(yùn)用于精確控制空投系統(tǒng)的軌跡控制.BERGERON 等[14-15]在十字型傘的其中一片傘臂上添加了擾流片,通過(guò)擾流片的關(guān)閉和打開(kāi)改變傘衣內(nèi)部的流場(chǎng)結(jié)構(gòu),使傘衣變形,達(dá)到改變降落傘氣動(dòng)性能的目的.FIELDS 等[16]通過(guò)在傘面安裝驅(qū)動(dòng)裝置改變傘繩長(zhǎng)度,控制降落傘的旋轉(zhuǎn)或滑翔,根據(jù)空投實(shí)驗(yàn)數(shù)據(jù)開(kāi)發(fā)四自由度及PID 控制模型,將著陸位置圓誤差降低2~4 倍;李龍恩[17]通過(guò)動(dòng)力學(xué)建模的方法,研究收放局部傘繩實(shí)現(xiàn)對(duì)圓形降落傘控制的可行性,發(fā)現(xiàn)收放局部傘繩能夠?qū)崿F(xiàn)對(duì)降落傘的主動(dòng)控制.
有關(guān)學(xué)者對(duì)于降落傘主動(dòng)控制的研究主要集中在實(shí)現(xiàn)傘繩收放的控制方法與控制裝置等方面,但對(duì)于降落傘收縮傘繩產(chǎn)生的氣動(dòng)力規(guī)律沒(méi)有系統(tǒng)地認(rèn)識(shí).掌握降落傘收縮傘繩之后的運(yùn)動(dòng)規(guī)律和氣動(dòng)力變化規(guī)律是實(shí)現(xiàn)降落傘主動(dòng)控制應(yīng)用的前提條件,同時(shí)也影響著發(fā)送調(diào)控指令的時(shí)間節(jié)點(diǎn).不同傘繩收縮長(zhǎng)度時(shí)降落傘的氣動(dòng)特性,尤其是側(cè)向控制力的規(guī)律是降落傘主動(dòng)控制成功應(yīng)用的重要因素.本文以十字型降落傘為研究對(duì)象,建立收縮不同傘繩長(zhǎng)度的降落傘無(wú)限質(zhì)量充氣的有限元模型,采用ALE流固耦合方法,對(duì)收縮不同傘繩長(zhǎng)度時(shí)十字型傘的氣動(dòng)特性進(jìn)行研究,為用十字型降落傘進(jìn)行主動(dòng)控制提供參考.
柔性降落傘在運(yùn)動(dòng)過(guò)程中和空氣的相互作用具有強(qiáng)非線性時(shí)變的特點(diǎn),兩者的耦合是同時(shí)涉及幾何非線性和材料非線性的瞬間大變形雙向耦合問(wèn)題.ALE 算法是一種基于罰函數(shù)的接觸耦合算法,結(jié)合了拉格朗日算法(主要用于求解固體力學(xué))和歐拉算法(主要用于求解流體力學(xué))的優(yōu)點(diǎn),可以實(shí)現(xiàn)結(jié)構(gòu)和流體的雙向耦合.
ALE 描述下的流體模型遵循Navier-Stokes 方程,其中包括連續(xù)性方程、動(dòng)量守恒方程、能量守恒方程.
ALE 描述下的連續(xù)性(質(zhì)量守恒)方程為
ALE 描述下的動(dòng)量守恒方程為
ALE 描述下的能量守恒方程為
式中: ρ為流體密度;t為時(shí)間;vi為流體速度;wi為參考構(gòu)型下網(wǎng)格點(diǎn)移動(dòng)速度;bi為單位質(zhì)量體力;e為單位質(zhì)量?jī)?nèi)能; σij為應(yīng)力張量,表達(dá)式為:
其中:p為流體壓力; δij為克羅尼克函數(shù); μ為流體動(dòng)力黏度系數(shù).
ALE 方法對(duì)固體域采用Lagrange 描述.降落傘的結(jié)構(gòu)控制方程為
式中: ρs為降落傘織物密度;u為單元位移; σs為柯西應(yīng)力張量;f為單元所受體積力;g為重力加速度.
降落傘傘衣采用不考慮厚度的二維薄殼結(jié)構(gòu),材料選擇織物,傘衣織物的透氣量與透過(guò)傘衣的流體速度呈二次多項(xiàng)式關(guān)系[18]:
式中: Δp為織物兩側(cè)壓力差;e為織物厚度;a、b分別為黏性系數(shù)和慣性系數(shù),均是關(guān)于多孔介質(zhì)相對(duì)透氣量 ε和空氣動(dòng)力黏度 μ的多項(xiàng)式.
建立傘衣結(jié)構(gòu)網(wǎng)格和流體網(wǎng)格重疊的有限元模型,使用罰函數(shù)接觸耦合算法實(shí)現(xiàn)力學(xué)參量的傳遞,在任意tn時(shí)間步內(nèi),追蹤拉格朗日節(jié)點(diǎn)(從節(jié)點(diǎn))穿透歐拉流體物質(zhì)(主節(jié)點(diǎn))的相對(duì)位移dn,則下一時(shí)刻迭代結(jié)果為:
式中:vrel為主從節(jié)點(diǎn)相對(duì)速度;主從節(jié)點(diǎn)速度為等參坐標(biāo)變化后得到的某一流體質(zhì)點(diǎn)速度vf;從節(jié)點(diǎn)速度為結(jié)構(gòu)節(jié)點(diǎn)速度vs; Δt為時(shí)間步長(zhǎng).
十字型降落傘主要包括傘衣和傘繩兩部分,其中,傘衣部分是由5 片尺寸相同、材質(zhì)相同的織物縫制而成,十字型傘的幾何結(jié)構(gòu)如圖1 所示.L為傘長(zhǎng)臂,W為傘短臂,Ls為傘繩長(zhǎng)度.十字型傘的氣動(dòng)特性主要受其幾何參數(shù)的影響,其中主要幾何參數(shù)包括臂長(zhǎng)比A=L/W和傘繩比Rs=Ls/L.
如圖2 所示建立坐標(biāo)系,傘繩節(jié)點(diǎn)作為坐標(biāo)系原點(diǎn)O,來(lái)流方向?yàn)樽鴺?biāo)系的+Z軸,與Z軸垂直的平面為XY平面,建立風(fēng)軸坐標(biāo)系OXYZ;以壓力中心O′作為坐標(biāo)系原點(diǎn),傘軸所在直線OO′作為Z′軸,與Z′軸垂直的平面為X′Y′平面,建立體軸坐標(biāo)系O′X′Y′Z′.
風(fēng)軸坐標(biāo)系Z軸與體軸坐標(biāo)系Z′軸之間的夾角為 α.在風(fēng)洞實(shí)驗(yàn)中,通過(guò)固定降落傘傘頂進(jìn)行降落傘的靜態(tài)氣動(dòng)特性研究,該夾角就是降落傘靜態(tài)實(shí)驗(yàn)中的攻角.
圖3 所示為降落傘的受力示意圖.風(fēng)軸坐標(biāo)系中降落傘的氣動(dòng)力主要包括阻力FD和升力FL.體軸坐標(biāo)系中,降落傘的主要受力為軸向力FT和法向力FN.
計(jì)算氣動(dòng)力系數(shù)時(shí)所用到的參考面積為十字型傘的名義面積S0,表達(dá)式[19]為:
為驗(yàn)證ALE 方法在十字型降落傘研究中的適用性,以文獻(xiàn)[20]中風(fēng)洞試驗(yàn)的十字型降落傘為例,采用ALE 算法對(duì)對(duì)稱十字型降落傘進(jìn)行數(shù)值模擬,十字型降落傘的具體尺寸和計(jì)算條件如表1 所示.
表1 十字型降落傘的尺寸和來(lái)流條件[20]Tab.1 Parameters of cross parachute and flow condition [20]
采用上述降落傘尺寸,模擬了攻角為0°和15°時(shí)十字型傘由平面狀態(tài)充氣至穩(wěn)定運(yùn)動(dòng)的過(guò)程,取計(jì)算穩(wěn)定后的1~4 s 的軸向力數(shù)據(jù),計(jì)算時(shí)間平均值,結(jié)果如圖4 所示.
圖4 風(fēng)洞試驗(yàn)數(shù)據(jù)[20]和模擬數(shù)據(jù)對(duì)比Fig.4 Comparison of wind tunnel data [20] and simulation data
攻角為0°時(shí)采用數(shù)值方法計(jì)算出十字型傘的軸向力系數(shù)大小為0.822,風(fēng)洞試驗(yàn)所得軸向力系數(shù)為0.8,兩者誤差約為2.75%;攻角為15°時(shí)數(shù)值計(jì)算出軸向力系數(shù)為0.765,風(fēng)洞試驗(yàn)所得軸向力系數(shù)為0.712,兩者誤差約為7.4%.
綜上所述,數(shù)值模擬計(jì)算與風(fēng)洞試驗(yàn)數(shù)據(jù)誤差小于10%,認(rèn)為本文所使用的數(shù)值模擬方法適用于十字型降落傘模擬計(jì)算.
本文基于某火箭彈小角度迅捷轉(zhuǎn)彎制定了一種十字型降落傘,具體幾何尺寸如表2 所示.
表2 十字型降落傘的幾何尺寸Tab.2 Geometric parameters of cross parachute
每個(gè)傘臂連接2 根傘繩,整個(gè)降落傘系統(tǒng)共有8 根傘繩.本文目的是研究十字型降落傘在穩(wěn)定狀態(tài)下的氣動(dòng)特性,對(duì)開(kāi)傘過(guò)程不做精確要求,所以將降落傘的初始狀態(tài)設(shè)置為十字形平面狀態(tài),見(jiàn)圖1.
傘面結(jié)構(gòu)為薄殼單元,對(duì)傘面進(jìn)行正四邊形網(wǎng)格劃分,保證每個(gè)傘面網(wǎng)格的質(zhì)量均等.傘繩結(jié)構(gòu)采取離散梁?jiǎn)卧鲌?chǎng)域采取固體單元.十字型傘的流場(chǎng)域采用直徑為5L、高度為6L的圓柱形計(jì)算域,采用結(jié)構(gòu)性網(wǎng)格劃分,相關(guān)的有限元模型參數(shù)如表3所示.
表3 降落傘的有限元模型參數(shù)Tab.3 Finite element model parameters of parachute
流場(chǎng)域和傘面的網(wǎng)格劃分如圖5 所示.
圖5 網(wǎng)格示意圖Fig.5 Diagram of mesh
如圖6 所示,收縮傘臂A所連接的2 根傘繩,收縮長(zhǎng)度定義為Δl=Ls-,單位為m,收縮比定義為δ=Δl/Ls.
圖6 傘繩收縮示意圖Fig.6 Diagram of suspension line retracting
來(lái)流速度為68 m/s,來(lái)流密度為1.205 kg/m3,當(dāng)收縮長(zhǎng)度 Δl為0、0.16、0.20、0.24、0.32、0.40 m 時(shí),所對(duì)應(yīng)的收縮比δ及其他計(jì)算條件如表4 所示.
表4 數(shù)值模擬計(jì)算條件Tab.4 Numerical simulation conditions
通過(guò)傘繩的收縮,十字型降落傘由對(duì)稱性幾何構(gòu)型變成非對(duì)稱性幾何構(gòu)型,同時(shí)降落傘的氣動(dòng)力參數(shù)發(fā)生改變.
圖7 所示為0°、5°、10°攻角時(shí),十字型傘的瞬時(shí)阻力系數(shù)隨時(shí)間的變化.
圖7 不同收縮比十字型傘的阻力系數(shù)Fig.7 Drag coefficients of cross parachutes with different retraction ratios
選擇計(jì)算穩(wěn)定后的1~4 s 的數(shù)據(jù),取時(shí)間平均值,計(jì)算得到平均阻力系數(shù),如圖8 所示.
圖8 平均阻力系數(shù)隨收縮比 δ變化圖Fig.8 Retraction ratio vs mean drag coefficient
從圖中可以看出,十字型傘的平均阻力系數(shù)隨傘繩收縮比 δ的增大而減小.降落傘的阻力系數(shù)主要受阻力面積的影響,而十字型傘的傘面構(gòu)型相較于其他傘型,受同等數(shù)量的傘繩長(zhǎng)度縮短的影響程度較大,所以傘繩收縮時(shí),平均阻力系數(shù)的減小較為顯著.
從圖7(b)、7(c)可以看出,在5°、10°攻角時(shí),傘繩收縮的十字型傘的瞬時(shí)阻力系數(shù)隨時(shí)間有較為明顯的波動(dòng).這是因?yàn)楫?dāng)十字傘處于非零攻角且傘繩同時(shí)收縮時(shí),由于傘面的柔性,迎風(fēng)面積會(huì)隨著十字傘自旋發(fā)生變化,圖9 中展示了當(dāng)攻角為10°、收縮比δ=20.0%時(shí),十字型傘發(fā)生自旋運(yùn)動(dòng)時(shí),t1、t2時(shí)刻(如圖7(c)中收縮比20%、攻角10°工況,選取某個(gè)周期中阻力系數(shù)最大和最小的2 個(gè)典型時(shí)刻)下,十字型傘的傘面狀態(tài).由圖中可以明顯看出,t1時(shí)刻傘面的迎風(fēng)面積明顯大于t2時(shí)刻的迎風(fēng)面積.
圖9 t1 時(shí)刻和t2 時(shí)刻的傘面狀態(tài) (收縮比δ=20.0%)δ=20.0%Fig.9 State of canopy at t1 and t2 (retracting ratio )
圖10 所示為0°、5°、10°攻角時(shí),十字型傘的瞬時(shí)升力系數(shù)隨時(shí)間的變化.
圖10 不同收縮比 δ時(shí)十字型傘的升力系數(shù)Fig.10 Lift coefficient of cross parachute vs retraction ratio δ
選擇計(jì)算穩(wěn)定后的1~4 s,取時(shí)間平均值,計(jì)算平均升力系數(shù),如圖11 所示.傘的平均升力系數(shù)基本沒(méi)有影響.在5°和10°攻角時(shí),傘繩收縮會(huì)引起平均升力系數(shù)的降低.
圖11 平均升力系數(shù)隨收縮比 δ變化圖Fig.11 Retraction ratio δ vs mean lift coefficient
根據(jù)圖10(a)可以看出,0°攻角時(shí),瞬時(shí)升力隨著自旋運(yùn)動(dòng)呈現(xiàn)出明顯的周期性,升力的方向隨著下拉傘臂A旋轉(zhuǎn)到不同位置而發(fā)生變化;同時(shí)從不同收縮比的工況可以看出,十字型傘瞬時(shí)升力系數(shù)的峰值與傘繩收縮比 δ密切相關(guān).
由圖10(b)、10(c)所示,攻角為5°、10°時(shí),升力系數(shù)受到攻角和收縮比 δ這2 個(gè)因素的共同影響.當(dāng)傘繩收縮比 δ為0 時(shí),升力主要受攻角影響,方向始終指向+X方向,并且升力為正值;當(dāng) δ不為0 時(shí),升力除了受攻角影響之外,還包括來(lái)自傘繩收縮的影響.隨著下拉傘臂A的旋轉(zhuǎn),如圖12(b)、12(c)所示,傘臂A的位置從傘軸左側(cè)運(yùn)動(dòng)至傘軸右側(cè),導(dǎo)致升力時(shí)刻發(fā)生變化.隨著傘繩收縮比 δ的增加,瞬時(shí)升力系數(shù)的峰值增大,如t1時(shí)刻(在收縮比為20%、攻角
圖12 收縮比δ=0 與δ=20.0%的傘面壓力云圖 (10°攻角)Fig.12 Pressure contour of canopy with δ=0 and δ=20.0% (angle of attack was 10°)
從圖中可以看出,在0°攻角時(shí),十字型傘的平均升力系數(shù)基本為0,傘繩收縮比 δ對(duì)0°攻角的十字型為10°工況,選取某個(gè)周期升力最大的典型時(shí)刻,如圖10(c)所示).隨著傘繩收縮比 δ的增加,瞬時(shí)升力系數(shù)的最小值減小.當(dāng) δ較大時(shí),收縮傘繩對(duì)升力的影響大于攻角對(duì)升力的影響,瞬時(shí)升力系數(shù)最小值逐漸出現(xiàn)負(fù)值,如t2時(shí)刻(在收縮比為20%、攻角為10°工況,選取某個(gè)周期升力最小的典型時(shí)刻,如圖10(c)所示).
當(dāng)攻角為10°、收縮比δ=0 和δ=20.0%的t1、t2時(shí)刻,Y=0 平面的壓力云圖如圖12 所示,速度矢量圖如圖13 所示.
圖13 收縮比δ=0 與δ=20.0%的速度矢量圖 (10°攻角)Fig.13 Velocity vector of canopy with δ=0 and δ=20.0% (angle of attack was 10°)
由圖12 可以看出,在10°攻角、收縮比δ=0工況時(shí),傘衣內(nèi)部的壓強(qiáng)關(guān)于傘軸非對(duì)稱分布.在收縮比δ=20.0%工況時(shí),t1時(shí)刻,收縮傘繩在傘軸左側(cè),加劇了傘面幾何構(gòu)型的非對(duì)稱性,導(dǎo)致高壓區(qū)向傘軸右側(cè)產(chǎn)生了進(jìn)一步的偏移,此時(shí)傘面受到的升力最大;t2時(shí)刻,收縮傘繩在傘軸的右側(cè),下拉傘臂A在來(lái)流的沖擊下,出現(xiàn)傘面底邊向內(nèi)凹陷的現(xiàn)象,由于內(nèi)凹傘面的阻擋,來(lái)流在A處發(fā)生繞流,如圖13(c)所示,傘臂A內(nèi)側(cè)出現(xiàn)漩渦區(qū),出現(xiàn)較為明顯的低壓區(qū),而高壓區(qū)出現(xiàn)在傘臂C內(nèi)側(cè),如圖12(c)所示,此時(shí),傘面受到的升力最小,為負(fù)值.
本文對(duì)不同傘繩收縮比 δ的十字型降落傘進(jìn)行了數(shù)值模擬研究,結(jié)果表明:
① 十字型降落傘的阻力系數(shù)隨收縮比 δ的增大而減小.
② 由于十字傘存在自旋運(yùn)動(dòng),下拉傘臂A的位置隨著自旋發(fā)生改變,導(dǎo)致升力具有較強(qiáng)的時(shí)變特性.十字型傘瞬時(shí)升力峰值隨收縮比 δ增大而增大.
本文的研究?jī)?nèi)容是十字型傘作為主動(dòng)控制裝置的基礎(chǔ)性研究工作,后續(xù)將對(duì)收縮多根(超過(guò)2 根)傘繩以及不同傘繩收縮組合下引起的其他非對(duì)稱構(gòu)型對(duì)側(cè)向氣動(dòng)力的影響進(jìn)行研究,制定出合理的傘繩收縮策略,實(shí)現(xiàn)十字型傘對(duì)彈箭系統(tǒng)運(yùn)動(dòng)姿態(tài)的調(diào)控.