**4. Flow Simulation Examples: Application to Permeability Computation**

#### *4.1. Flow Simulation*

The resulting mesh from the reconstruction process can be used to simulate various physical phenomena, such as those involved in fluid-structure interaction problems. Generally, for composite flow applications, incompressible Stokes flow around the fibers is considered. By considering a stationary regime and neglecting the volume forces, the variational form of the Stokes problem for velocity field, **u**, and pressure field, *p*, is written:

$$\begin{cases} (\mathbf{v}, q) \in \mathcal{V}\_0 \times \mathcal{Q} \\ (2\eta \boldsymbol{\varepsilon}(\mathbf{u}) : \boldsymbol{\varepsilon}(\mathbf{v}))\_{\Omega} - (p\_\prime \nabla \cdot \mathbf{v})\_{\Omega} = 0 \\ (\nabla \cdot \mathbf{u}, q)\_{\Omega} = 0 \end{cases} \tag{4}$$

where *ε* is the strain rate tensor.

A monolithic approach is used, i.e., the flow Equation (4) are solved on the single mesh defined over the whole computational domain, Ω, regardless of the type of phase it contains. The different phases are distinguished by their physical properties which are taken into account through a mixing law. A linear mixture relation is used for the viscosity, *η*, and described by the Equation (5).

$$
\eta = \eta\_f H\_\mathfrak{c} + \eta\_s (1 - H\_\mathfrak{c}) \tag{5}
$$

*η<sup>f</sup>* and *η<sup>s</sup>* are, respectively, the viscosities of the liquid and solid phases. *η<sup>s</sup>* acts as a penalty parameter: when it is high enough, shear rate in the penalized phase becomes close to zero and we find a rigid body motion. This is a simple way to obtain results similar to those provided by an augmented Lagrangian method where a Lagrange multiplier is used to impose a constraint on the solid phase to avoid its deformation [17]. To solve the system (4) using a finite element method, a stabilized approach of VMS type is employed [12]. The used software in this work is ICI-tech, developed at the High Performance Computing Institute (ICI) of Centrale Nantes and implemented for massively parallel context.

#### *4.2. Permeability Computation Procedure*

Predicting permeability is a very important issue in the field of composite forming process. However, it is tricky and complex to obtain experimentally and numerically reliable results, because most simulations are carried out in small periodic representative elementary volumes, under a lot of simplifying assumptions that idealize the real media. Here, we chose to rise to the challenge to numerically determine the permeability tensor of a large virtual sample of fibrous media that imitates sophisticated real media. In threedimensional cases, permeability is characterized by a symmetric second-order tensor **K**. This tensor relates the average fluid velocity **u** to the average pressure gradient on the fluid domain ∇*p<sup>f</sup>* , as shown by the Darcy law below:

$$
\langle \mathbf{u} \rangle = -\frac{\mathbf{K}}{\eta} \langle \nabla p \rangle^f \tag{6}
$$

Using a monolithic approach with finite element discretization, the homogenized velocity and pressure fields are written as the sum of their integration on each mesh element Ω*<sup>e</sup>* of the simulation domain Ω:

$$\langle \mathbf{u} \rangle = \frac{1}{V\_{\Omega}} \sum\_{\mathfrak{c}} \int\_{\Omega\_{\mathfrak{c}}} (1 - H\_{\mathfrak{c}}(\mathfrak{a})) \mathbf{u} \, d\Omega\_{\mathfrak{c}} \tag{7}$$

$$
\langle \nabla p \rangle^f = \frac{1}{V\_{\Omega\_f}} \sum\_{\mathfrak{c}} \int\_{\Omega\_{\mathfrak{c}}} (1 - H\_{\mathfrak{c}}(\mathfrak{a})) \nabla p \, d\Omega\_{\mathfrak{c}} \tag{8}
$$

where *V*<sup>Ω</sup> is the volume of the total domain and *V*Ω*<sup>f</sup>* is the volume of the fluid domain.

