考虑初始压密阶段的千枚岩变形破坏全过程统计损伤本构模型

李占甫，莘子健 ，付聪，李玉豪，张家铭

(1.安徽省交控建设管理有限公司，安徽 合肥 230088；2.中国地质大学（武汉）工程学院，湖北 武汉 430074)

1 引言

岩石是由具有一定结构构造的矿物集合体组成的非均质材料，因其生成条件、矿物成分、赋存环境等不同，使得岩石内部存在大量随机分布的微观裂隙。岩石的宏观破坏可看作其内部微观裂隙的扩展和累积，该过程称之为损伤[1]。损伤理论以材料内部微观裂隙的萌生、扩展和相互作用为基础，研究试件由微观直至发生宏观破坏的全过程，是目前研究岩石等含有天然微裂隙材料的有效方法[2]。基于岩石内部微观裂隙的分布及演化具有随机性，Krajcinovic等人[3]将统计强度理论与连续损伤理论相结合，开辟了一条能较好反映岩石损伤破坏过程和应力—应变关系的本构模型建立途径，即统计损伤本构模型。

自统计损伤理论提出以来，国内外学者利用其对岩石变形破坏过程的模拟方法进行了大量探索，并取得较为丰硕的研究成果。曹文贵[4]先基于Lemaitre应变等价性假设[5]定义了传统的岩石损伤变量，石崇[6]、曹瑞琅[7]引入不同修正系数对其修正，但因该定义忽略了岩石破坏后的残余强度，与实际不符，随后曹文贵[8]将岩石骨架部分抽象为损伤与未损伤两部分，对岩石的损伤重新定义，并建立可以反映岩石残余强度的损伤模型，尹杰[2]在此基础上基于不同准则进行计算和修正。周永强[9]和刘冬桥[10]分别考虑损伤阈值和微元强度的度量方式建立统计损伤模型。尽管当前理论模型能够较好描述岩石部分阶段的变形特征，但对于岩石变形破坏全过程的模拟仍存在一定不足，特别是在反映岩石初始压密阶段。究其原因，当前沿用较广的统计损伤模型忽略了岩石内部在初始状态下存在的大量初始缺陷，并在模型推导过程中将岩石弹性模量视为常数。事实上岩石在受力过程中，因内部微结构面的存在其岩石宏观变形与力学特性会产生较大变化，如千枚岩，因其特殊的千枚状构造，其微结构面发育，该类岩石在初始压密阶段的应力应变关系会呈现较强的非线性变形特征。此外，大量试验表明，岩石的弹性模量与初始缺陷的变形及闭合情况密切相关，并在一定范围内随所受围压水平的变化而改变，将其视为常数并不恰当。

为此，曹文贵[11]将岩石抽象为空隙与骨架两部分，采用宏观与微观相结合的方法，建立了可以反映脆性岩石初始宏观非线性变形力学行为的模拟方法，但其模型在模拟煤岩在围压为8MPa以下时拟合效果不理想；张超[12]、李修磊[13]在文献[11]的基础上对模型进行了改进，虽较好地拟合了砂岩在围压40MPa下的应力应变关系，但该岩石变形破坏过程中初始压密阶段并不明显，不能很好地验证该类模型对于反映含大量微结构面岩石初始压密阶段的变形特征。此外，前述三位学者均采用Mohr-Coulomb准则建立相关模型，该准则在描述岩石处于低应力状态下会存在一定偏差[2]。

基于上述分析，尽管现有岩石统计损伤模型对于岩石部分变形阶段特征可以较好描述，但能够反映考虑岩石初始压密阶段的全过程变形特征的本构模型较少。此外，多数模型所针对的岩石为砂岩、花岗岩以及粉砂质泥岩等，而对于含有大量初始微结构面的千枚岩研究鲜有报道。针对上述不足，本文在前人的研究基础上，基于文献[11]提出的考虑岩石初始宏观非线性变形力学行为的模拟方法，选用Drucker-Prager准则针对千枚岩在较低围压状态下的全过程变形特征进行模拟，以期建立出能反映岩石全过程变形特征的统计损伤模型。

2 岩石变形过程

岩石受外部荷载作用下的变形可划分为岩石骨架和内部微结构面两部分，且岩石骨架在屈服破坏前的变形为线弹性。低荷载条件下岩石的骨架变形与微裂隙闭合同时发生，应力应变曲线呈明显非线性；当荷载增加到一定水平时，微裂隙闭合完成，此时随荷载的增加至岩石屈服破坏前，岩石仅发生线弹性的骨架变形；荷载进一步增加达到岩石屈服强度后，岩石的骨架变形呈非线性，并依次出现屈服、应变软化和完全破坏的现象。因此，岩石在受荷变形直至破坏过程中具有阶段性变形特征：初始压密(OA)、线弹性变形(AB)、屈服(BC)、应变软化(CD)和完全破坏阶段（D点之后），如图1所示。

图1 千枚岩变形破坏全过程

3 岩石变形分析模型

