地下换热管安装施工过程中,由于受各种不确定因素的影响,地下换热管很容易出现弯曲,这种初始弯曲状态因施工条件等实际情况影响而具有随机性。为了便于计算分析,假设两种具有典型初始弯曲态的管体几何模型,分别为少弯的三弯曲和多弯的七弯曲管体模型Ⅰ、Ⅱ,最大弯曲挠度均为 P=0.1 m,模型Ⅲ为管体无弯曲理想模型,深度方向总长为50 m,如图1所示。

建立的模型基于以下假设:①换热管的初始弯曲为平面弯曲;②认为管体轴向温度均匀分布;③材料物性各向同性,回填土与围土的物理性质相近;④忽略初始应力、渗流及孔隙水对应力、应变和温度的影响。

图1

Fig.1

	Figure Option ViewDownloadNew Window
	图1 换热管模型Fig.1 Heat exchange pipe model

管体及回填土热形变弹性部分计算基于热弹性力学模型,根据线性热应力理论,物体内部微元体总应变由应力和温度变化共同作用产生,从而把基于广义虎克定律的弹性模型推广到包括热应力与应变在内^{[ 10]}的弹性模型;回填土及周围土壤热形变中的塑性形变部分采用摩尔库伦(Mohr-Coulomb)模型^{[ 11]}计算。结构热形变过程中,换热管与回填土之间存在高度非线性的复杂接触行为,需要添加接触单元^{[ 12]},接触单元力学模型采用各向同性的扩展库伦摩擦模型,建立接触面上的剪切应力与法向应力函数关系^{[ 13]}。

由于几何对称性,取模型一半结构用于计算,计算域半径为1 m,如图2所示。管体外径为0.032 m,壁厚为3.5 mm,钻孔直径为0.2 m,管中心间距为0.08 m。温度场计算中模型初始温度为10℃,管体内壁温度随时间变化,分为升温和降温两个阶段,如图3所示。形变位移计算中,换热管与围土间采用罚函数摩擦算法,摩擦因数为0.1,法向接触形式采用“硬接触”^{[ 14]};模型的对称面施加 X方向对称约束条件(XSYMM),顶、底面为完全固定约束,其他外侧面施加 X、 Y方向固定位移约束,以防止模型在任意方向上的刚体运动,确保计算不会产生数值奇异。采用六面体三维网格单元离散化几何模型,分别设置网格单元属性DC3D8及C3D8R用于温度场及位移场计算。Ⅰ、Ⅱ、Ⅲ三种结构模型的网格单元总个数分别为112 000、119 500、110 000,轴向分布500个网格节点,在靠近换热管附近区域加密网格,以增加单元数量来提高接触面附近区域参数的计算精度。

图2

Fig.2

	Figure Option ViewDownloadNew Window
	图2 计算域及边界Fig.2 Computational domain and boundary

图3

Fig.3

	Figure Option ViewDownloadNew Window
	图3 边界温度Fig.3 Temperature of boundary

首先用ABAQUS传热计算模块进行非稳态热传导数值计算,然后将温度场分布计算结果作为边界条件,输入到模型形变位移计算中进行热力耦合计算;有限元动态显式算法(ABAQUS/Explicit)可以模拟接触条件和其他不连续情况,且不需要迭代和收敛准则,能够有效计算结构形变过程中管体与回填土接触面间复杂的接触关系,所以采用有限元动态显式算法进行模型的热形变位移计算,该计算方法的有效性已经在文献[15]中得到了实验验证。

图5

Fig.5

	Figure Option ViewDownloadNew Window
	图5 受力状态Fig.5 Stress state

在升温段(储热工况),热膨胀形变、固定边界和回填土的约束作用使管体产生轴向、径向的热应力。由图5可知,管体的径向热应力由轴向热应力 σ_S的径向分力 σ_y与管体径向热应力 σ_R的合力构成,共同挤压回填土出现弹塑性应变,管间区域回填土同时受到进、出水管形变的挤压和剪切作用,导致此区域回填土 Y、 X向塑性应变呈负、正值,且绝对值大于其他区域,以特征面1塑性应变分布为例,如图6所示。

图6

Fig.6

	Figure Option ViewDownloadNew Window
	图6 回填土塑性应变区域Fig.6 Plastic strain area of backfill

