饱和冲量及其等效方法在舱室内爆炸中的应用

2021-04-27朱凌田岚仁李德聪王佳月张霄

中国舰船研究 2021年2期

朱凌，田岚仁，李德聪，王佳月，张霄

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

2 中国舰船研究设计中心，湖北武汉 430064

0 引言

船舶与海洋工程结构物在运行过程中会受到脉冲载荷的冲击作用，例如爆炸、砰击等。在结构响应研究和实际结构设计中，通常将爆炸载荷简化为具有初速度的冲击载荷或者脉冲载荷。而随着各种反舰武器性能的提高，其能够穿透船体外板在舱室内部发生爆炸，使得爆炸载荷对船体结构的威胁日益增加，因此研究约束空间（有限空间）内的抗爆、抗冲击影响越来越重要。

然而，相比敞开环境下的爆炸，约束空间的内爆炸载荷远比自由场的空爆载荷复杂，原因在于有限空间内的爆炸冲击波的壁面反射会造成多峰效应，经反射波叠加，其爆炸峰值会高于敞开环境下的峰值，作用时间也显著变长，且通常会存在准静态超压的现象[1-3]。迄今，针对舱室内爆炸冲击波的多峰效应及其存在的准静态超压等现象的研究成果颇丰，但尚未完全探明内爆炸冲击载荷作用下结构响应的力学机理。

20 世纪50 年代以来，许多学者研究了脉冲载荷作用下的结构动态塑性响应特性，所采用的大部分理论方法也都是基于理想的刚塑性（R-PP）假设，这是因为在忽略弹性变形时可显著简化理论模型，从而得到塑性变形的解析表达式。目前，在大挠度下的结构动态响应刚塑性理论求解方面应用得较广泛的是模态近似分析法[4]，运用该方法避免了考虑瞬态响应阶段的移行铰模型在求解中塑性铰两侧物理量不连续性和非线性的问题，使理论分析变得更简单。之后，国外许多学者又进一步发展了模态近似分析法[5-12]。在国内，黄震球[13]和颜丰等[14]基于动量和动量矩守恒发展了加筋板结构塑性变形的解析方法，所得结果与试验结果吻合较好。吴有生等[15]和牟金磊等[16]采用能量法考虑大变形时的应变关系及中面膜力的影响，求解了船体板结构塑性大变形。

上述研究爆炸冲击载荷作用下的结构响应特性针对的脉冲载荷持续时间较短，即结构变形将在载荷结束之前结束。然而，如果载荷持续时间较长，结构动力响应中将出现饱和冲量现象，从而影响采用理论方法计算所得结果的合理性和准确性。因此，本文将系统地总结饱和冲量概念的提出及研究发展成果，并以舱室内爆炸为典型算例，分析内爆炸载荷曲线的特性和结构响应特征，基于饱和等效方法将复杂的内爆炸载荷等效为矩形脉冲载荷，采用理论及数值方法计算等效载荷，给出相关的设计曲线和计算公式。通过研究舱室内爆炸作用下的饱和冲量现象及饱和等效方法，以更准确地得到结构塑性动力响应结果，并预测结构变形的终止时间与载荷间的关系。运用上述方法还可减少频繁的复杂非线性数值计算，更有效地开展船体结构的抗冲击设计优化。

1 饱和冲量概念的提出与研究发展

20 世纪90 年代，Zhao 等[17]定义了刚塑性结构在中等强度脉冲载荷作用下的一种特征现象—“饱和冲量”，其反映的是板在塑性动力响应中的一个特性，即板受到强烈的横向压力脉冲载荷作用时会产生较大的变形，而大变形诱发的膜力对板的承载能力起到了增强作用。若板受到足够长时间的矩形压力脉冲作用，就只有脉冲载荷达到饱和时间前的冲量对最大及最终变形有影响，之后的加载脉冲不会进一步增加变形量。此后，Zhao 等[18]又将饱和冲量概念拓展到基于理想的刚塑性模型的简支圆板、简支和固支方板以及圆柱壳中，利用模态分析法得到了封闭形式的解析解。然而，“饱和冲量”概念只对应于理想的刚塑性结构最大变形情况，未考虑弹性影响。鉴此，Zhu[19]提出一种运用有限差分方法来分析固支方板的弹塑性响应数值程序，基于该程序，Zhu等[20]提出了分别对应于“最大挠度”和“最终挠度”的“饱和冲量”，从而完善了基于刚塑性和弹塑性模型的饱和冲量现象的定义。

