中深层地埋管群的传热分析

方肇洪贾林瑞张方方崔萍

(山东建筑大学 热能工程学院，山东 济南250101)

0 引言

地埋管地源热泵供热空调系统是以大地为冷、热源，使作为中间介质的热载体在埋设在岩土中的封闭环路中循环流动，与大地进行热量交换，并进而通过热泵实现对建筑物的供热和供冷。 地源热泵技术对建筑节能、CO2减排和防治大气污染都起到了重要作用。 地源热泵供热空调系统由地埋管换热器、热泵主机以及建筑物内的末端系统组成。 地埋管换热器是地埋管地源热泵空调系统区别于其他供热空调系统的最大特点，其中的传热是管内流体与周围岩土之间的换热，可以看作是一种蓄热式换热器[1]。 中深层地埋管技术是近年来在国内外地源热泵领域出现的一种新技术，埋管深度达到数千米[2]。 在地源热泵和蓄热技术的应用中，常常由若干地埋管组成地埋管群。 中深层地埋管换热器的布置方式必须因地制宜，可以是多种多样的，空间尺度为数千米，涉及多个地层，因此其几何条件很复杂。地埋管换热器中的传热过程又是一个复杂的、非稳态的传热过程，所涉及的时间尺度很长(至少为数月至数年)，特别是由多个地埋管换热器组成的地埋管群中的传热问题都必须按三维问题来考虑。 因此，中深层地埋管换热器传热的研究是地源热泵空调系统的技术难点，同时也是该项技术研究的核心和应用的基础。

地埋管换热器的传热分析方法大体上可以分为两大类，即解析解方法和数值解方法，其各有长处和局限性[3-4]。 解析解方法的基本出发点是研究无限大或半无限大介质中单个钻孔埋管的传热规律，此时地埋管换热器简化为线热源或圆柱面热源，并得到其在恒定热流作用下温度响应的解析解。 20 世纪80 年代，ESKILSON[5]提出了一种叠加原理的新思路，基于描述无限大介质中单个钻孔在恒定热流加热条件下的温度响应，再利用叠加原理得到多个钻孔组成的地埋管换热器在变化负荷作用下的实际温度响应。 这种方法的物理概念清晰、计算简捷，考虑到了地埋管换热器的复杂几何配置和负荷随时间的变化，已经广泛应用于实际的工程设计计算和能耗分析[6-8]。 这种解析解方法的一个不足是需要假设地下岩土为均匀介质并有均匀的初始温度。 中深层地埋管换热器深入地下数千米，地温梯度成为影响传热过程的主要因素，因此继续假定地下岩土整体为均匀介质并忽略地温梯度明显是不适合的。 此外，在地埋管换热器传热的解析解分析方法中，需要把地埋管换热器的传热分成两个区域分别进行处理，在传统的解析解方法中需要在两个区域的界面，也就是钻孔壁上，通过简化假定将两者联系起来。通常的做法是对于钻孔内的问题假设钻孔壁面温度均匀，对于钻孔外的问题假设钻孔壁上热流均匀。这样的简化假定明显是不合理的，但是由于地埋管换热器传热计算复杂而耗时，为了能够采用简明的解析解公式，不得不采用这种不准确的界面连接条件。 LUO 等[9-10]研究了初始温度不均匀条件下单孔地埋管换热器的传热，为了摒弃连接界面上的简化假定同时保留在钻孔外区域采用解析解公式的求解路线，对钻孔内的“热源”做了空间(深度方向)和时间上的离散化处理。 然而，此研究仅涉及单孔的地埋管换热器，而且对于深孔和长时间的问题，这一方法的计算工作量巨大，失去了传统解析解方法计算简捷的优势。 此外，由于采用解析解，不可避免还要有整体均匀介质的假定，对中深层地埋管而言误差较大。

