基于误差马尾图量化爆轰数值模拟结果的置信度＊

王瑞利，梁 霄，林 忠

基于误差马尾图量化爆轰数值模拟结果的置信度＊

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

（1.北京应用物理与计算数学研究所，北京100094；2.山东科技大学数学学院，山东 青岛266590）

针对爆轰流体力学数值模拟过程中输入参数的不确定性，通过抽样技术，形成确定性爆轰流体力学程序的各种输入和数值求解，建立输入参数与输出响应量的样本，再通过概率框架下的误差累积分布函数与马尾图，给出了爆轰数值模拟过程中输入参数不确定度对模拟结果影响的置信度量化方法。通过一维黎曼问题、平面爆轰问题计算了误差马尾图，给出了二维爆轰拉氏自适应流体动力学LAD2D程序计算网格与模拟结果置信度的关系，对多物理爆轰过程发展高置信度数值模拟软件有很好的借鉴作用。

累积分布函数；误差马尾图；置信度量化；爆轰数值模拟

LAD2D程序［1］是一个非结构网格拉氏自适应的二维爆轰流体动力学软件，此程序已通过软件质量保证（SQA）和大量测试模型的考核，验证了程序的正确性［2－4］，但由于爆轰流体力学物理过程的复杂性和人们认识的缺陷，在建模过程中含有抽象、简化和近似，逼真建模很难，有时只能唯象建模和逐渐逼近真实情况，建模过程含有不确定性。加之描述其过程的数学物理模型是高度非线性的偏微分方程组，很难解析求解［5－6］。在数值求解过程中，由于连续到离散，存在计算模型误差、离散误差、计算机舍入误差和分析误差等，数值模拟过程始终是一种近似，含有不确定性。为此，数值模拟结果的置信度一直缺乏科学的论述［7－8］，不确定度量化是置信度评估的核心，不确定度量化（uncertainty quantification，UQ）方法分为概率框架和非概率框架下的多种量化方法［9］。本文中采用概率框架下的误差累积分布函数（cumulative distribution function，CDF）和互补累积分布函数（complementary cumulative distribution function，CCDF），进行统计分析，给出了CDF和CCDF曲线，形成爆轰计算中输入不确定性参数、网格尺度与输出响应量的误差马尾状图，从而可量化分析爆轰模拟结果的置信度，对多物理爆轰过程发展高置信度数值模拟软件有很好的借鉴作用。

1 数学模型

1.1 CDF方法及马尾图

假定在模型计算区间 ［a，b］上有一组标准解：X＝（x1，x2，…，xn），Y＝ （y1，y2，…，yn）。输入参数为λ，通过数值计算得到一组数值解：CDF方法的步骤为：

（1）计算数值点上的标准解。在区间［a，b］上，由标准解插值得到点上的标准解

（3）统计分析。假设误差分为K段，统计小于每段误差的个数，然后计算出平均值，即概率：

由式（3）就可以得到马尾图，建立输入参数λ对输出响应量误差的关系图，以量化输入参数λ对输出响应量的置信度，判断输入参数计算条件的合理性。

1.2 爆轰流体力学模型

数值模拟使用的基本方程是不定常可压缩理想流体力学方程和化学动力学方程的藕合方程组：

呈静止状态的凝固炸药P＝0。为了进行数值计算，用一条光滑曲线将它们连接起来，通常引进所谓的燃烧函数F来表征炸药反应程度，这里考虑如下燃烧函数：Wilkins函数（时间燃烧函数＋C－J比容燃烧函数）：

式中：VJ＝γV0／γ＋1（）是C－J比容，V0＝1／ρ0，tb是爆轰波刚到达计算网格的时刻（开始燃烧），即起爆时间［10］，t是当前计算时刻，ΔL＝rbΔR／DJ，ΔR 是网格宽度，DJ是 C－J爆轰速度，nb、rb是可调参数。计算时，考虑人为黏性。

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

式中：lNR为具有长度量纲的量，l2NR＝a2NRA；aNR为N－R人为黏性系数；A为计算网格面积。

（2）Landshoff人为黏性（一次黏性或线性黏性）：

