发展能源地下结构是实现建筑节能的有效手段. 相较于传统的钻孔埋管，将换热管路埋设在桩基和隧道等地下结构中，具有节省空间、降低成本和提高换热效率等优点．管路中水温的季节性变化，会引起地基土温度场的改变，导致能源结构产生附加温度应力^[1]和沉降^[2]，影响土-结构界面力学性质^[3]．式（1）为排水条件下侧摩阻极限承载力τ's的设计公式^[4]，

τ's=Kσ'z tan δ（1）

其中，K为水平土压力系数；σ'z为有效自重应力；摩擦角δ介于残余摩擦角和峰值摩擦角之间. 温度变化能够引起土压力Kσ'z和摩擦角δ的改变，对摩擦型桩的承载力有较大影响．因此，在承载力的设计中，亟需考虑温度对地基土力学性质的影响.

土体热力耦合行为研究是分析界面力学性质的基础，通过加入温度和孔压传感器的常规三轴仪、固结仪和真三轴装置，探索不同温度条件下超固结比和非轴对称应力状态等对温度应变的影响．研究表明，在对称三轴应力状态下，对于饱和黏土，一定温升范围内随着超固结比增加，体积应变由收缩转为膨胀，但非常规三轴应力状态下温度改变引起的应变分量的变化还鲜有研究．剪应力的存在使土体产生热剪切变形^[5-6]．在不排水温度循环条件下，存在发生热破坏的临界剪应力水平^[5，7]．温度对于临界状态摩擦角的影响不大^[8-9]．土体等温加卸载过程中的水平侧压力系数( K₀ )，一般参考MAYNE等^[10]根据大量实验总结的经验公式，正常固结土的K₀满足JAKY公式．目前没有变温条件下的K₀测定试验．

热弹塑性本构模型大多采用温度修正屈服面以模拟体积变化，为模拟循环变温的应变积累，需引入边界面^[11]．这些模型较少模拟应力状态对温度响应的影响．对温度引起剪应变缺少验证，正确模拟温度剪应变的关键在于定义合理塑性势函数，而且一般缺少对K₀的验证．弹塑性本构如修正剑桥模型（modified Cam-Clay model，MCC），其规定的正常固结下水平土压力系数K ^nc0 与M，压缩系数λ、κ和卸载泊松比ν'有关，通过临界状态确定的M往往导致K ^nc0 偏大^[12]．卸载过程为弹性，K ^nc0 与超固结比（over consolidated ratio，OCR）成正比，且和卸载泊松比有关，难以准确模拟K₀^[12]．

考虑温度对土-结构界面性质影响的模型或现场试验较少．目前有直剪试验^[13]、离心机试验^[14]、原位测试水平土压力和侧摩阻力等^[3]．现场试验发现饱和砂土地基中对能源桩循环温度加载后水平土压力降低^[3]．清华热力学岩土（Tsinghua thermo‐dynamical soil，TTS）本构模型可以合理描述多种热力耦合路径下的应力和应变响应^[15]．ZYMNIS等^[2]已经用TTS预测不同布置方式和换热效率下钻孔换热器群的沉降及差异沉降，并给出设计参数建议．本研究基于TTS本构模型进行单元试验研究，模拟平面应变条件下的温度体应变及应变分量，解释温度剪应变产生的机制，预测了温度循环对K₀的影响，研究结果对能源桩侧摩阻承载力设计有参考价值．

清华热力学岩土模型以颗粒固体流体动力学（granular solid hydrodynamics，GSH）为理论框架，给出了适合于岩土材料的弹性势能密度、迁移系数、耗散流函数的形式．本研究采用的TTS本构模型，适用于饱和黏土和砂土，在本构上通过两个机制实现应力场和温度场耦合：①热弹耦合，考虑温度对弹性应变的影响；②热塑耦合，考虑升温引起颗粒温度增加．

TTS通过弹性势能函数ω_e规定弹性应力和弹性应变的关系．ω_e表达式为

其中，B、θ、ξ和c为模型参数；ε^e_v为弹性体积应变；ε^e_s为弹性剪应变；∆t=t-t₀，t为当前温度，t₀为参考温度．B为

B=B₀e^B₁ρ_d（3）

其中，B₀为系数；B₁为模型参数；ρ_d为干密度，满足ρ_d=ρ_s φ_s，ρ_s和φ_s分别为固相本征密度和固相体积分数．割线弹性体积模量K为，

β_T为土骨架的体积热膨胀系数，由固相热膨胀系数β_s和结合水的热膨胀系数β_l决定^{[1 5]}，为

其中，φ_b为结合水体积分数.

