二維任意凸四邊形網(wǎng)格的MoF界面重構(gòu)解析算法

王成， 孔澤宇， 李濤， 賈祖朋

(1.北京理工大學(xué) 爆炸科學(xué)與技術(shù)國家重點(diǎn)實(shí)驗(yàn)室，北京 100081；2.北京應(yīng)用物理與計(jì)算數(shù)學(xué)研究所，北京 100088)

流體矩(moment of fluid，MoF)方法是一種分段線性網(wǎng)格重構(gòu)方法，主要用于多物質(zhì)流體數(shù)值模擬中進(jìn)行界面重構(gòu)，是目前VOF方法群中最為先進(jìn)的一種. 該方法可以在單元尺度上進(jìn)行界面重構(gòu)，只依賴于單元內(nèi)部數(shù)據(jù)，精度有著顯著的優(yōu)勢[1]. 在前人的成果中，MoF界面重構(gòu)方法應(yīng)用在很多情況中. 例如，在自適應(yīng)網(wǎng)格問題中利用MoF方法的形心數(shù)據(jù)判斷高曲率區(qū)域進(jìn)行界面細(xì)分[2]；在任意拉格朗日-歐拉模型中采用MoF方法的形心數(shù)據(jù)，對非結(jié)構(gòu)網(wǎng)格上的二維可壓縮流體進(jìn)行界面重構(gòu)[3]，以及直接采用MoF方法進(jìn)行多物質(zhì)條件下的ALE界面重構(gòu)[4]. 另外的一些針對方法的研究中，在軸對稱坐標(biāo)系下進(jìn)行了MoF方法界面重構(gòu)[5]，進(jìn)一步在柱坐標(biāo)的ALE方法的MoF方法的界面重構(gòu)實(shí)現(xiàn)[6]；通過使用對稱重構(gòu)方法，使多物質(zhì)MoF過程中實(shí)現(xiàn)兩材料的矩量對稱計(jì)算，同時(shí)提高計(jì)算精度和計(jì)算效率[7]；使用絲狀網(wǎng)格界面捕捉方法，在二維薄網(wǎng)格難以進(jìn)行常規(guī)界面重建處理的情況下實(shí)現(xiàn)了MoF方法，進(jìn)而實(shí)現(xiàn)了網(wǎng)格尖端處的界面重構(gòu)精度提升[8]. 針對MoF的界面求解精度和求解效率，一些研究中提出了求解一階偏導(dǎo)數(shù)以提高迭代效率和魯棒性的方法[9]，證明對目標(biāo)函數(shù)求解解析一階偏導(dǎo)數(shù)對混合單元問題界面重構(gòu)問題的解決效率有著較大的提升. Antoine Lemoine等[10]探究了正交網(wǎng)格的界面解析重構(gòu)，在笛卡爾坐標(biāo)下甚至所有矩形網(wǎng)格中實(shí)現(xiàn)了解析解的界面重構(gòu)工作，結(jié)合王成等[11]的重映改進(jìn)方法，可以大幅提升界面重構(gòu)的計(jì)算精度和計(jì)算效率并應(yīng)用在MMALE方法中. 本文對這一思路進(jìn)行了更進(jìn)一步的嘗試，在更為廣泛的凸四邊形網(wǎng)格中實(shí)現(xiàn)解析解的界面重構(gòu).

1 傳統(tǒng)MoF方法中的界面重構(gòu)

1.1 MoF方法簡介

在傳統(tǒng)的MoF方法處理多物質(zhì)網(wǎng)格界面重構(gòu)問題時(shí)，一般采用的是數(shù)值計(jì)算迭代尋找曲面的最優(yōu)近似的方法. MoF方法重構(gòu)界面主要依靠多邊形網(wǎng)格Ω中的各物質(zhì)子集ω的零階矩M0(ω)和一階矩M1(ω)來計(jì)算[1]，

(1)

