温度作用下岩石热弹塑性模型及其数值模拟

贾善坡， 林建品， 刘团辉， 崔丽华， 汤小琪， 蒋素萍

(1.长江大学 城市建设学院，湖北 荆州 434023； 2.华北油田 勘探开发研究院，河北 任丘 062552)

温度作用下岩石热弹塑性模型及其数值模拟

贾善坡1,2， 林建品2， 刘团辉2， 崔丽华2， 汤小琪2， 蒋素萍2

(1.长江大学 城市建设学院，湖北 荆州 434023； 2.华北油田 勘探开发研究院，河北 任丘 062552)

摘要：深入研究温度影响下的岩石受力、变形破坏机制及力学模型，对于保证岩石地下工程的安全和稳定具有十分重要的意义.从温度对岩石材料力学性质的影响出发，以修正Mohr-Coulomb准则为基础，建立了温度作用下的岩石热-弹塑性耦合本构模型，探讨了模拟岩石应变软化行为的方法；论述了隐式本构积分算法的主要内容，以ABAQUS软件为平台，编制了材料本构计算子程序.对三轴压缩试验进行了一系列数值模拟，其应力-应变曲线与已有实验规律表现出较好的一致性，所建立的模型能够较好地反映出岩石由于温度影响所表现的脆-塑性过渡，证明了模型的有效性.

关键词：岩石力学；热弹塑性；应变软化；修正Mohr-Coulomb准则

0引言

岩石材料除了遭受外荷载的作用，在很大程度上受环境变化的影响，尤其因温度变化而产生的热应力对其力学性能的影响更为突出.岩石热-力耦合理论在油气勘探开发、地热开发、深部资源开采、核废料地下储存等领域中有着重要的应用，因此多年来它一直是岩石力学的热点和前沿问题之一[1-2].

国内外学者围绕温度对岩石力学性质的影响开展了广泛研究，取得了丰硕成果.文献[3]研究了岩石的杨氏模量随温度升高而变化的情况，安山岩、花岗岩等的杨氏模量在300 ℃以下随温度升高急剧减小，超过300 ℃后杨氏模量变化不大.文献[4]发现花岗岩的弹性模型随温度升高逐渐降低.许锡昌等[5]研究了花岗岩在单轴压缩(20～600 ℃)状态下主要力学参数随温度的变化规律，指出75 ℃和200 ℃分别是花岗岩弹性模量和单轴抗压强度的门槛温度.WONG[6]总结了岩石强度与温度的关系，认为多数岩石随温度升高强度有所下降，下降趋势与岩石的种类相关.在实验室内常温和低围压或单轴压缩情况下，岩石的破坏主要表现为劈裂或剪切破坏形式，这类破裂机制主要受到岩石内部微裂纹的控制.文献[7]研究了温度对材料韧脆转变的影响，发现韧脆转变温度随着加载率升高而升高.桑祖南等[8]在高温高压条件下开展了辉长岩脆-塑性转化及其影响因素的实验研究，辉长岩在600 ℃时以脆性破裂为主，700～850 ℃时为半脆性变形(含微破裂)，900 ℃以上为塑性变形阶段.王颖轶等[9]研究了大理岩在常温至800 ℃ 的应力-应变全过程特性，试验结果表明：随着温度升高，岩石刚度、强度降低，表现出明显的软化特性，峰后特性及残余强度表现出由脆性向塑性的演化.

相关研究表明，在一定的温度和压力作用下，岩石的破坏形式多表现为由脆性破裂转变为塑性流动的现象.然而，现有的岩石热-弹塑性模型尚不能全面地反映岩石的刚度及强度弱化、应变软化以及脆-塑性转换的力学特征.笔者是在过去研究的基础上，建立基于修正Mohr-Coulomb准则的岩石热-弹塑性本构模型，提出模拟岩石应变软化行为的方法，并以ABAQUS软件为平台，编制计算程序，通过一系列三轴压缩数值试验验证程序的有效性.

1岩石热-弹塑性本构模型

1.1热-弹塑性耦合关系

岩石是由不同矿物组成的非均质体，受热后各种矿物的热膨胀效应各不相同，变形大的受压缩，变形小的受拉伸，由此在岩石中形成了一种由温度变化引起的热应力.如果温度由T0变化为T(T>T0)，则岩石将发生膨胀，设拉伸为正，则热应变可表示为

εT=β(T-T0)I，

(1)

式中：β为膨胀系数；I为单位张量.

温度影响下的岩石应变增量可表示为[2]

dε=dεe+dεp+dεT，

(2)

式中:dε、dεe、dεp、dεT分别为总应变增量、弹性应变增量、塑性应变增量、温度应变增量.

岩石弹性模量和泊松比均为温度的函数，应力与弹性应变满足胡克定律σ=De:εe，增量形式为

dσ=dDeT:εe+DeT:dεe,

(3)

式中:DeT为弹性刚度矩阵，是温度的函数,DeT=2GI+λI,G为剪切模量,λ为拉梅系数.

