自由约束条件下能源桩的离散元研究

王雪松， 龚晓南

（浙江大学建筑工程学院滨海中心，杭州 310012）

0 引言

气候变化和能源危机是如今面临的两大巨大挑战，当前形势下新能源以其具有可再生、无污染、形式种类丰富等优点深受各国青睐。

最早由日本的Morino［1］提出桩埋管换热器，而能源桩是一种新型桩基础形式，该桩基础结合地埋管换热器的塑料换热管与建筑物的混凝土桩基以实现结构物承载及换热功能。在国外方面，1948年Kelvin［2，3］提出线热源理论，该理论假定竖直埋管换热器为一恒定热流密度的无限长线热源。Lee等［4］根据半离散化网络的研究方法，提出三维有限差分的数值模型。Yoon等［5］基于有限元模拟分析，对比研究了不同埋管形式的能量桩在3个月内通过间歇运行和连续运行下的换热性能。Bezyaii等［6］通过三维数值模拟研究了不同换热管的能量桩的换热性能，发现螺旋型埋管换热性能优于单U型和W型埋管。

在国内方面，石磊［7］推导出了热源模型的温度分布解析解，并进一步提出了线圈模型。桑宏伟［8］尝试使用二维离散元模型对能源桩的传热特性进行研究，采用MatDEM离散元软件对能源桩建立模型并通过与现场工程监测结果对比，验证了离散元法在能源桩换热性能研究方面的有效性。

文中提出能源桩的离散元模型，并通过MatDEM软件对能源桩的传热性能和承载能力特性进行研究。

1 模型基本理论

离散元法由Cundall［9］首先提出，并将其主要用于计算分析非连续介质的各种力学行为。离散元受限于实际工程颗粒数量需求巨大导致的计算量较大等问题，在岩土工程实际应用中相对较少。而MatDEM软件采用全新的算法较好地解决了计算量问题，目前已经可实现应力场-位移场-温度场的多场耦合，给离散元的工程实际应用带来了可能。

1.1 离散元基本理论

在文中离散元模型中采用线弹性赫兹接触模型如下：

式中，F为接触力合力，Fn；Fs分别为法向接触力和切向接触力。

颗粒间的法向弹簧可以表示法向力Fn和法向变形Xn的关系，见式（2）。

式中，Kn为法向刚度；Xn是法向相对位移；Xb为断裂位移。

同理切向弹簧可以表示颗粒间的剪切力Fs和剪切变形Xs的关系，见式（3）。

式中，Ks为切向刚度；Xs为切向位移。

1.2 传热模型基本理论

1.2.1 传热学相关理论

文中中主要考虑颗粒间的热传导，热传导是介质内无宏观运动时的一种传热方式。根据傅里叶定律可推导出传热方程：

式中，T为温度场向量；ρ为质量密度；c为单位体积内的热量；t为时间；λ为导热系数。

1.2.2 离散元颗粒传热理论

根据有限差分法可以推导出离散元颗粒的传热方程［10］。

推导中做如下假设：

两相邻热源间的温度差与热阻间的关系。

式中，ΔT为温度改变；η为单位长度的热阻；L为热管长度(此处简化为两颗粒中心距离)；n→i为热管方向向量；q→i为两颗粒间的热通量。

在模型中，假设温度的变化只影响颗粒体积，且颗粒根据温度变化线性膨胀，如式（6）所示。

式中，ΔR为颗粒因温度改变的半径；α为热膨胀系数，m/℃；R为颗粒初始半径；ΔT为温度改变量。

在离散元体系中将颗粒视为独立的热源和邻近颗粒的热管构成，热管作为热量的传导媒介，以颗粒是否接触来判断热量是否传递，若颗粒间无接触则不传到热量，若两颗粒重叠则视为颗粒间的热管进行传热。如图1所示，T0颗粒与T1、T2、T3颗粒分别接触，颗粒形心之间形成3个热管，热管长度为形心间的距离Li=|R0Ri|，热流流经的截面Ai视为接触面积，ni为接触面指向低温物体的法向方向，T0颗粒每次通过3个热管与周围颗粒实现热传导，如图1所示。

图1 颗粒间传热

根据以上的假设，推导出颗粒的传热方程：

式中，Q（m）为颗粒周围热源流入该颗粒的热量；Qv为内热源热量；m为质量；cv为比热容。

2 能源桩几何模型建立

