黏性邊界層網格自動生成1)

甘洋科 劉劍飛

(北京大學工學院力學與工程科學系，北京100871)

黏性邊界層網格自動生成1)

甘洋科 劉劍飛2)

(北京大學工學院力學與工程科學系，北京100871)

高雷諾數黏性流動在壁面附近存在邊界層，在計算模擬中自動生成可靠且有效的計算網格仍然是計算流體力學存在的瓶頸.三棱柱/四面體混合網格技術在一定程度上緩解了這個困難.然而，對于復雜外形的情況，在邊界層內自動高效生成高質量的三棱柱單元仍然十分困難.常用的層推進法在凹凸區域及角點處生成的邊界層網格單元質量較差，邊界層網格最外層尺寸不均勻.針對這些問題，發展了一種黏性邊界層網格自動生成方法，通過頂點周圍邊的二面角識別物面網格特征確定多生長方向，預估并調整生長高度處理相交情況.同時提出一種三維前沿尺寸調節方式，提高了邊界層網格單元的正交性，保證了邊界層網格與遠場網格尺寸的光滑過渡.通過ONERA M6翼型以及帶發動機短艙的DLR-F6翼身組合體等外形的網格生成實例及繞流數值模擬，將計算值與標準實驗值進行對比，結果表明：該方法能夠自動高效地生成滿足數值計算需求的混合網格.

邊界層網格，生長法向，特征邊，相交，尺寸調節

引言

近年來，計算流體力學成為學術界研究的熱點，在航空航天、汽車、熱能動力等工程領域都有著重要的地位，廣泛用于設計、分析和優化.在計算流體力學模擬中，需要借助計算機數值求解物理量的控制方程.當前主要采用的計算方法，包括有限差分法、有限體積法、有限元法等，都需要對計算區域進行網格劃分，進而在網格單元節點上離散求解偏微分方程.30多年來，盡管計算流體力學取得了很大的成就，網格生成依然是需要解決的關鍵問題之一.大批商業網格生成軟件，包括Gambit，Gridgen以及ANSYS ICEM-CFD等，雖然具有功能全面、使用方便、技術服務好等優點，但在生成含邊界層網格時的自動性及可靠性較差，對于復雜外形需要大量手工操作，有時甚至無法生成.高雷諾數黏性流動中不能自動生成可靠且有效的網格仍然是計算流體力學中的瓶頸問題.

邊界層是流體運動中貼近壁面的薄層，其厚度與物體特征尺寸的比值為一個小量，通常小于千分之一.如果采用各向同性網格，則生成的邊界層網格單元數量巨大，大大提高了計算的代價，對于大規模的復雜問題，其計算代價是目前工程上不能接受的.而如果計算域都采用大長寬比的各向異性單元則會得到質量很差的數值解[1].因此，邊界層網格生成的難點在于如何用盡可能少的網格單元達到捕捉層內流動特性的目的.

三維黏性混合網格策略，即在邊界層內生成高質量的半結構化三棱柱單元而在其他遠場區域生成各向同性的四面體網格單元，是目前最常用的方法，這種方法兼具結構網格和非結構網格的優點.Pirzadeh[2]提出的層推進法(advancing layer method)是其中最具影響的一種.層推進方法從需要生長邊界層的表面三角網格出發，沿通過表面頂點的某一個方向布置各向異性網格的節點，然后根據表面網格的拓撲關系連接這些布置的節點來生成一層一層的三棱柱單元，有時由于求解器的原因，三棱柱單元可以轉化為四面體，如圖1所示.

圖1 三棱柱單元生長示意圖Fig.1 Prism elements growing example

圍繞層推進法，很多學者做了不同程度的改進，Kallingderis等[3]改進了推進向量的求法，提出了可見區域的概念，利用該方法求得的向量可以保證生成的單元都有正的體積.L¨ohner[4]提出了連接模型邊界上生成的半結構化網格與其他區域生成的各向同性網格的方法，在半結構化邊界層網格生成過程中，檢測形狀差,反轉、相交、尺寸不適合的單元，將相應的單元刪除.Aubry等[5]提出了一種求節點法向的方法，并將其用于表面各向異性網格的生成[67]，最近又利用Voronoi圖來確定非流形角點處的生長法向，進而更好地生成邊界層網格，并集成到軟件中[8].Connell等[9]實現了全四面體黏性網格的生成，在確保相容性的前提下將三棱柱單元全部轉化為四面體單元，通過刪除單元的方式處理相交問題.Garimella等[10]引進多生長方向，提高了非流形角點處的單元正交性及單元質量，提出采用收縮、平滑、裁剪等操作處理層推進過程中的單元自交及前沿相交情況.Athanasiadis等[1112]則提出了構造虛擬點來實現多法向生長及法向融合操作，提高了凹凸區域黏性網格的正交性及結構化特性.Sharov等[13]通過生成各向異性的表面網格來適應表面凹凸區域的半結構化黏性網格生成.王剛等[1415]在這方面做了一些工作，將層推進方法應用到黏性網格生成.陳建軍等[16]也在已有的研究基礎上，集成各種處理方法，發展了一套混合網格生成程序.

