1 引言
对于水面舰船而言，各种穿甲、半穿甲战斗部均可穿过船体外板，进入舱室内部爆炸。由于舱室的密封性，内部爆炸对舱室结构造成的毁伤效应更大。舱内爆炸对舱室结构造成的毁伤模式直接影响着舰船能否继续作战和生存[1-2]。
目前关于内爆毁伤模式的研究主要靠理论分析、数值仿真和模型试验等手段。内爆载荷是非线性载荷，直接使用基础物理理论分析较为困难，使用量纲分析能够定性研究舱室内爆的损伤[3]。姚术健[4]提出内爆分析的无量纲损伤数以及内爆载荷下箱型结构毁伤模式的快速预测方法。郑成[5]通过研究得出了舱壁变形与板厚比和无量纲损伤数之间存在线性关系。焦晓龙[6]在前两者的研究基础上，使用数值仿真与量纲分析结合的方式，提出了单次爆炸下多舱室结构的毁伤等级计算公式。
内爆的药量、炸点位置、加强结构和泄压口的大小等加载条件对于舱壁的毁伤程度有重要影响，有不少学者对其进行了研究。徐维铮等[7]研究了泄压过程中准静态压力的计算方法，并发现炸药对等熵指数的影响要大于泄压口的影响。汪维等[8]对多舱室的内爆进行试验，得出随着装药量的增加，刚箱体内变形程度逐渐加剧，箱体顶部出现破口。李营等[9]针对边缘拉伸失效的舱壁破坏形式进行研究，并发现设置加强筋和边缘变形协调装置能够较好地限制爆炸破口。以上对于舱室内爆的研究均以承受单发弹药攻击为背景，而实际上舱室可能承受多发弹药的攻击。
根据相关研究[10-11]多次爆炸对于舱壁的毁伤规律有明显的影响。因此有必要对舱室承受多发弹药攻击时舱室结构的毁伤规律进行研究。本文主要通过数值仿真和量纲分析的方法，对两次内爆加载下结构的毁伤规律进行研究。
2 有限元模型及计算方法
2.1 有限元模型
本文设计的舱室大小为2.4 m×2.4 m×2.4 m，为对称结构，本章研究药量、爆距、起爆时间间隔对舱室的破坏效应，有限元模型取1/4全尺寸模型进行计算，以节省计算资源。有限元模型如图1所示，共包含3种空气、Q235钢、TNT 3种物质，舱壁边界处采用固支边界条件，空气域边界处设置流出边界，整个模型的网格数量约为140万。
图1 有限元模型示意图
Fig.1 Finite element model
在反舰武器中，配备半穿甲战斗部的弹药主要有导弹和炮弹，炮弹的攻击距离近、火力密集，更有可能造成多次爆炸。目前，国内外发展比较成熟的舰炮口径主要在76～130 mm，各种炮弹的等效TNT当量约为0.5～10 kg[12]，因此本文的药量设置为2～10 kg。取药量为2、4、6、8、10 kg，起爆点距离舱壁中心为0.4、0.6、0.8、1、1.2 m的球形装药，2次爆炸均取上述各值，共记125种工况，进行数值仿真计算，仿真方法在2.1节中有叙述。
2.2 材料模型选择
在有限元仿真中，选取合适的本构方程和准确的材料参数，才可能得出正确的仿真结果。根据文献[6]本文选取的材料模型如下
空气采用Ideal Gas状态方程，其表达式如下：
P=(γ-1)ρe
(1)
式中：P为压强；γ为绝热指数；ρ为密度；e为比内能，具体参数如表1所示。
表1 空气状态方程参数
Table 1 Parameters of the air state equation
ρ/(kg·m-3)γe/J1.2251.42.068×105
舱壁钢材采用Q235钢，根据AUTODYN材料库中的STEEL 4340钢修改而得到，本构模型采用Johnson-Cook方程，适用于爆炸载荷下结构的变形，其表达式如下式。
(2)
式中：A为屈服静应力；B为应变硬化指数；C为应变率敏感系数；n为应变硬化指数；m为温度相关指数；
为有效塑性应变；
为等效塑性应变；T*=(T-T0)/(T1-T0)为相对温度，T0为参考温度，T1为材料熔点，T为某一时刻的温度，具体参数如表2。
表2 Q235钢材料参数
Table 2 Q235 steel material parameters
A/MPaB/MPaCnmT1/K3705100.260.0141.031 793
TNT炸药采用JWL状态方程描述，如下式。
(3)
式中：ν=ρ0/ρ为爆轰产物的相对比容，ρ0为炸药初始密度，ρ为炸药爆轰后某一时刻的密度。E为单位体积爆轰产物的内能；A，B，R1，R2，ω为相关参数，具体参数如表3。
表3 TNT材料参数
Table 3 TNT material parameters
ρ0/(g·cm-3)E/(kJ·g-1)A/kPaB/kPa1.634.83.73×1083.74×106R1R2ω4.150.90.35
3 数值仿真验证
3.1 数值仿真方法
本文通过2次数值仿真的方法进行2次内爆的加载。为了验证AUTODYN模拟多次内爆的可行性，将内爆实验与数值仿真结果进行对比。姚术健[4]对多次内部爆炸作用下箱型结构的累积损伤进行了实验，其结果可用于验证AUTODYN计算多次爆炸下累计损伤的可行性。本文选取SB-Ⅲ-1进行验证。该箱型结构边长为600 mm，边界板长度为120 mm，壁厚4 mm，第1次炸药质量为98.44 g，顶部开有圆形小孔用于放置炸药。若将箱型结构的形心视为原点，忽略顶部小孔，则该结构关于3个坐标平面是对称的，为了节省计算资源，可选择1/8模型进行计算，空气域模型为边长为450 mm的正方体，炸药为正方体形状，起爆点为坐标原点。对于流固耦合中的接触来说，流体网格不大于固体网格，对上述空气域模型分别采用5、6、7、8 mm的网格大小进行计算固体网格为流体网格的2倍。
空气初始内能设置为2.068×105 J，在空气域的最外层网格设置Flow_Out流出边界。空气采用欧拉网格，炸药填充在欧拉域中，舱壁选择拉格朗日网格，箱壁交界处采用刚性连接，欧拉和拉格朗日选择全自动接触，耦合类型选择全耦合，耦合方式为手动耦合，厚度为2倍欧拉网格尺寸，箱壁中心处和舱室内部设置高斯测点。
第一次爆炸计算时间取为20 ms，根据图2可知，计算进行到2 ms左右时，舱壁中心挠度已经不再大幅度增加，而是震荡收敛与某一固定值。为保证爆炸计算数据的连续性，第2次爆炸应该在第1次爆炸的结果文件的基础上进行，首先将结果文件中的欧拉域重新用空气材料填充，之后导入炸药的映射文件或者在欧拉域中重新填充炸药，上述操作会导致系统内部的能量不守恒，这是由于外部第二次重新填充炸药导致的系统不孤立造成的，并非计算出错导致，在AUTODYN 控制面板中能量检查选项中将循环次数设置为一个较大的值即可进行后续计算。
3.2 结构变形
首先对比不同网格尺寸模型的箱壁中心位移时程曲线，验证了该模型满足收敛性原则，之后使用误差最小的网格尺寸进行第2次内爆。2次内爆的结构变形程度相似，箱壁中心挠度误差较小，说明使用AUTODYN进行多次内爆计算时可行。
使用1/8模型进行计算时，仿真结果的每个箱壁中心位移曲线相同，取其中一面的位移曲线与实验值进行对比，为方便分析取前10 ms的曲线，图2表示不同网格尺寸下某一个侧面箱壁中心高斯测点位移值与实验值。从图2可以看出，大约2.5 ms时刻之后箱壁挠度值已经不再增加，而是震荡收敛与某一固定值。随着网格的减小，箱壁中心挠度的极大值和极小值均随着网格增大而增大。即收敛值逐渐增大，网格大小为5 mm时网格震荡幅度最小，其收敛趋势最靠近实验值，网格大小为6 mm和5 mm时的挠度曲线变化规律比较相似，震荡幅度相差也较小。网格大小为的7 mm挠度曲线震荡幅度明显大于5 mm和6 mm时的幅度，网格为8 mm的挠度曲线震荡幅度最大，距离实验值相差最远，这表明随着网格的减小，仿真结果是收敛的。
图2 不同网格尺寸的箱壁中心挠度时程曲线
Fig.2 Time-history curve of center deflection of box wall with different mesh sizes
根据图2，挠度值先大幅度增长，而后在一定区间内不断震荡，选取趋于稳定时的挠度数据，统计其平均值，用平均值来代表仿真的收敛值，各次仿真结果及其误差列于表4，其中“+”代表仿真值大于实验值，“-”代表仿真值小于实验值。从表4可知，随着网格尺寸的逐渐减小，误差也逐渐减小，但4种网格尺寸的计算误差均在8%以内，仿真结果的可靠性较高。
表4 不同网格尺寸的挠度
Table 4 Deflection/mm for different mesh sizes
网格仿真值实验值误差534.2332.93+3.9634.5532.93+4.9735.2932.93+7.2835.5332.93+7.9
图3表示两次爆炸的箱壁中心位移曲线，其中20 ms之后的位移曲线代表第2次爆炸的位移。可以看出，第2次爆炸箱壁震荡的幅度小于第1次，选取曲线挠度处于震荡收敛阶段的位移平均值作为仿真挠度值，如图所示，第2次爆炸挠度值在第1次的基础上有较小幅度的增加。
图3 两次内爆箱壁中心挠度时程曲线
Fig.3 Time-history curves of the center deflection of the chamber wall for two implosions
对仿真结果进行整理，将2次爆炸的挠度误差列于表5，可以得出第1次爆炸误差为3.9%，第2次爆炸误差为1.8%，2次仿真与实验有很高的相似度。
表5 两次内爆实验与仿真挠度
Table 5 Two implosion experiments and simulated deflections
爆炸次数TNT质量/g仿真值/mm实验值/mm误差198.4434.2332.93+3.9296.1435.7236.36-1.8
上述的箱型结构内部爆炸的实验与仿真整体变形结果，使用AUTODYN计算舱内爆炸和多次爆炸的累计毁伤的误差较小，仿真的可靠性和精度较高，说明使用AUTODYN进行舱室结构的累积毁伤是可行的。
4 内爆量纲分析与毁伤模式分布
舱内爆炸是一个高度非线性的过程，目前，针对舱内爆炸的理论分析主要从量纲分析的角度进行研究。此处的量纲分析在献中[13]的基础上进行。
内爆载荷作用下船体板架有挠曲变形、边界撕裂、局部破口和整体剪切等不同的破坏模式，根据计算结果从两方面来衡量舱壁的破坏情况：舱壁最大挠度，用字母δ表示；舱壁最大破口面积，用字母S表示。由于多次内爆作用下舱室结构的破坏模式问题比较复杂，影响δ、S的因素众多，本文主要考虑单位质量TNT的爆热Qv、爆炸TNT质量m，炸药距离舱壁的最短距离R，舱壁材料屈服强度σ0、密度ρs、声速Cs、厚度H，正方形舱室的边长L，初始时刻空气的密度ρ0、声速C0。
根据上述问题将记因变量为{δ，S}，记自变量为Qv，m，R，σ0，ρs，Cs，L，H，ρ0，C0，函数定性关系式记为，
{δ，S}=(Qv，m，R，σ0，ρs，Cs，L，H，ρ0，C0)
(4)
选取Qv，σ0，H为参考物理量，使用文献[14]中基于线性代数的量纲分析方法，得到的无量纲物理量为9个，包括2个无量纲因变量和7个无量纲自变量。利用得到的无量纲量可得到无量纲函数关系式如下，
{π1，π2}=ψ(π4，π5，π7，π8，π9，π11，π12)
(5)
(6)
对于确定种类炸药(本文中为TNT)在稳定的空气中爆炸，(Qv，ρ0，C0)=const，若认为舱壁材料的密度、声速和屈服极限在爆炸响应过程中保持不变，(ρ0，Cs，σ0)=const，可进一步简化为
(7)
上述函数关系表明
与
这些自变量相关，将各个自变量以指数形式耦合后形成的新自变量
即内爆时舱壁的最大毁伤程度可以用以Dm为自变量的函数表示。上述无量纲数Dm中包含炸药的种类Qv，炸药质量m，舱壁材料特性σs，舱壁厚度δ，正方形舱室的边长L，炸药和舱壁之间的最短距离R。值得注意的是，此处进行的量纲分析是针对内爆问题进行的研究，与爆炸次数无关，下文中的第1次内爆和第2次内爆均使用相同的表达式，区别在于两次爆炸式表达式中m和R的数值不同。
对仿真工况进行分析，发现2次爆炸的药量为2～10 kg时，单舱室结构的破坏模式有2种，即挠曲变形(记为模式一)和局部破口(记为模式二)。对于2次爆炸来讲，舱壁的最终变形结果与第1次和第2次爆炸都密切相关。
受文献[2]中的启发，为了初步研究2次爆炸对舱壁的毁伤程度。分别计算第1次爆炸和第2次爆炸的Dm值，分别将其记为Dm1和Dm2。从图4可以看出，左下部分区域为模式1，右上部分区域为模式2，将模式1与模式2有较为明显的分界线，使用多项式拟合可以较好的显示出两者的分界线。
图4 2次爆炸的无量纲损伤数关系
Fig.4 Dimensionless damage number map of two explosions
模式一与模式二的分界线的二次多项式关系如下：
(10)
利用破坏模式分布图7和式(10)，可以通过两次爆炸的无量纲损伤数来判断单舱室结构的毁伤模式。
5 结论
1)使用AUTODYN进行多次内爆加载的数值仿真方法是可行的。
2)2次爆炸时，随着爆炸当量和起爆位置的变化，舱壁出现挠曲变形和局部破口2种不同的毁伤模式。
3)分析2次爆炸的无量纲损伤数(Dm1，Dm2)关系可知，2种毁伤模式之间有较为明显的分界线，该分界线可用多项式函数拟合。
[1] 李旭东.内爆准静态压力载荷对舱壁结构的毁伤效应研究[D].太原：中北大学，2020.
Li X D.Research on the damage effect of implosion quasi-static pressure load on bulkhead structure[D].Taiyuan:North University of China，2020.
[2] Fast prediction method of failure modes for steelbox structures under internal blast loading[J].Engineering Failure Analysis,2021(120):104919.
[3] 孙博华.量纲分析以及应用[J].物理与工程，2016，26(06):11-20.
Sun B H.Dimensional Analysis and Application[J].Physics and Engineering，2016，26(06):11-20.
[4] 姚术健.箱形结构内部爆炸等效缩比实验方法及破坏特性研究[D].长沙：国防科学技术大学，2016.
Yao S J.Experimental method and failure characteristics of internal explosion of box-shaped structures[D].Changsha:National University of Defense Technology，2016.
[5] 郑成，孔祥韶，周沪，等.全封闭舱内爆炸载荷作用下薄板变形研究[J].兵工学报.2018，39(8):1582-1589.
Zheng C，Kong X S，Zhou H，et al.Research on the deformation of the thin plate under the explosion load in the fully enclosed cabin[J].Journal of Ordnance Engineering.2018，39(8):1582-1589.
[6] 焦晓龙，赵鹏铎，姚养无，等.基于仿真与量纲分析的不同药量TNT内爆下多舱室结构毁伤规律研究[J].爆炸与冲击，2020，40(08):124-133.
Jiao X L，Zhao P D，Yao Y W，et al.Research on the damage law of multi-compartment structures under different doses of TNT implosion based on simulation and dimensional analysis[J].Explosion and Shock，2020，40(08):124-133.
[7] 徐维铮，吴卫国.基于伯努利方程舱室内爆炸准静态泄压过程研究[J].应用力学学报，2019，36(06):1504-1508，1530.
Xu W Z，Wu W G.Research on quasi-static pressure relief process of cabin explosion based on Bernoulli equation[J].Chinese Journal of Applied Mechanics，2019，36(06):1504-1508，1530.
[8] 汪维，刘光昆，汪琴，等.钢结构箱体内爆作用下毁伤破坏试验研究[J].北京理工大学学报，2019，39(12):1225-1231.
Wang W，Liu G K，Wang Q，et al.Experimental study on damage and damage of steel structure box under internal explosion[J].Journal of Beijing Institute of Technology，2019，39(12):1225-1231.
[9] 李营，张磊，杜志鹏，等.舱内爆炸作用下舰船舱壁失效机理与抗破损设计[J].中国造船，2019，60(03):27-34.
Li Y，Zhang L，Du Z P，et al.Failure mechanism and anti-damage design of ship bulkhead under the action of explosion in cabin[J].China Shipbuilding，2019，60(03):27-34.
[10] 李旭东，尹建平，杜志鹏，等.多次水下爆炸钢制圆板应变与挠度增长规律分析[J].振动与冲击，2020，39(05):131-136.
Li X D，Yin J P，Du Z P，et al.Analysis of strain and deflection growth law of steel circular plate made of multiple underwater explosion[J].Vibration and Shock，2020，39(05):131-136.
[11] 张斐.水下爆炸作用下焊接钢板变形及破损规律研究[D].太原：中北大学，2019.
Zhang F.Study on deformation and damage law of welded steel plates under underwater explosion[D].Taiyuan:North University of China，2019.
[12] 杜志鹏，张磊，赵鹏铎.实尺度舰船舱内炮弹静爆试验方法[J].中国舰船研究，2018，13(03):103-109.
Du Z P，Zhang L，Zhao P D.Test method for static explosion of shells in the cabin of a full-scale ship[J].Chinese Ship Research，2018，13(03):103-109.
[13] 姚术健，张舵，郑监，等.内部爆炸作用下钢箱结构变形规律性实验[J].爆炸与冲击，2017，37(05):964-968.
Yao S J，Zhang D，Zheng J，et al.Experiment on deformation regularity of steel box structure under internal explosion[J].Explosion and Shock，2017，37(05):964-968.
[14] 钟巍，田宙，寿列枫.基于线性代数的大规模快速量纲分析算法及其在爆炸与冲击工程研究中的应用[J].计算数学，2020，42(02):170-195.
Zhong W，Tian Z，Shou L F.Large-scale fast dimensional analysis algorithm based on linear algebra and its application in explosion and impact engineering research[J].Computational Mathematics，2020，42(02):170-195.