一种用于同轴孔洞换热器的完全耦合数值方法，其中工作流体为（此处应具体说明工作流体的类型）

《Geoenergy Science and Engineering》：A fully coupled numerical method for coaxial borehole heat exchangers using as a working fluid

　　马格努斯·瓦根
能源技术研究所，Instituttveien 18，Kjeller，N-2027，挪威

摘要
我们提出了一种新的完全耦合和隐式的数值方法，用于求解使用可压缩工作流体（如超临界流体）的同轴钻孔热交换器（CBHE）的压力和温度问题。该方法的目的是展示一种具有良好收敛特性的稳健模拟方法，以用于设计和操作利用可压缩流体的CBHE。模型基于圆柱对称性，岩石中的瞬态热传导采用有限元方法（FEM）进行处理，而井内的质量、动量和能量守恒则采用有限差分方法（FDM）处理。该模型的新特点是所有未知数在规则的FE网格中都有对应的节点。利用FE网格中的节点进行岩石中的FEM和井内的FDM，简化了两种不同数值方法的集成以及系统矩阵的分配。通过与不可压缩流体恒定质量流量的情况的一维解析解和一维数值实现进行基准测试，结果表明基于均方根（RMS）误差的测试结果非常优秀。该模型用于比较在不同质量流量条件下使用液体/超临界流体的热产量。比较结果表明，对于所有流量范围，超临界流体的热产量均高于传统流体。由于其低压缩性，超临界流体在低流量时具有更高的出口温度，而由于其较低的热容量和粘度，在高流量时具有更高的出口温度。

1. 引言
地热能是一种重要的、几乎无排放的化石燃料替代能源，可用于供暖（IEA，2022年），并且可以将用于供暖的电力释放出来用于其他用途。浅层地热能已经成功应用了几十年（Lund和Boyd，2016年；Figueira等人，2024年；Roka等人，2025年）。下一步是进一步开发中深层和深层地热能（Welsch等人，2015年；2016年；Xu等人，2026年），因为随着井深的增加，岩石温度的升高可以带来更多的热量产出。浅层、中深层和深层地热系统之间的界限在深度或温度上并不明确。这里，中深层系统的初始井底温度范围为30°C至80°C（Berger等人，2025年）。浅层和深层地热系统的温度分别低于和高于这一范围。浅层地热能通常基于在150米至300米深度区间内的钻孔中的闭合回路系统，使用水作为工作流体。与开放式地热系统相比，闭合回路系统具有多个优势（Chen和Tomac，2023年）。开放式系统至少需要两口井，井之间有含水层或裂隙岩体。裂隙岩体生成的一个可能问题是水力压裂过程可能会引发地震活动（Schultz等人，2020年）。此外，很难确保压裂过程产生所需的裂隙网络，并且难以量化裂隙网络的水力和热性能（Viswanathan等人，2022年）。开放式系统的所有部分经常遭受矿物沉淀（结垢）问题——包括裂隙网络、井和地面设备（Kamal等人，2018年）。闭合回路系统只需要一口井，不会引发地震，结垢风险极小，并且水力和热性能通常更容易估算。如果用超临界流体替代水，则可以解决结垢问题。然而，水和超临界流体的组合具有腐蚀性，这可能是一个问题。超临界流体是浅层到中深层地热应用中常用的工作流体选择，其中井底温度低于100°C。一段时间以来，人们提出用超临界流体代替水作为工作流体。其中一个原因是超临界流体具有显著的热虹吸效应，可以无需泵就能循环流体。多位作者比较了超临界流体和水的适用性，结果表明两种流体的最佳选择取决于地热应用和操作条件（Hu等人，2020年；Scherer等人，2020年；Ungar等人，2024年；Evangeline和Darwin，2025年）。浅层闭合回路地热井通常基于U型管热交换器（Harris等人，2024年），它由两根平行管组成：一根管用于向下加热流体，另一根绝缘管用于将加热后的流体返回地面。研究表明，同轴钻孔热交换器（CBHE）的性能优于U型管热交换器（Brown等人，2024年）。使用超临界流体的CBHE比使用近似不可压缩流体（如水）的CBHE更复杂。对于不可压缩流体，温度方程与压力方程解耦，压力仅由井内流体柱的重量决定。而对于可压缩工作流体，通常的方法是在迭代循环中依次求解温度和压力（Tonkin等人，2021年；Wangen，2025年）。问题是，在使用可压缩流体的情况下，我们是否可以期望获得关于井内压力和温度以及岩石温度的可靠结果。在这里，我们提出了一种完全耦合和隐式的数值方案，用于求解井压力和温度以及岩石温度。该方法旨在展示一种具有良好收敛特性的稳健数值方法，以用于模拟利用可压缩流体的CBHE。模拟的目的可以是优化设计和操作，例如，通过所需的热产量和/或出口温度来限制井深、管径、材料属性和质量流量。模型假设井周围具有圆柱对称性，这是单井的常见假设。它将笛卡尔坐标系中的完整三维问题简化为圆柱坐标系中的二维问题。与完整的三维解相比，圆柱坐标系下的公式化在计算上更为轻量。假设流体处于准稳态，这对于井内停留时间超过几倍的时段是一个合理的假设（Wangen，2025年）。岩石的瞬态冷却通过热传导实现。井内的压力和温度采用有限差分方法（FDM）离散化，岩石中的瞬态温度通过有限元方法（FEM）离散化。

