基于單元相交的混合網(wǎng)格精確守恒插值方法*

徐春光，董海波，劉 君

(大連理工大學(xué)航空航天學(xué)院，遼寧 大連 116024)

基于單元相交的混合網(wǎng)格精確守恒插值方法*

徐春光，董海波，劉 君

(大連理工大學(xué)航空航天學(xué)院，遼寧 大連 116024)

基于網(wǎng)格切割思想，發(fā)展了二維/三維混合網(wǎng)格條件下的單元相交算法，可精確計(jì)算任意兩個(gè)多邊形/多面體的交集。在此基礎(chǔ)上，實(shí)現(xiàn)了基于單元相交(CIB/DC)的精確守恒插值算法。二維和三維驗(yàn)證算例表明，該方法能夠保證插值過程中計(jì)算域內(nèi)物理量的嚴(yán)格守恒，且具有比常規(guī)二階插值更高的精度。

流體力學(xué);守恒插值;單元相交;混合網(wǎng)格;局部超網(wǎng)格

在包含運(yùn)動(dòng)邊界的非定常流動(dòng)模擬中，流場物理空間隨時(shí)間變化，邊界運(yùn)動(dòng)導(dǎo)致網(wǎng)格也隨時(shí)間變化，典型應(yīng)用包括多體分離、機(jī)翼抖振顫振、艙門啟閉等[1]。CFD計(jì)算一般采用運(yùn)動(dòng)嵌套網(wǎng)格、笛卡爾自適應(yīng)網(wǎng)格和變形動(dòng)網(wǎng)格等技術(shù)處理網(wǎng)格變化。在采用上述技術(shù)處理網(wǎng)格變化時(shí)，均涉及流場信息的插值問題[1-2]，即利用舊網(wǎng)格上的流動(dòng)信息插值獲得重構(gòu)區(qū)域流動(dòng)參數(shù)。此外，在爆轟流場模擬中，爆轟計(jì)算網(wǎng)格與沖擊波傳播計(jì)算網(wǎng)格差兩個(gè)量級(jí)[3-4]，為提高計(jì)算效率，往往需要將物理量由爆轟計(jì)算時(shí)使用的密網(wǎng)格插值傳遞到?jīng)_擊波傳播模擬中的疏網(wǎng)格。

插值過程一般應(yīng)具備守恒性和單調(diào)性等基本特性。對(duì)爆轟計(jì)算等守恒性要求高的問題，一般的插值不能保證計(jì)算區(qū)域內(nèi)物理量的守恒，會(huì)降低計(jì)算精度，且無法捕捉激波的正確位置[5]，需要發(fā)展守恒插值方法。守恒插值方法主要可分為兩類，基于面通量(SFB/DC， simplified face-based donor-cell)的方法和基于單元相交(CIB/DC, cell-intersection-based donor-cell)的方法[6]。

SFB/DC方法[6]效率較高，但要求新、舊網(wǎng)格具有相同的拓?fù)浣Y(jié)構(gòu)，不適用于網(wǎng)格重構(gòu)等拓?fù)浒l(fā)生變化的情況。CIB/DC方法思想簡單直接，對(duì)計(jì)算網(wǎng)格的類型、結(jié)構(gòu)不進(jìn)行任何假設(shè)，適用范圍很廣，并且具有保號(hào)性，插值過程中不會(huì)產(chǎn)生新的極值[6]。但該方法推廣至三維情況時(shí)，由于三維網(wǎng)格單元重疊區(qū)域切割計(jì)算非常復(fù)雜，給這種方法的應(yīng)用帶來了很大困難。為了實(shí)現(xiàn)嚴(yán)格意義下的三維CIB/DC算法，P.E.Farrell等[7]進(jìn)行了一系列卓有成效的工作，在2011年進(jìn)一步發(fā)展為“局部超網(wǎng)格”(Local Supermesh)方法，成功實(shí)現(xiàn)了非結(jié)構(gòu)網(wǎng)格下的三維CIB/DC算法[8]。