地埋管换热器传热分析的第二类方法是以离散化数值计算为基础的传热模型，用有限元或有限差分法求解地下和流体中的温度响应并进行传热分析。 随着计算机技术的进步，数值计算方法以其适应性强的特点成为传热分析的基本手段，也已成为地埋管换热器理论研究的重要工具。 这种方法可以考虑比较接近现实的情况，如地下介质的不均匀性和地温梯度，在处理钻孔内外两个区域传热的耦合时，也可以不必采用简化假定而遵循比较真实的边界连接条件。 但是，地埋管换热器传热问题涉及空间范围的跨度大，从管壁厚度的毫米量级到深度的千米量级；同时问题涉及的时间范围的跨度也很大，负荷随时间的变化以及流动的特征时间都可以是分钟的量级，而系统运行的时间跨度长达十年以上。因此，对这一类传热问题的数值模拟，需要把区域划分为巨量的单元，耗费大量的计算时间。 单孔地埋管换热器的传热可以用柱坐标中的二维问题来描述，但求解由两个以上的钻孔组成的地埋管群的传热时，都将涉及三维非稳态传热问题，计算工作量往往大得无法忍受，进行基本的传热分析或直接求解工程问题都有很大的困难。 因此，现有的有关地埋管换热器(特别是中深层地埋管换热器)传热的数值分析研究主要涉及单孔地埋管换热器[11-13]。

文章是在尽可能真实的初始条件、边界条件(包括钻孔内外的连接条件)下，利用叠加原理把多孔地埋管换热器的三维传热问题分解成多个单孔的二维问题的叠加，实现降维的数值计算，从而使计算效率得到质的提升。 在研究单孔中深层地埋管换热器的传热时，开发了基于有限差分法的数值计算模型[14]。 在处理每一个时间步长上的二维瞬态导热问题时，由于采用了追赶法求解线性代数方程组，计算效率得到极大提高；同时在径向上进行了变步长的离散化处理，大大减少了节点的数量。 在计算机上，模拟单孔中深层地埋管换热器15 年的传热过程用时＜30 s，这就为采用叠加原理方法模拟多孔地埋管换热器的传热分析提供了基础模型。

1 地埋管传热问题的描述

中深层地埋管换热器的钻孔通常深入地下1 500～2 500 m，结构多采用套管式，如图1 所示。 换热器的外管采用钢管，内管应采用低导热系数的材料制成，现在通常采用耐热高密度聚乙烯管。 管内循环介质常采用水，循环水的流动方向有内进外出和外进内出两种。 从地下取热时，宜采用外进内出的方式；向地下蓄热时，宜采用内进外出的方式[13]。

图1 套管式地埋管换热器示意图

2 地埋管换热器传热数学模型

地埋管换热器的传热过程分成钻孔内、外两部分考虑，其各有不同的传热规律和特点，因此应分别建立数学模型。 两个区域的连接条件是界面上两者的温度和法向热流相等。 文章采用的重要方法是利用叠加原理把多孔地埋管换热器的传热问题分解为若干个单孔地埋管换热器传热的叠加，同时利用叠加原理处理带有温度梯度的初始温度分布。

2.1 钻孔内的传热模型

地埋管换热器钻孔内传热的基本特点是封闭管道(U 形管或套管)内流体在流动方向的对流传热与流体通过管壁和钻孔内材料与钻孔壁间径向传热的平衡。 通常需要给定换热器传热问题的定解条件，最常见的是给定流体的进口温度或换热量随时间的变化。在已知钻孔壁的温度分布以后，钻孔内的传热实际上就可以独立计算。 对U 形管或套管式的地埋管换热器的传热模拟已经开发了很多模型[15-17]，从最简单的利用热阻概念估算的方法到对管道内的流动和对流换热进行详细模拟的计算流体动力学(Computational Fluid Dynamics，CFD)软件，这些精细程度和准确度不同的各种模型应该都可以成为多孔地埋管换热器传热模型的一个组成部分。不过，对于地埋管换热器所涉及的超长直管道内的流动和传热，一维模型的精度已经可以满足实际需要，因此对于管内的流动和对流换热不建议采用复杂而费时的三维CFD 模拟。

