王述红, 邱伟, 高红岩, 张紫杉
东北大学 资源与土木工程学院, 辽宁 沈阳 110819
收稿日期：2017-12-25
基金项目：国家自然科学基金资助项目(51474050，U1602232);地质灾害防治与地质环境保护国家重点实验室项目(SKLGP2014K011);辽宁省高等学校优秀人才支持计划项目(LN2014006)。
作者简介：王述红(1969-)，男，江苏泰州人，东北大学教授，博士生导师。

摘要：基于三角形网格, 对裂缝扩展过程中流形单元变化情况进行了深入研究, 从几何网格的角度对数值流形方法的连续与非连续统一处理方式进行解读.采用一阶覆盖函数, 推导出数值流形算法的权函数表达式, 建立局部位移函数.通过数值流形计算程序, 得出裂缝尖端位移, 并计算尖端应力强度因子.通过经典的中心裂纹板模型, 对数值流形位移法求得的尖端应力强度因子进行验证, 算例的数值解和解析解吻合度较高, 证明数值流形法计算裂缝扩展的准确性, 为裂纹扩展过程中尖端应力强度因子的求解提供了新的数值解法.
关键词：裂缝扩展几何网格局部位移函数强度因子数值流形位移法
Application of Numerical Manifold Displacement Method in Crack Propagation of Rock Mass
WANG Shu-hong, QIU Wei, GAO Hong-yan, ZHANG Zi-shan
School of Resources & Civil Engineering, Northeastern University, Shenyang 110819, China
Corresponding author: WANG Shu-hong, E-mail: wangshuhong@mail.neu.edu.cn
Abstract: Based on the finite element triangle mesh, the changing of manifold elements in crack propagation process was studied in depth, and the continuous and discontinuous unified processing of the numerical manifold method are interpreted from the perspective of geometric mesh. The first-order coverage function was used to derive the weight function expression of numerical manifold algorithm, so that the local displacement function can be established. Through the numerical manifold calculation program, the crack tip displacement was obtained, and the tip stress intensity factor was calculated and verified by the classical cracked plate model. The numerical result is in good agreement with the analytical solution, which proves the numerical manifold method is accuracy and can be used as a new numerical solution for solving the tip stress intensity factor in the crack propagation process.
Key words: crack propagationgeometric meshlocal displacement functionintensity factornumerical manifold displacement method
当岩体中有未贯穿岩体的结构面存在时, 在受到工程扰动的情况下, 很容易造成结构面扩展开裂, 形成次生结构面, 进一步切割岩体, 形成规模较大的非连续块体, 该过程为连续到非连续的变形破坏^[1].由于其复杂性, 现有的主流数值分析软件, 例如有限元在对裂缝扩展问题分析时, 需要不断更新网格, 使得前处理工作复杂.为了能够统一求解连续和非连续变形问题, Shi^[2-4]提出了数值流形方法, 利用数学覆盖和物理覆盖两套覆盖系统构建单元网格, 并独立定义位移函数, 各覆盖位移函数通过权函数逼近求解整体场函数^[5].国内外众多学者对数值流形位移法及相关程序不断改进, 将其成功应用于岩体大变形破坏分析以及裂缝扩展模拟等相关研究中.
本文利用数值流形方法覆盖技术, 构建裂缝尖端位移函数, 用于求解尖端位移, 结合尖端位移场求解应力强度因子, 实现裂缝扩展的判据表征.
1 裂缝尖端位移求解采用数值流形的方法对岩体裂纹扩展进行分析, 首先同样需要进行网格划分, 建立两套覆盖系统^[6].在裂缝扩展过程中, 数学网格保持不变, 裂缝扩展不断剖分数学网格，从而不断形成新的物理网格^[7-8].本文采用三角形数学网格构建数学覆盖, 对岩体平面单裂纹扩展进行计算分析.裂缝扩展过程中流形单元的变化见图 1.
图 1(Fig. 1)

图 1 裂缝扩展过程中流形单元的变化Fig.1 Changing of manifold units during crack propagation (a)—裂缝未扩展；(b)—裂缝开始扩展；(c)—裂缝进一步扩展；(d)—裂缝贯穿.

