Cu/WCp 叠层功能梯度材料疲劳裂纹扩展速率的数值模拟

郭荣鑫，索玉霞，牛治亮，夏海廷，颜 峰，张玉波，李一博

(1.昆明理工大学 建筑工程学院 云南省土木工程防灾重点实验室(筹)，云南 昆明 650500；2.昆明理工大学 云南省先进材料力学行为与微结构设计高校重点试验室，云南 昆明 650500；3.云南农业大学 建筑工程学院，云南 昆明 650201；4.黄河勘测规划设计有限公司，河南 郑州 450000)

1 引 言

功能梯度材料(Functionally Gradient Materials，FGMs)是为了满足在极限温度环境(超高温、大温度落差)下正常工作而开发的一种新型复合材料[1]，Cu/WCpFGMs具备高强度、高硬度、高导电性和优异的热稳定性等优异力学性能，以及较好的耐电弧侵蚀性、抗熔焊性和耐高温抗氧化性，被广泛应用于高压断路器、继电器和开关柜等设备中的高压电触头[2-3]，在电阻焊电极及电火花加工电极等领域也有着广阔的应用前景[2-3]。

随着计算机技术的发展，以及无网格伽辽金法(EFG)[4]和扩展有限元法(XFEM)[5]等方法的提出，国内外学者对疲劳裂纹扩展行为的研究取得了重要进展，尤其在一些断裂参数(如应力强度因子、J积分等)的计算[6-7]和裂纹扩展路径的观测上，但对于裂纹扩展速率方面的研究仍然存在很多问题。杨利等[8]运用分子动力学方法，研究了微观尺度下铌元素对单晶γ-TiAl裂纹扩展过程的影响，但该方法针对的是合金材料，仅得到了γ-Ti Al裂纹扩展的轨迹图和能量演变图，并未得到裂纹扩展速率的结果。Harrison[9]提出了一种运用ABAQUS软件计算普通试件疲劳裂纹扩展速率的方法，该方法通过ABAQUS脚本界面进行模拟，提取出裂纹尖端相关的应力、应变值，用来计算各积分点的疲劳寿命，对裂纹扩展进行定义后，结合Paris公式，得到了α/β钛合金的疲劳裂纹扩展速率曲线。然而，该方法需要通过确定裂纹尖端的应力场来确定疲劳寿命，且需要自行定义裂纹尖端网格，通过相关公式计算每一个分析步中的裂纹增量，虽然结果精确但整个裂纹扩展的分析比较复杂。Nittur等[10]提出了一种采用有限元方法和塑性耗能标准相结合的数值计算方法，得到了铬镍铁合金和钛合金试件在循环载荷下的疲劳裂纹扩展速率，但该研究仅仅是针对合金等均质材料，对 FGMs 的研究还有欠缺。Bhattacharya[11]和王志勇[12]等运用扩展有限元方法(XFEM)，研究了FGMs的疲劳裂纹扩展行为，利用相互作用积分法计算应力强度因子，得到了FGMs的疲劳寿命及裂纹扩展路径。然而，这种方法在辅助应力场和位移场的确定上较复杂，并且未得出材料的疲劳 裂 纹 扩 展 速 率。JEONG-HO KIM 等[13-14]对FRANC2D 程序进行了改进，对FGMs的裂纹扩展行为进行数值模拟，但只得出了FGMs的应力强度因子值和裂纹路径的变化情况，该种方法研究的是材料参数梯度呈指数变化的FGMs，对于叠层FGMs而言，在模型的建立、材料参数的输入上较复杂。

因此，如何快速建立叠层FGMs的有限元模型，快速计算出叠层FGMs在裂纹扩展中的应力强度因子，从而得出疲劳裂纹扩展速率曲线，是亟待解决的问题之一。基于FRANC2D(Fracture Analysis Code in 2 Dimensions)有限元断裂分析软件，建立了一种FRANC2D 与ABAQUS相结合的FGMs裂纹扩展模拟方法，对Cu/WCp双层及多层FGMs标准拉伸试样的疲劳裂纹扩展行为进行数值模拟，该方法可以快速地建立叠层FGMs的有限元模型，迅速地计算出裂纹扩展中的应力强度因子，从而绘制出疲劳裂纹扩展速率曲线，在比较复杂的叠层FGMs裂纹扩展模拟研究中有很大优势。