在忽略黏滞应力的情况下，ω_e对弹性应变的微分π_ij即为有效应力，

其中，Kβ_T∆t为弹性应变被完全约束时的温度应力．

ρ_d作为模型的状态变量，可由质量守恒方程式得到，

其中， 为体积应变率；ρ̇d为干密度应变率.

ρ_s0和ρ_l0分别为固相、液相本征密度初值. 固相和结合水的本征密度分别为

ρ_s=ρ_s0[1-β_s(t-t₀) ]（8）

ρ_l=ρ_l0[1-β_l(t-t₀) ]（9）

由于试验数据较少，β_l暂时取水的体积热膨胀系数．水的密度和温度的关系如式（10），由此得到水的割线热膨胀系数表达式（11）^[5].

ρ=1 000. 07+0. 020 922 90t-0. 006 021 37t²+0. 000 016 30t³（10）

式（12）为颗粒熵增方程．颗粒温度T_g衡量颗粒无规则运动的程度．~T_g是为了简化方程形式采用的等价变量^[5]．方程包含升温时结合水向自由水转化引起颗粒熵增的机制．在外部激励下（如加载或升温，即体积应变率ε̇v、剪切应变率ε̇s或温度变化率ṫ都不等于0 ），颗粒温度增加；由于颗粒熵不断向真实熵转化，在无激励时颗粒温度降低，体系趋向于回归暂态弹性状态^[16]．

其中， m ₂、m ₃、m ₄、m ₆、α _bf和k _g为模型参数； φ _b由式（14）^[19]计算；φ_b0为初始结合水体积分数.

总应变由弹性应变ε^e和非弹性应变ε^D两部分组成.

非弹性应变率取决于~T_g和当前应力状态

（ε^e_ij）. 其中，e _ij为偏应变；m₁和a为模型参数．

为了考虑固结应力历史的影响，引入修正应变ε ^L_ij^[5]．修正后应变的计算方法为

其中，h和w为模型参数.

TTS模型参数包括等温和非等温参数两部分，共15个（热弹部分5个材料参数）．所有参数能够通过室内单元试验进行标定．标定两种土等温参数所用的试验为一维压缩、K ^oc0 经验规律和温度蠕变试验等温应力历史决定的蠕变方向．标定soft Bang‐kok clay非等温参数采用一维固结温度循环试验，标定Bonny silt非等温参数采用平面应变排水温度循环试验．

2. 1 等温参数

等温参数可通过一维压缩曲线、K₀随OCR的定性关系和热蠕变试验进行标定．ZYMNIS等^[17]给出的soft Bangkok clay参数，不符合K₀随OCR增加的经验规律，原因是符合一维压缩曲线的TTS参数不止一组，需要考虑更多应力路径得到最优解^[18].

参数θ取1. 5，a=0. 5^[5]，m₄取一个较大值．文献[9]中给出的初始干密度为789 kg/m³，先期固结应力约为70 kPa．状态变量初值ε ^e_ij=0，~T _g=0， ρ_d0取一个较小值，如600 kg/m³，相当于从泥浆状态开始固结．ρ_d0的取值仅影响正常固结线的起始段．图1为soft Bangkok clay^[9]和Bonny silt^[19]一维压缩模拟结果．其中，文献[6]给出Bonny silt压缩指数C_c=0. 04可能有误，故采用文献[19]的曲线标定．MAYNE等^[10]根据大量室内试验提出，正常固结下水平土压力系数与卸载、重加载时的水平土压力系数分别为

K ^nc0 =1-sinϕ'（19）

图1 （a）Bangkok黏土和（b）Bonny粉土压缩指数随压力的变化，数据点来源[9，19] Fig. 1 Variation of compression index with pressure，test points are from[9，19]，（a）Soft Bangkok clay and（b）Bonny silt.

其中，ϕ'为有效摩擦；OCR_max为卸载中最大超固结比.

图2是TTS和式（19）至式（21）预测的卸载、重加载K₀．式（20）和式（21）中OCR最大值由被动土压力系数K_p=(1+sinϕ')/(1-sinϕ')限制．TTS给出各点K₀的表达式，如式(22)，K₀由ξ、c、θ、ε_v^e和ε^e_s/ε^e_v决定．卸载时K ^oc0 随OCR的增加由参数m₁、m₂、m₃和h控制．参数B₀和B₁对K₀无影响，分别根据正常固结线（normal consolidation line，NCL）的位置和斜率^[17]标定.

图2 （a）Bangkok软黏土和（b）Bonny粉土由TTS模型和式（20）、（21）预测的K0与OCR关系Fig. 2 K0-OCR predicted by TTS and empirical formula (20) and (21). (a) Soft Bangkok clay and (b) Bonny silt.

