重疊網(wǎng)格隱式挖洞方法的改進(jìn)與應(yīng)用

江維青，楊愛明

(復(fù)旦大學(xué) 航空航天系，上海 200433)

在計算流體力學(xué)領(lǐng)域，針對結(jié)構(gòu)復(fù)雜的幾何體生成高質(zhì)量的單塊網(wǎng)格比較困難，而重疊網(wǎng)格技術(shù)[1]則是針對幾何體的各個部件分別生成單塊網(wǎng)格，各部件網(wǎng)格之間可以任意相交，重疊區(qū)域的網(wǎng)格單元經(jīng)過預(yù)處理后再從其他網(wǎng)格塊中插值獲得流場信息.因?yàn)槊恳痪W(wǎng)格塊的拓?fù)浣Y(jié)構(gòu)簡單，重疊網(wǎng)格方法大幅降低了復(fù)雜外形網(wǎng)格生成的難度.該方法在過去的30年里被廣泛應(yīng)用于航空航天飛行器流場的數(shù)值模擬，并在近些年延伸到用于模擬風(fēng)力發(fā)電機(jī)[2]、水下潛航器[3]和直升機(jī)全機(jī)[4]等復(fù)雜流場.

在重疊網(wǎng)格的預(yù)處理過程中，需要將重疊區(qū)域的網(wǎng)格單元分為計算單元、插值單元與非計算單元，這一過程也被稱為“挖洞”.傳統(tǒng)的挖洞方法需要識別出落在物面內(nèi)部的網(wǎng)格單元，將它們標(biāo)記為非計算單元從而形成洞區(qū)，位于洞區(qū)邊界上的點(diǎn)則作為插值單元從其他網(wǎng)格塊獲得流場信息，再通過割補(bǔ)法[5]、陣面推進(jìn)法[6]等方法優(yōu)化洞區(qū)邊界以獲得更好的插值效果.

相比傳統(tǒng)方法，由Baeder等最早提出的隱式挖洞(Implicit Hole Cutting，IHC)方法[7]則是根據(jù)網(wǎng)格單元質(zhì)量來判斷單元類型，有著自動化程度高、插值區(qū)域網(wǎng)格匹配度好等優(yōu)點(diǎn)，是目前較受關(guān)注的重疊網(wǎng)格挖洞方法之一.本文針對隱式挖洞方法的特點(diǎn)進(jìn)行了研究與改進(jìn)，選擇物面距與網(wǎng)格體積的乘積作為網(wǎng)格單元質(zhì)量判斷標(biāo)準(zhǔn)，增加識別物面內(nèi)網(wǎng)格單元與去除多余插值單元的步驟，通過建立網(wǎng)格塊挖洞優(yōu)先級矩陣以及基于模板跳躍法的搜索方法提高了貢獻(xiàn)單元搜索效率，并以復(fù)雜飛行器構(gòu)型的重疊網(wǎng)格生成與流場計算對本文發(fā)展的方法進(jìn)行了驗(yàn)證.

1 隱式挖洞方法簡介

隱式挖洞方法的核心是一個基于網(wǎng)格單元質(zhì)量做選擇的過程，除了給定網(wǎng)格本身信息之外不需要過多的人工干預(yù).其過程可以分為兩步.(1) 首先是在網(wǎng)格重疊區(qū)域計算所有網(wǎng)格單元的質(zhì)量，該參數(shù)表征了一個網(wǎng)格單元適合成為貢獻(xiàn)單元的程度，可以是網(wǎng)格單元的體積、物面距、長寬比等.(2) 然后為所有落在重疊區(qū)域內(nèi)的網(wǎng)格搜索貢獻(xiàn)單元，若存在質(zhì)量更好的貢獻(xiàn)單元，則該單元成為插值單元，從質(zhì)量最好的貢獻(xiàn)單元處通過三線性插值獲取流場信息；若貢獻(xiàn)單元不存在或貢獻(xiàn)單元的質(zhì)量均比該單元差，則該單元成為計算單元正常參與流場計算，這樣就避免了傳統(tǒng)方法中繁瑣的洞區(qū)邊界確定與調(diào)整過程.

