X80钢裂纹扩展与失效过程的多尺度模拟与实验验证

2022-01-15王汝恒钱家钰邓淞文

西南科技大学学报 2021年4期

盛鹰贾彬王汝恒钱家钰邓淞文

(1. 西南科技大学土木工程与建筑学院四川绵阳 621010；2. 西南科技大学工程材料与结构冲击振动四川省重点实验室四川绵阳 621010)

为满足长输油气管道向大口径、高压力发展的要求，针对高钢级管线钢的组织形态、力学性能、变形破坏机理等开展研究具有重要意义[1-3]。X80管线钢是当前应用最广泛的高强度钢之一，其在复杂的受力状况和工作环境中可能面临疲劳裂纹扩展、腐蚀损伤、局部屈曲失稳等多种破坏或失效风险[4-5]。

不少学者开展了X80钢在不同工况下宏观力学性能的研究。徐杰等[6]研究了不同的焊接温度对X80钢的断裂韧度的影响。刘英义等[7]研究了X80钢在土壤环境中的腐蚀行为。安腾[8]研究了高压氢对X80钢疲劳裂纹损伤行为的影响。然而，大部分学者重点研究X80钢的宏观尺度力学性能，而没有特别关注材料的微观缺陷的扩展对其宏观力学性能的影响。

为研究材料在原子尺度的微缺陷对连续介质尺度的宏观力学性能的影响，需要建立合理的多尺度耦合模型[9]。多尺度耦合计算方法的难点在于两个方面[10]：一是如何实现各物理量在不同尺度下的传递和过渡；二是如何实现边界条件与计算结果在不同尺度下的耦合。为解决这两个难题，不少学者在这个领域开展了大量的研究工作，建立了多种在原子尺度与连续介质尺度之间进行数据传递的方法，如：递进式原子-连续关联方法[11]、并发式原子-连续关联方法[12]、多尺度内聚力模型[13]等等。其中，多尺度内聚力模型的基本思想是：首先采用分子动力学方法获得材料的界面层黏结力与张开位移的关系(T-S曲线)，再将该T-S曲线的关键信息输入到有限元模型中用于进一步研究材料在宏观尺度上的力学行为。

为研究带微缺陷的X80钢的裂纹扩展与失效演化机制，本文在原子尺度建立了多组包含不同微观缺陷的X80钢拉伸微观计算模型，获得了各条包含裂纹尖端演化微观信息的T-S曲线，根据T-S曲线上关键转折点对应的原子构型图，从原子尺度观察分析了X80钢中不同微观缺陷扩展的现象和规律，并以钝裂纹模型为例，获得了其微裂纹的扩展轨迹；然后，在连续介质尺度建立了包含内聚力单元的单向拉伸CT试件有限元模型，将T-S曲线中反映出的材料核心微观特征引入内聚力模型中作多尺度分析，研究了微观缺陷对材料宏观断裂参数的影响程度。通过实验验证了多尺度数值模拟结果的可靠性。

1 基于内聚力模型的多尺度方法

1.1 分子动力学模型及势函数验证

X80钢的金相组织分析结果显示，X80钢主要由针状铁素体和少量珠光体构成：铁素体为体心立方(BCC)晶格结构，其强度、硬度不高，但具有良好的塑性与韧性；珠光体是铁素体与渗碳体的共析体，其力学性能介于铁素体与渗碳体之间，强度较高，硬度适中，塑性和韧性较好；渗碳体是X80钢中主要的强化相，是铁碳合金按亚稳定平衡系统凝固和冷却转变时析出的Fe3C，具有复杂斜方晶体结构，其硬度很高，但塑性和冲击韧性较差，易发生脆断。

建立α-Fe/Fe3C的晶体结构，其中，α-Fe为BCC结构，晶格常数a0=b0=c0=0.28664 nm，α=β=γ=90°；Fe3C为斜方晶体结构，晶格常数a0=0.506 nm，b0=0.674 nm，c0=0.451 nm，α=β=γ=90°，其每个单胞由12个Fe原子和4个C原子构成，如图1所示。

