Influence of pore structure parameters on flow characteristics based on a digital rock and the pore network model

Journal of Natural Gas Science and Engineering 31 (2016) 156e163

Contents lists available at ScienceDirect

Journal of Natural Gas Science and Engineering journal homepage: www.elsevier.com/locate/jngse

Influence of pore structure parameters on flow characteristics based on a digital rock and the pore network model Senyou An, Jun Yao, Yongfei Yang*, Lei Zhang, Jianlin Zhao, Ying Gao School of Petroleum Engineering, China University of Petroleum, Qingdao 266580, China

a r t i c l e i n f o

a b s t r a c t

Article history: Received 4 November 2015 Received in revised form 24 February 2016 Accepted 2 March 2016 Available online 4 March 2016

Many approaches have been developed to study the role that pore structure parameters play in flow characteristics. In this paper, regular network models are used to study the influence of the porosity, particle size, pore throat ratio, coordination number, shape factor and pore throat orientation on the absolute permeability. Among them, every parameter has obvious effect except shape factor when single-phase fluid flows in the models. Then the random pore network model is developed to analyze the relative permeability curve with a changing pore throat ratio, coordination number and shape factor. When the water saturation is constant, the oil relative permeability is higher for models that have a higher coordination number or have a lower pore throat ratio, and the wet phase has low permeability because of the corners. In order to study the dominant factor of the pore structure for two-phase flow further, a digital core, based on CT scanning, is constructed to firstly research the capillary influence of internal structure that determines whether the fluid can flow at a certain pressure difference. © 2016 Elsevier B.V. All rights reserved.

Keywords: Pore structure Flow simulation Pore network model Digital core

1. Introduction The structure parameters of a porous medium define its absolute permeability, relative permeability and whether flow occurs at a certain pressure difference. For a tight reservoir, in particular, the structures are various and the flow limit is of particular concern (Yao et al., 2013). To clarify the flow mechanism with different transmission parameters, as early as 1856, French engineer Darcy did massive experiments, and obtained relations between the pressure gradient and the seepage velocity, namely for Darcy Law(Bell, 1983), where the permeability k contains the effect of pore structures for flow. With the deepening of the study and the level of research to improve, some scholars have begun to consider the internal structure of porous media. Firstly, the single pipe model was used to study this phenomenon at the pore scale (Carman, 1937; Kozeny, 1927). Next, the classical capillary bundle model was used to accurately represent the porous media more (Purcell, 1949); the model can be used to explain wetting hysteresis. Other researchers improved this model further by defining the concept of the pseudo mean radius and tortuosity. In 1952, famous Ergun equation was proposed by Ergun (1952) and this equation

* Corresponding author. E-mail address: [email protected] (Y. Yang). http://dx.doi.org/10.1016/j.jngse.2016.03.009 1875-5100/© 2016 Elsevier B.V. All rights reserved.

have considered the relations between permeability and porous medium porosity and equivalent mean sphere diameter. Subsequently, statistical models and simple regular pore networks were gradually proposed (Childs and Collis-George, 1950; Dullien, 1975). The relations between permeability and pore structures have been proposed quantitatively. Currently, there are two main approaches for modeling porescale transport: one method is a direct method that is usually solved by NeS equation or Boltzmann equation based on scanning images or reconstructed images; the other method is the pore network model method, which is an approximate way for real structure. In direct modeling, CT is a higher accuracy approach, that was first applied to describe reservoir characteristics by Dunsmuir et al. (1991), compared with the numerical reconstruction method and the 2D scanning electron microscope method. The pore network model is similar to the single pipe and capillary bundle model in principle (Fatt, 1956). The network of pores and throats with simple geometric shapes can be extracted from the 3D image of a core (Blunt et al., 2013); methods that accomplish this extraction include the grain-based method (Bakke and Øren, 1997), the medial-axis based method (Sheppard et al, 2005) and the maximal-ball based method (Dong and Blunt, 2009). In addition, a random algorithm can be used to simulate flow by changing the structure parameters. In direct modeling, the modeled structure is very close to the real core, but most of the direct methods are very

S. An et al. / Journal of Natural Gas Science and Engineering 31 (2016) 156e163

