Computing Effective Permeability of Porous Media with FEM and Micro ‐ CT: An Educational Approach

: Permeability is a parameter that measures the resistance that fluid faces when flowing through a porous medium. Usually, this parameter is determined in routine laboratory tests by applying Darcy’s law. Those tests can be complex and time ‐ demanding, and they do not offer a deep understanding of the material internal microstructure. Currently, with the development of new computational technologies, it is possible to simulate fluid flow experiments in computational labs. Determining permeability with this strategy implies solving a homogenization problem, where the determination of the macro parameter relies on the simulation of a fluid flowing through channels created by connected pores present in the material’s internal microstructure. This is a powerful example of the application of fluid mechanics to solve important industrial problems (e.g., material characterization), in which the students can learn basic concepts of fluid flow while practicing the implementation of computer simulations. In addition, it gives the students a concrete opportunity to work with a problem that associates two different scales. In this work, we present an educational code to compute absolute permeability of heterogeneous materials. The program simulates a Stokes flow in the porous media modeled with periodic boundary conditions using finite elements. Lastly, the permeability of a real sample of sandstone, modeled by microcomputed tomography (micro ‐ CT), is obtained.


Introduction
Permeability is an important macro parameter usually used to characterize porous media. This parameter is a measure of the resistance that fluid faces when flowing through a porous medium. Therefore, this parameter plays an essential role in a great variety of fields, such as soil mechanics, oil prospection, petrophysics, material characterization, etc. The permeability is usually determined in routine laboratory tests by simply applying Darcy's law. However, laboratory tests can be complex, expensive, time-demanding, and sometimes destructive, and they do not offer a deep understanding of the internal microstructure of the material.
Computer simulations are increasingly applied as an alternative to laboratory tests due to their capability to overcome many difficulties presented in physical testing. Some of the advantages of numerical simulations are their non-destructive characteristic, low cost, ease of performing

