劉田田， 董慶利， 楊 洋， 冷玨琳， 鄭澎

(1.中物院高性能數(shù)值模擬軟件中心，北京 100088； 2.北京應(yīng)用物理與計(jì)算數(shù)學(xué)研究所,北京 100094；3.中國(guó)工程物理研究院計(jì)算機(jī)應(yīng)用研究所,綿陽(yáng) 621900)

1 引 言

高雷諾數(shù)繞流仿真模擬中存在緊貼壁面的粘性力不可忽略的流動(dòng)薄層，即邊界層。在邊界層內(nèi)，物理量梯度變化很大，為了提升模擬精度，要求邊界層內(nèi)的網(wǎng)格盡量垂直于壁面。目前，分塊結(jié)構(gòu)網(wǎng)格和混合網(wǎng)格是粘性計(jì)算中常用的兩類網(wǎng)格。分塊結(jié)構(gòu)網(wǎng)格采用分而治之的策略，將模型分解成若干個(gè)子塊，每個(gè)塊內(nèi)生成結(jié)構(gòu)化網(wǎng)格。這類網(wǎng)格的正交性好，但對(duì)于復(fù)雜外形，需要花費(fèi)大量的時(shí)間進(jìn)行子塊分解，操作非常繁瑣?；旌暇W(wǎng)格在邊界層采用半結(jié)構(gòu)化的三棱柱單元，剩余區(qū)域采用四面體單元填充。半結(jié)構(gòu)化的邊界層網(wǎng)格在垂直于壁面方向具備結(jié)構(gòu)化特征，能很好捕捉邊界層的粘性流動(dòng)特征；在平行于壁面方向具備非結(jié)構(gòu)化特性，能靈活地適應(yīng)復(fù)雜幾何外形。這兩類網(wǎng)格的邊界層網(wǎng)格生成都是粘性模擬過(guò)程中的難點(diǎn)問(wèn)題，吸引了眾多學(xué)者的關(guān)注和研究。

分塊結(jié)構(gòu)網(wǎng)格的邊界層網(wǎng)格生成方法有兩類，一類是基于表面網(wǎng)格，利用雙曲方程或橢圓方程等生成邊界層網(wǎng)格，如商業(yè)軟件Pointwise等，這類方法無(wú)法自動(dòng)處理邊界層網(wǎng)格相交的情況，需人工一層一層生成。另一類是在拓?fù)浞謪^(qū)時(shí)預(yù)留出邊界層的子區(qū)域，利用超限映射等方法生成初始網(wǎng)格，再對(duì)網(wǎng)格進(jìn)行正交性優(yōu)化，如商業(yè)軟件ICEM等，這類方法面向復(fù)雜幾何外形時(shí)，增加了人工分區(qū)的難度和工作量。盧風(fēng)順等[1]提出一種適用于飛機(jī)的分塊結(jié)構(gòu)網(wǎng)格邊界層拓?fù)渥詣?dòng)生成方法，該方法對(duì)表面網(wǎng)格進(jìn)行幾何特征提取，根據(jù)幾何特征構(gòu)造邊界層網(wǎng)格框架，基于超限插值方法在邊界層網(wǎng)格框架內(nèi)生成邊界層網(wǎng)格。該算法在一定程度上提高了分塊網(wǎng)格生成的自動(dòng)化程度，減少了拓?fù)浞謪^(qū)的工作量。

經(jīng)過(guò)數(shù)十年的發(fā)展，適用于混合網(wǎng)格的邊界層生成算法研究已經(jīng)取得了很多成果[2-17]。然而，對(duì)于復(fù)雜幾何外形的邊界層網(wǎng)格生成，至今沒(méi)有形成完善的解決方案。目前，混合網(wǎng)格邊界層網(wǎng)格生成方法有如下三類。

(1) 局部前沿層進(jìn)方法。局部計(jì)算前沿點(diǎn)的法向和步長(zhǎng)，逐層生成邊界層網(wǎng)格。這類方法是目前研究最廣泛的一類方法，發(fā)展了可視錐[2]、最優(yōu)法向[7,18]、多法向[3]、法向合并[5]、步長(zhǎng)縮減[19,20]和自動(dòng)縮層[15,21]等多種技術(shù)處理局部細(xì)小特征和凹凸特征。

