火箭空中爆炸冲击波峰值超压预测方法

2020-07-25王岩王华崔村燕段永胜赵蓓蕾

北京航空航天大学学报 2020年7期

王岩，王华，崔村燕，段永胜，赵蓓蕾

（1.航天工程大学研究生院，北京101416； 2.中国人民解放军96901部队22分队，北京100095；3.航天工程大学宇航科学与技术系，北京101416）

载人航天在空间探索及利用中发挥着越来越重要的作用。然而，对载人航天发射过程中安全性的研究相对匮乏。运载火箭系统是一个技术密集型的复杂大系统，其在载人航天任务中起重要运输作用。火箭在正常飞行过程中，因箭体结构或部件失效发生推进剂泄漏可能会引发爆炸。火箭爆炸事故发生后，乘员舱将受到冲击波超压作用，进而威胁宇航员的生命安全。目前，关于火箭飞行高度对冲击波参数影响的研究十分匮乏。［1］。因此，研究火箭空中爆炸冲击波参数预测方法对于乘员舱的安全评估具有重要意义。

国外研究人员通过开展大量试验以及对试验数据进行分析研究了液体推进剂的爆炸特性［2-4］。Farber［5］通过液体推进剂爆炸试验研究了爆炸当量；Blackwood等［6］建立了基于经验的火箭爆炸模型，对爆炸环境参数进行了预测研究。液体推进剂爆炸试验具有风险大、成本高、难度大、重复性差等特点［7-9］。随着计算机仿真技术的发展，数值仿真已成为解决这类问题的主要手段［10-14］。

谢雪腾［15］利用仿真手段比较了不同海拔高度下的冲击波参数变化规律，定量给出爆炸冲击波传播特性的高原影响程度，并通过海拔200 m和4 500m的静爆试验，证实了高原效应对爆炸冲击波传播特性的重要性。庞春桥等［16］提出了一种高原冲击波参数预测方法，并通过平原和海拔4 500m高原的试验进行了验证，得到了高原环境对冲击波参数的影响。目前，关于火箭飞行高度对爆炸冲击波参数影响的研究尚不足，如何给出一个科学的结论，是一个迫切需要解决的问题。

为了探究火箭空中爆炸时飞行高度对峰值超压的影响，以及找到作用于乘员舱的冲击波预测方法，本文利用显式非线性动态有限元分析软件ANSYS／LS-DYNA，建立了火箭飞行高度为0～20 km时空中爆炸有限元模型。通过改变大气环境参数，对不同飞行高度下火箭空中爆炸进行了数值模拟，定量给出了冲击波峰值超压随飞行高度的变化规律，在此基础上提出了火箭空中爆炸冲击波参数预测公式。研究结果对量化冲击波对乘员舱的危害作用、乘员舱快速危害性评估以及防护研究提供了一定参考。

1 火箭空中爆炸模式当量确定

假设火箭飞行时发生故障引发爆炸，本文选取8种典型爆炸场景建立仿真模型，研究火箭在0～20 km不同飞行高度发生爆炸时，冲击波对乘员舱的危害作用。与以往研究不同的是，本文考虑了不同飞行高度大气参数对冲击波参数的影响。

TNT当量模型是根据等效能量原则将液体推进剂转换成TNT炸药，应用TNT炸药的爆炸规律来预测液体推进剂爆炸规律［17］。由于液体推进剂与固体推进剂的化学反应机理不同，受多种因素的影响，因此，从理论上很难准确获得不同爆炸模式下的当量系数Y。通过对试验数据的统计分析［2-4］，可以估算液体推进剂的爆炸当量。目前，爆炸当量的计算方法主要采用查表法和读图法2种方法。因此，首先分别采用查表法和读图法计算了当量系数，然后根据最小原理［18］确定最终当量系数。根据上述方法确定的火箭在不同飞行高度爆炸时的当量系数如表1所示。表中：H为火箭飞行高度；M1、M2和Mzt分别为一子级、二子级和助推级推进剂质量。

表1 不同飞行高度下的推进剂质量及当量系数Tab le 1 Propellant m ass and equivalent coefficien t atdifferen t flight altitudes

2 火箭空中爆炸有限元模型

建立0～20 km不同飞行高度下火箭空中爆炸有限元模型，通过获取冲击波在乘员舱所在处产生的峰值超压，以量化研究大气参数对作用于乘员舱的冲击波的影响。

2.1 模型建立及网格划分

利用有限元分析软件ANSYS／LS-DYNA建立了包括火箭箭体（包括一子级、二子级和助推器）、等效TNT炸药、空气域3个部分在内的1∶1有限元模型。考虑到模型具有对称性，建模时采用四分之一模型。火箭空中爆炸模式下的有限元模型如图1所示。

