考虑损伤修正的岩石统计损伤本构模型研究

杨宇生，尹前锋，丰丛杰，贾洪彪

(1.中国地质大学工程学院,湖北 武汉 430074; 2.中国科学院大学地球与行星科学学院,北京 100049)

岩石是一种复杂的天然地质材料，由于岩石内部常包含大量随机分布的微裂隙等缺陷，使得岩石材料具有明显的各向异性和渐进破坏性[1]。近年来，基于连续损伤理论和统计强度理论的岩石损伤本构模型研究取得了很大的进展，成为研究岩石本构关系的主流。

文献[2]首次将连续损伤理论与统计强度理论结合，提出了损伤模型，开启了在岩石损伤本构模型方向上的探究之旅。文献[3]在前人的基础上，提出了将轴向应变作为统计分布变量的观点，对岩石损伤破裂过程进行了研究，为后续的研究奠定了基础。文献[4]首次提出岩石微元强度的概念，并以岩石微元强度服从Weibull分布为出发点，基于屈服准则建立了岩石损伤本构模型。文献[5]提出基于正态分布的岩石微元强度，以此建立了新的岩石损伤本构模型。文献[6]在前人研究的基础上提出了新的衰减型函数——Harris函数，作为岩石微元强度概率密度分布函数，并以此建立新的本构模型，较好的反应了岩石在三维应力状态下的应力——应变关系和破坏全过程。本文在现有的研究基础上，引入改进的Harris函数作为岩石微元强度概率密度分布函数，并对损伤变量进行修正，引入损伤变量修正系数，采用较为简洁的定义岩石微元强度随机分布变量的方式，提出了新的岩石损伤统计本构模型，运用MATLAB软件将实验数据点与本构模型进行拟合，求取相关模型参数，最终确定了岩石损伤本构模型。

1 岩石损伤本构模型的建立

假定理想条件下岩石内部不存在缺陷，为无损状态，在三维应力作用下其本构关系服从胡克定律，则有

(1)

式中：σ1、σ2、σ3为岩石所承受的3个主应力；ε1、ε2、ε3为岩石在3个主应力方向上所产生的主应变；G为岩石剪切模量，λ为拉梅常数。

在等围压条件下，无损岩石材料三轴压缩轴向应力——应变关系根据式(1)可得

σ1=Eε1+2νσ3

(2)

式中：E为岩石弹性模量；ν为泊松比。

在天然状态下，岩石内部往往分布有各种各样的微裂隙和结构面，这些均为岩石材料的初始缺陷，导致岩石材料的物理力学性质的非均质性[7]，为此我们引入损伤变量D，设在某荷载作用下已破坏的岩石微元体数目为c，岩石微元体总数为N，定义岩石损伤变量为D=C/N[8]，当D=0时，说明岩石为无损状态，D=1时说明岩石为完全损伤状态，D的取值反映了岩石材料内部的受损程度。在单一受力条件下，由连续损伤理论[9]，可得

σ=Eε(1-D)

(3)

由此可推导在三维应力作用下受损岩石的本构关系为

σ1=Eε1(1-D)+2νσ3

(4)

对于岩石材料而言，其单轴压缩曲线和三轴压缩曲线在累进性破裂阶段达到峰值点后，内部结构被完全破坏，裂隙迅速发展，岩石的承载力随变形的继续增大而下降，但并不会完全消失，即岩石的残余强度。而根据损伤参量的定义，当达到峰值点时岩石完全破坏，此时D=1，代入到式(3)中得出应力为零，此时表示岩石无残余强度，与实际情况不符。

为此，引入损伤变量修正系数δ∈(0，1)，通过与损伤变量D相乘的方式对其进行修正，使得到的岩石损伤统计本构模型能反映岩石材料的残余强度，以及岩石材料各向异性的特点，使得到的模型更接近实际情况，更好的模拟出岩石在达到峰值点后的应力——应变特征。通过分析，在引入损伤变量修正系数δ后得到的新的岩石损伤统计本构模型为

σ1=Eε1(1-δD)+2νσ3

(5)

根据前文所提，目前大多数研究中均假设岩石微元强度服从Weibull分布和正态分布，但两者均有其缺点和不足，为此本文引入著名的衰减型函数模型Harris函数[10]，其表达式为

(6)

式中：F为自变量，S为因变量，a、b为模型参数，均大于0。

在Harris函数中，S随F的递增而减小，其单调性与D相反，且S与D的变化范围均为[0，1]。为此，本文对Harris函数进行修正，使之作为岩石损伤变量的发展方程，即

(7)

式中：将F重新定义为岩石微元强度分布变量。

将式(7)两边对F求导，即得到岩石微元强度概率密度分布函数Φ(F)

