应力波与缺陷相互作用的宏观微观数值模拟＊

郭昭亮，任国武，汤铁钢，刘仓理

（中国工程物理研究院流体物理研究所冲击波物理与爆轰物理重点实验室，四川 绵阳621999）

在材料断裂机理研究以及工程应用中，裂纹诱致的断裂行为，一直都是力学家和工程师关注的焦点。尤其对金属及其合金等延性材料，裂纹尖端或刻槽尖端附近的塑性变形区，在裂纹的起裂、扩展、止裂过程中起着重要的作用［1－4］。对塑性区形状大小的估计历来已久，传统的做法从宏观断裂力学角度按照各种屈服准则确定其形状。G.R.Irwin估计了理想裂纹尖端塑性区的大小、形状，并作为弹塑性等效裂纹长度的修正［5］。张亚等［6］针对复合型裂纹，基于小范围屈服准则，采用俞茂宏统一强度理论，给出了裂尖塑性区的统一解析解。X.Gao等［4］基于小范围屈服条件，采用Hill屈服准则计算了不同加载情况下的裂尖塑性区状，发现塑性区呈现蝴蝶状，且弹塑性区边界光滑。Y.Huang等［7］利用最大裂纹张开位移估计了塑性区尺寸大小。然而对于真实材料，由于不存在理想的数学裂纹，则可以使用一定几何形状的缺陷近似裂纹，缺陷局域的塑性区大小、形状可以通过计算缺陷处的应力集中，并依据强度理论进行估计。S.Q.Shi等［8］使用简单的模型计算了刻槽尖端的最大正应力与塑性区的大小。张培源等［9］探讨了裂纹尖端为圆弧形的钝化模型，并依据Tresca屈服条件，得到了钝化裂纹前缘塑性区的线场分析解，对于含径向裂纹和圆弧钝化区的圆盘外围，得到其线弹性解，形成了一套估计裂纹前端塑性区尺寸的方法。现有的研究指出，对于平面应变与平面应力情况，塑性区的大小形状以及尖端的正应力分布不同［10］。

在处理裂纹局域的塑性区时，理论大都基于准静态假定，并且认为塑性区首先出现在裂尖（或缺陷）顶端的边界。钱才富等［11－12］从微观断裂力学的角度更精细地计算了裂纹前缘塑性区和无位错区，发现应力最大的位置出现在远离裂尖一定距离处，并指出它们的存在及性质决定着裂纹的扩展行为。

然而，对于由动态载荷诱致的缺陷局域的塑性区形成与演化，以及由此诱发的裂纹萌生、扩展等现象的研究相对较少，研究主要集中在孔洞的长大与演化过程。R.E.Rudd等［13－14］、E.T.Seppälä等［15］采用分子动力学程序，模拟了韧性金属中的孔洞成核、长大，以及孔洞之间的相互作用、贯通。祁美兰等［16］、王永刚等［17］采用二维LS－DYNA程序，研究了平板撞击加载下含初始杂质的纯铝样品中微孔洞的成核、长大与闭合现象。在快速载荷作用下，由于应力波与材料缺陷的相互作用，将出现与准静态加载时完全不同的物理现象。本文中，利用LS－DYNA 3D有限元程序以及分子动力学方法，从宏观微观两个层次分别模拟含有预置圆形（椭圆形）缺陷的薄板，在动态拉伸载荷作用下，塑性区的形成、演化过程以及随之而来的裂纹动态扩展过程。

1 计算模型

1.1 宏观模型

宏观尺度上，采用LS－DYNA有限元程序，模拟含有预置缺陷的薄板（见图1）在动态拉伸载荷（见图2）作用下，缺陷与应力波相互作用对塑性区形成演化的影响，以及由此诱发的裂纹萌生、裂纹动态扩展行为，并通过调节预置缺陷的长轴与短轴比例实现不同的缺陷预置。假定材料动态应力应变曲线如图3所示。

图1 建模示意图Fig.1 Scheme of model

图2 拉伸载荷Fig.2 History of load

图3 应力应变曲线Fig.3 Stress－strain curve

