卿 芳, 蔚喜軍, 趙曉龍, 鄒世俊, 賈祖朋,*

(1. 北京應(yīng)用物理與計算數(shù)學(xué)研究所, 計算物理實(shí)驗(yàn)室, 北京 100088;2. 中國工程物理研究院研究生院, 北京 100088)

0 引 言

多介質(zhì)可壓縮流體動力學(xué)問題的數(shù)值計算方法主要分為兩類：歐拉方法(Eulerian method)和拉氏方法(Lagrangian method)。歐拉方法采用固定的空間網(wǎng)格進(jìn)行計算。拉氏方法的網(wǎng)格跟隨流體一起運(yùn)動。這兩種方法各有利弊。歐拉方法易于處理流場大變形，但是難以精確地捕捉物質(zhì)界面。拉氏方法能清晰地刻畫、描述物質(zhì)界面，但是當(dāng)流場中發(fā)生大變形時，網(wǎng)格會嚴(yán)重扭曲，導(dǎo)致計算精度下降甚至計算終止。為了結(jié)合這兩種方法的優(yōu)點(diǎn)，Hirt等[1]提出了ALE(Arbitrary Lagrangian-Eulerian)方法，其中網(wǎng)格可以以任意指定的方式移動。大多數(shù)ALE算法[2-7]包括三步：拉氏步、重分步和重映步。傳統(tǒng)的ALE方法用單元邊描述物質(zhì)界面，并且要求每個單元只有一種介質(zhì)。當(dāng)物質(zhì)界面變形嚴(yán)重時，很難生成質(zhì)量良好的新網(wǎng)格。為了解決這個問題，Peer等[8]將傳統(tǒng)的ALE方法與界面捕捉方法相結(jié)合，提出了MMALE(Multi-material ALE)方法。這種方法引入了混合單元(一個單元中存在多種介質(zhì))，在混合單元中重構(gòu)物質(zhì)界面，重分網(wǎng)格時可以跨過界面生成新的網(wǎng)格，因此比傳統(tǒng)的ALE方法提供了額外的靈活性。MMALE方法的計算流程與傳統(tǒng)的ALE方法類似。首先，在拉氏步更新每種介質(zhì)的物理量以及網(wǎng)格坐標(biāo)，其中對混合單元需要用熱力學(xué)封閉模型來更新各種介質(zhì)的物理量。當(dāng)網(wǎng)格變形嚴(yán)重時，需要重構(gòu)混合單元的物質(zhì)界面，為后面重映步提供界面信息。然后，在重分步生成一個幾何品質(zhì)較好的新網(wǎng)格。最后,在重映步將舊網(wǎng)格的物理量映射到新網(wǎng)格上。

在拉氏步，常用的離散方法有兩種。一種是交錯型格式[9-12]，指的是速度定義在節(jié)點(diǎn)上，而其他變量(密度、壓力和內(nèi)能等)定義在單元中心。另一種是中心型格式[13-17]，所有變量都定義在單元中心。中心型格式比交錯型格式有一些優(yōu)點(diǎn)。例如，交錯型格式對變量使用不同的控制體,很難對所有變量構(gòu)建一致的高階精度。本文采用中心型格式。中心型格式可以分為兩類: 中心型有限體積方法[16-17]和中心型間斷有限元方法[18-20]。有限體積方法需要較大的單元模板來構(gòu)造高階格式；而間斷有限元是利用單元內(nèi)的高階插值多項(xiàng)式來構(gòu)造高階格式。因此，間斷有限元方法比有限體積方法更具有發(fā)展高階格式的前景。重映步通常采用兩種方法。一種是對流重映[21-25]，它要求新舊網(wǎng)格離得比較近，計算效率比較高.另一種是積分重映[26-31],它對新舊網(wǎng)格的關(guān)聯(lián)要求不高，但是需要計算新舊網(wǎng)格的交點(diǎn)，計算結(jié)果比較精確。近年來，在MMALE方法的研究上取得了一些進(jìn)展[32-35]。Galera和Maire等[33-34]提出了一種單元中心型二維MMAIE方法，界面重構(gòu)采用VOF(Volume of Fluid)或MOF方法。賈祖朋等[29]提出了一種基于MOF界面重構(gòu)的交錯型非結(jié)構(gòu)MMALE方法，模擬涉及強(qiáng)剪切變形的多介質(zhì)可壓縮流動。后來，賈祖朋等[36]發(fā)展了一種基于改進(jìn)的MOF方法的中心型MMALE方法。上面提到的MMALE方法在拉氏步都是用有限體積方法。

