三维与二维结合地源热泵地下换热器复合数值模型

T.M.Eikevik

(1 上海交通大学机械与动力工程学院 上海 200240； 2 挪威科技大学 特隆赫姆 7491)

鉴于我国面临的严峻能源和碳排放形势[1]，地源热泵技术以低能耗特点如今成为传统空调系统的替代方案之一。地埋管换热器作为地源热泵系统最主要的部件之一备受研究者关注并由此演化出诸多模型。L. R.Ingersoll等[2]提出的无限线热源模型属于最早期的地埋管模型。该模型计算简便但只考虑了土壤中径向的热交换而忽略了轴向的热交换。而这一点被D. Marcotte 等[3-5]证明了重要性。此外，无限线源模型鉴于其积分式特点只适用于相对管径较小的场合[6]。另一种广泛使用的解析模型是L.Ingersioll等[7-9]提出的无限圆柱热源模型。该模型由无限线热源模型发展而来，并继承了无限线热源模型纯热传导假设和不考虑轴向热交换的特点，因而同样不适用于长期模拟[10]。为了弥补这一缺陷，P. Eskilson[11]引入了有限线热源模型。P. Eskilson的模型考虑了轴向换热及多管叠加的热效应，但其缺陷在于计算相对耗时，且和其他解析模型一样设定沿管长的热流密度为一定值。

相较于解析模型，数值模型通常有高精度和对应实际工程场景可定制的特点。然而，数值模型尤其是三维数值模型对计算量的要求极高[12-13]。因此多篇论文采用了不同方法来降低计算耗时，这一尝试主要分成普通三维模型和准三维模型两类：一类是以T.Y.Ozudogru 等[14]为代表的普通三维模型，该模型通过一维线性单元来简化模拟过程却刻意忽略了管内液体的流动。不考虑流动的影响在流速较大、管长较长的情况下尤为明显。另一类是以Li Zhongjian[15]为代表的准三维模型，该模型考虑到模型的对称性因而只计算了其中一边。这一类准三维模型[16-18]的主要特点是不考虑地表和埋管底部的影响直接把模型简化为水平二维模型，同样不能反映轴向热交换，因而不适用于长期模拟。

综上所述，这些为了追求减少计算耗时的简化削弱了原三维模型的准确性。很多简化只是局部性的，不能从根本上改变三维模型耗时的特点，无法满足不规则管群温度分布的模拟。由此可见，单独一个的三维模型或二维模型不能很好地模拟真实气候条件下的实际问题。

本文提出了一个三维与二维结合地源热泵地下换热器复合数值模型，能引入实际气候数据的影响，获得不规则管群周边温度场分布情况，且提高计算效率。在新建的复合模型中，三维模型导出的模拟结果被用作二维模型输入的边界条件。结合两者的复合模型既发挥了三维模型准确性高的优点，又通过二维模型显著降低了管群模拟的运算量。

1 换热器模型的建立

复合模型由一个三维模型和一个二维模型组成。三维模型由两组U形管、管内流体、回填材料、管外及表层土壤组成。U形管外土壤的计算区域半径为2 m，计算区域深度距U形管底部2 m，具体结构如图1所示，换热模型参数如表1所示。二维模型可根据实际管群分布建立。实际管群分布为呈直角分布的两组矩形埋管阵列，管外土壤的计算半径为4 m，如图2所示。选择阵列中最具代表性的共27根地埋管组成的直角区域进行建模。为进一步简化，二维模型采用单位管化，相应的参数也做了适当调整[19]，相应的网格如图3所示，其中离U形管较近的网格较为紧凑，离U形管较远的网格较为稀疏。复合模型中，三维模型在运行状态下不同深度的钻孔壁温度模拟结果可视为二维模型的钻孔壁温度；其导出的钻孔壁暂态温度序列可被植入二维模型的边界条件函数。因此，在计算管群埋管区域某一个深度范围内的温度分布时可以针对性的选取三维模型该范围内等间距的输出结果作为输入边界条件带入二维模拟，可避免运行三维管群模型需要计算每一层网格所负担的庞大计算量。换热器模型建立过程中所作假设如下：1)地表温度随大气温度变化；2)管内换热流体温度仅随流动方向一维变化；3)土壤及回填材料为物性均匀；4)材料物性参数在模拟过程中保持不变；5)管间的相互影响仅限于垂直于U形方向。

图1 三维模型示意图Fig.1 Schematic diagram for the 3D model

系统在运行状态下流体暂态的总换热量等于因轴向流动引起的热量变化与径向导热引起的换热量之和，管壁与流体边界上的能量平衡方程为：

(1)

