王

（武漢紡織大學(xué) 數(shù)學(xué)與計(jì)算機(jī)學(xué)院，湖北 武漢 430073）

隨機(jī)多孔介質(zhì)中的不混溶驅(qū)替

（武漢紡織大學(xué) 數(shù)學(xué)與計(jì)算機(jī)學(xué)院，湖北 武漢 430073）

對(duì)隨機(jī)多孔介質(zhì)中粘滯指進(jìn)的分形性質(zhì)進(jìn)行研究，建立正方網(wǎng)格來(lái)模擬多孔介質(zhì)中潤(rùn)濕流體的侵入和流動(dòng)，運(yùn)用決定論模型進(jìn)行研究。結(jié)果表明，多孔介質(zhì)的幾何拓?fù)鋸?qiáng)烈地影響了粘滯指進(jìn)的結(jié)構(gòu)和驅(qū)替過(guò)程。粘滯指進(jìn)圖像掃及面積隨著迭代次數(shù)n的增加而增大。增加迭代次數(shù)n和網(wǎng)格尺寸，會(huì)導(dǎo)致驅(qū)掃效率E的增大。

多孔介質(zhì)；逾滲集團(tuán)；粘滯指進(jìn)；驅(qū)掃效率

隨機(jī)多孔介質(zhì)中流體的輸運(yùn)在許多工程領(lǐng)域都有重要的應(yīng)用。例如土壤中的有害廢物的擴(kuò)散，石油工程中油的移動(dòng)，色譜的分離和催化過(guò)程等。關(guān)于粘滯指進(jìn)問(wèn)題的第一個(gè)理論分析是由 Saffman 和 Taylor[1]針對(duì)Hele-Shaw[2]裝置進(jìn)行的。均勻網(wǎng)格中的兩相流已有大量研究結(jié)果[3]。本文把多孔介質(zhì)用正方網(wǎng)格中逾滲集團(tuán)表示，運(yùn)用決定論模型[3]對(duì)隨機(jī)多孔介質(zhì)中的粘滯指進(jìn)分形問(wèn)題進(jìn)行研究。

1 模擬方法

我們視多孔介質(zhì)為正方網(wǎng)格，格點(diǎn)表示孔隙，喉管表示通道，稱為鍵。如圖1所示，當(dāng)格點(diǎn)取值為0時(shí)(座(i+1,j))，該點(diǎn)不連通，此時(shí)與之相連的所有的鍵均完全封閉(圖1中與座(i,j)相連的右鍵，或圖2中的a或f鍵)，即導(dǎo)流率g(i,j),(i-1,j)=0；當(dāng)格點(diǎn)取值為1時(shí)(圖1中的座(i,j),(i-1,j),(i,j-1)和(i,j+1))，除其近鄰為0的座外，與之相連的所有的鍵均完全導(dǎo)通，g(i,j),(i-1,j)≠0等等。

對(duì)于后一類鍵又有三種情況：已完全被驅(qū)替流體所占據(jù)(圖1中與座(i,j)相連的左鍵，或圖2中的d鍵)，該鍵兩端點(diǎn)的壓力p均為1；已部分被驅(qū)體流體所占據(jù)(圖1中與座(i,j)相連的上鍵，或圖2中的c鍵)，該鍵一端點(diǎn)的壓力p為1，另一端點(diǎn)的壓力p小于1；完全沒(méi)有被驅(qū)體流體所占據(jù)(圖1中與座(i,j)相連的下鍵，或圖2中的b鍵)，該鍵上完全充滿被驅(qū)替流體，且一端點(diǎn)的壓力p為1，另一端點(diǎn)的壓力p小于1。還有一種鍵的狀態(tài)是，該鍵本身導(dǎo)通，但其兩端的格點(diǎn)的壓力均小于1，即驅(qū)替流體尚未流到這些格點(diǎn)(或圖2中的e鍵)。

圖1 網(wǎng)格中流體流動(dòng)狀態(tài)示意圖

圖2 網(wǎng)格中座鍵狀態(tài)示意圖

