泄爆薄板厚度对舱室破坏及舱内载荷影响的数值仿真

赵鹏铎，黄松*，，尹建平，张磊，黄阳洋，，徐豫新

1海军研究院，北京100161

2中北大学机电工程学院，太原030051

3北京理工大学爆炸科学与技术国家重点实验室，北京100081

0 引 言

目前，各种高性能的反舰武器成为水面舰船面临的主要威胁［1］。反舰武器战斗部进入舱内爆炸后形成的冲击波超压和准静态压力，都会导致舱室结构发生严重破坏，而采用泄爆技术可将舱内爆炸冲击波和爆轰产物能量通过泄爆口引出船体外，使舱内冲击波在还未形成积聚效应前（十几到几十毫秒）就将其压力降低，从而起到保护舰内指挥室和弹药库等重要舱室的作用。传统的弹药库通常布置在不易受敌方武器攻击的舰船中心位置，而DDG 1000驱逐舰出于隐身性考虑，将其设计成干舷内倾式船型，导致主甲板变窄而无法在舰船中心布置导弹垂直发射系统。为此，美国海军专门研发了沿船体周边布置的外围导弹垂直发射系统，如图1所示。该垂直发射系统弹药舱的外侧舱壁较为薄弱，因此对内侧舱壁进行了加强［2-5］。当导弹在发射单元内意外发生爆炸时，外侧舱壁将瞬间飞出，通过这种方式，可达到舱室泄爆的效果，这一效果己在美国海军实验中得到验证（图2）［5］。DDG 1000布置于舷侧的垂直发射系统颠覆了传统弹药库要隐藏在最不易受攻击的核心位置的设计原则，突破了泄爆防护关键技术，为驱逐舰抗反舰导弹爆炸攻击探索出了一种新的途径。

图1 DDG 1000驱逐舰外围垂直发射系统位置图Fig.1 Peripheral VLS cells located in DDG 1000 destroyer

图2 薄板泄爆模型爆炸实验及结果Fig.2 The explosion test and its results of the thin plate explosion venting model

从现有的文献来看，国内对于舱室泄爆结构的研究还主要集中在预开泄爆孔大小、布置方式对舱内冲击波特性及舱室破坏的影响等方面［6-9］，而对于DDG 1000这种新型舷侧薄板泄爆结构的研究还鲜有报道。本文将参考DDG 1000采用的舷侧薄板泄爆结构，通过AUTODYN有限元软件建立有限元模型进行数值仿真，对不同泄爆薄板厚度下舱室的破坏情况进行对比分析，研究内爆后舱内压力、比冲量的变化情况，用以为我国舰船舷侧舱室等结构开展泄爆设计提供一定的参考。

1 数值计算方法验证

1.1 数值仿真模型的建立

Edri等［10］针对一系列长方体装药下的舱内爆炸进行了实验（图3），并对不同装药下的载荷特性进行了对比研究。本文取该实验中装药量为4 kg TNT的舱室内爆工况，采用数值仿真手段模拟该工况下的内爆实验，并通过实验数据与仿真数据的对比来验证仿真方法的可靠性。实验所用舱室内部空间尺寸为290 cm×290 cm×270 cm，舱壁材料由2块钢板之间的增强高强度混凝土浇铸而成，整体厚度为35 cm（实验中忽略舱壁的动态变形，可视为刚体）。实验用TNT炸药为长方形，尺寸为208 mm×134 mm×104 mm（L1×L2×L3），位于舱室中心，其中L1与L3确定的平面平行于安装传感器的舱壁面。起爆点位于炸药上表面的几何中心。在实验中，共布置了9处压力传感器，其中，在测点RP1，RP7和RP8处未采集到数据，本文仿真将取实验中采集到了数据的6个位置（RP2，RP3，RP4，RP5，RP6和RP9）进行压力验证。

图3 文献［10］中的实验模型、测点布置及炸药示意图Fig.3 Test model，measuring points and TNT of Reference［10］

