高温花岗岩热冲击后力学特性及损伤演化规律研究

郭奇峰 ，钱志海 ，潘继良 ，席 迅 ✉，蔡美峰

1) 北京科技大学土木与资源工程学院，北京 100083 2) 北京科技大学金属矿山高效开采与安全教育部重点实验室，北京 100083

地热资源作为一种绿色清洁能源，近年来得到了世界各国的广泛关注. 传统的水热型地热资源已经得到了大面积的利用，而干热型地热资源的开发尚处于试验探索阶段[1-3]. 干热型地热资源是指赋存于地下3～10 km范围内温度在180 ℃以上的高温岩体，具有低孔隙度或低渗透性的特点，需要建造人工储层等手段进行热能提取[4-5].在传统的增强型地热系统（Enhanced geothermal system，EGS）工程中，通常从注入井将高压冷水介质灌输进深部干热岩储层，通过水力压裂、热刺激和化学刺激等方式形成人工裂隙网络，最终目的在于实现流体介质在注入井与生产井之间的有效流通和热量交换[6-8]. 在此期间，注入冷水介质诱发的热冲击可能会导致干热岩储层产生热破裂，这有助于在注入井附近形成裂缝. 然而，该过程也为高温岩体带来了诸如井壁塌孔和局部微震等一系列岩石损伤问题[9-13]. 因此，开展高温花岗岩的物理力学特性研究具有重要现实意义.

当前，学者们对高温花岗岩物理力学特性的研究主要集中在两个方面[14-16]：一是高温花岗岩在冷却作用后的物理力学特性研究；二是实时高温下花岗岩的物理力学特性研究. 受限于试验条件和设备，大多数研究集中于花岗岩在达到指定温度后冷却至室温下的力学性质. 冷却方式包括自然冷却、遇水冷却和其他条件冷却（如液氮等），其中以强烈温差形成的热冲击对岩体力学性能劣化最为明显[17-20]. 为了合理描述高温岩体在力学性能劣化后抵抗外力变形的能力，可通过建立本构模型来反映热损伤岩体在外力作用下的应力-应变关系. 例如，赵齐等[21]提出了由声发射能量表征的力损伤变量和由弹性模量表征的热损伤变量，建立了基于声发射技术的花岗岩热-力耦合损伤本构模型；朱振南等[22]基于正态分布函数和Mohr-Coulomb准则，通过引入热损伤变量建立了岩石统计热损伤本构模型；闵明[23]和蒋浩鹏等[24]基于Weibull分布函数，分别引入Drucker-Prager准则和Mohr-Coulomb准则来描述加载期间试件微元体的破坏，建立了热力耦合作用下的岩石统计损伤本构模型. 以上关于高温岩体的本构模型均考虑了高温热损伤作用，且理论计算获得的应力-应变曲线在整体上均能够较好地贴合试验曲线，但在初始非线性压密阶段的贴合效果较差，无法合理反映强烈热冲击引起的岩体孔隙结构劣化效应.

本文以EGS工程注入井近场岩体为研究对象，考虑注入冷水介质产生的热冲击效应，通过室内试验对花岗岩试件进行不同温度水平的热处理，然后采用快速高温水冷的方式完成热冲击. 首先，通过对热冲击花岗岩试件开展单轴压缩试验获得应力-应变曲线，为后续本构模型的验证提供数据支撑；然后，基于损伤力学理论，提出一种考虑高温初始损伤与加载期间试件微元破裂损伤相结合的热-力耦合损伤本构模型，并对本构模型的相关参数进行理论求解；之后，考虑热冲击损伤引起的孔隙结构劣化效应，引入压密系数对热冲击花岗岩的本构关系进行修正；最后，将提出的本构模型拟合曲线与试验应力-应变曲线进行对比和验证，讨论不同温度热冲击花岗岩试件在单轴压缩荷载下的损伤演化过程. 研究成果不仅有助于了解热冲击花岗岩宏观力学性能的损伤劣化过程，而且可为构建更加准确的数值计算模型和工程方案论证提供重要的理论指导.

