張進平，王煜曦，劉桂宏，何鐵柱，袁利娟，徐浩然，孔祥軍

(1. 北京市地熱調查研究所，北京 102218；2. 清華大學土木工程系，北京 100084)

熱儲工程數值模擬技術已在地熱資源評估、開采回灌方案優(yōu)化設計以及地熱井開采壽命預測等方面發(fā)揮重要作用。在地熱開發(fā)早期，由于已知數據較少，采用高度理想化、概念清晰、計算量相對較小的解析模型[1]或集中參數模型[2]對熱儲層特性演化特征進行粗略預估。隨著計算機技術的發(fā)展，數值模擬可以精細刻畫熱儲層非均質性及其熱儲物理、化學、熱力學等在開采回灌過程中的演化規(guī)律，數值模擬已成為現在熱儲工程中多場耦合效應研究與評價的主要方法[3]。熱儲工程熱-水-力-化(THMC)多場耦合模型主要有兩類：連續(xù)介質力學模型和離散裂隙網絡模型[4-6]?；谶B續(xù)介質力學模型的熱儲工程數值模擬方法發(fā)展很快，其中以有限元方法應用最為廣泛。由于地熱井(直徑為0.1 m～0.2 m)和地熱田(面積大于100 km2)之間的尺度相差巨大，群井系統(tǒng)數值模擬通常無法在日常工作計算機上完成，而城市尺度地熱田群井系統(tǒng)高效模擬方法解決了這一難題[7]。目前熱儲模型大多是對井系統(tǒng)，李騰風等[8]考慮熱源作用下非飽和土體水熱耦合作用機制，基于格子Boltzmann方法，采用雙分布函數分別描述溫度場及水分場的演化過程，建立了相應的水熱耦合模型。徐浩然等[9]提出了一種模擬工程尺度碳酸鹽巖熱儲酸化壓裂過程的數值方法，考慮碳酸鹽巖熱儲酸化壓裂過程中熱、水、力、化四場之間的耦合作用。地熱田含有數十口以上的地熱井，故考慮群井效應對地熱田可持續(xù)開發(fā)有著重要意義。通州區(qū)為北京市行政副中心，其地熱勘探數據較為完整，為了確保其地熱資源優(yōu)勢在本區(qū)內得到最大程度的發(fā)揮，構建科學高效的地熱資源開發(fā)利用體系，文章結合了實際地層和斷裂帶構建了通州區(qū)地質模型，開展深層熱儲數值模擬，探尋北京市通州區(qū)深部地熱資源的優(yōu)化開采模式，為通州區(qū)規(guī)?；責豳Y源開發(fā)利用提供依據。

1 地質背景及熱儲模型建立

1.1 地質背景和地熱井分布

野外實地調查及資料顯示[10-11]，通州全區(qū)地勢平坦而略有起伏，主要地貌類型為沖積平原，由永定河、潮白河、溫榆河等沖洪積而成，可劃分為河流階地、河床及漫灘、砂質決口扇及砂丘、洼地、采土和砂石坑等5 種類型。

工作區(qū)基巖地質構造如圖1，工作區(qū)新生界下伏地層從老到新包括中上元古界、古生界和中生界地層，主要發(fā)育夏墊斷裂帶、西集斷裂、姚辛莊斷裂等多條斷裂帶，其中夏墊斷裂帶對工作區(qū)鉆井影響最大, 其余斷裂規(guī)模相對較小，對地層的錯動、切割等影響較小，對地層蓋層條件、深部儲層埋深、儲層溫度等影響相對較小。

圖1 北京市通州區(qū)基巖地質構造圖Fig. 1 Geological structure of bedrock in Beijing Tongzhou District

區(qū)內主要地層從新到老依次為新生界、石炭-二疊系、寒武-奧陶系、青白口系和薊縣系。新生界地層巖性以砂質粘土、粉細砂、泥巖、頁巖為主，總厚度巨大，超過2000 m。新生界之下為古生界石炭-二疊系、奧陶系、寒武系，主要巖性為泥巖、頁巖、砂巖、灰?guī)r夾頁巖。古生界之下為元古界地層，包括青白口系和薊縣系等；青白口系巖性以泥灰?guī)r、砂巖、頁巖為主；薊縣系包括鐵嶺組、洪水莊組、霧迷山組和楊莊組，鐵嶺組、霧迷山組和楊莊組巖性都以白云巖為主，洪水莊組巖性為頁巖。根據地層巖性，工作區(qū)內主要以薊縣系霧迷山組熱儲為主。

