0.   引言

由于化石能源的大量使用，全球气候变化和能源短缺已经成为21世纪人类面临的最大挑战之一。随着中国“碳达峰、碳中和”目标的确定，为减少对化石能源的消耗和温室气体的排放，各种清洁能源将得到广泛应用。地热能作为一种绿色环保、储量充足、稳定可靠的可再生能源，具有较强的竞争力和广阔的应用前景。

地热能源桩简称能源桩、能量桩或地热桩，依靠地源热泵系统，在能源桩和桩周土体间进行热量交换。在夏季，通过吸收建筑物的热量制冷；在冬季，摄取地下热量用于建筑物供暖。过去20 a的相关研究证明了能源桩在技术和经济上的可行性。利用混凝土良好的储热性能，能源桩可用于建筑物以及道路、隧道等的制热与制冷。相比于传统钻孔地源热泵技术，能源桩将桩基础与热交换管结合，节约了用地面积和钻孔费用，因此该技术已逐渐成为岩土工程研究的热点。

通过原位试验可以研究能源桩在现场条件下的多场耦合行为。Laloui等^[1]和Bourne-Webb等^[2]分别开展了能源桩热-力响应的原位试验，试验场地为砂土和黏土，分别得到了桩体的应力-应变分布及其受力特性。桂树强等^[3]、路宏伟等^[4]也通过原位试验研究了能源桩的结构响应和承载特性。原位试验通常耗时长、成本高，因此缩尺模型试验也常用于研究能源桩的多场耦合特性。Stewart等^[5]、Goode Ⅲ等^[6]、Ng等^[7]对黏土和砂土层内端承型和摩擦型能源桩开展了离心模型试验，获得了能源桩应力-应变分布规律和承载特性。黄旭等^[8]和孔纲强等^[9]利用模型试验研究了温度循环对能源桩热-力行为的影响。考虑到场地和模型试验普遍具有复杂性和不易重复等局限性，一些学者也采用了数值分析方法，Laloui等^[1]、Wang等^[10-11]、Loria等^[12-13]、Saggu等^[14-15]、di Donna等^[16]、Rui等^[17]、王成龙等^[18]、李翔宇等^[19]验证了原位试验和模型试验结果，并进行了参数敏感性分析。

能源桩温度变化导致桩体产生附加应力与变形，桩周土体出现热-水-力耦合现象，显著影响桩土相互作用及桩承载力学特性，并对建筑物的长期稳定和地热能的利用效果产生影响。Laloui等^[1]未考虑温度梯度与达西流的耦合作用，而本文基于热流固三场全耦合理论，同时考虑了流体性质随温度的变化，研究能源桩的热-力耦合行为，并且比较了桩土间绝对粗糙和桩土间库仑摩擦的模拟结果。研究成果可为能源桩设计、施工及使用过程中的多场耦合问题提供参考。

1.   THM耦合理论及其验证

1.1   THM耦合分析简介

多孔介质的热-水-力耦合分析方法已广泛应用于岩土工程的诸多领域，包括地热利用、二氧化碳地质封存、核废料处置、地下管线冻胀、油藏模拟等。

Coussy^[20]从开放系统和不可逆过程热力学出发，建立了饱和多孔材料热-孔隙-弹塑性的一般理论。Booker等^[21]通过将温度场与位移和压力的计算解耦，提出了饱和热弹性体中点热源和球热源固结效应的解析解。Bai等^[22]给出了一维线性准静态弹性饱和体系的解析解，并比较了各种耦合情况和各耦合项的影响。Coussy^[23]从多孔介质力学基本理论出发推导了一维多孔弹性介质在温度荷载作用下的解析解。Olivella等^[24]开发了用于多孔介质中流动和热分析的模拟程序CODE_BRIGHT。Coussy^[23]、De Boer^[25]、Lewis等^[26]开展了饱和与非饱和多孔介质弹性和弹塑性热-孔隙-力学模型的有限元分析。

1.2   THM耦合理论

从动量、质量和能量等3个守恒条件出发，可推导THM三场全耦合的强形式方程组，该方程组给出了孔隙压力、固体位移和热流等在孔隙介质中的演化，详见Wang^[11]、Coussy^[23]和De Boer^[25]。

