干熱巖電學(xué)數(shù)值模擬方法與微觀特性分析

謝關(guān)寶

1 頁巖油氣富集機(jī)理與有效開發(fā)國家重點(diǎn)實(shí)驗(yàn)室， 北京 102206 2 中石化石油工程技術(shù)研究院有限公司， 北京 102206

0 引言

干熱巖能源是一種高效低碳清潔能源，目前全球發(fā)現(xiàn)的干熱巖儲量是所有油氣總儲量的數(shù)十倍(冉恒謙等，2010；楊方等，2012；陸川和王貴玲，2015；許天福等，2012，2016；付亞榮等，2018)，其內(nèi)部不存在流體或僅有少量流體(或致密不透水)的高溫巖體，儲層溫度一般在180 ℃以上，裂縫或斷層較為發(fā)育.作為一種重要的非常規(guī)資源(或能源)，研究其高溫條件下巖石物理性性質(zhì)，可指導(dǎo)干熱巖熱儲資源的高效開發(fā)；但在一些特殊情況下，如裂縫巖石、高溫超高溫條件，常規(guī)巖石物理實(shí)驗(yàn)難以開展，巖石物理響應(yīng)及特征研究面臨極大挑戰(zhàn).

近年來，數(shù)字巖心技術(shù)得到迅猛發(fā)展，通過構(gòu)建數(shù)字巖心模型，分析巖石物理性質(zhì)，已成為常規(guī)巖石物理實(shí)驗(yàn)的有益補(bǔ)充.考慮巖石粒度分布，構(gòu)建多孔隙數(shù)字巖心模型，可以精細(xì)表征巖石內(nèi)部微觀結(jié)構(gòu)，模擬實(shí)際巖石的導(dǎo)電特性(Bakke and ?ren，1997； ?ren and Bakke，2003).基于X-CT掃描構(gòu)建碳酸鹽巖模型，可分析碳酸鹽巖力學(xué)特性(Arns et al.，2002，2004，2005a,b).對于頁巖及其他致密巖性，可以采用Nano-CT(Sisk et al.，2010；Song et al.，2018)、FIB-SEM(Sok et al.，2009；Curtis，2010；聶昕等，2016；Salvati and Korsunsky，2017)、MAPS(李潮流等，2016，王民等，2018)等多種技術(shù)手段，構(gòu)建納米數(shù)字巖心，分析頁巖孔隙結(jié)構(gòu)特征.干熱巖熱儲層位的巖石物理性質(zhì)，在高溫高壓環(huán)境下會發(fā)生變化，尤其是高溫引起的電學(xué)性質(zhì)變化，已成為識別干熱巖資源的重要特征之一(Spichak et al.，2007；Chen，2015).許多學(xué)者研究了電阻率與溫度之間的關(guān)系(Winsauer et al.，1952；Ucok et al.，1979；Kariya and Shankland，1983；孟銀生等，2010；余琪祥等，2011；童小龍等，2020).蘇聯(lián)科學(xué)家(Cheremensky，1982)，總結(jié)了巖石電阻率與溫度之間的經(jīng)驗(yàn)公式，但忽略了巖石熱開裂引起的電阻率異常變化.基于數(shù)字巖心分析技術(shù)，改變巖石微觀參數(shù)，可分析巖石的微觀物理性質(zhì)(Adler et al.，1990；Blunt et al.，2002；Yue et al.，2011；Faisal et al.，2017；Tan et al.，2021).基于楓丹白露砂巖孔隙網(wǎng)絡(luò)模型，不考慮水濕巖石水膜傳導(dǎo)作用，在低含水飽和度情況下，電阻率指數(shù)的模擬結(jié)果遠(yuǎn)大于實(shí)驗(yàn)結(jié)果(Knackstedt et al.，2007)；如考慮水膜傳導(dǎo)作用，巖石電阻率模擬結(jié)果與實(shí)驗(yàn)結(jié)果準(zhǔn)確吻合(Liu et al.，2009).基于格子玻爾茲曼法，可模擬氣水分離過程，分析天然氣在地層水中的潤濕性和溶解度，分析其對儲層巖石電性的影響(Jiang et al.，2011).基于有限元法，采用消息傳遞機(jī)制(MPI)并行算法，可模擬超過2000×2000×2000個(gè)像素的巖石電性參數(shù)(Amabeoku et al.，2013).目前，基于數(shù)字巖心的巖石物理特性數(shù)值模擬研究面臨一個(gè)重大挑戰(zhàn)，即圖像分辨率與巖心尺寸(即圖像分辨率與儲層非均質(zhì)性)之間的矛盾(孔強(qiáng)夫等，2015).

