About this Journal Submit a Manuscript Table of Contents
Advances in Mechanical Engineering
Volume 2013 (2013), Article ID 137086, 11 pages
Research Article

On Full-Tensor Permeabilities of Porous Media from Numerical Solutions of the Navier-Stokes Equation

Yi Wang,1,2 Shuyu Sun,2 and Bo Yu1

1National Engineering Laboratory for Pipeline Safety, Beijing Key Laboratory of Urban Oil and Gas Distribution Technology, China University of Petroleum, Beijing 102249, China
2Computational Transport Phenomena Laboratory, Division of Physical Science and Engineering, King Abdullah University of Science and Technology, Thuwal 23955-6900, Saudi Arabia

Received 25 January 2013; Revised 28 May 2013; Accepted 3 June 2013

Academic Editor: Toshio Tagawa

Copyright © 2013 Yi Wang et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.


A numerical method is proposed to compute full-tensor permeability of porous media without artificial simplification. Navier-Stokes (N-S) equation and Darcy’s law are combined to design these numerical experiments. This method can successfully detect the permeability values in principle directions of the porous media and the anisotropic degrees. It is found that the same configuration of porous media may possess isotropic features at lower Reynolds numbers while manifesting anisotropic features at higher Reynolds numbers due to the nonlinearity from convection. Anisotropy becomes pronounced especially when convection is dominant.

1. Introduction

Permeability is a key parameter of porous media because it relates with many parameters, such as the infiltration, dielectric strength and thermophysics. It may strongly affect the design and production process of fibrous reinforcements, cement paste and so forth and related scientific researches. The permeability may vary over several orders of magnitude, and it plays an important role in petroleum migration and reservoir performance [1]. Therefore, the accurate measurement of permeability is always a hot topic in both industrial and academic fields. Many of these researches adopt experimental approaches, such as free-space methods, open-ended coaxial probe techniques, cavity resonators, transmission-line techniques, and gas-dynamic methods [2, 3].

Jensen and Heriot-Watt [4] proposed a statistical model for small-scale permeability using minipermeameter and core plug measurements. They suggested the minipermeameter measurement as a better choice. Wang et al. [5] discussed two experimental methods to determine the absolute values of in-plane permeability. They concluded that both the radial flow measurement method and the unidirectional flow measurement method were recommended to obtain reliable permeability data. Ferland et al. [6] proposed a concurrent method to estimate permeability at low cost. However, the experimental procedures can introduce uncertainty significantly on the estimated permeability. They suggested some special treatments to increase the reliability of experimental data. Some other researchers [710] also realized that permeability is difficult to measure although its definition is simple. Many factors, such as flow rate, pressure, fluid properties, handling process by human factors, edge effect, wall effect, single-equipment reproducibility, and between-equipment repeatability, could strongly influence the results of measurements. Therefore, the ability to obtain consistent permeability data depends on skilled and experienced experimental design, reproducible preparation of specimens, operation of equipment, and evaluation of measured raw data. This is the reason many published data observed by different persons often have significant differences for the same material. Weitzenböck et al. [11] pointed out that numerous practical problems caused three-dimensional (3D) permeability measurements to be very difficult. Then, they proposed an approach to measure the permeability by a two-dimensional (2D) radial flow method, which allows the experimental axes not to align with the principal direction of permeability [12, 13].

Some researchers made efforts on measuring transport properties (including permeability) by numerical simulation. Keehm et al. obtained good results in estimating permeability and electric conductivity of complicated pore geometries using Lattice-Boltzmann method (LBM) [14]. Kameda et al. applied LBM to estimate 3D permeability through 2D images of small fragments of rocks. They obtained a valid permeability-porosity trend by using a significant number of such small fragments in statistical sense [15]. Saenger et al. estimated permeability through LBM flow simulation and compared mechanical and transport properties for the same digital rocks [16].

