龔小權，賈洪印，趙 輝，唐 靜，張 健，*，付云峰

(1. 中國空氣動力研究與發(fā)展中心 計算空氣動力研究所，綿陽 621000；2. 空氣動力學國家重點實驗室，綿陽 621000)

0 引 言

重復使用[1]是天地往返飛行器大幅降低發(fā)射成本、減小發(fā)射周期、提高航天應用的重要措施。一方面，可重復使用航天飛機及空間軌道機動飛行器X-37B 的成功試飛使得火箭基兩級入軌（two stage to orbit, TSTO）重復使用成為較為現(xiàn)實的飛行模式。另一方面，隨著吸氣式?jīng)_壓發(fā)動機的技術進步，使用渦輪/吸氣沖壓組合動力的一級飛行器水平起降TSTO重復使用飛行方案也逐步取得突破性進展。

馮毅等[2]、唐偉等[3]提出基于TBCC 動力的TSTO 氣動布局概念設計，綜合考慮兩級分離安全性，將雙垂尾外移并適當偏轉為V 尾，同時適當增大下折翼梢以更大程度地利用超聲速、高超聲速激波升力，并增加航向穩(wěn)定性；通過不斷迭代優(yōu)化，最后獲得一級飛行器垂尾下反方案，同時調(diào)整二級飛行器的質心位置至64%附近，以進一步拓寬分離邊界。

TSTO 并聯(lián)飛行器級間分離是其設計時必須考慮的關鍵問題之一。兩級飛行器外形尺寸巨大且重量為數(shù)百噸量級。以唐偉設計的TSTO 為例，一級飛行器長度約85 m，二級飛行器長35 m，重量都在百噸以上。一方面，分離中如果兩級接觸，接觸力極有可能破壞飛行器結構或防熱層。另一方面，分離初始時刻兩級級間距離相當近，只有300 mm，相對兩級的尺寸來說太小，略微俯仰過快就會帶來頭部或者體襟翼的碰撞，這對分離初始階段兩級飛行器的姿態(tài)控制提出了較大要求。另外TSTO 兩級飛行器都是升力體布局，在一定的攻角范圍內(nèi)兩級都有較大的升力，且升力與重力同一量級，這將導致兩級不易分開，進一步增加了兩級碰撞的風險。TSTO 并聯(lián)分離兩級完全分開的時間很長，從能量管理角度考慮，設計師希望分離過程平穩(wěn)，盡量減少高度及能量損失。因此研究TSTO 并聯(lián)分離過程對TSTO 外形設計、分離點選擇及分離方案設計至關重要。

Liu 等[4]研究了TSTO 在不同速度下的分離，分析了組合體在超聲速以及不同襟翼偏角下的氣動特性，研究了不同馬赫數(shù)下無側滑角的縱向分離。采用國家數(shù)值風洞工程通用CFD 軟件NNW-FlowStar[5]，本文研究了TSTO 飛行器在有側滑角時的并聯(lián)級間分離過程，給出了分離過程一二級的姿態(tài)角及質心位移，并與試驗值進行了對比；分析了分離過程的安全性，驗證了FlowStar 軟件數(shù)值模擬TSTO 類并聯(lián)分離的能力和精度。文章主要由四個部分組成：首先是數(shù)值模擬方法簡介，介紹數(shù)值模擬軟件NNW-FlowStar、數(shù)值模擬用到的技術以及六自由度求解方法；然后在第二部分重點介紹TSTO 并聯(lián)分離的計算數(shù)模、網(wǎng)格，并采用WPFS 標準算例驗證了軟件的非定常計算精度；第三部分是分離安全性分析，首先分析了TSTO 組合體氣動特性，給出了TSTO 可能安全分離的分離點，然后數(shù)值模擬了在馬赫數(shù)Ma= 6、來流攻角α= ?2°、側滑角β= 2°、高度H= 30 Km 時，TSTO的并聯(lián)分離過程，分析了分離中兩級飛行器的姿態(tài)角變化及質心位移。最后對全文進行了總結。

1 數(shù)值方法

1.1 軟件介紹