1 试件制备流程

将花岗岩试样研磨成粉末状后，通过X射线衍射分析对样品成分进行测试，得到花岗岩试件的主要成分为石英和长石（钠长石、钙长石、钾长石）矿物. 将花岗岩加工为力学试验所需的标准圆柱体试件，试件的高度和直径分别约为100 mm和50 mm，高径比为2∶1，通过对试样的端面进行抛光处理，确保2个端面的粗糙度小于0.2 mm，不平行度在±0.1%以内，使加工后的试件满足国际岩石力学学会（ISRM）建议标准.

综合考虑干热型地热资源开发的经济性和现有技术条件，本文以600 ℃以内温度区间的干热岩为研究对象，选取的温度水平包括25、150、300、450和600 ℃. 注入高压冷水介质后，流体首先与高温干热岩之间产生快速的热交换，该工程环境在室内试验中即可简化为高温岩体快速遇水冷却形成的热冲击. 因此，在试验设计中，以温度梯度作为试验变量，通过高温热处理和快速遇水冷却制备热冲击花岗岩试件，之后开展相应的岩石力学试验.

为了降低加温过程中由不均匀温度分布引起的冲击热应力对岩石热膨胀行为的影响，采用马弗炉对试件进行缓慢升温处理，升温速率控制在5 ℃·min-1，当达到各个预设温度时，保持恒温状态2 h，其目的在于保证花岗岩试件内部受到均匀的热应力作用，同时可以确保在高温加热期间发生充分的物理化学反应. 之后，将热处理后的高温试件从马弗炉里取出，快速放入20 ℃的冷水中进行热冲击，水冷时间1 h，待试件完全冷却后取出擦干，然后移入105 ℃的恒温箱中进行12 h的烘干处理，直至试件充分干燥，避免不同含水率对试件宏观力学性质的影响.

不同温度热冲击作用后花岗岩试件的表面形貌变化如图1所示. 可以看出，150 ℃热冲击后花岗岩的表面形貌与常温状态下的花岗岩区别不大；300 ℃及以上热冲击试件表面均出现了大量肉眼可见的龟裂状微裂纹；随着温度的升高，由于脱水及部分矿物的氧化分解作用，试件表面的浅灰色逐渐转变为灰白色，在超过300 ℃以后，浅灰色花纹逐渐消失，只保留部分深黑色花纹，说明在高温热冲击作用后，部分矿物成分发生了转变.

图1 不同温度热冲击后花岗岩试件表面形貌变化Fig.1 Surface morphology changes in granite specimens after thermal shock at different temperatures

2 应力-应变曲线分析

使用如图2所示的MTS815电液伺服岩石力学试验系统开展单轴压缩试验，对热冲击花岗岩试件的宏观力学性质开展研究. 首先对试件进行预加载，加载速率为0.05 kN·s-1，预压力为0.5 MPa，保证试验机压头与试件完全接触；之后，以应力控制方式施加轴向压力，加载速率为0.25 MPa·s-1；当轴向应力加载至体积应变接近转折点时，转换为环向位移控制方式继续加载，加载速率为0.015 mm·min-1，直至试件发生破坏.

图2 MTS电液伺服岩石力学试验系统Fig.2 MTS electrohydraulic servo rock mechanics test system

试验获得的热冲击花岗岩试件单轴压缩应力-应变曲线如图3所示. 可以看出，随着温度的升高，峰值应力逐渐降低，峰值应变逐渐增大（150 ℃除外）. 完整的应力-应变曲线可分为初始压密阶段、弹性变形阶段、非稳定破裂发展阶段和峰后破裂阶段，与常温花岗岩试件的区别在于，热冲击试件在初始压密阶段的应力-应变曲线更多的呈下凹型，主要与初始热损伤微裂纹或微孔洞的闭合压密有关. 由于常温花岗岩本身较为致密，在25～150 ℃温度区间内热损伤尚未形成，此时的初始压密阶段几乎可以忽略不计. 当热冲击温度达到300 ℃及以上时，由于热应力的作用、脱水作用和部分矿物成分的氧化分解，花岗岩试件出现较为严重的热损伤现象，且随着热冲击温度的升高损伤持续增加，初始压密阶段愈发明显. 因此，在构建热冲击花岗岩的单轴压缩本构模型时，应当充分考虑加载初期由于孔隙结构劣化引起的非线性压密现象.

