矿井瓦斯爆炸是煤矿开采过程中六大灾害之一，爆炸后产生的高温、高压及有毒有害气体对井下人员生命和各种设备造成重大灾难^[1]。当煤矿发生瓦斯爆炸事故时，矿用救生舱可为无法及时逃离或等待救援的矿工提供一个安全的紧急避难所^[2]，救生舱已成为煤矿安全风险预防设施的六大系统之一^[3]。救生舱必须具备良好的抗爆性能和密封性能，才能保证救生舱在瓦斯爆炸事故中具有较好的稳定性。因此，研发抗爆性能强、稳定性好的救生舱，对于提高矿井爆炸事故中逃生人员的生存率具有重要作用。

从安全角度讲，救生舱整体结构必须稳定、灵活。救生舱的关键结构如外壳、舱门、法兰、加强筋的变形和扭曲必须在合理范围。因此，对救生舱的关键部位展开抗爆性分析就显得异常重要。相关学者对救生舱抗爆性能^[4-7]、隔热性能^[8-9]已展开了诸多研究。虽然也探讨了爆炸载荷类型、迎爆面形状、舱体结构及配筋率对救生舱动力响应的影响，然而，关于瓦斯爆炸作用下舱体抗爆性分析的研究较少。而且，大都采用了三角波^[10-13]和梯形方波^[14]甚至脉冲载荷^[15]及TNT当量法^[16]对瓦斯爆炸荷载进行简化，尽管这些研究大大节省了计算成本，但忽略了爆炸冲击波与结构间的流固耦合过程，导致模拟结果与真实值偏离较大；另外，未考虑舱体结构的尺寸效应对舱体抗爆性能的影响。

基于此，本文中先建立救生舱的物理模型、模拟区域及边界条件，建立瓦斯爆炸作用下救生舱的动态数学模型及条件，采用完全流固耦合算法揭示瓦斯爆炸流场与结构的耦合作用，计算并优化该救生舱的结构参数，分析思路拟为矿用救生舱的安全设计提供一种优化方法。

1.   有限元模型建立

1.1   救生舱模型

参考《煤矿井下救生舱抗爆性能数值模拟检测简要规范(建议稿)》及矿用救生舱的实际参数建立物理模型。实体结构如法兰、铰链、加强筋、把手、舱门、逃生门采用实体单元，蒙皮等采用壳体单元。图 1为救生舱的形状和几何结构。救生舱总长度7.2 m，直径1.8 m；救生舱壳体内部纵向均匀分布10根环向加强筋，6根轴向加强筋，各单元通过钢板和螺栓、法兰焊接在一起构成救生舱。

图  1  救生舱网格模型及配筋

Figure  1.  Meshes of rescue capsule and reinforcement

下载: 全尺寸图片 幻灯片

1.2   计算区域及边界条件

目前，主要通过入射冲击波作用在舱体上的压力时程曲线来确定救生舱爆炸载荷^[11-13]。依照《煤矿可移动式硬体救生舱通用技术条件(征求意见稿)》，巷道模型采用半圆拱形，巷道总长度为148 m，高2.6 m、宽3.2 m，为等截面直巷道，救生舱放置于巷道水平128 m处，爆源段长28 m，冲击波传播到舱体前的距离为100 m，巷道物理模型如图 2所示。

图  2  大型实验巷道瓦斯爆炸物理模型

Figure  2.  Physical model of gas explosion in a large tunnel

下载: 全尺寸图片 幻灯片

在瓦斯空气混合区域设置起爆点，能量为10 mJ。救生舱的底部焊接有槽钢，槽钢与滑轨配合，救生舱可沿着滑轨方向轻微滑动，救生舱底部支架固定在巷道底板且初始速度为零，巷道壁面平整、坚固，不考虑巷道壁面的吸能作用。

1.3   状态方程与材料参数

混合气体的主要成分是甲烷。模拟中选用浓度为9.5%的瓦斯/空气混合气体。采用空物质模型及线性多项式状态方程^[17]描述空气和瓦斯的流动状态，线性多项式状态方程如下：

式中：ρ₀为初始密度, ρ为密度, E为比内能；C₀~C₆为参数, μ=ρ/ρ₀－1。空物质材料模型和线性多项式方程中的参数值，见表 1。

表  1  线性多项式方程中的参数

Table  1.  Parameters used in linear polynomial equation of state