本文提出了一種二維中心型間斷有限元MMALE方法。在拉氏步，使用間斷有限元方法對流體力學(xué)方程組進(jìn)行離散，選取文獻(xiàn)[19]中的基函數(shù)。為了構(gòu)造二階積分重映算法，對拉氏步獲得的物理量單元均值做分片線性重構(gòu)。經(jīng)典的限制方法是對重構(gòu)多項(xiàng)式的梯度進(jìn)行限制[26,37]。近些年,Blanchard和Loubere[38]提出了基于MOOD (Multi-dimensional Optimal Order Detection)限制策略[39-40]的后驗(yàn)校正。本文采用后驗(yàn)校正，并做了一些改動以適應(yīng)多介質(zhì)計算。對新舊網(wǎng)格相交的計算，借鑒文獻(xiàn)[31]的三維“剪裁投影”算法的基本思想，設(shè)計了二維“剪裁投影”算法，顯著降低了多邊形相交算法的復(fù)雜度。

1 MMALE方法的計算步驟

本文的MMALE方法的計算步驟如圖1所示。首先是初始化步：定義初始網(wǎng)格上所有物理量的分布，使用文獻(xiàn)[41]中2.3.1節(jié)的方法定義初始網(wǎng)格中物質(zhì)體積分?jǐn)?shù)和中心坐標(biāo)。在拉氏步，采用單元中心型間斷有限元方法求解流體動力學(xué)方程，并將物理量和網(wǎng)格從時間tn更新到時間tn+1=tn+Δt上?；旌蠁卧姆忾]模型采用Tipton壓力松弛模型[42-43]。使用文獻(xiàn)[44]中等參坐標(biāo)方法更新物質(zhì)的中心坐標(biāo)。界面重構(gòu)采用卿芳等提出的一種健壯的MOF算法[45]。在網(wǎng)格重分步，利用Knupp算法[46]對拉氏網(wǎng)格進(jìn)行光滑處理；或者采用初始網(wǎng)格作為重分網(wǎng)格。最后，在重映步，使用二階積分守恒重映方法，將時間tn+1的拉氏網(wǎng)格中的所有物理量映射到重分后的新網(wǎng)格上。

圖1 MMALE 方法流程圖

2 間斷有限元拉氏方法

歐拉框架下的二維無粘可壓縮Euler方程組的向量形式為：

(1)

其中U=(ρ,ρu,ρv,ρE)T，F(xiàn)=(0,p,0,pu)T，G=(0,0,p,pv)T，?·是散度算子，即?·=?/?x+?/?y。ρ是流體密度,u=(u,v)是流體速度，p是壓力，E是總能量，相應(yīng)的比內(nèi)能為e=E-u2/2，狀態(tài)方程為p=p(ρ,e)。

下面給出式(1)的弱形式。式(1)兩邊同時乘以函數(shù)φ，并經(jīng)過微分運(yùn)算，得：

+?·(φ(F,G))-(F,G)·?φ=0

(2)

對上式第一項(xiàng)用Reynolds輸運(yùn)定理以及對第二項(xiàng)用散度定理，最終得到：

其中n=(nx,ny) 表示單元邊界S的單位外法向量。

?Wk,k=1,…,N}

(5)

其中Pn(Wk)為單元Wk上不超過n次多項(xiàng)式的全體集合。那么，對任意的單元Wk，下面的方程成立：

定義標(biāo)準(zhǔn)單元Ωst=[-1,1]2，它的坐標(biāo)記為(ξ,η)，ξ,η∈[-1,1]。對于任意的時間t，存在一個從Ωst到Wk的雙線性映射Fk:(ξ,η)∈Ωst→(x,y)∈Wk:

