地源热泵系统作为一种充分利用地热能资源为建筑供冷供热的高效系统，受到了广泛应用．世界范围内，地源热泵的安装容量从1995年的1 854 MW增长到2015年的50 258 MW，且这一趋势还在不断上升[1]．其中，基于土壤源的垂直地埋管换热器的系统因为占地面积小、换热面积大、运行效率高等特点，在工程上得到广泛的应用．在地埋管内一般常将单U或者双U型水管用于承载闭式循环中的热量输运，而大量实验及数值模拟研究[2-3]均已表明双U型地埋管的换热性能优于单U型地埋管，但是双U型地埋管换热器内部同时也存在着更加复杂的传热过程．地下换热器的传热模型是地源热泵系统设计、运行、优化、故障诊断等的重要基础，也是地源热泵研究的难点．基于通用求解器CFD的数值模型与基于热源法的解析模型分别是两大典型模型分支．其中，基于热源法衍生出的线热源[4]、圆柱面热源[5]、圆柱体热源[6]等模型方法为地源热泵的模拟计算提供了更高效与易用的工具．绝大多数基于热源法解析模型均基于二区域划分与建模思路[7]，即将地埋管内外分别看做两个建模子区域分别处理，在埋管外的土壤传热部分，可以将其看成线热源或者圆柱热源在半无穷大域内的热传导问题，而在埋管内部则是采用简化热阻传热的方式对水温进行模拟．更进一步看，不论单U或者双U型地埋管，其内部结构均不是完全中心对称的，因此无法像处理同轴套管型地埋管传热一样采用径向上单一热阻进行等效传热描述，这一问题触发了一系列的模型改进研究．代表模型包括将U型管结构进行等效处理的等效直径模型[8-10]、容积热阻模型(CaRM)[11]及热阻热容模型(TRCM)[12]，这些模型方法或改变地埋管内部几何特征，或用简化热网络替代内部复杂的非均匀非稳态传热机理．这种建模思路虽然可以避繁就简的实现模拟计算地埋管传热的目标，但是存在三个方面缺陷：a. 等效的直径或是热网络模型中的等效热阻与热容值难以准确确定，而且一旦回填材料热物性发生变化则很容易导致之前辨识的模型参数值不再适用；b. 一般在等效热阻热容模型中均忽略钻孔内部非均匀温度场规律，仅从等效的角度为水管和钻孔壁各自选取一个等效代表计算点，并采用热阻热容对各个计算点进行连接，在基本原理上该模型并不适合用于系统的参数分析研究，因为模型中的热阻热容值仅与某一特定条件一一对应，本身不具备广延性；c. 无法获取钻孔内真实的温度场、热流分布等情况，而采用CFD等数值模型去获取这一信息又存在计算资源需求过大的问题．本研究直接将被求解问题看成是U型管在孔内回填土与孔外土壤形成的二元复合介质中的热传导问题，建模中同时考虑地埋管在垂直方向上的热流不均匀特征，从解析建模角度提供一种计算精度高、容易操作、不依赖等效热阻热容假设的地埋管传热模型．1 模型方法利用线热源在二元复合介质中的非稳态传热解析模型函数[13-15]可以直接计算钻孔内外任一点温度的变化，但是该解析函数中存在着十分复杂的贝塞尔函数、无穷限定积分及无穷级数计算式，这为模型的应用设置了较大的阻碍，也使复合介质模型难以广泛应用．采用文献[16]中提出的模型假设：对于被计算点与热源所处同一地层的情况，直接采用该地层的物性参数进行相关传热计算；对于被计算点与热源不处于同一地层的情况，通过两者之间的连线穿越不同地层时的长度比例与不同地层的物性参数进行加权平均处理．具体如下ρs=(s1ρg+s2ρsoil)/(s1+s2),ks=(s1kg-1+s2ksoil-1)/(s1+s2)-1,Cs=(s1Cg+s2Csoil)/(s1+s2);as=ks/(ρsCs),式中：ρ，k，C和a分别为介质密度、导热系数、比热容及热扩散系数；下标s，g和soil分别表示复合介质、钻孔内回填材料及钻孔外土壤材料；s1与s2分别为计算点与热源连线在两者介质中行走的距离长度．从本质上说，二元复合介质模型也是一种二区域模型，只不过两个区域的分界面选在U型管的管壁而不是钻孔壁位置．在传统二区域方法中，区域划分为钻孔内的流体与回填土组成的部分以及钻孔外的土壤部分，其中由于钻孔内几何结构不是中心对称，难以获取完全解析解，于是只能采用基于热阻热容的简化模型进行替代．但是在复合介质思想指导下，计算区域被分为了管内流体与管外复合介质固体，因此无须采用等效处理的热阻热容进行建模，这不仅免去了热阻热容的参数辨识，而且还可以对钻孔内温度场进行模拟．在计算中，为了获取地埋管换热过程中水温的不断变化过程，须以U型管壁以及入口水温为边界条件，采用准动态传热方程对管内流体进行描述，即-mCpwdTfi1(z)dz=ki1Tfi1(z)-Twi1(z),mCpwdTfo1(z)dz=ko1Tfo1(z)-Two1(z); (1)ki1=(πdinhci1)-1+lndout/din/(2πλp)-1,ko1=(πdouthco1)-1+lndout/din/(2πλp)-1,式中：m为流体质量流量；Cpw为流体的比热容；T为温度；ki1和ko1分别为进水管中流体与管壁及出水管中流体与管壁的综合传热系数；hci1和hco1分别为水管内侧面和外侧面对流换热系数；din和dout分别为水管内直径与外直径；下标fi1和fo1分别为双U管中第一支U型管的进水管与出水管；下标wi1和wo1分别为双U管中第一支U型管的进水管管壁与出水管管壁．为了求解方便，对下水与上水管的壁温沿深度方向假设服从三次幂函数关系式，因此对方程(1)中的两个方程分别采用常微分方程求解方法获得一般解，然后利用入口水温边界条件，即可以计算出U型管内的水温函数Tf1(z)=f10+f11z+f12z2+f13z3+cx1e-a1z;Tf2(z)=f20+f21z+f22z2+f23z3+cx2e-a2z⋅式中：f10=c10-c11/a1+2c12/a12-6c13/a13；f11=c11-2c12a1+6c13a12；f12=c12-3c13/a1；f13=c13；f20=c20+c21a2+2c22a22-6c23a23；f21=c21+2c22a2+6c23a22；f22=c22+3c23/a2；   f23=c23；cx1=Tin-f10；cx2=[Tf1(H)-f20+f21H+f22H2+f23H3]/ea2H，其中，Tf1和Tf2分别为U型管中进水管与出水管的水温；z为沿程深度；H为U型管总深度；a1，a2，c10，c11，c12，c13，c20，c21，c22和c23均为公式拟合参数，无单位．在垂直地埋管换热计算中，埋管深度一般为50~150 m，所处区域内周围土壤环境温度比较均匀，因此大多数模型均采用钻孔壁热流沿深度方向均匀分布的假设．为了进一步提升模型的准确性，充分考虑了U型管内水温沿程变化导致管壁热流沿深度方向不均匀的问题，引入有限长分段圆柱热源SFCS模型[17]对该问题进行处理．SFCS模型的主要建模思路是将圆柱热源进行分段划分，每一段具有不同的热流密度，基于热源格林函数，对不同分段热源在空间与时间上积分，得到某一个分段热源强度作用下的温度场公式θir,z,t=ql(i)8πksπas∫0t1(t-t')3/2⋅I0rrb2ast-t'exp-r2+rb24ast-t'⋅∫(i-1)z0iz0exp-z-z'24ast-t'-exp-z+z'24ast-t'dz'dt',式中：θi为分段热源中第i段热源在土壤中任意位置(r，z)任意时间点t的过余温度；ql(i)为第i段热源的强度；rb为地埋管井半径；z0为分段热源的长度；z'和t'均为积分变量；I0(⋅)为零阶修正贝塞尔函数．然后利用该热源问题线性可叠加原理及动态热源强度可叠加原理，推导得θr,z,t=∑j=1Ns∑i=1nql(i,t/Δt-j)-ql(i,t/Δt-j-1)8πks∫0jΔt1jΔt-t'⋅I0rrb2asjΔt-t'exp-r2+rb24asjΔt-t'⋅Eiz0-z4asjΔt-t'-E(i-1)z0-z4asjΔt-t'-Eiz0+z4asjΔt-t'+E(i-1)z0+z4asjΔt-t'dt',式中：Ns为叠加时间段数；n为分段热源个数；Δt为时间步长；E(⋅)为误差函数．最后将二元复合材料处理方法、U型管流体传热解析函数、非均匀管壁热流处理方法及模型进行耦合，建立一种描述垂直双U型地下换热器的解析模型．该模型能够实现高效计算水流动态温度及钻孔内外温度场在时间尺度上的演变全过程．2 模型验证为了对新建立的系统模型准确性进行验证，选取文献[18]中提供的垂直双U型地下换热器热响应测试(TRT)数据作为支撑．采用Al-Khoury与Focaccia研究中提供的TRT数据以及地下换热器尺寸、物性参数等数据．地埋管深度为100 m，两管中心距为0.08 m，钻孔孔径为0.127 m，水管外壁直径为0.032 m，水管内壁直径为0.024 m，水管导热系数为0.42 W/(m⋅K)，流体导热系数为0.56 W/(m⋅K)，流体动力黏度为0.001 Pa⋅s，流体流速为0.42 m/s，流体比热容为4 180 J/(kg⋅K)，回填土密度为1 420 kg/m3，回填土导热系数为0.6 W/(m⋅K)，回填土比热容为1 197 J/(kg⋅K)，大地土壤导热系数为2.15 W/(m⋅K)，大地土壤比热容为400 J/(kg⋅K)，大地土壤密度为1 680 kg/m3，大地土壤初始温度为13℃．设置实验系统参数，在给定的条件下进行模型计算，短期3 h和中期4 h输出出口模拟水温与实验值对比如图1和2所示，图中θ为水温．采用数据误差分析指标进行计算，发现本模型的计算准确度在短期3 h和中期40 h的TRT数据中的平均相对误差分别仅为2.25%与0.17%；若采取均方根误差表示，则这两部分的均方根误差分别为0.51和0.21℃．此外，通过将本模型与基于等效热阻热容模型[18]对比，发现两者在同样的边界条件下，短期3 h和中期40 h的TRT数据误差分别是0.39和0.18℃，而本模型的优势在于保持高效计算的同时可以实现孔内外三维动态温度场的刻画．提出的基于二元复合介质的解析模型具备很好的计算精度．10.13245/j.hust.210121.F001图1短期3 h模拟水温与实验值对比10.13245/j.hust.210121.F002图2中期40 h模拟水温与实验值对比3 动态分析选取地埋管中间深度位置即z=50 m处，进行横截面的动态温度分布研究，系统参数与模型输入均沿用第2部分中的数据．运用本模型可以快速获取钻孔内外流体与固体构成的复杂动态温度分布数据，在时间线上选择第1，10，30 h的横截面温度分布，如图3所示，得出了传统基于热阻热容模型难以获取的内部机理与传热细节，同时无须消耗大量计算资源．本部分计算在普通配有CPU-i5-6300HQ的PC上，以0.002 m的网格密度、60 s的时间步长计算40 hTRT过程仅需10 min即可获取TRT全过程每一时刻的横截面温度分布情况．10.13245/j.hust.210121.F003图3垂直双U型地埋管中段横截面温度动态分布演变(色标单位：℃)3.1　横截面温度动态演化如图3所示，在TRT开始阶段的1 h，进口与出口水温相差较小，且周围的回填土与钻孔外土壤温度上升幅度较小，但是在两个并排的进口之间的区域已经出现少量温升大于其余区域的现象．进行到第10 h时，在U型管中间靠近进口水管的部分出现了热堆积，这主要是因为进口水温比较高，而当前此横截面的回填土温度同时受到水管在不同深度方向热流的综合作用，虽然在z=50 m截面水管温度低于回填土温，但是在此截面的上游部分水温较高，而正是上游部分较高水温导致了此截面处内部热堆积．这种内部热堆积现象及进口水温的不断变化，会使水温不断上升，这也是U型管的热短路现象．从图3(b)与(c)可以看到：内部热堆积的区域在不断扩张，该截面回水管水温也在上升．此外，该垂直双U型地埋管在工作30 h后，仍然具有充足的能力对水温进行降温处理，因为从温度分布上看，在背向中心区域的方向上及四周的回填土和钻孔外周围土壤温度仍然没有大幅度上升，这其中还蕴含大量的自然冷源可以被利用．图4给出了地埋管在25与75 m深度上的横截面温度场分布情况．首先，在30 h之后，地埋管各个深度(25，50，75)m位置的回填土温度场分布具有相似性，均呈现出从两个入口管内侧向两个出口管内侧温度扩散的现象．其次，温度场直观地反映了入口管水温与出口管水温随深度增加而温差减小的现象．10.13245/j.hust.210121.F004图4垂直双U型地埋管25与75 m深度横截面温度场分布(色标单位：℃)为了更进一步分析钻孔内详细的温度变化过程，对图3中的横截面再做A-A与B-B沿线取样分析，结果如图5所示．10.13245/j.hust.210121.F005图5横截面上A-A与B-B沿线温度取样分析可以看出：在前5 h内，A-A与B-B沿线温度逐渐上升且上升的速度在逐渐下降，并且进口水管与出口水管之间的热交换不可避免，而在温度分布动态演变过程中，A-A中两者之间的温差在缩小．在B-B沿线上，因为只涉及回填土内传热，所以温度变化过程比较平滑．双U型在孔内采用的一对U型管并联方式安装，使得温度曲线在A-A与B-B上呈现非对称分布，但是总体温度曲线在时间尺度上存在极大相似性．3.2　水管沿程水温变化过程为了更进一步分析TRT过程中水管内温度沿深度方向的空间变化及在时间线上的动态变化过程，分别对进出水管水温在1，10，20，30 h进行分析，结果如图6所示，图中：红色系为进水管水温；蓝色系为出水管水温．两者组合成V字型水温曲线，并且随着时间推移，该V型曲线逐渐上升，上升速度递减，但是可以看出不同时刻的V型曲线保持着相对稳定的形状结构．10.13245/j.hust.210121.F006图6进出水管1，10，20，30 h的温度分布4 结论提出了一种垂直双U型地下换热器的动态传热解析模型，并通过该模型对垂直双U型地下换热器的传热机理进行了分析，主要结论如下．a. 提出的解析模型采用二元复合介质建模思路，结合了分段圆柱热源法及流体准动态解析模型，为垂直双U型地下换热器提供了一种计算高效、操作容易、精度较高的实用模型．该建模方法同样适应于普通单U型地埋管传热分析与优化．b. 短期3 h和中期40 h的TRT模型数据与实验数据的平均相对误差分别仅为2.25%和0.17%，均方根误差分别为0.513和0.207℃．c. 得出了内部热堆积现象及U型管之间热短路机制．分析了横截面上两条代表性沿线上温度变化过程，发现非对称温度分布和温度分布随时间演化的现象具有相似性．