0 引言

岩石是一种天然地质材料, 由固体骨架、节理裂隙或孔隙及各类填充物组成, 赋存于一定地质环境(包括地应力、地温、地下水以及化学环境等) 中, 具有多组元、非均质、跨尺度的复杂结构^[1]。三轴压缩条件下, 岩石除线弹性和残余强度变形阶段外, 主要表现为非线性变形, 包括空隙压密、屈服硬化和应变软化等阶段。为了能够模拟出岩石变形破坏过程, 很多学者基于不同的力学理论如固体力学、损伤力学等, 建立了不同的本构模型如弹塑性本构模型^[2-3]、损伤本构模型^[4-5]等, 但其仅能模拟岩石变形阶段中的一个或几个, 无法模拟岩石变形破坏全过程, 也不能反映弹性模量随围压变化而变化, 其根本原因在于岩石屈服破坏前后的非线性变形机理不同, 亟需一个合理可行的非线性分析模型将两者联系起来, 并在此基础上提出相应的计算方法。

岩石内部含有随机分布的初始空隙, 当空隙率较小时, 空隙对岩石变形的影响几乎可以忽略, 可视为连续固体材料, 认为其在屈服破坏之前服从广义胡克定律, 弹性模量不随围压变化而变化。然而, 当空隙率较大时, 岩石在外荷载作用初期表现出明显的非线性压密变形过程, 弹性模量也会随围压变化而变化, 则必须考虑空隙对岩石变形的影响。然而, 很多学者在运用损伤力学和强度统计理论建立岩石变形破坏过程模拟方法^[6-7]时, 忽视或不考虑空隙压密非线性变形阶段以及弹性模量的变化, 因而该类模拟方法仅适用于空隙率较低的岩石。为此, 部分学者基于岩石组构特征将岩石类材料视为由空隙部分和骨架部分构成, 以空隙率来反映空隙对岩石变形的影响, 提出了新型岩石变形破坏过程模拟方法^[8-9], 但没有考虑空隙形态(孔隙和裂隙) 对岩石变形影响的差异性, 以及空隙压密引起岩石变形性质的不同, 导致该类方法虽能反映空隙压密引起的岩石非线性变形, 但模拟效果仍然较差。于是, 有学者将岩石初始压密变形区分为由可恢复变形和不可恢复变形2部分组成^[10], 基于真实应力-真实应变分析方法建立不可恢复应变与应力之间的关系, 建立能够模拟三轴压缩条件下岩石偏应力-应变关系曲线的本构模型, 但该模型仍不能反映岩石弹性模量变化特征, 也未反映出岩石残余强度变形阶段特征。由此可以看出, 目前研究尚未探究岩石屈服破坏前后的非线性变形机理以及空隙对岩石变形产生的影响, 虽能反映空隙岩石变形的部分特征, 但始终无法在同一模型中既能反映空隙岩石变形破坏全过程, 又能反映空隙岩石弹性模量变化特征, 因此, 仍须在此基础上继续对空隙岩石变形全过程的模拟方法进行研究。

理想的空隙岩石变形破坏过程模拟方法不仅能反映应变软硬化和残余强度等特性, 还能反映空隙压密阶段的非线性变形特性以及弹性模量随围压发生变化的特性。对此, 本研究将在探讨空隙岩石屈服破坏前后非线性变形机理与特征分析的基础上, 考虑空隙对岩石变形产生的影响, 基于细观唯象分析方法提出空隙岩石非线性变形细观唯象分析模型, 并在此基础上汲取现有岩石统计损伤模拟方法的优越性, 对不同围压作用下空隙岩石变形破坏全过程的模拟方法进行研究, 以期完善岩石损伤本构模型的研究内容与方法。

1 空隙岩石非线性变形机理与特征分析

空隙岩石非线性变形机理与特征分析是建立空隙岩石非线性变形细观唯象分析模型的基础, 因此, 根据相关研究成果^[11-18]对其进一步研究并总结如下:

