脆性岩石微裂纹压密段本构模型研究*

乔　彤，章　光，胡少华,3，赵顺利

(1. 武汉理工大学 资源与环境工程学院，湖北 武汉 430070；2. 安全预警与应急联动技术 湖北省协同创新中心，湖北 武汉 430070；3. 武汉大学 水资源与水电工程科学国家重点试 验室，湖北 武汉 430072；4. 黄河勘测规划设计有限公司，河南 郑州 450003)

0　引言

当前经济社会的迅速发展使得我国岩石力学与工程开始向深部进军，未来深地下空间利用、深部资源开采将成为常态。例如，我国矿产资源开采已逐步进入1 000～2 000 m深水平、西南地区大型水利水电工程地下硐室埋深普遍在高地应力地区、地下核废料处置库预选处置库大多埋深在500 m以上，伴随而来的是岩体面临高地应力、高地温、高渗透水压的复杂深地质赋存环境。因此，准确掌握岩石在复杂环境作用下的变形与强度特性是深部岩体工程安全迫切需要解决的基础性问题之一。

深部岩石力学工程围岩以脆性岩石为主，最常见的有花岗岩、石英闪长岩等，因此脆性岩石的应力-应变关系引起了广大学者的浓厚兴趣，是近二十年来深部岩石力学研究最活跃的课题。E. Eberhardt等[1]在总结Brace和Bieniawski[2]关于脆性岩石应力-应变关系研究成果基础上，从裂纹发展的角度提出图1所示的应力-应变关系曲线，其定性的描述了不同阶段裂纹行为状态与脆性岩石应力-应变之间的关系，主要包括微裂纹压密段、线弹性阶段、微裂纹稳定扩展阶段、微裂纹不稳定扩展阶段以及峰后阶段，从机理角度阐释岩石的本构关系曲线的特征，但是限于岩石裂纹行为的复杂性，并未建立对应的本构模型。

为了使脆性岩石的本构模型更加贴近实际，一些基于岩石断裂力学与裂纹行为状态的本构模型先后提出，Y.L. Lu等[3]基于传统的连续介质模拟方法和流体力学耦合理论，提出了双尺度模型，用以解决岩石断裂耦合过程中水力耦合问题；周小平等[4]则根据岩石不同阶段的应变特征，将围压状态分为低压、中压、高压3种状态，分别建立基于裂纹扩展贯通机制的本构模型；Hui-Hai Liu等[5]受自然应变和工程应变的启发，将岩石分为裂隙为主的软部和基质为主的硬部，得到两者串联的本构模型；袁小清等[6]认为脆性岩石的原始状态的损伤由于宏微观裂纹的存在，其值并不为零，从而将初始损伤的概念加入了本构模型的研究中；朱其志等[7]则进一步发展了热动力理论，将脆性岩石进行一定的均质化处理，并发展出了一种适用于脆性岩石非线性力学行为模拟的新型本构模型；区别于大多数的研究，刘红岩等[8]的研究重点主要是脆性岩石的动态损伤过程，更进一步的推广了断裂力学的应用范围。

虽然脆性岩石本构模型的研究引入断裂力学和裂纹行为状态的研究，但针对岩石微裂纹压密段的解析方程更多的停留在经验公式[9]。陈益峰等[10]提出的THM本构模型中考虑了微裂纹压密段，但仅仅是单裂隙结构的简单推广；荣冠等[11]基于柔度变形法推导出同时适合耦合和非耦合节理法向变形的g-λ模型；彭俊等[12]根据经验提出微裂纹压密段的本构模型，但是参数的物理意义指向不明，限制了进一步推广。为此，本文引入自然应变的概念，从单裂隙在法向荷载作用下的应力-应变关系出发，推导并建立了物理意义明确的微裂纹压密段的本构模型，然后将其和线弹性阶段的本构模型进行了统一。在此基础上，通过花岗岩的热-力耦合试验对本文提出的脆性岩石微裂纹压密段本构模型的合理与适用性进行了验证。

