基于過程法構建二維數(shù)字巖心

盛金昌,何光宇,周治榮,姚德生,詹美禮,羅玉龍

(1.河海大學水利水電學院,江蘇南京 210098;2.中國水電顧問集團成都勘測設計研究院,四川成都 610072;3.廣西電力工業(yè)勘察設計研究院,廣西南寧 530023;4.中水北方勘測設計研究有限責任公司,天津 300222)

達西定律是研究多孔介質(zhì)滲流的基礎,它是建立在宏觀試驗的基礎上,無法對細觀孔隙尺度的滲流進行描述,而石油開采等領域迫切需要從細觀孔隙角度出發(fā)研究滲流問題.目前多孔介質(zhì)滲流的一些細觀機理可以通過物理試驗來定性研究,但缺少科學的理論和技術來進行定量描述.目前國內(nèi)外已有不少學者開展了這方面的研究,而構建細觀孔隙網(wǎng)絡模型,即數(shù)字巖心是其中一項基礎性的工作.

自1956年Fatt[1]首次引入真正意義的孔隙網(wǎng)絡模型以來,隨著計算機和圖像分析技術的發(fā)展,很多學者[2-6]對孔隙網(wǎng)絡模型進行了研究.構建數(shù)字巖心有3類方法,即切片組合法、基于切片分析的圖像重建法和X射線立體成像法.基于切片分析的圖像重建法是基于二維巖心切片掃描圖像的構建方法,只需少量的巖心切片,而巖心切片在石油開采中常被用來進行地質(zhì)評價,較易獲得,故該方法被大多研究所采納.研究[7-8]表明,將傳統(tǒng)的兩點相關函數(shù)(隨機法)用來描述多孔介質(zhì)結構,尤其是孔隙之間的連通關系是不夠充分的.本文在前人研究的基礎上,基于過程法編程實現(xiàn)了二維數(shù)字巖心的構建,并對結果進行了可視化處理.

1 過程法構建數(shù)字巖心原理

過程法是Bakke等[4]于1997年提出的,它是在對二維巖心切片分析的基礎上,結合巖心的顆粒粒徑分布特性,通過對沉積類巖石3個形成過程(沉積、壓實和成巖)的模擬來構建數(shù)字巖心的.文獻[8]指出:建立準確、可靠的數(shù)字巖心所需的輸入數(shù)據(jù)可以從標準的巖石切片中獲得,這些輸入數(shù)據(jù)包括粒度分布、孔隙度、壓實程度的視覺估計、黏土材料、膠結物的類型和數(shù)量等.

1.1 模擬沉積過程的算法

基于過程法構建數(shù)字巖心的沉積模擬主要采用顆粒隨機堆積模型(分為連續(xù)堆積法和重整化群法2種方法)[9],本文采用連續(xù)堆積法.

1.1.1 基本假設

巖石顆粒沉積過程十分復雜,為簡化算法,在不影響沉積規(guī)律的前提下,進行如下假設[4,10-11]:(a)所有顆粒均為球體;(b)顆粒初始下落位置是隨機的,且各點幾率均等;(c)考慮沉積結構的穩(wěn)定性,忽略由于重力壓實引起的顆粒形變和結構重組效應;(d)顆粒下落時沿著重力勢能梯度下降最大的方向運動;(e)顆粒在下降過程中只受重力作用,忽略側向作用力;(f)顆粒下落達到相對平衡位置以后,其位置保持不變;(g)下降顆粒與已沉積穩(wěn)定顆粒發(fā)生碰撞時不發(fā)生彈跳;(h)在下落顆粒達到平衡位置以后,下一顆粒才開始下落;(i)二維模擬時沉積區(qū)域為一矩形區(qū)域.

1.1.2 沉積過程模擬

沉積過程模擬得到的是巖心顆粒的松散堆積體.模擬步驟為:(a)隨機生成第1層,作為表層顆粒鏈.(b)在粒度分布曲線上隨機抽取1個顆粒,半徑為Rd,在沉積區(qū)域(0≤x≤Lxmax,y≥0),根據(jù)式

