樊冬艷, 孫 海, 姚 軍, 嚴(yán) 俠, 張 凱, 張建光, 張 林

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

能源和環(huán)境是可持續(xù)發(fā)展戰(zhàn)略的核心與關(guān)鍵,常規(guī)的化石能源已十分有限且不可再生,同時(shí)隨著全球氣候變暖,節(jié)能減排的呼聲日益高漲,開發(fā)利用清潔的、可再生的新能源已成為各國(guó)21世紀(jì)的重要戰(zhàn)略選擇。地?zé)崮茏鳛橐环N經(jīng)濟(jì)、可再生、環(huán)保的新型能源,在國(guó)際能源結(jié)構(gòu)中所占比重越來越大[1]。中國(guó)地質(zhì)科學(xué)院和中國(guó)科學(xué)院分別進(jìn)行了中國(guó)干熱巖資源量的估算,分別為2.5×1025J[2]和2.1×1025J[3],按照2%的回收率計(jì)算,為2010年能源消耗的4 400倍,因此中國(guó)對(duì)于干熱巖地?zé)衢_采和利用給予了極大重視[4-6]。干熱巖儲(chǔ)層的開采一般將低溫水通過注入井注入熱儲(chǔ)層(人工裂縫和天然裂縫),經(jīng)加熱變成熱水,再?gòu)纳a(chǎn)井排出,實(shí)現(xiàn)地?zé)岬睦肹7]。目前,研究干熱巖裂縫內(nèi)流動(dòng)模型主要包括單裂縫模型[8]、連續(xù)介質(zhì)模型[9-10]以及離散裂縫網(wǎng)絡(luò)模型[11-13]等。其中單裂縫模型[8]假設(shè)水通過一個(gè)裂縫在注水井和生產(chǎn)井之間流動(dòng),實(shí)際儲(chǔ)層存在大量裂縫因此與實(shí)際不太相符。連續(xù)介質(zhì)模型[9-10]假設(shè)水通過均勻多孔介質(zhì)在注水井和生產(chǎn)井之間流動(dòng),無(wú)法準(zhǔn)確描述裂縫性質(zhì),因此存在一定的局限性。離散裂縫網(wǎng)絡(luò)模型[11-13]假定儲(chǔ)層由大量裂縫網(wǎng)絡(luò)相互構(gòu)成,更加符合實(shí)際情況。離散裂縫網(wǎng)絡(luò)滲流模型的研究始于20世紀(jì)60年代對(duì)巖石水力學(xué)的研究,因此關(guān)于流動(dòng)及流固耦合方面的研究已經(jīng)非常完善[14-15]。目前地?zé)岫鄨?chǎng)耦合模型主要基于數(shù)值模擬方法,采用商業(yè)軟件如TOUGH2進(jìn)行數(shù)值模擬[16-18],而解析方法較少[19-20],主要以單條裂縫模型為主。筆者在考慮流體黏度及裂縫開度變化的情況下,建立基于離散裂縫網(wǎng)絡(luò)模型干熱巖熱流固耦合開采數(shù)學(xué)模型,采用Laplace變換得到模型的解析解,以期為干熱巖儲(chǔ)層的開采提供技術(shù)指導(dǎo)及理論支持。

1 干熱巖儲(chǔ)層的離散裂縫網(wǎng)絡(luò)模型

1.1 物理模型