图1　岩石应力-应变关系曲线Fig.1　Stress strain relationship of rock

1　脆性岩石微裂纹压密段本构模型研究

1.1　单裂隙法向应变

脆性岩石在应力加载初期，随着应力增加，应力-应变曲线呈下凹型，表现出明显的微裂纹压密行为，为阐明微裂纹压密阶段脆性岩石的力学特征，本文首先研究单裂隙在法向荷载作用下的变形过程。

图2　岩石单裂隙Fig.2　Rock structure with single fracture

对于图2所示岩石单裂隙，在上下表面受均布法向应力σ作用，并假定岩石为均质材料，根据Hook定律可得:

dσ=Kdεv,t

(1)

式中：K为法向刚度；εv,t为法向自然应变或真实应变，其定义式如下：

(2)

式中：L为岩石结构在当前应力状态下的长度，在式(1)、(2)中，文中以压缩为正。

目前多数研究成果中，在应用Hook定律时多采用如下的形式[13]：

(3)

式中：L0为没有应力作用时岩石结构的初始法向长度；εv,e为法向工程应变，也是目前最多采用的表达方式。当针对法向工程应变εv,e与法向应力σ采用公式(3)所示的Hook定律时，其初始条件满足σ=0，L=L0。据此可得：

(4)

同理，当自然应变εv,t采用公式(2)中的Hook定律时，可得法向长度L与法向应力σ的函数关系如式(5):

(5)

基于公式(4)和(5)，定义σ/K为应力刚度比，分别绘制应力刚度比σ/K与法向长度L关系图(见图3)。通过图中不难发现，当应力刚度比较小时(σ/K<0.1)，公式(4)和公式(5)基本上是等价的。对于脆性完整的岩石而言，其刚度量级一般可以达到10～50 GPa，对应的单轴抗压强度一般为100～300 MPa，应力刚度比均满足σ/K<0.1的条件，因此选用形式更简单的公式(4)可以较好的表达其变形特征。随着σ/K的增大，公式(4)与公式(5)出现了明显的偏离。Freed[14]在回顾前人自然应变研究基础上，提出自然应变较工程应变能更精确的描述材料的变形特性，因此对于法向刚度较小的材料(流体、断层带、大尺度裂隙、软土等)不适宜采用公式(4)。

在法向应力作用下，岩块的变形特征可用公式(4)描述，但是裂隙结构的变形特征需要通过公式(5)进行精确描述。可得如下公式:

L0=Lc0+Lr0

(6)

dL=dLc+dLr

(7)

式中：L0为岩石结构的整体初始法向长度；Lr0为岩块初始法向长度；Lc0为裂隙初始法向长度；L为某一应力状态下的岩石结构的整体初始法向长度；Lr为某一应力状态下的岩块初始法向长度；Lc为某一应力状态下的裂隙初始法向长度。

将公式(4)、(5)代入(6)、(7)中，可得整体法向长度L的表达式：

(8)

图3　应力刚度比与法向长度Fig.3　Stress stiffness ratio versus normal length

通过式(8)可进一步求出岩石整体的法向变形量ΔL和采用常用的工程应变定义的εe：

(9)

(10)

式中：Kr和Kc分别为岩块和裂隙的法向刚度；nr和nc为工程应变中岩块工程应变和裂隙工程应变所占的比例因子，其定义式如式(11)所示：

(11)

通过式(9)、(10)可以发现，传统意义上的变形量与应变实际上由2部分组成：刚度较大的岩块部分，其应力-应变关系可以用线性关系较好的表达；刚度较小的裂隙部分，其应力-应变关系可以用负指数模型表达。

1.2　标准试样压密段应力-应变

