水泥圆柱试件断裂破坏的数值模拟

罗新荣， 李梦坤， 李亚伟， 杨 威

(中国矿业大学 安全工程学院， 江苏 徐州 221116)



水泥圆柱试件断裂破坏的数值模拟

罗新荣， 李梦坤， 李亚伟， 杨 威

(中国矿业大学 安全工程学院， 江苏 徐州 221116)

煤岩瓦斯动力灾害大多表现为瞬间断裂破坏，其破坏机理及过程还不十分清楚。以不同配比的水泥圆柱试件为研究对象，借助FLAC3D软件建立三维模型，对实验现象进行模拟和机理分析。结果表明：模型试件的轴向拉应力达到抗拉强度时，试件中心出现局部拉伸破坏区域，此时的围压与实验加载破断围压完全一致；继续增大围压，最大拉应力达到抗拉强度后不再增加。误差性分析证明，当σc=10σt时，模拟的临界围压与实验断裂破坏吻合最好。该结果对于研究三维岩石拉伸断裂瞬间破坏和煤岩瓦斯动力灾害瞬间断裂破坏机理有积极的意义。

水泥圆柱； 断裂破坏； 围压； 拉应力； 拉伸破坏； 误差分析

0 引 言

瞬间断裂破坏是岩石破坏的基本形式之一。在煤炭资源的开采过程中，冲击地压、煤与瓦斯突出等动力灾害事故经常发生，威胁着矿井的安全生产。因此，研究岩石瞬间断裂破坏具有重要的理论意义和工程应用价值。

在岩石力学领域，岩石断裂实验研究主要集中在压应力条件下的破坏特性方面，包括单轴、双轴和三轴压缩实验。很多学者深入研究了压应力条件下岩石的裂纹扩展、声发射特征以及围压对单轴抗压强度的影响等，并取得了丰富的成果[1-4]。拉应力作用下岩石断裂破坏特性的研究文献相对较少，分析主要有两方面的原因：一是实际岩土工程中岩石大多是处于受压应力状态；二是实验条件限制。以往的力学实验机和测试装置无法对拉应力作用下的裂缝扩展和变形破坏规律进行定量研究，测得的抗拉强度差别比较大。目前，岩石的抗拉强度通常采用巴西圆盘劈裂实验间接测得，实验过程中应力状态复杂，不能真实反映岩石在拉应力条件下的力学性质[5]。部分学者改进实验装置，运用直接拉伸方法对多种岩石以及煤的抗拉力学特性展开研究，不同装置下实验结果差别依然很大[6-12]。在岩石断裂数值模拟方面，很多学者也作了大量研究，梁正召等[13-16]采用RFPA3D软件对岩石断裂过程及裂纹扩展进行了模拟分析；赵延林等[17-18]采用Ansys软件实现了翼形裂纹尖端关键点随翼形裂纹扩展的自动生成；周斌等[19]采用FLAC3D模拟了层状岩体的三轴压缩实验；张春会等[20]采用数值模拟的方法对岩石应变软化及渗透率演化模型作了细致研究。目前，岩石断裂破坏的理论研究仍存在许多不足，导致实验与模拟结果无法匹配，不能准确反映岩石断裂破坏的力学特性。

基于此，对于岩石的断裂破坏，笔者试图改变已有的直接拉伸或者压缩的实验方法，以不同围压作用下的岩石瞬间断裂破坏实验为基础，采用FLAC3D软件对实验结果进行数值模拟，从细观层次上分析岩石瞬间断裂的宏观破坏行为。

1 围压加载实验及数值模型

1.1 试件制备与力学参数测定

实验中采用不同配比的水泥和沙子分别制作了A、B两组混凝土试件，其中，A组试件水泥、沙子质量比为1∶3.5，B组为1∶6。根据实验设备要求以及力学参数测试标准，制作两组共计八个混凝土试件。第一组为直径50 mm、高度100 mm、高径比为2的圆柱体试件六个。第二组制作两个直径110 mm、高150 mm的圆柱体试件，对圆柱体试件端面进行了打磨加工，自然条件下风干20 d后进行后续测试。

1.2 试件力学参数测试结果

