﻿ 基于扩展有限元的疲劳裂纹扩展分析
 舰船科学技术  2022, Vol. 44 Issue (20): 29-34    DOI: 10.3404/j.issn.1672-7649.2022.20.005 PDF

Fatigue crack propagation analysis based on extended finite element method
ZHANG Gong-xue, SHI Shao-bin, JIN Hui-jin, LIU Miao
School of Mechanical and Electrical Engineering, Shanxi University of Science and Technology, Xi′an 710021, China
Abstract: For the crack propagation problem of structural plate, based on the extended finite element method of Abaqus, three methods were used to calculate the static stress intensity factor of the specimens with unilateral crack and the central crack respectively. The results showed that the interaction integral method was more accurate and the error was within 2%. The interaction integral method was adopted to calculate the dynamic stress intensity factor during the crack propagation process, and the error was less than 3%. Combined with the low cycle fatigue analysis, the residual life of the specimen with hole edge crack was predicted, and the error was 5.9% by comparing with the test results. In addition, the test results in the stable stage of fatigue crack propagation are basically consistent with the finite element analysis, which proved the accuracy and feasibility of this method. Therefore, this method could provide reliable reference for calculating dynamic stress intensity factor and predicting residual life during the process of fatigue crack propagation.
Key words: extended finite element     stress intensity factor     interaction integration method     residual life
0 引　言

1999年Belytschko等[2]提出扩展有限元法，通过将富集函数引入常规有限元来反映计算区域的不连续性，使裂纹尖端的间断与网格独立，降低了裂纹尖端网格的划分难度，可以有效解决疲劳断裂问题。而且相较于疲劳试验，扩展有限元法应用灵活，周期短，能够模拟多种工况下的疲劳断裂，大大降低了试验成本，已发展成为解决疲劳裂纹扩展问题最主要的方法之一[3-5]。Shi[6]用扩展有限元法模拟岩石变形，使用考虑高阶项的位移外推法计算应力强度因子。Zou[7]分析了等效积分区域法不同积分路径对应力强度因子的影响。马梓聪[8]基于扩展有限元使用等效区域积分法计算了不同裂宽比的应力强度因子。余天堂[9]用相互作用积分法计算应力强度因子，之后，相互作用积分法得到广泛应用[10-12]。还有很多学者基于扩展有限元，对疲劳裂纹的剩余寿命进行预测[13-14]，为以后更好地理解和预测疲劳过程及寿命提供依据。

1 裂纹体的扩展有限元方法

