地下水源熱泵抽灌井優化布置及參數靈敏度

摘要：為了確定地下水源熱泵的最佳布井方式，探明含水層溫度場變化特征以及水文地質參數靈敏性，以石家莊地下水源熱泵適宜區肖家營—東兆通為研究區，基于熱儲層的地熱地質條件和熱物性參數，利用FEFLOW 軟件構建了三維水熱耦合數值模擬模型.在２０００m３/d的抽灌量下對不同井間距、不同布井模式含水層溫度場變化進行模擬，探討發生熱突破的可能性，并探討溫度場關于熱導率、滲透系數、孔隙率等參數的靈敏度.研究結果表明：１）在單抽雙灌模式下，方案C（兩眼回灌井連線與地下水流向垂直，一眼回灌井位于抽水井正下游，另一眼回灌井位于沿地下水流向４５°方向，即兩回灌井與抽水井組成等腰直角三角形，抽水井正下游的回灌井為直角頂點）的布井方案對含水層溫度場影響最小，為最優方案;２）單抽雙灌模式下，抽灌量為２０００m３/d時，方案C抽灌井間距設置在４０～５０m 較為合理;３）抽灌過程中，溫度場對于滲透系數的改變靈敏度較高，而對孔隙率以及熱導率的改變靈敏度較低.

關鍵詞：地下水源熱泵;抽灌井優化;溫度場變化;參數靈敏度;石家莊

doi：１０．１３２７８/j．cnki．jjuese．２０２１０２２２ 中圖分類號：TP３１ 文獻標志碼：A

０ 引言

地下水源熱泵作為一種新型的淺層地熱能開發應用技術，以其綠色、環保、高效、節能的特性被世界各國應用及推廣.其通過向熱泵機組輸入少量電能，以高比熱容的地下水為媒介進行抽灌換熱，進而達到將分散的、不可直接利用的地熱能轉化成為可直接利用熱能的目的[１２].該技術的推廣應用對于環境治理和改善具有重要意義，但是在系統的使用過程中也會產生許多問題，例如：系統建設缺乏規劃可能會導致后期抽灌量不足、管井堵塞、機組受腐蝕等一系列問題，從而降低運行效率[３５];抽灌井布局不合理可能會導致回灌冷水的冷鋒面向開采井靠攏，從而發生熱突破的現象，進而降低系統的換熱功效，污染地下水環境.其中熱突破相關問題尚處于相對缺乏系統研究的局面.所以，為提高機組的運行效率，避免引發次生地質災害，有必要進行地下水源熱泵系統運行后對地溫場產生影響的研究.國內外對于水源熱泵以及熱突破問題有了一定的基礎研究.Park等[６]研究了水源熱泵系統中的熱擴散現象，提出地下水含水層熱擴散對地下溫度分布有很大的影響.Smajil等[７]將水源熱泵系統放入新的能源系統優化模型組，介紹了新型的集成形式，提高了地下水源熱泵的應用效率.Park等[８]將地下水流動與傳熱模擬模型與遺傳算法聯系起來作為優化技術，該模型可以共同或分開確定最優井位和抽注量.文韜等[９]在二維層面上研究了布井形式和井間距對于熱貫通的影響.肖銳等[１０]指出在地下水源熱泵設計中，在抽灌水井連線間布設熱屏障井可改變地下水流場，降低熱量在抽灌井間的運移速度，有利于延長熱貫通發生時間并緩解熱貫通程度，但是熱貫通問題沒有從根本上得到解決.

雖然國內外相關研究人員在熱突破問題的研究方面已經積累了豐富的經驗，掌握了一定的理論基礎;但是從已有的研究中可以看出，地下水源熱泵研究主要集中在地面以上應用工程以及小規模區域的熱貫通模擬研究中，缺乏深層的理論研究.其中，國內水源熱泵工程多集中在對單井和雙井循環系統的研究方面，而在較大區域內地下水源熱泵采能對溫度場的影響，以及含水層儲能的熱突破問題方面的理論研究都較少.本文在充分考慮石家莊市城區的地質、水文地質條件、地下水特征和巖石熱物性參數等因素的基礎上，建立了地下水源熱泵適宜區水熱耦合模擬模型，研究熱突破問題與抽灌水量、井間距以及井布局方案的關系，并對相關參數進行了靈敏度評價，以期提出較為合理的抽灌模式，對揭示地下水源熱泵熱運移規律，以及為石家莊市地下水源熱泵系統技術的推廣和可持續開發應用提供參考依據.

１ 研究區概況

