地下群井改變徑流條件下的三維實驗模擬研究

王煜曦

（北京市地熱研究院，北京 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）在地下徑流速度一定的情況下，加大抽水井抽水量地埋管周邊溫度恢復并不明顯，但地溫下降的變化率增大。恢復到原始地溫需增加地溫場的恢復時間。加大抽水量即加大了強化傳熱過程，溫度變化率加大，但在抽水后期熱傳導成為溫度恢復的主要因素，換熱管和周地層溫度梯度減小，導熱變慢，溫度場恢復所需時間變長。