基于XFEM的主次裂纹间应力强度因子相互作用

陈小翠，杜成斌，江守燕

(河海大学力学与材料学院，江苏南京 210098)

尽管有限元法已存在数十年，但是有限元法仍不能有效地模拟工程中的裂纹及裂纹生长问题。对于裂纹问题，有限元法需布置高密度网格，并且在裂纹生长计算时有限元网格需要重剖分［1］。扩展有限元法(extended finite element method，XFEM)就是基于有限元处理不连续问题的弊端而提出的一种新的数值计算方法［2］。自1999年美国西北大学Belytschko教授及其研究组［3］提出扩展有限元法后，因XFEM在计算不连续问题时的优越性，十多年内得到迅速发展:Belytschko等［4］在XFEM中引入新的开裂准则，来判断裂纹生长路径和速度;Song等［5］增强了含裂纹单元的描述;Fries［6］、Cheng等［7］在常规 XFEM 基础上提出了Corrected-XFEM，以提高数值求解精度;Liu等［8］通过谱单元与XFEM结合，改善动态裂纹扩展的数值扰动问题。在发展过程中，XFEM的理论基础在不断更新进步，以提高数值求解的精度。

笔者采用文献［9］提出的新型裂尖改进函数，计算分析次裂纹的位置及长度对主裂纹应力强度因子的影响。改进的裂尖改进函数在保留了裂纹尖端场的应力奇异性和裂纹上、下表面位移不连续性的基础上，减少了裂尖改进单元的附加自由度。

1 Reduced XFEM的位移模式简介

式中:x——考察点的坐标;I——求解域中所有结点的集合;——被裂纹完全贯穿的单元结点集合;——裂尖改进结点集合;Ni(x)——结点i处的常规有限元插值形函数;(x)——单位分解函数，其形式可以与Ni(x)相同，也可以不同，文中取二者相同;H(x)——Heaviside改进函数;ui——常规有限元部分在结点i处的未知量;ai——与Heaviside改进相关的结点未知量;(x)——裂尖单元结点改进函数——与裂尖改进相关的结点未知量。

二维裂纹体的位移场用常规XFEM可表示为［10］

常规扩展有限元方法根据线弹性断裂力学裂尖解析位移场的基本形式，裂尖单元结点的改进函数可用式(2)表示［10］:

式中:r、θ——裂尖极坐标。

式(3)中π/4前的符号与θ保持一致。

采用文献［9］中裂尖改进函数，裂尖改进节点的自由度个数从8个减少为4个，并且该裂尖改进函数仍保留了裂纹尖端区应力奇异性项)和位移不连续性(sin(*)函数项)。XFEM位移模式可表示为

2 水平集法

2.1 裂纹的水平集法描述

裂纹用水平集函数 ψ(x，t)和 φk(x，t)(k=1，2)描述［11］，裂纹面的位置通过 ψ(x，t)的零水平集函数ψ(x，t)=0来描述。用符号距离函数来构造水平集函数，ψ(x，t)可表示为

波前水平集 φk(x，t)(k=1，2)与 ψ(x，t)正交，可表示为

式中:x*——考察点P在裂纹面上的投影点坐标;xk——第k个裂缝尖端的坐标;n——裂纹面的单位外法向;t——第k个裂纹尖端处的单位切向矢量;sign(x)——符号函数，x＞0时sign(x)=1，x=0时sign(x)=0，x＜0时sign(x)= －1。

2.2 改进单元类型的判断

应用XFEM计算裂纹问题时，共有裂尖改进单元和Heaviside改进单元这2类改进单元。文中分析时先计算出每个单元内各节点到裂尖的值(φk(x，t))i及 ψi(x，t)，再比较得到各个单元的 φmax、φmin和 ψmax、ψmin。裂尖改进单元为单元节点的水平集函数满足如下条件:

Heaviside改进单元为单元节点的水平集函数满足如下条件:

3 应力强度因子的计算

应力强度因子是衡量裂纹尖端区应力场强度的重要参数，也是断裂力学中裂纹的失效判据。应力强度因子计算方法主要是解析解法和数值解法。文中采用相互作用积分法计算应力强度因子，裂纹尖端相互作用能量积分公式为［12］

