引　　言

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

裂纹面受荷时，其裂尖附近应力场和位移场较复杂，难以定量描述，通常仍以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单元需满足奇异区内裂纹面自由的边界条件，即${\sigma _\varphi } = 0$，${\tau _{\rho \varphi }} = 0\left( {\varphi = \pm {\text{π }}} \right)$，而裂纹面受荷时不满足该边界条件，仍需改进。

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

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

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

图  1  裂尖局部面荷载等效

注　为了解释图中的颜色，读者可以参考本文的电子网页版本，后同。

Figure  1.  Crack tip local surface load equivalence

下载: 全尺寸图片 幻灯片

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

式中，$ {K_{\rm I}} $，$ {K'_{\rm I}} $分别表示仅裂尖局部受分布压力和裂纹面受一对集中压力的Ⅰ型SIFs；p，p′分别表示裂尖局部加载和裂纹面作用集中力时的荷载系数，详见文献[14]。

假定$ {K_{\rm I}} = {K'_{\rm I}} $，则集中力Q与裂尖局部荷载$ {\sigma _0} $有如下关系：

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

式中，P表示等效荷载系数，$P = \dfrac{p}{{p'}}{\text{π}}$。

综上所述，只需根据荷载类型即可得出$\displaystyle \int_{a - c}^a {\sigma \left( X \right){\text{d}}X}$与$ {\sigma _0}a $的关系式，并根据《应力强度因子手册》^[14]查询荷载系数p与p′，确定等效荷载系数P，进而获得等效集中力Q。

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

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

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

图  2  裂尖奇异区网格划分及分布荷载等效

Figure  2.  Discretization and distributed load equivalence in the crack tip singular region

下载: 全尺寸图片 幻灯片

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

式中，u，v分别表示裂尖局部直角坐标系下x，y轴对应的位移分量；$\rho $，$\varphi $分别表示裂尖局部极坐标系下的极径和极角分量，$ {f_{l,11}}\left( \varphi \right) $，$ {f_{l,12}}\left( \varphi \right) $，$ {f_{l,21}}\left( \varphi \right) $，$ {f_{l,22}}\left( \varphi \right) $为三角函数，具体取值可参考文献[12]。

2.2   奇异区刚度方程集成

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

式中，${\bar {\boldsymbol{K}}^{\left( k \right)}}$为第k层子单元刚度矩阵；${{\boldsymbol{w}}^{\left( k \right)}}$为第k层子单元节点位移列阵，平面问题每个节点具有2个自由度，即${{\boldsymbol{w}}^{\left( k \right)}} = {[ {\begin{array}{*{20}{c}} {u_1^{\left( k \right)}}&{v_1^{\left( k \right)}}& \cdots &{u_8^{\left( k \right)}}&{v_8^{\left( k \right)}} \end{array}} ]^{\text{T}}}$；${\bar {\boldsymbol{p}}^{\left( k \right)}}$为第k层子单元节点荷载列阵，其中奇异区裂纹面的荷载需等效至边界，即奇异区内所有虚节点荷载为0，故当k=2~n时，${\bar {\boldsymbol{p}}^{\left( k \right)}} = {\boldsymbol{0}} {\text{.}}$

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

式中，${{\boldsymbol{T}}^{\left( k \right)}}$表示第k层子单元转换矩阵；${\boldsymbol{\phi}}$表示广义参数列阵。

将式(5)等号两边同时左乘${{\boldsymbol{T}}^{\left( k \right){\text{T}}}}$得

并可改写为

式中，${{\boldsymbol{K}}^{\left( k \right)}}$，${{\boldsymbol{p}}^{\left( k \right)}}$分别为第k层子单元的广义刚度矩阵和广义荷载列阵，且当k=2 ~ n时，${{\boldsymbol{p}}^{\left( k \right)}} = {\mathbf{0}} {\text{.}}$

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

即${{\boldsymbol{K}}_{{\text{ss}}}}{\boldsymbol{\phi}} {\text{ = }}{\mathbf{0}}$，${{\boldsymbol{K}}_{{\text{ss}}}}$表示第2~n层子单元广义刚度矩阵之和。

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

式中，${\boldsymbol{w}}_{\text{b}}^{\left( 1 \right)}$表示过渡单元的实节点位移列阵；${\boldsymbol{w}}_{\text{s}}^{\left( 1 \right)}$表示过渡单元的虚节点位移列阵；$\boldsymbol{T}_{10 \times 2(m + 1)}^{\prime(1)}$表示过渡单元的转换矩阵$\boldsymbol{T}_{16 \times 2(m + 1)}^{(1)}$的后10行，下角标$10 \times 2\left( {m + 1} \right)$表示该矩阵大小。

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

式中，$\boldsymbol{K}_{\mathrm{bs}}^{(1)}=\overline{\boldsymbol{K}}_{\mathrm{bs}}^{(1)} \boldsymbol{T}_{10 \times 2(m + 1)}^{\prime(1)}, \; \boldsymbol{K}_{\mathrm{sb}}^{(1)}=\boldsymbol{T}_{10 \times 2(m + 1)}^{\prime(1)} \overline{\boldsymbol{K}}_{\mathrm{sb}}^{(1)}, \; \boldsymbol{K}_{\mathrm{ss}}^{(1)}=\boldsymbol{T}_{10 \times 2(m + 1)}^{\prime(1) \mathrm{T}} \overline{\boldsymbol{K}}_{\mathrm{ss}}^{(1)} \boldsymbol{T}_{10 \times 2(m + 1)}^{\prime(1)}$。由于过渡单元实节点无外加荷载，故其荷载列阵为0。

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

式中，${\hat {\boldsymbol{K}}_{{\text{bb}}}}$表示奇异区边界实节点的刚度矩阵；${\hat {\boldsymbol{K}}_{{\text{ss}}}}$表示奇异区内部虚节点的广义刚度矩阵；${\hat {\boldsymbol{K}}_{{\text{bs}}}}$，${\hat {\boldsymbol{K}}_{{\text{sb}}}}$为奇异区边界实节点与内部虚节点的关联刚度矩阵；${{\boldsymbol{w}}_{\text{b}}}$表示奇异区边界实节点的位移列阵；${\bar {\boldsymbol{p}}_Q}$表示奇异区边界实节点等效荷载Q的荷载列阵。

2.3   整体刚度方程的集成

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

式中，${\bar {\boldsymbol{K}}_{{\text{rr}}}}$表示常规区节点的刚度矩阵；${\bar {\boldsymbol{K}}_{{\text{bb}}}}$表示常规区与奇异区交界节点的刚度矩阵；${\bar {\boldsymbol{K}}_{{\text{rb}}}}$，${\bar {\boldsymbol{K}}_{{\text{br}}}}$表示常规区与奇异区交界节点的关联刚度矩阵；${{\boldsymbol{w}}_{\text{r}}}$，${\bar {\boldsymbol{p}}_{\text{r}}}$分别表示常规区节点的位移列阵与荷载列阵。

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

2.4   SIFs求解

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

3.   算 例 分 析

3.1   裂纹面受均布压力

例1　如图1所示的边界裂纹矩形板，其板高为2H，宽为W，且有$W = 2H = 100{\text{ cm}}$，裂纹长度为a，板厚t取单位厚度，奇异区范围取边长为2c的正方形。裂纹面施加均布压力（$\lambda = 0$），${\sigma _0} = 1{\text{ kN/c}}{{\text{m}}^2}$，弹性模量$E = 3 \times {10^4}{\text{ MPa}}$，Poisson比$\nu = 0.3$。因篇幅所限，本文先研究此方法在Ⅰ型裂纹中的适用性，Ⅱ型裂纹另做研究。

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

图  3  常规区网格离散

Figure  3.  Discretization in the regular region

下载: 全尺寸图片 幻灯片

根据《应力强度因子手册》（P303、P304）的裂尖局部面荷载作用下荷载系数p，p′随a/W变化曲线图^[14]，通过GetData软件提取荷载系数p，p′代表性点的坐标，利用Excel软件拟合即可获得P与a/W的关系式：$P = 1.94 + 0.57\left( {{a \mathord{\left/ {\vphantom {a W}} \right. } W}} \right) - 1.51{\left( {{a \mathord{\left/ {\vphantom {a W}} \right. } W}} \right)^2}$。本模型参数的具体取值如表1所示。

表  1  模型参数（λ=0）

Table  1.  Model parameters (λ=0)

a/W	$\displaystyle \int_{a - c}^a {\sigma \left( X \right){\text{d} }X}$	P	Q
0.05	0.5	1.96	0.982
0.10	0.5	1.98	0.991
0.20	0.5	1.99	0.997
0.30	0.5	1.98	0.988
0.40	0.5	1.93	0.963
0.50	0.5	1.85	0.924
0.60	0.5	1.74	0.869
0.70	0.5	1.60	0.800

下载: 导出CSV 

| 显示表格

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

W单元的3个重要参数的建议值^[12]为$\alpha = 0.9$，$m = 20$，$n = 300$。将W单元计算结果与文献[14]的解对比，为直观体现各计算结果的区别，将SIFs无量纲化$( {{{{K_{\rm I}}} \mathord{/ {\vphantom {{{K_{\rm I}}} {( {{\sigma _0}\sqrt {{\text{π }}a} } )}}} } {( {{\sigma _0}\sqrt {{\text{π}}a\;} })}}})$，如图4所示。

图  4  边界裂纹面受均布压力计算结果

Figure  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 = 100{\text{ cm}}$，矩形板厚度t为单位厚度，裂纹长度为2a，奇异区边长2c，材料参数同本算例的边界裂纹板。

图  5  中心裂纹面受均布荷载

Figure  5.  The uniformly distributed load on the central crack surface

下载: 全尺寸图片 幻灯片

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

图  6  常规区网格离散

Figure  6.  Discretization in the regular region

下载: 全尺寸图片 幻灯片

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

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

图  7  中心裂纹面受均布压力计算结果

Figure  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，仅荷载形式不一样，此时施加于裂纹面的分布压力为$\sigma \left( X \right) = {\sigma _0} {\left| {1 - {X \mathord{\left/ {\vphantom {X a}} \right. } a}} \right|^\lambda }$，${\sigma _0} = 1{\text{ kN/c}}{{\text{m}}^2}$，$\lambda = 1$。结合本文的荷载等效处理方法和W单元求解算例中的SIFs。

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

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

图  8  裂纹面受线性压力计算结果

Figure  8.  Calculation results for linear pressure on the crack surface

下载: 全尺寸图片 幻灯片

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

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

图  9  无量纲SIFs随尺寸变化的计算结果

Figure  9.  Calculation results for dimensionless SIFs varying with sizes

下载: 全尺寸图片 幻灯片

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

4.   结　　论

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

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

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

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