地埋管换热器是充分利用浅层地热能的技术之一，在节能和环保方面的优势使其得到广泛应用^[1-2]，但其占地面积大、钻孔费用高等缺陷在一定程度上限制了它的发展。能源桩作为一种新的地埋管换热器形式应运而生，它将不同形式的地埋管嵌于建筑物的桩基之中^[3]，通过管内循环工质的流动实现了热量在土壤和流体之间的传递，这大大降低了占地面积和钻孔费用^[4]，从而减少了资金投入，提升了工程经济性。然而，实际工程运用中能源桩的运行会受到诸多外界因素的干扰，地下水渗流就是其中一个非常重要的因素。地下水渗流不仅提高了能源桩的换热性能，也削弱了桩体位移和侧摩阻力的变化幅度等力学性能, 但目前对其影响规律尚不完全清楚。

Zhang等^[5-6]建立了考虑地下水流动的三维能源桩模型，得到了热传导和地下水渗流联合作用下温度响应分析解，表明地下水渗流可以有效提升能源桩的换热性能。Park等^[7]对目前的固体圆柱组合热源能源桩模型进行了适用性评价，发现当钻孔直径与钻孔长度之比较大时，模拟出的解析解高估了地下水流动对能源桩换热性能的影响。Wang等^[8]用热-多孔-力学全耦合有限元方法分析了能源桩与桩周土壤的相互作用，并借助连续介质力学基本定律，推导出了非线性控制方程来解释多孔介质力学响应、气体流动、水分流动和热-流耦合过程。Han等^[9]借助COMSOL软件建立了有地下水渗流和桩-土相互作用的能源桩三维耦合模型，并模拟计算了渗流流速对能源桩热-力耦合特性的影响，发现较快的渗流流速有利于热量在土壤中的传递，但对能源桩的桩轴剪切应力大小的影响微乎其微。Yang等^[10]数值模拟探究了地下水渗流的有无以及渗流流速对能源桩桩顶位移的影响，发现地下水渗流流速越大，水带走的包含在桩体和土壤中的热量越多，桩顶位移变化越小。Lou等^[11]针对居住建筑中常见的4种能源桩进行建模研究，发现地源热泵机组的性能系数(coefficient of performance, COP)随着渗流流速的增大而增大，且上下游桩之间的温差也与渗流流速呈正相关。You等^[12-13]利用三维有限元模型模拟研究了不同参数对能源桩桩群换热性能的影响，结果表明：地下水渗流流速对桩群换热性能影响最大，渗流流速越快，越有利于热量向桩周土壤的“迁移”，从而缓解土壤因取放热不平衡而导致的温度失衡状况。

综上所述，地下水渗流不仅可以使堆积在土壤中的热量充分扩散，还可以降低桩体温度的变化幅度，进而影响能源桩的力学性能。因此，能源桩在有地下水渗流的情况下运行过程是复杂的热-流-力耦合过程。然而，现有研究没有充分揭示地下水渗流对能源桩热-力耦合特性的作用机理。为此，本文建立了较为完善的数学模型，对有无地下水渗流工况下能源桩的运行情况进行数值模拟，探讨了地下水渗流对能源桩换热和热-力耦合作用的影响，从而为日后能源桩在实际工程中的运用提供借鉴。

1 计算模型 1.1 物理模型

能源桩与桩周土壤间的热-流-力耦合是一个非常复杂的过程，为简化分析作出如下假设:

1) 能源桩为均匀一致各向同性材料，物性参数保持不变, 忽略水分迁移导致的热量传递。

2) 不考虑深度方向传热，且桩体和桩周土壤温度初始值相同，土壤远边界温度不变。

3) 桩体周围土壤是多孔且充满水分的介质，热量在其中的传递方式包括热传导和热对流。

4) 埋管从管壁到桩周土壤的传热过程不考虑接触热阻。

5) 渗流流速和方向保持不变，且不考虑渗流流速沿深度方向的变化。

6) 桩体为弹性变形，桩周土壤发生弹塑性变形，桩土间摩擦因数保持不变。

基于以上假设，建立如图 1所示的内置并联双U形埋管能源桩的换热及热-力耦合物理模型。模型参数见表 1。

1.2 数学模型 1.2.1 U形管内流体控制方程

