水源热泵直线型异侧抽灌系统渗流场影响分析

李文溢， 周维博， 赖光东

(1.长安大学 环境科学与工程学院， 陕西 西安 710054； 2.长安大学 旱区地下水文与生态效应教育部重点实验室， 陕西 西安 710054)

1 研究背景

随着我国对绿色能源开发的大力支持，地热资源的研究与应用取得了巨大的进展，而地下水源热泵作为一种绿色环保、清洁高效的地热能开发利用技术，在国内乃至世界范围都有着广泛的应用前景[1-2]。然而该技术的推广往往会引发一些负面效应，如水源热泵抽回灌井群会引起含水层天然流场的变化，在水位降深超过一定限度时就可能诱发地裂缝、地面沉降等一系列环境地质问题[3-4]。国内外学者对此类课题进行了相关研究,倪龙等[5]分析了地下水源热泵同井回灌的动力学特征，并运用叠加原理推导了同井回灌的地下水降深理论解；Nam等[6]基于原型实验和数值模拟方法，对水源热泵系统中含水层水流运移规律进行研究；单明等[7]研究了水源热泵抽灌水对砂土地基的影响，指出抽回灌行为一定程度上会引起地基沉降和变形现象；何国锋等[8]研究了地下水源热泵周期性开采及回灌条件下浅层含水层水位变化特征，并指出周期交替抽灌方法和间歇式抽灌方式能有效地控制地面沉降。

国内外相关研究多集中在对水源热泵同井、对井抽灌模式的含水层流场特征和地面沉降量分析方面，而在实际工程应用中，一抽两灌模式的水源热泵机组(直线型异侧抽灌系统)数量居多[9]，但涉及的相应的渗流场特征、布井设计研究却较少。因此，本文以水源热泵直线型异侧抽灌系统为例，分别采用解析法、数值模拟法探讨渗流场中水位降深分布与变化特征，推导渗流场影响范围、安全距离与抽灌井距的函数关系，并得到了不同井距条件下的合理规划布井区，以期对水源热泵系统的合理布局和周边地质环境保护提供参考。

2 渗流场影响范围的分析

直线型异侧抽灌是指布设于承压含水层上的具有一眼抽水井、两眼回灌井的抽灌系统，抽水井流量为Q，回灌井流量为Q/2，两回灌井与抽水井的距离均为a，其布井示意图如图1所示。

图1 直线型异侧抽灌布井示意图

在单井抽水或回灌作用下，渗流场中会形成以抽水井为中心的水位降落漏斗或以回灌井为中心的水位抬升漏斗，当抽回灌系统共同作用时，含水层渗流场中会形成水位下降、抬升的叠加区。对于均质各向同性、水平等厚埋藏、侧向无限延伸、初始水力坡度为零的承压含水层，进行单井定流量抽水或回灌后所形成的渗流场影响范围是以井轴为圆心的圆域，其半径会随着时间不断增大，当抽水或回灌时间足够长时，含水层水流运动达到似稳定状态，渗流场影响范围(影响半径)会逐渐趋于一个稳定值，此时影响半径可看作为一个常数[10]。由上可知：直线型异侧抽灌系统中含水层渗流场的影响范围可看作各单井抽、回灌作用下水力影响范围的叠加区域，其平面示意图如图2所示。

图2中，D点为抽水井井位，A、E两点为回灌井井位，下降区(抬升区)表示抽水井(回灌井)单独作用下达到似稳定状态时所形成的降落漏斗(抬升漏斗)的范围，抽水井和回灌井的水力影响半径分别为R(BD)、r(AB)，抽灌井距为a(AD)，线段OF、OG、BC长分别为h1、h2、h3。

水源热泵工程建设中由于场地条件的制约，抽灌井距往往较小[12]，渗流场叠加区域的现象较为普遍，为与实际情形相符，a的取值范围为(R-r,R+r)，在相对误差不超过8‰的情况下[11]，可得弓形BFC的面积S1为：

(1)

