采用混合DEM-FDM方法对热机械耦合效应对井筒稳定性的微观分析

《Computers and Geotechnics》：Microscopic analysis of thermo-mechanical coupled effects on wellbore stability using a hybrid DEM-FDM approach

　　德浩·孟|费·王|海因茨·科尼茨基|何鹏国|陈林 中国沈阳110819，东北大学智能深部采矿与装备国家重点实验室 摘要 高效且经济高效的地下热能储存系统常常面临由热流体注入引起的井筒不稳定和破坏问题。离散元方法非常适合模拟高温条件下的井筒开裂与稳定性问题。然而，将离散元方法应用于深部地下系统时，会受到模型尺度要求以及大量元素带来的计算成本限制。为克服这些难题，本研究提出了一种离散-连续耦合方法，将离散元方法与有限差分框架相结合。通过建立新的传热方程和嵌入算法，实现了离散颗粒与周围连续体之间的热交换，并通过解析解对该方法进行了验证。在此基础上，利用三轴压缩试验数据构建了工程尺度的井筒模型。该耦合模型成功再现了热流体注入过程中的井筒传热与热裂解行为。进一步分析了三种不同工况下的井筒稳定性：（1）全水泥套管支撑；（2）无套管支撑；（3）流体柱压力作用。研究结果表明，实验室尺度与实际工程尺度之间的损伤演化过程既有共性，也存在尺度依赖性。通过将微观热力学过程与宏观岩体响应相联系，本研究深化了对晶质岩石中温度诱发损伤机制的理解，为评估和管理深部地热及能源储存系统的稳定性提供了有力工具。引言 地下结构稳定性及长期可持续性所受到的热效应影响日益受到重视，因为温度变化在能源生产和储存相关的地下作业中普遍存在（Vilarrasa和Rutqvist，2017；Im等人，2021）。尤其是地下热能储存系统，作为一种通过将多余热量或可再生能源储存在地下以供日后提取的方式，有望解决季节性能源供需不平衡的问题（Lee，2013）。地下热能储存系统的运行通常涉及在地质储层中循环注入、储存和抽取热流体，这些储层包括高孔隙度的沉积岩层以及裂隙发育的晶质岩石。这类系统的长期效率与安全性在很大程度上取决于井筒及其周围岩体在运行过程中所承受的耦合热力学载荷的机械完整性（Mutlu等人，2023）。在地下热能储存和地热应用中，井筒不稳定问题一直较为突出，尤其是在注入段附近，那里剧烈的温度梯度与较高的原位应力相互作用（Schmitt等人，2012）。井筒壁的快速加热或冷却可能会产生巨大的热应力，这些应力可能超过岩石的抗拉强度或剪切强度，从而导致岩石开裂、剥落或破碎（Yan等人，2014），进而增加非生产时间（Zhang等人，2009），并造成巨大的经济损失（Zeynali，2012）。因此，了解和预测井筒周围由热力学作用引发的损伤问题在地下工程中越来越重要，这对地下能源系统的设计、运行和风险管理具有直接意义。花岗岩作为地热和深部能源储存项目中的常见围岩，其热力学行为已得到了广泛研究。实验研究表明，热处理会改变花岗岩的关键力学性能，包括泊松比（Kumari等人，2017a）、断裂角度（Chen等人，2018）以及破坏应变（Singh等人，2015；Wang等人，2019a；Wang和Konietzky，2019；Meng等人，2021）。这些宏观响应源于微观结构的变化，如矿物膨胀不匹配、微裂纹萌生以及裂纹聚合（Srinivasan等人，2022）。然而在深部地质环境中，升高的温度与原位应力会对花岗岩的力学性能产生相互影响的作用（Gao等人，2022；Kumari等人，2017b；Ma等人，2023；Wang等人，2015；Wong，1982；Yi等人，2024；Yin等人，2021）。尤其是温度与围压之间通过控制微裂纹的形成、张开和闭合而以复杂方式相互作用（Batzle等人，1980；Yong和Wang，1980；Heard和Page，1982；Griffiths等人，2018）。虽然较高的围压可以抑制热裂解现象，但高温仍可能导致岩石强度显著下降。例如，即使在400 MPa的围压作用下，将温度升高到680℃时，花岗岩的差异强度也会从1660 MPa降至1030 MPa（Wong，1982）。更多近期研究指出，在高温（800℃）与高围压（150 MPa）共同作用下，岩石会出现更为分布广泛的整体损伤，因为岩石开始向延性转变（Meyer等人，2024）。同样，Wang等人（2024a）也报告了在高温高压条件下岩石出现的空间上不均匀的热裂解模式。由于花岗岩属于多晶且性质不均一，若没有将矿物尺度上的损伤过程与宏观岩体响应联系起来的理论框架，实验室观测结果就无法直接推广到实际工程规模的地质条件下。数值建模为连接实验室尺度观测结果与工程尺度下的岩体响应提供了一种有效途径。包括有限元方法（Guo等人，2024）、有限差分方法（Wang等人，2024a）以及相场方法（Yue等人，2024）在内的连续介质方法已被广泛用于岩石和地质材料的热力学分析。这些方法在处理大规模边界值问题时效率较高，但对于多晶质岩石中微裂纹的萌生与扩展过程则难以精确描述。通常情况下，这些方法是通过预定义的弱化本构律或沿预定义平面进行网格分离来描述损伤与破坏过程的（Enayatpour等人，2018；Cao等人，2025；Liu和Wang，2025；Liu和Wong，2025；Liu等人，2026；Wang等人，2026a；Wang等人，2026b）。相比之下，离散元方法则将岩石视为由相互作用的单元组成的集合体。花岗岩的宏观行为由其矿物晶体的集体相互作用决定，而离散元方法能够很好地模拟这种相互作用（Cundall等人，1996）。在离散元方法中，各个单元是在其界面处相互作用的，而非被视为连续体，这使得该方法能够模拟岩石中不同矿物的行为及其相互作用。因此，离散元方法是探究温度变化导致力学性能变化的微观机制的重要工具（Wang等人，2024b；Zhao，2016）。例如，一项离散元研究显示，除了人们熟知的温度弱化效应外，冷却过程中矿物碎片之间的新式互锁作用反而可能增强试样的强度（Wang等人，2024b）。离散元方法的另一个优点在于它能够直观地表示裂纹情况（Cai等人，2007）。当单元之间发生破坏时，其界面会失去抗拉强度和黏聚剪强度，从而允许发生滑动。这使得离散元方法能够有效模拟裂纹的萌生与扩展过程（Peng等人，2017；Tomac和Gutierrez，2020；De Silva等人，2023a；De Silva等人，2023b）。因此，离散元方法能够再现不同力学与热学条件下的井筒不稳定模式（Lin等人，2020；Xiang等人，2023a；Xiang等人，2023b）。需要指出的是，离散元方法既包括基于颗粒的框架，也包括基于块的框架。离散元方法最初是基于颗粒的框架发展起来的，其宏观行为是由刚性颗粒之间的接触相互作用所决定的。而在基于块的离散元方法中，材料则由相互作用的块体来表示，这种方法能更好地描述复杂的颗粒几何形态（Wang等人，2026a）。一些基于块的方法还将单个块体进一步离散为内部网格，以便在每个块体内实现弹性变形，由此形成了组合型有限离散元方法（FDEM）（Aboayanah等人，2022；Sun等人，2022；Yan等人，2022；Liu和Wong，2025；Liu和Wang，2025）。尽管这些与离散元方法相关的框架在单元表示方式上存在差异，但在应用于工程尺度的井筒问题时，它们都面临着类似的挑战，那就是计算成本问题（Vaezi等人，2025）。具体而言，模型尺寸与井筒直径的比值往往较低，这会导致边界效应不真实，且无法准确反映应力重分布情况，从而影响模拟的精度。这一限制源于扩大模型规模所带来的计算成本。随着模型尺寸的增大，需要计算的单元数量及相互作用数量也会迅速增加，进而导致计算效率大幅下降（Breugnot等人，2016）。虽然从理论上来说，大规模模型能够捕捉到井筒周围复杂的应力扰动与变形机制，但由于上述计算限制，其实际应用受到了阻碍。模型规模与计算可行性之间的这种权衡，凸显出有必要采用更先进的建模方法或优化技术来克服这些固有缺陷。为了解决工程尺度问题模拟中的计算难题，人们常常提出离散-连续耦合算法（Guo和Zhao，2014；Potyondy和Cundall，2004）。这类方法旨在在降低计算成本的同时，保留离散元方法在模拟微观行为方面的优势（Bai等人，2022；Zhang等人，2019）。典型的耦合算法会将建模域分解为多个子域，对特殊区域应用离散元方法，而对其他区域则采用连续介质方法，以此降低计算负荷（Cai等人，2007）。在那些需要提前识别潜在破坏区域的实际应用中，这种方法尤为有用（Fish等人，2021）。离散-连续耦合方法已被成功用于模拟纯力学耦合条件下的隧道诱导剥落现象（Potyondy和Mas Ivars，2020）以及井筒破裂现象（Wang等人，2020）。为了将这一方法应用于更复杂的多物理场条件，人们还开发出了动态耦合（Li等人，2015）、水力耦合（Zhang等人，2017；Zhang等人，2019）以及斜向入射（Gu等人，2024）等算法。这些方法依赖于有限差分代码FLAC（Itasca 2020）与颗粒流代码PFC（Itasca 2021）之间的TCP/IP套接字连接，从而实现将FLAC中连续域的节点位移传递给离散元方法中的离散颗粒。然而，尽管地下工程中对热力学过程的研究兴趣日益浓厚，但仅有少数研究探讨了离散元框架与连续介质框架之间的热耦合问题（Liu和Wong，2025；Liu和Wang，2025），而且由于这两种数值方法之间缺乏热交换功能，目前仍然没有一种完全整合的离散元-有限差分法热力学框架。在本研究中，我们针对这一不足，开发了一种基于颗粒的离散元-有限差分法框架，用于进行耦合热力学模拟。第2节中，我们提出了一个传热公式与嵌入算法，使得离散颗粒与周围连续体之间能够实现显式的热交换，从而可以将微观尺度的热裂解过程纳入工程尺度的井筒模型之中。第3节中，我们通过解析解对耦合方案进行了验证，并利用三轴压缩试验数据对模型进行了标定。第4节中，我们将经过验证的模型应用于研究在不同支撑条件下的井筒传热、热裂解及稳定性问题，这些条件代表了实际的工程场景。通过将微观热力学机制与宏观岩体响应相联系，本研究为理解晶质岩石在热作用下的损伤演化规律提供了新见解，同时也为评估地热能储存系统、地热工程以及其他深部地下工程中的地质稳定性及其他与热相关问题提供了可靠的建模工具。章节节选 结合颗粒模型 本研究中的离散元模拟采用了三维颗粒流代码PFC3D（Itasca 2021）。该代码将宏观行为视为由各个颗粒的集体相互作用所产生，而这些相互作用则是通过各种接触模型来定义的。结合颗粒模型（Potyondy和Cundall，2004）常被用于模拟岩石这类完整脆性材料的变形与损伤行为。本研究选择了线性平行键合模型，因为它能够根据实际情况调整键合强度。代表性基本体积设置与热参数确定 本研究选择来自中国某正在开发中的地热储层的公和花岗岩作为模型参数标定的目标岩石。此前已有研究探讨过公和花岗岩的物理性质（Cheng等人，2021）以及力学性质（Wang等人，2024c），这是一种典型的多晶质岩石，由35%的石英、30%的斜长石、23%的钾长石以及12%的其他矿物组成。公和花岗岩的各项性质详见表1。边长为……的立方形BPM模型 考虑热力学耦合因素的井筒稳定性 在向地下能源储存系统注入热流体时，井筒底部周围的岩石会经历显著的温度变化（见图6a）。热流体注入储层会在井筒底部附近引发热冲击。这就导致了井筒附近区域内的原位应力集中与热应力之间产生复杂的相互作用，进而增加了岩石受损的风险。此外，不同材料的线膨胀系数也存在差异……实验室尺度与井筒尺度响应的比较 通过实验室尺度下对花岗岩进行的受围压高温压缩试验结果，来模拟井筒尺度下的行为。如图14所示，花岗岩在较低围压条件下的热弱化现象更为明显。例如，在高围压环境下，处于220℃的样品产生的热裂纹明显较少（见图14a）。应力-应变曲线也证实了这一趋势：在45 MPa的围压下，200℃时的强度降低了21.4 MPa，而在没有围压的情况下则下降了51.2 MPa。

结论
本文提出了一种用于将基于BPM的离散元方法与基于连续体的有限元方法进行TM耦合的热交换算法，并将其应用于地热井筒模型中，重点研究地热能储存过程中热效应对井筒稳定性的影响。同时研究了花岗岩在不同尺度下对热效应的反应。主要结论如下：
(1)所提出的TM耦合算法能够有效整合离散元方法和有限元方法。这种混合模型减少了边界效应并降低了计算量。

作者贡献声明
孟德豪：撰写——初稿、可视化、验证、软件、方法论、研究、正式分析、数据整理。王飞：撰写——初稿、验证、监督、软件、资源、项目管理、概念设计。海因茨·科尼茨基：撰写——审阅与编辑、软件。何本国：撰写——审阅与编辑、资金获取。陈琳：撰写——审阅与编辑。

利益冲突声明
作者声明没有已知的可能影响本文所述工作的财务利益或个人关系。

致谢
本研究部分得到了中国国家自然科学基金的支持（项目编号：W2431040、52222810、52525808）。