针对基于刚塑性−弹塑性饱和冲量现象的问题，席丰等[21]利用最小加速度原理，通过建立数值方法，分析了脉冲载荷作用下简支圆板的动力响应，指出在高载荷范围内也存在脉冲载荷作用下的“饱和冲量”现象。此后，席丰等[22]又分析了脉冲载荷作用下的钢梁动力响应及反常行为的应变率效应，指出在脉冲载荷作用下发生塑性变形的钢梁总是存在“饱和冲量”现象，且发生时与载荷强度及载荷作用时间相关。

近年来，武汉理工大学朱凌教授的研究团队针对饱和冲量问题开展了进一步研究，并更系统地予以了分析，例如：尺度效应对方板饱和冲量的影响[23]；长宽比和边界条件对矩形板饱和冲量的影响[24]；材料应变率敏感性和应变硬化对饱和冲量的影响[25]；不同脉冲载荷作用下的方板饱和冲量[26-27]；考虑移行铰的方板饱和冲量[27]；同时考虑瞬态响应阶段和准确屈服面的梁饱和冲量[28]。此外，还基于饱和冲量的研究，提出了物理意义更明确且计算结果更准确的脉冲载荷等效方法[26-27]，改进了Youngdahl[29]于70 年代提出的半经验等效方法。然而，上述研究尚未在饱和冲量的应用层面开展更深入的探讨。

2 约束空间爆炸下船体板饱和冲量计算

2.1 有限元模型

结构响应的理论分析通常是在假定载荷已知的基础上开展的，不同于敞开环境下的爆炸，约束空间内爆炸载荷曲线通常很难被写成某种较为通用的函数形式。因此，可以针对需要分析的模型，利用有限元软件计算出爆炸所产生的脉冲载荷。此外，过去几十年以来，学者们开展了大量的爆炸载荷实验和数值研究，若针对的是某一系列结构，可通过前人的研究成果确定较为合理的载荷曲线形式，并直接对结构进行动力响应的理论分析。本文采用AUTODYN 有限元软件计算舱室内爆炸载荷及结构动态响应，选取图1 所示舰船右舷典型舱室的结构计算模型[1]，模型的尺寸为5 000 mm×3 000 mm×2 500 mm。

图1 简化后的舱室结构计算模型Fig. 1 Simplified calculation model of cabin structure

为简化研究和定量分析，仅采用无加筋板进行原理性探讨，并参考文献[30-31]将舱室内爆炸简化为炸药在舱室中心爆炸。如图2 所示，可变形壁面（横舱壁）为3 000 mm×2 500 mm×14 mm（2L（长）×2B（宽）×H（厚））的矩形板，其他5 个壁面为刚性壁面，设定舱室为完全封闭空间。矩形板的4 个端部为固支，空气域网格划分为110×60×50，网格数330 000，采用高阶单物质Euler-FCT 求解器。矩形板网格划分为120×100，网格数12 000，板与空气采用完全耦合实现相互作用。本文仿真计算选取4 组TNT 球状炸药。

图2 舱室内爆炸有限元模型Fig. 2 Finite element model of cabin explosion

2.2 计算材料模型

根据文献[1]中的材料参数，选取Q235 钢作为矩形板材料，密度为7 800 kg/m3，杨氏模量E=210 GPa，泊松比υ=0.3。采用双线性弹塑性本构模型，材料的应变率效应由Cowper-Symonds模型描述，相关参数见表1。

自由空气采用的是γ 律状态方程描述。TNT炸药爆轰产物的压力则由JWL 状态方程描述:

式中：P为爆轰产物的压力；V=ρ0/ρ，为相对体积，其中ρ0=1 630 kg/m3，为炸药初始密度，ρ 为炸药爆轰产物的密度；E0为单位体积炸药的初始内能；A1，A2，R1，R2，w为与炸药性质相关的材料常数；e 为自然常数。其他参数见表2。

表1 Cowper-Symonds 强度模型参数Table 1 Parameters of Cowper-Symonds strength model

表2 TNT 炸药JWL 状态方程参数Table 2 Parameters of TNT explosives JWL state equation

2.3 仿真计算结果

