基于NACA0012翼型的超燃沖壓高超聲速飛行器的氣動熱預(yù)示

汪 鈺,李 芃

（銅陵學(xué)院，安徽 銅陵 244000）

0 引言

高速氣流在穿過激波時會被壓縮，由此導(dǎo)致減速和溫度急劇升高，產(chǎn)生的高溫環(huán)境會引起飛行器燒蝕、干擾內(nèi)部設(shè)備的正常工作。因此，對飛行器氣動熱的準(zhǔn)確預(yù)示至關(guān)重要。風(fēng)洞試驗(yàn)、工程估算和計(jì)算流體力學(xué)（CFD）是氣動熱預(yù)示的常用方法。第一種方法耗時長、成本高，第二種方法的估算精度有限。隨著計(jì)算機(jī)技術(shù)的迅速發(fā)展，CFD已成為氣動熱預(yù)示研究的主要方法［1-3］。NA‐CA0012是研究飛行器外流場特性的典型翼型，它是由美國國家航空咨詢委員會（National Advisory Committee for Aeronautics簡稱NACA）設(shè)計(jì)的。目前基于NACA0012的氣動熱研究主要集中在兩個方面：（1）翼型的幾何變化對溫度或溫度系數(shù)的影響［4-5］；（2）溫度對翼型表面結(jié)冰特性的影響［6-7］。后緣形狀是影響預(yù)示精度的重要因素［8］，根據(jù)建模方法的不同，NACA0012包括尖后緣和鈍后緣兩種。本研究設(shè)計(jì)了兩種翼型外流場，通過與參考數(shù)據(jù)的對比，分析了不同參數(shù)對仿真預(yù)示的影響，給出了高超聲速飛行條件下的最佳參數(shù)配置，并進(jìn)行了相關(guān)的氣動熱預(yù)示，進(jìn)一步豐富和完善了有關(guān)NA‐CA0012的研究成果。

1 仿真基礎(chǔ)

1.1 超燃沖壓高超聲速飛行器

超燃沖壓高超聲速飛行器是指以超燃沖壓發(fā)動機(jī)為主要動力的巡航飛行器，此類發(fā)動機(jī)只有在飛行速度達(dá)到一定值時才能啟動，因此通常采用組合推進(jìn)方案。如以火箭為基礎(chǔ)的聯(lián)合循環(huán)（RBCC）、以渦輪為基礎(chǔ)的聯(lián)合循環(huán)（TBCC）和固體火箭雙模推進(jìn)。其中，RBCC充分發(fā)揮了火箭發(fā)動機(jī)大推重比和沖壓發(fā)動機(jī)高比沖的優(yōu)勢，是未來高超聲速飛行器和單級再利用飛行器的理想動力系統(tǒng)。

為了正確地預(yù)示RBCC高超聲速飛行器的氣動熱，需要分析飛行器的飛行路徑并確定自由來流的初始條件。如圖1所示，RBCC飛行器的飛行路徑由爬升段、巡航段和攻擊段組成。飛行器首先由載機(jī)運(yùn)送到指定地點(diǎn)，啟動固體火箭發(fā)動機(jī)將飛行器從低速模態(tài)帶入亞燃沖壓模態(tài)。該階段推進(jìn)效率最高，高度和速度的增量最大。最后進(jìn)入超燃沖壓模態(tài)，當(dāng)飛行器平穩(wěn)爬升到目標(biāo)高度和馬赫數(shù)后，進(jìn)入保持等高等速飛行的巡航段，其飛行距離約占總航程的90%，是飛行路徑的主要階段。到達(dá)目標(biāo)空域后進(jìn)入攻擊段，對目標(biāo)實(shí)施精準(zhǔn)打擊。本文提出了適用于巡航段的氣動熱預(yù)示方案，巡航段的飛行高度約為30km，自由來流的初始條件如下：馬赫數(shù)M為6.5，大氣壓力Pt為1179Pa，大氣溫度T t為226.65K，大氣密度ρ為0.0181 kg/m3，大氣動力黏度μ為1.475e-05 kg/（m·s）。

圖1 RBCC超燃沖壓高超聲速飛行器的飛行軌跡

1.2 計(jì)算域外流場