准静态条件下的动量守恒方程为

式中：$\mathit{\boldsymbol{\sigma }}$为总应力张量，受拉为正；$\rho$为混合物质量密度；$\mathit{\boldsymbol{b}}$为加速度矢量。总应力张量可写为

式中：$\mathit{\boldsymbol{\sigma }} '$为有效应力张量；b为Biot系数；p为孔隙流体压力；$\mathit{\boldsymbol{\delta }}$为Kronecker符号，表示单位二阶张量。

值得一提的是，弹性应变可分为两部分：

式中：${\mathit{\boldsymbol{\varepsilon }} ^\mathit{\boldsymbol{e}}}$为弹性应变张量；E为刚度矩阵；${\beta _{\text{s}}}$为固体的体积热膨胀系数；$\theta$为温度。

质量守恒方程为

式中：${\mathit{\boldsymbol{v}}_{\text{s}}}$为固体骨架的速度矢量；${c_{\text{m}}}$为混合物的压缩系数；${\beta _{\text{m}}}$为混合物的体积热膨胀系数；${\mathit{\boldsymbol{v}}^{\bf{D}}}$为达西速度矢量。其中${c_{\text{m}}}$，${\beta _{\text{m}}}$可由各组分定义，${\mathit{\boldsymbol{v}}^{\bf{D}}}$可由达西定律定义：

式中：n为多孔介质的孔隙率；${c_{\text{f}}}$，${c_{\text{s}}}$分别为流体和固体的压缩系数；${\beta _{\text{f}}}$，${\beta _{\text{s}}}$分别为流体和固体的体积热膨胀系数；k为多孔介质的渗透率；${\mu _{\rm{f}}}$为流体的动力黏度；${\rho _{\text{f}}}$为流体密度；g为重力加速度。

能量守恒方程为

式中：${(\rho C)_{\text{m}}}$为混合物的热容密度；K为固体骨架的体积模量；q为热流矢量；$\rho r$为热源项。其中${(\rho C)_{\text{m}}}$可由各组分定义，q可由傅里叶定律定义：

式中：${C_{\text{f}}}$，${C_{\text{s}}}$分别为流体和固体的比热容；${\lambda _{\text{m}}}$为混合物的平均导热系数；${\lambda _f}$，${\lambda _{\text{s}}}$分别为流体和固体的导热系数。

将式（1），（4），（8）联立可得THM全耦合强形式控制方程组。静力平衡、流体流动和热流的3个方程含有固体位移、孔隙压力和温度3因变量，且因变量间相互关联，因此该耦合是完全的，并非单向的。

另外，Adinolfi等^[27]认为热-水-力耦合中温度场对水力场的影响除了流体的热膨胀，还有流体性质会随温度的变化而变化。本文通过多项式拟合来考虑此效应^[28]：

1.3   THM耦合理论的数值验证

采用多物理场仿真软件COMSOL Multiphysics^[28]建立饱和介质THM耦合有限元模型。为了验证模型的正确性，将COMSOL软件的模拟结果、Wang^[11]的模拟结果与Bai等^[22]的解析解进行对比分析。

该算例为饱和非等温固结问题，不计重力影响，几何模型及边界条件如图 1所示。土柱高度为0.3 m，半径为0.1 m。该模型初始温度场为50℃。力学边界条件为：顶部施加竖向荷载F=-1×10⁷ Pa，底部固定边界，r为0，0.1 m处侧面边界约束其在r方向上的位移。水力边界条件为：顶部为排水边界（p=0），其它边界为不排水边界。热边界条件为：顶部设为恒温边界（θ=0℃），其它边界为绝热边界。在软件中使用了固体力学、达西定律和多孔介质传热3个模块，场变量分别为位移u、水压p和温度θ。非结构三角形网格（646个单元和360个节点）的模拟结果与图 1所示四边形网格（3个单元和8个节点）的模拟结果基本相同。表 1给出了模型的材料参数。

图  1  饱和非等温固结问题的几何模型及边界条件

Figure  1.  Geometric model and boundary conditions for saturated non-isothermal consolidation

下载: 全尺寸图片 幻灯片

表  1  饱和非等温固结问题的材料参数^{[11, 22]}

