由于受成岩过程中不同环境与作用力的影响, 岩石成为一种内部包含微裂隙、空隙等微缺陷的多晶材料^[1-2].从细观尺度观察到岩石主要由矿物晶粒和微结构缺陷组成, 晶粒的力学特性和晶界结构共同决定了岩石的宏观力学性质.当有外力作用时, 首先在岩石局部晶界出现微裂纹成核、萌发、扩展, 然后岩石局部发生破裂并出现应力集中, 使已有裂纹进一步扩展和新裂纹生成, 最终大量微裂纹联合形成宏观裂纹导致岩石破坏^[3].

近年来, 国内外研究人员对应力作用下岩石裂隙发育机制做了大量的研究工作.例如:Ge等^[4]进行了三轴和单轴荷载作用下煤岩破坏全过程的细观损伤演化规律的实时动态CT观察；杨永明等^[3]开展致密砂岩三轴压缩试验及CT扫描试验, 分析了不同围压对岩石试样内部破坏裂纹分布的影响规律；苏承东等^[5]对煤样进行常规单轴、三轴和三轴卸围压声发射(AE)试验, 研究了煤样在不同应力路径下变形破坏过程中声发射特征.由于室内试验很难从细观尺度获得试样内裂纹萌生—扩展—贯通的演化过程, 许多学者借助数值软件模拟研究了岩石宏-细观力学响应.例如:Wang等^[6]采用近场动力学对单轴压缩条件下岩石破裂过程进行了三维数值模拟；穆康等^[7]基于颗粒流程序PFC从裂纹类型的数目演化和角度分布等细观层面分析了不同围压作用下岩石的宏观力学响应；Chen等^[8]和Gao等^[9]利用二维块体离散元(UDEC)建立颗粒离散元模型(UDEC-GBM)，分析了岩石细观损伤演化过程.

基于上述研究现状, 本文首先利用开源软件Neper生成Voronoi晶粒镶嵌体, 将其导入3DEC建立考虑矿物颗粒成分的三维多晶离散元模型, 然后基于砂岩室内三轴压缩试验结果对模型的细观参数进行标定, 最后运用三维多晶离散元模型从细观尺度对非均质砂岩在不同围压下的损伤变形演化进行数值分析, 从试样内部裂纹数目演化等方面探讨砂岩宏观力学响应和破坏的细观机理.

天然砂岩是不同形状不同种类的矿物晶粒黏结形成的沉积岩, 如图 1a所示.为了真实体现砂岩的这种细观结构, 本文利用不同形状和大小Voronoi晶粒镶嵌建立模型，如图 1b所示.Voronoi晶粒由开源软件Neper(一款在Linux系统中运行, 用于多晶体生成和网格划分的软件)生成^[10], 导入3DEC建立三维多晶离散元模型, 模型中每一个颗粒都是凸多面体, 2个相邻颗粒之间共用一个平面, 3个相邻颗粒之间共用一条边, 4个或更多颗粒之间共用一个顶点.该模型的主要特点为:①更真实反映砂岩细观物理结构；②由于晶粒几何形状的不规则性, 外力作用下形成的微裂纹方向随机, 最终形成的裂纹曲折不规则, 与实际裂纹形态更接近.

图 1(Fig. 1)

图 1 砂岩扫描电镜图片与数值模型 Fig.1 SEM photo and numerical model of sandstone (a)—砂岩扫描电镜图；(b)—三维Voronoi镶嵌模型.

本文以砂岩为对象, 通过实验室试验和数值模拟研究细观尺度下非均质岩石的三轴压缩损伤破裂特征.该砂岩产自云南昆明, 主要矿物成分(质量分数)为:石英60 %, 长石19 %, 方解石8 %及黏结物13 %.按砂岩真实的矿物成分含量比例, 利用3DEC内嵌的FISH语言编写程序, 建立如图 2所示直径50 mm、高度100 mm的非均质多晶体离散元模型, 不同颜色代表不同种类的矿物晶粒.由于计算能力的限制和生成太小晶粒导致不良几何关系的原因, 建立与岩石真实矿物晶粒尺寸相同的数值模型是不现实的, 为平衡计算精度和计算成本, 本文所建立的模型包含5 300个颗粒, 忽略孔隙度对模型力学性质的影响^[11].本文通过监测模型上端面8个节点z轴方向上位移和计算其平均值得到轴向应变；通过监测模型中间平面上6个节点向平面圆心的位移和计算其平均值得到侧向应变；通过监测模型中心20 mm×20 mm×40 mm范围内的所有单元z轴向应力, 计算其平均值得到轴向应力.

图 2(Fig. 2)

图 2 非均质砂岩几何模型 Fig.2 Geometric model of heterogeneous sandstone

