高超聲速飛行器自適應減阻盤動態減阻機理

韓榮，劉偉，楊小亮

國防科技大學 空天科學學院，長沙 410073

為了改善頭部駐點處氣動加熱效應，高超聲速飛行器頭部通常采用鈍頭體外形。在高超聲速來流中，鈍頭體前方會形成較強烈的弓形激波，機體前后壓差陡增。研究表明，激波阻力占高超聲速飛行器整機阻力的80%以上［1］，導致飛行器有效載荷顯著降低，同時增加了控制難度［2］，給飛行安全帶來了極大的隱患。目前，高超聲速鈍體飛行器的減阻方法主要包括前緣凹腔、沿駐點線能量沉積、頭部加裝減阻桿以及它們的組合構型［3］等，其中減阻桿（頂端可加裝圓盤等不同外形，即減阻盤）是最為簡單有效的一種方法。

同時，未來飛行器對機動性和敏捷性提出了更高的要求。在機動飛行過程中，隨著攻角增加，傳統固定式減阻盤的效果迅速下降［4-5］。自適應減阻盤在一定程度上可以避免上述問題，即根據來流攻角自適應調整桿件方向，使減阻盤軸線始終和來流保持平行或較小偏角狀態。自適應減阻盤可以較好地解決有攻角狀態下減阻率下降的問題，對未來飛行器減阻系統可重復利用以及全攻角來流范圍適用具有重要價值。國內外相關研究人員對此進行了一系列風洞試驗及數值模擬研究。

Nicholson 等［6］在1968 年首次提出了自適應減阻桿的設計方案，借助肩部的皮托管監測飛行器迎風面與背風面的壓差，然后利用壓差驅動活塞、平行四邊形支架等機構，進而使減阻桿對準來流，防止飛行器頭部受到來流直接沖擊，以保護其內部結構以及所裝載的傳感器等設備。Schülein 等［7-12］采用被動控制方式，參考風向標工作原理設計了一種自適應減阻盤，其自適應過程主要由飛行器兩側翼面產生的氣動力矩驅動。研究發現，自適應減阻盤可以維持飛行器頭部回流區的對稱性，同時顯著降低再附激波強度。此外，耿云飛等［13］進一步設計了一種自適應無燒蝕減阻系統，通過在桿件前端噴射冷流，可以解決有攻角狀態下的減阻防熱以及桿件尖端燒蝕等問題。Deng 等［14］在2017 年將自適應的概念應用到了雙減阻盤結構上，對單減阻盤、雙減阻盤以及是否存在安裝角情況下的流場分別開展了數值仿真研究。結果顯示，當安裝角和攻角均為8°時，飛行器的升阻比可以達到3.634，相對原外形提高了9.1%。Huang 等［15］探討了安裝角對雙減阻盤防熱性能的影響，發現當桿件的安裝角比來流攻角大1.5°時，鈍頭體-減阻盤系統的防熱性能最佳。

此外，在高超聲速機動飛行過程中，飛行攻角、角速度以及角加速度等運動參數發生著劇烈變化，周圍流場呈現強烈的非定常效應。現有減阻方法研究主要集中在對高超聲速飛行器的靜態減阻性能及機理分析上，自適應減阻盤在動態過程中的現實效果以及流場特征必然和定常結果存在差異。高超聲速鈍頭體強迫俯仰振蕩問題有較多的試驗和工程計算結果［16-19］，而關于加裝自適應減阻盤后飛行器動態特性的研究則相對較少。

針對加裝自適應減阻盤的高超聲速飛行器開展動態研究是飛行器減阻設計的關鍵環節之一，具有重要的科學價值和工程意義。因此，本文通過非定常數值模擬方法，研究自適應減阻盤在高超聲速飛行器機動過程中的減阻效果和減阻機理，對比分析自適應減阻盤和傳統固定式減阻盤的差異，探討高超聲速飛行器氣動力特性以及流場特性隨減阻盤參數的變化規律，為自適應減阻盤的設計和選型提供參考。

1 物理模型和數值方法

1.1 物理模型

