崔鵬程，唐靜，李彬，馬明生，鄧有奇

中國(guó)空氣動(dòng)力研究與發(fā)展中心 計(jì)算空氣動(dòng)力研究所， 綿陽(yáng) 621000

在多體分離、艙門開啟、返回艙降落與螺旋槳滑流等非定常流動(dòng)模擬中，常采用重疊網(wǎng)格方法，將復(fù)雜的流動(dòng)區(qū)域劃分為多個(gè)簡(jiǎn)單的子區(qū)域，各子區(qū)域內(nèi)的網(wǎng)格獨(dú)立生成，數(shù)值計(jì)算在各個(gè)子網(wǎng)格上獨(dú)立進(jìn)行，各子網(wǎng)格之間存在重疊關(guān)系，流場(chǎng)信息通過插值在重疊區(qū)域進(jìn)行匹配和耦合[1]。重疊網(wǎng)格方法降低了網(wǎng)格生成難度，在多體相對(duì)運(yùn)動(dòng)數(shù)值計(jì)算方面具有很多優(yōu)勢(shì)，但各子網(wǎng)格間的數(shù)據(jù)插值很難保證守恒性[2]。

重疊網(wǎng)格間高精度插值是保證流場(chǎng)計(jì)算正確的重要基礎(chǔ)，目前主要有兩種插值方法：守恒插值和非守恒插值。非守恒插值方法實(shí)現(xiàn)簡(jiǎn)單，處理無(wú)間斷流動(dòng)時(shí)插值精度較高，對(duì)此國(guó)內(nèi)外很多學(xué)者開展了大量的研究工作。田書玲等[3-4]對(duì)物理空間的曲面單元進(jìn)行了標(biāo)準(zhǔn)化處理，實(shí)現(xiàn)了基于非結(jié)構(gòu)動(dòng)態(tài)網(wǎng)格的Lagrange插值方法；Marstin和Mcconnaughey[5]利用有限元思想構(gòu)造空間標(biāo)準(zhǔn)單元的線性插值函數(shù)，研究了雙線性、三線性重疊網(wǎng)格插值方法；康忠良等[6]基于線性重構(gòu)思想，發(fā)展了一種適用于任意單元的重疊網(wǎng)格插值策略，黃宇等[7]按單元尺度區(qū)分和擴(kuò)充貢獻(xiàn)單元模板，改善了網(wǎng)格尺寸相差較大時(shí)的插值精度；周乃春等[8]采用逆向距離權(quán)方法對(duì)重疊網(wǎng)格邊界的物理量進(jìn)行插值。

但非守恒插值方法對(duì)網(wǎng)格的相對(duì)尺寸依賴性較強(qiáng)，在網(wǎng)格尺寸相差較大的區(qū)域插值精度較低，且不能實(shí)現(xiàn)物理量在插值區(qū)域的守恒性，在非定常數(shù)值計(jì)算中可能會(huì)引起非物理擾動(dòng)，此外，非守恒方法即使計(jì)算網(wǎng)格質(zhì)量很好，也會(huì)使速度的計(jì)算值在重疊網(wǎng)格邊界之間產(chǎn)生錯(cuò)位[1,9]。

守恒插值方法精度較高，能實(shí)現(xiàn)物理量在插值區(qū)域的守恒性，對(duì)網(wǎng)格相對(duì)尺寸依賴性較小，但過程處理復(fù)雜，實(shí)現(xiàn)相對(duì)困難。常見的守恒插值方法有嵌入網(wǎng)格方法[10]、分段交界面方法[11]、ENO(Essentially Non Oscillatory)守恒重映方法[12]和總體守恒隱式插值方法[13]等。其中嵌入網(wǎng)格方法和分段交界面方法都是在幾何上將不重合邊界轉(zhuǎn)化為重合邊界，在三維情況下較難實(shí)現(xiàn)，ENO守恒重映方法需分解幾何單元，構(gòu)造幾何單元之間的面相交問題比較復(fù)雜，總體守恒隱式插值方法雖然計(jì)算精度較高，但需要求解附加方程組，給并行計(jì)算增添了難度。

Farrell和Maddison[14]發(fā)展了一種新的插值方法，提出“超網(wǎng)格”思想，構(gòu)造了二維網(wǎng)格的積分守恒插值，2011年文獻(xiàn)[15]進(jìn)一步發(fā)展了三維多面體的切割算法。徐春光等[16]發(fā)展了二維/三維混合網(wǎng)格下的單元相交算法，實(shí)現(xiàn)了二階精度的守恒插值方法。

本文在Farrell[14]、Menon[15]和徐春光[16]等的基礎(chǔ)上，基于“超網(wǎng)格”思想與單元相交算法，發(fā)展了一種適用于非定常數(shù)值計(jì)算的隱式并行混合重疊網(wǎng)格守恒插值方法，在每一個(gè)非定常時(shí)間步內(nèi)構(gòu)造局部超網(wǎng)格，利用局部超網(wǎng)格實(shí)現(xiàn)背景網(wǎng)格向目標(biāo)網(wǎng)格的守恒插值。經(jīng)過數(shù)值試驗(yàn)驗(yàn)證，對(duì)于二階分布的流場(chǎng)變量，該方法在重疊網(wǎng)格區(qū)域可以保證嚴(yán)格守恒插值，加快殘差收斂速度，提高流場(chǎng)變量的插值精度。

1 數(shù)值方法

1.1 控制方程

雷諾平均Navier-Stokes方程守恒形式的積分方程為

(1)