本文中在P.E.Farrell等和S.Menon等[9]工作的基礎(chǔ)上，構(gòu)造了應(yīng)用于二維/三維混合網(wǎng)格的CIB/DC算法，實(shí)現(xiàn)了基于精確網(wǎng)格切割的積分守恒插值。驗(yàn)證結(jié)果表明，該方法能夠保證插值過程中計(jì)算域內(nèi)物理量的嚴(yán)格守恒，具有比常規(guī)二階插值更高的精度。

1 守恒插值基本思想

網(wǎng)格Tb是計(jì)算域Ω的另一種形式的劃分，劃分形式與要求同Ta，則同樣有Tb上的物理量分布φb(x)。如果對(duì)于Ω的任意子區(qū)域D，滿足如下條件：

(1)

則稱φb(x)是φa(x)的守恒插值結(jié)果。

對(duì)于一般的情形，要使式(1)滿足，必須有Tb≡Ta，φb(x)≡φa(x)，但這就失去了在兩套網(wǎng)格間進(jìn)行插值的意義。為此，對(duì)D進(jìn)行限制，不能是Ω的任意子區(qū)域，而必須是網(wǎng)格Tb中的某個(gè)網(wǎng)格：

(2)

上式就是一般意義上的守恒插值關(guān)系式。

對(duì)式(2)，按照基于單元相交的思想，對(duì)D進(jìn)行分解：

(3)

式(3)就是基于單元相交守恒插值方法(CIB/DC)的基本依據(jù)。從式(2)到式(3)未引入任何近似。因此，在給定分布φ(x)的情況下，CIB/DC方法是嚴(yán)格守恒的。

從式(3)還可以看出，如果已知網(wǎng)格Ta上的物理量分布φa(x)，如何計(jì)算網(wǎng)格單元間的相交區(qū)域Vij便成為CIB/DC方法的關(guān)鍵。

2 “超網(wǎng)格”基本思想

圖1為文獻(xiàn)[8]給出的超網(wǎng)格示意圖，超網(wǎng)格是由Ta和Tb所有結(jié)點(diǎn)以及它們邊界交點(diǎn)構(gòu)成的網(wǎng)格，具有如下特性：

(1)每一個(gè)超網(wǎng)格不再和Ta或Tb中任何網(wǎng)格單元相交，即它總可以在Ta中找到一個(gè)網(wǎng)格包含或等于它，也總可以在Tb中找到一個(gè)網(wǎng)格包含或等于它；

(2)Ta的每一個(gè)網(wǎng)格總可以明確表示為一個(gè)或多個(gè)超網(wǎng)格的集合，Tb的每一個(gè)網(wǎng)格總可以明確表示為一個(gè)或多個(gè)超網(wǎng)格的集合。

圖1 超網(wǎng)格示意圖Fig.1 Schematic of the supermesh

3 單元相交算法

由上節(jié)可知，獲得超網(wǎng)格是實(shí)現(xiàn)精確守恒插值的關(guān)鍵。下面給出計(jì)算平面凸多邊形和三維凸多面體交集的具體算法。

3.1 平面凸多邊形相交算法

如圖2所示，將目標(biāo)三角形A1A2A3視為切割多邊形，使用切割三角形B1B2B3對(duì)目標(biāo)三角形進(jìn)行切割，每次切割時(shí)只保留目標(biāo)多邊形中與切割多邊形在切割線同側(cè)的部分。圖中依次使用△B1B2B3的3條邊B1B2、B2B3、B3B1切割△A1A2A3，切割后保留的多邊形依次為CA3A2G、CDEA2G、CDEFB1，多邊形CDEFB1即為△A1A2A3和△B1B2B3的交集。

圖2 二維網(wǎng)格切割示意圖Fig.2 Schematic of the two-dimensional grid cutting

圖3 直線切割平面凸多邊形Fig.3 Schematic of a straight line cutting plane convex polygon

在使用直線去切割平面凸多邊形時(shí)，以圖3為例，切割直線經(jīng)過點(diǎn)O，其法向單位矢量為n，切割直線將平面劃分為兩部分，定義法向n所指的半平面為正半平面，反方向的半平面為負(fù)半平面。如果切割之后目標(biāo)多邊形僅保留正半平面的部分，那么切割操作等效于求目標(biāo)多邊形和正半平面的交集。定義平面上任意一點(diǎn)P到切割直線的距離為：

