岩石力学试验是研究岩石力学特性，获取力学参数的重要途径，其应用涉及土木建筑、采矿、水利水电、地下工程等众多工程领域。同时，随着现代力学、数学、计算机技术的不断发展，数值仿真以其模拟分析的有效性、计算结果的三维可视化等优点，为工程问题研究提供了有效手段^[1-2]。

近年来，关于岩石力学试验的数值模拟已成为岩石力学特性分析领域的研究热点。贾善坡等^[3]利用数值模拟方法建立泥岩试样二维简化模型，并根据非排水条件下泥岩室内三轴试验结果，反演研究岩石本构模型参数。田洪铭等^[4]结合泥质红砂岩的三轴蠕变试验结果，建立非线性蠕变损伤模型，并利用有限元软件建立与室内试验标准岩样尺寸一致的模型，对试验过程进行模拟。王杰等^[5]通过单轴压缩试验, 对岩石的破坏过程进行分析, 并通过数值模拟获得了裂纹萌生、扩展、贯通直至形成宏观裂缝的过程。王青元等^[6]采用所提出的模型对7个不同尺寸的绿砂岩岩样进行了单轴压缩蠕变数值模拟，并对岩石长期强度尺寸效应进行了分析。赵志涛等^[7]建立了砂卵石-砂土复合试样有限元数值模型，通过不同围压、不同砂层厚度的复合试样数值试验，探究砂层对复合试样强度和破坏规律的影响。刘洪磊等^[8]模拟了花岗岩单轴压缩试验过程，对其破裂机理进行研究。

通过以上文献可知，目前针对岩石力学试验的模拟主要侧重于岩石材料本构模型的二次开发或针对特定尺寸、特定加载方式模型的分析，而针对试验过程尤其是模拟分析的前处理过程(模型尺寸、加载参数的确定等)和后处理过程(模拟结果的输出、数据的采集)等方面的研究成果较少。文章利用Python语言对ABAQUS软件进行二次开发，建立岩石力学试验过程的参数化数值模拟仿真系统，将考虑剪切效应的Drucker-Prager破坏准则编入vusdfld子程序中，引入单元生死算法研究荷载作用下岩石试块变形破坏过程，实现了岩石力学试验过程的数值仿真，最后，以花岗岩试块三轴压缩试验为例，验证了仿真系统的有效性。

一、岩石力学试验数值仿真系统的开发 (一)基于Python语言的ABAQUS二次开发流程

ABAQUS大型有限元软件具有强大的工程问题分析能力，拥有丰富的模型材料库和单元库，能模拟复杂的非线性问题。同时，ABAQUS软件具备很强的通用性和扩展性，允许用户通过自行编程实现对软件前、后处理的二次开发。Python语言作为ABAQUS软件的标准设计语言，为用户自行编程，实现对ABAQUS的操控。

通过基于Python语言的ABAQUS脚本接口，可以直接调用软件丰富的库函数，并绕过用户界面直接操纵内核，实现工程建模、材料属性定义、模型网格划分、边界条件设置、结果数据查看和分析等功能，还可对用户界面进行修改和二次开发^[9]。

利用Python语言编程对ABAQUS进行前后处理的二次开发，实现如下功能：

(1) 前处理。通过编程实现用户自定义交互式模型尺寸、材料属性、边界条件等参数的输入，并通过有限元计算，输出不同参数下的计算结果文件。

(2) 后处理。通过编程实现对计算结果文件的读取、数据分析和用户指定结果的显示。

具体流程见图 1。

(二)系统开发

适用于模拟不同模型尺寸、材料属性、边界条件的岩石试件试验破坏及过程分析数值仿真系统，实现了不重复有限元分析时繁琐的前后处理操作，仅通过程序调用，在交互窗口中输入不同计算参数就可以自动进行建模、计算，并形成相对应的计算结果文件。岩石试件有限元模型如图 2所示。

岩石力学试验过程中，通常将试块设计为圆柱形，然后将其放置在加载系统上，通过电脑操作对试块进行加载，研究其力学特性。通常，试验参数主要包括试样几何参数(试样的尺寸)和加载参数(围压、轴压的量值)。为与实验室试验过程一致，岩石力学试验数值仿真系统开发的具体实施步骤包括：

(1) 将图 2中加载参数(轴压P₁、围压P₂)、试样几何参数(高H、直径D)以及试块力学参数(密度、弹性模量、泊松比、粘聚力、摩擦角)设定为可变化的参数，利用Python语言编程进行参数化建模^[10]。

(2) 在参数化建模程序中加上如图 3所示的交互式模块程序代码，使程序代码中的参数与试验参数相匹配，当程序调用后首先进入用户自定义交互式窗口，即可对数值模拟试验参数进行修改，如图 4所示。

