﻿ 基于单元破裂的岩石裂纹扩展模拟方法
«上一篇
 文章快速检索 高级检索

 力学学报  2015, Vol. 47 Issue (1): 105-118  DOI: 10.6052/0459-1879-14-239 0

### 引用本文 [复制中英文]

[复制中文]
Wang Jie, Li Sihai, Zhang Qingbo. SIMULATION OF CRACK PROPAGATION OF ROCK BASED ON SPLITTING ELEMENTS[J]. Chinese Journal of Ship Research, 2015, 47(1): 105-118. DOI: 10.6052/0459-1879-14-239.
[复制英文]

### 文章历史

2014-08-18 收稿
2014-10-07 录用
2014–11–18 网络版发表

1. 中国科学院力学研究所, 北京 100190;
2. 铁道第三勘察设计院集团有限公司, 天津 300142

1 连续单元模型构建

 图 1 三棱柱单元的离散弹簧系统 Fig.1 Spring system of triangular prism element

 $\label{eq1}\left.\begin{array}{l} u_{7i} = \dfrac{u_{1i} + u_{4i} }{2},\ \ u_{8i} = \dfrac{u_{2i} + u_{5i}}{2},\ \ u_{9i} = \dfrac{u_{3i} + u_{6i} }{2} \\[2mm] u_{Mi} = \dfrac{u_{1i} + u_{2i} + u_{3i} }{3},\ \ u_{Ni} = \dfrac{u_{4i} + u_{5i}+ u_{6i} }{3} \\ \end{array} \right\}$ (1)

 $\label{eq2}\left.\begin{array}{l} \Delta u_{xi} = u_{8i} - u_{7i} \\\Delta u_{yi} = u_{9i} - mu_{7i} -nu_{8i} \\ \Delta u_{zi} = u_{Mi} - u_{Ni} \\ \end{array} \right\}$ (2)

 $\label{eq3}m = \vert 8P\vert / \vert 78\vert ,\ \ n = \vert 7P\vert / \vert78\vert$ (3)

 $\label{eq4} a^e = [u_{7x} u_{7y} u_{8x} u_{8y} u_{9x} u_{9y}]$ (4)

 $u=[{{N}_{7}}{{N}_{8}}{{N}_{9}}]{{a}^{e}}=N{{a}^{e}}$ (5)

 $\label{eq6}\left.\begin{array}{l} a_7 = x_8 y_9 - x_9 y_8 ,b_7 = y_8 - y_9 ,c_7 = - x_8 + x_9 \\ a_8 = x_9 y_7 - x_7 y_9 ,b_8 = y_9 - y_7 ,c_8 = - x_9 + x_7 \\ a_9 = x_7 y_8 - x_8 y_7 ,b_9 = y_7 - y_8 ,c_9 = - x_7 + x_8 \\ \end{array} \right\}$ (6)

 $\label{eq7}\left.\begin{array}{l} a_1 = ab,b_1 = - b,c_1 = - a + na \\ a_2 = 0,b_2 = b,c_2 = - na \\ a_3 = 0,b_3 = 0,c_3 = a \\ \end{array} \right\}$ (7)

 $\label{eq8} \varepsilon = {Lu} = {LNa}^e = [B_7 B_8 B_9] a^e$ (8)

 $\label{eq9} B_i = \dfrac{1}{2A}\left[{\begin{array}{cc} b_i &0 \\ 0&c_i \\ c_i &b_i \\ \end{array}} \right]\ \ (i = 7,8,9)$ (9)

 $\label{eq10}\Delta \varepsilon _{xx} = \dfrac{\Delta u_{xx} }{a},\Delta \varepsilon_{yy} = \dfrac{\Delta u_{yy} }{b},\gamma _{xy} = \dfrac{\Delta u_{xy}}{a} + \dfrac{\Delta u_{yx} }{b}$ (10)

 $\label{eq11}\Delta \varepsilon_{xi} = \dfrac{\Delta u_{xi }}{a},\Delta \varepsilon _{yi} = \dfrac{\Delta u_{yi}}{b},\Delta \varepsilon _{zi} = \dfrac{\Delta u_{zi }}{c}$ (11)

 $\label{eq12}\Delta \sigma_{ij} = \lambda \Delta \varepsilon _v \delta _{ij} + 2\mu \Delta \varepsilon _{ij}$ (12)

 $\label{eq13}F_{xi} = \dfrac{bc}{2}\sigma _{xi} ,\ \ F_{yi} = \dfrac{ac}{2}\sigma _{yi},\ \ F_{zi} = \dfrac{ab}{2}\sigma _{zi}$ (13)

 $\label{eq14}\left.\begin{array}{l} F_1^i = - \dfrac{F_{xi} }{2} - \dfrac{mF_{yi} }{2} - \dfrac{F_{zi} }{3}\\[2mm] F_2^i = \dfrac{F_{xi} }{2} - \dfrac{nF_{yi} }{2} - \dfrac{F_{zi} }{3}\\[2mm]F_3^i = \dfrac{F_{yi} }{2} - \dfrac{F_{zi} }{3}\\[2mm]F_4^i= - \dfrac{F_{xi} }{2} - \dfrac{mF_{yi} }{2} + \dfrac{F_{zi} }{3} \\[2mm]F_5^i =\dfrac{F_{xi} }{2} - \dfrac{nF_{yi} }{2} + \dfrac{F_{zi} }{3} \\[2mm] F_6^i = \dfrac{F_{yi} }{2} +\dfrac{F_{zi} }{3} \\ \end{array} \right\}$ (14)