(1) 岩石内部空隙形态各异, 大致可分为孔隙和裂隙2种, 其对岩石变形力学性质的影响存在显著的差异性。一般而言, 仅含孔隙的岩石在应力作用下的变形关系服从广义胡克定律, 其在屈服破坏前仅产生可恢复变形而不会产生不可恢复变形, 因此, 可将仅含孔隙的岩石等同于连续固体材料来研究屈服破坏前后的应力应变关系。仅含裂隙的岩石在应力作用下裂隙将发生闭合, 使空隙岩石在屈服破坏前产生不可恢复变形, 该变形的持续产生必须伴随持续增加的应力作用, 因此其应力应变关系的建立需基于瞬时应变分析方法。

(2) 裂隙闭合引起的不可恢复变形完成后, 空隙岩石在应力持续作用下的变形服从广义胡克定律, 但会在屈服破坏前产生因裂隙闭合而引起的可恢复变形, 虽然此时空隙岩石变形服从广义胡克定律, 但还会因裂隙闭合而引起弹性模量发生变化。当所受应力超过屈服强度后, 空隙岩石开始发生损伤破坏, 其内部微裂隙开始不断孕育、扩展并相互贯通, 产生因岩石屈服破坏而引起的非线性变形, 包括应变硬化和应变软化。当空隙岩石完全损伤破坏后, 其只能承受残余应力即残余强度, 此时空隙岩石在几乎不变的残余应力作用下产生持续增长的不可恢复变形。

因此, 概化后特定围压作用下的空隙岩石典型偏应力-轴应变关系曲线具有显著的阶段性特征, 见图 1, 概括分析如下:

(1) 空隙压密非线性上凹变形阶段(OA段): 基于三轴压缩试验规程可知, 荷载作用初期已避免了试验机与试件端部接触压密的影响, 因此该阶段主要为空隙受到压密的表现, 其变形包括不可恢复变形和可恢复弹性变形, 其中, 不可恢复变形是因裂隙闭合引起, 当到达A点后, 不可恢复变形已基本完成。

(2) 线弹性变形阶段(AB段): 低应力水平作用下的空隙岩石变形仅服从广义胡克定律, 但由于空隙岩石内部裂隙还会发生可恢复的闭合变形, 从而导致不同围压水平下岩石弹性模量发生变化。

(3) 应变硬化和应变软化变形阶段(BC段和CD段): 当所受应力达到屈服强度后, 空隙岩石内部新裂隙开始孕育, 若要使其继续变形, 需相应地增加应力, 即表现为应变硬化特征。当所受应力达到峰值强度后, 空隙岩石已孕育产生的大量裂隙开始迅速扩展、聚集, 甚至贯通形成宏观破裂面, 与此同时, 岩石承载能力迅速丧失, 直至稳定在残余强度数值附近, 表现为应变软化特性。

(4) 残余强度变形特征阶段(D点以后): 当所受应力下降接近至残余强度时, 裂隙扩展变缓并逐渐趋于停止, 岩石所能承受的应力为摩擦力提供, 基本保持不变, 而变形却持续增加。

值得注意的是, 尽管不同围压作用下空隙岩石变形破坏全过程表现出相似的阶段性特征, 但还会因围压水平不同而存在显著差异, 见图 2。不同围压作用下的空隙岩石应力特征值即空隙岩石裂隙闭合应力(A点)、屈服应力(B点)、峰值应力(C点) 和残余应力(D点) 随围压变化而变化。

基于常规三轴压缩试验规程可知, 特定围压或静水压力作用下空隙岩石变形完成后, 再受到轴向应力的作用, 并将岩石由静水压力状态向偏应力状态的转折处作为三轴压缩试验曲线的记录始点, 对于线弹性变形阶段, 相同围压下空隙岩石弹性模量在轴向应力增大过程中保持不变, 不同围压下空隙岩石弹性模量在轴向应力增大过程中发生变化, 即随着围压的增大, 弹性模量不断增大, 当围压增至一定程度时即空隙岩石内部裂隙完全闭合, 此时弹性模量既不因围压变化而变化, 也不因轴向应力的增大而变化, 而是基本保持为常量, 由此表明空隙岩石弹性模量主要取决于围压水平即三轴压缩试验初期所施加的静水压力。

2 空隙岩石非线性变形细观唯象分析模型