Table  1.  Material parameters for saturated non-isothermal consolidation ^{[11, 22]}

参数	取值	单位
Lamé第一常数λskel	7.05×107	Pa
Lamé第二常数μskel	4.7×107	Pa
Biot系数b	1.0	1
孔隙率n	0.42	1
渗透率k	1.0×10-14	m2
流体动力黏度${\mu _{\text{f}}}$	0.001	Pa·s
固体密度${\rho _{\text{s}}}$	3696	kg/m3
流体密度${\rho _{\text{f}}}$	1000	kg/m3
固体压缩系数${c_{\text{s}}}$	0	Pa-1
流体压缩系数${c_{\text{f}}}$	0	Pa-1
固体体积热膨胀系数${\beta _{\text{s}}}$	4.95×10-6	K-1
流体体积热膨胀系数${\beta _{\text{f}}}$	6.21×10-4	K-1
固体比热容${C_{\text{s}}}$	703	J/(K·kg)
流体比热容${C_{\text{f}}}$	4180	J/(K·kg)
固体导热系数${\lambda _{\text{s}}}$	1.38	W/(m·K)
流体导热系数${\lambda _{\text{f}}}$	0.6	W/(m·K)

下载: 导出CSV 

| 显示表格

图 2给出了土柱非等温固结过程中温度、孔隙流体压力和顶部位移的响应。土柱底部的温度变化如图 2（a）所示，底部温度在10~100 h期间变化较大，在100 h后基本达到稳定状态。图 2（b）的结果是无量纲孔隙流体压力（p/F）随时间的变化，可以看出100 s时超静孔隙压力基本消散，由于未考虑重力影响，孔隙压力在稳定状态时为0。固结也基本在100 s时完成，顶部沉降稳定值约为18 mm，如图 2（c）所示。通过比较可以看出模拟结果和解析解吻合较好，验证了本文所建立的饱和介质THM全耦合有限元模型有效性。

图  2  本文模拟结果与解析解的对比

Figure  2.  Comparison between simulated results and analytical solutions

下载: 全尺寸图片 幻灯片

2.   能源桩的THM全耦合有限元分析

2.1   模型配置与简化

实际工程中的能源桩通常是群桩的一部分，本文主要考虑单桩，但所建立的计算方法同样适用于群桩的计算分析。能源桩横截面通常为圆形，垂直放置在土中，假设土层和地下水位水平，因此可以建立轴对称模型。桩和岩层假定不透水，各土层渗透率可通过现场试验确定，并假设各向同性渗透。此外，由于本文主要研究能源桩与周围土体的耦合行为，无需对热交换管和热流体进行精细建模，因此将热交换管模拟为线热源或将整个桩体作为热源。

对于桩体本构模型，绝大部分文献[1, 10~13, 15~18, 29, 30]采用了热弹性模型，仅Saggu等^[14]采用了混凝土损伤塑性模型。对于土体本构模型，大部分文献采用了热弹性模型^{[10-11, 14, 18, 29-30]}、Drucker- Prager准则^[1]、Mohr-Coulomb准则^{[12-13, 16]}、剑桥黏土模型^[11]等常用本构模型，也有Saggu等^[15]、Rui等^[17]采用了CASM模型和热弹塑性剑桥黏土模型。土体采用弹塑性模型时塑性区基本只出现在桩顶和桩底附近，对重点关注的桩的应力应变分布影响较小；另外Laloui等^[1]和Bourne-Webb等^[2]的原位试验结果中桩和土也都表现出较明显的弹性响应特征，因此对于桩体和土体本文均假定为热弹性材料，土层力学参数应由三轴试验确定。

桩和土之间的接触假定有绝对粗糙^{[1, 12, 17]}、库仑摩擦^[14-15]和荷载传递模型^[31-32]3种。桩土界面绝对粗糙即桩和土的有限元节点之间没有相对运动，库仑摩擦通过下式来模拟桩土相互作用：

式中：${T_{\text{t}}}$为摩擦力；$\mu$为摩擦系数；${T_{\text{n}}}$为接触压力；${T_{{\text{cohe}}}}$为内聚抗滑阻力；${T_{{\text{t, max}}}}$为最大摩擦力；$\mathit{\boldsymbol{\delta }}$为界面摩擦角。

