0.   引言

在岩石压缩条件下的应变软化阶段, 尤其是在低围压时, 体积应变有相对大的增长的现象已多有报道^[1~5]。这一反常的现象常称之为“体积膨胀”或“剪切扩容”。其原因一方面是由于微裂纹的张开已经超过了闭合^{[1, 6]}; 另一方面是由于滑动块体在凹凸面上的抬升^[6~7]。

在实验室及现场可以探测到压缩条件下岩石体积膨胀产生的多种效应。由于许多宏观现象都与剪切扩容密切相关, 因此, 对剪切扩容的深入理解是极其重要的。Scholz等提出了著名的膨胀扩散假说, 可以用于解释断裂来临之前近场的行为^[8]。体积膨胀被视为一种重要的地震前兆, 前兆明显程度依赖于扩容量的大小^{[2, 9~10]}。与膨胀有关的非弹性变形对于流体传输过程有重要的含义^[4]。在与微裂纹有关的扩容过程中, 受到低围压作用的初始低孔隙度岩石的渗透率可以增加0.3至5个量级^[7]。井眼、隧洞、孔洞、巷道和井筒附近的岩石涌向空腔的大量变形部分是由于剪切膨胀而引起的^[11~12]。

遗憾的是, 在试验中我们不能确保若干岩样仅有不同的剪切扩容特性, 而有相同的物理、力学参数。因此, 剪切扩容的单独影响难于从试验角度进行研究。这一不足可由有限元法或有限差分法的数值试验避免。

考虑非均质性, 一些研究人员已采用数值方法模拟了岩石的破坏过程和宏观力学行为^[13~14]。通常采用威布尔分布函数描述单元强度的随机变化。在数值模拟中, 广泛采用峰后脆-理想塑性的本构关系。但上述研究并未涉及剪切扩容的影响。

拉格朗日元法(FLAC)是一个显式的有限差分程序, 可以有效地模拟地质体材料的屈服极限被达到之后的塑性流动行为^[15]。FLAC的一个突出特点是包含了一个强有力的内嵌的编程语言FISH。FISH允许用户定义新的变量、函数和本构模型。为了适合特殊的需要, 用户开发的函数不仅扩展了FLAC的功能, 而且添加了用户自定义的特色。

最近, 文献[16~17]采用FLAC的内嵌编程语言FISH编制了计算平面应变压缩岩样全部变形特征(轴向应变、侧向应变、体积应变、由侧向应变及轴向应变计算得到的泊松比)的FISH函数。文献[16]研究了扩容角对含单一材料缺陷岩样的破坏过程、模式、前兆及全部变形特征的单独影响。文献[18]编写了于试样内部预置初始缺陷的FISH函数, 研究了扩容角对内部含缺陷试样破坏过程及应力-轴向应变曲线的单独影响。

在文献[18]的基础上, 本文采用FLAC模拟了扩容角对具有初始随机材料缺陷的岩石试样的破坏前兆及变形特征的单独影响。测量了13条剪切带在不同位置的倾角, 并与3个经典理论进行了比较。

1.   本构关系及计算模型

计算模型的几何尺寸、单元划分及边界条件见图 1, 黑色单元代表材料缺陷, 白色单元代表密实的岩石。试样的高度及宽度分别为10cm及5cm。试样两端面是光滑的, 在试样的上端面施加常速度, v₀ =5×10^-10m/时间步, 计算在小变形模式及平面应变状态下进行。将试样划分为若干正方形单元, 单元边长为0.00125m。

图  1  含缺陷的岩样的边界条件

Figure  1.  Model geometry and boundary conditions of a rock specimen with imperfections

下载: 全尺寸图片 幻灯片

在塑性理论中, 由于塑性变形依赖于加载历史, 是不可恢复的, 因而塑性应变及应力之间没有直接的联系^[6]。为此, 引入了两个塑性势函数g_s及g_t用于分别控制剪切流动和拉伸流动。势函数g_s属于非关联流动法则, 依赖于主应力及扩容角Ψ。势函数g_t属于关联流动法则, 依赖于第3主应力。通过一个标量的塑性乘子, 塑性势函数与应变增量及应力建立了联系。

