姜黎明，孫建孟，劉學(xué)鋒，王海濤

（1.中國石油大學(xué)地球科學(xué)與技術(shù)學(xué)院，山東青島266555；2.中國石油大學(xué)理學(xué)院，山東青島266555）

天然氣飽和度對巖石彈性參數(shù)影響的數(shù)值研究

姜黎明1，孫建孟1，劉學(xué)鋒2，王海濤1

（1.中國石油大學(xué)地球科學(xué)與技術(shù)學(xué)院，山東青島266555；2.中國石油大學(xué)理學(xué)院，山東青島266555）

基于三維數(shù)字巖心，對巖石的彈性模量、縱橫波速度、拉梅常數(shù)、泊松比等參數(shù)隨不同含氣飽和度變化規(guī)律進(jìn)行數(shù)值研究。研究結(jié)果表明巖石的含氣飽和度對橫波速度影響不大，含氣飽和度的增加將引起橫波速度的緩慢增加。巖石的拉梅常數(shù)、泊松比、縱波速度、彈性模量和縱橫波速度比隨含氣飽和度的增加而減小，當(dāng)含氣飽和度較低時，它們隨含氣飽和度增加急劇降低。將各彈性參數(shù)隨含氣飽和度的變化率進(jìn)行比較，發(fā)現(xiàn)對天然氣飽和度變化最敏感的參數(shù)是拉梅常數(shù)，其次是泊松比、體積模量和縱橫波速度比等。該研究有助于選取對天然氣飽和度變化敏感的參數(shù)計算儲層的含氣飽和度，從而提高計算精度。

三維數(shù)字巖心；彈性參數(shù)；有限元方法；Gassmann－Wood理論；橫波速度；拉梅常數(shù)

0 引 言

地球物理資料的定性解釋需要了解巖石物理性質(zhì)和巖石微觀結(jié)構(gòu)的關(guān)系。隨著計算機(jī)技術(shù)的發(fā)展，可以根據(jù)巖石微觀結(jié)構(gòu)信息重建反映巖石真實(shí)孔隙空間的三維數(shù)字巖心，通過在三維數(shù)字巖心上進(jìn)行巖石物理實(shí)驗(yàn)?zāi)M，從而計算巖石的宏觀物理性質(zhì)（如彈性模量、電導(dǎo)率等）。這種方法被稱為數(shù)字巖石物理或計算巖石物理，可以用于模擬孔隙結(jié)構(gòu)、裂縫和流體對巖石聲學(xué)性質(zhì)和電學(xué)性質(zhì)的影響［1－3］。

多孔巖石的彈性性質(zhì)主要受巖石的孔隙結(jié)構(gòu)、固體相以及流體性質(zhì)的影響。巖石孔隙中的流體及其變化對巖石的彈性性質(zhì)產(chǎn)生一定的影響。定量研究巖石彈性力學(xué)性質(zhì)隨所含流體變化的特征對于深入了解巖石物理性質(zhì)，特別是對油氣田開發(fā)具有重要意義［4］。許多學(xué)者通過實(shí)驗(yàn)觀測的方法研究了孔隙度、泥質(zhì)含量、流體性質(zhì)等對巖石聲波速度的影響［5－7］，由于飽和氣實(shí)驗(yàn)對實(shí)驗(yàn)設(shè)備要求比較高，一般很難實(shí)現(xiàn)。本文基于三維數(shù)字巖心，對巖石的彈性模量、縱橫波速度、拉梅常數(shù)、泊松比等參數(shù)隨不同含氣飽和度變化規(guī)律進(jìn)行了數(shù)值實(shí)驗(yàn)研究。

1 數(shù)字巖心的建立

建立數(shù)字巖心的方法主要有3種：切片組合法、X射線CT法和基于巖石二維圖像分析的重建法。其中，切片組合法需要花費(fèi)很長時間制備大量的巖心切片，因而極少被采用。X射線CT法是建立數(shù)字巖心最直接、最準(zhǔn)確的方法，但由于實(shí)驗(yàn)成本過高，限制了該項技術(shù)的廣泛應(yīng)用。X射線CT法重建三維圖像的主要步驟：首先通過X射線CT掃描得到與樣品X射線吸收系數(shù)有關(guān)的投影數(shù)據(jù)，然后通過重建算法（如卷積濾波法和迭代法），將其轉(zhuǎn)換為一系列的巖心橫截面二維圖像，將巖心的二維橫截面圖像組合便得到巖心的三維灰度圖像，最后通過對灰度圖像進(jìn)行圖像分割，得到研究所需的三維二值圖像即三維數(shù)組?；趲r石二維圖像分析的重建算法只需要極少量巖石切片的掃描圖像，其獲取方便、經(jīng)濟(jì)，因此該方法具有重要的應(yīng)用價值［8］。

