具体实施方式

需要说明的是，本发明中：

所有物理量的单位均采用本领域常用单位即可。

CJ指的是爆轰波的CJ理论，即Chapman和Jouget在20世纪初分别提出的关于爆轰波的平面一维流体动力学理论。

针对现有技术存在的不足，本发明公开了一种混合炸药爆轰参数计算方法，以多相多组分热力学基本原理为基础，结合混合炸药配方组成和总体化学反应方程，以理想炸药经典CJ爆轰理论为前提，结合常用程序设计语言和数值计算方法，对混合炸药包括爆速、爆压在内的爆轰性能参数进行计算和分析。

所述混合炸药爆轰参数计算方法，从具体实现过程来看，将计算方法概括为计算数据读写、爆轰参数计算和计算结果处理三个主模块，根据所实现的具体功能的不同，每个主模块下又包含若干不同的子模块。

其中，计算数据读写模块，包括：通用模型参数设置子模块、 Becker-Kistiakowsky-Wilson(以下简写为BKW)状态方程参数读取子模块、周围环境初始物理状态读取子模块、混合炸药组分含量及热力学状态方程参数读取子模块和计算结果数据存储子模块，这些子模块组合在一起，用于实现的功能包括：混合炸药爆轰参数计算时通用计算模型参数、混合炸药元素构成及摩尔含量参数、爆轰产物状态方程参数、周围环境热力学初始状态参数和爆轰产物组分含量计算结果参数的输入、输出、读取和存储。

爆轰参数计算模块，包括：爆轰产物压强-温度-体积及自由能计算子模块、爆轰产物组分含量计算子模块、CJ状态计算子模块、Hugoniot关系计算子模块、绝热等熵膨胀计算子模块、和数值迭代计算子模块。用于实现的功能包括：混合炸药爆轰参数计算时爆轰产物热力学状态及组分含量变化的计算，混合炸药 CJ状态下爆速及爆压的计算，和爆轰产物冲击状态及等熵膨胀状态的计算。

计算结果处理模块，包括：计算结果数据提取子模块和图形化显示子模块。用于实现的功能包括：混合炸药爆轰参数计算结果数据的提取利用、分析处理和图形化显示。

所述混合炸药爆轰参数计算方法，从具体计算流程来看，将计算方法归纳为前处理、数值计算和后处理三个主流程，根据具体操作方式的不同，每个主流程下又包含若干不同的计算步骤。

其中，前处理流程，用于实现的核心功能包括：混合炸药爆轰参数计算时数据输入文件的编辑和构造。

前处理流程，如图1所示，包括以下步骤：

步骤100，数据输入文件编辑：

通过运行辅助程序根据页面提示构造数据输入文件，或者，使用文本编辑工具根据计算要求构造数据输入文件，包括以下步骤：

步骤110，通用计算模型参数编辑：选中并激活进行CJ状态计算、 Hugoniot关系计算、和绝热等熵膨胀计算的计算参数选项。

步骤120，输入BKW状态方程参数，BKW状态方程表示形式如式1所示：

式1中，p为混合炸药爆轰产物的压强，T为混合炸药爆轰产物的温度，V_g为混合炸药爆轰产物中气相组分的体积，R为普适气体常数，x的表达式如式2 所示：

式2中，T和V_g的含义同式1，x_i为第i种气相爆轰产物的摩尔分数，k_i为第i种气相爆轰产物的几何余容，求和只在气相爆轰产物之中进行。式1和式2 中，α、β、θ和κ为BKW状态方程参数。

混合炸药发生爆轰反应后生成组成结构复杂的爆轰产物体系。爆轰产物体系在整体上具有一致的压强p和温度T。从相态角度看，爆轰产物可大致划分为气相组分和凝聚相组分。其中爆轰产物气相组分的物态变化采用BKW状态方程进行描述。BKW状态方程建立在将分子间排斥势应用于Virial(维里)方程基础之上。进行数值计算时，需要知道BKW状态方程参数。经典的BKW状态方程包含两套参数，一组适用于碳原子含量较多的TNT类型的混合炸药，称之为TNT 参数；另一组适用于碳原子含量较少的RDX类型的混合炸药，称之为RDX参数。近些年来，为了得到更加准确的爆轰参数计算结果，以爆轰产物Hugoniot 冲击测试结果和爆轰参数测试结果为对标，参考爆轰产物气相组分的几何余容值，对BKW状态方程参数进行校正和新的标定。

步骤130，输入环境初始状态参数，包括初始温度、初始压强。

步骤140，混合炸药构成参数编辑：输入混合炸药初始密度，输入混合炸药的化学元素构成，输入混合炸药爆轰产物组成组分，根据爆轰产物组成组分的化学元素构成和质量比例分数折合计算出每种化学元素的摩尔含量，输入爆轰产物组成组分的热力学状态方程参数。