三维多晶离散元模型的宏观力学行为主要由颗粒和颗粒之间接触面的力学性质共同决定.本文中颗粒设置为可变形体, 颗粒本构模型采用线弹性模型, 接触面本构模型采用库伦滑移模型, 通过最大拉伸准则和摩尔-库伦滑移准则判断接触面的拉伸和剪切破坏, 如图 3所示.

图 3(Fig. 3)

图 3 颗粒接触面本构模型 Fig.3 Constitutive model of grain contacts

文献[12]在花岗岩压缩试验和现场圆形隧道开挖监测过程中观察到岩石脆性破裂中存在拉伸裂纹形成先于剪切破裂现象, 涉及到岩石内聚力弱化、摩擦强度增强, 且两者不同时进行的过程, 这种现象称为“内聚力弱化-摩擦强化”(cohesion weakening and frictional strengthening，CWFS).本文的接触面本构模型在库伦滑移本构模型基础上也采用CWFS概念, 设置残余内聚力和峰值摩擦角为0.

外力作用下岩石局部产生微破裂的同时以弹性波的形式释放一部分应变能的现象称为声发射, 通过监测声发射信号可以掌握随时间变化的岩石破裂空间结构和损伤劣化程度.唐春安等^[13]将外力作用下岩石材料体积单元中产生的微破裂(声发射事件)比率定义为损伤量, 并证明了定义的合理性.鉴于此, 三维多晶离散元模型中考虑岩石的细观破裂状态, 将损伤量D定义为在外力作用下模型内部的晶粒之间子接触面(模型中颗粒划分单元时每个接触面会自动划分为多个子接触面)破裂数量的比率, 用于定量描述岩石宏观性能的劣化过程.

(1)

式中:m为模型运行的时间步；n_i为模型运行的第i步破裂的子接触面数量；M为模型内不再生成新裂纹时运行的总时间步.

对于颗粒或块体离散元模型(PFC, 3DEC等), 岩石材料的宏观力学行为主要由细观尺度上颗粒和颗粒之间接触面的力学参数控制, 而这些力学参数无法直接从实验室获得, 根据岩石三轴压缩等试验中获得的宏观物理量(峰值强度、弹性模量等), 通过“反复试验法”^[9]可得到合适的细观力学参数, 参数标定过程如图 4所示.标定结果得到各种矿物晶粒的细观力学参数如表 1和表 2所示.不同矿物晶粒黏结形成岩石, 如果直接将各种矿物的细观力学参数赋值给模型, 模拟得到的砂岩试样强度远远大于真实岩石的强度值, 说明真实岩石内不同种类矿物胶结时晶粒之间接触面强度在一定程度上被弱化, 低于理想条件下纯矿物体内晶粒之间接触面强度, 目前还没有合适的方法估计两种相邻矿物晶粒之间接触面的强度值, 为了得到合理的晶粒接触面的强度值, 通过一系列不同强度的细观标定, 本文最终将不同种类矿物晶粒之间接触面强度值折减为0.5倍.

图 4(Fig. 4)

图 4 模型细观参数标定流程 Fig.4 Calibration process of model microscopic parameters

表 1(Table 1)

表 1 矿物颗粒细观参数 Table 1 Microscopic parameters of mineral grains 矿物 泊松比 石英 37.0 44.0 0.19 长石 37.5 15.0 0.27 方解石 73.3 32.0 0.28 黏结物 9.9 4.4 0.30

表 1 矿物颗粒细观参数 Table 1 Microscopic parameters of mineral grains

表 2(Table 2)

表 2 颗粒接触面细观参数 Table 2 Microscopic parameters of grain contacts 接触面 Kn/Ks φr/(°) 石英/石英 8e14 2 48.7 198 6 长石/长石 3e14 2 40.0 160 7 方解石/方解石 4e14 2 25.0 100 6 黏结物/黏结物 4e13 2 11.2 44.2 6 注：Kn为法向刚度；Ks为剪切刚度；JT为拉伸强度；JC为内聚力；φr为残余内摩擦角.

表 2 颗粒接触面细观参数 Table 2 Microscopic parameters of grain contacts

由于数值软件3DEC中没有施加法向力的边界条件命令, 导致模拟过程中无法对圆柱模型施加围压^[14].为此, 本文通过FISH语言编程完成围压的施加, 施加方法主要包括以下步骤:1)搜索模型侧面边界所有节点并提取每个节点坐标值；2)设置每个侧边界节点上径向合力值为预设围压值.将合力值分解为X, Y轴向力, 并通过边界命令施加于侧边界；3)施加过程中, 每个侧边界节点上施加的X, Y轴向力随时间步按比例增大, 直至达到预定压力值, 以保证整个模拟过程中模型处于准静态平衡状态；4)提取每个侧边界节点上X, Y轴向力并求矢量合力.所有侧边界节点上的合力相加求和, 并除以侧边界面积得到围压值, 实时监测围压值大小.根据监测得到的围压值实时微调整侧边界节点上施加的X, Y轴向力, 以保证围压值恒定.

