一种用于模拟含裂缝多孔储层的全耦合热-水-力学双孔隙度方法

《Computers and Geotechnics》：A fully coupled thermo-hydro-mechanical dual-porosity method for modeling fractured porous reservoirs

　　

**摘要**

fractured多孔介质的热-水-力学（THM）分析仍然是精确油藏模拟中的关键挑战，尤其是针对非均质域。本研究开发了一套完全耦合的计算框架，以增强裂缝系统的双重孔隙度建模，考虑固体变形、流体流动和热量传递之间的相互作用。推导了质量和热量传递形状因子的通用形式，并将渗透基质块中的导热和对流机制均纳入其中。与采用恒定传递参数的常规模型不同，所提出的方法通过对具有相裂缝构型和物理特性的代表性单元进行直接分析，分别评估质量和热量传递的瞬态形状因子。这些形状因子能够捕捉不断演化的传递过程，并在每个时间步动态地集成到双重孔隙度框架中。该方法在复杂问题上展示了其准确性和高效性，包括具有多个区域和不同块体几何形状的非均质地热油藏。结果表明，瞬态质量和热量形状因子表现出独特的、随时间变化的行为，而静态假设会在瞬态条件下导致巨大的建模误差。研究还探讨了在预测地热油藏行为时考虑实际裂缝构型的重要性。结果表明，所提出的基于有限元法的THM双重孔隙度方法是一种有前景的工具，可用于地热油藏热量开采性能的实际估算，特别是在地质非均质性影响流动和热力状态的区域。

**1. 引言**

有各种工程问题涉及含裂缝多孔域的热-水-力学建模，例如油气开采、冷热水注入以及热油藏优化（Sun et al., 2021, Ma et al., 2022, Mortazavi et al., 2022）。这些问题出现在地热能量开采、天然裂缝性油藏的烃类回收、地质碳封存、发热废物的深地质处置以及裂隙岩体中的热能储存等应用中（Koh et al., 2011, Saceanu et al., 2024, Zimmerman and Paluszny, 2024, Taghvaei et al., 2026）。在这些应用中，裂缝充当流体流动和热量传输的优势通道，而周围的基质块则作为流体和热能储存区域。因此，基质-裂缝传递项的准确性直接影响预测的压力演化、温度分布、热突破行为、产出液温度和整体热提取性能，使得可靠估算传递系数成为裂隙多孔介质工程尺度模拟的关键。在含有大量互连裂缝的域中，裂缝的存在显著改变了固体变形、流体流动和热量传递过程，从而增加了数值模拟的复杂性（Koh et al., 2011, Mortazavi et al., 2023）。目前提出了多种方法来建模非均质性，如裂缝，包括离散裂缝法（Fu et al., 2016, Gallyamov et al., 2018）、相场技术（Khoei et al., 2025, Khoei and Tabatabaei, 2026, Miehe et al., 2010, Yin et al., 2022）、扩展有限元法（Mohammadnejad and Khoei, 2013, Khoei et al., 2018, Wang et al., 2023, Mortazavi and Khoei, 2024）、近场动力学计算方法（Ni et al., 2022, Ni et al., 2023）、等几何分析框架（Saeedmonir and Sadeghirad, 2026）、等效连续介质法（Lu et al., 2013, Khoei et al., 2016）、计算多尺度技术（Nguyen et al., 2012, Khoei and Hajiabadi, 2018, Khoei et al., 2023, Saeedmonir et al., 2024）以及双重孔隙度法（Wang et al., 2018a, Liu et al., 2022, Hu et al., 2025）。在直接方法中，如上文中提到的前三种方法，裂缝被显式建模，从而产生更准确的结果（Zeng et al., 2020, Cusini et al., 2021, Liu et al., 2024）。然而，当模拟具有大量裂缝的大尺度域时，这些方法计算量巨大。相比之下，等效模型不显式建模裂缝，因此计算效率高。然而，它们通常会导致近似结果，不适合精确模拟。在多尺度方法中，大尺度行为通过研究微观层次提取（Khoei and Saeedmonir, 2021, Saeedmonir and Khoei, 2022, Svenning et al., 2016）。这些方法能够以适当的精度模拟非均质性；然而，当存在多个裂缝和非均质性时，它们的计算成本很高。

建模裂隙域的另一种方法是基于双重孔隙度法。该方法最初由Barenblatt等人（1960）提出，其中裂隙介质被分解为两个重叠的子域，即裂缝连续介质和基质连续介质。在此方法中，单个裂缝不被显式离散化；相反，其效应通过等效裂缝连续介质来表示。这使得双重孔隙度法成为一种成本效益较高的框架，因此被广泛用于模拟大规模裂隙域，如地热油藏（Liu et al., 2019, Hajiabadi and Khoei, 2019）。在双重孔隙度概念中，基质和裂缝连续介质通过表示两个域之间质量和/或热量传递的交换项进行耦合。这些项通常使用称为"形状因子"的传递参数来量化。双重孔隙度模型的准确性在很大程度上取决于该参数的正确评估。因此，过去几十年间，文献中开发了各种形状因子估计模型。Warren和Root（1963）使用双重孔隙度法对不可变形多孔介质中的单相流体流动进行建模。他们根据基质块的尺寸提出了一个形状因子的关系式。Kazemi等人（1976）采用有限元法提出了一个水力形状因子的几何关系式。尽管揭示出该模型会产生低估结果，但它仍被文献中广泛使用（Jerbi et al., 2017, Cordero et al., 2019）。Zimmerman等人（1993）进行了解析研究以求解一个非线性微分方程。他们的结果表明，基质-裂缝传递随时间演化，不应始终由伪稳态交换系数来表示。这项研究随后推动了瞬态和边界条件依赖形状因子公式的进一步发展。在他们的研究中，开发了一个水力形状因子的关系式，该关系式后来由Lim和Aziz（1995）在时域中验证，并由Mathias和 Zimmerman（2003）在Laplace域中验证。Kazemi和Gilman（1993）考虑了阶梯边界条件，并提出了水力形状因子的不同关系式。他们的模型与Zimmerman等人（1993）针对恒定裂缝压力条件所提出的模型相同。Hassanzadeh和Darvish（2006）首次考虑了裂缝压力制度连续变化的情况。他们的研究结果表明，水力形状因子具有随时间变化的行为，在整个传递过程中使用恒定值会导致瞬态阶段的显著误差。他们还证明，在特定块体中形状因子的值在不同边界条件下可能会发生变化。Rostami等人（2020）采用有限元法使用不规则块体的插值特征长度对水力形状因子进行数值计算。在以上所有研究中，固体变形的效应被忽略。Wang等人（2018b）考虑了基质系统的应力敏感性，以提出单相流体流动中形状因子的公式。Khoei和Taghvaei（2024）将可变形介质中的等温流体流动应用于推导水-力学形状因子的关系式，其中考虑了固体变形的影响。

双重孔隙度概念最适合对流场景，即裂缝区域的流速远大于基质域中的流速（Ding, 2022）。然而，在热量传递问题中，导热机制也起着重要作用。因此，使用双重孔隙度法模拟热问题通常伴随误差，特别是在瞬态阶段。另一方面，在地热模拟中，温度扩散的瞬态持续时间远长于压力扩散（van Heel et al., 2008）。在文献中，水力形状因子被广泛用于计算热量传递形状因子（Heinze and Hamidi, 2017）。然而，Abdallah等人（1995）指出，需要一种关系来考虑导热（扩散）和对流（平流）两种机制的热形状因子，这具有重大意义。Shaik等人（2011）使用双重孔隙度法数值研究了地热油藏的热量提取。他们的结果表明，常用的形状因子会导致温度被低估。van Heel等人（2008）针对一维热问题开发了瞬态热量传递形状因子，该值在恒定裂缝温度条件下的稳态时与Lim和Aziz（1995）所提出的值一致。然而，他们的结果表明，该值而非瞬态形状因子会在瞬态阶段导致误差。Heinze等人（2017）在一维稳态假设下推导了一维问题的热形状因子解析公式。他们指出，瞬态一维热形状因子取决于裂缝开口度、流速和热物性。Abbasi等人（2017, Abbasi et al., 2018）为两个不透水基质块之间的单条裂缝提出了热形状因子的公式。他们表明，基质块的尺寸影响一维热形状因子的值。Taghvaei和Khoei（2025）研究了裂缝与不透水基质之间的热量传递，针对方形和非方形块体进行了分析。根据他们的结果，在瞬态阶段使用恒定值的热形状因子可能导致显著误差。此外，热形状因子的量级受块体几何形状、边界条件和域的热物性影响。

尽管双重孔隙度法已应用于油气和热油藏的建模，但该方法尚未广泛用于耦合热-水-力学问题。Bai等人（1995）在热-水-力学问题中使用了双重孔隙度法，考虑了恒定裂缝压力和恒定通量边界条件。在他们的研究中，动量平衡方程和连续性方程进行了耦合，但能量平衡方程未耦合。Masters等人（2000）将热过程整合到Ghafouri和Lewis（1996）开发的水-力学模型中；然而，他们忽略了温度对固相的影响。Khalili和Selvadurai（2003）提出了一个耦合热-水-力学模型，用于模拟完全饱和的弹性域。在他们的研究中，质量和热量传递形状因子被假定为相等且恒定，其值基于Warren和Root（1963）模型确定。Nair等人（2004）也使用有限元法进行了热-水-力学研究；然而，连续性方程和能量平衡方程未耦合。Hosking等人（2020）提出了一个耦合双重孔隙度热-水-力学模型，用于模拟煤层CO2封存过程中的非等温气体流动。他们考虑了固体参数对连续性方程的影响；然而，他们忽略了动量平衡方程和连续性方程之间的耦合。

