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

姜黎明，孫建孟，劉學鋒，王海濤

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

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

姜黎明1，孫建孟1，劉學鋒2，王海濤1

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

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

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

0 引 言

地球物理資料的定性解釋需要了解巖石物理性質和巖石微觀結構的關系。隨著計算機技術的發展，可以根據巖石微觀結構信息重建反映巖石真實孔隙空間的三維數字巖心，通過在三維數字巖心上進行巖石物理實驗模擬，從而計算巖石的宏觀物理性質（如彈性模量、電導率等）。這種方法被稱為數字巖石物理或計算巖石物理，可以用于模擬孔隙結構、裂縫和流體對巖石聲學性質和電學性質的影響［1－3］。

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

1 數字巖心的建立

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

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

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

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

2 理論和數值模擬方法

2.1 流體替換：Gassmann理論

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

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

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

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

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

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

2.2 有限元法

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

3 數值模擬結果和討論

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

表1 巖石的彈性參數＊

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

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

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

圖4分別給出了巖石拉梅常數、泊松比、縱橫波速度、彈性模量及縱橫波速度比隨含氣飽和度的變化關系曲線。巖石的含氣飽和度對橫波速度影響不大，它的增加將引起橫波速度緩慢增加，這主要是因為含氣飽和度增加，導致了巖石密度降低。巖石的拉梅常數、泊松比、縱波速度、彈性模量和縱橫波速度比隨含氣飽和度的增加而減小，當含氣飽和度較低時，以上彈性參數隨含氣飽和度增加急劇降低。也就是說這些彈性參數在含水飽和度較高時，對于氣體含量非常敏感，當含有較多氣體時，它們對氣體含量不很敏感。

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

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

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

圖5 巖石彈性參數的相對變化率與含氣飽和度的關系

4 結 論

（1）把數值模擬結果與理論計算結果進行了比較，發現在整個含氣飽和度范圍內，有限元的數值模擬結果與GW理論計算結果基本吻合，證明了數值研究方法的可行性。

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

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

［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］.石油與天然氣地質，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］ 李生杰.巖性、孔隙及其流體變化對巖石彈性性質的影響［J］.石油與天然氣地質，2005，26（6）：760－764.

［8］ 姚軍，趙秀才，衣艷靜，等.數字巖心技術現狀及展望［J］.油氣地質與采收率，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高技術研究發展計劃項目（2006AA060105）

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