对于三维多晶离散元模型, 外力可以直接施加于侧边界，控制模型的围压值, 但模拟过程中发现模型局部容易发生破坏, 无法得到有效模拟结果.出现这种问题主要是由于直接施加于侧边界节点的力是非均匀的(图 5a), 导致局部出现应力集中, 某些颗粒急速变形破坏, 导致模型整体失效.为此, 本文考虑侧边界增加厚度为0.001 m的围压板, 围压板直接与模型侧边界接触, 如图 5b所示, 40个板块组成围压板, 板块之间相互接触, 但无黏结.在围压板外侧施加外力保证了模型受到均匀围压, 确保模拟有效进行.通过单轴压缩模拟试验证明了围压板的存在不影响模拟过程中模型的物理性质.另外, 本文的位移加载速率设置为0.05 m/s.

图 5(Fig. 5)

图 5 模型围压加载方式 Fig.5 Loading method of confining pressure in model (a)—围压直接加载示意图；(b)—围压板结构示意图.

采用TAW-2000 kN微机控制电液伺服岩石三轴试验系统开展了常规三轴压缩试验, 其围压为0, 5, 10, 20, 30和40 MPa.每个相应的围压值及单轴压缩试验各进行3次, 以减少岩样的离散性对试验结果的影响.试验前, 先将试件用聚烯烃热缩管包裹, 然后试件与应变测量仪一起放入三轴腔, 每个横向和纵向变形的测量值均为4个测点的平均值.轴向采用位移静态加载方式, 加载速率为0.002 mm/s, 侧向围压采用电液伺服控制, 加载速率为0.1 MPa/s.不同围压下应力-应变曲线如图 6所示.从图 6可知, 当施加围压从5 MPa增加到40 MPa时, 砂岩试件的峰值强度从99.1 MPa增加到212 MPa, 峰值应变从0.99 %增加到1.81 %, 随着围压增大，试样压密阶段减弱.

图 6(Fig. 6)

图 6 不同围压下砂岩试件三轴应力-应变曲线 Fig.6 Triaxial stress-strain curves of sandstones under different confining pressures

开展与实验室常规三轴压缩试验相对应的三维多晶离散元模型三轴压缩数值模拟试验, 模拟结果如图 7~图 10所示.自然状态下的砂岩试件内部存在大量不同方向的微裂纹和微孔洞, 在室内试验加载初期试样的应力-应变曲线出现下凹(即压密阶段)(图 6).因为建立的模型没有考虑岩石存在的初始缺陷条件, 所以数值模拟结果中不存在初始压密阶段.对于岩石而言, 考虑到试验中通常以弹性阶段斜率计算求得其弹性模量, 因此采用向左平移试验曲线的方式消去初始裂隙压密阶段^[15-16], 从图 7可以观察到单轴压缩条件下的模拟结果能够较好地拟合室内试验结果.图中σ_c和σ_cd分别为峰值应力和屈服应力.

图 7(Fig. 7)

图 7 单轴应力-应变曲线 Fig.7 Uniaxial stress-strain curves

图 8(Fig. 8)

图 8 砂岩试件的莫尔包络线 Fig.8 Mohr' s envelope of sandstone specimen

图 9(Fig. 9)

图 9 不同围压下应力-应变曲线及微裂纹数量-应变曲线 Fig.9 Stress-strain curves and microcrack numbers-strain curves under different confining pressures (a)—围压0；(b)—围压5 MPa；(c)—围压10 MPa；(d)—围压20 MPa；(e)—围压30 MPa；(f)—围压40 MPa.

图 10(Fig. 10)

图 10 不同围压下试件损伤演化特性 Fig.10 Damage evolution of specimens under different confining pressures (a)—损伤-应变曲线；(b)—损伤演化示意图.

根据5组不同围压下砂岩三轴压缩数值模拟结果, 可以拟合得到其Mohr-Coulomb强度包络线, 如图 8所示, 得到砂岩试样的内聚力为20 MPa, 内摩擦角为32.55°, 这与室内试验结果较为接近.另外, 本数值模型在不同围压下的抗压强度与室内试验结果拟合较好(见表 3), 说明模型能够模拟砂岩的基本力学性质.

表 3(Table 3)

