A Novel Finite Element Based Method for Predicting Permeability of Heterogeneous and Anisotropic Porous Microstructures

1. Introduction

Permeability is an average medium property that measures the ability of the porous medium to transport the fluid. Permeability is defined as a scalar value for isotropic porous media, but can be defined as a tensor to reflect its dependence on the orientation of the porous structure relative to the coordinate system.

In science and engineering (e.g. rocks, soils, synthetic and natural porous materials, etc...), permeability is a very important physical property of porous media because it controls the directional movement and flow rate of fluids in the porous structure. Permeability is influenced by factors such as the pore structure, through its size distribution and tortuosity, and the type of permeating fluid.

The applications of permeability are diverse and include various fields such as civil engineering, environmental engineering, petroleum engineering, geotechnical engineering, tissue engineering, agricultural engineering and materials engineering to name a few. For example, permeability is required for the design and construction of various civil engineering projects such as underground structures, dams and reservoirs. It helps engineers to predict how soil behaves when subjected to various loads and water pressures. Permeability is simple in concept, but has some very complex aspects in practice, especially when trying to make accurate measurements or predictions.

Manufacturing of composite materials, whether for semi-finished products such as prepregs or finished parts, relies on the impregnation by the resin of a dense network of fibers. The latter can be viewed as a porous medium. Controlling the transfer of a liquid resin within this medium is essential to obtain defect-free materials and parts. Knowledge of the permeability of the fiber network is generally required to adjust manufacturing parameters such as pressure, flow rate, and temperature.

The slender form factor of the fibers and the preferred orientations they assume in many engineered parts result in a fairly pronounced directionality of the flows. In addition, optimizing part design for the highest mechanical performance-to-mass ratio means varying the orientations and even the volume fraction of the fibers in the part. The fiber network becomes heterogeneous and anisotropic throughout the part. Measuring the components of the permeability tensor becomes difficult when characterizing fibrous preforms where the arrangement of the fibers changes spatially, such as materials produced by automated fiber placement. Due to the high cost and long duration of such experiments, it is attractive to determine the permeability by numerical simulation.

Compared to granular porous media whose permeability is often isotropic, the permeability of fibrous networks is almost never isotropic. Anisotropy is also observed in fractured porous rocks, where a dual-scale flow occur in the network of interconnected fractures and in the porous rock itself. The permeability of such a medium results from the interaction of flows within the two interconnected networks. In this paper, we focus on porous media with heterogeneous but monodisperse porosity, which can lead to anisotropic permeabilities. The single phase flow in such material is investigated here.

The general approach of computing the permeability of a porous material is to first solve the stationary Stokes equations to obtain the pressure and/or velocity fields in the computational domain, and then to average the computed fields by an appropriate permeability identification technique. Different numerical methodologies for the solution of this problem with different discretisation techniques, physical variables formulation, boundary conditions etc. exist in the literature. A recent benchmark on image-based permeability determination of fibrous materials has illustrated this diversity of numerical approaches [1]. The most widely used numerical approximation methods are the Finite Volume Method (FVM) and the Finite Element Method (FEM). Other methods, such as the Finite Difference Method (FDM), Pore Networks Modeling (PNM), the Lattice Boltzmann Method, or the Smoothed Particle Hydrodynamics (SPH) method, can also be used for flow field calculation [2]. However, as it was shown in the benchmark study [1], the choice of numerical discretisation method itself is not a dominant factor influencing the permeability prediction. The application of appropriate boundary conditions and the used permeability identification technique are much more important. Due to the heterogeneity of numerous materials, permeability calculations should be performed on a representative volume, which should be large enough to account for all geometric characteristics of the medium. Challenges then arise regarding the efficiency of numerical problem solving and the amount of memory required.

2. Objectives and Content of the Study

The main objectives of this work are:

• To propose and investigate a finite element based numerical methodology that can predict the homogenized permeability of heterogeneous and anisotropic porous microstructures where the flow in the microstructure is induced by a body force.
• To demonstrate the necessity to calculate the full permeability tensor including its off-diagonal components in order to capture the influence of the local anisotropy.

3. Materials and Methods

3.1. Modeling of Single-Phase Steady-State Porous Media Flow at Micro Scale

In this work, a fully saturated porous medium occupying the space $\Omega \in {\mathbb{R}}^n$ is considered to be made up of up to two non-overlapping phases, one representing a rigid, non-permeable solid (${\Omega }_s$) and the other representing an incompressible linear viscous fluid (${\Omega }_f$). The steady-state problem of flow [3] through such porous media decribed by the Stokes equation can be formulated as follows,

$$\begin{array}{cccc}\hfill \nabla ·\mathit{\sigma }& =\mathit{b}\hfill & \hfill \mathrm{in}{\Omega }_f& (\mathrm{Balance}\mathrm{of}\mathrm{linear}\mathrm{momentum})\hfill \end{array}$$

(1)

$$\begin{array}{cccc}\hfill \nabla ·\mathit{v}& =0\hfill & \hfill \mathrm{in}{\Omega }_f& (\mathrm{Incompressibility}\mathrm{constraint})\hfill \end{array}$$

(2)

$$\begin{array}{cccc}\hfill \mathit{\sigma }& =-p\mathit{I}+2\mu {\nabla }^s\mathit{v}\hfill & \hfill \mathrm{in}{\Omega }_f& (\mathrm{Constitutive}\mathrm{model}\mathrm{for}\mathrm{linear}\mathrm{viscous}\mathrm{fluid})\hfill \end{array}$$

