# A robust monolithic approach for resin infusion based process modelling

Lara Abouorm, Maxime Blais, Nicolas Moulin, Julien Bruchon, Sylvain Drapier

► To cite this version:

Lara Abouorm, Maxime Blais, Nicolas Moulin, Julien Bruchon, Sylvain Drapier. A robust monolithic approach for resin infusion based process modelling. Key Engineering Materials, 2014, 611-612, pp.306-315. 10.4028/www.scientific.net/KEM.611-612.306. emse-01144951

HAL Id: emse-01144951

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

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
# A robust monolithic approach for resin infusion based process modelling

Lara ABOUORM¹,ᵃ, Maxime BLAIS¹,ᵇ, Nicolas MOULIN¹,ᶜ, Julien BRUCHON¹,ᵈ and Sylvain DRAPIER¹,ᵉ

LGF UMR CNRS 5307, École Nationale Supérieure des Mines de Saint-Étienne, 158 cours Fauriel, 42023 Saint-Étienne cedex 2

ᵃabouorm@emse.fr, ᵇmaxime.blais@emse.fr, ᶜnicolas.moulin@emse.fr, ᵈbruchon@emse.fr, ᶜdrapier@emse.fr

**Keywords:** Stokes-Darcy coupling, ASGS method, Finite strains, Level-set method, Composites.

**Abstract.** The aim of this work is to focus on the Stokes-Darcy coupled problem in order to propose a robust monolithic approach to simulate composite manufacturing process based on liquid resin infusion. The computational domain can be divided into two non-miscible sub-domains: a purely fluid domain and a porous medium. In the purely fluid domain, the fluid flows according to the Stokes' equations, while the fluid flows into the preforms according to the Darcy's equations. Specific conditions have to be considered on the fluid/porous medium interface. Under the effect of a mechanical pressure applied on the high deformable preform/resin stacking, the resin flows and infuses through the preform which permeability is very low, down to $10^{-15}m^2$. Flows are solved using finite element method stabilized with a sub-grid scale stabilization technique (ASGS). A special attention is paid to the interface conditions, namely normal stress and velocity continuity and tangential velocity constraint similar to a Beaver-Joseph-Saffman's condition. The originality of the model consists in using one single mesh to represent the Stokes and the Darcy sub-domains (monolithic approach). A level set context is used to represent Stokes-Darcy interface and to capture the moving flow front. This monolithic approach is now perfectly robust and leads to perform complex shapes for manufacturing process by resin infusion.

## Introduction

Resin infusion-based processes are used for manufacturing large structures of good quality (low void content) with high fiber volume fraction, especially for large pieces in aeronautics. Resin infusion-based processes consist of infusing liquid resin through the thickness of the reinforcement rather than in the plane (Fig. 1 (a)).

Under the effect of a mechanical pressure applied on the whole stacking, the resin flows into the preforms which can be seen as a porous medium subject to large deformations. However, despite numerous advantages, the control of these processes is difficult, especially for the most critical properties related to the final piece like its dimensions and its fibre volume fraction. To control this process, we developed a model based on the coupling between the resin flow within the porous domain (Darcy), and the purely fluid domain (Stokes).

The Stokes-Darcy coupled problem has been studied by many researchers in many fields of engineering using both a decoupled approach [7], [10] and monolithic approach as proposed by [12]. The decoupled strategy consists of using two different element spaces to solve the Stokes and Darcy equations, whereas the unified strategy consists in using the same finite element space. In the monolithic approach proposed by Pacquaut et al. [12] P1+/P1 elements in Stokes and HVM (Hughes Variational Multiscale) Method for stabilization in Darcy are used. Due to consistency errors and spurious oscillations that appear in this previous approach, we use a robust approach which yields improvements compared to the latest one. The robustness of the approach which is assessed in this paper, is ensured by using ASGS method (Algebraic Subgrid Scale) to stabilize velocity and pressure approximated by
linear and continuous elements in Stokes and Darcy domains. Signed functions are used to represent the Stokes-Darcy interface and to capture the moving flow front.

