李 立，白 文，梁益華

(1.中國(guó)航空工業(yè)集團(tuán)公司第六三一研究所，陜西 西安 710068;2.中國(guó)航空研究院航空數(shù)值模擬技術(shù)研究應(yīng)用中心，北京 100012)

0 引言

隨著技術(shù)進(jìn)步，計(jì)算流體力學(xué)(CFD)已成為飛行器氣動(dòng)力設(shè)計(jì)的例行工具和日常應(yīng)用手段。但計(jì)算網(wǎng)格生成一直是制約CFD技術(shù)發(fā)展的一個(gè)瓶頸。通常對(duì)于復(fù)雜的計(jì)算問(wèn)題，整個(gè)數(shù)據(jù)準(zhǔn)備的70%以上工作量都被用于生成計(jì)算網(wǎng)格。近年來(lái)，非結(jié)構(gòu)網(wǎng)格(尤其是非結(jié)構(gòu)/結(jié)構(gòu)混合網(wǎng)格)因其對(duì)復(fù)雜幾何外形的適應(yīng)能力和良好的自動(dòng)生成能力，在復(fù)雜外形流場(chǎng)模擬中的應(yīng)用受到越來(lái)越多的重視［1］。本文圍繞非結(jié)構(gòu)自適應(yīng)網(wǎng)格生成技術(shù)開(kāi)展研究，其目標(biāo)是改進(jìn)非結(jié)構(gòu)網(wǎng)格技術(shù)在復(fù)雜流動(dòng)物理流場(chǎng)的模擬能力。

網(wǎng)格自適應(yīng)技術(shù)的核心是構(gòu)造各種自適應(yīng)探測(cè)器。不同的探測(cè)器對(duì)應(yīng)不同的網(wǎng)格自適應(yīng)方法［2］。目前，主要有兩類構(gòu)造自適應(yīng)探測(cè)器的典型方法。一類是基于流動(dòng)特征的傳統(tǒng)方法，如采用流場(chǎng)變量的差量、梯度等作為探測(cè)器。該類方法普遍認(rèn)為［3］，流場(chǎng)中梯度變化大的位置對(duì)流場(chǎng)計(jì)算精度的影響越大。主流商業(yè)軟件CFX、Fluent等所采用的策略都屬此類。一類是基于敏感性分析的新興方法，如本文研究的基于流場(chǎng)伴隨方程的網(wǎng)格自適應(yīng)探測(cè)器。與傳統(tǒng)方法相比，基于伴隨方程的網(wǎng)格自適應(yīng)方法是一種基于全局后驗(yàn)誤差估計(jì)的方法，直接建立了流場(chǎng)誤差與目標(biāo)輸出(如升力、阻力和力矩等)誤差的關(guān)聯(lián)關(guān)系，可有效避免因錯(cuò)誤判斷加密區(qū)域?qū)е碌木W(wǎng)格過(guò)度加密［3］，提高計(jì)算效率和目標(biāo)輸出的計(jì)算精度［4-5］，從而有效削減氣動(dòng)數(shù)值模擬對(duì)計(jì)算網(wǎng)格的敏感性。

本文對(duì)基于伴隨方程的網(wǎng)格自適應(yīng)技術(shù)的基本原理進(jìn)行了研究，構(gòu)造了基于流場(chǎng)伴隨方程的網(wǎng)格自適應(yīng)探測(cè)器(伴隨自適應(yīng)探測(cè)器)，利用局部加密網(wǎng)格的辦法，建立并實(shí)現(xiàn)了基于Euler方程的無(wú)粘伴隨自適應(yīng)能力。對(duì)于粘性流動(dòng)問(wèn)題，提出首先應(yīng)用無(wú)粘計(jì)算和伴隨自適應(yīng)從基準(zhǔn)網(wǎng)格自動(dòng)生成具有較好網(wǎng)格分辨率的非結(jié)構(gòu)自適應(yīng)網(wǎng)格，然后采用層推進(jìn)法自動(dòng)生成粘性網(wǎng)格的求解策略。通過(guò)對(duì)具體算例開(kāi)展的數(shù)值試驗(yàn)，表明本文建立的方法是有效的。

