竖向力-温度荷载-扭矩加载路径下单桩承载变形研究

2023-11-22陈秋怡黄中正陈朝棋欧孝夺

工程力学 2023年11期

江杰，陈秋怡，黄中正，陈朝棋，欧孝夺

(1.广西大学工程防灾与结构安全教育部重点实验室，广西，南宁 530004；2.广西大学土木建筑工程学院，广西，南宁 530004；3.广西大学防灾减灾与工程安全重点实验室，广西，南宁 530004)

近年来，基于传统地源热泵技术发展起来的能量桩逐渐受到国内外学者的广泛关注[1-3]。对于高层建筑等大型结构物使用的能量桩而言，除了受到竖向荷载、温度荷载作用外，可能还会承受水平荷载、扭矩荷载和动荷载。鉴于问题的复杂性，本文先讨论竖向力、温度和扭矩作用下能量桩的承载变形研究，以期更好地推广能量桩的应用。

对于竖向荷载下能量桩的受力变形特性，国内外学者进行了相关的试验。LALOUI 等[4]和BOURNEWEBB 等[5]开展了能量桩现场原位实验，研究温度变化对桩身力学行为的影响；路宏伟等[6]和蒋刚等[7]开展竖向力和温度荷载作用下的现场试验，揭示了桩体温度变化和桩顶竖向荷载会引起能量桩荷载传递机制和桩身内力的变化；方鹏飞等[8]结合现场试验探究升温工况下能量桩的工作机制，重点研究了能量桩桩身内力的变化规律；王成龙等[9]认为桩顶和桩端约束受温度荷载的影响，并通过开展相应的模型试验，分析不同约束条件下桩身位移和应力变化规律，并进一步研究了温度零点随约束不同的变化规律。鉴于试验费用高昂，试验条件难以控制，许多学者开始用理论方法研究能量桩的受力变形特性。KNELLWOLF等[10]将荷载传递法应用于能量桩，采用弹簧模拟温度荷载下上部结构与桩的相互作用，构造了竖向力和温度荷载作用下的荷载传递方程；PASTEN等[11]进一步将温度变形的影响考虑在桩段压缩量中，分析了温度荷载作用下能量桩位移变化规律；徐新丽等[12]基于弹性有限单元分析模型，通过对桩身和桩周土体进行受力分析，将温度荷载简化为单元内力中，最后结合静力平衡，建立了竖向力和温度荷载作用下能量桩桩身平衡方程，并结合London[4]试验和OUYANG 等[13]提出的混合传递法进行了对比验证。费康等[14]引入双曲线荷载传递函数，将温度荷载作用考虑到单元压缩量中，结合单元增量平衡方程，揭示了任意竖向力和温度荷载作用下的桩身变形、桩身内力分布。

综上，针对竖向力和温度荷载作用下能量桩承载变形特性，在试验和理论方面研究取得较多成果。但目前针对多向荷载作用下能量桩的研究仍较少，缺乏足够的理论研究，限制其进一步的推广应用。因此，本文对竖向力→温度荷载→扭矩加载路径下能量桩的承载变形展开研究。根据能量桩的工作条件，依次施加桩顶竖向力、温度荷载和桩顶扭矩。基于荷载传递法和边界单元法提出竖向力→温度荷载→扭矩加载路径下的能量桩单桩计算方法，并通过与已有试验和ABAQUS有限元结果进行验证，最后分析竖向荷载、长径比和温度对能量桩承载变形的影响。

1 荷载传递模型及假定

本文分析基于以下假定：

1) 竖向力和温度荷载作用下桩的径向变形较小[15]，故计算分析中只考虑桩的轴向变形。

2) 桩为线弹性等截面桩，桩-土界面处的荷载传递函数均符合双曲线形式，如图1 所示，桩侧和桩端荷载传递函数分别表示为：

图1 能量桩桩-土界面模型Fig.1 Pile-soil interface model of energy pile

