冲击载荷作用下柱壳链中的弹性波传播简化模型及其解析解*

彭克锋，崔世堂，潘 昊，郑志军，虞吉林

（1. 中国科学技术大学近代力学系中国科学院材料力学行为和设计重点实验室，安徽 合肥 230027；2. 北京应用物理与计算数学研究所，北京 100088）

由颗粒、柱壳或球壳组成的链结构能够实现应力波的传输，并具备较好的波形操控能力，如波形的放大和衰减[1-3]，因此可作为声学透镜[4]、减震材料[5-6]、声学开关[7]等，在工程领域有广泛的应用前景。链结构中的波传播过程因离散特性而显得十分复杂，一维链系统的波传播行为可能呈现线性、弱/强非线性的特征[8]，不同的波传播特征主要受组成链的单个胞元基体材料性质和胞元结构（如胞元形状和胞元尺寸）的影响[9-12]。其中，由柱壳组成的链结构，因胞元中空且壁厚可控，在操控波形方面表现出优异的性能[9-10,13-14]。因此，了解柱壳链中应力波的传播规律及弥散机理，对设计新型高效的波形控制器具有重要意义。

建立链结构中相邻胞元之间的接触关系，是研究柱壳链中波传播过程的基础。在一维均匀颗粒链中，Hertz 定律常用来表征相邻颗粒之间接触力F和相对位移δ 的关系，如F∝δ1.5[15]。对于空心壳，接触力和相对位移的关系为F∝δn，其中n为表征非线性程度的接触指数，取值与空心壳壁厚及形状相关。对于球壳，n通常为1.2～1.5[10,16]。Kim 等[9]指出，对于椭圆柱壳，可能n＜1，从而表现出应变软化的特征。

基于链结构中相邻胞元的接触关系，可以用离散元方法求解链结构中的波传播过程，即对每个胞元的力学守恒方程进行单独求解[9,17]。为了方便分析以获得解析解，也可将离散方程连续化后采用应力波理论来研究[18-22]。Pochhammer-Chree 精确理论可给出谐波激励下弹性圆杆中应力波传播的精确解，但其控制方程形式复杂，尤其对于瞬态冲击载荷，方程解析求解往往十分困难[23]。Rayleigh-Love 近似理论中，假设杆在受到加载时，其纵向运动会带来杆的横向涨缩，并且假设杆中质点的横向位移沿径向线性分布，从而推导出计及横向惯性效应的修正的波动方程[22]。Rayleigh-Love 理论具有较高精度，并成功地用于研究瞬态加载边界条件的波传播问题，如阶跃应力脉冲[24]、梯形应力脉冲[25]和指数衰减的边界速度载荷[26]。质量块初速冲击是一种生活中常见的冲击加载情形，但因其边界载荷随时间变化且需要根据牛顿运动定律进行反解，导致方程求解困难，尚未见相应的解析解。

本文中，在理论上分析冲击载荷下柱壳链中的弹性波传播过程及其几何弥散特性，将柱壳链等效为弹性均质方杆，并采用细观有限元方法构建柱壳链的细观结构对理论模型进行验证。首先，基于Rayleigh-Love 横向惯性修正理论和量纲分析方法，获得柱壳链在质量块初速冲击下的无量纲控制方程。再进一步通过Laplace 积分变换方法，对该控制方程进行求解。然后，将理论预测结果与细观有限元结果进行比较，并讨论该理论结果对于不同壁厚柱壳链的适用性。最后，基于理论模型的解析解，讨论柱壳链在初速冲击载荷下的波传播过程，和泊松比、惯性半径、冲击质量和速度对波形弥散的影响。

1 理论分析和数值模拟

1.1 问题描述和模型简化

考虑一个刚性质量块M以初速度v0沿X方向撞击一个静止的柱壳链的一端，柱壳链的另一端固定，如图1(a)所示。所有柱壳排布成一列，相邻柱壳相互绑定，链总长度为L，含有N个柱壳。柱壳的壁厚为h，外径为D=L/N，垂直方向厚度为a。柱壳的基体材料为弹性材料，密度为ρs，弹性模量为Es，泊松比为νs。在X方向压缩时，柱壳在Y方向将产生不可忽略的变形，即Y方向的泊松比较大，而在面外Z方向的变形量很小，可以忽略，即Z方向的泊松比为0。

