我国是化石燃料消耗大国，向国际社会承诺的双碳战略目标也给我国新能源和可再生能源的开发利用带来了机遇与挑战。地热能是一种储量丰富的可再生能源，近年来我国浅层地热能的开发利用增速较高，总装机容量2020年底已达到2.6万MW，连续20多年处于世界领先地位^[1]。与空气源热泵相比，地源热泵系统（ground source heat pump system, GSHPs）性能系数较高且运行稳定。然而，由于地源热泵的设计单位参差不齐，系统运行过程中出现了一些问题，如长期运行后系统冷热负荷不平衡导致的热/冷堆积现象。地源热泵系统达不到预期效果有多种原因，如地层热物性参数测量不准确、系统冷热负荷不匹配、运行管理不善等。一些设计能力较好的企业通常采用国外商用软件，如 ANSYS^[2]、TRNSY ^[3]等进行设计，这些软件采用的数值计算方法通常为有限体积法^[4-5]和有限差分法^[6]。但这些软件的使用需要掌握一定的专业知识，且模拟过程比较耗费机时。

地埋管换热器（borehole heat exchanger, BHE）是地源热泵系统中的关键部件，较早的单根埋管设计采用的理论方法是无限长线热源^[7]，此后出现了考虑井筒直径影响的无限长柱热源模型^[8-9]，以及考虑埋管深度影响的有限长线热源^[10]或柱热源模型^[11]等。为研究地埋管群之间的传热影响，Eskilson^[12]提出了无量纲g函数，Bernier等^[13-14]利用该方法及线性叠加原理，提出了无限长柱源模型下的井群热干扰计算方法。Marcotte 等^[15]结合线性叠加原理、时间卷积与快速傅里叶变换方法，提出了可以耦合地面“小时”级热负荷变化的地层温度计算方法。Man等^[16]同样基于短时间热负荷阶跃变化温度响应及线性叠加原理，采用有限长线源模型对埋管周围的地层温度变化及出水温度进行了计算。此外，Zhang等^[17]、Zhang等^[18]采用g函数与线性叠加原理，计算了有地下水渗流及地埋管间歇运行情况下的地埋管温度场。Wang等^[19]研究了地埋管冷热负荷不平衡对地层温度场及系统经济性的影响。杨露梅等^[20]在南京典型地埋管热泵工程布设监测孔拟合了地温场的变化方程。刘爱华等^[21]布设了观测孔以研究地埋管换热区地温垂向深度的影响。杨卫波等^[22]建立了能量桩热-力耦合模型研究了地下水渗流对桩埋管的影响。陈宝义等^[23]研究了土层挤密条件对换热器换热效果的影响。张延军等^[24]利用2 m测温法对厦门东山某热泉地区进行实地探测。

目前为止，针对有动态热冷负荷变化的地埋管地层温度场的解析解大部分都是基于Eskilson^[12]提出的线热源g函数结合线性叠加原理。该法是将随时间变化的实际热流曲线分解为多个不同的恒定热流及不同启动时间的g函数解的叠加。这一方法虽较数值方法计算简单，但也存在以下缺点：（1）g函数针对的是线源模型，未考虑井筒内径带来的影响；（2）当热负荷随时间变化较剧烈时，需要更多或间隔时间更短的不同启动时间下的阶跃g函数，以满足计算精度，但计算量也随之增加。本文提出了一种更为简便且可快速计算BHE长期年周期运行下的地层温度的解析算法。该法基于无限长圆柱模型，首先求解井筒壁面有周期性热流密度下的地层温度分布的解析解，再将年周期内地埋管热冷负荷随时间任意变化的曲线进行傅里叶或切比雪夫变换，最后将每个函数对应的解析解与其波幅相乘叠加，即可得到任意时刻的地层温度。由于该模型直接设定井筒壁面边界条件为周期性热流边界，对地埋管运行初期的地层温度计算会有一定误差，但随着地埋管的长期运行，该法快速计算得到的地层温度与数值解吻合良好，而用户更为关心的设计指标也正是地埋管长期运行后周围地层的温度场和热干扰半径。故与现有解析解和数值解相比，此法分别具有计算精度高和快速便捷的优点，可推广应用于工程实际。

