曲占慶, 張 偉, 郭天魁, 孫 江, 鞏法成, 田 雨, 李小龍

(中國石油大學(xué)(華東)石油工程學(xué)院,山東青島 266580)

干熱巖中地?zé)豳Y源豐富,但低孔低滲、埋深大等特點使開發(fā)難度增大。利用壓裂技術(shù)形成增強(qiáng)型地?zé)嵯到y(tǒng)為從干熱巖中高效提取熱能提供了可能?，F(xiàn)場施工試驗[1-2]、數(shù)值模擬研究[3]、物理模擬試驗[4]均表明低溫壓裂液注入高溫干熱巖儲層時的冷沖擊作用促進(jìn)了多裂縫的形成與裂縫網(wǎng)絡(luò)的構(gòu)建。針對干熱巖儲層改造后形成增強(qiáng)型地?zé)嵯到y(tǒng)的采熱性能,國內(nèi)外學(xué)者開展了溫度-滲流耦合[5-6]與溫度-滲流-應(yīng)力(thermo-hydro-mechanical,THM)耦合[7-9]采熱性能的相關(guān)研究。然而大部分?jǐn)?shù)值模擬研究基于熱平衡理論,忽略了熱儲層中巖體與流體間的熱交換,Sun等[10]在前期研究的基礎(chǔ)上基于局部熱非平衡理論開展了相關(guān)研究,但其中采熱性能僅對產(chǎn)出流體溫度進(jìn)行描述評價指標(biāo)較單一。 筆者以產(chǎn)出液質(zhì)量流量、產(chǎn)出液溫度、采熱速率為評價指標(biāo),基于局部熱非平衡理論,將儲層視為由基質(zhì)巖體與離散裂縫組成的雙重介質(zhì)模型,建立THM耦合模擬采熱過程中流體流動、熱量傳遞、巖體變形的相互作用,借助有限元軟件COMSOL Multiphysics實現(xiàn)所建立模型的全耦合求解,進(jìn)而對不同注采井網(wǎng)和模型參數(shù)下采熱性能進(jìn)行探討。

1 耦合模型的建立與求解

1.1 模型假設(shè)條件

假設(shè):含裂縫網(wǎng)絡(luò)干熱巖儲層為由基質(zhì)巖塊與離散裂縫構(gòu)成的雙重介質(zhì)模型,改造后形成的裂縫網(wǎng)絡(luò)與基質(zhì)相比滲透率較高,因此構(gòu)成了攜熱流體的主要流動通道,基質(zhì)和裂縫中的滲流規(guī)律符合達(dá)西定律;流體和巖石間不發(fā)生化學(xué)反應(yīng);導(dǎo)熱符合傅里葉定律,不考慮熱輻射的影響;儲層壓力高,不存在相變?yōu)閱蜗嘁后w流動;基于小變形理論認(rèn)為基巖為熱彈性。

1.2 控制方程

根據(jù)以上假設(shè),基于局部熱非平衡理論,含裂縫網(wǎng)絡(luò)干熱巖開發(fā)過程中THM耦合數(shù)學(xué)模型控制方程有質(zhì)量守恒方程、能量守恒方程和應(yīng)力平衡方程。

1.2.1 質(zhì)量守恒方程

巖體內(nèi)質(zhì)量守恒方程為

(1)

其中

式中,ρ為巖體密度,kg/m3;S為巖體的儲水系數(shù);t為時間,a;Qm為滲流源匯項,kg/(m3·s);v為巖塊中水流速,m/s;k為巖體滲透率,m2;μ為流體動力黏度,Pa·s;p為壓力,Pa。

裂縫內(nèi)質(zhì)量守恒方程為

(2)

其中

式中,df為裂縫寬度,m;ρf為裂縫密度,kg/m3;Sf為裂縫儲水系數(shù),Pa-1;Qf為巖塊與裂隙面的流量交換;kf為裂縫滲透率,m2;τ沿裂縫切向梯度。

1.2.2 能量守恒方程

對于增強(qiáng)型地?zé)嵯到y(tǒng)中流體的換熱,工作流體與固體骨架間的溫差不容忽視,因此局部熱平衡模型將導(dǎo)致較大的計算誤差。采用兩個不同的能量方程分別對流體和固體骨架中的熱量傳遞進(jìn)行描述實現(xiàn)非平衡傳熱,從而更加真實準(zhǔn)確地模擬基質(zhì)巖體與裂縫內(nèi)流體的能量傳遞過程。

流體能量守恒方程[11]為

(3)

固體骨架能量守恒方程為