舱室有限元模型直接在AUTODYN软件中建立。考虑到模型具有一定的对称性，为简化计算并节省时间，采用1/2对称模型进行仿真计算。舱室有限元模型由空气域、舱壁板和炸药组成，炸药位于舱室中心，如图4所示。由于不考虑舱壁的动态响应，舱壁材料选用4340钢，并在舱壁上设置刚性边界条件使其为刚性壁面，以模拟舱室模型的主体结构。

图4 仿真模型图Fig.4 Simulation model

1.2 仿真材料参数的选择

仿真中，空气采用理想气体状态方程：

式中：p为流体压力；γ为绝热指数；ρa为空气密度；e为比内能。本文中，取γ=1.4，ρ=1.225×10-3g·cm-3，e=2.068×105kJ/kg。

炸药采用TNT装药，爆轰产物用JWL状态方程描述。该状态方程本质上是一个经验公式，其数据来源于一系列的物理实验，能很好地反映产物的体积、压力和能量特性，适用于大多数高能炸药。爆轰产物的气体压强由下式给出：

式中：a，b，R1，R2，ω为实验拟合参数；E0为单位体积爆轰产物的内能；v=ρ0/ρt，为爆轰产物的相对比容，其中ρ0为炸药初始密度，ρt为爆轰后某一时刻的密度。TNT参数采用AUTODYN软件材料库中的默认参数［11］。

舱壁材料采用4340钢，将用Linear状态方程和Johnson-Cook强度模型来描述。Linear状态方程假设压力与内能无关，材料密度变化小，变化过程是可逆的，通常用于固体。

式中：K为材料的体积模量；μ为压缩比，μ=(ρ/ρ0)-1。

Johnson-Cook强度模型通常用来描述大应变、大应变率和高温下的应力应变关系。

式中：A，B，C，n，m均为材料常数，其中A为屈服应力，B为硬化常数，C为应变率常数，n为硬化指数，m为相应的温度指数；εp为有效塑性应变；εp*为有效塑性应变率；Tmelt为材料融化温度；Troom为室温。AUTODYN软件中对应的参数如表1和表2所示［11］。

表1 4340钢Linear状态方程参数Table 1 Parameters of Linear state equation of 4340 steel

表2 4340钢Johnson-Cook强度模型参数Table 2 Parameters of Johnson-Cook strength model of 4340 steel

1.3 仿真与实验结果对比

从计算结果中提取6个测点的压力时程曲线（图5）与实验中相对应测点的实验结果进行对比。文献［10］没有给出各测点的压力时程曲线，仅给出了反射超压峰值和准静态压力值。本文仿真结果与实验结果［10］的对比情况如表3所示。

由图5可以看到，同位于两壁面角隅位置的测点PR2，PR6及PR9的反射超压值相差较大。这是因为偏移起爆点会影响舱内压力的分布，根据文献［12］的研究，发现起爆点位置对结构内壁面各测点的反射冲击波超压峰值有显著影响。炸药中心起爆后，舱内等压线分布均匀，而当炸药偏中心起爆时，舱内等压线分布则较紊乱，偏向起爆位置一侧的等压线分布更为密集。此外，舱内压力会在不同角度的壁面测点处发生复杂的反射波叠加作用而导致压力大小有差异。

图5 仿真测点压力时程曲线Fig.5 Timehistorycurvesofmeasuringpointpressureinsimulation

表3 仿真结果与实验结果[10]对比Table 3 Comparison of simulation results and experimental results[10]

从表3中可以看到，除测点PR2处的误差较大外（27.3%），其余测点的反射超压及准静态压力的误差值均在允许范围内。表3中所示准静态压力值取6个测点在25 ms时刻测得压力的平均值（各测点的准静态压力值取波动基线处的压力值）。可见该数值计算方法具有一定的可靠性，可为下一步泄爆舱室模型的仿真计算提供保证。

2 泄爆舱室模型建立及结果分析

2.1 泄爆舱室仿真模型建立