根据以上讨论，尽管已有若干研究将双重孔隙度方法扩展至非等温或热-水-力学（THM）问题，但一个完全耦合的THM双重孔隙度框架——其中动量平衡、流体连续性和能量守恒方程被一致地集成——尚未得到充分发展。此外，在大多数现有公式中，基质-裂缝传递项使用恒定形状因子评估，质量传递和热量传递的形状因子要么被假定为相同，要么未被明确区分。为了解决这些局限性，本研究提出了一个用于裂隙多孔介质完全耦合THM分析的计算双重孔隙度框架。其主要目标是提高瞬态基质-裂缝质量与热量交换的预测精度，同时保持工程尺度裂隙油藏模拟所需的计算效率。在该框架中，动量平衡、流体连续性和能量守恒方程被完全耦合，并针对基质连续介质与裂缝连续介质之间的质量传递和热量传递分别引入了不同的时变形状因子。所提出的方法首先基于裂缝几何和材料特性在一般域中识别出不同的区域。然后，从每个区域中指定一个具有相同几何和物理特性的适当单元。该单元在一致的边界条件下进行分析。在每个时间步，使用所提出的瞬态关系评估质量和热量传递形状因子，其中考虑了热-水-力学过程的耦合。

随后，所计算得到的形状因子被应用于主域的双重孔隙度分析中。通过引入针对当前问题的几何条件和物理条件量身定制的瞬态且耦合的传热与传质形状因子，该方法能够有效提高双重孔隙度方法在裂隙域THM建模中的精度。此外，该方法能够处理包含具有任意裂隙构型的不同区域的非均质域。本文的组织结构如下：第2节介绍了完全耦合热-水-力分析的支配方程，并推导了相应的瞬态传热和传质形状因子。第3节介绍了双重孔隙度方法的有限元离散化。第4节描述了所提出的耦合THM双重孔隙度建模的计算框架。第5节研究了瞬态传质和传热形状因子在方形和非方形块体中的性能和精度。此外，为了研究所提出的THM算法的性能和效率，分析了三个数值算例。最后，第6节给出了一些结论性评述。

**2. 物理模型**

在双重孔隙度模型中，裂隙域被分解为两个重叠的子域：裂隙域和基质域，如图1所示。基本上，裂隙子域的渗透率远高于基质子域，因此输运机制在该域中起着关键作用。另一方面，基质域具有较高的储集能力。这两个子域通过适当的交换项进行耦合。相应地，物理模型基于以下关键假设：

- 假设基质域和裂隙域均为完全饱和状态
- 基质是可渗透的，允许在裂隙域和基质域之间发生质量交换
- 假设在多孔基质中，孔隙水和岩石瞬间达到局部热平衡（LTE）
- 由于基质和裂隙子域内存在流体流动，对流（平流）和传导（扩散）两种机制均参与传热过程
- 基质块是可变形的，固相与热之间的耦合是单向的；即热应变影响固体变形，但固体变形不影响传热过程


**图1. 双重孔隙度框架的概念示意——用于裂隙多孔介质完全耦合热-水-力（THM）分析**

需要注意的是，超出本假设范围的过程，例如局部热非平衡、非饱和或多相流条件、裂隙扩展以及损伤导致的裂隙开度演化等，均不在本研究考虑范围内。图2展示了本文采用的完全耦合热-水-力过程示意图。基质连续介质被建模为完全耦合的THM介质，其中力学、水力与热场通过有效应力、体积应变率、热应变、流体储集以及对流传热机制相互关联。相比之下，裂隙连续介质由耦合的热场和水力场表示，而力学变形则假定主要发生在基质中。两个连续介质通过基质-裂隙质量交换和热交换项进行联系。


**图2. 完全耦合热-水-力（THM）过程示意图**

**2.1. 支配方程**

为了模拟饱和可变形介质中的非等温流体流动，耦合方程的支配系统包括力学平衡方程、流体连续性方程和能量平衡方程。在u-p-T公式体系中，温度T、流体压力p和位移场u被作为基本变量。基于双重孔隙度概念，裂隙介质中的问题被分解为两个相互作用的连续介质，分别代表基质域和裂隙域。这些重叠的连续介质通过质量交换项和热交换项进行耦合。在该框架中，基质与裂隙之间的压力差和/或温度差驱动了流体和热量在两个域之间的传递。此外，在每个连续介质内部，热量传递通过对流和传导两种机制进行。孔隙流体与固体基质之间也采用了局部热平衡（LTE）假设。

在以下公式中，标量用斜体符号表示，向量和矩阵量用粗体符号表示。每个向量和矩阵的物理意义和维度将在其首次出现时予以说明。下标"m"和"f"分别表示基质连续介质和裂隙连续介质。

遵循饱和可变形多孔介质THM的标准公式体系（Lewis和Schrefler, 1998; Diersch, 2014; Khoei, 2015; Chan等, 2022），基质连续介质的平衡方程强形式可写为：

(1) ?·σ ? ρü + ρb = 0

其中σ为总应力。在本研究中，拉应力和拉应变定义为正值，而孔隙压力在受压时为正值。相应地，Biot有效应力关系式表达为σ = σ′ ? αpm_m，其中σ′代表有效应力，m为单位向量，m = [1 1 0]?。在存在热应变时，σ′可表达为：

σ′ = D:(ε ? ε^T)

其中D为弹性张量。在上式中，ρ、b和pm分别表示混合物的平均密度、体力向量和基质压力。此外，α为Biot系数，其计算式为α = 1 ? KT/KS，其中KT和KS分别为土样和固体颗粒的体积模量。本质边界条件和自然边界条件分别指定为：在Γu上u = ū，在Γt上σ·nΓ = t?，其中nΓ表示外部边界的单位外法向量。

遵循饱和可变形多孔介质的双重孔隙度公式体系（Rahman和Lewis, 1999; Hajiabadi和Khoei, 2019），双重孔隙度方法中的连续性方程可表达为：

(2) ?·v_w,m + α?·u? + Cm?m?m ? βw?m + q = 0（针对基质域）

(3) ?·v_w,f + Cf?f?f ? βw?f ? q = 0（针对裂隙域）

在上述方程中，k为渗透率，C表示总压缩性，?表示域的孔隙率，βw对应体积热膨胀系数。流体速度v_w根据广义达西定律（Lewis和Schrefler, 1998; Diersch, 2014; Khoei, 2015; Chan等, 2022）获得：

(4a) v_w,m = k_m(??pm + ρb ? ρü)

(4b) v_w,f = k_f(??pf + ρb)

此外，q为表示基质域和裂隙域之间交换的流体通量的耦合项。遵循Warren和Root（1963）的研究，该微分方程组可简化，流体交换项可评估为：

(5) q = s_M·k_m(p?m ? p?f)

其中s_M为传质形状因子。另外，p?m和p?f分别为平均基质压力和平均裂隙压力，其表达式为：

(6a) p?m = (1/Vm)∫_{V_m} p_m dV

(6b) p?f = (1/Vf)∫_{V_f} p_f dV

方程(2)和(3)可通过考虑本质边界条件和自然边界条件来求解，即在Γp上pf = p?，在Γ_qw上v_w·nΓ = q?。

遵循Heinze和Hamidi（2017）的热的双重孔隙度公式体系，并应用傅里叶热传导定律，基质域和裂隙域的温度由能量平衡方程控制：

(7) ρcm?m + (ρc)_w·v_w,m·?Tm ? λm?2Tm ? Q ? Q_HW,m = 0（针对基质域）

(8) ρcf?f + (ρc)_w·v_w,f·?Tf ? λf?2Tf + Q ? Q_HW,f = 0（针对裂隙域）

在上述方程中，c表示比热容，λ为导热系数，Tm和Tf分别为基质温度和裂隙温度。此外，ρcm、ρcf和(ρc)_w分别表示基质连续介质、裂隙连续介质和孔隙水的体积热容。另外，Q为基质域和裂隙域之间交换的热通量，其计算式为：

(9) Q = s_H·λ_m(T?m ? T?f)

其中s_H为传热形状因子，平均基质温度和平均裂隙温度T?m和T?f定义为：

(10a) T?m = (1/Vm)∫_{V_m} T_m dV

(10b) T?f = (1/Vf)∫_{V_f} T_f dV

通过质量传递在基质域和裂隙域之间传输的热量由Q_HW,m和Q_HW,f项表示，表达为：

(11a) Q_HW,m = ?·cw·(T_m ? T_f)·H(?)

(11b) Q_HW,f = ?·cw·(T_m ? T_f)·H(??)

其中H为海维赛德阶跃函数，对于正值参数其值为1，否则为零。在上述表达式中，?表示从裂隙域到基质域的单位时间内单位体积的传质速率。需要注意的是，每个域均为完全饱和状态，温度变化仅由裂隙连续介质和基质连续介质之间的流体交换引起。此外，?可计算为：

(12) ? = ρw·s_M·k_m(p?m ? p?f)