(8)

目前对岩石损伤本构模型的研究，均引入不同的岩石强度准(Mohr-Coulomb准则[11-12]、Hoek-Brown准则[13-14]、Druckre-Prager准则[15]等)作为岩石微元强度随机分布变量，通过研究发现，若以岩石屈服准则作为微元强度随机分布变量，推导出的表达式较为复杂，而以岩石轴向应变作为微元强度随机分布变量则简单许多，因此，本文同样采取此种简化方法，将岩石轴向应变作为微元强度随机分布变量，此时F=ε1。

将式(7)代入式(8)可得

(9)

2 模型的验证

为验证所建岩石损伤统计本构模型的合理性，采用INSTRON伺服机对微风化石英片岩进行了单轴压缩和三轴压缩实验(见图1～图3)，得到石英片岩的单轴抗压强度以及在2MPa、4MPa和6MPa三种不同围压下的全应力——应变曲线，得到了岩石弹性模量E及泊松比ν，其常规力学参数及不同围压下的特征参量如表1～表2所示。

图1 INSTRON伺服压力机及数据采集系统

图2 三轴试验所用压力室

图3 岩样实验前后照片

岩石种类泊松比ν弹性模量E/GPa单轴抗压强度σc/MPa内摩擦角φ/(°)石英片岩0.23521.3399.3955.276

表2 不同围压下的特征参量

通过获得的实验数据，运用曲线拟合法，利用MATLAB软件进行曲线拟合，得到模型参数a、b和修正系数δ的值，最终确定岩石损伤统计本构模型，并将模型曲线与实验数据点进行拟合，得到如图3所示的不同围压下模型曲线与实验数据点的吻合对比图，通过对比模型曲线与实验数据点的吻合情况可以发现，本文推导出的岩石统计损伤本构模型曲线能够很好的反映岩石应力应变全过程，实验数据点和模型曲线吻合良好，不仅在岩石峰前阶段吻合度良好，由于引入了修正系数δ，也很好的反映了岩石峰后破坏阶段。由实验可知，随着所加围压的增加，岩石的脆性降低，相对应的延性提高，岩石的残余强度表现的更加明显，本文推导出的岩石统计损伤本构模型对岩石峰后阶段有着很高的拟合精度，较传统的线性拟合方法更加优越，在某种程度也上体现了岩石的各向异性和渐进破坏性，符合已知的岩石材料的实际情况，具有一定的实用价值。

图4 不同围压下模型曲线与实验数据点拟合情况

3 相关参数的探讨

本文推导的本构模型中创新的引入了岩石损伤修正系数δ，在此来探讨δ对本构模型的影响，如图4所示，引用前文模型验证中INSTRON伺服机对微风化石英片岩进行的在4MPa围压下的三轴压缩实验数据，考虑修正系数δ的不同取值对模型曲线的影响。从图5中可以看出，当δ取1时，即不引入岩石损伤修正系数时，模型曲线显然与实际相差很大，很难反映出岩石破坏后的应力应变情况，说明修正系数δ的引入很好的表现了微风化石英片岩峰值点后的破坏情况，同时也反映了它的残余强度，优化了以往研究中本构模型的不足，具有很好的实际意义。从图5中还可以看出，修正系数δ的不同取值，对模型曲线的前半段影响不大，即石英片岩峰值点前影响不大，对曲线峰值点后影响较大，在实际工程应用中，δ的取值会随着岩性、围压以及初始缺陷的不同而变化，所以δ的准确取值是至关重要的，运用MATLAB软件进行曲线拟合确定δ的取值不失为一种好方法。

围压2MPa图5 修正系数δ的取值对模型曲线的影响

4 结论

本文提出了一种新的岩石统计损伤本构模型，以微风化石英片岩为例进行了实验论证，主要创新点和所得到的结论如下：

(1)根据以往的研究，提出运用改进的Harris 函数作为岩石微元强度概率密度分布函数，基于连续损伤理论和统计强度理论建立了岩石在三维应力作用下的统计损伤本构模型，同时在求取模型参数时，摒弃了以往研究中的多元函数求极值法，运用MATLAB软件，采用精度更高的曲线拟合法确定了模型参数。

(2) 由于对拟合精度产生的影响十分微小，本文将岩石轴向应变作为微元强度随机分布变量来建立本构模型，简化了模型，方便了运算。

(3)创新的引入岩石损伤修正系数δ，进一步提高了模型曲线在岩石峰值后阶段的拟合精度，很好的反映了岩石在三维应力状态下的渐进破坏性，通过相同受力条件下δ的不同取值对模型曲线与实验数据拟合精度的影响，表明引入δ的必要性，为理论研究和实际工程提供参考，具有很好的现实意义。