扩容角定义为塑性体积应变增量与塑性剪切应变增量之间的比率, 随着塑性剪切应变的增加, 体积增加时扩容角为正值。扩容角可以通过在剪应力作用下体积膨胀的岩石的剪切试验测量^[6]。对于上地壳的地质体材料, 其塑性行为与应力及剪胀有关, 扩容角的值通常小于内摩擦角^[19~20]。

本文采用了5个计算方案。从方案1至5, 扩容角分别为0°、10°、20°、30°及40°。方案1至5有相同的初始随机材料缺陷分布, 见图 1。初始随机材料缺陷的预置方法见文献[18]。

在弹性阶段, 密实岩石和材料缺陷有相同的本构关系, 弹性模量取为26.6GPa, 泊松比取为0.21。

在FLAC中, 粘结力和内摩擦角与塑性应变之间的关系可以任意指定。这样, 可以方便地模拟岩石在峰后应变软化阶段及残余阶段的力学行为。若粘结力和内摩擦角均被规定为常数, 则可以模拟理想塑性行为。

材料缺陷发生破坏之后经历理想塑性行为。材料缺陷的粘结力和内摩擦角分别为0.1MPa及38°。密实岩石一旦发生破坏先是经历线性应变软化行为(对应应变软化阶段), 然后是理想塑性行为(对应残余阶段)。密实岩石的破坏准则选取莫尔库仑剪破坏与拉破坏复合的模型, 抗拉强度取为0.2MPa。在峰后, 粘结力、内摩擦角与塑性应变的关系见文献[16~18, 21]。

2.   计算结果及分析

2.1   扩容角对破坏单元数目及侧向应变的影响

图 2给出了扩容角对试样最终破坏形态(12000时步)的影响。黑色单元表示这些单元已经发生了剪切或拉伸破坏; 白色单元表示这些单元尚处于弹性状态。图 3给出了扩容角对破坏单元数目-轴向应变曲线的影响。图 3中的黑点代表应变软化行为的开始(即试样上端面压缩应力开始下降的点, 之前是应力峰前阶段, 之后是应力峰后阶段, 应力始降点是由程序计算得到的)。

图  2  扩容角对试样最终破坏形态的影响

Figure  2.  Effects of the dilation angle on the final failure patterns of a specimen

下载: 全尺寸图片 幻灯片

图  3  扩容角对破坏单元数-轴向应变曲线的影响

Figure  3.  Number of failed elements-axial strain curves at different dilation angles

下载: 全尺寸图片 幻灯片

图 3表明, 当轴向应变小于一定值时, 不同扩容角时试样内部的破坏单元数目没有差别, 即:破坏单元数目-轴向应变曲线完全重合。随后, 破坏单元数目-轴向应变曲线发生了分离, 当轴向应变相等时, 高扩容角试样内部的破坏单元数目较多, 即:高扩容角时破坏单元数目-轴向应变曲线比较陡峭。当应力的峰值(对应各个黑点)被达到之后, 试样内部的破坏单元数目继续增加(对应应变软化阶段), 直到试样内部的破坏单元数目不再增加, 达到一个常数, 该常数随着扩容角的增加而增加。破坏单元数目恒定表明试样已处于残余变形阶段。

当扩容角为零时, 破坏单元的数目在应力峰前增加很快, 然而, 当扩容角较高时, 在应力峰前破坏单元数目的增加趋势已经放慢。

图 4给出了扩容角对侧向应变-轴向应变曲线的影响。图 4表明, 随着扩容角的增加, 试样在应力峰后侧向膨胀变得明显。

图  4  扩容角对侧向应变-轴向应变曲线的影响

Figure  4.  Lateral strain-axial strain curves at different dilation angles

