多區(qū)域網(wǎng)格尾跡積分技術(shù)在復(fù)雜構(gòu)型氣動力計算中的應(yīng)用研究＊

葉 露，王永虎，王元元

（1中國民航飛行學(xué)院，四川廣漢 618307；2西北工業(yè)大學(xué)航空學(xué)院，西安 710072）

0 引言

傳統(tǒng)的氣動力計算主要分為近場法和遠(yuǎn)場法兩種。近場法的準(zhǔn)確度受表面曲率影響較大，特別對表面外形變化復(fù)雜的物體，可能帶來較大誤差。因而它需要數(shù)量很大的高質(zhì)量網(wǎng)格來精確的擬合飛行器實際外表面。而基于經(jīng)典動量定理的傳統(tǒng)遠(yuǎn)場法，由于遠(yuǎn)場網(wǎng)格稀疏及流場求解的耗散等原因，并不能滿足計算要求。

尾跡積分法是對傳統(tǒng)的遠(yuǎn)場積分法的一種改進(jìn)，可以避免近場法和傳統(tǒng)遠(yuǎn)場法的不足并且能將阻力按產(chǎn)生的物理機(jī)理分解，這樣有利于設(shè)計師進(jìn)行有目的的減阻設(shè)計。近年來，國內(nèi)外一些研究人員對尾跡積分法進(jìn)行了相關(guān)的研究[1－5]。研究結(jié)果證明了尾跡積分方法的優(yōu)勢，同時也得到了很多有益的結(jié)論，包括影響尾跡積分效率和精度的因素以及影響規(guī)律等。但從發(fā)表的文章看，大部分的算例都是針對簡單外形（翼型和單獨機(jī)翼）或者較復(fù)雜外形（單獨翼身組合體）進(jìn)行的，原因是這些外形可以生成單塊計算網(wǎng)格，網(wǎng)格有統(tǒng)一的I、J、K編號，在進(jìn)行尾跡插值區(qū)的搜索時很方便。而復(fù)雜外形的多區(qū)域網(wǎng)格在進(jìn)行插值區(qū)域的搜索時很不方便，工作量很大，且存在不同區(qū)域之間數(shù)據(jù)傳遞問題，因此很少有人研究。

文中采用“對角線判斷”方法提高了尾跡積分技術(shù)在多區(qū)域網(wǎng)格應(yīng)用中的搜索效率，以繞DLR－F6翼身、翼／身／架／艙組合體等復(fù)雜構(gòu)型的流場為算例，對尾跡面法氣動力積分技術(shù)進(jìn)行了探討。

1 主控方程和數(shù)值方法

主控方程采用Navie－Stokes方程，在慣性Cartesian坐標(biāo)下，三維可壓縮非定常Navier－Stokes方程積分形式的表達(dá)式為：

其中：Q是守恒變量；F為無粘通量；G為粘性項；?Ω為控制體的邊界；n為控制體邊界單位外法向矢量。時間推進(jìn)采用LU－SGS隱式推進(jìn)方法，空間離散格式采用二階Roe方法[6]，通量限制器采用Vanbada通量限制器。邊界條件物面采用無滑移邊界條件，遠(yuǎn)場采用一維法向黎曼不變量，實現(xiàn)遠(yuǎn)場無反射邊界條件。采用S－A湍流模型。

2 尾跡面法氣動力積分原理

文中采用的尾跡積分的計算公式為：

式中：ρ∞、U∞、P∞分別是無窮遠(yuǎn)來流的密度、x方向的速度、靜壓；S代表尾跡面；Δs、ψ、ξ分別是尾跡面上的熵增、流函數(shù)和旋度。具體的推導(dǎo)過程詳見參考文獻(xiàn)[3]和文獻(xiàn)[7]。

3 尾跡面流動參數(shù)插值

由于尾跡積分方法是基于氣流軸系推導(dǎo)而來，因此文中尾跡面的選取方法為：先確定尾跡面與X軸交點，過該點作垂直于來流的平面并視為“尾跡積分面”，利用幾何原理確定了尾跡面位置后，進(jìn)行尾跡面積分的關(guān)鍵點就是尾跡面上流動參數(shù)的插值，文中采用以下方法：

1）尾跡面與網(wǎng)格各區(qū)域的位置判斷

尾跡面與網(wǎng)格各區(qū)域的位置判斷采用對角線判斷方法，如圖1所示。

圖1 尾跡面與網(wǎng)格各區(qū)域的位置判斷

根據(jù)來流方向在確定了尾跡面方程后，再確定每個網(wǎng)格區(qū)域的體對角線及其端點，根據(jù)幾何原理判斷任何對角線與尾跡面是否有交點，若存在交點則相應(yīng)區(qū)域確定為尾跡插值搜索區(qū)，進(jìn)而在確定的搜索區(qū)進(jìn)行網(wǎng)格單元邊界與尾跡面幾何關(guān)系的判斷。這樣可以大大提高插值搜索區(qū)的搜索效率。

表1為逐點確定尾跡面網(wǎng)格的方法以及文中所采用的搜索方法在網(wǎng)格量為500萬中的應(yīng)用效率，可以看出，文中采用的搜索算法大大提高了計算效率，從而能夠快速的進(jìn)行尾跡積分。

2）基于格心距離加權(quán)的流動參數(shù)插值

與異側(cè)相鄰網(wǎng)格格心流動參數(shù)中心插值不同的是，基于格心距離加權(quán)的流動參數(shù)插值采用循環(huán)距離判斷方法尋找與O點最近的4個點A、B、C、D，且至少有兩個點在尾跡面異側(cè)，分別求出4點與插值點的距離dA、dB、dC、dD，權(quán)值系數(shù)以及在尾跡面插值由以下公式計算：

