能源地下结构是利用建筑基础构件以及地下结构物与周围岩土体进行热交换的新型岩土结构^[1]。主要包括能源桩、能源地下连续墙、能源隧道等。时至今日，能源地下结构工程应用遍布世界，瑞士苏黎世机场新航站楼E建筑物采用440根现浇混凝土灌注桩支撑上部结构荷载，其中约300根桩为能源桩，能够满足整栋建筑物65%的供暖与70%的制冷需求^[2]。英国宝格丽酒店是第一个应用能源地下连续墙的工程案例^[3]，利用高36 m，壁厚800 mm并带有护板墙板的地下连续墙体来支撑一个开挖6层、深24 m的地下室空间，同时承受上部结构荷载。该工程计划将长100 m的地下连续墙体(总共3 600 m²)改造为能源地下连续墙，使之满足150 kW的供暖/制冷需求。奥地利延巴赫地区直径12 m的双轨铁路隧道通过安装热交换系统，将传统隧道改造为能源隧道，并为市政提供了15 kW的热负荷^[4-5]。国内方面，南京朗诗国际街区在1 200根桩基础中埋设了U型或W型传热管，桩长约为30 m, 为建筑面积10万m²、高21.6 m的建筑物节省了40%以上的耗能。上海自然博物馆新馆地下2层展厅下方为地铁十三号线区间，两者整体建构，并采用国内首创灌注桩与连续墙地埋管方案，承担了夏季冷负荷1 366 kW，冬季热负荷1 600 kW^[6]。上海世博会汉堡馆在43根钻孔灌注桩(桩长30 m，直径600 mm)中埋设3U型传热管，为面积3 113 m²、高18 m的建筑物提供能量。为了综合考虑换热需求，在建筑物外围布设了302口深60 m、间距5 m的换热井，以提高换热功率。此外，天津市梅江综合办公楼、同济大学旭日楼、天津市塘沽凯华商业广场等项目也广泛应用了能源桩技术^[6]。

上述能源地下结构工程中，能源桩最具代表性。作为地源热泵技术的一种新型地下换热器，将原本通过钻孔地埋的换热管路埋设于桩基础内，利用桩身混凝土较好的热传导性能和桩周岩土体的更大换热面积，能够提高地下换热器的换热性能，并节省大量钻孔费用和地下空间资源^[7]。英国与瑞士已相继提出了能源桩的技术标准^[8-9]。由清华大学主编、住房和城乡建设部颁布的“桩基地热能利用技术标准”^[10]填补了中国在能源桩技术领域的空白，使得中国桩基地热能利用工程的勘察、设计、施工、验收、运行和维护有了技术指导。

能源桩换热过程中，温度升降与冷热循环将导致桩体、尤其是桩周岩土体温度改变，桩土界面性质发生变化，进而影响其变形和承载能力。因而，桩周土体的热力耦合行为与能源桩变形和承载力紧密相关。岩土体热力特性的试验与本构模型研究也是能源桩热力特性研究的基础。因此，有必要建立一个能够合理描述岩土材料热力耦合特性的本构模型，并实现模型在地下能源结构工程设计计算中的应用。

为此，本文首先在室内单元试验的研究成果基础上，归纳了温度改变对饱和黏土固结蠕变体积变化的影响规律，同时总结了既有热力耦合本构模型中存在的问题。然后介绍了清华热力学岩土本构模型(TTS模型)理论框架与特点，简要给出模型参数的物理意义和标定方法。最后开展了饱和黏土升温排水试验的有限元分析和能源桩温度-应力耦合有限元分析，为能源地下结构的设计和安全服役提供理论依据。

1 饱和黏土热力耦合的体变规律与本构模拟

现有热固结试验大多在一维固结仪或三轴压缩仪中开展。对于一维压缩条件的温度固结排水试验，通常在传统固结仪基础上，增设水浴加热装置以实现对土样的温度控制，同时采用百分表量测土样的高度变化^[11]或用位移传感器量测土样的径向与轴向变形^[12]，以此获取土样的体积变化规律。对于等向应力状态的变温固结排水试验，通过在压力室外布设加热并在压力室内安置温度传感器^[13]，以实现对土样的温度控制。通过量测土样排水量来获取土样体积变化的方式简单易行，但忽略了土颗粒和水的热膨胀，从而导致试验误差。采用监测压力室内部空间的体积变化，进而反算土样体积变化则更为精确^[13]。

