深地能源工程熱水力多場(chǎng)耦合效應(yīng)高效模擬方法

趙志宏，劉桂宏，徐浩然

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

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

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

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

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

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

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

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

1 模擬方法

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

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

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

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

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

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

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

1.1 儲(chǔ)層

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

式中：t/s 為時(shí)間；fρ/(kg·m-3)為流體密度；φ為孔隙率；Qm/(kg·m-3·s-1)為流體質(zhì)量源項(xiàng)；達(dá)西流速u(mài)f為：

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

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

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

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

儲(chǔ)層中熱對(duì)流-傳導(dǎo)過(guò)程可用能量守恒方程描述：

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

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

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

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

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

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

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

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

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

1.2 蓋層

對(duì)于蓋層，其傳熱過(guò)程只考慮熱傳導(dǎo)：

1.3 斷層

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

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

斷層中熱對(duì)流-傳導(dǎo)過(guò)程可用能量守恒方程來(lái)描述：

1.4 一維井筒簡(jiǎn)化模型

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

式中：Aw/m2為地?zé)峋慕孛娣e；/(m/s)為沿井筒軸向的平均流速；Qw為通過(guò)地?zé)峋诎l(fā)生流體與外部圍巖的熱交換；fD為達(dá)西摩擦因子，與雷諾數(shù)(Re)、表面粗糙度e、地?zé)峋畯絛i有關(guān)：

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

對(duì)于湍流(Re＞2000)，Df可用Haaland 公式計(jì)算[32]： 定義等效傳熱系數(shù)(hZ)eq來(lái)近似描述地?zé)峋畠?nèi)流體與周?chē)鷰r石之間的傳熱過(guò)程：

式中：Tf/K 為地?zé)峋辛黧w的溫度；Ts/K 為巖石溫度?？紤]單位長(zhǎng)度內(nèi)從流體通過(guò)井筒的熱流相等，可推導(dǎo)出等效傳熱系數(shù)為：

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

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

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

式中，Pr表示普朗特?cái)?shù)。

1.5 初始條件與邊界條件

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

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

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

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

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

式中，lo/m 為地?zé)峋穆阊鄱伍L(zhǎng)度。

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

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

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

1.6 流體參數(shù)

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

1.7 模型求解過(guò)程

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

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

2 案例研究

2.1 北京東南城區(qū)地?zé)崽?/h3>

2.1.1 地?zé)岬刭|(zhì)條件

北京東南城區(qū)地?zé)崽镂挥谔彀查T(mén)廣場(chǎng)以東至東四環(huán)附近(圖4)，屬典型坳陷盆地，熱儲(chǔ)層包括薊縣系鐵嶺組(埋深約578 m～1200 m)和霧迷山組(埋深約517 m～2456 m)，均為水熱型熱儲(chǔ)，巖性為白云巖。由于長(zhǎng)期地質(zhì)構(gòu)造作用，上述兩熱儲(chǔ)層由洪水莊組頁(yè)巖隔開(kāi)，鐵嶺組熱儲(chǔ)由下馬嶺組頁(yè)巖覆蓋，受崇文門(mén)—胡家溝大斷裂帶的影響，地?zé)崽镏胁考捌狈较蜩F嶺組缺失(圖5)。

圖4 地?zé)峋植技捌拭鍵-I’位置示意圖 Fig.4 Locations of geothermal wells and geologic section I-I’

20 世紀(jì)70 年代北京東南城區(qū)地?zé)崽镩_(kāi)采以區(qū)域供暖為主，之后隨著開(kāi)采井?dāng)?shù)量的增加，地?zé)崴_(kāi)采量急劇增加，導(dǎo)致熱儲(chǔ)壓力迅速下降(水位下降0.83 m/a～2.36 m/a)，嚴(yán)重威脅地?zé)豳Y源的可持續(xù)開(kāi)采。自2000 年以來(lái)，通過(guò)地?zé)嵛菜毓嗖糠謱?shí)現(xiàn)了地?zé)豳Y源的可持續(xù)開(kāi)采。

2.1.2 模型建立與校正

數(shù)值模型尺寸為長(zhǎng)11.6 km，寬8.5 km，高2 km，共包含2 個(gè)熱儲(chǔ)層、9 個(gè)不同地層、4 條斷層以及39 口地?zé)峋?其中8 口回灌井)(圖6)。兼顧計(jì)算精度與效率，地?zé)峋c斷層周邊區(qū)域網(wǎng)格約尺寸約為2.5 m，而遠(yuǎn)離地?zé)峋驍鄬拥膮^(qū)域則使用較大尺寸的網(wǎng)格，約為 430 m。模型總計(jì)包含662817 個(gè)四面體單元，39 口地?zé)峋?291 個(gè)一維線單元代表(圖6)。

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