材料	ρ/(kg·m-3)	C0/Pa	C1/Pa	C2/Pa	C3/Pa	C4	C5	C6	E/(kJ·m-3)	V0
空气	1.29	-1×105	0	0	0	0.4	0.4	0	250	1
瓦斯	1.23	0	0	0	0	0.274	0.274	0	340	1

下载: 导出CSV 

| 显示表格

采用非线性塑性材料模型Plastic-Kinematic来模拟舱体各部位受力变形，该模型假定的材料屈服应力为^[18]：

式中：C、P为Cowper-Symonds常数，C=40.4，P=5，β为硬化参数，当β=0时为随动硬化，当β=1时为各向同性硬化，σ₀为屈服应力，ε′为应变率, ε_p^eff为有效塑性应变，σ_y为屈服强度，E_p为塑性硬化模量。Q460钢的状态方程参数分别为：σ_y=460 MPa，E_p=1 180 MPa，E=201 GPa，ν=0.3，β=0，ε_f=0.32。

2.   结果与分析

2.1   爆炸冲击波超压验证

吴斌^[19]和田志敏等^[20]在全尺寸巷道中展开了瓦斯浓度9.5%、体积200 m³甲烷/空气预混气体的爆炸实验，从理论上分析了瓦斯爆炸超压与距离、隧道断面面积和瓦斯体积(初始爆炸能量)之间的非线性关系。冲击波超压计算公式^[21]如下：

式中：c₀为空气传播速度，k为绝热压缩系数，x为爆炸中心与冲击波面之间的距离，S为巷道的横截面积, E_ρ0为瓦斯爆炸释放的总能量。

本文中采用ANSYS/LS DYNA软件的模拟结果与文献[2, 19, 21]中结果的对比，如图 3所示。

图  3  瓦斯爆炸实验值、理论值及模拟值对比

Figure  3.  Comparison of experimental, theoreticaland numerical data

下载: 全尺寸图片 幻灯片

由图 3通过比较，本文的爆炸冲击波超压模拟数据与参考文献报道的结果基本吻合。由于巷道壁面粗糙且具有吸能作用，因此数值模拟结果高于理论值和实验值。

2.2   爆炸流场与结构的流固耦合效应

图 4为瓦斯爆炸冲击波流场与救生舱的相互作用，表明爆炸流场与舱体结构之间相互影响。在瓦斯爆炸瞬间，爆炸产物迅速发生膨胀，致使高温高压气体高速向外扩散，同时周围气体在高压的作用下从巷道一端传播至另一端。在287.22 ms时刻，冲击波与救生舱相遇前，入射波超压值约为0.29 MPa，见图 4(a)。在293.68 ms时刻，冲击波与救生舱接触后，在救生舱对爆炸冲击波的激励作用及入射波与反射波的叠加综合作用下，救生舱表面及周围产生一个速度梯度场，冲击波超压迅速增大，达到0.67 MPa，见图 4(b)。由于救生舱的吸能作用及冲击波的多次反射和衍射，冲击波超压迅速衰减。在304.32 ms时刻，压力急剧下降至0.19 MPa，而救生舱的压力不超过0.06 MPa，见图 4(c)。在334.34 ms时刻，冲击波压力峰值衰减至0.038 MPa，见图 4(d)。

图  4  爆炸流场与结构的耦合作用

Figure  4.  Coupling interaction between flow field and structures

下载: 全尺寸图片 幻灯片

2.3   原模型动力响应

2.3.1   应力场

图 5是原模型救生舱的应力云图，从图中可以看出，随着爆炸冲击波与救生舱的相互作用，舱体的等效应力逐渐增大。在287.22 ms时刻，冲击波到达救生舱的表面并开始对救生舱施加压力。在冲击波与整个救生舱作用过程中，舱门的等效应力较高，而侧壁和顶部的应力较小。高应力出现在前舱门与正面的铰链处，最大应力单元719 321的应力为340 MPa，未超过材料的屈服强度，然而舱体的其他部位如法兰等结构单元应力没有超过材料的屈服强度，如图 5(d)所示。当瓦斯爆炸超压达到0.69 MPa，整体结构和主要部件的强度在弹性范围内。由于爆炸冲击波与救生舱的相互作用时间较短，应力集中区不会导致整体失效。

图  5  不同时刻救生舱的等效应力分布

Figure  5.  Distribution of equivalent stress of rescue capsule at different moments

下载: 全尺寸图片 幻灯片

2.3.2   位移场