则含水层渗流场中叠加区域的面积S2为：

(2)

根据几何原理可得抽回灌系统对整个渗流场的影响范围大小S3为：

(3)

在水源热泵工程中，含水层的水文地质参数以及抽、回灌的水力影响半径R、r可以通过抽灌试验获得，且α、β的值可通过抽灌井距a和R、r得出，因此在抽回灌条件确定的情况下，公式(3)中的g1、g2、g3、g4均为常数，直线型异侧抽灌系统的渗流场影响范围S3是一个仅与井距a有关的函数。为保证抽回灌系统的正常运行，水源热泵的布井工作中应当结合抽灌试验和井距值得出渗流场的影响范围，并明确在影响范围内是否有其他开采利用地下水的工程，从而根据周围建筑分布情况合理调整抽回灌井位布设和井距，以免对抽回灌系统造成干扰，引发地质环境问题。

图2 渗流场影响范围平面示意图

3 安全距离分析

在地下水源热泵的实际工程中，由于建设场地等条件的制约，抽灌井距往往较小且区域地下水流速较大，含水层的细颗粒介质更容易被水流携带而出，造成含水层介质的塌陷重组而导致地面沉降[13]，对建筑物地基稳定性构成潜在威胁，故抽回灌井应当与建筑保持一定的距离，即安全距离。再者地下水位的下降是威胁建筑物地基的主要原因，因此确定建筑安全距离前需要分析抽回灌系统作用下的渗流场水位降深分布情况。

3.1 渗流场叠加降深

根据渗流理论可知，承压含水层非稳定流单井抽水的Theis公式为[14]：

(4)

(5)

式(4)～(5)中：s为渗流场某一点水位降深，m；Q为抽水或灌水流量，m3/h；T为导水系数，m2/d；W(u)为井函数；S为贮水系数；t为抽、回灌时间，d;r为某点距抽水井的距离，m。

因为回灌是抽水的逆过程，所以适用于抽水的Theis井流公式同样适用于回灌[15-16]，且在实际工程中，水源热泵的抽回灌时间、抽回灌量等因素均能够满足Jacob公式的前提假设。

水源热泵系统运行条件下，抽、回灌井群共同作用于含水层渗流场，渗流场中任一点P(x，y)的水位符合降深叠加原理，渗流场中建立的平面直角坐标系如图3所示。其中，D点为抽水井位，A、E两点为回灌井位，抽灌井距为a，渗流场中P点距抽水井距离为r1，距回灌井距离分别为r2、r3。

当抽、回灌流量为定值时，t时刻渗流区内某点P的叠加降深可表示为：

sp=s1(r1,t)-s2(r2,t)-s3(r3,t)

(6)

式中：若sp>0，则表示P点水位下降；若sp<0，则表示P点水位抬升。

将公式(5)代入公式(6)可得P点的水位降深：

(7)

图3 叠加降深计算示意图

3.2 安全距离的计算

渗流场影响范围由水位抬升区和下降区组成，两区域内的地下水流运动状态存在差异，对建筑物地基造成的潜在威胁不同，因此分析水位抬升区和下降区的分界线，对于确定安全距离(即确定合理布井区)和布井位置就尤为重要。

根据水位叠加原理，水位抬升区和下降区的分界线上的各点水位降深为零，即sp满足：

(8)

在实际工程建设中，抽回灌流量Q和导水系数T均大于零，所以由式(8)可得：

(9)

(x2+y2)2=[(x+a)2+y2][(x-a)2+y2]

(10)

通过MATLAB软件绘制出公式(10)的曲线，如图4所示。图4中D为抽水井井位，A、E为回灌井井位，点M、N分别为曲线与x轴的交点。可以看出，分界曲线关于抽水井中心对称，在抽水井影响范围(外圆D)确定的情况下，分界线取值为定值(实线部分)。内圆D是以D(抽水井)为圆心，|DM|为半径的圆，Ⅰ、Ⅱ、Ⅲ、Ⅳ分别表示由分界线和内外圆所围成的封闭区域。