圖6 北京東南城區(qū)地?zé)崽飻?shù)值模型 Fig.6 Numerical model of the Beijing City geothermal field

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

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

2.1.3 計(jì)算結(jié)果分析

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

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

本文考慮兩種工況：一種為開(kāi)采率和回灌率相等，即地?zé)崽镒⒉杀葹?.26；另一種工況回灌率為 開(kāi)采率的2.875 倍，即地?zé)崽镒⒉杀葹?.0。通過(guò)對(duì)比分析這兩種模擬工況，揭示注采比對(duì)熱儲(chǔ)長(zhǎng)期性能的影響。注采比為1.0 是地?zé)豳Y源可持續(xù)開(kāi)發(fā)的趨勢(shì)，但我國(guó)實(shí)際地?zé)衢_(kāi)采中注采比仍遠(yuǎn)小于1.0。該地?zé)崽锘毓嗑恢梅植疾痪鶆?，選擇兩個(gè)包含采灌井的典型區(qū)域來(lái)解釋模擬結(jié)果(圖4中虛線方框)。在模擬過(guò)程中對(duì)開(kāi)采溫度進(jìn)行監(jiān)測(cè)，并將開(kāi)采溫度降低1 ℃的時(shí)間定義為開(kāi)采井的熱突破時(shí)間。兩種工況在普通臺(tái)式計(jì)算機(jī)(CPU i7-4790K，內(nèi)存16 GB)上都需要約9.5 h 的計(jì)算時(shí)間，這對(duì)于地?zé)犴?xiàng)目的設(shè)計(jì)和管理是合理的。

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

(續(xù)表)

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

圖10 北京東南城區(qū)地?zé)崽飻?shù)值模型初始 溫度場(chǎng)分布 /(℃) Figure 10 Initial temperature field in the numerical model of the Beijing City geothermal field

1) 工況I：注采比0.26

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

2) 工況II：注采比1.0

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

3) 熱儲(chǔ)變形分析

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

圖11 北京東南城區(qū)地?zé)崽飻?shù)值模型熱儲(chǔ)的溫度場(chǎng)分布 /(℃) Fig.11 Temperature distribution in the geothermal reservoirs of the Beijing City geothermal field

圖12 北京東南城區(qū)地?zé)崽飻?shù)值模型典型地?zé)峋臒嵬黄魄€及溫度分布 /(℃) 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ū)地?zé)崽飻?shù)值模型熱儲(chǔ)位移場(chǎng)分布 /mm Fig.13 Displacement distribution in the geothermal reservoirs of the Beijing City geothermal field

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

2.2 山東德州城區(qū)地?zé)崽?/h3>

2.2.1 地?zé)岬刭|(zhì)條件

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

圖15 德城區(qū)地?zé)峋恢梅植紙D Fig.15 Distribution of geothermal wells in the Decheng district

2.2.2 模型建立與校正

數(shù)值模型區(qū)域面積約為310 km2，深度2.5 km，共包含1 個(gè)熱儲(chǔ)層、5 個(gè)不同地層以及86 口地?zé)峋?圖16)。兼顧計(jì)算精度與效率，地?zé)峋爸車(chē)鷧^(qū)域的網(wǎng)格細(xì)化，最大單元尺寸2.5 m，熱儲(chǔ)層區(qū)域網(wǎng)格細(xì)化，最大單元尺寸1 km，其余蓋層網(wǎng)格粗化，最大單元尺寸12 km。模型總計(jì)包含308163 個(gè)四面體單元，86 口地?zé)峋?216 個(gè)一維線單元代表(圖16)。

圖16 德州城區(qū)地?zé)崽飻?shù)值模型 Fig.16 Numerical model of the Decheng district geothermal field

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

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

圖17 德城區(qū)年開(kāi)采量及地?zé)峋當(dāng)?shù)量 Fig.17 Evolution of geothermal well number and production volume in Decheng district

2.2.3 計(jì)算結(jié)果分析

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

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

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

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

3 討論

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

3.1 地?zé)崛壕?yīng)

