馬子涵，邢會(huì)林，靳國(guó)棟，譚玉陽(yáng)，閆偉超，李四海

（1.深海圈層與地球系統(tǒng)前沿科學(xué)中心，海底科學(xué)與探測(cè)技術(shù)教育部重點(diǎn)實(shí)驗(yàn)室，中國(guó)海洋大學(xué)海洋地球科學(xué)學(xué)院, 山東 青島 266100；2.青島海洋科學(xué)與技術(shù)試點(diǎn)國(guó)家實(shí)驗(yàn)室，山東 青島 266237；3.中國(guó)海洋大學(xué)海底科學(xué)與工程計(jì)算國(guó)際中心，山東 青島 266100）

地?zé)崮茏鳛橐环N清潔、穩(wěn)定、可再生的新型能源，符合綠色環(huán)?？沙掷m(xù)的發(fā)展要求，被譽(yù)為21世紀(jì)最具潛力的能源之一[1]。干熱巖型地?zé)豳Y源是指儲(chǔ)存在地下低滲透性高溫巖體中的熱能[2]，約占地?zé)豳Y源總量的30%，相當(dāng)于全球石油、天然氣和煤炭能量總和的30 倍[3]。但是，干熱巖儲(chǔ)層具有埋藏深、巖石致密、滲透率極低等特點(diǎn)，使用傳統(tǒng)技術(shù)難以有效開發(fā)其中的熱能[4]。利用增強(qiáng)型地?zé)嵯到y(tǒng)（Enhanced Geothermal Systems, EGS），即使用水力壓裂技術(shù)將原有儲(chǔ)層改造為高滲透性的人工儲(chǔ)層，而后通過(guò)水、二氧化碳等工作流體的循環(huán)注采，可實(shí)現(xiàn)干熱巖中熱能的高效開采[5]。EGS 技術(shù)作為實(shí)現(xiàn)深層干熱巖型地?zé)豳Y源綜合開發(fā)利用的關(guān)鍵技術(shù)，已經(jīng)成為各國(guó)新能源發(fā)展的重點(diǎn)關(guān)注方向之一[6]。