图2 二维管群模型示意图Fig.2 Schematic diagram of the boreholes distribution of the site in Hangzhou

图3 三维及二维网格模型Fig.3 The numerical mesh used to represent 3D and 2D model

系统停机状态下流体暂态的总换热量仅受到径向的导热的影响，其平衡方程为：

(2)

在U形管，回填材料和土壤内部的能量平衡方程为：

(3)

土壤的热力学能初值设为零，变化方程为：

(4)

式中：T为温度，K；τ为时间，s；m为质量，kg；r为半径，m；z为深度，m；u为流速，m/s；α为热扩散系数，m2/s；ρ为流体的密度，kg/m3；c为比热容，J/(kg·K)；h为表面传热系数，W/(m2·K)；下标w，p，f，s分别指代流体、U形管、回填材料以及土壤；下标pi，po，b，in分别指代U形管内、外管壁、钻孔壁及入口坐标位置。

假定各接触面紧密结合，土壤表面的温度，进水口的温度和速度均分别由温度函数序列Tudf(τ)和速度函数序列uudf(τ)给出。边界上的温度关系及输入条件为：

(5)

(6)

历史气候数据的获取有两种途径，国家气象局的网站或使用气象站数据。如果都没有，可以考虑建立气温的理想化周期波动曲线：

T(τ)=0.5Tmcos(2×10-7τ+φm)+

0.5Tdcos(2×10-7τ+φd)+Tavg

(7)

式中：τ为时间，s；Tm，Td分别为年周期及日周期温差,K；φm，φd分别为年周期及日周期的相位；Tavg为当地年平均气温，K。

二维模型传热方程的极坐标表达式为：

(8)

式中：θ为角度，(°)；k为导热系数，W/(m·K)。三维模型在运行状态下不同深度的孔壁温度模拟结果可视为二维模型的钻孔壁温度，因而其导出的钻孔壁暂态温度序列Tudf′(τ)可构成二维模型的钻孔壁输入边界条件：

Ts(τ)|r=rb=Tudf′(τ)

(9)

同时三维模型导出的不同深度稳定的土壤初始温度构成不同深度二维模型土壤初始温度。

三维及二维模型模拟流程如图4所示。

图4 三维及二维模型模拟流程Fig.4 The simulation process of 3D and 2D model

2 换热器模型的验证

为了验证本文采纳的三维模型暂态模拟方法的准确性，此处比照B.A.Beier等[20]提出的沙盒实验建立文献中同尺寸的地下换热器暂态数值模型。模型采用的具体几何与热物性参数见表2，同时模型采用了与实验相同的入口流速和温度作为输入边界条件。模拟预测的出水口暂态温度文与文献中传感器实测的出水口温度对比如图5所示。其平均绝对比例误差(MAPE)为3.13%，平均绝对误差(MAE)为0.416，相较传感器在测量时随机产生的标准差0.391 K，可认为所建立的三维模型能很好地反映实际情况。

图5 模型预测出水温度与实测出水温度对比Fig.5 Comparison for predicted and measured outlet temperature

3 模拟结果与分析

论文采纳了杭州市内的气象站数据，跨度为一个完整的供冷季、过渡季和供暖季。模型的模拟时间从2015年6月12日—2016年3月17日，期间分为供冷季、过渡季及供暖季。地源热泵在供冷季和供暖季中分别向地下释放热量和冷量，在过渡季长时间处于停机状态。从15年6月12日开始模拟直至16年3月17日的热力学能与气温数据如图6所示。地下换热器模型采纳的初始化条件为流体、U形管、回填材料和土壤各节点采用统一温度设定为291.8 K，通过将第一步的边界条件作为定值重复模拟50个假想日，可以有效缓冲初始化带来的误差影响。在模拟中第一个假想日的热力学能增量为0.039 4 MJ，供冷季开始前最后一个假想日的热力学能增量仅有0.001 9 MJ，热力学能随时间的斜率之比降至初始值的4.9%，近似于6月大气条件下的土壤热力学能自然增量的速度。表明经过50个假想日的前置稳定期后，土壤已基本达到适宜后续模拟运算的温度分布。

表2 沙盒实验参数

图6 土壤热力学能与气温对照数据Fig.6 Data of soil internal energy and air temperature