式中：τ和τ' 分别为桩侧剪应力和桩端阻力； Δ和Δb分别为桩-土相对位移和桩端位移；a、b、a＇和b＇均为模型参数。对桩侧，根据CAO 等[16]建议确定，a=Δu/τu，b=Rsf/τu，其中： Δu为桩体达极限摩阻力对应的极限位移，FLEMING 等[17]建议可取0.5%～2.0%D，D为桩的直径；Rsf为桩侧破坏比，根据费康等[11]建议取为1.0；τu为桩-土界面的破坏剪应力，τu=K0σtanφ，其中K0为静止土压力系数，K0=1-sinφ，φ为桩-土界面摩擦角；σ 为竖向应力。对桩端，根据RANDOLPH[18]建议确定，a'=其中：为桩端土的剪切模量；ν' 为桩端土的泊松比；为桩端破坏比，取为1.0；τ'u为桩端界面阻力破坏值，根据极限平衡理论确定，τ'u=ζcNc+ζqγhNq，其中 ζc和ζq为桩端为方形、圆形时的形状系数，Nc和Nq为承载力系数。

3) 桩周土体性质、桩侧和桩端荷载传递函数不受温度影响。

4) 取温度增加、桩体位移向下和桩侧阻力向上为正。

2 理论分析及推导

2.1 竖向力和温度荷载作用下单桩计算

将桩长为L的能量桩划分为N个单元，对单元和节点自上向下编号，建立竖向力和温度荷载作用下单桩承载变形特性的简化分析方法，如图2所示。

图2 能量桩沿桩长离散示意图Fig.2 Discrete diagram of energy pile along pile length

取任意桩体单元分析，由增量平衡条件，有：

式中： ΔPi为单元i上截面的轴力增量； ΔPi+1为单元i下截面的轴力增量； ΔQi为单元i桩侧阻力的增量；ksi为单元i中点处桩侧抗剪切刚度，即荷载传递函数的切线斜率； ΔLi为单元i长度；D为桩体直径； ΔXi为节点截面i的位移增量； ΔXi+1为节点截面i+1的位移增量。

式中， Δsi为单元的压缩量。考虑竖向力和温度荷载的相互作用，有：

式中：E为桩身弹性模量；A为单元截面面积；为自由膨胀量，根据线热膨胀理论，=αT·ΔT·ΔLi， αT为桩身混凝土线膨胀系数， ΔT为温度增量。

联立式(3)～式(6)，分离变量，有：

写成矩阵形式有：

式中：ki为桩身单元i刚度矩阵。其中：

结合单元上、下截面连续性，将每个桩身单元的刚度矩阵进行组装，得到竖向力和温度荷载作用下桩身整体矩阵通式为：

接着考虑桩顶边界第一个单元的轴力增量及桩端边界第N+1 节点反力得到式(11)、式(12)。

桩顶节点1：

桩端节点N+1：

式中，kb为桩端土体的抗压刚度系数。将式(12)代入式(8)，有：

综上，可得到竖向力和温度荷载作用下能量桩的控制方程：

其中，桩的总刚度矩阵K为：

对方程组求解可得到桩体节点位移增量 ΔX。求解过程如下：① 将荷载分为多级荷载增量，从第一个荷载增量 ΔF1，先施加荷载增量的1/2，即ΔF1/2，根据桩侧、桩端初始刚度建立桩身刚度矩阵，分别求出 ΔF1/2的节点位移；② 根据上一增量结束阶段的节点位移、应力建立刚度矩阵，求出增量中点位移；③ 用新求出的桩身节点节点位移求出增量 ΔF1的中点刚度矩阵，接着由式(14)求出第一级增量下的位移、应力；④ 重复步骤②～步骤③，直到最后一个增量，求出桩体节点位移。

根据能量桩的实际运行情况，依次进行竖向力加载和温度加载。竖向力加载中桩顶视为自由，温度荷载加载过程中上部结构对桩顶的约束作用，采用刚度系数为kt对节点1 的刚度矩阵进行修正，得到式(15)：

2.2 竖向力→温度荷载→扭矩加载路径下单桩计算

2.2.1 环向极限摩阻力的确定

为求解竖向力→温度荷载→扭矩加载路径能量桩承载变形特性，根据江杰等[19]和邹新军等[20]，通过边界元法，单元i表面作用有竖向力、温度产生的桩侧竖向摩阻力τv(i)和扭矩作用下的环向摩阻力τt(i)，如图3 所示，假设桩身摩阻力满足式(16)，有：

图3 桩周土体剪切作用Fig.3 Shear action of soil around piles

竖向力→温度加载结束后，扭矩作用下的环向极限摩阻力为：

式中：τf(i) 为单元i的极限摩阻力；τv(i) 为单元i的竖向摩阻力；τtf(i) 为单元i的环向极限摩阻力。

