土木与环境工程学报  2021, Vol. 43 Issue (5): 58-65   doi: 10.11835/j.issn.2096-6717.2020.142   PDF    
基于新型流形法的三维应力强度因子求解
祁勇峰 , 苏海东 , 龚亚琦     
长江科学院 材料与结构研究所; 水利部水工程安全和病害防治工程技术研究中心, 武汉 430010
收稿日期：2020-04-28
基金项目：国家自然科学基金（51779019、51809017）；中央级公益性科研院基本科研业务费（CKSF2020065/CL）
作者简介：祁勇峰(1978-), 男, 高级工程师, 主要从事水工结构数值分析和计算方法研究, E-mail: qiyf@mail.crsri.cn.
摘要：为保证在裂纹尖端具有最佳的逼近，裂纹尖端的解析解与其周边数值解联合应用的新型流形法用来求解应力强度因子，但仅限于平面问题的Ⅰ型和Ⅱ型裂纹。沿用解析解与数值解联合应用的思路，以三维穿透直裂纹为研究对象，在裂纹尖端引入Williams解析解级数，应用解析覆盖与周边数值覆盖联合求解三维应力强度因子，相对于其他数值方法而言计算精度高。推导相应的刚度矩阵和应变矩阵的具体表达式，通过典型算例验证了解决三维应力强度因子精确求解的有效性。
关键词：三维应力强度因子    数值流形方法    裂纹尖端Williams解析级数    解析覆盖    数值覆盖    
Computing 3D stress intensity factors based on new manifold method
Qi Yongfeng , Su Haidong , Gong Yaqi     
Material and Engineering Structure Department; Research Center on Water Engineering Safety and Disaster Prevention of Ministry of Water Resources, Yangtze River Scientific Research Institute, Wuhan 430010, P. R. China
Received: 2020-04-28
Foundation item: National Natural Science Foundation of China (No. 51779019, 51809017); Fundamental Research Funds for Central Public Welfare Research Institutes (No. CKSF2020065/CL)
Author brief: QI Yongfeng (1978-), senior engineer, main research interests: numerical analysis and calculation method of hydraulic structure, E-mail: qiyf@mail.crsri.cn.
Abstract: A new method is presented to compute the Stress Intensity Factors(SIF) by combining analytical solutions with numerical solutions in previous studies based on numerical manifold method, which obtains optimal approximation at the crack tip but only for the type Ⅰ and typeⅡ cracks in planar. Following the idea of combined application, Williams series is introduced at the crack tip and 3D SIF for three-dimensional straight through cracks is completed via combining analytical covers with numerical covers, which obtains higher computational precision than other numerical methods. Meanwhile, the corresponding formulas of the stiffness matrix and the strain matrix are deduced, and numerical examples shows the validity of the method in solving the exact solution of 3D stress intensity factor.
Keywords: 3D stress intensity factor (SIF)    numerical manifold method    Williams series for crack tip    analytical covers    numerical covers    

应力强度因子准则[1]是目前最常用的结构裂缝扩展准则之一，基于线弹性断裂力学的应力强度因子(Stress Intensity Factor，简称SIF)是用来表征裂纹尖端附近应力场和应变场强度，是控制裂纹尖端应力场和应变场强度的关键参数，在裂纹扩展分析中有着极其重要的地位。由于应力强度因子取决于外力的大小和分布、结构的几何条件以及裂缝的形状和位置，实际上只有少数问题存在解析解，对于复杂几何形状和加载条件的问题，只能通过数值方法来计算。目前裂缝扩展分析的主流数值计算方法有有限元法、扩展有限元法、数值流形法等。

有限元法[2-4]是目前分析裂缝等不连续问题的主要方法，为体现裂纹尖端(下简称裂尖)周围的应力集中和奇异性，往往需要在裂尖附近的复杂应力区布置高密度的单元网格，导致单元数目非常庞大；另外，在模拟裂缝扩展过程时，需不断重构网格，因此，有限元法对网格的要求和依赖性极大地影响了计算效率。扩展有限元法[5-7]通过在单元插值函数中引入能够反映裂缝面特性的不连续阶跃函数及反映裂尖局部特性的裂尖渐进位移场函数，裂缝可以穿过单元内部，裂缝扩展以后无需重新划分单元网格，采用同一网格就可以分析任意位置的裂缝问题，克服了常规有限元进行裂缝扩展分析的缺点，极大地简化了前处理工作。数值流形方法[8-11]引入不连续覆盖模拟裂缝，裂缝可以在网格内部穿过，巧妙地解决了常规有限元法裂缝面必须与单元边一致、裂缝扩展后需要重新划分网格的问题。相对于扩展有限元法设立不连续阶跃函数的方式而言，这种方式效果更好，在裂缝非常靠近单元边界时不会产生后者容易出现的数值误差。但无论是常规有限元法在裂尖布置细密网格的方式，还是扩展有限元法引入裂尖渐近位移场的方式，其主要目的都是为了提高裂尖附近的求解精度，从而提高应力强度因子的计算精度，这些方法都有改进的余地。即使是目前认为最适合于裂缝扩展分析的扩展有限元法，由于其只是使用了裂尖渐近位移场的部分特征函数，严格地说，还不能构成对裂尖位移场的最佳逼近。文献[12]提出的新型数值流形方法，实现了裂纹尖端的解析解与其周边数值解联合应用以求解应力强度因子，能够采用裂尖真实位移场的最佳逼近，并直接得到应力强度因子，计算精度高，但仅限于平面问题的Ⅰ型和Ⅱ型裂纹。

