地下群井改变径流条件下的三维实验模拟研究

王煜曦

（北京市地热研究院，北京 102218）

0 前言

能源和环境的可持续发展是我国近几年关注的焦点。对可再生能源与常规能源的智能耦合运行是实现绿色低碳、安全高效、城乡一体的现代能源城市的关键。据调查全世界范围内37 %的总能耗是建筑能耗，其中空调和供暖能耗占建筑能耗的50 %（Schmidt et al,1997），地源热泵系统在建筑节能源和环境可持续发展方面起着重要作用，作为清洁能源北京多采用U型地埋管地源热泵系统，但多年运行条件下地埋管周边温度场的热堆积、热泵系统的地热不平衡等问题较多，地埋管与周边岩土体包括地下水影响下的换热对地埋管换热器的运行性能尤为重要。关于垂直U型地埋管换热器的传热模型均是基于Kelvin的线热源理论（Ingetsoll et al,1954）、圆柱热源理论的纯导热模型（Bandos et al，2009；杨卫波等，2006），而在实际现场测试中发现土壤导热系数比实际值偏高（Chiasson，1992）、夏季管内流体实测温度比模拟值偏低（Witte et al,2002）等问题，范蕊等（2007）在室内搭建砂箱实验台模拟无渗流土壤、饱和土壤、有渗流土壤情况下地埋管换热器热负荷对其周边土壤温度场的影响，程金明等（2017）等基于内热源法建立了热渗耦合下的数学模型研究地下水渗流对换热管周围温度场的影响，表明地下水流越大对换热管周围温度场影响越大，龙会等（2007）等建立了考虑地下水渗流的单个钻孔三维模型，利用编程模拟过余温度场，当渗流速度一定是能有效缓解热堆积，转孔中部深度过余温度最大。以上都是针对单个线热源、圆柱热源三维温度场模拟或是结合渗流的一维模拟，模拟后均无实验验证对比。

本文研究改变流场条件下的群井换热孔周围的温度场的变化，对照实验基地实测数据进行了三维数值模拟，对比实测数据，验证模型的可靠性，同时对下一步工况进行了预测。

1 实验内容及理论研究

本次实验和数值模拟，主要针对人为改变地下流场的条件下（在抽水运行期间内）对地温场恢复情况的研究。实验共运行15天（9月8日—9月22日），前7天启动地源热泵系统，U形管开始换热，后8天关闭地源热泵系统，打开两个抽水井进行抽水试验，实验场地如图1所示。

1.1 渗流场特性

地下水在岩土体空隙中的运动称为渗流，地下含水层中三维渗流方程为：

式中：H为总水头，k ij为各个方向的渗透系数，Ss为弹性给水度或贮水率。

它表示在达西流动条件下单位体积、单位时间的水均衡关系。达西定律（或其他基本渗流定律）所描述的是水相对于固体骨架的流动。模拟地下渗流场特性时采用COMSOL里面的达西定律物理模块，包含流体参数和多孔介质参数（如有效孔隙度等），公式如下：

式中：ρ为流体密度，εP为孔隙率，k为多孔介质渗透率，P为压力，u为达西流速，Qm为质量源项（单位：kg/m3s）。

图1 实验场地平面图Fig.1 Layout of the experimental site

1.2 温度场特性

热能传递过程分为稳态和非稳态，温度场内任一点的温度不随时间变化而变化为稳态温度场，物体中个点的温度是时间的函数的热传递是非稳态热传递，即某一时刻的温度场，用数学描述为温度对时间的函数，本次模拟采用非稳态温度场。热传递分为热对流、热辅射和热传导。本次模拟只涉及到热传导和热对流，模拟埋管与不同地层进行地下换热的过程，采用COMSOL软件中的非等温管道流和多孔介质传热这两个物理场进行模拟，COMSOL软件里的非等温管道流模型流动方程如下：

式中：A为管道横截面积，u为管道里延管道曲线切线方向的流体流速，F为体积力，如重力，f D为达西摩擦因子，ρ为流体密度，d h为水力直径。

