A flux limiter strategy for solving the saturation equation in RTM process simulation

SUMMARY : The main aim of this work is the proposal of a numerical procedure for solving the saturation equation in RTM process simulation. In order to analyze in more detail the progressive impregnation of a fibrous preform by a fluid resin, the numerical model here proposed considers the flow through a partially saturated medium, including the dependence of permeability on the saturation degree. The model consists of an elliptic equation governing the pressure distribution and a transport hyperbolic equation governing the evolution of the saturation. A global flux limiter fixed mesh strategy is proposed for solving the transport equation with source term which describes the saturation evolution in RTM process. The flux limiter method has the ability to limit the extra numerical diffusion introduced by standard first-order schemes. This formulation could be an appealing choice for improving RTM flow simulations and optimize injection processes. Some preliminary numerical results are presented in order to validate the proposed numerical strategy.


INTRODUCTION
Liquid Composites Molding (LCM) processes are based on a proper resin impregnation of the reinforcement. Modeling and simulation play an important role in the development and optimization of molds production and in devising appropriate resin injection strategies. Minimization of mold filling time without losing part quality is an important issue in Resin Transfer Molding (RTM). Inadequate injection strategies tend to create macro and microvoids in the part, the formation of which depends on fluid flow velocity.
Over the last decades, different numerical techniques have been used to model the mold filling process in RTM. In general, RTM process simulation requires an accurate treatment of the advection equation which governs the evolution of different fluid properties: fluid presence function, incubation time, temperature, concentration of reacting components, degree of saturation, etc. In some of our former works [1,2] we proposed a new procedure to integrate accurately the transport of different fluid properties. In this work, that numerical procedure is extended for solving the transport equation governing the evolution of the saturation. In this model the permeability is assumed to be a function of saturation [3], and then the continuity equation that governs the pressure distribution, includes a source term which depends on the saturation. The presence of this source term explains the quadratic spatial pressure distribution in an unsaturated porous medium noticed in unidirectional flows.
To derive a closed model the saturation equation is considered whose source term depends on the micro and macro voids content.
The technique here used is based on a flux limiter fixed mesh strategy for solving the transport equation which governs the evolution of degree of saturation. It is well-known that spurious oscillations in the computed solution can be avoided by using first-order upwind discretization schemes, but unfortunately these first order schemes introduce an excessive numerical diffusion. One alternative for avoiding non-physical oscillations, reducing the extra numerical diffusion, lies in the use of second-order numerical fluxes in the regions where the solution is smooth enough that reduce to first-order in presence of discontinuities. This technique was successfully used for computing the volume fraction evolution in [1].

GOVERNING EQUATIONS
The model which describe the mould filling process is given by Darcy´s law (1) combined with the continuity equation, which should account for the loss of resin due to the perform impregnation. This can be done by incorporating a source term in the mass balance equation (2) where v is the velocity field, K is the preform permeability tensor, µ is the fluid viscosity, φ is the porosity, p is the pressure and S the degree of saturation. In order to evaluate the permeability, we make use of the model proposed by Breard et al. [3] for predicting void formation in LCM. Then, the permeability in Darcy's law should be considered as the product of the geometrical permeability, K Sat , and a relative permeability K rel (S), where the relative permeability depends on the saturation degree (4) where β is a fitting factor whose usual values range in the interval [3]. Combining Eqs. (1) and (2) yields (5) Finally, to make a closed description of the flow, we assume the saturation equation proposed in [3,5], which governs the transport of the variable S: (6) where the source term can take different forms. According to the model proposed by Breard et al. [3], the source terms writes: where the model proposed by Ruiz et al. [5] is assumed. In this last model, α M and α m represent the dispersive coefficients of the macro and micro voids, respectively, and Q S is a function of the modified capillary number [3].
The saturation S takes a unit value in the saturated domain, a zero value in the empty region and varies between 0 and 1 in the partially saturated region. The whole domain has been designated by Ω and its boundary by ∂Ω. Initially, we assume the condition (9) where the saturation function on the inflow boundary is assumed unchanged during the filling process, that is: The filling process simulation involves at each time step: 1. The calculation of saturation dependent permeability as well as the calculation of the source term of Eq. (5). 2. The calculation of the pressure distribution by applying a standard finite element discretization of Eq. (5). 3. The calculation of the velocity field from Darcy´s law (1). 4. The updating of the saturation distribution by integrating Eq. (6) using a flux limiter technique.
The boundary conditions are given by: • The pressure gradient in the normal direction to the mold walls vanishes.
• The pressure or the flow rate is specified on the inflow boundary (injection nozzle).
• The pressure is zero in the empty part of mold.

