1 引 言

高战略价值的重要军事目标的地下化，促进了侵彻武器与钻地武器的发展，而对于舰船、机库、防御工事等带腔室结构进行有效毁伤，内爆炸是最为有效的打击手段^[1]。与开放环境相比，炸药在密闭及半密闭的有限空间内爆炸过程具有冲击波效应显著增强、爆轰产物膨胀做功能力成为重要影响因素、温度效应明显等特征^[2]。一般将有限空间内爆炸的作用过程分为首个冲击波和准静态压力两部分，准静态压力由于其作用时间长、冲量大等特点，成为有限空间内爆炸研究的一个重要考虑因素，产生的效应倍受关注。因此，对有限空间内爆炸导致结构变形、失效的规律进行深入研究，建立合理有效的有限空间内爆炸威力评估方法，对弹药毁伤和防护结构设计具有指导性作用。目前主要通过小尺寸、小当量缩比实验对有限空间内爆炸进行研究，实验常常在爆炸容器、爆炸塔^[3,4,5,6]等设备中进行。爆炸容器作为一类特殊密封压力容器，通常是为限制一定量的爆炸物爆炸带来的破坏和污染等负面效应而设计，可以对容器内部的爆炸过程进行有效约束，将冲击波、爆炸破片等限制在其内部。目前常见的爆炸容器的承载主体普遍为球形或圆柱形，其壳体在受内部径向冲击载荷作用下的动态响应过程一直是爆炸容器研究领域的关注点之一。爆炸容器按照设计准则可分为单次使用与多次使用两种类型。其中，单次使用型爆炸容器允许材料发生一定的塑性变形，但不发生破坏；多次使用型爆炸容器的防护当量必须限制于材料弹性极限内，不能产生塑性变形。

1960年Baker^[7]曾针对薄球壳在内部冲击载荷作用下的动态冲击响应过程，提出了单自由度（Single degree of freedom, SDoF）运动方程，获得薄球壳受冲击载荷后的位移响应时程曲线的解析解。SDoF方程为研究球壳和圆柱壳在一维应变状态下的受内部冲击载荷动态响应过程提供了一个简明有效的解析方法，并出现了很多具有代表性的应用，包括Ko^[8]提出的多层式弹塑性金属薄球壳模型和赵士达^[9]提出的动力系数设计方法等，对金属材料爆炸容器的研究、设计和改进有着极为重要的指导性影响。相关实验和理论研究^[3,4,5,6,10]中发现，炸药在爆炸容器内部密闭空间内爆炸后，在首个冲击波后会出现一个比冲击波超压峰值小但作用时间较长、变化缓慢的准静态压力。之前普遍认为准静态压力因作用时间长、冲量大，对结构的毁伤效应更为显著，但还缺乏力学分析和明确认识^[11,12,13]。因此，针对爆炸容器在考虑准静态压力在有限空间内爆炸下的动态响应过程进行研究，对于爆炸效应评估和防护结构设计显得尤为重要。本文主要针对多次使用的球形爆炸容器，进行内爆炸准静压效应的弹性结构响应分析工作。

2 力学分析模型

2.1 研究模型

美国洛斯阿拉莫斯实验室获得的典型球形爆炸容器内爆载荷曲线^[11]如图1所示。前人研究中^{[7,8,9,10,11,12]}，常常将爆炸容器受内部爆炸冲击的载荷曲线简化为如图2a中的三角脉冲形式，忽略准静态压力带来的影响。因此，结合图1的爆炸载荷曲线，本研究将考虑准静态压力的爆炸载荷曲线简化为图2b^[14]，并建立简化载荷的表达式如式(1)。式中，p_m1为首个脉冲压力峰值，t₁为首个个脉冲作用时间，p_m2为准静态压力。当不考虑准静态压力时，令p_m2=0即可。

图1 球形爆炸容器内爆载荷曲线^[11]

Fig.1 Internal blast loading in an spherical explosion containment vessels^[11]

a. Loading without quasi‑static pressure

b. Loading with quasi‑static pressure

图2 简化内爆冲击载荷

Fig.2 Simplified internal blast loading

在承受来自容器中心的爆炸冲击载荷时，建立如图3所示的球壳单元来分析动态响应过程。此情况下，由于冲击载荷沿径向分布作用于球壳内表面，球壳的动态响应过程简化为在径向位移上的单一自由度受迫振动问题，因此可利用单自由度运动模型对动态响应过程进行推导分析，得出位移随时间变化的解析解。对于本研究中所用到的球壳，r为球壳中线半径，41 mm；h为壳体厚度，2 mm。球壳材料为弹性材料，E为弹性模量，200 GPa；v为泊松比，0.3；ρ为密度，7830 kg·m^-3。此外，图3中σ_θ和σ_φ为经向和纬向应力，u_r为径向位移。

