3D robust iterative coupling of Stokes, Darcy and solid mechanics for low permeability media undergoing finite strains

Arnaud Dereims, Sylvain Drapier, Jean-Michel Bergheau, Patrick de Luca

► To cite this version:

Arnaud Dereims, Sylvain Drapier, Jean-Michel Bergheau, Patrick de Luca. 3D robust iterative coupling of Stokes, Darcy and solid mechanics for low permeability media undergoing finite strains. Finite Elements in Analysis & Design, 2015, 94, pp.1-15. 10.1016/j.finel.2014.09.003. emse-01145419

HAL Id: emse-01145419

https://hal-emse.ccsd.cnrs.fr/emse-01145419

Submitted on 10 Sep 2021

**HAL** is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# 3D robust iterative coupling of Stokes, Darcy and solid mechanics for low permeability media undergoing finite strains

A. Dereims<sup>a,*</sup>, S. Drapier<sup>b</sup>, J.-M. Bergheau<sup>c</sup>, P. de Luca<sup>a</sup>

Stokes, Darcy and solid mechanics coupling is a matter of interest in many domains of engineering such as soil mechanics, bio-mechanics, and composites. The aim of this paper is to present a robust iterative method to deal with this coupling for low permeability media within the framework of industrial simulation, and especially for composite manufacturing processes. Stokes and Darcy problems are solved using a mixed velocity–pressure finite element using a mini-element formulation, and coupled together by the so-called Beavers–Joseph–Saffman conditions through their interface. This fluid formulation is then coupled to a non-linear solid mechanics formulation in finite deformations using Terzaghi's law at the pore level, and an explicit dependence of permeability with respect to porosity that is exactly computed from the solid mechanics kinematics. Then, those formulations are validated with test-cases and by the Method of the Manufactured Exact Solution (MMES) (Knupp and Salari, 2003[1]). Finally, a 3D curved transient example of application is presented.

## 1. Introduction

This paper presents an iterative method to strongly couple Stokes and Darcy flows, in low permeability media, with finite strain solid mechanics in order to deal with fluid seepage into and out from low permeability media undergoing finite strains. It is applied here to complex manufacturing processes for composites in an industrial framework. It represents an improved extension of the work presented by Celle in [2,3] to 3D and complex shapes with low permeability down to $10^{-15} \text{ m}^2$, opening the possibility to deal with industrial cases. In this approach two main problems have to be tackled, first the coupling between Stokes (pure fluid) and Darcy (flow through low permeability porous medium) flows and second, the fluid/structure interaction. In the case of low permeability media, the couplings are quite tricky to handle numerically especially due to the system conditioning where, for example, penalty coefficients will combine with very low contributions. Special attention will be paid, in this paper, to validate the numerical development in this framework.

Stokes and Darcy coupling has recently been a matter of interest for many researchers in various domains of application such as bio-mechanics [4], environment (river flow) [5,6] or composite processing

[2,3,7]. From those works, one can identify two strategies to deal with Stokes–Darcy coupling, a unified approach [4,7,8] relying on a monolithic formulation of the so-called Brinkman problem (Stokes + Darcy), and an iterative one [3,9] based on two formulations interacting through boundary/interface conditions.

As for the fluid structure interaction, it is usually considered at a scale of observation where the porous medium can be seen as homogeneous. Then, the corresponding solid–fluid mechanical formulation is seen as ‘smeared-out’ and corresponds to poroelasticity which has been widely used and therefore developed in soil mechanics [10–12] and bio-mechanics [13,4]. In this case as well, different ways are presented in the literature to deal with fluid–structure interaction. First, simplified fluid–solid mechanics frameworks can be considered, introducing in 1D models some relationships relating explicitly interstitial fluid pressure with thickness variation [14–18], therefore limiting the domain of interest to simplistic processes. Second, in a more general framework, some researchers chose to solve a fully three dimensional solid–fluid mechanical formulation either through heavy monolithic approaches [19] or using an iterative method based on permeability and stress updates depending on displacement–velocity considerations [2,3]. In any case, all these coupling methods between fluid flow and solid mechanics within porous media rely on Terzaghi's principle [10].

The originality of the proposed method is that it remains robust and accurate even for real severe physical parameters such as low

<footnote><sup>a</sup> ESI Group, Aero Business Center - Aeroparc, 25 Rue Marcel Issartier, BP 20005, 33700 Mérignac, France</footnote>

<footnote><sup>b</sup> École Nationale Supérieure des Mines de Saint-Étienne, Structure and Material Science Division, Georges Friedel Laboratory (UMR CNRS 5307), 158 cours Fauriel, 42023 Saint-Étienne Cedex 2, France</footnote>

<footnote><sup>c</sup> University of Lyon, École Nationale d'Ingénieur de Saint-Étienne, LTDS (UMR CNRS 5513), 58 rue Jean Parrot, 42023 Saint-Étienne Cedex 2, France</footnote>

<footnote>* Corresponding author. Tel.: +33 535540510.
E-mail address: Arnaud.Dereims@esi-group.com (A. Dereims).</footnote>
permeability values, thin pure fluid layer, large deformations and complex geometries which are often discarded in recent publications but definitively more relevant with respect to industrial applications. Second, the proposed industry-targeted method feature stands in its global approach which yields well-founded mechanistic relationship between flows and deformations. All the phenomena (pure and porous flows, and solid mechanics) are fully simulated and strongly coupled iteratively together. This permits to minimize the number of physical parameters which are required for describing the coupling mechanisms, somehow tricky and expensive to characterize within industrial context. Indeed, this approach does not require to perform compression test on wet fiber beds, only the constitutive law of the dry preforms is required. Moreover, the use of three independently formulated formulations makes it easy to use one, two of them, or the whole package to simulate all kinds of processes based on the same physical phenomena and to add, in the future, other formulations such as thermo-physico-chemical aspects, to simulate more complex coupled phenomena and processes.

This paper will present first the FE formulation for the fluid problem, in Stokes, Darcy, and coupled regimes, in Section 2. This coupled flow formulation will be then validated using test cases along with the Method of the Manufactured Exact Solution (MMES) [1] in Section 3. Second, the finite strain solid mechanics formulation will be exposed (Section 4) and coupled with the fluid flow problem (Section 5) and validated. Eventually, examples of fluid infiltration into 3D compacting orthotropic porous media will be given.

## 2. Formulation of the fluid mechanics problems

Let us formulate first the stabilized Stokes, Darcy, and Stokes-Darcy coupled problems. As depicted in Fig. 1, for the fluid flow one may split the whole domain $\Omega$ to be studied, into two areas.

The purely fluid region $\Omega_s$ is governed by the Stokes equations system:

$$ \begin{cases} \nabla \cdot (2\eta D(\underline{v}_s)) - \nabla p_s = 0 \\ \nabla \cdot \underline{v}_s = 0 \end{cases} \quad (1) $$

where $\eta$ is the fluid dynamic viscosity, $v_s$ the fluid velocity, $p$ the fluid pressure and $D(v) = \frac{1}{2}(\nabla v + \nabla^T v)$ the Eulerian strain rate associated with the velocity field $\underline{v}$. The porous medium area $\Omega_D$

is governed by the Darcy equations system:

$$ \begin{cases} \underline{v}_D = \frac{1}{\eta} K (-\nabla p_D + \rho g) \\ \nabla \cdot \underline{v}_D = 0 \end{cases} \quad (2) $$

where $\underline{v}_D$ is the Darcy velocity or the macroscopic mean velocity of the fluid with respect to the preforms, $p$ the pore pressure, $K$ the permeability tensor assumed diagonal in the element frame, $\eta$ the fluid viscosity and $\rho$ the volumetric mass of the fluid. The two domains are separated by an interface $\Gamma = \Omega_s \cap \Omega_D$ through both mass conservation and normal stress continuity must be fulfilled (see Fig. 1). In addition, it is chosen to control the tangential velocity jump with a drag coefficient $\alpha$ leading to the Beavers-Joseph-Saffman [20] condition. Since Darcy's law assumptions imply that shear is neglected ($\sigma_D = -p_D I$), the normal stress continuity is equivalent to a Dirichlet boundary condition on the hydrostatic pressure when using the primal formulation of Darcy system (5), leading to the following system of conditions on $\Gamma$:

$$ \begin{cases} \underline{v}_S \cdot \underline{n} = \underline{v}_D \cdot \underline{n} \\ 2 \underline{n} \cdot D(\underline{v}_S) \cdot \underline{\tau}_i = \frac{-\alpha}{\sqrt{K_{ii}}} (\underline{v}_S - \underline{v}_D) \cdot \underline{\tau}_i \\ p_D = p_S \end{cases} \quad (3) $$

where $\underline{n} = \underline{n}_S = -\underline{n}_D$ is the outward normal vector of Stokes/Darcy interface, $\underline{\tau}_i$ refers to the unit tangential vectors to the interface and $K_{ii}$ are the corresponding permeabilities.