室内进行岩石力学试验，多采用圆柱形试样，包含有多种裂隙等低刚度结构和岩石颗粒、孔隙等高刚度结构，如图4中的岩石试样模型所示。对于常规三轴试验中的岩石，其应力状态可用图5来表示，真实状态下的应力状态可分解为静水压力状态和偏应力状态。在具体的室内三轴压缩试验过程中，静水压力状态被认为是应变的起始测量点，其岩石试样产生的应变主要为偏应力状态下产生的应变。

图4　岩石试样结构Fig.4　Rock sample structure diagram

图5　应力状态分解Fig.5　Stress state decomposition

任取岩石中任意角度的单裂隙为分析对象，如图6所示。裂纹上下表面面积相等，为Si；裂纹长度为2ci；裂纹面方向与最大主应力方向夹角为βi。将裂纹投影到与最大主应力垂直的方向上和与最大主应力平行方向上，因此裂纹可用图7简化的表示。因裂纹为低刚度结构，因此应用公式(5)可得与最大主应力垂直的方向上投影后裂纹法向方向变形量如下：

(12)

式中：L0i为投影裂隙的初始长度。通过式(12)可进一步得出，图6中裂隙在最大主应力方向上对试样整体应变的贡献。同时将裂隙的变形在岩样的尺度上进行均布化处理，可得公式(13)：

图6　岩石结构中单裂隙分析对象Fig.6　Analysis object of single fracture in rock structure

图7　岩石结构原生裂纹简化模型Fig.7　Simplified model of rock structure primary crack

(13)

式中：L0为试样的初始长度；S0为试样的初始截面积；V0为试样的初始体积；Vi为裂纹投影后的体积；εi为均布化后的裂隙对试样整体应变贡献量。

通过式(13)可以发现，对于具有相同刚度Kj的裂纹系，其对试样的整体的应变可用式(14)表示：

(14)

式中：εj为裂纹系对试样整体应变的贡献量；N为具有相同刚度Kj的裂纹系中裂纹数目；Vj为裂纹系的中裂纹投影体积之和。进一步定义比例因子nj=Vj/V0，可以发现nj是一个与裂纹体积密切相关的量。此时式(14)可用式(15)表示：

(15)

最大主应力平行方向上的投影裂隙难以被压缩，其性质类似于岩石中孔洞。因此可认为是高刚度结构。

岩石试样中同样存在有其他刚度的裂纹系以及非裂纹低刚度结构。其在最大主应力方向上的应变均可以用式(15)表示。对于高刚度的岩石成分，如岩石颗粒，其应变可用公式(5)求得，岩石试样中高刚度成分在最大主应力方向上的应变表达式如式(16)：

(16)

设岩石中共有a种高刚度成分，b种低刚度成分。则岩石的应变ε可用式(17)表示：

(17)

(18)

考虑式(17)中的高刚度项和低刚度项。高刚度项的应变与偏应力呈线性关系，主要与岩石的线弹性阶段有关；低刚度项的应变主要表现为随着偏应力增加，其迅速减小。同样根据式(18)可以得出岩石试样的弹性模量E如下：

(19)

通过式(19)可以发现，低刚度部分对弹性模量的影响已经很小，并且对于一般脆性岩石而言，低刚度成分远小于高刚度成分的含量，对于岩石的应力-应变曲线而言，其线弹性阶段的应力-应变曲线的弹性模量可近似用下式表示：

(20)

通过式(20)不难看出，岩石高刚度的成分的初始体积含量和刚度的数值直接决定了弹性模量的大小。

1.3　压密段本构模型的简化

公式(17)用以计算岩石压密段应变时，具有较大的困难，不同岩石的低刚度成分和高刚度成分具有难以确定的难题。因此需要对公式(17)进行适当的简化，简化重点是负指数项，岩石中低刚度成分较为复杂，完全表达出来则会导致公式(17)适用性大大降低。对于脆性岩石，如果其低刚度成分中某一种成分占据主导地位，如热损伤造成的裂纹，则公式(17)可简化为一种低刚度材料。简化后的公式如下：

