趙志宏，劉桂宏，徐浩然

(清華大學土木工程系，北京 100084)

當今世界面臨的土地短缺、環(huán)境污染和能源枯竭等問題，使得深部地下空間開發(fā)利用已成必然趨勢[1―2]，如水電工程引水隧道最大埋深突破2400 m，南水北調(diào)輸水隧道最大埋深1150 m，高放廢物地質(zhì)處置庫埋深500 m～1000 m，地熱、礦產(chǎn)、油氣資源開采深度已達3000 m～8000 m(圖1)。深地能源工程旨在將地球深部蘊藏的能源提取出來或?qū)⑷祟惿a(chǎn)生活產(chǎn)生的廢物永久封存于地球深部，都涉及到儲層巖體復雜的熱(T)-水(H)-力(M)-化(C)多場耦合效應(yīng)及其動態(tài)調(diào)控。深部巖體在高地應(yīng)力、高地溫、高滲透壓等極端條件的耦合作用下，不僅巖體本身的物理力學性狀較淺部發(fā)生了顯著改變，而且多物理、化學場的耦合效應(yīng)更為明顯。深地能源工程導致地震、地下水位下降等災(zāi)害事故頻發(fā)，通常難以有效預(yù)測與防治，故全面深入理解深部巖體的多場耦合效應(yīng)是深地能源工程成功建設(shè)和安全運行的重要理論基礎(chǔ)。

圖1 典型深地能源工程示意圖 Fig.1 Schematic diagram for the typical deep geo-energy projects

由于深地能源工程的隱蔽性，僅通過室內(nèi)試驗和現(xiàn)場實測很難完全掌握儲層巖體復雜的耦合機理及其長期耦合效應(yīng)，主要存在的技術(shù)困難有：如何準確獲取深部巖體的物理力學參數(shù)；如何準確施加反映深部巖體賦存環(huán)境的邊界條件與初始條件；如何實現(xiàn)試驗數(shù)據(jù)的實時監(jiān)測與動態(tài)傳輸。準確高效的數(shù)值模擬則可以克服試驗研究遇到的上述難題，是研究深地能源工程多場耦合效應(yīng)的另一重要手段。自20 世紀80 年代初，國內(nèi)外學者已開始關(guān)注深部巖體多場耦合效應(yīng)，尤其在DECOVALEX (development of coupled models and their validation against experiments)[3]、 GTO-CCS (geothermal technology office’s code comparison study)[4]等國際合作項目的推動下，多孔介質(zhì)多場耦合作用的理論框架與模擬算法已日臻成熟，已可滿足近場模擬的需求。但是，對于遠場尺度且考慮圍巖-結(jié)構(gòu)相互作用的數(shù)值模擬，在精度和效率方面尚需提高。

已在深地能源工程多場耦合效應(yīng)評價中廣泛應(yīng)用的數(shù)值模擬方法主要包括有限元、邊界元、有限差分、有限體積等，主要的代表性多場耦合模擬程序如表1 所示[5―6]，其中，多數(shù)程序都可進行滲流傳熱過程的耦合分析，在此基礎(chǔ)上部分程序可進行熱-水-力三場、甚至熱-水-力-化四場的耦合模擬。Rutqvist[21]采用TOUGH-FLAC 研究了阿爾及利亞InSalah CO2地質(zhì)封存、美國Geysers 地熱田等深地項目中的巖體變形與壓力變化規(guī)律。朱萬成等[22]提出了巖體損傷過程中的熱-水-力多場耦合模型，并在COMSOL平臺實現(xiàn)了耦合方程的有限元求解。Kong 等[23]基于OpenGeoSys 建立了熱儲井間距優(yōu)化方法。THYME3D[24]、ROLG 等[25]程序考慮了液氣轉(zhuǎn)換與蒸氣傳輸。Pan 等[26]將細胞自動機的思想引入到巖石力學中，利用簡單的弱化元胞單元來模擬復雜裂隙網(wǎng)絡(luò)模型的力學行為，并應(yīng)用到熱-水-力耦合模擬中?？梢姡谶B續(xù)介質(zhì)力學的深部巖體多場耦合模擬方法發(fā)展很快。

表1 深地能源工程多場耦合程序[5―6] Table 1 The codes for modeling the coupled processes in deep geo-energy engineering[5―6]

