复合型裂纹的扩展路径模拟及疲劳寿命预测

王建明，刘伟，吕鹤婷

(1.山东大学机械工程学院，山东济南250061;2.山东大学 高效洁净机械制造教育部重点实验室，山东 济南250061)

构件的固有缺陷或加工损伤会引起初始裂纹，初始微裂纹在交变载荷作用下会发生扩展最终导致结构的断裂破坏。故研究裂纹扩展规律尤其是稳态扩展过程中的裂纹扩展行为具有重要意义。

目前国内外针对裂纹扩展的数值研究方法包括有限元法、边界元法、无网格法、扩展有限元法。而有限元法因其成熟性在线弹性断裂问题中被广泛采用，其关键技术涉及裂纹尖端应力奇异场的处理以及适应裂纹扩展的网格重构，已有众多学者做过相关研究。Bittencourt等［1］基于Franc2D 软件，针对线弹性材料采用递归剖分技术进行网格重构。Xu和Needleman［2］提出了单元间裂纹传播的方法，其裂纹被假设沿着单元边界进行扩展，虽然不再需要重新划分网格但裂纹的扩展方向和扩展步长均依赖初始网格的质量，计算精度不高。Bouchard等［3］采用节点释放技术模拟了平面裂纹扩展，并对最大周向应力(MCSC)、最小应变能密度(MSEDC)、最大能量释放率(MSERRC)3种裂纹扩展准则进行了分析对比，结果表明MSEDC精确度不如另外2种高，在有限元代码方面MCSC较为简单但裂尖网格要求较密，MSERRC最复杂但所得结果较为理想。Liu等［4］基于单元子域离散与应变光滑技术提出了光滑有限元法，创建了5节点三角形奇异单元以模拟裂纹尖端应力场，其计算时间较长。Khoei等［5］提出一种在裂纹扩展过程中进行误差估计并优化网格尺寸的自适应网格划分方法，尽管提高了计算精度但算法较繁琐。王慰军［6］基于ABAQUS进行裂纹扩展二次开发，采用Deluanay三角剖分算法进行网格自动划分，每一步的裂纹扩展增量Δa的大小不对裂纹扩展方向产生影响。王永伟［7］对拉伸载荷作用下的三维表面疲劳裂纹扩展进行了模拟，获得了推导裂纹形状变化规律的有效方法。

现有研究裂纹扩展数值模型大多针对Ⅰ型裂纹，其裂纹扩展方向可事先预知。本文针对Ⅰ、Ⅱ复合型裂纹开展裂纹扩展路径数值模拟研究。基于ANSYS进行二次开发，采用位移外推法和最大周向应力准则建立了裂纹扩展路径模拟及疲劳寿命预测数值计算模型，取得了理想的仿真结果。

1 建模相关理论

1.1 应力强度因子计算

应力强度因子是描述裂尖区域位移场和应力场的重要参数，其数值计算方法主要有J积分法、虚拟裂纹闭合法、虚拟裂纹扩展法、位移外推法等。本文采用节点位移外推法计算应力强度因子。

建立裂尖坐标系如图1所示，根据线弹性断裂理论，复合型裂纹尖端附近位移场可表示为［8］

式中:u、ν分别为平行于裂纹方向(x方向)、垂直于裂纹方向(y方向)的位移分量;r、θ为以裂尖为坐标系原点的极坐标;E为弹性模量，μ为泊松比;KI、KⅡ分别为Ι型、Ⅱ型裂纹的应力强度因子;k为膨胀模量，平面应力问题取k=3-μ/1+μ，平面应变问题取k=3-4μ。

图1 裂纹尖端坐标系Fig.1 The coordinate system at the crack tip

计算KΙ时，令θ=π，由式(2)可得

由式(3)可以看出，KΙ是裂纹尖端r→0( )处的极限值，无法直接由该式得到。在此通过有限元法，利用其节点位移进行外推，结合最小二乘法数据拟合，可形成KΙ的数值分析方法。在有限元模型中，由于裂尖处应力具有奇异性，故围绕裂纹尖端设置1/4奇异单元以提高数值计算精度。