文中根据以上离散元传热模型，并根据郭易木等［11］论文中PHC能源桩原型试验建立能源桩几何模型并验证模型有效性。

2.1 现场试验概况

现场勘查显示地下水位于1.0m深左右，于0～3.6m为硬壳层，土层分界线位于地下11.1m左右，分界线以上土层薄弱，其下土层强度较高，土层的热导率在1.15～1.17W/（m·K）范围内。现场试验中PHC能源桩桩长为24m，桩外径为500mm，内径为280mm，桩体材料为C80高强混凝土，换热管为单U型埋管。桩体弹性模量4.57×104MPa，热膨胀系数αc=11.57×10-6/℃。

2.2 离散元模型

2.2.1 土体模型

文中建立能源桩的三维离散元模型，模型盒长、宽均为4m，高度为28m，为保证计算的精确度并降低模型计算量，文中采用从中心到模型盒边缘粒径扩大的方式建立土体模型，颗粒半径分布如表1所示。

表1 颗粒半径分布 m

模型总颗粒数为571940个，其中活动颗粒为499739个，边界单元72201个，在模型中将边界单元定义为隔温光滑边界，不计算热传导。

为通过与现场试验对比验证模型有效性，参照现场土体条件设置模型的力学参数，并经过多次试算确定模型的微观力学参数见表2。

表2 颗粒力学性质

2.2.2 能源桩桩体模型

能源桩直径D1=500mm，桩长L=24m，为保证桩体光滑，桩体采用颗粒半径为0.03m的颗粒以一定系数重叠堆积而成，并设置为clump团簇，能源桩的U型传热管在离散元中以颗粒竖向堆叠为两根直管并紧贴桩体设置为clump团簇，以此进行传热模拟，能源桩模型如图2所示。

图2 能源桩模型

2.3 模型建模步骤

（1） 初始模型堆积，先建立1.1m×1.1m×28m的模型盒B1，随机生成颗粒半径为0.03m的颗粒经过重力沉积和压缩。再建立2.1m×2.1m×28m的模型盒B2，随机生成颗粒半径为0.04m的颗粒经过重力沉积和压缩，再切割B1模型盒生成半径0.5m的圆柱体结构体插入B2模型盒并平衡模型。同理依次生成不同颗粒半径、不同尺寸的模型盒，并不断切割结构体插入新的模型盒。

（2） 根据现场试验条件对地层进行简单划分，以便模型建立完成后赋予不同地层不同颗粒力学性质，地层划分如图3（a）所示。

图3 能源桩建模过程

（3） 建立能源桩及U型换热管并插入堆积好的土体模型中，给桩体及土体赋予不用的力学性质，先断开土体胶结，让土体在重力作用下填充间隙，模型充分平衡后再胶结地层与能源桩并最终实现模型平衡，如图3（b）所示。

2.4 传热过程模拟

离散元的传热过程是通过给不同颗粒设置不同的温度，以形成颗粒间的温度差来进行热传导模拟。能源桩的传热过程通过以下几步进行模拟：

（1） 先赋予颗粒初始温度，再给U型换热管设置相应的加热温度。

（2） 在每个时间步内根据有限差分法思想计算每个颗粒不同时间的温度值。

（3） 根据颗粒温度的改变，通过考虑颗粒线性膨胀更新颗粒半径。

（4） 根据每一个时间步的颗粒半径的改变，对模型进行堆积和平衡。

3 模拟结果和分析

3.1 桩体温度与分布规律

模型试验为对比现场试验，通过设置U型换热管与桩土颗粒的温度差进行168h的传热模拟，将U型换热管颗粒设为恒温40℃，依据现场试验中土体及桩体温度对模型中桩土颗粒初始温度进行设置，温度设置见表3。

表3 土体及桩体颗粒初始温度

图4给出了数值模拟试验和郭易木等现场试验0、24、72、168h的桩身受热后的温度对比图。根据模拟结果可得，能源桩加热初始阶段能源桩内土体温度及能源桩桩体温度上升较快，但桩内土体温度分布具有一定差异，与桑宏伟论文中观测现象一致，其原因一方面是靠近U型换热管颗粒的土体颗粒温度上升较快，另一方面是传热桩内的土体颗粒分布不均匀，无法像连续介质一样完全均匀传热。当桩体温度上升超过30℃后传热趋于平稳，桩身温度增长减缓，桩内土体温度分布趋于稳定，且分布较传热初始阶段更为均匀。

