基于誤差馬尾圖量化爆轟數(shù)值模擬結(jié)果的置信度＊

王瑞利，梁 霄，林 忠

基于誤差馬尾圖量化爆轟數(shù)值模擬結(jié)果的置信度＊

王瑞利1，梁 霄2，林 忠1

（1.北京應(yīng)用物理與計算數(shù)學(xué)研究所，北京100094；2.山東科技大學(xué)數(shù)學(xué)學(xué)院，山東 青島266590）

針對爆轟流體力學(xué)數(shù)值模擬過程中輸入?yún)?shù)的不確定性，通過抽樣技術(shù)，形成確定性爆轟流體力學(xué)程序的各種輸入和數(shù)值求解，建立輸入?yún)?shù)與輸出響應(yīng)量的樣本，再通過概率框架下的誤差累積分布函數(shù)與馬尾圖，給出了爆轟數(shù)值模擬過程中輸入?yún)?shù)不確定度對模擬結(jié)果影響的置信度量化方法。通過一維黎曼問題、平面爆轟問題計算了誤差馬尾圖，給出了二維爆轟拉氏自適應(yīng)流體動力學(xué)LAD2D程序計算網(wǎng)格與模擬結(jié)果置信度的關(guān)系，對多物理爆轟過程發(fā)展高置信度數(shù)值模擬軟件有很好的借鑒作用。

累積分布函數(shù)；誤差馬尾圖；置信度量化；爆轟數(shù)值模擬

LAD2D程序［1］是一個非結(jié)構(gòu)網(wǎng)格拉氏自適應(yīng)的二維爆轟流體動力學(xué)軟件，此程序已通過軟件質(zhì)量保證（SQA）和大量測試模型的考核，驗證了程序的正確性［2－4］，但由于爆轟流體力學(xué)物理過程的復(fù)雜性和人們認(rèn)識的缺陷，在建模過程中含有抽象、簡化和近似，逼真建模很難，有時只能唯象建模和逐漸逼近真實情況，建模過程含有不確定性。加之描述其過程的數(shù)學(xué)物理模型是高度非線性的偏微分方程組，很難解析求解［5－6］。在數(shù)值求解過程中，由于連續(xù)到離散，存在計算模型誤差、離散誤差、計算機舍入誤差和分析誤差等，數(shù)值模擬過程始終是一種近似，含有不確定性。為此，數(shù)值模擬結(jié)果的置信度一直缺乏科學(xué)的論述［7－8］，不確定度量化是置信度評估的核心，不確定度量化（uncertainty quantification，UQ）方法分為概率框架和非概率框架下的多種量化方法［9］。本文中采用概率框架下的誤差累積分布函數(shù)（cumulative distribution function，CDF）和互補累積分布函數(shù)（complementary cumulative distribution function，CCDF），進行統(tǒng)計分析，給出了CDF和CCDF曲線，形成爆轟計算中輸入不確定性參數(shù)、網(wǎng)格尺度與輸出響應(yīng)量的誤差馬尾狀圖，從而可量化分析爆轟模擬結(jié)果的置信度，對多物理爆轟過程發(fā)展高置信度數(shù)值模擬軟件有很好的借鑒作用。

1 數(shù)學(xué)模型

1.1 CDF方法及馬尾圖

假定在模型計算區(qū)間 ［a，b］上有一組標(biāo)準(zhǔn)解：X＝（x1，x2，…，xn），Y＝ （y1，y2，…，yn）。輸入?yún)?shù)為λ，通過數(shù)值計算得到一組數(shù)值解：CDF方法的步驟為：

（1）計算數(shù)值點上的標(biāo)準(zhǔn)解。在區(qū)間［a，b］上，由標(biāo)準(zhǔn)解插值得到點上的標(biāo)準(zhǔn)解

（3）統(tǒng)計分析。假設(shè)誤差分為K段，統(tǒng)計小于每段誤差的個數(shù)，然后計算出平均值，即概率：

