趙 鐘，張來平，2，赫 新，2

（1．中國空氣動力研究與發(fā)展中心 計算空氣動力研究所，四川 綿陽 621000；2．空氣動力學(xué)國家重點實驗室，四川綿陽621000）

0 引 言

在CFD工業(yè)化應(yīng)用中，生成計算網(wǎng)格一直占據(jù)整個計算流程的60%～70%時間，這一點對于復(fù)雜外形高質(zhì)量粘性流動計算網(wǎng)格的生成更為突出［1］。

依照網(wǎng)格的拓?fù)浣Y(jié)構(gòu)，計算網(wǎng)格分為結(jié)構(gòu)網(wǎng)格和非結(jié)構(gòu)網(wǎng)格。結(jié)構(gòu)網(wǎng)格的優(yōu)點是數(shù)據(jù)結(jié)構(gòu)簡單，存儲方便，計算簡單快捷，計算結(jié)果精度高；其缺點是難以處理復(fù)雜外形，對一些外形上的奇點需要進(jìn)行簡化或特殊處理。非結(jié)構(gòu)網(wǎng)格最早來源于固體力學(xué)中的有限元計算，主要指二維的三角形和三維的四面體網(wǎng)格，其優(yōu)點是易于處理復(fù)雜外形，數(shù)據(jù)結(jié)構(gòu)的隨機性使得其易于進(jìn)行網(wǎng)格自適應(yīng)；缺點是存儲量大，計算效率低。非結(jié)構(gòu)網(wǎng)格的缺點在高雷諾數(shù)計算時表現(xiàn)得尤為突出：高雷諾數(shù)邊界層模擬要求在物面法向有足夠的網(wǎng)格分辨率，這一點使得非結(jié)構(gòu)網(wǎng)格數(shù)量急劇增長。為了在一定程度上減少網(wǎng)格量，可以使用法向高度壓縮的“各向異性”四面體單元。然而，法向相鄰兩個各向異性四面體單元的中心連線和物面法向的偏離使得計算在某些情況下不夠精確（如利用格心型有限體積法求解單元內(nèi)的物理量梯度時）。

基于結(jié)構(gòu)網(wǎng)格和非結(jié)構(gòu)網(wǎng)格各自的優(yōu)缺點發(fā)展起來的混合網(wǎng)格技術(shù)結(jié)合了二者的優(yōu)勢：邊界層內(nèi)的網(wǎng)格一般采用半結(jié)構(gòu)的三棱柱單元（在物面法向是結(jié)構(gòu)的，而在貼體方向是非結(jié)構(gòu)的），邊界層以外單元是各向同性的四面體網(wǎng)格或Cartesian網(wǎng)格。與結(jié)構(gòu)網(wǎng)格相比，生成復(fù)雜外形的混合網(wǎng)格更加容易；和純非結(jié)構(gòu)網(wǎng)格相比，邊界層法線上有很高的網(wǎng)格分辨率和很好的網(wǎng)格正交性，又極大地減少了網(wǎng)格量（邊界層內(nèi)比非結(jié)構(gòu)網(wǎng)格減少大約三分之二）。這無疑對提高復(fù)雜外形的數(shù)值模擬精度和效率都是非常有利的。

混合網(wǎng)格的出現(xiàn)，在提高邊界層網(wǎng)格質(zhì)量的同時，也在一定程度上緩解了生成復(fù)雜外形網(wǎng)格帶來的困難。目前常用的三棱柱網(wǎng)格生成方法主要有幾何層推進(jìn)方法［2－3］和 求 解 雙 曲 型 方 程 的 方 法［4－5］。 這 些方法在實際工程應(yīng)用中得到了成功的應(yīng)用。但是，對于工業(yè)應(yīng)用中很多極端復(fù)雜的實際外形而言，要生成邊界層的三棱柱網(wǎng)格并非易事，往往會在幾何曲率變化劇烈的區(qū)域出現(xiàn)網(wǎng)格相交的情況，為此需要付出很多人工勞動調(diào)整局部網(wǎng)格分布，極大地制約了CFD的應(yīng)用。

