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

（北京理工大學(xué)爆炸科學(xué)與技術(shù)國家重點實驗室, 北京 100081）

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

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

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

1 偽弧長算法

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

1.1 離散格式

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

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

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

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

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

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

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

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

1.2 坐標變換

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

引入曲線坐標

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

1.3 網(wǎng)格移動算法

偽弧長控制函數(shù)取為

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

二維情況采用下述公式

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

因此方程(16)改寫為

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

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

式中：?t為時間步長。

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

圖1 新舊網(wǎng)格轉(zhuǎn)化示意圖Fig. 1 Diagram of old grid transforming to new grid

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

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

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

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

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

2 一維數(shù)值算例

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

2.1 一維Euler 方程激波管問題

初值條件：

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

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

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

2.2 一維爆炸問題

初值條件：

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

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

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

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

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

3 數(shù)值精度驗證

3.1 一維精度驗證

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

3.2 二維精度驗證

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

采取構(gòu)造一組人為解進行數(shù)值驗證。相應(yīng)的源項需要根據(jù)人為解進行修正。其中構(gòu)造的人為解詳細如下：

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

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

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

4 爆轟化學(xué)反應(yīng)流的計算

4.1 一維化學(xué)反應(yīng)流

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

初值條件：

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

圖5 一維化學(xué)反應(yīng)流密度Fig. 5 Density of one dimensional chemical reaction flow

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

4.2 二維化學(xué)反應(yīng)流前臺階爆轟問題

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

圖6 一維化學(xué)反應(yīng)流壓力Fig. 6 Pressure of one dimensional chemical reaction flow

圖7 一維化學(xué)反應(yīng)流的網(wǎng)格軌跡Fig. 7 Mesh trajectory of one dimensional chemical reaction flow

圖8 二維化學(xué)反應(yīng)流圖Fig. 8 Results of two dimensional chemical reaction flow

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

5 結(jié) 語

本文將高精度數(shù)值算法和偽弧長算法進行結(jié)合，較為詳細的介紹了高階偽弧長算法從推導(dǎo)到應(yīng)用的過程。針對高階偽弧長算法在處理通量計算和網(wǎng)格變換之后的物理量重構(gòu)中的難題，通過引入計算坐標，經(jīng)過坐標變換將計算過程轉(zhuǎn)換到曲線坐標系下進行計算。通過構(gòu)造的人為解，驗證了高階偽弧長算法的數(shù)值精度。

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

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