式中:W=[ρρuρvρwρE]T為守恒流場(chǎng)變量,ρ為密度,u、v、w為速度分量，E為能量；V為體積分;S為面積分；Ω為重疊區(qū)域;Hc為對(duì)流通量；Hv為黏性通量；n為流動(dòng)區(qū)域邊界上指向外側(cè)的單位法向量。

本文使用的流場(chǎng)解算器為自主研發(fā)的MFlow軟件[17]，對(duì)流通量和黏性通量的具體內(nèi)容可以參考文獻(xiàn)[18]，控制方程的離散采用格心格式的有限體積法，時(shí)間項(xiàng)采用LU-SGS(Lower-Upper Symmetric Gauss-Seidel)迭代[19]，湍流模型為一方程Spalart-Allmaras(S-A)模型[20]。黏性項(xiàng)采用中心格式，對(duì)流項(xiàng)采用Roe格式[21]，非定常計(jì)算過程采用多重網(wǎng)格加速收斂[22]。

1.2 積分型守恒插值條件

本文以兩套重疊子網(wǎng)格為例，簡(jiǎn)要說明重疊網(wǎng)格守恒插值方法。在守恒插值過程中要求守恒量的積分值保持恒定。

重疊子網(wǎng)格為ΓA、ΓB，重疊區(qū)域?yàn)棣?，設(shè)ΓA為背景網(wǎng)格，ΓB為插值網(wǎng)格，qA(x)為守恒變量在ΓA的分布函數(shù)，qB(x)為守恒變量在ΓB的分布函數(shù)，ΠAB為ΓA插值到ΓB的插值算子，如果滿足方程：

(2)

則稱qB(x)是qA(x)的守恒插值，其中qB(x)=qA(x)ΠAB。

2 超網(wǎng)格構(gòu)造方法

2.1 超網(wǎng)格定義

“超網(wǎng)格”概念最早由Farrell和Maddison提出[14]。對(duì)于任意兩套網(wǎng)格ΓA、ΓB，重疊區(qū)域?yàn)棣?，NA、NB分別為其網(wǎng)格節(jié)點(diǎn)，EA、EB分別為網(wǎng)格邊，K∈Γ為網(wǎng)格單元。由{ΓA,ΓB}組成的超網(wǎng)格ΓC滿足以下條件：

1)NC?NA∪NB。

2)V(KC∩K)∈{0,V(KC)}?KC∈ΓC,

K∈Γ,?！蕒ΓA,ΓB}。

即在ΓA、ΓB重疊區(qū)域內(nèi)，超網(wǎng)格的網(wǎng)格節(jié)點(diǎn)必須包含ΓA、ΓB所有的節(jié)點(diǎn)，且超網(wǎng)格ΓC的每一個(gè)網(wǎng)格單元KC與ΓA、ΓB網(wǎng)格單元的相交體體積為零或?yàn)镵C的體積，超網(wǎng)格的每個(gè)網(wǎng)格單元總是包含于ΓA或ΓB的網(wǎng)格單元。

2.2 局部超網(wǎng)格構(gòu)造

按照2.1節(jié)中超網(wǎng)格定義，任意兩套網(wǎng)格ΓA、ΓB的所有網(wǎng)格節(jié)點(diǎn)及邊可組成它們的超網(wǎng)格ΓC[15]。圖1為文獻(xiàn)[15]給出的超網(wǎng)格示意圖，圖中a～e為網(wǎng)格節(jié)點(diǎn)。將目標(biāo)網(wǎng)格ΓB的每一個(gè)網(wǎng)格單元與背景網(wǎng)格ΓA的網(wǎng)格單元求交，求交得到的網(wǎng)格單元可構(gòu)成ΓA、ΓB的超網(wǎng)格，母網(wǎng)格ΓA、ΓB的每一個(gè)網(wǎng)格單元總是可以拆分為一個(gè)或多個(gè)超網(wǎng)格單元的集合。

文獻(xiàn)[15]給出了全局超網(wǎng)格的構(gòu)造方法，即背景網(wǎng)格ΓA與目標(biāo)網(wǎng)格ΓB完全重合，所得超網(wǎng)格ΓC也與ΓA、ΓB完全重合。在重疊網(wǎng)格中各子網(wǎng)格區(qū)域不可能完全重合，且在非定常數(shù)值模擬中各子網(wǎng)格位置會(huì)相對(duì)變化。本文提出一種適用于重疊網(wǎng)格非定常數(shù)值模擬的局部超網(wǎng)格構(gòu)造方法，在每一個(gè)非定常時(shí)間步內(nèi)構(gòu)造局部超網(wǎng)格，利用局部超網(wǎng)格實(shí)現(xiàn)背景網(wǎng)格向目標(biāo)網(wǎng)格的守恒插值。局部超網(wǎng)格構(gòu)造方法如下：

1) 在重疊網(wǎng)格各子網(wǎng)格之間進(jìn)行挖洞操作，標(biāo)識(shí)出各子網(wǎng)內(nèi)不需要計(jì)算的網(wǎng)格、計(jì)算網(wǎng)格和插值網(wǎng)格等。

2) 用ADT(Alternating Digital Tree)搜索方法求出子網(wǎng)格每一個(gè)插值單元的所有貢獻(xiàn)單元[23]。普通的插值方法中一個(gè)插值單元只需用到一個(gè)貢獻(xiàn)單元信息，但在基于超網(wǎng)格的守恒插值方法中，一個(gè)插值單元需用到多個(gè)貢獻(xiàn)單元的信息。

