断裂力学是力学的一个重要分支,也是计算力学中的一大难点.边界元法由于只需在边界上离散和基本解的奇异性特性，相比有限元法，在解决断裂问题上具有天然的优势.但是，对于含裂纹的弹性体，传统的位移边界积分方程是不适定的^[1]，需要补充额外的方程，对此，最常见的方法是多域边界元法，在弹性体内沿裂纹面，将原问题分为多个不含裂纹的子区域.Cruse等^[2]最早使用了多域边界元法解决三维裂纹问题，Luchi等^[3]发展了这种方法，开发了奇异单元提高了三维应力强度因子的计算精度. Gao等^[4]用多域边界元法求解了变系数问题的断裂力学问题.后来，Hong等^[5]提出了双边界积分方程法，在非裂纹边界和裂纹面的一面采用位移边界积分方程，在裂纹面的另一面采用面力积分方程，这种方法不需要人为的将计算体划分成多个子区域，是一种计算精度高的方法.文献^{[6, 7]}发展了双边界积分方程法，并将其成功的用于边裂纹的疲劳扩展问题的研究.Xie等^[8]提出了一种直接面力边界积分方程法，在非裂纹边界和裂纹面都采用面力边界积分方程，实现了三维有限域和无限域裂纹问题的求解.

由于双边界积分方程法引入了面力积分方程，因此包含于其中的强奇异和超强奇异边界积分的精确计算是其解决裂纹问题的关键.为此，许多学者提出了各种方法试图消除奇异性求解奇异积分.Pan等^[9]采用库特(Kutt)积分准则解决面力边界积分方程中的超奇异积分；陈梦成等^[10]利用完备函数系，将超奇异积分展开为截断基数，然后做极限运算获取近似的解析解，处理了断裂问题中的超奇异积分；文献^[11]中提出了一种在等参坐标系下近似展开的方法来消除积分奇异性，Mi等^[12]将其成功地用于断裂力学问题中的超奇异分解计算；高效伟^{[13, 14, 15]}将非奇异部分展开成关于等参坐标下距离的幂级数，提出了一种能处理二维、三维高阶奇异曲线、曲面积分的方法；胡金秀和高效伟等^{[16, 17]}将被积函数的非奇异部分展开成B样条基函数，提高了计算结果的稳定性和收敛性.高效伟和冯伟哲等^{[18, 19, 20, 21, 22]}提出了一种在全局尺寸下的投影线、面幂级数展开法，可以计算任意高阶奇异曲线、曲面积分.

本文采用一种改进的双边界积分方程法解决裂纹问题，在非裂纹边界采用位移边界积分方程，在裂纹面中的一面采用面力边界积分方程，以非裂纹边界的位移和面力以及裂纹面的间断位移为未知量.这样，既比传统的双边界积分方程的计算量小，所得到的裂纹间断位移又可直接用于应力强度因子的计算.本文采用投影线、面幂级数展开法直接处理面力边界积分方程中的超奇异积分，并通过算例3验证了此方法对于C$^{0}$连续性的单元超奇异积分计算的有效性. 此外，本文还构造了两种裂纹尖端单元的间断位移插值函数，用于模拟裂纹尖端附近的位移场.

如图1所示，考虑一带裂纹的有限体，$S$为非裂纹边界，$\Gamma^+$为裂纹上边界，$\Gamma^{-}$为裂纹下边界，对于各向同性均质问题，不考虑体力时，其位移边界积分方程为 $$ c_{ij} u_j (x^P) = \int_{S + \Gamma ^ + + \Gamma ^ - } U_{ij} (x,x^P)t_j (x) d\Gamma (x)-\\ \int_{S + \Gamma ^ + + \Gamma ^ - } T_{ij} (x,x^P)u_j (x)d\Gamma (x) =\\ \int_S \left( {U_{ij} (x,x^P)t_j (x) - T_{ij} (x,x^P)u_j (x) } \right) d\Gamma (x) -\\ \int_{\Gamma ^ + + \Gamma ^ - } T_{ij} (x,x^P)u_j (x)d\Gamma (x) \tag{1}$$ 其中$U_{ij} (x,x^P)$和$T_{ij} (x,x^P)$为弹性力学问题的基本解^[23]，$x^P$表示源点坐标，$x$表示场点坐标，$c_{ij} $是依赖于源点$P$处几何形状的系数，对于光滑边界取值为0.5. 由基本解的特性和裂纹面面力平衡条件可得 $$ \left. \begin{array}{l} {U_{ij}}({x^ + },{x^P}) = {U_{ij}}({x^ - },{x^P})\\ {T_{ij}}({x^ + },{x^P}) = - {T_{ij}}({x^ - },{x^P})\\ {t_j}({x^ + }) = - {t_j}({x^ - }) \end{array} \right\} \tag{2}$$ 其中$x^ + \in \Gamma ^ + $,$x^ - \in \Gamma ^ - $. 将式(2)代入式(1)整理可得 $$ c_{ij} u_j (x^P) = \int_S U_{ij} (x,x^P)t_j (x)d\Gamma (x) -\\ \int_S T_{ij} (x,x^P)u_j (x)d\Gamma (x) -\\ \int_{\Gamma ^ + } T_{ij} (x,x^P)\Delta u_j (x)d\Gamma (x)\tag{3}$$ 其中，$\Delta u_j (x) = u_j (x^ + ) - u_j (x^ - )$，表示裂纹面的张开和滑移位移.

图1(Fig.1)

图1 带裂纹的有限域模型 Fig.1 Finite field model with cracks

由位移边界积分方程经过一定的推导可得到源点位移光滑边界下应力边界积分方程 $$ \dfrac{1}{2}\sigma _{ij} (x^P) = \int_\Gamma U_{ijk} (x,x^P)t_k (x)d\Gamma (x) -\\ \int_\Gamma T_{ijk} (x,x^P)u_k (x)d\Gamma (x) \tag{4}$$ 式中

