刘 冕，王根伟，宋 辉，王 彬

（1. 太原理工大学机械与运载工程学院应用力学研究所，山西 太原 030024；2. 材料强度与结构冲击山西省重点实验室，山西 太原 030024；3. 伦敦布鲁内尔大学机械航空工程系，伦敦 UB8 3PH，英国）

多胞金属兼具质轻和优异的吸能特性，已被广泛用于高速火车、汽车和航空航天等领域。在准静态下，多胞材料的应力-应变曲线呈现弹性、平台和致密化3 个不同的变形阶段[1-2]。在高加载速率下，多胞材料的动态行为可以通过应力增强和变形局部化来表征[3-5]。迄今为止，人们利用冲击波理论和有限元方法来研究具有不同密度梯度的多胞材料，发现密度梯度对多胞动态响应下的力学性能具有重要影响[6-7]，不仅可以影响冲击过程中应力波的传播[8-9]，同时特定的梯度多胞金属可以提高材料的能量吸收能力和抗冲击性能[10-12]。

Reid 等[13]提出用一维刚性-塑性-锁定(Rigid, perfectly-plastic, locking，R-PP-L)冲击波模型来描述多胞材料中应力波的动态响应；随后Tan 等[14]利用R-PP-L 模型对冲击载荷下泡沫材料的动态压缩性能进行了分析；Wang 等[15]通过R-PP-L 模型分析了线性密度梯度下多胞金属的耐撞性，发现正梯度多胞金属和负梯度多胞金属分别对冲击端和支撑端具有良好的保护效果；Shen 等[16-17]基于R-PP-L 模型研究了梯度多胞材料在高速压缩时的冲击响应，发现多胞材料沿正梯度方向压缩时仅在加载端产生一个波阵面（单波模型），沿负梯度方向压缩时则会在两端各产生一个波阵面（双波模型）。由于R-PP-L 模型仅包含平台应力和锁定应变两个材料参数，只能实现材料应力-应变曲线的一阶近似，为了更好地描述多胞材料的力学性能，Zheng 等[18]提出了一种更为精确的、率无关的刚性-塑性硬化（Rigid, plastichardening，R-PH）模型，用以表征梯度多胞材料的抗爆炸和抗冲击特性[19-22]。

有限元模型常被用来研究梯度多胞材料的力学性能[23-25]。Ajdari 等[23]发现，功能梯度蜂窝结构存在3 种不同的压溃模式：准静态、过渡态和动态，沿冲击方向降低相对密度可以增强蜂窝在压溃初期的能量吸收；Fan 等[24]研究了均匀、梯度和随机胞壁厚度的金属空心球泡沫材料的动态破碎响应，发现负密度梯度泡沫具有最大的能量吸收能力，传递到受保护结构的力最小。Zhang 等[25]对受恒速冲击影响下梯度Voronoi 多胞模型的动力学行为进行了模拟分析，定义了梯度多胞材料的第一、第二临界速度。

梯度泡沫金属在R-PH 模型下的应力波传播研究大多集中于正梯度泡沫，对负梯度泡沫的冲击波分析较少，负梯度泡沫冲击波模型结果与有限元结果之间的联系还有待完善。同时，负梯度泡沫材料虽然具有防护支撑端被保护物体且能量吸收能力较强的特点，但支撑端局部密实化现象的发生影响了负梯度泡沫优异的力学性能，因此，局部密实化的影响因素值得分析。

本研究基于R-PH 模型，建立恒速冲击荷载作用下负梯度泡沫材料的一维冲击波模型，给出冲击波传播的基本控制方程；利用随机Voronoi 技术构建梯度泡沫金属材料的三维细观有限元模型，使用LSDYNA 有限元软件对泡沫金属的动态压溃过程进行模拟分析，用于验证冲击波理论结果；通过分析冲击模型的响应历程，定义梯度泡沫材料的局部密实化应变和第二临界速度，并计算不同密度梯度、相对密度下局部密实化应变、第二临界速度的变化规律；最后讨论负梯度泡沫中支撑端局部密实化现象对被保护结构的影响，为工程防护提供参考。

1 理论与方法

1.1 冲击波模型

负梯度泡沫的一维冲击波模型如图1 所示，质量为M 的物块以恒定的中等冲击速度v 撞击闭孔梯度泡沫金属试件，试件的两端同时产生冲击波，且冲击波朝着相反的方向传播，此为双波模型。将冲击端处产生的冲击波定义为前冲击波，支撑端处产生的冲击波定义为后冲击波，前、后冲击波的物理量分别用下标f 和b 表示，在t 时刻，冲击波阵面的拉格朗日坐标分别为Xf(t)、Xb(t)。在波阵面处，波阵面前、后的物理量分别为{vA(t)，εA(t)，σA(t)}、{vB(t)，εB(t)，σB(t)}，两个冲击波阵面之间未变形区域的速度为vu(t)。