computationally demanding when addressing multi-phase fluid flow in porous media (Meakin and Tartakovsky, 2009; Raeini et al, 2012), such as the finite element method (FEM) and the LatticeBoltzmann method (LBM) (Gunstensen et al., 1991). The pore network model is a better approach in computational cost for multi-phase flow, so pore network model is usually constructed based on real pore structure to analyze relative permeability. Arns et al. (2004) introduced the importance of incorporating realistic 3D topologies in network models for predicting multiphase flow properties. The parameters of pore network model obtained by this approach cannot be adjusted. Thus, random three-dimensional pore network model is used in this paper to research the influence of pore structure parameters. After modeling, a suitable equation should be chosen to control the fluid flow, which is usually based on Knudsen number (Kn) for gas reservoir and is defined as the ratio of the molecular mean free path to the characteristic length. The Knudsen number is useful for determining whether statistical mechanics or the continuum mechanics formulation of fluid dynamics should be used: If the Knudsen number is near or greater than one, the mean free path of a molecule is comparable to a length scale of the problem, and the continuum assumption of fluid mechanics is no longer a good approximation. In this case statistical methods must be used(Colin, 2005). For conventional oil reservoir, the size of the pores is much greater than the size of fluid particle or element, which is far more than fluid molecules (Adler, 1992), so the fluid in a porous media can be regarded as a continuous medium and the NeS equation is valid. Due to the limitation of the number of computer calculations, many researchers simplify the NeS equation into the Poiseuille formula in the simulation. In this paper, we directly study the influence of structure parameters on the flow characteristics of a transmission based on the NeS equation for simple pore network models and real digital model. In addition, three methods with different strengths are used to study the absolute permeability, the relative permeability and the flow limit. The simple regular network models are designed to research absolute permeability based on NeS equation and random pore network is developed to analyze relative permeability. Besides, the real digital core model is firstly used to study the capillary influence of internal pore structure.

2. Mathematical model and validation 2.1. Flow equation

157

chose to be weight function of continuity equation and momentum equation respectively. Galerkin weaker form of NeS equation is derived firstly based on Green equation. And then, total finite element equation is deduced as following: ðeÞ

ðeÞ

ðeÞ

ðeÞ

Anr uar þ Bnbrs uar ubs þ Cant Pt þ Danbs ubs ¼ Ean ðeÞ

ðeÞ

Fmbr ubr ¼ Gm

(3) (4)

where, a ¼ 1,2; n ¼ 1,2,,,,,,Nu; m ¼ 1,2,,,,,,Np; b is 1, 2and take the sum; r,s is 1,2,,,,,,Nu and take the sum; t is 1,2,,,,,,Np and take the sum. Nu and Np are total number of velocity nodes and pressure nodes respectively. The symbols in above equation are listed as following:

Z

ðeÞ

Anr ¼

rFn Fr ,dU UðeÞ

Z

ðeÞ

Bnbrs ¼

rFn UðeÞ

Z

ðeÞ

Cant ¼

dab UðeÞ

Z

ðeÞ

Danbs ¼

m UðeÞ

ðeÞ


 vFn vFs vFs dab ,dU þ drb vxr vxr vxa Z

rfa Fn ,dU þ U

ðeÞ

Z

Fmbr ¼ UðeÞ ðeÞ

vFn Jt ,dU vxb

Z

Ean ¼ ðeÞ

vJr Fs ,dU vxb

pam Jm Fn ,dG ðeÞ G2

vJm Fr ,dU vxb

Z

Gk ¼

unr Fr Jk ,dG ðeÞ

G1

where, Fj is unit interpolatory basic function of velocity; Jk is unit interpolatory basic function of pressure; Fj is one order higher than Jk; d represents weight function; G is natural boundary. When the flow is steady flow, finite element equation can be expressed as following. This equation is solved by linearized iterative method.

32 3 2 3 u1s Bn1rs u1r þ D1n1s Bn1rs u1r þ D1n1s C1nt E1n 4 Bn1rs u2r þ D2n1s Bn2rs u2r þ D2n2s C2nt 54 u2s 5 ¼ 4 E2n 5 Fm1 Fm2 0 pt Gm 2

The fluid flowing in a porous media can be regarded as a continuous fluid. Considering the real condition, the liquid can be idealized as an uncompressed fluid. When the gravity is neglected, NeS equation can be deduced as follows:

(5)

2.2. Model validation


h  i vu þ ðu,VÞu þ Vp  V, m Vu þ VuT ¼0 r vt

(1)

V,u ¼ 0

(2)