標(biāo)準(zhǔn)單元Ωst的四個基函數(shù)為σ1=1,σ2=ξ,σ3=η,σ4=ξη。由于Fk是光滑的可逆映射，因此：

i=1,…,4

(8)

數(shù)值解在單元的交界面上可以是間斷的。因此在單元邊界上定義一個數(shù)值通量函數(shù)F*和G*來代替原來的通量函數(shù)F和G[20]。邊[Mr,Mr+1]上的數(shù)值通量定義為：

圖2 單元的記號

將(13)-(16)代入(12)：

m=1,2,3,4

(18)

m=1,2,3,4

(19)

最后，如果令：

uk1=[ρk1,ρk2,ρk3,ρk4]

(22)

uk2=[ρuk1,ρuk2,ρuk3,ρuk4]

(23)

uk3=[ρvk1,ρvk2,ρvk3,ρvk4]

(24)

uk4=[ρEk1,ρEk2,ρEk3,ρEk4]

(25)

其中Lkj(j=1,2,3,4)是(17)～(20)式空間離散算子得到的向量。

半離散格式(26)的時間離散采用二階顯式Runge-Kutta方法[48]。時間步長的估計遵循標(biāo)準(zhǔn)的CFL準(zhǔn)則和每個時間步上限制單元面積的變化，詳情請參考文獻(xiàn)[16-17]?？臻g二階重構(gòu)采用最小二乘算法[17]，并采用Barth-Jesperson限制器[49]抑制非物理的數(shù)值振蕩。

3 界面重構(gòu)

圖3 MOF方法示意圖

其中xh(θ)表示直線截取的多邊形的中心點(diǎn)。文獻(xiàn)[50]用最優(yōu)化方法求解上述極小化問題，有可能收斂到局部極小值。因此，近年有學(xué)者對經(jīng)典的MOF方法做了一些改進(jìn)。

注意到f(θ)的一階導(dǎo)數(shù)[50-51]，

文獻(xiàn)[51]中提到f′(θ)是關(guān)于θ的連續(xù)函數(shù)。因此，f(θ)的最小值點(diǎn)一定是f′(θ)的零點(diǎn)或者區(qū)間端點(diǎn)。文獻(xiàn)[45]提出了一種健壯的算法計算f′(θ)的零點(diǎn)。本文采用這種健壯的MOF方法[45]。

4 重映

4.1 多項(xiàng)式重構(gòu)

4.2 新舊網(wǎng)格相交

圖4 新單元和拉氏網(wǎng)格{W}的相交部分。拉氏網(wǎng)格和新網(wǎng)格分別用藍(lán)線和黑線表示。所要求的相交部分分別用不同顏色(a-f)標(biāo)記

圖5 單元內(nèi)物質(zhì)用物質(zhì)界面(紅線)分隔開。單介質(zhì)子多邊形用Pgreen，Pcyan表示

下面介紹求相交子多邊形的一個算法。這個算法遵循陳享等[31]的三維“剪裁投影”算法基本思想，并改為退化的二維“剪裁投影”算法。由于新單元是四邊形，可以先把它分解成兩個三角形。然后用這兩個三角形分別和舊網(wǎng)格求相交子多邊形，因而，新單元的面積等于這兩個三角形分別得到的相交子多邊形的面積之和。這個新舊網(wǎng)格相交的算法歸結(jié)為求三角形與多邊形的交點(diǎn)。算法分為以下三步：

1) 仿射變換

2) “剪裁投影”算法

圖6 單位三角形和多邊形的交點(diǎn)

3) 相交多邊形的面積和中心坐標(biāo)

最后，通過坐標(biāo)變換，將參考坐標(biāo)系中相交多邊形的面積和中心坐標(biāo)轉(zhuǎn)換為物理坐標(biāo)系中面積和中心坐標(biāo)：

4.3 積分

4.4 后驗(yàn)校正

