一种能源桩热力耦合作用下轴力计算方法

摘要

能源桩是一种高效新型的地源热泵系统，结构设计时需要考虑温度荷载引发的附加温度应力，为了解决这个问题，本发明提出了一种能源桩轴向内力计算方法。基于能源桩温度荷载变形特点—以中性点为界发生膨胀或收缩，以及温度应力大小由约束度决定的特点，结合荷载传递法的思想，构造了桩在热力耦合作用下一维有限差分计算方法。本发明可以计算桩在不同温度荷载用下的附加温度应力，以及热力耦合作用下的桩身内力大小及分布；考虑了温度对桩周土和桩端土对力学强度的影响以及不同土层参数的变化，并考虑了上部结构类型对温度应力的影响。

说明书

技术领域

本发明涉及浅层地温能利用新技术—能源地下结构技术领域，具体涉及一种能源桩热力耦合作用下轴向内力计算方法。

背景技术

能源地下结构是利用与地层接触的地下永久或临时围护结构来提取地温能的新技术，该技术是将换热管埋设在地下结构中，无需额外钻孔，相对传统钻孔埋管系统具有成本低、换热性能好、不额外占用地下空间等优点。能源桩是长细的垂直地下结构，能充分利用地下恒温层(浅层地温能利用价值最高)，换热性能更好，所以在实际工程应用最为普遍。能源桩自1984年首次在欧洲试用成功以后，安装数量出现指数增长，但是温度荷载可能带来的附加应力值得关注。Laloui团队(1999)首次对一根25m长能源桩进行原位试验，发现在+13.4℃温度荷载作用下，桩端产生了高达2MPa的附加温度应力，改变了传统桩身内力分布，后期有限的原位试验都证明温度荷载对桩身应力的大小和分布的影响不可忽视，所以能源桩的设计必须要考虑温度应力。

但是，由于当前国际上还缺乏能源桩设计规范和合理的内力计算方法，设计时主要依靠经验，采用放大安全系数的办法，比如瑞士和英国分别采用3和3.5的安全系数。显然这种粗糙的设计的方法会带来过大的桩基尺寸，并增加造价，影响这项绿色建筑节能技术的推广应用。

能源桩机理研究表明温度应力的大小取决于桩周围的约束度，比如上部结构对桩顶的约束，桩周土以及桩端持力层对桩的约束；桩以中性点为界发生变形当桩身温度升高，桩从中性点向两端膨胀，当桩身温度降低，桩的两侧向中性点收缩。所谓中性点，就是在温度荷载作用下，在沿桩身方向上不发生桩土相对位移的位置。Seed和Reese(1957)提出的荷载传递法被成功应用于桩基沉降和荷载计算，荷载传递位移协调法可以考虑不同土层力学参数的变化，在长期的实践和应用中，很多学者从不同角度，提出了多种荷载传递函数，但是荷载传递位移协调法要求在荷载作用下桩土相对位移的方向需要保持一致。根据能源桩的变形特点，可以以中性点为分界点，将桩分成两上下两部分，分别利用荷载传法的思想来计算能源桩内力和位移，还可以考虑温度对侧阻和端阻的影响。

发明内容

发明目的：能源桩设计需要考虑温度荷载引发的附加温度应力，当前能源桩设计依靠经验，采用简单放大安全系数的设计方法，造成桩基尺寸过大和成本增加等问题。为了合理设计能源桩，提出了一种能源桩热力耦合内力计算方法。该能源桩内力计算方法包括三个主要步骤：第一步、计算桩顶荷载P作用下桩身内力及桩土相对位移；第二步、确定在温度荷载△T作用下桩身上中性点的位置；第三步，以中性点为分界点将桩划分成上下两部分，然后以第一步的计算的结果为初始条件(力荷载作用下的位移)采用荷载传递法分别计算上下两部分在热力耦合(P+ΔT)作用下桩身的轴向内力。该方法结合了能源桩的结构变形特点和诱发温度应力的关键因素，能够合理计算能源桩热力耦合下的温度应力，为能源桩设计提供依据。

技术方案：

一种能源桩热力耦合作用下轴力计算方法，包括步骤：