1.   地埋管周期性热负荷边界下的无限长圆柱模型

1.1   周期性热流边界下无限长圆柱地层温度解析解

一维径向轴对称导热模型如图1所示，假设：

图  1  无限长柱源导热模型示意图

注：q为井内壁面的热流密度；R_w为井筒半径；r为地层半径；R_∞为地层远端半径；T_∞为地层远端温度。

Figure  1.  Schematic diagram showing the heat transfer model of BHE

下载: 全尺寸图片 幻灯片

（1）无限长圆柱径向导热轴对称；

（2）初始温度是常数且均匀；

（3）忽略地下水渗流，地层物性各向同性且为常数；

（4）井内壁面的热流密度（q）随时间（t）呈正弦或余弦周期性变化，即$q\left( t \right) = A\sin \left( {wt + \varphi } \right) + B$, 其中：A是波幅，B是偏离量或周期内平均值，w是年周期下的频率（w=2π/t_m），t_m为年周期（12个月），$\varphi$为相角。

井筒周围地层的能量控制方程为：

式中：角标“s”——地层；

$\rho$——密度/（kg∙m⁻³）；

C_p——比热/（J∙kg⁻¹∙K⁻¹）；

T——温度/°C；

λ——导热系数/（W∙m⁻¹∙K⁻¹）。

一般偏微分方程的解析解可通过分离变量法、积分变换法等得到。由于本方程包含有井筒壁面的非齐次边界条件，分离变量法不适用，因此，本文采用拉普拉斯积分变换（拉氏变换），先消去时间项得到修正的贝塞尔方程，再通过Gaver-Stehfest 进行拉氏逆变换。进行拉氏变换之前，将温度变量变换为过余温度（θ），即定义$\theta = ({T_{\text{s}}} - {T_\infty })$，拉式变换后的方程形式为^[9]：

式中：$\bar \theta$——拉氏变换后的过余温度/℃；

s、a——系数，s为复数，$a = {{{\lambda _{\text{s}}}} /{(\rho {C_p})_{\rm{s}}}}$，且${s / a}= {\varepsilon ^2}$。

式（3）为贝塞尔方程，其解的形式为：

式中：I₀——修正后的第一类贝塞尔函数；

K₀——修正后的第二类贝塞尔函数；

C₁、C₂——系数，由边界条件式（4）计算。

代入边界条件式（4）后得到经拉氏变换后的过余温度（$\bar\theta$），再通过拉氏逆变换可得原始过余温度（θ），如式（5）所示。式（5）的解可以分解为三部分，分别为θ_A1、θ_A2和θ_B，如式（6）—（8）所示。波幅A的乘积系数包括两项，对应的解分别设为θ_A1和θ_A2；常数B项的乘积系数只有一项，对应的解设为θ_B，即定常热流密度下的解析解^[9]。

式中：i——虚数；

K₁——一阶第二类修正贝塞尔函数；

$\gamma$、f(u)——自定义函数：$\gamma = \sqrt {{w /a}}$， $f(u)=$$\dfrac{{{J_0}\left( {ur} \right){Y_1}\left( {u{R_{\text{w}}}} \right) - {Y_0}\left( {ur} \right){J_1}\left( {u{R_{\text{w}}}} \right)}}{{{J_1}^2\left( {ur} \right) + {Y_1}^2\left( {u{R_{\text{w}}}} \right)}}$；

J₀、J₁——零阶、一阶第一类贝塞尔函数；

Y₀、Y₁——零阶、一阶第二类贝塞尔函数；

M——零阶第二类修正贝塞尔函数的模：$M = \sqrt {{{\ker }_0}{{\left( {\gamma r} \right)}^2} + {\text{ke}}{{\text{i}}_0}{{\left( {\gamma r} \right)}^2}}$；

N——一阶第二类修正贝塞尔函数的模：$N = \sqrt {{{\ker }_1}{{\left( {\gamma r} \right)}^2} + {\text{ke}}{{\text{i}}_1}{{\left( {\gamma r} \right)}^2}}$；