利用加工好的六个直径50 mm、高度100 mm的圆柱体试件，测定不同配比混凝土试件的单轴抗压强度、泊松比及弹性模量；对直径110 mm、高150 mm的第二组试件进行围压实验。利用围压加载仪分别对两组不同硬度的试块加载围压，测量加载过程的波速变化及断裂破坏。加载过程的全应力-应变曲线，如图1所示。

a A组试件

b B组试件

结合直线段全应力、应变，分析A、B两组全应力-应变曲线，得到两组试件基本力学参数，见表1。

从数据结果来看，混凝土试件具有一定的离散性，但随着沙子比例的增高，试件的抗压强度以及弹性模量都有明显降低，泊松比随试件配比变化较小，两组试件的泊松比都集中在0.17～0.20之间。

表1 两组试件力学参数

Table 1 Table of mechanical parameters of two sets of test pieces

水泥与沙子质量比试样σc/MPaE/MPaμ实验值均值实验值均值实验值均值1∶3.5111.36213.50310.2311.700.981.050.740.920.1750.1860.1740.1781∶614.3225.3934.554.750.340.530.490.450.1840.1960.2150.198

注：以三个试样参数均值作为模拟参数。

1.3 数值模型与参数确定

采用FLAC3D5.0软件分别模拟两组配比的水泥柱试件在不同围压下的拉应力及塑性区域变化。试件尺寸为φ110 mm×150 mm。由于轴对称，所以数值模型选择四分之一圆柱体，网格采用渐进式分布，在边缘加密，网格划分22 500个，如图2所示。为了消除试件两端的端头效应，围压施加时呈梯形分布。

a 围压分布 b 模型网格划分

岩体的破坏在宏观上表现为所受的应力达到极限强度而发生破坏[21]。大量实验证明，在低围压下，岩石发生破坏时的剪切应力与主应力为线性关系，当围压增大后，呈现出明显的非线性关系，是各种破坏条件下莫尔应力圆的弧形包络线[22]。针对岩性较脆的岩石，莫尔强度包络线近似于二次抛物线，如图3所示。表达式为

τ2=n(σ+σt),

(1)

式中：σt——岩石的单轴抗拉强度；

n——待定系数。

图3 二次抛物线形强度包络线

利用图3中的关系，有

(2)

联合式(1)和式(2)，得到二次抛物线形包络线的主应力表达式:

(σ1-σ3)2=2n(σ1+σ3)+4nσt-n2。

(3)

单轴压缩条件下，有σ3=0，σ1=σt，则式(3)变为

n2-2(σc+2σt)n+σc2=0。

(4)

由式(4)可解得

(5)

联立式(2)、(3)、(5)得出A组、B组试件的参数黏聚力c和抗拉强度σt，如表2所示。

表2 部分模拟参数

2 数值模拟与实验结果

2.1 数值模拟结果

数值模拟过程中，取抗压强度是抗拉强度的8～12倍(见表2)，对试件进行围压加载数值模拟。数值模拟结果表明，抗压强度是抗拉强度的10倍时模拟结果最符合实验结果。A、B组试件的加载围压分别为3.0～4.0、1.5～2.5 MPa，步进0.1 MPa，分别模拟10组结果；选择典型的三组数值模拟结果，对比分析圆柱中心侧轴向应力随围压的变化，如图4所示。从图4中可以看出，模型试件在围压的作用下，轴向方向有拉应力存在，拉应力大小呈抛物线形，A组试件在围压加载到3.6 MPa时，轴向拉应力恰好达到抗拉强度；B组试件的轴向拉应力达到抗拉强度时，围压临界值为2.2 MPa，此时模型试件的临界围压恰好与实验破断围压一致。随着围压的增大，模型试件中心轴向最大拉应力达到抗拉强度后，拉应力不再增大，并由中心方向向两端扩展，与实验不符合。实验结果表现为轴向应力(σx)达到抗拉强度时，试件瞬间拉伸破断。

a A组试件

b B组试件

