JWL状态方程及其等效多方状态方程在内爆炸计算中的应用分析

徐维铮，吴卫国

1武汉理工大学高性能舰船技术教育部重点实验室，湖北武汉430063

2武汉理工大学交通学院，湖北武汉430063

0 引 言

从战例统计和试验资料来看，反舰导弹已成为舰船生命力的主要威胁，舰船在舱内爆炸冲击载荷作用下的毁伤效应分析及防护结构设计是目前舰船工程领域的研究热点。对舱内爆炸冲击载荷的准确分析是结构响应评估和防护结构设计载荷输入的前提，因此起着至关重要的作用。数值模拟作为一种研究内爆炸载荷的有效方式，得到了国内外学者的普遍重视。

针对约束空间内爆炸过程的数值模拟，国内外学者开展了许多研究。曹玉忠等［1］采用瞬时爆轰后的高温高压气体作为爆炸初场，其与周围空气均采用理想气体状态方程来描述，数值模拟研究了半球顶圆柱筒密闭式抗爆容器的内部爆炸流场。Benselama等［2］基于瞬时爆轰假定，采用理想气体状态方程或JWL状态方程（Jones-Wilkins-Lee Equation of State，JWL-EOS）描述爆轰产物，对比分析了数值模拟结果与空中爆炸（空爆）的超压实测结果，结果表明，爆轰产物的状态方程对于空爆超压载荷的影响较小。任会兰等［3］采用基于消息传递接口（MPI）模式开发的内爆炸过程三维并行程序，对不同药量的炸药在典型工事中的爆炸过程进行了数值研究。在模拟过程中，假定炸药瞬时爆轰，其爆轰产物采用可变指数多方气体状态方程进行描述。Ma等［4］基于瞬时爆轰假定，采用可变指数多方气体状态方程来描述爆轰产物，通过求解非守恒形式的可压缩欧拉方程，数值模拟研究了爆炸波在地铁站内的传播过程。Sugiyama等［5］基于瞬时爆轰假定，采用理想的气体状态方程描述了爆轰产物及周围空气，全场具有统一的气体绝热指数1.4，在研究中，对二维轴对称可压缩欧拉方程进行数值离散，研究了房间内爆炸泄出过程中超压峰值随传播距离衰减的特性。近期，徐维铮等［6-8］基于瞬时爆轰假定，同样采用理想气体状态方程描述爆轰产物及周围空气，使用自主开发的三维内爆炸波高精度数值计算程序，探讨了泄压口大小、WENO格式精度、炸药质量等因素对约束空间内爆炸准静态超压载荷的影响规律。

从上述研究现状可知，描述炸药瞬时爆轰后的爆轰产物状态的方程主要有3种：JWL状态方程、可变指数多方气体状态方程和理想气体状态方程。其中，JWL状态方程既可用来描述爆炸冲击载荷的高压段，也可描述低压段，并且已被目前的商用程序（如AUTODYN，LSDYNA等）所采用；可变指数多方气体状态方程是以一种随密度改变绝热指数的方式来近似描述爆轰产物状态的方程，目前已被相关学者采用［3-4］；虽然上述研究并未给出确定其中参数的方法，且理想气体状态方程在描述高压段时也有一定的误差，但作为一种简单的状态方程，它便于程序编写，因此有学者也将其应用到了内爆炸场的数值研究中。

在真实内爆炸过程中，针对不同炸药的JWL状态方程，其ω参数存在一定的差异，并对爆炸过程会产生一定的影响。当爆炸发生在不同气体环境时，初始气体绝热指数γ0的差异也会对爆炸过程产生影响。因此，可变指数多方气体状态方程作为描述爆轰产物状态的一种等效方式，具备理想气体状态方程的简洁性，且便于程序编写及数值计算，是一种较有优势的状态方程。

鉴此，本文将基于Fortran平台，采用五阶WENO有限差分格式，利用自主开发的舱室内炸药爆炸过程二维数值计算程序，研究参数ω及γ0在不同取值下对内爆炸参数计算的影响规律，以此提出用于确定等效多方气体状态方程中参数的方法。

1 控制方程及数值求解

本文将整个爆炸场视作由爆轰产物及空气两种成分组成，将爆轰产物质量分数方程与可压缩欧拉方程进行耦合求解，具体形式如下［9］：