需要注意的是，Q_HW,m和Q_HW,f项代表与基质-裂隙质量交换相关的热量贡献。遵循Heinze和Hamidi（2017）的研究，采用海维赛德函数根据质量传递方向将此对流项分配给接收连续介质。相应地，这一贡献并不代表基质连续介质与裂隙连续介质之间的总热交换，因为传导热传递通过方程(9)中的热交换项单独考虑。

还需进一步强调，方程(5)和(9)中的交换项遵循经典双重孔隙度公式体系，其中基质-裂隙的传质和传热以两个连续介质之间的平均压力和温度差表示。在常规方法中，相关的形状因子通常被视为恒定的几何参数。相比之下，本研究在下一节中将推广这些量，以捕捉基质块的瞬态THM响应。

最后，为求解方程(7)和(8)，施加本质边界条件和自然边界条件，即在ΓT上Tf = T?，在Γ_qT上q_T = q?_T。

**2.2. 传质和传热形状因子**

传质和传热函数在双重孔隙度公式体系中起着核心作用，因此准确估计这些参数至关重要。在完全耦合热-水-力系统中，多种相互作用的工艺影响传递形状因子。所需瞬态形状因子的推导遵循经典双重孔隙度模型的一般框架，其中基质-裂隙交换项表示为基质块与周围裂隙连续介质之间平均压力或温度差的函数。然而，在本公式体系中，传递系数不被假定为恒定的几何参数。相反，它们通过对代表性质基质块的瞬态响应，将支配方程在基质体积上进行积分，并将所得边界通量与相应的双重孔隙度交换表达式关联来获得。因此，所推导的传质和传热形状因子能够捕捉基质块的瞬态状态。

为在THM框架中开发传质形状因子的广义表达式，考虑一个代表性基质块。通过将基质连续性方程(2)与广义达西定律(4)结合，块体内的压力场由以下方程控制：

(13) k_m?2pm = α?·u? + Cm?m?m ? βw?m

对基质块边界对方程(13)进行积分，应用散度定理，并将所得关系式除以Vm，得到以下关系：

(14) (1/Vm)∫_{?Vm} k_m(?pm/?n)dA = (1/Vm)∫_{V_m} α?·u? dV + Cm?m(?p?m/?t) ? βw(?T?m/?t)

其中?Vm为基质块的边界。

此外，流出通量q可沿基质块边界积分为：

(15) q = ?(1/Vm)∫_{?Vm} k_m(?pm/?n)dA

由方程(14)和(15)，可得到以下关系：

(16) ?q = (1/Vm)∫_{V_m} α?·u? dV + Cm?m(?p?m/?t) ? βw(?T?m/?t)

将上述关系与Warren和Root（1963）提出的方程(5)进行比较，可得到传质形状因子的广义形式：

(17) s_M = [1/(k_m(p?f ? p?m))]·[(1/Vm)∫_{V_m} α?·u? dV + ?mCm(?p?m/?t) ? βw(?T?m/?t)]

其中第一项和第三项分别考虑了固相效应和热变化效应。

为了推导广义传热形状因子，考虑基质中的传导和对流两种机制，采用类似的步骤。应用基质能量平衡方程(7)，代表性基质块内的温度场由以下方程控制：

(18) ρcm?m + (ρc)_w·v_w,m·?Tm = λm?2Tm

对方程(18)进行积分，将所得关系除以Vm，并应用散度定理，得到以下关系：

(19) ρcm(?T?m/?t) + (1/Vm)∫_{V_m} (ρc)_w·v_w,m·?Tm dV = (1/Vm)∫_{?Vm} λm(?Tm/?n)dA

流出通量Q定义为：

(20) Q = ?(1/Vm)∫_{?Vm} λm(?Tm/?n)dA

由此可得：

(21) ?Q = ρcm(?T?m/?t) + (1/Vm)∫_{V_m} (ρc)_w·v_w,m·?Tm dV

对比方程(21)和(9)，可得：

(22) s_H = [1/(λm(T?f ? T?m))]·[ρcm(?T?m/?t) + (1/Vm)∫_{V_m} (ρc)_w·v_w,m·?Tm dV]

上述方程提供了可渗透基质中传热形状因子的广义评估关系。在该公式体系中，第一项代表"传导"传热，第二项考虑了"对流"效应。值得注意的是，随着流体速度的降低，对流贡献逐渐变得微不足道。

需要注意的是，在推导传质和传热形状因子的广义表达式时，未对基质块的边界条件、材料属性或几何形状施加任何特定限制。相反，这些因素通过代表性单元问题的数值解来影响所计算的形状因子。因此，所得形状因子不是预设的常数，而是根据系统在真实几何构型、材料属性和边界条件下的单元瞬态响应来确定的。在简化条件下，所提出的公式体系可退化为常规表达式。例如，如果在方程(17)中忽略变形和热耦合效应，传质形状因子将退化为经典的水力形状因子。

同样地，当基质块内部的流体速度足够小时，方程(22)中的对流项消失，热传递形状因子变为以导热为主导。在足够晚的时间，随着系统趋近于伪稳态阶段，计算所得的形状因子可能收敛为与经典形状因子模型一致的常数值，具体取决于基质块的几何形状和边界条件。尽管如此，所提出的形状因子关系式在其推导中所采用的假设条件下是有效的。特别是，假设基质和裂隙连续体均处于完全饱和状态，基质被视为可渗透介质，并且在孔隙流体与固体骨架之间假定满足局部热平衡。此外，基质–裂隙交换项以两个连续体之间的平均压力差和温度差来表示。因此，对于涉及局部热非平衡、非饱和或多相流条件，或储层网格块内存在强烈非均匀压力场或温度场的应用，该公式可能需要修改。在这种情况下，单一的平均基质–裂隙压力差或温度差可能不足以准确表征局部交换过程。

**3. 控制方程的有限元离散化**

为了推导控制方程(1)、(2)、(3)、(7)和(8)的弱形式，采用Khoei（2015）提出的标准有限元Galerkin方法。有限元离散化通过将未知场变量表示为节点值与整体形状函数的形式来实现：

$$\mathbf{u}=\mathbf{N}_u\mathbf{U},\; p_m=\mathbf{N}_p\mathbf{p}_m,\; p_f=\mathbf{N}_p\mathbf{p}_f \tag{23}$$

$$T_m=\mathbf{N}_T\mathbf{T}_m,\; T_f=\mathbf{N}_T\mathbf{T}_f$$

其中$\mathbf{N}_u$、$\mathbf{N}_p$和$\mathbf{N}_T$分别表示位移、压力和温度形状函数矩阵。此外，$\mathbf{U}$、$\mathbf{p}_m$、$\mathbf{p}_f$、$\mathbf{T}_m$和$\mathbf{T}_f$分别表示节点位移、基质压力、裂隙压力、基质温度和裂隙温度的向量。通过施加有限元离散化，控制方程的弱形式可表示为：

$$\boldsymbol{\psi}_u=\mathbf{M}\ddot{\mathbf{U}}+\mathbf{K}\mathbf{U}-\mathbf{Q}\mathbf{p}_m-\mathbf{W}\mathbf{T}_m-\mathbf{f}_{\text{ext}}=\mathbf{0}$$

$$\boldsymbol{\psi}_{p_m}=\mathbf{Q}^T\dot{\mathbf{U}}+\mathbf{H}_m\mathbf{p}_m+\mathbf{S}_m\dot{\mathbf{p}}_m+\mathbf{H}'(\mathbf{p}_m-\mathbf{p}_f)-\mathbf{R}_m\dot{\mathbf{T}}_m=\mathbf{0}$$

$$\boldsymbol{\psi}_{p_f}=\mathbf{H}_f\mathbf{p}_f+\mathbf{S}_f\dot{\mathbf{p}}_f-\mathbf{H}'(\mathbf{p}_m-\mathbf{p}_f)-\mathbf{R}_f\dot{\mathbf{T}}_f-\mathbf{q}_{p,\text{ext}}=\mathbf{0}$$

$$\boldsymbol{\psi}_{T_m}=\mathbf{C}_m\mathbf{T}_m+\mathbf{L}_m\dot{\mathbf{T}}_m+\mathbf{H}_S(\mathbf{T}_m-\mathbf{T}_f)=\mathbf{0} \tag{24}$$

$$\boldsymbol{\psi}_{T_f}=\mathbf{C}_f\mathbf{T}_f+\mathbf{L}_f\dot{\mathbf{T}}_f-\mathbf{H}_S(\mathbf{T}_m-\mathbf{T}_f)-\mathbf{q}_{T,\text{ext}}=\mathbf{0}$$

其中$\mathbf{f}_{\text{ext}}$、$\mathbf{q}_{p,\text{ext}}$和$\mathbf{q}_{T,\text{ext}}$分别表示外力、流体通量和热通量向量。上述方程组表示双重孔隙度框架下热–水–力耦合问题的耦合方程组，其中各向量和矩阵定义如下：

$$\mathbf{M}=\int_V\mathbf{N}_u^T\rho\mathbf{N}_u\,dV,\quad \mathbf{K}=\int_V\mathbf{B}_u^T\mathbf{D}\mathbf{B}_u\,dV$$

$$\mathbf{Q}=\int_V\mathbf{B}_u^T\alpha_m\mathbf{N}_p\,dV,\quad \mathbf{W}=\int_V\mathbf{B}_u^T\frac{1}{3}\beta_m\mathbf{N}_T\,dV$$