以上所述的方法，在出現相交情況的時候采用刪除單元或停止推進的方式進行處理，這樣會導致各向同性網格生成器的初始表面網格尺寸不均勻，存在長寬比很大的狹長單元，不利于各向同性網格的生成;同時邊界層網格的層狀結構也會破壞.另外，該方法在凹凸脊及角點等幾何特征處需要特殊處理，現有的處理方法都會產生一定的問題，尚未合理地解決通用性的問題.

除了層推進方法，一些學者也提出了其他方法.Dyedov等[17]將偏移面(face o ff setting)方法[18]應用到黏性混合網格生成中，考慮面追蹤，并直接推進表面三角形單元，事先確定推進距離的上限，避免相交檢測，并提出采用變分優化的方法對三棱柱網格進行質量改善，該方法通過求解拉格朗日演化方程的方式來生成三棱柱邊界層網格，并將該混合網格生成方法應用到生物力學計算中，取得了一定效果.

與層推進方法不同，Ito等[19]以及Alauzet等[20]發展了一種黏性混合網格的生成方法，先生成全域的非結構各向同性網格，然后采用擠壓內部網格的方式生長邊界層半結構化網格.前者不改變各向同性網格的拓撲結構，后者結合拓撲可變網格變形方法.該方法可以避免網格生長過程中的相交檢測，但受限于網格變形技術，由于網格變形能力有限，在凹角點區域可以出現單元反轉的現象，同時由于擠壓，內部非結構網格的質量下降，影響計算的精度.

張來平等[2123]在生成壁面三棱柱網格方面發展了一種聚合的方法，將物面附近各向異性的四面體單元聚合生成三棱柱單元，從而得到半結構化的邊界層網格，并將生成的混合網格應用到非定常數值計算中[24-25].

本文在前人研究成果的基礎上，基于前沿層推進方法，發展了一套黏性邊界層網格自動生成程序，其中的細節處理方法，包括前沿多法向的確定以及生長高度的調節等，很好地解決了凹凸區域及角點附近的邊界層網格生成穩定性問題，能夠自動高效地生成可靠且有效的邊界層網格.同時發展了一種三維尺寸調節方式，可以為各向同性網格生成器提供均勻尺寸的初始表面網格，進而得到高質量的遠場各向同性網格，有利于提高流體計算解的質量.

1 邊界層網格生成流程

本文算法步驟如下:

輸入：計算區域物面三角形網格，邊界層首層厚度hf，厚度增長比例hr，層數nl.

輸出：邊界層網格，及最外層三角形網格前沿.

Step1讀入物面三角形網格數據，建立網格數據結構及拓撲關系.

Step2調整各節點總生長高度，此時生長高度為用戶給定總生長高度

Step 2.1提取前沿特征信息，包括特征邊及角點等信息，確定節點生長線數量及方向.

Step 2.2根據特征信息生長對應的單元，采用收縮高度的方式處理單元自交(反轉)問題，得到各點的合法生長高度.

Step 2.3對新前沿進行全局相交檢測，收縮高度避免全局相交，得到各點新的生長高度.采用拉普拉斯方法進行高度光滑處理得到最終的總生長高度.

Step3根據調整后的總生長厚度，得到每層生長高度，依次重復Step2.1—Step2.2，直到達到邊界層數nl，對最外層前沿進行全局相交檢測，采用收縮高度的方式處理.

2 前沿節點生長法向確定

2.1 定義

邊界層網格生成過程中的前沿為表面三角形網格.物面初始三角形網格剖分為初始前沿.前沿上節點 V，其周圍相鄰節點記為 Pi(i=1,2,···)，對應的鄰邊記為 ei(i=1,2,···).

V的鄰邊ei兩側的三角面片之間的有向二面角記為 θi，θi較大 (大于 240°)的邊定義為特征邊. 鄰邊中存在三條及以上特征邊的節點稱為角點.如圖2標出了邊e3及對應的二面角θ3.

2.2 生長法向計算

