Improvement of predictions of petrophysical transport behavior using three-dimensional finite volume element model with micro-CT images＊

LIU Jian-jun (劉建軍), SONG Rui (宋睿), CUI Meng-meng (崔夢夢)

1. State Key Laboratory of Oil and Gas Reservoir Geology and Exploitation, Southwest Petroleum University,Chengdu 610500, China 2. School of Geoscience and Technology, Southwest Petroleum University, Chengdu 610500, China,E-mail: liujj0906@sina.com 3. School of Petroleum and Natural Gas Engineering, Southwest Petroleum University, Chengdu 610500, China

Introduction

An accurate prediction of the material parameters for disordered systems such as rocks, soils, micro emulsions, composites, ceramics, papers, or complex fluids, requires the identification of the pore-scale geometries of these porous media[1-3]. Thus, the multiscale modeling was regarded as a starting point, as emphasized by many scholars[4,5]. Meanwhile, it is possible to conduct experiments to obtain such macroscopic properties. However, the full range of the displacement process is difficult to be determined. For example, experiments on three-phase flow at laboratory scale usually cost much energy and time, but the results are not reliable in the low saturation region.

The pore-scale modeling refers to a numerical simulation method representing the pore space with parameterized geometries, through which the multi flow can be simulated. In 1956, Fatt firstly proposed a twodimensinal regular lattice model and applied it to the prediction of the capillary pressure and the relative permeability of drainage in porous media. After that,many realistic models of porous media were proposed by assuming that the pore space is composed of tube bundles or the solid particles are represented by spheres[6]. These models can reproduce the spatial interconnectivity of pore systems, but cannot reflect the real shapes and distributions of natural porous media.

Developments of the imaging technology help to advance the modeling of porous media greatly, in providing two-dimensinal and three-dimensional imagesof porous media at the resolution of several micrometers. Many algorithms were developed for extracting the pore network from these images, such as the multi-orientation scanning method, the medial axis based method, the Voronoi diagram based method, the maximal ball method and the dual pore network approach[6-13]to construct pore-scale models. In addition,the lattice Boltzmann method is used to simulate the multiphase for predictions of the flow properties. Regular geometries (such as spheres and cylinders) are used to represent pores and throats, therefore, the complex shapes of micro structures in porous media can not be reproduced. In some cases, the pore-scale models were constructed through segmenting the images and extracting the outlines of pore spaces[14,15].These models can reflect the origin porous media images properly and can be used in the commercial simulation code after meshing. But most of these models are two-dimensinal, and can not reflect the spatial interconnectivity of the natural porous media.

To better understand the fundamentals of the multiphase flow in porous media, this paper proposes a novel approach of reconstructing the three-dimensional finite volume elements model from natural porous media images using the commercial software Mimics and ICEM, prior to the pore network model based on some basic assumptions. The tetra finite volume elements of different sizes are employed to represent the origin micro structures of porous media instead of some regular shapes (sphere or cylinder, etc.). The model is tested on two Berea sandstones, four sandstone samples, two carbonate samples, and one Synthetic Silica, which are generated on the micro-CT images scanned by the State Key Laboratory of Oil and Gas Reservoir Geology and Exploitation in the Southwest Petroleum University. Furthermore, single–and two phase flow simulations are conducted based on the Navier-Stokes equations in the Fluent software. Considering that the almost universal practice in the oil industry is to measure the relative permeability for the flow of only two-phases at one initial condition and to use empirical models of limited accuracy to predict the behavior in other cases[16], absolute and relative permeabilities of the models are predicted in the simulation for oil/gas reservoir engineering applications.

1. Pore scale finite volume element modeling

1.1 Micro-CT images preparation

Nine micro-CT images are used in this paper.Considering that the resolution of the Micro-CT images is dependent on the pixel size of the detector and the sample size, the mounted core sample must be small and well-shaped to make the best use of the detector and to acquire the highest possible resolution.Thus the small cylinder samples are drilled out from the core samples to acquire a better resolution, some images of which are shown in the Fig.1. However, it is difficult to construct the surface finite volume element for a cylinder with a complex micro structure. A cube of voxels is extracted from the origin micro-CT images. For example, 450×450×450 pixels of the synthe-tic Silica are selected and used to construct the pore scale model. The basic parameters of the selected micro-CT images used to construct the finite volume element model are listed in the Table 1.

Fig.1 Comparison of the denoised and segmented images