一般而言，重構(gòu)多項(xiàng)式(28)需要一個限制器來抑制間斷附近的振蕩。本文采用文獻(xiàn)[38]中提出的基于MOOD限制策略的后驗(yàn)校正。后驗(yàn)校正的原理是檢測出“壞”單元，然后降低“壞”單元上重構(gòu)多項(xiàng)式的次數(shù)，重復(fù)之前的步驟，直到這個單元是“好”單元為止。值得注意的是，文獻(xiàn)[38]中提出的方法只適用于單介質(zhì)流。為了使其適應(yīng)于多介質(zhì)流，我們對它做了一些小改動。對于純(單介質(zhì))新單元，對新單元進(jìn)行檢測；對于混合(多介質(zhì))新單元，對所有單介質(zhì)子多邊形進(jìn)行檢測，而不是整個新單元。

檢測的目的是過濾“好”單元，這些單元已經(jīng)使用無限制的重構(gòu)多項(xiàng)式進(jìn)行精確積分守恒重映；而把存在一些問題(振蕩，非物理等)的“壞”單元標(biāo)識出來。檢測準(zhǔn)則有兩個：

1) 物理相容性檢驗(yàn)準(zhǔn)則(PAD)

該準(zhǔn)則的依據(jù)是密度、內(nèi)能和壓強(qiáng)的必須為正：

2) 數(shù)值相容性檢驗(yàn)準(zhǔn)則(NAD)

該準(zhǔn)則是為了檢驗(yàn)數(shù)值解是否局部光滑：

DNAD=

圖7 基于后驗(yàn)校正的重映算法示意圖

5 數(shù)值算例

5.1 二維周期渦問題[36]

表1 誤差及收斂階(拉氏方法)

表2 誤差及收斂階(MMALE方法)

5.2 圓柱內(nèi)爆中的Rayleigh-Taylor不穩(wěn)定性問題[53]

這是一個類似ICF的測試問題。初始狀態(tài)如圖8所示，一團(tuán)輕氣體(R∈[0,1])被一團(tuán)重氣體包圍(R∈[0,1.2]))。兩種氣體均為理想氣體(絕熱指數(shù)γl=γh=1.5)。輕氣體的初始狀態(tài)為(ρl,pl,Ul)=(0.05,0.1,0),重氣體為(ρh,ph,Uh)=(1,0.1,0)。在重氣體外表面給定隨時間變化的壓力邊界條件，

圖8 初始區(qū)域的幾何數(shù)據(jù)

圖9 初始網(wǎng)格

對于MMALE方法，在開始時初始網(wǎng)格沒有混合單元，因此僅需執(zhí)行拉氏步。當(dāng)網(wǎng)格由于界面不穩(wěn)定而發(fā)生嚴(yán)重變形時，需要進(jìn)行重分和重映。如果網(wǎng)格中存在混合單元，采用Tipton壓力松弛模型進(jìn)行計算，并在重映步之前對這些混合單元執(zhí)行MOF界面重構(gòu)算法。

當(dāng)取a=10-3時，圖10給出t=0.3時用拉氏方法計算的網(wǎng)格和密度圖。注意到最終的網(wǎng)格質(zhì)量保持很好，沒有翻轉(zhuǎn)單元。圖11給出的是t=0.3時用MMALE方法計算網(wǎng)格與界面和密度等值線圖。從圖10和圖11可見，拉氏方法和MMALE方法的計算結(jié)果是非常相似的。

圖10 a=2×10-4，t=0.3時用間斷有限元拉氏方法計算的網(wǎng)格(左)和密度圖(右)

圖11 a=2×10-4，t=0.3時用MMALE方法計算的網(wǎng)格(左)和密度(右)

當(dāng)取a=10-3時，圖12給出了t=0.25時用拉氏方法計算的網(wǎng)格和局部放大圖，此時網(wǎng)格是沒有翻轉(zhuǎn)單元的。當(dāng)t=0.3時，用拉氏方法計算的結(jié)果是不準(zhǔn)確的，此時網(wǎng)格扭曲嚴(yán)重而且存在翻轉(zhuǎn)單元，如圖13所示。然而，用MMALE方法運(yùn)行到最終時刻沒有任何問題。圖14給出的是t=0.3時用MMALE方法計算的網(wǎng)格與界面和密度等值線圖。該算例的計算結(jié)果表明了MMALE方法具有較好的靈活性和魯棒性。

