唐靜，崔鵬程，賈洪印,*，李彬，李歡

1. 中國(guó)空氣動(dòng)力研究與發(fā)展中心 計(jì)算空氣動(dòng)力研究所， 綿陽(yáng) 621000

2. 西北工業(yè)大學(xué) 航空學(xué)院，西安 710072

近20年來(lái)，CFD已成為研究空氣動(dòng)力學(xué)的主要手段之一，并獲得了工業(yè)部門(mén)的廣泛認(rèn)可。然而，CFD在飛行器跨聲速氣動(dòng)特性計(jì)算、大迎角失速特性預(yù)測(cè)等方面仍面臨著較大的挑戰(zhàn)。CFD計(jì)算結(jié)果受數(shù)值格式、湍流模型和計(jì)算網(wǎng)格的影響很大。工程上常用的數(shù)值格式和湍流模型相對(duì)比較固定，且通過(guò)大量的研究和廣泛的應(yīng)用獲得了較為詳細(xì)的參數(shù)影響規(guī)律。然而，由于針對(duì)復(fù)雜外形的網(wǎng)格生成數(shù)學(xué)理論的欠缺，造成網(wǎng)格生成算法不完善，網(wǎng)格生成仍需要較多的人工干預(yù)。計(jì)算網(wǎng)格已成為影響CFD模擬精準(zhǔn)度的關(guān)鍵因素[1]，顯著影響著飛行器氣動(dòng)、控制和動(dòng)力等性能的預(yù)測(cè)精度。

網(wǎng)格生成的單元疏密分布依賴(lài)于對(duì)流動(dòng)特征的先驗(yàn)估計(jì)。工程問(wèn)題中流動(dòng)狀態(tài)通常很多，針對(duì)每個(gè)流動(dòng)狀態(tài)生成對(duì)應(yīng)網(wǎng)格的工作量太大。因此，與流場(chǎng)耦合、在計(jì)算過(guò)程中自動(dòng)優(yōu)化網(wǎng)格疏密分布的網(wǎng)格自適應(yīng)技術(shù)得到了國(guó)內(nèi)外研究人員的廣泛研究。美國(guó)NASA的報(bào)告《2030年CFD展望研究》明確將網(wǎng)格自適應(yīng)技術(shù)列為亟待發(fā)展的一項(xiàng)關(guān)鍵技術(shù)[2]。

網(wǎng)格自適應(yīng)技術(shù)的研究大約起步于20世紀(jì)80年代[3]，經(jīng)過(guò)近40年的發(fā)展，在網(wǎng)格單元誤差估計(jì)[4-5]、網(wǎng)格分布優(yōu)化[6]、物面網(wǎng)格幾何投影等方面都取得了較大的進(jìn)展[7]。然而，由于魯棒性不足(尤其是附面層內(nèi)半結(jié)構(gòu)化網(wǎng)格)和幾何保形實(shí)現(xiàn)難度大等問(wèn)題，網(wǎng)格自適應(yīng)技術(shù)還未能廣泛應(yīng)用于工程問(wèn)題中[7-8]。

網(wǎng)格分布優(yōu)化通常包括r型自適應(yīng)和h型自適應(yīng)兩種類(lèi)型。r型自適應(yīng)僅移動(dòng)網(wǎng)格點(diǎn)，網(wǎng)格的拓?fù)浣Y(jié)構(gòu)和單元連接關(guān)系在自適應(yīng)過(guò)程中保持不變。由于不增加網(wǎng)格單元的總量，r型自適應(yīng)改變網(wǎng)格分布的能力有限。h型自適應(yīng)則是通過(guò)局部區(qū)域網(wǎng)格單元的剖分加密或合并稀疏的方式優(yōu)化網(wǎng)格分布，以提高關(guān)鍵區(qū)域流動(dòng)模擬的精度。

四面體網(wǎng)格生成算法相對(duì)成熟[9-10]，因而四面體網(wǎng)格h自適應(yīng)方法的成熟度較高[11]，還可以給定網(wǎng)格尺寸場(chǎng)開(kāi)展各向異性的自適應(yīng)[12]。四面體網(wǎng)格主要應(yīng)用于無(wú)黏流動(dòng)，當(dāng)擴(kuò)展到黏性流動(dòng)時(shí)，物面附近需要采用層結(jié)構(gòu)的各向異性四面體網(wǎng)格單元[13]。黏性流動(dòng)的模擬更多地采用非結(jié)構(gòu)混合網(wǎng)格，其h型自適應(yīng)技術(shù)的難點(diǎn)集中在附面層內(nèi)層結(jié)構(gòu)網(wǎng)格的處理。

Alrutz[14]采用固定層結(jié)構(gòu)網(wǎng)格的策略，只針對(duì)四面體網(wǎng)格開(kāi)展自適應(yīng)，不能改變附面層網(wǎng)格的分辨率。Joubarne[15]和Lepage[16]等對(duì)四面體采用h自適應(yīng)，對(duì)附面層網(wǎng)格則采用y+自適應(yīng)，只能在法向高度方向改變網(wǎng)格的分辨率。Park[17]將棱柱和金字塔單元都轉(zhuǎn)化成四面體進(jìn)行處理，李立等[18]在純四面體網(wǎng)格上開(kāi)展h自適應(yīng)后再重新生成層結(jié)構(gòu)網(wǎng)格，這兩種方法可改變層結(jié)構(gòu)網(wǎng)格的流向尺寸分布，但前者失去了六面體或三棱柱網(wǎng)格模擬黏性流動(dòng)的精度優(yōu)勢(shì)，后者層結(jié)構(gòu)網(wǎng)格的生成原本就是網(wǎng)格生成的難點(diǎn)。