1999年Moes[15]用阶跃函数加强单元节点，体现了裂纹的不连续特征。近似的广义位移场形式为：

 $\begin{split}{u^h}(x) =& \sum\limits_{i \in {N_s}} {{u_i}{N_i}(x)} + \sum\limits_{j \in {N_{cut}}}{{a_j}{N_j}(x)H(x)} + \\ &\sum\limits_{k \in {N_{tip}}} {{N_k}(x)\sum\limits_{\alpha = 1}^4 {{b_{k\alpha }}{B_\alpha }(x)} }。\end{split}$ (1)

 $H(x) = \left\{ \begin{array}{ll} +1 & \psi (x,t) > 0 ，\\ - 1 & \varphi (x,t) < 0 。\end{array} \right.$ (2)

Bα(x)计算公式为：

 $\begin{split} [B_{\alpha} (x),\alpha =&1,\cdots \text{，}4]=\\ &\Biggr[\sqrt{r}\text{sin}\frac{\theta }{2}\text{,}\sqrt{r}\text{cos}\frac{\theta }{2}\text{,}\sqrt{r}\text{sin}\frac{\theta }{2}\mathrm{sin}\theta ,\sqrt{r}\text{cos}\frac{\theta }{2}\text{sin}\theta \Biggr]。\end{split}$ (3)

 图 1 加强节点分布 Fig. 1 Strengthen node distribution
 $\begin{split}{{{u}}^{{h}}}(x) = &\sum\limits_{i \in {N_s}} {{u_i}{N_i}(x)} + \sum\limits_{j \in {N_{cut}}}{N_j}(x)(H(x) - {H_j}(x)){a_j} + \\ &\sum\limits_{k \in {N_{tip}}} {{N_k}(x)\sum\limits_{\alpha = 1}^4 {({B_\alpha }(x) - {B_{k\alpha }}(x)){b_{k\alpha }}} }。\end{split}$ (4)
2 应力强度因子计算方法及结果分析 2.1 基于节点位移的外推法

 ${K_{\rm I}} = \mathop {\lim }\limits_{r \to 0} [{\sigma _y}(r,\theta = 0)\sqrt {2\text{π} r} ] 。$ (5)

 ${K_{Ii}} = \frac{{2\mu }}{{\kappa + 1}}{v_i}\sqrt {\frac{{2\text{π} }}{{{r_i}}}}，$ (6)

 图 2 基于位移的外推法计算应力强度因子 Fig. 2 Displacement based extrapolation to calculate stress intensity factor
 ${K_I} \approx \frac{{\displaystyle\sum {{r_i}} \sum {{K_{Ii}}} - N\sum {{r_i}{K_{Ii}}} }}{{{{\left( {\displaystyle\sum {{r_i}} } \right)}^2} - N{{\displaystyle\sum {{r_i}} }^2}}}。$ (7)
2.2 等效积分区域法

 $J = \int\limits_\varGamma {W{\rm{d}}y - {p_i}\frac{{\partial {u_i}}}{{\partial x}}} {\rm{d}}s。$ (8)

 $W = \frac{1}{2}{\sigma _{{{ij}}}}{\varepsilon _{{{ij}}}} 。$ (9)

 图 3 积分区域 Fig. 3 Integral area
 $J = {\smallint _s}\left[\left({\sigma _{xx}}\frac{{\partial u}}{{\partial x}} + {\tau _{xy}}\frac{{\partial v}}{{\partial x}} - W\right)\frac{{\partial q}}{{\partial x}} + \left({\tau _{xy}}\frac{{\partial u}}{{\partial x}} + {\sigma _{yy}}\frac{{\partial v}}{{\partial x}}\right)\frac{{\partial q}}{{\partial y}}\right]{\rm{d}}S 。$ (10)

2.3 相互作用积分法

 ${I^{(1,2)}} = \int_\varGamma {\left[{W^{(1,2)}}{\delta _{1,j}} - {\sigma _{ij}}^{(1)}\frac{{\partial {u_i}^2}}{{\partial {x_1}}} - {\sigma _{ij}}^{(2)}\frac{{\partial {u_i}^1}}{{\partial {x_1}}}\right]} {n_j}{\rm{d}}\varGamma 。$ (11)

 ${W^{(1,2)}}{\text{ = }}\frac{1}{2}({\sigma _{ij}}^{(1)}{\varepsilon _{ij}}^{(2)} + {\sigma _{ij}}^{(2)}{\varepsilon _{ij}}^{(1)})。$ (12)

 $I = \frac{{2({K_I}^{(1)}{K_I}^{(2)} + {K_{{\rm I}{\rm I}}}^{(1)}{K_{{\rm I}{\rm I}}}^{(2)})}}{{E{\text{*}}}}。$ (13)

 $E^*=\left\{\begin{array}{cc}E& 平面应力，\\ \dfrac{E}{1-{\text{v}}^{\text{2}}}&平面应变。\end{array}\right.$ (14)

 ${K_I} = \frac{{E^*}}{2}{I^{(1,\bmod el{\kern 1pt} I)}}，$ (15)
 ${K_{II}} = \frac{{E^*}}{2}{I^{(1,\bmod el{\kern 1pt} II)}} 。$ (16)
2.4 静态应力强度因子计算结果

 图 4 应力强度因子计算模型 Fig. 4 Stress intensity factor calculation models
 $\begin{split} K =& \sigma \sqrt {\text{π} a} \Biggr[1.12 - 0.23\frac{a}{B} + 10.6{\Biggr(\frac{a}{B}\Biggr)^2} -\\ &21.7{\Biggr(\frac{a}{B}\Biggr)^3} + 30.4{\Biggr(\frac{a}{B}\Biggr)^4}\Biggr](a/B < 0.6)，\end{split}$ (17)
 $K=\dfrac{\Biggr(\sigma \sqrt{\text{π} a}\Biggr[1-0.5\Biggr(\dfrac{2a}{B}\Biggr)+0.326{\Biggr(\dfrac{2a}{B}\Biggr)}^{2}\Biggr]\Biggr)}{\sqrt{\Biggr(1-\dfrac{2a}{B}}\Biggr)}(2a/B \leqslant 0.7) 。$ (18)

2.5 裂纹扩展过程动态应力强度因子计算结果

3 疲劳裂纹扩展过程及剩余寿命分析

 $\frac{{{\rm{d}}a}}{{{\rm{d}}N}} = C{(\vartriangle K)^m}。$ (19)

 $f = \frac{N}{{{c_1}\vartriangle {G^{{{\text{c}}_2}}}}} \geqslant 1.0。$ (20)

 $\frac{{{\rm{d}}a}}{{{\rm{d}}N}} = {c_3}{(\vartriangle G)^{{{\text{c}}_{\text{4}}}}}。$ (21)

 图 5 疲劳裂纹扩展试样模型 Fig. 5 Fatigue crack growth specimen model

 图 6 疲劳裂纹扩展过程应力云图 Fig. 6 Stress cloud diagram of fatigue crack propagation process

 图 7 疲劳裂纹扩展速率预测 Fig. 7 Prediction of fatigue crack growth rate

 图 8 疲劳裂纹扩展a/N曲线 Fig. 8 Fatigue crack growth a/N curve
4 结　语

1）在Abaqus中基于扩展有限元法建立了疲劳裂纹扩展模型，通过Python脚本实现了静态裂纹和裂纹扩展的动态应力强度因子计算。

2）分别使用位移外推法、等效积分区域法和相互作用积分法计算单边裂纹和中心裂纹的应力强度因子，发现相互作用积分法计算结果最稳定、相对误差最小，在2%以内。

3）使用相互作用积分法计算了单边裂纹扩展过程的动态应力强度因子，计算误差在3%以内。

4）计算了疲劳裂纹扩展过程的动态应力强度因子，结合Paris公式，使用低周疲劳分析，对孔边裂纹的剩余寿命进行预测，与试验结果相比，误差为3.8%，并且疲劳裂纹扩展稳定阶段的试验结果与扩展有限元分析基本一致。

 [1] 底月兰, 王海斗, 董丽虹, 等. 扩展有限元法在裂纹扩展问题中的应用[J]. 材料导报, 2017, 31(3): 70-74+85. DI Yuelan, WANG Haidou, DONG Lihong, et al. Application of the extended finite element method in crack propagation[J]. Materials Reports, 2017, 31(3): 70-74+85. DOI:10.11896/j.issn.1005-023X.2017.03.012 [2] T, Belytschko, T, Black. Elastic crack growth in finite elements with minimal remeshing[J]. International Journal for Numerical Methods in Engineering, 1999, 45: 602-620. [3] FENG S Z, LI W. An accurate and efficient algorithm for the simulation of fatigue crack growth based on XFEM and combined approximations[J]. Applied Mathematical Modelling, 2017, 55(MAR.): 600-615. [4] NAGASHIMA T, WANG C. XFEM analyses using two-dimensional quadrilateral elements enriched with only the heaviside step function[J]. International Journal of Computational Methods, 2022, 19(2): .2150063. DOI:10.1142/S0219876221500638 [5] 薛河, 王双, 王正, 等. 基于扩展有限元法的三维裂纹前缘应力强度因子计算方法[J]. 舰船科学技术, 2022, 44(3): 1-5. XU He, WANG Shuang, WANG Zheng, et al. XFEM-based calculation mothod for stress intensity factor of three-dimensional crack front[J]. Ship Science And Technology, 2022, 44(3): 1-5. DOI:10.3404/j.issn.1672-7649.2022.03.001 [6] SHI F, LIU J. A fully coupled hydromechanical XFEM model for the simulation of 3D non-planar fluid-driven fracture propagation[J]. Computers and Geotechnics, 2021, 132(5): 103971. [7] ZOU G, HE C. Path-dependent J-integrals under mixed-mode loads of mode I and mode II[J]. Theoretical and Applied Fracture Mechanics, 2018, 96: 380-386. DOI:10.1016/j.tafmec.2018.05.014 [8] 马梓聪. 含裂纹海洋平台结构剩余疲劳强度评估方法研究[D]. 大连: 大连理工大学, 2018. [9] 余天堂. 含裂纹体的数值模拟[J]. 岩石力学与工程学报, 2005, 24(24): 4434–4439. [10] 汪必升, 李毅波, 廖雅诗, 等. 基于扩展有限元模型的动态应力强度因子计算[J]. 中国机械工程, 2019, 30(11): 1294–1301. [11] MH A, BNN A, MTA B, et al. A new numerical approach in the seismic failure analysis of concrete gravity dams using extended finite element method. Engineering Failure Analysis. 2021. 132: 105835. [12] OT A, HO B, YY C. A New approach of the interaction integral method with correction terms for cracks in three-dimensional functionally graded materials using the tetrahedral finite element[J]. Theoretical and Applied Fracture Mechanics, 2021, 103065. [13] GRBOVIC A, KASTRATOVIC G, SEDMAK A, et al. Fatigue crack paths in light aircraft wing spars[J]. International Journal of Fatigue, 2019, 123(JUN.): 96-104. [14] LI Qiang, CHEN Zhixiang, LUO Shaohui. Reflection crack analysis of asphalt overlay on cement concrete pavement based on XFEM[J]. Advances in Civil Engineering. 2021, 1230447. [15] MOES N. A finite element method for crack growth without remeshing[J]. International Journal for Numerical Methods in Engineering, 1999, 46: 131-150. DOI:10.1002/(SICI)1097-0207(19990910)46:1<131::AID-NME726>3.0.CO;2-J [16] 夏世林, 磨季云. 基于J积分的含裂纹加筋板在单轴拉伸载荷下的极限强度计算[J]. 舰船科学技术, 2021, 43(7): 34–39, 59. XIA Shi-lin. MO Ji-yun. Ultimate strength calculation of cracked stiffened-panel under loading of uniaxial tension using J-integral theory[J]. Ship Science and Technology. 2021, 43(7): 34–39, 59. [17] BENVENUTI E. An effective XFEM with equivalent eigenstrain for stress intensity factors of homogeneous plates[J]. Computer Methods in Applied Mechanics & Engineering, 2017, 321(1): 427-454. [18] PARIS P, ERDOGAN F. A critical analysis of crack propagation laws[J]. Journal of Basic Engineering, 1963, 85(4). [19] 马颖涵, 张振, 郭孟雨, 等. F690超高强钢的腐蚀疲劳裂纹扩展行为及其有限元模拟[J]. 机械工程材料, 2021(4). 65–17, 80. [20] JIN Huijin, WU Sujun. A new driving force parameter for fatigue growth of multiple cracks[J]. International Journal of Fatigue, 2017, 10-16.