在上述研究的基础上，沿用解析解与数值解联合应用的思路，在裂纹尖端直接引入Williams解析解作为数学覆盖，应用裂纹尖端的解析解与周边数值解的三维流形覆盖联系技术，可直接得出裂纹尖端的三维应力强度因子，精度高且计算收敛快。

1 裂纹尖端位移场的Williams级数解

如图 1所示的裂纹尖端区域，Williams位移级数解

$ \begin{aligned} u &=\sum\limits_{i=0}^{m} \frac{1}{2} r^{i / 2}\left[a_{i} f_{i 1}^{x}(\theta)+b_{i} f_{i 2}^{x}(\theta)\right] ;\\ v &=\sum\limits_{i=0}^{m} \frac{1}{2} r^{i / 2}\left[a_{i} f_{i 1}^{y}(\theta)+b_{i} f_{i 2}^{y}(\theta)\right] ; \\ w &=\left\{\begin{array}{lr} c_{i} & i=0\\ \sum\limits_{i=1}^{m} r^{i-1 / 2}\left[c_{i} f_{i 3}^{z}(\theta)\right] & i=1, \cdots, m \end{array}\right. \end{aligned} $ (1)
图 1 裂纹尖端坐标系 Fig. 1 Coordinate system of crack tip

式中：$f_{i 1}^{x}(\theta)=c_{i 11}^{x} \cos \frac{i}{2} \theta+c_{i 12}^{x} \cos \left(\frac{i}{2}-2\right) \theta; c_{i 11}^{x}=$ $k + \frac{i}{2} + {(- 1)^i}; c_{i12}^x = - \frac{i}{2};f_{i2}^x(\theta) = c_{i21}^x\sin \frac{i}{2}\theta + $ $c_{i22}^x\sin \left({\frac{i}{2} - 2} \right)\theta; c_{i21}^x = - k - \frac{i}{2} + {(- 1)^i}; c_{i22}^x = \frac{i}{2}$；$f_{i1}^y(\theta) = c_{i11}^y\sin \frac{i}{2}\theta + c_{i12}^y\sin \left({\frac{i}{2} - 2} \right)\theta; c_{i11}^y = k - \frac{i}{2} - $ ${(- 1)^i}; c_{i12}^y = \frac{i}{2};f_{i2}^y(\theta) = c_{i21}^y\cos \frac{i}{2}\theta + $ $c_{i22}^y\cos \left({\frac{i}{2} - 2} \right)\theta; c_{i21}^y = k - \frac{i}{2} + {(- 1)^i}; c_{i22}^y = \frac{i}{2}$；$f_{i3}^z(\theta) = c_i^z\sin \left({i - \frac{1}{2}} \right)\theta; c_i^z = - \frac{2}{{i - 1/2}}(i = 1, \cdots, m)$。ai、bi、ci为待求系数; a0、b0、c0代表刚体位移；SIF与ai、bi、ci的关系为

$ K_{\mathrm{I}}=\sqrt{2 {\rm{ \mathsf{ π} }}} G a_{1}, K_{\mathrm{I I}}=\sqrt{2 {\rm{ \mathsf{ π} }}} G b_{1}, K_{\mathrm{II I}}=\sqrt{2 {\rm{ \mathsf{ π} }}} G{c_{1}} $ (2)
2 裂纹尖端位移场的解析解覆盖

在六面体网格内，将式(1)写成矩阵形式。

对于i=0:

$ \left\{\begin{array}{c} u \\ v \\ w \end{array}\right\}=\sum\limits_{j=1}^{8} w_{j}\left[\begin{array}{ccc} (k+1) / 2 & 0 & 0 \\ 0 & (k+1) / 2 & 0 \\ 0 & 0 & 1 \end{array}\right]\left\{\begin{array}{l} a_{0} \\ b_{0} \\ c_{0} \end{array}\right\} $ (3)

对于i≥1:

$ \begin{array}{c} \left\{\begin{array}{c} u \\ v \\ w \end{array}\right\}=\sum\limits_{j=1}^{8} w_{j} \sum\limits_{i=1}^{m} \frac{1}{2}\left[\begin{array}{lll} r^{\frac{i}{2}} & & \\ & r^{\frac{i}{2}} & \\ & & 2 r^{i-\frac{1}{2}} \end{array}\right] \rightarrow\\ \leftarrow\left[\begin{array}{ccc} f_{i 1}^{x}(\theta) & f_{i 2}^{x}(\theta) & 0 \\ f_{i 1}^{y}(\theta) & f_{i 2}^{y}(\theta) & 0 \\ 0 & 0 & f_{i 3}^{z}(\theta) \end{array}\right]\left\{\begin{array}{l} a_{i} \\ b_{i} \\ c_{i} \end{array}\right\} \end{array} $ (4)

将二维数值覆盖强制约束方法[13]推至三维，形成独立的三维解析覆盖，权函数wj均为结点1的权函数w1。其中，m为覆盖函数项数，下同。

对于第i项(i≥1)，其应变子矩阵为

$ \boldsymbol{B}_{i}=\left[\begin{array}{ccc} \frac{\partial}{\partial x} & & \\ & \frac{\partial}{\partial y} & \\ & & \frac{\partial}{\partial z} \\ \frac{\partial}{\partial y} & \frac{\partial}{\partial x} & \\ & \frac{\partial}{\partial z} & \frac{\partial}{\partial y} \\ \frac{\partial}{\partial z} & & \frac{\partial}{\partial x} \end{array}\right] \left(\frac{1}{2} w_{j} r^{i / 2}\left[\begin{array}{ccc} f_{i 1}^{x}(\theta) & f_{i 2}^{x}(\theta) & 0 \\ f_{i 1}^{y}(\theta) & f_{i 2}^{y}(\theta) & 0 \\ 0 & 0 & 2 r^{i / 2-1 / 2} f_{i 3}^{z}(\theta) \end{array}\right]\right) $ (5)

以f(θ)代表fi1x(θ)、fi2x(θ)、fi1y(θ)和fi2y(θ)、fi3z(θ)，由于r、f(θ)与坐标z无关，因此，式(5)中有关r·f(θ)多项式对z的偏导数为0，则式(5)为

$ \boldsymbol{B}_{i}=\left[\begin{array}{ccc} \frac{\partial}{\partial x}\left(\frac{1}{2} w_{j} r^{i / 2} f_{i 1}^{x}(\theta)\right) &\frac{\partial}{\partial x}\left(\frac{1}{2} w_{j} r^{i / 2} f_{i 2}^{x}(\theta)\right)&0\\ \frac{\partial}{\partial y}\left(\frac{1}{2} w_{j} r^{i / 2} f_{i 1}^{y}(\theta)\right) &\frac{\partial}{\partial y}\left(\frac{1}{2} w_{j} r^{i / 2} f_{i 2}^{y}(\theta)\right)&0\\ 0&0&r^{i+1 / 2} f_{i 3}^{z}(\theta) \frac{\partial}{\partial z}\left(w_{j}\right)\\ \frac{\partial}{\partial y}\left(\frac{1}{2} w_{j} r^{i / 2} f_{i 1}^{x}(\theta)\right)+\frac{\partial}{\partial x}\left(\frac{1}{2} w_{j} r^{i / 2} f_{i 1}^{y}(\theta)\right) & \frac{\partial}{\partial y}\left(\frac{1}{2} w_{j} r^{i / 2} f_{i 2}^{x}(\theta)\right)+\frac{\partial}{\partial x}\left(\frac{1}{2} w_{j} r^{i / 2} f_{i 2}^{y}(\theta)\right)&0\\ \frac{1}{2} r^{j / 2} f_{i 1}^{y}(\theta) \frac{\partial}{\partial z}\left(w_{j}\right) & \frac{1}{2} r^{i / 2} f_{i 2}^{y}(\theta) \frac{\partial}{\partial z}\left(w_{j}\right) & \frac{\partial}{\partial y}\left(w_{j} r^{i-1 / 2} f_{i 3}^{z}(\theta)\right) \\ \frac{1}{2} r^{j / 2} f_{i 1}^{x}(\theta) \frac{\partial}{\partial z}\left(w_{j}\right) & \frac{1}{2} r^{i / 2} f_{i 2}^{x}(\theta) \frac{\partial}{\partial z}\left(w_{j}\right) & \frac{\partial}{\partial x}\left(w_{j} r^{i-1 / 2} f_{i 3}^{z}(\theta)\right) \end{array}\right] $ (6)

将各项偏导数展开

$ \begin{aligned} &\frac{\partial}{\partial x}\left(\frac{1}{2} w_{j} r^{i / 2} f(\theta)\right)=\frac{\partial w_{j}}{\partial x}\left(\frac{1}{2} r^{i / 2} f(\theta)\right)+\\ &w_{j} \frac{\partial}{\partial x}\left(\frac{1}{2} r^{i / 2} f(\theta)\right) ;\\ &\frac{\partial}{\partial y}\left(\frac{1}{2} w_{j} r^{i / 2} f(\theta)\right)=\frac{\partial w_{j}}{\partial y}\left(\frac{1}{2} r^{i / 2} f(\theta)\right)+\\ &w_{j} \frac{\partial}{\partial y}\left(\frac{1}{2} r^{i / 2} f(\theta)\right) ;\\ &\frac{\partial}{\partial x}\left(w_{j} r^{i-1 / 2} f(\theta)\right)=\frac{\partial w_{j}}{\partial x}\left(r^{i-1 / 2} f(\theta)\right)+\\ &w_{j} \frac{\partial}{\partial x}\left(r^{i-1 / 2} f(\theta)\right) ;\\ &\frac{\partial}{\partial y}\left(w_{j} r^{i-1 / 2} f(\theta)\right)=\frac{\partial w_{j}}{\partial y}\left(r^{i-1 / 2} f(\theta)\right)+\\ &w_{j} \frac{\partial}{\partial y}\left(r^{i-1 / 2} f(\theta)\right)。\end{aligned} $ (7)

式中：

$ \begin{aligned} &\frac{\partial}{\partial x}\left(\frac{1}{2} r^{i / 2} f(\theta)\right)=\frac{\partial}{\partial r}\left(\frac{1}{2} r^{i / 2} f(\theta)\right) \frac{\partial r}{\partial x}+\frac{\partial}{\partial \theta} ; \\ &\left(\frac{1}{2} r^{i / 2} f(\theta)\right) \frac{\partial \theta}{\partial x}=\frac{i}{4} r^{(i / 2-1)} f(\theta) \frac{\partial r}{\partial x}+\frac{1}{2} r^{i / 2} \frac{\partial f(\theta)}{\partial \theta} \frac{\partial \theta}{\partial x} ;\\ &\frac{\partial}{\partial y}\left(\frac{1}{2} r^{i / 2} f(\theta)\right)=\frac{\partial}{\partial r}\left(\frac{1}{2} r^{i / 2} f(\theta)\right) \frac{\partial r}{\partial y}+\frac{\partial}{\partial \theta} ; \\ &\left(\frac{1}{2} r^{i / 2} f(\theta)\right)+\frac{\partial \theta}{\partial y}=\frac{i}{4} r^{(i / 2-1)} f(\theta) \frac{\partial r}{\partial y}+ \\ &\frac{1}{2} r^{i / 2} \frac{\partial f(\theta)}{\partial \theta} \frac{\partial \theta}{\partial y};\\ &\frac{\partial}{\partial x}\left(r^{i-1 / 2} f(\theta)\right)=\frac{\partial}{\partial r}\left(r^{i-1 / 2} f(\theta)\right) \frac{\partial r}{\partial x}+\frac{\partial}{\partial \theta}\left(r^{i-1 / 2} f(\theta)\right) ; \\ &\frac{\partial \theta}{\partial x}=(i-1 / 2) r^{(i-3 / 2)} f(\theta) \frac{\partial r}{\partial x}+r^{i-1 / 2} \frac{\partial f(\theta)}{\partial \theta} \frac{\partial \theta}{\partial x}; \\ &\frac{\partial}{\partial y}\left(r^{i-1 / 2} f(\theta)\right)=\frac{\partial}{\partial r}\left(r^{i-1 / 2} f(\theta)\right) \frac{\partial r}{\partial y}+\frac{\partial}{\partial \theta}\left(r^{-1 / 2} f(\theta)\right) ; \\ &\frac{\partial \theta}{\partial y}=(i-1 / 2) r^{(i-3 / 2)} f(\theta) \frac{\partial r}{\partial y}+r^{j-1 / 2} \frac{\partial f(\theta)}{\partial \theta} \frac{\partial \theta}{\partial y}。\end{aligned} $

具体到fi1x(θ)、fi2x(θ)、fi1y(θ)和fi2y(θ)，关于θ的偏导数为

$ \begin{aligned} &\qquad\frac{\partial f_{i 1}^{x}(\theta)}{\partial \theta}=-c_{i 11}^{x} \frac{i}{2} \sin \frac{i}{2} \theta+c_{i 12}^{x}\left(2-\frac{i}{2}\right) \cdot \\ &\sin \left(\frac{i}{2}-2\right) \theta ; \\ &\qquad\frac{\partial f_{i 2}^{x}(\theta)}{\partial \theta}=c_{i 21}^{x} \frac{i}{2} \cos \frac{i}{2} \theta+c_{i 22}^{x}\left(\frac{i}{2}-2\right)\cdot \\ &\cos \left(\frac{i}{2}-2\right) \theta; \\ &\qquad\frac{\partial f_{i 1}^{y}(\theta)}{\partial \theta}=c_{i 11}^{y} \frac{i}{2} \cos \frac{i}{2} \theta+c_{i 12}^{y}\left(\frac{i}{2}-2\right) \cdot \\ &\cos \left(\frac{i}{2}-2\right) \theta ; \\ &\qquad\frac{\partial f_{i 2}^{y}(\theta)}{\partial \theta}=-c_{i 21}^{y} \frac{i}{2} \sin \frac{i}{2} \theta+c_{i 22}^{y}\left(2-\frac{i}{2}\right)\cdot \\ &\sin \left(\frac{i}{2}-2\right) \theta ; \\ &\qquad\frac{\partial f_{i 3}^{z}(\theta)}{\partial \theta}=-2 \cos \left(i-\frac{1}{2}\right) \theta。\end{aligned} $

另外，

$ \begin{aligned} & \frac{\partial r}{\partial x}=\frac{x-x_{0}}{r} ; \frac{\partial r}{\partial y}=\frac{y-y_{0}}{r} ; \frac{\partial \theta}{\partial x}=\frac{-\left(y-y_{0}\right)}{r^{2}} ; \\ \frac{\partial \theta}{\partial y}=& \frac{x-x_{0}}{r^{2}}。\end{aligned} $

刚度子矩阵

$ \boldsymbol{K}_{\mathrm{ij}}=\int_{V} \boldsymbol{B}_{i}^{\mathrm{T}} \boldsymbol{D} \boldsymbol{B}_{j} \mathrm{~d} V \quad i, j=0,1, \cdots \cdots, m $ (8)

式中：D为弹性矩阵；V为流形元体积。

3 裂纹周边网格的解析解与数值解覆盖

在裂纹周边的各网格内，位移统一表示为式(9)，其中i≥1。

$ \begin{aligned} &\left\{\begin{array}{l} u \\ v \\ w \end{array}\right\}=\sum\limits_{j=1}^{J} w_{j 1} \sum\limits_{i=1}^{m} \frac{1}{2} r^{i / 2}\left[\begin{array}{ccc} f_{i 1}^{x}(\theta) & f_{i 2}^{x}(\theta) & 0 \\ f_{i 1}^{y}(\theta) & f_{i 2}^{y}(\theta) & 0 \\ 0 & 0 & 2 r^{i / 2-1 / 2} f_{i 3}^{z}(\theta) \end{array}\right]\left\{\begin{array}{l} a_{i} \\ b_{i} \\ c_{i} \end{array}\right\}+\\ &\sum\limits_{j 2=1}^{L} w_{j 2} \sum\limits_{n=0}^{p} \sum\limits_{k=0}^{n} \sum\limits_{l=0}^{n-k}\left[\begin{array}{llll} x^{n-k-l} y^{l} z^{k} & \\ & x^{n-k-l} y^{l} z^{k} & \\ & & x^{n-k-l} y^{l} z^{k} \end{array}\right]\left\{\begin{array}{l} d_{n k l} \\ e_{n k l} \\ f_{n k l} \end{array}\right\} \end{aligned} $ (9)

式中：dnkl、enkl、fnkl为待求的多项式系数；n、k、l为多项式阶次(p为多项式的最高阶次)；wj1、wj2为相应于各结点的权函数；J和L的个数根据网格而定，但保持J+L=8。显然，取L=0则退化到解析解覆盖式(4)。

如图 2所示(图中窄条为部分重叠区域，长方体区域为独立覆盖，阴影为裂纹)，采用强制约束的方法，将网格结点28、40、39、33、34、45、46约束到结点27，即令

$ \begin{aligned} \left\{\begin{array}{l} u \\ v \\ w \end{array}\right\}_{27}=&\left\{\begin{array}{l} u \\ v \\ w \end{array}\right\}_{28}=\left\{\begin{array}{l} u \\ v \\ w \end{array}\right\}_{33}=\left\{\begin{array}{l} u \\ v \\ w \end{array}\right\}_{34}=\left\{\begin{array}{l} u \\ v \\ w \end{array}\right\}= \\ &\left\{\begin{array}{l} u \\ v \\ w \end{array}\right\}_{40}=\left\{\begin{array}{l} u \\ v \\ w \end{array}\right\}_{45}=\left\{\begin{array}{c} u \\ v \\ w \end{array}\right\}_{46} \end{aligned} $
图 2 六面体数学网格中的解析解覆盖和数值解覆盖 Fig. 2 Analytical covers and numerical covers in hexahedral mathematical mesh

根据形函数$\sum\limits_{j=1}^{8} w_{j}=1$的关系，从而在该网格内形成包含裂纹尖端的独立的解析覆盖。

对于周边的其他网格，比如：网格1-2-14-13-7-8-20-19中的结点均为数值解的独立覆盖；解析覆盖的相邻网格39-40-52-51-45-46-58-57中的结点39、40、45、46为解析结点，其他4个结点为数值结点；而裂纹自由端所在网格25-26-38-37-31-32-44-43在竖直方向采用了两个独立覆盖，实现裂纹两边的独立运动。

对于解析覆盖结点j处的第i项，其应变矩阵为

$ \boldsymbol{B}_{ji}=\left[\begin{array}{ccc} \frac{\partial}{\partial x} & & \\ & \frac{\partial}{\partial y} & \\ & & \frac{\partial}{\partial z} \\ \frac{\partial}{\partial y} & \frac{\partial}{\partial x} & \\ & \frac{\partial}{\partial z} & \frac{\partial}{\partial y} \\ \frac{\partial}{\partial z} & & \frac{\partial}{\partial x} \end{array}\right] \left(\frac{1}{2} w_{j} r^{i / 2}\left[\begin{array}{ccc} f_{i 1}^{x}(\theta) & f_{i 2}^{x}(\theta) & 0 \\ f_{i 1}^{y}(\theta) & f_{i 2}^{y}(\theta) & 0 \\ 0 & 0 & 2 r^{i / 2-1 / 2} f_{i 3}^{z}(\theta) \end{array}\right]\right) $ (10)

而对于数值覆盖结点t处，设多项式的总项数为q(为$\sum\limits_{n=0}^{p} \sum\limits_{k=0}^{n} \sum\limits_{l}^{n-k}$的项数之和)，其应变矩阵为

$ \boldsymbol{B}_{tq}=\left[\begin{array}{ccc} \frac{\partial}{\partial x} & & \\ & \frac{\partial}{\partial y} & \\ & & \frac{\partial}{\partial z} \\ \frac{\partial}{\partial y} & \frac{\partial}{\partial x} & \\ & \frac{\partial}{\partial z} & \frac{\partial}{\partial y} \\ \frac{\partial}{\partial z} & & \frac{\partial}{\partial x} \end{array}\right] \left(W_{t} \sum\limits_{n=0}^{p} \sum\limits_{k=0}^{n} \sum\limits_{l=0}^{n-k}\left[\begin{array}{lll} x^{n-k-l} y^{l} z^{k} & & \\ & & x^{n-k-l} y^{l} z^{k} & \\ & & & x^{n-k-l} y^{l} z^{k} \end{array}\right]\right) $ (11)

令f(x, y, z)=xn-k-lylzk，则式(11)可写为

$ \boldsymbol{B}_{t q}=\sum\limits_{n=0}^{p} \sum\limits_{k=0}^{n} \sum\limits_{l=0}^{n-k}\left[\begin{array}{ccc} \frac{\partial}{\partial x}\left(W_{t} f(x, y, z)\right) & 0 & 0 \\ 0 & \frac{\partial}{\partial y}\left(W_{t} f(x, y, z)\right) & 0 \\ 0&0&\frac{\partial}{\partial z}\left(W_{t} f(x, y, z)\right) \\ \frac{\partial}{\partial y}\left(W_{t} f(x, y, z)\right) & \frac{\partial}{\partial x}\left(W_{t} f(x, y, z)\right) & 0 \\ 0 & \frac{\partial}{\partial z}\left(W_{t} f(x, y, z)\right) & \frac{\partial}{\partial y}\left(W_{t} f(x, y, z)\right) \\ \frac{\partial}{\partial z}\left(W_{t} f(x, y, z)\right) & 0 & \frac{\partial}{\partial x}\left(W_{t} f(x, y, z)\right) \end{array}\right] $ (12)

解析覆盖与数值覆盖相关的刚度子矩阵为

$ \boldsymbol{K}_{j i t q}=\int_{\mathrm{V}} \boldsymbol{B}_{j i} \boldsymbol{D} \boldsymbol{B}_{t q} \mathrm{~d} V \text { 或 } \boldsymbol{K}_{tqj i}=\int_{\mathrm{V}} \boldsymbol{B}_{i q} \boldsymbol{D} \boldsymbol{B}_{j i} \mathrm{~d} V $ (13)

由式(10)和式(12)可见，刚度矩阵积分中具有非多项式的函数，基于多形式函数的单纯形积分公式无法应用，必须采用数值积分。因此，采用四面体区域的Hammer积分[14]，将流形元边界上的三角片与其形心相连形成四面体，然后再进行积分计算。

4 算例分析

以含边界裂纹的无限长柱体为例，考虑两种不同类型荷载。

4.1 两端受均布拉力

图 3所示为包含边界裂缝的两端受均布拉力的无限长柱体，柱体的宽度为w，应力强度因子理论值为KI=$\sqrt {\pi a} f\left({a/w} \right)$[15]，其中f(a/w)=1.12-0.23(a/w)+10.6(a/w)2-21.7(a/w)3+30.4(a/w)4，取裂纹长度a=0.05 m，柱体宽度w=0.4 m，模型的柱体长度可取h=2 m，柱体厚度t取单位宽度，弹性模量E=200 000 kPa，泊松比为0.167，均布荷载P为30 kN/m2。KI理论值为1.45 kPa/m1/2。

图 3 两端受均布拉力的无限长柱体内的边界裂缝 Fig. 3 Boundary crack in infinite long column with both ends bearing uniform force

整体的流形元网格如图 4所示，共划分6个独立覆盖和9个部分重叠覆盖(窄条区域)，每个独立覆盖的大小基本相同。裂纹所在的独立覆盖区域大小为0.65 m×0.2 m×1 m(长×宽×厚)，柱体底面施加法向约束。

图 4 部分重叠覆盖的流形元网格 Fig. 4 Manifold method mesh based on partially overlapping covers

应力强度因子的计算结果见表 1，考虑了裂纹尖端所在的独立覆盖取Williams级数的不同项数，独立覆盖周边的数值解覆盖可取多项式的不同阶数。

表 1 应力强度因子 Table 1 Stress intensity factors

随着覆盖函数阶数的增加，三维应力强度因子计算值基本接近理论值。Williams级数的阶数以及周边数值覆盖阶数对计算值有较大影响。

1) Williams级数的阶数(m)影响最大。当m≤3，则KI与理论值相差较大；当m≥4时，KI接近理论值。