图4 分界线示意图

由图4可以看到，在抽回灌井群系统的作用下，区域Ⅲ、Ⅳ及内圆D内各点降深均大于零，形成水位下降区；区域Ⅰ、Ⅱ各点降深小于零，形成水位抬升区。地下水位抬升对地基沉降的潜在威胁很小[17]，而区域Ⅲ和区域Ⅳ虽然位于下降漏斗区，但距离抽水井较远，属于水力坡度较小、水流平缓的外围区域，其造成地基沉降的可能性很小或者沉降量在建筑安全标准之内[7]；而内圆D属于水位降落漏斗的中心区域，地下水动力过程较为活跃，渗流速度较快，区域内的含水层细颗粒介质更容易在地下水流运动中塌陷重组，进而诱发地面沉降现象，对建筑物造成的潜在威胁也就更大。通过以上分析可知：抽水井与周边建筑的距离应当大于DM(DN)，即抽水井与建筑物之间的安全距离为DM(DN)，内圆D以外的区域为建筑安全区。令公式(11)中y=0，可求得分界线与x轴交点M、N横坐标分别为：

(11)

因此安全距离的大小为：

(12)

在水源热泵的布井工作中，抽回灌井群系统对建筑的威胁主要来自于抽水井，抽水井位置的选取关系到周边建筑物地基的安全，而安全距离的推导同时考虑到场地因素和地面沉降因素。由公式(12)可以看出：在水文地质参数和抽回灌量参数确定的条件下，安全距离与抽灌井距a呈线性正相关关系，且会随着井距的增大而增大，其原因为抽灌井距增大时，抽水井和回灌井的相互干扰作用减小，渗流场中的叠加区域也随之缩减，抽水井运行所形成的水头下降漏斗的范围增大，这符合地下水流运动的一般规律。对比图2和图4可以看出：抽、回灌井影响范围的交点B、C位于渗流场叠加区域的边界，叠加降深始终为零，即分界曲线恒过B、C两点。随着井距的增大，分界曲线的左右两支的距离增大，地下水动力活跃区(内圆D)增大，抽水井应与建筑保持的安全距离也就增大。这与公式(12)的理论分析结果相符。

4 渗流场数值分析

4.1 渗流场数学模型的建立

本文以西安市北郊某水源热泵试验场地为例，建立地下水流数值模型。试验场地地形基本平坦，地貌单元属于渭河南岸河漫滩。试验场地为典型的承压含水层取水，以地面为零基准面，含水层顶底板标高分别为-40、-80 m，岩性主要为中粗砂和砂砾石。场地内有生产井3眼，为一抽两灌直线型异侧布井，抽灌井距27 m。根据试验场地的水文地质资料，承压含水层可概化为：水平等厚、均质各向同性、侧向无限延伸的二维非稳定流系统，侧向边界为定水头边界[12]，顶底板无垂向入渗补给，其二维地下渗流数学模型如下：

(13)

式中：H为承压含水层水头，m；T为导水系数，m2/d；S为贮水系数；t为抽灌时间，d；D为模拟区范围；H0为含水层初始水头，m；ri表示某点距抽回灌井的距离；Qi为抽回灌井的流量，m3/d；Γ为模拟区域四周边界。

4.2 模拟参数

根据研究区抽水试验，抽水井的影响半径约为200 m，为获得稳定的水头边界，选取300 m×300 m×40 m的区域为模拟区，垂向上岩层厚度取50 m，其中0～5 m为黏性土层，5～45 m为承压含水层，45～50 m为黏性土层，抽回灌井均为完整井，井径为0.167 m，数值模拟参数如表1所示。为提高计算精度，采用不等距剖分网格，网格间距4m，抽水井流量为-80 m3/h，两个回灌井流量均为40 m3/h。

