U型地源热泵系统热-流-力学耦合模型研究

赵 军，刘泉声，张程远

（1. 安徽理工大学 土木建筑学院，安徽 淮南 232001；2. 中国科学院武汉岩土力学研究所 岩土力学与工程国家重点实验室，武汉 430071）

1 研究背景

地下水源热泵技术是一种采集浅层低温地能，同时满足供暖和制冷的需求，并且实现零污染排放的能源利用方式。2005年该技术被中国建设部列为建筑业十项新技术之一，在建筑物中的推广应用是国家列为节约资源节约工作重点之一，同时，许多地方都把发展地下水源热泵作为发展本地经济的一个契机，对“节能减排”和“两型社会”建设具有建设性的意义。在地下水源丰富的地区，如大江大河流域，地下水源热泵是可以普遍采用的一种地源热泵形式。地下水源热泵利用了地下水，必须涉及取用水的回灌，不回灌可以避免地下热积累对热泵系统效能的影响，但是，只取水不进行有效回灌或回灌不慎，都会造成地面沉降和已有地下管线的破坏。当前水资源税的开征，表明国家正着力改变长期以来实行资源和环境无价制度导致资源和环境恶化的现状，通过征收资源与环境税，体现资源和环境的价值。

关于单井系统（见图 1），国内外不少专家和学者[1－8]有相当多的贡献。在渗流场研究方面，李旻等[9]对单井回灌的含水层渗流场给出了解析解，为数值模拟分析奠定了基础。何满潮等[10－12]针对地热水，对井回灌渗流场中的渗透系数进行了研究，得出受地下水温度的影响，井体周围产生的物理堵塞，使得井体周围的渗透系数减小，增加了回灌的难度，并得出 K=K0e-λt（K0为初始孔隙系数；t为时间变量；λ为常数）这一衰减方程。虽然国内的不少学者在地源热泵井的回灌问题上做了大量的工作，但关于地下水源热泵系统热-温度-力全耦合研究文献鲜见，仅有少量以将热泵井群与含水层作为一个整体系统研究其整体储热性能这样的理念为基础的研究。

图1 单井回灌系统简示意图Fig.1 Stetch of standing well system

本文从能量守恒定理、质量守恒定理等原理出发，分别建立热传递方程、渗流方程、渗流方程、连续性方程以及 THM 耦合本构方程，将水源热泵单井和地下水环境作为整个系统来进行多场的耦合分析研究，利用数值软件对 THM 耦合模型给出定性的评价方法。

基本假定：①岩（砂）体是均质的，各向同性材料；②地下水流动服从Brinkman方程，Darcy定律；③岩体的变形为小变形，连续。

2 控制方程的建立

2.1 热传递方程的建立

单井热泵系统是利用地下水和地表温差提取能量。在夏季，由于地下水水温低于地表水，可以利用它来降温，而冬季则恰好相反。被利用过的冷（热）水往往都是要回灌到地下，它与原有地下水一般来说存在着温度差异，将引起地下水形成温度峰面，势必对再次抽取地下水将产生影响。

引出参数有效热传导率keff用以定性的描述热传递过程：

式中：n为岩体的孔隙率；kl、ks分别为流体和固体的热传导率。控制体的能量平衡方程为

式中：ρl为流体的密度；Cl为流体的比热容；Cs为固体的比热容；T为温度； vx、vy、vz分别为流体在x、y、z方向的流体速度。

以M和N来分别替代等式（2）左端系数，即

式中：Q为汇源项。

利用坐标矢量关系后可进一步简化为

在单井系统中，往往更多的时候是关心径向的热传递，由于考虑了岩体的各向同性性，式（5）可进一步简化为

式中：ρs为固体骨架的密度；V为流体的速度。

式（6）就是单井系统中热传递的能量表达式。

2.2 流动方程

2.2.1 区域流动方程

在利用水源热泵井的过程中，抽取的地下水的流动速度的快慢对利用能量的过程起到关键的作用，接近井体壁处，流体速度较大，而远离井体处流体的速度相对较小，因此，对井壁处和远离井壁处分别建立不同的流动方程加以描述：

式中：∇为拉普拉斯算子；η为黏度系数（kg/(ms)）；u为流体对固体颗粒的相对速度矢量（m/s）；ui,j、uj,i分别为I、j向的位移分量；k为渗透率（m2）；p为水压力（kPa）。

2.2.2 连续性方程

地下水运动的连续性方程可以从质量守恒原理出发来考虑，即渗流场中水在某一单元体内的增减速率等于进出该单元流量速率之差，分别得到流体和固体的质量守恒方程[13－14]为

式中：s为饱和度；ρf为流体密度；vf为流体速度，其他符号意义同前。