為了同時(shí)改變層結(jié)構(gòu)網(wǎng)格物面法向和切向的尺寸分布，更有效的方法是采用網(wǎng)格剖分方法。按自適應(yīng)后網(wǎng)格單元類(lèi)型，自適應(yīng)可以分為標(biāo)準(zhǔn)單元類(lèi)型、懸空節(jié)點(diǎn)類(lèi)型和多面體單元類(lèi)型。標(biāo)準(zhǔn)單元類(lèi)型方法需要剖分臨近單元以消除懸空節(jié)點(diǎn)，層結(jié)構(gòu)內(nèi)的網(wǎng)格通常只能沿著物面法向方向剖分并“貫穿”整個(gè)層結(jié)構(gòu)以改變物面切向方向的網(wǎng)格尺寸[19-20]，改變法向尺寸多采用y+方法[21]。該類(lèi)方法自適應(yīng)后網(wǎng)格保持了完整的層結(jié)構(gòu)，但要求初始網(wǎng)格的層結(jié)構(gòu)清晰，且在處理金字塔單元時(shí)可能出現(xiàn)剖分方式不能自洽的問(wèn)題。

懸空節(jié)點(diǎn)類(lèi)方法對(duì)網(wǎng)格單元剖分方向不作任何限制，但要求流場(chǎng)解算器能處理懸空節(jié)點(diǎn)。該類(lèi)方法可分為基于網(wǎng)格面和基于網(wǎng)格單元兩類(lèi)?；诰W(wǎng)格面的方法[22-23]可在一定程度上實(shí)現(xiàn)網(wǎng)格各向異性分布，但通常網(wǎng)格單元的幾何質(zhì)量隨著自適應(yīng)迭代下降很快?；诰W(wǎng)格單元的方法既可用于各向異性的網(wǎng)格自適應(yīng)[24]，也可用于各向同性的網(wǎng)格自適應(yīng)[25]，后者網(wǎng)格單元剖分類(lèi)型很少，對(duì)數(shù)據(jù)結(jié)構(gòu)的依賴(lài)較輕，軟件實(shí)現(xiàn)更容易。

多面體單元類(lèi)型的自適應(yīng)[26-27]將帶有懸空節(jié)點(diǎn)的標(biāo)準(zhǔn)單元轉(zhuǎn)換成多面體單元，不需要對(duì)相鄰單元進(jìn)行額外的處理，消除了網(wǎng)格單元之間的相互影響，魯棒性更高。多面體類(lèi)型也可類(lèi)似地分成基于網(wǎng)格面[26]和基于網(wǎng)格單元兩種子類(lèi)型[28]。

表面網(wǎng)格加密后，新增網(wǎng)格點(diǎn)通常不在幾何曲面上，需要將其投影到物面上。物面幾何信息可通過(guò)幾何系統(tǒng)(CAD)直接獲得[28]，也可通過(guò)初始表面網(wǎng)格擬合成的近似曲面得到[29-30]。前者能將新增網(wǎng)格點(diǎn)精確投影到真實(shí)幾何上，但對(duì)于復(fù)雜外形容易出現(xiàn)投影多解或網(wǎng)格交叉等問(wèn)題。擬合曲面與真實(shí)幾何有一定的偏差，但消除了自適應(yīng)方法對(duì)CAD系統(tǒng)的依賴(lài)。依據(jù)采用網(wǎng)格單元的模板大小，擬合方法可分為局部曲面擬合[29]和非局部曲面擬合[30]兩種類(lèi)型。非局部曲面擬合需要使用較大區(qū)域內(nèi)的表面網(wǎng)格信息和采用分片高次擬合。局部曲面擬合只使用當(dāng)前表面網(wǎng)格單元的頂點(diǎn)信息，如局部Coons曲面[29]。

一般情況下，物面附近的層結(jié)構(gòu)網(wǎng)格單元法向尺寸遠(yuǎn)小于切向尺寸，表面網(wǎng)格點(diǎn)發(fā)生微小的法向移動(dòng)都可能造成物面附近的網(wǎng)格單元出現(xiàn)重疊交錯(cuò)等嚴(yán)重問(wèn)題。因此，表面網(wǎng)格投影后需要對(duì)體網(wǎng)格單元進(jìn)行重構(gòu)或移動(dòng)變形，確保計(jì)算網(wǎng)格有效。復(fù)雜外形的棱柱網(wǎng)格重構(gòu)難度很大，且魯棒性通常不高。通過(guò)移動(dòng)體網(wǎng)格頂點(diǎn)的網(wǎng)格變形技術(shù)更加高效，但要求采用的網(wǎng)格變形技術(shù)能有效處理棱柱層網(wǎng)格[31]。

