解锁中仿Rocscience RS2 与 RS3 的力量：利用Voronoi镶嵌模型(VTM)对硬岩破坏的真实模拟
预测硬岩在应力作用下的行为一直是一个长期挑战。与延性材料不同，硬岩通常以脆性方式破坏，在低约束条件下会突然且不可预测地断裂。
传统岩石建模的局限性
传统建模方法将岩石视为均质体，忽略了真实岩石中存在的颗粒尺度非均质性和微观结构缺陷。
这些微小缺陷会充当薄弱点，在突然破坏发生之前就已引发微裂纹。
因此，传统模型难以模拟脆性破坏的真实发展过程。
Voronoi镶嵌模型(VTM)：颗粒尺度的方法
为克服脆性破坏模拟的局限性，VTM用于捕捉硬脆性岩石中的颗粒尺度非均质性。该方法将岩石数字化划分为不重叠的多边形（或多面体）块体，每个块体代表矿物颗粒。块体之间的边界被建模为节理，模拟天然颗粒边界，从而允许在颗粒尺度上真实地模拟裂纹的萌生和扩展。
VTM类型：
非连续体VTM：能够显式模拟裂纹萌生与张开，但计算量大，尤其在三维中。
连续体VTM：计算效率更高，但过去局限于二维模拟。
在一次具有里程碑意义的合作中，Rocscience公司与达尔豪斯大学共同开发了基于RS2和RS3的三维连续体VTM。其目标是兼顾计算效率，并将颗粒尺度分析扩展到三维，从而克服二维连续体模型长期存在的局限性。
方法：模拟Lacdu Bonnet(LdB)花岗岩的脆性破坏
研究首次将三维连续体VTM应用于模拟LdB花岗岩的破坏过程。方法被总结为以下系统工作流程：
第一步——在RS2中构建二维Voronoi模型
RS2内置的Voronoi节理网络工具将模型域划分为代表矿物颗粒的多边形块体，同时在颗粒边界上定义节理单元。该过程天然引入了非均质性：不规则块体形状、不同节理长度与取向，以及块体与节理单元在刚度和强度上的差异。
RS2-VTM包含约515个中等规则的多边形，平均边长为3.5mm，并采用6节点三角形单元网格（见图1）。
通过迭代校准16个微观参数（块体与节理单元的强度与变形参数），使其与LdB花岗岩的宏观特性（峰值强度、杨氏模量）和宏观行为（破坏模式与峰后响应）一致。具体参数及范围详见表1。
无约束压缩试验模拟的边界条件：
下边界：竖向固定，中点在水平方向与竖直方向完全约束，以防止刚体运动并确保数值稳定。
侧边界：不约束（自由位移）。
上边界：施加竖向位移以模拟轴向加载。
表1：二维与三维VTM的微观参数（引自Amiri和Bahrani(2023)）
微观参数 | 块体（固体单元） | 块体边界（节理单元） |
E(弹性模量,GPa) | 78 | - |
ν(泊松比)υ | 0.22 | - |
kn(法向刚度,GPa/m) | - | 145,000 |
kn/ks(刚度比) | - | 2.5 |
σt(抗拉强度,MPa) | 22 | 9.8 |
c(黏聚力,MPa) | 82 | 45 |
φ(摩擦角,°) | 59 | 59 |
σtres(残余抗拉强度,MPa) | 0.1 | 0.1 |
cres(残余黏聚力,MPa) | 0.1 | 0.1 |
φres(残余摩擦角,°) | 59 | 59 |
第二步——开发平面应变RS3-VTM
经校准的RS2-VTM Voronoi节理网络被导入RS3，并沿y轴拉伸6mm以建立平面应变模型。施加的边界条件如下：
前/后边界：在y方向固定，施加平面应变约束。
上边界：施加竖向位移，模拟单轴加载。
下边界：竖向固定，中点全方向约束以防止刚体运动。
图1：用于模拟无约束压缩试验的Voronoi结构、网格单元与边界条件：（a）RS2-VTM；（b）平面应变RS3-VTM
通过在三维平面应变模拟中保持二维模型的微观参数与几何特征，本方法允许直接比较RS2与RS3模型的力学行为。
第三步——实现完整三维VTM
完整的三维RS3-VTM（60mm×60mm×120mm）由Neper的Voronoi镶嵌算法生成。该结构被导入RS3（也支持DFN节理网络导入），随后进行几何清理与节理面赋值。
图2：（a）Neper生成的三维Voronoi结构；（b）导入RS3后生成的完整三维VTM；（c）RS3-VTM的正视图，显示Voronoi块体
关键规格：
块体几何：约2000个Voronoi块体
网格配置：1,454,270个四节点四面体单元
运行优化：平均节理尺寸6mm（大于平面应变模型）
材料参数：第一步校准的微观参数
边界条件：下边界竖向固定且中心点全约束；上边界施加竖向位移以模拟单轴加载。
结果：对比分析与敏感性研究
二维与三维模型对比分析
应力–应变行为：RS2-VTM与平面应变RS3-VTM的应力–应变曲线几乎一致（见图3a），均表现为脆性响应。峰值强度分别为158MPa（RS2-VTM）与143MPa（RS3-VTM），差异源于网格单元排列不同。
破坏模式：两种模型均表现为以拉伸破坏为主：节理单元在约30%峰值应力时发生屈服，随后块体单元在接近峰值时屈服。在RS3-VTM中，峰值/峰后阶段出现少量剪切或混合模式屈服。
图3：RS2-VTM与平面应变RS3-VTM结果对比：（a）应力–应变曲线；（b）破坏模式
完整三维RS3-VTM性能
无约束压缩结果：
该模型的UCS为151MPa，破坏由峰前拉伸节理屈服引发，并在峰值/峰后阶段转变为剪切屈服。
图4：完整三维RS3-VTM无约束压缩试验：（a）应力–应变曲线；（b）不同加载阶段的破坏模式（节理与网格单元屈服）
敏感性分析
表2：网格单元密度对UCS与运行时间的影响
网格单元数 | 单轴抗压强度(MPa) | 运行时间(小时) |
~135,000 | 260 | 3 |
~240,000 | 216 | 4 |
~1,500,000 | 151 | 11 |
表3：块体抗拉强度对UCS的影响
主要发现：
网格细化显著影响UCS。随着网格密度增加，UCS降低42%（见图5a）。杨氏模量保持在70GPa，破坏模式变化较小（见图5b）。
块体抗拉强度是峰值强度的主要控制因素（见图6a）。
破坏模式始终以拉伸为主（见图6b）。
平面应变与完整三维模型均成功再现了物理试验中的颗粒尺度破坏过程。
图5：完整三维RS3-VTM无约束压缩试验：网格单元数量的影响：（a）应力–应变曲线与峰值强度；（b）破坏模式
图6：完整三维RS3-VTM无约束压缩试验：块体（网格单元）抗拉强度的影响：（a）应力–应变曲线与峰值强度；（b）破坏模式
结论
在RS2与RS3中集成VTM技术，代表了硬岩脆性破坏模拟的重大进展，并通过Lacdu Bonnet花岗岩案例得以验证。欲了解RS2与RS3的建模能力与应用，请联系我们。
联系我们
中仿智能科技（上海）股份有限公司
公司电话：021-80399555
地址：上海市松江区九新公路1005号临港松江科技城中仿大厦(201615)
电子邮箱：info@cntech.com
客服热线：400-888-5100
技术交流（QQ群）：534132040（GeoStudio），346550371（Rocscience），870585197（Leapfrog）
购买咨询（微信）：请用微信，扫一扫