图3 所示为采用20 kg 炸药在内爆炸作用下横舱壁的冲击波压力时历曲线。由图可见：壁面反射作用使冲击波压力含有多个峰值；因舱室模型未设置泄爆口，内爆炸后产生的准静态超压Ps始终作用在结构上，板上不同位置的压力峰值有所差别，但准静态压力值基本一致。由于受到复杂的反射波影响，约束空间内爆炸下的最大压力峰值Pm尚无可广泛使用的经验公式。对于舱室内爆炸准静态超压，文献[2]总结了4 种经验公式。图4 为本文有限元计算得到的准静态超压值（带圆形黑色实线）与这4 种经验公式的对比。图中红色叉点在横坐标上的时间即对应于炸药量下的饱和时间tsat。由图可见：本文计算结果与Moir和Carlson 经验公式得到的结果较为接近。同时，还可发现最大压力峰值（带矩形黑色实线）和准静态超压随着炸药量的增加都大致呈线性上升的趋势。

图3 横舱壁冲击波压力时历曲线（20 kg TNT）Fig. 3 Pressure time histories of blast wave on the transverse bulkhead (20 kg TNT)

图5 所示为3 组不同炸药量下的爆炸载荷冲量I随时间的变化曲线。由图可见，在反射波作用下，冲量经过初始时间段的波动后大致呈线性上升的趋势，这主要是后期准静态超压所致。图6所示为3 组不同炸药量下板中心挠度W0随时间的变化曲线。结合图3 可发现，板在饱和时间tsat（蓝色点划线所示）后仍受到大于板的静态塑性极限载荷的压力作用，冲量也始终处于上升的趋势，但板的变形已达到最大值且未再增加。这说明只有早期压力脉冲载荷所形成的冲量对板的实际变形有效，后期的冲击压力并未引起板的变形值进一步增加，此即为“饱和冲量”现象，对应时间tsat为饱和时间，对应冲量Isat为饱和冲量。

图4 压力峰值和准静态超压随炸药量的变化曲线Fig. 4 Variation of pressure peak values and quasi-static overpressure with explosive mass

图5 不同炸药量下冲量时历曲线Fig. 5 Time histories of impulse with different explosive mass

图6 不同炸药量下板中心挠度时历曲线Fig. 6 Deflection time histories of plate center with different explosive mass

上述分析表明，舱室内爆炸下结构存在“饱和冲量”现象，若预测由此导致的结构大变形应基于“饱和冲量”的相关方法。但是，在弹性效应的影响下，由永久挠度来确定tsat的精确值较难，故可选择最大挠度进行分析。需要指出的是，基于最大挠度分析的tsat值通常偏小，但随着载荷幅值的增大，即结构塑性变形的增加，2 种挠度确定的tsat的差值会逐渐缩小。

图7 和图8 分别给出了饱和时间tsat、饱和冲量Isat和板中心点最大挠度W0sat随炸药量m的变化关系。由图可知：饱和时间随炸药量的增加逐渐减少，但减少的趋势逐渐变缓，这与刚塑性的理论分析中饱和时间不随压力峰值变化有所出入，主要是由弹性效应及载荷曲线发生变化所致。此外，最大挠度和饱和冲量则随炸药量的增加而大致呈线性上升的趋势。尽管本文计算算例中的载荷曲线较复杂，但两者的趋势与文献中采用的矩形脉冲[24]、线性衰减脉冲[26]或者线性上升指数衰减脉冲[27]的图线趋势一致。

图7 饱和时间随炸药量的变化曲线Fig. 7 Variation of saturated time with explosive mass

图8 饱和冲量与板的中心最大挠度随炸药量的变化曲线Fig. 8 Variation of saturated impulse and maximum deflection of plate center with explosive mass

3 饱和等效方法的应用与验证

3.1 载荷等效处理过程

在工程实际中，脉冲形状通常十分复杂，例如舱室内爆炸的压力时历曲线具有多个峰值且形状不规则，这对于求解脉冲载荷作用下的结构响应和结构设计造成了极大困难。而工程中有很多结构需优化设计和反复计算，若每次都采用非线性有限元计算，计算资源消耗太大。因此，基于有限元计算或实验测量等方法得到载荷时历曲线，可将其引入到载荷等效方法中，结合理论分析快速计算结构响应，以减少繁琐的复杂非线性数值计算。文献[29] 提出的消除脉冲载荷形状影响的等效方法，可将一个任意形状的脉冲由有效冲量Ie和有效载荷Pe这2 个量来表征：