其中

裂纹尖端相互作用积分与应力强度因子的关系为

取XFEM计算的数值解作为真实场，线弹性断裂力学裂尖解析场为辅助场(应力场、位移场)，由式(11)可得真实场的Ⅰ型应力强度因子:

4 数值算例

4.1 改进的XF EM计算KⅠ

本算例考察有限尺寸拉伸板的裂纹问题，如图1所示。拉伸板宽b=1 m，高h=2 m，弹性模量E=1 MPa，泊松比υ=0.3，轴向拉伸应力σ=1 kPa。网格划分为19×39(共741个)均匀网格，如图1(c)所示。计算时分别取缝长 a为0.1 m、0.2 m、0.3 m、0.4 m、0.5 m。

边缘裂纹问题应力强度因子的解析解为［13］

中心裂纹问题应力强度因子的解析解为［14］

图1 有限尺寸拉伸板Fig.1 A finite stretched plate

采用改进的XFEM来计算图1(a)和图1(b)所示裂纹问题的KⅠ，并与解析解的计算结果进行比较，如图2所示。由图2可知，改进的XFEM方法的计算结果与解析解结果吻合较好。

图2 解析解与XFEM数值解比较Fig.2 Comparison of analytical and XFEM numerical results

4.2 有限尺寸拉伸板同侧裂纹问题

本算例分析有限尺寸拉伸板两同侧裂纹之间的影响，计算的拉伸板宽b=1 m，高h=2 m，弹性模量E=1 MPa，泊松比υ=0.3，轴向拉伸应力σ=1 kPa。板边缘设2条裂纹C1和C2，位置如图3所示。为研究C2的位置和长度对C1的应力强度因子的影响，保持C1长度a1=0.5 m不变，分别改变C2的长度a2及2条裂纹之间的距离d，用改进的XFEM来计算C1相应的应力强度因子值。

计算结果如图4和图5所示，C1单裂纹情况下KⅠ由文献［15］公式计算为3.545 kPa·m1/2。由图4和图5可知，由于C2存在，C1的KⅠ明显比单裂纹状态结果(3.545 kPa·m1/2)小;由图4可知，C1长度a1不变，随着C2的长度a2从0.1 m递增到0.6 m，C1的KⅠ呈明显递减趋势;当a2/a1=0.9左右时，C2长度增量对C1的KⅠ影响最大;2条裂纹间的距离d对C1的KⅠ也有一定影响，如图5所示，C1的应力强度因子随裂纹间距离增加呈“勺”形分布，且这种趋势在a2/a1≥0.86时消失，呈递增趋势;从图5左右两端差异可知，右端无限向单裂纹状态下应力强度因子接近，这主要是由于当C2无限远离C1时，C1更接近单裂纹状态;左端裂纹越长，C1的应力强度因子相对单裂纹的减少量越大。

图3 拉伸板Fig.3 A stretched plate

图4 C2长度对C1的KⅠ影响Fig.4 Effect of length of C2 on KⅠ of C1

图5 裂纹间距离对C1的KⅠ影响Fig.5 Effect of crack distance on KⅠ of C1

4.3 有限尺寸拉伸板异侧裂纹问题

本算例分析有限尺寸拉伸板的异侧2条裂纹之间的影响，仍采用4.2节的算例，并在板的左右设2条裂纹，位置如图6所示。为研究C2的位置和长度对C1的应力强度因子的影响，保持C1长度a1=0.4 m不变，分别改变C2的长度a2及2条裂纹之间距离d，再用改进的XFEM来计算C1相应的应力强度因子值。

计算结果如图7所示，C1长度a1不变，随着裂纹间距离d增大，C2长度变化对C1的KⅠ影响减小，而随着C2的长度a2从0.1 m递增到0.4 m，C1的KⅠ值呈明显递减趋势。

图6 有限尺寸拉伸板(2条裂纹)Fig.6 A finite stretched plate(two cracks)

5 结 语