利用EGS 技術(shù)提取熱能是一個(gè)包含傳熱、滲流等多物理場(chǎng)耦合的復(fù)雜過(guò)程[7]，因其影響參數(shù)眾多[8]，物理性質(zhì)相差懸殊，時(shí)空演化不均衡等特點(diǎn)，通常利用數(shù)值模擬技術(shù)分析EGS 采熱過(guò)程[9-10]。近年來(lái)，國(guó)內(nèi)外學(xué)者提出了不同的數(shù)值模型，其中最常用的方法為等效多孔介質(zhì)模型。Jiang 等[6]、Ayling 等[11]建立了用于研究EGS 采熱的熱流耦合模型，將熱儲(chǔ)層等效為多孔介質(zhì)，進(jìn)行模擬并與實(shí)際情況進(jìn)行對(duì)比，驗(yàn)證了模型的有效性。Llanos 等[12]為庫(kù)珀盆地Habanero EGS 開發(fā)了一個(gè)三維熱流模型，考慮了滲透率各向異性，結(jié)果表明滲透率各向異性是控制儲(chǔ)層流體流動(dòng)的關(guān)鍵，同時(shí)為干熱巖儲(chǔ)層的開發(fā)設(shè)計(jì)了不同的井布局，并得出東西向交錯(cuò)的布局在該模型中為最優(yōu)井布局。Hogarth等[13]在此模型基礎(chǔ)上進(jìn)行了單井注入和2 口井聯(lián)通、閉環(huán)測(cè)試的模擬，將模型計(jì)算的生產(chǎn)井溫度曲線與實(shí)際記錄溫度進(jìn)行對(duì)比，驗(yàn)證了模擬結(jié)果的準(zhǔn)確性。Samin 等[14]提出孔隙度均一的三維多孔介質(zhì)THM 模型，儲(chǔ)層分為基巖和裂縫，裂縫區(qū)的滲透率具有各向異性，以此模型優(yōu)化裂縫區(qū)滲透率、流體注入壓力、注采井距等，實(shí)現(xiàn)熱發(fā)電性能與成本控制的綜合改善。然而，均勻多孔介質(zhì)模型沒(méi)有考慮儲(chǔ)層空間分布上的非均質(zhì)性，過(guò)于理想化，不能精確地表征壓裂改造后的儲(chǔ)層滲流特征[15]。此外，另一種比較常用的方法為規(guī)則裂縫或隨機(jī)裂縫模型，Sun 等[1]提出二維模型中儲(chǔ)層由巖石基質(zhì)和規(guī)則裂縫組成，并以此模型對(duì)地?zé)醿?chǔ)層中的流體流動(dòng)、傳熱等特征進(jìn)行了綜合模擬，通過(guò)敏感性分析研究了EGS 采熱的主控參數(shù)，結(jié)果表明THM 耦合效應(yīng)對(duì)采熱效率有重要影響。Xu等[2]提出由隨機(jī)裂縫和巖石基質(zhì)組成的模型，以此研究EGS 的熱流耦合系統(tǒng)，將該模型應(yīng)用于Habanero 熱儲(chǔ)層，對(duì)其采熱過(guò)程進(jìn)行模擬，應(yīng)用結(jié)果表明該模型是模擬工業(yè)規(guī)模地?zé)醿?chǔ)層傳熱的有效方法。Wang等[4]提出的THM 模型由巖石基質(zhì)塊體和隨機(jī)裂縫組成，模擬包含三維隨機(jī)裂縫的地?zé)醿?chǔ)層的長(zhǎng)期地?zé)嵘a(chǎn)過(guò)程，并將其應(yīng)用于澳大利亞的Habanero EGS 項(xiàng)目，評(píng)估了20 a 期間不同注入和生產(chǎn)壓力下的采熱量，結(jié)果表明，通過(guò)適當(dāng)控制注采壓差，可以實(shí)現(xiàn)持續(xù)穩(wěn)定的供電。另外有學(xué)者利用單裂縫概念化模型研究采熱過(guò)程，即認(rèn)為熱儲(chǔ)層是由基質(zhì)巖體和水平單裂縫組成或?qū)⒘芽p簡(jiǎn)化為一維線段嵌入巖石基質(zhì)網(wǎng)絡(luò)中，Yao 等[16]在三維熱流單裂縫模型基礎(chǔ)上，研究區(qū)域流動(dòng)方向?qū)Σ蔁岬挠绊?，結(jié)果表明，采熱可持續(xù)性隨著區(qū)域流動(dòng)方向和2 井中心連線的夾角的增加而增加。Li 等[17]在THM 耦合單裂縫模型基礎(chǔ)上進(jìn)行開發(fā)過(guò)程研究，發(fā)現(xiàn)高滲透率或大孔徑會(huì)導(dǎo)致過(guò)早發(fā)生熱突破，因此要在流動(dòng)和傳熱之間進(jìn)行權(quán)衡。然而，在這些模型中單一或者規(guī)則的裂縫分布過(guò)于理想，不能體現(xiàn)儲(chǔ)層中裂縫的真實(shí)情況；隨機(jī)裂縫相對(duì)靈活，但自由變量過(guò)多、模型數(shù)據(jù)量大、試驗(yàn)數(shù)據(jù)缺乏。

盡管已經(jīng)開展了大量干熱巖EGS 多物理場(chǎng)耦合的數(shù)值模擬工作，但由于改造后的滲透率通常極為復(fù)雜，常規(guī)方法難以描述，因此準(zhǔn)確地估計(jì)儲(chǔ)層中的熱量提取過(guò)程仍然存在大量挑戰(zhàn)。近年來(lái)，Xing 等[8-9]利用微震數(shù)據(jù)計(jì)算出地?zé)醿?chǔ)層的滲透率分布，為干熱巖開采設(shè)計(jì)了8 井布局，模擬結(jié)果展示了地?zé)醿?chǔ)層中流場(chǎng)和溫度場(chǎng)隨時(shí)間的演化。李庭樑等[15]結(jié)合地質(zhì)背景和壓裂數(shù)據(jù)，計(jì)算出滲透率等水動(dòng)力學(xué)參數(shù)的分布，并對(duì)熱儲(chǔ)進(jìn)行數(shù)值模擬，將其結(jié)果與分塊均質(zhì)熱儲(chǔ)進(jìn)行比較，結(jié)果表明熱儲(chǔ)的非均質(zhì)性會(huì)導(dǎo)致EGS 采熱性能顯著下降。Xu 等[18]利用水力壓裂過(guò)程中記錄的微震事件生成了裂縫模型，并將該方法應(yīng)用于Habanero 儲(chǔ)層，證明了 Habanero 井之間的連通性。Fang 等[19]分別利用基于立方定律的等效多孔介質(zhì)方法和離散裂縫網(wǎng)絡(luò)方法估算壓裂改造后的儲(chǔ)層滲透率，結(jié)果表明儲(chǔ)層中滲透率最高的區(qū)域即微震最密集的區(qū)域。以上研究表明水力壓裂誘發(fā)微地震的分布可以較為準(zhǔn)確地表征儲(chǔ)層改造后的滲透率，微地震分布優(yōu)勢(shì)的方向也即儲(chǔ)層中流體流動(dòng)的優(yōu)勢(shì)路徑。目前利用微地震數(shù)據(jù)評(píng)估滲透率也存在不足，如網(wǎng)格分布和震源精度的限制等。可通過(guò)優(yōu)化網(wǎng)格劃分、改善定位方法等提高滲透率評(píng)估的準(zhǔn)確性。

