岩石是采矿、基坑和隧道等各类工程中的常见材料，对于各类工程的稳定性具有重要影响，因此研究其力学行为显得十分必要。岩石是一种非均匀脆性材料，在对其进行应力加载至峰值强度后，岩石材料的强度会迅速降低，但依旧会保有一定强度，这种由于变形引起岩石材料强度降低的现象称为岩石的应变软化（Wang et al，2020）。为准确模拟岩石的应变软化现象，数值软件通过特定的本构模型（刘志祥等，2023；Chen et al，2024）控制岩石力学参数，从而间接控制岩石的变形与破坏。然而，研究人员对岩体的本构关系和变形破坏机理的理解不够透彻，加之输入数值模型的数据参数不够准确，导致模拟结果的可靠性和准确性难以保证（刁心宏等，1999；唐春安，1999），因此构建精确岩石力学参数模型尤为关键。
在数值模拟中，构建真实的岩石应变软化模型是岩体工程计算分析的关键问题。数值模拟技术被广泛应用于揭示岩石峰后力学行为等方面，其突破主要体现在参数关联模型的构建和数值实现方法的创新上。通过在FLAC3D平台引入动态参数修正机制，实现了岩石本构行为与实验数据的精准映射。首先，建立了黏聚力和内摩擦角等关键参数与塑性应变梯度的非线性耦合方程，揭示了不同岩石在应变软化阶段的动态响应规律（张帆等，2008；王卫华等，2021）；其次，提出了广义力学参数框架，通过状态参数重构屈服面模型，显著提升了软弱岩体在复杂围压下的模拟精度（陆银龙等，2010；孙闯等，2015）。同时，基于Voronoi多边形离散模型（裴书锋等，2024）和预制裂隙试件试验（蒲成志等，2010；Liu et al，2016），研究发现矿物接触参数、裂隙倾角与裂隙密度三者共同控制岩石强度劣化路径，在高围压条件下，闭合裂隙更易诱发剪切破坏模式（Sun et al，2012），而张开裂隙主要控制拉张破坏模式（Eftekhari et al，2016；王雄等，2025）。在此背景下，通过改进FLAC3D弹脆性本构与超细网格划分技术，双裂隙岩石破坏模拟精度提升至90%以上（付金伟等，2012）。三轴压缩试验进一步表明，围压增大能够有效抑制岩石脆性破坏，且残余强度和初始残余应变随围压呈指数增长（唐海燕，2016；刘冬桥等，2017），据此建立的损伤本构模型（Chen et al，2023）通过损伤变量量化了裂隙发展对强度退化的贡献。在理论建模方面，部分学者创新性地引入包含退化角等特征量的新本构方程，定量阐明了围压对峰后脆性退化的影响（于永江等，2012），以莫尔—库伦准则为框架，采用内摩擦角和黏聚力的应变相关函数构建峰后应力—应变关系（李文婷，2012），结合位移法解析软化段斜率（陈康等，2013），最终实现了从理论推导到有限元程序计算的闭环验证（王水林等，2014；Zhang et al，2022）。总体而言，关于岩石力学行为的研究范式已从单纯参数调整转向本构方程迭代优化，基于Fish语言的二次开发极大地增强了参数修正系统的自适应性（王嵩等，2017），使得不同尺度模拟结果与试验数据的匹配误差普遍控制在5%以内（王小平等，2018）。
1 数值软件
1.1 FLAC3D软件
FLAC3D软件具备三维复杂地质建模能力，支持弹性、塑性、流变多相材料模型及参数自定义，通过显式拉格朗日算法与混合离散技术，高效模拟岩土材料的塑性破坏、流动及接触滑移过程，为岩土工程提供高精度数值分析，因此被广泛用于岩土工程和采矿工程等领域的数值模拟过程中。其中，早期应用主要是简单的地下开挖和边坡稳定性分析，后来扩展至多场耦合分析，再到现在的智能算法结合。本研究基于FLAC3D软件三维力学响应特征与高效解算优势，开展岩石单轴压缩数值模拟研究，精准解析岩石的应变软化与结构破坏特性。
1.2 Python软件
Python是一种高级的编程语言，其语法清晰、简洁且易于学习，采用缩进表示代码块，使得代码结构更加直观。此外，Python不需要预先声明变量类型，增强了灵活性和易用性，无论是初学者还是专业开发者，Python都是完成项目的理想选择，其应用领域极其广泛，包括Web开发、数据科学、人工智能和科学计算等。因此，选择Python作为数据处理和程序循环的工具。
2 多参数逆向优化理论框架
2.1 莫尔—库伦准则
在FLAC3D中，莫尔—库伦准则与内摩擦角、黏聚力以及最大主应力和最小主应力相关，如下为破坏准则，当 =0时，受力岩石处于临界状态。表示为
式中： 为最小主应力； 为最大主应力； 为内摩擦角；c为黏聚力； 。
2.2 塑形流动法则
在FLAC3D中，莫尔—库伦模型最后计算的应力，是以流动法则为基础，最后经过塑形应力修正得到，修正函数可表示为
式中： 为塑形应力修正值； 为流动法则系数； 为内摩擦角； ； ； ； 为体积模量； 为剪切模量； 为关于剪胀角的常数； ， 为材料常数。
2.3 公式推导
通过联立式（1） 、式（2） 和式（3） ，得到一个关于内摩擦角、黏聚力和剪胀角的多元函数，该函数表示FLAC3D中应力塑形修正后的值，通过对该函数的各元素分别求偏导，最后得出各元素对函数值的影响，以便后续写入程序。计算公式为
令 ，可得剪胀角相关值：
由上式可知，剪胀角与塑形修正值呈负相关关系，即可以通过减小剪胀角来增大塑形修正函数值，或通过增大剪胀角来减小塑形修正函数值。同理可以得出：
由以上分析可知，岩石峰后的强度与剪涨角、内摩擦力和黏聚力密切相关，FLAC3D可通过岩石的本构关系，以多参数逆向优化理论为基础，控制岩石这3个力学参数，从而间接控制岩石的变形与破坏。
3 应变梯度自适应模型实现
3.1 模型构建
本文模拟岩石单轴压缩试验，利用FLAC3D中的Extrusion模块进行建模，先建立平面图形，再添加厚度，形成实体单元。然后，采用实体单元模拟受压岩体，模型尺寸长度为50 mm，宽度为50 mm，高度为100 mm。整个模型有2 000个实体单元，2 541个网格节点，模型结构如图1所示。
3.2 边界条件的设定
对于模拟单轴压缩试验，只需固定z=0平面所有节点z向速度，四周不设定约束条件，为自由面，而顶部z=100平面设定向下恒定的速率，以模拟试验中的位移加载。
3.3 采用的本构模型与计算参数
由于岩石的结构和力学特征非常复杂，在外力作用下，岩石不仅会发生弹性变形，还会产生塑性变形。为了分析本模型，采取以下假设：（1）连续性：假设岩石为连续的介质，不考虑内部裂缝或缺陷的影响；（2）均质性：假设岩石在整体上具有均匀性，各部分的物理和力学特性相同；（3）各向同性：假设岩石在各个方向上的力学特性相同，不考虑材料的方向效应；（4）忽略构造应力场，不考虑由于地质构造导致的应力场变化，以简化模型分析。本研究采用莫尔—库伦强度准则，这是岩体工程中常用的强度准则，能够有效预测岩石在受力情况下的破坏行为。通过上述假设和准则，可以清晰地分析岩土体在外力作用下的力学响应。
表1 岩石力学参数Table 1 Rock mechanics parameters |
| 参数名称 | 参数值 | 参数名称 | 参数值 |
|---|---|---|---|
| 弹性模量/GPa | 1.4 | 黏聚力/MPa | 2.5 |
| 泊松比 | 0.28 | 剪胀角/（°） | 10 |
| 内摩擦角/（°） | 40 |
3.4 模型加载
4 实时动态修正算法构建
4.1 Python与FLAC3D的交互
在FLAC3D中已有与Python联合使用的例子，FLAC3D中的Fish语言在Python中是以字符串（string）的形式写入，最后导入FLAC3D中，从而实现Python语言与Fish语言的联合使用。
4.2 修正模型的程序实现
在运行程序前需要对已有试验数据进行提取并处理。本研究利用Python强大的数据处理能力，将试验数据提取出来，并将数据放入list容器，命名为lines，随后调用FLAC3D关于各参数的内嵌命名，构建实时动态修正算法，最后将整个程序输入FLAC3D内置的Python的平台并运行，实现对模型参数的动态修正。单轴压缩修正模型程序实现过程如图3所示，以运算到第n步为例，对流程图加以说明。
（1）运算到第n步，通过函数s.stress（）［2］［2］提取此时模型中的应力，且与试验数据lines［a］进行作差，将得到的值进行绝对值处理，如果这个值小于0.01，说明模拟值与实际值相差不大，不用对参数进行修正，则直接进入下一个循环，即n=n+10。
（2）如果s.stress（）［2］［2］与lines［a］二者相差的绝对值大于0.01，说明模拟值与实际值存在较大差异，需要进行参数修正。
（3）如果s.stress（）［2］［2］>lines［a］，则模拟应力值过高，根据前文计算结果，需要增大剪胀角，减小黏聚力和内摩擦角，修正系数为0.01。如果s.stress（）［2］［2］<lines［a］，则模拟应力值过低，此时根据前文计算结果，需要减小剪胀角，增大黏聚力和内摩擦角，修正系数为0.01（修正系数是根据试验数据精度选取，试验数据精度提取到百分位，因此设置为0.01）。进入下一循环，即n=n+10。
（4）每10个时步进行一次循环参数，直至规定时步结束。
表2 部分内嵌函数与命令Table 2 Partially embedded functions and commands |
| 名称 | 用途 |
|---|---|
| s=it.zone.find（） | 获得指定id号的单元 |
| s.prop（'dilation'） | 获得单元剪胀角 |
| s.prop（'friction'） | 获得单元内摩擦角 |
| s.prop（'cohesion'） | 获得单元黏聚力 |
5 单轴压缩数值试验分析
利用Python建立试验数据驱动的实时反向传播算法，将单轴压缩修正模型程序嵌入到FLAC3D内置的本构模型中，实现莫尔—库伦本构参数的多线程动态反馈迭代。
5.1 时步应力分析
通过FLAC3D内嵌命令，提取修正模型各时步瞬时应力，具体数据见表3。
表3 时步瞬时应力Table 3 Instantaneous stress at time step |
| 时步 /（×10³） | 原模型 应力/MPa | 修正模型 应力/MPa | 实验应力/MPa | 原模型相对误差/% | 修正模型相对误差/% |
|---|---|---|---|---|---|
| 平均误差 | 84.18 | 3.89 | |||
| 1.20 | 11.90 | 8.51 | 8.88 | 34.01 | 4.19 |
| 1.40 | 11.70 | 7.86 | 7.39 | 58.32 | 6.40 |
| 1.60 | 11.00 | 6.56 | 6.50 | 69.23 | 1.00 |
| 1.80 | 11.60 | 4.98 | 5.25 | 120.95 | 5.08 |
| 2.00 | 11.30 | 4.61 | 4.74 | 138.40 | 2.80 |
由表3可知，在各时步修正模型的相对误差均小于原模型。尤其当时步为2×103时，原模型的相对误差达到138.4%，而修正模型的相对误差仅为2.8%，后者较前者精度提升了49.4倍。此外，原模型平均误差为84.18%，而修正模型为3.89%，后者较前者精度提升了21.6倍，说明该修正模型较好地调整了岩石的力学参数，更好地描述了岩石峰后的力学行为。
5.2 应力—位移曲线分析
岩石试件在单轴压缩荷载下的变形全过程可通过图4中的试验曲线表示。根据全应力—位移曲线，岩石的变形可划分为4个阶段。（1）裂隙压密阶段：试件内部的裂隙逐渐闭合，岩石被逐渐压密，表现为早期非线性变形，应力—位移曲线呈上凹型。（2）弹性变形至微弹性裂隙稳定发展阶段：试件处于弹性变化阶段，岩石展现出线性弹性特性，应力与应变成正比，符合胡克定律，应力—应变曲线呈近似直线。（3）非稳定破裂发展阶段：岩石从弹性转向塑性的关键点，岩石承载能力达到最大值，出现非稳定裂隙，并逐渐扩展。（4）破裂后阶段：岩块承载力达到峰值强度后，其内部结构遭到破坏，但试件基本保持整体状。到本阶段，裂隙快速发展，交叉且相互联合形成宏观断裂面。此后，岩块变形主要表现为沿宏观断裂面的块体滑移，试件承载力随变形的增大而迅速下降，但并不会降至0，说明破裂的岩石仍有一定的承载力。
原模型的力学参数为定值，采用FLAC3D中内置的莫尔—库伦本构模型。由图4（a）可知，修正模型模拟曲线与原模型模拟曲线在峰值前的曲线基本一致，但是峰后曲线相差较大。原模型在达到峰值强度后，其应力—应变曲线未能反映岩石内部损伤扩展与应力下降的趋势，而是持续维持峰值强度水平。这主要是由于模型中的关键岩石参数（如剪胀角、内摩擦角和黏聚力）在峰值后阶段未发生相应演化，仍保持为初始定值，从而导致模拟结果与岩石单轴试验中观察到的峰后应变软化行为存在显著差异，难以符合实际岩石的力学响应。因此，有必要对该本构模型进行合理修正。
修正模型是以FLAC3D中内置的莫尔—库伦本构模型为基础，利用Python将多参数逆向优化理论写入程序，对岩石峰后的力学参数进行动态调整，从而实现对原本构模型的修正。由图4（b）可知，经过程序修正的模型，岩石峰后的应力随位移的增加而逐渐降低，且峰后应力—位移曲线与试验曲线基本一致，较为准确地还原了单轴压缩试验应力—位移曲线，较好地描述单轴压缩条件下岩石峰后力学特性。峰后曲线的相似性表明，岩石在破坏后的残余强度和应变软化等行为得到合理反映。通过精准模拟岩石峰后力学行为，可以在确保工程安全的情况下，充分利用岩石的峰后强度，从而创造更多价值，也可以进行风险评估，为工程决策提供科学依据，降低灾害的发生概率。
6 结论
（1）创新性地提出利用Python建立试验数据驱动的实时反向传播算法，通过FLAC3D内嵌接口实现莫尔—库伦本构参数的多线程动态反馈迭代，得到模拟曲线与试验曲线较为吻合，高估趋势完全消除。
（2）以FLAC3D内置的莫尔—库伦本构模型为基础，结合Python进行程序开发，较为准确地修正了岩石的力学参数，且易于实现，能够广泛应用在工程实践中。
（3）提出了一种岩石力学参数优化方法，准确再现了岩体峰后阶段的力学行为演化过程，为岩石关键参数精度缺失的优化提供了新的思路。
甘公网安备 62010202000672号
