羅劍波，栗莉，方明恩，羅帥，綦龍，張巖

（航空工業(yè)洪都，江西南昌330024）

笛卡爾網(wǎng)格在氣動設(shè)計中的應(yīng)用研究

羅劍波，栗莉，方明恩，羅帥，綦龍，張巖

（航空工業(yè)洪都，江西南昌330024）

采用非貼體笛卡爾網(wǎng)格及有限體積法求解Euler方程的方法對兩個模型的氣動特性進行了計算，并考慮粘性項對阻力結(jié)果進行了修正，將修正后的計算結(jié)果與風洞試驗值進行對比分析，結(jié)果表明：修正后的CFD計算值與風洞試驗值吻合較好，驗證了在氣動設(shè)計中笛卡爾網(wǎng)格方法具有計算精度高、求解速度快的特點，適用于氣動方案初步設(shè)計階段的實際工程應(yīng)用。

笛卡爾網(wǎng)格；Euler方程；數(shù)值模擬；粘性阻力修正

0 引言

在飛行器氣動設(shè)計過程中，尤其在氣動方案初步設(shè)計階段，往往需要進行大量CFD仿真，以獲得各狀態(tài)下飛行器的氣動特性，用于氣動方案的選型。以往采用的方法是：首先對設(shè)計的氣動外形生成結(jié)構(gòu)/非結(jié)構(gòu)網(wǎng)格，再求解N-S方程或Euler方程，最后對求解結(jié)果進行處理分析。采用這種方法通常能獲得較為精確的解，但網(wǎng)格劃分和流場求解過程需花費大量的人力和計算時間，限制了其在氣動方案設(shè)計階段的應(yīng)用。

近年來，基于笛卡爾網(wǎng)格技術(shù)求解Euler方程的方法越發(fā)成熟，與傳統(tǒng)CFD仿真方法不同的是：該方法首先在全流場域生成各向尺寸一致的粗糙網(wǎng)格，再根據(jù)模型結(jié)構(gòu)在物面附近自動逐步加密，得到尺寸合適的流場網(wǎng)格，最后求解Euler方程得到流場結(jié)果，該方法網(wǎng)格生成效率高、流場求解速度快，能大大縮短方案階段的氣動外形設(shè)計時間。

1 Cart3D程序簡介

本文采用的笛卡爾網(wǎng)格劃分及求解Euler方程的過程在Cart3D程序中完成。Cart3D程序能通過定義網(wǎng)格區(qū)域及網(wǎng)格密度，自動捕捉模型的幾何特征，快速生成非貼體笛卡爾網(wǎng)格（圖1），極大地壓縮網(wǎng)格生成時間；求解模塊基于有限體積法求解Euler方程，采用顯式龍格庫塔法時間推進，空間離散采用迎風格式，并采用通量限制器提高精度，具有TVD高階格式特性，在求解過程中，可采用多重網(wǎng)格、自適應(yīng)網(wǎng)格等措施提高計算效率和計算精度，適用于亞、跨、超音速復雜飛行器的外流場分析。

2 控制方程與數(shù)值算法

三維可壓縮流動的Euler方程可以表示為：

式中：U為守恒變量；F、G、H為無粘通量，分別表示為：

式中：ρ、p、e分別表示密度、壓力和單位質(zhì)量總能；u、v、w分別為x、y、z三個方向的速度分量。假設(shè)來流為完全氣體，則方程組可以使用下式封閉：

對于空氣，γ取1.4。

3 算例驗證

3.1 ONERA M6機翼

ONERA M6機翼因其外形簡單，在跨音速環(huán)境中機翼表面繞流呈現(xiàn)出局部超音速流動、激波、邊界層分離等復雜的流動狀態(tài)，因而成為驗證CFD方法的典型算例。對ONERA M6機翼劃分的笛卡爾網(wǎng)格，如圖2所示，采用Euler方程求解，并將計算的壓力分布與NASA標準計算結(jié)果及ONERA S2MA的風洞試驗值進行對比分析，驗證采用笛卡爾網(wǎng)格離散并求解Euler方程方法的計算精度。計算條件為：馬赫數(shù)Ma=0.84，攻角α=3.06°，側(cè)滑角β=0°，Re=11.72×106（基于平均氣動弦長bA=0.64607m），參考面積Sref=0.7532m2。

圖3左側(cè)是采用笛卡爾網(wǎng)格求解Euler方程得到的ONERA M6機翼上表面壓力云圖，右側(cè)是NASA求解N-S方程得到的標準計算結(jié)果，可以看出兩種計算方法得到的上翼面壓力分布具有相當高的一致性，壓力分布結(jié)果基本吻合，采用笛卡爾網(wǎng)格的方法很好地捕捉到了彈翼上表面的λ型激波和波后的流動。

提取ONERA M6機翼展向z/B=0.20、0.44、0.65、0.80、0.90、0.95處6個典型截面位置的壓力分布，并與ONERA S2MA風洞的試驗值進行對比。如圖4所示，圖例中Exp表示風洞試驗值，Euler為采用笛卡爾網(wǎng)格方法的計算值。結(jié)果表明：彈翼下表面流動相對簡單，笛卡爾網(wǎng)格方法精確地模擬了彈翼下表面的流動，各截面處計算的壓力值與試驗的壓力值吻合較好；而彈翼上表面由于出現(xiàn)激波、邊界層分離等復雜流動現(xiàn)象，笛卡爾網(wǎng)格方法的計算結(jié)果與試驗值之間有一定的偏差，彈翼上表面前緣處的壓力峰值的位置及壓力峰值的模擬較為準確，但對彈翼上表面中部激波位置及激波強度的模擬不大準確；總體上看計算結(jié)果與試驗值吻合較好，可驗證笛卡爾網(wǎng)格方法的計算精度。

