類X-37B航天器氣動力天地相關性數值模擬

馬率，張露，劉釩，孫俊峰，崔興達

中國空氣動力研究與發展中心 計算空氣動力研究所，綿陽 621000

跨大氣層飛行器氣動特性的復雜性在于飛行走廊(即飛行軌道)的參數非常復雜，包括離地高度、馬赫數、迎角、側滑角等流動參數和副翼、升降舵、方向舵、襟翼偏角等控制參數，單純依靠風洞和飛行試驗的手段從耗費、時間等角度來說遠不能滿足需求。而航天飛行器相比航空飛行器的風洞試驗更加困難、外形較簡單更適合CFD計算，CFD高效批量生產數據的特點可以對試驗數據進行補充，尤其是對地面數據向真實飛行條件的修正與外推的天地相關研究工作具有重要意義[1]。

X-37B是波音公司為競標美國空軍“軌道飛行器項目”而研制的無人且可重復使用的跨大氣層在軌飛行器，該機由火箭發射進入太空，是第一種既能在環繞地球衛星軌道上飛行又能自主重返大氣層并最終著陸的航天飛行器，截至2014年10月17日，X-37B(OTV3)已完成連續飛行超過674天[2]。

波音公司設計的X-37B特點有：機翼由細長邊條翼和后部短小三角翼組合而成雙三角翼，機翼后緣設計有全展長的襟副翼、機體上部有一對集方向舵和升降舵作用的V型尾翼、V尾中間有減速板，機體后帶有體襟翼，飛行器腹部為平面的翼身組合體布局(也稱升力體布局)[3]，見圖1。X-37B的氣動特性研究主要包括風洞試驗(見圖2)、數值模擬和飛行試驗3個方面，根據波音公司的報告，X-37的風洞試驗數據只占了研發與創建用于再入飛行的氣動數據庫的19%，剩下的絕大部分采用數值計算方法，只有小部分是用飛行試驗來驗證地面試驗和數值計算結果，彌補地面試驗的不足[4]。鑒于X-37B項目的成功，國內也先后采用了CFD技術對此類飛行器開展了布局概念設計[5]及氣動特性研究[6]。

圖1 X-37B的氣動控制面[3]Fig.1 Aero-surfaces of X-37B[3]

圖2 AEDC風洞中的X-37B氣動試驗Fig.2 Aerodynamics tests of X-37B in AEDC wind tunnels

為了研究類X-37B布局航天器風洞試驗與飛行數據的天地相關性問題，建立了類X-37B航天器相關模型，然后開展了計算模型網格的規模影響及修正研究，在此基礎上對比分析了雷諾數和試驗狀態支架干擾的影響，完成了類X-37B航天器氣動力天地差異性分析，以期為今后開展相關研究積累技術基礎。

1 計算模型及數值方法

參照已公布的X-37B外形信息及風洞試驗圖片，生成了類似X-37B的飛行器模型和帶支桿1∶24的縮比模型，見圖3。計算網格分為真實構型和帶尾支撐桿的縮比模型,真實構型的第1層網格距離物面的距離為5×10-6m，而縮比模型的第1層網格距離物面的距離相應為2.3×10-7m，這兩個距離在馬赫數范圍(Ma=0.4～8.0)、離地高度55 km以下都能達到0.3≤y+≤5,圍繞飛行器生成的是C型對接網格，在保持相同物面距離的條件下，網格規模在研究網格無關性時為1 000萬 ～8 000萬。在不同雷諾數和試驗狀態支架干擾的影響研究中采用的網格規模為1 000萬，有關網格規模對計算結果的影響將在下文做詳細說明。圖4給出了飛行器的相關計算網格。

圖3 類X-37B基本模型及帶支桿縮比模型Fig.3 X-37B-like base model and model with wind tunnel strut

圖4 類X-37B基本模型及帶支桿縮比模型計算網格Fig.4 Computational grid of X-37B-like base model and model with wind tunnel strut

計算采用筆者所在課題組自主研發的并行多塊結構化網格的數值計算PMB3D求解器[7-8]完成，該程序目前已經升級到了9.0版本，在國內航空航天領域多個復雜流動工程項目的計算模擬中得到了應用與驗證[9-12]。本文數值模擬應用了有限體積法求解Navier-Stokes方程，在慣性坐標系下的積分形式為

(1)