θ₁——零阶第二类修正贝塞尔函数的相角：${\theta _1} = \arctan \left[ {{{{\text{ke}}{{\text{i}}_0}\left( {\gamma r} \right)} \mathord{\left/ {\vphantom {{{\text{ke}}{{\text{i}}_0}\left( {\gamma r} \right)} {{{\ker }_0}\left( {\gamma r} \right)}}} \right. } {{{\ker }_0}\left( {\gamma r} \right)}}} \right]$；

θ₂——一阶第二类修正贝塞尔函数相角补角：${\theta _2} = \arctan \left[ {{{{{\ker }_1}\left( {\gamma r} \right)} \mathord{\left/ {\vphantom {{{{\ker }_1}\left( {\gamma r} \right)} {{\text{ke}}{{\text{i}}_1}\left( {\gamma r} \right)}}} \right. } {{\text{ke}}{{\text{i}}_1}\left( {\gamma r} \right)}}} \right]$；

ker₀——零阶第二类修正贝塞尔函数的实数部： ${\ker _0}\left( {\gamma r} \right) = {Re} \left[ {{K_0}\left( {{{\text{e}}^{{\text{i}}{\text{π}}/4}}\gamma r} \right)} \right]$；

kei₀——零阶第二类修正贝塞尔函数的虚数部： ${\ker _0}\left( {\gamma r} \right) = {Im} \left[ {{K_0}\left( {{{\text{e}}^{{\text{iπ }}/4}}\gamma r} \right)} \right]$；

ker₁——一阶第二类修正贝塞尔函数的实数部： ${\ker _1}\left( {\gamma r} \right) = {Re} \left[ {{{\text{e}}^{ - {\text{iπ }}/2}}{K_1}\left( {{{\text{e}}^{{\text{iπ }}/4}}\gamma r} \right)} \right]$；

kei₁——一阶第二类修正贝塞尔函数的虚数部： ${\ker _1}\left( {\gamma r} \right) = {Im} \left[ {{{\text{e}}^{ - {\text{iπ }}/2}}{K_1}\left( {{{\text{e}}^{{\text{iπ }}/4}}\gamma r} \right)} \right]$。

因此，式（5）可以写为θ_A1、θ_A2和θ_B之和的形式：

1.2   周期性冷热负荷变化情况下的地温场求解

上文推导得到的式（9）是针对井壁面热流随时间满足完整周期性正弦函数情况下的解。但实际上，地埋管换热过程中管壁热流随时间的变化虽然是周期性函数，但并非一个完整的正弦函数。因此，如何将此周期性时间函数转换为正弦函数，利用上述推导得到的解是解决该问题的关键。本文利用傅里叶变换方法，将实际热负荷波动函数近似表示为有限个（n个）正余弦周期函数的和，同时利用导热微分方程解的可叠加性。具体表示如下：

式中：L——半个周期（6个月）；

n——展开的三角函数的个数。

将常数a₀带入式（8），将每一个${a_j}\cos \dfrac{{j{\text{π}}}}{L}t$与${b_j}\sin \dfrac{{j{\text{π}}}}{L}t$带入公式（9）进行计算累加得到地层温度

式中：${\theta _{\cos ,j}}(t)$、${\theta _{\sin ,j}}(t)$——对应j个余弦周期函数和正弦 周期函数的解。

2.   我国典型地区浅层地埋管周围地层温度场的快速求解

2.1   解析解模型验证

为了验证本文提出的解析解，本文以天津市为例，分析北方典型地区的地埋管换热器供热/冷负荷，并将数值模拟结果与本文解析解进行对比验证。选取的地埋管几何参数和地层热物性如表1所示。

表  1  地埋管几何与地层物性参数

Table  1.  Geometrical and physical parameters given in the simulation

参数名称	值	参数名称	值
Rw /mm	55	T0/°C	15
λs/（W·m−1·K−1）	2.2	(ρCp)s/（J·m−3·K−1）	2.16×106
tm （2L）/月	12	t/a	10
nm	20	nh	500
注：nm、nh为年、天周期傅里叶展开级数。

下载: 导出CSV 

| 显示表格