d=xOP·n=(xP-xO)·n

式中：xP和xP分別為點(diǎn)P和點(diǎn)O坐標(biāo)。

根據(jù)上述定義，正半平面的點(diǎn)到切割直線的距離均為正數(shù)，負(fù)半平面的點(diǎn)到切割直線的距離均為負(fù)數(shù)，切割直線上的點(diǎn)到切割直線的距離為零。如果凸多邊形各頂點(diǎn)到切割直線的距離均為非負(fù)數(shù)(圖3(a))，則切割結(jié)果即為原多邊形；如果凸多邊形各頂點(diǎn)到切割直線的距離均為非正數(shù)，則切割結(jié)果為空集；除上述兩種情況外，凸多邊形與切割直線必然有兩個(gè)交點(diǎn)(圖3(b))。具體如下：

(1)設(shè)目標(biāo)多邊形為N，其頂點(diǎn)依次為N1,N2,…,NN；切割直線為l，切割結(jié)果為R；

(2)計(jì)算N的各個(gè)頂點(diǎn)到l的距離，記頂點(diǎn)Ni到l的距離為di；

(3)初始化一個(gè)空的有序點(diǎn)集P；

(4)從頂點(diǎn)N1起，依次對(duì)多邊形N的邊N1N2,N2N3,…,NNN1進(jìn)行判定。

對(duì)線段NiNj的判定過程為：如果didj<0，則線段與切割直線有新的交點(diǎn)，記交點(diǎn)為Nnew；依次將Ni、Nnew、Nj中距離為非負(fù)數(shù)的頂點(diǎn)加入點(diǎn)集P；如果didj≥0，則依次將Ni、Nj中距離為非負(fù)數(shù)的頂點(diǎn)加入點(diǎn)集P；

(5)完成所有邊的判定后，如果P為空集，則R為空集；若P為非空集，則依次連接P中的點(diǎn)形成的多邊形就是所求的切割結(jié)果R。

3.2 三維凸多面體相交算法

在三維流體計(jì)算中，常用的網(wǎng)格包括四面體、四棱錐、三棱柱、六面體等多種類型。其中，四棱錐、三棱柱、六面體這3種網(wǎng)格都具有四邊形表面，網(wǎng)格生成過程中不能嚴(yán)格保證四邊形表面的四個(gè)頂點(diǎn)位于同一平面內(nèi)，會(huì)給多面體表面的定義帶來歧義。

為了簡化算法的復(fù)雜性，同時(shí)消除歧義，在計(jì)算兩個(gè)任意凸多面體的交集時(shí)，首先將凸多面體分解為多個(gè)小四面體，如圖4所示，其中點(diǎn)N6為點(diǎn)N2、N3、N4、N5的重心；然后通過計(jì)算這些小四面體的交集，來獲得多面體的交集。計(jì)算多面體Pa和Pb的交集的具體步驟為：

圖4 三維網(wǎng)格分解示意圖Fig.4 Schematic of the three dimensional mesh decomposition

(3)多面體Pa和Pb交集可表示為：

交集的體積和重心分別是：

4 相交單元的搜索算法

相對(duì)于通常的流場插值方法，守恒插值的計(jì)算效率較低，主要有兩個(gè)原因：(1)計(jì)算兩個(gè)單元交集的算法較為復(fù)雜；(2)對(duì)同一個(gè)新網(wǎng)格單元，普通插值方法只需用到一個(gè)舊網(wǎng)格單元(稱為貢獻(xiàn)單元)的信息，而守恒插值需要用到多個(gè)貢獻(xiàn)單元的信息。前文已對(duì)單元相交算法進(jìn)行了詳細(xì)的描述，本節(jié)將介紹守恒插值中多個(gè)貢獻(xiàn)單元的搜索算法。

