豫南燕山期花岗岩蠕变特性及非线性蠕变损伤模型

2022-07-10唐志强吉锋许汉华冯文凯何萧

科学技术与工程 2022年16期

唐志强，吉锋*，许汉华，冯文凯，何萧

(1.成都理工大学地质灾害防治与地质环境保护国家重点实验室，成都 610059； 2.中国有色金属工业昆明勘察设计研究院有限公司，昆明 650051； 3.云南省岩土工程与地质灾害重点实验室，昆明 650051)

秦岭-大别造山带作为世界上规模最大的超高压变质带，其后造山阶段的壳幔作用以及造山带的折返、隆升和剥蚀是个复杂的动力学过程；具体的演化主要分为俯冲、同碰撞和后碰撞三个阶段，后碰撞阶段主要出现走滑变形和伸展断裂活动，出现大规模的花岗岩浆活动[1]。信阳鸡公山一带花岗岩属于燕山期花岗岩体的一部分，主要出露在桐商断裂和龟梅断裂以南，北大别构造带上，以高钾钙碱性I型花岗岩为主。伴随着豫南地区众多抽水蓄能电站的建设施工，该地区深埋地下洞室的围岩长期稳定性问题亟待解决。许多研究表明，深埋地下洞室的围岩变形失稳破坏往往与时间密切相关[2-3]。因此，研究该地区特有高钾钙碱性 I型花岗岩的蠕变力学性质，对于豫南燕山期花岗岩地区深埋地下洞室围岩长期稳定性评价及支护设计具有极其重要的意义。

目前针对蠕变模型的研究方法主要是通过对试样的蠕变试验结果进行分析，根据试验蠕变变形表现出的特征，对模型中可能存在的基本元件进行辨识，并将各基本元件进行串、并联等组合，预测在不同应力作用史及变化的物理环境下岩体的变形及破坏与时间的关系[4]。元件模型通常是线性模型，而许多工程实践和试验研究表明，岩石的蠕变往往是一个非线性的过程。鉴于元件模型存在缺陷，国内外学者针对如何更吻合的描述岩石蠕变全过程展开了大量研究。Zhao等[5]通过连接胡克体、胡克体与塑性滑块并联、Kelvin体和广义Bingham体，建立了包含多种变形分量的非线性弹塑性流变模型；Liu等[6]采用分数阶黏壶替代西原正夫模型中的Newton黏壶，从而得到岩石的非线性蠕变本构模型；韩阳等[7]提出了一种非线性黏壶元件替换Burgers模型中的两个线性黏壶元件，建立了一种非定常参数Burgers模型；李任杰等[8]用分数阶微积分软体元件替换Kelvin模型中整数阶黏壶元件,同时引入非定常化非线性黏性元件,采用元件组合方式,构建了岩石硬性结构面蠕变模型。这些研究主要采用非线性元件与传统元件串并联或替代传统元件来建立非线性模型。但非线性模型组合似乎只是为了是理论曲线与试验曲线相吻合，非线性元件的参数却缺乏具体的物理意义[9-10]。

岩石的蠕变破坏本质是岩石在长期荷载作用下，内部裂纹不断产生与发展，岩石材料力学性质逐渐劣化损伤直至破坏的过程[11]。在蠕变模型中考虑材料力学参数随应力状态以及时间的损伤变化，建立岩石的蠕变损伤模型更能反映蠕变的真实状态。李杰林等[12]通过定义基于核磁共振孔隙度的岩石冻融损伤因子与基于体应变的蠕变损伤因子，并将损伤因子引入基本流变元件中，建立了冻融岩石的非线性蠕变损伤本构模型及方程；Lu等[13]通过对缺陷砂岩进行冻融循环试验和三轴压缩试验，建立了冻融砂岩在长期荷载作用下的损伤演化方程；钟祖良等[14]基于试验数据将传统Burgers模型改进为全损伤Burgers蠕变模型，该模型能较好地描述花岗质片麻岩的蠕变特性。综上，中外学者针对特定地质环境下的不同岩石构建蠕变损伤本构模型，并应用于计算模拟，显著提升了岩石变形预测精度。然而，目前针对豫南地区燕山期高钾钙碱性 I型花岗岩的蠕变力学性质的研究仍较薄弱。随着豫南某水电站大型地下厂房洞室开挖，对该地区燕山期花岗岩的蠕变力学特性进行研究显得极为重要。

为此，以该电站地下厂房内的燕山晚期中粒二长花岗岩为研究对象，采用YSJ-01-00型岩石三轴蠕变试验机开展了三轴压缩蠕变实验，分析了花岗岩在不同围压下长期强度的变化规律。构建了能够描述加速蠕变阶段的蠕变损伤模型，并对改进后的西原模型进行了参数辨识及拟合。研究成果可为豫南燕山期花岗岩地区深埋地下洞室围岩长期稳定性评价提供试验依据和理论指导。