(4)

式中,C為比熱容,J/(kg·K);λ為導(dǎo)熱系數(shù),W/(m·K);ρ為密度,kg/m3;hsf為固體骨架與液體間換熱系數(shù),W/(m2·K);αsf為兩相界面的比表面積,m2/g;Ts為巖石骨架溫度,K;Tf為流體溫度,K;下標(biāo)s、f分別對應(yīng)巖石和流體;ε為巖石孔隙率。

裂縫內(nèi)熱傳導(dǎo)與熱對流能量平衡方程[10]為

(5)

其中

(ρC)eff=ρfCf(1-φ)+ρfrCfrφ,keff=λf(1-φ)+λfrφ.

式中,dfr為裂隙寬度,m;(ρC)eff為有效比熱容,J/(kg·K);Keff為有效熱傳導(dǎo)系數(shù),W/(m·K);v為巖塊中水流速,m/s,流體遵從達(dá)西定律。

1.2.3 應(yīng)力平衡方程

(6)

其中

G=E/2(1+ν).

裂隙近似由一對法向和切向位移表面構(gòu)成,巖體裂隙變形方程為

(7)

(8)

式中,u為位移,m;σ為總應(yīng)力,Pa;σ′為有效應(yīng)力,Pa;K為剛度,Pa/m;下標(biāo)n、s分別指示裂隙面法向和切向。通過對裂縫內(nèi)部邊界進(jìn)行弱形式描述代表裂縫的剛度[5,10]。

1.3 THM耦合相互作用關(guān)系

(1)滲流場對溫度場的影響:裂縫內(nèi)冷流體與高溫巖體進(jìn)行熱交換,巖體中原有平衡狀態(tài)的溫度場發(fā)生改變。

(2)溫度場對滲流場的影響:高溫高壓巖體內(nèi),裂隙內(nèi)水滲流物理特性(如密度、黏度等)不再是常數(shù),隨溫度變化發(fā)生變化,從而影響巖體內(nèi)流體流動。

(3)溫度場對應(yīng)力場的影響:巖體溫度的改變影響巖體固有的物理力學(xué)性質(zhì),同時溫度效應(yīng)產(chǎn)生的熱應(yīng)力會導(dǎo)致原有應(yīng)力場分布改變。Zhu等[12]通過實驗獲得了不同溫度下彈性模量。

(4)應(yīng)力場對溫度場的影響:應(yīng)力場變化引起巖體結(jié)構(gòu)改變時的變形生熱,同時巖體的結(jié)構(gòu)變形會影響孔隙內(nèi)部的導(dǎo)熱性能。

(5)滲流場對應(yīng)力場的影響:孔隙壓力的存在影響巖石有效應(yīng)力。

(6)應(yīng)力場對滲流場的影響:應(yīng)力場變化導(dǎo)致了孔隙和裂隙變化,進(jìn)而改變了巖體與裂縫的滲透性能。

與裂縫的滲透率相比巖石基質(zhì)的滲透率極低,因此忽略基質(zhì)滲透率的變化。裂縫滲透率的變化進(jìn)行描述為

(9)

1.4 模型的求解

應(yīng)用有限元軟件COMSOL Multiphysics可將各物理場聯(lián)立求解,實現(xiàn)耦合分析中各物理場的全面耦合,從而反映真實的模擬過程。采用有限元方法對裂縫網(wǎng)絡(luò)的模擬通常使用具有一定厚度的單元層刻畫裂縫,但當(dāng)裂縫網(wǎng)絡(luò)規(guī)模較大時,往往受限于計算能力而無法求解。借鑒文獻(xiàn)[10]和[13],將二維離散裂縫網(wǎng)絡(luò)設(shè)置為線單元,并借助軟件中內(nèi)置的裂隙模塊進(jìn)行求解。所建立耦合模型與借助COMSOL實現(xiàn)的數(shù)值求解方法的可行性與準(zhǔn)確性通過對比白冰[14]針對一維條件下飽和土熱彈性固結(jié)問題獲得的解析解得到驗證。

2 計算模型

2.1 模型描述