基于上述研究，本文以澳大利亞庫(kù)珀盆地的Habanero 熱儲(chǔ)為研究區(qū)，創(chuàng)新性地構(gòu)建了微地震監(jiān)測(cè)數(shù)據(jù)與儲(chǔ)層滲透率之間的等效關(guān)系并建立了不同的熱儲(chǔ)模型，包括分區(qū)均質(zhì)模型和非均質(zhì)模型；利用自主研發(fā)的多物理場(chǎng)耦合軟件PANDAS 對(duì)不同模型中滲流場(chǎng)、溫度場(chǎng)和生產(chǎn)井井底溫度的變化等進(jìn)行了數(shù)值模擬、分析比較，探究采熱性能的差異。

1 理論與方法

本文使用的模擬工具為自主研發(fā)的PANDAS，這是一種基于有限元和格子玻爾茲曼方法（Lattice Boltzmann Method，LBM）的數(shù)值模擬程序；其中熱-流動(dòng)耦合數(shù)值模擬模塊[20]將應(yīng)用于Habanero 地?zé)醿?chǔ)層開發(fā)中的有限元模擬及評(píng)估研究。該模型基本假設(shè)為：

（1）干熱巖儲(chǔ)層中充滿單相的飽和流體水，且流體在儲(chǔ)層中的流動(dòng)遵循達(dá)西定律；

（2）流體和巖石間的換熱通過(guò)對(duì)流和傳導(dǎo)實(shí)現(xiàn)；

（3）不考慮流體與礦物的物理化學(xué)反應(yīng)。

1.1 連續(xù)性方程

流體在多孔介質(zhì)中的流動(dòng)可以表示為：

式中： φ—孔隙度；

βm—綜合壓縮系數(shù)/Pa-1；

βf、βr—流體和基巖的壓縮系數(shù)/Pa-1；

P—流體壓力/Pa；

ρ—流體密度/（kg·m-3）；

v—流體流速/（m·s-1）。

遵循達(dá)西定律的單相流體的流動(dòng)控制方程為：

式中：k—多孔介質(zhì)的滲透率/m2；

μ—流體的動(dòng)力黏度/（Pa·s）；

g—重力加速度/（m·s-2）；

D—深度/m。

1.2 多孔介質(zhì)中的對(duì)流-傳導(dǎo)換熱方程

多孔介質(zhì)中基質(zhì)和流體能量平衡方程分別為：

式中：C—比熱/（J·kg-1·K-1）；

λ—導(dǎo)熱系數(shù)/（W·m-1·K-1）；

T—溫度/K；

r—巖石基質(zhì)；

f—流體。

從宏觀角度出發(fā)，假定EGS 各階段中局部處于熱平衡狀態(tài)，即Tr=Tf=T。聯(lián)立式（4）和式（5），多孔介質(zhì)中的傳導(dǎo)-對(duì)流換熱控制方程為：

式中：m—基質(zhì)和流體的混合系統(tǒng)。

2 研究區(qū)概況

