非結(jié)構(gòu)混合網(wǎng)格魯棒自適應(yīng)技術(shù)

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

1. 中國空氣動力研究與發(fā)展中心 計算空氣動力研究所， 綿陽 621000

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

(1)

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

(2)

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

(3)

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

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

gr=gmaxtr,gc=gmaxtc

(4)

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

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

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

(5)

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

(6)

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

(7)

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

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

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

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

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

1) 三角形剖分為4個三角形，見圖1(a)。

2) 四邊形剖分為4個四邊形，見圖1(b)。

3) 四面體剖分為8個四面體，見圖1(c)。

4) 三棱柱剖分為8個三棱柱，見圖1(d)。

5) 六面體剖分為8個六面體，見圖1(e)。

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

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

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

相鄰兩個單元層級不同時，在兩個單元的交界面或交界線上將出現(xiàn)圖2所示的懸空節(jié)點。由于常規(guī)流場模擬軟件不支持懸空節(jié)點，為了提高網(wǎng)格自適應(yīng)方法對流場模擬軟件的通用性，采用多面體轉(zhuǎn)換方法間接消除懸空節(jié)點。以圖2左側(cè)單元為例，轉(zhuǎn)化為多面體后共包括9個面：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個面為標準的四邊形，后4個面為包含5個頂點的多邊形。一般情況下，基于非結(jié)構(gòu)網(wǎng)格的流場模擬軟件只支持三角形和四邊形兩種標準網(wǎng)格面單元。因此，為了進一步提高自適應(yīng)系統(tǒng)的通用性，本文將多邊形面轉(zhuǎn)換成多個三角形和四邊形。對于三角形面，標準化轉(zhuǎn)換共有3種拓撲類型，如圖3所示；對于四邊形面，共有5種拓撲類型，如圖4所示。

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

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

圖3 非標準三角形網(wǎng)格面拓撲類型

圖4 非標準四邊形網(wǎng)格面拓撲類型

圖5 多面體單元多邊形網(wǎng)格面標準化示意圖

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

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

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

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

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

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

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

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

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

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

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

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

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

(8)

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

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

本文采用三次Hermit插值直接構(gòu)造三維邊界曲線。如圖6所示，表面網(wǎng)格邊兩個頂點的坐標為x1和x2，相應(yīng)法向為n1和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)

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

對于三角形網(wǎng)格面，如圖7所示，3個網(wǎng)格點為x1～x3，相應(yīng)法向為n1～n3。首先將三角形x1x2x3轉(zhuǎ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)格邊，實線：重構(gòu)的邊界曲線)

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

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

對于四邊形網(wǎng)格面，如圖8所示，4個網(wǎng)格點為x1～x4，相應(yīng)法向為n1～n4。類似地，首先將四邊形x1x2x3x4轉(zhuǎ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)格點投影

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

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

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

(13)

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

(14)

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

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

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

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

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

(15)

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

6 結(jié) 論

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

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

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

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

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

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