本文中網(wǎng)格單元的質(zhì)量選取為網(wǎng)格體積與物面距的乘積，其數(shù)值越小則網(wǎng)格質(zhì)量越好.引入物面距作為質(zhì)量判斷標(biāo)準(zhǔn)可以有效地減少挖洞結(jié)束后孤點(diǎn)的數(shù)量[8]，而且靠近物面處的網(wǎng)格單元一般來說體積與物面距均較小，因此會被判別為高質(zhì)量的計算單元，使得插值邊界自然地遠(yuǎn)離了物面處流場變化梯度較大的區(qū)域，從而確保了更好的插值精度.

2 隱式挖洞方法的改進(jìn)

2.1 減少計算單元與插值單元數(shù)量

原始的隱式挖洞方法僅僅通過比較網(wǎng)格質(zhì)量區(qū)分插值單元與計算單元，不會挖去落在物體內(nèi)部的非計算單元.如圖1所示，繞圓形的曲線網(wǎng)格與均勻的笛卡爾背景網(wǎng)格互相重疊，以網(wǎng)格面積作為質(zhì)量判斷標(biāo)準(zhǔn)則曲線網(wǎng)格質(zhì)量較高，背景網(wǎng)格相對而言質(zhì)量較差.從圖1可見隱式挖洞后落在圓形物面內(nèi)部的背景網(wǎng)格被標(biāo)記為計算單元，不過因?yàn)檫@些物面內(nèi)部單元被插值單元所包圍，所以其計算得到的錯誤信息不會傳播到外界，這也是隱式挖洞方法不需要像傳統(tǒng)方法那樣識別出落在物面內(nèi)部網(wǎng)格的原因[9].這些多余的計算單元雖然不影響計算結(jié)果的準(zhǔn)確性，但卻不利于最終的流場顯示，也增加了不必要的計算量，應(yīng)當(dāng)采用適當(dāng)?shù)姆椒ǔ?

圖1 原始的隱式挖洞方法結(jié)果示例Fig.1 Illustration of original implicit hole cutting method

對于落在物面內(nèi)部的單元，可以在隱式挖洞開始前使用洞映射方法[10]識別并標(biāo)記為非計算單元，在之后的隱式挖洞過程中即可跳過這些單元，這樣雖然增加了網(wǎng)格預(yù)處理時的工作量，但確保了最終形成的網(wǎng)格不包含非必要的計算單元，也有利于流場計算結(jié)果的正確顯示.

由于重疊區(qū)域所有質(zhì)量較差的網(wǎng)格都從質(zhì)量較好的網(wǎng)格處插值獲得流場信息，原始的隱式挖洞方法最終得到的插值單元數(shù)量較大，而其中的許多插值單元并非必要.針對這一問題可以采用如下措施進(jìn)行改進(jìn)： 網(wǎng)格質(zhì)量較差的單元不再全部被標(biāo)記為插值單元，而是先統(tǒng)一標(biāo)記為非計算單元，再根據(jù)計算方法所需要的插值邊界層數(shù)進(jìn)行修正，將與計算單元相鄰的若干層非計算單元轉(zhuǎn)變?yōu)椴逯祮卧纯?圖2展示了應(yīng)用上述優(yōu)化步驟后的繞圓形網(wǎng)格與背景網(wǎng)格隱式挖洞結(jié)果，從圖可見重疊區(qū)域中質(zhì)量較差的的背景網(wǎng)格僅在邊界上根據(jù)求解器需要保留兩層插值單元，其余均作為非計算單元刪去，再加上刪去了落在物面內(nèi)部的計算單元，流場計算過程中計算效率會得到相應(yīng)的提升.

