地下水源热泵运行模式对地温场效能影响研究

周念清，孔令熙，王小清，王 洋

(1. 同济大学水利工程系，上海 200092； 2. 上海浅层地热能工程技术研究中心，上海 200072)

0引 言

地下水源热泵从含水层中抽取地下水，通过热交换器实现向使用对象供暖或制冷，然后再将抽取的地下水重新注入到含水层中，从而保持地下水量的平衡[1-2]。在系统运行期间，沿地下水流方向引起的热对流要比热传导传递的热能更多[3]，因此在水文地质条件良好的含水层中，系统运行具有出色的供暖/制冷性能和效率[4]。经过热交换后的地下水在含水层中产生热羽流，当其沿着抽灌行为引起的局部地下水流到达抽水井时，就会出现热贯通现象[5]。抽水井温度的变化会降低系统性能和运行效率，甚至可能导致系统故障[6-7]。国内外学者关于地下水源热泵对地温场的影响做了大量的研究，但多数以室内试验[8-9]和数值模拟[10-11]为主，只有结合实际工程监测数据进行研究才真正具有应用价值[12-13]。COMSOL作为一款多物理场耦合分析软件，具有精度高、误差小的特点，近年来在地热能开发利用的研究当中得到了越来越广泛的应用[14-15]。

井群布局和运行模式是地下水源热泵设计重点考虑的问题，关系到系统的可持续运行[16-17]。在方案设计阶段，通过比较不同井群布局下热泵运行的地温场变化，可以优化抽灌井的布置方式[18-19]。对已投入运行的工程项目，一般在原有设计的基础上评价热泵系统的运行效果并提出改进策略。骆祖江等[20]对地下水源热泵工程热平衡发展趋势进行了分析和预测，通过增大温差的方式来进行调节，可有效缓解工程中出现的热贯通现象。对于实际工程而言，情形比较复杂，多数研究是基于假设简化模型进行探讨[21-22]。曾锦雯[23]建立了一个“三抽三灌”的对井模型，模拟一个供暖制冷周期不同情形下含水层温度场的分布情况，以此来探讨地下水源热泵系统运行模式的优化问题。此类研究虽具有普适性但缺少实例验证，且研究时间尺度较短。业界普遍认为“大温差、小流量”的运行模式能减小对地温场的影响，但结合实际工程运行策略优化方案时缺少相应的量化指标。

本研究以上海某地下水源热泵工程为例，采用COMSOL Multiphysics有限元分析工具建立地下采能区含水层热渗耦合数值模型，模拟系统长期运行时地温场的演变过程，从冷、热负荷设计和抽灌模式两个方面进行探讨，获取相关的技术参数，为地下水源热泵系统的长期稳定运行提供技术支撑。

1场地工程概况

表1场地地层结构Table 1Site Stratigraphic Texture

地下水源热泵系统主要用于农业种植温室的温度调节，温室面积约21 000 m2。该工程建设项目在2010年3月进行了水资源论证，地下水源热泵系统于2014年11月建成后开始试运行，于2015年11月正式投入使用，至今仍在稳定运行中。系统供暖季为12月～次年5月，制冷季为6～10月，全天24 h运行。供暖季日均热负荷为350 kW，制冷季日均冷负荷为262.5 kW。

工程采用“两抽两灌”的L型布井方案，井1、井2为冬灌夏抽井(冷井)，井3、井4为夏灌冬抽井(热井)，抽、灌井之间距离为124 m。水源井井深104.5 m，滤水管位于井深77.6～102.6 m处，总长度16.7 m，为承压非完整井。系统设计单井最大地下水抽灌量为30 m3·h-1。井管、滤水管均采用φ250 PVC管材。抽取的地下水100%回灌入相同目标含水层。另外，场地井2和井3之间设置了3口地温监测井，工程平面布置见图1。

图1工程平面布置Fig.1Project Plane Layout

2地下水渗流与热量运移耦合模型

2.1基本假设

地下含水层热量运移是一个极其复杂的传热传导过程，为便于分析和计算，对抽灌井计算模型做了如下假设：①计算区岩土体为饱和多孔介质，同一地层岩土体均质、各向同性，产状水平，厚度不变；②地下水流为水平径向流，忽略竖直方向上的渗流，地下水流运动服从达西定律；③不考虑上部弱透水层的贮水性和弹性释水量；④流体与介质的比热容、导热系数、密度、渗透率等物性参数不随温度和压力变化而变化；⑤忽略自然对流和辐射换热的影响；⑥忽略井身传热带来的热损失，不考虑地面环境和大地热流对采能区的温度影响。

2.2地下水渗流与热量运移数学模型

