Emmanuelle Abisset-Chavanne
email: emmanuelle.abisset-chavanne@ec-nantes.fr

Francisco Chinesta

Toward an optimisation of the reactive resin transfer molding process: thermo-chemico-mechanical coupled simulations

Keywords: RTM process, Multiphysic coupled model, Numerical simulations

In a view of reactive RTM process optimization, a numerical strategy is proposed to perform simulation of the filling stage of such a process, including the numerous thermo-chemico-mechanical couplings involved. Numerical simulations are then performed illustrating the capabilities of the strategy and underlying the numerical features that have to be taken into account. Finally, in order to build virtual charts for processability windows definition, the possibility to combine this numerical strategy with a multiparametric solver is addressed.

Introduction

In recent years, thermoplastic composite materials are more and more considered in industrial sectors, as the aerospace and automotive industries, because they present high mechanical performances, an environmental interest and a potential to reduce the forming process time. In parallel, the liquid composites moulding (LCM) processes and, particularly, the resin transfer moulding (RTM), are being increasingly used in the manufacture of fiber-reinforced composite materials because of its potential to become a high-speed process for the manufacturing of components at a low cost. The objective of this work is to combine these material and process advantages by developing a RTM process modelling adapted to thermoplastic materials.

As for the classical RTM process, it consists in injecting a low viscosity fluid into a mould cavity containing a preshaped dry fiber reinforcement. During this injection phase, the fluid impregnates and wets out the fibers. Once the mould is filled, the system is cooled down. In our case, because polymerized thermoplastic resins have a high viscosity, the prepolymers are mixed just before the injection into the mould and the polymerization reaction occurs during the injection phase. In other words, the resin simultaneously fills the mould and polymerizes-and can continue after-to become a very viscous thermoplastic, and then crystallization can follows in the case of semicristalline polymers. The mould in this case acts also as a sort of chemical reactor. Therefore, this leads to numerous interactions between the two mechanisms that modify the final properties of the part.

The aim of the work proposed here is to develop a numerical tool that allows to predict the final properties of the cured material by taking into account the numerous couplings and phenomena occurring during the reactive RTM process, and particularly during the mould filling stage. Therefore, the proposed model integrates the different physics occurring during the injection: the filling process, the polymerization reaction, the associated temperature increase and the viscosity variation with polymerization degree and temperature during the injection stage.

The modelling of such phenomena have already been addressed separatly: the mould filling process in one side [START_REF] Trochu | Numerical analysis of the resin transfer molding process by the finite element method[END_REF][START_REF] Advani | Flow and rheology in polymer composites manufacturing[END_REF][START_REF] Chinesta | Some difficulties in the flow front treatment in fixed mesh simulations of composites forming processes[END_REF] and the chemical reactions and couplings in the others (e.g. in [START_REF] Luisier | Time-temperaturetransformation diagram for reactive processing of polyamide 12[END_REF]). Then, the idea here is to assemble and couple these different models in a complete one and to introduce it in a numerical tool to simulate all the physics occurring in the mould.

The numerical resolution of mould filling problems has already been investigated in the literature but technical points still have to be circumvented in order to build robust numerical tools. In particular, a key-point is the numerical treatment of the transient free surface associated with the liquid flow front. Indeed, in such injection problems, the properties of the flow (namely here the temperature and the polymerization degree) are not defined in the empty part of the mould. A particular care has then to be taken concerning the treatment of the transport problems and the initial conditions when discretized in an Eulerian framework, for all the different quantities considered.

In the numerical tool here proposed, the resolution of the transport problem of the different quantities is based on the use of a fluid-fraction function within a volumeof-fluid scheme [START_REF] Chinesta | Some difficulties in the flow front treatment in fixed mesh simulations of composites forming processes[END_REF]. The coupling between the different physics is introduced by means of adequate source terms in the different transport equations.

In the following, we begin by introducing all the equations that allow us to model the plysics and couplings occurring during the mould filling stage of a reactive RTM process. Then, the numerical integration scheme of the model equations is developed.

In a third part, some simulations of the reactive RTM process are performed to demonstrate the capability of the model and the necessity to take into account all the couplings involved, and in particular the coupling of the viscosity increase with the polymerization degree in order to ensure a complete mould filling. We will analyse carefully the numerical troubles introduced by the time discretization on the quantities transported, and in particular on the temperature evolution.

Finally, the possibility of using this simulation facility in a view of performing process optimization is analysed within the model reduction framework.

Mechanical and chemical modelling

In order to simulate the reactive RTM process, two main points have to be addressed:

-the resolution of the injection problem: modelling the filling process and the transport of the quantities that characterise the resin state (temperature, polymerization degree ...) -the introduction of the couplings between all the physics that occur in the mould during the filling: polymerization kinetic, temperature evolution and evolution of the viscosity (η) with the polymerization degree (α) of the resin.

The first point has been solved by using the strategy proposed in [START_REF] Sanchez | Towards an efficient numerical treatment of the transport problems in the resin transfer molding simulation[END_REF] and whose main ideas are revisited and generalized in what follows. The second point has been treated following [START_REF] Luisier | Time-temperaturetransformation diagram for reactive processing of polyamide 12[END_REF] that considers PA12, by introducing some source terms in the transport equations governing the evolution of the quantities of interest.

Fluid flow treatment

The flow problem is defined in the mould cavity which can be divided into two subdomains

= f (t) e (t)
f (t) representing the part of the mould saturated with the resin and e (t) the empty part at time t. The common boundary of these two parts is the flow front which is the most delicate entity to treat in such a flow modelling. Indeed, the flow front position has to be updated at each time step, and with it, the boundary conditions of each mould parts e and f .

The main difficulty when simulating this problem by using standard fixed mesh methods is the dependency of the front position with the mesh [START_REF] Scardovelli | Direct numerical simulation of free-surface and interfacial flow[END_REF][START_REF] Sethian | Level set methods for fluid interfaces[END_REF]. Moreover, even if advanced techniques are used to deal with the front updating [START_REF] Hirt | Volume-of-fluid (VOF) method for the dynamics of free boundaries[END_REF][START_REF] Osher | Fronts propagating with curvature dependent speed: algorithms based on Hamilton Jacobi formulations[END_REF][START_REF] Cruchaga | Movinginterface computations with the edge-tracked interface locator technique (ETILT)[END_REF][START_REF] Pilliot | Second-order accurate volume-offluid algorithms for tracking material interfaces[END_REF], a tricky problem remains: the definition of the initial conditions to be applied to the transported quantities in the empty part of the mould (temperature or polymerization degree) as they are physically not defined in absence of fluid.

To overcome these difficulties, in [START_REF] Sanchez | Towards an efficient numerical treatment of the transport problems in the resin transfer molding simulation[END_REF], authors propose the use of a volume-of-fluid method (derived from the one defined in [START_REF] Sweby | High resolution schemes using flux limiters for hyperbolic conservation laws[END_REF]) coupled with the introduction of a fluid fraction function I (x, t) which takes values between 0 and 1 and represents the resin saturation level at each point:

I (x, t) = 1 if x ∈ f (t) 0 if x ∈ e (t) (1)
Thus, using this function, it is not necessary anymore to divide the mould into two parts, we just have to calculate the evolution of the function I in the whole volume to represent the evolution of the flow domain during the filling.

By considering the quantities of interest (unknown fields) multiplied by function I, the calculation of the different quantities can then be done over the whole domain, avoiding the problem of quantities definition in the empty part of the mould. Indeed, let consider for example the calculation of the pressure in the mould. The pressure distribution is governed in f by the Darcy's law combined with the fluid incompressibility, and is equal to zero in e :

∇(-K/μ • ∇p) = 0 ∀x ∈ f (t) p(x) = 0 ∀x ∈ e (t) (2)
Thanks to the use of the I function, we can combine this two behaviours to define an extended weak form of the pressure problem over the whole domain :

∇p * (K/μ) -1 ∇p • I + p * p • (1 -I ) d = 0 (3)
and likewise for all the quantities involved in the filling simulation. By this way, the calculation is made over the whole domain. We do not have to consider anymore the flow front tracking, because it is captured from the knowledge of the fluid fraction I. Thanks to this strategy, we then avoid the first tricky problem of injection simulations: the flow front position update.

Concerning the second problem-the initialization of the transported quantities in the empty part of the mould-it is developed in section "Transport and coupling: updating the other fields".

Thermal-chemical-rheological coupling Numerous works have already been published on the numerical treatment of injection processes. But, the principle of a reactive RTM process is that the polymerization reaction occurs during the filling stage leading to the viscosity (η) increase as the degree of polymerization (α) increases and the resin temperature (T) evolves. Then, we have introduced in the model the transport of the characteristic quantities of the resin (namely the temperature and the polymerization degree) and their couplings.

The chemical kinetics and the couplings laws considered in our model are issued from results of identification performed on the PA12.

In [START_REF] Luisier | Time-temperaturetransformation diagram for reactive processing of polyamide 12[END_REF], the following polymerization kinetic has been identified (Fig. 1):

dα dt = K • exp -E act RT • α m • (1 -α) n (4)
Fig. 1 Evolution of the degree of polymerization with time calculated with Eq. 4 (see [START_REF] Luisier | Time-temperaturetransformation diagram for reactive processing of polyamide 12[END_REF])

with

K = 4, 574 • 10 6 (5)
This polymerization process affects the temperature evolution by introducing a source term in the energy balance (heat equation):