where u is the flow velocity for one point and r is the corresponding density. This equation with highly nonlinear variables is typically a second-order non-self-adjoint equation (Yuan and Lin, 2007; Jian and Chao, 2010). It is difficult to find the appropriate variational principles for this equation, so the Galerkin finite element method is used to solve this function (Ye and Liu, 2008). The interpolatory basic function of pressure p and velocity ui are

To validate the accuracy of the above-described numerical model that describes the control of the fluid through porous media, the physical model is designed as shown in Fig. 1. The parallel circles of infinite extent represent solid cylinders, with fluid flowing in the gray area from the left boundary to the right boundary. Constant pressures are assigned for the inlet and outlet boundaries. The symmetrical boundary condition is used for the top and bottom sides. This two-dimensional model is equivalent to a threedimensional cylinder structure due to the identical characteristics that are found in every slice. In adjusting the porosities of models, the radius is set as 5 mm, 10 mm, 15 mm, 20 mm, 25 mm, and 30 mm respectively, and each

158

S. An et al. / Journal of Natural Gas Science and Engineering 31 (2016) 156e163

Fig. 1. Schematic diagram of the validation model.

cylinder is assigned a number of spaces. Fig. 2 and Fig. 3 show the velocity and pressure distributions when the radius is 10 mm and the distance is 10 mm; these distributions are found to be quite symmetric. Both the simulation result and the analytical solution are presented in Fig. 4, which shows the change of the normalized permeability with porosity. From Gebart (1992), the analytical solution for this model is the following equation:

sffiffiffiffiffiffiffiffiffiffiffiffiffiffi !5=2 K 1  fc ¼C 1 1f a2

(6)

where K is the absolute permeability; a is the radius of the solid particle; C is the geometric factor; and fc is critical value of porosity below which there is no permeating pffiffiffi flow. For the square array of cylinders, the value of C is 16 9 2p, and fc is equal to 1  p 4 (Gebart, 1992). The limit of this analytical equation is f < 0.65. Fig. 4 indicates that the simulation results are highly consistent with the analytical solution. In other words, the flow equations and code are accurate, especially when 0.3 < f < 0.65. =

=

3. The influence of the pore structure on the absolute permeability According to the validating model, we have obtained the change of the absolute permeability with porosity. Next, the impact of the pore throat ratio, coordination number and shape factor on the absolute permeability are analyzed as follows. 3.1. The impact of the pore throat ratio The pore throat ratio is defined as the ratio of the pore diameter to the throat diameter, as denoted by Rpt. lp is used to describe the greatest length of the pore slice, and the shortest length of the throat slice is lt. Thus, the pore throat ratio of Fig. 5 is given by the following equation:

Fig. 3. The pressure image.

Rpt ¼

lp lt

(7)

To analyze the influence, the three-dimensional mode is designed as shown in Fig. 6. The porosity of different models is set an approximately constant value, being about 14%. Changing the diameters of the pores and the throats at the same porosity, we obtain Fig. 7 through flow simulations. The absolute permeability decreases exponentially as the pore throat ratio increases. The difference in size between the pore and the throat becomes prominent when the pore throat ratio is higher. In other words, the heterogeneity of the model is more obvious. In addition, the ratio of small pores relatively increases when the porosity is maintained at a constant value. All of these factors cause the flow capacity to decrease.

3.2. The impact of the coordination number and the pore throat orientation The coordination number is the average number of throats linked with one pore in the model, which is a topology parameter. The coordination number of the isolated pore is zero, and for a blind side pore, only one side is connected with a throat, so the number is 1. The range of the coordination number is usually 2e14 and mostly fluctuates between 3 and 6. In this section, we construct the models as shown in Fig. 8, all of

10 10

Analytical result Simulation result

10

K/a

2

10 10 10 10 10 10 0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1.0

Porosity Fig. 2. The velocity image for the validation model.

Fig. 4. Comparison between the simulation result and the analytical result.

S. An et al. / Journal of Natural Gas Science and Engineering 31 (2016) 156e163

159

Fig. 5. The connection of the pore and throat.

Fig. 8. Physical models for studying the coordination number.

Fig. 6. Three-dimensional model for studying the pore throat ratio.