(21)

式中：Kc为岩石微裂纹法向刚度；n为微裂纹在与最大主应力方向垂直方向上的投影体积比，对于各向同性的岩石材料，n为微裂纹体积比的一半。

为进一步研究裂纹刚度对脆性岩石应力-应变曲线的影响，取弹性模量E=20 Pa，投影体积比n=1%，通过式(21)绘制不同裂纹法向刚度下的应力-应变曲线(如图8所示)。

图8　裂纹法向刚度对应力-应变曲线的影响Fig.8　The influence of normal stiffness of crack on stress-strain curve

通过图8可以发现，随着裂纹法向刚度的增加，应力-应变曲线的压密段逐渐缩小，当裂纹法向刚度达到200 MPa时，已难以辨别微裂纹压密段，目前部分研究人员采用微裂纹闭合应力表示微裂纹完全压缩，但是微裂纹闭合应力往往难以确定，尤其是当微裂纹法向刚度较大时常存在较大误差[1,15-16]。从图8中还可以发现，当微裂纹刚度较小时，弹性模量呈现出随着法向刚度增加减小的趋势，当微裂纹的法向刚度进一步增加时，弹性模量呈现出增加的趋势，事实上岩石真正的高刚度成分的弹性模量是基本保持不变的，微裂纹法向刚度继续增加时，其在一定程度上可以归为高刚度成分，此时脆性岩石将不存在微裂纹压密段。

对于一般的脆性岩石，微裂纹法向刚度通常情况下为相对小值，随着微裂纹刚度增加，弹性模量会出现一定程度的减小。

式(21)可同时表达压密段和线弹性段，因此用式(21)去表达岩石应力-应变曲线，能够更为准确体现出脆性岩石高刚度材料特性，微裂纹性质也能得到较好的表达，更为重要是避免了微裂纹闭合应力的确定难题。对于具体的岩石试验，式(21)表达的是轴向变形，对于环向应变而言，因不存在微裂纹的压密过程，其应力-应变特征可以通过线性去表达。用ε1c和ε2c分别表达微裂纹压密段及线弹性段的轴向应变和环向应变，则其应力-应变可用式(22)表示：

(22)

式中：υ为岩石高刚度材料泊松比。

2　花岗岩热-力耦合室内试验研究

2.1　试样制备及试验过程

为了验证式(22)所示本构模型的适用性，设计花岗岩热-力耦合室内试验。试验选用φ50 mm×100 mm的圆柱形标准试件。试验处理温度为室温，200，400，600，800，900，1 000 ℃ 7种，围压为0，5，10，15 MPa。试验过程如下：首先将花岗岩试样放入到马弗炉中开始热处理，设定升温速率为5℃/min，当温度升高到设定温度时，设定热处理温度恒温4 h，保证花岗岩试样充分受热；停止加热后令试样自然冷却至室温，然后将其封装至三轴压缩试验仪上，给试样施加围压至所需值；然后给试样施加轴压，当试样发生明显应力降时，停止加载，记录并整理相关试验数据。

2.2　试验结果

通过热-力耦合试验，获得不同围压、不同热处理温度下的应力-应变曲线如图9所示。由图可知，随着温度的升高，岩石试样微裂纹压密段逐渐显现，围压会在一定程度上抑制微裂纹压密段的发展；花岗岩在低温低围压条件下，表现出明显的脆性破坏特征，其破坏规律呈现出明显的高温劣化特性，破坏和峰值强度迅速下降。随着热处理温度和围压的升高，花岗岩试样逐渐向延塑性转化。

图9　不同围压条件下的应力-应变曲线Fig.9　Stress versus strain curves under different confining pressures

3　脆性岩石微裂纹压密段本构模型验证

3.1　模拟参数确定