进一步的，爆轰产物组成组分的热力学状态方程参数表示为以温度为自变量的多项式拟合函数。进行拟合时，查阅热化学数据表，如JANAF表，以爆轰产物组成组分的热力学熵作为基准，根据热力学熵随温度的变化关系，进行多项式数据拟合。兼顾计算准确度和拟合复杂度，通常可选择带常数的四阶多项式拟合。之后，再在热力学熵的基础上，应用热力学基本关系式进行推导，从而获得包括热力学焓、吉布斯自由能在内的其他热力学量的多项式拟合函数关系。

进一步的，爆轰产物组成组分的热力学状态方程参数还包括组分的生成焓以及气相组分的几何余容。

数值计算流程，用于实现的核心功能包括：爆轰产物组分含量的计算，混合炸药CJ爆轰状态的计算，爆轰产物Hugoniot关系的计算，和爆轰产物绝热等熵膨胀过程的计算。数值计算流程，包括以下步骤：

步骤200，爆轰参数计算，基于多相多组分化学平衡热力学基本原理和CJ 爆轰理论，对混合炸药爆轰性能参数进行热力学计算，包括以下步骤：

步骤210，爆轰产物气相组分自由能的计算，基于BKW状态方程，推导得出的计算表达式如式3所示：

式3中，μ_g为爆轰产物气相组分的化学势，G为爆轰产物组分的吉布斯自由能，ΔH⁰ _f为爆轰产物组分的生成焓，p₀为初始状态压强，下标i表示第i种气相组分，x_i、p、f(x)、β、k_i和k的含义与式1和式2相同。

步骤220，爆轰产物凝聚相组分自由能的计算，爆轰产物凝聚相组分的物态变化采用Cowan状态方程进行描述，如式4所示：

p＝p₁(V_s)+a(V_s)T+b(V_s)T² 式4；

基于Cowan状态方程，推导得出的计算表达式如式5所示：

式4和式5中，p₁、a和b为Cowan状态方程参数，V_s为爆轰产物凝聚相组分体积，μ_s为爆轰产物凝聚相组分的化学势，G和ΔH⁰ _f的含义同式3，p、T 的含义同式1，下标s表示凝聚相组分，上标0表示初始状态。

步骤230，爆轰产物组分含量的计算，根据化学平衡自由能最小原理，化学元素原子数守恒原理和爆轰产物质量守恒原理，计算表达式如式6所示：

式6中，x_i的含义同式1，x_T为当前计算步时爆轰产物组分的总含量，y_i为下一个计算步时第i种气相爆轰产物的摩尔分数，y_T为下一个计算步时爆轰产物组分的总含量，求和符号上标N表示对N种气相爆轰产物求和，λ_k为第k种化学元素的拉格朗日乘子符号，α_ik为系数矩阵，表示第i种爆轰产物中第k种元素的原子数含量，求和符号上标M表示对M种化学元素求和，NT为所有种类的爆轰产物，b_k为混合炸药中第k种元素的原子数含量。

通过步骤210和步骤220，分别得到当前计算步的爆轰产物气相组分和凝聚相组分各自的自由能。根据混合法则，以气相组分和凝聚相组分各自的含量作为相应的权重值，将各自部分的自由能进行线性加权，即可得到当前计算步的混合体系的总自由能。根据式6，以当前计算步的爆轰产物组分含量为初始值，使用数值迭代方法进行计算，得到下一个计算步的爆轰产物组分含量。

步骤240，爆轰产物CJ状态的计算，对混合炸药爆轰产物处于CJ状态时的爆速和爆压进行计算，计算表达式如式7和式8所示：

式7和式8中，D_CJ为CJ状态下的爆速，p_CJ为CJ状态下的压强，V_CJ为 CJ状态下的体积，偏导数符号下标S表示等熵状态，V₀为初始体积，p₀为初始压强，V_g的含义同式1。

进一步的，所述爆轰产物CJ状态的计算包括CJ状态下混合炸药爆速和爆压的计算。

步骤250，爆轰产物Hugoniot关系的计算，结合混合炸药爆轰产物流场的能量守恒关系，得到的计算表达式如式9所示：

式9中，e为爆轰产物总的能量，e₀为初始能量，p、V_g含义同式1，p₀、 V₀含义同式8。

步骤260，计算爆轰产物绝热等熵膨胀过程，计算表达式如式10所示：

S_Total-S_CJ＝0 式10；

式10中，S_Total为膨胀过程中的熵，S_CJ为CJ状态下的熵。S_Total的计算表达式如式11所示：

