基于三维裂隙网络的岩体剪切特性尺寸效应分析

宋盛渊,黄 迪,隋佳轩,陶 勇,马牧野,李豪杰

(1.吉林大学 建设工程学院,长春 130000; 2.吉林省地质环境监测总站,长春 130000)

裂隙作为一种重要的不连续面广泛存在于岩体中,对岩体的结构及力学性质起重要作用,影响着岩体的强度和破坏模式。不同尺寸岩体试样包含的裂隙存在差异性,导致裂隙岩体的几何和力学参数存在明显的尺寸效应,当岩体参数达到基本稳定时的试样尺寸即为表征单元体(REV)[1]。表征单元体(REV)是能够反映材料统计平均性质的最小体积,从岩石材料的角度讲,它是表征岩体尺寸效应的参量[2-3]。在工程实践中,研究大尺度岩体所需要的力学参数来源主要是通过实验室测定小尺度岩体参数推算的,但是由于岩石材料的不均匀性和裂隙发育的随机性,这种方法的可靠性难以明确[4-5]。岩体尺寸效应的研究为获取大尺度岩体力学参数提供了一个新方法,它将岩块的实验数据与岩体的力学性质关联起来,在评估岩体力学性质方面具有重要意义[6],是获取岩体的等效参数的重要手段。

首先,当达到几何REV时,裂隙岩体在结构特性上就具有了统计学意义,结构稳定后,才能逐步达到能够用于工程计算和模拟的力学REV。众多学者用不同的几何参数探讨裂隙岩体的尺寸效应,可分为裂隙网络的几何特性和完整岩体的几何特性。前者包括裂隙的密度[7]、张开度[8]以及三维裂隙连通率[9-10]等;后者包括能够体现岩石完整性的块体百分比[11],表征岩体质量好坏和连续程度的岩石质量指标(RQD)[12],以及岩体结构发育程度的地质强度指标(GSI)[13]等。基于此,文献[14]进一步对块体百分比加以修正,提出块体体积综合百分比用以研究裂隙岩体的REV尺寸。以上研究均从几何角度对裂隙岩体的尺寸效应进行研究。

但在岩体工程计算中,掌握岩体的力学性质是合理计算的关键,离开材料力学性质谈REV没有任何意义[15]。因此,确定岩体的力学REV也至关重要。众多学者通过试验法、解析法和数值模拟法研究裂隙岩体的力学参数尺寸效应。文献[16]对现场岩体进行承压板试验和承载力试验进行变形模量的计算,求其REV尺寸,但受到试验设备及成本的限制,所测得的岩石面积分别仅有2 000 cm2和500 cm2。文献[17]通过3D打印技术制作裂隙网络类岩体试件,对120 cm以内的岩体进行单轴压缩试验研究尺寸效应。可见,试验法受到室内和现场条件的影响较大,在大尺度节理岩体的尺寸效应研究上受到限制。解析法将岩体的力学参数与岩石块体间的力学参数通过公式联系起来,提供了一种计算REV尺寸的途径。文献[18]依据中心极限定理推导了REV的数学表达式,其中涉及到了众多参数,例如,应力场、弹性形变势能、岩块及结构面的柔度矩阵、结构面的厚度等。解析法的公式多设计到的参数众多,且获取困难,这在一定程度上制约了其发展。随着计算机的发展,多种数值模拟技术被应用到裂隙岩体的尺寸效应研究中。文献[19-20]通过有限元方法对不同裂隙岩体进行REV特征分析,将岩体视为岩石和裂隙的二元结构。相对于上述有限元方法,离散元更加适用于非连续介质模型的模拟,也被应用到裂隙岩体的研究中。文献[21]采用Goodman节理单元模拟随机节理,通过数值模拟计算变形参数,最终将REV所对应的力学参数应用到的边坡的稳定性计算中。文献[22]运用颗粒流软件开展岩体单轴压缩试验,探究了直线型和分形节理模型峰值应力的尺寸效应,结果表明,二者在尺寸效应规律上较为一致。以上研究多从变形和抗压强度的角度对裂隙岩体的尺寸效应进行诸多研究。值得注意的是,剪切破坏也是边坡失稳的重要因素,当下滑力大于抗滑力时,就会产生滑移面导致边坡失稳,而滑动面的抗滑力主要由岩体的抗剪强度来提供。除此之外,在岩体工程连续介质模型的计算中抗剪强度指标也是必不可少的参数,尤其是在运用强度折减法时,主要是对内聚力内摩擦角进行折减,因此,研究岩体的剪切特性尺寸至关重要。然而目前多数抗剪强度尺寸效应的研究多集中在结构面上,关于岩体的相关研究较少。文献[23]通过室内和现场直剪试验对柱状节理玄武岩以抗剪强度参数为指标进行研究。由于受到试验条件的限制,仅对小尺度岩块进行了研究,研究的岩块最大尺寸为100 cm。数值模拟为大尺度岩体研究提供了途径。