which have the same diameters of pores and throats (5 mm and 15 mm, respectively), and the porosities are similar. The coordination numbers are set differently for model (1) and model (2). In addition, model (3) is designed to be different regarding the orientation of pores and throats. The flow of the unit thickness can be obtained by using the integral for the outlet end, which is 15509.29 mm2/s, 16812.10 mm2/s and 12207.02 mm2/s. For the fluid flow in the space from the left boundary to the right boundary under a pressure difference of 1000 Pa, we can calculate the permeability, being 0.0217 mm2, 0.0236 mm2, 0.0172 mm2 respectively. Analyzing the results of model (1) and model (2), we arrive at the conclusion that when the orientation of pores and throats is similar, the higher the coordination number, the better the flow capacity. The cause of this condition is the improved connection for models with higher coordination numbers. However, if the pore throat orientation is different, the result will be not comparable, as is clarified by model (2) and model (3). Their coordination numbers

and porosities are the same, but the differences of their flow velocities are obvious. Model (2) has the best condition, whereas model (3) is the worst one. Thus, when we evaluate the geological quality of a reservoir, the orientation cannot be neglected. 3.3. The impact of the shape factor The shapes of the pores and throats are various because of the complex diagenetic environment. Mason and Morrow defined a shape factor to describe different sectional shapes. The shape factor is defined as follows:

G¼

A P2

(8)

Where A is the sectional area, m2; P is the section perimeter, m; the shape factor is 1/16 for a square; and G ¼ 1/(4p ) for a cycle. The range for a triangle is from 0 to 0.0481, and the maximum value is

300

Permeability/mD

250

200

k = 520.52e-0.608Rpt

150

100 50 0 0

1

2

3 4 Pore throat ratio

5

6

7

Rpt

Fig. 7. The change of permeability with the pore throat ratio.

Fig. 9. The models with different shape factors.

160

S. An et al. / Journal of Natural Gas Science and Engineering 31 (2016) 156e163

1.0 0.9 0.8

0.0796

Relative permeability

0.7

0.0625

0.6

kro krw

0.0481

0.5 0.4 0.3

0.2 0.1 0.0 0

Fig. 10. Three-dimensional stochastic pore network model.

0.2

0.4 0.6 Water Saturation(Sw)

0.8

1

Fig. 13. Relative permeability curves of models with different shape factors.

1.0 0.9

Relative permeability

4. The influence of the pore structure on the relative permeability

4

0.8

3

0.7

kro

2.5

0.6

1.5

0.5 0.4

krw

0.3 0.2 0.1 0.0 0

0.2

0.4 0.6 Water saturation(Sw)

0.8

1

To some extent, the absolute permeability can be used to evaluate the flow capacity of the porous media. However, the proportion of a single flow is small in actual production. We study the influence of some parameters on the relative permeability in this chapter. The impact of the pore throat ratio, coordination number and the shape factor are simulated by a random three-dimensional pore network model(Mohammad and Blunt, 2005), as shown in Fig. 10; these three parameters are changed in the codes. The pore network model with simple regular geometric shapes is a simplification of the actual pore geometry, and the simulation of flow is easier to be calculated.

Fig. 11. Relative permeability curves of different pore throat ratio.

4.1. The impact of the pore throat ratio obtained for an equilateral triangle. Three types of models are designed to obtain the law of absolute permeability versus the shape factor, as shown in Fig. 9. All of these models have the same coordination number (4), porosity, and pore throat orientation. Only the shape factor is changed. The flow simulation is operated at a pressure difference of 1000 Pa. We obtain corresponding permeability of 0.0236 mm2, 0.0235 mm2 and 0.0229 mm2 from the models, which tells us that the influence of the shape factor on absolute permeability can be neglected. 1.0 0.9

2.5

Relative permeability

0.8

kro

3

0.7

Firstly, the pore throat ratio of the pore network models is designed to have values of 4, 3, 2.5 and 1.5 for the flow simulations. All of the models have the same absolute size (500 mm  500 mm), porosity (28%), pore throat ratio, coordination number and wetting angle. In the simulation, the models saturated with oil are flooded by water. The relative permeability curve can be determined by monitoring the changes of oil and water saturations, as shown in Fig. 11. When the water saturation is constant, the oil relative permeability is higher for models that have a lower pore throat ratio, and the oil relative permeability decreases more slowly. That means the common area (water and oil relative permeability are all greater than zero in relative permeability curve) is greater correspondingly. The more serious heterogeneity results that more oil cannot be moved at same pressure difference. To observe inner phenomenon,

3.5

0.6