1 數(shù)值方法

1.1 流場(chǎng)主控制方程的計(jì)算方法

出發(fā)方程為Euler方程或雷諾平均的N-S方程，應(yīng)用格點(diǎn)格式的有限體積方法離散，控制體采用對(duì)偶方法生成，空間離散采用標(biāo)準(zhǔn)的Jameson中心差分格式，時(shí)間推進(jìn)采用三階顯式Runge-Kutta格式。另外，在實(shí)際計(jì)算中，均采用三重多重網(wǎng)格W循環(huán)的聚合多重網(wǎng)格算法進(jìn)行計(jì)算加速。對(duì)于粘性計(jì)算，湍流模型采用Menter SST k-ω 模型。

1.2 流場(chǎng)伴隨方程的計(jì)算方法

流場(chǎng)伴隨方程的定義式如下［4-5］:

其中，λ即伴隨變量，R為流場(chǎng)控制方程殘差，Q為流場(chǎng)解變量，f為目標(biāo)輸出函數(shù)(如升力、阻力、力矩等)。顯然，當(dāng)流場(chǎng)解已知，伴隨變量就可以通過(guò)求解式(1)得到。為了解的穩(wěn)定性和計(jì)算過(guò)程的魯棒性，本文采用求解流場(chǎng)主控制方程的類似辦法，通過(guò)在伴隨方程中添加虛擬時(shí)間項(xiàng)進(jìn)行時(shí)間推進(jìn)求解［7］:

這里，Rλ(λ)為伴隨方程的殘差，

實(shí)際求解時(shí)，對(duì)式(2)，與流場(chǎng)主控制方程一樣，采用格點(diǎn)格式的有限體積方法離散，空間離散采用標(biāo)準(zhǔn)Jameson中心差分格式，時(shí)間推進(jìn)采用三階Runge-Kutta顯式格式，并采用多重網(wǎng)格加速。

1.3 伴隨自適應(yīng)探測(cè)器

通過(guò)伴隨變量可直接建立目標(biāo)輸出殘差與流場(chǎng)殘差之間的關(guān)系［4］:

式中，(λ0，Q0)表示初始網(wǎng)格上的伴隨解和流場(chǎng)解;(λ*，Q*)表示理論解或非常密網(wǎng)格上的解，一般無(wú)法直接給出;R(Q0)表示粗網(wǎng)格上流場(chǎng)控制方程的殘差;表示在非常密網(wǎng)格上伴隨方程的殘差。

由式(4)不難分析，影響目標(biāo)輸出的誤差項(xiàng)可分為兩部分。一部分是可利用現(xiàn)有網(wǎng)格直接計(jì)算的絕對(duì)誤差項(xiàng)(λ0)TR(Q0)，而剩下的一部分是需在更密網(wǎng)格上計(jì)算才能得到的修正誤差項(xiàng)。這表明，如對(duì)修正誤差項(xiàng)加以估計(jì)，就可判斷計(jì)算網(wǎng)格中哪些位置對(duì)目標(biāo)輸出有更大的影響，進(jìn)而可用于建立新型的網(wǎng)格自適應(yīng)探測(cè)器，即伴隨自適應(yīng)探測(cè)器。本文中，伴隨自適應(yīng)探測(cè)器定義如下:

通過(guò)估計(jì)或判斷該值大小就可標(biāo)記出流場(chǎng)中哪些位置的網(wǎng)格最需要加密。

1.4 自適應(yīng)網(wǎng)格生成策略