本研究的含水层类型为承压含水层，根据假设条件，含水层饱和多孔介质地下水流动遵循达西定律，并且满足质量守恒定律。地下水渗流过程和定解条件可用方程(1)描述。

(1)

(2)

热量运移与传输的机制主要通过对流和传导实现，忽略辐射换热的影响。地下水源热泵抽灌井热量运移主要发生在含水层水与水之间、水与介质之间，在此过程中地下多孔介质基体和流体成为了热传递的载体，二者构成了渗流场和温度场的耦合关系。通过达西渗流速度可将渗流方程与热传导方程进行耦合。含水层多孔介质热量运移过程和定解条件可以用方程(3)描述。

(3)

q=-keffT

(4)

2.3数值模型与网格剖分

地下水源热泵的运行过程是复杂的三维多场耦合问题，许多研究常将其简化为二维模型来计算。陈必光[24]讨论了对井系统的二维简化问题，结果表明，横向剖切模型对于模拟热储层中间深度区域的温度场效果良好。项目选定的热储层即采能区埋深为74.0～104.2 m，含水层岩性为粉细砂、中砂互层，选择渗透性较好的中砂层(91.5～94.3 m)进行模拟，取中间深度93 m地层进行横切获得二维简化模型。

采用COMSOL Multiphysics有限元软件进行模拟，研究热渗耦合作用下的地温变化特征，分别选用达西定律模块和多孔介质传热模块模拟抽灌系统的渗流和传热过程。为了确定模型的几何尺寸以消除边界的影响，先进行了多种尺寸的试算，判定依据是研究区边界上的水头和温度条件不受抽灌井引起的地下渗流场和温度场改变带来的影响，试算结果得到抽灌井的影响半径约为300 m。为了网格剖分和计算的便捷，几何形状通常采用正方形，故研究区尺寸最终确定为800 m×800 m。模型网格采用细化尺寸的自由三角形网格，最小网格单元边长为0.1 m，共包含6 496个域单元和132个边界单元，计算区域单元剖分参见图2。

图2抽灌井模型网格(单位:m)Fig.2Pumping and Recharge Well Model Grid (Unit:m)

2.4初始条件和边界条件

为了评价含水层的出水能力，并计算含水层的渗透系数，在方案阶段进行了现场抽水试验。此次采用承压水非完整井稳定流单孔抽水试验，根据经验公式(5)可计算含水层的渗透系数。

(5)

式中：Q′为抽水井流量；l为滤水管长度；rw为抽水井半径；sw为抽水井水位降深。

利用J2监测井进行单孔抽水试验，根据抽水试验期间观测数据，取平均值进行计算，最终得到目标含水层渗透系数K为26.8 m·d-1。结合室内试验结果，得到含水层主要物理和热物性参数，见表2。

表2含水层主要参数Table 2Main Parameters of Aquife

根据水文地质资料，计算区的初始水头为-1.63～-0.27 m，地下水径流方向由西向东，平均水力坡度为1.26×10-3，地下水平均渗流速度为3.89×10-7m·s-1。含水层初始温度为19.67 ℃，假定同一深度初始地温分布一致。初始水头和渗流场如图3所示。

图3初始水头和渗流场Fig.3Initial Water Head and Seepage Field

在渗流场中，含水层四周设定为第一类(水头)边界；抽灌井设定为第二类(流量)边界。为方便计算，抽灌井流量取供暖(制冷)期的平均流量，且抽灌量相等。在温度场中，含水层四周设为开放边界，即认为四周为恒温边界，上游温度设为初始地温；抽灌井设为第二类(热通量)边界，回灌水温度恒定。

2.5模型的识别与验证

模型初步建立后，需要对模型参数进行识别，并验证模型的准确性，取J1、J2、J3监测井处一个供暖制冷周期地温变化的实测值和模拟值进行对比。根据系统运行3年的监测数据，供暖期和制冷期抽灌水量变化呈正态分布，最大值分别约为500 m3·d-1和1 000 m3·d-1，平均值约为300 m3·d-1；冷井供暖期回灌水温度分别在8 ℃上下变化，制冷期抽水温度则由8 ℃逐渐升高至14 ℃；热井供暖期抽水温度基本稳定在20 ℃左右，制冷期回灌水温度受热负荷变化的影响在17～30 ℃动态变化，平均在21 ℃左右。因此，模型的单井抽灌流量设为300 m3·d-1，供暖期和制冷期回灌水温度分别设为8 ℃和21 ℃。结果显示，J1、J2、J3处地温实测值与模拟值的变化曲线较为吻合，见图4。计算可得J1、J2、J3处地温实测值和模拟值的平均误差分别为5.22%、2.60%和7.28%，说明该数值模型能较好地拟合实际情况，满足参数反演要求，可进行下一步计算分析工作。