4

0.5

krw

0.4 0.3 0.2 0.1

0.0 0

0.2

0.4

0.6

0.8

1

Water saturation(Sw) Fig. 12. Relative permeability curves of different coordination numbers.

Fig. 14. Gray image of a sandstone core via a CT scan.

S. An et al. / Journal of Natural Gas Science and Engineering 31 (2016) 156e163

161

corresponding simulations. After flooding, the relative permeability curve can be obtained, as shown in Fig. 12. When the water saturation is constant, the oil relative permeability is higher for models that have a higher coordination number, and the oil relative permeability decreases more slowly. The corresponding situation is that the residual oil saturation is smaller for a higher coordination number model because the water phase prefers to stay in the corners for water wet cores, thereby improving the chance to enter more space. The common area of water and oil is greater in the curve when the connection is better. Another obvious phenomenon is that the impact of the change in coordination number is greater for the oil phase as water is attached to the solid surface. Fig. 15. The distribution of particles after binary segmentation.

Fig. 16. Meshing result of the geometry.

single throat model based digital core is studied in later content.

4.3. The impact of the shape factor In addition, the complex shapes of pores and throats have different impacts on the flow capacity of fluid for rocks with diverse wettability because of the corners, whereas the absolute permeability is almost independent of the shape factor. The change of the number of the shape factor in the models has an obvious influence on the relative permeability curve, as shown in Fig. 13. Analyzing the curve, we can see that the model that has cycle pores has smallest residual oil saturation; at the same time, the space of triangle pores are primarily occupied by remaining oil. When the water saturation is the same, the relative permeability of the wet phase, that is, water, is greater with a larger shape factor in the pore network model, whereas the capacity of the non-wetting phase flow is weaker because the corner is more prominent in the model with a small shape factor; in addition, the solid wall is water wet, which results in water tending to be continuous in the corners(Xu et al., 2014). Oil wrapped in the water cannot easily flow because of the Jamin effect, so the residual oil saturation will increase.

4.2. The impact of the coordination number 5. The capillary influence of internal pore structure Coordination number is an average number, which is equal to total number of throats over the number of pores in the random network model. Then, we set the coordination number of the pore network models to have values of 2.5, 3, 3.5, and 4 for the

As just mentioned, the structure of the core is quite complex, even for one flow channel, with the sectional size changing in different locations. In section 4.1, we just know that the relative

Fig. 17. The map of the volume fraction after displacement at a differential pressure of 0.4 MPa.

162

S. An et al. / Journal of Natural Gas Science and Engineering 31 (2016) 156e163

To determine the dominant factor, single-phase flow is simulated because the velocity field can reveal the pore structure. In the simulation, the inlet and outlet are designed diagonally to guarantee that fluid can flow throughout the space. The distribution of velocity is shown in Fig. 18, from which we can find a high velocity channel, the diameter of which can be calculated as follows.

r¼

rffiffiffiffiffiffi rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Q 193:12 ¼ ¼ 1:368mm pv 3:14  37:24

(9)

where Q is the flow, which can be obtained by integration for the outlet, and v is the highest velocity in the whole space. Thus, the greatest capillary force can be obtained when the interface force is known.

Pc ¼

Fig. 18. The velocity distribution of single-phase displacement.

permeability changes with pore throat ratio, but what is internal factors? In this chapter, we study the dominant factor of the pore structure for two-phase flow by using phase-field method(Bao et al., 2012) based on a digital core. First, a gray image of a sandstone core is obtained via CT scanning, as shown in Fig. 14; the gray image can approximately reflect the internal structure of the core (Dunsmuir et al., 1991). Next, the reconstructed cylinder is intercepted and segmented by binarization, as shown in Fig. 15 (0 represents the skeleton, and 1 represents a pore). A smaller volume with various diameters is further selected. To simulate flow, the data are translated into geometry (Zhao, 2009) and then are meshed, as shown in Fig. 16 (Fabri and Pion, 2009). The model that has a wall that is absolutely lipophilic is occupied by oil. Water floods this model at a pressure difference of 0.4 MPa from t. The result is shown in Fig. 17, which tells us that only part of the space is filled with water.

2s,cos q 2  0:031  1 ¼ ¼ 45321:6Pa r 1:368  106

(10)

