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

(1. 同濟大學(xué)水利工程系，上海 200092； 2. 上海淺層地?zé)崮芄こ碳夹g(shù)研究中心，上海 200072)

0引 言

地下水源熱泵從含水層中抽取地下水，通過熱交換器實現(xiàn)向使用對象供暖或制冷，然后再將抽取的地下水重新注入到含水層中，從而保持地下水量的平衡[1-2]。在系統(tǒng)運行期間，沿地下水流方向引起的熱對流要比熱傳導(dǎo)傳遞的熱能更多[3]，因此在水文地質(zhì)條件良好的含水層中，系統(tǒng)運行具有出色的供暖/制冷性能和效率[4]。經(jīng)過熱交換后的地下水在含水層中產(chǎn)生熱羽流，當(dāng)其沿著抽灌行為引起的局部地下水流到達抽水井時，就會出現(xiàn)熱貫通現(xiàn)象[5]。抽水井溫度的變化會降低系統(tǒng)性能和運行效率，甚至可能導(dǎo)致系統(tǒng)故障[6-7]。國內(nèi)外學(xué)者關(guān)于地下水源熱泵對地溫場的影響做了大量的研究，但多數(shù)以室內(nèi)試驗[8-9]和數(shù)值模擬[10-11]為主，只有結(jié)合實際工程監(jiān)測數(shù)據(jù)進行研究才真正具有應(yīng)用價值[12-13]。COMSOL作為一款多物理場耦合分析軟件，具有精度高、誤差小的特點，近年來在地?zé)崮荛_發(fā)利用的研究當(dāng)中得到了越來越廣泛的應(yīng)用[14-15]。

井群布局和運行模式是地下水源熱泵設(shè)計重點考慮的問題，關(guān)系到系統(tǒng)的可持續(xù)運行[16-17]。在方案設(shè)計階段，通過比較不同井群布局下熱泵運行的地溫場變化，可以優(yōu)化抽灌井的布置方式[18-19]。對已投入運行的工程項目，一般在原有設(shè)計的基礎(chǔ)上評價熱泵系統(tǒng)的運行效果并提出改進策略。駱祖江等[20]對地下水源熱泵工程熱平衡發(fā)展趨勢進行了分析和預(yù)測，通過增大溫差的方式來進行調(diào)節(jié)，可有效緩解工程中出現(xiàn)的熱貫通現(xiàn)象。對于實際工程而言，情形比較復(fù)雜，多數(shù)研究是基于假設(shè)簡化模型進行探討[21-22]。曾錦雯[23]建立了一個“三抽三灌”的對井模型，模擬一個供暖制冷周期不同情形下含水層溫度場的分布情況，以此來探討地下水源熱泵系統(tǒng)運行模式的優(yōu)化問題。此類研究雖具有普適性但缺少實例驗證，且研究時間尺度較短。業(yè)界普遍認(rèn)為“大溫差、小流量”的運行模式能減小對地溫場的影響，但結(jié)合實際工程運行策略優(yōu)化方案時缺少相應(yīng)的量化指標(biāo)。

本研究以上海某地下水源熱泵工程為例，采用COMSOL Multiphysics有限元分析工具建立地下采能區(qū)含水層熱滲耦合數(shù)值模型，模擬系統(tǒng)長期運行時地溫場的演變過程，從冷、熱負(fù)荷設(shè)計和抽灌模式兩個方面進行探討，獲取相關(guān)的技術(shù)參數(shù)，為地下水源熱泵系統(tǒng)的長期穩(wěn)定運行提供技術(shù)支撐。

1場地工程概況

表1場地地層結(jié)構(gòu)Table 1Site Stratigraphic Texture