应力持续作用下的空隙岩石必然产生不可恢复和可恢复2部分变形, 前者变形由空隙永久性闭合和空隙岩石损伤破坏过程引起, 后者变形由空隙暂时性闭合和空隙岩石弹性力学行为引起。

岩石属于多缺陷性介质, 除孔隙外, 在地质作用过程中更是发育各种成因的裂隙, 导致岩石的空隙性较为复杂^[19], 其在应力作用下可能发生闭合而引起空隙岩石变形, 但是, 由于裂隙和孔隙的形态不同, 其闭合特点存在较大差别, 有的裂隙和孔隙在应力作用下发生永久闭合即应力卸除后仍处于闭合状态, 使空隙岩石产生永久不可恢复变形; 有的裂隙和孔隙在应力作用下虽能够完全或部分闭合, 但应力卸除后能够恢复为张开状态, 使空隙岩石产生可恢复变形。因此, 本研究在空隙岩石非线性变形机理与特征分析的基础上, 采用细观唯象分析方法将空隙岩石细观化为2种材料, 一种材料仅产生不可恢复变形, 而另一种材料既能产生可恢复变形又能产生因其屈服破坏引起的不可恢复变形。于是, 将这2种材料进行串联以构成空隙岩石非线性变形细观唯象分析模型, 根据串联组合体的应力应变特性可知, 这2种材料产生的变形之和等于空隙岩石产生的变形。

由于空隙闭合产生的不可恢复变形主要是因裂隙永久闭合引起的, 当然其中少部分也可能是由孔隙永久闭合引起的, 如图 3所示, 本研究将上述因空隙永久闭合产生不可恢复变形的部分材料定义为裂隙材料, 而将另一部分材料定义为基质材料, 其仍然包含裂隙和孔隙, 只是其中裂隙在其屈服破坏前仅发生可恢复的闭合变形, 并引起基质材料弹性模量的变化, 而孔隙体积与基质材料弹性模量两者并无直接关系。于是, 对裂隙材料和基质材料的变形关系进行分析, 建立空隙岩石非线性变形分析模型, 并在此基础上分别探讨2部分材料的变形计算方法, 这就是本研究中空隙岩石非线性变形细观唯象分析模型研究的思路。设空隙岩石代表性柱状单元长度为l, 其中, 基质材料累计长度为l^r, 裂隙材料累计长度为l^p, 在应力作用下, 空隙岩石代表性柱状单元的变形量为Δl, 基质材料变形量为Δl^r, 裂隙材料变形量为Δl^p, 则:

(1)

空隙岩石沿代表性柱状单元长度i方向的宏观应变ε_i可表示为:

(2)

式中, 和分别为基质材料和裂隙材料沿代表性柱状单元长度i方向的应变; γ为裂隙材料累计长度在空隙岩石代表性柱状单元中的占比, 可分别表示为:

(3)

(4)

(5)

由于为裂隙材料的应变, 是空隙岩石在屈服破坏前产生的不可恢复非线性变形, 因而根据真实应力-真实应变分析方法可得出裂隙材料变形计算方法, 即:

(6)

式中, E_p和μ_p分别为裂隙材料的弹性模量和泊松比; σ_i, σ_j和σ_k为空隙岩石的3个主应力分量; (i, j, k) = (1, 2, 3) 或(2, 3, 1) 或(3, 1, 2)。由此可以看出, 若要建立空隙岩石变形破坏全过程模拟方法, 关键在于解决基质材料的分析方法。

3 基质材料变形计算方法

应力作用下基质材料既产生屈服破坏前可恢复变形, 也产生屈服破坏后不可恢复变形, 可基于岩石变形统计损伤分析方法^[20-21]进行计算。然而, 空隙岩石在屈服破坏前不可恢复变形完成后, 空隙岩石实际上即可视为基质材料, 其中仍存在未发生闭合的裂隙, 其在应力作用下还会继续发生闭合, 从而引起空隙岩石弹性模量即基质材料弹性模量发生变化, 因此, 在计算基质材料变形之前, 需要首先建立基质材料弹性模量分析模型。

