Bubble Dynamics Characteristics Near Double-Layer Cylindrical Structures with a Hole Based on PIV Experiments

颗粒金属材料具有多尺度、孔隙度、流动性及一定强度等特性，已被广泛应用于许多重要工业领域，例如航空航天、自动化、电子^[1–2]。颗粒金属经过冲击压缩可以制备出不同性能的材料，因此金属粉末压实越来越受到工业制造领域的关注^[3]。相较于传统的压实方法，冲击压缩可以将粉末快速熔化、快速固化，压缩得更加密实^[4–5]。国内外学者从20世纪80年代开始对颗粒金属材料的冲击压缩行为进行研究，但其冲击压缩过程的力学响应尚未得到很好的理解，而这恰恰是颗粒金属材料设计和应用的关键问题，尤为必要进行深入研究。

颗粒金属材料的冲击压缩过程历时短、反应迅速，且粒子直径很小，准确捕捉其响应特征较为困难^[6–8]。一些学者以单质金属粉末、金属粉末混合物作为研究对象，借助温度压力传感器等瞬态测试工具，通过高速飞片冲击试验对材料冲击前、后的初态和终态进行测定^[9–10]。然而，实验测试难以获得压缩过程的响应信息。近来，大量数值模拟方法从细观角度出发研究颗粒材料的冲击压缩行为，如欧拉流体动力学方法、有限元法、离散元法，主要研究的不均匀颗粒材料包括铝^[11]、铜^[12]、多晶铁^[13]、铍^[14]等。连续介质有限元方法可以表征局部相对密度分布、应力分布和粉末流动，但由于粉末是不均匀的，不能从粒子尺度表征颗粒接触和相互作用。离散元方法能弥补上述缺点，可从颗粒水平观察力学响应，但由于简化了粒子接触算法，所以不适用于压实过程。多颗粒有限元方法^[15–17]集成以上两种方法的优点，适用于完全压实过程且能很好地阐明颗粒局部变形及力链传递。

已有的一些理论模型可以用来描述颗粒金属材料的动态压缩行为。压力-孔隙比(P-$\alpha$)模型^[18]仅需要定义参数$\alpha$，即孔隙材料与基体材料体积之比，就可以捕捉到压力随孔隙比的变化关系。P-$\alpha$模型通过引入$\alpha$的方式成功分离孔穴压溃和基体压缩的影响，但是$\alpha$的演化关系需要额外补充。经验性的球壳压溃模型^[6]随后被发展，该模型提出包含惯性项的空穴压溃方程，成功应用于静水压加载，压力为准静态和惯性相关项之和。基于球壳压溃模型，许多学者对其进行修正，分别引入偏应力^[7]、温度^[19]、黏性^[20]、微惯性^[21]的影响。然而，大部分模型都是基于经验提出的，直接得到动态的应力-应变关系在实验上是困难的，而细观有限元模型可以得到局部应力-应变信息。因此，本研究采用细观有限元模拟冲击压缩过程，研究颗粒金属材料的冲击压缩模型。

本研究采用多颗粒有限元方法，从粒子层次建立颗粒金属材料的二维细观有限元模型，研究在恒速冲击载荷作用下铜粉末的力学行为，分析冲击波的传播特性，观察变形及速度分布特征，表征冲击速度与冲击波速度的关系及波后动态应力-应变状态。

1.   数值模拟

1.1   几何模型

模型简化为二维等径圆形颗粒，颗粒直径为30 ${\text{μ}}$m，初始颗粒堆积方式大致分为4类：随机堆积、简单立方堆积、密排六方堆积和蜂窝堆积^[22–24]。其中密排六方堆积及随机堆积的准静态压缩测试结果表明，压力-平均相对密度关系的变化趋势与实验结果吻合较好，随机堆积可以较好地描述颗粒间的摩擦作用机制。为此，采用随机堆积的建模方式，颗粒在空间随机分布，颗粒之间不重叠。

采用Matlab软件编写颗粒堆积模型，建立细观分布模型的步骤为：（1）计算出在一定空间内设定孔隙率$\phi$（孔隙部分所占比重）下等径颗粒的个数，如表1所示；（2）随机撒点的坐标作为颗粒核点的初始位置，通过调整核点位置确保任意两点之间的距离大于或等于2倍颗粒半径；（3）实时更新颗粒位置坐标，直到所有颗粒都不重叠为止。例如图1颗粒材料生成试件的孔隙率分别为0.50、0.35、0.25。

