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

徐春光，董海波，劉 君

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

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

徐春光，董海波，劉 君

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

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

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

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

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

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

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

1 守恒插值基本思想

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

(1)

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

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

(2)

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

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

(3)

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

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

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

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

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

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

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

3 單元相交算法

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

3.1 平面凸多邊形相交算法

如圖2所示，將目標三角形A1A2A3視為切割多邊形，使用切割三角形B1B2B3對目標三角形進行切割，每次切割時只保留目標多邊形中與切割多邊形在切割線同側(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

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

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

式中：xP和xP分別為點P和點O坐標。

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

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

(2)計算N的各個頂點到l的距離，記頂點Ni到l的距離為di；

(3)初始化一個空的有序點集P；

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

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

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

3.2 三維凸多面體相交算法

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

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

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

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

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

4 相交單元的搜索算法

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

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

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

5 守恒插值驗證算例

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

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

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

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

算例1

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

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

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

式中：α(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ù)的準確分布，然后依次插值到網(wǎng)格M1、M2、…、MN上，共進行N次插值，最后比較網(wǎng)格MN上的函數(shù)分布與準確值之間的誤差。在本文計算中，取imax=jmax=65，N=200，T=20。

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

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

算例2

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

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

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

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

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

圖10是計算誤差與網(wǎng)格步長之間的關(guān)系，其中橫坐標是表1中4個不同結(jié)構(gòu)網(wǎng)格計算單元的邊長，縱坐標代表計算區(qū)域內(nèi)各單元函數(shù)分布插值后計算值與理論值之差的1范數(shù)?？梢钥闯鲈谙嗤木W(wǎng)格步長下，守恒型插值的誤差更小，而且誤差與網(wǎng)格步長之間的變化斜率更大，說明其具有更高階的精度。計算結(jié)果表明，常規(guī)二階插值精度約為1.4階，而守恒型二階插值精度約為1.87階。表2是不同算例的守恒誤差比較，守恒插值方法的守恒誤差接近機器零，明顯優(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)

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

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

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

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

[1] 郭正,劉君,瞿章華.非結(jié)構(gòu)動網(wǎng)格在三維可動邊界問題中的應用[J].力學學報,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] 白曉征.包含運動界面的爆炸流場數(shù)值模擬方法及其應用[D].長沙:國防科技大學,2009.

[4] 徐春光,白曉征,劉瑜,等.爆炸近區(qū)流場的數(shù)值模擬方法研究[J].兵工學報,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] 郭正.包含運動邊界的多體非定常流場數(shù)值模擬方法研究[D].長沙:國防科技大學,2002.

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

(責任編輯 曾月蓉)

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

國家自然科學基金項目(11272074);遼寧省自然科學基金項目(201202033)

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

劉 君,liujun65@dlut.edu.cn。

O354 <國標學科代碼：13025 class="emphasis_bold"> 國標學科代碼：13025 文獻標志碼：A國標學科代碼：13025

A