本文以工程實(shí)用為目標(biāo)，發(fā)展了網(wǎng)格自適應(yīng)系統(tǒng)涉及的網(wǎng)格分布優(yōu)化、網(wǎng)格幾何投影和空間網(wǎng)格匹配3項(xiàng)關(guān)鍵技術(shù)，建立了魯棒的非結(jié)構(gòu)混合網(wǎng)格自適應(yīng)方法。首先，在多面體網(wǎng)格自適應(yīng)方法的基礎(chǔ)上，發(fā)展了標(biāo)準(zhǔn)網(wǎng)格面的網(wǎng)格分布優(yōu)化方法，提高了方法的通用性，解決了物面網(wǎng)格投影造成的“縫隙”問(wèn)題。其次，通過(guò)Hermite插值方法構(gòu)造了局部Coons曲面，采用參數(shù)化投影方法實(shí)現(xiàn)了表面網(wǎng)格的快速投影。最后，采用改進(jìn)的距離函數(shù)動(dòng)網(wǎng)格方法實(shí)現(xiàn)了空間網(wǎng)格與投影后表面網(wǎng)格的協(xié)調(diào)匹配。

1 網(wǎng)格單元加密或稀疏的自適應(yīng)判據(jù)

流場(chǎng)是構(gòu)造網(wǎng)格單元自適應(yīng)判據(jù)的基礎(chǔ)，本文流場(chǎng)求解采用自主研發(fā)的MFlow軟件[32]。該軟件采用二階格式的有限體積方法求解積分型的雷諾平均Navier-Stokes方程。文中對(duì)流通量采用Roe格式[33]，黏性項(xiàng)離散采用中心差分，時(shí)間項(xiàng)離散采用LU-SGS方法[34]，湍流模擬采用一方程SA模型[35]或RC修正形式的SA模型[36]。

1.1 擬梯度自適應(yīng)判據(jù)

擬梯度自適應(yīng)判據(jù)的理論依據(jù)假定流場(chǎng)梯度大的區(qū)域計(jì)算誤差更大。本文采用密度ρ、壓力p、速率|v|和溫度T4個(gè)流場(chǎng)變量加權(quán)的形式。為了避免流場(chǎng)變量數(shù)值量級(jí)差別大造成權(quán)值取值困難的問(wèn)題，本文采用遠(yuǎn)場(chǎng)流場(chǎng)值(下標(biāo)∞表示)進(jìn)行無(wú)量綱化，則加權(quán)流動(dòng)變量可表示為

(1)

式中：|v|為流動(dòng)速度矢量；w1～w4分別為4個(gè)變量的權(quán)值。進(jìn)一步考慮網(wǎng)格尺寸因素，本文擬梯度計(jì)算采用的形式為

(2)

式中:下標(biāo)1、2分別指網(wǎng)格面兩側(cè)的單元;d12為體心間的距離; 上標(biāo)s為距離的影響因子，s值越大距離影響越大，一般取值范圍為-1.0～1.3。特別地，當(dāng)s取-1時(shí)，類(lèi)似于采用近似一階數(shù)值方法計(jì)算梯度。則相鄰網(wǎng)格單元中最大梯度值gm為

(3)

式中:NF為面相鄰網(wǎng)格單元的總數(shù)。

遍歷所有網(wǎng)格單元得到gm的最大值gmax，給定自適應(yīng)加密和粗化的閾值參數(shù)tr和tc，則加密和粗化的臨界梯度值為

gr=gmaxtr,gc=gmaxtc

(4)

當(dāng)某網(wǎng)格單元的gm大于gr時(shí)標(biāo)記為加密，當(dāng)網(wǎng)格單元的gm小于gc時(shí),標(biāo)記為粗化，其他的網(wǎng)格單元保持不變。

1.2 剪切應(yīng)變率自適應(yīng)判據(jù)

剪切應(yīng)變率自適應(yīng)判據(jù)[37]主要用來(lái)識(shí)別流場(chǎng)中的分離旋渦，由Q判據(jù)進(jìn)行無(wú)量綱化得到,即

(5)

式中：ω和S分別為速度梯度的旋轉(zhuǎn)分量和剪切分量，定義為

(6)

由于自由來(lái)流的S值相對(duì)很小，剪切應(yīng)變率 計(jì)算可能會(huì)產(chǎn)生數(shù)值噪聲，造成旋渦區(qū)域的誤判。本文采用應(yīng)變過(guò)濾方法消除數(shù)值噪聲，則判定網(wǎng)格單元處于旋渦區(qū)域還需要滿(mǎn)足如下條件：

(7)

式中：ε為過(guò)濾閾值；N為網(wǎng)格單元總數(shù)；h取值0.1～0.5。

2 網(wǎng)格單元分布優(yōu)化方法

本文采用剖分或合并網(wǎng)格單元的方法實(shí)現(xiàn)網(wǎng)格分布的優(yōu)化。為提高網(wǎng)格加密或粗化操作的運(yùn)行效率，應(yīng)建立合適的數(shù)據(jù)結(jié)構(gòu)。

2.1 網(wǎng)格細(xì)化方法