The above calculation can explain the flooding result with a pressure difference 0.4 MPa. Next, we adjust the pressure difference to a larger degree to analyze the critical flow condition. When the value arrives at 0.7 MPa, the water spreads almost throughout the entire space as shown in Fig. 19, without blind-side pores, that is, the greatest capillary force is overcome. Considering the pressure decrease in the process of flow (slightly larger than 0.2 MPa from inlet to the smallest throat), the pressure difference across the smallest space just exceeds the capillary force. We can arrive at the conclusion that the capillary pressure of the smallest throat is the dominant factor. Besides, water or oil exists when a gas reservoir is produced and the surface tension is quite large. Thus, we must design a reasonable minimum pressure difference to ensure costeffective production. 6. Conclusion Three approaches were used in this paper to help us understand the relationship between the flow capacity and the pore structure parameters. The absolute permeability was analyzed using a regular pore network model because the flow characteristics can be simulated clearly in this way. In addition, the NeS equation was chosen due to its applicability and accuracy, and the solution was validated by fitting Gebart formula. We found that the absolute

Fig. 19. The map of the volume fraction after displacement at a differential pressure of 0.7 MPa.

S. An et al. / Journal of Natural Gas Science and Engineering 31 (2016) 156e163

permeability decreases exponentially with the increase of the pore throat ratio. When the orientation of pores and throats is similar, the higher the coordination number, the better the flow capacity, and the pore throat distribution is also a critical factor; however, the influence of the shape factor to absolute permeability can be ignored. The direct method has the advantage of high accuracy, but it needs great computational cost for multi-phase flow, which is difficult to achieve. And we cannot change structure parameters of real core model as our design. Thus, a random pore network model was adopted to simulate two-phase flow in this paper. The pore throat ratio, coordinate number and shape factor are main structure parameters for the relative permeability curve. When the water saturation is constant, the oil relative permeability is higher for models that have a higher coordination number, and the wet phase has the lowest permeability because of the corners. By observing relative permeability curve, the oil relative permeability is higher for models that have a lower pore throat ratio when the water saturation is constant, and the oil relative permeability decreases more slowly. That means the common area greater correspondingly. The heterogeneity results that more oil cannot be moved at same pressure difference. In real production, engineers pay attention to whether oil can flow at a certain pressure difference and how large of a pressure difference should be set for effective development. Pore network model just tells us that the relative permeability changes with pore throat ratio. The digital core is a best choice for studying internal structure, because the real pore structure can be selected and we can clearly observe the internal structure. After conducting simulations at various pressure differences, it was found that the capillary force of the minimum connected diameter is the dominant parameter considering the flow pressure drop. The results reveal the influence of the pore structure parameters, which can be referenced for oil reservoir production. In addition, the minimum flow pressure theory is a principle to guide the production of water flooding oil. Acknowledgments We express appreciation for the following financial support: the National Natural Science Foundation of China (Nos. 51304232, 51490654, 51234007), the Major Programs of Ministry of Education of China (No. 311009), the Specialized Research Fund for the Doctoral Program of Higher Education (20120133120017), the Fundamental Research Funds for the Central Universities (No. 14CX05026A), the China Postdoctoral Science Foundation (No. 2015M580621), the Introducing Talents of Discipline to Universities (B08028), and the Program for Changjiang Scholars and Innovative Research Team in University (IRT1294). References Adler, P.M., 1992. Porous Media: Geometry and Transports. Arns, J.-Y., Robins, V., Sheppard, A.P., Sok, R.M., Pinczewski, W.V., Knackstedt, M.A.,

163