鉴于此,本文基于三维裂隙网络通过颗粒流PFC3D软件进行大尺度裂隙岩体直剪试验数值模拟,从剪切特性的角度定量研究裂隙岩体的尺寸效应,探讨抗剪强度和相关参数的尺寸效应并确定REV尺寸,拟合力学参数和试样尺寸的函数关系。

1 坝肩岩体三维裂隙网络模拟

现场结构面多以二维迹线的形式出露,baecher圆盘模型将三维空间内的裂隙视为薄圆盘,裂隙的特征参数可由圆盘模型简便且准确地表征,圆盘直径、方向分别表征裂隙的直径和产状。三维裂隙网络模拟即是在地质结构面现场采样的基础上,以概率理论为依据,利用迹长校正、直径估算等方法将二维迹线扩展为三维空间中的薄圆盘[24]。三维裂隙网络模拟基于蒙特卡洛原理,认为裂隙的特征参数在空间内服从某种特定的分布,将现场采样得到的随机性裂隙进行分布概率拟合,蒙特卡洛抽取随机数以逼近这种分布函数,即可得到与现场结构面分布一致的裂隙网络。

1.1 研究区概况

怒江松塔水电站(图1(a))位于青藏高原东南缘,区内地势总体北高南低,呈阶梯式下降,且具由西北向南东掀斜特征。坝址河段属中山-高山地貌单元,河谷主要呈较对称的“V”形,宽度为80～100 m,两岸岸坡坡度为30°～60°(图1(b))。坝址区地质构造复杂,断裂较发育,为“入”字型构造体系,节理裂隙发育,可观察到岩脉侵入,如图1(c)和1(d)所示。坝址区岩性主要为燕山晚期黑云二长花岗岩,粗粒变晶结构,整体块状构造。岩体受构造影响,轻微变质,岩石、矿物普遍碎裂,局部糜棱岩化,矿物定向分布形成片麻状和片理状构造等现象。大型水电站对坝肩岩体的承载能力及抗滑能力要求较高,其力学参数特别是抗剪能力的研究对工程的安全性至关重要,因此本文以西藏怒江松塔水电站右岸坝肩PDC3平硐内的裂隙岩体为例,开展三维裂隙网络模拟和直剪试验数值模拟以确定其剪切特性REV大小。现场结构面多以二维迹线的形式出露,采用矩形窗口法在右岸坝肩的平硐PDC3人工测量裂隙,共采集到裂隙247条。以二维迹线图的方式描述裂隙,如图2所示。

图1 研究区概况

图2 裂隙二维迹线图

1.2 优势分组及产状概率分布

采用文献[25]提出的QPSO-FCM全局优化算法对结构面产状进行优势分组,如表1和图3所示。划分优势组数为3组,3组裂隙数分别为31、37、179,第1组裂隙陡倾且与坡面大角度相交,第2、3组裂隙分别陡倾坡内和缓倾坡外。平均产状分别为180.3°∠89.8°、278.2°∠85°、118.1°∠32.2°。Fisher分布是表征裂隙产状最常用的概率分布。其概率密度函数为

