含原生隐微裂隙岩石颗粒流模型构建及细观参数标定方法研究*

李志强 刘国锋 晏长根 董 凯 崔 洁

(①长安大学,公路学院,西安 710064,中国)(②海南大学,土木建筑工程学院,海口 570228,中国)

0 引 言

原生隐微裂隙是火成岩、沉积岩成岩过程的伴生物,在玄武岩、大理岩等硬、脆性岩石中十分常见。与工程尺度的断层或断裂、岩体统计尺度的节理相比,隐微裂隙尺度小,且分布隐蔽,无规律。以往建立岩石体力学模型或进行数值仿真时,较少考虑隐微裂隙的影响效应,实际上其对岩石的力学特性和破坏模式有着显著影响,具体表现为:工程扰动下,裂隙尖端局部应力急剧增大,加剧次生裂隙的扩展演化过程,导致岩石的宏观特征强度大幅降低(魏元龙等,2015; 唐礼忠等,2019; 张传庆等,2019); 且岩石的破坏模式与裂隙的分布特征如位置、角度、长度等密切相关(Xie et al.,2011; 袁广祥等,2019; 赵建军等,2019; 张杰等,2021)。在岩体工程领域,含裂隙岩石力学特性、破裂机理及裂隙作用机制的理论和应用研究一直是热点问题之一,目前常用的研究手段主要为室内试验与数值仿真,室内力学试验多以预制含裂隙试样为研究对象(Liu et al.,2018; Jiang et al.,2019; Song et al.,2019; 陈伟等,2021),然而预制裂隙工艺复杂、制作困难且难以准确反映真实裂隙的物理力学特性。以颗粒流分析程序PFC为代表的离散元仿真作为一种简单、便捷、行之有效的研究手段愈发受到研究者的青睐。

PFC能够从细观结构角度进行介质材料的力学特性的描述和受力变形的分析,可直观显现荷载作用下岩石内次生裂隙萌生、扩展直至贯通的全过程,在研究含裂隙岩石的力学特性和裂纹扩展特征方面具有明显优势(周喻等,2013; Itasca Consulting Group,Inc.2014)。陈鹏宇(2018)从接触模型选取、细观参数标定等方面详尽介绍了PFC2D模拟裂隙岩石的主流方法。陈秀云(2015)基于PFC构建了单一中心闭合裂隙岩石颗粒流模型,开展多工况数值实验,探究闭合裂隙特征参数,如宽度、长度、倾角等对裂隙岩石破坏特征的影响。黄丹等(2020),于辉等(2020)分别针对性地研究了闭合裂隙倾角与含裂隙大理岩、砂岩强度和变形特性的响应关系,其中于辉等(2020)同时考虑了围压对裂隙砂岩力学特性和破坏特征的影响。黄达等(2020)开展了岩石拉-压数值模拟的试验,实现了岩石高压拉强度比的离散元模拟,并分析了其破坏机制。邓清海等(2017)以直切槽式圆盘试样为研究对象,基于PFC进行多组数值巴西试验,探究裂隙倾角、长度对直切槽式圆盘试样裂纹扩展规律的影响。唐世斌等(2021)重复开展多次数值巴西劈裂试验,解析裂隙类型(张开、闭合、填充)、裂隙倾角等参数对次生裂隙演化过程的影响。李玉成等(2019),张晗等(2021)利用PFC开展含平行双裂隙岩样的数值压缩试验,揭示单轴加载条件下含平行双裂隙岩样的破坏规律及裂隙分布对该规律的影响。黄彦华等(2014)考虑岩桥倾角的差异构建多种含断续节理岩样颗粒流模型,并在不同围压下进行离散元仿真,分析围压和岩桥倾角对红砂岩力学特性和破坏模式的影响。何栋梁等(2019)建立了多个具有不同锚固角度的锚固节理颗粒流模型,并在不同法向荷载下进行直剪,探究锚固节理的力学特性及节理面的破坏特征。方前程等(2014)通过类岩石材料模型试验将研究对象尺度扩展为含断续节理岩体,基于PFC仿真,研究不同围压下含不同倾角断续节理岩体的力学响应。周喻等(2016)通过PFC构建断续节理岩质边坡颗粒流模型,分析工程尺度岩质边坡的细观破坏机制。上述有关裂隙岩石体的研究中,无论对象为小尺度的岩块或是大尺度的岩体,PFC仿真均可作为一项重要的研究手段。