基于上述试验仪器与体积量测方法，国内外学者开展了大量饱和黏土升温-降温固结排水试验^[14-16]，从中可总结得出以下试验现象(图 1)：

1) 正常固结土在升温过程中均发生体积收缩，温度体应变关系曲线具有非线性，此后的降温过程仍发生体积收缩，温度体应变关系曲线近似线性。

2) 轻微超固结土升温体积先膨胀后收缩，降温体积收缩；重超固结土升温体积膨胀，超固结比(overconsolidation ratio, OCR)越大，体积从膨胀转变到收缩的温度也越高。

3) 正常固结土多次温度循环过程中有显著的累积体缩变形产生，而超固结土温度循环后累积体积膨胀变形，累积的体缩/体胀变形随循环次数的增加而减少，多次循环后累积变形基本稳定。

上述试验现象与升温过程中吸附结合水和颗粒间应力的改变有关。一方面，温度的升高导致结合水弱化，土颗粒变得更紧密；另一方面，土颗粒和水的热膨胀也会引起Van der Waals引力和双电子层斥力发生改变，导致颗粒接触力发生调整^[17]，最终导致颗粒间接触力随土体温度的升高逐渐降低，并引起土骨架发生结构坍塌，孔隙比减小。土体完成结构坍塌后，颗粒间产生了附加接触力，并使颗粒间接触在更高的温度下达到新的平衡。对于正常固结土，升温过程颗粒重排列引起的塑性变形占主导。而对于超固结土，因为土体较为密实，颗粒重排列受到限制，土颗粒和水的热膨胀弹性变形占主导^[18]。

尽管现有试验获得了较为丰富的研究成果，但仍存在一些问题值得深入探讨：1)体积变形量测方法多样，没有比较量测方法之间的差异性；2)黏性土温度固结排水试验时间往往较长，土体在恒定有效应力条件下会发生蠕变变形，缺乏相关蠕变试验与变温试验对比，导致对现有试验结果的分析隐含了蠕变的影响；3)有学者认为土的温度效应与塑性、液限有关^[11]，但目前土性对蠕变变形与温度变形的影响研究相对匮乏；4)缺乏细致的热弹耦合、热塑耦合、温度蠕变机制的研究。

现有描述热力耦合的本构模型大多以临界状态硬化弹塑性理论框架为基础，针对温度变化对土体弹性变形和塑性变形的影响，本构模型可分为仅考虑热弹耦合的弹塑性模型和热弹耦合-热塑耦合模型。前者代表模型有日本名古屋工业大学提出的热弹-黏塑性本构模型^[19]。该模型假定塑性部分与温度无关，但与应变率相关，通过引入热弹性线膨胀系数描述温度改变引起的土体弹性体积变化。但对于正常固结土的升温体缩试验现象，熊勇林等^[20]认为是热-水-力耦合边值问题的结果而非材料的本构行为。尽管模型能够在有限元模拟中实现正常固结土升温体缩特点，但对试验结果定量描述的能力不足。热弹耦合-热塑耦合模型的代表有变温剑桥模型、ACMEG-T模型、UH模型^[21-23]，这些模型主要既考虑了温度对弹性模量的影响，也能描述温度对土体屈服特性的影响，可以定性描述正常固结土与超固结土升温降温过程的塑性体积变化，但对温度-体积应变曲线的非线性描述不足，无法描述温度循环的累积变形，也无法描述与温变导致的体积变形同时发生的蠕变。在有限元分析方面，现有工作大多采用热弹耦合的本构模型进行^[24-26]，由于上述热塑本构模型的复杂性，有限元分析工作鲜有报道，其模拟能力没有得到很好的体现。

2 TTS模型

TTS模型是以非平衡态热力学为基础的新一代岩土本构模型，能够为岩土材料力学响应提供更本质的物理基础。较弹塑性本构模型而言，TTS模型摒弃了经典塑性力学中的若干概念(屈服函数、流动法则、硬化规律等)，以颗粒固体流体动力学(granular solid hydrodynamics, GSH)为理论框架，能够统一描述岩土材料温度效应、临界状态、率相关性、各向异性、循环滞回、多场耦合、小应变刚度等重要性质^[27-29]。限于篇幅，本文仅给出温度-应力耦合的TTS模型框架。