TTS可以正确模拟正常固结下的K ^nc0 ，在NCL上，常规应力范围内K ^nc0 基本不变，如图3．在定量上，TTS预测的 K ^oc0 小于经验公式，模拟结果可以通过减小迁移系数m₁、m₂、m₃或h改善，但在超固结比较大时会导致温度应变方向不正确．修改弹性弛豫方程形式可能改善模拟情况，如将m₁还原为系数λ_v和λ_s^[15]．

2. 2 非等温参数

非等温参数用于定量描述温度应变及其非线性．尽管a取0. 5，本研究模型仍是温度率相关的，因此m₄的取值对模型有影响．由于缺少等温蠕变试验数据，m₄和a没有标定，通过调整m₆尽量拟合．固相热膨胀系数暂与ZHANG等^[15]标定的Kaolin土一致．热塑耦合参数及φ_b0通过试验曲线拟合得到．

图4（a）是soft Bangkok clay一维固结温度体积应变模拟．首先将竖向应力加载到200 kPa，然后卸载到超固结比分别为1、2、4和8，经过1次温度循环：22℃→90℃→22℃，试验中升温速率约为 0. 002 8℃/min（1 次升温时长 24 400 min）． WANG等^[20-21]通过有限元验证此实验样品内部及边缘孔压积累最大不超过0. 1 kPa，可以认为完全排水．本研究采用的方程形式与文献[17，21]不同，且满足K₀规律，因此参数取值不同．如图4（a）所示，模拟略高估OCR=2的应变，其余情况定量较好．图4（b）是50次温度循环后的应变积累，不同应力状态下应变趋于稳定，正常固结下的最终应变约27%．

图3 TTS模拟NCL上各点K nc0 随加载压力p的变化Fig. 3 Variation of K nc0 with p load-on（blue line），load-off （red line），and re-loaded（yellow line）during normal consolidation simulated by TTS.

Bangkok黏土不同OCR下K₀随温度循环的变化如图4（c）至（f）．OCR=1，第1个循环升温时K₀先略增加再略减小，降温时K₀减小；OCR>1，第1个循环升温时K₀减小，降温时K₀减小．

图4 Soft Bangkok clay温度应变及温度循环下K0预测（a）不同超固结比下温度体积应变；（b）温度循环下应变积累；（c）OCR=1、（d）OCR=2、（e）OCR=4、（f）OCR=8条件下K0随温度循环变化的预测Fig. 4 Simulation of thermal volume strains and prediction of K0 under thermal cycles of soft Bangkok clay. （a）thermal volume strains under different OCR；（b）thermal strain accumulation under thermal cycling；prediction of K0 under thermal cycles with （c）OCR=1，（d）OCR=2，（e）OCR=4，and（f）OCR=8．

OCR决定K₀随着循环次数的变化，与OCR决定温度体积应变方向的机制一致．当1次升降温体积收缩时，K₀随着循环减少，当1次升降温体积膨胀时，K₀随着循环增加．如图5，弹性剪切应变ε^e_s和修正剪切应变ε^L_s、修正体积应变ε^L_v的变化很小， ε^e _v最终等于ε^L_v . ε^e _v>ε^L_v时，ε̇Dv =m₁~T_g ^a( ε^e _v-ε^L_v )>0 ，升温体积收缩，ε^e_v随循环减少，由式(22)，K₀随着循环减少；ε^e_v<ε^L_v时，升温体积膨胀，K₀随着循环增加．循环后的K₀最终值由等温条件下的弹性应变ε^e和修正应变ε^L决定，升降温速率增加，或升降温幅值减小，只是增加达到稳定需要的循环次数．

COCCIA等^[6]改进的平面应变试验装置中，立方体土样的上下两个平面为刚性平面，以限制土样的位移，并控制温度和孔压；其余4个面用气囊施加应力，试样处于平面应变状态．

图5 （a）OCR=1、（b）OCR=2、（c）OCR=4、（d）OCR=8条件下εev、εes、εLv和εLs在温度循环中的变化Fig. 5 Variation of εev，εes，εLv，εLs in thermal cycles with（a） OCR=1，（b）OCR=2，（c）OCR=4，and（d）OCR=8.

表1为各试样温度循环前的应力状态．其中， σ'x、σ'y和σ'z分别为x、y和z方向的应力. 试验材料为饱和Bonny粉土．首先在4个试样的x和y方向施加固结应力580 kPa，然后两个方向均卸载到440 kPa，得到K1. 0组的应力状态．然后对其余3个试样的y方向继续卸载，分别得到σ'y/σ'x=0. 9、0. 7和0. 5．z方向应力σ'z和OCR为TTS模拟结果，可以看出各试样OCR接近，均在1～2．OCR采用平均应力定义，2. 2节模拟一维固结状态温度蠕变时OCR采用竖向应力定义．