在201×201正方網(wǎng)格中，給定一個(gè)逾滲概率P，若P≥0.59，則存在一個(gè)跨越集團(tuán)；當(dāng)P=1時(shí)，格點(diǎn)的占有率為100%，對(duì)于不同的P對(duì)應(yīng)于不同連通性的孔隙介質(zhì)。對(duì)于不屬于跨越集團(tuán)的格點(diǎn)和集團(tuán)，它對(duì)應(yīng)沉積巖的不連通區(qū)域。我們的模擬僅在跨越集團(tuán)中進(jìn)行。為此，我們建立孔隙介質(zhì)和兩相流的模型如下：不同逾滲概率P的逾滲集團(tuán)對(duì)應(yīng)不同的孔隙介質(zhì)，視格點(diǎn)為孔隙，鍵為喉管，且喉管的半徑是非均勻分布的：r→(r-rmin)/(rmax-rmin)，rmin和rmax分別是喉管半徑的最小和最大值。在參考文獻(xiàn)[4,5,6]曾運(yùn)用過(guò)類似的模型。

因此，對(duì)于任意一個(gè)喉管i，其半徑ri在0-1之間。例如圖3(a)是一個(gè)5×5的正方網(wǎng)格，假設(shè)網(wǎng)格中心為流體注入點(diǎn)。驅(qū)替規(guī)則如下：每一步，依照逐次超松弛技術(shù)只有一個(gè)喉管充滿驅(qū)替流體與被驅(qū)替流體。圖 3(b)顯示了圖 3(a)第一次驅(qū)替的過(guò)程。假設(shè)當(dāng)驅(qū)替流體到達(dá)逾滲集團(tuán)的邊界時(shí)凝集。為了模擬凝集的破壞性，我們假設(shè)每個(gè)被注入流體的喉管的尺寸隨著下面的規(guī)則增加[7]: ri→ri+∈(1-ri)，其中∈表示取0-1之間的隨機(jī)數(shù)。這樣做使得ri(t+1)總是在[ri(t), 1]之間變化。在某種意義上也就是說(shuō)，喉管的半徑隨著驅(qū)替流體凝集的密度而變化。

例如圖 3(c)表明了 3(b)驅(qū)替后被破壞的多孔介質(zhì)。充滿凝集流體的喉管半徑遵循下面的公式ri→ri+∈(1-ri)，并且介質(zhì)完全干縮，這時(shí)候新的驅(qū)替流體進(jìn)入。依照對(duì)r的一般假設(shè)，ri(t)的最大值是1，跟rmax的取值無(wú)關(guān)。比如，如果喉管的半徑為最大值rmax，那么經(jīng)過(guò)凝集破壞后的半徑增加到rmax+ ε (ε ＞ 0)。這時(shí)，r→[(rmax+ε )-rmin)]/[ (rmax+ε )-rmin)。我們研究當(dāng)P從0.6到1變化時(shí)201×201網(wǎng)格中的逾滲，迭代次數(shù)直到n=20。

圖3 驅(qū)替凝集過(guò)程示意圖

由Poiseuille定律，通過(guò)兩個(gè)相鄰格點(diǎn)i，j之間鍵中的流體流量為

我們采用逐次超松弛技術(shù)從方程(1)和(2)中解出壓力場(chǎng)

這里，超松弛系數(shù)ω為1.66。邊界條件為：在注入點(diǎn)p = 1，在網(wǎng)格外周界p = 0。在確定了壓力之后，驅(qū)替流體的前鋒移動(dòng)一段位移

在一個(gè)界面相鄰的格點(diǎn)，我們選擇各鍵中驅(qū)替前沿最先到達(dá)下一個(gè)格點(diǎn)所需時(shí)間（即各鍵驅(qū)替時(shí)間中最短的那一個(gè)）為時(shí)間步長(zhǎng)Δt。然后，按照相應(yīng)的時(shí)間步長(zhǎng)Δt移動(dòng)兩相流的交界面并引入新的交界面，計(jì)算壓力場(chǎng)，重復(fù)整個(gè)過(guò)程，直到到達(dá)外部邊界為止。這種方法的主要優(yōu)點(diǎn)是我們可以對(duì)驅(qū)替過(guò)程進(jìn)行研究，也被稱為決定論模型[8]。

2 多孔介質(zhì)的粘滯指進(jìn)生長(zhǎng)