图6中大气温度曲线在8月上旬(8月4日14∶00)达到峰值后开始下行并在一月下旬(1月16日6∶00)达到谷值。地埋管在供冷季向地下释放热量并在供暖季从地下抽取热量。对应的土壤热力学能曲线在供冷季上升，在过渡季中基本保持不变，并在供暖季开始后逐步下降。与仅受气温影响的热力学能对比可见系统运转下的内能曲线变化在供冷季和供暖季明显更加陡峭。在供冷季结束后整个三维计算区域的热力学能上升了3.96 MJ，平均温升为1.09 K。一个供冷季、过渡季、供暖季结束后，整个三维计算区域的热力学能上升了0.292 MJ，与仅受气温影响的情况相比，土壤的相对热力学能增量为1.037 MJ，相对温度增量为0.28 K，如图7所示。这反映了长江三角洲地区每年对冷量的需求高于对热量的需求，因而在年复一年的换热过程后，将有大量的热量被堆积在土壤中，进而影响地埋管换热器的换热效率。

图7 不同深度观测点的土壤温度Fig.7 The soil temperature at at different depth observation points

指定-0.5 m及-1.5 m的进出水管连线中心轴上的观察点，考察其全年的温度变化情况。-1.5 m的观测点的平均温度及和仅受气温影响情况下的相对温升峰值滞后于-0.5 m的观测点约2周的时间。计算区域内土壤的平均温度和仅受气温影响情况下的相对温升峰值滞后于-0.5 m的观测点约9周的时间。说明沿埋管深度深的绝对温度和相对仅受气温影响情况下的温升曲线均在相位上滞后于沿埋管深度较浅的点。

图8 进出水口温度对比Fig.8 Comparison of the inlet fluid temperature and outlet fluid temperature

杭州市内的地源热泵系统在仅在需要供冷或供热季节的工作时间开机，在多数周末、节假日停机。尤其在9月末到11月初这段气温适宜的过渡季期间，系统有较长的停机时段。在系统停机期间，土壤各点的温度及进出口处的水温仅受大气及环境温度的影响，显示出较为平滑的曲线。由图8整理得知，当系统运行时，整个供冷季地埋管的出水口水温与同时刻气温温差均值为-8.313 K，整个供暖季地埋管的出水口水温与同时刻气温温差均值为9.077 K。

图9所示为在完整的供冷季、过渡季、供暖季结束后埋管周围的温度分布情况，埋管区域一侧及转角处受到热堆积的最大影响半径见表3。由表3可知，沿埋管深度越深受热堆积影响的范围越小，埋管区域在转角处的热堆积影响半径显著大于埋管区域一侧的影响半径，其影响半径的比值随深度减小，反映出不规则排列地埋管群的热堆积在较深的地层的影响范围更局限于其自身的排列形状。通过用系统运转条件下的截面平均温度减去仅受大气影响条件下的平均温度可得不同深度截面的平均温升，如图10所示。-1 m深和-10 m以下的截面有不同的均温变化规律，并且-1 m深度在整个供冷季的均温峰值低于-10 m的均温峰值。这一现象表现出与大气的换热在带走土壤热堆积尤其是表层土壤热堆积方面的显著影响。结合三维模型的模拟结果可知，系统运转所造成的热堆积效应主要集中在7 m以下稍深的土壤层而非与大气频繁换热的表层土壤(最大温差0.5 K)。

表3 不同深度处的最大影响半径

图9 供暖季结束时不同深度的二维温度分布Fig.9 Temperature distribution at different depth at March 17th, 2016

图10 截面平均温升随不同深度的变化Fig.10 Cross section average increment of temperatureat changes with different depth

3 结论

本文建立了基于实际气候条件下的三维与二维结合地源热泵地下换热器复合数值模型。实现了用较小的计算耗时处理模拟管群换热的问题。相较于现有模型，新建立的复合模型利用三维模型的暂态模拟结果作为二维模型输入的边界条件，进而在降低管群运算复杂度的同时保证了精确性。案例中，采用了杭州市内的实测数据，模拟时间为2015年6月12日—2016年3月17日，跨度为一个完整的供冷季、过渡季和供暖季，主要结论如下：

1)当计算三维换热时，不需要采用设置初始温度的方法，仅需设定初始气温，将第一步的边界条件作为定值重复模拟50个假想日的时长作为缓冲，可将热力学能随时间的日均斜率降至初始值的4.9%，有效冲抵了初始化时全部节点设置同一个初始温度带来的误差干扰。

2)影响半径为2 m以内的土壤在模拟一个完整的供冷季、过渡季和供暖季后的热力学能增量为0.292 MJ，相较仅受气温影响的土壤的相对热力学能增量为1.037 MJ，相对温度增量为0.28 K。该增量反映出长三角地区用户对冷量的需求高于对热量的需求，因此释放到土壤中多余的热量将随时间逐步堆积，对换热效率产生不利影响。

