多胞材料在高速冲击下具有变形局部化的典型特征。Reid等^[1-2]研究了由圆环组成的一维系统的动态冲击行为，提出了“结构冲击波”的概念，并在对木材的动态压溃实验中证实了压溃带的冲击波传播特征。Zheng等^[3]、Liu等^[4]在对随机蜂窝的动态冲击中观察到3种变形模式，即均匀模式、过渡模式和冲击模式。张新春等^[5]、胡玲玲等^[6]探讨了胞元构型对蜂窝结构动态冲击性能的影响。Zou等^[7]认为蜂窝中传播的冲击波波阵面在高速冲击下具有胞元尺寸的宽度。Liao等^[8]采用了局部应变场计算方法表征了蜂窝铝中一维冲击波的传播并得到了冲击波速度的大小。一系列的一维冲击波模型被发展来表征多胞材料中塑性压溃行为，如R-PP-L模型(率无关、刚性-理想塑性-锁定假设)^{[2, 9-10]}和R-PH模型(率无关，刚性-塑性硬化假设)^[11-13]。然而，大部分的冲击波模型都是基于准静态名义应力-应变曲线得到的，直接得到动态的应力应变关系在实验上是困难的。

多胞材料在动态冲击下存在变形局部化和应力增强等现象，名义应力应变曲线失去了物理意义^[4]。Zheng等^[13]通过对泡沫铝的三维细观有限元模型的分析得到了波后应力应变状态点，并提出了D-R-PH(动态，刚性-塑性硬化)模型来表征动态应力-应变行为。Barnes等^[14]、Sun等^[15]分别对开孔泡沫和闭孔泡沫进行实验和数值研究，也得到了多胞材料的动态应力应变状态。但这些研究都是针对冲击模式下的应力应变状态，对过渡模式下的应力应变状态的认识仍然不清楚。此外，文献[13-15]中波后应力都是由撞击端的名义应力得到的，缺少对材料内部局部应力信息的直接表征。

认识多胞材料在中速过渡模式下的应力应变行为，对于完整地理解多胞材料的动态本构关系有重要意义。本文中将采用截面应力计算方法研究基于细观有限元模型的随机蜂窝结构在恒速冲击下的应力分布，分析冲击波在随机蜂窝中的传播特性，表征随机蜂窝的动态应力应变关系。

1.   数值模型

1.1   细观有限元模型

采用二维Voronoi技术构造随机蜂窝结构，详细建模过程参见文献[3]。如图 1，蜂窝试件的长和宽均为50 mm，试件的不规则度为0.3，由500个胞元组成。试件相对密度ρ₀/ρ_s= 0.1，其中ρ₀为蜂窝试样变形前的初始密度，ρ_s为基体材料的密度。蜂窝试件的平均胞元尺寸约为2.5 mm，其值定义为与平均胞元面积相等的圆的直径。

图  1  Voronoi蜂窝试件细观有限元模型

Figure  1.  A cell-based finite element model of a Voronoi honeycomb specimen

下载: 全尺寸图片 幻灯片

采用ABAQUS/Explicit有限元软件对蜂窝材料的恒速压溃过程进行数值模拟。胞壁材料由S4R壳单元(4节点，减缩积分、沙漏控制、有限膜应变)进行模拟。通过网格收敛性分析将模型划分了10 480个壳单元，平均单元长度约为0.2 mm。基体材料属性设置为弹性-理想塑性，材料参数设置为泊松比ν = 0.3，弹性模量E= 66 GPa，屈服应力σ_y = 175 MPa，基体密度ρ_s = 2 700 kg/m³。所有接触面均定义为通用接触，摩擦因数取为0.02。为了模拟面内变形，模型中所有节点的面外变形均被限制。

蜂窝试件的左右两端设置两个刚性面，左侧刚性面以恒定不变的速度v沿x方向向右压缩试件。在压缩过程中右侧刚性板为固支端，如图 1所示。