下载: 全尺寸图片 幻灯片

2.2   扩容角对计算得到的泊松比的影响

图 5给出了扩容角对计算得到的泊松比-轴向应变曲线的影响。图 5表明, 计算得到的泊松比-轴向应变曲线可以划分为3个阶段。

图  5  扩容角对计算得到的泊松比-轴向应变曲线的影响

Figure  5.  Calculated Poisson' s ratio-axial strain curves at different dilation angles

下载: 全尺寸图片 幻灯片

第1阶段对应初始加载阶段, 计算得到的泊松比随轴向应变的增加而增加, 对此解释见文献[16~18, 21]。在第2阶段, 计算得到的泊松比为常量, 这反映了试样内部的变形为均匀变形。在第3阶段, 计算得到的泊松比快速增加, 这是由于试样发生了显著的侧向膨胀。在第2阶段, 不同扩容角时计算得到的泊松比总是低于0.5 (各向同性材料弹性阶段泊松比的上限)。扩容角为0时, 在第3阶段, 计算得到的泊松比可以超过0.5。然而, 当扩容角较高时(40°), 计算得到的泊松比在应力峰前就可以超过0.5。

在各黑点之右, 计算得到的泊松比可以称之为应力峰后泊松比。对于一些脆性材料, 在后破坏阶段, 通过试验测量得到的泊松比可以大于0.5^[22~26]。

Gold由刚度计算了冰板的泊松比在0.29至0.60之间, 依赖于冰的类型、温度及颗粒的尺寸^[22]。刘新荣等研究了一种典型软岩的变形特征, 他们发现, 泊松比超过了0.5^[23]。徐志伟及殷宗泽给出了粉煤灰在真三轴试验条件下的实验结果, 泊松比可以超过0.5, 在最小主应力方向甚至可达到1^[24]。朱建明等测试了14个岩样, 得出了破坏之后岩石具有摩擦-滑动特征的结论, 并且测量得到的泊松比通常大于0.5^[25]。从这一实验结果可以发现, 泊松比的最大值为1.75。Min及Jing使用离散元法研究了断裂岩体依赖于应力的力学特性及泊松比的边界, 研究发现泊松比高于0.5^[26]。他们建议, 对于断裂的岩石, 仔细重新评估泊松比取值(在0.2至0.3之间)的习惯做法是必要的。

这些反常的现象似乎与泊松比应该小于0.5的常识不一致。实际上, 剪切局部化开始之后, 尤其是在应变软化阶段, 整个试样内部的变形不再保持均匀分布。穿透试样的倾斜的剪切带极大地增加了侧向应变的值。这样, 试样的应力峰后侧向行为将明显有别于弹性阶段。因此, 上述观察到的实验现象是合理的。

当单轴压缩岩样遭受到贯通试样的剪切破坏时, 文献[27]推导了应力峰后泊松比的解析解, 从理论上证实了上述反常的实验现象。对于常规岩样, 应力峰后泊松比的预测值可以达到1.4。

2.3   扩容角对体积应变的影响

图 6给出了扩容角对体积应变-轴向应变曲线的影响。图 6表明, 当扩容角为0时, 体积应变的峰值出现于应力峰值附近。而且, 体积应变总为正。这意味着, 变形后的试样的体积不能大于初始体积。对于扩容角为0时的岩石试样, 应力峰后的体积扩容源于被剪切带所分割的毗邻块体之间的相对滑动。

图  6  扩容角对体积应变-轴向应变曲线的影响

Figure  6.  Volumetric strain-axial strain curves at different dilation angles

下载: 全尺寸图片 幻灯片

单轴压缩岩样遭受到贯通试样的剪切破坏时, 为了准确地估计岩样的应力峰后体积应变, 文献[28]将处于后破坏阶段的试样划分为3部分:两个靠近试样端部的弹性区域和一个中心区域(包括倾斜的剪切带, 剪切带经历塑性变形)。应力峰后总体积应变ε_v可以被划分为弹性和塑性部分:

式中ε_v^e及ε_v^p分别称之为弹性和塑性, 它们的表达式为^[28] :

式中L_e是两个弹性区域的总高度, σ_c是单轴压缩强度, L是试样的高度, E是弹性模量, σ是应力, α是剪切带与应力之间所夹的锐角, w是剪切带宽度, c是描述岩石峰后脆性的软化模量, υ是泊松比, B是试样的宽度。

应当指出, 尽管在上述解析解中扩容角为0, 但是, 在应力峰后的变形阶段岩样的体积增加。文献[28]的不考虑扩容的理论分析支持目前的数值结果(对于扩容角为0时的岩石试样, 试样于应力峰后发生体积扩容)。

当扩容角达到40°时, 在应力峰前剪切扩容发生。最终体积应变变为负值。这意味着, 在压缩条件下, 变形后的试样的体积已大于初始体积。非零扩容角时, 试样应力峰后体积膨胀的原因为:剪切局部化(导致了毗邻块体之间的相对滑动)及剪切扩容(发生于剪切带内部)。

2.4   扩容角对破坏前兆的影响

图 3至图 6中的黑点代表应变软化行为的开始(应力开始下降)。由图 3至图 6可以发现, 在峰前破坏的单元数目-轴向应变曲线、侧向应变-轴向应变曲线、计算得到的泊松比-轴向应变曲线及体积应变-轴向应变曲线明显偏离了线性状态。这些偏离现象可以被视为破坏的一种前兆。偏离越大, 破坏的前兆越明显。因此, 扩容角越高时, 前兆越明显。

当应变软化行为刚发生时, 方案1至5的剪切应变增量等值线图见图 7。可以发现, 当应变软化行为刚发生时, 在试样内部可以观察到多条倾斜的剪切带。这些剪切带可以看作破坏的另一种前兆。另外, 在图 7 (a)中可以看到一个更明显的剪切带网络。在图 7 (c~d)中, 观察不到剪切破坏网络, 而是几条占优势的剪切带。当扩容角为0时, 在峰前, 试样的侧向膨胀相对比较均匀, 然而, 当扩容角达到40°时, 不均匀的侧向膨胀尤其显著。

图  7  方案1至5的剪切应变增量等值线图(应变软化行为的开始时)

Figure  7.  Shear strain contour diagrams for schemes 1 to 5 when strain-softening behavior just occurs

下载: 全尺寸图片 幻灯片

对于方案1至5, 当应变软化行为刚发生时, 预测的最大剪切应变增量分别为5.47×10^-5、8.63×10^-5、2.10×10^-4、2.74×10^-4及3.00×10^-4。因此, 对于高扩容角时的岩样, 破坏的前兆(剪切带)比较明显。这一数值结果与以前针对均质岩样的数值结果一致^[16]。

2.5   扩容角对剪切带倾角的影响

Coulomb理论认为剪切带倾角与剪切应力与法向应力之比达到最大值时的平面的倾角相一致^[29]。这样, Coulomb倾角θ_C可以表示为

式中ϕ是启动内摩擦角。

Roscoe强调应变的重要性, 建议剪切带与零拉伸线的方向一致^[29]。Roscoe倾角θ_R为

式中Ψ是扩容角。

Arthur等通过试验发现, 剪切带倾角等于内摩擦角和扩容角的平均值^[29]。这样, Arthur倾角θ_A可以表示为

这一表达式被Vardoulakis的理论分析^[29]所证实。

当时间步达到12000时, 方案1至5的变形后的剪切应变增量等值线图及体积应变增量等值线图见图 8。为了便于肉眼观察, 在图 8中每个节点的水平及垂直位移均被放大了1.0 ×10³倍。对于方案1及2, 最终于试样内部可以观察到3条充分发展的剪切带, 见图 8 (a~b)。对于方案3至5, 有4或5条剪切带得到了充分的发展, 见图 8 (c~e)。