Table 1 Scaled Micro-CT images of core sample

Since the scaled images should be segmented to distinguish the pore space and the matrix in the images, the smoothing and denoising processes are usually employed to reduce sharp angles and small holes in the finite volume element model. As a cornerstone of the modern image processing, the median filter is used here to replace the gray scale value of a voxel with the median value of the nearest 26 surrounding cells. In this way, the images are despeckled from the gray scales without introducing new gray values into the origin. In Fig.1, it is evident that the small holes can be removed effectively by the median filter processing,as witnessed by a comparison between (c) to (d). A larger value of the median filter threshold will remove larger miscellaneous point sets.

1.2 Finite volume element reconstruction

Just like (c) and (d) in Fig.1, the white part with the RGB value of 0 represents the rock matrix and the black part with the RGB value of 255 is the pore space.The rock matrix and the pore space can be distinguished by the threshold value after importing into the Mimics software. The reconstructed models of some samples are shown in Fig.2 (left part), in which the yellow part is the rock matrix and the green part is the pore space. However, these images are completed by a simple surface modeling, to provide only a visual perception and cannot be used in the simulation directly.Traditionally, these models will be transferred into polygon geometries, and then meshed to establish the finite volume elements model. Yet, this traditional process faces two big challenges:

(1) The reconstructed polygon geometry model requires tremendously large storage space in the computer. Taking B2 as an example, the file of the reconstructed geometry is nearly 1GB in size, which will take a huge amount of computer memory and video memory to scan or edit.

Fig.2 Images of reconstructed porous samples in Mimics and finite volume elements of pore space

(2) The meshing process might fail under the existing conditions of the personal computers and the meshing method. The required computer memory for the meshing process will grow exponentially along with the file size, eventually, the existing meshing methods cannot handle the complex disorder microstructures of the porous media.

In view of these difficulties, a new method is proposed, which is explained concisely as follows:

(1) The surface finite mesh modeling is employed instead of the traditional surface geometry modeling, that is, the triangle surfaces meshes are used to reconstruct the surface geometry of the origin porous sample. And the total number of the elements in the models mainly depends on the complexity of the pore shape. Under the circumstance that the skewness (Sk)must not exceed the highest tolerance value of 0.95 for the Fluent software in the remesh process in the Mimics 10.0 software, it is possible to subdivide one triangle of bad quality into more and better triangles.Though unstructured meshes are used in this study,the skewness of these models, shown in Fig.3, is small enough to ensure the reliability for the convergence in the simulation.

針對(duì)不合理信念的調(diào)查，本研究采用王玉（2009）[6]編制的大學(xué)生不合理信念問卷，共包括4個(gè)因子：絕對(duì)化要求、糟糕至極、概括化評(píng)價(jià)和低挫折忍耐；針對(duì)專業(yè)滿意度的調(diào)查，本研究采用柳會(huì)（2014）[7]編制的專業(yè)滿意度調(diào)查問卷，均采用Likert五級(jí)尺度。

Fig.3 Mesh skewness of reconstructed FVM models

Fig.4 Radial pore distribution of origin sandstone core and reconstructed finite volume models. So refers to the origin sandstone core sample

(2) The tetrahedral volume mesh is generated from the original surface mesh by importing the surface meshes into the ICEM software. The volume mesh quality is ensured by a subdivision process. Meanwhile, the boundary meshes are selected to form various parts, preparing for applying the boundary conditions in the numerical simulation. Some of the completed volume meshes are shown in Fig.2 (right part). And the comparative curves of the radial distribution of pores between the reconstructed finite volume models of S1-S4 and the origin sandstone core sample (SO)are shown in Fig.4. It is found that the radius of the pores is close to that of the origin sample, which verifies the feasibility of this method.

2. Numerical simulation

2.1 Mathematical model

For theithphase, the continuity equation is[17]

whereαiis the volume fraction of theithfluid in the cell,ρiis theithfluid’s density,Sαiis the source item. Whenn=1, the equation is just the continuity equation for the single flow.

Fig.5 Boundary conditions in the simulation when the pressure drop is imposed along z direction

The phases volume fraction will be determined based on the following constraint[17]

The properties in the transport equations are determined by the volume fraction of the component phases in each cell. In a multiphase system, the density of the fluid in each cell is given by[17]

Table 2 Absolute permeability of simulation and experiment

All other properties (e.g., viscosity) are computed in this manner.