我國青海省共和盆地干熱巖發(fā)育，地?zé)豳Y源豐富.本研究以該地區(qū)熱儲層巖心為研究對象，采用X-CT、QEMSCAN、Maps等多種實(shí)驗(yàn)方法，構(gòu)建了多尺度、多組分、高分辨率數(shù)字巖心；采用分形布朗運(yùn)動(dòng)方法(fBm)，構(gòu)建了裂縫性數(shù)字巖心；考慮干熱巖儲層高溫的特點(diǎn)，在數(shù)值模擬中引入了溫度對巖石骨架的影響，并利用有限元法計(jì)算了巖石電阻率，分析了干熱巖電學(xué)特性微觀響應(yīng)機(jī)理，為干熱巖熱儲層位評價(jià)和地?zé)豳Y源開發(fā)提供了微觀理論基礎(chǔ).

1 多尺度數(shù)字巖心構(gòu)建

研究所用巖心取自共和盆地GR1井，其熱儲層位主要巖性為黑云母片麻巖、花崗巖和花崗閃長巖，孔隙度、滲透率都極低，裂縫和斷層較為發(fā)育(陳梓慧等，2015).研究鉆取了三顆巖心，如圖1所示，分別編號為No.1、No.2、No.3，對應(yīng)取心深度為2980 m、3380 m和3600 m，巖性均為花崗巖.采用X射線衍射(XRD)及覆壓孔滲獲實(shí)驗(yàn)取了巖心主要礦物成分與孔隙度等參數(shù)，如表1所示.實(shí)驗(yàn)結(jié)果表明，研究區(qū)干熱巖主要造巖礦物為石英、斜長石和鉀長石，含少量黏土礦物；孔隙度非常小(最小為0.44%)，巖心較為致密，其中φ3>φ2>φ1.為了構(gòu)建高分辨率數(shù)字巖心，采用QEMSCAN進(jìn)行定量礦物分析，采用Maps進(jìn)行微觀孔隙分析，采用X-CT掃描構(gòu)建三維數(shù)字模型.

首先，對樣品進(jìn)行二維大面積掃描電鏡成像(Maps，分辨率為250 nm，圖2a)和QEMSCAN礦物成分分析(圖2b)，確定巖心孔隙尺寸分布特征及礦物分布情況.根據(jù)巖心孔徑分布特征，選擇合適的分辨率進(jìn)行X-CT掃描. 然后，根據(jù)礦物分布特征來確定納米級別的掃描區(qū)域(圖2b中的黃圈)，鉆取子樣品進(jìn)行納米CT掃描，1號樣品子樣品記為GRY1，2號樣品子樣品記為GRY2，3號樣品子樣品記為GRY3-1.最后，由Maps圖像提取干熱巖樣品孔隙半徑特征，發(fā)現(xiàn)其孔隙尺寸很小，主要分布范圍為0.25～4.0 μm，后續(xù)采用的X-CT掃描分辨率為500 nm，其三維重構(gòu)圖如圖2c所示，但此分辨率下孔隙識別度很低，孔隙度僅0.12%.因此，將X-CT掃描分辨率提高至65 nm，子樣品名稱為GRY3-2，其三維重構(gòu)圖如圖2d所示，此時(shí)孔隙較為容易識別，孔隙度為1.39%.