庫(kù)珀盆地位于南澳大利亞?wèn)|北部，如圖1 所示。據(jù)資料顯示，庫(kù)珀盆地下方3 500～4 500 m 范圍內(nèi)埋藏有大量的高溫花崗巖，溫度約為250 °C[21]，是其他地區(qū)相同深度花崗巖溫度的2 倍?；◢弾r巖體上方覆蓋著約3 500 m 厚的沉積蓋層，其主要由低孔、低滲的碎屑巖和煤系夾層組成[12]，是花崗巖基底天然的“絕熱屏障”。

圖1 庫(kù)珀盆地?zé)醿?chǔ)層研究區(qū)域（修改自文獻(xiàn)[12，24]）Fig.1 Geothermal reservoir region to be studied in Cooper basin (modified from [12, 24])

Habanero 3 號(hào)井（H03）位于1 號(hào)井東北方向約560 m處，井深約4 221 m。在2008年對(duì)3 號(hào)井進(jìn)行了水力壓裂[22]。并對(duì)1 號(hào)井和3 號(hào)井進(jìn)行了時(shí)長(zhǎng)為42 d 的閉環(huán)流動(dòng)測(cè)試。實(shí)測(cè)數(shù)據(jù)表明，1 號(hào)井井口壓力為44.5 MPa，3 號(hào)井井口壓力為33.5 MPa，注采流量為15.5 kg/s。

水力壓裂誘發(fā)的微地震的空間分布表明存在呈近水平結(jié)構(gòu)的Habanero 斷層[23]。測(cè)井資料表明，該斷層是構(gòu)造成因，水力壓裂重新激活并擴(kuò)展了之前存在的裂縫，擴(kuò)大了儲(chǔ)層高滲區(qū)域[25]。本文將以庫(kù)珀盆地現(xiàn)場(chǎng)水力壓裂及相關(guān)聯(lián)通試驗(yàn)為基礎(chǔ)，開展數(shù)值模擬評(píng)估研究。

3 Habanero EGS 數(shù)值模型

3.1 三維數(shù)值模型

基于Habanero 地?zé)醿?chǔ)層的地質(zhì)構(gòu)造特征及水力壓裂改造效果，建立了三維數(shù)值模型。該模型在x，y和z方向的模型尺寸分別為1 200，1 500，500 m。模型在水力壓裂區(qū)劃分較密，網(wǎng)格尺寸為2.5 m，其他區(qū)域的網(wǎng)格尺寸由壓裂區(qū)向外逐漸變大，在15～20 m 之間。模型計(jì)算深度為地下4 000～4 500 m，如圖2（a）所示。

圖2 分區(qū)均質(zhì)模型和剖面圖Fig.2 Homogeneous model (a) and homogeneous model profile view (b)

3.1.1 分區(qū)均質(zhì)模型

根據(jù)庫(kù)珀盆地現(xiàn)場(chǎng)資料及研究成果，將研究區(qū)簡(jiǎn)化為從上至下分為3 層的三維模型，如圖2（b）所示。上層的深度為4 000～4 245 m，為未壓裂的致密花崗巖；中間層深度4 245～4 255 m，厚度為10 m，為孔隙度和滲透率較高的改造區(qū)；下層的深度為4 255～4 500 m，為未壓裂的致密花崗巖。

3.1.2 基于微地震數(shù)據(jù)的非均質(zhì)模型

微震可以反映改造過(guò)程中巖石破裂發(fā)生的位置及擴(kuò)展范圍，故本文利用已知地質(zhì)信息、微震監(jiān)測(cè)數(shù)據(jù)構(gòu)建Habanero 非均質(zhì)熱儲(chǔ)模型。本文的數(shù)據(jù)為Geodynamics 公司在Habanero 1 號(hào)井水力壓裂過(guò)程中記錄的誘發(fā)微震事件。研究表明，微震點(diǎn)的密集程度可以用來(lái)表征儲(chǔ)層內(nèi)的滲透率分布，一般認(rèn)為，微震分布密集區(qū)域的滲透率要高于微震分布稀疏區(qū)域[15]。因此，本文基于Habanero 熱儲(chǔ)區(qū)巖石破裂的機(jī)理，根據(jù)微震點(diǎn)到模型中各個(gè)節(jié)點(diǎn)的距離建立了三維非均質(zhì)儲(chǔ)層滲透率模型。