图3 热冲击花岗岩试件单轴压缩应力-应变曲线Fig.3 Uniaxial compressive stress–strain curve of granite specimens after thermal shock

3 热-力耦合损伤本构模型

3.1 损伤本构方程

Kachanov[25]将连续度定义为有效承载面积与无损状态下的横截面面积之比，即：

式中：A为 无损状态时的横截面面积，m2；为损伤后的有效承载面积，m2；ψ为连续度，是一个量纲一的标量场变量， ψ=1对应于不存在任何损伤的理想材料状态， ψ=0对应于完全破坏的没有任何承载能力的材料状态.

将外加载荷F与有效承载面积A˜之比定义为有效应力σ˜，即：

式中，σ为Cauchy应力，是外加载荷F与无损面积A的比值.

基于连续度ψ，Rabotnov[26]提出了损伤因子D的概念，将其定义为：

对于完全无损状态，D=0；对于完全丧失承载能力的状态，D=1. 于是可以得到有效应力σ˜与损伤因子D之间的关系为：

根据Lemaitre应变等价性假设，认为损伤单元在名义应力作用下的应变响应与无损单元在有效应力作用下的应变响应相同，在外力作用下受损材料的本构关系可采用无损时的形式，只需要把其中的Cauchy应力替换成有效应力即可，即：

式中，ε为 应变，E和分别为无损材料和受损材料的弹性模量.

因此，可以建立基本损伤本构方程：

张全胜等[27]将岩石的天然损伤定义为基准损伤状态，提出推广后的应变等价原理：材料受到力F的作用，损伤产生扩展，任取其中的2种损伤状态，则材料在第1种损伤状态下的有效应力作用于第2种损伤状态引起的应变等价于材料在第2种损伤状态下的有效应力作用于第1种损伤状态引起的应变，即：

式中： σ(1)和 σ(2)分别为第1种和第2种损伤状态下的有效应力；E(1)和E(2)分别对应为第1种和第2种损伤状态下的弹性模量.

根据推广后的应变等价原理，不妨将岩石的基准损伤状态作为第1种损伤状态，热冲击损伤后的状态作为第2种损伤状态，有：

式中：σ0和σT分别为基准损伤状态和热损伤状态的有效应力；E0和ET分别对应为基准损伤状态和热损伤状态的弹性模量；DT为热损伤变量.

于是得到热冲击损伤状态的基本损伤本构方程为：

弹性模量E0和ET的关系式为：

同理，将热损伤后的状态作为第1种损伤状态，热损伤后加载引起的损伤状态作为第2种损伤状态，再次应用推广后的应变等价原理可得材料内部热损伤与加载受荷损伤共同作用下的本构关系为：

式中，σ*为 热损伤后加载期间的有效应力；DM为 加载引起的损伤变量.

联立式(10)和式(11)，可得用热损伤变量和受荷损伤变量表示的花岗岩受荷应力-应变关系为：

式中，DT-M为热损伤花岗岩受载期间的总损伤变量，DT-M=DT+DM-DTDM.

可以看出，花岗岩在受荷期间的总损伤DT-M是初始热损伤DT和 加载损伤DM相互叠加或耦合的结果，DT+DM为 叠加项，DTDM为耦合项.

根据损伤力学中的宏观唯象学方法，岩石宏观物理性能的响应能够代表材料内部的劣化程度[28].材料的弹性模量在热作用前后更便于分析和测量，可将岩石热损伤变量定义为：

加载损伤主要由加载期间岩石材料微元的失效破坏引起的，将加载损伤变量DM定义为失效的微元数量Nd与 微元总数Nt的比值，可得：

