位移外推法和扩展有限元法计算应力强度因子的比较

2024-03-20张亚洲钟红王立强

科学技术与工程 2024年5期

张亚洲, 钟红, 王立强

(中国水利水电科学研究院, 北京 100048)

混凝土作为非均质材料,其表面和内部存在大量微裂缝甚至宏观缺陷[1]。混凝土坝中的裂缝可能损害大坝的整体性,改变其受力状态,降低其耐腐蚀性,从而缩短其使用寿命[2-4]。对混凝土抗裂性能的研究有利于延缓、减轻混凝土的开裂以及对已开裂混凝土采取止裂措施[5],从而提高混凝土结构的耐久性和安全性[6]。其中混凝土材料中存在的微裂缝可能影响裂尖应力强度因子。作为判断结构失效和裂缝的扩展,高效准确地计算裂尖应力强度因子尤为重要。

自1957年Irwin提出应力强度因子(stress intensity factor,SIF)概念后,该值成为断裂力学的重要参量,引起广大学者的研究[7-8]。常用的计算方法有根据定义求解的外推法和间接法,间接法利用应力强度因子的相关参数求解,例如由Rybicki等[9]提出的虚拟裂缝闭合法,该法通过计算虚拟裂缝后的节点力和节点位移从而避免求解裂缝尖端复杂的应力场和位移场。Rice[10]提出的J积分法,该法与积分路径无关,因此对于复杂的裂缝尖端可以转换成其他路径进行求解,大大降低难度,并且该方法不仅适用线弹性模型[11],在弹塑性断裂分析中也可广泛运用,Abaqus有限元软件已集成基于J积分的扩展有限元法(extended finite element method,XFEM)计算应力强度因子。Song等[12]基于比例边界有限元法进行应力强度因子的计算,该方法精度较高,是由于其相似中心会尽可能趋近于裂缝尖端位置,并且该方法不仅适用于均质材料,对于界面材料同样适用。文献[13-14]也通过有限元法对裂尖应力强度因子进行计算。此外,学者们采用声发射[15]、光弹性贴片[16]试验与有限元法相结合的方式计算裂尖应力强度因子。

上述研究对裂尖应力强度因子的求解做出了突出贡献,但是目前尚未有使用两种不同的方法对同一模型进行求解,并与解析解或文献解进行对比的研究。基于此,现通过4个经典案例,分别使用位移外推法和扩展有限元法计算裂尖应力强度因子,并检验两种方法的精度。

1 位移外推法计算应力强度因子

外推法最先由Chan等[17]利用应力场和位移场中奇异项进行展开,基于最小二乘法绘制K与r曲线拟合所得应力强度因子,公式如下。

(1)

(2)

2 扩展有限元法计算应力强度因子

相比于传统有限元法,XFEM以单位分解为理论基础并在模型中裂纹是可以穿过网格,如图1所示。XFEM常通过J积分法来计算应力强度因子,J积分在裂缝尖端周围的曲线曲面积分存在路径无关情况,其断裂韧度指标以能量形式呈现,而相互积分法作为其扩展可以合理、自然地描述总能量的释放[18],因此可以运用在扩展有限元法中。通过引入辅助场[式(3)～式(5)中用上标(2)表示]与真实场[式(3)～式(5)中用上标(1)表示]进行叠加得到J积分分量,称为相互作用积分,其表示式如式(3)所示。而J为J(1)、J(2)和M的和,因此J的表达式如式(5)所示。

图1 扩展有限元法的各单元示意图

(3)

(4)

(5)

3 实例计算

3.1 三点弯曲梁模型

含有预制裂缝的三点弯曲梁模型示意图如图2所示,其中长L=800 mm、跨度S=600 mm、高度H=200 mm、厚度B=60 mm,预制裂缝长度a0=20、40、60、80和100 mm,即缝高比α范围为0.1～0.5,跨高比β=3。模型上表面中间施加F=4 kN压力荷载,下部左端施加水平向和竖直向边界条件,右端施加竖直向边界条件。材料属性为:弹性模量E=30 GPa,泊松比υ=0.167。

图2 三点弯曲梁模型

应力强度因子的计算公式[19]如下。

(6)

(7)

(8)

3.1.1 位移外推法

在Abaqus中建立二维模型用于位移外推法计算应力强度因子,裂尖网格加密处理,考虑平面应变状态设置网格为CEP4单元。