与其他能够指定模型介质的宏观本构关系并赋予物理力学参数的数值软件不同,PFC采用颗粒局部接触和演化特征来反映宏观问题(李东东,2018),颗粒流模型的准确性及细观参数的合理性直接决定着模拟成功与否,然而多数研究者利用其进行裂隙岩石数值分析时,主流思路仍是以预设规则裂隙颗粒流岩石模型为研究对象,考虑原生隐微裂隙真实分布形态构建真实裂隙岩石颗粒流模型的研究不多见,此外含原生裂隙颗粒流模型所涉及的多个接触本构,其细观参数的标定流程同样尚不够明确。鉴于此,本文以含原生隐微裂隙岩石为研究对象,对其真实颗粒流模型的构建方法及仿真所需细观参数的标定流程展开系统研究,研究成果不仅极大增强了含原生隐微裂隙岩石颗粒流模型构建的准确性及细观参数取值的透明性和合理性,同时可为内置多接触本构颗粒流模型细观参数的标定提供重要参考。

1 隐微裂隙岩石建模及参数标定方法

根据适用维度的不同,离散元仿真软件PFC又细分为2D及3D版(Itasca Consulting Group,Inc.,2014)。相较于3D版,2D版计算效率高,硬件要求低,结果清晰、直观便于分析,在模拟裂隙岩石过程中更为常用,本文建模及参数确定方法主要针对PFC2D,但并不局限于2D版,结果同样可为3D版提供参考。

1.1 隐微裂隙岩石颗粒流模型构建

隐微裂隙岩石颗粒流模型的构建遵循“先构建完整颗粒模型,后导入数字化裂隙”的原则,首先建立与岩样同尺寸的闭合“墙体”,作为模型边界; 其后在“墙体”内部填充颗粒单元,运算平衡,完整颗粒流模型即构建完成。在此,需要明确颗粒最小粒径、最大最小粒径比、孔隙率、密度等颗粒几何参数。众多学者(Jensen et al.,1999; 周喻等,2011; 石崇等,2015; Chen,2016; 黄彦华等,2016; 陈鹏宇,2018)针对如何确定颗粒单元几何参数开展了大量工作,参照其研究结果同时结合室内物理力学测试即可确定上述参数。

完整模型建立成功并验证无误后,后续工序便是将隐微裂隙准确地导入其中,如何更精确对岩样内隐微裂隙的分布形态进行刻画显得尤为关键。朱珍德等(2007)综合利用电镜扫描、数字图像处理及数理统计等分析手段对大理岩内原生微裂隙的分布规律进行研究。钟志彬等(2017),汪权明等(2020)通过普通低倍放大镜分别对流纹岩、白云岩侧表面的微裂隙进行观察和描绘,借助CAD实现裂隙几何分布信息的量化。袁广祥等(2019)凭借偏光显微镜观测、拍照及数字化处理获取花岗岩内原生裂隙的分布规律。Jing et al.(2021)基于CT扫描和数字图像处理技术实现煤岩三维裂隙网络模型的重构及裂隙结构的定量化提取。以上诸多方法中,CT扫描能够在无损条件下获得岩样内隐微裂隙更为准确、全面的分布特征。利用该特性,提出图1所示含隐微裂隙岩石颗粒流模型构建方法。该方法包括:试样CT扫描,裂隙识别、切片导入CAD,坐标系(与完整模型相同)建立,裂隙临摹,几何坐标获取,fish裂隙生成、导入完整模型等步序。相对于拍照、临摹等方式所建立的含裂隙岩石离散元模型,上述方法所得模型更为准确。

图1 隐微裂隙岩石颗粒流模型构建方法

1.2 岩石及隐微裂隙细观参数标定