为分析岩石在初始压密阶段内部微裂隙产生的变形，在岩石内取一代表性柱体单元，该柱体单元由岩石骨架和微裂隙两部分组成，如图2所示。由前述岩石变形机理可知，在初始压密阶段微裂隙变形与骨架变形并不协调，故设柱体单元在加载前初始高度为h0，其中微裂隙部分高度为，骨架部分高度为该柱体单元在某荷载应力σi下总变形量为Δh，岩石内部微裂隙和骨架部分变形量为Δhw和Δhr，则有如下关系：

图2 千枚岩变形分析模型

故岩石在荷载应力σi的作用下的总应变（宏观应变）εi、微裂隙部分应变和骨架部分应变可表示为：

为建立岩石总应变与其组成部分应变间的关系，特定义柱体单元微裂隙部分占比为b0,即利用式（1）～（5）可得：

分析式（6）可看出岩石在初始压密阶段（b0≠0）的宏观应变由微裂隙与骨架两部分组成，且呈非线性，当微裂隙完全闭合（b0=0）后，岩石宏观应变仅由骨架部分提供。由岩石变形破坏过程知，在岩石未达到屈服强度前骨架部分材料的变形线弹性状态，进一步证明在初始压密阶段的非线性变形是由微裂隙闭合导致的，显然式（6）作为岩石变形分析模型是符合实际的。

3.1 微裂隙部分变形分析方法

按前述分析思路，若想获得岩石微裂隙在荷载应力σi作用下的应变需先求得相应微裂隙部分（hw）的变形量Δhw，故将荷载σi划分为n个等级组成的应力增量并逐级施加，故有：

岩石微裂隙部分的闭合变形具有不可恢复的特点，采用真应变能够更准确地描述其变形特点，故将变形量Δhw可看作在某一级应力增量作用下变形量的累加，设微裂隙部分（hw）在应力增量作用后的高度变为采用文献[14]方法可得此时微裂隙部分的应变可表示为：

由于微裂隙部分材料的累计变形高度无法通过试验测量，则假定服从广义胡克定律，微裂隙部分材料在第s级应力增量作用下产生的应变增量可表示为：

联立式（8）和式（9）并变形有：

由式（11）可知，若想得到空隙部分在荷载应力σi下的应变得先求出在各级荷载作用下相应的变形模量和泊松比，但常规试验无法获得，为便于研究，假设微裂隙部分材料的力学参数不随岩石变形而改变，从而式（11）可改写为：

式中E1和μ1为岩石微裂隙部分材料的弹性模量和泊松比，相关参数的确定方法详见第4.2。

3.2 岩石骨架部分变形分析方法

岩石骨架部分的变形可根据荷载是否达到屈服强度分为两部分，未达到屈服强度时呈线弹性变形，超过屈服强度岩石骨架发生损伤，致使变形呈非线性，并依次出现应变硬化、应变软化和完全破坏的现象，对此可采用基于连续介质的统计损伤理论进行分析。目前常用的统计损伤模型基于Lemaitre应变等效假设，并不能准确反映岩石骨架部分损伤的力学本质，尤其是认为岩石完全破坏后不具备任何承载能力，即忽略了岩石的残余强度，为此，曹文贵等[8]对此类模型进行修正，以解决该类本构模型的不足，即：

式中，σi为名义应力为有效应力为岩石的轴向残余强度；D为损伤变量，是材料损伤程度的度量。

将岩石骨架部分抽象为未损伤和损伤两部分，骨架整体受力为σi，其作用面积为A，则未损伤材料面积为A1，损伤部分面积为A2，两部分所受微观应力分别为如图3所示。则损伤变量D可定义为损伤部分面积与岩石骨架总面积之比，即D=A2/(A1+A2)。

图3 千枚岩骨架部分变形分析模型

基于Lemaitre应变等效假设，可认为岩石骨架部分的宏观应变与未损伤部分微观应变相等。为骨架未损伤部分变形的微观应力，且骨架未损伤部分为线弹性体，服从广义胡克定律，则有：

式中，E2和μ2为岩石骨架部分弹性模量和泊松比为骨架未损伤部分所受有效应力。假设骨架部分材料的损伤仅沿轴向发生，即在侧向方向上名义应力等于有效应力，则有：

由式（13）、（14）及上述假设可得：

故岩石骨架部分变形可表示为：

4 岩石统计损伤本构模型及参数确定方法

4.1 岩石统计损伤本构模型

对于岩石骨架部分的变形破坏为一连续损伤过程，可采用统计损伤理论进行模拟，模拟关键在于对岩石骨架部分微元强度的合理度量，经对比不同强度准则的优劣，本文选用Drucker-Prager准则进行相关计算。该准则表达式为：

式中F为岩石骨架部分微元强度，取值分两种情况：当骨架部分达到损伤阈值时，F≥0未达到时，变形呈线弹性，F＜0；I1和J2分别为应力张量的第一不变量和应力偏量的第二不变量；a和k分别为与岩石骨架部分粘聚力（c）、内摩擦角（φ）相关的常数。上述参量可表示为：