研究區（圖１）位于石家莊肖家營—東兆通一帶，靠近滹沱河沖洪積扇中部地帶，面積為１９．７４km２，屬于地下水源熱泵適宜區[１１１２].區內屬北溫帶半濕潤、半干旱大陸性季風氣候，平均氣溫約１３．４９ ℃（圖２），１９５６—２０１５ 年平均日照時長數約為２４４０．４２h，全年無霜期１８０～２１０d，受自然因素與人為影響，年際氣溫呈逐年升高、日照時長呈降低趨勢[１３].地層主要由基底地層與第四系松散沉積層組成，基底為白堊系和第三系，呈臺階狀展布;第四系沉積物的成因主要是河流沖積、洪積，顆粒上粗下細，沉積厚度由山前向東部逐漸變厚，含水層巖性以砂層和卵礫石層為主.

２ 模型構建

２．１ 水文地質概念模型

研究區含水層巖性以砂、卵礫石層為主，厚度１１５m，其中含水層頂板設為７０ m，底板設置為－４５m，水力坡降為０．１２０％.根據研究區水文地質條件構建水文地質概念模型，將該模型分為９個亞層，其中，第一亞層為粗砂層（７０～６５m），第二亞層為粉砂、中砂互層（６５～５０m），第三至九亞層為砂、卵礫石（５０～－４５m）.地下水流向為自西北向東南.為測試單抽雙灌模式的最佳布井方式，將抽回灌井濾管放置在埋深為６０～１１０m 的位置.根據石家莊地區實際的供暖周期，將每年的１１月１５日至翌年３月１５日設置為供暖期，其余時間為間歇期，機組２４h 運轉，根據研究區單位涌水量為５０～１００m３/（hm），抽灌井水量在此區間內較為適宜，同時選取一系列梯度的抽灌量模擬進行試算，得出抽灌量為２０００m３/d具有較好的代表性.將回灌井的初始溫度設置為８℃.整個地下水流動系統可概化為非均質各向異性三維非穩定流系統.

２．２ 數值模型

研究區的水文地質概念模型可以概化為水平結構、非均質、各向異性的潛水三維非穩定流系統.地下水運動的基本微分方程[１４１５]為：

２．３ 模型初始條件和參數設定

２．３．１ 初始條件

１）滲流場的邊界條件

考慮到研究區水力坡降很小，可以近似忽略其對區域流場的影響，將四周設置為定水頭邊界，故將模型西北側邊界設置為定水位邊界，水頭值為２６．４５m，將東南側邊界設置為定水位邊界，水頭值為２５．２５m.將垂直于地下水流向的方向設置為隔水邊界.

２）溫度場的邊界條件

將模型的頂板溫度設置為變溫邊界即大氣溫度，溫度變化根據近年平均大氣氣溫變化給出，以時間曲線的形式進行賦值.模型底部邊界作已知溫度邊界處理，設為定溫邊界，已有研究[１７]計算得出研究區恒溫層的深度在３０～３５m，平均深度為３２m，恒溫層的溫度為１４．４０００ ℃，平均地溫梯度為２．８１℃/hm.根據恒溫層地溫梯度、深度以及溫度計算，可將底板溫度設置為１６．７３２３℃.

３）初始溫度場

頂板溫度為平均年內大氣溫度，根據研究區概況可以得到埋深在３０～３５m 時所得恒溫層溫度為１４．４０００℃，將底板溫度設置為１６．７３２３℃，內插值法得出其他深度的溫度值，見圖３.

２．３．２ 參數設定與模型校正

研究區的各地層的熱物性參數根據含水層巖性給出.表１為含水層水文地質參數及巖土熱物性參數值[１８１９].

在不設置地下水源熱泵抽灌井的情況下，對研究區水熱耦合模型進行模擬.經運算，模型運行１０a后，溫度場基本達到穩態，以穩態后的溫度場作為初始溫度場.為觀測不同深度溫度場變化規律，在模型中心不同深度同一節點處設置觀測井，設置深度分別為０、５、２０、３０、４０m.模型繼續運行１０a后，不同深度溫度場變化曲線見圖４.

從圖４結合圖２可以看出：觀測點溫度變化趨勢與大氣溫度基本一致，呈周期性變化，這是由于受大氣熱交換的影響;地表處觀測點溫度變幅最大，與大氣溫度變化保持一致，隨著深度的增加，溫度變幅逐漸變小，深度達到４０m 時，溫度基本不變;地表以下觀測點溫度與大氣溫度存在滯后期，且滯后性隨著深度的增加愈發嚴重，這是由于受熱傳導的影響.此外，３０ m 與４０ m 觀測點溫度變幅較小，３０m深度溫度在１３．３～１４．６ ℃變化，４０m 深度溫度在１４．３～１４．４℃之間變化，基本不受大氣溫度影響，該深度與恒溫層深度３２m 相近;說明對于模型的概化以及含水層水文地質參數與巖土熱物性參數的賦值基本符合石家莊市地下巖土分布規律.

