0.   引言

能量桩是一项结合桩基埋管换热工艺与地源热泵技术的节能减排技术^[1],已逐渐受到国内外工程与学术界的广泛关注。在高层建筑工程及具有桥面融雪系统的桥梁工程中,一方面,能量桩需承担上部结构传递的倾斜荷载（包括上部结构自重引起的竖向荷载与风荷载、车辆制动荷载等水平荷载）；另一方面,还要承受桩–土换热引起的温度荷载。因此,热（温度荷载）-力（结构荷载）耦合作用下的能量桩承载机理与受力变形特性较为复杂,不同于常规工程桩,是目前影响能量桩这一新型技术推广应用的关键问题之一。

近年来,能量桩相关研究已有了较大的进展,开始从热响应测试与传热理论研究^[2-4]逐渐转向结构响应规律及机理分析^[5]。然而,针对倾斜荷载下能量桩受力变形特性的研究似未见报道,现有研究多聚焦于竖向荷载下能量桩的沉降与应力分析。现场试验方面,Laloui等^[6]探究了热-力耦合作用下能量桩的附加应力分布规律；Bourne-Webb等^[7]提出了降温工况下能量桩将产生附加温度拉应力及附加沉降,再度升温后,仍存在部分不可恢复的沉降；桂树强等^[8]通过单次冷热循环,揭示了附加温度应力与冷热循环下不可恢复沉降的形成原因；路宏伟等^[9]分析了摩擦型能源桩的荷载传递特性,升温或降温工况下,分别在桩身上部或下部出现负摩阻力,且桩顶竖向荷载可弱化升温导致的负摩阻力效应；方金城等^[10]着重分析能量桩热力学特性及其对桩周土体的影响规律；任连伟等^[11]初步探讨了冷热循环下能量桩的附加温度应力与桩侧摩阻力演化规律。模型试验方面,黄旭等^[12],孔纲强等^[13]提出了附加温度应力受桩周约束条件影响,降温工况下桩端甚至可能出现较大拉应力,且多次冷热循环下能量桩的累计沉降将影响上部结构的安全性；随后,陆浩杰等^[14]以长期冷热循环下摩擦型桩为研究对象,初步探讨了桩顶累积沉降的形成机理；此外,刘干斌等^[15]探究了在常温、升温与降温3种工况下桩周土体的热固结过程及桩侧负摩阻力的发挥机制。理论研究方面,Cesar等^[16]在桩单元压缩量计算中考虑温度影响,然其荷载传递函数假定为折线型,故难以考虑冷热循环过程中塑性变形的累积；Suryatriyastuti等^[17]提出了一个包含硬化-软化机制的侧摩阻力-深度循环函数,以反映冷热循环过程中能量桩的承载力变化；此外,费康等^[18]引入双曲线型荷载传递函数与Modjoin法则模拟冷热循环过程,并通过刚度折减近似考虑塑性变形的积累。相比而言,数值模拟方面的研究工作相对较少^[19-21],且多为基于竖向荷载下单桩的轴对称性质的二维数值分析。综上,针对竖向荷载下的能量桩沉降与应力分析,已取得了较多的成果；但实际工程中,桩基往往不止承受单一竖向荷载,水平荷载（风荷载、车辆制动荷载等）的影响也是不容忽视的。因而,对倾斜荷载下能量桩的受力变形特性展开三维数值分析是有必要的。

鉴于此,本文基于已有试验及数值分析,借助ABAQUS有限元软件,建立了倾斜荷载工况下考虑热-力耦合作用的能量桩三维数值模型：①通过与模型试验及二维数值结果的对比分析,验证本文三维建模方法及参数选取的合理性；②在此基础上,探究桩顶竖向荷载和温差对沉降及桩侧摩阻力变化规律的影响,并分析桩顶沉降的变化原因及负摩阻力的发挥机制；③通过桩身水平位移与桩身最大弯矩演化规律,分析升温及降温工况下,能量桩内力变形发生变化的原因；④定量分析桩顶沉降、桩顶水平位移及桩身最大弯矩值随温差变化的关系,并定义各自的增大系数以便求解其他工况下的桩身响应；⑤基于上述分析,提出相应的能量桩工程设计建议。

1.   三维有限元模型建立