式中：lL为具有长度量纲的量，lL＝aL槡A，aL为Landshoff人为黏性系数；A 为计算网格面积；c为当地声速。

1.3 CDF法与爆轰流体力学数值模拟结果置信度量化

1.2 节中爆轰模型计算涉及很多输入参数，包括产物JWL状态方程中A、B、R1、R2、w，反应程度及燃烧函数中nb、rb和流体计算时人工黏性系数aNR、aL，这些输入参数对输出计算结果有多大影响，需要进行分析评估。

采用CDF法时，首先凭经验给出一套计算参数，针对其中一个参数进行扰动，其他参数固定不变，通过确定性程序，进行系列模拟计算，得到数值解。其次根据基准解，可以是解析解，也可以是高精度数值解，计算数值解与基准解的误差。然后形成误差的CDF及马尾图，以量化数值模拟结果的置信度。

2 计算结果与分析

2.1 一维黎曼问题

SOD问题常称为激波管问题。许多间断解方法的设计和构造，都利用这种激波管问题进行可靠性和准确度的数值实验，从而判断和检验方法、格式和程序的优劣。SOD问题是一个非常柔和的例子，精确解包含了一个左稀疏波、接触间断和一个激波。下面采用SOD问题作为数值例子测试考核LAD2D程序的正确性和数值格式的精度。

SOD问题一般计算区域为［－1，1］，初值为：

式中：ρL、UL、PL分别指左端初始状态的密度、速度与压力，ρR、UR、PR分别指右端初始状态的密度、速度与压力。左右边界采用连续边界条件，状态方程采用理想气体：P＝γ－1（）ρe，γ＝1.4。计算采用CFL＝0.5，一次黏性系数取0.06，二次黏性系数取0.2。图1给出了20～800个网格单元间隔10的16套网格计算到t＝0.5时计算结果与解析解的比较。图2给出了16套网格计算到t＝0.5时计算结果与解析解的比较误差的累积分布函数CDF的统计分析马尾图。从图中可以清楚地看出计算误差与网格尺度的关系。

从图2可以清楚地查出LAD2D程序计算SOD问题的置信度情况。假设要求计算误差在0.01（1%），那么LAD2D程序计算SOD问题只有在网格尺度为0.01（200个网格单元）下，置信度才能超过90%。假设要求计算误差在0.02（2%），那么LAD2D程序计算SOD问题只有在网格尺度为0.025（80个网格单元）下，置信度就能超过90%。

爆炸波问题也是一个一维黎曼问题，由于有一边状态的压力达到爆轰的压力，所以称为爆炸波问题。同SOD问题相比，此问题在激波处有一个很窄的强激波区，对程序格式健壮性考核有重要意义。爆炸波问题一般计算区域为［－1，1］，初值为：

式中：ρL、UL、PL、γL分别指左端状态初始的密度、速度、压力和理想气体状态方程系数，ρR、UR、PR、γR分别指右端状态初始的密度、速度、压力和理想气体状态方程系数。左右边界采用连续边界条件，状态方程采用理想气体：P＝γ－1（）ρe。计算采用CFL＝0.5，一次黏性系数取0.06，二次黏性系数取0.2。图3给出了14套网格计算到t＝0.1时计算结果与解析解的比较误差的累积分布函数CDF的统计分析马尾图。从图中可以清楚地看出计算误差与网格尺度的关系。假设要求计算误差在0.01（1%），那么LAD2D程序计算爆炸波问题只有在网格尺度为0.003 3（600个网格单元）下，置信度才能超过90%。假设要求计算误差在0.02（2%），那么LAD2D程序计算爆炸波问题只有在网格尺度为0.005（400个网格单元）下，置信度能超过90%。从SOD和爆炸波的统计分析看，网格计算尺度在0.04～0.025之间，计算误差如果要求在0.01～0.05之间时，SOD问题的置信度为0.6～0.99，爆炸波问题的置信度为0.4～0.7，置信范围太大，很难把握，需缩小置信范围。

2.2 平面爆轰问题