PFC软件内置多个接触模型用以描述介质颗粒间接触的力学行为,多位学者(Lee et al.,2011; Mas et al.,2011; 胡波等,2019)研究表明:平行黏结模型、平直节理模型可以很好地模拟完整岩石的力学特性,节理裂隙用光滑节理模型模拟较为合适,影响上述3种接触力学行为的主要细观参数分别如表1、表2和表3所示。

表1 平行黏结模型主要细观参数

表2 平直节理模型主要细观参数

表3 光滑节理模型主要细观参数

表中各细观参数并没有实际的物理意义(Mehranpour et al.,2017),标定其数值时常采用大量试算建立模拟结果(宏观力学参数)与所用细观参数的关系,依据试验所得宏观参数反推细观参数的方式。通过归纳、总结前人研究,提出了图2所示完整颗粒流模型细观参数的标定方法,平行黏结模型及平直节理模型均可用于描述完整模型颗粒单元接触间的力学行为,对于不同类型裂隙岩石,“平行黏结模型+光滑节理模型”和“平直节理模型+光滑节理模型”两种组合在描述裂隙与岩石的交互影响时各有千秋,可通过手动试算确定最佳接触模型,也可参考他人研究进行选择。整体上,完整模型细观参数的标定流程可概括为:模型选择、参数组合、批量模拟、敏感性分析、参数关系构建、参数反算和结果微调。

图2 完整颗粒流模型参数确定流程图

光滑节理模型模拟节理裂隙时,节理经过的颗粒,其原本的接触本构会被光滑节理接触替换,采用该方式模拟节理裂隙远比移除节理位置的接触或降低黏结强度的方式更为合理。结合前人研究,确定裂隙与颗粒单元间接触的细观参数时,首先需在考虑节理裂隙影响的条件下,分析光滑节理模型各细观参数与模拟所得宏观参数的响应关系。为此应先预建含不同倾角裂隙岩石的颗粒流模型,完整模型及参数上节已定保持不变,指定裂隙接触为光滑节理模型,设计多种参数组合方案,批量模拟构建结果样本,定性分析模拟结果与细观参数间的映射关系,随后根据CT扫描结果,建立能够反映裂隙真实分布形态的隐微裂隙岩石颗粒流模型,结合试验所得裂隙岩石力学参数及宏观结果与细观参数的关系,反复调算,即可确定裂隙颗粒流模型所需细观参数,详细流程见图3。

图3 隐微裂隙颗粒流模型参数确定流程图

2 方法应用与验证

隐晶质玄武岩是我国西南山区白鹤滩水电站地下厂房施工期间最广泛揭露的围岩之一,新鲜、完整、质地坚硬、岩性极好,见图4。然而厂房施工期间,隐晶质玄武岩洞壁却出现了大量片帮剥落、卸荷开裂等破坏现象,如图5所示。室内实验结果显示,隐晶质玄武岩的单轴强度呈现出极强的离散性,最高可达351MPa,最低仅为107MPa(江权等,2017)。刘国锋等(2016)发现其内部随机发育有大量分布隐蔽、无规律,细小闭合的隐裂隙。多位学者(江权等,2017; 张传庆等,2019)均认为正是由于原生隐裂隙的存在才导致部分隐晶质玄武岩试样强度极低,进而导致宏观强度表现出较高的离散性。

图4 隐裂隙发育的隐晶质玄武岩

图5 洞室施工期间隐晶质玄武岩典型破坏(片帮剥落)

以原生隐微裂隙发育的隐晶质玄武岩为例进行裂隙岩石颗粒流模型构建及参数标定方法的应用和验证。实验所用试样为国际岩石力学学会ISRM建议的标准圆柱样(直径50mm,高100mm),利用CT扫描测定内部隐裂隙分布形态,判定试样为完整试样或含隐裂隙试样,其后在RMT150-C刚性试验机上进行力学特性实验,采用位移控制模式,以0.002mm·s-1的速度进行加载,直至试样破坏。对于完整试样而言,当加载强度即将达到峰值强度之前,试样会断续发出“噼啪”声响,一旦加载达到试样破坏所需的峰值强度,试样会随即产生脆性爆裂破坏,伴随着巨大声响,炸裂成块状或片状并弹射出去。而对于隐含裂隙的试样而言,其破坏过程与完整试样则明显不同,加载过程较为平静,峰值破坏后,试样仍可保持相对较好的完整性。典型试样宏观力学参数及内部隐裂隙数量如表4所示。

