基于温度损伤的红砂岩本构模型

A Constitutive Model of Red Sandstone Based on Temperature Damage

摘要: 本文针对深部地下岩体工程中岩石处于高温环境，且所处的三向应力状态一般为σ1 > σ2 > σ3。考虑了岩石高温下的损伤和真三轴应力条件下的岩石的应力损伤。基于Lemaitre应变等效假设理论，通过引入能够反映岩石在热力耦合条件下的温度损伤变量和力损伤变量来描述岩石材料的损伤程度，根据热力耦合损伤变量的指数分解形式引入热力耦合损伤变量，基于D-P准则建立岩石高温–荷载耦合损伤的真三轴瞬时本构模型。

Abstract: This study focuses on the high-temperature environment encountered in deep underground rock engineering, where the triaxial stress state is generally σ1 > σ2 > σ3. Considering both the damage to rocks under high temperatures and the stress damage under true triaxial stress conditions, based on Lemaitre’s strain equivalence hypothesis theory, the temperature damage variable and mechanical damage variable, which can reflect the damage degree of rock materials under thermo-mechanical coupling conditions, are introduced to describe the damage extent of rock materials. The thermo-mechanical coupling damage variable is introduced according to the exponential decomposition form of the thermo-mechanical coupling damage variable. True triaxial instantaneous constitutive models for high temperature and load coupling damage of rocks are established based on the D-P criterion.

1. 引言

岩石是由一种或多种矿物组成的具有不同结构构造的非均质天然材料。由于矿物成分、地层条件、胶结材料及地质构造强度的差异性，岩石内部存在大量随机分布的缺陷，使得构成岩石细观的微单元强度有极大差异。岩石内部结构存在大量微裂缝和初始孔隙，处于随机分布且数量众多，这些均质的细观破坏只能用统计方法来研究。曹文贵等[1]假定岩石所受的轴向荷载由损伤和未损伤材料共同承担，利用统计损伤理论，把岩石损伤阈值影响纳入考虑范围，提出一种反映岩石统计损伤全过程的本构模型；孙梦成等[2]假设岩石各项同性，基于连续损伤理论利用最小耗能原理推导了复杂条件下的岩石线性损伤演化方程；石崇等[3]假定岩石微元强度服从Weibull分布，基于H-K准则建立三维应力作用下的岩石损伤模型，结合岩石应变软化特点引入修正参数q；朱振南等[4]引入热损伤变量，岩石基于Normal分布和Mohr-Coulomb准则条件下建立热力作用下岩石损伤本构模型；高玮等[5]通过得到岩石破裂程度变量来建立一种能反映深部工程破裂区破裂岩体的本构模型；张卫强[6]开展了不同高温处理后的岩石孔隙率测试，实验表明岩石的热损伤主要是由于岩石自身矿物结构的改变以及热应力造成的孔洞缝隙的发育而造成的；杨文东[7]开展了实时高温和热处理后的三轴压缩试验，结果表明岩石的弹性模量随温度的增大而减小，但是泊松比呈现不同的变化规律；Chen [8]研究了岩石物理力学性质显著改变时的阈值温度，发现随着加热温度的升高，受热后岩石峰值应力和弹性模量减小，峰值应变增大；Krajcinovic和Silva [9]从岩石材料内部缺陷的随机分布出发，利用岩石的微量元素强度遵循Weibull分布，建立了岩石断裂过程的统计损伤本构方程；Xu [10]基于Weibull分布和Lemaitre应变等效原理，为了考虑热力耦合对岩石的损伤，引入一个非线性的热力耦合损伤参数变量；郤保平等[11]通过对高温花岗岩进行不同冷却，发现岩石的强度特性不仅与温度有关，加热方式、升温速度、降温速度、冷却方式等因素都会使岩石的强度特性发生改变；成泽鹏等[12]通过对热冲击过程中的岩石力学参数、微观结构、流动性和破裂过程数值分析，发现花岗岩的渗透率变大，突变存在阈值，随着温度的升高，花岗岩的渗透率和孔隙率先缓慢增加后急剧增加。本文在前人基础上，基于损伤变量指数分解和D-P准则，考虑热力共同作用建立三轴条件下岩石统计损伤本构模型。带入数据进行验证，以期进一步得到岩石温度和围压共同作用下的损伤破坏规律。