基于巖石二維圖像分析的重建算法分為2類：隨機(jī)法和過程法。隨機(jī)法是在巖石二維圖像統(tǒng)計函數(shù)的約束下調(diào)整隨機(jī)場中孔隙和骨架的位置，使隨機(jī)場和二維圖像具有相同的統(tǒng)計特性，從而構(gòu)建三維數(shù)字巖心。高斯場法［9］和模擬退火算法［10］是應(yīng)用較為普遍的隨機(jī)法。隨機(jī)法重建數(shù)字巖心孔隙連通性較差，但適用性強(qiáng)。過程法重建數(shù)字巖心是結(jié)合巖石顆粒的粒徑分布，通過對沉積類巖石的形成過程（包括沉積、壓實(shí)和成巖作用）的模擬建立數(shù)字巖心［11－12］。過程法重建的數(shù)字巖心更接近真實(shí)巖心。

本文選用4塊通過X射線CT構(gòu)建的數(shù)字巖心，孔隙度分別為8%、13%、19%、21%，滲透率范圍是80～3 000mD＊＊非法定計量單位，1mD＝9.87×10－4μm2，下同（見圖1）。圖1中藍(lán)色代表巖石骨架，紅色代表孔隙。由于巖石物理性質(zhì)受巖石尺寸的影響，因此在研究中選了代表元體積（REV），通過多次測試驗(yàn)證表明當(dāng)數(shù)字巖心大小為200×200×200體素時，巖心的物理性質(zhì)（比如孔隙度、彈性模量）不再受巖心尺寸的影響。所以在研究中，選代表元體積為200×200×200體素。

圖1 X射線CT獲取的楓丹白露砂巖三維數(shù)字巖心

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

2.1 流體替換：Gassmann理論

Gassmann理論是研究流體飽和對巖石聲波速度影響的最常用理論。該理論把飽和流體孔隙巖石的彈性模量與干巖石骨架的彈性模量、固體顆粒彈性模量和孔隙流體的彈性模量聯(lián)系起來。流體飽和巖石的體積模量K為

如果孔隙空間被n種流體的混合物所填充，則混合流體的體積模量可以由Wood方程［5］給出

式中，xi和Ki分別為第i種流體相的體積含量和體積模量。對于氣和水2種復(fù)合流體，式（2）可寫為

式中，Sw為含水飽和度；Kw和Kg分別為水和氣的體積模量。把式（3）的流體模量Kf代入式（1）便可得到計算兩相流體飽和巖石體積模量的Gassmann－Wood（GW）方程

由此可計算單相和多相流體飽和巖石的有效彈性模量，進(jìn)而得到巖石的縱波速度vp、橫波速度vs、縱波阻抗Ip、拉梅常數(shù)λ及泊松比σ等。飽和流體巖石密度與孔隙度的關(guān)系可以表示為

式中，ρ0為固體顆粒的密度；ρf為孔隙流體的密度。巖石彈性參數(shù)可由以下彈性方程得到

2.2 有限元法

有限元法非常適合基于多孔巖石的微觀結(jié)構(gòu)預(yù)測巖石的有效彈性性質(zhì)［14］，并且可以用來處理具有任何體素的微觀結(jié)構(gòu)。給定固體組分和流體相的體積和剪切模量，就可以基于多孔巖石的三維微觀結(jié)構(gòu)用有限元的方法研究巖石的宏觀物理性質(zhì)，比如巖石的有效彈性模量。在三維數(shù)字巖心的基礎(chǔ)上，利用有限元法計算巖石彈性模量的理論基礎(chǔ)是變分原理。對于一個給定的數(shù)字巖心，施加一個宏觀應(yīng)變，利用快速共軛梯度的方法通過把體系的彈性自由能最小化來求取由這個外加應(yīng)變引起的平均應(yīng)力，從而求取巖石的有效彈性模量。首先把三維數(shù)字巖心離散化，離散后的每一個三維的體素都看成是1個三線性的有限元。實(shí)際計算的精確性取決于離散的誤差，分割的三維體素越多，誤差越小。然而，體素的多少又取決于計算機(jī)的內(nèi)存和計算速度（每個體素占內(nèi)存230B）?？紤]到計算機(jī)存儲和運(yùn)算能力的限制，本文研究中取200×200×200體素。

3 數(shù)值模擬結(jié)果和討論

本文基于三維數(shù)字巖心利用有限元的方法模擬了含氣飽和度對巖石彈性性質(zhì)的影響，在數(shù)值模擬中，巖石骨架和孔隙流體的彈性參數(shù)取值見表1。

表1 巖石的彈性參數(shù)＊

