江維青，楊愛明

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

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

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

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

1 隱式挖洞方法簡(jiǎn)介

隱式挖洞方法的核心是一個(gè)基于網(wǎng)格單元質(zhì)量做選擇的過程，除了給定網(wǎng)格本身信息之外不需要過多的人工干預(yù).其過程可以分為兩步.(1) 首先是在網(wǎng)格重疊區(qū)域計(jì)算所有網(wǎng)格單元的質(zhì)量，該參數(shù)表征了一個(gè)網(wǎng)格單元適合成為貢獻(xiàn)單元的程度，可以是網(wǎng)格單元的體積、物面距、長(zhǎng)寬比等.(2) 然后為所有落在重疊區(qū)域內(nèi)的網(wǎng)格搜索貢獻(xiàn)單元，若存在質(zhì)量更好的貢獻(xiàn)單元，則該單元成為插值單元，從質(zhì)量最好的貢獻(xiàn)單元處通過三線性插值獲取流場(chǎng)信息；若貢獻(xiàn)單元不存在或貢獻(xiàn)單元的質(zhì)量均比該單元差，則該單元成為計(jì)算單元正常參與流場(chǎng)計(jì)算，這樣就避免了傳統(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)格單元一般來說體積與物面距均較小，因此會(huì)被判別為高質(zhì)量的計(jì)算單元，使得插值邊界自然地遠(yuǎn)離了物面處流場(chǎng)變化梯度較大的區(qū)域，從而確保了更好的插值精度.

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

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

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

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

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

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

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

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

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

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

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

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

3 計(jì)算方法

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

(1)

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

(2)

(2) 式中：

(3)

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

(4)

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

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

4 算 例

4.1 DLR-F4翼身組合體模型

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

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

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

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

為檢驗(yàn)計(jì)算結(jié)果的精確性，我們將本文計(jì)算結(jié)果與風(fēng)洞實(shí)驗(yàn)數(shù)據(jù)及其他學(xué)者的計(jì)算結(jié)果[12]進(jìn)行了對(duì)比.圖6為機(jī)翼展向典型站位上的壓力系數(shù)(CP)分布曲線，圖7為升力系數(shù)曲線(CL)和阻力系數(shù)(CD)曲線，圖中y/b為等y截面的y坐標(biāo)與機(jī)翼展長(zhǎng)b的比值，表明了該截面與機(jī)翼的相對(duì)位置關(guān)系，x/c為x坐標(biāo)與該截面上的機(jī)翼弦長(zhǎng)c的比值，表明了該點(diǎn)在機(jī)翼弦向的相對(duì)位置.從圖可見本文計(jì)算結(jié)果與實(shí)驗(yàn)數(shù)據(jù)及其他計(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)型計(jì)算能力的標(biāo)準(zhǔn)模型，其幾何構(gòu)型與流動(dòng)機(jī)理較為復(fù)雜，有著豐富的風(fēng)洞實(shí)驗(yàn)數(shù)據(jù)，被AIAA組織的第一屆高升力預(yù)測(cè)會(huì)議選作標(biāo)準(zhǔn)研究模型.本文中使用的的是全展襟翼構(gòu)型，計(jì)算模型為半模，共劃分為20個(gè)網(wǎng)格塊，總網(wǎng)格單元數(shù)量約為6.5×106個(gè).計(jì)算條件為： 來流馬赫數(shù)Ma=0.2，雷諾數(shù)Re=4.63×106，迎角α=13.0°.

該高升力構(gòu)型主翼與前緣縫翼、后緣襟翼之間的縫道狹窄，流動(dòng)現(xiàn)象復(fù)雜，應(yīng)用本文發(fā)展的隱式挖洞方法能較好地處理縫道內(nèi)各網(wǎng)格塊交界處區(qū)域的嵌套關(guān)系，確保插值區(qū)域的網(wǎng)格匹配良好.圖8展示了機(jī)翼展向截面y/b=0.65處的計(jì)算網(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)計(jì)算結(jié)果的精確性，我們同樣將本文計(jì)算結(jié)果與NASA蘭利研究中心的風(fēng)洞實(shí)驗(yàn)數(shù)據(jù)及波音公司研究人員的計(jì)算結(jié)果[13]進(jìn)行了對(duì)比.圖9為機(jī)翼展向典型站位上的壓力系數(shù)分布曲線，圖10為升力系數(shù)曲線和阻力系數(shù)曲線.從圖可見本文計(jì)算值與實(shí)驗(yàn)數(shù)據(jù)及其他計(jì)算結(jié)果吻合較好，說明本文發(fā)展的隱式挖洞方法同樣在該高升力構(gòu)型計(jì)算中得到了成功應(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é) 論

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

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

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

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