由于干熱巖儲(chǔ)層地?zé)豳Y源的開采系統(tǒng)一般由熱儲(chǔ)層、注水井、生產(chǎn)井和壓裂裂縫以及天然裂縫組成。冷水通過注水井注入熱儲(chǔ)層,基于熱傳導(dǎo)和對(duì)流等機(jī)制加熱變成熱水,再?gòu)纳a(chǎn)井排出,實(shí)現(xiàn)地?zé)豳Y源的開發(fā)利用。為了準(zhǔn)確模擬地?zé)醿?chǔ)層的熱、流、固等復(fù)雜機(jī)制,針對(duì)干熱巖儲(chǔ)層采用離散裂縫網(wǎng)絡(luò)模型(圖1),其中采用一維的線單元代表二維的矩形裂縫,流體通過裂縫網(wǎng)絡(luò)流動(dòng),基巖中無(wú)流體流動(dòng)只進(jìn)行熱傳導(dǎo)過程;離散裂縫內(nèi)的流動(dòng)為穩(wěn)定滲流,應(yīng)力通過改變裂縫的開度對(duì)流動(dòng)產(chǎn)生影響,裂縫內(nèi)主要以熱對(duì)流為主;其次離散裂縫網(wǎng)絡(luò)模型主要由注水井、生產(chǎn)井、邊、結(jié)點(diǎn)構(gòu)成,考慮了單井注入單井開采情況,只有相互溝通的裂縫參與流體流動(dòng),只與一個(gè)結(jié)點(diǎn)相連的裂縫無(wú)流體流動(dòng)。

圖1 干熱巖儲(chǔ)層離散裂縫網(wǎng)絡(luò)模型Fig.1 Discrete fracture network model of dry-hot rock reservoir

1.2 離散裂縫模型的流動(dòng)模型

為了研究裂縫中水的流動(dòng),基于質(zhì)量守恒方程以任意一個(gè)結(jié)點(diǎn)n作為研究對(duì)象(圖2),則流入與流出的質(zhì)量應(yīng)相等,當(dāng)該結(jié)點(diǎn)為注水井時(shí),定注入量生產(chǎn),因此水的質(zhì)量變化,即注水井的注入量minj已知,當(dāng)結(jié)點(diǎn)為生產(chǎn)井時(shí),定生產(chǎn)壓力pw,該結(jié)點(diǎn)的壓力已知。因此對(duì)于結(jié)點(diǎn)n可建立質(zhì)量守恒方程為

(1)

其中

裂縫中的流動(dòng)服從Poiseuille 流動(dòng)規(guī)律,因此

(2)

μw(T)=ATB.

(3)

其中

A=109.574-8.405 64S+0.313 34S2+8.722 13×

10-3S3,

B=-1.121 66+2.639 51×10-2S-6.794 61×10-4S2-5.471 19×10-5S3+1.555 86×10-6S4.

式中,S為水中礦物質(zhì)質(zhì)量分?jǐn)?shù),%;T為溫度,1 ℉=[(T-32)×5/9]℃。

1.3 離散裂縫模型的開度模型

大量的實(shí)驗(yàn)和礦場(chǎng)實(shí)踐表明裂縫開度在開采過程中具有動(dòng)態(tài)性質(zhì)[23-25],例如裂縫的開度受孔隙彈性和熱彈性的影響,目前大部分模擬都沒有考慮這些動(dòng)態(tài)變化,筆者基于Ahmad Ghassemi 模型[26]在不考慮濾失的情況下給出裂縫開度在不同時(shí)間和位置隨巖石特性和溫度的變化關(guān)系為

(4)

其中

χ=αT(1+v)/(1-v),TΔ=T(x,z,0)-T(0,0,t).

式中,w為裂縫開度,m;w0為初始裂縫開度;z為沿裂縫方向的長(zhǎng)度坐標(biāo),m;t為時(shí)間,s;erfc(x)為互補(bǔ)誤差函數(shù);χ為熱彈性應(yīng)力系數(shù);αT為基巖的線性熱膨脹系數(shù),℃-1;v為排水泊松比;TΔ為初始時(shí)刻溫度與端點(diǎn)的溫度差,℃;Kr為巖石導(dǎo)熱系數(shù),W/(m·℃);ρr和ρw分別為基巖和流體水的密度,kg/(m3);Cr和Cw分別為基巖和流體水的比熱容,J/(kg·℃);q0為注入流體速度,m2/s。

1.4 離散裂縫模型的溫度模型

為了求解離散裂縫網(wǎng)絡(luò)模型的流動(dòng),取任意一條裂縫i作為研究對(duì)象,由于假設(shè)流體只在裂縫中流動(dòng),沿著裂縫建立相應(yīng)的坐標(biāo),如圖2所示。