1 花岗岩三轴蠕变试验

1.1 研究区花岗岩矿物成分及分布规律

豫南某水电站工程区位于秦岭～大别造山带上，所在区域划分为3个一级大地构造单元，包括中朝准地台、秦岭地槽褶皱系和扬子准地台。进一步可划分为9个二级大地构造单元，即中朝准地台的华熊台缘坳陷、华北坳陷，秦岭褶皱系的北秦岭褶皱带、南秦岭褶皱带、桐柏—大别褶皱带、南阳—襄樊坳陷、潢川坳陷和扬子准地台的上扬子台隆、两湖断坳，如图1所示。

图1 区域大地构造单元划分

该电站的上水库和地下厂房区均位于出露在龟山-梅山断裂(产状为N85°W/SW60°～80°)和桐柏-商城断裂带(产状为N60°W/SW50°～80°)之间，北大别构造带上的燕山期花岗岩株内，为高钾钙碱性I型花岗岩。为揭示此种花岗岩的矿物成分及分布规律，沿燕山期花岗岩株横剖面方向取6个样品(图2)进行磨片测试。

图2 研究区岩性分布及取样位置

根据岩石薄片鉴定结果，研究区花岗岩的矿物成分主要为碱性长石、斜长石、石英，次要矿物见黑云母、白云母，副矿物为锆石和磁铁矿。其中碱性长石含量为33%～39%，呈半自形-他形短柱状；斜长石含量为32%～40%，多为自形-半自形的厚板状；石英含量为21%～25%，无色，他形粒状；黑云母含量为2%～3%，呈褐黄色-褐色多色性，鳞片状，具一组极完全解理；白云母含量为1%～3%，无色，具闪突起，干涉色鲜艳。各类岩石矿物颗粒粒径沿花岗岩株横剖面方向呈先减小后增大的分布规律(图3)。

图3 颗粒粒径分布规律

1.2 蠕变试验设备

本次试验采用的试验设备是由成都理工大学地质灾害防治与地质环境保护国家重点试验室自主研发的YSJ-01-00型岩石三轴蠕变试验机(图4)，该试验机的围压范围为0～30 MPa，轴向荷载范围为0～600 kN，轴向荷载测试精度和围压测试精度均为0.5%F.S，轴向荷载和围压恒定时间可以保持6个月以上，满足本次力学试验的要求。

图4 试验加载系统

1.3 试样制备

本次试验的岩石试样取自河南某水电站大型地下厂房内，岩性为花岗岩；此种花岗岩呈灰白色，中～细粒结构，块状构造，矿物粒径在1～4 mm。参照《水利水电工程岩石试验规程》的试件加工精度要求，制作成Φ50 mm×100 mm的标准圆柱样。试样的高度、直径允许误差为±0.3 mm；试件两端面的不平行度允许偏差为±0.05 mm；端面应垂直于试件轴线，允许偏差为±0.25°。

该水电站地下厂房区岩体应力主要为自重应力场，最大主应力σ1为8.59 MPa，岩石单轴抗压强度平均值为101.17 MPa。为模拟研究区地应力条件下花岗岩的力学特征及蠕变特性，分别在围压σ3=5 MPa和σ3=7.5 MPa条件下进行及常规三轴压缩试验，物理力学特性指标如表1所示。

表1 花岗岩物理力学参数

1.4 蠕变试验方案

基于单轴抗压试验和常规三轴压缩试验结果，分别设置围压σ3=5 MPa和σ3=7.5 MPa两个试验组，加载方式采用单试件逐级加载试验方式“陈氏加载法”[15-16]，根据单轴抗压强度以及不同围压下三轴最大偏应力进行应力水平分级，具体加载应力路径如图5所示。采用荷载控制加载速率方式，加载速率为2.5 MPa/min，考虑到蠕变试验周期较长，每一级荷载的持续时间为24 h，然后重复操作直至试样破坏，停止加载。整个实验过程中需保持围压恒定不变。

图5 不同围压下加载应力路径

2 试验结果分析

2.1 蠕变特征分析

根据花岗岩三轴蠕变力学试验曲线，可分析不同应力状态下花岗岩的蠕变变形随时间变化的规律。图6为花岗岩试样在围压σ3=5、7.5 MPa时的三轴压缩蠕变轴向应变-时间曲线。表2为花岗岩试样在不同围压下的蠕变试验数据。

图6 试样轴向应变-时间