計算每個節點周圍的有向二面角，確定其鄰邊中的特征邊，根據特征邊將節點周圍的三角面片分組，如圖 2，V點周圍存在三條特征邊 VP1,VP3,VP4，則將周圍的三角面片分為3組.每一組確定一條生長方向ni(i=1,2,3)，由該組三角面片的外法向加權得到，Ni為第i組三角面片數，mj為第 j個三角片外法向，權值wj可簡單取為1或者采用文獻[26]中提出的角度權值

圖2 生長法向確定Fig.2 Calculate growth directions

以上的法向計算策略能夠很好地滿足推進單元正交性及可視化條件，而且三條法向就足夠了.因此本文設定節點的生長方向最多有三條.當節點周圍存在三條以上特征邊時，則從中挑出有向二面角最大的三條特征邊作為最后的分組標準.

節點生長法向確定后，周圍三角片在該節點處有一個對應的生成法向；特征邊在該點處有一到兩個生長方向，可由邊兩側的三角片的對應情況確定.當節點存在三個生長法向時，該節點會單獨生長出一個四面體單元.

3 生長高度確定

初始生長高度為用戶給定，即程序的輸入中包含的第一層網格高度hf，生長高度增長比率hr及層數nl.由于節點的生長法向不一定都垂直于物面，因此有必要對初始生長高度進行適當的調整.按照調整后的高度信息逐層推進生長邊界層，在推進過程中，由于邊界形狀的復雜性，可能出現單元自交及前沿全局相交的問題.因此也需要對當前的生長高度進行調節.

3.1 初始生長高度

在用戶給定初始生長高度等信息作為輸入后，由于上一節確定的節點生長方向不一定垂直于節點周圍所有的三角面片，因此為了盡量保證推進前沿與當前前沿的距離符合用戶給定的信息，需要對初始生長高度進行適當的調整.

對于二維情況，由于生長法向是節點外角的N等分線 (N為生長法向數目)，因此只需要考慮法向與物面垂線夾角即可.圖3為二維示意圖，前沿推進高度為h，節點V1有一個法向，則生長高度為h/sinθ，點V2處θ=90°，所以生長高度為h.

圖3 二維初始高度調節Fig.3 2D initial height setting

對于三維情況，由于采用多生長方向，同時法向的確定為分組三角面片法向的加權平均，因此為了保證推進前沿與當前前沿盡量平行，則需要采用加權法調整每個節點的初始生長高度.如圖4，節點V周圍的一個三角片Ti，其單元法向為nT，它對應的分組節點生長法向為ni(節點V可能有多個法向，圖4只顯示了對應于Ti的法向ni)，則初始高度調整后為

遍歷節點V周圍所有三角片加權平均得到節點V調整后的初始節點生長高度hV，記T(V)為節點V周圍的三角片集合，NT(V)為集合元素個數，則有

圖4 三維初始高度調節Fig.4 3D initial height setting

3.2 自交

自交即單元出現反轉，如圖5(a)所示.當出現自交時，收縮自交單元節點的生長高度，二維單元則收縮至交點處，即法向融合，這樣可以避免前沿小尺寸出現，同時也可以減少之后處理自交的頻次.

圖5 自交處理二維示意圖Fig.5 Dealing with 2D self-intersection

對于三維情況，當三棱柱單元出現反轉時，如圖6，記三角形ABC推進到三角形A′B′C′，面積由S變為S′，則采取收縮高度的方式處理，收縮因子γ為

自交處理完以后，需要采用拉普拉斯方法對前沿各點的生長高度進行光滑處理，這樣可以避免三棱柱單元鄰邊突變的情形.

圖6 自交處理三維示意圖Fig.6 Dealing with 3D self-intersection

3.3 全局相交

對于復雜的計算域，前沿層推進的過程中，不同物面推進的前沿之間可能存在相交的情況，對于前沿規模非常大的問題，全局相交檢測是一個耗時的過程.為了減少相交檢測的范圍，可以建立八叉樹(二維四叉樹)結構使得全局檢測局部化.盡管八叉樹結構提高了全局相交檢測的效率，但以往的層推進方法每一層都需要檢測，使效率下降.為了提高邊界層網格生成的效率，本文的方法先對全局預測一個總生長高度，即先采用總生長高度生長一層邊界層，在單元生長的過程中進行自相交處理，得到每個表面節點的生長高度，然后進行拉普拉斯光滑處理，最后進行一次全局相交檢測，采用收縮的方式得到每個節點的總生長高度，進而得到合適的每層生長高度.采用預測得到的生長高度進行邊界層網格推進則可以減少全局相交檢測的次數，提高效率.

