具有恒定冲击载荷的梯度泡沫金属材料设计

多胞材料具有高孔隙率、可变形量大和轻量化特征以及超强的能量吸收能力，广泛地应用于汽车、航空航天和铁路等领域。在动态压缩情形下，多胞材料具有明显的变形局部化特征和应力增强现象[1-3]，可以很好地提高结构的耐撞性。如何优化设计多胞材料以满足结构耐撞性要求，是近几年结构冲击领域的研究热点。

多胞材料的准静态名义应力应变关系呈现明显的三段式特征：弹性段、长而稳定的平台段和压实段。基于此特征，Reid等[1]提出了刚性-塑性-锁定（R-PP-L）模型，并发展了一维冲击波模型以表征多胞材料的动态力学行为。Gibson等[4]定量化给出了R-PP-L模型中多胞材料的平台应力和锁定应变与相对密度和基体材料屈服应力的关系。该理论已被应用于研究不同相对密度和强度分布下梯度多胞材料的吸能能力和峰值应力[5-10]。然而，R-PP-L模型只满足对多胞材料力学性能的一阶近似，Zheng等[11]提出了更为精确的刚性-塑性硬化（R-PH）模型，其应力应变关系写作

式中：σ0和C分别为初始压溃应力和应变硬化参数，这两个材料参数与多胞材料的相对密度相关[12-13]。Ding等[14]发现基于R-PH模型可以更好地指导爆炸牺牲层的设计。

随着更为精确的多胞材料本构模型的提出，梯度多胞材料的动态变形机理的研究也趋于完善。但多数研究停留在讨论不同梯度多胞材料的动态力学性能上，分析具有不同密度和强度分布的多胞模型的动态力学性能[15-21]。只通过观察比较其力学响应，选择满足耐撞性的要求的梯度模型，无法直接达到最优化的目的。Yang等[13]首次提出了梯度多胞材料的耐撞性反向设计理论，直接设定物体冲击力要求，反向设计多胞材料的密度分布，并利用三维细观有限元模拟验证其有效性。不过，该耐撞性反向设计理论无显式解，Yang等[13]运用四阶龙格库塔法获得了数值解。常白雪等[22]利用该理论对二维Voronoi模型进行耐撞性反向设计，并提出了该理论的简单而方便的渐近解。由于二维随机蜂窝与三维闭孔胞元的结构差异，针对于二维模型提出的渐近解，无法适用于三维多胞构型的耐撞性设计需求，而对于闭孔多胞构型的耐撞性设计理论的近似解尚未提出。

本文中进一步研究了质量块以特定的初速度撞击梯度闭孔泡沫金属杆的理论模型和有限元模拟：通过讨论物体冲击过程中受载恒定的耐撞性要求，对梯度多胞材料耐撞性反向设计理论进行简化。针对闭孔泡沫金属，拟合R-PH模型中初始压溃应力和应变硬化参数与相对密度的关系；运用级数法，求解耐撞性反向设计理论，获得梯度多胞材料的密度分布；利用变胞元尺寸法和三维随机Voronoi技术，生成特定密度分布的梯度闭孔多胞模型；利用有限元软件ABAQUS/Explicit模块计算，获取冲击力结果验证设计的有效性，并与均匀密度分布的多胞材料结果进行对比，讨论合理的密度分布设计的必要性。

1 理论与方法

1.1 R-PH模型的材料参数

本文中采用的均匀密度的闭孔泡沫模型与文献[11]中的一致。在10 m/s的恒速压缩下，利用ABAQUS/Explicit模块进行有限元模拟，获取了不同相对密度的Voronoi模型准静态压缩下的力学性能，其名义应力应变曲线如图1(a)所示。进而，利用Zheng等 [11]提出的R-PH模型对准静态名义应力应变曲线进行最小二乘法拟合，如图1(a)所示。所得的初始压溃应力和应变硬化参数与相对密度呈幂律关系[13]，如图1(b)所示。为方便理论分析，本文采用幂次为3/2的幂律关系：

式中：ρ为相对密度，σys为基体材料的屈服应力，k1和k2为拟合参数，分别为1.166和0.115。

图1(a) 不同相对密度下闭孔泡沫模型的准静态应力应变关系
Fig.1(a) Nominal stress-strain relations of closed-cell foam models with different relative densities

图1(b) 无量纲初始压溃应力和应变硬化参数与多胞模型相对密度的幂律关系[13]
Fig.1(b) Power-law fitting of material parameters of the R-PH idealization with relative density[13]

1.2 耐撞性设计模型及其渐近解