2.3 应力本构方程的建立

2.3.1 有效应力原理

岩体中饱和-非饱和岩体的有效应力原理[8]可表示为式中：εp为孔隙水压力引起的应变；εw为岩体裂隙吸水后膨胀引起的应变；βs为固相线热膨胀系数；βw为吸水线膨胀系数；p为孔隙水压力和空气压力（kPa）；C为岩体的弹性刚度张量；m为法向应力单位矩阵；Ks为固相的压缩系数。

2.3.2 几何方程

根据假定③可得几何方程和体积应变：

式中：u为岩体骨架的位移量；εv为体应变。

2.3.3 岩/砂体中力学平衡方程

单井回灌系统中，由于冷（热）水的抽取和回灌都将对岩体的内部应力（ Fi）状态产生影响，而对其主要贡献的是水头压力的变化和温度的变化。

将式（10）中第二、三式代入式（12）可得

因孔隙水压力是时间的函数，对其求偏微分可得

对式（14）取对时间的偏微分，忽略饱和度对时间变化的影响：

将式（14）代入式（15），可得

式中：h为水头高度；γ为流体重度；c为常变量；pk为孔隙水压力盒空气压力（kPa）。

对式（17）求关于时间t的偏微分：

根据假定③可得

式中：β1为热膨胀系数。

将式（17）代入式（19），可得

再将式（20）代入式（16），经过组合整理后，可得

式（21）中只含有温度T、水头压力h和位移分量u以及内力Fi的本构方程。

为便于说明数值模拟的需求，下面给出单井系统岩土体破坏准则，考虑到实际工程，特给出三维状态下的库仑破坏准则：

式（23）中N和Q可分别用下式来表达：

式中：so、φ分别为黏聚力和内摩擦角。

通过式（23）可以看出，破坏准则是通过 fail来表达的，具体的判断准则可表述为：当 fail = 0时，单井周围的岩体开始破坏；当fail＞0时，单井周围岩体处于稳定的状态；当fail＜0时，单井周围岩体处于破坏状态。

3 工程实例

以武汉某小区作为工程为例，施工前对场地的地质状况进行了解，特别注意是否有地下管线及其准确位置。对地面进行清理，铲除地面杂草、杂物和浮土，平整地面，确定钻孔的具体位置，图2、3分别为现场钻井图和测试井管道连接示意图。整理现场监测到的数据后，得到地下水温与时间的变化关系，图4为理论值与实测值的关系。从图4可以得出，理论值和实测值的误差控制在10%内，说明文中所提出的传热模型对预测水源热泵井的井周温度场的变化是有效的。

图2 钻孔Fig.2 Drill hole

图3 测试管道连接示意图Fig.3 Sketch of test tube

图4 理论值与实测值的比较Fig.4 Comparison between theoretical values and actual measurement values

4 数值模拟

数值模拟前，先对地质资料图进行适当的简化，经过简化后的地址柱状图依次为①杂填土层，厚6 m；②黏土层，厚24 m；③细砂层，厚45 m；④粗砂砾层，厚25 m；地下稳定水位位于6 m水头处。建立在第四系地层上，呈砂、砂砾石层的交互叠置，为了便于模型的模拟，取图5的模型，厚4 m，长5 m，宽1 m，为完全饱和砂土，各向同性弹性材料，井体贯穿地层100 m，井体半径为0.25 m，材料的力学参数见表1。

图5 井体模型（单位：m）Fig.5 The model of well (unit: m)

表1 井体材料参数Table1 Parameters of well

经过参数取值后，首先得到井体的水压力梯度同轴向距离的关系图，如图6所示。从图中可以看出，当井体半径r ＞ 3时，曲线较平缓；当0.5 ＜ r ＜1.5时，曲线将发生较大的变化，此时可以用Darcy-Brinkman方程来描述；当0 ＜ r ＜ 0.5时，曲线呈“陡峭”型，可以用N-S方程来描述；水压力随着井体半径的增加而增加，直至达到一个稳定的值。因此，本文所提到的流动耦合方程是适合本模型的。

图6 水压力与井体半径关系图Fig.6 Relationship between r and P

在实际工程中必须要考虑到井体成形后井体的整体变形图，为了便于说明问题，得到了单井的整体位移图，如图7所示，图中x、y、z分别为模拟整体位移图，单位为m。

图7 单井整体位移及变形图Fig.7 The whole displacement and deformation diagram of single-well

模型在计算过程中只考虑了半对称结构，从图7可以看出，井体的最大位移发生在井体的边缘处，最大变形量达到3.309×10-3m，在井壁处变形量接近为0。