该模型的新特点是所有未知数在规则的FE网格中都有对应的节点。这意味着岩石温度以及环形空间和内管中的流体压力和温度都由同一网格中的节点表示。然后，FE网格中节点的系统矩阵覆盖了所有未知数，一个涉及系统矩阵的解可以得出所有未知数。无需分配专门的系统矩阵来表示井与岩石之间的相互作用（Diersch等人，2010年）。使用规则网格的计算优势在于系统矩阵是一个带宽有限的带状矩阵。使用规则网格简化了井中节点的管理，例如找到环形空间中与节点的相邻岩石界面节点。

本文的组织结构如下：首先简要回顾了一些常用的和最新的使用可压缩工作流体的BHE方案。在介绍FE网格之前，给出了单相管流的一维控制方程。接下来的部分介绍了岩石中瞬态热传导的FEM以及完全耦合和隐式的压力和温度方案。然后，该方法首先与不可压缩流体的一维解析解进行基准测试，然后再与使用超临界流体的一维数值实现进行基准测试。最后，在相同操作条件下比较了超临界流体和水作为工作流体的情况。

2. 使用超临界流体模拟CBHE的数值方法
在地质热井建模中，已经采用了多种方法。从早期的单相流分析模型到最近的多相现象的计算流体动力学模拟，管道流动的流体力学有着悠久的历史（Chung，2002年；Ishii，2011年；Ferziger和Peri?，2012年；Hammer和Morin，2014年；Lecheler，2022年；Leontidis等人，2023年）。早期的多相流体流动模型包括与周围环境的热交换（Ishii，2011年）。Ishii（2011年）阐述了一相和两相流动的基本流体力学。Tonkin等人（2021年）回顾了瞬态两相流的数学基础及其在地热模拟器中的实现方式。这些模拟器被应用于具有独立注入井和生产井的开放式地热系统中的蒸汽和水两相流动。商业的管道瞬态多相流模拟器（如Olga和LedaFlow）也被应用于地热井（J?rgensen等人，2020年；Fard等人，2006年）。