對于兩相流體飽和巖石，利用格子玻爾茲曼方法［15］確定三維數(shù)字巖心孔隙空間的流體分布。圖2是孔隙度為19%的楓丹白露砂巖（親水巖石）在不同含氣飽和度下流體分布圖，其中藍(lán)色代表巖石骨架，紅色代表天然氣，綠色代表水。

圖2 不同含氣飽和度下的流體分布（親水巖石）

在不同含氣飽和度條件下，孔隙空間的流體分布不同，巖石的彈性模量會受到流體分布和流體性質(zhì)的影響。首先利用格子玻爾茲曼方法確定了不同含氣飽和度下孔隙空間的流體分布，然后基于三維數(shù)字巖心利用有限元方法計算了巖石的體積模量和剪切模量，研究了含氣飽和度和彈性模量的關(guān)系并和GW理論進(jìn)行了比較，具體計算結(jié)果見圖3。從圖3可見，在整個含氣飽和度范圍內(nèi)，有限元的數(shù)值模擬結(jié)果與GW理論計算結(jié)果基本吻合，這說明基于三維數(shù)字巖心利用有限元方法計算兩相流體飽和巖石的彈性模量是可行的。從圖3還可以看出隨含氣飽和度的增加，巖石的體積模量減小，這是因?yàn)闅夂康脑黾訉?dǎo)致了巖石的壓縮性增大，巖石的剪切模量不受流體變化的影響。

圖4分別給出了巖石拉梅常數(shù)、泊松比、縱橫波速度、彈性模量及縱橫波速度比隨含氣飽和度的變化關(guān)系曲線。巖石的含氣飽和度對橫波速度影響不大，它的增加將引起橫波速度緩慢增加，這主要是因?yàn)楹瑲怙柡投仍黾?，?dǎo)致了巖石密度降低。巖石的拉梅常數(shù)、泊松比、縱波速度、彈性模量和縱橫波速度比隨含氣飽和度的增加而減小，當(dāng)含氣飽和度較低時，以上彈性參數(shù)隨含氣飽和度增加急劇降低。也就是說這些彈性參數(shù)在含水飽和度較高時，對于氣體含量非常敏感，當(dāng)含有較多氣體時，它們對氣體含量不很敏感。

對于天然氣儲層，氣體的儲集會引起整個儲集層的彈性模量、拉梅常數(shù)、泊松比、縱波速度、彈性模量及縱橫波速度比明顯降低，而橫波速度略有升高。為了研究各彈性參數(shù)對含氣飽和度變化的敏感性，定義彈性參數(shù)隨含氣飽和度相對變化率

式中，ypg和yg分別表示部分飽和氣和完全飽和氣巖石的彈性參數(shù)。

圖5將各彈性參數(shù)隨含氣飽和度的變化率進(jìn)行了比較。對含氣飽和度變化最敏感的參數(shù)是拉梅常數(shù)，其次是泊松比、體積模量和縱橫波速度比，特別是當(dāng)含氣飽和度較小的時候，含氣飽和度的變化能引起這些彈性參數(shù)的明顯變化。在結(jié)合測井資料計算巖石含氣飽和度時，可以選擇對含氣飽和度變化反應(yīng)敏感的參數(shù)計算。首先要選擇拉梅常數(shù)，其次是泊松比、體積模量和縱橫波速度比，這樣可以減小計算誤差，提高計算精度。

圖5 巖石彈性參數(shù)的相對變化率與含氣飽和度的關(guān)系

4 結(jié) 論

（1）把數(shù)值模擬結(jié)果與理論計算結(jié)果進(jìn)行了比較，發(fā)現(xiàn)在整個含氣飽和度范圍內(nèi)，有限元的數(shù)值模擬結(jié)果與GW理論計算結(jié)果基本吻合，證明了數(shù)值研究方法的可行性。

（2）巖石的含氣飽和度對橫波速度影響不大，它的增加將引起橫波速度緩慢增加；巖石的拉梅常數(shù)、泊松比、縱波速度、彈性模量和縱橫波速度比隨含氣飽和度的增加而減小，當(dāng)含氣飽和度較低時，它們隨含氣飽和度增加急劇降低，說明這些彈性參數(shù)在含水飽和度較高時，對于氣體含量非常敏感，當(dāng)含有較多氣體時，它們對氣體含量不很敏感。

（3）將各彈性參數(shù)隨不同含氣飽和度的變化率進(jìn)行了比較，發(fā)現(xiàn)對含氣飽和度變化最敏感的參數(shù)是拉梅常數(shù)，其次是泊松比、體積模量和縱橫波速度比。該研究有助于選取對含氣飽和度變化敏感的參數(shù)計算儲層的含氣飽和度，從而提高計算精度。

［1］ Meille S，Garboczi E J.Linear Elastic Properties of 2D and 3DModels of Porous Materials Made from Elongated Objects［J］.Modelling Simul.Mater.Sci.Eng，2001，9（1）：371－390.