ρc p (T ) DT dt = k(T ) T + ρw H ∂α ∂t with D• dt = ∂• ∂t + v • ∇• ( 6
)
where ρ is the density of the resin, c p and k respectively its specific heat and its thermal conductivity, w the resin volume fraction and H the total heat generated in the reaction. A classical linear evolution of c p and k with the temperature is used:

c p (T ) = a c p • T + b c p k(T ) = a k • T + b k (7)
Concerning the relation between the viscosity and the degree of polymerization, the following law has been identified (also in [START_REF] Luisier | Time-temperaturetransformation diagram for reactive processing of polyamide 12[END_REF], Fig. 2):

η(T , α) = 2, 14 • exp (-0, 028.T + (12, 08 + 0, 016T ) • α) (8)
It can be noticed that an increase of the temperature leads to a decrease of the viscosity whereas the evolution of the polymerization degree increases the viscosity.

Let us now combine this two models to develop the coupled thermo-chemico-mechanical model of he reactive RTM injection process.

Injection process modelling

The resin impregnation is modelled as a flow in a porous media in which we introduce the different couplings Fig. 2 Evolution of the viscosity with the polymerization degree calculated with Eq. 8 described previously. The model equations which describe this process in the domain are given by:

• the Darcy's law linking the pressure gradient and the velocity:

v = -(K/η) • ∇p (9)
• the incompressibility of the fluid:

∇v = 0 (10)
• the evolution of the resin fraction function I which is governed by the advection equation:

∂I ∂t + v • ∇I = 0 ( 11
)
• the transport and the evolution of the polymerization degree which is governed by the advection equation:

∂α ∂t + v • ∇α = K • exp -E act RT • α m • (1 -α) n (12)
• the transport and the evolution of the temperature which is governed by the advection equation:

ρc p (T ) ∂T ∂t + v • ∇T = k(T ) T + ρw H ∂α ∂t (13)
• the evolution of the viscosity which is governed by the chemo-rheological law defined previously:

η(T , α) = 2, 14 • exp (-0, 028.T + (12, 08 + 0, 016T ) • α) (14)
Concerning the initial and boundary conditions, we assume that:

• for the flow:

-the pressure gradient in the normal direction to the mould walls is null, which is equivalent to a null fluid flow through these walls; -the injection pressure or the flow rate is prescribed on the inflow boundary (injection nozzle) ∂ -; -the flow is saturated at the nozzle::

∀t, I (x ∈ ∂ -, t) = 1 (15)
-at time t = 0, the mould is empty:

∀x ∈ , I (x, t = 0) = 0 (16)
• for the temperature:

-the mould, and all the domain boundaries (excluding the injection nozzle) is at the imposed temperature T mould ; -the fluid is injected at the temperature T inj ;

• for the polymerization degree:

-no polymerization occurs before the resin enter the mould:

∀t, α(x ∈ ∂ -, t) = 0 (17)
We are now going to develop the numerical scheme used to solve the filling process model.

Numerical treatment

To perform the mould filling simulation, including the different physical/chemical couplings, we use the following four-step incremental explicit numerical algorithm. For each time step:

• step 1: compute the pressure field by solving the extended variational formulation on the whole domain [START_REF] Sanchez | Towards an efficient numerical treatment of the transport problems in the resin transfer molding simulation[END_REF]:

∇p * (K/μ) -1 ∇p • I + p * p • (1 -I ) d = 0 (18)
• step 2: compute the resin velocity at each element by applying Darcy's law; • step 3: update the function I by solving the advection equation governing its evolution; • step 4: update all the other transported quantities as the temperature or the polymerisation degree, taking into account the different couplings.

In what follows and for the sake of simplicity, we do not develop the two first steps as we use a classical finite element method to compute the pressure and the calculation of the pressure gradient at each element.

Simulation of the filling process: updating the fluid fraction I

The entire development of the advection equation integration is developed in [START_REF] Sanchez | Towards an efficient numerical treatment of the transport problems in the resin transfer molding simulation[END_REF]. Then, here we just recall the main features of the technique.

The evolution of the fluid fraction function I is given by the advection equation:

∂I ∂t + v • ∇I = 0 (19)
The initial and boundary conditions are given respectively by Eqs. 15 and 16. We now discretize the whole domain into N nonoverlapping cells or elements (triangles in our computer implementation). We integrate the weak formulation of Eq. 19 in each element e, use the divergence theorem and the incompressibility of the fluid and it results:

e ∂I ∂t d + ∂ e + I (v • n) ds + ∂ e - I (v • n) ds = 0 (20)
where ∂ e -denotes the inflow cell boundary and ∂ e + the outflow one, n is the unit outward vector defined on the cell boundary ∂ e .

To avoid problems related to the definition of function I in the element boundaries when using constant values of I in each element, we consider a piecewise constant approximation of I. Then, on the cell outflow boundaries, the value of I is taken equal to the one existing in the considered element e, and on the cell inflow boundaries, the value of I is taken equal to its value in the upstream element.

Finally, using a first-order explicit approximation of the time derivative, Eq. 20 can be written on each element as:

From a physical point of view these parameters ensure that an element must be almost fully filled (I e 1) before the downstream one begins to fill.

Transport and coupling: updating the other fields

The equations governing the transport and the evolution of the other quantities, namely the temperature and the polymerization degree, are given by the general linear advection equation:

∂J ∂t + v • ∇J = S ( 25
)
where J stands for the transported quantities (T or α) and S for the different source terms introduced by the sources and the different couplings.

As for the fluid fraction, we can develop the discretization of Eq. 25 that leads: However, it has been demonstrated in [START_REF] Sanchez | Towards an efficient numerical treatment of the transport problems in the resin transfer molding simulation[END_REF] that this equation involves some difficulties: we must use specific time and space discretization to avoid non-physical situations and the issue related to the definition of the quantities in the empty part of the mould where the quantities are not defined physically remains.

Indeed, lets consider an element that just starts its filling stage. If we assume that J is zero in the empty cells (hypothesis commonly done to initialize the different fields in the empty domain), then when the element e starts its filling process Eq. 26 gives:

J n+1 e = δ e-J n e- e- | e | + S n+1 t (27)
This equation shows that, for e-small enough (value that depend on the imposed flow rate and the time and space discretizations), we will have a value of J related to the fluid just introduced in element e smaller than the one associated with the upstream fluid. This situation is totally unphysical and, moreover, it shows that the initialization of the quantity in the element that just begins its filling process determines the corresponding field evolution.

To overcome these problems, the strategy proposed in [START_REF] Chinesta | Some difficulties in the flow front treatment in fixed mesh simulations of composites forming processes[END_REF][START_REF] Sanchez | Towards an efficient numerical treatment of the transport problems in the resin transfer molding simulation[END_REF] consists in considering the quantity I • J instead of J, that allows rewrite Eq. 26 as: Now, if we consider again the element that just begins its filling process at time t n+1 and that consequently has J n e = 0, then by combining the Eqs. 24 and 28, we obtain:

J n+1 e = J n e-+ S n+1 t (29)
We here ensure a right physical updating. And more, if we consider a pure convection situation (S n+1 = 0), we prove the exact transport:

J n+1 e = J n e- ( 30
)
These scheme and equations have been implemented in a numerical code than we are now applying for simulating some reactive RTM scenarios.

Numerical results

In this part, we first consider simulations to demonstrate the ability of the model to describe the process, and then, the influence of the introduction of the couplings in the mould filling evolution.

Simulation of reactive RTM process

We consider a rectangular 0, 4 × 1, 5 m mould and a fibrous preshaped reinforcement with permeability (homogeneous and isotropic) K. The injection is performed with a constant flow rate V inj . The value of all the parameters used in the simulations are reported in Table 1.

It can be noticed that in our simulations the mould is heated at the injection temperature. In future works, depending on the noticed influence of the injection temperature, alternatives could be considered in order to optimize the processing efficiency.

The simulations presented hereafter aim at evidencing the importance of taking into account the different couplings occurring in the mould during the injection stage, and in particular the coupling between the viscosity and the polymerization degree.

To do so, we perform a first simulation without introducing the coupling between the viscosity and the couple polymerization degree/temperature. For a total time equal to 8 s, Fig. 3 shows the prediction of the filling state and the polymerization degree of the resin in the mould.

It can be noticed that in this case the mould is almost filled, and if we let the calculation going forward, the mould ends up to fill with a complete polymerization of the entire volume of resin. In order to verify the simulation code we check that the filling time coincides with the theoretical value | |/V inj .

We now introduce the chemo-rheological law giving the evolution of the viscosity with the polymerization degree/ temperature and report the prediction in the same conditions in Fig. 4.

It can be notice that the filling process has been completely modified because the polymerization of the resin and the increase of viscosity associated leads to the formation of a blockage that avoids the complete filling of the mould.

Results prove that the couplings have to be introduced in the simulation in order to be able to predict the resin state during the mould filling stage. They also show that the numerical tool is a precious tool in order to define the conditions that allow the complete mould filling. But to do so, the identified thermo-chemical laws have to be precise enough. Indeed, in our example we do not analyze the effect of polymerization on the temperature evolution because the law considered here and that was identified in [START_REF] Luisier | Time-temperaturetransformation diagram for reactive processing of polyamide 12[END_REF] allows that assumption.

However, these first simulations allow to highlight some specific numerical features that we are addressing in what follows.