圖7所示的是一個二維全局收縮的示意圖，收縮后的總生長高度能夠較好地避免內部邊界層網格前沿的全局相交，因此，在很大程度上能夠避免多次的全局相交檢測及收縮操作.

圖7 全局相交處理二維示意圖Fig.7 Dealing with global intersection(2D example)

4 單元生長及前沿尺寸調節

4.1 單元生長及細分

根據第2節可知，前沿存在三種特征信息:角點，特征邊以及三角片.因此，本文方法的單元生長方式共有如圖8四種情況，粗點及粗線在當前前沿上，箭頭表示生長方向:角點生長出一個四面體單元，三角片生長出一個三棱柱單元，特征邊則可以生長出一個三棱柱單元或者一個退化的三棱柱單元.

圖8 單元生長示意圖Fig.8 Element growing types

由于大多數流體求解器局限于純四面體網格[7]，為了適應大多數流體求解程序的要求，本文實現了將得到的邊界層網格轉化全四面體網格的過程:連接三棱柱或退化三棱柱的側面對角線，將三棱柱細分為三個各向異性的四面體單元，如圖9.

為了保持單元之間的相容性，即單元相鄰當且僅當單元具有共同的節點、邊或面.然而，不是所有的連接方式都能對三棱柱單元進行有效的剖分，如圖9(a)的S型連接方式則不能對三棱柱單元進行剖分得到三個四面體單元.因此在連接三棱柱側邊四邊形對角線的時候需要進行特殊的處理.

排除以上的S型連接方式，可以得到合法連接的充要條件是存在兩條有公共端點的對角線.如圖9(b)連接方式將三棱柱剖分為ABCF,ABFE,AEFD三個四面體.

圖9 三棱柱側面對角線連接方式Fig.9 Connecting types of diagonals for lateral faces of a prism

因此本文處理連接對角線的方式為：對角線統一由編號小的點連向編號大的生長點，如圖10編號i1＞i0＞i2，則對角線由Pi0和Pi2連向編號大的點Pi1的生長點，Pi2連向Pi0的生長點.對于全部的三棱柱和退化的三棱柱的側面都進行如此連接處理，則可以得到全局相容的四面體網格.

圖10 節點拓撲連接關系示意圖Fig.10 Diagonal choice for prism tetrahedronization

4.2 前沿尺寸調節

邊界層網格生成結束之后，邊界層最外層前沿將作為各向同性網格生成器的表面網格，因此最外層前沿三角網格的質量直接關系到生成的各向同性網格的質量.所以，很有必要在生成邊界層網格的同時進行前沿尺寸調節，為之后的各向同性網格生成提供高質量的三角形網格輸入.

在邊界層生長推進的過程中，其前沿的尺寸可能不斷發生變化.二維情況下，在生成邊界層網格時，表面線段在推進后的尺寸可能會變大(變小的情況按自交情況處理)，如圖11中的BC推進到AD，AD的長度變大，記

文獻[27]在生成二維混合網格時對這種情況進行了處理，提出當α大于給定值(本文取1.5)時，則取AD中點E，將四邊形單元ABCD細分為三個三角形.

本文將其推廣應用到三維尺寸控制中.由于三維邊界前沿在推進過程中只出現四種情形，角點產生四面體，特征邊和三角面片產生三棱柱，不管哪種情形，其表面只存在兩種形狀:三角形或四邊形.首先標記尺寸大的邊，三角形的處理方式為標記邊的中點與對應的頂點相連，然后消除標記，依次進行，直至所有標記邊處理結束;而對于四邊形，則采用圖11的細分方式，同時也加入三角形的操作方式;參照圖12，紅色粗線為標記的尺寸較大的邊，藍色線為新加入線.對應于三維中的調節操作可參見附錄A,另外有兩種情況需要進行一次邊置換操作，如圖13.

圖12 三角形及四邊形處理方式Fig.12 Subdivision of triangles and quads

圖13 邊置換操作(藍色線)Fig.13 Edge swap operation(blue line)

5 邊界層網格生成實例及數值驗證

程序在輸入物面三角形網格之后，能夠根據給定的邊界層網格第一層高度、層數及總高度，自動生成邊界層網格.以下實例的遠場網格均采用前沿推進法[28]生成.同時，為了驗證本文生成的混合網格的計算有效性，采用本文方法生成的混合網格對二維NACA0012翼型、三維ONERA M6翼型、帶發動機短艙F6翼身組合體進行黏性繞流計算.