$$\mathbf{f}_{\text{ext}}=\int_V\mathbf{N}_u^T\rho\mathbf{b}\,dV+\int_{\Gamma_t}\mathbf{N}_u^T\mathbf{t}\,d\Gamma$$

$$\mathbf{H}_m=\int_V\mathbf{B}_p^T\mathbf{k}_m\mathbf{B}_p\,dV,\quad \mathbf{S}_m=\int_V\mathbf{N}_p^TC_m\phi_m\mathbf{N}_p\,dV$$

$$\mathbf{H}'=\int_V\mathbf{N}_p^T s\mathbf{k}_m\mathbf{N}_p\,dV,\quad \mathbf{R}_m=\int_V\mathbf{N}_p^T\beta_w\mathbf{B}_T\,dV$$

$$\mathbf{H}_f=\int_V\mathbf{B}_p^T\mathbf{k}_f\mathbf{B}_p\,dV,\quad \mathbf{S}_f=\int_V\mathbf{N}_p^TC_f\phi_f\mathbf{N}_p\,dV$$

$$\mathbf{R}_f=\int_V\mathbf{N}_p^T\beta_w\mathbf{B}_T\,dV,\quad \mathbf{q}_{p,\text{ext}}=\int_{\Gamma_t}\mathbf{N}_p^T\bar{q}_p\,d\Gamma \tag{25}$$

$$\mathbf{H}_S=\int_V\mathbf{N}_T^Ts\lambda_m\mathbf{N}_T\,dV,\quad \mathbf{q}_{T,\text{ext}}=\int_{\Gamma_t}\mathbf{N}_T^T\bar{q}_T\,d\Gamma$$

其中$\mathbf{B}_u$表示应变–位移矩阵，$\mathbf{B}_p$表示由压力形状函数的空间梯度获得的导数矩阵。在上述关系式中，矩阵$\mathbf{M}$、$\mathbf{K}$、$\mathbf{Q}$和$\mathbf{W}$分别与力学惯性、刚度、水–力耦合和热–力耦合项相关。此外，矩阵$\mathbf{H}_m$和$\mathbf{H}_f$表示水力渗透率矩阵，$\mathbf{S}_m$和$\mathbf{S}_f$表示水力储容矩阵，$\mathbf{R}_m$和$\mathbf{R}_f$对应与热膨胀效应相关的热–水力耦合矩阵。矩阵$\mathbf{H}'$表示基质–裂隙流体交换耦合项。下标'm'和'f'分别指代基质和裂隙连续体。

此外，当对流占主导时，传统的Galerkin格式可能变得不稳定。在这种情况下，温度场中可能出现虚假的空间振荡。在热传递分析中，描述平流与扩散热传递之比的爱克哈特数可表示为：

$$\text{Pé}=\frac{(\rho c)_w v_w h}{\lambda} \tag{26}$$

其中$h$表示沿流动方向的单元长度。爱克哈特数用于衡量平流热传递的重要性。爱克哈特数大于一表明以平流为主导，而爱克哈特数小于一则表明以导热为主导的热流动。为了消除高度对流流动中的数值振荡并获得可靠结果，此处采用Petrov-Galerkin迎风格式。该方法提供了一种鲁棒的手段来有效处理数值不稳定性。因此，权函数被选取为不同于单元形状函数的形式（$W_I^e \neq N_I^e$），引入相对于流场的方向性不对称。由此，对于方程组(24)定义以下矩阵：

$$\mathbf{C}_m=\int_V\frac{d\mathbf{W}_T}{d\mathbf{x}}\lambda_m\mathbf{B}_T^T\,dV+\int_V\mathbf{W}_T^T(\rho c)_w\mathbf{v}_{w,m}\mathbf{B}_T^T\,dV$$

$$\mathbf{L}_m=\int_V\mathbf{W}_T^T(\rho c)_m\mathbf{N}_T\,dV$$

$$\mathbf{C}_f=\int_V\frac{d\mathbf{W}_T}{d\mathbf{x}}\lambda_f\mathbf{B}_T^T\,dV+\int_V\mathbf{W}_T^T(\rho c)_w\mathbf{v}_{w,f}\mathbf{B}_T^T\,dV \tag{27}$$

$$\mathbf{L}_f=\int_V\mathbf{W}_T^T(\rho c)_f\mathbf{N}_T\,dV$$

其中矩阵$\mathbf{C}_m$和$\mathbf{C}_f$分别表示基质和裂隙连续体的导热和 convective heat transfer 贡献，而$\mathbf{L}_m$和$\mathbf{L}_f$表示相应的热储容矩阵。此外，$\mathbf{B}_T$是由温度形状函数的空间梯度构建的导数矩阵，权函数$\mathbf{W}_T$基于Petrov-Galerkin格式使用最优迎风参数定义，表示为：

$$\alpha_{\text{opt}}=\coth\frac{\text{Pé}}{2}-\frac{2}{\text{Pé}}$$

关于Petrov-Galerkin迎风技术、权函数的推导以及最优加权因子的定义，详见附录A，参考Brooks和Hughes（1982）、Cui等人（2018）以及Taghvaei和Khoei（2025）的研究。

此外，方程组(24)采用已建立的Newmark格式（Khoei, 2015）进行时间离散化，将$t_{n+1}=t_n+\Delta t$时刻的未知场变量与$t_n$时刻的已知变量通过以下关系联系起来：

$$\ddot{\mathbf{U}}_{n+1}=a_0(\mathbf{U}_{n+1}-\mathbf{U}_n)-a_2\dot{\mathbf{U}}_n-a_4\ddot{\mathbf{U}}_n$$

$$\dot{\mathbf{U}}_{n+1}=a_1(\mathbf{U}_{n+1}-\mathbf{U}_n)-a_3\dot{\mathbf{U}}_n-a_5\ddot{\mathbf{U}}_n$$

$$\dot{\mathbf{p}}_{m,n+1}=a_6(\mathbf{p}_{m,n+1}-\mathbf{p}_{m,n})-a_7\dot{\mathbf{p}}_{m,n}$$

$$\dot{\mathbf{p}}_{f,n+1}=a_6(\mathbf{p}_{f,n+1}-\mathbf{p}_{f,n})-a_7\dot{\mathbf{p}}_{f,n}$$

$$\dot{\mathbf{T}}_{m,n+1}=a_6(\mathbf{T}_{m,n+1}-\mathbf{T}_{m,n})-a_7\dot{\mathbf{T}}_{m,n} \tag{28}$$

$$\dot{\mathbf{T}}_{f,n+1}=a_6(\mathbf{T}_{f,n+1}-\mathbf{T}_{f,n})-a_7\dot{\mathbf{T}}_{f,n}$$

其中

$$a_0=\frac{1}{\beta\Delta t^2},\quad a_1=\frac{\gamma}{\beta\Delta t},\quad a_2=\frac{1}{\beta\Delta t},\quad a_3=\frac{\gamma}{\beta}-1 \tag{29}$$

$$a_4=\frac{1}{2}\beta-1,\quad a_5=\frac{\Delta t}{\gamma}\left(\frac{1}{2\beta}-1\right),\quad a_6=\frac{1}{\theta\Delta t},\quad a_7=\frac{1}{\theta}-1$$

其中$\beta$、$\gamma$和$\theta$表示Newmark参数。最后，方程组可通过Newton-Raphson方法线性化，得到以下代数方程组：

$$\boldsymbol{\Psi}_{n}^{i+1}=\boldsymbol{\Psi}_{n}^{i}+\mathbf{J}_n^i\cdot\Delta\mathbf{X}_n^i$$

$$\boldsymbol{\Psi}=\begin{bmatrix}\boldsymbol{\psi}_u\\ \boldsymbol{\psi}_{p_m}\\ \boldsymbol{\psi}_{p_f}\\ \boldsymbol{\psi}_{T_m}\\ \boldsymbol{\psi}_{T_f}\end{bmatrix},\quad \mathbf{X}=\begin{bmatrix}\mathbf{U}\\ \mathbf{p}_m\\ \mathbf{p}_f\\ \mathbf{T}_m\\ \mathbf{T}_f\end{bmatrix} \tag{30}$$

其中$\boldsymbol{\Psi}$和$\mathbf{X}$分别表示残差向量和主变量的节点值向量。Jacobian矩阵$\mathbf{J}$在第$i$次迭代中可表示为：

$$\mathbf{J}=\begin{bmatrix}\dfrac{\partial\boldsymbol{\psi}_u}{\partial\mathbf{U}} & \dfrac{\partial\boldsymbol{\psi}_{p_m}}{\partial\mathbf{U}} & \dfrac{\partial\boldsymbol{\psi}_u}{\partial\mathbf{p}_m} & \dfrac{\partial\boldsymbol{\psi}_{p_m}}{\partial\mathbf{p}_m} & \cdots\\[10pt] \vdots & & & & \ddots\end{bmatrix} \tag{31}$$

将方程(28)中的变量代入空间离散化关系式(24)后，Jacobian矩阵可表示为：