地下水源熱泵系統(tǒng)主要用于農(nóng)業(yè)種植溫室的溫度調(diào)節(jié)，溫室面積約21 000 m2。該工程建設(shè)項目在2010年3月進行了水資源論證，地下水源熱泵系統(tǒng)于2014年11月建成后開始試運行，于2015年11月正式投入使用，至今仍在穩(wěn)定運行中。系統(tǒng)供暖季為12月～次年5月，制冷季為6～10月，全天24 h運行。供暖季日均熱負(fù)荷為350 kW，制冷季日均冷負(fù)荷為262.5 kW。

工程采用“兩抽兩灌”的L型布井方案，井1、井2為冬灌夏抽井(冷井)，井3、井4為夏灌冬抽井(熱井)，抽、灌井之間距離為124 m。水源井井深104.5 m，濾水管位于井深77.6～102.6 m處，總長度16.7 m，為承壓非完整井。系統(tǒng)設(shè)計單井最大地下水抽灌量為30 m3·h-1。井管、濾水管均采用φ250 PVC管材。抽取的地下水100%回灌入相同目標(biāo)含水層。另外，場地井2和井3之間設(shè)置了3口地溫監(jiān)測井，工程平面布置見圖1。

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

2地下水滲流與熱量運移耦合模型

2.1基本假設(shè)

地下含水層熱量運移是一個極其復(fù)雜的傳熱傳導(dǎo)過程，為便于分析和計算，對抽灌井計算模型做了如下假設(shè)：①計算區(qū)巖土體為飽和多孔介質(zhì)，同一地層巖土體均質(zhì)、各向同性，產(chǎn)狀水平，厚度不變；②地下水流為水平徑向流，忽略豎直方向上的滲流，地下水流運動服從達西定律；③不考慮上部弱透水層的貯水性和彈性釋水量；④流體與介質(zhì)的比熱容、導(dǎo)熱系數(shù)、密度、滲透率等物性參數(shù)不隨溫度和壓力變化而變化；⑤忽略自然對流和輻射換熱的影響；⑥忽略井身傳熱帶來的熱損失，不考慮地面環(huán)境和大地?zé)崃鲗Σ赡軈^(qū)的溫度影響。

2.2地下水滲流與熱量運移數(shù)學(xué)模型

本研究的含水層類型為承壓含水層，根據(jù)假設(shè)條件，含水層飽和多孔介質(zhì)地下水流動遵循達西定律，并且滿足質(zhì)量守恒定律。地下水滲流過程和定解條件可用方程(1)描述。

(1)

(2)

熱量運移與傳輸?shù)臋C制主要通過對流和傳導(dǎo)實現(xiàn)，忽略輻射換熱的影響。地下水源熱泵抽灌井熱量運移主要發(fā)生在含水層水與水之間、水與介質(zhì)之間，在此過程中地下多孔介質(zhì)基體和流體成為了熱傳遞的載體，二者構(gòu)成了滲流場和溫度場的耦合關(guān)系。通過達西滲流速度可將滲流方程與熱傳導(dǎo)方程進行耦合。含水層多孔介質(zhì)熱量運移過程和定解條件可以用方程(3)描述。

(3)

q=-keffT

(4)

2.3數(shù)值模型與網(wǎng)格剖分

地下水源熱泵的運行過程是復(fù)雜的三維多場耦合問題，許多研究常將其簡化為二維模型來計算。陳必光[24]討論了對井系統(tǒng)的二維簡化問題，結(jié)果表明，橫向剖切模型對于模擬熱儲層中間深度區(qū)域的溫度場效果良好。項目選定的熱儲層即采能區(qū)埋深為74.0～104.2 m，含水層巖性為粉細(xì)砂、中砂互層，選擇滲透性較好的中砂層(91.5～94.3 m)進行模擬，取中間深度93 m地層進行橫切獲得二維簡化模型。