对于同种材料，初始载荷σ0及缺陷的形状位置，决定了薄板中塑性区的形成位置与变形失效行为。

（1）当σ0＜Y／2时，薄板加载处发射弹性波，在远离缺陷处的应力波相互作用不会形成超过屈服应力的区域，只有在缺陷边界处的反射会形成塑性区。

（2）当Y／2＜σ0＜Y时，薄板加载处发射弹性波，在远离缺陷处的应力波相互作用区由于应力超过屈服应力将形成带状塑性区，此外由于缺陷边界处的反射也会形成塑性区。

（3）当σ0＞Y时，加载边界首先发射弹性前驱波，其后跟随塑性波，塑性波扫过的区域全部进入塑性，弹性前驱波由于相互作用及在缺陷附近的反射也会在短期内形成塑性区，当塑性波扫过后，将整体进入塑性。

本文中的计算主要侧重于第2种载荷情况。薄板材料为20钢，密度ρ0＝7.85g／cm3，剪切模量G＝82GPa，弹性模量E＝210GPa，泊松比ν＝0.286［18］。采用流体弹塑性本构模型，Grüneisen状态方程［19］；屈服应力取490MPa，为准静态屈服应力的2倍，塑性硬化模量取2GPa，为弹性模量的约1%［20］。

1.2 微观模型

微观尺度上，采用分子动力学方法模拟椭圆微缺陷在拉伸下的动态裂纹扩展。通过数值求解牛顿方程，获得整个体系粒子的运动轨迹。采用Verlet速度算法，多粒子相互作用势为Lennard－Jones势，并采用多项式截断［21］：

式中：η、ξ分别是能量尺度和长度尺度，参数a2、a3及rmax可通过式（1）中两个函数在rspl处的连续、包括其一阶和二阶导数连续，来共同决定，分别是a2＝0.542 449 4、a3＝0.093 505 27、rspl＝1.244 455和rmax＝1.711 238。由此势函数，可知其能量最低点对应的平衡位置a0＝21／6ξ。

计算中，采用恒定位移加载方法，即在每个时间步，都对最外层上下边界的两个原子层分别施加一个固定的位移增量Δu，这样可避免边界复杂的应力状态，其中Δu＝0.000 6。同时为避免外边界带来的波反射，在垂直于裂纹传播方向考虑加入线性速度梯度，减少模拟的时间。整个模拟研究将限制于拉伸型裂纹，其加载为沿x方向。

依据维里定理，在原子尺度上对i原子，其局域应力状态可定义为：

式中：α和β是方向，V是原子的体积，对于二维系统而言，是原子所占据的面积。

2 LS－DYNA计算结果

2.1 含偏心圆孔缺陷非对称薄板塑性区演化过程

图4为含有偏心圆孔缺陷的薄板，薄板尺寸为50mm×75mm×1mm，缺陷中心位于（25mm，25mm），圆孔直径10mm。在薄板两侧边界同时施加如图1～2所示的拉伸型阶跃载荷400MPa，并持续到计算结束。算例中，所施加的载荷满足Y／2＜σ0＜Y。

图4 含圆孔缺陷非对称薄板网格划分Fig.4 Mesh of an asymmetry sheet with a circularity hole

关注物理量等效塑性应变，可以得到应力波与缺陷相互作用时，塑性区随时间的演化过程，如图5所示。在4.5μs时，由于左侧应力波与缺陷作用，而右侧应力波尚未到达缺陷边界，所以只形成2个小的塑性变形区。由此可见，缺陷边界处的塑性区形成与整体应力是否平衡无关，它们是应力波与缺陷相互作用、反射塑性波导致的结果。随着时间演化，在7.3μs时，左右两侧的应力波在空间远离缺陷的位置发生相互作用，形成塑性带，这是应力波在空间相互作用导致的结果。随后，此塑性带逐步演化，最终与圆孔附近的塑性区汇合，在空间形成整体的塑性区域。

图5 含偏心圆孔缺陷非对称薄板塑性区演化Fig.5 Plastic zone evolution in an asymmetry sheet with a circularity hole

