爆轟波強間斷問題的偽弧長算法及其人為解驗證

馬天寶，陳建良，寧建國，原新鵬

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

如何能夠精確捕捉和追蹤爆轟波和沖擊波波陣面的傳播過程一直是爆炸與沖擊問題數(shù)值模擬的難點。爆炸與沖擊問題通常由Euler方程組進行描述，求解這些方程組得到的解具有奇異性，包括稀疏波、激波間斷等。為了解決這些氣相爆轟問題，人們提出了很多數(shù)值格式，如TVD格式，ENO格式，WENO[1]格式等。中國科研人員張德良等[2]、王成等[3]、王昌建等[4]在氣相爆轟領(lǐng)域取得了豐碩的成果。

在數(shù)值模擬氣相爆轟問題時，巨大的計算量是經(jīng)常遇到的挑戰(zhàn)，自適應(yīng)網(wǎng)格方法是有效的解決途徑之一。它通過合理地分布網(wǎng)格，高效精確地捕捉到強間斷，極大地節(jié)省了計算資源。自適應(yīng)網(wǎng)格包括三類典型的方法：p方法、h方法和r方法。p方法根據(jù)一定的誤差評估方法或指標(biāo)改變插值多項式近似階數(shù)。h方法通過增加或刪除網(wǎng)格節(jié)點來生成新網(wǎng)格，例如將網(wǎng)格點添加到數(shù)值梯度很大的區(qū)域，或者在結(jié)果很光滑的區(qū)域移除節(jié)點。r方法中通過改變節(jié)點的位置達到反映結(jié)構(gòu)特征的目的，而不改變節(jié)點總數(shù)目。文獻[5]中研究了基于自適應(yīng)網(wǎng)格r方法的偽弧長算法并應(yīng)用于模擬氣相爆轟問題。

雖然國內(nèi)外科研人員在數(shù)值模擬氣相爆轟問題的數(shù)值格式上取得了很多成果，但對相關(guān)工作的驗證與確認(rèn)卻較少。驗證與確認(rèn)的概念有美國計算機協(xié)會在1979年首次正式提出，并得到歐美國家的高度重視。在高性能計算機支撐下的大規(guī)模工程仿真，其可靠性一直是亟待解決的重難點，細(xì)微的偏差都可能對最終結(jié)果造成不可估量的損失。在美國開展的“先進模擬與計算計劃”中，為了模擬預(yù)測庫存核武器的安全性和可靠性，各大實驗室投入大量人力財力促進驗證與確認(rèn)的發(fā)展。在處理復(fù)雜的工程、物理問題時，人們通常會對所研究的問題適當(dāng)簡化然后進行數(shù)值模擬，但模型簡化、邊界條件的近似等過程都會對數(shù)值模擬結(jié)果帶來一定的誤差，而程序邏輯結(jié)構(gòu)、迭代格式的選取等因素也會對模擬結(jié)果造成影響[6]。對數(shù)值程序進行驗證與確認(rèn)，就是通過一定的方法衡量模擬結(jié)果與真實的物理現(xiàn)象的相關(guān)程度。

Roache[7]在文章中提出關(guān)于數(shù)值結(jié)果精確度控制的相關(guān)問題，隨后發(fā)展為不確定度的量化，并包含程序驗證、計算驗證和確認(rèn)三部分。Oberkampf等[8]深入探討了關(guān)于人為解、解析解、高精度數(shù)值解的驗證與確認(rèn)基準(zhǔn)建立的問題。Roy[9]總結(jié)了人為解方法的實現(xiàn)步驟，并指出解驗證包括截斷誤差、迭代收斂誤差和離散誤差的評估。鄧小剛等[10]梳理了計算流體力學(xué)中驗證與確認(rèn)的原理、方法，并用實例說明經(jīng)過驗證和確認(rèn)的模型才具有高可信度。王瑞利等[11]對二維流體力學(xué)方程組構(gòu)造了一組普適的人為解，并應(yīng)用PPM程序的數(shù)值結(jié)果驗證其正確性。

