复合型疲劳裂纹扩展角度计算方法

1. 简介

复合型疲劳裂纹扩展断裂是机械构件主要的断裂形式之一，由于其在实际结构中存在的广泛性和理论分析的复杂性而一直受到学者的广泛关注。而以I型（张开型）和II型（滑开型）并存的复合型断裂模式是导致疲劳断裂的主要原因，因此研究I-II复合型疲劳裂纹的扩展机理对构件剩余寿命的精确评估具有重大意义^[1]。复合型疲劳裂纹的扩展需要确定两个参数：裂纹尖端的扩展角度以及扩展长度。裂纹扩展角度一般需要结合相应的断裂准则来确定，而裂纹扩展长度需要首先计算裂尖等效应力强度因子K_eq，随后结合相应的裂纹扩展速率模型计算裂纹扩展长度。无论是裂纹扩展角度还是扩展长度，都需要预先计算复合型疲劳裂纹的I型应力强度因子和II型应力强度因子。

本文以一个含中心斜裂纹的均匀拉伸平板和一个无限大平面中受滚动接触载荷的斜裂纹为例，结合相应的裂纹扩展准则，介绍如何利用Abaqus计算复合型疲劳裂纹扩展角度，并重点讨论当I型应力强度因子和II型应力强度因子在单个循环周次内非等比例变化时，裂纹扩展角度的确定方法。

2. 复合型断裂准则

目前国内外已经提出了大量用于计算复合型裂纹扩展角度的断裂准则，这些准则一般都是从以下三个方面进行推导得到的：（1）以应力作为参数；（2）以位移作为参数；（3）以能量作为参数。目前常见的几种复合型断裂准则有：最大切向应力断裂准则，最大能量释放率断裂准则，最大应变能密度因子准则^[2]。这些准则由于考虑的角度和观点有所不同，因此所得结果有一定差异，但都能得到较为相近的结果。

这里本文使用最大切向应力断裂准则（Maximum Tangential Stress criteria, MTS）来评估复合型疲劳裂纹的断裂角度，该准则由Erdogan和G. C. Sih^[3]根据中心斜裂纹承受均匀拉伸的树脂玻璃板的试验结果提出。

在I-II复合型断裂问题中，裂纹尖端附近的应力场可通过下式以极坐标形式给出：

为使得裂纹沿垂直于最大环向切应力的方向扩展（尽管根据弹性理论，裂纹尖端的应力应该是无穷大的，若假设裂纹尖端的真实应力是有限的，并且弹性理论仍然成立，则还是可以使用这样的准则），应使得最大环向切应力（注意，这里的切向应力（Tangential stress）不是剪切应力（Shearstress））取得极大值，根据多元函数求极值有：

求解可得到：

因此有：

上式中第一个解θ=±π显然对应裂纹沿自由面扩展的情况，这并非我们想要的解。而第二个解中θ为随K_I和K_II变化的函数。从上式中还可以看出，该表达式与裂纹尖端剪切应力的表达式是相同的，这意味着当剪切应力为0时环向切向应力达到最大值。

为了便于计算，可以将上式中的扩展角度θ解出，有：

其中当K_I/K_II>0时上式取负号，反之取正号。

裂纹扩展角度θ也可以采用下面的等价形式进行计算：

裂纹扩展角度θ应从逆时针方向量起，当K_II=0时，θ=0；当K_II>0时，θ<0；当K_II<0时，θ>0。

3. 计算算例

3.1 含中心斜裂纹的均匀拉伸平板

假定有一块50mm×100mm×6mm的平板，平板中心含有一条长度2c为10mm的斜裂纹，裂纹与水平面的倾角β为30°。在平板远端承受有载荷比R=0、应力幅值Δσ=200MPa的均匀循环拉伸载荷，如图3.1所示，试计算斜裂纹的扩展角度。

图3.1 含中心斜裂纹平板

为了计算复合型疲劳裂纹的扩展角度，首先需要获得裂纹尖端I型应力强度因子和II型应力强度因子在循环载荷作用下的变化规律。在Abaqus中可利用奇异单元建立裂纹体，并通过围线积分法计算复合型裂纹的应力强度因子，相应的计算流程可参考如何在ABAQUS中构建奇异单元。通过计算可得到循环载荷作用下含中心斜裂纹平板的Von Mises应力分布如图3.2所示。

图3.2 含中心斜裂纹平板的Von Mises应力分布

从图3.2中可以看出，复合型裂纹尖端的应力分布与纯I型裂纹尖端的应力分布非常类似，都呈现蝴蝶形状。由于II型断裂模型下的应力是反对称分布的，可以看到上下裂纹面的应力分布是非对称的。

图3.3给出了在一个循环周次内裂纹尖端的应力强度因子变化规律。

图3.3 中心斜裂纹平板应力强度因子变化规律