由式（3）就可以得到馬尾圖，建立輸入?yún)?shù)λ對輸出響應(yīng)量誤差的關(guān)系圖，以量化輸入?yún)?shù)λ對輸出響應(yīng)量的置信度，判斷輸入?yún)?shù)計算條件的合理性。

1.2 爆轟流體力學(xué)模型

數(shù)值模擬使用的基本方程是不定?？蓧嚎s理想流體力學(xué)方程和化學(xué)動力學(xué)方程的藕合方程組：

呈靜止?fàn)顟B(tài)的凝固炸藥P＝0。為了進行數(shù)值計算，用一條光滑曲線將它們連接起來，通常引進所謂的燃燒函數(shù)F來表征炸藥反應(yīng)程度，這里考慮如下燃燒函數(shù)：Wilkins函數(shù)（時間燃燒函數(shù)＋C－J比容燃燒函數(shù)）：

式中：VJ＝γV0／γ＋1（）是C－J比容，V0＝1／ρ0，tb是爆轟波剛到達計算網(wǎng)格的時刻（開始燃燒），即起爆時間［10］，t是當(dāng)前計算時刻，ΔL＝rbΔR／DJ，ΔR 是網(wǎng)格寬度，DJ是 C－J爆轟速度，nb、rb是可調(diào)參數(shù)。計算時，考慮人為黏性。

（1）von Neumann－Richtmyer人為黏性（二次黏性）：

式中：lNR為具有長度量綱的量，l2NR＝a2NRA；aNR為N－R人為黏性系數(shù)；A為計算網(wǎng)格面積。

（2）Landshoff人為黏性（一次黏性或線性黏性）：

式中：lL為具有長度量綱的量，lL＝aL槡A，aL為Landshoff人為黏性系數(shù)；A 為計算網(wǎng)格面積；c為當(dāng)?shù)芈曀佟?/p>

1.3 CDF法與爆轟流體力學(xué)數(shù)值模擬結(jié)果置信度量化

1.2 節(jié)中爆轟模型計算涉及很多輸入?yún)?shù)，包括產(chǎn)物JWL狀態(tài)方程中A、B、R1、R2、w，反應(yīng)程度及燃燒函數(shù)中nb、rb和流體計算時人工黏性系數(shù)aNR、aL，這些輸入?yún)?shù)對輸出計算結(jié)果有多大影響，需要進行分析評估。

采用CDF法時，首先憑經(jīng)驗給出一套計算參數(shù)，針對其中一個參數(shù)進行擾動，其他參數(shù)固定不變，通過確定性程序，進行系列模擬計算，得到數(shù)值解。其次根據(jù)基準(zhǔn)解，可以是解析解，也可以是高精度數(shù)值解，計算數(shù)值解與基準(zhǔn)解的誤差。然后形成誤差的CDF及馬尾圖，以量化數(shù)值模擬結(jié)果的置信度。

2 計算結(jié)果與分析

2.1 一維黎曼問題

SOD問題常稱為激波管問題。許多間斷解方法的設(shè)計和構(gòu)造，都利用這種激波管問題進行可靠性和準(zhǔn)確度的數(shù)值實驗，從而判斷和檢驗方法、格式和程序的優(yōu)劣。SOD問題是一個非常柔和的例子，精確解包含了一個左稀疏波、接觸間斷和一個激波。下面采用SOD問題作為數(shù)值例子測試考核LAD2D程序的正確性和數(shù)值格式的精度。

SOD問題一般計算區(qū)域為［－1，1］，初值為：