U形管内流体处于湍流形式，使用k-ε湍流模型^[14]能更加贴切地描述湍流过程:

$ \begin{gathered} \frac{\partial}{\partial t}\left(\rho_{\mathrm{f}} k\right)+\frac{\partial}{\partial x_{j}}\left(\rho_{\mathrm{f}} k u_{j}\right)= \\ \frac{\partial}{\partial x_{j}}\left[\left(\mu+\frac{\mu_{\mathrm{t}}}{\sigma_{k}}\right) \frac{\partial k}{\partial x_{j}}\right]+G_{k}-\rho_{\mathrm{f}} \varepsilon. \end{gathered} $

(1)

$ \begin{gathered} \frac{\partial \rho}{\partial t}\left(\rho_{\mathrm{f}} \varepsilon\right)+\frac{\partial \rho}{\partial x_{j}}\left(\rho_{\mathrm{f}} \varepsilon u_{j}\right)=\frac{\partial}{\partial x_{j}}\left[\left(\mu+\frac{\mu_{\mathrm{t}}}{\sigma_{\varepsilon}}\right) \frac{\partial \varepsilon}{\partial x_{j}}\right]+ \\ \rho_{\mathrm{f}} C_{1} S \varepsilon-\rho_{\mathrm{f}} C_{2} \frac{\varepsilon^{2}}{k+\sqrt{\nu \varepsilon}} . \end{gathered} $

(2)

式中: t表示时间，s；ρ_f为桩埋管内流体密度，kg/m³；k为湍流动能，J；u_j为x_j坐标上的速度矢量分量，m²/s；μ为流体的动力黏度，Pa·s；ν为运动黏度，m²/s；μ_t为湍流黏度，Pa·s；ε为湍流耗散率；G_k是由平均速度梯度产生的k的迭代项；湍动能k和ε的Prandtl数σ_k和σ_ε分别取1.0和1.2；C₁和C₂是经验常数。

1.2.2 能源桩桩体传热控制方程

能源桩桩体内部传递热量的方式为热传导，控制微分方程可写成

$ \rho_{\mathrm{c}} C_{\mathrm{c}} \frac{\partial T_{\mathrm{c}}}{\partial t}=\nabla \cdot\left(\lambda_{c} \nabla T_{\mathrm{c}}\right)+q_{\mathrm{c}}. $

(3)

式中: C_c为桩体比热容，J/(kg·K)；T_c为桩体温度，K；ρ_c为桩体材料密度，kg/m³；q_c为桩体体积热源，W/m³；λ_c为桩体材料导热系数，W/(m·K)。

1.2.3 能源桩周围土壤渗流控制方程^[15]

能量方程：

$ \begin{gathered} \frac{\partial}{\partial t}\left(\gamma\rho_{\mathrm{sf}} E_{\mathrm{sf}}+(1-\gamma) \rho_{\mathrm{s}} E_{\mathrm{s}}\right)+\nabla \cdot\left(\boldsymbol{v}_{\mathrm{sf}}\left(\rho_{\mathrm{sf}} E_{\mathrm{sf}}+p\right)\right)= \\ S_{\mathrm{sf}}+\nabla \cdot\left[\lambda_{\mathrm{eff}} \nabla T-\left(\sum\limits_{i} h_{i} \boldsymbol{J}_{i}\right)+\left(\bar{\tau} \cdot \boldsymbol{v}_{\mathrm{sf}}\right)\right]. \end{gathered} $

(4)

式中: γ为土壤孔隙率；ρ_s和ρ_sf分别为土壤中固体密度和渗流水密度，kg/m³；E_s和E_sf分别为多孔介质中固体和渗流流体单位质量的内能，J；p为压强，Pa；S_sf为多孔介质中渗流流体的熵，J/K；h_i为渗流流体的焓，J；J_i为渗流流体扩散通量，kg/(m²·s)；v_sf为渗流流速，m/s；λ_eff为土壤的有效热传导系数，W/(m·K)，

$ \lambda_{\text {eff }}=\gamma \lambda_{\mathrm{sf}}+(1-\gamma) \lambda_{\mathrm{s}}. $

(5)

式中: λ_sf和λ_s分别为渗流水和土壤中固体物质的导热系数，W/(m·K)。

动量方程:

由于地下水渗流流速非常缓慢，因此将其视为层流流动，采用经过修正的动量方程描述，