2 FGMs裂纹扩展速率的数值模拟方法

FRANC2D 断裂有限元分析软件是由美国Cornell大学断裂力学小组开发的，是用于裂纹分析的专业软件，其最大优点在于在模拟自动扩展裂纹的同时，可以自动重置网格，并计算裂纹尖端的应力强度因子(SIF)，节约了时间又提高了计算精度[15-16]，在模拟裂纹的起裂和扩展方面有很大优势。

2.1 应力强度因子

在FRANC2D 中，对裂纹扩展过程中应力强度因子的计算主要基于三种方法：位移相关法(DCT)、虚裂纹闭合法(MD CRK-CLOS)和J 积分法(JINTEGRAL)。Ⅰ型裂纹的应力强度因子可以使用位移相关法来进行计算[17]，这种方法基于节点位移来计算用应力强度因子，且在有限元方法计算中，基于位移公式的计算方法比基于应力的计算方法更为精确。

采用DCT 方法对层间应力强度因子进行计算，当设置裂纹和裂纹扩展增量Δa后，FRANC2D 将对裂纹尖端和裂纹尖端附近的单元重新进行单元结构划分，由裂尖附近裂纹面的节点位移，可以实现每一步裂纹扩展后应力强度因子的计算。通过DCT 方法得到的应力强度因子为[18]：

式中，G=E/(1+ν)，为剪切模量；ν为泊松比；对于平面应变状态，к=3-4ν；对于平面应力状态，к=(3-4ν)/(1+ν)；a、b、c、d分别为裂纹尖端有限元单元上的节点，如图1所示；V a、V b、V c、V d分别为裂纹尖端相对于节点a、b、c、d处的y方向位移；U a、U b、U c、U d分别为裂纹尖端相对于节点a、b、c、d处的x方向位移；Δa为裂纹尖端的单元长度。由DCT 方法计算出的应力强度因子将被代入所选定的断裂判据中，断裂判据用来判断局部裂纹的稳定性及预测裂纹扩展角度，新裂纹尖端的位置则完全取决于所设定的裂纹扩展增量Δa。

图1 裂纹尖端附近单元划分Fig.1 Element division near crack tip

2.2 裂纹扩展速率

如果在应力循环ΔN次后，裂纹扩展量为Δa，则应力每循环一周，裂纹扩展量为Δa/ΔN(mm/周)，这称为裂纹扩展速率。对于疲劳裂纹扩展速率的研究，主要在于寻求裂纹扩展速率与有关各种力学参量之间的数学表达式。断裂力学认为裂纹的扩展速率同裂纹尖端应力强度因子幅值ΔK有关，P.C.Paris提出的描述这种关系的裂纹扩展公式为：

式中，a为裂纹长度，N为循环次数，C和m是与试验条件(环境、频率、温度和应力比R 等)有关的材料常数。ΔK称为应力强度因子幅，是每个周期内的最大与最小应力强度因子之差：

对于现有的有限元软件来说，一般很难直接输出疲劳裂纹扩展速率曲线，同样，FRANC2D 也不能直接给出疲劳裂纹扩展速率曲线。因此，本研究根据Cornell大学断裂力学小组所提供的表格，通过FRANC2D 计算出每一步裂纹扩展增量所对应的的应力强度因子值，从而确定出应力强度因子幅(ΔK)，结合Paris公式以及由试验确定的C和m值，可以绘制出疲劳裂纹扩展速率曲线(da/d N-ΔK)；再根据裂纹尖端距界面的距离，绘制出随裂尖距界面距离变化的裂纹扩展速率曲线(da/d N-l)。

本研究结合ABAQUS和FRANC2D 两个有限元软件，对叠层FGMs的裂纹扩展行为进行数值模拟，最终得到疲劳裂纹扩展速率曲线，过程流程图如图2所示。

3 FGMs疲劳裂纹扩展行为数值模拟过程