(1) 每个迭代步由单元节点位移增量,求得弹簧变形增量$\Delta u_{ij}(i,j=x,y,z)$;

(2) 由弹簧变形增量$\Delta u_{ij}$及弹簧特征长度求得单元应变增量$\Delta \varepsilon _{ij}$和应力增量$\Delta \sigma _{ij}$;

(3) 由应力全量$\sigma_{ij}$及弹簧特征面积求得弹簧力$F_{ij}$ ($i,j=x,y,z)$,通过坐标变换,插值求得单元节点力;

(4) 根据节点动力平衡方程,求得新的节点位移、速度,如此迭代求解直至收敛.

 图 2 三点弯曲梁模型 Fig.2 Numerical model of three-point beam

 图 3 有限元(a)和弹簧元(b)计算的$y$向位移云图(单位: m) Fig.3 Displacement-contour of three-point beam in y-direction by (a) FEM and(b) SEM (unit: m)
 图 4 两种数值方法计算的节点位移 Fig.4 Node displacement of two numerical methods
 图 5 两种数值方法计算的单元应力 Fig.5 Element stress of two numerical methods

 图 6 转动计算模型 Fig.6 Numerical model of rotation
 $\label{eq15}y = - r\sin (wt)$ (15)

 图 7 转动计算的数值结果和理论解对比 Fig.7 Comparison between numerical result and theory resolution of rotation
2 单元内裂纹的形成与扩展

2.1 单元的破裂

 图 8 块体破裂的莫尔-库伦准则 Fig.8 The Mohr--Coulomb criterion for the block
 $\label{eq16} - \sigma _3 \geqslant \sigma _{\rm t}$ (16)

 $\label{eq17}\sigma _{\rm n} = \dfrac{1}{2}(\sigma _1 + \sigma _3 ) + \dfrac{1}{2}(\sigma _1 -\sigma _3 )\cos 2\theta$ (17)

 $\label{eq18}\tau = \dfrac{1}{2}(\sigma _1 - \sigma _3 )\sin 2\theta$ (18)

 $\label{eq19}s = c + \sigma _{\rm n} \tan \varphi$ (19)

 $\label{eq20} - \sigma _3 < \sigma _{\rm t} ,R \leqslant r$ (20)

 $\label{eq21}\sigma _1 \geqslant \sigma _3 \tan ^2(\pi / 4 + \varphi / 2) + 2c\tan (\pi / 4 + \varphi/ 2)$ (21)

 图 9 三棱柱单元内部破裂示意图 Fig.9 Internal failure modes of triangular prism element

 图 10 节点处裂纹扩展的三种方式 Fig.10 Three modes of crack propagation on node