(2) 四面體聚合方法。在邊界層附近生成各向異性四面體網(wǎng)格，然后合并四面體得到三棱柱網(wǎng)格[16,22]。

(3) 全局前沿層進(jìn)方法。整體計(jì)算新前沿的位置。一種是將前沿的推進(jìn)過(guò)程轉(zhuǎn)化成物理量傳播的過(guò)程，通過(guò)求解PDE方程得到下一層前沿的位置[10]。另一種是利用水平集方法整體推進(jìn)前沿[9,23]。

目前混合網(wǎng)格的邊界層網(wǎng)格大都是利用三角形生成三棱柱單元，實(shí)際上這類方法同樣適用于四邊形生成六面體單元。分塊結(jié)構(gòu)網(wǎng)格的表面網(wǎng)格拓?fù)浞謪^(qū)可以看成是一個(gè)非結(jié)構(gòu)四邊形網(wǎng)格，因此，混合網(wǎng)格邊界層的網(wǎng)格生成方法可以用于分塊結(jié)構(gòu)網(wǎng)格邊界層拓?fù)浞謪^(qū)的生成。分塊結(jié)構(gòu)網(wǎng)格的一類邊界層生成方法是在拓?fù)浞謪^(qū)內(nèi)進(jìn)行結(jié)構(gòu)網(wǎng)格細(xì)分，這種思想也可以用于混合網(wǎng)格邊界層的生成。因此，混合網(wǎng)格的邊界層網(wǎng)格生成方法和分塊結(jié)構(gòu)網(wǎng)格的邊界層網(wǎng)格生成方法可以互相借鑒，形成優(yōu)勢(shì)互補(bǔ)。受文獻(xiàn)[1，20]方法的啟發(fā)，在已有研究成果的基礎(chǔ)上，本文提出了一種雙前沿推進(jìn)思想，綜合采用多種技術(shù)處理局部幾何特征，形成了一種邊界層網(wǎng)格生成算法。

2 算法概述

雙前沿推進(jìn)過(guò)程中，存在標(biāo)架前沿和參考前沿兩套數(shù)據(jù)。標(biāo)架前沿是需要生成邊界層的前沿，參考前沿可以是更密的網(wǎng)格也可以是幾何，用來(lái)輔助標(biāo)架前沿在推進(jìn)過(guò)程中進(jìn)行相交檢測(cè)。如圖1所示，在分塊結(jié)構(gòu)網(wǎng)格中，標(biāo)架前沿是表面的拓?fù)浞謪^(qū)，參考前沿是表面網(wǎng)格或幾何數(shù)據(jù)。在混合網(wǎng)格中，標(biāo)架前沿是待生成邊界層的表面網(wǎng)格，參考前沿可以是加密的表面網(wǎng)格，如果參考前沿和標(biāo)架前沿使用同一套網(wǎng)格，則等同于以往的混合網(wǎng)格邊界層前沿推進(jìn)算法。因此，雙前沿推進(jìn)算法可以看成是局部前沿層進(jìn)算法的擴(kuò)展。

基于雙前沿推進(jìn)的邊界層網(wǎng)格生成算法流程如圖2所示，算法的輸入、輸出和具體流程如下。

圖2 雙前沿推進(jìn)邊界層網(wǎng)格生成流程

輸入：標(biāo)架前沿、參考前沿、邊界層的第一層高度、相鄰層之間的過(guò)渡比以及總的邊界層層數(shù)等。

輸出：邊界層網(wǎng)格。