岩石在受外力挤压时多表现为断裂破坏。通常情况下，在FLAC3D数值模拟中塑性云图反应了不同围压下模型圆柱试件单元的塑性损伤破坏情况，由于圆柱试件具有对称性，文中取模型试件圆柱的1/4来观察破坏单元的塑性情况。塑性破坏云图如图5所示。从图5中可以看出，A、B两组试件在围压分别为3.3、2.0 MPa的情况下，试件本身没有发现拉伸破坏区域，当围压分别增加到临界围压3.6、2.2 MPa时，模型试件中心出现局部拉伸破坏区域；随着围压的增大，破坏区域在向试件两端方向扩展的同时，也慢慢由试件的中心向侧面扩展，破坏区域形状呈“椭球形”，而没有瞬时拉伸破断。原因是现有的数值模拟软件，均无法统一岩石应力-应变的本构方程和瞬间拉伸断裂现象。因此，数值模拟过程中轴向拉应力达到临界值以后，模型拉伸破坏并没有瞬时断开，这一点与实验结果不符。

a A组试件

b B组试件

Fig. 5 Failure of plastic damage under different confining pressures

2.2 实验结果验证

在实验过程中，由于横波首波位置不明显，纵波容易观察，采用纵波来分析实验加载过程。随着围压的增大，纵波波速变化较大，在围压增大到一定程度时，试件发生断裂，波速降为0。加载后试件断裂如图6所示，围压-波速的变化如图7所示。

a A组试件 b B组试件

由图7可见，实验过程中波速随围压(p)的变化而变化，A组试件加压至3.5 MPa的过程中，波速随围压的增大而基本保持不变;围压由3.5 MPa加压至4.0 MPa时，波速迅速下降，有明显的响声，试件发生瞬时断裂，断裂处在试件本身的1/2位置。B组试件围压低于2.1 MPa时，波速随围压增大保持不变，围压由2.1 MPa加压至2.5 MPa的过程中，有明显响声，试件开始断裂，断裂后波速下降为0，在1/3处出现断裂痕迹。两组试件断裂方式多表现为由轴心向外的水平断裂，即为拉伸破坏。模拟的试件临界围压分别为3.6和2.2 MPa，恰好在实验结果范围之内，可以验证该模拟具有一定的指导意义。

图7 波速-围压变化

3 分析讨论

3.1 试件微观裂纹扩展的基本原理