隨機生成下落顆粒的x坐標(ξ為隨機產(chǎn)生的數(shù)),其y坐標取為常數(shù)y0.(c)根據(jù)xd與表層顆粒鏈中顆粒的x坐標之間的關系,確定與下落顆粒碰撞的表層鏈顆粒.(d)計算下落顆粒中心坐標,判斷其穩(wěn)定性,如果不穩(wěn)定,向左或者向右滾動1個顆粒,重復第(c)步,反復操作,直到顆粒達到穩(wěn)定位置.二維情況下顆粒最終平衡位置有4種,如圖1所示.(e)記錄下落顆粒的中心坐標及其半徑.(f)根據(jù)下落顆粒的中心位置與目前表層鏈顆粒中心位置的關系,更新表層顆粒鏈.(g)判斷所需顆粒數(shù)是否已沉降完畢,如果是,則結束沉積模擬,反之繼續(xù)第(b)～(g)步,直到所有顆粒沉降完成.

圖1 顆粒二維堆積過程中的4種可能平衡位置Fig.1 Four possible equilibrium positions of 2D particle accumu lation

1.2 模擬壓實過程的算法

巖石的壓實過程是指受到來自上層巖石的重力以及周圍巖石、流體的接觸壓力和摩擦力等作用后,巖石顆粒發(fā)生平動或轉(zhuǎn)動2種運動,是一個非常復雜的過程.文獻[4,7]認為,下覆巖層受到上覆巖層壓實后,其總體積和孔隙度下降,因此在算法實現(xiàn)時,應引進一個壓實因子來描述孔隙減少程度;Jin等[10]則從動力過程來考慮巖石的壓實過程.本文進行了如下假設:顆粒的幾何形狀不發(fā)生變化,無側向位移,顆粒不發(fā)生轉(zhuǎn)動,顆粒無彈性變形,只通過改變各個顆粒的垂直坐標來實現(xiàn)巖石的壓實過程.對于不同的巖石壓實程度,通過引入壓實因子 λ,使構建的數(shù)字巖心與真實巖心相符合.λ的取值范圍是[0,1],λ越大表示巖石受壓縮程度越大,反之則越小.

二維空間壓實前后顆粒相對位置的關系為

式中:y——壓實后的坐標;y0——壓實前的坐標.

1.3 模擬成巖過程的算法

成巖過程是指未固結的疏松沉積物變?yōu)閳杂矌r石的過程.這一過程十分復雜,包含了物理、化學、力學和生物化學等各種作用,這些作用從沉積到變質(zhì)一直影響著沉積物的各種特性.本文對成巖過程中的主要部分——礦物質(zhì)的膠結作用進行研究.膠結作用是礦物質(zhì)沿巖石顆粒表面結晶和沉淀的過程,最普通的礦物質(zhì)膠結類型是石英膠結物.

Roberts等[12]首次提出了一致增長顆粒固結模型.為了提高該模型的適用性,Schwartz等[13]提出了一個模擬顆粒膠結作用的交替增長算法:

交替增長算法認為增長速率應該依賴于增長方向.對于每一個增長方向r,定義 Δ為初始球(R)表面與Wigner-Seitz方塊邊界之間的距離,通過引入增生指數(shù)β來控制顆粒中心與孔隙和顆粒界面之間距離的增長速率,參數(shù) α用于控制增長數(shù)量.

Jin等[10]提出了一種改進的顆粒生長算法,該算法與文獻[12-13]的算法不同,考慮了硅石膠結增生率和巖石顆粒大小對巖石顆粒增長的影響:

式中:Δ(r)——礦物質(zhì)膠結增生從巖石顆粒中心沿r方向的增量;ˉR——所有顆粒的平均半徑;R0——原始顆粒的半徑;l(r)——沿r方向原球形顆粒表面到該顆粒對應的多邊形表面或者與其相鄰顆粒表面之間的距離;L(r)——從顆粒中心到孔隙膠結表面的距離;κ——參數(shù),控制孔隙度下降數(shù)量.

根據(jù)文獻[12-13]的原理,筆者編制程序模擬了巖石的成巖過程.圖2給出了α分別為-2,0和2時巖石顆粒石英膠結增生圖(圖中原巖石顆粒為白色,石英膠結物為灰色,孔隙空間為黑色).對比圖2可以看出:(a)若α＜0,膠結物優(yōu)先在大顆粒的表面生長,多數(shù)的膠結物沉淀在大顆粒的表面;(b)若 α=0,膠結物均勻地沉淀于所有巖石顆粒的表面;(c)若α＞0,膠結物優(yōu)先在小顆粒表面生長,多數(shù)的膠結物沉淀在小顆粒的表面.從物理觀點分析來看,參數(shù)α取正值似乎更合理并與現(xiàn)實相符,因為小顆粒比大顆粒具有更大的比表面積.

圖2 不同α條件下巖石顆粒的石英膠結Fig.2 Quartz cementation around rock particles under different values of α

圖3為β不同時巖石顆粒石英膠結增生圖.β可以取任意數(shù)(負數(shù)、零或正數(shù)),本文β分別取-1,0,1來代表負數(shù)、零或正數(shù)對巖石顆粒石英膠結增生的影響.從圖3可以看出:(a)若 β＜0,石英膠結物優(yōu)先在喉道中生長,優(yōu)先沉淀在狹窄的巖石顆粒接觸區(qū)域,降低了孔隙的面容比.當孔隙空間變得不連貫時,這種生長傾向于保持大孔隙的連貫性.(b)若β=0,石英膠結物在孔隙和喉道中均勻生長,均勻地沉淀于所有顆粒表面.(c)若β＞0,膠結物優(yōu)先在孔隙中生長,多數(shù)膠結物沿著 l(r)最大方向沉淀.這種類型的膠結傾向于增加原孔隙空間的孔隙面容比,并且當孔隙度非常低時,會使孔隙空間相互連接.總之,β影響著膠結增生的幾何趨勢.

圖3 不同β時巖石顆粒的石英膠結Fig.3 Quartz cementation around rock particles under different values of β

2 構建二維數(shù)字巖心實例及可視化

2.1 顆粒半徑分布曲線

應用計算機圖像處理技術對高分辨率的二維巖心薄切片圖像進行分析處理,可得試驗樣本建模所需參數(shù).假設所有的顆粒都是圓(球)形顆粒,通過光散射原理計算可構造出粒度分布曲線.本文選取240個巖石顆粒,假設最小巖石顆粒的半徑為10μ m,最大巖石顆粒半徑為600μ m,經(jīng)過處理,得到如圖4所示的一條巖石顆粒半徑及其編號分布曲線.

2.2 沉積模擬

根據(jù)顆粒隨機堆積過程基本思想,筆者采用MATLAB軟件編制了巖石顆粒在二維空間中隨機沉積模擬程序,并用MATLAB軟件實現(xiàn)了隨機沉積結果的可視化.在圖4所示的曲線上選取81個顆粒進行1次沉積模擬(每個編號對應1個顆粒半徑),可視化結果如圖5所示.因為隨機模擬,每次模擬都會呈現(xiàn)不同的結果.

2.3 壓實模擬結果

將沉積模擬結果作為初始輸入數(shù)據(jù)進行壓實算法模擬,得到的巖石壓實模擬結果如圖6所示,其中λ分別為0.05,0.1,0.20,隨著λ的增大,孔隙不斷縮小.

圖4 巖心顆粒粒度與編號分布關系Fig.4 Relationship between rock particle size and particle number

圖5 二維沉積模擬結果(白色為巖石顆粒)Fig.5 Results of 2D deposition simu lation(where white denotes rock particles)