图1 Fe3C单胞结构Fig.1 Fe3C crystal structure

建立α-Fe/Fe3C分子动力学计算模型，模型尺寸为28.664 nm×14.332 nm×1.146 nm，共计39 788个原子。选择嵌入式原子势函数EAM/FS[14]，其第i个原子的势能如式(1)所示：

(1)

式中Ak，Rk，ak，rk为常数，可由第一性原理计算得到[15]。

势函数通常是经验性的，在分子动力学仿真中往往存在争议。为验证α-Fe/Fe3C的EAM/FS势函数及其参数，采用LAMMPS获得了α-Fe/Fe3C的熔点。

设置初始构型的时间步长为0.01 ps，采用三维周期性边界条件，在2.5 K的温度下弛豫500 000步。利用Nose-Hover方法将压力保持在0附近，并将系统温度由T=2.5 K逐渐升高至2 500 K。在仿真过程中，每1 000步输出一次热动力学结果。升温过程中原子平均体积随模拟温度的变化如图2所示。

图2 升温过程中原子平均体积随模拟温度的变化Fig.2 Variation of average atomic volume with simulated temperature during heating process

升温过程中，平均原子体积几乎随温度线性增加。然而，当温度增加到某一特定值时，平均原子体积突然增加。这个现象说明，系统正经历着某种相变。α-Fe/Fe3C的熔点约为 1 873 K，仿真获得的相变时的温度值为 1 990 K。两者之间的误差约为 6.2%。验证了势函数的合理性。

1.2 T-S曲线与内聚力模型

考虑到裂纹尖端的应力不可能达到无穷大，学者们发展了内聚力模型。材料的开裂、裂纹扩展和破坏是由界面单元决定的，但是在内聚力模型中，裂纹扩展区域被简化为一个零厚度区域，该区域包含两个相重叠的面，其中一个面对任意材料都视为无损伤连续介质，另一个面为连续介质单元之间的黏聚界面，用于判定材料的损伤。在加载过程中，两个面彼此分离，但内聚力单元并不会被破坏。当材料损坏发生时，黏聚界面张开并丧失刚度，这就意味着该黏聚界面连接的两端连续介质不再相连，即材料在该处断裂。基于该原因，在内聚力模型中，裂纹只能沿黏聚界面的单元传播。

在使用黏聚单元模拟裂纹扩展之前，必须首先获得黏聚界面的黏结力T与裂纹张开位移S之间的关系(简称T-S关系)。黏聚界面的分离程度是由位移间断值δ来衡量的，即等于相邻连续单元的位移u+和u-的差值，如式(2)所示：

δ=u+-u-

(2)

黏聚界面的分离量δ取决于分别作用于黏聚界面表面的正应力和剪应力值，而裂纹张开位移S即为各处黏聚界面的分离量δ之总和。判断黏聚界面是否分开的判据有两个：内聚力单元的临界分离位移δ0和临界分离应力T0。

更具体地说，在分子动力学模型中，取裂纹面上、下一定范围的原子区域为T-S区域，分别计算这两个区域内原子的平均位移，将两者之差作为裂纹的张开位移S；将T-S区域内原子应力的平均作为材料的界面黏结力T。在计算过程中，只采用了T-S区域内的原子来获得界面黏结力与裂纹张开位移，这样能有效地避免数据的分散性。将不同时刻得到的S，T的值绘制于同一坐标系中即可得到T-S曲线。

在内聚力模型中，它们之间的黏结力和张开位移的关系服从事先指定的T-S曲线。而T-S曲线的形状、极值、分布等信息充分反映了与裂纹尖端演化相关的微观信息，如位错发射、裂尖钝化、孪晶、相变等等。运用内聚力单元模拟裂纹扩展过程必须确定3个特征量：内聚力单元的初始刚度、失效阈值、失效演化参量。而这3个特征量均可从T-S曲线中获得。