式中：狀態變量Q=[ρ,ρu,ρv,ρw,e]T用守恒變量表示，分別代表流體的密度、動量分量和總能；HI表示控制體邊界面上的對流通量；HV表示黏性通量；Ω表示控制體的體積；S表示控制體邊界面的面積；n為邊界面的外法向量。計算中黏性項采用中心差分，無黏項采用Roe平均迎風通量差分分裂格式離散，選取Condiff限制器來保證將網格單元格心物理量插值到界面處的精度及魯棒性，時間格式采用了隱式LU-SGS方法，湍流模型選取了考慮可壓縮修正的Menter’sk-ωSST兩方程模型，并運用了多重網格、殘值平均和局部時間步長、并行加速等高效計算方法。

2 網格規模影響及修正

2.1 網格規模差異對計算結果影響

保證足夠大的網格規模對提高計算結果可信度是十分重要的，這里采用從1 000萬、2 000萬、3 000萬、 5 000萬直至8 000萬網格規模，對類X-37B布局的飛行器網格無關性計算做深入研究。圖5給出了1 000萬和8 000萬網格規模的表面網格比較。以Ma=0.4，迎角α=0°，4°，8°為例，研究網格對縱向氣動特性的影響。圖6給出了氣動系數隨網格規模的變化趨勢，其中，N為網格格心點數，1/N2/3為AIAA阻力預測會議(Drag Prediction Workshop, DPW)[13-14]歸納的研究氣動力預測精度與網格變化關系的參數。

圖5 不同網格規模的物面網格比較Fig.5 Comparison of different sizes of surface meshes

由圖6可見，軸向力系數Cx隨1/N2/3減小基本上線性下降，而法向力系數CN、俯仰力矩系數Cmz隨1/N2/3減小變化很小。圖6中縱坐標數值是軸向力系數隨1/N2/3的線性關系推出的采用PMB3D軟件計算Ma=0.4，α=0°、4°、8°狀態得到的軸向力系數Cx理論網格收斂值。為進一步研究亞跨超范圍內采用1 000萬網格和8 000萬網格對計算結果的影響，圖7給出了1 000萬網格相對8 000萬網格計算0°迎角時軸向力系數Cx和摩阻系數Cf隨馬赫數變化曲線，可見在亞跨聲速時，兩套網格的軸向力計算結果差別較大，在馬赫數0.4時誤差最大為32.4%，在馬赫數1.05以后誤差不到5%，而兩套網格的摩阻系數計算結果差別較小，最大誤差不到6%。另外兩套網格在亞跨超范圍內法向力系數的誤差最大為2.05%， 俯仰力矩系數的誤差最大為1.95%，縱向壓心的誤差不超過0.7%機身長度。

圖6 氣動系數隨網格數的變化趨勢(Ma=0.4，β=0°,Re=7.91×107)Fig.6 Varying trends of aerodynamic coefficients with grids (Ma=0.4，β=0°,Re=7.91×107)

圖7 α=0°時網格數對Cx和Cf的影響及相對誤差Fig.7 Influence of grid numbers on Cx and Cf and relative error at α=0°

綜上所述，亞聲速來流計算狀態網格規模對由壓差產生的軸向力影響較大，對法向力系數、俯仰力矩系數和縱向壓心影響很小。

2.2 網格影響修正

2.1節研究表明，采用密網格計算的軸向力系數明顯比稀網格計算的合理可信，但采用密網格耗費的計算機資源龐大，比如8 000萬網格比1 000萬 網格的計算機資源占用大出近一個數量級。因此，研究稀網格和密網格計算結果差異的規律性，并運用規律將稀網格的計算結果修正到理論網格收斂值上來是十分有意義的。

研究軸向力系數與迎角的變化趨勢可以發現1 000萬網格計算的軸向力系數與8 000萬網格計算的軸向力系數之間的差異基本上與迎角無關，前體軸向力系數和底部軸向力系數也存在類似的關系。因此可以用表1中0°迎角8 000萬網格與1 000萬網格計算的軸向力系數之差、前體軸向力系數之差和底部軸向力系數之差作為修正量，來修正1 000萬網格計算結果中的軸向力系數、前體軸向力系數和底部軸向力系數，升力系數、阻力系數和升阻比根據軸向力系數和法向力系數對不同迎角的投影而重新計算，其他氣動系數由于隨網格規模變化基本上在2%以內，不做修正。

由于篇幅所限，圖8只給出了在亞跨聲速范圍1 000萬網格、8 000萬網格和1 000萬網格修正計算得到的軸向力系數Cx、底部軸向力系數Cxb。由圖可見，1 000萬網格修正后的結果與8 000萬網格的計算結果基本重合，說明以上修正方法是可信的，在工程上也是非常實用的。