在质量块初速度冲击梯度泡沫金属材料的过程中，为保护物体安全，且平稳减速，须控制物体在冲击过程中受冲击载荷恒定且低于物体耐受载荷。因而，梯度泡沫金属杆需根据耐撞性要求进行特定的相对密度分布设计。针对二维随机蜂窝模型的耐撞性反向设计中，常白雪等[22]提出了控制质量块冲击梯度多胞金属杆过程中受载恒定的相对密度分布渐近解，通过二维细观有限元模拟验证了设计的有效性，并指出特定的相对密度分布只适用于特定的冲击情形。对于闭孔梯度泡沫材料，Yang等[13]提出了耐撞性反向设计理论，通过四阶龙格库塔法，针对于控制冲击过程中质量块受载恒定、线性递增和线性递减的耐撞性要求，获得了相对密度分布的数值解，但未给出可方便实际工程应用的公式解。本文中针对物体受载恒定的耐撞性要求，对耐撞性反向设计理论进行简化，运用级数法，求得相对密度分布的渐近解。

考虑一质量为M的物块以冲击速度ν0撞击横截面面积为A0的闭孔泡沫金属杆，被保护质量块在整个冲击过程中受载恒定，额定受载力为F0。冲击过程中质量块的速度为ν=ν0-a0t，其中a0=F0/M为冲击物体的加速度，t为冲击时间。冲击的总时间T0=Mν0/F0，特征长度L0=Mν02/F0。根据Yang等[13]的耐撞性反向设计理论，当冲击载荷为恒定值时，该理论可以简化为

式中：ρ为泡沫金属材料相对密度，ρs为基体材料密度，Φ为冲击波阵面位置，撇号表示对ρ的导数，考虑到相对密度、长度、冲击速度和时间的特征量，定义无量纲量：

式中：ρ0为冲击端的相对密度。此时，式(3)可以进一步写作：

式中利用级数法，可以求得渐近解：

式中：一阶系数和二阶系数分别为：

和

梯度泡沫金属杆冲击端的初始相对密度ρ0可以由一维塑性冲击波模型给出：

在质量块冲击过程中，相对密度单调递增变化的梯度泡沫金属杆中塑性压溃波的传播满足一维冲击波假设，且质量块随波后压实区同步运动，即速度一致[13,22]，波后压实区应力记为σB。根据应力波理论和牛顿运动定律，如文献[22]，可以获得：

根据渐近求解的相对密度分布和对应的冲击情形代入上式，可以获得冲击力等参量随时间的变化。

1.3 梯度泡沫金属材料的数值建模

变胞元尺寸法可用于生成壁厚均一且密度分布连续变化的二维或三维的随机Voronoi细观有限元模型[6-7,13]。与基于六角蜂窝构型的二维多胞模型的撒点规则不同，三维多胞结构的撒点规则是基于卡尔文十四面体提出的。Yang等 [13]基于均匀三维Voronoi结构引入了新的撒点规则，在一定区域内随机撒点，通过已知的相对密度分布和壁厚，控制任意两个相邻成核点i和j之间的间距不小于当前位置的最小许可间距，记为

式中：h为胞元胞壁厚度，k为胞元不规则度，ρ（xij，yij，zij）为i和j两核点连线中点位置对应的相对密度。本文中只考虑x方向上的相对密度分布，xij=（xi+xj）/2，其他两个方向相对密度均匀分布。

梯度泡沫金属的基体材料模型和参数与文献[13]一致。基体材料假设为弹性-理想塑性，其密度ρs=2 700 kg/m3，杨氏模量Es=69 GPa，泊松比ν=0.3，屈服应力σys=165 MPa。采用壳单元S3R划分胞壁网格，通过网格收敛性分析[11]，其壳单元的特征尺寸大小为约0.3 mm。梯度泡沫材料细观有限元模型的横截面为正方形，其面积A0=30×30 mm2。模型中所有可能的接触设置为摩擦因数为0.02的通用接触，冲击加载示例见图2。

图2 质量块冲击梯度多胞杆
Fig.2 A diagram of 3D graded Voronoi models under mass impact

2 结果与讨论

2.1 梯度泡沫金属材料相对密度分布的设计

上文已经给出了恒定冲击力耐撞性反向设计理论的一阶和二阶近似解，本文中采用与Yang等[13]冲击情形一致：冲击物体质量M=50 g，初始冲击速度ν0=100 m/s，额定受力为F0=7.5 kN。由此可得，特征时间T0=66.7 ms，特征长度L0=66.7 mm。冲击端相对密度为0.086 9。参数则有渐近解中的系数为a1=0.570和a2=0.191。