图 1a中为连续的块体材料, 块体边界对数学覆盖进行切割，形成物理网格, 各物理网格相互叠加形成流形单元^[9].块体中产生一条单裂缝, 当发生扩展时, 对应的流形单元也将发生变化, 如图 1b~图 1d.块体中产生的裂缝对块体进行了部分切割, 裂缝路径范围内的流形单元产生了新的节点, 因此裂缝切割形成的流形元可以产生不同的位移如图 1b、图 1c所示.当裂缝完全贯穿块体之后, 块体被分割成不连续的两部分, 裂隙路径范围内的任一流形单元分别具有不同的节点, 形成的两个新的块体可以独立自由移动, 如图 1d所示.
将三角形单元各节点位移场定义如下:

(1)

用权函数W(x, y)表示流形元单元的权函数, 则各节点对应的权函数为

(2)

以图 1a中流形单元1₁, 2₁, 3₁为例, 该流形单元具有三个数学覆盖V₁, V₂, V₃, 同时具有三个物理覆盖1₁, 2₁, 3₁.设物理覆盖上的位移函数分别为常函数(u₁, v₁), (u₂, v₂), (u₃, v₃), 对应权函数分别为W₁(x, y), W₂(x, y), W₃(x, y), 权函数由坐标系可直接求出.则该流形单元的位移函数可表示为

(3)

式中:
2 应力强度因子算法应力强度因子是用来表征裂纹端部应力场特征的物理量, 也是用来描述应力在裂纹端部的奇异性.断裂力学认为, 当裂纹尖端的应力强度因子达到材料的断裂韧度时, 裂纹将发生扩展^[10-11].因此对块体裂纹扩展的研究, 首先应计算裂纹尖端的强度因子.本文采用计算方法求应力强度因子.
对于平面裂纹扩展问题, 一般情况下, Ⅰ, Ⅱ两种裂纹同时存在, 二者对应的裂纹尖端的应力强度因子分别为K_Ⅰ, K_Ⅱ.综合相关文献, 给出K_Ⅰ, K_Ⅱ两种混合型裂纹扩展问题的位移场计算公式:

(4)

式中:, G为剪切模量; ，μ为泊松比.
利用数值流形方法求得裂纹尖端位移场的分布, 结合线弹性断裂力学应力强度因子的解析公式, 从而求出应力强度因子.
3 算例分析简单的有限大平板中心裂纹模型如图 2所示.根据应力强度因子手册, 应力强度因子解析解计算式为

(5)

图 2(Fig. 2)

图 2 有限大中心裂纹模型及裂缝尖端坐标系Fig.2 Finite central crack model and crack tip coordinate system (a)—有限大中心裂纹模型；(b)—裂缝尖端坐标系.

式中, C为相对无限大平板计算公式的修正系数,
以有限大平板中心裂纹模型作为边界条件, 创建数值流形法计算模型.采用三角形单元, 模型上下均布荷载σ=2 MN/m(采用无限密集的集中荷载来实现), 取材料的弹性模量E=2× 10⁵ MPa, 泊松比μ=0.25.模型宽度2b=6 m, 裂纹长度2a=2 m, 材料的断裂韧度K_IC=1 MN·m^1/2.模型建立完成后, 需要求解裂缝尖端位移.首先将平面直角坐标系转化成极坐标系, 引入以裂纹端点为原点的坐标系, 也称裂纹前缘坐标系, 如图 2所示.对于长度为2a的裂纹, 对于指定一点(x, y), 对应的极坐标是(r, θ), 如图 2所示, 坐标转化公式如下:

(6)

长度为2a的裂纹在数值流形方法中的位置如图 3所示, 以裂纹右侧为研究对象, 经过极坐标转化, 计算出随裂缝扩展时的尖端位移, 见表 1.
图 3(Fig. 3)

图 3 数值流形模型Fig.3 Numerical manifold model (a)—数学网格覆盖；(b)—流形元划分.

表 1(Table 1)

