基于改进三维刚体弹簧法的横观各向同性岩石力学行为模拟及多尺度各向异性机理研究

本文针对横观各向同性岩石力学行为模拟的挑战，推荐一项结合三维改进刚体弹簧法（3D mRBSM）与新型各向异性工具箱的研究。该研究通过引入各向异性块体几何生成、基于三次贝塞尔曲线的弹簧参数定向分配及定向弹簧组模型（DSM），系统揭示了块体几何各向异性在控制材料力学响应中的关键作用。研究成果成功复现了Callovo-Oxfordian粘土岩和Tournemire页岩的实验观测特征，为复杂地质材料的多尺度力学建模提供了稳健且可扩展的离散方法框架。

　　

在地球表面广泛分布的沉积岩和变质岩，是地下工程、能源开采与储存、放射性废物处置等领域经常遇到的材料。这些岩石的力学行为往往表现出明显的各向异性，在大多数情况下可以简化为横观各向同性。横观各向同性岩石存在一个优选方向，沿该方向的力学性能呈现旋转对称性，无法用经典的各向同性力学理论完全解释。然而，用于模拟各向异性岩石变形和破裂的数值建模技术的发展仍值得进一步探索。岩石的横观各向同性力学行为源于多种尺度的地质特征，包括定向矿物组构、差异矿物沉积以及原生和次生结构面的存在。这种多尺度和非均质的各向异性起源给数值建模带来了巨大挑战，需要不同的各向异性技术来有效捕捉特定尺度的特征。

离散方法是模拟岩石变形和破裂过程的有效工具。这类方法的核心思想是将岩体离散为相互作用的颗粒或块体，从而能够模拟不连续现象，如裂纹萌生和扩展。基于不同的岩体离散化方法，常见的离散方法大致可分为：基于球体的方法（如PFC、YADE）和基于块体的方法（如UDEC、DDA、FDEM、RBSM）。基于块体的方法似乎特别适合表示粘结性岩石，而基于球体的方法通常用于颗粒材料。因此，本研究采用了基于块体的方法，具体是改进的刚体弹簧法（mRBSM）。mRBSM已被证明在模拟类岩石材料时表现出良好的特性，例如能够明确描述界面以及破裂和分离机制。

为了应对横观各向同性岩石形成的多样性，研究人员开发了一个用于三维改进刚体弹簧法（3D mRBSM）的各向异性工具箱，包含三项技术：各向异性块体几何（基于空间变换生成具有可控纵横比的块体）、各向异性弹簧物理参数（使用三次贝塞尔曲线捕捉弹性性能和破坏强度的方向变异性）以及定向弹簧组模型（DSM，用于纳入微观结构弱点）。通过系统的参数分析，研究了块体几何各向异性、参数取向各向异性和局部岩石组构对离散块体-弹簧系统力学行为的影响。开发了一个专用的建模工具箱来捕捉由方向相关机制引起的各向异性变形、强度和破坏模式。还提出了相应的参数标定程序。通过对两种代表性各向异性岩石的模拟，验证了模型的有效性和适用性。

为开展研究，作者主要采用了以下关键技术方法：基于空间线性变换的各向异性Voronoi镶嵌技术，用于生成具有可控几何特征的块体系统；基于三次贝塞尔曲线的弹簧参数方向分配函数，用于定义弹性刚度、强度等参数随界面方向的变化；定向弹簧组模型（DSM），通过调整选定界面上弹簧组的取向和参数来模拟微观结构弱点；以及改进的刚体弹簧法（3D mRBSM）计算框架，结合k-means聚类策略布置弹簧组，并采用基于内聚力模型（CZM）的渐进破坏算法模拟裂纹扩展。

2.1. 改进刚体弹簧法

该方法将研究域离散为块体系统的镶嵌结构，相邻块体通过弹簧组连接。每个弹簧组包含三个相互垂直的弹簧：一个垂直于界面，另外两个平行于界面。这些弹簧代表了相邻刚性块体之间界面的力学行为，能够模拟材料的变形、断裂起始和扩展。本研究采用Yao等人首先提出的改进刚体弹簧模型（mRBSM），该方法使用分布式弹簧组来捕捉由界面非均匀应力引起的渐进破坏过程。基于刚体假设，三维mRBSM模型为每个块体分配六个独立的自由度（DoF），包括三个平移未知量（u_x, u_y, u_z）和三个旋转未知量（ω_x, ω_y, ω_z）。块体的位移通过基于增量法的一系列线性方程求解。采用块雅可比预处理共轭梯度法（BJ-PCG）求解线性方程组。当发生非线性行为（如岩石破裂）时，采用改进的牛顿-拉夫森法作为收敛算法。

2.2. 3D mRBSM的各向异性技术

2.2.1. 各向异性块体几何

