馬天寶，陳建良，寧建國(guó)，原新鵬

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

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

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

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

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

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

1 偽弧長(zhǎng)算法

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

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為邊界長(zhǎng)度微元，ni,j是邊界?Ki,j的單位外法向量，δ(w)=(F(w),G(w))。定義Lax-Friedrichs通量格式為：

(4)

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

(5)

對(duì)時(shí)間的離散，采用三階Runge-Kutta格式：

(6)

1.2 偽弧長(zhǎng)移動(dòng)算法

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

(7)

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

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

(8)

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

(9)

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

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

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

(10)

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

(11)

(12)

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

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

(13)

(14)

(15)

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

(16)

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

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

2 程序驗(yàn)證

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

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

通常，為了更好地使用人為解方法驗(yàn)證程序，人為解的選取需要遵循以下規(guī)定[14]：(1)人為解應(yīng)當(dāng)由光滑解析函數(shù)組成，例如多項(xiàng)式、指數(shù)函數(shù)、三角函數(shù)或者它們的組合，這樣方便求導(dǎo)運(yùn)算，對(duì)精度保證也極其重要；(2)人為解能夠使方程組的每一項(xiàng)都成立，即應(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ù)性。

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

(17)

狀態(tài)方程為：

(18)

(19)

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

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

(20)

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

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

收斂階計(jì)算公式為[15]：

(21)

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

(22)

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

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

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

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

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

4 結(jié) 論

將偽弧長(zhǎng)算法應(yīng)用在處理爆轟波的強(qiáng)間斷問(wèn)題，說(shuō)明將計(jì)算空間由物理空間轉(zhuǎn)換到弧長(zhǎng)空間的方法，詳細(xì)說(shuō)明了程序的人為解構(gòu)造過(guò)程，對(duì)比分析偽弧長(zhǎng)算法對(duì)計(jì)算精度的提高，驗(yàn)證偽弧長(zhǎng)算法求解非線性偏微分方程組的有效性；通過(guò)對(duì)比分析氣相爆轟波傳播過(guò)程中波陣面處網(wǎng)格的自適應(yīng)加密，說(shuō)明偽弧長(zhǎng)算法相比于直接有限體積方法在捕捉波陣面效果上的優(yōu)勢(shì)，計(jì)算過(guò)程中可以實(shí)現(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] 張德良,謝巍,郭長(zhǎng)銘,等.氣相爆轟胞格結(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ìn)展[J].科學(xué)通報(bào),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] 王星,馬天寶,寧建國(guó).雙曲偏微分方程的局部偽弧長(zhǎng)方法研究[J].計(jì)算力學(xué)學(xué)報(bào),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] 王瑞利,劉全,劉希強(qiáng),等.人為解方法及其在流體力學(xué)程序驗(yàn)證中的應(yīng)用[J].計(jì)算機(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] 鄧小剛,宗文剛,張來(lái)平,等.計(jì)算流體力學(xué)中的驗(yàn)證與確認(rèn)[J].力學(xué)進(jìn)展,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] 王瑞利,林忠,袁國(guó)興.科學(xué)計(jì)算程序的驗(yàn)證和確認(rèn)[J].北京理工大學(xué)學(xué)報(bào),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)洋,劉睿,劉希強(qiáng).人為解與人為解方法[J].聊城大學(xué)學(xué)報(bào)(自然科學(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.