如圖2所示，基于NACA4 airfoil generator和air‐foil tools設(shè)計(jì)了尖后緣和鈍后緣兩種NACA0012翼型并繪制了計(jì)算域外流場，其中翼型的特征長度L=1m。對外流場的數(shù)值模擬，一般建議遠(yuǎn)場距離為翼型特征長度的12-20倍。建立遠(yuǎn)場距離為20m的C型網(wǎng)格拓?fù)湟宰钚』诰W(wǎng)格的偏度。尖后緣采用ICEM CFD進(jìn)行網(wǎng)格劃分時可直接折疊后緣處的塊，如圖2（a）所示。鈍后緣保留后緣處的塊，如圖2（b）所示。

圖2 （a）尖后緣計(jì)算域外流場；（b）鈍后緣計(jì)算域外流場

1.3 湍流模型

對于涉及可壓縮和熱傳遞的流動，需要求解質(zhì)量、動量和能量守恒方程。此外，還要求解湍流帶來的輸運(yùn)方程?？紤]到湍流方程無法通過直接數(shù)值模擬（DNS）求解，采用了目前應(yīng)用最為廣泛的雷諾平均可壓縮NS方程（RANS）。K-epsilon、SST komega和Spalart-Allmaras是被廣泛應(yīng)用且相對準(zhǔn)確的RANS湍流模型［9-11］。K-epsilon模型主要用于工業(yè)CFD仿真，該模型存在對邊界層分離和逆壓梯度不敏感的問題，所以在翼型外流場的應(yīng)用并不廣泛。Spalart-Allmaras是為航空應(yīng)用設(shè)計(jì)的單方程模型，能夠?yàn)槭苣鎵禾荻扔绊懙倪吔鐚臃抡嫣峁┝己玫慕Y(jié)果。SST k-omega是基于標(biāo)準(zhǔn)k-omega模型設(shè)計(jì)的雙方程模型，其對自由來流的初始條件不敏感，且經(jīng)過校準(zhǔn)，能夠準(zhǔn)確計(jì)算光滑表面的流動分離。因此采用雙方程的SST k-omega湍流模型。

2 驗(yàn)證測試

2.1 網(wǎng)格策略

因?yàn)檠埠蕉蔚鸟R赫數(shù)M大于5，此時的計(jì)算域外流場處于高超聲速階段，所以采用參考文獻(xiàn)［12］中的數(shù)據(jù)來驗(yàn)證仿真參數(shù)。計(jì)算域外流場如圖3（a）所示，其中粗黑點(diǎn)為坐標(biāo)原點(diǎn)，a、b和c是輸入邊界INPUT boundary，e、f和g是壁面邊界WALL boundary，d是輸出邊界OUTPUT boundary。湍流流動受翼型壁面的影響，因此，近壁面區(qū)域的網(wǎng)格劃分十分重要，需要保證邊界層的第一層網(wǎng)格處于sub viscosity layer，關(guān)鍵在于計(jì)算第一層網(wǎng)格高度y H。相關(guān)計(jì)算公式如式（1）：

其中R e是雷諾數(shù)，ρ、U、μ和L分別是自由來流的密度、速度、動力黏度和翼型的特征長度。

其中y P是第一層網(wǎng)格中間位置到壁面的距離，y H=2y P。氣動熱預(yù)示對壁面邊界層的精度要求較高，一般建議近壁面網(wǎng)格y+≤1。驗(yàn)證測試的參數(shù)取值如下：R e=10e+06，L=1m，μ=4.46082e-5kg/（m·s），ρ=0.0247kg/m3，Pt=576Pa，Tt=81.2K，U=1806m/s，近壁面網(wǎng)格加密區(qū)如圖3（b）所示。為了保證數(shù)值仿真結(jié)果的可靠性進(jìn)行了網(wǎng)格無關(guān)性測試，如表1所示，測試了420 000、590 000和830 000三種網(wǎng)格數(shù)，P和T是翼型表面不同位置的壓力和溫度，830 000網(wǎng)格T/Tt的平均誤差率優(yōu)于590 000網(wǎng)格，但是后者的P/Pt的平均誤差率數(shù)值和總平均誤差率均優(yōu)于830 000網(wǎng)格，因此采用590 000網(wǎng)格。

圖3 （a）計(jì)算域外流場；（b）近壁面網(wǎng)格加密區(qū)