2.2 含中心椭圆孔缺陷对称薄板塑性区演化及失效过程

图6为含有椭圆孔缺陷的薄板，薄板尺寸为50mm×50mm×1mm，椭圆孔长轴10mm，短轴4mm，中心位置在（25mm，25mm），加载形式如图1～2所示。

图6 含椭圆孔缺陷对称薄板网格划分Fig.6 Mesh of a symmetry sheet with an elliptic hole

当塑性应变累积到一定程度，局部将会出现失效。为此，在程序中设置了失效应变0.04。应力波与缺陷相互作用后，塑性区域随时间的演化过程以及由应变累积导致的局部失效行为，如图7所示。

在4.4μs时，椭圆边缘同时形成4个小的塑性变形区，随着应力波的相互作用，塑性区迅速形成，在7.0μs时，塑性应变最大值的位置已经上移到椭圆边界外。由此诱致的裂纹萌生首先出现在远离缺陷边界，然后裂纹与缺陷贯通并向外扩展。有趣的是，随着裂纹的扩展，由于卸载波的相互作用，在与加载平行的方向，将出现裂纹的萌生、扩展，最终形成十字架形的断裂形式。塑性区的形成演化、裂纹的萌生扩展都与应力波和缺陷的相互作用密切相关，它展现了动态载荷作用下独特的变形与失效模式。

图7 含中心椭圆孔缺陷对称薄板塑性区演化及失效过程Fig.7 Plastic zone evolution and failure process in a symmetry sheet with an elliptic hole

3 分子动力学计算结果

对于微观情况，所选用模拟体系的几何形状是一个包含椭圆裂纹的二维密堆积体系，其长度沿x方向是lx＝600a0，沿y方向是ly＝693a0，而椭圆裂纹尺寸是沿x方向的长轴为30ξ，沿y方向的短轴为10ξ。总模拟原子数目为约480000。模拟采用微正则系综（粒子数、体积和总能守恒），此系综适合研究动态裂纹的非平衡过程。初始的模拟体系温度趋近于零。在实际模拟中，将长度参数ξ、能量参数η和粒子质量m都设置为参考量，那么所有量的量纲均为1。由德拜频率，时间步长选择为Δt＝0.003t0，t0＝ξ－1（η／m）1／2，每隔1000Δt获得一幅图像。

计算结果如图8所示。在计算至99t0时，远离缺陷上下边缘的位置出现位错堆积，在随后的6t0内，局部空洞化，由此诱致裂纹的萌生。裂纹沿着纵向扩展，与缺陷边缘贯穿，最终导致整体的断裂失效。在微观模拟中，即使不考虑无位错区，含有预置缺陷的薄板在动态拉伸下应力最大的位置也将出现在远离缺陷边缘，这与宏观模拟（见图7）类似。分子动力学模拟中，缺陷横向边缘附近已经存在了相当数量的位错，并出现了应变局域化现象。但由于系统提供的能量不足以使这些位错发展为裂纹，因而没有得到图7中所示十字架形的断裂形式。

图8 含有椭圆缺陷的二维密堆积体系的分子动力学模拟Fig.8 Molecular dynamic simulation in 2Dclosed－pack system with an elliptic hole

4 讨 论

宏观和微观两个尺度的计算表明，在动态载荷作用下探讨含缺陷结构的变形与失效行为时，必须考虑应力波与缺陷相互作用所带来的特有现象。它与含缺陷结构在准静态载荷作用下的变形与失效行为有着本质的区别。

在准静态加载下，不考虑结构中的应力波传播所带来的影响，认为材料内部任何时刻都处在应力平衡态。

图9 材料内部缺陷分布示意图Fig.9 Flaw distribution within material

如图9所示，Ω为关注的材料区域，∂Ω为该区域的外边界，内部存在分布缺陷，可用内边界∂Ωi、∂Ωi＋1、∂Ωi＋2等效。准静态加载下，应力将实时在区域Ω内达到平衡，并形成相应的应力分布。考虑缺陷的几何形状，那么就不可能出现所谓理想的尖端，也就不会出现应力奇点，此时应力分布为解析函数。根据最大模原理［22］，应力绝对值的最大值只可能出现在边界上，也即是说，如果结构内部存在应力极大值，则必然分布于缺陷的边界上。