本文綜合考慮加密方法的魯棒性和效率，選擇各向同性加密網(wǎng)格單元，則每類(lèi)單元對(duì)應(yīng)著唯一的加密模式。并且，除金字塔單元外，其他類(lèi)型網(wǎng)格單元加密后的子網(wǎng)格單元具有自相似性，即子單元類(lèi)型與父級(jí)單元相同。各類(lèi)型的網(wǎng)格單元加密方法如下：

1) 三角形剖分為4個(gè)三角形，見(jiàn)圖1(a)。

2) 四邊形剖分為4個(gè)四邊形，見(jiàn)圖1(b)。

3) 四面體剖分為8個(gè)四面體，見(jiàn)圖1(c)。

4) 三棱柱剖分為8個(gè)三棱柱，見(jiàn)圖1(d)。

5) 六面體剖分為8個(gè)六面體，見(jiàn)圖1(e)。

6) 金字塔剖分為6個(gè)金字塔和4個(gè)四面體，見(jiàn)圖1(f)。

圖1 非結(jié)構(gòu)混合網(wǎng)格單元加密模式

定義1自適應(yīng)層級(jí)某網(wǎng)格單元的自適應(yīng)層級(jí)是指由初始網(wǎng)格單元得到該網(wǎng)格單元需要加密的次數(shù)。

相鄰兩個(gè)單元層級(jí)不同時(shí)，在兩個(gè)單元的交界面或交界線上將出現(xiàn)圖2所示的懸空節(jié)點(diǎn)。由于常規(guī)流場(chǎng)模擬軟件不支持懸空節(jié)點(diǎn)，為了提高網(wǎng)格自適應(yīng)方法對(duì)流場(chǎng)模擬軟件的通用性，采用多面體轉(zhuǎn)換方法間接消除懸空節(jié)點(diǎn)。以圖2左側(cè)單元為例，轉(zhuǎn)化為多面體后共包括9個(gè)面：F1:{1,2,6,5},F2:{4,a,c,b},F3:{a,3,d,c},F4:{b,c,e,8},F5:{c,d,7,e},F6:{1,2,3,a,4},F7:{2,3,d,7,6},F8:{6,7,e,8,5},F9:{1,4,b,8,5}。其中，前5個(gè)面為標(biāo)準(zhǔn)的四邊形，后4個(gè)面為包含5個(gè)頂點(diǎn)的多邊形。一般情況下，基于非結(jié)構(gòu)網(wǎng)格的流場(chǎng)模擬軟件只支持三角形和四邊形兩種標(biāo)準(zhǔn)網(wǎng)格面單元。因此，為了進(jìn)一步提高自適應(yīng)系統(tǒng)的通用性，本文將多邊形面轉(zhuǎn)換成多個(gè)三角形和四邊形。對(duì)于三角形面，標(biāo)準(zhǔn)化轉(zhuǎn)換共有3種拓?fù)漕?lèi)型，如圖3所示；對(duì)于四邊形面，共有5種拓?fù)漕?lèi)型，如圖4所示。

圖2中定義多面體的后4個(gè)多邊形網(wǎng)格面都包括5個(gè)網(wǎng)格頂點(diǎn)，每個(gè)多邊形標(biāo)準(zhǔn)化后剖分成3個(gè)三角形單元，因此該多面體單元標(biāo)準(zhǔn)化后最終包括17個(gè)面，如圖5所示。

圖2 懸空節(jié)點(diǎn)示意圖

圖3 非標(biāo)準(zhǔn)三角形網(wǎng)格面拓?fù)漕?lèi)型

圖4 非標(biāo)準(zhǔn)四邊形網(wǎng)格面拓?fù)漕?lèi)型

圖5 多面體單元多邊形網(wǎng)格面標(biāo)準(zhǔn)化示意圖

2.2 網(wǎng)格粗化方法

本文網(wǎng)格粗化方法采用“回退”方法，即網(wǎng)格粗化是網(wǎng)格加密的逆向過(guò)程。因而只有當(dāng)某網(wǎng)格單元經(jīng)過(guò)加密后生成的所有子網(wǎng)格單元都被探測(cè)器標(biāo)記為粗化時(shí)，該網(wǎng)格才回退到加密前狀態(tài)。

由于相鄰網(wǎng)格單元體積的比值過(guò)大會(huì)降低CFD求解的收斂速度和穩(wěn)定性，因此通常需要對(duì)相鄰單元的層級(jí)差進(jìn)行限制。另外，同時(shí)具備網(wǎng)格細(xì)化和粗化功能時(shí)，自適應(yīng)屬性的標(biāo)記可能出現(xiàn)矛盾情況，即同一網(wǎng)格單元被探測(cè)器標(biāo)記為粗化，卻由于層級(jí)差限制需要加密。因此，本文自適應(yīng)方法增加了自適應(yīng)加密/粗化準(zhǔn)則。

準(zhǔn)則1自適應(yīng)加密/粗化準(zhǔn)則某網(wǎng)格單元自適應(yīng)加密或粗化屬性的標(biāo)記滿(mǎn)足如下3條規(guī)則：