圖14是自交處理和尺寸調節結合的二維實例圖.從圖中可以看出，調整生長高度融合生長法向之后，單元自交得到很好的處理(凹角區域).

圖14 自交處理和尺寸調節二維實例Fig.14 2D example dealing with self-intersection and size controlling

圖15和圖16為一個凹槽黏性網格圖，展示了凹角區域收縮生長高度的效果.

圖15 凹槽黏性網格Fig.15 Viscous mesh example of a cavity

圖16 凹槽黏性網格內部放大圖Fig.16 Enlarged inside view of the cavity mesh

5.1 NACA0012翼型

圖17是一個二維NACA0012翼型混合網格生成實例，生成20層邊界層網格，含8394個四邊形，102個三角形;遠場包括11302個三角形.圖18為翼型尾部放大圖.從二維NACA0012實例可以看出，采用尺寸調節后，最外層前沿線段長度比較均勻，因此邊界層網格與遠場網格銜接較好.

圖17 NACA0012二維混合網格Fig.17 2D hybrid mesh for NACA0012

圖18 NACA0012二維混合網格尾部放大圖Fig.18 Enlarged view of the 2D hybrid mesh of NACA0012

圖19是一個三維NACA0012混合網格生成實例，生成10層邊界層網格，包括36819個四面體；遠場包含22787個四面體.圖20為內部視圖，可以看出多生長方向提高了邊界層網格的正交性.

圖19 NACA0012三維網格Fig.19 3D hybrid mesh of NACA0012

圖20 NACA0012三維網格內部視圖Fig.20 Inside view of the 3D hybrid mesh of NACA0012

二維NACA0012翼型繞流的入口邊界馬赫數為0.8，攻角為0°，采用Spalart-Allmaras湍流模型，得到的壓力分布云圖如圖21所示，其上表面壁面附近的速度矢量圖如圖22所示，從結果可以看出本文采用的混合網格能夠捕捉壁面附近及流場內部的信息，網格是有效的.

圖21 二維NACA0012壓力分布云圖Fig.21 Pressure contours of NACA0012 wing

圖22 機翼上表面附近速度矢量圖Fig.22 Velocity vectors near upper wing surface

5.2 三維ONERA M6翼型

采用本文方法生成M6的黏性混合網格如圖23所示，圖24為內部視圖.

圖24 ONERA M6混合網格內部視圖Fig.24 Inside view of hybrid mesh of M6 wing

為了與現有實驗數據進行對比，三維M6翼型繞流采用標準實驗工況，Ma=0.8395，攻角為3.06°，采用k-ω湍流模型，機翼不同站位上的壓力系數分布計算值與實驗數據[29]對比如圖25～圖29.選取的站位y/b分別為0.2，0.44，0.65，0.8，0.9.

圖25 y/b=0.2壓力系數分布Fig.25 Pressure coeffients on the wing surface at y/b=0.2

圖26 y/b=0.44壓力系數分布Fig.26 Pressure coeffients on the wing surface at y/b=0.44

圖27 y/b=0.65壓力系數分布Fig.27 Pressure coeffients on the wing surface at y/b=0.65

圖28 y/b=0.8壓力系數分布Fig.28 Pressure coeffients on the wing surface at y/b=0.8

圖29 y/b=0.9壓力系數分布Fig.29 Pressure coeffients on the wing surface at y/b=0.9

從圖中對比可以看出，壓力系數分布總體上吻合得較好，各站位上的激波位置捕捉較為理想.總之，從對比結果可以看出本文生成的黏性混合網格能夠有效應用到黏性繞流計算中.

5.3 帶發動機短艙的F6翼身組合體

采用本文方法對F6翼身組合體進行黏性混合網格生成，圖30為內部網格剖分圖.

圖30 DLR-F6混和網格內部視圖Fig.30 Inside view of hybrid mesh of DLR-F6

計算中，為了與實驗數據進行對比，采用Ma=0.75，Re=3×106，α=1.738的計算條件，文獻[30]表明計算結果與湍流模型的選取有關，本文的湍流模型采用k-ω模型.將得到的機翼不同站位上的壓力系數與實驗值[31]進行對比.圖31～圖38展示了計算結果與所有8個實驗站位上的對比圖.從圖中可以看出，計算結果與實驗得到的物面壓力系數分布較為符合，同時激波位置的捕捉也較為準確，因此本文生成的混合網格運用到F6翼身組合體的繞流計算是有效的.

圖31 y/b=0.15壓力系數分布Fig.31 Pressure coeffients on the wing surface at y/b=0.15