将底部的竖向位移和侧面的水平位移约束，作为模型的力学边界条件。每次计算前需要进行地应力平衡，即首先计算桩和土在重力作用下产生的初始应力，然后将初始应力作为计算初始值，以保证初始时刻合理的应力状态和零位移。水力边界条件是设定顶部为排水边界，沿对称轴、岩层和桩土界面为不排水边界。地应力平衡时也需要考虑地下水，初始孔隙水压力对应静水压力分布，水位位于地表。需要说明的是，自然条件下模型右侧是可渗透的，但结构-温度荷载的影响范围有限，对模型最右侧几乎没有影响，因此为了提高计算效率，将右侧也设为不排水边界。类似地，热边界条件包括对称轴、底部和右侧的绝热边界条件以及顶部的恒温边界条件。

2.2   离心机模型的有限元分析

本节采用轴对称THM耦合有限元模型模拟饱和介质中能源桩的热-力行为，并与Wang等^[10]的模拟结果进行对比分析。模型假设线弹性和准静态加载过程。该模型考虑了土骨架和水两相，分别满足各向同性可压缩及不可压缩条件。有限元模型采用混合四边形单元模拟耦合效应。其中，位移场采用双二次插值，孔隙压力和温度场则采用双线性插值。

该算例的几何模型、边界条件及网格如图 3所示，包含35个单元和48个节点。根据Wang等^[10]，此模型的原型是饱和粉土中能源桩缩尺离心机试验，离心加速度为24g（g为重力加速度）。模型桩半径为0.025 m，长度为0.38 m；模型箱半径为0.25 m，高度为0.54 m。轴向荷载F=-268 kPa作用于能源桩顶，方向竖直向下，在模拟过程中保持恒定。整个模型的温度均为20℃，初始时刻起将桩体的对称轴设定为50℃，模型顶部的温度在室温下保持为恒定20℃。对位移、水压和温度3个因变量，通过向后欧拉法时间积分，并采用全耦合牛顿迭代法求解。为与Wang等^[10]的数值模型保持一致，本模型桩土界面假定为绝对粗糙。混凝土桩和土的材料参数如表 2所示。

图  3  能源桩离心机模型的几何模型、边界条件及网格

Figure  3.  Geometric model, boundary conditions and mesh for centrifuge model for energy piles

下载: 全尺寸图片 幻灯片

表  2  离心机模型中桩和土的材料参数^[10]

Table  2.  Material parameters for energy piles and soils in centrifuge tests ^[10]

参数	取值	单位
桩的弹性模量Ec	3.452×1010	Pa
桩的泊松比vc	0.18	1
土的弹性模量Esoil	1.222×104	Pa
土的泊松比vsoil	0.3	1
饱和水力传导率ksoil	5×10-7	m/s
桩的密度${\rho _{\text{c}}}$	2564	kg/m3
土的密度${\rho _{{\text{soil}}}}$	1916	kg/m3
水的密度${\rho _{\text{w}}}$	ρf(θ)	kg/m3
桩的线热膨胀系数${\alpha _{\text{c}}}$	8×10-6	K-1
土的线热膨胀系数${\alpha _{{\text{soil}}}}$	1.17×10-5	K-1
桩的比热容${C_{\text{c}}}$	855	J/(K·kg)
土的比热容${C_{{\text{soil}}}}$	1190	J/(K·kg)
水的比热容${C_{\text{w}}}$	Cf(θ)	J/(K·kg)
桩的导热系数${\lambda _{\text{c}}}$	1.92	W/(m·K)
土的导热系数${\lambda _{{\text{soil}}}}$	1.2	W/(m·K)
水的导热系数${\lambda _{\text{w}}}$	λf(θ)	W/(m·K)

下载: 导出CSV 

| 显示表格