2.1 理论假定

TTS模型主要有3点理论假定：

1) 假设饱和多孔介质由固相和液相组成，液相又由结合水及自由水组成。随着温度的变化，结合水与自由水之间存在质量交换^[30]。

2) 假设各相在空间上均为连续分布，整体方程则可表示为各相方程的叠加。

3) 假设不同相的温度相同，液相不产生固化和汽化。

2.2 弹性势能与有效应力

TTS模型采用超弹性理论来描述岩土材料的非线性弹性行为。在超弹性理论中，材料的应力应变关系通过自由能函数(弹性势能函数或应变余能函数)得到，完全满足热力学基本定律。为了合理描述岩土材料压硬性、压剪耦合与非等向固结剪切刚度，TTS模型的弹性势能函数为

$ \omega_{\mathrm{e}}=B\left(\varepsilon_{\mathrm{v}}^{\mathrm{e}}+c\right)^{n}\left[\left(\varepsilon_{\mathrm{v}}^{\mathrm{e}}\right)^{2}+\xi\left(\varepsilon_{\mathrm{s}}^{\mathrm{e}}\right)^{2}\right]. $

(1)

其中: ω_e为弹性势能函数，ε_v^e和ε_s^e分别为弹性体积应变和弹性剪切应变，B、c、n和ξ为模型参数。

温度变化会引起土骨架发生弹性热膨胀，进而引起弹性势能的改变。为了描述温度对弹性势能函数的贡献，TTS模型在式(1)的基础上引入了一个热弹耦合项，提出了考虑温度的弹性势能函数：

$ \omega_{\mathrm{e}}=B\left(\varepsilon_{\mathrm{v}}^{\mathrm{e}}+c\right)^{n}\left[\left(\varepsilon_{\mathrm{v}}^{\mathrm{e}}\right)^{2}+\xi\left(\varepsilon_{\mathrm{s}}^{\mathrm{e}}\right)^{2}\right]+\int K \beta_{\mathrm{T}} \Delta T \mathrm{d} \varepsilon_{\mathrm{v}}^{\mathrm{e}}. $

(2)

其中: K为弹性体积模量，β_T为热膨胀系数，ΔT为温度增量, Kβ_TΔT为完全约束热弹性膨胀所需的附加应力，∫Kβ_TΔTdε_v^e为该附加应力在弹性体积变形上所作的功。通过式(2)可得到考虑温度的饱和黏土有效应力公式，即

$ \begin{array}{c} \sigma_{i j}^{\prime}=B\left(\varepsilon_{\mathrm{v}}^{\mathrm{e}}+c\right)^{n-1}\left[(n+2)\left(\varepsilon_{\mathrm{v}}^{\mathrm{e}}\right)^{2}+2 c \varepsilon_{\mathrm{v}}^{\mathrm{e}}+\right. \\ \left.n \xi\left(\varepsilon_{\mathrm{s}}^{\mathrm{e}}\right)^{2}\right] \delta_{i j}+K \beta_{\mathrm{T}} \Delta T \delta_{i j}+2 B \xi\left(\varepsilon_{\mathrm{v}}^{\mathrm{e}}+c\right)^{n} \varepsilon_{i j}^{\mathrm{e}}. \end{array} $

(3)

式(3)中的变量、参数和式(1)、(2)相同。

2.3 非弹性与能量耗散

1) 颗粒熵增方程。

岩土颗粒材料受外界因素(应力、温度)扰动时，除了土骨架的弹性变形，颗粒之间的摩擦、碰撞、滑移、滚动等非弹性相互作用迫使单一颗粒发生无规则的涨落运动(图 2)。涨落运动的强弱程度反映了土体颗粒间非弹性相互作用的大小，是颗粒固体引起其不可恢复变形与能量耗散的主要物理机制。为了合理描述颗粒固体的涨落运动，Jiang和Liu^[31]提出了一个新的独立状态变量——颗粒熵，取代了弹塑性理论中的塑性应变、等效塑性应变或塑性功，甚至黏滞应变等变量，以非平衡态热力学的视角，俯瞰岩土的非弹性变形等宏观耗散行为关系，即颗粒熵增方程：

