吸氣式高超聲速飛行器助推分離過程數(shù)值仿真*

聞 訊,柳 軍,夏智勛

(國防科技大學(xué) 空天科學(xué)學(xué)院, 湖南 長沙 410073)

吸氣式高超聲速飛行器是典型的高動態(tài)臨近空間飛行器，具有飛行速度快、機(jī)動性能好等優(yōu)點(diǎn)[1-3]，是近年來國際軍事研究熱點(diǎn)，其相關(guān)研究已列入多國戰(zhàn)略技術(shù)部署計(jì)劃[4-5]。為了達(dá)到啟動條件，吸氣式高超聲速飛行器必須借助運(yùn)載器提速，并在指定彈道高度實(shí)現(xiàn)助推分離，助推分離過程屬于典型的多體分離問題[6-9]。

針對多體分離問題，目前已開展多項(xiàng)相關(guān)研究并取得了一定進(jìn)展。張海瑞等[10]采用重疊網(wǎng)格技術(shù)對子母彈分離過程進(jìn)行了數(shù)值研究，研究發(fā)現(xiàn)，滿足俯仰角速度分離指標(biāo)比滿足俯仰角條件更苛刻。Chamberlain等[11]對整流罩分離問題進(jìn)行了高精度數(shù)值仿真研究，研究發(fā)現(xiàn)，在拋罩分離指令下達(dá)后的一段時間內(nèi)，罩體的尾流區(qū)域仍影響后體的流場分布。多體分離過程通常存在多體干擾問題[12-15]，王元靖等[16]在超聲速多體分離實(shí)驗(yàn)中發(fā)現(xiàn)載荷模型的氣動特性受分離位置的影響而變化顯著。Togashi 等[17]采用非結(jié)構(gòu)網(wǎng)格對投彈問題進(jìn)行了數(shù)值仿真研究，得到了分離過程中清晰的激波干擾圖像，且與風(fēng)洞紋影實(shí)驗(yàn)結(jié)果吻合。Li等[18]通過改變分離條件對火箭轉(zhuǎn)級過程進(jìn)行了研究，得到不同分離距離下的級間區(qū)流場分布差異。然而，關(guān)于多體分離問題的快速非定常數(shù)值計(jì)算方面的研究較少，尤其針對前體飛行器通流快速非定常數(shù)值計(jì)算的研究，對其分離規(guī)律和機(jī)理仍沒有明確的結(jié)論。

隨著網(wǎng)格構(gòu)建技術(shù)的逐漸發(fā)展，多種網(wǎng)格模型可應(yīng)用于多體分離問題的求解，主要包括整體網(wǎng)格變形技術(shù)、動網(wǎng)格技術(shù)、局部網(wǎng)格重構(gòu)技術(shù)以及嵌套網(wǎng)格技術(shù)等[19]。這使得復(fù)雜構(gòu)型飛行器的非定常數(shù)值研究成為可能。

本文針對帶圓頂泄流以及減速板布局的試飛器構(gòu)型，采用非結(jié)構(gòu)網(wǎng)格下的局部網(wǎng)格重構(gòu)技術(shù)，提出適用于求解復(fù)雜構(gòu)型多體分離問題的六自由度非定常方法。并采用重疊網(wǎng)格技術(shù)加快收斂，對該構(gòu)型在四組弱干擾冷態(tài)分離條件下進(jìn)行數(shù)值仿真，得到試飛器在不同分離條件下的流場布局和氣動特性。通過對仿真結(jié)果的分析得到最優(yōu)來流分離條件。

1 物理模型

1.1 幾何模型

帶級間泄流結(jié)構(gòu)以及減速板構(gòu)造的轉(zhuǎn)級階段試飛器氣動模型如圖1所示，該試飛器主要由吸氣式高超聲速飛行器、級間段和助推器三部分組成，氣動模型總長約8 m。飛行器采用頭部進(jìn)氣、多模塊進(jìn)氣道以及中心燃燒的固體沖壓式發(fā)動機(jī)構(gòu)型，尾部共有四片舵面和導(dǎo)流塊，呈“×”字布局。級間段采用圓頂泄流結(jié)構(gòu)，助推器主體為直徑0.6 m圓柱結(jié)構(gòu)，沿助推器左右兩側(cè)對稱布置四片減速板(見圖2)，全部張開時相對于飛行器縱向?qū)ΨQ面呈迎風(fēng)70°，助推器尾部共有四片尾翼，呈“×”字布局，翼展約2 m。

圖1 試飛器氣動外形圖Fig.1 Vehicle geometric shape