圖1 干熱巖柱塞樣(a) 1號樣(25 mm); (b) 2號樣(25 mm); (c) 3號樣(25 mm).Fig.1 HDR plug cores(a) No.1 25 mm plug sample; (b) No.2 25 mm plug sample; (c) No.3 25 mm plug sample.

表1 干熱巖XRD和孔隙度實(shí)驗(yàn)分析Table 1 XRD and porosity measurement experiment of HDR

Maps掃描區(qū)域與QEMSCAN掃描區(qū)域一致，將兩者圖像重疊對比，可以確定Maps圖像中每個(gè)像素點(diǎn)對應(yīng)的灰度、礦物類型和微孔隙發(fā)育程度.根據(jù)Maps中礦物灰度值大致確定礦物密度排序，顏色越黑代表礦物密度越低，顏色越亮白代表礦物密度越高.通過3個(gè)Maps圖像分析可得到礦物的密度，密度大小排序?yàn)椋害咽ⅰ堞研遍L石<ρ鉀長石<ρ綠泥石+方解石<ρ黃鐵礦.這個(gè)密度排序?yàn)楹罄m(xù)X-CT重建三維灰度圖像的多閾值分割提供了參考依據(jù).

數(shù)字巖心組分劃分通常采用單閾值進(jìn)行分割，0表示孔隙格架，1表示骨架，但單一閾值分割方式會喪失某些礦物成分信息.因此，在X-CT掃描圖像的基礎(chǔ)上，結(jié)合QEMSCAN礦物識別技術(shù)和XRD礦物組分分析，依據(jù)前述Maps獲取的礦物密度排序，對干熱巖數(shù)字巖心進(jìn)行多閾值分割.數(shù)字巖心礦物組分分割結(jié)果顯示如圖2e、圖2f所示.可以看出，干熱巖樣品礦物以石英、斜長石和鉀長石為主，含少量綠泥石、蒙脫石等黏土礦物，含極少量鐵礦物，且?guī)r心孔隙度均很小.從表2可以看到，從數(shù)字巖心得到的各礦物組分與XRD礦物分析結(jié)果接近，因此，可將數(shù)字巖心各礦物組分的電學(xué)綜合響應(yīng)定義為骨架電阻率，根據(jù)理論方法和巖石物理實(shí)驗(yàn)，構(gòu)建骨架電阻率與溫度關(guān)系.通過多次掃描實(shí)驗(yàn)構(gòu)建的數(shù)字核心分別命名為GRS1、GRS2、GRS3-1、GRS3-2(65 nm).

2 電學(xué)特性微觀數(shù)值模擬方法

電學(xué)特性的數(shù)值模擬通常采用有限元方法(FEM)進(jìn)行.考慮溫度對巖心電阻率的影響，將溫度對巖心電學(xué)特性的影響，轉(zhuǎn)換為其對巖心骨架的影響；將不同溫度下的骨架電導(dǎo)率作為有限元輸入，模擬不同溫度下整體數(shù)字巖心電流，得到整個(gè)數(shù)字巖心電阻率.

表2 干熱巖礦物成分對比分析Table 2 Companion of the mineral composition of HDR from XRD and constructed digital cores

圖2 多種掃描實(shí)驗(yàn)及多組分高分辨三維數(shù)字巖心(a) Maps掃描結(jié)果； (b) QEMSCAN掃描結(jié)果； (c) 500 nmX-CT掃描結(jié)果； (d) 65 nmX-CT掃描結(jié)果； (e) 三維多組分模型； (f) 三維孔隙模型.Fig.2 Multiple scanning experiments and multi-component high-resolution 3D digital cores(a) High-resolution scanning by Maps; (b) mineral analysis by QEMSCAN; (c) Reconstruction of 3D grayscale model with 500 nm X-CT; (d) reconstruction of 3D grayscale model with 65 nm X-CT; (e) 3D multi component digital core, including potash (light purple), mixture of quartz and plagioclase (white), chlorite (red), and ferriferous minerals (yellow) ; (f) 3D pore space.