Numerical features related to the numerical integration of the transport equations

The numerical features that we are going to address now concern the integration of the quantities transported and in particular the treatment of the convection term. In the following, we will focus on pure convection so we consider that S n+1 = 0 and then:

Let us consider now and element that is being filled. At each time step, this element receives an amount of fluid, and then its fluid fraction increases of I depending on the flow rate and the time step.

Because we introduce the parameters δ e to avoid numerical diffusion, the fluid fraction in the element has to reach this value before it can transfer the fluid that it contains to the elements located downstream.

Then, if the fluid fraction is such that: . This error comes from the fact that the inflow volume of fluid coming to element e in the interval [t n , t n+1 ] is not compensated by the flow leaving the element because, as I n e < δ e , the fluid in element e cannot be transferred to the neighbour elements located downstream.

I n e <
This issue remains in the integration of the equation governing the transport of all the other quantities. For the quantity I, this numerical error is corrected in the code by enforcing I = 1 when the computed fluid fraction becomes greater than 1. However, as we can't impose a maximum value for some quantities (temperature, residence time, ..), the integration process will introduce some inevitable errors in the field updating.

For the polymerization degree, this error does not affect too much the calculation because its maximum value is one. But, for the temperature, this error can become more critical because the temperature order of magnitude is of hundred and it has not a maximum physical value.

Indeed, errors in the temperature updating lead to a nonnegligible spurious increase of the temperature in the resin. As our polymerization law is very sensitive to the temperature, this error may lead to a misrepresentation of the real process simulation.

One solution to avoid this problem would be to consider lower values of δ e . Then, whatever the elementary fluid fraction I , the element will transfer fluid to the downstream elements before overpassing the value I = 1. However such a choice has an incidence on the flow front thickness, that spreads artificially through some elements, and even if in reality such a flow front representation could be quite realistic, an appropriate physically based choice of δ e should be proposed. By this way, we ensure no excess. However, as pointed out in the last section, changing the time step depending of the process conditions (e.g. injection flow rate) makes difficult the use of standard model reduction strategies. Then, finally, other solution to avoid troubles is to use very small time steps when integrating the transport equations. And, in order to prevent too large calculation costs, an idea would be to use different time steps depending on the equations solved:

-for the integration of the Darcy's law over the whole domain, a quite large time step can be used; -lower time steps should be considered for the integration of the transport equations.

Now that we have shown the capabilities of our numerical strategy and underlined the tricky numerical points, in the following we try to go further and envisage the possibility to use this model to develop a functional software able to optimize the reactive RTM process.

Toward efficient RTM process optimisation

The objective from the point of view of process optimization is to define correct parameters that ensure a complete mould filling. In our case, the two main process parameters that play a key role in the filling process and that can be easily controlled, are the flow rate imposed at the injection nozzle and the injection temperature (which in our case coincides with the one applied on the mould).

Indeed, both parameters influence the filling:

-a too low flow rate can lead to an incomplete filling of the mould due to the formation of blockages created by the increase of the resin viscosity during polymerization; -on the opposite, a too high flow rate, can originate voids formation due to a bad impregnation of the preform (Fig. 5); -a high temperature helps in maintaining a low viscosity of the resin but increase the processing cost if the mould has to be maintained at the same temperature.

The idea is then to determine the optimal process parameters (V inj , T inj ) that ensure a complete filling of the mould without implying too high processing costs (a too high mould temperature or injection velocity).

Thus, we propose to create a virtual chart, having V inj and T inj as variables so that, depending on the material used and the mould size, the operator can operate within a processability window that ensure a correct filling of the mould. Fig. 6 SVD method applied to the separation time/space and T inj of the calculated polymerization degree α Such work has already been performed for forming process optimization in [START_REF] Ghnatios | Methodological approach to efficient modeling and optimization of thermal processes taking place in a die: application to pultrusion[END_REF]. The idea is to use the PGD method [START_REF] Chinesta | A short review on model order reduction based on proper generalized decomposition[END_REF] to solve a multiparametric problem where process parameters are considered as extra-coordinates. In the following, the objective is to determine whether such method can be applied for the reactive RTM process underlying the key point that will have to be addressed to reach this objective.

To do so, we have to check whether the model solution can easily be expressed in a separated form.

We performed several simulations, with different values of the parameters (V inj , T inj ) and determine if these solutions can be written as a finite sums decompositions involving functions of space and time, functions depending on the injection rate and finally functions involving the mould temperature.

Separation of the temperature T inj First, we focus on the mould and injection temperature. Several simulations are performed with values of the mould temperature from 200 • C to 500 • C. In order to check the possible separated form, a SVD is performed on the solution. As space and time decompositions in evolving domain problems involve too many terms, we decided to check the separated representation that follows:

J (x, t, T mould ) = F i (x, t) • G i (T mould ) (35)
The result of the SVD applied on the evolution of the polymerization degree is presented on Fig. 6. It can be notice that the solution can be represented with a few number of space/time and temperature modes.

However, one has to keep in mind that this test has been performed with a mould temperature equal to the injection temperature. In this case, even if the polymerization process introduces a small temperature variation in the resin, one can consider that the temperature is quite homogeneous in the whole volume of resin and so, the separability space/time and temperature is better.

Separation of the injection flow rate V inj Now we analyse the possible separation of the solution in space/time and injection flow rate modes:

J (x, t, V inj ) = F i (x, t) • G i (V inj ) ( 36
)
In this case, the separability issue is trickier as the injection flow rate has a huge influence on the numerical integration of the model. Indeed, as said in the previous part, this parameter play a crucial role in the calculation and the integration of the quantities transported during the filling simulation.

The dependency of the injection time on the injection flow rate makes also difficult the separability.

Then in order to be able to develop virtual charts to be used as process optimization support, a crucial point will be to focus on the V inj influence to make the separability possible. These points constitute a work in progress.

Conclusion

In this work, a complete model has been developed in order to take into account the different thermo-chemicomechanical couplings that occurs during the mould filling in reactive RTM processes. It includes the evolution of the temperature with the polymerization but more importantly, the evolution of the resin viscosity with the polymerization degree. This last phenomenon is crucial in the simulations as it allows to check the complete mould filling or the blockage originated by a too high increase of the viscosity.