假定壓裂改造后含裂縫網(wǎng)絡(luò)干熱巖計算模型為700 m×700 m的二維模型。改造后干熱巖中裂縫網(wǎng)絡(luò)認(rèn)為是在注入水壓力、熱應(yīng)力與儲層中原有裂縫的共同作用下所形成的[1-4,10],模型中采用隨機(jī)裂縫進(jìn)行描述。隨機(jī)生成的裂縫網(wǎng)絡(luò)由兩組裂隙構(gòu)成,裂縫數(shù)量為245條,傾角分別為0°和90°,裂縫長度均值為70 m,方差為30 m,由于增強(qiáng)型地?zé)嵯到y(tǒng)中的主要導(dǎo)水通道為聯(lián)通裂縫,在構(gòu)建裂縫網(wǎng)絡(luò)模型時去掉不聯(lián)通的裂縫[13]。模型采用兩注一采的注采方式,計算模型如圖1所示,藍(lán)色圓點代表注入井,紅色圓點代表采出井,注入井與采出井距離均為350 m。

圖1 含裂縫網(wǎng)絡(luò)干熱巖模型Fig.1 Hot dry rock model with fracture network

2.2 計算條件

利用上述模型模擬含裂縫網(wǎng)絡(luò)干熱巖采熱全過程,耦合分析從注入水的時刻開始,時間步長為1 d,總時間段為40 a。邊界條件和初始條件如下:

(1)滲流場。為確保系統(tǒng)中的水循環(huán)給定壓力邊界條件。假設(shè)注入井給定壓力為44 MPa,生產(chǎn)井給定壓力為40 MPa,注采井直徑設(shè)置為0.2 m。其他外部條件為不可滲透邊界。在攜熱流體注入前,儲層中的初始水壓為40 MPa。

(2)溫度場。注入井的水溫為20 ℃。為真實模擬EGS周邊基巖對EGS區(qū)域的傳熱將邊界設(shè)為開邊界。對于干熱巖基質(zhì)巖體初始溫度為200 ℃。

(3)應(yīng)力場。模型四周邊界設(shè)置為位移約束邊界。為便于明確低溫流體注入過程對儲層應(yīng)力場的擾動,不考慮地應(yīng)力場的影響。

應(yīng)用COMSOL Multiphysics實現(xiàn)了固體形變、流體滲流與熱量傳遞的全耦合,并對所建立的二維模型進(jìn)行求解。采用三角形網(wǎng)格進(jìn)行網(wǎng)格剖分,形成三角形單元49 725個與邊界單元5 932個。

模型中采用的主要材料參數(shù)如下:基巖密度為2 700 kg/m3,滲透率為1×10-17m2,孔隙度為1%,導(dǎo)熱系數(shù)為3.1 W/(m·K),彈性模量為40 GPa,泊松比為0.25,比熱容為950 J/(kg·K),儲水率為1×10-8Pa-1,熱膨脹系數(shù)為6×10-6K-1;裂縫密度為1 200 kg/m3,寬度為0.002 m,滲透率為3×10-11m2,孔隙度為30%,儲水率為1×10-9Pa-1,比熱容為850 J/(kg·K),法向剛度為100 GPa/m,切向剛度為30 GPa/m;Biot系數(shù)為1。

3 模擬結(jié)果分析

3.1 滲流場、溫度場、應(yīng)力場演化

向含裂縫網(wǎng)絡(luò)干熱巖儲層中注入低溫攜熱流體循環(huán)采熱過程伴隨著滲流場、溫度場、應(yīng)力場的互相影響與制約。由于基巖滲透率極低,少量攜熱流體濾失到基巖中,如圖2所示。流體主要沿裂縫滲流,聯(lián)通的裂縫網(wǎng)絡(luò)成為EGS中主要的導(dǎo)水通道。同時攜熱流體在流動過程中與裂縫周圍基巖發(fā)生熱量交換,如圖3所示。低溫流體被不斷加熱,裂縫周圍基巖中熱量被不斷提取,攜熱流體流動一段距離后被加熱至儲層溫度,此時的流體與基巖達(dá)到熱平衡不再發(fā)生熱交換。不論是流體滲流產(chǎn)生的水壓作用還是注入流體與儲層溫差導(dǎo)致的熱應(yīng)力作用都將對應(yīng)力場分布形成擾動,應(yīng)力場的變化將直接影響儲層中裂縫的滲流能力。在裂縫網(wǎng)絡(luò)中設(shè)置裂縫滲透率探針,如圖4所示。低溫流體的注入使高溫基巖遇冷收縮,同時流體壓力也導(dǎo)致裂縫和基巖受力發(fā)生形變,在兩者的共同作用下裂縫滲透率提高。由于a、d兩點處離注入井距離近,開發(fā)初期即受到水壓力與熱應(yīng)力的影響,裂縫滲透率升高較早,而b、c處距離注入井遠(yuǎn),水壓力與冷鋒面擴(kuò)散至b、c處需一段時間,因此系統(tǒng)開發(fā)約4 a后滲透率才有所提升,且此處流體已被加熱,熱應(yīng)力的作用不及a、d處劇烈,導(dǎo)致提升幅度不大。