1) 層級(jí)差限制規(guī)則。共面相鄰的網(wǎng)格單元層級(jí)差不超過(guò)1，否則，加密相鄰層級(jí)較低的網(wǎng)格單元。

2) 網(wǎng)格光滑性規(guī)則。當(dāng)與某個(gè)單元共面的所有單元都標(biāo)記為加密屬性時(shí)，該單元也標(biāo)記為加密屬性。

3) 加密優(yōu)先規(guī)則。當(dāng)某個(gè)網(wǎng)格單元標(biāo)記為粗化卻因規(guī)則1)或規(guī)則2)需要標(biāo)記為加密時(shí)，加密標(biāo)記優(yōu)先于粗化。

2.3 網(wǎng)格自適應(yīng)的數(shù)據(jù)結(jié)構(gòu)

由于本文網(wǎng)格粗化方法采用網(wǎng)格加密的逆向操作，因此從初始網(wǎng)格到最細(xì)網(wǎng)格，包括中間過(guò)程的網(wǎng)格，都需要記錄和存儲(chǔ)。各層級(jí)的網(wǎng)格還需要建立父級(jí)-子級(jí)的對(duì)應(yīng)關(guān)系，以便實(shí)現(xiàn)粗化操作。多級(jí)的粗、細(xì)網(wǎng)格形成自然的樹(shù)形數(shù)據(jù)關(guān)系，并且由各類(lèi)單元加密模式可知父-子級(jí)網(wǎng)格單元是一對(duì)多的關(guān)系，因此可以采用多級(jí)多路的樹(shù)形數(shù)據(jù)結(jié)構(gòu)處理由最粗網(wǎng)格單元到最細(xì)網(wǎng)格單元之間的層級(jí)關(guān)系。

3 表面網(wǎng)格幾何投影方法

實(shí)際應(yīng)用中物面通常為曲面，因此表面網(wǎng)格加密后新增加的網(wǎng)格點(diǎn)通常不在物面上，需要將其投影到物面或重構(gòu)的近似物面上。

3.1 物面幾何重構(gòu)

本文物面幾何重構(gòu)方法采用局部Coons曲面擬合，該方法只依賴(lài)于表面網(wǎng)格點(diǎn)坐標(biāo)和相應(yīng)的法向矢量。法向單位矢量通常不是標(biāo)準(zhǔn)的網(wǎng)格數(shù)據(jù)，本文通過(guò)與網(wǎng)格點(diǎn)直接相連的面網(wǎng)格的法向加權(quán)得到，即網(wǎng)格點(diǎn)i的法向向量ni為

(8)

式中：nk和Sk分別為網(wǎng)格面k的單位法向量和面積；Nk為與網(wǎng)格點(diǎn)i相連的網(wǎng)格面總數(shù)。

局部Coons曲面擬合的基本思想是先通過(guò)表面網(wǎng)格點(diǎn)和法向信息重構(gòu)網(wǎng)格邊對(duì)應(yīng)的三維邊界曲線，然后通過(guò)邊界曲線的融合構(gòu)造三維曲面。

本文采用三次Hermit插值直接構(gòu)造三維邊界曲線。如圖6所示，表面網(wǎng)格邊兩個(gè)頂點(diǎn)的坐標(biāo)為x1和x2，相應(yīng)法向?yàn)閚1和n2，參數(shù)化的邊界曲線c(t)為

c(t)=f1x1+f2x2+g1τ1+g2τ2

(9)

式中：f和g為權(quán)值；τ為相應(yīng)切向修正量，即

d1=x2-x1,d2=x1-x2

(10)

針對(duì)三角形和四邊形兩種表面網(wǎng)格單元，邊界曲線構(gòu)造方法相同，但曲線融合成曲面的方法有較大差別，本文基于文獻(xiàn)[29]中的方法，針對(duì)三角形和四邊形分別采用如下的融合方法。

對(duì)于三角形網(wǎng)格面，如圖7所示，3個(gè)網(wǎng)格點(diǎn)為x1～x3，相應(yīng)法向?yàn)閚1～n3。首先將三角形x1x2x3轉(zhuǎn)化到標(biāo)準(zhǔn)三角形參數(shù)平面V1V2V3，V1V2方向參數(shù)為S′，V1V3方向參數(shù)為S″。網(wǎng)格邊x1x2構(gòu)成邊界曲線c1，網(wǎng)格邊x2x3構(gòu)成邊界曲線c2，網(wǎng)格邊x1x3構(gòu)成邊界曲線c3，則三維擬合曲面S(S′,S″)表示為

S(S′,S″)=0.5(Xa+Xb+Xc-Xt)

(11)

圖6 三維邊界曲線重構(gòu)示意圖(虛線：網(wǎng)格邊，實(shí)線：重構(gòu)的邊界曲線)

圖7 三角形網(wǎng)格曲面重構(gòu)示意圖

圖8 四邊形網(wǎng)格曲面重構(gòu)示意圖

