強(qiáng)間斷多介質(zhì)流的高精度偽弧長(zhǎng)方法

王晨濤，馬天寶，李坤

（北京理工大學(xué) 爆炸科學(xué)與技術(shù)國(guó)家重點(diǎn)實(shí)驗(yàn)室，北京 100081）

基于連續(xù)介質(zhì)力學(xué)求解雙曲守恒律方程得到的解通常具有奇異性，即解存在強(qiáng)間斷，例如激波、接觸間斷，或者稀疏波之類的弱間斷. 針對(duì)求解間斷解發(fā)展了很多高分辨率和高精度的數(shù)值算法. 高分辨率的算法可以追溯到早期的TVD 總變差不增算法，但是算法精度不高于二階. 隨后，發(fā)展了一系列的高精度數(shù)值格式，例如ENO，RKDG，WENO[1-2]，SD，SV等. 但是高精度的數(shù)值格式數(shù)值色散較強(qiáng)，使得解的震蕩較強(qiáng). 通常需要額外的限制震蕩的方法來抑制震蕩. 因此，為了減少數(shù)值震蕩，提高間斷分辨率，從削弱方程奇異性的角度發(fā)展了弧長(zhǎng)算法. RISK[3]通過引入弧長(zhǎng)參數(shù)，增加了一個(gè)弧長(zhǎng)約束方程，使得求解過程的奇異性得到減弱或消除. CHAN[4]通過引入偽弧長(zhǎng)控制參數(shù)將其引申為偽弧長(zhǎng)算法，并應(yīng)用于求解非線性方程的的奇異性問題. 武際可[5]等人應(yīng)用偽弧長(zhǎng)算法解決了包含奇異性的Burgers 方程.TANG[6]提出了包含偽弧長(zhǎng)參數(shù)的移動(dòng)網(wǎng)格控制方程，并給出了在二維情況下網(wǎng)格最優(yōu)分布的自適應(yīng)函數(shù)表達(dá)式. 寧建國(guó)等發(fā)展了二階偽弧長(zhǎng)算法，建立了偽弧長(zhǎng)算法穩(wěn)定性[7]理論，得到了偽弧長(zhǎng)算法的收斂性結(jié)論，并成功將偽弧長(zhǎng)算法運(yùn)用到爆炸沖擊波強(qiáng)間斷問題的求解. 由于在變形的物理空間不易構(gòu)造高階數(shù)值格式，因此上述的弧長(zhǎng)算法限制在了二階精度.

在多介質(zhì)流體的求解過程中，圍繞物質(zhì)界面的追蹤發(fā)展了經(jīng)典的VOF[8-9]方法，Level Set[10]法，處理多介質(zhì)的五方程模型[11]方法，多介質(zhì)ALE[12]算法，以及基于界面形心重構(gòu)物質(zhì)界面的MOF[13]方法.Level Set 通過求解一個(gè)對(duì)流方程實(shí)現(xiàn)界面的隱式追蹤，無需對(duì)界面進(jìn)行繁瑣的重構(gòu)，因此具有計(jì)算簡(jiǎn)單，界面分辨率高等優(yōu)點(diǎn). 由于界面兩側(cè)流體狀態(tài)方程的差異，容易在界面附近產(chǎn)生需要虛假振蕩，因此需要合理定義界面的邊界條件. RGFM[14]通過求解一個(gè)黎曼問題，修改界面附近的物理量，有效降低了界面處的非物理反射引起的數(shù)值振蕩.

本文將偽弧長(zhǎng)算法和WENO 格式以及Level-Set和RGFM 相結(jié)合，發(fā)展了多介質(zhì)流的高精度偽弧長(zhǎng)算法，為了構(gòu)造高精度格式，采用坐標(biāo)變換將控制方程映射到正交均勻的弧長(zhǎng)坐標(biāo)系中，計(jì)算過程在均勻正交的弧長(zhǎng)坐標(biāo)下進(jìn)行計(jì)算. 為了消除變形網(wǎng)格幾何誤差，采用滿足幾何守恒律[15-17]的表達(dá)式計(jì)算. 針對(duì)網(wǎng)格移動(dòng)以后距離函數(shù)的插值提出了基于泰勒插值的三階非守恒格式. 計(jì)算結(jié)果表明，本文的算法在保持較高的精度的同時(shí)也削弱了間斷處的數(shù)值振蕩，大幅提升了間斷的分辨率.

1 控制方程

1.1 物理空間下的控制方程

可壓縮歐拉方程在物理空間的張量形式為

1.2 弧長(zhǎng)空間下的控制方程