为研究Cu/WCp叠层FGMs的疲劳裂纹扩展行为，验证本研究中所采用数值模拟方法的准确性及可靠性，文中基于FRANC2D 断裂分析软件，分别对疲劳试验所采用的Cu/WCp叠层FGMs拉伸试件进行疲劳裂纹扩展数值模拟。

3.1 有限元模型的建立

疲劳试验研究所采用的试件为双层及多层Cu/WCpFGMs哑铃状单边缺口拉伸试件，经粉末冶金热压烧结制备而成，试件平行段尺为30×12×5mm，厚度为2mm，如图3所示。初始裂纹分别在不同WCp含量方向。试件分层布置为：双层哑铃状单边缺口拉伸试件中，梯度层的厚度从左到右分别为3和9mm，WCp体积分数变化为3%～15%；多层哑铃状单边缺口拉伸试件中，梯度层的厚度从左到右依次为2、1和9mm，WCp体积分数变化为3%～9%～15%。

因此，本研究使用ABAQUS 有限元软件，完成FGMs几何模型的建立、材料参数的输入(弹性模量E和泊松比ν)以及有限元网格的划分，再经MATLAB程序将所建模型转化为适用于FRANC2D 断裂模拟的输入程序。这种建模方法简单方便，能够得到已经赋予了材料属性并进行了网格划分的有限元模型，在FRAN2D 中只需对模型的边界条件(荷载和约束)进行设置，即可进行断裂模拟。

在ABAQUS有限元软件中建立试件的结构模型和有限元模型，如图4所示。采用六节点三角形单元进行自由网格划分，在模型上下两边和弧形部分所布的种子长度为0.2mm，考虑到应力应变在裂纹尖端具有奇异性，将裂纹附近的网格进行加密，所布种子长度为0.1mm；将相应的材料参数E和ν输入有限元软件中，具体数值见表1。随后，在FRANC2D中对所建模型施加荷载和约束，通过在模型下边界定义约束边界的起始点、起始点的相邻点和约束边界的截止点，能够对施加水平和竖直两个方向的约束；同时，在模型上边界施加一个均匀拉力，参考试验，控制最大力不变，输入的拉力大小分别为40 MPa(当裂纹从WCp体积分数为3%方向开始扩展)和22 MPa(当裂纹从WCp体积分数为15%方向开始扩展)。

图2 模拟计算疲劳裂纹扩展速率曲线的流程图Fig.2 Flowchart of simulating fatigue crack growth rate curves

图3 双层3%～15%(a)、三层3%～9%～15%(b)、双层15%～3%(c)和三层15%～9%～3%(d)的Cu/WCp FGMs试验疲劳试样Fig.3 Material specimen in experiment of 3%～15%bi-layered(a)，3%～9%～15% multilayer(b)，15%～3%bi-layered(c)，and 15%～9%～3% multilayer(d)

图4 双层(a)、三层(b)FGMs几何模型和Cu/WCp FGMs几何模型及有限元模型(c)Fig.4 Geometric model of bi-layered(a)，geometric model of multilayer(b)and finite element model(c)in Cu/WCp FGMs

表1 材料参数Table 1 Material parameters

3.2 分析过程

FRANC2D 对裂纹扩展行为的分析主要分为前处理、模拟计算和后处理三个步骤，有平面应力(Plane Stress)、平 面 应 变 (Plane Strain)、轴 对 称(Axisymmetric)和弯曲(Bending)四种分析类型，本研究所模拟的对象为薄板状拉伸试件，因此选择的分析类型为平面应力状态。

FRANC2D 在普通有限元的等参数单元基础上，引入奇异单元，采用标准的具有二次形函数的六节点三角形单元(Bilinear3 side)和八节点四边形单元(Bilinear4 side)。在进行网格划分时，采用裂纹尖端的K主导区适当加密网格，对裂纹尖端径向和环向的网格进行细化。这种特殊梯度单元，对各向异性材料(如FGMs)裂纹尖端的网格划分具有很大的优势，即使采用较稀疏的网格划分也能得到满足模拟要求的结果。由用户定义的裂纹增量确定新的裂纹尖端位置后，沿着增量上裂纹扩展路径的单元将被消除，并自动进行局部网格的重构。通过这种网格局部重新划分方法以及特殊网格的移动，可以进行试件疲劳裂纹扩展的模拟，如图5所示。