通过解耦模型的不同部分，可以简化数值方法：（1）可以将压力与温度解耦，（2）可以将环形空间和内管解耦，（3）可以将环形空间与岩石解耦。这些选项可以以不同的方式组合。例如，可以将井内流体流动的完全耦合模型与岩石解耦。Wangen（2025年）测试了四种不同的压力和温度耦合策略。Dendys等人（2015年）提供了用于模拟不同类型地热系统的软件概述。Sun等人（2019年）和Zhao等人（2023年）将环形空间中的温度与内管中的温度解耦。Tonkin等人（2021年）通过数值积分考虑重力和摩擦力的静压力梯度来求解压力。Bai等人（2022年）通过求解井内和岩石中的瞬态温度方程来模拟超临界流体作为工作流体，采用中心差分法。这种方法导致了井和岩石温度的三对角线性方程系统。解耦后的压力是通过数值积分密度依赖的压力梯度得到的。作者采用了小时间步长而不是迭代方法。

几种用于模拟BHE的软件工具基于现有的管道流动模型，例如T2WELL（Vasini等人，2018年；Pan和Oldenburg，2014年）和COMSOL多物理场（Beckers等人，2024年）。使用超临界流体的BHE模拟软件（White等人，2024年；Sandia，2020年；Ransom，1985年）是基于开源多物理场库MOOSE（Permann等人，2020年）构建的。Diersch等人（2010年）使用基于开源库FEFLOW（Diersch，2013年）的3D FE模型完全耦合地求解了与岩石的温度。OpenGeoSys是一个替代的开源库，用于构建地热井的三维模拟工具（Naumov等人，2025年）。OpenGeoSys拥有针对不同类型BHE的模块。当质量流量恒定时，通常使用Ramey的热通量近似方法。Ramey的近似方法将问题简化为一维问题，其中唯一的未知数是井内的流体压力和温度（Yu等人，2019年；Sun等人，2019年；Zhao等人，2023年）。除了短时间外，Ramey的近似方法是准确的（Hagoort，2004年）。Sun等人（2019年在使用超临界流体求解CBHE的质量、动量和能量守恒瞬态方程时假设了Ramey的热通量。环形空间和内管中的压力和温度是在一个主循环中依次解决的，其中收敛标准是对模拟误差的限制。赵等人（2023年）提出了一个使用作为工作流体的模型，并假设了来自岩石的Ramey热流。他们通过使流体沿着内管向下流动并沿着环形空间向上流动来模拟反向模式下的CBHE（环形套管换热器）。作者使用有限差分方法解决了井的瞬态温度方程。稳态压力方程通过数值积分来解决。温度和压力方程在另一个外部循环中解耦，直到温差小于给定的限制。Wangen（2025年）比较了基于Ramey热流的三种CBHE方案：（1）完全耦合的隐式压力和温度，（2）隐式温度和显式压力，以及（3）显式温度和显式压力。第二种方案，即隐式温度和显式压力，产生的结果与完全耦合方案相似。

3. 同轴换热器中流体流动的控制方程
管道中流体流动的基本物理定律是质量守恒、动量守恒和能量守恒（Ishii, 2011）。一维流动的质量守恒方程为
（1）
其中t是时间，x是沿井的距离，ρ是流体密度，v是横截面积上的平均流体速度。类似地，动量守恒方程为
（2）
Darcy-Weisbach方程（Incorpera等人, 2011）给出了由于摩擦而产生的压力降的梯度，μ是摩擦系数，d是水力直径
（3）
而tw是湿润周长。对于内管，水力直径等于管径；对于环形空间，则是外径与内径之差。在层流状态下，摩擦系数为
（4）
而在湍流状态下，我们使用Colebrook摩擦系数（Incorpera等人, 2011）
（5）
其中k是管道的粗糙度。沿井的单位体积流体的重力为
（6）
其中g是重力加速度，h是井的坡度。当深度h随沿井的位置增加而增加时，坡度h为正。一维流体流动的能量守恒表示为
（7）
其中U内是内管的流体能量，Q是能量的源/汇项。它代表从内管到环形空间的热量泄漏以及来自周围岩石的热量输入。Q的单位是J（焦耳）。
假设管道流动处于准稳态，这意味着质量守恒、动量守恒和能量守恒方程中的瞬态项可近似为零。能量方程中唯一的时间依赖项是通过管道壁的热流。例如，由于质量流量变化或输入温度变化引起的瞬态状态只会持续几个停留时间。停留时间是从一小批流体进入环形空间直到它再次出现在内管出口的时间。对于长井来说，停留时间可能是几小时，这与几周、几个月或几年的恒定运行条件相比是非常短的时间。Wangen（2025年）在水作为工作流体的情况下测试了能量方程中的准稳态近似，发现其与完全瞬态解几乎没有差异。当岩石的冷却或加热时间尺度远大于流体的停留时间时，这种近似是合理的。在准稳态流动的情况下，质量守恒、动量守恒和能量守恒方程的离散形式在附录A中推导出来。