多數(shù)深地能源工程都用到鉆井進行能量或物質(zhì)的轉(zhuǎn)移，但常規(guī)井筒結(jié)構(gòu)(直徑約0.1 m～0.2 m)與儲層本身(＞100 km2)存在3 個到4 個數(shù)量級以上的尺度差異，導致大尺度數(shù)值模型中精細刻畫井筒結(jié)構(gòu)的計算難度很大。Al-Khoury等[27―28]和Saeid等[29]將三維井筒結(jié)構(gòu)簡化為一維線單元模型，顯著提高了深地工程數(shù)值模擬的效率。

本文基于三維井筒結(jié)構(gòu)一維線單元假設(shè)，提出一種適合于深地能源工程遠場尺度的熱-水-力多場耦合效應(yīng)高效模擬方法，并依托我國京津冀及周邊地區(qū)典型水熱型地熱田群井系統(tǒng)開展案例研究，揭示深部熱儲溫度場、滲流場、變形場的長期演化 特征。

1 模擬方法

深地能源工程包含若干儲層與蓋層，以及大量開采井和回灌井，其典型井筒結(jié)構(gòu)如圖2 所示，通常采用多開成井模式，包括泵室段、井壁段、濾水段等，井段數(shù)量根據(jù)井深與儲層性質(zhì)設(shè)計。蓋層以上部分放入套管，并采用水泥砂漿完井。儲層部分放入濾水管，提供流體交換通道。

圖2 深地能源工程典型井筒結(jié)構(gòu) Fig.2 Schematic diagram of a typical well completion in deep geo-energy engineering

本節(jié)引入一維線單元對該三維井筒結(jié)構(gòu)進行簡化，考慮沿井筒軸向的滲流傳熱過程，井筒內(nèi)流體與蓋層的熱交換通過等效換熱系數(shù)來近似考慮，其它主要假設(shè)條件包括：

1) 儲層處于完全飽和狀態(tài)，且不考慮儲層內(nèi)氣液相變過程；

2) 蓋層處于完全干燥狀態(tài)，且不考慮不同儲層之間的水力聯(lián)系；

3) 井筒內(nèi)流體沿軸向流動，且同一深度流速沿徑向處處相等；

4) 考慮井筒內(nèi)流體性質(zhì)如密度、粘度、熱傳導系數(shù)、熱容等與溫度的相關(guān)性，但不考慮井筒內(nèi)氣液相變過程。

1.1 儲層

儲層中滲流過程可用質(zhì)量守恒方程描述：

式中：t/s 為時間；fρ/(kg·m-3)為流體密度；φ為孔隙率；Qm/(kg·m-3·s-1)為流體質(zhì)量源項；達西流速uf為：

式中：μ/(Pa·s)為流體動力粘度；κ/m2為儲層滲透率；pf/Pa 為流體壓力；g/(m·s-2)為重力加速度；z/m為豎向坐標，向上為正。

根據(jù)多孔介質(zhì)彈性理論有以下關(guān)系[30]：

式中：S/Pa-1為儲水系數(shù)，可表示為孔隙率φ、Biot系數(shù)Bα、流體體積模量Kf和儲層巖體體積模量Kd的函數(shù)：

聯(lián)立式(1)～式(3)可得：

儲層中熱對流-傳導過程可用能量守恒方程描述：

式中：Cp,f/(J·kg-1·K-1)為恒壓下流體的比熱容；T/K為溫度；Q/(W·m-3)為熱源； (ρCp)eff/(J·kg-1·K-1)為儲層巖體等效體積熱容，可表示為：

式中，effk/(W·m-1·K-1)為儲層巖體有效導熱系數(shù)，是儲層導熱系數(shù)sk和流體導熱系數(shù)fk的加權(quán)平均值：

多孔介質(zhì)彈性理論可用于描述熱儲中流體、溫度、變形之間的相互作用關(guān)系，即應(yīng)力張量σ、應(yīng)變張量ε、熱應(yīng)變Tε與孔隙水壓力fp的關(guān)系為[30]：

式中，C/(N/m)為儲層巖體剛度張量，應(yīng)在排水且恒定孔壓條件下測量應(yīng)力引起的應(yīng)變。

溫度應(yīng)變可以表示為：