在貢獻(xiàn)單元的搜索算法中，首先需要快速定位出新網(wǎng)格單元在舊網(wǎng)格中的位置。對(duì)此，可采用四叉樹[11]或八叉樹[12](三維)方法加速搜尋，如圖5(a)所示，完全填充網(wǎng)格為舊網(wǎng)格，其上的非填充網(wǎng)格為新網(wǎng)格中的某個(gè)單元，通過四叉樹方法，可快速定位出新網(wǎng)格單元的格心所在的舊網(wǎng)格單元(圖中填充單元)，這是新網(wǎng)格單元的第1個(gè)貢獻(xiàn)單元；然后，逐個(gè)搜索各個(gè)貢獻(xiàn)單元的直接相鄰單元，如果相鄰單元與新網(wǎng)格單元相交，則將該單元加入貢獻(xiàn)單元的列表，如圖5(b)中淺色填充單元所示；當(dāng)所有貢獻(xiàn)單元的相鄰單元(圖5(c)與(b)中不同的填充單元)都不與新網(wǎng)格單元相交時(shí)，搜索結(jié)束。

圖5 相交單元的搜索方法Fig.5 Search method of intersection elements

5 守恒插值驗(yàn)證算例

圖6 數(shù)值計(jì)算建模圖Fig.6 Geometry of the specimen used for simulation

圖7 計(jì)算結(jié)果的誤差比較Fig.7 Different errors of computational results

圖8 計(jì)算結(jié)果的積分比較Fig.8 Different integrals of computational results

設(shè)計(jì)了3個(gè)算例來驗(yàn)證本文發(fā)展的守恒插值方法。第1個(gè)算例主要考核相同拓?fù)浣Y(jié)構(gòu)下網(wǎng)格點(diǎn)位置變化的影響；第2個(gè)算例主要考核網(wǎng)格形狀的影響；第3個(gè)算例主要考核三維情況下的應(yīng)用情況。

算例1

考慮如下的測試函數(shù)：

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

計(jì)算區(qū)域?yàn)閇0,1]×[0,1]，計(jì)算網(wǎng)格為一組拓?fù)浣Y(jié)構(gòu)相同的非結(jié)構(gòu)三角形網(wǎng)格，分別表示為M0、M1、…、MN，其中M0如圖6所示。網(wǎng)格MN中的第(i,j)個(gè)網(wǎng)格點(diǎn)的坐標(biāo)由下式給出：

式中：α(t)=0.5sin(4πt),ξ=(i-1)/(imax-1),η=(j-1)/(jmax-1)t=n/T,其中i=1,2,…,imax;j=1,2,…,jmax;n=0,1,…,N。

測試過程中，在網(wǎng)格M0上給定測試函數(shù)的準(zhǔn)確分布，然后依次插值到網(wǎng)格M1、M2、…、MN上，共進(jìn)行N次插值，最后比較網(wǎng)格MN上的函數(shù)分布與準(zhǔn)確值之間的誤差。在本文計(jì)算中，取imax=jmax=65，N=200，T=20。

圖7給出了計(jì)算誤差隨插值次數(shù)的變化情況，其中縱坐標(biāo)代表計(jì)算區(qū)域內(nèi)各單元函數(shù)分布插值后計(jì)算值與理論值之差的1范數(shù)。兩種插值方法中，插值誤差都隨插值次數(shù)的增加而增大，但守恒型插值方法的誤差始終較小，不到常規(guī)二階插值誤差的一半。

圖8為兩種方法的物理量積分值比較，其中縱坐標(biāo)代表每個(gè)單元計(jì)算值總和與理論值的比，從圖中明顯可以看出，守恒插值保證了計(jì)算域內(nèi)物理量的守恒。

算例2

考慮如下的測試函數(shù)：

φ(r)=2+cos(πr/L)