本文中，重點研究使用偽弧長算法處理爆轟波的強間斷奇異性問題，將物理空間的坐標(biāo)和物理量推廣到弧長空間，計算空間的轉(zhuǎn)換會巧妙地繞過物理空間的奇異點；針對偽弧長算法編寫了二維程序，用人為解方法對二維程序進行驗證，通過計算數(shù)值解的精度，驗證程序的可靠性；最后將程序應(yīng)用于爆轟波傳播等物理問題的數(shù)值模擬，模擬結(jié)果顯示該方法能較好地捕捉和追蹤爆轟波波陣面的傳播。

1 偽弧長算法

本文采用的偽弧長算法包含2個部分：第1部分是有限體積法，給出物理空間控制方程在時間和空間的離散方式；第2部分是偽弧長移動算法，具體說明物理空間的坐標(biāo)和物理量轉(zhuǎn)換到弧長空間的方法。

1.1 數(shù)值離散方法

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

(1)

式中：w為任意物理量，F(xiàn)(w)、G(w)、S(w)為w的函數(shù)。

?Ki,jS(w)dσ

(2)

式中：dσ為面積微元。

利用Green-Gauss定理，可以將上式轉(zhuǎn)化為：

(3)

式中：|Ki,j|為網(wǎng)格Ki,j的面積；?Ki,j為網(wǎng)格Ki,j的邊界；ds為邊界長度微元，ni,j是邊界?Ki,j的單位外法向量，δ(w)=(F(w),G(w))。定義Lax-Friedrichs通量格式為：

(4)

式中：u、v為任意物理量，n為單位法向量，c=max(δ′(u)·n,δ′(v)·n)。式(3)可以寫成半離散格式：

(5)

對時間的離散，采用三階Runge-Kutta格式：

(6)

1.2 偽弧長移動算法

首先說明偽弧長算法在計算空間的轉(zhuǎn)換關(guān)系。物理空間Ωp的坐標(biāo)采用x=(x,y)表示，弧長空間Ωc的坐標(biāo)使用ζ=(ξ,η)來表示，〈xi-1,j-1,xi-1,j,xi,j,xi,j-1〉為網(wǎng)格Ki,j的4個頂點。將物理空間Ωp和弧長空間Ωc坐標(biāo)的一一對應(yīng)關(guān)系表示為(x,y)=(x(ξ,η),y(ξ,η))。通過變分方法，可知需要求得泛函E(x,y)的最小值[12]，E(x,y)表達式為：

(7)

式中：G1和G2是控制函數(shù)，=(?ξ,?η)T，?ξ=?/?ξ，?η=?/?η。則相應(yīng)的Euler-Lagrange方程是：

?ξ(G1?ξx)+?η(G1?ηx)=0, ?ξ(G2?ξy)+?η(G2?ηy)=0

(8)

定義偽弧長函數(shù)為：

(9)

通過調(diào)整弧長參數(shù)α1、α2的大小，可以控制網(wǎng)格的加密區(qū)域，通常ψ>0。

將式(8)中G1和G2都取為ψ，可以得到：

·(ψx)=0,·(ψy)=0

(10)

在計算中，使用Gauss-Seidel方法進行迭代。對上面 (8) 式，可以寫成下面形式：

(11)

(12)

式中：1≤i≤Nξ，1≤j≤Nη，Nξ,Nη為偽弧長空間ξ,η方向的網(wǎng)格總數(shù)；上角標(biāo)[z]和[z+1]分別表示舊網(wǎng)格和新網(wǎng)格。

基于Han等[13]的研究工作，采用幾何插值方法，對新網(wǎng)格x[z+1]和舊網(wǎng)格x[z]之間物理量進行更新，使新舊網(wǎng)格保證通量守恒：

(13)

(14)

(15)

式中：m、n為任意函數(shù)，D1由式(16)計算，D2、D3、D4以此類推：

(16)

結(jié)合前文算法應(yīng)用，將偽弧長算法應(yīng)用在二維雙曲守恒系統(tǒng)，計算步驟總結(jié)如下。