通过简化DDG 1000舷侧垂直发射系统弹药舱的结构，参考其整体尺寸建立了1∶1三维舱室模型，如图6所示。其中，舱室长5 m，宽3 m，高5 m。舷侧舱壁外板的中间部分为泄爆薄板结构，居中设计，其尺寸大小为2.5 m×2.5 m。舱室四周（上、下甲板及前、后横舱壁）为刚性舱壁，不考虑其动态响应。两条纵向舱壁为正常舱壁，其中左侧纵向舱壁（靠近舷侧）具有泄爆薄板结构，右侧纵向舱壁为防护舱壁，舱壁材料选用4340钢。除泄爆薄板厚度外，其余舱壁的等效厚度均为10 mm。建立的泄爆舱1/2对称有限元模型如图7所示，其中欧拉空气域边界设置流出边界条件，以避免冲击波在空气域边界发生反射而造成计算误差。炸药为135 kg TNT方形裸炸药，位于舱室中心。4340钢、TNT炸药及空气的材料参数见1.2节。图7中，“1”为在防护舱壁中心处设置的冲击波压力测点，“2”为在舱壁中心处设置的挠度变化测点。

图6 舱室的三维模型Fig.6 3D model of cabin

图7 1/2有限元模型Fig.7 1/2 finite element model

2.2 典型薄板泄爆舱室泄爆过程分析

在前期的仿真中，发现左、右舱壁厚度均为10 mm时，在舱室中心处放置的135 kg方形TNT炸药起爆后防护舱壁刚好发生破损（图9中的3#舱室）。为便于进行对比分析，仿真中所使用的TNT药量也为135 kg。这里取泄爆薄板厚度为4 mm，对135 kg TNT炸药在舱室中心的爆炸进行分析。不同时刻舱内的压力变化云图如图8所示。

从图8中可以看到，炸药在舱内起爆后，冲击波最先到达离炸点最近的左、右舱壁面。0.4 ms时，初始冲击波到达左、右舱壁并和舱壁发生相互作用，作用范围逐渐扩大。0.8 ms时，垂直方向的初始冲击波与上、下舱壁及左、右舱壁的相互作用导致在舱壁附近出现了高压区。1.2 ms时，舱室壁面反射的冲击波与向角隅结构方向传播的冲击波在角隅处发生汇聚，此时，在压力载荷的作用下，左侧泄爆薄板开始发生塑性变形，右侧防护舱壁变形不明显。3.2 ms时，初始冲击波在传播至上、下舱壁与左、右舱壁相交的位置处时形成冲击波汇聚，与紧随其后而来的上、下舱壁的反射冲击波叠加，在舱室中部位置形成了汇聚高压区域，由于右侧舱壁的反射冲击波强于左侧，导致高压区逐渐向左舱壁一侧移动（如4.0 ms时）［13］。4.8 ms时，泄爆薄板中心处在压力载荷的作用下发生失效，出现破口，舱内压力开始外泄，右侧防护舱壁仅发生了弯曲变形。8.0 ms时，左侧舱壁反射的冲击波到达右侧舱壁，在右侧舱壁附近形成高压区。随着时间的推移，在压力泄出的过程中，左侧舱壁的破坏沿着泄爆薄板与正常舱壁的连接处及边角逐渐加大，破损的舱壁形成破片飞出。由于破孔较大，压力泄出导致舱内压力逐渐降低。左侧舱壁完全破坏时的舱室的破损情况如图8（h）所示，可见右侧的防护舱壁仅发生了塑性变形，未发生破损。

图8 不同时刻舱内压力云图Fig.8 The pressure contours at different times in cabin

2.3 泄爆薄板厚度对舱室破坏方式的影响

本节对泄爆薄板厚度分别为2，6，10 mm的3种舱室进行仿真分析，以研究薄板厚度对舱室破坏方式的影响。为便于后文的表述，这里将2，6和10 mm薄板舱室分别定义为1#舱室、2#舱室和3#舱室。炸药均位于舱室中心处，3种情况下舱室典型时刻的破坏情况如图9所示。