(3)

$$\begin{array}{cccc}\hfill \mathit{v}& =0\hfill & \hfill \mathrm{in}{\Omega }_s& (\mathrm{Impermeable}\mathrm{solid})\hfill \end{array}$$

(4)

$$\begin{array}{cccc}\hfill \mathit{v}& ={\mathit{v}}_D\hfill & \hfill \mathrm{in}{\Gamma }_f^D& (\mathrm{Dirichlet}\mathrm{boundary}\mathrm{conditions}\mathrm{for}\mathrm{fluid})\hfill \end{array}$$

(5)

$$\begin{array}{cccc}\hfill \mathit{n}·\mathit{\sigma }& =\mathit{t}\hfill & \hfill \mathrm{in}{\Gamma }_f^N& (\mathrm{Neumann}\mathrm{boundary}\mathrm{conditions}\mathrm{for}\mathrm{fluid})\hfill \end{array}$$

(6)

$$\begin{array}{cccc}\hfill {\mathit{v}}_D& =0\hfill & \hfill \mathrm{in}{\Gamma }_s& (\mathrm{No}\text{-}\mathrm{slip}\mathrm{boundary}\mathrm{conditions}\mathrm{on}\mathrm{solid}\mathrm{interfaces})\hfill \end{array}$$

(7)

where $\mu$ is the dynamic viscosity of the fluid, p is the mechanical definition of the pressure, $\mathit{t}$ is the imposed boundary traction and ${\nabla }^s\mathit{v}$ is the strain-rate defined by the symmetric part of the gradient of the velocity as,

$$\begin{array}{c}\hfill {\nabla }^s\mathit{v}={\displaystyle \frac{1}{2}}\left[\nabla \mathit{v}+{(\nabla \mathit{v})}^T\right]\end{array}$$

(8)

In the literature, several studies have demonstrated the use of a penalty based approach [4,5]. It involves a relaxation of the incompressibility condition with a pseudo-compressibility condition as follows,

$$\begin{array}{cccc}\hfill \nabla ·\mathit{v}& =-{\displaystyle \frac{p}{\lambda }}\hfill & \hfill \mathrm{in}{\Omega }_f& (\mathrm{Pseudo}\text{-}\mathrm{compressibility}\mathrm{constraint})\hfill \end{array}$$

(9)

where $\lambda$ is a penalty parameter. The convergence of the numerical solution using a penalty based approach has been demonstrated in [6]. The benefit of using the penalty based approach is that it uncouples pressure and velocity. With this, one can solve the momentum equations in terms of only the velocity field without considering additional degrees of freedom (DoFs) associated with the pressure.

Despite the limitations and challenges associated with the choice of value for the penalty parameter itself as discussed in [4]; the reduction in the total number of degrees of freedom to solve is attractive for computationally costly problems such as the one in the Virtual permeability benchmark [1] discussed further in Section 4.1. Thus, with a computational cost in mind, a modeling choice was made to use the penalty based approach in this work.

3.2. PoroS 1.0: An Image-based Stokes Flow Solver

