中深层地热井换热数学模型研究

Research of Mathematical Models for Heat Exchange in Medium-Deep Geothermal Wells

摘要: 本文聚焦中深层地热井换热数学模型的研究，在传热过程分析的基础上做出合理假设以简化模型，建立了井孔内和井孔外的传热数学模型，其中井孔内模型考虑不同流体流动方式，井孔外模型则提出有限长圆柱热源模型，以更准确地反映地热井传热实际，为地热井传热性能预测和优化设计提供了理论支撑。

Abstract: This paper focuses on the research of heat exchange mathematical models for Medium-deep geothermal wells. Based on the analysis of the heat transfer process, reasonable assumptions are made to simplify the model. Heat transfer mathematical models for both inside and outside the wellbore are established. The model for inside the wellbore considers different fluid flow patterns, while the model for outside the wellbore proposes a finite-length cylindrical heat source model to more accurately reflect the actual heat transfer in geothermal wells. This study provides theoretical support for the prediction and optimization design of geothermal well heat transfer performance.

1. 引言

随着全球能源需求的不断增加，可再生能源的开发利用日益受到重视。地热资源作为一种清洁、可再生的能源，具有巨大的开发潜力。中深层地热井换热技术作为地热能利用的重要方式，其换热效率的提升对于地源热泵系统的性能优化具有重要意义[1]-[4]。建立准确的中深层地热井换热数学模型，对于深入研究地热井的传热机理[5]-[9]、优化地热井设计[10]-[12]、提高地源热泵系统的整体效率具有重要的理论和实用价值。

2. 中深层地热井传热过程分析与简化

2.1. 地热井传热过程

地热传热主要分为两部分：一是同心管换热器内部结构的功能性传热及其与地热井围岩之间的传热过程；二是地下岩土热储的传热过程。岩土热储内部的传热包括固体岩石骨架的热传导、孔隙中流体的热传导和对流换热等。

中深层地热井下换热器系统主要由管内换热流体、井下同轴套管、钻孔回填材料以及地热井外的岩土构成[13] [14]。该系统的传热过程主要有：(1) 井孔外部岩土的导热过程、(2) 回填材料和井孔壁面的换热过程、(3) 回填材料内部导热过程、(4) 同轴套管外管壁面与回填材料的换热过程、(5) 同轴套管外管管壁的导热过程、(6) 同轴套管环腔内流体与外管内壁的对流换热过程、(7) 同轴套管环腔内流体与内管外壁的对流换热过程、(8) 同轴套管内管管壁的导热过程、(9) 同轴套管内管流体与内管内壁的对流换热过程。其中，地热井下同轴换热器外管壁面与回填材料和回填材料与井孔壁面之间的接触间隙可视为无限小，且温度一致，因此该传热过程可忽略不计。

2.2. 传热过程简化假设

为了简化地热井传热过程的分析，本文做出了一系列合理假设，包括：1、中深层同轴换热器周围的土壤和岩石被认为是一个或几个水平层的均质介质，其热物性不随温度变化；2、换热器的各传热构件(绝热内管、套管)、封固水泥以及岩土微元体内认为是匀质且具有各向同性，传热介质的导热系数等热性质稳定，不随温度变化；3、忽略了大气的温度波动及其对土壤顶层的影响，因此地表以上的气温保持恒定；4、地下换热区域被认为是半无限的；5、计算域内的地热通量(地热密度)均匀；6、同心管换热器水平截面流体的温度、速度沿切向相等；7、采热循环工质为液态纯水，在其工作环境(压力0~20 MPa，温度20℃~60℃)中不被蒸发；8、岩石含水饱和度为1；9、热储及其蕴含流体的黏性和密度与同等条件下实验室测定结果相同；10、地下水的渗流方向可被查明，且仅存在唯一的流动方向；11、忽略封固水泥与外管、热储岩土接触面的余隙热阻；12、忽略流体的粘性摩擦；13、循环工质主要通过热对流进行能量传递，忽略其轴向的热传导；14、渗流场中仅考虑单一组分“水及水蒸气”。这些假设不仅不会降低模型的计算精度，还能简化模型求解，提高计算速度。

3. 地热井孔内传热模型

3.1. 物理模型建立

根据生产资料边界及描述，建立了中深层地热井传热的物理模型。坐标原点设置在内管顶部(图1)，整个传热区域以井孔壁面为界分为井孔内和井孔外两部分进行传热分析。