4. 数值网格
该模型假设井的中心轴具有圆柱对称性，网格在圆柱坐标系中变为二维。水平轴是从井中心出发的径向距离，记为r，垂直轴是沿井的距离，记为z（见图1）。岩石中元素的径向厚度按几何级数增加。这使得靠近井的地方网格更加细腻，因为那里的温度梯度最大。

图1. 数值网格分为三个部分：岩石、井以及井底部下方的四个空隙单元。沿着网格线1和3的蓝色井节点分别用于表示内管和环形空间的压力，而沿着网格线2和4的蓝色节点分别用于表示内管中的温度。
网格的径向范围假设足够大，以便在模拟过程中井的热扰动能够传播到外部边界。岩石中的热传导是一个缓慢的过程，通常m的径向范围对于深度达3公里的井，在典型热扩散率范围内的岩石中，25年的时间跨度是足够的。在模拟结束时，可以直接检查半径接近r的元素温度是否没有变化。如果m不够，则将m增加到更大的值并不会显著增加计算时间，因为元素的径向长度是按几何级数增加的。
岩石温度在位于元素角处的节点处表示。因此，在井与流体之间的界面、表面以及模型的其他边界处都有节点。每个元素都具有均匀的性质。FEM（有限元）网格中的节点代表岩石中的温度，同时也代表井中的温度和压力。图1显示了网格节点如何布置以同时表示岩石和井。
每个节点有一个未知数。图1中的黑色节点位于岩石中，蓝色节点位于井中。井中的蓝色节点分为四列：第1列和第2列是内管中的压力和温度节点，第3列和第4列分别是环形空间中的压力和温度节点。第五列中的红色节点位于套管/岩石与环形空间中的流体之间的界面处，这些节点由FE（有限元）和FD（有限差分）方法共享。这些界面节点位于环形空间外表面的半径处。这些节点确保了界面两侧的温度连续性。热流在界面处也是连续的，因为从岩石通过界面的热流与进入环形空间另一侧流体的热流相同。必须在界面的两侧施加热流连续性的边界条件。
井节点引入的节点耦合比FE网格的最近邻连接更多。内管和环形空间的温度节点是耦合的，但它们不是最近邻节点。带状矩阵为这些额外的耦合提供了空间。对于 Regular Grid，如果垂直于井的节点数量远少于沿井的节点数量，那么带状矩阵的带宽通常较低。此外，Gauss带状求解器是一种直接、简单且高效的求解器。
在井的底部下方有四个空隙单元（见图1）。空隙单元意味着井底部的节点与井正下方的岩石节点之间没有耦合。相反，要求内管和环形空间中的井底部的压力和温度相同。使用四个空隙单元将井底部与下方的岩石解耦，可以保持网格在元素编号方面的规则性。
岩石中的节点表示实际的网格位置的温度。井中的温度和压力节点位于实际的z位置，而不是实际的径向位置。相反，在FD离散化中直接使用了内管和环形空间的直径。
交错网格是瞬态流体流动的常见数值方法，其中未知的速度在单元接口处求解，而压力和温度在单元中心获得（Gerya, 2019）。这里假设了一维和准稳态流动，直接通过质量流量除以流体密度和管道的横截面积来计算流体速度（见方程（A.1）。尽管没有数学证明这种方案是一致的、稳定的和收敛的，但数值测试表明它的表现良好。例如，测试表明即使在极端流量下，该方案也能很好地工作，除非关闭了流体摩擦。
尽管圆柱对称性不适用于曲线井路径，但它通常是一个很好的近似。即使在井路径有弯曲或井是水平的情况下也是如此。如果温度梯度不大，那么围绕水平井的圆柱计算域内的温差也很小。例如，如果垂直梯度为35°C/km，直径为100米的水平计算域的基部和顶部之间的初始温差为3.5°C。在轴周围10米的直径范围内，向井的热梯度较陡时，温差仅为0.35°C。圆柱对称性在岩石中的流体流动和相关的热对流方面存在局限性。唯一具有圆柱对称性的流动模式是平行于井的流体流动。

