焦开拓 孙成珍 白博峰 杨瑞涛

1. 动力工程多相流国家重点实验室·西安交通大学 2.陕西延长石油国际勘探开发工程有限公司

0 引言

国家能源局颁发的《关于促进地热能开发利用的若干意见》[1]指出，到2025年全国地热能供暖（制冷）面积比2020年增加50%，地热能发电装机容量比2020年翻一番，可见中深层地热能是我国重点发展的可再生能源。目前中深层地热能的开采方式主要有水热型和地埋管型，水热型要求有地热流体，对回灌技术要求高，且需面临地下水化学和生物化学方面的风险[2]。当地质条件和回灌技术无法满足时，“只取热不取水”的闭式中深层地埋管型就成为一种替代技术[3-4]。中深层地埋管深度一般在1 000 m以上，常见类型有同轴套管型、U型和双U管型，在波兰[5]、美国[6]、德国[7]等国家都进行了先导试验项目。

地下水渗流使得地层中热传递方式由纯导热转变为导热加渗流换热，有利于将上游远处地热能传递至埋管周围[8-9]，而埋管周围低温区在流体带动下向下游迁移形成低温热羽流区域[10-11]，埋管周围地层温度升高。地下水渗流加强了埋管与地层之间的换热，因此会对中深层地埋管的取热产生明显影响。Kang等[12]计算结果表明，在1.16×10－7m/s地下水渗流速度下可使7天运行时长中埋管沿线热流提升10.08%，也使得取热总功率更快趋于稳定。Jia等[12]对2 505 m水平对接埋管进行研究，发现240 m厚含水层可使得埋管取热总功率提升近5.73 kW。Angelotti等[11]计算结果表明，60天中的埋管取热量随着渗流速度的增加而增加，但是增长斜率会随着渗流速度的增加而慢慢降低。虽然前人研究涉及了地下水渗流对中深层地埋管取热性能的影响，但是对跨度几百天的长时间运行或者几十年的多周期运行情况下的取热性能和取热可持续性影响规律尚待探讨。

中深层地层相较于浅层地层需考虑到更多的地质因素，例如地温梯度[13]、地层物性波动[14]、地层表面温度[15]等。上述几种因素会导致埋管沿线热流随着深度和时间变化明显，进而也需考虑地层和埋管内部流体之间的流固耦合传热。地层温度响应的两种计算方法为数值方法和解析方法。数值方法采用传统离散方法，可精细描述上述的复杂地质因素，但是计算耗时，不同埋管计算时需重新建模，通用性较差[16-17]。较之而言，解析方法计算快捷，在长期取热中具有较好的准确度[18]。传统解析模型通过线热源理论[19]和圆柱热源理论[20]推导得出，近期也有学者发展了混合模型[21]。需要指出的是，解析模型多用于浅层地层，无法直接描述上述复杂地质因素，计算中深层地层需结合其他传热理论进行推导和改进[22]。

本文旨在定量探究不同地下水渗流速度对中深层地埋管长时间、多周期取热性能的影响，进而基于线热源理论引入了自主开发的中深层埋管与地层双向耦合传热解析模型。通过一原位试验井参数，研究了埋管在单个取热期、多周期取热两个情况下的取热性能，并通过变化率、降低程度等自定义参数定量评估了地下水渗流的影响，总结了埋管取热性能随地下水渗流速度的变化规律，以期助推我国地热能开发。

1 物理模型设计

中深层地埋管的传热过程分为两部分：地层渗流传热和井筒内传热。基本建模策略为先单独求解此两部分，之后在井孔外壁处将此两部分耦合。自主开发的埋管与地层双向耦合传热解析模型[23-24]为：在地层渗流传热中应用了改进的解析模型，综合考虑了地下水渗流、分层热物性、地温梯度、地表温度等影响，并且以埋管沿线取热功率随时间和深度变化为前提计算地温扰动；在井筒内传热中，研究了典型的同轴套管式结构，其从内到外依次为内管、内管壁、外环、外管壁、填充材料（通常为水泥），循环水从外环注入，由外管壁加热后从内管采出，此种注入方式相较于从内管注入可有效降低热损失[13]。物理模型如图1所示。图1中的ρi表示第i层地层密度，kg/m；λi表示第i层地层热导率，W/(m·℃)；Cρi表示第i层地层的比热容，J/(kg·℃ )；Ui表示第i层地层达西速度，m/s；qs(i)表示第i层线热源，W/m；zi表示空间坐标，m；Tin和Tout分别表示流进和流出管段的水温，℃；下标i表示地层，i=1, 2, 3, …, N。