图3 球壳单元

Fig.3 Element of the spherical shell

2.2 响应过程分析

对于球壳受内部冲击载荷作用情况，对如图3所示的球壳微元，在dθ方向上，故σ_θ在径向上的分量可表示为：

同理可得σ_φr≈σ_φdφ，根据牛顿第二定律

由于球壳表面有σ_θ=σ_φ=σ，且根据胡克定律有

故可知，球壳受内部冲击载荷时，其单自由度运动方程^[7]可表示为式(5)

式中，ω为振动频率，振动周期T=2π/ω。

将考虑准静态压力的爆炸载荷简化表达式(1)带入方程后可以求解得到位移响应的解，如式(6)

令p_m2=0时可得到不考虑准静态压力的特解，和Baker提出的薄球壳在内爆载荷下的弹性响应解^[7]是一致的。

对式(6)中u_r求导得到式(7)，当du_r/dt=0时，可得达到最大位移的时间与最大位移值。

当满足t=t₁时du_r/dt=0，即最大位移出现时刻恰好在首个脉冲冲量结束，则有式(8)

因此首个个脉冲作用时间t₁需满足

在赵士达^[9]的工作中将其约等为t_c=0.375T，本研究中获得的值更为精确。

当t₁＜t_c时，首个脉冲作用时间相对较短，最大位移出现于准静态压力作用期间；当t₁＝t_c时，最大位移出现时刻恰好为首个脉冲冲量结束时刻；当t₁＞t_c时，最大位移会随准静态压力p_m2大小的改变，可能出现在准静态压力作用期间，也可能出现在首个脉冲作用期间，需要进一步比较准静态压力阶段位移最大值与首个脉冲段位移最大值从而判断最大位移的发生阶段。

通过式(10)可计算得到一个p_m2临界值，当满足t₁＞t_c与p_m2＝p_m2c时，后续准静态压力段的位移最大值与首个脉冲作用期间的位移最大值恰好相等。

满足t₁＞t_c且0 MPa≤p_m2＜p_m2c时，后续准静态压力段的位移最大值与小于首个脉冲作用期间的位移最大值，因此最大位移发生于首个脉冲作用期间，不随准静态压力的增大而变化；满足t₁＞t_c且p_m2＞p_m2c时，后续准静态压力段的位移最大值与大于首个脉冲作用期间的位移最大值，最大位移发生于准静态压力段，受准静态压力p_m2的影响，随着准静态压力的增大而增大。

由式(6)可知，进入准静态压力作用段后的后续响应过程可表示为正弦函数形式，因此后文中所提及的振幅值，均指t>t₁后准静态压力阶段的弹性等幅度振动振幅值，计算方式如式(11)

3 数值模拟校验

使用LS‑DYNA^[15]对球壳内部冲击载荷下的弹性动态响应过程进行数值模拟，并与理论解析结果进行对比。建立一个如图4所示的1/8球壳模型，球壳中线半径41 mm，厚度2 mm，进行有限元网格划分时使用六面体SOLID单元，沿经线方向与赤道线方向上网格数量150个，径向网格数量6个。在使用LS‑DYNA求解器进行计算时，材料属性设置为MAT_ELASTIC，弹性模量200 GPa，泊松比0.3，密度7830 kg·m^-3。冲击载荷按照图2b所示的简化载荷曲线及式(1)进行定义，径向均布加载于球壳内表面。计算时长150 μs，计算步长0.1 μs。

图4 1/8球壳有限元模型

Fig.4 Finite element model of one‑eighth spherical shell

如表1所示，假定多种不同的考虑准静态压力的载荷Load 1~Load 3，带入式(6)进行计算，并将公式计算得到的响应的相关特征值如最大位移、振幅和振动周期，与LS‑DYNA数值模拟结果进行对比验证。载荷Load 1~Load 3的对比结果如图5和表2所示，图5中u_r为径向位移，mm；表2中u_{r max}为径向最大位移，mm；Amplitude为后续阶段振幅，mm；T为振动周期，μs；ε_r为理论解与数值模拟结果的相对误差。可以发现理论解析公式与LS‑DYNA数值模拟结果符合情况很好，公式的准确性得到了数值模拟的验证。

图5 位移响应理论解与LS‑DYNA数值模拟结果对比

