基于FLAC3D二次开发的含预制缺陷岩体数值模拟

周廷强，李江华

（1．重庆三峡学院 土木工程学院，重庆404100；2．河南工业职业技术学院 城市建设学院，河南 南阳473000）

经过漫长的地质作用以及各种工程扰动，岩体内部存在着各种节理、裂隙、孔隙、空洞等缺陷，在外部载荷的作用下，岩体的变形与破坏明显受控于岩体中的缺陷，呈现出明显的结构特征。同时缺陷的存在使得岩体的强度特征和变形破坏规律变得更加复杂，研究含缺陷岩体的失稳断裂规律对于保障缺陷结构工程的稳定性具有重要的理论价值和实践意义[1]。对此，国内外学者从室内试验到数值模拟做了大量的研究工作，其中通过室内试验还原缺陷岩体的结构特征、地质条件及工程力学环境模拟岩体的变形破坏过程和机制已经得到广泛应用。无论是单条[2]、双条[3]还是多条[4]预制缺陷的单轴抗压[5]、直剪[6]以及不同围压[7]条件下的室内试验研究均已得到开展。但室内试验建立的含缺陷试件比较复杂，且不可控因素较多，取样技术难，成本高，重复性差，这些都成为研究的瓶颈，且工程尺度岩体往往存在尺寸效应，室内试验尺度不能直接应用于工程尺度岩体，利用数值方法研究节理岩体力学性质成为目前很普遍的手段。Reyes 等[8]运用有限元方法研究了含2 条预制缺陷岩体试样的起裂、扩展和贯通模式，并利用物理试验结果对数值试验结果正确性进行了初步验证。指出适当的选用数值模型可以得到与模型试验相同的裂纹的起裂、扩展以及贯通模式。Bobet[9]等利用FROCK 程序建立了1 个适用于拉、剪复合破坏模式的准则的模型，并且研究了预制缺陷与外荷载形成不同角度时裂纹尖端的应力场。唐春安[10]基于细观损伤理论自行开发的RFPA数值计算软件，在细观层次上对岩石、混凝土等非均匀脆性材料的裂纹扩展规律进行了系统研究。梁正召等[11]采用细观损伤数值模拟方法模拟了单轴压缩条件下岩石试样的破坏过程，数值模拟得到了表面裂隙内部扩展、贯通过程，动态再现翼型裂纹、壳体裂纹的形态，探讨三维裂纹内部的受力机制，推测可能发生的断裂类型。上述研究从实验和数值模拟角度完善了缺陷岩体中强度变形规律以及裂纹扩展贯通原理，取得了积极有益的进展。但上述的试验与模拟研究中很少涉及到工程尺度上缺陷岩体，限制了研究成果的应用。为此提出了1 种弹脆性损伤本构模型，并利用C++语言进行二次开发将该模型嵌入到FLAC3D软件中，该方法不仅可以表征裂纹扩展规律，也可以得到裂纹扩展时微观上的拉剪破裂类型，而且基于FLAC3D软件的二次开发使得该方法适用于大变形下的岩石力学问题。作为案例分析，利用该模型对含不同倾角的双缺陷岩体进行数值模拟得到其强度及微裂纹扩展演化规律。

1 数值计算模型

1.1 库伦破坏准则与弹脆性本构关系

假定岩石的破坏主要是剪切破坏，岩石的强度，即抗摩擦强度等于岩石本身抗剪切摩擦的黏聚力和剪切面上法向应力产生的摩擦力。即平面上的剪切强度准则为：

式中：τ 为剪切面上的剪应力（剪切强度）；σ 为剪切面上的正应力；c 为黏聚力；φ 为内摩擦角。

在此描述的基础上，可以讨论岩石的破裂条件及其加载应力场的关系。在平面应力状态下，库伦准则可以用莫尔极限圆直观的图解描述，Mohr-Coulomb 破坏准则如图1。

确切的准则由强度曲线AL 表示，其斜率为f=tanφ，在τ 轴上截距为c。平面应力状态下应力σ 与τ由最主应力σ1与最小主应力σ3确定的应力圆所确定。

图1 Mohr-Coulomb 破坏准则