采用有限差分软件Visual MODflow进行模拟计算，假设模拟区各层初始承压水头标高一致，均为48 m，承压含水层顶底板边界为隔水边界，四周边界为定水头边界(48 m)，对于布设在承压含水层区域的一抽两灌系统，渗流场达到稳定所需时间一般为50～60 h，且理论最佳布井间距一般为40～80 m[12]，在实际工程中，由于受到场地因素的制约，实际布井间距往往小于最佳间距，故本文选30～60 m之间的4组井间距进行模拟(间隔10 m)；同时为确保渗流场达到稳定，模拟时长设为72 h。

表1 数值模拟参数

4.3 模拟结果及分析

利用Visual MODflow模拟了不同井距条件下含水层渗流场的水头变化特征并得到了稳定的渗流场分布图，不同井距条件下的水位降深等值线如图5所示。

图5 不同井距下直线型异侧抽灌系统运行72 h水位降深等值线图(单位：m)

从图5中井距30、40、50、60 m条件下水位降深场等值线可以看出：等值线关于抽回灌井所在直线呈轴对称分布，且在抽回灌井附近区域等值线较为密集，说明该区域内水力坡度较大，地下水渗流速度较快；而距离抽回灌井较远的区域等值线较为稀疏，说明这些区域水力坡度较小，地下渗流受抽回灌井的影响也较小。从图5可以看出不同井距条件下得到的等值线图中均出现了叠加降深为0的曲线(即抬升区和下降区的分界曲线)，根据公式(15)可以得到井距30、40、50、60 m的条件下对应的安全距离分别为：21.2、28.2、35.3、42.4 m，该结果与图5中各井距下的对应值基本一致。

由图5还可以看出：不同井距条件下降深等值线并不是以抽水井或回灌井为中心的圆，而是沿着抽回灌井连线方向发生偏移，且井距越大，产生的偏移距离也就越明显，这与图2所示的水头抬升区(下降区)存在差异。其原因为图2基于Theis井流假设(含水层均质各向同性，初始水力坡度为零)，所得的单井影响范围是以抽水井(回灌井)为中心的圆，而在数值模拟中抽回灌井井之间互相干扰，造成了渗流场中地下水流速分布的不均匀，使得降深场发生偏移。对于因抽回灌井互相干扰所引起的影响范围大小的差异，作出以下分析：

根据数值模拟结果，得出不同井距下渗流场达到稳定状态时，抽回灌井的水位变化值(由于两侧回灌井的水位抬升值差异很小，故选用其平均值)。对于单井水力影响半径R、r的计算，本文采用了一般(1～10 m/d)渗透性的承压含水层影响半径的计算公式即奚哈德公式[18]：

(14)

式中：R为水力影响半径，m；sw为抽水井水位降深(或回灌井水位抬升)，m；k为渗透系数，m/d。

根据公式(14)计算抽、回灌井的影响半径，代入到公式(3)中得出不同井距下渗流场影响范围的解析结果，利用Surfer 12得到不同井距所得渗流场影响范围的数值模拟结果，如表2所示。

表2 不同井距条件下抽灌井水位变化及渗流场影响范围

表3 差异性显著检验结果

采用科学试验中广泛应用的单因素方差分析法，对两种不同方法下得到的渗流场影响范围的差异性进行显著检验。检验结果如表3所示。由表3可以看出，在置信度为95%的条件下，显著性水平P大于0.05，检验统计量F远远小于临界值Fcrit，这也就说明解析法和数值模拟法所得到的影响范围没有显著差异。结合图2、5和表3，可以得出：抽回灌井的相互影响虽然使得抽回灌系统的渗流场影响范围发生了偏移，但所引起的影响范围的差异并不显著。其原因为：在水文地质参数和抽灌量参数确定的条件下，虽然地下水流速分布不均匀造成了渗流场偏移，但渗流场影响范围的大小仅与抽回灌井距有关，所以相同井距条件下解析法和数值法所得的影响范围也就没有显著差异。结合表2和图5可以看出：随着井距增大，抽、回灌井之间的等水位线变得稀疏，抽、回灌井中的水位差和渗流场影响范围增大，其原因为：随着井距的增大，渗流场叠加区域减小，抽回灌之间的干扰作用减弱。这说明在水文地质参数和抽回灌参数确定的条件下，井距是决定渗流场影响范围的重要因素。在实际工程中，可以通过减小井距来减小渗流场影响范围，但同时也要考虑到井距减小会使得抽灌井之间的区域地下水动力循环更为活跃，对周边建筑物地基的威胁变大。