裂縫分布存在非均勻性,裂縫密集區(qū)域的熱對流更明顯,冷鋒面范圍更大,如圖3所示。裂縫稀疏區(qū)域,周圍基巖中未提取熱量多,基巖對攜熱流體加熱能力更強(qiáng),冷鋒面范圍較小。裂縫分布的非均勻性也同樣導(dǎo)致了儲層有效應(yīng)力場的非均勻性,如圖5所示。

隨著開發(fā)的進(jìn)行,冷鋒面逐漸向采出井推進(jìn),直至該EGS系統(tǒng)出現(xiàn)熱突破,如圖3(c)所示。裂縫周圍基巖溫度的降低也使裂縫內(nèi)流體黏度增加流速下降(圖2)。冷鋒向采出井的不斷推進(jìn)也導(dǎo)致了有效應(yīng)力場的不斷變化。因此在EGS開發(fā)中流體滲流、熱量傳遞、固體形變的相互作用伴隨著整個開發(fā)過程。

圖2 不同時刻裂縫網(wǎng)絡(luò)內(nèi)速度場分布Fig.2 Distribution of velocity field in fracture network at different time

圖3 不同時刻干熱巖儲層溫度場分布Fig.3 Distribution of temperature field in hot reservoir at different time

圖4 裂縫滲透率探針分布與不同時刻滲透率演化Fig.4 Location of fracture permeability probes and permeability evolution

圖5 不同時刻干熱巖儲層有效應(yīng)力場等值線圖Fig.5 Contour map of effective stress field in hot reservoir at different time

3.2 注采井網(wǎng)對采熱性能的影響

注采井網(wǎng)的變化必然對滲流場、溫度場、應(yīng)力場產(chǎn)生影響,最終改變干熱巖儲層的采熱性能。采用圖1中裂縫網(wǎng)絡(luò)模型,分別設(shè)置一注一采、兩注一采、四注一采3種不同的注采井網(wǎng),注采井井距均設(shè)置為350 m。以采出井的產(chǎn)出液溫度、產(chǎn)出液質(zhì)量流量、采熱速率[15-17]作為采熱性能評價指標(biāo)，產(chǎn)熱速率為

Wh=q(hpro-hinj).

(10)

式中,Wh為產(chǎn)熱速率,MW;q為產(chǎn)出液質(zhì)量流量,kg/s;hpro為產(chǎn)出液比焓,kJ/kg;hinj為注入液比焓,kJ/kg。

圖6 不同注采井網(wǎng)下采熱性能Fig.6 Heat mining performance with different inject-product well patterns

井網(wǎng)中注入井?dāng)?shù)量的增加可為EGS注入更多的攜熱流體,在裂縫網(wǎng)絡(luò)中進(jìn)行熱交換后,更多熱流體在采出井被采出。如圖6(a)所示,一注一采開發(fā)5 a時產(chǎn)出液質(zhì)量流量為32 kg/s,井網(wǎng)改為兩注一采后開發(fā)5 a時產(chǎn)出液質(zhì)量流量為55 kg/s,四注一采開發(fā)5 a時達(dá)到了75 kg/s。不同注采井網(wǎng)干熱巖儲層溫度場分布如圖7所示。注入井的加入,使干熱巖中更多區(qū)域的地?zé)豳Y源得到動用,同時合理的注采井布局并未造成冷鋒面擴(kuò)大導(dǎo)致的熱突破加快。如圖8(b)所示,注入井的加入反而緩解了產(chǎn)出液溫度的降低。分析圖8可見,裂縫網(wǎng)絡(luò)中同一位置處,采用一注一采時裂縫內(nèi)流速最快,增加注入井?dāng)?shù)量后,在采出井附近產(chǎn)生流體累加,裂縫內(nèi)流速減緩,增加了其在干熱巖中的換熱從而緩解了產(chǎn)出液溫度的降低。注入井?dāng)?shù)量的增加使產(chǎn)出液質(zhì)量流量大幅提高,同時產(chǎn)出液溫度未出現(xiàn)下降,因此EGS采熱速率提升明顯,四注一采開發(fā)5 a時采熱速率達(dá)到60 MW。可見在井網(wǎng)中加入注入井對于系統(tǒng)采熱速率提高有效。

