SCIENCE CHINA Technological Sciences • Article •
August 2015 Vol.58 No.8: 1375–1384 doi: 10.1007/s11431-015-5842-z
Pore scale simulation of liquid and gas two-phase flow based on digital core technology ZHANG Lei1,2, KANG QinJun2, YAO Jun1*, GAO Ying1, SUN ZhiXue1, LIU HaiHu3 & VALOCCHI Albert J.4 1 School of Petroleum Engineering, China University of Petroleum, Qingdao 266580, China; Earth and Environmental Sciences Division, Los Alamos National Laboratory, Los Alamos 87545, USA; 3 School of Energy and Power Engineering, Xi’an Jiaotong University, Xi’an 710049, China; 4 Department of Civil & Environmental Engineering, University of Illinois at Urbana-Champaign, Urbana IL 61801, USA 2
Received March 1, 2015; accepted May 21, 2015; published online July 7, 2015
Two-phase flow in two digital cores is simulated by the color-gradient lattice Boltzmann method. This model can be applied to two-phase flow with high-density ratio (on order of 1000). The first digital core is an artificial sandstone core, and its three-dimensional gray model is obtained by Micro-CT scanning. The gray scale images are segmented into discrete phases (solid particles and pore space) by the Otsu algorithm. The second one is a digital core of shale, which is reconstructed using Markov Chain Monte Carlo method with segmented SEM scanning image as input. The wettability of solid wall and relative permeability of a cylindrical tube are simulated to verify the model. In the simulations of liquid and gas two phase flow in digital cores, density ratios of 100, 200, 500 and 1000 between liquid and gas are chosen. Based on the gas distribution in the digital core at different times, it is found that the fingering phenomenon is more salient at high density ratio. With the density ratio increasing, the displacement efficiency decreases. Besides, due to numerous small pores in the shale, the displacement efficiency is over 20% less than that in the artificial sandstone and the difference is even about 30% when density ratio is greater than 500. As the density ratio increases, the gas saturation decreases in big pores, and even reaches zero in some small pores or big pores with small throats. Residual liquid mainly distributes in the small pores and the edge of big pores due to the wettability of liquid. Liquid recovery can be enhanced effectively by decreasing its viscosity. pore scale, digital core, liquid and gas two-phase, lattice Boltzmann method, shale Citation:
Zhang L, Kang Q J, Yao J, et al. Pore scale simulation of liquid and gas two-phase flow based on digital core technology. Sci China Tech Sci, 2015, 58: 13751384, doi: 10.1007/s11431-015-5842-z
1 Introduction Under reservoir conditions, it is hard to get the inside information of the rock by using physical experiment, and the experimental measurements could be both difficult and time consuming. An alternative or complementary technique to get the inside information of rocks is pore scale simulation
*Corresponding author (email:
[email protected]) © Science China Press and Springer-Verlag Berlin Heidelberg 2015
based on digital geometry of rocks. Recently, pore scale imaging and modelling, also called “digital core analysis” [1], has been developed rapidly in oil and gas industry. In order to study the flow mechanism in porous media, the digital pore structure should be obtained or reconstructed first. Normally, there are two methods to construct digital rocks: Physical experiments and numerical reconstruction. The physical experiments include X-ray computed tomography (CT) [2], focused ion beams-scanning electron microscopy (FIB-SEM) [3], etc. X-ray CT is tech.scichina.com link.springer.com
1376
Zhang L, et al.
Sci China Tech Sci
an advanced imaging technique that allows nondestructive and noninvasive imaging of specimens to depict crosssectional and three-dimensional internal structures, and yields high-resolution, three-dimensional representations of pore space and fluid distribution within porous materials. An overview of recent advances and developments in the application of X-ray CT to subsurface flow and transport problems can be found in Wildenschild and Sheppard [4], with a particular emphasis on multi-phase flow. Usually FIB-SEM can produce images down to resolutions of tens of nanometer, which is sufficient for low permeability sandstones and unconventional shale gas and shale oil reservoirs. However, FIB-SEM method is destructive. 3D rock images can be acquired by X-ray CT or FIB-SEM directly, but both the methods are expensive and the latter is destructive. Numerical reconstruction methods, based on analysis of highresolution two-dimensional images, can be used to construct three-dimensional representations of the pore space, including Gaussian simulation [5–7], simulated annealing [8,9], processed-based [10,11], multiple-point statistic [12], etc. In this work, the Markov Chain Monte Carlo (MCMC) method [13] will be adopted, which is fast and only requires thin section images or 2D images as input. Several pore-scale modelling methods have been applied for multiphase flow simulation in porous media, such as pore-network model [14], level set method [15], the volume-of-fluid (VOF) method [16,17] and some grid based computational fluid dynamics methods. In recent years, the lattice Boltzmann (LB) method is becoming a promising method and has been successfully applied to a variety of fields, such as flows in cavities [18,19], flows in porous media [20], chemical reactions [21,22], heat transfer [23], particle suspensions [24], dissolution and precipitation [25–27], and multiphase and multicomponent flows [28–31]. The LB method can be efficiently implemented on parallel computers and easily deal with complex boundaries. A number of LB models have been developed to simulate multiphase flow problems, like the color model proposed by Gunstensen et al. [32], the pseudo-potential model introduced by Shan and Chen [33], the free-energy model presented by Swift et al. [34], and some other models [35]. The density ratio of liquid & gas varies in different reservoirs, and is on the order of 1000 under some conditions. To the best knowledge of the authors, there is no literature on LB simulation of liquid & gas flow (density ratio on order of 1000) in digital core of real rock has been reported yet. In this work, two digital cores are constructed, including artificial core scanned by X-ray computed microtomography (Micro-CT) and shale digital core reconstructed by MCMC with SEM scanning image as input. A color model developed by Liu et al. [36] is adopted,which used a perturbation step to generate the interfacial tension and a recoloring step to promote phase segregation and maintain surfaces. Liu’s model can deal with immiscible two-phase fluids problem with high viscosity and density ratio. In the
August (2015) Vol.58 No.8
simulations, density ratios of liquid and gas are set as 100, 200, 500 and 1000, in addition, four cases for investigating the influence of capillary number are simulated. The digital cores are saturated by liquid phase, and the gas is injected at the inlet. The displacement efficiency and residual liquid distribution are studied.
2 Digital cores In this section, two digital cores are constructed. The first one is artificial sandstone core, and its three-dimensional gray model is obtained by X-ray Micro-CT scanning. The second one is digital core of shale, which is reconstructed by MCMC method with segmented SEM scanning image as input. 2.1 X-ray Micro-CT scanning digital core of artificial sandstone Three-dimensional digital core of artificial sandstone is made by stacking sequences of two-dimensional slice images. The slices are tomographic reconstructions from a large number of projection images obtained by the scanning system. The scanner used for this study is Micro-XCT 400, and 3600 projection images are obtained. Figure 1 is a slice of tomographic reconstructions. Unfortunately, X-ray CT is not free of artifacts, which complicates quantitative image analysis owing to obscuration of significant features or misinterpretation of attenuation values of a single material in different image sections [37]. Commonly experienced problems include high-frequency noise, beam hardening, defective detector pixels, scattered X rays or poorly centered samples. To partially alleviate these problems, CT images are typically preprocessed with smoothing filters, varying from simple low-pass (blurring) to advanced edge-preserving smoothing algorithms. A review of filtering methods applicable for preprocessing of X-ray CT images can be found in the work of Kaestner et al. [38]. However, image segmentation is the most crucial step affecting all subsequent quantitative analysis and modeling efforts [39]. Many binarization approaches produce poor results when applied directly because CT images commonly suffer from high noise levels. In this study, the non-local means algorithm [40] is adopted to alleviate noise, and then, the digital core is segmented by Otsu [41] algorithm. A core with size of 160×160×160 is cropped for flow simulation, as shown in Figure 2, and its resolution is 4.08 microns. The porosity of the digital core is 61.8%. 2.2
Reconstructing digital core of shale
The shale rock is located in Sichuan Basin of China. Figure 3 is the SEM scanning image of the shale core, with resolution of 2.8 nm per pixel. Firstly, the SEM image is segmented
Zhang L, et al.
Figure 1
Figure 2
Sci China Tech Sci
Gray scale CT image of artificial sandstone.
August (2015) Vol.58 No.8
1377
are labeled as orange, with a total of 0.6% of the whole volume or 3.5% of the pore volume, and connected pores are labeled as blue, with porosity decreasing to 16.2%. The pore size distribution of the shale sample is calculated by using the thirteen directions approach in ref. [42], which is based three dimensions method by extending the two dimensions method in modelling Knudsen flow through porous media. For each pore cell, an average pore diameter is computed. A pore length is computed by moving in a given direction (both backward and forward) until a nonpore cell is reached and a “length” of pore cells in the given direction can be computed. The x, y, and z directions are used. Two diagonal directions for each of the x−y, x−z, and y−z planes are used along with four diagonal directions which traverse through 3D space. These 13 different lengths are averaged to compute an “effective pore diameter” for a given pore cell. The pore size distribution is shown in Figure 6. It can be seen that over 92% of the pore diameters are
Three-dimensional digital core after segmentation.
into binary image, only contains black and white pixels. Then, the three-dimensional digital core model is reconstructed by MCMC method, with the binary image as input. In the MCMC method, Markov Chain is applied into Monte Carlo method to realize the dynamic simulation. The essence of MCMC method is to construct a Markov Chain with stable distributions. Based on the 2D binary images, the transition probabilities between voids and solids controlling the Markov Chain are determined. The MCMC method takes spatial structural information into account that identifies all the transition probabilities for a given local training lattice stencil. A multiple-voxel interaction scheme is developed in order to generate individual realizations. Such reconstruction approach and the structures it generates are referred to as “pore architecture models (PAMs)”. Generally speaking, four procedures are used to reconstruct the 3D digital rock, the details of the PAMs and reconstruction procedures are described in Wu et al. [13]. The generated three-dimensional digital core model is shown in Figure 4, with porosity of 16.8% and size of 300× 100×100. The blue part is pore structure, and the rest part is solid matrix. As the unconnected pores have no contribution to the flow, they are marked as solid matrix to save computational resources. As shown in Figure 5, unconnected pores
Figure 3
SEM scanning image of shale in Sichuan Basin.
Figure 4
Three-dimensional digital core of shale (the blue part is pores).
Figure 5 (blue).
Illustration of disconnected pores (orange) and connected pores
1378
Zhang L, et al.
Sci China Tech Sci
August (2015) Vol.58 No.8
9 3 2 3 f i k ,eq k i k i 2 ei u 4 ei u 2 u2 , (4) 2 2 c c c where i k is related to the density ratio , i is the weight coefficients and u is the local fluid velocity. And i k is given as k , i 1 k / 12, 1 k / 24,
i 0, i 1, 2, , 6,
k
Figure 6
in the range of 8.4 to 25.6 nm and 66% are between 11.2 and 19.6 nm (4 to 7 lattices).
k
2
i
3 Lattice Boltzmann method The Liu’s LB model is introduced briefly in this section. For more details one can refer to ref. [36]. In Liu’s model, the density ratio can be set as high as 1000 and the viscosity ratio can be set to 100, which can cover most situations of liquid and gas. There are two phases (red and blue) defined in color model [27]. Let f i k x, t represents the k phase particle distribution function (PDF) in i-th velocity direction at position x and time t, where k = R or B denotes the red or blue color phase. The total PDF is defined as f i f i R f i B . The evolution equation for each phase is (1)
time step, and i is the collision operator. The collision operator contains three parts:
(i )
i i
k
i
1
3
k i
1
i
k
2
.
(2)
k
i
1
1
k
f
k i
fi k ,eq
,
where k is the relaxation time of k phase and f i k
k , eq
the equilibrium distribution function (EDF) of fi . The EDF is chosen to consider different densities for red and blue phases, and is defined as follows:
2
Bi ,
(7)
R x, t B x, t , R x, t B x, t
(8)
We assume that AR AB A , and the interfacial tension satisfies 4 / 9A . To account for different viscosities of red and blue phase, the viscosity of fluid mixture is determined by 1 1 N 1 N , 0.5 2 R 0.5 2 B 0.5
(9)
which can ensure the continuity of viscosity flux across the interface [43]. The relation between viscosity and dimensionless relaxation time is k k 0.5 t c 2 / 3 . k
is
and Bi is set to B0 1 / 3 , B1 6 1 / 18 , B7 18 1 / 36 .
i
(3)
e N A i k N N 2 2
N x, t
is the BGK collision operator, and defined as
(6)
where the phase field N is defined as
k
k
R 1 B . B 1 R
is the perturbation operator, and is defined as
k (2)
where ei is the lattice velocity in i-th direction, t is the
k
i 7,8, ,18,
k is a free parameter and should satisfy 0 k 1 . In the two phase interface, there is a relation between R and B :
Pore size distribution of reconstructed digital core.
f i k x ei t , t t f i k x, t i k x, t ,
(5)
3
is the recoloring operator, and is defined by
R * fi (i k )(3) fi B B f i* (i k )(3) fi R
R B cos i f i eq 2 R B cos i f i eq 2
u= 0
,
(10) u= 0 ,
where fi * denotes the postperturbation, presegregation value of the total PDF in i-th direction, i is the angle between
Zhang L, et al.
Sci China Tech Sci
the phase field gradient N and the lattice vector ei , is the segregation parameter related to the interface thickness, and its value must be between 0 and 1 to ensure positive particle distribution functions. And cos i is defined by cos i
ei N ei N
.
(11)
The different densities and viscosities are considered in EDF and total viscosity, respectively.
4 Model validation This section uses two benchmark cases to validate the color-fluid model. Subsection 4.1 demonstrates that different contact angles can be obtained by adjusting the value of the phase field at the solid wall. Subsection 4.2 compares the relative permeability of co-current annular flow in a cylindrical tube with analytic predictions. All the computations in this study are implemented in Los Alamos National Laboratory’s high performance computer Mustang. Mustang is an Appro Xtreme-X cluster with 1600 compute nodes. Each node has two 12-Core AMD Opteron model 6176 processors and a total of 24 processor cores. The CPUs run at 2.3 GHz. Each node has 64 GB DDR3 memory (1333 MHz), or 2.6 GB memory per core. The L3 cache is 12 MB, shared among the 12 cores of a socket. Therefore, it has 102.4 TB memory for 38400 cores. Its aggregate peak cluster performance is (2.3 GHz×4 ops/clock×24 cores/node×1600 nodes) 353 TF/s. 4.1
Contact angle
In the present model, we follow the assumption of Rowlinson and Widom [44] that the solid wall is a mixture of two fluids, thus having a certain value of the phase field N. In the LB community, similar approaches were also adopted in the free-energy models and mean-field theory model as well. In the simulations, the initial condition is a hemisphere stationary droplet (blue fluid) with radius R = 20 lattices sitting along the center of the bottom wall. The size of computational area is 100 × 100 × 100. The top wall is assumed to be neutral wetting. The top and bottom boundaries are solid wall and standard bounce back condition is adopted. The periodic boundary condition is implemented on the other boundaries. The densities of both phases are fixed at 1.0, the kinetic viscosities are 0.166, and the surface tension is 0.03. The segregation parameter is set to 0.7. The contact angle () of blue phase is set to 30°, 60°, 90°, 120° and 150° for each case. At the same time, the densities of red and blue phase are set as 1 cos and 1 cos on the solid wall, and the value of the phase field is
August (2015) Vol.58 No.8
1379
cos( ) . We run the simulations until the droplet shape does not change, i.e., reaching an equilibrium state. The final steady states are shown in the Figure 7. Figure 7 demonstrates that different contact angles can be obtained by adjusting the value of the phase field at the solid wall.
4.2 Relative permeability
Immiscible two-phase flow in porous media commonly described in terms of a generalized Darcy’s law using the concept of relative permeability. Here, the co-current annular flow in a cylindrical tube is considered. The wetting fluid flows along the tube walls with non-wetting fluid in the center. The relative permeability is given as follows, at viscosity ratio equals 1 [45]: kr , w S w 2 , 2 kr , nw 2 1 S w 1 S w ,
(12)
where subscripts w and nw represent wetting phase and non-wetting phase separately. The length of the tube is 100, and the radius R is 30 (both in lattice unit). The periodic boundary condition is implemented on the inlet and outlet. The standard bounce back condition is adopted at the solid wall. The densities of both phases are fixed at 1.0, the kinetic viscosities are 0.03, and the surface tension is 0.03. The segregation parameter is set to 0.7. Driving of the fluids is done by applying a uniform body force to the PDFs. The details about the forcing term are given by Guo et al. [46]. To measure the relative permeability, we first calculate the total momentum of each phase at a given wetting phase saturation S w
S w u .
(13)
The relative permeability in each phase is given as kr ,
S w P S 1 . S 1 P S w
(14)
Figure 8 compares the simulated relative permeability
Figure 7
Different contact angles (30°, 60°, 90°, 120°, 150°) at the wall.
1380
Zhang L, et al.
Sci China Tech Sci
August (2015) Vol.58 No.8
as Ca
Figure 8
Relative permeability for co-current flow in a cylindrical tube.
with the analytic results from eq. (12). The simulated results are in excellent agreement with the analytic results. The relative permeability of both phases is less than one and we notice that kr , w kr , nw 1 .
5 Liquid & gas two-phase flow simulation In this section, liquid & gas two-phase flow simulations are conducted in digital core of artificial sandstone and reconstructed digital core of shale rock, which are introduced in section 2. The gas phase distribution in digital core at final equilibrium state will be clearly illustrated at different conditions. 5.1
un n , where un and n are the velocity cos
and dynamic viscosity of the advancing non-wetting fluid, respectively. The capillary number is fixed at 4×105. The surface tension is set to 0.03. The contact angle of liquid phase is 30 . The dynamic viscosity of liquid varies considerably under different pressures or temperatures. Here the kinetic viscosities of liquid and gas phases are fixed at 0.03, the density of gas is fixed at 1 for simplicity. Four cases with liquid densities of 100, 200, 500 and 1000 are simulated. In these four cases, the drainage velocity of gas are kept with same value, due to the same capillary number, gas density and gas kinetic viscosity. Figure 9 shows the gas distribution in the digital core at t = 5×105. The gas saturations are 12.4%, 12.2%, 12.1% and 11.3% when liquid densities are 100, 200, 500 and 1000. The gas saturations are almost the same for all cases, due to the same drainage velocity at the inlet. As shown in Figure 9, the higher the liquid density, the deeper the gas invades into the domain. The fingering phenomenon is more salient at higher density ratio. The dynamic viscosity of liquid phase increases with the density increasing, and the velocity will decrease with the density increasing under the same pressure gradient. Figure 10 shows the gas phase distribution in digital core at final equilibrium state, where the gas saturations are 59.2%, 54.7%, 43.2% and 33.3% when liquid densities are 100, 200, 500 and 1000. As the density becomes higher, the dynamic viscosity of liquid is bigger and the liquid is more difficult to be displaced, resulting in lower displacement efficiency, and higher residual liquid saturation. Hence, liquid recovery can be enhanced effectively by decreasing its viscosity.
Two-phase flow in artificial sandstone
The simulation is conducted in the three-dimensional cropped digital core of artificial sandstone, shown in section 2.1. Firstly, we need check the connectivity of the pores and mark the unconnected pores as solid matrix. This model only contains one connected pore. The slice of x = 0 is the inlet and the slice of x = 160 is the outlet. Twenty slices of pore are added to the inlet and outlet, in order to eliminate the effects of the solid phase. The final size of the simulated domain is 200 × 160 × 160 lattices. Initially, the digital core is fully saturated with the wetting phase (liquid). The non-wetting phase (gas) is injected continuously at the left inlet, with a constant flow rate, while a constant pressure is imposed at the right outlet. Therefore, the liquid and gas phase are the displaced and advancing phase, respectively. This is a typical drainage process. The boundaries on y- and z-directions are assumed to be solid walls. The standard bounce back condition is applied on all the solid boundaries. The injecting velocity of the gas is determined by capillary number, which is defined
Figure 9 The gas phase distribution in digital core at t = 5×105 (The densities of liquid phase are, in order, 100, 200, 500 and 1000.).
Zhang L, et al.
Sci China Tech Sci
August (2015) Vol.58 No.8
1381
sional digital core model is 300×100×100. There are 375 cores adopted in parallel computing for the current core. The size of each small block is 20×20×20.
Figure 10 The gas phase distribution in digital core at equilibrium state (The densities of liquid phase are, in order, 100, 200, 500 and 1000).
5.2
Two-phase flow in shale digital core
In this section, four cases with different liquid densities are simulated on the shale digital core that reconstructed in section 2.2. The slice of x = 0 is the inlet and the slice of x = 300 is the outlet. In order to eliminate the effect of the solid phase, the first and last twenty slices near the inlet and outlet are set as pore space. The boundary conditions are set as in section 5.1. Velocity boundary is set at the inlet, and pressure boundary is set at the outlet and standard bounce back condition is set on all the solid walls. Initially, the digital core is fully saturated with the liquid phase. The gas phase is injected continuously at the inlet with a constant flow rate. As mentioned before, the size of three-dimen-
Figure 11
5.2.1 The influence of density ratio The capillary number is set to 4×104, and other parameters are kept the same as in section 5.1 for comparison. The capillary number for these simulations is relatively high, Ca = 4×104. The reason for choosing a higher Ca compared with that in section 5.1, is to obtain a complete invasion within a relatively short time. The densities for liquid are 100, 200, 500 and 1000 as well. In these four cases, the drainage velocity of gas is kept with the same value, due to the same capillary number, gas density and gas kinetic viscosity. Figure 11 shows the gas phase distribution in shale digital core at the final equilibrium state, and the gas saturations are 43.9%, 34.2%, 13.3% and 6.9% when liquid densities are 100, 200, 500 and 1000. The same tendency as section 5.1 is obtained, the displacement efficiency decreases as the liquid density increasing. However, the displacement efficiency is over 20% less than that in artificial sandstone. The main reason is that numerous small size pores are distributed in the shale, as shown in Figure 6. The two phases distribution slices at y = 50 based on the phase field N are shown in Figure 12. Three parts with
N 0 , 1 N 0 and N 1 are gas phase, liquid phase and solid matrix, respectively. As the liquid density increases, the gas saturation in big pores decreases, and is even zero in some small size pores as well as big pores with small size throat. We can see that the residual liquid mainly distribute in the small size pores and the edge of big pores due to the wettability of liquid. Comparison of the gas saturation between artificial sandstone and shale is shown in Figure 13. The x coordinate is density ratio ( g /o ). The difference is about 30% when density ratio is less than 0.02 (the liquid density is bigger
The gas phase distribution in shale at equilibrium state (The densities of liquid phase are, in order, 100, 200, 500 and 1000).
1382
Zhang L, et al.
Sci China Tech Sci
August (2015) Vol.58 No.8
Figure 12 The slice of two phases distribution at y = 50 based on the phase field N (The densities of liquid phase are, in order, 100, 200, 500 and 1000. The blue part is solid matrix where N<1).
tively. Figure 14 shows the gas phase saturation as a function of capillary number at breakthrough, and the gas saturations are 60.7%, 42.5%, 35.3% and 37.2% as capillary numbers decreases. The displacement efficiency increases as the capillary number increasing. According to the definition of capillary number, the velocity of advancing phase fluid increases with the capillary number increasing and the pressure gradient increases. Therefore, fluid in smaller size pores can be displaced, which results in high displacement efficiency. In addition, the rate of increase is nearly constant, and the same near linear increase tendency was reported by Cottin et al. [47] and Zhang et al. [48].
Figure 13
The gas phase saturation versus density ratio.
than 500). The effect of pore size is more significant at low density ratio. For the two-phase flow in tight rock (such as shale), the phase with larger density is more likely to be trapped in the rock, when the density ratio is lower. For instance, when the liquid density is 1000, the gas saturation is only 6.9% in shale. 5.2.2 The influence of capillary number The density ratio has significant influence in the liquid & gas flow; in addition, the influence of capillary number can not be neglected. In this section, the displacement efficiency at breakthrough is studied with different capillary numbers. The kinetic viscosities of gas and liquid phases are fixed at 0.03 and 0.3 respectively. The density ratio is set to 1.0. The contact angle and the surface tension are set as previous section. Four capillary numbers are chosen in the simulations, which are 4×103, 8×104, 4×104, 8×105, respec-
6
Conclusion
In this paper, two kinds of digital core, artificial sandstone
Figure 14
The gas phase saturation versus capillary number.
Zhang L, et al.
Sci China Tech Sci
core and shale digital core, are constructed using different methods. The three-dimensional gray scale model of artificial sandstone is obtained by Micro-CT scanning. The gray scale images are segmented into discrete phases (solid particles and pore space) by Otsu algorithm. The shale digital core is reconstructed by Markov Chain Monte Carlo method with segmented SEM scanning image as input. Both digital cores are simulated by lattice Boltzmann color model, which can be applied to two-phase flow with high-density ratio. It is proved that the model can accurately reproduce wetting effects on solid wall. Besides, the relative permeability of co-current annular flow in a cylindrical tube is compared with analytic predictions. The simulated results show that they are in excellent agreement with the analytic results. The liquid & gas two-phase flow simulation is implemented in both constructed digital cores. During the simulation, four different liquid densities (100, 200, 500 and 1000) are chosen, with the density of gas fixed at 1 for simplicity. Based on the gas distribution in the digital core at different times, the fingering phenomenon is more salient at higher density ratio. The displacement efficiency decreases with the liquid density increases, which results in higher residual liquid saturation. Besides, due to numerous small size pores distributed in shale, the displacement efficiency is over 20% less than that in artificial sandstone and the difference is even about 30% when density ratio is less than 0.02 (the liquid density is bigger than 500). As the liquid density increases, the gas saturation in big pores decreases, and the gas saturation is even zero in some small size pores and big pores with small size throat. We can see that the residual liquid mainly distribute in the small size pores and the edge of big pores due to the wettability of liquid. Liquid recovery can be enhanced effectively by decreasing its viscosity. Besides, four capillary numbers are chosen for studying its influence. The displacement efficiency increases as the capillary number increasing, due to the velocity of advancing phase fluid increases and the pressure gradient increases with the capillary number increasing. Fluid in smaller size pores can be displaced, which results in high displacement efficiency. However, not only the density ratio and capillary number, but also different rock microstructures have effects on displacement efficiency. We will investigate it systematically in our future work. This work was supported by the National Natural Science Foundation of China (Grant No. 51234007, 51404291), Program for Changjiang Scholars and Innovative Research Team in University (Grant No. IRT1294), and Introducing Talents of Discipline to Universities (Grant No. B08028). The author ZHANG Lei would like to give appreciation to Chinese Scholarship Council for supporting the one-year study in Los Alamos National Laboratory. Valocchi A J and Liu H gratefully acknowledge the support of the International Institute for Carbon Neutral Energy Research (WPI-I2CNER), sponsored by the Japanese Ministry of Education, Culture, Sports, Science and Technology. Kang Q acknowledges the support from
August (2015) Vol.58 No.8
1383
the LDRD Program and Institutional Computing Program of Los Alamos National Laboratory. 1 2
3
4
5 6 7 8 9 10 11
12
13
14
15
16 17 18
19
20
21
22
23
Blunt M J, Bijeljic B, Dong H, et al. Pore-scale imaging and modelling. Adv Water Resour, 2013, 51: 197–216 Dunsmuir J H, Ferguson S, D'amico K, et al. X-ray microtomography: A new tool for the characterization of porous media. SPE-22860-MS, 1991 Tomutsa L, Silin D, Radmilovic V. Analysis of chalk petrophysical properties by means of submicron-scale pore imaging and modeling. SPE-99558-PA, 2007 Wildenschild D, Sheppard A P. X-ray imaging and analysis techniques for quantifying pore-scale structure and processes in subsurface porous medium systems. Adv Water Resour, 2013, 51: 217–246 Joshi M Y. A class of stochastic models for porous media. Dissertation for the Doctor Degree. Lawrence: University of Kansas, 1974 Quiblier J A. A new three-dimensional modeling technique for studying porous media. J Colloid Interf Sci, 1984, 98: 84–102 Adler P, Jacquin C, Quiblier J. Flow in simulated porous media. Int J Multiphase Flow, 1990, 16: 691–712 Hazlett R. Statistical characterization and stochastic modeling of pore networks in relation to fluid flow. Math Geol, 1997, 29: 801–822 Yeong C, Torquato S. Reconstructing random media. Phys Rev E, 1998, 57: 495 Bryant S, Blunt M. Prediction of relative permeability in simple porous media. Phys Rev A, 1992, 46: 2004 Øren P-E, Bakke S. Process based reconstruction of sandstones and prediction of transport properties. Transport Porous Med, 2002, 46: 311–343 Okabe H, Blunt M J. Prediction of permeability for porous media reconstructed using multiple-point statistics. Phys Rev E, 2004, 70: 066135 Wu K, Van Dijke M I, Couples G D, et al. 3D stochastic modelling of heterogeneous porous media–applications to reservoir rocks. Transport Porous Med, 2006, 65: 443–467 Blunt M, King P. Relative permeabilities from two-and threedimensional pore-scale network modelling. Transport Porous Med, 1991, 6: 407–433 Osher S, Sethian J A. Fronts propagating with curvature-dependent speed: algorithms based on Hamilton-Jacobi formulations. J Comput Phys, 1988, 79: 12–49 Hirt C, Nichols B. Volume of fluid (VOF) method for the dynamics of free boundaries. J Comput Phys, 1981, 39, 201–225 Jiang C B, Deng B, Hu S X, et al. Numerical investigation of swash zone hydrodynamics. Sci China Tech Sci, 2013, 56: 3093–3103 Tang X L, Su Y W, Wang F J, et al. Numerical research on lid-driven cavity flows using a three-dimensional lattice Boltzmann model on non-uniform meshes. Sci China Tech Sci, 2013, 56: 2178–2187 Huang B, Wu Q, Wang G Y. Numerical simulation of unsteady cavitating flows around a transient pitching hydrofoil. Sci China Tech Sci, 2014, 57: 101–116 Martys N S, Chen H. Simulation of multicomponent fluids in complex three-dimensional geometries by the lattice Boltzmann method. Phys Rev E, 1996, 53: 743 Chen L, Kang Q, He Y L, et al. Pore-scale simulation of coupled multiple physicochemical thermal processes in micro reactor for hydrogen production using lattice Boltzmann method. Int J Hydrogen Energy, 2012, 37: 13943–13957 Chen L, He Y L, Kang Q, et al. Coupled numerical approach combining finite volume and lattice Boltzmann methods for multi-scale multi-physicochemical processes. J Comput Phys, 2013, 255: 83–105 Biferale L, Perlekar P, Sbragaglia M, et al. Convection in multiphase
1384
24 25 26 27
28
29
30
31 32 33 34 35
36
Zhang L, et al.
Sci China Tech Sci
fluid flows using lattice Boltzmann methods. Phys Rev Lett, 2012, 108: 104502 Joshi A S, Sun Y. Multiphase lattice Boltzmann method for particle suspensions. Phys Rev E, 2009, 79: 066703 Kang Q, Zhang D, Chen S, et al. Lattice Boltzmann simulation of chemical dissolution in porous media. Phys Rev E, 2002, 65: 036318 Kang Q, Zhang D, Chen S. Simulation of dissolution and precipitation in porous media. J Geophys Res, 2003, 108: 2505 Chen L, Kang Q, Robinson B A, et al. Pore-scale modeling of multiphase reactive transport with phase transitions and dissolutionprecipitation processes in closed systems. Phys Rev E, 2013, 87: 043306 Chen L, Luan H B, He Y L, et al. Pore-scale flow and mass transport in gas diffusion layer of proton exchange membrane fuel cell with interdigitated flow fields. Int J Therm Sci, 2012, 51: 132–144 Kang Q, Lichtner P C, Zhang D. Lattice Boltzmann pore-scale model for multicomponent reactive transport in porous media. J Geophys Res, 2006, 111: doi: 10.1029/2005JB003951 Chen L, Kang Q, Mu Y, et al. A critical review of the pseudopotential multiphase lattice Boltzmann model: Methods and applications. Int J Heat Mass Tran, 2014, 76: 210–236 Gui N, Xu W K, Ge L, et al. LBE-DEM coupled simulation of gassolid two-phase cross jets. Sci China Tech Sci, 2013, 56: 1377–1386 Gunstensen A K, Rothman D H, Zaleski S, et al. Lattice Boltzmann model of immiscible fluids. Phys Rev A, 1991, 43: 4320–4327 Shan X, Chen H. Lattice Boltzmann model for simulating flows with multiple phases and components. Phys Rev E, 1993, 47: 1815–1819 Swift M R, Osborn W R, Yeomans J M. Lattice Boltzmann Simulation of Nonideal Fluids. Phys Rev Lett, 1995, 75: 830–833 He X, Chen S, Zhang R. A Lattice Boltzmann Scheme for Incompressible Multiphase Flow and Its Application in Simulation of Rayleigh–Taylor Instability. J Comput Phys, 1999, 152: 642–663 Liu H, Valocchi A J, Kang Q. Three-dimensional lattice Boltzmann model for immiscible two-phase flow simulations. Phys Rev E, 2012,
August (2015) Vol.58 No.8
37
38
39
40 41 42
43
44 45 46 47
48
85: 046309 Ketcham R A, Carlson W D. Acquisition, optimization and interpretation of X-ray computed tomographic imagery: applications to the geosciences. Comput Geosci, 2001, 27: 381–400 Kaestner A, Lehmann E, Stampanoni M. Imaging and image processing in porous media research. Adv Water Resour, 2008, 31: 1174–1187 Iassonov P, Gebrenegus T, Tuller M. Segmentation of X-ray computed tomography images of porous materials: A crucial step for characterization and quantitative analysis of pore structures. Water Resour Res, 2009, 45: W09415 Buades A, Coll B, Morel J M. A non-local algorithm for image denoising. IEEE Computer Society Conf, 2005.2: 60–65 Otsu N. A threshold selection method from gray-level histograms. Automatica, 1975, 11: 23–27 Lange K J, Sui P C, Djilali N. Pore scale simulation of transport and electrochemical reactions in reconstructed PEMFC catalyst layers. J Electrochem Soc, 2010, 157: B1434–B1442 Zu Y, He S. Phase-field-based lattice Boltzmann model for incompressible binary fluid systems with density and viscosity contrasts. Phys Rev E, 2013, 87: 043301 Rowlinson J S, Widom B. Molecular Theory of Capillarity. Oxford: Clarendon Press, 1982 Goldsmith H L, Mason S G. The flow of suspensions through tubes. II. Single large bubbles. J Colloid Sci, 1963, 18: 237–261 Guo Z, Zheng C, Shi B. Discrete lattice effects on the forcing term in the lattice Boltzmann method. Phys Rev E, 2002, 65: 046308 Cottin C, Bodiguel H, Colin A. Drainage in two-dimensional porous media: From capillary fingering to viscous flow. Phys Rev E, 2010, 82: 046315 Zhang C, Oostrom M, Wietsma T W, et al. Influence of viscous and capillary forces on immiscible fluid displacement: pore-scale experimental study in a water-wet micromodel demonstrating viscous and capillary fingering. Energ Fuel, 2011, 25: 3493–3505