图1 物理模型示意图

1.1 地层渗流传热特征

为了简化地层渗流传热计算，引入以下假设：①各地层沿深度方向排布，在水平方向无限长，埋管直径相较而言很小，因此埋管取热在地层中可近似为一根垂向的线热源；②地下水仅沿水平方向流动，且同一地层中的渗流速度相同，由于渗流速度较小，多孔介质中为局部热平衡状态；③地表温度（T0）和埋管底部向下100 m处无扰动层的温度（TH+100）不受埋管取热影响且不随时间变化。定义地层中(x,y,z)位置处τ时刻的过余温度Θ为该处温度与埋管底部向下100 m处无扰动层的温差，即Θ由两部分组成，长度为H的埋管取热引起的温度扰动 和上下恒温层的影响Θs(z)。

Θb与地层初始和边界温度无关，当线热源和渗流速度沿深度方向为常数且均不随时间变化、各地层热物性相同时，可应用移动有限长线热源模型求解[23]，表示为：

式中Θb表示过余温度，℃；qs表示线热源，W/m；H表示埋管长度，m；x、y、z表示空间坐标，m；τ表示时间，s；g表示单位长度线热源的积分函数。

式（1）中引入了从－H至0积分的虚拟热源的影响，使得地表的温度扰动为零；定义则函数g可表示为：

式中λ表示地层热导率，W/(m·℃)；a表示地层热扩散系数，m2/s；v表示等效渗流速度，m·s－1；U表示达西速度，m/s；ρf、ρm分别表示流体和地层的密度，kg/m；Cρf、Cρm分别表示流体和地层的比热容，J/(kg·℃)。

实际运行中，深层地层温度更高，与管内流体温差更大，因而qs值更大。为在模型中考虑qs随深度的变化，需将各地层段不同值的线热源以及虚拟线热源叠加。此外，实际运行中qs随时间也是变化的：在取热期，中深层地埋管处于运行状态，qs随运行时间增长指数降低；在恢复期，埋管停止运行，qs维持为0；在恢复期转变至下一取热期时，qs有一阶跃增长，开始下一周期的取热。为在模型中考虑qs随时间的变化，引入Duhamel定理，即将qs随时间的连续变化近似等价为离散的阶梯变化，并记录相邻阶梯间的变化量，进一步将式（1）和（2）修正为[23-24]：

其中

式中N表示地层层数；i表示第i层地层的参数；ω表示总的离散时间步数；上标R、V分别表示真实热源和虚拟热源；下标comp.表示复合地层的热物性，用于表征地层热物性沿深度的变化。

Θs(z)反映了地表和无扰动层之间地温梯度的影响，其求取可被视为无限大复合平板中的一维稳态非齐次导热问题，相关的控制方程和求解公式参见本文参考文献[24]。总的过余温度由Θb和Θs相加得到：

1.2 井筒内传热特征

在井筒传热中，由于沿线水温发展至稳态的时间远小于地层中离散的一个时间步，可采用稳态传热模型来计算井内传热。内管和外环内流体的传热如式（7）所示[23]：

式中下标1、2分别表示内管和外环中的水；Tf表示管段内平均水温，℃；Tw表示管壁温度，℃；qv表示体积流量，m3/s；Tin、Tout分别表示流进和流出管段的水温，℃；R12表示内管中循环水和外环中循环水之间的换热热阻，(m·℃)/W；R2g表示外环中循环水与井孔外壁之间的换热热阻，(m·℃)/W。

为耦合地层渗流传热和井筒内传热两部分，井孔外壁处需满足热流相等、温度连续条件，记为：

式中rb表示井孔外壁半径，m。