圖32 y/b=0.239壓力系數分布Fig.32 Pressure coeffients on the wing surface at y/b=0.239

圖33 y/b=0.331壓力系數分布Fig.33 Pressure coeffients on the wing surface at y/b=0.331

圖34 y/b=0.377壓力系數分布Fig.34 Pressure coeffients on the wing surface at y/b=0.377

圖35 y/b=0.411壓力系數分布Fig.35 Pressure coeffients on the wing surface at y/b=0.411

圖36 y/b=0.514壓力系數分布Fig.36 Pressure coeffients on the wing surface at y/b=0.514

圖37 y/b=0.638壓力系數分布Fig.37 Pressure coeffients on the wing surface at y/b=0.638

圖38 y/b=0.847壓力系數分布Fig.38 Pressure coeffients on the wing surface at y/b=0.847

6 結論

本文提出的黏性邊界層網格自動生成策略可以很好地解決物面凹凸區域及角點處邊界層網格的生長問題，避免產生非法單元，同時減少全局相交檢測的次數，提高了網格生成的效率;另外還發展了一種三維尺寸調節方法，保證了邊界層網格最外層三角形網格的尺寸統一，有利于生成高質量的各向同性遠場網格.最后將本文生成的二維及三維的黏性混合網格應用到高雷諾數黏性繞流計算中，得到的計算結果表明本方法生成的黏性混合網格能夠滿足實際計算的要求，得到合理的計算結果.結合網格變形方法，如彈簧類方法，代數插值類方法[32]等,本文提出的混合網格生成策略也將很快應用到非定常流場計算[3334]中.

本文在處理自交或全局相交時只采用收縮因子的形式，在容易產生自交的凹角區域，其法向推進距離小于其他區域，前沿尺寸也相應減小，如圖39所示為凹角區域收縮效果，在該區域的未推進空間過小，生成各向同性網格會有一定困難.目前的處理方式是增加該區域表面網格的密度，減小表面網格尺寸，即生成的初始表面網格在凹凸脊及角點區域的單元會加密.同樣，在一些狹縫區域，兩側壁面同時推進的情況下，由于空間狹小，會造成未推進區域過小，給各向同性網格生成造成困難.本文實例中尚未涉及此種特別狹窄的區域，全局的前沿合并策略是解決這類情況的一種好的方式，將在今后的工作中進一步考慮和完善.

圖39 凹角網格收縮效果Fig.39 Mesh shrinking at concave ridges

另外，本文的程序將表面網格直接作為輸入，計算域邊界的表面三角剖分對于邊界層網格的順利進行具有重要意義[6,35]，因此，在生成邊界層網格之前，如何獲得適用于邊界層推進的高質量的表面網格值得進一步研究.

1 Sahni O,Jansen KE,Shephard MS,et al.Adaptive boundary layer meshing for viscous fl ow simulations.Engineering with Computers,2008,24(3):267-285

2 Pirzadeh S.Unstructured viscous grid generation by the advancinglayers method.AIAA Journal,1994,32(8):1735-1737

3 Kallinderis Y,Ward S.Prismatic grid generation for threedimensional complex geometries.AIAA Journal,1993,31(10):1850-1856

4 Lohner R.Matching semi-structured and unstructured grids for Navier–Stokes calculations.AIAA Paper 1993-3348,Orlando,FL,U.S.A.1993

5 Aubry R,L¨ohner R.On the‘most normal’normal.Communications in Numerical Methods in Engineering,2008,24(12):1641-1652

6 Aubry R,Dey S,Mestreau E,et al.An entropy satisfying boundary layer surface mesh generation.SIAM Journal on Scienti fi c Computing,2015,37(4):A1957-A1974

7 Aubry R,Dey S,Karamete K,et al.Smooth anisotropic sources with application to three-dimensional surface mesh generation.Engineering with Computers,2016,32(2):313-330

8 Dey S,Aubry R,Karamete B,et al.Capstone:A geometry-centric platform to enable physics-based simulation and design of systems.Computing in Science&Engineering,2016,18(1):32-39

9 Connell SD,Braaten ME.Semistructured mesh generation for threedimensional Navier-Stokes calculations.AIAA Journal,1995,33(6):1017-1024

10 Garimella RV,Shephard MS.Boundary layer mesh generation for viscous fl ow simulations. International Journal for Numerical Methods in Engineering,2000,49(1):193-218

11 Athanasiadis AN,Deconinck H.Object-oriented three-dimensional hybrid grid generation.International Journal for Numerical Methods in Engineering,2003,58(2):301-318

