爆炸下球壳变形空间周期分布的理论计算方法*

刘文祥，张德志，钟方平，程 帅，张庆明

（1. 西北核技术研究院强动载与效应实验室，陕西 西安 710024；2. 北京理工大学爆炸科学与技术国家重点实验室，北京 100081）

1976 年Buzukov[1]实验发现了爆炸容器弹性变形范围内的应变增长现象，在这种现象中容器壳体的最大变形没有出现在第一个振动周期内，而是出现在后续周期内。应变增长现象对容器安全来说是不利的，因为传统的容器设计方法一般以第一个周期内的变形峰值为依据。近几年，应变增长现象的研究有不少新发现。Dong 等[2]、Li 等[3]发现壳体不稳定性导致弯曲振动进而形成应变增长现象，并研究了相关机理；刘文祥等[4-5]实验观察到壳体塑性变形下的应变增长现象，获得了应变增长系数超过6 的数据，甚至推测存在系数超过10 的情况。

振动叠加是形成应变增长现象的重要原因[2,6-7]。Duffey 等[7]基于振动力学中“拍”形成机理，结合壳体振动频率分析结果以及出现应变增长现象的应变曲线的特征，推测频率相近的振动叠加形成应变增长现象。Liu 等[8]提出了对应变增长现象解剖式分析的方法，针对带扰动源球壳上的应变增长现象，分离出弯曲应变（弯曲波引起的壳体弯曲变形导致的应变）和膜应变（壳体拉伸变形导致的应变），直接证明了频率相近的膜振动和弯曲振动叠加是导致应变增长现象的主要原因，还发现了球壳变形呈空间周期分布。

解剖式分析是研究振动叠加形成应变增长现象的新方法，早期研究虽提出了壳体弯曲波、变形空间周期性分布的相关机理，但缺乏理论的支撑。本文中参考Timoshenko 梁的弯曲理论[9-10]，理论分析球壳上弯曲波的波速以及壳体变形空间分布的周期，并与数值仿真结果进行对比，以验证理论计算的可靠性。

1 球壳上的弯曲波和变形空间周期分布

球形爆炸容器往往存在人孔门、观察窗等功能部件，在爆炸作用下这些部件与球壳的运动存在差异，由扰动源处产生弯曲波并在球壳上传播，这些部件被称为扰动源。完全约束扰动源是指绝对静止状态的扰动源，其导致的弯曲波比其他约束条件扰动源的更强，导致的应变增长现象也更严重。Liu 等[8]利用数值仿真分析了内径523 mm、厚度3 mm 的完全约束扰动源下球壳的响应情况。壳体模型见图1，其中扰动源半径为L，应变输出位置与旋转对称轴的距离为r′，相应的夹角为α，应变输出位置与扰动源正对位置的弧线长度为x。

图2 为L=10 mm 时球壳上扰动源正对位置外壁和内壁的应变曲线，可见该位置出现了明显的应变增长现象。Liu 等[8]的研究表明，由于扰动源处产生的弯曲波在扰动源正对位置处汇聚，导致此位置出现了球壳上最严重的应变增长现象。

Liu 等[8]根据薄壳变形机理（见图3），提出了由壳体外壁应变和内壁应变计算得到弯曲应变和膜应变的方法，其中弯曲应变为弯曲波引起壳体弯曲振动导致的应变，膜应变是壳体膜振动引起的壳体切向拉压变形导致的应变。

图1 球壳数值模型Fig. 1 The numerical model for the spherical shell

图2 扰动源半径为10 mm 时球壳上扰动源正对位置的应变曲线Fig. 2 Total strain-time curves of spherical shellat the pole opposite to the site of perturbation when the disturbance source radius is 10 mm

图3 薄壳变形机理Fig. 3 Deformation mechanism of thin shell

图4 为分离出来的L=10 mm 时壳体上不同位置的弯曲应变曲线，可以观察到弯曲波由扰动源产生并在球壳上的传播。其中弯曲波A 为频率与球壳膜振动频率相近的第一个弯曲波，弯曲波波群的头部为波长最短的弯曲波。由于弯曲波频率越高，波速越快，波长短的波将逐渐超过波长长的波。根据弯曲波到达的壳体位置以及相应时刻，可计算出弯曲波A 以及最短弯曲波的传播速度，如表1所示，其中弯曲波A 的传播速度为420～450 m/s，最短弯曲波的传播速度为2 900～3 200 m/s。

图4 壳体上不同位置的弯曲应变曲线和不同时刻的壳体上弯曲波Fig. 4 Bending strain-time curves at different positions on spherical shell and bend waves at different moments

表1 弯曲波的速度Table 1 Velocities of bending waves propagating along shell.

Liu 等[8]研究发现，弯曲波运动到扰动源正对位置后反射，而反射返回的弯曲波与入射的弯曲波由于频率相近，两者相互叠加将形成“驻波”现象，具备“驻波”的弯曲振动与膜振动叠加导致球壳上应变呈周期性空间分布。图5 为L=10，62.5 mm 时球壳上应变峰值的空间分布情况，可见两种情况下离扰动源正对位置较近时空间分布周期约为41 mm，离扰动源正对位置较远时空间分布周期约为82 mm。

2 理论分析

Timoshenko 在研究梁弯曲时考虑了剪切效应，推导了梁弯曲波波长与波速之间的关系[9-10]。本文中将参考相关理论，推导球壳弯曲波波长和波速的关系。取图1 中球壳的微元段dx 为对象，其弯曲变形如图6 所示，Q 表示作用在任一截面的切力，M 表示相应的弯矩，α 为球壳的微元段转角。需要说明的是，球壳微元除了受到外力Q、M，还受到切向拉压力以及垂直于纸面的拉压力的作用，但这两个外力对微元弯曲变形没有贡献，因此未显示在图中。