图4J1、J2、J3处实际地温与模拟地温对比Fig.4Comparison of Actual and Simulated Ground Temperature at J1, J2 and J3

3计算结果与分析

3.1地温场特征

为了分析热泵系统运行期间地温变化规律和换热性能，共模拟了系统运行5个供暖和制冷周期含水层的地温场变化，选取3个供暖和制冷周期进行分析，见图5。由图5可知，冷、热井交替抽灌地下水的运行模式使得冷、热影响区呈现交替扩散收缩的趋势，且热羽形状受地下水渗流和抽灌作用影响呈现水滴状。“冬灌夏抽”和“夏灌冬抽”的地下水抽灌模式使得供暖期和制冷期灌入含水层的冷量和热量被储存起来，作为下一运行周期的冷源和热源，实现废热再利用，从而避免了热冗余导致热堆积区的产生。由于供暖期换热量大于制冷期，且供暖期回灌水与初始地温的温差较大，故冷井附近的冷影响区相对于热影响区其扩散速度更快、范围更大，这就导致每一运行周期冷影响区相比上一周期都在不断地扩大，因此需要考察冷影响区对于系统运行的影响。

图5供暖和制冷周期含水层地温场变化特征Fig.5Geothermal Field Change Characteristics of Aquifer During Heating and Cooling Period

在系统冷、热负荷满足温室实际需求的情况下，本研究考虑是否会发生热贯通来评价系统换热性能的优劣。取井2和井3之间的截线，设定低于初始地温0.5 ℃的19.17 ℃等温线作为冷影响区边界，当冷影响区扩散至热井，即影响距离大于井2和井3间距124 m时，认为出现热贯通现象。由模拟结果可以得出，第1、第3和第5个制冷期末的影响距离分别为55.5、80.2、100 m，平均扩散速度为11.1 m·a-1。由此可见，该热泵系统的短期运行效果良好，但预计到第7个供暖期末冷影响区将扩散到热井，发生热贯通，导致热泵换热性能降低。

3.2冷、热负荷设计的影响

在热泵系统实际运行时，日均冷、热负荷会根据用户需求进行调整。通过模拟冷、热负荷变化时含水层热量运移过程，评价其对热泵系统换热性能的影响。为了简便计算，在满足分析要素的基础上设置了“同时增加”“一增一减”“同时减小”等不同工况。根据系统实际运行时冷、热负荷变化情况，确定各工况设定的负荷增减幅度为20%。因为冷、热负荷与系统循环水量及抽灌水温差有关，所以通过同时改变日均抽灌水流量和温差，分别比较冷、热负荷增加和减小20%时的运行工况模拟得到的结果，各工况设定见表3。

表3不同冷、热负荷的工况设定Table 3Working Condition Setting of Different Cooling and Heating Loads

图6不同冷、热负荷时的冷影响区半径Fig.6Radius of Cold Affected Zone Under Different Cooling and Heating Loads

从模拟结果可以得出，在供暖期冷影响区在不同工况下仍在不断向热井扩散，如图6所示。工况1～5冷影响区平均扩散速度分别为10.8、11.9、16.9、5.8、9.9 m·a-1，预计出现热贯通现象的时间分别为2 355、2 004、1 594、4 527、2 704 d。增加热负荷会加剧热贯通现象的发生。与基础工况相比，工况3的冷影响区扩散速度提高了56.5%，预计提前2个周期发生热贯通；增加冷负荷可以明显延缓这种趋势，工况2比工况3推后1个周期达到热贯通。工况2、工况5与工况3相比，冷影响区扩散速度分别降低了29.6%和41.4%，因此冷、热负荷设计对冷影响区扩散速度是敏感的。与增加冷负荷相比，减小热负荷更能有效延缓冷影响区的扩散，工况5预计比基础工况推后1个周期发生热贯通。如果减小热负荷同时增加冷负荷，那么这种缓解作用将得到显著的提升，工况4的扩散速度相比基础工况降低了46.3%，在模拟的9个运行周期内未出现热贯通现象。另外，在第9个供暖期末，工况1～5冷影响区场地外扩散距离分别为135、140、150、118、133 m，因此，减小热负荷、增加冷负荷均能减小冷扩散羽对场地周边的影响，且有利于系统的可持续运行。

此外，为了评价冷、热负荷变化对热扩散羽的影响，设定20.17 ℃等温线作为热影响区边界，考察井2、井3截线上的热影响区在制冷期运行时的变化趋势。在不同工况下，制冷期末的热影响区随着热泵运行周期总体呈现缓慢收缩的趋势见图7。由于冷负荷增加使得热影响区扩散距离增加，工况2和工况4分别比基础工况平均增加了4.4%和26.6%。由模拟结果可知，热影响区距离在20～42 m之间，不会对下一周期冷井的运行产生影响。