流程： (1) 構(gòu)建標(biāo)架前沿和參考前沿；(2) 計(jì)算標(biāo)架前沿點(diǎn)的法向；(3) 計(jì)算標(biāo)架前沿點(diǎn)的步長(zhǎng)；(4) 光滑化標(biāo)架前沿點(diǎn)的法向和步長(zhǎng)；(5) 判斷標(biāo)架前沿中是否存在相交，若存在相交，調(diào)整標(biāo)架前沿點(diǎn)的法向或步長(zhǎng)，若不存在，執(zhí)行下一步；(6) 插值計(jì)算參考前沿點(diǎn)的法向和步長(zhǎng)；(7) 判斷參考前沿中是否存在相交，若存在相交，調(diào)整標(biāo)架前沿的法向或步長(zhǎng)，執(zhí)行(6)，若不存在，執(zhí)行下一步；(8) 生成一層邊界層單元，判斷邊界層單元是否有效，如網(wǎng)格單元質(zhì)量、單元是否反轉(zhuǎn)等；(9) 根據(jù)需要細(xì)化邊界層單元；(10) 優(yōu)化邊界層網(wǎng)格單元質(zhì)量。

與傳統(tǒng)的局部前沿層進(jìn)方法相比，雙前沿推進(jìn)算法不僅考慮當(dāng)前待推進(jìn)前沿的特征，還考慮了更精細(xì)的幾何特征，其優(yōu)勢(shì)主要包括兩個(gè)，(1) 將局部前沿層進(jìn)算法從混合網(wǎng)格推廣至分塊結(jié)構(gòu)網(wǎng)格，可同時(shí)適用于分塊結(jié)構(gòu)網(wǎng)格的邊界層生成和混合網(wǎng)格的邊界層生成。將結(jié)構(gòu)網(wǎng)格的表面分區(qū)看成是非結(jié)構(gòu)網(wǎng)格，利用雙前沿推進(jìn)算法可以有效提升生成結(jié)構(gòu)網(wǎng)格分區(qū)拓?fù)渲羞吔鐚臃謪^(qū)構(gòu)建的自動(dòng)化程度。(2) 由于考慮了更精細(xì)的特征，能有效避免因網(wǎng)格離散精度不夠?qū)е碌倪吔鐚游恢貌缓侠怼?/p>

3 邊界層網(wǎng)格生成

3.1 標(biāo)架前沿點(diǎn)的法向計(jì)算

最傳統(tǒng)的方法采用相鄰單元的法向平均作為前沿點(diǎn)的生長(zhǎng)方向，但在尖銳特征附近極易產(chǎn)生反轉(zhuǎn)單元。文獻(xiàn)[2]提到，不產(chǎn)生反轉(zhuǎn)單元的一個(gè)必要條件是前沿點(diǎn)的生長(zhǎng)法向與點(diǎn)連接的所有單元都是可見(jiàn)的，即法向在可視錐范圍內(nèi)，如圖3所示。設(shè)點(diǎn)P周圍的面單元集合為Fp，面單元的法向分別為Nj(j∈Fp)。三角形單元的法向取其外法向方向，對(duì)于四邊形單元，首先將其劃分成兩個(gè)三角形，單元的法向取兩個(gè)三角形法向的平均值。對(duì)每個(gè)前沿點(diǎn)，可視錐定義為[2]

{M∈R3|PM·Nj>0,?j∈Fp}

(1)

圖3 可視錐

可視錐本身是一個(gè)多面錐，為了簡(jiǎn)化，可以用圓錐來(lái)替代。文獻(xiàn)[7]的最優(yōu)法向技術(shù)可以保證計(jì)算出的法向在可視錐范圍內(nèi)。最優(yōu)法向的計(jì)算方法是，極小化法向與周圍單元法向的最大夾角。該問(wèn)題的求解可以轉(zhuǎn)化為點(diǎn)集的最小覆蓋圓問(wèn)題，其中兩點(diǎn)之間的距離定義為球面距離。有定理表明，最小覆蓋圓一定會(huì)穿過(guò)點(diǎn)集的兩個(gè)或至多三個(gè)點(diǎn)。由于每個(gè)前沿點(diǎn)周圍的面?zhèn)€數(shù)一般在4～8個(gè)左右，可以采用枚舉的方法求解該問(wèn)題。計(jì)算所有過(guò)任意兩點(diǎn)或三點(diǎn)的覆蓋圓，篩選其中半徑最小的圓，該圓心與P點(diǎn)的連線方向就是前沿點(diǎn)的法向。