Figure 1. Physical heat transfer model for medium-deep geothermal wells

图1. 中深层地热井传热物理模型

3.2. 传热过程的热阻分析

中深层地热井换热器中的流体与地热井周围岩土的传热是一个地热井内外耦合的传热过程，与地下岩土相比，井孔内的热容量可忽略不计，所以井孔内流体与套管间的传热过程可以迅速达到平稳阶段，因此可以将井孔内的传热视为稳态传热，井孔外的传热可视为一个线热源或者是柱热源在无限大均匀介质中的非稳态传热过程。

整个传热过程中存在的传热热阻分别为环腔中的流体与地热井壁面的传热热阻和环腔内流体与内管流体之间的传热热阻。将环腔中流体与地热井壁面的传热热阻记为$R_1$ ，将环腔中流体与内管流体间的传热热阻记为$R_2$ 。根据传统传热学原理，可以得到$R_1$ 和$R_2$ 的表达式分别为：

$R_1=\frac{1}{2\pi k_b}\ln \left(\frac{r_b}{r_{1o}}\right)+\frac{1}{2\pi r_{1i}{h}_1}+\frac{1}{2\pi k_{p1}}\ln \left(\frac{r_{10}}{r_{1i}}\right)$ 式(1)

$R_2=\frac{1}{2\pi r_{2i}{h}_2}+\frac{1}{2\pi k_{p2}}\ln \left(\frac{r_{2o}}{r_{2i}}\right)+\frac{1}{2\pi r_{2o}{h}_1}$ 式(2)

式中：$r_b$ ——地热井孔半径，m；$r_{1i}、r_{1o}$ ——分别为外套管的内、外半径，m；$r_{\text{2}i}、r_{\text{2}o}$ ——分别为内套管的内、外半径，m；$k_b、k_{p1}、k_{p\text{2}}$ ——分别为地热井孔的回填材料、外套管和内管的导热系数，W/(m∙K)；$h_1、h_2$ ——分别为环腔流体与环腔管壁、内管流体与内管内壁的对流换热系数，W/(m²∙K)。

影响对流换热系数的因素有很多[15]-[18]，由于流动动力的不同、流动状态的区别、流体是否存在相变以及换热表面的几何形状的差别构成了多种类型的对流传热现象。在地热井内传热过程中所涉及到的对流换热系数，可以根据经典的传热学关联式计算得到。表面对流换热系数可表示为：

$h=f\left(u_f,l,{\rho }_f,v_f,k_f,c_f\right)=\frac{\text{Nu}\cdot k_f}{l}$ 式(3)

式中：$l$ ——换热表面的特征长度，m；$k_f$ ——换热流体的导热系数，W/(m²∙K)。

式(3)中Nu为努赛尔数，是一个无量纲量。其物理意义是传热壁面上流体的无量纲温度梯度。应用Dittus-Boelter公式：

$\text{Nu}=0.023{\mathrm{Re}}^{0.8}{\mathrm{Pr}}_f^n$ 式(4)

式中：$\mathrm{Re}$ ——流体的雷诺数，无量纲；$P_{rf}$ ——流体的定性温度下的普朗特数，无量纲。

当$\mathrm{Re}={10}^4~1.2\times {10}^5$ ，$\text{Pr}=0.7~120,l/d\ge 60$ 。当加热流体时，$n=0.4$ ，冷却流体时，$n=0.3$ 。

根据式(3)可以计算出内管和环管内的对流换热系数，进而求得环腔中流体与地热井壁面的传热热阻$R_1$ 以及环腔中流体与内管流体间的传热热阻$R_2$ 。

3.3. 井孔内传热数学模型的建立

基于热平衡理论和热传导微分方程，本文建立了单个地热井下换热器的传热数学模型。通过分析换热器内的具体传热过程，研究了地热井下换热器系统的传热机理以及影响换热效率的各个因素。

模型假设井孔壁面在轴向的温度均匀、钻孔内部的传热视为稳态过程等。根据循环流体的不同流动方式(外进内出和内进外出)，分别建立传热数学模型。取地热井下同轴换热器在z方向上的一个微元段dz进行分析。管内传热过程可看作为稳态，所以内管中微元段流体通过内管壁的传热量为：

$\delta Q_2=\frac{T_{f2}-T_{f1}}{R_2}{d}_z$ 式(5)

