0.   序言

随着现代焊接技术、缺陷检测技术以及焊后修复技术的发展，焊接区域内存在的初始缺陷的尺度被逐渐缩小，微裂纹萌生以及短裂纹逐渐发展至宏观裂纹，断裂力学相关理论在此过程中不再适用，因此，对焊接区域裂纹萌生及短裂纹扩展的仿真模拟成为亟待解决的问题. 对此，邓彩艳等人^[1]、刘小刚等人^[2]基于Tanaka-Mura模型对ABAQUS进行二次开发，并结合Voronoi图法建立了对钢材焊接区域的裂纹萌生仿真方法. 成立夫等人^[3]运用基于晶界不可逆损伤理论的疲劳指标参量(FIP)及Voronoi图法完成了对焊接头焊趾裂纹萌生过程的仿真. 与上述两种方法对比，循环内聚力模型(cycle cohesive zone model，CCZM)在设置裂纹萌生路径时，将具有累积损伤特性的Cohesive单元预制在金属材料晶界及晶内，因而能够得到更加直观的结果.

目前基于CCZM和Voronoi图法进行的多晶金属裂纹萌生仿真主要集中在滚动接触疲劳领域. 杨静等人^[4]将赋予疲劳损伤本构的Cohesive单元插入Voronoi微观晶粒间以模拟晶界特性，实现在滚动接触疲劳下裂纹萌生及扩展的模拟.Ghodrati等人^[5]基于Neper程序开发了包括晶粒和晶界微观结构的Voronoi模型，进行了滚动钢轮接触疲劳损坏的研究.Benedetti等人^[6]基于不可逆损伤内聚力模型模拟循环荷载下微观晶粒晶界的退化过程，建立了一种用于分析多晶材料疲劳性能的方法.Sun等人^[7]实现了相邻实体单元与Cohesive单元的损伤信息传递，结合循环内聚力模型与晶体塑性本构实现了微观尺度下的疲劳裂纹萌生仿真. 上述研究表明，结合CCZM和Voronoi图法对多晶金属材料裂纹萌生进行模拟具备可行性，而焊接区域与母材同样由微观晶粒组成，且同样包含穿晶断裂和沿晶断裂.因此，基于CCZM和Voronoi图法形成具有疲劳累积损伤特性的焊接区域、热影响区和焊缝微观晶粒模型，并与宏观对接焊缝模型合并，进行多尺度裂纹萌生仿真，得到裂纹萌生位置、短裂纹扩展路径、临界循环次数及微观晶粒组织力学响应特征，形成焊缝焊趾疲劳裂纹萌生行为仿真方法.

1.   内聚力模型

1.1   单调内聚力模型

如图1所示，裂纹尖端塑性区内微孔萌生、发展并与裂纹尖端结合，促使裂纹向前扩展，内聚力模型将裂纹视为上下彻底分离的自由面，两自由面间作用着的内聚力σ使得裂纹张开后不立即断开，将内聚力σ与张开位移δ之间演化方程定义为牵引分离定理(traction-separation law，TSL)，而在开裂过程中释放的能量定义为断裂能.

图  1  内聚力模型

Figure  1.  Cohesive zone model

下载: 全尺寸图片 幻灯片

针对不同的应用环境，研究学者提出不同类型TSL，其中双线性TSL运用的最为广泛，如图2所示，双线性TSL可以分为线弹性阶段和损伤阶段，处于线弹性阶段的内聚力区域在承载后，内聚力($ {\sigma}_{\text{n}} $为法向，$ {\sigma}_{\text{s}} $为切向)随张开位移($ {\delta }_{\text{n}} $为法向，$ {\delta }_{\text{s}} $为切向)增加而增加，并呈现线性关系，其斜率为初始刚度($ {{k}}_{\text{0,n}} $为法向，${{k}}_{\text{0,s}}$为切向)，随之达到一个应力最大值为临界应力($ {\sigma}_{\text{max,n}} $为法向，$ {\sigma}_{\text{max,s}} $为切向)，此时对应的张开位移为临界位移($ {\delta }_{\text{0,n}} $为法向，$ {\delta }_{\text{0,s}} $为切向)，随着张开位移的增加，内聚力逐渐降低直至内聚力区域完全破裂，此时对应的张开位移为破坏位移($ {\delta }_{\text{f,n}} $为法向，$ {\sigma}_{\text{f,s}} $为切向)，内聚力与牵引位移变化关系曲线围成的面积为断裂能($ {{G}}_{\text{IC}} $为法向，$ {{G}}_{\text{I}\text{I}\text{C}} $为切向).