采用COMSOL Multiphysics有限元軟件進行模擬，研究熱滲耦合作用下的地溫變化特征，分別選用達西定律模塊和多孔介質(zhì)傳熱模塊模擬抽灌系統(tǒng)的滲流和傳熱過程。為了確定模型的幾何尺寸以消除邊界的影響，先進行了多種尺寸的試算，判定依據(jù)是研究區(qū)邊界上的水頭和溫度條件不受抽灌井引起的地下滲流場和溫度場改變帶來的影響，試算結(jié)果得到抽灌井的影響半徑約為300 m。為了網(wǎng)格剖分和計算的便捷，幾何形狀通常采用正方形，故研究區(qū)尺寸最終確定為800 m×800 m。模型網(wǎng)格采用細(xì)化尺寸的自由三角形網(wǎng)格，最小網(wǎng)格單元邊長為0.1 m，共包含6 496個域單元和132個邊界單元，計算區(qū)域單元剖分參見圖2。

圖2抽灌井模型網(wǎng)格(單位:m)Fig.2Pumping and Recharge Well Model Grid (Unit:m)

2.4初始條件和邊界條件

為了評價含水層的出水能力，并計算含水層的滲透系數(shù)，在方案階段進行了現(xiàn)場抽水試驗。此次采用承壓水非完整井穩(wěn)定流單孔抽水試驗，根據(jù)經(jīng)驗公式(5)可計算含水層的滲透系數(shù)。

(5)

式中：Q′為抽水井流量；l為濾水管長度；rw為抽水井半徑；sw為抽水井水位降深。

利用J2監(jiān)測井進行單孔抽水試驗，根據(jù)抽水試驗期間觀測數(shù)據(jù)，取平均值進行計算，最終得到目標(biāo)含水層滲透系數(shù)K為26.8 m·d-1。結(jié)合室內(nèi)試驗結(jié)果，得到含水層主要物理和熱物性參數(shù)，見表2。

表2含水層主要參數(shù)Table 2Main Parameters of Aquife

根據(jù)水文地質(zhì)資料，計算區(qū)的初始水頭為-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

在滲流場中，含水層四周設(shè)定為第一類(水頭)邊界；抽灌井設(shè)定為第二類(流量)邊界。為方便計算，抽灌井流量取供暖(制冷)期的平均流量，且抽灌量相等。在溫度場中，含水層四周設(shè)為開放邊界，即認(rèn)為四周為恒溫邊界，上游溫度設(shè)為初始地溫；抽灌井設(shè)為第二類(熱通量)邊界，回灌水溫度恒定。

2.5模型的識別與驗證

模型初步建立后，需要對模型參數(shù)進行識別，并驗證模型的準(zhǔn)確性，取J1、J2、J3監(jiān)測井處一個供暖制冷周期地溫變化的實測值和模擬值進行對比。根據(jù)系統(tǒng)運行3年的監(jiān)測數(shù)據(jù)，供暖期和制冷期抽灌水量變化呈正態(tài)分布，最大值分別約為500 m3·d-1和1 000 m3·d-1，平均值約為300 m3·d-1；冷井供暖期回灌水溫度分別在8 ℃上下變化，制冷期抽水溫度則由8 ℃逐漸升高至14 ℃；熱井供暖期抽水溫度基本穩(wěn)定在20 ℃左右，制冷期回灌水溫度受熱負(fù)荷變化的影響在17～30 ℃動態(tài)變化，平均在21 ℃左右。因此，模型的單井抽灌流量設(shè)為300 m3·d-1，供暖期和制冷期回灌水溫度分別設(shè)為8 ℃和21 ℃。結(jié)果顯示，J1、J2、J3處地溫實測值與模擬值的變化曲線較為吻合，見圖4。計算可得J1、J2、J3處地溫實測值和模擬值的平均誤差分別為5.22%、2.60%和7.28%，說明該數(shù)值模型能較好地擬合實際情況，滿足參數(shù)反演要求，可進行下一步計算分析工作。

圖4J1、J2、J3處實際地溫與模擬地溫對比Fig.4Comparison of Actual and Simulated Ground Temperature at J1, J2 and J3

3計算結(jié)果與分析

3.1地溫場特征