式中：$T_{f1}$ ——外管流体微元段流体的温度，℃；$T_{f2}$ ——内管流体微元段流体的温度，℃。

换热器环管中流体与内管管壁、外管管壁、回填材料以及地下岩土微元段之间的传热量为：

$\delta Q_1=\frac{T_{f1}-T_b}{R_1}{d}_z$ 式(6)

式中：$T_b$ ——地热井壁面温度，℃。

根据能量守恒定律可得流体流动方式为外进内出时流体的能量平衡方程为：

$M{c}_f\frac{\text{d}{T}_{f2}\left(z\right)}{\text{d}z}=\frac{T_{f2}\left(z\right)-T_{f1}\left(z\right)}{R_2}$ 式(7)

通过系列拉普拉斯变换的方法求解，可计算出地热井下同轴换热器内管流体出口的解析表达式为：

${\Theta }_f^{*}={\Theta }_2\left(0\right)=\frac{f_1\left(1\right)+f_2\left(1\right)}{f_3\left(1\right)-f_2\left(1\right)}=\frac{2\beta k_1^{\ast }\cosh \left(\beta \right)-\sinh \left(\beta \right)}{2\beta R_1^{\ast }\cosh \left(\beta \right)+\sinh \left(\beta \right)}$ 式(8)

写为有因次的形式，可得到内管流体出口温度的解析表达式为：

$T_{f2}^{*}=T_b+\left({T^{\prime }}_f-T_b\right)\frac{2\beta R_1^{\ast }\cosh \left(\beta \right)-\sinh \left(\beta \right)}{2\beta R_1^{\ast }\cosh \left(\beta \right)-\sinh \left(\beta \right)}$ 式(9)

考虑到循环流体入口处的无量纲温度为1，可以得到地热井下同轴换热器中流体流动方向为外进内出的布置方式时，流体工质的无量纲温度沿井孔深度方向上分布的解析表达式为：

${\Theta }_2\left(Z\right)=\exp \left(-\frac{Z}{2R_2{}^{\ast }}\right)\left[-f_2\left(Z\right)+f_3\left(Z\right)\cdot \frac{2\beta R_1{}^{\ast }\cosh \left(\beta \right)-\sinh \left(\beta \right)}{2\beta R_1{}^{\ast }\cosh \left(\beta \right)+\sinh \left(\beta \right)}\right]$ 式(10)

同理当循环流体布置为内管进入、环腔流出时，根据能量守恒可以得到流体的能量平衡方程。

$\begin{array}{l}M{c}_f\frac{\text{d}{T}_{f1}\left(z\right)}{\text{d}z}=\frac{T_{f1}\left(z\right)-T_{f2}\left(z\right)}{R_2}+\frac{T_{f1}\left(z\right)-T_b}{R_1}\\ -M{c}_f\frac{\text{d}{T}_{f2}\left(z\right)}{\text{d}z}=\frac{T_{f2}\left(z\right)-T_{f1}\left(z\right)}{R_2}\end{array}$ 式(11)

最终可以得到地热井下同轴换热器在内进外出的流动方式下流体工质的无量纲温度沿井孔深度方向上分布的解析表达式为：

${\Theta }_2\left(Z\right)=\exp \left(\frac{Z}{2R_1{}^{\ast }}\right)\left[f_1\left(Z\right)+f_2\left(Z\right)\cdot \frac{2\beta R_1{}^{\ast }\cosh \left(\beta \right)-\sinh \left(\beta \right)}{2\beta R_1{}^{\ast }\cosh \left(\beta \right)+\sinh \left(\beta \right)}\right]$ 式(12)

通过上述推导得到了地热井下同轴换热器内流体在外进内出和内进外出两种不同的流体流动布置方式下管内流体温度沿井孔深度方向上分布的解析解。

对比式(10)和式(12)可以看出，地热井下同轴换热器内流体在外进内出和内进外出两种不同的布置方式下，换热器内外管流体的温度沿井孔深度方向上的分布式不同的，另外可以看出，在使用均匀的井孔壁温的假设条件下，地热井下同轴换热器内外管的流体在两种不同的流动方式下得到的出口温度却是相等的，说明两种不同的流动方式不影响地热井下同轴换热器的取热性能。

4. 地热井孔外传热模型

