H-M-C条件下模拟溶质运移和矿物溶解的LBM-FDEM耦合框架

《Water Resources Research》：A Coupled LBM-FDEM Framework for Simulating Solute Transport and Mineral Dissolution Under H-M-C Conditions

　　

将CO2注入地质储层会使储层环境酸化，导致矿物颗粒溶解，从而增加CO2泄漏的风险。因此，深入理解这些微尺度溶解机制对于地质CO2封存的安全性至关重要。为应对这一挑战，研究人员开发了一种多场耦合模拟框架，该框架将格子玻尔兹曼方法（LBM）与有限-离散元方法（FDEM）相结合，用于模拟溶质运移和矿物微观结构溶解。该框架引入了几何映射检测方案和改进的动量交换算法，以捕捉流体-固体相互作用。此外，还提出了一种新的溶解速率和溶质浓度更新方案，明确考虑了质量守恒和复杂矿物形态。同时，还提出了一种动态几何拓扑更新算法，用于追踪矿物颗粒的形态演化。通过多个数值基准测试（包括颗粒沉降、运动颗粒的溶质释放、多晶型颗粒溶解以及非均匀通道内的反应性运移）验证了LBM-FDEM框架的准确性和能力，结果与现有解决方案高度一致。最后，将该方法应用于研究不同Péclet数下裂缝溶解的时空演化。数值结果表明，在扩散主导的低Péclet数条件下，产物积累将溶解限制在上游入口处，呈现出典型的表面溶解模式。相反，在高Péclet数条件下，强对流驱动反应物快速渗透，引发沿整个裂缝的均匀通道化拓宽。最终，该框架在模拟地下工程中复杂的流体-力学-化学耦合方面展现出显著潜力。

随着全球碳中和战略的推进，CO2地质封存（碳捕集、利用与封存，CCUS）被视为减缓全球变暖的关键技术之一。然而，CO2注入储层后会酸化储层环境，导致矿物胶结物溶解，进而改变岩石微观结构。这些微观结构变化可能引发储层内应力重分布，诱发微裂纹的萌生或扩展，最终降低储层完整性并增加CO2泄漏的风险。这一渐进劣化过程本质上涉及高度复杂的孔隙尺度流体-力学-化学（H-M-C）耦合过程，涵盖反应性溶质运移和矿物溶解。虽然大量实验研究已显著推进了对这些过程的理解，但受限于当前实验观测技术在实时性和视场方面的局限，仅通过实验难以完整揭示其演化机理。

为突破上述瓶颈，数值模拟方法凭借其可视性、可重复性和经济性优势，已成为探索复杂H-M-C耦合环境下矿物溶解和微观结构演化的重要替代手段。在现有数值方法中，连续-非连续耦合方法（如有限-离散元方法，FDEM）在模拟岩石变形和破裂方面展现出最大潜力，因为它无需额外假设或重新网格化即可同时模拟岩石破坏过程中的连续-非连续转变。然而，现有基于FDEM的流体-力学模型依赖岩石裂缝中理想化的层流假设，难以模拟复杂H-M-C耦合环境下的矿物溶解和微观结构演化。为弥补这一缺陷，需要引入更先进的流体模拟方法替代当前的FDEM流体-力学耦合模块。格子玻尔兹曼方法（LBM）通过将流体域离散为规则正交格点并利用粒子分布函数求解玻尔兹曼方程，无需追踪单个流体粒子，显著提高了计算效率，且其局部的碰撞和迁移操作天然适应复杂流体边界。尽管已有研究将LBM与FEM或DEM耦合以捕捉流体-结构相互作用，但在大多数现有反应性运移模型中，流体-固体界面在固定正交网格上持续更新，固有地产生数值“阶梯”效应，损害了溶解过程中多角度矿物边界连续动态后退的捕捉精度。因此，开发一种完全耦合的LBM-FDEM框架以消除阶梯效应并显式追踪反应性运移过程中复杂边界的形态演化，具有迫切的必要性，但目前尚无相关报道。

针对上述问题，研究人员开发了一种耦合LBM与FDEM的多场耦合模拟方法，用于模拟复杂孔隙-裂缝岩石介质中的矿物溶解和溶质运移。在构建该耦合框架时，首先提出了一种有效的几何映射检测方法，用于识别被三角形单元覆盖的LBM格点并建立二者间的几何映射关系。在此基础上，引入改进的动量交换算法计算流体与矿物颗粒边界间的相互作用力。为模拟化学反应诱导的流体-矿物相互作用（即矿物溶解），进一步开发了H-M-C耦合框架，其中基于质量守恒方程开发了溶解速率更新算法，以提升溶质浓度增量和矿物颗粒质量损失的预测精度。为更精确地模拟溶解对矿物颗粒几何形状的影响，基于更新的质量损失提出了FDEM三角形单元的拓扑更新算法。通过上述方法，LBM-FDEM耦合方法被应用于多个数值基准测试：首先模拟粘性流体中的颗粒沉降过程，并与解析解对比以验证H-M框架的准确性；随后模拟流动流体中圆形颗粒的溶质释放过程以验证H-C框架的准确性；为突出所提出的H-M-C框架在捕捉复杂形态颗粒非均匀形状变化方面的优势，还进行了不同多晶型颗粒的完整溶解过程模拟；此外，通过在非均匀通道内模拟可溶颗粒的反应性运移，并与已建立的LBM-FVM解进行定量对比，验证了颗粒溶解模拟的准确性。在此基础上，以相同几何尺寸的孔隙尺度微流控实验为参照，模拟了单一裂缝通道内的反应性运移并将所得形态模式与实验结果直接对比，从而验证了模型的宏观能力。基于上述验证，该框架被应用于系统研究不同Péclet数下裂缝溶解和离子运移的时空演化。该项研究的成果发表于《Water Resources Research》，为量化地下工程中复杂的H-M-C过程提供了稳健的数值工具，对保障CO2封存安全性具有重要意义。

