Sébastien Comas-Cardona

Paul H L Groenenboom

Christophe Binétruy
email: binetruy@ensm-douai.fr

Patricia Krawczak

Simulation of liquid composite molding processes using a generic mixed FE-SPH method

Keywords: Smoothed Particle Hydrodynamics, finite elements, liquid composite molding, explicit numerical schemes. MOTS-CLES : Smoothed Particle Hydrodynamics, éléments finis, liquid composite molding, schémas d'intégration explicite

Composite manufacturing processes involve multi-scale phenomena. Equations governing mechanics, heat transfer and fluid flow dynamics can be derived in the scale of the problem to solve. This article focuses on addressing fluid flow simulation needs for Liquid Composite Molding processes using a generic mixed FE-SPH method. The SPH method is Lagrangian and models fluids as particles. The method has been proven to be suitable to simulate fluid flows. Special solutions have been developed for flow through porous media and non-Newtonian fluid flow. Numerical schemes for such solutions are also given. The implementation of the SPH algorithm within a structural finite element software facilitates simulation of fluid-structure for LCM processes. Several applications are presented and discussed. RÉSUMÉ. La fabrication de matériaux composites met en jeu des phénomènes multiéchelles. Les équations qui gouvernent la dynamique des fluides, la mécanique et le transfert thermique peuvent être développées à l'échelle du phénomène à observer. Cet article se focalise sur les besoins de simulation des moulages composites en utilisant une méthode générique FE-SPH. La méthode SPH est lagrangienne et modélise le fluide à l'aide de particules. Il a été prouvé que cette méthode est appropriée pour simuler et modéliser l'écoulement de fluides. Des solutions spéciales ont été développées pour l'écoulement dans un milieu poreux et pour un fluide non newtonien. Les schémas d'intégration numérique sont également présentés. Les particules représentant le fluide peuvent interagir avec des éléments finis, ce qui permet de simuler des procédés de fabrication requérant des physiques couplées. Plusieurs applications sont présentées et discutées.

Introduction

Composite materials are covering a wide range of applications for aeronautics, automotive and marine structures. Such composite materials are usually made of fibers enclosed in a thermosetting or thermoplastic matrix. Liquid Composite Molding (LCM) processes, such as compression, Resin Transfer Molding (RTM) or Vacuum Assisted Resin Transfer Molding (VARTM), have become of increasing interest to manufacture composite parts. Reproducibility, quality, performances, environment friendliness and cost reduction are goals a composite manufacturing company seeks. Those goals will inevitably require the use of reliable and powerful simulation tools that enhances predictions and reduces trial and error procedures.

It is well established that micro-, meso-and macroscale media are present in composite manufacturing processes. For instance, for continuous fiber reinforced composites, the microscale corresponds to the inner tow region, whereas the inter tow region would be mesoscale, and the macroscale represents the overall composite structure. Equations governing mechanics, heat transfer or fluid flow can be derived for the three scales depending on the problem to be solved or the phenomenon looked at. More precisely for fluid flow simulation, the Navier-Stokes equations, averaged Navier-Stokes for porous media [START_REF] Tucker | Governing equations for flow and heat transfer in stationary fibre beds[END_REF], or Darcy's law will be generally used. The equations are usually discretized in an Eulerian reference frame using Finite Difference Method (FDM) or Finite Element Method (FEM) [START_REF] Trochu | Numerical analysis of the Resin Transfer Molding process by the Finite Element Method[END_REF]. However, when the geometry of the flow domain changes, because of interfaces with moving structures (such as compression, injectioncompression or VARTM) or the presence of a free surface, special solutions are required. Also when coupling is present, e.g. hydro-mechanical coupling, [START_REF] Han | Fibre mat deformation in liquid composite molding, II: Modelling[END_REF] showed that mixed Control Volume/FDM can be used. (Also [START_REF] Kang | Analysis of vacuum bag resin transfer molding process[END_REF][START_REF] Acheson | The implications of fiber compaction and saturation on fully coupled VARTM simulation[END_REF] coupled resin flow, fiber preform mechanics and tow saturation during VARTM process through Control Volume/FEM or FDM algorithms. Those Eulerian numerical schemes are limited to small deformations for instance or do not allow to mix domains with porous and non-porous media. Lately, Lagrangian numerical schemes have been of interest with [START_REF] Farina | Flow in deformable porous media: Modelling and simulations of compression moulding processes[END_REF][START_REF] Belov | Modelling of permeability of textile reinforcements: lattice Boltzmann method[END_REF][START_REF] Sawley | Modelling of Flow in Porous Media and Resin Transfer Moulding using Smoothed Particle Hydrodynamics[END_REF]. The later qualitatively reproduced the flow front shape of an isothermal RTM filling using the Smoothed Particle Hydrodynamics (SPH) method for mesoscale simulations. This paper particularly focuses on the use of the generic mixed FE-SPH method to address the simulation needs of LCM processes. First, fundamentals on fluid flow dynamics are given for porous and non-porous regions. Then, a detailed description of the mixed FE-SPH method is presented. Finally, several examples of applications for LCM simulations will be proposed and discussed.

Fluid flow dynamics

In the following sections, a bold lower case letter represents a vector and a bold upper case letter denotes a tensor. [START_REF] Tucker | Governing equations for flow and heat transfer in stationary fibre beds[END_REF] and [START_REF] Pillai | Governing equations for unsaturated flow through woven fibre mats, Part 1. Isothermal flows[END_REF] derived the conservation of mass and momentum equations for flow in porous media from averaging the conservation equations for a two-phase medium in which one of the materials is fixed in space. In that case, the volume fraction available for the fluid remains fixed and is also called porosity. Considering laminar flow, the suitable forms of the equations to be solved (conservation of mass and momentum) are respectively:

Fluid flow through porous media

0 u t [1] u K g uu u 1 2 p t [2]
where is the fluid volume fraction or porosity, u the fluid velocity, the viscosity, the viscous stress tensor, p the pressure, g the gravitational acceleration, K the permeability (tensor or scalar) and the apparent density where is the fluid density. For Equation 2, it has been assumed that the flow velocity in the porous region is small enough to assume that the viscous drag depends linearly on the velocity.

A scaling analysis can show that the viscous term (i.e., the second term of the RHS) can be, in most cases, neglected for LCM. With this assumption and using Equation 1, Equation 2 becomes:

u K g u u u 1 2 p t [3]
Most of the time, the computational domain can be divided in porous (e.g., fiber reinforcements) and non-porous (e.g., channels or tubing) regions. In order to simulate the flow for all regions, Equation 3 provides the basis of the solution with the viscous drag term only being used within the porous domains. As a limiting case, within the porous region (i.e., where the flow is dominated by the viscous drag) the convective and inertia terms (LHS of Equation 3) becomes negligible and therefore Darcy's law for the transport (or filtration) velocity u v may be verified.

Non-Newtonian fluid flow

For the flow of a liquid around individual and immersed solid objects, the following conservation of mass and Navier-Stokes equations should be solved:

0 u t [4] g u u u p t [5]
For composite molding conditions, turbulence is irrelevant, but non-Newtonian laminar flow should be considered. The components of the viscous stress tensor may be written as:

ij k k i j j i ij x u x u x u 3 2 [6]
where the indices i, j, and k refer to the three Cartesian coordinates and ij is the Kronecker delta. Note that the divergence term on the RHS of Equation 6 is maintained since limited compressibility will be assumed in the SPH formulation (Equation 15), but maintain as low as possible (see Sections 2.3 and 3.3). Among several non-Newtonian rheological models, the Carreau one will be chosen because it can model fairly well most polymeric fluid behaviors [START_REF] Mccabe | Predicting the Newtonian viscosity of complex fluids from high strain rate molecular simulations[END_REF][START_REF] Achilleos | On numerical simulations of polymer extrusion instabilities[END_REF]:

1 0 1 ) ( E [7]
where E is the total strain rate and and are all parameters of the fluid to be fitted to experimental data. The asymptotic behavior at large E is Newtonian with a viscosity of , for small E the effective viscosity tends to , and at intermediate regime it displays shear-thinning behavior rather like a power-law fluid. Since this model does not contain a yield stress, nor infinities at any limit, no numerical instabilities are expected. Note that if } (TXDWLRQ 7 models a Newtonian fluid.

Compressible fluid formulation

During composite molding, the flow of resin may be considered as incompressible. In that case Equation 1 and 4 would reduce to .u = 0 and the pressure should be obtained from a Poisson type of equation. Such a system of equations can be solved efficiently by an implicit method which allows for a relatively large time step. In the present study, however, the solution time step needs to be small. The first reason is that the motion of the flow front requires an accurate solution in time and space when the resin passes the border of the porous region. A second reason is that the weak coupling of the flow with moving structures for instance, requires an explicit solution with a sufficiently small time step. If the time step has to be small, an explicit solution in combination with compressibility of the liquid is more suitable.

Smoothed particle hydrodynamics fundamentals

The fully Lagrangian and meshless SPH method has been originally developed by [START_REF] Gingold | Smoothed particle hydrodynamics: Theory and application to non-spherical stars[END_REF] and [START_REF] Lucy | A Numerical Approach to the Testing of Fusion Process[END_REF] to solve astrophysical problems. Later, [START_REF] Monaghan | Shock Simulation by the Particle Method SPH[END_REF][START_REF] Monaghan | Simulating free surface flows with SPH[END_REF] applied the SPH method to the flow of compressible fluids and also showed that when pressure wave propagation is not of prime interest (as in LCM simulations), the use of an artificial equation of state is quite suitable to model incompressible fluid flows. SPH solutions for incompressible flow have also been developed, but such solutions have not been pursued for LCM simulations.

Integral interpolants and kernel functions

Mathematically, an integral interpolant of an arbitrary function (r) can be defined as:

' d ) , ' ( ) ' ( ) ( r r r r r h W [8]
where W is a kernel function having the following properties:

when 0 ) , ( 0 when ) ( ) , ( 1 d ) , ( h h W h W h W r r r r r r r [9]
where r is the vector position, h a scalar called the smoothing length, a constant parameter taken as 2, the Dirac function and the domain. The integral interpolant is an estimate of the function . Moreover, it can be shown that the integral interpolant of the gradient of a function, if the boundary integrals are neglected, may be written as:

' d ) , ' ( ) ' ( r r r r h W [10]
From a choice of kernel functions available in the software, the frequently used cubic spline function has been chosen for this study. As quoted by [START_REF] Shao | Incompressible SPH method for simulating Newtonian and non-Newtonian flows with a free surface[END_REF], this kernel has compact support, a continuous second derivative and is not very sensitive to particle disorder.

Summation interpolants

For numerical simulation of continuum dynamics, the material is modeled by a = 1, 2, … N particles of mass m a and density a . A particle a interacts with the neighboring particles b belonging to the sphere of influence of a. A sphere of influence (Figure 1), is defined by a kernel function that weighs interactions and a smoothing length that limits interactions in space. Therefore, the integral interpolant of a function (defined in section 3.1.) and its gradient can be approximated by a summation interpolant over a collection of particles:

) , ( ) ( h W m b a b b b b a r r r [11] ) , ( ) ( h W m b a b a b b b a r r r [12]
Equation 12 shows that no mesh is required to numerically evaluate the gradient of the unknown functions but only the gradient of the kernel function. The role of grid or mesh size in FDM or FEM is played by the particle size in SPH: a smaller particle size implies a better spatial resolution of the flow, but requires more CPU. The smoothing length is defined to be proportional to the particle size.

The fact that SPH is a meshless method explains the attractiveness of the method to the dynamics of materials undergoing large relative displacements, while maintaining a Lagrangian frame of reference.

Figure 1. Sphere of influence of a particle a

SPH formulation for fluid flow through porous media

Within the SPH method, particles represent a fixed amount of material. If the fluid enters a porous region, the mass has to be distributed over a larger volume in space. Since the fluid particles should be able to model the progress of the flow front, the proper choice for the SPH formulation in porous media is to let the particles move with the microscopic fluid velocity. Thus, Equations 1 and 3 are amenable to the regular SPH discrete equations, which involves the apparent density and the microscopic velocity:

ab a b b a b a W m t u u d d [13] a a a a ab a b ab a a b b b a a W p p m t u K g u 1 2 2 2 d d [14]
The artificial viscosity contribution ab from [START_REF] Monaghan | Shock Simulation by the Particle Method SPH[END_REF] has been included in the momentum equation in order to stabilize the solution scheme. This type of artificial viscosity is superior to the standard finite element viscosity for SPH. For flow at low speed, the strength of the artificial viscosity may be taken quite small. For fluid flow through porous media, the most relevant physical viscosity contribution is directly entered as the viscous drag term of the momentum equation. Conservation of energy may also be solved with a similar SPH formulation, but since the pressure does not depend on the internal energy and since it will be assumed here that thermal effects on the viscosity may be neglected, the energy equation will not be considered.

Since Equations 13 and 14 allow some compressibility, they will be solved in combination with an artificial equation of state (Murnaghan model) for nearly incompressible liquid:

1 0 0 B p p [15]
with p 0 a reference pressure, the reference density, B the (artificial) bulk modulus and the exponent. The parameters B and are adjusted to keep the sound speed for Equation 15 greater, by at least a factor of 10, than the maximum fluid velocity. Therefore, tuning properly those parameters will keep the compressibility lower than a few percent.

SPH formulation for non-Newtonian fluid flow

A similar SPH formulation of Equations 4 and 5 can be given for non-Newtonian fluid flow: [18] where is the total stress tensor.

W m t 2 2 d d [17] ij a ij a ij a p , ,
To evaluate the shear stresses, the following SPH formulation for the velocity gradients will be adopted [START_REF] Libersky | Smooth Particle Hydrodynamics with strength of materials[END_REF]:

j a ab b i a i b b b j i a x W u u m x u , , , , [19]
This method has been validated for materials with strength in PAM-CRASH by [START_REF] Groenenboom | Numerical simulation of 2D and 3D hypervelocity impact using the SPH option in PAM-SHOCK[END_REF][START_REF] Groenenboom | Copper cylinder impact at high velocity: Numerical Simulation using SPH and Finite Elements in PAM-SHOCK[END_REF]. The current model is similar to that used by [START_REF] Ellero | Viscoelastic flows studied by Smoothed Particle Dynamic[END_REF], with the visco-elastic model replaced by the Carreau model, but it differs from the implementation for incompressible flow by [START_REF] Shao | Incompressible SPH method for simulating Newtonian and non-Newtonian flows with a free surface[END_REF], who arrive at an expression for the viscosity similar to the usually adopted artificial viscosity.

To overcome numerical instabilities during the transient period when the flow is started, the artificial viscosity is linearly reduced and the Carreau model is linearly activated during a user-specified time interval.

Finally, the Murnaghan model (Equation 15) is also used as a constitutive equation. The parameters of this model are adjusted so that the compressibility of the fluid remains lower than few percents.

Mixed FE-SPH

Some advanced LCM modeling and simulations involve hydro-mechanical coupling when a foam core is compressed while manufacturing sandwich materials or when processes such as injection-compression or VARTM are used, for instance. In those cases, either flow induced deformations or forced deformations are present. Such interactions between resin and the moving structures (foam core, mold…) are modeled by one of the sliding interface algorithms available in PAM-CRASH. Sliding interfaces model the interaction between structures and parts of structures that are not permanently connected by standard finite element connectivity conditions. The employed sliding interfaces is based on the well known penalty formulation, where geometrical interpenetrations between so-called slave nodes and contacting faces are penalized by counteracting forces that are in essence proportional to the penetration depth. The contact algorithm will automatically detect when a particle (slave) penetrates any segments (master) of the outer surface of the finite element of the foam. The contact thickness indicates the distance away from a contact face where physical contact is established. In this case the contact thickness should be representative of the particle size. The use of such interaction between the SPH and finite elements has been validated for a range of applications such as sloshing [START_REF] Meywerk | Fluid-structure interaction in crash simulation[END_REF][START_REF] Cha | Industrial applications of PAM-SHOCK using SPH[END_REF], heart valve opening [START_REF] Haack | On the use of a Particle Method for Analyis of Fluid-structure Interaction[END_REF] and impact of aeronautical structures on water [START_REF] Pentecote | Simulation of water impact problems using the Smoothed Particle Hydrodynamics Method[END_REF].

Numerical schemes

Explicit integration

Temporal integration of Equations 13 and 14, is performed by a similar explicit method as for PAM-CRASH finite element nodes:

K g u u 1 2 2 2 2 1 2 1 [21] 2 1 2 1 1 n n a n a n a t u r r [22]
where the superscript denotes the time step, the tilde on u a (Equation 21) the time step to be discussed below and r a the particle coordinate vector. Note that all velocities are defined at intermediate (or half-integer) time steps whereas densities, pressure and positions are defined at integer time steps, rendering second-order accuracy. A numerically stable time step, accounting for finite elements, particles and any sliding interface, is determined automatically by the software during each solution step. Within the porous regions, the flow is usually dominated by the drag force contribution of Equation 21, which in equilibrium should be balanced by the pressure gradient, combined with gravity, if relevant. Hence, the drag term should be integrated in an implicit manner to guarantee numerical stability at an acceptable computational time step. Therefore Ê a should be taken as (u a ) n+1/2 . Since the expression for the drag force only depends on the properties of particle a, the unknowns at the new time step remain uncoupled from the neighboring particles.

Similarly, in the case of the Non-Newtonian formulation, Equations 16 and 17 will be integrated as:

b n ab a n b n a b n n a n a W m t 2 1 2 1 2 1 1 u u [23] g I u u n ab a b n ab b n b a n a b n n a n a W m t 2 2 2 1 2 1 [24] 2 1 2 1 1 n n a n a n a t u r r [25]
where the total stress is the sum of the hydrostatic pressure and the viscous stress (Equation 18).

Transition to the porous regions

When a particle travels from a non-porous to a porous region, it will suddenly experience a high deceleration due to the drag force. To avoid the numerical instabilities that may arise from this situation, interpolation is used. The range corresponding to the representative smoothing length has been found to provide an adequate solution.

Additional SPH features

A well-known but undesirable feature of SPH is the tendency of the particle distribution of moving particles to become less homogeneous. A particle motion correction, introduced as XSPH by [START_REF] Monaghan | Simulating free surface flows with SPH[END_REF], but adapted to the apparent density , is defined by:

b ab b a a b b a a W m t 2 1 d d u u u r [26]
The relative strength of this correction is given by which is usually taken close to 0.3. Other relevant features of the PAM-CRASH/SPH method are planes of symmetry, a variable smoothing length and non-interacting regions in space. The option of mirror particles [START_REF] Morris | Modelling low Reynolds number incompressible flows using SPH[END_REF] allows to introduce mathematically exact planes of symmetry. Non-interacting regions are also defined by parts of space in which particles do not experience mutual interaction. This option allows to define inflow regions of particles moving under the influence of inertia and external forces only. Those options have been used for the applications discussed below.

Examples of applications

Isothermal isotropic filling

Following is an example of mold filling simulation with isotropic permeability. A Newtonian fluid is injected in a complex 2D mold containing 50% of fiber reinforcement. The permeability of the fiber reinforcement is isotropic and constant in time (K=10 -10 m 2 ). The fluid viscosity is 0.001 Pa.s and the fluid inflow velocity is 0.05 m.s -1 at both entries. Comparisons were made on flow front positions and pressure distributions with the well-established FE-based commercial simulation tool PAM-RTM (ESI, France) to model the RTM process including filling and curing. For this case, only the filling part based on Darcy's law, including flow front propagation has been used. Figure 2 shows the results of the simulations after 38 s (a and b) and 54 s (c and d) of injection. The comparison of the flow front positions using the SPH method (a and c) to the PAM-RTM ones (b and d) gives a good match. Minor differences in the filling sequence can be explained by the numerically enhanced compressibility used in the SPH simulations. Increasing the compressibility in SPH provides better agreement but at higher computational cost. In terms of injection pressure responses, the injection pressure obtained at the end of the injection is 0.07 MPa and 0.08 MPa from SPH and PAM-RTM simulations respectively. It is concluded that there is a quite reasonable agreement between the SPH method and PAM-RTM for isotropic porous media flow.

Figure 2. Comparison of an isothermal isotropic RTM filling using the SPH method and PAM-RTM

Isothermal anisotropic filling

Orthotropic effects of the porous medium are easily incorporated in the method by considering the permeability K, to be a tensor. A Newtonian fluid is injected in a rectangular 2D mold containing 50% of fiber reinforcement. The permeability of the upper half of the mold is isotropic (K=1.25 10 -9 m 2 ) while the bottom half is anisotropic (K xx =10 -9 m 2 , K xy =K yx =0.5 10 -9 m 2 and K yy =4.0 10 -9 m 2 ). Note that the off-diagonal terms (K xy and K yx ) in the permeability tensor represent the angle between the simulation domain axes and the main flow directions of the fiber reinforcements. The simulation results shown in Fig. 3 demonstrate that the SPH method (a) is as accurate as PAM-RTM (b), and that the circular (upper half domain) and ellipsoidal (lower half domain) flow front profiles are respectively representative of an isotropic and anisotropic filling.

Hydro-mechanical coupling

RTM can be a very attractive single-step process for sandwich structures manufacturing. The two fiber-reinforced skins and the impermeable core are laid up in a mold cavity. Hydro-mechanical coupling is present when the pressure field developing while saturating the fiber-reinforced skins leads to shifting or compression of the core [START_REF] Wirth | Experimental analysis of core crushing and core movement in RTM and SRIM foam cored composite parts[END_REF][START_REF] Binetruy | Foam core deformation during liquid molding of sandwich: Modelling and experimental analysis[END_REF]. For such process, the skins are modeled as porous regions in which Lagrangian SPH particles are evolving. Figure 4 shows particles and finite elements modeling resin and sandwich core respectively.

Figure 4. Snapshot of the local foam deformation when hydro-mechanical coupling is present

The injection is performed at constant flow rate. The equations for the deformation and stresses in the foam are solved by the explicit finite element algorithm. The choice has been made to employ the first-order, iso-parametric and under-integrated hexahedron elements with hourglass control of PAM-CRASH. Although specific material models for various types of foam are available within this software, a linear-elastic material behavior with initial internal compressive stress is assumed since the compression remains limited. When the foam is deformed because of resin pressure built-up, more volume becomes available to the resin. Both porosity and permeability increase leading to a discontinuity in the pressure response. Userdefined subroutines have been developed to include porosity and permeability variations as a function of the local foam core displacement. This behavior simulated with the presented method agrees well with the physics and the experimental data. Full details, pressure and foam deformation results on hydro-mechanical coupling simulation using the mixed FE-SPH method are given in [START_REF] Comas-Cardona | A generic mixed FE-SPH method to address hydro-mechanical coupling in liquid composite moulding processes[END_REF].

Saturated permeability of fabrics

In order to perform accurate LCM filling simulations, physical parameters such as fabric permeability are needed. It is well known that the experimental measurements of that fabric property is very delicate. Usually a fluid is injected at constant pressure (or constant flow rate) through several layers of fabric of cross section A and length L. Then pressure loss p and flow rate Q are measured, and the permeability K is calculated using Darcy's law:

p A QL K [27]
The simulation of flow through a periodic cell should provide a reliable solution and avoid that experimental procedure. Numerically, an injection at constant flow rate could be performed. Then the pressure loss can be measured and the permeability is also given by Equation 27. However, another solution would be to start from a filled cavity, apply periodic boundary conditions and drive the fluid through the unit cell by gravity. For such flows, Darcy's law can be written as:

g u K x xx ˆ [28]
where û x is the steady state average fluid velocity in the flow direction. A validation of this technique has been done for a Newtonian fluid flow through a rectangular pipe. Results of the validation for two viscosities are given in Table 1 using a 5 x 10 x 20 particle lattice. It can be seen that the theoretical permeability of a rectangular pipe (K th ) matches well the permeability obtained from the simulations (K SPH ). Calculations on a full size fabric cell (Figure 5) is currently under investigation.

Discussion and conclusions

This paper presents the mathematics and numerical schemes for a unified approach to LCM processes at various scales of discretization, using a mixed FE-SPH approach. It has been demonstrated that the SPH method allows to simulate mold filling in both isotropic and anisotropic porous media, and viscous flow to determine permeability. In combination with Finite Elements, hydro-mechanical coupling during an RTM process may be simulated as well.

Flow front and pressure levels obtained with the FE-SPH method may slightly differ from experiments because of the limited compressibility allowed in the solution (Equation 15). Therefore a compromise need to be found between reducing compressibility and obtaining results at a reasonable CPU. The level of these differences is, however, surpassed by the variations in the results due to the inaccuracies in the measured physical input parameters. For regular filling, these inaccuracies propagate in a linear manner, but with hydro-mechanical coupling, the variations in the results may become much bigger.

The generic FE-SPH method is not expected to prevail for regular resin filling simulations since FDM or FEM filling tools are faster to use. However, when coupling is of interest, the SPH-FE method does provide an attractive alternative. A major advantage of this method with respect to PAM-RTM and similar software is that it allows the interactive coupling with the structural deformation leading to a completely different flow because of permeability variations and the formation of channels (preferential flow) in between the fiber reinforcement and the foam for instance. Therefore, the potential of the method relies on its Lagrangian formulation as well as its ability to interact with finite elements. Also, it would be of interest for mesoscopic and microscopic flows in unit fabric cells. More work could involve the hydro-mechanical coupling between the flow and the fibers, and the permeability measurement for sheared fabrics for instance.

Future developments for a more powerful LCM simulation tool would be to incorporate the equation of energy to solve for non-isothermal processes. This original method provides a suitable solution for various coupled phenomena since the integration of many physical effects is rather straightforward in the SPH formulation and implementation. This feature is an important asset of the method

Figure 3 .

 3 Figure 3. Comparison of an isothermal anisotropic RTM filling using the SPH method and PAM-RTM

Figure 5 .

 5 Figure 5. Snapshot of the filling of a periodic fabric unit cell

Table 1 .

 1 Flow through a rectangular pipe (Re is the Reynolds number)

Acknowledgements

The authors would like to thank S.V. Lomov and I. Verpoest (KUL, Belgium) for providing the fabric unit cell.


S Comas-Cardona

P Groenenboom

C Binetruy
email: binetruy@ensm-douai.fr

P Krawczak

A generic mixed FE-SPH method to address hydro-mechanical coupling in liquid composite moulding processes

Keywords: A. Layered structures, C. Computational modelling, C. Numerical analysis, E. Resin transfer moulding (RTM)

Resin Transfer Moulding (RTM) can be a very attractive single-step process for sandwich structure manufacturing. During an RTM injection, large deformation of the foam core material is present due to the pressure field developing while saturating the fibre-reinforced skins. Such effects induce poor dimensional accuracy of skins and core thickness. Moreover, the foam core crushing increases filling time and injected resin amount. In order to model this hydro-mechanical coupling, the use of a mixed Finite Element and Smoothed Particle Hydrodynamics (FE-SPH) method is presented in this study. This method combines Lagrangian particles and finite elements, modelling, respectively, the resin and the sandwich core. Computed pressure responses during the filling stage are compared to symmetric sandwich manufacturing test results. The potential of this original numerical method is also discussed.

Introduction

Resin Transfer Moulding (RTM) can be a very attractive single-step process for sandwich structures manufacturing. The two fibre-reinforced skins and the core are laid up in a mould cavity. Both skins are saturated and bonded to the core either using a pressure bucket to perform an injection at constant pressure or using a piston to perform an injection at constant flow rate. There is no need for an extra adhesive bonding step compared to other sandwich manufacturing processes. During sandwich manufacturing using RTM, the pressure field developing while saturating the fibrereinforced skins can lead to large shifting [START_REF] Al-Hamdan | Dynamic core movements during liquid moulding of sandwich structures[END_REF][START_REF] Wirth | Experimental analysis of core crushing and core movement in RTM and SRIM foam cored composite parts[END_REF] or compression [START_REF] Wirth | Experimental analysis of core crushing and core movement in RTM and SRIM foam cored composite parts[END_REF][START_REF] Binetruy | Foam core deformation during liquid molding of sandwich: modelling and experimental analysis[END_REF] of the core as found in experiments. When such injections are performed at constant flow rate, it has been shown that measuring pressure at the resin inlet can be an appropriate sensor for the detection of hydro-mechanical coupling during sandwich manufacturing [START_REF] Binetruy | Foam core deformation during liquid molding of sandwich: modelling and experimental analysis[END_REF]. Moreover, during resin injection stage, fibre reinforcements can also be submitted to compression [START_REF] Al-Hamdan | Dynamic core movements during liquid moulding of sandwich structures[END_REF][START_REF] Han | Fibre mat deformation in liquid composite molding. I: experimental analysis[END_REF]. Such effects are damaging to the part because they induce poor dimensional accuracy of both skin and core thickness. Also, filling time and injected resin amount can be increased due to the core crushing.

Previous works focused on addressing those manufacturing issues by modelling the mechanical behaviour of fibre reinforcement under compression [START_REF] Saunders | Compression in the processing of polymer composites 1. A mechanical and microstructural study for different glass fabrics and resins[END_REF][START_REF] Batch | Compaction of fibre reinforcements[END_REF] and the hydromechanical coupling involved in composite manufacturing processes using finite element methods (FEM) [START_REF] Wirth | Experimental analysis of core crushing and core movement in RTM and SRIM foam cored composite parts[END_REF][START_REF] Pham | Simulation of compression resin transfer molding with displacement control[END_REF], mixed control volume/finite difference methods (FDM) [START_REF] Han | Fibre mat deformation in liquid composite molding. II: modelling[END_REF] or a onedimensional Lagrangian formulation [START_REF] Farina | Flow in deformable porous media: modelling and simulations of compression moulding processes[END_REF]. Numerical simulation of fluid flow is most frequently performed using FDM or FEM within an Eulerian reference frame. When the geometry of the flow domain changes, due to a free surface or interfaces with moving structures, or for mixing of multiple fluids, special solutions are required within these numerical schemes. There also exist many numerical simulation tools for flow through porous media, but it is often not possible to combine them with flow in non-porous regions.

In this paper, it will be shown how the Smoothed Particle Hydrodynamics (SPH) method provides an attractive alternative solution to simulate flow-induced foam core compression during symmetric sandwich manufacturing based on Binetruy and Advani's former work [START_REF] Binetruy | Foam core deformation during liquid molding of sandwich: modelling and experimental analysis[END_REF]. Foam core shifting will not be taken into consideration. The SPH method has previously been used to simulate RTM filling [START_REF] Sawley | Modelling of flow in porous media and resin transfer moulding using smoothed particle hydrodynamics[END_REF] and gave fairly good qualitative results in terms of resin flow front profiles. In order to model this coupling, a mixed Finite Element and Smoothed Particle Hydrodynamics (FE-SPH) method has been chosen. An SPH package is available within the commercial simulation tool PAM-CRASH (ESI, France) and has been further developed to solve for viscous flow in porous media. The following section presents the fundamentals involved with flow within porous media and a brief description of the mixed FE-SPH method. A comparison with another simulation tool for isotropic flow during RTM filling will be presented. Then, pressure responses computed from the filling simulations will be compared to symmetric sandwich manufacturing test results where hydromechanical coupling is present. Finally, the potential of this alternative and generic numerical method will be discussed.

Fundamentals

Fluid flow dynamics within a porous medium

For transport of a liquid within porous media, the flow is dominated by the viscous drag the liquid experiences. Therefore, these effects need to be included whenever the liquid is inside a porous region. It is possible to derive the conservation of mass and momentum equations for flow in porous media using the volume averaging method on the conservation equations for a two-phase medium in which one of the materials is fixed in space [START_REF] Tucker | Governing equations for flow and heat transfer in stationary fibre beds[END_REF][START_REF] Pillai | Governing equations for unsaturated flow through woven fibre mats. Part 1. Isothermal flows[END_REF][START_REF] Groenenboom | Application of the mixed SPH-FE approach in PAM-CRASH to deformable porous media[END_REF]. Considering that the volume fraction available for the fluid remains fixed and is called porosity, constant fluid viscosity, Newtonian fluid and laminar flow, the general forms of the equations to be solved (conservation of mass and momentum) are, respectively:

vh vt C V$ðhũÞ Z 0 ( 1
)
vhũ vt C V$hũũ Z KfVp C fmVðV$ũÞ C hg K f f d ( 2
)
where f is the fluid volume fraction or porosity, ũ the fluid velocity, m the viscosity, p the pressure, g the gravitational acceleration, f d the drag force and h the apparent density defined as:

h Z m f V t Z fr (3)
where m f is the fluid mass within a control volume V t , and r the fluid density. Even though the flow of resin during composite manufacturing may be considered as incompressible, the governing equations (Eqs. ( 1) and ( 2)) are written for a compressible fluid. The reason for this will be given in the following section.

According to Tucker et al. [START_REF] Tucker | Governing equations for flow and heat transfer in stationary fibre beds[END_REF], for composite moulding processes the viscous term in Eq. ( 2) is negligible. These authors also mentioned that the drag force might be considered as linearly dependent on the fluid velocity when inertia is negligible

f d Z fmK K1 $ũ ( 4
)
where K is the permeability. Taking into account the previous assumptions and combining Eqs. ( 1) and (2) leads to:

h vũ vt C ũðV$ũÞ Z KfVp C hg K f 2 mK K1 $ũ (5)
In this study, in order to simulate the flow for all regions of the computational domain, Eq. ( 5) provides the basis of the solution with the viscous drag term only being used within the porous domain. As a limiting case, within the porous region (i.e. where the flow is dominated by the viscous drag and the Reynolds number based on the pore size is much lower than 1) the LHS of Eq. ( 5) may be neglected and hence Darcy's equation may be derived

ṽ Z K K m ðVp K rgÞ (6)
where ṽ Z fũ is the transport velocity.

SPH formulation

The SPH is a Lagrangian and mesh less method which was originally developed to solve astrophysical problems [START_REF] Lucy | A numerical approach to the testing of fusion process[END_REF][START_REF] Gingold | Smoothed particle hydrodynamics: theory and application to non-spherical stars[END_REF] and has later been applied also to the flow of compressible fluids [START_REF] Monaghan | Simulating free surface flows with SPH[END_REF][START_REF] Schlatter | A pedagogical tool using smoothed particle hydrodynamics to model fluid flow past a system of cylinders[END_REF][START_REF] Hyncik | Lagrangean approach to flow simulations by smoothed particle hydrodynamics[END_REF][START_REF] Shao | Incompressible SPH method for simulating Newtonian and non-Newtonian flows with a free surface[END_REF][START_REF] Pentecote | Simulation of water impact problems using the smoothed particle hydrodymanics method[END_REF]. In SPH, the fluid is modelled by N particles of mass m i . Each particle i will interact with all neighbouring particles j that belong to the sphere of influence of the particle i (Fig. 1). The sphere of influence is defined by: -a kernel function W in order to weigh interactions according to the distance between the particles; -a smoothing length h to limit interactions in space. Several kernel functions are available in the literature and in the used software. For this study, the cubic spline kernel has been chosen [START_REF]PAM-CRASH Notes Manual[END_REF][START_REF] Fulk | A numerical analysis of smoothed particle hydrodynamics[END_REF]. Note that the kernel functions have properties similar to the Dirac function. Another specificity of the SPH method is the particle size. It plays a similar role as the grid or mesh size in FDM or FEM: a smaller particle size implies a better spatial resolution of the flow, but more CPU. The smoothing length is defined to be proportional to the particle size.

Within the SPH method, particles represent a fixed amount of fluid. If the fluid enters a porous region, the mass has to be distributed over a larger volume in space. Since the fluid particles should be able to model the progress of the flow front, the proper choice for the SPH formulation in porous media is to let the particles move with the microscopic fluid velocity. Therefore, the advancement of the front is a direct consequence of the motion of the particles representing the resin, both inside and outside the porous region. This is one advantage of using the Lagrangian SPH description.

It is straightforward to derive the SPH formulation corresponding to Eqs. ( 1) and ( 5) [START_REF] Lucy | A numerical approach to the testing of fusion process[END_REF][START_REF] Gingold | Smoothed particle hydrodynamics: theory and application to non-spherical stars[END_REF][START_REF] Monaghan | Simulating free surface flows with SPH[END_REF][START_REF] Schlatter | A pedagogical tool using smoothed particle hydrodynamics to model fluid flow past a system of cylinders[END_REF][START_REF] Hyncik | Lagrangean approach to flow simulations by smoothed particle hydrodynamics[END_REF][START_REF] Shao | Incompressible SPH method for simulating Newtonian and non-Newtonian flows with a free surface[END_REF][START_REF] Monaghan | Shock simulation by the particle method SPH[END_REF] which in this case involves the apparent density and the microscopic velocity:

dh i dt Z X j m j ðũ i K ũj Þ$V i W ij (7) dũ i dt ZKf i X j m j p j h 2 j C p i h 2 i CP ij ! V i W ij Cg Kf 2 i mK K1 i ũi h i (8)
where W ij is the SPH smoothing kernel, m j the mass of the particle j. In order to stabilize the solution scheme, the standard Monaghan-Gingold artificial viscosity contribution P ij has to be added to the momentum equation [START_REF] Monaghan | Shock simulation by the particle method SPH[END_REF]. It has been shown that when numerical shock phenomena are no longer dominating the solution, this artificial viscosity may be used to mimic the physical viscosity, but it may be difficult to identify the proper parameters [START_REF] Schlatter | A pedagogical tool using smoothed particle hydrodynamics to model fluid flow past a system of cylinders[END_REF]. The possibility to include non-Newtonian viscosity in the fluid flow is currently being investigated.

For the current study, only the most relevant contribution of viscosity to Eq. ( 8), the viscous drag term, is included. Numerical solution of Eqs. ( 7) and ( 8) allows to model the flow in both porous and non-porous regions of a domain. When the particle is in the porous region, the RHS of Eq. ( 8) becomes very large and the equation becomes equivalent to Darcy's law. However, when the fluid particle is in a non-porous region, the viscous drag term is dropped. Orthotropic effects of the porous medium are easily incorporated by considering the permeability K to be a tensor in Eq. ( 8). Usually, the flow of resin during moulding is considered as that of an incompressible fluid. In that case, Eq. ( 1) reduces to V$ũ Z 0 and the pressure should be obtained from a Poisson type of equation. Such a system of equations can be solved efficiently by an implicit method which allows for a relatively large time step. In the present study, however, the solution time step needs to be small. The first reason is that the motion of the flow front requires an accurate solution in time and space when the resin passes the border of the porous region. A second reason is that the weak coupling of the flow with the foam, as discussed below, requires an explicit solution with a sufficiently small time step. Hence, the original Eqs. ( 7) and (8) (corresponding to Eqs. ( 1) and ( 2)) will be solved in combination with an artificial equation of state (Murnaghan model) for nearly incompressible liquid [START_REF] Sawley | Modelling of flow in porous media and resin transfer moulding using smoothed particle hydrodynamics[END_REF][START_REF] Monaghan | Simulating free surface flows with SPH[END_REF][START_REF] Schlatter | A pedagogical tool using smoothed particle hydrodynamics to model fluid flow past a system of cylinders[END_REF][START_REF] Hyncik | Lagrangean approach to flow simulations by smoothed particle hydrodynamics[END_REF][START_REF] Shao | Incompressible SPH method for simulating Newtonian and non-Newtonian flows with a free surface[END_REF][START_REF] Pentecote | Simulation of water impact problems using the smoothed particle hydrodymanics method[END_REF]:

p Z p 0 C B r r 0 g K 1 ( 9
)
with p 0 a reference pressure, r 0 the reference density, B the (artificial) bulk modulus and g the exponent. The parameters B and g of Eq. ( 9) are adjusted to keep density variations lower than a few percent. SPH solutions for incompressible flow have also been developed, but such solutions have not been pursued for this study.

Interaction of SPH particles with finite elements

The interaction of the particles with the finite elements, representing, respectively, the resin and the foam in this study, is modelled by the sliding interface algorithms available within PAM-CRASH [START_REF]PAM-CRASH Notes Manual[END_REF]. Sliding interfaces model the interaction between structures and parts of structures that are not permanently connected by standard finite element connectivity conditions. Most of these sliding interfaces are based on the well-known penalty formulation, where geometrical interpenetrations between contacting faces are penalized by counteracting forces that are in essence proportional to the penetration depth. The contact algorithm will automatically detect when a particle (slave) penetrates any segments (master) of the outer surface of the finite element of the foam. The contact thickness indicates the distance away from a contact face where physical contact is established. In this case the contact thickness should be representative of the particle size. The use of such interaction between the SPH and finite elements has been validated for a range of applications such as sloshing [START_REF] Meywerk | Fluid-structure interaction in crash simulation[END_REF][START_REF] Cha | Industrial applications of PAM-SHOCK using SPH[END_REF], heart valve opening [START_REF] Haack | On the use of a particle method for analyis of fluidstructure interaction[END_REF] and impact of aeronautical structures on water [START_REF] Pentecote | Simulation of water impact problems using the smoothed particle hydrodymanics method[END_REF].

Numerical schemes

Due to the motion of a free surface and of the interface with the deformable structural parts (e.g. the foam core during sandwich manufacturing) the computational time step will be required to be small and an explicit solution scheme in combination with the nearly-incompressible fluid model will provide a suitable solution. Pressure, gravity and viscous drag (through the viscosity to permeability ratio) contribute to the momentum equation (Eq. ( 8)). For low permeability, the drag force may become so high that the direct explicit solution becomes unstable. Therefore, the drag force contribution to the SPH momentum is integrated in an implicit manner. Furthermore, a special solution has been implemented for particles very close to the interface between the porous and non-porous regions in order to avoid numerical oscillations.

Comparison between the predictions of FE-SPH method and PAM-RTM

Before using the mixed FE-SPH method to solve for hydro-mechanical coupling during sandwich manufacturing, comparisons were made on flow front positions and pressure distributions with the commercial simulation tool PAM-RTM (ESI, France). The mould filling part of this software provides the solution for transient flow based on Darcy's law using linear non-conforming finite elements [START_REF] Gauvin | Permeability measurement and flow simulation through fiber reinforcement[END_REF] and incorporates tracking of the flow front during injection. Following is an example of mould filling simulation with isotropic permeability. A Newtonian fluid is injected in a complex 2D mould containing 50% of fibre reinforcement (Fig. 2). The permeability K of the fibre reinforcement is isotropic and constant in time (KZ10 K10 m 2 ). The fluid viscosity is 0.001 Pa s and the fluid inflow velocity is 0.05 m s K1 . Fig. 2 shows the results of the simulations after 38 s (a and b) and 54 s (c and d) of injection. The comparison of the flow front positions using the SPH method (a and c) to the PAM-RTM ones (b and d) gives a good match. Minor differences in the filling sequence can be explained by the numerically enhanced compressibility used in the SPH simulations. Increasing the compressibility in SPH provides better agreement but requires higher computational effort. In terms of injection pressure responses, the maximal injection pressure obtained at the end of the injection is 0.07 and 0.08 MPa from SPH and PAM-RTM simulations, respectively. It is concluded that there is a quite reasonable agreement between the SPH method and PAM-RTM for porous media flow.

Experimental

Material characterization

An exhaustive materials characterization is performed at room temperature in order to feed the code with material input data. A non-reactive Newtonian fluid (Di-Octyl Phtalate, DOP) is used. Its viscosity, measured using a Brookfield viscometer, is 0.07 Pa s. The permeability of the woven fabric (1500S3, Chomarat Composites) is measured for fibre volume fractions varying from 0.43 to 0.57. The variation can be fitted to a second order polynomial as:

K Z 1:08 !10 K8 V 2 f K 1:42 !10 K8 V f C 4:6 !10 K9 (
10) where V f is the fibre volume fraction and the permeability K is in m 2 . The areal weight of the fabrics is 1500 g m K2 . The foam core (NA1100, Alve ´o) has also been characterized in compression. For the pressure levels used in this study, the foam can be considered as a linear elastic material and its bulk modulus is 1.2 MPa.

Sandwich manufacturing

Sandwich panels are manufactured at room temperature using a mould which cavity dimensions are 500!200! 16 mm 3 (Fig. 3). The mould cavity is filled up with a 10 mm-thick soft closed-cell polyethylene foam core and two skins each consisting of 3 plies of glass woven fabric. The fibre volume fraction reached in both skins is close to 53%. Once the mould is closed, a fluid (DOP) is injected at constant flow rate (15 cl min K1 ). Pressure at the injection gate (P1) and 29 cm away from the gate (P2) are monitored and recorded throughout the injection (Fig. 4). The pressure response at the gate location (P1 curve) can be separated into two regions: a rapid initial pressure increase during the first 25 s followed by a slower pressure increase. In the first region, the linearity between pressure and time verifies Darcy's law for unidirectional flow while the porous medium (fibre reinforcement) is not deformed. After 25 s of injection, foam core compression occurs due to a local pressure built-up in the mould cavity. In the neighbourhood of the gate, both porosity and permeability of the porous media increase, creating a discontinuity of slope of the pressure response-the resistance of the fluid to flow through the porous media decreases [START_REF] Binetruy | Foam core deformation during liquid molding of sandwich: modelling and experimental analysis[END_REF].

Simulation

In order to simulate the hydro-mechanical coupling observed experimentally, the mixed FE-SPH method is used. The foam core is modelled with FE as an elastic material, which is initially loaded under compressive stress. This initial compressive load is calculated from the equilibrium between the fibre reinforcement and the foam compression at the closing of the mould [START_REF] Binetruy | Foam core deformation during liquid molding of sandwich: modelling and experimental analysis[END_REF]. The fibrereinforced skins are modelled as porous regions in which Lagrangian particles are evolving.

Special routines have been developed and implemented within PAM-CRASH to model as closely as possible the hydro-mechanical coupling. When the foam is being compressed, more volume becomes available to the flow of resin leading to expansion of the porous material (Fig. 5). When the foam deformations exceed a specific minimum value, a thin layer of free volume between the foam core and the porous material is created. Practical limitations do not allow to include this effect directly in the geometry of the computational model, but simple mathematical expressions are derived for the porosity and permeability as a function of the local core displacement. Therefore the user-defined routines (for porosity and permeability) are implemented and both parameters are calculated at each time step based on the local foam core compression (Fig. 6). As mentioned earlier, the fluid pressure response is transferred to the foam core using contact algorithms. Thus, no extra boundary conditions in between the porous and solid regions are required.

Fig. 7 shows the simulation results in terms of liquid pressure and foam core deformation at the injection gate. The pressure response increases with respect to time until the foam starts to deform (tZ28 s). After reaching a foam deformation of few tenths of millimetres, the pressure levels off. This behaviour simulated with the mixed FE-SPH method corresponds to the physics described experimentally. A comparison of the liquid pressure from the simulation with the experimental one also shows a quantitative good agreement. Fig. 5 is a close-up look at the foam deformation in the neighbourhood of the injection gate after 35 s of injection. Because of technical issues it was not possible to measure it experimentally, using a contact-less displacement sensor for instance, since the experiments were carried out in an industrial mould. However, the foam core deformation (few tenths of millimetres) and the pressure response effects obtained in this study from the simulation, are in the same range as those obtained by Binetruy et al. [START_REF] Binetruy | Foam core deformation during liquid molding of sandwich: modelling and experimental analysis[END_REF] and Wirth et al. [START_REF] Wirth | Experimental analysis of core crushing and core movement in RTM and SRIM foam cored composite parts[END_REF] for RTM using very similar experimental conditions. Therefore, the accuracy of the method allows us to properly predict pressure and core crushing responses during the injection. Finally, the results of the study show that the simulation tool should be very sensitive. As a matter of fact, few tenths of millimetres of foam core crushing have a tremendous effect on the pressure response at the injection gate.

Discussion and conclusion

The pressure curve responses obtained from experiments and from simulation may differ because of inaccuracies of the physical input parameters (viscosity, foam bulk modulus, permeability.). For porous media flow according to Darcy's law, a typical 10% error on the measurement of the permeability of the fibre network will induce a 10% error on the pressure response for injections at constant flow rate. However, when hydro-mechanical coupling is relevant, the errors will propagate in a non-linear manner and will be quite difficult to estimate. Errors can also be present on the experimental curve since race-tracking and edgeeffects may be present. The most important result of the simulation is that the liquid pressure response increases with respect to time until the foam starts to deform, then because of the permeability increase of the porous media, the liquid pressure reaches a plateau, as found in experiments.

In view of these uncertainties it was not found meaningful to conduct extensive testing of the sensitivity of the results on some numerical parameters (artificial viscosity, penalty coefficients of contacts, and coefficients of the Murnaghan equation). From previous studies on SPH and the coupling to finite elements, it is known that the influence on the results remains limited. The coupled SPH-FE solution is not expected to become a practical alternative to the existing Finite Element of Finite Difference solution for regular filling. However when coupling needs to be accounted, the SPH-FE method does provide an alternative. The main advantage of this method with respect to PAM-RTM and similar software is that it allows the interactive coupling with the structural deformation leading to a completely different flow because of permeability variations and the formation of channels (preferential flow) in between the fibre reinforcement and the foam for instance.

The mixed FE-SPH method has been used to solve for an example of hydro-mechanical coupling. A step toward the simulation of hydro-mechanical coupling has been reached using a very generic and versatile method. Another direct application would be to simulate the displacement of a foam core during sandwich manufacturing with RTM for instance. Also, the mixed FE-SPH method should have the potential to solve for non-isothermal fluid flows. Future developments would focus on the implementation of the equation of energy to simulate thermo-hydro-mechanical coupling in liquid composite moulding processes or include mechanical behaviour of fabrics for foam core shifting modelling for instance.

Fig. 1 .

 1 Fig.1. Sphere of influence of a particle i.

Fig. 2 .Fig. 3 .

 23 Fig. 2. Flow front position comparison for isotropic permeability filling simulations. Fig. 3. Experimental setup.

Fig. 5 .

 5 Fig. 5. Part of the domain used for the FE/SPH simulation. Grey scale of the mesh represents the stress levels.

Fig. 6 .

 6 Fig. 6. Flow chart of the permeability K and porosity f adjustments.

Fig. 4 .

 4 Fig. 4. Experimental pressure responses (P 1 : at the gate. P 2 : 29 cm from the gate).

Acknowledgements

The authors would like to thank DGA (De ´le ´gation Ge ´ne ´rale pour l'Armement, France) for financial support for this study.