構(gòu)造局部超網(wǎng)格需搜索出所有與插值單元相交的貢獻(xiàn)單元，圖2給出了搜索貢獻(xiàn)單元的示意圖，首先用ADT搜索方法快速定位插值單元的體心在背景網(wǎng)格的位置，體心所在的單元為1個(gè)貢獻(xiàn)單元，如圖2(a)所示；然后以此貢獻(xiàn)單元為基礎(chǔ)，檢索所有與此單元共點(diǎn)的背景網(wǎng)格，搜索出與插值網(wǎng)格相交的網(wǎng)格單元，這些單元也是貢獻(xiàn)單元；重復(fù)檢索與貢獻(xiàn)單元共點(diǎn)的網(wǎng)格單元，當(dāng)所有與貢獻(xiàn)單元共點(diǎn)的網(wǎng)格都不與插值單元相交時(shí)，終止搜索，如圖2(b)所示。

3) 求出重疊網(wǎng)格子網(wǎng)格的插值單元與其每一個(gè)貢獻(xiàn)單元的相交多面體，相交多面體可構(gòu)成局部超網(wǎng)格。

4) 以局部超網(wǎng)格為媒介，在每一步內(nèi)迭代計(jì)算中實(shí)現(xiàn)貢獻(xiàn)單元向插值單元的守恒插值。當(dāng)前時(shí)間步計(jì)算結(jié)束時(shí)，釋放超網(wǎng)格。下一個(gè)時(shí)間步開始計(jì)算時(shí)，根據(jù)當(dāng)前子網(wǎng)格的相對(duì)位置重新建立超網(wǎng)格。

2.3 網(wǎng)格求交算法

由2.2節(jié)可知，插值單元與貢獻(xiàn)單元求交是獲取超網(wǎng)格的關(guān)鍵，不同的網(wǎng)格單元有不同的網(wǎng)格求交算法。對(duì)凸多邊形或多面體而言，二維和三維網(wǎng)格求交方法已經(jīng)較為成熟[24]，本文采用快捷的Sutherland-Hodgman 逐邊切割方法[16,25]。

Sutherland-Hodgman方法的基本思想為用切割多邊形的各邊界線逐次對(duì)目標(biāo)多邊形進(jìn)行切割，從而得到一個(gè)新的多邊形。圖3為Suther-land-Hodgman方法示意圖，以兩個(gè)二維多邊形求交為例，三角形A1A2A3與B1B2B3求交。把A1A2A3當(dāng)做目標(biāo)三角形，B1B2B3當(dāng)做切割三角形。用B1B2B3的各條邊依次對(duì)A1A2A3進(jìn)行切割，得到相交多邊形B1CDEF。

三維情況下Sutherland-Hodgman方法的思想為用切割多面體的各個(gè)面逐次對(duì)目標(biāo)多面體進(jìn)行切割，從而形成一個(gè)新的多面體。在三維數(shù)值模擬中，常見的網(wǎng)格單元包括四面體、三棱柱、金字塔和六面體等單元。其中三棱柱、金字塔和六面體單元具有四邊形表面。在用常見的Pointwise、Icem等商業(yè)軟件生成網(wǎng)格的過程中，為了提高物面網(wǎng)格投影的精度和空間網(wǎng)格的光滑性，會(huì)對(duì)網(wǎng)格進(jìn)行扭轉(zhuǎn)或彎曲變形，從而導(dǎo)致四邊形的4個(gè)頂點(diǎn)不能嚴(yán)格地分布在同一平面內(nèi)，在多面體切割的過程中會(huì)產(chǎn)生歧義[26]。為了消除歧義，在多面體求交時(shí)，首先將三棱柱、金字塔和六面體單元分解為多個(gè)小四面體單元。通過小四面體相交體的組合來(lái)實(shí)現(xiàn)多面體求交。圖4所示為四面體網(wǎng)格求交算法，四面體P1與P2求交，將P2當(dāng)做目標(biāo)四面體，P1當(dāng)做切割四面體，用P1的各個(gè)面對(duì)P2進(jìn)行切割，保留切割面內(nèi)側(cè)的交點(diǎn)與交線，從而獲得P1、P2的相交體P3。相比平面多邊形而言，四面體求交得到的多面體構(gòu)造比較復(fù)雜，圖4中四面體P1、P2的相交體P3是具有10個(gè)頂點(diǎn)的7面體。

圖5所示為用網(wǎng)格求交方法構(gòu)造三維超網(wǎng)格的方法，具體步驟如下：

1) 在獲取重疊網(wǎng)格局部超網(wǎng)格時(shí)，先用ADT搜索方法求出子網(wǎng)格每一個(gè)插值單元的所有貢獻(xiàn)單元。在三維情況下，一個(gè)插值單元可能有多個(gè)對(duì)應(yīng)的貢獻(xiàn)單元，圖5(a)為一個(gè)插值單元與其相交的背景單元，這個(gè)插值單元有2個(gè)貢獻(xiàn)單元。

(3)

(4)

則插值單元PC與貢獻(xiàn)單元PD的相交體PI可表示為小四面體交集的并集：

(5)

其體積VI與體心XI分別為

(6)

3) 插值單元與貢獻(xiàn)單元的相交體可組成局部超網(wǎng)格，圖5(b)為一個(gè)插值單元與其超網(wǎng)格，每一個(gè)超網(wǎng)格單元都由插值單元與其貢獻(xiàn)單元相交獲得。

3 基于超網(wǎng)格的守恒插值方法

3.1 積分轉(zhuǎn)換

本文使用的流場(chǎng)計(jì)算器MFlow軟件控制方程的離散采用格心格式的有限體積法，流場(chǎng)數(shù)據(jù)儲(chǔ)存在網(wǎng)格體心處，網(wǎng)格K的體積為VK，其體心XK可以表示為

(7)