2) 周边数值覆盖阶数升高有利于提高解的精度。当周边数值解取1阶时，KI的计算值普遍小于取2阶的情况，当m≥4时，KI与理论值接近，但仍有差距，仅当m取为7时才达到1.42，与理论值1.45最为接近。而当周边数值覆盖取2阶，m≥7时，计算值与理论值基本一致。

前期平面问题研究表明，在大的覆盖中单纯依靠提高覆盖函数阶次的方法往往会带来计算结果的振荡跳跃。反之，如果采用较小的覆盖而用相对简单的低阶多项式，则可以更好地逼近实际复杂的分布情况。表 2的计算结果也说明，三维问题中，基于大覆盖，仅仅采用提高级数解及相邻覆盖函数的阶数的做法来提高计算精度，计算也表现出一定的不稳定性，要取得满意的计算精度，裂纹尖端及其周边的覆盖函数的阶数不小于7阶。

表 2 应力强度因子(n=2) Table 2 Stress intensity factors

考虑到裂纹所在的独立覆盖区域较大，因此，将裂纹所在的独立覆盖进行局部加密，将裂纹尖端覆盖分别加密1倍及2倍，采用局部覆盖加密技术[16]，重新计算应力强度因子，如表 2所示。

表 2结果表明，当裂纹尖端独立覆盖加密1倍后、Williams级数的阶数≥4或当裂纹尖端独立覆盖加密2倍后，Williams级数的阶数≥3时，KI与理论解十分接近，且随着Williams级数阶数的提高，计算结果趋于稳定。