為了分析熱泵系統(tǒng)運行期間地溫變化規(guī)律和換熱性能，共模擬了系統(tǒng)運行5個供暖和制冷周期含水層的地溫場變化，選取3個供暖和制冷周期進行分析，見圖5。由圖5可知，冷、熱井交替抽灌地下水的運行模式使得冷、熱影響區(qū)呈現(xiàn)交替擴散收縮的趨勢，且熱羽形狀受地下水滲流和抽灌作用影響呈現(xiàn)水滴狀?！岸嘞某椤焙汀跋墓喽椤钡牡叵滤楣嗄Ｊ绞沟霉┡诤椭评淦诠嗳牒畬拥睦淞亢蜔崃勘粌Υ嫫饋?，作為下一運行周期的冷源和熱源，實現(xiàn)廢熱再利用，從而避免了熱冗余導(dǎo)致熱堆積區(qū)的產(chǎn)生。由于供暖期換熱量大于制冷期，且供暖期回灌水與初始地溫的溫差較大，故冷井附近的冷影響區(qū)相對于熱影響區(qū)其擴散速度更快、范圍更大，這就導(dǎo)致每一運行周期冷影響區(qū)相比上一周期都在不斷地擴大，因此需要考察冷影響區(qū)對于系統(tǒng)運行的影響。

圖5供暖和制冷周期含水層地溫場變化特征Fig.5Geothermal Field Change Characteristics of Aquifer During Heating and Cooling Period

在系統(tǒng)冷、熱負(fù)荷滿足溫室實際需求的情況下，本研究考慮是否會發(fā)生熱貫通來評價系統(tǒng)換熱性能的優(yōu)劣。取井2和井3之間的截線，設(shè)定低于初始地溫0.5 ℃的19.17 ℃等溫線作為冷影響區(qū)邊界，當(dāng)冷影響區(qū)擴散至熱井，即影響距離大于井2和井3間距124 m時，認(rèn)為出現(xiàn)熱貫通現(xiàn)象。由模擬結(jié)果可以得出，第1、第3和第5個制冷期末的影響距離分別為55.5、80.2、100 m，平均擴散速度為11.1 m·a-1。由此可見，該熱泵系統(tǒng)的短期運行效果良好，但預(yù)計到第7個供暖期末冷影響區(qū)將擴散到熱井，發(fā)生熱貫通，導(dǎo)致熱泵換熱性能降低。

3.2冷、熱負(fù)荷設(shè)計的影響

在熱泵系統(tǒng)實際運行時，日均冷、熱負(fù)荷會根據(jù)用戶需求進行調(diào)整。通過模擬冷、熱負(fù)荷變化時含水層熱量運移過程，評價其對熱泵系統(tǒng)換熱性能的影響。為了簡便計算，在滿足分析要素的基礎(chǔ)上設(shè)置了“同時增加”“一增一減”“同時減小”等不同工況。根據(jù)系統(tǒng)實際運行時冷、熱負(fù)荷變化情況，確定各工況設(shè)定的負(fù)荷增減幅度為20%。因為冷、熱負(fù)荷與系統(tǒng)循環(huán)水量及抽灌水溫差有關(guān)，所以通過同時改變?nèi)站楣嗨髁亢蜏夭睿謩e比較冷、熱負(fù)荷增加和減小20%時的運行工況模擬得到的結(jié)果，各工況設(shè)定見表3。

表3不同冷、熱負(fù)荷的工況設(shè)定Table 3Working Condition Setting of Different Cooling and Heating Loads

圖6不同冷、熱負(fù)荷時的冷影響區(qū)半徑Fig.6Radius of Cold Affected Zone Under Different Cooling and Heating Loads