式中：Tα為儲層巖體熱膨脹系數(shù)；refT/K 為儲層初始溫度。

根據(jù)Biot 定理，流體孔隙壓力與多孔基質(zhì)的膨脹和流體含量有如下關(guān)系：

式中，M/Pa 表示Biot 模量，是儲水系數(shù)S的倒數(shù)。 在純重力荷載下，處于平衡狀態(tài)的儲層可用Navier 方程來表示：

式中，ρav=ρfφ+ (1 -φ)ρs/(kg·m-3)為儲層巖體等效密度。

1.2 蓋層

對于蓋層，其傳熱過程只考慮熱傳導：

1.3 斷層

儲層中常含有斷層等不連續(xù)體，其滲流過程可用質(zhì)量守恒方程描述：

式中：dfr/m 為斷層開度；qfr為斷層中的體積流量，可表示為：

斷層中熱對流-傳導過程可用能量守恒方程來描述：

1.4 一維井筒簡化模型

一維井筒單元中不可壓縮流體的能量守恒方程為[31]：

式中：Aw/m2為地熱井的截面積；/(m/s)為沿井筒軸向的平均流速；Qw為通過地熱井壁發(fā)生流體與外部圍巖的熱交換；fD為達西摩擦因子，與雷諾數(shù)(Re)、表面粗糙度e、地熱井徑di有關(guān)：

對于層流(Re＜2000)，Df與井筒的表面粗糙度e無關(guān)，可用Stokes 公式計算：

對于湍流(Re＞2000)，Df可用Haaland 公式計算[32]： 定義等效傳熱系數(shù)(hZ)eq來近似描述地熱井內(nèi)流體與周圍巖石之間的傳熱過程：

式中：Tf/K 為地熱井中流體的溫度；Ts/K 為巖石溫度?？紤]單位長度內(nèi)從流體通過井筒的熱流相等，可推導出等效傳熱系數(shù)為：

式中：r0/m、r1/m 和r2/m 分別為套管內(nèi)徑、外徑和水泥砂漿外徑；kp/(W·m-1·K-1)和kc/(W·m-1·K-1)分別為套管和水泥砂漿的導熱系數(shù)。

流體井筒內(nèi)產(chǎn)生的熱膜阻hint可通過下式計算：

式中：Nu為努塞爾數(shù)。對于層流時，Nu為常數(shù)3.66。對于湍流(3000＜Re＜6×106，0.5＜Pr＜2000)，努塞爾數(shù)可由下式確定：

式中，Pr表示普朗特數(shù)。

1.5 初始條件與邊界條件

儲層初始水壓按靜水壓力考慮，Dirichlet 或Neumann 邊界條件均可應(yīng)用于模型頂部、側(cè)面和底部邊界。井筒內(nèi)初始水壓也按靜水壓力考慮，運行開始后(t＞ 0)，井口設(shè)為流速邊界條件。

模型初始溫度場可根據(jù)地溫梯度TΔ 確定：

式中：Ts/K 為地面溫度。Dirichlet 或Neumann 邊界條件均可應(yīng)用于熱儲模型頂部、側(cè)面和底部邊界。井同內(nèi)初始溫度等于圍巖初始溫度，運行開始后(t＞ 0)，回灌井口溫度固定為：

開采井井口的熱通量設(shè)為：

在井筒與儲層圍巖界面，將質(zhì)量通量(Mt)應(yīng)用于儲層內(nèi)邊界：

式中，lo/m 為地熱井的裸眼段長度。

根據(jù)回灌井井底溫度(Tb)，在儲層內(nèi)邊界設(shè)置線熱源：

開采井井底溫度取為儲層內(nèi)邊界溫度均值[29]：

模型頂部為自由邊界，底面設(shè)為固定位移邊界。模型側(cè)邊水平位移為0，只允許豎向位移。

1.6 流體參數(shù)

深部儲層中流體密度、粘度、導熱系數(shù)和比熱容隨溫度變化按如下經(jīng)驗公式考慮[33―34]：

1.7 模型求解過程