4.2 两端受剪切荷载作用

考虑三维荷载作用下的算例。如图 5所示，无限长柱沿z方向施加均布荷载Q=100 kN/m2，裂纹长度a=0.02 m，柱体尺寸以及材料参数同上，其应力强度因子理论值为KⅢ=$Q\sqrt {\pi c} $[15]，其中$c = \frac{a}{w}$。KⅢ的理论值为2.51。

图 5 两端受剪力的无限长柱体内的边界裂缝 Fig. 5 Boundary crack in infinite long column with both ends bearing shear force

流形元网格如图 6所示，共划分18个独立覆盖(方块区域)和25个部分重叠覆盖(窄条区域)，每个独立覆盖的大小基本相同。裂纹所在的独立覆盖区域大小为0.2 m×0.2 m×1 m(长×宽×厚)，柱体底面施加法向约束。应力强度因子的计算结果见表 3。

图 6 部分重叠覆盖的流形元网格 Fig. 6 Manifold method mesh based on partially overlapping covers

表 3 应力强度因子 Table 3 Stress intensity factors

表 3计算结果表明：采用图 6所示计算网格，当m≥7时，周边数值覆盖阶数取2、3的多项式阶数时，计算结果与理论值比较符合。当m≤7，计算结果与理论值有一定差别，局部数值出现跳跃，表明裂纹附近的网格还没有达到足够的密度。当网格加密一倍后，周边数值覆盖阶数均取2阶，当m≥7时，计算值与理论值十分接近，且随着阶数的提高，计算结果趋于稳定。