5. 岩石中热传导的FE（有限元）方法
岩石中的瞬态三维热传导由温度方程控制
（10）
其中C oslo是恒压下的岩石热容量，ρ oslo是岩石密度，λ oslo是岩石热导率。将（10）与FE基函数相乘，然后在岩石域上积分，得到标准的弱形式
（11）
其中温度用基函数表示为φ，n oslo是井与岩石之间界面的外向单位法向量，q oslo是穿过界面的热流。外垂直边界处的热流为零，井底部中心的岩石也是如此。顶部（z=0）和底部（z=H）的边界具有用户给定的温度。FE网格使用双线性四边形基函数。
围绕中心轴（z=0）的完全对称性给出了del算子
（12）
其中e_z是z方向上的单位向量，e_x是x方向上的单位向量。在圆柱坐标系中对体积的积分包括了这个因子。这个因子可以在对元素的数值积分中包含。这个因子在每个积分点计算并乘以等参映射的雅可比矩阵（Langtangen, 2003）。这样，FE公式在圆柱坐标系中几乎与笛卡尔坐标系中的公式相同。
位于井界面处的节点确保了岩石和环形空间之间的温度连续性。为了确保热流也连续，我们必须施加条件，即界面岩石侧的热流与井侧的热流相同。方程（11）中的表面积分中的热流为
（13）
其中n interfaces是界面节点，n adj是沿井相同z位置的环形空间中的相邻温度节点，k oslo是热传递系数。表达式（13）中的热流与环形空间中的热传递FD方案中的表达式相同（见附录B）。在井侧和岩石侧使用相同的热流表达式确保了界面两侧的热流连续性。完全耦合的压力和温度
完全耦合方法提供了一个系统的解，其中包含了岩石中的温度以及井中的压力和温度。系统矩阵由覆盖岩石的有限元（FE）部分（11）和用于井的有限差分法（FDM）部分组成。附录A提供了稳态管道流动守恒定律的离散化方法，附录B包含了热传递系数，附录C给出了井压和井温的FDM方程，附录D则提供了井压和井温的矩阵元素、热传递系数以及FDM方程。在一般情况下，流体属性（如密度和粘度）取决于流体压力和温度，这使得系统变得非线性。为了求解耦合系统的压力和温度，采用了连续替换法（也称为皮卡德迭代法）。迭代循环开始时使用上一个时间步长的压力和温度来计算流体属性，在后续迭代中则使用上一次迭代的解来计算流体属性。这种方法的收敛性非常好，表现为连续两次压力和温度解之间的最大差异非常小。即使在初始点远离实际解的情况下，收敛速度也很快，通常每次迭代都能减少一个数量级。经过3到5次迭代后，最大温差通常会小于C。

连续替换法是一种适用于弱非线性问题的良好迭代方法。另一种方法是牛顿-拉夫森法，它在接近解时的收敛速度更快。然而，牛顿-拉夫森法在计算上通常要求更高，因此可能收敛速度较慢（Paniconi和Putti，1994年）。牛顿-拉夫森法的另一个缺点是流体属性（如压缩性和热膨胀性）是通过压力和温度的数值导数来计算的。这里的压缩性和热膨胀性是基于线性插值的密度表得出的，而线性插值会导致密度的二阶导数为零。