基于有限元計算平臺COMSOL，建立遠場尺度深地能源工程三維數(shù)值模型，采用四面體單元離散數(shù)值模型，一維線單元代表井筒。對于儲層巖體中的熱-水-力耦合控制方程(式(5)、式(6)、式(12))、蓋層巖體中傳熱控制方程(式(13))、井筒中滲流傳熱控制方程(式(16)和式(21))，采用COMSOL 中的并行直接稀疏求解器接口求解非線性系統(tǒng)，模型求解流程如圖3 所示。

圖3 計算流程示意圖 Fig.3 Calculation procedure for the reservoir-well model

2 案例研究

2.1 北京東南城區(qū)地熱田

2.1.1 地熱地質(zhì)條件

北京東南城區(qū)地熱田位于天安門廣場以東至東四環(huán)附近(圖4)，屬典型坳陷盆地，熱儲層包括薊縣系鐵嶺組(埋深約578 m～1200 m)和霧迷山組(埋深約517 m～2456 m)，均為水熱型熱儲，巖性為白云巖。由于長期地質(zhì)構(gòu)造作用，上述兩熱儲層由洪水莊組頁巖隔開，鐵嶺組熱儲由下馬嶺組頁巖覆蓋，受崇文門—胡家溝大斷裂帶的影響，地熱田中部及偏北方向鐵嶺組缺失(圖5)。

圖4 地熱井分布及剖面I-I’位置示意圖 Fig.4 Locations of geothermal wells and geologic section I-I’

20 世紀70 年代北京東南城區(qū)地熱田開采以區(qū)域供暖為主，之后隨著開采井數(shù)量的增加，地熱水開采量急劇增加，導致熱儲壓力迅速下降(水位下降0.83 m/a～2.36 m/a)，嚴重威脅地熱資源的可持續(xù)開采。自2000 年以來，通過地熱尾水回灌部分實現(xiàn)了地熱資源的可持續(xù)開采。

2.1.2 模型建立與校正

數(shù)值模型尺寸為長11.6 km，寬8.5 km，高2 km，共包含2 個熱儲層、9 個不同地層、4 條斷層以及39 口地熱井(其中8 口回灌井)(圖6)。兼顧計算精度與效率，地熱井與斷層周邊區(qū)域網(wǎng)格約尺寸約為2.5 m，而遠離地熱井或斷層的區(qū)域則使用較大尺寸的網(wǎng)格，約為 430 m。模型總計包含662817 個四面體單元，39 口地熱井由1291 個一維線單元代表(圖6)。

圖5 I-I’剖面示意圖 Fig.5 Geologic section I-I′ of geothermal field

圖6 北京東南城區(qū)地熱田數(shù)值模型 Fig.6 Numerical model of the Beijing City geothermal field

模型內(nèi)包含4 口水位監(jiān)測井，即井2、井3、 井5 和井7，監(jiān)測歷時11 年(1971 年―1981 年)，開采井流量與監(jiān)測井水位變化如圖7 所示[35]。通過對水位監(jiān)測數(shù)據(jù)的擬合來反演熱儲滲透率及標定模型邊界水位，所求得的熱儲滲透率及其余輸入?yún)?shù)如表2 所示。計算結(jié)果表明，4 口監(jiān)測井的模擬水位與實測水位基本吻合，可在此基礎(chǔ)上模擬兩個熱儲對回灌的響應(yīng)以及預(yù)測地熱井在未來50 年使用壽命內(nèi)的熱突破。此外，由于鐵嶺組熱儲的厚度小于霧迷山組熱儲(霧迷山組未揭穿)且整體水位在1971 年更低，隨著開采井數(shù)量及開采量的增加，經(jīng)過11 年的開采，到了1981 年，鐵嶺組熱儲出現(xiàn)了大范圍的抽水漏斗(圖8)。

東南城區(qū)地熱田溫度約為90 ℃，而回灌水溫度為30 ℃，故應(yīng)按式(31)～式(34)考慮流體性質(zhì)(密度、粘度、導熱系數(shù)和比熱容)隨溫度的變化。Yang 等[35]提供了地溫梯度分布(圖9)，結(jié)合模型輸入?yún)?shù)(表2)，計算得到模型的初始溫度場分布如圖10 所示。

2.1.3 計算結(jié)果分析

圖7 北京東南城區(qū)地熱田數(shù)值模型水位校正 Fig.7 Calibration of water level for numerical model of the Beijing City geothermal field