图2为2020年天津市的年平均气温变化。可以看出，天津一年的平均气温可近似成一个三角函数，日最高气温与最低气温也可近似为三角函数。由于室外气温与地埋管冷/热负荷变化的线性关系较强^[25]，故可根据室外气温变化得到地埋管井壁处的径向热流通量，可表示为年周期函数${q_{\text{m}}}({t_{\text{m}}})$和日周期函数${q_{\text{h}}}({t_{\text{h}}})$之和，即

图  2  天津市气温与地埋管延米功率的变化

Figure  2.  Ambient air temperature and heat flux at the borehole wall in Tianjin City

下载: 全尺寸图片 幻灯片

天津市集中供暖时间为每年11月15日至次年3月15日，供冷时间为每年5月15日至9月15日；日热量峰谷值可近似为以天（t_h）为周期的三角函数${q_{\text{h}}}({t_{\text{h}}})$。因此，从供冷时间开始，即开始运行年份的5月15日，t_m=0，全年冷、热负荷时长相同，冷、热负荷延米功率相同，冷、热负荷不平衡率较小故可认为冷、热负荷平衡。假设地埋管延米功率月负荷峰值为30 W/m，日负荷峰谷值为5 W/m，由此可近似得到天津市井壁处的总周期延米功率q_TJ的函数形式，如式（13）所示：

通过式 （13），可以得到解析解与数值解的对比结果，如图3所示。可以看到，在前3个周期，解析解和数值解之间偏差较大，这是因为解析解的求解忽略了初始条件；随着地埋管运行时间的增大，解析解和数值解之间的偏差逐渐减少，直至忽略不计。另外可以看到，在井壁处解析解和数值解之间的偏差波动较别处更为剧烈。这是因为此偏差受限于动态热负荷的傅里叶展开精度，增加傅里叶展开级数的项数可降低此偏差。此外，解析解和数值解之间的偏差会随着距管壁距离的增加而减小，距管壁较远处地层温度的2种解的偏差可忽略不计。式（13）中，总周期负荷可分为年周期负荷与日周期负荷两部分。对比总周期和年周期的解的偏差可以看出，靠近井壁处的总周期偏差大于年周期偏差，这同样受限于日周期负荷的傅里叶展开精度，可通过增加傅里叶级数展开的项数降低偏差。另外可以看到，距井壁较远处地层温度的2种解之间的偏差可忽略不计。由此可知，总周期负荷下的解析解适用于计算地埋管长期运行时的地层温度变化，运行3 a后的解析解与数值解的偏差小于0.1 °C。

图  3  不同径向距离下地层温度解析解与数值解的对比

注：$\theta_{\rm{ana}}$为解析解，$\theta_{\rm{num}}$为数值解。

Figure  3.  Comparison of analytical and numerical solutions of temperature changes and temperature deviations

下载: 全尺寸图片 幻灯片

2.2   天津市浅层地埋管周围的地层温度场

分别计算式（12）中的年周期函数q_m(t_m)与日周期函数q_h(t_h)对地温场的影响，结果如图4所示。可以看出，在q_h(t_h)峰谷值对称的情况下，q_h(t_h)的负荷波动只影响近井壁处的地层温度，且q_m(t_m)波动引起的地层温升为相同时间段内q_TJ的平均值，日周期延米功率q_h(t_h)引起的地层温升在q_m(t_m)的基础上对称波动，近井壁处的波动范围在0.4 °C以内，而距井壁较远处的地层波动趋于0，故较远处q_m(t_m)与q_TJ波动引起的地层温升相同。此外还可看出，地层温度主要由年周期负荷确定，日周期负荷对于地层温度的影响很小，相比于年周期函数的傅里叶展开级数，为得到更好的拟合精度，日周期函数需要更多的傅里叶展开项数，这会影响计算效率。因此计算过程中，可酌情忽略日周期负荷的影响，从而简化解析解。

图  4  不同半径处的地层温度随时间的变化

Figure  4.  Temperature changes at different radii of BHEs with time

下载: 全尺寸图片 幻灯片

2.3   其他城市浅层地埋管周围地层温度场

为研究地埋管冷/热负荷的均衡程度和不同运行时间下地层温度场的变化规律及热影响半径，本文同时选取哈尔滨市为研究对象。该地区以冬季热需求为主，即冬季热负荷运行时间长，夏季冷负荷运行时间短。而夏季冷负荷时间长，冬季热负荷时间短的广州市可作为哈尔滨市的反向类型。