模型参数取值如表1所示，其中E取试样三轴压缩应力-应变曲线弹性阶段的弹性模量，其余参数裂纹体积比n、裂纹刚度Kc通过最佳曲线拟合得到。由表中结果可知，在室温条件下，弹性模量E在37.74～51.61 GPa，裂纹体积比n在0.71～1.50%，裂纹刚度Kc在19.56～38.87 GPa之间，说明在高温预处理之前，花岗岩试样以高刚度成分为主；随着热处理温度升高，花岗岩试样的裂纹增多、刚度劣化，开始向低刚度成分转化，在1 000 ℃时，弹性模量E减小到3.82～21.69 GPa，裂纹体积比n增加到6.45～7.78%，裂纹刚度Kc减小到4.10～28.79 GPa之间。

表1　热-力耦合作用下本构模型验证参数

续表1

3.2　模型验证结果

将第2节试验中压密段和线弹性段的轴压和围压数据代入确定参数后的模型，计算得到相应的应变数据模拟值。图10给出了热处理温度从室温至1 000 ℃、围压为0，5，10，15 MPa下花岗岩三轴压缩试验与数值结果对比。由图可见，本文提出的脆性岩石微裂纹压密段本构模型能够非常好的描述不同热处理温度、不同围压下花岗岩的压缩压密段、线弹性段试验结果，反映了初始阶段花岗岩岩石试样中低刚度成分的裂隙被压缩这一显著变形特征，同时对于脆性岩石的高温劣化性和围压对压密段的抑制作用也有很好的体现。相较于Hui-Hai Liu等[5]、袁小清等[6]、彭俊等[12]提出的脆性岩石本构模型，本文模型中各参数均是在经典的裂纹模型基础上直接推导而来，物理意义明确，对岩石试样三轴压缩应力-应变过程压密段的模拟更为准确，且只有弹性模量E、裂纹体积比n和裂纹刚度Kc3个参数，容易从试验曲线获取。本文提出的脆性岩石微裂纹压密段本构模型，对深入掌握脆性岩石的强度规律具有一定的指导意义。

图10　热-力耦合试验模拟成果验证Fig.10　Validation of thermal mechanical coupling test simulation results

4　结论

1)引入自然应变的概念，建立物理意义明确的脆性岩石微裂纹压密段本构模型，并将其简化后与线弹性阶段的本构模型进行统一。模型共有3个参数，均可通过三轴压缩试验获取。

2)将脆性岩石微裂纹压密段本构模型模拟结果与花岗岩热-力耦合试验结果对比发现，模型能够描述不同热处理温度、不同围压下花岗岩的压密段、线弹性段试验结果，验证了模型的合理性。

3)花岗岩热-力耦合试验表明脆性岩石具有明显的高温劣化性，围压对压密段具有抑制作用，本文所提出的模型能体现这些特性。

4)今后可在脆性岩石微裂纹压密段本构模型的基础上，进一步探讨微裂纹稳定扩展阶段和不稳定扩展阶段，建立应力-应变全过程本构模型。

[1]Eberhart E, Stead D, Stimpson B, et al. Identifying crack initiation and propagation thresholds in brittle rock[J]. Canadian Geotechnical Journal, 1998, 35(2):222-233.

[2]Brace W F, Bombolakis E G. A note on brittle crack growth in compression [J]. Journal of Geophysical Research, 1963,68(12): 3709-3713.

[3]Lu Y L, Elsworth D, Wang L G. Microcrack-based coupled damage and flow modeling of fracturing evolution in permeable brittle rocks[J]. Computers and Geotechnics, 2013, 49:226-244.

[4]周小平, 张永兴. 卸载岩体本构理论及其应用[M]. 北京：科学出版社, 2007:28-89.

[5]Liu H.H, Doughty C, Bodvarsson G S. An active fracture model for unsaturated flow and transport in fractured rocks[J]. Water Resour. Res, 1998, 34(10), 2633-2646.

[6]袁小清，刘红岩，刘京平. 基于宏细观损伤耦合的非贯通裂隙岩体本构模型[J]. 岩土力学, 2015, 36(10): 2804-2814.