以Menezes 等［20］試驗中 的鈍頭 體外形 作為基準模型展開研究，幾何模型如圖1 所示。其中，鈍頭體頂角為120°，底部直徑D=100 mm。數值計算中，參考長度取鈍體底部直徑D，參考面積即底部面積πD2/4。

自由來流參數與試驗條件一致，其中來流馬赫數為5.75，來流靜溫為140 K，壁面設置為等溫壁，溫度為300 K，其余來流條件見表1。數值模擬基于此來流條件開展。

表1 試驗來流條件Table 1 Experimental conditions of incoming flow

1.2 非定常流動/運動耦合計算方法

1.2.1 非定常流體力學控制方程

飛行器機動過程涉及到網格邊界運動，為精確描述動邊界流動問題，基于任意拉格朗日歐拉（Arbitrary Lagrangian-Eulerian，ALE）方法建立坐標系。在ALE 坐標系中，積分形式的三維可壓縮非定常Navier-Stokes（N-S）方程表示為

式中：各量均為有量綱量；Q表示守恒變量，具體表達式為

采用來流密度ρ∞、來流速度U∞、來流動力黏度μ∞以及物體的特征長度L0作為參考量，將變量和方程進行無量綱化，轉換得到無量綱形式的控制方程

同時，采用AUSM 混合格式進行空間離散，時間離散使用雙時間步隱式LU-SGS 方法，詳情可參照文獻［21］。

1.2.2 運動學方程

采用強迫俯仰振蕩運動模型對高超聲速飛行器開展動態氣動特性研究，俯仰振蕩運動可以模擬飛行中的拉升、俯沖等常用機動方式，是典型的動態非定常問題。在振蕩過程中飛行器攻角按照如下規律變化：

式中：α0為初始攻角；αm表示振蕩幅值；t為無量綱時間；k是減縮頻率。具體定義如下：

本文主要研究高超聲速飛行器在大攻角機動中的動態減阻問題，因此初始攻角取0°，振蕩幅值αm=20°，振蕩周期為125 s。

1.2.3 流動/運動耦合算法

高超聲速飛行器機動問題需要涉及流體流動、飛行器運動及變形等多個學科領域，其中流體流動主要是求解N-S 方程，飛行器運動則通過求解運動學方程實現。本文通過圖2 所示的流動/運動耦合求解方案推進求解。

圖2 流動/運動耦合求解方案Fig. 2 Hydrodynamics/kinematics coupling method

在計算開始之后，首先需要對計算網格初始化，求解N-S 方程得到初始狀態下的流場信息。然后通過計算運動學方程，得到飛行器在下一時刻位置和角度等參數信息；進而執行網格運動或變形，更新計算網格。接下來便采用更新后的網格進行新一輪的流場模擬。計算完成后，對流場進行后處理，輸出飛行器氣動力、氣動力矩以及流場相關計算結果，完成計算。

為保證數值模擬精度和計算的穩定性，采用一階隱式歐拉格式的緊耦合方法［22］，通過隱式方式在每一時間步下構造多個子迭代，并在每一子迭代步驟中進行數據交換。

1.2.4 動態混合網格生成技術

高超聲速飛行器動態機動問題涉及物理邊界運動，其中需要處理的關鍵問題是動態網格的生成，本文使用的動態網格技術主要為基于徑向基函數（Radius Basis Function，RBF）的網格變形技術。

基于徑向基函數的網格變形方法主要分為求解和更新兩個模塊，其中求解模塊主要涉及插值系數的確定。

在RBF 插值過程中，“中心點”表示d維歐氏空間給 定的 一組位置不同點X=｛xc1，xc2，…，xcn｝? RD，各“中心點”上對應的一系列標量值為gc1，gc2，…，gcn。RBF 插值指的是當基函數φ（x）給定時，尋找合適的連續函數：

使得條件

得到滿足的過程。

顯然，插值函數f（x）需要通過所有的“中心點”，因此RBF 插值的關鍵就是要確定插值系數γi。具體到動網格問題，插值“中心點”表示物體運動邊界上的點，而相應的標量值指的是物體在3 個方向上的位移。插值系數通過式（9）確定：

式中：向量γ=[γ1γ2…γn]T；g=｛g（xci）｝ni=1；插值矩陣M為n階方陣，其元素為