2.2.2 竖向力→温度荷载→扭矩加载路径下单桩计算

扭矩作用下，桩单元的静力平衡条件为：

式中：T(z)为z处扭矩；τt(i) 为深度z处桩侧环向摩阻力。

扭矩作用下桩的变形条件为：

式中： θ(z) 为扭转角；Gp为桩单元的剪切模量；Jp为桩单元的极惯性矩。

根据江杰等[21]，扭转角与扭转位移的关系可以表示为：

式中，st(z) 为深度z处的扭转位移。

桩顶、桩端边界条件分别如下：

桩顶：

桩端边界条件引入POULOS[22]的建议，有：

式中：T1为桩顶扭矩； θi为第i单元的桩身扭转角为桩端扭矩。

采用有限差分进行求解：① 假设桩身扭转角θ为任意非零矩阵，根据式(17)、式(20)和式(22)求得桩身扭转矩阵；② 由式(25)求得到桩身扭转角 θk；③ 用 θk求出新的桩身扭转刚度矩阵，由式(25)求得到新的 θk+1。当 |θk-θk+1|小于限定值时，可得到桩身扭转角θ ，若 |θk-θk+1|不满足要求则重复步骤②～步骤③直至满足限定值。

2.3 数值实现

本文计算依次进行竖向力加载、温度荷载加载和扭矩加载，有：

1) 竖向力加载

令 ΔT=0，将竖向力划分为若干增量步，根据土层参数和桩身参数得到初始刚度矩阵，并采用中点增量法对式(14)进行求解，得到桩体节点位移增量。

2) 温度荷载加载

以竖向力加载结束的桩身节点位移作为温度加载的初始加载状态，将温度荷载划分为若干增量步，求解过程与竖向力加载大致相同，最后得到竖向力→温度加载下的桩基响应。

3) 扭矩加载

计算竖向力加载和温度荷载加载结束后的桩侧竖向摩阻力，通过式(17)计算桩侧环向极限摩阻力，取深度z处的微段桩段 dz进行分析，建立桩身扭矩角整体方程，利用有限差分法，对式 (25)进行求解，通过控制 |θk-θk+1|迭代误差得到桩身节点扭转角，最后输出竖向力→温度荷载→扭矩加载路径下能量桩的桩基响应。

3 方法的验证

3.1 纯温度荷载作用下单桩计算方法验证

1) 模型简介

本文以文献[23]中的能量桩为例，桩径D=0.88 m，桩长L=19.6 m。通过ABAQUS 有限元软件进行模拟，在桩-土界面设置接触对模拟实际桩-土界面的相互作用。桩-土界面的摩擦角根据试验取为25°，桩身和桩周砂土分别采用弹性本构模型和Mohr-Coulomb 理想弹塑性模型。桩-土其他材料参数见表1。

表1 桩-土材料参数Table 1 Pile and soil parameters

2) 边界条件及网格划分

根据工程实际，土体侧面约束径向位移，底部和顶部视为固定边界和自由面边界条件。

通过自动平衡法来实现地应力平衡，在地应力平衡分析步，仅对完整的土体模型进行地应力平衡；在静力分析步，将桩单元激活，从而实现地应力平衡。在静力分析步施加竖向力，竖向力加载结束后，导入温度场模拟能量桩受热状态，最后达到指定温度后在桩顶施加扭矩。桩体和土体均采用C3D8 网格类型，桩体模型通过结构化网格技术进行划分，总共5760 个桩体网格单元；土体模型首先对土体定义Partition，最后结合扫掠技术划分，总共45 000 个土体网格单元，如图4 所示。

图4 三维数值模型Fig.4 3D numerical model

3) 桩身轴力分布

纯升温工况下，桩身的轴力在数值上沿深度呈现先增大、后减小的规律，大致在桩身0.6L处达到最大值。图5 为不同温度增量下沿桩长方向的桩身轴力图，随着温度增量的增大，桩身轴力逐渐增大，其中桩顶和桩端的轴力较小，桩身中下部轴力较大。可以看出，本文方法计算结果、有限元结果与NG 等[23]的试验数据都具有较好的一致性。

图5 温度荷载下桩身轴力分布Fig.5 Axial force distribution of pile under thermal load

4) 桩侧竖向摩阻力分布

纯升温工况下，桩身体积变化导致与桩周土体将产生相对位移，从而桩侧竖向摩阻力发生改变。图6 为不同温度增量作用下沿桩长方向的桩身竖向摩阻力分布图。