YUAN Xiaoqing, LIU Hongyan, LIU Jingping. Constitutive model of rock mass with non-persistent joints based on coupling macroscopic and mesoscopic damages[J]. Rock and Soil Mechanics, 2015, 36(10):2804-2814.

[7]朱其志, 刘海旭, 王伟, 等．北山花岗岩细观损伤力学本构模型研究[J]．岩石力学与工程学报, 2015, 34(3): 433-439.

ZHU Qizhi, LIU Haixu, WANG Wei. A Micromechanical Constitutive Damage Model for Beishan Granite[J]. Chinese Journal of Rock Mechanics and Engineering, 2015, 34(3): 433-439.

[8]刘红岩, 王新生, 张力民, 等. 非贯通节理岩体单轴压缩动态损伤本构模型[J]. 岩土工程学报, 2016, 38(3): 426-436.

LIU Hongyan, WANG Xinsheng, ZHANG Limin, et al. A dynamic damage constitutive model for rock mass with non-persistent joints under uniaxial compression[J]. Chinese Journal of Geotechnical Engineering, 2016, 38(3): 426-436.

[9]朱维申, 赵成龙, 周浩, 等. 当前岩石力学研究中若干关键问题的思考与认识[J]. 岩石力学与工程学报, 2015, 34(4): 649-658.

ZHU Weishen, ZHAO Chenglong, ZHOU Hao, et al. Discussion on several key issues in current rock mechanics[J]. Chinese Journal of Rock Mechanics and Engineering, 2015, 34(4): 649-658.

[10]陈益峰, 胡冉, 周创兵, 等. 热-水-力耦合作用下结晶岩渗透特性演化模型[J]. 岩石力学与工程学报, 2013, 32(11): 2185-2195.

CHEN Yifeng, HU Ran, ZHOU Chuangbing, et al. A permeability evolution model for crystalline rocks subjected to coupled thermo-hydro-mechanical loading[J]. Chinese Journal of Rock Mechanics and Engineering, 2013, 32(11): 2185-2195.

[11]荣冠, 黄凯, 周创兵, 等. 岩石节理法向加载非线性变形本构模型研究[J]. 中国科学: 技术科学, 2012, 42(4): 402-414.

RONG Guan, HUANG Kai, ZHOU Chuangbing, et al. A new constitutive law for the nonlinear normal deformation of rock joints under normal load[J]. Science China Press, 2012, 42(4): 402-414.

[12]彭俊, 荣冠, 周创兵, 等. 岩石裂纹闭合效应及其定量模型研究[J]. 岩土力学, 2016, 37(1): 126-132.

PENG Jun, RONG Guan, ZHOU Chuangbing, et al. A study of crack closure effect of rocks and its quantitative model[J]. Rock and Soil Mechanics, 2016, 37(1): 126-132.

[13]Liu H H, Rutqvist J, Berryman J G. On the relationship between stress and elastic strain for porous and fractured rock[J]. International Journal of Rock Mechanics and Mining Sciences, 2009, 46(2): 289-296.

[14]Freed A D. Natural strain[J]. Journal of Engineering Materials and Technology, 1995, 117(4): 379-385.

[15]Martin C D, Chandler N A. The progressive fracture of Lac Dubonnet granite[J]. International Journal of Rock Mechanics and Mining Sciences and Geomechanics Abstracts, 1994, 31(6), 643-659.

[16]彭俊, 蔡明, 荣冠, 等. 裂纹闭合应力及其岩石微裂纹损伤评价[J]. 岩石力学与工程学报, 2015, 34(6): 1091-1100.

PENG Jun, CAI Ming, RONG Guan, et al. Stresses for crack closure and its application to assessing stress-induced microcrack damage[J]. Chinese Journal of Rock Mechanics and Engineering, 2015, 34(6): 1091-1100.