(1)

表1 三维裂隙网络模拟参数

图3 裂隙优势分组结果

式中:θ为裂隙倾角,φ为裂隙倾向,k为Fisher常数。

Fisher常数k用于描述裂隙产状的离散程度,k值越大,裂隙产状越集中,反之,裂隙产状越离散。常数k计算公式为

(2)

1.3 直径分布拟合及计算

由于采样窗口的限制,很难直接获取裂隙真实迹长,即裂隙迹长存在取样偏差,因此首先要对现场取样迹长进行校正[26]。迹长校正公式[27]为

(3)

式中:w为采样窗口宽度,h为采样窗口高度,R0为两端相交型裂隙的占比,R2为两端可见型裂隙的占比,θ0为裂隙真倾角。

卡方检验用于对校正后的样本数据分布和标准分布之间的差异性进行判断,卡方值越小,表示越服从某种分布。对常见的迹长分布类型进行卡方检验,组1和组2的迹长服从正态分布,组3服从对数正态分布。基于校正后得到的迹长值,采用文献[28]提出的方法进行直径大小的估算。计算结果见表1。

1.4 密度计算

采用张量法计算裂隙在三维空间内的体积密度P32,计算公式[27]为

(4)

1.5 三维裂隙网络生成及验证

利用PFC3D软件中DFN模块编制fish语言,基于蒙特卡洛原理将表2中确定的裂隙参数组合,生成尺寸40 m×40 m×40 m的三维裂隙网络,包含34 250条裂隙,3组裂隙数分别为2 527、2 179以及29 544(图4)。

表2 数据对比

图4 三维裂隙网络

蒙特卡洛模拟按照裂隙参数随机地生成三维裂隙网络,因此模型具有一定随机性和误差,准确度还需要进一步检验。将三维裂隙网络中裂隙信息导出,进行直径分布拟合计算以及产状、密度的计算,与现场裂隙数据对比,以此来检验模型的准确度,见表2。

2 等效岩体模型

2.1 等效岩体技术

颗粒流中以平行黏结和光滑节理两大模型为基础构建等效岩体。颗粒被假设为刚性体,在颗粒间赋予平行黏结模型,在外力作用下相互运动,传递力和力矩,当颗粒间的接触被破坏时,单个颗粒本身不会发生破坏,而是在两个颗粒之间的平行黏结断裂处产生微裂纹,微裂纹积累导致破裂面产生,引起模型整体强度降低,这符合岩石材料的力学特性和破坏特征,可有效表征完整岩石的力学行为[29];光滑节理模型可反映结构面摩擦性质,摩擦行为通过将光滑节理模型分配给两侧颗粒之间所有的接触实现,允许颗粒穿过和滑动,更好地反映试样的脆性破坏特征[30-31]。等效岩体技术即将裂隙网络嵌入到完整岩石中,在完整岩石被裂隙交切部位的平行黏结模型替换为光滑节理模型,其原理见图5。

图5 等效岩体原理

基于上述原理,将平行黏结和光滑节理模型中所涉及的细观参数与实验室所得的完整岩石及结构面的宏观力学参数对标,保证模型能够表征真实岩体。从裂隙网络中每隔1 m切割出一个正方体裂隙网络(图6),同时,处于边界上的裂隙也切割处理,根据等效岩体技术,嵌入至相应尺寸的颗粒体模型中。

图6 裂隙网络二维切割示意图

2.2 细观参数标定

细观力学参数表征模型的力学性质,当试样宏观力学参数与室内试验结果一致时,便可将该组参数用于实际计算模型根据试验结果统计,坝址区花岗岩单轴抗压强度平均值89 MPa,弹性模量30 GPa。结构面内摩擦系数0.70,内聚力0.10 MPa。

2.2.1 平行黏结模型参数标定