综上分析可得，本文建立的基于内聚力模型分析复合微观缺陷扩展与失效的多尺度方法的基本步骤如下：

(1)建立相应的带裂纹的分子动力学模型，确定T-S区域；

(2)统计T-S区域内的原子平均应力与裂纹张开位移，并绘制T-S曲线；

(3)根据T-S曲线获得内聚力单元材料参数；

(4)建立标准紧凑拉伸试件和三点弯曲试件的有限元模型，并将通过T-S曲线获得的参数赋予内聚力单元；

(5)在内聚力模型中模拟裂纹扩展，获得相应的材料宏观力学性能参数。

需要特别说明的是，在分子动力学模型中，T-S区域的尺寸选择对结果有较大影响。在裂纹的有效影响区域内，T-S曲线应该由上升段和下降段构成。其中，在T-S曲线的上升阶段，该区域的变形处于弹性阶段，基本没有金属键的断裂；而在T-S曲线的下降阶段，裂纹尖端会出现明显钝化和材料损伤。受裂纹影响不大的区域，其T-S曲线只有应力上升段而没有下降段，该部分区域的变形始终处于弹性阶段，不应纳入T-S曲线的综合统计范围，否则计算得到的界面黏结力的值将比真实值偏小，是不合理的。研究发现，裂纹的影响区域为裂纹长度约6倍的范围内，即合理的T-S区域应取裂纹长度约6倍的范围。

2 计算结果及分析

2.1 含不同微缺陷力学模型的T-S曲线

油气管道在冶金制造过程中会产生微裂纹和微空洞等缺陷，在服役过程中，也会因为输送介质具有很强的腐蚀性而造成管道内表面的腐蚀开裂，从而形成微观缺陷。这些缺陷在管道内的高压作用下可能逐渐扩展形成宏观裂纹，存在使管道破裂的风险。因此，建立如图3所示的分子动力学微观尺度计算模型：模型Ⅰ不含微缺陷；模型Ⅱ为中心钝裂纹(裂纹长2.5 nm，宽0.52 nm)；模型Ⅲ是为中心空洞(空洞直径0.6 nm)。

图3 微观尺度计算模型Fig.3 Microscale computation models

设置[100]和[010]方向为自由边界条件(s)，在[001]方向采用周期性边界条件(p)。分子模拟的初始温度为300 K，加载方式为沿[100]方向施加速率为0.025 nm/ps的拉伸应变。模拟中取时间步长t=1 ps，求解牛顿运动方程的积分方法选用Velocity-Verlet算法。模拟过程分两步：在微正则系综(NVE)下，首先对初始构型弛豫使其达到位能最低状态，然后对弛豫后的纳米单晶体沿[100]方向均匀拉伸20 000步。加载过程中利用Nose-Hoover方法调节温度使其保持恒定。模拟过程中每200步记录原子的位置、应力、总能量，并输出系统的几何构型。取距裂纹上、下表面1.0 nm范围内的区域为T-S数据区。

图4、图5、图6分别为通过分子动力学模拟获得的模型Ⅰ、模型Ⅱ、模型Ⅲ的T-S曲线。在裂纹张开的任一时刻，T-S曲线中都能得到一个离散点并对应一种原子构型。因此，裂纹尖端的复杂变形机制能综合反映到T-S曲线中。为了观察微观变化(如位错、裂尖钝化、二次裂纹萌生与扩展等等)对T-S曲线的影响，在T-S曲线中标出了关键转折点所对应的原子构型图。

图4 模型Ⅰ在拉伸载荷下的T-S关系与变形机制Fig.4 T-S relationship and deformation mechanism of model Ⅰ under tensile load