图5 裂纹尖端局部网格重新划分 (a)用户在初始几何构型中自定义裂纹几何形状；(b)删除新裂纹尖端位置附近的网格并进行网格更新与重构；(c)在裂纹尖端建六节点三角形单元；(d)增加裂纹尖端附近径向和环向单元数量并进行局部网格细化Fig.5 Mesh local area around the crack tip (a)user defines the crack geometry in the initial geometry；(b)mesh near the tip of the new crack is removed and the mesh is updated and reconstructed；(c)bilinear3 side are constructed at the crack tip；(d)number of radial and ring elements near the crack tip is increased and local mesh is refined

根据试验中试件的情况，在FRANC2D 中先引入缺口，长度为1mm；通过选取裂纹的中心点坐标和裂纹尖端位置，定义一个线弹性裂纹，长度为0.5mm；设置裂纹扩展量Δa为65μm，分别对不同的Cu/WCp叠层FGMs模型进行裂纹扩展，图6为双层FGMs拉伸试件的裂纹扩展路径。

图6 裂纹扩展 (a)裂纹开始扩展；(b)裂纹尖端到达WCp 含量3%(体积分数)和15%(体积分数)界面；(c)裂纹扩展结束；(d)试件变形情况Fig.6 Crack propagation (a)start propagating；(b)crack tip reaches the interfaces between two contents of WCp 3%(Vol)and 15%(Vol)；(c)stop propagating；(d)deformation of specimen

4 结果分析与对比验证

采用FRANC2D 与ABAQUS结合的数值模拟方法，计算出了Cu/WCp叠层FGMs单边缺口拉伸试件的应力强度因子值，结合Paris公式，得出了疲劳裂纹扩展速率曲线。为了对该数值模拟方法的可行性进行验证，在MTS810疲劳实验机上进行疲劳裂纹扩展试验。用电火花线切割机分别在WCp体积分数3%(或15%)一侧预制1mm 的缺口，根据规范要求[19]，预制长度约为0.5mm 的疲劳裂纹，采用长焦距显微镜M75进行裂纹长度监测。疲劳裂纹扩展过程中保持最大循环载荷Pmax不变，采用正弦波加载，频率为10 Hz，应力比为0.1，同时记录裂纹扩展量Δa i以及相邻两个裂纹尖端所对应的疲劳载荷循环次数，利用传统均质材料的计算公式来确定应力强度因子值，绘制出裂纹扩展速率曲线。研究结果显示，该方法可以快速地计算出Cu/WCp叠层FGMs的应力强度因子值，并可用于模拟叠层FGMs的疲劳裂纹扩展速率。

4.1 裂纹扩展速率和应力强度因子幅度的关系(da/d N-ΔK)

图7为数值模拟方法及试验方法下，双层及多层Cu/WCpFGMs裂纹扩展速率da/d N和应力强度因子幅度ΔK的关系对比曲线图。从图可见，梯度层界面两端的da/d N有很大的差异。由7(a)～(b)可知，双层Cu/WCpFGMs中，当控制裂纹从低WCp体积分数(3%)层向高WCp体积分数(15%)层扩展时，da/d N随着ΔK的增大而增大，在靠近界面前变化趋势有所减慢，穿过界面后迅速增大直至试件断裂；当控制裂纹从高WCp体积分数(15%)层开始扩展时，在靠近界面前da/d N变化趋势减缓甚至下降，穿过界面后迅速增大直至试件断裂。由图7(c)～(d)可知，多层Cu/WCpFGMs中裂纹扩展速率也存在和双层Cu/WCpFGMs类似的变化规律，当控制裂纹从3%(体积分数)WCp层开始扩展时，随着ΔK的增大，在靠近第一个界面(WCp体积分数为3%与WCp层与体积分数为9%)时，裂纹扩展速率减小，穿过界面后速率迅速上升；当靠近第二个界面(WCp体积分数为9%与 WCp层与体积分数为15%)时，裂纹扩展速率减小，穿过界面后迅速上升直至试件断裂；当控制裂纹从15%(体积分数)WCp层开始扩展时，随着ΔK的增大，裂纹扩展速率分别在靠近第一个界面(WCp体积分数为15%与WCp层与体积分数为9%)和第二个界面(WCp体积分数为9%与WCp层与体积分数为3%)时，都有降低的趋势。从数值模拟方法的结果中同样可以发现，无论裂纹从高WCp体积分数向低WCp体积分数扩展还是从低WCp体积分数向高WCp体积分数扩展，梯度层界面两端的da/d N都有很大的差异，裂纹扩展速率在靠近界面时均减缓或降低。因此，梯度层间界面对功能梯度材料裂纹扩展的整个阶段都存在很大影响，它可以起到减缓裂纹扩展的作用。