图4 桩身温度分布

168h时桩身平均温度增长18.8℃，桩体最高温度达到38℃，桩身温度总体呈现两端低中部高的形式，原因分析为一方面由于桩体初始温度分布呈现两端低中间高现象；另一方面由于桩两端散热条件良好。但随着传热，顶部桩体温度升温较快，原因分析为传热模型未考虑顶部土体及桩体向上方空气中散热的影响，在后续的建模中可对此进行修正，模拟结果的桩体温度与分布规律与现场试验结果基本一致。

桩外土体每隔5m依次取靠近桩体0.1m的土体平均温度，不同时刻温度分布如图5所示，可见土体颗粒温度在0～48h桩外土体温度上升较小，说明初始阶段能源桩传热对桩外土体影响较小。

图5 桩周土体温度分布

桩外土体在48～168h持续稳定上升，在168h时桩外土体最高温度达到26.6℃，分布规律整体呈现中间高两端低，是由于桩端以下土体散热空间更大，没有周围热量互相影响，因此桩底土体散热较快，而模型顶部因该离散元模型未考虑土体及桩体与空气交换热量而导致顶部土体温度上升较快。

3.2 桩周土体位移

能源桩模型中通过U型换热管传热，使桩、土体温度发生改变并根据式（7）产生相应的颗粒热膨胀，取桩体及土体温度变化处即埋深16.3m水平面的土体颗粒分析桩周土体位移变化规律。

如图6、图7所示，桩体内部回填土体因受桩体限制颗粒位移较小，而随着热量传递可见桩外土体产生相应的位移。桩外土体沿X、Y方向均产生对称、相反的位移分布，越靠近桩体的土体相对位移越大。而模型中U型管紧贴桩体分布，靠近U型管的桩体、土体温度上升较快，邻近土体相对其他位置土体沿着桩体径向产生更大的位移，即桩外土体沿X方向位移相对沿Y方向位移绝对值较大。

图6 桩周土体沿X方向位移

图7 桩周土体沿Y方向位移

3.3 桩身纵向应力变化规律

桩身受热产生由中间向两端的线性热膨胀，而桩周土体及桩底土体对桩身产生约束，能源桩以钻孔放入土体颗粒中进行平衡的形式，故可视为桩体除自身重力作用外无初始应力，能源桩的纵向附加热应力如图8所示。模拟结果为与郭易木等现场试验对照，取现场试验应变片位置的颗粒进行附加应力计算，可见能源桩随着加热产生轴向压应力，传热结束后桩体应力总体呈现两端小、中间大的趋势，但相对现场试验测得结果较小，原因分析为与现场试验相比模型颗粒孔隙度较大、土体给桩体约束较小。在几个分析点中桩体的最大附加应力达4.12MPa，位于16.3m处，最小附加应力为桩顶部，与现场试验现象基本一致。

图8 桩身应力分布

4 结语

文中使用南京大学自主研发的MATDEM软件建立能源桩三维离散元传热模型模拟PHC能源桩（U型换热管）传热过程，并与郭易木等论文中PHC能源桩原型试验的模拟结果进行验证对比，得出以下结论：

（1） 在保证计算精度的前提下，可利用粒径优化方法，从模型中心向边缘逐层堆积不同粒径离散元颗粒来降低离散元模拟的计算量，文中验证了该方法在能源桩的传热模型中的可行性。

（2） 模型中桩体温度上升趋势与郭易木等现场试验结果基本一致，桩体温度呈现两端低中间高的变化，但由于离散元模型中无法考虑土体及桩体向空气中散热，因此模拟试验中桩体温度相对现场试验结果较高，在后续的离散元模型中可以考虑对模型改进以使结果更准确。

（3） 在能源桩传热过程中，能源桩与桩周土体随着传热的进行温度发生改变，且桩周土体随着能源桩桩体膨胀而沿着远离能源桩方向产生对称、相反的位移分布。

（4） 将离散元传热模型应用到实际工程的研究中具有可行性，文中已验证U型换热管能源桩的传热模拟与现场试验结果具有较高的重合度，并且根据文中的建模思想调整换热管颗粒模型可以实现不同形状换热管的能源桩的传热模拟，如单U型、双U型、W型、螺旋型换热管的模拟均可实现。