式中：ρL、UL、PL分別指左端初始狀態(tài)的密度、速度與壓力，ρR、UR、PR分別指右端初始狀態(tài)的密度、速度與壓力。左右邊界采用連續(xù)邊界條件，狀態(tài)方程采用理想氣體：P＝γ－1（）ρe，γ＝1.4。計算采用CFL＝0.5，一次黏性系數(shù)取0.06，二次黏性系數(shù)取0.2。圖1給出了20～800個網(wǎng)格單元間隔10的16套網(wǎng)格計算到t＝0.5時計算結(jié)果與解析解的比較。圖2給出了16套網(wǎng)格計算到t＝0.5時計算結(jié)果與解析解的比較誤差的累積分布函數(shù)CDF的統(tǒng)計分析馬尾圖。從圖中可以清楚地看出計算誤差與網(wǎng)格尺度的關(guān)系。

從圖2可以清楚地查出LAD2D程序計算SOD問題的置信度情況。假設(shè)要求計算誤差在0.01（1%），那么LAD2D程序計算SOD問題只有在網(wǎng)格尺度為0.01（200個網(wǎng)格單元）下，置信度才能超過90%。假設(shè)要求計算誤差在0.02（2%），那么LAD2D程序計算SOD問題只有在網(wǎng)格尺度為0.025（80個網(wǎng)格單元）下，置信度就能超過90%。

爆炸波問題也是一個一維黎曼問題，由于有一邊狀態(tài)的壓力達到爆轟的壓力，所以稱為爆炸波問題。同SOD問題相比，此問題在激波處有一個很窄的強激波區(qū)，對程序格式健壯性考核有重要意義。爆炸波問題一般計算區(qū)域為［－1，1］，初值為：

式中：ρL、UL、PL、γL分別指左端狀態(tài)初始的密度、速度、壓力和理想氣體狀態(tài)方程系數(shù)，ρR、UR、PR、γR分別指右端狀態(tài)初始的密度、速度、壓力和理想氣體狀態(tài)方程系數(shù)。左右邊界采用連續(xù)邊界條件，狀態(tài)方程采用理想氣體：P＝γ－1（）ρe。計算采用CFL＝0.5，一次黏性系數(shù)取0.06，二次黏性系數(shù)取0.2。圖3給出了14套網(wǎng)格計算到t＝0.1時計算結(jié)果與解析解的比較誤差的累積分布函數(shù)CDF的統(tǒng)計分析馬尾圖。從圖中可以清楚地看出計算誤差與網(wǎng)格尺度的關(guān)系。假設(shè)要求計算誤差在0.01（1%），那么LAD2D程序計算爆炸波問題只有在網(wǎng)格尺度為0.003 3（600個網(wǎng)格單元）下，置信度才能超過90%。假設(shè)要求計算誤差在0.02（2%），那么LAD2D程序計算爆炸波問題只有在網(wǎng)格尺度為0.005（400個網(wǎng)格單元）下，置信度能超過90%。從SOD和爆炸波的統(tǒng)計分析看，網(wǎng)格計算尺度在0.04～0.025之間，計算誤差如果要求在0.01～0.05之間時，SOD問題的置信度為0.6～0.99，爆炸波問題的置信度為0.4～0.7，置信范圍太大，很難把握，需縮小置信范圍。

2.2 平面爆轟問題

眾所周知，爆轟波達到定態(tài)時，數(shù)值計算一定要符合Chapman－Jouguet理論。也就是說，在聲速點的波后流場與CJ模型相一致，即可達到CJ狀態(tài)。炸藥取PBX－9404炸藥，參數(shù)為：K＝2.827、DJ＝8.88km／s，ρ0＝1.842g／cm3，由CJ理論可以解析推出PBX－9404炸藥CJ狀態(tài)［11］為：pCJ＝37.27GPa。計算區(qū)域為［0，10］，初始在左邊起爆。此問題起爆后經(jīng)過一段時間可以達到穩(wěn)定的爆轟波及CJ狀態(tài)，以檢驗程序能否達到CJ狀態(tài)，是爆轟數(shù)值模擬軟件適應(yīng)性考核的最基本問題。

（1）網(wǎng)格尺度置信度分析