$$\mathbf{J}=\begin{bmatrix} a_0\mathbf{M}+\mathbf{K}-\mathbf{Q}\mathbf{0} & -\mathbf{W}\mathbf{0} & \mathbf{0} & \mathbf{0} & \mathbf{0}\\[6pt] a_1\mathbf{Q}^T & \mathbf{H}_m+\mathbf{H}'+a_6\mathbf{S}_m & \mathbf{H}' & -a_6\mathbf{R}_m & \mathbf{0}\\[6pt] \mathbf{0} & -\mathbf{H}' & \mathbf{H}_f+\mathbf{H}'+a_6\mathbf{S}_f & \mathbf{0} & -a_6\mathbf{R}_f\\[6pt] \mathbf{0} & \mathbf{0} & \mathbf{0} & \mathbf{C}_m+\mathbf{H}_S+a_6\mathbf{L}_m & -\mathbf{H}_S\\[6pt] \mathbf{0} & \mathbf{0} & \mathbf{0} & -\mathbf{H}_S & \mathbf{C}_f+\mathbf{H}_S+a_6\mathbf{L}_f \end{bmatrix} \tag{32}$$

显然，问题的不同物理场之间存在强耦合关系。需要注意的是，方程(32)中的Jacobian矩阵通常是非对称的。这种非对称性源于控制方程的耦合特性、非对称的非对角耦合项、能量平衡方程中对流项的存在，以及Petrov-Galerkin权函数的使用。因此，需要非对称矩阵求解器来求解所得的线性化方程组。在本实现中，采用了一种适用于非对称矩阵的求解器。

**4. 耦合THM双重孔隙度模型的计算框架**

传统的双重孔隙度模型通常采用恒定的热传递和质量传递形状因子，忽略了物性参数、基质构型、边界条件以及问题瞬态动态的影响。这种简化可能导致显著的不准确，特别是在流体流动、热传递和力学变形动态交互的全耦合热–水–力系统中。此外，热传递中的瞬态期在地热系统等应用中至关重要。为解决这些问题，为裂隙域的全耦合THM分析引入了一种两步计算算法。该算法采用随时间变化的热传递和质量传递形状因子，这些形状因子以耦合方式计算，考虑了当前域的特征，捕捉了质量传递（流体流动）、热传递和力学变形之间的相互依存关系。因此，所得的质量传递和热传递形状因子能够更准确地表征裂隙系统。两步程序概述如下：

- **单元分析**：选取一个代表性单元（代表性体积元，RVE），反映主域的基质构型和THM特性。这些特性包括基质和裂隙的渗透率、压缩性、热导率、比热以及其他相关参数。采用有限元方法，在一致的边界条件下求解该单元的全耦合THM边界值问题。在每个时间步，利用方程(22)计算瞬态热传递形状因子，考虑可渗透基质中的导热和对流机制。同时，利用方程(17)在每个时间步计算瞬态质量传递形状因子，考虑水力、热力与力学过程之间的相互作用。

- **主问题求解**：在整个域上使用相同的时间步求解主要的THM边界值问题。在先前步骤中计算的瞬态质量传递和热传递形状因子用于表征基质–裂隙之间的相互作用。该方法确保了耦合系统物理一致性的表征。

这种两步程序的主要优势在于，代表性单元问题远小于全裂隙分辨的直接数值模拟（DNS）模型。因此，所提出的程序能够在保持双重孔隙度框架计算效率的同时，在分析过程中更新传递系数。该计算算法的示意图和相关流程图分别见图3和图4。



**5. 数值模拟结果**

为了展示两步THM计算算法在全耦合THM建模裂隙域中的准确性和性能，本文数值求解了若干算例。在第一组算例中，分析了方形和非方形块体，并施加了各种边界条件。通过将瞬态质量传递和热传递形状因子所得结果与细网格单孔隙度结果进行比较，评估了结果的准确性。下一个算例涉及包含方形块体的裂隙土柱，旨在展示THM算法在建模较大域时的高效性和准确性。最后，建模了两个场地尺度地热储层问题，以展示所提算法在处理涉及非均质区域的实际应用方面的能力。在本节中展示的所有数值算例中，均未观察到与方程(11)中Heaviside切换项相关的收敛困难。尽管如此，对于涉及基质–裂隙质量传递方向频繁反转的问题，可采用平滑Heaviside近似以增强数值鲁棒性。此外，在每个时间步监控了Newton-Raphson收敛情况，未观察到与Jacobian矩阵非对称性相关的收敛问题。

**5.1. 方形和非方形块体的THM分析**

第一组算例选取用于验证方形和非方形块体THM分析中瞬态质量传递和热传递形状因子的准确性。为此，选取了正方形和等腰直角三角形几何形状的块体，如图5所示。热–水–力分析的材料属性列于表1中。这些问题采用双重孔隙度方法模拟，所得结果与细网格单孔隙度模型（作为直接数值模拟，DNS）的结果进行比较。在双重孔隙度模型中，裂隙和基质域均假定为与总域面积相同。双重孔隙度模型的等效属性遵循Oda（1986）的工作，基于各域的体积变化计算。方形和三角形块体中相应的裂隙面积分别为0.059 m2和0.05 m2。双重孔隙度分析采用双线性有限元进行。每个块体所使用的DNS网格也见图5。如前所述，文献中通常采用仅依赖于基质块几何尺寸的恒定质量传递和热传递形状因子。此类模型未考虑基质构型或本研究所考虑的THM物性参数。在本研究中，采用了Warren和Root（1963）、Kazemi等人（1976）以及Lim和Aziz（1995）提出的三种广泛使用的恒定形状因子模型。表2总结了这些模型的无量纲形式。


**表1.** 方形和三角形块体THM问题中的材料属性（DNS模型）。

| 材料属性 | 基质 | 裂隙 |
|---------|------|------|
| 固有渗透率$k$(m2) | $9\times10^{-17}$ | $9\times10^{-13}$ |
| 弹性模量$E$(Pa) | $3\times10^7$ | $3\times10^7$ |
| 泊松比$\nu$ | 0.2 | 0.2 |
| 压缩性系数$C_\phi$(1/Pa) | $1.4286\times10^{-8}$ | $1.4286\times10^{-8}$ |
| 流体动力粘度$\mu_f$(Pa·s) | $10^{-3}$ | $10^{-3}$ |
| Biot系数$\alpha$ | 1 | 1 |
| 密度$\rho$(kg/m3) | 2000 | 1000 |
| 比热$c$(kJ/kg·°C) | 1 | 4 |
| 热导率$\lambda$(W/m·°C) | 3 | 0.6 |

**表2.** 文献中提出的无量纲恒定形状因子($sA_m$)。

| 研究者 | 二维流动 |
|--------|---------|
| Warren和Root（1963） | 32 |
| Kazemi等人（1976） | 8 |
| Lim和Aziz（1995） | 18.17* 或 $2\pi^2$* |

*圆柱形块体。

第一个算例考虑一个被裂隙区域包围的等腰直角三角形基质块。初始基质压力和温度设为零，而周围裂隙区域中的压力和温度分别保持在$p=20$ kPa和$T=100$°C恒定。在这些条件下，由于基质–裂隙之间的质量和热交换，基质块内的压力和温度随时间增加，最终趋近平衡。利用方程(17)和(22)计算的质量传递和热传递形状因子绘制于图6中。可以看出，质量传递形状因子与热传递形状因子不同，且在瞬态阶段和稳态阶段均表现出更大的值。图7展示了该三角形块体中基质平均温度和压力随时间的变化。如图所示，使用双重孔隙度方法所得结果在趋势和数值上均与细网格模拟（DNS）一致。