7. 使用不可压缩流体的基准测试
数值解需要通过与精确解进行比较来验证。不幸的是，在一般情况下，对于可压缩流体，似乎没有精确解。不过，对于流体属性恒定的情况，存在解析解。Kabir等人（1996年）提出了一个在流体属性恒定、流量恒定，并假设使用Ramey近似（Ramey，1962年）来计算井内热流量的情况下的CBHE（井-岩-热交换器）的解。Kabir等人的解后来被扩展到包含水平段的CBHE的情况（Sharma等人，2020年）。Wangen（2024年）提出了一种半解析的多段解方法。该方法对每个段与段之间的两个参数进行了数值求解。

图2(a)–(c)分别展示了本文提出的数值方法与Kabir等人（1996年）的解之间的比较，对应的质量流量分别为0.1 kg/s、1 kg/s和10 kg/s。表1列出了该案例中的参数。在1天、10天、100天和1000天时比较了这些解。结果显示，在超过10天的时间跨度内，两者非常吻合。表2列出了图2中温度解之间的偏差，以均方根误差（RMS error）表示。Ramey近似在一天时间尺度上不够准确，但随着时间的延长，其准确性逐渐提高。这也是表2中偏差逐渐减少的原因之一。尽管图2中的测试是针对不可压缩流体的，但它验证了FDM和FEM在流体与岩石界面处的耦合是准确的。同时，这也证明了在圆柱坐标系中的数值瞬态FE解的准确性。

下载：下载高分辨率图片（590KB）
下载：下载全尺寸图片

图2. 三种质量流量（a）0.1 kg/s、（b）1 kg/s 和 (c) 10 kg/s 下的数值解与解析解的比较。环形空间中的输入温度为25摄氏度。

表1. 不可压缩流体基准测试案例中的参数。井内沿程的温度差是指井表面与底部之间的初始温差。

- r方向上的单元数 [-] 49
- s方向上的单元数 [-] 100
- 周围岩石半径 [m]
- 质量流量 [kg/s]
- 进口温度 [C]
- 表面温度 [C]
- 井长度 [m]
- 井内沿程温度差 [C]
- 内管 radial 厚度 [m]
- 保温层 radial 厚度 [m]
- 环形空间 radial 厚度 [m]
- 流体热容量 ()
- 流体粘度 [Pas]
- 流体热导率 ()
- 岩石热导率 () 2.5
- 保温层热导率 () 0.01

表2. 图2中曲线之间的偏差，以均方根误差（RMS error）表示。单位为摄氏度。

- 空单元
- 流量 0.1 kg/s
- 1天 0.13
- 10天 0.12
- 20天 0.15
- 70天 0.10
- 80天 0.22

- 10天 0.125
- 11天 0.116
- 0.084
- 100天 0.116
- 0.345

图2中的三种质量流量展示了CBHE运行的三种不同状态。图2(a)中的环形空间温度接近初始岩石温度，说明流量足够低，使得环形空间中的流体与岩石达到热平衡。内管中的温度向表面逐渐升高，表明低流量会导致内管向环形空间的热量泄漏。环形空间中的温度在接近表面的地方也呈曲线状，这是因为进口温度高于井上部的岩石温度。因此，流体在加热上部岩石之前，先加热了这部分岩石。图2(c)显示，环形空间中的流体从周围岩石中吸收热量的时间较短，因此其温度明显低于井底部的初始岩石温度。当流体在内管中返回表面时，几乎呈直线，这表明其在内管中停留的时间不足以导致显著的热量泄漏。图2(b)展示了中等流量下的情况，在运行1000天后，井底温度仍保持在50摄氏度以上。

8. 使用可压缩流体的基准测试
为了测试可压缩流体情况下的数值方法，采用了另一种数值解。同样，该测试解基于Ramey的热流量近似，允许使用一维FDM方法来求解CBHE中的压力和温度（Wangen，2025年）。流体属性（如密度、热容量、热导率和粘度）从二维表格中查找。数据来自NIST（2024年），并在数据点之间进行线性插值。两种数值方法（一维参考解和完全耦合解）都使用了相同的流体属性表格。