計(jì)算區(qū)域是邊長為1的正方形，L為正方形的對(duì)角線長度，r是與正方形中心的距離。測試中使用兩種不同的網(wǎng)格，網(wǎng)格1是均勻分布的結(jié)構(gòu)網(wǎng)格，網(wǎng)格2是非結(jié)構(gòu)網(wǎng)格，并與網(wǎng)格1具有相同的邊界網(wǎng)格點(diǎn)分布，如圖9所示。測試中，首先在網(wǎng)格1上給出測試函數(shù)的準(zhǔn)確分布，然后將其插值到網(wǎng)格2上，再插值回網(wǎng)格1，如此反復(fù)，進(jìn)行200次插值后，將物理量分布與初始的準(zhǔn)確分布進(jìn)行比較，統(tǒng)計(jì)插值過程中的誤差。分別采用常規(guī)二階插值和守恒型二階插值進(jìn)行測試，并使用不同疏密的網(wǎng)格以分析插值方法的精度。測試中使用的4組網(wǎng)格的網(wǎng)格量如表1所示。

表1 計(jì)算網(wǎng)格參數(shù)Table 1 Parameters of the computational grid

圖9 數(shù)值計(jì)算建模圖Fig.9 Geometry of the specimen used for simulation

圖10是計(jì)算誤差與網(wǎng)格步長之間的關(guān)系，其中橫坐標(biāo)是表1中4個(gè)不同結(jié)構(gòu)網(wǎng)格計(jì)算單元的邊長，縱坐標(biāo)代表計(jì)算區(qū)域內(nèi)各單元函數(shù)分布插值后計(jì)算值與理論值之差的1范數(shù)。可以看出在相同的網(wǎng)格步長下，守恒型插值的誤差更小，而且誤差與網(wǎng)格步長之間的變化斜率更大，說明其具有更高階的精度。計(jì)算結(jié)果表明，常規(guī)二階插值精度約為1.4階，而守恒型二階插值精度約為1.87階。表2是不同算例的守恒誤差比較，守恒插值方法的守恒誤差接近機(jī)器零，明顯優(yōu)于普通二階插值。

表2 不同算例的守恒誤差Table 2 Conservation errors indifferent numerical cases

圖10 插值方法的精度比較Fig.10 Comparison of precision with different interpolation methods

算例3

考慮如下的測試函數(shù)：

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

計(jì)算區(qū)域是邊長為1的正方體。每次測試使用兩套網(wǎng)格，首先在網(wǎng)格1上給出測試函數(shù)的準(zhǔn)確分布，然后將其插值到網(wǎng)格2上，再插值回網(wǎng)格1，如此反復(fù)，進(jìn)行10次插值后，將物理量分布與初始的準(zhǔn)確分布進(jìn)行比較，統(tǒng)計(jì)插值過程中的誤差。分別采用普通二階插值和守恒型二階插值進(jìn)行測試，并使用不同疏密的網(wǎng)格以分析插值方法的精度。測試中使用的4組網(wǎng)格的網(wǎng)格量如表3所示，其中網(wǎng)格步長由下式得到：

式中：NE1和NE2分別為兩套網(wǎng)格的單元個(gè)數(shù)，ND=3為網(wǎng)格的維數(shù)。

圖11是測試函數(shù)在三維空間中的分布云圖，圖12是計(jì)算誤差與網(wǎng)格步長之間的關(guān)系，其中橫坐標(biāo)是表3中dx的值，縱坐標(biāo)代表計(jì)算區(qū)域內(nèi)各單元函數(shù)分布插值后計(jì)算值與理論值之差的1范數(shù)。在三維條件下，守恒型插值依然具有更小的計(jì)算誤差和更高的計(jì)算精度。

表3 計(jì)算網(wǎng)格參數(shù)Table 3 Parameters of the computational grid

圖11 測試函數(shù)分布云圖Fig.11 Nephogram of different test functions

圖12 插值方法的精度比較Fig.12 Comparison of precision with different interpolation methods

6 結(jié) 論

采用在新舊網(wǎng)格相交單元剖分構(gòu)建出來的超級(jí)網(wǎng)格作為網(wǎng)格之間信息傳遞的基礎(chǔ)，采用四叉樹(二維)或八叉樹(三維)方法搜尋算法提高計(jì)算效率，基于以上方法提出了網(wǎng)格之間流場物理量的守恒插值算法。該算法理論上可以實(shí)現(xiàn)網(wǎng)格之間物理量的精確守恒，二維和三維算例也表明產(chǎn)生的誤差接近機(jī)器零，明顯優(yōu)于常用的二階插值。本文中提出的方法可以解決爆轟流場模擬中網(wǎng)格尺度和物理量相差較大情況下的插值難題。