圖2 分離組件局部放大圖Fig.2 Enlargement of separation components

1.2 基準(zhǔn)網(wǎng)格

計(jì)算所使用的基準(zhǔn)網(wǎng)格是借助ICEM軟件生成的，根據(jù)分離模型特點(diǎn)對氣動模型表面做了適當(dāng)加密。生成模型的表面網(wǎng)格后，采用八叉樹方法生成計(jì)算域的體網(wǎng)格，并進(jìn)行網(wǎng)格質(zhì)量優(yōu)化。最終生成四面體網(wǎng)格規(guī)模約為112萬，多重網(wǎng)格層數(shù)為5。圖3和圖4分別給出了氣動模型縱向?qū)ΨQ面網(wǎng)格和體軸水平截面網(wǎng)格示意圖。

圖3 網(wǎng)格外形縱向?qū)ΨQ面網(wǎng)格Fig.3 Grid profile on the longitudinal symmetrical surface

圖4 網(wǎng)格外形體軸水平截面網(wǎng)格Fig.4 Grid profile on the horizontal center surface

2 數(shù)值方法及分離條件

2.1 六自由度運(yùn)動方程及求解

建立與前體固連的慣性坐標(biāo)系和與后體(分離剛體)固連的本體坐標(biāo)系。在慣性坐標(biāo)系中求解有運(yùn)動邊界條件的流動問題，三維N-S 方程如式(1)所示。

(1)

式中：Q為守恒變量；Fc為無黏通量；Fv為黏性通量；xc為網(wǎng)格移動速度。

本體坐標(biāo)系下，針對分離剛體可以得到式(2)～(4)。

(2)

(3)

(4)

其中：M為外力對剛體質(zhì)心的合力矩；I為剛體慣性矩張量；ωx，ωy，ωz分別為剛體的角速度ω在本體坐標(biāo)系下三個方向的分量；φ為滾轉(zhuǎn)角；θ為俯仰角；φ為偏航角。

2.2 分離條件

非定常六自由度計(jì)算的初始條件是通過定常計(jì)算得到的，數(shù)值模擬的來流條件為：馬赫數(shù)4.3，高度14 km，助推器分離的初始角速度和初始線速度均為零，其他初始參數(shù)見表1。

表1 助推器分離初始參數(shù)

基于上述基礎(chǔ)工況再參考不同來流攻角和側(cè)滑角條件設(shè)計(jì)四組算例，用于考察來流條件波動對安全分離結(jié)果的判斷和影響，算例條件及編號說明在表2 中給出。

采用四核2.67 GHz主頻、12 GB內(nèi)存的 Inter Core i7 CPU 920處理器開展本算例的非定常計(jì)算工作，完成每個時間狀態(tài)點(diǎn)的循環(huán)迭代計(jì)算約需要6.5 h，若取用30 ms作為非定常計(jì)算的時間間隔，則完成300 ms的單狀態(tài)計(jì)算周期為3 d。

表2 算例來流參數(shù)說明

3 結(jié)果與討論

采用非定常的計(jì)算方法得到試飛器在四種分離條件下飛行器YZJ和助推器ZTQ的氣動阻力和分離速度，定義分離力Δ為ZTQ和YZJ兩者阻力的差值，并將其作為一項(xiàng)分離指標(biāo)，Δ為正值時即可實(shí)現(xiàn)分離。另外，定義分離加速度比Rat為YZJ加速度比與ZTQ加速度的商值，不大于1的分離加速度比被認(rèn)為對分離有利。

3.1 分離過程

圖5給出了各算例在分離過程中的機(jī)體阻力及氣動分離力，圖6給出了各算例在分離過程中的分離加速度及分離加速度比。其中，帶攻角狀態(tài)結(jié)果用實(shí)心符號表示、帶側(cè)滑角狀態(tài)結(jié)果用實(shí)線表示。由圖5可見，在分離過程中YZJ氣動阻力在分離初始時刻稍有波動，約在100 ms以后流場穩(wěn)定，阻力基本不再變化。從狀態(tài)對比的數(shù)值上看，A5B0和A5B2狀態(tài)阻力較A0B2和A0B0狀態(tài)阻力略大。分離力曲線變化趨勢主要受ZTQ阻力發(fā)展規(guī)律影響，總體上先稍有減小后迅速增大，且?guī)Чソ菭顟B(tài)的增長趨勢更強(qiáng)。不難發(fā)現(xiàn)，圖6中的加速度與圖5中的阻力之間存在正比例關(guān)系。另外從圖6中右側(cè)坐標(biāo)值所對應(yīng)的分離加速度比數(shù)值來看，分離加速度比在整個分離過程中均保持在安全值0.8以內(nèi)，無攻角狀態(tài)最大值不超過0.2，其余狀態(tài)最大值約0.3，滿足安全分離需求。