2.2 裂缝的力学描述

(1) 假定起裂时裂缝面的法向应力为$\sigma _{\rm t}$,即式(16)中的单元抗拉强度. 采用文献[21]建议的线性软化曲线(图11),最大裂缝张开度

 图 11 内聚力模型线性软化曲线[21] Fig.11 Linear softening curve of cohesive model[21]
 $\label{eq22}w_{\rm f} = 2G_{\rm f} / \sigma _{\rm t}$ (22)

 $\label{eq23}\sigma _{\rm n} = \left\{ {\begin{array}{ll} \sigma _{\rm t} (1 - w / w_{\rm f} ) ,\ \ & 0 \leqslant w \leqslant w_{\rm f} \\ 0,& w > w_{\rm f} \\ \end{array}} \right.$ (23)

(2)根据文献[22]假定,定义裂缝面切向刚度按法向张开度进行折减,折减系数为$\beta$,裂缝面黏聚力$c$与裂缝面抗拉强度具有相同软化模式,则有

 $\beta ={{\sigma }_{\text{n}}}/{{\sigma }_{\text{t}}},{c}'=\beta c$ (24)
 $\tau =\left\{ \begin{array}{*{35}{l}} \beta {{k}_{\text{s}}}{{u}_{\text{s}}} & (\left| \tau \right|\le -{{\sigma }_{\text{n}}}\tan \varphi +{c}') \\ \text{sign}({{{\dot{u}}}_{\text{s}}}){{\sigma }_{\text{n}}}\tan \varphi & (\left| \tau \right|>-{{\sigma }_{\text{n}}}\tan \varphi +{c}') \\ \end{array} \right.$ (25)

 图 12 裂缝面两侧的子块体发生分离或滑动 Fig.12 Sub-blocks separate or slide along the crack face

 图 13 子六面体单元的离散方式 Fig.13 Discrete mode of the Sub-hexahedral element
2.3 裂纹扩展方式

 图 14 单元内部的初始裂纹 Fig.14 Initial crack in element

 图 15 裂纹尖端的扩展路径 Fig.15 Propagation path on the crack tip

 $\label{eq26}\bar {\sigma }_{ij} = \dfrac{1}{N}\sum\limits_{k = 1}^N {\sigma _{ij}^k }$ (26)

 图 16 裂纹扩展方式示意 Fig.16 A diagram of crack propagation mode

(1) 根据破裂准则,判断单元123沿路径1---11开裂,同时在裂纹尖端处建立新节点11;

(2) 通过尖端节点11的应力值,按强度准则获得裂纹扩展方向,进而确定开裂路径为11---12;

(3) 在节点12处判断,按此方式,可以依次建立开裂路径12---13和13---14等;

3 数值验证及结果分析

3.1 三点弯曲梁试验模拟

 $\label{eq27}S_{\rm E} = \dfrac{G_{\rm f} }{\sigma _{\rm t} b}$ (27)

 图 17 三点弯曲梁试样的无量纲荷载-位移曲线 Fig.17 Non-dimensional load--deflection curves for the three-point bending
 图 18 不同加载点位移下的试样$x$向位移云图 Fig.18 Displacement-contour of three-point bending specimen in $x$-direction
3.2 拉伸平板中的裂纹扩展

 图 19 几何模型和网格划分 Fig.19 Geometric model and mesh configuration

 图 20 不同切口张开位移时平板$y$向位移云图 Fig.20 Displacement-contours of $y$-direction in different CMOD
 图 21 加载前后网格变化局部放大图 Fig.21 Amplification of local mesh before and after loading
3.3 双切口试样中的裂纹扩展模拟

 图 22 双切口试样几何模型和网格划分 Fig.22 Geometric model and mesh configuration of double-notched specimen

 图 23 双切口试样加载过程的裂纹扩展情况 Fig.23 The crack propagation of the double-notched specimen at differentvertical displacement
 图 24 双切口试样裂纹模式的两种数值方法对比 Fig.24 The crack pattern of double edge-notch specimen; a for continuousdiscontinuous element method; b for Lagrange multiplier method
3.4 单轴压缩试验模拟

 图 25 几何模型和网格划分 Fig.25 Geometric model and mesh configuration

 图 26 单轴压缩试样加载过程的裂纹扩展情况 Fig.26 The crack propagation of the uniaxial compression specimen atdifferent stage

 图 27 试样轴向应力-应变曲线和声发射特征 Fig.27 The axial stress--strain curve and AE rate feature of the uniaxial compression specimen

 图 28 单轴压缩试样加载过程中数值网格变化图 Fig.28 The change of mesh during loading process of uniaxial compressionspecimen
4 结 论

(1) 将块体单元离散为具有物理意义的弹簧系统,在局部坐标系下求解单元的变形和应力,连续问题计算结果与有限元一致,引入最大拉应力与莫尔-库伦复合准则确定单元破坏状态和破裂方向,采用局部切割块体的方式实现单元内部开裂,形成初始裂纹,进而引入内聚力模型描述裂缝特征,扩展路径可穿过单元内部和边界,可以显示模拟裂纹的形成和扩展过程.

(2) 通过单切口平板拉伸、三点弯曲梁、双切口试样等典型岩石裂纹扩展试验对所提方法进行验证分析,结果表明: 该方法可较好地解决拉伸、拉剪以及压剪复合等应力状态下的裂纹形成和扩展问题,裂纹扩展模式能够与试验以及已有数值方法吻合.

(3) 通过单轴压缩试验,对岩石的破坏过程进行分析,数值模拟结果获得了裂纹形成、声发射与应力-应变曲线各阶段之间的对应关系: 弹性阶段后微裂纹开始萌生,声发射速率增加,表现为曲线的非线性段; 峰值点后,局部带不断演化,开始出现宏观裂缝,声发射速率达到最大值,应力-应变曲线表现为软化段,数值模拟揭示的破坏机理与试验一致,说明所提方法模拟岩石裂纹扩展以及破坏问题的有效性.

(4) 本文的连续-非连续单元法是在三棱柱单元上实现,裂纹扩展分析均局限于拟三维情形,后续工作中,可进一步对真三维状态进行研究.

5 讨 论

 [1] 张楚汉. 论岩石、混凝土离散-接触-断裂分析. 岩石力学与工程学报, 2008, 27(2): 217-235 (Zhang Chuhan. Discrete-contact-fracture analysis of rock and concrete. Chinese Journal of Rock Mechanics and Engineering, 2008, 27(2): 217-235 (in Chinese)) [2] Swenson DV, Ingraffea AR. Modeling mixed-mode dynamic crack propagation using finite elements: Theory and applications. Computational Mechanics, 1988, 3: 381-397 [3] Belytschko T, Black T. Elastic crack growth in finite elements with minimal remeshing. International Journal for Numerical Methods in Engineering, 1999, 45(5): 601-620 [4] Ke TC. Simulated testing of two dimensional heterogeneous and discontinuous rock masses using discontinuous deformation analysis. [PhD Thesis]. Berkeley: University of California, 1993 [5] Shi GH. Manifold method. In: Salami MR and Banks D eds. Proceedings of the First International Forum on Discontinuous Deformation Analysis (DDA) and Simulation of Discontinuous Media. Berkeley, California, 1996. 52-204 [6] Shi GH. Numerical manifold method. In: Li JC, Wang CY, Sheng J eds. Proceedings of the First International Conference on Analysis of Discontinuous Deformation (ICADD-1). Chungli, Taiwan, China, 1995. 187-222 [7] 刘丰, 郑宏, 李春光. 基于NMM的EFG方法及其裂纹扩展模拟. 力学学报, 2014, 46(4): 582-590 (Liu Feng, Zheng Hong, Li Chunguang. The NMM-based EFG method and simulation of crack propagation. Chinese Journal of Theoretical and Applied Mechanics, 2014, 46(4): 582-590 (in Chinese)) [8] 杨永涛, 徐栋栋, 郑宏. 动载下裂纹应力强度因子计算的数值流形元法. 力学学报, 2014, 46(5): 730-738 (Yang Yongtao, Xu Dongdong, Zheng Hong. Evaluation on stress intensity factor of crack under dynamic load using numerical manifold method. Chinese Journal of Theoretical and Applied Mechanics, 2014, 46(5): 730-738 (in Chinese)) [9] Munjiza A, Owen DRJ, Bicanic N. A combined finite-discrete element method in transient dynamics of fracturing solids. Engineering Computations, 1995, 12(2): 145-174 [10] Munjiza A. The Combined Finite-Discrete Element Method. New York: John Wiley and Sons, 2004: 277-290 [11] Owen DRJ, Feng YT, De Souza Neto EA, et al. The modeling of multi-fracturing solids and particulate media. International Journal for Numerical Methods in Engineering, 2004, 60(efeq1): 317-339 [12] Li SH, Zhang YN. Feng C. A spring system equivalent to continuum model. In: The Proc. of The Fifth International Conference on Discrete Element Methods. London, 2010. 75-85 [13] 冯春, 李世海, 姚再兴. 基于连续介质力学的块体单元离散弹簧法研究. 岩石力学与工程学报, 2010, 29(增刊1): 2690-2704 (Feng Chun, Li Shihai, Yao Zaixing. Study of block-discrete-spring method based on continuum mechanics. Chinese Journal of Rock Mechanics and Engineering, 2010, 29(Supp.1): 2690-2704 (in Chinese)) [14] 王杰, 李世海, 周东, 等. 模拟岩石破裂过程的块体单元离散弹簧模型. 岩土力学, 2013, 34(8): 2355-2362 (Wang Jie, Li Shihai, Zhou Dong, et al. A block-discrete-spring model to simulate failure process of rock. Rock and Soil Mechanics, 2013, 34(8): 2355-2362 (in Chinese)) [15] 常晓林, 胡超, 马刚, 等. 模拟岩石失效全过程的连续-非连续变形体离散元方法及应用. 岩石力学与工程学报, 2011, 30(10): 2004-2011 (Chang Xiaolin, Hu Chao, Ma Gang, et al. Continuous-discontinuous deformable discrete element method to simulate the whole failure process of rock masses and application. Chinese Journal of Rock Mechanics and Engineering, 2011, 30(10): 2004-2011 (in Chinese)) [16] 马刚, 周创兵, 常晓林, 等. 岩石破坏全过程的连续-离散耦合分析方法. 岩石力学与工程学报, 2011, 30(12): 2444-2455 (Ma Gang, Hu Chao, Chang Xiaolin, et al. Continuous-discontinuous coupling analysis for whole failure process of rock. Chinese Journal of Rock Mechanics and Engineering, 2011, 30(12): 2444-2455 (in Chinese)) [17] 侯艳丽, 周元德, 张楚汉. 用3D离散元实现I/II型拉剪混凝土断裂的模拟. 工程力学, 2007, 24(3): 1-7 (Hou Yanli, Zhou Yuande, Zhang Chuhan. I/II tensile shear mixed mode fracture simulation by 3D discrete element method. Engineering Mechanics, 2007, 24(3): 1-7 (in Chinese)) [18] 张青波, 李世海, 冯春. 四节点矩形弹簧元及其特性研究. 岩土力学, 2012, 33(11): 3497-3502 (Zhang Qingbo, Li Shihai, Feng Chun. Study of four-node rectangular spring element and its properties. Rock and Soil Mechanics, 2012, 33(11): 3497-3502 (in Chinese)) [19] 张青波, 李世海, 冯春, 等. 基于SEM的可变形块体离散元法研究. 岩土力学, 2013, 34(8): 2385-2392 (Zhang Qingbo, Li Shihai, Feng Chun, et al. Study of deformable block discrete element method based on SEM. Rock and Soil Mechanics, 2013, 34(8): 2385-2392 (in Chinese)) [20] 张楚汉, 金峰. 岩石和混凝土离散-接触-断裂分析. 北京: 清华大学出版社, 2008 (Zhang Chuhan, Jin Feng. Discrete-Contact-Fracture Analysis of Rock and Concrete. Beijing: Tsinghua University Press, 2008 (in Chinese)) [21] Hillerborg A, Modeer M, Petersson PE. Analysis of crack formation and crack growth in concrete by means of fracture mechanics and finite elements. Cement and Concrete Research, 1976, 6(6): 773-782 [22] Wells GN, Sluys LJ. A new method for modeling cohesive cracks using finite elements. International Journal for Numerical Methods in Engineering, 2001, 50(12): 2667-2682 [23] Carpinteri A, Colombo G. Numerical analysis of catastrophic softening behavior (snap-back instability). Computer & Structures, 1989, 31(4): 607-636 [24] Ning YJ, An XM, Ma GW. Footwall slope stability analysis with the numerical manifold method. International Journal of Rock Mechanics & Mining Sciences, 2001, 48(6): 964-975 [25] Nooru-Mohamed MB. Mixed-mode fracture of concrete: An experimental approach. [PhD Thesis]. Delft: University of Technology, 1992 [26] Oliver J, Huespe AE, Samaniego E, et al. Continuum approach to the numerical simulation of material failure in concrete. International Journal for Numerical and Analytical Methods in Geomechanics, 2004: 28(7-8): 609-632 [27] Zi G, Rabczuk T, Wall W. Extended meshfree methods without branch enrichment for cohesive cracks. Computational Mechanics, 2007, 40(2): 367-382
SIMULATION OF CRACK PROPAGATION OF ROCK BASED ON SPLITTING ELEMENTS
Wang Jie, Li Sihai, Zhang Qingbo
1. Institute of Mechanics, Chinese Academy of Sciences, Beijing 100190, China;
2. The Third Railway Survey and Design Institute Group Corporarion, Tianjin 300142, China
Fund: The project was supported by the Strategic Priority Research Program (B) of Chinese Academy of Sciences (XDB10030303) and Youth Science Fund of the National Natural Science Foundation of China (11302230, 11302229).
Abstract: In conventional discrete element methods, fracture is judged by criterion of interface and cracks can only propagate along the boundary of elements. However, criterion of interface can only be used rationally on the condition that macro or micro fractures exist in physical problems. The path and direction of crack will be limited severely by the initial mesh when crack propagates along the boundary. Given these two limitations, a continuous-discontinuous element method is proposed and applied to simulate the progressing cracking problem of rocks. Specifically, criterion is applied on element and intra-element fracture will form. In continuous calculation, element is denoted by a discrete spring system which has specific physical meaning and its deformation and stress are calculated by the characteristic length and area of springs in local coordinate system. The continuous calculation results demonstrate a satisfactory agreement with the traditional finite element method. By updating spring information and local coordinate system, large displacement and rotation of elements can be calculated directly. In addition, Mohr-Coulomb criterion is implemented into the new model to specify the failure state and fracture direction, and intact element will be divided into two elements by means of cutting block. In this way, fracture may be inserted along the boundary of elements or within intact element. A cohesive zone model is employed to simulate the fracture and the elements on two sides of the crack are set to two different nodes at the same time, causing the displacement to be discontinuous. Finally, from numerical results of several intense examples with crack propagation, this method can satisfactorily simulate the progressing cracking problems under tensile, compressive and shear conditions, and its rationality is approved. The continuous-discontinuous element method has been shown to be insensitive to quality of mesh and thus has the potential to simulate crack initiation and propagation.
Key words: continuous-discontinuous    element fracture    crack propagation    rock failure