现在对井体的破坏准则进行讨论，利用式（23）作为理论依据，对fail取不同的值后得到结果如图8所示，为无量纲。从图中可以看出，考虑到模型的对称性，因此只得出井体的两个角处产生了热位移，随着 fail值越来越小，破坏影响范围将越大，施工前对岩体的性质了解对井体是否失效起到关键作用。

图8 井体破坏准则模型图（单位：m）Fig.8 The model of the failure criterion for well (unit: m)

5 结 论

（1）通过井体的破坏准则可以得出，井体随着破坏值的变化而变化，在实际工程中应充分了解岩体的性质后才能施工。

（2）单井系统中水压力随着井体半径的增加而增加，直至达到一个稳定的值。

（3）现场温度实测值和理论值具有较好的吻合性，本文的模型对预测水源热泵井温度的变化是有效的。

[1]MILLARD A,REJIB A,CHIJIMATSU M,et al.Numerrical study of the THM effects on the near-field safety of a hypothetical nuclear waste repository—BMT1 of the DECOVALEX Ⅲ project. Part 2: Effects of THM coupling in continuous and homogeneous rocks[J]. Int. J.Rock Mech. Min. Sci.,2005,42: 731－744.

[2]MASSEI N,LACROIX M,WANG H Q,et al. Transport of particulate material and dissolved tracer in a highly permeable porous medium:comparision of the transfer parameters[J]. Journal of Contaminant Hydrology,2002,57: 21－39.

[3]STEPHEN E S. Particle transport through twodimensional,saturated porous media:influence of physical structure of the medium[J]. Journal of Hydrology,1995,167: 79－98.

[4]STEPHEN E S. The importance of the third dimension on transport through saturated porous media:case study based on transport of particles[J]. Journal of Hydrology,1996,179: 181－195.

[5]AHFIR N D,BENAMAR A,ALEMA,et al. Influence of internal structure and medium length on transport and deposition of suspended particles:a laboratary study[J].Transport Porous Media,2009,76: 589－307.

[6]BARTELDS G A,BRUINING J.,AND MOLENAAR J.The modeling of velocity enhancement in polymer flooding[J]. Transport Porous Media,1997,26: 75－88.

[7]NORIO T,KASUMI Y,GEORGE Z. Model study of the thermal storage system[J]. Geothermics,2003,32: 603－607.

[8]MASSEI N,LACROIX M,WANG H Q,et al. Transport of particulate material and dissolved tracer in a highly permeable porous medium: Comparison of the transfer parameters[J]. Journal of Contaminat Hydrology,2002,57: 21－39.

[9]李旻,刁乃仁,方肇洪. 单井回灌地源热泵承压含水层渗流解析解[J]. 山东建筑工程学院学报,2006,21(1): 1－5.LI Min,DIAO Nai-ren,FANG Zhao-hong. Analytical solution of seepage flow in a confined aquifer with a standing column well[J]. Journal of Shandong Institute of Architecture and Engineering,2006,21(1): 1－5.

[10]何满潮,李启民. 地热资源在移民小区可持续发展应用研究[J]. 太阳能学报,2005,25(2): 223－226.HE Man-chao,LI Qi-ming. Reserch of sustainable development of geothermal resources in immigrant communities[J]. Acta Energiae Solaris Sinica,2005,25(2): 223－226

[11]何满潮,刘斌,姚磊华,等. 地热水对井回灌渗流场理论研究[J]. 中国矿业大学学报,2004,33(3): 245－248.HE Man-chao,LIU Bin,YAO Lei-hua,et al. Study on theory of seepage field a round geothemal productionreinfection doublets wells[J]. Journal of China University of Mining &Technology,2004,33(3): 245－248.

[12]何满潮,刘斌,姚磊华,等. 地下热水回灌过程中渗透系数研究[J]. 吉林大学学报(自然科学版),2002,32(4):374－377.HE Man-chao,LIU Bin,YAO Lei-hua,et al. Study on hydraulic conductivity during geothermal reinjection[J].Journal of Jilin University(Earth Science Edition),2002,32(4): 374－377.

[13]孔祥言,李道伦,徐献芝,等. 热-流-固耦合渗流的数学模型研究[J]. 水动力学研究与进展,2005,20(2): 270－275.KONG Xiang-yan,LI Dao-lun,XU Xiang-zhi,et al.Study on the mathematical models of coupled therma hydrological mechanical process[J]. Journal of Hydrodynamics,2005,20(2): 270－275.

[14]张强林,王媛. 岩体 THM 耦合模型控制方程建立[J].西安石油大学学报[J]. 2007,22(2): 139－145.ZHANG Qiang-lin,WANG Yuan. Establishment of THM coupling governing equations in rock masses[J]. Journal of Xi′an Shiyou University,2007,22(3): 139－145.