圖4表明了在5種迭代次數(shù)下（n=1,2,5,10,20），逾滲集團(tuán)中粘滯指進(jìn)的圖像（P = 0.8，M = 10）。n = 1時(shí)的粘滯指進(jìn)圖像已經(jīng)在參考文獻(xiàn)[9]中研究過(guò)了。當(dāng)n ＞ 1時(shí)，驅(qū)替流體掃及范圍隨著參數(shù)n的增加而增大，而當(dāng)n ≥ 10時(shí)，粘滯指進(jìn)圖像的形狀幾乎沒(méi)有改變了。而且，改變迭代次數(shù)不會(huì)影響逾滲集團(tuán)中粘滯指進(jìn)的各向異性。這個(gè)結(jié)果說(shuō)明多孔介質(zhì)的幾何拓?fù)浣Y(jié)構(gòu)很強(qiáng)地影響了粘滯指進(jìn)的形狀和結(jié)構(gòu)。

圖4 逾滲集團(tuán)（201×201網(wǎng)格）中粘滯指進(jìn)圖像

圖5 分形維數(shù)D與逾滲概率P的函數(shù)關(guān)系圖

為了研究多孔介質(zhì)的幾何拓?fù)浣Y(jié)構(gòu)對(duì)粘滯指進(jìn)圖像的形狀和結(jié)構(gòu)的影響，本文得到了粘滯比 M=10，不同逾滲概率P及不同迭代參數(shù)n下的粘滯指進(jìn)圖像，并計(jì)算出其分形維數(shù)D。首先測(cè)量不同網(wǎng)格L×L ，不同迭代次數(shù)n，以及不同概率P下粘滯指進(jìn)圖像中驅(qū)替流體掃及面積S，假定[4]: S ～ LD。在圖5顯示了粘滯指進(jìn)圖像的分形維數(shù)D與逾滲概率P之間的關(guān)系。當(dāng)P=1時(shí)，無(wú)論迭代幾次，分形維數(shù)D=2。也就是說(shuō)，在二維正方網(wǎng)格中粘滯指進(jìn)圖像的結(jié)構(gòu)依賴于有限的粘滯比M。當(dāng)n=10及n=20時(shí)，粘滯指進(jìn)圖像的分形維數(shù)D和P無(wú)關(guān)，也就是說(shuō)，當(dāng)n ≥ 10時(shí)，粘滯指進(jìn)的結(jié)構(gòu)是穩(wěn)定的。由圖像可以看出，分形維數(shù)D隨著P的增加而增加，與迭代次數(shù)n無(wú)關(guān)。這表明，多孔介質(zhì)的幾何拓?fù)浣Y(jié)構(gòu)對(duì)粘滯指進(jìn)有很大影響。

3 驅(qū)掃效率

研究多孔介質(zhì)粘滯指進(jìn)現(xiàn)象的一個(gè)特別有實(shí)用價(jià)值的方面就是驅(qū)掃效率E，它是預(yù)測(cè)驅(qū)替性能的一個(gè)重要參數(shù)。在現(xiàn)場(chǎng)驅(qū)替過(guò)程中，體積驅(qū)掃效率定義為：E=驅(qū)掃體積/總體積。由于現(xiàn)在我們研究的是隨機(jī)多孔介質(zhì)中的粘滯指進(jìn)，因此，此時(shí)的效率應(yīng)為面積驅(qū)掃效率，設(shè)掃及面積為Ae掃及區(qū)域內(nèi)多孔介質(zhì)的面積為As，于是有[8]

式中L為網(wǎng)格尺寸。該式表明分維D能夠被合理地作為評(píng)價(jià)驅(qū)掃效率和石油儲(chǔ)量的有用參數(shù)。

圖6顯示了驅(qū)掃效率E與網(wǎng)格尺寸L之間的函數(shù)關(guān)系。由圖可以看出增加迭代次數(shù)n和網(wǎng)格尺寸L，驅(qū)掃效率E隨之減小，當(dāng)L達(dá)到最大值時(shí)E有最小值。圖中每條E～L曲線都存在有兩個(gè)線性區(qū)域，分別表示粘滯指進(jìn)的凝集區(qū)域和驅(qū)替區(qū)域[10]。