2.1 電阻率與溫度關(guān)系

干熱巖大多屬于酸性侵入巖，骨架導(dǎo)電能力很差，可看作絕緣體.當(dāng)溫度T升高時(shí)，巖心電阻率ρcore會降低，巖石電阻率和溫度有如下關(guān)系(Cheremensky，1982)：

(1)

式中，α是經(jīng)驗(yàn)系數(shù)；ρ0初始電阻率(Ωm)；T0為初始溫度(℃)；ρcore是巖心電阻率(Ωm)；T為溫度(℃).

選取1號、2號、3號巖心子樣本測量其不同溫度下電阻率，樣品分別命名為No.1-1、No.1-2、No.1-3、No.2、No.3-1、No.3-2；測量樣品單軸加壓(12 MPa)，利用電阻絲對樣品加熱，在測量溫度點(diǎn)穩(wěn)定15 min；測量溫度從室溫開始，至500 ℃，每50 ℃測量一次，每個(gè)測量點(diǎn)測量5次，取平均值作為該測量點(diǎn)數(shù)值，結(jié)果如圖3所示.電阻率的對數(shù)與溫度大致成線性相關(guān)特征，與式(1)的描述一致.擬合可以得式(1)中系數(shù)α(本地區(qū)為0.02)，ρ0在150 ℃時(shí)為7.0×108Ωm.

圖3 不同溫度干熱巖電阻率Fig.3 The measured HDR resistivity results at different temperatures

干熱巖體可視作一種多孔介質(zhì)復(fù)合材料，包括骨架和孔隙兩部分.假定孔隙不導(dǎo)電，則巖心電阻率ρrock與骨架電阻率ρskeleton有如下關(guān)系式：

(2)

式中，φ為孔隙度(小數(shù))；ρskeleton為骨架電阻率(Ωm)；ρrock為巖石電阻率(Ωm).

通過巖心測量的巖心電阻率，可求解骨架電阻率，將其代入有限元算法中，即可模擬不同溫度下干熱巖巖心電學(xué)參數(shù).

2.2 有限元算法

有限元基本思想是變分原理，即將求解每個(gè)像素上的電壓分布問題，轉(zhuǎn)化為求解系統(tǒng)線性能量極值問題，從而確定數(shù)字巖心有效電導(dǎo)率.為使能量En取極小值，需滿足En對變量um(結(jié)點(diǎn)電壓)的偏導(dǎo)數(shù)均為零，即：

(3)

在三維數(shù)字巖心中，由于其自身就是由離散的像素構(gòu)成，選定每一個(gè)像素為一個(gè)單元，單元的形狀為立方體，每個(gè)像素通過8個(gè)頂點(diǎn)與相鄰單元連接，像素的能量由各個(gè)節(jié)點(diǎn)上的電壓確定，則像素內(nèi)的能量為

(4)

式中，σpq電導(dǎo)率張量；e(ex，ey，ez)是像素內(nèi)點(diǎn)(x，y，z)處的局部場.

在有限元方法中每個(gè)像素的能量可以表示為

(5)

式中，矩陣Drs是三維模型中的8×8矩陣，也稱為“剛度矩陣”，ur是像素上第r(1，2,…，8)個(gè)節(jié)點(diǎn)處的電壓；us是像素上第s(1，2,…，8) 個(gè)節(jié)點(diǎn)處的電壓.

將所有像素能量相加，可以計(jì)算整個(gè)三維數(shù)字巖心的總能量，將系統(tǒng)總能量取最小值，計(jì)算所有節(jié)點(diǎn)電壓.如果已知每個(gè)像素單元的電導(dǎo)率張量，可以求得每個(gè)單元內(nèi)的電流.然后，對三維數(shù)字巖心中所有單元電流取體積平均，得到三維數(shù)字巖心的等效電流.最后，根據(jù)電流與電壓的關(guān)系，得到巖樣等效電阻率.

3 電學(xué)微觀響應(yīng)特征分析