COMSOL多孔介质传热模块，依据流体扩散方程和热力学性质来说明固体骨架中有流体流动是的热力学现象。适用于流体处于稳定流状态的多孔介质传热，公式如下：

式中：ρ为流体密度，CP为流体恒压热容，q为热通量，u为流体流速，keff为有效热导率，Q为热源。

地下水的流动是一个复杂的过程，它伴随着能量的迁移和转化，地下水的渗流运动促成了岩体内热能以对流的方式发生转移，使岩体内温度场得以重新分布。当换热管与地层换热时，含水层实质上是复杂的渗流场、温度场的动力耦合，处于一种复杂的动态变化过程之中。

2 数值模型的建立

考虑到地下流场的方向和边界问题，模型区域为长方体，长124 m，宽135.5 m，高150 m，囊括张家湾实验基地。

2.1 几何模型的建立

模型简化假设：

（1）将土壤视为一个均匀、各向同性的饱和多孔介质，热物性不变。（2）忽略地下水纵向渗流影响，仅考虑水平方向层流流动。（3）抽水井壁与原状土壤相同，物性参数一致。（4）忽略U型换热孔周边回填土，假设换热孔直接与地层土壤接触。（5）考虑U型管壁层（PE管）导热系数。（6）忽略渗流路径迂回曲折，仅分析主要流向，土粒空为渗流所充满。

模型基于钻孔资料概化为5层具有不同热性能和水力属性的地质层，相关含水层已完全饱和，顶部和底部均为弱透水层，中间有一层弱透水层（隔水层）。地层分布及参数设定如表1所示，含水层孔隙率均为0.39，水力传导率2.94 ×10-4。

表1 地层分布及参数设置Tab.1 Stratigraphic distribution and parameter settings

2.2 参数设定

U型地埋管参数设置：换热孔区域并排设置有9组U型换热管，管径3.2 cm，壁厚0.3 cm，间距5m，筒壁导热系数k=0.41 W/(m·K)。初始条件：管中循环液为水，入口温度22.8℃，入口流速取实际流量1.4m3/h，入口流速设为分段函数，前7天流量为1.4 m3/h，后8天关闭热泵机组，系统停止运行，流量为0。管内壁温度与循环液体温度相同，管外壁温度与地层温度相同。

抽水井参数、初始条件设置：两口抽水井位置与工区位置相同，如图2所示，距离地埋管换热孔近的为7号抽水井，内径为22.5 cm深30 m，抽水量为10 m3/h，较远处为6号抽水井，内径为58 cm深120 m，抽水量为24 m3/h，前7天不抽水（抽水量为0），后8天开始抽水。

水的物性系数为：动力粘度μ=1×10-3N·s/m2，密度ρ=998 kg/m3，比定压热容C=4200 J/(kg·℃),导热系数k=0.58 W/(m·K)。

地层初始条件：温度采用取换热孔H5（中间换热孔）9月8日的各层温度值作为初始地温场分布函数带入模型，初始地温随深度函数如图2所示。

模型边界条件：图2中A、B面为地下水入口边界，C、D面是地下水出口边界，模型边界为隔热边界。

网格划分：为了精确计算，对U形管进行单独划分，网格最大单元尺寸为0.3 m，剩余实体采用较细化的自由四面体网格划分，共996836个域单元。

图2 模型和初始地温场Fig.2 The physical model and the initial geothermal field

3 数值模拟和实验验证

模型模拟时间是15天，前7天启动地源热泵系统，U形管开始换热，之后关闭地源热泵系统，打开两个抽水井进行抽水试验，抽水运行8天，实验共运行15天。初始地温取换热孔H5处的温度，文中所用数据都基于换热孔H5的计算结果与实验值。图3是初始地温场、初始U型地埋管道温度场，初始地埋管温度为22.8 ℃，没有与地温场进行换热。

3.1 数值模拟和实验验证

地源热泵系统运行第6天时（0～6天）地下5m、48m、89m处U形换热管附近温度和地下温度场如图4，由图可见：地层温度较原始地温场略有提高，总体上看越靠近地埋管附近，温度越高，温度场以地埋管为中心，环形递减分布。

图3 模型初始场Fig.3 The initial temperature field of the model