$ \frac{\partial}{\partial t}\left(\rho_{\mathrm{sf}} u\right)+\frac{\partial}{\partial t}\left(\rho_{\mathrm{sf}} u u\right)=\mu \frac{\partial^{2} u}{\partial x^{2}}+S_{i}. $

(6)

S_i是附加源项，由黏性损失项和惯性损失项两部分组成。假定土壤为均匀多孔介质，则动量源项可以简化为

$ S_{i}=\frac{\mu}{\alpha} v_{i}+C_{3} \frac{1}{2} \rho\left|\boldsymbol{v}_{\mathrm{sf}}\right| \boldsymbol{v}_{i}. $

(7)

式中: α是土壤渗透率，C₃是惯性阻力因子。

连续性方程:

$ \frac{\partial\left(\gamma \rho_{\mathrm{sf}}\right)}{\partial t}+\nabla \cdot\left(\gamma \rho_{\mathrm{sf}} \boldsymbol{v}_{\mathrm{sf}}\right)=0. $

(8)

1.2.4 桩土本构模型

由于桩体周围土壤的强度远远小于能源桩本身，随着时间的推移，土壤的强度慢慢超过其弹性极限，进入弹塑性阶段，而能源桩依旧保持弹性变形的状态。因此，对于能源桩周围土壤采用Mohr-Coulomb弹塑性模型^[16]，见式(9)；而对于能源桩采用弹性本构模型^[17]，见式(10)和(11)。

$ \boldsymbol{A}=\boldsymbol{D}\left(\varepsilon_{\mathrm{obs}}+\alpha_{\mathrm{c}} \Delta T\right). $

(9)

$ R_{\mathrm{mc}} \sigma_{q}-\sigma_{p} \tan \varphi-c=0. $

(10)

$ \begin{gathered} R_{\mathrm{mc}}=\frac{1}{\sqrt{3}} \cos \varphi \sin \left(\varTheta+\frac{\pi}{3}\right)+ \\ \frac{1}{3} \cos \left(\varTheta+\frac{{\rm{ \mathsf{ π}}}}{3}\right) \tan \varphi . \end{gathered} $

(11)

式中: A为应力矩阵；D为弹性矩阵，与弹性模量和Poisson比相关联；α_c为桩体材料的热膨胀系数，1/℃；ε_obs为桩体因为温升产生的观测应变；ΔT为桩体温度的变化值，℃; R_mc为π平面上屈服面形状的一个度量；φ和Θ分别为桩周土壤的内摩擦角和偏极角，(°)；c表示能源桩周围土壤黏聚力，σ_q为偏应力，σ_p为平均有效应力，kPa。

1.3 定解条件 1.3.1 温度场

1) 初始条件。

$ \begin{aligned} &T_{\mathrm{f}}(z, t)=T_{\mathrm{c}}(x, y, z, t)= \\ &T_{\mathrm{s}}(x, y, z, t)=T_{0}(t=0) . \end{aligned} $

(12)

2) 边界条件。

土壤远边界为

$ \left.T_{\mathrm{s}}(x, y, z, t)\right|_{x, y=\pm 7.5}=T_{0} . $

(13)

上部边界条件为

$ \left.\frac{\partial T}{\partial z}\right|_{z=0}=0. $

(14)

下部边界条件为

$ \left.\frac{\partial T}{\partial z}\right|_{z=25}=0. $

(15)

式中: T_f(z, t)为t时刻能源桩埋管z处的温度，℃；T_c(x, y, z, t)是t时刻(x，y，z)处的桩基温度，℃；T_s(x, y, z, t)是t时刻(x，y，z)处的土壤温度，℃；T₀为土壤初温，℃。

1.3.2 渗流场

地下水渗流入口采用速度边界，规定渗流方向和流速为定值(x轴方向为正)，进口温度为17.5 ℃，将出口设置为压力出口，从而防止回流效应的干扰。

1.3.3 结构场

对于结构场的控制约束主要为: 1) 土体的底部设为固定端约束，即x、y、z 3个方向的位移均为零。2) 土体顶部不进行约束，作自由面处理。3) 初始地应力会造成桩基的承载力下降，因此在数值模型中应进行考虑，即z方向自由移动、x和y两个方向设置位移为零。4) 桩顶施加向下的力荷载，同时桩体受热(受冷)可以自由膨胀(收缩)，即桩体在x、y、z 3个方向上均为自由移动状态。

