多胞材料的动态应力应变状态及其一致近似关系*

王 鹏，朱长锋，郑志军，虞吉林

(1.中国科学技术大学近代力学系中国科学院材料力学行为和设计重点实验室,安徽 合肥 230026;2.中国工程物理研究院总体工程研究所,四川 绵阳 621999)

多胞材料在高速冲击下具有变形局部化的典型特征。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个胞元组成。试件相对密度ρ0/ρs= 0.1，其中ρ0为蜂窝试样变形前的初始密度，ρs为基体材料的密度。蜂窝试件的平均胞元尺寸约为2.5 mm，其值定义为与平均胞元面积相等的圆的直径。

图1 Voronoi蜂窝试件细观有限元模型Fig.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/m3。所有接触面均定义为通用接触，摩擦因数取为0.02。为了模拟面内变形，模型中所有节点的面外变形均被限制。

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

1.2 截面应力计算方法

(1)

(2)

2 结 果

2.1 一维应力分布

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

图2 两种冲击速度下的应力分布图Fig.2 One-dimensional stress distributions at two impact velocities

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

图3 一维的应力分布及其应力梯度分布Fig.3 One-dimensional stress distributions and the corresponding stress gradients

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

图4 不同冲击速度下冲击波位置随时间的关系Fig.4 Variation of shock front position with impact time at different impact velocities

2.3 一维冲击波理论

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

vB(t)-vA(t)=vs(t)(εB(t)-εA(t))

(3)

σB(t)-σA(t)=ρ0vs(t)(vB(t)-vA(t))

(4)

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

(5)

图5 一维应力和应变分布及其理想化Fig.5 One-dimensional stress and strain distributions and the idealizations

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

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

2.4 冲击波速度的比较

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

vs=v/εL

(6)

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

(7)

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

图6 不同方法得到的冲击波速度的比较Fig.6 Comparison of shock wave speeds obtained using different methods

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

3 分析和讨论

3.1 冲击波速度的一致近似

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

(8)

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

(9)

图7 冲击波速度与冲击速度的关系Fig.7 Variations of the shock wave speed with the impact velocity

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

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

(10)

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

(11)

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

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

(12)

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

(13)

图8 波后应变和波后应力随冲击速度的变化Fig.8 Variations of strain and stress behind shock front with impact velocity

(14)

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

图9 不同冲击速度下的速度分布Fig.9 Velocity distributions at different impact velocities

图1 0 准静态及动态应力应变关系Fig.10 Quasi-static and dynamic stress-strain relations

4 结 论

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