1. Introduction

2. Objectives and Content of the Study

3. Materials and Methods

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.

4. Results and Discussion

5. Conclusions

Abbreviations

References

1. Syerko, E.; Schmidt, T.; May, D.; Binetruy, C.; Advani, S.; Lomov, S.; Silva, L.; Abaimov, S.; Aissa, N.; Akhatov, I.; others. Benchmark exercise on image-based permeability determination of engineering textiles: Microscale predictions. Composites Part A: Applied Science and Manufacturing 2023, 167, 107397. [CrossRef]
2. Wagner, A.; Eggenweiler, E.; Weinhardt, F.; Trivedi, Z.; Krach, D.; Lohrmann, C.; Jain, K.; Karadimitriou, N.; Bringedal, C.; Voland, P.; others. Permeability estimation of regular porous structures: A benchmark for comparison of methods. Transport in porous media 2021, 138, 1–23. [CrossRef]
3. Donea, J.; Huerta, A. Finite element methods for flow problems; John Wiley & Sons, 2003.
4. Hughes, T.J.; Liu, W.K.; Brooks, A. Finite element analysis of incompressible viscous flows by the penalty function formulation. Journal of computational physics 1979, 30, 1–60. [CrossRef]
5. Heinrich, J.; Marshall, R. Viscous incompressible flow by a penalty function finite element method. Computers & Fluids 1981, 9, 73–83.
6. Temam, R. Navier-Stokes equations: theory and numerical analysis; Vol. 343, American Mathematical Soc., 2001.
7. Long, A.; Brown, L. Modelling the geometry of textile reinforcements for composites: TexGen. In Composite reinforcements for optimum performance; Elsevier, 2011; pp. 239–264.
8. Brown, L.; mike matveev.; georgespackman. louisepb/TexGen: TexGen v3.13.1, 2023. [CrossRef]
9. Taylor, C.; Hood, P. A numerical solution of the Navier-Stokes equations using the finite element technique. Computers & Fluids 1973, 1, 73–100.
10. Sánchez-Palencia, E. Non-homogeneous media and vibration theory. Lecture Note in Physics, Springer-Verlag 1980, 320, 57–65.
11. Whitaker, S. Flow in porous media I: A theoretical derivation of Darcy’s law. Transport in porous media 1986, 1, 3–25. [CrossRef]
12. Lopez, E.; Abisset-Chavanne, E.; Lebel, F.; Upadhyay, R.; Comas, S.; Binetruy, C.; Chinesta, F. Flow modeling of linear and nonlinear fluids in two and three scale fibrous fabrics. International Journal of Material Forming 2016, 9, 215–27. [CrossRef]
13. Guibert, R.; Horgue, P.; Debenest, G.; Quintard, M. A Comparison of Various Methods for the Numerical Evaluation of Porous Media Permeability Tensors from Pore-Scale Geometry. Mathematical Geosciences 2016, 48, 329–347. [CrossRef]
14. Wen, X.; Durlofsky, L.; Edwards, M. Use of border regions for improved permeability upscaling. Mathematical Geology 2003, 35, 521–547. [CrossRef]
15. Nabovati, A.; Llewellin, E.W.; Sousa, A.C. A general model for the permeability of fibrous porous media based on fluid flow simulations using the lattice Boltzmann method. Composites Part A: Applied Science and Manufacturing 2009, 40, 860–869. [CrossRef]
16. Zakirov, T.; Khramchenkov, M. Study of the pore space heterogeneity effect on the absolute permeability tensors calculated under different boundary conditions and driving forces using a “Computational Rock Physics” technology. Journal of Petroleum Science and Engineering 2022, 216, 110750. [CrossRef]
17. White, F.M. Fluid mechanics; Tata McGraw-Hill Education, 1979.
18. Gebart, B.R. Permeability of unidirectional reinforcements for RTM. Journal of composite materials 1992, 26, 1100–1133. [CrossRef]
19. Saxena, N.; Hofmann, R.; Alpak, F.O.; Berg, S.; Dietderich, J.; Agarwal, U.; Tandon, K.; Hunter, S.; Freeman, J.; Wilson, O.B. References and benchmarks for pore-scale flow simulated using micro-CT images of porous media and digital rocks. Advances in Water Resources 2017, 109, 211–235. [CrossRef]
20. Mavko, G.; Mukerji, T.; Dvorkin, J. The rock physics handbook; Cambridge university press, 2020.
21. Syerko, E. International Virtual Permeability Benchmark 3D Image Dataset of the Fiber Tow Microscopic Sample, 2022. [CrossRef]
22. Yun, M.; Sas, H.; Simacek, P.; Advani, S.G. Characterization of 3D fabric permeability with skew terms. Composites Part A: Applied Science and Manufacturing 2017, 97, 51–59. [CrossRef]