$ \begin{array}{c} \rho_{\mathrm{s}} \dot{T}_{\mathrm{g}}^{\prime}=m_{2} m_{4} \dot{e}_{i j} \dot{e}_{i j}+m_{3} \dot{\varepsilon}_{\mathrm{v}}^{2}+ \\ \frac{m_{5} \pi_{k k} \alpha_{\mathrm{bf}} \phi_{\mathrm{b}}}{3 \phi_{\mathrm{s}}} \dot{T}-m_{4} T_{\mathrm{g}}^{\prime}. \end{array} $

(4)

其中: ρ_s为颗粒密度，π_kk为平均有效应力，T_g为颗粒熵温度，m₂、m₃、m₄、m₅为迁移系数，α_bf表示单位质量结合水升高单位温度向自由水的质量转移，ϕ_b和ϕ_s分别为结合水体积分数和固相体积分数。

2) 弹性弛豫方程。

在外界因素作用下，岩土材料除了发生颗粒涨落，颗粒之间的碰撞、滑移和滚动等相互作用还会迫使土体偏离弹性接触状态。土体偏离弹性接触的状态也称之为暂态弹性状态。而土体由暂态弹性状态向稳态弹性状态过渡时，往往伴随着颗粒间持续的相互作用，并引起土体弹性能量的减少，非弹性应变的增加，这个过程称为“弹性弛豫”(elastic relaxation)。弹性弛豫方程为

$ \dot{\varepsilon}_{i j}^{\mathrm{e}}=\dot{\varepsilon}_{i j}-\left(m_{1} T_{\mathrm{g}}^{\prime} \alpha_{\varepsilon}^{\mathrm{e}} \delta_{i j}+T_{\mathrm{g}}^{\prime} \alpha_{i j}^{\mathrm{e}}\right). $

(5)

式(5)中，m₁为迁移系数。

2.4 模型参数

TTS模型参数包括超弹性参数4个、能量耗散参数或迁移系数4个、温度相关参数4个，共12个，如表 1所示。

1) 参数n、B、c、ξ。

参数n为常数，取值范围为0.5~1.5，为了避免数值奇异，通常取1.5。参数B为干密度的函数，表达式为B=B₀exp(B₁ρ_d)。其中，B₀和B₁分别与正常固结线在e-logp平面上的截距和斜率有关(或者说与剑桥模型中参考孔隙比和压缩指数有关)。参考孔隙比越大，B₀越大；压缩指数越大，B₁越小。通过等向固结试验的压缩曲线和回弹曲线确定。参数c与土的黏聚力有关，参数ξ与土的临界状态应力(p_cs，q_cs)有关，可通过不排水三轴剪切试验确定。

2) 迁移系数相关参数m₁、m₂、m₃、m₄。

参数m₁控制土体非弹性应变的变化率，由不排水剪切试验确定；参数m₂反映体积应变率和剪切应变率对颗粒熵温度变化速率的影响，与临界状态应力比有关；参数m₃由等向固结试验确定；参数m₄控制颗粒熵温度变化程度，通过应力松弛试验确定。

3) 温度相关参数β_s、α_bf、k_g、m₅。

热膨胀系数β_s与矿物成分有关，可取为各矿物热膨胀系数的加权平均。结合水转化自由水速率参数α_bf，可根据结合水含量与温度关系确定。参数k_g和m₅通过恒定有效应力的热固结试验标定。

3 基于TTS模型的有限元分析 3.1 高岭土升温排水试验有限元分析

为了研究温度-应力耦合对饱和土体变化规律的影响，本文在自行编制的轴对称有限元程序基础上，对饱和高岭土的升温排水试验进行模拟。试验详细过程参见相关文献[32]。基于上述参数标定方法，得到了高岭土的模型参数，如表 1所示。

本文采用3种模拟方式：方式一将本构模型退化为热弹-黏塑性模型，即忽略温度增量对颗粒熵增方程和弹性弛豫方程的影响，以此描述饱和土的热弹耦合性质。方式二本构模型考虑为热弹-热塑模型，但假定结合水体积热膨胀系数为常数。方式三的本构模型为热弹-热塑率相关模型，并且能够全面考虑升温过程和结合水向自由水的转化机制。