圖2 改進(jìn)后的隱式挖洞方法結(jié)果示例Fig.2 Illustration of enhanced implicit hole cutting method

2.2 提高貢獻(xiàn)單元搜索效率

對于傳統(tǒng)的挖洞方法，在洞邊界確定后只需要為所有插值單元搜索貢獻(xiàn)單元即可，而對于隱式挖洞方法，需要對所有位于重疊區(qū)域的網(wǎng)格單元執(zhí)行搜索貢獻(xiàn)單元的流程，再與可能的貢獻(xiàn)單元比較網(wǎng)格質(zhì)量以確定網(wǎng)格單元的類型.因?yàn)樗阉髫暙I(xiàn)單元的次數(shù)遠(yuǎn)高于傳統(tǒng)挖洞方法，隱式挖洞方法在自動化程度較高的同時也有著計算量較大的缺陷.

圖3 貢獻(xiàn)單元搜索流程示例Fig.3 Illustration of donor cells searching process

為了提高貢獻(xiàn)單元搜索效率，首先根據(jù)幾何體特征建立一個關(guān)于各個網(wǎng)格塊互相挖洞優(yōu)先級的N維矩陣，N為網(wǎng)格塊的數(shù)量.該矩陣中第n行第m列的元素表示第n個網(wǎng)格塊是否可能為第m個網(wǎng)格塊中網(wǎng)格單元提供貢獻(xiàn)單元，1表示可能，0表示不可能.若兩個網(wǎng)格塊互相沒有重疊區(qū)域，該值即為0，若兩個網(wǎng)格塊雖然互相重疊但第m個網(wǎng)格塊的網(wǎng)格質(zhì)量顯然優(yōu)于第n個網(wǎng)格塊，例如以網(wǎng)格體積與物面距乘積作為質(zhì)量判斷標(biāo)準(zhǔn)時，機(jī)翼端部(Cap)網(wǎng)格質(zhì)量明顯優(yōu)于同區(qū)域的繞機(jī)翼網(wǎng)格，則該值也取為0.反之若第n個網(wǎng)格塊存在為第m個網(wǎng)格塊中網(wǎng)格單元提供貢獻(xiàn)單元的可能性，該值即取為1.搜索貢獻(xiàn)單元之前先根據(jù)該矩陣的取值進(jìn)行篩選，只在有可能存在貢獻(xiàn)單元的網(wǎng)格塊中進(jìn)行貢獻(xiàn)單元搜索，其余網(wǎng)格塊直接跳過，這樣可以減少許多不必要的計算量.

貢獻(xiàn)單元的搜索方法則是基于模板跳躍法[11]發(fā)展而來.如圖3所示，對需要搜索貢獻(xiàn)單元的網(wǎng)格P，從某一初始單元開始，計算網(wǎng)格P中心點(diǎn)在該網(wǎng)格單元中的三線性插值信息，再根據(jù)等參變換獲得的局部坐標(biāo)決定初始單元需要向哪個方向移動以及移動的步數(shù)，如此不斷循環(huán)，直到找到能包圍網(wǎng)格P中心點(diǎn)的貢獻(xiàn)單元為止，該貢獻(xiàn)單元同時也作為下一網(wǎng)格搜索時的初始單元.

3 計算方法

本文中的算例均采用格心格式的有限體積法求解三維雷諾平均Navier-Stokes控制方程，其在絕對坐標(biāo)系下的積分形式為：

(1)

式中：W為守恒變量；Hc和Hv為對流通量和粘性通量.其表達(dá)式為：

(2)

(2) 式中：

(3)

根據(jù)理想氣體的性質(zhì)，再引入以下3個方程封閉原方程組：

(4)

(2)～(4)式中：ρ為流體密度；u,v,w為流體速度在3個坐標(biāo)方向上的分量；p為壓強(qiáng)；E為單位體積的總能；H為單位體積的總焓；k為熱傳導(dǎo)系數(shù)；T和τ為溫度和粘性應(yīng)力；V為流體速度矢量；γ為等壓比熱與等容比熱的比值；R為理想氣體常數(shù).

