• ISSN 1001-1455  CN 51-1148/O3
  • EI Compendex、CA收录
  • 力学类中文核心期刊
  • 中国科技核心期刊、CSCD统计源期刊
高级检索 E-mail Alert

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

王瑞利 梁霄 林忠

引用本文:
Citation:

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

    作者简介: 王瑞利(1964—), 男, 研究员, wang_ruili@iapcm.ac.cn;
  • 基金项目: 国家自然科学基金项目 91630312
    国防基础科学研究计划项目 C1520110002
    中国工程物理研究院科学技术发展基金项目 2015B0202045
    国家自然科学基金项目 11372051
    国家自然科学基金项目 11475029

  • 中图分类号: O381;O241

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

  • CLC number: O381;O241

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

    Figure 1.  Comparison between numerical results and analytic solutionsfor different mesh scales of 0.1~0.025(16 sets mesh)

    图 2  网格尺度0.1~0.025(16套)计算结果与解析解的对比误差统计分析的马尾图

    Figure 2.  Horsetail figure of error between numerical results and analytic solutions for different mesh scales of 0.1~0.025(16 sets mesh)

    图 3  网格尺度0.04~0.025(14套)时计算结果与解析解的对比误差统计分析的马尾图

    Figure 3.  Horsetail figure of error between numerical results and analytic solutions for different mesh scales of 0.1~0.025(14 sets mesh)

    图 4  网格尺度0.1~0.01(10套)时不同时刻计算结果误差统计分析的马尾图

    Figure 4.  Horsetail figure of error between numerical results and analytic solutions at different times

    图 5  4套参数在不同网格尺度下的计算结果误差统计分析的马尾图(t=8 μs)

    Figure 5.  Horsetail figure of error with four sets parameters for different mesh scales(t=8 μs)

    表 1  JWL状态方程4套参数

    Table 1.  Four sets of parameters for JWL state eqution

    参数代号 A/GPa B/GPa R1 R2 ω
    PAP1 765.788 14.249 4.3 1.45 0.28
    PAP2 947.774 59.206 5.0 2.00 0.28
    PAP2 1 502.851 134.913 6.3 2.40 0.28
    PAP3 3 250.433 233.271 8.2 2.80 0.28
    下载: 导出CSV
  • [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 2D hydrodynamics 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.
  • [1] 章冠人 . 爆轰的数值模拟. 爆炸与冲击, 1981, 1(1): 99-99.
    [2] 肖作智 . 球面爆轰波的数值模拟. 爆炸与冲击, 1985, 5(2): 51-58.
    [3] 梁霄王瑞利 . 基于自适应和投影Wiener混沌的圆筒实验不确定度量化. 爆炸与冲击, 2019, 39(4): 041408-1-041408-12. doi: 10.11883/bzycj-2018-0253
    [4] 孙锦山曹菊珍张铁桥龚静号 . 钝感炸药爆轰性能的数值模拟. 爆炸与冲击, 1995, 15(1): 28-36.
    [5] 代淑兰许厚谦 . 驻定斜爆轰波并行数值模拟. 爆炸与冲击, 2006, 26(6): 572-576. doi: 10.11883/1001-1455(2006)06-0572-05
    [6] 袁帅柏劲松张汉钊李平 . 滑移爆轰密封的初步数值模拟. 爆炸与冲击, 2015, 35(4): 496-500. doi: 10.11883/1001-1455(2015)04-0496-05
    [7] 焦清介魏继锋周钢严楠蔡瑞娇 . 爆轰波拐角传播三维数值模拟. 爆炸与冲击, 2003, 23(6): 534-538.
    [8] 袁帅文尚刚李平董玉斌 . 强爆轰驱动飞片的数值模拟研究. 爆炸与冲击, 2015, 35(2): 197-202. doi: 10.11883/1001-1455(2015)02-0197-06
    [9] 金柯李平吴强金孝刚 . 爆轰产物驱动飞片运动数值模拟研究. 爆炸与冲击, 2004, 24(5): 419-424.
    [10] 朱号锋靳平肖卫国 . 硬岩中地下爆炸震源函数的数值模拟. 爆炸与冲击, 2009, 29(6): 648-653. doi: 10.11883/1001-1455(2009)06-0648-06
    [11] 洪滔秦承森 . 爆轰波管中铝粉尘爆轰的数值模拟. 爆炸与冲击, 2004, 24(3): 193-200.
    [12] 龙仁荣付跃升张庆明 . 水下爆炸爆源定位方法与误差分析. 爆炸与冲击, 2013, 33(2): 181-185. doi: 10.11883/1001-1455(2013)02-0181-05
    [13] 梁霄王瑞利 . 混合不确定度量化方法及其在计算流体动力学迎风格式中的应用. 爆炸与冲击, 2016, 36(4): 509-515. doi: 10.11883/1001-1455(2016)04-0509-07
    [14] 刘云峰王健平 . 有限谱ENO格式在爆轰波数值模拟中的应用. 爆炸与冲击, 2003, 23(4): 343-348.
    [15] 刘尔岩刘邦第王元书 . 应用阵面追踪法对散心爆轰波传播的数值模拟. 爆炸与冲击, 2001, 21(1): 57-61.
    [16] 刘晓利李鸿志 . 玉米粉-氧气混合物中爆轰波的数值模拟. 爆炸与冲击, 1994, 14(3): 208-216.
    [17] 姜羲王荪源 . 含能材料密实床燃烧转爆轰的数值模拟. 爆炸与冲击, 1992, 12(2): 97-105.
    [18] 朱建士魏振典周德忠 . 定常爆轰数值模拟中人为粘性与人为反应率的选取. 爆炸与冲击, 1983, 3(1): 21-27.
    [19] 孙晓晖陈志华张焕好 . 激波绕射碰撞加速诱导爆轰的数值模拟. 爆炸与冲击, 2011, 31(4): 407-412. doi: 10.11883/1001-1455(2011)04-0407-06
    [20] 潘振华范宝春归明月张旭东 . 流动系统中爆轰波传播特性的数值模拟. 爆炸与冲击, 2010, 30(6): 593-597. doi: 10.11883/1001-1455(2010)06-0593-05
  • 加载中
图(5)表(1)
计量
  • 文章访问数:  504
  • HTML全文浏览量:  41
  • PDF下载量:  407
  • 被引次数: 0
出版历程
  • 收稿日期:  2016-05-03
  • 录用日期:  2016-07-02
  • 刊出日期:  2017-11-25

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

    作者简介:王瑞利(1964—), 男, 研究员, wang_ruili@iapcm.ac.cn
  • 1. 北京应用物理与计算数学研究所, 北京 100094
  • 2. 山东科技大学数学学院, 山东 青岛 266590
基金项目:  国家自然科学基金项目 91630312国防基础科学研究计划项目 C1520110002中国工程物理研究院科学技术发展基金项目 2015B0202045国家自然科学基金项目 11372051国家自然科学基金项目 11475029

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

English Abstract

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

    • 假定在模型计算区间[a, b]上有一组标准解:X=x1, x2, …, xn, Y=y1, y2, …, yn。输入参数为λ, 通过数值计算得到一组数值解:$\tilde X = \left( {{{\tilde x}_1}, {{\tilde x}_2}, \cdots , {{\tilde x}_\mathit{m}}} \right), \mathit{\tilde Y = }\left( {{{\tilde y}_1}, {{\tilde y}_2}, \cdots , {{\tilde y}_\mathit{m}}} \right)$。CDF方法的步骤为:

      (1) 计算数值点上的标准解。在区间a, b上, 由标准解插值得到$\tilde X = \left( {{{\tilde x}_1}, {{\tilde x}_2}, \cdots , {{\tilde x}_\mathit{m}}} \right)$点上的标准解$Y' = \left( {{{y'}_1}, {{y'}_2}, \cdots , {{y'}_m}} \right)$:

      $ {{y'}_i} = \frac{{{x_j}_{ + 1} - {{\tilde x}_i}}}{{{x_j}_{ + 1} - {{\tilde x}_j}}}{y_j}\frac{{{{\tilde x}_i} - {x_j}}}{{{x_j}_{ + 1} - {x_j}}}{y_{j + 1}}\;\;\;\;\;\;\;\;\;{{\tilde x}_i} \in \left[ {{x_j}, {x_j}_{ + 1}} \right] $

      这样通过式(1), 就可得到$\tilde X = \left( {{{\tilde x}_1}, {{\tilde x}_2}, \cdots , {{\tilde x}_\mathit{m}}} \right)$对应的标准解$Y' = \left( {{{y'}_1}, {{y'}_2}, \cdots , {{y'}_m}} \right)$。

      (2) 计算误差。由${\tilde X}$的标准解$Y' = \left( {{{y'}_1}, {{y'}_2}, \cdots , {{y'}_m}} \right)$和数值解$\mathit{\tilde Y = }\left( {{{\tilde y}_1}, {{\tilde y}_2}, \cdots , {{\tilde y}_\mathit{m}}} \right)$, 就可以得到计算误差。即

      $ \Delta {y_i} = {{\tilde y}_i} - {{y'}_i}\;\;\;\;\;\;\;\;\;\;i = 1, 2, \cdots , m $

      (3) 统计分析。假设误差分为K段, 统计小于每段误差的个数, 然后计算出平均值, 即概率:

      $ F\left( {{\rm{erro}}{{\rm{r}}_l}} \right) = p\left( {\Delta {y_i} \le {\rm{erro}}{{\rm{r}}_l}} \right)\;\;\;\;\;\;\;\;\;l = 1, 2, \cdots , K $

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

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

      $ 质量方程\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\frac{{\partial \rho }}{{\partial t}} + \nabla \cdot \left( {\rho \;\mathit{\boldsymbol{u}}} \right) = 0 $

      $ 动量方程\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\frac{{\partial \left( {\rho \mathit{\boldsymbol{u}}} \right)}}{{\partial t}} + \nabla \cdot \left( {\rho \;\mathit{\boldsymbol{uu}}} \right) + \nabla P = 0 $

      $ 能量方程\;\;\;\;\;\;\;\;\;\;\;\;\;\frac{{\partial \left( {\rho \mathit{E}} \right)}}{{\partial t}} + \nabla \cdot \left( {\rho E\;\mathit{\boldsymbol{u}}} \right) + \nabla \cdot \left( {P\;\mathit{\boldsymbol{u}}} \right) = 0 $

      $ 状态方程\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;P = P\left( {\rho , e, F} \right) $

      式中:E=e+12u2, 其中ρuEeP分别表示密度、速度、单位质量的总能量、单位质量的内能与压力。方程组(4)~(7)的主要问题是状态方程(7)与一般的流体力学计算不同, 其中包含了燃烧函数F, 燃烧函数F要能正确反映化学反应的特性。爆炸产物的状态方程可以是JWL形式的状态方程, 也可以是理想气体状态方程P=(γ-1)ρe, 这里选用JWL形式状态方程:

      $ P = A\left( {1 - \frac{w}{{{R_1}V}}} \right){{\rm{e}}^{ - {R_1}V}} + B\left( {1 - \frac{w}{{{R_2}V}}} \right){{\rm{e}}^{ - {R_2}V}} + \frac{{wE}}{V} $

      式中:$V = \frac{v}{{{v_0}}}$, E为比热容力学能, ABR1R2w是常数, vv0是比容。

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

      $ F = {\left[ {{\rm{max}}\left( {{F_1}, {F_2}} \right)} \right]^{{n_{\rm{b}}}}} $

      $ {F_1} = \left\{ \begin{array}{l} 0\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;V \ge V\\ \frac{{{V_0} - V}}{{{V_0} - {V_{\rm{J}}}}}\;\;\;\;\;\;\;\;\;{V_{\rm{J}}} < V < {V_0}\;\\ 1\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;V \le \;{V_{\rm{J}}} \end{array} \right. $

      $ {F_2} = \left\{ \begin{array}{l} 0\;\;\;\;\;\;\;\;\;\;\;\;\;t \le {t_{\rm{b}}}\\ \frac{{t - {t_{\rm{b}}}}}{{\Delta L}}\;\;\;\;\;\;\;\;{t_{\rm{b}}} < t\; < \;{t_{\rm{b}}} + \Delta L\\ 1\;\;\;\;\;\;\;\;\;\;\;\;\;t \ge {t_{\rm{b}}} + \Delta L \end{array} \right. $

      式中:VJ=γV0/(γ+1)是C-J比容, V0=1/ρ0, tb是爆轰波刚到达计算网格的时刻(开始燃烧), 即起爆时间[10], t是当前计算时刻, ΔL=rbΔR/DJ, ΔR是网格宽度, DJ是C-J爆轰速度, nbrb是可调参数。计算时, 考虑人为黏性。

      (1) von Neumann-Richtmyer人为黏性(二次黏性):

      $ {q_{{\rm{NR}}}} = \left\{ \begin{array}{l} l_{{\rm{NR}}\rho }^2{\left( {\frac{{\mathop V\limits^. }}{V}} \right)^2}\;\;\;\;\;\;\;\mathop V\limits^. < 0\\ 0;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\mathop V\limits^. \ge 0 \end{array} \right. $

      式中:lNR为具有长度量纲的量, $l_{{\rm{NR}}}^2 = a_{{\rm{NR}}}^2A$;aNR为N-R人为黏性系数;A为计算网格面积。

      (2) Landshoff人为黏性(一次黏性或线性黏性):

      $ {q_{\rm{L}}} = \left\{ \begin{array}{l} {l_{\rm{L}}}\rho \;c\left( {\frac{{\mathop V\limits^. }}{V}} \right)\;\;\;\;\;\;\;\;\;\;\mathop V\limits^. < 0\\ 0;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\mathop V\limits^. \ge 0 \end{array} \right. $

      式中:lL为具有长度量纲的量, ${l_{\rm{L}}} = {a_\rm{L}}\sqrt A $, aL为Landshoff人为黏性系数;A为计算网格面积;c为当地声速。

    • 1.2节中爆轰模型计算涉及很多输入参数, 包括产物JWL状态方程中ABR1R2w, 反应程度及燃烧函数中nbrb和流体计算时人工黏性系数aNRaL这些输入参数对输出计算结果有多大影响, 需要进行分析评估。

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

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

      SOD问题一般计算区域为[-1, 1], 初值为:

      $ \left\{ \begin{array}{l} {\rho _{\rm{L}}} = 1.0\\ {U_{\rm{L}}} = 0.0\\ {P_{\rm{L}}} = 1.0 \end{array} \right.\;\;\;\;\;\;\;\;\;\;\;x \le 0, \;\;\;\;\;\;\;\;\;\;\;\left\{ \begin{array}{l} {\rho _{\rm{R}}} = 0.125\\ {U_{\rm{R}}} = 0.0\\ {P_{\rm{R}}} = 0.1 \end{array} \right.\;\;\;\;\;\;\;\;\;\;\;x > 0 $

      式中:ρLULPL分别指左端初始状态的密度、速度与压力, ρRURPR分别指右端初始状态的密度、速度与压力。左右边界采用连续边界条件, 状态方程采用理想气体: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的统计分析马尾图。从图中可以清楚地看出计算误差与网格尺度的关系。

      图  1  网格尺度0.1~0.025(16套)时计算结果与解析解的对比

      Figure 1.  Comparison between numerical results and analytic solutionsfor different mesh scales of 0.1~0.025(16 sets mesh)

      图  2  网格尺度0.1~0.025(16套)计算结果与解析解的对比误差统计分析的马尾图

      Figure 2.  Horsetail figure of error between numerical results and analytic solutions for different mesh scales of 0.1~0.025(16 sets mesh)

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

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

      $ \left\{ \begin{array}{l} {\rho _{\rm{L}}} = 1.4\\ {U_{\rm{L}}} = - 2.6\\ {P_{\rm{L}}} = 37.176\;5\\ {\gamma _{\rm{L}}} = 1.4 \end{array} \right.\;\;\;\;\;\;\;\;x \le 0, \;\;\;\;\;\;\;\;\;\;\left\{ \begin{array}{l} {\rho _{\rm{R}}} = 15.293\\ {U_{\rm{R}}} = 0.0\\ {P_{\rm{R}}} = 1.0\\ {\gamma _{\rm{R}}} = 3.1 \end{array} \right.\;\;\;\;\;\;\;\;\;x > 0 $

      式中:ρLULPLγL分别指左端状态初始的密度、速度、压力和理想气体状态方程系数, ρRURPRγ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, 置信范围太大, 很难把握, 需缩小置信范围。

      图  3  网格尺度0.04~0.025(14套)时计算结果与解析解的对比误差统计分析的马尾图

      Figure 3.  Horsetail figure of error between numerical results and analytic solutions for different mesh scales of 0.1~0.025(14 sets mesh)

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

      (1) 网格尺度置信度分析

      在爆轰计算条件(主要是输入参数)固定的情况下, 改变网格尺度, 统计分析程序的置信度。起爆采用压缩比σ=1.03, 燃烧函数采用Wilkins函数, 参数nb=1.0, rb=2.1。状态方程采用JWL形式, 输入参数A=765.788 GPa、B=14.249 GPa、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。

      图  4  网格尺度0.1~0.01(10套)时不同时刻计算结果误差统计分析的马尾图

      Figure 4.  Horsetail figure of error between numerical results and analytic solutions at different times

      针对一维平面爆轰, 假设要求计算误差在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)时, 其参数ABR1R2w敏感性分析对爆轰模型的使用有重要意义, 其参数取值大小对实际问题模拟结果的影响是爆轰模型及模拟结果用于决策的重要依据。针对爆轰产物JWL状态方程参数随机选取了4套(见表 1)参数, 其他输入参数同(1)的计算条件, 在网格尺度分别为0.1和0.01情况下进行了数值模拟, 将网格尺度为0.001、JWL状态方程参数取第一组时的计算结果作为基准解。图 5为4套参数计算结果的误差统计的累积分布函数CDF, 其中基准解参数采用PAP1, 计算分点采用10 000个点。

      参数代号 A/GPa B/GPa R1 R2 ω
      PAP1 765.788 14.249 4.3 1.45 0.28
      PAP2 947.774 59.206 5.0 2.00 0.28
      PAP2 1 502.851 134.913 6.3 2.40 0.28
      PAP3 3 250.433 233.271 8.2 2.80 0.28

      表 1  JWL状态方程4套参数

      Table 1.  Four sets of parameters for JWL state eqution

      图  5  4套参数在不同网格尺度下的计算结果误差统计分析的马尾图(t=8 μs)

      Figure 5.  Horsetail figure of error with four sets parameters for different mesh scales(t=8 μs)

      假设基准解选取精度很高时, 从图 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状态方程参数应该引起重视。但选取哪一套参数合理, 置信度高, 需要建立合理、精度高的基准解, 才能应用本文方法合理选取参数。

    • (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) 研究结果表明, 基于误差马尾图量化参数敏感性和置信度的思想是可行的, 为非线性、强间断、多物理耦合问题的数值模拟结果置信度评估提供了一种行之有效的方法。

参考文献 (11)

目录

    /

    返回文章
    返回