圖8 北京東南城區(qū)地熱田數(shù)值模型水位分布圖 /m Fig.8 Distribution of water level in the numerical model of the Beijing City geothermal field

本文考慮兩種工況：一種為開采率和回灌率相等，即地熱田注采比為0.26；另一種工況回灌率為 開采率的2.875 倍，即地熱田注采比為1.0。通過對比分析這兩種模擬工況，揭示注采比對熱儲長期性能的影響。注采比為1.0 是地熱資源可持續(xù)開發(fā)的趨勢，但我國實際地熱開采中注采比仍遠小于1.0。該地熱田回灌井位置分布不均勻，選擇兩個包含采灌井的典型區(qū)域來解釋模擬結(jié)果(圖4中虛線方框)。在模擬過程中對開采溫度進行監(jiān)測，并將開采溫度降低1 ℃的時間定義為開采井的熱突破時間。兩種工況在普通臺式計算機(CPU i7-4790K，內(nèi)存16 GB)上都需要約9.5 h 的計算時間，這對于地熱項目的設(shè)計和管理是合理的。

表2 模型擬合及輸入?yún)?shù)列表[36―38] Table 2 Fitting and input parameters for the numerical model[36―38]

(續(xù)表)

圖9 北京東南城區(qū)地熱田數(shù)值模型地溫 梯度分布 /(℃/100 m) Fig.9 Distribution of temperature gradient in the numerical model of the Beijing City geothermal field

圖10 北京東南城區(qū)地熱田數(shù)值模型初始 溫度場分布 /(℃) Figure 10 Initial temperature field in the numerical model of the Beijing City geothermal field

1) 工況I：注采比0.26

圖11 給出兩個熱儲溫度場隨時間的演變過程。隨著開采時間的延長，冷鋒面由回灌井向開采區(qū)移動，部分靠近回灌井的開采井發(fā)生了熱突破(圖12)，對于霧迷山組熱儲，由于井29 和井35 距離回灌井R2 只有240 m 和460 m，因此井29 在8.8 a 就發(fā)生了熱突破，井35 在22.5 a 發(fā)生了熱突破。對于鐵嶺組熱儲，井31 和井R3 形成地熱對井系統(tǒng)，井31在23.3 a 發(fā)生熱突破，而回灌井井G 的冷鋒面在 50 a 內(nèi)尚未到達井32、井15、井9 和井27。

2) 工況II：注采比1.0

與工況I 相比，由于增大了回灌量，冷鋒面在工況II 中運移得更快，導致有更多的開采井發(fā)生了熱突破(圖12)，對于霧迷山組熱儲，井29、井35和井5 分別在2.1 a、8.0 a 和14.8 a 時就發(fā)生了熱突破，但回灌沒有影響井16 的溫度變化，因為井16距離回灌井R2 和井E 約為883 m 和611 m。對于鐵嶺組熱儲，除了井27 以外的所有開采井都發(fā)生了熱突破，考慮到地熱系統(tǒng)會運行很長時間，井27在未來也有可能發(fā)生熱突破。

3) 熱儲變形分析

從回灌井中注入冷水到高溫高壓的熱儲中，必然會導致熱儲在回灌壓力及熱應(yīng)力的作用下發(fā)生變形，在長期回灌條件下，回灌井周圍的熱儲逐漸被冷卻，沉降量較大的區(qū)域均出現(xiàn)在回灌井周圍(圖13)，但這并不一定意味著回灌一開始就導致熱儲產(chǎn)生較大的變形，事實上，在回灌開始的較短時間內(nèi)，熱儲因回灌壓力的作用導致熱儲發(fā)生了膨脹，即產(chǎn)生了正向位移，并反映在地表沉降量變化中(圖14)，說明在回灌初期，回灌壓力對熱儲的變形起主導作 用。但隨著回灌時間的增加，熱儲逐漸被冷卻，這時熱應(yīng)力對熱儲的變形作用大于回灌壓力的作用，并發(fā)生了“熱收縮”現(xiàn)象，導致熱儲沉降量逐漸增大。對于開采井，由于開采高溫地熱水，導致開采井周圍的巖石遇熱發(fā)生了膨脹，在地表開采井周圍均出現(xiàn)了正向位移(5 mm～12 mm)，但隨著冷鋒面逐漸向開采井運移，開采井井底周圍巖石遇冷收縮，導致沉降量逐漸增大，如井29(圖14)。