鑒于復(fù)雜外形三棱柱／四面體混合網(wǎng)格生成的困難，本文發(fā)展了一種自動化的基于“各向異性”四面體網(wǎng)格聚合的混合網(wǎng)格生成方法（為了敘述方便，下文稱“聚合法”），即首先用已經(jīng)發(fā)展成熟的非結(jié)構(gòu)網(wǎng)格生成方法生成粘性層內(nèi)的各向異性四面體單元和無粘區(qū)的各向同性四面體單元，然后將邊界層內(nèi)的各向異性四面體單元聚合為三棱柱，而外場的各向同性四面體單元保持不變。這種方法生成的網(wǎng)格在整個計算域是由四面體、三棱柱和少數(shù)多面體單元組成的混合網(wǎng)格。通過對DLR－F6翼身組合體和復(fù)雜戰(zhàn)斗機外形的混合網(wǎng)格生成及數(shù)值模擬表明，該方法生成的混合網(wǎng)格具有網(wǎng)格質(zhì)量好、自動化程度高等特點，對提高復(fù)雜外形飛行器關(guān)鍵氣動參數(shù)的模擬精度具有重要意義。

1 非結(jié)構(gòu)四面體網(wǎng)格生成方法簡介

目前，非結(jié)構(gòu)的四面體網(wǎng)格生成方法已經(jīng)發(fā)展的比較成熟［6－7］，多種商業(yè)軟件業(yè)已集成了這種方法，比如Gridgen、ICEM－CFD等。這里僅做簡要介紹。非結(jié)構(gòu)的四面體網(wǎng)格包括粘性層內(nèi)的各向異性四面體網(wǎng)格和無粘區(qū)域的各向同性四面體網(wǎng)格。

生成各向異性四面體網(wǎng)格采用的是由陣面推進(jìn)法（Advancing Front Method）改進(jìn)而來的層推進(jìn)法（Advancing Layer Method）。在生成好物面的三角形單元后，通過確定三角形頂點的法向和推進(jìn)距離生成出拉伸的四面體單元，生成步驟如下：

（1）計算物面三角形單元頂點的推進(jìn)方向。通過面積加權(quán)點p周圍的面法向得到點p的推進(jìn)法向，并做一定次數(shù)迭代進(jìn)行光滑：

這里，Np表示點p周圍的頂點數(shù)目，Nf表示和點p相連的三角形數(shù)目表示p點的法向，nf表示面的法向，n表示第n次迭代。

（2）在確定了推進(jìn)方向后，需要確定推進(jìn)距離，一般采用幾何指數(shù)增長模式：

這里，δ0和δn分別表示第一層的網(wǎng)格高度和第n層的網(wǎng)格高度，r是法向網(wǎng)格距離變化率。

（3）在確定了推進(jìn)方向和推進(jìn)距離后，從每個頂點向計算域推出一層網(wǎng)格點。對于每個三角形而言，三個頂點向外推出三個點。分別將推出的三點和前一層的三個點相連，形成三個四面體單元。

在推進(jìn)過程中，為了避免單元之間相交，保證交界面的相容性，需要兩兩判斷是否相交，這在三維情況下會帶來網(wǎng)格生成效率急劇下降。為此，可以采用“彈簧法”來判斷相交性以提高效率，具體過程可參考文獻(xiàn)［7］。

在推出給定層數(shù)的各向異性四面體網(wǎng)格后，停止各向異性四面體網(wǎng)格生成過程，形成各向異性四面體網(wǎng)格外邊界面。各向異性四面體網(wǎng)格推進(jìn)過程中，在局部會由于可能相交而停止推進(jìn)，這些局部停止推進(jìn)而形成的邊界面連同其他邊界面（如遠(yuǎn)場、對稱面）一起，形成了各向同性四面體網(wǎng)格生成的初始陣面，從而進(jìn)入各向同性四面體網(wǎng)格生成階段。生成空間各向同性四面體單元一般采用Delaunay法或陣面推進(jìn)法。在現(xiàn)有的商業(yè)軟件中，一般將二者結(jié)合起來使用。這樣既可以提高四面體網(wǎng)格生成效率，同時可以保證初始陣面的相容性。具體方法請參見文獻(xiàn)［6－7］，這里不再詳述。