图6 温度荷载下桩侧摩阻力分布Fig.6 Distribution of pile shaft resistance under thermal load

由图6 可知，计算值与模拟值有很好的一致性，桩侧摩阻力沿深度的变化规律与龚建清等[24]的模拟趋势一致。中性点(桩侧摩阻力为0 的点)位置大约在桩长0.6L处，与BOURNE-WEBB 等[5]现场实测的中性点位置相近，桩侧竖向摩阻力数值上沿桩长方向呈现上部分先增大后减小至0，下部分从0 开始反向增大的分布，这与LAOUI 等[4]研究成果一致。本文算例工况采用NG 等[23]的进行计算，中密砂土层，弹性模量较小，对桩端的约束较小，故中性点位置在0.6L处，而LAOUI等[4]的桩端持力层刚度大，可提供更大的桩端约束，因此其中性点更接近桩端。此外，可以看出，纯温度荷载作用下 (ΔT>0)，导致桩身上半部分的摩阻力为负值， ΔT= 3 0 ℃相对于 ΔT=15 ℃，能量桩最大侧摩阻力(绝对值)有所增大。

将纯温度加载工况下的计算值与相关试验数据和有限元对比，证明了方法的准确性。

3.2 竖向力→温度荷载加载路径下单桩计算方法验证

本文的计算理论不仅可以分析竖向力→温度荷载→扭矩加载路径下的能量桩承载变形特性，也可以用来分析竖向力→温度荷载加载路径下能量桩计算。根据蒋刚等[7]试验参数，重新建立数值模型，将竖向力→温度荷载加载路径下单桩计算值、模拟值与蒋刚等[7]试验值进行对比。图7为桩身的轴力随深度变化曲线，对比单一竖向力加载、竖向力→温度荷载工况可知，温度荷载会改变桩身轴力分布。由图7 可知，本文计算值、模拟值与试验值接近，因此可证明本文计算方法的正确性。

图7 竖向力→温度荷载加载路径桩身轴力分布曲线Fig.7 Axial force distribution curve under vertical force →thermal loading path

3.3 竖向力→温度荷载→扭矩加载路径下单桩计算方法验证

本文基于NG 等[23]的试验，对温度增量为30 ℃的单桩施加扭矩，通过与有限元进行对比，验证竖向力→温度荷载→扭矩加载路径下单桩承载变形计算方法的正确性。

3.3.1 桩侧环向摩阻力分布

为研究竖向力→温度荷载→扭矩加载路径下，桩顶扭矩对能量桩环向摩阻力的影响。对竖向力→温度荷载结束后施加扭矩后，为了验证本文理论方法的正确性，通过有限元分别对桩顶施加扭矩100 kN·m、200 kN·m、300 kN·m、400 kN·m，提取扭矩作用后的桩身环向摩阻力后与计算模拟结果对比如图8 所示，4 组计算值与数值模拟结果均比较接近，证明了本文方法在计算竖向力→温度荷载→扭矩加载路径下桩侧环向摩阻力的正确性。

图8 桩侧环向摩阻力分布Fig.8 Distribution of circumferential friction resistance on pile side

由图8 可知，当桩顶扭矩分别为100 kN·m、200 kN·m、300 kN·m、400 kN·m 时，桩身环向摩阻力随深度呈现线性增大的趋势，总体上，竖向力→温度荷载→扭矩加载路径下桩身环向摩阻力会随着桩顶扭矩的增大而增大，当桩顶扭矩增大到一定范围后，环向摩阻力的增长较之前缓慢。原因是，扭矩较小时，桩-土界面进入弹性阶段，温度荷载引起土体的抵抗作用发挥明显，从而导致环向摩阻力随深度增大；当扭矩增大到一定数值时，环向摩阻力趋势趋于一致，表明桩顶扭矩起主导作用。

3.3.2 桩顶扭矩-扭转角变化曲线

对于竖向力→温度荷载→扭矩加载路径下单桩计算，通过确定桩侧极限环向摩阻力，并通过式(24)考虑桩端边界条件，最后得到T-θ 曲线。取竖向荷载P=1000 kN 时，研究扭转角随不同桩顶扭矩作用下的变化。整体上来看，本文计算得到的曲线和有限元结果较为吻合，证明了竖向力→温度荷载→扭矩加载路径下单桩承载变形计算方法具有较好的准确性。从图9 可看出，ΔT=30 ℃工况下单桩的极限扭矩[25]为645 kN·m，对比ΔT=0 ℃工况可知，能量桩的极限扭矩降低了20.37%，说明温度荷载对能量桩的极限扭矩有影响。