众所周知，爆轰波达到定态时，数值计算一定要符合Chapman－Jouguet理论。也就是说，在声速点的波后流场与CJ模型相一致，即可达到CJ状态。炸药取PBX－9404炸药，参数为：K＝2.827、DJ＝8.88km／s，ρ0＝1.842g／cm3，由CJ理论可以解析推出PBX－9404炸药CJ状态［11］为：pCJ＝37.27GPa。计算区域为［0，10］，初始在左边起爆。此问题起爆后经过一段时间可以达到稳定的爆轰波及CJ状态，以检验程序能否达到CJ状态，是爆轰数值模拟软件适应性考核的最基本问题。

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

在爆轰计算条件（主要是输入参数）固定的情况下，改变网格尺度，统计分析程序的置信度。起爆采用压缩比σ＝1.03，燃烧函数采用 Wilkins函数，参数nb＝1.0，rb＝2.1。状态方程采用JWL形式，输入参数A＝765.788GPa、B＝14.249GPa、R1＝4.3、R2＝1.45、ω＝0.28。计算CFL＝0.7，一次黏性系数为0.06，二次黏性系数为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个网格的数值解作为基准解）对比误差的累积分布函数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状态方程中输入参数敏感性分析

爆轰计算中采用JWL状态方程式（8）时，其参数A、B、R1、R2和w敏感性分析对爆轰模型的使用有重要意义，其参数取值大小对实际问题模拟结果的影响是爆轰模型及模拟结果用于决策的重要依据。针对爆轰产物JWL状态方程参数随机选取了4套（见表1）参数，其他输入参数同（1）的计算条件，在网格尺度分别为0.1和0.01情况下进行了数值模拟，将网格尺度为0.001、JWL状态方程参数取第一组时的计算结果作为基准解。图5为4套参数计算结果的误差统计的累积分布函数CDF，其中基准解参数采用PAP1，计算分点采用10 000个点。

表1 JWL状态方程4套参数Table 1Four sets of parameters for JWL state eqution

假设基准解选取精度很高时，从图5可以看出，当误差要求为0.01（1%）时，4套参数对计算结果影响很大，在网格尺度为0.01时，密度分布置信度从0.9下降到0.3，压力分布置信度从0.45下降到0.28。当误差要求为0.04（4%）时，4套参数对计算结果密度影响不大，置信度均大于0.9，压力仍影响很大，置信度仍从0.45下降到0.28。从这个计算结果说明，选取爆轰产物JWL状态方程参数应该引起重视。但选取哪一套参数合理，置信度高，需要建立合理、精度高的基准解，才能应用本文方法合理选取参数。

3 结 论

（1）通过标准解与数值解之间的误差累积分布函数（CDF）及马尾图表征方法，建立数值计算输入参数与网格尺度对输出响应量影响的置信度的表征方法。

（2）基于自主开发的爆轰弹塑性流体力学软件LAD2D，针对一维黎曼问题，给出了网格尺度变化下数值模拟结果与解析解之间误差的马尾表征图。对于经典的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状态方程参数取值对计算结果的影响，给出了4套参数在两种网格尺度下误差的马尾表征图。当误差要达到0.01（1%）时，在网格尺度为0.01时，密度分布置信度从0.9下降到0.3，压力分布置信度从0.45下降到0.28。当误差要求达到0.04（4%）时，即降低要求时，密度分布置信度均大于0.9，压力分布置信度仍从0.45下降到0.28。从这个分析结果说明，选取爆轰产物JWL状态方程参数应该引起重视。

（4）研究结果表明，基于误差马尾图量化参数敏感性和置信度的思想是可行的，为非线性、强间断、多物理耦合问题的数值模拟结果置信度评估提供了一种行之有效的方法。

［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］.北京：国防工业出版社，1991.

［6］ 孙锦山.凝聚炸药非理想爆轰的数值模拟［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］ 王瑞利，江松.多物理耦合非线性偏微分方程与数值解不确定性量化数学方法［J］.中国科学（数学），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参数评估中的应用［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］ 汤涛，周涛.不确定性量化的高精度数值方法和理论［J］.中国科学（数学），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］.北京：国防工业出版社，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）；

中国工程物理研究院科学技术发展基金项目（2015B0202045）

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

（责任编辑 曾月蓉）