2. 岩石损伤演化方程

2.1. 微元强度分布

在研究地下工程建设的同时，岩体的热损伤问题一直是国内外学者研究的重点，岩石材料的物理力学特性和其高温损伤机理非常重要。高温岩石材料在热损伤作用下，其表现主要为内部的原生裂隙进一步发育，新的微裂纹出现生成、扩展、贯通，最终形成肉眼可见的裂纹趋势，随着温度升高，内部的裂纹越多，内部损伤情况也越严重，其物理力学特性也就发生相应改变。由于岩石内部缺陷的分布演化存在随机性，因此，采用单一、固定的常量往往难以对其内部矿物颗粒的力学性质进行描述。因此本文基于统计强度理论的损伤力学定义应力损伤变量，来描述岩石在加载过程中岩石内部的多尺度无序损伤演化和破坏。定义应力损伤变量D_Q为某一应力水平下已经破坏的微元数目N_F与初始状态下微元总数目N的比值。

$D_Q=\frac{N_F}{N}$ (1)

引入概率密度函数来描述微元体的破坏数目，设岩石微元体的破坏概率服从Weibull分布。

$P\left(F\right)=\frac{w}{F_0}{\left(\frac{F}{F_0}\right)}^{w-1}\exp \left[-{\left(\frac{F}{F_0}\right)}^w\right]$ (2)

式中，F表示Weibull分布变量，w和F₀为Weibull分布的参数。

在任意区间内$\left[F,\text{d}F\right]$ 已破坏的微元数目为$NP\left(F\right)\text{d}F$ ，当加载到某一水平F时，破坏的微元体数目为：

$N_F={\displaystyle {\int }_0^{F}NP\left(F\right)\text{d}F}=N\left\{1-\exp \left[-{\left(\frac{F}{F_0}\right)}^w\right]\right\}$ (3)

由式(1)和(3)得到应力损伤变量为概率密度函数$\varphi \left(F\right)$ 的累计分布函数。

$D_Q={\displaystyle {\int }_0^{F}P\left(x\right)\text{d}x}=1-\exp \left[-{\left(\frac{F}{F_0}\right)}^w\right]$ (4)

2.2. 微元强度确定

利用统计损伤理论建立岩石的统计损伤演化模型，关键在于对岩石微单元强度的合理度量。在岩土领域中，最早出现且使用最广泛的准则之一是Mohr-Coulomb准则。但它忽视了中间主应力对岩石破坏的重要影响。而Drucker-Prager (D-P)屈服准则是M-C准则的一种近似表达。Drucker-Prager准则考虑了静水和中主应力的影响，克服了M-C的不足。D-P准则的屈服面在材料逐渐屈服的过程中保持恒定，不会发生改变，屈服强度会随着侧限压力的增大而相应增强。D-P考虑了由于屈服过程而引起的体积膨胀现象，有助于更准确地预测岩石在屈服状态下的变形行为。本文使用D-P破坏准则表示岩石的微观单元强度准则，引进文献[13]。

$f\left(\sigma \right)=\alpha I_1+\sqrt{J_2}-C$ (5)