3.2 某飛行器

采用笛卡爾網(wǎng)格方法對某飛行器算例的升阻特性隨攻角和馬赫數(shù)的變化趨勢進行大量數(shù)值計算，并將計算與風洞試驗結(jié)果進行比較。飛行器頭部附近網(wǎng)格如圖5所示。

利用單機計算了某飛行器在亞音速、超音速下共24個飛行狀態(tài)的升阻特性，完成整個網(wǎng)格劃分及流場求解過程總計耗時不超過3H。由于笛卡爾網(wǎng)格方法僅求解Euler方程，不能預測物面的摩擦阻力，故在計算完成后對阻力系數(shù)計算結(jié)果進行了粘性修正。

本文對阻力系數(shù)采用的粘性修正方法為：

1）利用Fluent等求解N-S方程，得到飛行器的零升阻力系數(shù)Cd0NS；

2）利用笛卡爾網(wǎng)格方法求解Euler方程，得到零升阻力系數(shù)Cd0E；

3）計算零升阻力系數(shù)差量ΔCd0=Cd0NS-Cd0E；

4）將零升阻力系數(shù)差量ΔCd0疊加到笛卡爾網(wǎng)格方法的阻力系數(shù)求解結(jié)果，得到修正后的阻力系數(shù)。

圖6～圖8為飛行器升阻力系數(shù)的計算值與試驗值的對比，其中：圖例Exp表示風洞試驗值，Euler為笛卡爾網(wǎng)格方法的CFD計算值，CFD阻力修正為對笛卡爾網(wǎng)格方法的CFD阻力計算值進行粘性修正的結(jié)果?？梢钥闯觯荷ο禂?shù)的CFD計算值與試驗值基本重合；而未修正的阻力系數(shù)計算值與試驗值存在差量，但粘性修正后的阻力系數(shù)曲線與風洞試驗值基本重合，且對飛行器的阻力發(fā)散馬赫數(shù)的預測較為準確，說明阻力修正方法可行。由此可見，笛卡爾網(wǎng)格方法能較為準確的預測飛行器的升力系數(shù)，阻力系數(shù)也能通過粘性修正獲得較為精度的結(jié)果。

與傳統(tǒng)采用貼體結(jié)構(gòu)/非結(jié)構(gòu)網(wǎng)格進行空間離散，并采用Fluent、CFX等求解器求解N-S方程的方法比較，采用笛卡爾網(wǎng)格求解Euler方程的方法不僅具有較高的計算精度，而且因其網(wǎng)格生成迅速、求解效率高，能大幅度縮短設(shè)計周期。就上述飛行器模型而言，利用單機對24個飛行狀態(tài)從建模、求解到后處理，采用笛卡爾網(wǎng)格求解Euler方程的方法共耗時不超過3H，而相同條件下若利用傳統(tǒng)方法處理耗時則不少于4天。

4 結(jié)論

通過對ONERA M6機翼及某飛行器模型的計算驗證可知：

1)笛卡爾網(wǎng)格離散方法能較為準確地模擬物面的流動，包括激波位置、激波強度及物面流動分離現(xiàn)象，具有較高的精度，可滿足方案階段氣動特性快速設(shè)計的工程需求；

2)利用笛卡爾網(wǎng)格離散方法求解Euler方程，配合多重算法的加速方法，網(wǎng)格生成迅速、求解速度快，計算結(jié)果可靠，即使在單機情況下，也比常規(guī)采用N-S方程求解器的效率提高至少10倍，從而大大縮短設(shè)計周期。

[1]肖涵山，陳作斌，劉剛，江雄。基于Euler方程的三維自適應(yīng)笛卡爾網(wǎng)格應(yīng)用研究[J]。空氣動力學學報，2003；21（2）：202-210.

[2]戚姝妮，張祖庚，董軍。基于自適應(yīng)笛卡爾網(wǎng)絡(luò)之三維標模氣動特性的數(shù)值模擬[J]。氣動研究與試驗，2005；22（3）：1-8.

Application Research of Cartesian Grid in Aerodynamic Design

Luo Jianbo,Li Li,Fang Mingen,Luo Shuai,Qi Long,Zhang Yan
(AVIC-HONGDU,Nanchang,Jiangxi 330024)

This paper adopts the way of non-body-fitted Cartesian grid and finite volume method for solving Euler equation to calculate the aerodynamic characteristics,and corrects drag result with giving consideration to viscosity term.The corrected calculation result is analyzed by comparing with wind tunnel test value.The result shows:the corrected CFD calculation value fits wind tunnel test value well,so it demonstrates that Cartesian grid method features highcalculation accuracy and fast solving speed,and is applicable to the practical engineering application of aerodynamic scheme in initial design stage.

Cartesian grid；Euler equation；numerical simulation；correction of viscous drag

2017-08-29）

>>>作者簡介 羅劍波，1989年5月出生，2011年畢業(yè)于南京航空航天大學，工程師，現(xiàn)主要從事飛行器氣動設(shè)計工作。