由图9可以看到，在40 ms时刻，1#舱室和2#舱室左侧舱壁均出现了大破口，且为薄板部分整体飞出（1#舱室在8.4和12 ms时，由于薄板飞出的距离较远，故未给出已飞出的薄板图），右侧防护舱壁仅发生了塑性变形，未发生破坏，而3#舱室的左、右舱壁均发生了破口。通过对比1#舱室和2#舱室可以看到，1#舱室在1.6 ms时在薄板与舱壁连接边的中心处最先发生失效。这是因为当板厚较薄时，在载荷的作用下发生了较大的塑性变形，连接处因受到拉伸作用而发生了拉伸失效。2#舱室在4.8 ms时在薄板与舱壁的连接边角处最先发生失效。这是因为当板厚较厚时，在相同的载荷作用下塑性变形困难，而在连接边角处发生了应力集中，其因受到较大的作用力而发生失效。随后，这2种舱室的破坏情况都将沿着连接边界逐渐扩大，最后在冲击波压力的作用下，整块薄板飞出舱体，形成泄爆口。最后，破损的左舱壁在压力的作用下向外翻出（例如，1#舱室和2#舱室的12和40 ms时刻）。

通过对以上3种舱室破坏情况的比较可以发现，1#，2#舱室的右侧防护舱壁仅发生了塑性变形，而3#舱室的右侧防护舱壁则发生了破口，这说明泄爆薄板结构可以有效防护右侧舱壁。舱壁的破坏失效最先发生在薄板与舱壁的连接处，并逐渐向舱壁边角处扩大。因此，薄板厚度越小，在薄板和舱壁连接处就越容易发生失效，也就可以更早地形成泄爆口，其形成表现为薄板整体飞出舱室。

图9 不同薄板厚度下不同时刻舱室典型破坏情况Fig.9 Typical damage conditons of cabin with different thin plate thickness at different times

2.4 泄爆薄板厚度对舱内载荷的影响

图10和图11分别为内爆后3种舱室内测点位置的压力变化及比冲量变化时程曲线。

由图10可以看到，在7.5 ms之前，3种舱室在测点处的压力曲线完全重合。1#，2#舱室在7.5 ms之前左侧舱壁均发生了破损，薄板飞出；对比3#舱室可以发现，泄爆结构对舱内前两次的压力峰值几乎没有影响。从图10中也可以看到7.5 ms之前，测点处的比冲量也基本相同。

图10 测点位置压力变化时程曲线Fig.10 Pressure change curves of measuring point location

图11 测点位置比冲量变化时程曲线Fig.11 Specific impulse time curves of determination location

在 7.5～11.5 ms时刻，2#，3#舱室在测点处的压力曲线重合度较高，其中2#舱室测点处的压力略低于3#舱室，1#舱室在测点处的压力最低。这是因为1#舱室在2.4 ms时薄板就发生失效而形成了泄爆口，在7.5 ms后舱内大部分压力已泄出，无法形成压力很高的第3次反射波。而2#舱室在7.5 ms时刚形成泄爆口不久，舱内压力还未来得及大量泄出，故此时间段内测点处的压力与3#舱室相比压力略低。比冲量为压力对时间的积分，在该段时间，3种舱室的比冲量开始发生变化，如图11所示。

由图11可看到，在11 ms之后，由于3#舱室还未出现破口，舱内冲击波在多次复杂的反射和相互作用下，逐渐衰减为作用时间很长的准静态压力，其值约为2.11×103kPa（取波动基线处的压力值）。而1#舱室和2#舱室由于泄爆的作用，舱内压力逐渐降低，最终没有形成准静态压力。从图11所示的比冲量时程图中还可以看到，在11 ms之后，1#，2#舱室的比冲量变化逐渐趋于平缓，最终1#舱室维持在3.53×104Pa·s，2#舱室维持在 4.14×104Pa·s，降比约14.7%。由于3#舱室内存在作用时间很长的准静态压力，导致舱内的比冲量随时间继续上升，在40 ms时达到 8.95×104Pa·s，为 1#，2#舱室的 2倍以上（由于发生了破口，随着时间的推移，其值还会继续增大，但增长速率有所降低）。可见，3#舱室防护舱壁失效破坏发生的主要原因是在比冲量的作用下舱壁中心处变形达到拉伸极限。