表 1 坐标及位移统计表Table 1 Table of coordinates and displacement 直角坐标 极坐标 y方向尖端位移/m σ=1 MN/m σ=2 MN/m (3.5, 5) (6.10, 0.96) 0.000 801 0.001 601 (3.6, 5) (6.16, 0.94) 0.000 801 0.001 601 (3.7, 5) (6.19, 0.93) 0.000 806 0.001 609 (3.8, 5) (6.28, 0.92) 0.000 806 0.001 609 (3.9, 5) (6.34, 0.90) 0.000 805 0.001 609 (4.0, 5) (6.40, 0.89) 0.000 805 0.001 609 (4.1, 5) (6.46, 0.88) 0.000 806 0.001 609 (4.2, 5) (6.52, 0.87) 0.000 805 0.001 608 (4.3, 5) (6.59, 0.86) 0.000 832 0.001 661 (4.4, 5) (6.66, 0.84) 0.000 832 0.001 661 (4.5, 5) (6.72, 0.83) 0.000 832 0.001 661

表 1 坐标及位移统计表 Table 1 Table of coordinates and displacement

在得到尖端位移后, 根据式(4), 计算裂缝扩展过程, 此过程通过Matlab程序计算实现, 各裂缝宽度下的K_Ⅰ值, 对于纯Ⅰ型裂纹, K_Ⅰ=0, G=8×10⁴ MPa, k=0.78.
令
则

(7)

根据式(7)分别求得裂纹扩展过程中的裂纹尖端强度因子, 见表 2.
表 2(Table 2)

表 2 Ⅰ型裂纹应力强度因子随裂缝宽度变化(σ=2 MN/m)Table 2 Variation of stress intensity factor of Ⅰ crack with crack width(σ=2 MN/m) 0.5 2.724 9 2.550 4 6.84 0.6 2.899 4 2.815 6 2.97 0.7 3.089 3 3.069 5 0.64 0.8 3.296 4 3.317 3 0.63 0.9 3.523 0 3.562 7 0.39 1.0 3.771 7 3.809 3 1.11 1.1 4.045 6 4.059 8 0.35 1.2 4.348 4 4.317 3 0.72 1.3 4.684 6 4.584 8 2.17 1.4 5.059 7 4.865 5 3.99 1.5 5.480 4 5.163 0 6.14

表 2 Ⅰ型裂纹应力强度因子随裂缝宽度变化(σ=2 MN/m) Table 2 Variation of stress intensity factor of Ⅰ crack with crack width(σ=2 MN/m)

由表 2可知, 随着裂缝不断扩展, 宽度不断增加, 尖端应力强度因子随之不断增大, 裂缝稳定扩展, 最终裂纹贯通.保持数学覆盖不变, 对表征裂缝的物理网格不断更新, 实现裂缝扩展, 见图 4.
图 4(Fig. 4)

图 4 中心裂纹扩展示意图Fig.4 Diagram of central crack propagation (a)—裂纹未扩展；(b)—裂纹开始扩展；(c)—裂纹进一步扩展；(d)—裂纹贯穿.

在拉伸荷载固定的条件下, 随着裂缝宽度a的变化, 采用数值解法求得的K_Ⅰ值与解析解求得的K_Ⅰ值误差成规律性变化, 即在裂缝宽度较小和较大时, 二者误差较大, 可能是因为裂缝较小或较大时, 数值流形网格大小(数学网格精度)与裂缝宽度不匹配.
当裂纹宽度a=1 m时, 其直角坐标为(4, 5), 转化为极坐标(6.403 1, 0.896 1).结合数值流形法计算得到的尖端裂纹位移, 通过设置测量点, 求得尖端位移, 在不同荷载情况下求得K_Ⅰ数值解, 如表 3所示.
表 3(Table 3)

表 3 Ⅰ型裂纹应力强度因子随拉伸荷载变化Table 3 Variation of stress intensity factor of Ⅰ crack with tensile load 0.5 0.952 3 0.952 3 0.00 1.0 1.897 6 1.904 6 0.37 1.5 2.845 2 2.856 9 0.41 2.0 3.792 8 3.809 3 0.43 2.5 4.740 4 4.761 5 0.44

表 3 Ⅰ型裂纹应力强度因子随拉伸荷载变化 Table 3 Variation of stress intensity factor of Ⅰ crack with tensile load