$\{\begin{array}{l}\alpha =\frac{2\sin \phi }{\sqrt{9+3{\sin }^2\phi }}\\ C=\frac{C_0\cos \phi }{\sqrt{1+\frac{\sin {\phi }^2}{3}}}\end{array}$ (6)

式中$\alpha$ ，C为材料强度参数可由摩擦角$\phi$ 和C₀确定。I₁为有效应力第一不变量，J₂为有效偏应力第二不变量。

$I_1={\sigma }_1^{*}+{\sigma }_2^{*}+{\sigma }_3^{*}$ (7)

$J_2=\frac{{\left({\sigma }_1^{*}-{\sigma }_2^{*}\right)}^2+{\left({\sigma }_2^{*}-{\sigma }_3^{*}\right)}^2+{\left({\sigma }_3^{*}-{\sigma }_1^{*}\right)}^2}{6}$ (8)

2.3. 损伤变量建立

当岩石材料经历高温处理后，其内部会产生显著的热应力，导致矿物颗粒膨胀并发生挤压错动。这一过程加剧了岩石内部裂隙的扩展与连通，从而显著影响其力学性质。为了更准确地描述这一现象，引入一个与温度紧密相关的函数变量。使用弹性模量作为关键参数，来定义热损伤变量：

$D_T=1-\frac{E_T}{E_0}$ (9)

式中E_T是岩石经过温度T作用后的弹性模量，E₀为岩石常温条件下的弹性模量。

由于岩石在温度与载荷作用下，岩石表现出不同的损伤特性，温度–压力耦合作用下的岩石总损伤不是简单的损伤叠加。文献[14]基于应变等价原理对总损伤变量进行了推导，表达式为：

$D=D_Q+D_T-D_Q\cdot D_T$ (10)

在式(10)右边减去$D_T^2$ 可得下式[15]：

$D=D_Q+D_T-D_Q\cdot D_T-D_T^2$ (11)

整理得：

$D=\left(D_Q+D_T\right)-\frac{1}{2}{\left(D_Q+D_T\right)}^2$ (12)

使用泰勒展开，得到热力耦合下损伤变量为：

$D=1-{\text{e}}^{-\left(D_Q+D_T\right)}$ (13)

联立式(5)、式(13)可得弹性模量在热损伤作用下岩石的损伤演化方程表达式：

$D=1-\exp \left\{\exp \left[-{\left(\frac{\alpha I_1+\sqrt{J_2}-C}{F_0}\right)}^w\right]+\frac{E_T}{E_0}-2\right\}$ (14)

2.4. 热–力作用下的岩石损伤本构模型

基于Lemaitre应变等价性原理和有效应力的概念[16]，三维各项同性损伤的名义应力$\sigma$ 和有效应力${\sigma }^{*}$ 的关系为：

$\left[{\sigma }^{*}\right]=\frac{\left[\sigma \right]}{1-\left[D\right]}=\frac{\left[E\right]\left[\epsilon \right]}{1-\left[D\right]}$ (15)

根据广义胡克定律，有：

${\epsilon }_i=\frac{{\sigma }_i^{*}-\mu \left({\sigma }_j^{*}+{\sigma }_k^{*}\right)}{E}\left(i,j,k=1,2,3\right)$ (16)

根据应变协调原理和应变等价原理，假定各个方向上的应变不变，则：

${\sigma }_i=E{\epsilon }_i\left(1-D\right)+\mu \left({\sigma }_j+{\sigma }_k\right)$ (17)

将式(14)代入式(17)可得本构方程：

${\sigma }_i=E_T{\epsilon }_i\exp \left\{\exp \left[-{\left(\frac{\alpha I_1+\sqrt{J_2}-C}{F_0}\right)}^w\right]+\frac{E_T}{E_0}-2\right\}+\mu \left({\sigma }_j+{\sigma }_k\right)$ (18)

3. 模型参数确定

在岩石三轴试验中，可以得到应变${\epsilon }_1$ 和应力${\sigma }_1,{\sigma }_2,{\sigma }_3$ 根据胡克定律和应变等价原理可知：

$\{\begin{array}{l}{\epsilon }_1=\frac{1}{E}\left({\sigma }_1^{*}-\mu {\sigma }_2^{*}-\mu {\sigma }_3^{*}\right)\\ {\sigma }_i^{*}=\frac{{\sigma }_i}{1-D} \left(i=1,2,3\right)\end{array}$ (19)

则I₁和J₂可表示为：

$\{\begin{array}{l}{I}_1=\frac{E{\epsilon }_1\left({\sigma }_1+{\sigma }_2+{\sigma }_3\right)}{{\sigma }_1-\mu {\sigma }_2-\mu {\sigma }_3}\\ J_2=\frac{{\left(E{\epsilon }_1\right)}^2\left[{\left({\sigma }_1-{\sigma }_3\right)}^2+{\left({\sigma }_1-{\sigma }_2\right)}^2+{\left({\sigma }_2-{\sigma }_3\right)}^2\right]}{6{\left({\sigma }_1-\mu {\sigma }_2-\mu {\sigma }_3\right)}^2}\end{array}$ (20)

则有：

${\sigma }_1=E_T{\epsilon }_1\exp \left\{\exp \left[-{\left(\frac{\alpha I_1+\sqrt{J_2}-C}{F_0}\right)}^w\right]+\frac{E_T}{E_0}-2\right\}+\mu \left({\sigma }_2+{\sigma }_3\right)$ (21)

结合边界条件：

${\epsilon }_1={\epsilon }_p,{\sigma }_1={\sigma }_p$ (22)

${\sigma }_1={\sigma }_p,\frac{\partial {\sigma }_p}{\partial {\epsilon }_p}=0$ (23)

$\frac{\partial \sigma }{\partial \epsilon }=E_T\exp \left[-{\left(\frac{F_{sc}}{F_0}\right)}^w\right]-E_T{\epsilon }_p\exp \left[-{\left(\frac{F_{sc}}{F_0}\right)}^w\right]\left[\frac{w}{F_0}{\left(\frac{F_{sc}}{F_0}\right)}^{w-1}\right]\frac{\partial F_{sc}}{\partial \epsilon }=0$ (24)

${\sigma }_p=E_T{\epsilon }_p\exp \left[-{\left(\frac{F_{sc}}{F_0}\right)}^w\right]+\mu \left({\sigma }_2+{\sigma }_3\right)$ (25)

联立式(24)、式(25)后可得参数w，F₀：

$\{\begin{array}{l}m=\frac{-1}{\ln \left(D_T+1-\ln \frac{E_T{\epsilon }_p}{{\sigma }_p-\mu \left({\sigma }_2+{\sigma }_3\right)}\right)\left(D_T+1-\ln \frac{E_T{\epsilon }_p}{{\sigma }_p-\mu \left({\sigma }_2+{\sigma }_3\right)}\right)}\\ F_0=\frac{F}{{\left(-\ln \left(D_T+1-\ln \frac{E_T{\epsilon }_p}{{\sigma }_p-\mu \left({\sigma }_2+{\sigma }_3\right)}\right)\right)}^{\frac{1}{w}}}\end{array}$ (26)

4. 本构模型验证

4.1. 不同温度、不同围压下红砂岩的理论应力–应变曲线

在对三轴压缩条件下热–力耦合作用下的岩石损伤本构模型进行验证时，本文选取高温后红砂岩三轴压缩试验的相关试验数据进行验证，对25~900℃高温后红砂岩进行围压10 MPa的三轴压缩试验，见图1。试验所得不同温度下红砂岩的三轴压缩相关强度变形参数，根据相关参数带入理论模型，得到了理论应力–应变曲线。

Figure 1. Theoretical stress-strain curves of granite at different temperatures

图1. 不同温度下的红砂岩理论应力–应变曲线

4.2. 结果分析

由图可以看出，理论曲线和实验曲线的峰值应力和峰值应变基本相同，但峰值应力后，温度的升高，导致岩石内部的损伤缓慢累积，使得应变也逐渐增大，引起残余应力不同。损伤变量的指数分解使得温度损伤变量对岩石强度的影响呈现出指数阶的特性较为显著。且随着温度的升高，岩石内部的损伤累积过程发生变化，导致其在破坏后仍能表现出较高的残余承载能力。

5. 总结

本文在分析真三轴条件下岩石变形破坏全过程的基础上，通过对岩石力学变形机理和特征的深入研究，基于损伤理论，根据广义应变理论，假设岩石的微单元强度符合Weibull函数的分布特征，基于D-P准则并对损伤变量进行分解建立三维应力作用下的岩石损伤模型，主要结论如下。

1) 对比三轴压缩条件下理论和实验的峰值应力、峰值应变，发现拟合效果较好，说明构造的统计损伤本构模型能够反映红砂岩热损伤后的应力–应变峰值特性，这验证了模型的合理性。

2) 拟合曲线能够较好模拟实验曲线的走势，且损伤变量的指数分解使得拟合曲线的延性更强。