图3显示了在恒定质量流量下，垂直井的井压和温度随时间变化的情况，时间点分别为1天、10天、100天、1000天和10000天。在这种情况下，输入压力为10 MPa，以确保流体处于液态或超临界状态。液体在C温度下注入，并在m到m的深度区间内变为超临界状态，此时其温度超过了临界温度C。即使是在短短的一天时间内，数值结果也非常吻合。这种优异的吻合再次证明了在圆柱坐标系中实现FE公式以及井与岩石界面之间的热传递是正确的。两种不同的CBHE数值解之间的良好吻合表明数值解是准确的。表3列出了图3中压力和温度解之间的偏差，以均方根误差（RMS error）表示。

图3. 基于Ramey近似热流的数值解（a）与压力（b）和温度（b）的比较。环形空间中的流体输入压力为10 MPa，输入温度为C。

图3(a)展示了热虹吸效应，其中输出压力大于输入压力。这种情况下的压力差为 MPa。热虹吸效应是由于环形空间和内管中流体的密度差异造成的。对于可热膨胀的流体（如水），内管的平均密度显著低于环形空间，因为内管的温度较高。热虹吸效应产生的功率估计为 kW，其中 为出口和入口之间的平均密度。图3(b)显示了流体压力和温度之间的耦合关系。流体在向上流动过程中，随着压力降低，温度也随之下降。由于能量守恒（A.11），压力下降会导致温度下降。在没有热输入的情况下，压力变化会产生温度变化（14），其中温度 以开尔文为单位。使用数据C、1/K和 ，可以计算出m长度内温度差为C。这一项解释了从井底到表面温度下降的大部分原因。

图4显示了10,000时间跨度内输出压力、输出温度和热量产生的比较。图4(a)中输出压力的差异小于7 kPa，对于输出压力为10 MPa的情况来说这是一个很好的吻合。图4(b)中输出温度的差异小于0.15摄氏度，对于输出温度超过30摄氏度的情况也是如此。图4(c)中热量产生的差异小于0.5 kW，对应于大约100 kW的热量产生。

在迭代过程中，使用上一个时间步长的温度解作为当前时间步长的猜测值。初始猜测值与实际温度解之间的差异在第一个时间步长最大，随后连续两次解之间的差异逐渐减小。图5展示了前五个和最后五个时间步长的收敛特性。收敛性通过连续两次迭代解之间的最大温度和压力差来衡量。图5中的收敛标准为C。收敛性向零趋近表明每次迭代都在改进解的质量。虽然标准C相当严格，但它展示了收敛如何趋近于零差异。收敛标准可以放宽至少一个数量级，而不会显著影响结果的质量。

图5. 环形空间和内管中温度节点所需的迭代次数。图表显示了（a）初始差异最大的前五个时间步长，以及（b）最后五个时间步长。最后五个迭代模式有所重叠。还应提到的是，即使对于取心井，井周围岩石的类型和属性也存在不确定性。数值实验表明，在3公里深的井中，将sc-CO2的热导率从2提高至2.2（提高10%），5年后输出温度会升高超过0.5摄氏度，25年后则升高超过1摄氏度。在相同的垂直井中，对于不同的质量流量条件，对两种流体进行了对比。表4列出了相关参数：入口压力为10兆帕（MPa），该压力足够高，使得流体可以处于液相或超临界相。液态sc-CO2在某个温度下被注入井中，当其在环空中被加热时，一旦温度超过其临界温度，就会变成超临界状态。在本节中，“sc-CO2”既可表示液态也可表示超临界状态，因为这两种状态都可能同时存在于环空和内管中。

