裂纹面分布加载裂尖SIFs分析的广义参数Williams单元*

徐 华， 曹 政， 邹云鹏， 杨绿峰

（1.广西大学 土木建筑工程学院 工程防灾与结构安全教育部重点实验室，南宁 530004；2.广西大学 土木建筑工程学院 广西防灾减灾与工程安全重点实验室，南宁 530004；3.广西桂能工程咨询集团有限公司，南宁 530015）

引 言

工程结构通常是带裂缝工作的，如水文环境中的混凝土坝、海洋工程、地下工程和压力容器等，其结构表面的裂缝易遭受流体侵入，有可能加速裂缝扩展与连通，引起结构的局部强度和整体承载力降低，甚至引发结构失效破坏，这些重大工程一旦破坏将引起难以估量的损失.因此，对裂纹面受荷的工程结构断裂行为研究具有十分重要的现实意义.

裂纹面受荷时，其裂尖附近应力场和位移场较复杂，难以定量描述，通常仍以SIFs作为裂尖附近应力-应变场强弱程度的度量，并可据此进一步判断结构的安全性.基于此，国内外学者围绕裂纹面加载的裂尖SIFs开展了大量研究工作.刘钧玉等[1-2]、Zhong等[3]和陈白斌等[4]基于比例边界有限元法，分别对裂纹面受拉伸荷载的单边裂纹、多裂纹和界面裂纹的裂尖SIFs进行了分析，该方法相较于普通有限单元法，需引入多个系数矩阵，增加了计算工作量.李亚等[5]采用线场分析方法，对无限大板中心裂纹面受均布荷载的应力场进行修正，进而得到了有限宽板中心裂纹的裂尖SIFs，但不同的裂纹模型均需相对应的应力函数.Fett等[6]利用ABAQUS有限元软件拟合半无限大板边界斜裂纹面上受集中力的权函数，以此求解了裂尖复合型SIFs.文献[5-6]针对不同的裂纹模型，均需找到相对应的函数才能求解SIFs，不具有通用性.Li等[7]采用柔度法求解有限宽板边界裂纹面受分布荷载的裂尖SIFs，需通过叠加法对原始模型进行拆分处理，不能直接求解.Walters等[8]针对三维模型中裂纹面受荷的平面裂纹SIFs，基于相互作用积分法分离出独立的裂纹面荷载积分项，提高了计算精度.Muthu等[9]提出了一种基于扩展无网格Galerkin法的裂纹闭合积分理论，通过求解应变能释放率获得裂尖SIFs的方法.在实际工程方面，贾金生等[10]基于断裂力学理论估算了重力坝坝踵水平裂缝中含高压水的裂尖SIFs，但基本假设较为粗略，与实际结果存在较大误差.唐世斌等[11]建立了在土压力、裂纹面和井筒内水压耦合作用下的射孔裂尖SIFs计算方程.综上所述，当前的裂纹面受荷问题仍以有限单元法为主，但该方法本身仍存在较大的改进空间.广义参数有限元法[12]是有限单元法的一个分支，该方法通过整体控制方程求解出与SIFs直接相关的广义参数，避免了线性外推拟合等繁琐的后处理过程引起的二次误差，从而具有很高的精度.但该方法建立的W单元需满足奇异区内裂纹面自由的边界条件，即 σφ=0，τρφ=0(φ=±π)，而裂纹面受荷时不满足该边界条件，仍需改进.

针对同样不满足边界条件的曲线裂纹[13]，课题组在裂尖建立等效区，将等效区内的曲线转换为折线以满足W单元的边界条件，受此思想启发，对裂纹面直接受荷模型，本文提出了一种基于SIFs互等的裂纹面裂尖奇异区荷载等效处理方法，即将奇异区内的分布荷载等效为奇异区边界的集中荷载，避免奇异区裂纹面受荷而无法使用W单元求解，充分发挥了W单元求解SIFs的高精高效优势.同时，通过算例分析，给定荷载类型等效系数建议值，为裂纹面受荷求解裂尖SIFs提供了新的思路.

1 裂尖局部分布荷载的等效处理

以一含边界裂纹的矩形板为例，其板高为2H，板宽为W，裂纹长度为a，在裂纹张口处建立整体坐标系XOY，取裂尖局部等效区边长为2c，裂纹面受到分布压力 σ(X)=σ0|1-X/a|λ， σ0为O点处压力值，如图1所示.

图1 裂尖局部面荷载等效Fig.1 Crack tip local surface load equivalence

