爆炸沖擊波強間斷問題的高階偽弧長算法研究*

馬天寶，王晨濤，趙金慶，寧建國

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

雙曲守恒律方程有著廣泛的應用，在內爆動力學，慣性約束聚變，計算天文學，磁流體力學，計算爆炸力學，計算生物學等領域都和它密切相關。它有一個明顯的特征，無論初值條件多么的光滑，隨著時間的演化，最終都可能會演化為帶有強間斷的解。爆炸與沖擊問題是強非線性的雙曲方程，這些方程的解通常帶有奇異性，也就是在解的附近存在較大的變化，包括激波，稀疏波，接觸間斷等。為了克服上述難題，一些高精度的數值格式被提出，如MUSCL、WENO[1-2]等。但采用高階重構時，解在間斷附近會產生虛假震蕩。因此構造有效的高精度數值算法，不僅要求能捕捉到間斷區(qū)域，同時又可以消除虛假震蕩。針對虛假震蕩的非物理現象，可以考慮將原始守恒方程映射到特征空間上[3]，在特征空間上求得數值通量之后[4]，再將通量映射回原始空間。這種通過對方程組多變量的耦合進行解耦的過程，可以巧妙地避開虛假震蕩。

對于常規(guī)的有限體積方法，計算網格是均勻網格，若提高計算精度，那么相應的策略是增加網格數目。但是對于物理量梯度變化較小的區(qū)域，粗糙的網格便可以得到滿意的結果，對光滑區(qū)域求解增加的計算量是沒有必要的，因此這無疑在一定程度上降低了計算效率?；诖耍鶕獾男再|，對網格進行合理的節(jié)點分布，對于高效、精確地計算雙曲守恒方程有著極其重要的作用。偽弧長算法是解決這種問題的有效方法之一。在偽弧長算法中，網格節(jié)點會隨著時間的變化而變化，在保證網格節(jié)點數目不變的情況下，在物理量梯度較大的區(qū)域自動加密，而在解較為平滑的區(qū)域自動稀疏。Tang 等[5]通過求解一個基于變分原理的網格方程得到了一種移動網格方法，2006 年Tang[6]將移動網格算法拓展到二維情況，近年來，自適應網格算法得到了較大的發(fā)展[7-9]。

本文將高精度的WENO 格式同偽弧長算法相結合，針對網格移動之后造成的網格非均勻和非正交使得通量的計算和網格物理量重構的過程較為復雜的現象，提出采用坐標變換，將坐標映射到均勻正交的計算坐標系下進行計算。結果表明這樣可以提高數值精度。通過一維激波管問題結果的比較，可以看出偽弧長算法相較于傳統(tǒng)的固定網格算法可以更好地捕捉到強間端，而高階偽弧長算法對間斷捕捉的分辨率最高。結果表明高階偽弧長算法相較于有限體積以及低階格式的偽弧長算法，它更適合處理處理爆炸與沖擊強間斷問題。此外，可以改變弧長參數使得網格移動更加明顯，進而提高對間斷捕捉的分辨率。

1 偽弧長算法

本文采用的偽弧長算法包含3 個部分：第1 部分給出物理空間的控制方程在時間和空間上的離散格式；第2 部分是坐標變換，將不均勻的物理空間變換到均勻的計算空間中進行求解；第3 部分則是偽弧長算法。

1.1 離散格式

考慮如下的雙曲守恒系統(tǒng)：

對式(1)在網格 單元上進行積分，得到：

對上式采用Green 公式進行變換，則有：

數值通量采用局部Lax-Friedrichs 格式，詳細定義如下：

式(3) 可以寫成半離散格式：

以下僅對x方向進行離散，y方向僅需要將變量i替換為j即可。

對于時間的離散，忽略掉網格的索引i,j，則采用如下的四階TVD-RK 格式：

式(8)和式(9)的重構中，重構的物理量可以選擇原始變量，守恒變量以及特征變量。由于歐拉格式的非線性性質，采用高階格式容易引起非物理震蕩，因此可以考慮采用將變量映射到特征空間中，在特征空間中進行重構，重構完成之后再映射回原始空間，具體計算如下：

1.2 坐標變換