图7

Fig.7

	Figure Option ViewDownloadNew Window
	图7 特征点相对位移Fig.7Relative displacement of characteristic points

由图7可知,随着边界温度的不断升高, A、 B特征点的相对位移持续增加, A_s与 A_p、 B_s与 B_p的相对位移值一致,管体和回填土在升温过程中保持接触状态。事实上,由于钻孔尺寸的限制,管体最大弯曲挠度较小,导致 σ_y相对于 σ_R为一个微小值,所以管体径向膨胀形变占主要地位。两种模型与理想模型(模型Ⅲ)相比较,其特征点相对位移值趋于一致。在第7200、14 400、21 600 s时, A特征点相对位移量分别是 B特征点相对位移量的2 .87、1 .45、1 .61倍,这是由于特征点 B位于管间区域,其 Y轴负方向的形变位移受到出水管体 Y轴正方向热形变位移的抑制作用,随着温度的升高,抑制作用整体呈减弱趋势。

在降温段(取热工况),特征点相对位移减小。降温段中期(28 800 s),位于管间区域的特征点 B相对位移值较升温段中期(14 400 s)减小,说明管间区域回填土热形变回弹能力减弱。降温段末期,不同模型管-土特征点的径向相对位移出现差异,其差值即为间隙宽度值,简称间隙值。

特征面1处于管体 Y正向弯曲最大挠度位置。由图7(a)可知,管体弯曲凸侧的特征点 A_p、 A_s相对位移几乎无差异,而弯曲凹侧特征点 B_p与 B_s的相对位移方向相反,表明 B_s未能跟随 B_p回弹而导致管与回填土接触面分离,在管体弯曲凹侧出现值为4.0×10^-4m的间隙。这是由于管体径向热应力大幅减小使其迅速回弹,而管体轴向尺度远大于径向,轴向残留较大的热应力,在这种情况下 σ_y相对于 σ_R不再是一个微小值。在管体自身径向回弹形变的同时, σ_y促使管体沿弯曲方向继续移动,含有塑性应变的回填土未能跟随管体回弹,从而导致管体弯曲凹侧出现较大间隙。特征面2位于管体弯曲挠度较小的位置,其特征点相对位移状态与理想模型接近,如图7(b)(e)所示,弯曲凹侧特征点间隙值为2.9×10^-4m,小于特征面1弯曲凹侧的特征点间隙值,相应的弯曲凸侧特征点仍然存在1.0×10^-4m的间隙,说明弯曲挠度越小的位置 σ_y越弱,即管体弯曲凸侧向管体弯曲方向形变位移的能力越差。

特征面1'、2'分别位于管体 Y轴正、负向弯曲最大挠度位置,由图7(c)(d)可知,随管体弯曲方向的不同, B_p、 A_p特征点均明显出现与 B_s、 A_s相反方向的相对位移,分别在特征面1'、2'中管体弯曲凹侧特征点出现5.8×10^-4m、6.1×10^-4m的间隙,均大于特征面1、2中的特征点间隙值,这说明同一钻孔中相同最大弯曲挠度的管体弯曲段越多,管体弯曲凹侧的间隙越大。

以特征面2'、3的间隙分布为例进行分析,限于篇幅其他模型特征面间隙分布未列出。由图8(a)可知,对于 Y正向弯曲的管体,沿特征点 A-1 -8 -B、 A'-1 '-8 '-B'间隙值呈现递增趋势,分布在进、出水管体弯曲凹侧的特征点4 -B、4 -B'间隙值大于平均值(3.6×10^-4m、3.9×10^-4m),在约占1/2管体周长范围内分布,出水管间隙值整体大于进水管间隙值,特征面1'间隙分布规律与之相反。对于特征面3,间隙值沿特征点 A-1 -8 -B、 B'-8 '-1 '-A'呈现递增趋势,进、出水管体大于平均间隙值(3.4×10^-4m、3.1×10^-4m)的特征点分布范围主要集中在管间区域,如图8(b)所示。

图8

Fig.8

	Figure Option ViewDownloadNew Window
	图8 特征面间隙分布(×10-4m)Fig.8 Gap distribution on characteristic face(×10-4m)