2 基于各向異性四面體聚合的混合網(wǎng)格生成方法

在四面體網(wǎng)格生成好后，就可以用本文發(fā)展的聚合法來生成混合網(wǎng)格。聚合法的基本思想是，將各向異性四面體聚合為三棱柱，而各向同性四面體保持不變。這種方法一方面充分利用了非結(jié)構(gòu)網(wǎng)格填充復(fù)雜外形計算域的自動化特點，另一方面又能聚合得到高質(zhì)量的混合網(wǎng)格以適應(yīng)粘性流動計算的需要。聚合過程主要包括四面體單元聚合、棱柱側(cè)面聚合和邊界條件處理。

2．1 四面體單元聚合

聚合四面體網(wǎng)格是整個算法的核心部分，其原理是充分利用層推進(jìn)法生成各向異性四面體算法的特點，將三個各向異性四面體單元聚合為一個三棱柱。聚合方法直接影響到聚合后的網(wǎng)格質(zhì)量，具體包括幾個步驟：

首先，提取單元的面幾何特征。對初始四面體網(wǎng)格，計算每個四面體的所有表面（三角形）的面積，并找出每個單元的面積最大面、次大面和最小面。

然后，根據(jù)四面體單元面的幾何特征判斷是否聚合，具體算法如下：

（1）根據(jù)單元的最小面和最大面的面積之比確定單元性質(zhì)：

若r＜α，則該單元為各向異性，否則為各向同性。這里，Smin和Smax分別表示最小面面積和最大面面積。α是一經(jīng)驗參數(shù)，實踐表明，α一般取值小于0．4，α取值越大則三棱柱層數(shù)越多，反之則邊界層網(wǎng)格和無粘區(qū)域網(wǎng)格過渡越光滑；

（2）第一次聚合。遍歷初始網(wǎng)格的所有面，設(shè)面Face兩側(cè)的單元分別為C1和C2，若Face是C1或者C2的最大面，并且C1或者C2是各向異性網(wǎng)格，則將這兩個網(wǎng)格聚合為一個粗網(wǎng)格CC。定義每個粗網(wǎng)格里含有的子單元個數(shù)為粗網(wǎng)格聚合率，則CC的聚合率為2；

（3）第二次聚合。遍歷初始網(wǎng)格的所有面，找出其兩側(cè)單元，若面Face兩側(cè)單元中有一個（C1）已被聚合而另一個（C2）尚未聚合，其中C1所在粗網(wǎng)格為CC，F(xiàn)ace是C1或C2的最大面，C1或C2是各向異性單元，并且CC的聚合率小于3，則將C2聚合到CC，CC聚合率加1。

（4）第三次聚合。遍歷每個尚未聚合的各向異性單元C1，找出其最大面以及最大面另一側(cè)單元C2，和C2所對應(yīng)的粗網(wǎng)格為CC。若CC的聚合率小于等于3，則將C1聚合到CC。

（5）遍歷所有初始四面體單元Cell，若Cell尚未被聚合，則將之單獨聚合為一個粗網(wǎng)格。

聚合過程如圖1所示。第一次聚合和第二次聚合的主要目的是將三個四面體單元聚合為一個三棱柱單元。但是對于外形變化劇烈的局部網(wǎng)格幾何特征不明顯，導(dǎo)致其不能被聚合而成為孤立單元。如果放任之，則該單元和周圍單元的體積比約為1∶3，網(wǎng)格過渡不光滑。為了避免網(wǎng)格過渡不光滑，第三次聚合將孤立單元聚合到周圍的粗網(wǎng)格，使之和周圍單元體積比約為4∶3，從而使網(wǎng)格光滑過渡。

圖1 聚合網(wǎng)格過程Fig．1 Procedure of cell agglomeration

2．2 側(cè)面的聚合