對(duì)于四邊形網(wǎng)格面，如圖8所示，4個(gè)網(wǎng)格點(diǎn)為x1～x4，相應(yīng)法向?yàn)閚1～n4。類(lèi)似地，首先將四邊形x1x2x3x4轉(zhuǎn)化到標(biāo)準(zhǔn)四邊形參數(shù)平面V1V2V3V4，V1V2方向參數(shù)為S′，V1V4方向參數(shù)為S″。網(wǎng)格邊x1x2構(gòu)成邊界曲線c1，網(wǎng)格邊x2x3構(gòu)成邊界曲線c2，網(wǎng)格邊x4x3構(gòu)成邊界曲線c3，網(wǎng)格邊x1x4構(gòu)成邊界曲線c4，則三維的擬合曲面S(S′,S″)表示為

S(S′,S″)=Xa+Xb-Xt

(12)

3.2 新增網(wǎng)格點(diǎn)投影

本文擬合曲面為參數(shù)曲面，可以采用參數(shù)點(diǎn)映射方法實(shí)現(xiàn)新增表面網(wǎng)格點(diǎn)的幾何投影，即先在標(biāo)準(zhǔn)參數(shù)面上取投影點(diǎn)的參數(shù)坐標(biāo)，然后代入?yún)?shù)曲面方程求得投影點(diǎn)在曲面上的坐標(biāo)點(diǎn)，最后將新增網(wǎng)格點(diǎn)移動(dòng)到目標(biāo)點(diǎn)位置。

如圖9所示，對(duì)于三角形網(wǎng)格面，常用的新增表面網(wǎng)格點(diǎn)的位置為N1、N2和N3，對(duì)應(yīng)參數(shù)坐標(biāo)為

圖9 網(wǎng)格參數(shù)曲面和新增點(diǎn)參數(shù)

(13)

對(duì)于四邊形網(wǎng)格面，常用的新增點(diǎn)位置為N1、N2、N3、N4和N5，對(duì)應(yīng)參數(shù)坐標(biāo)為

(14)

4 空間網(wǎng)格匹配方法

因模擬黏性流動(dòng)的需要，附面層內(nèi)緊鄰物面網(wǎng)格的法向尺寸遠(yuǎn)小于其他兩個(gè)方向的尺寸，新增表面網(wǎng)格點(diǎn)投影后，極有可能造成物面附近棱柱網(wǎng)格單元出現(xiàn)重疊交錯(cuò)的問(wèn)題，如圖10所示，因此必須對(duì)空間網(wǎng)格點(diǎn)進(jìn)行處理。網(wǎng)格變形技術(shù)不改變網(wǎng)格拓?fù)浜筒辉鰷p網(wǎng)格點(diǎn)數(shù)，通用性好，其中計(jì)算速度最快的一類(lèi)方法是基于距離函數(shù)的動(dòng)網(wǎng)格方法。

圖10 表面網(wǎng)格投影造成網(wǎng)格單元重疊交錯(cuò)示意圖

本文采用改進(jìn)型的基于距離函數(shù)的動(dòng)網(wǎng)格技術(shù)實(shí)現(xiàn)表面網(wǎng)格投影后體網(wǎng)格單元的快速匹配。該方法通過(guò)引入自調(diào)節(jié)的影響半徑大大提高了附面層內(nèi)半結(jié)構(gòu)化網(wǎng)格的變形能力，同時(shí)通過(guò)包圍盒加速技術(shù)提高了變形方法的計(jì)算效率。

給定表面網(wǎng)格點(diǎn)變形量Δxi,S，體網(wǎng)格點(diǎn)坐標(biāo)變化量Δxj,V表示為

(15)

式中：NS為表面網(wǎng)格點(diǎn)總數(shù)；N為體網(wǎng)格點(diǎn)總數(shù)；η為衰減函數(shù)；ωi為表面網(wǎng)格點(diǎn)i的權(quán)值。衰減函數(shù)和權(quán)值的取值方法以及更多細(xì)節(jié)信息可參考文獻(xiàn)[31]。

5 網(wǎng)格自適應(yīng)方法數(shù)值驗(yàn)證

5.1 30P30N三段翼構(gòu)型繞流自適應(yīng)模擬

30P30N三段翼構(gòu)型廣泛應(yīng)用于CFD數(shù)值方法的驗(yàn)證，主要考察數(shù)值方法對(duì)尾跡剪切層的模擬能力。圖11給出了其幾何構(gòu)型[38]和風(fēng)洞試驗(yàn)選擇的6個(gè)速度型測(cè)量位置。計(jì)算來(lái)流狀態(tài)為：?jiǎn)挝焕字Z數(shù)Re=1.61×107，馬赫數(shù)Ma∞=0.2，靜溫T∞=288.15 K，迎角α=19°，采用的湍流模型為標(biāo)準(zhǔn)一方程SA模型。計(jì)算采用六面體為主的非結(jié)構(gòu)混合網(wǎng)格，如圖12所示，其翼型主翼周向共251個(gè)網(wǎng)格點(diǎn)，前緣縫翼和后緣襟翼分別有131和151個(gè)網(wǎng)格點(diǎn)。法向第1層網(wǎng)格高度為2×10-5m，網(wǎng)格增長(zhǎng)比為1.12，展向共2個(gè)網(wǎng)格點(diǎn)。

