基于數(shù)字巖心研究流體性質(zhì)對(duì)裂縫性低滲透儲(chǔ)層彈性參數(shù)的影響規(guī)律

孫建孟,閆國亮,姜黎明,崔利凱,趙建鵬,崔紅珠

(1.中國石油大學(xué)地球科學(xué)與技術(shù)學(xué)院,山東青島266580;2.中國石油勘探開發(fā)研究院西北分院,甘肅蘭州730020; 3.中國石油集團(tuán)測(cè)井有限公司技術(shù)中心,陜西西安710077;4.中國石油集團(tuán)測(cè)井有限公司長(zhǎng)慶事業(yè)部,陜西西安710077)

基于數(shù)字巖心研究流體性質(zhì)對(duì)裂縫性低滲透儲(chǔ)層彈性參數(shù)的影響規(guī)律

孫建孟1,閆國亮2,姜黎明3,崔利凱1,趙建鵬1,崔紅珠4

(1.中國石油大學(xué)地球科學(xué)與技術(shù)學(xué)院,山東青島266580;2.中國石油勘探開發(fā)研究院西北分院,甘肅蘭州730020; 3.中國石油集團(tuán)測(cè)井有限公司技術(shù)中心,陜西西安710077;4.中國石油集團(tuán)測(cè)井有限公司長(zhǎng)慶事業(yè)部,陜西西安710077)

為研究裂縫及流體性質(zhì)對(duì)低滲透儲(chǔ)層彈性參數(shù)的影響規(guī)律,采用X射線CT掃描技術(shù)構(gòu)建低滲透儲(chǔ)層巖石的3維數(shù)字巖心,應(yīng)用圖像處理算法加入定向排列的平行便士狀裂縫,形成橫向各向同性數(shù)字巖心,采用擴(kuò)展的有限元方法計(jì)算含有裂縫的數(shù)字巖心的彈性模量,分析裂縫及流體性質(zhì)對(duì)其彈性模量的影響規(guī)律。結(jié)果表明:裂縫縱橫比和密度保持不變,隨流體體積模量的增加,縱波各向異性參數(shù)呈線性趨勢(shì)減小;裂縫密度和流體性質(zhì)不變,隨裂縫縱橫比增加,儲(chǔ)層巖石的彈性模量均呈線性增加;裂縫縱橫比和流體性質(zhì)不變,隨裂縫密度增加,儲(chǔ)層巖石的彈性模量也呈線性減小。

測(cè)井;裂縫性低滲透儲(chǔ)層;巖石性質(zhì);彈性模量;數(shù)字巖心;有限元方法

由于定向應(yīng)力的作用,低滲透儲(chǔ)層存在的自然裂縫或壓裂縫通常表現(xiàn)為具有一定取向且相互平行的裂縫模式,形成裂縫性低滲透儲(chǔ)層[1]。裂縫的存在及裂縫中填充流體的性質(zhì)對(duì)裂縫巖石的彈性性質(zhì)有重要的影響,會(huì)表現(xiàn)出大的地震各向異性。定量研究裂縫巖石彈性性質(zhì)隨所含流體變化的特征對(duì)于深入了解巖石物理性質(zhì),特別是對(duì)油氣田開發(fā)具有重要意義。由于缺少表征巖石全面結(jié)構(gòu)的資料,巖石物理模型主要局限于經(jīng)驗(yàn)公式[2]、上下邊界理論[3]、有效介質(zhì)理論[4-9],這些都不完全令人滿意。隨著計(jì)算機(jī)技術(shù)的發(fā)展,可以基于微觀結(jié)構(gòu)的數(shù)字模型來直接求解線性彈性方程,從而得到巖石的彈性參數(shù)。Knackstedt[10]利用有限元方法計(jì)算了楓丹白露砂巖的體積模量與剪切模量,數(shù)值計(jì)算結(jié)果與Han[2]的實(shí)驗(yàn)結(jié)果一致。Roberts[11]利用有限元方法計(jì)算了隨機(jī)多孔材料的彈性性質(zhì),結(jié)果表明在整個(gè)孔隙度范圍內(nèi),多孔材料的楊氏模量與固體相的泊松比無關(guān)。Dina[12]等基于數(shù)字巖心用有限元的方法研究了部分飽和巖石的彈性性質(zhì),并把計(jì)算結(jié)果與低頻Gassmann-Wood方程和高頻Gassmann-Hill方程進(jìn)行了比較,發(fā)現(xiàn)巖石均勻飽和情況下計(jì)算結(jié)果與低頻Gassmann-Wood方程計(jì)算結(jié)果一致。但是,上述研究都沒有進(jìn)行裂縫巖石彈性性質(zhì)隨所含流體變化特征的研究。筆者基于重建的裂縫性低滲透儲(chǔ)層數(shù)字巖心,采用有限元方法計(jì)算不同裂縫參數(shù)和流體性質(zhì)下彈性模量的變化規(guī)律并與沒有裂縫發(fā)育的巖石進(jìn)行對(duì)比,分析裂縫開度、長(zhǎng)度及流體性質(zhì)對(duì)裂縫性低滲透儲(chǔ)層彈性參數(shù)的影響規(guī)律。