图4 地下5 m、48 m、89 m温度场分布图Fig.4 The temperature distribution at 5，48，89 m underground

由图5可见，随着地源热泵的开启，换热孔H5周边的温度逐渐上升，在地下5m深的换热孔周边0.5m处温度比原始地温升高1.928℃，周边1m处比原始地温升高1.1℃。地下48m深的换热孔周边0.5m处温度比原始地温升高2.386℃，周边1 m处比原始地温升高1.042℃。地下89m深的换热孔周边0.5m处温度比原始地温升高2.189℃，周边1m处比原始地温升高1.14℃。

图5 H5换热孔5 m、48 m、89 m深周边温度变化图Fig.5 The temperature around H5 buried pipes at 5, 48, 89 m underground

地源热泵系统关闭，抽水井运行8天，地下5m、48m、89m 处U形换热管附近温度和地下温度场恢复如图6。当开启两个抽水井后地下水的渗流速度增大，由于受到原有地下流场方向的影响，距离7#抽水井最近的换热孔H7，H8周围的温度场微偏向抽水井方向，其余温度场向地下水流动的下游扩散，扩散方向与地下水径流方向一致。渗流速度越大，热量传播速度就越快，温度场的热影响半径范围就越广。

由于6#抽水井较远，且抽水量不大，换热孔周围的温度场沿抽水井方向和渗流方向的和速度场方向扩散，温度场被拉伸，H7，H8，H9温度场偏向抽水井方向。

由图7可见，在抽水井运行过程中地埋管换热孔周边的温度逐步恢复，在地下5m深（原始地温15.8℃）的换热孔周边0.5m处温度恢复到16.793℃，周边1m处温度恢复到16.842℃，与原始地温分别相差0.993℃、1.042℃；地下48 m深（原始地温13.844℃）的换热孔周边0.5m处温度恢复到14.628℃，周边1 m处温度恢复到14.624℃，与原始地温分别相差0.784℃、0.78℃；地下89m深（原始地温14.43℃）的换热孔周边0.5 m处温度恢复到15.094℃，周边1m处温度恢复到15.113℃，与原始地温分别相差0.664℃、0.683℃。

实际监测中发现在抽水井运行8天后地埋管周边温度难以恢复到原始地温场，均有1℃左右的余热，根据导热微分方程可知温差变小后温度梯度减小，导热变慢，达到热平衡时间加长，模拟值与实际值相似，验证了模型的准确性。

图8为换热管工作期间U形管出口温度计算值与实际值对比图（9月11日实验区断电），由图可见出口温度的计算值基本在稳定在20℃左右，地埋管出口温度实际值在9月10日，12日较高，由于环境温差，地埋管间歇式关停工作，导致温度与模拟值有差别但与系统运行稳定的实际温度值相差不到1℃，个别天几乎完全拟合，模拟值比较精确。

图9为H5不同深度温度计算值和实测值对比图，孔深10m处15天后计算值与比实际数据略高0.95℃（实际值14.4℃，计算值15.35℃）；孔深60m处计算值与比实际数据略高2 ℃，（计算值15.921℃，实际值13.9℃）；孔深80 m处计算值与比实际数据略高1.8℃，（计算值15.834℃，实际值14℃）；孔深120 m处实验值和计算值对比图，计算值与实际数据吻合较好（计算值15.269℃，实际值15.3℃）。

综 上，H5换 热 孔10m、60m、80m、120m处的温度场计算值与实际值相比，地埋管实际温度和计算值变化趋势相同，关闭地源热泵后实测地温较计算值恢复快。在10m处和120 m处两者吻合较好。60m和80m处的温度场相差2℃以内，由于60m和80m均处于模型的隔水层，实际工作区里，在50～100 m存在细砂中砂多层砂结构，单层厚度10～14m，模型中简化后各参数取平均值，导致温度较实际值较高。U型换热孔处的回填料是膨润土和原浆，回填料的导热性与地层不同，导致地埋管换热结束后模型模拟值的恢复较实际慢。总体地埋管与周围土体换热趋势和温差与实际值差别不大，可以对实际进行预测。