图  2  双线性TSL

Figure  2.  Bilinear traction-separation law. (a) normal direction; (b) tangent

下载: 全尺寸图片 幻灯片

实际工程结构往往处于复杂的受力情形下，混合模式下的双线性内聚力模型可用于表征复杂应力下的断裂情形，混合模式下TSL如图3所示，图中混合模式下$ {\sigma}_{\text{max,}\text{m}}\mathrm{为} $临界应力；$ {{G}}_{\text{I}\text{-II}\text{C}} $为断裂能；$ {{k}}_{\text{0,}\text{m}} $为初始刚度；$ {\delta }_{\text{0,}\text{m}} $为临界位移；$ {\delta }_{\text{f}\text{,}\text{m}} $为破坏位移为$ {\delta }_{\text{f}\text{,}\text{m}} $.

图  3  混合模式下的TSL

Figure  3.  Traction-separation law of mixed-mode

下载: 全尺寸图片 幻灯片

切向、法向与混合模式的牵引−位移本构关系之间存在转化关系，混合模式下的位移变化量$ {\delta }_{\text{t}\text{,}\text{m}} $、初始临界位移$ {\delta }_{\text{0}\text{,}\text{m}} $、初始破坏位移$ {\delta }_{\text{f}\text{,}\text{m}} $计算方程分别为

式中：$ {\delta }_{\text{t}\text{,}\text{n}} $、$ {\delta }_{\text{t}\text{,}\text{s}} $分别为法向、切向位移变化量；β为混合加载比；η为B-K准则系数.

1.2   循环内聚力模型

循环内聚力模型本构关系如图4所示. OA阶段为单调加载条件下的线弹性阶段，此阶段刚度不发生变化，损伤并未发生；当张开位移$ {\delta }_{\text{m}} $大于临界位移$ {\delta }_{\text{0,m}} $后进入单调加载下的软化阶段(AB段)，此时材料发生不可逆损伤，刚度逐渐下降，D_t为当前时刻的累积损伤量，应力随着张开位移$ {\delta }_{\text{m}} $的增加而下降；BO段为卸载阶段，此阶段损伤不进行累积；当再次加载时(OC段)，由于不可逆损伤的产生，应力不会沿着上一次卸载时的斜率增加，而是沿着曲线OC上升至点C，此过程中单元的刚度逐渐减少，而在之后的卸载过程(CO段)中刚度不发生变化，循环反复，损伤逐渐累积，刚度逐渐下降，直至内聚力区域完全失效，裂纹向前发展.

图  4  循环内聚力模型本构关系

Figure  4.  Constitutive relation of cyclic cohesive zone model

下载: 全尺寸图片 幻灯片

Roe等人^[8]认为内聚力区域的疲劳破坏应遵循以下3个原则：①只有当变形量大于临界值时损伤才会开始累积；②损伤的增量与当前加载水平变形增量有关；③存在临界应力水平，低于这个水平循环可以无限进行而不失效. 并在此基础上提出循环荷载作用下的Cohesive单元累积损伤本构演化方程式和单调荷载作用下的软化损伤${\Delta } $D_m通过式分别为

式中：$ \Delta $D_C为疲劳累积损伤增量；A为损伤放大因子；$ \Delta\dot{\bar {{u}}} $为位移变形增量，在计算时需对其取绝对值；$ {\delta }_{\Sigma } $为内聚力特征长度，可视为初始临界位移的倍数并根据实测数据进行调整，即${\delta }_{{ \Sigma }}\text{=}\text{n}{\delta }_{\text{0,m}}$；$ \bar {{T}} $为当前时刻有效界面应力；$ {\sigma}_{\text{max}\text{,}\text{m}} $为初始临界应力；σ_f为材料的疲劳极限；H(x)为Heaviside阶跃函数，当x小于0时，H(x)为0，当x大于等于0时，H(x)为1；$ \Delta \bar {{u}} $为当前应变值，即累积的位移变形增量；$ {\delta }_{\text{0,m}} $为初始张开位移，也可根据实际情况进行调整.

1.3   疲劳累积损伤子程序

采用ABAQUS用户自定义材料本构子程序VUMAT将Roe-Siegmund损伤演化模型进行程序化，并将VUMAT定义的本构模型赋予到指定的Cohesive 单元中，以实现单元刚度退化失效. 图5为疲劳累积损伤子程序流程.