通过对哈尔滨市2020年气温变化的统计，得到周期冷/热负荷周期函数，如图5所示，哈尔滨市法定供暖时间为每年10月20日至次年4月20日，供冷时间为每年6月20日至8月20日。自供冷时间开始计算，即开始运行年份的6月20日为t_m=0，由2.2节可知，日周期负荷对计算结果影响不大，因此此处忽略日周期负荷，可近似得到哈尔滨市井壁处的总周期延米功率（q_HB）。类似可以得到广州市的总周期延米功率（ q_GZ）。公式为：

图  5  哈尔滨市年气温与地埋管延米功率变化

Figure  5.  Air temperature and heat flux at the borehole wall in Harbin City

下载: 全尺寸图片 幻灯片

分别对式（14）与式（13）的年周期函数q_m(t_m)进行计算，研究年周期负荷的影响，结果如图6所示。可以看到，由于哈尔滨的热负荷远大于冷负荷，因此地埋管从地层提供热量大于提供冷量，运行过程中会造成冷堆积现象，导致地层温度逐年下降，系统运行10 a后，距地埋管5 m处的地层温度会降低2.2 °C，影响地埋管取热效率。反之，由于广州的冷负荷远大于热负荷，造成热负荷堆积，地层温度逐年上升，系统运行10 a后，距地埋管5 m处的地层温度上升2.1 °C。

图  6  不同地区相同年周期负荷下地层温度随时间的变化

Figure  6.  Comparison of temperature changes at different radii of BHEs in different area

下载: 全尺寸图片 幻灯片

2.4   总周期延米功率解析解对地层温度的影响分析

通过傅里叶级数展开，我们对天津市与哈尔滨市拟合的周期负荷进行展开，如表2所示。通过式（10）可以看到，其傅里叶展开系数有一定关联，a_j的结果为cos函数，与a_j相比，b_j的展开式中除A₁项系数项少了第二项外，其余系数项只需将a_j公式中的cos函数改为sin函数即可。

表  2  不同地区的傅里叶展开级数

Table  2.  Fourier expansion series in different areas

地区	${a_0}$	${a_{_j} } (j=1,\cdots,n)$	${b_{_j} } (j=1,\cdots,n)$
天津	$\dfrac{1}{6}\left(\dfrac{8}{\text{π} }{A}_{1}-\dfrac{8}{\text{π} }{A}_{2}\right)$	$\dfrac{ {6{A_1}\left[ {\cos \left(\dfrac{ {2j{\text{π } } } }{3}\right) + 1} \right] - 6{A_2}\left[ \cos \left(\dfrac{ {j{\text{π } } } }{3}\right) + \cos (j{\text{π } }) \right]} }{ {(3 + 2j)(3 - 2j){\text{π } } } }$	$\dfrac{ {6{A_1}\sin \left(\dfrac{ {2j{\text{π } } } }{3}\right) + 6{A_2}\left[ {\sin \left(\dfrac{ {j{\text{π } } } }{3}\right) + \sin (j{\text{π } })} \right]} }{ {(3 + 2j)(3 - 2j){\text{π } } } }$
哈尔滨	$\dfrac{1}{6}\left(\dfrac{4}{\text{π} }{A}_{1}-\dfrac{12}{\text{π} }{A}_{2}\right)$	$\begin{array}{l}\dfrac{12{A}_{1}\left[\mathrm{cos}\left(\dfrac{j\text{π} }{3}\right)+1\right]}{(6+2j)(6-2j)\text{π} }(j=3,0)+\\ \dfrac{4{A}_{2}\left[\mathrm{cos}(\dfrac{5j\text{π} }{3})+\mathrm{cos}\left(\dfrac{2j\text{π} }{3}\right)\right]}{(2+2j)(2-2j)\text{π} }(j=1,-\dfrac{\sqrt{3} }{4})\end{array}$	$\begin{array}{l}\dfrac{12{A}_{1}\left[\mathrm{sin}\left(\dfrac{j\text{π} }{3}\right)\right]}{(6+2j)(6-2j)\text{π} }(j=3,\dfrac{1}{6})+\\ \dfrac{4{A}_{2}\left[\mathrm{sin}(\dfrac{5jn\text{π} }{3})+\mathrm{sin}\left(\dfrac{2j\text{π} }{3}\right)\right]}{(2+2j)(2-2j)\text{π} }(j=1,-\dfrac{1}{4})\end{array}$
天津年周期负荷：${q}_{\text{TJ} }(t)={A}_{1}\mathrm{sin}\left(\dfrac{\text{π} }{4}{t}_{\text{m} }\right)\left(0\leqslant {t}_{\text{m} }\leqslant 4\;\text{月}\right)+{A}_{2}\mathrm{sin}\left(\dfrac{\text{π} }{\text{4} }{t}_{\text{m} }-\dfrac{\text{π} }{2}\right)(6\leqslant {t}_{\text{m} }\leqslant 10\;\text{月})$ 哈尔滨年周期负荷：${q}_{\text{HB} }(t)={A}_{1}\mathrm{sin}\left(\dfrac{\text{π} }{2}{t}_{\text{m} }\right)(0\leqslant {t}_{\text{m} }\leqslant 2\;\text{月})+{A}_{2}\mathrm{sin}\left(\dfrac{\text{π} }{6}{t}_{\text{m} }+\dfrac{\text{π} }{\text{3} }\right)(4\leqslant {t}_{\text{m} }\leqslant 10\;\text{月})$ 傅里叶展开公式：$q\left( t \right) \approx \displaystyle\sum\limits_{j = 0}^n { {q_j}(t) = \;} {a_0} + \displaystyle\sum\limits_{j = 1}^n {\left( { {a_j}\cos \dfrac{ {j{\text{π } } } }{L}t + {b_j}\sin \dfrac{ {j{\text{π } } } }{L}t} \right)}$