本文拟通过优选几何尺寸、桩-土材料参数、边界条件、分析步与网格类型等,借助ABAQUS有限元软件进行数值试验,揭示热-力耦合作用下能量桩的受力变形特性,并通过MATLAB平台对数值结果进行后处理。

1.1   几何模型建立及参数选取

以某能量桩^[19]为例,桩径D=1.5 m,桩长L=8.2 m,由于桩轴线2倍桩径范围外,土体温度变化基本可以忽略^[20-21],并考虑边界效应的影响,将数值模型中桩周土体尺寸取为半径为10 m,深度为20 m的圆柱体,即桩端与土体底面的垂直距离为11.8 m。桩-土界面设置接触,将桩-土界面摩擦角考虑为δ=24.5°。桩体混凝土采用弹性本构模型,混凝土密度ρ_p为2200 kg/m³,孔隙率n_p为0.15；桩周砂土选择Mohr-Coulomb理想弹塑性模型,内摩擦角、黏聚力及剪胀角分别为φ=35°,c=1 kPa和ψ=9.40°,土颗粒密度${\rho }_{\text{s}}$为2650 kg/m³,孔隙率n_s为0.42；此外,弹性模量E、泊松比$\nu$、线膨胀系数$\alpha$、导热系数$\lambda$、比热容c_s等桩-土材料热力学参数如表1所示。

表  1  桩-土材料参数

Table  1.  Material parameters of pile and soil

材料	E/MPa	$\nu$	α/(με·°C-1)	λ/(W·(m·°C-1))	cs/(J·(kg·°C-1))
砂土	78	0.30	10	0.25	930
混凝土	33	0.25	16	1.7	930

下载: 导出CSV 

| 显示表格

1.2   边界、荷载条件及网格划分

能量桩三维数值模型的边界条件：①力学边界条件。土体底面视为固定边界条件,土体侧面约束径向位移,顶面自由；②温度边界条件。除桩顶外,模型外边界的温度均始终与初始温度T₀保持一致。

能量桩的荷载施加考虑采用在Edit keywords中材料属性后加上“initial condition, type=stress, geostatic”的方法对桩-土模型进行地应力平衡,同时考虑桩-土自重引起的初始地应力,将土压力系数取为K₀=1-sinφ=0.43；在施加重力荷载后,再以均布应力的形式在桩顶施加竖向荷载P；已知桩体与桩周土体的初始温度为T₀=25℃,分别将桩体升温至32℃,37℃,43℃,令温差ΔT=T₁-T₀,则温差为7℃,12℃与18℃,此时,桩顶应力分别为396,436,476 kPa^[19]。该算例的温度变化路径如图1所示,可通过自定义幅值曲线amplitude实现。

图  1  桩体温度变化曲线

Figure  1.  Variation of temperature of energy piles

下载: 全尺寸图片 幻灯片

为进行完全热-力耦合分析,桩体和土体均采用C3D8T（coupled temperature-displacement）特定网格类型,以保证网格的质量与分析的精确度。此外,为兼顾计算精度与计算效率,采用结构化网格技术对桩体进行网格划分,沿径向与桩长方向均匀划分,共960个桩体网格单元；对于较为复杂的土体模型,通过定义partition并结合扫掠网格技术对其进行网格划分,靠近桩体部分网格密集,远离桩体部分网格疏松,共5240个土体网格单元,如图2所示。

图  2  三维数值模型

Figure  2.  Three-dimensional numerical modelling

下载: 全尺寸图片 幻灯片

2.   数值模型验证

将桩顶附加沉降曲线与Rotta等^[19]的模型试验及二维数值结果进行对比,并规定竖直向下为正方向,如图3所示。

图  3  桩顶附加沉降曲线与试验及二维数值分析结果对比

Figure  3.  Comparison among additional pile head settlement curves, experimental and two-dimension numerical results

下载: 全尺寸图片 幻灯片

由图3可知,本文三维有限元模型计算结果与Rotta等^[19]的试验及二维数值分析结果吻合较好,尤其是初始阶段与最终阶段；而中间阶段模拟结果偏大的原因主要是试验中桩身温度荷载不均匀^[19]。此外,本文值与Rotta等^[19]二维数值结果的误差可能是桩-土界面接触、模型尺寸、网格划分等导致的,其中,Rotta等^[19]在桩-土之间单独设置一薄层模拟桩-土界面,而本文则采用设置接触对的方法进行模型。由此可见,本文三维有限元分析可以较好地模拟热-力耦合作用下能量桩的受力变形特性。