以上算例验证了三维裂缝计算公式和程序的正确性，表明裂纹尖端解析解覆盖和周边数值解覆盖联合应用求解三维线弹性断裂力学问题可行。与常规有限元方法相比，无需在裂纹尖端布置细密的网格，计算精度高，收敛相对较快。

裂纹尖端独立覆盖的密度、解析覆盖的级数以及相邻数值覆盖的阶数是影响应力强度因子计算精度的重要因素，但在保证独立覆盖有一定密度的情况下，提高与独立覆盖相邻数值覆盖的阶数可以得到应力强度因子的精确解。

裂纹尖端独立覆盖的合理布置对应力强度因子的计算精度及稳定性有一定的影响，因此，下一步要重点研究裂纹尖端附近的覆盖自动布置及密度问题，以保证方法的收敛性，便于开展三维裂缝扩展的动态模拟研究。

5 结论

将裂纹尖端解析解覆盖和周边数值解覆盖联合应用，分析三维线弹性断裂力学问题，得到以下主要结论：

1) 在包含裂纹尖端的解析覆盖中，应用裂纹尖端附近的Williams位移解析解作为覆盖函数，并采用高阶多项式三维覆盖函数与解析覆盖的条形连接技术，实现了在解析覆盖中直接求得裂纹尖端的三维应力强度因子。