在爆轟計算條件（主要是輸入?yún)?shù)）固定的情況下，改變網(wǎng)格尺度，統(tǒng)計分析程序的置信度。起爆采用壓縮比σ＝1.03，燃燒函數(shù)采用 Wilkins函數(shù)，參數(shù)nb＝1.0，rb＝2.1。狀態(tài)方程采用JWL形式，輸入?yún)?shù)A＝765.788GPa、B＝14.249GPa、R1＝4.3、R2＝1.45、ω＝0.28。計算CFL＝0.7，一次黏性系數(shù)為0.06，二次黏性系數(shù)為2.0。圖4給出了1～10μs時10套（網(wǎng)格尺度分別為0.1、0.09、0.08、0.07、0.06、0.05、0.04、0.03、0.02、0.01）網(wǎng)格計算結(jié)果不同時刻與基準(zhǔn)解（選取10 000個網(wǎng)格的數(shù)值解作為基準(zhǔn)解）對比誤差的累積分布函數(shù)CDF。

針對一維平面爆轟，假設(shè)要求計算誤差在0.01（1%），從圖4中的t＝1μs時計算結(jié)果與網(wǎng)格尺度變化來看，網(wǎng)格尺度從0.1～0.01之間時，置信度都可超過90%。但從t＝4、8μs的計算結(jié)果可以看出，假設(shè)要求計算誤差在0.01（1%），網(wǎng)格尺度在0.1～0.01之間時，t＝4μs置信度為0.6～0.8，t＝8μs僅0.3～0.4。從這個分析可以看出，隨著計算時間的推進，置信度逐漸下降，甚至有的置信度很差。說明目前爆轟計算的建模很不適應(yīng)，需要大大改進計算模型。但無論從t＝1、4μs，還是從t＝8μs計算結(jié)果看，網(wǎng)格尺度達到0.01時，置信度大于0.6，所以在目前建模情況下，爆轟計算的網(wǎng)格尺度要盡量小，網(wǎng)格尺度最好小于0.01，甚至更小，才能大大提高爆轟計算的置信度。

（2）JWL狀態(tài)方程中輸入?yún)?shù)敏感性分析

爆轟計算中采用JWL狀態(tài)方程式（8）時，其參數(shù)A、B、R1、R2和w敏感性分析對爆轟模型的使用有重要意義，其參數(shù)取值大小對實際問題模擬結(jié)果的影響是爆轟模型及模擬結(jié)果用于決策的重要依據(jù)。針對爆轟產(chǎn)物JWL狀態(tài)方程參數(shù)隨機選取了4套（見表1）參數(shù)，其他輸入?yún)?shù)同（1）的計算條件，在網(wǎng)格尺度分別為0.1和0.01情況下進行了數(shù)值模擬，將網(wǎng)格尺度為0.001、JWL狀態(tài)方程參數(shù)取第一組時的計算結(jié)果作為基準(zhǔn)解。圖5為4套參數(shù)計算結(jié)果的誤差統(tǒng)計的累積分布函數(shù)CDF，其中基準(zhǔn)解參數(shù)采用PAP1，計算分點采用10 000個點。

表1 JWL狀態(tài)方程4套參數(shù)Table 1Four sets of parameters for JWL state eqution

假設(shè)基準(zhǔn)解選取精度很高時，從圖5可以看出，當(dāng)誤差要求為0.01（1%）時，4套參數(shù)對計算結(jié)果影響很大，在網(wǎng)格尺度為0.01時，密度分布置信度從0.9下降到0.3，壓力分布置信度從0.45下降到0.28。當(dāng)誤差要求為0.04（4%）時，4套參數(shù)對計算結(jié)果密度影響不大，置信度均大于0.9，壓力仍影響很大，置信度仍從0.45下降到0.28。從這個計算結(jié)果說明，選取爆轟產(chǎn)物JWL狀態(tài)方程參數(shù)應(yīng)該引起重視。但選取哪一套參數(shù)合理，置信度高，需要建立合理、精度高的基準(zhǔn)解，才能應(yīng)用本文方法合理選取參數(shù)。