图  5  疲劳累计损伤子程序流程

Figure  5.  Fatigue cumulative damage subroutine flow

下载: 全尺寸图片 幻灯片

2.   多尺度对接焊缝疲劳断裂分析

2.1   Q345 焊接区域内聚力参数

文中以金属拉伸颈缩断裂时的真实应力^[9]作为材料的法向临界应力值$ {\sigma}_{\text{max}\text{,}\text{n}} $. 文献[10]对Q345热影响区和焊缝带缺口圆棒试件进行轴拉试验，图6为缺口半径R为1.5 mm热影响区及焊缝圆棒试件在位移−加载过程中缺口处的名义轴向应力与真实轴向应力变化图. 如图6所示，名义轴向应力在到达临界伸长量时急剧下降，此时缺口处发生颈缩断裂现象，对应的真实轴向应力即为法向临界应力，Q345热影响区、焊缝的法向临界应力$ {\sigma}_{\text{max}\text{,}\text{n}} $分别为1294.92、1296.01 MPa. 切向临界应力$ {\sigma}_{\text{max}\text{,}\text{s}} $取值为35％的法向临界应力值^[11]，即0.35$ {\sigma}_{\text{max}\text{,}\text{n}} $，分别为453.22、453.60 MPa.

图  6  热影响区与焊缝临界应力

Figure  6.  Heat affected zone and weld critical stress. (a) heat affect zone; (b) weld

下载: 全尺寸图片 幻灯片

断裂力学理论中的断裂韧性J_IC与内聚力模型中的断裂能G_C的概念相同，取值相等^[12]. 文中将Q345热影响区的断裂韧性设置为235.8 N/mm^[13]，焊缝的断裂韧性设置为781.47 N/mm^[14]，初始刚度设置为10⁸ MPa/mm^[15]. 除内聚力参数以外，还需要确定疲劳强度σ_f，热影响区与焊缝疲劳强度的取值统一取值为150 MPa^[16]. 需要说明的是，文中疲劳强度为限制裂纹萌生行为仿真进行的拟定阈值，未考虑实际过程加载中循环次数与应力比对疲劳强度的影响.

2.2   多尺度有限元模型

对接焊缝存在焊趾附近萌生裂纹沿板厚方向扩展的开裂模式，焊趾附近为焊缝、热影响区的交界处. 后续仿真模拟中，将Q345热影响区晶粒尺寸设定为10 μm，将焊缝区晶粒设定为30 μm^[17]，同时考虑晶粒的形态，将热影响区的晶粒长宽比设定接近为1∶1，而将焊缝区晶粒的长宽比设定为2.5∶1. 裂纹在金属组织中的萌生与扩展，通常包括穿晶断裂和沿晶断裂的两种断裂模式，而焊接区域通常以穿晶断裂为主，在后续模拟中考虑以上两种断裂模式，并以晶内与晶界内聚力参数赋值的差异体现以穿晶断裂为主导的断裂模式.

建立图7所示的多尺度对接焊缝有限元模型，采用单面V形焊，保留焊缝余高，焊缝宽度为6 mm，余高为1.2 mm，将对接焊缝有限元模型的一端设置为固定边界，在另一端施加应力比为0.1，最大荷载为111.11 MPa的均布荷载. 基于Voronoi图法在焊趾位置处建立半径为0.3 mm的焊缝与热影响区微观晶粒区域. 在微观晶粒晶界及晶粒内部插入0 mm厚度的Cohesive单元，赋予其Roe-Siegmund本构模型和相应的内聚力参数. 考虑微观晶粒组织内部晶粒的各向异性特征，运用Python脚本将随机的局部坐标系赋予各个晶粒，并将C₁₁ = 228 GPa，C₁₂ = 132.0 GPa，C₄₄ = 116.5 GPa各向异性线弹性本构赋予各个晶粒^[18]. Q345钢焊接接头焊趾微观晶粒组织的实际金相观察结果如图8所示^[19]，与建立的微观晶粒模型具有基本一致的分布特征. 除微观晶粒区域以外的模型部分采用各向同性线弹性本构，弹性模量取值为206 GPa，泊松比为0.3. 二维有限元模型假定取自对接焊缝中部截面，因此，微观晶粒组织模型与剩余宏观对接焊缝模型均设置为平面应变单元.