1 裂縫性低滲透儲(chǔ)層巖石數(shù)字巖心建模

裂縫性低滲透儲(chǔ)層巖石一般具有各向異性特征,各向異性的存在使巖石的彈性性質(zhì)研究變得尤為復(fù)雜。一般情況下,巖石的各向異性特征可以用相對(duì)簡(jiǎn)單的橫向各向同性來模擬。為了研究橫向各向同性巖石的彈性性質(zhì),構(gòu)建了具有橫向各向同性的裂縫性低滲透儲(chǔ)層巖石的數(shù)字巖心。建模步驟分為兩步:首先獲取具有各向同性特征的低滲透儲(chǔ)層巖石基質(zhì)的3維數(shù)字巖心;然后在巖石基質(zhì)中加入一系列定向排列的裂縫,得到具有橫向各向同性特征的裂縫性數(shù)字巖心。

目前,建立巖石基質(zhì)3維數(shù)字巖心的方法主要有兩大類:一類是數(shù)學(xué)方法,通過巖石2維薄片圖像獲得的統(tǒng)計(jì)信息來重建數(shù)字巖心,包括隨機(jī)法和過程模擬法兩種;另一類是物理方法,通過實(shí)驗(yàn)儀器掃描直接構(gòu)建數(shù)字巖心,包括切片組合掃描法和X-CT掃描法兩種。與數(shù)學(xué)方法相比,物理方法構(gòu)建的3維數(shù)字巖心具有直觀、準(zhǔn)確的優(yōu)點(diǎn)。

切片組合掃描法通過逐層刨光和切割巖心并進(jìn)行高精度成像得到,不僅會(huì)破壞巖心的完整性,而且非常耗時(shí)。X-CT掃描法可以無損地檢測(cè)非透明物體的組成及結(jié)構(gòu)[13-14],因此采用該方法構(gòu)建低滲透儲(chǔ)層巖石的數(shù)字巖心。X射線CT掃描系統(tǒng)得到的最原始的圖像是反映巖心徑向X射線吸收系數(shù)的一系列投影圖像,應(yīng)用濾波(卷積)反投影技術(shù)處理后得到巖心的3維灰度圖像。

圖1(a)為某低滲透儲(chǔ)層砂巖巖心(巖心編號(hào)為CJ157,實(shí)驗(yàn)測(cè)量孔隙度為9.2%)掃描得到的3維灰度圖像的一個(gè)切片,掃描分辨率為2μm/像素。圖像中像素點(diǎn)灰度值反映了巖石各組分密度,像素點(diǎn)越亮,即灰度值越大,表示密度越大。從圖中可以直接識(shí)別出孔隙的位置(圖像上比較暗的像素點(diǎn)),但孔隙與骨架的邊緣比較模糊,需要對(duì)圖像進(jìn)行濾波處理,從而提高孔隙和骨架的對(duì)比度[15]。圖1 (b)為采用中值濾波算法對(duì)圖1(a)濾波處理后的結(jié)果。最后選取合適的閾值將灰度圖像分割為二值圖像,分割結(jié)果見圖1(c)。將二值化后的圖像逐層疊合,得到CJ157巖心的3維數(shù)字巖心,如圖1(d)所示。圖中紅色部分為孔隙,藍(lán)色部分為骨架,孔隙度為8.8%,與CJ157巖心實(shí)驗(yàn)測(cè)量的孔隙度9.2%接近。