图5 模型Ⅱ在拉伸载荷下的T-S关系与变形机制Fig.5 T-S relationship and deformation mechanism of model Ⅱ under tensile load

图6 模型Ⅲ在拉伸载荷下的T-S关系与变形机制Fig.6 T-S relationship and deformation mechanism of model Ⅲ under tensile load

为直观比较模型Ⅰ、模型Ⅱ、模型Ⅲ的T-S曲线的差异，将3组T-S曲线放置在一张图中，如图7所示。

图7 模型Ⅰ、模型Ⅱ、模型Ⅲ的T-S曲线对比Fig.7 T-S curve comparison of model Ⅰ, Ⅱ and Ⅲ

由图4-图7可得如下结果：

(1)图4显示了无缺陷模型Ⅰ在拉伸荷载作用下原子键断裂、微裂纹萌生和扩展、模型断裂的过程，其T-S曲线走势与宏观尺度试件拉伸的应力-应变曲线走势相似。T-S曲线大致分为4个阶段：第一阶段的T-S曲线近似呈线性关系，最高点①处的原子应力为12.75 GPa；第二阶段随着裂纹张开位移的不断增大，原子应力增长缓慢，由①和②的原子构型对比可见，处于模型中线附近的原子排列方式由穿插排列变为并行排列，裂纹附近原子呈无序排列状态；第三阶段原子应力随裂纹张开位移的增大而继续增大，局部出现裂纹，当应力峰值达到19.89 GPa时，裂纹完全贯通；第四阶段原子应力急剧下降，模型完全断裂。

(2)图5显示了模型Ⅱ的裂尖变形过程，其裂尖变形及失效机制为钝化和位错发射。模型Ⅱ的T-S曲线出现了两次峰值。当T-S曲线第一次达到峰值点时，裂纹尖端应力集中严重，裂纹两端与裂纹长度方向呈45° 的方位出现位错，原子出现无序排列状态，随后应力小幅下降；随着加载的继续，裂纹钝化程度进一步加大，位错区域不断扩大，裂纹前方不远处出现二次裂纹，在主裂纹与二次裂纹连通前，应力有一次小幅上升，当主裂纹与二次裂纹间的大量连接原子都呈无序排列状态时，应力达到第二个峰值9.52 GPa；随着进一步加载，应力迅速下降，直至裂纹完全连通。

(3)图6显示了模型 Ⅲ 的裂尖变形过程，其裂尖变形及失效机制仍为钝化和位错发射。模型Ⅲ的T-S曲线出现了3次峰值。当T-S曲线第一次达到峰值点10.08 GPa时，空洞形状由圆形拉伸为椭圆形，空洞附近的原子间距增大；随着加载的继续，空洞两端沿垂直于加载的方向出现位错，位错前方原子呈无序排列状态，当位错导致的呈无序排列状态的原子大量增多时，空洞裂纹扩展为条状并逐渐钝化，出现应力松弛，应力先小幅下降再小幅上升，并出现第二次峰值；继续加载，裂纹继续扩展，扩展路线呈非线性，应力小幅下降后上升至10.56 GPa，此时模型完全断裂。

由图4-图7可进一步发现：

(1)由图5与图4对比可知，当模型完全断裂时，模型Ⅱ裂纹张开位移约为模型Ⅰ裂纹张开位移的50%，模型Ⅱ最大原子应力也约为模型Ⅰ最大原子应力的50%。表明在微观尺度，钝裂纹对模型的强度具有不可忽略的影响。

(2)由图5与图6对比可知，当模型完全断裂时，模型 Ⅲ 裂纹张开位移比模型Ⅱ裂纹张开位移约大20%，模型 Ⅲ 最大原子应力也比模型Ⅱ最大原子应力略大。表明在微观尺度，空洞对模型的强度仍具有不可忽略的影响，但空洞对模型强度的影响比钝裂纹对模型强度的影响更小。