A numerical strategy has been then proposed to integrate the model equations in order to develop a code able to perform reactive RTM process simulations. It uses a volume-of-fluid scheme and introduces a fluid fraction function I that allows to avoid the issue related to the initialization of variables associated with the fluid in the empty regions of the mould. A particular focus has been done on the discretization which is a tricky point in order to avoid numerical errors. Some results of reactive RTM process simulations have then been presented to illustrate the capability of the model and the numerical tool. Finally, a discussion has been addressed on the possibility of using model order reduction based on a separated representation of the multiparametric solution. This question requires a deeper analysis that will be addressed in future works.

  refers to the considered element and the superscript to the time step, | e | denotes the surface of the element e, I e-defines the value of I in the upstream neighbour, e+ and e-are the inflow and outflow fluid volumes defined by: scheme can introduce too diffuse flow fronts and certain instabilities related to small values of I n e . To overcome theses effects, we introduce the parameters δ e and δ e-defined as:δ e =1if I e ≥ I th 0 if I e < I th and δ e-= 1 if I e-≥ I th 0 if I e-< I th (

Fig. 3

 3 Fig. 3 Fluid fraction function and polymerization state of the resin at t = 8 s for a simulation without chemo-rheological coupling

  time n + 1, we will have an extra unphysical fluid fraction in the considered element (because physically the fluid fraction cannot be greater than 1) whose value is equal to 1 -I n+1 e

Fig. 4 Fig. 5

 45 Fig. 4 Fluid fraction function and polymerization state of the resin at t = 8 s for a simulation with chemo-rheological coupling

Table 1

 1 Parameters values used to performed the simulationsChemical parameters (with respect to Eq. 4)

	E act (kJ/mol)	m	n	R	H (J/g)
	80.2	0.3	1.05	8.314	-53
	Resin properties (with respect to Eq. 8)
	ρ (kg/m 3 )	w	a k	b k	a cp	b cp
	1560	0,4	0,5836	-0, 0012	1205,3	1,6
	Process parameters
	Preform permeability ()	V inj (m.s -1 )	T inj (K)	T mould (K)
	10 -7	0,05	493	493

Acknowledgments This work has been partially supported by the IUF-Institut Universitaire de France-and the french ANR RTMPlast project.


R Ata

A Soulaïmani

F Chinesta

The natural volume method (NVM): Presentation and application to shallow water inviscid flows

Keywords: NEM, shallow water equations, stabilization, Lagrangian description, collocation, inviscid flows 1

In this paper a fully Lagrangian formulation is used to simulate 2D shallow water inviscid flows. The natural element method (NEM), which has been used successfully with several solid and fluid mechanics applications, is used to approximate the fluxes over Voronoi cells. This particle-based method has shown huge potential in terms of handling problems involving large deformations. Its main advantage lies in the interpolant character of its shape function and consequently the ease it allows with respect to the imposition of Dirichlet boundary conditions.

In this paper, we use the NEM collocationally, and in a Lagrangian kinematic description, in order to simulate shallow water flows that are boundary moving problems. This formulation is ultimately shown to constitute a finite-volume methodology requiring a flux computation on Voronoi cells rather than the standard elements, in a triangular or quadrilateral mesh. St Venant equations are used as the mathematical model. These equations have discontinuous solutions that physically represent the existence of shock waves, meaning that stabilization issues have thus been considered. An artificial viscosity deduced from an analogy with Riemann solvers is introduced to upwind the scheme and therefore stabilize the method. Some inviscid bidimensional flows were used as preliminary benchmark tests, which produced decent results, leading to well-founded hopes for the future of this method in real applications.

INTRODUCTION

Free surface flows are often simulated in shallow domains using Eulerian descriptions. Most software applications use the finite-element or the finite-volume method in Eulerian formulations. However, several difficulties arise in simulating some phenomena related to the Lagrangian behavior of the flow, such as in moving boundary problems or moving shocks. The Lagrangian description is thus preferred for overcoming most such difficulties. Fully Lagrangian applications are few in the literature and we can mostly find applications that use the so-called moving mesh methods. The improvement in accuracy is achieved by adapting a fixed mesh so that nodes are concentrated in the areas where the solution presents high gradients (such as shocks or boundary layer areas). The mesh is thus moved using equations that involve node velocities. Normally, the partial differential equation (PDE) and the mesh equations are solved simultaneously. Examples of moving mesh methods include the moving finite-difference methods (for instance, see Dorfi and Drury [START_REF] Dorfi | Simple adaptive grids for 1-D initial value problems[END_REF]) and the moving finite-element method (FEM) [START_REF] Miller | Moving finite elements I[END_REF][START_REF] Miller | Moving finite elements II[END_REF]. More recently, the meshless FEM was introduced by Idelsohn et al. in [START_REF] Idelsohn | The meshless finite element method[END_REF] and applied in fluid mechanics with free surface in [START_REF] Oñate | The particle finite element method. An overview[END_REF]. In the same family of methods, Huang and Russell [START_REF] Huang | A moving collocation method for solving time dependent partial differential equations[END_REF] used a cell averaging cubic Hermite collocation discretization for the physical PDEs combined with a finite-difference discretization for the PDE to solve problems with discontinuous solutions.

Most works covering these methods employ conventional numerical approaches. Actually, in the last decades, meshless-or mesh-free methods have surfaced as a serious mean of handling problems that could not be handled conventionally; these include problems presenting large distortion and fragmentation, and problems with unconventional behaviors. These methods are used mostly in the Galerkin formulations and some of them in a collocation manner. The family of mesh-free methods are used on a wide variety of mechanical and electrical applications. In hydraulics, Monaghan has explored the smooth particle hydrodynamics (SPH) method for free surface flows using an isentropic state law [START_REF] Monaghan | Simulating free surface flows with SPH[END_REF]. Many others have used the same law to simulate several problems in fluid mechanics, including Bonet et al. [START_REF] Bonet | Variational formulation for the smooth particle hydrodynamics (SPH) simulation of fluid and solid problems[END_REF] who covered the use of SPH in fluids. The volume of fluid technique applied to the Navier-Stokes equations gave spectacular results when simulating free surface flows with breaking waves and for sloshing problems [START_REF] Lin | A numerical study of breaking waves in the surf zone[END_REF]. Ata and Soulaïmani presented a stabilized SPH method for inviscid shallow water flows using the St Venant equations [START_REF] Ata | A stabilized SPH method for inviscid shallow water flows[END_REF], where they showed that the stabilization problem can be solved using an artificial viscosity arising from the upwinding of the flux in Riemann solvers. The treatment of boundary conditions, nonetheless, remains a problem that is yet to be solved. The non-interpolant nature of the shape function of almost all the mesh-free methods, including SPH, represents the main handicap when imposing Dirichlet boundary conditions. These methods are therefore only well adapted to simulating unboundeddomain problems in fluid mechanics. The natural element method (NEM) represents a serious mean of overcoming such problems; the NEM has the advantage of having an interpolant shape function.

The NEM was originally a Galerkin-based method. It is considered to fall under the meshless method family even though it uses the Delaunay triangulation or the Voronoi diagram to build its shape function. The NEM uses the notion of natural neighbor interpolation, an interpolation that has very striking properties such as a strictly interpolating character, the ability to exactly interpolate piece-wise linear boundary conditions on convex domain, and a well-defined and robust approximation with no user-defined parameters on non-uniform grids [START_REF] Cueto | Overview and recent advances in natural neighbor Galerkin methods[END_REF]. For a non-convex domain, the NEM does not reproduce linear boundary conditions. A counterexample is shown in [START_REF] Cueto | On the imposition of essential boundary conditions in natural neighbor Galerkin methods[END_REF]. This problem could be overcome by the use of -shape technique [START_REF] Cueto | Overview and recent advances in natural neighbor Galerkin methods[END_REF] or a constrained Voronoi diagram [START_REF] Yvonnet | Interpolation naturelle sur les domaines non-convexes par l'utilisation du diagrame de Voronoi contraint[END_REF].

This method was first used for the interpolation of data from a set of scattered points [START_REF] Sibson | A brief description of natural neighbor interpolation[END_REF]. Recently, it has been used for some applications in solid mechanics [START_REF] Sukumar | The natural elements method in solid mechanics[END_REF], such as those involving multiscale problems [START_REF] Missoum-Benziane | A new fully coupled two-scales modelling for mechanical problems involving microstructure: the 95/5 technique[END_REF] and crack simulation [START_REF] Yvonnet | Interpolation naturelle sur les domaines non-convexes par l'utilisation du diagrame de Voronoi contraint[END_REF]. Liang et al. [START_REF] Liang | Simulation of dam-and dyke-break hydrodynamics on dynamically adaptive quadtree grids[END_REF] used the Sibson NEM for the interpolation of data on new inserted nodes when applying the finite-volume method (FVM) on a dynamically adaptive quadtree grid. Del Pin [START_REF] Pin | A Lagrangian approach based on the natural neighbor interpolation[END_REF] used the NEM shape function as a kernel to implement the moving-particle semi-implicit method in a Lagrangian formulation; the latter is a moving finite-difference method.

In this paper, the NEM is investigated in order to explore its capacity to simulate free surface flows in a fully Lagrangian formulation. The NEM's shape function is defined on Voronoi cells, giving rise to a non-standard Lagrangian finite-volume scheme called the natural volume method (NVM).

We begin by giving a brief description of the NEM in Section 2 and we recall the C-NEM method in Section 3. In Section 4, we present the shallow water equations (SWE) and their discretization in the NVM framework. Finally, results on 2D applications are shown.

THE NATURAL ELEMENT METHOD

In the NEM, interpolation functions (the trial and test functions for Galerkin formulation) are constructed using the notion of natural neighbors. Natural neighbor co-ordinates are based on wellknown geometric concepts such as the Voronoi diagram and the Delaunay tessellation [START_REF] Sukumar | Natural neighbour Galerkin methods[END_REF]. It was Sibson [START_REF] Sibson | A brief description of natural neighbor interpolation[END_REF] who first presented the notion of natural neighbor interpolation for data approximation. Further works on their construction and properties were carried out later by Farin [START_REF] Farin | Surfaces over Dirichlet tesselations[END_REF] and others. Recently, a new type of interpolation (non-Sibsonian), which is similar to but less costly than Sibson's, has been introduced by Belikov et al. [START_REF] Belikov | The non-Sibsonian interpolation: a new method for interpolation of the values of a function on an arbitrary set of points[END_REF]. In fact, the non-Sibsonian interpolant requires the computation of Lebesgue measures of order D -1 in R D to compute the shape functions, whereas for the Sibson interpolant, a computation of areas (in 2D) or volumes (in 3D) is necessary.

Voronoi diagram and Delaunay tessellation

The Voronoi (or Dirichlet) diagram is a unique geometric construct obtained from a set of distinct nodes in a Euclidian space of dimension D.

Let be a set of M distinct nodes of a bounded domain : = {n 1 , n 2 , . . . , n M }. The Voronoi diagram associated with the set is the subdivision of the domain into regions V (n I ) such that each point in V (n I ) is the closest to node n I than to any other node n J ∈ (I = J ). The region V (n I ) (first-order Voronoi cell) of node n I is the convex polygon (polyhedron):

V (n I ) = {x ∈ R D : d(x, x I )<d(x, x J ) ∀J = I } (1)
where d is a Lebesgue measure defined on . The dual of the Voronoi diagram is the Delaunay tessellation. The Delaunay triangles are constructed by connecting nodes whose Voronoi cells have common boundaries. It is important to note that the Voronoi diagram is unique, even if the Delaunay tessellation is not. Figure 1 shows the triangulation and its dual Voronoi diagram for the case of 9 nodes scattered in a square in 2D.

For the sake of simplicity, from now on we will consider the space dimension D = 2.

Sibsonian interpolants [14]

The definition of the first-order Voronoi diagram was introduced by Equation ( 1). Similarly, we can extend this definition to a higher order k>1. Then, we define the Voronoi diagram of second order as

V I J = {x ∈ R 2 : d(x, x I )<d(x, x J )<d(x, x K ) ∀K = I, J } (2)
The second-order Voronoi diagram then becomes the subdivision of the plane (or the space) into domains V I J such that V I J is the locus of all points that have I as the closest node and J as the second closest node.

Sibson has used the second-order Voronoi diagram notion to define the natural neighbor-type interpolation function. This function can be used as a trial or test function in Galerkin-type formulations.

If a point P with coordinates x is to be inserted into the Delaunay tessellation, natural neighbortype shape function associated with P is defined as the ratio of the area of the second-order Voronoi cell A I (P) to the total area of the first-order Voronoi cell A(P) associated with P (see Figure 2):

I (x) = A I (P) A(P) , A(P) = n J =1 A J (P) ( 3
)
where n is the number of natural neighbors of P. The derivatives of the Sibson-type shape function could be obtained by deriving Equation (3) as follows:

I, j (x) = A I, j (P)-I (x) A , j (P) A(P) , j = 1, 2 (4)
For more details and to consult exhaustive studies on all characteristics of the Sibson interpolation, we refer to References [START_REF] Cueto | Overview and recent advances in natural neighbor Galerkin methods[END_REF][START_REF] Yvonnet | A new extension of the natural element method for nonconvex and discontinuous problems: the constrained natural element method (C-NEM)[END_REF].

Non-Sibsonian or Laplace interpolation

Let be a set of nodes to which are associated their Voronoi cells defined by Equation (1). Let t I J be the facet associated with cells V I and V J (of dimension D -1: segment in 2D or polygon in 3D) and d(t I J ) the measure of Lebesgue of t I J . If any point P with coordinates x ∈ R 2 is introduced into the tessellation and if this point has n natural neighbors, then the Laplace shape function associated with the node I , in the 2D case, is [START_REF] Sukumar | Natural neighbour Galerkin methods[END_REF][START_REF] Belikov | The non-Sibsonian interpolation: a new method for interpolation of the values of a function on an arbitrary set of points[END_REF][START_REF] Christ | Weights of links and plaquettes in a random lattice[END_REF]

I (x) = I (x) n J =1 J (x) , J (x) = s J (x) h J (x) (5)
where J (x) is the Laplace weight function, s I (x) is the length of the Voronoi edge associated with P and node I , and h I (x) is the Euclidian distance between points P and I (see Figure 2). The derivative of the Laplace interpolation function is

I, j (x) = I, j (x)-I (x) , j (x) (x) (6)
with (x) = J J (x). In 2D, the Laplace shape function uses a fraction of length whereas Sibson's function uses a fraction of areas, and as a result the computational cost is lower for the Laplace interpolation. This advantage is more remarkable in the 3D case.

Properties

In addition to positivity, the most important properties of the natural element-type interpolation given by Equation (3) or ( 5) are listed below

• The interpolation is defined from the nodal values according to

u h (x) = n I I (x)u I (7)
• Nodal interpolation: This can be expressed as

I (x J ) = I J (8)
with I J being the delta Kronecker function. Thus, as for the cases of FEM or FVM, the NEM is an interpolant method. This property is very important as it simplifies the imposition of Dirichlet boundary conditions. This property is lost for most of the meshless method family, which therefore represents a major drawback for them. • The partition of unity: Sibson and Laplace interpolations possess zero-order consistency or the partition of unity property

N I I (x) = 1 (9)
Consequently, all NEM interpolations reproduce exactly constant functions. • Linear completeness: The NEM-type interpolations reproduce exactly linear functions, which yield their linear completeness

x = N I I (x)x I ( 10
)
See the demonstration in [START_REF] Cueto | Overview and recent advances in natural neighbor Galerkin methods[END_REF]. • Strict linearity on the border: At the boundary of a convex domain, the shape function is strictly linear between two successive nodes; a proof can be found in [START_REF] Sukumar | The natural elements method in solid mechanics[END_REF]. • The support nature: The support of NEM-type interpolation functions is the union of Delaunay circumcircles of node I . Figure 3 presents a comparison between the NEM support and the moving least-squares (MLS) support [START_REF] Sukumar | Natural neighbour Galerkin methods[END_REF]. For Galerkin-type formulations, this shape of the support makes the numerical integration a relatively tough task when compared with finite elements, for instance. Besides, another difficulty lies on the fact that the shape function is rational, which complicates the integration task more. For our work, the collocation formulation avoids all these drawbacks; for further details, see [START_REF] González | Numerical integration in natural neighbour Galerkin methods[END_REF]. • In 2D, the implementation of Bowyer-Watson algorithm [START_REF] Bowyer | Computing Dirichlet tesselations[END_REF][START_REF] Watson | Computing the n-dimensional Delaunay tesselation with application to Voronoi polytopes[END_REF] can be adopted to compute the Sibson-and Laplace-type shape functions. The Lasserre algorithm has been used also in [START_REF] Sambridge | Geophysical parametrization and interpolation of irregular data using natural neighbours[END_REF]. The Watson algorithm is not extensible in 3D. Besides, it fails on the nodes located on the Delaunay edges. Thus, the Bowyer-Watson algorithm is efficient to evaluate the Laplace shape function in 2D [START_REF] Cueto | Overview and recent advances in natural neighbor Galerkin methods[END_REF]. For a full description and for more details on this algorithm, we refer, for instance, to Cueto et al. [START_REF] Cueto | Overview and recent advances in natural neighbor Galerkin methods[END_REF].

THE CONSTRAINED NEM OR C-NEM

To overcome the problems seen in the NEM in its standard formulation when handling the nonconvex character of a domain, Cueto et al. [START_REF] Cueto | Imposing essential boundary conditions in the natural elements method by mean of density-scaled alpha-shapes[END_REF] first proposed a new technique called 'alpha shape,' which is a generalization of the concept of the convex hull of a cloud of points. Yvonnet et al. [START_REF] Yvonnet | A new extension of the natural element method for nonconvex and discontinuous problems: the constrained natural element method (C-NEM)[END_REF] suggested an original approach based on the use of the constrained Voronoi diagram, which was introduced by Klein and Lingas [START_REF] Klein | A linear-time randomized algorithm for the bounded Voronoi diagram of a simple polygon[END_REF]. Giving a constrained Voronoi diagram, the C-NEM-type interpolation is given as: An unknown function u is approximated by u h (x) given by

u h (x) = N I =1 C I (x)u I ( 11
)
where N is the number of natural neighbors visible from x and C I are the constrained natural neighbor shape functions. u I are the nodal unknown values of natural neighbors visible from the point x [START_REF] Yvonnet | A new extension of the natural element method for nonconvex and discontinuous problems: the constrained natural element method (C-NEM)[END_REF].

It is very important to note that when a point x comes close to an edge or a face of , the contour of the concerned Voronoi cells is extended to infinity to the outside of such that the contribution in area (in 2D) or in volume (in 3D) is infinite. These infinite contributions allow the retrieval of the linearity of the interpolation on the boundary of the domain. Figure 4 shows the C-NEM-type shape functions (Sibson and Laplace).

PROBLEM SETTING

The Shallow water equations (SWE)

We attempt to find a numerical solution to the SWE using the NEM. The NEM was used to simulate some fluid mechanics problems especially relating to metal forming. The updated Lagrangian formulation was mostly used [START_REF] Martinez | Updated Lagrangian free surface flow simulations with natural neighbour Galerkin methods[END_REF]. That is the case of the FEM when applied to treat SWE [START_REF] Soulaïmani | An arbitrary Lagrangian-Eulerian finite element method for solving three dimensional free surface flows[END_REF].

González et al. [START_REF] González | A natural element updated Lagrangian strategy for free-surface fluid dynamics[END_REF] recently presented an updated Lagrangian application to some free surface problems in fluid dynamics using Navier-stokes equations with moderate Reynolds numbers. No fully Lagrangian application of the NEM seems to be available in the literature. The rapid distortion of used mesh, and thus the frequent re-meshing, represents the biggest challenge to overcome for scientists. The paper by Birknes and Pedersen [START_REF] Birknes | A particle finite element method applied to long wave run-up[END_REF] could be taken as a typical work. It introduces the use of the 'particle FEM (PFEM)' to the SWE to simulate the long wave run-up. The originality of this work lies in its use of a mixed Eulerian-Lagrangian formulation of the PFEM. The inviscid SWE in the non-conservative form when neglecting the bed slope and friction terms are given by

Dh Dt +h∇ •u = 0 (12) Du Dt + g∇h = 0 (13)
where h, u, and g are, respectively, water height, depth-average velocity, and gravity. D/Dt refers to the total derivative. Cases that introduce bed slope and friction terms will be covered in future works.

Mass conservation

Usually, continuity equation ( 12) is used for mass conservation. However, although we use a Lagrangian formulation of the NEM, the mass is automatically conserved. In fact, the total mass of the fluid is divided into N particles, where N is the total number of nodes. Each particle is assigned to a node such that, in any given time step, the mass (the volume for an incompressible medium) of each particle is conserved. The volume of a particle I is given by V I = A I h I , where A I is the Voronoi cell area and h I is the averaged water depth over the area. The new value of the depth h n+1 I at the time t n+1 = t n + t is computed as follows:

V n+1 I = h n+1 I A n+1 I = V n I = h n I A n I ⇒ h n+1 I = h n I A n I A n+1 I ( 14
)
A n+1 I
is the area of the deformed cell associated with the node I at the time t n+1 . A n+1 I is obtained after time integration of the coordinates of all the vertices of the Voronoi cell associated with node I . It is important to note that the geometric Voronoi cell, i.e. the cell related to node I corresponding to the dual of the Delaunay triangulation, does not correspond to the one we use to compute the new water depth. In fact, to ensure the local mass conservation, we have to use the initial Voronoi diagram (at t = 0 s) and update its coordinates at each time step (Figure 5).

The coordinates updating of the Voronoi diagram needs the computation of the acceleration and the velocity values at the vertices of each cell. Such computations require that first we stock the coordinates and the old values of acceleration and velocity. Furthermore, we need to 'interpolate' such functions for each time step. Interpolation [START_REF] Cueto | Overview and recent advances in natural neighbor Galerkin methods[END_REF] can be used to carry out that purpose. For instance, to compute the velocity at a vertex s of the cell of a node I , s is first introduced in the tessellation, and the velocity is given by

u h (s) = N J =1 C J (s)u J ( 15
)
where J = 1, . . . , N are the visible natural neighbors of s.

Conservation of momentum

We look for the application of the NEM in a collocation form to simulate shallow water flows using St Venant equations. Although, most NEM applications are in a Galerkin formulation, a finite-difference scheme of the method was first proposed by Sukumar in [START_REF] Sukumar | Voronoi cell finite difference method for the diffusion operator on arbitrary unconstructed grids[END_REF], and then discussed extensively in [START_REF] Cueto | Overview and recent advances in natural neighbor Galerkin methods[END_REF]. Unlike Sukumar and Cueto's use of the terminology of finite difference, we prefer to use the term of finite volumes while it consists of the use of averaged values of the flux on all Voronoi diagram edges. For a node I located at x, the discrete nodal gradient can be approximated directly by deriving Equation [START_REF] Oñate | The particle finite element method. An overview[END_REF] in the coordinate directions. This approach is adopted in the most meshless methods when used in a collocation manner. As a second way to approximate the nodal gradient, we can consider the finite-volume averaging for the strain operator, which was proposed by Chen et al. [START_REF] Chen | A stabilized conforming nodal integration for Galerkin meshfree methods[END_REF]. Generally, the jth component of the gradient of a field w is given by

*w h *x j I = A I w h , j (x) dA A I = I w h n j d A I ( 16
)
where w is the field to be derived, A I is the area, and I is the boundary of the Voronoi cell related to node I . n j is the jth component of the unit outward normal vector to I . The Gauss divergence theorem has been used in Equation ( 16).

For the special case of a regular grid, we can approximate the integral using an average weighted value of w h on each Voronoi edge [START_REF] Cueto | Overview and recent advances in natural neighbor Galerkin methods[END_REF]. The choice w h (x m ) = (w I +w J )/2 renders an exact result for a linear field w on a regular grid (x m is the middle point of the segment [I J]):

*w h *x j I = N J =1 (w I +w J )n I J j s I J 2A I , A I = 1 4 N J =1 s I J h I J ( 17
)
where n I J j is the jth component of the unit outward normal to the Voronoi edge common to nodes I and J and s I J is the length of this edge, and h I J is the length of the segment [I J]. We have a linear finite-element discretization here. In fact, the finite-element approximation is linear on the edges of the triangle and in the case of a regular nodal distribution, the middle point of the Voronoi edge relative to the segment [I J] is located on the edge of the triangle with I and J as vertices.

With the approximation of the gradient on hand, one can proceed to write a discrete version of the SWE. For a particle I , we apply Equation [START_REF] Liang | Simulation of dam-and dyke-break hydrodynamics on dynamically adaptive quadtree grids[END_REF] to evaluate the gradient of the pressure P = gh, following which Equation ( 13) is expressed as

Du I Dt = - 1 2A I N J =1 [(gh I + gh J )n J s I J ] ( 18
)
where u is the velocity field, g is the gravity acceleration, h I is the water depth of the particle I , and N is the number of edges of the cell related to the particle I . Note that N is not always equal to the number of natural neighbors of I , and that it is only in the case of a regular distribution of nodes and for an internal node (far from the boundary) that N corresponds to the number of visible natural neighbors. The approximation given by Equation ( 18) is the easiest and lowest cost choice. However, it is limited to cases where the nodal distribution is either mostly regular or where the solution is linear or quasi-linear. For more complicated cases, any other approximation type may be used (finite element, finite difference, natural element, or any Godunov-type scheme, etc.). For our case, the use of the Voronoi diagram favors the trivial choice of the natural element approximation.

As mentioned above, Equation (18) may only be used for cases with regular nodal distribution or simple geometries and for linear or quasi-linear solutions. In our case, the Lagrangian simulation induces a relative chaotic displacement of nodes. Therefore, the regularity of the initial nodal distribution is quickly lost. Figure 6 shows a very simple example where, for three nodes {1, 2, 3}, Equation ( 18) cannot be applied. In fact, for the Voronoi edge relative to nodes 1 and 2, the interpolation of any function w in the middle point of this edge using w h (x m ) = (w I +w J )/2 is far from accurate. To overcome this problem, an interpolation of the flux on the Voronoi edges is achieved using the NEM. The easiest way to quantify Equation ( 16), i.e. to compute the flux on the Voronoi edge, is to integrate this equation on a unique integration point (the middle point of the edge). For the case of a convection-diffusion problem under some regularity conditions, Mishev demonstrated, in [START_REF] Mishev | Finite volume methods on Voronoi meshes[END_REF], that this quadrature (as well as the one given by Equation ( 18)) ensures a O(h s-1 ) convergence, where s is the order of the Sobolev space (H s ) in which the solution is defined and regular (see the work of Mishev [START_REF] Mishev | Finite volume methods on Voronoi meshes[END_REF] for a full mathematical study of this quadrature scheme)

*w h *x j I = A I w h , j (x) d A I = *A I w h n j d A I = 1 A I N J ( wI J s J n I J j ) ( 19
)
where wI J is the NEM approximation of w estimated at M J , the middle point of the edge, s J is the edge length, and N is the number of edges of the Voronoi cell of I . wI J is obtained by applying Equation [START_REF] Cueto | Overview and recent advances in natural neighbor Galerkin methods[END_REF], which is a 'richer' approximation than the finite-element approximation used by Chen in Equation [START_REF] Liang | Simulation of dam-and dyke-break hydrodynamics on dynamically adaptive quadtree grids[END_REF] as it involves all the natural neighbors of M J rather than only I and J for the other approximations. Applying Equation ( 19) on the non-conservative form of the momentum equation of the SWE (Equation ( 13)) gives

u = Du I Dt = - 1 A I N J =1 (g h J )n I J s J ( 20
)
where h J = N J K =1 K (M J )h K and N J is the number of edges in the Voronoi cell related to M J . Using the identity given by n ds = 0, or in a discretized form N J =1 n I J s J = 0, Equation [START_REF] Farin | Surfaces over Dirichlet tesselations[END_REF] Figure 6. The case frequently obtained in Lagrangian simulations where Equation [START_REF] Pin | A Lagrangian approach based on the natural neighbor interpolation[END_REF] does not produce good approximation (especially for neighbor 2).

becomes

uI = Du I Dt = - 1 A I N J =1 (g h J + g h I )n I J s J (21)
Equation ( 21) is the final approximation of the momentum equation, which is preferred to Equation [START_REF] Pin | A Lagrangian approach based on the natural neighbor interpolation[END_REF]. This equation will be used for all tests that will be discussed in Section 10. However, owing to the symmetry of the scheme given by this equation, an upwinding technique must be used for stabilization purposes.

Remark

As described in the previous paragraphs, the Voronoi I ) that should be treated implicitly. However, an implicit treatment becomes very expensive in CPU time. For this reason, we opted for an explicit linearization, which assumes that h n+1 A n+1 h n A n .