單元雷諾數(shù)［13］和近激波面的aspect ratio［14］是影響仿真準(zhǔn)確度的兩個重要網(wǎng)格參數(shù)。文獻(xiàn)［13］建議對二維圓柱體，其單元雷諾數(shù)應(yīng)不小于8。針對NACA0012翼型，對單元雷諾數(shù)分別取值4、8和16來尋找理想的數(shù)值，對應(yīng)的y H和y+分別為（3.5e-6m 0.0825）（8e-6m 0.175）和（1e-5m 0.35）?；诰W(wǎng)格總數(shù)和單元雷諾數(shù)不變的前提，再比較不同as‐pect ratio數(shù)值下的誤差率，來確認(rèn)最合適的網(wǎng)格參數(shù)，見表1。

表1 網(wǎng)格無關(guān)性

2.2 數(shù)值方法

表2給出了在高超聲速階段采用的數(shù)值仿真參數(shù)。求解器選用專為高速可壓縮流設(shè)計(jì)的密度基，通量采用ROE和AUSM類型，其仿真結(jié)果與試驗(yàn)數(shù)據(jù)較為接近。SST k-omega湍流模型應(yīng)用intensity and viscosity ratio turbulence method，其中的湍流強(qiáng)度（%）和湍流黏度比均為1。計(jì)算域的邊界條件如下：輸入邊界INPUT boundary類型為pressure far field，velocity specification method選 擇magnitude and direction，攻角為0°，自由來流X方向分量為cos（0°）=1、Y分量為sin（0°）=0。輸出邊界OUTPUT boundary類型為pressure outlet，壁面邊界WALL boundary類型為no-slip，等溫壁，初始值311K。梯度計(jì)算的常用方法有cell-based和node-based，第二種方法的準(zhǔn)確度更高，尤其是針對具有傾斜度的非結(jié)構(gòu)化網(wǎng)格，但是其對于硬件資源的要求更高?？紤]到生成的結(jié)構(gòu)化網(wǎng)格和有限計(jì)算資源，采用least-squares cell-based方案計(jì)算梯度，該方法可以提供與node-based相當(dāng)?shù)木惹矣?jì)算開銷較低。對于壓力，采用基于中心差分的second-order方案。

表2 數(shù)值仿真參數(shù)

2.3 數(shù)值結(jié)果分析

基于一種湍流模型、兩個通量類型和三個單元雷諾數(shù)，對每一種后緣翼型進(jìn)行如下6組參數(shù)配置的 仿 真：SST+AUSM+yplus0.35、SST+ROE+yplus0.35、SST+AUSM+yplus0.175、SST+ROE+yplus0.175、SST+AUSM+yplus0.0825、SST+ROE+yplus0.0825。

對尖后緣，上述參數(shù)配置下的P/Pt分布如圖4（a）所示，其中X軸為x/L，Y軸為P/Pt的數(shù)值。T/Tt的分布如圖4（b）所示。P/Pt和T/Tt的平均誤差率分別為（4.01%4.17%）（19.91%14.30%）（5.84%3.86%）（11.64%9.21%）（4.09%4.10%）和（12.96%12.88%）。如圖4（c）所示，SST+ROE+yplus0.35配置下的總平均仿真誤差率最低。接著基于該配置，對近激波面處網(wǎng)格的aspect ratio分別取值30、60、120、180、240和300再執(zhí)行6組仿真，如圖4（d）所示，對應(yīng)的P/Pt和T/Tt的平均誤差率分別為（5.06%10.05%）（4.01%4.17%）（4.74%3.98%）（5.47%3.10%）（4.71%3.14%）和（10.48%7.06%）。當(dāng)as‐pect ratio取值較小時，壓力和溫度的平均誤差率較大。之后隨著aspect ratio的增大，平均誤差率開始下降。當(dāng)aspect ratio等于240時可獲得最佳結(jié)果，總的平均誤差率約3.93%。

圖4 （a）尖后緣P/Pt的分布圖；（b）尖后緣T/Tt的分布圖；（c）尖后緣在不同參數(shù)配置下的平均誤差率；（d）尖后緣在不同aspect ratios下的平均誤差率