定义负梯度泡沫在拉格朗日坐标下的密度分布为

图 1 双波模型Fig. 1 Double shock model

式中：ρ0为平均相对密度，γ 为密度梯度，X 为x 方向的拉格朗日坐标，H 为梯度试件x 方向的总长度。

根据应力波理论[26]，冲击波阵面上的运动学相容条件为

同时，动力学相容条件为

式中：ρs为基体材料密度。联立式(2)、式(3)可得波后应力为

前波波阵面的初始条件为{vfA(t) = vu(t)，εfA(t) = 0，σfA(t)}，{vfB(t) = v，εfB(t)，σfB(t)}，后波波阵面上的初始条件为{vbA(t) = vu(t)，εbA(t) = 0，σbA(t)}，{vbB(t) = 0，εbB(t)，σbB(t)}。

将初始条件代入式(4)，可得前、后波阵面处波后应力为

基于R-PH 模型的刚性假设

近年来，人工造林已成为造林的主要方式，忽视了天然林自我更新和人工促进森林自我更新的重要互补作用。随着国家生产和建设对木材资源的需求日益增加，森林采伐遵循“消费少于增长”的原则，导致人工造林的更新不能满足林业良性发展和管理的需要。追求人工造林成果，不注意森林的自然更新，是造林方式的一大缺陷。

式中：σ0为泡沫材料准静态的初始压溃应力，C 为应变硬化参数，ε 为应变。波阵面的前方区域接近塑性压溃临界状态，且处于应力平衡场，因此波阵面处两种波前应力分别为

由此推出，波后应变为

以中间未变形区域为研究对象，由牛顿第二定律可得

积分后，未变形区域速度为

因此前、后波阵面的坐标为

在冲击波传播过程中，当vu= v，即中间未变形区域的速度增大到冲击速度时，前冲击波停止传播，定义此时刻为t*。当t ＜ t*时，定义为第一阶段，冲击波在该阶段的响应结果可由上述所得方程求解。当t ＞ t*时，第二阶段开始，只有单一的后冲击波向左端冲击端继续传播。当后波阵面与前波阵面相遇时，响应结束，定义此时刻为t2，则第二阶段后波阵面的坐标为

波后应力为

结合上述方程，双波模型中冲击端应力σimp与支撑端应力σsta分别为

1.2 有限元模型

3D-Voronoi 模型是由N 个形核点完全随机分布在体积为V0的立方体区域内生成的，为了防止模型中产生过小胞元，立方体内任意两个相邻形核点之间的距离不小于δ[27]

本研究中模型的平均相对密度ρ0分别为0.06、0.09、0.12，密度梯度γ 分别为-0.4、-0.6、-0.8，尺寸为40 mm × 40 mm × 60 mm，不规则度为0.5，模型胞元个数分别为700（γ = -0.4）、732（γ = -0.6）和768（γ = -0.8），胞壁厚度由给定的相对密度确定。经有限元收敛分析后，使用S3R 和S4R 壳单元混合网格对泡沫模型进行划分，网格尺寸设置为0.2 mm[28]。以ρ0= 0.09、γ = -0.8 为例，负梯度泡沫金属的3D-Voronoi 有限元模型如图2 所示，运用LS-DYNA 软件进行数值模拟。在右端固定一个刚性板，左端刚性板以恒定的冲击速度撞击泡沫模型。基体材料为铝，选定双线性应变强化本构模型，模型参数见表1。其中，ρs为密度，E 为杨氏模量，ν 为泊松比，σys为屈服应力，Et为切线模量。定义泡沫模型为自接触，两端刚性板和泡沫模型之间为面面接触，设定动摩擦系数为0.2。

图 2 有限元模型Fig. 2 Finite element model

表 1 本构模型参数Table 1 Constitutive model parameters

1.3 R-PH 模型的材料参数

R-PH 模型中初始压溃应力σ0、硬化参数C 与相对密度之间存在幂律关系[29]

式中：k1、k2为拟合参数。如图3 所示，在5 m/s 恒速压缩条件下，用最小二乘法拟合不同相对密度的三维均匀Voronoi 泡沫模型的名义应力-应变曲线，得到k1、k2分别为1.25、0.12。

图 3 不同相对密度的准静态应力-应变曲线和参数拟合结果Fig. 3 Quasi-static stress-strain curve and parameter fitting results with different relative densities

冲击波控制方程的初始条件为vu(0) = 0，Xf(0) = 0，Xb(0) = H，这是非线性微分方程组，没有显式解，因此采用四阶Runge-Kutta 法进行求解。

2 结果和讨论

2.1 冲击波历程