CONSISTENCY AND STABILIZATION ISSUES

The study of the consistency of the method on non-uniform grid poses interesting challenges. The consistency basically ensures that in the limit, when the grid spacing tends to zero, the approximated solution tends to the analytical one. Using Taylor series, Cueto et al. [START_REF] Cueto | Overview and recent advances in natural neighbor Galerkin methods[END_REF] have demonstrated that using the finite-difference method, the approximation of the Laplacian operator is O(h m ) for a non-uniform grid, with h m being the mean of the internodal distance. At the very least, a first-order consistency is thus obtained on non-uniform grids. In the case of a regular grid, the consistency order raises to the second order. For more details, see the paper by Sukumar [START_REF] Sukumar | Voronoi cell finite difference method for the diffusion operator on arbitrary unconstructed grids[END_REF]. In our applications, the regularity of the nodal distribution cannot be conserved in the time because of the Lagrangian character of the simulations. That is why, we cannot guarantee more than a first-order consistency.

While the solutions for the SWE are discontinuous, shock waves are encountered when simulating shallow flows. On the other hand, Equations ( 18) and ( 21) offer a symmetric centereddifference-like scheme that is are mostly unstable. That is why a stabilization technique must be used in order to avoid spurious spatial oscillations that may grow in amplitude especially near discontinuities. As is the case in Ata and Soulaïmani [START_REF] Ata | A stabilized SPH method for inviscid shallow water flows[END_REF], we can adopt a stabilization technique by the analogy of what is often adopted with Riemann solvers. Doing so will lead to an upwinding of the scheme given by Equations ( 18) and ( 21) as it will introduce an artificial viscosity.