从图3.3中可以看出，由于计算模型中认为材料满足线弹性、小变形假设，并且不存在其他的非线性因素，因此I型应力强度因子和II型应力强度因子在循环周次内完全是比例变化的。利用节2中的公式计算裂纹尖端不同时刻的扩展角度，结果如图3.4所示，为了对比计算结果的准确性，图中还给出了利用Abaqus内置的基于最大切向应力断裂准则的算法输出的裂纹扩展角度(Crackpropagation direction, Cpd)。

图3.4 中心斜裂纹平板裂纹扩展角度

可以看到，由于I型应力强度因子和II型应力强度因子是等比例变化的，因此在整个循环周次内，疲劳裂纹的扩展角度将为恒定值，并且这个角度约为46°，如图3.4所示。

图3.4 中心斜裂纹扩展角度示意图

尽管利用I型应力强度因子和II型应力强度因子的变化规律我们计算得到了复合型疲劳裂纹的扩展角度，但上述方法并未解决裂纹在何时扩展的问题。并且，在某些特殊情况下，I型应力强度因子和II型应力强度因子在循环周次内可能是非等比例变化的，这意味着对于不同时刻t，我们将计算得到不同的裂纹扩展角度θ，显然裂纹只可能沿其中的某一个角度扩展。

Sih和Erdogan^[3]假设当裂纹尖端的等效应力强度因子（Resultant stressintensity factor）超过材料的门槛值时疲劳裂纹将会扩展，而拉伸断裂模式下的等效应力强度因子K_σ定义为：

显然在不同时刻t，不同的K_I、K_II和θ都将得到不同的等效应力强度因子K_σ，为了获得裂纹实际扩展的角度，我们可以假设裂纹将会在等效应力强度强度因子K_σ达到最大的时刻扩展^[4]。利用图3.3中的应力强度因子和图3.4中的扩展角度，通过上式可计算得到不同时刻的等效应力强度因子K_σ变化如图3.5所示。

图3.5 中心斜裂纹平板等效应力强度因子变化规律

从图3.5中可以看出，尽管当I型和II型应力强度因子等比例变化时裂纹扩展角度为恒定值，但不同时刻的等效应力强度因子K_σ却是变化的，并且当t=2.5s时达到最大值。该时刻恰好对应于远端拉伸载荷达到最大值的时刻的，因此认为裂纹在该时刻扩展是比较合理的。

3.2 无限大平面中承受滚动接触载荷的斜裂纹

假定无限大平面中存在一条长度为1mm斜裂纹，裂纹与水平面呈30°的倾角，平面上作用有滚动接触载荷，将以匀速v滚过裂纹，试计算滚动接触载荷穿过裂纹前后的应力强度因子变化，以及裂纹的扩展角度。

图3.6 无限大平面中承受滚动接触载荷的斜裂纹

在本例中，将作用平面的滚动接触简化为Hertz接触。因此，作用在平面上的接触压力p(x)和切应力τ(x)可表示为：

其中，p₀为接触中心的最大接触压力，本例中取为100MPa，a为接触斑的半径，取为0.5mm，µ为摩擦系数，取为0.15。假定接触斑沿平面以0.1mm/s的速度移动，为了在Abaqus中模拟移动接触斑，需要分别使用用户子程序DLOAD和UTRACLOAD，其中DLOAD用于施加移动的接触压力载荷，UTRACLOAD用于施加切应力载荷。

在Abaqus中建立的有限元模型如图3.7所示，整个模型采用平面应变单元划分。尽管在接触压力的作用下，斜裂纹的裂纹面可能会处于相互接触的状态，出于简化考虑，本文不考虑裂纹面之间的接触效应。并且，似乎当考虑裂纹面接触之后，即使在裂纹面使用非常细密的网格，Abaqus的围线积分法也无法给出收敛的应力强度因子。

图3.7 含斜裂纹的无限大平板有限元模型

图3.8给出了在滚动接触载荷作用的不同阶段，斜裂纹附近的Von Mises应力分布，图中滚动接触载荷的运动起始点为距离裂纹5mm的位置。

图3.8 无限大平面斜裂纹区域Von Mises应力分布

从图中可以看出，当滚动接触载荷远离裂纹体时，Von Mises应力的峰值出现在接触平面的次表面区域，并且分布规律与Hertz接触下的应力分布规律一致。随着接触载荷逐渐靠近裂纹体，裂纹尖端的应力将会增大，使得峰值应力出现在裂纹尖端，并且靠近接触面的裂纹表面出现了较大的变形，使得裂纹面可能处于闭合状态。

图3.9给出了滚动接触载荷作用下裂纹尖端的应力强度因子变化规律。

图3.9 滚动接触载荷作用下裂纹尖端应力强度因子变化规律

可以看到，由于未考虑裂纹面之间可能存在的接触，当接触载荷靠近裂纹区域时出现了I型应力强度因子为负的情况，在真实环境中这显然是不可能出现的。利用节2中的裂纹扩展角度计算公式，可以计算得到不同时刻裂纹尖端的扩展角度如图3.10所示。

图3.10 滚动接触载荷作用下裂纹扩展角度变化规律

利用图3.10中的裂纹扩展角度，可以计算得到裂纹尖端的等效应力强度因子K_σ随时间的变化规律如图3.11所示。

