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

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

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

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

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

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

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

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

1 偽弧長算法

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

1.1 數值離散方法

考慮如下雙曲守恒系統：

(1)

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

?Ki,jS(w)dσ

(2)

式中：dσ為面積微元。

利用Green-Gauss定理，可以將上式轉化為：

(3)

式中：|Ki,j|為網格Ki,j的面積；?Ki,j為網格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 偽弧長移動算法

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

(7)

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

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

(8)

定義偽弧長函數為：

(9)

通過調整弧長參數α1、α2的大小，可以控制網格的加密區域，通常ψ>0。

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

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

(10)

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

(11)

(12)

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

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

(13)

(14)

(15)

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

(16)

結合前文算法應用，將偽弧長算法應用在二維雙曲守恒系統，計算步驟總結如下。

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

2 程序驗證

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

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

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

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

(17)

狀態方程為：

(18)

(19)

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

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

(20)

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

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

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

(21)

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

(22)

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

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

3 氣相爆轟波傳播問題的數值模擬

模擬二維管道中氣相爆轟波傳播問題(所有物理量均為無量綱形式)，取計算域為[0,300]×[0,50]，測試偽弧長算法對爆轟波陣面的捕捉效果。固定網格的有限體積法初始網格數設定1 500×250，偽弧長算法網格數取為750×250。選取y=25.0的切面，研究網格節點沿x方向的分布情況。圖3為T=2.0時密度與壓力的計算結果。可以看出，兩種方法的波陣面基本重合，也就是說雖然偽弧長算法在x方向的網格數只有有限體積法的一半，但波陣面傳播位置保持一致，而且在波陣面處，偽弧長算法由于網格的自適應移動，依然可以高效地捕捉到間斷面。圖4為有限體積法和偽弧長算法在波陣面附近壓力云圖的局部放大圖。對比結果說明偽弧長算法網格分布會隨波陣面的傳播而移動，在波陣面附近網格分布和有限體積法網格分布效果相同，在未反應區和已反應區網格分布相對稀疏，這說明偽弧長算法用有限體積法一半數目的網格依然可以捕捉和追蹤爆轟波陣面的傳播，同時有效減小計算量。

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

4 結 論

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

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

[2] 張德良,謝巍,郭長銘,等.氣相爆轟胞格結構和馬赫反射數值模擬[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. 爆炸力學高精度數值模擬研究進展[J].科學通報,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] 王昌建,徐勝利.直管內胞格爆轟的基元反應數值研究[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].計算力學學報,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] 王瑞利,劉全,劉希強,等.人為解方法及其在流體力學程序驗證中的應用[J].計算機應用與軟件,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] 鄧小剛,宗文剛,張來平,等.計算流體力學中的驗證與確認[J].力學進展,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] 王瑞利,林忠,袁國興.科學計算程序的驗證和確認[J].北京理工大學學報,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] 曾現洋,劉睿,劉希強.人為解與人為解方法[J].聊城大學學報(自然科學版),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.