式（9）计算的是填充材料与地层的交界面处地层温度的平均值，选取渗流上游、下游和两个垂直方向共四个交点。联立求解各层对应的式（7）～（9）所构成的方程组，可得各层的六个未知量之后将回代至式（6）可得任意时间步和任意位置处的地层温度。

1.3 原位试验参数

基于一真实试验井进行模型验证和案例研究，该井为同轴套管式，深2 500 m，位于渭河盆地中南部的地热田。原位试验井地埋管参数如表1所示，利用该原位试验井的测完井曲线可得地层热物性参数如热导率、密度和比热容随深度的分布[23-24]，如图2所示。根据地质调查结果，该地地温梯度为33 ℃/km，在850 ～1 750 m处存在含水层，达西速度为9.59×10－8m/s。该试验系统不仅记录埋管进出口水温，还采用分布式光纤系统对沿线外管壁温度数据进行采集，进而可利用外管壁温度随时间变化的试验数据对模型的非稳态传热计算进行验证，模型验证结果详见笔者之前成果[23]，此处不再赘述。

表1 原位试验地埋管系统参数表

图2 地层热物性参数随深度的变化图

2 地下水渗流对取热性能影响

根据上述原位试验井参数，分析了地下水渗流对中深层地埋管在单个取热期、多个周期两个情况下的取热性能。根据前人研究[25-26]，埋管热流密度随着深度的增加而增加，深处地层渗流传热情况对整根埋管的取热性能影响更大，因此以深层地层作为研究对象更具代表性。在本节分析中，假设含水层仅存在于埋深2 000～2 500 m处，达西速度（U）介于0～8×10－7m/s，埋管入口工况维持在Tin=5 ℃、Qin=25 m3/h。在模型计算方面，每层地层考虑的厚度为50 m，单个取热期持续120天，离散时间步长6小时，而多周期取热持续50年，一个周期包括取热期120天、恢复期245天，离散时间步长120小时。

2.1 单个取热期

在含水层中的各个渗流速度下，120天内埋管的出口温度和取热功率如图3所示。由图3-a可知，取热初始时刻埋管的出口温度和取热功率达到最大值，随着运行时间累积而衰减。伴随着达西速度(U)从0 m/s升高至8×10－7m/s，取热期内同一时刻的埋管出口温度和取热功率均有所增加，第120天时埋管的出口温度从16.1 ℃增加至17.8 ℃，取热功率从323.5 kW增加至371.6 kW。图3-b中还显示了出口温度增幅和取热功率增幅，代表某达西速度下的值与同一时刻U=0 m/s时的值之差。可以看出在运行初始阶段，地下水渗流引起的出口温度增幅和取热功率增幅较小，各个达西速度下均近似为0 m/s，说明运行初期埋管的取热性能不如地下水渗流速度敏感。随着运行时间的累积，出口温度增幅和取热功率增幅的值渐渐增大且增长斜率先大后小，各个达西速度下增幅值之间的差异也越来越大。总体来讲，地下水渗流速度越大，埋管取热性能提高越明显，但是地下水渗流对于取热性能的促进作用主要体现在取热后期。

图3 各达西速度下埋管出口温度和取热功率随时间的变化图

为了便于定量分析埋管的取热性能，定义如下两个参数[23]：

由上式定义可知，变化率是指当前时步埋管取热总功率相对于上一时步埋管取热总功率的变化量，而降低程度指的是当前时步埋管取热总功率相对于第一个时步埋管取热总功率（此时功率最大）的降低量。此外，将式中求和符号去掉后，可分别表示第i层地层的局部变化率和局部降低程度。

图4反映了U=0 m/s时埋管取热性能降低程度和变化率随时间的变化。可看出降低程度随时间增长，最大值为54.7%，而变化率随时间降低，最小值为0.10%。两个参数均在取热初始阶段变化明显。本文定义当变化率小于0.5%时为稳定取热阶段，否则为瞬态取热阶段。因此，图4中达到稳定取热的时间为27.7天。由图4中饼图可知，在第27.7至120天的稳态阶段的取热量占总取热量的71.87%，说明稳态阶段是衡量埋管取热性能的关键。

图4 取热性能降低程度和变化率随时间的变化图