表1 確定尾跡面網(wǎng)格的搜索效率比較

其中q可以是原始變量也可以是守恒變量，如圖2所示。

圖2 距離加權(quán)插值

4 應(yīng)用算例

算例1：繞DLR－F6翼身組合體跨音速流動

計算狀態(tài)：自由來流Ma∞＝0.75，α＝1°，Re＝3.0×106

采用多塊點對接網(wǎng)格，共16塊，300萬單元。

圖3為翼身組合體表面網(wǎng)格以及壓力云圖分布，圖4、圖5為典型剖面壓力分布計算與實驗對比，結(jié)果顯示采用的流場求解器是可靠的。

圖3 表面網(wǎng)格和壓力云圖

圖4 15%半翼展處表面壓力系數(shù)分布

圖6 尾跡面渦量分布

圖6為不同位置尾跡面的渦量分布圖（圖中x／c代表尾跡面位置距機(jī)身尾部的距離與平均氣動弦長的比值）。圖7為尾跡面位置對尾跡積分結(jié)果的影響。尾跡面位置選取與網(wǎng)格分布有關(guān)，文中網(wǎng)格分布下可以看出尾跡面位置距機(jī)尾0.05倍 ～0.5倍的氣動弦長距離最為合理，此區(qū)間的阻力增長段由于數(shù)值耗散熵增形成，而下降段是渦量耗散引起。表2列出了F6翼身組合體氣動力計算結(jié)果。

圖7 氣動力積分與尾跡面位置的關(guān)系

表2 F6翼身組合體氣動力計算和實驗結(jié)果

算例2：繞DLR－F6機(jī)翼／機(jī)身／掛架／短艙組合體跨音速流動

計算狀態(tài)：自由來流

采用多塊點對接網(wǎng)格，共43塊，430萬單元。

圖8為F6翼／身／架／艙組合體表面網(wǎng)格以及壓力云圖分布，圖9、圖10為典型剖面壓力分布，計算與實驗數(shù)據(jù)吻合良好。

表3給出了F6翼身架艙組合體氣動力計算結(jié)果。從結(jié)果可以看出，隨尾跡面位置向下游移動，熵增阻力逐漸增大，升力系數(shù)逐漸減小，這是因為越遠(yuǎn)離機(jī)翼，由于方程數(shù)值耗散引起的“熵增”也越大，從而覆蓋了真實的熵增阻力，同時渦量耗散使得升力系數(shù)減小。同樣的，可以看出在文中網(wǎng)格分布下尾跡面位置距機(jī)尾0.05倍～0.5倍的氣動弦長距離最為合理。

圖8 表面網(wǎng)格分布以及壓力云圖

圖9 15%半翼展處表面壓力系數(shù)分布

圖10 33.1%半翼展處表面壓力系數(shù)分布

表3 F6翼身架艙組合體氣動力計算和實驗結(jié)果

圖11給出了各種方法計算出的極曲線的對比。從圖中看出，采用文中的多區(qū)域網(wǎng)格尾跡積分技術(shù)得到的極曲線與表面積分結(jié)果和實驗數(shù)據(jù)比較吻合，尤其是距離機(jī)身尾部0.1倍平均氣動弦長的尾跡面積分精度最高。

圖11 各種方法計算出的極曲線對比

5 結(jié)論

1）文中很好的運(yùn)用了多區(qū)域網(wǎng)格尾跡積分方法求解了DLR－F6復(fù)雜構(gòu)型的氣動力，數(shù)值試驗表明，尾跡面的位置處于x／c∈（0.05，0.5）范圍內(nèi)時，具備較好的氣動力積分精度。

2）采用“對角線判斷”方法進(jìn)行多區(qū)域網(wǎng)格尾跡面插值搜索區(qū)的搜索，可以大大減少搜索工作量，提高尾跡積分的計算效率。

3）采用尾跡面積分法計算復(fù)雜外形的氣動力可以在保證和物面積分相當(dāng)精度的前提下將阻力按產(chǎn)生機(jī)理分解，對于設(shè)計師有針對性的減阻設(shè)計是很有利的。

4）特別地，對于復(fù)雜構(gòu)型來說尾跡面位置選擇經(jīng)驗性較強(qiáng)，從文中研究結(jié)果看尾跡面要盡量靠近機(jī)身，從而降低機(jī)身長度對機(jī)翼產(chǎn)生的流場變量的影響，諸如渦量、數(shù)值耗散等。

[1]Cummings R M，Giles M B，Shrinivas G N.Analysis of the elements of drag in three－dimensional viscous and inviscid flows，AIAA－96－2482[R].1996.

[2]Hunt D L，Cummings R M，Giles M B.Determination of drag from three－dimensional viscous and inviscid flowfield computation，AIAA－97－2257[R].1997.

[3]Zhu Ziqiang，Wang Xialu，Liu Jie，et al.Comparison of predicting drag methods using computational fluid dynamics in 2d／3dviscous flow[J].Science in China Series E：Technological Science，2007，50（5）：534－549.

[4]Van Dam C P，Nikfetrat K.Drag prediction at subsonic and transonic speeds using Euler methods[J].Journal of Aircraft，1995，32（4）：839－845.

[5]Paparone L，Tognaccini R.Computational fluid dynamics－based drag prediction and decomposition[J].AIAA Journal，2003，41（9）：1647－1657.

[6]Roe P L.Approximate Riemann solvers，parameter vector and difference schemes[J].Journal of Computational Physics 1997，135（2）：250－258.

[7]J Mark Janus，Animesh Chatterjee.Use of a wake－integral method for computational drag analysis[J].AIAA Journal，1996，34（1）：188－190.