數(shù)值模擬軟件為中國空氣動力研究與發(fā)展中心基于國家數(shù)值風洞工程研發(fā)的通用CFD 軟件NNWFlowStar。該軟件內(nèi)核MFlow 是基于非結構混合網(wǎng)格、用于求解亞跨超聲速流場的專業(yè)CFD 求解器[6-9]。NNW-FlowStar 能處理任意形狀的網(wǎng)格單元，具有較強的靈活性。采用有限體積法離散控制方程，求解變量存儲于網(wǎng)格單元的體心。在本文數(shù)值模擬中，對流通量采用Roe 通量差分分裂格式進行離散，該格式具有很高的間斷和黏性分辨率[10-11]。單元內(nèi)使用線性重構使得空間離散具有二階精度，變量梯度求解使用節(jié)點型Gauss 方法[12]，限制器采用Venkatakrishnan[13]限制器，該限制器具有較高的精度和良好的殘差收斂性。粘性通量采用中心差分離散。假定流場為全湍流，采用Spalart-Allmaras 一方程湍流模型[14]，湍流模型方程空間離散采用一階迎風格式。

守恒形式的非定常Navier-Stokes（N-S）方程為：其中，為守恒變量，H為無黏通量，Hv為黏性通量，n表示控制體單元邊界面的外法矢。

采用雙時間步方法求解整個非定常過程，由于非定常計算工作周期較長，本文內(nèi)迭代采用LU-SGS（lower upper-symmetric Gauss-Seidel）隱式時間離散[15]，并在迭代中引入局部時間步長加速收斂技術加速收斂。采用基于MPI （message passing interface） 的大規(guī)模并行計算[16]縮短計算周期。重疊挖洞及插值相關技術見文獻[17]。

1.2 六自由度求解

FlowStar 軟件在初始狀態(tài)的機體體軸系中求解分離物體的質心運動方程，在分離物體的體軸系（非慣性系）中求解分離物體的繞質心運動方程。采用四階R-K 方法求解分離物體的運動方程（2）～方程（4），得到質心的速度以及物體的角速度，從而得到物體的姿態(tài)角。

其中，φ為滾轉角，θ為俯仰角，ψ為偏航角，MA表示氣動力力矩，MT表示除氣動力外的其他力的力矩（外力力矩）。

2 計算模型、網(wǎng)格及軟件驗證

2.1 計算模型

本文研究模型是大型TSTO 并聯(lián)飛行器，由唐偉設計提供。圖1 給出了該模型的三視圖，下面級長度為80 m，重心距離下面級飛行器實際頭部尖點（55.45，0，?1.07）約65%自身長度位置；上面級長度為35 m，重心距離下面級飛行器實際頭部尖點（55.45，0，4.00）約64%自身長度位置。兩級間最短距離0.3 m。詳細尺寸及質量特性見文獻[3-4]及表1。

表1 TSTO 尺寸及質量特性Table 1 Geometry size and mass characteristics of TSTO

圖1 TSTO 外形三視圖Fig. 1 Three views of the TSTO

2.2 計算網(wǎng)格

本文的計算網(wǎng)格為非結構混合網(wǎng)格，包括六面體、三棱柱、金字塔、四面體網(wǎng)格。采用六面體和三棱柱模擬附面層，四面體模擬空間流場各向同性區(qū)域，通過金字塔過渡三棱柱和四面體。為節(jié)省網(wǎng)格量并較好地模擬前緣，在兩級飛行器機翼、垂尾表面前后緣盡可能采用四邊形網(wǎng)格，并將四邊形網(wǎng)格的長寬比控制在100 以內(nèi)。對于一些規(guī)則且獨立的部件，例如體襟翼、副翼、升降舵表面，全部采用結構四邊形網(wǎng)格布置，如圖2、圖3。圖2 給出了TSTO 一級的表面網(wǎng)格，圖3 給出了二級的表面網(wǎng)格。從圖中看到，在兩級表面外形變化較大的區(qū)域都有加密網(wǎng)格，以保證計算網(wǎng)格與數(shù)模的一致性。針對二級翼身組合體外形，在機翼機身連接處加密網(wǎng)格。

圖2 TSTO 一級表面網(wǎng)格Fig. 2 Surface grid of the first stage of TSTO

圖3 TSTO 二級表面網(wǎng)格Fig. 3 Surface grid of the second stage of TSTO

圖4 給出了TSTO 一二級分離初始時刻空間網(wǎng)格的y= 0 截面。從圖中看到，為保證重疊網(wǎng)格插值精度，將一級飛行器空間區(qū)域中二級飛行器可能的運動區(qū)域的空間網(wǎng)格加密（圖中外層較大的加密區(qū)）。

圖4 一二級分離初始時刻對稱面網(wǎng)格Fig. 4 Grid distribution in the symmetry plane at the initial time instance of separation

TSTO 一二級距離非常近，僅0.3 m，且分離的馬赫數(shù)Ma= 6，在一級機身上表面和二級機身下表面存在嚴重的激波邊界層干擾以及激波反射。因此針對TSTO 一級機身上表面和二級頭部及機身下表面進行網(wǎng)格局部加密，如圖5、圖6，以期更好地捕捉兩級間的激波反射，提高對氣動特性的模擬精度。