3.   受力变形特性

为进一步分析倾斜荷载作用下能量桩的受力变形特性,仍基于上述算例,探究桩顶沉降、桩侧摩阻力、水平位移及桩身弯矩的变化规律。在此分析中,考虑竖向荷载在0～4000 kN变化,水平荷载在0～800 kN变化^[22],温差ΔT工况主要考虑为0,±7℃,±12℃,±18℃。

3.1   桩顶沉降变化规律

考虑到桩顶沉降受水平荷载影响较小,故固定H=200 kN条件不变,通过改变桩顶竖向荷载P和温差ΔT,探讨桩顶沉降的变化规律,如图4所示。图中,规定竖向位移向下为正,向上为负。

图  4  桩顶沉降–竖向荷载曲线

Figure  4.  Curves of pile-head settlement-vertical load

下载: 全尺寸图片 幻灯片

由图4可知,当ΔT<0℃时,桩顶沉降w₀明显大于ΔT≥0℃工况,其主要原因为降温工况下桩顶收缩方向与竖向荷载方向一致,其沉降量为桩顶竖向荷载与温度荷载共同作用造成的^[9]；当ΔT=-18℃时,无荷载与P=1000 kN工况下的桩顶沉降值分别为2.077,4.975 mm；而当T=0℃时,无荷载与P=1000 kN工况下的桩顶沉降值分别为0,2.230 mm；显然,桩体降温收缩引起的桩顶附加沉降是不容忽视的。反之,当ΔT>0℃时,桩体将产生热膨胀,在P=0 kN与P=500 kN时,甚至会出现向上的位移（-2.230,-1.103 mm）,此时桩顶膨胀方向与竖向荷载方向相反,桩侧摩阻力变化较为复杂,但随着竖向荷载的增大,在桩顶荷载与桩端持力层的共同约束下,升温对桩顶沉降的影响逐渐减弱,w₀–P曲线逐渐趋近ΔT=0℃工况。在实际工程中,桩顶与桩端往往存在约束,沉降随温差的变化受到限制,故能量桩在正常使用状态下,将在桩顶与桩端产生不可忽视的附加温度应力（压应力或拉应力）。因此,在能量桩设计时,应对附加温度应力（尤其是拉应力）进行初步估算,使其低于桩体混凝土抗压及抗拉强度,并适当提高桩顶与桩端钢筋配筋率。

为定量分析温差ΔT对桩顶沉降w₀的影响,定义α(w₀)=Δw₀/w_0(ΔT=0℃)为桩顶沉降增大系数,其中,Δw₀为竖向荷载相同时,某温差工况下桩顶沉降与ΔT=0℃工况下桩顶沉降w_0(ΔT=0℃)的差值,探究桩顶沉降增大系数随温差的变化规律,如图5所示。由图5可知,桩顶沉降增大系数与温差的关系近似呈线性增大,其斜率随竖向荷载增大而减小,规律较为清晰,故对图中4条曲线进行拟合,可得其统一表达如下所示：

图  5  桩顶沉降增大系数–温差曲线

Figure  5.  Curves of increment coefficient of pile-head settlement-temperature difference

下载: 全尺寸图片 幻灯片

其余温差与荷载工况下桩顶沉降可通过式（1）计算出桩顶沉降增大系数α(w₀),并结合ΔT=0℃工况下桩顶沉降进行求解。

3.2   桩侧摩阻力变化规律

能量桩升温或降温时,桩体与桩周土体将产生相对位移,引起桩侧摩阻力的变化,不同深度土体的侧摩阻力可通过轴力及桩体自重求解。图6为不同温差工况下(P=0)桩侧摩阻力q_s随深度分布的情况。图中,规定桩侧摩阻力向上为正,向下为负。

图  6  温度荷载下桩侧摩阻力分布

Figure  6.  Distribution of pile-side frictions under temperature loads

下载: 全尺寸图片 幻灯片