COMSOL模拟结果和Wang等^[10]模拟结果的对比如图 4所示。图 4（a）是对称轴上桩顶部、中部、底部以及模型中部表层土的竖向位移变化，测点的具体位置在图中以（r，z）坐标的形式给出。由图可知能源桩顶部（0，0）处竖向位移向上，而底部（0，-0.38）处竖向位移向下，这主要是土骨架变形、土固结以及桩和土的热膨胀共同作用的结果。同时可以看出能源桩和表层土在初始时刻发生瞬时沉降，这是由于初始时刻轴向荷载瞬间施加于桩顶引起。随着温度逐渐升高，桩和土都相应地膨胀。如图 4（b）所示，剪应力σ_rz沿桩土界面随深度的变化也呈现类似现象。桩土界面上半部分剪应力为负，而下半部分剪应力为正，中部附近剪应力较小。如图 4（c），竖向应变包括弹性应变和热应变两部分，能源桩顶部的竖向应变较小，而中部和底部的竖向应变较大。由于顶部设置了恒温边界，因此顶部附近的温度更低，热应变较小，所以总竖向应变也较小。同时能源桩底部的竖向应变相对于中部也有降低，这是由于桩底应变受桩下土的作用。另外，桩土界面超静孔隙水压力在加载瞬间形成，尤其在桩底（z=-0.38 m）附近，但逐渐消散，在100 s时基本达到静水压力，可见相对于热传导和热对流，孔压消散较快。

图  4  本文结果与Wang等^[10]模拟结果的对比

Figure  4.  Comparison between simulated results and those from Wang et al^[10]

下载: 全尺寸图片 幻灯片

2.3   原位试验的有限元分析

本节通过THM耦合有限元模拟Laloui等^{[1, 33]}开展的能源桩现场原位试验，并进行对比分析。该原位试验数据完整，结果广泛用于数值模拟校核验证。共开展了8组能源桩热循环原位试验（试验1~8），每组试验为期28 d，其中加热持续12 d冷却持续16 d。能源桩同时承受结构荷载和热荷载，前者主要是在建建筑物的自重，桩顶上最大结构荷载约为1088 kN；后者则通过控制热交换管中的水温来施加，在每层施工结束时施加一次完整的热循环。进行试验1时施工尚未开始，桩顶无结构荷载；进行试验7时施工完成，桩顶结构荷载达到最大。因此本文选取这两组代表性试验进行模拟。试验1的最大热荷载Δθ=21℃，而其它试验的最大热荷载Δθ=15℃。

能源桩的几何模型和边界条件如图 5所示，该轴对称模型的半径为7 m，深度为52 m。能源桩半径为0.5 m，长度为26 m，桩周土是5层成层土。边界条件设置与2.1节类似，模型初始温度和顶部恒温边界均为15℃。将整个桩体作为热源，试验1中热荷载随时间变化如图 6所示^[1]。能源桩和土均假设为热弹性材料，桩土接触模型假定为绝对粗糙或库仑摩擦，界面摩擦角根据相关文献[16, 34, 35]确定，表 3给出了材料参数。网格采用非结构三角形单元，包含6399个单元和3455个节点，在桩土界面附近较密，远处则相对稀疏。

图  5  原位试验有限元模型的几何模型与边界条件

Figure  5.  Geometric model and boundary conditions for finite element analysis of field tests

下载: 全尺寸图片 幻灯片

图  6  试验1中施加在桩上的温度荷载时程曲线

Figure  6.  Time histories of temperature imposed on piles in Test 1

下载: 全尺寸图片 幻灯片

表  3  原位试验中桩和土的材料参数^{[1, 16]}

Table  3.  Material parameters for energy piles and soils in field tests ^{[1, 16]}

参数	土层A1	土层A2	土层B	土层C	土层D	能源桩
体积模量K/MPa	122	122	59	83	1860	17381
剪切模量G/MPa	113	113	1000	1400	1675	14313
密度ρ/(kg·m-3)	2000	1950	2000	2200	2550	2500
孔隙率n	0.1	0.1	0.35	0.3	—	—
水力传导率k/(m·s-1)	2×10-6	7×10-7	1×10-6	1×10-6	—	—
导热系数λ/(W/(m·K)-1)	1.8	1.8	1.8	1.8	1.1	2.1
比热容C/(J/(K·kg)-1)	1200	1231	1200	1091	784	800
体积热膨胀系数β/K-1	1×10-5	1×10-5	1×10-4	1×10-4	1×10-6	1×10-5
界面摩擦角δ/(°)	19.1	18.5	16.1	18.9	24.5	—
界面摩擦系数μ	0.346	0.335	0.289	0.342	0.456	—