其中：||x-xci||表示歐氏距離，在三維空間中，可以直接由r'替代：

上述插值系數的求解過程即RBF 方法的求解模塊，在求解得到位移插值系數之后，將計算域內各個網格點的坐標直接代入RBF 插值函數，便可求出內場每個網格點的位移，從而完成整個計算域內的網格移動，如式（12）所示。此過程即RBF 方法的更新模塊。

2 數值方法驗證及網格無關性分析

2.1 數值計算方法驗證

作為高超聲速導彈的標準模型，HBS（Hyper Ballistic Shape）動態特性既有試驗結果，也有半經驗理論預測數據［23］，可以較好地驗證程序的可靠性。HBS 外形的幾何參數如圖3 所示，其中：r=22.5 mm，θ1=5°，θ2=15°。計算來流條件與風洞試驗一致，其中來流馬赫數Ma∞=6.85，雷諾數Re=0.72×106（參考長度取頭部直徑），定姿態攻角及動態驗證初始攻角分別為0°，4°，8°，13.5°，16°和20°。

圖3 HBS 模型幾何外形示意圖Fig. 3 Geometry of HBS model

HBS 外形靜態氣動力系數隨攻角變化曲線如圖4 所示，圖中同時給出了文獻［24-25］的數值結果。該模型的升阻力系數均隨攻角逐漸增加，其中阻力系數先緩慢上升，而后在10°攻角附近增速加快。升力系數曲線斜率變化較小，基本上隨攻角線性變化。計算所得升阻力系數與參考結果沒有明顯差異。

圖4 HBS 外形靜態氣動力系數隨攻角變化曲線Fig. 4 HBS static aerodynamic coefficients varying with angle of attack

以定姿態計算所得流場作為非定常模擬的初場，分別模擬0°，4°，8°，13.5°，16°和20°初始攻角條件下的強迫俯仰振蕩運動，振幅均取1°，減縮頻率設置為0.05。俯仰阻尼導數的計算使用基于Ektin 非定常氣動力模型的強迫振蕩動導數辨識方法［26］。圖5（a）為不同初始攻角下的靜導數計算結果，數值模擬結果與參考結果的變化趨勢保持一致，隨著來流攻角增加，HBS 鈍錐體模型的靜穩定性參數從正變為負，基本在8°攻角之后轉變為靜穩定狀態。隨著攻角繼續增加，靜導數符號在17°左右再度改變，靜穩定性再次發生變化。數值結果與試驗值［23］基本吻合；與文獻［24］中的計算結果之間存在散布，靜導數在小攻角和大攻角處吻合較好。圖5（b）為俯仰阻尼動導數隨初始攻角的變化情況，模擬結果與參考結果基本一致，HBS 標模均處于俯仰動態穩定狀態。當來流攻角小于13.5°時，動導數變化較小，此后動導數迅速增加。在13.5°攻角附近結果與試驗值［23］存在一定差異，此處流動復雜需要進一步研究。

圖5 靜導數及動導數隨初始攻角變化曲線Fig. 5 Static and dynamic derivatives varying with initial angle of attack

綜上所述，程序模擬結果基本和試驗數據［23］以及文獻中的計算結果［24］相符合，說明所采用程序具備了非定常問題的求解能力，可以較好地模擬高超聲速飛行器動態運動問題。

2.2 網格無關性分析

針對圖1 所示物理模型，選取3 套不同密度網格進行網格無關性驗證，3 套網格所對應的壁面第1 層網格高度及網格雷諾數見表2。來流馬赫數為5.75，其余條件見1.1 節。

表2 飛行器壁面第1 層網格高度設置Table 2 First layer height of aircraft wall

壁面網格雷諾數的定義如下：

式中：ρ∞、v∞、μ∞分別表示來流密度、速度及黏性系數；Δx為壁面第1 層網格高度。

分別選取一方程SA 和二方程SST 湍流模型，對3 套不同密度網格進行數值模擬，所得阻力系數值與對應試驗值的對比見表3。結果表明，采用2 種不同的湍流模型時，計算結果與試驗值均相差很小，其中一方程SA 模型的計算結果和試驗值更為接近。此外，3 套網格的計算結果相差不大，并且隨著壁面第1 層網格高度的減小，阻力系數逐漸向試驗結果靠近。綜合考慮了計算量及計算精度之后，選取中等密度網格（即Δx=0.003 mm）及一方程SA 湍流模型開展數值模擬。