(2)求解弧長空間控制方程：

2 程序驗證

通常對程序和解的驗證方法有精確解方法、人為解方法、軟件質(zhì)量保證方法、標(biāo)準(zhǔn)數(shù)值解對照方法、專家判斷法和代碼對比等方法[14]，其中精確解方法和人為解方法是最常用的驗證方法。但描述物理模型大多是非線性偏微分方程(組)，能夠得到精確解的方程類型極其有限，特別是對于二維及多維方程組，人為解相對更容易獲取，因此人為解方法在程序驗證領(lǐng)域的應(yīng)用更為廣泛。本文中，采用人為解方法對二維程序進行驗證。

人為解構(gòu)造方法的基本思路是：針對需要求解的偏微分方程(組)，先假設(shè)一個(組)可達解，將這個(組)解帶入原方程(組)，逆向推導(dǎo)出使方程(組)成立所必須添加的源項，并且確定對應(yīng)的邊界條件，然后調(diào)整源項和邊界等對應(yīng)部分的程序代碼，得到相應(yīng)的數(shù)值解[6]。

通常，為了更好地使用人為解方法驗證程序，人為解的選取需要遵循以下規(guī)定[14]：(1)人為解應(yīng)當(dāng)由光滑解析函數(shù)組成，例如多項式、指數(shù)函數(shù)、三角函數(shù)或者它們的組合，這樣方便求導(dǎo)運算，對精度保證也極其重要；(2)人為解能夠使方程組的每一項都成立，即應(yīng)具備普適性；(3)根據(jù)具體的控制方程組，人為解應(yīng)當(dāng)有一定數(shù)量的非平凡導(dǎo)數(shù)；(4)人為解在求解區(qū)域不能具有奇異性；(5)人為解應(yīng)當(dāng)滿足控制方程的基本假設(shè)，例如如果程序需要保正性，那選取的人為解應(yīng)當(dāng)滿足非負(fù)性。

針對式(1)給出的二維雙曲系統(tǒng)，令：

(17)

狀態(tài)方程為：

(18)

(19)

式中：k1為指前因子，Ea為單位質(zhì)量反應(yīng)物的活化能，R為氣體常數(shù)，T0=p/(ρR)為溫度。計算中取k1=2566.4，Ea=50，γ=1.2。

根據(jù)人為解構(gòu)造準(zhǔn)則，我們構(gòu)造出上述方程組的一組人為解(所有物理量均為量綱一形式)：

(20)

這樣選取的人為解形式簡單，有足夠階的非無效導(dǎo)數(shù)，而且可以使方程組源項為0，程序更改方便。

計算域取為[-1,1]×[-1,1]，采用周期性邊界條件，網(wǎng)格數(shù)分別取為(40×40)，(80×80)，(160×160)，(320×320)，計算終止時間取為T=1.0以保證經(jīng)歷一個完整的周期。

收斂階計算公式為[15]：

(21)

式中：Δsk-1和Δsk為網(wǎng)格步長，εk-1和εk分別為網(wǎng)格步長為Δsk-1和Δsk時數(shù)值解與精確解之間的L1范數(shù),用以表征兩者之間的額誤差。范數(shù)ε具體的表達式為[26]：

(22)

通過計算得到直接的有限體積法以及偽弧長算法的密度的數(shù)值解誤差ε及精度O如表1所示。在T=0.48時刻直接有限體積法和偽弧長算法的網(wǎng)格分布示意圖和密度云圖分別如圖2所示，圖中顯示的是使用40×40個網(wǎng)格的情況。通過對比可以看出，格式計算精度為2階，而且當(dāng)使用偽弧長算法時，程序精度得到適當(dāng)提高，因為偽弧長算法能夠使網(wǎng)格在數(shù)值梯度較大的區(qū)域自適應(yīng)加密，使網(wǎng)格點分布與物理解耦合進而提高解的精度和分辨率。

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

3 氣相爆轟波傳播問題的數(shù)值模擬