在质量块的冲击下，柱壳链中的应力波从近端传向远端。柱壳链的几何非线性将导致复杂的非线性应力波传播问题。为简化分析，假设柱壳链为由均匀的连续介质组成的杆，杆在X方向的长度为L，在Y方向的高度为D，在Z方向的宽度为a，如图1(b) 所示。杆的横截面面积为A=aD，等效密度为ρ0=πρs(D-h)h/D2，等效弹性模量为E，在Y方向的等效泊松比为ν，而在Z方向的等效泊松比为0。

图 1 柱壳链简化模型Fig. 1 Simplified model of a cylindrical shell chain

1.2 理论分析

1.2.1 计及横向惯性的波传播模型

进一步得出应变和无量纲速度：

值得注意，式(14)～(19)均为无穷级数，其数值计算的结果受到求和项数目k的影响，k越大，结果越精确。因此，后续的分析为保证结果的准确性，k需足够大，本文中均取k＞50。

1.3 数值模拟

1.3.1 细观有限元模型

采用细观有限元方法，建立质量块初速冲击柱壳链的有限元模型。对每个柱壳单独划分网格，再导入有限元软件ABAQUS/Explicit 进行求解。本文中，柱壳链几何参数取值分别为：单个柱壳的壁厚h=0.56 mm，外径D=10 mm，厚度a=4 mm，柱壳链中共含30 个柱壳，链总长L=300 mm。通过网格收敛性分析，网格尺寸取为0.3 mm，单个柱壳含1 920 个六面体单元(C3D8R)。模拟中，所有的面接触方式均采用硬接触，设置摩擦因数为0.2。柱壳基体材料密度ρs=966 kg/m3，弹性模量Es=1 600 MPa，泊松比νs=0.3。柱壳链近端受到刚性质量块（M=1 g）初始速度v0=1 m/s 的冲击，远端固定。不同时刻柱壳链中的应力云图如图2 所示。结果表明，应力波从冲击端往支撑端沿着柱壳逐个传播，且应力主要集中在柱壳之间的接触区域。

图 2 不同时刻柱壳链中的应力云图Fig. 2 Von Mises stress distributions in the cylindrical shell chain at different times

图 3 名义应力应变曲线和等效泊松比Fig. 3 Nominal stress-strain curves and equivalent Poissons ratios

2 结 果

2.1 位移和载荷曲线

质量块（M=1 g）以初速v0=1 m/s 冲击壁厚0.56 mm 的柱壳链，其冲击端位移随着加载时间的变化，如图4(a)所示。结果表明，冲击端位移随着时间增大，然后逐渐趋于平缓，理论分析与数值模拟结果基本一致。理论分析与细观有限元模拟所得的支撑端和冲击端的载荷历史曲线如图4(b)所示。观察发现，理论分析能够较好地捕捉载荷曲线的整体特征，但曲线振荡存在一些差异，理论计算的载荷曲线呈现低频振荡，但在数值模拟中，载荷曲线呈现高频振荡，这是由于罚函数接触算法所致。由应力波理论，对于不考虑几何弥散的均质弹性杆，支撑端与冲击端峰值载荷的比为2。而该算例中，由理论分析得到的支撑端和冲击端的峰值载荷分别为2.1、1.3 N，支撑端与冲击端峰值载荷的比小于2。因此，柱壳链可作为一种有效的冲击缓冲器。

2.2 应变和速度分布

对于第2.1 节的算例，理论分析和数值模拟给出的应变场和速度场如图5 所示。需要说明的是，在数值模拟中，难以计算沿冲击方向连续变化的宏观应变场，因此定义一些局部点上的应变以检验理论结果。在过柱壳中心的YZ横截面上，取位于其表面的所有点作为该柱壳的参考点，提取这些点在冲击过程中的平均位移和速度，进一步通过相邻两个柱壳参考点的相对位移ΔX计算应变，即 ε =ΔX/D。为方便比较，理论上也可以计算与数值模拟中对应参考点处的速度和应变。图5(a)为不同时刻的应变分布，结果表明，随着波往远端传播，应变峰值略有减小，且随着时间的推移，波后振荡的区域越来越大。理论计算与细观有限元模拟所得结果整体趋势一致，波后的波形振荡和波形前沿宽度也大致相同。不同时刻的速度分布如图5(b)所示。从速度分布可以发现，峰值速度随着波的传播逐渐降低，理论计算的粒子速度与细观有限元结果较好吻合。

图 4 冲击端的位移、冲击端和支撑端的载荷Fig. 4 Displacements at the support end, forces at the impact and support ends

图 5 不同时刻杆中的应变和速度分布Fig. 5 Strain and velocity distributions in rod at different times