建立基巖中一維熱傳導(dǎo)數(shù)學(xué)模型為

(5)

其中

Ψ=(T-Tw)/(Tr-Tw)i.

式中,Tr為初始儲(chǔ)層溫度,℃;Tw為注入水的溫度,℃;ar為基巖的熱擴(kuò)散系數(shù),ar=Kr/(ρrCr),故無(wú)因次溫度的范圍為[0,1]。

圖2 單條裂縫坐標(biāo)示意圖Fig.2 Schematic diagram of single fracture

其次,裂縫中對(duì)流的流體在邊界處,通過熱傳導(dǎo)的方式獲取基巖中熱量,忽略巖石與水界面阻力,假定水的溫度與裂縫表面溫度相等,即

(6)

式中,wi為第i條裂縫的開度,m;vi為第i條裂縫的流速,m/s。

2 地?zé)醿?chǔ)層離散裂縫網(wǎng)絡(luò)模型的求解

2.1 離散裂縫網(wǎng)絡(luò)模型溫度場(chǎng)求解

將一維熱傳導(dǎo)方程(5)、(6)變換到Laplace空間,則方程化為

(7)

常微分方程的通解為

(8)

代入初始和邊界條件方程(7)可得單條裂縫的溫度分布為

(9)

其中

當(dāng)裂縫為離散裂縫網(wǎng)絡(luò)時(shí),第i條裂縫的初始溫度由上一級(jí)裂縫j的溫度決定,故

代入到方程(9)得:

(10)

(11)

其中sk表示到達(dá)裂縫i的所有路徑。

圖3 注入井與生產(chǎn)井的多級(jí)裂縫網(wǎng)絡(luò)示意圖Fig.3 Schematic diagram for a multi-stage fracture network between injection wells and production wells

最后,將Laplace空間溫度分布進(jìn)行Laplace反演,得到實(shí)空間下裂縫網(wǎng)絡(luò)中,任意一條裂縫的溫度分布為

(12)

其中

2.2 離散裂縫網(wǎng)絡(luò)模型耦合求解過程

為了說明離散裂縫網(wǎng)絡(luò)模型的耦合求解過程,結(jié)合簡(jiǎn)單縫網(wǎng)進(jìn)行說明,如圖4所示。地?zé)醿?chǔ)層有一口注水井和一口生產(chǎn)井,注水井定注入量,生產(chǎn)井定生產(chǎn)壓力,注入井與生產(chǎn)井通過6條裂縫、5個(gè)結(jié)點(diǎn)相連,其中裂縫1、2、5、6為水平縫,裂縫3、4為垂直縫,所有裂縫與地?zé)醿?chǔ)層垂直,簡(jiǎn)化為二維處理,把離散裂縫網(wǎng)絡(luò)模型耦合求解過程分為以下3步。

(1)求解各邊的流量分布?；陔x散裂縫網(wǎng)絡(luò)的流動(dòng)模型,對(duì)質(zhì)量守恒方程進(jìn)行求解,其中黏度采用地層溫度下流體黏度,開度取初始裂縫開度;計(jì)算得到各結(jié)點(diǎn)的壓力值,通過壓差計(jì)算各裂縫段的流量分布。如圖4所示,包括注水井和生產(chǎn)井共有6個(gè)結(jié)點(diǎn),通過方程(1)建立方程組DP=F,D為6×6的系數(shù)矩陣,P為各結(jié)點(diǎn)的壓力值,F為流入流出量。

圖4 離散裂縫網(wǎng)絡(luò)實(shí)例1Fig.4 Discrete fracture network example 1

(2)求解各條邊的溫度分布。如果某條裂縫只和單個(gè)結(jié)點(diǎn)相連,則認(rèn)為該裂縫中流體不流動(dòng),裂縫兩個(gè)端點(diǎn)無(wú)壓差;分析流入某裂縫的網(wǎng)絡(luò)路徑,如圖4從注水井到生產(chǎn)井可能流動(dòng)的路徑為s1={1,3,5},s2={1,2,4},存在兩條路徑,代入式(12),則裂縫6的溫度分布可表示為

(13)