(1)桩周围约束的定义：定义上部结构与桩顶的约束刚度K_h，桩侧土层对桩单元i的约束定义为一个非线性的弹簧作用在单元中部，持力层对桩端的约束定义为一个非线性的弹簧作用在桩端；获得桩侧与桩端在常温和变温条件下的荷载传递函数；

(2)根据步骤(1)计算在桩顶外力载荷P作用下的桩身轴力以及桩土相对位移；

(3)假定中性点单元，以中性点为界将桩划分成上下两部分；所述中性点为桩在温度荷载ΔT作用下，在轴向不发生桩土相对附加位移的位置；

(4)根据步骤(2)得到的在桩顶外力载荷P作用下的桩身轴力及桩土相对位移以及步骤(3)确定的中性点，采用荷载传递法分别计算中性点上下两部分在热力耦合作用下桩身的轴力。

所述步骤(2)中计算在桩顶外力载荷P作用下的桩身轴力以及桩土相对位移具体如下：

符号说明：Q、F、Z分别代表桩身荷载、侧摩阻力和位移；下标H、S和B分别代表单元顶部、中部和底部位置；下标M、T和MT分别代表力学荷载、温度荷载单独作用和热力耦合作用；

(21)根据桩的尺寸和精度将桩离散成N个线弹性单元，每个单元i的刚度K_i＝E·A_i/L_i，其中，E为桩的弹性模量，A_i和L_i分别表示单元i的横截面积和长度；根据胡克弹性定律，单元i在外力荷载作用下的压缩量Δ_i：

(22)则单元i顶部位移Z_H,M,i和中部位移Z_S,M,i：

Z_H,M,i＝Z_B,M,i+Δ_i

(23)根据单元i的中部位移和荷载传递函数，计算单元i的侧阻F_S,M,i：

F_S,M,i＝f(Z_S,M,i)

(24)根据单元i静力平衡，计算单元i新的顶部荷载Q_H,M,i,new，：

Q_H,M,i,new＝Q_B,M,i+Q_S,M,i＝Q_B,M,i+f(Z_S,M,i)

(25)验证|Q_H,M,i,new-Q_H,M,i|＜10^-6是否成立；若成立，则满足精度要求，计算下一个单元；否则，假设单元i的顶部荷载Q_H,M,i＝Q_H,M,i,new，重复步骤(21)～(24)，直到满足|Q_H,M,i,new-Q_H,M,i|＜10^-6；

(26)根据计算得到的每个单元的轴力和位移计算得到整个桩身的轴力以及桩土相对位移。

步骤(3)中确定中性点的位置具体为：

(31)假设在单元NP的顶端为中性点位置，那么Z_T,H,NP＝0或者Z_T,B,NP-1＝0；上部分桩单元序号为1～NP-1，下部分桩单元NP～N；

(32)对于桩下部单元：

1)假设单元i的温度变形量Δ_i为自由变形，即

Δ_i＝L_iαΔT

其中α为桩的线膨胀系数。

2)计算单元i位移：

顶部位移：Z_T,H,i＝Z_T,B,i-1；

中部位移：其中+表示ΔT＞0条件下，-表示ΔT＜0；

底部位移：Z_T,B,i＝Z_T,H,i±Δ_i，其中+表示ΔT＞0条件下，-表示ΔT＜0；

其中单元NP，

顶部位移：Z_T,H,NP＝0；

中部位移：其中+表示ΔT＞0条件下，-表示ΔT＜0；

底部位移：Z_T,B,i＝±Δ_NP，其中+表示ΔT＞0条件下，-表示ΔT＜0；

3)计算单元i轴力：

底部轴力：Q_T,B,i＝Q_T,H,i+1；

侧阻：F_T,S,i＝f_T,S,i(Z_T,S,i)；

顶部轴力：Q_T,H,i＝Q_T,B,i+F_T,S,i；

其中，单元N的底部轴力：Q_T,B,N＝f_T,B,N(Z_T,B,N)；

4)计算单元i的变形

5)验证|Δ_i,new-Δ_i|＜10^-6是否成立：

如果成立计算结束；如果不成立，则将Δ_i＝Δ_i,new，重复1)～4)迭代，直到满足|Δ_i,new-Δ_i|＜10^-6；

(33)对于桩上部单元：

1)假设单元i的温度变形量Δ_i为自由变形，即

Δ_i＝L_iαΔT

2)计算单元i位移：

