用广义扩展有限元计算界面裂纹应力强度因子*

双材料被广泛应用于如航空、航天和交通运输等领域,双材料界面处往往会出现裂纹等缺陷,这些缺陷会导致材料强度的降低,而且损伤通常开始于界面或界面附近区域,因此对双材料界面裂纹的研究至关重要。

研究人员已使用如有限元法、边界元法[1]、扩展有限元法(Extended Finite Element Method,XFEM)[2-6]等数值方法计算双材料界面裂纹的应力强度因子(Stress Intensity Factor,SIF)。有限元法在处理裂纹问题时需要在裂纹尖端划分很细的网格,不可避免地会出现几何形状比较差的单元,且随着裂纹增长网格需要重构,从而影响计算精度和计算效率。XFEM[7-8]通过在单元位移场函数中增加富集函数来描述不连续问题,使得网格与结构内部细节或物理界面无关,从而克服了裂纹尖端等高应力区网格划分的困难,同时在模拟裂纹增长过程中无需对网格重构。但是XFEM在计算精度上并没有显著提高。文献[9]结合广义有限元法(Generalized Finite Element Method,GFEM)和XFEM的优点提出广义扩展有限元法(Generalized Extended Finite Element Method,GXFEM)概念,并应用于各向同性材料的裂纹分析。文献[10]分别用无网格伽辽金法和XFEM计算了双材料界面裂纹的应力强度因子,无网格伽辽金法因在其求形函数及其导数时常涉及到矩阵求逆及多个矩阵相乘,使得计算时间很长,并且其在XFEM中所用的界面富集函数往往会因为混合单元的存在,影响解的收敛性。文献[11]基于XFEM提出一种新的界面富集函数,有效提高了解的收敛性。文献[12]利用三角变换将裂纹尖端位移场富集函数由4项减少到2项,减少了计算量,但是其没有涉及到比较复杂的双材料裂纹尖端富集函数。文献[13]用常规XFEM计算双材料界面裂纹的应力强度因子,没有引入界面的富集函数,计算精度并不高。

本文首先建立了分析双材料界面裂纹的GXFEM,其中引入界面富集函数改进界面单元,应用三角变换等手段将界面裂纹的裂纹尖端富集函数由12个减少到6个,然后采用偏移富集函数对混合单元的问题进行处理,对裂纹尖端和裂纹面的单元引入附加的结点广义参数来提高计算精度。最后通过数值算例验证本文GXFEM的合理性和有效性。

对于双材料界面裂纹,GXFEM的位移场由4个部分构成：常规有限元法的位移和反映裂纹面、裂纹尖端、材料界面存在的加强位移。

式中：i、 j、 k和m为单元结点编号；G、G0、G1和G2分别为常规单元、裂纹贯穿单元(空心圆)、裂纹尖端单元(实心圆)和材料界面单元(三角形)结点的集合；g为裂纹尖端富集函数的数目；x为高斯点的坐标向量；uh为结点位移,h为空间维数；ui为结点i的常规位移；aj、bkg和cm分别为贯穿裂纹结点j、裂纹尖端结点k和材料界面结点m的附加自由度向量；H(x)、φg(x)和ψ(x)分别为裂纹面、裂纹尖端和材料界面的富集函数；Nu、 Nja、Nkb和Nmc分别为常规单元、裂纹贯穿单元、裂纹尖端单元和材料界面单元的广义形函数,上标u为常规位移,a、b和c分别为裂纹贯穿结点、裂纹尖端结点、材料界面结点的附加自由度向量。

式中：Fiu、Fja、Fbk和Fcm分别为常规单元、裂纹贯穿单元、裂纹尖端单元和材料界面单元的结点位移插值函数，由于不同结点位移的插值函数之间的协调性没有强制要求,在理论上Fui、Faj、Fbk和Fcm可以随意选取；φm、φm、φm和φm为常规有限元的形函数。