表4 隐晶质玄武岩典型试样宏观力学参数及隐裂隙数量

2.1 完整模型构建及参数标定

通过指定“墙体”的端点和法线方向生成长100mm,宽50mm的矩形闭合“墙体”,采用“半径扩大法”在“墙体”内填充颗粒单元,调整内应力,删除悬浮颗粒,完整模型即构建完成。颗粒单元几何参数取值见表5。

表5 颗粒单元几何参数

颗粒分布满足预定要求后,即可按图2所示流程进行参数标定。首先指定颗粒间的接触本构,这里以平行黏结模型为例。关于平行黏结模型,众多学者(阿比尔的等,2018; 胡训健等,2021; Wu et al.,2021)已经开展了大量研究,其细观参数与宏观计算结果间的关系较为明晰,主要结论有:宏观弹性模量(E)对平行黏结有效模量(pb＿emod)最为敏感; 法向刚度(pb＿kn)和切向刚度(pb＿ks)的比值(pb＿kratio)对泊松比(ν)的影响最大; 而法向和切向黏结强度的比值(pb-ten/pb＿shear)则主要影响试样的破坏模式。保证结果可靠性的前提下,为简化平行黏结模型细观参数的标定过程,预先做出两个基本假设,即颗粒接触模量等于平行黏结有效模量,颗粒法向和切向的刚度比等于平行黏结刚度比,同时固定对计算结果影响不大的参数值,颗粒间的摩擦系数(fric)取0.5,黏结半径系数(pb＿rad)为1.0。

2.1.1 平行黏结有效模量(pb＿emod)的确定

分别取黏结有效模量为5GPa、10GPa、15GPa、20GPa、25GPa、30GPa和35GPa,固定其他参数的取值不变,重复进行7次单轴压缩模拟,模拟所得弹性模量与黏结有效模量的关系如图6所示。

图6 pb＿emod与宏观弹性模量的关系

弹性模量随黏结有效模量取值的增加而增大,根据表4,完整岩石的弹性模量为53.45GPa,初步确定平行黏结有效模量等于21.87GPa。

2.1.2 平行黏结刚度比(pb＿kratio)的确定

固定颗粒接触模量和黏结有效模量为21.87GPa,设置黏结刚度比为1.0、1.5、2.0、2.5、3.0、3.5和4.0,其他参数保持不变,7次重复试验所得泊松比随黏结刚度比的变化见图7。

图7 pb＿kratio与泊松比的关系

随着黏结刚度比的增加,泊松比也在不断增大,但是增大的速率在逐渐变缓。将表4所示完整试样的泊松比(0.23)带入图中对数函数拟合公式,所得黏结刚度比为2.06。

2.1.3 法向黏结强度(pb＿ten)和切向黏结强度(pb＿shear)的确定

黏结强度比(pb＿ten/pb＿shear)影响模型最终的破坏模式,确定黏结强度前,应先根据试样的破坏模式,确定黏结强度比。人为设定切向黏结强度为10MPa,取黏结强度比为0.2、0.6、1.0、1.4、1.8和2.2,其他参数不变,重复进行6次单轴压缩数值模拟,典型破坏形式见图8。

图8 不同黏结强度比pb＿ten/pb＿shear对应的破坏形式(蓝色为张拉裂隙,红色为剪切裂隙)

图8a为黏结强度比取0.2的模拟结果,典型的脆性劈裂破坏,破坏产生的裂隙主要为张拉裂隙,未出现剪切裂隙。随着黏结强度比的增大,剪切破坏产生的裂隙逐渐增多,试样的破坏形式也在由劈裂破坏向剪切破坏转换。玄武岩属于硬脆性岩石,完整试样破坏形式以脆性劈裂破坏为主,反复试算后,取黏结强度比为0.29,则对应的法向黏结强度为2.9MPa,设置黏结强度变化系数k为0.5、1、2、4、6、8和10,得到7组比值均为0.29的法向黏结强度和切向黏结强度,展开模拟,峰值强度与黏结强度变化系数的关系如图9所示。

