基于高精度格式的翼型失速特性脫體渦數(shù)值模擬

王 瑞，鐘伯文*

（1.南昌航空大學(xué) 飛行器工程學(xué)院，南昌 330063；2.江西省飛行器設(shè)計(jì)與氣動仿真重點(diǎn)實(shí)驗(yàn)室，南昌 330063）

機(jī)翼在大攻角狀態(tài)下的氣動特性因湍流分離流動的存在呈現(xiàn)高度的非定常特征。當(dāng)攻角超過某臨界值時，機(jī)翼的升力會急劇減小，阻力會大幅度增大，導(dǎo)致飛機(jī)處于失速狀態(tài)而引發(fā)飛行事故。因此，精確評估飛機(jī)的失速特性對飛機(jī)的設(shè)計(jì)和飛機(jī)的性能優(yōu)化具有重要意義。

通常，機(jī)翼翼型在小攻角狀態(tài)的流動接近定常狀態(tài)，可以使用傳統(tǒng)的雷諾平均方程（Reynolds-averaged Navier Stokes Equation，RANS）進(jìn)行較為精確的數(shù)值模擬。然而，當(dāng)攻角一旦增大至接近甚至超過某臨界值，翼型外流場會出現(xiàn)大幅度的分離，傳統(tǒng)的RANS 方法無法精確模擬這類湍流的大分離非定常流動。

除了RANS 方法，還有直接數(shù)值模擬（Direct Numerical Simualtion，DNS） 和大渦模擬（Large Eddy Simulation，LES）[1]等湍流模擬方法。從理論上講，DNS方法能夠準(zhǔn)確模擬這類分離流動，但由于其對網(wǎng)格量的要求很高，計(jì)算耗費(fèi)巨大，無法滿足飛機(jī)工程型號的應(yīng)用需求。LES 方法采用空間過濾的方式過濾湍流的小尺度脈動，直接模擬大尺度渦結(jié)構(gòu)，而小尺度湍流脈動的影響通過亞格子模型來模擬，LES 方法相對于DNS可以大幅度降低計(jì)算量，但在近壁面仍然對網(wǎng)格密度有很高的要求[2]，使得LES 難于獲得工程型號的實(shí)際應(yīng)用。波音公司的Spalart 等[3]根據(jù)LES 和RANS 兩種方法的特點(diǎn)和優(yōu)點(diǎn)，提出了RANS/LES 混合方法，即脫體渦模擬（Detached Eddy Simulation ，DES），在近壁區(qū)采用RANS 方法，在遠(yuǎn)離壁面區(qū)域采用LES 方法，使得其方法在飛機(jī)工程型號的應(yīng)用成為可能。其后，還形成了不同的改進(jìn)方法，如延遲脫體渦模擬方法[4]（Delayed-DES，DDES）、DDES 增強(qiáng)版本[5]（Improved-DDES，IDDES）及分區(qū)域脫體渦模擬方法[6]（Zonal-DES，ZDES）等。

目前DES 類方法已經(jīng)對翼型大分離流動問題進(jìn)行了廣泛研究和應(yīng)用[7-8]。Shur 等[9]將DES 方法應(yīng)用于NACA0012 翼型低速狀態(tài)的失速特性預(yù)測，并且與URANS 進(jìn)行了對比，得到了不錯的模擬效果。本文作者Zhong 等[10]將DES 方法與URANS、隱式大渦模擬（Implicit LES，ILES）進(jìn)行了對比，體現(xiàn)了DES 方法對預(yù)測NACA0012 翼型失速特性的較強(qiáng)能力。Li[11]將BL 模型（Baldwin -Lomax）與SA 模型（Spalart-Allmaras）結(jié)合構(gòu)造了基于SA+BL 模型的DES 方法，并對NACA 系列翼型進(jìn)行了失速特性的數(shù)值模擬。

分離流動的精確模擬不僅受到湍流模擬方法的影響，同時也與數(shù)值格式密切相關(guān)。Ghosal[12]和Chow 等[13]評估了數(shù)值誤差對LES 方法的影響，結(jié)果都發(fā)現(xiàn)數(shù)值格式對LES 計(jì)算結(jié)果的影響占比甚至超過了亞格子模型。Kawai 等[14]將基于B-L 模型（Baldwin-Lomax）的RANS/LES 混合方法與高階緊致差分格式（Compact Differencing Scheme，CMPT）結(jié)合模擬了NACA64A-006 薄翼型失速特性，模擬結(jié)果表明了相同網(wǎng)格密度下高階CMPT與傳統(tǒng)迎風(fēng)格式相比可以得到更可靠的計(jì)算結(jié)果。