表1 Bonny silt平面应变试验温度循环前的应力状态Table 1 Stress state before thermal loading of Bonny silt samples in plane strain test

每个试样经过1次温度循环（25℃→65℃→25℃），试样中心的传感器显示内部升温速率为1℃/h．

平面应变下Bonny silt的温度应变如图6．其中，虚线是实测结果，实线是TTS的模拟结果．试验结果表明，升温时，应力比σ'y/σ'x越小，x方向（即应力保持方向）的温度应变越大，y方向的温度应变越小，并在σ'y/σ'x=0. 7和0. 5时出现收缩，而样品的体积应变几乎相同．

降温时，体积应变呈现出较强的非线性变化，试验土的降温应变幅值与应力状态无显著关系， K1. 0试样出现降温体变先收缩后膨胀．

TTS对于K0. 7和K0. 9的体应变模拟不准，大主应力方向应变偏大，试验数据呈现的规律较复杂，具体机制有待研究．K1. 0和K0. 5体应变及两个方向的应变模拟较好．

两个主应力方向的应变分量不同，主要是温度引起蠕变的贡献．温度加快蠕变速率，升温前各个方向的近期加载历史不同，导致各个方向蠕变的大小、收缩或膨胀性质不同，产生温度剪应变．TTS模型对此的解释为：随着卸载，ε^e_v减小，y方向的弹性应变ε^e_y减少，修正体积应变ε^L_v及其分量基本不变．K0. 7和K0. 5情况下，升温时非弹性应变率 ，因此出现膨胀．由于各试样OCR接近，弹性体积应变接近，升温时 ，因此均出现升温体积收缩．

图6 平面应变下Bonny silt温度应变（a）体积应变；（b）x方向应变；（c）y方向应变Fig. 6 Thermal strain of Bonny silt under plain strain situation for (a) volume strain, (b) strain in x direction, and (c) strain in y direction.

图7是根据试验绘制的温度剪应变发展曲线．其中， ε _s= 2 3 ( ε_x-ε_y )，箭头所指位置为降温开始点．剪应力水平越高，升温引起的剪应变越大．降温的规律不明显，K0. 5和K1. 0几乎无剪应变增长，K0. 7剪应变显著增加．TTS模拟的升温剪胀关系与试验符合较好，降温中几乎无剪应变发展，但随卸载后应力比略有增加．热弹塑性本构，若要合理模拟温度剪应变，需要定义合理的流动准则．温度剪胀曲线与三轴剪胀曲线的关系尚未得到确证．

表2 Soft Bangkok clay和Bonny silt的TTS模型参数取值Table 2 TTS parameters for soft Bangkok clay and Bonny silt

图7 平面应变下Bonny silt温度剪应变与体应变的关系Fig. 7 Thermal shear strain versus thermal volumetric strain of Bonny silt under plain strain situation.

采用与第3节相同的试验条件，将材料改为Bangkok黏土，进行平面应变下Bangkok黏土的温度应变预测，如图8．等温条件下，TTS给出各个应力比下的OCR与Bonny粉土接近，如表3．仍具有升温时应力比σ'y/σ'x越小，x方向的温度应变越大，y方向的温度应变越小，并在σ'y/σ'x=0. 7和0. 5时出现收缩，且体积应变几乎相同的性质．降温阶段与第2. 2节一维固结温度试验类似，几乎无应变发展．与Bonny silt相比，Bangkok软黏土在相同温升速率下收缩、膨胀蠕变均更大，由于结合水含量高，升温对蠕变速率的提高更明显．

表3 假想Bangkok软黏土平面应变试验温度循环前的应力状态Table 3 Supposed stress state before thermal loading of soft Bangkok clay samples in plane strain test

图8 平面应变下Bangkok黏土温度应变预测（a）体积应变；（b）x方向应变；（c）y方向应变Fig. 8 Prediction of thermal creep corresponding to soft Bangkok clay in plain strain situation for (a) volume strain, (b) strain in x direc‐tion, and (c) strain in y direction.

基于TTS模型模拟了一维压缩和平面应变状态下饱和黏土和粉土的热蠕变，分析了温度循环过程中K₀的变化规律，可知TTS可以定性模拟K ^oc0 与OCR的关系．真三轴试验各主应力方向上的热蠕变大小不同，产生温度剪应变．温度引起总体积应变收缩时，某些应变分量可能出现膨胀．基于TTS的本构模型可以合理描述这一性质．