基于數字巖心研究流體性質對裂縫性低滲透儲層彈性參數的影響規律

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

(1.中國石油大學地球科學與技術學院,山東青島266580;2.中國石油勘探開發研究院西北分院,甘肅蘭州730020; 3.中國石油集團測井有限公司技術中心,陜西西安710077;4.中國石油集團測井有限公司長慶事業部,陜西西安710077)

基于數字巖心研究流體性質對裂縫性低滲透儲層彈性參數的影響規律

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

(1.中國石油大學地球科學與技術學院,山東青島266580;2.中國石油勘探開發研究院西北分院,甘肅蘭州730020; 3.中國石油集團測井有限公司技術中心,陜西西安710077;4.中國石油集團測井有限公司長慶事業部,陜西西安710077)

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

測井;裂縫性低滲透儲層;巖石性質;彈性模量;數字巖心;有限元方法

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

1 裂縫性低滲透儲層巖石數字巖心建模

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

目前,建立巖石基質3維數字巖心的方法主要有兩大類:一類是數學方法,通過巖石2維薄片圖像獲得的統計信息來重建數字巖心,包括隨機法和過程模擬法兩種;另一類是物理方法,通過實驗儀器掃描直接構建數字巖心,包括切片組合掃描法和X-CT掃描法兩種。與數學方法相比,物理方法構建的3維數字巖心具有直觀、準確的優點。

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

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

圖1 CJ157砂巖的數字巖心Fig.1 Digital core of CJ157

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

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

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

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

2 理論和數值模擬方法

2.1 各向異性Gassmann理論

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

干巖石的彈性模量和流體飽和巖石彈性模量的關系滿足:

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

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

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

2.2 有限元方法

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

在數值求解過程中,當能量En對第m個結點彈性位移的偏導數構成的梯度矢量的平方和小于某一給定允許誤差時,可近似認為等式(7)成立,即確定了3維數字巖心中的應力分布和有效彈性參數。

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

圖4 6個宏觀應變基矢Fig.4 Six basis vectors of macro-strain

3 裂縫性低滲透儲層彈性性質的影響因素

3.1 數值算法的驗證

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

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

3.2 流體性質對儲層彈性性質的影響

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

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

3.3 裂縫對儲層彈性性質的影響

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

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

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

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

4 結 論

(1)采用擴展的有限元方法計算的裂縫性低滲透儲層彈性參數和理論計算結果基本吻合,驗證了所擴展的數值模擬算法用于橫向各向同性介質彈性性質研究的可行性。

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

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

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

[1] 曾聯波.低滲透砂巖儲層裂縫的形成與分布[M].北京:科學出版社,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] 康志勤,趙陽升,孟巧榮,等.油頁巖熱破裂規律顯微CT實驗研究[J].地球物理學報,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

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

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