图9 T -θ曲线对比Fig.9 Comparison of T -θ curves

4 参数分析

为研究竖向力→温度荷载→扭矩加载路径下单桩的影响因素，分别对竖向力、长径比和温度增量进行分析。

4.1 竖向力的影响

为研究不同竖向荷载(根据NG 等[23]的试验，取温度荷载增量为30 ℃)对能量桩荷载传递特性的影响，分别施加25%Pu、50%Pu、75%Pu的竖向力(Pu为竖向极限荷载)。图10 给出3 种荷载工况下，桩顶扭矩角随扭矩的发展情况。可见，温度增量和桩顶扭矩为定值，随着竖向力的增加，25%Pu、50%Pu、75%Pu对应的极限扭矩分别为680 kN·m、590 kN·m 和515 kN·m，极限扭矩降低幅度为26.2%，说明了竖向力→温度荷载→扭矩加载路径下，在达到竖向极限荷载的过程中，由竖向力产生的桩侧竖向摩阻力会不断逼近极限值，进而影响桩侧极限环向摩阻力，最终导致单桩极限扭矩受到影响，此结果验证了式(17)的正确性。

图10 三种竖向力工况下 T-θ曲线对比Fig.10 Comparison of T -θ curves under three vertical load conditions

图11 为25%Pu、50%Pu、75%Pu(温度增量30 ℃，扭矩350 kN·m)三种工况下的桩身扭转角随深度的分布情况。由图11 可知，三种竖向力工况下桩身扭转角存在较大的差异。25%Pu、50%Pu和75%Pu作用下桩顶桩身扭转角分别为0.020、0.023和0.027，而桩顶扭转角最大增幅约35%。因此，对于竖向力→温度荷载→扭矩加载路径下，随着竖向力的增大，桩身变形会受到较大的影响。因此工程桩基中的能量桩需适当控制桩顶竖向力水平。

图11 三种竖向力工况下桩身扭转角分布情况Fig.11 Distribution of pile torque angle under three vertical load conditions

4.2 长径比的影响

为探究桩身长径比对承载力的影响规律，选取3 组不同长径比的桩进行分析，获得相应的桩身承载力(桩顶所能承受竖向力与扭矩极值)包络图如图12 所示。

图12 不同长径比L/D 对桩顶承载力包络线的影响Fig.12 Effect of different ratio L/D on the envelopes of pile top bearing capacity

由图12 可以看出，相同温度增量下，当长径比L/D=10 时（保持桩径不变，改变桩身长度），桩身承载能力较小；当长径比L/D=30 时，对应的扭矩极限承载力相对于L/D=20 增加约 1.1 倍，表明随着桩身长径比的增加，桩身承载力得到提高，承载力包络线有所外扩。

4.3 温度增量的影响

为研究温度增量对能量桩承载变形特性的影响。图13 给出了温度增量分别为10 ℃、20 ℃、30 ℃、40 ℃工况下，桩顶扭矩-扭矩角分布情况。可以看出，增大温度增量，单桩抗扭能力会受到影响。

图13 不同温度增量下 T -θ曲线对比Fig.13 Comparison of T -θ curves under different temperature increments

图14 为不同温度增量(扭矩为500 kN·m)工况下的桩身扭转角分布情况。由图14 可知，四种温度工况下的沿深度分布的桩身扭转角存在较大的差异。在竖向力→温度荷载→扭矩加载路径下，随着温度增量的增大，桩身变形会受到影响。对比温度增量10 ℃、20 ℃、30 ℃、40 ℃的桩身扭转角沿深度的分布规律，可知桩端扭转角受温度影响较大，说明升温引起的桩身体积变化，会导致桩端阻力受影响。此外还可以看出，大约0.6L以上扭转角较大，0.6L以下扭转角较小，原因是随着温度荷载的增大，桩端约束有所提高，故靠近桩端附近的变形较小。因此，竖向力→温度荷载→扭矩加载路径下的单桩变形主要发生在桩身0.6L以上，故在能量桩施工过程中要注意对浅层地基的保护。