(3)黏度及開度的耦合迭代求解。在上一步中得到各裂縫段的溫度分布,通過式(3)得到相應(yīng)水的黏度分布,代入式(4)計(jì)算開度變化,利用平均黏度和開度循環(huán)進(jìn)行步驟(1)、(2)計(jì)算溫度分布,直至平均黏度與開度的值穩(wěn)定。

3 模型正確性驗(yàn)證及參數(shù)敏感性分析

為了驗(yàn)證本模型解析方法的正確性,在不考慮黏度及開度變化時(shí),與TOUGH2數(shù)值計(jì)算結(jié)果進(jìn)行對(duì)比,TOUGH2軟件是由Berkeley 國(guó)家實(shí)驗(yàn)室研發(fā)的一款通用數(shù)值模擬器,針對(duì)裂縫及多孔介質(zhì)中的多維、多相、多組分非等溫流動(dòng)過程,對(duì)離散區(qū)域采用有限差分,時(shí)間采用一階向下全隱式方法[27],分析黏度及開度的變化對(duì)產(chǎn)水量的影響,以及熱提取過程注入量和注入溫度的影響。根據(jù)上述求解過程,對(duì)實(shí)例1(圖4)進(jìn)行求解,基本參數(shù):裂縫1、6長(zhǎng)度為100 m,裂縫2、3、4、5長(zhǎng)度為500 m,裂縫初始開度為0.003 m,裂縫高度為1 m;水密度為980 kg/m3,水的比熱容為4 300 J/(kg·℃);巖石礦化度為2 000 mg/L,巖石比熱容為1 050 J/(kg·℃),巖石密度為2 700 kg/m3,巖石導(dǎo)熱系數(shù)為2.7 W/(m·℃),泊松比為0.185;基巖的線性熱膨脹系數(shù)為8×10-6℃-1。邊界條件:注入量為0.1 kg/s;生產(chǎn)井定壓力為10 MPa;注入水溫度為10 ℃;儲(chǔ)層初始溫度為70 ℃。

3.1 模型的正確性驗(yàn)證

基于離散裂縫網(wǎng)絡(luò)模型,在不考慮水的黏度及裂縫開度變化的情況下,計(jì)算實(shí)例1,將裂縫網(wǎng)絡(luò)中各裂縫的溫度分布與TOUGH2數(shù)值計(jì)算結(jié)果進(jìn)行對(duì)比,如圖5所示。本模型解析解結(jié)果與數(shù)值解TOUGH2計(jì)算結(jié)果吻合,生產(chǎn)井出口端早期誤差約為3%,但隨著時(shí)間增加誤差減小,結(jié)果幾乎一致,說明本解析模型的可行性及正確性,當(dāng)然隨著裂縫復(fù)雜程度的增加,誤差稍有變化;其次各裂縫出口端的溫度變化,隨著時(shí)間的增加,出口端的溫度逐漸降低,其中注水井裂縫1由于冷水的不斷注入,溫度下降最快,其次為裂縫2、3出口端,裂縫4、5出口端,裂縫6出口端即生產(chǎn)井處。

圖6為實(shí)例1離散裂縫網(wǎng)絡(luò)各點(diǎn)處的溫度分布,包括生產(chǎn)1 a和生產(chǎn)10 a的兩種情況。由圖6可見,生產(chǎn)1 a時(shí),生產(chǎn)井的溫度還沒有突破,產(chǎn)出井水的溫度為70 ℃,而到了生產(chǎn)10 a整個(gè)裂縫網(wǎng)絡(luò)系統(tǒng)的溫度下降,由圖6可見生產(chǎn)井熱突破時(shí)間大約為3 a。

圖5 解析模型與TOUGH2結(jié)果對(duì)比Fig.5 Comparison between analytic model and TOUGH2 results

圖6 離散裂縫網(wǎng)絡(luò)溫度分布Fig.6 Temperature distribution of discrete fracture network

3.2 離散裂縫網(wǎng)絡(luò)黏度及開度的影響

3.2.1 離散裂縫網(wǎng)絡(luò)黏度

