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

汪 鈺,李 芃

（銅陵學院，安徽 銅陵 244000）

0 引言

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

1 仿真基礎

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

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

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

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

1.2 計算域外流場

如圖2所示，基于NACA4 airfoil generator和air‐foil tools設計了尖后緣和鈍后緣兩種NACA0012翼型并繪制了計算域外流場，其中翼型的特征長度L=1m。對外流場的數值模擬，一般建議遠場距離為翼型特征長度的12-20倍。建立遠場距離為20m的C型網格拓撲以最小化近壁網格的偏度。尖后緣采用ICEM CFD進行網格劃分時可直接折疊后緣處的塊，如圖2（a）所示。鈍后緣保留后緣處的塊，如圖2（b）所示。

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

1.3 湍流模型

對于涉及可壓縮和熱傳遞的流動，需要求解質量、動量和能量守恒方程。此外，還要求解湍流帶來的輸運方程。考慮到湍流方程無法通過直接數值模擬（DNS）求解，采用了目前應用最為廣泛的雷諾平均可壓縮NS方程（RANS）。K-epsilon、SST komega和Spalart-Allmaras是被廣泛應用且相對準確的RANS湍流模型［9-11］。K-epsilon模型主要用于工業CFD仿真，該模型存在對邊界層分離和逆壓梯度不敏感的問題，所以在翼型外流場的應用并不廣泛。Spalart-Allmaras是為航空應用設計的單方程模型，能夠為受逆壓梯度影響的邊界層仿真提供良好的結果。SST k-omega是基于標準k-omega模型設計的雙方程模型，其對自由來流的初始條件不敏感，且經過校準，能夠準確計算光滑表面的流動分離。因此采用雙方程的SST k-omega湍流模型。

2 驗證測試

2.1 網格策略

因為巡航段的馬赫數M大于5，此時的計算域外流場處于高超聲速階段，所以采用參考文獻［12］中的數據來驗證仿真參數。計算域外流場如圖3（a）所示，其中粗黑點為坐標原點，a、b和c是輸入邊界INPUT boundary，e、f和g是壁面邊界WALL boundary，d是輸出邊界OUTPUT boundary。湍流流動受翼型壁面的影響，因此，近壁面區域的網格劃分十分重要，需要保證邊界層的第一層網格處于sub viscosity layer，關鍵在于計算第一層網格高度y H。相關計算公式如式（1）：

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

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

圖3 （a）計算域外流場；（b）近壁面網格加密區

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

表1 網格無關性

2.2 數值方法

表2給出了在高超聲速階段采用的數值仿真參數。求解器選用專為高速可壓縮流設計的密度基，通量采用ROE和AUSM類型，其仿真結果與試驗數據較為接近。SST k-omega湍流模型應用intensity and viscosity ratio turbulence method，其中的湍流強度（%）和湍流黏度比均為1。計算域的邊界條件如下：輸入邊界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。梯度計算的常用方法有cell-based和node-based，第二種方法的準確度更高，尤其是針對具有傾斜度的非結構化網格，但是其對于硬件資源的要求更高。考慮到生成的結構化網格和有限計算資源，采用least-squares cell-based方案計算梯度，該方法可以提供與node-based相當的精度且計算開銷較低。對于壓力，采用基于中心差分的second-order方案。

表2 數值仿真參數

2.3 數值結果分析

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

對尖后緣，上述參數配置下的P/Pt分布如圖4（a）所示，其中X軸為x/L，Y軸為P/Pt的數值。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配置下的總平均仿真誤差率最低。接著基于該配置，對近激波面處網格的aspect ratio分別取值30、60、120、180、240和300再執行6組仿真，如圖4（d）所示，對應的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%）。當as‐pect ratio取值較小時，壓力和溫度的平均誤差率較大。之后隨著aspect ratio的增大，平均誤差率開始下降。當aspect ratio等于240時可獲得最佳結果，總的平均誤差率約3.93%。

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

對鈍后緣，上述參數配置下的P/Pt分布如圖5（a）所示，其中X軸為x/L，Y軸為P/Pt的數值。類似的，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配置下的平均誤差率最低。接著基于該配置，對近激波面處網格的aspect ratio分別取值40、80、160、240、320和400再執行6組仿真，如圖5（d）所示，對應的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%）。當aspect ratio取值較小時，壓力和溫度的平均誤差率較大。之后隨著aspect ratio的增大，平均誤差率開始下降。當as‐pect ratio等于320時可獲得最佳結果，總的平均誤差率約2.95%。

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

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

3 氣動熱預示

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

表3 數值仿真參數

4 結論

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

（1）鈍后緣的仿真性能優于尖后緣。

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

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