图  1  颗粒金属材料的几何模型

Figure  1.  Geometric models of granular materials

下载: 全尺寸图片 幻灯片

表  1  不同孔隙率下细观模型所含颗粒数

Table  1.  Particle numbers of mesoscale model at different porosities

Porosity ($\phi$)	Number of particles (N)
0.60	142
0.50	177
0.45	195
0.40	213
0.35	230
0.30	248
0.28	255
0.25	265

下载: 导出CSV 

| 显示表格

1.2   细观有限元模型

数值模拟采用ABAQUS/Explicit 6.11，基体材料选择实验数据丰富的铜颗粒。对每个颗粒单独划分网格后导入有限元软件，颗粒网格划分如图2(b)所示，边界处网格进行加密处理，保证内部颗粒相互接触作用的精度，内部采用粗网格能有效提高计算效率，以节约计算时间。经过网格收敛性分析，每个颗粒划分成100个六面体单元，包含226个节点，采用体单元可以增加计算精度，但需约束节点在厚度方向的位移，限制粒子在二维平面内的变形。

图  2  铜粉末恒速压缩的细观有限元模型(a)以及颗粒初始网格划分(b)

Figure  2.  A finite element model of copper powders under constant-velocity compression (a) and the initial mesh of a particle (b)

下载: 全尺寸图片 幻灯片

铜粉末在高速冲击加载时表现出应变率硬化、应变硬化、温度软化等特性，在此采用典型的应变率、温度相关的Johnson-Cook模型表征^[25]，写作

式中：屈服应力$\sigma$₀为90 MPa，应变硬化系数n为0.31，材料参数B为292 MPa，C为0.025，参考应变率${\dot \varepsilon _0}$取为1 s^–1，熔点温度T_m为1360 K，参考温度T_r为300 K，温度软化系数m为1.09。

本研究采用的基体材料模型可以较好地表征真实铜粉末在冲击加载时的温度响应特征，粉末在压缩变形后，部分塑性功转化为热能，转化系数为0.9，从而带来温度升高的现象^[26]。

颗粒试件位于长和宽均为0.5 mm的正方形容器内，上端、下端及支撑端刚性板均固定不动，冲击端刚性板以速度V向右恒速压缩，如图2(a)所示。颗粒之间的摩擦系数设为0.2，颗粒与刚性板无摩擦^[15]。颗粒之间的所有接触面定义为通用接触。法向接触属性采用硬接触，切向选用罚函数接触。

2.   数值结果

2.1   变形特征

通过对初始孔隙率为0.28的样本进行冲击压缩模拟，观察发现在动态冲击下，颗粒金属材料呈现出高度局部化的变形特征。随着冲击速度的提高，变形局部化的特征更加明显，波阵面变得更加平整，波阵面后方的颗粒也随着冲击速度的增大而变形更加急剧，由均匀变形渐渐转变为不均匀变形。应力波沿加载方向以近似平面的形式逐层传播，波阵面宽度大约为2～4个颗粒直径。当冲击速度为95 m/s时，波阵面后方还未完全压实，仍存在一些孔隙，孔隙率的减小依赖于颗粒的平动、转动和塑性变形。当速度增加到一定值时，波后几乎完全压实，没有孔隙，此时波后应变为最大压实应变，如图3所示。

图  3  不同时刻试件变形图

Figure  3.  Deformation patterns at different times

下载: 全尺寸图片 幻灯片

2.2   一维速度场分布

单轴冲击压缩下，应力波以近似一维应变平面波的形式从冲击端向支撑端传播，因此颗粒材料的冲击压缩问题可以简化为一维问题进行分析。引入一维速度场分布定量化地表征冲击波的传播特征。从有限元结果可以直接得到在Lagrange坐标系下每个节点的速度信息，将各截面上的节点在加载方向的速度分量进行平均，得到铜粉末在不同时刻t的一维速度场分布。速度场结果很好地捕捉到高速冲击下铜粉末呈冲击波模式的变形特征，表明速度急剧变化集中在一条狭窄的区域内。冲击速度越高，此区域越窄；而冲击速度较低时，波阵面越宽，如图4所示。