設(shè)置模擬溫度10～500 ℃，每50 ℃模擬一次，數(shù)字巖心代表體積元設(shè)置為400體素，模擬結(jié)果以Z-Axis為準(zhǔn)，計(jì)算不同溫度下的數(shù)字巖心電阻率，模擬結(jié)果如圖4所示.隨溫度升高，干熱巖電阻率呈減小變化趨勢，與式(1)的描述是一致的，且GRS1、GRS2數(shù)值模擬結(jié)果與實(shí)驗(yàn)結(jié)果也相吻合；但GRS3-1、GRS3-2的數(shù)值模擬結(jié)果與實(shí)驗(yàn)結(jié)果不一致，在其物理實(shí)驗(yàn)中間段，變化范圍較大且呈現(xiàn)非線性特征，與式(1)相差較大.

是什么原因?qū)е?號巖心電阻率變化趨勢與1號、2號不同呢？與1號、2號樣品相比，3號樣品埋藏深度更深，儲層溫度更高，熱應(yīng)力超過巖心抗壓強(qiáng)度臨界值，導(dǎo)致熱開裂或巖心破裂，如圖5所示.此外，在3號樣品的Maps和X-CT圖像中可以清楚地看到一些微裂縫.實(shí)驗(yàn)室加熱和原生地下高溫條件的共同作用導(dǎo)致3號樣品的裂縫更加發(fā)育，甚至破裂.

圖4 不同溫度下干熱巖電阻率數(shù)值模擬結(jié)果(a) GRS1； (b) GRS2； (c) GRS3-1； (d) GRS3-2.Fig.4 Resistivity numerical simulation results of HDR at different temperatures(a) GRS1； (b) GRS2； (c) GRS3-1； (d) GRS3-2.

圖5 3號巖樣裂縫顯示(a) Maps圖像； (b) 65 nm X-CT圖像； (c) 巖心開裂照片.Fig.5 Fractures in maps image; X-CT image and sample of No. 3 HDR(a) Maps image; (b) 65 nm X-CT image; (c) Rock sample after high temperature heating.

圖6為3號樣品的2個(gè)子樣品不同溫度下的電阻率數(shù)值模擬結(jié)果.由圖可見，隨著溫度升高，其電阻率在300 ℃后發(fā)生較大差異變化.圖6中虛線指示的是兩個(gè)子樣品的平均電阻率，當(dāng)溫度達(dá)到300 ℃時(shí)，電阻率變化趨勢發(fā)生了變化，故推測此時(shí)巖樣產(chǎn)生了裂縫.因此，加熱過程中產(chǎn)生裂縫的溫度邊界約為300 ℃.

圖6 3號樣品不同溫度下電阻率： 熱開裂溫度邊界約為300 ℃Fig.6 Resistivity measurement results of two subsamples from the No.3 sample at different temperatures, showing the temperature boundary of thermal cracking is about 300 ℃.

為進(jìn)一步分析熱裂解產(chǎn)生裂縫對電阻率的影響，研究采用分形布朗運(yùn)動(dòng)(fBm)方法在3號數(shù)字巖心中人為添加裂縫(Voss，1988； Kim， 2007)，形成帶裂縫數(shù)字巖心模型，命名為GRS3-3.因研究區(qū)裂縫多為低傾角，因此模擬中添加的裂縫均為水平裂縫.為了模擬加熱和熱開裂實(shí)際情況，模擬中隨著溫度的升高，生成的裂縫數(shù)量也相應(yīng)增加.具體模擬過程及參數(shù)如下：當(dāng)300 ℃時(shí)，添加第一條裂縫，其長度200體素，開度為2體素；當(dāng)400 ℃時(shí)，添加第二條裂縫，450 ℃時(shí)添加第三條裂縫；隨著溫度升高，先前生成裂縫的尺度會增大(假設(shè)長度每升高50 ℃增長100體素)；當(dāng)溫度升至450 ℃(或第一條裂縫貫穿巖心)時(shí)，則不再增加裂縫條數(shù).整個(gè)添加裂縫的過程如圖7a—圖7c所示，最后構(gòu)建的裂縫性數(shù)字巖心模型如圖7d所示.

