# Towards an efficient numerical treatment of the transport problems in the resin transfer molding simulation

Fernando Sánchez, Juan Antonio Garcia, Llanos Gascón, Francisco Chinesta

## ► To cite this version:

Fernando Sánchez, Juan Antonio Garcia, Llanos Gascón, Francisco Chinesta. Towards an efficient numerical treatment of the transport problems in the resin transfer molding simulation. Computer Methods in Applied Mechanics and Engineering, 2007, 196 (21-24), pp.2300-2312. 10.1016/j.cma.2006.07.018. hal-01007392

HAL Id: hal-01007392

https://hal.science/hal-01007392

Submitted on 2 Apr 2017

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

Public Domain
# Towards an efficient numerical treatment of the transport problems in the resin transfer molding simulation

F. Sánchez<sup>a</sup>, J.A. García<sup>b</sup>, Ll. Gascón<sup>c,*</sup>, F. Chinesta<sup>d</sup>

Efficient liquid composites molding (LCM) processes simulation requires an efficient treatment of the advection equation governing the evolution of different process variables (volume of fluid, heat transfer, incubation time, etc.). In a previous work [J.A. García, Ll. Gascón, F. Chinesta, A fixed mesh numerical method for modelling the flow in liquid composites molding processes using a volume-of-fluid technique, Comput. Methods Appl. Mech. Engrg. 192 (7-8) (2003) 877-893], a second-order scheme with flux limiters was proposed for the integration of the advection equation which governs the volume fraction evolution. Due to the fact that other properties, like the incubation time, are not defined in the empty part of the mold, some numerical difficulties are found during the process updating [F. Chinesta, T. Mabrouki, A. Ramón, Some difficulties in the flow front treatment in fixed mesh simulations of composites forming processes, in: 5th ESAFORM Conference on Material Forming, Kracow, 2002, pp. 295-298]. Then the scheme described in García et al. (2003) must be modified to circumvent the difficulty just refereed. This paper describes a new flux limiter fixed mesh technique for the calculation of the incubation time and the fluid fraction in mold filling simulation of thin cavities containing a fiber reinforcement. Numerical results in two dimensions are presented in order to validate the proposed numerical strategy.

**Keywords:** Composite forming processes; Liquid injection molding; Transport problems; Fixed mesh simulation; Flux limiters

## 1. Introduction

Resin transfer molding (RTM) is a process for the manufacture of fiber-reinforced composites. In recent years, liquid composites molding processes (LCM) and, particularly, resin transfer molding (RTM), are considered more and more in industrial sectors like the aerospace and automotive industries. An interesting feature of the RTM process is its potential to become a high-speed process for the manufacturing of components at a low cost, giving cost

savings and performance improvements over traditional methods.

In the RTM process dry fiber reinforcement is preshaped and oriented into a skeleton known as the preform, which is inserted into a mold cavity which has the shape of the desired part. The mold is then closed and a pressurized low viscosity reactive fluid is injected into the mold cavity. During this stage, known as the injection phase, the resin impregnates and wets out the fibers. The curing reaction begins when the prepolymers are mixed just before the injection into the mold. During the curing reaction the resin polymerizes to become rigid plastic. The filling process in RTM can be investigated by using numerical simulation and/or experimental analysis. Because the experimental analysis is often expensive and time-consuming, the

<footnote><sup>a</sup> Escuela Superior de Enseñanzas Técnicas, Universidad Cardenal Herrera CEU, C/San Bartolomé 55, E-46115 Alfara del Patriarca, Valencia, Spain</footnote>

<footnote><sup>b</sup> Depto. Ing. Mecánica y Materiales, Universidad Politécnica de Valencia, Camino Vera s/n, 46022 Valencia, Spain</footnote>

<footnote><sup>c</sup> Depto. Matemática Aplicada, Universidad Politécnica de Valencia, Camino Vera s/n, 46022 Valencia, Spain</footnote>

<footnote><sup>d</sup> Laboratoire de Mécanique des Systèmes et des Procédés, UMR CNRS-ENSAM-ESEM, 151 Boulevard de l'Hôptital, F-75013 Paris, France</footnote>
numerical simulation seems to be an appealing and efficient tool to analyze the mold filling stage, and it allows to reduce significantly the number of experimental tests needed.

An important issue in the simulation of the mold filling process is the numerical treatment of the transient free surface (or moving boundary) defined by the liquid resin flow front. The fluid flow described in this work represents a typical case of a free surface flow. Fluid flows are classified as free surface flows when a fluid has one or more boundaries in the form of interfaces with another fluid and the external fluid has a negligible influence on the motion of the internal fluid. This assumption is valid, typically, in situations where the viscosity and density of the internal fluid are much larger than the viscosity and density of the external fluid. Two different types of numerical procedures are being generally used to simulate it: moving mesh techniques and fixed mesh methods. The moving mesh techniques [3] are based on remeshing the resin saturated part of the mold at each time step and they have the drawback related to the dependence of its accuracy on the frequency of this remeshing with the associated CPU cost. Some methods have been developed to reduce the frequency of the remeshing (see the works of Tezduyar [4–7] and its review paper [8] and the references therein). However, in some problems it is not be possible to reduce the remeshing frequency to a practical level. Two alternative Eulerian numerical techniques using a fixed mesh have been developed for accounting moving interfaces (see [9,10] for a recent review). These include the Volume-of-Fluid techniques (since the pioneer work of Hirt and Nichols [11] to other more recent alternatives introduced in [12–14], for instance) and the level set methods [15]. Volume-of-Fluid techniques describe the interface by discretizing the computational domain into a collection of cells, where a variable referring the fluid volume fraction is assigned to each cell. Finally, let us mention that level set methods have also been successfully used for solving free surface flows, however, the mathematical model is different because they describe the interface from the zero level of a distance function.

Some fixed mesh methods in the RTM context [16] use two mold discretizations, one to solve motion equations by using the finite element method, and the other one to calculate the flow front evolution by using a Volume-of-Fluid method. In other techniques the control volumes coincide with the finite elements. In this case, Trochu et al. [17] propose the use of a non-conforming linear finite element approximation of the pressure field, in which, in spite of that the pressure is discontinuous on the element boundaries (except at the middle point), the velocity remains continuous between elements, yielding a numerical model that satisfies locally the resin mass conservation. In general, equations of motions are solved in the part of the mold filled by the resin. The scheme considered in this work to simulate the RTM process is based on an extended formulation to the whole domain (mold). Moreover, it uses the same fixed mesh for solving momentum equations and for updating the volume of fluid.