圖1 CJ157砂巖的數(shù)字巖心Fig.1 Digital core of CJ157

通過數(shù)字圖像處理技術(shù)把一系列定向排列的便士形狀裂縫嵌入巖石基質(zhì)中,從而得到橫向各向同性的數(shù)字巖心,如圖2所示。圖中便士形狀裂縫的3個(gè)半軸滿足a2=a3?a1。圖3為施加裂縫后CJ157砂巖數(shù)字巖心的橫截面。裂縫密度和縱橫比是定量描述裂縫的兩個(gè)重要參數(shù)。其中縱橫比a定義為:α=a1/a2;裂縫密度e計(jì)算公式為

式中,N為裂縫總的數(shù)目;V為數(shù)字巖心的體積;φc為裂縫孔隙度。在數(shù)值模擬中設(shè)定便士形狀裂縫厚度為3個(gè)像素點(diǎn),直徑為24個(gè)像素點(diǎn),則縱橫比為0.125。

圖2 施加裂縫后CJ157砂巖的數(shù)字巖心Fig.2 Digital core of sandstone CJ157 with cracks

圖3 施加裂縫后CJ157砂巖數(shù)字巖心的水平橫截面Fig.3 Horizontal cross-section of digital core of sandstone CJ157 with cracks

2 理論和數(shù)值模擬方法

2.1 各向異性Gassmann理論

各向同性Gassmann流體替換理論適用于各向同性巖石而不適用于各向異性巖石流體替換的研究,因此不能用于裂縫中飽和流體巖石彈性性質(zhì)的研究。Gassmann[16]隨后給出了適用于各向異性介質(zhì)的流體替換理論,稱為各向異性Gassmann理論。該理論假定在巖石孔隙空間中流體壓力是均衡的,所以只適用于低頻條件下各向異性巖石彈性性質(zhì)的研究。

干巖石的彈性模量和流體飽和巖石彈性模量的關(guān)系滿足:

對(duì)于橫向各向同性巖石,其彈性剛度矩陣中含有5個(gè)獨(dú)立的彈性參數(shù)。式(2)中比例系數(shù)αi、αj可以表示為

Thomsen各向異性參數(shù)ε、γ、δ可以用來描述橫向各向同性介質(zhì)的各向異性程度[17],具體定義為

式中,ε為度量準(zhǔn)縱波各向異性的參數(shù),ε越大,介質(zhì)的縱波各向異性越大,ε=0時(shí),縱波無各向異性; γ可以看作是度量準(zhǔn)橫波各向異性強(qiáng)度,或橫波分裂強(qiáng)度的參數(shù);δ主要控制水平垂直正交方向上的各向異性強(qiáng)弱。

2.2 有限元方法

由于數(shù)字巖心是由離散的體素組成,每一個(gè)體素可以視為一個(gè)立方體單元,因此應(yīng)用有限元方法計(jì)算巖石物理屬性時(shí)無需進(jìn)行網(wǎng)格劃分。在3維數(shù)字巖心的基礎(chǔ)上,利用有限元法計(jì)算巖石彈性模量的理論基礎(chǔ)是變分原理。對(duì)于給定的數(shù)字巖心,沿主應(yīng)力和切應(yīng)力方向分別施加一個(gè)宏觀應(yīng)變,通過使系統(tǒng)的彈性自由能En最小來確定每個(gè)像素點(diǎn)上的最終彈性位移分布。根據(jù)變分原理,求解每個(gè)像素點(diǎn)上的位移分布問題轉(zhuǎn)化為求解系統(tǒng)線性彈性自由能極值的問題,并最終確定數(shù)字巖心的有效彈性模量。為使能量En取極小值,需滿足能量對(duì)變量um(結(jié)點(diǎn)彈性位移)的偏導(dǎo)數(shù)均為零,即