为了保证计算精度，经过多次试算，确定出计算空气域为Φ50 m×62 m（半径为25 m，高为-2～60m）的圆柱体计算域。运载火箭贮箱中的液体推进剂等效为圆柱形TNT炸药装药，其尺寸随火箭爆炸场景的不同而变化。定义TNT圆柱体的半径为r0、高度为h0（单位均为m），且长径比设为h0／（2r0）＝1。不同飞行高度火箭空中爆炸下TNT炸药的初始参数如表2所示。表中：mT和z0分别为TNT的质量和质心位置。若以火箭底面的中心为坐标原点，那么乘员舱位于Z＝46.47m平面内。

图1 火箭空中爆炸有限元模型Fig.1 Finite elementmodel of rocket explosion in air

表2 不同飞行高度数值仿真TNT初始参数Tab le 2 Initial param eters of TNT at differen t flight altitudes in num erical sim u lations

为尽可能减少计算量，在保证精度的前提下采用网格渐变法划分网格。为了进行网格无关性验证，选取了6种不同尺寸的网格，对同一数值模型进行了仿真分析，网格1～网格6的网格尺寸逐渐减小。为了观察冲击波峰值超压与网格尺寸的关系，得到了冲击波峰值超压对比曲线如图2所示，P为冲击波峰值超压，R为比例距离。从图2中可以看出，当网格尺寸增大到“网格4”时，冲击波峰值超压增加不再明显，因此选择“网格4”对模型进行网格划分。将空气域划分为3个部分，分别为0～2.7m，2.7～9m，9～25 m，对应的网格尺寸分别为0.2m，0.8m和1.6m。TNT和火箭箭体的网格尺寸分别为0.2 m和0.4 m。爆心位于TNT炸药质心位置，并采用点起爆方式［19］。

图2 不同网格尺寸模型的冲击波峰值超压对比曲线Fig.2 Comparison of shock wave peak overpressure among models with differentmesh sizes

2.2 算法及边界条件

单元类型采用实体164八节点单元。在大变形和非线性结构分析中，为了防止单元变形，TNT和空气域模型采用ALE多物质算法，而火箭模型采用拉格朗日算法，2种网格之间采用流固耦合算法［20］。

此外，为了模拟爆炸冲击波在对称面的传播效果，在对称平面内施加法线方向的位移约束；在空气域的顶面、底面和侧面施加无反射约束，以模拟无限空气域［21］。

2.3 观测点设置

在乘员舱所在平面Z＝46.47m处等间距选取3个观测点A、B、C，以显示不同位置的冲击波超压随时间变化过程以及峰值大小。3个观测点坐标分别是A（0，0，46.47）m、B（1.25，0，46.47）m、C（2.5，0，46.47）m，后处理中取3个观测点的平均值作为最终仿真结果。乘员舱观测点分布如图3所示。

2.4 材料参数和状态方程

图3 乘员舱观测点分布Fig.3 Observation points distribution of crew module

空气［22］采用MAT_NULL材料模型和线性多项式状态方程EOS_LNIEAR_POLYNOMIAL进行描述，TNT［23］采用 MAT-HIGH-EXPLOSLVEBURN材料模型，爆轰压力Pb和单位体积内能E及相对体积V的关系采用JWL状态方程进行描述。箭体铝合金［24］采用Johnson-Cook模型，状态方程采用Gruneisen进行描述。

火箭在飞行中爆炸，由爆轰产物膨胀强烈压缩空气形成爆炸冲击波，因此空气介质参数变化将对冲击波传播产生一定影响。在不同飞行高度下，环境气体的物理特性（密度、压力等）会发生变化，环境气体参数总是随着飞行高度的增加而递减。高空环境与平原有较大不同，空气相对稀薄，通过改变仿真中的大气参数，定量研究飞行高度对爆炸冲击波参数的影响规律。不同飞行高度下的大气参数如表3所示［25］。表中：ρ为大气密度，E0为初始的单位体积内能，其他参数不随飞行高度变化而变化。

表3 不同飞行高度下的大气参数［25］Tab le 3 A tm ospheric param eters at different flight altitudes［25］

3 火箭空中爆炸冲击波传播规律

3.1 模型有效性验证

图4 仿真结果与经验公式平均值对比曲线Fig.4 Comparison of simulation results with average of several empirical formulas

图4给出了火箭在海平面高度爆炸时冲击波峰值超压仿真值与多个经验公式［26-28］平均值的对比曲线。从图中可以看出，两者吻合较好，说明仿真模型适用于火箭爆炸冲击波的研究。

仿真值与经验公式平均值之间存在一定偏差，主要是因为经验公式是根据试验结果拟合而得到的，而仿真模型是通过求解空气动力学方程而得到的相对较为理想的结果。仿真结果与试验条件存在一定的差别，但仿真结果仍在工程允许范围内。

3.2 空中爆炸火箭解体过程