图6 抽水结束后地下5 m、48 m、89 m温度场分布图Fig.6 The temperature distributions at 5 m, 48 m and 89 m underground after pumping

图7 地埋管不同深度周边温度恢复图Fig.7 The temperature around buried pipes with different depths

图8 地埋管出口温度计算值与实际值对比图Fig.8 Comparison of the simulation temperature values and the experimental values at the outlet of the U-vertical buried pipes

图9 地埋管不同深度温度计算值与实际值对比图Fig.9 Comparison of the simulation temperature values and the experimental values at different depths of buried pipes

3.2 预测加大抽水情况下地温场恢复

假设7号抽水井，抽水量为60m3/h，较远处6号抽水井抽水量为120m3/h时，地埋管换热区温度漂移如图10所示，浅层地温场由于受7#抽水井的影响，温度扩散明显偏向抽水井方向。在地下89m处，由于受6#抽水井的强烈影响，抽水运行初期地埋管温度场偏向抽水井方向，随后温度漂移较快，地温场恢复明显，有温度云飘向抽水井方向。

由图11可见，加大抽水量后，5m深的换热孔（原始地温15.8℃）周边0.5m处温度恢复到16.703℃，周边1 m处温度恢复到16.702℃，与原始地温分别相差0.903℃、0.902℃，与3.1节模拟值相比，地温场降低0.09℃、0.14℃。48 m深的换热孔周边0.5m处温度恢复到14.422℃，周边1m处温度恢复到14.418℃，与原始地温（13.844℃）相差0.578℃、0.574℃，与3.1节模拟值相比，地温场降低0.206℃。89m深的换热孔周边0.5m处温度恢复到14.954℃，周边1m处温度恢复到14.95℃，与原始地温（13.844℃）分别相差0.524 ℃、0.52℃。与3.1节模拟值相比，地温场降低0.14℃、0.163℃。

根据模拟结果得出：模拟加大抽水量后地埋管周边温度恢复并不明显，但根据图7和图11的黑色趋势线可以看出：加大抽水量后地温下降的变化率变大，但若要地温场恢复到原始状态，在地下径流速度一定的情况下，再次加大径流速度效果并不明显，需增加地温场的恢复时间。

温度恢复初期抽水井在地层抽水可假设成为管内强制对流传热（含水层为块体，孔隙为管道），根据对流传热Dittus-Boelter公式可知加大抽水量即加大了强化传热过程，所以温度变化率加大，在抽水后期热传导成为温度恢复的主要因素，换热管和周地层温度变化率减小，导热变慢，温度场恢复所需时间变长。

图10 抽水结束后地下5 m、48 m、89 m温度场分布图Fig.10 The temperature distribution at 5 m, 48 m and 89 m underground after pumping

图11 换热孔不同深度周边温度恢复图Fig.11 The recovery of temperatures around the buried pipes at different depth

4 结论

（1）U型地埋管工作期间越靠近地埋管附近，温度越高，温度场以地埋管为中心，环形递减分布，抽水井开启后地温恢复，换热孔温度场沿抽水井方向和渗流方向的和速度场方向扩散，扩散方向与地下水径流方向一致。渗流速度越大，热量传播速度就越快，温度场的热影响半径范围就越广。

（2）COMSOL三维模拟抽水井抽水后期U型管周边地温场有1℃左右的余热没有恢复与实测值类似；地埋管不同深度的温度实测值和计算模拟值变化趋势相同，但关闭地源热泵后地埋管实测温度比模拟计算值恢复快，初步分析原因：模型相关地层概化较粗、忽略考虑回填料的热导率导致；地埋管模拟出口温度模拟值与实际值相似。总体上模型计算值与实测值误差小，可用来做进一步的预测研究。

（3）在地下径流速度一定的情况下，加大抽水井抽水量地埋管周边温度恢复并不明显，但地温下降的变化率增大。恢复到原始地温需增加地温场的恢复时间。加大抽水量即加大了强化传热过程，温度变化率加大，但在抽水后期热传导成为温度恢复的主要因素，换热管和周地层温度梯度减小，导热变慢，温度场恢复所需时间变长。