圖5 二級頭部附近對稱面網(wǎng)格Fig. 5 Grid distribution around the nose of the second stage in the symmetry plane

圖6 兩級機身縫隙間的對稱面網(wǎng)格Fig. 6 Grid distribution around the two-stage fuselage gaps in the symmetry plane

2.3 WPFS 模型分離數(shù)值模擬

機翼/掛架/帶舵外掛物模型WPFS（wing/pylon/finned-store）[18]是美國發(fā)起第一次分離投放CFD 驗證時使用的一個機翼和導彈的簡化模型。該模型具有比較翔實的風洞數(shù)據(jù)，因此被各種CFD 軟件用做驗證分離計算的標準算例。圖7 給出了WPFS 模型視圖。為確保外掛物和機翼/掛架的安全分離，試驗過程中使用了虛擬的彈力作用，在分離距離大于給出的安全距離后，虛擬彈力的作用消失。試驗條件、外掛物質量特性和虛擬彈力的相關參數(shù)見文獻[18]。

圖7 WPFS 模型Fig. 7 WPFS model

包裹機翼和掛架的背景網(wǎng)格有網(wǎng)格單元239 萬，其中四面體183 萬、三棱柱56 萬；包裹外掛物的子網(wǎng)格有網(wǎng)格單元208 萬，其中四面體101 萬、三棱柱107 萬。采用非定常雷諾平均N-S 方程、SA 湍流方程來模擬計算該模型的外掛物分離歷程。計算條件和CTS （captive trajectory simulation）試驗條件相同，計算中真實時間步長取0.01 s。通過計算得到了外掛物質心位置、外掛物姿態(tài)角隨時間變化曲線。

圖8 給出了計算得到的帶舵外掛物質心位置與試驗結果的對比曲線。計算結果和試驗值吻合很好。

圖8 計算和試驗的外掛物質心隨時間變化對比Fig. 8 Displacement variation with time for the centroid of the external store (computation VS experiment)

圖9 給出了計算得到的帶舵外掛物姿態(tài)角與試驗結果的對比曲線。計算得到的外掛物姿態(tài)角中的偏航角和滾轉角與試驗結果吻合得較好；但俯仰角卻比試驗值偏大，在0.3 s 時最大偏差值有0.8°（約19%）。分析其原因可能是數(shù)值模擬直接選擇了在0.054 s 時取消彈射力，和試驗中嚴格地根據(jù)距離取消彈射力有一定差別，這種累計效應隨著時間的增長顯現(xiàn)了出來。

圖9 計算和試驗的外掛物姿態(tài)角隨時間變化對比Fig. 9 Attitude angle variation with time for the external store(computation VS experiment)

通過計算結果和試驗結果的比較分析可以看出，計算結果與試驗數(shù)據(jù)吻合得較好，表明FlowStar 軟件的數(shù)值方法和非結構重疊網(wǎng)格方法的可靠性較高，具有較強的實用性，能夠較好地模擬多體干擾與分離復雜流場及物體分離投放軌跡。

3 分離安全性分析

3.1 組合體氣動特性分析

TSTO 兩級飛行器安全分離對于整個方案設計至關重要。一方面，TSTO 兩級飛行器都是升力體布局，且分離初始時刻相距很近，這對分離方式提出了較高的要求；另一方面，從能量管理角度，設計師期望正攻角下分離，從而避免TSTO 高度及能量損失。因此研究TSTO 設計分離工況點（Ma= 6，H= 30 Km）的分離特性異常重要。

本節(jié)首先分析了TSTO 組合飛行器在設計分離工況的定常氣動特性，用于選擇分離點。圖10 給出了組合體一二級升力系數(shù)隨來流攻角變化的曲線，從圖中可以看到，一級飛行器升力系數(shù)隨著攻角增大而增大，攻角大于1°之后升力系數(shù)為正。二級飛行器在計算的攻角內(nèi)升力都為正值，二級飛行器完全處于一級飛行器的背風區(qū)，其升力系數(shù)隨著攻角增大而減小。圖11 給出了組合體一二級飛行器俯仰力矩系數(shù)隨來流攻角變化的曲線，在計算的攻角范圍內(nèi)一級飛行器的俯仰力矩都是負值，二級飛行器的俯仰力矩都是正值，這代表分離初始階段兩級會相對旋轉，而不是同向旋轉，這對分離極為不利。從圖10、圖11 看到，在計算的攻角范圍內(nèi)都是不利于分離的。二級具有較大的抬頭力矩，由于兩級間距離很小，初始距離為0.3 m，過大的抬頭力矩容易導致一二級相撞。