3 結(jié) 論

（1）通過標(biāo)準(zhǔn)解與數(shù)值解之間的誤差累積分布函數(shù)（CDF）及馬尾圖表征方法，建立數(shù)值計算輸入?yún)?shù)與網(wǎng)格尺度對輸出響應(yīng)量影響的置信度的表征方法。

（2）基于自主開發(fā)的爆轟彈塑性流體力學(xué)軟件LAD2D，針對一維黎曼問題，給出了網(wǎng)格尺度變化下數(shù)值模擬結(jié)果與解析解之間誤差的馬尾表征圖。對于經(jīng)典的SOD問題，假設(shè)要求模擬誤差達到0.01（1%），那么計算網(wǎng)格尺度必須小于0.01，置信度才能超過90%。假設(shè)要求計算誤差達到0.02（2%），即降低要求時，計算網(wǎng)格尺度只要小于0.025，置信度就會超過90%。對于強激波的爆炸波問題，若要求計算誤差達到0.01（1%），那么計算網(wǎng)格尺度必須小于0.003 3，置信度才能超過90%。假設(shè)要求計算誤差達到0.02（2%），即降低要求時，計算網(wǎng)格尺度只要小于0.005，置信度就會超過90%。但從爆炸波與SOD問題相比，計算網(wǎng)格尺度與激波強弱關(guān)系很大，激波越強，需要計算網(wǎng)格尺度越小，這與理論分析是一致的。

（3）針對炸藥爆轟產(chǎn)物JWL狀態(tài)方程參數(shù)取值對計算結(jié)果的影響，給出了4套參數(shù)在兩種網(wǎng)格尺度下誤差的馬尾表征圖。當(dāng)誤差要達到0.01（1%）時，在網(wǎng)格尺度為0.01時，密度分布置信度從0.9下降到0.3，壓力分布置信度從0.45下降到0.28。當(dāng)誤差要求達到0.04（4%）時，即降低要求時，密度分布置信度均大于0.9，壓力分布置信度仍從0.45下降到0.28。從這個分析結(jié)果說明，選取爆轟產(chǎn)物JWL狀態(tài)方程參數(shù)應(yīng)該引起重視。

（4）研究結(jié)果表明，基于誤差馬尾圖量化參數(shù)敏感性和置信度的思想是可行的，為非線性、強間斷、多物理耦合問題的數(shù)值模擬結(jié)果置信度評估提供了一種行之有效的方法。

［1］ 王瑞利，林忠，溫萬治，等.多介質(zhì)拉氏自適應(yīng)流體動力學(xué)軟件LAD2D研制及其應(yīng)用［J］.計算機輔助工程，2014，23（2）：1－7.Wang Ruili，Lin Zhong，Wen Wanzhi，et al.Development and application of adaptive multi－media Lagrangian fluid dynamics software LAD2D［J］.Computer Aided Engineering，2014，23（2）：1－7.

［2］ 王瑞利，溫萬治.復(fù)雜工程建模與模擬的驗證與確認(rèn)［J］.計算機輔助工程，2014，23（4）：61－68.Wang Ruili，Wen Wanzhi.Advances in verification and validation of modeling and simulation of the complex engineering［J］.Computer Aided Engineering，2014，23（4）：61－68.

［3］ Wang Ruili，Liang Xiao，Lin Wenzhou，et al.Verification and validation of the detonation computational fluid dynamics model［J］.Defect and Diffusion Forum，2016，366（3）：40－46.

［4］ Wang Ruili，Zhang Shudao，Liu Quan.One method of constructing manufactured solutions to 2Dhydrodynamics Euler equations［C］∥Proceedings of the International Conference on Numerical Analysis and Applied Mathematics 2014（ICNAAM－2014），2015，1648（1）：599－620.