綜上所述，翼型的失速現(xiàn)象所表現(xiàn)的分離流動十分復(fù)雜，精確高效的模擬對湍流模擬方法和數(shù)值格式具有很高的要求。盡管DES 方法已經(jīng)在翼型失速特性的預(yù)測中得到了應(yīng)用，但是在數(shù)值格式、方法的適用性及模擬精度等方面，仍然存在進(jìn)一步的提升空間，需要深入地探索研究。

本文采用脫體渦模擬（DES）方法，并將HLLC 低耗散近似黎曼求解器與五階WENO 格式結(jié)合對NACA0015 翼型的失速特性和相應(yīng)的流場進(jìn)行了數(shù)值模擬研究，并與二階Jameson 中心差分格式（Centered Differencing Scheme，CDS）[15]的數(shù)值模擬結(jié)果分別從失速特性、壓力系數(shù)分布及非定常湍流結(jié)構(gòu)3 個方面進(jìn)行對比分析。

1 數(shù)值模擬方法

1.1 HLLC 格式

Toro 等[16]為解決HLL 格式在間斷面會產(chǎn)生較大的數(shù)值耗散現(xiàn)象提出了HLLC 格式。HLLC 格式本身是一階精度并且構(gòu)造簡單，具有低耗散、高分辨率的特性，與Roe 格式相比不需要進(jìn)行熵的修正來保持物理數(shù)值解，適合跨音速和低速計(jì)算。

HLLC 通量和相關(guān)的黎曼求解器都包含了黎曼問題的所有波，不需要將方程線性化，而且不通過固定熵值的方法來解決低密度問題。圖1 給出了黎曼問題示意圖，Toro 等[16]改進(jìn)的兩態(tài)HLLC 格式可以準(zhǔn)確地解決孤立激波、間斷波和剪切波，并且在黎曼扇形內(nèi)部結(jié)構(gòu)中保留了初始的正密度和壓力。

圖1 黎曼問題示意圖

HLLC 格式求解的格式如下

1.2 WENO 格式

WENO 格式（Weighted Essentially Non-oscillatory）是Liu 等[17]在ENO 算法（Essentially Non-oscillatory）的基礎(chǔ)上發(fā)展的。近年來WENO 算法在復(fù)雜流動的計(jì)算得到了廣泛應(yīng)用[18]。Jiang 等[19]建議采用IS 光滑因子度量函數(shù)的光滑性，構(gòu)成了WENO-JS 格式。

五階WENO-JS 格式的三階數(shù)值通量在3 個模板中可以表示為

光滑因子表示為

非線性權(quán)重因子定義如下

式中：k 的取值范圍為[0，1，2]；ε 為防止分母為0 的充分小正數(shù)，本文取10-6；P 為用于控制ENO 的行為，本文取4。

最終的重構(gòu)函數(shù)形式如下

1.3 Jameson 中心差分格式

Jameson 中心差分格式用于空間離散，具體形式如下

該格式允許雜波振蕩，為了解決這個問題，使用耗散通量模擬了人工耗散，得到的最終Jameson 中心差分格式如下

1.4 時間推進(jìn)格式

以偏微分方程y′=f（x，y），y（x，y）=y0為例，在點(diǎn)x=x0+h 處，y 和h 為矢量。五階Runge-Kutta 格式的具體形式如下

式中：ν 為Runge-Kutta 階數(shù)，本文取5。

2 數(shù)值模擬分析

本文采用NACA0015 翼型作為研究算例，展長設(shè)置為0.2c，c 為翼型弦長。計(jì)算網(wǎng)格采用多塊結(jié)構(gòu)化網(wǎng)格，網(wǎng)格塊劃分如圖2 所示，流場計(jì)算域?yàn)?5c，網(wǎng)格點(diǎn)數(shù)分布為：流向點(diǎn)數(shù)647 個、法向點(diǎn)數(shù)115 個、展向點(diǎn)數(shù)81 個，網(wǎng)格單元總數(shù)為6×106。實(shí)驗(yàn)數(shù)據(jù)來源LM Glasfiber 風(fēng)洞[20]，計(jì)算的模擬流動狀態(tài)為：自由來流馬赫數(shù)為0.078，基于弦長的雷諾數(shù)為1.6×106。計(jì)算程序中，為了加速收斂和減少計(jì)算耗費(fèi)，非定常時間推進(jìn)采用雙時間步隱式迭代法。時間步長采用無量綱形式，時間步長Δt=0.01，對應(yīng)的物理時間步長為3.5×10-5s。

圖2 多塊網(wǎng)格劃分

2.1 失速特性預(yù)測分析