根據以上資料建立通州區(qū)三維地質模型，區(qū)域面積 906 km2，深度4.5 km，地層分為 7 層，包含20 條斷裂帶，如圖2(a)。地層由新到老依次為新生界(Kz)、石炭—二疊系(C-P)、寒武—奧陶系(?-O)、青白口(Qn)、薊縣系鐵嶺組和洪水莊組(Jxt+Jxh)、薊縣系霧迷山組(Jxw)和太古界(Ar)如圖2(b)。

圖2 通州區(qū)三維地質模型示意圖Fig. 2 Three-dimensional geological model of Tongzhou District

新生界全區(qū)范圍內均有分布，通州西南部的地層厚度超過了4500 m。石炭—二疊系只在通州南部小部分區(qū)域存在，其余地區(qū)缺失。寒武—奧陶系在通州區(qū)的東部和南部分布較完整，北部大部分地層缺失，只在副中心附近小區(qū)域存在。青白口系、薊縣系鐵嶺組和洪水莊組在工作區(qū)中部及西南部有不同程度的缺失，根據之前的鉆井資料，洪水莊組厚度約70 m 左右，為了簡化模型，將鐵嶺組和洪水莊組合并成一層。薊縣系霧迷山組為通州區(qū)主要開采的熱儲層，在全區(qū)范圍內均有分布，在南部及東部有部分地區(qū)缺失。太古界為基巖，全區(qū)范圍內均有分布，由于鉆井并未揭穿霧迷山組，東南部地區(qū)的太古界地層埋深可能超過4500 m，地層厚度分布及模型物理力學參數如表1 所示。

表1 地層分布及物理力學參數Table 1 Formation distribution and mechanical parameters

通州鉆鑿完成地熱井 25 眼，其中開采井 22 眼，回灌井3 眼，主要集中在副中心區(qū)域，南部及西部處于未開發(fā)狀態(tài)，沒有地熱井，地熱井分布如圖3。地熱井成井深度、出水溫度和出水量差異較大，現有數據統(tǒng)計：成井深度約642 m～3589 m，出水溫度約28 ℃～92 ℃，出水量約603 m3/d～3073 m3/d。

圖3 地熱井分布圖Fig. 3 Distribution of geothermal wells

1.2 多尺度模型的建立

COMSOL Multiphysics 是一款大型多物理場仿真軟件，以有限元法為基礎，通過求解偏微分方程、偏微分方程組來解決單物理場、多物理場問題，實現真實物理系統(tǒng)的模擬。多物理場問題本質上就是對偏微分方程組的求解，因此只要可以用偏微分方程組進行描述的物理現象，都可以用COMSOL Multiphysics 多物理場耦合來進行模擬分析。強大的多物理場耦合功能是該軟件相比于其他有限元軟件的主要優(yōu)勢所在。

1.2.1 地熱井模型簡化

李馨馨等[12]利用線單元模擬地熱對井提出了裂隙巖體三維熱流耦合的等效模擬方法，李敬元等[13]將地層巖石視為多孔介質，考慮了孔隙壓力和滲流作用的影響，根據彈性理論和摩爾-庫侖屈服條件以及巖石報傷后的軟化特性，對井筒周圍巖石進行了彈塑性力學分析，給出了井壁穩(wěn)定條件。應用COMSOL 軟件建立遠場尺度熱儲三維數值模型，對于細長的地熱井結構，其長度與半徑之比可達104以上，是建模中挑戰(zhàn)性的數值問題，Al-Khoury 等[14-15]和Saeid 等[16]提出了一種用一維線單元模擬地熱井中傳熱過程的偽三維井筒模型。

井筒結構示意圖如圖4(a)所示，蓋層以上部分放入套管，并采用水泥砂漿完井。本文進一步擴展了此模型，引入一維線單元對該三維井筒結構進行簡化[17-18]，圖4(a)為地熱井筒結構示意圖，圖4(b)為簡化的一維井筒模型，考慮沿井筒軸向的滲流傳熱過程，井筒內流體與蓋層的熱交換通過等效換熱系數來近似考慮，儲層部分放入濾水管，提供流體交換通道，其它主要假設條件包括：