初始網(wǎng)格在經(jīng)過聚合后形成的“三棱柱”側(cè)面還是兩個三角形。為了減少存儲開銷和計算量，需要將每個側(cè)面的兩個三角形合并為一個四邊形，如圖2。聚合面的關(guān)鍵在于判斷哪些面需要聚合為一個面，可以通過面的左右單元屬性和面法向夾角判斷。值得注意的是，在合并邊界面時，只有邊界類型相同的兩個面才能合并，比如對稱面，否則會破壞幾何外形和邊界條件。

圖2 合并三棱柱側(cè)面Fig．2 Procedure of face agglomeration

邊界條件的處理比較簡單，將初始網(wǎng)格的邊界條件類型直接賦予聚合后的粗網(wǎng)格即可。在處理好邊界條件后就可以直接輸出網(wǎng)格并進(jìn)行計算。

2．3 復(fù)雜外形混合網(wǎng)格生成實例

為了驗證本文的混合網(wǎng)格生成方法對復(fù)雜外形的適應(yīng)性，本文首先選取了第二屆AIAA CFD阻力預(yù)測會議（DPW－II）［8］中的 DLR－F6翼身組合體進(jìn)行驗證。如圖3所示為該外形數(shù)模和局部表面網(wǎng)格。初始網(wǎng)格是全四面體網(wǎng)格，包括各向異性四面體網(wǎng)格和各向同性四面體網(wǎng)格，共1743萬個四面體單元，3496萬個面，295萬個點。在聚合過程中，α取為0．4，聚合后得到781萬個網(wǎng)格單元，2080萬個面，網(wǎng)格單元總量減少大約一半，面的總量大約減少三分之一。如圖4所示為聚合后的混合網(wǎng)格在x＝0．13截面和z＝0．3截面的網(wǎng)格視圖。聚合后邊界層有40層三棱柱網(wǎng)格，外場是四面體網(wǎng)格，從邊界層到外場網(wǎng)格均勻過渡。通過聚合過程不僅網(wǎng)格質(zhì)量得到提高，網(wǎng)格數(shù)量也大幅減少，計算效率也會有很大提高。

圖3 DLR－F6數(shù)模和局部表面網(wǎng)格Fig．3 Configuration and surface grids of DLR－F6

為了說明聚合法的魯棒性，本文生成了具有復(fù)雜外形的類F16戰(zhàn)斗機全機外形的混合網(wǎng)格。如圖5（a）所示為該外形的表面網(wǎng)格，其幾何構(gòu)型包括機身、機翼、平尾、立尾、腹鰭。我們采用聚合法生成該外形的混合網(wǎng)格量約為全場500萬。圖5（b）－（d）分別顯示了典型截面的混合網(wǎng)格。通過該實例表明，聚合法有很好的魯棒性，生成的混合網(wǎng)格可以保證邊界層內(nèi)絕大部分是三棱柱，而且從三棱柱網(wǎng)格光滑地過渡到外場四面體網(wǎng)格。

3 混合網(wǎng)格生成方法的應(yīng)用

為了驗證本文的混合網(wǎng)格生成方法的實用性，我們對DLR－F6（帶發(fā)動機和翼吊）外形進(jìn)行數(shù)值計算，并與風(fēng)洞試驗結(jié)果和他人計算結(jié)果進(jìn)行了比較。本文采用的流場解算器是作者所在團隊開發(fā)的CFD解算軟件 USTAR［9－10］。USTAR 是一款基于非結(jié)構(gòu)／混合網(wǎng)格的格心型有限體積雷諾平均N－S方程定常／非定常解算器。軟件采用迎風(fēng)型格式進(jìn)行空間無粘項離散，粘性項采用中心格式離散。時間推進(jìn)采用隱式Block LU－SGS算法［10］。為了提高計算效率，采用多重網(wǎng)格技術(shù)加速；對于大規(guī)模計算問題，采用了基于METIS分區(qū)的并行計算方法。軟件集成了典型湍流模型，包括：一方程SA模型，二方程k－ε模型、k－ω模型和SST模型。該軟件已在多種復(fù)雜外形的數(shù)值模擬中得到成功的應(yīng)用。其中的具體算法請參見文獻(xiàn)［9－10］。