如果应力圆上的点落在强度曲线AL 之下，则说明该点表示的应力还没有达到材料的强度值，材料不发生破坏；如果应力圆上的点落在强度曲线AL之上，则说明该点的应力已超过材料的强度并发生破坏；如果应力圆上的点正好与强度曲线AL 相切（图1 中D 点），则说明材料处于极限平衡状态，材料所产生的剪切破坏将可能在该点所对应的剪切面上发生。

事实上，细观单元的破裂不仅仅只有剪切破坏，将细观单元破裂分为2 种，即拉伸破裂和剪切破裂。当应变状态达到了最大拉伸应变准则引起拉伸破裂；压应力或者剪应力导致应力状态满足剪切破坏准则，则发生剪切破裂。带拉伸破坏准则的库伦准则的表达式如下：

式中：σc为岩石的单轴抗压强度；σt为岩石的单轴抗拉强度

假定岩石单元在细观上的本构关系服从弹-脆性损伤模型（图2），在阶段I，岩石单元具有线弹性的简单应力-应变特性，该阶段为弹性段。在阶段II，岩石单元达到一定的强度准则后发生破裂转换为弱介质，则时候岩石单元具有残余强度特性。图中的εc和εr分别代表岩石单元达到强度准则时对应的应变，σcr与σtr分别代表岩石单元破裂后残余抗压与抗拉强度。

1.2 FLAC3D自定义本构开发

FLAC3D提供了1 种相对便利的开发环境，允许用户按照自己的意愿开发自定义本构。该本构利用C++开发后编译成DLL（动态链接库）文件，用户在需要使用该本构时只需调用即可。自定义本构的主要功能是返回产生相应应变增量后的新应力值，同时该模型必须提供其他的信息，如模型的名字及该模型中用到的材料参数，同时需要制定模型同软件交互的方式。文献中[12-13]已经提及基于C++的FLAC3D自定义本构开发的方法与流程，这里只对其2 个关键进行介绍。

图2 弹-脆性损伤本构关系

1）本构模型的注册。每1 个用户自定义本构模型被编译进1 个DLL 文件，该DLL 文件在FLAC3D中必须被实例化。按照惯例，在被用作FLAC3D插件的DLL 文件中包含4 个输出函数：getName、getMajorVersion、getMinorVersion 和creatInstance。此外还需提供1 个被称为DLLMain 的存根函数，当函数库从系统中加载或写在是会被调用。其中DLLMain 文件始终是相同的；输出的getName 始终返回1 个以“model”开头的字符串，且该字符串必须是独立的，以同其他本构模型区分开来。getMajorVersion 函数不应当被改变，表示本构模型的主要版本；getMinorVersion 表示的本构模型的小版本；creatInstance 函数实际上创建并返，1 个类别实体，它被存储于注册表中，用以创建其他的实体。

2）循环过程中本构与FLAC3D间信息的传递。FLAC3D与用户自定义本构间最重要的联系是被称为run 的成员函数，它计算出循环过程中模型的力学响应。1 个被称为State 的构架被用于与模型进行信息交流。成员函数run 的主要任务是计算新应变增量对应的新应力。在非线性模型中，与模型的内部状态交流信息也是非常重要的，以方便状态的输出。开发的弹脆性损伤模型需要输入的参数有弹性模量E，泊松比μ，密度ρ，黏聚力c，内摩擦角φ，抗压强度σc，抗拉强度σt，残余抗压强度σcr，残余抗拉强度σtr共9 个参数。

2 含预制缺陷脆性材料数值模拟

2.1 数值计算模型

利用弹脆性损伤本构模型（图2），研究不同缺陷角度对岩石的强度特性、应力变化与裂纹扩展规律及破坏模式的影响。数值模拟的模型示意图如图3，其中模型为50 mm×50 mm×100 mm 的长方体试件；缺陷角度为缺陷面与水平面的夹角，分别为0°、15°、30°、45°、60°、75°、90°共7 中工况；缺陷长度为预制缺陷与试样自由表面间的交线长度，为20 mm；裂纹厚度为预制缺陷与试样自由表面间的交线厚度，为2 mm；岩桥角度为两缺陷间最短连线与水平面的夹角，为45°；岩桥长度为两缺陷间最短连线的长度，为14 mm。

图3 数值计算模型

利用文献[14]中提及的网格划分方法，采用TCL脚本语言对SURPAC 进行二次开发，利用FLAC3D模型构建的数据转换方法，能够有效避免产生重复节点，从而实现FLAC3D模型的快速构建。数值计算网格如图4，划分网格均为1 mm×1 mm×1 mm 的正方体网格。