２．４ 網格剖分

模型采用三角網格剖分模式，并在抽水井和回灌井處進行網格加密，共剖分成２９００個網格.為便于觀測不同深度的溫度變化，模型垂向上共剖分為９層結點網絡面，分為８個網絡層，層面高程自上而下分別為７０、６５、５０、４０、３０、２０、１０、－４０、－４５m.

２．５ 方案設定

研究區抽灌模式設置為單抽雙灌模式，單抽雙灌模式下井的間距設置為３０，４０，５０m;井的抽水量設置為２０００m３/d;結合實際模擬時產生的溫度場情況將熱突破現象產生的溫度差限定為０．１０００℃，即如果抽水井處的溫度降低超過０．１０００ ℃，則認為出現熱突破[２０].傳輸模式以熱對流為主時，水力梯度的大小和方向為影響熱傳輸的主要因素;當水力梯度大小不變時，為探究單抽雙灌模式下不同井間距、不同布井方式對于熱突破的影響，將水力梯度方向對熱突破現象的抑制作用最大化，分別以平行于地下水流向和垂直于地下水流向布置抽水井和回灌井.同時考慮到回灌井的占地面積和影響半徑，以及對抽水井溫度的影響等現實問題，以回灌井布置在抽水井的正下游為原則，在確定一眼回灌井位置后，分別在抽水井對角線和垂向方向上布置另一眼回灌井.共設置４種布井方案（圖５）.

方案A：兩眼回灌井位于抽水井的正下游，連線與地下水流向平行，抽水井與回灌井位于同一條流線上（圖５a）.

方案B：兩眼回灌井的連線與地下水流向垂直，抽水井與回灌井位于同一條等水頭線上（圖５b）.

方案C：兩眼回灌井連線與地下水流向垂直，一眼回灌井位于抽水井正下游，另一眼回灌井位于沿地下水流向４５°方向，即兩回灌井與抽水井組成等腰直角三角形，抽水井正下游的回灌井為直角頂點（圖５c）.

方案D：一眼回灌井位于抽水井正下游，另一眼回灌井位于抽水井的垂向方向，兩回灌井與抽水井組成等腰直角三角形，抽水井為直角頂點（圖５d）.

３ 結果與討論

３．１ 抽灌井布井方式優化

為確定最佳布井方式和其對溫度場的影響，文中采取兩種對比形式：橫向對比（定流量、定井間距下，不同排布方式對天然溫度場的影響）;縱向對比（同一種排列方式、不同井間距對天然溫度場的影響）.在定流量、不同井間距下，通過試算得知，比較符合現實狀況的井間距為４０m，因此，選擇井間距為４０m、流量為２０００m３/d，根據圖５中A、B、C、D４種不同布井方案進行模擬.從中選取對于天然溫度場影響最小的布井方案，溫度場變化結果見圖６.

由圖６可見，A、B、C、D４種布井方案在２０００m３/d、４０m 井間距下均產生熱突破現象，其中方案B、方案D布井的方式熱突破現象更為明顯.

其中：方案A 為兩口回灌井交互影響，對于天然溫度場影響較大并產生疊加現象，使抽水井處的溫度產生較強波動，產生了一定的熱突破現象（圖６a）;方案B為兩回灌井連線垂直于地下水流向，對流產生的熱交換較為明顯，熱突破現象較方案A 更為明顯（圖６b）;方案C 為兩回灌井位于抽水井下游，其中一回灌井沿地下水流向４５°方向，因側方向的回灌井對于抽水井影響較小，產生熱突破的現象不明顯（圖６c）;方案D產生了較明顯的熱突破現象（圖６d）.通過設置的溫度監測孔提取了每種方案抽水井周圍溫度的變化，結果見圖７和表２.可知：方案B和方案D溫度差值分別為０．７５７７，０．７２７５ ℃（表２），溫度曲線波動非常明顯（圖７），二者均產生了較明顯的熱突破現象;方案A 和方案C溫度差值分別為０．３０７２，０．１８３８℃（表２），溫度曲線波動較平緩（圖７），溫度差較小，但都大于０．１０００ ℃，也產生了熱突破.方案A產生的熱突破值大于方案C，這是因為方案A 中二回灌井互相干擾較大，增大了水動力場的擾動，使溫度變化加劇，進而導致抽水井處溫度變化大.通過４類布井方案的比較可知，單抽雙灌模式下井的間距為４０m 時，方案C布井最為合適.