圖6 驅(qū)掃效率E與網(wǎng)格尺寸L的函數(shù)關(guān)系圖

4 結(jié)論

本文中，我們把隨機(jī)多孔介質(zhì)看做是正方網(wǎng)格，借助于決定論方法在喉管半徑為非均勻分布的假定下，模擬了粘滯指進(jìn)在逾滲集團(tuán)中的一些新的性質(zhì)。逾滲集團(tuán)中粘滯指進(jìn)掃及面積的大小隨著迭代次數(shù)n的增加而增大。當(dāng)n ≥ 10時(shí)，粘滯指進(jìn)圖像幾乎不再變化。增加逾滲概率P導(dǎo)致分形維數(shù)D的增加，與迭代次數(shù)n無(wú)關(guān)。這表明多孔介質(zhì)的幾何拓?fù)浣Y(jié)構(gòu)嚴(yán)重影響了粘滯指進(jìn)。驅(qū)掃效率E隨著迭代次數(shù)n和網(wǎng)格尺寸L的增加而減小，當(dāng)網(wǎng)格尺寸L取最大值時(shí)，E達(dá)到最小。并且在網(wǎng)格尺寸L的整個(gè)變化過(guò)程中，E～L存在兩種線性關(guān)系，分別表示粘滯指進(jìn)在逾滲集團(tuán)中的凝集區(qū)域和驅(qū)替區(qū)域。

[1] P. G. Saffman, G. I. Taylor. The Penetration of A Fluid into A Porous Medium or Hele-show Cell Containing A More Visous Fluid[J].Proc. Royal Society A, 1958, 245: 312-329.

[2] H. S. Hele-Show. The Flow of Water[J]. Nature, 1898, 58: 34-36.

[3] M. Sahimi. Flow Phenomena in Rocks[J]. Reviews of Msdern Physics, 1993,65: 1393-1534.

[4] E. Salmon, M. Ausloos, N. Vandewalle. Aging of Prous Mdia Fllowing Fuid Ivasion, Feezing and Tawing[J]. Physical Review E,1997, 55: R6348-52.

[5] K. T. Tallakstad, H. A. Knudsen, T. Ramstad. Steady-State Two-Phase Flow in Porous Media: Statistics and Transport Properties [J].Physical Review Letter, 2009, 102: 074502.

[6] Ch. Cottin, H. Bodiguel, A. Colin. Drainage in two-dimensional porous media: From capillary fingering to viscous flow [J]. Physical Review E, 2010, 82: 046315.

[7] J. P. Tian, K. L. Yao. Ageing of Random Porouse Media Following Fluid Deterministic Displacement, Freeing, Thawing[J]. The European Physical Journal B, 2000, 14:543-549.

[8] H. Siddiqui, M. Sahimi. Computer Simulations of Miscible Displacement Processes in Disordered Porous Media[J]. Chemical Engineering Science, 1990, 45: 163-182.

[9] J. P. Tian, K. L. Yao. Computer Simulations of Two-phase Flow with Surface Tension in The Percolation Cluster[J]. Physics Letters A, 1999, 251: 259-261.

[10] H. Boularot, G. Albinet. Frozen and Active Regions in Diffusion Limited Aggregation Clusters[J]. Physical Review E, 1996, 53:5106-5110.

The Immiscible Displacement in Random Porous Media

WANG Man
(College of Mathematics and Computer Science, Wuhan Textile University, Wuhan Hubei 430073, China)

In this paper, the fractal nature of Viscous fingeringin (VF) random porous media has been studied. The fluid transport is based on the deterministic method. The result shows that the topology and the geometry of the porous media have a strong effect on displacement processes. The cluster size of viscous fingering pattern in percolation cluster increases with the increase of iteration parameter n. When iteration parameter n ≥ 10, VF pattern does not change with n. The distribution of velocities normal to the interface of VF in percolation cluster is also studied. The sweep efficiency E increases along with the increasing of iteration parameter n and decreases with the network size L.

Porous Media; Percolation Cluster; Viscous Fingering; Sweep Efficiency

O189.11

A

1009－5160(2011)06－0090－05

王嫚（1986-），女，碩士研究生，研究方向：無(wú)序空間中的分形生長(zhǎng).