笔者尝试推导整个裂纹面加载时裂尖应力场和位移场的Williams级数表达式，即裂纹面为非自由面，此时，因边界条件较难处理，其表达式异常复杂.根据《应力强度因子手册》[14]，裂尖局部加载和裂纹面上作用集中力情况均有高精度的参考解.基于此，现提出利用SIFs互等，即将裂尖局部分布力等效为裂纹面上一对集中力Q，如图1所示.以λ=0为例，裂尖局部均布荷载σ0和裂纹面集中荷载Q作用时，其裂尖Ⅰ型SIFs表达式为

式中，KI，分别表示仅裂尖局部受分布压力和裂纹面受一对集中压力的Ⅰ型SIFs；p，p′分别表示裂尖局部加载和裂纹面作用集中力时的荷载系数，详见文献[14].

假定KI=则集中力Q与裂尖局部荷载σ0有如下关系：

由此，可得到裂纹面局部受均布荷载的等效公式.同理，可以推导裂纹面裂尖局部长度c上受分布荷载σ (X)的等效公式，即

式中，P表示等效荷载系数，

2 基于W单元确定裂纹面受荷的SIFs

2.1 裂纹面奇异区受荷的W单元模型

如图1所示，绿色正方形区域为裂尖奇异区，其外围为常规区，因W单元对奇异区尺寸不敏感，假定本文选取的奇异区尺寸等于上节中的裂尖局部尺寸.奇异区内采用W单元，外围常规区选用四边形8节点等参单元；常规区网格可通过ANSYS有限元软件进行自由离散，奇异区的网格划分规则如下：将奇异区均分为8个三角形条元，各条元离散为n个基于径向离散比例因子α自相似的梯形子单元和一个裂尖三角形微单元，因裂尖三角形微单元极小，可忽略其刚度贡献，则将内部n个位移场由Williams级数控制的子单元集成的梯形条元定义为W单元.W单元最外层子单元因其一条边（图2中的绿色边界）为奇异区与常规区的交界，保留边界上3个实节点，内部5个节点转化为虚节点，称该单元为过渡单元.奇异区内分布荷载σ (X)可通过等效处理为奇异区边界一对集中力Q，如图2所示.

图2 裂尖奇异区网格划分及分布荷载等效Fig.2 Discretization and distributed load equivalence in the crack tip singular region

以裂尖为局部坐标原点o，建立局部直角坐标系xoy，如图2所示.以Williams级数表示裂尖奇异区位移场，并取级数的前m+ 1项：

式中，u，v分别表示裂尖局部直角坐标系下x，y轴对应的位移分量；ρ，φ分别表示裂尖局部极坐标系下的极径和极角分量，fl,11(φ)，fl,12(φ)，fl,21(φ)，fl,22(φ)为三角函数，具体取值可参考文献[12].

2.2 奇异区刚度方程集成

以任一W单元为例，根据四边形8节点等参单元理论，该条元内任意第k层子单元的刚度方程可表示为

结合式(4)，将任意第k层子单元位移列阵表示为矩阵形式：

式中，T(k)表示第k层子单元转换矩阵；ϕ表示广义参数列阵.

将式(5)等号两边同时左乘T(k)T得

并可改写为

式中，K(k)，p(k)分 别为第k层子单元的广义刚度矩阵和广义荷载列阵，且当k=2 ～n时，p(k)=0.

根据式(11)可知，各层子单元刚度矩阵大小相同，且存在相同乘数ϕ，故可将条元内第2 ～n层子单元刚度方程进行叠加：

即Kssϕ=0，Kss表示第2～n层子单元广义刚度矩阵之和.

奇异区最外层的过渡单元，即W单元中第1层子单元，有3个实节点位于常规区与奇异区交界处，5个虚节点位于奇异区内部.根据节点所在区域不同，将过渡单元的位移向量分块表示为

将过渡单元的广义刚度方程按照虚、实节点分块表示，并与第2～n层子单元的刚度方程进行叠加，因此可得一个W单元的刚度方程：

根据式(14)，将8个W单元刚度方程进行集成，并将奇异区边界的裂纹面等效荷载Q叠加至相应节点可得到奇异区的整体刚度方程：

2.3 整体刚度方程的集成

将常规区所有单元的刚度方程进行集成，并以奇异区外围节点为界，分块表示为

计算模型分为常规区和奇异区两个部分，即将式(15)和(16)的刚度方程进行集成，得到整个模型的刚度方程：

2.4 SIFs求解