從模擬結(jié)果可以得出，在供暖期冷影響區(qū)在不同工況下仍在不斷向熱井?dāng)U散，如圖6所示。工況1～5冷影響區(qū)平均擴散速度分別為10.8、11.9、16.9、5.8、9.9 m·a-1，預(yù)計出現(xiàn)熱貫通現(xiàn)象的時間分別為2 355、2 004、1 594、4 527、2 704 d。增加熱負(fù)荷會加劇熱貫通現(xiàn)象的發(fā)生。與基礎(chǔ)工況相比，工況3的冷影響區(qū)擴散速度提高了56.5%，預(yù)計提前2個周期發(fā)生熱貫通；增加冷負(fù)荷可以明顯延緩這種趨勢，工況2比工況3推后1個周期達到熱貫通。工況2、工況5與工況3相比，冷影響區(qū)擴散速度分別降低了29.6%和41.4%，因此冷、熱負(fù)荷設(shè)計對冷影響區(qū)擴散速度是敏感的。與增加冷負(fù)荷相比，減小熱負(fù)荷更能有效延緩冷影響區(qū)的擴散，工況5預(yù)計比基礎(chǔ)工況推后1個周期發(fā)生熱貫通。如果減小熱負(fù)荷同時增加冷負(fù)荷，那么這種緩解作用將得到顯著的提升，工況4的擴散速度相比基礎(chǔ)工況降低了46.3%，在模擬的9個運行周期內(nèi)未出現(xiàn)熱貫通現(xiàn)象。另外，在第9個供暖期末，工況1～5冷影響區(qū)場地外擴散距離分別為135、140、150、118、133 m，因此，減小熱負(fù)荷、增加冷負(fù)荷均能減小冷擴散羽對場地周邊的影響，且有利于系統(tǒng)的可持續(xù)運行。

此外，為了評價冷、熱負(fù)荷變化對熱擴散羽的影響，設(shè)定20.17 ℃等溫線作為熱影響區(qū)邊界，考察井2、井3截線上的熱影響區(qū)在制冷期運行時的變化趨勢。在不同工況下，制冷期末的熱影響區(qū)隨著熱泵運行周期總體呈現(xiàn)緩慢收縮的趨勢見圖7。由于冷負(fù)荷增加使得熱影響區(qū)擴散距離增加，工況2和工況4分別比基礎(chǔ)工況平均增加了4.4%和26.6%。由模擬結(jié)果可知，熱影響區(qū)距離在20～42 m之間，不會對下一周期冷井的運行產(chǎn)生影響。

圖7不同冷、熱負(fù)荷時的熱影響區(qū)半徑Fig.7Radius of Heat Affected Zone Under Different Cooling and Heating Loads

需要說明的是，如根據(jù)用戶負(fù)荷需求本身對系統(tǒng)冷、熱負(fù)荷工況進行調(diào)整，則這種調(diào)控措施完全能夠滿足實際工程運行的需求；如為了緩解熱貫通現(xiàn)象進行調(diào)整，則需要其他能源的運用或項目本身運營策略的調(diào)整作為補充，以滿足用戶對供暖制冷的需求。

3.3抽灌模式的影響

當(dāng)熱泵系統(tǒng)冷、熱負(fù)荷設(shè)計恒定時，可以通過改變系統(tǒng)循環(huán)水量和抽灌水溫差，分析熱泵抽灌模式調(diào)整對換熱性能的影響。下面分別對增大循環(huán)水量、減小抽灌水溫差以及減小循環(huán)水量、增大抽灌水溫差的運行工況進行模擬，并比較其結(jié)果。根據(jù)監(jiān)測數(shù)據(jù)可知，熱泵運行時實際的單井涌水量最高可達500～1 000 m3·d-1，全年平均為300 m3·d-1，說明該含水層開采能力和水泵功率有較大的可調(diào)整余量，但由于冷、熱負(fù)荷恒定，循環(huán)水流量與抽灌水溫差呈反比關(guān)系，溫差的限制使得流量不能過大或過小。如果實際運行中井1、井2供暖期回灌水在8 ℃上下變化，抽水溫度基本穩(wěn)定在20 ℃，抽灌水溫差在10～15 ℃之間。因此，循環(huán)水流量的調(diào)整范圍不宜超過±20%，各工況設(shè)定見表4。

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