2) 典型的张开型和撕开型的裂纹算例表明，应力强度因子的计算精度较高。鉴于三维裂缝扩展问题的复杂性，裂纹尖端周边数值覆盖阶数以及独立覆盖网格密度对应力强度因子计算精度的影响较二维问题更大。因此，协调独立覆盖密度、阶数与周边三维数值覆盖阶数的关系，来保证高精度求解收敛性的快速、稳定是下一步研究的重点。

考虑到解析级数是裂尖附近真实位移场的最佳逼近，相比其他方法而言，可以认为该方法在应力强度因子求解方面逼近效果更好、收敛更快，同时，由于网格布置根据不同区域的精度要求，只在裂尖附近进行覆盖加密，因此，相比采用均匀网格的扩展有限元而言，计算效率将有所提高，可以实现大规模计算。另外，应力强度因子SIF本身就是裂尖解析级数的未知数，在求解系统方程组时一并得到，而不需要像其它方法那样通过所谓的“直接”法或“间接”法来推求，不仅方便，而且不会引入额外误差，这也是该方法的优势所在。

该方法可以同时求解Ⅰ型、Ⅱ型、Ⅲ型(撕开型)裂纹的应力强度因子，应用复合型裂纹扩展准则就可以判断其是否继续开裂，因此，该方法在三维裂缝扩展的动态模拟方面极具应用前景。

