Single- and two-phase flow simulation based on equivalent pore network extracted from micro-CT images of sandstone core.

Due to the intricate structure of porous rocks, relationships between porosity or saturation and petrophysical transport properties classically used for reservoir evaluation and recovery strategies are either very complex or nonexistent. Thus, the pore network model extracted from the natural porous media is emphasized as a breakthrough to predict the fluid transport properties in the complex micro pore structure. This paper presents a modified method of extracting the equivalent pore network model from the three-dimensional micro computed tomography images based on the maximum ball algorithm. The partition of pore and throat are improved to avoid tremendous memory usage when extracting the equivalent pore network model. The porosity calculated by the extracted pore network model agrees well with the original sandstone sample. Instead of the Poiseuille's law used in the original work, the Lattice-Boltzmann method is employed to simulate the single- and two- phase flow in the extracted pore network. Good agreements are acquired on relative permeability saturation curves of the simulation against the experiment results.

Background

Accurate acquisition of the micro structure and flow properties of porous media is of great importance in petroleum engineering, biomedical science, micro-electronics, and composites (Hilfer and Zauner 2011; Sanya et al. 2014; Tagar et al. 2016). Thus, the pore-scale modeling, where fluid displacement is simulated in a model of porous medium, has been used as a platform to study multi-phase flow at the pore scale.

The pore network was described by Fatt (1956a, b, c) for the first time. He randomly assigned radii to the two-dimensional (2D) regular lattice and predicted the capillary pressure and relative permeability of drainage. Since then, extensive studies on pore network models have been carried out to develop the predictive models of multi-phase flow focusing on the topology of the structure, pore size distribution, throats and their spatial correlation (Chatzis and Dullien 1977; Jerauld and Salter 1990; Lowry and Miller 1995). However, these realistic pore networks are consequently unable to reflect the real topology and complex geometry of natural porous media. In most of these studies, flow governing equations are obtained from the Poiseuille’s law. To accurately describe the flow conductance, the Poiseuille radius and length should correctly reflect the real pore and throat shape (Oren et al. 1998; Sholokhova et al. 2009). With the development of imaging technology, the three-dimensional (3D) images with resolution at the micron scale can be acquired by various approaches, such as micro computed tomography scanning technique (micro-CT), scanning electron microscopy (SEM), serial sectioning, confocal laser scanning microscopy, and reconstructed porous media by mathematical methods. Benefited from this, the reconstruction of pore network representing real rock structures has been extended from 2D thin sections (Laroche and Vizika 2005; Liu and Song 2015) or numerical reconstructions (Blunt 1998, 2001) to three-dimensional (3D) images at the micron scale (Bauer et al. 2012; Liu et al. 2015) currently.

Based on the maximum ball algorithm, a modified method of extracting the equivalent pore network model from the three-dimensional micro computed tomography images is presented in this paper. The equivalent pore network model is able to reproduce the intricate micro structure of the natural porous media sample properly. Based on the equivalent pore network model, simulation codes on single- and two-phase flow using the Lattice-Boltzmann method are developed to predict the flow properties. Meanwhile, the effects of wettability on oil recovery are studied.

Pore network model extraction algorithm

A improved method of extracting the equivalent pore network model from the three-dimensional micro computed tomography (micro-CT) images based on the maximum ball algorithm (Silin and Patzek 2006) is proposed in this paper. In the original work, the maximal ball algorithm starts from each voxel in the pore space to find the largest inscribed spheres that just touch the grain or the boundary, and the pore and throat are segmented by ranking the local radii of the spheres (Silin and Patzek 2006), which consume tremendous memory to find the maximum spheres and remove the smaller balls included by larger ones. Here, some improvements are made to simplify the process of building the maximum spheres and the partition of pores and throats. The workflow of this algorithm can be summarized in the following steps.

Image acquire and segmentation. In this paper, two sandstone samples named by ST1 and ST2 from Shengli Oilfield in China are imaged in the State Key Laboratory of Oil and Gas Reservoir Geology and Exploitation in Southwest Petroleum University and used an input to the equivalent pore network model reconstruction. Image segmentation is a crucial step in image analysis to turn the grey scales of the original micro-CT image into black and white representing the void and solid. Meanwhile, the median filter algorithm is used to remove noise, especially the salt and pepper noise which causes tiny and isolate pores in the model. The 3D segmented and filtered images of ST1 and ST2 containing 4003 pixels are shown in Fig. 1, in which the black part represents the pore and the cyan part is the matrix.

Fig. 1

3D segmented image of ST1 and ST2. The black part is the pore and the cyan part is the matrix