图5给出了火箭空中爆炸不同时刻等效应变云图，再现了火箭箭体在冲击波超压作用下的解体过程［29］。炸药引爆后发生强烈的化学反应，产生的高温高压爆炸产物向四周急剧膨胀扩散，冲击波以近似球面波的形式向外传播。当t＝0.9ms时，爆炸冲击波传播至箭体结构处，使一子级箭体结构开始出现变形和破坏。当t＝2.1ms时，助推火箭结构出现破坏，碎片呈球形向四周飞散；随着冲击波的传播，箭体结构破坏逐渐向上下两侧延伸，且破坏面积逐渐增大。当t＝6.2 ms时，一子级和助推火箭结构几乎全部破坏，且二子级结构开始出现破坏；t＝8.4ms时，冲击波传播至乘员舱位置，并开始作用于乘员舱。由材料力学基本理论可知，冲击波作用于火箭箭体，使箭体所受应力超过了箭体材料的屈服极限，发生了塑性变形和破坏［30］。

3.3 空气中冲击波传播过程

图5 火箭空中爆炸解体过程等效应变云图Fig.5 Effective strain contour of rocket disintegration process during rocket explosion in air

图6 火箭空中爆炸时的冲击波超压云图演化过程Fig.6 Evolution process of pressure contour of shock wave during rocket explosion in air

图6给出了火箭在飞行中爆炸时空气冲击波的超压云图演化过程，再现了冲击波在空气中的传播规律。不同云图代表冲击波传播的典型时刻。当t＝0.9ms时，炸药产生的高温高压爆炸产物向四周急剧膨胀扩散，冲击波以近似球面波的形式向外传播［31］；当t＝2.1ms时，空气冲击波的传播由于受到箭体的阻碍，出现前后冲击波的叠加，球形冲击波开始出现局部变形；当t＝6.2ms时，箭体内部出现冲击波聚焦现象，而已破坏区域形成稳定的冲击波向外传播。最终波形上下呈锥形状，四周呈球形状，高压区出现在上下冲击波锥形段处。当t＝8.4ms时，冲击波传播至乘员舱，将对其结构构成破坏。

从图6中可以看出，流场的发展过程与实际物理过程基本一致。最初冲击波不是球面波，在一定距离处逐渐接近于球面波。空气冲击波与箭体接触瞬间，冲击波阵面处的气流质点受到阻碍作用速度降低，且方向发生一定改变，与后来的冲击波形成叠加，在箭体内部形成增强的冲击波。

4 结果与分析

4.1 飞行高度对火箭空中爆炸冲击波参数的影响

火箭在不同飞行高度都有可能发生爆炸，大气参数是影响爆炸冲击波参数的重要因素。通过有限元仿真获取作用于乘员舱的冲击波峰值超压，以量化研究大气参数对作用于乘员舱的冲击波参数的影响。

图7给出了冲击波峰值超压P与飞行高度H关系曲线，显示了火箭在不同飞行高度发生爆炸时产生的作用于乘员舱的峰值超压。其中，空爆模型考虑了大气参数的变化，而参考模型未考虑大气参数变化。通过对2类模型进行对比，以研究不同飞行高度大气参数对冲击波峰值超压的影响规律。

图7 火箭空中爆炸冲击波峰值超压随飞行高度变化关系Fig.7 Relation between flight altitude and peak overpressure of shock wave generated by rocket exp losion in air

从图7可以看出，2条曲线均随着飞行高度的上升呈下降趋势，但下降幅度有所不同，且飞行高度越高，2条曲线压力值差异越大。参考模型的峰值超压曲线下降幅度较为平缓，其压力值的下降主要是受到燃料消耗的影响。因为随着火箭飞行高度的上升，火箭燃料在不断消耗，导致总爆炸当量的减小，比例距离逐渐增大。根据峰值超压理论公式可知，冲击波峰值超压随着比例距离的增大而减小。而空爆模型的峰值超压曲线下降较为陡峭，这是由于受到燃料消耗、大气参数变化双重作用的结果。一方面，如上所述，燃料消耗使比例距离增大进而使峰值超压减小；另一方面，大气密度等参数随着飞行高度的上升而减小，加速了冲击波在大气中的传播衰减，进一步导致冲击波峰值超压减小。因此，冲击波峰值超压随着飞行高度的上升而下降，且飞行高度越高，下降趋势越显著。

4.2 火箭空中爆炸冲击波参数预测公式修正

由爆炸力学相关理论可知，冲击波峰值超压与比例距离之间的关系遵循式（1）。将多个经验公式取平均值进行拟合，得到的拟合曲线如式（2）所示。式（2）主要用于计算不考虑大气参数影响的冲击波峰值超压。

为了获得火箭空中爆炸冲击波预测公式，定义压强衰减系数为σ＝Pa／P＝aH2+bH+c。Pa为考虑大气参数的冲击波峰值超压，并假设σ与H遵循二次多项式规律，a、b、c为拟合参数。那么，考虑高度效应的冲击波峰值超压预测公式便可得到，即Pa＝P（aH2+bH+c）。

对于火箭空中爆炸，利用仿真结果可以拟合得到压力衰减因子α随飞行高度H的方程，拟合公式中的参数及其标准差见表4。从表中可以看出，3个参数的标准差均足够小，说明方程拟合质量较高。拟合方程的函数图像如图8所示。从图中可以看出，压力衰减系数随飞行高度的上升而逐渐减小。最终得到的火箭空中爆炸冲击波峰值超压预测公式为