基于微地震數(shù)據(jù)的非均質(zhì)模型使用的網(wǎng)格模型與分區(qū)均質(zhì)模型使用的網(wǎng)格相同。在本文的計(jì)算方法中，滲透率是網(wǎng)格節(jié)點(diǎn)到微震點(diǎn)距離的指數(shù)函數(shù)，各個(gè)節(jié)點(diǎn)處的滲透率由2 部分組成：第1 部分為初始滲透率，第2 部分通過(guò)計(jì)算各個(gè)微震事件到該節(jié)點(diǎn)的距離進(jìn)行評(píng)估。具體計(jì)算方法引自Xing 等[8]。在有限元計(jì)算中對(duì)微震事件和滲透率計(jì)算節(jié)點(diǎn)的加權(quán)距離評(píng)估：

式中：K—節(jié)點(diǎn)處的滲透率/m2；

K0—初始滲透率/m2；

A—震級(jí)對(duì)滲透率的影響系數(shù)（此處A=1.0×10-15）；

dx、dy、dz—考慮破裂面傾角影響的微震事件與網(wǎng)格節(jié)點(diǎn)的距離；

B1、B2、B3—微震事件在x、y、z方向上影響范圍的參數(shù)。

3.2 模型參數(shù)

根據(jù)測(cè)井?dāng)?shù)據(jù)，庫(kù)珀盆地在3 500～5 000 m 深度范圍內(nèi)，初始溫度和壓力隨深度而變化，溫度梯度為31 °C/km，壓力梯度為8.6 MPa/km。將模型上下表面的壓力和溫度設(shè)置為恒定，上表面壓力為71.8 MPa，溫度為240.1 °C，下表面壓力為76.1 MPa，溫度為255.6 °C。根據(jù)Llanos 等[12]和O’Sullivan 等[26]的研究，模型基于所有側(cè)面沒(méi)有熱量或質(zhì)量的傳遞建立，因此模型四周設(shè)置為封閉邊界。

分區(qū)均質(zhì)模型中除滲透率以外，其他的巖石物理性質(zhì)參考Llanos 等[12]的研究，如表1 所示，在基于微地震數(shù)據(jù)的非均質(zhì)模型中，滲透率通過(guò)微震監(jiān)測(cè)數(shù)據(jù)進(jìn)行計(jì)算，其他巖石物理性質(zhì)參數(shù)與分區(qū)均質(zhì)模型保持一致。

表1 分區(qū)均質(zhì)模型中巖石物理性質(zhì)參數(shù)Table 1 Parameters of rock in the zonal homogenization model

在2 個(gè)模型中，均使用水作為循環(huán)工作流體。在模型初始狀態(tài)時(shí)，假設(shè)儲(chǔ)層中水的溫度與其所處的圍巖溫度相同，注入工作流體溫度為90 ℃。表2 為循環(huán)工作流體的物理性質(zhì)參數(shù)。

表2 循環(huán)工作流體的物理性質(zhì)參數(shù)Table 2 Parameters of the circulating fluid

3.3 模型滲透率的確定

利用2008年8月進(jìn)行的Habanero 1 號(hào)井（H01，注入井）與Habanero 3 號(hào)井（H03，生產(chǎn)井）（圖3）的閉環(huán)流動(dòng)現(xiàn)場(chǎng)測(cè)試數(shù)據(jù)，對(duì)研究區(qū)滲透率進(jìn)行校正。以現(xiàn)場(chǎng)實(shí)測(cè)壓力、流量數(shù)據(jù)為基礎(chǔ)，確定模型中的井底壓力，調(diào)整模型中改造區(qū)的滲透率，并對(duì)比模擬流量與實(shí)際流量。

圖3 模型中H01、H03 位置Fig.3 Locations of H01 and H03 in the model

3.3.1 分區(qū)均質(zhì)模型滲透率

Habanero 地?zé)醿?chǔ)層是致密花崗巖體，地層初始滲透率低，約為1.0×10-20m2。根據(jù)水力壓裂產(chǎn)生的微震點(diǎn)的空間分布，結(jié)合Llanos 等[12]的研究，Habanero 地?zé)醿?chǔ)層壓裂后滲透率存在各向異性，y方向和x方向的滲透率之比為2∶1，z方向滲透率較低，在此基礎(chǔ)上利用實(shí)測(cè)流量對(duì)滲透率進(jìn)行調(diào)整。表3 為在不同滲透率下，生產(chǎn)井中的模擬流量與實(shí)際流量的對(duì)比，由此可得，壓裂后改造區(qū)滲透率的最佳擬合值為x方向1.3×10-13m2，y方向2.6×10-13m2，z方向4.0×10-14m2，10 m 厚的改造區(qū)以外滲透率均為1.0×10-20m2。