圖8不同抽灌模式時的冷影響區(qū)半徑Fig.8Radius of Cold Affected Zone Under Different Pumping and Injection Modes

從模擬結(jié)果可知，冷、熱負(fù)荷設(shè)計恒定時各工況冷影響區(qū)范圍仍會到達熱井處，但在程度上存在差異，如圖8所示。工況6～9冷影響區(qū)擴散速度分別為12.0、13.2、9.8、9.1 m·a-1，預(yù)計發(fā)生熱貫通的時間分別為2 000、1 977、2 717、3 083 d。工況6、7冷影響區(qū)擴散速度分別比基礎(chǔ)工況提高了11.1%和22.2%，而工況8、9則分別降低了9.3%和15.7%?？梢缘玫窖h(huán)水量對冷影響區(qū)擴散速度的靈敏度為：每增加1%循環(huán)水量，冷影響區(qū)擴散速度提高1.11%；每減小1%循環(huán)水量，冷影響區(qū)擴散速度提高0.86%。由此可見，抽灌水流量增大、溫差減小時，熱貫通現(xiàn)象加?。怀楣嗨髁繙p小、溫差增大時，有利于緩解熱貫通趨勢。另外，減小抽灌流量還可以減小冷羽場地外方向的擴散距離，第9個供暖期末工況9的冷羽半徑比基礎(chǔ)工況減小了3.6%。因此，采用“大溫差、小流量”抽灌模式能減小系統(tǒng)運行對地溫場的影響。

隨著熱泵運行，各工況制冷期末的熱影響區(qū)距離均在20～40 m之間，呈現(xiàn)先擴散后收縮的趨勢(圖9)，平均收縮速率為1.2 m·a-1。熱影響區(qū)的發(fā)展與抽灌模式的選擇無明顯關(guān)系，對系統(tǒng)運行的影響也較小。值得注意的是，在冷、熱負(fù)荷恒定時，無論循環(huán)水量和抽灌水溫差如何調(diào)整，經(jīng)換熱器換熱后提供給溫室的熱(冷)量同樣是恒定的，所以調(diào)整水量和溫差后仍能符合溫室對溫度的要求。

圖9不同抽灌模式時的熱影響區(qū)半徑Fig.9Radius of Heat Affected Zone Under Different Pumping and Injection Modes

4結(jié)語

(1)采用COMSOL Multiphysics建立的熱滲耦合數(shù)值模型能較為準(zhǔn)確地模擬地下水源熱泵運行地溫場的演變特征?，F(xiàn)行工況的模擬結(jié)果顯示，熱泵系統(tǒng)的短期運行效果良好，但由于全年熱負(fù)荷大于冷負(fù)荷，隨著熱泵運行，冷影響區(qū)擴散范圍將逐漸增大，冷、熱源井之間會在第7個供暖期末出現(xiàn)熱貫通現(xiàn)象，長期運行將不利于熱井的取熱。

(2)減小熱負(fù)荷、增加冷負(fù)荷可以顯著延緩熱貫通發(fā)生的趨勢。若減小20%熱負(fù)荷的同時增加20%冷負(fù)荷，冷影響區(qū)的擴散速度比原來降低46.3%，在模擬周期內(nèi)并未發(fā)生熱貫通。雖然冷負(fù)荷增加使得熱影響區(qū)擴散距離增加，但不會對冷井的運行造成影響。因此，合理調(diào)節(jié)冷、熱負(fù)荷有利于熱泵系統(tǒng)的長期穩(wěn)定運行。

(3)當(dāng)冷、熱負(fù)荷設(shè)計恒定時，分別增大10%和20%循環(huán)水量，冷影響區(qū)擴散速度將分別比原工況提高11.1%和22.2%；如果分別減小10%和20%循環(huán)水量，則分別降低了9.3%和15.7%。說明減小抽灌水流量、增大溫差的抽灌模式能有效緩解熱貫通現(xiàn)象的發(fā)生。