12 Athanasiadis AN,Deconinck H.A folding/unfolding algorithm for the construction of semi-structured layers in hybrid grid generation.Computer Methods in Applied Mechanics and Engineering,2005,194(48):5051-5067

13 Sharov D,Luo H,Baum JD,et al.Unstructured Navier-Stokes grid generation at corners and ridges.International Journal for Numerical Methods in Fluids,2003,43(6-7):717-728

14 王剛,葉正寅,陳迎春.三維非結構黏性網格生成方法.計算物理,2001,18(5):402-406(Wang Gang,Ye Zhengyin,Chen Yingchun.Generation of three-dimensional unstructured viscous grids.Chinese Journal of Computational Physics,2001,18(5):402-406(in Chinese))

15 王剛,葉正寅.三維非結構混合網格生成與N-S方程求解.航空學報,2003,24(5):385-390(Wang Gang,Ye Zhengyin.Generation of three dimensional mixed and unstructured grids and its application in solving Navier-Stokes equations.Acta Aeronautica et Astronautica Sinica,2003,24(5):385-390(in Chinese))

16 陳建軍,曹建,徐彥等.適應復雜外形的三維黏性混合網格生成算法.計算力學學報,2014,31(3):363-370(Chen Jianjun,Cao Jian,Xu Yan,et al.Hybrid mesh generation algorithm for viscous computations of complex aerodynamics con fi gurations.Chinese Journal of Computational Mechanics,2014,31(3):363-370(in Chinese))

17 Dyedov V,Einstein DR,Jiao X,et al.Variational generation of prismatic boundary-layer meshes for biomedical computing.International Journal for Numerical Methods in Engineering,2009,79(8):907-945

18 Jiao X.Face o ff setting:A uni fi ed approach for explicit moving interfaces.Journal of Computational Physics,2007,220(2):612-625

19 Ito Y,Nakahashi K.Improvements in the reliability and quality of unstructured hybrid mesh generation.International Journal for Numerical Methods in Fluids,2004,45(1):79-108

20 Alauzet F,Marcum D.A closed advancing-layer method with connectivity optimization-based mesh movement for viscous mesh generation.Engineering with Computers,2015,31(3):545-560

21 Zhang L,Zhao Z,Chang X,et al.A 3D hybrid grid generation technique and a multigrid/parallel algorithm based on anisotropic agglomeration approach.Chinese Journal of Aeronautics,2013,26(1):47-62

22 Zhong Z,Zhang L P,Xin H E.Hybrid grid generation technique for complex geometries based on agglomeration of anisotropic tetrahedrons.Acta Aerodynamica Sinica,2013,31(1):34-39

23 張來平,張涵信.復雜無粘流場數值模擬的矩形/三角形混合網格技術.力學學報,1998,30(1):104-108(Zhang Laiping,Zhang Hanxin.A Cartesian/triangular hybrid grid solver for complex inviscid fl ow fi elds.Acta Mechanica Sinica,1998,30(1):104-108(in Chinese))

24 張來平,王志堅,張涵信.動態混合網格生成及隱式非定常計算方法.力學學報,2004,36(6):664-672(Zhang Laiping,Wang Zhijian,Zhang Hanxin.Hybrid moving grid generation and implicit algorithm for unsteady fl ows.Acta Mechanica Sinica,2004,36(6):664-672(in Chinese))

25 張來平,常興華,段旭鵬,等.基于動態混合網格的不可壓非定常流計算方法.力學學報,2007,39(5):577-586(Zhang Laiping,Chang Xinghua,Duan Xupeng,et al.Implicit algorithm based on dynamic hybrid mesh for incompressible unsteady fl ows.Acta Mechanica Sinica,2007,39(5):577-586(in Chinese))

26 Th¨urrner G,W¨uthrich C A.Computing vertex normals from polygonal facets.Journal of Graphics Tools,1998,3(1):43-46

27 Dussin D,Fossati M,Guardone A,et al.Hybrid grid generation for two-dimensional high-Reynolds fl ows.Computers&Fluids,2009,38(10):1863-1875

28 Geuzaine C,Remacle JF.Gmsh:A 3-D fi nite element mesh generator with built-in pre-and post-processing facilities.International JournalforNumericalMethodsinEngineering,2009,79(11):1309-1331

29 Schmitt V,Charpin F.Pressure distributions on the ONERA-M6-Wing at transonic Mach numbers,experimental data base for computer program assessment.Report of the Fluid Dynamics Panel Working Group 04,AGARD AR-138,1979.