Medial axis extraction and the pore and throat partition. The medial axis can be obtained by a pore space burning algorithm (Lindquist et al. 1996; Lindquist and Venkatarangan 1999), and the burn number is recorded for each medial axis, which defines the radius of maximum spheres. Then the maximum spheres are centered along the medial axis. The intersection of greater than or equal to three median axes and the balls on the boundary are defined as pores, and the others are throats. The throats are replaced by cylinders with the average radius of the biggest balls in the pore-throat chains linking two pores. The minimum size of a pore is defined in codes. For example, a pore space containing two or less voxels is defined as a throat, but cannot be pore in numerical simulation in this paper. The median axis and initial partition of ST1 is illustrated in Fig. 2, in which the line represents the throat and the small ball represents the pore. The extracted equivalent pore network model of ST1 and ST2 are illustrated in Fig. 3.

Fig. 2

Median axis and initial partition of pores and throats. The lines represents the throat and the small balls represent pores

Fig. 3

EPNMs of ST1 and ST2. Different color in the image represents different radii

LBGK mathematical model

In recent years, the Lattice-Boltzmann method (LBM) has been developed into an alternative and promising numerical scheme for the simulation of fluid flows and physics modeling in fluids. The scheme is particularly effective in fluid flow applications involving interfacial dynamics and complex boundaries (porous media, e.g.) (Chen and Doolen 1998). Several algorithms have been developed since the Lattice-Boltzmann method was proposed by Lugwig Boltzmann in 1872, among which, the Bhatnagar–Gross–Krook (BGK) collision term (Qian 1993; Chen et al. 1995) has been treated as a popular algorithm in field of shock formation, multi-phase flow, and porous material property test (Crouse et al. 2002; Luan et al. 2011; Gooneratne et al. 2013). The Lattice Bhatnagar–Gross–Krook (LBGK) model used in this paper can be expressed as:where is the particle distribution function at location x and time t along the ith direction (i = 0, 1, 2…18) using 3D nineteen velocity model; is the collision term; is the velocity vector.

The two phase flow mathematical model is the color-gradient model proposed by Gunstensen et al. (1991). For a two-phase flow, the local fluid particle distribution f is defined as two parts (Gunstensen et al. 1991),where and refer to the red and blue fluid in the two-phase system; ρ and ρ are the fluid density of red and blue color; ν is the angle of local color gradient versus the lattice calculated direction; is a quadratic expansion of the Maxwell–Boltzmann distribution; β represents the separation trends of the two immiscible fluids between 0 and 1.

By introducing the surface tension term, the collision term can be calculated by (Gunstensen et al. 1991),where τ is the single time relaxation parameter; σ is the surface tension; ν is the angle of the ith direction of the lattice.

The conserved quantities of mass and momentum are calculated by,where ρ and u is the density and the local velocity, respectively.

The interaction between different fluid phases is defined as (Ramstad et al. 2010),where G is a parameter that controls the strength of the inter-particle force.

When we assume the density of rock is ρ, the wettability of the rock is defined by the interaction between the solid and fluids, and can be described as (Ramstad et al. 2010),where ω is the weight coefficient; is an indicator function that is equal to 1 or 0 for a solid or a fluid domain node, respectively; ρ is used to represents the different wall properties for different contact angles.

Then the macroscopic flow rate Q can be calculated to determine the intrinsic permeability. The intrinsic permeability K of the network is obtained from Darcy’s law at complete saturation. The absolute permeability K of the network is derived from Darcy’s law,where the network is fully saturated with a single phase j of viscosity μ; Q is the total single phase flow rate through the pore network of length L with the potential drop ΔP; A is the cross-sectional area of the model.

Then relative permeability iswhere Q is the total flow rate of phase j in multiphase conditions with the same imposed pressure drop ΔP.

Results

Some basic parameters of the pore network models are listed in Table 1. It is shown that the porosity of the equivalent pore network models is close to experimental data of the original sandstone samples in the case of that the spheres volume is equal to the original pore space. Based on the extracted equivalent pore network models and LBM, the simulation codes for single- and two- phase flow are developed. As is listed in Table 2, the absolute permeability of the equivalent pore network models are computed and compared to the experimental data, in which good agreement is acquired between the absolute permeability of the equivalent pore network models and the experimental benchmark data.

Table 1

Basic parameters of the models

	Size (mm)	Number of pores	Number of throats	Average connection number	Porosity of EPNM (%)	Experimental porosity (%)
ST1	2.14	6298	12,558	3.92	19.27	18.73
ST2	1.46	810	1563	3.72	13.71	13.11

Table 2