Since all physical parameters can be easily and exactly controlled in numerical simulations, the shortcomings of laboratory experiments mentioned previously can be avoided. In particular, numerical methods have the advantage of separately studying different physical processes coexisting in nature, which are uneasily separated in laboratory. It is easier to form standardization of measurement methods and obtain unchangeable results. It will be a very efficient and economical way for measuring properties of porous media combining with the digital rock physics [17, 18]. In the present paper, we studied how to establish this digital laboratory method for measuring permeability of porous media through directly solving the N-S equation, other than reconstructing it by LBM method. To the best knowledge of the authors, all studies of permeability prediction have been concerned only with diagonal permeability tensor actually. Full tensor of permeability has not been studied extensively since it has more general meanings for practical applications. Thus, we expand our research object from the simple diagonal permeability tensor to more general full permeability tensor. We select gravity in periodic domain as the driven force, instead of pressure widely used in previous studies [118]. This can avoid the edge effect encountered in experimental measurements. As a first trial, we assume incompressible single-phase flow to pass through the porous media in rectangular geometry, but the methodology can be extended to more complex applications, for example two-phase flow using diffuse interface models, which will be pursued in our future researches. The basic principles and proposed methods in this study will be presented next.

2. Principles and Methods

2.1. General Principles

The basic law that is used by experimental and numerical studies for measuring permeability is the Darcy’s law: where is Darcy velocity (m/s), is the permeability tensor (m2), is the gravitational acceleration (m/s2), is the pressure gradient (Pa/m), is the density of fluid (kg/m3), and is the dynamic viscosity of fluid .

The N-S equation (2) and the continuity equation (3) that describe fluid flow in the continuum sense are as follows ( is the Laplace operator): Conventional methods for measuring permeability neglect gravity and drive the fluid flow by pressure because it is easier to control, especially for laboratory experiments. Here, we select to use the full expressions of (1)~(3) to determine the permeability. Periodic boundary condition is adopted so that boundary pressure needs not to be considered. It is pretty easier for numerical implementation and more sensible in nature because subsurface formation is in fact a periodic system or statistically periodic one. The general procedure is stated below.

Step 1. Solve (2) and (3) for incompressible single-phase flow in porous media and obtain several velocity fields as samples.

Step 2. Obtain volumetric velocity of the whole domain according to the samples obtained in Step 1.

Step 3. Solve permeability for the whole domain by directly setting the volumetric velocity as the Darcy velocity in (1).

The detailed numerical methods used in these three steps will be introduced in the next section.

2.2. Numerical Methods

As stated previously, we only consider simple 2D rectangular cases in the present study. The incompressible single-phase flow in porous media is modeled as the flow passing through a square barrier in a square domain with periodic boundary conditions. The numerical model is shown in Figure 1(a). is the domain side length while is the side length of the barrier.  means periodic boundary.

Figure 1: Calculation methodology: (a) numerical model; (b) sampling; and (c) coordinate rotation.

On this kind of domain, the vectors and tensors in (1)~(3) have the following forms: Here , , , , and are the diagonal components of the permeability tensor while and are the off-diagonal components of the permeability tensor. These four variables are all independent of each other. Three typical directions of gravity are chosen to calculate the full-tensor permeability. They are represented by , , and , respectively, and named “Sample a,” and “Sample b,” and “Sample c” as shown in Figure 1(b).

The numerical algorithm for solving (2)~(3) in “Step 1” is projection method with fully explicit spatial discretization (but solve pressure implicitly) using second-order finite difference central type scheme. Cyclic tridiagonal matrix algorithm (CTDMA) is applied to accelerate the convergence of pressure iteration. The mesh size is set as . Flow is considered to reach the steady state as the temporal truncation averaged over the whole domain is lower than .

The volumetric velocity in “Step 2” is defined as follows: where is the volumetric velocity with the components and in the and directions. is the volume of the whole domain. Equation (5) has a discrete expression in this study: where and are the total number of and on grid cells, respectively. The values of and can be determined by solving the N-S equation so that and can be determined.

The volumetric velocity in (5) is actually the Darcy velocity in (1); that is, . In the Darcy scale averaged in the whole domain, pressure gradient should be zero because of the periodic boundary condition; that is, . Thus, (1) can be rewritten as For Sample a , (7) becomes

For Sample b , (7) becomes