$$ U_{ijk} = \dfrac{(1 - 2\nu )(\delta _{ki} r_{,j} + \delta _{kj} r_{,i} - \delta _{ij} r_{,k} ) + \beta r_{,i} r_{,j} r_{,k} }{4\pi \alpha (1 - \nu )r^\alpha } \tag{5}$$

$$ T_{ijk} = \dfrac{1}{2\pi \alpha (1 - \nu )r^\beta }\Big\{ (1 - 2\nu )(\beta n_k r_{,i} r_{,j} + \\ n_j \delta _{ik} + n_i \delta _{jk} ) + \beta \nu (n_i r_{,j} r_{,k} + n_j r_{,i} r_{,k} ) - \\ (1 - 4\nu )n_k \delta _{ij} + \beta r_{,m} n_m [(1 - 2\nu )\delta _{ij} r_{,k} + \\ \nu (\delta_{ik} r_{,j} + \delta _{jk} r_{,i} ) - (\beta + 2) r_{,i} r_{,j} r_{,k}] \Big\} \tag{6}$$ 其中$\left\{ \begin{array}{l} \alpha = 1,\beta = 2,{\rm{二维问题}}\\ \alpha = 2,\beta = 3,{\rm{三维问题}} \end{array} \right.$

将式(4)乘以源点处外法线方向余弦$n_j (x^P)$可得面力边界积分方程 $$ \dfrac{1}{2}t_j (x^P) = n_i (x^P)\int_\Gamma U_{ijk} (x,x^P)t_k (x)d\Gamma (x) -\\ n_i (x^P)\int_\Gamma T_{ijk} (x,x^P)u_k (x)d\Gamma (x) \tag{7}$$ 如图1所示，将面力积分方程配置在裂纹面$\Gamma ^ + $上 $$ \dfrac{1}{2}\sigma _{ij} (x_{\Gamma ^ + }^P ) + \dfrac{1}{2}\sigma _{ij} (x_{\Gamma ^ - }^P ) = \\ \int_{S + \Gamma ^ + + \Gamma ^ - } U_{ijk} (x,x^P)t_k (x)d\Gamma (x) -\\ \int_{S + \Gamma ^ + + \Gamma ^ - } T_{ijk} (x,x^P)u_k (x)d\Gamma (x) \tag{8}$$ $$ \dfrac{1}{2}n_j (x_{\Gamma ^ + }^P )\left( {\sigma _{ij} (x_{\Gamma ^ + }^P )+\dfrac 12 \sigma _{ij} (x_{\Gamma ^ - }^P )} \right) = \\ n_j (x_{\Gamma ^ + }^P )\int_{S + \Gamma ^ + + \Gamma ^ - } {U_{ijk} (x,x^P)t_k (x)d\Gamma (x)} -\\ n_j (x_{\Gamma ^ + }^P )\int_{S + \Gamma ^ + + \Gamma ^ - } {T_{ijk} (x,x^P)u_k (x)d\Gamma (x)} \tag{9}$$ 其中$x_{\Gamma ^ + }^P \in \Gamma ^ + ,x_{\Gamma ^ - }^P \in \Gamma ^ - $，由于$n_j (x_{\Gamma ^ + }^P ) = - n_j (x_{\Gamma ^ - }^P )$,式(9)可写为 $$ \dfrac{1}{2}\left( {t_j (x_{\Gamma ^ + }^P ) - t_j (x_{\Gamma ^ - }^P )} \right) = \\ n_j (x_{\Gamma ^ + }^P )\int_{S + \Gamma ^ + + \Gamma ^ - } {U_{ijk} (x,x^P)t_k (x)d\Gamma (x)} -\\ n_j (x_{\Gamma ^ + }^P )\int_{S + \Gamma ^ + + \Gamma ^ - } {T_{ijk} (x,x^P)u_k (x)d\Gamma (x)} \tag{10}$$ 再由基本解的特性和裂纹面上面力平衡条件 $$ \left. \begin{align} & {{U}_{ijk}}({{x}^{+}},{{x}^{P}})={{U}_{ijk}}({{x}^{-}},{{x}^{P}}) \\ & {{T}_{ijk}}({{x}^{+}},{{x}^{P}})=-{{T}_{ijk}}({{x}^{-}},{{x}^{P}}) \\ & {{t}_{j}}({{x}^{+}})=-{{t}_{j}}({{x}^{-}}) \\ \end{align} \right\} \tag{11}$$ 将式(11)代人式(10)，整理可得 $$ \begin{array}{l} {t_j}(x_{{\Gamma ^ + }}^P) = {n_j}(x_{{\Gamma ^ + }}^P)\int_S {{U_{ijk}}(x,{x^P}){t_k}(x)d\Gamma (x)} - \\ {n_j}(x_{{\Gamma ^ + }}^P)\int_S {{T_{ijk}}(x,{x^P}){u_k}(x)d\Gamma (x)} - \\ {n_j}(x_{{\Gamma ^ + }}^P)\int_{{\Gamma ^ + }} {{T_{ijk}}(x,{x^P})\Delta {u_k}(x)d\Gamma (x)} \end{array} \tag{12}$$ 式(12)就是源点配置在裂纹面$\Gamma ^+ $上的面力边界积分方程，再与源点配置在非裂纹面$S$上的位移边界积分方程(3)组成双边界积分方程用于求解断裂问题. 当由其求出未知量后，如果还需要裂纹面$\Gamma ^+$上的位移，直接用位移边界积分方程(3)对源点在$\Gamma ^+ $上的节点配点即可. 对于无限域问题，则不需要位移边界积分方程，只需要源点配置在裂纹面$\Gamma ^+$上的面力边界积分方程，且式(12)退化为 $$ t_j (x_{\Gamma ^ + }^P ) = - n_j (x_{\Gamma ^ + }^P )\int_{\Gamma ^ + } {T_{ijk} (x,x^P)\Delta u{ }_k(x)d\Gamma } (x) \tag{13}$$