图5反映了第120天时埋管取热性能降低程度、达到稳定取热阶段的时间和稳定取热占比随达西速度的变化，也考虑了埋管沿线均为无含水层时的情况（全线地层不含水）。在有含水层但内部流动是静止的情况下（即U=0 m/s），其120天时降低程度和稳态时间均较无含水层时更大，分别高1.1%和0.6天，而稳定取热占比相较无含水层时低0.7%。原因为含水层内孔隙水为饱和状态，孔隙水热扩散系数与固体骨架相比较小，进而导致含水层中的地层热扩散系数比无含水层中的小；较小的热扩散系数意味着埋管周围的低温响应区更难扩散出去，致使周围土壤温度随着时间累积降低得更明显，且更难达到稳定取热状态。另一方面，随着达西速度从0 m/s升至8×10－7m/s，120天时降低程度和稳态时间降低（分别从54.7 %降至48.2 %，从27.7天降至21.5天），稳定取热占比增大（从71.9 %升至78.6 %）。这说明地下水渗流增强有助于取热稳定和减弱取热性能随时间的降低，原因为地下水渗流促进了埋管取热周围低温区的耗散，提高了埋管周围土壤温度。综合来看，仅有地下水的存在并不能使得埋管快速达到取热稳定，而是在有明显的地下水流动时才会。

图5 降低程度、稳态时间和稳定取热占比随达西速度的变化图

为研究埋管沿线局部取热性能，图6展示了无含水层以及各个渗流速度下局部热流和局部降低程度随深度的变化。由图6-a可知，由于地温梯度的影响，随着深度增加，0～2 000 m处埋管局部热流近似线性增长，其中也存在着由分层热物性所致的波动。在遇到含水层时，由于孔隙水增加了地层热导率，局部热流相较于无含水层有一阶跃增长。在第1天取热时这一阶跃增长在各个达西速度下相同，均约为103 W/m。然而，在已达到稳定取热的第120天时，阶跃增长随着达西速度的增大而增大，例如在U=0 m/s时阶跃增长仅为35.2 W/m，而在U=6×10－7m/s时阶跃增长为81.5 W/m。由图6-b可知，0～2 000 m处埋管的局部降低程度随着深度的增加而减小，说明深处地层的取热可持续性优于浅层地层。在2 000 m处遇到含水层时，不同达西速度的局部降低程度呈现不同的变化：当U为0～2×10－7m/s较小值时，局部降低程度有一阶跃增长，其值比无含水层中的大，说明此状况下地层的局部可持续性低于无含水层。当U＞3×10－7m/s时，局部降低程度有一阶跃减小，说明此状况下地层的局部可持续性高于无含水层，且越大的达西速度具有越明显的阶跃减小，从而在U=8×10－7m/s时局部降低程度仅为38.3%。总体而言，虽然地下水渗流可以增长取热稳定期埋管的局部热流，但是仅在地下水U≥3×10－7m/s时才能增强埋管的局部可持续性。

图6 第120天时局部热流和局部降低程度随深度的变化图

2.2 多周期取热

在多周期取热中，变化率和降低程度的定义变为如下[23]：

由上式定义可知，变化率是指第N年第m时步埋管取热功率相对于第（N－1）年第m时步的变化量，而降低程度为第N年第m时步埋管取热功率相对于第1年第m时步的降低量。同样的，将（1）式中求和符号去掉后，可表示第i层的局部变化率和局部降低程度。

图7-a显示了U=0 m/s时前5年周期运行下的出口温度和取热功率。可以看出取热期开始时刻和结束时刻的出口温度和取热功率均随着运行年限的增加而衰减，取热期开始时刻可视为每个取热季的瞬态阶段，而取热期结束时刻可视为取热季的稳定阶段。图7-b展示了U=0 m/s时两个时刻的降低程度和变化率随运行年限的变化情况。在变化率方面，运行前10年取热期开始时刻的变化率高于取热期结束时刻，之后两者基本重合，说明10年之后地层中取热期取走的热量和恢复期恢复的热量达到基本平衡。而对于降低程度来说，取热期开始时刻和结束时刻之间的差距在前15年逐渐增大，之后两者差值基本不变，约为2.2%。这说明取热季瞬态阶段取热性能随运行年限的变化比取热季稳定阶段更明显，降低幅度更高。