For Sample c , (7) becomes

Theoretically, (8a), (8b), (9a), and (9b) are adequate to determine the four components of the full-tensor permeability. However, and are always very close to zeros since they are orthogonal to the bulk flow directions and large numerical errors may exist. Thus, and determined by (8a) and (9b) are always very close to zeros. This apparently violates the fact that and may be far from zeros for some cases. The application of (8a) and (9b) actually makes the full-tensor permeability decay to a diagonal permeability so that they should be canceled. To predict the off-diagonal components accurately, (10a) and (10b) are utilized since the bulk flow in Sample c has equivalent projections in the and directions so that and are both accurate. Therefore, the final equations to determine the four components of the permeability tensor should be the four equations (8b), (9a), (10a), and (10b), which are used in “Step 3.” The procedure is calculating by (9a) and by (8b) and substituting them to (10a) and (10b) to obtain and , respectively.

In practice, it is of interest to detect the maximum permeability, the minimum permeability, principle direction, and anisotropy so that essential features of reservoir, which are independent of coordinate systems and samplings, can be described. Therefore, the original permeability tensor can be transformed by rotating the original coordinate to a new coordinate according to tensors’ characteristics [19]. The angle between the coordinate system and the coordinate system is represented by the angle (Figure 1(c)). Therefore, a series of transformed permeability tensors can be obtained by continuously rotating the coordinate system, that is, continuously changing the value of angle . Once , is achieved at a certain angle ; it represents the case that fluid flows no longer in any tangential directions but only in the two normal directions and so that maximum and minimum components of permeability among all transformed permeability tensors can be examined (, ). In this case, represents the principle direction. To examine the principle direction precisely, ~360° with interval 1° is adopted. Another important parameter is the anisotropy of the permeability in the porous media. Here, we define the ratio to determine the anisotropic degree of permeability: . The tensor transformation formula [19] is shown as follows: where , .

3. Results and Discussion

Two important factors are discussed here: porosity and Reynolds number. To study them separately, two groups of numerical cases are designed. One is changing porosity by changing the size of barrier but keeping constant length of domain. The other one is changing the length of domain but keeping constant porosity. These two factors are discussed in the following two sections. Fluid passing through the porous media is water with density    and dynamic viscosity .

3.1. Discussion on Porosity

The domain size is set as . The size of barrier is set as ~ with interval so that 9 cases named Case 1~Case 9 are generated. According to the size of the computational domain and spatial element sizes on the fixed mesh, time step is set as its maximum value available: .

The flow simulation results by directly solving the Navier-Stokes equation are shown in Figure 2. The contours represent the magnitudes of velocity (i.e., the characteristic speed to be defined later), while the vectors represent the directions of velocity. It can be seen that the flow fields are well represented in the 9 cases. The square barriers are clearly identified by the vectors. We define characteristic speed as , characteristic length as the hydraulic radius , Reynolds number as , and Darcy number as . The characteristic parameters are listed in Table 1 (results for Samples a, b, and c are the same so that the data in Table 1 are any of them). The very low Reynolds number (order of magnitude ~) indicates that the flows are very slow laminar flows. It can be seen that the dimensionless number ReDa1/2 has the order of magnitude ~, which satisfies the suggested restriction of Darcy’s law by Bear and Cheng [20].

Table 1: Characteristic parameters.
Figure 2: Flow fields for different porosities from Navier-Stokes simulation: columns a, b, and c represent Sample a, Sample b, and Sample c; numbers 1~9 represent Case 1~Case 9; the contour shows total speed while the vector shows directions of total velocity.