為了縱向對比探究同一種排列方式下不同井間距對天然溫度場的影響，選擇方案C 布井，保持抽灌水量為２０００m３/d，分別對３０m 和５０m 的間距也進行了模擬.在３０m 井間距下，溫度曲線下凹明顯（圖８），溫度差值達到０．９４１４℃（表３），在水動力場的作用下整體呈現２個相連的橢球形，與抽水井產生強烈的干擾，可見發生了強烈的熱突破現象（圖９a），說明在回灌井自身的干擾下，增大了低溫的影響范圍，同時由于井間距的減小，回灌井的冷流大股匯入抽水井，對抽水井工程影響較大，同時熱交換效率降低，工程壽命縮短.在５０m 間距下，溫度曲線較為平直（圖８），溫度場則呈現２個涇渭分明的橢球，水動力場的作用不足以將回灌井的冷流帶入５０m 外的抽水井處（圖９b），５０m 井間距下回灌井溫度變化僅僅為０．０５３１℃（表３），尚未發生熱突破現象.但抽灌井間距增大勢必增大占地面積，亦會增加工程成本;因此綜合考慮，在２０００m３/d的抽灌量下，采用方案C 布井，井間距設置在４０～５０m較為合理.

３．２ 參數靈敏度分析

為了探究不同參數對于溫度場的靈敏度，選擇方案C布井模式，在４０m 井間距、２０００m３/d流量的情況下，探究滲透系數、孔隙率、熱導率三者對溫度場的影響.本文采取將各項數值分別提高至１５０％（方案一、三、五）和減少至８０％（方案二、四、六）的方式來探究參數靈敏度大小.具體的模擬方案如表４所示，抽水井溫度變化曲線見圖１０.

３．２．１ 孔隙率

選用表４中方案一與方案二對孔隙率參數重新進行賦值，將孔隙率分別提高至１５０％或減小至８０％，根據熱突破溫度差的大小得到孔隙率對于溫度場的影響強度.由圖１０a和表５中可以發現：孔隙率改變對熱突破溫度差的影響不大，當孔隙率增大至１５０％（方案一）時，熱突破溫度差為０．１６１０℃;而孔隙率降低至８０％（方案二）時，熱突破溫度差為０．１９７０℃，兩種方案熱突破溫度差變幅僅為０．０３６０℃.可見溫度場對于孔隙率的變化靈敏度較低.

３．２．２ 滲透系數

選用表４中方案三與方案四對滲透系數參數重新進行賦值，通過滲流系數變化前后產生的熱突破溫度差大小來分析溫度場對于滲透系數的靈敏度.由圖１０b和表６可知：方案三滲透系數增大至１５０％時，回灌井處的冷水未對抽水井處溫度造成影響，熱突破溫度差僅為０．０３４１ ℃，未發生熱突破現象;而方案四滲透系數減小至８０％，熱突破溫度差為０．４５８７ ℃，發生了熱突破現象，兩種方案熱突破溫度差變幅高達０．４２４６℃.綜上，滲透系數的變化對溫度場的影響較大.

３．２．３ 熱導率

選用表４中方案五與方案六對熱導率參數重新進行賦值，通過熱導率變化前后產生的熱突破溫度差大小來分析溫度場對于滲透系數的靈敏度.由圖１０c和表７可知：熱導率的改變對溫度場的影響變化不大，當熱導率增大至１５０％（方案五）時，熱突破溫度差為０．１９７１℃;而熱導率降低至８０％（方案六）時，熱突破溫度差為０．１７９４ ℃，兩種方案熱突破溫度差變幅僅為０．０１７７ ℃.可見熱導率對于由對流傳熱方式控制的溫度場所能造成的影響較小，當以對流為主要傳熱方式時，溫度場對巖石傳熱效率的變化靈敏度較低.

４ 結論

１）對于單抽雙灌模式下，抽灌井距離４０m 以及２０００m３/d的抽灌量下采用方案C（兩眼回灌井連線與地下水流向垂直，一眼回灌井位于抽水井正下游，另一眼回灌井位于沿地下水流向４５°方向，即兩回灌井與抽水井組成等腰直角三角形，且抽水井正下游的回灌井為直角頂點）布井方式，對于含水層溫度場影響最小，其熱突破溫度差為０．１８３８℃.