表3 分區(qū)均質(zhì)模型不同滲透率下生產(chǎn)井中的模擬流量與實(shí)際流量Table 3 Simulated flow and actual flow of production wells in the homogeneous model with different permeabilities

3.3.2 非均質(zhì)模型滲透率

在非均質(zhì)模型中，通過(guò)調(diào)整相關(guān)參數(shù)，得到不同滲透率下生產(chǎn)井中的模擬流量（表4），并與實(shí)際流量進(jìn)行對(duì)比。模型中滲透率的最佳擬合相關(guān)參數(shù)為B1=5.0×10-5、B2=2.0×10-5、B3=2.5（B1、B2、B3是微震事件在x、y、z方向上影響范圍的權(quán)重參數(shù)，庫(kù)珀盆地水力壓裂微震事件以水平分布為主，因此設(shè)置Habanero地?zé)醿?chǔ)層壓裂效果對(duì)節(jié)點(diǎn)至微地震點(diǎn)x、y方向的距離較為敏感）。該情況下，滲透率分布見圖4 和圖5，圖4 顯示整體模型的滲透率的分布，圖5 為滲透率分布的等值面圖，主要顯示水力壓裂區(qū)的滲透率分布。

表4 非均質(zhì)模型不同滲透率相關(guān)參數(shù)下生產(chǎn)井中模擬流量與實(shí)際流量Table 4 Simulated flow and actual flow of production wells in the heterogeneous model with different permeability related parameters

圖4 （a）非均質(zhì)模型的滲透率分布圖和（b）滲透率剖面圖Fig.4 Permeability distribution for the heterogeneous model (a)and permeability profile (b)

圖5 （a）滲透率分布等值面圖和（b）滲透率等值面的剖面圖Fig.5 Contour surface view of permeability distribution (a) and isosurface profile of permeability distribution (b)

4 結(jié)果與討論

4.1 流場(chǎng)

圖6（a）（d）是2 種模型中深度為4 250 m 截面的流速分布，可以看出，2 種滲透率模型中流速較大的是2 口井之間的連線及其附近的區(qū)域，說(shuō)明儲(chǔ)層內(nèi)流體流動(dòng)的優(yōu)勢(shì)路徑為兩井之間的連通區(qū)域。對(duì)比2 種滲透率模型的平面流速，分區(qū)均質(zhì)模型中該深度截面流速較大區(qū)域的范圍要比非均質(zhì)模型更大。原因在于2 種模型滲透率的不同，分區(qū)均質(zhì)模型10 m 厚改造區(qū)內(nèi)是具有各向異性的均一分布滲透率，x方向1.30×10-13m2，y方向2.60×10-13m2，z方向4.0×10-14m2。非均質(zhì)模型中滲透率由微震事件及破裂面決定，如圖7 所示，注入井附近滲透率最大，約為9.4 ×10-13m2，從注入井向生產(chǎn)井方向滲透率逐漸變小，在生產(chǎn)井附近滲透率約為1.0×10-14m2，這是因?yàn)槲⒄鹗录l(fā)生的集中區(qū)域?yàn)樽⑷刖浇?，并由此向外延伸，在微震事件發(fā)生較多的區(qū)域滲透率較高，儲(chǔ)層中2 井連線及其附近的區(qū)域表現(xiàn)出流動(dòng)優(yōu)勢(shì)，遠(yuǎn)離注入井的區(qū)域由于微震事件較少，所以裂隙較少、滲透率較低，流速也小。非均質(zhì)模型中除注入井附近，其他位置滲透率要比分區(qū)均質(zhì)模型中低，因此分區(qū)均質(zhì)模型井間流速大于非均質(zhì)模型。

圖6 不同模型的流速分布（上為分區(qū)均質(zhì)模型，下為非均質(zhì)模型）Fig.6 Velocity distribution of different models (the upper is the zonal homogeneous model and the bottom is the heterogeneous model)

圖7 非均質(zhì)模型中z=4 250 m 剖面滲透率分布Fig.7 Permeability distribution at the section z=4 250 m in the heterogeneous model

