固体火箭发动机结构意外爆炸的数值模拟与风险评估❋

樊凯旋 彭 旭 任家帆 饶国宁

南京理工大学化工学院(江苏南京，210094)

引言

近年来，随着我国航天航空技术的不断快速发展，对于航天航空产品作业过程中的安全要求越来越高。 作为航天航空产品的动力装置，固体火箭发动机是不利用周围介质，而是利用自身携带的固体推进剂来生成工质的直接反作用发动机。 固体推进剂是化学推进剂的一种，它利用化学反应所释放出的热能作为推进系统的热源[1]。 随着现代固体推进剂能量的提高和火炸药技术的融合，固体推进剂中会含有高能炸药成分，在外界刺激作用下一旦发生爆炸，将会导致灾难性事故的发生[2-4]。

从20 世纪60 年代起，国外对固体火箭发动机的意外爆炸及其后果开展了深入的研究[5-7]，不仅进行了大量的模拟仿真和理论分析工作，还进行过多次的固体发动机爆炸试验，积累了大量的试验数据，并且依此制定了安全防范措施。目前，国内的研究工作主要集中在几个方面：1)针对固体发动机在外来冲击下是否安全[8]，并未涉及到爆炸所产生的危害；2)固体火箭发动机自毁的危险性研究[9-10]；3)高能固体火箭发动机爆炸冲击波毁伤效应研究[11]，这里只是涉及在无遮蔽的情况下固体发动机爆炸冲击波对人的杀伤作用，并没有考虑到固体发动机在工房中爆炸的情况。 而针对航天企业装配作业时工房内的固体发动机爆炸事故的模拟和后果评估方面鲜有研究报道。

本文中，以航天产品装配工房作业为研究对象，采用显式动力学软件AUTODYN 模拟固体发动机作业过程中的意外爆炸，分析工房内的冲击波流场分布规律以及对人员伤害的影响，研究结果可以为装配作业的风险评估、安全对策与防范提供理论支撑。

1 有限元模型的建立

以某型号固体火箭发动机装配作业为研究对象，只考虑冲击波对人的伤害作用，所以可以将固体发动机等效为相应质量的TNT。 有限元模型如图1所示。

1.1 材料方程及参数的确定

1)空气。 空气采用Ideal Gas 理想气体状态方程，其表达式为

式中：p为压力；ρ为空气密度，ρ ＝1.225 kg/m3；e为气体单位质量内能，e ＝206.8 kJ/g；γ为绝热指数，γ ＝1.4。

2)炸药。 根据固体发动机的TNT 当量，采用JWL 状态方程来描述其爆轰过程，具体表达式为

式中：p为爆压；V为爆轰产物的相对体积；E为初始比内能；A1、B1、R1、R2和ω为材料常数，其具体参数取值如表1 所示。 表1 中，ρ为TNT 的密度；D为炸药的爆速[12-13]。

3)水泥混凝土。 CONC-35MPA 混凝土的状态方程为p-alpha，其表达式为

表1 炸药的JWL 状态方程参数Tab.1 Parameters of JWL equation of state for explosives

式中：ρ0为材料初始密度；ρ为压缩过程中材料的密度；e为初始内能；A1、A2、A3、B0、B1、T1、T2为材料常数，其具体参数取值如表2 所示[14]。

表2 CONC-35MPA 混凝土的材料方程参数Tab.2 Material equation parameters of CONC-35MPA concrete

强度模型则选择适用于描述爆炸过程大变形、高应变的RHT Concrete 本构模型。 RHT 本构模型引入了3 个极限面， 即弹性极限面、失效面和残余强度面，分别描述混凝土的初始屈服强度、失效强度及残余强度的变化规律。 失效面的方程为

表3 RHT 本构模型强度参数Tab.3 Strength parameters of RHT constitutive model

1.2 数值计算模型

AUTODYN 里建立的模型为空气模型、产品模型与工房模型。 气域为长方体结构，将产品和工房包括在内，将三者用完全耦合的方式连接，空气域的边界条件为Flow-out (无限域、无反射边界条件)。

1.2.1 空气模型

空气模型采用Euler，3D-Multi-Material 算法，空气域为100 m×100 m×15 m 的长方体结构，采用六面体结构化网格，将网格设置为渐变网格，最小网格尺寸为50 mm。 网格划分如图2 所示。

1.2.2 产品模型

根据《军工燃烧爆炸品工程安全技术规范征求意见稿2019》可以得出，此型号固体火箭发动机的TNT 当量为0.4[17]，该数值是通过实验得出的。

此产品的TNT 当量比为0. 4，TNT 的密度为1.63 g/cm3，密度体积公式为