Porosity is defined as the ratio of void space volume to the total volume as follows: Thus, the 9 cases relate to 9 different values of porosity. Full permeability tensors obtained from flow simulation results of the 9 cases are listed in Table 2. It can be seen that the diagonal components of the permeability tensor all decrease with decreasing porosity. Their values show that the studied porous media fall within the ranges of oil rocks ~, sandstone ~, and good limestone ~ [20]. The diagonal components of transformed permeability tensor are shown in Figure 3 for the 9 porosities; there always exist and , . These mean that the two diagonal components of the full permeability tensor are the same and are always equal to each other along with the rotating coordinate. This indicates that the permeable properties are all the same for any direction so that principle direction does not exist. Therefore, the full-tensor permeabilities are always diagonal tensors with equivalent components for the 9 porosities. The off-diagonal components, which are 3~6 orders of magnitudes lower than the diagonal components (Table 2), should be considered as numerical errors. In this case, the off-diagonal components should be considered as zeros. Revealed from the mathematical meaning, physical meaning is that the porous media in the 9 cases are all isotropic. Corresponding characteristic parameters are listed in Table 3.

Table 2: Full-tensor permeabilities at different porosities.
Table 3: Characteristic permeabilities and anisotropies at different porosities.
Figure 3: Diagonal components of transformed permeability tensor at different porosities.
3.2. Discussion on Reynolds Number

Besides the porosity, Reynolds number is also a key factor which may strongly affect velocity fields so that the results of permeability may be affected. Thus, it is important to check whether the Reynolds number, actually controlled by the domain size, may change the conclusions in Section 3.1. Therefore, the domain size is gradually scaled up from  m by keeping porosity at 0.64 (Case 6 in Tables 2 and 3). The results are shown in Table 4.

Table 4: Full-tensor permeabilities at different Reynolds numbers.

From  m to  m, the diagonal components of the full-tensor permeability just increase two orders of magnitudes once the domain size increases one order of magnitude. The off-diagonal components should still be considered as numerical errors and set to be zeros as stated previously. These three cases are quite similar, and their difference is only the magnitude. It can be verified by the flow fields in Figures 4(a1), 4(b1), 4(c1)~Figures 4(a3), 4(b3), 4(c3), Figures 5(a)~5(c), and Table 5 that the domain size scaling up from  m to  m just leads the velocity to increase proportionally. They are all isotropic. However, the situation is different when the domain size increases up to  m (see line 5 in Table 4). The off-diagonal components are no longer several orders of magnitudes smaller than the diagonal components so that the full-tensor permeability is no longer diagonal. Figures 4(a4)~4(c4) show that the fluid tends to depart the solid walls and the vortex expands. Figure 5(d) shows that the components of the transformed permeability tensor are no longer constant. Maximum and minimum values of the diagonal components can always be obtained at the same time once the off-diagonal components become zeros. Their values (Table 5) show that the porous media at the domain length of 10−3 m enter the range of clean stand (~) which is a better aquifer than the previous ones [20]. Correspondingly, four angles (, , , and ) can be obtained, which are shown in Table 6. Surprisingly from Table 5, the porous media start to show anisotropic property () when the domain length reaches 10−3 m. It can be summarized from Table 6 that is always equal to angle or when   equals . These two angles are on the same line across the second and the fourth quadrants in the original coordinate (see Figure 1(c)). This line is the direction where permeability achieves maximum value, that is, the easiest direction where fluid can pass through, so that this is the principle direction. Accordingly, the direction fluid receives the most resistance (i.e., the direction for ) is the line represented by angles  and , so that this line is the direction orthogonal to the principle direction.

Table 5: Characteristic permeabilities and anisotropies at different Reynolds numbers.
Table 6: Angles at zero off-diagonal components of transformed permeability for  m.
Figure 4: Flow fields for different Reynolds numbers from Navier-Stokes simulation: columns a, b, and c represent Sample a, Sample b, and Sample c; numbers 1~4 represent domain size , , , and , respectively; the contours show total speed while the vectors show directions of total velocity.
Figure 5: Components of transformed permeability tensor at different Reynolds numbers: (a) , (b) , (c) , and (d) .

When the length of domain is greater than or equal to 10−2 m, the flow simulation is divergent. This is probably because the mesh size is not dense enough for larger domain. Larger domains with finer grid, which are very time consuming, will be studied in the future.