3.1 基质材料弹性模量分析模型

假定基质材料内部裂隙均为硬币形椭球裂隙, 其半轴长在x, y, z方向分别等于a, a和λa, 其中, a和λ为裂隙的半长度和纵横比(即裂隙最小直径与最大直径之比)。所有硬币形椭球裂隙的半长度数值保持不变, 但其纵横比各不相同^{[16, 19]}, 则:

(7)

式中, β和V分别为基质材料的体积柔度和体积; N为基质材料内部张开裂隙总数; β_s和μ_s分别为岩石骨架的体积柔度和泊松比。由材料力学可知, 基质材料体积模量K = E/[3(1-2μ)] = 1/β, 岩石骨架体积模量K_s = E_s/[3(1-2μ_s)] = 1/β_s, 其中, E和μ为基质材料的弹性模量和泊松比, E_s为岩石骨架的弹性模量。对于空隙率较大的岩石, 其基质材料和骨架材料的泊松比存在差异, 但目前其参数之间定量关系的可靠性尚未得到充分论证, 而对于空隙率较低的岩石, 可近似认为其基质材料和骨架材料的泊松比相等。因此, 为便于开展后续研究, 本研究暂忽略岩石材料泊松比发生变化的影响(μ =μ_s), 可得:

(8)

由于应力作用下基质材料内部初始裂隙的长度不发生变化^[16], 因此, 如果忽略应力作用下基质材料体积变化, 即认为V近似为常数, 则基质材料弹性模量仅与其未闭合的裂隙数直接相关。为了建立应力作用下基质材料弹性模量分析模型, 就必须探讨裂隙数与应力之间的关系。静水压力p作用下空隙岩石内部硬币形椭球裂隙发生闭合而引起弹性模量增大^{[16, 19]}, 即:

(9)

式中, λ₀和λ分别为静水压力作用前后裂隙纵横比。于是, 将p视为围压σ₃, 则:

(10)

如果裂隙发生闭合(即λ = 0), 由式(10) 可得裂隙闭合应力为:

(11)

由于基质材料内部裂隙的纵横比不同, 所以不同裂隙的闭合应力也不同, 因此, 随着σ₃的增大, 基质材料内部裂隙持续发生闭合, 其中未闭合裂隙数N不断减小, 直至所有裂隙完成闭合, 此时N = 0。于是, 由式(8) 可知E = E_s, 即基质材料弹性模量为常数, 不再随σ₃的增大而发生变化。

假设基质材料内部存在N₀条裂隙, 将裂隙纵横比按由小到大的顺序排列即λ₁, λ₂, …, λ_i, …, λ_N0, 并设裂隙纵横比在[0, λ_N0] 且服从均匀分布, 则第i条裂隙纵横比可表示为:

(12)

设第i条裂隙闭合应力为, 则由式(11) 可得:

(13)

而且:

(14)

因此, 当空隙岩石应力σ₃达到 (即) 时, 则有i条裂隙已发生闭合, 而未闭合裂隙数N为:

(15)

利用式(12) 和式(13) 可得:

(16)

令, 将式(16) 代入式(15) 可得:

(17)

于是, 将式(17) 代入式(8) 可得:

(18)

其中:

(19)

将σ₃ = 0代入式(18) 可得出单轴压缩条件下空隙岩石弹性模量为E₀, 即:

(20)

也可将式(20) 改写为:

(21)

将式(21) 代入式(18) 可得出σ₃作用下的E, 即:

(22)

显然, 式(22) 仅适用于σ₃ ≤ A的情况, 如果σ₃ > A, 则基质材料内部所有裂隙已完成闭合, 此时基质材料弹性模量E与岩石骨架弹性模量E_s相同, 不再随σ₃的增大而发生变化, 因此, E随σ₃的变化关系可表示为:

(23)

式中, E_s可视为高围压作用下(基质材料内部裂隙已全部闭合) 空隙岩石的弹性模量; A可利用低围压作用下(基质材料内部裂隙未全部闭合) 空隙岩石弹性模量与围压之间的关系进行拟合分析确定。

3.2 基质材料损伤本构模型