圖12 a=10-3，t=0.25時用間斷有限元拉氏方法計算的網(wǎng)格(左)和局部放大圖(右)

圖13 a=10-3，t=0.3時用間斷有限元拉氏方法計算的網(wǎng)格(左)和局部放大圖(右)

圖14 a=10-3，t=0.3時用MMALE方法計算的網(wǎng)格(左)和密度(右)

5.3 Sedov 問題[54]

這是一個點(diǎn)爆炸問題。在原點(diǎn)處給定一個能量。隨著時間的推進(jìn)，一個圓柱型激波會向外傳播。初始密度為ρ=1，速度為(u,v)=(0,0)。在原點(diǎn)處，存在單位內(nèi)能，并且在其它點(diǎn)處初始壓力為p=0。這個問題的精確解由Sedov在文獻(xiàn)[54]中給出：在時間t=1時，激波運(yùn)動到離原點(diǎn)距離為1處，密度峰值達(dá)到6。使用本文的中心型MMALE方法進(jìn)行計算。計算區(qū)域?yàn)棣?[0,1.125]×[0,1.125]，初始剖分為45×45的一致網(wǎng)格。在靠近原點(diǎn)的第一個單元上給定內(nèi)能e0=400。在初始時刻，界面位于以原點(diǎn)為圓心，半徑為0.43的圓上。注意到在混合單元中，這兩種介質(zhì)都是理想氣體，具有相同的絕熱指數(shù)γ=1.4。我們將它們視為混合單元，以便將數(shù)值解與精確解進(jìn)行比較。圖15 給出了t=0.5,0,75,1時的網(wǎng)格和界面。圖16是t=1時刻的密度等值線圖和以半徑為變量的密度圖。從圖16可知。計算得到的激波位置是準(zhǔn)確的，數(shù)值解與精確解很接近。

圖15 t=0.5,0,75,1時的網(wǎng)格和界面(從左到右)

圖16 t=1時的密度等值線圖(左)和以半徑為變量的密度圖(右)

5.4 二維Rayleigh-Taylor 不穩(wěn)定性問題[6，36]

計算區(qū)域?yàn)榫匦蝃0,1/3]×[0,1].初始時刻將它剖分為35×100的均勻網(wǎng)格。初始時刻，密度為ρh=2的重氣體和密度為ρl=1的輕氣體被擾動界面分隔開，界面方程為yi(x)=0.5+0.01cos(6πx)。重氣體在輕氣體之上。這兩種氣體都有相同的絕熱指數(shù)γ=1.4。重力場垂直向下，重力加速度為g=(gx,gy)T=(0,-0.1)T。初始時刻兩種氣體都處于靜止?fàn)顟B(tài)。初始壓力分布通過靜壓平衡條件推導(dǎo)出來:

眾所周知，這種狀態(tài)是不穩(wěn)定的。由于正弦界面的存在，在界面附近會產(chǎn)生旋渦。隨著時間的推移，較重的氣體會落入較輕的氣體中形成一個尖釘；較輕的氣體會上升到較重的氣體中形成一個氣泡。圖17給出了t=5,7,9,11,13時的網(wǎng)格和界面。圖18給出了t=5,7,9,11,13時的密度等值線圖。從圖17和圖18可以看出，我們的結(jié)果與文獻(xiàn)[6,29]中的結(jié)果非常接近，這里的計算結(jié)果對稱性好，界面清晰。

圖17 t=5,7,9,11,13時的網(wǎng)格和界面(從左到右)

圖18 t=5,7,9,11,13時的密度等值線圖(從左到右)

5.5 激波與氦氣泡相互作用問題[6，36]