式中：P(t)为脉冲载荷；ty和tf分别为塑性变形开始和结束时刻；2tmean为等效后的矩形脉冲长度。通过式（3）找出该有效脉冲的重心，作为等效矩形脉冲的中心，再由式（4）确定等效矩形脉冲的幅值，从而对原脉冲载荷曲线进行“掐头去尾”后得到有效冲量Ie。

在实际应用中，通常很难确定tf，文献[29]提出采用经验公式（5）进行“去尾”处理。

式中，Py为结构的塑性极限载荷。

在过去几十年内，Youngdahl 这种半经验等效方法被广泛用来简化处理实际的工程问题。然而，该等效方法并未给出理论依据。近几年，文献[24-25]在饱和冲量研究的基础上，提出了物理意义更明确，且计算结果更准确的解析脉冲载荷等效方法，从而改进了Youngdahl 半经验等效方法。改进后的等效方法由饱和时间tsat代替式（5）中的tf，划分出脉冲中的有效部分，得到饱和冲量Isat，然后利用文献[29]提出的“重心等效”方法获得等效的矩形脉冲载荷。

在本文算例中，采用饱和等效方法对实际舱室的内爆炸载荷进行了等效，原理如图9 所示。图中，对于固支矩形板，静态塑性极限载荷由式（6）计算[32]，

将塑性极限弯矩M0=σ0H2/4和宽长比β=B/L代入，可以计算得到本文模型的静态塑性极限载荷为Py=74.6 kPa。

图9 饱和等效方法示意图Fig. 9 Schematic diagram of saturated equivalent method

由图9 可见，载荷曲线初始段上升得很快。为简化分析，ty可以取为0，计算结果的误差将随着载荷幅值（炸药量）的增大而减小。Youngdahl 等效方法建议采用式（5）近似地确定塑性变形结束时刻tf，但从本文算例可见，在存在准静态超压的情况下，整条载荷曲线的值长时间处于静态塑性极限载荷Py之上，使得式（5）无法成立。这说明Youngdahl 等效方法不适用于舱室内爆炸载荷的等效计算。表3 给出了等效后的矩形脉冲载荷参数。需要注意的是，基于等效后的载荷采用理论方法或者有限元法计算响应时，应考虑是否有大气压的影响。表3 中，P0为实际压力峰值，Isat'为减去大气压影响后的饱和冲量。

表3 等效后的矩形脉冲载荷参数Table 3 Parameters of equivalent rectangular pulse loading

3.2 等效方法的理论与仿真验证

3.2.1 刚塑性理论方法

得到矩形等效载荷后，可非常方便地对结构响应进行理论分析及求解。白雪玉等[24]针对矩形脉冲载荷作用下矩形板的饱和冲量现象进行刚塑性理论分析，给出了板的无量纲饱和挠度及饱和时间的计算公式：

式（7）和式（8）适用于2tmean≥tsat的情况，但等效后的矩形脉冲长度2tmean

然而，在爆炸载荷作用下，由于船体钢结构响应通常受应变率的影响很大，所以为得到更准确的结构大变形的解，分析时计及应变率效应具有工程实际意义。文献[33]假设应变率服从Cowper-Symonds 模型，从对刚体撞击下的矩形板不同区域的应变率进行计算，给出了矩形板的整体应变率加强因子。为量化应变率效应对饱和冲量现象的影响，文献[25]针对脉冲加载的方板进一步开展了研究。这些方法也可扩展到脉冲加载的矩形板研究中。

对于2tmean≥tsat的情况，假设应变率服从Cowper-Symonds 模型，空间平均应变率为

对于2tmean

则考虑了应变率效应的无量纲饱和挠度为

因此，可得等效后的矩形脉冲载荷作用下船体结构的最大挠度理论结果。由于文献[20]前期研究中指出饱和时间与板厚无关，所以可将等效载荷代入到上述理论公式中进行结构的板厚设计，从而可减少复杂的非线性数值计算。

3.2.2 弹塑性有限元仿真

鉴于AUTODYN 软件中不便于在板结构上准确施加均布载荷，本文采用Abaqus/Explicit 软件对板结构施加等效后的矩形脉冲载荷进行仿真计算。计算时，采用与AUTODYN 中相同的模型，但无需设置空气域。矩形板网格取60×50，网格类型为S4R 四节点减缩积分壳单元，如图10 所示。矩形板材料设置为AUTODYN 软件计算中相同的材料参数，计及应变率效应。