在传统的地埋管换热器解析解模型中，钻孔内的传热通常采用所谓的“准三维”模型，也就是在考虑钻孔内横截面上的二维(或一维)导热的基础上，考虑流体温度在深度方向的变化以及轴向的对流换热。 由于钻孔换热器结构特别细长的特点，为保持模型的简明，可忽略钻孔内固体部分的轴向导热；同时，钻孔内材料(包括回填材料、管材和管道中的水)的热容量也可忽略不计。 后一个简化假定实际上把钻孔内的传热当作稳态问题处理，主要是为了使解析求解成为可能。 但这个假定明显是对实际情况的偏离，特别是在系统启动和负荷波动时会引起较大的偏差。 文章所提出的方法总体上将采用数值计算方法求解，在钻孔内模型的方程中加入流体温度随时间的变化，同时采用简化假定，忽略内管壁与内管流体以及外管壁、回填材料与外管流体各部分的温度差异。 这样的假定在不增加变量的数量的条件下，提高了模型的保真度，同时对计算工作量几乎没有影响。 因此，在这一传热模型中，钻孔内的传热采用考虑流体温度随时间变化和材料热容量的修正的准三维模型。 文章给出套管式地埋管换热器钻孔内传热的数学模型，就是兼顾了模型的精度和计算效率的模型，根据已知的上一时刻钻孔壁的温度分布，确定内外管中流体在该时刻的温度分布。 套管式地埋管换热器可以有外进内出和内进外出两种流动形式，并且不同的流动形式将导致不同的温度分布和传热量。

根据能量的平衡，可以写出流动形式为内进外出时的能量方程，由式(1)表示为

式中C1、C2分别为单位长度的外管和内管管段的热容量，J/(m·K)；z为深度方向坐标，m；τ为时间，s；Tf1和Tf2分别为外管和内管中的流体温度，℃；c为流体的比热，J/(kg·K)；M为管内流体的质量流量，kg/s；H为钻孔深度，m；R1和R12分别是外管流体与钻孔壁之间、内外管流体之间的热阻，m·K/W；且R1和R12由式(2)表示为

式中k1、k2分别为外、内管的导热系数，W/(m·K)；rb为钻孔半径，m；r1o、r1i分别为外管外径和内径，m；r2o、r2i分别为内管外径和内径，m。

该微分方程组的定解条件是内、外管在钻孔底部贯通，温度相同，由式(3)表示为

进口的条件可以给定温度，由式(4)表示为

或给定传热量(这里规定从地下取热为正)，由式(5)表示为

式中Q为取热量，W。

当流动形式改为外进内出时，方程式(1)中右边最后一项，即对流项要改变符号，定解条件式(4)和(5)也作相应的改变。

对于这个微分方程组，可以对坐标和时间作离散化处理后直接数值求解。 对于钻孔内的传热过程，采用更复杂或更简单的不同的数学模型进行模拟并不会影响以下所述的用叠加原理方法求解多孔地埋管换热器钻孔外传热的有效性。

2.2 钻孔外的传热模型

钻孔外的传热为空心圆柱形域内(如图1 所示)的导热问题，对钻孔外区域的传热采用以下基本假设:

(1) 钻孔外区域岩土中的传热机理仅考虑导热，而不考虑流体渗流等引起的其他传热；

(2) 钻孔外区域由水平分层的若干地层组成，各层均为均匀介质，其热物理性质(导热系数ki、热扩散率ai)不随温度变化；

(3) 在所研究的区域内有均匀的大地热流qg，W/m2；

(4) 在水平方向的外边界远离地埋管换热器的钻孔位置，为绝热边界；

(5) 在初始时刻域内各处未受扰动，保持稳态温度分布。

研究多孔地埋管换热器的传热必须先分析单孔的情况。 在以下的数学描述中没有涉及岩土分层的问题，以保持表达式的简洁。 实际上，虽然这种分层将对解析求解造成极大的困难，但对于数值求解几乎没有影响。

单孔地埋管换热器的传热问题可采用柱坐标描述，其数学描述包括导热微分方程与边界条件。 导热微分方程由式(6)表示为

式中T为温度，℃；r为径向坐标，m；rboundary为径向边界的半径，m；zboundary为深度方向边界的坐标，m。