方式一的模拟结果如图 3a所示，升温初期土的孔压开始累积，有效应力减小，弹性应变增加，土体膨胀。随着升温的进行，主固结过程加剧，宏观上土体由体胀转变为体缩。超固结比对饱和土升温过程的体胀与体缩均有影响。因此，热弹耦合的本构模型能够定性描述升温体变的OCR相关性。

考虑热塑机制后(方式二)，模型能够描述升温导致的塑性应变累积，即热蠕变。由于热蠕变在正常固结土和轻超固结土中表现得更为显著，因此模拟结果与试验结果也更加吻合(图 3b)。也说明考虑热弹-热塑机制的TTS模型同时具备描述热胀冷缩的孔压增长、主固结与热次固结的能力。但难以描述重超固结土升温过程的体积变化规律。

方式三考虑了温度对结合水热膨胀系数的影响，模拟了OCR=1~12高岭土的升温排水试验。结果表明：通过引入升温过程中结合水向自由水转化这一物理机制，TTS能够很好地描述重超固结土升温先膨胀后收缩的现象，也能模拟体积随温度显著的非线性变化(图 3c)。

3.2 能源结构温度-应力耦合平面应变有限元分析

本节应用TTS模型开展能源结构温度-应力耦合的有限元分析。图 4为能源基础有限元模型，能源基础尺寸为长10 m，宽1.2 m，总单元数为520。土体地基土干密度为1 550 kg/m³，TTS模型参数如表 1所示。具体计算过程为：首先进行重力加载(地基温度22℃)，再对基础中心线施加温度荷载，从初始温度22℃加热至82℃，随后降温至初始温度，最后对基础顶部施加竖向荷载。

图 5为排水和不排水2种条件下不同温度的荷载-位移曲线。给定参数条件下，温度升高竖向荷载有微小降低，通常温度对竖向荷载的影响主要与土体的热力特性和界面性质有关。一方面，升温导致土体剪切强度和刚度发生弱化，降低了土体对桩的约束能力。另一方面，升温过程中，由于桩周土体的升温体缩，界面正应力发生衰减，并导致界面切应力的减小。随着基础深度的增大，界面切应力随温度升高呈非线性增加，而界面正应力变化较小。

图 6为能源结构不排水条件下，初始温度(情况一)、升温后(情况二)、先升温再降温(情况三)3种情况下的荷载-位移曲线。对比发现，升温后(情况二)能源基础的位移增大，后续的降温(情况三)中位移继续增大。原因在于，建筑基础的热应变基本恢复，而土体刚度和强度发生弱化，且产生累积变形。由于结构和土体在变温条件下响应的差异性导致了界面接触的减弱，增加了基础的沉降。此外，排水条件和土体密实度对结构沉降的影响显著，排水条件和土体低密实度情况下，基础沉降量更大。

4 结论

本文总结了饱和黏土温度-应力耦合的体积变化规律，指出了现有热力耦合本构模型存在的问题，阐述了TTS模型的理论框架，开展了饱和高岭土升温排水试验和能源桩的有限元分析，主要结论如下：

1) 饱和黏土的热应变与OCR有关，正常固结土升温降温过程均发生体积收缩；超固结土升温膨胀降温收缩；累积变形随循环次数的增加而减少，多次循环后累积变形基本稳定。热塑性本构模型能够较好描述黏土的热力特征，因本构模型的复杂性，尚未开展相关的有限元分析工作。

2) TTS模型是热弹-热塑率相关模型，可以反映岩土材料热胀冷缩的孔压增长与主固结、温度蠕变和热塑率相关性。模型特定参数8~10个。

3) 基于TTS模型，在给定参数条件下，在升温过程中，能源基础的承载能力有微小降低。界面正应力、切应力均发生衰减，随着基础深度的增大，界面切应力随温度升高呈非线性增加，而界面正应力变化较小。升温降温后，桩土界面的接触减弱，基础沉降有所增大、承载能力有一定降低。

4) 现有模拟中，仍以热弹耦合的本构模型为主，有必要进一步开展以热塑耦合本构模型为基础的有限元模拟；值得注意的是，现有热塑耦合本构模型没有考虑温度蠕变的影响，开展热弹-热黏塑性本构模型研究也十分重要。