(3)由图7可知，带微缺陷的X80钢的微观模型的破坏情况存在如下共同点：T-S曲线中的最高峰值点表明材料中开始出现损伤(裂纹起裂)，裂尖开始出现明显的位错发射；随着加载的继续，应力集中导致原子间金属键被破坏，当应力集中区域附近的大量原子呈无序排列状态时，裂尖明显钝化，应力出现松弛，应力小幅下降。

(4)在微观尺度，微缺陷会降低模型强度，不同形状的微缺陷对模型强度的影响不同，在外荷载作用下不同形状缺陷扩展的方向也不同。因此，研究不同微缺陷对微观模型强度的影响是有意义的。

(5)T-S曲线综合反映了裂纹尖端复杂的变形机制(如位错发射、裂尖钝化等等)，因此，将分子动力学得到的反映微观特性的T-S曲线作为搭建多尺度计算的桥梁是可行的。

2.2 微裂纹扩展轨迹与实验验证

由上述分析可知，钝裂纹对微观模型强度的影响比空洞更大，对X80钢而言，钝裂纹更容易导致材料沿滑移面破坏。故以钝裂纹模型Ⅱ为算例，分析其微裂纹的扩展轨迹。

裂纹尖端原子的扩散常数可以通过Einstein布朗动力学模型得到：

(3)

式中：N为系统原子总数；S为模型维数，当前系统采用三维模型，故S=3；t为时间。

通过原子位置及最大应力，即可确定裂尖原子位置及位移，进而计算裂尖原子的扩散系数D0。

微观模型下，裂纹的扩展路径由位移函数(4)控制：

(4)

式中：(X(t),Y(t))，(X(t+Δt),Y(t+Δt))分别表示当前时刻、下一时刻的裂尖坐标；Vm为裂纹扩展临界速度，本文取Vm=605.74 m/s，Δt=0.01 μs；G(t)是由C++程序生成的一个标准正态随机变量。

于是，由式(4)可迭代得到裂尖原子下一时刻的位置。图8为微观尺度下裂纹尖端的轨迹。

图8 微观尺度下裂纹尖端的轨迹Fig.8 Trace of crack tip at microscale

由图8可知，裂尖原子在微观尺度下是剧烈震荡的，而在宏观尺度作有限元计算得到的结果近似为一条直线。由此可以看出，从微观尺度研究裂纹的扩展是有必要的。

下面采用X80管件母材三点弯曲(SEB)试样的断裂韧性实验验证裂纹扩展轨迹的非线性[16]。图9所示为试样CTOD试验结果的δ-Δa试验关系及断口照片。

观察图9可知，X80的冲击断口呈非镜像对称，断口沿裂纹延长线方向的扩展呈锯齿状，并非为严格的直线。这与分子动力学计算得到“裂纹尖端的轨迹在微观尺度下是剧烈震荡的”结论一致，这也间接定性验证了本文分子动力学计算结果的合理性。

图9 CTOD试样的δ-Δa试验关系及断口照片Fig.9 δ-Δa experiment curve of Sample CTOD and the photographs of fracture surface

2.3 T-S关系对宏观断裂的多尺度模拟

为分析原子尺度下微结构的演化对较高尺度下材料失效行为的影响，本节建立了包含内聚力单元的CT试件模型，并模拟了CT试件在位移加载条件下的裂纹扩展过程。通过对CT试件裂纹扩展的模拟，能检验分子动力学获得的黏聚单元法则(T-S曲线)能否满足线弹性断裂力学的理论预期结果。CT试件的几何尺寸和网格模型如图10所示。CT试件的宽和高分别为1.250W和1.200W，其中W为销孔中心到不含裂纹边界的距离，取W=200 μm。初始裂纹长度为a=0.500W，销孔半径为 0.250W。在预定义裂纹扩展路径上设置宽为0.1 nm的内聚力单元。在参考点RP-1和RP-2上施加位移载荷。