岩石固有各向异性的一个主要原因是矿物颗粒的择优取向。通常，RBSM模型使用泊松盘采样生成Voronoi网格，用于模拟各向同性类岩石材料的断裂过程。为了在刚体弹簧系统中产生各向异性，研究人员引入了一种基于空间变换的各向异性块体生成策略。该策略能够生成具有可控横观各向同性几何特征的块体系统，理论上可以创建具有任意纵横比的块体几何形状，其大小近似由M_ani量化。

2.2.2. 各向异性弹簧物理参数

在mRBSM模型中，所有物理参数都分配给相邻块体之间的弹簧组。对于准静态问题，这些参数通常分为两类：弹性相关参数和破坏准则相关参数。为了模拟各向异性材料，这些参数通常依赖于界面的方向。假设弹簧组的法向量与岩石的横观各向同性对称轴形成角度θ（范围从0°到90°），定义一个函数γ(θ)来描述不同θ值对应的弹簧参数变化。在三维mRBSM模型中，基于三次贝塞尔曲线定义了一个S形函数，该函数具有平滑性和易于进一步修改的理想特性。

2.2.3. 定向弹簧组模型

除了矿物几何各向异性，定向微观结构或弱点（如层理、叶理或带状结构）也是导致岩石各向异性的关键因素。当前的mRBSM模型在表示定向微观结构或弱点方面仍有局限性。受Cundall等人提出的平滑节理模型（SJM）的启发，研究人员通过人为调整弹簧组的方向，将一种各向异性方法引入3D mRBSM。由于与SJM不同，本研究将该方法称为定向弹簧组模型（DSM）以区别于前者。DSM的实现包括基于预定义规则选择需要应用DSM的界面，调整这些界面上弹簧组的取向，以及修改弹簧组参数以反映定向弹簧组的力学性能。

3. 各向异性工具箱的参数研究

3.1. 杨氏模量相关参数

基于初步模拟，发现影响模型杨氏模量的参数主要包括M_ani, E_0,ani, r_ani, ρ^D, E₀^D和r^D。通过一系列控制变量试验，研究了每个参数对杨氏模量的具体影响。结果表明，增加M_ani或降低E_0,ani都会降低E(0°)/E(90°)的比值，并使E(β)曲线呈现S形。然而，M_ani产生的影响更具中心对称性。当E₀^D < E(90°)时，增加DSC密度ρ^D会降低E(0°)/E(90°)的比值。此外，如果ρ^D固定为10%，降低DS弹簧刚度参数E₀^D会进一步降低E(0°)/E(90°)。r^D对E(β)曲线的形状有显著影响，但几乎不影响E(0°)/E(90°)的比值。

3.2. 强度相关参数

基于初步模拟，发现所有参数都影响横观各向同性mRBSM试样的强度。通过进行不同参数组合试样的单轴压缩试验和直接拉伸试验，分析了不同各向同性面倾角β下的单轴抗压强度（UCS）和直接抗拉强度（DTS）。结果表明，M_ani, E_0,ani和r_ani对UCS(β)曲线的特征有显著影响，即使它们不直接改变弹簧组的强度。增加M_ani, E_0,ani或r_ani会导致UCS(0°)升高。但它们对UCS(90°)的影响不同：增加E_0,ani导致UCS(90°)降低；增加r_ani导致UCS(90°)升高；而改变M_ani对UCS(90°)的影响很小。降低T_ani会导致UCS(90°)降低。降低S_ani会降低整个UCS(β)值，在β≈45°附近出现更显著的UCS下降。增加定向弹簧的密度（ρ^D）和弹簧刚度参数（E₀^D）显著增强了UCS(0°)，同时导致UCS(90°)略有下降。调整r^D可以使UCS(β)曲线从UCS(30°) < UCS(60°)转变为UCS(30°) > UCS(60°)。修改S^D可以影响UCS各向异性的程度，同时几乎保持UCS(0°)和UCS(90°)不变。T_ani和T^D对β接近0°时的DTS有显著影响，但对β接近90°时的DTS影响有限。

3.3. 破坏模式相关参数

横观各向同性岩石通常表现出与各向同性面取向和应力状态相关的破坏模式。因此，也有必要验证模型的宏观破坏模式。为了检验单个工具箱参数对数值试样破坏模式的影响，对单轴压缩试验（试验3.2-1至试验3.2-15）的峰后破坏模式进行了比较分析。结果表明，块体几何各向异性参数M_ani和剪切强度参数S_ani（用于一般弹簧组）和S^D（用于定向弹簧组）对试样的破裂面特征有显著影响。增加M_ani使得主破坏面倾向于与各向同性面对齐，并且在β介于30°和75°之间时，与M_ani=1.0的情况相比，观察到更明显的剪切破坏趋势。降低S_ani可以显著改善破坏面倾角与各向同性面倾角之间的相关性。S^D的效果与S_ani相似。然而，定向弹簧组不能完全取代块体几何各向异性的作用。