为了建立基质材料变形计算方法, 拟基于岩石统计损伤理论进行研究, 并做如下假设:

(1) 基质材料视为由损伤和未损伤2部分材料并联组成, 并联组合体总应力等于这2部分材料的应力之和, 且损伤材料、未损伤材料与基质材料变形协调。

(2) 轴向承受荷载并发生变形破裂, 近似认为损伤仅发生在轴向, 损伤材料轴向应力等于残余应力; 侧向约束试件变形并提高承载能力, 未损伤材料变形服从广义胡克定律。

基于上述假定, 根据文献[6] 的研究思路, 可得:

(24)

(25)

(26)

(27)

(28)

式中, D为基质材料沿最大主应力方向的损伤变量; 为损伤材料沿最大主应力方向的残余主应力; σ_i和分别为基质材料和未损伤材料沿i方向的主应力; 和分别为未损伤材料和损伤材料沿i方向的主应变。于是, 利用式(24) ~ 式(28) 可得:

(29)

(30)

(31)

根据三轴压缩试验可获得空隙岩石残余抗剪强度指标(残余黏聚力c_r和残余内摩擦角φ_r), 由M-C强度准则可得出残余应力, 即:

(32)

根据韦伯分布统计损伤理论^[6], D可表示为:

(33)

式中, F为基质材料微元体强度; m和F₀为韦伯分布参数。根据M-C强度准则将F表示为:

(34)

其中:

(35)

(36)

式中, c_y和φ_y分别为临界损伤状态下空隙岩石的黏聚力和内摩擦角。将式(2) 和式(29) 联立, D也可表示为:

(37)

再将式(37) 和式(24) ~ 式(26) 代入式(34), F也可表示为:

(38)

于是, 将式(38) 代入式(33), 再将式(33) 代入式(29) ~ 式(31), 可得出基质材料变形计算方法, 即:

(39)

其中:

(40)

(41)

(42)

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

前述已获得裂隙材料和基质材料的应变(和) 分析方法, 于是将式(6) 和式(39) 代入式(2) 可得出空隙岩石损伤本构模型, 即:

(43)

其中:

(44)

(45)

(46)

由该模型可看出, E和F分别由式(23) 和式(38) 确定, 除了包含屈服和残余强度参数之外, 还包含一些非常规参数(如γ, E_p, μ_p, E_s, μ, E₀, A, m和F₀), 须对这些参数的确定方法问题给予解决。鉴于E_s, E₀, μ和A的确定方法已在前述部分介绍, 而γ, E_p和μ_p可采用文献[10] 中的方法确定, 还需给出参数m和F₀的确定方法。

由于三轴压缩试验条件下空隙岩石应力应变曲线具有极值特性, 因此采用文献[10] 中的方法可确定参数m和F₀, 但其中关键在于如何利用式(43) 获得∂σ₁/∂ε₁的解析表达式。如果直接利用式(43) 进行求解, 将使其计算过程及由此获得的参数m和F₀的表达式十分复杂, 不利于应用, 因此需依据空隙岩石非线性变形机理及特征分析进行简化处理。

当空隙岩石进入AB段时, 该应力水平下空隙闭合引起的初始宏观非线性变形已经结束, 式(43) 中表达式γ [1-exp(-σ′/E_p)] 趋于最大值γ。由于三轴压缩试验应力组合方式为σ₁ >σ₂ = σ₃, 则对式(43) 进行化简以计算OA段以后的空隙岩石理论偏应力-应变曲线, 即:

(47)

于是, 利用文献[10] 中的方法可得出参数m和F₀的确定方法, 可分别表示为:

(48)

(49)

其中:

(50)

(51)

由于峰值点处最大主应力σ_sc及其相应的主应变ε_sc需从空隙岩石三轴压缩试验中获得, 而其结果一般为偏应力应变试验曲线, σ_sc和ε_sc不能由此直接获取, 因此需解决其确定方法。设和分别为三轴压缩条件下空隙岩石轴向极限强度及其相应应变, 则:

(52)

(53)

式中, ε_c为仅在σ₃作用下空隙岩石的初始应变。于是, 基于式(43) 中第2式可得:

(54)