可以看出，所有的表達(dá)式均是守恒格式，因此方便用中心差分格式計(jì)算，從而可以消除因網(wǎng)格變形造成的幾何誤差.

2 偽弧長(zhǎng)算法

2.1 偽弧長(zhǎng)算法幾何意義

為了削弱方程的奇異性，通過引入一個(gè)弧長(zhǎng)約束方程，即

其中s為弧長(zhǎng)變量. 可以發(fā)現(xiàn)，當(dāng)弧長(zhǎng)變量為常值時(shí)，dw越大，則 dx越小，它的物理意義就是物理空間下變量的梯度越大，網(wǎng)格的尺寸越小，也就是網(wǎng)格朝著間斷點(diǎn)移動(dòng). 相關(guān)的幾何意義參見圖1，從圖中可以看出，經(jīng)過引入弧長(zhǎng)約束方程，使得網(wǎng)格自適應(yīng)朝著間斷移動(dòng)，原始物理空間下陡峭的物理量在弧長(zhǎng)空間中得到了平滑的過渡，因此削弱了方程的奇異性.

圖1 1D 偽弧長(zhǎng)算法的幾何意義Fig. 1 Geometric significance of one-dimensional PALM

二維空間和一維空間類似，基于維度解耦的思想，引入兩個(gè)弧長(zhǎng)約束方程，把它寫成向量的形式，即

相關(guān)的幾何示意圖如圖2 所示. 它和一維情況下的物理意義類似，都是引入約束方程使得網(wǎng)格朝著間斷移動(dòng)，在物理空間的表現(xiàn)形式就是網(wǎng)格匯聚到了間斷附近，然后再弧長(zhǎng)空間下網(wǎng)格依舊是正交均勻的.

圖2 2D 偽弧長(zhǎng)算法的幾何意義Fig. 2 Geometric significance of two-dimensional PALM

2.2 偽弧長(zhǎng)算法網(wǎng)格移動(dòng)

偽弧長(zhǎng)算法的網(wǎng)格滿足一個(gè)能量極值函數(shù)：

然后可以給出滿足弧長(zhǎng)控制函數(shù)的網(wǎng)格迭代方程(8)，即

至此，網(wǎng)格便可以采用數(shù)值方法迭代計(jì)算，例如采用收斂較快的高斯-賽德爾迭代，格式如下.

3 數(shù)值格式

3.1 控制方程的離散

3.2 守恒變量的更新

網(wǎng)格移動(dòng)完以后，需要更新在新網(wǎng)格上的物理量，本文采用守恒插值[6]來計(jì)算：

4 Level-Set 技術(shù)

4.1 控制方程的離散

4.2 距離函數(shù)的非守恒插值

網(wǎng)格移動(dòng)完以后，需要更新新網(wǎng)格點(diǎn)上的Level-Set 值，本文基于泰勒插值，給出三階插值的表達(dá)式，由于新舊網(wǎng)格間距較小，因此將新網(wǎng)格點(diǎn)在舊網(wǎng)格點(diǎn)作泰勒展開，有：

由于計(jì)算過程是在弧長(zhǎng)坐標(biāo)系下，因此將上式映射至弧長(zhǎng)空間，可得：

4.3 多介質(zhì)流高精度偽弧長(zhǎng)算法流程

⑤根據(jù)式(10)和(18)進(jìn)行控制方程以及符號(hào)距離函數(shù)在時(shí)間步上的更新. 更新完以后根據(jù)距離函數(shù)的符號(hào)來重新分配兩種介質(zhì).

⑥判定計(jì)算時(shí)間是否到達(dá). 若到達(dá)，則終止計(jì)算，否則轉(zhuǎn)到步驟①.

5 數(shù)值算例

為便于計(jì)算，以下算例計(jì)算均采用無量綱計(jì)算.

5.1 精度測(cè)試

用一個(gè)包含精確解的一維的Euler 方程來驗(yàn)證高階偽弧長(zhǎng)算法PALM-5 的收斂精度. 初值條件和精確解為

從表1 的計(jì)算結(jié)果可以看出，弧長(zhǎng)坐標(biāo)系下的高階偽弧長(zhǎng)算法保證了較高的精度，沒有因?yàn)榫W(wǎng)格的變形而損失精度. 這說明本文的算法滿足了高精度的性質(zhì).

表1 固定網(wǎng)格下WENO 和PALM-5 的收斂階Tab. 1 The accuracy of WENO, PALM-5 schemes

5.2 一維爆炸波

