Study on bearing capacity characteristics of energy pile based on three-dimensional modeling*
随着时代的发展，人们逐渐重视清洁能源的利用.地源热泵是一种通过在地下埋设换热管的方式，利用浅层地热资源为建筑供热或制冷的环保型空调系统.能量桩是基于传统地源热泵技术引申出的新型建筑节能技术，既能承担上部建筑物荷载，又因其内置换热管兼具地源热泵换热器的作用.目前，大部分采用地源热泵的工程是建筑物外钻孔埋管技术，占用场地较大，成本较高［1］，而能量桩技术具有节约建筑用地、工期短、造价低及换热效率高等优点.1984年，能量桩技术在奥地利被首次提出，之后很快在德国、瑞士等西方国家推广，并在工程应用中取得众多成功案例［2-3］.近年来，能量桩技术也被应用在我国南京市朗诗国际街区和天津市梅江综合办公楼等众多工程项目中.
能量桩在运行过程中，桩-土之间的热交换不仅会引起桩体周围土体温度场的变化，还会导致桩身产生附加热应力和变形，从而影响桩基承载性能.针对这一问题，相关学者开展了系列研究：Knellwolf等［4］基于荷载传递法，提出一种土工数值分析方法，用以评估温度循环对能量桩的影响，但未考虑温度场对桩体周围土的影响；Yavari等［5］开展模型试验，对饱和黏土中能量桩在热-力耦合作用下的力学性能进行研究，表明温度循环后桩顶产生的不可逆沉降随工作载荷的增加而逐渐增大；桂树强和程晓辉［6］针对能量桩运行过程中的结构响应开展现场原位试验，得出对桩体施加温度荷载后，桩身产生的附加温度荷载大小受桩端约束的影响等结论；Wu等［7］研究了饱和黏土中能量桩与邻桩和承台之间的相互作用，显示邻桩和承台对能量桩的沉降起到约束作用；方鹏飞等［8］通过室内模型试验和现场原位试验，研究了温度变化对能量桩承载性能的影响，表明温度循环在桩身产生附加热应力，进而对桩体侧摩阻力有较大影响；Stewart和 Mccartney［9］基于离心模块化方法，研究了能量桩在温度循环时的瞬态热力响应，显示桩端附近的热应力更大，多次温度循环后桩顶产生沉降.
在开展模型试验或现场原位试验过程中，会受到试验成本、工期和试验仪器等条件的约束，研究能量桩在运行过程中的热力学特性有一定局限性.随着计算机技术的进步，数值模拟更加便捷化，因其不受季节和地域等因素的限制被广泛运用于科学研究中.近年来，数值模拟的方法频繁地被运用到能量桩的研究中：Morino和 Oka［10］率先在数值模拟中引入有限差分法分析能量桩的传热性能，并定义了桩基埋管换热器的概念；Gashti等［11］基于有限元数值模拟对复合能量桩基础的结构和岩土阻力进行研究，认为当降温时桩身会产生附加应力；郝耀虎等［12］基于室内模型试验与数值模拟，开展冷-热温度循环荷载作用下桩基在不同桩端约束条件下的桩身应力和位移特性研究，表明桩身位移变化规律受温度与桩端约束共同影响；王成龙等［13］基于室内模型试验和数值模拟，对饱和砂土中能量桩承载特性和传热特性进行研究，显示桩顶沉降随温度循环次数的增加不断累积；刘干斌等［14］开展模型试验和数值模拟研究，分析饱和黏土中不同温度工况下能量桩承载力特性，得出升温后地基发生热固结现象表现为沉降变形，能量桩的单桩极限承载力随温度的升高而增大的结论；费康等［15］基于ABAQUS数值分析软件建立有限元模型，分析正常固结黏土地基中能量桩长期运行的性能，显示温度循环次数越多，桩顶累积沉降越大；王中豪［16］运用COMSOL软件建立能量桩的三维数值模型，分析其换热性能的影响因素，显示螺旋型埋管能量桩换热能力最强；翁凯文［17］通过试验和数值模拟针对不同界面温度和循环模式对能量桩桩-土界面摩擦性状的影响，表明能量桩桩-土界面的温度与界面抗剪强度呈正相关趋势.
上述成果对研究能量桩的变形特性和承载特性起到了指导作用，但针对温度荷载对能量桩及其周围土体影响的研究仍相对较少.为了进一步确定能量桩工作时的承载力特性，需要进行大量的研究，数值模拟能够较好地反映能量桩的各项参数而不受物理条件的约束，其对于能量桩换热系统的研究有重大意义.因此，本文开展饱和黏土中能量桩在不同工况下的数值模拟研究.
COMSOL Multiphysics软件基于偏微分方程（partial differential equation，PDE）建模，内置结构力学、流体流动、传热和数学等多种模块，可快捷定义和求解多种物理场耦合的问题.COMSOL软件架构开放，操作界面简洁明了，初始条件、边界条件等可自由定义，有丰富的数据处理方式［18-19］.运用COMSOL软件建立能量桩数值模型，可以有效分析其运行过程中的热-力响应特性.
为优化模型求解的计算过程，在建立能量桩数值模型时需要引入的假设条件如下：（1）桩身材料和桩体周围土体各向同性，模型各部分的热物性参数保持恒定；（2）不考虑土壤中渗流影响，无地下水流动换热；（3）将U型换热管及其中的换热液体简化为管道流；（4）建立理想线弹性能量桩模型；（5）模型各处的初始温度相同且均匀分布.
模型传热过程主要包括换热管与混凝土的耦合传热以及桩身与土体传热，由于不考虑土壤中渗流及地下水流动换热的影响，桩体周围土体主要考虑热传导和部分热辐射.传热控制方程为：
式中：ρ为材料密度（kg/m3）；CP为定压比热容（J·kg-1·℃-1），i=p 时，CPi表示桩身的比热容，i=s时，CPi表示饱和黏土的比热容；u为流体速率（m/s）；∇为哈密顿算子符号；T为温度（℃）；q为热通量；Q表示热源；Qted表示热弹性阻尼；keq为热传导率（W·m-1·℃-1）.
换热管的导热效率主要与管内液体的流速有关，换热管的材料沿管身的换热能力在换热系数计算时应当同时考虑，通过 Dittus-Boelter公式［20］换算可得
式中：Kpipe为换热管管壁的导热系数；
分别为换热管的外径和内径；k为接触面对流传热系数，其数值主要与换热液体的流速及换热液体的物理性质有关，其表达式为
式中：Re为换热液体的雷诺系数；Pr为普朗特数；Kf为流体热导率.
桩内换热管管壁处的耦合传热控制方程表示为
式中：p为管内换热液体的压强（N/m2）；A为换热管截面面积（m2）；CP为换热液体的比热容（J·kg-1·℃-1）；fD为临界流动的达西摩擦系数；dh为水力直径（m）；u为截面平均速度（m/s）；et为t时刻水流的切向速度（m/s）；Q为热强度（W/m3）.
1.3.1 基本参数
模型桩直径D1=84 mm，桩长L=500 mm，桩周围土体直径D2=550 mm，高度H=500 mm，桩顶高于土体表面100 mm，模型尺寸如图1所示.
图1 模型示意
桩内埋设单U形换热管，管中换热液体为水，设置流速度为22 L/min，其各项参数从COMSOL Multiphysics材料库中导入.通过桩-土主要物理参数计算得饱和黏土地基的单桩竖向承载力特征值为0.5 kN，桩、土基本材料参数如表1所示.数值模拟设置2种工况，分别为桩顶零荷载和桩顶设置0.5 kN的工作荷载，均将换热液体加热至70℃维持24 h，之后降温至5℃，维持5 h，重复循环3次，共87 h.在COMSOL软件中将温度循环设定为分段函数，区间设置0～24 h，对应函数值是70℃；24～29 h对应函数值是5℃.选择外推为周期性，无平滑处理.
表1 材料主要物理参数
images/BZ_36_236_414_2244_471.png混凝土饱和黏土2 400 1 900 1.28 1.50 1×10-5 1×10-5 0.35 0.50 3×104 15
1.3.2 边界条件及网格划分
设定模型各处为线弹性材料，土体表面可自由膨胀和收缩，四周采用辊支撑，底部采用固定约束，模型四周和底部设置为热绝缘，其余部分设置热通量.模拟过程考虑地应力影响，对桩体和土体施加重力荷载.三维数值模型的网格最小单元越小，计算所得的各项数据越可靠，但计算所需的时间更长，对运行软件的设备性能要求更高.因此，需要根据模拟的具体情况选择适当的网格密度.本模型对换热管、桩身和土体分别采取不同的网格，对换热管采取精度较高的边单元，最小单元长度为0.12 mm；桩身和土体选用自由四面体网格，最小单元长度分别为0.9和2.4 mm，划分结果如图2所示.
图2 模型网格划分
（a）土体网格划分；（b）桩体网格划分
选取零荷载工况下，温度循环过程中不同时刻桩-土温度场分布如图3所示.升温阶段，由于换热管作为热源，其温度始终最高，高温区以换热管为中心向外扩散，桩周围土体高温区大致沿桩身呈长轴为400 mm、短轴为200 mm的椭圆形对称分布.由于土壤热阻的影响，土体温度随径向距离增大而减小.降温阶段，桩体温度骤降，土体表面温度较为恒定.随着循环次数的增加，桩-土温度逐渐降低，考虑是降温阶段使桩身温度低于初始温度所致.
图3 不同时刻模型温度场分布
（a）t=24 h；（b）t=29 h；（c）t=53 h；（d）t=58 h；（e）t=82 h；（f）t=87 h
温度循环为82和87 h的桩顶和土体表面位移变化规律如图4所示，考虑到位移量较小，故选择1∶100比例因子放大变形.由图4可知，加热阶段桩体膨胀，有、无工作荷载作用时桩顶产生不同程度的上升位移.零荷载作用下，桩顶最大上升位移为0.503 2 mm；工作荷载作用下，桩顶最大上升位移为0.151 3 mm，桩顶位移的差别表明工作荷载的约束作用会抑制桩顶膨胀.制冷阶段，零荷载作用下桩体受冷收缩产生的沉降未将受热膨胀产生的上升位移完全抵消；工作荷载作用下桩顶表现为少量沉降，沉降量为0.038 1 mm，可能是由于降温时桩体受冷收缩，其与土体之间的水平土压力减小.加热过程中，受桩-土温度场的影响，土体膨胀量沿深度递减；降温时桩周围土由膨胀转为沉降，考虑是由于降温导致超静孔隙水压力的消散，使土体产生固结沉降所致.
图4 桩身和土体位移
（a）零荷载加热；（b）零荷载制冷；（c）工作荷载加热；（d）工作荷载制冷
综上表明，当桩体受热膨胀时，工作荷载会抑制桩顶的上升位移；当桩体制冷收缩时，工作荷载会引起桩顶沉降.能量桩在实际应用中需长期运行并保证其稳定性，而桩顶位移尤其是沉降位移是影响结构安全的重要因素.因此，应重视温度循环导致的桩顶位移，并且设计阶段需考虑到黏土产生变形对桩基承载力的影响.
零荷载和工作荷载2种工况下，桩端压应力随时间变化的曲线如图5所示.加热时，桩体受热膨胀，桩端压应力在4 h内迅速升高，可能是由于桩体温度未达到稳定状态；制冷时，桩体受冷收缩，桩端压应力逐渐减小.加热阶段结束后，零荷载和工作荷载作用下桩端的最大压应力分别为0.128 7和0.222 5 MPa；制冷阶段结束后，桩端最小压应力分别为0.017 6和0.119 0 MPa.工作荷载作用下，加热结束时的最大压应力约为零荷载作用下的1.73倍，制冷结束时的桩端压应力约为零荷载作用下的6.8倍.随着循环次数的增加，桩端最大压应力逐渐减小，考虑是由于桩周围土体在多次冷热交替后发生热固结作用所致，土体强度增大，在实际工程运用中，应着重关注桩侧土的变化及桩端土的性质.
图5 桩端压应力随时间变化曲线
COMSOL Multiphysics岩土力学模块是一个特殊附加模块，主要用于模拟岩土工程应用，模块中设置专用接口以研究岩土与人造结构之间的交互作用问题.沿桩身土体表面至桩端每50 mm深度取1个测点，共9个测点，运用岩土工程模块计算各测点的桩身应力和桩体侧摩阻力.
选取2种工况24和29 h的桩体侧摩阻力绘制其沿深度分布如图6所示.定义侧摩阻力向上为正，向下为负.通过零荷载和工作荷载作用下的桩侧摩阻力对比，分析工作荷载作用对桩侧摩阻力的影响.桩体受热时，桩体上部侧摩阻力为负，下部侧摩阻力为正，制冷时恰好相反.第24 h时刻，工作荷载和零荷载作用下最大侧摩阻力值分别为0.016 3和0.025 8 MPa；第29 h时刻，二者最大侧摩阻力值分别为0.029 9和0.047 3 MPa.零荷载作用下，位移零点位于桩体中部偏下位置；工作荷载作用下，位移零点上移至桩中部偏上位置.由于工作荷载作用使桩身整体下沉，导致产生负摩阻力的区域变小.零荷载作用下，加热时桩顶无约束可自由膨胀，在桩身上部产生的负摩阻力最大值为0.025 8 MPa.由于工作荷载的约束作用，此时桩体上部最大负摩阻力仅为0.012 6 MPa.
图6 侧摩阻力沿深度分布
选取零荷载作用下3次温度循环加热和制冷结束时的桩侧摩阻力数据绘制其沿深度分布图（图7），以分析温度循环对桩体侧摩阻力的影响.多次温度循环位移零点始终保持在桩体中部偏下位置.3次循环加热时，最大侧摩阻力值依次为0.025 8、0.028 8和0.030 8 MPa；制冷时最大侧摩阻力值依次为0.010 4、0.016 0和0.016 5 MPa.随着循环次数的增加，桩侧摩阻力逐渐增大，增量逐渐减小，可能是多次温度循环使桩体周围土产生固结变形，对桩体的约束作用增强.
图7 零荷载3次温度循环侧摩阻力沿深度分布
运用COMSOL Multiphysics软件建立饱和黏土中能量桩的三维数值模型，分析桩-土温度场、桩顶位移和桩基承载力，得到以下几点结论：
（1）升温阶段，桩-土温度场以换热管为中心向外扩散，土体温度随径向距离增大而减小；降温阶段，桩体温度骤降，土体表面温度较为恒定.加热时，零荷载作用下桩顶膨胀量较大，工作荷载的约束作用会限制桩顶膨胀；制冷时，工作荷载会引起桩顶0.038 1 mm（0.46% D1）的少量沉降，靠近桩身的土体由加热时的膨胀转为收缩.
（2）桩端压应力在加热时增大，制冷时减小.工作荷载作用会引起较大的桩端阻力，加热结束时的最大压应力约为零荷载作用下的1.73倍，制冷结束时约为零荷载作用下的6.8倍.随循环次数的增加，桩端最大压应力逐渐减小.
（3）加热后桩体上部产生负侧摩阻力，下部产生正侧摩阻力，而制冷时恰好相反.工作荷载作用导致产生负摩阻力的区域逐渐变小，位移零点上移至桩中部偏上位置.随着循环次数的增加，桩侧摩阻力逐渐增大，增量逐渐减小.
［1］李高建，胡玉叶.地源热泵技术的研究与应用现状［J］.节能技术，2007，25（2）：176-178.
［2］BRANDL H.Energy foundations and other thermoactive ground structures［J］.Géotechnique，2006，56（2）：81-122.
［3］MONIQUE M，BACH P M，BOUAZZA A，et al.Technological advances and applications of geothermal energy pile foundations and their feasibility in Australia［J］.Renewable and Sustainable Energy Reviews，2010，14（9）：2683-2696.
［4］KNELLWOLFC，PÉRON H，LALOUIL，etal.Geotechnical analysis of heat exchanger piles［J］.Journal of Geotechnical and Geoenviron Engineering，2011，137（10）：890-902.
［5］YAVARIN，TANGAM，PEREIRAJM，etal.Mechanical behaviour of a small-scale energy pile in saturated clay［J］.Géotechnique，2016，66（11）：878-887.
［6］桂树强，程晓辉.能源桩换热过程中结构响应原位试验研究［J］.岩土工程学报，2014，36（6）：1087-1094.
［7］WU D，LIU H L，KONG G Q，et al.Displacement response ofan energy pile in saturated clay［J］.Proceedings of the Institution of Civil Engineers-Geotechnical Engineering，2018，171（4）：285-294.
［8］方鹏飞，ABDELMALEK B，WANG B，等.温度对地热能源桩承载性状的影响［J］.防灾减灾工程学报，2017，37（4）：540-545.
［9］STEWARTM A，MCCARTNEY JS.Centrifuge modeling of soil-structure interaction in energy foundations ［J］. Journal of Geotechnical and Geoenvironmental Engineering， 2014， 140 （4） ：04013044.
［10］MORINO K，OKA T.Study on heat exchanged in soil by circulating water in a steel pile［J］.Energy and Buildings，1994，21（1）：65-78.
［11］GASHTIEH N，MALASKAM，KUJALAK.Evaluation of thermo-mechanical behaviour of composite energy piles during heating/cooling operations［J］.Engineering Structures，2014，75：363-373.
［12］郝耀虎，孔纲强，彭怀风，等.桩端约束对桩身热力学特性影响的模拟分析［J］.防灾减灾工程学报，2017，37（4）：532-539.
［13］王成龙，刘汉龙，孔纲强，等.工作荷载下温度循环对桩基变形与应力的影响分析［J］.岩土力学，2016，37（增1）：317-322.
［14］刘干斌，谢琦峰，范高飞，等.饱和黏土中热交换桩承载力特性模型试验研究［J］.岩石力学与工程学报，2017，36（10）：2535-2543.
［15］费康，钱健，洪伟，等.黏土地基中能量桩力学特性数值分析［J］.岩土力学，2018，39（7）：2651-2661.
［16］王中豪.地源热泵桩基埋管的三维数值模拟与应用分析［D］.武汉：华中科技大学，2016.
［17］翁凯文.软土地基能源桩桩-土界面摩擦性状研究［D］.杭州：浙江工业大学，2019.
［18］毛笑笑，周开发.COMSOL在土木工程多物理场问题中的应用前景［J］.智能城市，2017，3（1）：28-29+51.
［19］孔为，肖建昆，陈代芬.Comsol Multiphysics在“节能技术”教学中的应用［J］.科技资讯，2013（6）：209.
［20］BABUS’HAQ R F.Forced-convection heat transfer from a pipe to air flowing turbulently inside it［J］.Experimental Heat Transfer，1992，5（3）：161-173.