在偽弧長算法中，移動之后的網格是非均勻網格，而傳統(tǒng)的數值格式均是基于均勻網格給出的。通過限制網格的移動距離，可以近似用傳統(tǒng)的格式進行插值運算，然而僅僅在一階或者二階情況下，這種插值引起的誤差不大。若直接運用到高階格式，將會引起較大的誤差。因此，對于高階偽弧長算法的通量計算，必須考慮網格尺度的影響。為了解決這種問題，通常有兩種策略，第一種策略是將網格尺度引入WENO 格式的非線性權重，然后直接進行構造高階格式[11-13]。這種構造清晰明了，但是格式較為復雜，而且不易推廣到二維或者三維情況。另外一種策略是轉換坐標系，將當前坐標系變換到曲線坐標系下，因為曲線坐標系中計算網格是均勻正交的，因此可以直接用傳統(tǒng)的格式進行通量計算。詳細變換如下。

引入曲線坐標

則式(1)可以由原始坐標轉換到曲線坐標下進行計算，即

1.3 網格移動算法

偽弧長控制函數取為

為了保證網格的光滑性，可以用低通濾波器進行光滑處理，一維情況表達式為

二維情況采用下述公式

具體光滑處理次數應該視情況而定，通常情況下處理2～3 次即可，若非線性情況較嚴重，可以處理10 次左右。

因此方程(16)改寫為

通過對式(20)的計算，可以求得網格的分布。式(20)的計算可以采用Jacobian 迭代進行計算，計算步驟如下

因此根據式(21)可以得到式(13)的網格移動速度為

式中：?t為時間步長。

網格移動完成之后，接下來的工作就是需要得到物理量在新網格上的值，這一重要步驟又稱物理量重映。常用的物理量重映方法分為2 類，一是基于通量的物理量重映方法，二是基于網格相交的精確積分守恒重映方法?？紤]到在相鄰的時間間隔要求，網格移動距離較小，而且網格始終處于保凸性，因此本文采用的映射方法為第一種。

圖1 新舊網格轉化示意圖Fig. 1 Diagram of old grid transforming to new grid

綜上，偽弧長算法的詳細步驟如下。

第2 步，求解網格控制函數式(17)，并用式(18)和式(19)進行光滑打磨控制函數。

第4 步，物理量映射，根據式(23)更新新網格下的物理量。

第5 步，求解物理量控制方程：(a)如果采用低階偽弧長算法，則直接在物理空間進行求解。通過等式(10)更新時刻的物理量；(b)如果采用高階偽弧長算法，則利用式(12)、式(13)和式(14)將物理空間的物理量轉換到均勻計算坐標系下進行求解。通過式(10)更新時刻的物理量，求解完成之后根據逆變換再把物理量映射回物理空間。

第6 步，如果求解達到終點時間，則程序終止，否則轉到第2 步。

2 一維數值算例

下面基于偽弧長數值算法進行一些算例測試，由于上面的理論是基于二維情況推導的，對于一維情況直接進行降維處理，然后再計算即可。限于篇幅這里不再詳細說明一維的理論。此外，以下所有數值算例均采用無量綱進行計算，偽弧長控制函數均采用式(17)計算。

2.1 一維Euler 方程激波管問題

初值條件：

計算域為[?5，5]，終止時間為tmax=，偽弧長控制函數(17)取為=(1, 5)，邊界條件為自由輸出邊界條件，計算結果如圖2 所示。

上面的計算均是在150 個網格數目下求得的，其中參照解是采用20 000 個固定網格下結合五階特征變換進行計算的。其中MUSCL,WENO-5,WENOCV-5 分別代表的是固定網格下用二階MUSCL，五階WENO，五階特征變量WENO 格式計算；而PALM-2, PALM-5 代表的是在采用二階偽弧長和五階偽弧長的格式計算。

