基于FRANC3D的CT试样疲劳裂纹扩展计算
在断裂力学分析中一个重要的计算参量就是应力强度因子，对于一些含有裂纹的复杂结构，计算其裂纹尖端的应力强度因子只能依赖于数值方法，如通过有限元方法或边界元方法进行求解，常用的有限元分析软件有ABAQUS和ANSYS等。但这些软件并非专业的断裂分析软件，在进行裂纹网格的前处理时较为复杂。FRANC3D是由美国康奈尔大学断裂力学工作组开发的一套具有建模、应力分析、应力强度因子分析和裂纹扩展分析的计算软件，可以用来计算工程结构在任意复杂的几何形状、载荷条件和裂纹形态下的三维裂纹扩展形态和疲劳裂纹扩展寿命，并且与ANSYS和ABAQUS等软件具有接口。
本文基于有限元分析软件ABAQUS建立了CT试样(compact tension specimen,紧凑拉伸试样)的有限元模型，并通过断裂分析软件FRANC3D插入裂纹体进行应力强度因子分析，计算疲劳裂纹扩展寿命，并与解析解进行对比。
1. 应力强度因子计算
计算采用的标准CT试样的尺寸如图1所示，尺寸为W=100mm，B=25mm，初始裂纹长度a设置为30mm，计算过程取载荷P为5000N。
图1 标准CT试样
FRANC3D无法完成有限元模型的前处理过程，因此需要在ABAQUS中建立CT试样的有限元模型，并设置相应的材料属性和边界条件，但无需在ABAQUS中建立含裂纹体的网格，而是由FRANC3D自动完成裂纹体的插入及裂纹网格的划分，在ABAQUS中建立的有限元模型如图2所示。
图2 CT试样有限元模型
在ABAQUS中完成有限元模型后导出生成输出文件(.inp)，并将其导入到FRANC3D中，导入的模式可以选择两种：一种是导入完整的有限元模型(complete model)，一种为将导入的模型分割为全局模型和局部模型(globalmodel and local model)，FRANC3D将会在局部模型上插入裂纹并划分裂纹网格，这将极大地缩短划分裂纹网格所需的时间，因此非常适用于复杂结构的裂纹扩展分析。需要注意的是，与ABAQUS的子模型技术不同，尽管FRANC3D将模型分割为了全局模型和局部模型，但在进行应力分析时仍然采用的是完整的全局模型，在全局模型中，局部模型将通过绑定连接(tie)或者接触以及节点相连的形式(三种连接方式可以在FRANC3D中选择，默认为绑定连接)与局部模型之外的区域连接在一起，因此这样的处理模式只是简化了局部网格划分的规模，并没有减小模型的计算规模。
在FRANC3D中划分的局部模型如图3所示，图中的红色 区域即为局部模型，在划分裂纹网格时，网格重划分只会在红色 区域内完成，而对白色 区域没有影响，红色和白色 区域的交界面处默认以绑定接触的形式连接，因此在交界面处的网格形式将会被保留，不会被重划分，对于存在边界条件的区域(retained BC surfaces)网格也会被保留下来，由于本模型中红色 区域不存在边界条件作用的区域，因此只有交界面处的网格会被保留，其余区域的网格均会随着裂纹网格的划分而发生改变。
图3 FRANC3D分割局部模型
分割完成的局部模型如图4所示。
图4 CT试样局部模型
从图4中可以看出，交界面处的网格被保留了下来，因此在划分裂纹网格的过程中这些区域的网格不会发生改变。
通过FRANC3D中的Cracks→New Flaw Wizard插入裂纹，FRANC3D支持多种初始裂纹类型，对于CT试样可以设置为穿透裂纹，初始裂纹长度a为30mm，如图5所示。
图5 插入初始裂纹
在FRANC3D中划分裂纹网格时，为了保证能够精确计算应力强度因子，与采用奇异单元进行应力强度因子分析类似，需要裂纹前缘划分一圈环状的楔形单元和六面体单元，如图6所示。
图6 裂纹前缘网格构成
裂纹尖端由一圈楔形的奇异单元构成，这些奇异单元将会被两圈或更多的六面体单元包围，而在六面体单元的外围将会采用四面体单元进行过渡。FRANC3D在裂纹尖端默认采用8个楔形奇异单元和2圈六面体单元，这样的网格布置方式将会模板的形式(templatemesh)应用到每一个裂纹尖端的网格划分中。
由于默认的网格布置仿真已经能够获得非常高的计算精度，因此一般在分析中只需要设置裂纹前缘的半径(template radius)即可，若默认的裂尖网格布置方式不满足计算要求，也可以在FRANC3D中进行更改，如图7所示。
图7 裂纹前缘半径的定义
为了比较裂纹前缘半径对应力强度因子的影响，本文将选取不同的裂纹前缘半径进行计算。
通过FRANC3D中的Analysis→Static Crack Analysis进行应力强度因子计算，FRANC3D将自动调用ABAQUS的求解器进行应力分析以及应力强度因子计算，为了提高计算速度，可以在FRANC3D中开启ABAQUS的并行计算功能，选择Abaqus command→View/Edit command，在输入文件名的后面添加命令:cpus=n，其中n为计算参与的CPU数量，如图8所示。
图8 设置并行计算
计算得到CT试样的应力分布如图9所示。
图9 CT试样应力分布
采用裂纹扩展分析软件FRANC3D对应力强度因子进行计算，FRANC3D提供三种方法进行应力强度因子的计算：交互积分法(Interaction Integral, M-Integral)、位移外推法(Displacement Correlation)和虚拟裂纹闭合法(Virtual Crack Closure Technique，VCCT)，计算过程中默认采用精度较高的交互积分法(M积分法)进行计算，计算结果可以通过与位移外推法的结果进行比较来验证计算结果的精度。
分别取裂纹前缘半径为0.1、0.2、0.3和0.4mm，计算得到的结果与采用解析公式计算得到的结果对比如图10所示。标准GB/T 6398-2000中给出的CT试样的应力强度因子的解析计算表达式为：
其中α=a/W，上式仅在α大于等于0.2时有效。
从图10中可以看出，应力强度因子在沿试样厚度方向的分布并不是一致的，而是在试样两侧的应力强度因子较小，而在试样中心的应力强度因子达到最大值，这是由于试样在两侧和中心位置所处的应力状态不同而造成的。这种应力强度因子上的差异也将导致裂纹在扩展时并非一条直线，而是呈现出一条拇指状的曲线。同时，本文选用的三种网格尺寸计算结果差异不大，与解析结果也较为接近。
图10 应力强度因子计算结果
2. 裂纹扩展仿真
2.1解析计算方法
对于处于中速率的稳定疲劳裂纹扩展，裂纹扩展速率da/dN与应力强度因子范围ΔK之间近似满足对数线性关系，如图11所示，即Paris公式：
其中裂纹扩展参数C和m是描述材料疲劳裂纹扩展特性的基本参数，由试验确定。计算中取C为1.14×10-15，m为3.85。
图11 裂纹扩展速率曲线
上式中ΔK为应力强度因子范围，可表示为：
其中Kmax为循环载荷中最大载荷对应的应力强度因子，Kmin为循环载荷中最小载荷对应的应力强度因子，R为循环载荷的应力比，计算中取循环载荷为脉动循环载荷，即有R=0，ΔK=Kmax。
当使用Paris公式进行裂纹扩展寿命计算时，可采用数值方法进行计算，在指定的循环次数ΔN内，裂纹扩展长度Δa可近似表示为：
假定裂纹扩展之前的长度为ai-1，则当前裂纹长度ai可表示为：
通过迭代计算，即可计算得到循环次数N与裂纹扩展长度a之间的关系。上述计算过程中假设裂纹为I型裂纹，沿直线扩展，因此不存在裂纹扩展过程中的角度变化的问题。
2.2 裂纹扩展仿真计算
下面采用FRANC3D进行裂纹扩展模拟，由于在进行裂纹扩展分析时需要已知裂纹前缘的应力强度因子，因此首先需要采用Static Crack Analysis进行应力强度因子计算，在获得裂纹前缘的应力强度因子后再进行裂纹扩展仿真。
在Analysis中选择Crack Growth Analysis，FRANC3D提供多种裂纹扩展模型(Crack Extension Model)，这里选择给定应力比(Given Stress Ratio, R)的恒幅疲劳裂纹扩展分析(Fatigue, Constant Amplitude)，由于裂纹扩展仿真中不涉及角度变化，因此在裂纹扩展角度模型(Kink Angle Model)中选择平面裂纹(Planar)。
FRANC3D提供了多种裂纹扩展速率模型以及考虑应力比影响的模型用于疲劳裂纹扩展的模拟，如图12所示，这里选择不考虑应力比效应的Paris模型，除了裂纹扩展速率参数C和n，FRANC3D还要求提供裂纹门槛值ΔKth和材料断裂韧度Kc，当应力强度因子范围ΔK小于裂纹门槛值ΔKth，裂纹将不会扩展，当应力强度因子的最大值Kmax大于材料断裂韧度Kc时，材料将发生断裂，因此FRANC3D将终止裂纹扩展的计算。本文的计算过程并不具备任何工程背景，因此可以将裂纹门槛值ΔKth设定为非常小的值，将材料断裂韧度Kc为非常大的值防止计算终止。
图12 裂纹扩展速率模型设置
在计算应力强度因子范围ΔK时，通常需要计算出最大载荷对应的应力强度因子Kmax和最小载荷对应的应力强度因子Kmin，由于在线弹性断裂力学中对于给定的裂纹体，当裂纹长度确定时，载荷大小与应力强度因子的取值存在比例关系，因此FRANC3D只计算了最大载荷对应的应力强度因子Kmax，而最小载荷对应的应力强度因子Kmin直接通过应力比计算得到，如图13所示，在ABAQUS中给定载荷下计算得到的应力强度因子将作为Kmax用于裂纹扩展分析，如果ABAQUS中的载荷值对应的应力强度因子并非裂纹扩展使用的Kmax，还可以通过Load Multiplier对Kmax进行比例调整。
图13 指定循环载荷
随着疲劳裂纹的扩展，裂纹前缘的轮廓也会发生变化，FRANC3D使用多项式对裂纹前缘轮廓进行平滑拟合，避免在应力强度因子计算中产生不必要的数值噪声。FRANC3D在自动裂纹扩展模拟中默认采用固定阶次的多项式对裂纹前缘轮廓进行拟合，拟合多项式的最佳阶次由程序自动计算给出，如图14所示。
图14 裂纹前缘轮廓拟合
通常程序给出的多项式已经足够拟合裂纹前缘轮廓，只通过调整extrapolate的取值使得拟合的多项式能够外延处模型一段距离即可，如果拟合多项式完全位于网格内部，则将导致裂纹网格划分失败。
FRANC3D提供了两种裂纹扩展类型用于疲劳裂纹扩展：固定长度扩展(medianextension)和固定循环次数扩展(cyclesextension)，采用固定长度扩展时，则在每一次计算分析中裂纹扩展的长度Δa是给定的，其对应的循环次数ΔN通过裂纹扩展速率模型计算得到；采用固定循环次数ΔN扩展时，则每一计算分析时的裂纹扩展长度Δa通过裂纹扩展速率模型计算得到。当采用固定长度扩展时，用户指定的长度为裂纹前缘上的代表裂纹前缘应力强度因子平均值(median K)的节点处的扩展量，如图15所示。
图15 裂纹前缘扩展方式
若裂纹前缘应力强度应力平均值(median KI)处的节点的扩展量(由用户指定)为Δam，则在其他节点处的扩展量Δai可以表示为：
通过分析可知上式中的n即为Paris公式中的裂纹扩展速率参数m。
本文中取每个计算步中裂纹扩展增量为2mm，裂纹扩展计算步的数量为20，已知初始裂纹长度为30mm，则计算终止时裂纹长度约为70mm。
计算终止时CT试样的应力分布如图16所示。
图16 CT试样应力分布
图17 裂纹扩展前缘轮廓变化
将FRANC3D计算得到的裂纹扩展寿命与解析方法计算得到的寿命进行对比，如图18所示。
图18 疲劳裂纹扩展寿命对比
从图18中可以看出，通过解析方法计算得到的裂纹扩展寿命要高于通过FRANC3D计算得到的寿命，这主要是由于解析方法计算得到的应力强度因子要低于裂纹尖端实际的应力强度因子，这一点从图10中的结果也可以看出，因此解析方法给出的寿命要更长。
为了验证解析方法与FRANC3D的计算误差是由于计算的应力强度因子偏低造成的，本文提取了图17中虚线所示路径(裂纹中心线)上应力强度因子范围ΔK与裂纹长度a之间的关系(a-ΔK曲线)，如图19所示，根据Paris公式可以得到：
因此裂纹扩展寿命可以表示为：
其中a0和af为裂纹初始长度和最终长度，应力强度因子范围ΔK是一个与裂纹长度a有关的函数。令函数f(a)满足下述关系：
可以得到f(a)和应力强度因子范围ΔK与裂纹长度a之间的关系如图19所示。
图19 f(a)和应力强度因子范围ΔK与裂纹长度a之间的关系
从图19中可以看出，随着裂纹长度a的增加，应力强度因子范围ΔK的值也随之增加，而函数f(a)的取值则逐渐下降，由于函数f(a)与横坐标轴围成的面积的实质为该段裂纹长度内的扩展寿命，因此可以看出在裂纹扩展的后期阶段，随着循环次数的增加，裂纹长度a将迅速增长。
通过将函数f(a)进行数值积分即可得到疲劳裂纹的扩展寿命，其扩展寿命与解析计算方法和FRANC3D给出的扩展寿命的对比如图20所示。
图20 疲劳裂纹扩展寿命对比
图20中FRANC3D Integration代表通过FRANC3D提取的a-ΔK曲线积分计算得到的扩展寿命，可以看出该方法计算得到的寿命与FRANC3D给出的寿命非常接近，误差仅为2.82%，因此可以确定解析公式与FRANC3D计算得到的误差较大是由于解析公式计算的应力强度因子偏低引起的。
本文主要讨论了恒幅载荷下CT试样的疲劳裂纹扩展寿命的计算，对于变幅载荷情况下的疲劳裂纹扩展寿命计算，当不考虑载荷次序的影响时，同样可以通过Cycle-by-Cycle的方式逐个循环累计计算扩展寿命，而通过a-ΔK曲线计算将极大地简化计算流程，此时由于每个载荷谱块内的应力水平有所不同，因此可以计算单位应力水平下的a-ΔK曲线，再通过线弹性比例关系获得每个载荷谱块对应的a-ΔK曲线。