值得注意的是, 上述是根据特定围压作用下空隙岩石偏应力应变试验曲线确定的参数m和F₀。为使其能够适应任意围压情况, 需分别建立σ_sc和ε_sc与σ₃之间的关系, 可表示为:

(55)

(56)

式中, c_f和φ_f为峰值应力状态下空隙岩石的黏聚力和内摩擦角; a和b为待定常数, 可利用式(56) 进行线性拟合计算。于是, 将式(56) 代入式(53), 然后将式(53) 和式(55) 代入式(50) 和式(51), 再将式(50) 和式(51) 代入式(48) 和式(49) 可得到适用于任意σ₃下m和F₀的计算表达式。

5 实例分析与验证

前述已建立空隙岩石本构模型并给出了参数的确定方法, 但其合理性和可行性仍需验证, 为此根据红砂岩和煤岩2类岩石的试验资料进行验证。

5.1 红砂岩三轴压缩试验验证

基于林卓英等^[13]对红砂岩进行的围压σ₃分别为0, 15, 40 MPa和60 MPa的三轴压缩试验曲线, 采用本研究方法即可求得相关本构模型参数, 见表 1。然后将其代入式(43) 获得三轴压缩条件下红砂岩偏应力-应变理论曲线, 见图 4, 并将其与红砂岩试验曲线以及文献[6] 中的理论模型曲线进行对比分析, 见图 5。由此可以看出, 本研究模型和方法不仅能反映弹性模量随σ₃增大而变化的特点, 也能反映空隙压密非线性上凹变形特征阶段和残余强度变形特征阶段, 较文献[6] 中的理论模型曲线与试验曲线更加吻合, 表明本研究模型和方法更具有合理性与可行性。

5.2 煤岩三轴压缩试验验证

根据Yumlu等^[18]对煤岩进行的σ₃分别为0, 3, 5 MPa和8 MPa的三轴压缩试验曲线可知, 不同σ₃条件下煤岩弹性模量不发生变化, 可认为经过空隙压密变形阶段后的煤岩不再包含可进一步发生闭合的裂隙, 故取参数A = 0, 并根据煤岩试验曲线利用本研究方法确定其余模型参数, 见表 1, 在此基础上获得了本研究模型理论曲线, 并将其与文献[6] 和文献[10] 中的理论模型曲线以及试验曲线进行比较, 见图 6。由此可以看出, 本研究理论模型曲线较文献[6] 和文献[10] 中的理论模型曲线更加接近于试验曲线; 文献[6] 中的理论模型曲线虽然能反映空隙岩石峰后残余强度变形特征, 但不能反映空隙压密阶段的非线性变形特征; 文献[10] 中的理论模型曲线虽然能反映空隙压密非线性上凹变形特征, 但无法反映岩石峰后残余强度变形特征; 本研究模型兼具文献[6] 和文献[10] 中模型的优越性, 更具有合理性与可行性。

6 结论

本研究考虑三轴压缩条件下岩石空隙压密非线性上凹变形特征、弹性模量随围压变化而变化的特点以及峰值后区变形特征, 在探讨空隙岩石非线性变形机理与特征分析的基础上, 借助细观唯象分析方法和损伤力学理论, 对三轴压缩条件下空隙岩石偏应力-应变过程的模拟方法进行了研究, 得出的结论如下:

(1) 基于空隙岩石临界损伤破坏前后的非线性变形机理, 采用细观唯象分析方法将空隙岩石视为由基质材料和裂隙材料2部分串联组成, 建立了空隙岩石非线性变形细观唯象分析模型。

(2) 基于体积柔度与空隙岩石内部张开裂隙数量密切相关的特点, 引入体积柔度与弹性模量之间以及空隙岩石所受静水压力与裂隙纵横比之间的关系, 建立了基质材料弹性模量分析模型。

(3) 建立了空隙岩石统计损伤本构模型并给出了模型参数的确定方法。通过红砂岩和煤岩的实例分析表明, 本研究模型和方法能够反映出常规三轴压缩条件下空隙岩石空隙压密非线性上凹变形特征、弹性模量随围压变化而变化的特点以及峰值后区变形特征。