與式(6)類似，對(duì)于任意的多面體網(wǎng)格，其體積和體心都可以分解為四面體來(lái)計(jì)算，即

VK=∑Vi

(8)

式中：Vi和Xi分別為四面體的體積與體心。

一般流體計(jì)算中認(rèn)為流場(chǎng)變量呈二階分布，將流場(chǎng)變量Taylor展開：

(9)

可以證明，流場(chǎng)變量在K單元的平均值就等于體心值，流場(chǎng)變量在K單元的積分值等于體心值與體積的乘積,即

qK

(10)

3.2 守恒插值方法

以超網(wǎng)格為媒介，可以實(shí)現(xiàn)重疊網(wǎng)格子網(wǎng)格之間的守恒插值，具體方法如下。

式(2)是重疊網(wǎng)格守恒插值條件，將其在背景網(wǎng)格ΓA、插值網(wǎng)格ΓB上離散化，即

(11)

根據(jù)超網(wǎng)格的特點(diǎn)，插值網(wǎng)格ΓA的每一個(gè)網(wǎng)格單元KB總是可以拆分為一個(gè)或多個(gè)超網(wǎng)格單元的集合，設(shè)ΓC為其超網(wǎng)格，網(wǎng)格單元為KC，則式(11)可寫為[15]

(12)

這樣就把ΓA到ΓB的守恒插值轉(zhuǎn)換到ΓA到ΓC的守恒插值。根據(jù)超網(wǎng)格的特點(diǎn)，背景網(wǎng)格ΓA的每一個(gè)網(wǎng)格單元KA總是可以拆分為一個(gè)或多個(gè)超網(wǎng)格單元的集合。超網(wǎng)格單元KC的體積為VKC，則其體心XKC可表示為

(13)

設(shè)KC為貢獻(xiàn)單元KA所包含的一個(gè)超網(wǎng)格單元，可以利用二階Taylor重構(gòu)[6,27-28]將任意變量插值從貢獻(xiàn)單元體心XKA插值到超網(wǎng)格體心XKC處，即

(14)

設(shè)貢獻(xiàn)單元KA包含m個(gè)超網(wǎng)格，可以證明通過式(9)的插值方式，變量q在KA上的積分可寫為在KA所包含的m個(gè)超網(wǎng)格上的積分之和，即

(15)

綜合式(11)、式(12) 、式(15)，背景網(wǎng)格ΓA的流場(chǎng)變量可通過超網(wǎng)格ΓC插值到插值網(wǎng)格ΓB，即

(16)

式(15)整個(gè)插值過程對(duì)二階分布的流場(chǎng)變量嚴(yán)格守恒，具體實(shí)現(xiàn)過程總結(jié)如下：

1) 將流場(chǎng)變量從貢獻(xiàn)單元KA插值到它的超網(wǎng)格單元KC，即

(17)

2) 將超網(wǎng)格ΓC的流場(chǎng)變量采用體積權(quán)插值到插值單元KB，設(shè)KB包含n個(gè)超網(wǎng)格單元KC，則其流場(chǎng)變量q(XKB)可寫為

(18)

這樣以超網(wǎng)格為媒介，將流場(chǎng)變量從貢獻(xiàn)單元KA插值到插值單元KB。如果流場(chǎng)變量q(x)二階分布，則整個(gè)插值過程嚴(yán)格守恒，這種插值方法不會(huì)改變?cè)攘鲌?chǎng)的二階分布，實(shí)現(xiàn)了超網(wǎng)格單元和插值單元流場(chǎng)變量的精確守恒插值。

3.3 隱式并行插值方法

在重疊網(wǎng)格插值過程中，根據(jù)貢獻(xiàn)單元是否同時(shí)也是插值單元，插值方法可分為隱式插值和顯式插值。顯式插值方法需區(qū)分背景網(wǎng)格與插值網(wǎng)格，插值速度較慢，不適合復(fù)雜構(gòu)型的大規(guī)模重疊網(wǎng)格。本文發(fā)展了一種隱式并行插值方法，不嚴(yán)格區(qū)分背景網(wǎng)格與插值網(wǎng)格，如圖6所示，每一套網(wǎng)格既向其他網(wǎng)格發(fā)送數(shù)據(jù)，同時(shí)又從其他網(wǎng)格接收數(shù)據(jù)。

這種隱式并行的重疊網(wǎng)格插值方法插值速度較快，適用于多個(gè)相對(duì)運(yùn)動(dòng)的子網(wǎng)格間的流場(chǎng)數(shù)據(jù)插值，能處理復(fù)雜構(gòu)型的大規(guī)模重疊網(wǎng)格。

4 算例分析

本文首先采用不同類型的流場(chǎng)方程對(duì)所建立的插值方法的守恒性和插值精度進(jìn)行驗(yàn)證。然后將該方法用于CFD重疊網(wǎng)格標(biāo)模的數(shù)值模擬，通過跟試驗(yàn)數(shù)據(jù)對(duì)比，考核該方法的魯棒性和工程實(shí)用性。

4.1 守恒性驗(yàn)證

4.1.1 線性方程

假設(shè)物理量φ在網(wǎng)格區(qū)域線性分布，滿足方程：

φ(x,y,z)=2x+3y+z

(19)

如圖7所示，ΓA、ΓB分別為結(jié)構(gòu)與非結(jié)構(gòu)四面體網(wǎng)格，網(wǎng)格區(qū)域一致。在網(wǎng)格上設(shè)定物理量φ的分布滿足式(19)，然后ΓA、ΓB間進(jìn)行100次相互插值，統(tǒng)計(jì)φ在ΓB網(wǎng)格區(qū)域內(nèi)的積分值與每個(gè)單元的平均插值誤差絕對(duì)值。