图3.11 滚动接触载荷作用下裂纹尖端等效应力强度因子变化规律

从图3.11中可以看出，当t到达约为55s时等效应力强度因子达到最大值，然而在图3.10中，当t大于45s之后，通过公式计算得到的裂纹扩展角度与Abaqus给出的扩展角度有很大的差异。这主要是由于本文建立的模型未考虑裂纹面之间的接触，当滚动接触载荷穿过裂纹时，将计算得到负的I型应力强度因子，由于真实裂纹在接触之后不会存在负的I型应力强度因子，因此Abaqus会强制将裂纹扩展角度设定为0，导致公式计算得到的角度与Abaqus给出的结果有很大的差距。

尽管通过在裂纹面之间建立接触可以避免由于裂纹面之间的穿透产生的负的I型应力强度因子，但目前似乎当考虑裂纹面接触时，利用Abaqus的围线积分法无法获得收敛的应力强度因子（在每个围道上获得的应力强度因子分散性很大），或许只能使用诸如虚拟裂纹闭合法以及位移外推法来评估裂尖应力强度因子。

总结

本文以一个含中心斜裂纹的均匀拉伸平板和一个无限大平面中受滚动接触载荷的斜裂纹为例，结合最大切向应力断裂准则，讨论了复合型疲劳裂纹扩展角度的计算方法，并通过引入等效应力强度因子，来解决复合型裂纹中应力强度因子非等比例变化时裂纹扩展角度的计算问题。在承受滚动接触载荷的斜裂纹中，若不考虑裂纹面之间的接触，将得到不真实的负值I型应力强度因子，从而计算得到错误的裂纹扩展角度。目前，Abaqus的围线积分法似乎无法在考虑裂纹面接触时给出收敛的应力强度因子，必须使用诸如虚拟裂纹闭合法以及位移外推法来评估裂尖应力强度因子。

参考文献

[1]  贺红凯. 基于扩展有限元法的复合(Ⅰ-Ⅱ)加载条件下CTS疲劳裂纹扩展行为研究[D]. 南昌大学, 2014.

[2]  程靳，赵树山. 断裂力学[M]. 北京：科学出版社，2006

[3]  F. Erdogan, G. C. Sih. On the Crack Entensionin Plates Under Plane Loading and Tranverse Shear[J]. J. of Basic. Eng. ASME,1963, 85:519.

[4]  David I. Fletcher. Numerical simulation ofnear surface rail cracks subject to thermal contact stress[J]. Wear 314(2014)96-103.

附录A. 滚动接触载荷计算子程序

      !程序用于模拟二维滚动接触载荷      SUBROUTINE DLOAD(F,KSTEP,KINC,TIME,NOEL,NPT,LAYER,KSPT,     1 COORDS,JLTYP,SNAME)C      INCLUDE 'ABA_PARAM.INC'C      DIMENSION TIME(2), COORDS (3)      CHARACTER*80 SNAME      REAL::P0=100         !最大接触应力      REAL::X_INI=-5        !接触斑初始位置      REAL::RADIUS=0.5      !接触斑半径      REAL::VELOCITY=0.1    !移动速度      REAL::X_CENTER        !接触斑中心坐标      REAL::X_LOCAL         !接触斑局部坐标系中的坐标            !计算当前积分点的局部坐标值      X_CENTER=VELOCITY*TIME(1)+X_INI      X_LOCAL=COORDS(1)-X_CENTER      IF(ABS(X_LOCAL).LE.RADIUS)THEN      !积分点位于接触斑内部          F=P0*SQRT(1-X_LOCAL**2/RADIUS**2)      ELSE          F=0      ENDIF            RETURN      END            SUBROUTINE UTRACLOAD(ALPHA,T_USER,KSTEP,KINC,TIME,NOEL,NPT,     1 COORDS,DIRCOS,JLTYP,SNAME)C      INCLUDE 'ABA_PARAM.INC'C      DIMENSION T_USER(3), TIME(2), COORDS(3), DIRCOS(3,3)      CHARACTER*80 SNAME            REAL::P0=100         !最大接触应力      REAL::X_INI=-5        !接触斑初始位置      REAL::RADIUS=0.5      !接触斑半径      REAL::VELOCITY=0.1    !移动速度      REAL::FRIC_COEF=0.15  !摩擦系数      REAL::X_CENTER        !接触斑中心坐标      REAL::X_LOCAL         !接触斑局部坐标系中的坐标            !计算当前积分点的局部坐标值      X_CENTER=VELOCITY*TIME(1)+X_INI      X_LOCAL=COORDS(1)-X_CENTER      !指定切向力作用方向      T_USER(1)=-1.0      T_USER(2)=0.0      T_USER(3)=0.0      IF(ABS(X_LOCAL).LE.RADIUS)THEN      !积分点位于接触斑内部          F=P0*SQRT(1-X_LOCAL**2/RADIUS**2)          ALPHA=FRIC_COEF*F               !切向力      ELSE          ALPHA=0      ENDIF      RETURN      END