式11中，x_ig为爆轰产物气相组分的含量，S_g为爆轰产物气相组分的热力学熵，x_is为爆轰产物凝聚相组分的含量，S_s为爆轰产物凝聚相组分的热力学熵。

绝热等熵膨胀过程以CJ状态为起始点，爆轰产物以CJ状态下的热力学熵值保持不变而向周围环境膨胀。

后处理流程，用于实现的核心功能包括：混合炸药爆轰参数计算结果的可视化处理。后处理流程，包括以下步骤：

步骤300，计算结果可视化，包括混合炸药爆轰产物CJ状态、Hugoniot关系、绝热等熵膨胀的图形化显示。

本发明所解决的技术问题涉及到混合炸药爆轰参数的计算方法。发明构思是以化学反应平衡状态自由能最小原理和CJ理想爆轰理论为基础，结合混合炸药配方组成以及化学反应方程式，以BKW状态方程和Cowan状态方程分别描述爆轰产物气相组分和凝聚相组分的物态变化规律，从热力学角度出发，应用数值迭代方法，对包括爆速、爆压在内的混合炸药爆轰参数进行计算。本发明所公开的技术方案中包含了一系列通用化的计算模块和次序化的计算流程，适合于采用常用程序设计语言开发形成的标准化计算程序，以供各种不同种类、不同含量的混合炸药配方爆轰性能设计和预估使用。

本发明的混合炸药爆轰参数计算方法，分别使用BKW状态方程和Cowan 状态方程描述爆轰产物气相组分和凝聚相组分的物态变化规律，基于化学平衡时自由能最小原理，对包括爆速、爆压在内的混合炸药爆轰参数进行数值计算，包括数据输入文件编辑、爆轰参数计算及计算结果可视化步骤。本发明可对不同元素构成及不同组分含量的混合炸药进行反复多次的爆轰参数计算，计算精度满足实际需求，方便用户使用。

以下给出本发明的具体实施例，需要说明的是本发明并不局限于以下具体实施例，凡在本申请技术方案基础上做的等同变换均落入本发明的保护范围。

实施例1：

遵从上述技术方案，本实施例使用混合炸药爆轰参数计算方法对混合炸药 EX1的爆轰参数进行计算。

混合炸药EX1，以质量份数计，配方组成为：黑索今(RDX)：聚异丁烯(PIB)：癸二酸二辛酯(DOP)：油(Oil)＝91：5.3：2.1：1.6。

步骤100，数据输入文件编辑：使用文本编辑工具根据计算要求构造数据输入文件。

步骤110，通用计算模型参数编辑：选中并激活进行CJ状态计算、 Hugoniot关系计算、和绝热等熵膨胀计算的选项。

步骤120，数组BKW状态方程参数，其中α＝0.5，β＝0.16、θ＝400和κ＝10.91。

步骤130，输入环境初始状态参数，初始温度为300K、初始压强为 0.1MPa。

步骤140，混合炸药构成参数编辑：混合炸药的初始密度为1.68g/cm³，按比例折算后，混合炸药EX1的化学式表示为C_1.6234H_3.2468N_2.5659O_2.5800，生成焓为1.4054×10⁵cal/g。混合炸药爆轰产物组分包括气相组分H₂O、CO₂、N₂、H₂、 O₂、CO、NH₃、H、OH和CH₄以及凝聚相组分固相碳(Sol C)，分别从爆轰产物数据库中提取并输入爆轰产物组成组分的热力学状态方程参数。

步骤200，爆轰参数计算。

步骤210，根据式3计算爆轰产物气相组分自由能。

步骤220，根据式4和式5计算爆轰产物凝聚相组分自由能。

步骤230，应用数值迭代算法，根据式6，计算爆轰产物组分含量。

步骤240，爆轰产物CJ状态计算，根据式7计算混合炸药CJ状态下的爆速，D_CJ＝8266m/s；根据式8计算混合炸药CJ状态下的爆压，p_CJ＝29.2GPa。

步骤250，根据式9计算爆轰产物的Hugoniot关系。

步骤260，根据式10计算爆轰产物绝热等熵膨胀过程。

步骤300，计算结果可视化，爆轰产物CJ状态图形显示如附图2所示， Hugoniot关系图形显示如附图3所示，绝热等熵膨胀图形显示如附图4所示。

本实施例能够对C、H、N、O化学元素构成的混合炸药爆轰参数开展计算，计算精度满足工程需求，有利于节约成本、缩短时间、提高效率。

实施例2：

遵从上述技术方案，本实施例使用混合炸药爆轰参数计算方法对混合炸药 EX2的爆轰参数进行计算。

混合炸药EX2，以质量份数计，配方组成为：黑索今(RDX)：梯恩梯(TNT)：铝粉(Al)：蜡(Wax)＝40：38：17：5。