根据式(7)可以得到TNT 的尺寸为0.55 m ×0.40 m ×0.40 m。 TNT 采用Euler 算法，单发产品炸点为炸药的中心，双发产品的炸点有两个，分别布置在两个炸药的中心，起爆时间为同时起爆，双发产品间距为1 m，炸药嵌在空气域里。 炸药采用的是六面体结构化网格，网格尺寸为50 mm，单发和双发产品网格划分如图3 和图4 所示。

1.2.3 工房模型

工房为长方体状，尺寸为30 m×15 m×9 m，工房采用六面体结构化网格，网格尺寸为200 mm。 单发产品放置在工房地面的中心处，双发产品在工房地面的中轴线上并排放置。 网格划分如图5 所示。

为了定性和定量地分析工房内部压力分布规律，在空气域内设置了一系列监测点，监测点的位置均为距工房地面1 m 高。 为了便于比较，设置单发产品和双发产品爆炸时监测点的位置相同，如图6所示。

2 数值模拟结果分析

产品在工房地面爆炸时，形成的爆炸冲击波向四周传播，当遇到壁面时反射形成反射冲击波，这些反射波及入射波相互叠加，使空间内的压力增大。TNT 在密闭空间爆轰后，能量的释放主要分为两个过程：爆轰和后燃烧过程。 爆轰形成了冲击波压力，后燃烧形成了密闭空间中的准静态压力，高频冲击波压力波与低频准静态压力波叠加[18]。 这里只研究爆轰产生的冲击波超压的影响，以此来进行定性和定量分析。

2.1 不同时刻的压力分布

单发产品爆炸过程中，不同时刻工房的压力云如图7 所示。

当t＝0 ms 时，由于爆炸还没有开始，工房处在空气域内，所以此时工房内的压力应该与环境压力相同，压力为101 kPa。

当t为10 ms时，窗口处的压力大于工房内其他位置的压力，窗口处的最大压力为1 018 kPa。这是由于：炸药爆炸冲击波以球形的方式向外传播，炸药距离窗口比距离墙面更近，所以冲击波到达得更早，窗口的压力要比其余位置大；泄压时，工房内部的压力从窗口泄放，就会导致窗口附近的压力急剧升高。

当t为20 ms 时，工房顶部的最大压力可以达到2 800 kPa，工房内大部分区域的压力为负值。

当t为30 ms 时，工房内的压力值维持在450 kPa 左右，工房顶部的压力开始下降。

当t为40 ms 时，工房内大部分区域压力在620 kPa 左右，最大压力为5 600 kPa。

t从50 ms 到70 ms，窗口处的压力为工房内的最大压力，最大压力从5 700 kPa 上升到8 000 kPa，工房内出现负压区域。

双发产品爆炸过程中，不同时刻工房的压力云如图8 所示。

模拟未开始时，整个工房的压力与环境压力相同；当t为10 ms 时，工房顶部和窗口处的压力率先到达最大值，最大压力为2 600 kPa；当t为20 ms时，工房顶部的最大压力可以达到3 380 kPa，工房内大部分区域的压力为负值；当t为30 ms 时，工房内的压力为580 kPa，工房顶部的压力开始下降；当t为40 ms 时，窗口处的压力依然是工房内的最大压力，最大压力为9 100 kPa；t从50 ms 到70 ms，窗口处的压力为工房内的最大压力，最大压力从10 730 kPa 上升到13 200 kPa。

2.2 冲击波超压与伤害半径

在整个空气域内要对冲击波超压进行定量分析，就必须在空气域的不同位置添加监测点，得到实时压力变化情况，根据压力变化曲线来判断冲击波超压的大小，再根据超压对人的伤害作用，从而求出伤害半径。 考虑窗口对压力的变化会有影响，所以分两种情况来对单发产品和双发产品爆炸的伤害范围进行分析。 其中，X方向与窗口平行，Y方向与窗垂直。 查阅文献得到冲击波超压对人的伤害作用如表4 所示[19]。

2.2.1X方向上的压力分析

表4 冲击波超压对人的伤害作用Tab.4 Damage effect of shock wave overpressure on humans

X方向上监测点的分布为：沿着工房的中心向X轴分布，所有的点均在与X轴平行的直线上。 测点的位置如表5 所示。

表5 X 方向上测点的位置Tab.5 Location of points on the X direction

各测点的压力曲线见图9、图10。

可以看出，单发产品爆炸时，测点1＃～5＃的压力均从101 kPa开始变化，不同位置的冲击波超压到达时间也不相同，到炸点的距离越远，超压到达的时间也越大。 并且随着到炸点距离的增大，冲击波超压也越来越小。 这也符合冲击波超压作用的原则。

