王瑞利，梁 霄，林 忠

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

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

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

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

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

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

1 數(shù)學模型

1.1 CDF方法及馬尾圖

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

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

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

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

1.2 爆轟流體力學模型

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

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

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

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

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

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

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

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

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

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

2 計算結果與分析

2.1 一維黎曼問題

SOD問題常稱為激波管問題。許多間斷解方法的設計和構造，都利用這種激波管問題進行可靠性和準確度的數(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個網格單元間隔10的16套網格計算到t＝0.5時計算結果與解析解的比較。圖2給出了16套網格計算到t＝0.5時計算結果與解析解的比較誤差的累積分布函數(shù)CDF的統(tǒng)計分析馬尾圖。從圖中可以清楚地看出計算誤差與網格尺度的關系。

從圖2可以清楚地查出LAD2D程序計算SOD問題的置信度情況。假設要求計算誤差在0.01（1%），那么LAD2D程序計算SOD問題只有在網格尺度為0.01（200個網格單元）下，置信度才能超過90%。假設要求計算誤差在0.02（2%），那么LAD2D程序計算SOD問題只有在網格尺度為0.025（80個網格單元）下，置信度就能超過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套網格計算到t＝0.1時計算結果與解析解的比較誤差的累積分布函數(shù)CDF的統(tǒng)計分析馬尾圖。從圖中可以清楚地看出計算誤差與網格尺度的關系。假設要求計算誤差在0.01（1%），那么LAD2D程序計算爆炸波問題只有在網格尺度為0.003 3（600個網格單元）下，置信度才能超過90%。假設要求計算誤差在0.02（2%），那么LAD2D程序計算爆炸波問題只有在網格尺度為0.005（400個網格單元）下，置信度能超過90%。從SOD和爆炸波的統(tǒ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］，初始在左邊起爆。此問題起爆后經過一段時間可以達到穩(wěn)定的爆轟波及CJ狀態(tài)，以檢驗程序能否達到CJ狀態(tài)，是爆轟數(shù)值模擬軟件適應性考核的最基本問題。

（1）網格尺度置信度分析

在爆轟計算條件（主要是輸入參數(shù)）固定的情況下，改變網格尺度，統(tǒng)計分析程序的置信度。起爆采用壓縮比σ＝1.03，燃燒函數(shù)采用 Wilkins函數(shù)，參數(shù)nb＝1.0，rb＝2.1。狀態(tài)方程采用JWL形式，輸入參數(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套（網格尺度分別為0.1、0.09、0.08、0.07、0.06、0.05、0.04、0.03、0.02、0.01）網格計算結果不同時刻與基準解（選取10 000個網格的數(shù)值解作為基準解）對比誤差的累積分布函數(shù)CDF。

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

（2）JWL狀態(tài)方程中輸入參數(shù)敏感性分析

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

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

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

3 結 論

（1）通過標準解與數(shù)值解之間的誤差累積分布函數(shù)（CDF）及馬尾圖表征方法，建立數(shù)值計算輸入參數(shù)與網格尺度對輸出響應量影響的置信度的表征方法。

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

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

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

［1］ 王瑞利，林忠，溫萬治，等.多介質拉氏自適應流體動力學軟件LAD2D研制及其應用［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］ 王瑞利，溫萬治.復雜工程建模與模擬的驗證與確認［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］ 張冠人，陳大年.凝聚炸藥起爆動力學［M］.北京：國防工業(yè)出版社，1991.

［6］ 孫錦山.凝聚炸藥非理想爆轟的數(shù)值模擬［J］.力學進展，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ù)學方法［J］.中國科學（數(shù)學），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］ 王瑞利，劉全，溫萬治.非嵌入式多項式混沌法在爆轟產物JWL參數(shù)評估中的應用［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］.中國科學（數(shù)學），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］ 王瑞利，林忠，魏蘭.利用前沿推進法計算復雜炸藥結構的起爆時間［J］.高壓物理學報，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］ 惲壽榕，趙衡陽.爆炸力學［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 國標學科代碼：13035

A

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

2016－05－03；

2016－07－02

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

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

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

（責任編輯 曾月蓉）