文献[32-33]对颗粒尺寸和试样模型的大小关系作出了研究,分别认为模型高度与颗粒尺寸的比例大于100、125、200时对宏观力学参数和破坏模式影响较小,但二者只考虑了模型高径比为2∶1的情况。针对高径比1∶1时,考虑到颗粒数量对模拟速度的影响,将试样高度与颗粒尺寸的比例设置为62.5,把颗粒数量维持在一定范围内,保证模拟的效果和计算速率。由于不同尺寸模型中颗粒尺寸不一致,因此需要每一个模型都进行细观参数标定,采用试错法进行多次单轴压缩数值模拟试验不断地调整细观参数直到所得到的宏观力学参数与室内试验结果匹配。单轴压缩数值模拟试验如图7所示。参数标定结果见表3,标定的宏观参数与实验室力学试验结果误差均在可控范围内,可以用于表征完整花岗岩的力学性质。

表3 平行黏结模型细观参数

图7 单轴压缩数值模拟试验

2.2.2 光滑节理模型参数标定

构建尺寸为0.1 m×0.1 m×0.1 m且包含贯通裂隙面的立方体试样,裂隙设置在试样的中间。剪切盒由8面“墙体”组成,其中,左右两侧均有两面“墙体”分别组成上下剪切盒,为了方便剪切试验的模拟,上下剪切盒之间设定微小间距。在进行直剪试验时,移动上剪切盒,同时对顶墙施加速度使法向压应力分别维持在1、2、3、4 MPa直至试样被剪坏,根据峰值抗剪强度画出抗剪强度包络线。通过不断地调整细观参数直至与室内结构面直剪试验结果匹配,求出内摩擦系数为0.693,内聚力为0.1 MPa,符合实验室试验结果,其中,法向刚度和切向刚度均为8×109MPa,摩擦系数为0.8。另外,对于有胶结填充的黏结型裂隙,内聚力由摩擦系数、法向黏结强度和切向黏结强度共同提供,现场裂隙多为非黏结裂隙,其内聚力仅有摩擦系数提供,法向及切向黏结强度设为0。抗剪强度包络线如图8所示。

图8 裂隙抗剪强度包络线

3 REV尺寸确定

通过PFC3D对20个不同尺寸的等效岩体依次进行2、3、4 MPa正应力下的直剪试验,做出其抗剪强度包络线,求出内摩擦系数和内聚力。根据不同正应力下的抗剪强度,以及内摩擦系数和内聚力判断REV尺寸。

3.1 抗剪强度尺寸效应

由不同尺寸的直剪试验求出其抗剪强度,得出抗剪强度τ与试样长度D的关系图,见图9。从图9可以看出,抗剪强度τ随试样长度D的变化而波动,说明试样的抗剪强度存在尺寸效应。在试样长度D为0～7 m时降低,7 m以后逐渐趋于稳定。计算结果局部出现波动,还需要具体判断指标来量化波动的大小以确定REV尺寸。设定差值变化率为量化指标,当差值变化率小于10%时确定为岩体的REV尺寸[34]。差值变化率计算公式为

(5)

图9 不同正应力条件下抗剪强度τ与试样长度D的关系曲线

式中:Ri为差值变化率,Pi为第i个试样力学参数值,Pi-1为第i-1个试样力学参数值。

抗剪强度的差值变化率见图10。由图10可以看出,正应力σ为2、3、4 MPa时,试样尺寸在11 m×11 m×11 m之后差值变化率均在10%以下,因此,确定裂隙岩体REV尺寸为11 m×11 m×11 m。

图10 差值变化率

3.2 抗剪强度参数尺寸效应

内摩擦系数和内聚力是表征岩体抗剪强度的重要参数。基于上述不同正应力下的抗剪强度,依据莫尔库伦理论求出不同尺寸试样抗剪强度参数,即内摩擦系数和内聚力。根据求得的内摩擦系数及内聚力分别做出其与试样尺寸的关系图,见图11。