在地?zé)醿?chǔ)層隨著冷水的注入,水的黏度隨溫度和壓力的變化發(fā)生改變,特別是溫度對(duì)黏度的影響較大,本文中假定水的黏度滿足McCain公式,主要受礦化度及溫度的影響。當(dāng)不考慮水的黏度變化時(shí),水的黏度通過注入水初始溫度由式(3)計(jì)算得到,在實(shí)例1的基礎(chǔ)上不考慮開度變化,計(jì)算各條裂縫的平均黏度,如圖7所示。

由圖7可見,隨著時(shí)間的增加,各條裂縫內(nèi)水的黏度增加,滲流速度減小,故出口端流量減小,趨于不考慮黏度變化時(shí)流量值;當(dāng)不考慮黏度變化時(shí),黏度由注入水初始溫度(10 ℃)計(jì)算得到,黏度較大為1.368 mPa·s,隨著冷水的注入,沿著離散裂縫分布,裂縫1最靠近注入端,因此裂縫1黏度與不考慮黏度變化時(shí)幾乎重合,裂縫2、3與裂縫4、5位置對(duì)稱,因此計(jì)算黏度相同,同時(shí)說明算法的正確性,沿著冷水由注入端至出口端,裂縫中水的溫度逐漸上升,因此沿著水的流動(dòng)方向黏度逐漸減小,但隨著時(shí)間的增加,裂縫中水的溫度逐漸降低,黏度逐漸增大,增加幅度越來越小,特別是早期黏度增加比較明顯,晚期幾乎不變,對(duì)于出口端流量在開采8 a后幾乎沒有太大的影響。

3.2.2 離散裂縫網(wǎng)絡(luò)開度

地?zé)醿?chǔ)層基于離散裂縫模型,考慮孔隙壓力及熱彈性的影響,建立開度隨溫度壓力的變化過程。在不考慮黏度變化的基礎(chǔ)上,計(jì)算實(shí)例1得到各裂縫的開度變化,圖8(a)為各裂縫在模擬20 a的過程中的平均開度隨時(shí)間變化,圖8(b)為考慮與不考慮裂縫開度的產(chǎn)量變化。

圖7 考慮與不考慮黏度變化時(shí)曲線對(duì)比Fig.7 Comparison curves of considered viscosity change or not

圖8 考慮與不考慮開度變化時(shí)曲線對(duì)比Fig.8 Comparison curves of considered fracture aperture change or not

由圖8可見,隨著時(shí)間的增加,各條裂縫的開度逐漸增大,同時(shí)出口端流量增加迅速,主要是由于裂縫的開度對(duì)流量的影響最大,裂縫開度的3次方為滲透率的函數(shù),而流量直接取決于滲透率;不考慮開度變化時(shí),設(shè)置開度為0.003 m,出口端流量較小;由于冷水的注入,受應(yīng)力和熱彈性的影響導(dǎo)致裂縫開度的增加,且沿著水流動(dòng)的方向裂縫開度的增加幅度越來越小,即在注入端裂縫開度的增加最大,生產(chǎn)端裂縫開度變化最小。

3.2.3 同時(shí)考慮黏度及開度變化

通過以上分析可見,考慮黏度后出口端的流量變化不大,而考慮開度變化后出口端的流量增加明顯,在此同時(shí)考慮黏度和開度變化時(shí),出口端流量的變化規(guī)律。在實(shí)例1的基礎(chǔ)上同時(shí)考慮黏度隨溫度的變化,開度隨溫度和地層性質(zhì)的影響得到圖9。

由圖9可見,在不考慮黏度和開度變化時(shí),出口端的流量為定值;考慮黏度變化后,早期流量變化明顯,逐漸減小并趨于不考慮黏度變化時(shí)流量;考慮開度變化時(shí)流量迅速增加,主要因?yàn)榱芽p的開度與流量成3次方關(guān)系,當(dāng)同時(shí)考慮裂縫的開度變化和水的黏度變化時(shí),產(chǎn)量依舊呈上升趨勢(shì),由圖9可見相對(duì)而言裂縫開度對(duì)出口端流量的影響更加明顯,因此在模擬過程中有必要考慮裂縫開度對(duì)地?zé)醿?chǔ)層的影響。