The quantity F h (x I ) = gh I can be considered as the flux through the Voronoi edge. Thus, Equation ( 18) could be expressed as

uI = Du Dt I = - 1 A I N J =1 (F h (x I )+ F h (x J )) 2 n I J s I J ( 22
)
Copyright q 2008 John Wiley & Sons, Ltd.

For the purpose of stabilization, the term (F h (x I )+ F h (x J ))/2 in Equation ( 22) is replaced by a simple Lax-Friedrich scheme given by

1 2 ((F h (x I )+ F h (x J ))-(u J -u I )•n I J ) ( 23
)
We can thus define an artificial viscosity that will stabilize the system and will especially not allow interpenetration between particles. This artificial viscosity is given by (we refer to [START_REF] Ata | A stabilized SPH method for inviscid shallow water flows[END_REF] for more details)

I J = -u I J •r I J r 2 I J + 2 = -[(v J -v I )(y J -y I )+(u J -u I )(x J -x I )] r 2 I J + 2 (24)
where

= 1 2 (c I +c J ) = 1 2 ( √ gh I + √ gh J )
is the averaged wave speed, u I J = u I -u J the difference between the velocities of I and J . u and v are the components of the velocity field in the x and y directions, respectively, r I J = r I -r J , is a small constant used to avoid the division by zero, and is a free parameter. The final momentum equation will be expressed as

