賈艷雨,王玉星,李亞輝,張海雄,王 茜
(1.中石化新星(北京)新能源研究院有限公司 鄭州分公司,河南 鄭州 450000; 2.河南省地質(zhì)研究院,河南 鄭州 450000; 3.地下清潔能源勘查開發(fā)產(chǎn)業(yè)技術(shù)創(chuàng)新戰(zhàn)略聯(lián)盟,河南 鄭州 450000)
我國目前正處在經(jīng)濟結(jié)構(gòu)轉(zhuǎn)型階段,亟需加快綠色能源的開發(fā)利用與提高當(dāng)前能源資源開發(fā)效率[1]。地?zé)崮芨鶕?jù)埋藏深度分為3類:200 m以淺的為淺層地?zé)崮?200~3 000 m的為中深層地?zé)崮?3 000 m以下的為深層地?zé)崮躘2-3]。其中,中深層地?zé)崮苤兴疅嵝偷責(zé)崮苁俏覈壳伴_發(fā)利用最為廣泛的地?zé)崮茴愋汀D壳?我國正大力推進以地?zé)豳Y源代替?zhèn)鹘y(tǒng)資源為建筑物供暖的資源利用方式以減輕資源供給壓力,超過20個省市區(qū)應(yīng)用水熱型地?zé)豳Y源進行供暖,地?zé)峁┡偯娣e位居世界第一。水熱型地?zé)崮苤械膫鳠峤橘|(zhì)為地?zé)崃黧w,具有埋藏深,天然補給緩慢等特點,屬于相對不可再生資源[4]。隨著水熱型地?zé)崮艿拇笠?guī)模開發(fā)利用,地?zé)崃黧w補采失衡,造成了儲層壓力下降、水化學(xué)污染及熱污染等一系列地質(zhì)環(huán)境問題。而回灌是解決上述問題,保證地?zé)豳Y源可持續(xù)開發(fā)利用的最有效途徑[5]。但在近幾年的回灌研究中發(fā)現(xiàn),地?zé)嵛菜毓嘤謳砹诵聠栴},即回灌井周邊熱儲溫度降低,在下一個供暖季之前遠不能恢復(fù)到最初的熱儲溫度[6]?;毓嘁矔股a(chǎn)井周邊溫度發(fā)生下降,這種現(xiàn)象稱為熱突破[7]。熱突破時間直接決定了地?zé)嵯到y(tǒng)的使用壽命,一旦發(fā)生了熱突破,熱儲的產(chǎn)熱能力將逐漸下降甚至導(dǎo)致地?zé)嵯到y(tǒng)完全停運,因此地?zé)峄毓嘌芯康闹饕康氖穷A(yù)測開采井的熱突破時間,并通過各種手段來降低開采井的熱突破風(fēng)險,從而延長地?zé)嵯到y(tǒng)的使用壽命,確保地?zé)嵯到y(tǒng)高效經(jīng)濟地運行。 除現(xiàn)場回灌試驗與在運營期間對井內(nèi)溫度進行監(jiān)測等現(xiàn)場手段外,使用數(shù)值模擬建立可反映熱儲層溫度與流場隨時間的演化規(guī)律對預(yù)測生產(chǎn)井熱突破的時間非常重要[8-10]。
許勇[11]通過數(shù)值模擬對采灌對井之間的水力聯(lián)系,合理井間距、回灌過程中化學(xué)場、滲流場及溫度場的演化過程進行模擬,預(yù)測多年后不同采灌方案下地下熱水水位的動態(tài)變化規(guī)律影響熱突破的因素有地?zé)峋髁?、采灌井間距、熱儲層巖性、熱儲層地質(zhì)構(gòu)造等因素。朱家玲等[12]以天津濱海新區(qū)塘沽地區(qū)館陶組孔隙性地層為研究重點,利用TOUGH2軟件擬合研究了地?zé)峄毓嗑g距壓差補償對回灌效率的影響。趙志宏等[13]應(yīng)用等效滲流通道模型理論框架進行了示蹤試驗反演及開采井熱突破預(yù)測。魏凱等[14]建立了含裂縫地?zé)醿拥臐B流—傳熱弱耦合模型,分析了裂縫的傾角、長度、寬度等特征參數(shù)對地?zé)醿訚B流場、溫度場的影響。李靜巖等[15]研究了熱儲上下巖層熱補償作用對CO2羽流地?zé)嵯到y(tǒng)性能的影響。崔翰博等[16]結(jié)合流固耦合傳熱理論并運用Comsol軟件,建立了離散型裂隙巖體流體傳熱模型。Lederer C等[17]利用傳熱傳質(zhì)有限元(FEHM)的模擬器進行了水—熱—力學(xué)(THM)耦合模擬,并通過單一裂縫將注入井和生產(chǎn)井連接起來,研究發(fā)現(xiàn),THM 效應(yīng)引導(dǎo)注入井和生產(chǎn)井之間的流動,導(dǎo)致溫度下降更快,并減少了熱能生產(chǎn)。Vallier B等[18]構(gòu)建了一套模擬水—熱—力學(xué)的耦合模型來模擬法國 Rittershoffen 深層地?zé)釋到y(tǒng) GRT-1 和 GRT-2 的水熱運移機制與力學(xué)響應(yīng)機制。Wang Shihao等[19]對裂隙儲層中水、熱、力學(xué)耦合過程中的半解析解及其在增強型地?zé)嵯到y(tǒng)中的應(yīng)用進行了研究。
由于中深層地?zé)崧裆钶^大,熱儲層地質(zhì)情況復(fù)雜,精確構(gòu)建與實際情形相同的物理模型非常困難。在構(gòu)建大尺度地?zé)釘?shù)值模型時,當(dāng)前的地?zé)崮P涂煞譃?類:單孔隙率模型,雙孔隙率模型和多孔隙率模型[16]。其中多孔隙率模型考慮了實際地層中的巖性幾何分布特征與裂隙的空間分布,具有較高的準確性,但往往由于需要大量現(xiàn)場勘測資料與試驗數(shù)據(jù)支持,造成構(gòu)建成本較高,當(dāng)前未見廣泛應(yīng)用。
本文依托河南蘭考縣館陶組砂巖熱儲資料,基于簡化的多孔隙率模型,采用COMSOL有限元多場耦合數(shù)值模擬軟件,在考慮基巖和蓋層熱補償與滲流對地?zé)醿觽鳠徇^程的影響的情況下,建立了含斷層地?zé)醿拥臐B流—傳熱耦合模型。采用周期“一采一灌”對井注采模式,對不同井間距與不同滲透率地層的地?zé)峋L期運行過程進行模擬,通過對熱儲層采熱運行期間溫度場與流場變化進行分析與對比,研究長時間尺度下回灌對采熱過程的熱突破現(xiàn)象發(fā)生規(guī)律。
蘭考縣隸屬于河南省開封市,地理位置處于開封市城區(qū)東部;地質(zhì)構(gòu)造上處于南華北盆地開封凹陷。開封凹陷構(gòu)造上處于南華北盆地與渤海灣盆地交叉疊合部位,北東向構(gòu)造帶的交叉疊合部位,北部緊鄰渤海灣盆地東濮凹陷,南部緊鄰南華北盆地太康隆起,5組斷裂構(gòu)成開封凹陷邊界。開封凹陷印支期抬升剝蝕,缺失白堊系地層,侏羅系僅局部發(fā)育,厚度薄,但保留了較厚的三疊系地層,喜山期再次沉降,新生界地層穩(wěn)定沉積,厚度較大。開封凹陷地?zé)崴a給來源為西部山區(qū)地質(zhì)歷史時期的大氣降水[20-21]。開封凹陷平均大地?zé)崃髦禐?6~60 mW/m2,地溫梯度2.8~3.6 ℃/hm,整體按凹陷延伸方向自西向東增加,開封—蘭考一帶地溫梯度3.4~3.6 ℃/hm。研究熱儲層選擇為館陶組砂巖,以長石巖屑砂巖、巖屑長石砂巖為主,多為孔隙型膠結(jié),表現(xiàn)出近源沉積的特征。館陶組頂面埋深為1 250~1 350 m,底面埋深1 900~2 000 m,地層厚600 m左右,砂巖密度1.70~1.96 g/cm3,孔隙度15%~35%,平均22%,滲透率1×10-13m2。蘭考館陶組地?zé)崴疄镹aCL型水,礦化度14 000~24 000 mg/L,Ca2+含量390~880 mg/L,Mg2+含量80~120 mg/L,pH值在7以上,弱堿性水。蘭考區(qū)域單井涌水量為100~155 m3/h,出水溫度為68~76 ℃[22-23]。
假定熱儲層、基巖與蓋層為勻質(zhì)的多孔介質(zhì),回灌水的滲流過程滿足Darcy定律。根據(jù)流體的質(zhì)量守恒定律,式(1)為基巖在穩(wěn)態(tài)下的滲流控制方程,式(2)為基巖在瞬態(tài)下的滲流控制方程。對斷層內(nèi)部滲流場的控制,采用等效連續(xù)模型進行描述,式(3)為斷層在穩(wěn)態(tài)下的滲流控制方程,式(4)為斷層在瞬態(tài)下的滲流控制方程。
(1)
(2)
(3)
(4)
(5)
(6)
式中,p為熱儲層基巖內(nèi)部孔隙壓力;pf為斷層內(nèi)部的孔隙壓力;u為基巖內(nèi)流體的滲流速度;uf為斷層內(nèi)流體的滲流速度;Km為基巖滲透率;Kf為斷層的等效滲透率;μ為流體的動力黏度;Qm為單位厚度基巖內(nèi)流體的質(zhì)量流量;Qf為單位高度斷層內(nèi)流體的質(zhì)量流量;df為斷層的等效直徑;λ為流體密度。
對于多孔介質(zhì),熱傳遞方式主要為巖石顆粒間的熱傳導(dǎo)、孔隙流體間的熱傳導(dǎo)和熱對流。若基巖為勻質(zhì)多孔介質(zhì),根據(jù)能量守恒原理,穩(wěn)態(tài)與瞬態(tài)下基巖中的熱傳遞控制方程分別為式(7)與(8),穩(wěn)態(tài)與瞬態(tài)下斷層中的熱傳遞控制方程分別為式(9)與(10)。
ρCpuT+q=Q
(7)
(8)
dsQs-dsρCputT-tdsqs=Qfac
(9)
(10)
q=-keffT
(11)
qs=-kefftT
(12)
keff=θpkp+(1-θp)k
(13)
(ρCp)eff=θpρpCp,p+(1-θp)ρCp
(14)
式中,Cp為流體的恒壓熱容;Cp,p為基巖的恒壓熱容;q為總傳導(dǎo)熱通量;qs為斷層內(nèi)壁傳導(dǎo)熱通量;Q為熱源總熱量;ds為斷層壁厚;Qs為斷層壁熱儲量;Qfac為斷層內(nèi)熱儲量;keff為有效導(dǎo)熱系數(shù);kp為基巖的導(dǎo)熱系數(shù);k為斷層內(nèi)流體的導(dǎo)熱系數(shù);(ρCp)eff為有效體積熱容;ρp為基巖的密度;θp為體積分數(shù)。
2.2 算例設(shè)置
本文根據(jù)蘭考熱儲地層分布幾何特征,建立簡化的多孔隙率模型如圖1所示。模型幾何尺寸為5 000 m×5 000 m×3 000 m。蓋層、熱儲層與基巖深度分別為0~1 300 m、1 300~1 900 m和1 900~3 000 m。在本文模型中對采熱井與回灌井進行簡化處理,僅考慮裸井部分進行生產(chǎn)活動,忽略了其他井段在地層中的滲流與換熱作用,在幾何上使用與z方向垂直的直線段表示裸井。采熱井與回灌井位長度均為200 m,在模型中心呈對稱分布。取地溫梯度為0.034 ℃/m,地表溫度為20 ℃。熱突破的定義為采熱井的平均溫度降低1 ℃。
圖1 幾何模型示意Fig.1 Geometric model schematic
為了研究采灌井間距與滲透率對溫度與滲流場的影響,設(shè)置4個井間距:300、400、500、600 m,并在每個井間距內(nèi)設(shè)置4個滲透率Kn(n=1,2,3,4),取值詳情見表1,共計16種工況。
表1 滲透率變量設(shè)置Tab.1 Permeability variable setting
設(shè)置1年為1個運行周期,每個運行周期內(nèi)采熱井與回灌井實際運行時間設(shè)為0.4年,回灌率設(shè)置為100%,即回灌流速等于采熱井的抽水流速。設(shè)置時間周期函數(shù)A,循環(huán)周期為1年,具體見式(15)。式(16)為回灌井與采熱井的水流控制方程,式(17)為回灌井的熱源方程:
(15)
Mout=Min=βρA(t)
(16)
(17)
式中,t為運行時間;Min為回灌井質(zhì)量流率;Mout為采熱井質(zhì)量流率;β為采熱井的抽水流速;Qw為回灌井的流失熱量;Tinj為回灌井的回灌水溫,本文取值為20 ℃。
設(shè)置穩(wěn)定的壓力邊界條件15 MPa。首先計算采熱井水位在第10年下降100 m,計算在每種工況下第10年采熱井孔隙壓力下降2 MPa的抽水流速,并將其帶入對應(yīng)的后續(xù)計算中,分析對比不同井間距與滲透率下溫度場與滲流場的在100年運營期內(nèi)的時空演化規(guī)律。斷層厚度為0.01 m,孔隙度為0.2,滲透率為1×10-19m2。設(shè)置垂直與水平分析平面對計算結(jié)果進行分析(圖2)。本文模型材料參數(shù)見表2。
圖2 分析平面Fig.2 Analysis plane
表2 模型材料參數(shù)Tab.2 Model material parameters
經(jīng)過對模型的抽水流速進行試算,得到了每種工況下使10年內(nèi)水位降低20 m的采熱井抽水流速(圖3)。由圖3可知,在井間距為300 m時有最大抽水流速,隨著井間距減小,在井間距500 m抽水流速降為最低,而井間距為600 m較500 m下抽水流速小幅增加。
不同井間距下采熱井井底溫度歷史演化曲線如圖4所示。由圖4可知,由于每年的時間周期內(nèi),采熱運營僅為0.4年,非運營期間的井底溫度得到了熱量補償,因此井底溫度隨時間的變化為鋸齒狀,其中齒狀起伏越劇烈,熱補償效果越好。相較于其他井間距,當(dāng)井間距為300 m時采熱井底溫度隨時間下降最快。當(dāng)滲透率為K1和K2時,井間距為500 m與400 m下采熱井底溫度隨時間的變化不大。隨著滲透率增加,井間距對采熱井底溫度隨時間的變化影響逐漸增加。不同井間距下回灌井井底溫度歷史演化曲線如圖5所示。由圖5可知,相較于采熱井,回灌井的井底溫度在10年內(nèi)急劇下降,與回灌水溫(20 ℃)快速靠近,并且由于低溫區(qū)域隨著時間推移范圍擴大,熱補償對回灌井的井底溫度提升幅度也迅速降低。隨著時間推移,回灌井井底溫度逐漸與回灌水溫重合。
圖3 抽水流速隨井間距的變化Fig.3 Change of pumping flow rate with well spacing
圖4 不同井間距下采熱井井底溫度歷史演化曲線Fig.4 Historical evolution curves of bottom hole temperature of heat recovery wells under different well spacings
圖5 不同井間距下回灌井井底溫度歷史演化曲線Fig.5 Historical evolution curves of bottom hole temperature of recharge wells under different well spacings
第50年與第100年的采熱井底溫度隨井間距變化如圖6所示。分析圖6發(fā)現(xiàn),采熱井底溫度與井間距具有很強的相關(guān)性,隨著井間距的增加,回灌造成的地層溫度降低對采熱井底溫度的影響越小,但隨著時間增加而增大。500、600 m井間距下的采熱井底溫度隨時間變化較小,井間距300、400 m下井底溫度發(fā)生了較大變化。以熱儲層滲透率K4為例,分析運行第100年采熱周期流線隨井間距的變化(圖7)。由圖7可知,井間距對工作區(qū)域內(nèi)的滲流場的變化影響較小,所有的井間距下滲流場的特征較為一致。當(dāng)水通過回灌井進入地下后,由于采熱井與回灌井間存在壓力梯度造成的吸力,部分回灌水向采熱井流動,被吸向采熱井的回灌水運動軌跡呈現(xiàn)雙橢圓形。
圖6 第50與100年采熱井井底溫度隨井間距的變化Fig.6 Variation of bottom hole temperature of heat recovery wells with well spacing in the 50th and 100th years
圖7 運行第100年采熱周期滲流場隨井間距的變化Fig.7 The variation of seepage field with well spacing during the 100th year of operation of thermal recovery cycle
運行第100年采熱周期垂直平面內(nèi)不同井間距滲流場與溫度場如圖8所示。由圖8可知,在300 m井間距下,采熱井被整體包裹在低溫區(qū)域內(nèi),右側(cè)沒有形成明顯的冷鋒凸出面,說明300 m井間距下有熱突破現(xiàn)象發(fā)生,回灌水的利用效率最高。在400 m井間距下,僅有小部分的冷鋒凸出面與采熱井接觸,而500、600 m井間距下,采熱井沒有被冷鋒凸出面覆蓋。通過將回灌井布設(shè)在遠離生產(chǎn)井的位置,可以最大限度地減少該現(xiàn)象的產(chǎn)生,而回灌井靠近生產(chǎn)井的位置可以最大限度地提高回灌井的效益,因此必須在這2個相互矛盾的布局之間找到適當(dāng)?shù)钠胶狻2蔁峋臒嵬黄茣r間見表3。
圖8 運行第100年采熱周期垂直平面內(nèi)不同井間距 滲流場與溫度場Fig.8 The seepage field and temperature field of different well spacings in the vertical plane of the 100th year of the heat recovery cycle
表3 采熱井熱突破時間Tab.3 Thermal breakthrough time of production well
由表3可知,300、400 m井間距下均會產(chǎn)生熱突破;在500 m井間距下,僅在K4滲透率時產(chǎn)生熱突破;在600 m井間距下未產(chǎn)生熱突破。這說明在蘭考縣進行地?zé)豳Y源開發(fā)利用時,井間距不易小于500 m,這樣既可以保證在運行壽命期間不產(chǎn)生熱突破,又可以保證回灌水利用最大化。由圖9可知,熱儲層滲透率越高,抽水流速越高,即生產(chǎn)效率越高,這意味著在滿足回灌率為100%的情況下熱儲層滲透率越高冷鋒凸出面越大,越容易發(fā)生熱突破現(xiàn)象。
圖9 運行第100年采熱周期垂直平面內(nèi)不同滲透率下 滲流場與溫度場(500 m)Fig.9 The seepage field and temperature field of different permeability in the vertical plane of the 100th year of the heat recovery cycle(500 m)
為研究中深層地?zé)峄毓鄬Σ蔁徇^程的影響,依托河南蘭考縣館陶組砂巖熱儲資料,采用COMSOL有限元多場耦合數(shù)值模擬軟件,在考慮基巖和蓋層熱補償與滲流對地?zé)醿觽鳠徇^程的影響的情況下,在三維多孔隙率模型的基礎(chǔ)上利用Darcy定理和多孔介質(zhì)傳熱公式建立了含裂縫地?zé)醿訚B流—傳熱耦合模型,分析了蘭考熱儲層采熱運行期間的溫度場與流場變化,得出以下結(jié)論。
(1)相較于采熱井,回灌井的井底溫度在10年內(nèi)急劇下降,與回灌水溫(50 ℃)快速靠近,隨著時間推移,回灌井井底溫度逐漸與回灌水溫重合。
(2)井間距對工作區(qū)域內(nèi)滲流場變化的影響較小。當(dāng)水通過回灌井進入地下后,由于采熱井與回灌井間存在壓力梯度,回灌水運動軌跡呈現(xiàn)雙橢圓形。
(3)隨著井間距的增加,回灌造成的地層溫度降低對采熱井底溫度的影響越小,但隨著時間增加而增大。
(4)在蘭考縣進行地?zé)豳Y源開發(fā)利用時,井間距不易小于500 m,這樣既可以保證在運行壽命期間不產(chǎn)生熱突破,又可以保證回灌水利用最大化。