圖4 井筒結構和一維地熱井模型Fig. 4 Well structure and one-dimensional model

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

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

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

4)考慮井筒內流體性質如密度、粘度、熱傳導系數、比熱容等與溫度的相關性，但不考慮井筒內氣液相變過程。

1.2.2 模型網格劃分與時間離散

兼顧計算精度與效率，地熱井及周圍區(qū)域的網格細化，最大單元尺寸2.5 m；熱儲層區(qū)域網格細化，最大單元尺寸1 km；其余蓋層網格粗化，最大單元尺寸6 km，如圖5 所示。模型總計包含299 959 個四面體單元，25 眼地熱井由798 個一維線單元代表。

圖5 模型網格剖分示意圖Fig. 5 Mesh generation diagram of the model

地熱系統(tǒng)只在每年的11 月15 日至來年的3 月15 日共4 個月的供暖季里處于開采或回灌運行狀態(tài)，在非供暖季的8 個月里地熱井處于停運狀態(tài)，圖6 為地熱井運行時間周期。

圖6 時間周期示意圖Fig. 6 Time period diagram

1.2.3 模型求解及邊界條件

將三維地熱井簡化成了一維線單元，考慮井筒內流體沿井筒軸向的滲流傳熱，井筒內流體與周圍巖體的換熱過程采用等效換熱系數近似考慮，該系數考慮了井徑、套管和水泥砂漿的導熱性、流體性質和流量的影響。蓋層傳熱過程可用熱傳導方程：

模型初始條件和邊界條件：熱儲初始溫度60 ℃，井口流速0.09 m/s，回灌溫度20 ℃，儲層初始水壓為靜水壓力，井筒內初始水壓也按靜水壓力考慮，運行開始后(t＞ 0)，井口設為流速邊界條件，井筒內初始溫度等于圍巖初始溫度。

模型的初始溫度場可根據地溫梯度 ΔT來確定，頂面溫度邊界條件設為全年平均氣溫，側面溫度開邊界和水位邊界，底邊溫度開邊界和不透水邊界，如圖7 所示。模型所用熱儲巖石熱物理參數如表1 所示。

圖7 模型邊界條件示意圖Fig. 7 Model boundary conditions

2 模型驗證及當前模式分析

2.1 計算模型與參數

根據不同地層的增溫率，Surfer 軟件插值得到不同地層的地溫梯度分布，霧迷山組下伏巖層地溫梯度數據未掌握，根據前人研究基礎，取1.5 ℃/100 m。選擇6 眼溫度監(jiān)測井進行擬合，調整模型的熱物性參數使得溫度監(jiān)測井中的模擬溫度與出水溫度保持一致。結果表示溫度與監(jiān)測溫度吻合結果較好，圖8 為溫度監(jiān)測井和溫度擬合數據，計算得到模型的初始溫度場分布如圖9 所示。

圖8 溫度擬合曲線Fig. 8 The temperature fitting curves

圖9 初始溫度場分布 /(℃)Fig. 9 The initial temperature distribution

利用區(qū)內5 眼地熱井的抽水試驗數據進行水位擬合，用裘布依公式和吉哈爾特降壓影響半徑經驗公式迭代計算出監(jiān)測井附近滲透系數的參考值。模型取大落程的滲透系數作為參考值，用Surfer軟件插值得到儲層的滲透系數分布，基于抽水試驗(5 眼監(jiān)測井)和回灌試驗(京通-4)的抽水量、降深數據與水位監(jiān)測數據進一步調整熱儲的滲透系數分布(表2)，最終使模擬水位與實測水位較好吻合，如圖10 所示。

表2 校正前后滲透系數對比表 /(m/d)Table 2 Comparison of permeability coefficient before and after calibration

圖10 抽水、回灌試驗水位擬合曲線Fig. 10 Fitting curves of water level of exploitation and recharge tests

2.2 真實采灌模式評價

區(qū)域內共25 眼地熱井，3 眼回灌井(通熱灌-9、京通灌-4、京通灌-6)。結合地熱井現有開采能力以及北京冬季供暖需求，假設回灌水溫25 ℃，當單井采灌量為1000 m3/d、1500 m3/d 和2000 m3/d 時，計算開采100 年后熱儲的溫度場和水位場的演化。