图  4  不同冲击速度下试件内的一维速度(v)分布

Figure  4.  One-dimensional velocity distribution in specimens under different impact velocities

下载: 全尺寸图片 幻灯片

2.3   冲击波速度

塑性冲击波波阵面的位置可由局部应变率最大值的位置确定。基于一维应力波理论，连续性方程为

可知速度梯度与应变率相等。一维速度场分布特征表明，高速冲击时，铜粉末的速度分布可以分为3个区域：波前的未变形区域、波后的压实区和状态发生突变的区域。依据速度梯度的变化趋势可以看出速度的变化程度，速度梯度最大处，速度变化程度最剧烈，将该位置作为此时刻冲击波波阵面的位置，如图5所示。

图  5  一维速度分布及其速度梯度 (V=225 m/s, t=0.4 μs, $\phi$=0.28)

Figure  5.  One-dimensional velocity distribution and the corresponding velocity gradient (V=225 m/s, t=0.4 μs, $\phi$=0.28)

下载: 全尺寸图片 幻灯片

不同时刻的波阵面位置与时间近似呈线性关系，此直线的斜率即为冲击波波速。冲击速度越高，斜率越大，冲击波波速越大，如图6所示。

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

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

下载: 全尺寸图片 幻灯片

2.4   与实验结果的比较

将多颗粒细观有限元模拟所得结果与文献中气枪冲击实验结果^[27]及欧拉有限元数值模拟结果^[28]进行比较，如图7所示。Benson^[28]所采用的欧拉有限元方法适用于速度更高的情形，而在较低的冲击速度下，该方法对波阵面位置的捕捉不准确，导致冲击波波速的结果偏高。而采用多颗粒有限元方法所得结果与实验结果吻合较好，二维简化模型所得冲击波波速略低于实验结果。Borg等^[29]通过对陶瓷颗粒进行二维简化，数值模拟分析其冲击Hugoniot关系，所得结果类似。原因是相同孔隙率下，二维模型的堆积配位数比三维少，较真实粉末有更多的自由空间，粒子连接链较短，在一定程度上妨碍冲击波的传播，从而导致冲击波波速较实验结果偏低。此外，在较低的冲击速度下，颗粒材料存在孔隙和分布不均匀性，波阵面不平整，选取波阵面位置和计算波速时会带来一定的误差。随着冲击速度的增大，波阵面越来越狭窄，位置较易确定，波速计算精度提高，误差逐渐减小。

图  7  冲击波速度-粒子速度关系与已有文献结果对比

Figure  7.  Comparison of the shock wave speed-particle velocity relations with the data in the literature

下载: 全尺寸图片 幻灯片

3.   理论分析

3.1   塑性冲击波模型

由第2节可知，在较高速度冲击下，应力波在颗粒试件中近似以塑性冲击波的形式沿加载方向逐层传播，基于应力波理论^[30]分析，冲击波波阵面上的质量守恒和动量守恒方程分别为

式中：V_s为冲击波速度，{v_A, $\varepsilon$_A, $\sigma$_A}、{v_B, $\varepsilon$_B, $\sigma$_B}分别为波阵面前后的速度、应变、应力状态。

波前颗粒材料由于处于较疏松的状态，颗粒间存在间隙，弹性波不能向前传播，因此应力、应变、速度都假设为零，即波前状态为v_A=0、$\varepsilon$_A=0、$\sigma$_A=0，与冲击速度无关。

在动态冲击载荷作用下，波后应变与冲击速度是有关的，会随着冲击速度的增大而逐渐增大，直到达到理论上的最大压实应变。当冲击速度增大到一定程度时，波后应变逐渐逼近理论上能达到的最大压实应变，写作

式中：$\varepsilon$_m=$\phi$，为最大压实应变。

恒速加载时波后速度为v_B=V，因此可求得冲击波波速

将（6）式代入（4）式可得波后应力