［2］ Knackstedt M A，Arns C H，Sheppard A P，et al.Archie’s Exponents in Complex Lithologies Derived from 3DDigital Core Analysis［C］∥The SPWLA 48th Annual Logging Symposium，paper UU，Austin：2007.

［3］ Jiang Liming，Sun Jianmeng，Liu Xuefeng，et al.Study of Different Factors Affecting the Electrical Properties of Natural Gas Reservoir Rocks Based on Digital Cores［J］.Journal of Geophysics and Engineering，2011，8（2）：366－371.

［4］ 杜韞華，錢凱，張守鵬，等.中國天然氣儲層的巖石、古地理類型與勘探方向［J］.石油與天然氣地質(zhì)，1999，20（2）：140－143.

［5］ King M S，Marsden J R，Dennis J W.Biot Dispersion for P－and S－wave Velocity in Partially and Fully Saturated Sandstone［J］.Geophysical Prospecting，2000，48（6）：1075－1089.

［6］ Clark V A.The Effect of Oil under In－situ Conditions on the Seismic Properties of Rocks［J］.Geophysics，1992，57（7）：894－901.

［7］ 李生杰.巖性、孔隙及其流體變化對巖石彈性性質(zhì)的影響［J］.石油與天然氣地質(zhì)，2005，26（6）：760－764.

［8］ 姚軍，趙秀才，衣艷靜，等.數(shù)字巖心技術(shù)現(xiàn)狀及展望［J］.油氣地質(zhì)與采收率，2005，12（6）：52－54.

［9］ Adler P M，J acquin C G，Quiblier J A.Flow in Simulated Porous Media［J］.International Journal of Multiphase Flow，1990，16（4）：691－712.

［10］Hazlett R D.Statistical Characterization and Stochastic Modeling of Pore Networks in Relation to Fluid Flow［J］.Mathematical Geology，1997，29（6）：801－822.

［11］Oren P E，Bakke S.Process Based Reconstruction of Sandstones and Prediction of Transport Properties［J］.Transport in Porous Media，2002，46：311－343.

［12］Oren P E，Bakke S.Reconstruction of Berea Sandstone and Pore－scale Modelling of Wettability Effects［J］.Journal of Petroleum Science and Engineering，2003，39：177－199.

［13］Krief M，Garat J，Stellingwerff J，et al.A Petrophysical Interpretation Using the Velocities of P and S Waves（Full－wave Sonic）［C］∥The 12th SPWLA Logging Symposium，Paper HH，Paris：1989.

［14］Garboczi E J，Day A R.Algorithm for Computing the Effective Linear Elastic Properties of Heterogeneous Materials：Three－dimensional Results for Composites with Equal Phase Poisson Ratios［J］.Journal of the Mechanics and Physics of Solids，1995，43（9）：1349－1362.

［15］Pan C，Hilpert M，Miller C T.Lattice－Boltzmann Simulation of Two－phase Flow in Porous Media［J］.Water Resources Research，2004，40：W01501.

Numerical Study of the Effect of Natural Gas Saturation on the Reservoir Rocks’Elastic Parameters

JIANG Liming1，SUN Jianmeng1，LIU Xuefeng2，WANG Haitao1
（1.School of Geo－sciences，China University of Petroleum，Qingdao，Shandong 266555，China；2.College of Science，China University of Petroleum，Qingdao，Shandong 266555，China）

The variation law of the effect of natural gas saturation on the reservoir rocks’elastic parameters，such as elastic moduli，P－wave and S－wave velocities，Laméconstant，Poisson’s ratio，is studied based on 3Ddigital cores.The results show that gas saturation of reservoir rocks has little influence on S－wave velocity and the S－wave velocity increases slowly with the increase of gas saturation.The elastic parameters，such as Laméconstant，Poisson’s ratio，P－wave velocity，the ratio of P－wave velocity to S－wave velocity，decrease as the gas saturation increases.Moreover，they are changed sharply when the sandstone contains very small amount of gas.After comparing the change ratios of the elastic parameters with the variation of gas saturation，it is found the Laméconstant is most sensitive to the variation of gas saturation，followed by Poisson’s ratio，bulk modulus and the ratio of P－wave velocity to S－wave velocity.Through this study，the parameter which is sensitive to change of gas saturation can be used to calculate gas saturation of reservoir rocks in order to improve the calculation accuracy.

3Ddigital core，elastic parameter，finite element method，Gassmann－Wood theory，S－wave velocity，Laméconstant

TE19

A

2011－05－13 本文編輯 王小寧）

1004－1338（2012）03－0239－05

國家863高技術(shù)研究發(fā)展計劃項目（2006AA060105）

姜黎明，男，1983年生，博士研究生，從事巖石物理數(shù)值模擬及測井解釋研究。