**主要技术方法**：为开展上述研究，研究人员采用了以下关键技术方法：（1）几何映射检测方法，通过三角形面积坐标判据精确识别流体、边界和固体节点，并采用加权外推方案重构新生成边界节点的分布函数；（2）改进的动量交换算法，基于流固交界面的精确交点计算流体动力节点力，并等效映射至三角形单元节点；（3）溶解速率与溶质浓度更新算法，基于质量守恒方程将颗粒质量损失与边界溶质浓度相关联，采用有限差分格式结合广义浓度边界条件（可涵盖Dirichlet、Neumann及Robin边界）求解边界浓度；（4）FDEM三角形单元拓扑更新算法，根据元素边缘暴露于流体的程度将边界单元划分为四种拓扑类型，分别采用共享节点平均策略、质心收缩策略及各向同性缩放方法更新节点坐标。数值验证中，LBM-FDEM模拟的颗粒沉降过程与Wan和Turek（2006）及Jiao等（2021）的基准结果吻合，移动颗粒溶质释放与Xia等（2023）的基准结果一致，多晶型颗粒溶解形态演化与Derksen等（2015）、Bouissonnié等（2021）及Adobes-Vidal等（2017）的实验观察相符，非均匀通道内可溶颗粒溶解与Gray等（2021）的LBM-FVM解高度一致。

**数值基准验证与结果分析**：在颗粒沉降验证（第3.1节）中，通过模拟粘性流体中圆形颗粒的沉降过程，成功捕捉了典型的冯·卡门涡街尾流结构，颗粒垂向速度和位置随时间的变化与基准结果吻合良好，验证了H-M耦合框架捕捉FDEM固体颗粒与LBM流体格点间相互作用的能力。在运动颗粒溶质扩散验证（第3.2节）中，模拟了平移圆形颗粒周围的溶质扩散过程，捕捉到颗粒周围明显的高浓度区域及显著的上游扩展，浓度场分布及沿通道中心线的浓度剖面与基准数据吻合，验证了H-C耦合框架及几何边界检测算法的有效性。在多晶型颗粒溶解验证（第3.3节）中，对比了圆形、方形和六边形矿物颗粒的溶解过程，结果表明圆形颗粒经历均匀各向同性溶解，保持原有形态；而角状颗粒表现出明显更快的初始边界后退，高曲率尖角溶解速度远快于平坦边缘，呈现典型的角圆化现象（corner-rounding），后期逐渐趋近于平滑圆形的动力学行为，与实验观察一致，验证了拓扑更新算法捕捉复杂几何与化学耦合的能力。在反应性流动及动态界面演化验证（第3.4节）中，模拟了流道内圆形颗粒的反应性溶解过程，连续流体注入产生了非对称浓度分布，上游正面溶解最快，下游尾迹区溶解速率最低，颗粒逐渐由圆形演化为流线型几何形态，颗粒面积、周长及平均溶解速率的时序演化与基准结果高度一致，验证了LBM-FDEM框架准确再现对流-反应流驱动的固液界面复杂动态演化的能力。

**裂缝溶解数值模拟与机制分析**：在裂缝溶解模式复现（第4.1节）中，通过调节边界压力和初始孔径匹配无量纲Péclet数（Pe）和几何比（λ），成功复现了微流控实验中两种截然不同的溶解模式：Regime I（Pe=34.1，λ=0.27）表现为高度局域化的入口溶解，Regime II（Pe=341，λ=1.02）表现为均匀的通道扩宽，模拟结果与实验在最终溶解形态和归一化反应表面积（A/A0）演化趋势上均高度吻合。模拟提供的高分辨率浓度场揭示了局域化溶解源于低注入速度下反应物在入口缓冲区内被完全消耗，而均匀扩宽则归因于高流速下反应物被强对流输运至裂缝深处。在不同压力梯度下溶质运移与通道溶解时空演化研究（第4.2节）中，通过固定几何参数（λ=0.44）和初始孔径（d2=0.132mm），系统增大注入压力以生成Pe=3.41、34.1和102.3三种工况。结果表明：在低Pe条件下，溶质运移受扩散控制，反应物在接触上游矿物表面时即被完全消耗，呈现典型的表面溶解模式（face dissolution）；在中等Pe条件下，对流与扩散高度竞争，反应物进一步渗入裂缝但沿程逐渐消耗，导致溶解速率出现显著空间梯度，最终将均匀裂缝塑造成锥形形态（tapered shape）；在高Pe条件下，强对流确保未反应的反应物沿整个流动路径持续供给，同时生成的产物被快速冲刷出下游边界，系统呈现明显的沟槽化扩宽模式（channelized widening pattern），驱动裂缝壁面快速均匀地向侧向扩展。这些不同的溶解和溶质迁移特征本质上是由于控制传质机制从扩散主导向对流主导的根本转变。

**讨论与结论**：研究结果表明，所提出的LBM-FDEM耦合框架不仅能忠实再现实验中观察到的形态转变，还能通过提供空间浓度场明确揭示其驱动机制，即不同Pe下的不同溶解过程从根本上取决于对流通量与局部界面反应速率之间的竞争。该方法为量化矿物溶解和溶质运移的H-M-C过程提供了一个稳健且通用的数值工具。当前的几何映射和拓扑更新算法因计算复杂度限制仅限于二维，未来工作将把LBM-FDEM框架扩展至三维。该研究为地质CO2封存安全性评估、裂缝性储层反应性运移模拟以及地下工程中复杂多场耦合问题提供了一种高保真度的数值模拟手段。