30 石磊,楊云軍,周偉江等.兩種湍流模型在高速旋轉翼身組合彈箭中的對比研究.力學學報,2017,49(1):84-92(Shi Lei,Yang Yunjun,Zhou Weijiang.A comparative study of two turbulence models for Magnus e ff ect in spinning projectile.Chinese Journal of Theoretical and Applied Mechanics,2017,49(1):84-92(in Chinese))

31 La fl in KR,Brodersen O.Summary of data from the Second AIAA CFD Drag Prediction Workshop.AIAA2004-0555,2004

32 劉中玉,張明鋒,聶雪媛等.一種基于徑向基函數的兩步法網格變形策略.力學學報,2015,47(3):534-538(Liu Zhongyu,Zhang Mingfeng,Wei Xueyuan,et al.A two-step mesh deformation strategy based on radial basis function.Chinese Journal of Theoretical and Applied Mechanics,2015,47(3):534-538(in Chinese))

33 常興華,馬戎,張來平等.基于計算流體力學的“虛擬飛行”技術及初步應用.力學學報,2015,47(4):596-604(Chang Xinghua,Ma Rong,Zhang Laiping,et al.Study on CFD-based numerical virtual fl ight technology and preliminary application.Chinese Journal of Theoretical and Applied Mechanics,2015,47(4):596-604(in Chinese))

34 童福林,李新亮,唐志共.激波與轉捩邊界層干擾非定常特性數值分析.力學學報,2017,49(1):93-104(Tong Fulin,Li Xinliang,Tang Zhigong.Numerical analysis of unsteady motion in shock wave/transitional boundary layer interaction.Chinese Journal of Theoretical and Applied Mechanics,2017,49(1):93-104(in Chinese))

35 Aubry R.Ensuring a smooth transition from semi-structured surface boundary layer mesh to fully unstructured anisotropic surface mesh//53rd AIAA Aerospace Sciences Meeting and Exhibit,AIAA,Reston,VA,2015,AIAA-2015-1507

AUTOMATIC VISCOUS BOUNDARY LAYER MESH GENERATION1)

Gan Yangke Liu Jianfei2)
(Department of Mechanics and Engineering Science，College of Engineering，Peking University，Beijing 100871,China)

High Reynolds number fl uid fl ows have boundary layers at the wall.To automatically generate robust and valid boundary layer mesh for the simulations is still the bottleneck problem of computational fl uid dynamics.Prisms/Tetrahedra hybrid mesh leads to signi fi cant savings in mesh size and solution costs.However,it’s still difficult to generate prismatic elements of high quality within boundary layers of complex models.Previous advancing layers techniques sometimes lead to invalid meshes and poor quality elements at concave/convex ridges and sharp corners.To improve these situations,we present a strategy for automatically generating viscous boundary layer mesh.In this method,multiple growth directions at ridges and corners are well de fi ned by the dihedral angles around the vertices and the growth heights are adjusted appropriately.Therefore,boundary layer mesh grows well at sharp corners,convex and concave ridges of the domain.We also decrease the number of global intersection checks by prede fi ning the total growth heights before generating elements through one global check,which improves the efficiency of mesh generation.At the same time,we develop a 3D strategy of mesh size control to get a size uniform triangular mesh of the outer boundary of the boundary layer mesh,which is bene fi cial to generate far- fi eld isotropic mesh of high quality.Finally,mesh examples and the viscous fl ow simulations including 2D and 3D are presented.In 3D,the hybrid mesh over the standard ONERA M6 and DLR-F6 con fi gurations are generated with the present method.The numerical results agree very well with experiment data which indicates that the hybrid viscous meshes generated by the proposed method are e ff ective and efficient.

boundary layer mesh，growth normal，feature edge，intersection，size control

O242.21

A

10.6052/0459-1879-17-154

2017–05–04收稿，2017–07–05 錄用,2017–07–07 網絡版發表.

1)國家自然科學基金(11172005)和國家重點基礎研究發展計劃(2010CB832701)資助項目.

2)劉劍飛，副教授，主要研究方向：網格生成，計算幾何，計算機圖形學.E-mail：liujianfei@pku.edu.cn

甘洋科,劉劍飛.黏性邊界層網格自動生成.力學學報,2017,49(5):1029-1041

Gan Yangke,Liu Jianfei.Automatic viscous boundary layer mesh generation.Chinese Journal of Theoretical and Applied Mechanics,2017,49(5):1029-1041

附錄A三維邊界層網格生成中單元尺寸調節操作

角點生成的四面體Tetraheron generated from one corner

三角片生成的三棱柱Prism generated from one triangle

特征邊生成的三棱柱Prism generated from one feature edge