底部位移：Z_T,B,i＝Z_T,H,i+1；

中部位移：其中+表示ΔT＞0条件下，-表示ΔT＜0；

顶部位移：Z_T,H,i＝Z_T,B,i±Δ_i，其中+表示ΔT＞0条件下，-表示ΔT＜0；

其中单元NP-1，

底部位移：Z_T,B,NP-1＝0；

中部位移：其中+表示ΔT＞0条件下，-表示ΔT＜0；

顶部位移：Z_T,H,NP-1＝±Δ_NP-1，其中+表示ΔT＞0条件下，-表示ΔT＜0；

3)计算单元i轴力：

顶部轴力：Q_T,H,i＝Q_T,B,i+1；

侧阻：F_T,S,i＝f_T,S,i(Z_T,S,i)；

底部轴力：Q_T,B,i＝Q_T,H,i+F_T,S,i；

其中，单元1的Q_T,H,1＝K_h·Z_T,H,1；

4)计算单元i的变形

5)验证|Δ_i,new-Δ_i|＜10^-6是否成立：

如果成立计算结束；如果不成立，则将Δ_i＝Δ_i,new，重复1)～4)迭代，直到满足|Δ_i,new-Δ_i|＜10^-6；

(34)能源桩的上下两部分保持静力平衡，即：

在计算时令NP＝1～N,计算出f_T(1)，…，f_T(N)，其中最小的f_T(NP)中NP即为中性点所在的单元，单元顶部的位置为中性点的位置。

所述步骤(4)中采用荷载传递法分别计算中性点上下两部分在热力耦合作用下桩身的轴向内力具体如下：

(41)对于桩下部单元：

1)假设单元i的温度变形量Δ_i为自由变形，即

Δ_i＝L_iαΔT

2)计算单元i位移：

顶部位移：Z_MT,H,i＝Z_MT,B,i-1；

中部位移：其中+表示ΔT＞0条件下，-表示ΔT＜0；

底部位移：Z_MT,B,i＝Z_MT,H,i±Δ_i，其中+表示ΔT＞0条件下，-表示ΔT＜0；

对于单元NP：

顶部位移：Z_MT,H,NP＝Z_M,H,NP；

中部位移：其中+表示ΔT＞0条件下，-表示ΔT＜0；

底部位移：Z_MT,B,NP＝Z_M,H,NP±Δ_NP，其中+表示ΔT＞0条件下，-表示ΔT＜0；

3)计算单元i轴力：

底部轴力：Q_MT,B,i＝Q_MT,H,i+1；

侧阻：F_MT,S,i＝f_T,S,i(Z_MT,S,i)；

顶部轴力：Q_MT,H,i＝Q_MT,B,i+F_MT,S,i；

4)计算单元i的变形：

(42)对于桩上部单元：

1)假设单元i的温度变形量Δ_i为自由变形，即

Δ_i＝L_iαΔT

2)计算单元i位移：

底部位移：Z_MT,B,i＝Z_MT,H,i+1；

中部位移：其中+表示ΔT＞0条件下，-表示ΔT＜0；

顶部位移：Z_MT,H,i＝Z_MT,B,i±Δ_i，其中+表示ΔT＞0条件下，-表示ΔT＜0；

对于单元NP-1：

底部位移：Z_MT,B,NP-1＝Z_M,B,NP-1；

中部位移：其中+表示ΔT＞0条件下，-表示ΔT＜0；顶部位移：Z_MT,H,NP-1＝Z_M,B,NP-1±Δ_NP-1，其中+表示ΔT＞0条件下，-表示ΔT＜0；

3)计算单元i轴力：

顶部轴力：Q_MT,H,i＝Q_MT,B,i+1；

侧阻：F_MT,S,i＝f_t,S,i(Z_MT,S,i)；

底部轴力：Q_MT,B,i＝Q_MT,H,i+F_MT,S,i；

4)计算单元i的变形：

验证|Δ_i,new-Δ_i|＜10^-6是否成立，如果成立计算结束；如果不成立，则将Δ_i＝Δ_i,new，重复中性点上部和下部单元的迭代，直到满足|Δ_i,new-Δ_i|＜10^-6，得到不同深度能源桩桩身轴力的大小。

有益效果：