在數(shù)值求解過程中,當(dāng)能量En對(duì)第m個(gè)結(jié)點(diǎn)彈性位移的偏導(dǎo)數(shù)構(gòu)成的梯度矢量的平方和小于某一給定允許誤差時(shí),可近似認(rèn)為等式(7)成立,即確定了3維數(shù)字巖心中的應(yīng)力分布和有效彈性參數(shù)。

對(duì)于橫向各向同性巖石,具有5個(gè)獨(dú)立參數(shù),執(zhí)行一次數(shù)值模擬難以求取所有的彈性參數(shù),因此可以利用彈性方程σ=Cε的線性關(guān)系,引入6個(gè)正交應(yīng)變基矢,其中α表示基矢的數(shù)目,i(i=1,…,6)表示Voigt指數(shù),應(yīng)用這6個(gè)正交基矢足以求出剛度矩陣的36個(gè)元素。圖4為6個(gè)應(yīng)變基矢的示意圖,其中(a)～(c)為軸向應(yīng)變,(d)～(f)為切向應(yīng)變,具體表達(dá)式為ε1=(1,0,0,0,0,0),ε2=(0,1,0,0, 0,0),ε3=(0,0,1,0,0,0),ε4=(0,0,0,1,0,0),ε5= (0,0,0,0,1,0),ε6=(0,0,0,0,0,1)。通過選其中一個(gè)應(yīng)變基矢作為宏觀應(yīng)變條件執(zhí)行模擬,可以得到6個(gè)宏觀應(yīng)力。根據(jù)彈性方程的線性關(guān)系,這6個(gè)宏觀應(yīng)力定義了剛度矩陣Cαβ的一列,因此分別獨(dú)立執(zhí)行6次模擬便可以得到剛度矩陣的所有元素。由于橫向各向同性巖石含有5個(gè)獨(dú)立參數(shù),所以僅需執(zhí)行兩次模擬,就可以得到所求參數(shù)。擴(kuò)展的數(shù)值模擬算法不僅能計(jì)算各向同性巖石的彈性參數(shù),也能計(jì)算具有橫向各向同性巖石的彈性參數(shù)。

圖4 6個(gè)宏觀應(yīng)變基矢Fig.4 Six basis vectors of macro-strain

3 裂縫性低滲透儲(chǔ)層彈性性質(zhì)的影響因素

3.1 數(shù)值算法的驗(yàn)證

圖5 剛度矩陣元素與流體模量的關(guān)系Fig.5 Relationship between stiffness matrix and fluid modulus

基質(zhì)巖石的剛度矩陣與流體模量的關(guān)系如圖5 (a)所示。從圖中可以看出,所選的基質(zhì)巖石是各向同性的,c11與c33相等,并隨流體模量的增加而增加; c44與c66相等,即巖石的剪切模量μ,這兩個(gè)量與流體性質(zhì)無關(guān),這與各向同性Gassmann理論也是相符的。通過對(duì)比數(shù)值模擬結(jié)果(圖中散點(diǎn))與利用各向同性流體替換Gassmann理論計(jì)算的結(jié)果(圖中實(shí)線),發(fā)現(xiàn)數(shù)值模擬結(jié)果和理論計(jì)算結(jié)果吻合非常好,證實(shí)了數(shù)值算法的正確性。

3.2 流體性質(zhì)對(duì)儲(chǔ)層彈性性質(zhì)的影響

基于構(gòu)建的裂縫性低滲透儲(chǔ)層橫向各向同性數(shù)字巖心,利用有限元方法計(jì)算了巖石的彈性參數(shù)并與各向異性Gassmann理論計(jì)算結(jié)果進(jìn)行了比較,如圖5 (b)所示。圖中給出了含裂縫巖石的剛度矩陣與流體模量的關(guān)系,其中裂縫密度為0.05。從圖上可以看出嵌入基質(zhì)巖石的平行裂縫引起了顯著的各向異性。剛度矩陣元素c11和c33發(fā)生偏離,并且在平行于對(duì)稱軸方向(即c33方向)巖石變得更軟。c44與c66隨流體模量的增大產(chǎn)生了很小的偏差,這說明流體對(duì)切向彈性參數(shù)影響比較小。通過對(duì)比數(shù)值模擬結(jié)果(圖中散點(diǎn))與利用各向異性流體替換Gassmann理論計(jì)算的結(jié)果(圖中實(shí)線),發(fā)現(xiàn)數(shù)值模擬結(jié)果和理論計(jì)算結(jié)果基本吻合,說明所擴(kuò)展的數(shù)值模擬算法用于氣層巖石各向異性研究是可行性的。基于擴(kuò)展的數(shù)值模擬算法進(jìn)一步研究了各向異性參數(shù)隨流體模量的變化關(guān)系,如圖6所示。從圖中可以看出,當(dāng)流體模量為0時(shí),即干巖石,縱波各向異性參數(shù)ε最大,隨流體體積模量的增加而減小。