圖6 二維壓實模擬不同λ的壓實對比(白色為巖石顆粒)Fig.6 Results of 2D compaction simu lation with different values of λ(where white denotes rock particles)

2.4 成巖模擬結果

成巖模擬采用改進的交替增長算法,在壓實模擬數(shù)據(jù)的基礎上進行,在此取 λ=0,α=1,β=1,κ=0.2來模擬巖心顆粒石英膠結增生,結果如圖7所示.

圖7 巖石顆粒的石英膠結增生結果Fig.7 Overgrow th results of quartz cementation

3 結 語

巖石類多孔介質(zhì)的滲流特性非常復雜,以往的研究多集中于宏觀意義上的滲流問題.目前,盡管可以通過實驗來研究多孔介質(zhì)中細觀滲流的一些機理,但還缺少必要的手段和科學合理的方法對細觀滲流進行定量描述.在孔隙級別上研究多孔介質(zhì)細觀滲流問題,建立一套可以準確描述細觀滲流機理的理論體系,具有非常重要的理論意義和應用價值,而數(shù)字巖心則是開展細觀滲流機理研究的基礎.

本文在前人工作的基礎上,基于過程法編程完成了二維數(shù)字巖心的構建工作,并實現(xiàn)了二維巖心的可視化.在巖心構建過程中,參數(shù)α控制著石英膠結物沉淀的優(yōu)先位置,增生指數(shù) β控制著石英膠結物優(yōu)先增生的幾何趨勢,κ則控制著孔隙度下降的數(shù)量,壓實因子λ反映上覆巖層壓力,壓實因子越大,孔隙壓縮程度越高.

[1]FATT I.The network model of porous media:I capillary pressure characteristics[J].Transactions of AIME,1956,207(1):144-159.

[2]JOSH I M.A class of stochastic models for porous media[D].Law rence Kansas:University of Kansas,1974.

[3]HAZLETT R D.Statistical characterization and stochastic modeling of pore networks in relation to fluid flow[J].Mathematical Geology,1997,29(6):801-822.

[4]BAKKE S,?REN P E.3-D pore-scale modeling of sandstones and flow simulations in the pore networks[J].SPE Journal,1997,2:136-149.

[5]STREBELLE S.Sequential simulation drawing structures from training images[D].Palo Alto:Stanford University,2000.

[6]OKABE H,BLUNT M J.Prediction of permeability for porous media reconstructed using multip le-point statistics[J].Physical Review E,2004,70:106-135.

[7]BISWAL B,MANSWARTH C,H ILFER R,et al.Quantitative analysis of experimental and synthetic microstructures for sedimentary rocks[J].Physics A,1999,273:452-475.

[8]?REN P E,BAKKE S.Process based reconstruction of sandstones and predictions of transport properties[J].Transport in Porous Media,2002,46(2):311-343.

[9]HE D,EKERE N N.Computer simulation of powder compaction of spherical particles[J].Journal of MaterialsScience Litters,1998,17:1723-1725.

[10]JIN G,PATZEK T W,SILIN D Y B.Physics-based reconstruction of sedimentary rocks[R].Richardson:SPE,2003:305-318.

[11]薄啟煒,董長銀,張琪,等.礫石充填層孔喉結構可視化模擬[J].石油勘探與開發(fā),2003,30(4):108-110.(BO Qi-wei,DONG Chang-yin,ZHANG Qi,et al.Pore filling gravel layer structure visualization simu lation[J].Petroleum Exploration and Development,2003,30(4):108-110.(in Chinese))

[12]ROBERTS J N,SCHWARTZ L M.Grain consolidation and electrical conductivity in porous media[J].Physical Rev B,1985,31(9):5990-5997.

[13]SCHWAR TZ L M,KIMMINAU S.Analysis of electrical conduction in the grain consolidation model[J].Geophysics,1987,52(10):1402-1411.