1.4 网格划分及其独立性验证

利用Gambit软件进行能源桩模型的网格划分^[15]。由于深度方向温度变化小，U形管深度方向采用较为稀疏的网格，而在U形管弯管处流场急剧变化，因此选择较为紧密的网格。为保证各接触面处模拟结果的可靠性，接触面处网格疏密程度保持一致。网格划分如图 2所示。

为了验证所划分的网格的有效性，进行相关独立性验证。选取网格数量为70×10⁴、116×10⁴、225×10⁴和342×10⁴的桩体数值模型，采用不同网格数的模型在同样的参数条件下不间歇运行24 h。图 3为出口温度随网格数量的变化。可以发现，4种网格数量计算出的埋管出口温度相差甚微。因此，综合计算时间和精度，网格数定为116×10⁴。同时，选取6种不同时间步长，验证其对运行结果的影响。图 4给出了出口温度随时间步长的变化。可以看出，时间步长大于600 s时，随着时间步长的增加，出口温度下降显著，但时间步长小于600 s时，计算结果与600 s时模拟出的结果接近，因此选择600 s作为模拟时间步长。

2 模拟结果与分析

采用Fluent和Abaqus软件对渗流场下能源桩的传热及热-力耦合特性进行求解，模拟运行时间为14 d，桩内埋管入口水温为35 ℃，流速为0.6 m/s，桩体及周围土壤的初始温度设定为17.5 ℃，地下水渗流流速为60 m/a，其余相关参数设置见表 2^[15]。表 2中混凝土相关热力学参数根据型号为C30混凝土对应的参数进行选取，土壤相关热力学参数根据黏土对应的参数进行选取。

2.1 换热性能 2.1.1 换热量

由图 5和6可知，在地下水渗流的影响下，能源桩换热量经过一段时间逐渐趋于稳定，而无渗流工况下，桩体换热量持续减少且换热量明显小于有地下水渗流工况；随着时间推移，两者换热量的差值愈来愈大，如第4天两者的日换热量差值为9.17 MJ，而到运行的最后一天，差值增加至30.97 MJ。这主要是由于地下水的存在使热量随着渗流方向迁移，从而降低了桩土温度变化幅度，使得能源桩换热性能增强。

2.1.2 热影响范围

图 7给出了有无地下水渗流工况下运行14 d后桩深7.5 m处桩周土壤温度分布云图。可以看出，无地下水渗流时，桩体向周围土壤对称散热，热影响范围约为1.88 m；有地下水渗流时，温度场沿地下水流动方向偏移，其中地下水上游、下游和垂直方向热影响范围分别是1.06、4.03、1.85 m。进一步对比分析可以看出，在地下水流动的影响下，地下水上游和下游的热影响范围是无地下水时的43.62%和114.36%。这表明渗流场的存在，使得桩群中前排桩排入土壤中的热/冷量会沿渗流方向传递，进而影响后排桩的散热。

2.1.3 桩体温升

图 8给出了有无地下水渗流工况下桩截面平均温升沿深度方向的变化。显然，渗流导致桩截面平均温升降低，如运行结束时有无渗流工况下分别为12.95和14.25 ℃。渗流导致桩截面平均温升降低了9.12%。进一步分析图 8还可以发现，随着运行时间增加，渗流工况下的桩截面平均温升变化幅度明显低于无渗流工况，如在有地下水渗流工况下，运行结束时相较于运行3 d时的桩截面平均温度仅相差0.60 ℃，而在无地下水渗流工况下为1.50 ℃。这表明地下水渗流不仅降低桩体温升，也缩短了桩体温度场达到稳定状态所需时间。

2.2 桩体位移

图 9给出有无地下水渗流工况下单独建筑荷载和热-力耦合荷载下桩体位移沿深度的变化曲线。可以看出，由于桩体受到周围土壤向上的侧摩阻力束缚，且阻力随着深度的增加而增大，导致在单独建筑荷载下能源桩沿深度方向位移量逐渐变小。但在热-力耦合作用下，桩体受热向桩两端膨胀，桩体中下部存在位移零点，且有地下水渗流工况下桩体位移明显小于无地下水渗流工况。