uI = Du Dt I = - 1 2A I N J =1 ((gh I + gh J )+ I J )n I J s I J ( 25
)
The stabilization-upwinding term I J is introduced in the same manner as for Equation [START_REF] Bowyer | Computing Dirichlet tesselations[END_REF]. The final form of the discretized moment equation that we used is then

uI = Du Dt I = - 1 A I N J =1 (g h J + g h I + I J )n I J s J ( 26
)
Remarks 1. We are using two different cells: one for the mass and one for momentum conservation. The Voronoi cell is the one used for the momentum equation; the cell used for the computation of the water depth is a 'squeezed' Voronoi cell and then it does not verify any of the known properties of the Voronoi cell. 2. As for all Lagrangian formulations with large deformations, we could possibly see distorted cells appear. The distortion could appear for cells used for mass conservation (which are Voronoi cells at t = 0 s and only at t = 0 s). It cannot affect Voronoi cells which are always convex and which can never be distorted in a finite-element manner. In such a case, the reference configuration would be updated and even a mass redistribution could be required. For all the cases presented in Section 10, the cells did not attain a high level of distortion, and so we did not update any of them.

BOUNDARY CONDITIONS

The treatment of the boundary conditions is greatly simplified for the NEM as compared with other meshless methods. The interpolant character of the method simplifies the Dirichlet boundary imposition. For our cases, we need a perfect slipping condition on the wall, meaning that the normal component of the velocity and acceleration fields of particles situated on the boundary have to be zero u n I = 0 and un I = 0 ∀I on the domain boundary [START_REF] Sambridge | Geophysical parametrization and interpolation of irregular data using natural neighbours[END_REF] where u n I and un I are the normal components of the velocity and the acceleration of the particle I . To that end, velocity and acceleration vectors can be projected onto the tangent vector.

The Lagrangian character of the formulation induces several difficulties in the treatment of boundary conditions especially when the geometry is complex. In fact, the boundary of the domain (the geometry of the channel) is defined by a set of geometrical nodes. A fluid particle that moves at every time step is attached to every node. Therefore, the geometry could be lost or at least modified when the particles of the boundary move. Thus, to maintain a well-defined geometry, the boundary particles are fixed to the geometric locations of their corresponding nodes. This procedure must physically ensure a perfect slipping condition. To compute the value of acceleration, velocity, or water depth at the old position of the particle (the geometric position), a linear interpolation is done using the two closest particles. Figure 7 illustrates such a procedure. We must pay particular attention to a particle located in a corner or at an angular point in general. In this case, a direction is chosen (privileged) between the two directions that form the corner. However, the direction that is 'parallel' to the flow direction matches better with the flow's nature and naturally captures the physics of the problem. Such a problem is well known in all finite-volumes and finite-elements-type approaches. It is one of the points to be studied in further works.

We can also let boundary particles move freely at the edges of the domain. Such a procedure guarantees better results whereas it does not require accelerations' or velocities' interpolation. However, it causes more difficulties when following a particle over the boundary. The angular points and corners are handled in a similar manner as in the first procedure.

TIME INTEGRATION

The time integration is treated using an explicit scheme. This avoids the high resolution of systems but requires specific time step selection in order to stabilize time integration. In this work, we used an explicit Euler scheme in some cases, but mostly it was Newmark explicit scheme that we used. The velocity and the coordinates of a particle I are given by the following:

u n+1 I = u n I + t ((1-) un I + un+1 I ) (28)
x n+1 I = x n I + tu n I + t 2 (( 1 2 -) un I + un+1 I ) ( 29
)
where t is the time step, u is the velocity, u is the acceleration, and the notation u n means u(t n ).

By choosing the parameters and as, respectively, 1 2 and 1 4 , the scheme corresponds to the trapezoidal rule and is second-order accurate [START_REF] Bardella | Newmark's time integration method from the discretization of extended functionals[END_REF].

Loukili and Soulaïmani suggested in [START_REF] Loukili | Numerical tracking of shallow water waves by the unstructured finite volume WAF approximation[END_REF] that according to a linearized stability analysis applied to the SWE, the following Courant-Friedrichs-Lewy (CFL) restriction must be respected:

CFL = t max I ( √ gh + u 2 x +u 2 y ) min(d I J ) <1 ∀I, J (30)
For most finite-volume applications in Eulerian formulation, the CFL value is often greater than 0.5. In our cases, when particles are moving, and the distances between them change with every time step, we reduce the CFL values to a maximum of 0.5.

ADAPTIVITY AND NODES REFINEMENT

Given the fact that the computational procedure is fully Lagrangian, the nodes' distribution changes with every time step. Therefore, in certain cases, we could obtain strongly irregular particle concentration in which a big number of nodes is located in a small part of the domain, and the remaining few nodes are scattered across the other part. Actually, the introduction of an artificial viscosity given by Equation ( 24) ensures a minimum regularity for the node distribution, and as long as there is no interpenetration between particles, the Lagrangian aspect of the formulation may dominate the viscosity effect. In such cases, an adaptation or a node redistribution on the entire domain will become necessary. That is what is known as 'adaptivity' or 'mesh refinement.'

The adaptivity procedure must be carried out with respect to the mass conservation. The total mass of the fluid flowing in the channel, before and after nodal refinement, must be conserved. In our case, this condition is easily respected. To begin the refinement procedure, we need to use an indicator that will provide information regarding whether the concentration of some particles is increasing in a particular area, which could lead to two or more particles sticking together in short order. The indicator we use in this work is based on the distance between the nodes. For every node I , the following entity is computed:

R = min J d I J min J d 0 I J , J is a natural neighbor of I (31)
where d I J is the distance between nodes I and J at time t and d 0 I J is the distance between nodes I and J at time t = 0 s. When R is less than 10%, it means that particle I is too close to one of its neighbors and this kicks in the adaptivity procedure. The two particles must be merged into one particle located at an equal distance between them. The new mass is now equal to the sum of the two masses, and the acceleration and velocity are interpolated linearly.

If the indicator R increases up to 190%, then an insertion of a node is achieved. The new particle will be located at an equal distance between the particle and its furthest neighbor. The mass is distributed between them avoiding the occurrence of any numerical perturbance. To that end, few iterations are completed in order to stabilize the simulation before going on to the next time step.

It should be noted here that for the benchmark tests conducted in this work and which are presented in the following paragraphs, the adaptivity procedure was used only when adding new nodes in corners.

COMPUTATIONAL ALGORITHM

The global algorithm of the method is given as follows:

1. Carry out geometrical initialization and construct the Voronoi diagram. 2. Stock all Voronoi cell lattices. They will be used for water depth computing.

For every time step:

• For every node I = 1, 2, . . . , M the acceleration is computed using Equation ( 26).

• If I is a boundary node, treat it as described in Section 6.

• The velocity is computed in every node using time integration.

• Using Equation [START_REF] Cueto | Overview and recent advances in natural neighbor Galerkin methods[END_REF], compute the accelerations at the initial Voronoi lattices (for mass conservation). • Update the velocities and the coordinates of the initial Voronoi diagram lattices using Equations ( 28) and ( 29). • Update the particle coordinates using Equation [START_REF] Klein | A linear-time randomized algorithm for the bounded Voronoi diagram of a simple polygon[END_REF].

• With the new node positions, re-mesh the entire domain: Delaunay triangulation and Voronoi subdivision. • Compute the new water depths using Equation ( 14).

• Compute the new time step using Equation [START_REF] Martinez | Updated Lagrangian free surface flow simulations with natural neighbour Galerkin methods[END_REF].

It should be noted that Voronoi cell reconstruction in every time step is simply a 'luxury' we accorded ourselves, while we possess a performing and fast Voronoi constructor. However, it does ensure better accuracy and consequently good results.

NUMERICAL EXAMPLES

In the following, we will present some applications that represent a validation step of our new approach. As seen in the preceding paragraphs, the formulation we present is applied to inviscid shallow water flows modelized with St Venant equations ( 12) and ( 13) by neglecting source terms (friction and bed slope terms). The introduction of these terms will be covered in future publications.

Rectangular channel

The first benchmark test is obviously the rectangular section channel. This case is a generalization of the 1D case in which the transversal component does not-or only slightly-influence the behavior of the flow. This flow could be seen as a set of 1D flows set one next to the other transversally. An analytical solution for the 1D case was given by Stoker in [START_REF] Stoker | Water Waves[END_REF]. The initial conditions are presented in Figure 8. The water depth is 1 m upstream and 0.25 m downstream. The channel length is 2000 m and its width is 80 m. The dam is located at x = 1000 m. The depth discontinuity is linearized on a distance of 15 m. The number of particles used for the simulation was between 1000 and 30 000. However, a minimum of 1000 in the longitudinal direction is necessary to well capture the shape of the solution. In Figure 8, 10 000 particles were used as 10×1000 particles.