5 结 论

(1)在抽回灌量确定的条件下，井距的减小会使得抽回灌井中的水位差减小，但同时也会使得抽回灌井之间的地下水动力循环更为活跃，对周边建筑物地基的威胁变大。

(2)在水文地质参数、抽回灌量参数确定的条件下，直线型异侧抽灌系统对渗流场影响范围、建筑安全距离均与抽灌井井距呈正相关关系，渗流场影响范围会随着井距的增大而显著增大，可以通过减小井距的方法来减小影响范围。

(3)直线型异侧抽灌系统对建筑地基的威胁主要来自于抽水井，在水源热泵布设工作中应重点注意抽水井井位的选择，保证抽水井与周边建筑物保持一定安全距离。

参考文献:

[1] NISHMURA T. “Heat pumps- status and trends” in Asia and the Pacific[J]. International Journal of Refrigeration, 2002, 25(4):405-413.

[2] 武 强,孙海洲,徐生恒,等. 多流场耦合的单井循环换热地能采集井物理场分析[J].太阳能学报, 2015, 36(12):3015-3022.

[3] 刘 洋.某水源热泵抽灌试验及建筑物沉降影响分析[J].发电与空调, 2015，36(2):48-51.

[4] 马最良,姚 杨,姜益强,等.地下水源热泵若不能100%回灌地下水将是子孙后代的灾难[J].制冷技术，2007(4):5-8.

[5] 倪 龙,马最良.同井回灌地下水源热泵地下水渗流理论研究[J]. 太阳能学报, 2006, 27(12):1219-1224.

[6] NAM Y, OOKA R. Numerical simulation of ground heat and water transfer for groundwater heat pump system based on real-scale experiment[J]. Energy & Buildings, 2010, 42(1):69-75.

[7] 单 明,刘忠昌,舒昭然,等.水源热泵抽取和回灌地下水对建筑物沉降的影响[J].建筑结构, 2010，40(1):62-64.

[8] 何国峰,张 云,孙 铁,等.周期性开采和回灌条件下浅层地下水位变化特征研究——以天津滨海新区为例[J].水文地质工程地质, 2016, 43(6):27-34.

[9] 田良河.地下水源热泵系统布井方案分析[J].地质与勘探,2013,49(2):367-372.

[10] 薛禹群.关于稳定井流模型和Dupuit公式的讨论——对陈崇希教授“商榷”一文的答复[J]. 水利学报, 2011, 42(10):1252-1256.

[11] 张二艳.利用Taylor级数计算弓形面积的近似公式[J].北京印刷学院学报,2009, 17(2):72-74.

[12] 周维博,赖光东,姚炳光,等.西安地区水源热泵抽回灌参数优化试验与适用性规划研究[R].西安:长安大学,西安市水务局,2015.

[13] 廖 荣,丁跃元,刘立才.浅析水源热泵对区域水文地质条件的影响[J].山西建筑,2009, 35(4):213-214.

[14] 薛禹群,吴吉春.地下水动力学[M].第3版.北京：地质出版社,2010.

[15] 张蔚榛.地下水非稳定流计算和地下水资源评价[M].北京：科学出版社,1983.

[16] 马 聪,周维博,云 涛,等.地下水源热泵对井抽灌系统回灌能力研究[J].太阳能学报,2015, 36(3):593-598.

[17] 杨红禹,田春德.未来上海地区海平面上升对土压缩模量的影响分析[J].水文地质工程地质, 2001，28(2):25-27.

[18] 胡 月,刘国东,暴伟越,等.地下水环境影响评价中的“影响半径”及其确定[J].山西科技, 2015，30(6):59-62.