[1] 郭正,劉君,瞿章華.非結(jié)構(gòu)動(dòng)網(wǎng)格在三維可動(dòng)邊界問題中的應(yīng)用[J].力學(xué)學(xué)報(bào),2003,35(2):140-146. Guo Zheng, Liu Jun, Qu Zhanghua. Dynamic unstructured grid method with applications to 3d unsteady flows involving moving boundaries[J]. Acta Mechanica Sinica, 2003,35(2):140-146.

[2] Zeeuw D D, Powell K G. An adaptively refined Cartesian mesh solver for the Euler equations[J]. Journal of Computational Physics, 1993, 104(1),104:56-68.

[3] 白曉征.包含運(yùn)動(dòng)界面的爆炸流場數(shù)值模擬方法及其應(yīng)用[D].長沙:國防科技大學(xué),2009.

[4] 徐春光,白曉征,劉瑜,等.爆炸近區(qū)流場的數(shù)值模擬方法研究[J].兵工學(xué)報(bào),2012,33(5):565-573. Xu Chunguang, Bai Xiaozheng, Liu Yu, et al. Research on numerical simulation method of near-field flows in air blast[J]. Acta Armmamentarii, 2012,33(5):565-573.

[5] P?rt-Enander E, Sj?green B. Conservative and non-conservative interpolation between overlapping grids for finite volume solutions of hyperbolic problems[J]. Computers and Fluids, 1994,23(3):551-574.

[6] Margolin L G, Shashkov M. Second-order sign-preserving conservative interpolation (remapping) on general grids[J]. Journal of Computational Physics, 2002,184(1):266-298.

[7] 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/34/35/36):2632-2642.

[8] 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/2/3/4):89-100.

[9] 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/44):2797-2804.

[10] 郭正.包含運(yùn)動(dòng)邊界的多體非定常流場數(shù)值模擬方法研究[D].長沙:國防科技大學(xué),2002.

[11] 劉君,白曉征,郭正.非結(jié)構(gòu)動(dòng)網(wǎng)格計(jì)算方法----及其在包含運(yùn)動(dòng)界面的流場模擬中的應(yīng)用[M].長沙:國防科技大學(xué)出版社,2009.

(責(zé)任編輯 曾月蓉)

An accurate conservative interpolation method for the mixed grid based on the intersection of grid cells

Xu Chunguang, Dong Haibo, Liu Jun

(SchoolofAeronauticsandAstronautics,DalianUniversityofTechnology,Dalian116024,Liaoning,China)

A method is presented for conservatively transferring cell-centered physical quantities from one mesh to another with second-order accuracy, which is well-suited for finite-volume methods that rely on cell-centered variables. The proposed methodology implements the cell-intersection algorithm of 2D/3D mixed grid based on the local supermesh idea, and is able to accurately calculate the intersection of any two polygons or polyhedrons, thus providing a basis for accurate conservative interpolation. The accuracy and the efficacy of the new method are demonstrated with multiple 2D and 3D numerical experiments. The test results show that this method ensures strict conservation of physical quantities in the interpolation process, and achieves a higher accuracy than that achieved by conventional second-order interpolation methods.

fluid mechanics; conservative interpolation; cell-intersection; mixed grid; local supermesh

10.11883/1001-1455(2016)03-0305-08

2014-11-10； < class="emphasis_bold">修回日期：2015-03-30

2015-03-30

國家自然科學(xué)基金項(xiàng)目(11272074);遼寧省自然科學(xué)基金項(xiàng)目(201202033)

徐春光(1977— )，男，博士,高級(jí)工程師;

劉 君,liujun65@dlut.edu.cn。

O354 <國標(biāo)學(xué)科代碼：13025 class="emphasis_bold"> 國標(biāo)學(xué)科代碼：13025 文獻(xiàn)標(biāo)志碼：A國標(biāo)學(xué)科代碼：13025

A