在北京東南城區(qū)地?zé)崽镏校瑢?duì)比兩種模擬工況可以發(fā)現(xiàn)，在回灌率較大和長(zhǎng)時(shí)間的運(yùn)行條件下，更多的開(kāi)采井可能受到回灌的影響(圖12)，例如，回灌井G 周?chē)娜诘責(zé)峋?2、井15 和井9 組成了地?zé)崛壕到y(tǒng)，在德州城區(qū)地?zé)崽镏?，DZ48井、DZ54 井及其回灌井同樣組成了地?zé)崛壕到y(tǒng)，這意味著傳統(tǒng)的地?zé)釋?duì)井?dāng)?shù)值模型可能無(wú)法準(zhǔn)確地預(yù)測(cè)開(kāi)采井的熱突破，而本文提出的深地能源工程熱水力多場(chǎng)耦合高效數(shù)值模擬方法可為地?zé)嵯到y(tǒng)的長(zhǎng)期管理提供一種有效考慮地?zé)崛壕?yīng)的模擬 方案。

3.2 熱儲(chǔ)之間的相互作用

北京東南城區(qū)地?zé)崽锬Ｐ椭邪藘蓚€(gè)熱儲(chǔ)層，即霧迷山組和鐵嶺組熱儲(chǔ)層，本文重點(diǎn)研究了回灌井R1 和開(kāi)采井30 之間的相互作用，它們彼此接近但在不同的儲(chǔ)層中(圖23)，回灌井R1 位于深層的霧迷山組熱儲(chǔ)層中，而開(kāi)采井30 則位于淺層的鐵嶺組熱儲(chǔ)中，由于這兩個(gè)熱儲(chǔ)層被洪水莊組蓋層所 隔開(kāi)，R1 中的回灌冷水不會(huì)影響開(kāi)采井30 的溫度，由于蓋層不透水，冷鋒面只能通過(guò)熱傳導(dǎo)在蓋層中移動(dòng)，在地?zé)崛壕到y(tǒng)運(yùn)行50 a 后，霧迷山組熱儲(chǔ)的回灌不會(huì)影響鐵嶺組熱儲(chǔ)溫度。

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

圖21 德城區(qū)地?zé)崽锊晒嗑季桨?Fig.21 Well layout scheme for production and injection wells in the Decheng district geothermal field

圖22 德城區(qū)地?zé)崽锏湫偷責(zé)峋疅嵬黄魄€ Fig.22 Thermal breakthrough curves for the representative geothermal wells in the Decheng district geothermal field

表4 德城區(qū)地?zé)崽餆嵬黄茣r(shí)間統(tǒng)計(jì)表 Table 4 Thermal breakthrough times for the geothermal wells in the Decheng district geothermal field

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

3.3 斷層的影響

北京東南城區(qū)地?zé)崽锬Ｐ椭羞€考慮了幾條斷層對(duì)地?zé)崛壕到y(tǒng)運(yùn)行的影響(圖24～圖25)，例如，回灌井T3 和T4 的深度基本相同，但井T4 穿過(guò)了 斷層，回灌冷水可以快速流過(guò)斷層，并使靠近斷層區(qū)域的溫度降低。對(duì)于穿過(guò)斷層的開(kāi)采井22，由于深層地?zé)崴ㄟ^(guò)斷層流向井22，使得該井的開(kāi)采溫度隨時(shí)間緩慢增加(圖25)。

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

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

4 結(jié)論

本文提出了深地能源工程熱-水-力多場(chǎng)耦合高效模擬方法，并應(yīng)用于北京市東南城區(qū)地?zé)崽锖蜕綎|德城區(qū)地?zé)崽锏壬顚拥責(zé)衢_(kāi)采工程，研究了地?zé)嵯到y(tǒng)的群井效應(yīng)、采灌方案優(yōu)化設(shè)計(jì)、以及深部熱儲(chǔ)溫度場(chǎng)、滲流場(chǎng)、變形場(chǎng)的長(zhǎng)期演化特征。主要結(jié)論如下：

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

(2) 基于北京市東南城區(qū)地?zé)崽?，?shù)值模擬結(jié)果表明，地?zé)峋奈恢?、深度和采灌量?huì)對(duì)熱儲(chǔ)的長(zhǎng)期性能產(chǎn)生重大影響，且不可忽略地?zé)崛壕?yīng)，應(yīng)在實(shí)際工程中優(yōu)化地?zé)峋荚O(shè)位置以及開(kāi)采量與回灌量，以避免開(kāi)采井快速發(fā)生熱突破。如果熱儲(chǔ)層之間有蓋層隔開(kāi)，則不同熱儲(chǔ)層之間并無(wú)顯著影響。斷層可為流體流動(dòng)和熱量傳遞提供快速通道。熱儲(chǔ)變形同時(shí)受到回灌壓力和熱應(yīng)力的作用，這兩種應(yīng)力的主次作用決定了熱儲(chǔ)變形量的大小。

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