該方程的空間離散采用簡單高效的Jameson中心格式，加入人工耗散項(xiàng)抑制數(shù)值振蕩，時間離散采用魯棒性較強(qiáng)的隱式LU-SSOR方法完成.湍流模型采用計算量適中的Spalart-Allmaras方程模型.物面邊界為無滑移絕熱邊界條件，遠(yuǎn)場邊界采用Riemann不變量處理.采用當(dāng)?shù)貢r間步長、隱式殘值光順等措施加速收斂.

4 算 例

4.1 DLR-F4翼身組合體模型

DLR-F4翼身組合體是廣泛應(yīng)用于現(xiàn)代運(yùn)輸飛機(jī)跨音速氣動特性研究的經(jīng)典模型.該模型有著豐富的風(fēng)洞實(shí)驗(yàn)數(shù)據(jù)，被AIAA組織的第一屆阻力預(yù)測工作會議選作標(biāo)準(zhǔn)研究模型.計算模型為半模，共劃分為6個網(wǎng)格塊，總網(wǎng)格單元數(shù)量約為5.1×106個.計算條件為： 來流馬赫數(shù)Ma=0.75，雷諾數(shù)Re=3.0×106，迎角α=1.0°.

應(yīng)用本文發(fā)展的隱式挖洞方法得到的重疊網(wǎng)格如圖4所示，可見機(jī)翼網(wǎng)格與背景網(wǎng)格之間的插值邊界被有效地推離了機(jī)翼表面.圖5對隱式挖洞方法與傳統(tǒng)方法進(jìn)行了對比，比較兩套重疊網(wǎng)格進(jìn)行流場計算后得到的機(jī)翼展向截面壓力等值線可以發(fā)現(xiàn)，由于隱式挖洞方法將插值邊界推離了流場梯度變化較大的物面以及插值邊界處網(wǎng)格匹配性更好，該方法得到的壓力等值線過渡更光滑.

圖4 應(yīng)用隱式挖洞方法得到的重疊網(wǎng)格Fig.4 Illustration of the overset grids after implicit hole cutting method

圖5 不同挖洞方法得到的壓力等值線對比Fig.5 Comparison of pressure contours between different hole cutting methods

為檢驗(yàn)計算結(jié)果的精確性，我們將本文計算結(jié)果與風(fēng)洞實(shí)驗(yàn)數(shù)據(jù)及其他學(xué)者的計算結(jié)果[12]進(jìn)行了對比.圖6為機(jī)翼展向典型站位上的壓力系數(shù)(CP)分布曲線，圖7為升力系數(shù)曲線(CL)和阻力系數(shù)(CD)曲線，圖中y/b為等y截面的y坐標(biāo)與機(jī)翼展長b的比值，表明了該截面與機(jī)翼的相對位置關(guān)系，x/c為x坐標(biāo)與該截面上的機(jī)翼弦長c的比值，表明了該點(diǎn)在機(jī)翼弦向的相對位置.從圖可見本文計算結(jié)果與實(shí)驗(yàn)數(shù)據(jù)及其他計算結(jié)果吻合較好，本文發(fā)展的隱式挖洞方法在該算例中得到了成功應(yīng)用.

圖6 DLR-F4翼身組合體模型典型站位壓力系數(shù)分布曲線Fig.6 CP on typical stations of DLR-F4 wing-body configuration

圖7 DLR-F4翼身組合體模型升力系數(shù)曲線和阻力系數(shù)曲線Fig.7 Lift and drag coefficients curves of DLR-F4 wing-body configuration

4.2 NASA高升力梯形機(jī)翼模型