表1 α=0°時網格影響修正量Table 1 Grid influence correction at α=0°

圖8 修正后氣動特性與采用1 000萬、8 000萬網格計算結果比較(β=0°，δ=0°)Fig.8 Comparison between correction aerodynamic coefficients and calculation results with 10 million and 80 million grids (β=0°，δ=0°)

3 高空飛行與風洞試驗氣動差異

3.1 雷諾數對基本外形的氣動特性影響

分別按飛行雷諾數和風洞試驗雷諾數計算類X-37B飛行器無舵偏無側滑的氣動特性，這里氣動特性計算仍然使用1 000萬網格并采用了第2節中的網格修正方法來修正相關的氣動力結果，高空飛行與風洞試驗雷諾數的比較見圖9,可見高空飛行雷諾數都要大于風洞雷諾數[15]。

飛行雷諾數、風洞試驗雷諾數和風洞試驗雷諾數帶支桿三者的氣動特性部分計算結果見圖10。 結果表明在亞跨聲速來流時，雷諾數增大使法向力系數CN、升力系數CL和俯仰力矩系數Cmz的絕對值略有增大，隨著馬赫數增加，雷諾數影響逐漸減小，在高超聲速時雷諾數影響量又開始微幅增大。在亞跨聲速來流時雷諾數增大使軸向力系數Cx和阻力系數CD明顯減小，隨著馬赫數增加雷諾數的影響量減小，在高超聲速時雷諾數影響量又開始略有增大。在所計算的馬赫數范圍內除0°迎角外雷諾數增大使壓心略有后移，但最大后移量不超過機身長度的0.5%。在亞跨聲速范圍內由于雷諾數對軸向力和阻力系數影響較大，因此雷諾數對此范圍內的升阻比影響較大，雷諾數增大使升阻比增大，最大變化在Ma=0.4、α=10° 處，增大量為0.586 22。

圖9 高空飛行雷諾數與地面風洞試驗雷諾數比較Fig.9 Comparison of Reynolds number of high altitude flight and ground wind tunnel test

在真實風洞試驗中，測量得到的模型阻力是很難將摩擦阻力和壓差阻力分開并區分準確的，而CFD的優勢之一就是能將計算得到的阻力分解為壓差阻力和摩擦阻力，這就為飛行器的阻力特性分析提供定性和定量的優勢。為研究雷諾數對阻力系數的影響，采用CFD方法對真實飛行狀態和試驗工況狀態的0°迎角阻力系數CD進行了數值模擬，并將其分解為壓差阻力系數CDp和摩擦阻力系數CDf，以各自飛行雷諾數的Cx作為歸一化系數，它們的構成見表2，可以看到在亞聲速階段，摩阻能占到阻力的近30%，隨著馬赫數的增大所占比重急劇減小，到Ma=7.0時已不超過阻力的5%。另外飛行和風洞試驗的雷諾數差異對壓差阻力和摩擦阻力都有一定的影響。根據附體流動的一般規律，雷諾數增大使摩擦阻力系數減小，由于Ma=0.4 飛行雷諾數與風洞試驗雷諾數差異最大，因此在Ma=0.4時雷諾數對摩擦阻力系數影響最大。但在Ma=5.0，7.0時，雷諾數的增大并沒有使摩擦阻力系數減小，這是因為按風洞雷諾數計算時物面分離導致的，由于風洞雷諾數的邊界層厚度相比飛行雷諾數的變厚，機身上的流動頂著逆壓梯度前進的能力變弱，風洞雷諾數相比飛行雷諾數在機身上更容易分離，而分離區內靠近壁面的流動實際上變為回流或逆流，導致物面的摩擦力相比未分離時反向，與飛行器總摩擦力方向相反，產生了類似附加推力的作用。圖11給出的是Ma=3.0、Ma=7.0按飛行雷諾數和風洞試驗雷諾數計算的對稱面流場等馬赫數圖，在Ma=3.0時按兩個雷諾數計算的流場沒有分離，而在Ma=7.0時可以看到按風洞雷諾數計算的流場機身前體已經明顯分離。

進一步把壓差阻力CDp分解為前體壓差阻力CDfp和底阻CDd，由表2可見，在亞跨聲速時，雷諾數對前體壓差阻力系數有一定的影響，一方面雷諾數對底部分離流動的影響會傳到上游流動，從而改變上游流動的壓強分布；另一方面雷諾數對附面層厚度的影響會導致跨聲速激波位置的變化從而引起波阻改變。而在超聲速時，由于雷諾數對底部分離流動的影響不會傳到上游流動，所以雷諾數對前體壓差阻力系數的影響非常小。在高超聲速時，由于飛行器前體流動分離，使雷諾數的影響分析更加復雜。

