高应变率下延性金属中微孔洞贯通行为的数值分析

王永刚，徐东明

(宁波大学 力学与材料科学研究中心，浙江 宁波315211)

0 引言

延性金属材料在高应变率加载下发生拉伸断裂主要包括3 个物理过程:微孔洞的成核、长大和贯通［1］。迄今，人们对于微孔洞的成核和长大规律进行了大量的实验和理论分析工作［2-4］，但微孔洞相互贯通行为的研究工作相对很少。已有学者研究指出微孔洞的贯通在动态拉伸断裂过程中起着很关键的作用［5］，因此，研究微孔洞贯通机理对揭示材料动态拉伸断裂的整个物理过程和本质有重要意义。

Curran 等［6］基于实验观测结果认为有两类微孔洞贯通机理，即微孔洞之间直接贯穿和微孔洞之间的变形局部化。事实上，微孔洞相互作用是通过它们周围的塑性应变场来实现的。当两个微孔洞之间距离达到一个临界值时，它们周围的塑性应变场发生交接，从而改变微孔洞的长大速率。一些学者曾对微孔洞之间开始贯通临界条件进行了理论研究。Needleman 和Tvergarrd［7］认为材料内部微孔洞体积分数f 达到一个临界值fc时微孔洞贯通开始。Brock和Sun［8］研究指出当材料中初始微孔洞体积分数较小时，临界微孔洞体积分数fc对应力三轴度和加载应变率不敏感。而Pardoen 和Hutchinson［9］认为fc依赖于微孔洞初始体积分数、加载应力三轴度及微孔洞形状。Brown 和Embury［10］提出了一个微孔洞开始贯通的临界韧带距离(ILDc)准则，指出当微孔洞间韧带长度等于微孔洞的直径时，微孔洞从独立长大开始向贯通转变。Horstemeyer 等［11］采用微力学有限元计算方法研究准静态加载条件下304L不绣钢和6061T6 铝中两个微孔洞的贯通行为，结果表明:ILDc的值依赖于边界条件和材料的特性，约为2～8 倍微孔洞直径。最近，Seppälä 等［12-13］采用分子动力学方法对金属铜动态拉伸断裂的微孔洞长大和聚集行为进行了数值模拟，计算结果显示:当ILDc等于0.5 倍微孔洞直径时，微孔洞周边塑性形变区相遇，微孔洞长大的速率突然增加，最终导致了微孔洞之间发生了直接贯穿。

无论是采用fc还是ILDc作为微孔洞贯通开始的准则，人们对微孔洞贯通行为的认识还是比较初步的，存在一些争议。本文拟采用动力学有限元计算模型，对高应变率加载条件下纯铝材料中双孔洞贯通行为进行数值分析，重点讨论初始孔洞间韧带距离和加载应变速率对双孔洞贯通行为的影响规律，为合理建立孔洞贯通的理论模型奠定一定的基础，从而进一步揭示延性金属拉伸型断裂的物理机制。

1 计算模型建立

高应变率拉伸载荷下延性金属中出现的损伤主要是大量的随机分布的球形微孔洞。因此，在一定假设条件下，可以把含损伤材料离散成大量含孔洞的圆柱体胞的集合体。这里，为了研究孔洞之间的贯通行为，在数值模拟中设定每一个圆柱体胞含有两个大小相同的球形孔洞，二分之一模型见图1(a)所示。在对称载荷条件下，我们采用轴对称的二维单元来进行计算分析，有限元模型见图1(b)所示。微孔洞的直径d0为50 μm，圆柱体胞的直径D 和高度H 相等，均取为1 mm.圆柱体胞的直径和高度是球形微孔洞直径的20 倍。

图1 二分之一体胞模型及轴对称有限元分析模型Fig.1 Half of the axisymmetric unit cell and FEM model

计算中，边界条件采用速度边界，如图1所示，vz、vr分别表示在轴向和径向施加的速度函数。速度从0 开始经历0.1 μs 线性增加到恒值，随后保持不变。选取的模拟材料为纯铝，在高应变率加载下，采用考虑塑性硬化的流体弹塑性本构关系和Grüneisen 状态方程进行数值分析。

式中:σij为应力张量;p 为压力;sij为偏应力张量。

式中:σy为屈服应力考虑塑性硬化时，

式中:σ0为初始屈服应力;为等效塑性应变;Eh塑性硬化模量。