在一般情況下，前兩階矩的等價(jià)量有著實(shí)際的物理意義，即體積分?jǐn)?shù)μ(ω)和形心xc(ω):

(2)

因此，MoF方法一般使用物質(zhì)的體積分?jǐn)?shù)和形心來進(jìn)行界面重構(gòu).

MoF界面重構(gòu)的核心問題在于找到實(shí)際的物質(zhì)子集ω*的一個(gè)最優(yōu)多邊形子網(wǎng)格近似ω*，這個(gè)近似子網(wǎng)格的邊界Γ*是實(shí)際界面Γ的一個(gè)最優(yōu)近似，如圖1所示.

由于MoF方法是體積守恒的，因此更進(jìn)一步地說，確定這個(gè)近似子網(wǎng)格ω*是一個(gè)求最小值的問題，即找到一個(gè)多邊形近似子網(wǎng)格ω*，使得實(shí)際物質(zhì)子集ω的形心xc(ω)和這個(gè)子網(wǎng)格的形心xc(ω*)距離最小.

(3)

圖1 實(shí)際物質(zhì)子集和重構(gòu)的多邊形近似子網(wǎng)格

1.2 傳統(tǒng)MoF界面重構(gòu)方法

在傳統(tǒng)的MoF方法中，求解最小值問題通過迭代圖1(b)顯示的界面法向量角度來實(shí)現(xiàn)，在保證體積分?jǐn)?shù)μ(ω)不變的條件下，每選定一個(gè)法向量n，就可以對應(yīng)確定一個(gè)線性子網(wǎng)格界面Γ*，進(jìn)而求解出多邊形子網(wǎng)格的形心xc(ω*)，進(jìn)一步找出一個(gè)法向量n0使得這個(gè)情況下的多邊形形心xc(ω*)與實(shí)際形心xc(ω)之間的距離最近. 求解過程中法向量n迭代旋轉(zhuǎn)角度決定了迭代求解的精度. 為了簡化這一過程，賈祖朋等[1]采用了Brent迭代、限定旋轉(zhuǎn)范圍等各種方法進(jìn)行迭代，在此不再一一說明.

1.3 解析解求解凸四邊形MoF界面重構(gòu)方法

對任意的一個(gè)凸四邊形網(wǎng)格，進(jìn)行線性界面重構(gòu)后形成的兩個(gè)子網(wǎng)格形狀可能性是有限的，針對不同形狀的每一個(gè)子網(wǎng)格都可以求解出這個(gè)子網(wǎng)格情況下的一個(gè)物質(zhì)形心(x*,y*)，(x*,y*)隨著界面的變動而運(yùn)動，因?yàn)榻缑娴倪\(yùn)動是連續(xù)的，所以(x*,y*)的運(yùn)動也是連續(xù)的，即所有的子網(wǎng)格形心(x*,y*)可以組成一條曲線F(x*,y*).

當(dāng)物質(zhì)的實(shí)際子網(wǎng)格ω形心(x,y)已知時(shí)，MoF方法的求最小值問題簡化為計(jì)算曲線F(x*,y*)上距離形心(x,y)最近的一個(gè)點(diǎn). 當(dāng)曲線復(fù)雜度不高時(shí)，可以直接計(jì)算出該點(diǎn).

2 對任意凸四邊形網(wǎng)格的切分

2.1 對任意四邊形網(wǎng)格的旋轉(zhuǎn)簡化

x′y′坐標(biāo)系內(nèi)任意的一個(gè)四邊形A′B′C′D′，首先將其平移至xy坐標(biāo)系中的ABCD位置，令點(diǎn)A與xy坐標(biāo)系的原點(diǎn)重合；然后旋轉(zhuǎn)坐標(biāo)系，使得向量AB在xy坐標(biāo)系中為(1,0). 旋轉(zhuǎn)方程為：

(4)