微元的破坏同样满足强度准则，可以表示为：

其中：f(σ*)为有效应力的某种组合函数，可用微元强度分布变量S表示；K0为常数，代表岩石材料的强度. 当微元的应力组合达到其强度时，即发生失效破坏.

岩石本身即是一种非均质性材料，破坏微元的分布具有随机性，因此采用统计损伤力学的方法研究加载期间微元破坏引起的损伤更为合理.在任意载荷水平作用下，累计失效微单元的总数量为加载区间内所有随机分布的破坏微单元的积分，即：

式中：x为分布变量；P(x)为概率分布函数.

于是可以得到加载损伤变量DM的表达式为：

假定花岗岩微元的强度失效概率满足Weibull分布函数：

其中，m和S0是分别反映微元集中程度和强度大小的Weibull参数.

将分布函数代入式(17)中，积分得到加载损伤变量DM的表达式为：

下面将采用Drucker-Prager强度准则确定岩石微元强度分布变量S，并对Weibull分布参数m和S0的值进行求解.

Drucker-Prager强度准则表达式为[29]

式中：σ*1、 σ*2和σ*3分别为不同主应力方向的应力分量；I1为 应力张量第一不变量；J2为应力偏张量第二不变量；α0和k为 材料参数；φ为岩石材料的内摩擦角.

引入应变等价性假设，由有效应力σ*替换线弹性Hooke定律中的表观主应力，在单轴压缩状态下，可以得到如下本构关系：

式中：σ1和 ε1分别为轴向应力和轴向应变.

可以得到：

于是可以求出岩石微元强度分布变量S的表达式为：

通过对式(24)进行变换形式可以得到加载损伤变量DM的表达式

式(28)与式(19)是等价的，即：

将式(27)的微元强度分布变量S代入式(29)中，化简可得：

式(30)即为单轴压缩条件下基于Weibull分布的热冲击花岗岩统计损伤本构方程.

3.2 统计参数的确定

岩石材料的强度参数（内聚力c和内摩擦角φ）可以通过真实的岩石力学试验获得，Weibull分布参数m和S0可以通过线性回归法、灾变理论和峰值点法等进行确定[30]. 对于岩石的压缩应力-应变曲线，其中最重要的特征点为峰值点，可以反映岩石的抗压强度和峰值应变. 通过峰值点法对统计损伤本构模型的参数进行求解，在应力-应变曲线的峰值点处，存在以下关系：

式中，σp和εp分别为峰值轴向应力和峰值轴向应变.

对式(30)的本构模型取相对于所述轴向应变ε1的导数，有：

由式(27)求出在峰值点处的微元强度分布变量Sp，可得：

代入到峰值条件关系式(31)和(32)中，可得：

于是可以得到：

峰值点处的应力-应变关系可以表示为：

将式(37)代入式(35)，求得Weibull参数m的表达式为：

将m代入式(36)即可得到Weibull参数S0为：

从应力-应变曲线可以看出，热损伤试件在压缩变形初期均会出现更为显著的非线性压密阶段，因此，通过引入压密系数K对热损伤试件的本构关系式(30)进行修正，修正后的单轴压缩损伤本构方程如下：

压密系数K可按下式取值[31]：

式中，εcd为屈服应力对应的应变；n为常数，与应力-应变曲线的曲率有关，n>0且n≠1.

3.3 损伤本构模型验证