在材料处于压缩状态下，Grüneisen 状态方程的压力计算:

当材料处于拉伸状态下，

计算参数见表1.

表1 纯铝的本构和状态方程参数Tab.1 Constitutive and equation of state parameters for pure aluminium

2 计算结果讨论

在轴向速度场和径向速度场相同条件下，图2显示轴对称体胞模型边界单元上宏观轴向应力和径向应力几乎完全相等，应力开始时基本呈线性增加，随后由于微孔洞长大的影响，出现非线性增长，很快达到极值点，随后快速下降。由此可见，体胞模型中微孔洞承受的是一个准静水拉伸载荷条件，这与一维应变层裂实验中样品中微孔洞所承受了应力状态基本一致。因此，后续的有关微孔洞之间贯通行为的研究结果可以直接应用于层裂实验的理论分析和数值模拟中。另外，在准静水拉伸载荷下，图3显示不同时刻两个微孔洞周围的塑性应变场的演化云图。图3(a)显示两个微孔洞各自形成几乎均匀的球对称塑性应变场，微孔洞呈现相互独立球形快速长大，而随着微孔洞不断长大，两个微孔洞之间韧带区域塑性应变场相互交叠，如图3(b)所示，当此区域的应变场达到一个临界值时，两个微孔洞即发生贯通行为，直至直接贯穿，如图3(c)和图3(d)所示。

图3给出了微孔洞贯通行为的一些定性描述，为了定量讨论微孔洞之间的贯通行为，这里采用模拟计算中实时跟踪微孔洞长大过程中形状改变来判断，即通过比较图1中标识的微孔洞上顶点A 和下顶点B 的轴向位移。如果微孔洞处于独立长大阶段，A 点和B 点轴向位移是一致的，如果A 和B 点轴向位移发生明显的改变，表明微孔洞之间韧带区域塑性应变场交叠改变了B 点的轴向位移增长速率，此时刻即可判断微孔洞贯通开始，所对应的图1中标识的B 点和C 点之间距离即为微孔洞开始贯通的临界韧带距离(ILDc).下面首先讨论初始孔洞间韧带距离ILD0对微孔洞贯通动力学行为的影响。

图2 宏观轴向应力和径向应力的比较Fig.2 Comparison of global axial stress and radial stress

图3 不同时刻微孔洞周围等效塑性应变云图Fig.3 Effective plastic strain contours around void

2.1 ILD0的影响

在加载应变率104s-1条件下，考察不同初始微孔洞间韧带距离对孔洞贯通动力学行为的影响，选择ILD0/d0为2、3、4 和5，即ILD0分别为100 μm、150 μm、200 μm、250 μm(上面图2和图3的数据对应的ILD0/d0=3).图4(a)～图4(d)左边纵轴给出了不同ILD0/d0条件下实时的微孔洞间韧带距离ILD 与实时微孔洞直径d 之比随时间演化过程，而其右边纵轴同时给出了A 点和B 点的轴向位移时程。这里需要说明一下，为什么没有直接来讨论临界ILD，而是讨论它与实时微孔洞直径d 之比，这主要是因为微孔洞贯通是通过微孔洞周围的塑性应变场来实现的，而微孔洞的直径直接影响着其周围的塑性应变场。对比A 点和B 点的轴向位移时程，发现:微孔洞贯通开始前A 点和B 点轴向位移是相同的，两个微孔洞之间没有发生相互作用，微孔洞处于独立长大阶段，随着微孔洞不断长大，它们周围的塑性应变场发生相互影响，导致随后A 点和B 点轴向位移明显发生偏离，利用此点恰好可以判断微孔洞开始贯通时刻，因此，这里的计算结果也检验了上面提出的通过微孔洞形状改变来判断微孔洞贯通开始的有效性。而从ILD/d 的时程曲线上来看，在0.3 μs 之前，ILD/d 减小速率缓慢，此阶段微孔洞基本还处于弹性变形阶段，随后进入微孔洞独立快速长大阶段，这时ILD/d 迅速减小，紧接着微孔洞间发生贯通行为，导致ILD/d 减小速率再次变缓，最后慢慢趋于0，即微孔洞发生了直接贯穿，如上面图3(d)所示。对比图4(a)～图4(d)，由A 点和B 点轴向位移明显发生偏离点容易确定微孔洞间贯通开始的临界韧带距离，在图中用黑点标识，结果显示:微孔洞从独立长大向微孔洞贯通转变开始的临界韧带距离约为0.5 倍(ILDc=0.5 ±0.1d)的微孔洞直径，并且它对初始的ILD0/d0不敏感。