3)沿埋管深度深的绝对温度和相对仅受气温影响情况下的温升曲线均在相位上滞后于沿埋管深度较浅的点。在进出水口连线中心轴上-1.5 m的观测点的平均温度及和仅受气温影响情况下的相对温升峰值滞后于-0.5 m的观测点约2周的时间。计算区域内土壤的平均温度和仅受气温影响情况下的相对温升峰值滞后于-0.5 m的观测点约9周的时间。

4)系统运行条件下，整个供冷季地埋管的出水口水温与同时刻气温温差均值为-8.313 K，整个供暖季地源热泵的出水口水温与同时刻气温温差均值为9.077 K。

5)系统运转造成的热堆积效应主要集中在7 m以下稍深的土壤层而非与大气频繁换热的表层土壤。埋管区域在转角处的热堆积影响半径大于其一侧的影响半径，且其绝对值与比值均随深度的增大而减少，说明不规则排列地埋管群的热堆积的影响范围在越深的地层越局限于自身的排列形状。

[1] 中华人民共和国国家统计局. 2014年中国统计年鉴[M].北京: 中国统计出版社,2014：54-56. (National Bureau of Statics of China. China statistical yearbook in 2014[M]. Beijing: China Statistic Press, 2014：54-56.)

[2] INGERSOLL L R, PLASS H J. Theory of the ground pipe heat source for the heat pump[J].Heating Piping & Air Conditioning, 1948, 54(7): 339-348.

[3] MARCOTTE D, PASQUIER P, SHERIFF F, et al. The importance of axial effects for borehole design of geothermal heat-pump systems[J]. Renewable Energy, 2010, 35(4): 763-770.

[4] ZENG H Y, DIAO N R, FANG Z H. A finite line-source model for boreholes in geothermal heat exchangers[J]. Heat Transfer—Asian Research, 2002, 31(7): 558-567.

[5] MOLINA-GIRALDO N, BLUM P, ZHU K, et al. A moving finite line source model to simulate borehole heat exchangers with groundwater advection[J]. International Journal of Thermal Sciences, 2011, 50(12): 2506-2513.

[6] YANG H, CUI P, FANG Z. Vertical-borehole ground-coupled heat pumps: A review of models and systems[J]. Applied Energy, 2010, 87(1): 16-27.

[7] INGERSIOLL L, ZOBEL O J, INGERSOLL A C. Heat conduction: with engineering geological and other applications[M]. US: University of Wisconsin Press, 1954.

[8] KAVANAUGH S P. Simulation and experimental verification of vertical ground-coupled heat pump systems[D].US: Oklahoma State University Stillwater, 1985.

[9] CARSLAW H S, JAEGER J C. Conduction of heat in solids[M]. Oxford: Clarendon Press, 1960.

[10] REES S J, HE Miaomiao. A three-dimensional numerical model of borehole heat exchanger heat transfer and fluid flow[J]. Geothermics, 2013, 46: 1-13.

[11] ESKILSON P. Thermal analysis of heat extraction boreholes[M]. 1987.

[12] HE Miaomiao, REES S J, SHAO Li. Simulation of a domestic ground source heat pump system using a three-dimensional numerical borehole heat exchanger model[J]. Journal of Building Performance Simulation, 2011, 4(2): 141-155.

[13] HE Miaomiao. Numerical modelling of geothermal borehole heat exchanger systems[D].Leicester: De Montfort University, 2012.

[14] OZUDOGRU T Y, OLGUN C G, SENOL A. 3D numerical modeling of vertical geothermal heat exchangers[J]. Geothermics, 2014, 51: 312-324.

[15] LI Zhongjian. A new constant heat flux model for vertical U-tube ground heat exchangers[J]. Energy and Buildings, 2012, 45: 311-316.

[16] ZHANG Linfeng, ZHANG Quan, HUANG Gongsheng. A transient quasi-3D entire time scale line source model for the fluid and ground temperature prediction of vertical ground heat exchangers (GHEs)[J]. Applied Energy, 2016, 170: 65-75.

[17] MA Weiwu, LI Min, LI Ping, et al.New quasi-3D model for heat transfer in U-shaped GHEs (ground heat exchangers): Effective overall thermal resistance[J]. Energy, 2015, 90: 578-587.

[18] WETTER M, HUBER A. Vertical borehole heat exchanger EWS Model[J]. Trnsys Type, 1997, 451.

[19] GU Y, O′NEAL D L. Development of an equivalent diameter expression for vertical U-tubes used in ground-coupled heat pumps[J]. Ashrae Transactions, 1998, 104(2): 347-355.

[20] BEIER R A, SMITH M D, SPITLER J D. Reference data sets for vertical borehole ground heat exchanger models and thermal response test analysis[J]. Geothermics, 2011, 40(1): 79-85.