距裂尖ri处，裂纹表面张开位移νi为

式中:νa和νb分别是裂尖后端位于r=ri处的上裂纹面节点a和对应的下裂纹面节点b的y方向位移如图2所示。在裂纹面上规律选取一系列节点组成数据对:

利用最小二乘法对数据对进行拟合处理，假定KΙi和ri之间用线性关系近似:

当r=0时，。根据最小二乘法原理，其最佳拟合应满足:

类似地，

图2 裂尖后端裂纹表面张开位移示意图Fig.2 The deformation of crack surfaces behind the crack tip

1.2 裂纹扩展方向判定

由于裂纹方位不对称或者载荷分布不对称，裂纹并非简单地沿着直线向前扩展，而是根据KΙ与KⅡ的关系表现为曲线轨迹。目前裂纹扩展方向的判断主要有3种方法:1)最大周向应力准则;2)最大能量释放率准则;3)最小应变能密度准则。本文采用最大周向应力准则计算裂纹扩展方向角θ，即认为裂纹沿周向正应力σθθ最大的方向扩展。在裂尖极坐标系中，裂尖附近的应力场可表示为［8］

令 ∂σθθ/∂θ=0，τrθ=0 可得

求解上述方程可得裂纹扩展角:

Keq与材料断裂韧性Kc比较可作为裂纹是否发生失稳扩展的判据。

1.3 Paris准则

Paris用应力强度因子幅定量描述稳态区疲劳裂纹扩展速率，其表达式为

式中:C和m是实验确定的材料常数;ΔKeq是等效应力强度因子幅值，可用 ΔKΙ、ΔKⅡ代替方程(13)中的求得。利用Paris公式模拟疲劳裂纹扩展有2种处理方法:

1)在各裂纹扩展步中设定载荷循环次数ΔN，计算相应的裂纹扩展增量Δa:

2)在各裂纹扩展步中设定裂纹扩展增量Δa，计算相应的载荷循环次数ΔN:

本文采用方法2)，即在各裂纹扩展步中给定裂纹扩展增量Δa，计算相应的载荷循环次数ΔN。

2 算法流程

本文利用ANSYS平台，基于上述建模理论，通过APDL编程形成裂纹扩展路径模拟和疲劳寿命预测专用分析模块，其计算流程可分为以下4个步骤:

1)建立含裂纹的参数化几何模型;

2)网格自动划分，其中围绕裂纹尖端采用1/4奇异单元，其余部分采用二次三角形单元，参数化施加载荷及位移边界条件;

3)进行有限元计算，得到裂尖后端上下表面对应节点位移，分别利用式 ( 7)、式 ( 8)计算应力强度因子KΙ、KΙΙ，利用式(13)计算等效应力强度因子Keq;

4)通过设置多个裂纹扩展步以跟踪裂纹扩展路径。在每个裂纹扩展步中设定裂纹扩展增量Δa，计算相应的Keq(或ΔKeq)，利用式 16( )计算该扩展步对应的载荷循环次数ΔN，当Keq≥Kc时停止计算，视为裂纹稳态扩展阶段结束;否则根据裂纹扩展方向角θ及裂纹扩展增量Δa形成新的裂纹扩展步。如此循环计算直至满足Keq≥Kc。裂纹扩展增量△a取值应与成反比，Δa越小路径模拟结果越精确，但计算时间加长，一般可取初始裂纹长度的 10%～20%［9］。

图3给出该裂纹扩展路径模拟和疲劳寿命预测专用分析模块的分析流程，其中在计算裂尖位置并形成新的裂纹扩展步模块中，首先根据式(12)确定在新扩展步中的裂纹方向，再根据给定的扩展增量确定新的裂尖位置，根据新的裂尖位置调整计算网格，完成相应计算。