Another important variable in the RTM process is the elapsed time since the reactives are mixed, also known as *incubation time* [2], which for each fluid particle depends on the path travelled throughout the mold, and is directly related with the resin curing conversion. In the case of complex molds, with multiple injection points, racetracking effects, zones with different permeability, etc., the incubation time at the flow front can be significantly different to the injection time. Injection strategy is a practical issue in RTM in which the resin can be injected into the mold using multiple gates. Traditionally the selection of gate and vent locations in mold design has been based on the designer knowledge and trial and error attempts, but many recent research studies have been conducted to optimize the process by using computer simulation [18,19]. In [20] the incubation time has been employed as a variable related with an optimal performance by means of the definition of a process performance index based on gate-distance and incubation-time for liquid composite molding processes. This index allows the consideration of both resin flow and curing issues in process design optimization. Consequently, an accurate numerical simulation able to compute the *age* of the fluid located at each point of the mold is very useful for simulating the curing reaction.

In general, efficient 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 reactive, etc. Our goal in the present work is the development of an accurate, efficient and robust numerical method for the numerical integration of the advection equation governing the evolution of a generic variable *J* associated with the fluid.

Because of the fluid incompressibility, if a function satisfies an advection equation, then it can be expressed as a conservation law using the incompressibility of the fluid, and therefore one can apply the extensive knowledge for high-resolution numerical methods for conservation laws to devise a method for integrating the advection equation (see, for instance [21–23]). These finite difference methods are designed for accounting accurately the presence of discontinuities or high gradients in the problem solution by dropping order (upwinding) in the vicinity of a discontinuity, while in smooth regions the high-order of accuracy remains unchanged. Furthermore, these methods have a rigorous mathematical formulation that admits precise numerical analysis in special cases. An important kind of these schemes are the flux limiters, introduced by Sweby [24], and based on the definition of hybrid schemes that use second-order fluxes in smooth regions of the solution, with a limitation of the numerical fluxes to first-order, to avoid numerical oscillations in the vicinity of the discontinuities.

In a former work [1], the authors described a 2D numerical model for RTM simulations in which the volume fraction was updated by using a new second-order scheme with flux limiters, that provide excellent results regardless of the
position of the injection points. That Volume-of-Fluid scheme was based on a simple heuristic 2D extension of the one proposed in the 1D case in [24] to construct finite difference schemes for modelling solutions of scalar conservation laws. However, the updating of certain variables, such as the incubation time, introduces numerical difficulties related to the fact that they are not defined in the empty part of the mold. Thus, it is necessary to develop a scheme able to be applied to the updating of any variable related to the fluid in the framework of fixed mesh simulations.

In this work, we propose a general strategy able to carry out an accurate field transport when an Eulerian technique, using a fixed mesh, is adopted. We redefine the flux limiters proposed in [1] in order to account for the numerical difficulties related to the integration of variables related to the fluid, such as the incubation time, in fixed mesh simulations, when they cannot be defined properly in the empty part of the mold. Thus, we propose a new second-order technique for the numerical integration of the volume fraction, which is employed to determine the new position of the flow front, and for updating the incubation time in the RTM process, both variables being governed by their corresponding linear advection equations. The proposed technique will be applied to solve some problems related to the resin transfer molding.

The organization of this paper is as follows. In the next section, the governing equations and boundary conditions for the RTM flow model are presented. In order to compute the pressure and velocity fields we use an extended variational formulation that combines the governing equations associated with the filled and empty domains. In Section 3, we present a new strategy to define schemes for integrating the advection equation governing the evolution of some variables related to the fluid. Thus, a new version of the flux limiters technique has been developed in the two-dimensional case in order to update the volume fraction and the incubation time in the context of the RTM processes. Section 4 presents some numerical results. The paper finishes with some concluding remarks in Section 5.

## 2. Mechanical modelling

The fluid flow problem is defined in a volume $\Omega$ (whole domain or mold cavity), being $\partial\Omega$ its boundary, which is divided into two subdomains

$$ \Omega = \Omega_f(t) \cup \Omega_e(t), $$

where $\Omega_f(t)$ denotes the resin saturated part of the mold at time $t$, and $\Omega_e(t)$ is the empty part of the domain at that time. The boundaries of $\Omega_f(t)$ and $\Omega_e(t)$ are $\partial\Omega_f(t)$ and $\partial\Omega_e(t)$, respectively. The filled and empty mold regions have a common boundary, which corresponds to the fluid flow front $\partial\Omega_{ad}(t)$ (see Fig. 1). Therefore, the presence at time $t$ of the fluid in $\Omega$ is defined by the fluid fraction function $I(x, t)$ which takes values between 0 and 1 and represents the saturation level at each point.

<figure><img src="image_1.png" /><figcaption>Fig. 1. Two-dimensional model.</figcaption></figure>

The resin impregnation is usually modelled as a flow in porous media. The model equations which describe this process are given by

* Darcy's law, which states that the flow velocity, $v$, is proportional to the pressure gradient, according to

$$ v = -(K/\mu)\nabla p, \quad (1) $$

where $K$ is the preform permeability tensor, $\mu$ is the viscosity of the fluid and $p$ is the pressure. For the sake of simplicity, from now on, we assume constant permeability and viscosity.

* The incompressibility of the fluid, expressed as

$$ \nabla \cdot v = 0. \quad (2) $$

The pressure distribution is found in the RTM simulations by replacing Darcy's law (1) into Eq. (2). Then, the flow velocity derives from the pressure gradient, according to the Darcy's law (1).

* The evolution of the resin volume fraction $I$ (also known as fluid presence function), which is governed by the following advection equation:

$$ \frac{\partial I}{\partial t} + v \cdot \nabla I = 0. \quad (3) $$

The value of the variable $I$ measures the degree of saturation and will be used to know the position of the flow front at each time step. The volume fraction $I$, is defined as

$$ I(x, t) = \begin{cases} 1 & \text{if } x \in \Omega_f(t), \\ 0 & \text{if } x \notin \Omega_f(t), \end{cases} $$

where $\Omega_f(t)$ denotes the fluid domain at time $t$.

* The evolution of the incubation time $E$, governed by the following transport equation:

$$ \frac{\partial E}{\partial t} + v \cdot \nabla E = 1. \quad (4) $$

The value of the incubation time $E$ is zero in the injection nozzle and varies throughout the filled part of the mold, but it is not defined in the empty mold region, i.e.

$$ E(x, t) = \begin{cases} E(x, t) & \text{if } x \in \Omega_f(t), \\ \text{Not defined} & \text{if } x \notin \Omega_f(t). \end{cases} $$
The initial and boundary conditions are

* The pressure gradient in the normal direction to the mold walls is zero. Physically this means that material cannot leaves the mold cavity through the mold walls.

* The pressure or the flow rate are prescribed on the inflow boundary (injection nozzle), $\partial\Omega^-$, i.e.

$$p(x \in \partial\Omega^-) = P_i \quad \text{or} \quad \underline{v}(x \in \partial\Omega^-) = \underline{v}_i,$$