通过引入模型的应力、位移边界条件，对整体刚度方程式（17）进行求解，即可获得裂尖奇异区的广义参数列阵ϕ.将列阵中的广义参数a1，b1分别代入式(18)，计算获得裂尖应力强度因子：

3 算例分析

3.1 裂纹面受均布压力

例1如图1所示的边界裂纹矩形板，其板高为2H，宽为W，且有 W=2H=100cm，裂纹长度为a，板厚t取单位厚度，奇异区范围取边长为2c的正方形.裂纹面施加均布压力（ λ=0）， σ0=1kN/cm2， 弹性模量E=3×104MPa ，Poisson比ν =0.3.因篇幅所限，本文先研究此方法在Ⅰ型裂纹中的适用性，Ⅱ型裂纹另做研究.

取 a /W=0.3的 模型为例，根据前期研究成果建议奇异区尺寸取 c =a/20，使用ANSYS有限元软件对模型常规区进行网格自由离散，此时离散为146个单元，501个节点，如图3所示.将常规区的单元、节点信息以及约束等导入自编程的FORTRAN程序，通过求解模型整体刚度方程以获得裂尖SIFs.P=1.94+0.57(a/W)-1.51(a/W)2.本模型参数的具体取值如表1所示.

表1 模型参数（λ=0）Table 1 Model parameters (λ=0)

图3 常规区网格离散Fig.3 Discretization in the regular region

根据《应力强度因子手册》（P303、P304）的裂尖局部面荷载作用下荷载系数p，p′随a/W变化曲线图[14]，通过GetData软件提取荷载系数p，p′代表性点的坐标，利用Excel软件拟合即可获得P与a/W的关系式：

由表1可以看出，随着a/W的变化，等效荷载系数P在1.60～2.00范围取值.鉴于此，本算例适当放大P的取值范围，研究在P取不同值时，a/W的变化对SIFs的影响.

W单元的3个重要参数的建议值[12]为 α =0.9， m =20， n =300.将W单元计算结果与文献[14]的解对比，为直观体现各计算结果的区别，将SIFs无量纲化如图4所示.

图4 边界裂纹面受均布压力计算结果Fig.4 Calculation results for uniformly distributed pressure on the crack surface

均布压力作用于边界裂纹面时，若忽略奇异区裂纹面的荷载（P=0），计算误差为4%～17%，且随着a/W的增大而减小，表明随着裂纹长度的增加，裂尖局部荷载对SIFs的影响逐渐变小.对奇异区裂纹面荷载等效处理后，现选取1%为误差限，根据叠加法可知，当P的取值范围为[2.0，2.1]时，任意a/W均满足该误差限.整体上，随着a/W的增加，无量纲SIFs呈增大的趋势，且增大的趋势愈加显著.

例2当均布压力作用于中心裂纹上，如图5所示，其板高为2H，宽为2W，且有W=H=100cm，矩形板厚度t为单位厚度，裂纹长度为2a，奇异区边长2c，材料参数同本算例的边界裂纹板.

图5 中心裂纹面受均布荷载Fig.5 The uniformly distributed load on the central

根据对称性，取半模型进行分析，以a/W=0.3的 模型为例，根据前期研究成果建议奇异区尺寸取c=a/20，使用ANSYS有限元软件对模型常规区进行网格自由离散，此时离散为181个单元，602个节点，如图6所示.将常规区的单元、节点信息以及约束等导入自编程的FORTRAN程序，通过求解模型整体刚度方程以获得裂尖SIFs.

图6 常规区网格离散Fig.6 Discretization in the regular region

根据《应力强度因子手册》[14]拟合出中心裂纹的等效荷载公式P=2.03+0.11(a/W)+0.11(a/W)2.当a/W取值在[0.01，0.8]区间时，P的取值范围为[2.03，2.18].适当放大P的取值范围，研究P取不同值的情况下，a/W的变化对SIFs的影响.

W单元的3个重要参数的建议值[12]为 α =0.9，m=20，n=300.将W单元计算结果与文献[14]的解对比，如图7所示.

图7 中心裂纹面受均布压力计算结果Fig.7 Calculation results for uniformly distributed pressure on the central crack surface

均布压力作用于中心裂纹面时，若忽略奇异区裂纹面的荷载（P=0），计算误差均在11%以上，最大达到20.2%，且该误差随着a/W的增大而逐渐减小.同样以1%为误差限，根据叠加原理可以推断出，对任意a/W，均满足该误差限要求的P取值范围为[1.95, 2.05].随着中心裂纹板的裂纹长度增加，无量纲SIFs呈递增的趋势，且递增幅度逐渐变大.

