一種中心型間斷有限元MMALE方法

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

(1. 北京應用物理與計算數學研究所, 計算物理實驗室, 北京 100088;2. 中國工程物理研究院研究生院, 北京 100088)

0 引 言

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

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

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

1 MMALE方法的計算步驟

本文的MMALE方法的計算步驟如圖1所示。首先是初始化步：定義初始網格上所有物理量的分布，使用文獻[41]中2.3.1節的方法定義初始網格中物質體積分數和中心坐標。在拉氏步，采用單元中心型間斷有限元方法求解流體動力學方程，并將物理量和網格從時間tn更新到時間tn+1=tn+Δt上。混合單元的封閉模型采用Tipton壓力松弛模型[42-43]。使用文獻[44]中等參坐標方法更新物質的中心坐標。界面重構采用卿芳等提出的一種健壯的MOF算法[45]。在網格重分步，利用Knupp算法[46]對拉氏網格進行光滑處理；或者采用初始網格作為重分網格。最后，在重映步，使用二階積分守恒重映方法，將時間tn+1的拉氏網格中的所有物理量映射到重分后的新網格上。

圖1 MMALE 方法流程圖

2 間斷有限元拉氏方法

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

(1)

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

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

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

(2)

對上式第一項用Reynolds輸運定理以及對第二項用散度定理，最終得到：

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

?Wk,k=1,…,N}

(5)

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

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

標準單元Ωst的四個基函數為σ1=1,σ2=ξ,σ3=η,σ4=ξη。由于Fk是光滑的可逆映射，因此：

i=1,…,4

(8)

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

圖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]。時間步長的估計遵循標準的CFL準則和每個時間步上限制單元面積的變化，詳情請參考文獻[16-17]。空間二階重構采用最小二乘算法[17]，并采用Barth-Jesperson限制器[49]抑制非物理的數值振蕩。

3 界面重構

圖3 MOF方法示意圖

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

注意到f(θ)的一階導數[50-51]，

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

4 重映

4.1 多項式重構

4.2 新舊網格相交

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

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

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

1) 仿射變換

2) “剪裁投影”算法

圖6 單位三角形和多邊形的交點

3) 相交多邊形的面積和中心坐標

最后，通過坐標變換，將參考坐標系中相交多邊形的面積和中心坐標轉換為物理坐標系中面積和中心坐標：

4.3 積分

4.4 后驗校正

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

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

1) 物理相容性檢驗準則(PAD)

該準則的依據是密度、內能和壓強的必須為正：

2) 數值相容性檢驗準則(NAD)

該準則是為了檢驗數值解是否局部光滑：

DNAD=

圖7 基于后驗校正的重映算法示意圖

5 數值算例

5.1 二維周期渦問題[36]

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

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

5.2 圓柱內爆中的Rayleigh-Taylor不穩定性問題[53]

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

圖8 初始區域的幾何數據

圖9 初始網格

對于MMALE方法，在開始時初始網格沒有混合單元，因此僅需執行拉氏步。當網格由于界面不穩定而發生嚴重變形時，需要進行重分和重映。如果網格中存在混合單元，采用Tipton壓力松弛模型進行計算，并在重映步之前對這些混合單元執行MOF界面重構算法。

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

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

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

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

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

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

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

5.3 Sedov 問題[54]

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

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

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

5.4 二維Rayleigh-Taylor 不穩定性問題[6，36]

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

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

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

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

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

計算區域為[0,0.65]×[-0.089,0.089]，初始時刻將它剖分成260×72個均勻單元。如圖19所示，氣泡是由圓心(xc,yc)=(0.32,0)和半徑r=0.025定義的圓盤。在右邊界(初始時刻位于x=0.65)上給出了一個類似活塞的邊界條件，它以速度(u*,0)向左移動。在其他邊界上給固壁邊界條件。入射激波的馬赫數為Ms=1.22。氦氣泡和空氣在初始時刻處于靜止狀態。氦氣泡和空氣的初始數據分別為(ρ1,p1)=(0.182,105),(ρ2,p2)=(1,105)；摩爾質量和絕熱指數分別為(M1,γ1)=(5.269×10-3,1.648)，(M2,γ2)=(28.963×10-3,1.4)。由Rankine-Hugoniot關系可以得到活塞在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三個時刻的網格和界面以及這三個時刻的試驗紋影圖[55]。圖21 給出了tfinal時刻的密度等值線圖以及局部放大圖。從圖20和圖21可以看出，這里的計算結果與文獻[55]中的試驗紋影圖符合的很好，與文獻[6，36]的數值結果非常接近。

圖19 初始區域的幾何數據

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

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

5.6 三重點問題[6]

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

圖22 初始區域的幾何數據

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

6 結 論

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