下载: 导出CSV 

| 显示表格

图 7以试验1为例呈现了加热和冷却结束时的温度场模拟结果，可以看出受温度荷载影响的部分主要限于桩的深度，因此后续的部分结果仅给出模型上半部分桩长范围的模拟结果。

图  7  试验1加热和冷却结束时温度变化Δθ的分布

Figure  7.  Distribution of temperature variation at end of heating and cooling periods in Test 1

下载: 全尺寸图片 幻灯片

如图 8~12可以看出数值模拟与试验结果及其它已发表的数值结果吻合较好。试验1中桩顶无结构荷载，能源桩可自由上移，仅热荷载在桩中产生应力和应变。图 8将加热和冷却期间的竖向应变模拟值与测量值进行了对比，可以看出加热期间竖向应变不均匀且受界面摩擦力的影响。桩土接触为库仑摩擦时的模拟结果比绝对粗糙更合理，因为绝对粗糙大大高估了界面摩擦力，使得竖向应变分布特别不均匀。图 9中模拟的16 m深度处径向应变变化与试验结果基本一致，说明热循环期间桩土始终保持着径向接触。图 10给出了整个热循环过程桩顶竖向位移的时程曲线，可以看出位移呈明显弹性特征。土层A1、A2，土层B、C和土层D中的应力水平差异明显，如图 11所示。由于土层D模量很大且导热系数和热膨胀系数较小，其受荷载影响较小。试验7桩同时承受结构荷载和热荷载，模拟结果较好反映了能源桩的热-力行为。由于桩侧土的摩擦和桩底土的端承作用，能源桩不能自由移动，因此桩体会产生附加应力。Δθ=13.4℃时桩中竖向应力的模拟结果如图 12所示，与Saggu等^[14-15]的模拟结果较为一致。在土层B、C段桩中竖向应力较大，这是由于这两层土模量相对较大，对桩的约束更大，因此产生的热应力也更大。

图  8  试验1中热循环期间桩的竖向应变分布

Figure  8.  Vertical strain distribution during a heating-cooling cycle of the pile in Test 1

下载: 全尺寸图片 幻灯片

图  9  试验1中桩的径向应变时程曲线

Figure  9.  Time histories of radial strain for piles in Test 1

下载: 全尺寸图片 幻灯片

图  10  试验1中桩顶竖向位移时程曲线

Figure  10.  Variation of vertical displacement at pile head in Test 1

下载: 全尺寸图片 幻灯片

图  11  试验1中热循环期间土的竖向有效应力分布

Figure  11.  Distribution of vertical effective stress during a heating- cooling cycle of soils in Test 1

下载: 全尺寸图片 幻灯片

图  12  试验7中Δθ=13.4℃时桩的竖向应力分布

Figure  12.  Distribution of vertical stress along pile at Δθ=13.4℃ in Test 7

下载: 全尺寸图片 幻灯片

3.   结论与展望

本文基于热-水-力三场耦合理论，建立能源桩与土相互作用的THM全耦合有限元分析模型，研究能源桩的热-力响应，得到以下3点结论。

（1）利用多孔介质混合理论推导了THM全耦合控制方程，通过饱和非等温固结问题验证了该模型的有效性，表明采用全耦合有限元模型分析THM多场相互作用的可行性。

（2）考虑流体性质随温度的变化和不同的桩土接触模型，建立了基于THM全耦合的能源桩轴对称有限元分析模型，并与离心数值模型和现场试验结果进行了对比，表明该模型能够准确模拟能源桩复杂的热-力行为，且能较好地再现试验结果，可为能源桩相关多场耦合问题的数值分析提供参考。

（3）模拟结果表明，相对于热传导和热对流，超静孔压的消散较快，热对流对分析结果的影响较大，孔隙压力的影响也不可忽略，因此考虑THM耦合效应非常重要。

土的本构模型和长期的热循环都会显著影响能源桩和土的承载特性，并对地源热泵系统的工作性能产生影响。因此，后续研究中将考虑引入土的弹塑性本构关系，分析能源桩及周围土体在长期热循环下的各种耦合响应。