The paper is organized as follows. The first section presents the mathematical modelling for the Stokes-Darcy coupled problem. The next section introduces both the velocity-pressure mixed formulation for the Stokes-Darcy problem and the variational multiscale method used for stabilization. Then, a specific section is dedicated to the interface treatment and capturing. The last section shows numerical validation and results in severe regimes (low permeability, down to $10^{-15}m^2$, complex geometries) to illustrate the capability of modelling manufacturing processes by resin infusion.

## Mathematical model

Let us define $\Omega \subset \mathbb{R}^m$ ($m = 2$ or $m=3$) as a bounded domain made up by two non overlapping subdomains $\Omega_s$ and $\Omega_d$ separated by a surface $\Gamma = \partial\Omega_s \cap \partial\Omega_d$ with $\Gamma_{s,D} \cap \Gamma_{s,N} = \emptyset$, corresponding to two different kinds of boundary conditions: Dirichlet or Neumann conditions (Fig. 1 (a)). Index $s$ is used to denote everything that concerns the purely fluid domain (governed by the Stokes' equations) and index $d$ for the porous medium (governed by the Darcy's equations). $\Gamma$ is the interface between the Stokes and Darcy domains. In this section, we present the modelling of the resin flow, considered as an incompressible fluid, into the fluid distribution medium and into the preforms. Both Stokes' and Darcy's equations and their weak formulation are then introduced.

The Stokes' equations, which express momentum and mass balances when inertia is neglected, is written as: Find the velocity $v_s$ and the pressure $p_s$ fields defined by

$$ 
\begin{aligned}
-\nabla \cdot (2\eta\dot{\varepsilon}(\mathbf{v}_s)) + \nabla p_s &= \mathbf{f}_s && \text{in } \Omega_s \\
\nabla \cdot \mathbf{v}_s &= 0 && \text{in } \Omega_s \\
\mathbf{v}_s &= \mathbf{v}^1 && \text{on } \Gamma_{s,D} \\
\boldsymbol{\sigma}_{n,s} &= -p_{\text{ext},s}\mathbf{n}_s && \text{on } \Gamma_{s,N}
\end{aligned}
\quad (1) $$

where $\dot{\varepsilon}$ is the strain rate tensor defined by $\dot{\varepsilon}(\mathbf{v}_s) = \frac{1}{2}(\nabla\mathbf{v}_s + \nabla\mathbf{v}_s^T)$, $\mathbf{f}_s$ is the volume force, $\mathbf{v}^1$ is the velocity prescribed on the boundary $\Gamma_{s,D}$, $\mathbf{n}_s$ is the unit vector normal to the boundary of $\Omega_s$, $\boldsymbol{\sigma}_{n,s}$ is the normal stress prescribed on $\Gamma_{s,N}$ to be equal to $-p_{\text{ext},s}\mathbf{n}_s$ and $\eta$ is the fluid viscosity assumed to be constant (Newtonian fluid).

Respectively, the Darcy's equations are then expressed as: Find the velocity $v_d$ and pressure $p_d$ such as

$$ 
\begin{aligned}
\frac{\eta}{K} \mathbf{v}_d + \nabla p_d &= \mathbf{f}_d && \text{in } \Omega_d \\
\nabla \cdot \mathbf{v}_d &= 0 && \text{in } \Omega_d \\
\mathbf{v}_d \cdot \mathbf{n}_d &= g_d && \text{on } \Gamma_{d,D} \\
p_d &= p_{\text{ext},d} && \text{on } \Gamma_{d,N}
\end{aligned}
\quad (2) $$

where $K$ is the permeability tensor which can be reduced to a scalar in the isotropic case, $p_{\text{ext},d}$ is a pressure to be prescribed on $\Gamma_{d,N}$, $\mathbf{n}_d$ is the outward unit vector normal to the boundary of $\Omega_d$ and $\mathbf{f}_d$ is the volume force.

Conditions (mass conservation and continuity of the stress) has to be considered on the interface $\Gamma$ of normal $n = n_s = -n_d$ (Fig. 1 (b)).

**Continuity of the normal velocity.** The mass conservation through the interface $\Gamma$ is expressed by the continuity of the normal velocity field $\mathbf{v}$ across $\Gamma$

$$ \mathbf{v}_s \cdot \mathbf{n}_s + \mathbf{v}_d \cdot \mathbf{n}_d = 0 \quad \text{on } \Gamma \quad (3) $$

**Continuity of the fluid normal stress.** The continuity of normal stress over the interface $\Gamma$ is expressed by