3 計算結果與分析

本節研究基于C 構型（不同減阻盤構型如圖6 所示）、桿件長度與鈍頭體直徑之比L/D=1.0、減阻盤直徑與鈍頭體直徑之比d/D=0.25 的自適應減阻盤，高超聲速飛行器運動方程見式（5），初始攻角取0°，振蕩幅值為20°，振蕩周期為125 s。

圖6 5 種不同構型減阻盤示意圖Fig. 6 Schematic diagrams of aerodisks of five different configurations

3.1 自適應減阻盤作用機理及動態流場

不同于固定式減阻盤，在鈍頭體運動過程中，自適應減阻盤始終保持對準來流方向，二者作用方式對比如圖7 所示。

圖7 固定式減阻盤及自適應減阻盤作用方式示意圖Fig. 7 Principle of fixed and self-aligned aerodisks

為分析自適應減阻盤對動態機動過程中飛行器周圍流場演化情況的影響，圖8 給出了高超聲速飛行器在1/2 周期內的對稱面瞬態密度云圖以及流線示意圖，其中流線采用流場馬赫數進行染色，可以發現，回流區、再附激波以及弓形激波共同決定著飛行器的阻力特性。如圖8（a）所示，減阻桿下方回流區（以下簡稱“下回流區”）主要受再附激波和弓形激波相對位置的影響。隨著俯仰角增加，再附激波向外偏折，在14.1°俯仰角附近已經和弓形激波共線，此時下回流區不再受到激波的限制，流動向鈍頭體后方擴散。盡管如此，在弓形激波的保護作用下，壁面免于來流直接沖擊，因此壓力依然處于較低水平。

圖8 加裝自適應減阻盤飛行器對稱面瞬態密度云圖及流線示意圖Fig. 8 Transient density contours and streamlines of aircraft symmetry plane with self-aligned aerodisks

對于減阻桿上方回流區（以下簡稱“上回流區”），在自適應減阻盤的作用下，其面積變化較小。0°俯仰角和7.7°俯仰角時，上回流區面積較大，鈍體上表面整體處于上回流區包覆之中，壁面壓力水平基本不變。在14.1°俯仰角時，受上壁面推動，上回流區跨過減阻桿向下運動；此時弓形激波向下偏折，導致再附點前移，鈍頭體肩部附近壁面壓力上升。自適應減阻盤的存在使得弓形激波形狀得以保持，在20°俯仰角時，減阻盤仍然可以將再附點控制在鈍體肩部附近。20°俯仰角之后，回流區又開始逐漸恢復原狀。

3.2 自適應減阻盤和固定式減阻盤對比

當采用2 種不同的減阻方式時，高超聲速飛行器在機動過程中的阻力系數時間歷程曲線如圖9 所示。在圖示的2 個振蕩周期中，采用自適應減阻盤方法之后，80%以上的大部分時間范圍內，飛行器阻力系數大幅下降。而在0°俯仰角附近，自適應減阻盤的阻力系數略大于固定式減阻盤，此時二者相對于基準飛行器（圖1）的減阻率均比較高。

圖9 阻力系數時間歷程曲線對比Fig. 9 Comparison of drag coefficients varying with time

圖10和圖11 分別為飛行器在20°俯仰角瞬時的回流區示意圖和對稱面流線圖，其中圖10 沿軸向對流場做剖面處理，各剖面均給出了馬赫數云圖，從中可以觀察到回流區的大致形狀。

圖10 20°俯仰角瞬時回流區示意圖Fig. 10 Recirculation zone at pitch angle of 20°

結合圖10 和圖11 可以看出，當高超聲速飛行器機動至20°俯仰角時，由于來流可以直接越過桿件，導致固定式減阻盤重構流場作用下降。而自適應減阻盤始終對準來流，因此桿件上下兩側流場對稱性更強，在回流區的保護下，有效降低了飛行器所承受的壓差阻力。

圖11 20°俯仰角瞬時對稱面流線圖Fig. 11 Streamlines of symmetry plane at pitch angle of 20°