步骤100，数据输入文件编辑：使用文本编辑工具根据计算要求构造数据输入文件。

步骤110，通用计算模型参数编辑：选中并激活进行CJ状态计算、Hugoniot关系计算、和绝热等熵膨胀计算的选项。

步骤120，数组BKW状态方程参数，其中α＝0.5，β＝0.16、θ＝400和κ＝10.91。

步骤130，输入环境初始状态参数，初始温度为300K、初始压强为 0.1MPa。

步骤140，混合炸药构成参数编辑：混合炸药的初始密度为1.72g/cm³，按比例折算后，混合炸药EX2的化学式表示为C_2.0678H_2.6299N_1.5863O_2.0843Al_0.6303，生成焓为4.926×10³cal/g。混合炸药爆轰产物组分包括气相组分H₂O、CO₂、N₂、 H₂、O₂、CO、NH₃、H、OH和CH₄以及凝聚相组分Al₂O₃、固相铝(Al(s))和固相碳(Sol C)，分别从爆轰产物数据库中提取并输入爆轰产物组成组分的热力学状态方程参数。

步骤200，爆轰参数计算。

步骤210，根据式3计算爆轰产物气相组分自由能。

步骤220，根据式4和式5计算爆轰产物凝聚相组分自由能。

步骤230，应用数值迭代算法，根据式6，计算爆轰产物组分含量。

步骤240，爆轰产物CJ状态计算，根据式7计算混合炸药CJ状态下的爆速，D_CJ＝7353m/s；根据式8计算混合炸药CJ状态下的爆压，p_CJ＝23.4GPa。

步骤250，根据式9计算爆轰产物的Hugoniot关系。

步骤260，根据式10计算爆轰产物绝热等熵膨胀过程。

步骤300，计算结果可视化，爆轰产物CJ状态图形显示如附图5所示， Hugoniot关系图形显示如附图6所示，绝热等熵膨胀图形显示如附图7所示。

本实施例能够对C、H、N、O、Al化学元素构成的混合炸药爆轰参数开展计算，计算精度满足工程需求，有利于节约成本、缩短时间、提高效率。

实施例3：

遵从上述技术方案，本实施例使用混合炸药爆轰参数计算方法对混合炸药 EX3的爆轰参数进行计算。

混合炸药EX3，以质量份数计，配方组成为：奥克托今(HMX)：氟橡胶 (F-Rubber)＝95：5。

步骤100，数据输入文件编辑：使用文本编辑工具根据计算要求构造数据输入文件。

步骤110，通用计算模型参数编辑：选中并激活进行CJ状态计算、 Hugoniot关系计算、和绝热等熵膨胀计算的选项。

步骤120，数组BKW状态方程参数，其中α＝0.5，β＝0.16、θ＝400和κ＝10.91。

步骤130，输入环境初始状态参数，初始温度为300K、初始压强为 0.1MPa。

步骤140，混合炸药构成参数编辑：混合炸药的初始密度为1.83g/cm³，按比例折算后，混合炸药EX3的化学式表示为C_1.5108H_2.9213N_2.5659O_2.5659F_0.1003，生成焓为9.221×10³cal/g。混合炸药爆轰产物组分包括气相组分H₂O、CO₂、N₂、 H₂、O₂、CO、NH₃、H、OH、CH₄、HF、CF₄、F₂、COF₂以及凝聚相组分固相碳(Sol C)，分别从爆轰产物数据库中提取并输入爆轰产物组成组分的热力学状态方程参数。

步骤200，爆轰参数计算。

步骤210，根据式3计算爆轰产物气相组分自由能。

步骤220，根据式4和式5计算爆轰产物凝聚相组分自由能。

步骤230，应用数值迭代算法，根据式6，计算爆轰产物组分含量。

步骤240，爆轰产物CJ状态计算，根据式7计算混合炸药CJ状态下的爆速，D_CJ＝8805m/s；根据式8计算混合炸药CJ状态下的爆压，p_CJ＝35.5GPa。

步骤250，根据式9计算爆轰产物的Hugoniot关系。

步骤260，根据式10计算爆轰产物绝热等熵膨胀过程。

步骤300，计算结果可视化，爆轰产物CJ状态图形显示如附图8所示， Hugoniot关系图形显示如附图9所示，绝热等熵膨胀图形显示如附图10所示。

本实施例能够对C、H、N、O、F化学元素构成的混合炸药爆轰参数开展计算，计算精度满足工程需求，有利于节约成本、缩短时间、提高效率。

尽管这里参照本发明的多个解释性实施例对本发明进行了描述，但是，应该理解，本领域技术人员可以设计出很多其他的修改和实施方式，这些修改和实施方式将落在本申请公开的原则范围和精神之内。