表 3 单轴和三轴压缩试验与模拟结果 Table 3 Uniaxial and triaxial compression testing and simulation results 围压/MPa 抗压强度/MPa 误差/% 试验值 模拟值 0 52.5 52.3 0.2 5 98.3 91.5 6.9 10 115.4 113.9 1.3 20 160.2 149.1 6.9 30 182.1 182.8 0.4 40 212.0 217.7 2.7

表 3 单轴和三轴压缩试验与模拟结果 Table 3 Uniaxial and triaxial compression testing and simulation results

由图 9可见, 砂岩试件三轴压缩应力-应变曲线中峰值应力之前试样轴向和侧向应变都随围压增大而增加, 且随围压增大试样塑性变形也不断增大, 屈服阶段逐渐明显.试样由脆性逐渐转化为延性, 峰值后应力跌落幅度变缓, 残余强度增大.

试样宏观上的变形破裂与其内部微裂纹的扩展演化密切相关, 在弹性阶段仅有少量微裂纹, 但一旦达到屈服应力点, 微裂纹开始持续快速生成, 试样发生不可逆转的塑性变形.不同围压下试样内总微裂纹数量演化规律类似, 但生成拉裂纹和剪裂纹的演化规律不同.当围压小于10 MPa时, 微裂纹以拉裂纹为主, 且弹性阶段几乎不生成剪裂纹；当围压增加到20 MPa时, 弹性阶段生成剪裂纹但依然以拉裂纹为主, 最终拉裂纹主导了试样的破坏；当围压增加到40 MPa时, 试样中首先生成剪裂纹, 出现了剪裂纹数量大于拉裂纹的情况.说明随着围压增大, 试样内剪裂纹数量所占百分比增大, 低围压时试件以张拉破坏为主, 随着围压的增大, 剪切滑移破坏区域逐渐增大, 张拉破坏区域逐渐减少, 即高围压时试件以剪切滑移破坏为主.这从某种程度上反映了围压对岩石内部孔裂隙的拉伸破裂发展具有一定抑制作用.这种裂纹演化规律的变化也反映了围压效应下砂岩试样脆性转延性的特性, 即以拉微裂纹为主的细观演化体现出砂岩脆性破坏特征, 而剪微裂纹大量发育体现出试样脆性转延性特征.另外, 随着围压增大, 砂岩试件扩容点应力占强度百分比增大, 例如围压从5 MPa增加到30 MPa, 相应的扩容点应力占强度百分比从70 %增加到82 %, 这说明围压越大扩容现象出现越滞后, 围压对砂岩试样内部微裂纹的发展具有一定抑制作用.图中ε_d代表体积扩容点.

由图 10可以看出, 不同围压条件下砂岩试样的损伤破坏演化规律是一致的, 基本可以划分为4个阶段:初始损伤阶段(Ⅰ), 有很少量的微裂纹生成, 损伤量几乎为零, 砂岩岩样处于弹性阶段前期；损伤发展阶段(Ⅱ), 砂岩试样处于弹性阶段后期, 其内部一些局部应力集中位置陆续发生破裂生成微裂纹, 微裂纹萌生速率加快, 损伤演化加快, 但总体并不剧烈；损伤加速阶段(Ⅲ), 几乎出现在屈服应力之后, 砂岩试样内部微裂纹生成持续在一个较高数值, 微裂纹不断生成、汇合和贯通, 试样最终形成宏观破裂面；损伤回落阶段(Ⅳ), 试样内微裂纹生成减慢, 损伤量趋于平缓, 试样基本处于破坏后阶段.砂岩试件的整个损伤演化过程很好地反映了砂岩内部微裂纹萌生、扩展、贯通直至试样破坏的渐进演化过程.这与文献[7]的研究结果是一致的.

利用开源软件Neper生成Voronoi晶粒镶嵌模型导入3DEC，建立更能真实反映岩石材料结构(例如:矿粒胶结形成的砂岩)的三维多晶离散元模型, 考虑砂岩矿物成分, 模型中引入强度非均质性.考虑围压效应得到砂岩压缩变形过程中的细观演化规律.模拟结果表明:低围压时, 模型内生成拉裂纹, 且拉裂纹的数量远远多于剪裂纹, 说明拉裂纹主导岩石的破裂, 岩石偏向于发生脆性破坏；但随着围压增大, 剪裂纹数量增多, 且试样所受围压在大于20 MPa时, 其内部生成的剪裂纹数量大于拉裂纹, 岩石呈现出延性变形, 扩容现象也随围压增大而更明显(相应的扩容点应力占强度百分比增大), 最终发生剪切破坏.砂岩的压缩变形过程是一个渐进损伤演化过程, 基本可以分为4个阶段:初始损伤阶段、损伤发展阶段、损伤加速阶段和损伤回落阶段.