坐標(biāo)變換后存在4個(gè)互相獨(dú)立的未知變量xC、yC、xD、yD，其中點(diǎn)D在AD邊上，可以用AD邊斜率倒數(shù)k1和點(diǎn)D的y坐標(biāo)yD來表示點(diǎn)D位置：

xD=k1yD

(5)

同樣地可以使用BC邊斜率倒數(shù)k2和點(diǎn)C的y坐標(biāo)yC來表示點(diǎn)C的位置

xC=k2yC+1

(6)

從而可以利用4個(gè)獨(dú)立變量率k1、yC、k2、yD確定四邊形的位置.

2.2 凸四邊形網(wǎng)格的限制條件

(7)

而若k1≤k2，則四邊形網(wǎng)格必然存在. 綜合兩種情況，得到四邊形網(wǎng)格的限制條件為

(8)

2.3 凸四邊形網(wǎng)格拆分情況

采用一個(gè)線性界面拆分凸四邊形網(wǎng)格，存在如圖2所示的以下4種情況：

圖2 四邊形網(wǎng)格一次拆分可能性

① 2個(gè)交點(diǎn)位于相鄰的2條邊上，線性界面將四邊形網(wǎng)格拆分成1個(gè)三角形和1個(gè)五邊形；

② 2個(gè)交點(diǎn)位于一組對邊上，線性界面將四邊形拆分為2個(gè)四邊形；

③ 1個(gè)交點(diǎn)與四邊形的1個(gè)節(jié)點(diǎn)重合，另一個(gè)交點(diǎn)在這個(gè)節(jié)點(diǎn)的2條對邊之一上，線性界面將四邊形拆分成一個(gè)三角形和一個(gè)四邊形；

④ 2個(gè)交點(diǎn)分別是對角的2個(gè)四邊形節(jié)點(diǎn)，線性界面將四邊形網(wǎng)格拆分成2個(gè)三角形.

其中情況③和情況④可以視作情況①和情況②的一種特殊情況. 因此，在進(jìn)行解析解計(jì)算時(shí)，只需要考慮三角形子網(wǎng)格和四邊形子網(wǎng)格的情況，對五邊形子網(wǎng)格的情況可以通過對側(cè)的三角形的形心x(1)與體積V(1)和整個(gè)網(wǎng)格的形心xcell與體積Vcell來推導(dǎo)出五邊形子網(wǎng)格的形心x(2)，

(9)

3 任意凸四邊形的解析解求法

在本節(jié)中，首先考慮的是交點(diǎn)位于AD邊和AB邊時(shí)拆分出三角形的情況，其他拆分出三角形的情況可以通過軸變換直接轉(zhuǎn)變到這種情況；其次考慮交點(diǎn)位于AD邊和BC邊時(shí)的拆分出2個(gè)四邊形的情況，同樣的，其他拆分出四邊形的情況也可以通過軸變換的方法來轉(zhuǎn)變成這種情況.

3.1 拆分出三角形的情況

圖3 四邊形網(wǎng)格拆分出三角形網(wǎng)格的示意圖

因此，對情況①中的三角形，其形心和體積可以給定為

(10)

可以推出x*和y*之間的關(guān)系為

(11)

根據(jù)二次曲線的性質(zhì)[12]，三角形形心(x*,y*)軌跡曲線為一條雙曲線. 對于給定形心(Px,Py)，在這條雙曲線上距離給定形心最近的一點(diǎn)(x*,y*)與給定形心(Px,Py)的連線一定垂直于雙曲線在點(diǎn)(x*,y*)處的切線，求解得：

(12)

圖4 四邊形網(wǎng)格拆分出三角形網(wǎng)格的可能范圍

3.2 拆分出四邊形的情況

對于四邊形網(wǎng)格拆分出2個(gè)四邊形的情況，類似于3.1中的情況，首先進(jìn)行旋轉(zhuǎn)變換，使2個(gè)交點(diǎn)R、S分別位于AD邊和BC邊上；然后給定各點(diǎn)坐標(biāo)數(shù)據(jù)，如圖5所示.

