MCTS试件的三维有限元计算断裂分析

李庆芬，齐桂营，朱莉，何水清

(哈尔滨工程大学机电工程学院，黑龙江哈尔滨150001)

工程构件和机器零部件中往往含有三维复杂裂纹和其他缺陷，这些裂纹和缺陷造成受力构件的应力集中、强度下降、并会导致灾难性的断裂事故.断裂力学作为研究含裂纹结构强度的一门学科，可以正确评价带裂纹构件的力学特性，其中应力强度因子作为表征外力作用下弹性体裂纹尖端附近应力场强度的一个重要参量，是评价断裂破坏和裂纹扩展规律的重要指标.因此，求解裂纹体的应力强度因子是线弹性断裂力学中很重要的一项工作，对评估含裂纹构件剩余强度和疲劳裂纹扩展寿命等都具有重要的意义.

在断裂力学发展的初期，为了数学处理上的简便，往往将实际裂纹简化为平面二维模型.目前，二维断裂问题的研究已十分成熟，已建立了一系列应力强度因子的理论计算公式［1］，并基于不同理论确定了不同的断裂判据［2-6］.然而工程实际构件中的裂纹，由于构件几何形状及承载方式的多样化，通常呈现的是三维复合型断裂模式，理想的二维裂纹几乎是不存在的.为了更好地解决工程实际中的断裂问题，研究三维复合型裂纹的断裂特性就变得尤为重要.而三维断裂问题虽已受到广大研究者的关注，但有关三维复合型断裂问题的研究还很不成熟.由于结构和裂纹的复杂性，目前只有少数典型裂纹体的应力强度因子存在理论解析解［7］，大多数三维复合型断裂问题在求解裂纹前沿应力强度因子时存在较大困难，尤其是在裂纹平面与三维构件自由表面的交汇处，迄今尚无精确的解析解.现有的三维断裂判据，如 S 判据［3］、POOK 判据［8］和最大主应力判据［9］等，也只给出了平面纯拉伸(Ⅰ型)、平面纯剪切(Ⅱ型)和反平面剪切(Ⅲ型)加载条件下的裂纹起始角和断裂极限值，而对于三维复合型裂纹，尤其是Ⅰ+Ⅱ+Ⅲ型复合裂纹的断裂行为，则迄今还少有研究［10-11］.

本文将采用由Richard设计的特殊加载设备［12］和带倾斜裂纹面(倾角γ为45°)的修正紧凑拉剪试件，通过改变裂纹的几何形状方位构建三维复合型裂纹，应用修正的虚拟裂纹闭合积分方法，对三维裂纹体复合型裂纹的应力强度因子进行分析计算，探讨其裂纹尖端的断裂行为.该研究将为深入探讨三维复合型断裂问题做出有益的贡献，并为解决实际工程中存在类似裂纹构件的安全可靠性评估提供参考依据.

1 数学模型与计算方法简介

1.1 裂纹尖端应力场及应力场强度因子

裂纹尖端应力场是建立断裂判据的理论依据.直角坐标系下裂纹尖端应力场如下:

式中:σx、σy、σz、τxy、τxz、τyz为直角坐标系下各方向的应力分量，r、θ为裂纹尖端附近点的极坐标.其中，应力强度因子Ki(i=Ⅰ、Ⅱ、Ⅲ)唯一地描述了裂纹尖端应力场的强度，是评价断裂行为的重要参数.因此，分析计算裂纹前沿的应力强度因子的分布情况是研究裂纹体断裂特性的基础.

1.2 修正的虚拟裂纹闭合积分方法

Rybicki和Kanninen在虚拟裂纹闭合积分方法(virtual crack closure integral method，VCCIM)［13］的基础上提出了修正的虚拟裂纹闭合积分方法(modified virtualcrack closure integralmethod，MVCCI)［14］.该方法用裂纹长度为a的原始裂纹尖端节点位移ΔUy，i-1(a)代替虚拟裂纹闭合积分方法所要求扩展裂纹的节点位移ΔUy，j-1(a+Δa)，避免了当裂纹扩展为a+Δa时，二次有限元分析的附加影响(见图1).具体表达式如下:

