非結構混合網格魯棒自適應技術

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

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

2. 西北工業大學 航空學院，西安 710072

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

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

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

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

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

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

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

懸空節點類方法對網格單元剖分方向不作任何限制，但要求流場解算器能處理懸空節點。該類方法可分為基于網格面和基于網格單元兩類。基于網格面的方法[22-23]可在一定程度上實現網格各向異性分布，但通常網格單元的幾何質量隨著自適應迭代下降很快。基于網格單元的方法既可用于各向異性的網格自適應[24]，也可用于各向同性的網格自適應[25]，后者網格單元剖分類型很少，對數據結構的依賴較輕，軟件實現更容易。

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

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

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

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

1 網格單元加密或稀疏的自適應判據

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

1.1 擬梯度自適應判據

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

(1)

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

(2)

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

(3)

式中:NF為面相鄰網格單元的總數。

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

gr=gmaxtr,gc=gmaxtc

(4)

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

1.2 剪切應變率自適應判據

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

(5)

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

(6)

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

(7)

式中：ε為過濾閾值；N為網格單元總數；h取值0.1～0.5。

2 網格單元分布優化方法

本文采用剖分或合并網格單元的方法實現網格分布的優化。為提高網格加密或粗化操作的運行效率，應建立合適的數據結構。

2.1 網格細化方法

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

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 非結構混合網格單元加密模式

定義1自適應層級某網格單元的自適應層級是指由初始網格單元得到該網格單元需要加密的次數。

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

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

圖2 懸空節點示意圖

圖3 非標準三角形網格面拓撲類型

圖4 非標準四邊形網格面拓撲類型

圖5 多面體單元多邊形網格面標準化示意圖

2.2 網格粗化方法

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

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

準則1自適應加密/粗化準則某網格單元自適應加密或粗化屬性的標記滿足如下3條規則：

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

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

3) 加密優先規則。當某個網格單元標記為粗化卻因規則1)或規則2)需要標記為加密時，加密標記優先于粗化。

2.3 網格自適應的數據結構

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

3 表面網格幾何投影方法

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

3.1 物面幾何重構

本文物面幾何重構方法采用局部Coons曲面擬合，該方法只依賴于表面網格點坐標和相應的法向矢量。法向單位矢量通常不是標準的網格數據，本文通過與網格點直接相連的面網格的法向加權得到，即網格點i的法向向量ni為

(8)

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

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

本文采用三次Hermit插值直接構造三維邊界曲線。如圖6所示，表面網格邊兩個頂點的坐標為x1和x2，相應法向為n1和n2，參數化的邊界曲線c(t)為

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

(9)

式中：f和g為權值；τ為相應切向修正量，即

d1=x2-x1,d2=x1-x2

(10)

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

對于三角形網格面，如圖7所示，3個網格點為x1～x3，相應法向為n1～n3。首先將三角形x1x2x3轉化到標準三角形參數平面V1V2V3，V1V2方向參數為S′，V1V3方向參數為S″。網格邊x1x2構成邊界曲線c1，網格邊x2x3構成邊界曲線c2，網格邊x1x3構成邊界曲線c3，則三維擬合曲面S(S′,S″)表示為

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

(11)

圖6 三維邊界曲線重構示意圖(虛線：網格邊，實線：重構的邊界曲線)

圖7 三角形網格曲面重構示意圖

圖8 四邊形網格曲面重構示意圖

對于四邊形網格面，如圖8所示，4個網格點為x1～x4，相應法向為n1～n4。類似地，首先將四邊形x1x2x3x4轉化到標準四邊形參數平面V1V2V3V4，V1V2方向參數為S′，V1V4方向參數為S″。網格邊x1x2構成邊界曲線c1，網格邊x2x3構成邊界曲線c2，網格邊x4x3構成邊界曲線c3，網格邊x1x4構成邊界曲線c4，則三維的擬合曲面S(S′,S″)表示為

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

(12)

3.2 新增網格點投影

本文擬合曲面為參數曲面，可以采用參數點映射方法實現新增表面網格點的幾何投影，即先在標準參數面上取投影點的參數坐標，然后代入參數曲面方程求得投影點在曲面上的坐標點，最后將新增網格點移動到目標點位置。

如圖9所示，對于三角形網格面，常用的新增表面網格點的位置為N1、N2和N3，對應參數坐標為

圖9 網格參數曲面和新增點參數

(13)

對于四邊形網格面，常用的新增點位置為N1、N2、N3、N4和N5，對應參數坐標為

(14)

4 空間網格匹配方法

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

圖10 表面網格投影造成網格單元重疊交錯示意圖

本文采用改進型的基于距離函數的動網格技術實現表面網格投影后體網格單元的快速匹配。該方法通過引入自調節的影響半徑大大提高了附面層內半結構化網格的變形能力，同時通過包圍盒加速技術提高了變形方法的計算效率。

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

(15)

式中：NS為表面網格點總數；N為體網格點總數；η為衰減函數；ωi為表面網格點i的權值。衰減函數和權值的取值方法以及更多細節信息可參考文獻[31]。

5 網格自適應方法數值驗證

5.1 30P30N三段翼構型繞流自適應模擬

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

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

圖11 30P30N三段翼幾何構型和6個速度型站位

圖12 30P30N三段翼原始基準網格

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

圖13 30P30N三段翼自適應后的網格

圖14 網格自適應迭代過程密度殘差收斂曲線

圖15 網格自適應前后翼型壓力分布比較

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

圖16 網格自適應前后第1個站位處速度型分布比較

圖18 網格自適應前后第5個站位處速度型分布比較

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

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

圖19 網格自適應前后主翼前段區域速度分布

5.2 三角翼大迎角繞流自適應模擬

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

圖20 三角翼幾何構型和試驗測壓站位

圖21 三角翼流場模擬的原始基準網格

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

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

圖22 三角翼網格自適應過程密度殘差收斂曲線

圖23 自適應前后x/c=0.6截面處網格對比

圖24 自適應前后x/c=0.95截面處網格對比

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

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

圖25 三角翼x/c=0.6截面上表面前緣處網格自適應前后對比

圖26 三角翼網格自適應前后空間流線對比

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

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

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

6 結 論

本文針對非結構混合網格，發展了網格單元分布優化、網格幾何投影和空間網格匹配3項關鍵技術，建立了高魯棒性、幾何保真的網格自適應系統，并采用30P30N三段翼型和三角翼兩個標準算例對自適應技術進行了驗證。

1) 結合基于流場特征的自適應探測器，基于網格加密和粗化方法建立的網格單元分布優化技術能有效地實現計算網格的動態優化。

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

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

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

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