图3 裂纹扩展模拟流程图Fig.3 Flow chart of crack growth path prediction

3 数值算例

本节将从应力强度因子计算、复合裂纹扩展路径模拟和疲劳寿命预测三方面通过典型算例验证前文提出的算法及程序设计的正确性。

3.1 单边裂纹平板应力强度因子计算

单边裂纹平板几何模型如图4所示，裂纹长度a=0.4 mm，平板宽度b=1 mm，高度h=2 mm，轴向载荷σ=1 MPa。材料参数:弹性模量E=30 MPa，泊松比μ=0.3，按平面应变问题计算，其Ι型裂纹应力强度因子解析解为

式中:修正系数Cf取:

图5为Von Mises等效应力云图。裂尖后端5点对应力强度因子数值结果如表1所示。

将表1数据代入式(7)，得到位移外推法应力强度因子数值解为2.370 MPa·mm1/2。由式(17)可得其解析解为2.358 MPa·mm1/2，两者相对误差仅为0.5%，说明本文所给出的应力强度因子数值计算方法具有较高的精度。

图4 单边裂纹平板几何模型Fig.4 A plate with single-edge crack under tension

图5 Von Mises等效应力云图Fig.5 The contours of Von Mises stress distribution

表1 裂尖点对KΙi数值解Table 1 The KIinumerical solutions at crack tip points

3.2 有机玻璃试样裂纹扩展路径模拟

本算例相关模型参数均来源于文献［1］，以便于将仿真结果与文献［1］中的实验结果进行对照分析。选用有机玻璃PMMA试样模拟Ⅰ、Ⅱ型复合裂纹的扩展过程，同时考虑孔对裂纹扩展的影响。试样长度L=20 mm，宽度W=8 mm，3个孔的半径r=0.25 mm，位置如图6所示，将试样在底边两点处固定并在顶边中点施加集中载荷P=1 N。材料参数为:弹性模量E=3×104MPa，泊松比μ=0.3。裂纹初始长度a=1 mm、距试样左端距离b=4 mm。

图7给出了工况Ⅰ下裂纹扩展第2步、第6步时的扩展路径，分析中取裂纹扩展增量Δa=0.5 mm。由模拟结果可知初始裂纹经过8个裂纹扩展步后最终从左边趋向于中间孔，模拟路径与文献［1］实验结果对比如图8所示，二者非常接近，同时与文献［9-10］中的数值结果吻合较好。图9为工况ⅠVon Mises等效应力云图。图10、图11为裂纹扩展过程中Ⅰ、Ⅱ型应力强度因子。

图6PMMA试样几何模型Fig.6 Geometric model of the PMMA specimen

图7PMMA试样裂纹扩展路径Fig.7 Crack growth path of the PMMA specimen

图8 裂纹扩展路径模拟结果与实验结果对比Fig.8 Comparison of the crack growth path between numerical results and experimental results

图9 Von Mises等效应力云图Fig.9 The contours of Von Mises stress distribution

图10 Ⅰ型应力强度因子Fig.10 Type I stress intensity factor

图11 Ⅱ型应力强度因子Fig.11 TypeⅡstress intensity factor

3.3 单边斜裂纹疲劳寿命预测

单边斜裂纹几何模型如图12所示，长度b=100 mm，宽度h=200 mm，裂纹倾斜角度θ=40°，初始长度a0=20 mm，将模型视为平面应变问题，试样承受脉冲循环载荷，Δσ=40 MPa(σmax=40 MPa，σmin=0)，材料参数为:E=74 000 MPa，泊松比μ=0.3。Paris模型参数C=2.087 136 × 10-13，m=3.32，平面应变断裂韧性Kc=1 897 MPa·mm1/2。

图13所示为应力强度因子随裂纹长度的变化曲线图，各步裂纹扩展增量取Δa=1 mm，从图13可以看出应力强度因子KⅡ经过微幅上升之后迅速变为0，之后裂纹由Ι-Ⅱ复合型变为纯Ι型。裂纹扩展轨迹如图14所示，裂纹开裂时其方向发生突变，随后保持直线扩展直至试样失稳断裂。