模擬二維管道中氣相爆轟波傳播問題(所有物理量均為無量綱形式)，取計算域為[0,300]×[0,50]，測試偽弧長算法對爆轟波陣面的捕捉效果。固定網(wǎng)格的有限體積法初始網(wǎng)格數(shù)設(shè)定1 500×250，偽弧長算法網(wǎng)格數(shù)取為750×250。選取y=25.0的切面，研究網(wǎng)格節(jié)點沿x方向的分布情況。圖3為T=2.0時密度與壓力的計算結(jié)果?？梢钥闯?，兩種方法的波陣面基本重合，也就是說雖然偽弧長算法在x方向的網(wǎng)格數(shù)只有有限體積法的一半，但波陣面?zhèn)鞑ノ恢帽３忠恢拢以诓嚸嫣?，偽弧長算法由于網(wǎng)格的自適應(yīng)移動，依然可以高效地捕捉到間斷面。圖4為有限體積法和偽弧長算法在波陣面附近壓力云圖的局部放大圖。對比結(jié)果說明偽弧長算法網(wǎng)格分布會隨波陣面的傳播而移動，在波陣面附近網(wǎng)格分布和有限體積法網(wǎng)格分布效果相同，在未反應(yīng)區(qū)和已反應(yīng)區(qū)網(wǎng)格分布相對稀疏，這說明偽弧長算法用有限體積法一半數(shù)目的網(wǎng)格依然可以捕捉和追蹤爆轟波陣面的傳播，同時有效減小計算量。

為進一步驗證程序?qū)庀啾Z問題的適應(yīng)性，選擇采用Ar稀釋、氫氣、氧氣體積分?jǐn)?shù)之和為70%的氫氧混合氣體，采用一維定常解作為初始條件，對二維管道中的氣相爆轟問題進行數(shù)值模擬研究。在波陣面前方設(shè)置密度擾動，密度擾動并不會影響胞格爆轟的三波點數(shù)的多少也不會影響胞格的形狀，只會加快胞格的形成時間，提高計算效率。計算中，取q=19.713 2，γ=1.44，CFL系數(shù)為0.8，入射爆轟波的馬赫數(shù)Ma=5.6，管道縱向為100個單位長度，橫向初始賦值占100個單位長度。圖5給出了模擬結(jié)果圖像。從圖5(a)看出當(dāng)爆轟波傳播到約x=360時開始形成完整的胞格結(jié)構(gòu)，但這些胞格是因為擾動才快速形成的。隨著爆轟波繼續(xù)推進，部分三波點的壓力值降低，而另外一些壓力則增大，并最終形成圖5(b)所示的規(guī)則穩(wěn)定的胞格結(jié)構(gòu)，這與文獻[17]的計算結(jié)果保持一致。

4 結(jié) 論

將偽弧長算法應(yīng)用在處理爆轟波的強間斷問題，說明將計算空間由物理空間轉(zhuǎn)換到弧長空間的方法，詳細(xì)說明了程序的人為解構(gòu)造過程，對比分析偽弧長算法對計算精度的提高，驗證偽弧長算法求解非線性偏微分方程組的有效性；通過對比分析氣相爆轟波傳播過程中波陣面處網(wǎng)格的自適應(yīng)加密，說明偽弧長算法相比于直接有限體積方法在捕捉波陣面效果上的優(yōu)勢，計算過程中可以實現(xiàn)在保證精度的前提下適當(dāng)減少網(wǎng)格數(shù)量，從而提高求解效率。

[1] LIU X D, OSHER S, CHAN T. Weighted essentially non-oscillatory schemes[J]. Journal of Computational Physics, 1994,115(1):200-212.

[2] 張德良,謝巍,郭長銘,等.氣相爆轟胞格結(jié)構(gòu)和馬赫反射數(shù)值模擬[J].爆炸與沖擊,2001,21(3):161-167.

ZHANG Deliang, XIE Wei, GUO Changming, et al. Numerical simulation of cellar structures and machreflection of gaseous detonation waves[J]. Explosion and Shock Waves, 2001,21(3):161-167.