原模型救生舱的位移云图如图 6所示，由于迎爆面先受到爆炸冲击波和动压作用，因此在这个区域的变形最明显。初始时刻287.22 ms时，救生舱的舱门边缘变形较严重，最大位移达到25 mm，直接影响舱体密闭性能(参见《煤矿可移动式硬体救生舱通用技术条件(征求意见稿)》)，因此需要对迎爆面进行优化；301.66 ms以后，由于冲击波的传播和衍射作用，受力区域逐渐扩大，舱体正面的边缘区域变形较严重，达到25.1 mm。因此，为了保证舱体的密闭性能，需要对迎爆面进行加固处理，以防止瓦斯爆炸产生的有毒气体入侵。变形相对较小的部分，如逃生门、底座支架的位移均小于4 mm，且没有出现局部的脆性断裂和裂缝，这说明了所设计的救生舱其他部位的安全性和整体刚度符合安全要求。

图  6  不同时刻救生舱的位移分布

Figure  6.  Distribution of displacement of rescue capsule at different moments

下载: 全尺寸图片 幻灯片

结合上述分析，冲击波对救生舱各部位的影响差异较大。迎爆面边缘区域是保证救生舱气密性的关键部位。因此，需要对迎爆面进行优化如增加加强筋数量，保证舱体结构在质量较轻的情况下迎爆面的强度和刚度满足要求，提高该区域的抗冲击性能。

3.   结构设计优化与数值分析

3.1   模型优化

结合上述分析可知，选用Hyper Optistruct模块对迎爆面结构进行拓扑优化，获得加强筋的最佳布设位置，结合原模型的变形特征，对原模型迎爆面进行配筋优化后的模型如图 7所示。

图  7  优化模型网格及配筋

Figure  7.  Meshes of optimization model and its reinforcement

下载: 全尺寸图片 幻灯片

3.2   新模型结构响应

3.2.1   应力场

图 8是优化模型的应力云图，从图中可以看出，随着爆炸冲击波与救生舱的相互作用，舱体的等效应力逐渐增大。在287.22 ms时刻，冲击波到达救生舱的表面并开始对救生舱施加压力。在冲击波与整个救生舱作用过程中，舱门的等效应力较高，而侧壁、底部和顶部的应力较小。高应力出现在迎爆面加强筋处，最大应力为436 MPa，但未超过材料的屈服强度，如图 8(d)所示，而舱体的其他部位如法兰等结构的单元应力未超过材料的屈服强度。整体结构和主要部件的强度在弹性范围内，且由于爆炸冲击波与救生舱的相互作用时间较短，应力集中区不会导致整体失效。

图  8  不同时刻优化模型等效应力分布

Figure  8.  Distribution of equivalent stress of optimization model at different moments

下载: 全尺寸图片 幻灯片

3.2.2   位移场

优化模型救生舱的位移云图如图 9所示，由于迎爆面首先受到爆炸冲击波和动压作用，因此在这个区域的变形最明显。初始时刻287.22 ms时，救生舱的舱门边缘变形较严重，最大位移达到18.3 mm，此后随着冲击波的传播，该值逐渐减小，并在《煤矿可移动式硬体救生舱通用技术条件(征求意见稿)》的合理范围内。在301.66 ms时刻，由于冲击波的传播和衍射作用，受力区域逐渐扩大，舱体正面的边缘区域变形较严重，最大位移达到18.6 mm(小于极限变形值20 mm)，还能够有效避免瓦斯爆炸产生的有毒气体入侵。变形相对较小的部分，如逃生门、底座支架的位移均小于4 mm，未出现局部的脆性断裂和裂缝，表明优化后的救生舱的安全性和整体刚度符合安全要求。

图  9  不同时刻优化模型位移分布

Figure  9.  Distribution of displacement of optimal model at different moments

下载: 全尺寸图片 幻灯片

3.3   原、新模型对比

3.3.1   救生舱关键结构动态响应