相反，将恒定且相等的形状因子应用于热量和质量传递函数的计算会导致显著误差，尤其是在瞬态阶段。下载：下载高分辨率图像（133KB）下载：下载全尺寸图像图6. 在恒压和恒温边界条件下，三角形岩块的质量与热交换形状因子随时间的变化。下载：下载高分辨率图像（299KB）下载：下载全尺寸图像图7. 在恒压和恒温边界条件下，三角形岩块平均基质温度和压力随时间的变化；不同模型之间的对比。在下一个示例中，考虑如图5所示的方形基质岩块。初始基质压力和温度设为零，而周围裂隙区域的压力和温度则规定为随时间线性增加，分别为p=20×101+10??tkPa和T=100×101+10??t°C。相应的瞬态质量和热传递形状因子绘制在图8中。显然，这些形状因子并不相等，且它们稳定在比三角形岩块更低的值。图9对比了使用这些瞬态形状因子获得的双重孔隙度结果与DNS（直接数值模拟）及恒定双重孔隙度模型的结果。可以观察到，图8中的瞬态质量和热形状因子能够准确模拟裂隙与基质区域之间的质量和热交换。另一方面，传统模型在瞬态阶段对压力场和温度场均表现出误差。下载：下载高分辨率图像（137KB）下载：下载全尺寸图像图8. 在线性压力和线性温度边界条件下，方形岩块的质量与热交换形状因子随时间的变化。下载：下载高分辨率图像（294KB）下载：下载全尺寸图像图9. 在线性压力和线性温度边界条件下，方形岩块平均基质温度和压力随时间的变化；不同模型之间的对比。方形岩块也在阶梯式加载边界条件下进行研究。载荷施加于岩块顶部，公式为：（33）f=0.15tkN，当t<100s时；f=15kN，当t≥100s时。此外，裂隙温度设为：（34）T=50(1+0.01t)°C，当t<100s时；T=100°C，当t≥100s时。瞬态耦合的质量与热传递形状因子根据公式（17）、（22）计算，如图10所示。在这种情况下，瞬态期间流体形状因子和热形状因子之间的差异更为明显。然而，它们在稳态时的常数值相似。使用不同模型获得的基质压力和温度随时间的变化绘制在图11中。显然，与传统模型不同，采用计算所得瞬态形状因子的双重孔隙度方法的结果与DNS模型吻合良好。下载：下载高分辨率图像（119KB）下载：下载全尺寸图像图10. 在梯形加载和线性温度边界条件下，方形岩块的质量与热交换形状因子随时间的变化。下载：下载高分辨率图像（294KB）下载：下载全尺寸图像图11. 在梯形加载和线性温度边界条件下，方形岩块平均基质温度和压力随时间的变化；不同模型之间的对比。根据上述结果，可以推断，公式（17）、（22）中引入的耦合形状因子能够准确捕捉不同边界条件下方形和非方形岩块中基质与裂隙区域之间质量和热交换的瞬态行为。此外，还观察到质量与热传递形状因子的值在瞬态和稳态阶段通常不同。因此，对两个形状因子应用相同的值会导致显著误差。在下一组示例中，将评估两步THM算法在大规模问题建模中的准确性和性能。5.2 含裂隙土柱的THM分析为了评估所提出的两步热-水-力学框架的准确性和性能，研究了含裂隙土柱问题。该土柱的几何形状、代表性单元体以及边界条件如图12所示。土柱的热-水-力学性质列于表3中。双重孔隙度模型的等效热和液压性质也根据Oda（1986）进行了计算。含裂隙土柱的边界条件如图12所示。在土柱顶部和底部边缘，规定Dirichlet型压力条件，分别为p=0和p=2MPa。裂隙区域维持恒温100°C，而侧边界假设为热绝缘和液压绝缘。对于力学边界条件，底部边缘在竖直方向固定，侧边界在水平方向约束。本例采用1×1m2的单元体，如图12(d)所示。基于这些边界条件，计算框架的第一步是在单元体上施加等效条件进行定义。相应地，温度设为100°C，裂隙区域压力指定为2/7MPa。基质区域具有渗透性，因此导热和对流机制均对基质和裂隙域中的传热过程有贡献。为验证所提方法的准确性，结果与DNS结果进行了对比。为此，DNS模型采用了5789个双线性单元体和6143个节点，如图12(b)所示，而双重孔隙度模拟则使用140个单元体和213个节点。所选单元体的相应瞬态质量和热传递形状因子绘制在图13中。下载：下载高分辨率图像（337KB）下载：下载全尺寸图像图12. 土柱问题；(a) 几何形状和边界条件，(b) 采用5789个双线性单元体的直接数值模拟（DNS），(c) 采用140个双线性单元体的双重孔隙度模型，(d) 两步THM算法单元体的有限元网格。表3. 含裂隙土柱的材料性质（DNS模型）。材料性质基质裂隙固有渗透率k(m2)低：9×10?1?高：9×10?1?9×10?13弹性模量E(Pa)3×10?3×10?泊松比ν0.20.2压缩系数Cφ(1/Pa)1.4286×10??1.4286×10??流体动力黏度μf(Pa·s)10?310?3Biot系数α11密度ρ(kg/m3)20001000比热c(kJ/kg·°C)14导热系数λ(W/m·°C)30.6下载：下载高分辨率图像（115KB）下载：下载全尺寸图像图13. 土柱单元体的质量与热交换形状因子随时间的变化。图14、图15展示了如图12中所示的B点和C点处基质压力和温度场以及A点处位移随时间的变化，采用不同方法进行分析。可以观察到，使用当前计算方法获得的结果与DNS模型结果吻合良好。相反，传统双重孔隙度模型在瞬态阶段的压力场和温度场表现出偏差。另一个重要发现是，包含QHW项——该术语考虑了非等温水进入基质或裂隙域的效果——可以提高结果的准确性。相反，忽略该热源会导致瞬态阶段温度的高估。为进一步检验该术语对温度预测的影响，考虑了如表3所述的低基质渗透率和高基质渗透率两种情况。图16对比了来自DNS解和当前THM算法在C点的温度演化结果，分别包含和排除QHW。结果清楚地表明，纳入QHW可以提高两种情况下热传递预测的准确性。此外，随着基质渗透率的增加，其效果更加显著，此时非等温水的交换增强。下载：下载高分辨率图像（700KB）下载：下载全尺寸图像图14. B点和C点的基质压力和温度随时间的变化；不同模型的对比。下载：下载高分辨率图像（184KB）下载：下载全尺寸图像图15. A点的位移随时间的变化；不同模型的对比。下载：下载高分辨率图像（361KB）下载：下载全尺寸图像图16. 不同基质渗透率下C点的温度演化；当前模型包含和排除Qhw热耦合的对比。为定量评估不同双重孔隙度方法在瞬态响应中的准确性，计算了从t=0到t=8×10?s时间区间内的平均百分比误差，并报告在表4中。在此时间段之后，系统接近稳态条件。为评估准确性与计算效率之间的权衡，不同方法的CPU时间也在表5中进行了对比。报告的CPU时间在配备两个AMD EPYC 9334 32核处理器（运行频率为2.70 GHz）和1.00 TB内存的工作站上获得。CPU时间对比表明，所提出的具有瞬态形状因子的双重孔隙度公式所需的计算成本高于具有恒定形状因子的传统双重孔隙度模型，这是因为需要额外的计算来评估瞬态质量和热传递系数。然而，对于所考虑的含裂隙土柱问题，所提出的方法在计算效率上仍优于细网格单一孔隙度（DNS）模型。误差分析进一步表明，在瞬态期间，所提出的方法对监测的压力、温度和位移场的误差小于传统的恒定形状因子模型。总体而言，这些结果表明，所提出的公式在相对于经典双重孔隙度模型计算成本适度增加的情况下，提高了双重孔隙度方法的瞬态准确性，同时相对于本示例中的DNS仍保持明显的计算优势。表4. 含裂隙土柱问题中不同方法的平均百分比误差对比。双重孔隙度（DP）方法平均百分比误差（0-8×10?s）pm@Bpm@CTm@BTm@Cuy@ADP – 本研究（包含Qhw）1.53%3.18%3.57%7.09%4.48%DP – 本研究（不包含Qhw）1.53%3.18%9.26%11.41%4.48%DP – （Warren-Root模型）6.65%6.09%12.83%13.22%6.15%DP – （Lim-Aziz模型）12.57%12.36%11.77%10.04%13.23%表5. 含裂隙土柱问题中不同方法的CPU时间对比。方法CPU时间（小时）直接数值模拟（DNS）方法13.57具有瞬态形状因子的双重孔隙度方法（本研究）5.18具有恒定形状因子的双重孔隙度方法（传统模型）1.085.3 地热储层采热的THM模拟本节考虑了两个地热储层案例研究，以进一步评估所提出的两步THM双重孔隙度框架在建模具有不同裂隙构型的非均匀裂隙域方面的效率和能力。在通过单一岩块和含裂隙土柱示例与DNS解进行验证后，这些储层尺度模拟用于展示所提方法在具有复杂裂隙网络、多样化基质岩块几何形状和非均匀液压条件的工程尺度系统中的应用。此类大规模储层的直接细网格（DNS）模拟计算成本高昂，因为需要显式离散化大量裂隙和基质岩块。因此，第一个案例研究重点关注裂隙构型的影响，而第二个案例研究则考察液压和热非均匀性对储层行为的影响。5.3.1 裂隙构型对地热储层行为的影响第一个案例研究旨在评估所提框架捕捉裂隙构型对地热储层中耦合THM响应影响的能力。考虑的储层几何形状如图17所示，尺寸为L=200m和l=5m。相应的材料性质总结在表6中。计算域采用由20×20双线性有限元组成的均匀网格进行离散化。施加的边界条件如图17所示。冷水通过注入井以Tin=20°C的温度和pin=13 MPa的压力注入，而生产井压力规定为pout=3MPa。所有其余外边界假设为不渗透和热绝缘，阻止流体和热通量穿过域边界。对于力学问题，左右边界在水平方向约束（ux=0），而顶部和底部边界在竖直方向约束（uy=0）。考虑了两种裂隙构型情景以研究裂隙几何形状对储层性能和THM行为的影响：•均匀裂隙构型：储层由均匀的方形基质岩块组成，•非均匀裂隙构型：储层被分为四个具有不同裂隙构型的区域，如图17所示。区域I包含与均匀情况相同的方形基质岩块。区域II和III分别由梯形和三角形岩块组成。区域IV由比其它区域裂隙间距更大的较小方形岩块组成。下载：下载高分辨率图像（182KB）下载：下载全尺寸图像图17. 五点法地热储层的四分之一模型：构型、不同的储层区域及相应的单元体，以及由400个双线性单元体组成的双重孔隙度网格。表6.

# 断裂构型地热储层中不同区域单元格的物性参数

| 物性参数 | 基质 | 裂缝 |
|---|---|---|
| 固有渗透率 k (m2) | 低：8.5×10?1? | 区域I、II、III：1.9×10?13；区域IV：5.28×10?13 |
| 弹性模量 E (Pa) | 4×101? | 1×101? |
| 泊松比 ν | 0.2 | 0.2 |
| 可压缩系数 C_φ (1/Pa) | 1.5×10?? | 1.5×10?? |
| 流体动力粘度 μ_f (Pa·s) | 10?3 | 10?3 |
| Biot 系数 α | 1 | 1 |
| 密度 ρ (kg/m3) | 2000 | 1000 |
| 比热 c (kJ/kg·°C) | 1.25 | 4.18 |
| 导热系数 λ (W/m·°C) | 2.85 | 0.6 |