当源点在非裂纹面时本文采用位移边界积分方程，在裂纹面上时本文采用面力边界积分方程. 对于非裂纹面的离散与单元形函数的选择，与一般的边界元方法一样，不需特殊处理. 由于面力边界积分方程(12)中的$\int_{\Gamma^+} {T_{ijk} (x,x^P)\Delta u_k(x)d\Gamma } $积分核具有超强奇异性，其奇异积分存在的充分条件为源点处边界位移导数连续(C$^{1}$连续性)^[1]，而普通的拉格朗日形函数不能保证单元交点处的C$^{1}$连续性，所以针对裂纹面需要根据奇异积分的求解方法采用相应的离散方式和单元形函数.此外，由于裂纹尖端附近的位移场具有$\sqrt r $特征，应力场具有$\dfrac{1}{\sqrt r}$的奇异的奇异性，普通的拉格朗日形函数难以精确模拟裂尖附近的弹性场，所以对于裂纹尖端单元的形函数也需要特殊处理.

对于非裂纹面，通常是离散成采用拉格朗日形函数的线性单元或二次单元，本文也不例外，离散成二次单元. 对于裂纹面的离散，最普遍的做法是采用非连续元，即源点都取在单元的内部，因为拉格朗日形函数能保证单元内部位移导数的连续性，也就能满足哈达玛主值积分源点处的光滑性要求.但是，非连续元的使用不仅不能保证单元交点处位移导数的连续性，甚至不能保证单元交点处位移的连续性，这将造成计算精度的损失^[1].如果将与单元交点相连的所有单元结合在一起考虑，从物理上讲各单元积分的发散部分应互相抵消，如此便可不采用非连续元，本文使用的奇异积分方法正是考虑了这点，所以无需采用非连续元便可得到相对精确的解，这点在算例3中给出了验证.

对于裂纹尖端单元，由于其周围有一边无相邻单元，无法抵奇异积分发散项，所以对裂纹尖端单元，我们采用非连续元. 而对于裂纹面上非裂尖单元，我们采用连续元. 离散后的位移边界积分方程和面力边界积分方程分别如下 $$ c_{ij} t_j (x_{\Gamma ^ + }^P ) = \sum_{e = 1}^{NE1} {\left\{ {\sum_{\alpha = 1}^{M1} {t_j^\alpha (x^Q)} \int_{S_e } {U_{ij} [x,x^P(\xi ,\eta )]N_\alpha (\xi ,\eta )Jd\xi d\eta } } \right\}}-$$ $$ \sum_{e = 1}^{NE1} {\left\{ {\sum_{\alpha = 1}^{M1} {u_j^\alpha (x^Q)} \int_{S_e } {T_{ij} [x,x^P(\xi ,\eta )]N_\alpha (\xi ,\eta )Jd\xi d\eta } } \right\}} +$$ $$ \sum_{e = 1}^{NE2} {\left\{ {\sum_{\alpha = 1}^{M2} {\Delta u_j^\alpha (x^Q)} \int_{\Gamma _e^ + } {T_{ij} [x,x^P(\xi ,\eta )]N_\alpha (\xi ,\eta )Jd\xi d\eta } } \right\}} \tag{14}$$ $$ t_j (x_{\Gamma ^ + }^P ) = n_i (x_{\Gamma ^ + }^P )\sum_{e = 1}^{NE1} {\left\{ {\sum_{\alpha = 1}^{M1} {t_j^\alpha (x^Q)} \int_{S_e } {U_{ijk} [x,x^P(\xi ,\eta )]N_\alpha (\xi ,\eta )J d\xi d\eta } } \right\}} -$$ $$ n_i (x_{\Gamma ^ + }^P )\sum_{e = 1}^{NE1} {\left\{ {\sum_{\alpha = 1}^{M1} {u_j^\alpha (x^Q)} \int_{S_e } {T_{ijk} [x,x^P(\xi ,\eta )]N_\alpha (\xi ,\eta )Jd\xi d\eta } } \right\}} +$$ $$ n_i (x_{\Gamma ^ + }^P )\sum_{e = 1}^{NE2} {\left\{ {\sum_{\alpha = 1}^{M2} {\Delta u_j^\alpha (x^Q)} \int_{\Gamma _e^ + } {T_{ijk} [x,x^P(\xi ,\eta )]N_\alpha (\xi ,\eta )Jd\xi d\eta } } \right\}} \tag{15}$$

对于非裂纹面单元和裂纹面非裂尖单元，本文采用拉格朗日二次形函数.对于裂尖单元，在求解裂纹问题的子区域法中通常采用 1/4节点奇异单元；本文采用双积分方程解法，裂纹面的未知量只有裂纹间断位移，没有面力，只需要裂纹间断位移的插值函数.对于裂尖单元，本文提供了3种形式的间断位移插值函数，其中裂尖单元1采用现有的形函数，裂尖单元2和裂尖单元3采用本文构造的插值函数.

裂尖单元1： 裂纹间断位移插值函数与单元几何形函数相同，二维问题采用3节点二次单元形函数，三维问题采用8节点或9节点二次单元形函数.