由圖6（b）（e）可知，在垂直方向上，分區(qū)均質(zhì)模型的流動(dòng)區(qū)域限于10 m 厚的改造區(qū)內(nèi)，而非均質(zhì)模型在深度方向上流體流動(dòng)區(qū)域集中于該區(qū)域，但不限于此，這是因?yàn)槲⒄鹗录诖瓜蛏嫌休^大延伸。

4.2 溫度場(chǎng)

溫度分布等值面圖，范圍為90～230 °C（圖8），總的來(lái)看，在初始運(yùn)行階段，冷流體自注入井進(jìn)入壓裂改造區(qū)，由于與周圍高溫巖體進(jìn)行了熱交換，注入井附近巖體溫度接近注入溫度90 °C。隨著注入流體的增多，低溫區(qū)（90～115 °C）逐漸擴(kuò)大，并沿優(yōu)勢(shì)流動(dòng)通道向生產(chǎn)井?dāng)U散。

圖8 不同模型溫度隨時(shí)間變化的等值面圖Fig.8 Isosurface of temperature variation with time in different models

對(duì)比來(lái)看，分區(qū)均質(zhì)模型中的冷流體從注入井向改造區(qū)中流動(dòng)擴(kuò)散，逐漸形成一個(gè)近似扁橢圓體的低溫區(qū)，近y方向上的低溫區(qū)影響范圍大于近x方向的影響范圍。在深度方向上，低溫區(qū)范圍基本限制在深度4 245～4 255 m 的10 m 厚改造區(qū)內(nèi)。

在非均質(zhì)模型中，注入井兩側(cè)冷流體流動(dòng)的影響范圍與均質(zhì)模型相比，冷流體向著生產(chǎn)井流動(dòng)的趨勢(shì)更加明顯。近y方向上的低溫區(qū)影響范圍同樣大于近x方向的影響范圍。但深度方向?qū)Ρ葋?lái)看，低溫區(qū)影響范圍不再局限于10 m 厚的范圍。這是因?yàn)樗毫研纬傻奈⒄鹗录谏疃壬系姆植挤秶^廣，微震發(fā)生的區(qū)域內(nèi)儲(chǔ)層滲透率也隨之變大。

4.3 采熱性能

2 種模型的生產(chǎn)井井底溫度在30 a 間的變化情況見圖9（a）。分區(qū)均質(zhì)模型30 a 間溫度變化趨勢(shì)為前期（0～7.9 a）保持不變，之后（7.9～30 a）緩慢下降，非均質(zhì)模型30 a 間溫度變化趨勢(shì)為前期（0～8.5 a）保持不變，之后（8.5～30 a）快速下降，且下降速度逐漸變慢。系統(tǒng)運(yùn)行到30 a，分區(qū)均質(zhì)模型溫度下降到220.9 °C，下降了26.9 °C，非均質(zhì)模型溫度下降到170.0 °C，下降了77.8 °C。

圖9 分區(qū)均質(zhì)和非均質(zhì)模型生產(chǎn)井井底溫度和采熱量對(duì)比Fig.9 Comparison of the bottom hole temperature of production wells and heat extraction comparison in the two models

對(duì)比來(lái)看，2 種模型生產(chǎn)井初始溫度保持不變的時(shí)間十分接近；后期分區(qū)均質(zhì)模型溫度下降少，非均質(zhì)模型下降多。由此分析，分區(qū)均質(zhì)模型略早下降0.6 a，即熱突破時(shí)間稍微早一些，這是分區(qū)均質(zhì)模型生產(chǎn)井附近滲透率比非均質(zhì)模型大導(dǎo)致的。非均質(zhì)模型溫度下降幅度比分區(qū)均質(zhì)模型大，下降速度比非均質(zhì)模型快，這是由于非均質(zhì)模型的高滲區(qū)集中于2 口井連通區(qū)域內(nèi)，使得非均質(zhì)模型中優(yōu)勢(shì)流動(dòng)通道顯著，所以溫度下降更快。