温度对岩石力学性质的影响包括两个方面：①温度的变化对岩石物理力学参数有影响，温度升高导致岩石劣化，弹性模量、强度等均随着温度的升高而下降；②温度的变化会导致热应力问题.综合利用式(1)～式(3)，可得岩石热-弹塑性增量形式的本构关系为

(4)

1.2修正Mohr-Coulomb模型

常规的Mohr-Coulomb是一个压剪型的屈服准则，不适用拉剪型应力状态，用它预测的锥顶处的拉力值大于实验值(抗拉强度)[10].对常用的Mohr-Coulomb 准则进行修正，采用双曲线方程对拉伸型Mohr-Coulomb屈服准则和剪切型Mohr-Coulomb屈服准则进行拟合(图1)，即

(5)

根据式(5)，构造基于修正Mohr-Coulomb准则的塑性势函数，即

(6)

式中:φ为膨胀角.

1.3应变软化模型

岩石材料在常温和低围压或单轴压缩情况下多表现为峰后区弱化，如果采用强度准则描述这种过程，岩石应力达到峰值强度后，随着变形的进一步增加，岩石强度由峰值(或快或慢)下降至强度的残余值，强度参数发生弱化[11].用ω来代表岩石强度参数(如岩石黏聚力c、内摩擦角φ等)，通过建立应变软化模型来模拟岩石随着温度的升高由脆性破坏到延性破坏的机制转变.图1为强质参数随等效塑性应变的溶化曲线，其溶化方程可定义为

(7)

峰值强度随温度的升高有逐渐降低的趋势，为了简化计算，假设峰值强度ωmT与温度T之间满足线性关系，即

ωmT(T)=ωm0-B·T，

(8)

式中:ωm0、B均为模型参数，B>0.

在一定的温度和压力作用下，岩石的主要破坏形式会由脆性破裂向塑性流动转变.大量试验结果表明，残余强度随温度的升高有逐渐增大的趋势[2]，定义残余强度ωrT为温度的函数(图2)，

(9)

2数值计算积分方案

温度对岩石力学的影响属于非线性热弹塑性问题，由温度引起的热应力、塑性流动以及随温度变化的屈服准则等结合在一起，屈服面随着温度升高而逐渐减小.笔者采用隐式向后欧拉积分算法实现计算过程中的应力更新，如图3所示.其基本原理是：首先是弹性预测阶段，假设增量步的应变增量全部为弹性应变，计算出应力更新后的数值；其次是塑性修正阶段，检查弹性预测阶段的应力更新数值是否满足屈服条件，若不满足屈服条件，则无需进行修正，反之则对应力进行修正，使应力返回到屈服面[12].

在数值求解的过程中，通常以时间为增量步寻求在离散时间点时刻的平衡，在每个时间步内，屈服函数随温度升高而逐渐减小，屈服函数不仅与时间有关，还与等效塑性应变和温度有关.

(1)根据温度场计算结果确定tn+1时刻的温度，根据计算公式确定变温下的材料力学特性(弹性刚度矩阵、黏聚力、内摩擦角等).

(2)在增量步开始时，主程序提供应力张量σn、总应变εn+1、总应变增量Δε、及时间增量Δt，计算热应变增量ΔεT和应变增量Δεep=Δε-ΔεT.

(3)根据弹性关系计算试探应力的增量值Δσ*和试探应力值σ*n+1.

(10)

式中:DeT为弹性矩阵.

3数值算例

数值实验选用标准圆柱形岩样(见图4)，几何模型尺寸为50mm×100mm(半径×高)，侧面施加围压，垂直方向采用位移控制，上面施加非零位移，下面施加位移约束，与岩样实验中采用的应变控制加载过程一致.选取某深埋隧洞围岩大理岩作为试样[13]，单轴压缩试验的弹性模量、泊松比、峰值应力与温度之间的关系式为

(11)

式中:弹性模量的单位为GPa;峰值应力的单位为MPa.

根据Mohr-Coulomb准则，可建立黏聚力与单轴抗压强度的关系式σc=2c·cosφ/(1-sinφ)， 因此，峰值黏聚力cmT与温度的关系式为

(12)

针对如下3种工况进行分析：①单轴压缩，采用理想弹塑性模型；②单轴压缩情况，采用应变软化模型；③温度T=200 ℃，三轴压缩情况，采用应变软化模型.

以ABAQUS软件为平台，对上述模型进行二次开发，对圆柱形试样模型进行有限元计算，图5～图7为岩样在不同工况下的应力-应变曲线.

工况1情况下应力-应变曲线规律与大量试验结果的规律一致，随着温度的升高，屈服强度逐渐减小，弹性模量逐渐减小；温度越高，屈服流动的效应增长，数值计算较易收敛.

对于工况2而言，当岩石进入屈服阶段后，随着塑性应变的增加，黏聚力逐渐减小直至残余强度值，在温度较低时表现为弹脆塑性性质，随着温度的升高，由脆性向韧性过渡，当温度达到300 ℃时，岩石已开始变现为应变软化特性，当达到500 ℃时，岩石的塑性有理想塑性流动的趋势.由此可以看出，笔者提出的热-弹塑性应变软化模型可以较好地模拟由于温度的影响岩石所表现的弹-脆塑性到弹塑性的过渡.