計算的來流條件為DPW－II中的CASE2，即Ma＝0．75，Re＝3．0×106，攻角從－3．0°到1．5°。計算過程中采用SST湍流模型，沒有考慮轉(zhuǎn)捩影響。圖6所示為聚合后混合網(wǎng)格和聚合前非結(jié)構(gòu)網(wǎng)格的殘差、俯仰力矩收斂曲線對比，可以看到混合網(wǎng)格不管在收斂性還是計算效率上都有很大提高。

圖6 收斂歷程曲線Fig．6 Convergence history

DPW會議的重點是氣動力精確預(yù)測。對于氣動力預(yù)測來說，網(wǎng)格質(zhì)量尤其是邊界層內(nèi)網(wǎng)格質(zhì)量的好壞直接影響到計算精度。為了驗證聚合法生成的混合網(wǎng)格計算結(jié)果，本文不僅對比了風(fēng)洞實驗結(jié)果，還與DPW－II中一些參與者的計算結(jié)果進(jìn)行了對比［11－14］，包 括 NASA Langley 的 USM3D 和FUN3D，Cessna的非結(jié)構(gòu)解算器NSU3D和日本Kawasaki Heavy Industries的UG3。圖7是聚合后的混合網(wǎng)格（USTAR）、聚合前的非結(jié)構(gòu)四面體網(wǎng)格（Unstructured）、USM3D、FUN3D、NSU3D 計算結(jié)果對比。由圖可見，聚合前的非結(jié)構(gòu)網(wǎng)格計算結(jié)果和其他計算結(jié)果普遍偏差較大，證明了聚合后有效提高了網(wǎng)格質(zhì)量。

圖7（a）是計算得到的升力系數(shù)曲線。由圖可見，USTAR計算結(jié)果在攻角大于0．5°時和實驗結(jié)果符合較好，在攻角小于0．5°時與實驗值有一定偏差，但是升力曲線和FUN3D、USM3D的結(jié)果符合較好。DPW－II會議中，大部分提交的計算結(jié)果升力曲線也符合這種趨勢，即升力計算值比實驗值偏大，斜率比實驗值偏小［15］，這種現(xiàn)象在DPW－I中也普遍存在［15］。

圖7（b）是升阻力極曲線。第一、二屆阻力預(yù)測會議中，多數(shù)參與者提交的計算結(jié)果中阻力都比實驗值偏大，如圖所示，NSU3D和FUN3D的阻力都比實驗值偏大，但USTAR計算的極曲線和實驗值符合的較為一致。

圖7（c）是俯仰力矩曲線。圖中計算結(jié)果都和實驗值符合的不是很理想，USTAR的俯仰力矩曲線和USM3D符合較好，處于FUN3D和NSU3D之間。和阻力預(yù)測會議收到的大部分提交結(jié)果一樣，除FUN3D以外，俯仰力矩都比實驗值偏大［15］。

除了氣動力之外，壓力分布也是本文觀察的一個重點。圖8是機翼沿展向不同站位處的壓力分布，圖中為USTAR在混合網(wǎng)格上的計算結(jié)果與實驗值、USM3D、FUN3D的計算結(jié)果對比。在z／b＝0．150處，由于機翼上表面后緣翼身結(jié)合處分離渦的影響，計算結(jié)果和實驗值都有一定偏差；在z／b＝0．331站位處，由于機翼下表面的發(fā)動機內(nèi)側(cè)有分離渦，計算結(jié)果也和實驗有差異，但USTAR和USM3D明顯要比FUN3D符合的好；在z／b＝0．638站位，該處外形簡單流場結(jié)構(gòu)單一，因此三者都和實驗值均符合較好。

圖7 DLR－F6翼身組合體氣動特性計算Fig．7 Aerodynamic coefficients of DLR－F6

圖8 機翼不同站位處壓力系數(shù)分布Fig．8 Pressure coefficient distributions at z／b＝0．150，0．331，0．638