Numerically, there are several ways to formulate the coupled problem. Indeed, the form of the Darcy problem allows us to separately solve for pressure and velocity fields [17]. Those methods consist in finding the pressure field and then posttreating the velocity with Darcy's law. However, those approaches imply velocity discontinuities if permeability changes within the domain, leading to mass conservation issues. This is even more drastic for Stokes/Darcy coupling where velocities may strongly differ in the vicinity of the Stokes/Darcy interface proportionally to the permeability jump. As for the Stokes equations, they are ideally solved using a mixed velocity/pressure formulation to ensure mass conservation. For those reasons, a mixed velocity/pressure formulation is chosen for both Stokes and Darcy problems. These formulations can take two different forms in which velocity and pressure boundary conditions appear as Dirichlet or Neumann conditions. A convenient choice for our approach is the dual mixed formulation for Stokes and primal mixed formulation for Darcy [2,3] presented below.

<figure><img src="image_2.png" /><figcaption>Fig. 1. Stokes-Darcy domains and 2D coupling conditions.</figcaption></figure>
The dual mixed velocity/pressure formulation of Stokes equation writes as follows:

$$
\left\{
\begin{array}{l}
\text{Find } (v, p) \in E_{\underline{v}} \times L^2(\Omega_S) \mid \forall (\delta v, \delta p) \in E_{\delta v} \times L^2(\Omega_S), \\
\int_{\Omega_S} 2\eta \underline{\underline{\cdot}}(\delta v) : \underline{\underline{D}}(v) dv - \int_{\partial\Omega_S} p \underline{\underline{v}} \cdot \delta v \, dv = \underbrace{\int_{\partial\Omega_{S,i}} 2\eta \delta v \underline{\underline{\cdot}} \underline{\underline{D}}(\delta v) \cdot \underline{\underline{n}} \, ds}_{\int_{\partial\Omega_{S,i}} \delta v \underline{\underline{\cdot}} \underline{\underline{n}} \, ds} - \int_{\partial\Omega_{S,i}} p \delta v \cdot \underline{\underline{n}} \, ds \\
\int_{\Omega_S} \delta p \underline{\underline{\nabla}} \cdot \underline{\underline{v}} \, dv = 0 \\
E_{\delta v} = \{ \delta v \in H^1(\Omega_S) \mid \delta v = 0 \text{ on } \partial\Omega_{S,v} \} \\
E_{\underline{v}} = \{ \underline{\underline{v}} \in H^1(\Omega_S) \mid \underline{\underline{v}} = \underline{\underline{v}}^d \text{ on } \partial\Omega_{S,v} \}
\end{array}
\right.
\qquad (4)
$$

where $\partial\Omega_{S,i}$ is a boundary of the Stokes domain submitted to normal stress $\underline{t} = \underline{\sigma} \underline{n}$ and $\partial\Omega_{S,v}$ is the complementary boundary submitted to a prescribed velocity $\underline{v}^d$.

The primal mixed velocity/pressure formulation of Darcy equation writes as follows:

$$
\left\{
\begin{array}{l}
\text{Find } (v, p) \in L^2(\Omega_D) \times E_p \mid \forall (\delta v, \delta p) \in L^2(\Omega_D) \times E_{\delta p}, \\
\int_{\Omega_D} \delta \underline{v} \frac{\eta}{K} \underline{v} \, dv + \int_{\Omega_D} \delta \underline{v} \underline{\nabla} p \, dv = \int_{\Omega_D} \delta v \rho \underline{g} \, dv \\
\int_{\Omega_D} \underline{v} \underline{\nabla} \delta p \, dv = \int_{\partial \Omega_{D,v}} \delta p \underline{v} \underline{n} \, ds \qquad (5) \\
E_{\delta p} = \{ \delta p \in H^1(\Omega_D) \mid \delta p = 0 \text{ on } \partial \Omega_{D,p} \} \\
E_p = \{ p \in H^1(\Omega_D) \mid p = p^d \text{ on } \partial \Omega_{D,p} \}
\end{array}
\right.
$$

where $\partial\Omega_{D,p}$ is a boundary of Darcy domain submitted to an imposed pressure $p^d$ and $\partial\Omega_{D,v}$ is the complementary boundary submitted to a prescribed normal velocity $\underline{v}^d = \underline{v} \cdot \underline{n}$.

When discretizing those formulations with a velocity-based finite element method, stability problems appear due to the fluid incompressibility. In order to avoid non-physical pressure modes due to this instability, one must consider to choose elements that fulfill the Brezzi-Babuska condition [21]. A good way to respect this condition is to ensure that the degree of approximation for the velocity field is larger than the pressure one [22]. In the literature two families of elements can be found. The first one is the Taylor-Hoods elements type (P2/P1 or P1/P0) [23,24]. The second type corresponds to standard finite element formulations stabilized with the so-called multi-scale stabilization method that consists in formulating the problem on a coarse scale (nodes) and a fine scale (center of the element), the fine scale being possibly condensed in a stabilization term. Among those approaches one can find the Hughes variational multi-scale (HVM) method [7,8] and the mini-element [2,3,25]. In the present work, the mini-element is selected for its robustness and simplicity of implementation. Indeed, it is based on the same degree of interpolation for both pressure and velocity and there is no constant to introduce and to calibrate in the stabilization term contrary to HVM [8] or ASGS methods [26,27].

## 2.1. Mini-element stabilization method for Stokes formulation

In what follows, is described the practice of the stabilization for the dual mixed formulation of Stokes problem, the same procedure is applied to stabilize the primal mixed Darcy formulation. The mini-element stabilization method stands on the enrichment of the velocity field approximation by the introduction of an extra node in the center of the element (see Fig. 2). Since the velocity field degree of interpolation is higher than the pressure one, the Brezzi-Babuska condition is fulfilled avoiding pressure parasitic modes [21].

Let us introduce the superscript h to distinguish the approximated fields from the exact ones. The extra velocity degree of freedom leads to the following decomposition for the velocity approximation $(\underline{v}^h(x))$ and the corresponding weighting function

$(\delta \underline{v}^h(x))$, taken in the functional spaces described in Eq. (8):

$$
\begin{equation}
\begin{aligned}
\underline{v}^h(x) &= \underline{v}_I^h(x) + \underline{v}_B^h(x) = \sum_{i=1}^{n_e} [N_i(x)] \{\underline{v}_I^h\} + [b_e(x)] \{\underline{v}_B^h\} \\
\delta \underline{v}^h(x) &= \delta \underline{v}_I^h(x) + \delta \underline{v}_B^h(x) = \sum_{i=1}^{n_e} [N_i(x)] \{\delta \underline{v}_I^h\} + [b_e(x)] \{\delta \underline{v}_B^h\}
\end{aligned}
\tag{6}
\end{equation}
$$

where $N_i(x)$ are the shape functions of the element equal to 1 at node i and 0 at the other nodes, $n^e$ is the number of nodes in the element and $b_e$ is the bubble shape function that vanishes on the edge of the element (see Fig. 2). The following decomposition for the pressure approximation $(p^h(x))$ and the corresponding weighting function $(\delta p^h(x))$ is proposed (functional spaces are described in Eq. (8)):

$$
\begin{align}
p^h(x) &= \sum_{i=1}^{n_e} N_i(x) p_i^h \nonumber \\
\delta p^h(x) &= \sum_{i=1}^{n_e} N_i(x) \delta p_i^h \tag{7}
\end{align}
$$

Introducing the previous discretization in system (4), the weak form of the Stokes problem writes as follows:

$$
\left\{
\begin{aligned}
& \text{Find } (\underline{v}_I^h, \underline{v}_B^h, p^h) \in E_{\underline{v}_I^h} \times E_{\underline{v}_B^h} \times L^2(\Omega_S) \mid \forall (\delta \underline{v}_I^h, \delta \underline{v}_B^h, \delta p^h) \in E_{\delta \underline{v}_I^h} \times E_{\delta \underline{v}_B^h} \times L^2(\Omega_S), \\
& \int_{\Omega_S} 2\eta \underline{\underline{D}}(\delta \underline{v}_I^h) : \underline{\underline{D}}(\underline{v}_I^h + \underline{v}_B^h) \, dv - \int_{\Omega} p^h \underline{\underline{\nabla}} \cdot \delta \underline{v}_I^h \, dv = \int_{\partial \Omega_{S,i}} \delta \underline{v}_I^h \cdot \underline{\underline{\sigma}} \cdot \underline{n} \, ds \\
& \int_{\Omega} 2\eta \underline{\underline{D}}(\delta \underline{v}_B^h) : \underline{\underline{D}}(\underline{v}_I^h + \underline{v}_B^h) \, dv - \int_{\Omega} p^h \underline{\underline{\nabla}} \cdot \delta \underline{v}_B^h \, dv = 0 \\
& \int_{\Omega} \delta p^h \underline{\underline{\nabla}} \cdot \underline{\underline{v}}^h \, dv + \int_{\partial \Omega} \delta p^h \underline{\underline{\nabla}} \cdot \underline{\underline{v}}_B^h \, dv = 0 \\
& E_{\delta \underline{v}_I^h} = \{ \delta \underline{v}_I^h \in H^1(\Omega_S) \mid \delta \underline{v}_I^h = 0 \text{ on } \partial \Omega_{S, \underline{v}_I^h} \} \\
& E_{\underline{v}_I^h} = \{ \underline{v}_I^h \in H^1(\Omega_S) \mid \underline{v}_I^h = \underline{v}_I^{hd} \text{ on } \partial \Omega_{S, \underline{v}_I^h} \} \\
& E_{\delta \underline{v}_B^h} = E_{\underline{v}_B^h} = \{ \underline{v}_B^h \in H^1(\Omega_S) \mid \underline{v}_B^h = 0 \text{ on } \partial \Omega_{S, \underline{v}_B^h} \}
\end{aligned}
\right.
\tag{8}
$$

After few mathematical steps described in [21], the discrete stabilized Stokes system can be written under the following form:

$$
\left[ \begin{bmatrix} K_{ll} & K_{lp} \\ K_{pl} & C_{pp} \end{bmatrix}^T \right] \left\{ \begin{matrix} \{ \underline{v}_I \} \\ \{ p \} \end{matrix} \right\} = \left\{ \begin{matrix} \underline{F} \\ \underline{0} \end{matrix} \right\} \qquad (9)
$$

where $[C_{pp}] = -[K_{pb}][K_{bb}]^{-1}[K_{bp}]^T$ is the stabilization term. Finally, one can summarize the method to the introduction of a stabilization term $[C_{pp}]$ inversely proportional to the viscosity $\eta$ that relaxes the fluid incompressibility constrain. In the same way, one can stabilize the Darcy formulation with a term proportional to $K/\eta$ [2,3].

This stabilization method has been validated on test cases in [2,3] for hierarchical bubble functions (P1+, see Fig. 2) and its good rate of convergence verified with the method of the manufactured exact solution (MMES) proposed in [28]. In Section 3, are presented an analytical validation of the coupled problem and an illustration of the MMES validation.

## 2.2. Numerical aspects of Stokes/Darcy coupling

The solution selected for coupling Stokes and Darcy flows consists in dealing with two sub-meshes, respectively for Stokes and Darcy domains, non-necessarily structured but matching at the interface. The flows are solved individually, in an iterative scheme, until equilibrium is reached at their interface. Practically, normal velocity taken from Darcy's domain, $\underline{v}_D \cdot \underline{n}$, is imposed by a penalty method into Stokes equations through a Dirichlet boundary condition in velocity, and the hydrostatic pressure taken in the Stokes' domain, $p_s$, is prescribed into Darcy equations as a Dirichlet boundary condition in pressure (see points 2 and 3 in Fig. 12). Convergence is reached when both pressure and velocity corrections tend toward zero, lower than a given admissible bound.

Numerically, the coupling methods introduced imply to pay attention in the construction of the global system to be solved, let
<figure><img src="image_3.png" /><figcaption>Fig. 2. Mini-element P1+/P1: (a) triangle element (2D), (b) graphical representation of the hierarchical bubble shape function [2,3].</figcaption></figure>

us say $[\kappa]\{u\} = \{F\}$. First of all, the kinematic boundary conditions (tangential and normal velocity) prescribed by penalty method lead to the following system for a general interface defined by its normal components ($\underline{n} = n_x \underline{x} + n_y \underline{y} + n_z \underline{z}$):

$$ \begin{pmatrix} \cdots & \cdots & \cdots & \cdots & \cdots \\ \cdots & A & B & C & \cdots \\ \cdots & D+Qn_x & E+Qn_y & F+Qn_z & \cdots \\ \cdots & G & H & I & \cdots \\ \cdots & \cdots & \cdots & \cdots & \cdots \end{pmatrix} \begin{Bmatrix} v_x^i \\ v_y^i \\ v_z^i \\ \vdots \end{Bmatrix}_{\text{global}} = \begin{Bmatrix} \vdots \\ L \\ M+Qv_D \cdot \underline{n} \\ N \\ \vdots \end{Bmatrix} \quad (10) $$

where $Q$ is a penalty coefficient usually taken as $10^{8-12}\max(\kappa_{ij})$. However, on one hand high values of $Q$ ensure that the condition is well fulfilled, while on the other hand they will lead to an ill conditioned system due to extra out-of diagonal terms for which iterative solvers poorly behave. To ensure a good conditioning, the Dirichlet conditions can be prescribed by penalty but in a local reference system attached to the nodes where normals to the frontier (or the interface) have to be considered as one of the axes of this local reference system (see Fig. 3). In practice, the velocity degrees of freedom are locally rotated. For an interface of normal $\underline{n} = n_x\underline{x} + n_y\underline{y}$ the system at the elementary level will become

$$ \begin{bmatrix} [I] & [0] \\ [0] & [R] \end{bmatrix}^T \begin{bmatrix} [\kappa^e] \\ [0] & [R] \end{bmatrix} \begin{Bmatrix} \{V_g\} \\ \{V_l\} \end{Bmatrix} = \begin{bmatrix} [I] & [0] \\ [0] & [R] \end{bmatrix}^T \{f^e\} \quad (11) $$

where $[R]$ is the rotation matrix from the global to the local reference system, $\{V_g\}$ is the unknown nodal velocity vector expressed in the global reference system and $\{V_l\}$ is the unknown nodal velocity vector expressed in the local reference system. The local problem thus writes now in a dominant-diagonal form

$$ \begin{pmatrix} \cdots & \cdots & \cdots & \cdots & \cdots \\ \cdots & An_y - Bn_x & An_x + Bn_y & C & \cdots \\ \cdots & Dn_y - En_x & Dn_x + En_y + Q & F & \cdots \\ \cdots & Gn_y - Hn_x & Gn_x + Hn_y & I & \cdots \\ \cdots & \cdots & \cdots & \cdots & \cdots \end{pmatrix} \begin{Bmatrix} \vdots \\ v_x^i \\ v_y^i \\ v_z^i \\ \vdots \end{Bmatrix}_{\text{local}} = \begin{Bmatrix} \vdots \\ L \\ M+Qv_D \cdot \underline{n} \\ N \\ \vdots \end{Bmatrix} \quad (12) $$

After solving the system, the velocity needs to be post-treated to be expressed in the global reference system. It must be noted that due to the relatively small permeabilities considered in this study, the Darcy velocity to impose as Dirichlet boundary conditions in Stokes formulation is small. So in order to penalize properly the system, the penalty coefficient $Q$ should be chosen regarding permeability to ensure $M \ll Qv_D \cdot n$. It was verified for permeabilities down to $10^{-15} \text{ m}^2$ that $Q = 10^5/K$ for $K < 10^{-6} \text{ m}^2$ leads to satisfactory results [27,28].

Secondly, an important issue concerns the use of the mini-element P1+/P1. Indeed, as shown in Section 2.1, the stabilization method relies on the relaxation of the fluid incompressibility

<figure><img src="image_4.png" /><figcaption>Fig. 3. Kinematic boundary condition.</figcaption></figure>

constrain equivalent to

$$ \underline{\nabla} \cdot \underline{v} + \beta p = 0 \quad (13) $$

with $\beta$ the stabilization term proportional to $1/\eta$ for Stokes and to $K/\eta$ for Darcy. When considering a regular Stokes flow $\beta$ is sized by the viscosity ($\eta$) and remains sufficiently low to consider $\beta p$ negligible, ensuring the fluid incompressibility ($\nabla \cdot \underline{v} \approx 0$). However when this Stokes' flow is coupled to a porous flow with low permeability, the Darcian normal velocity to be imposed into the Stokes formulation can be extremely low, down to $10^{-13} \text{ m s}^{-1}$, implying small velocities over the whole domain. In such a case, typical in composite process simulation, the stabilization term $\beta p$ may not be negligible at all, leading to a non-zero velocity divergence in Stokes domain. The incompressibility constrain cannot be enforced anymore. To ensure an acceptable value of the velocity divergence, it may be necessary to introduce an artificial viscosity in the computation of the stabilization term $\beta$ in Eq. (13) or $[C_{pp}]$ in the global system (9). It has been found that an artificial viscosity equal to 10 times the real one guarantees good results down to permeability values of $10^{-15} \text{ m}^2$. This is another improvement proposed in this paper, since it allows us to deal with real severe physical applications as shown in validation cases and in the final example of Section 7.

In order to demonstrate the benefit of the artificial viscosity a comparison was made on a simple test case. This test case is described in Fig. 4 and consists of a flow perpendicular to the Stokes/Darcy interface. Given the boundary conditions, this case has a simple analytic solution $\underline{v} = -1.76 \times 10^6 K_{yy}$. Fig. 5 shows the ratio of simulated velocity magnitude over analytical velocity magnitude without and with artificial viscosity for different permeability values.

### 3. Validation of the Stokes/Darcy coupling method

#### 3.1. Analytical validation: the parallel flow test case

To validate the Stokes/Darcy coupling, we investigate the test case of a flow parallel to the interface in a steady state, presented in Fig. 6.
<figure><img src="image_5.png" /><figcaption>Fig. 4. Transverse flow test case.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 5. Comparison of the ratio of simulated velocity magnitude over analytical velocity magnitude in the case of a transverse flow: (a) without artificial viscosity, (b) with artificial viscosity.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Fig. 6. Parallel flow test case for Stokes/Darcy coupling.</figcaption></figure>

Let us consider a porous medium of permeability $K = 10^{-15}$ m² and dimensions $H_D = 2$ m and $L=5$ m crossed by a fluid of viscosity $\eta=1$ Pa s submitted to a pressure differential of $p_1 - p_0 = 10$ Pa in contact with a pure fluid layer of thickness $H_S = 1$ m and submitted to the same pressure differential. The remaining boundaries are considered as impervious walls $v \cdot n = 0$. This test case has an analytical solution [7]:

$$ 
\begin{cases}
v_x = -\frac{K}{2\eta} \frac{\lambda^2 + 2a\lambda}{1+a\lambda} \right) dP \left(1 + \frac{\alpha}{\sqrt{K}y}\right) + \frac{1}{2\eta} (y^2 + 2ay\sqrt{K}) \frac{dP}{dx} & \text{for } 0 \le y \le H_S \\
v_x = -\frac{K}{\eta} \left(\frac{p_1 - p_0}{L}\right) & \text{for } y < 0 \\
v_y = 0 & \text{for any y}
\end{cases}
\quad (14) $$

with $\lambda = H_S/\sqrt{K}$, $dP/dx = (p_1 - p_0)/L$ and $\alpha$ being the drag coefficient.

Fig. 7 shows the results obtained for four values of $\alpha$, the coefficient controlling the tangential velocity on the Stokes/Darcy interface (Eq. (3)). It confirms the good implementations of the coupling methods. One can verify the effect of this coefficient leading to Poiseuille's flow in the extreme case ($\alpha$ large). Finally, one can remark that $\alpha$ should be chosen regarding permeability, in the present case we chose $\alpha=0$, $\alpha = \sqrt{K}$, $\alpha = 10^2\sqrt{K}$ and $\alpha = 10^4\sqrt{K}$. In Fig. 8, the pressure and velocity fields computed for $\alpha=0$, corresponding to Fig. 7(a) case, are presented.

The numerical simulations were performed over a mesh of 3094 triangles made up of two non-necessarily structured meshes matching at the interface. The convergence, corresponding to Stokes/Darcy equilibration, was reached (for the four cases) in three iterations with a relative error between the two last iterations of the order of $10^{-9}$.
<figure><img src="image_8.png" /><figcaption>Fig. 7. Stokes/Darcy coupling analytical (continuous line)/numerical (dots) comparison: (a) $\alpha = 0$, (b) $\alpha = 3.16 \times 10^{-8}$, (c) $\alpha = 3.16 \times 10^{-6}$, (d) $\alpha = 3.16 \times 10^{-4}$.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Fig. 8. Simulation results for $\alpha = 0$: (a) pressure field, (b) velocity field.</figcaption></figure>

## 3.2. Method of the manufactured exact solutions

The proper implementation of the coupling algorithm was verified with the method of the manufactured exact solutions (MMESs) [1]. This method consists in building an analytical solution, that is then introduced into the analytical system of PDE to be solved, to deduce the corresponding source term. Finally, the finite element problem is loaded with this source term, and solved for the discrete solution that can be compared with the exact analytical one, to assess the accuracy.

Let us consider a domain $\Omega = (0, 1) \times (0, 2)$ m² divided into two subdomains $\Omega_S = (0, 1) \times (1, 2)$ m² governed by the Stokes' equation and $\Omega_D = (0, 1) \times (0, 1)$ m² governed by Darcy's law in contact along the interface $\Gamma = \Omega_S \cap \Omega_D$. The pressure and velocity fields proposed by Discacciati in [6] are considered:

$$ \begin{cases} \nabla_S \cdot \mathbf{x} = - \cos(\pi/2y) \sin(\pi/2x) \\ \nabla_S \cdot \mathbf{y} = \cos(\pi/2x) \sin(\pi/2y) - 1 + x \\ p_S = 1 - x \end{cases} \tag{15} $$
a

<figure><img src="image_10.png" /><figcaption>Fig. 9. MMES results with η=1 Pa s, K=1 m², h=0.0125 m: (a) pressure field, (b) velocity field corresponding to Eqs. (15) and (16).</figcaption></figure>

b

<figure><img src="image_11.png" /></figure>

$$
\begin{aligned}
\left| \frac{\mathbf{v}_D \cdot \mathbf{x}}{\mathbf{v}_D \cdot \mathbf{y}} \right| &= \sin(\pi/2x) \cos(\pi/2y) + y \\
\left| \frac{\mathbf{v}_D \cdot \mathbf{y}}{p_d} \right| &= \cos(\pi/2x) \sin(\pi/2y) + x - 1 \\
\left| \frac{p_d}{2/\pi} \right| &= \frac{1}{\pi} \cos(\pi/2x) \cos(\pi/2y) - y(x - 1)
\end{aligned}
\qquad (16)
$$

Let us consider regular meshes constructed with squares of size $h \times h$ split into two triangles. Each computation was made over a more refined mesh ($h$ is divided by 2) than the previous one starting with $h=0.1$ down to $h=0.0125$. For the sake of simplicity we chose $\eta=1$ Pa s and $K=1$ m² as proposed in [6]. Fig. 9 shows isovalues of the numerical solution for pressure and for velocity, for $h=0.0125$.

For every mesh and each numerical solution $v_h$ and $p_h$, the errors are calculated using the $L^2$ norm:

$$ \| u \|_{0,\Omega} = \left( \int_{\Omega} u^2 d\Omega \right)^{1/2} \qquad (17) $$

and $H^1$ norm:

$$ \| u \|_{1,\Omega} = \| u \|_{0,\Omega}^2 + \sum_{j=1}^{d} \| \frac{\partial u}{\partial x_j} \|_{0,\Omega}^2 \qquad (18) $$

Errors thus calculated permit to compare the numerical rates of convergence to the theoretical ones. Stokes rate of convergence should verify [29]:

$$ \| \underline{v} - \underline{v}_h \|_{1,\Omega_s} + \| p - p_h \|_{0,\Omega_s} \le C_1 h (\| \underline{v} \|_{2,\Omega_s} + \| p \|_{1,\Omega_s}) \qquad (19) $$

and Darcy rate of convergence should verify [29]:

$$ \| \underline{v} - \underline{v}_h \|_{1,\Omega_D} + \| p - p_h \|_{0,\Omega_D} \le C_2 h (\| \underline{v} \|_{2,\Omega_D} + \| p \|_{1,\Omega_D}) \qquad (20) $$

where $C_1$ and $C_2$ are some constants. Convergence of the solution was studied and computed errors for pressure and velocity for the different mesh sizes are presented in Tables 1 and 2 for Stokes and Darcy domain respectively. Fig. 10 shows that pressure and velocity converge respectively with rates of 1.8 and 2.2 in Stokes domain and with rates of 2 and 1.4 in Darcy domain. These results correspond to the theoretical and numerical convergence rates determined in [29] with Eqs. (19) and (20).

<caption>**Table 1**</caption>

Errors on the $L^2$ norm for the pressure and velocity, in Stokes' domain for coupled problem.

<table><thead><tr><th>Mesh</th><th>h (m)</th><th>||p - p<sub>h</sub>||<sub>0,Ω</sub></th><th>||v - v<sub>h</sub>||<sub>0,Ω</sub></th></tr></thead><tbody><tr><td>10 × 10</td><td>0.1</td><td>0.0017176</td><td>4.93044e-5</td></tr><tr><td>20 × 20</td><td>0.05</td><td>0.00045979</td><td>1.01198e-5</td></tr><tr><td>40 × 40</td><td>0.025</td><td>0.000121</td><td>2.19269e-6</td></tr><tr><td>80 × 80</td><td>0.0125</td><td>3.27068e-5</td><td>5.00105e-7</td></tr></tbody></table>

Finally, the good convergence rates and the good agreement with analytical results allow us to validate the implementation of the coupling algorithm and to confirm its robustness.

## 4. Formulation of the solid mechanics problem

The solid mechanics problem is non-linear on two counts, first due to the non-linear response of the porous medium and second due to the finite deformations. An updated Lagrangian formulation in displacement is adopted here [2,3] which writes

$$
\begin{align*}
& \text{find } \underline{u} \in E_u \mid \forall \delta \underline{u} \in E_{\delta u}, \\
& \int_{t_1 \le t \le \Omega} \underline{u}^{t+\Delta t} \underline{G} : \delta \underline{u}^{t+\Delta t} \underline{E} \underline{u}^{t+\Delta t} dv - \int_{t_1 \le t \le \Omega} \underline{u}^{t+\Delta t} \underline{f}^{u} \cdot \delta \underline{u}^{t+\Delta t} dv - \int_{\partial u^{t+\Delta t} \Omega_{\underline{g}}} \underline{u}^{t+\Delta t} \underline{f}^{s} \cdot \delta \underline{u}^{t+\Delta t} ds = 0 \\
& \underline{u} = \underline{u}^{d}, \quad \forall x \in \partial \Omega_{\underline{u}} \\
& E_{\underline{u}} = \{ \underline{u} \in H^1(\Omega) \mid \underline{u} = \underline{u}^d \quad \forall x \in \partial \Omega_{\underline{u}} \} \\
& E_{\delta \underline{u}} = \{ \delta \underline{u} \in H^1(\Omega) \mid \delta \underline{u} = 0 \quad \forall x \in \partial \Omega_{\underline{u}} \}
\end{align*}
\qquad (21)
$$

where $S$ is the second Piola–Kirchoff stress tensor, $\delta E$ is the virtual Green–Lagrange tensor, $f^s$ are the contact forces acting on the solid surface and $\delta u$ is the virtual displacement field. This formulation relies on a Newton–Raphson scheme to deal with those nonlinearities. Actually, the hyper-elastic orthotropic constitutive relationship is solved in a local reference system for each element to
take into account the material orientations. The procedure consists of a first elastic prediction:

$$[M_T((0))](u_0) = \{R\} \qquad (22)$$

where $[M_T((0))]$ is the initial tangent matrix computed with the stiffness taken at the origin of the constitutive law and $\{u_0\}$ is the unknown displacement at the prediction iteration. Then convergence is checked with the relative residual as follows:

$$\frac{\|[M_S((u_{i-1}))]\{u_{i-1}\} - \{R\}\|}{\|\{R\}\|} \le r \quad \text{for } i \ge 1 \qquad (23)$$

where $[M_S((u_{i-1}))]$ is the current secant matrix computed from the constitutive response, $\{u_{i-1}\}$ the displacement at previous iteration and $r$ the target component residual value. If convergence is not reached, a correction increment is computed:

$$\begin{cases} [M_T((u_{i-1}))]\{\Delta u_i\} = [M_S((u_{i-1}))]\{u_{i-1}\} - \{R\} \\ \{u_i\} = \{u_{i-1}\} + \{\Delta u_i\} \end{cases} \qquad (24)$$

where $\{\Delta u_i\}$ is the correction increment at iteration $i$, $\{u_i\}$ is the displacement at iteration $i$. This procedure is encapsulated in the block 1 and 4 of the global algorithm (see Fig. 12) used to solve the entire fluid-solid problem.

In addition, the porous medium is represented by an equivalent homogeneous medium composed of rigid fibers, so that each macroscopic strain is reflected by fiber rearrangements at the microscopic scale (see Fig. 11). Therefore the porosity can be directly computed from the mechanical rate of deformations.

Indeed, mass conservation of the homogeneous porous medium writes as follows:

$$\frac{D}{Dt}(m(\Omega)) = \frac{D}{Dt}\int_{\Omega_t} \rho \, dv = 0 \qquad (25)$$

where $D/Dt$ is the particular derivative with respect to time and $\rho$ is the density of the porous medium. This porous medium is actually a biphasic medium constituted of fibers and fluid. Using the mixture low and noting that there is no mass transfer between

<caption>**Table 2**</caption>

Errors on the $L^2$ norm for the pressure and velocity in Darcy's domain for coupled problem.

<table><thead><tr><th>Mesh</th><th>h (m)</th><th>||p - p_h||<sub>0,Ω</sub></th><th>||v - v_h||<sub>0,Ω</sub></th></tr></thead><tbody><tr><td>10 × 10</td><td>0.1</td><td>0.00130628</td><td>0.00499473</td></tr><tr><td>20 × 20</td><td>0.05</td><td>0.00032102</td><td>0.00129536</td></tr><tr><td>40 × 40</td><td>0.025</td><td>6.90808e-5</td><td>0.00033242</td></tr><tr><td>80 × 80</td><td>0.0125</td><td>1.39862e-5</td><td>0.00011783</td></tr></tbody></table>

<figure><img src="image_12.png" /></figure>

<figure><img src="image_13.png" /><figcaption>**Fig. 10.** Error convergence in Stokes and Darcy domain: (a) velocity, (b) pressure.</figcaption></figure>

the fluid and the fibers Eq. (25) leads to

$$\left| \frac{D}{Dt} \int_{\Omega_t} \rho_s (1 - \phi) \, dv \right| = 0 \\ \left| \frac{D}{Dt} \int_{\Omega_t} \rho_r \phi \, dv \right| = 0 \qquad (26)$$

with $\rho_s$ being the fiber density, $\rho_r$ the fluid density and $\phi$ the porosity or fluid fraction corresponding to micro- and macro-voids in Fig. 11, i. e., $\phi = 1 - V_f$ with $V_f$ being the fiber fraction. The first equation of system (26) is actually the mass conservation of the solid phase, i.e., the fibers. This equation will be used to update the porosity ($\phi$) from the displacement field. Using a Lagrangian approach, this equation can be written in the reference configuration as follows:

$$\frac{D}{Dt} \int_{\Omega_{ref}} \rho_s (1 - \phi) J \, dV_{ref} = 0 \qquad (27)$$

where $J$ denotes the Jacobian of the transformation between the reference configuration $\Omega_{ref}$ and the current configuration $\Omega_t$ and $dV_{ref}$ is the elementary volume defined on $\Omega_{ref}$. Since, the mechanical problem is also treated with a Lagrangian approach, the reference domain is the domain of the material particles, so, the convective velocity is zero. In addition, we can note that $dV_{ref}$ is constant over time. So Eq. (27) becomes

$$\int_{\Omega_{ref}} \frac{\partial}{\partial t}(\rho_s(1-\phi))J \, dV_{ref} = 0 \qquad (28)$$

Furthermore, we consider that the fibers are incompressible, so

$$\frac{\partial \rho_s}{\partial t} = 0 \qquad (29)$$

leading to

$$\int_{\Omega_{ref}} \frac{\partial}{\partial t}((1-\phi)J) \, dV_{ref} = 0 \qquad (30)$$

Finally, Eq. (30) means that the product $(1-\phi)J$ is constant over time, so it can be reduced to the following relation:

$$J(\bar{x}, t + \Delta t)(1 - \phi(\bar{x}, t + \Delta t)) = J(\bar{x}, t)(1 - \phi(\bar{x}, t)) \qquad (31)$$

where $J$ is the Jacobian of the transformation, $\bar{x}$ the position vector and $\phi$ the porosity of the porous medium, i.e., $\phi = 1 - V_f$ with $V_f$ being the fiber volume fraction. This approach proposed by Celle et al. in [2,3] and also used later in [19] does not restrict the analysis to one dimension, as it is usually the case in the literature when using semi-empirical laws expressing porosity or thickness variation as a function of pressure [14–16]. Moreover, it allows us to use any constitutive law for the porous medium (elastic or hyper-elastic and isotropic or orthotropic). This is one of the major feature of the presented global approach, opening new prospects in terms of industrial outcomes.
This formulation has been validated on analytical test cases and by comparison with the commercial solver Abaqus [2,3]. In the present paper, we focus on the fluid/structure validation (see Section 5.1).

## 5. Fluid/structure coupling method

The fluid/structure strong coupling is treated iteratively in a quasi-static context (see Fig. 12). Indeed, the Stokes/Darcy coupled

problem (see points 2 and 3 in Fig. 12) is solved on a fixed grid (i.e. rigid porous medium), the fluid influence on the porous medium being taken into account through pore (interstitial) pressure by Terzaghi's law [10]:

$$ \sigma_{tot}(\underline{\epsilon}) = \sigma_{eff}(\underline{\epsilon}) + SpI \qquad (32) $$

where $\sigma_{tot}$ is the overall stress applied to the system, $\sigma_{eff}$ is the effective stress in the fibers, $I$ is the unity tensor, $p$ is the

<figure><img src="image_14.png" /><figcaption>Fig. 11. Illustration of the compression of a small section of a fiber bed at microscopic level.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Fig. 12. Stokes/Darcy/solid mechanics coupling algorithm.</figcaption></figure>
hydrostatic pressure and S is the saturation factor of the pore. Let us notice that the updated Lagrangian scheme adopted for the solid mechanics problem permits us to use the quasi-Eulerian stress field to compute the corresponding Terzaghi stress tensor [2,3]. The non-linear finite strain mechanical problem (see point 4 in Fig. 12) is then solved considering a given pore pressure (i.e. for a Terzaghi equivalent behavior) acting as an internal load. Then, the deformation is reflected by the porosity variation modifying the permeability $K$ of the medium. The permeability can be computed with any porosity dependent measurement or model. Here, for the example, the Carman-Kozeny's law [30] can be considered:

$$K = \frac{d_f^2}{16h_k(1-\phi)^2} \phi^3 \quad (33)$$

where $d_f$ is the fiber diameter, $h_k$ the Kozeny constant and $\phi$ the porosity. Finally, convergence is verified by a relative error computation on both fluid velocity and pressure fields of fluid, and solid displacement between two iterations.

### 5.1. Analytical validation of the fluid/structure interaction

Using an analytical approach based on Airy's functions (see Appendix A), it is possible to build a test case for Terzaghi's principle which has an analytical solution in the case of linear elastic preforms. The test case represents the compaction of a porous medium submitted simultaneously to an internal loading due to the pressure of the fluid present in the pores (see Fig. 22) [2,3] and an external loading acting on the upper side.

We consider an isotropic material of Young's modulus $E=266$ kPa and Poisson's ratio $\nu = 0.3$. The plate is submitted to an external load $P_{imp} = 10^4$ Pa, a linear internal load $P_{int}(y) = -5 \times 10^5 y + 10^4$ Pa. Dimensions of the domain are $H = 2 \times 10^{-2}$ m and $L = 4 \times 10^2$ m.

The analytical solution is

$$ \begin{cases} u_x(x,y) = \frac{1+\nu}{E} [(1-2\nu)(ay+b)x+\nu Gx] \\ u_y(x,y) = \frac{1+\nu}{E} \left[ (1-2\nu)\left(\frac{a}{2}y+b\right)y + (\nu-1)Gy - \frac{1-2\nu}{2}ax^2 \right] \end{cases} \quad (34) $$

where $G = aH+b+P_{imp}$, $a = -5 \times 10^5$ and $b = 10^4$.

Fig. 13 shows the perfect agreement between analytical and numerical results.

## 6. Flow front tracking

There are several methods for tracking the flow front, such that the level-set method and the FE-CV method [31]. As part of this work, we did not study these methods but have used the filling algorithm already present in PAM-RTM™ [32,33]. This algorithm is based on a division of the transient regime in a series of quasi-static states. New elements filled between two time steps are determined from the flow rate calculated at the last known time step. A filling factor varying from 0 to 1 is associated with each element. In the case of flow in deformable porous medium, porosity changes over time. Because the time increment is determined before the calculation, it is necessary to correct it so that the total amount of resin present in the system corresponds to the volume of saturated pore. The following correction is proposed:

$$ \Delta t = \Delta t_{algo} \frac{V_{t+\Delta t}}{V_t} \quad (35) $$

where $\Delta t$ is the corrected time increment, $\Delta t_{algo}$ is the time increment given by the filling algorithm, $V_t$ is the volume of resin that would have been injected in the system during the time step with no deformation and $V_{t+\Delta t}$ is the actual injected volume of resin during the time step.

<figure><img src="image_16.png" /><figcaption>Fig. 13. Terzaghi's test case results: (a) analytical and numerical final shape of the plate, (b) displacement magnitude field over the deformed domain.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>Fig. 14. Stiffened panel geometry (Daher-Socata).</figcaption></figure>
## 7. Example of Liquid Resin Infusion simulation

In order to demonstrate the ability of the approach to handle industrial cases, the complete infusion simulation of a "T" stiffened panel is presented. The geometry, proposed by Daher-Socata within the INFUCOMP project, consists of the assembly of a 380 × 280 × 3 mm³ plate and a "T" stiffener which is made of a two 3 mm thick "L" preforms as described in Fig. 14.

Infusion condition is presented in Fig. 15. One distribution medium is located between the panel and the mould and a second one is placed on the top of the "T" under the vacuum bag. Resin

inlet is located in the center of the plate under the mould while the vent is in the center of the second distribution medium on the top of the "T".

Preforms are made of 48,302 woven fabric from Hexcel Reinforcements constituted of T700 12k carbon fibers. Orthotropic permeability measured by Hexcel Reinforcements is

$$ 
\begin{array}{ll}
|K_p| = 9.76 \times 10^{-12} \text{ m}^2 & \text{in the plane of the fabrics} \\
K_t = 6.75 \times 10^{-13} \exp(-9.8V_f) \text{ m}^2 & \text{in the transverse direction}
\end{array}
\qquad (36) $$

where $V_f$ is the fiber volume fraction. The resin used for this infusion is the RTM6, epoxy from Hexcel whose viscosity is 0.033 Pa s at 120 °C.

Due to the symmetry of the problem, we simulate only one-quarter of the geometry as described in Fig. 17(a). The mesh used is composed of 8100 tetrahedrons and 1884 nodes. The distribution medium at the bottom of the preform is represented by a pure resin flow (Stokes area) and is 1 mm thick. It is considered as non-deformable. The preform is considered as an equivalent porous homogeneous medium (Darcy area), whose constitutive law has not been studied. However, for the example, we will assume that it obeys to the same behavior as fabric G1157 (UD plain) from Hexcel Reinforcement that has been extensively studied in [34–36]. The constitutive law is presented in Fig. 16. The behavior of a stack of UD, not preformed (G1157 tested by P. Wang), is much more flexible than a preform held in shape by an epoxy powder (industrial preforms), so deformations presented in the results are larger than the ones usually observed. Injection is represented by a pressure boundary condition of 1 bar and the vent by a pressure boundary condition of 0 bar. Finally, the other boundary conditions are impervious walls ($v \cdot n = 0$). The boundary conditions are summarized in Fig. 17(b) and (c).

Figs. 18 and 19 show the filling ratio at times $t=63$ s, $t=273$ s and $t=550$ s. The results reflect a fast filling of distribution medium (1 s), then a transverse infusion of the plate. When the resin starts to infuse the plate, beginning of the swelling can be identified (see Fig. 20). After about one minute, the plate is completely filled, the pressure is balanced in its extremities and the preform swells back to its initial thickness (see Fig. 19 at $t=69$ s). Then, the resin begins to migrate in the base of the "T",

<figure><img src="image_18.png" /><figcaption>Fig. 15. Infusion condition description.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>Fig. 16. Compression response of a dry stack made of plain UD G1157.</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>Fig. 17. Stiffened panel infusion: (a) problem geometry (mesh used in blue), (b) fluid boundary conditions, (c) solid boundary conditions. (For interpretation of the references to color in this figure caption, the reader is referred to the web version of this paper.)</figcaption></figure>
first in the transverse direction, and then, when it reaches the
curvature, in the plane of the plies (see Fig. 18 at t=273 s). Finally,
resin finishes impregnating the upper part of the “T” and the filling
ends after 550 s.

Fig. 20 shows the variation of the thickness during the infusion. It
may be noted that the areas of the part that have been impregnated
in the in-plane direction (curvature and upper part of the T) present
a linear swelling along time, while the plate and the base of the T,
which have been impregnated transversely, have three distinct
phases. Indeed, one can observe a first phase of filling with a small
swelling. Then, when resin reaches the vacuum bag (top of the plate),
a sharp increase of thickness occurs. Finally, the preform continues to
swell linearly during the end of the infusion.

8. Conclusion

An innovative and robust iterative method to fully couple Stokes,
Darcy and solid mechanics has been presented, applied here to
composite process simulation. Based on the previous work by Celle
et al., the model has been extended to three dimensional cases and
new coupling conditions have been introduced to improve the global
performance, stability and robustness of the code. All the numerical
methods have been detailed and validated with analytical
approaches and the convergence of the fully coupled Stokes/Darcy
formulation was checked with the method of manufactured exact
solutions. Finally, a numerical example in 3D on a complex shape
showed the capability of the tool to simulate industrial cases.

<figure><img src="image_21.png" /><figcaption>Fig. 18. Filling ratio evolution during the infusion of the stiffened panel (isometrical view).</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>Fig. 19. Filling ratio evolution during the infusion of the stiffened panel (side view), black lines represent the initial geometry.</figcaption></figure>
<figure><img src="image_23.png" /><figcaption>Fig. 20. Thickness variation during the infusion of the stiffened panel.</figcaption></figure>

Concerning forthcoming developments, a full experimental validation is currently conducted with industrial partners in the frame of a FP7 funded project and a transient formulation is studied for the simulation of three-dimensional dynamic consolidation of the porous medium. Moreover, to fully meet industrial needs, thermo-chemical couplings must be introduced allowing us to simulate non-isothermal transient processes and an elastoplastic behavior for the fabric that may be more realistic with respect to compaction experiments.

## Acknowledgment

The authors would like to thank the support of the European Union in the INFUCOMP project (European Community's Seventh Framework Program FP7/2009-2013 under Grant agreement no. 233926) financing partially these researches.

## Appendix A. Terzaghi validation: analytical solution with Airy functions

### A.1. Principle of the Airy functions

Neglecting the inertial forces, the momentum conservation equation can be written as follows:

$$ \nabla \cdot \sigma + f_v = 0 \qquad (A.1) $$

Let us admit that volume forces derive from a potential V as follows:

$$ \mathbf{f}^{\nu} = \begin{cases} f_x \\ f_y \end{cases} = \text{grad } V(x,y) \qquad (A.2) $$

leading to the following form of the momentum equation:

$$ \left| \begin{array}{l} \frac{\partial(\sigma_{xx} + V)}{\partial x} + \frac{\sigma_{xy}}{\partial y} = 0 \\ \frac{\sigma_{xy}}{\partial x} + \frac{\partial(\sigma_{yy} + V)}{\partial y} = 0 \end{array} \right. \qquad (A.3) $$

### A.1.1. Stress field determination

To determine the stress field, we need to find a scalar biharmonic function $\phi(x, y)$ (known as the Airy function) which satisfies

$$ \nabla^4 \phi = \frac{\partial^4 \phi}{\partial x^4} + 2 \frac{\partial^4 \phi}{\partial x^2 \partial y^2} + \frac{\partial^4 \phi}{\partial y^4} = C(\nu) \left( \frac{\partial f_x}{\partial x} + \frac{\partial f_y}{\partial y} \right) \qquad (A.4) $$

with

$$ C(\nu) = \begin{cases} \frac{1-\nu}{1-2\nu} & (\text{plane strains}) \\ \frac{1}{1-\nu} & (\text{plane stresses}) \end{cases} $$

the $\phi$ function needs also to satisfy the stress boundary conditions on the frontiers of the domain:

$$ \begin{aligned} \frac{\partial^2 \phi}{\partial y^2} n_x - \frac{\partial^2 \phi}{\partial x \partial y} n_y &= t_x \\ \frac{\partial^2 \phi}{\partial x^2} n_x - \frac{\partial^2 \phi}{\partial y \partial x} n_x &= t_y \end{aligned} \qquad (A.5) $$

In a 2D case, the plane strains stress field is given by the following relations:

$$ \begin{aligned} \sigma_{xx} &= \frac{\partial^2 \phi}{\partial y^2} - V \\ \sigma_{yy} &= \frac{\partial^2 \phi}{\partial x^2} - V \\ \sigma_{xy} &= \sigma_{yx} = -\frac{\partial^2 \phi}{\partial x \partial y} \\ \sigma_{zz} &= \nu (\sigma_{xx} + \sigma_{yy}) \\ \sigma_{xz} &= \sigma_{yz} = 0 \end{aligned} \qquad (A.6) $$

### A.1.2. Displacement field determination

Knowing the stress field, we can compute the strain field with Hooke's law. To evaluate the displacements the relation $\epsilon_{ij} = \frac{1}{2}(u_{i,j} + u_{j,i})$ is used. Let us consider a point $x_0$ in the solid, where displacement and rotation are equal to zero. We can compute the displacement of any point in the solid by integrating the strains along the most proper way (see Fig. 21):

$$ u_i(\mathbf{x}) = \int_{x_0}^{x} U_{ij}(\mathbf{x}, \xi) d\xi_j \qquad (A.7) $$

with

$$ U_{ij}(\mathbf{x}, \xi) = \varepsilon(\xi) + (x_k - \xi_k) \left[ \frac{\partial \varepsilon_{ij}(\xi)}{\partial \xi_k} - \frac{\partial \varepsilon_{kj}(\xi)}{\partial \xi_l} \right] \qquad (A.8) $$

where $x_k$ are the position vector components of the point in which we are computing the displacements and $\xi_j$ are the position vector components of a point along the integration path.
<figure><img src="image_24.png" /><figcaption>Fig. 21. Displacement computation procedure: integration of the strains along the proper path.</figcaption></figure>

<figure><img src="image_25.png" /><figcaption>Fig. 22. Test case for Terzaghi: preforms compaction with free edge.</figcaption></figure>

## A.2. Application to Terzaghi's test case

The problem is summarized in Fig. 22 and boundary conditions of the problem are presented below:

$$ egin{cases} \sigma_{ij}h_j = 0 \Rightarrow \sigma_{xx} = \sigma_{yx} = 0 & \text{on } x = \pm L \text{ for any } y \\ \sigma_{xy} = 0, & \sigma_{yy} = -P_{\text{imp}} \text{ on } y = H \text{ for any } x \\ u_x(0, 0) = u_y(0, 0) = 0 & \end{cases} $$

Using Terzaghi's law ($\bar{\sigma} = \bar{\sigma}_{\bar{x}} - p_r I$) with $p_r(x,y) = V = -(ay+b)$ we can write the momentum conservation:

$$ \underline{\underline{\underline{\sigma^e}}} + \underline{\underline{\underline{f_p}}} = 0 $$

with

$$ \underline{\underline{\underline{f_p}}} = \begin{cases} 0 \\ -a \end{cases} $$

So we can find an Airy function that satisfies the boundary conditions whose form is

$$ \phi(x, y) = A_1 x^2 + A_2 xy + A_3 y^2 $$

This polynomial function verifies Eq. (A.4) for any value of $A_1$, $A_2$ and $A_3$. From Eq. (A.6), we obtain the following relations:

$$ \sigma_{xx} = 2A_3 + (ay+b) $$

$$ \sigma_{yy} = 2A_1 + (ay+b) $$

$$ \sigma_{xy} = -A_2 $$

Left and right frontiers being free, $\sigma_{xx} = \sigma_{yx} = 0$ on $x = \pm L$, implying $A_3 = A_2 = 0$ and leading to

$$ \sigma_{xx} = (ay+b) $$

$$ \sigma_{yy} = 2A_1 + (ay+b) $$

$$ \sigma_{xy} = \sigma_{yx} = 0 $$

The last constant $A_1$ can be determined with the boundary conditions on $y=H$ ($\sigma_{yy} = -P_{\text{imp}}$):

$$ \sigma_{yy} = 2A_1 + (aH+b) = -P_{\text{imp}} $$

$$ \Rightarrow A_1 = - \frac{P_{\text{imp}} + aH + b}{2} $$

Finally the stress field in the 2D plate is:

$$ egin{cases} \sigma_{xx} = ay + b \\ \sigma_{yy} = -P_{\text{imp}} + a(y - H) \\ \sigma_{xy} = \sigma_{yx} = 0 \end{cases} \qquad (A.9) $$

Then, using Hooke's law we can compute the strain field:

$$ egin{cases} \epsilon_{xx} = \frac{1+\nu}{E} [(1-2\nu)(ay+b)+\nu G] \\ \epsilon_{yy} = \frac{1+\nu}{E} [(1-2\nu)(ay+b)+(\nu-1)G] \\ \epsilon_{xy} = 0 \end{cases} \qquad (A.10) $$

with $G = aH+b+P_{\text{imp}}$.

Using Eq. (A.8), we have:

$$ U_{xx} = \frac{1+\nu}{E} [(1-2\nu)(a\eta+b)+\nu G] + (y-\eta) \frac{1+\nu}{E} [(1-2\nu)a] $$

$$ U_{yy} = \frac{1+\nu}{E} [(1-2\nu)(a\eta+b)+(\nu-1)G] $$

$$ U_{xy} = 0 $$

$$ U_{yx} = - \frac{1+\nu}{E} (x-\xi) [(1-2\nu)a] $$

We can finally determine the displacements $u_x$ and $u_y$ by resolving the following integral:

$$ egin{cases} u_x(x, y) = \int_{(0,0)}^{(x,y)} U_{xx} d\xi + U_{xy} d\eta \\ u_y(x, y) = \int_{(0,0)}^{(x,y)} U_{yx} d\xi + U_{yy} d\eta \end{cases} $$

So, the displacement field is:

$$ egin{cases} u_x(x, y) = \frac{1+\nu}{E} [(1-2\nu)(ay+b)x+\nu Gx] \\ u_y(x, y) = \frac{1+\nu}{E} [(1-2\nu)(\frac{a}{2}y+b)y+(\nu-1)Gy - \frac{1-2\nu}{2}ax^2] \end{cases} $$

with $G = aH+b+P_{\text{imp}}$.

## References

[1] P. Knupp, K. Salari, Verification of Computer Codes in Computational Science and Engineering, Chapman & Hall/CRC, 2003.

[2] P. Celle, S. Drapier, J.-M. Bergheau, Numerical modelling of liquid infusion into fibrous media undergoing compaction, Eur. J. Mech. A/Solids 27 (4) (2008) 647–661. http://dx.doi.org/10.1016/j.euromechsol.2007.11.002, URL http://linkinghub.elsevier.com/retrieve/pii/S0997753807001003.

[3] P. Celle, S. Drapier, J.-M. Bergheau, Numerical aspects of fluid infusion inside a compressible porous medium undergoing large strains, Eur. J. Comput. Mech. 17 (5-7) (2008) 819–827. http://dx.doi.org/10.3166/remn.17.819-827, URL http://remn.revuesonline.com/article.jsp?articleId=12264.

[4] S. Badia, A. Quaini, A. Quarteroni, Coupling Biot and Navier-Stokes equations for modelling fluid-poroelastic media interaction, J. Comput. Phys. 228 (21) (2009) 7986–8014. http://dx.doi.org/10.1016/j.jcp.2009.07.019, URL http://linkinghub.elsevier.com/retrieve/pii/S0021999109004069.

[5] M. Discacciati, Mathematical and numerical models for coupling surface and groundwater flows, Appl. Numer. Math. 43 (1–2) (2002) 57–74. http://dx.doi.org/10.1016/S0168-9274(02)00125-3, URL http://linkinghub.elsevier.com/retrieve/pii/S0168927402001253.

[6] M. Discacciati, Domain Decomposition Method for the Coupling of Surface and Groundwater Flows (Ph.D.), Ecole Polytechnique Fédéral de Lausane (EPFL), 2004.

[7] G. Pacquaut, J. Bruchon, N. Moulin, S. Drapier, Combining a level-set method and a mixed stabilized P1/P1 formulation for coupling Stokes-Darcy flows, Int. J. Numer. Methods Fluids 69 (2) (2012) 459–480. http://dx.doi.org/10.1002/fld.2569, URL http://doi.wiley.com/10.1002/fld.2569).
[8] A. Masud, A stabilized mixed finite element method for Darcy-Stokes flow, Int. J. Numer. Methods Fluids 54 (6-8) (2007) 665-681. http://dx.doi.org/10.1002/fld.1508, URL (http://doi.wiley.com/10.1002/fld.1508).

[9] M. Discacciati, Iterative methods for Stokes/Darcy coupling, in: Domain Decomposition Methods in Science and Engineering, vol. 40, 2004, pp. 563-570. doi: http://dx.doi.org/10.1007/3-540-26825-1_59. URL (http://www.springerlink.com/index/K22M2P65L1021317.pdf).

[10] K. Terzaghi, R.B. Peck, G. Mesri, Soil Mechanics in Engineering Practice, Wiley-IEEE, 1956. URL (http://books.google.com/books?id=bAwVv071FXoC&pgis=1) (http://eu.wiley.com/WileyCDA/WileyTitle/productCd-0471086584.html).

[11] M.A. Biot, General theory of three-dimensional consolidation, J. Appl. Phys. 12 (2) (1941) 155-164.

[12] M.A. Biot, General Solutions of the Equations of Elasticity and Consolidation for a Porous Material, J. Appl. Mech. (1956) 91-96.

[13] K.R. Rajagopal, A. Wineman, M. Gandhi, On boundary conditions for a certain class of problems in mixture theory, Int. J. Eng. Sci. 24 (8) (1986) 1453-1463. http://dx.doi.org/10.1016/0020-7225(86)90074-1, URL (http://linkinghub.elsevier.com/retrieve/pii/0020722586900741).

[14] T.G. Gutowski, T. Morigaki, The consolidation of laminate composites, J. Compos. Mater. 21 (2) (1987) 172-188. http://dx.doi.org/10.1177/002199838702100207, URL (http://jcm.sagepub.com/cgi/doi/10.1177/00219838702100207).

[15] T. Ouahbi, A. Saouab, J. Bréard, P. Ouagne, S. Chatel, Modelling of hydro-mechanical coupling in infusion processes, Composites Part A: Appl. Sci. Manuf. 38 (7) (2007) 1646-1654. http://dx.doi.org/10.1016/j.compositesa.2007.03.002, URL (http://linkinghub.elsevier.com/retrieve/pii/S1359835X07000346).

[16] A. Farina, L. Preziosi, Non-isothermal injection molding with resin cure and preform deformability, Composites Part A: Appl. Sci. Manuf. 31 (12) (2000) 1355-1372. http://dx.doi.org/10.1016/S1359-835X(00)00040-3, URL (http://linkinghub.elsevier.com/retrieve/pii/S1359835X09000403) (http://www.ncbi.nlm.nih.gov/pubmed/21435355).

[17] Q. Govignon, S. Bickerton, P.A. Kelly, Simulation of the reinforcement compaction and resin flow during the complete resin infusion process, Composites Part A: Appl. Sci. Manuf. 41 (1) (2010) 45-57. http://dx.doi.org/10.1016/j.compositesa.2009.07.007, URL (http://linkinghub.elsevier.com/retrieve/pii/S1359835X09002292).

[18] P. Šimáček, O. Eksik, D. Heider, J.W. Gillespie Jr., S.G. Advani, Post-filling flow in vacuum assisted resin transfer molding processes: theoretical analysis, Composites Part A: Appl. Sci. Manuf. 40 (6-7) (2009) 913-924. http://dx.doi.org/10.1016/j.compositesa.2009.04.018, URL (http://linkinghub.elsevier.com/retrieve/pii/S1359835X09001183) (http://linkinghub.elsevier.com/retrieve/pii/S1359835X11003277).

[19] R. Larsson, M. Rouhi, M. Wysocki, Free surface flow and preform deformation in composites manufacturing based on porous media theory, Eur. J. Mech. A/Solids 31 (1) (2012) 1-12. http://dx.doi.org/10.1016/j.euromechsol.2011.06.015, URL (http://linkinghub.elsevier.com/retrieve/pii/S0997753811001021).

[20] G.S. Beavers, D.D. Joseph, Boundary conditions at a naturally permeable wall, J. Fluid Mech. 30 (1) (1967) 197. http://dx.doi.org/10.1017/S0022112067001375, URL (http://www.journals.cambridge.org/abstract_S0022112067001375).

[21] D.N. Arnold, F. Brezzi, M. Fortin, A stable finite element for the Stokes equations, Calcolo 21 (4) (1984) 337-344. http://dx.doi.org/10.1007/BF02576171, URL (http://www.springerlink.com/index/10.1007/BF02576171).

[22] J.N. Reddy, D.K. Gartling, The Finite Element Method in Heat Transfer and Fluid Dynamics, third edition, CRC Press, 2010.

[23] C. Taylor, P. Hood, A numerical solution of the Navier-Stokes equations using the finite element technique, Comput. Fluids 1 (1) (1973) 73-100. http://dx.doi.org/10.1016/0045-7930(73)90027-3, URL (http://linkinghub.elsevier.com/retrieve/pii/0045793073900273).

[24] A. Hannukainen, M. Juntunen, R. Stenberg, Computations with finite element methods for the Brinkman problem, Comput. Geosci. 15 (1) (2010) 155-166. http://dx.doi.org/10.1007/s10596-010-9204-4, URL (http://www.springerlink.com/index/10.1007/s10596-010-9204-4).

[25] E. Feulvarch, J.-C. Roux, J.-M. Bergheau, A simple and robust moving mesh technique for the finite element simulation of Friction Stir Welding, J. Comput. Appl. Math. (2012) 1-9. http://dx.doi.org/10.1016/j.cam.2012.07.013, URL (http://dx.doi.org/10.1016/j.cam.2012.07.013)(http://linkinghub.elsevier.com/retrieve/pii/S0377042712003032).

[26] R. Codina, Stabilized finite element approximation of transient incompressible flows using orthogonal subscales, Comput. Methods Appl. Mech. Eng. 191 (39-40) (2002) 4295-4321. http://dx.doi.org/10.1016/S0045-7825(02)00337-7, URL (http://linkinghub.elsevier.com/retrieve/pii/S0045782502003377).

[27] L. Abouorm, N. Moulin, J. Bruchon, S. Drapier, Monolithic approach of Stokes-Darcy coupling for LCM process modelling, Key Eng. Mater. 554-557 (2013) 447-455. http://dx.doi.org/10.4028/www.scientific.net/KEM.554-557.447.

[28] L. Abouorm, R. Troian, S. Drapier, J. Bruchon, N. Moulin, Stokes/Darcy coupling in severe regimes using multiscale stabilization for mixed finite elements: monolithic approach versus decoupled approach, Eur. J. Comput. Mech. 23 (3-4) (2014) 1113-1137.

[29] L. Abouorm, R. Troian, S. Drapier, J. Bruchon, N. Moulin, Stokes/Darcy coupling in severe regimes using multiscale stabilization for mixed finite elements: monolithic approach versus decoupled approach, Eur. J. Comput. Mech., 23 (3-4), 2014, 1113-1137, http://dx.doi.org/10.1002/num.20349, URL (http://doi.wiley.com/10.1002/num.20349).

[30] P.C. Carman, Fluid flow through granular beds, Trans. Inst. Chem. Eng. 15 (1937) 150-166.

[31] D. Pino Muñoz, J. Bruchon, S. Drapier, F. Valdivieso, A finite element-based level set method for fluid-elastic solid interaction with surface tension, Int. J. Numer. Methods Eng. 93 (9) (2013) 919-941. http://dx.doi.org/10.1002/nme.4415, URL (http://doi.wiley.com/10.1002/nme.4415).

[32] F. Trochu, R. Gauvin, D.-M. Gao, Numerical analysis of the resin transfer molding process by the finite element method, Adv. Polym. Technol. 12 (4) (1993) 329-342. http://dx.doi.org/10.1002/adv.1993.060120401, URL (http://doi.wiley.com/10.1002/adv.1993.060120401).

[33] F. Trochu, R. Gauvin, J.-F. Boudreault, RTMFLOT-an integrated software environment for the computer simulation of the resin transfer molding process, J. Reinf. Plast. Compos. 13 (3) (1994) 262-270. http://dx.doi.org/10.1177/073168449401300306, URL (http://jrp.sagepub.com/cgi/doi/10.1177/073168449401300306).

[34] P. Wang, S. Drapier, J. Molimard, A. Vautrin, J.-C. Minni, Characterization of Liquid Resin Infusion (LRI) filling by fringe pattern projection and in situ thermocouples, Composites Part A: Appl. Sci. Manuf. 41 (1) (2010) 36-44. http://dx.doi.org/10.1016/j.compositesa.2009.09.007, URL (http://linkinghub.elsevier.com/retrieve/pii/S1359835X09002826).

[35] P. Wang, S. Drapier, J. Molimard, A. Vautrin, J.-C. Minni, Numerical and experimental analyses of resin infusion manufacturing processes of composite materials, J. Compos. Mater. 46 (13) (2011) 1617-1631. http://dx.doi.org/10.1177/0021998311421990, URL (http://jcm.sagepub.com/cgi/doi/10.1177/0021998311421990).

[36] P. Wang, J. Molimard, S. Drapier, A. Vautrin, J.-C. Minni, Monitoring the resin infusion manufacturing process under industrial environment using distributed sensors, J. Compos. Mater. 46 (6) (2011) 691-706. http://dx.doi.org/10.1177/0021998311410479, URL (http://jcm.sagepub.com/cgi/doi/10.1177/0021998311410479).