通过已获得的上述参数表征相对密度分布，并与Yang等[13]的结果进行比较验证。根据初始压溃应力与相对密度的关系，获得相对密度梯度分布的渐近解，将Yang等[13]的数值结果作为精确解进行对比，如图3(a)所示。渐近解与精确解在冲击端的相对密度有小量差距。渐近解的结果0.086 9与精确解的结果0.085 0相比，相对误差约为2%，可以忽略。比较渐近解与精确解给出的相对密度分布表明，一阶近似解只是线性近似，冲击端附近差异较小，越靠近远端与精确解差异越大；二阶近似解的整体差异相比较小，冲击端附近稍大于精确解，远端附近稍小于精确解。

结合应力波理论和牛顿第二运动定律，理论预测特定密度分布的梯度泡沫金属材料的力学行为。将相对密度分布和特定冲击情形代入公式进行理论预测，运用四阶龙格库塔法求解，可以获得冲击力、速度和波阵面位置随时间的演化过程。通过观察冲击力随时间的演化过程，二阶近似和精确解模型的结果与耐撞性额定冲击力值非常吻合，而一阶近似模型结果初期吻合较好，后期差别较大，如图3(b)所示。通过冲击速度降为0时确定波阵面传播停止的位置，获得梯度泡沫金属杆的长度分别为：69.4 mm（一阶近似模型）和65.9 mm（二阶近似模型）。精确解给出的梯度泡沫金属杆的长度为66.3 mm，如图4(b)所示。通过理论预测，二阶近似解具有足够的精度，可用于设计具有恒定冲击载荷的梯度泡沫金属材料细观有限元模型的相对密度分布。

图3 相对密度分布的渐近解和冲击载荷历程的理论预测值
Fig.3 Theoretical predictions of asymptotic solutions of relative density distribution and the history curves of impact force

图4 冲击速度与冲击波波阵面位置的历史曲线
Fig.4 Evolution history of impact velocity and location of shock front

2.2 有限元模拟结果

基于已获得的相对密度梯度分布，利用变胞元尺寸法，生成梯度泡沫金属杆的细观有限元模型，并生成均匀密度分布的泡沫金属细观有限元模型作对比。均匀密度模型与二阶近似密度梯度模型的整体尺寸一致，平均相对密度一致，相对密度为0.117，三维Voronoi模型如图5所示。取中轴面x=15mm处剖面图观察其变形过程，可以发现一阶近似和二阶近似模型在刚性板的冲击过程中，以逐层压溃变形模式为主，在密度较小端胞元开始发生压溃，并形成一道较为稳定的塑性冲击波向远端传播，如图6(a)～(b)。而对于均匀密度分布的细观有限元模型，其撞击初始变形以冲击端胞元压溃为主，形成一道向远端传播的冲击波。随着刚性板的继续推进，由于端部效应，杆内应力波在固支端反射加强，其应力超过了支撑端部分胞元的初始压溃应力，随即发生压溃变形。由于冲击情形为质量块初速度撞击，冲击能量有限，端部反射波只能使得少量胞元发生压溃，并不能使得该压溃波传播很远，最终压溃带宽度约为3～4个胞元直径，如图6(c)。比较密度均匀和梯度分布多胞模型变形图，发现合理地设计相对密度分布，可以有效地控制压溃波的形成和传播。

图5 密度梯度多胞杆细观有限元模型
Fig.5 Cell-based finite element models of density gradient cellular rods

图6 梯度多胞细观有限元模型中轴剖面变形图
Fig.6 Deformation patterns of cell-based finite element models in the middle section of the density gradient cellular rods perpendicular to x-axis

梯度泡沫金属杆冲击端发生局部化压溃，产生塑性冲击波向远端传播。波前未压实区应力和速度分别为为σ0和0，其中初始压溃应力为σ0(ρ(Φ(t)))=s0ρ3/2，即当前t时刻波阵面位置处对应的初始压溃应力，则支撑端力Fsup=A0σ0。有限元计算的支撑端和冲击端力的结果与相对密度分布精确解模型的理论预测结果比较见图7。有限元计算的结果出现震荡，其原因为：刚性板初始接触阶段，试件和刚性板间的罚函数接触算法会导致震荡；后期频率较高的震荡为胞壁内部的弹性波来回传播及交互作用导致的；后期低频且幅值较大的震荡是由于胞元的失稳压溃产生的。一阶近似下的有限元模型的计算结果与理论预测值比较：初始阶段有限元计算的冲击力结果偏高于额定冲击力，后逐渐趋于额定值并重合；有限元计算的支撑端力结果与理论预测结果初期呈递增趋势，初始阶段有限元计算的结果高于理论预测值，上升斜率较小，后逐渐低于理论预测值，整体相比差异明显。二阶近似下的有限元模型的计算结果与理论预测值比较表明，冲击力结果整体维持较为稳定的状态，在某一恒定值附近震荡，其值稍大于额定设计值；支撑端力结果与理论预测值初期都呈现递增趋势，初始阶段有限元计算结果稍高于理论预测值，后逐渐趋于理论预测值，整体偏高。二阶近似模型相比于一阶近似模型可以很好地控制冲击力在某一恒定值附近震荡，其值偏高于额定值，与Yang等[13]的结果基本一致，如图8(a)所示。这表明二阶近似解可以很好地表征精确解，而一阶近似解的精度不够。初始压溃应力与相对密度3/2次幂的关系可以应用于指导梯度泡沫材料的耐撞性设计。