地熱回灌可能導致熱突破，將開采溫度下降2 ℃定義為熱突破，關注回灌井分布密集區(qū)域溫度場的變化。距離通熱灌-9 最近的井是通熱-15，井間距為182 m，冷鋒面很快到達了通熱-15，當采灌量分別為1000 m3/d、1500 m3/d、2000 m3/d，對應的熱突破時間為8 a、5 a、4 a，開采井中溫度逐漸下降并穩(wěn)定在38 ℃左右，如圖11 所示。距離京通灌-6 最近的井是京通-3，在采灌量較大的情況下(2000 m3/d)出現了熱突破，熱突破時間為82 a。京通灌-4 位置相對孤立，距離最近的京通-5 的井間距超過2000 m，不受到回灌影響。

圖11 開采溫度變化曲線Fig. 11 The temperature variation curves

回灌井分布較稀疏的區(qū)域，回灌不足可能會帶來水位場的變化。在長期開采的情況下，開采井集中區(qū)，形成大面積的抽水漏斗，并隨著采灌量的增加而增大。

即使回灌井稀少，但對維持熱儲壓力、穩(wěn)定開采水位起到了相當大的作用。隨著回灌量的增大、采灌井間距的減小，其作用也越明顯。但是不合理的回灌井布置，可能會導致開采井過早發(fā)生熱突破(如通熱灌-9 導致通熱-15 過早發(fā)生了熱突破)，因此，合理的采灌井布置才是實現研究區(qū)地熱資源長期開發(fā)利用的有效途徑。

3 優(yōu)化開采模型研究

針對目前采灌方案的不足，在100%回灌的前提下，結合工作區(qū)內的開采井、回灌井提出三種布井方案， “一采一灌”、“一采兩灌”和“兩采三灌”。

結合區(qū)域內地熱井的開采能力及薊縣系霧迷山組熱儲開采經驗，單井的采灌量分別考慮1000 m3/d、1500 m3/d 和2000 m3/d，采灌井間距分別考慮300 m、500 m 和800 m 三種工況，回灌溫度為供暖尾水溫度25 ℃?，F有采灌模式下水位場分布如圖12 所示。

圖12 現有采灌模式下水位場分布 /mFig. 12 Water level distribution with current exploitation and recharging mode

“一采一灌”布井方案中，采灌量大于1500 m3/d時，冷鋒面運移到了開采井處，開采井在78 a 發(fā)生了熱突破，采灌量增大為2000 m3/d，熱突破時間縮短至56 a；采灌井間距為500 m 和800 m 時，無論采灌量為1000 m3/d、1500 m3/d、2000 m3/d，回灌冷鋒面均未運移到開采井處(圖13)，開采井溫度也沒有下降，當采灌井間距大于500 m 時，開采井的熱突破風險小?！耙徊蓛晒唷辈季桨钢校捎诨毓嗑贾迷陂_采井下游位置，在三種采灌量下，回灌冷鋒面均未運移到開采井處，未造成開采井發(fā)生熱突破，開采井的熱突破風險最小(圖13)。“兩采三灌”的布井方案中，位于下游的開采井受到回灌冷鋒面的影響較大，熱突破風險高。

圖13 熱儲溫度場分布 /(℃)Fig. 13 Temperature distribution of heat storage

在三種采灌量下，下游的開采井均發(fā)生了熱突破，熱突破時間分別為96 a、72 a 和58 a；位于上游的開采井在采灌量小于1500 m3/d 沒有發(fā)生熱突破，當采灌量增大到2000 m3/d 時在85 a 發(fā)生了熱突破，位于下游的開采井熱突破風險極高，位于上游的開采井也面臨著一定程度的熱突破風險。

從水位計算結果來看：不同布井方案下，熱儲的水位場分布有所不同，在開采井附近會形成低壓區(qū)，回灌井附近形成高壓區(qū)(圖14)。在三種布井方案下，隨著采灌量和井間距的增大，開采井的水位逐漸下降，回灌井的水位逐漸上升。采灌量的大小直接影響了采灌井水位的高低，而采灌井間距越小，回灌對穩(wěn)定開采水位的作用越明顯。