图7不同冷、热负荷时的热影响区半径Fig.7Radius of Heat Affected Zone Under Different Cooling and Heating Loads

需要说明的是，如根据用户负荷需求本身对系统冷、热负荷工况进行调整，则这种调控措施完全能够满足实际工程运行的需求；如为了缓解热贯通现象进行调整，则需要其他能源的运用或项目本身运营策略的调整作为补充，以满足用户对供暖制冷的需求。

3.3抽灌模式的影响

当热泵系统冷、热负荷设计恒定时，可以通过改变系统循环水量和抽灌水温差，分析热泵抽灌模式调整对换热性能的影响。下面分别对增大循环水量、减小抽灌水温差以及减小循环水量、增大抽灌水温差的运行工况进行模拟，并比较其结果。根据监测数据可知，热泵运行时实际的单井涌水量最高可达500～1 000 m3·d-1，全年平均为300 m3·d-1，说明该含水层开采能力和水泵功率有较大的可调整余量，但由于冷、热负荷恒定，循环水流量与抽灌水温差呈反比关系，温差的限制使得流量不能过大或过小。如果实际运行中井1、井2供暖期回灌水在8 ℃上下变化，抽水温度基本稳定在20 ℃，抽灌水温差在10～15 ℃之间。因此，循环水流量的调整范围不宜超过±20%，各工况设定见表4。

表4不同抽灌模式工况Table 4Working Condition of Different Pumping and Injection Modes

图8不同抽灌模式时的冷影响区半径Fig.8Radius of Cold Affected Zone Under Different Pumping and Injection Modes

从模拟结果可知，冷、热负荷设计恒定时各工况冷影响区范围仍会到达热井处，但在程度上存在差异，如图8所示。工况6～9冷影响区扩散速度分别为12.0、13.2、9.8、9.1 m·a-1，预计发生热贯通的时间分别为2 000、1 977、2 717、3 083 d。工况6、7冷影响区扩散速度分别比基础工况提高了11.1%和22.2%，而工况8、9则分别降低了9.3%和15.7%。可以得到循环水量对冷影响区扩散速度的灵敏度为：每增加1%循环水量，冷影响区扩散速度提高1.11%；每减小1%循环水量，冷影响区扩散速度提高0.86%。由此可见，抽灌水流量增大、温差减小时，热贯通现象加剧；抽灌水流量减小、温差增大时，有利于缓解热贯通趋势。另外，减小抽灌流量还可以减小冷羽场地外方向的扩散距离，第9个供暖期末工况9的冷羽半径比基础工况减小了3.6%。因此，采用“大温差、小流量”抽灌模式能减小系统运行对地温场的影响。

随着热泵运行，各工况制冷期末的热影响区距离均在20～40 m之间，呈现先扩散后收缩的趋势(图9)，平均收缩速率为1.2 m·a-1。热影响区的发展与抽灌模式的选择无明显关系，对系统运行的影响也较小。值得注意的是，在冷、热负荷恒定时，无论循环水量和抽灌水温差如何调整，经换热器换热后提供给温室的热(冷)量同样是恒定的，所以调整水量和温差后仍能符合温室对温度的要求。

图9不同抽灌模式时的热影响区半径Fig.9Radius of Heat Affected Zone Under Different Pumping and Injection Modes

4结语

(1)采用COMSOL Multiphysics建立的热渗耦合数值模型能较为准确地模拟地下水源热泵运行地温场的演变特征。现行工况的模拟结果显示，热泵系统的短期运行效果良好，但由于全年热负荷大于冷负荷，随着热泵运行，冷影响区扩散范围将逐渐增大，冷、热源井之间会在第7个供暖期末出现热贯通现象，长期运行将不利于热井的取热。

(2)减小热负荷、增加冷负荷可以显著延缓热贯通发生的趋势。若减小20%热负荷的同时增加20%冷负荷，冷影响区的扩散速度比原来降低46.3%，在模拟周期内并未发生热贯通。虽然冷负荷增加使得热影响区扩散距离增加，但不会对冷井的运行造成影响。因此，合理调节冷、热负荷有利于热泵系统的长期稳定运行。

(3)当冷、热负荷设计恒定时，分别增大10%和20%循环水量，冷影响区扩散速度将分别比原工况提高11.1%和22.2%；如果分别减小10%和20%循环水量，则分别降低了9.3%和15.7%。说明减小抽灌水流量、增大温差的抽灌模式能有效缓解热贯通现象的发生。