图9 k与峰值强度的关系

峰值强度σ随着黏结强度变换系数k近似线性变化,依据表4所示完整试样的峰值强度(331.81MPa),黏结强度变化系数计算结果为13.69,则法向黏结强度和切向黏结强度分别为39.7MPa和136.9MPa。

以上述参数进行数值单轴压缩试验,所得力学参数与力学试验结果存在较大误差,见表6。

表6 初始参数模拟结果

上述问题是由于模型各细观参数对宏观计算结果存在交叉影响(石崇等,2018),即同一个细观参数影响多个宏观参数的计算结果,同一个宏观参数的计算结果受到多个细观参数的影响,需要结合已建立的宏、细观参数间的关系按照式(1)对已确定的初始细观参数进行微调。

(1)

式中:plt表示室内实验所得宏观参数;psm为数值模拟的结果;a代表宏、细观参数函数关系的比例系数;pbi指第一次确定的初始细观参数;pbc表示调整后的细观参数。

杜朗的徒弟陈洋接替了他在罗恬身边的位置——不光是要替代杜朗在富人圈子里表演魔术，还要替代杜朗做罗恬的情人。

调整后的参数如表7所示,该参数模拟结果见表8,此时模拟与实验结果的误差皆已控制在5%以内。数值模拟和室内试验所得应力-应变曲线及破坏形态的对比见图10,两者吻合度高,模拟较为成功。

表7 PFC数值仿真所用最终参数

表8 最终参数模拟结果

图10 数值模拟和室内试验结果的对比

2.2 原生隐微裂隙岩石模型构建及参数确定

裂隙试样中的裂隙“导入”过程如图11所示,该裂隙试样的力学参数及裂隙特征见表4。细观参数按图3所示流程进行标定,预设规则裂隙岩石颗粒流模型,其中预设裂隙长为1.5cm,倾角φf分别为0°、15°、30°、45°、60°、75°和90°,且预设裂隙位于模型中心位置。部分模型如图11所示,指定裂隙与颗粒单元间接触的数值本构为光滑节理模型,细观参数组合通过折减平行黏结模型对应细观参数及人为给定的方式生成(具体方案见表9),控制变量重复进行119组数值单轴压缩实验,模拟所得力学参数、设定的节理参数、裂隙倾角三者间的关系见图12。

表9 光滑节理模型参数组合方案

图11 裂隙导入流程

图12 含预设规则裂隙岩石颗粒流模型

对于裂隙倾角为90°的模型试样而言,改变法向、切向刚度和摩擦系数,模拟所得力学参数均未发生明显改变,表明其对图示3个细观参数并不敏感。其他裂隙试样的力学参数则在细观参数改变时均发生相应变化。

图13a、图13b、图13c所示分别为节理法向刚度对裂隙颗粒流模型峰值强度、弹性模量和泊松比的影响。当裂隙倾角为0°～30°时,随着节理法向刚度的减小,裂隙模型峰值强度明显逐渐降低。倾角为45°时,其峰值强度变化趋势同上,但变化值大幅减小。倾角为60°时,随法向刚度减小,峰值强度先增大后减小往复变化,倾角为75°时,模型试样的峰值强度恢复随法向刚度折减而减小的变化规律,但减小量较小。对于0°和15°倾角的裂隙模型,法向刚度依次折减,弹性模量则不断降低。倾角为30°时,弹性模量则先增加后降低。倾角为45°的试样,弹性模量随法向刚度减小而增加,且只有法向刚度减小到一定程度后,弹性模量才开始有明显增加。对于60°倾角裂隙模型,其弹性模量随法向刚度的变化规律恰好与45°模型相反,随着法向刚度的减小,弹性模量先减小后变大。75°倾角模型,弹性模量并未随法向刚度改变产生明显变化。对于0°～30°倾角的裂隙模型,节理法向刚度减小,其泊松比对应增加,且只有法向刚度减小到一定程度后,泊松比开始有明显增加。裂隙倾角为45°和60°时,随着法向刚度的减小,泊松比分别呈现先减小后增大与先增大后减小的相反规律。75°倾角试样,其泊松比不随法向刚度的变化而变化,这一点与弹性模量相似。