两种冲击速度下vu的变化曲线如图4 所示，对应的波阵面位置如图5 所示。由图4 和图5 可以看出，当冲击速度较低(v = 50 m/s)时，t ＜ t*（第一阶段），随着vu逐渐增大，前波阵面坐标Xf由零逐渐增加，后波阵面坐标Xb由H = 60 mm 处开始逐渐减小；当vu增大到50 m/s 时，t = t*，前冲击波阵面停止传播；当t*＜ t ＜ t2时，后冲击波阵面继续传播；在t2时刻，后波阵面到达前波阵面位置，响应结束。当冲击速度较高(v = 150 m/s)时，前冲击波始终存在，且传播距离变大，t*不再存在；在t2时刻，前后波阵面相遇，此时未变形区域的速度vu小于150 m/s。

图 4 不同冲击速度下的vu 曲线Fig. 4 Curves ofvu at different impact velocities

图 5 不同冲击速度下波阵面的位置曲线Fig. 5 Location of wave front at different impact velocities

2.2 局部密实化应变与第二临界速度

在负梯度泡沫的变形过程中会发生支撑端的局部变形，局部变形结束时的现象称之为支撑端的局部密实化，该时刻所产生的应变定义为局部密实化应变[30]。局部密实化应变可以通过能量吸收效率曲线来定义。能量吸收效率[31]是指材料压缩至某一名义应变时所吸收的能量与对应的名义应力的比值，其表达式为

式中：η 为能量吸收效率，εy为初始峰值应力时的应变，ε 为材料压缩时的名义应变，σ(ε)为材料压缩至名义应变时所对应的名义应力。

负梯度泡沫的名义应力-应变曲线与能量吸收效率曲线如图6 所示。由图6 可知，能量吸收效率曲线呈现先上升后下降的趋势，曲线到达最高点时对应的应变为局部密实化应变，此刻冲击端曲线到达屈服应力后的最低点。到达局部密实化应变后，冲击端应力值开始上升，能量吸收效率曲线开始下降。

在不同的冲击速度下，负梯度泡沫存在准静态模态、过渡模态和冲击模态3 种变形模态[32]，其中过渡模态和冲击模态分别对应两种冲击波传播历程。第二临界速度vcr2[33]是指过渡模态向冲击模态转变的临界速度。图7、图8 分别是相对密度为0.09、密度梯度为-0.8 的梯度泡沫在两种模态下的应变云图与名义应力-应变曲线。

图 6 负梯度泡沫的名义应力-应变曲线与能量吸收效率曲线Fig. 6 Nominal stress-strain curve and energy absorption efficiency curve of negative graded foam

图 7 负梯度泡沫的应变云图Fig. 7 Strain nephograms of negative graded foam

图 8 两端应力的FE 结果和理论预测对比Fig. 8 Comparison of FE results and theoretical predictions at both sides stress

由图7、图8 可知，冲击波模型推导出的理论应力值与有限元模型的模拟结果吻合较好，理论模型能够较好地预测负梯度泡沫在冲击载荷作用下的应力变化趋势。若泡沫处于过渡模态(v = 50 m/s ＜vcr2)：在冲击波第一阶段，负梯度泡沫两端同时发生变形；当t = t*时，冲击端停止变形，应力突然下降，小于此处的准静态屈服应力，此时支撑端局部变形结束，t*时刻所对应的应变为局部密实化应变；随后支撑端继续向中间区域变形，两端应力逐渐上升。若泡沫处于冲击模态(v = 150 m/s ＞ vcr2)：惯性效应改变了泡沫的变形特性，在压缩过程中，泡沫主要在冲击端发生变形；当t = t2时，冲击波阵面相遇，支撑端的局部变形结束，t2所对应的应变为局部密实化应变。当v = vcr2时，未变形区域速度vu(t)达到冲击速度v 的时刻刚好等于两波阵面相遇时刻，故t*= t2时所对应的冲击速度为第二临界速度vcr2。

2.3 参数分析

冲击波理论模型与有限元模型得到的局部密实化应变随密度梯度和相对密度变化的曲线如图9所示。由图9 可知，两种模型计算结果的最大误差小于15%，说明冲击波理论模型能够较好地预测负梯度泡沫中的局部密实化应变。在不同冲击速度下，局部密实化应变存在3 个增长趋势：当冲击速度较小时，局部密实化应变随冲击速度的增加缓慢增大；当冲击速度中等时，局部密实化应变增长幅度上升；当冲击速度较高且超过第二临界速度时，局部密实化应变增长逐渐平缓且无限接近压实应变。通过比较相同冲击速度下密度梯度与相对密度对局部密实化应变的影响发现，随着密度梯度绝对值和相对密度的增大，局部密实化应变逐渐减小，且密度梯度对局部密实化应变的影响比相对密度更大。