其中，

式中：U为守恒变量；G，H分别为x,y方向的数值通量；ρ为气体密度；u,v分别为x,y方向气体的速度分量；p为气体压力；E为单位体积流体的总能量；e为气体内能；λ为爆轰产物质量分数。

由于式（1）中含有ρ，u，v，p，e，λ这6个未知变量，所以还需补充状态方程才能将方程组封闭。在本文计算中，爆轰产物的状态采用JWL状态方程［2］描述：

式中：A，B，R1，R2，ω均为JWL状态方程中的参数，其中A=3.712×1011Pa，B=3.232×109Pa，R1=4.15，R2=0.95，ω=0.3，ρE为炸药密度，ρE=1 630 kg/m3。

此外，周围气体采用理想气体状态方程描述：

式中：γ0=1.4。此时初始气体密度ρ0=1.0 kg/m3；初始大气压力p0=1.0×105Pa。

本文分别采用五阶WENO有限差分格式［10］和三阶TVD龙格—库塔（TVD-RK）法对上述方程中的空间项及时间项进行数值离散［7］。

2 内爆炸数值计算

2.1 舱室尺寸及测点布置

本节主要讨论状态方程（4）及式（5）中ω和γ0参数的不同取值对内爆炸参数计算的影响规律。图1（a）所示为用于计算的长方形舱室尺寸，其中舱内设置有3个测点，通过测点输出爆炸参数的时程曲线。舱室壁面边界条件设为刚性边界，不考虑舱室结构的变形［11］。均匀网格尺寸取为10 mm，如图1（b）所示。

图1 舱室尺寸及网格划分Fig.1 Sizes of cabin and mesh partition

2.2 爆炸初场

基于瞬时爆轰假定［5］，本文将初始方形炸药（边长l=40 mm）等效为均匀的高压气团，如图2所示。具体参数按照JWL状态方程计算如下：，其中气体初始内能e0等于炸药爆热值Q，e0=4.69×106J/kg。炸药周围气体参数由理想气体状态方程计算，即ρ=ρ0,p=p0。

图2 爆炸初场（ω=0.3,γ0=1.4）Fig.2 Initial condition for the explosion（ω =0.3,γ0=1.4）

2.3 参数ω的影响

由于炸药类型不同，ω的值也会有所不同，所以探讨ω的取值对爆炸过程计算的影响规律具有一定的工程意义。本文将式（4）右端前两个指数高压项简记为A+B项，其随气体密度的变化曲线如图3所示。

图3 指数高压项随气体密度衰减曲线Fig.3 Curves of the exponential high pressure term attenuated with gas density

由图3可知，随着爆轰产物密度的降低，高压项数值迅速减小，当爆轰产物密度膨胀为初始密度的1/12时，高压项数值趋于0。即爆轰产物膨胀到一定程度后，爆轰产物JWL状态方程演化为理想气体状态方程形式：

且等效气体绝热指数γe=ω+1。

由式（4）及式（6）可知，在 JWL状态方程中A，B，R1，R2，ρE这5个参数存在于前两个指数高压项中，因在爆炸过程中高压项迅速衰减，故本文不对这5个参数的影响进行探讨。鉴于参数ω既可影响爆炸冲击载荷高压段的计算，也可影响其低压段的计算，且不同炸药类型有不同的ω值［12］，因此，为探讨参数ω对整个爆炸过程计算的影响规律，本文算例选取了3种ω值（ω=0.2，0.3，0.4）进行计算，同时保持γ0=1.4 。

图4所示为不同ω取值下炸药中心测点3处的超压衰减时程曲线。由图可知，随着ω值的增大，炸药内部的压力衰减速率减缓。

图5所示为不同ω取值下测点1和测点2处爆轰产物质量分数λ的时程曲线。由图可知，ω的取值对于不同测点处的爆轰产物质量分数时程曲线有较大影响，且随着ω值的增大，爆轰产物到达舱壁面的初始时间提前，亦即ω值的增大加快了爆炸初期爆轰产物的膨胀速率。

图5 不同ω取值下测点处的爆轰产物质量分数时程曲线Fig.5 Time history curves of mass fraction of detonation products at gauging points for different values ofω