分析表2中数据可知：同一围压条件下，每一级应力增量相等时，随着轴向偏应力的增加，瞬时轴向应变增量逐渐减小，表明试样的弹性模量逐渐增大，试样呈现应变硬化特性；两种围压条件下，轴向蠕变总量与瞬时弹性应变的比值εc/ε0均不超过10%(其中，ε0为瞬时弹性应变，εc为轴向蠕变总量)，表明花岗岩试样在蠕变全过程中轴向总应变以瞬时弹性应变为主；围压7.5 MPa，偏应力42.5 MPa时，轴向蠕变总量εc为0%，表明此时偏应力水平并未达到花岗岩试样的蠕变起始应力，验证了坚硬岩也存在蠕变变形阈值现象[17]；当围压升高时，在相当的偏应力水平作用下，轴向蠕变总量会减小，如围压5 MPa，偏应力150 MPa时，轴向蠕变总量为0.084%，围压7.5 MPa，偏应力170 MPa时，轴向蠕变总量为0.04%；后者的偏应力水平比前者高20 MPa，但轴向蠕变总量却小于前者，表明围压提高，可以制约花岗岩的蠕变变形。

表2 不同围压作用下花岗岩试样三轴压缩蠕变实验数据

2.2 长期强度分析

岩石蠕变长期强度对评价工程长期稳定性而言是极其重要的。通过三轴压缩蠕变试验可知，附加应力水平较小时，岩石的蠕变变形速率会逐渐减小，最终趋于稳定，蠕变变形也会趋于收敛，岩石不会发生破坏。当附加应力水平大于或等于某一定值时，试样的蠕变变形会随时间逐渐增大，且岩石的蠕变速率会在某一阶段迅速增加，试样最终发生破坏。这个应力水平的临界值就认为是岩石的长期强度。

常用的岩石蠕变长期强度的确定方法有3种：一是等时应力-应变曲线法；二是过渡蠕变法(主要根据广义Kelvin模型和Burgers模型对数据的辨识度进行判断)；三是蠕变曲线第一拐点法。

采用等时应力-应变曲线法来确定长期强度。试验时每一级荷载作用持续时间为24 h，分别采取0、2、6、12、18、24 h对应的应变量，绘制六簇等时应力-应变关系曲线，如图7所示。

图7 不同围压下等时应力-应变图

观察六簇等时应力-应变曲线(图7)可知，曲线主要由两部分组成，一部分由左侧斜率较大线段组成，另一部分由斜率较小线段组成，中间存在一个转折点，其对应的轴向应变突然增大，故该转折点所对应的轴向应力即为所求的长期强度值。

当围压为5 MPa时，岩石长期强度为180.23 MPa；当围压为7.5 MPa时，长期强度为168.9 MPa；与常规三轴试验相比，在5 MPa和7.5 MPa条件下，长期强度分别下降了6.37%和21.0%。

2.3 破坏特征分析

花岗岩试件在围压5 MPa、7.5 MPa时，三轴压缩蠕变试验的破坏形态如图8所示。从破坏时间来看，由图6(b)的最后一级蠕变破坏曲线可知，围压5 MPa下试样蠕变总时长约为1.69 h，发生加速蠕变至试样破坏的总历时为0.11 h，其加速蠕变持续时间较为短暂，仅占总时长的6.5%。花岗岩试样的蠕变破坏呈现出脆性破坏特征。

图8 不同围压下试样蠕变破坏形态及素描图

从破坏形式来看，两种不同围压条件下，试样均表现为双剪切面“Y”形剪切破坏，且都发育有一条主贯通裂纹，主裂纹附近存在破坏较彻底的情况，产生大量粉末状固体颗粒，分析其原因是花岗岩在破坏过程中，由于附加应力克服了主裂纹面上的法向应力所产生的摩擦力以及岩石本身的内聚力，导致试样沿主裂纹处发生了强烈的剪切滑移，破坏了剪切面上岩石的晶体结构，所以产生了大量粉末状的固体颗粒。围压7.5 MPa的试样表现出整体湿润，这是由于试样受力破坏时将包裹试样的热缩管炸裂，围压缸内的液压油浸入所致。

由图8可知，围压5 MPa下，试样破坏时2条主裂纹的宏观破坏面夹角约36°，主裂纹交汇处位于试样底部；围压7.5 MPa下，试样破坏时2条主裂纹的宏观破坏面夹角约26°，主裂纹交汇处位于试样中部。对比可知，围压7.5 MPa时，主裂纹沿轴线方向的扩展长度及沿径向的扩展角度均小于围压为5 MPa时，但破坏时的轴向偏应力却比前者大22.5 MPa。因此，相比于轴向偏应力，围压对剪切裂纹的扩展具有更为明显的控制作用。

综上，花岗岩试样在三轴压缩蠕变试验下的蠕变破坏形式为脆性剪切破坏。

3 非线性损伤蠕变模型研究