以缝高比为0.2的模型对裂缝周围进行4种数量网格细化,分别包括50、100、150和200个,得到的结果如表1所示,计算相对误差,即解析解KI减去位移外推计算得到KI的绝对值与解析解KI的比值,发现随网格密度的增加,其相对误差越小。根据解德等[20]实例中计算应力强度因子时需要将裂缝尖端奇异点去掉可以提高精确度,以200个网格为例,绘制Ki-r曲线,根据公式法计算出的KI=0.512 MPa·m1/2。图3中没有去除奇异点所得到的KI=0.502 MPa·m1/2,相对误差为1.95%。发现裂缝尖端和末端拟合效果较差,因此去除裂尖奇异点所得到的KI=0.511 MPa·m1/2,相对误差为0.16%,其决定系数R2=0.999,拟合直线效果更好。因此使用位移外推法时可通过进行对模型细化网格和在结果处理中去除裂缝尖端和末端的部分数据来提高结果的准确性。

表1 缝高比为0.2的模型不同网格数量结果

图3 Ki-r曲线图

基于上述要求,计算不同缝高比下运用位移外推法计算的应力强度因子如表2所示,发现应力强度因子随缝高比的增大而增加,且每组数据的相对误差均小于1%。

表2 位移外推法下不同缝高比应力强度因子结果

3.1.2 扩展有限元法

在Abaqus中建立三维模型,并制定裂缝类型为XFEM裂缝,设置裂缝不发生扩展并且围道积分数目设置为5个,全局采用的网格单元类型为C3D8R。

以缝高比为0.2的模型为例,讨论网格密度对应力强度因子的影响:分别进行了尺寸为10、5和2.5 mm正方形网格,得到的结果如表3所示,发现结果相对误差均在0.5%以内且网格细化可以略提高结果精度,但是对于三维模型,过于细化网格会导致计算速度降低,网格尺寸为2.5 mm计算时间是5 mm的14倍,因此综合精度和时间考虑,在计算三维模型时可以对网格适当加密。

表3 缝高比为0.2的三点弯曲梁不同网格尺寸结果

不同缝高比下运用扩展有限元法计算出的应力强度因子如表4所示,发现该方法计算应力强度因子的相对误差保持在1%内。

表4 XFEM下不同缝高比应力强度因子结果

3.2 矩形板中心裂缝

含有中心预制裂缝的矩形薄板模型如图4所示,其中高2H=20 m、宽度2W=10 m、厚度B=0.5 m,预制裂缝长度2a=1、2、3、4、5和6 m,即缝宽比α范围为0.1～0.6,模型上、下表面施加σ=1 MPa拉应力,材料属性为:弹性模量E=30 GPa,泊松比υ=0.167。

图4 矩形板中心裂缝模型

根据《应力强度因子手册》[21],得到该模型下左、右裂尖的计算公式如下。

(9)

(10)

表5 不同缝宽比下无量纲应力强度因子结果

表6 不同倾角下无量纲应力强度因子结果

3.3 矩形板中心斜裂缝

含有中心预制裂缝的矩形薄板模型示意图如图5所示,其中模型高2H=20 m,宽度2W=10 m,厚度B=0.5 m,为尽量减小边界对裂缝的影响,取预制裂缝长度2a=1 m,即缝宽比为0.1,裂缝与水平方向倾角β,分别取20°、30°、45°、60°和70°。模型上、下表面施加σ=1 MPa拉应力,材料属性为:弹性模量E=30 GPa,泊松比υ=0.167。

图5 矩形板中心斜裂缝模型

图6 矩形板中心界面裂缝模型示意图

根据《应力强度因子手册》[21],得到该模型上、下裂尖的计算公式如下。

(11)

式(11)中:KI和KII分别为Ⅰ型和Ⅱ型应力强度因子,MPa·m1/2。

3.4 矩形板中心界面裂缝

含有中心预制界面裂缝的矩形薄板模型示意图如图11所示,其中模型高2H=20 m、宽度2W=10 m、厚度B=10 m,预制裂缝长度2a=1.5和5 m,即缝宽比为0.15和0.5。模型上、下表面施加σ=1 MPa拉应力。将模型上下平分给定不同弹性模量E1和E2来形成两个材料,其中的E1/E2分别取2、5和10,泊松比υ1=υ2=0.3。

界面材料的应力强度因子目前没有解析解,比例边界有限元法在计算应力强度因子时得到广泛运用[22-25]。通过文献[22]的模型和结果来与位移外推法和XFEM得到的应力强度因子进行对比。

表7 位移外推法下界面应力强度因子结果

使用位移外推法法得到的结果如表7所示,发现K1的值远大于K2,表明张开型应力强度因子占主导因素。因此运用扩展有限元法计算应力强度因子时只提取K1。将K1无量纲化处理结果如表8所示,发现当两种材料的弹性模量相差较小时,K1与文献[22]结果较为接近。但是随模量比值的增大,其相对误差有增大的趋势。这可能是由于Abaqus集成的XFEM计算应力强度因子是根据M积分,而M积分在计算过程中依赖弹性模量。当周围的弹性模量变化过大可能导致计算出现较大误差,从而与文献和外推方法的结果相差较大。与位移外推法相比,使用扩展有限元法计算得到的应力强度因子略大。