圖4(a)給出了利用最優(yōu)法向計(jì)算得到的前沿點(diǎn)的初始方向。這種方法對(duì)一般的幾何外形可以滿足，但是對(duì)于有凹角的外形，在凹角處很容易網(wǎng)格相交。為了減少凹角處的網(wǎng)格相交，采用加權(quán)的Laplacian平滑算法對(duì)法向進(jìn)行平滑，

(2)

圖4 標(biāo)架前沿點(diǎn)法向

表1 平滑前后凹角點(diǎn)的最大步長(zhǎng)

3.2 標(biāo)架前沿點(diǎn)的總步長(zhǎng)計(jì)算

邊界層的第一層高度記為h1，相鄰層之間的過(guò)渡比為t，總的邊界層層數(shù)為n，則計(jì)算出邊界層的初始總高度為

(3)

為減少前沿相交，首先利用射線法計(jì)算每個(gè)前沿點(diǎn)的最大可行步長(zhǎng)。為節(jié)省計(jì)算量，采用R-Tree對(duì)前沿單元進(jìn)行空間位置分割?；赗-Tree和單元的碰撞檢測(cè)技術(shù)，計(jì)算從前沿點(diǎn)Pi出發(fā)沿法向方向Ni長(zhǎng)度為2h的線段是否與其他前沿相交，若相交，則計(jì)算第一個(gè)交點(diǎn)Qi與Pi之間的距離，記為前沿點(diǎn)Pi的最大可行步長(zhǎng)，hmax(i)=|PiQi|，否則，hmax(i)=2h。根據(jù)最大可行步長(zhǎng)，重新計(jì)算每個(gè)前沿點(diǎn)的步長(zhǎng)為

(4)

前沿面以波陣方式向前推進(jìn)時(shí)更有利于凹凸特征附近生成光滑的前沿面，從而提高網(wǎng)格單元質(zhì)量。因此，在凹角點(diǎn)附近減小步長(zhǎng)，凸角點(diǎn)附近增加步長(zhǎng)，像吹氣球一樣使步長(zhǎng)平滑過(guò)渡。凹凸角點(diǎn)的判定方法參見(jiàn)文獻(xiàn)[15]。根據(jù)凹凸性計(jì)算步長(zhǎng)的公式為

(5)

式中βi為前沿點(diǎn)的特征角，取Pi鄰接的所有邊二面角的平均值，

(6)

dihedral(lj)為Pi鄰接的第j條邊相鄰兩個(gè)面的二面角。

圖5 特征角

圖6給出了一個(gè)立方體分別向內(nèi)和向外生長(zhǎng)邊界層時(shí)的步長(zhǎng)。向外生長(zhǎng)時(shí)立方體邊上的點(diǎn)是凸點(diǎn)，圖中步長(zhǎng)明顯縮減，向內(nèi)生長(zhǎng)時(shí)立方體邊上的點(diǎn)是凹點(diǎn)，圖中步長(zhǎng)明顯增大。為了避免相鄰點(diǎn)之間的步長(zhǎng)相差過(guò)大導(dǎo)致單元畸形，對(duì)前沿點(diǎn)的步長(zhǎng)進(jìn)行Laplacian平滑，

(7)

圖6 凹凸特征步長(zhǎng)

3.3 標(biāo)架前沿相交性檢測(cè)

為了避免標(biāo)架前沿生成網(wǎng)格時(shí)產(chǎn)生相交問(wèn)題，需檢測(cè)前沿推進(jìn)后在空間中是否相交，并據(jù)此調(diào)整標(biāo)架前沿點(diǎn)的步長(zhǎng)。采用R-Tree空間搜索樹(shù)對(duì)數(shù)據(jù)進(jìn)行空間分割，從而提升數(shù)據(jù)搜索效率。具體算法流程如下。

(1) 將標(biāo)架前沿的原始前沿面和推進(jìn)后的前沿面進(jìn)行三角化。

(2) 構(gòu)建所有三角形單元的空間包圍盒搜索樹(shù)。