需要指出的是，这里的结论只适合于外加载荷在材料内部尚未形成塑性区的情况，因为一旦塑性区形成，应力极大值就会向材料内部移动。于是在准静态假设下，塑性区只可能首先产生在边界，并随着载荷的增加由边界向材料内部扩展。

动态载荷作用下的情况则不同，区域内的应力并非实时达到平衡，考虑应力波的传播，应力波前缘应力值会出现跳跃，这就破坏了区域内的应力可导性，此时应力函数非解析，不能使用最大模原理理解应力极大值的位置。在发生应力波相互作用后，以及应力波与缺陷相互作用后，局域应力非平衡时，材料内部的应力状态非常复杂，很难用解析方式对应力极大值的位置给出合适的描述，尤其是对于存在不规则几何的缺陷，因此选择LS－DYNA程序以及分子动力学方法通过数值模拟给出其发展过程。

应力波与缺陷相互作用时，可以将内部缺陷视为内边界，应力波在内边界上进行反射，反射波与加载波的相互作用形成塑性区。此塑性区必然形成于远离缺陷边缘的位置，并由此诱导裂纹的萌生，引发裂纹的动态扩展。考虑应力波与缺陷的相互作用对塑性区形成演变以及材料失效行为的影响，虽然目前只进行了初期的数值模拟，却得到了完全不同于准静态分析的有趣现象。这可以为动态载荷作用下结构的防护提供相应的支持。

5 结 论

在宏观上利用LS－DYNA有限元程序，在微观上使用分子动力学程序，分别模拟了含有预置缺陷的薄板在动态拉伸载荷作用下的塑性区的形成演化以及裂纹的萌生与扩展过程。主要结论如下：

（1）动载作用下，从本质上讲塑性区是塑性波扫过的区域，它可以由应力波在缺陷处的反射以及应力波在空间的相互作用形成。

（2）宏观LS－DYNA程序模拟与微观分子动力学模拟共同指出：在动态载荷作用下，即使不考虑无位错区的存在，裂纹萌生的位置也可能出现在远离缺陷边缘的位置。从数学的角度看，这个现象的出现是由于动载作用下，应力波与缺陷的相互作用改变了应力函数的解析性。这个改变直接导致应力极大值空间分布不再保持准静态载荷作用下的边界分布特征。对于真实材料在动态载荷作用下缺陷附近的裂纹萌生，建议同时考虑无位错区的存在以及应力波效应带来的影响。

（3）由于裂纹扩展导致的卸载波相互作用，将在与加载平行的方向形成新的塑性区，诱导新的裂纹萌生，产生类十字架形的失效形式。这体现了应力波相互作用在裂纹萌生过程中所产生的重要影响。

［1］Fields R J，de Wit R.Plastic zone formation around an arresting crack［C］∥Knauss W G，Rosakis A J.Non－linear fracture：Recent advances.Springer，1990：231－238.

［2］Khan S M A，Khraisheh M K.A new criterion for mixed mode fracture initiation based on the crack tip plastic core region［J］.International Journal of Plasticity，2004，20（1）：55－84.

［3］Bian Li－chun，Kim K S.The minimum plastic zone radius criterion for crack initiation direction applied to surface cracks and through－cracks under mixed model loading［J］.International Journal of Fatigue，2004，26（11）：1169－1178.

［4］Gao Xin，Wang Han－gong，Kang Xing－wu，et al.Analysis solutions to crack tip plastic zone under various loading conditions［J］.European Journal of Mechanics A：Solids，2010，29（4）：738－745.

［5］Janssen M，Zuidema J，Wanhill R J H.Fracture Mechanics［M］.Spon Press，2004.

［6］张亚，强洪夫，杨月诚.复合型裂纹小范围屈服下裂尖塑性区统一解［J］.机械工程学报，2007，43（2）：50－54.Zhang Ya，Qiang Hong－fu，Yang Yue－cheng.Unified solutions to mixed mode crack tip under small scale yielding［J］.Chinese Journal of Mechanical Engineering，2007，43（2）：50－54.