原模型及优化模型的关键结构如舱门、逃生门、加强筋、连接法兰的最大单元的应力和位移曲线对比结果，如图 10所示。图中，d_{m, or}、σ_{m, or}为原模型的最大位移、最大应力，d_{m, op}、σ_{m, op}为优化模型的最大位移、最大应力。在瓦斯爆炸冲击作用下，迎爆面的动力响应更强烈，因此具有最大的等效应力和位移，原模型与优化模型等效应力均未超过材料的屈服强度，但优化后最大单元位移减小了近61%且未超过20 mm，见图 10(a)。由于爆炸冲击载荷衰减较快且作用时间极短，在救生舱末端处的逃生门的应力和位移相对较小，在280~305 ms内，配筋优化对救生舱的最大单元应力和最大单元位移基本无影响，在305 ms后，二者在最大单元应力和最大单元位移的差异逐渐扩大，这是由于经过配筋优化后的新模型对冲击波超压具有更强的承受能力。此外，优化模型的逃生门最大单元应力和最大单元位移出现时刻均晚于原模型，见图 10(b)。优化模型的加强筋的最大单元应力和最大单元位移均大于原模型，由于爆炸冲击波先作用在舱体迎爆面，而原模型由于迎爆面产生应力集中致变形较大(超过20 mm)，优化后舱体的迎爆面的应力转移至该区段加强筋，且该区段加强筋吸收了迎爆面外壳的冲击波能量，因此优化后的救生舱迎爆面加强筋的应力和位移均大于原模型，见图 10(c)。而优化模型与原模型的法兰的残余位移均小于3 mm，见图 10(d)，应力远小于材料的屈服强度。综合分析结果表明，对原模型迎爆面壳体进行配筋处理对该区段壳体具有一定保护作用，且优化后的舱体强度和刚度均优于原模型。

图  10  原模型与优化模型关键结构动态响应对比

Figure  10.  Comparison of dynamic response of key components on original and optimized models

下载: 全尺寸图片 幻灯片

3.3.2   舱体能量

图 11为原模型与优化模型的动能、内能及总能量能量变化曲线。

图  11  原模型与优化模型的动能、内能及总能量对比

Figure  11.  Comparison of kinetic energy, internal energyand total energy of original and optimized models

下载: 全尺寸图片 幻灯片

能量变化是模型在载荷作用下的一个明显特征，内能包括弹性应变能、黏弹性或者蠕变过程的能量耗散以及伪应变能等，动能主要包括舱体振动能量，而总能量能为动能与内能之和。从图中可以看出，内能与动能变化趋势相反，二者在数值上相差1个数量级，其中的差包括能量耗散及伪应变能等。当爆炸冲击载荷作用于舱体，动能迅速增大，287 ms时达到最大值，随着冲击波的传播，动能逐渐转化为内能，在290 ms以后，动能急剧减小，随着冲击波的传播，舱体动能逐渐减小并趋于平衡。当冲击波完全淹没救生舱舱体时，舱体几乎处于静止状态，动能最小，内能达到最大值。原模型与优化模型的动能、内能及总能量变化趋势基本一致。二者的动能差异基本不大，内能相差近280 kJ。原模型的总能量峰值为1 400 kJ，而新模型总能量峰值为1 150 kJ，比原模型有所衰减。这是由于：原模型迎爆面主要承受爆炸冲击波作用，主要吸能构件是主舱门和舱门，导致这些构件的破坏较严重，而增加了加强筋的改进舱的主要吸能构件是加强筋，加强筋破坏较严重，具有很好的缓冲吸能作用，主舱门吸收能量较小。从上面的能量分析进一步证明，增加迎爆面的加强筋对爆炸冲击波具有一定缓冲作用。

4.   结论

采用ANSYS/LS-DYNA考察了圆柱壳救生舱在瓦斯爆炸载荷下救生舱的动态响应，并进行结构参数优化，结论如下。

(1) 采用ALE算法揭示了爆炸冲击波流场与救生舱的耦合作用，瓦斯爆炸冲击波的数值模拟结果与文献数据基本吻合。

(2) 冲击波与救生舱相互作用过程中，应力集中区位于前舱门与迎爆面的铰链处，原模型最大绝对位移处于救生舱迎爆面边缘。舱门的等效应力较高，而侧壁和顶部的应力较小。舱体的其他部位如法兰等结构单元应力没有超过材料的屈服强度且没有塑性变形，整个舱体处于弹性状态。但迎爆面变形较大，需要对原模型进行优化处理。

(3) 采用Hyper Optistruct对原模型迎爆面结构进行拓扑优化，优化模型抗爆性模拟结果表明：舱体正面的边缘区域变形较严重，最大位移达到18.6 mm(小于极限变形值20 mm)，还能够有效避免瓦斯爆炸产生的有毒气体入侵。其他部位未出现局部的脆性断裂和裂缝，表明优化后的救生舱的安全性和整体刚度符合安全要求。

(4) 从能量转换角度阐释了救生舱优化前后能量变化规律，表明舱体能量遵循一定规律，动能与内能变化趋势基本相反；原模型与优化模型的动能、内能及总能量变化基本一致。原模型与优化模型的动能差异基本不大，内能相差近280 kJ，表明对原模型迎爆面进行配筋优化可以减缓冲击波对迎爆面和舱门的破坏作用。