由式（14）代入式（18）和（19），并在常规三轴压缩试验条件下（σ2=σ3），可得：

为得到损伤变量D，采用可以较好表征骨架材料微元强度变形破坏随机性的Weibull概率密度函数，则损伤变量D与微元强度F间的关系可表示为：

式中m与F0为岩石骨架部分微元强度Weibull分布参数，当岩石处于完全破坏阶段时D=1，岩石未发生损伤时D=0。将式（12）、（13）、（16）代入式（6）可得到考虑岩石内部微裂隙及残余强度的岩石变形破坏全过程模型，即：

4.2 模型参数确定方法

①b0的确定方法

由前述岩石变形破坏过程可知，岩石宏观变形由微裂隙闭合与骨架变形两部分组成，当微裂隙完全闭合时且未达到屈服强度前，骨架部分材料处于线弹性状态，此时E1趋近于0，式（25）中第1式中可变为：

若σ3=0时，即岩石处于单轴压缩状态，式（26）可变为：

式（27）经转换可得：

由式（28）可知，b0为岩石单轴压缩试验曲线线弹性段的延长线在应变轴上的截距。

②E2和μ2的确定方法

μ2为岩石骨架材料的泊松比，故可通过对已完成微裂隙压密阶段的岩石进行泊松比试验直接测定。E2为岩石骨架部分弹性模量，可通过式（26）变换并在三轴试验中弹性变形阶段取值代入求得，即：

③E1和μ1的确定方法

在已知岩石在初始压密阶段满足式（25）中第1式和已求得模型参数b0、E2和μ2的基础上，在初始压密阶段任意取值，采用基于最小二乘法原理的曲线拟合方法可得E1和μ1。

④微元强度随机分布参数的确定

微元强度随机分布参数m与F0为模型建立的关键，由岩石具有应变软化特性可知，其应力应变关系曲线在峰值点（σ1=σsc，ε1=εsc）处具有极值特性，即：

式中σsc和εsc分别为岩石三轴试验中峰值点的应力和应变。

同时式（30）满足式（25）中第2式，联立方程即可确定参数m与F0：

至此，已建立考虑岩石内部微裂隙及残余强度的岩石变形破坏全过程模拟方法，并给出相关参数的确定方法，现对其合理性进行验证分析。

5 实例分析与验证

引入文献[15]的试验曲线对前述模型进行验证分析。利用该文献针对千枚岩在结构面夹角为0°时，围压分别为0、2、4、6和8MPa下的岩石三轴试验曲线，如图4所示。由该文献可知黏聚力C=11.66MPa，内摩擦角φ=50.35°，按照前述参数确定方法可知：b0=0.0017，不同围压下相关计算参数见表1所示。将相关计算参数代入式（31）～式（34）可求得微元强度随机分布参数m与F0。将试验数据及各计算参数代入本文所提出的本构模型，得到不同围压作用下千枚岩全过程变形应力应变理论曲线，并将其与试验曲线相比较，如图5所示。

模型参数 表1

图4 千枚岩三轴试验曲线

图5 本文理论曲线与试验曲线的比较

为进一步验证本文理论模型和方法的合理性与可行性，分别根据文献[11]、文献[12]和文献[13]所提出的本构模型，以文献[15]围压4MPa试验数据为基础，分别获取各模型理论曲线，并与本文模型所获理论曲线相对比，如图6所示。

图6 围压4MPa下不同模型理论曲线与试验曲线的比较

由图5可以看出，本文所提出的统计损伤本构模型，针对千枚岩在不同围压条件下的计算结果与试验曲线均有较高的吻合程度，很好地反映千枚岩变形破坏全过程应力应变特征。以围压4MPa试验数据为例，比较已有模型和本文模型对其模拟效果，可以看出，由于文献[11]所建立模型并未考虑岩石破坏后的残余强度，导致其残余强度变形阶段与试验曲线差别较大。此外，文献[11]、文献[12]和文献[13]所建立模型均基于Mohr-Coulomb准则，致使其在低围压状态下与试验曲线差别较大，特别是在岩石屈服和应变软化阶段。由此可见，本文所构建模型能够较好地模拟含大量微结构面的千枚岩在低围压状态下的全过程变形破坏特征。

6 结语

本文利用统计损伤理论，对能反映岩石在初始压密、屈服破坏及残余变形破坏等阶段的岩石变形全过程模拟方法进行了研究。得出如下结论：

①将岩石抽象为骨架和内部空隙两部分，分别建立了空隙和骨架部分变形分析方法，进而建立了可模拟空隙岩石变形破坏全过程的统计损伤本构模型；

②由于Mohr-Coulomb准则，并不能很好描述岩石处于低应力状态下的变形破坏特征，本文选用Drucker-Prager准则，并以含有大量初始微结构面的千枚岩为对象，对其在低围压状态下的全过程变形特征进行模拟研究；

③通过本文模型、现有相关模型和试验曲线的比较，验证了本文所建岩石损伤模型的有效性和合理性。