临界韧带距离作为微孔洞贯通开始的准则，比较直观，但相对缺乏物理本质的认识和提炼。分子动力学研究成果［12-13］已经在微观层次上揭示微孔洞间韧带区域位错密度迅速增长导致微孔洞发生贯通。位错运动宏观上对应于塑性应变，因此为了进一步理解微孔洞贯通动力学行为，图5给出了不同ILD0/d0条件下微孔洞间韧带距离中点处典型单元的等效塑性应变时程曲线。利用微孔洞开始贯通的时间对应条件，可以在图5中找出微孔洞贯通时刻所对应的塑性应变，图中用黑点标识，结果显示:临界塑性应变阈值约为0.5，随着ILD0/d0增大，临界塑性应变阈值缓慢减小。

2.2 加载应变速率的影响

通过改变边界速度函数，可以调整加载应变速率。这里以ILD0/d0=3 模型为例，讨论加载应变速率对微孔洞贯通行为的影响。图6给出了不同加载应变速率下ILD/d 时程曲线和A 和B 点轴向位移时程曲线，结果显示:随着加载应变速率增大，微孔洞长大速率提高，微孔洞间贯通时间提前，但对微孔洞间开始贯通临界韧带距离没有影响，还保持在0.5 倍微孔洞直径。图7给出了不同加载应变速率条件下微孔洞间韧带距离中点处典型单元的等效塑性应变时程曲线，结果也同样显示微孔洞开始贯通的临界塑性应变阈值为0.5，图中黑点标识。图8给出了不同应变速率条件下体胞模型边界单元上的宏观平均应力随时间的变化曲线，从中可见:随着加载应变速率增大，在微孔洞开始贯通时宏观平均应力(黑点标识处)在增大，不对应于应力峰值，也就是说微孔洞间贯通开始后，宏观平均应力场还在增大，材料仍然可以继续承载。以前有很多学者为了回避对微孔洞贯通过程的理论描述，常常假定微孔洞贯通开始后材料迅速失去承载能力，这个观点是值得商榷的。以后在建立合理的微孔洞贯通理论模型时必须考虑微孔洞贯通也是需要一个过程的，它不是瞬时事件。

图4 ILD/d 演化时程曲线及A 点和B 点轴向位移时程曲线Fig.4 Evolution of ILD/d profiles and A and B points’axial displacement profiles

图5 微孔洞间韧带中心处典型单元的等效塑性应变时程Fig.5 Effective plastic strain profiles of representative element in the center of intervoid ligament distance

图6 ILD/d 时程曲线和A 和B 点轴向位移时程曲线Fig.6 Evolution of ILD/d profiles and A and B points’axial displacement profiles

图7 微孔洞间韧带中心处典型单元的等效塑性应变时程Fig.7 Effective plastic strain profiles of representative element in the center of intervoid ligament distance

图8 微孔洞间韧带中心处典型单元的等效塑性应变时程Fig.8 Effective plastic strain profiles of representative element in the center of intervoid ligament distance

3 关于ILDc和fc的讨论

目前，关于微孔洞开始贯通的准则主要有两种观点，有些学者采用微孔洞间临界韧带距离ILDc，还有些学者常采用临界微孔洞体积分数fc.如果基于微孔洞在材料中统计均匀分布的假定，设用边长为L 的立方体来代表受损伤材料的一个特征体积，假设有8 个直径为d 的微孔洞均匀分布于立方体的8 个角上，如图9所示。每个微孔洞只有八分之一属于此立方体，因此微孔洞开始贯通时的临界微孔洞体积分数:

由此可见，ILDc和fc是可以相互转换的，在物理本质也是一致的。这里近似假定图9中微孔洞两两之间的相互贯通行为是相互独立的，当然这个假定具有一定的合理性，因为微孔洞之间相互贯通是通过塑性应变场交叠来实现的，图9中其他6 个微孔洞由于距离相对比较远，对两个微孔洞之间的塑性应变场的影响相对较弱。因此，在微孔洞两两之间的相互贯通行为是相互独立的假定下，可以通过含两个微孔洞圆柱体胞组合排列形成空间的8 个微孔洞，同时体胞模型中微孔洞间临界韧带距离ILDc在物理本质上可以近似替代含8 个微孔洞的立方体损伤单元中ILDc.对于纯铝，上述数值模拟得到ILDc为0.5 倍微孔洞直径，通过(7)式，计算得到的fc约为0.15.文献［14］曾对纯铝开展过系列层裂实验研究，基于对回收样品中含有的微孔洞宏观统计分析，获得了微孔洞开始贯通的微孔洞体积分数fc约为0.11，与数值分析结果非常接近，这表明通过数值分析获得微孔洞开始贯通临界韧带距离或临界孔洞体积分数具有一定的合理性。