圖10 一二級升力系數(shù)隨攻角變化曲線Fig. 10 Lift coefficient variation with the angle of attack for the two stages

圖11 一二級俯仰力矩系數(shù)隨攻角變化曲線Fig. 11 Pitching moment coefficient variation with the angle of attack for the two stages

表2 給出了一二級舵面不偏轉時的定常氣動特性。從表中看到，在來流攻角為?2°、0°、2°時，一級的總的法向力（氣動力加重力）都為負，二級的總的法向力（氣動力加重力）為負。因此在這三個攻角下分離一、二級都要損失高度和能量。從表2 中還可以看到，在攻角為?2°時，二級法向力最大，一級法向力最小，單從這方面考慮，此狀態(tài)是最有利于分離的；但是此時二級的俯仰力矩較大，從這個方面看，該狀態(tài)又不利于分離。

表2 無舵偏定常氣動力Table 2 Steady aerodynamic forces without rudder deflection

3.2 分離過程模擬

從3.1 節(jié)的分析看到，在計算攻角內(nèi)分離都不是特別理想，但是這只是初始狀態(tài)。隨著分離進行，兩級間距離逐漸增大，兩級間復雜的激波干擾以及兩級相對位置、姿態(tài)的變化也導致兩級的氣動特性急劇變化。

本節(jié)選取來流攻角?2°，側滑角2°來分析兩級的級間分離過程。計算采用兩套網(wǎng)格，如圖4 所示，包括包裹一級飛行器的網(wǎng)格和包裹二級飛行器的網(wǎng)格。包裹一級飛行器的網(wǎng)格量約為2122 萬，其中，四面體網(wǎng)格單元1254 萬個、三棱柱單元472 萬個、六面體單元342 萬個、金字塔單元54 萬個。包裹二級飛行器的網(wǎng)格量約為1288 萬，其中四面體網(wǎng)格單元623 萬個、三棱柱單元448 萬個、六面體單元179 萬個、金字塔單元38 萬個。CTS 試驗在中國空氣動力研究與發(fā)展中心Φ1 m量級高超聲速風洞FD-30 中進行[19]。該風洞在2019 年建成雙軌跡捕獲試驗系統(tǒng)，具備開展Ma= 6 并聯(lián)分離雙分離軌跡捕獲試驗的能力。為減少支撐干擾影響，試驗采用雙尾支撐結構，如圖12。圖13 給出了CTS 試驗中的典型流場紋影圖，從圖中能夠看到，級間分離過程中兩級間存在復雜的激波干擾以及激波邊界層干擾。試驗在求解六自由度方程時，將滾轉力矩Mx置零，也就是將方程（3）中的滾轉力矩置零。為與試驗一致，計算過程也采用相同的處理。計算來流條件為Ma= 6、α= ?2°、β= 2°、H= 30 Km，采用全尺寸外形，在中國空氣動力研究與發(fā)展中心計算集群上采用1024 核完成整個非定常計算。

圖13 CTS 典型流場紋影圖[19]Fig. 13 Schlieren image of the typical flowfield in CTS[19]

圖14 、圖15 給出了TSTO 兩級在分離過程中兩級姿態(tài)角隨分離時間變化的曲線。在計算的分離時間內(nèi)兩級能夠安全分離，計算和試驗的俯仰角和偏航角隨時間變化具有很好的一致性，進一步驗證了NNW-FlowStar 的數(shù)值模擬精度。從圖14 一級飛行器的姿態(tài)角變化曲線看，計算和試驗獲得的偏航角和滾轉角變化規(guī)律完全一致。在2°側滑角下，一級飛行器具有正偏航角，頭部逐漸轉向來流側滑方向，減小側滑角；由于滾轉力矩置零，滾轉角很小。一級飛行器在分離初始階段具有較小的低頭力矩，隨著兩級的移動，尤其是一二級法向方向的位移，導致二級的頭激波及其反射激波在一級上表面的位置后移，一級飛行器俯仰力矩逐漸增大，并逐步抬頭。

圖14 Ma = 6、α = ?2°、β = 2°，一級飛行器姿態(tài)角隨時間變化曲線Fig. 14 Attitude angle variation with time for the first stages at Ma = 6，α = ?2°，β = 2°

圖15 Ma = 6、α = ?2°、β = 2°，二級飛行器姿態(tài)角隨時間變化曲線Fig. 15 Attitude angle variation with time for the second stages at Ma = 6，α = ?2°，β = 2°