下载: 导出CSV 

| 显示表格

式（9）表明，地层温度变化的解析解由三部分组成。为研究各部分对地层温度的影响，本文分别计算了哈尔滨市（式14）与天津市（式13）的解析解的三部分，结果如图7所示。可以看出，在地埋管长期运行条件下，θ_A2对于地层温度的影响有限，可以忽略不计，因此可认为周期性负荷引起的地层温度波动主要由θ_A1与θ_B组成，即：$\theta = {\theta _{{\text{A1}}}} + {\theta _{\text{B}}}$。其中，θ_A1由式（6）可得，θ_B可近似为式（15）：

图  7  不同地区地埋管不同半径处地层温度总周期延米功率的解析解

Figure  7.  Temperature changes at different radii of BHEs in different parts of analytical solution in Harbin and Tianjin

下载: 全尺寸图片 幻灯片

其中：

当$r/{R_{\rm{w}}} \gg 1{\kern 1pt}$时，无线长线热源与无线长柱热源的解高度吻合，θ_B可进一步近似为：

壁面温升可以近似为^[8]：

因此，耦合实际建筑物供热/冷负荷年周期性的地埋管换热器周围地层温度场的解析解可由式（9）简化为式（16）：

3.   结论

（1）提出了耦合实际建筑物供热/冷负荷的地埋管换热器周围地层温度场的解析解，可以快速求解地埋管换热器周围地层温度场的变化，较以往解析解计算精度高，且与数值解相比具有计算快捷、简单方便的特点，适宜工程实际的推广应用。

（2）与数值解的对比验证结果表明：本文提出的方法虽未考虑地层初始条件对温度场的影响，但仅影响前3 年的计算精度；随着时间的推移，解析解和数值解的偏差逐渐减小，3 a以后偏差率基本保持在1%以内，工程上可以忽略不计。

（3）地层温度的变化主要由年周期冷/热负荷函数控制，日周期负荷对于地层温度的影响可忽略。

（4）本文给出了中国3个地理典型地区单根地埋管运行工况下的地层地温场算例。通过修改冷/热负荷函数波幅，可以获取适合不同建筑物类型的地埋管周围地层的温度变化规律。若用户有多年实测的冷热负荷数据，也可通过实测数据得到的实际冷热负荷函数计算地埋管周围的地温场。

（5）结合线性叠加原理，本文所提出的解析解也可用于分析地埋管群的地层温度场变化。