NASA高升力梯形機(jī)翼(NASA TrapWing)模型是專門用于驗(yàn)證機(jī)翼高升力構(gòu)型計算能力的標(biāo)準(zhǔn)模型，其幾何構(gòu)型與流動機(jī)理較為復(fù)雜，有著豐富的風(fēng)洞實(shí)驗(yàn)數(shù)據(jù)，被AIAA組織的第一屆高升力預(yù)測會議選作標(biāo)準(zhǔn)研究模型.本文中使用的的是全展襟翼構(gòu)型，計算模型為半模，共劃分為20個網(wǎng)格塊，總網(wǎng)格單元數(shù)量約為6.5×106個.計算條件為： 來流馬赫數(shù)Ma=0.2，雷諾數(shù)Re=4.63×106，迎角α=13.0°.

該高升力構(gòu)型主翼與前緣縫翼、后緣襟翼之間的縫道狹窄，流動現(xiàn)象復(fù)雜，應(yīng)用本文發(fā)展的隱式挖洞方法能較好地處理縫道內(nèi)各網(wǎng)格塊交界處區(qū)域的嵌套關(guān)系，確保插值區(qū)域的網(wǎng)格匹配良好.圖8展示了機(jī)翼展向截面y/b=0.65處的計算網(wǎng)格以及壓力等值線，可見縫道內(nèi)的壓力等值線基本能光滑連續(xù)地穿越插值區(qū)域.

圖8 機(jī)翼展向截面網(wǎng)格與壓力等值線圖Fig.8 Illustration of the overset grid and pressure contours at spanwise section

為了檢驗(yàn)計算結(jié)果的精確性，我們同樣將本文計算結(jié)果與NASA蘭利研究中心的風(fēng)洞實(shí)驗(yàn)數(shù)據(jù)及波音公司研究人員的計算結(jié)果[13]進(jìn)行了對比.圖9為機(jī)翼展向典型站位上的壓力系數(shù)分布曲線，圖10為升力系數(shù)曲線和阻力系數(shù)曲線.從圖可見本文計算值與實(shí)驗(yàn)數(shù)據(jù)及其他計算結(jié)果吻合較好，說明本文發(fā)展的隱式挖洞方法同樣在該高升力構(gòu)型計算中得到了成功應(yīng)用.

圖9 高升力梯形機(jī)翼模型典型站位壓力系數(shù)分布曲線Fig.9 CP on typical stations of TrapWing configuration

圖10 高升力梯形機(jī)翼模型升力系數(shù)和阻力系數(shù)曲線Fig.10 Lift and drag coefficients curves of TrapWing configuration

5 結(jié) 論

本文對重疊網(wǎng)格隱式挖洞方法進(jìn)行了改進(jìn)研究，并將改進(jìn)后的隱式挖洞方法應(yīng)用于復(fù)雜飛行器外形重疊網(wǎng)格生成與流場計算，得到了以下結(jié)論：

(1) 隱式挖洞方法的核心是網(wǎng)格單元質(zhì)量的選擇，以網(wǎng)格體積與物面距的乘積作為質(zhì)量判斷標(biāo)準(zhǔn)可以有效地將插值邊界推離流場變化梯度較大的物面，從而獲得更好的插值效果；

(2) 原始的隱式挖洞方法自動化程度高，但由于需要不斷重復(fù)貢獻(xiàn)單元搜索流程而計算量較大，通過適當(dāng)改進(jìn)后可以減少非必要的插值單元與計算單元數(shù)量并提高貢獻(xiàn)單元搜索效率，從而在保留原方法優(yōu)點(diǎn)的基礎(chǔ)上提高計算效率；

(3) 基于本文發(fā)展的改進(jìn)的隱式重疊網(wǎng)格方法，我們對翼身組合體和高升力構(gòu)型流場進(jìn)行了數(shù)值模擬，計算結(jié)果與實(shí)驗(yàn)值吻合較好，網(wǎng)格插值邊界處壓力等值線光滑，說明該方法在復(fù)雜外形網(wǎng)格生成工作中有較高的工程應(yīng)用價值.