The shock-capturing term [START_REF] González | Numerical integration in natural neighbour Galerkin methods[END_REF] was used with a free parameter chosen at between 0.5 and 2. In this case = 2 is used. We can note, therefore, that the solution of the 1D dam break case is properly reproduced. Figure 9 shows a tri-dimensional view of the depth at t = 15 s, whereas Figure 10 makes a comparison with the analytic solution. We can see the presence of the diffusive effect of the artificial viscosity especially at the rarefaction wave location. This could be explained by the fact that the Lax-Friedrich scheme does not include any MUSCL-type reconstruction to reduce the diffusive effect commonly associated with first-order Godunov-type schemes [START_REF] Ata | A stabilized SPH method for inviscid shallow water flows[END_REF]. Therefore, the shock capturing is very successful and the stabilization precess works well when there is no numerical oscillation. The use of a free parameter as a flux limiter that depends on the depth gradient (such as = 0 |h i -h j |/h i with 0 = 2) gives a better result from a stabilization point of view, but has more diffusive effect.

Figure 11 shows the evolution of L 2 error in time. We can see that this latter remains acceptable even though the Lagrangian character of the flow implies frequent relative displacement between particles, which could affect the solution accuracy and thus could increase the global error.

Convergent and divergent channels

For now, we will try to explore the effect of the geometry on the solution. We will first examine the influence of a transversal component on the behavior of the solution. To that end, we will present the cases of a convergent channel and of a divergent channel. Figure 12 presents the divergent and convergent channel cases. The water depth is 1 m upstream and 0.25 m downstream. The artificial viscosity parameters are the same as in the last test. To achieve this simulation 10 760 particles were used. Here the domain length was still 2000 m and the dam was located at x = 1000 m, whereas the discontinuity was linearized over a distance of 15 m.

As we can see in these figures, the introduction of the transversal component of the flow due to the change in the geometry in this direction does not affect the performance of the solution. The shocks are well captured and the stabilization tools perform as in the last case. Figure 13 gives us an idea about the velocity distribution on a divergent channel. We can see that the velocity vectors follow the geometry of the domain and the perfect slipping condition is verified.

Cylindrical dam break

In his book [START_REF] Toro | Shock-capturing Methods for Free Surface Flows[END_REF], Toro presents the cylindrical dam break test. It consists of a cylindrical dam with a radius R = 5 m located in an unbounded domain. For the numerical implementation, the domain is taken as measuring 100 m×100 m; the water depth is 1 m in the dam and 0.95 m outside; the number of particles used here is 2490 and it was increased up to about 10 000. Figures 14 and15 illustrate this case. The behavior of the method with respect to the number of particles and the variation of the artificial viscosity parameter is quite similar as illustrated in Ata and Soulaïmani [START_REF] Ata | A stabilized SPH method for inviscid shallow water flows[END_REF] to which we therefore refer the reader for more details. In our case, we will focus on the importance of the 'nature' of the mesh or the initial nodal distribution. Figure 14 shows the case of a structured mesh of 80×80. This figure gives the initial nodal distribution, the initial mesh, the contour of the solution, the velocity distribution and, finally, a 3D view of the water depth. As we can see, the axisymmetry of the problem is lost once the wave starts to move. This problem is overcome in the second case. In fact, Figure 15 illustrates the case of a pseudo-circular nodal distribution. The quality of the solution is strongly improved and, in particular, the axisymmetry of the solution is almost revealed. The caption of shocks is well done and the birth of a second wave in the middle of the domain is shown clearly in this figure. Actually, we should note the presence of numerical oscillations in both cases and especially in Figure 15. This could be explained by the low nodal density in the vicinity of the central point of the dam. Another probable reason is the non-compatibility of the Voronoi subdivision of the domain (which is a polygonal subdivision) with the circular geometry and the circular nature of the flow.

Rectangular channel with a slight contraction

In the second step of the investigation of the sensitivity of the method to the geometry used, a slight contraction is introduced in the case of a rectangular channel dam break. Thus, two angular points are introduced into both sides of the channel. The inclination angle is 2.29 • (10 m/100 m). The initial conditions are presented in Figure 16. The channel length is 200 m; its width is 80 m in its broad part and 60 m in its narrow part. The dam is located at x = 85 m; water depth is 1 m upstream and 0.95 m downstream. The CFL is taken as 0.3. The discontinuity in the water depth was linearized across a distance of 15 m.

Figure 17 gives a 3D view of the water depth obtained at t = 5 s. We can see that the sudden geometry change induces the birth of a wave perturbance, which is expected and is physically explained. In fact, when the shock (front) wave reaches the geometry deviation, it submits to a reflection, which is concretized by a sudden water-height elevation that vanishes as the wave continues to move downstream. We can, consequently, affirm that, the method overcomes problems generated by the presence of a slight geometrical deviation. More complicated geometrical applications will be covered in the future works.

Rectangular channel with a sudden contraction

The ultimate test is a rectangular channel with a sudden contraction. It is a test with a severe change in the geometry. This case will show the capacity of the scheme to handle wide variations in the geometry. The length of the domain is 200 m; its width is 80 m upstream the dam and 20 m downstream. The dam is located exactly at x = 100 m. The water depth is 1 m upstream and 0.95 m downstream. Figure 18 shows the initial conditions for this problem. The water discontinuity is linearized over a distance of 10 m. of 0.5 was fixed for the time integration, a value that induces a time step in the neighborhood of 0.123 s in the beginning, which decreases down to 0.121 s at the end of the simulation. We can conclude that the solution is generally well reproduced and that the waves are well captured. The shape and the smoothness of the solution is slightly perturbed in the vicinity of sharp angles of the geometry. This could be naturally explained by the fact that the linear interpolation discussed in Section 6 cannot be achieved for nodes located in the corner. In this case, the application of the continuity equation ( 12) is used to determine the water depth in this node. Few iterations are needed in order to reach the optimal value of h.

CONCLUSION

For the numerical approximation of the inviscid SWE, a new finite-volume method based on a flux computation using a natural neighbor approximation on the Voronoi diagram is applied. Moreover, the application of this method is achieved within a purely Lagrangian collocation scheme. A perfect slipping condition is assumed at the boundary. The stabilization of the obtained scheme is achieved using an upwinding artificial viscosity that is obtained from an analogy with Riemann solver techniques. The first benchmark tests performed and shown in this paper reveal good results and hold out the promise for more interesting ones in the future. The presence of a complicated geometry and the introduction of source terms (bathymetry and bed friction) will be covered in the following steps before we go on to handle real cases of river or flood simulations.

Figure 1 .

 1 Figure 1. Delaunay triangulation (left) and its relative Voronoi diagram (right) for a 9-node domain.

Figure 2 .

 2 Figure 2. Sibson vs Laplace interpolations [19].

Figure 3 .

 3 Figure 3. Support of the NEM-type interpolation function vs those of MLS and EFG-figure from [19].

Figure 4 .

 4 Figure 4. Sibson and Laplace shape functions: (a) nodal distribution; (b) Sibson shape function; and (c) Laplace shape function-figure from [22].

Figure 5 .

 5 Figure 5. The mass conservation strategy (•: node, : 'mass' cell vertex, →: velocity). Up: domain decomposition at time t n . Down right: real Voronoi cell at time t n+1 used for momentum computation. Down left: mass cell or sub-domain used for mass conservation and water depth computation at time t n+1 .

Figure 7 .

 7 Figure 7. Boundary treatment: a particle on the boundary is fixed and its physical parameters ( u, u, and h) are interpolated from parameters of boundary particles at time t n+1 .

Figure 8 .

 8 Figure 8. Initial conditions for the rectangular channel dam break problem.

Figure 9 .Figure 10 .

 910 Figure 9. Water depth at time t = 15 s.

Figure 11 .

 11 Figure 11. Evolution of L 2 error in time.

Figure 12 .

 12 Figure 12. Divergent channel at time t = 15 s (left) and convergent channel at t = 30 s (right).

Figure 13 .

 13 Figure 13. Velocity distribution for the case of divergent channel (zoom).

Figure 16 .Figure 17 .

 1617 Figure 16. Geometry and contour lines for the rectangular channel with a slight contraction.

Figure 18 .Figure 19 .

 1819 Figure 18. Initial conditions for the problem of rectangular channel with a sudden contraction.

1 HFigure 20 .

 120 Figure 20. Water depth for the rectangular channel with a sudden contraction at time t = 5s.

  diagram is constructed at each time step. This procedure induces a geometrical reconstruction of the cells associated with each particle at every time step. Rigorously, the time derivative of the momentum should be discretized, for instance, as

	Dm I u I /Dt (h n+1 I A n+1 I depend on u n+1 A n+1 I u n+1 I	-h n I A n I u n I )/ t, which defines a non-linear problem (while h n+1 I	and

ACKNOWLEDGEMENTS

The authors wish to acknowledge the contribution of Dr L. A. Illoul for his invaluable help in the realization of this work. This work is partially funded by the CRSNG and Hydro-Quebec http://www.interscience. wiley.com/jpages/0271-2091/.