圖9 對(duì)比不同因素出口端流量Fig.9 Export flow rates with different factors

3.3 注入量的影響

為了分析注入量對(duì)離散裂縫網(wǎng)絡(luò)溫度的影響,在其他參數(shù)不變的情況下,分別設(shè)置注入量為0.1、0.2和0.3 kg/s時(shí),得到裂縫1出口端和生產(chǎn)井出口端處的溫度變化,如圖10所示。由圖10可見,隨著注入量的增加,裂縫1出口端和生產(chǎn)井出口端中水的溫度逐漸減小,且隨著注入量的線性增加,溫度降低幅度越來越小,主要由于注入量增加,流體從地層帶走的熱量越多,裂縫中熱量通過基巖塊熱傳導(dǎo)越難及時(shí)補(bǔ)充。裂縫1出口端早期溫度變化迅速,后期變緩;同時(shí)隨著注入量的增加,生產(chǎn)井出口端熱突破時(shí)間越來越早,因此在考慮熱提取的過程中,不僅要考慮生產(chǎn)井流出水量,還需要綜合考慮熱水的溫度及其他參數(shù)的影響,在各參數(shù)優(yōu)化的情況下,選擇合適的注入量,使產(chǎn)出熱值實(shí)現(xiàn)最大化。

圖10 不同注入量下生產(chǎn)井出口端溫度變化Fig.10 Variation of temperature with injection volumeat production well outlet

3.4 注入溫度的影響

為了分析注入水的溫度對(duì)離散裂縫網(wǎng)絡(luò)溫度的影響,在其他參數(shù)不變的情況下,得到裂縫1出口端和生產(chǎn)井出口端處的溫度變化,如圖11所示。

圖11 不同注入溫度下生產(chǎn)井出口端溫度變化Fig.11 Variation of temperature with injection temperature at production well outlet

由圖11可見,隨著注入水溫度的增加,生產(chǎn)井出口端處溫度早期相同,到后期注入溫度越低,出口端生產(chǎn)井溫度下降越快,溫度越低。在注入水量和裂縫網(wǎng)絡(luò)結(jié)構(gòu)相同的情況下,熱損失近似，且由圖11可見在不同注入溫度下熱突破時(shí)間幾乎相同,受注入水溫的影響不大,而主要受地層及裂縫網(wǎng)絡(luò)的影響;其次隨著注入水溫度的線性增加,裂縫1出口端的溫度早期迅速降低,并趨于相應(yīng)的注入溫度。

4 結(jié) 論

(1)隨著時(shí)間的增加,各裂縫內(nèi)水的黏度增加,同時(shí)出口端流量減小;當(dāng)冷水由注入端流經(jīng)采出端采出,裂縫中水的溫度逐漸上升,沿著水的流動(dòng)方向黏度逐漸減小,但隨著時(shí)間的增加,裂縫中水的溫度逐漸降低,黏度逐漸增大,增加幅度減小;黏度對(duì)流量的影響主要在早期,晚期影響較小。

(2)隨著時(shí)間的增加,各條裂縫的開度逐漸增大,同時(shí)出口端流量增加迅速,主要是由于裂縫的開度對(duì)流量的影響較大;在同時(shí)考慮裂縫的開度變化和水的黏度變化時(shí),產(chǎn)量依舊呈明顯的上升趨勢(shì),故在干熱巖模擬過程中不應(yīng)忽略開度的變化。

(3)隨著時(shí)間的增加,出口端的溫度逐漸降低;隨著注入量的增加,出口端中水的溫度逐漸減小,且隨著注入量的線性增加,生產(chǎn)井出口端中溫度降低幅度越來越小;隨著注入水溫度的增加,出口端生產(chǎn)井處的溫度早期相同,到后期注入溫度越低,出口端生產(chǎn)井溫度下降越快,溫度越低,且在不同注入溫度下熱突破時(shí)間幾乎相同,受注入水溫度的影響不大。