对于地热井孔外传热许多学者都采用一维无线长线热源理论和一维无线长圆柱面热源模型，而这种模型相对来说比较简单，计算速度快。但由于其忽略了井孔的有限深度，计算得到的地下岩土温度分布误差较大，并且当时间趋于无穷大时，温度场不会趋于稳定。因此不能用来讨论长时间的地热井下传热问题。针对有限长线热源模型的弊端，本文提出了有限长圆柱热源模型。

地热井下同轴换热器在地下区域中占有一定的空间，有限长线热源模型中忽略了地热井径向的传热会对计算结果产生一定的影响[15] [16]。本节针对有限长线热源的弊端进行改进，在所建地热井外有限长面热源模型中考虑了径向传热，忽略了地热井下岩土中的渗流情况。该模型更符合中深层地源热泵系统的热源井的传热情况[17] [18]。考虑地面边界条件、井孔的有限深度，同时考虑热源在井孔径向的空间分布情况，对井孔外的热源模型进行以下假设：(1) 井孔周围岩土为半无限大的均匀介质，初始温度分布均匀，且热物性恒定；(2) 地面温度恒等于初始温度T₀；(3) 将热源考虑为在井孔壁圆柱面上的均匀分布；(4) 忽略地下水渗流对换热性能的影响。忽略热对流、热辐射，仅考虑导热，在直角坐标系下的热传导微分方程为：

$\frac{1}{a}\frac{\partial T}{\partial t}=\frac{{\partial }^{2}T}{\partial x^2}+\frac{{\partial }^{2}T}{\partial y^2}+\frac{{\partial }^{2}T}{\partial z^2}$ 式(13)

式(13)中，T为温度，t为时间，a为热扩散率。结合相应的边界条件可以确定地源热泵系统作用下的温度场。

在圆柱面热源模型中，对于地面恒温的边界条件，同有限长线热源模型一样采用虚拟热汇的方式来处理，假设在地面上存在一个热物性和地面以下岩土相同的半无限大的均匀介质，虚拟介质中存在与实际热源对称、强度相等的热汇。在t'时刻开始，圆柱面热源以强度为$q_s$ 开始放热，地面上虚拟热汇强度为−$q_s$ ，如图2所示。在热源和热汇的共同作用下，实际介质内的温度场发生改变。则t时刻井孔外某一点的实际温升为Ta。

Figure 2. Finite-length cylindrical surface heat source model

图2. 有限长圆柱面热源模型示意图

在t时刻h处截面上的微元体$\text{d}\theta \text{d}h$ 至圆柱坐标中计算的点$M\left(R,z\right)$ 的水平距离为：

$d=\sqrt{{\left(R-r_b\cos \theta \right)}^2-{\left(r_b\sin \theta \right)}^2}=\sqrt{R^2+r_b-2R{r}_b\cos \theta }$ 式(14)

式中：$R$ ——计算点至井孔中心的水平距离，m；$r_b$ ——热源井半径，m。

考虑瞬时点热源函数，可得t时刻点$M\left(R,z\right)$ 在微元体$\text{d}\theta \text{d}h$ 作用下的温度响应：

$\text{d}{T}_a=\frac{1}{8{\left[\pi \sqrt{\left(at-t^{\prime }\right)}\right]}^3}{e}^A\text{d}\theta \text{d}h$ 式(15)

其中：

$A=-\frac{R^2+r_b^2-2R{r}_b\cos \theta +{\left(h-z\right)}^2}{4a\left(t-t^{\prime }\right)}$ 式(16)

根据叠加原理，式(15)的温度响应$\text{d}{T}_a$ 对$h,\theta ,t$ 进行积分，可以得到：

$T_a=\frac{q_s}{{\rho }_g{c}_g}{\displaystyle {\int }_0^H\text{d}}h{\displaystyle {\int }_0^{2\pi }\text{d}}\theta {\displaystyle {\int }_0^t\exp }\left(-\frac{R^2+r_b^2-2R_{r_b}\cos \theta +{\left(h-z\right)}^2}{4a\left(t-t^{\prime }\right)}\right)\text{d}{t}^{\prime }$ 式(17)

引进变量u，进行变量替换：

$u=\sqrt{\frac{R^2+r_b^2-2R{r}_b\cos \theta +{\left(h-z\right)}^2}{4a\left(t-t^{\prime }\right)}}$