To predict permeability, the proposed simulation procedure relies on microstructure generation, phase reconstruction, mesh adaptation, and resolution of the Stokes equations, considering that fibers are static and impermeable. In fact, to determine all components of **K**, three flows in the three directions *x*, *y*, and *z* are successively simulated, an exponent {1, 2, 3} is referred to each one. The flow is induced by an imposed pressure gradient. Depending on the direction where the flow is desired, a constant pressure field on the input face of the simulation domain against a null field on the output face is imposed. For the other faces of the domain, only the normal component of the velocity field is imposed as null. Assuming that the permeability tensor is symmetric and positive definite, its components can be calculated by the resolution of the overdetermined linear system given by:

⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ∇*px*<sup>1</sup> ∇*py*<sup>1</sup> ∇*pz*<sup>1</sup> <sup>000000</sup> <sup>000</sup> ∇*px*<sup>1</sup> ∇*py*<sup>1</sup> ∇*pz*<sup>1</sup> <sup>000</sup> <sup>000000</sup> ∇*px*<sup>1</sup> ∇*py*<sup>1</sup> ∇*pz*<sup>1</sup> ∇*px*<sup>2</sup> ∇*py*<sup>2</sup> ∇*pz*<sup>2</sup> <sup>000000</sup> <sup>000</sup> ∇*px*<sup>2</sup> ∇*py*<sup>2</sup> ∇*pz*<sup>2</sup> <sup>000</sup> <sup>000000</sup> ∇*px*<sup>2</sup> ∇*py*<sup>2</sup> ∇*pz*<sup>2</sup> ∇*px*<sup>3</sup> ∇*py*<sup>3</sup> ∇*pz*<sup>3</sup> <sup>000000</sup> <sup>000</sup> ∇*px*<sup>2</sup> ∇*py*<sup>2</sup> ∇*pz*<sup>3</sup> <sup>000</sup> <sup>000000</sup> ∇*px*<sup>3</sup> ∇*py*<sup>3</sup> ∇*pz*<sup>3</sup> 010 −10 0 0 0 0 001000 −10 0 0000010 −1 0 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ *Kxx Kxy Kxz Kyx Kyy Kyz Kzx Kzy Kzz* ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ = −*η* ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ *ux* 1 *uy* 1 *uz* 1 *ux* 2 *uy* 2 *uz* 2 *ux* 3 *uy* 3 *uz* 3 0 0 0 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ (9)

The solution obtained from the resolution of this matrix system (9) is, obviously, an approximate solution. To ensure a perfect symmetry of **K**, if necessary, the following modification to the extra diagonal terms is made:

$$\mathcal{K}\_{ij}^{final} = \mathcal{K}\_{ji}^{final} = \frac{\mathcal{K}\_{ij} + \mathcal{K}\_{ji}}{2} \tag{10}$$

#### *4.3. Permeability Computation Validation*

To validate permeability computation, the whole procedure was applied to a parallel square packing of fibers having an identical diameter. Rigidity of fibers was ensured by imposing *η<sup>s</sup>* = 500*η<sup>f</sup>* and a zero velocity condition was imposed upon them. Figure 9a shows the used geometry configuration for *Vf* = 25.65%. Equation (11) represents its calculated permeability tensor adimensionalized by the square of fiber radius which respect a transverse isotropic form as expected from the symmetry of the packing.

$$\mathbf{K} = \begin{bmatrix} 0.21 & 0 & 0 \\ 0 & 0.21 & 0 \\ 0 & 0 & 0.28 \end{bmatrix} \tag{11}$$

Permeability evolution according to fiber volume fraction was studied by varying fiber diameter and keeping same the domain size for all simulations. The obtained results of normalized transverse permeability are reported in Figure 9b and compared to the model of [18–20]. The observed permeability values through this graph are in the same order than the one obtained from analytical laws which is relevant to our approach.

**Figure 9.** Comparison of computed permeability with some analytical models: (**a**) Simulated parallel square packing configuration. (**b**) Normalized transverse permeability results.

#### *4.4. Application for 10,000 Fibers*

#### 4.4.1. Microstructure Generation