從圖15 的二級姿態(tài)角變化曲線看，計算和試驗的二級俯仰角和偏航角變化規(guī)律一致。在上一節(jié)中分析組合體的氣動特性時看到，二級飛行器具有較大的抬頭力矩，因此在分離過程中二級飛行器一直保持抬頭，且抬頭較快；隨著兩級的分離，二級飛行器的頭激波的反射激波逐步后移且強度逐步減弱，影響二級飛行器機身和機翼下表面壓力分布，導致俯仰力矩斜率變化。二級飛行器具有負偏航角，頭部逐漸轉向來流側滑的反方向，增大側滑角。

圖16、圖17 給出了一二級飛行器質心位移隨時間變化曲線。在2°側滑角下，兩級飛行器均向飛行員左側移動，計算與試驗的兩級側向移動位移變化規(guī)律及大小一致。從圖16、圖17 看到，一二級飛行器均向下運動，計算與試驗的兩級質心法向位移變化規(guī)律及量值一致。由于一級飛行器為負升力，二級飛行器為正升力，因此相同時間時，一級向下運動距離更大。試驗和計算的軸向位移有一定差別，隨著分離進行，兩者的差量增大。在1.6 s 時刻一級飛行器計算和試驗的軸向位移相差約1.5 m，這與阻力計算精度相關；而阻力計算精度又與計算網(wǎng)格尺寸、湍流模型、兩級后體分離等多種因素相關。但總的來說，相比飛行器尺寸（85 m），這一差量較小。

圖16 Ma = 6，α = ?2°，β = 2°，一級質心位移隨時間變化曲線Fig. 16 Displacement variation with time for the centroids of the first stages at Ma = 6，α = ?2°，β = 2°

圖17 Ma = 6，α = ?2°，β = 2°，二級質心位移隨時間變化曲線Fig. 17 Displacement variation with time for the centroids of the second stages at Ma = 6，α = ?2°，β = 2°

圖18 給出了Ma= 6、α= ?2°、β= 2°，TSTO 兩級并聯(lián)分離不同分離時刻兩級對稱面的壓力云圖以及表面俯視圖。分離過程兩級間流動復雜，二級飛行器的頭部激波打到一級飛行器上表面并在二級飛行器機身下表面和一級飛行器上表面之間反射。在分離的初始階段，一二級飛行器的最小距離出現(xiàn)在二級體襟翼附近，二級飛行器體襟翼下表面出現(xiàn)高壓。隨著分離時間增大，二級頭激波在一級上表面的位置逐步后移，一級上表面的高壓區(qū)逐步后移。在計算的時間范圍內(nèi)兩級仍未相互脫離干擾區(qū)。

圖18 Ma = 6、α = ?2°、β = 2°，不同時刻兩級對稱面及表面壓力云圖Fig. 18 Pressure distributions in the symmetry plane and on the surface of the two stages at different time instances for Ma = 6，α = ?2°，β = 2°

4 結 論

本文基于國家數(shù)值風洞工程通用CFD 軟件NNW-FlowStar，采用各向異性的非結構混合網(wǎng)格以及動態(tài)重疊網(wǎng)格數(shù)值模擬了TSTO 一二級飛行器并聯(lián)分離。通過與CTS 試驗值對比一二級飛行器質心位移以及姿態(tài)角變化曲線，驗證了軟件對升力體并聯(lián)級間分離類問題的數(shù)值模擬能力。通過數(shù)值模擬，可以得到以下結論：

1）NNW-FlowStar 軟件在重疊插值算法、高分辨率Roe 格式、高斯節(jié)點型梯度計算方法、大規(guī)模并行等一系列數(shù)值算法方面的開發(fā)和改進保證了高超聲速流場的精細模擬。

2）NNW-FlowStar 軟件具備開展TSTO 升力體并聯(lián)級間分離類數(shù)值模擬能力，模擬結果與CTS 試驗值具有較高的一致性，驗證了軟件在高超飛行器并聯(lián)分離數(shù)值模擬的精度。

3）TSTO 并聯(lián)級間分離關系到兩級飛行器設計，目前的分離工況能夠安全分離，但是安全余量非常小，分離中兩級最小距離出現(xiàn)在二級體襟翼處，如果考慮到兩級飛行器尺寸、氣動彈性以及其他因素，這種最小距離是不安全的。

4）需要進一步研究TSTO 并聯(lián)級間分離方案，包括來流攻角、預置舵偏等，以期能夠增加分離過程兩級的安全距離，在減少TSTO 能量損失的同時，實現(xiàn)兩級快速分離。

致謝：感謝中國空氣動力研究與發(fā)展中心林敬周研究員提供的CTS 試驗數(shù)據(jù)。