A software suite named `PoroS’ consisting of new voxel- (3D) or pixel- (2D) based porous media permeability solvers has been developed that uses finite element approach. The numerical architecture of the solver has been specifically designed so that:

• It can directly read the segmented images of the microstructure obtained from the micro-CT scan and perform the flow and permeability computations on them.
• It can also directly read a voxelized geometry of a digital twin of a material generated using the widely used open-source software TexGen^® [7,8].
• It can handle flow problems with a large number of degrees of freedom (order of billions of DoFs).

Figure 1. Overview and scope of PoroS: The input can be either a segmented 3D scan of the material or a digital twin of the material generated using TexGen^®. The scope of PoroS consists of the Stokes flow problem solver and a full-field homogenization subroutine. Thus the output of the solver consists of flow fields of the flow problems and a full 3D permeability tensor

Figure 1. Overview and scope of PoroS: The input can be either a segmented 3D scan of the material or a digital twin of the material generated using TexGen^®. The scope of PoroS consists of the Stokes flow problem solver and a full-field homogenization subroutine. Thus the output of the solver consists of flow fields of the flow problems and a full 3D permeability tensor

Alhough PoroS contains several memory efficient solvers to solve the linear system of equations, the results presented in this work use in particular one of the solvers which is a global-matrix free preconditioned conjugate gradient solver.

3.2.1. Weak Form

Following the finite element approach with penalty; the weak form for the problem described in Section 3.1 reads:

$$\begin{array}{c}\hfill \underset{\mathrm{Viscous}\mathrm{Term}}{\underbrace{\underset{\Omega }{\int }\mu (\nabla \mathit{w}:\nabla \mathit{v})\mathrm{d}\Omega }}+\underset{\mathrm{Penalty}\mathrm{Term}}{\underbrace{\underset{\Omega }{\int }\lambda (\nabla ·\mathit{w})(\nabla ·\mathit{v})\mathrm{d}\Omega }}=\underset{\mathrm{Body}\mathrm{Force}\mathrm{Term}}{\underbrace{\underset{\Omega }{\int }\mathit{w}·\mathit{b}\mathrm{d}\Omega }}+\underset{\mathrm{Neumann}\mathrm{Boundary}\mathrm{Term}}{\underbrace{\underset{{\Gamma }_f^N}{\int }\mathit{w}·\mathit{t}\mathrm{d}\Gamma }}\end{array}$$

(10)

where $\mathit{w}$ is the test function for velocity such that,

$$\begin{array}{c}\hfill \mathit{w}\in \mathcal{V}:={\mathcal{H}}_{\mathrm{d}{\Omega }_f^D}^1(\Omega )=\left\{\mathit{w}\in {\mathcal{H}}^1(\Omega )|\mathit{w}=\mathbf{0}\mathrm{on}\mathrm{d}{\Omega }_f^D\right\}\end{array}$$

(11)

3.2.2. Discretization

In PoroS, a voxel-based approach (pixel-based in 2D) is followed. For 3D, a 27-node Taylor-Hood hexahedral element (HEX27) is implemented for the resolution of the velocity field. Similarly for 2D, a 9-node quadrilateral element (QUAD9) is implemented [3,4,9]. From here onwards, the derivations are presented in 3D formulation; however, equivalent development is done for 2D formulation as well.

The velocity vector at an arbitrary point $\mathit{q}\in {\Omega }_f$ is then given by,

$$\begin{array}{c}\hfill \left[\begin{array}{c}\underset{(3\times 1)}{{\mathit{V}}_{\mathit{q}}}\end{array}\right]=\sum _{i=1}^{27}\left[\begin{array}{ccc}{N}_i& 0& 0\\ 0& N_i& 0\\ 0& 0& N_i\end{array}\right]\left[\begin{array}{c}{V}_i^x\\ V_i^y\\ V_i^z\end{array}\right]=\left[\begin{array}{c}\underset{(3\times 81)}{\mathit{N}}\end{array}\right]\left[\begin{array}{c}\underset{(81\times 1)}{{\mathit{V}}_{\mathit{e}}}\end{array}\right]\end{array}$$

(12)

where $N_i$ are the classical shape functions corresponding to the HEX27 element evaluated at point $\mathit{q}$, and ${\mathit{V}}_{\mathit{e}}$ is the nodal velocity vector of the element.

3.2.3. Linear System of Equations

Using 12, the Equation (10) can be transformed into a linear system of equations as follows,

$$\begin{array}{c}\hfill \left[{\mathit{A}}^{\mathit{v}}+{\mathit{A}}^{\mathit{\lambda }}\right]\left[\mathit{V}\right]=\left[\mathit{F}\right]\end{array}$$

(13)

where ${\mathit{A}}^{\mathit{v}}$ is the global viscous stiffness matrix and ${\mathit{A}}^{\mathit{\lambda }}$ is the global penalty stiffness matrix assembled from the elemental stiffness matrices ${\mathit{A}}_{\mathit{e}}^{\mathit{v}}$ and ${\mathit{A}}_{\mathit{e}}^{\mathit{\lambda }}$, respectively. They are defined as,

$$\begin{array}{cc}\hfill {\mathit{A}}_{\mathit{e}}^{\mathit{v}}& =\underset{{\Omega }_e^f}{\int }{{\mathit{B}}_{\mathit{v}}}^T\mathit{D}{\mathit{B}}_{\mathit{v}}\mathrm{d}{\Omega }_e^f\hfill \end{array}$$

(14)

$$\begin{array}{cc}\hfill {\mathit{A}}_{\mathit{e}}^{\mathit{\lambda }}& =\underset{{\Omega }_e^f}{\int }{{\mathit{B}}_{\mathit{\lambda }}}^T\lambda {\mathit{B}}_{\mathit{\lambda }}\mathrm{d}{\Omega }_e^f\hfill \end{array}$$

(15)

In the latter, the strain-rate ($\mathit{B}$) matrices are defined as,

$$\begin{array}{cccc}\hfill \underset{(6,81)}{{\mathit{B}}_{\mathit{v}}}& =\left[\begin{array}{cccc}\underset{(6,3)}{{\mathit{B}}_{\mathit{v}}^{\mathbf{1}}}& \underset{(6,3)}{{\mathit{B}}_{\mathit{v}}^{\mathbf{2}}}& ....& \underset{(6,3)}{{\mathit{B}}_{\mathit{v}}^{\mathbf{27}}}\end{array}\right]\hfill & \hfill \underset{(6,81)}{{\mathit{B}}_{\mathit{\lambda }}}& =\left[\begin{array}{cccc}\underset{(6,3)}{{\mathit{B}}_{\mathit{\lambda }}^{\mathbf{1}}}& \underset{(6,3)}{{\mathit{B}}_{\mathit{\lambda }}^{\mathbf{2}}}& ....& \underset{(6,3)}{{\mathit{B}}_{\mathit{\lambda }}^{\mathbf{27}}}\end{array}\right]\hfill \end{array}$$

(16)

$$\begin{array}{cccc}\hfill \underset{(6,3)}{{\mathit{B}}_{\mathit{v}}^{\mathit{i}}}& =\left[\begin{array}{ccc}{\displaystyle \frac{\partial N_i}{\partial x}}& 0& 0\\ 0& {\displaystyle \frac{\partial N_i}{\partial y}}& 0\\ 0& 0& {\displaystyle \frac{\partial N_i}{\partial z}}\\ {\displaystyle \frac{\partial N_i}{\partial y}}& {\displaystyle \frac{\partial N_i}{\partial x}}& 0\\ 0& {\displaystyle \frac{\partial N_i}{\partial z}}& {\displaystyle \frac{\partial N_i}{\partial y}}\\ {\displaystyle \frac{\partial N_i}{\partial z}}& 0& {\displaystyle \frac{\partial N_i}{\partial x}}\end{array}\right]\hfill & \hfill \underset{(6,3)}{{\mathit{B}}_{\mathit{\lambda }}^{\mathit{i}}}& =\left[\begin{array}{ccc}{\displaystyle \frac{\partial N_i}{\partial x}}& {\displaystyle \frac{\partial N_i}{\partial y}}& {\displaystyle \frac{\partial N_i}{\partial z}}\\ {\displaystyle \frac{\partial N_i}{\partial x}}& {\displaystyle \frac{\partial N_i}{\partial y}}& {\displaystyle \frac{\partial N_i}{\partial z}}\\ {\displaystyle \frac{\partial N_i}{\partial x}}& {\displaystyle \frac{\partial N_i}{\partial y}}& {\displaystyle \frac{\partial N_i}{\partial z}}\\ 0& 0& 0\\ 0& 0& 0\\ 0& 0& 0\end{array}\right]\hfill \end{array}$$

(17)

and the viscous 3D constitutive matrix $\mathit{D}$ is a diagonal matrix given by,

$$\begin{array}{cc}\hfill \underset{(6,6)}{\mathit{D}}& =\mu \mathrm{diag}(2,2,2,1,1,1)\hfill \end{array}$$

(18)

In 13 the term $\left[\mathit{V}\right]$ corresponds to the global velocity nodal vector, whereas the term $\left[\mathit{F}\right]$ corresponds to the external force, which in the scope of this work is applied in the form of a body force similar to gravity. It is given by,

$$\begin{array}{c}\hfill \underset{(81,1)}{\mathit{F}}=\underset{{\Omega }_e^f}{\int }{\left[\begin{array}{c}\underset{(3\times 81)}{\mathit{N}}\end{array}\right]}^T\left[\begin{array}{c}\underset{(3\times 1)}{\mathit{b}}\end{array}\right]\mathrm{d}{\Omega }_e^f\end{array}$$

(19)

where $\left[\begin{array}{c}\mathit{b}\end{array}\right]={\left[\begin{array}{ccc}{b}_x& b_y& b_z\end{array}\right]}^T$ is a vector denoting the body force.

3.2.4. Permeability Determination Procedure

The permeability tensor of a saturated porous medium ($\mathit{K}$) is associated with Darcy’s law, which relates the pressure gradient ($\nabla P$) to the volume-averaged fluid velocity ($\mathit{v}$) by,

$$\begin{array}{c}\hfill \mathit{v}=-\frac{\mathit{K}}{\mu }\nabla P\end{array}$$

(20)

$\mathit{K}$ is a second order symmetric, positive-definite tensor, it is an intrinsic property of the porous medium. The emergence of the Darcy equation from the Stokes equation is well justified by multi-scale asymptotic homogenization [10] or averaging techniques [11]. In full-field numerical methods, the permeability tensor is obtained by solving a series of incompressible single-phase flow problems in a domain of the porous medium using a numerical method. Typically, for a three-dimensional problem, one has to solve 3 independent flow problems with an averaged flow imposed by the boundary conditions in the way to ensure the independence of these 3 solutions. A first approach consists then in extracting the averaged velocity and pressure gradients. Then one has to solve the linear system with 9 equations and 9 unknowns,

$$\begin{array}{c}\hfill \begin{array}{c}\left[\begin{array}{c}{K}_{xx}\\ K_{xy}\\ K_{xz}\\ K_{yy}\\ K_{yz}\\ K_{zx}\\ K_{zy}\\ K_{zz}\end{array}\right]=\mu {\left[\begin{array}{ccccccccc}\nabla P_x^1& \nabla P_y^1& \nabla P_z^1& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& \nabla p_x^1& \nabla p_y^1& \nabla p_z^1& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& \nabla p_x^1& \nabla p_y^1& \nabla p_z^1\\ \nabla P_x^2& \nabla P_y^2& \nabla P_z^2& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& \nabla p_x^2& \nabla p_y^2& \nabla p_z^2& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& \nabla p_x^2& \nabla p_y^2& \nabla p_z^2\\ \nabla P_x^3& \nabla P_y^3& \nabla P_z^3& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& \nabla p_x^3& \nabla p_y^3& \nabla p_z^3& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& \nabla p_x^3& \nabla p_y^3& \nabla p_z^3\end{array}\right]}^{-1}\end{array}\left[\begin{array}{c}{v}_x^1\\ v_y^1\\ v_z^1\\ v_x^2\\ v_y^2\\ v_z^2\\ v_x^3\\ v_y^3\\ v_z^3\end{array}\right]\end{array}$$

(21)

where the subscripts x, y, z denote directions, and the subscripts 1, 2, 3 denote a numerical simulation run in each direction. The solution obtained by solving this matrix system is an approximate solution. For this reason, the symmetry of the permeability tensor identified at this stage is not always respected. In most cases, the calculated $K_{ij}$ components are not identical to the $K_{ji}$ components. In such cases, to ensure the symmetry, it is common to make the following modification to the off-diagonal terms:

$$\begin{array}{c}\hfill K_{ij}^{final}=K_{ji}^{final}=\frac{K_{ij}+K_{ji}}{2}\end{array}$$

(22)

The degree of asymmetry is also induced by the type of boundary conditions used to perform the calculation. The asymmetry reduces as the sample size increases because boundary effects contribute much less.

Another approach relies on the $\mathit{homogenization}$ to connect the pore-scale flow and the continuum scale flow [12]. Similar to the homogenization in micromechanics, where the Hill-Mandel lemma is used as an equivalence criterion between physical quantities at both scales, here the total power dissipated by the flow at both scales is used. When the flow occurs through porous media, there is a dissipation of energy equal to the hydraulic pressure loss along the flow field. This energy loss is equal to the work of the viscous forces resisting the flow.

The powers at micro- and macro-scales are given by,

$$\begin{array}{c}\hfill P_{micro}=\underset{RV}{\int }\mathit{\sigma }:\nabla \mathrm{v}\mathrm{d}\mathit{x}\end{array}$$

(23)

$$\begin{array}{c}\hfill P_{macro}=\nabla P·\mathit{V}\left|w\left(\mathit{X}\right)\right|\end{array}$$

(24)

$\left|w\right(\mathit{X}\left)\right|$ is called as the measure of $w\left(\mathit{X}\right)$. Following the derivation in [12] one obtains,

$$\begin{array}{cc}\hfill 2\mu \langle \mathit{D}:\mathit{D}\rangle & =\nabla P·\mathit{V}\hfill \end{array}$$

(25)

$$\begin{array}{cc}\hfill 2\mu \langle \mathit{D}:\mathit{D}\rangle & ={\mathit{V}}^T·{\mathit{K}}^{-1}·\mathit{V}\hfill \end{array}$$

(26)

where $\mathit{D}$ is the rate of strain tensor. Thus, knowing the velocity fields (and subsequently the viscous dissipation) from the flow solutions, one can obtain the full permeability tensor from Equation (26). This method does not require to calculate the pressure field in order to compute the permeability.

3.2.5. Boundary Conditions

Solving the equations of the heterogeneous pore flow model and those of the upscaled macroscopic flow model requires in both cases the definition of a set of boundary conditions (BC).

A commonly used condition is to reproduce those used in the laboratory to measure the permeability of a material. In the case of composite fiber reinforcements, the domain is a rectangular parallelepiped (rectangular cuboid). A constant pressure is applied to one face and a different pressure is applied to the opposite face, creating a constant pressure gradient, the other surfaces are impermeable (i.e. zero pressure-gradient normal to their plane, no-slip velocity). This boundary condition is sometimes called permeameter boundary condition. An alternative is to impose a constant macroscale velocity on the two permeable boundaries. The main issue with this BC is that it suppresses important transverse flows and can result in unphysical flows. This problem is overcome by using periodic boundary conditions. This have been used widely in the framework of homogenization theory. However, this condition rarely corresponds to the situation encountered with fibrous reinforcement stacks, where the bounding box of the RVE intersects the solid and fluid phases non-periodically.

$\mathit{Symmetry}$ is an intermediate condition between the permeameter and periodic BCs. The velocities tangential to the bounding walls are allowed to evolve as if there were no walls, but all velocities normal to the walls are zero. This BC also suppresses transversal flow within the material, but partially circumvents the no-flow boundary problem along all four faces of the computational domain in terms of diagonal permeability tensor terms by allowing non-zero velocities along the boundaries. Periodization of non-periodic porous media can be achieved by the translation or symmetry procedures [13], the first one, however, creating a risk to provide a non-percolating geometry, while the second one increasing multiple times the size of a computational domain.

It has been observed that near the boundary of the computational domain, flow characteristics may better reflect prescribed physical flow constraints than the effects of local permeability. This can bias the expected directionality of fluid flow and therefore the predicted permeability. To mitigate this effect, it has been proposed to add a boundary region to impose boundary conditions further away from the boundaries by immersing the domain of interest in a larger domain with the same type of heterogeneity [13,14]. An equivalent approach is to average the calculated fields over a restricted sub-volume of the flow away from the boundaries. This approach of $\mathit{sub}\text{-}\mathit{volume}\mathit{immersion}$ involves additional computational cost due to the presence of the border region.

3.2.6. Body Force Driven Flow

A common element in minimizing the bias introduced by the non-representativeness of the processed digital samples and boundary conditions is to address the largest possible porous media in 3D, avoiding a periodization strategy. The computational efficiency of the numerical technique and the ability to apply appropriate boundary conditions are central in this challenge. Similar to LBM, where a uniform body force can be used instead of a pressure gradient to produce the same flow rate as pressure-driven flow [15], body force driven flow is used in PoroS. It has been shown that the permeability tensor asymmetry problem for very heterogeneous samples is reduced by forcing a flow with a body force into the sample and applying periodic boundary conditions [16].

In this work for the calculation of permeability, a uniform unit body force is applied in an appropriate direction depending on the flow problem under consideration, i.e. for flow problem in X direction, $b_x=1,b_y=0,b_z=0$ is imposed, similarly, for flow problem in Y and Z directions, $b_x=0,b_y=1,b_z=0$ and $b_x=0,b_y=0,b_z=1$ are imposed, respectively.

3.3. Validation of PoroS

This section discusses three test cases from the literature that serve as the validation of the PoroS solver for the Stokes flow solution. The first test case consists of a Poiseuille flow in a 3D pipe with Dirichlet-type boundary condition without considering the body force (Section 3.3.1), whereas the second and third test cases (Section 3.3.2 and Section 3.3.3) involve the calculation of the permeability using body force approach.

3.3.1. Poiseuille Flow in a 3D Pipe with Circular Cross Section

In order to validate the Stokes flow solver, a voxelized version of the classical test case of Poiseuille flow in a 3D pipe of a circular cross section was considered. Even though the cross sectional geometry was considered to be rectangular, $\mathit{v}=0$ was imposed on all nodes on the voxelized boundary of the cross section and all elements outside the boundary were set to solid regions thereby restricting the flow outside the circular pipe (Figure 2a). The pipe has a length of 5 m and a diameter of 1 m. 20 voxels were used across the diameter for discretization. A uniform velocity (along X) of 1 m/s was applied at the circular inlet. After performing the simulation, the fully developed $V_x$ profile at the center of the cross section (Figure 2b) was taken for comparison with the analytical solution discussed in [17]. A good match was observed between the numerical solution calculated using PoroS and the analytical solution (Figure 2c).

3.3.2. Validation of Prediction of Transverse Permeability of a Fiber Array

To validate the numerical prediction of permeability by PoroS of a microstructure in the form of array of aligned cylinders (fibers) transversely to their direction, a test from [18] has been considered. It involves a flow through a RVE of this microstructure with a quadratic packing, i.e. around a cylindrical obstacle depicted in Figure 3a. A voxelized geometry of $101\times 101\times 11$ [$\mu .m$] was considered with a resolution of 1 [$\mu .m$/voxel]. Several cases with different fibre volume fractions ($V_f$) were considered by varying the diameter of the cylindrical obstacle. Simulations were carried out considering (i) Dirichlet condition driven flow and (ii) body force driven flow without additional fluid buffer zone, both with a convergence criteria of 0.1% for the permeability. The $K_{xx}$ component of the permeability tensor for all the cases is plotted in Figure 3b that shows a good match with the analytical solution discussed in [18].

3.3.3. Validation of Prediction of Longitudinal Permeability

To validate the numerical prediction by PoroS of the longitudinal permeability, a test inspired from [19] has been considered. It involves a flow in the domain containing a channel with a constant square cross section, and the volume fraction of this channel within the domain is being varied (Figure 4). A geometry of dimensions $100\times 100\times 100$ [$\mu .m$] was considered with a resolution of 1 [$\mu .m$/voxel]. Four sub-cases of this test were considered where a different diameter of the inner square channel was considered: 20, 40, 60 and 80 [$\mu .m$] as shown in Figure 4.

The problem was solved for a fluid of viscosity 1 [$Pa.s$], and a penalty number 1000 was used for imposing the incompressibility constraint. The simulations were carried out considering the body-force driven loading. The numerical results obtained using PoroS were compared to the empirical results presented in [20], where the equivalent longitudinal permeability of the domain with channels with a square cross section is given by:

$$\begin{array}{c}\hfill K_{\mathrm{empirical}}=0.035144{\displaystyle \frac{l^4}{A}}\end{array}$$

(27)

where l is the channel cross section edge length and A is the full cross section area including the solid region.

It is evident from the comparison (Table 1), that the numerical prediction of the permeability is in very good agreement with the corresponding empirical values.

The validation using these three conventional test cases served for testing not only the solver but also the new finite element based approach using volumetric forcing condition to predict the permeability. This first step of validation is essential before proceeding to compare the permeability prediction results on a much bigger RVE from the international virtual permeability benchmark. This is discussed in the next section.

4. Results and Discussion

This section consists of two parts. In Section 4.1, the results computed by PoroS were compared to the results from the international virtual permeability benchmark on fibrous microstructures. In Section 4.2, the importance of a correct computation of a full permeability tensor including the off-diagonal terms is discussed illustrating it, firstly, with a simple example of a 2D inclined channel network, and, secondly, with a real 3D microstructure.

4.1. Comparison of Permeability Values Obtained by PoroS with the Results from the International Virtual Permeability Benchmark

The first international virtual permeability benchmark conducted for the engineering textiles [1] has discussed and compared various numerical modeling approaches in detail. The benchmark has provided 50 values computed by 16 participants for the same input geometry (a 3D segmented image of the representative volume element), which can be found at the repository [21]. It was obtained from the 3D X-ray microscope scan of a twill-weave glass fiber reinforced composite (HexForcce 01102 from HEXCEL) (Figure 5a). The input was a 3D voxelized image with dimensions of 1003x124x973 voxels and a scan resolution of 0.521 [$\mu m/voxel$] in each direction. This problem with a total of about 120 million voxels resulted in a total of about 3 billion DoFs.

In the following sections, three comparative studies were undertaken where the results were compared for: the full geometry, the geometry divided into 10 sub-volumes across the fiber direction, and finally the geometry divided into 10 sub-volumes along the fiber direction.

4.1.1. Comparison with the Results on the Full Geometry of the RVE

The results obtained with PoroS 1.0 using the body forcing fall within the cluster of the results of the participants of the benchmark (Figure 5b). Additionally, the results are very close to the mean values of the benchmark results which was obtained after eliminating the outliers (refer to [1] for more details on how the means were obtained)

This indicates that the approach of using body-force driven flow in a microstructure provides good results for the permeability of a heterogeneous microstructure, which are consistent with the results obtained using different modeling techniques reported by various participants in the benchmark.

4.1.2. Comparison with the Results on the Sub-Volumes: Transverse Cuts

Further analysis was conducted to ensure the validity of this approach by comparing the permeability results computed on 10 sub-volumes of the RVE of the benchmark, where the RVE was subdivided by slicing the geometry transverse to the fiber direction as shown in Figure 6a. With such a subdivision a monotonic evolution of the fiber volume fraction along the fiber direction can be observed within this sample: it increases from 0.54 to 0.59 along the fiber direction. This evolution is due to the woven geometry of the fabric where warp and weft tows intersect and thus can be locally more compacted at these crossings. The natural decrease in permeability with increasing $V_f$ in this case was not always correctly captured by the benchmark participants, who made this study on sub-volumes [1]. For the longitudinal component of the permeability $K_{yy}$ in the direction of the fibers, i.e. when the flow has straighter and less complex paths through the geometry, this trend was correctly captured by all calculations; PoroS also showed a correct trend (Figure 6c). However, predicting consistent values of $K_{xx}$ with respect to $V_f$ was challenging for the benchmark participants (Figure 6b), while PoroS succeeded in predicting the decreasing permeability trend. The only exception to the decreasing trend was a higher $K_{xx}$ value calculated on sub-volume 3 with $V_f=0.547$, which was also detected by Software-Q and Software-R (Figure 6b). This value can be explained by the fact that with a $V_f$ close to the previously calculated $V_f=0.544$, this subvolume had a locally very open channel with the dominant effect on permeability.

It was noted in [1] that a correct dependence of the Software-R result on the fiber volume fraction was obtained thanks to the permeability identification procedure, which allows to compute the full permeability tensor, which is also the case for PoroS (this point will be discussed in more detail in Section 4.2). Quantitatively, however, the $K_{xx}$ result of Software-R is too far away from the cluster of consistent benchmark results, while the values predicted by PoroS correlate well with the cluster. The results of Software-T also yielded the full tensor, but were classified as non-consistent in the benchmark due to the numerical methodology used.

For the transverse permeability component $K_{zz}$, more participants were able to predict a correct decreasing trend (Figure 6d), because numerically the flow had much simpler paths through the small thickness of the domain in this direction compared to the case of $K_{xx}$. PoroS also predicted very close to other participants’ values.

4.1.3. Comparison with the Results on the Sub-Volumes: Longitudinal Cuts

The main motivation for the participants to conduct the study on the sub-volumes of the RVE of the benchmark was the significant computational cost associated with handling the whole RVE. From the computations on the transverse sections considered in the previous section, the participants derived the equivalent transverse permeability components of the entire RVE. Similarly, they used longitudinal slices of the RVE (parallel to the fiber direction: Figure 7a) to calculate the equivalent axial permeability component $K_{yy}$ considered in this section.

On the contrary, as shown previously, PoroS was not limited by memory capacity in predicting the permeability of the entire RVE. Its study on the cross sectional volumes was performed to prove its ability to correctly capture the permeability dependence on the fiber volume fraction compared to the approaches used in the benchmark. Contrary to the direction along the fibers with their gradual compaction and thus monotonic change of $V_f$, the longitudinal sub-volumes did not have such trend of $V_f$ change since they were distributed along the width of the tow. However, locally the fiber volume fraction still varied from sub-volume to sub-volume (Figure 7a) due to the misalignment and twisting of the fibers within the tow. The results of the predicted permeability $K_{yy}$ on 10 sub-volumes are plotted against increasing $V_f$ in Figure 7b. The general trend of decreasing permeability is mainly observed in the results of the indicated software and in the results of PoroS. It should be noted that, unlike the transverse slices of the previous section, the longitudinal slices each contain a very small number of fibers that are not statistically representative of the entire microstructure. This can explain the fluctuations in the decreasing trend of the permeability, especially for close values of $V_f$. Moreover, all results in Figure 7b follow the same trend and are close to each other. The exception is the $K_{yy}$ value predicted by PoroS at $V_f=0.516$. This small $V_f$ corresponds to the boundary subvolume 10, which has a peculiar geometry with some fibers strongly misaligned with respect to the axial Y direction. The determination of the full permeability tensor by homogenization in PoroS allowed to capture these effects. This point is further discussed in Section 4.2.2.

The benchmark [1] demonstrated the importance of boundary conditions applied tangentially to the flow direction. For this tow microstructure, it was shown that the application of the symmetric BC resulted in higher axial permeability $K_{yy}$ values than the periodic BC. The same is confirmed for the sub-volume study: the results of Software-A also used symmetric BC in Figure 7b. This tendency can be explained by the geometry of this microstructure due to some channels at the boundaries, which are external to the whole RVE, and which double their size in the case of symmetry application, thus opening large flow paths in the axial direction. These channels do not increase in size when periodicity is applied. Since PoroS is a body force driven flow, as opposed to conventional boundary conditions, the majority of the values in Figure 7b computed by PoroS lie between the values of periodic BC and symmetric BC. This fact underlines another advantage of the method implemented in PoroS: in domains with strong influence of boundary conditions, the body force approach solves the problem of their adequate choice.

In conclusion, the observations discussed in these comparative sections indicate that the approach of using body-force driven flow through a microstructure results in a good prediction of its permeability, and this was demonstrated starting with academic test cases from the literature and then further by comparisons with the results from the virtual permeability benchmark. This completes the first objective of this article. The second objective, which is to demonstrate the need to calculate the full permeability tensor including its off-diagonal components, is discussed next.

4.2. Importance of Predicting the Full Permeability Tensor

To demonstrate the importance of predicting the full permeability tensor of the structure, two test cases are presented. Section 4.2.1 discusses a representative test case that has been specially constructed. Section 4.2.2 revisits the longitudinal subchannel 10 to discuss the peculiarities of the flow field computed within it.

4.2.1. A Case with a 2D Inclined Channel Network

This representative test case consists of a 2D orthogonal channel network where wider channels are inclined at an angle of ${60}^{\circ }$ with respect to X axis (Figure 8a). PoroS uses a pixelized image in 2D ($500\times 285$ pixels in X and Y directions for this test case with the resolution of 1 [$\mu m/px$]) as an input where one pixel corresponds to one element as shown in Figure 8b.

The steady-state Stokes flow problems in X and Y directions were solved using PoroS solver with the body force approach as described in Section 3.2. The field of magnitude of velocity for X and Y flow problems is shown in Figure 9.

After performing the homogenization of the structure, the macroscopic permeability tensor was calculated to be,

$$\begin{array}{c}\hfill \left[\begin{array}{cc}{K}_{11}& K_{12}\\ K_{21}& K_{22}\end{array}\right]=\left[\begin{array}{cc}4.41\times {10}^{-13}& 6.66\times {10}^{-13}\\ 6.66\times {10}^{-13}& 1.20\times {10}^{-12}\end{array}\right]\end{array}$$

(28)

It is evident from the relative magnitudes of the off-diagonal terms, that they cannot be neglected for such a structure. Going further, the principal flow directions for this case can be obtained by calculating eigenvalues (${\lambda }_i$) and eigenvectors ($V_i$) of the permeability tensor in Equation 28 which are,

$$\begin{array}{c}\hfill {\lambda }_1=5.45e-14,V_1=\left[\begin{array}{c}-0.87\\ 0.50\end{array}\right]=149.{91}^{\circ };{\lambda }_2=1.59e-12,V_2=\left[\begin{array}{c}0.50\\ 0.87\end{array}\right]=59.{91}^{\circ }\end{array}$$

(29)

Angles calculated in Equation 29 provide the validation that the main flow direction (corresponding to the highest eigenvalue i.e. ${\lambda }_2$) is along $59.{91}^{\circ }$ which matches with the orientation of the wider channels (${60}^{\circ }$).

4.2.2. The Importance of Identifying the Dominant Flow Direction

The calculation of the skew terms of the permeability tensor becomes essential for the materials with the architecture for which the flow principal directions are unknown. For instance, the 3D woven fabric studied in [22] exhibited an asymmetric flow in the thickness direction due to its complex weaving pattern. This effect was clearly reflected on the non-negligible order of magnitude of the measured off-diagonal components of its permeability tensor.

For a typical quasi-isotropic or orthotropic material, it might be sufficient to just evaluate the diagonal components of the permeability. However, the permeability of the anisotropic materials will be determined by the principal flow paths and their directions. This principal flow path will often be driven by the orientations of the material itself (as demonstrated in Section 4.2.1) but also it depends on the dominant flow direction which may or may not be the same as that of the material orientation. For example, an interesting flow pattern was observed while considering the longitudinally cut subvolume 10 from Section 4.1.3 while looking at the flow problem Z in the axial direction of fibers. Here even though the majority of the fibers are indeed oriented along the Z direction, the dominant flow direction is different from Z direction (indicated by the dotted line in Figure 10). This is due to a strong misalignment of some fibers.

It was observed that the ratio $K_{zz}/K_{xz}$ was high for most of the longitudinally cut subvolumes of the benchmark (going as high as 1600 for subvolume 6). Whereas, in the case of subvolume 10, of it was found to be to 44. Thus, even though one can know the overall material direction of a microstructure in consideration, it is not obvious that the dominant flow direction would be along this direction of the material. This emphasizes the need for evaluating the full permeability tensor including its off-diagonal components. These effects are visible on a large scale during the manufacture of composite parts, but not necessarily on too small digital samples. It is therefore important to ensure the representativeness of the latter in order to predict these directional effects and the resulting cross-flows.

5. Conclusions

In this work, a novel finite element based approach for predicting the permeability of heterogeneous and anisotropic porous microstructures was proposed where the flow in the microstructure is induced by a body force. In order to investigate the applicability of this approach, a new solver suite called PoroS 1.0 was developed that consists of Stokes equation flow solvers and the permeability calculation program. It has image based FEM solvers that work directly on the segmented images of the microstructure obtained from either the micro-CT scan of the material or using the digital twin of the material generated using the widely used TexGen^® software. The architecture of PoroS has been specifically designed to work with a large number of degrees of freedom (order of billions), which is often the case for such voxelized models generated from the scans.

PoroS was first validated with various standard flow and permeability problems from the literature and then a detailed comparison was made with respect to the results of the international numerical permeability benchmark, thereby validating both the solver and the novel approach. Finally, the importance of predicting the full permeability tensor was demonstrated especially in the case of heterogeneous and anisotropic porous microstructures where the dominant flow direction might be different from the dominant material structure orientation.

As a perspective, a version 2.0 of PoroS has been under development that will address the dual-scale flow problems using the resolution of Stokes-Brinkman equations while also taking into account the local orientations of the porous microstructure. It will retain all advantages of PoroS 1.0, while introducing new features in order to handle more advanced multi-scale problems in numerical permeability estimation.