本发明提供一种能源桩内力计算方法，可以计算桩在不同温度荷载用下的附加温度应力，以及热力耦合作用下的桩身内力大小及分布；考虑了温度对桩周土和桩端土对力学强度的影响以及不同土层参数的变化，并考虑了上部结构类型对温度应力的影响。该计算方法原理清楚，通过计算语言(如Matlab等)编制简单的迭代程序就可以求解，为能源桩的合理设计提供计算方法和理论依据。

附图说明

图1能源地下结构示意图。

图2为能源桩有限差分模型示意图。

图3为计算外力荷载P作用下桩身位移和内力的示意图。

图4为中性点确定示意图。

图5为热力耦合作用下桩身内力计算示意图。

图6为力荷载单元计算示意图。

图7为能源桩内力计算示意图。

具体实施方式

下面结合附图对本发明作更进一步的说明，下面的符合和上面“技术方案”中的符号一直。

图1为能源地下结构示意图。如图1所示，由于能源桩内力计算模型为一维有限差分模型，模型涉及到复杂热力耦合荷载和复杂桩周土和上部结构的非线性边界条件，需要采用位移协调的方法编写计算机程序迭代求解。

(1)桩周围约束的定义

能源桩附加温度应力的产生原因是由于在发生温度变化的时候不能自由膨胀或收缩。本发明借鉴传统荷载传递法的思想，结合能源桩结构和变形特点，将桩离散成N个线弹性单元，从桩顶到桩端的单元序号依次为1，…，N。为了单位统一，单元i的刚度可以表示为K_i＝E·A_i/L_i(kN/m)(其中，E为桩的弹性模量，A_i和L_i分别表示单元i的横截面积和长度)。

符号说明：Q、F、Z分别代表桩身荷载、侧摩阻力和位移；下标H、S和B分别代表单元顶部、中部和底部位置；下标M、T和MT分别代表力学荷载、温度荷载单独作用和热力耦合作用；i表示单元的编号。比如：

Q_M,H,i表示在桩顶荷载单独作用下单元i顶部的轴力；Q_T,S,i表示在温度荷载单独作用下单元i中部的轴力；Q_MT,B,i表示在热力耦合作用下单元i底部的轴力。

F_M,S,i表示在桩顶荷载单独作用下单元i的侧阻；F_T,S,i表示在温度荷载单独作用下单元i的侧阻；F_MT,S,i表示在热力耦合作用下单元i的侧阻。

Z_M,H,i表示在桩顶荷载单独作用下单元i顶部的位移；Z_T,S,i表示在温度荷载单独作用下单元i中部的位移；Z_MT,B,i表示在热力耦合作用下单元i底部的总位移。

单元i在外力荷载作用下的压缩量Δ_i(m)可以通过下面的公式计算：

上部结构对桩顶的约束被考虑，桩在温度荷载作用下，桩顶发生位移会产生由于上部结构的约束会产生附加桩顶荷载。由于上部结构主要为钢筋混凝土材料，所以将上部结构对桩的约束定义为一个刚度为K_h(kN/m)的线性弹簧约束，桩顶的位移(上升或下降)将改变桩顶的荷载，刚度系数的大小与上部结构类型有关；温度荷载产生的附加桩顶荷载可以表示为：

Q_T,H,1＝K_h·Z_T,H,1

桩侧土层对桩单元i的约束定义为一个非线性的弹簧作用在单元中部，同样，持力层对桩端的约束定义为一个非线性的弹簧作用在桩端，即单元N的底部。桩周土对桩侧(侧阻)以及持力层对桩端(端阻)的约束，可以采用荷载传递函数来定义，桩土荷载传递函数可以通过实测来确定。自从Seed和Reese(1957)提出荷载传递法的位移协调法以来，多种荷载传递函数被提出用于描述不同桩型和不同桩周围的地质条件的桩土作用。由于每个单元都可以对应一个荷载传递函数，所以可以考虑不同土层参数的变化，甚至能够考虑温度变化对桩土作用的影响。桩侧与桩端在常温和变温条件下的荷载传递函数：

在常温下，桩侧阻和端阻分别可以表示为：

F_M,S,i＝f_M,S,i(Z_M,S,i)

Q_M,B,N＝f_M,B,N(Z_M,B,N)

在温度T作用下，桩侧阻和端阻可以表示为：