Fig.5 Comparison between analytical solution and LS‑DYNA simulation

4 准静态压力的影响

基于前期研究经验^[16]，图6~图7、图9,10,11的算例中，使首个脉冲冲量大小一致。图中u_{r max}表示整个响应阶段最大位移值，mm；amplitude表示后续弹性等幅振动阶段的振幅，mm；u_{r min}表示后续弹性振动阶段最小位移值即反向最大位移，mm。针对本研究中的球壳算例，t_c=11 μs。

a.

b.

图9 位移响应曲线和不同准静压下的关键计算参数（p_m1=364 MPa，t₁=11 μs）

Fig.9 Dynamic response and calculated key parameters with different quasi‑static pressure (p_m1=364 MPa，t₁=11 μs)

a. Dynamic response

b. Parameters related to p_m2

图10 位移响应曲线和不同准静压下的关键计算参数（p_m1=200 MPa，t₁=20 μs）

Fig.10 Dynamic response and calculated key parameters with different quasi‑static pressure (p_m1=200 MPa，t₁=20 μs)

a. Dynamic response

b. Parameters related to p_m2

图11 位移响应曲线和不同准静压下的关键计算参数（p_m1=100 MPa，t₁=40 μs）

Fig.11 Dynamic response and calculated key parameters with different quasi‑static pressure (p_m1=100 MPa，t₁=40 μs)

4.1 t₁<t_c的情况

在保证首个脉冲冲量不变的情况下，针对t₁<11 μs的情况，即最大位移必然发生于准静态压力作用期间情况，进行了如图6,7中p_m1=800 MPa，t₁=5 μs和p_m1=500 MPa，t₁=8 μs等情况在不同准静态压力p_m2水平下的分析。图6a和图7a分别为两种情况在p_m2=0,30,60 MPa三种准静态压力下计算得到的位移响应曲线，图6b和图7b则是将p_m2作为变量，得到最大位移、最小位移和振幅值关于p_m2变化的关系曲线。

a. Dynamic response

b. Parameters related to p_m2

图6 位移响应曲线和不同准静压下的关键计算参数（p_m1=800 MPa，t₁=5 μs）

Fig.6 Dynamic response and calculated key parameters with different quasi‑static pressure (p_m1=800 MPa，t₁=5 μs)

a. Dynamic response

b. Parameters related to p_m2

图7 位移响应曲线和不同准静压下的关键计算参数（p_m1=500 MPa，t₁=8 μs）

Fig.7 Dynamic response and calculated key parameters with different quasi‑static pressure (p_m1=500 MPa，t₁=8 μs)

因频率仅受结构尺寸与材料属性的影响，在不同准静态压力水平下，最大位移和最小位移出现的时间有一定的差异。在式(6)的计算中也发现因p_m1、t₁与p_m2的不同，导致C₁、C₂的区别，从而在计算中最大位移出现时刻存在差异；此外，最大位移出现在准静态压力段，受到准静态压力的影响，随着准静态压力水平的增大而增大，当首个脉冲冲量不变、t₁接近临界值11 μs时，由于最大位移受到准静态压力作用时间的减少，增长趋势放缓；最小位移随准静态压力增大而增大，受到准静态压力的影响相较于最大位移更为明显；振幅随准静态压力增大呈减小趋势。

在图8首个脉冲冲量不等情况的算例计算中，p_m1=200 MPa，t₁=5 μs。对比图6和图8中情况可以发现，当首个脉冲峰值减小，而作用时间不变时，首个脉冲冲量降低，p_m2对最大位移、最小位移和振幅的影响程度更为明显。

a. Dynamic response

b. Parameters related to p_m2

图8 位移响应曲线和不同准静压下的关键计算参数（p_m1=200 MPa，t₁=5 μs）

Fig.8 Dynamic response and calculated key parameters with different quasi‑static pressure (p_m1=200 MPa，t₁=5 μs)

4.2 t₁＝t_c的情况

在保证首个脉冲冲量不变的情况下，针对t₁=11 μs，即最大位移出现时刻恰好为首个脉冲冲量结束时刻情况，进行如图9中p_m1=364 MPa，t₁=11 μs情况在不同准静态压力p_m2水平下的分析。图9a为该情况在p_m2=0,30,60 MPa三种准静态压力下计算得到的位移响应曲线，图9b为最大位移、最小位移和振幅值关于p_m2变化的关系曲线。结果显示在t₁＝t_c情况下，准静态压力对最小位移和振幅产生影响，最小位移随准静态压力增大而增大，振幅随准静态压力增大而减小，但最大位移值不会发生改变。