2 種模型單位質(zhì)量流體采熱量的對(duì)比關(guān)系，即采出熱量與采出質(zhì)量的比值[27]，見圖9（b）。分區(qū)均質(zhì)模型和非均質(zhì)模型在儲(chǔ)層模型運(yùn)行初期單位質(zhì)量流體采熱量比較接近，約為1.04 MJ/kg，分區(qū)均質(zhì)模型初始時(shí)間略高。之后隨著生產(chǎn)井溫度的快速下降，2 種模型的單位質(zhì)量流體采熱量也迅速降低，對(duì)比來(lái)看，兩種模型單位質(zhì)量流體采熱量下降發(fā)生的時(shí)間接近，分區(qū)均質(zhì)模型采熱量下降速度很慢，在模型計(jì)算到30 a時(shí)，單位質(zhì)量流體采熱量降低至1.003 MJ/kg，降低了3.74%；非均質(zhì)模型采熱量下降速度較快，在模型計(jì)算到30 a 時(shí)，采熱量降低至0.906 MJ/kg，降低了12.72%。從采熱量對(duì)比看，分區(qū)均質(zhì)模型在熱儲(chǔ)層運(yùn)行期間采熱量變化相對(duì)更穩(wěn)定，非均質(zhì)模型采熱量降幅較大。

對(duì)比生產(chǎn)井溫度曲線和采熱量曲線，分區(qū)均質(zhì)模型的模擬結(jié)果相比于非均質(zhì)模型結(jié)果來(lái)說(shuō)，溫度下降幅度小、單位質(zhì)量流體采熱量高，但由微震分布可知實(shí)際儲(chǔ)層中的滲透率分布并不均勻，2 井之間優(yōu)勢(shì)流動(dòng)通道效應(yīng)顯著，因此分區(qū)均質(zhì)模型的模擬單位質(zhì)量流體采熱量高于實(shí)際單位質(zhì)量流體采熱量。

原始地?zé)醿?chǔ)層壓裂后，激活并擴(kuò)展了先存裂縫，在水壓足夠高時(shí)產(chǎn)生新的裂縫，提高了儲(chǔ)層中尤其是近井區(qū)域的滲透率，形成優(yōu)勢(shì)流動(dòng)通道，改善了儲(chǔ)層連通性。此外，水力壓裂后儲(chǔ)層滲透率具有顯著的滲透率各向異性，是影響干熱巖開發(fā)動(dòng)態(tài)的重要因素之一[28]。相比分區(qū)均質(zhì)模型，非均質(zhì)模型中的優(yōu)勢(shì)流動(dòng)通道、滲透率各向異性更加明顯，更加符合壓裂后的儲(chǔ)層實(shí)際情況，因此在實(shí)際應(yīng)用中，非均質(zhì)模型的模擬結(jié)果對(duì)實(shí)際工程更具參考意義。

5 結(jié)論

（1）水力壓裂誘發(fā)的微震事件集中于近井區(qū)域，導(dǎo)致該區(qū)滲透率增高，因此非均質(zhì)模型中優(yōu)勢(shì)流動(dòng)通道顯著，冷流體從注入井更快流向生產(chǎn)井，溫度下降速度快；分區(qū)均質(zhì)模型中優(yōu)勢(shì)流動(dòng)通道效應(yīng)較弱，溫度下降速度慢。

（2）深度方向上，分區(qū)均質(zhì)模型僅限于10 m 厚的壓裂改造區(qū)內(nèi)流動(dòng)；非均質(zhì)模型由于微震事件在深度方向上有較大延伸，流體流動(dòng)區(qū)域集中于該層但不限于該范圍。

（3）分區(qū)均質(zhì)模型在熱儲(chǔ)層開采模擬期間采熱量變化相對(duì)穩(wěn)定，降幅為3.74%；非均質(zhì)模型單位質(zhì)量流體采熱量降幅較大，為12.72%。

（4）對(duì)比2 種模型的模擬結(jié)果，分區(qū)均質(zhì)模型的生產(chǎn)井溫度下降幅度小、單位質(zhì)量流體采熱量高。但由微震監(jiān)測(cè)數(shù)據(jù)可知，實(shí)際儲(chǔ)層中的滲透率分布不均，優(yōu)勢(shì)流動(dòng)通道效應(yīng)顯著，分區(qū)均質(zhì)模型的模擬采熱量高于實(shí)際采熱量，因此非均質(zhì)模型的模擬結(jié)果在干熱巖資源實(shí)際開發(fā)過(guò)程中更具指導(dǎo)意義。