对比分析数值计算结果和应力强度因子手册计算得到的K_Ⅰ解析解, 二者基本接近，误差均在1 %范围内.随着荷载的不断增大误差逐渐稳定在0.45 %左右, 充分体现出基于数值流形位移法计算尖端应力强度因子的精确度和稳定性.
针对误差分析结果, 增加网格精度至12, 同时提高拉荷载至5 MN/m, 得到图 5所示的模拟结果.通过结果可以看出, 提高荷载和数学网格精度对裂缝扩展的连贯性和程序的收敛性均有所提高.
图 5(Fig. 5)

图 5 右侧裂缝扩展开裂示意图Fig.5 Diagram of crack propagation on the right (a)—右侧单裂缝模型；(b)—裂缝扩展开裂.

4 结论本文在数值流形法的基础上, 对单裂缝扩展过程进行模拟, 展现出裂缝扩展过程中的流形单元变化, 并建立单元的位移函数, 求解尖端位移.通过经典的中心裂纹板模型, 对数值流形位移法求得尖端应力强度因子进行验证, 发现随着裂缝宽度不断增加, 尖端应力强度因子随之不断增大; 宽度一定时采用数值解法求得的K_Ⅰ值与解析解求得的K_Ⅰ值误差成规律性变化, 在不同荷载情况下求得K_Ⅰ数值解和应力强度因子手册计算得到的K_Ⅰ解析解的误差均在1 %范围内, 数值解和解析解吻合度较高, 说明以数值流形法为基础求解裂纹尖端应力强度因子是一种新型的有效方法.
参考文献

[1]	Wang S H, Ni P P. Application of block theory modeling on spatial block topological identification to rock slope stability analysis[J].International Journal of Computational Methods, 2014, 11(1): 7–49.
[2]	Shi G H.Manifold method of material analysis[C]// The 9th Army Conference on Applied Mathematics and Computing.Minneapolis, 1991: 57-76.
[3]	Shi G H.Manifold method[C]// Proceedings of the First International Forum on DDA and Simulations of Discontinuous Media.Berkeley, 1996: 52-204.
[4]	Shi G H.Numerical manifold method[C]//Proceedings of the Second International Conference on Analysis of Discontinuous Deformation.Kyoto, 1997: 1-35.
[5]	王水林, 葛修润. 流形元方法在模拟裂纹扩展中的应用[J].岩石力学与工程学报, 1997, 16(5): 405–410. ( Wang Shui-lin, Ge Xiu-run. Application of manifold method in simulating crack propagation[J].Chinese Journal of Rock Mechanics and Engineering, 1997, 16(5): 405–410.DOI:10.3321/j.issn:1000-6915.1997.05.002)
[6]	Liu Z J, Zheng H. Two-dimensional numerical manifold method with multilayer covers[J].Science China Technological Sciences, 2016, 59(4): 515–530.DOI:10.1007/s11431-015-5907-z
[7]	Zheng H, Xu D D. New strategies for some issues of numerical manifold method in simulation of crack propagation[J].International Journal for Numerical Methods in Engineering, 2014, 97(13): 986–1010.DOI:10.1002/nme.v97.13
[8]	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–157.DOI:10.1016/j.cma.2015.07.001
[9]	Yang Y T, Zheng H. A three-node triangular element fitted to numerical manifold method with continuous nodal stress for crack analysis[J].Engineering Fracture Mechanics, 2016, 162: 51–75.DOI:10.1016/j.engfracmech.2016.05.007
[10]	任中俊, 陈跃欣. 有限大平板中心裂纹应力强度因子分析[J].矿业工程研究, 2013(28): 1–3. ( Ren Zhong-jun, Chen Yue-xin. Investigation of stress intensity factor of centre crack in finite plate[J].Mineral Engineering Research, 2013(28): 1–3.)
[11]	徐慧, 伍晓赞, 程仕平, 等. 复合裂纹的应力强度因子有限元分析[J].中南大学学报(自然科学版), 2007, 38(1): 79–83. ( Xu Hui, Wu Xiao-zan, Chen Shi-ping, et al. Finite element analysis of stress intensity factor in composite mode crack[J].Journal of Central South University(Science and Technology), 2007, 38(1): 79–83.DOI:10.3969/j.issn.1672-7207.2007.01.016)