基于图像的低应力水平下岩石非线性变形模拟

《INTERNATIONAL JOURNAL OF MECHANICAL SCIENCES》：Image-based modelling of rock non-linear deformation under low stress levels

　　

摘要翻译：
本研究提出了一种基于图像的方法，用于模拟在围压增加条件下由孔隙闭合效应引起的岩石非线性变形，该方法考虑了低分辨率数字岩石图像中固有的部分体积效应。该方法首先重建数字岩石模型，其中每个体素被赋予介于0到1之间的孔隙分数值。随后，每个部分体积体素被视为一个两部分弹簧，其中柔顺的孔隙相先于较坚硬的固相发生变形。连续的岩石非线性变形过程被离散化为一系列线性有限元步骤，每一步由固定的宏观应变驱动。该方法的预测结果与三个砂岩样品上实验测量的应力-应变曲线和波速吻合良好。这些结果表明，所提出的方法不仅能够捕捉宏观非线性变形，还能直接提供孔隙尺度的洞察，揭示围压条件下控制岩石力学行为的微结构演化过程。

论文解读：

**一、研究背景与问题**

岩石的力学变形行为对于许多地下工程应用至关重要，尤其是在低应力水平（通常低于10-20 MPa）条件下，包括浅层油气开采、地热能源提取、地下二氧化碳封存以及地下土木工程设计等。从岩石压缩试验的应力-应变关系分析中可以发现，多孔岩石，特别是颗粒状砂岩，在达到线弹性阶段之前往往经历非线性变形阶段。然而，这一非线性阶段由于线性弹性理论的简化处理而常被忽略。岩石地球物理性质在这一非线性变形过程中可能发生显著变化，例如孔隙度和渗透率急剧下降，弹性模量和波速显著增加。这些变化主要归因于应力增加时预先存在的微裂纹和固体颗粒间隙的闭合。简单地以恒定弹性模量处理岩石无法真实反映其孔隙闭合行为及由此引起的地球物理性质变化。

**二、研究内容与结论**

研究人员提出了一种基于图像的方法来模拟低应力水平下由孔隙闭合效应引起的岩石非线性变形。该方法的核心创新在于将数字岩石图像中的部分体积效应纳入考虑，通过累积Beta分布函数重建包含亚分辨率孔隙空间的数字岩石模型，将每个部分体积体素视为包含柔顺孔隙相和坚硬固相的两部分弹簧模型，并通过变形正则化循环实现逐步加载模拟。研究使用两个Bentheimer砂岩样品和一个Clashach砂岩样品的实验数据对方法进行了验证。结果表明，该方法能够有效再现岩石在围压条件下的非线性变形特征，模拟的应力-应变曲线和波速与实验测量结果吻合良好。更重要的是，该方法能够提供孔隙尺度的微观结构演化洞察，揭示围压作用下微裂纹和微孔隙的渐进闭合过程，这一信息是宏观数学模型无法获得的。

**三、主要关键技术方法**

研究人员采用了以下关键技术方法：(1) 基于累积Beta分布函数的孔隙分数重建方法，用于从岩石CT扫描图像中计算亚分辨率孔隙分数；(2) 局部两部分弹簧模型，源自两部分胡克模型（TPHM），将每个部分体积体素中的孔隙空间视为先于固体相变形的柔顺相；(3) 变形正则化循环，包括变形后孔隙度重新计算、变形网格重新映射到规则网格以及更新模型供下一步加载使用；(4) 基于体素的三维有限元方法（FEM）求解逐步线性弹性问题；(5) 改进的Hashin–Shtrikman等效介质理论（EMT）用于计算波速。样本来源包括Bentheimer砂岩图像和实验波速数据，以及Clashach砂岩样品的压缩实验数据。

**四、研究结果**

**4.1 重建岩石模型**

通过Beta分布方法分割后，部分体积体素主要位于固体颗粒边界处，形成固体颗粒与纯孔隙相之间的过渡区域。这一区域物理上代表未分辨的微孔隙和微裂纹，在低应力下会闭合。与传统的多相阈值方法相比，Beta分布方法能够保留窄的颗粒界面，而传统方法会丢失许多这样的界面特征。研究将三个砂岩模型的部分体积相划分为500个子相，以实现从孔隙到固体的平滑过渡。