自適應(yīng)計(jì)算共迭代3次，自適應(yīng)判據(jù)采用基于速率的擬梯度方法，對(duì)應(yīng)的參數(shù)tr=7×10-5，tc=0，s=-1.0，加權(quán)流場(chǎng)變量速率項(xiàng)權(quán)值為1，其他項(xiàng)為0。圖13給出了自適應(yīng)后的網(wǎng)格，由于流動(dòng)速度在翼型附近和尾跡剪切層變化劇烈，因此網(wǎng)格在翼型附近和尾跡區(qū)都進(jìn)行了加密。初始網(wǎng)格單元總數(shù)約為4萬(wàn)，由于流場(chǎng)模擬采用了三維網(wǎng)格，本文自適應(yīng)直接采用三維網(wǎng)格加密方法，因而自適應(yīng)后網(wǎng)格增量很大，約為970萬(wàn)。

圖11 30P30N三段翼幾何構(gòu)型和6個(gè)速度型站位

圖12 30P30N三段翼原始基準(zhǔn)網(wǎng)格

圖14給出了自適應(yīng)模擬過(guò)程中密度殘差收斂曲線，可以看出，自適應(yīng)后殘差收斂性更好。圖15給出了翼型壓力分布與試驗(yàn)值[38]及參考值[39]的比較，初始網(wǎng)格和自適應(yīng)后計(jì)算的壓力系數(shù)Cp分布與試驗(yàn)值吻合得都很好，且自適應(yīng)后主翼和襟翼上表面的壓力略低。

圖13 30P30N三段翼自適應(yīng)后的網(wǎng)格

圖14 網(wǎng)格自適應(yīng)迭代過(guò)程密度殘差收斂曲線

圖15 網(wǎng)格自適應(yīng)前后翼型壓力分布比較

圖16～圖18給出了圖11所示第1、第3和第5這3個(gè)站位的速率分布，其橫軸為采用自由來(lái)流值進(jìn)行無(wú)量綱化的當(dāng)?shù)厮俾蕌/ainf，縱軸為點(diǎn)到翼型表面的無(wú)量綱法向距離(參考值為平均氣動(dòng)弦長(zhǎng))。綜合來(lái)看，隨著自適應(yīng)迭代和網(wǎng)格在翼型附近加密，網(wǎng)格對(duì)流場(chǎng)的分辨率提高。與初始網(wǎng)格相比，自適應(yīng)后速度分布與試驗(yàn)值吻合度明顯提高，且本文計(jì)算結(jié)果比文獻(xiàn)[39]中采用高階格式計(jì)算的結(jié)果略好。

圖16 網(wǎng)格自適應(yīng)前后第1個(gè)站位處速度型分布比較

圖18 網(wǎng)格自適應(yīng)前后第5個(gè)站位處速度型分布比較

圖19給出了自適應(yīng)前后主翼前段區(qū)域的速率流場(chǎng)分布，前緣縫翼上表面產(chǎn)生的附面層低速氣流和流過(guò)縫隙的高速氣流在主翼前緣區(qū)域發(fā)生摻混。相比初始網(wǎng)格，自適應(yīng)后網(wǎng)格尺寸減小，分辨率增加，因此網(wǎng)格的數(shù)值耗散小，前緣縫翼上表面附面層發(fā)展更充分，附面層更厚，摻混的兩股氣流剪切作用更強(qiáng)，因此能夠更好地模擬出圖16所示的速度型中向左的“凹坑”。由于第1個(gè)站位離主翼前緣不遠(yuǎn)，因此主翼的附面層相對(duì)很薄，流動(dòng)速度在法向很小的距離內(nèi)迅速增加，對(duì)應(yīng)于圖16中第1個(gè)峰值。

對(duì)于襟翼上表面區(qū)域，由主翼產(chǎn)生的附面層氣流與流經(jīng)襟翼前部縫隙的高速氣流發(fā)生摻混，同時(shí)主翼上部的摻混區(qū)繼續(xù)向后發(fā)展，因此襟翼上表面的速度型出現(xiàn)了兩處向左的“凹坑”，如圖18所示。類(lèi)似地，通過(guò)網(wǎng)格自適應(yīng)對(duì)網(wǎng)格進(jìn)行加密后，速度型的模擬與試驗(yàn)吻合得更好。

圖19 網(wǎng)格自適應(yīng)前后主翼前段區(qū)域速度分布

5.2 三角翼大迎角繞流自適應(yīng)模擬

后掠角65°尖前緣三角翼(VFE-2)常用于驗(yàn)證和評(píng)估CFD方法預(yù)測(cè)復(fù)雜分離渦的能力，其幾何模型和試驗(yàn)測(cè)壓站位[40]如圖20所示。計(jì)算來(lái)流條件為：?jiǎn)挝焕字Z數(shù)Re=1.377×107，靜溫T∞=322.04 K，馬赫數(shù)Ma∞=0.4，側(cè)滑角β=0°，迎角α=20.17°，湍流模型為采用RC修正的一方程SA模型。

圖20 三角翼幾何構(gòu)型和試驗(yàn)測(cè)壓站位

圖21 三角翼流場(chǎng)模擬的原始基準(zhǔn)網(wǎng)格