為了觀察流場結(jié)構(gòu)，圖9畫出了物面流線。圖9（a）為機翼上翼面流線，圖9（b）為翼身結(jié)合處流線，可看到有明顯的分離現(xiàn)象。圖9（c）為USTAR計算的下翼面發(fā)動機短艙內(nèi)側(cè)的流線，圖9（d）為UG3計算結(jié)果，二者現(xiàn)象一致，分離區(qū)大小也差不多，但是和實驗結(jié)果（圖9e）相比，分離區(qū)明顯偏大，在第二屆阻力預(yù)測會議的總結(jié)文獻(xiàn)里也指出大多計算結(jié)果具有分離區(qū)偏大的特點［15］。

圖9 USTAR、UG3計算的物面流線和實驗結(jié)果Fig．9 Streamlines over DLR－F6

4 結(jié)束語

本文針對復(fù)雜外形粘性流動網(wǎng)格生成費時費力的問題，發(fā)展了一種魯棒較好、自動化程度較高的混合網(wǎng)格生成方法。聚合法生成的網(wǎng)格集成三棱柱、四面體和少數(shù)的多面體，在減少網(wǎng)格數(shù)量的同時提高了網(wǎng)格質(zhì)量以增強對流場結(jié)構(gòu)的捕捉能力。通過采用第二屆阻力預(yù)測會議的標(biāo)準(zhǔn)算例對該方法進(jìn)行了驗證，結(jié)果表明該方法生成的混合網(wǎng)格質(zhì)量適應(yīng)計算的需要。通過對類F16復(fù)雜戰(zhàn)斗機全機構(gòu)型的混合網(wǎng)格生成，進(jìn)一步說明其魯棒性較好，表明其具有廣闊的應(yīng)用前景。

［1］ BAKER T J．Mesh generation：art or science［J］．Progress in Aerospace Sciences．2005，41：29－63．

［2］ PIRZADEH S．Three－dimensional unstructured viscous grids by the advancing－layers method［J］．AIAA Journal，1996，34：43－49．

［3］ KALLINDERIS Y，KHAWAJA A．Hybrid prismatic tetrahedral grid generation for complex geometries［J］．AIAA Journal，1996，34：291－298．

［4］ MATSUNO K．Hyperbolic upwind method for prismatic grid generation［R］．AIAA2000－1003，2000．

［5］ MATSUNO K．High order upwind method for hyperbolic grid generation［J］．Computers ＆ Fluids，1999，28：825－851．

［6］ LOHNER R，PARIKH P．Generation of three dimensional unstructured grids by the advancing front method［J］．Int．J．Numer．Meth．Fluids，1988，8：1135－49．

［7］ PIRZADEH S．Unstructured viscous grid generation by Advancing－Layers Method［R］．AIAA 93－3453，1993．

［8］ 2nd AIAA CFD Drag Prediction Workshop［EB］．http：／／aaac．larc．nasa．gov／tsab／cfdlarc／aiaa－dpw／，June 2003．

［9］ ZHANG L P，YANG Y J，ZHANG H X．Numerical simulations of 3Dinviscid／viscous flow fields on Cartesian／unstructured／prismatic hybrid grids［R］．In Proceedings of the 4th Asian CFD Conference［C］，Mianyang，Sichuan，China，September 2000．

［10］ZHANG L P，WANG Z J．A block LU－SGS implicit dual time－stepping algorithm for hybrid dynamic meshes［J］．Computers ＆ Fluids，2004，33：891－916．

［11］YAMAMOTO K，OCHI A．CFD sensitivity of drag prediction on DLR－F6configuration by structured method and unstructured method［R］．AIAA2004－398，2004．

［12］MAVRIPLIS D J．Drag prediction of DLR－F6using the turbulent Navier－Stokes calculations with multigrid［R］．AIAA2004－397，2004．

［13］LEE－RAUSCH E M，MAVRIPLIS D J．Transonic drag prediction on a DLR－F6transport configuration using unstructured solvers［R］．AIAA2004－554，2004．

［14］SCLAFANI A J，DEHAAN M A．OVERFLOW drag prediction for the DLR－F6transport configuration：a DPW－II case study［R］．AIAA 2004－393，2004．

［15］LAFLIN K R，BRODERSEN O．Summary of data from the second AIAA CFD Drag Prediction Workshop［R］．AIAA2004－0555，2004．