图4 数值计算网格

模拟砂岩选取力学参数[15]：弹性模量为2.0 GPa，泊松比2.5，密度2 500 kg/m3，黏聚力2.5 MPa，内摩擦角46°，单轴抗压强度为20 MPa，单轴抗拉强度为1.0 MPa，残余抗压强度为2.0 MPa，残余抗拉强度为0.1 MPa。采用两端位移加载方式，上下两端加载速率均为5×10-5mm/时步[16]。

2.2 数值计算结果

力学参数随缺陷角度的变化关系如图5，缺陷角度0°、15°、30°、45°、60°、75°、90°试件的峰值强度分别为13.75、13.45、12.74、11.23、9.45、9.73、12.27 MPa，很明显的可以看出随着倾角的增大，峰值强度成非线性变化趋势，从0°到60°,随着缺陷倾角的增大，峰值强度逐渐降低，从60°到90°,随着缺陷倾角的增大，峰值强度逐渐增加。于之相对应的弹性模量表现出与峰值强度类似的变化规律，随着缺陷倾角的增大，弹性模量先减小后增大，弹性模量分别为1.54、1.48、1.43、1.24、1.16、1.23、1.47 GPa，倾角 为60°的弹性模量最小。

图5 力学参数随缺陷角度的变化关系

最大主应力随缺陷角度的变化关系图如图6。当缺陷倾角为0°和15°时，首先在预制缺陷的两侧与岩桥区域产生一定的应力集中，岩桥迅速贯通，同时在缺陷尖端萌生拉应力引起的翼型裂纹，翼型裂纹沿着加载方向延伸一定程度不再扩展，在预制缺陷尖端产生新的次生裂纹，最终压剪裂纹贯通使试件发生破坏。缺陷倾角为30°和45°时破坏模式相似，2 条预制缺陷均以翼型裂纹形式扩展，应力集中出现在岩桥区域，使岩桥区域发生剪切线贯通破坏，两外端的裂纹继续沿加载方向扩展，最终由翼型裂纹贯通试件引起整体失稳破坏。缺陷倾角为60°和75°时，在预制缺陷尖端同时产生翼型和反翼型裂纹，岩桥区域出现1 条剪切裂纹，与两预制裂隙尖端的裂纹贯通使得试件发生宏观破坏。缺陷倾角为90°时，岩桥贯通不明显，最终由预制缺陷与试件上下表面贯通引发宏观破坏。

图6 最大主应力随缺陷角度的变化关系图

对试件破坏过程中的微破裂类型进行分析，塑性区随缺陷角度的变化关系如图7。

图7 塑性区随缺陷角度的变化关系

当缺陷角度为0°和15°时，拉伸破坏在整个破坏模式中相对其他角度有更大的占比，约40%，主要集中在缺陷的上下部分，该部分只有少量的剪切破坏掺杂。2 条缺陷之间区域为剪切贯通，拉破坏较少。当缺陷角度为30°和45°时，拉伸破坏相对于0°和15°时较为集中，主要分布于2 条缺陷上下侧，缺陷中间区域以剪切破坏为主，但相对于0°和15°时，缺陷内侧尖端区域以拉伸破坏为主，缺陷外侧尖端以混合型破坏和剪切破坏为主。整体上看，拉伸破坏占比相对0°和15°时有些微下降，约20%，和剪切破坏占比均有所增加。当缺陷角度为60°时，拉伸破坏区急剧缩减到5%，主要集中在缺陷内测靠尖端区域，整个试件明显开始以剪切破坏为主。当缺陷角度为75°和90°，拉破坏占比下降到1%以内，试件的破坏模式转变为剪切破坏。

3 结 语

1）通过C++平台FLAC3D进行二次开发，开发了1 种弹脆性损伤本构模型，该模型可以很好的反应岩石单元细观破坏。对其弹脆性本构进行详细描述，同时也展示FLAC3D二次开发的基本过程。

2）通过自主开发的FLAC3D弹脆性本构模型，对含不同倾角的预制缺陷岩体进行数值模拟计算，得到其强度随缺陷角度的增加先下降后上升的规律，同时其细观破裂由低缺陷角度时的拉伸-剪切型的混合破坏转变为高缺陷角度时的剪切破坏。