4.3 t₁＞t_c的情况

如图10和图11，针对t₁＞11 μs，即最大位移会随着准静态压力p_m2大小的改变，可能出现在准静态压力作用期间，也可能出现在首个脉冲作用期间的情况，进行了p_m1=200 MPa，t₁=20 μs及p_m1=100 MPa，t₁=40 μs两种状态在不同准静态压力p_m2水平下的分析。类似地，图10a和图11a分别为两种情况在p_m2=0,30,60 MPa三种准静态压力下计算得到的位移响应曲线，图10b和图11b为最大位移、最小位移和振幅值关于准静态压力p_m2变化的关系曲线。

最大位移发生于首个脉冲冲量作用期间时，因不受准静态压力的影响，最大位移随着准静态压力的增大而基本保持不变，最小位移受到准静态压力的影响更为明显。而随着p_m2的增大，准静态压力作用期间的弹性响应过程中会发生更大的最大位移，直至大于首个脉冲段的最大位移，此时表现为最大位移出现在准静态压力阶段，因此受到p_m2的影响，最大位移再次随着p_m2的增大而增大，振幅值随准静态压力增大而减小。

通过式(10)，可以求解出最大位移发生于准静态压力段的p_m2临界值，例如在p_m1=200 MPa，t₁=20 μs情况中计算可得p_m2c＝26 MPa，当p_m2>p_m2c后，最大位移出现在准静态压力期间并随p_m2的增大而增大，而p_m2＜p_m2c时，最大位移出现在首个脉冲作用期间，不随p_m2的改变而改变。这种变化在图10b中表现为最大位移随着p_m2的增大而出现了一个拐点。计算中还发现，对于p_m1=100 MPa，t₁=40 μs情况，p_m2需要超过86 MPa，才能使最大位移出现在准静态压力期间。

准静压p_m2临界值的计算结果受到振动频率、首个脉冲压力峰值、首个个脉冲作用时间、壳体厚度、材料杨氏模量和泊松比等多个因素的影响。此外，即使图6,7和图9,10,11中的首个冲击波冲量相同，但结构响应由于首个冲击波超压、脉宽不同而出现差异，特别是t₁＞t_c时首峰幅值相差较大。因此，在研究有限空间内爆炸压力对结构的效应时，还应该关注结构自身参数，比如密度、厚度、杨氏模量、固有频率等。采用多次使用型容器开展有限空间内爆炸效应威力评估时，不仅要考虑爆炸载荷的首个冲击波及后续准静态压力的具体细节，还要与结构动态响应历程及其自身参数的分析紧密结合。

5 能量分析方法

利用LS‑DYNA数值模拟输出响应过程中系统总能量变化曲线，并与位移响应曲线和载荷曲线进行比对，从能量分析方法的角度对准静压的影响进行机理分析与解释。

5.1 t₁<t_c的情况

当t₁＜t_c时，最大位移发生于准静态压力作用期间，以图12所示p_m1=800 MPa，t₁=5 μs，p_m2=60 MPa情况为例进行分析。

图12 系统总能量受冲量影响(p_m1=800 MPa，t₁=5 μs，p_m2=60 MPa)

Fig.12 Total energy affected by impulse(p_m1=800 MPa，t₁=5 μs，p_m2=60 MPa)

在位移响应首次达到最大值后，由于首个脉冲作用时间相对较短，此时处于准静态压力段。在阶段0中，首个脉冲载荷结束，进入准静态压力段，准静态压力与壳体运动方向相同，对系统做正功W₀，因此准静态压力p_m2的提高会导致W₀的增大，系统获得更多的能量，从而令最大位移值提高。

在阶段1位移由最大值运动至最小值，准静态压力与壳体运动方向相反，对系统做负功W₁，在阶段2位移由最小值运动至最大值，准静态压力与壳体运动方向相同，对系统做正功W₂。由于准静态压力p_m2在两个阶段中相同，因此W₁＝W₂，正负功互相抵消，系统最大总能量不变，因此后续响应过程中的位移最大值与首个位移最大值保持相等。准静态压力p_m2水平的提高，会导致位移由最大值运动至最小值的过程中负功W₁的增大，因此系统最小总能量减小，位移最小值减小，振幅减小。

5.2 t₁=t_c的情况

当t₁＝t_c时，最大位移出现时刻恰好为首个脉冲冲量结束时刻，以图13所示p_m1=364 MPa，t₁=11 μs情况为例进行分析。