F_T,S,i＝f_T,S,i(Z_T,S,i)

Q_T,B,N＝f_T,B,N(Z_T,B,N)

考虑温度T条件下桩卸载函数为：

F_T,S,i＝f_t,S,i(Z_T,S,i)

Q_T,B,N＝f_t,B,N(Z_T,B,N)

(2)能源桩内力计算过程：

①桩顶外力P作用下桩身位移和轴力

1、假设单元i的顶端轴力为Q_M,H,i

2、根据胡克弹性定律，计算单元i的压缩量Δ_i：

3、计算单元i顶部位移Z_H,M,i和中部位移Z_S,M,i：

Z_H,M,i＝Z_B,M,i+Δ_i

4、根据单元i的中部位移和荷载传递函数，计算单元i的侧阻F_S,M,i：

F_S,M,i＝f_S,M,i(Z_S,M,i)

5、根据单元i静力平衡，计算单元i新的顶部荷载Q_H,M,i,new，：

Q_H,M,i,new＝Q_B,M,i+Q_S,M,i＝Q_B,M,i+f(Z_S,M,i)

6、验证|Q_H,M,i,new-Q_H,M,i|＜10^-6是否成立

如果成立，那么计算满足精度要求，同理计算下一个单元(i+1)；否则，假设单元i的顶部荷载Q_H,M,i＝Q_H,M,i,new，重复步骤1～5，直到满足|Q_H,M,i,new-Q_H,M,i|＜10^-6。

在实际计算时，从桩端单元1开始，假设桩端位移Z_B,M,N＝0，得到桩端轴力Q_M,B,N＝f(Z_M,B,N)，然后按1-6的迭代方法一次求出单元N，N-1，…1的侧阻：F_S,M,N，…，F_S,M,1。

整个能源桩保持静力平衡，其控制方程如下：

如果平衡方程满足，则结算结束；如果平衡方程不等于0，则

Z_B,M,N＝Z_B,M,N+10^-6，重复上面的循环，直到满足整体平衡方程要求的精度。

②确定中性点的位置

所谓中性点，就是桩在温度荷载ΔT作用下，在轴向不发生桩土相对附加位移的位置。假设能源桩在取热和蓄热过程中，桩身温度分布均匀。假设在单元NP的顶端为中性点位置，那么Z_T,H,NP＝0或者Z_T,B,NP-1＝0。以中性点为界将桩分为两部分，上部分桩单元序号为1～NP-1，下部分桩单元NP～N。在迭代计算时，上部单元和下部分别利用荷载传递位移协调法计算，并满足精度要求。

对于桩下部单元(NP～N)：

1、假设单元i的温度变形量Δ_i为自由变形，即

Δ_i＝L_iαΔT

其中α为桩的线膨胀系数。

2、计算单元i位移

顶部位移：Z_T,H,i＝Z_T,B,i-1

中部位移：其中+表示ΔT＞0条件下，-表示ΔT＜0

底部位移：Z_T,B,i＝Z_T,H,i±Δ_i，其中+表示ΔT＞0条件下，-表示ΔT＜0

其中单元NP，

顶部位移：Z_T,H,NP＝0

中部位移：其中+表示ΔT＞0条件下，-表示ΔT＜0

底部位移：Z_T,B,NP＝±Δ_NP，其中+表示ΔT＞0条件下，-表示ΔT＜0。

3、计算单元i轴力

底部轴力：Q_T,B,i＝Q_T,H,i+1

侧阻：F_T,S,i＝f_T,S,i(Z_T,S,i)

顶部轴力：Q_T,H,i＝Q_T,B,i+F_T,S,i

其中，单元N

底部轴力：Q_T,B,N＝f_T,B,N(Z_T,B,N)

4、计算单元i的变形

5、验证|Δ_i,new-Δ_i|＜10^-6是否成立

如果成立计算结束；如果不成立，则将Δ_i＝Δ_i,new，重复1～4迭代，直到满足|Δ_i,new-Δ_i|＜10^-6。

对于桩上部单元(NP-1～1)：

1、假设单元i的温度变形量Δ_i为自由变形，即

Δ_i＝L_iαΔT

2、计算单元i位移

底部位移：Z_T,B,i＝Z_T,H,i+1

中部位移：其中+表示ΔT＞0条件下，-表示ΔT＜0

