笛卡爾網格在氣動設計中的應用研究

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

（航空工業洪都，江西南昌330024）

笛卡爾網格在氣動設計中的應用研究

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

（航空工業洪都，江西南昌330024）

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

笛卡爾網格；Euler方程；數值模擬；粘性阻力修正

0 引言

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

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

1 Cart3D程序簡介

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

2 控制方程與數值算法

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

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

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

對于空氣，γ取1.4。

3 算例驗證

3.1 ONERA M6機翼

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

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

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

3.2 某飛行器

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

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

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

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

2）利用笛卡爾網格方法求解Euler方程，得到零升阻力系數Cd0E；

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

4）將零升阻力系數差量ΔCd0疊加到笛卡爾網格方法的阻力系數求解結果，得到修正后的阻力系數。

圖6～圖8為飛行器升阻力系數的計算值與試驗值的對比，其中：圖例Exp表示風洞試驗值，Euler為笛卡爾網格方法的CFD計算值，CFD阻力修正為對笛卡爾網格方法的CFD阻力計算值進行粘性修正的結果。可以看出：升力系數的CFD計算值與試驗值基本重合；而未修正的阻力系數計算值與試驗值存在差量，但粘性修正后的阻力系數曲線與風洞試驗值基本重合，且對飛行器的阻力發散馬赫數的預測較為準確，說明阻力修正方法可行。由此可見，笛卡爾網格方法能較為準確的預測飛行器的升力系數，阻力系數也能通過粘性修正獲得較為精度的結果。

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

4 結論

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

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

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

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

[2]戚姝妮，張祖庚，董軍。基于自適應笛卡爾網絡之三維標模氣動特性的數值模擬[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年畢業于南京航空航天大學，工程師，現主要從事飛行器氣動設計工作。