圖10 雷諾數及有/無支桿對氣動特性影響比較(β=0°，δ=0°)Fig.10 Comparison of influence of Reynolds number on aerodynamic characteristics with and without support sting (β=0°，δ=0°)

表2 迎角0°時基于CFD計算結果的雷諾數對阻力系數的影響分析(歸一化)Table 2 Re influence correction of CD base on CFD results(normalization) at α=0°

圖11 雷諾數對流場的影響(α=0°)Fig.11 Effects of Reynolds number on flow fields (α=0°)

3.2 有/無支桿對基本外形的氣動特性影響

模型支桿的氣動干擾問題自風洞試驗誕生以來就成為阻礙模型風洞試驗精準度提高的難題。從更廣義的角度上看,由于模型支桿破壞了試驗模型的有效外形和支桿本身的繞流特性，干擾了試驗模型的流場，這兩方面的因素改變了按飛行器外形幾何相似縮小的試驗模型的空氣動力學特性[16]。

目前風洞試驗結果都是扣除了帶尾支桿影響在內的底部阻力或底部軸向力，這就是試驗中常被提及的“前體阻力”或“前體軸向力”，扣除了底阻的帶尾支桿支撐的模型風洞試驗,很難準確給出依賴于總阻力或總軸向力的試驗模型的配平迎角和相關力矩，因此對支桿影響的修正在風洞試驗中不可忽視。

本節按風洞試驗雷諾數分別計算了有/無支桿情況下的氣動特性[17-18]，并通過比較兩者之間的差異分析支桿的影響。支桿外形盡量模擬風洞試驗的支桿，見圖3，值得指出的是，為了使帶支桿外形的底阻計算與對應的無支桿外形保持一致，真實風洞中對應的風洞試驗底阻計算一般是通過支桿在底部附近布置的n個采樣測壓點的壓強平均獲得的，具體為

(2)

式中：CpA、CpB…CpN分別為支桿在底部1～n個測壓點的壓力系數；Sb為無支桿外形的完整底部面積；Sref為計算參考面積。這樣帶支桿外形軸向力系數計算值就修改為

(3)

由圖10可見，支桿對升力系數CL和俯仰力矩系數Cmz影響不明顯，對阻力系數CD影響較大，在亞跨聲速來流時，支桿的存在使軸向力系數和阻力系數減小，隨著馬赫數增大，支桿的影響量逐漸減小，特別是超聲速和高超聲速，由于下游流動不能向上傳遞，因此除0°迎角支桿的存在使升阻比增大外,其他狀態下支桿的影響量很小。

圖12 底部軸向力系數和前體軸向力系數在有/無支桿情況下比較 (β=0°，δ=0°)Fig.12 Comparison of Cxb and Cxf with and without support sting (β=0°，δ=0°)

圖13 底部測壓點示意圖Fig.13 Diagram of bottom pressure measuring points

圖14 有/無支桿外形底部壓強分布比較Fig.14 Comparison of bottom pressure distribution with and without support sting

4 結 論

本文采用數值模擬方法研究了網格規模對類X-37B布局航天器數值計算結果的影響，獲得了不同網格規模對縱向氣動特性的影響規律，探討了實用于工程的網格規模影響修正方法，通過對該布局飛行器高空飛行和風洞試驗狀態數值計算數據分析和研究，得到以下結論：

1) 不同網格規模對類X-37B氣動布局航天器的亞聲速軸向力系數影響較大，對法向力系數、俯仰力矩系數和縱向壓心影響不大。

2) 不同網格規模對軸向力系數的影響基本上與迎角無關，本文提出的網格規模影響修正方法對類X-37B氣動布局飛行器的計算結果修正是可信的，在工程上也是非常實用的。

3) 由于風洞試驗模型縮比較大，特別是在亞跨聲速風洞試驗雷諾數與高空飛行雷諾數差一個量級，因此雷諾數對類X-37B氣動布局飛行器氣動特性特別是軸向力系數、阻力系數和升阻比有較大的影響。

4) 由于風洞試驗狀態支桿的存在，在亞跨聲速來流時對類X-37B氣動布局飛行器底阻影響較大，因而對此類飛行器軸向力系數、阻力系數和升阻比有較大的影響。

5) 在本文研究范圍內，雷諾數和支桿對類X-37B氣動布局飛行器的軸向力系數影響趨勢相反，從而使高空飛行與風洞試驗狀態軸向力系數的差異縮小。