自適應(yīng)網(wǎng)格生成的一般性策略和方法在文獻(xiàn)［2，8］有較詳細(xì)的介紹。本文采用局部加密網(wǎng)格的辦法(又稱h-方法)來(lái)進(jìn)行自適應(yīng)網(wǎng)格生成。h-方法通過(guò)細(xì)劃網(wǎng)格單元來(lái)實(shí)現(xiàn)自適應(yīng)網(wǎng)格生成，通常分為點(diǎn)插入、粗網(wǎng)格細(xì)分和網(wǎng)格重構(gòu)等三類方法。本文采用的方法屬于粗網(wǎng)格細(xì)分法。該方法依賴于自適應(yīng)探測(cè)器所標(biāo)記的邊。根據(jù)邊的標(biāo)記情況，對(duì)于四面體，共有三種細(xì)分網(wǎng)格的模式，分別為二分模式、四分模式和八分模式。本文自適應(yīng)網(wǎng)格生成策略如下:首先通過(guò)伴隨自適應(yīng)探測(cè)器對(duì)網(wǎng)格單元進(jìn)行標(biāo)記，然后根據(jù)網(wǎng)格單元的標(biāo)記情況，采用二分模式進(jìn)行各向同性加密。對(duì)四面體單元，其方法是連接四面體最長(zhǎng)邊的中點(diǎn)和該邊所對(duì)的2個(gè)點(diǎn)，從而將四面體分成2個(gè)子四面體。為了利用式(5)定義的自適應(yīng)探測(cè)器對(duì)網(wǎng)格單元或邊進(jìn)行標(biāo)記，公式中的理論解或非常密網(wǎng)格上的解首先必須用從原始粗網(wǎng)格的解(λ0，Q0)往一個(gè)全局加密網(wǎng)格上插值得到的解(λh，Qh)代替。這里，全局加密網(wǎng)格采用八分模式進(jìn)行網(wǎng)格劃分來(lái)自動(dòng)生成。為了簡(jiǎn)單，文中采用分段線性插值來(lái)計(jì)算全局加密網(wǎng)格的解。通過(guò)循環(huán)掃描粗網(wǎng)格的每個(gè)單元所包含的全局加密網(wǎng)格的網(wǎng)格點(diǎn)并由式(5)來(lái)估計(jì)當(dāng)?shù)貙?shí)際的修正誤差:

自適應(yīng)準(zhǔn)則定義如下:用戶給定全局誤差閥值ε，按照誤差等分布的原則，判斷當(dāng)?shù)卣`差與平均全局誤差的比值η=N·εk/ε是否大于1，如該條件成立，則對(duì)所在的單元進(jìn)行標(biāo)記。這里，N表示現(xiàn)有粗網(wǎng)格的網(wǎng)格點(diǎn)數(shù)。

值得指出，采用二分模式對(duì)四面體網(wǎng)格進(jìn)行自適應(yīng)剖分可能會(huì)導(dǎo)致自適應(yīng)后的網(wǎng)格包含懸掛點(diǎn)或非四面體單元，并且還有可能使自適應(yīng)網(wǎng)格的網(wǎng)格質(zhì)量比較差，因此采取措施消除懸掛點(diǎn)和對(duì)網(wǎng)格進(jìn)行光順?lè)浅１匾?］。推薦將自適應(yīng)后的網(wǎng)格導(dǎo)入通用網(wǎng)格生成軟件進(jìn)行必要的處理。

1.5 粘性網(wǎng)格生成

為了捕捉邊界層內(nèi)的粘性流動(dòng)，必須在沿物面法向的一定厚度內(nèi)生成各向異性網(wǎng)格，即需將原始的純四面體非結(jié)構(gòu)網(wǎng)格轉(zhuǎn)換成混合網(wǎng)格?；旌暇W(wǎng)格中可能包含的網(wǎng)格單元類型有四面體、棱柱、六面體或金字塔單元。有多種方法可以實(shí)現(xiàn)從四面體網(wǎng)格自動(dòng)生成混合網(wǎng)格，本文推薦采用層推進(jìn)方法［10］。其思路是在物面附近按期望的邊界層厚度通過(guò)網(wǎng)格局部重構(gòu)產(chǎn)生貼體的棱柱單元，而在這一厚度之外的大部分網(wǎng)格單元保持不變，過(guò)渡單元采用四面體或金字塔單元填充。

2 程序校驗(yàn)和方法驗(yàn)證

2.1 RAE2822跨聲速翼型計(jì)算

RAE2822跨聲速翼型是一個(gè)超臨界翼型，在全世界得到廣泛研究，其試驗(yàn)結(jié)果被用于大量流場(chǎng)解算器的驗(yàn)證［11］。這里，選取的計(jì)算狀態(tài)為 Case 9:M∞=0.734，攻角α=2.79°。初始計(jì)算網(wǎng)格如圖1(a)所示，共有15852個(gè)網(wǎng)格點(diǎn)，31388個(gè)三角形網(wǎng)格單元。圖1(b)給出由初始網(wǎng)格得到的馬赫數(shù)云圖?？梢钥吹?，該初始網(wǎng)格較好地捕捉到流場(chǎng)的激波結(jié)構(gòu)。選用升力系數(shù)作為目標(biāo)輸出函數(shù)進(jìn)行自適應(yīng)，設(shè)定全局誤差閥值為0.1。自適應(yīng)后的網(wǎng)格如圖1(c)所示，共有20763個(gè)網(wǎng)格點(diǎn)，41141個(gè)三角形單元?？梢钥吹?，伴隨自適應(yīng)加密的區(qū)域主要在翼型的前緣、上表面和激波區(qū)，這和主要流動(dòng)的特征非常匹配，也正是我們希望加密的區(qū)域。按伴隨自適應(yīng)探測(cè)器的基本原理，這同時(shí)還表明，RAE2822翼型的前緣、上表面和激波區(qū)對(duì)于準(zhǔn)確計(jì)算升力有主要貢獻(xiàn)。

2.2 ONERA-M6跨聲速機(jī)翼計(jì)算

圖1 RAE2822跨聲速翼型的伴隨自適應(yīng)結(jié)果Fig.1 Effect of mesh adaptation for RAE2822 airfoil

圖2 ONERA-M6跨聲速機(jī)翼的伴隨自適應(yīng)結(jié)果Fig.2 Effect of mesh adaptation for ONERA-M6 wing

圖3 自適應(yīng)前后ONERA-M6機(jī)翼物面沿展向壓力系數(shù)分布的對(duì)比Fig.3 Comparison of spanwise pressure coefficient distributions between before and after adaptation with experimental data for ONERA-M6 wing

與RAE2822翼型一樣，ONERA-M6機(jī)翼也是一個(gè)經(jīng)典的驗(yàn)證算例［12］。這里，選取的計(jì)算狀態(tài)為:M∞=0.84，攻角α=3.06°。初始計(jì)算網(wǎng)格為純四面體網(wǎng)格，包含52144個(gè)網(wǎng)格點(diǎn)，296517個(gè)四面體單元，其中，機(jī)翼物面分布的網(wǎng)格點(diǎn)數(shù)為5459。圖2(a)給出了初始的機(jī)翼物面網(wǎng)格?？梢钥吹剑跏嫉奈锩婢W(wǎng)格點(diǎn)在機(jī)翼表面分布比較均勻，除在翼尖、機(jī)翼前緣位置從幾何保形的角度作了處理外，沒(méi)有從流動(dòng)特征方面作任何特殊考慮。同樣，選用升力系數(shù)作為目標(biāo)輸出函數(shù)進(jìn)行自適應(yīng)，設(shè)定全局誤差閥值為0.1。自適應(yīng)后的四面體網(wǎng)格共有128386個(gè)點(diǎn)，740366個(gè)四面體單元，其中機(jī)翼物面分布的網(wǎng)格點(diǎn)數(shù)為8041。圖2(b)給出自適應(yīng)后物面網(wǎng)格的示意圖。從自適應(yīng)的物面網(wǎng)格點(diǎn)分布可以很清晰看到一個(gè)呈現(xiàn)λ形的加密區(qū)域。這也正好對(duì)應(yīng)本計(jì)算狀態(tài)下ONERA-M6機(jī)翼表面應(yīng)呈現(xiàn)的激波結(jié)構(gòu)。從圖中可看出，通過(guò)自適應(yīng)，網(wǎng)格點(diǎn)對(duì)機(jī)翼前后緣也進(jìn)行了加密。圖2(c)給出由自適應(yīng)網(wǎng)格得到的ONERA-M6機(jī)翼物面馬赫數(shù)云圖，計(jì)算得到的機(jī)翼表面 激波非常清晰。圖3比較了自適應(yīng)前后計(jì)算得到的機(jī)翼物面沿展向壓力系數(shù)分布。從壓力系數(shù)(尤其在靠近翼根的位置)的比較，可以看到，對(duì)于ONERA-M6機(jī)翼，采用伴隨自適應(yīng)還可在一定程度上提高計(jì)算的精度。

3 復(fù)雜流場(chǎng)計(jì)算應(yīng)用

這一節(jié)討論伴隨自適應(yīng)技術(shù)在復(fù)雜流場(chǎng)計(jì)算中的應(yīng)用。選擇的考核算例是第二屆國(guó)際渦流試驗(yàn)項(xiàng)目(VFE-2)中的尖前緣65°三角翼模型［13］，主要考察伴隨自適應(yīng)技術(shù)在提高基于非結(jié)構(gòu)網(wǎng)格的CFD技術(shù)預(yù)測(cè)大攻角復(fù)雜渦流場(chǎng)的能力。VFE-2項(xiàng)目是美國(guó)和歐洲在2004年聯(lián)合發(fā)起的國(guó)際合作項(xiàng)目，采用NASA Langley研究中心提供的不同前緣鈍度的三角翼標(biāo)準(zhǔn)模型作為研究對(duì)象，旨在確認(rèn)和評(píng)估當(dāng)前CFD方法在復(fù)雜渦流場(chǎng)的預(yù)測(cè)能力和技術(shù)狀態(tài)。關(guān)于該模型的介紹和具體描述參見(jiàn)文獻(xiàn)［14］。本文選取的計(jì)算狀態(tài)是:M∞=0.4，攻角α=20.3°。對(duì)于粘性計(jì)算，基于平均氣動(dòng)弦長(zhǎng)的雷諾數(shù)Re=2.0×106。為了減小計(jì)算量，實(shí)際計(jì)算中，采用半翼展模型開(kāi)展計(jì)算研究。

3.1 無(wú)粘計(jì)算結(jié)果

初始網(wǎng)格采用純四面體網(wǎng)格，包含64835個(gè)網(wǎng)格點(diǎn)，365681個(gè)四面體單元，其中，物面上分布的的網(wǎng)格點(diǎn)為9464個(gè)。對(duì)于初始計(jì)算網(wǎng)格，除了從幾何保形方面考慮對(duì)三角翼前后緣適當(dāng)進(jìn)行加密控制外，沒(méi)有進(jìn)行任何特殊處理。圖4(a)上半部分給出了三角翼初始的物面網(wǎng)格，可以看到，物面上點(diǎn)的分布幾乎是均勻的。圖4(a)下半部分是由該初始計(jì)算網(wǎng)格采用無(wú)粘計(jì)算得到的物面壓力系數(shù)分布云圖，給出了一個(gè)很明顯的低壓區(qū)。這一低壓區(qū)正是三角翼上表面前緣渦掃過(guò)的痕跡，說(shuō)明由初始計(jì)算網(wǎng)格可捕捉到渦流場(chǎng)的存在。仍選用升力系數(shù)作為目標(biāo)輸出函數(shù)進(jìn)行伴隨自適應(yīng)，全局誤差閥值設(shè)為0.1。自適應(yīng)后的四面體網(wǎng)格共有139682個(gè)網(wǎng)格點(diǎn)，766948個(gè)四面體單元，其中，物面上分布的網(wǎng)格點(diǎn)為13333個(gè)。圖4(b)上半部分給出了自適應(yīng)后的三角翼物面網(wǎng)格示意圖。很明顯，自適應(yīng)加密的物面網(wǎng)格點(diǎn)主要分布在機(jī)翼翼根、機(jī)翼前后緣以及計(jì)算捕捉到的前緣渦所掃過(guò)的區(qū)域。圖4(b)下半部分給出了由自適應(yīng)網(wǎng)格計(jì)算得到的物面壓力系數(shù)分布云圖。與上圖比較，不難看到，自適應(yīng)后得到的低壓區(qū)更加清晰，這表明捕捉到的前緣渦的強(qiáng)度有所增加。

圖4 VFE-2尖前緣65°三角翼模型自適應(yīng)前后結(jié)果的對(duì)比Fig.4 Comparison of mesh and pressure coefficient distribution on the wall surface between before and after adaptation for VFE-2 delta wing with sharp leading-edge

圖5通過(guò)沿流向方向截取兩個(gè)典型站位x-c=0.4及x-c=0.6的空間網(wǎng)格，對(duì)比了自適應(yīng)前后空間網(wǎng)格的變化?？梢钥吹?，對(duì)于此算例，空間網(wǎng)格主要在靠近機(jī)翼上表面的位置(尤其是翼根和主渦的附近)進(jìn)行加密，而在機(jī)翼下表面附近網(wǎng)格幾乎沒(méi)有任何加密。并且，在流向方向的不同站位，網(wǎng)格的加密位置和密度隨主渦渦核的位置和強(qiáng)度有所變化，但并不完全集中在渦核周圍。這與基于流動(dòng)特征的渦自適應(yīng)或熵自適應(yīng)方法觀察到的有明顯差別［2］。

圖5 尖前緣65°三角翼模型自適應(yīng)前后兩個(gè)典型站位的空間網(wǎng)格對(duì)比Fig.5 Comparison of mesh resolution in space from two typical sections for VFE-2 delta wing with sharp leading-edge

圖6比較了自適應(yīng)前后采用無(wú)粘計(jì)算得到的尖前緣65°三角翼在典型站位的物面壓力系數(shù)分布?？梢钥吹?，自適應(yīng)前后，采用無(wú)粘計(jì)算在預(yù)測(cè)渦核的位置和主渦吸力峰強(qiáng)度方面都有很大改善，但與試驗(yàn)相比，整體上偏差都比較大。這也從側(cè)面印證了人們長(zhǎng)期形成的一個(gè)計(jì)算經(jīng)驗(yàn)，即采用Euler方程進(jìn)行無(wú)粘計(jì)算不足以很好地預(yù)測(cè)三角翼大攻角渦流場(chǎng)這類存在大的分離和以渦流為主導(dǎo)的復(fù)雜流場(chǎng)。

圖6 自適應(yīng)前后尖前緣65°三角翼模型在典型站位的物面壓力系數(shù)分布與試驗(yàn)值的對(duì)比(無(wú)粘計(jì)算)Fig.6 Comparison of pressure coefficient distributions between before and after adaptation with experimental data for VFE-2 delta wing with sharp leading-edge(in inviscid computation)

表1 兩套粘性計(jì)算網(wǎng)格的網(wǎng)格信息統(tǒng)計(jì)Table1 Information for two meshes in viscous computations

3.2 粘性計(jì)算結(jié)果

在粘性計(jì)算中，粘性計(jì)算網(wǎng)格分別在初始四面體網(wǎng)格和自適應(yīng)后的四面體網(wǎng)格基礎(chǔ)上生成。為了保證第一層網(wǎng)格的Y+≈1，邊界層內(nèi)第一層網(wǎng)格距離物面的距離按0.01進(jìn)行估計(jì)。在實(shí)際計(jì)算中，取第一層網(wǎng)格的實(shí)際間距為5.0×10-6，網(wǎng)格拉伸比為1.2，總的附面層網(wǎng)格層數(shù)為30層。兩套網(wǎng)格的詳細(xì)信息參見(jiàn)表1。

為了比較計(jì)算效果，圖7首先給出利用兩套粘性計(jì)算網(wǎng)格得到的物面極限流線和空間流場(chǎng)解的對(duì)比。計(jì)算中所有參數(shù)設(shè)置在前后保持一致。從物面極限流線的結(jié)果看，由自適應(yīng)局部加密的粘性網(wǎng)格得到的解不僅觀察到主再附線(說(shuō)明渦流場(chǎng)存在)，而且還清晰觀察到二次再附線，而由初始粘性網(wǎng)格得到的解僅能觀察到主再附線，二次再附線不明顯。并且，兩者主再附線的位置明顯有差別，由初始粘性網(wǎng)格得到的主再附線更靠近翼尖的位置。這表明，由自適應(yīng)局部加密的網(wǎng)格得到的主渦位置會(huì)明顯沿展向翼根靠近。從空間流場(chǎng)解的結(jié)果進(jìn)行分析，圖中采用沿流向方向不同站位的總壓比等值線和繞渦核的流線的方式清晰表現(xiàn)出渦流場(chǎng)的存在，但由自適應(yīng)局部加密的粘性網(wǎng)格得到的解明顯觀察到二次渦，并且主渦核的位置明顯更靠近翼根。因此，兩者的分析結(jié)論完全一致，即采用伴隨自適應(yīng)局部加密后，粘性計(jì)算進(jìn)行渦捕捉的能力明顯增強(qiáng)。另外，還容易注意到，由初始粘性網(wǎng)格計(jì)算得到的解還提前觀察到渦破裂。這與試驗(yàn)觀察到的情況不太吻合［15］。

圖7 尖前緣65°三角翼模型的物面極限流線和空間流場(chǎng)解對(duì)比(粘性計(jì)算)Fig.7 Comparison of flow solutions between before and after adaptation for VFE-2 delta wing with sharp leading-edge(in viscous computation)

圖8比較了粘性計(jì)算得到的物面壓力系數(shù)分布。與無(wú)粘計(jì)算僅給出三個(gè)典型站位的結(jié)果不同，圖8給出了在試驗(yàn)中所測(cè)量的所有站位數(shù)據(jù)的比較。其中，試驗(yàn)結(jié)果包含兩類，一類由NASA Langley中心的跨音速風(fēng)洞在雷諾數(shù)Re=6.0×106下進(jìn)行試驗(yàn)得到，一類由德國(guó)風(fēng)洞協(xié)會(huì)的跨音速風(fēng)洞在雷諾數(shù)Re=2.0×106下進(jìn)行試驗(yàn)得到，二者有較小的雷諾數(shù)效應(yīng)差別。可以看到，由初始粘性網(wǎng)格得到的各站位上表面壓力系數(shù)分布與試驗(yàn)值相差甚遠(yuǎn)，主要表現(xiàn)在壓力系數(shù)吸力峰位置明顯沿展向朝翼尖位置偏移，而吸力峰峰值比試驗(yàn)值偏小;與之相較，由自適應(yīng)局部加密的粘性網(wǎng)格除了在x/c=0.2這一站位外，在其他四個(gè)站位都給出與試驗(yàn)符合較好的結(jié)果。自適應(yīng)局部加密對(duì)計(jì)算結(jié)果的改進(jìn)非常明顯。在x/c=0.2這一站位，推測(cè)計(jì)算結(jié)果沒(méi)有明顯改善的原因是，兩套粘性網(wǎng)格在機(jī)翼前端附近的位置保持了基本一致的網(wǎng)格分辨率，而這一網(wǎng)格分辨率不足以保證在x/c=0.2站位給出與試驗(yàn)符合良好的計(jì)算結(jié)果?？赡艿慕鉀Q辦法是在這一區(qū)域附近進(jìn)行手動(dòng)加密以保證在機(jī)翼前端都有足夠的網(wǎng)格分布。

圖8 尖前緣65°三角翼模型在典型戰(zhàn)位的物面壓力系數(shù)分布與試驗(yàn)值的對(duì)比(粘性計(jì)算)Fig.8 Comparison of pressure coefficient distributions with experimental data for VFE-2 delta wing with sharp leading-edge(in viscous computation)

4 結(jié)論

對(duì)基于伴隨方程的網(wǎng)格自適應(yīng)技術(shù)的基本原理進(jìn)行了研究，構(gòu)造和實(shí)現(xiàn)了基于流場(chǎng)伴隨方程的非結(jié)構(gòu)網(wǎng)格自適應(yīng)探測(cè)器，采用局部加密網(wǎng)格的方法，建立了基于Euler方程的無(wú)粘伴隨自適應(yīng)能力。對(duì)粘性計(jì)算，提出首先采用無(wú)粘計(jì)算和伴隨自適應(yīng)生成網(wǎng)格分辨率好的初始網(wǎng)格，然后采用層推進(jìn)方法生成粘性計(jì)算網(wǎng)格的求解策略。采用二維和三維經(jīng)典算例對(duì)方法進(jìn)行了驗(yàn)證，并應(yīng)用本文提出的求解策略對(duì)VFE-2尖前緣三角翼大攻角渦流場(chǎng)進(jìn)行了數(shù)值模擬。計(jì)算結(jié)果表明，采用自適應(yīng)技術(shù)能提高計(jì)算精度，尤其對(duì)三角翼大攻角渦流場(chǎng)這類存在大的分離和以渦主導(dǎo)的復(fù)雜流場(chǎng)，采用本文策略，可明顯改善預(yù)測(cè)結(jié)果，有效提高預(yù)測(cè)準(zhǔn)確度。粘性計(jì)算時(shí)，在網(wǎng)格生成階段，采用Euler無(wú)粘計(jì)算進(jìn)行伴隨自適應(yīng)生成具有較好網(wǎng)格分辨率的初始四面體網(wǎng)格還可在一定程度上減小計(jì)算代價(jià)，提高計(jì)算效率。

［1］張來(lái)平，王志堅(jiān).三維復(fù)雜外形的非結(jié)構(gòu)網(wǎng)格自動(dòng)生成技術(shù)與應(yīng)用［J］.計(jì)算物理，1999，16(5):552-558.

［2］BAI W，QIU Z，LI L.Recent efforts to establish adaptive hybrid grid computing capability at ACTRI［J］.Computational Fluid Dynamics journal，2007，7(4):438-449.

［3］WARREN G P，ANDERSON J T，THOMAS J T，KRIST S L.Grid convergence for adaptive methods［R］.AIAA paper 1991-1592，1991.

［4］PARK M A.Adjoint-based three-dimensional error prediction and grid adaptation［R］.AIAA paper 2002-3286，2002.

［5］NIELSEN E J.Adjoint-based algorithms for adaptation and design optimization on unstructured grids［A］.The 3rd East-West High-Speed Flow-field Conference［C］.Beijing，China，2005.

［6］LEE-RAUSCH E A，PARK M A，JONES W T，HAMMOND D P，NIELSEN E J.Application ofparalleladjoint-based error estimation and anisotropic grid adaptation for three-dimensionalaerospace configurations［A］.The 23rd AIAA Applied Aerodynamics Conference［C］.Toronto，Ontario，2005.

［7］VENDITTI D A.Grid adaptation for functional outputs of compressible flow simulations［D］.Massachusetts Institute of Technology，Cambridge，MA，2002.

［8］PIRZADEH S Z.An adaptive unstructured grid method by grid subdivision，local remeshing，and grid movement［R］.AIAA paper 1999-3255，1999.

［9］SHAROV D，F(xiàn)UJII K.Three-dimensional adaptive bisection of unstructured grids for transientcompressible flow computations［R］.AIAA paper 1995-1708，1995.

［10］TIMOTHY J B.Mesh generation:art or science［J］.Progress in Aerospace Sciences，2005，41:29-63.

［11］http://www.grc.nasa.gov/www/wind/valid/raetaf/raetaf.html

［12］http://www.grc.nasa.gov/www/wind/valid/m6wing/m6wing.html

［13］UMMEL D，REDEKER G.A new vortex flow experimentfor computer code validation［A］.RTO Symposium on"Advanced Flow Management:Part A-Vortex Flows and High Angle of Attack for Military Vehicles"［C］.Loen，Norway，2001.

［14］CHU J，LUCKRING J M.Experimental surface pressure data obtained on 65°delta wing across Reynolds number and Mach number ranges［R］.NASA TM-1996-4645，1996.

［15］KONRATH R，SCH?DER A，KOMPENHAS J.Analysis of PIV results obtained for the VFE-265°delta wing configuration at sub-and transonic speeds［A］.24th Applied Aerodynamics Conference［C］.San Francisco，California，2006.