$$ \mathbf{n} \cdot \boldsymbol{\sigma}_s \cdot \mathbf{n} = \mathbf{n} \cdot \boldsymbol{\sigma}_d \cdot \mathbf{n} \quad \text{on } \Gamma \quad (4) $$
**Beaver-Joseph-Saffman condition.** The Beaver Joseph Saffman condition allows the tangential velocity to be specified on the interface Γ [4], it writes

$$2\mathbf{n} \cdot \dot{\boldsymbol{\varepsilon}}(\mathbf{v}_s) \cdot \boldsymbol{\tau}_i = -\frac{\alpha}{\sqrt{K_{ii}}}(\mathbf{v}_s \cdot \boldsymbol{\tau}_i), \quad i = 1, 2 \qquad (5)$$

where α is a dimensionless parameter, so-called slip coefficient, τᵢ are the tangential vectors on the interface and Kᵢᵢ are the permeabilities related to the tangential directions.

<figure><img src="image_2.png" /><figcaption>Fig. 1: Domain decomposition into three zones for the modelling of resin infusion process: fluid distribution medium with resin, dry and wet preforms (a), 2D representation of the Stokes-Darcy domain (b).</figcaption></figure>

## Weak formulation

The weak formulation of the Stokes-Darcy coupled problem is obtained by summing up the Stokes and Darcy's weak formulations and by taking into consideration interface conditions described in the *Mathematical model* section. All details concerning both the weak formulation of Stokes' and Darcy's systems separately are given in [1, 2, 3]. For a sake of simplicity, we choose to write the L² inner product in Ω<sub>d/s</sub> as <·, ·>. The following functional spaces are used as velocity, pressure and test function spaces

$$ 
\begin{align*}
L^2(\Omega_i) &= \{q : \Omega_i \to \mathbb{R} \mid \int_{\Omega_i} q^2 d\Omega_i < \infty\} \\
H^1(\Omega_s)^m &= \{\mathbf{u} \in L^2(\Omega_s)^m \mid \nabla \mathbf{u} \in L^2(\Omega_s)^{m \times m}\} \\
H_{\Gamma_{s,D}}^1(\Omega_s)^m &= \{\mathbf{u} \in H^1(\Omega_i)^m \mid \mathbf{u} = \mathbf{u}^1 \text{ on } \Gamma_{s,D}\} \\
H_{\Gamma_{s,D}}^{1,t}(\Omega_s)^m &= \{\mathbf{u} \in H^1(\Omega_s)^m \mid \mathbf{u} = 0 \text{ on } \Gamma_{s,D}\} \\
H(\mathrm{div}, \Omega_d) &= \{\mathbf{u} \in L^2(\Omega_d)^m \mid \nabla \cdot \mathbf{u} \in L^2(\Omega_d)\} \\
H_{\Gamma_{d,D}}(\mathrm{div}, \Omega_d) &= \{\mathbf{u} \in H(\mathrm{div}, \Omega_d) \mid \mathbf{u} \cdot \mathbf{n} = u^2 \text{ on } \Gamma_{d,D}\} \\
H_{\Gamma_{d,D}}^t(\mathrm{div}, \Omega_d) &= \{\mathbf{u} \in H(\mathrm{div}, \Omega_d) \mid \mathbf{u} \cdot \mathbf{n} = 0 \text{ on } \Gamma_{d,D}\}
\end{align*}
\qquad (6) $$
with $i = s$ or $i = d$ and $m$ is the dimension equal to 2 or 3.

The mixed formulation of the Stokes-Darcy problem is established by considering a velocity $v$ on $\Omega$ and a pressure field $p$ on $\Omega$ such as $v|_{\Omega_i} = v_i$ and $p|_{\Omega_i} = p_i$ with $i = s$ or $i = d$. The integrals over $\Omega_s$ and $\Omega_d$ must be re-defined on $\Omega$. This is achieved by introducing the Heaviside function $H_i$ equal to 1 in the domain $i$ and vanishing elsewhere. These functions allow us to write

$$ \int_{\Omega_i} (...) d\Omega = \int_{\Omega} (...) H_i d\Omega \quad \text{with } i = s \text{ or } i = d \qquad (7) $$