圖5 四邊形網(wǎng)格拆分出四邊形網(wǎng)格的示意圖

由于沒有一個(gè)確定的公式直接給出四邊形形心，使用將四邊形ABSR拆分成2個(gè)三角形ASR和ABS的方法，因此，子網(wǎng)格ABSR體積和形心分別為

(13)

求解方程組得到：

(14)

式中：

Δ=4V*2A，A=

符號±表示此時(shí)yS是正的，yS>yR，反之亦然.

將求出yS,yR代入式(13)可以得x*和y*之間的關(guān)系為

F4(x*,y*)=x*2-(k1+k2)x*y*+

(15)

根據(jù)二次曲線標(biāo)準(zhǔn)型不變量的性質(zhì)[12]，當(dāng)k1=k2時(shí)子網(wǎng)格形心(x*,y*)軌跡曲線為拋物線；當(dāng)k1≠k2時(shí)子網(wǎng)格形心(x*,y*)軌跡曲線為雙曲線.

當(dāng)k1≠k2時(shí)，實(shí)際形心(Px,Py)的最接近位置(x*,y*)應(yīng)滿足兩點(diǎn)連線垂直于曲線在(x*,y*)點(diǎn)處的切線條件. 即：

(x*-Px,y*-Py)×

(16)

求解得到x*和y*關(guān)系后代回到式(13)，化簡得到

(HF2-36C2)y*4+(2FGH+F2J-72CD)y*3+

(HG2+KF2+2FGJ-72CE-36D2)y*2+

(JG2+2FGK-72DE)y*+(KG2-36E2)=0

(17)

式中：

H=9(k2-k1)2；

K=-[8V*(k2-k1)+3].

當(dāng)k2≠k1時(shí),HF2-36C2=36(k2-k1)2(k22+1)(k12+1)≠0，式(17)恒為4次方程，可以采用文獻(xiàn)[13]的方法得到解析解，至多有4個(gè)實(shí)根. 又有式x*=f4(y*)，故對每一個(gè)求出的y*，都可以求解出唯一一個(gè)x*.

k1=k2時(shí)，式(17)可以修正三次方程，可以采用文獻(xiàn)[13]的方法得到解析解，至多有3個(gè)實(shí)根. 又有式x*=f4(y*)，故對每一個(gè)求出的y*，都可以求解出唯一一個(gè)x*.

3.3 一般情況拓展

在2.3節(jié)中，推導(dǎo)了式(9)來計(jì)算五邊形子網(wǎng)格的形心. 當(dāng)凸四邊形網(wǎng)格的A點(diǎn)位于五邊形子網(wǎng)格中時(shí)，可以使用式(9)對圖形進(jìn)行適當(dāng)?shù)男D(zhuǎn)以簡化計(jì)算.

從幾何角度考慮，當(dāng)子網(wǎng)格體積V*>SΔABC時(shí)，可以選擇對側(cè)三角形DAC作為拆分三角形計(jì)算，也就是將點(diǎn)D平移至原點(diǎn)，將點(diǎn)C旋轉(zhuǎn)并縮放至(1,0)，之后進(jìn)行計(jì)算. 總的來說，即選擇三角形的第一個(gè)點(diǎn)為原點(diǎn)，第二個(gè)點(diǎn)為(1,0)點(diǎn). 進(jìn)行該情況的計(jì)算. 對于所有可能的情況，選擇計(jì)算三角形如下：

V*>SΔABC時(shí)，選擇ΔDAC求解；V*>SΔABD時(shí)，選擇ΔBCD求解;

V*>SΔBCA時(shí)，選擇ΔCDA求解；V*>SΔCDB時(shí)，選擇ΔDAB求解.

4 界面重構(gòu)算例

4.1 隨機(jī)四邊形網(wǎng)格中進(jìn)行界面重構(gòu)