图  7  多尺度对接焊缝模型(mm)

Figure  7.  Multi-scale butt weld model

下载: 全尺寸图片 幻灯片

图  8  焊趾微观晶粒组织分布

Figure  8.  Microstructure distribution of weld toe

下载: 全尺寸图片 幻灯片

2.3   裂纹萌生仿真结果

由于晶粒组织的形态和晶粒各向异性的随机性，相同荷载循环加载作用下裂纹萌生至相同深度时的循环次数会存在于一定范围内. 文中设置10组具有不同晶粒组织和各向异性特征的仿真组，依据参考文献[20]将裂纹萌生临界深度取为0.3 mm. 其中一组仿真组(BUTT-1)裂纹萌生及短裂纹扩展过程中，图9为其x轴方向应力的变化. 由图9可知，裂纹在宏观应力最为集中的焊趾位置萌生，并随着短裂纹的扩展微观晶粒组织的应力分布处于不断更新的状态，裂纹尖端附近微观晶粒存在小范围的应力集中区域，加速了附近晶界和晶内Cohesive单元的累积损伤过程，不断更新的应力分布促使短裂纹沿板厚方向发展. 并且由于所设置的各向异性本构，微观晶粒区域的应力分布并非均匀过渡.

图  9  短裂纹扩展应力变化

Figure  9.  Stress variation of short crack propagation. (a) initial state; (b) 0.1 mm; (c) 0.2 mm; (d) 0.3 mm

下载: 全尺寸图片 幻灯片

短裂纹扩展路径如图10所示，短裂纹发展以穿晶断裂为主导，并呈现典型的“ZIGZAG”形态. 文献[21]对焊接热影响区进行微观断口分析得到如图11所示的金相观察结果，可知短裂纹在热影响区中的扩展路径特征及断裂模式与模拟结果基本一致. 因此，所建立的多尺度裂纹萌生仿真方法在进行裂纹萌生及短裂纹扩展行为模拟时具备可靠性.

图  10  短裂纹扩展路径

Figure  10.  Short crack propagation path

下载: 全尺寸图片 幻灯片

图  11  裂纹在热影响区中的扩展

Figure  11.  Crack propagation in heat affected zone

下载: 全尺寸图片 幻灯片

图12为不同仿真组中短裂纹扩展至临界深度0.3 mm时的循环加载次数. 由图12可知，由于微观晶粒组织分布及力学特性的随机性，不同裂纹萌生仿真组中短裂纹扩展至临界深度时的循环加载次数均不相同，但分布在一定范围内. 图12中标注了获得上述循环加载次数时的n值和A值，其中n值为由累积内聚力长度$ {\delta }_{\Sigma } $与临界初始位移$ {\delta }_{\text{0}\text{,}\text{m}} $之比，其值大小与材料本身相关，依靠试验数据拟合得到；A值为损伤放大因子，在子程序中设置以提高仿真效率. 真实的裂纹萌生寿命需要在保持A/n比值不变的情况下，根据试验获取的n值调整A值，将调整后的A值与循环加载次数相乘得到.

图  12  不同仿真组临界循环加载次数

Figure  12.  Critical cyclic times of different simulation groups

下载: 全尺寸图片 幻灯片

图13为不同仿真组中短裂纹在微观晶粒区域的扩展路径. 由图13可知，不同仿真组的裂纹萌生仿真结果中，裂纹均在焊趾处萌生，短裂纹扩展路的具体细节均不相同，而在短裂纹扩展的大致方向与形态上保持基本一致，此仿真结果与大量试验现象一致，一定程度上证明了其可靠性.

图  13  不同仿真组短裂纹扩展路径

Figure  13.  Short crack propagation paths of different simulation groups

下载: 全尺寸图片 幻灯片

3.   结论

(1) 基于Roe-Siegmund循环内聚力模型的多尺度疲劳裂纹萌生分析方法，在进行焊缝焊趾疲劳裂纹萌生仿真时，能够自发地选择符合实际情况的裂纹萌生位置以及短裂纹扩展路径.

(2) 基于Voronoi图法与内聚力单元法生成的焊接区域微观晶粒组织，在与宏观试件模型形成多尺度模型后，可用来完成对焊接区域材料微观断裂过程的模拟.

(3) 由多尺度疲劳裂纹萌生方法获取的临界循环次数存在于一定的分布范围内，通过试验拟合累积内聚力长度可获取对应的裂纹萌生寿命.