Hence, the variational formulation of the Stokes-Darcy coupled problem consists in finding $[v, p] \in \mathcal{V} \times \mathcal{Q}$ such that

$$ B_c[(v, p), (w, q)] = L_c([w, q], \quad \forall [w, q] \text{ defined in } \mathcal{V}_0 \times \mathcal{Q}) \qquad (8) $$

with

$$ \begin{aligned} \mathcal{V} &= H_{\Gamma_{s,D}}^{1} \times H_{\Gamma_{d,D}}(\mathrm{div}, \Omega_d) \\ \mathcal{V}_0 &= H_{\Gamma_{s,D}}^{1,t} \times H_{\Gamma_{d,D}}^{t}(\mathrm{div}, \Omega_d) \\ \mathcal{Q} &= L^2(\Omega_s) \times L^2(\Omega_d) \end{aligned} \qquad (9) $$

The bilinear form $B_c$ and the linear form $L_c$ are defined by

$$ \begin{aligned} B_c([v, q], [w, q]) &= \langle 2\eta\dot{\epsilon}(v) : \dot{\epsilon}(w) H_s \rangle_{\Omega} + \left\langle \frac{\mu}{K} v, w H_d \right\rangle_{\Omega} - \langle p, \nabla \cdot w \rangle_{\Omega} \\ &\quad + \langle q, \nabla \cdot v \rangle_{\Omega} + \left\langle \frac{\alpha\eta}{\sqrt{K}} v, w \right\rangle_{\Gamma} \end{aligned} \qquad (10) $$

$$ L_c([v, p], [w, q]) = \langle \mathbf{f}, \mathbf{w} \rangle_{\Omega} + \langle h, q \rangle_{\Omega} + \langle p_{\mathrm{ext},s} \mathbf{n}, \mathbf{w} \rangle_{\Gamma_{s,N}} + \langle p_{\mathrm{ext},d} \mathbf{n}, \mathbf{w} \rangle_{\Gamma_{d,N}} $$

($\mathbf{f}, h$) are defined by ($\mathbf{f}_i, h_i$) in $\Omega_i$.

## Finite element approximation with ASGS stabilization

We consider the bounded domain $\Omega \subset \mathbb{R}^m$ discretized into $n_{el}$ non-overlapping elements. This one single unstructured mesh is made up of triangles if $m = 2$ and of tetrahedrons if $m = 3$. The Galerkin approximation of both the Stokes and the Darcy problems requires the use of velocity-pressure interpolation that satisfy the adequate inf-sup condition. Different interpolation pairs are known to satisfy this condition for each problem independently, but the key issue is to find interpolations that satisfy both at the same time. In this paper, we choose the use of stabilized finite element methods. The philosophy of the stabilized methods is to strengthen classical variational formulations so that discrete approximation which would otherwise be unstable becomes stable and convergent. One of the stabilized finite element method approximation for Stokes-Darcy problem is VMS (Variational Multiscale Method) [8, 9]. An important feature is that the finite element space are conforming $\mathcal{V}_h \times \mathcal{Q}_h \subset \mathcal{V} \times \mathcal{Q}$ and that the finite element spaces in Stokes' and Darcy's domains are the same because we use a monolithic approach. Let $\mathcal{V}_h$ and $\mathcal{Q}_h$ be the finite element spaces of the piecewise linear and continuous functions, which contains the solutions $\mathbf{v}_{h,i}$ and $p_{h,i}$. The basic idea of this method is to approximate the effect of the component of the continuous solution which cannot be captured by the finite element solution. It consists in splitting the continuous solution for velocity and pressure into two components, one coarse corresponding to the finite element scale $[\mathbf{v}_{h,i}, p_{h,i}]$, and a finer component corresponding to lower scale $[\mathbf{v}'_i, p'_i]$ for resolutions. The velocity is approximated as $\mathbf{v}_i = \mathbf{v}_{h,i} + \mathbf{v}'_i$ and the pressure field is approximated as $p_i = p_{h,i} + p'_i$. We consider a subgrid space $\mathcal{V} \times \mathcal{Q} = (\mathcal{V}_h \times \mathcal{Q}_h) \oplus (\mathcal{V}' \times \mathcal{Q}')$.
Invoking this decomposition in the continuous problem (Eq. 10) for both solution and test functions, one gets the two-scales systems