地表通常采用第一类边界条件，即假设地表温度为定值t0，由式(7)表示为

区域的底部边界设置在远离钻孔底部的深处，采用第二类边界条件，由式(8)表示为

区域的径向边界设在离钻孔足够远处，其为绝热边界，由式(9)表示为

在钻孔壁上，是钻孔内、外两个区域的连接条件，分别由式(10)和(11)表示为

钻孔底部以下的岩土部分，简化为径向边界条件，由式(12)表示为

初始条件:岩土中的初始温度分布通常采用未受干扰时温度分布，满足稳态导热方程。 当假设岩土为均匀介质时，则初始温度分布由式(13)表示为

当假设岩土层由多个水平地层组成时，则初始温度分布由式(14)表示为

式中L为不同地层；zL为第L层地层的下边界的坐标。

3 叠加降维算法

对于以上数学模型所描述的单孔地埋管换热器的传热问题，利用线性叠加原理，钻孔外岩土区域的温度响应的解可以表示成两部分温度响应的叠加，由式(15)表示为

式中θ可以看作是岩土中任一点的温度相对于其自己的初始温度的过余温度。

把以上表达式代入式(6)～(14)的定解问题，可以看到θ(r，z，τ) 将满足以下定解问题，由式(16)表示为

即θ(r，z，τ) 将满足导热微分方程、零初始条件和上、下两个边界以及径向远端边界的齐次边界条件，同时也满足钻孔壁上的连接边界条件；这个二维导热问题可以结合钻孔内的数学模型，如式(1)所示，进行快速的数值求解[13]。

对于多孔地埋管换热器的传热问题，多孔在平面上的任意布置如图2 所示。 对于n个钻孔组成的地埋管换热器群，此时地下岩土中的温度响应必须用三维函数表示。

图2 多孔地埋管换热器的平面配置图

可以证明，多孔温度响应的叠加表达式(17)充分满足除式(11)以外的式(6)～(14)。 此时对于第i个钻孔壁的边界条件仍是在钻孔壁上热流连续，即在钻孔壁上通过导热进入岩土的热量等于钻孔内流体传给钻孔壁的热量，对于地埋管群中的第i个钻孔，式(11)应表示为

将式(17)代入以上边界条件得式(19)为

钻孔外岩土中的任一点(x，y，z)在任一时刻τ的温度响应可以由式(17)表示为

式中Гi是钻孔壁在横截面上投影的圆周边界。的物理意义是在单位长度上第j个钻孔中的换热(热源)所产生的热流进入第i个钻孔的净热量；第一项是第i个钻孔本身的内部传热产生的径向热流。 从物理概念上可以容易理解，从数学上也可以证明，在各钻孔间的距离远大于各钻孔半径的条件下， 即满足rij ＝1，2，…，n；i≠j) 时，式(19)中的第二项(求和项)中的各项都远小于第一项，可以忽略不计。 此时每一个单孔温度响应θi(ri，z，τ) 都满足其钻孔壁上如式(11)的边界条件。 对于几乎所有的多孔地埋管换热器的工程问题，钻孔之间的间距远大于钻孔半径的条件都是充分满足的，因此该多孔地埋管换热器的传热问题可以满足叠加原理的条件，也就是地埋管群传热问题的解可以用各单孔温度响应的解，即定解问题式(16)的解，叠加得到。 根据这个思路，多孔地埋管换热器数值计算流程如图3 所示。

图3 多孔地埋管换热器降维数值计算流程图

对于钻孔内的计算，主要关注钻孔壁的温度。 第j个钻孔的钻孔壁上的温度分布由式(20)表示为

式中rij为第i个钻孔与第j个钻孔中心之间的距离，

这样，多孔地埋管换热器的三维传热问题就可以通过求解若干个单孔地埋管换热器的二维传热问题得到解决。 得到的计算结果可以是任意时刻岩土中的温度分布，也可以给出工程上更关注的任意时刻任一地埋管换热器的出口温度。

4 算例