1.2   截面应力计算方法

对于蜂窝试件在承受匀速压缩的情况下，一维近似是适用的，因此可以通过截面工程应力来表征局部应力信息，详见文献[16]。不同的拉格朗日截面位置的工程应力由截面内力除以截面初始面积得到。截面的内力由两部分组成，一部分是通过基体材料的单元节点所传递的力，称为节点力，记为F_x⁽ⁿ⁾；另一部分是由变形中胞壁接触所引起的接触力，记为F_x^(c)。这两部分力由下式可以计算得到:

式中:M为所考察的截面所贯穿的单元在参考构型中左侧的节点数，K是从撞击端到所考察截面的总的节点数目。对于恒速压缩情形，K包含了刚性端的节点。接触力F_x^(c)的简化计算公式的算法得益于牛顿第三定律，完整的接触对的接触力大小相等方向相反，互相抵消之后剩余的部分为被截面所隔开的接触对的接触力。由于刚性面与试件撞击端位置重合，因此撞击端面上的接触力由刚性面输出接触力得到。F_xi⁽ⁿ⁾、F_xi^(cn)、F_xi^(ct)分别为任意节点编号i的节点力的x分量以及节点的法向和切向接触力的x分量，均可以直接从ABAQUS的输出文件提取得到。截面的工程应力可以由下式计算为：

式中:σ_x⁽ⁿ⁾=F_x⁽ⁿ⁾/A₀，表示由节点所传递的应力；σ_x^(c)=F_x^(c)/A₀，表示由胞壁接触所引发的应力。其中A₀为截面的初始面积。

2.   结果

2.1   一维应力分布

冲击速度分别为200和50 m/s时的应力分布如图 2所示。在任意名义应变时刻下，应力会在某个位置发生陡然减小，并且这个位置随着名义应变的增大往试件右侧移动。这说明在高速和中速撞击下试件内均存在冲击波的传播，应力发生急剧变化的位置即为冲击波波阵面的位置，中速撞击下的波阵面区域相较于高速撞击更加平缓。从图 2可以明显地观察到应力增强现象，跨过冲击波阵面，应力由波前应力迅速增加到波后应力并几乎保持不变。

图  2  两种冲击速度下的应力分布图

Figure  2.  One-dimensional stress distributions at two impact velocities

下载: 全尺寸图片 幻灯片

2.2   基于应力分布得到冲击波速度

可以由应力分布梯度来确定冲击波波阵面Φ的位置，如图 3所示。应力分布梯度反映了一维应力在加载方向上的变化，在应力梯度绝对值达到最大的拉格朗日坐标位置，应力变化最快，这个位置即为当前时刻冲击波波阵面的位置。中速和高速压溃下波阵面位置均可以采用此方法得到。不同冲击速度下得到的波阵面位置随冲击时间变化关系如图 4所示。在恒速压缩情形下，冲击波的波阵面位置与冲击时间关系为线性关系。可通过对线性关系拟合得到的斜率来预估冲击波速度。从图 4中可以看出，冲击速度越大，斜率越大，则相应的冲击波速度也越大。

图  3  一维的应力分布及其应力梯度分布

Figure  3.  One-dimensional stress distributions and the corresponding stress gradients

下载: 全尺寸图片 幻灯片

图  4  不同冲击速度下冲击波位置随时间的关系

Figure  4.  Variation of shock front position with impact time at different impact velocities

下载: 全尺寸图片 幻灯片

2.3   一维冲击波理论

除了由有限元结果直接得到的应力信息来得到冲击波速度之外，一维冲击波理论也可以对蜂窝中的冲击波传播行为进行描述。由应力波理论^[17]，跨波阵面的质量和动量守恒关系分别为:

式中:v_s(t)是冲击波速度, σ_B、ε_B、v_B为波后方的应力、应变和速度，σ_A、ε_A、v_A为波前方的应力、应变和速度。联立这两个守恒关系可以得到冲击波速度与波前波后应力、应变的关系为：