表1给出了热冲击花岗岩试件单轴压缩统计损伤本构模型对应的参数，图4分别给出了不同温度作用下热冲击花岗岩试件的单轴压缩应力-应变曲线与本构模型理论曲线拟合效果对比. 从拟合曲线的对比结果可以看出，对于常温25 ℃和150 ℃热冲击试件，由于没有形成明显的热损伤，初始非线性压密特征较弱，未修正的统计损伤本构模型与试验应力-应变曲线偏差不大. 随着温度的升高，热冲击效应不断增强，初始非线性压密阶段越来越明显，未修正的本构模型与试验应力-应变曲线的偏离程度逐渐加大，且温度越高偏离幅度越大，此时已无法通过该模型反映真实的岩石本构关系. 对比可以看出，在300 ℃及以上温度热冲击条件下，引入压密系数进行修正的统计损伤本构模型则能够很好地拟合真实应力-应变关系，即使是在没有明显热损伤的条件下，其拟合效果依然优于未修正的损伤本构模型. 由于试验中采用的是环向位移控制加载方式，峰后表现为ClassⅡ曲线的变形特征，因此无法对损伤本构模型的峰后阶段进行验证. 验证结果表明，本文提出的考虑初始非线性压密的统计损伤本构模型可以合理地反映热冲击花岗岩试件的峰前应力-应变关系.

图4 热冲击花岗岩单轴应力-应变曲线与本构模型理论曲线对比. (a) 25 ℃; （b）150 ℃; （c）300 ℃; （d）450 ℃; （e）600 ℃Fig.4 Comparison between the uniaxial stress–strain curve of granite and the theoretical curve of the statistical damage constitutive model: (a) 25 ℃;（b）150 ℃; （c）300 ℃; （d）450 ℃; （e）600 ℃

表1 热冲击花岗岩试件单轴压缩统计损伤本构模型参数Table 1 Parameters of the statistical damage constitutive model of granite specimens after thermal shock under uniaxial compression

4 讨论

单轴加载期间不同温度热冲击花岗岩试件的损伤变量演化特征如图5所示. 可以看出，在经过150 ℃的热冲击后，试件的初始热损伤略有下降，分析认为这是由于高温脱水和热膨胀作用增加了晶体的结合力和矿物颗粒间的摩擦力，导致弹性模量增大. 在150 ℃以上温度区间，随着温度的升高，初始热损伤不断增大. 在加载初期，总损伤变量演化曲线基本保持不变，此时主要以初始微裂纹或微孔隙的压密闭合为主，尚未形成新的损伤；在达到某一损伤阈值后，伴随着新裂纹的萌生和扩展，损伤变量随着轴向荷载的增加而快速上升，但不同温度水平热冲击试件的上升曲线形态存在较大差异. 在热冲击温度水平较低时，损伤变量演化曲线上升形态较为陡峭；随着温度的升高，初始热损伤程度不断加深，外载作用下的损伤阈值逐渐前移，在接近峰值时的损伤变量曲线上升速率也逐渐变缓. 例如在600 ℃时，损伤变量曲线随轴向应变的增加几乎呈线性上升. 曲线的发展形态也从侧面说明了热冲击花岗岩试件在低温条件下的破坏主要以脆性为主，而温度的增加使得其破坏形式逐渐向延性和塑性转变.

图5 单轴加载期间热冲击花岗岩试件损伤变量演化特征Fig.5 Evolution characteristics of the damage variables of granite specimens after thermal shock during uniaxial loading

5 结论

为了研究高温花岗岩在热冲击作用后的力学特性和损伤演化规律，开展了25～600 ℃范围内不同温度热冲击作用下花岗岩试件的单轴压缩试验，提出了一种考虑初始热冲击损伤与加载期间试件微元破裂损伤相结合的统计损伤本构模型，引入压密系数对本构模型进行了修正，并通过试验获得的应力-应变曲线对本构模型的有效性进行了对比和验证，讨论了温度对加载期间热冲击花岗岩试件损伤变量演化特征的影响. 主要结论如下：

（1）随着热冲击温度的升高，花岗岩试件的初始热损伤不断增大，应力-应变曲线表现出明显的非线性压密现象；

（2）引入压密系数修正后的统计损伤本构模型能够更加充分地反映热冲击花岗岩试件在初始加载阶段的非线性压密特征；

（3）在热冲击温度较低时，损伤变量演化曲线上升形态较为陡峭，随着热冲击温度的升高，曲线上升速率逐渐变缓并逐渐由非线性向线性转变.