根据文章提供的降维计算方法开发计算程序，采用有限差分法求解。 算例考虑的系统是由两个钻孔组成的套管式中深层地埋管换热器，两个钻孔的几何参数、物理条件完全相同，埋设在由3 层水平地层组成的岩土中。 系统的主要参数如下:钻孔直径d＝0.28 m、深度H＝2 000 m、钻孔间距D＝10 m、大地热流qg＝0.075 W/m2、单孔循环水流量M＝12 kg/s、单孔取热速率恒定Q＝200 kW。 在取热负荷不变的条件下模拟了系统运行8 000 h 的传热工况。

采用文章提出的降维算法，可以把三维问题分解为若干柱坐标中的二维问题求解。 柱坐标中的二维导热问题在径向采用变步长，纵向由于温度梯度很小，可以采用较大的步长。 整个问题(包括钻孔内)采用的节点数为10 752 个，完全可以保证计算这样的大型三维传热问题的精度。 经测试，在台式电脑英特尔酷睿(i7-8700CPU＠320GHz 3.19GHz)上完成以上运算用时10.412 s。

对算例进行计算的结果如图4 所示。 其中，系统初始时刻温度场温度分布如图4(a)所示，系统结束运行时(8 000 h)岩土中的温度场分布如图4(b)所示。 随着地埋管换热器不断取热，钻孔周边岩土层的温度逐渐降低。 靠近埋管区域的岩土温度降低得最为剧烈，远离地埋管换热器的岩土温度逐渐趋于初始地温。 当运行8 000 h 时，双孔之间的岩土温度明显低于外侧岩土温度，说明二者之间存在热干扰现象，这对埋管换热器的取热效率有不利的影响。

图4 双钻孔地埋管换热器的岩土温度响应图

单孔取热及双孔取热两个不同工况下，流体进出口温度随运行时间的变化如图5 所示。 在恒定取热工况的条件下，单孔及双孔的进出口水温均随着运行时间的增大而逐渐降低。 持续取热的时间小于约半年时，两者的进出水温度几乎完全相同，说明在运行时间足够短时，双孔之间的热影响可以忽略不计。 计算结果表明，在第1 年的运行期间，两工况的进水和出水温度基本一致，说明在运行时间足够短的情况下，可以忽略相邻钻孔间的热干扰。 这充分表明了降维算法的可信度。

随着运行时间的增加，两者之间的差别逐渐加大。 在两孔的间距为10 m 的情况下，持续运行10 年后，双孔和单孔换热器的进出口水温可以相差约6.1 ℃。 这表明两个相邻地埋管换热器之间的热干扰已经严重影响到了地埋管换热器的效能。

图5 双孔及单孔地埋管换热器进出口水温随时间的变化趋势图

5 结论

文章基于线性叠加原理提出了多孔地埋管换热器的降维计算方法及地埋管群的传热模型，摒弃了传统地源热泵传热模型中关于钻孔内外区域连接条件的均匀温度和均匀热流的简化假定，因而更准确地描绘了中深层地埋管中的传热过程。 主要结论如下:

(1) 地下岩层中初始温度分布的影响可以由单独的一项叠加到零初始温度的温度响应中。 因此，叠加原理可以扩展到考虑地温梯度的场合。

(2) 单个钻孔地埋管换热器钻孔外的导热通常可以简化为柱坐标中的二维问题；而多个钻孔地埋管换热器钻孔外的导热必须要按三维问题处理。 在中深层地埋管群的传热问题中，线性叠加原理可以适用于钻孔外区域的导热问题。 因此，地埋管群的传热计算可以转化为多个单孔地埋管换热器传热(二维问题)的叠加来求解，从而大大提高计算的效率。

(3) 基于有限差分法的单孔中深层地埋管换热器的二维数值计算模型，采用了追赶法求解相应的线性代数方程组，具有计算高效快捷的特点，特别适合在用叠加原理求解多孔地埋管换热器的场合使用。

(4) 相邻的地埋管换热器在长期运行过程中存在热干扰现象，其影响的大小与钻孔间距、运行时间以及取热速率有关。 关于中深层地埋管群的参数分析以及设计和优化运行将在后续的研究中给出。