各区域对应的单元格及其裂缝面积和基质面积如图17和表7所示。根据各区域的断裂构型，区域II和区域III的裂缝渗透率几乎相同，且高于区域I。区域IV由于裂缝数量更多、裂缝宽度更大，因此具有最高的裂缝渗透率。四个区域对应的瞬态质量传递和热量传递形状因子如图18所示。需要注意的是，计算得到的区域I的形状因子也适用于均匀案例。如图18所示，位于不同区域的单元格的质量传递和热量传递形状因子存在差异，表明基质-裂缝交换过程强烈受到局部基质块构型和演化的THM（热-水-力）条件的影响。这些结果表明，单一等效方形块表征或常数形状因子的使用可能不足以准确捕捉基质与裂缝连续体之间质量和热量交换的空间变异性。因此，考虑各区域特有的基质块几何形状和瞬态传递特性对于提高双重孔隙度THM模型的预测能力至关重要。

图19比较了利用图18中的形状因子，在所提出的两步THM方法中，均匀案例和非均匀案例在各时间步获得的压力等值线。显然，在非均匀案例中，由于区域II、III和IV具有较高的裂缝渗透率，压力在不同时间步上扩散速度更快。另一方面，在均匀案例中，压力在整个域中均匀传播。非均匀案例中区域I的压力耗散受到其低渗透率的限制。然而，在区域IV中，高渗透率促进了向生产井的流动，导致该区域压力降低更为有效，从而在该区域形成更宽的等值线间距。

**表7 与不同区域相关的各单元格的裂缝面积和基质面积**

| 块体形状 | 裂缝面积 | 基质面积 |
|---|---|---|
| 方形块体（区域I） | A_f = 1.475 m2 | A_m = 23.525 m2 |
| 方形块体（区域IV） | A_f = 2.4375 m2 | A_m = 22.56 m2 |
| 三角形块体 | A_f = 2.5 m2 | A_m = 22.5 m2 |
| 梯形块体 | A_f = 2.25 m2 | A_m = 22.75 m2 |

**图18.** 断裂构型地热储层中不同区域单元格的质量交换和热量交换形状因子随时间的变化（均匀和非均匀案例研究）。

**图19.** 断裂构型地热储层（均匀和非均匀案例研究）在不同时间步的基质压力等值线。

在图20中，展示了均匀案例和非均匀案例的基质温度等值线。从该图可以看出，具有均匀断裂构型的均匀案例在整个储层中表现出更为均匀分布的温度场。相比之下，非均匀案例在所有报告的时间步上均表现出更快的热锋传播。这种行为可归因于区域II、III和IV较高的渗透率，增强了流体流动，从而增加了这些区域内的对流热量传递。此外，非均匀案例中的温度等值线明显向区域IV（渗透率最高的区域）发生畸变。热锋的优先推进突出了渗透率非均匀性与热量传递之间的强耦合关系，表明了断裂网络构型对储层热行为演变的显著影响。这种行为可归因于区域IV中高渗透率的裂缝，它们在储层中充当优先流动通道。渗透率的增加降低了水力阻力，导致流体速度提高和对流热量传递增强。因此，热能通过该区域更高效地传递，显著影响了周围的温度场。由此产生的不对称和弯曲的温度等值线清楚地表明了渗透率非均匀性与热量传递过程之间的强耦合关系。这些结果表明，所提出的THM框架能够准确捕捉断裂地热储层中非均匀断裂特性、流体流动和热量传递之间的复杂相互作用。

**图20.** 断裂构型地热储层（均匀和非均匀案例研究）在不同时间步的基质温度等值线。

图21展示了均匀案例和非均匀案例的水平位移等值线。在瞬态过程的早期阶段，变形主要由孔隙压力变化控制。因此，非均匀案例由于压力在高渗透率区域传播更快，表现出更宽的空间变形分布和更大的位移量。随着模拟的进行，热锋通过储层推进，降温引起的热收缩开始影响力学响应，导致位移减小。在非均匀案例中，冷锋传播更快，并向区域IV呈现明显的弯曲，反映了与其高渗透率相关的优先流动和热量传递。这些结果突出了断裂地热储层中强烈的热-水-力耦合，并展示了所提出框架在捕获非均匀渗透率条件下流体流动、热量传递和固体变形之间复杂相互作用的能力。

**图21.** 断裂构型地热储层（均匀和非均匀案例研究）在不同时间步的水平位移等值线。

为了进一步评估均匀案例和非均匀案例在断裂构型储层场景中的性能，利用方程(35)计算的生产水温度和热效率曲线如图22所示。

$$\text{热效率} = \frac{\text{热量提取速率}}{\text{可用热量}} = \frac{T_{out} - T_{in}}{T_{ref} - T_{in}} \tag{35}$$

结果表明，均匀案例和非均匀案例在热行为上存在显著差异。在均匀案例中，裂缝和基质的热-水-力特性在整个域中均匀分布，生产井温度逐渐降低，40年运行后降至65.2°C。相比之下，非均匀案例的特征是具有不同渗透率的区域，表现出更早的热突破和出口温度的急剧下降，在相同时间内降至23.2°C。这种加速冷却主要是由区域IV中高渗透率的裂缝造成的，它们为注入流体创造了优先流动通道。由此产生的流体速度增加使冷注入水更快地到达生产井。这种通道效应加速了冷锋从注入井到生产井的对流传输，超过了周围基质通过传导补充热量的能力，从而导致更早的热突破。因此，储层的热性能受到渗透率分布的显著影响。在均匀案例中，热效率在整个40年运行期间保持在67%以上，表明持续的热量提取效率。相比之下，非均匀案例表现出热效率的快速下降，在模拟结束时降至约11%。这些结果突出了渗透率非均匀性对热突破和储层耗竭的主导影响，强调了在预测地热储层长期性能时准确表征非均匀断裂网络的重要性。

**图22.** 断裂构型地热储层中均匀和非均匀案例在40年注水期间的(a)出口水温和(b)热效率随时间的变化。

总体而言，压力、温度和位移场中观察到的差异表明，依赖于区域的断裂构型、等效渗透率和瞬态传递形状因子共同控制了储层尺度的THM响应。特别是，非均匀案例中观察到的倾斜等值线图案无法通过假设由单一等效方形基质块和均匀基质-裂缝传递特性所代表的均匀储层来再现。热量传递形状因子控制基质与裂缝之间热交换的速率；因此，它在不同区域之间的变化影响冷锋的传播，并可能显著改变热突破行为。同样，质量传递形状因子控制基质与裂缝之间的压力平衡和流体交换，从而影响流体再分配以及整个储层中由此产生的对流热量传递。这些发现强调了在双重孔隙度THM模拟中考虑空间变化的基质块几何形状和传递特性的重要性。非均匀渗透率和依赖于区域的传递过程的联合作用进一步反映在生产性能中，如图23所示的流出速率和热量提取结果所证实的。

**图23.** 断裂构型地热储层中均匀和非均匀案例在25年注水期间的(a)每宽度流出速率和(b)热量提取速率随时间的变化。

图23展示了均匀案例和非均匀案例在25年运行期间的生产流量时间演化。非均匀案例达到约0.08 L/s/m的稳态流出速率，比均匀案例（0.048 L/s/m）高出约66%。这种增强的生产能力源于区域IV中高渗透率的裂缝，它们形成了优先流动通道，促进了注入井和生产井之间更有效的流体输送。该图还展示了两个案例相应的热量提取速率，假设储层宽度为100m。这些结果进一步揭示了断裂网络非均匀性对储层生产力和长期热性能的影响。如图所示，非均匀案例在运行早期表现出优越的性能，在前五年内达到约2.3 MW的峰值热量提取速率。这一数值比均匀案例高约64%，反映了高渗透率区域所促进的增强流体循环。然而，在运行约15年后，非均匀案例经历了热性能的逐步下降，到第25年热量提取速率降至1.27 MW。这种急剧下降主要与高渗透率区域中强化的对流热量传递有关，特别是区域IV。在这些区域建立的优先流动通道加速了冷注入水从注入井到生产井的运动，导致过早的热突破。因此，热能的提取速度超过了周围岩石基质通过传导热量传递进行补充的速度，导致储层加速热耗竭。相比之下，均匀案例表现出更为平缓的热量提取速率下降，仅在约24年运行后才出现明显的热耗竭。在约第24年观察到性能交叉点，在此之后，均匀案例的热量提取速率（约1.4 MW）超过了非均匀案例。总体而言，这些结果强调了地热储层中短期生产力和长期可持续性之间的权衡。虽然高渗透率区域增强了对流热量传递并增加了早期能源生产，但它们也加速了热突破和储层耗竭。相反，更为均匀的储层特性虽然产生较低的热量提取速率，但支持更稳定和可持续的长期能源生产。该案例研究的结果突出了渗透率非均匀性在控制热突破、热量提取性能和长期储层耗竭中的关键作用。

总体而言，这些发现表明了所提出THM框架在模拟具有复杂断裂几何形状、非均匀渗透率分布和空间变化的基质-裂缝交换过程的断裂地热储层方面的能力和有效性。该框架成功捕捉了流体流动、热量传递和固体变形之间的耦合相互作用，为工程尺度地热系统的分析和设计提供了强大的工具。