3.3 不同參數對自適應減阻盤動態減阻效果的影響

本節采用控制變量的研究方法，通過非定常數值模擬，分別改變減阻盤構型、桿件長度以及減阻盤直徑等參數，探討自適應減阻盤在高超聲速飛行器動態機動過程中減阻效果的變化情況。飛行器的具體運動形式見1.2.2 節。

3.3.1 減阻盤構型

首先針對安裝不同形狀減阻盤的高超聲速飛行器，對比研究其在動態機動過程中的阻力特性。圖12 為不同減阻盤構型對應的飛行器阻力系數和升阻比隨攻角變化曲線。

如圖12（a）所示，構型A 在動態機動過程中阻力系數變化較大。其中，當飛行器“抬頭”（攻角：0°→20°、0°→-20°）時，阻力系數較小；“低頭”（攻角：20°→0°、-20°→0°）時，阻力系數相對較大，最大、最小阻力系數均出現在±6°俯仰角左右。另外4 種減阻盤構型對應的飛行器在“低頭”時的阻力系數同樣要高于“抬頭”過程，不過差別較小，最大阻力系數出現在±20°俯仰角附近。按照飛行器阻力系數大小，5種減阻盤構型可排序為：構型A＞構型D＞構型E＞構型B＞構型C。

圖12（b）為飛行器的升阻比變化曲線，其中，T表示飛行器的俯仰振蕩周期。在0°俯仰角（0、T/2、T）附近，由于飛行器上下對稱，其升力系數接近0，5 種構型的差別較小。不同形狀減阻盤的主要差別在20°俯仰角（T/4、3T/4）左右，構型C與構型E 的升阻比十分接近，并且小于其余三者，另外3 種構型按升阻比大小可排序為：構型B＜構型D＜構型A。

圖12 不同減阻盤構型對應的阻力系數及升阻比曲線Fig. 12 Drag coefficients and lift-to-drag ratios of different aerodisk configurations

5種不同減阻盤構型作用下0°俯仰角瞬時鈍頭體壁面壓力系數分布如圖13 所示，其中y＞0 部分代表軸線上方壁面，y＜0部分表示軸線下方壁面。

顯然，上下壁面存在一定的壓強差，在0°俯仰角時為飛行器提供了微弱的升力。圖13 中同時給出了基準飛行器在0°俯仰角時的壁面壓力分布，可以看出，在回流區覆蓋之下，鈍體頂點附近壁面壓強大幅度下降。在鈍頭體肩部，受弓形激波和流動再附等共同作用，壁面壓強遠高于駐點處。相對基準飛行器，構型C、構型B 以及構型E 依然能夠降低肩部壁面壓力，說明這3 種減阻盤構型的回流區面積較大，可以覆蓋到鈍體肩部，其中構型C 作用效果最佳。

圖13 0°俯仰角瞬時鈍頭體壁面壓力系數分布Fig. 13 Pressure coefficients of wall at pitch angle of 0°

隨著俯仰角增加，鈍頭體迎風面壓強逐漸上升，背風面由于處于弓形激波之后，壓力系數相對較低。20°俯仰角瞬時鈍頭體壁面壓力情況如圖14所示，其中構型A 壓力系數分布與基準飛行器差別較小，為0.4左右。在其余4種頭部構型作用下，背風面壓力系數可以下降至0.2 以下，下降率可以達到50%以上，按壓力系數大小排序為：構型C≈構型B＜構型E＜構型D。由于流動再附等原因，壁面壓力系數峰值出現在迎風面鈍頭體肩部附近；后體壓力系數相對較小，且基本不受減阻盤形狀影響。

圖14 20°俯仰角瞬時鈍頭體壁面壓力系數分布Fig. 14 Pressure coefficients of wall at pitch angle of 20°

3.3.2 減阻桿長度

圖15 給出的是加裝不同長度減阻桿的高超聲速飛行器在一個周期內的阻力系數變化情況。

圖15 阻力系數時間歷程曲線Fig. 15 Drag coefficients varying with time