图11 内摩擦系数f、内聚力c与试样长度D的关系

由图11可知,内摩擦系数f和内聚力c具有尺寸效应,且二者趋势不同:内摩擦系数f在试样长度D为0～4 m时逐渐增大,随后趋于稳定;在试样长度为0～6 m时,内聚力c逐渐减小,在9～11 m时亦有较大波动。同样采用差值变化率来量化波动的大小,以此作为确定REV的判断指标。内摩擦系数f及内聚力c的差值变化率见图10。由图10可知,内摩擦系数f在试样尺寸6 m×6 m×6 m以后的差值变化率均小于10%,因此,以内摩擦系数f为指标的REV尺寸为6 m×6 m×6 m。内聚力在试样尺寸11 m×11 m×11 m之后的试样的差值变化率均小于10%,判定以内聚力c为指标的REV尺寸为11 m×11 m×11 m。可见,以不同参数为指标确定的REV尺寸不尽相同,其中,依据内聚力和内摩擦系数求出的REV尺寸差别较大,内摩擦系数先达到稳定。

综合比较抗剪强度REV尺寸、内摩擦系数REV尺寸以及内聚力REV尺寸,选取11 m×11 m×11 m作为裂隙岩体的剪切特性REV尺寸。

3.3 剪切特性函数关系拟合

为了探究剪切特性力学参数之间与试样尺寸之间是否存在某种规律,对不同尺寸、不同正应力下的抗剪强度以及内聚力和内摩擦系数分别进行函数拟合。

非线性回归拟合迭代计算表明,抗剪强度、内聚力以及内摩擦系数与试样尺寸之间均存在指数函数关系。由于篇幅限制,在此对抗剪强度的分析仅以2 MPa正应力为例,抗剪强度、内聚力以及内摩擦系数与试样尺寸之间拟合的相关系数分别为0.85、0.89、0.89,拟合程度较好,如图12～14所示。拟合关系式分别为

图12 抗剪强度τ与试样长度D拟合

图13 内聚力c与试样长度D拟合

图14 内摩擦系数f与试样长度D拟合

τ=14.25+10.49×0.69D

(6)

c=11.39+12.25×0.64D

(7)

(8)

式中:τ为抗剪强度,c为内聚力,f为内摩擦系数,D为试样长度。

根据以上拟合关系式可知,当试样尺寸趋于无穷时,抗剪强度、内聚力和内摩擦系数均会趋于一个定值,符合尺寸效应最终会使力学参数达到稳定的特征,这也验证了剪切特性尺寸效应的存在和函数关系式的合理性。并且抗剪强度和内聚力与试样尺寸的关系均可由式y=a+b×cx来表达,其中,y为抗剪强度/内聚力,x为试样长度,a,b,c为常数。

4 结 论

本文以怒江松塔水电站坝肩岩体为研究对象,基于蒙特卡洛原理构建了岩体结构三维网络模型,利用PFC3D实现了不同尺寸等效结构岩体的直剪模拟试验,从剪切特性角度揭示了复杂结构岩体的尺寸效应。其主要结论如下:

1)研究区坝肩岩体的抗剪强度、内聚力、内摩擦角均存在明显的尺寸效应。试样尺寸在一定范围内时,抗剪强度、内聚力随试样尺寸的增大而增大,而内摩擦系数随试样尺寸的增大而减小。

2)基于差值变化率定量分析剪切力学参数的变化程度以确定REV尺寸,结果表明抗剪强度、内聚力、内摩擦系数的REV尺寸不完全一致,针对不同岩体工程应分别考虑其影响。为安全起见,本文综合确定研究区坝肩岩体的剪切特性REV尺寸为11 m×11 m×11 m。

3)通过非线性回归拟合分析,确定抗剪强度、内聚力、内摩擦系数与岩体试样尺寸之间近似存在指数函数关系,将为后续确定研究区坝基岩体抗剪强度参数提供参考。