圖11 北京東南城區(qū)地熱田數(shù)值模型熱儲的溫度場分布 /(℃) Fig.11 Temperature distribution in the geothermal reservoirs of the Beijing City geothermal field

圖12 北京東南城區(qū)地熱田數(shù)值模型典型地熱井的熱突破曲線及溫度分布 /(℃) Fig.12 Thermal breakthrough curves of the representative production wells and temperature distribution in the geothermal reservoirs of the Beijing City geothermal field

圖13 北京東南城區(qū)地熱田數(shù)值模型熱儲位移場分布 /mm Fig.13 Displacement distribution in the geothermal reservoirs of the Beijing City geothermal field

圖14 北京東南城區(qū)地熱田數(shù)值模型位移變化曲線 Fig.14 The displacement curve of the geothermal reservoirs and ground surface of the Beijing City geothermal field

2.2 山東德州城區(qū)地熱田

2.2.1 地熱地質(zhì)條件

研究區(qū)位于山東省德州市德城區(qū)，自中、新生代以來，受燕山運動和喜馬拉雅運動的影響，地殼運動總的趨勢是以下降為主，長期接受沉積，并沉積了巨厚的新生界地層，厚度超過3000 m，其下為中生界地層。鉆孔揭露的地層分布從上往下分別為：第四系平原組(Qp)、新近系明化鎮(zhèn)組(Nm)、新近系館陶組(Ng)、古近系東營組(Ed)和古近系沙河街組(Es)，熱儲層為新近系館陶組，區(qū)內(nèi)共有地熱井86 口，其中回灌井2 口，監(jiān)測井4 口(圖15)。

圖15 德城區(qū)地熱井位置分布圖 Fig.15 Distribution of geothermal wells in the Decheng district

2.2.2 模型建立與校正

數(shù)值模型區(qū)域面積約為310 km2，深度2.5 km，共包含1 個熱儲層、5 個不同地層以及86 口地熱井(圖16)。兼顧計算精度與效率，地熱井及周圍區(qū)域的網(wǎng)格細化，最大單元尺寸2.5 m，熱儲層區(qū)域網(wǎng)格細化，最大單元尺寸1 km，其余蓋層網(wǎng)格粗化，最大單元尺寸12 km。模型總計包含308163 個四面體單元，86 口地熱井由2216 個一維線單元代表(圖16)。

圖16 德州城區(qū)地熱田數(shù)值模型 Fig.16 Numerical model of the Decheng district geothermal field

研究區(qū)內(nèi)共有4 口水位監(jiān)測井，即DZ1 井、DZ28 井、DZ48 井和DZ53 井，監(jiān)測歷史20 年(1998年―2017 年)。從1998 年開始，陸續(xù)有新的地熱井投入使用且開采量逐年增加(圖17)，通過對水位監(jiān)測數(shù)據(jù)的擬合來反演熱儲滲透率及標定模型邊界水位(水位埋深約70 m)，所求得的熱儲滲透系數(shù)及其余輸入?yún)?shù)如表3 所示。計算結(jié)果表明，4 口監(jiān)測井的模擬水位與監(jiān)測水位基本吻合(圖18)，可在此基礎(chǔ)上利用模型進行德城區(qū)的采灌優(yōu)化設(shè)計，此外，隨著開采井數(shù)量的增加，研究區(qū)內(nèi)因抽水而形成的漏斗范圍在逐年擴大(圖19)。

研究區(qū)內(nèi)有1 口溫度監(jiān)測井，即DZ17 井，該井為回灌井，回灌期2016 年12 月14 日―2017 年4 月30 日，回灌結(jié)束后從7 月4 日開始，每隔1 個月監(jiān)測不同井深的溫度變化，持續(xù)4 個月，至11月3 日結(jié)束。通過調(diào)整模型的導熱系數(shù)和比熱容等參數(shù)，使模擬溫度與加測溫度基本吻合，最終得到的模型參數(shù)如表3 所示，計算結(jié)果表明，模擬溫度值與實測溫度值基本一致(圖20)，說明可用該模型進行后續(xù)計算。

圖17 德城區(qū)年開采量及地熱井數(shù)量 Fig.17 Evolution of geothermal well number and production volume in Decheng district