$$B_c([\mathbf{v}_h, p_h]), [\mathbf{w}_h, q_h]) + B_c([\mathbf{w}^\prime, p^\prime], [\mathbf{w}_h, q_h]) = L_c([\mathbf{w}_h, q_h]) \quad (11)$$

$$B_c([\mathbf{v}_h, p_h], [\mathbf{w}^\prime, p^\prime]) + B_c([\mathbf{v}^\prime, p^\prime], [\mathbf{w}^\prime, q^\prime]) = L_c([\mathbf{w}^\prime, q^\prime]) \quad (12)$$

for all $[\mathbf{w}_h, q_h] \in \mathcal{V}_h \times \mathcal{Q}_h$ and $[\mathbf{w}^\prime, q^\prime] \in \mathcal{V}' \times \mathcal{Q}'$. After approximating Eq. 12 with an algebraic formulation, by introducing an operator of projection onto $\mathcal{V}' \times \mathcal{Q}'$, the approximated fields $[\mathbf{v}^\prime, p^\prime]$ are taken into account in the finite element problem (Eq. (11) [1].

For Stokes and Darcy flows coupled through their interfaces, the stabilized problem with ASGS can be written as follows

$$B_{c,\text{stable}}([\mathbf{v}_h, p_h][\mathbf{w}_h, q_h]) = L_{c,\text{stable}}([\mathbf{w}_h, q_h]) \quad (13)$$

with

$$\begin{aligned} B_{c,\text{stable}}([\mathbf{v}_h, p_h], [\mathbf{w}_h, q_h]) = & 2\eta \langle H_s \dot{\mathbf{\epsilon}}(\mathbf{v}_h) : \dot{\mathbf{\epsilon}}(\mathbf{w}_h) \rangle_{\Omega} + \frac{\eta}{K} \langle H_d \mathbf{w}_h, \mathbf{v}_h \rangle_{\Omega} \\ & - \langle \nabla \cdot \mathbf{w}_h, p_h \rangle + \langle \nabla \cdot \mathbf{v}_h, q_h \rangle_{\Omega} + \tau_{p,c} \langle \nabla \cdot \mathbf{v}_h, \nabla \cdot \mathbf{w}_h \rangle_{\Omega} \\ & + \langle \alpha \frac{\eta}{\sqrt{K}} (\mathbf{v}_h \cdot \boldsymbol{\tau}), (\mathbf{w}_h \cdot \boldsymbol{\tau}) \rangle_{\Gamma} \\ & + \tau_{u,c} \langle \frac{\eta}{K} \mathbf{v}_h + \nabla p_h, -\frac{\eta}{K} \mathbf{w}_h - \nabla q_h \rangle_{\Omega} \end{aligned} \quad (14)$$

$$\begin{aligned} L_{c,\text{stable}}([\mathbf{w}_h, q_h]) = & \langle \mathbf{f}, \mathbf{w}_h \rangle_{\Omega} + \langle h, q_h \rangle_{\Omega} + \langle p_{\text{ext},s} \mathbf{n}, \mathbf{w}_h \rangle_{\Gamma_{s,N}} \\ & + \langle p_{\text{ext},d} \mathbf{n}, \mathbf{w}_h \rangle_{\Gamma_{s,D}} - \tau_{p,c} \langle h, \nabla \cdot \mathbf{w}_h \rangle_{\Omega} \\ & + \tau_{u,c} \langle \mathbf{f}, -\frac{\eta}{K} \mathbf{w}_h - \nabla q_h \rangle_{\Omega} \end{aligned}$$

$\tau_{p,c}$, $\tau_{u,c}$ are the stabilization parameters, that we compute in the case of a scalar permeability $K$ as

$$\begin{aligned} \tau_{p,c} &= c_p \frac{\eta}{K} l_p^2 H_d + c_1 \eta H_s \\ \tau_{u,c} &= \left( c_1 \eta H_s + c_u \frac{\eta}{K} l_u^2 H_d \right)^{-1} h_k^2 \end{aligned} \quad (15)$$

with $c_p$ and $c_u$ some algorithmic constants. $l_u$ and $l_p$ are length scales which we choose to take $(L_0 h_k)^{1/2}$, $L_0$ is a characteristic length of the domain and $h_k$ is the element size.

**Orthotropic permeabilities.** In the Liquid Resin Infusion processes, the permeability is orthotropic, for that we have to extend our stabilization method. Consequently, the stabilization parameters $\tau_{p,c}$ and $\tau_{u,c}$ must be updated accordingly. In the orthotropic case, the permeability tensor $K$ in the structural frame $\mathcal{R} = (0; X_I, X_{II}, X_{III})$ is given by

$$K = \begin{bmatrix} K_I & 0 & 0 \\ 0 & K_{II} & 0 \\ 0 & 0 & K_{III} \end{bmatrix}_{(0;X_I,X_{II},X_{III})} \quad (16)$$

In considering that the Darcy's equations are independent in each space direction, the stabilization parameter $\tau_{u,c}$ can be extended as a tensor of stabilization terms, the new stabilization parameters are defined by

$$\begin{aligned} \tau_{p,c} &= c_p \frac{\eta}{K} l_p^2 H_d + c_1 \eta H_s \\ \tau_{u,c} &= \left( \frac{1}{c_1 \eta} H_s \mathbf{I} + \frac{1}{c_u \eta l_u^2} H_d K \right) h_k^2 \end{aligned} \quad (17)$$
where $\overline{K}$ is a "representative" permeability defined here, such as
$$[max(K_I, K_{II}, K_{III}) + min(K_I, K_{II}, K_{III})]/2.$$

**Interface capturing.** In our monolithic approach, both interfaces $\Gamma$ separating the Stokes' and Darcy's domains and $\Gamma_f$ the moving flow front, are not described by a set of boundary elements. These interfaces go through the mesh elements. Consequently, two functions $\phi$ and $\phi_f$ has to be introduced to depict these interfaces. $\phi, \phi_f$ are chosen as a signed distance function, respectively to $\Gamma$ and $\Gamma_f$. $\Gamma$ and $\Gamma_f$ are then respectively described by the zero iso-surface of $\phi : \Gamma = \{\phi = 0\}$ and $\phi_f : \Gamma_f = \{\phi_f = 0\}$. When considering the discrete problem, $\phi$ and $\phi_f$ are approximated by continuous and piecewise linear $\phi_h$ and $\phi_{fh}$ on $\Omega_h$.

The drawback of the monolithic approach is that the surface integrals are not easily and directly computable. There is two ways to compute the surface integral $\langle \alpha \frac{\eta}{K} (\mathbf{v}_h \cdot \boldsymbol{\tau}), (\mathbf{w}_h \cdot \boldsymbol{\tau}) \rangle_{\Gamma}$ involved in Eq. 14. Either the surface integral is turned into a volume integral by introducing a Dirac delta function, or it is computed exactly by rebuilding a piecewise linear interface [13]. The exact computation of the surface integrals shows more accurate results than Dirac approximation [3]. Moreover, the Stokes' ans Darcy's contributions are averaged over the elements cut by the interface $\Gamma$ thanks to the Heaviside functions $H_s$ and $H_d$ introduced in Eq. 7, 10.

The specificity of our approach is the introduction of a moving flow front, depicted above by the level-set function $\phi_f$, separating the part of the domain which already filled with the resin from the part which is not filled yet. Because of our monolithic approach, the velocity and pressure fields have to be defined on the whole computational domain assuming that the empty part is filled with a Newtonian incompressible fluid (referred to the "air") having a very low viscosity $\eta_a \ll \eta_f$, where $\eta_f$ is the fluid viscosity. Taking into account this moving flow front into Eq. 14 consists of simply replacing the constant resin viscosity $\eta$ by the viscosity $\eta_\phi$, equal to $\eta_f$ in the filled domain (when $\phi_f > 0$) and equal to $\eta_a$ in the empty domain (when $\phi_f < 0$):

$$\eta_{\phi} = H(\phi)\eta_{f} + (1 - H(\phi))\eta_{a} \quad (18)$$

Once the velocity is computed, $\phi_f$ is transported by solving the following transport equation [11]:

$$\begin{aligned}
\frac{\partial \phi_f}{\partial t} + \mathbf{v} \cdot \nabla \phi_f &= 0 \quad \forall (x, t) \in \Omega \times (0, T) \\
\phi_f(x, t = 0) &= \phi_0(x) \quad \forall x \in \Omega \\
\phi_f(x, t) &= g(x, t) \quad \forall x \in \partial\Omega^-, \forall t \in (0, T)
\end{aligned} \quad (19)$$

where $\partial\Omega^- = \{x \in \partial\Omega; \mathbf{v} \cdot \mathbf{n} < 0\}$ is the inflow part of the boundary $\partial\Omega$.

In order to avoid spurious oscillations when applying the continuous piecewise linear function $\phi_{fh}$ approximation, the advection equation (Eq. 19) is solved by using a classical Streamline-Upwind Petrov-Galerkin method and by using an implicit Euler or Crank-Nicolson schemes for the time discretization. As stated above, $\phi_f$ is a signed distance function (i.e. $||\nabla\phi_f|| = 1$ which ensures the regularity of $\phi_f$) but this property is generally not conserved during the computation (in fact, it depends on the velocity field $\mathbf{v}$). Hence, a reinitialization step may be necessary to ensure that the solution is not deteriorated and to recover the property of the signed distance function (i.e., $||\nabla\phi_f|| = 1$) without disturbing the position of the flow front. Further details about the reinitialization can be found in [6].

## Numerical tests and simulations of resin infusion based process

**Bi-dimensional flow.** The first simulation consists in a bi-dimensional flow with an orthotropic permeability of the porous medium. The aim of the numerical test is to validate the Stokes-Darcy coupling presented in this paper for a transient flow with a very low permability. The simulation has been carried out by using the finite element software Z-set [5]. The computational domain is divided into a
purely fluid domain elliptical in shape (Stokes' domain at the center) and a porous medium (Fig. 2). The Stokes-Darcy interface and the resin flow front are depicted by both the zero iso-surface value coloured in white and the zero iso-surface value coloured in black on the Fig. 2 (a). The physical parameters for this simulation are $\eta_a = 1.10^{-5} Pa.s$ and $\eta_f = 3.10^{-2} Pa.s$. The permeability tensor $K$ is defined in the structural frame $\mathcal{R} = (0; X_I, X_{II})$ by the permeability $K_I = 5.10^{-15} m^2$ and $K_{II} = 1.10^{-15} m^2$. A pressure equals to $1.10^5 Pa$ is applied on the center of the computational domain and a pressure equals to $0 Pa$ is applied on the boundary of the whole domain. This case has been investigated in [14]. An analytical solution allows to describe the position of the flow front during the simulation. The Fig. 2 (a) shows the pressure iso-values provided by the simulation. A pressure gradient appears in the wet preform zone. Fig. 2 (b) compares both the simulation results to the analytical solution in the two main directions $X_I$ and $X_{II}$. One can verify the good correlation between the analytical and the numerical solutions.

<figure><img src="image_3.png" /><figcaption>Fig. 2: Pressure field in a bi-dimensional Stokes-Darcy flow with an orthotropic permeability (a), evolution of the flow front (b).</figcaption></figure>

**Numerical simulation of the manufacturing process by resin infusion with moving interface.**

The example proposed here represents the manufacturing of a simple piece in which resin is injected and then infuses through the preforms. The computational domain is depicted on the Fig. 3 (a) with an injection channel through which the resin is injected. The boundary conditions prescribed for this simulation are a normal stress on the inflow part of the domain and a zero normal velocity on the other edges of the domain. The permeability tensor $K$ is defined in the structural frame $\mathcal{R} = (0; X_I, X_{II})$ of each porous domain of the part, by the permeability $K_I = 1.10^{-13} m^2$ and $K_{II} = 1.10^{-15} m^2$. The compu-
tational domain is discretized with a fixed mesh of 9,592 triangular elements corresponding to 5,031 nodes. Fig. 3 shows the evolution of the flow front showing that the orthotropic permeability (Fig. 3 (b),(d),(f)) has a strong effect on the filling of the porous medium compared with the isotropic case (Fig. 3 (c),(c),(g)). The type of permeability also affects the pressure fields in the porous domain and then modify the flow in the tilted part. Moreover, the results seem to show that when the front reaches the top of the horizontal part of the domain, air is entrapped air, causing defects. These defects could be reduced by placing properly the vents in an optimization process. The robustness of our approach to simulate realistic geometries in the context of the manufacturing processes of composite materials has been demonstrated. These simulations were conducted with a Stokes-Darcy monolithic approach with a moving flow front and very low orthotropic permeabilities ($10^{-13} \sim 10^{-15}m^2$).

All the developments have been extended to 3D case and a numerical simulation of the manufacturing process by resin infusion has been performed but not included in this paper.

<figure><img src="image_4.png" /><figcaption>Fig. 3: Geometry and boundary conditions used for the simulation of the manufacturing process by resin infusion (a) Evolution of the resin flow front and pressure fields for three different times for orthotropic permeability $K_I = 1.10^{-13}m^2$ and $K_{II} = 1.10^{-15}m^2$ (b), (d), (f) and for an isotropic permeability $K_I = K_{II} = 1.10^{-15}m^2$ (c), (e), (g).</figcaption></figure>
## Conclusion

A unified strategy has been developed to solve the Stokes-Darcy coupled problem. A stabilized finite element method has been proposed to stabilize Stokes-Darcy coupled problem in the case where the Brezzi-Babuska stability condition is not fulfilled. This stabilization method is based on a variational multiscale technique called ASGS method. Convergence of this method was validated in [1] and compared with another approach [3] showing the accuracy of the results.

A numerical simulation of a bi-directional flow were presented to validate the Stokes-Darcy coupling with orthotropic permeabilities. At the end of the paper, 2D simulations of the manufacturing process by resin infusion were presented. In these simulations, two signed distance functions were used, one to represent the interface between the purely fluid domain and the porous medium and a second to capture the flow front. It was shown that our monolithic approach is relevant to simulate the resin infusion processes in severe regimes (very low orthotropic permeabilities...).

Regarding the prospects for the future work, we shall take into account the deformation of preforms by using an Updated Lagrangian scheme relying on displacement-based finite element. A special care will be paid to the interaction of the preform deformation and the resin infusion. Indeed, the resin pressure will be modified by the permeability change induced by preform compaction, and conversely the mechanical response of the porous medium will be represented via a Terzaghi's model modified according to the current resin pressure.

## Literature References

### References

[1] L. Abouorm: *VMS (Variational Multi Scale) stabilised Stokes-Darcy coupled flows in porous media undergoing finite deformations: application to infusion-based composite processing*, Ph.D. Thesis, Ecole Nationnale Supérieure des Mines de Saint-Etienne (2013).

[2] L. Abouorm, N. Moulin, J. Bruchon and S. Drapier: Key Engineering Materials Vol. 554-557 (2013) p. 447-455.

[3] L. Abouorm, R. Troian, S. Drapier, J. Bruchon and N. Moulin: European Journal of Computational Mechanics (2014) *accepted*.

[4] G. Beavers and D. Joseph: Journal of Fluid Mechanics Vol. 30 (1967) p. 197-207.

[5] J. Besson and R. Foerch: Computer Methods in Applied Mechanics and Engineering Vol. 142 (1997) p. 691-706.

[6] J. Bruchon, H. Digonnet and T. Coupez: International Journal for Numerical Methods in Engineering Vol. 78 (2009) p. 980-1008.

[7] P. Celle, S. Drapier, and J-M. Bergheau: European Journal of Computational Mechanics Vol. 17 (2008) p. 819-827.

[8] S. Badia and R. Codina: SIAM journal on Numerical Analysis, Vol. 47 (2008), p. 1971-2000.

[9] S. Badia and R. Codina: Computer Methods in Applied Mechanics and Engineering, Vol. 199 (2010), p. 654-667.

[10] M. Discacciati and A. Quarteroni, Revista Metematica Complutense Vol. 22 (2009) p. 315-426.

[11] S. Osher and J.A. Sethian, Journal of Computational Physics Vol. 79 (1988) p. 12-49.
[12] G. Pacquaut, J. Bruchon, N. Moulin, and S. Drapier: International Journal for Numerical Methods in Fluids, Vol. 69 (2012), p. 459-480.

[13] D. Pino Munoz, J. Bruchon, S. Drapier and F. Valdivieso, International Journal for Numerical Methods in Engineering Vol. 93 (2013) p. 919–941.

[14] J.R. Weitzenböck, R.A. Shenoi and P.A. Wilson: Composites Part A: Applied Science and Manufacturing Vol. 30 (1999) p. 781-796.