曲占慶, 張 偉, 郭天魁, 孫 江, 鞏法成, 田 雨, 李小龍
(中國(guó)石油大學(xué)(華東)石油工程學(xué)院,山東青島 266580)
干熱巖中地?zé)豳Y源豐富,但低孔低滲、埋深大等特點(diǎn)使開(kāi)發(fā)難度增大。利用壓裂技術(shù)形成增強(qiáng)型地?zé)嵯到y(tǒng)為從干熱巖中高效提取熱能提供了可能。現(xiàn)場(chǎng)施工試驗(yàn)[1-2]、數(shù)值模擬研究[3]、物理模擬試驗(yàn)[4]均表明低溫壓裂液注入高溫干熱巖儲(chǔ)層時(shí)的冷沖擊作用促進(jìn)了多裂縫的形成與裂縫網(wǎng)絡(luò)的構(gòu)建。針對(duì)干熱巖儲(chǔ)層改造后形成增強(qiáng)型地?zé)嵯到y(tǒng)的采熱性能,國(guó)內(nèi)外學(xué)者開(kāi)展了溫度-滲流耦合[5-6]與溫度-滲流-應(yīng)力(thermo-hydro-mechanical,THM)耦合[7-9]采熱性能的相關(guān)研究。然而大部分?jǐn)?shù)值模擬研究基于熱平衡理論,忽略了熱儲(chǔ)層中巖體與流體間的熱交換,Sun等[10]在前期研究的基礎(chǔ)上基于局部熱非平衡理論開(kāi)展了相關(guān)研究,但其中采熱性能僅對(duì)產(chǎn)出流體溫度進(jìn)行描述評(píng)價(jià)指標(biāo)較單一。 筆者以產(chǎn)出液質(zhì)量流量、產(chǎn)出液溫度、采熱速率為評(píng)價(jià)指標(biāo),基于局部熱非平衡理論,將儲(chǔ)層視為由基質(zhì)巖體與離散裂縫組成的雙重介質(zhì)模型,建立THM耦合模擬采熱過(guò)程中流體流動(dòng)、熱量傳遞、巖體變形的相互作用,借助有限元軟件COMSOL Multiphysics實(shí)現(xiàn)所建立模型的全耦合求解,進(jìn)而對(duì)不同注采井網(wǎng)和模型參數(shù)下采熱性能進(jìn)行探討。
假設(shè):含裂縫網(wǎng)絡(luò)干熱巖儲(chǔ)層為由基質(zhì)巖塊與離散裂縫構(gòu)成的雙重介質(zhì)模型,改造后形成的裂縫網(wǎng)絡(luò)與基質(zhì)相比滲透率較高,因此構(gòu)成了攜熱流體的主要流動(dòng)通道,基質(zhì)和裂縫中的滲流規(guī)律符合達(dá)西定律;流體和巖石間不發(fā)生化學(xué)反應(yīng);導(dǎo)熱符合傅里葉定律,不考慮熱輻射的影響;儲(chǔ)層壓力高,不存在相變?yōu)閱蜗嘁后w流動(dòng);基于小變形理論認(rèn)為基巖為熱彈性。
根據(jù)以上假設(shè),基于局部熱非平衡理論,含裂縫網(wǎng)絡(luò)干熱巖開(kāi)發(fā)過(guò)程中THM耦合數(shù)學(xué)模型控制方程有質(zhì)量守恒方程、能量守恒方程和應(yīng)力平衡方程。
1.2.1 質(zhì)量守恒方程
巖體內(nèi)質(zhì)量守恒方程為
(1)
其中
式中,ρ為巖體密度,kg/m3;S為巖體的儲(chǔ)水系數(shù);t為時(shí)間,a;Qm為滲流源匯項(xiàng),kg/(m3·s);v為巖塊中水流速,m/s;k為巖體滲透率,m2;μ為流體動(dòng)力黏度,Pa·s;p為壓力,Pa。
裂縫內(nèi)質(zhì)量守恒方程為
(2)
其中
式中,df為裂縫寬度,m;ρf為裂縫密度,kg/m3;Sf為裂縫儲(chǔ)水系數(shù),Pa-1;Qf為巖塊與裂隙面的流量交換;kf為裂縫滲透率,m2;τ沿裂縫切向梯度。
1.2.2 能量守恒方程
對(duì)于增強(qiáng)型地?zé)嵯到y(tǒng)中流體的換熱,工作流體與固體骨架間的溫差不容忽視,因此局部熱平衡模型將導(dǎo)致較大的計(jì)算誤差。采用兩個(gè)不同的能量方程分別對(duì)流體和固體骨架中的熱量傳遞進(jìn)行描述實(shí)現(xiàn)非平衡傳熱,從而更加真實(shí)準(zhǔn)確地模擬基質(zhì)巖體與裂縫內(nèi)流體的能量傳遞過(guò)程。
流體能量守恒方程[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分別對(duì)應(yīng)巖石和流體;ε為巖石孔隙率。
裂縫內(nèi)熱傳導(dǎo)與熱對(duì)流能量平衡方程[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+ν).
裂隙近似由一對(duì)法向和切向位移表面構(gòu)成,巖體裂隙變形方程為
(7)
(8)
式中,u為位移,m;σ為總應(yīng)力,Pa;σ′為有效應(yīng)力,Pa;K為剛度,Pa/m;下標(biāo)n、s分別指示裂隙面法向和切向。通過(guò)對(duì)裂縫內(nèi)部邊界進(jìn)行弱形式描述代表裂縫的剛度[5,10]。
(1)滲流場(chǎng)對(duì)溫度場(chǎng)的影響:裂縫內(nèi)冷流體與高溫巖體進(jìn)行熱交換,巖體中原有平衡狀態(tài)的溫度場(chǎng)發(fā)生改變。
(2)溫度場(chǎng)對(duì)滲流場(chǎng)的影響:高溫高壓巖體內(nèi),裂隙內(nèi)水滲流物理特性(如密度、黏度等)不再是常數(shù),隨溫度變化發(fā)生變化,從而影響巖體內(nèi)流體流動(dòng)。
(3)溫度場(chǎng)對(duì)應(yīng)力場(chǎng)的影響:巖體溫度的改變影響巖體固有的物理力學(xué)性質(zhì),同時(shí)溫度效應(yīng)產(chǎn)生的熱應(yīng)力會(huì)導(dǎo)致原有應(yīng)力場(chǎng)分布改變。Zhu等[12]通過(guò)實(shí)驗(yàn)獲得了不同溫度下彈性模量。
(4)應(yīng)力場(chǎng)對(duì)溫度場(chǎng)的影響:應(yīng)力場(chǎng)變化引起巖體結(jié)構(gòu)改變時(shí)的變形生熱,同時(shí)巖體的結(jié)構(gòu)變形會(huì)影響孔隙內(nèi)部的導(dǎo)熱性能。
(5)滲流場(chǎng)對(duì)應(yīng)力場(chǎng)的影響:孔隙壓力的存在影響巖石有效應(yīng)力。
(6)應(yīng)力場(chǎng)對(duì)滲流場(chǎng)的影響:應(yīng)力場(chǎng)變化導(dǎo)致了孔隙和裂隙變化,進(jìn)而改變了巖體與裂縫的滲透性能。
與裂縫的滲透率相比巖石基質(zhì)的滲透率極低,因此忽略基質(zhì)滲透率的變化。裂縫滲透率的變化進(jìn)行描述為
(9)
應(yīng)用有限元軟件COMSOL Multiphysics可將各物理場(chǎng)聯(lián)立求解,實(shí)現(xiàn)耦合分析中各物理場(chǎng)的全面耦合,從而反映真實(shí)的模擬過(guò)程。采用有限元方法對(duì)裂縫網(wǎng)絡(luò)的模擬通常使用具有一定厚度的單元層刻畫(huà)裂縫,但當(dāng)裂縫網(wǎng)絡(luò)規(guī)模較大時(shí),往往受限于計(jì)算能力而無(wú)法求解。借鑒文獻(xiàn)[10]和[13],將二維離散裂縫網(wǎng)絡(luò)設(shè)置為線單元,并借助軟件中內(nèi)置的裂隙模塊進(jìn)行求解。所建立耦合模型與借助COMSOL實(shí)現(xiàn)的數(shù)值求解方法的可行性與準(zhǔn)確性通過(guò)對(duì)比白冰[14]針對(duì)一維條件下飽和土熱彈性固結(jié)問(wèn)題獲得的解析解得到驗(yàn)證。
假定壓裂改造后含裂縫網(wǎng)絡(luò)干熱巖計(jì)算模型為700 m×700 m的二維模型。改造后干熱巖中裂縫網(wǎng)絡(luò)認(rèn)為是在注入水壓力、熱應(yīng)力與儲(chǔ)層中原有裂縫的共同作用下所形成的[1-4,10],模型中采用隨機(jī)裂縫進(jìn)行描述。隨機(jī)生成的裂縫網(wǎng)絡(luò)由兩組裂隙構(gòu)成,裂縫數(shù)量為245條,傾角分別為0°和90°,裂縫長(zhǎng)度均值為70 m,方差為30 m,由于增強(qiáng)型地?zé)嵯到y(tǒng)中的主要導(dǎo)水通道為聯(lián)通裂縫,在構(gòu)建裂縫網(wǎng)絡(luò)模型時(shí)去掉不聯(lián)通的裂縫[13]。模型采用兩注一采的注采方式,計(jì)算模型如圖1所示,藍(lán)色圓點(diǎn)代表注入井,紅色圓點(diǎn)代表采出井,注入井與采出井距離均為350 m。
圖1 含裂縫網(wǎng)絡(luò)干熱巖模型Fig.1 Hot dry rock model with fracture network
利用上述模型模擬含裂縫網(wǎng)絡(luò)干熱巖采熱全過(guò)程,耦合分析從注入水的時(shí)刻開(kāi)始,時(shí)間步長(zhǎng)為1 d,總時(shí)間段為40 a。邊界條件和初始條件如下:
(1)滲流場(chǎng)。為確保系統(tǒng)中的水循環(huán)給定壓力邊界條件。假設(shè)注入井給定壓力為44 MPa,生產(chǎn)井給定壓力為40 MPa,注采井直徑設(shè)置為0.2 m。其他外部條件為不可滲透邊界。在攜熱流體注入前,儲(chǔ)層中的初始水壓為40 MPa。
(2)溫度場(chǎng)。注入井的水溫為20 ℃。為真實(shí)模擬EGS周邊基巖對(duì)EGS區(qū)域的傳熱將邊界設(shè)為開(kāi)邊界。對(duì)于干熱巖基質(zhì)巖體初始溫度為200 ℃。
(3)應(yīng)力場(chǎng)。模型四周邊界設(shè)置為位移約束邊界。為便于明確低溫流體注入過(guò)程對(duì)儲(chǔ)層應(yīng)力場(chǎng)的擾動(dòng),不考慮地應(yīng)力場(chǎng)的影響。
應(yīng)用COMSOL Multiphysics實(shí)現(xiàn)了固體形變、流體滲流與熱量傳遞的全耦合,并對(duì)所建立的二維模型進(jìn)行求解。采用三角形網(wǎng)格進(jìn)行網(wǎng)格剖分,形成三角形單元49 725個(gè)與邊界單元5 932個(gè)。
模型中采用的主要材料參數(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),儲(chǔ)水率為1×10-8Pa-1,熱膨脹系數(shù)為6×10-6K-1;裂縫密度為1 200 kg/m3,寬度為0.002 m,滲透率為3×10-11m2,孔隙度為30%,儲(chǔ)水率為1×10-9Pa-1,比熱容為850 J/(kg·K),法向剛度為100 GPa/m,切向剛度為30 GPa/m;Biot系數(shù)為1。
向含裂縫網(wǎng)絡(luò)干熱巖儲(chǔ)層中注入低溫?cái)y熱流體循環(huán)采熱過(guò)程伴隨著滲流場(chǎng)、溫度場(chǎng)、應(yīng)力場(chǎng)的互相影響與制約。由于基巖滲透率極低,少量攜熱流體濾失到基巖中,如圖2所示。流體主要沿裂縫滲流,聯(lián)通的裂縫網(wǎng)絡(luò)成為EGS中主要的導(dǎo)水通道。同時(shí)攜熱流體在流動(dòng)過(guò)程中與裂縫周圍基巖發(fā)生熱量交換,如圖3所示。低溫流體被不斷加熱,裂縫周圍基巖中熱量被不斷提取,攜熱流體流動(dòng)一段距離后被加熱至儲(chǔ)層溫度,此時(shí)的流體與基巖達(dá)到熱平衡不再發(fā)生熱交換。不論是流體滲流產(chǎn)生的水壓作用還是注入流體與儲(chǔ)層溫差導(dǎo)致的熱應(yīng)力作用都將對(duì)應(yīng)力場(chǎng)分布形成擾動(dòng),應(yīng)力場(chǎng)的變化將直接影響儲(chǔ)層中裂縫的滲流能力。在裂縫網(wǎng)絡(luò)中設(shè)置裂縫滲透率探針,如圖4所示。低溫流體的注入使高溫基巖遇冷收縮,同時(shí)流體壓力也導(dǎo)致裂縫和基巖受力發(fā)生形變,在兩者的共同作用下裂縫滲透率提高。由于a、d兩點(diǎn)處離注入井距離近,開(kāi)發(fā)初期即受到水壓力與熱應(yīng)力的影響,裂縫滲透率升高較早,而b、c處距離注入井遠(yuǎn),水壓力與冷鋒面擴(kuò)散至b、c處需一段時(shí)間,因此系統(tǒng)開(kāi)發(fā)約4 a后滲透率才有所提升,且此處流體已被加熱,熱應(yīng)力的作用不及a、d處劇烈,導(dǎo)致提升幅度不大。
裂縫分布存在非均勻性,裂縫密集區(qū)域的熱對(duì)流更明顯,冷鋒面范圍更大,如圖3所示。裂縫稀疏區(qū)域,周圍基巖中未提取熱量多,基巖對(duì)攜熱流體加熱能力更強(qiáng),冷鋒面范圍較小。裂縫分布的非均勻性也同樣導(dǎo)致了儲(chǔ)層有效應(yīng)力場(chǎng)的非均勻性,如圖5所示。
隨著開(kāi)發(fā)的進(jìn)行,冷鋒面逐漸向采出井推進(jìn),直至該EGS系統(tǒng)出現(xiàn)熱突破,如圖3(c)所示。裂縫周圍基巖溫度的降低也使裂縫內(nèi)流體黏度增加流速下降(圖2)。冷鋒向采出井的不斷推進(jìn)也導(dǎo)致了有效應(yīng)力場(chǎng)的不斷變化。因此在EGS開(kāi)發(fā)中流體滲流、熱量傳遞、固體形變的相互作用伴隨著整個(gè)開(kāi)發(fā)過(guò)程。
圖2 不同時(shí)刻裂縫網(wǎng)絡(luò)內(nèi)速度場(chǎng)分布Fig.2 Distribution of velocity field in fracture network at different time
圖3 不同時(shí)刻干熱巖儲(chǔ)層溫度場(chǎng)分布Fig.3 Distribution of temperature field in hot reservoir at different time
圖4 裂縫滲透率探針?lè)植寂c不同時(shí)刻滲透率演化Fig.4 Location of fracture permeability probes and permeability evolution
圖5 不同時(shí)刻干熱巖儲(chǔ)層有效應(yīng)力場(chǎng)等值線圖Fig.5 Contour map of effective stress field in hot reservoir at different time
注采井網(wǎng)的變化必然對(duì)滲流場(chǎng)、溫度場(chǎng)、應(yīng)力場(chǎng)產(chǎn)生影響,最終改變干熱巖儲(chǔ)層的采熱性能。采用圖1中裂縫網(wǎng)絡(luò)模型,分別設(shè)置一注一采、兩注一采、四注一采3種不同的注采井網(wǎng),注采井井距均設(shè)置為350 m。以采出井的產(chǎn)出液溫度、產(chǎn)出液質(zhì)量流量、采熱速率[15-17]作為采熱性能評(píng)價(jià)指標(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)所示,一注一采開(kāi)發(fā)5 a時(shí)產(chǎn)出液質(zhì)量流量為32 kg/s,井網(wǎng)改為兩注一采后開(kāi)發(fā)5 a時(shí)產(chǎn)出液質(zhì)量流量為55 kg/s,四注一采開(kāi)發(fā)5 a時(shí)達(dá)到了75 kg/s。不同注采井網(wǎng)干熱巖儲(chǔ)層溫度場(chǎng)分布如圖7所示。注入井的加入,使干熱巖中更多區(qū)域的地?zé)豳Y源得到動(dòng)用,同時(shí)合理的注采井布局并未造成冷鋒面擴(kuò)大導(dǎo)致的熱突破加快。如圖8(b)所示,注入井的加入反而緩解了產(chǎn)出液溫度的降低。分析圖8可見(jiàn),裂縫網(wǎng)絡(luò)中同一位置處,采用一注一采時(shí)裂縫內(nèi)流速最快,增加注入井?dāng)?shù)量后,在采出井附近產(chǎn)生流體累加,裂縫內(nèi)流速減緩,增加了其在干熱巖中的換熱從而緩解了產(chǎn)出液溫度的降低。注入井?dāng)?shù)量的增加使產(chǎn)出液質(zhì)量流量大幅提高,同時(shí)產(chǎn)出液溫度未出現(xiàn)下降,因此EGS采熱速率提升明顯,四注一采開(kāi)發(fā)5 a時(shí)采熱速率達(dá)到60 MW??梢?jiàn)在井網(wǎng)中加入注入井對(duì)于系統(tǒng)采熱速率提高有效。
圖7 開(kāi)發(fā)20 a時(shí)不同注采井網(wǎng)干熱巖儲(chǔ)層溫度場(chǎng)分布Fig.7 Distribution of temperature field in hot reservoir with different inject-product well pattern after 20-year development
圖8 開(kāi)發(fā)20 a時(shí)不同注采井網(wǎng)干熱巖儲(chǔ)層速度場(chǎng)分布Fig.8 Distribution of velocity field in fracture network with different inject-product well pattern after 20-year development
采熱速率結(jié)合了產(chǎn)出液質(zhì)量流量與產(chǎn)出液溫度,是評(píng)價(jià)系統(tǒng)采熱性能的綜合評(píng)價(jià)指標(biāo)。鑒于對(duì)注采井網(wǎng)研究的認(rèn)識(shí),采用四注一采的開(kāi)發(fā)方式,針對(duì)模型關(guān)鍵參數(shù)開(kāi)展敏感性分析。
圖9 儲(chǔ)層物性對(duì)采熱速率的影響Fig.9 Effect of reservoir properties on heat mining rate
基質(zhì)滲透率的增加首先引起滲流場(chǎng)的改變,更多低溫?cái)y熱流體由裂縫流向基質(zhì),高溫基巖遇冷收縮加劇。因此裂縫開(kāi)度增加,應(yīng)力場(chǎng)受到影響,同時(shí)裂縫滲透率提高,產(chǎn)出液質(zhì)量流量增加采熱速率提高,如圖9(a)所示,基質(zhì)滲透率為1×10-17m2與1×10-18m2時(shí)相比采熱速率增幅較大,但與1×10-16m2時(shí)相比增幅較小,這是由于基巖滲透率增大到一定程度后,大量工作流體濾失所致。
基質(zhì)導(dǎo)熱系數(shù)通過(guò)改變溫度場(chǎng)從而影響滲流場(chǎng)與應(yīng)力場(chǎng)。基質(zhì)導(dǎo)熱系數(shù)對(duì)儲(chǔ)層溫度場(chǎng)影響較小,如圖9(b)所示。其對(duì)采熱速率的影響也較微弱,該模擬結(jié)果與Sun等[10]的研究結(jié)果一致。
基質(zhì)熱膨脹系數(shù)的不同首先影響應(yīng)力場(chǎng),熱膨脹系數(shù)越大,高溫基巖遇冷產(chǎn)生的熱應(yīng)力越強(qiáng),基巖收縮位移越大,裂縫滲透率提升越大,如圖9(c)所示。熱膨脹系數(shù)增加開(kāi)發(fā)0~30 a間采熱速率增加,但30 a后采熱速率降低。這是由于低溫流體注入產(chǎn)生的熱應(yīng)力作用在使裂縫滲透率提高的同時(shí)也導(dǎo)致了熱竄流[7],使熱突破更快,產(chǎn)出液溫度降低速度增加。
注采壓差增加儲(chǔ)層中流體流速增加,采熱速率提高,同時(shí)流速的增加將導(dǎo)致攜熱流體換熱不充分,如圖10所示,隨著開(kāi)發(fā)時(shí)間的延長(zhǎng)由注采壓差增加引起的采熱速率提升幅度減小。若壓差過(guò)大將導(dǎo)致EGS開(kāi)發(fā)壽命縮短,采熱速率快速下降。不同注入溫度對(duì)應(yīng)力場(chǎng)與滲流場(chǎng)均有影響。注入溫度越高低溫流體導(dǎo)致的熱應(yīng)力作用越小,高流體溫度也更有利于流體流動(dòng),因此注入流體溫度越高采熱速率越高。
(1)向含裂縫網(wǎng)絡(luò)干熱巖中注入攜熱流體采熱過(guò)程中涉及THM耦合,其中應(yīng)力場(chǎng)受水壓力與溫差所致熱應(yīng)力的共同作用,裂縫網(wǎng)絡(luò)分布的非均勻性對(duì)流體滲流、冷鋒面擴(kuò)散及應(yīng)力場(chǎng)分布均有影響。
(2)采用局部熱非平衡理論對(duì)干熱巖采熱過(guò)程進(jìn)行模擬,考慮換熱過(guò)程中工作流體與固體骨架間的溫度差異可更真實(shí)地描述巖石與流體間熱交換。
(3)井網(wǎng)中注入井的加入,使干熱巖中更多區(qū)域的地?zé)豳Y源得到動(dòng)用,產(chǎn)出液質(zhì)量流量大幅提升,合理的布局可避免熱突破加快,從而有效提高系統(tǒng)采熱速率,四注一采開(kāi)發(fā)5 a時(shí)采熱速率可達(dá)60 MW。
(4)儲(chǔ)層基質(zhì)滲透率對(duì)采熱速率的影響大于基質(zhì)導(dǎo)熱系數(shù)?;|(zhì)熱膨脹系數(shù)通過(guò)控制溫差引起的熱應(yīng)力影響采熱性能,其數(shù)值越大攜熱流體越容易發(fā)生熱竄流。注入溫度越高對(duì)干熱巖開(kāi)發(fā)越有利,注采壓差應(yīng)調(diào)節(jié)適當(dāng),過(guò)大的壓差雖然產(chǎn)出液質(zhì)量流量高,但會(huì)影響開(kāi)發(fā)壽命。