THE ADVECTION EQUATION
The saturation equation has a hyperbolic character, needing for appropriate stabilizations. In this section we propose a flux-limiter technique for its solution. For the sake of simplicity from now on we only consider one dimensional models.
In the one-dimensional case Eq. (6) reads: and Eq. (11) can be integrated by applying a second-order upwind scheme preserving the TVD property [6], whose discrete form writes: where and Here h represents the mesh size, Δt the time step and χ(r) is the flux limiter function. For the aproximation of the source term, and considering the model (7) or using the expression of the source term given by (8) The superscript UP is associated with the first-order upwind scheme and the superscript SW with the second-order scheme using a modified flux limiter technique (in our case the Sweby flux limiter). To define the coefficient r i+1/2 in Eq.(13), we propose to define its value comparing consecutive variations of the approximate numerical solution with respect to the flow direction [6]: The Sweby flux limiter reads [6]: In general, if what guarantees that the scheme will be of first-order in the neighbourhood of a discontinuity, since implies that the slopes of the solution have opposite signs. On the other hand, it is necessary that to obtain a second-order scheme in smooth regions of the solution. The limiter function should be chosen verifying both conditions and maximizing the antidiffusive flux. We remark in Eq.(13) that the choice gives rise to the original first-order upwind method, whereas the choice defines a fully second-order scheme. Hence, the proposed hybrid scheme, with in Eq. (13), is of second-order in smooth regions of the solution but it is a firstorder method when a discontinuity is present. This allows avoid the spurious numerical oscillations associated with the conventional second-order methods in the presence of discontinuities and reduces the excessive numerical diffusion introduced by the first-order upwind schemes.

NUMERICAL SIMULATIONS
In order to analyze the accuracy and efficiency of the proposed numerical scheme just described for the discretization of the transport equation governing the evolution of the saturation, where the source term is given by Eq. (8), we analyze some one dimensional described [3]. A mold of 0.5m length is considered. A constant injection pressure is prescribed (10 -5 Pa) being the saturated permeability K sat and the resin viscosity 5.10 -11 m 2 and 0.1Pa.s, respectively. For the numerical simulation, we consider a time step of 0.1 seconds, α M =α m = 10 -10 , and a constant value R S =0.8. The domain is assumed initially empty, except the first element that represents the injection nozzle that is assumed full-filled.
To analyze the influence of mesh size on the simulated results, three meshes with different nodal distributions are considered, consisting in 30, 60 and 150 nodes, respectively. The associated numerical solutions are depicted in Fig. 1. for a filling time of 700 seconds, using the first-order scheme (Eqs. (12)-(13) with χ(r)=0) and using the Sweby flux limiter scheme (Eq. (12)-(13) with χ(r) defined by Eq.(17)). We can notice that the convergence is faster when the flux limiter second-order scheme is considered.   1. proves also that the diffusivity of the scheme decreases with the mesh refinement. It can be noticed that the use of a first-order scheme introduces a significant numerical overdiffusion. Obviously, there are two terms that contribute to the front smoothing, the one related to the source term, and the one purely numerical introduced by the dicretization scheme. The last one can be reduced by using higher-order schemes and fine enough meshes.  To quantify the convergence different computations were performed using different mesh sizes (30, 60, 90, 120 and 150 nodes, respectively). The computed results are illustrated in Fig. 3. that represents the partially filled domain to the total length of the mold for both discretization schemes.
The pressure distribution varies linearly within each saturated element and is continuous between elements. The saturation and pressure distributions for a filling time of 700 seconds are represented in Fig. 4. It can be noticed the presence of three different behaviors: fully impregnation, unsaturated and empty domains. In the fully saturated regions the model reduces to the usual Darcy model, leading to a linear pressure distribution. In the partially saturated regions, the pressure distribution becomes parabolic, as announced in [3]. Finally, in the empty domain the pressure vanishes.

CONCLUSIONS
A new numerical procedure for simulating LCM processes has been presented. The numerical model is based on the consideration of partially saturated flows. For this purpose, the advection-diffusion equation describing the evolution of the saturation is solved by using a flux limiter upwind scheme. Numerical results confirm that first-order schemes exhibit an excessive and no realistic diffusion due to the numerical approximation of the advective term, while the flux-limiter scheme shows less extra-diffusive effects. Thus, the flux limiter proposed improves significantly the results (with respect to the first-order solutions). The results suggest that the present numerical method has a satisfactory capability of simulating the LCM process but further convergence analysis must be carried out.