图12 单边斜裂纹平板承受循环载荷Fig.12 A plate with edge-angled crack under a cyclic tension

图13 应力强度因子随裂纹长度变化曲线Fig.13 Stress intensity facors vs.crack length

图14 单边斜裂纹扩展路径Fig.14 Crack growth path of the edge-angled crack problem

图15为裂纹扩展的a-N曲线，当裂纹等效应力强度因子达到临界值Kc时，其裂纹长度为61.4 mm，载荷循环次数为130 016次，与文献［11］采用无网格法模拟得到的结果:裂纹长度61.3 mm，载荷循环次数131 411次比较，两者基本吻合，验证了本文所得方法的正确性与有效性。

图15 单边斜裂纹疲劳寿命曲线Fig.15 Fatigue life curves of the edge-angled crack problem

4 结束语

本文基于ANSYS编制了裂纹扩展路径模拟和疲劳寿命预测专用分析模块，采用位移外推法数值计算应力强度因子，并结合最大周向应力准则和Paris准则形成模拟复合型裂纹扩展路径及疲劳寿命预计的有效算法。该分析模块由APDL参数化编程，在裂纹扩展的每个计算步中采用全自动网格划分，其分析过程简便高效，为裂纹扩展模拟和疲劳寿命预测提供有利工具。

［1］BITTENCOURT T N，WAWRZYNEK P A，INGRAFFEA A R，et al.Quasi-automatic simulation of crack propagation for 2D LEFM problems［J］.Engineering Fracture Mechanics，1996，52(2):321-334.

［2］XU X P，NEEDLEMAN A.Numerical simulations of fast crack growth in brittle solids［J］.Journal of the Mechanics and Physics of Solids，1994，42(9):1397-1434.

［3］BOUCHARD P O，BAY F，CHASTEL Y.Numerical modelling of crack propagation:automatic remeshing and comparison of different criteria［J］.Computer Methods in Applied Mechanics and Engineering，2003，192(35/36):3887-3908.

［4］LIU G R，NOURBAKHSHNIA N，CHEN L，et al.A novel general formulation for singular stress field using the ESFEM method for the analysis of mixed-mode cracks［J］.International Journal of Computational Methods，2010，7(1):191-214.

［5］KHOEI A R，AZADI H，MOSLEMI H.Modeling of crack propagation via an automatic adaptive mesh refinement based on modified superconvergent patch recovery technique［J］.Engineering Fracture Mechanics，2008，75(10):2921-2945.

［6］王慰军.基于ABAQUS的裂纹扩展仿真软件及应用［D］.杭州:浙江大学，2006:17-44.WANG Weijun.The program and its applications for simulation of crack propagation based on ABAQUS［D］.Hangzhou:Zhejiang University，2006:17-44.

［7］王永伟.结构疲劳裂纹扩展的数值模拟［D］.大连:大连理工大学，2005:26-42.WANG Yongwei.Numerical simulations of fatigue crack growth of structure［D］.Dalian:Dalian University of Technology，2005:26-42.

［8］程靳，赵树山.断裂力学［M］.北京:科学出版社，2006:21-45.

［9］NOURBAKHSHNIA N，LIU G R.A quasi-static crack growth simulation based on the singular ES-FEM［J］.International Journal for Numerical Methods in Engineering，2011，88(5):473-492.

［10］XUAN H N，LIU G R，BORDAS S，et al.An adaptive singular ES-FEM for mechanics problems with singular field of arbitrary order［J］.Computer Methods in Applied Mechanics and Engineering，2013，253:252-273.

［11］DUFLOT M，DANG H N.A meshless method with enriched weight functions for fatigue crack growth［J］.International Journal for Numerical Methods in Engineering，2004，59(14):1945-1961.