圖22給出了自適應(yīng)迭代過(guò)程中密度殘差收斂曲線，由于初始網(wǎng)格不能有效捕捉空間分離流動(dòng)，密度殘差收斂較慢且共降低約2個(gè)量級(jí)。自適應(yīng)后密度殘差收斂性明顯改善，其收斂速度更快且下降更多。

圖23和圖24分別給出了x/c=0.6，0.95這2個(gè)站位等x截面上自適應(yīng)前后網(wǎng)格對(duì)比，自適應(yīng)后網(wǎng)格加密主要集中在主渦的渦核區(qū)域和旋渦脫離三角翼前緣后的強(qiáng)剪切區(qū)域。

圖22 三角翼網(wǎng)格自適應(yīng)過(guò)程密度殘差收斂曲線

圖23 自適應(yīng)前后x/c=0.6截面處網(wǎng)格對(duì)比

圖24 自適應(yīng)前后x/c=0.95截面處網(wǎng)格對(duì)比

圖25給出了三角翼x/c=0.6站位處上表面前緣附近自適應(yīng)前后由網(wǎng)格代表的幾何形狀對(duì)比，整體來(lái)看，自適應(yīng)加密并投影到擬合曲面后，前緣處的幾何形狀光順且連續(xù)性好。從局部放大圖可以看出，自適應(yīng)后新增網(wǎng)格點(diǎn)投影后的位置幾乎位于真實(shí)的幾何曲線上。由此表明，基于一定的網(wǎng)格密度，本文曲面擬合技術(shù)能高精度地重構(gòu)真實(shí)的幾何曲面。

自適應(yīng)前后分離渦附近的空間流線比較見(jiàn)圖26，初始網(wǎng)格計(jì)算的分離渦在接近三角翼后緣附近發(fā)生了破裂，二次渦在后緣附近與主渦相互混合；自適應(yīng)后主分離渦未發(fā)生破裂，主渦和二次渦向后運(yùn)動(dòng)的旋轉(zhuǎn)流態(tài)都保持得很好，且彼此獨(dú)立向后發(fā)展。

圖25 三角翼x/c=0.6截面上表面前緣處網(wǎng)格自適應(yīng)前后對(duì)比

圖26 三角翼網(wǎng)格自適應(yīng)前后空間流線對(duì)比

圖27給出了x/c=0.6站位處自適應(yīng)前后計(jì)算的壓力分布與試驗(yàn)值[40]的比較，自適應(yīng)后壓力分布與試驗(yàn)吻合更好。圖28給出了x/c=0.95站位處自適應(yīng)前后計(jì)算的壓力分布，由于自適應(yīng)后主渦未發(fā)生破裂，計(jì)算的壓力分布精度得到顯著提升。

圖27 自適應(yīng)前后x/c=0.6截面處壓力分布對(duì)比

圖28 自適應(yīng)前后x/c=0.95截面處壓力分布對(duì)比

6 結(jié) 論

本文針對(duì)非結(jié)構(gòu)混合網(wǎng)格，發(fā)展了網(wǎng)格單元分布優(yōu)化、網(wǎng)格幾何投影和空間網(wǎng)格匹配3項(xiàng)關(guān)鍵技術(shù)，建立了高魯棒性、幾何保真的網(wǎng)格自適應(yīng)系統(tǒng)，并采用30P30N三段翼型和三角翼兩個(gè)標(biāo)準(zhǔn)算例對(duì)自適應(yīng)技術(shù)進(jìn)行了驗(yàn)證。

1) 結(jié)合基于流場(chǎng)特征的自適應(yīng)探測(cè)器，基于網(wǎng)格加密和粗化方法建立的網(wǎng)格單元分布優(yōu)化技術(shù)能有效地實(shí)現(xiàn)計(jì)算網(wǎng)格的動(dòng)態(tài)優(yōu)化。

2) 基于僅依賴(lài)表面網(wǎng)格點(diǎn)信息的局部曲面擬合技術(shù)能高精度地重構(gòu)物面幾何信息，通過(guò)參數(shù)化目標(biāo)點(diǎn)可實(shí)現(xiàn)新增網(wǎng)格點(diǎn)的快速投影。

3) 基于改進(jìn)型的距離函數(shù)方法能有效解決因表面網(wǎng)格點(diǎn)投影造成的空間網(wǎng)格單元重疊交錯(cuò)的問(wèn)題，實(shí)現(xiàn)空間網(wǎng)格單元與投影后的表面網(wǎng)格協(xié)調(diào)匹配。

4) 本文建立的網(wǎng)格自適應(yīng)技術(shù)可以明顯改善流場(chǎng)求解的收斂性和提高流場(chǎng)細(xì)節(jié)的分辨能力，如壓力分布和速度分布等，從而提高飛行器氣動(dòng)特性的預(yù)測(cè)精度。

對(duì)于真實(shí)戰(zhàn)斗機(jī)或運(yùn)輸機(jī)等更復(fù)雜的幾何構(gòu)型，網(wǎng)格單元總量可達(dá)到千萬(wàn)量級(jí)。并行計(jì)算是減少自適應(yīng)系統(tǒng)運(yùn)行時(shí)間的有效方法，因而下一步將針對(duì)網(wǎng)格自適應(yīng)系統(tǒng)的并行化開(kāi)展研究。