## 5.3.2 水力非均匀性和热非均匀性对地热储层行为的影响

为了进一步评估所提出THM算法在模拟非均匀域中的性能和能力，考虑了一个尺寸为100×100 m2的第二地热储层，如图17所示（L=100m，l=2.5m）。该储层的物理特性列于表8中。对该储层考虑了两种情景，即水力非均匀性和热-水力非均匀性。在水力非均匀性案例中，储层包含四个区域，这些区域仅在断裂构型上有所不同，导致各区域之间的裂缝渗透率存在差异，而其他THM特性（如导热系数和热容）则保持均匀。

相比之下，在热-水力非均质性算例中，IV区的裂隙区域导热系数和热容量均被设为其他区域的30%；因此，该算例同时存在水力非均质性和热学非均质性。该储层的初始压力设为4.5 MPa，初始温度设为90°C。在注入点处，水压和温度分别规定为$P_{\text{in}}$=5 MPa和$T_{\text{in}}$=20°C。在生产点处，压力和温度分别规定为$P_{\text{out}}$=3 MPa和$T_{\text{out}}$=90°C。其余水力、热学和力学边界条件与前文储层算例中采用的相同，即：外部边界对热量和流体通量绝热封闭；左右边界在水平方向受约束$u_x$=0；上下边界在垂直方向受约束$u_y$=0。这两个算例中不同区域的单元质量传递和热传递形状因子随时间的变化如图24所示。

表8：水力与热学非均质储层中不同区域相关单元的材料属性

| 材料属性 | 基岩 | 裂隙 |
|---|---|---|
| 固有渗透率 $k$ (m2) | 7.2×10?1? | I、II、III区：8.5×10?13；IV区：2.36×10?12 |
| 弹性模量 $E$ (Pa) | 4×10? | 3×10? |
| 泊松比 $\nu$ | 0.2 | 0.2 |
| 压缩系数 $C_\varphi$ (1/MPa) | 1.5×10?? | 1.5×10?? |
| 流体动力黏度 $\mu_f$ (Pa·s) | 10?3 | 10?3 |
| Biot系数 $\alpha$ | 1 | 1 |
| 密度 $\rho$ (kg/m3) | 2000 | 1000 |
| 比热容 $c$ (kJ/(kg·°C)) | 基准模型：1.25；IV区（低值）：0.375 | 4.2 |
| 导热系数 $\lambda$ (W/(m·°C)) | 基准模型：2.9；IV区（低值）：0.87 | 0.65 |
| 热膨胀系数 $\beta_s$ | 8×10?? | 8×10?? |

图24：水力与热学非均质储层中不同区域单元的质量传递和热传递形状因子随时间的变化

图25和图26分别展示了水力非均质算例中基岩压力的等值线和流线图。这些结果清晰地揭示了水力非均质性对储层水力行为的影响。值得注意的是，由于I区渗透率低、阻碍流体流动并限制压力耗散，因此该区域压力等值线较平坦且压力值较高。相比之下，IV区等值线较稀疏，反映出其高渗透率特性有利于流体向生产井方向流动并更有效地降低压力。需要指出的是，由于该储层尺寸较小且整体渗透率较高，其瞬态持续时间比前节研究的裂隙构型储层更短。此外，从图26的结果可以看出，在较早的时间步，IV区的流线较密集。这反映了IV区的高渗透率创造了冷注入水从注入点到产出区域的优先通道。I区流线较稀疏，表明由于其低渗透率限制了流体渗透。在较晚的时间步，流线呈现出更为分散的分布模式，表明流动状态正在向稳态过渡，使流体能够更有效地渗入低渗透率区域。

图25：地热储层中水力非均质算例在不同时间步的基岩压力等值线

图26：地热储层中水力非均质算例在不同时间步的流线图

此外，图27展示了水力非均质算例的基岩温度等值线。在水力非均质算例中，四个区域具有均匀的热学特性，温度等值线向高渗透率的IV区偏移，反映了储层的水力-热学耦合行为。因此，冷注入水由于其较高的渗透率而更深入地渗入IV区，使流体和热量能够快速向II区运动。相比之下，I区受冷锋影响较小，由于其低渗透率限制了热传递，温度保持较高。最后，图28展示了热-水力非均质算例的温度等值线，其中IV区的裂隙热容量和导热系数降低。可以看出，在水力非均质算例中，渗透率差异控制着冷锋传播的主要路径和速率，而热-水力非均质算例中引入的热学非均质性则进一步加速并锐化了冷锋的传播。

图27：地热储层中水力非均质算例在不同时间步的基岩温度等值线

图28：地热储层中热-水力非均质算例采用所提出的THM双重孔隙模型在不同时间步的基岩温度等值线

相应地，冷锋深入渗透至IV区，这得益于该区域较低的热容量和导热系数。换言之，当冷锋到达IV区时，由于其低热容量，岩石能够被更迅速地冷却。此外，由于导热系数较低，热量难以从周围岩石扩散进入该区域。这一行为强调了在预测热突破和设计地热工程时同时考虑水力特性和热学特性变化的重要性。这些结果表明，所提出的两步THM算法能够捕捉地质非均质性对裂隙地热储层中流动模式和热行为的影响。

**5.4 数值结果的讨论**

数值结果展示了所提出的两步THM双重孔隙框架的关键特征和优势。首先，所提出的质量传递和热传递形状因子是瞬态耦合量，而非固定参数。其次，如公式(17)和(22)所示，质量传递形状因子综合考虑了流体储存、体积变形和热膨胀的影响，而热传递形状因子则考虑了热储存、热传导和热对流过程。结果进一步表明，质量传递形状因子和热传递形状因子本质上是不同的，在一般情况下不应假设二者相等或简单成正比。所提出的算法显著提高了双重孔隙模型的精度，特别是在瞬态条件下。在单块和裂隙土柱算例中，瞬态形状因子得出的预测结果比基于恒定形状因子的传统模型更接近细网格单一孔隙（DNS）解。这种性能的改善源于对形状因子进行依赖于问题的评估，该评估动态反映了当前边界条件、材料属性、基岩块几何形状以及瞬态THM行为。尽管不同的双重孔隙模型趋向于在后期时间收敛到相同的稳态解，但所提出的方法提供了对瞬态演化过程更为精确的描述，从而提高了双重孔隙模型对耦合THM过程的预测能力。

数值研究还展示了所提出框架的通用性和计算效率。由于形状因子以通用形式推导，其适用性不限于方形基岩块；可以方便地扩展到各种块体几何形状、边界条件和材料属性。对裂隙土柱问题的运行时-误差分析表明，与传统的双重孔隙模型相比，所提出的方法与DNS解的吻合度更好，同时计算量远低于细网格单一孔隙/DNS方法。对于大规模裂隙系统，所提出框架的计算优势预计将更加显著，因为在这些系统中直接进行细网格或DNS模拟通常代价过高。通过将精度提升与计算成本降低相结合，所提出的方法为模拟裂隙多孔介质中耦合THM过程提供了一种有效且实用的替代方案。

最后，采用Petrov-Galerkin迎风格式来增强对流主导热传递问题中的数值稳定性。尽管本研究未对标准Galerkin格式和Petrov-Galerkin格式进行直接比较，但该稳定化方法的有效性已在先前的研究中得到了充分验证，包括Cui等人（2018）和Taghvaei与Khoei（2025）的研究。与这些发现一致，本研究中进行的完全耦合THM模拟表现出稳定的温度解，未出现明显的数值振荡。

**6. 结论**

本文提出了一个完全耦合的热-水-力（THM）计算框架，以提高裂隙储层双重孔隙建模的精度。推导了质量传递形状因子的通用形式，考虑了固体变形、流体流动和热传递过程的耦合效应。通过在可渗透基岩块中考虑对流和传导机制，推导了瞬态热传递形状因子的解析表达式。此外，通过对方形和非方形块的分析，验证了这些瞬态质量传递和热传递形状因子的性能和精度。结果表明，热传递形状因子和质量传递形状因子均表现出随时间变化的特性，且不一定相等。同时表明，恒定形状因子在瞬态阶段可能导致显著误差。

在所提出的基于有限元的双重孔隙计算方法中，第一步是选取一个具有与主问题相同裂隙构型和热-水-力属性的单元。在一致的边界条件下对该单元进行分析，并在每个时间步以完全耦合的方式计算相应的热传递和质量传递形状因子。然后将计算所得的形状因子用于主域的双重孔隙分析中。因此，当基于当前域几何和物理条件采用瞬态质量传递和热传递形状因子时，裂隙与基岩之间的传递过程可以以更高的精度进行建模。通过多个代表性算例验证了所提出的两步THM框架的性能和效率，包括两个包含多个具有不同块体构型区域的非均质地热储层。总体而言，结果表明该框架能够捕捉在地质非均质性显著影响流动模式和热行为情况下的储层响应。

**CRediT作者贡献声明**

Mahtab Taghvaei：撰写——原稿，验证，软件，方法学，形式分析，概念构思。Amir R. Khoei：撰写——审阅与编辑，项目管理，指导，验证，概念构思，调查研究，方法学。

Mahtab Taghvaei | Amir R. Khoei

结构抗震工程技术中心，土木工程学院，德黑兰理工大学，11365-9313，德黑兰，伊朗