圖6 縱波各向異性參數(shù)ε與流體體積模量的關(guān)系Fig.6 Relationship between anisotropy parameter ε of longitudinal wave and fluid bulk modulus

3.3 裂縫對(duì)儲(chǔ)層彈性性質(zhì)的影響

圖7 剛度矩陣元素與裂縫縱橫比的關(guān)系Fig.7 Relationship between stiffness matrix and crack aspect ratios

便士形狀裂縫的縱向半軸a1的2倍對(duì)應(yīng)裂縫的開度h,橫向半軸a2的2倍對(duì)應(yīng)裂縫的長(zhǎng)度l。圖7為裂縫長(zhǎng)度為48 μm,裂縫中分別飽和天然氣、油和水,通過改變裂縫開度得到的不同裂縫縱橫比情況下CJ157數(shù)字巖心的剛度矩陣元素變化圖。計(jì)算過程中天然氣的體積模量為0.2 GPa,油的體積模量為1.02 GPa,水的體積模量為2.2 GPa。從圖中可以看出,隨著裂縫縱橫比的增大,剛度矩陣元素均呈線性增大。

圖8為裂縫開度為6 μm,裂縫長(zhǎng)度為48 μm,不同裂縫密度情況下CJ157數(shù)字巖心的剛度矩陣元素變化圖。從圖中可以看出,隨著裂縫密度的增大,剛度矩陣元素也呈線性減小。

圖8 剛度矩陣元素與裂縫密度的關(guān)系Fig.8 Relationship between stiffness matrix and crack density

4 結(jié) 論

(1)采用擴(kuò)展的有限元方法計(jì)算的裂縫性低滲透儲(chǔ)層彈性參數(shù)和理論計(jì)算結(jié)果基本吻合,驗(yàn)證了所擴(kuò)展的數(shù)值模擬算法用于橫向各向同性介質(zhì)彈性性質(zhì)研究的可行性。

(2)當(dāng)裂縫縱橫比和密度保持不變,流體模量為0時(shí),縱波各向異性參數(shù)ε最大;隨流體體積模量的增加,縱波各向異性參數(shù)ε呈線性減小。

(3)裂縫密度和流體性質(zhì)保持不變,隨裂縫縱橫比增加,裂縫性低滲透儲(chǔ)層巖石的彈性模量均呈線性增大。

(4)裂縫縱橫比和流體性質(zhì)保持不變,隨裂縫密度增加,裂縫性低滲透儲(chǔ)層巖石的彈性模量呈線性減小。

[1] 曾聯(lián)波.低滲透砂巖儲(chǔ)層裂縫的形成與分布[M].北京:科學(xué)出版社,2008.

[2] HAN D H.Effects of porosity and clay content on acoustic properties of sandstones and unconsolidated sediments [D].Standford:Standford University,1986.

[3] HASHIN Z,SHTRIKMAN S A.variational approach to the theory of the elastic behaviour of multiphase materials [J].J Mech Phys Solids,1963,11(2):127-140.

[4] BERRYMAN J G.Mixture theories for rock properties,in a Handbook of physical Constants[M].Washington:A-merican Geophysical Union,1995:205-228.

[5] BERRYMAN J G.Effective medium theories for multicomponent poroelastic composites[J].Journal of Engineering Mechanics,2006,132(5):519-531.

[6] CLEARY M P,CHEN I W,LEE S M.Self-consistent techniques for heterogeneous media[J].J Eng Mech, 1980,106(5):861-887.