Absolute permeability of EPNMs

	EPNM/mD	MB/mD
ST1	1210	1100
ST2	872	703

In the water flooding process, the inlet or outlet boundaries have been selected the same as the actual experimental water flooding process, which means, the top and bottom of the model is regarded as the inlet and outlet respectively, while the other boundaries are treated as impermeable. The model is initially saturated by oil and water is injected from the inlet. Fluid parameters used in the simulation are listed in Table 3. Water volume fraction distribution for different saturation in the water flooding process of the extracted pore network is shown in Fig. 4, from which it is found that the oil is displaced in the pore space along with the water invading from the inlet gradually. Oil and water move along the respective channel to reduce the flow resistance, and it is rare for oil and water to flow alternate in the same pore. The displacement lag is remarkable. This phenomenon is called viscous fingering caused by fluid flowing along the lower flow resistance channels. Meanwhile, oil in tiny throats is impossible to be displaced by water because of the tremendous capillary force. To show the displacement process in detail, the water volume fraction distribution is illustrated in this paper.

Table 3

Fluid parameters used in the simulation

Interfacial tension σ (mN/m)	Water viscosity μ w (cp)	Oil viscosity μ o (cp)	Water density ρ w (kg/m3)	Oil density ρ o (kg/m3)	p inlet (Pa)	p outlet (Pa)
48.0	0.991	5.23	1000	890	10	0

Fig. 4

Water volume fraction distribution of ST1 for different water saturation. Each point represents the center of pore

In the simulation, the intrinsic contact angle is assigned as a random distribution to pores and throats in the network in a distribution interval. By adopting the average distribution of intrinsic contact angle intervals, different wetting systems from water-wet to oil-wet can be obtained. Considering that the original sandstones are water-wettable, the contact angle used in the simulation is assigned as [30-60]°. Meanwhile, the permeability saturation curves of both the experiment and simulation results for different wettability are plotted in Fig. 5, which also verifies the feasibility of the equivalent pore network model and the two-phase code.

Fig. 5

Permeability saturation curves of both the EPNMs and the experimental results

Wettability of the reservoir rock directly affects the flow of fluids, and fluids distribution in a reservoir, which is important for transport properties determination. Many experimental investigations on the impact of wettability have been conducted. However, it is impossible to obtain two natural rock cores with the same micro structure and the opposite wettability, that is, it is difficult to realize the single variable control of wettability at the laboratory scale. Simulation of water flooding in the extracted equivalent pore network model can be an effective method to solve this problem.

In this paper, different wettability conditions of the same pore network (ST1 and ST2) are assigned to study the effects on oil recovery. As is shown in Fig. 6, the optimal oil recovery is found in mixed wettability sandstone, when the contact angle is in the interval [80, 100]°. This trend of the simulation results is identical to the experimental result by Morrow et al. (1986) and Jadhunandan et al. (1995). Combined with the permeability saturation curves in Fig. 6, it is found that the water relative permeability and the residual oil saturation increase with the contact angle in [0, 90]° and decrease with the increase of the contact angle in [90, 180]°. In water-wet media, water can flow readily through the wetting layers in the corners of the pore space, while oil is trapped in the larger pores. Then the increase of contact angle leads to a connected oil flowing channel and less trapping. When the medium becomes oil-wet, oil attaches to the solid surface and viscous fingering occurs for water in the pore space, which will reduce the effective channel radius of water and oil recovery.

Fig. 6

Recovery efficiency V.S. Average contact angle θ

Conclusion

This paper presents a modified method of extracting equivalent pore network model from the 3D micro-CT images based on the maximum ball algorithm. The improved partition methods are able to avoid tremendous memory usage when extracting the equivalent pore network model. Two types of micro sandstone images are used to display the pore network models and simulate for single- and two- phase flow. The porosity, calculated permeability of the pore network model agree well with experimental benchmark data of the original sandstone sample. Using the extracted pore network model and two-phase flow codes based on Lattice Boltzmann method, the simulation on water flooding mechanism is conducted to obtain the effects of wettability on oil recovery in the porous sandstone. Visualized water flooding process and the relative permeability saturation curves are obtained. Moreover, it is found that the optimal oil recovery would be realized in the mixed wettability reservoirs. Both of the two simulation results are identical to the experimental benchmark data, which verifies the feasibility of our pore network model and the simulation codes.

However, the assumptions used in the pore network extraction algorithm, which leads to idealized shape, radius and connectivity of pore and throat, result in a lager porosity and permeability of the reconstructed model compared to the original sample. A further research focus on the reproduction of the pore structure for the real shape is worth of studying to realize a wider application in other fields.