Effective Permeability through Numerical Homogenization
The determination of permeability through material microstructure analysis is a multiscale problem. The determination of this macro parameter depends on the determination of the velocity field, which is a function of the material's microstructure. This means that, to solve the permeability problem of a medium, it is necessary to solve the governing equations of the microscale. The connection between micro and macroscales must be done through homogenization techniques.
Computational homogenization techniques seek to determine relationships between the micro and the macrostructure of the material, which consist of the determination of the effective properties of representative elementary volumes (REV) by considering characteristics and properties of microscopic elements. For determining the effective permeability of a porous medium, a known pressure gradient or body force is applied to the material REV. Under this condition, a velocity field is induced in the domain. The homogenized effective property of the material is defined as the relationship between the applied pressure gradient 〈 〉 or body force vector 〈 〉 , the fluid viscosity μ, and the average velocities 〈 〉 at the macroscale, known as Darcyʹs law (Ω denotes representative volume).
where 〈 〉 〈 , 〉 is the average velocity vector resulted from the microscopic scale, i.e., ( Darcy's law, given by Equation (1), states that the total velocity (units of displacement by time, e.g., m/s) is equal to the product of the medium permeability (m 2 ), the pressure gradient (unit of pressure per length, e.g., Pa/m), and the inverse of the fluid viscosity (Pa•s).
Considering a two-dimensional problem, the matrix is obtained by two analyses, in which either a unit pressure gradient or a unit body force per direction is applied. The mean velocity vector determined in each analysis is equivalent to that column of matrix related to the respective direction of the pressure gradient (or body force) applied in each case. To determine the velocity field in the problem domain, the governing equations of the problem must be solved through numerical methods, as described in the next section.

Computational Fluid Dynamics with Finite Element Method (FEM)
The continuity equation and the Navier-Stokes equation are the equations that govern the fluid flow problem for most cases of transport in porous medias. As the scale of the problem to be modeled reduces (e.g., nano scales), it is important to be aware that some effects may emerge in the system (e.g., capillarity, chemical interaction). However, this work focuses on viscous fluids at a scale of micrometers; therefore, it can be considered that there is no other effect than those that the Navier-Stokes equation is capable of contemplating. As this work aims to evaluate low-speed incompressible viscous flows only, the convective effects may be neglected (due to low Reynold's number) and, therefore, the Navier-Stokes equation is reduced to the Stokes equation. Although this simplifies the problem, for the great majority of practical applications, these differential equations that govern the problem still have no analytical solution; thus, numerical methods such as the finite element method (FEM) are employed as an essential tool to obtain reliable and approximate solutions.
The FEM proposes to approximate the problem solution by a finite set of elements connected by a finite number of points (commonly called nodes). The collection of finite elements that represents the continuous space is commonly called the discretized domain. Unlike the continuous domain, in the discretized domain, the variables are calculated only in the nodes of the finite elements. The continuous domain of the problem is represented by interpolation functions, called shape functions, which interpolate the calculated variables at the element nodes.
One of the FEM basic requirements is to transform the differential form (strong form) of the governing equations of the problem into integral equations (weak form). This process can be done by applying the weighted residual method, which allows the exact solution of the problem to be approximated by one with null residuals if integrated in the domain of interest. Therefore, the boundary value problem of the continuous medium can be solved by the mean values of a defined range.
The governing equations of the incompressible Stokes flow problem arise from conservation of linear momentum and conservation of mass, which are expressed, in a closed domain Ω, by Equations (3) and (4), respectively.
where μ is the viscosity, , is the velocity vector, p represents the pressure, and , represents the body force. Note that Equation (4) represents the mass continuity equation simplified to incompressible flows, i.e., the density is assumed to be constant (independent of space and time) and the conservation of mass is simplified to the equation of conservation of volume. That assumption is possible because the permeability of porous media is mostly determined in steadystate flow, in which the compressibility effect can be ignored (Akanji and Mattha [4]). The solution of a Stokes flow problem requires the specification of boundary conditions in addition to Equations (3) and (4). As mentioned before, there are basically two ways to compute the average velocity field in order to determine the permeability. It is either by applying a known pressure gradient as boundary conditions or by simply applying known body forces in the domain. The latter approach was chosen because applying domain quantities (such as body forces) instead of applying boundary conditions is the simplest and the most natural choice when using a domainbased discretization procedure such as the FEM to compute the velocity fields.
By applying the weighted residual method, the governing equations are multiplied by continuous and differentiable weight functions  . The sum of residues generated by the approximation is considered null within a defined range, resulting in Using the divergence theorem in the first term of the integral in Equation (5) in combination with the compact support of ψ , i.e., the vanishing boundary integrals, yields Different weight functions can be used in the weighted residual method. One of the most common strategies is to assign to the weight functions the interpolation functions of the nodal values, i.e., the shape functions. This method is known as the Galerkin method. The matrix of shape functions is called N. In classical formulations of the FEM, there arises a matrix whose elements contain derivatives of the shape functions, known as the B matrix. These matrices establish the following relationships: p p = N p, (9)  B LN, (10) where the values and represent the nodal values of u and p, respectively, and L is the matrix of the differential operators.
According to Reddy [4], in order to make an analogy with the virtual powers principle, the weight function vector ψ presented in Equation (5) can be understood as virtual velocities taken by u; therefore, u u  ψ N ψ . The weight function  presented in Equation (6) is associated with the volume change of the fluid element; thus, it should be understood as a virtual hydrostatic pressure that causes this volume variation. That way, we define Applying Equations (8)-(10) and the chosen weight functions to Equations (6) and (7), the weak form of the governing equations can be represented as functions of matrices B and N. Equations (7) and (6) should be respectively rewritten, in the element domain, as follows: As the nodal weight values are arbitrary, the terms in the parentheses of Equations (11) and (12) should vanish, resulting in The terms on the right-hand side of the conservation of momentum equations are the result of the effect of the body forces in each element, called f.
Finally, solving the partial differential Equations (3) and (4) with the application of FEM is simplified to solve the following system of linear algebraic equations: where K is the viscosity matrix, G is the gradient matrix, and D is the velocity divergent matrix and is equal to the transpose of G matrix, The above model is known as mixed formulation. By mixed formulation, it explains that both pressure and velocity are unknown variables. Considering that, usually, fluid dynamic problems have boundary conditions in terms of both velocity and pressure, it becomes the most intuitive formulation to be used. However, the coupling of velocity and pressure variables brings new difficulties to the resolution of the system of algebraic equations.
According to Reddy [6], one way to interpret the coupling between velocity and pressure is to understand the continuity equation as a constraint of incompressibility to momentum conservation equations. Due to the fact that the incompressibility constraint does not depend on pressure variables, constructing a finite element model becomes complicated. The direct solution of the system cannot be performed since the coefficient matrix that multiplies the variable vector has a sub-matrix with a diagonal of null coefficients, which makes the coefficient matrix singular and non-inversible. Thus, it is necessary to use other strategies to solve the equation system.
Another common difficulty of mixed formulation is the choice of shape functions. The relationship between the second derivative of the velocity field and the first derivative of the pressure field stated in the Stokes equation shows that the velocity field must have a greater degree than the pressure field. The use of equal-order approximation functions for both pressure and velocity fields generates instability. This instability can be explained by the solvability analysis of the equation system. The solvability of a linear system is explained by several authors as the guarantee of the infsup condition, also known as the LBB condition, thanks to the works of Ladyzhenskaya [7], Babuska [8], and Brezzi [9]. According to Elman et al., the inf-sup condition is a sufficient condition for the unique solvability of linear systems of equations in the form of Equation (15). Satisfying the inf-sup condition for Equation (15) will guarantee a unique pressure in a closed domain flow.
The use of approximation functions of the same order does not guarantee the inf-sup condition. Therefore, it does not ensure a unique solution for the problem, and the generated results may be non-reliable Hence, finding a system with unique solution requires the use of other methods to approximate the velocity and pressure fields.
Other approximation methods that are recommended for mixed formulation are found in the literature. Elman et al. [10] and Bathe [11]cited the following approximations as satisfactory for the inf-sup condition: Q2-Q1 approximation (bi-quadratic approximation for velocity and bi-linear approximation for pressure), also known as the Taylor-Hood method, Q2-P1 approximation (biquadratic approximation for velocity and linear approximation for pressure in each element without guaranteeing continuity), and Q2-P0 approximation (bi-quadratic approximation for velocity and constant for pressure in each element).
However, using same-order shape functions to interpolate pressure and velocity fields has several advantages from a computational perspective. In order to use same-order shape functions for velocity and pressure fields, and to eliminate instability, stabilization techniques can be adopted. Stabilization techniques have the basic idea of relaxing the restriction of incompressibility, ensuring the inf-sup condition. Many authors proposed different methodologies that proved to be very efficient in stabilizing different approximation models.
One of the simplest discretization methodologies is the use of bi-linear form functions for the velocity and pressure domains. This type of approximation is known as Q1-Q1 approximation. The Q1-Q1 approximation, in addition to being easy to implement, produces lower computational cost for solving the system of equations. For this reason, many authors seek to use stabilization techniques to make use of this type of discretization. As an example, we can cite the works of Andreassen and Andreasen [1], Yang et al. [5], Karim et al. [12], Braack and Schieweck [13], Codina and Blasco [14], and Becker and Hansbo [15].
To construct a stable computational model of Q1-Q1 discretization, this work follows the same stabilization adopted by Andreassen and Andreasen [1]. According to Andreassen and Andreasen [1], the stabilization of the Q1-Q1 model can be obtained by adding a new term with a stabilization coefficient in the continuity equation. The weak form of the continuity equation with the stabilization term is defined by where  is the pressure stabilization parameter, defined as where h is the size of the largest diagonal of the finite element.
The integral with the pressure stabilization term results in a matrix that must be added to the equation system of Equation (15), resulting in where P represents the pressure stabilization matrix.
The introduction of the new term in the continuity equation causes the coefficient matrix to become positively defined. It allows the system to be solved in the direct form.

Implementation of Boundary Conditions
As the volume to be analyzed increases, the property measured for that volume approaches the materialʹs effective property. When the increasing analyzed volume reaches the effective property of the material, the volume can be considered representative. Using periodic boundary conditions, the convergence of the volume's effective properties to the REV's effective properties happens much faster (with smaller volumes). This type of boundary condition simulates the volume of interest within a periodic region, which means that the volume of interest repeats in all directions, as shown in Figure 1. Thus, it is possible to obtain homogeneous and periodic results by modeling smaller volumes, thereby reducing computational effort. The computational implementation of periodic boundary conditions can be performed by assigning the same degrees of freedom to correspondent nodes at the opposite boundaries. As a consequence, the opposite boundaries will have the same behavior, i.e., the velocities and pressure values will be the same at the opposite boundary nodes, enforcing that the volume of interest is indeed a periodic medium. The simulation of flows in a porous medium poses challenges to the problem modeling. The porous media modeling performed in this work is based on fluid modeling only. Ignoring solid elements represented in the image requires imposing the following condition to the solid-fluid interface: the degrees of freedom in the solid-fluid node interface relative to the velocity components are null.

Octave/Matlab Implementation
The code proposed in this paper is designed to receive as input an image (in special micro-CT images) of the problem to be modeled. It means that a pixel-based finite element mesh would be a natural choice of modeling strategy. In this approach, each pixel corresponds to a bi-linear finite element. The benefit of this strategy is that it results in a structured regular mesh (all the elements are squares of the same size) and, therefore, in an extremely low computational cost to generate the mesh. Taking full advantage of the fact that all elements in the mesh are of the same size, we compute the element matrices only once, working solely in the "element domain".
It is important to mention that obtaining accurate values of permeability of real samples through computational simulations relies on contemplating the narrowest constrictions in the porous medium. Covering the material's fine-scale features requires, firstly, images with high resolutions. The resolution has to be high enough to give good representation of the features in the image. On the other hand, the accuracy of the simulation also relies on modeling the image's fine features appropriately. It means that narrowest constrictions such as pore throats and cavities represented by few pixels must be modeled with a finer mesh. Therefore, convergence tests are necessary to compute permeability, taking into account the influence of image's fine-scale features.
The code described in this section is published in Zenodo, an open-source repository. The complete code is available at the following URL: https://doi.org/10.5281/zenodo.3612168.

Input and Initialization
The developed program performs simulation of Stokes flow in binary images. As initial data, the user should enter a binary image (line 03 of the code) representing the REV (as in Figure 1 %% Determining permeability of porous materials %% Reads the image z = imread(ʹfilename.tifʹ); %% Increases the image resolution by n times n = 1; z = repelem(z,n,n); %% Displays the image figure(1) imshow(z) %% Determines the effective permeability lx = 1; % horizontal dimension ly = 1; % vertical dimension KH = compute_permeability(lx,ly,z); The MatLab function "repelem" (code line 06) is used to refine the mesh by a factor n, dividing each image pixel in a group of n × n square finite elements. This is an important procedure to perform convergence tests.
Note that, in the initialization step, the matrix image is transformed into a vector that indicates the nature of each element in the mesh (whether it is solid or fluid).
It is convenient that input data (lines 01 to 13) be stored in their own file, thus decoupling the main program (which calculates the permeability) from the data inputs (model: image and dimensions). This way, we can save the data of several different problems. Therefore, the main program is defined in its own function, which is presented from line 14 of the code.
To determine permeability, we call the function "compute_permeability", which takes as parameters the dimensions of the model and the information of each element whether it is solid or fluid (Listing 2): Listing 2. Implementation of the main function to compute permeability. 14  It is also necessary to create a matrix containing the element connectivity, which indicates the nodal incidences of each finite element. The connectivity matrix is assembled so that the image is inserted in a periodic medium. For achieving that, the nodes of the opposite edges are numbered the same. Figure 2 illustrates how the image nodes are numbered using periodic boundary conditions. To implement this strategy of mesh node numbering (see lines 21 to 24) and generation of a connectivity matrix (see lines 25 to 29), which considers periodic boundary conditions, Listing 3 is used.

Coeficient Matrix and RHS Contribution of Each Element
Since the program proposes performing simulations from images, and the mesh generation uses a pixel-based method, all elements are equal; thus, the viscosity, gradient, and pressure stabilization matrices can be calculated only once. This considerably reduces the simulation time. These matrices are incorporated into the main program using Listing 4.   (8) The previous code, although presented in a didactic way, does not make explicitly state the required elements for the calculation of each matrix independently. However, since the shape functions used to interpolate the velocity and pressure fields are the same, the vector f and the matrices k, g, and pe can be calculated within the same function, making the algorithm more compact and elegant.

Boundary Conditions and Linear Equation System Solution
The implementation of zero velocity condition on the fluid-solid interface is done by reducing the coefficient matrix to solve only non-null variables. This process is performed on lines 113 through 123. The system is solved on line 123 (Listing 8).

Obtaining Permeability
After assembling the vector with a nodal variable, permeability is calculated relating the average velocities to the prescribed pressure gradient. In this step, the mean velocities of each element are calculated integrating the velocity field inside each element. The velocity field of each element is achieved using the nodal velocities, and the shape function matrix N is used to approximate the velocity field. This process can be performed using Listing 9. In the above algorithm, the element average velocity is calculated in each element during the for loop. The permeability is determined using the micro and macro relationships stated in Equations (1) and (2). However, for a more elegant and efficient way, taking advantage of the already created matrix of degrees of freedom, we use Listing 10.

Validation Test
In order to validate the program presented in this work, a Stokes flow simulation was performed in the porous medium represented in Figure 3. The permeability obtained by the code demonstrated in this work, which was based on the code of Andreassen and Andreasen [1], was compared with the analytical solution proposed by Drummond and Tahir [16]. According to Drummond and Tahir [16], the permeability of the medium shown in Figure 3 can be determined using Equation (22).
where permeability (k) is a function of the grain's radius (r) and the volume fraction of solid (c). In Reference [14], the proposed analytical equation is valid for small values of c. Therefore, for program validation, the REV of the medium was modeled with unit dimensions and grains with radius equal to 0.1. For these adopted values, the corresponding c can be registered, and the permeability k = 0.0281 × 10 −3 is obtained via the analytical solution. Figure 4 shows the results obtained through numerical modeling using the program proposed in this paper. As Figure 4 shows, the convergence of the permeability of the medium displayed in Figure 3 can be obtained from mesh size 400 × 400. From this mesh size, the difference between the numerical result and analytical result is 0.71%. This means that the program gives good results and fulfills its purpose.

Image Acquisition: Microcomputed Tomography
The convergence of numerical results to experimental data has direct dependence on the fidelity of the representation of the computational model. One way to create reliable computational models is by using X-ray microcomputed tomography (micro-CT). Micro-CT is a technique used for internal three-dimensional visualization of any material on the micrometer scale. This technique consists of generating multiple radiographic projections taken from various angles to produce a 3D image of an object's internal structure. The generation of a micro-CT involves several steps: sample preparation and assembly, projection acquisition, reconstruction, image processing and segmentation, mesh generation, and computer simulations, as illustrated in Figure 5.  In the micro-CT acquisition process, the sample to be digitalized is placed on a turntable inside the micro-CT scanner that rotates around its own axis while radiographies are taken. The radiographic projections are used by a reconstruction algorithm to create a 3D image of the object. The reconstruction is done by an inverse analysis that determines the level of X-ray absorption of a spatial point inside the material. This inverse analysis is made based on the different absorption levels on each radiographic projection. The result of a micro-CT is a 3D digital image of an object. A digital image is a visual representation of an object formed by small elements called voxels. Each voxel is given a specific value that is responsible for attributing it a color. Figure 6 illustrates the result of a micro-CT performed by the authors of a sandstone sample.

Image Processing and Segmentation
After scanning and reconstruction, a 3D virtual object is generated. The visualization of the materialʹs internal microstructure can be done by cutting planes, as illustrated in Figure 6. However, for quantitative analyzes, such as porosity, permeability, and other effective properties, images must be processed. The process of image processing includes selection of region of interest (ROI), filter application for noise reduction, and image segmentation.
To determine sandstone permeability from 2D models, a slice (2D image) that was considered to be representative was chosen. After that, image filters were used. Filters are tools that reduce interference and artefacts, and that enhance certain features of images. They are based on sophisticated algorithms that facilitate the extraction of specific information from an image by modifying the imageʹs gray histogram. As a last step, the image was segmented. Since the image obtained by micro-CT is a grayscale image, it is important to determine which regions of the material will receive their respective properties (e.g., which voxels will be defined as solid or pores). Segmentation consists of determining the material phases according to the gray tones of the micro-CT voxels. This process can be performed by using the image histogram. By the use of the image histogram, one can define the range of gray tones that represent the pores and the range of gray tones that represent the material. There are different strategies that can be used in the image segmentation process. This process can take a lot of time, and it demands the user's experience to produce good results. After basic phase separation of the sandstone virtual sample, a tool for grain separation was used to enhance the quality of the segmentation. Using the grain separation tool allows for a richer analysis of fluid flow through the grains, which is essential for permeability determination. The grain separation tool used is based on a distance map function that determine the boundaries of the grains considering the 3D geometry of pores and grains. Since this tool is performed on 3D images, the result of the grain separation in a 2D image slice can present connections to pores (as presented in Figure 7) that seems to artificially enforce pore connection. The use of this tool was necessary to ensure that the fluid simulated in this model would be able to flow inside the material. It is possible to make a comparison between before and after image processing and segmentation of the sandstone in Figure 7. In Figure 7, blue represents the solid region and black represents the porous connected region where the fluid flows. Note that different users and different tools would eventually generate different segmented images. It is important to reinforce that the segmentation of the 2D slice of sandstone is used here as an example of input for the program to compute absolute permeability. Different results of segmentation of the sandstone sample can be used as input for the program as any other 2D image.

Sandstone Absolute Permeability
After validation, the program was used to measure the permeability of a sandstone sample represented by a 2D image. The two-dimensional image of the sandstone to be analyzed is exactly the same image presented in Figure 6, which has dimensions of 1 mm × 1 mm. Since the 2D image of the sandstone sample is much smaller than the sample size tested experimentally in the laboratory, periodic boundary conditions must be applied to the sandstone virtual model; thus, no change in the proposed program is necessary. Table 1 presents the obtained values for different mesh sizes (number of elements to represent the image). Figure 8 shows the sandstone permeability convergence as the finite element mesh is refined. The permeability obtained for the 2D image of the sandstone sample analyzed with the most refined mesh was Kxx = 0.1019 μm 2 and Kyy = 0.1603 μm 2 . This means that the sample is anisotropic with respect to permeability. The mesh size of 2400 × 2400 was limited by the random-access memory (RAM) of 32 Gb. A computer with 16 Gb RAM would allow refinement up to 1600 × 1600. Note that the purpose of this paper was to present the methodology to compute permeability through computational tool in an educational way. The computation of permeability of the sample of sandstone here was used as an example of the application of the methodology. In order to obtain more accurate results for this sample through convergence tests with more refined meshes, we suggest the use of more powerful computers or the implementation of the presented methodology in higher-performance languages, such as C.

Conclusions
This paper presented an educational program capable of computing the permeability of porous materials by simulating fluid flows governed by the Stokes equation. The necessary theoretical basis for modeling the problem was explained, and the peculiarities of the FEM applied to computational fluid dynamics were discussed. Different ways for assembling the system of equations to be solved were explained, as one of the formulations capable of reducing significantly the processing time.
Because it is part of the digital petrophysical field, the program proposed here uses binary images as data input. This feature allows the permeability of real porous materials to be determined by means of micro-CT, which is capable of revealing the internal microstructure of the material through different voxels to represent both the fluid and the solid phases. In this work, the porous medium was modeled in a way such that only the fluid domain was considered. This methodology is extremely efficient for materials that have low porosities, as it allows the permeability to be obtained using less computational effort compared to solid and fluid domain modeling strategies.
The program described by the authors was validated through a simple problem that has a known analytical solution. Hereafter, the effective permeability of a sandstone sample was computed by a 2D analysis. Although the validation of the code confers credibility to the proposed computational program, the sandstone permeability value calculated in two dimensions may not be representative for a three-dimensional sandstone sample due to the three-dimensional geometric characteristics of its pores. However, the development of the 2D program presented here gives the reader the necessary knowledge to extend the program for 3D cases in a simple and intuitive way.