当质量块（M=1 g）冲击速度为1 m/s 时，柱壳链中的应变较小，增大冲击速度可以检验应变较大时理论模型的适用性。图6 为冲击速度为5、10 m/s 时柱壳链中的应变分布。结果表明，理论分析与数值模拟的结果较好吻合，冲击速度越大，峰值应变越大，表明理论能够很好地预测峰值较大应变(约0.1)时柱壳链中不同时刻的应变分布。

图 6 不同冲击速度下杆中的应变分布Fig. 6 Strain distributions in rod under different impact velocities

进一步比较不同壁厚柱壳链的理论分析与数值模拟的结果。通过理论分析与细观有限元模拟计算，壁厚0.33、1.07 mm 的柱壳链在受到质量块冲击(M=1 g,v0=1 m/s)时，链中的应变分布如图7 所示。结果表明，理论与数值模拟很好吻合，理论能够很好地预测不同壁厚柱壳链中不同时刻的应变分布。

图 7 不同壁厚时杆中的应变分布Fig. 7 Strain distributions in rod with different wall thicknesses

3 讨 论

3.1 泊松比和无量纲惯性半径对波传播过程的影响

在理论分析中，可以发现泊松比ν 和无量纲惯性半径 κ¯ 总是以组合的形式出现，因此可以直接讨论无量纲参数 α =ν2κ¯2的影响。但是，为了更直观地考察泊松比ν 和无量纲惯性半径 κ¯ 的影响，这里分别改变这两个参数。图8(a)为不同泊松比某个时刻杆中的应变分布，其中无量纲惯性半径为0.04，泊松比为0、0.25、0.50、0.75 和1.00。结果表明，泊松比为0 时，应变分布曲线光滑无振荡，应变波形前沿宽度很窄。当泊松比不为0 时，应变分布存在明显的振荡现象，且随着泊松比的增大，振荡幅值增大，应变波形前沿宽度增宽。此外，应变峰值也随着泊松比的增大而逐渐降低。以应变峰值处作为波阵面位置，泊松比越大，波速越小。类似地，固定泊松比为某个值，如ν=0.5，讨论不同惯性半径 κ¯ 对变波形的影响。 κ¯ 受柱壳半径的影响，图8(b)中 κ¯ 由小到大分别为0.02、0.04 和0.06。与泊松比的影响类似，惯性半径越大，应变峰值越小，波速越小，振荡幅值越大。综上可见，一维杆中横向惯性效应的影响越大，应变分布振荡越剧烈，并且波形前沿宽度增宽、应变峰值减小。

图 8 泊松比和惯性半径对应变分布的影响Fig. 8 Influences of Poisson’s ratio and inertia radius on strain distributions

3.2 无量纲速度和质量对波传播过程的影响

无量纲速度v¯0和质量m¯ 对应变波形的影响如图9 所示。可以看出，随着无量纲速度的增大，应变波形的振动波长不变，应变最大位置处对应的波速不变，波形前沿宽度不变，变化的仅是应变幅值。其实，由式(14)可以发现，无量纲速度仅出现在影响应变幅值的项。不同于无量纲速度，无量纲质量不仅影响应变幅值，也影响波形的波长。无量纲质量越大，应变幅值越大，但由图9(b)可见，质量在一定范围内变化时，对波形前沿宽度、应变最大位置处对应波速的影响均较小。

图 9 无量纲速度和无量纲质量对应变分布的影响Fig. 9 Influences of dimensionless velocity and mass on strain distributions

3.3 弥散效应的分析

4 结 论

建立柱壳链结构的等效连续介质模型，对冲击载荷作用下柱壳链中的应力波传播过程及其几何弥散特性进行了分析，并用细观有限元模型对理论模型进行验证。

（1）基于考虑横向惯性修正的Rayleigh-Love 波动方程，分析了几何弥散对应力波传播的影响，给出了柱壳链在质量块初速冲击下的无量纲控制方程。基于Laplace 变换及逆变换，获得了该控制方程的解析解。

（2）解析解的载荷和位移曲线、应变和速度分布均与细观有限元结果较好吻合，且适用于不同壁厚的柱壳链。

（3）讨论了泊松比、惯性半径、冲击质量和速度对冲击载荷作用下柱壳链中波形弥散的影响。结果表明，应变峰值、波形振荡幅度、波形前沿宽度和应变最大位置处对应的波速主要与泊松比和惯性半径相关，泊松比和惯性半径越大，应变峰值越小，波形振荡幅度越大，波形前沿宽度越宽，应变最大位置处对应的波速越小。

本文中建立的冲击载荷作用下柱壳链中的弹性波传播简化模型及解析解，可为揭示柱壳链的弥散机理和波形控制器的优化设计提供参考。