双发产品爆炸时，测点1＃～5＃的压力也是从环境压力开始变化的，压力变化趋势和单发产品爆炸时的压力变化基本相同，但是超压峰值的增长与爆炸产品数量的增长呈现非线性关系。 根据表4 冲击波超压对人的伤害作用，得到单发产品和双发产品爆炸的X方向上的伤害范围，如表6 所示。

表6 X 方向上的伤害半径Tab.6 Damage radius in the X direction m

所以，单发产品爆炸在X方向上无伤害半径为40.0 m 之外；轻伤半径＞29.0 ～40.0 m；重伤半径20.5 ～29.0 m；死亡半径为20.5 m 之内。 双发产品爆炸在X方向上无伤害半径为45.0 m 之外；轻伤半径＞32.0 ～45.0 m；重伤半径24.0 ～32.0 m；死亡半径为24.0 m 之内。

2.2.2Y方向上的压力分析

Y方向上监测点的分布为：沿着工房的中心向Y轴分布，所有测点均在与Y轴平行的直线上。 测点的位置如表7 所示。 各测点的压力曲线如图11、图12 所示。

表7 Y 方向上点的位置Tab.7 Location of the point in the Y direction

可以看出，单发产品爆炸时，测点7＃～12＃的压力也是从环境压力开始变化。 为了更好地比较各个点的冲击波超压的差异以及压力变化规律，这里将测点6＃的压力曲线舍去，因为测点6＃的冲击波超压过大，如果放在同一个坐标系中，很难根据压力曲线对其他点做出比较。

图11 反映出随着到炸点的距离变远，冲击波超压到达的时间随之增大，冲击波超压也在逐渐减小，同样符合冲击波超压作用的原则。

双发产品爆炸时，测点7＃～12＃压力增加的幅度明显要比测点1＃～5＃的压力增加幅度大，这是因为测点7＃～12＃分布于两个炸药连线的中垂线上，由于冲击波是以球形的方式向外扩散，所以测点7＃～12＃受到两股冲击波同时作用。 根据表4 冲击波超压对人的伤害作用，就可以得出在Y方向上的伤害半径，如表8 所示。

表8 Y 方向上的伤害半径Tab.8 Damage radius in the Y direction m

所以单发产品爆炸在Y方向上无伤害半径为43.0 m 之外，轻伤半径＞31.0 ～43.0 m；重伤半径22.0 ～31.0 m；死亡半径为22.0 m 之内。 双发产品爆炸在Y方向上无伤害半径为55.0 m 之外，轻伤半径＞42.5 ～55.0 m；重伤半径32.0 ～42.5 m；死亡半径在32.0 m 之内。

2.3 模拟结果与实验结果的比较

根据阳建红等[20]做的无遮蔽情况下高能推进剂爆炸实验，可以将实验结果与模拟结果作比较，从而印证上述模拟方法是正确的。

具体做法为：由实验中推进剂的质量和推进剂TNT 当量可以计算出TNT 质量。 由于双发产品之间的冲击波会相互影响，这里只对单发产品进行分析。 将文章中有限元模型的工房部分删除，产品质量与实验中换算得出的TNT 质量相同，然后再增加测点，模拟测点位置与试验测点位置相同，其余条件保持不变。 可以得出测点处冲击波超压的数据，如表9 所示。

从表9 可以看出，在距离爆源大于6 m 时，数值模拟结果与实验所测数据误差范围在10%以内；距离爆源越近，爆炸场内冲击波传播越复杂，马赫反射的影响越显著，从而造成了模拟数据与实验数据相差较大。 通过实验结果与模拟结果相比较，说明文中所建立的计算模型与方法是正确的。

3 结论

通过数值模拟的手段研究了航天产品装配工房内单发产品和双发产品的爆炸过程，得出如下结论：

表9 实验结果与模拟结果相比较Tab.9 Result comparison of experiment with simulation

1)根据不同时刻的压力云图可以得到，产品在工房内爆炸要比在空气中爆炸复杂得多，超压峰值的大小与工房的结构有很大的关系。

2)当单发产品爆炸时，冲击波超压在X方向上对人的杀伤距离为29.0 m，爆炸后人员的安全距离为40.0 m；在Y方向上对人的杀伤距离为31.0 m，爆炸后人员的安全距离为43.0 m。 当双发产品爆炸时，冲击波超压在X方向上对人的杀伤距离为32.0 m，爆炸后人员的安全距离为45.0 m；在Y方向上对人的杀伤距离为42.5 m，爆炸后人员的安全距离为55.0 m。