2.2.3 計算結(jié)果分析

為控制館陶組熱儲層地下水水位的持續(xù)下降和地熱尾水的排放污染，地熱供暖尾水必須100%回灌，以保證采灌均衡。研究區(qū)采用“一采一灌”對井模式，區(qū)內(nèi)目前DZ17 井與DZ1 井組成地熱對井系統(tǒng)，DZ31 井與DZ28 井組成地熱對井系統(tǒng)，需給其余82 口開采井匹配相應(yīng)的回灌井，在采灌量等于90 m3/h 的前提下，模擬出防止開采井熱突破的最優(yōu)采灌井間距，考慮到住宅小區(qū)的范圍與規(guī)模，采灌井間距初步設(shè)置方案為200 m、300 m、400 m 和500 m，布井方案如圖21 所示，回灌溫度為30 ℃，并考慮采灌周期(4 個月開采，8 個月停采)的影響。通過計算，得到了DZ48 井、DZ53 井和DZ56 井在不同采灌井間距下的熱突破曲線(圖22)，隨著采灌井間距的增大，開采井發(fā)生熱突破的時間越長，曲線的振蕩頻率與采灌周期有關(guān)，振幅與井深關(guān)，開采井井深越大，溫度曲線的振幅越大。將3 口監(jiān)測井在不同采灌井間距下的熱突破時間進行統(tǒng)計(表4)，當在井間距為200 m 時，3 口監(jiān)測井的熱突破時間均為13 a，熱突破時間較短，說明采灌井間距不應(yīng)小于200 m，DZ48 井在井間距小于500 m 時相比DZ53 井和DZ58 井來說，熱突破時間較短，這是由于DZ48 井受到DZ54 井回灌井的影響，造成DZ48 井的熱突破時間較短，由于DZ53井和DZ58 井位置比較孤立，不受周圍回灌井的影 位置，避免出現(xiàn)“一采多灌”的情況發(fā)生，建議在實際工程中，采灌井間距應(yīng)大于400 m 才能延長開采井的熱突破時間。

圖18 德城區(qū)地熱田數(shù)值模型水位校正 Fig.18 Calibration of water level for the numerical model of the Decheng district geothermal field

圖19 德城區(qū)地熱田數(shù)值模型不同年份的水位埋深分布云圖 /m Fig.19 Distribution of water levels in the numerical model of the Decheng district geothermal field

表3 德城區(qū)地熱田數(shù)值模型輸入?yún)?shù)列表 Table 3 Parameters for the numerical model of the Decheng district geothermal field

3 討論

響，這兩口監(jiān)測井的熱突破時間基本相同。以上結(jié)果說明，隨著采灌井間距增大，能有效延長熱突破時間，但同時也要考慮在集中開采區(qū)采灌井的相對

3.1 地熱群井效應(yīng)

在北京東南城區(qū)地熱田中，對比兩種模擬工況可以發(fā)現(xiàn)，在回灌率較大和長時間的運行條件下，更多的開采井可能受到回灌的影響(圖12)，例如，回灌井G 周圍的三口地熱井井32、井15 和井9 組成了地熱群井系統(tǒng)，在德州城區(qū)地熱田中，DZ48井、DZ54 井及其回灌井同樣組成了地熱群井系統(tǒng)，這意味著傳統(tǒng)的地熱對井數(shù)值模型可能無法準確地預(yù)測開采井的熱突破，而本文提出的深地能源工程熱水力多場耦合高效數(shù)值模擬方法可為地熱系統(tǒng)的長期管理提供一種有效考慮地熱群井效應(yīng)的模擬 方案。

3.2 熱儲之間的相互作用

北京東南城區(qū)地熱田模型中包括了兩個熱儲層，即霧迷山組和鐵嶺組熱儲層，本文重點研究了回灌井R1 和開采井30 之間的相互作用，它們彼此接近但在不同的儲層中(圖23)，回灌井R1 位于深層的霧迷山組熱儲層中，而開采井30 則位于淺層的鐵嶺組熱儲中，由于這兩個熱儲層被洪水莊組蓋層所 隔開，R1 中的回灌冷水不會影響開采井30 的溫度，由于蓋層不透水，冷鋒面只能通過熱傳導在蓋層中移動，在地熱群井系統(tǒng)運行50 a 后，霧迷山組熱儲的回灌不會影響鐵嶺組熱儲溫度。