岩石的破裂是一个演化问题，针对于脆性岩石，破坏机理表现出拉伸破坏。Griffith用能量平衡的观点确定了脆性材料裂纹扩展准则，微观上岩石内部存在错综复杂的细小裂纹，在外应力作用下，材料中的微小裂纹尖端存在拉应力集中且达到极限应力状态致使断裂开始扩展[23]。图8给出了二维构形下脆性介质水泥圆柱的微观开裂过程。图8a表示该二维介质在无应力状态下的长度为L；若对介质施加围压，介质在竖直方向上产生拉应力，位移增量为Δ，如图8b所示；图8c表示该介质中央切开一条长度为2a的水平走向中心裂纹(a<

a b

c d

Fig. 8 Schematic diagram of cracking process of brittle

medium

根据热力学第一定律有如下能量方程：

U+T+Γ=W+Q，

(6)

式中：U——储存在介质中的热力学能；

T——动能；

Γ——表面能；

W——外载在变形过程中所做的功；

Q——外界提供给物体的热量。

假设断裂过程是绝热过程和准静态过程(即裂纹扩展速率远小于应力波速)，则Q=0，T=0，式(6)为

U+Γ=W。

(7)

对于固定位移加载下的断裂，有dW=0。若介质为弹性，则U等同于弹性应变能Ue，对式(7)进行求导，得

(8)

由线弹性力学的叠加原理[24]可知：图8d中在形成2a中心裂纹时所释放的能量等价于裂纹表面作用均匀压应力σc而远场为自由的构形中所蕴藏的能量。则易于求出

Γ=4adγ，

(9)

式中：d——板厚；

γ——单位表面自由能。

断裂过程所释放的应变能(用-Ue表示)可由虚功原理得出

(10)

(11)

式中：κ——与泊松比有关的无量纲参数，κ=(3-v)/(1+v)。

将式(11)代入式(10)积分得出

(12)

(13)

E——非破坏材料的弹性模量；

μ——泊松比。

因此，试件微观断裂依据为

(14)

微观角度上试件在围压的加载下，当轴向拉应力σt大于临界拉应力σcr时，试件内部有裂纹出现并开始扩展[25]。针对脆性水泥柱试件，若内部裂纹尖端的应力场超过材料的抗拉强度，裂纹迅速扩展和贯通，最终导致脆性岩石的拉断裂破坏。

3.2 数值模型误差性分析

以B组试件为研究对象，在σc=8σt，σc=10σt，σc=12σt三种情况下进行数值模型误差性分析，如图9所示。从图9中得知：当σc=12σt(即σt=0.396 MPa)的情况下，B组试件在围压加载到1.9 MPa时，轴向拉应力达到0.396 MPa，当围压加大到2.2 MPa时，最大拉应力仍停滞在抗拉强度上，表明此种情况下B组试件断裂的临界围压为1.9 MPa。由此分析，当σc=8σt(即σt=0.594 MPa)的情况下，B组试件的断裂临界围压为2.8 MPa。实验结果表明，在围压为2.1～2.5 MPa时，B组试件开始断裂，而上述两种情况下的临界围压都在实验结果范围之外，因此，取抗压强度是抗拉强度的10倍进行数值模拟更能符合实验结果。

a σt=0.396 MPa

b σt=0.594 MPa

4 结 论

(1)针对两组不同配比的圆柱形水泥试件进行围压加载，A、B组试件分别在围压加载到3.5、2.2 MPa时，纵波波速开始迅速下降，试件发生瞬时断裂。A组试件的断裂处在其本身的1/2位置，B组试件在1/3处发生断裂，断裂形式为瞬时拉伸破坏。

(2)试件数值模拟结果表明，随着围压的增大，试件轴向拉应力达到抗拉强度后收敛，此时试件的临界围压分别为3.6、2.2 MPa。模拟结果与实验加载破断的围压基本一致。对数值模型进行误差性分析，证明σc=10σt时，临界围压最符合实验结果。

(3)实验中试件是瞬时拉伸破坏，而数值模拟结果是，随着围压的增大，拉伸破坏区域缓慢向四周扩展，未表现出拉伸破坏的瞬时性，这与软件的本构模型有关。对于脆性岩石的拉伸破坏，计算软件系统还需要进一步研究完善。

[1] 李守巨, 李 德, 武 力, 等. 非均质岩石单轴压缩试验破坏过程细观模拟及分形特性[J]. 煤炭学报, 2014, 39(5): 849-854.[2] 张社荣, 孙 博, 王 超, 等. 双轴压缩试验下岩石裂纹扩展的离散元分析[J]. 岩石力学与工程学报, 2013(S2): 3083-3091.[3] 杨永杰, 王德超, 郭明福, 等. 基于三轴压缩声发射试验的岩石损伤特征研究[J]. 岩石力学与工程学报, 2014, 33(1): 98-104.[4] 杨永杰, 宋 扬, 陈绍杰. 三轴压缩煤岩强度及变形特征的试验研究[J]. 煤炭学报, 2006, 31(2): 150-153.

[5] 中华人民共和国国家标准编写组. GB/T50266—99 工程岩体试验方法标准[S]. 北京： 中国计划出版社， 1999.

[6] 周小平, 张永兴, 哈秋舲, 等. 单轴拉伸条件下细观非均匀性岩石变形局部化分析及其应力-应变全过程研究[J]. 岩石力学与工程学报, 2004, 23(1): 1-6.

[7] 代高飞, 夏才初, 晏 成. 龙滩工程岩石试件在拉伸条件下的变形特性试验研究[J]. 岩石力学与工程学报, 2005, 24(3): 384-388.

[8] 余贤斌, 谢 强, 李心一, 等. 直接拉伸、劈裂及单轴压缩试验下岩石的声发射特性[J]. 岩石力学与工程学报, 2007, 26(1): 137-142.

[9] 刘世奇, 李海波, 李俊如, 等. 有侧向压力的岩石材料动态直接拉伸试验装置[J]. 岩土力学, 2007, 28(11): 2365-2368.

[10] 李地元, 李夕兵, CHARLIE LIC. 2种岩石直接拉压作用下的力学性能试验研究[J]. 岩石力学与工程学报, 2010, 29(3): 624-632.

[11] 张泽天, 刘建锋, 王 璐, 等. 煤的直接拉伸力学特性及声发射特征试验研究[J]. 煤炭学报, 2013, 38(6): 960-965.

[12] 张绪涛, 张强勇, 袁圣渤, 等. 岩石轴向直接拉伸试验装置的研制及应用[J]. 岩石力学与工程学报, 2014, 33(12): 2517-2523.[13] 梁正召, 唐春安, 张永彬, 等. 岩石直接拉伸破坏过程及其分形特征的三维数值模拟研究[J]. 岩石力学与工程学报, 2008, 27(7): 1403-1410.

[14] 梁正召, 李连崇, 唐世斌, 等. 岩石三维表面裂纹扩展机理数值模拟研究[J]. 岩土工程学报, 2011, 33(10): 1616-1622.[15] 梁正召, 唐春安, 张永彬, 等. 岩石三维破裂过程的数值模拟研究[J]. 岩土力学与工程学报, 2006, 25(5): 932-936.

[16] 梁正召, 肖东坤, 李聪聪, 等. 断续节理岩体强度与破坏特征的数值模拟研究[J]. 岩土工程学报, 2014, 36(11): 2087-2095.[17] 赵延林, 曹 平, 万 文, 等. 岩石裂纹亚临界扩展实验与压剪流变断裂模型[J]. 中南大学学报: 自然科学版, 2014, 45(1): 277-286.

[18] 赵延林, 万 文, 王卫军, 等. 类岩石材料有序多裂纹体单轴压缩破断试验与翼形断裂数值模拟[J]. 岩土工程学报, 2013, 35(11): 2098-2109.

[19] 周 斌, 张可能, 柳群义. 考虑应变软化特征的层状岩体三轴压缩数值试验分析[J]. 中南大学学报: 自然科学版, 2009, 40(3): 829-832.

[20] 张春会, 赵全胜, 王来贵, 等. 三轴压缩岩石应变软化及渗透率演化的试验和数值模拟[J]. 煤炭学报, 2015， 40(8): 1775-1782.

[21] 张 荣. 保护层开采上覆煤岩体走向边界卸压特征研究[D]. 北京： 中国矿业大学, 2016: 11-40.

[22] 蔡美峰, 何满潮, 刘东燕. 岩石力学与工程[M]. 北京： 科学出版社, 2002: 219-227.

[23] 杨 卫. 宏微观断裂力学[M]. 北京： 国防工业出版社, 1995: 15-21.

[24] 李世愚, 和泰名, 尹祥础. 断裂力学导论[M]. 北京： 中国科学技术大学出版社, 2010: 42-50.

[25] 王自强, 陈少华. 高等断裂力学[M]. 北京： 科学出版社, 2009: 14-30.

(编辑 晁晓筠 校对 荀海鑫)

Numerical simulation and analysis of mechanism behind fracture and failure of cement column

LuoXinrong,LiMengkun,LiYawei,YangWei

(School of Safety Engineering, China University of Mining & Technology, Xuzhou 221116, China)

The dynamic disasters of coal and gas occur in the typical form of instant fractures governed by an elusive mechanism and process. This paper introduces the simulation analysis and mechanism analysis of experimental phenomena, using different proportions of cement column specimens as the research object and drawing on a three-dimensional model developed by the aid of FLAC3Dsoftware. Numerical simulation analysis shows that, when the axial tensile stress of model specimen reaches the tensile strength, the local tensile failure zone occurs at the center of the specimen, during which time the confining pressure is completely consistent with that of the experiment, meaning that a continued increase in the confining pressure produces no more increase in the maximum tensile stress after the stress is up to tensile strength; and the error analysis of numerical model verifies a best agreement between the critical confining pressure and the experimental fracture, whenσc=10σt. The research results are of positive significance for a better insight into the mechanism underlying both the three-dimensional rock tensile fracture instant destruction and gas dynamic disaster instant fracture.

cement column; fracture failure； confining pressure； tensile stress； tensile failure； error analysis

2016-08-23

国家自然科学基金项目(U1361102)

罗新荣(1957-)，男，江西省樟树人，教授，博士生导师，研究方向：矿井瓦斯防治理论与技术，E-mail：aq204@cumt.edu.cn。

10.3969/j.issn.2095-7262.2016.05.006

TU452

2095-7262(2016)05-0496-06

A