**4.2 高应力下的岩石变形**

模拟结果显示，随着单轴应力增加，窄孔隙（微裂纹和颗粒边界间隙）逐渐收缩并最终消失，导致颗粒直接接触。该孔隙闭合效应在相对较低的应力水平（如30 MPa）下尤为显著。变形后的模型与实验X射线显微CT图像定性吻合良好。采用Beta分布方法重建的模型表现出明显的非线性变形行为，而传统多相模型则呈现一致的线性响应。模拟应力-应变曲线与多个文献报道的实验结果对比表明，该方法成功再现了非线性阶段的特征。与Clashach砂岩实验数据的定量比较显示，整体均方根误差（RMSE）为4.45 MPa，在应变超过0.1%后相对误差降至10%以下。

**4.3 岩石变形过程中的波速**

研究人员在每个加载步后使用改进的Hashin–Shtrikman（HS）等效介质理论计算纵波和横波速度。当围压超过5 MPa时，纵波和横波速度预测的相对误差均低于5%；在5 MPa以下非线性最强区域，相对误差为10%–15%。自洽（SC）模型分析表明，纵横比为0.01的孔隙最能代表亚分辨率孔隙空间的特征。值得注意的是，低应力下测量波速接近HS下界预测，高应力下移向上界预测，这与微结构演化一致：低应力时颗粒被微孔隙分离，高应力时颗粒直接接触。

**4.4 图像分辨率和细化效应**

图像分辨率对模拟结果有显著影响。降低分辨率会改变图像灰度直方图，导致关键参数选择的不确定性增加。然而，如果从高分辨率图像获得可靠的参数集，即使在下采样模型上使用这些参数，仍能获得满意的应力-应变曲线。网格细化因子（5倍和10倍）对宏观响应的影响可以忽略不计。总孔隙度的累积相对误差保持在0.013%以下，平均体素体积误差在60 MPa时约为8%，表明数值守恒误差在可接受范围内。

**4.5 黏土矿物效应**

参数化研究表明，随着部分体积相中黏土含量的增加，应力-应变曲线下移，波速系统降低，表明岩石变得更柔软。在0.4%轴向应变下，黏土含量从0%增加到20%时，应力从约37 MPa降至约33 MPa，纵波速度从3.9 km/s降至约3.5 km/s。同时，非线性变形阶段随黏土含量增加而变得不那么明显。这表明忽略黏土矿物即使含量很小也会引入不可忽略的误差。

**五、讨论与结论**

在讨论部分，研究人员指出模拟初始阶段的较大误差源于CT图像无法分辨在加载初期闭合的最细微孔隙空间。随着应变增加，这些缺失的孔隙空间逐渐被消耗，机械响应主要由可分辨孔隙的闭合和固体骨架的弹性变形主导，因此误差减小。此外，自洽模型分析揭示了亚分辨率孔隙空间具有纵横比分布而非单一值，这为未来开发更精确的波速预测方法提供了方向。

研究结论可翻译为：本研究提出了一种基于图像的方法，用于模拟增加围压下由孔隙闭合效应引起的岩石非线性变形。通过将部分体积体素视为过渡相，并假设其中的孔隙空间先于固相变形（两部分弹簧假设），该方法捕捉了围压下微孔隙和微裂纹的渐进闭合过程，这是非线性刚度增加的主要驱动力。通过与三个砂岩样品的实验测量结果直接对比，该方法有效再现了实验室压缩试验中观察到的特征非线性变形阶段。该方法的一个关键能力是提供宏观经验或解析模型无法获得的孔隙尺度洞察。通过计算每个变形步后的孔隙度变化并正则化变形网格，该方法有效追踪了应力增加下孔隙结构的演化。模型可视化揭示了孔隙渐进闭合过程：较大孔隙收缩、颗粒直接接触、先前连通的孔隙变得连通性降低或孤立。这项工作的可能扩展方向是与降尺度技术和微连续介质理论集成，以实现耦合流体流动模拟和孔隙空间拓扑统计分析。总体而言，该方法为基于图像的力学建模提供了新途径，明确考虑了亚分辨率孔隙特征，为预测岩土和能源应用中的岩石变形提供了更现实的基础。