圖20 DZ17 井溫度擬合曲線 Fig.20 Temperature fitting curve of the well DZ17 in the Decheng district geothermal field

圖21 德城區(qū)地熱田采灌井布井方案 Fig.21 Well layout scheme for production and injection wells in the Decheng district geothermal field

圖22 德城區(qū)地熱田典型地熱井熱突破曲線 Fig.22 Thermal breakthrough curves for the representative geothermal wells in the Decheng district geothermal field

表4 德城區(qū)地熱田熱突破時間統(tǒng)計表 Table 4 Thermal breakthrough times for the geothermal wells in the Decheng district geothermal field

圖23 北京東南城區(qū)地熱田鐵嶺組和霧迷山組熱儲的溫度分布圖 /(℃) Fig.23 Distribution of temperature in Tieling and Wumishan formation geothermal reservoirs

3.3 斷層的影響

北京東南城區(qū)地熱田模型中還考慮了幾條斷層對地熱群井系統(tǒng)運行的影響(圖24～圖25)，例如，回灌井T3 和T4 的深度基本相同，但井T4 穿過了 斷層，回灌冷水可以快速流過斷層，并使靠近斷層區(qū)域的溫度降低。對于穿過斷層的開采井22，由于深層地熱水通過斷層流向井22，使得該井的開采溫度隨時間緩慢增加(圖25)。

圖24 北京東南城區(qū)地熱田霧迷山組熱儲的溫度分布(包括T3 和T4) /(℃) Fig.24 Temperature distribution in the Wumishan reservoir including well T3 and T4

圖25 北京東南城區(qū)地熱田井22 的熱突破曲線 及熱儲溫度分布 /(℃) Fig.25 Thermal breakthrough curve for the well 22 and the temperature distribution in the Beijing City geothermal field

4 結(jié)論

本文提出了深地能源工程熱-水-力多場耦合高效模擬方法，并應(yīng)用于北京市東南城區(qū)地熱田和山東德城區(qū)地熱田等深層地熱開采工程，研究了地熱系統(tǒng)的群井效應(yīng)、采灌方案優(yōu)化設(shè)計、以及深部熱儲溫度場、滲流場、變形場的長期演化特征。主要結(jié)論如下：

(1) 針對深地能源工程井筒結(jié)構(gòu)獨特的細長型幾何特點，將井筒簡化為一維線單元，考慮井筒內(nèi)流體沿井筒軸向的滲流傳熱，而井筒內(nèi)流體與周圍巖體的換熱過程則采用等效換熱系數(shù)近似考慮，套管和砂漿層的影響包含在等效換熱系數(shù)中，同時，給出了描述熱儲熱-水-力多場耦合過程的表達式，并考慮了斷層的影響和流體性質(zhì)隨溫度的變化，利用該方法實現(xiàn)了遠場尺度深地能源工程熱-水-力多場耦合效應(yīng)的高效模擬。

(2) 基于北京市東南城區(qū)地熱田，數(shù)值模擬結(jié)果表明，地熱井的位置、深度和采灌量會對熱儲的長期性能產(chǎn)生重大影響，且不可忽略地熱群井效應(yīng)，應(yīng)在實際工程中優(yōu)化地熱井布設(shè)位置以及開采量與回灌量，以避免開采井快速發(fā)生熱突破。如果熱儲層之間有蓋層隔開，則不同熱儲層之間并無顯著影響。斷層可為流體流動和熱量傳遞提供快速通道。熱儲變形同時受到回灌壓力和熱應(yīng)力的作用，這兩種應(yīng)力的主次作用決定了熱儲變形量的大小。

(3) 基于山東省德州城區(qū)地熱田，根據(jù)“一采一灌”布井方案，在模型中給82 口開采井匹配了相對應(yīng)的回灌井，提出了一套布井方案，并對井間距進行了敏感性分析，計算結(jié)果表明，在實際工程中，采灌井間距大于400 m 才能避免開采井發(fā)生快速熱突破，且在集中開采區(qū)布置回灌井時應(yīng)考慮采灌井的相對位置，避免“一采多灌”的不利情況。