圖8對(duì)比了本文插值方法與三線性插值方法、逆向距離權(quán)插值方法在守恒性和插值精度上的異同。圖8(a)為φ在ΓB網(wǎng)格區(qū)域內(nèi)的積分值，可以看出用本文守恒插值方法計(jì)算得出的積分值在插值過程中保持恒定，按照式(11)守恒插值的定義，該方法在計(jì)算過程中可保持物理量積分的嚴(yán)格守恒。而三線性插值方法和逆向距離權(quán)插值方法不具備守恒性，且隨著插值次數(shù)的增加，物理量積分值與理論值的誤差越來(lái)越大。

圖8(b)為φ在ΓB網(wǎng)格區(qū)域內(nèi)每個(gè)單元的平均插值誤差絕對(duì)值，可以看出與三線性插值方法和逆向距離權(quán)插值方法相比，本文所建立的插值方法插值精度更高。

4.1.2 三角函數(shù)

假設(shè)物理量φ在網(wǎng)格區(qū)域呈三角函數(shù)分布，滿足方程：

φ(x,y,z)=1+sin(2πx)sin(2πy)sin(2πz)

(20)

插值網(wǎng)格與計(jì)算條件同4.1.1小節(jié)。圖9對(duì)比了三角函數(shù)的插值過程中，本文插值方法與三線性插值方法、逆向距離權(quán)插值方法在守恒性和插值精度的異同?？梢钥闯觯疚陌l(fā)展的守恒插值方法對(duì)三角函數(shù)的插值計(jì)算具有較強(qiáng)的守恒性(見圖9(a))，且插值精度明顯高于三線性插值方法與逆向距離權(quán)插值方法(見圖9(b))。

4.2 定常數(shù)值計(jì)算

4.2.1 ONERA-M6機(jī)翼

將重疊網(wǎng)格數(shù)值計(jì)算結(jié)果與單套整體網(wǎng)格數(shù)值計(jì)算結(jié)果對(duì)比是一種常用的檢驗(yàn)重疊網(wǎng)格技術(shù)優(yōu)劣的方法[7,23]。為了考察本文重疊網(wǎng)格插值方法的插值精度，對(duì)ONERA-M6機(jī)翼分別采用非結(jié)構(gòu)重疊網(wǎng)格與單套整體網(wǎng)格進(jìn)行數(shù)值模擬，計(jì)算來(lái)流條件為:馬赫數(shù)Ma=0.84, 迎角α=3.06°, 來(lái)流溫度T∞=255.56 K, 雷諾數(shù)Re=1.172×107。

重疊網(wǎng)格由兩部分組成：背景網(wǎng)格和機(jī)翼網(wǎng)格，機(jī)翼表面采用各向異性網(wǎng)格對(duì)機(jī)翼前后緣進(jìn)行加密，機(jī)翼近壁面采用三棱柱網(wǎng)格以較好地模擬黏性，遠(yuǎn)壁面采用四面體、金字塔和六面體的混合網(wǎng)格，機(jī)翼網(wǎng)格量約為60.3萬(wàn)。背景網(wǎng)格全部采用四面體網(wǎng)格單元，網(wǎng)格量約為6.1萬(wàn)。兩套子網(wǎng)格的重疊效果如圖10所示，總網(wǎng)格量為66.4萬(wàn)。

單套全機(jī)網(wǎng)格的近壁面網(wǎng)格分布與重疊網(wǎng)格完全一致，空間網(wǎng)格尺度和網(wǎng)格數(shù)量與兩套重疊網(wǎng)格相當(dāng)，總網(wǎng)格量為54.8萬(wàn)。

如圖11(a)和圖11(b)所示，本文將不同重疊網(wǎng)格插值方法的壓力系數(shù)Cp計(jì)算結(jié)果與相同網(wǎng)格尺度單套整體網(wǎng)格的計(jì)算結(jié)果進(jìn)行對(duì)比，可以看出相比三線性插值方法與逆向距離權(quán)插值方法，本文所提出的插值方法計(jì)算得出的壓力系數(shù)和單套整體網(wǎng)格的計(jì)算結(jié)果更為符合,其中L為機(jī)翼展長(zhǎng)。由于本文建立的插值方法提高了重疊網(wǎng)格交界面的插值精度，使得激波附近的壓力系數(shù)計(jì)算更為精確。

從圖11(a)和圖11(b)可以發(fā)現(xiàn)CFD計(jì)算結(jié)果與風(fēng)洞試驗(yàn)結(jié)果符合良好，前緣吸力峰和激波位置與風(fēng)洞試驗(yàn)數(shù)據(jù)較為一致。但仔細(xì)對(duì)比發(fā)現(xiàn)，整體網(wǎng)格與守恒插值方法計(jì)算得出的Cp與試驗(yàn)值略有偏差。ONERA-M6機(jī)翼是一個(gè)跨聲速機(jī)翼，如圖11(c)和圖11(d)所示，大量的數(shù)值計(jì)算[29-33]表明，對(duì)于ONERA-M6機(jī)翼跨聲速狀態(tài)，全局加密網(wǎng)格、S-A模型、DES和DDES方法的數(shù)值計(jì)算結(jié)果都無(wú)法與風(fēng)洞試驗(yàn)完全符合，這是數(shù)值計(jì)算和風(fēng)洞試驗(yàn)的固有偏差[29-30]。而三線性與距離權(quán)插值方法在激波后部與風(fēng)洞試驗(yàn)完全符合，這是完全不合理的，也恰好說明這兩種重疊網(wǎng)格插值方法數(shù)值誤差較大。本文建立的守恒插值方法計(jì)算結(jié)果與整體網(wǎng)格、文獻(xiàn)[29-33]一致，說明該方法插值精度較高，產(chǎn)生的數(shù)值誤差較小。