The first step of the process is the microstructure generation using the octree optimized algorithm described in Section 2. A sample of approximately 10,000 (exactly 10,062) collision-free fibers is created in a cubic domain with an edge length of 1.35 mm. The fibers have a common diameter of 15 μm and a length that follows a normal distribution of mean 0.2 mm and standard deviation 0.03 mm. The obtained volume fraction is **V***<sup>f</sup>* = 14%. The orientation state is nearly isotropic and is given by the following orientation tensor *a*<sup>2</sup> [21]:


Figure 10 shows the set of the generated fibers. Despite the fact that the generation is sequential, these fibers are created in only **1min44s** thanks to the octree contribution.

**Figure 10.** Studied generated microstructure: (**a**) 10,000 generated fibers. (**b**) A random slice showing no collisions.

4.4.2. Microstructure Reconstruction with Adaptative Mesh

The computation was performed on 384 cores. Starting from an initial mesh of ≈4.6 million nodes and ≈27 million elements, after 30 iterations, an adapted final mesh of ≈67 million nodes and ≈391 million elements is created by the methods described in Sections 3.1 and 3.2. For *H<sup>ε</sup>* with *ε* = 3.125 μm, the total immersion and adaptation process required **4h52min** for the 30 iterations. Figure 11 shows the evolution, in a number of elements for each iteration of the mesh adaptation, as well as the computational time. During the first iterations of immersion of the generated fibers in the initial mesh, the mesher adds a considerable number of elements until reaching a peak at the ninth iteration, in order to properly capture the geometries of all the fibers at first. Then, the mesher focuses its work on optimizing the mesh adaptation at the interfaces while respecting a criterion of mesh quality. Once an efficient mesh is achieved, the number of elements stabilizes. The time evolution curve naturally follows the evolution of the mesh size.

**Figure 11.** Evolution of the mesh number of elements (**left** axis) and calculation cost (**right** axis) during the 30 iterations of adaptation of anisotropic mesh, performed on 384 cores

#### 4.4.3. Flow Resolution and Permeability Tensor Computation

Three pressure gradients are applied to the constructed finite element mesh in order to generate the flows required for the identification of **K**. Figure 12 shows the pressure field and velocity vectors around the immersed fibers for the flow in the *x* direction. These results were obtained for a resolution time of the system (4) equal to approximately **7min** minutes on 384 CPUs.

**Figure 12.** Flow according to *x* direction: (**a**) velocity vector around the fibers. (**b**) Zoom around a zone of the figure.

The predicted full permeability tensor adimensionalized by the square of fiber radius for this media is as follows:

> ⎤ ⎦

$$\mathbf{K} = \begin{bmatrix} 0.7322 & -0.0013 & -0.0033 \\ -0.0013 & 0.7444 & -0.0025 \\ -0.0033 & -0.0025 & 0.7089 \end{bmatrix}$$

For isotropic material, only the three diagonal elements are non-null and they are equal. Here, the studied sample is nearly isotropic. For this reason, the obtained diagonal elements are quite similar and the off-diagonal elements are smaller by around two orders of magnitude.

#### **5. Conclusions**

Obtained results show our capability thanks to an octree implementation to deal with big data in terms of input of permeability simulation and to perform reliable finite element calculation on complex geometries. Through the proposed method, further studies can be conducted to better quantify the impact of the microstructural parameters on the permeability and, thus, avoiding problems related to the choice of the size of the simulation domains, which remains rather delicate to define, especially in the case of non-periodic geometries. We can also think about exploring the permeability of multiaxial tissues of the non-crimp fabric (NCF) or textile type. Thanks to the several numerical optimization, the permeability can thus be evaluated at the microscopic scale on several layers by representing the fibers inside the wicks.

**Author Contributions:** Conceptualization and methodology, N.A., L.D., L.S., H.D. and E.A.-C.; software, N.A., L.D., H.D. and P.L.; validation, N.A., L.S., H.D. and E.A.-C.; writing—original draft preparation, N.A. and L.D.; writing—review and editing, N.A., L.S. and E.A.-C. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Acknowledgments:** This work was performed by using HPC resources of the Centrale Nantes Supercomputing Centre on the cluster Liger. Additionally, we would like to warmly thank GENCI and TGCC French National Supercomputing Facility for CPU time under projects A00\*0607575 allowing us to test the approach presented in this article on Joliot-Curie.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