由图6可知,本文数值所得桩侧摩阻力沿深度方向的分布规律与Laloui等^[6]、Bourne-Webb等^[7]、孔纲强等^{[4, 13]}的研究成果一致；中性点位置约在0.8L处,接近Laloui等^[6]现场实测的中性点位置,而Bourne- Webb等^[7]与孔纲强等^{[4, 13]}现场实测的中性点位置约在0.6L处；中性点位置存在差异的主要原因在于桩端约束条件的不同,本文算例与Laloui等^[6]的桩端持力层为砂土或砂岩,弹性模量或单轴抗压强度较大,对能量桩提供的桩端约束远大于Bourne-Webb等^[7]与孔纲强等^{[4, 13]}的黏土持力层,故其中性点位置更靠近桩端。此外,在此分析中,由于桩顶无竖向荷载,故桩体上部的桩侧摩阻力比下部更为明显；当ΔT>0℃时,桩身上、中部桩侧摩阻力为负,负摩阻力最大值为-8 kPa,下部桩侧摩阻力为正；当ΔT<0℃时,桩身上、中部桩侧摩阻力为正,摩阻力最大值为3 kPa,而下部桩侧摩阻力为负,反映了桩体升温膨胀或降温收缩时,桩周土体对能量桩的约束^[13]。综上,能量桩桩侧摩阻力分布规律及大小受温差影响较大,在能量桩设计时,应充分考虑负摩阻力的影响,如适当提高通长钢筋配筋率、选址时尽量避开高压缩性土层、避免大面积堆载等。

3.3   水平位移变化规律

取P=2000 kN,通过改变桩顶水平荷载H和温差ΔT,探讨桩身水平位移y的变化规律,如图7所示（图中位移为桩边缘处水平位移）。由于图中深度–水平位移曲线与刚性桩类似,通过计算可得,$\alpha h$=2.73≥2.5（其中$\alpha$=(mb₀/EI)^0.2,m为弹性地基比例系数,b₀为基桩计算宽度,EI为抗弯刚度）,故可将其视为刚性桩。

图  7  深度–水平位移曲线

Figure  7.  Curves of depth-lateral displacement

下载: 全尺寸图片 幻灯片

由图7可知,当ΔT>0℃时,桩顶水平位移y₀大于ΔT=0℃的工况,其主要原因为：升温时,桩体沿径向将产生热膨胀,此外,桩端因温度升高而产生桩土共同膨胀,相当于桩周及桩端约束加强^[8],限制了桩端位移。反之,当ΔT<0℃时,桩体沿径向将产生收缩,然而,桩土共同收缩将导致桩侧及桩端约束减弱^[8],旋转中心上移引起桩顶水平位移随之增大,上述两种作用相互抵消,尤其在水平荷载较小时（H=200 kN）,但随着水平荷载的不断增大,桩体自身径向收缩的影响逐渐可以忽略。此外,还可以观察到桩端位移与旋转中心深度均随ΔT的增大而增大,表明了温度变化将引起桩周及桩端约束的变化。综上,升温和降温工况下桩顶水平位移均有所增大,但各自的原因不尽相同,工程实际中,若要弱化温差的影响,可考虑整体加强桩周、桩顶及桩端的侧向约束,如嵌入岩层或扩大下部桩身半径等。

为定量分析温差ΔT对桩顶水平位移y₀的影响,定义$\alpha$(y₀)=Δy₀/y_0(ΔT=0℃)为桩顶水平位移增大系数,其中,Δy₀为水平荷载相同时,某温差工况下桩顶水平位移与ΔT=0℃工况下桩顶水平位移y_0(ΔT=0℃)的差值,探究桩顶水平位移增大系数随温差的变化规律,如图8所示。由图8可知,ΔT>0℃和ΔT<0℃工况下桩顶水平位移均大于ΔT=0℃工况,然其原因却不尽相同。此外,桩顶水平位移增大系数与温差的关系近似呈线性变化,其斜率随水平荷载增大而减小,规律较为清晰,规律较为清晰,故对图中4条曲线进行拟合,可得其统一表达式为

图  8  桩顶水平位移增大系数–温差曲线

Figure  8.  Curves of increment coefficient of pile-head deflection-temperature difference

下载: 全尺寸图片 幻灯片

其余温差与荷载工况下桩顶沉降可通过式（2）计算出桩顶水平位移增大系数$\alpha$(y₀),并结合ΔT=0℃工况下桩顶水平位移进行求解。

3.4   桩身弯矩变化规律

同样地,取P=2000 kN,通过改变桩顶水平荷载H和温差ΔT,探讨桩身弯矩M的变化规律,如图9所示。

图  9  深度–弯矩曲线