圖12給出了用不同重疊網(wǎng)格插值方法計(jì)算的殘差收斂曲線，可以看出由于本文守恒插值方法在重疊網(wǎng)格區(qū)域插值精度較高，引入誤差較小，使得收斂速度較快。本文插值方法的收斂速度分別比逆向距離權(quán)方法和三線性方法快20.0%和10.8%。

4.2.2 NACA0012翼型激波干擾

取2個(gè)并列錯(cuò)位的NACA0012翼型進(jìn)行超聲速擾流計(jì)算，來(lái)流條件為:Ma=2,α=0°,Re=1.812×107。

在此計(jì)算條件下，由于翼型并列錯(cuò)位，2個(gè)翼型之間會(huì)產(chǎn)生強(qiáng)烈的激波干擾，流場(chǎng)結(jié)構(gòu)比較復(fù)雜，本文采用重疊網(wǎng)格對(duì)此算例進(jìn)行數(shù)值模擬，如圖13所示,重疊網(wǎng)格由兩部分組成：上翼型網(wǎng)格和下翼型網(wǎng)格。

為了考核插值方法在重疊網(wǎng)格邊界的插值精度，人為地將計(jì)算網(wǎng)格分布得較為稀疏，在重疊網(wǎng)格邊界,網(wǎng)格的尺度相差較大。

圖14(a)～圖14(c)給出了不同插值方法計(jì)算得到的馬赫數(shù)等值線分布，白色、黑色分別為上、下翼型網(wǎng)格的馬赫數(shù)等值線，圖14(d)為細(xì)網(wǎng)格計(jì)算結(jié)果，3種插值方法都正確模擬了流場(chǎng)結(jié)構(gòu)?？梢钥闯?，在上下翼型重疊網(wǎng)格交界面附近，用三線性和逆向距離權(quán)插值方法計(jì)算得到的馬赫數(shù)等值線有明顯的錯(cuò)位，而守恒插值方法得到的馬赫數(shù)等值線較為光滑。在下面翼型的尾流區(qū)域，人為設(shè)置了兩套子網(wǎng)格的尺度，使得在重疊網(wǎng)格邊界兩套網(wǎng)格的尺度相差較大，在這種情況下，三線性和逆向距離權(quán)插值方法計(jì)算得到的等值線出現(xiàn)了明顯的間斷，而守恒插值方法得到的馬赫數(shù)等值線仍然較為連續(xù)。此外，2個(gè)NACA0012翼型并列錯(cuò)位，產(chǎn)生了復(fù)雜的激波干擾流場(chǎng)，相比三線性和逆向距離權(quán)方法，用守恒插值方法計(jì)算得到的激波分布更為精細(xì)、光滑。

4.3 三維機(jī)翼分離投放模型

選擇機(jī)翼/掛架/帶舵外掛物(Wing/Pylon/Finned-Store，WPFS)模型，對(duì)本文所建立的重疊網(wǎng)格守恒插值方法進(jìn)行驗(yàn)證。該模型是美國(guó)進(jìn)行分離投放CFD數(shù)值驗(yàn)證時(shí)使用的一個(gè)機(jī)翼和導(dǎo)彈的簡(jiǎn)化模型，有詳實(shí)可靠的風(fēng)洞試驗(yàn)數(shù)據(jù)，因此常被當(dāng)做驗(yàn)證重疊網(wǎng)格分離投放計(jì)算的標(biāo)準(zhǔn)算例。外掛物共有4片尾翼，呈X形分布。在該試驗(yàn)中，還使用了彈射力，使外掛物在分離初期保持抬頭趨勢(shì)，避免其尾部和掛架碰撞，具體數(shù)據(jù)詳見文獻(xiàn)[34]。

計(jì)算狀態(tài)取和文獻(xiàn)[34]一致：來(lái)流馬赫數(shù)Ma=0.95，初始迎角α=0°，初始側(cè)滑角為0°，本文采用非定常雷諾平均Navier-Stokes方程、S-A湍流方程、6自由度(DOF)方程來(lái)模擬計(jì)算該模型的外掛物分離歷程。采用6套重疊子網(wǎng)格進(jìn)行計(jì)算，分別為：用于模擬計(jì)算的包裹機(jī)翼和掛架的背景網(wǎng)格、包裹帶掛物的子網(wǎng)格以及包裹每一個(gè)舵的子網(wǎng)格，總網(wǎng)格量為553萬(wàn)。WPFS模型子網(wǎng)格的重疊效果如圖15(a)所示，表面網(wǎng)格如圖15(b)所示。

圖16(a)是初始狀態(tài)在執(zhí)行“挖洞”操作后得到的重疊網(wǎng)格計(jì)算單元的剖面圖，物理量通過相互重疊的網(wǎng)格部分進(jìn)行插值傳遞。圖16(b)是時(shí)間t=0.3 s時(shí)各個(gè)子網(wǎng)格相互重疊示意圖，可以看到在計(jì)算過程中，子網(wǎng)格跟隨物體位置的變化而變化。圖17給出了t=0 s時(shí)刻與t=0.3 s時(shí)刻的實(shí)時(shí)壓力分布云圖，計(jì)算得出的流場(chǎng)結(jié)構(gòu)清晰合理。

