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

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

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

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

波音公司設(shè)計的X-37B特點有：機翼由細長邊條翼和后部短小三角翼組合而成雙三角翼，機翼后緣設(shè)計有全展長的襟副翼、機體上部有一對集方向舵和升降舵作用的V型尾翼、V尾中間有減速板，機體后帶有體襟翼，飛行器腹部為平面的翼身組合體布局(也稱升力體布局)[3]，見圖1。X-37B的氣動特性研究主要包括風洞試驗(見圖2)、數(shù)值模擬和飛行試驗3個方面，根據(jù)波音公司的報告，X-37的風洞試驗數(shù)據(jù)只占了研發(fā)與創(chuàng)建用于再入飛行的氣動數(shù)據(jù)庫的19%，剩下的絕大部分采用數(shù)值計算方法，只有小部分是用飛行試驗來驗證地面試驗和數(shù)值計算結(jié)果，彌補地面試驗的不足[4]。鑒于X-37B項目的成功，國內(nèi)也先后采用了CFD技術(shù)對此類飛行器開展了布局概念設(shè)計[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布局航天器風洞試驗與飛行數(shù)據(jù)的天地相關(guān)性問題，建立了類X-37B航天器相關(guān)模型，然后開展了計算模型網(wǎng)格的規(guī)模影響及修正研究，在此基礎(chǔ)上對比分析了雷諾數(shù)和試驗狀態(tài)支架干擾的影響，完成了類X-37B航天器氣動力天地差異性分析，以期為今后開展相關(guān)研究積累技術(shù)基礎(chǔ)。

1 計算模型及數(shù)值方法

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

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

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

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

(1)

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

2 網(wǎng)格規(guī)模影響及修正

2.1 網(wǎng)格規(guī)模差異對計算結(jié)果影響

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

圖5 不同網(wǎng)格規(guī)模的物面網(wǎng)格比較Fig.5 Comparison of different sizes of surface meshes

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

圖6 氣動系數(shù)隨網(wǎng)格數(shù)的變化趨勢(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°時網(wǎng)格數(shù)對Cx和Cf的影響及相對誤差Fig.7 Influence of grid numbers on Cx and Cf and relative error at α=0°

綜上所述，亞聲速來流計算狀態(tài)網(wǎng)格規(guī)模對由壓差產(chǎn)生的軸向力影響較大，對法向力系數(shù)、俯仰力矩系數(shù)和縱向壓心影響很小。

2.2 網(wǎng)格影響修正

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

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

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

表1 α=0°時網(wǎng)格影響修正量Table 1 Grid influence correction at α=0°

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

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

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

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

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

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

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

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

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

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

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

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

模型支桿的氣動干擾問題自風洞試驗誕生以來就成為阻礙模型風洞試驗精準度提高的難題。從更廣義的角度上看,由于模型支桿破壞了試驗?zāi)Ｐ偷挠行庑魏椭U本身的繞流特性，干擾了試驗?zāi)Ｐ偷牧鲌?，這兩方面的因素改變了按飛行器外形幾何相似縮小的試驗?zāi)Ｐ偷目諝鈩恿W特性[16]。

目前風洞試驗結(jié)果都是扣除了帶尾支桿影響在內(nèi)的底部阻力或底部軸向力，這就是試驗中常被提及的“前體阻力”或“前體軸向力”，扣除了底阻的帶尾支桿支撐的模型風洞試驗,很難準確給出依賴于總阻力或總軸向力的試驗?zāi)Ｐ偷呐淦接呛拖嚓P(guān)力矩，因此對支桿影響的修正在風洞試驗中不可忽視。

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

(2)

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

(3)

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

圖12 底部軸向力系數(shù)和前體軸向力系數(shù)在有/無支桿情況下比較 (β=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 結(jié) 論

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

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

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

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

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

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