對鈍后緣，上述參數(shù)配置下的P/Pt分布如圖5（a）所示，其中X軸為x/L，Y軸為P/Pt的數(shù)值。類似的，T/Tt的分布如圖5（b）所示。P/Pt和T/Tt的平均誤差率分別為（5.67%3.35%）（4.22%4.14%）（7.94%3.08%）（8.91%6.70%）（5.27%3.36%）和（9.47%6.02%）。如圖5（c）所示，SST+AUSM+yplus0.35配置下的平均誤差率最低。接著基于該配置，對近激波面處網(wǎng)格的aspect ratio分別取值40、80、160、240、320和400再執(zhí)行6組仿真，如圖5（d）所示，對應(yīng)的P/Pt和T/Tt的平均誤差率分別為（19.64%7.75%）（4.22%4.14%）（3.67%3.28%）（5.91%10.03%）（3.06%2.83%）和（2.74%9.58%）。當(dāng)aspect ratio取值較小時，壓力和溫度的平均誤差率較大。之后隨著aspect ratio的增大，平均誤差率開始下降。當(dāng)as‐pect ratio等于320時可獲得最佳結(jié)果，總的平均誤差率約2.95%。

因此，對NACA0012翼型的后緣形狀，鈍后緣的性能要優(yōu)于尖后緣，前者較后者帶來的仿真精度提升了約24.97%，后者較前者導(dǎo)致仿真精度下降了約33.28%，后緣翼型選擇不當(dāng)會導(dǎo)致較大的仿真誤差。所以在高超聲速階段，建議采用的參數(shù)配置如下：鈍后緣翼型+y+取值0.35+單元雷諾數(shù)為16+近激波面aspect ratio為320+密度基求解器應(yīng)用AUSM類型的通量。如果計(jì)算域外流場選用了尖后緣翼型，那么建議采用如下參數(shù)配置：y+取值0.35+單元雷諾數(shù)為16+近激波面aspect ratio為240+密度基求解器應(yīng)用ROE類型的通量。

圖5（a）鈍后緣P/Pt的分布圖；（b）鈍后緣T/Tt的分布圖；（c）鈍后緣在不同參數(shù)配置下的平均誤差率；（d）鈍后緣在不同aspect ratios下的平均誤差率

3 氣動熱預(yù)示

基于NACA0012鈍后緣翼型建立了遠(yuǎn)場距離為20m的C型計(jì)算域外流場，網(wǎng)格y+取值0.35（對應(yīng)y H=1.6e-5m），單元雷諾數(shù)為16，近激波面aspect ra‐tio取值320。高超聲速飛行器巡航段的氣動熱預(yù)示方案仿真參數(shù)如表3所示，在巡航段飛行器保持6.5馬赫的恒定速度和30km的恒定高度，自由來流的初始壓力和溫度分別為1179Pa和226.65K。初始化選用基于拉普拉斯方程的hybrid非均勻方法，穩(wěn)態(tài)分析，計(jì)算迭代5000次。氣動熱溫度最終達(dá)到了2158K（1885℃），預(yù)示結(jié)果表明RBCC高超聲速飛行器在飛行過程中將承受極端的氣動熱問題。

表3 數(shù)值仿真參數(shù)

4 結(jié)論

基于NACA0012翼型，采用文獻(xiàn)［12］的初始條件和數(shù)據(jù)進(jìn)行了驗(yàn)證測試，分析了不同后緣形狀和網(wǎng)格參數(shù)（單元雷諾數(shù)和近激波面aspect ratio）在高超聲速條件下對數(shù)值仿真精度的影響，給出了最優(yōu)參數(shù)配置并對RBCC高超聲速飛行器在巡航段的氣動熱進(jìn)行了預(yù)示，結(jié)論如下。

（1）鈍后緣的仿真性能優(yōu)于尖后緣。

（2）不同于二維圓柱體，對NACA0012翼型，單元雷諾數(shù)的建議取值應(yīng)當(dāng)不小于16，近激波面網(wǎng)格的aspect ratio數(shù)值取決于選用的后緣形狀。對尖后緣翼型，aspect ratio建議取值為240。對鈍后緣翼型，aspect ratio的建議取值為320。

（3）在20米遠(yuǎn)場距離條件下采用鈍后緣+SST k-omega+單元雷諾數(shù)16+近激波面aspect ratio320+AUSM通量類型的密度基求解器可獲得最佳仿真效果，P/Pt和T/Tt的平均誤差率分別為3.06%和2.83%。RBCC高超聲速飛行器在飛行過程中的氣動熱可達(dá)到1885℃，熱障問題十分突出。