顶部位移：Z_T,H,i＝Z_T,B,i±Δ_i，其中+表示ΔT＞0条件下，-表示ΔT＜0

其中单元NP-1

底部位移：Z_T,B,NP-1＝0

中部位移：其中+表示ΔT＞0条件下，-表示ΔT＜0

顶部位移：Z_T,H,NP-1＝±Δ_NP-1，其中+表示ΔT＞0条件下，-表示ΔT＜0

3、计算单元i轴力

顶部轴力：Q_T,H,i＝Q_T,B,i+1，

侧阻：F_T,S,i＝f_T,S,i(Z_T,S,i)

底部轴力：Q_T,B,i＝Q_T,H,i+F_T,S,i

其中，单元1

Q_T,H,1＝K_h·Z_T,H,1

4、计算单元i的变形

5、验证|Δ_i,new-Δ_i|＜10^-6是否成立

如果成立计算结束；如果不成立，则将Δ_i＝Δ_i,new，重复1～4迭代，直到满足|Δ_i,new-Δ_i|＜10^-6

能源桩的上下两部分要求总体保持静力平衡，即：

在计算时令NP＝1～N，可以计算出f_T(1)，…，f_T(N)，其中最小的f_T(NP)中NP即为中性点所在的单元，单元顶部的位置为中性点的位置。

③计算热力耦合作用下能源桩内力

能源热力耦合计算时，在ΔT＞0的条件下桩的上部发生向上的膨胀，以及在ΔT＜0条件下桩下部发生向上的收缩，导致负摩阻力，在荷载传递函数上是一个卸载过程。在计算热力耦合荷载作用时，仍然要将中性点上部和下部分开进行迭代计算，详细过程如下：

对于桩下部单元(NP～N)：

1、假设单元i的温度变形量Δ_i为自由变形，即

Δ_i＝L_iαΔT

2、计算单元i位移

顶部位移：Z_MT,H,i＝Z_MT,B,i-1

中部位移：其中+表示ΔT＞0条件下，-表示ΔT＜0

底部位移：Z_MT,B,i＝Z_MT,H,i±Δ_i，其中+表示ΔT＞0条件下，-表示ΔT＜0

对于单元NP

顶部位移：Z_MT,H,NP＝Z_M,H,NP

中部位移：其中+表示ΔT＞0条件下，-表示ΔT＜0

底部位移：Z_MT,B,NP＝Z_M,H,NP±Δ_NP，其中+表示ΔT＞0条件下，-表示ΔT＜0

3、计算单元i轴力

底部轴力：Q_MT,B,i＝Q_MT,H,i+1，

侧阻：F_MT,S,i＝f_T,S,i(Z_MT,S,i)

顶部轴力：Q_MT,H,i＝Q_MT,B,i+F_MT,S,i

4、计算单元i的变形

对于桩上部单元(NP-1～1)：

1、假设单元i的温度变形量Δ_i为自由变形，即

Δ_i＝L_iαΔT

2、计算单元i位移

底部位移：Z_MT,B,i＝Z_MT,H,i+1

中部位移：其中+表示ΔT＞0条件下，-表示ΔT＜0

顶部位移：Z_MT,H,i＝Z_MT,B,i±Δ_i，其中+表示ΔT＞0条件下，-表示ΔT＜0

对于单元NP-1

底部位移：Z_MT,B,NP-1＝Z_M,B,NP-1

中部位移：其中+表示ΔT＞0条件下，-表示ΔT＜0

顶部位移：Z_MT,H,NP-1＝Z_M,B,NP-1±Δ_NP-1，其中+表示ΔT＞0条件下，-表示ΔT＜0

3、计算单元i轴力

顶部轴力：Q_MT,H,i＝Q_MT,B,i+1，

侧阻：F_MT,S,i＝f_t,S,i(Z_MT,S,i)

底部轴力：Q_MT,B,i＝Q_MT,H,i+F_MT,S,i

4、计算单元i的变形

验证|Δ_i,new-Δ_i|＜10^-6是否成立，如果成立计算结束；如果不成立，则将Δ_i＝Δ_i,new，重复中性点上部和下部单元的迭代，直到满足|Δ_i,new-Δ_i|＜10^-6，得到不同深度能源桩桩身轴力的大小。