如果能直接从有限元结果中得到波前和波后应力/应变信息，就可以由式(5)计算得到冲击波速度的大小。为了简便，这里将这种结合冲击波理论和数值模拟结果的方法叫做一维冲击波理论方法。

通过对一维应力和应变分布求平均值，可以得到波前和波后的应力和应变，如图 5所示。图中一维应变分布采用文献[8]提出的局部应变计算方法得到。在求平均的过程中，为了消除冲击波波阵面平均效应的影响，将紧邻波阵面附近1.5倍胞元半径R_c内的数据排除掉。波后应力和应变信息由冲击端到Φ-1.5R_c区域内平均得到，波前应力应变信息由Φ+1.5R_c到支撑端区域平均得到。再由式(5)就可以得到不同冲击速度下的冲击波速度。

图  5  一维应力和应变分布及其理想化

Figure  5.  One-dimensional stress and strain distributions and the idealizations

下载: 全尺寸图片 幻灯片

2.4   冲击波速度的比较

通过截面应力和局部应变场等有限元方法以及一维冲击波理论方法都可以得到冲击波速度。此外，也可由已有的冲击波模型得到冲击波速度的大小。由R-PP-L模型^[2]得到的冲击波速度为:

式中:ε_L为多胞材料的锁定应变。由R-PH模型^[13]得到的冲击波速度为:

式中:C为多胞材料的塑性硬化参数。在本文中，ε_L=0.64，C=0.226 MPa。

由不同方法得到的冲击波速度如图 6所示。由一维冲击波理论方法以及R-PH冲击波模型得到的冲击波速度比较接近有限元方法得到的冲击波速度，并且与冲击速度的差值均为一个常值，即材料冲击参数。然而由R-PP-L模型预估得到的冲击波速度明显高于有限元结果。这是由于R-PP-L模型采用的是锁定不变的压实应变，实际上其值小于动态压溃的压实应变。由质量守恒关系可以发现，压实应变越小，冲击波速度越大。R-PH模型中的压实应变是随冲击速度变化的动态应变，因此更加适合描述多胞材料在动态压溃下的冲击波传播。

图  6  不同方法得到的冲击波速度的比较

Figure  6.  Comparison of shock wave speeds obtained using different methods

下载: 全尺寸图片 幻灯片

3.   分析和讨论

3.1   冲击波速度的一致近似

前文中通过应力分布得到了高速冲击下的冲击波速度与冲击速度的关系。在中等冲击速度下，仍然可以观察到冲击波传播的现象^{[4, 8]}，因此仍然采用应力梯度分布方法来确定冲击波速度的大小，结果见图 7。当冲击速度很高的时候，冲击波速度与冲击速度呈线性关系并且两者之间的差值几乎为一个常数^[13]。但随着冲击速度的减小，冲击波速度趋近于常数。因此，冲击波速度对冲击速度的一阶导数在低速和高速下分别趋于0和1，从低速到高速的完整曲线具有明显的S型特征。利用S函数(sigmoid function)，冲击波速度对冲击速度的一阶导数可以写作：

图  7  冲击波速度与冲击速度的关系

Figure  7.  Variations of the shock wave speed with the impact velocity

下载: 全尺寸图片 幻灯片

式中:a和b为待定参数。上式的积分形式可以写作：