將計(jì)算得到的外掛物分離動(dòng)態(tài)特性同試驗(yàn)數(shù)據(jù)進(jìn)行比較，圖18(a)和圖18(b)分別為外掛物的質(zhì)心位移隨時(shí)間的變化量以及對(duì)應(yīng)時(shí)刻的姿態(tài)角。可以看出，計(jì)算得到的軌跡同試驗(yàn)數(shù)據(jù)吻合較好，證明本文建立的重疊網(wǎng)格守恒插值方法準(zhǔn)確有效，魯棒性強(qiáng)，工程實(shí)用。

5 結(jié) 論

1) 針對(duì)二階線性分布的流場(chǎng)變量，采用局部超網(wǎng)格方法在重疊網(wǎng)格區(qū)域可實(shí)現(xiàn)精確守恒插值，對(duì)非線性分布的流場(chǎng)變量，本文插值方法守恒性更高。

2) 結(jié)合局部超網(wǎng)格與ADT搜索方法，重疊網(wǎng)格插值單元能夠快速搜索到多個(gè)貢獻(xiàn)單元，合理地?cái)U(kuò)大了貢獻(xiàn)單元模板，提高了插值精度和流場(chǎng)光滑性，加快了殘差收斂速度。

3) 本文重疊網(wǎng)格守恒插值方法適用于非結(jié)構(gòu)網(wǎng)格與結(jié)構(gòu)網(wǎng)格的各類型網(wǎng)格單元，工程實(shí)用、魯棒性強(qiáng)。插值時(shí)采用體積加權(quán)方法，當(dāng)子網(wǎng)格尺度相差較大時(shí)，仍然能夠保持較好的流場(chǎng)光滑性與連續(xù)性。

本文發(fā)展了一種重疊網(wǎng)格守恒插值新方法，下一步將開展重疊網(wǎng)格洞面優(yōu)化技術(shù)與嵌入網(wǎng)格方法，進(jìn)一步提高重疊子網(wǎng)格尺度相差較大時(shí)流場(chǎng)的光滑性。

參 考 文 獻(xiàn)

[1] 李鵬, 高振勛, 蔣崇文. 重疊網(wǎng)格方法的研究進(jìn)展[J]. 力學(xué)與實(shí)踐, 2014, 36(5):551-565.

LI P, GAO Z X, JIANG C W. The progress of the overlapping grid techniques[J]. Mechanics in Engineering, 2014, 36(5):551-565(in Chinese).

[2] SHIH T I P. Overset grids: Fundamentals and practical issues: AIAA-2002-3259[R]. Reston, VA: AIAA, 2002.

[3] 田書玲, 伍貽兆, 夏健. 用動(dòng)態(tài)非結(jié)構(gòu)重疊網(wǎng)格法模擬三維多體相對(duì)運(yùn)動(dòng)繞流[J]. 航空學(xué)報(bào), 2007, 28(1):46-51.

TIAN S L, WU Y Z, XIA J. Simulation of flows past multi-body in relative motion with dynamic unstructured overset grid method[J]. Acta Aeronautica et Astronautica Sinica, 2007, 28(1): 46-51 (in Chinese).

[4] 伍貽兆, 田書玲, 夏健. 基于非結(jié)構(gòu)動(dòng)網(wǎng)格的非定常流數(shù)值模擬方法[J]. 航空學(xué)報(bào), 2011, 32(1):15-26.

WU Y Z, TIAN S L, XIA J. Unstructured grid methods for unsteady flow simulation[J]. Acta Aeronautica et Astronautica Sinica, 2011, 32(1): 15-26 (in Chinese).

[5] MARSTIN C, MCCONNAUGHEY H. Computational problems on composite grids: AIAA-1984-1611[R]. Reston, VA: AIAA, 1984.

[6] KANG Z L, YAN C, YU J, et al. A fast and reliable overset unstructured grids approach[J]. Acta Mechanica Sinica, 2013, 29(2): 149-157.

[7] 黃宇, 閻超, 王文, 等. 混合重疊網(wǎng)格插值方法的改進(jìn)及應(yīng)用[J]. 北京航空航天大學(xué)學(xué)報(bào), 2017, 43(2):285-292.

HUANG Y, YAN C, WANG W, et al. An improved interpolation method for hybrid overset grid and its application[J]. Journal of Beijing University of Aeronautics and Astronautics, 2017, 43(2): 285-292 (in Chinese).

[8] 周乃春, 李彬, 鄭鳴, 等. 帶控制率導(dǎo)彈投放數(shù)值模擬[J]. 空氣動(dòng)力學(xué)學(xué)報(bào), 2013, 31(3): 107-115.

ZHOU N C, LI B, ZHENG M, et al. Missile separation simulation with control laws[J]. Acta Aerodynamica Sinica, 2013, 31(3): 107-115 (in Chinese).

[9] FARRELL P E, PIGGOTT M D, PAIN C C, et al. Conservative interpolation between unstructured meshes via supermesh construction[J]. Computer Methods in Applied Mechanics and Engineering, 2009, 198(33): 2632-2642.

[10] ZHENG Y, LIOU M S. Progress in the three-dimensional DRAGON grid scheme: AIAA-2001-2540[R]. Reston, VA: AIAA, 2001.

[11] XU K, SUN G, CAI J Y. On interface conservative attributions for computations of the complex flows of high-lift system based on chimera tech-nique[C]∥Computer and Automation Engineering (ICCAE), 2010 The 2nd International Conference on IEEE. Piscataway, NJ: IEEE Press, 2010: 279-283.

[12] 張宇飛, 陳海昕, 符松. 基于高階守恒重映對(duì)窗口嵌入技術(shù)的改進(jìn)[J]. 計(jì)算物理, 2011, 28(2): 167-173.