Figure  9.  Curves of depth-bending moment

下载: 全尺寸图片 幻灯片

由图9可知,各温差工况下桩身弯矩随深度的变化趋势大致相同。当ΔT>0℃时,桩身弯矩始终表现出增大的趋势,且其增量随ΔT的增大而增大,其主要原因是桩土升温时桩端及桩周约束加强。当ΔT<0℃时,同一深度桩身弯矩表现出：在水平荷载较低时,弯矩增大,而在水平荷载较高时,桩身弯矩反而小于ΔT=0℃的工况,其主要原因为：桩土共同收缩时,桩端及桩周约束减弱,基桩变形主要表现为绕旋转中心偏转。此外,不难看出最大弯矩位置均处于3 m深度（0.36L）附近,即温差对最大弯矩位置影响较小,工程设计时应考虑在该深度附近适当提高配筋率,以避免桩体发生破坏。

为定量分析温差ΔT对桩身最大弯矩M_max的影响,定义$\alpha$(M_max)=ΔM_max/M_{max(ΔT=0℃)}为桩身最大弯矩增大系数,其中,ΔM_max为水平荷载相同时,某温差工况下桩身最大弯矩与ΔT=0℃工况下桩身最大弯矩M_{max(ΔT=0℃)}的差值,探究桩身最大弯矩增大系数随温差的变化规律,如图10所示。由图10可知,当ΔT>0℃时,桩身最大弯矩增大系数与温差的关系近似呈线性变化,且其斜率随水平荷载增大而减小,其主要原因为水平荷载增大时,温度影响相对减弱,故对图中4条曲线ΔT>0℃部分进行拟合,可得其统一表达式为

图  10  桩身最大弯矩增大系数-温差曲线

Figure  10.  Curves of increment coefficient of maximum bending moment of energy pile-temperature difference

下载: 全尺寸图片 幻灯片

其余ΔT>0℃温差与荷载工况下桩顶沉降可通过式（3）计算出桩身最大弯矩增大系数$\alpha$(M_max),并结合ΔT=0℃工况下桩身最大弯矩进行求解,但ΔT<0℃工况下桩身最大弯矩受温差影响较小。

4.   结论

本文通过与模型试验及二维数值结果的对比,验证了考虑热-力耦合作用的能量桩三维有限元分析方法合理性,针对倾斜荷载下能量桩受力变形特性这一问题,探讨了考虑温差影响的桩顶沉降、桩侧摩阻力、水平位移及桩身弯矩变化规律,主要得出4点结论。

（1）当ΔT<0℃时,桩顶沉降w₀明显大于ΔT≥0℃工况；而当ΔT>0℃时,在竖向荷载较小时,桩顶甚至出现向上的位移,但其w₀–P曲线随着竖向荷载的增大逐渐趋近ΔT=0℃工况。能量桩在正常使用状态下,将在桩顶及桩端产生不可忽视的附加温度应力（压应力或拉应力）,因而,在能量桩设计时,应对附加温度应力（尤其是拉应力）进行初步估算,使其低于混凝土抗压及抗拉强度,并适当提高桩顶与桩端钢筋配筋率。

（2）当ΔT>0℃时,桩身上、中部桩侧摩阻力为负,下部桩侧摩阻力为正；而当ΔT<0℃时,桩身上、中部桩侧摩阻力为正,而下部桩侧摩阻力为负。能量桩桩侧摩阻力分布规律及大小受温差影响较大,在能量桩设计时,应充分考虑负摩阻力的影响,如适当提高通长钢筋配筋率、选址时尽量避开高压缩性土层、避免大面积堆载等。

（3）升温和降温工况下桩顶水平位移均有所增大,但其原因不尽相同；工程实际中,若要弱化温差的影响,可考虑整体加强桩周、桩顶及桩端的侧向约束,如嵌入岩层或扩大下部桩身半径等。

（4）各温差工况下桩身弯矩随深度的变化趋势大致相同。当ΔT>0℃时,桩身弯矩始终表现出增大的趋势；而当ΔT<0℃时,桩端及桩周约束减弱,基桩变形主要表现为绕旋转中心偏转,桩身弯矩小于ΔT=0℃的工况。此外,最大弯矩位置受温差影响较小,均处于3 m深度（0.36L）附近,工程设计时应考虑在该深度附近适当提高配筋率,以避免桩体发生破坏。