用文献［9］提出的改进XFEM进行次裂纹对主裂纹应力强度因子的影响研究，通过断裂力学中有限尺寸板的算例，得出该XFEM计算的应力强度因子精度较高，与解析解结果吻合。在此基础上研究次裂纹位置和长度对主裂纹应力强度因子的影响，结果表明:对于平行裂纹问题，只要结构有另一条裂纹存在，主裂纹的应力强度因子会明显比单裂纹状态下的小;同侧裂纹情况，主裂纹的应力强度因子随次裂纹长度增加有明显的递减趋势，当d/a1达到1.5时主裂纹的应力强度因子变化趋于平缓，主裂纹的应力强度因子随裂纹间距离增大呈“勺”形，且这种趋势在a2/a1≥0.86时消失，呈递增趋势;对于异侧裂纹情况，主裂纹的应力强度因子随次裂纹长度的增加呈整体递减趋势，但随着2条裂纹间距离的增大，次裂纹长度变化对主裂纹应力强度因子的影响减小。

图7 C2长度和位置对C1的影响Fig.7 Effects of length and position of C2 on C1

［1］李录贤，王铁军.扩展有限元法(XFEM)及其应用［J］.力学进展，2005(1):5-20.(LI Luxian，WANGTiejun.The extended finite element method and its applications［J］.Advances in Mechanics，2005(1):5-20.(in Chinese))

［2］应宗权，杜成斌，程丽.含夹杂非均质材料的扩展有限元数值模拟［J］.河海大学学报:自然科学版，2008，36(4):546-549.(YING Zongquan，DU Chengbin，CHENG Li.Application of extended finite element method in heterogeneous materials with inclusions［J］.Journal of Hohai University:Natural Sciences，2008，36(4):546-549.(in Chinese))

［3］BELYTSCHKO T，Black T.Elastic crack growth in finite elements with minimal remeshing［J］.Int J Numer Meth Engng，1999，45(5):601-620.

［4］ BELYTSCHKO T，CHEN Hao，XU Jingxiao，et al.Dynamic crack propagation based on loss of hyperbolicity and a new discontinuous enrichment［J］.Int J Numer Meth Engng，2003，58(12):1873-1905.

［5］SONGJ，BELYTSCHKOT.A method for dynamic crack and shear band propagation with phantom nodes［J］.Int JNumer Meth Engng，2006，67:863-893.

［6］FRIES T.A corrected XFEM approximation without problems in blending elements［J］.International Journal for Numerical Methods in Engineering，2008，75(5):503-532.

［7］CHENGK W，FRIEST.Higher-order XFEM for curved strong and weak discontinuities［J］.International Journal for Numerical Methods in Engineering，2010，82(5):564-590.

［8］LIU Z L，MENOUILLARD T，BELYTSCHKO T.An XFEM/Spectral element method for dynamic crack propagation［J］.Int J of Fracture，2011，169(2):183-198.

［9］江守燕，杜成斌.一种扩展有限元断裂分析的裂尖单元新型改进函数［J］.力学学报，2013，45(1):134-138.(JIANG Shouyan，DU Chengbin.A novel enriched function of elements containing crack tip for fracture analysis in the framework of extended finite element methods［J］.Chinese Journal of Theoretical and Applied Mechanics，2013，45(1):134-138.(in Chinese))

［10］NICOLAS M，JOHN D，BELYTSCHKO T.A finite element method for crack growth without remeshing［J］.International Journal for Numerical Methods in Engineering，1999，46(1):131-150.

［11］DUFLOT M.A study of the representation of cracks with level sets［J］.2007，70(11):1261-1302.

［12］TOSHIO N，YOUHEI O，SHUICHI T.Stress intensity factor analysis of interface cracks using XFEM［J］.International Journal for Numerical Methods in Engineering，2003，56:1151-1173.

［13］赵建生.断裂力学及断裂物理［M］.武汉:华中科技大学出版社，2003.

［14］ ZAMANI A，GRACIE R，ESLAMI M R.Cohesive and non-cohesive fracture by higher-order enrichment of XFEM［J］.International Journal for Numerical Methods in Engineering，2012，90:452-483.

［15］MOHAMMADI S.Extended finite element method for fracture analysis of structures［M］.Oxford，UK:Blackwell Publishing Ltd，2008.