图7 多周期运行下取热期初始时刻和结束时刻变化情况图

图8展示了各个达西速度下多周期运行第50年时埋管取热性能的降低程度，其中也考虑了全线为无含水层的情况。与单个取热期中的分析类似，无含水层时的降低程度在取热期开始时刻和结束时刻均比达西速度为0 m/s时小。随着U从0 m/s升高至8×10－7m/s，第50年时取热期开始时刻的降低程度和取热期结束时刻的降低程度均开始减小，分别从17.7%降至9.6%，从15.5%降至7.7%，降低速度均为由快变慢。取热期开始时刻的降低程度始终大于取热期结束时刻，两者的差值呈先减小后增大的趋势，当U=4×10－7m/s时达到最小值1.4%。总体而言，愈大的地下水渗流有助于减弱埋管取热性能随运行年限的降低，提高埋管的可持续运行能力，但是内部静止流动的含水层会起到相反效果。

图8 多周期运行第50年时取热期开始时刻和结束时刻的降低程度图

图9展示了多周期运行第50年时各渗流速度下埋管沿线的局部降低程度。可以看出无论是取热期初始时刻还是结束时刻，随着深度的增加局部降低程度总体上在减小，但在遇到含水层时不同达西速度对应着不同的变化。当U=0 m/s时，局部降低程度相对于无含水层有一微弱的阶跃增长，说明此状况下地层随运行年限增长的局部可持续性低于无含水层。当U≥1×10－7m/s时，在遇到含水层时有一阶跃减小，且减小幅度随着达西速度的增加而增加，以至于U≥4×10－7m/s时局部降低程度近似减小为0，这意味着埋管50年每个取热期的取热能力都维持在第1年的状况，有利于埋管的取热可持续性。值得指出的是，在U=0 m/s和无含水层状况时，埋管端部和底部的局部降低程度有一阶跃减小，这是由于恒温边界T0和TH+100的影响减弱了周围地层的地温波动，增强了附近地层的局部取热可持续性。

图9 多周期运行第50年时埋管沿线局部降低程度图

3 结论

1）对于120天的单个取热期，在运行初始时刻，地下水渗流对埋管取热性能影响较小，而随着运行时间累积，其引起的出口温度增幅和取热增幅开始显现。地下水渗流引起的增幅作用随着达西速度的增大而增加，当U从0 m/s升高至8×10－7m/s时，第120天时出口温度和取热功率分别提高了1.7 ℃和48.1 kW。然而，仅有含水层而没有地下水流动会使得稳态时间和降低程度高于无含水层，原因为孔隙水减小了地层整体的热扩散系数。

2）单个取热期下埋管的局部热流随着深度的增加而近似线性增长，且在遇到含水层时有一阶跃增长。在埋管取热的初始时刻，该阶跃增长值在各个达西速度下维持不变，而在取热稳定阶段，该阶跃增长值随着达西速度的增加而增长。当U≥3×10－7m/s时，埋管的局部降低程度低于无含水层情况，说明此时地下水渗流增强了埋管的局部取热可持续性。

3）对于运行50年的多周期取热，埋管在第10年时取热性能达到准平衡状态。取热期开始时刻（瞬态取热阶段）的取热性能随运行年限增长的变化率和降低幅度均要高于取热期结束时刻（稳定取热阶段）的情况。更大的达西速度有助于减小埋管取热性能随运行年限的降低，随着U从0 m/s增长至8×10－7m/s，取热期开始时刻和结束时刻的降低程度均减小，第50年时范围分别为9.6%～17.1%和7.7%～15.5%。

4）在多周期取热中，埋管沿线的局部降低程度总体上随着深度的增加而减小，说明深层地层的取热可持续性优于浅层地层。在遇到含水层时，若内部孔隙水不流动，埋管局部降低程度有一阶跃增长，不利于埋管在运行年限上的局部取热可持续性，而内部孔隙水有流动时对应着阶跃减小。此外，当U＞4×10－7m/s时，埋管每个周期的取热性能近似维持不变。