图13 模拟所得力学参数、设定的节理参数、预设裂隙倾角三者间的关系

节理切向刚度与裂隙模型峰值强度、弹性模量和泊松比的关系分别见图13d、图13e、图13f。当裂隙倾角为0°、60°和75°时,改变节理切向刚度,裂隙岩石的力学参数并未产生明显变化。对于15°和30°倾角的裂隙模型,降低节理切向刚度,峰值强度显著降低。含45°倾角裂隙模型,峰值强度同样随切向刚度降低,但降幅较小。对于倾角为15°的裂隙模型,改变节理切向刚度,模型弹性模量发生小范围变化。对于45°裂隙模型,其弹性模量随切向刚度减小呈“台阶形”增长,即切向刚度折减到一定程度后,弹性模量才会发生明显变化,泊松比变化规律类似。含30°倾角裂隙模型,当节理切向刚度折减1/50后,其泊松比显著增加。

图11所示为典型的含裂隙试样,其裂隙倾角分别为30°、45°和60°。对于30°倾角裂隙试样,降低节理法向、切向刚度,减小摩擦系数,可以降低其峰值强度。需要降低弹性模量时,则可通过大幅降低节理法向刚度和摩擦系数来实现。欲使泊松比增大时,则可采用减小节理法向、切向刚度及摩擦系数的方式。含45°倾角裂隙试样,增大节理法向刚度和摩擦系数,可以使试样的峰值强度增大; 减小节理法向、切向刚度可以增大试样弹性模量; 增大节理切向刚度和摩擦系数,可以减小裂隙试样的泊松比。对于含60°倾角裂隙试样,通过折减节理法向刚度的方式可以减小试样的峰值强度; 改变节理法向刚度和摩擦系数可以使试样的宏观弹性模量发生变化,但对应规律需要具体分析; 增加节理法向刚度可以增大试样的泊松比。

结合上述定性分析,经反复试算,适用于表4所示典型含裂隙试样的光滑节理接触细观参数如表10所示。考虑到隐裂隙为闭合裂隙,特增加了法向黏结强度和切向黏结强度两个参数,表征闭合裂隙的黏结作用。基于表10中参数及图11所示真实的裂隙岩石颗粒流模型进行数值单轴压缩实验,结果见表11,数值模拟与室内试验获得的力学参数十分接近。图14所示为模拟和试验所得应力-应变曲线的对比,两者十分吻合,可见,该方法所建模型及所确定的细观参数较为合理。

表10 仿真所用细观参数

表11 裂隙岩石颗粒流模型模拟结果

图14 模拟与试验结果的对比

3 主要结论

本文结合前人研究及CT扫描、图像数字化技术建立了原生隐微裂隙岩石颗粒流模型构建和细观参数标定的方法,并以隐晶质玄武岩为例,对方法进行应用和验证。主要结论如下:

(1)CT扫描配合图像数字化技术可实现试样内裂隙形态的精确刻画,基于此提出了能够反映裂隙真实分布形态的含隐微裂隙岩石颗粒流模型的精细化构建方法。

(2)平行黏结模型细观参数与宏观模拟结果(力学参数)的关系较为简明,但考虑各参数的交互影响,初始参数仍需微调后方可使用。光滑节理模型细观参数与宏观模拟结果的关系相对较为复杂,确定其参数时,应先根据已建模型定性分析不同倾角下宏观结果与细观参数间的响应。

(3)“平行黏结+光滑节理”的数值本构组合能够很好地描述含隐微裂隙隐晶质玄武岩颗粒流模型接触间的力学性质。

(4)通过典型案例,验证了本文所提出的含隐微裂隙岩石模型构建和参数标定方法,模拟与试验结果的对比分析证实了本文所提方法的合理性和可行性。