圖5 機(jī)體阻力及氣動分離力對比曲線Fig.5 Comparison curve of body forces and aerodynamic separation forces

圖6 分離加速度及分離加速度比對比曲線Fig.6 Comparative curve of separating acceleration and separating acceleration ratio

助推器線速度分布如圖7所示，助推器角速度分布如圖8所示。由圖7線速度變化曲線可知，X方向分離運(yùn)動最為明顯，隨時間發(fā)展近似呈二次曲線變化規(guī)律，助推器分離運(yùn)動由快至慢的工況順序依次是A5B2、A5B0、A0B2、A0B0。分析最佳分離來流條件，除了分離速度可作為一個參考指標(biāo)以外，分離的安全性、穩(wěn)定性、可靠性都很重要，本文從角速度以及周向線速度的角度來考核后三項(xiàng)重要指標(biāo)。從Y方向的線速度結(jié)果分析來看，無攻角對應(yīng)兩組工況的助推器分離線速度更小，即認(rèn)為工況A0B0和A0B2的Y方向分離更滿足上述指標(biāo)；而在Z方向無側(cè)滑角對應(yīng)的兩組工況在Z方向的線速度基本為零，假設(shè)控制系統(tǒng)完善，鑒于安全分離的考慮，工況A0B2是更合適的選擇。由圖8角速度分離曲線可知，A0B0和A5B0狀態(tài)Y方向角速度基本為零，Z方向角速度相對很小，可認(rèn)為其對應(yīng)翻轉(zhuǎn)運(yùn)動較弱。A5B2和A5B0在X方向運(yùn)動最為明顯，Y方向運(yùn)動和Z方向運(yùn)動對應(yīng)最明顯的狀態(tài)是A0B2和A5B2，且角速度增長較快，伴隨有非常明顯的滾轉(zhuǎn)運(yùn)動；同時，分離過程中A5B2狀態(tài)助推器的滾轉(zhuǎn)運(yùn)動較A5B0和A0B2更為明顯，A0B2其次，A0B0狀態(tài)飛行器滾轉(zhuǎn)效應(yīng)最弱。

圖7 助推器線速度分布Fig.7 Line velocity distribution for booster

圖8 助推器角速度分布Fig.8 Angular velocity distribution for booster

計(jì)算得到300 ms內(nèi)試飛器在A0B0狀態(tài)下縱向?qū)ΨQ面流場密度云圖如圖9～14所示，圖中Rho表示無量綱密度，是流場當(dāng)?shù)孛芏扰c來流密度的比值，該變量的流場等值線分布在一定程度上可反映激波的形態(tài)。在初始時刻，來流在前體飛行器頭部產(chǎn)生弓形激波，進(jìn)氣道捕捉大部分流量，激波在內(nèi)流道中經(jīng)內(nèi)壁面多次反射交叉后在噴管處膨脹，然后經(jīng)級間段泄流的氣體在后體助推器頭部圓頂位置再次產(chǎn)生明顯的弓形激波。在60 ms時刻，試飛器整體的激波形態(tài)接近初始時刻，從級間段附近網(wǎng)格變化可知后體助推器位置略后移。在90 ms時刻，分離距離繼續(xù)增大。分離過程對前體飛行器流場形態(tài)產(chǎn)生微弱影響，前體尾噴管出口低壓區(qū)面積增大。在此基礎(chǔ)上，分離距離在150 ms和210 ms時刻進(jìn)一步增大，且助推器姿態(tài)及位置發(fā)生較明顯的變化。至300 ms時刻已基本飛離前體尾流影響域。

圖9 0 ms時A0B0狀態(tài)流場密度云圖Fig.9 Flow field density contour for A0B0 on 0 ms

圖10 60 ms時A0B0狀態(tài)流場密度云圖Fig.10 Flow field density contour for A0B0 on 60 ms

3.2 氣動分析

圖11 90 ms時A0B0狀態(tài)流場密度云圖Fig.11 Flow field density contour for A0B0 on 90 ms

圖12 150 ms時A0B0狀態(tài)流場密度云圖Fig.12 Flow field density contour for A0B0 on 150 ms

圖13 210 ms時A0B0狀態(tài)流場密度云圖Fig.13 Flow field density contour for A0B0 on 210 ms