应用FRANC2D 与ABAQUS结合的数值模拟方法，可以更为快速的计算出FGMs在疲劳载荷下裂纹尖端的应力强度因子，并用于分析FGMs中梯度层的存在对疲劳裂纹扩展速率的影响规律。然而，数值模拟结果中梯度层界面的应力强度因子值较试验中的低，这主要是因为WCp在铜基体中的分散并不能保证是完全均匀的，试验中采用传统均质材料的计算方法得到的应力强度因子值与有限元模型建立中对材料的均质化假设，均会造成试验和模拟结果曲线的不完全吻合。

图7 裂纹扩展速率曲线(da/d N-ΔK)Fig.7 Crack growth rate curve(da/d N-ΔK)(a)3%～15%Bi-layered FGMS；(b)15%～3%Bi-layered FGMS；(c)3%～9%～15% Multilayer FGMS；(d)15%～9%～3% Multilayer FGMS

4.2 裂纹扩展速率和裂纹尖端距试件界面距离的关系(da/d N-l)

图8 裂纹扩展速率曲线(da/d N-l)Fig.8 Crack growth rate curve(da/d N-l)

图8 为数值模拟方法下，双层及多层Cu/WCpFGMs裂纹扩展速率da/d N和裂纹尖端距试件界面距离l的关系曲线。从图可见，双层梯度材料的裂纹扩展速率明显快于多层功能梯度材料的裂纹扩展速率；在双层功能梯度材料中，裂纹从高含量梯度(WCp体积分数为15%)层开始扩展，裂纹扩展速率快于裂纹从低含量梯度(WCp体积分数为3%)层开始扩展的裂纹扩展速率；在多层功能梯度材料中也存在相同的变化规律，裂纹从高含量梯度(WCp体积分数为15%)层开始扩展时，裂纹扩展速率也要稍快一些。说明梯度含量的变化对功能梯度材料的疲劳裂纹扩展影响很大；同时，梯度层的数量对功能梯度材料的疲劳裂纹扩展速率也有一定的影响。

5 结 论

本研究运用二维断裂分析有限元软件FRANC2D和ABAQUS相结合的方法，对Cu/WCp叠层FGMs的疲劳裂纹扩展行为进行数值模拟，分别对四种不同模型进行了数值分析，并与试验计算结果进行对比，得到以下结论：

1.所建立的FRANC2D 与ABAQUS 相结合的FGMs裂纹扩展模拟方法，可以较快速建立有限元模型，迅速计算出裂纹扩展中的应力强度因子，有利于比较复杂的叠层FGMs的裂纹扩展模拟研究。

2.应用该方法对Cu/WCp叠层FGMs单边缺口拉伸试件进行疲劳裂纹扩展数值模拟，得到了裂纹从不同WCp含量(体积分数为3%或15%)层开始扩展时的疲劳裂纹扩展速率曲线。发现曲线吻合度较高且变化规律及趋势基本相同，证明使用FRANC2D 计算FGMs疲劳裂纹扩展速率是可行的。

3.无论裂纹从低含量(WCp体积分数为3%)层还是从高含量层(WCp体积分数为15%)开始扩展，裂纹扩展速率曲线在界面处均明显减缓甚至下降，裂纹尖端穿过界面后，裂纹扩展速率迅速增大直至试件破坏，说明梯度变化及梯度层间界面的存在对功能梯度材料裂纹扩展的整个阶段都存在很大影响。

4.梯度含量的变化对功能梯度材料的疲劳裂纹扩展影响很大；同时，梯度层的数量对功能梯度材料的疲劳裂纹扩展速率也有一定影响。