計算區(qū)域?yàn)閇0,0.65]×[-0.089,0.089]，初始時刻將它剖分成260×72個均勻單元。如圖19所示，氣泡是由圓心(xc,yc)=(0.32,0)和半徑r=0.025定義的圓盤。在右邊界(初始時刻位于x=0.65)上給出了一個類似活塞的邊界條件，它以速度(u*,0)向左移動。在其他邊界上給固壁邊界條件。入射激波的馬赫數(shù)為Ms=1.22。氦氣泡和空氣在初始時刻處于靜止?fàn)顟B(tài)。氦氣泡和空氣的初始數(shù)據(jù)分別為(ρ1,p1)=(0.182,105),(ρ2,p2)=(1,105)；摩爾質(zhì)量和絕熱指數(shù)分別為(M1,γ1)=(5.269×10-3,1.648)，(M2,γ2)=(28.963×10-3,1.4)。由Rankine-Hugoniot關(guān)系可以得到活塞在x方向的速度u*=-124.824。入射激波的x方向速度是Dc=-456.482。入射激波在ti=668.153×10-6時刻撞擊氦氣泡。計算終止時間tfinal=ti+674×10-6=1342.153×10-6。圖20給出了t1=ti+245×10-6=913.153×10-6，t2=ti+4247×10-6= 1095.153×10-6和tfinal三個時刻的網(wǎng)格和界面以及這三個時刻的試驗(yàn)紋影圖[55]。圖21 給出了tfinal時刻的密度等值線圖以及局部放大圖。從圖20和圖21可以看出，這里的計算結(jié)果與文獻(xiàn)[55]中的試驗(yàn)紋影圖符合的很好，與文獻(xiàn)[6，36]的數(shù)值結(jié)果非常接近。

圖19 初始區(qū)域的幾何數(shù)據(jù)

圖20 t1，t2, tfinal時刻(從左到右)用MMALE方法計算的密度(上)和試驗(yàn)紋影圖(下)

圖21 tfinal時刻的密度等值線圖(左) 和局部放大圖(右)

5.6 三重點(diǎn)問題[6]

考慮如圖22所示矩形域中的三種狀態(tài)的二維黎曼問題。初始區(qū)域Ω=[0,7]×[0,3]分為三個子區(qū)域Ω1=[0,1]×[0,3]；Ω2=[1,7]×[0,1.5]和Ω3=[1,7]×[1.5,3]。Ω1中是初始狀態(tài)為(ρ1,p1,U1)=(1,1,0)的高壓高密度氣體。Ω2中是初始狀態(tài)為(ρ2,p2,U2)=(1,0.1,0)的低壓高密度氣體。Ω3中是初始狀態(tài)為(ρ3,p3,U3)=(0.125,0.1,0)的低壓低密度氣體。Ω1和Ω3中物質(zhì)的絕熱指數(shù)均為γ1=γ3=1.5；Ω2中物質(zhì)的絕熱指數(shù)為γ2=1.4。邊界條件均為固壁邊界條件,時間計算到tfinal=10。初始網(wǎng)格數(shù)為210×90。由于聲阻抗的不同，Ω2和Ω3中的兩個激波以不同的速度傳播。這在Ω2和Ω3的界面處產(chǎn)生強(qiáng)剪切力。這種剪切作用產(chǎn)生Kelvin-Helmholtz不穩(wěn)定性，并形成渦旋。用拉氏方法計算該問題時，在渦旋發(fā)展前由于網(wǎng)格扭曲纏結(jié)而導(dǎo)致計算終止。當(dāng)使用經(jīng)典的ALE方法時，捕獲渦流是困難的。我們使用新的中心型MMALE方法計算此問題。圖23 給出了t=5,6時的網(wǎng)格和界面。

圖22 初始區(qū)域的幾何數(shù)據(jù)

圖23 t=5,6的網(wǎng)格和界面(從上到下)

6 結(jié) 論

本文提出了一種具有二階精度的中心型間斷有限元MMALE方法。在該方法中，拉氏步采用間斷有限元方法求解可壓縮歐拉方程。對于混合網(wǎng)格，采用Tipton壓力松弛模型更新物理量。采用具有魯棒性的MOF方法進(jìn)行界面重構(gòu)。針對重映步的多邊形相交問題，提出了一種“裁剪投影”算法；在此基礎(chǔ)上構(gòu)建了一個基于后驗(yàn)校正的二階積分守恒重映方法。數(shù)值試驗(yàn)的結(jié)果表明，該方法具有較好的魯棒性和靈活性。