A Robust Monolithic Approach for Resin Infusion Based Process Modelling

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 m2. 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 represents 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 infusionbased 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 Ω ⊂ R m (m = 2 or m=3) as a bounded domain made up by two non overlapping subdomains Ω s and Ω d separated by a surface Γ = ∂Ω s ∩ ∂Ω d with Γ s,D ∩ Γ s,N = ∅, 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). Γ 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 whereε is the strain rate tensor defined byε(v s ) = 1 2 (∇v s + ∇v s T ), f s is the volume force, v 1 is the velocity prescribed on the boundary Γ s,D , n s is the unit vector normal to the boundary of Ω s , σ n,s is the normal stress prescribed on Γ s,N to be equal to −p ext,s n s and η 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 where K is the permeability tensor which can be reduced to a scalar in the isotropic case, p ext,d is a pressure to be prescribed on Γ d,N , n d is the outward unit vector normal to the boundary of Ω d and f d is the volume force.
Conditions (mass conservation and continuity of the stress) has to be considered on the interface Γ of normal n = n s = −n d ( Fig. 1 (b)). Continuity of the normal velocity. The mass conservation through the interface Γ is expressed by the continuity of the normal velocity field v across Γ v s · n s + v d · n d = 0 on Γ Continuity of the fluid normal stress. The continuity of normal stress over the interface Γ is expressed by n · σ s · n = n · σ d · n on Γ (4) Beaver-Joseph-Saffman condition. The Beaver Joseph Saffman condition allows the tangential velocity to be specified on the interface Γ [4], it writes where α is a dimensionless parameter, so-called slip coefficient, τ i are the tangential vectors on the interface and K ii are the permeabilities related to the tangential directions.

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 2 inner product in Ω d/s as < ·, · >. The following functional spaces are used as velocity, pressure and test function spaces

Material Forming ESAFORM 2014
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 Ω and a pressure field p on Ω such as v |Ω i = v i and p |Ω i = p i with i = s or i = d. The integrals over Ω s and Ω d must be re-defined on Ω. 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 ∫ Hence, the variational formulation of the Stokes-Darcy coupled problem consists in finding The bilinear form B c and the linear form L c are defined by

Finite element approximation with ASGS stabilization
We consider the bounded domain Ω ⊂ 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 Invoking this decomposition in the continuous problem (Eq. 10) for both solution and test functions, one gets the two-scales systems After approximating Eq. 12 with an algebraic formulation, by introducing an operator of projection onto V ′ × Q ′ , the approximated fields [v ′ , p ′ ] 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 with τ p,c , τ u,c are the stabilization parameters, that we compute in the case of a scalar permeability K as 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 τ p,c and τ u,c must be updated accordingly. In the orthotropic case, the permeability tensor K in the structural frame R = (0; X I , X II , X III ) is given by In considering that the Darcy's equations are independent in each space direction, the stabilization parameter τ u,c can be extended as a tensor of stabilization terms, the new stabilization parameters are defined by

Material Forming ESAFORM 2014
where 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 Γ separating the Stokes' and Darcy's domains and Γ f the moving flow front, are not described by a set of boundary elements. These interfaces go through the mesh elements. Consequently, two functions ϕ and ϕ f has to be introduced to depict these interfaces. ϕ, ϕ f are chosen as a signed distance function, respectively to Γ and Γ f . Γ and Γ f are then respectively described by the zero iso-surface of ϕ : Γ = {ϕ = 0} and ϕ f : When considering the discrete problem, ϕ and ϕ f are approximated by continuous and piecewise linear ϕ h and ϕ f h on Ω 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 < α η K (v h · τ ), (w h · τ ) > Γ 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 Γ 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 ϕ 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 η a ≪ η f , where η f is the fluid viscosity. Taking into account this moving flow front into Eq. 14 consists of simply replacing the constant resin viscosity η by the viscosity η ϕ , equal to η f in the filled domain (when ϕ f > 0) and equal to η a in the empty domain (when ϕ f < 0): Once the velocity is computed, ϕ f is transported by solving the following transport equation [11]: where ∂Ω − = {x ∈ ∂Ω; v · n < 0} is the inflow part of the boundary ∂Ω.
In order to avoid spurious oscillations when applying the continuous piecewise linear function ϕ f h 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, ϕ f is a signed distance function (i.e. ||∇ϕ f || = 1 which ensures the regularity of ϕ f ) but this property is generally not conserved during the computation (in fact, it depends on the velocity field 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., ||∇ϕ 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 η a = 1.10 −5 P a.s and η f = 3.10 −2 P a.s. The permeability tensor K is defined in the structural frame 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 P a is applied on the center of the computational domain and a pressure equals to 0P a 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. 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 R = (0; X I , X II ) of each porous domain of the part, by the permeability K I = 1.10 −13 m 2 K II = 1.10 −15 m 2 . The compu-

312
Material Forming ESAFORM 2014 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 with a Stokes-Darcy monolithic approach with a moving flow front and very low orthotropic permeabilities (10 −13 ∼ 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.

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.