图13 系统总能量受冲量影响(p_m1=364 MPa，t₁=11 μs，p_m2=60 MPa)

Fig.13 Total energy affected by impulse(p_m1=364 MPa，t₁=11 μs，p_m2=60 MPa)

与5.1中类似，同样存在位移由最大值运动至最小值的阶段1与位移由最小值运动至最大值的阶段2。由于准静态压力p_m2在两个阶段中相同，因此W₁＝W₂，正负功互相抵消，系统最大总能量不变，因此后续响应过程中的位移最大值与首个位移最大值保持相等。而准静态压力p_m2水平的提高同样会导致位移由最大值运动至最小值的过程中负功W₁的增大，系统最小总能量减小，位移最小值减小，振幅减小。

与t₁＜t_c情况的不同之处在于首个位移最大值出现时刻恰好首个脉冲冲量结束，对应W₀＝0，因而准静态压力p_m2的提高并不会导致最大位移值提高。

5.3 t₁>t_c的情况

当t₁＞t_c时，最大位移会随着准静态压力p_m2大小的改变，可能出现在准静态压力作用期间，也可能出现在首个脉冲作用期间，以p_m1=200 MPa，t₁=20 μs情况为例进行分析，通过式(10)计算得p_m2c＝26 MPa。因此如图14，分析p_m2=0,26,60 MPa三种情况的LS‑DYNA数值模拟得到的系统总能量变化曲线，可以发现在三种工况下，都可以分为三个阶段进行功和能量的分析，0阶段首个脉冲与壳体运动方向相反，对系统做负功W₀，在1阶段准静态压力与壳体运动方向同样相反，做负功W₁，而在2阶段，准静态压力与壳体运动方向相同，做正功W₂。

a. p_m2=0 MPa

b. p_m2=26 MPa

c. p_m2=60 MPa

图14 系统总能量受冲量影响(p_m1=200 MPa， t₁=20 μs)

Fig.14 Total energy affected by impulse (p_m1=200 MPa， t₁=20 μs)

(1) 当p_m2<26 MPa时，如图14a情况中p_m2=0 MPa的情况下，(W₀+W₁)>W₂，负功使得进入准静态压力段后系统最大总能量小于首个脉冲段，因此准静态压力段的最大位移小于首个脉冲段的最大位移，即最大位移出现在首个脉冲作用期间；

(2) 随着p_m2的增大，W_Ⅰ与W_Ⅱ也相应的增大，直到当p_m2=26 MPa时，即图14b情况下，(W₀+W₁)=W₂，负功与正功相互抵消，此时系统最大总能量不变，因此准静态压力段的最大位移与首个脉冲段的最大位移相等；

(3) p_m2>26 MPa时，如图14c中的p_m2=60 MPa情况，此时由于(W₀+W₁)<W₂，准静态压力段的正功超过了首个脉冲与部分准静态压力段的负功总和，增大了系统最大总能量，因此准静态压力段的最大位移大于首个脉冲段的最大位移，即最大位移出现在准静态压力段。

此外，图14b和图14c两种情况同时也表明在后续阶段即使准静态压力作用时间很长，输入冲量很大，但在振动响应过程中由于正负功相互抵消，实际上系统总能量并没有一直增加，而是同样维持着一个正弦变化的过程。

6 结论

（1）首个脉冲作用时间临界值t_c与结构固有频率有关，准静态压力临界值p_m2c受首个脉冲压力峰值，首个个脉冲作用时间、结构频率、壳体厚度、材料弹性模量和泊松比等多个因素的影响。首个脉冲作用时间和准静态压力与其临界值的大小关系对最大位移出现时刻具有决定性作用。

（2）首个脉冲冲量作用时间t₁<t_c时，最大位移出现于准静态压力期间；t₁＝t_c时最大位移出现时刻恰好为首个脉冲冲量结束时刻；t₁＞t_c时存在一个准静态压力临界值p_m2c，若p_m2＜p_m2c则最大位移出现在首个脉冲作用期间，若p_m2＞p_m2c则最大位移出现在准静态压力作用期间。最大位移出现在准静态压力作用期间时随准静态压力增大而增大，出现在首个脉冲作用期间时则不受准静态压力影响。

（3）最小位移随准静态压力增大而增大，后续等幅振动幅值随准静态压力增大而减小。

本研究为球形容器在考虑准静态压力的内爆炸载荷作用下弹性动态响应研究提供了有效的分析方法，可对爆炸效应评估、弹药设计和防护结构设计提供指导，后续还将针对弹塑性模型和更多复杂材料与结构开展深入研究。

参考文献