式(2)适用于承受平面纯拉伸载荷即二维张开型(Ⅰ型)裂纹模式，图1(a)为二维模型裂纹尖端的有限元网格.式中:t为试件厚度，F为节点力，ΔUy为节点的相对位移，Wy为节点力作用产生节点位移所做的功，G1c(a)为二维裂纹尖端的能量释放率.

图1 MVCCI计算方法的低阶单元有限元网格Fig.1 FE low order standard elements for MVCCI-method

MVCCI这种高精度的解析方法还可以直接结合低阶的实体单元来计算复杂的三维断裂问题.图1(b)为三维裂纹前沿的有限元网格.具体的三维MVCCI解析计算公式如下:

MVCCI方法是一种近似的方法，其值收敛于精确解.该方法已被用于计算典型的二维和三维断裂问题［15］，并取得了较好的计算精度.

2 带倾斜裂纹面的MCTS试件模型及有限元分析

2.1 MCTS试件及加载设备

本文选用带倾斜裂纹面(倾角γ为45°)的修正紧凑拉剪试件(modified compact tension shear，MCTS).试件及其特定的加载设备如图2所示.

裂纹平面的倾斜导致裂纹方位关于加载方向的不对称分布，形成复合型裂纹，裂纹在不同加载条件下将会产生不同的复合型断裂模式.随加载角度α的改变，试件将从承受纯拉伸(α=0°)加载，拉伸/剪切(0°＜α＜90°)复合型加载，变化到纯剪切(α=90°)加载.

图2 带倾斜裂纹面的MCTS试件及其特定加载设备Fig.2 MCTS-specimen with inclined crack plane and the special loading device

试件的主要几何参数如下:宽度W=90 mm，高度h=145 mm，厚度t=10 mm，裂纹长度a=50 mm，裂纹面倾角γ为45°.所选材料参数为:弹性模量E=80 kN/mm2，泊松比 ν =0.3，材料为铝合金，试件所承受的静态载荷F=40 kN.

2.2 有限元模型的建立

采用HyperMesh软件提供的模型变形技术对已有的紧凑拉剪试件(compact tension shear specimen，CTS)［16］的计算模型进行改进，通过倾斜其裂纹平面构建具有复合型裂纹的MCTS计算分析模型，揭示裂纹面方位对断裂行为的影响.为了形成倾斜裂纹平面并保持各个网格之间的关联关系，建模过程中将对CTS试样上下1/2有限元模型分别建立修改域，修改域中的操纵点沿总体坐标系y轴上下移动，使裂纹平面所在的整个横截面发生倾斜，然后通过节点耦合模拟倾斜裂纹平面的直裂纹前沿，并合并非裂纹面区域的双重节点，最终形成倾斜裂纹平面MCTS试样的有限元模型.由于上述数值计算方法无需用奇异单元来模拟裂纹尖端的奇异性，所以整个三维有限元模型均采用八节点的solid45单元，并对裂纹前沿附近进行了网格细化，以保证计算精度，最精细(Δa=0.25 mm)的三维有限元模型共4 104个单元，4 886个节点.

图3给出了带倾斜裂纹面MCTS试件的整体有限元模型及裂纹前沿区域网格.

图3 带倾斜裂纹面MCTS试件的有限元模型及网格Fig.3 FE-model and mesh of the CTS-specimen with inclined crack plane

2.3 有限元计算方法

本文采用商用ANSYS软件分别对Ⅰ型、Ⅱ型、及Ⅰ+Ⅲ复合型(α=45°)加载条件下的MCTS试样进行了有限元计算.

首先根据MVCCI方法，计算得出三维裂纹体沿裂纹前沿各个节点处的能量释放率(strain energy release rate，SERR)Gi(a，z/t)，i=Ⅰ、Ⅱ、Ⅲ.然后应用Irwin关系式［13］将能量释放率转换为应力强度因子如下:

将结果标准化如下:

式中:σ0=F/wt，σ0为裂纹面正应力.

最后，绘制出具有倾斜裂纹平面的MCTS试样直裂纹前沿的应力强度因子分布曲线.进而探讨倾斜裂纹平面MCTS试样在不同加载条件下的断裂行为.

3 计算结果与分析

3.1 Ⅰ型载荷作用下的计算结果与分析

Ⅰ型载荷作为脆性断裂最危险的加载方式一直是研究重点.图4给出了I型载荷(α=0°)作用下带倾斜裂纹面(γ=45°)和直裂纹面(γ=0°)试样裂纹前沿的应力强度因子的分布曲线.由图可知:直裂纹面CTS的裂纹前沿上只有Ⅰ型应力强度因子的存在;而MCTS试样的裂纹前沿上却同时存在着Ⅰ型、Ⅱ型和Ⅲ型应力强度因子.说明在纯张开型载荷作用下，直裂纹面试样只能发生张开型断裂，而倾斜裂纹面试样不仅发生了张开型断裂，同时还产生了滑开型和撕开型断裂，即产生了Ⅰ+Ⅱ+Ⅲ复合型断裂模式.说明I型载荷下，倾斜裂纹面试样的断裂模式不再满足单一载荷与断裂模式唯一对应的关系.

图4 Ⅰ型载荷作用下MCTS试样裂纹前沿的标准化SIFFig.4 Normalized SIF along crack front of MCTS-specimen under mode-Ⅰloading

MCTS试样在I型加载下的变形情况如图5(a)所示.由图可看出，试件整体变形以张开型断裂为主，但可以清楚地看到倾斜裂纹平面同时发生了滑开型和撕开型位移.CTS试样的外形尺寸及加载条件与MCTS试样完全相同，而MCTS试样裂纹平面的倾斜，导致裂纹方位关于加载方向的不对称分布，形成了复合型裂纹，即倾斜裂纹平面所形成的复杂裂纹几何形状是产生三维复合型断裂模式的主要原因.

由计算结果还可看出，Ⅱ型应力强度因子在z/t=0.5处具有最大值，由此可知试件将在两侧自由表面附近首先发生失稳扩展.

图5 MCTS试样在不同加载下的变形情况Fig.5 Deformed FE-model of the MCTS-specimen under different loading conditions

3.2 Ⅱ型载荷作用下的计算结果与分析

图6给出了Ⅱ型载荷(α =90°)作用下带倾斜裂纹面(γ=45°)MCTS试样裂纹前沿的应力强度因子的分布曲线.

由图6可知:带倾斜裂纹面MCTS试样的裂纹前沿上虽然仍以Ⅱ型应力强度因子为主，但在试样两侧自由表面附近区域可以看到少量Ⅰ型和Ⅲ型应力强度因子的存在.说明在纯滑开型载荷作用下，倾斜裂纹面试样不仅发生了滑开型断裂，同时还产生了张开型和撕开型断裂，即产生了Ⅰ+Ⅱ+Ⅲ复合型断裂模式.这与直裂纹面试样在纯Ⅱ型载荷下往往只发生滑开型断裂模式是不同的.由计算结果还可看出，应力强度因子在两侧自由表面附近具有最大值，由此可知试件将在两侧自由表面附近首先发生失稳扩展.

图6 Ⅱ型载荷作用下MCTS试样裂纹前沿的标准化SIFFig.6 Normalized SIF along crack front of MCTS-specimen under mode-Ⅱloading

MCTS试样在 II型加载下的变形情况，如图5(b)所示.由图可看出，试件整体变形以滑开型断裂为主，但也可以清楚的看到倾斜裂纹平面同时发生了张开型和撕开型位移.

3.3 复合型载荷作用下作用下的计算结果与分析