图 9 密度梯度与相对密度对局部密实化应变的影响Fig. 9 Influence of density gradient and relative density on local densification strain

不同密度梯度和相对密度下负梯度泡沫的变形模态如图10 所示，第二临界速度将负梯度泡沫的变形模态图分为两个区域，分别对应速度中等时的过渡模态和速度较高时的冲击模态。从图10 可以看出，更大的密度梯度和相对密度会延迟泡沫模态的转变，随着密度梯度绝对值和相对密度的增加，由过渡模态向冲击模态转变所需的速度越高，第二临界速度越大。

图 10 不同密度梯度和相对密度下的变形模态Fig. 10 Deformation modes under different density gradients and relative densities

2.4 局部密实化的影响

图11 给出了两种速度下密度梯度和相对密度不同时支撑端的应力-应变曲线。为了消除相对密度大小对应力增长幅度的影响，图11(b)和图11(d)中，用支撑端应力分别除以对应的相对密度。由图11可以看出，密度梯度和相对密度对支撑端具有显著影响。密度梯度绝对值越大，泡沫的初始屈服应力越小；相对密度越大，泡沫的应力值越大；密度梯度绝对值和相对密度越大的泡沫，支撑端应力值增长越早，且应力增长幅度越大。当v = 50 m/s 时，γ = -0.8 和ρ0= 0.12 的梯度泡沫最早发生应力增长，且拥有最大的应力增长值；当v = 150 m/s 时，支撑端平台阶段变长，应力逐渐平稳，γ = -0.4 的梯度泡沫应力增长现象消失，而其他梯度泡沫的后期应力仍有增加。

图 11 不同冲击速度下支撑端的名义应力-应变曲线Fig. 11 Nominal stress-strain curve of the support end under different impact velocities

图12 为相对密度为0.09、密度梯度为-0.6 的梯度泡沫有限元模型在不同冲击速度下冲击端与支撑端的名义应力-应变曲线。当冲击速度为30 m/s 时，局部密实化应变为0.11，当变形到达局部密实化应变后，负梯度泡沫的冲击端与支撑端应力曲线上升，其中支撑端应力值由2.1 MPa 持续上升到5.3 MPa；随着冲击速度的增加，局部密实化应变逐渐增大，冲击端应力上升值减小，而支撑端的应力增长明显，应力增长幅度加大；当冲击速度为80 m/s 时，在局部变形过程中，支撑端应力值逐渐上升，最高上升到6.6 MPa；当冲击速度超过第二临界速度，泡沫进入冲击模态，冲击端应力值在初始峰值应力后始终处于逐渐下降的变化过程，支撑端应力值逐渐平稳。

图 12 不同冲击速度下冲击端与支撑端的名义应力-应变曲线Fig. 12 Nominal stress-strain curves of the impact end and the support end under different impact velocities

支撑端的局部密实化现象导致被保护物体受到较高强度的载荷，影响了负梯度泡沫优异的抗冲击性能，对被保护物体不利。在保证能量吸收满足冲击载荷作用的前提下，适当地减小相对密度和密度梯度，能延缓支撑端应力值的增长；同时，冲击速度高于第二临界速度可使支撑端处于应力平台状态，确保被保护物体在压缩过程中承受较小的载荷。因此，可以利用由冲击波理论得到的局部密实化应变与第二临界速度，有效地预测负梯度泡沫在不同材料参数下两端应力值增长对应的变形范围和速度范围，改善不同工况下负梯度泡沫材料作为抗压材料的防护效果。

3 结 论

利用R-PH 理论设计了恒速冲击载荷下负梯度泡沫材料的一维冲击波模型，推导了冲击波传播的基本控制方程。利用LS-DYNA 有限元软件对三维Voronoi 有限元模型计算模拟得到的结果，对冲击波的理论分析进行验证。根据冲击波理论定义了局部密实化应变和第二临界速度，并探讨了冲击速度、密度梯度、相对密度参数的影响，得到如下结论。

(1) 冲击波模型的理论解与有限元模型的数值解吻合较好，基于R-PH 模型的冲击波理论能较好地预测负梯度泡沫金属的应力-应变曲线、局部密实化应变和第二临界速度。

(2) 局部密实化应变存在3 个增长过程：当速度较小时应变增长比较缓慢，当速度中等时应变快速增长，当速度较高时应变增长逐渐平缓；密度梯度绝对值与相对密度越大，局部密实化应变越小，且支撑端的应力增长现象越早发生，应力增长越多。

(3) 密度梯度绝对值与相对密度越大，第二临界速度越大；当冲击速度小于第二临界速度时，速度越大，支撑端的应力增长越多；当冲击速度大于第二临界速度时，速度越大，支撑端的曲线越平稳。