综上，对奇异区荷载进行等效处理是有必要的，在1%的误差限下，现给出P建议取值为2.0，不仅对裂纹面受均布压力的边界裂纹和中心裂纹都有良好的适用性，且保证了计算精度.均布荷载和线性荷载是工程中常见的荷载类型，故下文将以裂纹面受线性压力为例，验证P建议取值的通用性.

3.2 裂纹面受线性压力

算例各参数具体取值与模型网格划分同本文的例1和例2，仅荷载形式不一样，此时施加于裂纹面的分布压力为σ (X)=σ0|1-X/a|λ，σ0=1kN/cm2，λ =1.结合本文的荷载等效处理方法和W单元求解算例中的SIFs.

将常规区的单元、节点信息以及约束等导入自编程的FORTRAN程序，通过求解模型整体刚度方程以获得裂尖SIFs.本算例中，W单元重要参数建议值为α =0.9，m=20，n=300.

① 裂纹面受到线性压力时，由于文献[14]无荷载系数p，p′的具体取值，故无法拟合出等效荷载系数P的关系式.现尝试P=2.0的建议取值时，研究不同a/W模型下的边界裂纹和中心裂纹SIFs，并与文献[14]的参考解进行比较.无量纲SIFs的计算结果如图8所示.

由图8可知，当P=2.0时，计算结果误差均满足1%的误差限，证明在分析裂纹面受线性压力的中心裂纹或边界裂纹问题时，荷载等效的处理方法均有很好的通用性和高精性.当裂纹长度较小时，边界裂纹与中心裂纹的无量纲SIFs接近，而随着a/W的增大，无量纲SIFs呈增大的趋势，且边界裂纹的增长幅度大于中心裂纹.

图8 裂纹面受线性压力计算结果Fig.8 Calculation results for linear pressure on the crack surface

② 当P=2.0 时 ，现取边界裂纹和中心裂纹的矩形板宽为W=100cm，研究不同的H/W情况下，a/W的变化对SIFs的影响，并与文献[14]的参考解进行比较.无量纲SIFs结果如图9所示.

由图9可看出，H/W和a/W取不同值时，W单元计算结果均满足1%的误差限，这说明在奇异区尺寸取建议值的情况下，等效荷载系数P对裂纹位置或模型尺寸都不敏感，且适用于常见的裂纹面分布荷载类型，证明了本文的荷载等效处理方法具有很强的适用性，但由于篇幅所限，算例中只对均布荷载与线性荷载进行了论证.当H/W增大时，无量纲SIFs呈减小的趋势；当a/W增大时，无量纲SIFs呈增大的趋势，这表明裂纹模型是存在尺寸效应的.

图9 无量纲SIFs随尺寸变化的计算结果Fig.9 Calculation results for dimensionless SIFs varying with sizes

4 结 论

本文对裂纹面奇异区荷载进行等效处理，使得等效奇异区内裂纹面满足W单元裂纹面自由的边界条件，通过裂纹面受不同荷载的边界裂纹和中心裂纹奇异区荷载等效算例分析，验证了本文等效处理的正确性与通用性，得到了以下结论：

1) 确定等效奇异区尺寸后，结合算例简化了等效荷载系数P的关系式，并建议取值为2.0.算例结果表明，本文方法与《应力强度因子手册》对比，其SIFs误差均小于1%；而且《应力强度因子手册》所采用的权函数法，当模型、荷载变化时，均需重新确定权函数，相较于W单元缺乏灵活性.

2) 等效荷载系数P的取值不再是由模型尺寸、裂纹长度决定的关系式，无需通过《应力强度因子手册》拟合，且适用于裂纹面受均布或线性压力的边界裂纹和中心裂纹，证明了本文对奇异区荷载进行等效处理方法具有很强的简便性与通用性.

3) 含裂纹的矩形板，裂纹面受荷时，随着a/W的增大，Ⅰ型无量纲SIFs均呈增大的趋势，且增大幅度愈加明显；随着H/W的增大，无量纲SIFs呈减小的趋势，这些都表明结构内部缺陷的增多会加速其破坏进程.相同模型尺寸下，边界裂纹的无量纲SIFs整体大于中心裂纹，说明了裂纹的约束条件对SIFs的影响较大.

参考文献（ References ）：