All phases will share a single momentum equation in the fluid domain, as the result, the velocity field is shared among all phases. The momentum equation, shown below, is dependent on the volume fractions of all phases through the properties contained in Eqs.(3) and (4)[17].

The surface tension for two-phase flow can be regarded as volume force and added to the momentum equation. It has the following form[17]:

Considering the wall adhesion effect, the contact angle between the solid surface and the fluid is adopted to modify the unit normalof phase interface nearby the surface.

By the Fluent software, the outlet flow rate can be obtained. Then the absolute permeability can be calculated as[18]

Then the relative permeability is[18]

whereQiis the total single phase flow rate through the model,Qsiis the total flow rate of the phaseiunder multiphase conditions with the same imposed pressure drop. And bothQiandQsican be obtained by the Fluent software.

2.2 Single-phase flow simulation for permeability prediction

In the single-phase simulation, the pressure drop is imposed on the inlet and the outlet. Taking thezdirection as an example, the top is defined as the inlet and the bottom is the outlet. The other surrounding four faces are defined as impermeable boundaries.The same definitions are in thexandydirections.Since the origin core sample is cylindrical, the experimental permeability data are just acquired along thezdirection. Figure 5 shows the sketch of the boundary conditions in the simulation when the pressure drop is imposed along thezdirection.

By Eq.(9), the absolute permeability along thex,y,zdirections can be calculated. The experiment and simulation results are listed in Table 2. Though there are some differences between the simulation and experiment results for the same type porous sample,they are acceptable for the simulation model is just a part of the origin core sample. Meanwhile, the permeability differs from the simulated direction, which verifies the heterogeneity of the natural porous media.The fitting curves shown in Fig.6 are found to pass though the origin point, which verifies the classic Darcy’s law. The pressure field and the velocity field of the models B1, S1, C2 and SS along thezdirection are shown in Fig.7. It can be seen that the fluid flows mainly along some channels of good connectivity. This phenomenon is verified by the classical micro fluid experiment in porous media.

Fig.6 The variation of volume flow rate with the applied pressure gradient. Linear relationship observed in the graph verifies the Darcy’s law

2.3 Two-phase flow simulation

Based on the finite volume models and the volume of fraction (VOF) model in the Fluent software,the water flooding process is simulated. Here, only the top and bottom faces along thezdirection are selected as the fluid import and export boundaries. The other surrounding four faces are defined as impermeable boundaries. The gradient of 5 MPa/m is imposed between the inlet and outlet boundaries. The Fluid Physical Parameters in the simulation study are listed in Table 3. In the two-phase flow simulation, we choose the models S1, S2, S3 and S4 for the experimental benchmark data of the water flooding in these origin sandstone samples.

The oil volume fraction for different time steps in the water flooding process can be acquired by the Fluent software. The volume fraction images of the model S1 for four different time steps are shown in Fig.8, in which the water moves along the main channels shown in the velocity field of the single-flow simulation. Meanwhile, one may obviously see the viscous fingering phenomenon when the water invades the oil, which causes a hysteresis in the displacement process and some residual oil in some narrow pore spaces. And the oil in the micro fractures of the pore space is impossible to be displaced by the water because of the tremendously large capillary force.

Fig.7 Velocity and pressure fields of four samples along z direction (B1, S1, C2, and Silica)

The oil/water volume fraction and the flow rate can be obtained in every time step, as well as the water saturation and the relative permeability by Eq.(10). The traditional long core-flooding experimentis conducted for the origin sandstone sample in the laboratory. The permeability saturation curves of the simulation and experiment results are shown in Fig.9.Though different contact angles are applied to different models, it is found that the relative permeability of S4 is close to the experimental results of the origin sandstone sample, from which the sample is taken when the contact angle iso90. And the simulation results show that the water saturation of the relative permeability curve intersection changes from ＜0.5 to＞ 0.5 when the sandstone wettability changes from the water-wet to the oil-wet. It is consistent to the qualitative analysis of theory and experiment. All these verify the feasibility and applicability of the reconstructed finite volume model of the porous media and the numerical simulation.

Table 3 Fluid Physical Parameters in the simulation study

Fig.8 Oil volume fraction of model S1 for different steps. From left-top to right-bottom the time steps are 1, 22, 40 and 60

Fig.9 Relative permeability curves of models S1, S2 S3, S4 and experiment

3. Conclusions