由图6分析可知，当偏应力施加过程中，岩石试样会产生瞬时弹性变形，则蠕变模型中应存在弹性元件；当附加应力恒定后，岩石试样蠕变变形速率会逐渐减小趋于稳定，蠕变变形量会随时间逐渐减小直至收敛，则蠕变模型中应存在与时间相关的黏性元件；当附加偏应力高于长期强度值且蠕变应变达到某一值时，岩石试样会进入加速蠕变阶段，产生不可恢复的塑性变形并导致最终断裂破坏，则蠕变模型中应存在变形不可恢复的塑性元件。传统西原模型由胡克体元件、黏弹性体及黏塑性体(Newton体)通过串联组合而成，能够较为全面的反应岩石的黏弹塑性变形性质，但对于加速蠕变阶段的描述还不够准确。分析其原因在于传统西原模型中的Newton的黏性系数为定值，然而在实际蠕变试验过程中，岩石会发生损伤，岩石内部本身存在的一些随机分布的微裂纹、孔洞等缺陷在外界影响因素作用下会逐步扩大导致岩石性能劣化，岩石的黏性系数就会发生改变，尤其是加速蠕变阶段的黏性系数与蠕变损伤密切相关[18-19]。

因此，考虑蠕变过程中岩石发生的损伤累积，来描述黏性系数随时间的变化规律，得到一个考虑损伤的变黏性系数黏壶元件，如图9所示。

图9 考虑损伤的变黏性系数黏壶元件

图9中的变黏性系数黏壶元件实质上就是基于损伤累积的Newton体，根据Newton体的本构关系可得出该变黏性系数黏壶元件的本构关系为

(1)

考虑荷载长期作用对黏性系数的影响，可得

η1(D)=η1(1-DC)

(2)

式(2)中：η1为变黏性系数黏壶元件的初始黏性系数；DC为损伤变量，0≤DC<1。

根据文献[20-22]可知，岩石蠕变过程中的损伤变量与时间呈负指数函数关系，即

DC=1-e-αt

(3)

式(3)中：α为与花岗岩材料性质相关的系数；t为时间。

将式(3)代入式(2)中，可得到变黏性系数黏壶元件的黏性系数与时间关系为

η1(D)=η1e-αt

(4)

再将式(4)代入式(1)中，并对其积分，可得变黏性系数黏壶元件的一维本构方程为

(5)

在三维应力作用下，变黏性系数黏壶元件的三维蠕变本构方程为

(6)

式(6)中：ε为变黏性系数黏壶元件的应变值；σS为长期强度值。

通过将上述得到的变黏性系数黏壶元件替换传统西原模型中的Newton体，对于三轴压缩蠕变试验过程中发现的蠕变变形阈值现象，考虑引入一个开关元件来表征，该开关元件以蠕变起始应力σcr为界限；当σ<σcr时，模型表现出弹性元件性质，只能发生弹性变形，无蠕变变形；当σ≥σcr时，该开关元件失去作用，模型发生弹性变形后，在长期荷载的作用下可发生蠕变变形。得到了基于损伤累计的改进西原模型，用于描述岩石蠕变的全过程尤其是加速蠕变阶段的特性。改进的西原模型如图10所示。

E0、E1及η0分别为瞬时弹性模量、黏弹性模量和黏滞系数

基于传统西原模型的本构方程，并结合式(6)，可以推导得出改进西原模型的三维蠕变本构方程为

(7)

式(7)中：K0、G0及G1分别为体积模量、剪切模量和黏弹性剪切模量。

4 蠕变损伤模型验证及参数辨识

为了验证基于岩石损伤累积得到的改进西原模型的适用性和合理性。根据花岗岩三轴压缩蠕变试验结果，对模型进行参数辨识，并得到相应的流变力学参数。借助于七维高科有限公司独立开发的软件1-stOpt，采用Levenberg-Marquardt优化算法(即最小二乘算法)+通用全局优化算法对模型参数进行逐次线性迭代求解，并根据相关性系数反映模型参数与蠕变试验数据的拟合精度。

在常规三轴蠕变试验中，满足σ2=σ3，根据胡克定律有

(8)

限于篇幅，仅对围压5 MPa的三轴压缩蠕变试验数据进行参数辨识，得到各应力水平下的模型参数如表3所示。

表3 围压5 MPa时，改进西原模型参数辨识结果

将计算出的模型参数带入改进的西原模型蠕变方程，得到改进西原模型理论计算曲线和蠕变试验曲线的对比，如图11所示。

图11 围压5 MPa时蠕变试验曲线与改进西原模型理论曲线对比

根据图11可知，改进的西原模型能够对花岗岩的蠕变全过程特征进行较好的描述，准确反映了最后一级加载过程中出现的蠕变三阶段特性，尤其是对加速蠕变阶段的拟合效果也较好，拟合相关系数为0.973 2。因此，改进的西原模型能够较好地描述河南某水电站大型地下厂房洞室花岗岩的蠕变力学特征，可为豫南燕山期花岗岩地区地下洞室围岩长期稳定性评价及支护设计提供试验依据和理论指导。