本文對比了HLLC+五階WENO 格式與CDS 的計(jì)算結(jié)果。圖3 和圖4 分別給出了不同數(shù)值格式得到的升力系數(shù)曲線（CL）和阻力系數(shù)（CD）曲線。圖3 結(jié)果表明，2 種數(shù)值格式都能準(zhǔn)確預(yù)測失速的臨界攻角及小攻角的升力系數(shù)。但是，在近失速段（攻角在8.13～13°范圍內(nèi)），CDS 預(yù)測得到的近失速升力系數(shù)與實(shí)驗(yàn)值相比偏小，13°攻角的升力系數(shù)預(yù)測偏差最大，而HLLC+五階WENO 格式得到的升力系數(shù)誤差較小，增長趨勢和實(shí)驗(yàn)值相同。在過失速段（攻角在14.3～19°范圍內(nèi)），HLLC+五階WENO 格式能夠準(zhǔn)確預(yù)測升力系數(shù)，而CDS 的預(yù)測的整體數(shù)值偏小。

圖3 不同數(shù)值格式的升力系數(shù)曲線

圖4 不同數(shù)值格式的阻力系數(shù)曲線

根據(jù)圖4 阻力系數(shù)預(yù)測結(jié)果，當(dāng)攻角小于17°，2種格式得到的數(shù)值結(jié)果基本一致，在臨界攻角附近與實(shí)驗(yàn)值相比偏大。當(dāng)攻角為17°時，CDS 預(yù)測得到的阻力系數(shù)與實(shí)驗(yàn)相比偏大，HLLC+五階WENO 格式的結(jié)果更接近實(shí)驗(yàn)值，而在19°攻角狀態(tài)下，2 種差分格式的預(yù)測結(jié)果都偏大。

2.2 壓力系數(shù)分布預(yù)測

圖5 給出了在15.66°攻角下時間平均后的翼型表面壓力系數(shù)（CP）分布模擬結(jié)果。結(jié)果表明，在翼型上表面，兩種格式得到的上表面壓力分布系數(shù)轉(zhuǎn)折點(diǎn)基本一致，并且與實(shí)驗(yàn)值相同，因此2 種格式都準(zhǔn)確地預(yù)測了流動的分離點(diǎn)位置。CDS 預(yù)測翼型表面壓力系數(shù)分布的整體數(shù)值偏小，這導(dǎo)致了該攻角下升力系數(shù)的數(shù)值結(jié)果小于實(shí)驗(yàn)值。HLLC+五階WENO 格式預(yù)測得到的壓力峰值和分離區(qū)的壓力系數(shù)偏大，而其他位置的數(shù)值與實(shí)驗(yàn)值比較接近，因此在該攻角下預(yù)測的升力系數(shù)比實(shí)驗(yàn)值略大。

圖5 AOA=15.66°時間平均的翼型表面壓力系數(shù)分布

表1 給出了15.66°攻角的平均升力系數(shù)對比情況。根據(jù)升力系數(shù)相對誤差的對比結(jié)果，HLLC+五階WENO 格式得到的升力系數(shù)預(yù)測的相對誤差較小，因此該格式得到的大攻角壓力系數(shù)分布的數(shù)值模擬結(jié)果更接近實(shí)驗(yàn)結(jié)果。

表1 AOA=15.66°平均升力系數(shù)對比

2.3 大攻角湍流結(jié)構(gòu)分析

圖6 給出了在19°攻角狀態(tài)下的瞬時湍流結(jié)構(gòu)等值面圖（Q 準(zhǔn)則）。結(jié)果顯示，2 種格式都能捕捉到大尺度湍流渦結(jié)構(gòu)。HLLC+五階階WENO 格式的耗散性相對較低、模擬精度較高，能夠捕捉到小尺度湍流渦結(jié)構(gòu)，能夠更為準(zhǔn)確地模擬大分離湍流運(yùn)動，從而更準(zhǔn)確地預(yù)測翼型的失速特性。而CDS 的模擬精度低、耗散性大，無法捕捉較小尺度的湍流渦結(jié)構(gòu)。

圖6 AOA=19°瞬時湍流結(jié)構(gòu)渦量等值面（Q 準(zhǔn)則）

3 結(jié)論

本文將HLLC 格式與五階WENO 格式結(jié)合，通過對比二階Jameson 中心差分格式的數(shù)值模擬結(jié)果，研究了高階數(shù)值格式對DES 方法模擬翼型失速特性的影響。

1）首先對比了HLLC+五階WENO 格式與CDS 的升力系數(shù)和阻力系數(shù)的預(yù)測結(jié)果。2 種格式都能較為準(zhǔn)確地預(yù)測失速臨界攻角，而HLLC+五階WENO 格式能提供與實(shí)驗(yàn)結(jié)果相比更為準(zhǔn)確的升力系數(shù)和阻力系數(shù)。

2）HLLC+五階WENO 格式的耗散性較小、模擬精度高，與CDS 相比能夠捕捉更多小尺度湍流渦結(jié)構(gòu)，HLLC+五階WENO 格式能夠提高DES 方法對翼型失速狀態(tài)下湍流分離渦結(jié)構(gòu)的模擬精度，從而能更為準(zhǔn)確地預(yù)測翼型的失速特性。