圖14 熱儲水位場分布 /mFig. 14 Water level distribution of heat storage

在“一采一灌”布井方案中，當采灌量增大到2000 m3/d 時，回灌井的水位會超過0 m， “一采一灌”布井方案無法實現100%自然回灌。而在“一采兩灌”和“兩采三灌”布井方案中，由于回灌井數量的增加而減小了單眼回灌井的回灌量，使得回灌井水位均在0 m 以下，可實現100%自然回灌。

優(yōu)化開采模式研究表明區(qū)內回灌井布設需要滿足以下條件：

① 采灌井間距大于500 m，采灌量不超過2000 m3/d；

② 大地流場方向由西北流向東南，回灌井需布置在開采井下游；

③ “一采兩灌”布井方案，部分開采井集中的地區(qū)可采用“兩采三灌”；

④ 回灌井的布設位置不應對周圍其他地熱井造成干擾。

根據以上原則區(qū)內已有的距離其他井較遠的開采井匹配回灌井時，采用“一采兩灌”的布設方案(若不新增開采井，無法實現“兩采三灌”)；部分距離較近的開采井，如京通-1 和通熱-19、通熱-9 和通熱-15、京通-3 和京通6、園-1 和園-2 等采用“兩采三灌”的布設方案，以節(jié)約運維成本。具體的回灌井布設方案如圖15 所示，區(qū)內共布設回灌井46 眼。在該布井方案下，集中開采區(qū)的抽水漏斗范圍得到了有效的控制，穩(wěn)定了開采水位，在100 年的開采年限里，開采溫度沒有明顯的下降，說明開采井未發(fā)生熱突破(圖16)。對比目前采灌模式下監(jiān)測井的水位可知，優(yōu)化采灌模式下能有效恢復并穩(wěn)定開采水位不超過150 m(圖17)，說明優(yōu)化采灌模式能滿足研究區(qū)內地熱資源的長期開發(fā)利用。

圖15 回灌井布設方案Fig. 15 The distribution plan of recharge wells

圖16 優(yōu)化采灌模式模擬結果示意圖Fig. 16 The simulation results of optimized exploitation and recharging mode

圖17 兩種模式下監(jiān)測井水位對比曲線Fig. 17 The curves of water levels monitored under the two modes

“一采兩灌”布井方案，回灌井布置在下游可大幅降低開采井的熱突破風險，且能滿足100%自然回灌的需求，通過對采灌井間距的敏感性分析可知，最適合研究區(qū)的采灌井間距應大于500 m，通過對采灌量的敏感性分析可知，為確保開采水位不超過150 m，研究區(qū)的采灌量不宜超過2000 m3/d。優(yōu)化采灌模式下能有效恢復并穩(wěn)定開采水位不超過150 m，優(yōu)化采灌模式能滿足研究區(qū)內地熱資源的長期開發(fā)利用。

4 結論

本文基于城市尺度地熱群井高效數值模擬方法，建立了工作區(qū)三維熱儲數值模型，并通過對溫度和水位監(jiān)測數據的擬合，校正了模型參數，為后面地熱資源優(yōu)化開采模式研究奠定了基礎。主要結論如下：

(1)通過對目前采灌模式的模擬研究，回灌對維持熱儲壓力和穩(wěn)定開采水位十分重要，并隨著回灌量的增大以及采灌井間距的減小，作用越明顯，不合理的回灌井布置，可能會導致開采井過早發(fā)生熱突破，因此，合理的采灌井布置才是實現工作區(qū)地熱資源長期開發(fā)利用的有效途徑。

(2)通過研究不同布井方案下的井組模型得到了最適合工作區(qū)的布井方案為“一采兩灌”和“兩采三灌”，且當采灌井間距大于500 m，采灌量不超過2000 m3/d 時，能保證100%回灌且在100 年開采年限里，開采井不發(fā)生熱突破，開采水位不超過150 m 的需求。提出了工作區(qū)內的優(yōu)化采灌布井方案原則，在滿足布設原則的優(yōu)化采灌模式下，開采井不會發(fā)生熱突破，且能有效恢復并穩(wěn)定開采水位不超過150 m，實現了區(qū)內地熱資源的長期開發(fā)利用。