除圖3(g)和3(h)算例以外其余均采用的網(wǎng)格數(shù)目是300，參照解是在WENO 格式下用20 000 個(gè)網(wǎng)格計(jì)算得到的. 從圖3(a)和(3b)中可以看出，采用相同的格式偽弧長(zhǎng)算法比固定網(wǎng)格的算法分辨率要高很多，對(duì)間斷的捕捉能力很強(qiáng). 而且5 階偽弧長(zhǎng)算法比2 階偽弧長(zhǎng)算法分辨率更高，對(duì)間斷的捕捉能力要高很多，耗散很低. 對(duì)比圖3(c)和3(d)可以看出采用不同的偽弧長(zhǎng)參數(shù)得到的結(jié)果差異較大，本算例中第二套系數(shù)的分辨率要優(yōu)于第一套系數(shù)，可以從圖3(e)和3(f)中看出，第二套系數(shù)網(wǎng)格移動(dòng)的更加劇烈，因此對(duì)間斷捕捉的分辨率更高. 因此為了提高分辨率，可以適當(dāng)?shù)恼{(diào)整偽弧長(zhǎng)控制函數(shù). 圖3(g)和3(h)均采用了220 和440 個(gè)網(wǎng)格數(shù)進(jìn)行對(duì)比，可以看出，PALM-5 在220 個(gè)網(wǎng)格數(shù)目下和WENO 在440 個(gè)網(wǎng)格數(shù)目下的近似分辨率幾乎完全相同，而從表2中可以看出采用PALM-5 算法在近似分辨率下需要的時(shí)間大大縮短.

表2 近似分辨率下WENO 和PALM-5 的時(shí)間對(duì)比Tab. 2 Time comparison of WENO and PALM-5 at approximate resolution

圖3 爆炸波問題的結(jié)果Fig. 3 Results of blast wave problem

5.3 一維激波擊打氣-氣問題

這是一個(gè)激波擊打氣-氣界面問題的多介質(zhì)算例，初值條件為

5.4 二維水下爆炸

從計(jì)算結(jié)果可以看出，兩種計(jì)算結(jié)果的對(duì)稱性都比較好，但是固定網(wǎng)格下捕捉到的界面擴(kuò)散嚴(yán)重，而偽弧長(zhǎng)算法下界面銳利清楚，且解的光滑性非常好，由于狀態(tài)方程以及物理量的巨大差異，因此固定網(wǎng)格下的結(jié)果在界面附近有輕微的數(shù)值震蕩，而采用偽弧長(zhǎng)算法計(jì)算的結(jié)果非常光滑，完全消除了數(shù)值震蕩,如圖5 所示. 從圖5(c)中可以看出網(wǎng)格的自適應(yīng)性非常好，而且網(wǎng)格足夠光滑，沒有出現(xiàn)間斷的網(wǎng)格點(diǎn).

圖5 水下爆炸結(jié)果Fig. 5 Results of two-dimensional underwater explosion problem

6 結(jié) 論

將偽弧長(zhǎng)算法和高階WENO 格式相結(jié)合發(fā)展了高精度偽弧長(zhǎng)算法，結(jié)合RGFM 和Level-Set 將其應(yīng)用到多介質(zhì)的計(jì)算中. 為了構(gòu)造高精度的數(shù)值格式，所有的計(jì)算均在坐標(biāo)變換之后的弧長(zhǎng)坐標(biāo)系下計(jì)算，此外，在多介質(zhì)的計(jì)算中，為了使的距離函數(shù)保持物理性質(zhì)，采用了法向速度的表達(dá)式，針對(duì)網(wǎng)格移動(dòng)以后距離函數(shù)的插值，給出了3 階的非守恒插值形式.根據(jù)數(shù)值算例，得出以下結(jié)論：

①經(jīng)過坐標(biāo)變換以后，高精度偽弧長(zhǎng)算法仍舊可以保持較高的精度；

②爆炸波算例說明高階偽弧長(zhǎng)算法比低階偽弧長(zhǎng)算法捕捉界面的分辨率更高，而且相對(duì)固定網(wǎng)格下結(jié)果提升顯著，在近似分辨率下偽弧長(zhǎng)算法需要的時(shí)間大大縮短，計(jì)算效率較高；

③多介質(zhì)算例說明，高精度的偽弧長(zhǎng)算法結(jié)合RGFM 可以成功的計(jì)算多介質(zhì)問題，并且有效降低了界面以及其他間斷處的數(shù)值震蕩，使得解的光滑性非常好，網(wǎng)格自適應(yīng)朝著間斷移動(dòng)，使得解捕捉間斷的能力非常強(qiáng). 解決了低階格式耗散強(qiáng)，高階格式震蕩強(qiáng)的問題. 發(fā)展了高分辨率，低震蕩的算法.