对于工况3而言，随着围压的升高，岩石的屈服强度逐渐升高，在高围压下，岩石有脆性向韧性过渡的趋势，与已有的岩石力学试验规律一致.

4结论

(1)以热力学和弹塑性理论为基础，根据已有实验的规律性认识，对服从Mohr-Coulomb准则的岩石强度参数演化行为进行了研究，建立了温度作用下的岩石热弹塑性本构模型，并以ABAQUS软件为平台，提出了本构数值积分算法.

(2)对三轴压缩试验进行了一系列数值模拟，得到了不同温度和压力条件下的应力-应变曲线，模拟结果与实验的规律性认识一致，热-弹塑性应变软化模型可以较好地反应由于温度和压力的影响岩石所表现的脆-塑性过渡，证实了该模型的有效性.

参考文献：

[1]周宏伟, 谢和平, 左建平. 深部高地应力下岩石力学行为研究进展[J]. 力学进展, 2005, 35(1): 91-99.

[2]康健.随机介质固热耦合数学模型与岩石热破裂数值实验[D]. 阜新：辽宁工程技术大学力学与工程学院，2004.

[3]林睦曾，岩石热物理学及其工程应用[M]. 重庆：重庆大学出版社，1991：130-138.

[4]刘月妙，王驹，谭国焕，等. 热力耦合条件下北山深部花岗岩长期性能研究[C]// 第九届全国岩石力学与工程学术大会论文集. 北京：科学出版社，2006：31-36.

[5]许锡昌, 刘泉声. 高温下花岗岩基本力学性能初步研究[J]. 岩土工程学报, 2000, 22(3): 332-335.

[6]WONGTeng-fong.Effectsoftemperatureandpressureonfailureandpost-failurebehaviorofwesterleygranite[J].MechanicsofMaterials, 1982, (1): 3-17.

[7]BREDEM.Thebrittle-to-ductiletransitioninSilicon[J].ActaMetallMaster, 1993, 41(1): 211.

[8]桑祖南，周永胜，何昌容，等. 辉长岩脆-塑性转化及其影响因素的高温高压实验研究[J]. 地质力学学报，2001，7(2)：130-137.

[9]王颖轶，张宏君，黄醒春，等. 高温作用下大理岩应力-应变全过程的试验研究[J]. 岩石力学与工程学报，2002, 21(2)：2345-2349.

[10]贾善坡,陈卫忠,杨建平,等. 基于修正Mohr-Coulomb准则的弹塑性本构模型及其数值实施[J]. 岩土力学，2010,31(7):2051-2058.

[11]王水林, 王威，吴振君. 岩土材料峰值后区强度参数演化与应力-应变曲线关系研究[J].岩石力学与工程学报，2010，29(8)：1525-1529.

[12]贾善坡.BoomClay泥岩渗流应力损伤耦合流变模型、参数反演与工程应用[D]. 武汉: 中国科学院武汉岩土力学研究所, 2009.

[13]施毅,朱珍德,李志敬.考虑温度效应时深埋洞室围岩变形特性[J].水利水电科技进展,2008,28(3):33-36.

AThermo-elasto-plasticConstitutiveModeloftheRockundertheTemperatureEffectandItsNumericalImplementation

JIAShan-po1，2，LINJian-pin2，LIUTuan-hui2，CUILi-hua2，TANGXiao-qi2，JIANGSu-ping2

(1.SchoolofUrbanConstruction，YangtzeUniversity，Jingzhou434023，China； 2.ExplorationandDevelopmentResearchInstituteofHuabeiOilfield,Renqiu062552,China)

Abstract:Further study of the mechanism and law of rock deformation and breakage under the temperature effect has great significance to the safety and stability of the underground rock engineering. Based on the evolution characteristics of rock mechanics parameters under the temperature effect, a new thermo-elasto-plastic constitutive model with improved Mohr-Coulomb criterion is presented, which can depict the temperature effect and softening behaviours of the rock. The main contents of the implicit constitutive integration algorithm are discussed systematically; and the consistent stiffness matrix is deduced. Based on the implicit Euler stress integration algorithm, the UMAT subroutine of the model is developed in ABAQUS. The numerical simulation under various stress states and temperatures shows that the numerical simulations agree with the experimental rules and the proposed model is able to depict effectively the main features of brittle-plastic transition of rock.

Key words:rock mechanics; thermo-elasto-plastic; strain-softening; modified Mohr-Coulomb criterion

中图分类号：TU457

文献标志码：A

doi:10.3969/j.issn.1671-6833.2015.02.012

文章编号:1671-6833(2015)02-0052-05

作者简介：贾善坡(1980-)，男，长江大学副教授，博士后，主要从事岩石力学与工程方面的教学与研究工作,Email:jiashanporsm@163.com.

基金项目：中国博士后科学基金资助项目(2014M551055)； 岩土力学与工程国家重点实验室开放基金资助项目(Z013007);湖北省建设厅科技计划资助项目(2011-119)

收稿日期:2014-10-25；

修订日期：2014-12-03