参考文献
[1]
PEREZ N. Linear-elastic fracture mechanics[M]. Cham: Springer International Publishing, 2016: 79-130.
[2]
FAGEEHI Y A, ALSHOAIBI A M. Nonplanar crack growth simulation of multiple cracks using finite element method[J]. Advances in Materials Science and Engineering, 2020, 2020: 1-12.
[3]
ALSHOAIBI M. Finite element simulation of crack growth path and stress intensity factors evaluation in linear elastic materials[J]. Jounal of Computational and Applied Research in Mechanical Engineering, 2019.
[4]
HAJIAN M, MORADI M. An improved approach for computation of stress intensity factors using the finite element method[J]. Engineering Analysis with Boundary Elements, 2019(1): 54-63.
[5]
李录贤, 王铁军. 扩展有限元法(XFEM)及其应用[J]. 力学进展, 2005, 35(1): 5-20.
LI L X, WANG T J. The extended finite element method and its applications: a review[J]. Advances in Mechanics, 2005, 35(1): 5-20. (in Chinese) DOI:10.3321/j.issn:1000-0992.2005.01.002
[6]
GINER E, SUKUMAR N, DENIA F D, et al. Extended finite element method for fretting fatigue crack propagation[J]. International Journal of Solids and Structures, 2008, 45(22/23): 5675-5687.
[7]
ELGUEDJ T, DE SAINT MAURICE R P, COMBESCURE A, et al. Extended finite element modeling of 3D dynamic crack growth under impact loading[J]. Finite Elements in Analysis and Design, 2018, 151: 1-17. DOI:10.1016/j.finel.2018.08.001
[8]
石根华. 数值流形方法与非连续变形分析[M]. 裴觉民, 译. 北京: 清华大学出版社, 1997.
SHI G H. Numerical manifold method and discontinuous deformation analysis[M]. PEI J M, translated. Beijing: Tsinghua University Press, 1997.
[9]
YANG S K, CAO M S, REN X H, et al. 3D crack propagation by the numerical manifold method[J]. Computers & Structures, 2018, 194: 116-129.
[10]
ZHANG G X, LI X, LI H F. Simulation of hydraulic fracture utilizing numerical manifold method[J]. Science China Technological Sciences, 2015, 58(9): 1542-1557. DOI:10.1007/s11431-015-5901-5
[11]
ZHENG H, LIU F, DU X L. Complementarity problem arising from static growth of multiple cracks and MLS-based numerical manifold method[J]. Computer Methods in Applied Mechanics and Engineering, 2015, 295: 150-171. DOI:10.1016/j.cma.2015.07.001
[12]
苏海东, 祁勇峰, 龚亚琦. 裂纹尖端解析解与周边数值解联合求解应力强度因子[J]. 长江科学院院报, 2013, 30(6): 83-89.
SU H D, QI Y F, GONG Y Q. Compute stress intensity factors via combining analytical solutions around crack tips with surrounding numerical solutions[J]. Journal of Yangtze River Scientific Research Institute, 2013, 30(6): 83-89. (in Chinese) DOI:10.3969/j.issn.1001-5485.2013.06.019
[13]
祁勇峰, 苏海东, 崔建华. 部分重叠覆盖的数值流形方法初步研究[J]. 长江科学院院报, 2013, 30(1): 65-70.
QI Y F, SU H D, CUI J H. Preliminary study on numerical manifold method with partially overlapping covers[J]. Journal of Yangtze River Scientific Research Institute, 2013, 30(1): 65-70. (in Chinese) DOI:10.3969/j.issn.1001-5485.2013.01.013
[14]
王勖成. 有限单元法[M]. 北京: 清华大学出版社, 2003.
WANG X C. Finite element method[M]. Beijing: Tsinghua University Press, 2003.
[15]
中国航空研究院. 应力强度因子手册[M]. 北京: 科学出版社, 1993.
Chinese Aeronautical Establishment. Stress intensity factors handbook[M]. Beijing: Science Press, 1993. (in Chinese)
[16]
苏海东, 龚亚琦, 颉志强, 等. 基于矩形独立覆盖初步实现结构静力分析的自动计算[J]. 长江科学院院报, 2016, 33(2): 144-150.
SU H D, GONG Y Q, XIE Z Q, et al. Preliminary implementation of automatic computation for static analysis of structures using NMM based on independent rectangular covers[J]. Journal of Yangtze River Scientific Research Institute, 2016, 33(2): 144-150. (in Chinese)