为进一步分析爆轰产物质量分数的最终演变趋势，通过质量加权，可得到充分混合后的舱内爆轰产物的质量分数λf为

式中：m为炸药质量；mkong为空气质量；V为舱室体积。爆炸初期，爆轰产物的质量分数为1.0，周围环境气体对应的爆轰产物质量分数为0。

由式（7）可知，λf仅与炸药质量、舱室体积相关，而与ω的取值无关，即舱内爆轰产物的质量分数将最终趋于一致。从图5也可见这种趋势。

图6所示为不同ω取值下测点1和测点2处气体绝热指数γ的时程曲线。由图可知，ω的取值对气体绝热指数的时程有着显著影响，即随着ω的增大，最终混合物的气体绝热指数也增大。

图6 不同ω取值下测点处的气体绝热指数时程曲线Fig.6 Time history curves of gas specific heat ratio at gauging points for different values ofω

根据式（7）的分析思路，可得到舱内最终混合物的气体绝热指数γf为

式（8）的推导用到了式（6）的结论。由式（6）可知，当爆轰产物膨胀到后期，JWL状态方程演化为理想的气体状态方程形式，此时爆轰产物的绝热指数为γe，舱室内最终混合物的气体绝热指数γf为爆轰产物与空气绝热指数的质量加权之和。因此，由式（8）可知，当其他参数不变时，γf会随参数ω的增大而增大。

图7所示为不同ω取值下测点处的超压时程曲线。由图可知，随着ω的增大，冲击波峰值增大，冲击波到达舱壁面的初始时间提前，且最终形成的准静态压力峰值也随之增大。

图7 不同ω取值下测点处的超压时程曲线Fig.7 Time history curvs of overpressure at gauging points for different values ofω

为了说明准静态压力峰值增大的原因，本文采用了文献［6］中的密闭空间内爆炸准静态超压峰值Δps计算公式：

然后，进一步将式（8）代入式（9），可得

由式（10）可知，当其他参数不变时，Δps将随着参数ω的增大而增大。

2.4 参数 γ0的影响

当爆炸发生在不同气体环境（例如氮气、氧气、空气）中时，不同气体有着不同的初始绝热指数γ0［13］，并对爆炸过程产生一定的影响，故探讨γ0的取值对于研究其对爆炸参数计算的影响规律具有一定的工程意义。本文选取3种γ0值（γ0=1.2，1.3，1.4），同时保持ω=0.3。

图8所示为不同γ0取值下测点处的爆轰产物质量分数时程曲线。由图可知，γ0的取值对不同测点处爆轰产物质量分数时程曲线有着较大影响，且随着γ0值的减小，爆轰产物到达舱壁面的初始时间提前，即γ0值的增大可使爆炸初期爆轰产物的膨胀速率减缓。

图8 不同γ0取值下测点处的爆轰产物质量分数时程曲线Fig.8 Time history curves of mass fraction of detonation products at gauging points for different values ofγ0

图9所示为不同γ0取值下测点处的气体绝热指数时程曲线。由图可知，γ0的取值对气体绝热指数γ的时程有着显著影响，随着γ0的增大，最终混合物的气体绝热指数γf增大。由式（8）也可知，当其他参数不变时，γf随着参数γ0的增大而增大。

图9 不同γ0取值下测点处的气体绝热指数时程曲线Fig.9 Time history curves of gas specific heat ratio at gauging points for different values ofγ0

图10所示为不同γ0取值下测点处的超压时程曲线。由图可知，随着γ0的增大，冲击波峰值增大，冲击波到达舱室壁面的初始时间提前。由图还可看出，γ0的取值对最终形成的准静态压力峰值的影响很小，具体可由式（10）判断。由式（10）可知，当其他参数不变时，Δps与γ0呈非线性的关系，乘号两端均包含γ0且具有相反的变化特征，从而使得准静态压力峰值的变化并不显著。

3 等效多方气体状态方程

图10 不同γ0取值下测点处的超压时程曲线Fig.10 Time history curves of overpressure at gauging points for different values ofγ0

由于JWL状态方程中的参数较多且包含指数的计算，因此使程序变量增多、计算时间增加。为提高工程计算的效率，作为一种可靠的方式，可以采用一种等效的多方气体状态方程来描述爆轰产物的状态。例如，文献［3-4］在研究炸药爆炸数值计算时，采用了可变绝热指数多方气体状态方程来描述炸药爆轰产物，具体形式如下：