与二阶近似模型结果对比，均匀密度分布的泡沫金属模型的冲击力结果具有明显的峰值且高于额定载荷，后随时间衰减。由于冲击初期，质量块冲击速度最高，冲击端相对密度（0.117）大于二阶近似模型（0.086 9），其初始压溃应力和应变硬化参数较高，使得初期冲击力明显高于二阶近似模型，如图8(a)所示。随冲击速度的降低，其对应的冲击力也相应减小，出现明显的冲击力峰值后衰减，而二阶近似模型相对密度单调递增，补偿了因速度衰减对冲击力的影响，很好地控制了冲击力的稳定演化。观察物体冲击速度随时间的衰减过程发现，有限元模拟结果要比理论预测衰减偏快，用时稍短。以冲击速度降为0时表征冲击结束，一阶和二阶近似模型实际冲击时间分别为0.64 ms和0.62 ms，而均匀密度分布模型速度衰减得最快且不稳定，冲击时间为0.56 ms，耗时最短，如图8(b)所示。合理地优化梯度泡沫金属材料相对密度的分布，对于控制质量块冲击力是必要的。

图7 冲击端和支撑端载荷理论预测与有限元计算结果的比较
Fig.7 Comparisons of impact force and support force curves history between theoretical predictions and finite element (FE) results

图8 有限元计算结果对比
Fig.8 Comparisons of finite element (FE) results

图7和图8中有限元模拟的结果偏高于理论预测值，可能原因是本文中采用的冲击波模型未考虑多胞材料的率敏感性影响，初始压溃应力和应变硬化参数均使用准静态应力应变关系拟合得到，导致理论值低估有限元模拟的结果。已有研究表明，即使在未考虑基体材料应变率效应的影响下，多胞材料的动态应力应变关系与准静态的应力应变关系有明显的不同。Zheng等[11]发现了泡沫金属一种特有的率敏感性——“冲击速率敏感性”，即多胞材料动态压实区的应力应变状态对应于特定的冲击速度，不同冲击速度下的状态点可以用一条曲线描述，该曲线可用于表征多胞材料动态下的应力应变关系。通过比较动态应力应变关系和准静态的应力应变关系发现，动态初始压溃应力高于准静态初始压溃应力，而动态应变硬化参数小于准静态应变硬化参数。Ding等[23]和Wang等[24]进一步研究发现了泡沫金属的初始压溃应力具有应变率敏感性，在应变率高到一定程度后二者有幂率关系。因而，后期研究中可以考虑泡沫金属材料率敏感性的影响，以提高设计的准确度。在实验研究方面，He等[25]提出了通过控制水冷前发泡时间可以制备连续密度梯度变化的闭孔泡沫金属，将有助于后期耐撞性反向设计理论的实验验证。

3 结 论

本文中基于一维率无关刚性-塑性硬化（R-PH）模型，获得了泡沫金属的初始压溃应力和应变硬化参数与相对密度呈3/2次幂的关系。针对于恒定冲击力的耐撞性要求，对梯度泡沫金属的耐撞性反向设计理论进行简化，运用级数法，获得了闭孔梯度泡沫金属材料相对密度分布的渐近解，与Yang等[13]的精确解进行了比较。通过理论结果的比较发现：两者初始相对密度具有微小差异，但可以忽略；一阶近似解是对精确解的线性近似，误差较大；二阶近似与精确解吻合较好。通过冲击力理论预测判断，二阶近似密度分布可以达到恒定冲击力要求，而一阶近似结果在后期无法确保冲击力维持恒值。根据相对密度分布的渐近解，设计三维梯度细观有限元模型并进行计算验证。通过与数值模拟结果的比较发现：二阶近似模型冲击力一直维持较为稳定的恒定冲击力，很好地满足了恒定冲击载荷的耐撞性要求；一阶近似模型冲击力无法维持较为稳定的载荷；均匀密度分布的多胞模型冲击力具有明显的峰值，后随时间衰减。相比均匀密度多胞材料，单调递增的二阶近似密度梯度设计，合理地分配密度的强弱，控制载荷稳定演化，使得速度平稳降低，提高了撞击过程物体的稳定性。综上所述，耐撞性反向设计理论简化模型的渐近解对于梯度闭孔泡沫材料的耐撞性设计是有效的，所提出的耐撞性设计方法在控制冲击吸能过程和冲击物受载方面具有指导意义。