图9 临界孔洞间韧带距离和临界孔洞体积分数的几何说明Fig.9 Geometry relationship between ILDc and fc

4 结论

采用LS-DYNA 动力学有限元程序，对轴对称体胞模型中双孔洞的贯通行为进行了数值分析，讨论了初始孔洞间韧带距离和加载应变速率的影响，得到了如下几点认识:

1)提出了一个微孔洞开始贯通的判断标准，即通过比较微孔洞上顶点和下顶点的轴向位移曲线的明显偏离来判断微孔洞开始贯通，数值结果显示该方法具有一定的合理性和有效性。

2)微孔洞从独立长大向微孔洞贯通转变开始的临界韧带距离约为0.5 倍微孔洞实时直径，并且不依赖于初始微孔洞间韧带距离和加载应变速率。

3)微孔洞间贯通行为是通过塑性应变场相互作用实现的，数值计算结果显示:在孔洞间韧带区域的塑性应变场需要达到一个临界阈值才能驱动微孔洞开始贯通，该临界塑性应变阈值基本保持在0.5左右。

4)微孔洞在材料中统计均匀分布和微孔洞之间的相互贯通行为是相互独立的假定下，临界孔洞间韧带距离ILDc判据和临界孔洞体积分数判据在物理本质也是一致的，可以相互转换。

References)

［1］ Antoun T，Seaman L，Curran D R，et al.Spall fracture［M］.Berlin:Springer，2003.

［2］ Johnson J N.Dynamic fracture and spallation in ductile solids［J］.Journal of Applied Physics，1981，52:2812 -2825.

［3］ Molinari A，Wright T W.A physical model for nucleation and early growth of voids in ductile materials under dynamic loading［J］.Journal of the Mechanics and Physics of Solids，2005，53:1476-1504.

［4］ Chen D，Tan H，Yu Y，et al.A void coalescence-based spall model［J］.International Journal of Impact Engineering，2006，32(11):1752 -1767.

［5］ Czarnota C，Mercier S，Molinari A.Modeling of nucleation and void growth in dynamic pressure loading，application to spall test in tantalum［J］.International Journal of Fracture，2006，141(1 -2):177 -194.

［6］ Curran D R，Seaman L，Shockey D A.Dynamic failure of solids［J］.Physics Reports，1987，147:253 -388.

［7］ Tvergarrd V，Needleman A.Analysis of the cup-cone fracture in a round tensile bar［J］.Acta Metal.，1984，27:1 -15.

［8］ Brock W，Sun D Z，Honign A.Verification of the transferability of micromechanical parameters by cell model calculations with visco-plastic materials［J］.International Journal Plasticity，1995，11(8):971 -989.

［9］ Pardoen T，Hutchinson J W.An extended model for void growth and coalescence［J］.Journal of the Mechanics and Physics of Solids，2000，48:2467 -2512.

［10］ Brown L M，Embury J D.The initiation and growth of void at second phase particles［C］∥Proceedings 3rd International Conference on Strength of Metals and Alloys.London:Institute of Metals，1973.

［11］ Horstemeyer M F，Matalanis M M，Sieber A M，et al.Micromechanical finite element calculations of temperature and void configuration effects on void growth and coalescence［J］.International Journal Plasticity，2000，16:979 -1015.

［12］ Seppälä E T，Belak J，Rudd R E.Onset of void coalescence during dynamic fracture of ductile metals［J］.Physics Review Letters，2004，93(24):245503.

［13］ Seppälä E T，Belak J，Rudd R E.Three-dimensional molecular dynamic simulation of void coalescence during dynamic fracture of ductile metals［J］.Physics Review B，2005，71:064112.

［14］ Qi M L，He H L.Simulation of critical behavior on damage evolution［J］.Chinese Physics B，2010，19(3):036201-5.