２）采用方案C布井，在單抽雙灌模式下，抽灌水量為２０００m３/d時，３０，４０，５０m 井間距下熱突破溫度差分別為０．９４１４，０．１８３８，０．０５３１℃，故將抽灌井間距控制在４０～５０m 之間較為合理.３）在抽灌井運轉，溫度以對流方式傳輸為主時，溫度場對于滲透系數的改變靈敏度較高，將滲透系數提高至１５０％時，熱突破溫度差為０．０３４１ ℃，減小至８０％時，熱突破溫度差為０．４５８７ ℃;而溫度場對孔隙率以及熱導率的改變敏感度較小，溫度變幅均在０．０５００℃之內.

參考文獻（References）：

[１] ?地熱能開發利用“十三五”規劃?發布[J]．地質裝備，２０１７，１８（２）：３ ４．

[２] 孫明行，王瑞湖，管彥武，等．廣西陸域干熱型地熱資源潛力估算[J]．吉林大學學報（地球科學版），２０２２，５２（４）：１３０２ １３１３．

[３] 王海龍．北京平原區地下水源熱泵數值模擬研究[D]．北京：中國地質大學（北京），２００９．

[４] 馮波，崔振鵬，趙璞，等．基于T２WELL的U 型地熱井供暖潛力數值模擬[J]．吉林大學學報（地球科學版），２０２２，５２（２）：５６０ ５７０．

[５] 趙振，秦光雄，羅銀飛，等．西寧盆地地熱水特征及回灌結垢風險[J]．水文地質工程地質，２０２１，４８（５）：１９３ ２０４．

[６] Park B H，Bae G O，Lee K K．Importance of ThermalDispersityinDesigningGroundwaterHeat Pump（GWHP）System：FieldandNumericalStudy[J]．RenewableEnergy，２０１５，８３：２７０ ２７９．

[７] SmajilH，OderskyL，HamacherT．Integrationof Groundwater Heat Pumps into Energy System Optimization Models [J]．Energy，２０２２，２３８：１２１６０７．

[８] ParkD，LeeE，KaownD，etal．Determinationof Optimal Well Locations and Pumping/Injection Ratesfor Groundwater Heat Pump System [J]．Geothermics，２０２１，９２：１０２０５０．

[９] 文韜，崔先澤，范勇．地下水源熱泵抽灌井群優化布置研究[J/OL]．長江科學院院報，https：//kns．cnki．net/kcms/detail/４２．１１７１．tv．２０２１０２２６．１６１８．０１２．html．

[１０] 肖銳，黃堅，王小清．熱屏障井對地下水源熱泵換熱影響模擬[J]．水文地質工程地質，２０２１，４８（２）：１９０１９８．

[１１] 閆佰忠，孫劍，陳佳琦，等．基于自適應BPNN GIS耦合的地下水源熱泵適宜性分區研究[J]．現代地質，２０２１．DOI：１０．１９６５７/j．geoscience．１０００ ８５２７．２０２１．０６３．

[１２] 閆佰忠，孫劍，王昕洲，等．基于GIS FAHP的石家莊市地下水源熱泵適宜性分區[J]．吉林大學學報（地球科學版），２０２１，５１（４）：１１７２ １１８１．

[１３] 湯麗莎，張健挺．日照時數減少的原因分析及應對措施[J]．氣象研究與應用，２００９，３０（增刊１）：１０８１０９．

[１４] 竇明，張彥，米慶彬，等．地溫空調井布局方式對地下水流場和溫度場的影響分析[J]．鄭州大學學報（工學版），２０１４，３５（５）：１２４ １２８．

[１５] 馬子涵，邢會林，靳國棟，等．基于微地震數據的增強型地熱儲層參數及采熱的數值模擬研究[J]．水文地質工程地質，２０２２，４９（６）：１９０ １９９．

[１６] 李海嶺，周雨澤，曹亞新，等．基于地下水壓采目標的地下水源熱泵系統井群調控方案[J]．南水北調與水利科技，２０２０，１８（５）：１６５ １７３．

[１７] 孫劍．石家莊市地下水源熱泵適宜區及含水層溫度場特征變化研究[D]．石家莊：河北地質大學，２０２０．

[１８] 倪龍，馬最良．含水層參數對同井回灌地下水源熱泵的影響[J]．天津大學學報，２００６，３９（２）：２９ ２３４．

[１９] 趙朝兵，程偉，田自浩，等．河北平原區第四系（２００m 以淺）巖土體綜合熱導率特征綜述[J]．河北地質，２０１６（１）：８ １０．

[２０] 陳響亮．抽灌井群熱交互性及其布控特性研究[D]．長春：吉林大學，２０１１．