從圖2(a)中可以看出在固定網格下，高階格式相對于低階格式可以更加精確地捕捉到間斷，但是從局部放大中可以看出，五階固定網格在間斷附近產生了數值震蕩，這是高階格式帶來的弊端，因此可以通過特征變換進行求解，進行特征變換之后，高階格式消除了局部震蕩，整體解的光滑性較好，因此可以經過特征變換得到較為滿意的結果。從圖2(b)中可以看出，偽弧長算法可以很容易的捕捉到間斷面，而且二階偽弧長算法對間斷的捕捉幾乎接近于固定網格下五階的捕捉程度，但是弱于偽弧長五階算法。因此這說明了高階偽弧長算法是有效的。圖2(c)是網格的運動軌跡，說明偽弧長算法很好地控制了網格的自適應移動，使得網格在梯度較大的區(qū)域進行加密，進而捕捉到間斷解。

2.2 一維爆炸問題

初值條件：

該算例的計算域為[0,1]，終止時間tmax=0.038，兩套偽弧長控制函數分別取為(a1,a20,20)，(a1,a2)2=(10,80)，邊界條件為固壁反射條件，計算結果如圖3 所示，其中2-1, 2-2, 5-1, 5-2 分別代表二階和五階分別采用系數和系數時的計算結果，可以看出，當采用不同的系數時，逼近的分辨率也不同。

圖3 一維爆炸算例計算結果Fig. 3 Results of one dimensional explosion example

該算例是一個爆炸問題的算例，是檢驗高精度格式的經典算例。普通低階格式在網格數目較少的情況下很難捕捉到強間斷，而且在計算的過程中，由于稀疏波的傳播，壓力接近于零，使得計算容易產生負值，導致程序終止。因此，在計算的時候均采用了保正性條件，詳細的保正性理論可以參考文獻[9]。圖中的參照解是在固定網格下用五階特征變換在網格數為25 000 下進行計算得到的。

從上面的計算結果可以看出，二階固定網格耗散較強，對最高點的捕捉能力較弱；五階固定網格較二階提升很大，可以較為清晰的捕捉到多個間斷點，而且在間斷點附近耗散性非常??；二階偽弧長算法較固定網格算法有一定提升，固定網格對斷點的捕捉能力很弱，耗散性較強，而采用偽弧長算法以后，耗散能力減弱，對間斷點的捕捉能力增強；采用五階偽弧長算法最為滿意，各個間斷點都可以較為準確的捕捉，而且在間斷點鄰域的結果和參照解也最為接近。

通過網格的軌跡線圖，可以看出在t=0.027 時刻附近，兩個沖擊波相遇，同時伴隨著稀疏波的傳播，正是因為稀疏波的傳播，造成壓力過小，接近于零。這樣，采用普通的通量差分格式計算會因為計算機的浮點誤差和計算格式的數值誤差造成物理量為負，從而導致程序的終止。因此需要采用特殊的格式進行修正才能保證程序的進行，本文采用了保正性格式進行限制。

3 數值精度驗證

3.1 一維精度驗證

表1 有限體積法與偽弧長算法在不同網格數下的誤差和精度Table 1 Numerical errors and precision of FVM and PALM changing with grid numbers (example 3)

3.2 二維精度驗證

采用二維兩步化學反應流來測試偽弧長算法的范數誤差和收斂階。對應式(1)的詳細表達式如下：

采取構造一組人為解進行數值驗證。相應的源項需要根據人為解進行修正。其中構造的人為解詳細如下：

表2 給出了二維驗證計算結果，可以看出，隨著計算網格的加密，誤差在逐漸減小，而且收斂階和算法本身較為吻合。通過圖4 可以看出，高階有限體積算法對等密度圖的捕捉較低階有限體積更為精細，而且采用偽弧長算法之后，網格自動在物理量梯度較大的區(qū)域進行自適應加密，同時也沒有失去解的準確性。

表2 有限體積法與偽弧長算法在不同網格數下的誤差和精度Table 2 Numerical errors and precision of FVM and PALM changing with grid numbers (Example 4)

圖4 密度云圖（T =2π，網格40×40）Fig. 4 Density contours (T =2π,grid 40×40)

4 爆轟化學反應流的計算

4.1 一維化學反應流

一維爆炸沖擊問題，控制方程取為如下：

初值條件：

偽弧長控制函數取為 (a1,a2)=(0.5,5)，邊界條件為反射邊界，終止時間T=0.1。計算結果如圖5 所示。

圖5 一維化學反應流密度Fig. 5 Density of one dimensional chemical reaction flow