表4. 流体参数。井中的温差是指井表面与井底之间的初始温差。
- 参数名称 | 值 |
|-----------------|-----------------|
| 单元数（径向方向） | 49 |
| 单元数（轴向方向） | 100 |
| 周围岩石半径 | [m] |
| 入口压力 | 10 MPa |
| 入口温度 | [C] |
| 表面温度 | [C] |
| 井内温差 | [C] |
| 井长 | [m] |
| 内管径向厚度 | [m] |
| 保温层径向厚度 | [m] |
| 环空径向厚度 | [m] |
| 岩石热导率 | [W/(m·K) |
| 保温层热导率 | [W/(m·K) |

内管、井壁和环空的径向厚度是典型的数值，适用于2公里深的CBHE系统（参考文献：Pan et al., 2020; Gruben et al., 2021; Chen and Tomac, 2023; Lund, 2024）。深度超过2公里的井可能需要多层套管（Chen and Tomac, 2023）。为了将输入温度保持在10摄氏度左右，需要使用热泵；不过这些热泵并未包含在模拟过程中。由于输入温度接近地表温度，CBHE系统实际上已经从岩石表面开始吸收热量。

图6(a)显示了在不同质量流量下，两种流体的输出温度以及出口与入口之间的压力差。输出数据分别在10天、100天、1000天和10,000天时进行了记录。从图中可以看出，当质量流量较低时，内管中的流体几乎与周围岩石达到了热平衡；随着流量增加，流体通过内管时热量损失减少，输出温度随之升高，直至达到最大值；而当流量较高时，流体在环空中被加热的时间变短，输出温度则下降。

图6(b)比较了不同质量流量下两种流体的输出温度和压力差。结果表明，在质量流量较低的情况下，某种流体具有更高的输出温度；而在质量流量较高时，另一种流体表现出最高的输出温度。图7展示了两种流体在输出温度最高时的井内温度变化情况。较低的质量流量允许流体在下降过程中有更多时间被加热，因此两种流体在井底的最高温度接近初始温度。然而，在上升过程中，由于压力减小，某种流体的温度会下降。

图8显示，在整个质量流量范围内，某种流体的热输出量均高于另一种流体。不过，内管向环空的热量泄漏与输出温度的变化趋势一致——热量泄漏量与环空与内管之间的温差成正比。因此，输出端与进口端之间的最大温差导致了最大的热量泄漏。

在10兆帕的压力下，某种流体的密度和粘度等物理性质与大气条件下的性质差异不大，这使得该模型结果适用于接近大气压力的入口压力情况。

**结论**：我们提出了一个基于完全耦合和隐式方法的数值模型，用于研究使用可压缩工作流体（如sc-CO2）的CBHE系统。由于模型采用隐式求解方式，因此即使在较大的时间步长下也能保持稳定的数值解；同时，由于模型实现了完全耦合，因此具有出色的收敛性能。该模型基于圆柱对称性，采用二维规则网格进行计算，其中x轴表示从井中心的径向位置，y轴表示沿井方向的距离。假设工作流体的流动处于准静态状态，这一假设在停留时间较长的情况下是合理的。井内的压力和温度通过有限差分法（FDM）求解，而岩石的瞬态温度则通过有限元法（FEM）求解。模型的一个创新之处在于使用有限元网格来表示所有未知量（即井内的压力、温度和岩石温度）。我们采用了带状矩阵来处理模型中的耦合力，这种求解器简单、直接且高效。该模型通过与不可压缩流体和恒定质量流量情况下的解析解进行了对比验证，结果显示两者吻合度很高，证明了圆柱坐标系下的有限元方法和FDM与FEM之间的井/岩石耦合效果良好。此外，还通过与使用Ramey瞬态热通量的完全耦合一维数值模型进行了对比，结果同样显示两种方法的性能非常接近。在相同的操作条件和质量流量范围内（0.1 kg/s至15 kg/s），对2.5公里深的井进行了模拟，结果显示某种流体在所有情况下产生的热量都更多。尤其是在高质量流量下，该流体的输出温度也更高。