結果顯示，在整個振蕩周期內，飛行器的阻力系數均隨著減阻桿長度增加而下降。其中當L/D從0.5 上升至1.0 時，阻力系數下降幅度最大，而L/D=2.0 和2.5 兩種構型之間差距相對較小。

由圖16 可知，在其余參數相同的情況下，隨桿件長度增加，鈍頭體前方壁面壓力水平明顯下降，壓差阻力逐漸減小。當L/D從0.5 增加至1.0 時，鈍頭體前方回流區尺寸逐漸增加，因此圖15 中L/D=0.5 和L/D=1.0 構型對應飛行器的阻力系數差別較大。此外，隨著L/D繼續增大，減阻桿頭部附近開始出現二次小回流區（L/D=1.5，2.0，2.5），不過由于二次回流區距離鈍體壁面較遠，對壁面壓力影響較小。鈍體前方大回流區面積也有增加，但變化較小，其中L/D=2.0 和2.5 時較為接近，因此L/D=1.0 之后加裝不同長度桿件對飛行器阻力系數的影響不大。

圖16 0°攻角瞬時流線及壁面壓力云圖Fig. 16 Streamlines and wall pressure contours under angle of attack of 0°

3.3.3 減阻盤直徑

當減阻盤直徑不同時，高超聲速飛行器在俯仰振蕩過程中的阻力系數變化情況如圖17 所示。不難發現，增大減阻盤直徑有利于改善飛行器機動過程中的阻力特性，不過隨著減阻盤尺寸增加，這種改善作用逐漸趨于飽和，即當減阻盤直徑為0.312 5D和0.375D時，飛行器對應阻力系數曲線逐漸靠近。

圖17 阻力系數隨攻角變化曲線Fig. 17 Drag coefficients varying with angle of attack

此外，d/D=0.125 時，最大阻力系數出現在±10°俯仰角附近，最小阻力系數出現在±12°俯仰角左右。其余4 種尺寸減阻盤對應的阻力系數的最大值和最小值則分別對應±20°及0°俯仰角。

如圖18 所示，在減阻盤肩部，高速氣流的方向會發生偏折，其上方形成一個膨脹扇面，Kharati-Koopaee 等［27］詳細解釋了這一現象。隨著減阻盤直徑增加，膨脹扇尺寸以及傾角隨之增大，分離激波向外擴展，弓形激波波后影響范圍擴大，因此不難理解“隨減阻盤尺寸增加，阻力系數逐漸下降”的變化規律。

圖18 不同俯仰角瞬時壓力等值線圖Fig. 18 Transient pressure contours at different pitch angles

觀察0°攻角時的等壓力線云圖可以發現，由于俯仰振蕩過程中存在遲滯現象，流場關于減阻盤軸線并不完全對稱。當d/D=0.312 5 時，分離激波和再附激波的位置逐漸趨于共線，此時回流區的面積基本不再大幅變化，因此圖17中d/D=0.312 5和d/D=0.375 對應的減阻盤構型在0°俯仰角附近阻力系數曲線接近重合。

4 結 論

針對自適應減阻盤，通過非定常數值模擬，首先探討了其作用機理及對動態流場演化的影響，然后對比了固定式減阻盤以及自適應減阻盤的差異，并詳細分析了高超聲速飛行器的動態阻力特性隨減阻盤構型、桿件長度以及減阻盤直徑等參數的變化規律。主要結論如下：

1）當存在來流攻角時，自適應減阻盤依然可以發揮其流場重構的作用，有效地解決減阻率急劇下降的問題，有利于機動飛行。

2）相對于傳統固定式減阻盤，在強迫俯仰振蕩過程中，采用自適應方法后，80%以上時間中飛行器阻力系數大幅下降。并且隨著俯仰角增大，自適應減阻盤相對于固定式減阻盤的優勢逐漸增加。

3）在減阻盤構型、桿件長度以及減阻盤直徑等參數的研究中發現，構型C 減阻效果最佳，高超聲速飛行器的阻力系數隨著桿件長度增加而下降，并且增大減阻盤直徑有利于改善飛行器在機動過程中的阻力特性。不過隨著減阻盤尺寸增加，該改善作用逐漸趨于飽和。

致 謝

感謝國防科技創新研究院的張來平老師和常興華老師提供的指導和幫助。