圖7 開發(fā)20 a時不同注采井網(wǎng)干熱巖儲層溫度場分布Fig.7 Distribution of temperature field in hot reservoir with different inject-product well pattern after 20-year development

圖8 開發(fā)20 a時不同注采井網(wǎng)干熱巖儲層速度場分布Fig.8 Distribution of velocity field in fracture network with different inject-product well pattern after 20-year development

3.3 模型參數(shù)對采熱性能的影響

采熱速率結(jié)合了產(chǎn)出液質(zhì)量流量與產(chǎn)出液溫度,是評價系統(tǒng)采熱性能的綜合評價指標(biāo)。鑒于對注采井網(wǎng)研究的認(rèn)識,采用四注一采的開發(fā)方式,針對模型關(guān)鍵參數(shù)開展敏感性分析。

圖9 儲層物性對采熱速率的影響Fig.9 Effect of reservoir properties on heat mining rate

基質(zhì)滲透率的增加首先引起滲流場的改變,更多低溫攜熱流體由裂縫流向基質(zhì),高溫基巖遇冷收縮加劇。因此裂縫開度增加,應(yīng)力場受到影響,同時裂縫滲透率提高,產(chǎn)出液質(zhì)量流量增加采熱速率提高,如圖9(a)所示,基質(zhì)滲透率為1×10-17m2與1×10-18m2時相比采熱速率增幅較大,但與1×10-16m2時相比增幅較小,這是由于基巖滲透率增大到一定程度后,大量工作流體濾失所致。

基質(zhì)導(dǎo)熱系數(shù)通過改變溫度場從而影響滲流場與應(yīng)力場?；|(zhì)導(dǎo)熱系數(shù)對儲層溫度場影響較小,如圖9(b)所示。其對采熱速率的影響也較微弱,該模擬結(jié)果與Sun等[10]的研究結(jié)果一致。

基質(zhì)熱膨脹系數(shù)的不同首先影響應(yīng)力場,熱膨脹系數(shù)越大,高溫基巖遇冷產(chǎn)生的熱應(yīng)力越強(qiáng),基巖收縮位移越大,裂縫滲透率提升越大,如圖9(c)所示。熱膨脹系數(shù)增加開發(fā)0～30 a間采熱速率增加,但30 a后采熱速率降低。這是由于低溫流體注入產(chǎn)生的熱應(yīng)力作用在使裂縫滲透率提高的同時也導(dǎo)致了熱竄流[7],使熱突破更快,產(chǎn)出液溫度降低速度增加。

注采壓差增加儲層中流體流速增加,采熱速率提高,同時流速的增加將導(dǎo)致攜熱流體換熱不充分,如圖10所示,隨著開發(fā)時間的延長由注采壓差增加引起的采熱速率提升幅度減小。若壓差過大將導(dǎo)致EGS開發(fā)壽命縮短,采熱速率快速下降。不同注入溫度對應(yīng)力場與滲流場均有影響。注入溫度越高低溫流體導(dǎo)致的熱應(yīng)力作用越小,高流體溫度也更有利于流體流動,因此注入流體溫度越高采熱速率越高。

4 結(jié) 論

(1)向含裂縫網(wǎng)絡(luò)干熱巖中注入攜熱流體采熱過程中涉及THM耦合,其中應(yīng)力場受水壓力與溫差所致熱應(yīng)力的共同作用,裂縫網(wǎng)絡(luò)分布的非均勻性對流體滲流、冷鋒面擴(kuò)散及應(yīng)力場分布均有影響。

(2)采用局部熱非平衡理論對干熱巖采熱過程進(jìn)行模擬,考慮換熱過程中工作流體與固體骨架間的溫度差異可更真實地描述巖石與流體間熱交換。

(3)井網(wǎng)中注入井的加入,使干熱巖中更多區(qū)域的地?zé)豳Y源得到動用,產(chǎn)出液質(zhì)量流量大幅提升,合理的布局可避免熱突破加快,從而有效提高系統(tǒng)采熱速率,四注一采開發(fā)5 a時采熱速率可達(dá)60 MW。

(4)儲層基質(zhì)滲透率對采熱速率的影響大于基質(zhì)導(dǎo)熱系數(shù)。基質(zhì)熱膨脹系數(shù)通過控制溫差引起的熱應(yīng)力影響采熱性能,其數(shù)值越大攜熱流體越容易發(fā)生熱竄流。注入溫度越高對干熱巖開發(fā)越有利,注采壓差應(yīng)調(diào)節(jié)適當(dāng),過大的壓差雖然產(chǎn)出液質(zhì)量流量高,但會影響開發(fā)壽命。