图7给出了复合型载荷(α=45°)作用下带倾斜裂纹面(γ=45°)MCTS试样裂纹前沿的应力强度因子的分布曲线.由图可知:在Ⅰ+Ⅱ复合型载荷作用下，裂纹前沿产生了Ⅰ+Ⅱ+Ⅲ复合型断裂.裂纹前沿强度因子的分布曲线关于中心平面不对称，在z/t=0.5处应力强度因子具有最大值，由此可知自由表面附近最易发生裂纹的失稳扩展.

图7 复合型载荷作用下MCTS试样裂纹前沿的标准化SIFFig.7 Normalized SIF along crack front of MCTS-specimen under mixed-mode loading

MCTS试样在Ⅰ+Ⅱ复合型加载下的变形情况如图5(c)所示.从图7中也可看出，试件上下倾斜裂纹平面在Ⅰ+Ⅱ复合型载荷作用下同时产生了张开型、滑开型和撕开型位移.

上述数值分析结果表明:带倾斜裂纹平面的MCTS试样在不同加载条件下，其裂纹前沿均产生了Ⅰ+Ⅱ+Ⅲ复合型变形模式，即裂纹失稳扩展的临界条件将受到3种主要断裂模式共同作用的影响，其不再沿原裂纹平面扩展.

裂纹平面倾斜导致的裂纹方位不对称及断裂模式之间的耦合作用是产生复合型断裂模式的主要原因.根据Richard判据，该复合型裂纹在3种基本断裂模式的共同影响下，其裂纹前沿在初始扩展阶段发生偏折的同时也产生了扭转变形，断裂破坏的不稳定区域则随载荷的变化而变化，MCTS试样在I型或Ⅱ型载荷作用下将在两侧自由表面附近首先发生失稳扩展，而对于平面内任意复合型加载条件，z/t=0.5处自由表面附近区域是该复合型裂纹断裂破坏的危险区域.

纵观不同加载条件下裂纹前沿应力强度因子的分布曲线，I型载荷作用时裂纹前沿各节点应力强度因子值均大于其他载荷情况，由此可知，倾斜裂纹平面MCTS试样在平面纯拉伸载荷作用下最易发生失稳扩展，说明I型载荷对于三维复合型裂纹来说仍是最具破坏的加载模式.

4 实验结果

为了对数值计算结果进行验证，采用透明有机玻璃脆性材料(便于观察裂纹开裂情况)加工了带有45°倾斜平面裂纹的MCTS试样，并分别在Ⅰ型和Ⅱ型加载模式下进行了实验(见图8、9).

由实验结果可知，在I型加载模式下(图8)，裂纹在初始扩展阶段同时发生扭转和偏折变形，形成一个新的“螺旋”状的裂纹面，随着裂纹的不断扩展，扭转的裂纹面逐渐汇合成一个新的光滑平面，新裂纹面与加载方向垂直，裂纹沿这个新的裂纹表面方向扩展直至断裂，即倾斜裂纹平面MCTS试样最终发生了张开型断裂.这与倾斜平面裂纹在初始变形阶段，试件整体变形以张开型断裂为主，但同时发生了滑开型和撕开型断裂的计算结果是一致的.

图8 在Ⅰ型加载下的MCTS试样Fig.8 MCTS-specimen under mode Ⅰ loading

在Ⅱ型加载模式下(图9)，裂纹面在初始扩展阶段发生偏折的同时也产生了轻微的扭转变形，扭转变形在试样两侧自由表面附近变形幅度较大，然而随着扩展的深入，扭转现象逐渐消失，裂纹最终发生了滑开型断裂.这与裂纹在初始变形阶段，试件整体变形以滑开型断裂为主，但也同时发生了张开型和撕开型位移的计算结果是一致的.

图9 在Ⅱ型加载下的MCTS试样Fig.9 MCTS-specimen under mode Ⅱ loading

根据对比分析，实验中试样的裂纹开裂扩展特征与本文的计算结果是一致的，一定程度上验证了本文采用MVCCI数值计算方法所做的计算分析.

5 结论

本文对带有倾斜裂纹面的MCTS试件的应力强度因子进行了有限元分析及实验验证.研究结果表明:

1)试件在不同加载条件下均产生了Ⅰ+Ⅱ+Ⅲ复合型断裂，裂纹平面倾斜导致的裂纹方位不对称及断裂模式之间的耦合作用是产生三维复合型断裂模式的主要原因.

2)实验中试样的裂纹开裂扩展特征与有限元分析的结果一致，一定程度上验证了本文的计算分析.

本文的研究拓展了三维断裂研究的空间，为深入探讨三维复合型断裂问题做出了有益的贡献，并为解决实际工程中存在类似裂纹构件的安全可靠性评估提供了参考依据.

［1］中国航空研究院.应力强度因子手册［M］.北京:科学出版社，1993:7-438.Aeronautical Academy of China.Academy handbook of stress intensity factors［M］.Beijing:Science Press，1993:7-438.

［2］ERDOGAN F，SIH G C.On crack extension in plates under plane loading and transverse shear［J］.Trans ASME，Journal of Basic Engineering，1963，85:519-527.

［3］SIH G C.Strain energy density factor applied to mixed mode crack problems［J］.Int J Fracture，1974，10(3):305-321.

［4］陈增涛，王铎.复合型裂纹的混合开裂准则［J］.哈尔滨工业大学学报，1994，26(4):127-130.CHEN Zengtao，WANG Duo.A mixed fracture criterion for mixed mode crack propagation［J］.Journal of Harbin Institute of Technology，1994，26(4):127-130.

［5］崔德渝，吴绍富，张行.复合型裂纹的能量释放率计算与断裂判据［J］.北京航空航天大学学报，1990，3:43-48.CUI Deyu，WU Shaofu，ZHANG Xing.Energy release rate and failure criterion for mixed mode cracking［J］.Journal of Beijing University of Aeronautics and Astronautics，1990，3:43-48.

［6］蒋玉川，胡兴福.复合型裂纹扩展的最大塑性区尺度准则［J］.四川大学学报，2005，37(4):11-14.JIANG Yuchuan，HU Xingfu.The fracture criterion of the maximum size of plastic zone of mixed mode crack growth［J］.Journal of Sichuan University，2005，37(4):11-14.

［7］范天佑.断裂理论基础［M］.科学出版社，2003:138-186.FAN Tianyou.Foundation of fracture theory［M］.Beijing:Science Press，2003:138-186.

［8］POOK L P.A finite element analysis of the angle crack specimen［M］.London:Mechanical Engineering Publications，1993:285-302.

［9］SCHOLLMANN M，RICHARD H A，KULLMER G，FULLAND M.A new criterion for the prediction of crack development in multiaxially loaded structures［J］.Int J Fract，2002，117(2):129-141.

［10］QI Guiying，BUCHHOLZ F G，YAN Shengyuan，et al.Computational fracture analyses of a compact tension shear(CTS)specimen with an inclined crack front［J］.Key Engineering Materials，2010，417-418:629-632.

［11］CITARELLA R，CRICRÌ G.Comparison of DBEM and FEM crack path predictions in a notched shaft under torsion ［J］.Engineering Fracture Mechanics，2010，77(11):1730-1749.

［12］RICHARD H A，BENITZ K.A loading device for the creation of mixed mode in fracture mechanics［J］.Int J Fracture，22，1983:R55-R58.

［13］IRWIN G R.Analysis of stresses and strains near the end of a crack traversing a plate［J］.J of Applied Mechanics，1957，24:361-364.

［14］RYBICKI E F，KANNINEN M F.A finite element calculation of stress intensity factors by a modified crack closure integral［J］.Engng Fract Mech，1977，9:931-938.

［15］FRIEDRICH-G B，AZZEDINE C，HANS A R.Computational fracture analyses by means of virtual crack closure integral methods［C］//CD-ROM-Proceedings of the 6th Congreso Argentino de Mecanica Computational(MECOM99).Mendoza，Argentina，1999:29.1-29.19.

［16］RICHARD H A.A new compact shear specimen［J］.Int J Fracture，1981，17:R105-R107.