To discuss the previous phenomenon that different Reynolds numbers may generate essentially different permeability tensors, or in other words predict different types of porous media for the same configuration, four characteristic parameters are shown in Table 7. and are defined as the ratios of the mean convection effect to the mean diffusion effect in the momentum equation (2): where the superscript “—” represents the average value over the whole domain. It is clear in Table 7 that the fluid flows are actually diffusion dominated for the domain lengths smaller than 10−3 m since and are all much smaller than 1. However, and are greater than 1 at the domain length of 10−3 m. This means that the flow is dominated by the convection effect of fluid. It is well known from their intentions that the diffusion effect tends to transport variables to all directions uniformly while the convection effect tends to transport variables along specific directions. Therefore, the existence of convection dominant flow may be the main reason for the anisotropy. It is worth to point out that Table 7 also verifies the conclusion made by Bear and Cheng [20] that Darcy’s law is usually valid as long as Reynolds number is lower than 1 (occurred at the domain lengths ~ where ~) but sometimes as high as 10 (occurred at the domain length  m where and the suggested restriction is violated).

Table 7: Characteristic parameters at different Reynolds numbers.

4. Conclusion

A numerical method is proposed to compute the permeability in the form of full tensor. The flow simulation results show that flow fields can be well represented by solving the Navier-Stokes equation directly. Original and transformed permeability tensors are obtained so that maximum and minimum components of permeability with principle directions and anisotropies are detected successfully. With this information, the directions of largest and smallest resistances for fluid flow in porous media can be inferred easily.

Through the analyses on the porosity effect and the Reynolds number effect, it is found that porous media with the same porosity and the same configuration can manifest different levels of anisotropy at different Reynolds numbers. At the Reynolds numbers that diffusion dominates flow, isotropy is a good description. At the Reynolds number that convection dominates flow, anisotropy occurs. Thus, it is important to pay attention to Reynolds numbers for porous media applications. This is especially important for applications with large Darcy velocities, such as flow and transport in packed columns and in fractured geological media.

5. Future Issues to Be Addressed

The previous discussions show that the proposed method in the present study is an easy way to determine full-tensor permeability numerically and related characteristics. However, many issues still need to be addressed in future works. Only two of them are briefly listed below as an example.(1)Method for solving the Navier-Stokes equation in porous media should be improved to accelerate the computation. Current speed is not acceptable for engineering applications. Multigrid method needs to be developed for the flow around solids so that the iteration of pressure can be largely accelerated. High-resolution time integration scheme is expected to largely reduce the number of time steps needed. Since only steady-state results are needed for obtaining the full-tensor permeability, a method directly solving the steady-state Navier-Stokes equation is also expected so that a large amount of time integrations can be avoided.(2)Shape, number, position, and so forth of inner barriers should be intensively studied to know the application and limitation of this method and to improve precision. Denser meshes are also expected to study what will happen for larger Reynolds numbers.