4. 模型应用

为了验证带有各向异性工具箱的新型3D mRBSM模型，研究人员对两种典型的横观各向同性岩石进行了模拟：Callovo-Oxfordian（COx）粘土岩和Tournemire页岩。选择这两种材料的原因是COx粘土岩具有相对完整的微观结构，而Tournemire页岩与COx粘土岩相比含有更多的定向微观缺陷。这种对比可以全面评估模型在低和高定向弹簧组影响条件下的性能。

4.1. Callovo-Oxfordian粘土岩

用于模拟COx粘土岩的微观参数中，定向弹簧组的密度ρ^D定义为0.5%的低水平。基于标定的参数，对各向同性面倾角β为0°、15°、30°、45°、60°、75°和90°的情况进行了数值模拟。数值模拟结果与实验数据在β相关的杨氏模量和三轴压缩峰值强度方面表现出良好的一致性。数值试样表现出与实验相似的β依赖性破坏模式：在β=0°时，破坏表现为拉-剪混合模式；随着β增加到45°左右，破坏逐渐转变为剪切主导模式；在β=75°时，破坏演变为网状劈裂模式。微观裂纹方向的分布也显示出与宏观破坏模式一致的β依赖性变化。

4.2. Tournemire页岩

用于模拟Tournemire页岩的微观参数中，ρ^D定义为5.5%，显著高于COx粘土岩的模拟。基于标定的参数，对各向同性面倾角β为0°、15°、30°、45°、60°、75°和90°的情况进行了数值模拟。数值模拟结果与实验数据在β相关的杨氏模量和不同围压下的三轴压缩峰值强度方面进行了比较。模拟结果和实验数据在两个各向异性指标k₁和k₂上表现出一致的趋势。具体来说，随着围压σ₃的增加，k₁增加，而k₂减小。数值试样的破坏模式与各向同性面的倾角β密切相关。在低围压（σ₃=1 MPa）下，β=0°、15°和30°的试样表现出拉-剪混合破坏模式。β=45°和60°的试样主要显示剪切破坏，但破裂角高于各向同性面的倾角。随着β向90°增加，破坏模式转变回拉-剪混合破坏。在高围压（σ₃=20 MPa）下，数值试样表现出比低围压下更强的β依赖性破坏模式。微观裂纹的方向分布表明，在低围压下，破坏模式以拉伸裂纹为主。然而，当围压增加到20 MPa时，剪切裂纹的比例显著增加。

5. 关于I型断裂能的讨论

断裂能已被广泛认为是控制离散方法中裂纹扩展的关键参数。3D mRBSM代码结合了内聚力模型（CZM）来表征弹簧破坏后的渐进软化和能量耗散。通过直接拉伸试验确定断裂能，在试样的顶部和底部施加位移控制载荷直至完全破坏。I型断裂能对试样的抗拉强度有显著影响。弹簧组参数各向异性和定向弹簧组工具都能在mRBSM框架内合理地表征各向异性断裂能。其中，前者倾向于产生高于预测值的G_soft，而后者产生的结果与理论预期一致。

本研究提出了一个用于模拟横观各向同性岩石各向异性力学行为的新型三维改进刚体弹簧法（3D mRBSM）框架。这项工作的核心目标是系统探索细长颗粒方法在离散各向异性建模中的作用，并评估其与其他显式和隐式各向异性技术的兼容性和潜在协同作用。为此，开发了一个专用的各向异性工具箱，以灵活控制弹性性能、强度和破坏机制中的方向依赖性。该工具箱集成了三个关键组成部分：各向异性块体生成、各向异性弹簧组参数分配和定向弹簧组模型（DSM）。尽管上述技术在离散方法领域并非原创，但本研究首次在刚性块体-弹簧框架内对它们进行了集成和系统分析。通过参数分析，全面评估了这些工具对关键力学响应的影响，包括方向相关的杨氏模量、压缩/拉伸强度和破坏模式。值得注意的是，块体的几何各向异性对所有讨论的评价指标都有显著的非线性影响。为了解决高维设置中参数识别的复杂性，设计了一种迭代标定策略，以解耦弹性和强度相关效应，并降低计算成本。通过对两种代表性各向异性岩石（Callovo-Oxfordian粘土岩和Tournemire页岩）的验证，证明了所提出框架的有效性，在弹性各向异性、方向强度变化和破坏模式方面与实验数据高度吻合。

总的来说，这项研究建立了一个全面且可扩展的建模框架，将微观几何特征与宏观各向异性力学响应联系起来。通过将细长颗粒方法与方向参数化和弹簧组配置相结合，所提出的方法能够在统一的数值方案中协调控制多种各向异性机制。这不仅推进了离散方法在模拟各向异性岩石方面的现有能力，而且证明了将各向异性几何技术与广泛使用的显式和隐式各向异性建模技术杂交的可行性。