(3) 對(duì)每一個(gè)推進(jìn)后的前沿面，利用搜索樹(shù)查找與其包圍盒相交的單元，采用三角形碰撞檢測(cè)判斷是否存在實(shí)際相交。若相交，則縮短其包含的前沿點(diǎn)的步長(zhǎng)，直到不再相交,h″(i)=ε·h′(i)，其中ε為步長(zhǎng)收縮因子，取值[0.5,1.0)。對(duì)于無(wú)法通過(guò)步長(zhǎng)縮減避免相交的前沿單元，局部調(diào)整其包含的前沿點(diǎn)的法向，再進(jìn)行相交檢測(cè)，直到不再相交。

3.4 參考前沿相交性檢測(cè)

將計(jì)算參考前沿的每一個(gè)前沿點(diǎn)投影到標(biāo)架前沿，計(jì)算投影點(diǎn)相鄰的標(biāo)架前沿點(diǎn)，根據(jù)這些點(diǎn)的法向和步長(zhǎng)插值計(jì)算參考前沿點(diǎn)的法向和步長(zhǎng)。為了避免復(fù)雜特征處，由于網(wǎng)格離散精度不夠?qū)е碌倪吔鐚泳W(wǎng)格不合理的問(wèn)題，判斷參考前沿點(diǎn)推進(jìn)后是否會(huì)存在空間相交。相交性檢測(cè)的方法同 3.3 節(jié)相同，如果參考前沿單元相交，需縮短標(biāo)架前沿點(diǎn)的步長(zhǎng)，并重新計(jì)算參考前沿點(diǎn)的步長(zhǎng)，直到參考前沿中不存在相交。

3.5 單元有效性檢測(cè)和修正

確定前沿點(diǎn)的法向和步長(zhǎng)，可以獲得一層總的邊界層單元。為了避免產(chǎn)生質(zhì)量非常差的單元，判斷這一層總的邊界層單元中是否存在無(wú)效或質(zhì)量差的單元。判斷網(wǎng)格質(zhì)量的標(biāo)準(zhǔn)是計(jì)算單元的Jacobian值和扭曲度。Jacobian值若為負(fù)值，表示單元反轉(zhuǎn)。扭曲度在一定程度上能表征與理想單元的偏離程度，定義如下，

(8)

式中θe為面單元或體單元的等角，θmin和θmax為面單元或體單元的最小角和最大角。對(duì)于體單元，計(jì)算偏斜度時(shí)，采用skewness=max{Surface(skew),Volume(skew)}，其中Surface(skew)為體單元包含的所有面單元的偏斜度，Volume(skew)為利用體單元二面角計(jì)算出的偏斜度。

對(duì)于Jacobian異常或扭曲度差的單元，縮短其包含的前沿點(diǎn)的步長(zhǎng)或局部調(diào)整法向，直到所有單元質(zhì)量達(dá)到要求。

3.6 邊界層單元細(xì)化

對(duì)于分塊結(jié)構(gòu)網(wǎng)格，標(biāo)架前沿推進(jìn)后得到的每一個(gè)單元都是邊界層的一個(gè)拓?fù)渥訅K。在每個(gè)拓?fù)渥訅K內(nèi)，根據(jù)邊界層參數(shù)要求，采用超限映射方法(TFI)對(duì)邊界層網(wǎng)格進(jìn)行細(xì)化，再將網(wǎng)格點(diǎn)投影到真實(shí)的幾何。對(duì)于混合網(wǎng)格，得到一層總的邊界層網(wǎng)格后，根據(jù)邊界層的層數(shù)和過(guò)渡比，計(jì)算每個(gè)前沿點(diǎn)對(duì)應(yīng)的每層邊界層單元的高度，將總的邊界層細(xì)分成指定層的邊界層網(wǎng)格。記前沿點(diǎn)Pi的總步長(zhǎng)為h*(i)，則Pi對(duì)應(yīng)的每層邊界層單元高度為

(9)

3.7 邊界層網(wǎng)格優(yōu)化