式中：r和θ为以裂纹尖端为原点的极坐标；为间断面上离x距离最短的点；f(x) 为间断面函数；ξ(x)为最短距离函数；n为间断面在点处的外法向矢量；ε为裂纹尖端双材料参数。双材料界面裂纹裂纹尖端富集函数有12项。在GXFEM中裂纹尖端结点有74个自由度,约束式(1)的bk1=bk2,bk3=bk4,bk5=bk6,bk7=bk8,bk9=bk10,bk11=bk12,利用三角变换,且不考虑常数项系数,使式(5)变为

经过变换,裂纹尖端的富集函数变为6项。式中：π/4前的符号应与θ的符号保持一致。引用式(6)作为裂纹尖端富集函数,裂纹尖端结点自由度为38个,大大减少了运算量。当ε=0时,裂纹尖端富集函数φg(x)为

双材料间的弹性不匹配可用Dundurs参数α和β来描述[13] ：

式中：μi为材料i(i=1,2)的剪切模量；κi为Kolosov常数：

式中：νi为材料i(i=1,2)的泊松比。

可以选择问题的已知解,问题本身所具有的物理特性或者部分解析解作为富集函数。文献[14]为了模拟孔洞和夹杂问题,首次在XFEM中引入材料界面富集函数：

式中：ξm(x)为材料界面富集结点m在材料界面水平集函数中的值,其表达式见式(4)。但是由于混合单元存在使得解的收敛性不好。文献[11]提出了新的改进方法,表达式为

富集单元的结点位移并不是真实的结点位移,为了使所有的结点都是真实的结点位移,需要对式(1)进行转基处理,变换为

式中：H(xj)为单元结点j的Heaviside值；φg(xk)为单元结点k的裂纹尖端单元富集函数值；ψ(xm)为单元结点m的材料界面富集函数值；xj、xk和xm分别为裂纹贯穿单元结点、裂纹尖端结点和材料界面结点坐标向量。

将式(2)代入式(14),位移模式可表示为

本文在常规结点和材料界面结点的结点位移插值函数取0阶,对应每个结点有2个自由度。结点的广义形函数为

式中：Fi0为0阶结点位移插值函数。

裂纹的贯穿单元结点和裂纹尖端单元结点的结点位移插值函数取1阶,即对应每个结点具有6个广义自由度。结点的广义形函数为

式中：Fi1为1阶结点位移插值函数。

双材料界面裂纹模型的加强方式如图 1所示。

和常规有限元法一样,将式(15)代入弹性体的虚功方程,可导出GXFEM的支配方程为

式中：δ为结点位移列阵；K和F分别为结构总刚度矩阵和结点载荷列阵,K按常规步骤由单元刚度矩阵组集而成。单元层次上的k和f分别为

式中：keij为刚度矩阵；Bir为单元应变转换矩阵；D为材料矩阵；Bjs为单元应变转换矩阵；Ω为弹性体的区域；fie为所有结点荷载力矩阵；fui为常规结点荷载力矢量；fai为裂纹贯穿结点荷载力矢量；fib1,fib2,…,fibl为裂纹尖端结点荷载力矢量；fic为材料界面结点荷载力矢量。单元结点荷载力矢量分量表示为

式中：H为裂纹贯穿单元结点的富集函数；为表面力向量；b为体力矢量；Γt为表面力给定的表面；φg为裂纹尖端富集函数。

Biu、Bia、Bib和Bic分别为常规结点、裂纹贯穿单元结点附加、裂纹尖端单元结点附加和材料界面单元附加的应变转换矩阵,其定义如式 (26)～式(29)所示：

式中：(),x和(),y分别为形函数对x求偏导和形函数对y求偏导。

双材料界面裂纹尖端的应力场的形式为[13]

式中:σ12为剪切应力；σ22为拉应力；Kf=K1+iK2为复应力强度因子,i为虚数单位；lref为参考长度。