为更直观地反映渗流作用对桩体位移的影响，图 10给出了有无渗流工况下热-力耦合荷载下桩体位移增量(位移增量为热-力耦合荷载与单独建筑荷载下桩体位移的差值)的变化曲线。由图 10可知，运行14 d后，无渗流和有渗流工况下桩顶位移增量分别是1.04和0.95 mm，显然地下水渗流使得桩顶位移增量减少了5.77%。这主要是由于渗流作用使得桩体温升低于无渗流工况，导致其热膨胀变形也较小。进一步分析可发现，有地下水渗流工况下，能源桩运行14 d相比于运行第3天的位移增量仅增加0.054 mm，而无地下水渗流工况下为0.12 mm，增量变化量减少55.00%。这表明，地下水渗流降低了桩体的温升，进而导致桩体位移增量减少。

2.3 桩体轴力

由图 11可以看出，在单独建筑荷载下，有无地下水渗流工况下桩体轴力都表现出沿深度方向逐渐变小的趋势，其主要原因是在单独力荷载的传递过程中桩受到周围土壤向上侧摩阻力的削弱。进一步可以看出，由于在加热过程中产生的附加温度荷载对桩体产生的是压应力，因此热-力耦合荷载作用下产生的轴力大于单独建筑荷载。分析图 12可以看出，附加温度荷载随着深度先变大后减小，且最大值在6.5 m深处。这是因为桩体上部和下部受热，分别向桩顶和桩端膨胀，这两处受到土壤侧摩阻力的束缚，使得附加温度荷载向两端逐步降低。例如，在有地下水渗流工况下，能源桩运行3 d和运行14 d的最大附加温度荷载分别是63.4和64.7 kN，对应的无地下水渗流工况下分别为64.8和69.8 kN，即运行结束后地下水渗流使最大附加温度荷载减少了7.31%；而且，运行3 d至14 d期间，有地下水渗流工况下的附加温度荷载变化幅度只是无地下水渗流工况下的26.00%。这是因为附加温度荷载与桩体温升呈正相关变化。

2.4 侧摩阻力

图 13给出了单独建筑荷载和热-力耦合荷载下桩体侧摩阻力的变化。可以看出，与受单独建筑荷载相比，热-力耦合荷载下桩体下部侧摩阻力变大，而桩体上部变小。这主要是因为热-力耦合荷载下桩体受热向两端膨胀，桩的上部侧摩阻力向下，而下部的侧摩阻力向上，但单独荷载下桩所受侧摩阻力都向上，所以桩的上部侧摩阻力相互抵消，下部的侧摩阻力互相累加。进一步分析图 13，运行14 d后，有无地下水渗流工况对应的最大侧摩阻力分别为-6.26和-7.11 kPa，有地下水渗流工况下的最大负侧摩阻力比无地下水渗流工况下降低了11.95%。另外，有无地下水渗流工况下在第3 d至第14 d之间最大负侧摩阻力分别增加了0.11和0.74 kPa，这表明在此期间，有地下水渗流工况下负侧摩阻力变化幅度相对于无地下水渗流工况下减少了85.14%。可见，地下水的渗流降低了能源桩受热膨胀的程度，从而减小了侧摩阻力的变化幅度。

3 结论

1) 地下水渗流会带走能源桩释放的热量，降低桩土温度变化幅度，维持换热温差，强化能源桩的换热性能，使桩体温度场尽快稳定。运行结束时，有地下水渗流工况下能源桩日换热量比无地下水渗流工况下多出30.97 MJ。

2) 桩体周围土壤温度场随渗流流动方向发生偏移，这会影响后排桩的传热性能。在有地下水渗流工况下，桩下游的热影响范围是无地下水渗流工况下的114.36%，而桩上游的热影响范围是无地下水渗流工况下的43.62%。

3) 与无地下水渗流工况相比，地下水的存在一定程度上导致桩体温升幅度较小，降低了约9.12%，从而导致桩体位移、桩体轴力和侧摩阻力的变化量相较于无地下水渗流工况偏小，有利于能源桩的长期运行。运行结束后，有渗流和无渗流工况对应的桩顶位移分别是0.95和1.04 mm，7.5 m桩深处的最大附加温度荷载分别是64.7和69.8 kN，桩体上部所受最大负侧摩阻力分别是-6.26和-7.11 kPa。