$\text{d}u=\frac{a\sqrt{R^2+r_b^2-2R{r}_b\cos \theta +{\left(h-z\right)}^2}}{4{\left[\sqrt{4a\left(t-t^{\prime }\right)}\right]}^3}\text{d}{t}^{\prime }$ 式(18)

此时公式16化为：

$T_a=\frac{q_s}{p_g{c}_g}{\displaystyle {\int }_0^H\text{d}}h{\displaystyle {\int }_0^{2\pi }\text{d}}\theta {{\displaystyle \int }}_{u_1}^{\infty }\frac{1}{2a{\pi }^2\sqrt{R^2+r_b^2-2R{r}_b\cos \theta +{\left(h-z\right)}^2}}\exp \left(-u^2\right)\text{d}u$ 式(19)

其中$u_1$ 为$t^{\prime }=0$ 时$u$ 的值，引入余误差函数$erfc\left(z\right)$ ：

$erfc\left(z\right)=1-\frac{2}{\sqrt{\pi }}{\displaystyle {\int }_0^z\exp }\left(-u^2\right)\text{d}u$ 式(20)

则地下热源部分引起的温升为：

$T_a=\frac{q_s}{4\pi k_g}{\displaystyle {\int }_0^H\text{d}}h{{\displaystyle \int }}_0^{2\pi }\left\{\frac{erfc\left(\frac{\sqrt{R^2+r_b^2-2R{r}_b\cos \theta +{\left(h-z\right)}^2}}{2\sqrt{at}}\right)}{\sqrt{R^2+r_b^2-2R{r}_b\cos \theta +{\left(h-z\right)}^2}}\right\}\text{d}\theta$ 式(21)

根据对称性可得虚拟热汇部分引起的温升为：

$T_a=-\frac{q_s}{4\pi k_g}{\displaystyle {\int }_0^H\text{d}}h{{\displaystyle \int }}_0^{2\pi }\left\{\frac{erfc\left(\frac{\sqrt{R^2+r_b^2-2R{r}_b\cos \theta +{\left(h+z\right)}^2}}{2\sqrt{at}}\right)}{\sqrt{R^2+r_b^2-2R{r}_b\cos \theta +{\left(h+z\right)}^2}}\right\}\text{d}\theta$ 式(22)

叠加热源和热汇两部分温升的影响，得到有限长圆柱面热源模型：

$T_a=\frac{q_s}{4\pi k_g}{\displaystyle {\int }_0^H\text{d}}h{{\displaystyle \int }}_0^{2\pi }\left\{\frac{erfc\left(\frac{A_1}{2\sqrt{at}}\right)}{A_1}-\frac{erfc\left(\frac{A_2}{2\sqrt{at}}\right)}{A_2}\right\}\text{d}\theta$ 式(23)

式中：

$A_1=\sqrt{R^2+r_b^2-2R{r}_b\cos \theta +{\left(h-z\right)}^2}$ 式(24)

$A_2=\sqrt{R^2+r_b^2-2R_b\cos \theta +{\left(h+z\right)}^2}$ 式(25)

其中$q_s$ 与有限长线热源模型中的线热∙源单位长度热流$q_l$ 有以下的关系：

$q_s=\frac{q_l}{2\pi }$ 式(26)

该模型考虑了地热井径向的传热，更符合中深层地源热泵系统的热源井的传热情况。通过热传导微分方程和相应的边界条件，可以确定地源热泵系统作用下的温度场。

5. 结论

1、建立了中深层地热井传热的物理模型及井孔内、外的传热数学模型，为地热井的传热分析与优化提供了坚实基础。

2、推导出了关键的热阻表达式，并得到了流体温度沿井孔深度方向的解析解，为换热器设计和传热性能预测提供了有力工具。

3、提出的有限长圆柱热源模型提高了地热井传热计算的准确性，为地源热泵系统的性能优化提供了理论支持。

本文成功建立了中深层地热井的传热物理模型及井孔内、外的传热数学模型，为深入理解和优化地热井的传热性能提供了坚实的理论基础。通过推导关键的热阻表达式和流体温度沿井孔深度方向的解析解，为换热器的设计和传热性能预测提供了有力的工具。此外，提出的有限长圆柱热源模型显著提高了地热井传热计算的准确性，为地源热泵系统的性能优化和实际应用提供了更加可靠的依据。综上所述，本研究不仅在地热井换热模型的理论研究上取得了重要进展，还为地热能的高效开发和利用提供了有力的技术支撑。