通过对舱内载荷变化的分析发现，薄板泄爆结构对舱内初始冲击波超压几乎没有影响，但可以有效降低舱内的准静态压力及比冲量，且薄板越薄，泄压效果越明显。在本文所研究药量和结构下，造成防护舱壁破坏的主要因素不是初始冲击波超压，而是长时间存在的准静态压力。

2.5 泄爆薄板厚度对防护舱壁变形挠度的影响

图12所示为3种舱室防护舱壁在内爆炸载荷下的挠度变化时程曲线。由图可以看出，在7.5 ms之前，3种舱室的防护舱壁变形挠度变化趋势基本相同；在7.5 ms之后的准静态压力作用阶段，防护舱壁的挠度不再增大，变化波动较小。40 ms时，3种舱室防护舱壁的变形挠度依次为871.6，910.8，937.2 mm，最大和最小挠度相差约65 mm。可见在该药量和结构下，造成防护舱壁达到最大变形挠度的主要因素是前期的初始冲击波和反射冲击波，后期的准静态压力对该挠度增大的影响较小。图12中，3#舱室后期的挠度曲线为直线，是因为防护舱壁发生了破口而未能测出后续的挠度变化波动。

图12 测点位置挠度变化时程曲线Fig.12 Deflection curves of measuring point location

3 比冲量及挠度变化函数模型拟合

舱内发生爆炸后，为能快速预估不同薄板厚度下作用在防护舱壁上的比冲量和挠度变化大小，本节通过数据拟合，得到了计算比冲量和挠度的函数模型。除上述已有的3组仿真数据外，又增加了泄爆薄板厚度为4和8 mm的仿真数据。不同板厚下获得的比冲量和挠度大小如表4所示。由于板厚为10 mm的舱室在40 ms时仅发生小破口，比冲量还在继续增大，故拟合时舍去了其值。拟合结果如图13和图14所示，得到的函数模型及相关系数（R2）如表5所示。由表5可以看到，两者的函数模型均为二次函数，相关系数达到了0.99以上。该数值可以反映出数据的离散程度，越接近±1，表示数据的相关度越高，拟合越好。

表4 不同薄板厚度下防护舱壁处比冲量和挠度值Table 4 Specific impulse and deflection of the protective bulkhead at different thickness

图13 比冲量随薄板厚度变化曲线Fig.13 Variation of specific impulse with respect to thin plate thickness

图14 挠度随薄板厚度变化曲线Fig.14 Variation of deflection with respect to thin plate thickness

表5 拟合的函数模型及相关系数Table 5 Fitting function model and correlation coefficients

4 结 论

本文采用数值分析的方法模拟了内爆情况下新型薄板泄爆舱室的毁伤过程，分析了不同泄爆薄板厚度下舱室的破坏情况及舱内载荷的变化情况，得到以下主要结论：

1）舱内发生爆炸后，新型薄板泄爆结构能有效形成泄爆口。舱壁的破坏失效最先发生在薄板和舱壁连接处，并逐渐向舱壁边角处扩大。薄板厚度越小，在内爆载荷作用下，薄板和舱壁连接处越容易发生失效，从而也就越容易形成泄爆口，且泄爆口的形成表现为薄板整体飞出舱室。

2）造成防护舱壁达到最大变形挠度的主要因素是前期的初始冲击波和反射冲击波，而造成防护舱壁最终破坏的主要因素是长时间作用的准静态压力。泄爆结构对舱内的初始冲击波超压几乎没有影响，但可以明显降低舱内的准静态压力和比冲量。其中，具有泄爆薄板的舱室在一定时间后舱内准静态压力将降为0；比冲量随着薄板厚度的减小而降低。

3）通过数据拟合得到比冲量和挠度随泄爆薄板厚度变化的二次函数模型，其相关系数达0.99以上，可以很好地描述比冲量和挠度随泄爆薄板厚度的变化规律。