在复杂应力状态下,已经广泛应用从相互作用积分中提取应力强度因子。对于界面裂纹,相互作用积分可以表示为

式中：δ1j为kronecker记号；ui,1aux、εikaux和σijaux分别为辅助位移对x方向求偏导、辅助应变和辅助应力；A为积分区域的面积； q,j为权函数对x、y方向求偏导；σij和σik为真实应力。

相互作用积分与应力强度因子的关系可以表示为

式中：K1aux=1,K2aux=0,此时把I记为I1,或者K1aux=0,K2aux=1,此时把I记为I2,可以求得真实场作用下的K1和K2；E的定义为

式中：E1为相对应下半平面的弹性模量；E2为相对应的上半平面的弹性模量。

σij 和ui,1由广义扩展有限元算得,ui,1aux、εikaux和σijaux由辅助位移场求出：

当计算K1时,辅助位移场函数为

当计算K2时,辅助位移场函数为

式中：C和D定义参见文献[13]。

辅助应变为

由辅助位移函数可求得

式中：

利用胡克定律,辅助应力场σijaux可由辅助应变σijaux求得。

一个二维双材料正方形板上下边界受均匀拉应力σ=1MPa,模型中心有长度为2a=2mm的裂纹,板长为Hc=2L=60mm,E1/E2=2～1000,E2=1MPa,ν1=ν2=0.3。模型是对称的,采用一半模型来模拟,其板内的界面裂纹如图 2所示。假设板处于平面应变状态,需要采取合适的边界条件消除刚体位移的影响。本问题可以看作是无限大板双材料界面裂纹问题。式(30)中l=1。令,裂纹尖端的应力强度因子K1和K2的解[15]为

本算例网格密度为151×75,计算结果列于表 1。从表 1可以看出：本文的GXFEM计算得到的应力强度因子和式(39)的结果近似,F1和F2的相对误差分别在0.2%和3%之内,说明本文的GXFEM可以准确地计算界面裂纹的应力强度因子。随着E1/E2的增大,F1和F2的绝对值随着增加。

单边界面裂纹和副裂纹板如图 3所示。一个二维弹性双材料平板上边界和下边界受单向拉伸载荷σ=1MPa,在沿界面有长度为a1=40mm的边缘主裂纹,板上还有a2=10mm副裂纹。板尺寸为L×W,L=2W=200mm,E2/E1=2～100,E2=100MPa,ν1=ν2=0.3。副裂纹的中心和主裂纹的右裂纹尖端水平距离为0,垂直移动副裂纹,计算主裂纹尖端的应力强度因子,求得的应力强度因子均通过进行无量纲化,并与文献[10]的结果进行比较。

本算例网格密度为19×39,图 4在不同的弹性模量比下,随着裂纹之间相互距离变化,主裂纹尖端无量纲应力强度因子(KⅠ、KⅡ)的变化。在同等的裂纹之间相互距离时,E2/E1=100比E2/E1=2的无量纲应力强度因子大一些。图中的无网格伽辽金法(EFGM)结果[10]和本文的GXFEM结果符合的很好。

1) GXFEM结合了XFEM和GFEM的优点,不需要在裂纹尖端设置过密的网格,同时能够保证较高的计算精度。

2) 本文在运用GXFEM处理界面裂纹时,使用三角变换将裂纹尖端位移场逼近函数由12个减少到6个,在不损失精度的情况下大大减少了计算量；采用偏移富集函数对混合单元的问题进行了处理,提高了计算的收敛性；应用GFEM对裂纹尖端和裂纹面的元素通过增加结点广义位移参数,提高了计算精度。算例表明,采用本文的GXFEM计算双材料界面裂纹的应力强度因子是成功和有效的。

3) FEM、GFEM及XFEM均可以看作GXFEM在特定条件下的简化。在处理不连续问题时,只需要对不连续的区域增加富集函数、进行结点自由度广义化,对其他区域则采用FEM进行处理,从而可以减少建模的复杂性和工作量。