［1］刘钧玉, 林皋, 范书立, 等.裂纹面受荷载作用的应力强度因子的计算[J].计算力学学报, 2008, 25（5）: 621-626.（LIU Junyu, LIN Gao, FAN Shuli, et al.The calculation of stress intensity factor including the effects of surface tractions[J].Chinese Journal of Computational Mechanics, 2008, 25（5）: 621-626.（in Chinese））

［2］LIU J Y, LIN G, LI X C, et al.Evaluation of stress intensity factors for multiple cracked circular disks under crack surface tractions with SBFEM[J].China Ocean Engineering, 2013, 27（3）: 417-426.

［3］ZHONG H, LI C L, LI H J, et al.Stress intensity factors of interfacial crack with arbitrary crack tractions[J].IOP Conference Series:Earth and Environmental Science, 2019, 304（5）: 052111.

［4］陈白斌, 李建波, 林皋.无需裂尖增强函数的扩展比例边界有限元法[J].水利学报, 2015, 46（4）: 489-496, 504.（CHEN Baibin, LI Jianbo, LIN Gao.An extended scaled boundary finite element method without asymptotic enrichment of the crack tip[J].Journal of Hydraulic Engineering, 2015, 46（4）: 489-496, 504.（in Chinese））

［5］李亚, 易志坚, 王敏, 等.裂纹面局部均布荷载下Ⅰ型裂纹有限宽板应力强度因子[J].应用数学和力学, 2020, 41（10）:1083-1091.（LI Ya, YI Zhijian, WANG Min, et al.The stress intensity factor of a finite-width plate with a mode-Ⅰcenter crack subjected to uniform stress on the crack surface near the crack tip[J].Applied Mathematics and Mechanics, 2020, 41（10）: 1083-1091.（in Chinese））

［6］FETT T, RIZZI G.Weight functions for stress intensity factors and T-stress for oblique cracks in a halfspace[J].International Journal of Fracture, 2005, 132（1）: L9-L16.

［7］LI J, WANG X, TAN C L.Weight functions for the determination of stress intensity factor and T-stress for edge-cracked plates with built-in ends[J].International Journal of Pressure Vessels and Piping, 2004, 81（3）: 285-296.

［8］WALTERS M C, PAULINO G H, DODDS JR R H.Interaction integral procedures for 3-D curved cracks including surface tractions[J].Engineering Fracture Mechanics, 2005, 72（11）: 1635-1663.

［9］MUTHU N, MAITI S K, FALZON B G, et al.A comparison of stress intensity factors obtained through crack closure integral and other approaches using extended element-free Galerkin method[J].Computational Mechanics, 2013, 52（3）: 587-605.

［10］贾金生, 汪洋, 冯炜, 等.重力坝高压水劈裂模拟方法与特高重力坝设计准则初步探讨[J].水利学报, 2013, 44（2）:127-133.（JIA Jinsheng, WANG Yang, FENG Wei, et al.Simulation method of hydraulic fracturing and discussions on design criteria for super high gravity dams[J].Journal of Hydraulic Engineering, 2013, 44（2）: 127-133.（in Chinese））

［11］唐世斌, 刘向君, 罗江, 等.水压诱发裂缝拉伸与剪切破裂的理论模型研究[J].岩石力学与工程学报, 2017, 36（9）:2124-2135.（TANG Shibin, LIU Xiangjun, LUO Jiang, et al.Theoretical model for tensile and shear crack initiation at the crack tip in rock subjected to hydraulic pressure[J].Chinese Journal of Rock Mechanics and Engineering, 2017, 36（9）: 2124-2135.（in Chinese））

［12］杨绿峰, 徐华, 李冉, 等.广义参数有限元法计算应力强度因子[J].工程力学, 2009, 26（3）: 48-54.（YANG Lüfeng,XU Hua, LI Ran, et al.The finite element with generalized coefficients for stress intensity factor[J].Engineering Mechanics, 2009, 26（3）: 48-54.（in Chinese））

［13］徐华, 邓鹏, 蓝淞耀, 等.曲线裂纹裂尖SIFs等效分析的广义参数Williams单元确定方法[J].工程力学, 2020, 37（6）:34-41.（XU Hua, DENG Peng, LAN Songyao, et al.The determination method of Williams element with generalized degrees of freedom for equivalent analysis of SIFs at the curved crack tip[J].Engineering Mechanics, 2020,37（6）: 34-41.（in Chinese））

［14］中国航空研究院.应力强度因子手册[M].增订版.北京: 科学出版社, 1993.（Chinese Aeronautical Establishment.Handbook of Stress Intensity Factors[M].Revised ed.Beijing: Science Press, 1993.（in Chinese））