图  8  当时间步达到12000时方案1至5的剪切应变等值线图

Figure  8.  Shear strain contour diagrams for schemes 1 to 5 at 12000 timesteps

下载: 全尺寸图片 幻灯片

利用未变形的剪切应变增量等值线图(限于篇幅未给出), 测量了13条得到充分发展的剪切带倾角, 并与Coulomb理论、Roscoe理论及Arthur理论进行了比较, 见图 9。在此, ϕ被假定等于密实岩石的内摩擦角, 即:ϕ=44°。图 9 (a)给出了方案1至5靠近试样左下角的剪切带倾角的数值解。图 9 (b)给出了方案1至5靠近试样右下角的剪切带倾角的数值解。图 9 (c)给出了方案3至5与试样左下角剪切带平行的剪切带倾角的数值解。

图  9  剪切带倾角的计算结果与Coulomb理论、Roscoe理论及Arthur理论的比较

Figure  9.  Comparison of the numerical shear band inclination and the predicted results using classical theories

下载: 全尺寸图片 幻灯片

由于一条剪切带的两个边界可能不平直, 剪切带倾角可能发生改变。因此, 一条剪切带在不同位置有不同的倾角。目前关于曲线剪切带的数值模拟结果类似于采用X射线及光学测量方法观测到的实验结果^[30]。他们通过对试样逐一切片发现, 剪切带的边界不是直的, 而是稍微弯向最大主应力方向。

当充分发展的剪切带贯穿试样的侧面时, 剪切带的倾角在45°左右, 低于内部剪切带的倾角。图 9不包括这些(低倾角)剪切带的倾角的数值结果。Needleman及Ortiz认为驻波微分方程的通解(以瑞利表面波为例)的出现与存在自由边界时应变局部化的开始相对应^[31]。针对压力敏感摩擦材料的分析揭示, 固体内部及边界上的剪切带方向出现不匹配。他们发现, 对于塑性不可压缩固体, 在边界上对称的剪切带与加载轴成45°。然而, 内部的剪切带可以十分陡峭。这一现象已在其它的使用有限元法^[32]及FLAC的数值试验^[33]中观测到。

充分发展的剪切带的倾角的数值解与Arthur理论的预测结果更接近, 见图 9, 并且低于Coulomb理论的预测结果。通常, Roscoe理论低估了剪切带倾角, 除非扩容角较高时。

在低扩容角时, 靠近试样右下角和左下角的剪切带的倾角的数值解比较分散, 这是由于剪切带的两边界不是直的, 见图 9 (a~b)。弯曲的剪切带出现的原因或许是由于试样内部初始材料缺陷的随机分布。

3.   结论

1) 随着轴向应变的增加, 试样内部破坏的单元数目增加, 直到达到一个常数, 该常数随着扩容角的增加而增加。当扩容角较高时, 计算得到的泊松比在峰前就可以超过0.5, 剪切扩容于峰前发生; 最终, 试样的体积应变变为负值(体积大大膨胀)。剪切局部化(导致了毗邻块体之间的相对滑动)及剪切扩容(发生于剪切带内部)是非零扩容角试样峰后体积膨胀的原因。

2) 在峰前, 破坏的单元数目-轴向应变曲线、计算得到的泊松比-轴向应变曲线及体积应变-轴向应变曲线均发生了偏离线性状态的现象。通过观察上述曲线偏离线性状态的量及剪切应变增量可以发现, 扩容角越高, 试样破坏的前兆越明显。

3) 在低扩容角时, 对试样内部的充分发展的剪切带的倾角的测量表明, 剪切带倾角的数值解比较分散, 这是由于弯曲的剪切带边界。充分发展的剪切带的倾角更接近Arthur理论的预测结果, 低于Coulomb理论的预测结果。通常, Roscoe理论低估了剪切带倾角, 除非扩容角较高。