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

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

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

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

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

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

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

0 引言

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

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

1 Cart3D程序簡(jiǎn)介

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

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

三維可壓縮流動(dòng)的Euler方程可以表示為：

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

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

對(duì)于空氣，γ取1.4。

3 算例驗(yàn)證

3.1 ONERA M6機(jī)翼

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

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

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

3.2 某飛行器

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

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

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

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

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

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

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

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

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

4 結(jié)論

通過對(duì)ONERA M6機(jī)翼及某飛行器模型的計(jì)算驗(yàn)證可知：

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

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

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

[2]戚姝妮，張祖庚，董軍?；谧赃m應(yīng)笛卡爾網(wǎng)絡(luò)之三維標(biāo)模氣動(dòng)特性的數(shù)值模擬[J]。氣動(dòng)研究與試驗(yàn)，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）

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