式中:c为待定参数。可用式(9)对有限元结果进行拟合，从图 7可见拟合效果较好，其中拟合参数为a=3.603 m/s, b=40.68 m/s和c=54.91 m/s。若式(8)右边的值较小，对应于低速冲击，此时并不能观察到冲击波传播的现象，因此不妨取该值为0.05，可以估算从均匀模式向过渡模式转变的临界速度v_c1=b－aln19=30.07 m/s。若式(8)右边的值较大，对应于高速冲击，不妨取该值为0.95，可以估算从过渡模式向冲击模式转变的临界速度v_c2=b+aln19=51.29 m/s。可见参数a和b的线性组合反映了模式转变速度的大小。当0＜v≪b时，冲击波速度v_s近似等于c。当v≫b时，$1 + {{\rm{e}}^{\frac{{v - b}}{a}}} \approx {{\rm{e}}^{\frac{{v - b}}{a}}}$，此时冲击波速度v_s=v-b+c，和R-PH模型中得到的冲击波速度与冲击速度变化关系一致，即c₁=c-b。因此，式(9)表征了冲击波速度与冲击速度的一致近似关系。

3.2   动态应力应变状态的一致近似

在本文中考察的恒速冲击情形中，若取ε_A=0, v_A=0, v_B=v和σ_A=σ₀，由跨波阵面的守恒关系式(3)和(4)有:

将式(9)代入上式，可得：

与有限元结果的比较如图 8所示，图中虚线对应于式(11)。可见式(11)与有限元结果存在一定的误差。事实上，在恒速冲击下，冲击波波阵面前方的应变和粒子速度不能直接忽略。

图 9给出了不同冲击速度下试件中的速度分布情况，局部粒子速度在冲击波波阵面位置附近约两个胞元区间内急剧变小。当冲击速度为中等大小时，波阵面位置处的粒子速度变化更为平缓。整个波前区域处于低速状态，近似呈线性分布，与冲击速度的大小关系不大。冲击波前方的粒子速度和变形不能忽略，可以通过图 9估算v_A≈10 m/s，初始压溃应变ε_A取准静态实验中常用的ε₀ =0.02。因此，由跨波阵面的守恒关系式(3)和(4)有:

图  8  波后应变和波后应力随冲击速度的变化

Figure  8.  Variations of strain and stress behind shock front with impact velocity

下载: 全尺寸图片 幻灯片

图  9  不同冲击速度下的速度分布

Figure  9.  Velocity distributions at different impact velocities

下载: 全尺寸图片 幻灯片

将式(9)代入上式，可得:

如图 8所示，这个结果与有限元结果吻合较好。引入无量纲参数$\xi = {{\rm{e}}^{\frac{{v - b}}{a}}}$，式(13)可以进一步写作：

式中:α=v_A/a, β =b/a, γ=c/a和B=ρ₀a²。因此，式(14)给出恒速冲击下动态应力应变状态的一致近似关系。准静态和动态应力应变关系的比较如图 10所示。在图 10中，准静态应力应变曲线及动态应力应变状态点由有限元计算得到，动态应力应变曲线由式(14)得到。在动态应力应变关系的曲线上的每个点都对应着特定的一个冲击速度。在高应变情形下，动态应力应变关系位于准静态曲线右侧，这与已有文献中的结论一致。随着应变的减小，动态的应力应变关系逼近准静态应力应变曲线。

图  10  准静态及动态应力应变关系

Figure  10.  Quasi-static and dynamic stress-strain relations

下载: 全尺寸图片 幻灯片

4.   结论

采用细观有限元模型研究了随机蜂窝结构在恒速冲击下的动态性能，采用截面应力计算方法得到了随机蜂窝试件内的应力分布和冲击波速度，比较了高速冲击下由不同方法得到的冲击波速度与冲击速度的关系。结果表明，R-PP-L模型高估了冲击波速度，但R-PH模型以及一维冲击波理论得到的冲击波速度与有限元模拟得到的结果接近。在冲击速度很高的时候，冲击波速度与冲击速度的关系趋于线性，但随着冲击速度的减小，冲击波速度不断减少并趋于常数。最后，发展了可以表征冲击波速度与冲击速度的关系的一致近似模型，并基于一维冲击波理论发展了动态应力应变关系的一致近似模型。这个模型可以对多胞材料在过渡模式及冲击模式下的动态压溃行为进行表征。