[7] NORRIS A N.A differential scheme for the effective moduli of composites[J].Mech Mater,1985,4:1-16.

[8] ZIMMERMAN R W.Compressibility of sandstones[M]. New York:Elsevier,1991.

[9] BERRYMAN J G.Single-scattering approximations for coefficients in Biot?s equations of poroelasticity[J].J Acoust Soc Am,1992,91(2):551-571.

[10] KNACKSTEDT M A,ARNS C H,PINCZEWSKIZ W V, et al.Computation of linear elastic properties from microtomographic images:methodology and agreement between theory and experiment[J].Geophysics,2002,67 (5):1396-1405.

[11] ROBERTS A P,GARBOCZI E J.Computation of the linear elastic properties of random porous materials with a wide variety of microstructure[J].Proceedings of the Royal Society of London,2002,458:1033-1054.

[12] DINA M,BORIS G,RADIM C,et al.Finite element modelling of the effective elastic properties of partially saturated rocks[J].Computers&Geosciences,2008, 34(6):647-657.

[13] DUNSMOIR J H,FERGUSON S R,D?AMICO K L,et al.X-ray microtomography:a new tool for the characterization of porous media[R].SPE22860,1991.

[14] 康志勤,趙陽升,孟巧榮,等.油頁巖熱破裂規(guī)律顯微CT實(shí)驗(yàn)研究[J].地球物理學(xué)報(bào),2009,52(3): 842-848. KANG Zhi-Qin,ZHAO Yang-Sheng,MENG Qiao-Rong,et al.Micro-CT experimental research of oil shale thermal cracking laws[J].Chinese Journal of Geophysics,2009,52(3):842-848.

[15] DONG H,TOUATI M,MARTIN J B.Pore network modeling:analysis of pore size distribution of arabian core samples[R].SPE105156,2007.

[16] GASSMANN F.Uber die elastizitat porous media[J]. Vier der Natur Gesellschaft in Zurich,1951,96(1):1-23.

[17] THOMSEN L.Weak elastic anisotropy[J].Geophysics,1986,51:1954-1966.

(編輯 修榮榮)

Research of influence laws of fluid properties on elastic parameters of fractured low permeability reservoir rocks based on digital core

SUN Jian-meng1,YAN Guo-liang2,JIANG Li-ming3,CUI Li-kai1,ZHAO Jian-peng1,CUI Hong-zhu4
(1.School of Geosciences in China University of Petroleum,Qingdao 266580,China; 2.Research Institute of Petroleum Exploration&Development-Northeast,PetroChina,Lanzhou 730020,China; 3.Technology Center,China Petroleum Logging Company Limited,Xi?an 710077,China; 4.Changqing Division,China Petroleum Logging Company Limited,Xi?an 710077,China)

In order to investigate the effects of fracture and fluid properties on the elastic parameters of fractured low permeability reservoirs,the digital cores of low permeability reservoirs were constructed by X-ray CT.The transversely isotropic digital cores with parallel penny shape cracks were firstly constructed using digital image processing technology,then the elastic parameters of the cores were calculated by an extended finite element method,and lastly the effects of fracture and fluid properties on the elastic modulus were studied.The results show that with fixed crack aspect ratio and density,the longitudinal wave anisotropy linearly decreases with the increase of the fluid bulk modulus.With fixed crack density and fluid properties, the elastic moduli of fractured low permeability reservoir rock show a linear relationship with increasing crack aspect ratio. And with fixed crack aspect ratio and fluid properties,the elastic moduli of fractured low permeability reservoir rock show a linear relationship with increasing crack density.

well logging;fractured low permeability reservoir;rock properties;elastic modulus;digital core;finite element method

TE 19

:A

1673-5005(2014)03-0039-06

10.3969/j.issn.1673-5005.2014.03.006

2013-10-05

國家自然科學(xué)基金項(xiàng)目(41374124);國家重大科技專項(xiàng)(2011ZX05006-002)

孫建孟(1964-),男,教授,博士生導(dǎo)師,從事測(cè)井資料處理、解釋與巖石物理研究。E-mail:sunjm@upc.edu.cn。