图5 2 种扰动源半径下壳体应变峰值的分布Fig. 5 Total strain along spherical shell for two kinds of disturbance source radii

图6 壳体弯曲变形示意图Fig. 6 Schematic diagram of spherical shell bending deformation

进行了以下假设和限制，将问题简化为可以进行较易求解的数学模型：

(1)弯曲波理论分析是基于平截面假定的，即假设与图6 中虚线垂直的平截面在变形后仍为平面，并保持和变形后的虚线垂直。该假设可把平面弯曲问题简化为一维问题。

(2)壳体发生较小的弯曲变形，即tan α 远小于1。该限制可以将后续分析的某些方程简化为线性方程。

按照动量定理和动量矩定理，可得到两个动力学公式：

式中：ρ0为变形前壳体的密度，A0为变形前壳体截面面积，I 为截面对中性轴的转动惯量，V 为壳体的横向移动速度，ω 为截面转动的角速度。

在平截面假定下，可得到：式中：w 为壳体的横向位移，K 为曲率。

当壳体仅发生较小的弯曲变形时，即tan α 远小于1，可得到：

另外，弯矩M 按其定义是截面上法应力对中性轴的合力矩，为：

式中：E 为材料弹性模量。

Timoshenko 认为，当弯曲波波长较短时，必须考虑剪切效应。对弯曲变形考虑剪切应力时，原来的平截面将因为剪切应变而发生挠曲，即壳体的横向位移w 实际上是由于弯矩M 作用下的转角α 所对应的wM和由于切力Q 作用下的剪切应变γ 所对应的wQ两部分组成：

切力Q 与切应变γ 存在关系：

式中：G 为剪切弹性模量，η 为壳体截面的形状系数。矩形面η 的计算公式为[11]：

泊松比υ=0.29，则η=0.85。

把式（3）和（7）～(11）代入式（1）和式(2)，则得到：

两式消去含wM的项，得到四阶偏微分方程：

式（17）的通解形式为：

式中：D 为弯曲波振幅， ψ =2π/λ 为波数， p=ψC 为圆频率， λ 为波长，C 为波速。则得到：

式（19）存在两个解，即：

当λ 无限小时，式（20）中右边根号内取正号“+”时，得到C=C0；右边根号内取负号“-”时，得到C=CQ。显然，C=CQ才是真解，因此式（19）的真解取：

对于矩形截面，旋转半径计算公式为[12]：

式（23）给出了球壳上弯曲波波速和波长的关系。

3 理论计算结果

由公式（23）得到球壳上弯曲波波速与弯曲波波长的关系曲线，见图7。

图7 球壳上弯曲波波速与波长的关系Fig. 7 Relation between length and velocity of bending wave

可见弯曲波的波长越短，波速越快；当波长λ 无限短时，弯曲波波速趋于极限值。取弹性模量E=200 GPa，密度ρ=78 00 kg/m3，泊松比υ=0.29，其计算为：

可见，最短弯曲波的波速为声速的0.574 倍，约为2 895 m/s，与表1 中由数值仿真得到的结果2 965～3 188 m/s 相近。

另外，波的波速和波长还存在关系：

式中：f 为波的频率。

壳体膜振动主频的表达式为：

弯曲波A 为膜振动频率相近的第一个弯曲波，其频率近似取为壳体膜振动主频，取球壳半径a=261.5 mm，即可计算得到弯曲波A 的频率约为5 173 Hz。

由式（23）和式（25）均可绘制出弯曲波A 的波速和波长之间的关系曲线，如图8 所示，两曲线的交点可以得到弯曲波A 的波速。

图8 弯曲波A 的波长与波速Fig. 8 Length and velocity of bending wave A

由图8 可知，频率与呼吸振动频率相近的弯曲波A 的波速约为376 m/s。而表1 中数值仿真得到的波速为427～443 m/s，计算结果与数值仿真结果相差(11.9～15.1)%。

Liu 等[8]给出了壳体变形空间周期分布的计算公式，离扰动源正对位置较近时其周期为：

离扰动源正对位置较远时其周期为:

式中：C 为波速，取弯曲波A 的计算结果376 m/s；f 为相应的频率，取计算结果5 173 Hz。则式(27)～(28)分别计算得到在离扰动源正对位置较近时周期为36.3 mm，离扰动源正对位置较远时其周期为72.7 mm，与数值仿真结果41、82 mm 相差约11.5%。

可见，理论分析结果与数值仿真结果存在一定的差异，该差异主要来源于两个因素：一个因素是对理论分析的简化处理，理论分析的两个主要假设或者限制条件在数值仿真中并不完全成立；另一个因素是对数值仿真数据的读取方法，计算弯曲波波速要读取弯曲波在球壳上运动的距离和时间，壳体变形空间分布周期要读取球壳上每个周期的长度，这些数据的读取方法本身必然带来误差。

4 结 论

参考Timoshenko 梁的弯曲理论，基于平截面假定和壳体发生较小的弯曲变形的假设，推导出球壳上弯曲波波速和波长的关系，计算得到了最短弯曲波和与膜振动频率相近的弯曲波的波速，再结合早期研究提出的壳体变形分布周期与弯曲波波速的关系，计算得到了壳体变形空间分布的周期。主要结论有：

（1）理论计算结果与数值仿真结果基本吻合，其中弯曲波波速的计算结果与数值仿真结果相差在15%以内，壳体变形空间分布周期的计算结果与数值仿真结果相差在12%以内。

（2）弯曲波的波长越短，波速越快。当波长无限短时，波速趋于一极限值，约为声速的0.574 倍。

本文的理论研究为解剖式分析应变增长现象的新方法提供了一定的理论依据，进一步验证了该分析方法的合理性。