ZHANG Y F, CHEN H X, FU S. Improvement on window embedment technology with high order conservative remapping[J]. Chinese Journal of Computational Physics, 2011, 28(2): 167-173(in Chinese).

[13] ZHAO X, GUAN H W, YANG Z, et al. An implicit and globally conservative unstructured chimera grid method: AIAA-2011-0777[R]. Reston, VA: AIAA, 2011.

[14] FARRELL P E, MADDISON J R. Conservative interpolation between volume meshes by local Galerkin projection[J]. Computer Methods in Applied Mechanics and Engineering, 2011, 200(1): 89-100.

[15] MENON S, SCHMIDT D P. Conservative interpolation on unstructured polyhedral meshes: An extension of the supermesh approach to cell-centered finite-volume variables[J]. Computer Methods in Applied Mechanics and Engineering, 2011, 200(41): 2797-2804.

[16] 徐春光, 董海波, 劉君. 基于單元相交的混合網(wǎng)格精確守恒插值方法[J]. 爆炸與沖擊, 2016, 36(3): 305-312.

XU C G, DONG H B, LIU J. An accurate conservative interpolation method for mixed grid based on the intersection of grid cells[J]. Explosion and Shock Waves, 2016, 36(3): 305-312 (in Chinese).

[17] 崔鵬程, 鄧有奇, 唐靜, 等．基于伴隨方程的網(wǎng)格自適應(yīng)及誤差修正技術(shù)[J]．航空學(xué)報(bào), 2016,37(10): 2992-3002.

CUI P C, DENG Y Q, TANG J, et al. Adjoint-based grid adaptation and error correction[J]. Acta Aeronautica et Astronautica Sinica, 2016, 37(10): 2992-3002 (in Chinese).

[18] BLAZEK J. Computational fluid dynamics: principles and applications[M]. 3rd ed. Oxford: Elsevier, 2015:75-120.

[19] KIM J S, KWON O J. Improvement on block LU-SGS scheme for unstructured mesh Navier-Stokes computations: AIAA-2002-1061[R]. Reston, VA: AIAA, 2002.

[20] SPALART S R. A one-equation turbulence model for aerodynamic flows: AIAA-1992-0439[R]. Reston, VA: AIAA, 1992.

[21] ROE P L. Approximate Riemann solvers, parameter vectors, and difference schemes[J]. Journal of Computational Physics, 1981, 43(2): 357-372.

[22] 李彬, 唐靜, 鄧有奇, 等. 并行的多重網(wǎng)格方法在離散伴隨優(yōu)化中的應(yīng)用[J]. 航空學(xué)報(bào), 2014, 35(8): 2091-2101.

LI B, TANG J, DENG Y Q, et al. Application of parallel multigrid algorithm to discrete adjoint optimization[J]. Acta Aeronautica et Astronautica Sinica, 2014, 35(8): 2091-2101 (in Chinese).

[23] BEATRICE R, JAY S. Robust and scalable overset grid assembly for partitioned unstructured meshes: AIAA-2013-0797[R]. Reston, VA: AIAA, 2013.

[24] THEOHARIS T, IAN P. Two parallel methods for polygon clipping[J]. Computer Graphics Forum, 2010, 8(2):107-114.

[25] MARTINEZ F, RUEDA A J, FEITO F R. A new algorithm for computing Boolean operations on polygons[J]. Computers & Geosciences, 2009, 35(6): 1177-1185.

[26] 唐靜, 鄧有奇, 馬明生, 等. 飛翼氣動(dòng)優(yōu)化中參數(shù)化和網(wǎng)格變形技術(shù)研究[J]. 航空學(xué)報(bào), 2015, 36(5): 1480-1490.

TANG J, DENG Y Q, MA M S, et al. Parametrization and grid deformation techniques for fly-wing shape optimization[J]. Acta Aeronautica et Astronautica Sinica, 2015, 36(5): 1480-1490 (in Chinese).

[27] ALAUZET F, MEHRENBERGER M. P1-conservative solution interpolation on unstructured triangular meshes[J]. International Journal for Numerical Methods in Engineering, 2010, 84(13): 1552-1588.

[28] GARIMELLA R, KUCHARIK M, SHASHKOV M. An efficient linearity and bound preserving conservative interpolation (remapping) on polyhedral meshes[J]. Computers & Fluids, 2007, 36(2): 224-237.

[29] MAYEUR J, DUMONT A, DESTARAC D, et al. RANS simulations on TMR test cases and M6 wing with the ONERA elsA flow solver: AIAA-2015-1745[R]. Reston, VA: AIAA, 2015.

[30] DURRANI N, QIN N. Comparison of RANS, DES and DDES results for ONERA M6 wing at transonic flow speed using an in-house parallel code: AIAA-2011-0190[R]. Reston, VA: AIAA, 2011.

[31] DA SILVA R G, AZEVEDO J L F, BASSO E. Simulation of ONERA M6 wing flows for assessment of turbulence modeling capabilities: AIAA-2016-0549[R]. Reston, VA: AIAA, 2016.

[32] SAXENA S K, NAIR M T. Implementation and testing of Spalart-Allmaras model in a multi-block code: AIAA-2002-0835[R]. Reston, VA: AIAA, 2002.

[33] LUO H, BAUM J D, L?HNER R. A Hermite WENO-based limiter for discontinuous Galerkin method on unstructured grids[J]. Journal of Computational Physics, 2007, 225(1): 686-713.

[34] HEIM E R. CFD wing/pylon/finned store mutual interference wind tunnel experiment: AEDC-TSR-91-P4[R]. New York: AEDC, 1991.