$$\partial\Omega^- = \{\underline{x} \in \partial\Omega, \underline{v}(\underline{x}) \cdot \underline{n}(\underline{x}) < 0\},$$

where $\underline{n}(\underline{x})$ is the unit outwards vector, defined on the boundary at point $\underline{x}$.

* The pressure is zero on the flow front, that is

$$p(\underline{x} \in \partial\Omega_{\text{ad}}(t)) = 0.$$

* We assume that at time $t=0$, the mold is empty, i.e.

$$I(\underline{x}, t=0) = \begin{cases} 0 & \text{if } \underline{x} \in \Omega, \\ 1 & \text{if } \underline{x} \in \partial\Omega^- \end{cases} \qquad (5)$$

and that for any time

$$I(\underline{x} \in \partial\Omega^-, t) = 1. \qquad (6)$$

* We assume that the resin reactive are mixed on the inflow boundary. In this manner, for the incubation time it results

$$E(\underline{x}, t=0) = \begin{cases} \text{Not defined} & \text{if } \underline{x} \in \Omega, \\ 0 & \text{if } \underline{x} \in \partial\Omega^- \end{cases} \qquad (7)$$

being

$$E(\underline{x} \in \partial\Omega^-, t) = 0. \qquad (8)$$

The simplest incremental explicit numerical algorithm for performing the filling simulation solves, at each time step, until the complete filling of the mold (i.e. while $|\Omega_f(t)| < |\Omega|$, the next four steps:

Step 1. Compute the pressure field using a finite element discretisation of the variational formulation extended to the whole domain $\Omega$, which is equivalent to prescribe a null pressure at the nodes not connected with, at least, a fulfilled element.

Step 2. Once the pressure distribution is known at the nodes of the mesh, the velocities of the elements are obtained applying Darcy's law (Eq. (1)).

Step 3. Update the domain occupied by the fluid by solving the equation governing the fluid fraction evolution, Eq. (3).

Step 4. Update the incubation time $E$, solving Eq. (4).

An efficient technique for the resolution of the transport equations, defined by Eqs. (3) and (4), will be explained in detail in the next section.

### 3. Numerical treatment of the advection equations governing different transports

In the just described algorithm, once the pressure and velocity fields have been calculated at time $t$, the fluid fraction

can be updated in order to redefine the fill and empty domains and to locate ("capture") the new position of the flow front (Step 3). From the flow kinematics and the just computed fluid fraction, the incubation time can be updated (Step 4). The evolution of both, the volume fraction $J$ and the incubation time $E$, is given by the general linear advection equation

$$\frac{\partial J}{\partial t} + \underline{v} \cdot \nabla J = S, \qquad (9)$$

where for the incubation time $J=E$, $S=1$, and for the volume fraction evolution $J=I$, $S=0$. The initial conditions result from Eqs. (5) and (7), being the boundary conditions the ones given by Eqs. (6) and (8), respectively.

However, despite of this apparent similitude, we are going to show that there are important differences in the numerical treatment of both fields due to the fact that the incubation time is not defined in the empty domain, where the fluid does not exist.

When the domain $\Omega$ is discretized into $N$ nonoverlapping cells, i.e.,

$$\Omega = \bigcup_{e=1}^{N} \Omega_{e} \qquad (10)$$

an appropriate technique for the discretization of the advection equation (9) that governs the evolution of a generic field $J$ lies in integrating Eq. (9) in each cell $\Omega_e$

$$\int_{\Omega_e} \left( \frac{\partial J}{\partial t} + \nabla \cdot (J \underline{v}) - J(\nabla \cdot \underline{v}) \right) d\Omega = \int_{\Omega_e} S d\Omega. \qquad (11)$$

From Eq. (11), using the divergence theorem and taking into account the fluid incompressibility, results

$$\int_{\Omega_e} \frac{\partial J}{\partial t} d\Omega + \int_{\partial\Omega_{e+}} J(\underline{v} \cdot \underline{n}) ds + \int_{\partial\Omega_{e-}} J(\underline{v} \cdot \underline{n}) ds = \int_{\Omega_e} S d\Omega, \qquad (12)$$

where $\partial\Omega_{e+}$ denotes the outflow cell boundary and $\partial\Omega_{e-}$ the inflow cell boundary related to $\Omega_e$, and $\underline{n}$ is the unit outward normal defined on the cell boundary $\partial\Omega_e$.

Stabilized discretizations of Eq. (12) usually involve continuous field approximations into each cell, being discontinuous across the cell boundaries. In order to discretize Eq. (12) all integrals must be approximated. One of the main difficulties related to Eq. (12) is that the approximation of $J$ (that as just indicated is not continuous across the cell boundaries) is not defined in the element boundaries along which are defined the boundary integrals. If we consider a piecewise constant approximation of $J$, one could assume that on the cell outflow boundary the function $J$ takes the existing value inside the element $\Omega_e$, i.e. $J(\underline{x} \in \partial\Omega_{e+}) = J_e$ and that on the inflow cell boundary the function $J$ takes the value existing in the neighbor upstream element $\Omega_k$, i.e., $J(\underline{x} \in \partial\Omega_k^k) = J_k(\partial\Omega_k^k = \partial\Omega_{e-} \cap \partial\Omega_k)$. Thus, Eq. (12) can be rewritten in the equivalent form
$$ \frac{\partial J_e}{\partial t} |\Omega_e| = -J_e \int_{\partial \Omega_{e+}} (\underline{v} \cdot \underline{n}) ds - \sum_{k=1}^{k=N_{\text{up}}} J_k \int_{\partial \Omega_e^k} (\underline{v} \cdot \underline{n}) ds + S_e |\Omega_e|, \quad (13) $$

where $|\Omega_e|$ denotes the volume of $\Omega_e$ and $N_{\text{up}}$ the number of upstream neighbor elements, with $\bigcup_{k=1}^{k=N_{\text{up}}} \partial \Omega_e^k = \partial \Omega_{e-}$. For the sake of simplicity in the next paragraphs we assume (without loss of generality) that $N_{\text{up}} = 1$.

Considering a first-order explicit approximation of the time derivative, Eq. (13) can be written as

$$ J_e^{n+1} = J_e^n - J_e^n \frac{\Omega_{e+}}{|\Omega_e|} + J_e^n \frac{\Omega_{e-}}{|\Omega_e|} + (\Delta t) S_e, \quad (14) $$

where now $J_{e-}$ defines the value of $J$ in the single upstream neighbor cell, and we have also introduced the inflow and outflow fluid volumes defined, respectively, by

$$ \Omega_{e-} = \left| \int_{\partial \Omega_{e-}} (\underline{v}_e \cdot \underline{n}) ds \right| \Delta t, \quad (15) $$

$$ \Omega_{e+} = \left| \int_{\partial \Omega_{e+}} (\underline{v}_e \cdot \underline{n}) ds \right| \Delta t. \quad (16) $$

In Eq. (14) the subscript denotes the considered element and the superscript indicates the time step.

The discrete form of the equation governing the fluid fraction evolution is obtained from Eq. (14) taking $J = I$ and $S = 0$

$$ I_e^{n+1} = I_e^n - I_e^n \frac{\Omega_{e+}}{|\Omega_e|} + I_e^n \frac{\Omega_{e-}}{|\Omega_e|}. \quad (17) $$

The incubation time results from Eq. (14) taking $J = E$ and $S = 1$

$$ E_e^{n+1} = E_e^n - E_e^n \frac{\Omega_{e+}}{|\Omega_e|} + E_e^n \frac{\Omega_{e-}}{|\Omega_e|} + \Delta t. \quad (18) $$

We notice again that the value of the incubation time $E$ is set to zero at the injection nozzle, it takes positive values in the filled part of the mold, but it is not defined on the empty region. Thus, a first problem appears if we consider the time $t$ for which an element $\Omega_e$ starts its filling process from its upstream element. At that time, we can assume that the outflow fluid volume vanishes, i.e., $\Omega_{e+} = 0$, which leads to

$$ E_e^{n+1} = E_e^n + E_{e-}^n \frac{\Omega_{e-}}{|\Omega_e|} + \Delta t, \quad (19) $$

Eq. (19) establishes that solution at $t_{n+1}$ is dominated by solution existing in the element at the previous time step (for $\Delta t$ small enough $\Omega_{e-} \approx 0$). A possible solution, widely used, consists in convecting the field existing in the upstream element, and then apply Eq. (19). However, this solution implies the necessity of identify accurately the time at which each element starts its filling stage. We can notice that if we assume wrongly that the incubation time is zero in the empty cells, then Eq. (19) gives

$$ E_e^{n+1} = E_{e-}^n \frac{\Omega_{e-}}{|\Omega_e|} + \Delta t \quad (20) $$

which states that for $\Omega_{e-}$ and $\Delta t$ small enough the fluid located at the flow front is younger than the upstream fluid, which is absolutely wrong.

To overcome this difficulty, we apply the strategy proposed in [2], which lies in multiplying Eq. (4) by the fluid fraction $I$, Eq. (3) by the field $E$ and then combining both resulting equations to obtain

$$ \frac{\partial(IE)}{\partial t} + \underline{v} \cdot \nabla(IE) = I \quad (21) $$

with the initial condition

$$ (IE)(\underline{x}, t) = \begin{cases} E(\underline{x}, t) & \text{if } \underline{x} \in \Omega_f(t), \\ 0 & \text{if } \underline{x} \notin \Omega_f(t) \end{cases} \quad (22) $$

whose discretised form, according to Eq. (14), (with $J = IE$ and $S = I$) is given by

$$ I_e^{n+1} E_e^{n+1} = I_e^n E_e^n - I_e^n E_e^n \frac{\Omega_{e+}}{|\Omega_e|} + I_e^n E_e^n \frac{\Omega_{e-}}{|\Omega_e|} + (\Delta t) S_e, \quad (23) $$

where taking $S_e = I_e^{n+1}$, yields

$$ E_e^{n+1} = \frac{I_e^n E_e^n}{I_e^{n+1}} - \frac{I_e^n E_e^n}{I_e^{n+1}} \frac{\Omega_{e+}}{|\Omega_e|} + \frac{I_e^n E_e^n}{I_e^{n+1}} \frac{\Omega_{e-}}{|\Omega_e|} + \Delta t. \quad (24) $$

It is important to notice that one can compute $E_e^{n+1}$ after the fluid fraction updating which gives $I_e^{n+1}$. Taking into account Eq. (17), Eq. (24) becomes

$$ E_e^{n+1} = \frac{I_e^n E_e^n \left(1 - \frac{\Omega_{e+}}{|\Omega_e|}\right) + I_e^n E_e^n \frac{\Omega_{e-}}{|\Omega_e|}}{I_e^n \left(1 - \frac{\Omega_{e+}}{|\Omega_e|}\right) + I_e^n \frac{\Omega_{e-}}{|\Omega_e|}} + \Delta t. \quad (25) $$

For applying Eq. (25) with high accuracy (avoiding numerical instabilities related to small values of $I_e^n$) we assume that $\Omega_{e-} \neq 0$ if and only if $I_{e-} \approx 1$, and that $\Omega_{e+} \neq 0$ if and only if $I_e \approx 1$. This condition can be introduced in Eq. (25) from the definition of the $\delta$ parameter

$$ \delta_e = \begin{cases} 1 & \text{if } I_e \ge I_{\text{th}} \text{ and } \delta_{e-} = \begin{cases} 1 & \text{if } I_{e-} \ge I_{\text{th}}, \\ 0 & \text{if } I_{e-} < I_{\text{th}} \end{cases} \\ 0 & \text{if } I_e < I_{\text{th}} \end{cases} $$

with $I_{\text{th}}$ a threshold parameter very close to one ($I_{\text{th}} \approx 1$). Thus Eq. (24) can be finally written as

$$ E_e^{n+1} = \frac{I_e^n E_e^n}{I_e^{n+1}} - \delta_e^n \frac{I_e^n E_e^n}{I_e^{n+1}} \frac{\Omega_{e+}}{|\Omega_e|} + \delta_e^n \frac{I_e^n E_e^n}{I_e^{n+1}} \frac{\Omega_{e-}}{|\Omega_e|} + \Delta t \quad (26) $$

being

$$ I_e^{n+1} = I_e^n \left( 1 - \delta_e^n \frac{\Omega_{e+}}{|\Omega_e|} \right) + I_e^n \delta_e^n \frac{\Omega_{e-}}{|\Omega_e|}. $$

To illustrate the accuracy obtained by using Eq. (26) we will consider different particular situations:

(1) The beginning in the element filling. In this case, $I_{e-}^n = 1$, $I_e^n = 0$, $\delta_e^n = 0$, from which Eq. (26) gives $E_e^{n+1} = E_{e-}^n + \Delta t$
i.e., the incubation time takes the upstream element value increased with the time step. This result is in good agreement with the physical considerations.

(2) During the element filling $I_e^n < 1$ and $I_{e-}^n = 1$. In consequence, the fluid transfer from element $\Omega_e$ towards the downstream element vanishes and the updating results

$$E_e^{n+1} = \frac{I_e^n E_e^n}{I_e^{n+1}} + \frac{I_{e-}^n E_{e-}^n}{I_e^{n+1}} \frac{\Omega_{e-}}{|\Omega_e|} + \Delta t.$$

Here the incubation time is updated taking into account the age of the resin in $\Omega_e$, the age of the incoming resin and the time step.

(3) In the fulfilled elements, where the fluid presence function results $I_e^{n+1} = I_e^n \approx 1$ and the updating results

$$E_e^{n+1} = E_e^n - E_e^n \frac{\Omega_{e+}}{|\Omega_e|} + E_{e-}^n \frac{\Omega_{e-}}{|\Omega_e|} + \Delta t$$

that corresponds to the usual discontinuous finite element method for the incubation time, given by Eq. (18).

For the general case, the new basic scheme for the numerical discretization of the advection equation (14) can be rewritten as

$$J_e^{n+1} = J_e^n - \delta_e^n J_e^n \frac{\Omega_{e+}}{|\Omega_e|} + \delta_e^n J_{e-}^n \frac{\Omega_{e-}}{|\Omega_e|} + (\Delta t) S_e^{n+1}. \quad (27)$$

### 3.1. A flux limiter strategy

The numerical scheme described by (27) is only a first-order method. First-order schemes are regarded as too diffusive for many applications, they have a strong numerical dissipation in the neighborhood of discontinuities giving a low accuracy in smooth regions of the solution, which means in RTM a diffusive definition of the flow front. High-order techniques avoiding spurious oscillations, requires the use of particular nonlinear advection schemes. There are many examples of nonlinear nonoscillatory advection schemes in the literature, including schemes based on flux limiters or slope limiters (see [21,22] for a review).

The integration of the advection equation to update the fluid fraction in RTM was accomplished in the 2D case by the authors in [1] introducing an alternative second-order scheme combined with appropriate flux limiters. Since the resolution of the incubation time makes use of the fluid fraction, the scheme described in [1] must be modified for applying it for the incubation time updating.

For this purpose we will introduce in the following paragraphs a new flux limiter strategy for solving problems defined by the generic advection equation, given by Eq. (9) where $J$ is a field that can be related to the fluid, being not defined in the empty domain.

The considered meshes consist of a collection of non-overlapping triangles, that cover the domain $\Omega$. For the generic triangle $\Omega_e$, we denote by $A(e)$ its area, $\Omega_j$ represents a neighbor triangle, $\underline{n}_{ej}$ is the outward unit vector on the common edge of the triangles $e$ and $j$, which is assumed pointing towards $j$, $l_{ej}$ is the length of that edge and the velocity vector in that element is $\underline{v}_e$. In this case, Eq. (12) is modified as follows:

$$\frac{\partial J_e}{\partial t} A(e) + \sum_{j=1}^{3} \int_{\partial \Omega_{ej}} J(\underline{v} \cdot \underline{n}) ds = \int_{\Omega_e} S d\Omega,$$

where the field $J$ is assumed constant into each element, and $\partial\Omega_{ej}$ represents the $j$th face of the triangle $e$. Thus, in meshes composed of triangles, Eq. (14) can be written as

$$J_e^{n+1} = J_e^n - \frac{\Delta t}{A(e)} \underline{v}_e^n \cdot \sum_{j=1}^{3} l_{ej} \hat{f}_{ej}^n \underline{n}_{ej} + (\Delta t) S_e, \quad (28)$$

where $\hat{f}_{ej}^n$ is defined by

$$\hat{f}_{ej}^n = \frac{1}{2} \{(J_j^n + J_e^n) - \text{sign}(\underline{v}_e^n \cdot \underline{n}_{ej})(J_j^n - J_e^n)\} \quad (29)$$

which takes the existing value inside the element $\Omega_e$, i.e., $\hat{f}_{ej}^n = J_e^n$ on the outflow boundaries and the value of the neighbor element, i.e., $\hat{f}_{ej}^n = J_j^n$, on the inflow boundaries.

In order to obtain an accurate transport of the field $J$ we assume that the fluid comes to an element only from its fulfilled neighbor elements, as it was discussed in the previous section. This condition is introduced in the above formulation from the use of the parameter

$$\delta_k = \begin{cases} 1 & \text{if } I_k \ge I_{th} \\ 0 & \text{if } I_k < I_{th} \end{cases} \quad \text{for } k = e, j \quad (30)$$

which is introduced in Eq. (29) in a natural way

$$\hat{f}_{ej}^{\text{L}} = \frac{1}{2} \{(\delta_j^n J_j^n + \delta_e^n J_e^n) - \text{sign}(\underline{v}_e^n \cdot \underline{n}_{ej})(\delta_j^n J_j^n - \delta_e^n J_e^n)\}. \quad (31)$$

Then, if the source term in Eq. (28) is evaluated at the time step $(n+1)$, Eqs. (28) and (31) are equivalent to Eq. (14), for the mesh here considered.

As starting point to extend the flux limiter technique established in [1] we first propose a purely second-order scheme for the numerical integration of Eq. (9). For this end, we suggest that the estimation of $\hat{f}_{ej}$ in Eq. (28) could be expressed as

$$\hat{f}_{ej}^{\text{H}} = \frac{1}{2} \left\{ (\delta_j^n J_j^n + \delta_e^n J_e^n) - \left( \frac{\Delta t}{d_{ej}} \underline{v}_{ej}^n \cdot \underline{n}_{ej} \right) (\delta_j^n J_j^n - \delta_e^n J_e^n) \right\}, \quad (32)$$

where $\underline{v}_{ej}$ denotes an average velocity between the cells $e$ and $j$ (the arithmetic average, for example) and $d_{ej}$ represents the distance between the baricenters of the triangles $e$ and $j$. This approach is based on a heuristic 2D extension of the classical 1D Lax-Wendroff scheme and it includes the fact that neighbor elements can exchange fluid only when the upstream element is fulfilled.
In order to prevent the presence of spurious oscillations in the vicinity of high gradients of discontinuities taken place on the flow front, induced by the second-order scheme defined by Eqs. (28) and (32) we need to introduce a suitable generalization of the flux limiter technique introduced in [1] which takes into account the physics related to the introduction of the $\delta$ parameter given by Eq. (30).

Following the methodology employed in [1], the numerical flux $\hat{f}_{ej}$ in Eq. (28) can be constructed in a nonlinear way by combining high and lower order approximations

$$ \hat{f}_{ej} = \hat{f}_{ej}^{L} + \chi(r_{ej}) (\hat{f}_{ej}^{H} - \hat{f}_{ej}^{L}) $$

or, in an equivalent form

$$ \hat{f}_{ej} = \hat{f}_{ej}^{L} + \frac{1}{2}\chi(r_{ej}) \left( \text{sign}(\underline{v}_{e}^{n} \cdot \underline{n}_{ej}) - \frac{\Delta t}{d_{ej}} \underline{v}_{ej}^{n} \cdot \underline{n}_{ej} \right) \times (\delta_{j}^{n} J_{j}^{n} - \delta_{e}^{n} J_{e}^{n}). \quad (33) $$

Here $\chi$ is an adaptive correction called flux limiter. This function allows to add a proper amount of nonlinear antidiffusion to the lower-order flux approximation $\hat{f}_{ej}^{L}$ improving its accuracy without generating spurious oscillations in presence of high gradients. For this reason, the flux limiter should be equal to zero in the vicinity of high gradients, and very close to one in regions where the solution is smooth enough.

To determine the admissible antidiffusion correction, the flux limiter $\chi$ should depend on the local smoothness of the solution. The smoothness quantification is carried out by introducing a suitable smoothness indicator which uses the upstream and downstream solution slopes. For classical 1D flux limiters this quantification is made from the ratio between the upstream and the downstream slopes at node $i$. In the 2D case we propose the following quantification:

$$ r_{ej} = \begin{cases} \min_{k \text{ adjacent to } e \text{ and } \cos(\theta_{ek}) < 0} & \left\{ \frac{\delta_{e}J_{e} - \delta_{k}J_{k}}{\delta_{j}J_{j} - \delta_{e}J_{e}} \right\} & \text{if } \cos(\theta_{ej}) \ge 0, \\ \min_{k \text{ adjacent to } j \text{ and } \cos(\theta_{jk}) < 0} & \left\{ \frac{\delta_{k}J_{k} - \delta_{j}J_{j}}{\delta_{j}J_{j} - \delta_{e}J_{e}} \right\} & \text{if } \cos(\theta_{ej}) < 0 \end{cases} $$

in which $\theta_{ej}$ denotes the angle between $\underline{n}_{ej}$ and the velocity vector of element $e$.

The construction of the above smoothness indicator is based on the following considerations:

* In order to detect discontinuities in the numerical solution, when the flow goes from $e$ to $j$ we compare the gradients of the solution between the triangles $e$ and $j$ with the gradients of the solution between the triangle $e$ and their adjacent triangles from which the flow is coming to the triangle $e$.

* When the flow is coming from $j$ to $e$, we compare the gradients of the solution between the triangles $e$ and $j$ with the gradients of the solution between the triangle $j$ and their adjacent triangles from which the flow is coming to the triangle $j$.

Therefore, the estimation of $r_{ej}$ includes the parameter $\delta$, accounting the physics previously described. However, we observe that this condition can be ignored assuming

$$ \delta = 1 \quad \text{if } I \le 1. \quad (34) $$

In this case, the scheme allows that neighbor elements can exchange fluid even when they are not fulfilled. Eq. (34) is equivalent to Eq. (30) when we consider $I_{th} = 0$.

A variety of flux limiters have been proposed in the literature (we just mention the general reviews given in [21–23]). In the present work the Superbee limiter has been considered, which is the most compressive classical limiter

$$ \chi_{SB}(r) = \max\{0, \min(2r, 1), \min(r, 2)\}. $$

## 4. Numerical examples

In this section the accuracy and order of convergence of the proposed algorithms are analyzed using different numerical test.

### 4.1. Saturated flow

Firstly we analyze the accuracy of the different numerical schemes just described for the discretization of the advection equation, considering the saturated mold depicted in Fig. 2. In that case a saturated porous media using unidirectional flow with constant velocity has been considered (the preform permeability $K$ and the resin viscosity $\mu$ are $10^{-7}$ m$^2$ and 0.1 Pa s, respectively). In this case we avoid the resolution of the Darcy’s equation, focusing in the errors induced by the discretization of the advection equations.

Now, we are solving the equation governing the evolution of $I$, but due to the fact that the considered domain is fulfilled, $I$ is not more, in these examples, the fluid fraction, but a generic property of the fluid whose evolution is governed by Eq. (3).

Exact and numerical solution concerning the fluid property $I$ are depicted in Figs. 3–5 for a mesh of 1600 elements using, respectively, the first-order scheme (Eqs. (28) and (31)), the second-order one (Eqs. (28) and (32)), and finally using the Superbee flux limiter scheme just proposed (Eqs. (28) and (33)). In these three schemes, the parameter $\delta$ has

<figure><img src="image_2.png" /><figcaption>Fig. 2. Saturated flow conditions.</figcaption></figure>
<figure><img src="image_3.png" /><figcaption>Fig. 3. Exact solution of *I* at *t* = 3.6 s (left) and numerical prediction using the first-order scheme described by Eqs. (28) and (31) without including the *δ* parameter (right).</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Fig. 4. Exact solution of *I* at *t* = 3.6 s (left) and numerical prediction using the second-order scheme described by Eqs. (28) and (32) without including the *δ* parameter (right).</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 5. Exact solution of *I* at *t* = 3.6 s (left) and numerical prediction using the Superbee flux limiter scheme described by Eqs. (28) and (33) without including the *δ* parameter (right).</figcaption></figure>

been omitted (i.e. expression (34) is considered) allowing a
fluid exchange between partially filled elements.

It can be noticed that the use of the first-order or the
second-order schemes produces significant deviations with
respect the exact solution. In the second-order case
(Fig. 4) due to the numerical oscillations in the vicinity of
the flow front. The first-order technique (Fig. 3) avoids
oscillations, but it introduced an over-diffusion in the front
neighborhood. The best results has been obtained using the
Superbee flux limiter (Fig. 5), described by Eqs. (28) and
(33) even when the δ parameter is not taken into account.

In order to obtain a better representation of numerical
results Fig. 6 shows the profile of *I* in the flow direction
obtained using a mesh with 240 elements. Fig. 6 compares
the solutions obtained with the first-order and the Superbee
limiter schemes both computed without using the δ para-
meter, with the exact solution. A similar comparison of
results has been omitted for the purely second-order case
(Fig. 4) because of the oscillatory character of the numeri-
cal solution. In Fig. 6 we can notice the expected numerical
dissipation characteristic of the first-order schemes as well
as the accuracy obtained by using the Superbee scheme
whose solution describes perfectly the front position with-

out introducing the spurious oscillations characteristic of
purely high-order schemes.

Although the numerical results obtained for a generic
field advected with an equation like Eq. (3) without consid-
ering the δ parameter are accurate enough when an appro-
priate flux limiter technique is applied, this scheme is not
suitable for the resolution of the advection equation gov-
erning the evolution of a generic fluid property, when it
cannot be defined in the empty part of the mold, as it is
the case of the incubation time. Fig. 7 compares the exact
solution of the incubation time for a filling time of 4.35 s
with the numerical one obtained with the Superbee flux
limiter scheme described by Eqs. (28) and (33), without
considering the δ parameter. As can be noticed in Fig. 7
(right) the value of the incubation time is too high for
the elements that start its filling. The introduction of the
δ parameter seems to be then pertinent for the integration
of the incubation time.

Fig. 8 shows the numerical results for the incubation
time obtained using a mesh with 240 elements and both
the first-order scheme and the Superbee flux limiter tech-
nique considering the δ parameter defined by Eq. (30).
These numerical results are compared with the exact
<figure><img src="image_6.png" /><figcaption>Fig. 6. Numerical results of *I* obtained with the first-order scheme (left) and with the Superbee flux limiter (right) without taking into account the *δ* parameter. Both solutions are compared with the exact one.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Fig. 7. Exact solution of the incubation time for a filling time of 4.35 s (left) and numerical solution obtained with the Superbee flux limiter scheme described by Eqs. (28) and (33) without introducing the *δ* parameter (right).</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Fig. 8. Numerical results for the incubation time obtained with the first-order scheme (left) and with the Superbee flux limiter technique (right) using the *δ* parameter defined in Eq. (30), both compared with the exact solution.</figcaption></figure>

As noticed, the use of the adapted Superbee flux limiter allows to obtain a better accuracy in the definition of the incubation time not only in the front neighborhood, but also in the region characterized by smooth variations of the solution.

Fig. 9 shows the exact solution of *I* compared with the one that results from the first-order scheme assuming the value of *δ* given by Eq. (30). Fig. 10 shows similar simulation that Fig. 6 when the parameter *δ* is introduced in both formulations (first-order and Superbee). In this case, we remark that the consideration of the parameter *δ* in the Superbee scheme for the calculation of *I* yields the inabilitation of the flux limiter. In fact, we can observe that when the flow goes from triangle *e* to *j*, $\delta_e = I_e = 1$, $\delta_k = I_k = 1$ and $\delta_j = 0$ and then $r_{ej} = 0$; and when the flow is coming to the triangle *e*, $\delta_j = I_j = 1$, $\delta_k = I_k = 1$ and $\delta_e = 0$ and

then $r_{ej} = 0$ again. Then, from Eq. (33), $\hat{f}_{ej} = \hat{f}_{ej}^L$ and the limited second-order term vanishes. This conclusion can be appreciated in Fig. 10 which shows that the first-order solution (at the left) and the Superbee one (at the right) are equivalent for the calculation of the volume fraction $I$. However, we can notice the accuracy in the front description.

The convergence analysis for this problem has been carried out when the front reaches the outflow boundary. The errors quantified by the L2-norm for different mesh sizes have been plotted in Figs. 11 and 12. We can notice that the order of convergence is always higher when the Superbee flux limiter is used instead of the first-order upwind scheme, independently of the consideration of the *δ* parameter. However, when the parameter *δ* is included in the schemes from Eq. (30) (Fig. 11), the order of convergence
<figure><img src="image_9.png" /><figcaption>Fig. 9. Exact solution of $I$ at $t = 3.6$ s (left) and numerical prediction using the first-order scheme described by Eqs. (28) and (31) including the $\delta$ parameter defined by Eq. (30) (right).</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Fig. 10. Numerical results of $I$ obtained with the first-order scheme (left) and with the Superbee flux limiter (right) using the $\delta$ parameter. Both solutions are compared with the exact one.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Fig. 11. Evolution of the error in the incubation time solution with the mesh size for the numerical schemes defined by Eqs. (28) and (31) and by Eqs. (28) and (33) using the $\delta$ parameter defined by (30).</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Fig. 12. Evolution of the error in the incubation time solution with the mesh size for the numerical schemes defined by Eqs. (28) and (31) and by Eqs. (28) and (33) without considering the $\delta$ parameter.</figcaption></figure>

for the Superbee flux limiter technique is greater than the one obtained in absence of this parameter (Fig. 12).

Analogously, using as reference the incubation time solution, Fig. 13 illustrates the L2-errors for different mesh sizes obtained by the schemes described by Eqs. (28)–(31), and by Eqs. (28)–(33) using the $\delta$ parameter defined by (30) for the above problem when a nonsaturated flow is considered, coupling the Darcy’s equation resolution with the ones governing the evolution of $I$ and $E$. The results indicate similar tendencies with respect to the order of convergence, which is also higher when the Superbee flux limiter is used instead of the first-order upwind scheme, although in this case both schemes give lower accuracy and order of

<figure><img src="image_13.png" /><figcaption>Fig. 13. Evolution of the error in the incubation time solution with the mesh size for the numerical schemes defined by Eqs. (28) and (31) and by Eqs. (28) and (33) using the $\delta$ parameter defined by (30), for a similar flow problem involving non-saturated flow.</figcaption></figure>
convergence than the saturated flow problem because of
the errors introduced in the resolution of the kinematics.

4.2. *Mold with different permeabilities*

In order to illustrate the applicability of the scheme
defined by Eqs. (28) and (33) using the δ parameter defined
by (30) in more complex geometries, a particular filling case
of the mold shown in Fig. 14 has been considered. In this
case, a finite element mesh of 800 elements of a mold
1000 × 500 × 5 mm with lateral injection and a constant

flow rate of 5 cm³/s has been considered. The preform per-
meability was 10⁻⁸ m² everywhere except in two areas, one
with 10⁻⁷ m² and the other one with 10⁻⁹ m². Due to the
differences in the preform permeabilities the flow is signifi-
cantly perturbed. The resin viscosity was 0.1 Pa s in all
simulations.

Figs. 15–20 show the numerical resolution of the trans-
port problem. Fig. 15 shows the mold flow pattern, that is
the time at which each point is reached by the resin flow.
Figs. 16–18 illustrate the effects of the different permeabil-
ities on the flow kinematics and pressure distribution. It

<figure><img src="image_14.png" /><figcaption>Fig. 14. Finite element mesh, gate locations and zones of different permeability.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Fig. 15. Mold filling pattern (in s) for a filling time of 47 s.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Fig. 17. Velocity (m/s) averaged in nodes for a filling time of 37 s.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>Fig. 16. Pressure distribution (Pa) and node averaged velocities for a filling time of 37 s.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>Fig. 18. Pressure distribution (Pa) and node averaged velocities for a filling time of 47 s.</figcaption></figure>
<figure><img src="image_19.png" /><figcaption>Fig. 19. Incubation time (s) for a filling time of 37 s.</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>Fig. 20. Incubation time (s) for a filling time of 47 s.</figcaption></figure>

can be noticed that *young resin* moves faster through the high permeability area, mixing with the *old resin* located in the low permeability area (see Figs. 19 and 20).

## 5. Conclusions

This paper presents a new approach for solving different transport problems arising from the RTM simulation. Some of these problems involve fields related to the fluid, and in consequence cannot be properly defined in the empty domain, such as the incubation time. In this form, the application of standard discretization techniques can fail.

An efficient numerical technique which allows to solve the incubation time and the volume fraction transport for RTM process has been achieved. The incubation time calculation accuracy in a fixed mesh discretization depends directly on the previous resolution of the element volume fraction, which allows to locate the flow front position of the injected resin. The accurate resolution of both advection problems has been implemented by using similar strategies based on the original definition and implementation of the flux limiter techniques [24].

The development of numerical techniques to updating some fluid properties that introduce a particular difficulty related to the fact that it is not defined in the empty part of the mold is based on two basic ideas. On the one hand, we replace the correspondent transport equation by another one that combines those fields that are not defined in the empty part of the mold with the fluid presence function in order to integrate a new advection equation associated with

a field defined in the whole domain. On the other hand, in order to avoid numerical instabilities related to small values of the volume fraction at the beginning of an element filling we assume that the fluid comes to an element only from its fulfilled neighbor elements. This condition has been introduced in the formulation of the schemes from the use of a parameter $\delta$ that depends on the fluid fractions of the elements that exchange fluid and which has to be considered in all stages of the discretization.

With this aim in view, some different numerical techniques have been proposed, tested and compared in order to obtain an efficient treatment of the fluid properties transport through the whole domain in fixed mesh simulations. The aforementioned schemes have been implemented and validated in order to solve the flow front evolution and the incubation time updating through the whole domain using a fixed mesh in RTM processes. The schemes formulation has been generalized to any fluid property and allows an efficient treatment of its transport. A new version in 2D of the Superbee flux limiter, which includes the parameter $\delta$ in its formulation, has been found to be the most robust. The numerical results validate its use in the simulation of the advection in RTM process.

For a complete process modelling it will be necessary to introduce the energy equations as well as the resin viscosity evolution during the curing. This modelling can be easily achieved in the framework of the numerical strategy proposed in this paper. Since the code implemented treats only with bidimensional cases it can not deal with the energy equation resolution for which is necessary a 3D numerical resolution. Hence, important process variables as temperature and curing may be necessary to simulate as well. We think that the technique proposed in this work can be used to simulate efficiently these advection problems.

Our results about the resin flow front and incubation time profile throughout the mold can be an important simulation tool with low computation cost because can be applied to bidimensional isothermal simulations and allow an optimized RTM process design. Some related goals are: to decline incomplete curing, avoid non homogeneous curing, the use of genetic algorithms in optimization of the nozzle and vents location (which computational implementation cost is very high), implementation of online curing control strategies, etc.

## Acknowledgments

This research work is supported by a grant from Ministerio de Ciencia y Tecnología (Spain) and another one from Fondo Europeo de Desarrollo Regional (FEDER) through Project DPI2004-03152 and by Generalitat Valenciana under Project GVA05/209.

## References

[1] J.A. García, Ll. Gascón, F. Chinesta, A fixed mesh numerical method for modelling the flow in liquid composites molding processes using a
volume-of-fluid technique, Comput. Methods Appl. Mech. Engrg. 192 (7-8) (2003) 877-893.

[2] F. Chinesta, T. Mabrouki, A. Ramón, Some difficulties in the flow front treatment in fixed mesh simulations of composites forming processes, in: 5th ESAFORM Conference on Material Forming, Kracow, 2002, pp. 295-298.

[3] S. Li, R. Gauvin, Numerical analysis of the resin flow in resin transfer molding, J. Reinf. Plast. Compos. 10 (3) (1990) 314-327.

[4] T.E. Tezduyar, M. Behr, S. Mittal, A.A. Johnson, Computation of unsteady incompressible flows with the stabilized finite element methods: space-time formulations, iterative strategies and massively parallel implementations, New Methods Transient Analysis, ASME AMD-Vol. 143 (1992) 7-24.

[5] A.A. Johnson, T.E. Tezduyar, Mesh update strategies in parallel finite element computations of flow problems with moving boundaries and interfaces, Comput. Methods Appl. Mech. Engrg. 119 (1994) 73-94.

[6] T.E. Tezduyar, Computation of moving boundaries and interfaces and stabilization parameters, Int. J. Numer. Methods Fluids 43 (2003) 555-575.

[7] K. Stein, T.E. Tezduyar, R. Benney, Automatic mesh update with the solid-extension mesh moving technique, Comput. Methods Appl. Mech. Engrg. 193 (2004) 2019-2032.

[8] T.E. Tezduyar, Finite element methods for flow problems with moving boundaries and interfaces, Arch. Comput. Methods Engrg. 8 (2001) 83-130.

[9] R. Scardovelli, S. Zaleski, Direct numerical simulation of free-surface and interfacial flow, Annu. Rev. Fluid Mech. 31 (1999) 567-603.

[10] J.A. Sethian, P. Smereka, Level set methods for fluid interfaces, Annu. Rev. Fluid Mech. 35 (2003) 341-372.

[11] C.W. Hirt, B.D. Nichols, Volume-of-fluid (VOF) method for the dynamics of free boundaries, J. Comput. Phys. 39 (1981) 201-225.

[12] J. López, J. Hernández, P. Gómez, F. Faura, A volume of fluid based on multidimensional advection and spline interface reconstruction, J. Comput. Phys. 195 (2004) 718-742.

[13] M.A. Cruchaga, D.J. Celentano, T.E. Tezduyar, Moving-interface computations with the edge-tracked interface locator technique (ETILT), Int. J. Numer. Methods Fluids 47 (2005) 451-469.

[14] J.E. Pilliot, E.G. Puckett, Second-order accurate volume-of-fluid algorithms for tracking material interfaces, J. Comput. Phys. 199 (2004) 465-502.

[15] S. Osher, J.A. Sethian, Fronts propagating with curvature dependent speed: algorithms based on Hamilton-Jacobi formulations, J. Comput. Phys. 79 (1988) 12-49.

[16] S.G. Advani, Flow and Rheology in Polymer Composites Manufacturing, Elsevier, 1994.

[17] F. Trochu, R. Gauvin, D.M. Gao, Numerical analysis of the resin transfer molding process by the finite element method, Adv. Polym. Technol. 12 (4) (1993) 329-342.

[18] S. Jiang, C. Zhang, B. Wang, A process performance index and its applications to optimization of the RTM process, Polym. Compos. 22 (5) (2001) 690-701.

[19] J. Luo, Z. Liang, C. Zhang, B. Wang, Optimum tooling design for resin transfer molding with virtual manufacturing and artificial intelligence, Composites Part A 32 (2001) 877-888.

[20] F. Sánchez, J.A. García, F. Chinesta, Ll. Gascón, C. Zhang, Z. Liang, B. Wang, A process performance index based on gate-distance and incubation time for the optimization of gate locations in liquid composite molding processes, Composites Part A 37 (2006) 903-912.

[21] C. Hirsch, Numerical computation of internal and external flows, Computational Methods for Inviscid and Viscous Flows, vol. 2, Wiley, Chichester, 1990.

[22] R. Leveque, Numerical Methods for Conservation Laws, Birkhäuser, Basel/Boston/Berlin, 1990.

[23] E.F. Toro, Riemann Solvers and Numerical Methods for Fluid Dynamics, Springer-Verlag, New York, 1997.

[24] P.K. Sweby, High resolution schemes using flux limiters for hyperbolic conservation laws, SIAM J. Numer. Anal. (1984) 995-1011.