在尺寸為1×1，網(wǎng)格數(shù)量為160×160的隨機(jī)四邊形網(wǎng)格中構(gòu)建一個(gè)圓心為(0.5,0.5)，半徑為0.35的圓形，結(jié)果如圖6所示，圖6(a)為整體重構(gòu)界面圖，圖6(b)為局部界面圖.

圖6 圓心為(0.5,0.5),半徑為0.35的圓在隨機(jī)四邊形網(wǎng)格中的重構(gòu)

4.2 界面重構(gòu)效率比較

本節(jié)采用重構(gòu)同心圓的方法對傳統(tǒng)Brent迭代方法和解析解求法的計(jì)算效率進(jìn)行比較. 在[0,1]×[0,1]的網(wǎng)格上重構(gòu)10個(gè)同心圓，其圓心均為(0.5,0.5)，半徑為0.05至0.5(間隔0.05). 在網(wǎng)格密度為30×30的情況下重構(gòu)結(jié)果如圖7所示，傳統(tǒng)方法耗時(shí)47.95 s，解析解方法4.68×10-2s，計(jì)算效率提高了約1 000倍.

圖7 采用傳統(tǒng)及解析方法重構(gòu)10個(gè)同心圓

將網(wǎng)格密度增加到500×500后，傳統(tǒng)方法用時(shí)403.17 s，解析解方法用時(shí)0.36 s.

4.3 二維強(qiáng)激波沖擊水中氣泡

本節(jié)采用二維強(qiáng)激波沖擊水中氣泡算例[14]，驗(yàn)證算法在應(yīng)用中的穩(wěn)定性及計(jì)算效率.

空氣和水兩種介質(zhì)采用統(tǒng)一的狀態(tài)方程[15]，p=(γ-1)ρe-γB. 空氣、靜止的水及水中強(qiáng)激波的初始密度、速度、壓力及狀態(tài)方程參數(shù)取值為

式中：p為壓力；ρ為密度(g/cm3)；u、v分別為x、y方向的速度(cm/ms)；γ和B為狀態(tài)方程參數(shù).

計(jì)算域選擇為12.0 cm×15.0 cm，網(wǎng)格步長取0.1 cm；x<3 cm設(shè)置為強(qiáng)激波，參數(shù)取Ωs，設(shè)置圓心為(3 cm,6 cm)，半徑3 cm的圓形氣泡，參數(shù)取Ωa，其余部分為靜止水，參數(shù)取Ωw；計(jì)算域左側(cè)設(shè)為入流邊界條件，其余邊界設(shè)為固壁有滑移條件.

典型時(shí)刻的密度梯度圖像如圖8(左為解析方法計(jì)算結(jié)果，右為傳統(tǒng)方法計(jì)算結(jié)果)：

圖8 典型時(shí)刻的密度梯度圖像

本文還對比了解析方法與傳統(tǒng)方法的計(jì)算效率. 運(yùn)行1 000步，傳統(tǒng)方法MoF重構(gòu)部分用時(shí)34 289.36 s；解析方法MoF重構(gòu)部分用時(shí)26.83 s，計(jì)算效率提升顯著.

5 結(jié)束語

本文基于MoF界面重構(gòu)方法提出了一種適用于任意凸四邊形網(wǎng)格的解析界面重構(gòu)算法，通過分別討論不同子網(wǎng)格拆分情況下的子網(wǎng)格形心曲線，解析地求解出界面重構(gòu)的最小化問題，實(shí)現(xiàn)MoF方法界面重構(gòu)，大幅提高了界面重構(gòu)效率. 采用隨機(jī)四邊形網(wǎng)格內(nèi)界面重構(gòu)驗(yàn)證，證明了界面重構(gòu)的有效性和準(zhǔn)確性；采用激波沖擊氣泡的算例驗(yàn)證了算法在實(shí)際應(yīng)用中的有效性和對計(jì)算效率的提高.