馬天寶，王晨濤，趙金慶，寧建國(guó)

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

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

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

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

1 偽弧長(zhǎng)算法

本文采用的偽弧長(zhǎng)算法包含3 個(gè)部分：第1 部分給出物理空間的控制方程在時(shí)間和空間上的離散格式；第2 部分是坐標(biāo)變換，將不均勻的物理空間變換到均勻的計(jì)算空間中進(jìn)行求解；第3 部分則是偽弧長(zhǎng)算法。

1.1 離散格式

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

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

對(duì)上式采用Green 公式進(jìn)行變換，則有：

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

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

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

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

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

1.2 坐標(biāo)變換

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

引入曲線坐標(biāo)

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

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

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

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

二維情況采用下述公式

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

因此方程(16)改寫(xiě)為

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

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

式中：?t為時(shí)間步長(zhǎng)。

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

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

綜上，偽弧長(zhǎng)算法的詳細(xì)步驟如下。

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

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

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

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

2 一維數(shù)值算例

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

2.1 一維Euler 方程激波管問(wèn)題

初值條件：

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

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

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

2.2 一維爆炸問(wèn)題

初值條件：

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

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

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

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

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

3 數(shù)值精度驗(yàn)證

3.1 一維精度驗(yàn)證

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

3.2 二維精度驗(yàn)證

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

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

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

表2 有限體積法與偽弧長(zhǎng)算法在不同網(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)流的計(jì)算

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

一維爆炸沖擊問(wèn)題，控制方程取為如下：

初值條件：

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

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

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

4.2 二維化學(xué)反應(yīng)流前臺(tái)階爆轟問(wèn)題

計(jì)算域?yàn)?[0,3]×[0,1]。 障礙物的區(qū)域?yàn)?[1,3]×[0,0.4]，偽弧長(zhǎng)控制函數(shù)取為 (a1,a2)=(2,0.2)，初值條件取ZND[10]的解析解作為初值?？刂品匠滩捎檬?27)兩步化學(xué)反應(yīng)流。邊界條件：除了左右兩側(cè)，其他均設(shè)為反射邊界。左右兩側(cè)為流入流出邊界。終止時(shí)間T=0.22。圖6 為各方法給出的T=0.22 時(shí)刻的狀態(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 有限體積進(jìn)行計(jì)算的結(jié)果。其余計(jì)算均在 360×90網(wǎng)格下進(jìn)行計(jì)算。通過(guò)對(duì)比可以發(fā)現(xiàn)偽弧長(zhǎng)算法都較為準(zhǔn)確的捕捉到了爆轟反應(yīng)的細(xì)節(jié)，對(duì)爆炸沖擊波陣面都得到了較為精確的捕捉，而且用偽弧長(zhǎng)五階算法的結(jié)果和參照解的結(jié)果更為接近，并且在局部拐角區(qū)域，高階算法對(duì)細(xì)節(jié)的捕捉更為精確。這說(shuō)明了高階偽弧長(zhǎng)算法收斂速度較快，因此可以更好地處理爆炸與沖擊強(qiáng)間斷的問(wèn)題。

5 結(jié) 語(yǔ)

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

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

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