［7］Huang Yi，Chen Jing－jie，Liu Gang.A new method of plastic zone size determined based on maximum crack opening displacement［J］.Engineering Fracture Mechanics，2010，77（14）：2912－2918.

［8］Shi S Q，Puls M P.A simple method of estimating the maximum normal stress and plastic zone size at a shallow notch［J］.International Journal of Pressure Vessels and Piping，1995，64（1）：67－71.

［9］张培源，张晓敏，严波，等.裂尖曲率对裂纹前缘塑性区的影响［J］.应用力学学报，2004，21（4）：93－96.Zhang Pei－yuan，Zhang Xiao－min，Yan Bo，et al.Plastic zone affected by crack tip curvature［J］.Chinese Journal of Applied Mechanics，2004，21（4）：93－96.

［10］布洛克 D.工程断裂力学［M］.王克仁，译.北京：科学出版社，1980.

［11］钱才富，姜忠军，陈平，等，裂纹尖端塑性区和无位错区的微观模拟［J］.金属学报，2004，40（2）：159－162.Qian Cai－fu，Jiang Zhong－jun，Chen Ping，et al.Micro－simulation of crack tip plastic zone and dislocation－free zone［J］.Acta Metallurgica Sinica，2004，40（2）：159－162.

［12］钱才富，李慧芳，崔文勇.Ⅰ型裂纹尖端塑性区和无位错区及其对裂纹扩展的影响［J］.材料研究学报，2007，21（6）：599－603.Qian Cai－fu，Li Hui－fang，Cui Wen－yong.ModeⅠcrack tip plastic zone，dislocation－free zone and their effects on crack propagation［J］.Chinese Journal of Materials Research，2007，21（6）：599－603.

［13］Rudd R E，Belak J F.Void nucleation and associated plasticity in dynamic fracture of polycrystalline copper：An atomistic simulation［J］.Computational Materials Science，2002，24（1／2）：148－153.

［14］Rudd R E.Void growth in BCC metals simulated with molecular dynamics using the Finnis－Sinclair potential［J］.Philosophical Magazine，2009，89（34／35／36）：3133－3161.

［15］SeppäläE T，Belak J，Rudd R E.Onset of void coalescence during dynamic fracture of ductile metals［J］.Physical Review Letters，2004，93（24）：245503.

［16］祁美兰，贺宏亮，王永刚，等.动态冲击下纯铝中微损伤演化的仿真研究［J］.振动与冲击，2007，26（8）：133－135.Qi Mei－lan，He Hong－liang，Wang Yong－gang，et al.Simulation of micro void evolution in the pure aluminum under dynamic loading［J］.Journal of Vibration and Shock，2007，26（8）：133－135.

［17］王永刚，刘宏伟.强冲击载荷下含杂质的纯铝中微孔洞长大的动力学行为［J］.高压物理学报，2010，24（4）：248－254.Wang Yong－Gang，Liu Hong－wei.Dynamic behavior of void growth in aluminum with a preexisting flaw under intense impact loading［J］.Chinese Journal of High Pressure Physics，2010，24（4）：248－254.

［18］工程材料实用手册委员会.工程材料实用手册：结构钢、不锈钢［M］.北京：中国标准出版社，1988.

［19］LS－DYNA Keyword user’s manual－2003［M］.California：Livermore，2003.

［20］郭昭亮，刘仓理，汤铁钢.预置圆孔膨胀环动态断裂行为研究［J］.实验力学，2010，25（5）：546－552.Guo Zhao－liang，Liu Cang－li，Tang Tie－gang.On the expanding ring dynamic fracture behavior with a preset circular hole［J］.Journal of Experimental Mechanics，2010，25（5）：546－552.

［21］Wagner N J，Holian B L，Voter A F.Molecular－dynamics simulations of two－dimensional materials at high strain rates［J］.Physical Review A，1992，45（12）：8457－8470.

［22］潘永亮，汪琥庭，汪芳庭，等.复变函数［M］.北京：科学出版社，2004.