分塊結(jié)構(gòu)網(wǎng)格的邊界層網(wǎng)格可以采用帶源項(xiàng)的橢圓方程進(jìn)行網(wǎng)格正交性優(yōu)化?；旌暇W(wǎng)格的邊界層網(wǎng)格采用文獻(xiàn)[9]的能量函數(shù)法改善邊界層單元的質(zhì)量。三棱柱的能量函數(shù)定義可以推廣到六面體網(wǎng)格，利用采用棱邊單元的正交性能量函數(shù)E⊥(τ)和底面單元的形態(tài)能量函數(shù)Eθ(τ)，定義單元的質(zhì)量能量如下，

E(τ)=μEθ(τ)+(1-μ)E⊥(τ)

(10)

式中μ為底面單元形態(tài)與棱邊正交性之間的權(quán)重比重。優(yōu)化網(wǎng)格質(zhì)量的方法是極小化所有單元的總能量，

(11)

利用擬牛頓迭代求解該極小化問(wèn)題。

4 數(shù)值實(shí)驗(yàn)

輸入前沿后，給定邊界層第一層的高度、層高的過(guò)渡比和總的層數(shù)，可以自動(dòng)生成邊界層網(wǎng)格。以下實(shí)例中，混合網(wǎng)格的表面網(wǎng)格為三角形或四邊形單元，遠(yuǎn)場(chǎng)網(wǎng)格采用四面體單元。分塊結(jié)構(gòu)網(wǎng)格的邊界層內(nèi)采用TFI生成初始網(wǎng)格，再用橢圓優(yōu)化進(jìn)一步優(yōu)化網(wǎng)格質(zhì)量。圖7給出一個(gè)簡(jiǎn)單球體的網(wǎng)格生成效果，包含了混合網(wǎng)格和多塊結(jié)構(gòu)網(wǎng)格的邊界層網(wǎng)格，驗(yàn)證了雙前沿推進(jìn)算法對(duì)兩類網(wǎng)格的適用性。

圖7 球體模型邊界層網(wǎng)格

圖8給出了三維翼型模型的混合網(wǎng)格和多塊結(jié)構(gòu)網(wǎng)格結(jié)果。翼型區(qū)域內(nèi)表面設(shè)置為壁面。從圖8(a)可以看出混合網(wǎng)格的邊界層過(guò)渡非常平滑，單元形態(tài)規(guī)整。圖8(b)左圖是翼型模型的拓?fù)浞謪^(qū)，將該拓?fù)浞謪^(qū)映射到幾何模型后生成多塊結(jié)構(gòu)網(wǎng)格，如圖8(b)右圖所示。

圖8 翼型模型

為了進(jìn)一步驗(yàn)證算法的有效性，對(duì)F6模型進(jìn)行了混合網(wǎng)格剖分，如圖9所示。圖10對(duì)邊界層網(wǎng)格的扭曲度進(jìn)行了統(tǒng)計(jì)?？梢钥闯鼋^大部分邊界層網(wǎng)格質(zhì)量都很好，但包含非常少量的扭曲度大于0.9的單元，這也是后續(xù)要繼續(xù)研究的問(wèn)題。

圖9 F6模型混合網(wǎng)格

圖10 F6混合網(wǎng)格質(zhì)量統(tǒng)計(jì)

5 結(jié) 論

本文面向粘性流體模擬的邊界層網(wǎng)格生成問(wèn)題，提出了一種雙前沿推進(jìn)思想，并發(fā)展了一種同時(shí)適用于分塊結(jié)構(gòu)網(wǎng)格邊界層生成和混合網(wǎng)格邊界層生成的算法。為了適應(yīng)復(fù)雜幾何特征，算法綜合采用了最優(yōu)法向技術(shù)、步長(zhǎng)調(diào)整技術(shù)和極小化能量?jī)?yōu)化網(wǎng)格等關(guān)鍵技術(shù)。測(cè)試算例表明了該算法的有效性，但對(duì)復(fù)雜模型在凹凸特征附近會(huì)產(chǎn)生少量質(zhì)量差的單元，這是下一步研究工作的重點(diǎn)。