裂尖单元2： 设裂纹间断位移为 $$ \Delta u_j = \sqrt r \Delta \overline {u_j } = \sqrt r \sum_{\alpha = 1}^m {N_\alpha } \Delta \overline {u_j^\alpha } = \sum_{\alpha = 1}^m {{N}'_\alpha } \Delta \overline {u_j^\alpha } $$ 其中，$r$为源点到裂尖的距离，$\Delta \overline {u_j^\alpha } $表示消除$\sqrt r $特性后的裂纹间断位移.$$ {N'_\alpha } = \left\{ \begin{array}{l} {N_\alpha }(\xi ,\eta )\sqrt {r(\xi )} ,{\rm{二维问题}}\\ {N_\alpha }(\xi ,\eta )\sqrt {r(\xi ,\eta )} ,{\rm{三维问题}} \end{array} \right.$$ $r(\xi ) = l_\xi ^n (1 + \xi )$，$\xi = - 1$为裂尖, 三维单前沿单元($\eta = - 1$为前沿)$$ r(\xi ,\eta ) = l_\eta ^n (1 + \eta )$$ 三维双前沿单元($\xi = - 1$,$\eta = - 1$为前沿)$$ r(\xi ,\eta ) = l_\xi ^n l_\eta ^n (1 + \xi )(1 + \eta )$$ 其中$N_\alpha (\xi )$和$N_\alpha (\xi ,\eta )$为单元几何形函数，$l_\xi ^n $和$l_\eta ^n $分别表示裂尖单元$n$在$\xi $和$\eta $方向上的长度.

裂尖单元3： 根据裂纹尖端附近位移场特征^[24] $$ \Delta u = \lambda _1 \sqrt r + \lambda _2 r + \lambda _3 r^{\tfrac{3}{2}} + \cdots $$ 设张开位移插值函数为：

二维($\xi = - 1$为裂尖)$$ N_\alpha (\xi ) = a_{_1 }^\alpha \sqrt {1 + \xi } + a_2^\alpha (1 + \xi ) + a_3^\alpha (1 + \xi )^{\tfrac{3}{2}} \tag{16}$$ 三维单前沿单元($\eta = - 1$为前沿)$$ N_\alpha (\xi ,\eta ) = a_{_1 }^\alpha \sqrt {1 + \eta } + a_2^\alpha (1 + \eta ) + a_3^\alpha (1 + \eta )^{\tfrac{3}{2}} +\\ a_4^\alpha \xi \sqrt {1 + \eta } + a_5^\alpha \xi (1 + \eta ) + a_6^\alpha \xi (1 + \eta )^{\tfrac{3}{2}} +\\ a_7^\alpha \xi ^2\sqrt {1 + \eta } + a_8^\alpha \xi ^2(1 + \eta ) + a_9^\alpha \xi ^2(1 + \eta )^{\tfrac{3}{2}} \tag{17}$$ 三维双前沿单元($\xi = - 1,\eta = - 1$为前沿)$$ N_\alpha (\xi ,\eta ) = a_{_1 }^\alpha \sqrt {(1 + \xi )(1 + \eta )} + a_2^\alpha (1 + \xi )(1 + \eta ) +\\ a_3^\alpha [(1 + \xi )(1 + \eta )]^{\tfrac{3}{2}} + a_4^\alpha \xi \sqrt {(1 + \xi )(1 + \eta )} +\\ a_5^\alpha \xi (1 + \xi )(1 + \eta ) + a_6^\alpha \xi [(1 + \xi )(1 + \eta )]^{\tfrac{3}{2}} +\\ a_7^\alpha \eta \sqrt {(1 + \xi )(1 + \eta )} + a_8^\alpha \eta (1 + \xi )(1 + \eta ) +\\ a_9^\alpha \eta [(1 + \xi )(1 + \eta )]^{\tfrac{3}{2}} \tag{18}$$ 系数$a_\beta ^\alpha $的值，可根据插值函数的性质由程序自动求出.

在式(14)中，源点配置在非裂纹面$S$上，右端第1项积分存在强奇异性，采用刚体位移法间接求解；第2项存在弱奇异性，采用对数积分求解；第3项无积分奇异性，直接用高斯数值积分求解.

在式(15)中，源点配置在裂纹面$\Gamma ^ + $上，右端第1、第2项积分都不存在奇异性，直接采用高斯数值积分计算；右端第3项存在超强奇异性，无法采用刚体位移法进行间接计算，需采用高阶奇异积分的直接法进行计算.

对于面力边界积分方程中的超强奇异积分，将其写出如下形式 $$ \int_{\Gamma _e^ + } {T_{ijk} [{ x},{ x}^P(\xi ,\eta )]N_\alpha (\xi ,\eta )J d\xi d\eta } =\\ \int\limits_{\Gamma _e^ + } {\dfrac{f(x,x^P)}{r^\beta (x,x^P)}} d\Gamma (x) = I_e (x^P) \tag{19}$$ 式中，$f(x,x^P)$为核函数正则部分，$r^\beta (x,x^P)$代表奇异性程度，采用文献 ^[20] 的方法计算有限部分积分. 首先将原边界单元在等参坐标系原点 ($\xi = 0$和$\eta = 0$)处的切线或切平面进行投影，并在投影线或投影面上沿等参坐标建立局部坐标系如图2和图3.

图2(Fig.2)

图2 二维边界单元在投影线上的局部坐标系及投影线上变量的定义 Fig.2 Local coordinate on projection line of a two dimension element and variables defining on the projection line

图3(Fig.3)

图3 三维边界单元在投影平面上的局部坐标系 Fig.3 Local coordinate on projection plane of a three dimension boundary element

坐标间的转换关系为 $$ {x}'_i = L_{ij} (x_j - x_j^o ) \tag{20}$$ $$ x_i - x_i^o = L_{ji} {x}'_j \quad \tag{21}$$ 式中，$x_i^o $为新坐标系的原点在旧坐标系中的坐标，$L_{ij}$ 是新坐标系的坐标轴相对于旧坐标系坐标轴的方向余弦.

投影线、投影平面内的径向距离$\rho $及其对新坐标的偏导数$\rho _{,I} $(见图2和图4)，其定义为 $$\rho = \sqrt {\left( {{x}'_1 - {x}'^P_1 } \right)^2 + \left( {{x}'_2 - {x}'^P_2 } \right)^2} \tag{22}$$ $$\rho _{,I} = \dfrac{\partial \rho }{\partial {x}'_I } = \dfrac{{x}'_I - {x}'^P_I }{\rho } \tag{23}$$

图4(Fig.4)

图4 投影平面内各量的定义 Fig.4 Variables defining on the projection plane

通过一系列推导可得出原坐标系下各个量在投影线、面中局部坐标系下的表达式，详细推导过程见文献 ^[20]，最终式(19)可写为 $$ {I_e}({x^P}) = \left\{ \begin{array}{l} \mathop {\lim }\limits_{\varepsilon \to 0} \int_{{\rho _\varepsilon }}^{{\rho _E}} {\frac{{f({x^P},x)}}{{{g^\beta }(\rho ){\rho ^\beta }n_i^0{n_i}}}d\rho } = \\ \mathop {\lim }\limits_{\varepsilon \to 0} \int_{{\rho _\varepsilon }}^{{\rho _E}} {\frac{{\bar F(\rho )}}{{{\rho ^\beta }}}d\rho } ,{\rm{二维问题}}\\ \int_A {\frac{{f({x^P},x)}}{{{g^\beta }(\rho ){\rho ^\beta }n_i^0{n_i}}}} dA(x),{\rm{三维问题}} \end{array} \right. \tag{24}$$ 对于二维问题，$\bar {F}(\rho )$只是$\rho $的函数，将其通过幂级数展开 $$ \bar {F}(\rho ) = \sum_{n = 0}^N B^{(n)} \rho^n \tag{25}$$ 其中，$N$为幂级数阶数，$B^{(n)} $为待求系数. 在积分区间$(0,\rho )$布置$N + 1 $个点$\left( {0,\rho _1 ,\rho _2 ,\cdots ,\rho _N } \right)$，求解式(25)便可确定$B^{(n)} $.对于三维问题，将上式在投影平面内的积分通过采用径向积分法^[25]转换成投影平面边界线$L$ (见图4)的线积分 $$ I({ x}^P) = \int_L {\dfrac{1}{\rho }} \dfrac{\partial \rho }{\partial n^L}F d L \tag{26}$$ 式中，$\partial \rho / \partial n^L = \rho _{,I} n_I^L $，$n_I^L $为边界线$L$的外法线方向余弦在$I$方向的分量(见图4)；$F$为径向积分 $$ F = \mathop {\lim }\limits_{\varepsilon \to 0} \int_{\rho _\varepsilon }^{\rho _L } {\dfrac{\bar {F}(\rho )}{\rho ^{\beta - 1}} d\rho } \tag{27}$$ 式中 $$ \bar {F}(\rho ) = \dfrac{f({ x}^P,{ x})}{g^\beta (\rho )n_i^0 n_i } \tag{28}$$ 此时，$\bar {F}(\rho )$与二维情况一样只是$\rho $的函数，将其通过幂级数展开，见式(25). 系数$ B $确定之后，就可求出二维积分项${I_\varepsilon }({x^P})$和三维径向积分项$F$ $$ I_e (x^P) = \mathop {\lim }\limits_{\varepsilon \to 0} \int_{\rho _\varepsilon }^{\rho _E } {\dfrac{\bar {F}(\rho )}{\rho ^\beta }d\rho } =$$ $$ \sum_{n = 0}^N {B^{(n)} \mathop {\lim }\limits_{\varepsilon \to 0} \int_{\rho _\varepsilon }^{\rho _E } {\rho ^{n - \beta } d\rho = } } \sum_{n = 0}^N B^{(n)} E_n \tag{29}$$ $$F = \mathop {\lim }\limits_{\varepsilon \to 0} \int_{\rho _\varepsilon }^{\rho _L } {\dfrac{\bar {F}(\rho )}{\rho ^\alpha }d\rho } = \sum_{n = 0}^N B^{(n)} \mathop {\lim }\limits_{\varepsilon \to 0} \int_{\rho _\varepsilon }^{\rho _L } \rho ^{n - \alpha }d\rho =$$ $$ \sum_{n = 0}^N B^{(n)} {E}'_n \tag{30}$$ 式中，$\alpha = \beta - 1$ $$ {E_n} = \left\{ {\begin{array}{*{20}{l}} {\frac{1}{{n - \alpha + 1}}\left( {\rho _E^{n + 1 - \beta } - \mathop {\lim }\limits_{\varepsilon \to 0} \rho _\varepsilon ^{n + 1 - \beta }} \right),n \ne \beta - 1}\\ {\ln {\rho _E} - \mathop {\lim }\limits_{\varepsilon \to 0} \ln {\rho _\varepsilon }, n = \beta - 1} \end{array}} \right. \tag{31}$$ 对于确定的物理问题，上述积分结果应该是一有限值，也就是说，当考虑源点周围所有单元的结果后，式(31)发散项的最终贡献应为零. 但是，要消除该发散项，要求源点周围单元都采用全局尺寸下的无穷小尺度$\varepsilon $，为此，需将投影线上局部尺度$\rho _\varepsilon $转化为全局尺度$\varepsilon $，经过转化并消除发散项后可得 $$ E_n = \left\{ \!\!\begin{array}{ll} \dfrac{1}{n - \beta + 1}\left( {\dfrac{1}{\rho _E ^{\beta - n - 1}} - H_{\beta - n - 1} } \right) , & 0 ≤ n ≤ \beta - 2 \\ \ln \rho _E - \ln H_0 , & n = \beta - 1 \\ \dfrac{1}{n - \beta + 1}\rho _E ^{n - \beta + 1} , & n > \beta - 1 \end{array} \!\! \right. \tag{32}$$ 采用同样的处理方法，可得 $$ {E}'_n = \left\{ \!\! \begin{array}{ll} \dfrac{1}{n - \alpha + 1}\left( {\dfrac{1}{\rho _L ^{\alpha - n - 1}} - H_{\alpha - n - 1} } \right) , & 0 ≤ n ≤ \alpha - 2 \\ \ln \rho _L - \ln H_0 , & n = \alpha - 1 \\ \dfrac{1}{n - \alpha + 1}\rho _L ^{n - \alpha + 1} , & n > \alpha - 1 \end{array} \!\! \right. \tag{33}$$ 式(32)和式(33)中的$H$的值详见文献 ^[20]，式(33)中的$\rho _L $是面单元边界线在切平面上的投影曲线$L$上的取值.当源点$P$位于轮廓线$L$的一边上，$\rho \to 0$时，有$\partial \rho / \partial n^L \to 0$，无奇异性，因此该边不需要考虑，故式(26)中的积分可以直接利用常规高斯积分进行计算.

由于本文采用的双边界积分方程是以裂纹面的间断位移为未知量，所以直接取裂纹尖端附近点的间断位移，采用两点位移公式计算裂尖点应力强度因子^[1]，既简单又方便，其计算公式如下 $$ {K}_{i} = H^i\left[{\dfrac{r_B }{(r_B - r_A )\sqrt {r_A } }\Delta u_j (r_A ) - \dfrac{r_A }{(r_B - r_A )\sqrt {r_B } }\Delta u_j (r_B )} \right] \tag{34}$$ 式中，$i$表示Ⅰ,Ⅱ,Ⅲ,$j$表示2,1,3,$A$和$B$分别表示靠近裂尖的两点，$r_A $和$r_B $为$A$点和$B$点到裂尖的距离.$$ H^{\rm Ⅰ} = H^{\rm Ⅱ} = \dfrac{G}{2(1 - \nu )}\sqrt {\dfrac{\pi }{2}} ,\ H^{\rm Ⅲ} = \dfrac{G}{2}\sqrt {\dfrac{\pi }{2}} \tag{35}$$ 式中，$G$为材料剪切模量，$\nu $为材料的泊松比，平面应力时$\nu = \dfrac{\nu }{1 + \nu }$.

算例1： 二维矩形中心斜裂纹

如图5所示矩形板，矩形长$2h$，宽为$2w$，且$h / w = 2$，矩形中心有一斜裂纹，裂纹长为$2a$，倾斜角为$\theta $，两端受均布拉力$t$.采用上文所述的边界积分方程法，在非裂纹面布置32个二次连续元，在裂纹面其中一面上布置6个二次连续元和2个二次非连续裂尖单元如图6所示，其中裂尖单元分别采用两种插值函数，NTCES=1表示采用上文3.2节中的裂尖单元1，NTCES=3表示采用上文3.2节中的裂尖单元3. 首先，固定$\theta = 45^\circ $，使$a / w$分别取0.1, 0.2到0.8，讨论规则化应力强度因子$K / (t\sqrt {\pi a} )$随$a / w$的变化规律；然后，固定$a / w = 0.2$，使$\theta $分别取15$^{\circ}$,30$^{\circ }$,45$^{\circ }$,60$^{\circ }$和75$^{\circ }$，讨论规则化应力强度因子$K / (t\sqrt {\pi a} )$随$\theta $的变化规律.我们将所得结果与应力强度因子手册上^[26]的参考值作比较，分别如图7和图8所示. 从图中可看出：$K_{\rm Ⅰ}$和K_Ⅱ都随着$a / w$的增大而增大，且K_Ⅰ的增速要大于K_Ⅱ的增速；$K_{\rm Ⅰ}$随着$\theta $的增大而减小，K_Ⅱ随着$\theta $的增大先增大后减小,在$\theta $等于45$^{\circ }$时最大. 本文方法所得结果与应力强度因子手册^[26]给出的参考值基本一致，最大误差不超过1%.

图5(Fig.5)

图5 矩形中心斜裂纹模型 Fig.5 Model of rectangle with inclined crack in center

图6(Fig.6)

图6 二维矩形中心斜裂纹网格模型 Fig.6 Mesh model of rectangle with inclined crack in center

图7(Fig.7)

图7 应力强度因子随变化规律($\theta = 45^\circ $) Fig.7 Change rule of SIF along with $a / w$ ($\theta = 45^\circ $)

图8(Fig.8)

图8 应力强度因子随$\theta $变化规律($a / w = 0.2$) Fig.8 Change rule of SIF along with $\theta $ ($a / w = 0.2$)

详细结果数据与误差见表1 $\sim$表4，从表中可以看出采用两种裂尖单元都能得到精度不错的结果.其中采用了特别构造的裂尖单元3的计算结果有略微的改善，在$a / w$小于0.3时尤为明显；随着$a / w$的增大，相对误差也增大. 如表5所示，在$a / w$为0.8时，通过加密裂纹上的单元，采用裂尖单元3所得结果精度有所提高，逐渐收敛于参考值；而采用普通的裂尖单元1所得结果收敛情况并不太理想.由此算例可以验证本文方法对二维深埋裂纹计算的有效性，即使不采用特殊构造的裂尖单元也能得到比文献 ^[27]更高精度的结果.

表1(Table 1)

表1 $a / w$变化时规则化应力强度因子KⅠ计算结果与误差Table 1 Computed result and error of normalized SIF KⅠ along with $a / w$ $\theta = 45^\circ$ $a / w$ KⅠ /$(t\sqrt {\pi a} )$ Relative error NTCES=1 NTCES=3 reference NTCES=1 NTCES=3 0.1 0.503 1 0.503 9 0.504 6 $-0.303%$ $-0.138%$ 0.2 0.516 7 0.517 4 0.518 1 $-0.271%$ $-0.145%$ 0.3 0.539 3 0.539 4 0.540 6 $-0.237%$ $-0.215%$ 0.4 0.570 8 0.570 2 0.571 9 $-0.189%$ $-0.293%$ 0.5 0.611 4 0.609 8 0.611 9 $-0.089%$ $-0.350%$ 0.6 0.661 6 0.658 6 0.661 1 0.082% $-0.381%$ 0.7 0.723 4 0.718 1 0.721 0 0.334% $-0.401%$ 0.8 0.800 2 0.791 3 0.795 0 0.658% $-0.459%$

表1 $a / w$变化时规则化应力强度因子KⅠ计算结果与误差 Table 1 Computed result and error of normalized SIF KⅠ along with $a / w$

表2(Table 2)

表2 $a / w$变化时规则化应力强度因子KⅡ计算结果与误差Table 2 Computed result and error of normalized SIF KⅡ along with $a / w$ $\theta = 45^\circ$ $a / w$ KⅡ /$(t\sqrt {\pi a} )$ Relative error NTCES=1 NTCES=3 reference NTCES=1 NTCES=3 0.1 0.500 3 0.501 1 0.501 8 $-0.303%$ $-0.130%$ 0.2 0.505 6 0.506 4 0.507 2 $-0.320%$ $-0.160%$ 0.3 0.514 5 0.515 0 0.516 2 $-0.332%$ $-0.230%$ 0.4 0.527 1 0.527 3 0.529 0 $-0.367%$ $-0.315%$ 0.5 0.543 6 0.543 6 0.545 8 $-0.405%$ $-0.398%$ 0.6 0.564 7 0.564 5 0.567 4 $-0.479%$ $-0.510%$ 0.7 0.591 5 0.591 1 0.595 0 $-0.596%$ $-0.657%$ 0.8 0.625 9 0.625 4 0.630 0 $-0.655%$ $-0.735%$

表2 $a / w$变化时规则化应力强度因子KⅡ计算结果与误差 Table 2 Computed result and error of normalized SIF KⅡ along with $a / w$

表3(Table 3)

表3 $\theta $变化时规则化应力强度因子KⅠ计算结果与误差Table 3 Computed result and error of normalized SIF KⅠ along with $\theta $ $a / w = 0.2$ $\theta $ KⅠ /$(t\sqrt {\pi a} )$ Relative error NTCES=1 NTCES=3 reference NTCES=1 NTCES=3 15$^{\circ}$ 0.954 8 0.956 2 0.957 7 $-0.301%$ $-0.160%$ 30$^{\circ}$ 0.770 7 0.771 8 0.773 0 $-0.295%$ $-0.160%$ 45$^{\circ}$ 0.516 7 0.517 4 0.518 1 $-0.271%$ $-0.145%$ 60$^{\circ}$ 0.259 8 0.260 1 0.260 5 $-0.282%$ $-0.164%$ 75$^{\circ}$ 0.069 9 0.070 0 0.070 1 $-0.316%$ $-0.204%$

表3 $\theta $变化时规则化应力强度因子KⅠ计算结果与误差 Table 3 Computed result and error of normalized SIF KⅠ along with $\theta $

表4(Table 4)

表4 $\theta $变化时规则化应力强度因子KⅡ计算结果与误差Table 4 Computed result and error of normalized SIF KⅡ along with $\theta $ $a / w = 0.2$ $\theta $ KⅡ /$(t\sqrt {\pi a} )$ Relative error NTCES=1 NTCES=3 reference NTCES=1 NTCES=3 15$^{\circ}$ 0.250 2 0.250 7 0.251 0 $-0.307%$ $-0.132%$ 30$^{\circ}$ 0.435 3 0.436 0 0.436 7 $-0.324%$ $-0.156%$ 45$^{\circ}$ 0.505 6 0.506 4 0.507 2 $-0.320%$ $-0.160%$ 60$^{\circ}$ 0.440 4 0.441 0 0.441 7 $-0.301%$ $-0.150%$ 75$^{\circ}$ 0.255 3 0.255 7 0.256 0 $-0.273%$ $-0.128%$

表4 $\theta $变化时规则化应力强度因子KⅡ计算结果与误差 Table 4 Computed result and error of normalized SIF KⅡ along with $\theta $

表5(Table 5)

表5 $a / w = 0.8$时计算结果误差随裂纹单元数变化情况Table 5 Error of normalized SIF $K$ along with number of crack element $a / w = 0.8$ Element's number KⅠ /$(t\sqrt {\pi a} )$ Relative error KⅡ /$(t\sqrt {\pi a} )$ Relative error NTCES=1 {NTCES=3} NTCES=1 NTCES=3 8 0.658% $-0.459%$ $-0.655%$ $-0.735%$ 16 0.115% $-0.287%$ $-0.745%$ $-0.626%$ 24 $-0.063%$ $-0.237%$ $-0.795%$ $-0.586%$ 32 $-0.134%$ $-0.190%$ $-0.816%$ $-0.553%$

表5 $a / w = 0.8$时计算结果误差随裂纹单元数变化情况 Table 5 Error of normalized SIF $K$ along with number of crack element

算例2： 三维无限域圆形裂纹面受均布剪力

如图9所示无限域圆盘裂纹，半径$r = 1$，受面内均布剪力$q = 1$，杨氏模量$E = 1$，泊松比$\nu = 0.3$.网格模型如图10所示，采用24个8节点裂尖非连续元和108个8节点连续元，对于裂尖单元分别采用两种张开位移插值函数，NTCES=1表示采用上文3.2节中的裂尖单元1，NTCES=2表示采用上文3.2节中的裂尖单元2.

图9(Fig.9)

图9 三维圆形裂纹受均布剪力 Fig.9 A penny shape crack under shear

图10(Fig.10)

图10 三维圆形裂纹网格模型 Fig.10 Mesh model of penny shape crack

根据应力强度因子手册^[28]，其精确解为 $$ K_{{\rm Ⅰ}A} = 0 ,K_{{\rm Ⅱ}A} = \dfrac{4q\cos w}{\pi (2 - \nu )}\sqrt {\pi a}$$ $$ K_{{\rm Ⅲ}A} = \dfrac{4(1 - \nu )q\sin w}{\pi (2 - \nu )}\sqrt {\pi a} $$

根据本文方法分别计算出采用两种插值函数下的应力强度因子，并与精确解作比较，最大误差都在1%以内. 从图11和图12中可以看出，采用两种裂尖张开位移插值函数都能得到精度较高的解，一方面验证了本文构造的裂尖单元2的有效性，另一方面也验证了本文方法对于普通的二次裂尖张开位移插值函数也同样适用. 在本例中，相对于普通的裂尖单元1，采用构造的裂尖单元2后，并没有带来精度的提高；裂尖单元2在原插值函数上乘上了$\sqrt r $，虽然考虑了裂尖附近位移场特征，但也会导致裂尖单元与非裂尖单元公共点的值不连续，这会给计算结果带来一定的误差.

图11(Fig.11)

图11 三维圆形裂纹规则化应力强度因子KⅠ,KⅡ Fig.11 Normalized SIF KⅠ,KⅡ of penny shape crack

图12(Fig.12)

图12 三维圆形裂纹规则化应力强度因子KⅢ Fig.12 Normalized SIF KⅢ of penny shape crack

算例3：三维无限域、有限域深埋方形裂纹

如图13所示，在有限域长方体中有一方形裂纹，长方体高$2H$，长和宽都为$W$，方形裂纹长为$2a$，其中 $H / W = 1$，$a / W = 0.25$，在长方体底部和顶部施加垂直于裂纹面的单位拉应力$\sigma $.在本例中，非裂纹面采用8节点连续单元，对无限域问题，则只需裂纹面网格，不需要非裂纹面网格.裂纹面分别采用两种网格计算，一种裂尖单元采用9节点非连续元，非裂尖单元采用9节点连续元，一种整个裂纹面都采用9节点非连续元，如图14所示.

图13(Fig.13)

图13 三维有限域方形裂纹模型 Fig.13 Model of square crack in finite field

图14(Fig.14)

图14 裂纹面网格 Fig.14 Mesh of crack plane

选取上文3.2节的裂尖单元1的形函数进行计算，结果如图15和图16所示，其中$x$表示裂尖点离裂纹前沿一端的距离.从图中可以看出，裂纹前沿中点的应力强度因子最大，在裂纹前沿的两端应力强度因子趋于零，本文所得结果变化曲线与潘尔年等所得结果基本一致.对于无限域方形裂纹，只采用裂尖非连续元和裂纹面全部采用非连续元计算所得结果基本相等，最大规则化应力强度因子分别为0.754 5和0.752 2，与肖洪天等^[29]的结果0.756 4、Pan等^[9]的结果0.762 6、Weaver^[30]的结果0.74基本一致.这充分说明本文方法与现存方法非常吻合，也验证本文所采用的超奇异有限部分积分方法由于考虑了源点周围所有单元的影响可以不采用非连续元的观点.对于有限域方形裂纹，只采用裂尖非连续元的计算结果比裂纹面全部采用非连续元的计算结果稍大，最大规则化应力强度因子分别为0.818 1和0.805 1，这与肖洪天等^[29]的结果0.811 2、文献^[31]的结果0.81、Pan等^[9]的结果0.818 3都非常吻合.此外，采用本文构造的裂尖单元2和裂尖单元3同样能得到一致的结果，如表6所示.

图15(Fig.15)

图15 无限域方形裂纹前沿规则化应力强度因子 Fig.15 Normalized SIF of square crack in infinite field

图16(Fig.16)

图16 有限域方形裂纹前沿规则化应力强度因子 Fig.16 Normalized SIF of square crack in finite field

表6(Table 6)

表6 采用不同裂尖单元的计算结果比较Table 6 Comparison of computed result for different crack tip elements Problem type NTCES KⅠ /$(t\sqrt {\pi a} )$ Maximun Tip discontinuous All discontinuous 1 0.754 5 0.752 2 infinite 2 0.753 5 0.751 3 3 0.754 0 0.752 6 1 0.818 1 0.805 1 finite 2 0.818 9 0.806 5 3 0.815 4 0.805 7

表6 采用不同裂尖单元的计算结果比较 Table 6 Comparison of computed result for different crack tip elements

算例4: 有限厚度板半圆形表面裂纹

如图17所示有限厚度板半圆型裂纹，$a = b = 1$，$h = w = t = 5 $. 板的材料参数：泊松比为0.3，剪切模量为10. 板上下表面都受单位拉应力. 我们将本文方法所得规则化应力强度因子K_Ⅰ/$(\sigma \sqrt {\pi a} )$结果与文献^[32]采用有限元方法所得结果作比较，如图18所示. 从图中可以看出，$2\phi /\pi $小于0.1时，本文结果与文献^[32]的结果基本一致，在$ 2\phi /\pi $接近0时，本文结果与文献^[32]的结果差异较大，但最大相对差异在2%以内，详细差异见表7.由于文献^[32]的结果是基于有限元的计算结果，并非精确的理论解，至于两种计算结果的精度高低，尚无法作出判断.

图17(Fig.17)

图17 有限厚度板半圆型裂纹几何模型 Fig.17 A semi-elliptical surface crack in a finite-thickness plate

图18(Fig.18)

图18 有限厚度板半圆形裂纹规则化应力强度因子 Fig.18 Model-Ⅰ SIF of semi-elliptical surface crack in a finite-thickness plate

表7(Table 7)

表7 本文方法与文献[32]方法的比较Table 7 Comparisons between our method and Ref.[32] $ 2 \phi/\pi $ Ref.[32] Our method Relative difference 0 0.747 0.759 7 1.705% 0.125 0.729 0.7397 1.465% 0.25 0.703 0.7051 0.299% 0.375 0.689 0.6880 $-0.140%$ 0.5 0.679 0.6783 $-0.107%$ 0.625 0.674 0.6723 $-0.256%$ 0.75 0.67 0.6687 $-0.199%$ 0.875 0.668 0.6668 $-0.176%$ 1 0.668 0.6662 $-0.264%$

表7 本文方法与文献[32]方法的比较 Table 7 Comparisons between our method and Ref.[32]

本文将一种高阶奇异积分的直接计算法用于双边界积分方程，成功的解决了二维、三维断裂力学中应力强度因子的计算问题. 本文方法，不需要在整个裂纹面都采用非连续元，只在裂纹尖端采用非连续元就能得到精度较高的结果. 本文还构造了两种特殊的裂尖单元，用于体现裂尖附近的位移场特征. 无论是采用本文构造的裂尖单元，还是采用普通的裂尖单元，都能得到高精度的结果，本文方法对裂尖间断位移的插值函数依赖性不强.

本文方法主要是针对线弹性断裂问题，对于裂尖有塑性变形的情况，这将带来域积分项，传统边界单元法需要在塑性变形区域进行域内离散，这将部分丧失其降维的优势. 本文方法将考虑采用径向积分法将域积分转化为边界积分，从而不需要域内离散，但在源点处在裂纹边界时，出现的超强奇异积分的处理尚有难度. 这正是本文方法的发展方向，如若解决，不仅可适用于弹塑性断裂问题还可适用于非均质材料的断裂问题.

Li Jun, Feng Weizhe, Gao Xiaowei