观察发现，在速度较高时波前状态为v_A=0、$\varepsilon$_A=0、$\sigma$_A=0，波后状态为v_B=V、$\varepsilon$_B=$\phi$、$\sigma$_B=–$\rho$₀V²/$\phi$，仅含有一个材料参数初始孔隙率$\phi$。基于此特征，发展以动态锁定应变为唯一参数的冲击波模型，如图8所示。图8中箭头表示从波前状态沿Rayleigh线突跃到波后状态。该模型未考虑材料的率敏感性和温度的影响，仅以初始孔隙率$\phi$表征较高冲击速度下铜粉末的动态应力-应变关系。

图  8  动态应变锁定冲击波模型

Figure  8.  Dynamic shock model with strain locked

下载: 全尺寸图片 幻灯片

3.2   理论模型的应用

通过已有文献实验结果验证了多颗粒有限元方法数值模拟的有效性，继续通过此方法计算不同孔隙率的铜粉末在不同冲击速度下粒子速度与冲击波速的关系。孔隙率分别为0.25、0.28、0.30、0.35、0.40、0.45、0.50、0.60的V_s-v关系，如图9所示。而通过动态应变锁定的冲击波模型得出的V_s-v关系与细观有限元模拟结果吻合极好，由此证明动态应变锁定的冲击波模型可以很好地表征不同孔隙率铜粉末的Hugoniot关系。

图  9  不同孔隙率的铜粉末粒子速度与冲击波速度的关系

Figure  9.  Variation of the shock wave speed with the particle velocity at different porosity of copper powders

下载: 全尺寸图片 幻灯片

不同冲击速度下，冲击载荷不同；冲击速度越大，冲击端应力越大。高速恒速加载情形下，冲击端应力近似等于波后应力，由细观有限元模拟提取出的冲击端应力如图10所示。颗粒材料的不均匀性以及孔隙的存在导致冲击力振荡较大，但基本维持在某一特定值上下。由动态应变锁定的冲击波模型计算得出的波后应力可以很好地捕捉到这一水平。

图  10  不同冲击速度下冲击端应力与冲击波模型所得波后应力比较

Figure  10.  Comparison of the stress at the impact end obtained from finite element models and the stress behind the shock wave front predicted by the shock model at different impact velocities

下载: 全尺寸图片 幻灯片

对于不同孔隙率的颗粒材料，波后应力随着孔隙率的增大逐渐减小，在某一恒定速度加载情形下，对塑性冲击波传到支撑端前这一时间范围内的冲击端应力进行平均，将所得应力与通过理论模型计算得出的波后应力进行比较。结果表明，冲击端平均应力与冲击波模型所得波后应力相对误差较小（小于6%），如图11所示。相对误差由下式求得

图  11  不同孔隙率下冲击端应力与冲击波模型所得波后应力比较

Figure  11.  Comparison of the stress at the impact end and the stress behind the shock wave front obtained by the shock model under different porosities

下载: 全尺寸图片 幻灯片

式中：$\sigma$_avg为冲击端平均应力，$\sigma$_s为理论模型计算所得波后应力。

4.   结　论

采用多颗粒有限元方法建立二维颗粒金属材料的冲击压缩模型。基于不同孔隙率下等径圆形颗粒的随机堆积模型，采用恒速压缩的加载方式，得到铜粉末的一维速度场分布，确定了塑性冲击波波阵面位置，进而得到了粒子速度与冲击波波速的关系。波后应变随着冲击速度增大，最终保持为最大压实应变不变，而波前应力、应变、速度均为零。基于颗粒材料的这一特征，发展了以动态锁定应变为唯一参数的冲击波模型，该模型可以很好地表征V_s-v关系并确定波后应力振荡的平均值。同时需要指出的是，该冲击波模型仅在速度较大、波后应变为最大压实应变时有较好的适用性，而在速度较小时会有较大的误差，在中低速下还需要发展更加精确的模型。颗粒金属材料的细观特性对其宏观力学性能有重要影响，本研究虽然从细观尺度建立模型，可以分析其宏观和微观力学响应，但是并未研究细观参数颗粒形状、尺寸、分布等的影响。基体材料选用的是Johnson-Cook模型，考虑了基体材料的应变率效应，但是在高速冲击加载时，应变率效应与惯性效应难以解耦，理论模型并未考虑应变率效应，即使这样仍可以较好地捕捉其主要特征。