對裂縫性數(shù)字巖心進(jìn)行電學(xué)特性微觀數(shù)值模擬，結(jié)果如圖8所示.由圖可見，在300 ℃前，電阻率仍呈線性下降趨勢；加入裂縫后，由于裂縫內(nèi)空氣不導(dǎo)電，電阻率減少趨勢變緩，并保持穩(wěn)定(在巖心礦物相變前，如礦物脫水).此時(shí)，含裂縫數(shù)字巖心GRS3-3模擬的電阻率，與巖石實(shí)驗(yàn)測得的平均電阻率變化趨勢一致，這說明裂縫性數(shù)字巖心構(gòu)建與數(shù)值模擬方法是正確的.在獲取巖樣或進(jìn)行高溫巖石物理實(shí)驗(yàn)受限的情況下，利用上述數(shù)字巖心方法，能夠更方便地分析干熱巖高溫電學(xué)性質(zhì).

圖7 不同溫度增加裂縫模型示意圖(a) 350 ℃時(shí)裂縫模型； (b) 400 ℃時(shí)裂縫模型； (c) 450 ℃時(shí)裂縫模型 ； (d) 450 ℃時(shí)多裂縫模型.Fig.7 Schematic diagram of adding fractures at different temperatures(a) Adding the first fracture at 350 ℃；(b) Adding the second fracture at 400 ℃；(c) Adding the third fracture at 450 ℃； (d) The combined digital core with multipal fractures at 450 ℃.

圖8 3號巖樣電響應(yīng)數(shù)值模擬結(jié)果Fig.8 Resistivity simulation results of two sub-cores from No.3 sample after adding fractures

4 結(jié)論

為分析不同溫度下干熱巖電學(xué)特性與響應(yīng)機(jī)理，構(gòu)建了多尺度、多組分、高分辨率數(shù)字巖心三維模型，分析了多尺度孔隙結(jié)構(gòu)及和礦物成分特征，并研究了不同溫度下干熱巖電阻率變化規(guī)律.具體認(rèn)識如下：

(1)采用X-CT、Maps、QEMSCAN和XRD多種實(shí)驗(yàn)分析方法，構(gòu)建了三維多尺度、多組分、高分辨率數(shù)字巖心模型，分析了干熱巖儲層孔隙特征和礦物成分特征，為研究干熱巖電學(xué)特性提供了有效載體.

(2)干熱巖電阻率隨溫度升高呈線性下降趨勢.對于深埋層伴有熱開裂的巖石，300 ℃前電阻率仍呈線性下降，300 ℃后電阻率變化幅度減小并趨于穩(wěn)定.

(3)裂縫性數(shù)字巖心模擬與高溫巖石物理實(shí)驗(yàn)結(jié)果一致，表明巖心加熱過程中可加熱過程中可能出現(xiàn)熱開裂.

(4)數(shù)字巖心技術(shù)可快速研究干熱巖高溫石物理特性，可有效解決現(xiàn)有巖石物理實(shí)驗(yàn)無法定量分析的問題，為研究干熱巖熱儲層高溫巖石特性提供了一種新方法.

干熱巖高溫電學(xué)特性在300 ℃出現(xiàn)拐點(diǎn)現(xiàn)象，文中利用裂縫性數(shù)字巖心進(jìn)行了驗(yàn)證，物模結(jié)果與數(shù)模結(jié)果有較好的一致性；但溫度引起的干熱巖“裂紋到縫網(wǎng)再到破裂”這個(gè)過程是比較復(fù)雜的，還需綜合考慮干熱巖體的縫網(wǎng)拓展機(jī)理、巖心礦物相變、礦物會脫水、原位地層應(yīng)力條件等因素，才能夠?qū)Ω蔁釒r儲層巖石物理性質(zhì)有更深刻、更準(zhǔn)確的認(rèn)識.