(3) 对不同试验参数下的模型进行有限元计算。

二、岩石力学试验数值模拟 (一)单元破坏失效准则

岩石试样在压缩作用下，破坏模式主要包括剪切破坏、张拉破坏和拉剪混合破坏等^[11]。基于此，设置考虑剪切效应的Drucker-Prager准则(简称D-P准则)作为模型单元破坏的判断准则，以模型加载计算过程单元体是否达到以上状态作为判断其是否失效的依据^[12]。

D-P准则可表示为

$ f=\sqrt{J_{2}}+\alpha I_{1}-K=0 $

(1)

式中，I₁=σ₁+σ₂+σ₃，表示应力张量的第一不变量；$J_{2}=\frac{1}{6}\left[\left(\sigma_{1}-\sigma_{2}\right)^{2}+\left(\sigma_{2}-\sigma_{3}\right)^{2}+\left(\sigma_{3}-\sigma_{1}\right)^{2}\right] $，表示偏应力张量的第二不变量；α和K表示岩石的材料参数，可由以下公式计算：

$\alpha=\frac{2 \sin \varphi}{\sqrt{3}(3-\sin \varphi)}, K=\frac{6 C \cos \varphi}{\sqrt{3}(3-\sin \varphi)} $

式中，C和φ分别表示岩石的粘聚力和内摩擦角。

ABAQUS中vusdfld子程序可以读取模型迭代计算过程的单元状态(如应力、应变、等效塑性应变等)，将考虑剪切效应的Drucker-Prager破坏准则编入其中，当试块模型某单元应力状态满足D-P准则进入塑性状态时，认为该单元失去承载能力发生破坏，将该单元设定为失效单元，在下一加载步的迭代计算中将其删除。这样反复迭代计算，直至所有加载步计算结束。

(二)花岗岩试样室内试验分析

花岗岩质地坚硬致密、强度高，是一种较为常见的工程岩体。选取某工程现场取得的岩芯加工成φ50 mm×100 mm标准试样2块，在RMT多功能岩石试验机上进行三轴压缩试验(围压10 MPa)。由三轴压缩试验得出试样偏应力——轴向应变关系曲线，如图 5所示。花岗岩试样破坏如图 6所示。可以看出，三轴压缩条件下花岗岩试样的破坏形式较为一致，属于剪切破坏，可采用考虑剪切效应的D-P准则作为岩石破坏的判断准则。

对两个试块加载过程的试验数据进行整理，得到每个试样的物理力学参数如表 1所示。同时根据室内试验分析，得到岩石试块的粘聚力为27 MPa，摩擦角为45°。

(三)模拟结果分析

调用岩石力学试验数值仿真系统，利用图 4所示的用户自定义交互式窗口输入相应系数，对花岗岩试样破坏过程进行数值模拟。用户自定义交互式窗口参数输入如下：

(1) 加载参数设置。轴压P₁设定为25 MPa，围岩为10 MPa；

(2) 试样几何参数设置。试样高H设定为100 mm，直径D为50 mm；

(3) 试块力学参数设置。见表 1，粘聚力设定为27 MPa，摩擦角为45°。

通过岩石力学试验数值仿真系统对试块加载过程进行模拟，得到的模拟结果如图 7所示。可以看出，试块加载结束后，整体分为两部分，两者间存在较大的相对滑移，说明试块破坏形式为剪切破坏，这与试验得出的结论基本一致(图 6)。

同时在试样顶部取一关键点，通过自动读取该点在模拟过程的应力应变关系，得到模拟过程中试样应力-应变曲线，如图 8所示。可以看出，当模拟加载偏应力达到230 MPa左右时，试块失去承载能力，发生整体破坏，模拟过程中试样关键点应力-应变曲线与试验过程得到的曲线具有较好的一致性。

通过以上分析可知，岩石力学试验数值仿真系统计算结果与试验结果在试样破坏形式和应力-应变关系方面具有较好的一致性，这也验证了岩石力学试验数值仿真系统的合理性和有效性。

三、结语

(1) 利用Python语言对ABAQUS软件进行二次开发，对前处理过程(包括模型尺寸、加载参数的确定等)和后处理过程(模拟结果的输出、数据的采集)进行编程，建立岩石力学试验过程的参数化数值模拟仿真系统。

(2) 将考虑剪切效应的Drucker-Prager破坏准则编入vusdfld子程序中，引入失效单元删除算法研究荷载作用下岩石试块变形破坏过程，实现了岩石力学试验过程的数值仿真。以花岗岩试块三轴压缩试验为例，验证了仿真系统的合理性和有效性。

(3) 岩石力学试验数值仿真系统可根据用户需要与实验室加载平台进行数据匹配，极大提高了建模和分析效率，同时丰富了力学试验内容。