式中：γb0为爆轰产物的初始等效绝热指数；γbf为爆轰产物膨胀到最后的绝热指数；下标b为爆轰产物绝热指数随膨胀过程的衰减系数。

3.1 参数取值

尽管文献［3-4］提出了采用可变绝热指数多方气体状态方程来描述炸药爆轰产物，但却未对炸药爆炸过程影响较大的参数γb0，γbf及b的取值问题进行充分探讨。实际上，如何合理地确定这3个参数值十分重要。因此，本文拟基于JWL状态方程来探讨上述3个参数的取值方法，提出一种与JWL状态方程近似等效的多方气体状态方程，并通过相关数值算例来验证其可行性及可靠性。

本文研究时，γb0的取值由式（12）计算。假设炸药瞬时爆轰，则ρ=ρE，e=Q，将其代入式（4）和式（11）中，按照压力等效原则，可得

式中，γbf的取值由式（13）计算。当爆轰产物得到充分膨胀后，产物的状态方程退化为式（6），并按照压力等效原则，可得

式中，b的取值直接影响到气体绝热指数随密度变化的关系（主要在高压段起作用），目前，暂无一种解析方式可直接确定其数值，但可通过数值计算结果与JWL状态方程的计算结果进行对比，来近似确定其数值。

3.2 数值计算结果分析

本文选取2.2节的炸药参数（边长l=40 mm）进行数值计算。图11所示为不同b取值下测点处的超压时程曲线与JWL状态方程计算结果的对比。由图可知，当参数b取为2 000时，数值计算得到的超压载荷与JWL状态方程的计算结果吻合较好。

图11 不同b取值下测点处的超压时程曲线（l=40 mm）Fig.11 Time history curves of overpressure at gauging points for different values ofb（l=40 mm）

图12所示为不同b取值下测点处爆轰产物的质量分数时程曲线，图13所示为不同b取值下测点处气体绝热指数时程曲线（两图中方形炸药边长l=40 mm）。由图12和13图可知，b的取值对爆炸参数时程曲线影响较大，当b=2 000时，其在高压段（爆炸初期）与JWL状态方程的计算结果更吻合。因此，针对该爆炸工况，b=2 000是一种合理的取值。

图12 不同b取值下测点处爆轰产物的质量分数时程曲线Fig.12 Time history curves of mass fraction of detonation products at gauging points for different values ofb

图13 不同b取值下测点处的气体绝热指数时程曲线Fig.13 Time history curves of gas specific heat ratio at gauging points for different values ofb

为了验证参数b=2 000同样适用于其他药量，本文在其他参数不变的情况下，选取了边长l=160 mm的方形炸药再次进行数值计算。图14给出了在此炸药参数及不同b取值下两个测点处的超压时程曲线，图15给出了不同b取值下测点2处的爆轰产物质量分数时程曲线及气体绝热指数时程曲线。由对图14和图15的分析可知，针对大质量炸药内爆炸工况，参数b=2 000的数值计算结果与JWL状态方程的计算结果同样吻合较好。

图14 不同b取值下测点处的超压时程曲线（l=160 mm）Fig.14 Time history curves of overpressure at gauging points for different values ofb（l=160 mm）

图15 不同b取值下测点2处的爆轰产物质量分数及气体绝热指数时程曲线Fig.15 Time history cures of mass fraction of detonation products and gas specific heat ratio at gauging point-2 for different values ofb

4 结 论

通过本文的研究，得到如下主要结论：

1）在JWL状态方程中，参数ω值的增大可使爆炸初期爆轰产物的膨胀速率加快，最终增大了混合物气体绝热指数和冲击波峰值，使冲击波到达舱室壁面的初始时间提前，且最终形成的准静态压力峰值也增大。

2）初始气体绝热指数γ0的增大可使爆炸初期爆轰产物的膨胀速率减缓，增大了最终混合物气体绝热指数和冲击波峰值，且冲击波到达舱室壁面的初始时间提前，但对最终形成的准静态压力峰值影响很小。

3）本文用来确定等效多方气体状态方程参数的方法具有可行性和可靠性，通过数值模拟结果分析，初步确定了其中的参数b=2 000。