圖14 300 ms時A0B0狀態(tài)流場密度云圖Fig.14 Flow field density contour for A0B0 on 300 ms

圖15～18分別給出了飛行器在300 ms分離過程中的升力系數(shù)、阻力系數(shù)、側(cè)向力系數(shù)以及俯仰力矩系數(shù)隨時間的變化曲線。從結(jié)果分析得到，工況A5B2和A5B0氣動結(jié)果接近，但側(cè)向力系數(shù)有差異，這主要是側(cè)滑角方向相對更強(qiáng)的來流條件所致。同時，工況A0B2和A0B0狀態(tài)下的飛行器氣動結(jié)果接近，且相對受分離擾動影響更小、數(shù)據(jù)波動更小。分離初始時期阻力系數(shù)變化是分離過程中噴管處壓力場差異造成的。從俯仰力矩數(shù)據(jù)來看：工況A5B0和A5B2狀態(tài)下飛行器存在低頭力矩且有逐漸增大的趨勢；而工況A0B0和A0B2存在一個小的抬頭力矩且波動較小，這對于飛行器維持靜穩(wěn)定飛行狀態(tài)更有利。

圖15 升力系數(shù)隨分離時間變化曲線Fig.15 Lift coefficient curve on the separation time

圖16 阻力系數(shù)隨分離時間變化曲線Fig.16 Drag coefficient curve on the separation time

圖17 側(cè)向力系數(shù)隨時間變化曲線Fig.17 Lateral coefficient curve on the separation time

圖18 俯仰力矩系數(shù)隨時間變化曲線Fig.18 Pitch coefficient curve on the separation time

當(dāng)分離距離大于3倍前體直徑，即大于1.2 m時即可認(rèn)為助推器安全分離。圖19為四組算例的飛行器與助推器的分離矢量距離隨時間的變化曲線。由圖可知，0.3 s內(nèi)四組工況的分離距離隨時間的變化曲線十分接近，均呈拋物線型分布，這說明該方案在弱干擾冷態(tài)分離條件下的適應(yīng)性較好，同時0.3 s計(jì)算狀態(tài)下四組工況均實(shí)現(xiàn)了助推器的有效分離。另得到0.1 s時分離距離約為0.3 m，0.15 s時的分離距離在0.6～0.7 m范圍內(nèi)，0.3 s時的分離距離約為2.8 m。同時，由插值方法得到的分離距離對應(yīng)所需要的時間及對應(yīng)時刻的分離力如表3 所示。

圖19 300 ms分離距離隨時間變化曲線Fig.19 Separating curves on the 300 ms separation time

分離距離/mA5B0狀態(tài)A0B2狀態(tài)A5B2狀態(tài)A0B0狀態(tài)時間/ms分離力/N時間/ms分離力/N時間/ms分離力/N時間/ms分離力/N0018 104025 829017 793025 8790.5128.821 901128.422 872127.722 202128.723 4931.0183.427 796182.828 434181.924 938182.927 4221.5223.133 465222.931 460221.725 645223.030 4232.0255.236 148255.834 454255.139 150256.232 349

4 結(jié)論

本文采用動網(wǎng)格技術(shù)和非定常數(shù)值計(jì)算方法評估了帶減速板和泄流圓頂分離結(jié)構(gòu)的助推分離方案的可靠性，解決了助推安全分離問題，并得到各時間節(jié)點(diǎn)下的分離結(jié)構(gòu)的流場分布。通過對流場結(jié)構(gòu)的分析得到：在分離初期，助推器級間段的圓頂形成了流體滯止區(qū)，使前體飛行器尾噴管流場壓強(qiáng)增大，導(dǎo)致飛行器氣動阻力減小，進(jìn)而有利于助推分離；隨著分離距離的增加，泄流結(jié)構(gòu)對前體阻力干擾減弱，飛行器阻力增大對安全分離產(chǎn)生不利因素，而減速板加快實(shí)現(xiàn)了助推器的分離，能夠增加中后期的安全分離的可靠性。

通過對氣動力結(jié)果的分析發(fā)現(xiàn)，在初始分離中，前體飛行器的氣動阻力存在先減小再逐漸增大的變化趨勢，隨后保持在一個較初始值更高的穩(wěn)定值，即在整個分離過程中存在一個分離力極小值，這個結(jié)論與流場顯示結(jié)果吻合。通過對氣動參數(shù)的變化規(guī)律以及分離數(shù)據(jù)的分析發(fā)現(xiàn)，來流干擾因素越少對分離越有利，即攻角和側(cè)滑角均為0時為最佳分離條件。