The work presented in this paper has been supported in part by the project entitled “Simulation of Subsurface Geochemical Transport and Carbon Sequestration,” funded by the GRP-AEA Program at KAUST. The work has also been supported in part by National Science Foundation of China (no. 51206186, no. 51174206) and Science Foundation of China, University of Petroleum, Beijing (no. 2462012KYJJ0403, no. 2462012KYJJ0404).


  1. C. L. Dinwiddie, F. J. Molz III, and J. W. Castle, “A new small drill hole minipermeameter probe for in situ permeability measurement: fluid mechanics and geometrical factors,” Water Resources Research, vol. 39, no. 7, 2003. View at Publisher · View at Google Scholar · View at Scopus
  2. B. J. Wolfson and S. M. Wentworth, “Complex permeability and permeability measurement using a rectangular waveguide,” Microwave and Optical Technology Letters, vol. 27, no. 3, pp. 180–182, 2000.
  3. A. N. Ivanov, S. N. Kozlova, and A. V. Pechenov, “Permeability measurement,” Measurement Techniques, vol. 43, no. 12, pp. 1086–1088, 2000. View at Scopus
  4. J. L. Jensen and U. Heriot-Watt, “A model for small-scale permeability measurement with applications to reservoir characterization,” in Proceedings of the SPE/DOE 7th Symposium on Enhanced Oil Recovery, pp. 891–900, April 1990. View at Scopus
  5. T. J. Wang, C. H. Wu, and L. J. Lee, “In-plane permeability measurement and analysis in liquid composite molding,” Polymer Composites, vol. 15, no. 4, pp. 278–288, 1994. View at Scopus
  6. P. Ferland, D. Guittard, and F. Trochu, “Concurrent methods for permeability measurement in resin transfer molding,” Polymer Composites, vol. 17, no. 1, pp. 149–158, 1996. View at Scopus
  7. R. Gauvin, F. Trochu, Y. Lemenn, and L. Diallo, “Permeability measurement and flow simulation through fiber reinforcement,” Polymer Composites, vol. 17, no. 1, pp. 34–42, 1996. View at Scopus
  8. Y. Luo, I. Verpoest, K. Hoes, M. Vanheule, H. Sol, and A. Cardon, “Permeability measurement of textile reinforcements with several test fluids,” Composites A, vol. 32, no. 10, pp. 1497–1504, 2001. View at Publisher · View at Google Scholar · View at Scopus
  9. T. F. Fwa, S. A. Tan, C. T. Chuai Y, and K. Guwe, “Expedient permeability measurement for porous pavement surface,” The International Journal of Pavement Engineering, vol. 2, pp. 259–270, 2001.
  10. R. Arbter, J. M. Beraud, C. Binetruy et al., “Experimental determination of the permeability of textiles: a benchmark exercise,” Composites A, vol. 42, no. 9, pp. 1157–1168, 2011. View at Publisher · View at Google Scholar · View at Scopus
  11. J. R. Weitzenböck, R. A. Shenoi, and P. A. Wilson, “Measurement of principal permeability with the channel flow experiment,” Polymer Composites, vol. 20, no. 2, pp. 321–335, 1999. View at Scopus
  12. J. R. Weitzenböck, R. A. Shenoi, and P. A. Wilson, “Radial flow permeability measurement. Part A: theory,” Composites A, vol. 30, no. 6, pp. 781–796, 1999. View at Publisher · View at Google Scholar · View at Scopus
  13. J. R. Weitzenböck, R. A. Shenoi, and P. A. Wilson, “Radial flow permeability measurement. Part B: application,” Composites A, vol. 30, no. 6, pp. 797–813, 1999. View at Publisher · View at Google Scholar · View at Scopus
  14. Y. Keehm, T. Mukerji, and A. Nur, “Computational rock physics at the pore scale: transport properties and diagenesis in realistic pore geometries,” Leading Edge, vol. 20, no. 2, pp. 180–183, 2001. View at Scopus
  15. A. Kameda, J. Dvorkin, Y. Keehm, A. Nur, and W. Bosl, “Permeability-porosity transforms from small sandstone fragments,” Geophysics, vol. 71, no. 1, pp. N11–N19, 2006. View at Publisher · View at Google Scholar · View at Scopus
  16. E. H. Saenger, F. Enzmann, Y. Keehm, and H. Steeb, “Digital rock physics: effect of fluid viscosity on effective elastic properties,” Journal of Applied Geophysics, vol. 74, no. 4, pp. 236–241, 2011. View at Publisher · View at Google Scholar · View at Scopus
  17. B. Quintal, M. Frehner, C. Madonna, N. Tisato, M. Kuteynikova, and E. H. Saenger, “Integrated numerical and laboratory rock physics applied to seismic characterization of reservoir rocks,” Leading Edge, vol. 30, no. 12, pp. 1360–1367, 2011. View at Publisher · View at Google Scholar · View at Scopus
  18. Y. Keehm, T. Mukerji, and A. Nur, “Computational rock physics at the pore scale: transport properties and diagenesis in realistic pore geometries,” Leading Edge, vol. 20, no. 2, pp. 180–183, 2001. View at Scopus
  19. R. L. Bishop and S. I. Goldberg, Tensor Analysis on Manifolds, Dover Publications, 1968.
  20. J. Bear and A. H. D. Cheng, Modeling Groundwater Flow and Contaminant Transport, vol. 23 of Theory and Applications of Transport in Porous Media, Springer Science+Business Media B.V., 2010. View at Publisher · View at Google Scholar