2004. Effect of network topology on relative permeability. Transp. Porous media 55 (1), 21e46. Bakke, S., Øren, P.-E., 1997. 3-D pore-scale modelling of sandstones and flow simulations in the pore networks. Spe J. 2 (02), 136e149. Bao, K., Shi, Y., Sun, S., Wang, X.-P., 2012. A finite element method for the numerical solution of the coupled CahneHilliard and NaviereStokes system for moving contact line problems. J. Comput. Phys. 231 (24), 8083e8099. Bell, J., 1983. Hydrodynamics in Porous Media. Chen Chongxi translate. China: Building industry Publishing house, pp. 91e94. Blunt, M.J., Bijeljic, B., Dong, H., Gharbi, O., Iglauer, S., Mostaghimi, P., Paluszny, A., Pentland, C., 2013. Pore-scale imaging and modelling. Adv. Water Resour. 51, 197e216. http://dx.doi.org/10.1016/j.advwatres.2012.03.003. Carman, P., 1937. Review of literature, confirming validity of Kozeny's theory. Trans. Inst. Chem. Eng. 15, 150. Childs, E.C., Collis-George, N., 1950. The permeability of porous materials. Proc. R. Soc. Lond. Ser. A, Math. Phys. Sci. 201 (1066), 392e405. http://dx.doi.org/ 10.2307/98381. Colin, S., 2005. Rarefaction and compressibility effects on steady and transient gas flows in microchannels. Microfluid. Nanofluidics 1 (3), 268e279. Dong, H., Blunt, M.J., 2009. Pore-network extraction from micro-computerizedtomography images. Phys. Rev. E 80 (3), 036307. Dullien, F.A.L., 1975. New network permeability model of porous media. AIChE J. 21 (2), 299e307. http://dx.doi.org/10.1002/aic.690210211. Dunsmuir, J.H., Ferguson, S., D'Amico, K., Stokes, J., 1991. X-ray microtomography: a new tool for the characterization of porous media. In: SPE Annual Technical Conference and Exhibition. Society of Petroleum Engineers. Ergun, S., 1952. Fluid flow through packed columns. Chem. Eng. Prog. 48. Fabri, A., Pion, S., 2009. CGAL: the computational geometry algorithms library. In: Proceedings of the 17th ACM SIGSPATIAL International Conference on Advances in Geographic Information Systems. ACM, pp. 538e539. Fatt, I., 1956. The network model of porous Media. I, II and III. Petrol. Trans. AIME 207. Gebart, B., 1992. Permeability of unidirectional reinforcements for RTM. J. Compos. Mater. 26 (8), 1100e1133. Gunstensen, A.K., Rothman, D.H., Zaleski, S., Zanetti, G., 1991. Lattice Boltzmann model of immiscible fluids. Phys. Rev. A 43 (8), 4320. Jian, Y., Chao, Y., 2010. Study on discontinuous Galerkin method for NaviereStokes equation. Chin. J. Theor. Appl. Mech. 42 (5), 962e970. Kozeny, J., 1927. Über kapillare Leitung des Wassers im Boden:(Aufstieg, Versick€lder-Pichler-Tempsky 136, erung und Anwendung auf die Bew€ asserung). Ho 271e306. Meakin, P., Tartakovsky, A.M., 2009. Modeling and simulation of pore-scale multiphase fluid flow and reactive transport in fractured and porous media. Rev. Geophys. 47 (3). Mohammad, P., Blunt, M.J., 2005. Three-dimensional mixed-wet random pore-scale network modeling of two- and three-phase flow in porous media. I. Model description. Phys. Rev. E 71 (2), 199e201. Purcell, W.R., 1949. Capillary pressures-their measurement using mercury and the calculation of permeability. Trans. AIME 39 (1949). Raeini, A.Q., Blunt, M.J., Bijeljic, B., 2012. Modelling two-phase flow in porous media at the pore scale using the volume-of-fluid method. J. Comput. Phys. 231 (17), 5653e5668. http://dx.doi.org/10.1016/j.jcp.2012.04.011. Sheppard, A., Sok, R., Averdunk, H., 2005. Improved pore network extraction methods. In: International Symposium of the Society of Core Analysts, pp. 21e25. Xu, W., Ok, J.T., Xiao, F., Neeves, K.B., Yin, X., 2014. Effect of pore geometry and interfacial tension on water-oil displacement efficiency in oil-wet microfluidic porous media analogs. Phys. Fluids 26 (9), 093102. http://dx.doi.org/10.1063/ 1.4894071. Yao, J., Sun, H., Huang, Z., Zhang, L., Zeng, Q., Sui, Y., Fan, D., 2013. The key mechanical problems in the development of Shale Gas. Sci. CHINA Phys. Mech. Astronomy 1527e1547. Ye, L., Liu, J., 2008. Numerical Simulation of Microcosmic Seepage in Porous Media Baseed on N-s equation. Wuhan Polytechnic University. Yuan, S., Lin, Y., 2007. An EEP method for post-computation of super-convergent solutions in one-dimensional Galerkin FEM for second order non-self-adjoint boundary-value problem. Chin. J. Comput. Mech. (02), 142e147. Zhao, X., 2009. Numerical Rock Construction and Pore Network Extraction. China University of Petroleum, Dong Ying.