王晨濤，馬天寶，李坤

（北京理工大學 爆炸科學與技術國家重點實驗室，北京 100081）

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

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

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

1 控制方程

1.1 物理空間下的控制方程

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

1.2 弧長空間下的控制方程

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

2 偽弧長算法

2.1 偽弧長算法幾何意義

為了削弱方程的奇異性，通過引入一個弧長約束方程，即

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

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

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

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

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

2.2 偽弧長算法網格移動

偽弧長算法的網格滿足一個能量極值函數：

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

至此，網格便可以采用數值方法迭代計算，例如采用收斂較快的高斯-賽德爾迭代，格式如下.

3 數值格式

3.1 控制方程的離散

3.2 守恒變量的更新

網格移動完以后，需要更新在新網格上的物理量，本文采用守恒插值[6]來計算：

4 Level-Set 技術

4.1 控制方程的離散

4.2 距離函數的非守恒插值

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

由于計算過程是在弧長坐標系下，因此將上式映射至弧長空間，可得：

4.3 多介質流高精度偽弧長算法流程

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

⑥判定計算時間是否到達. 若到達，則終止計算，否則轉到步驟①.

5 數值算例

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

5.1 精度測試

用一個包含精確解的一維的Euler 方程來驗證高階偽弧長算法PALM-5 的收斂精度. 初值條件和精確解為

從表1 的計算結果可以看出，弧長坐標系下的高階偽弧長算法保證了較高的精度，沒有因為網格的變形而損失精度. 這說明本文的算法滿足了高精度的性質.

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

5.2 一維爆炸波

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

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

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

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

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

5.4 二維水下爆炸

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

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

6 結 論

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

①經過坐標變換以后，高精度偽弧長算法仍舊可以保持較高的精度；

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

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