[3] 王成, SHU C-W. 爆炸力學(xué)高精度數(shù)值模擬研究進展[J].科學(xué)通報,2015(10):882-898.

WANG Cheng, SHU C-W. Progress in high-resolution numerical simulation of explosion mechanics[J]. Chinese Science Bulletin, 2015(10):882-898.

[4] 王昌建,徐勝利.直管內(nèi)胞格爆轟的基元反應(yīng)數(shù)值研究[J].爆炸與沖擊,2005,25(5):405-416.

WANG Changjian, XU Shengli. Numerical study on cellular detonation in a straight tube based on detailed chemical reaction model[J]. Explosion and Shock Waves, 2005,25(5):405-416.

[5] 王星,馬天寶,寧建國.雙曲偏微分方程的局部偽弧長方法研究[J].計算力學(xué)學(xué)報,2014(3):384-389.

WANG Xing, MA Tianbao, NING Jianguo. Local pseudo arc-length method for hyperbolic partial differential equation[J]. Chinese Journal of Computational Mechanics, 2014(3):384-389.

[6] 王瑞利,劉全,劉希強,等.人為解方法及其在流體力學(xué)程序驗證中的應(yīng)用[J].計算機應(yīng)用與軟件,2012(11):4-7.

WANG Ruili, LIU Quan, LIU Xiqiang, et al. Artificial solution and its application in verifying hydrodynamics program[J]. Computer Applications and Software, 2012(11):4-7.

[7] ROACHE P J. Verification and validation in computational science and engineering[M]. Albuquerque, NM, USA: Hermosa Publishers, 1998:8-9.

[8] OBERKAMPF W L, TRUCANO T G. Verification and validation benchmarks[J]. Nuclear Engineering and Design, 2008,238(3):716-743.

[9] ROY C J. Review of code and solution verification procedures for computational simulation[J]. Journal of Computational Physics, 2005,205(1):131-156.

[10] 鄧小剛,宗文剛,張來平,等.計算流體力學(xué)中的驗證與確認(rèn)[J].力學(xué)進展,2007,37(2):279-288.

DENG Xiaogang, ZONG Wengang, ZHANG Laiping, et al. Verification and validation in computational fluid dynamics[J]. Advances in Mechanics, 2007,37(2):279-288.

[11] 王瑞利,林忠,袁國興.科學(xué)計算程序的驗證和確認(rèn)[J].北京理工大學(xué)學(xué)報,2010,30(3):353-356.

WANG Ruili, LIN Zhong, YUAN Guoxing. Verification and validation in scientific computing code[J]. Transactions of Beijing Institute of Technology, 2010,30(3):353-356.

[12] NING J G, YUAN X P, MA T B, et al. Positivity-preserving moving mesh scheme for two-step reaction model in two dimensions[J]. Computers and Fluids, 2015,123:72-86.

[13] HAN J Q, TANG H Z. An adaptive moving mesh method for two-dimensional ideal magnetohydrodynamics[J]. Journal of Computational Physics, 2007,220(2):791-812.

[14] 曾現(xiàn)洋,劉睿,劉希強.人為解與人為解方法[J].聊城大學(xué)學(xué)報(自然科學(xué)版),2010,23(1):71-75.

ZENG Xianyang, LIU Rui, LIU Xiqiang. Manufactured solution and the method of manufactured solution[J]. Journal of Liaocheng University (Natural Science Edition), 2010,23(1):71-75.

[15] NING J G, WANG X, MA T B, et al. Numerical simulation of shock wave interaction with a deformable particle based on the pseudo arc-length method[J]. Science China Technological Sciences, 2015,58(5):848-857.

[16] GOLLAN R J, JACOBS P A. About the formulation, verification and validation of the hypersonic flow solver Eilmer[J]. International Journal for Numerical Methods in Fluids, 2013,73(1):19-57.

[17] LI J, REN H L, NING J G, et al. Additive Runge-Kutta methods for H2/O2/Ar detonation with a detailed elementary chemical reaction model[J]. Chinese Science Bulletin, 2013,58(11):1216-1227.