一維化學反應流的參照解是基于5 000 個固定網格結合特征變換進行計算的，而二階格式采用的是120 個網格，五階格式采用的是120 個網格?？梢钥闯觯诠潭ňW格下，五階格式對間斷的捕捉要優(yōu)于二階格式，而經過偽弧長的變化之后，二階格式也可以較為準確是捕捉到強間斷。對解捕捉的最好結果是采用了五階偽弧長格式，它對網格捕捉的分辨率最高，圖中可以看出高階偽弧長和參照解幾乎重合。這充分說明了高階偽弧長格式的優(yōu)越性。在本算例中，由于壓力的初值變化跨越了10 個數量級以上，計算時采用的變量均為原始守恒變量，并沒有采用局部特征變換，從而導致在間斷附近存在微小震蕩。采用了偽弧長算法，并沒有完全消除震蕩，但是卻很明顯的減弱了數值震蕩，這說明偽弧長算法對于減弱方程的奇異性有著很好的效果。需要注意的是，由于爆炸沖擊問題的強烈的奇異性，常規(guī)的高階格式會因為計算物理量為負，而導致程序終止。在本文中，修正式(7)的 ε=10?15，并結合保正性算法[9]，這樣才能使得程序可以繼續(xù)計算。此外可以看到，網格的的軌跡得到了很好的自適應，它在解梯 度較大的區(qū)域自適應加密。

4.2 二維化學反應流前臺階爆轟問題

計算域為 [0,3]×[0,1]。 障礙物的區(qū)域為 [1,3]×[0,0.4]，偽弧長控制函數取為 (a1,a2)=(2,0.2)，初值條件取ZND[10]的解析解作為初值?？刂品匠滩捎檬?27)兩步化學反應流。邊界條件：除了左右兩側，其他均設為反射邊界。左右兩側為流入流出邊界。終止時間T=0.22。圖6 為各方法給出的T=0.22 時刻的狀態(tài)。

圖6 一維化學反應流壓力Fig. 6 Pressure of one dimensional chemical reaction flow

圖7 一維化學反應流的網格軌跡Fig. 7 Mesh trajectory of one dimensional chemical reaction flow

圖8 二維化學反應流圖Fig. 8 Results of two dimensional chemical reaction flow

圖6 的參照解是采用了1200×400網格數結合五階WENO 有限體積進行計算的結果。其余計算均在 360×90網格下進行計算。通過對比可以發(fā)現偽弧長算法都較為準確的捕捉到了爆轟反應的細節(jié)，對爆炸沖擊波陣面都得到了較為精確的捕捉，而且用偽弧長五階算法的結果和參照解的結果更為接近，并且在局部拐角區(qū)域，高階算法對細節(jié)的捕捉更為精確。這說明了高階偽弧長算法收斂速度較快，因此可以更好地處理爆炸與沖擊強間斷的問題。

5 結 語

本文將高精度數值算法和偽弧長算法進行結合，較為詳細的介紹了高階偽弧長算法從推導到應用的過程。針對高階偽弧長算法在處理通量計算和網格變換之后的物理量重構中的難題，通過引入計算坐標，經過坐標變換將計算過程轉換到曲線坐標系下進行計算。通過構造的人為解，驗證了高階偽弧長算法的數值精度。

各算例的計算驗證表明：偽弧長算法相較于有限體積算法可以更成功地捕捉到強間斷；而且高階偽弧長算法對強間斷的捕捉相對于傳統(tǒng)高階有限體積法和二階偽弧長算法有更大的優(yōu)勢，它可以以更少的網格數目獲得更高的分辨率。

本文采用偽弧長算法結合高階WENO 有限體積格式進行計算，有限體積格式在處理物理量重映過程中的過程較為繁瑣，而且對網格移動的限制較為嚴格，要求每兩步網格迭代變化的距離不能太大。因此為了提高計算效率，可以嘗試采取有限差分進行計算，有限差分結合動網格和坐標映射可以避免物理量的重映過程，因此可以適當提高計算效率。為了提高計算精度獲得高分辨率的結果，而且又不需要很高的網格數目，可以采用高階偽弧長算法進行計算，它可以以較少的網格數目得到較為滿意的結果，從而提高計算效率。