In this paper, a novel approach of reconstructing the three-dimensional finite volume elements model from natural porous media images is proposed. This model can reproduce the real shape of pores in porous media, which will better promote the prediction of the petrophysical transport behavior as compared to other equivalent pore network models. The model is tested on two Berea sandstones, four sandstone samples, two carbonate samples, and one Synthetic Silica, which are generated on micro-CT images. With the volume mesh reconstruction process in the Mimics and ICEM software, the mesh qualities of these models are well enough for numerical simulations. In order to validate the efficiency, the models are employed to study the single- and two- phase flows in the Fluent software.The absolute permeability and the relative permeability curve are predicted in the simulation study, which agree well with the experimental results.

It should be noted that this study provides some preliminary results for micro-fluids in porous media using finite volume methods. Though examples used in this study mainly come from the petroleum industry,it can be widely applied to micro fluid and structural analyses in other related areas.

[1] ZHOU Xiao-jun, YU Min-chao. The structural flow in pipe containing porous medium saturated with powerlaw fluid[J]. Journal of Hydrodynamics, 2012, 24(1):138-144.

[2] LI Pei-chao, KONG Xiang-yan and LU De-tang. Mathematical modeling of flow in saturated porous media[J]. Journal of Hydrodynamics, Ser. A, 2003,18(4): 419-426(in Chinese).

[3] SILIN D., TOMUTSA L. and BENSON S. M. et al. Microtomography and pore-scale modeling of two-phase fluid distribution[J]. Transport in Porous Media, 2011,86(2): 495-515.

[4] GRUCELSKI A., POZORSKI J. Lattice Boltzmann simulations of flow past a circular cylinder and in simple porous media[J]. Computers and Fluids, 2013, 71:406-416.

[5] JAMSHIDZADEH Z., TSAI F. T. C. and MIRBAGHERI S. A. et al. Fluid dispersion effects on density-driven thermohaline flow and transport in porous media[J].Advances in Water Resources, 2013, 61: 12-28.

[6] BAUER D., YOUSSEF S. and HAN M. et al. From computed microtomography images to resistivity index calculations of heterogeneous carbonates using a dualporosity pore-network approach: Influence of percolation on the electrical transport properties[J]. Physical Review E, 2011, 84(1): 011133.

[7] LIANG Z., IOANNIDIS M. A. and CHATZIS I. Geometric and topological analysis of three-dimensional porous media: Pore space partitioning based on morphological skeletonization[J]. Journal of Colloid and Interface Science, 2000, 221(1): 13-24.

[8] SHIN H., LINDQUIST W. B. and SAHAGIAN D. L. et al. Analysis of the vesicular structure of basalts[J].Computers and geosciences, 2005, 31(4): 473-487.

[9] BRYANT S., BLUNT M. Prediction of relative permeability in simple porous media[J]. Physical Review A,1992, 46(4): 2004-2011.

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

[11] ?REN P. E., BAKKE S. Reconstruction of Berea sandstone and pore-scale modelling of wettability effects[J].Journal of Petroleum Science and Engineering, 2003,39(3): 177-199.

[12] SILIN D., PATZEK T. Pore space morphology analysis using maximal inscribed spheres[J]. Physica A: Statistical Mechanics and its Applications, 2006, 371(2):336-360.

[13] SILIN D., JIN G. and PATZEK T. Robust determination of the pore space morphology in sedimentary rocks[C]. SPE Annual Technical Conference and Exhibition. Houston, USA, 2003.

[14] GUNDE A. C., BERA B. and MITRA S. K. Investigation of water and CO2(carbon dioxide) flooding using micro-CT (micro-computed tomography) images of Berea sandstone core using finite element simulations[J]. Energy, 2010, 35(12): 5209-5216.

[15] LIU J., SONG R. and ZHAO J. Numerical simulation research on seepage mechanism in pore-scale deformable porous media[J]. Disaster Advances, 2013, 6(S1):49-58.

[16] BLUNT M. J. Flow in porous media-pore-network models and multiphase flow[J]. Current Opinion in Colloid and Interface Science, 2001, 6(3): 197-207.

[17] YANG Z., PENG X. F. and YE P. Numerical and experimental investigation of two phase flow during boiling in a coiled tube[J]. International Journal of Heat and Mass Transfer, 2008, 51(5): 1003-1016.

[18] VALVATNE P. H., BLUNT M. J. Predictive pore-scale modeling of two-phase flow in mixed wet media[J].Water Resources Research, 2004, 40(7): 1-21.