［5］ 張冠人，陳大年.凝聚炸藥起爆動力學(xué)［M］.北京：國防工業(yè)出版社，1991.

［6］ 孫錦山.凝聚炸藥非理想爆轟的數(shù)值模擬［J］.力學(xué)進展，1995，25（1）：127－133.Sun Jinshan.Numerical modeling of non－ideal detonation in condensed explosives［J］.Advances in Mechanics，1995，25（1）：127－133.

［7］ 王瑞利，江松.多物理耦合非線性偏微分方程與數(shù)值解不確定性量化數(shù)學(xué)方法［J］.中國科學(xué)（數(shù)學(xué)），2015，45（6）：723－738.Wang R L，Jiang S.Mathematical methods for uncertainty quantification in nonlinear multi－physics systems and their numerical simulations［J］.Scientia Sinica（Mathematica），2015，45（6）：723－738.

［8］ 王瑞利，劉全，溫萬治.非嵌入式多項式混沌法在爆轟產(chǎn)物JWL參數(shù)評估中的應(yīng)用［J］.爆炸與沖擊，2015，35（1）：9－15.Wang Ruili，Liu Quan，Wen Wanzhi.Non－intrusive polynomial chaos methods and its application in the parameters assessment of explosion product JWL［J］.Explosion and Shock Waves，2015，35（1）：9－15.

［9］ 湯濤，周濤.不確定性量化的高精度數(shù)值方法和理論［J］.中國科學(xué)（數(shù)學(xué)），2015，45（7）：891－928.Tang Tao，Zhou Tao.Recent developments in high order numerical methods for uncertainty quantification［J］.Scientia Sinica（Mathematica），2015，45（7）：891－928.

［10］ 王瑞利，林忠，魏蘭.利用前沿推進法計算復(fù)雜炸藥結(jié)構(gòu)的起爆時間［J］.高壓物理學(xué)報，2015，29（4）：286－292.Wang Ruili，Lin Zhong，Wei Lan.Calculating the initiation time of the explosive with complex structure using the advancing front technique［J］.Chinese Journal of High Pressure Physics，2015，29（4）：286－292.

［11］ 惲壽榕，趙衡陽.爆炸力學(xué)［M］.北京：國防工業(yè)出版社，2005.

Confidence level of numerical simulation of detonation through quantifying the horsetail of error

Wang Ruili1，Liang Xiao2，Lin Zhong1

（1.Institute of Applied Physics and Computational Mathematics，Beijing100094，China；2.College of Mathematics，Shandong University of Science and Technology，Qingdao 266590，Shandong，China）

In the present study，sampling technique was used to deal with the input parameter uncertainty in the numerical simulation of detonation CFD（computational fluid dynamics）.Then the deterministic detonation CFD program was constructed with different input.The sample of the input parameter and system response quantity was obtained through the previous result.The cumulative distribution function and the horsetail of error was utilized to achieve the confidence level，which was then used to assess the influence of the input parameter uncertainty on the simulation result of detonation CFD.The horsetail graph of error in one dimensional Riemann problem and planar detonation problem were presented to analyze the relationship between the confidence level of simulation result and the mesh used in LAD2D.This method provides a reference for developing the software of multiphysics detonation process on high confidence level.

cumulative distribution function；horsetail of error；quantification of confidence level；numerical simulation of detonation

O381；O241 國標(biāo)學(xué)科代碼：13035

A

10.11883／1001－1455（2017）06－0893－08

2016－05－03；

2016－07－02

國家自然科學(xué)基金項目（11372051，91630312，11475029）；國防基礎(chǔ)科學(xué)研究計劃項目（C1520110002）；

中國工程物理研究院科學(xué)技術(shù)發(fā)展基金項目（2015B0202045）

王瑞利（1964— ），男，研究員，wang＿ruili＠iapcm.ac.cn。

（責(zé)任編輯 曾月蓉）