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

Emmanuelle Abisset-Chavanne, Francisco Chinesta

► To cite this version:

Emmanuelle Abisset-Chavanne, Francisco Chinesta. Toward an optimisation of the reactive resin transfer molding process: thermo-chemico-mechanical coupled simulations. International Journal of Material Forming, 2013, 7 (2), pp.249 - 258. 10.1007/s12289-013-1124-0 . hal-01007388

HAL Id: hal-01007388

https://hal.science/hal-01007388

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
# Toward an optimisation of the reactive resin transfer molding process: thermo-chemico-mechanical coupled simulations

Emmanuelle Abisset-Chavanne · Francisco Chinesta

**Abstract** 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.

**Keywords** RTM process · Multiphysic coupled model · Numerical simulations

## 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 [1–3] and the chemical reactions and couplings in the others (e.g. in [4]). 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.

E. Abisset-Chavanne (✉) · F. Chinesta
GeM - Institut de recherche en Génie Civil et Mécanique,
UMR CNRS-ECN-Université de Nantes,
1 rue de la Noe, BP 92101, 44321 Nantes Cedex 3, France
e-mail: emmanuelle.abisset-chavanne@ec-nantes.fr
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 volume-of-fluid scheme [3]. 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 physics 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 ($\eta$) with the polymerization degree ($\alpha$) of the resin.

The first point has been solved by using the strategy proposed in [5] and whose main ideas are revisited and

generalized in what follows. The second point has been treated following [4] 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 $\Omega$ which can be divided into two subdomains

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

$\Omega_f(t)$ representing the part of the mould saturated with the resin and $\Omega_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 $\Omega_e$ and $\Omega_f$.

The main difficulty when simulating this problem by using standard fixed mesh methods is the dependency of the front position with the mesh [6, 7]. Moreover, even if advanced techniques are used to deal with the front updating [8–11], 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 [5], authors propose the use of a volume-of-fluid method (derived from the one defined in [12]) coupled with the introduction of a fluid fraction function $I(\mathbf{x}, t)$ which takes values between 0 and 1 and represents the resin saturation level at each point:

$$ I(\mathbf{x}, t) = \begin{cases} 1 & \text{if } \mathbf{x} \in \Omega_f(t) \\ 0 & \text{if } \mathbf{x} \in \Omega_e(t) \end{cases} \qquad (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 $\Omega$ 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 $\Omega_f$ by the Darcy’s law combined with the fluid incompressibility, and is equal to zero in $\Omega_e$:

$$ \begin{cases} \nabla(-\mathbf{K}/\mu \cdot \nabla p) = 0 & \forall \mathbf{x} \in \Omega_f(t) \\ p(\mathbf{x}) = 0 & \forall \mathbf{x} \in \Omega_e(t) \end{cases} \qquad (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 Ω:

$$ \int_{\Omega} [\nabla p^* (\mathbf{K}/\mu)^{-1} \nabla p \cdot I + p^* p \cdot (1-I)] d\Omega = 0 \quad (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 [4], the following polymerization kinetic has been identified (Fig. 1):

$$ \frac{d\alpha}{dt} = K \cdot \exp\left(\frac{-E_{\text{act}}}{RT}\right) \cdot \alpha^m \cdot (1-\alpha)^n \quad (4) $$

<figure><img src="image_1.png" /><figcaption>Fig. 1 Evolution of the degree of polymerization with time calculated with Eq. 4 (see [4])</figcaption></figure>

with

$$ K = 4,574 \cdot 10^6 \quad (5) $$

This polymerization process affects the temperature evolution by introducing a source term in the energy balance (heat equation):

$$ \begin{gathered} \rho c_p(T) \frac{DT}{dt} = k(T) \Delta T + \rho w \Delta H \frac{\partial \alpha}{\partial t} \\ \text{with } \frac{D\bullet}{dt} = \frac{\partial\bullet}{\partial t} + \mathbf{v} \cdot \nabla\bullet \end{gathered} \quad (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:

$$ \begin{aligned} c_p(T) &= a_{c_p} \cdot T + b_{c_p} \\ k(T) &= a_k \cdot T + b_k \end{aligned} \quad (7) $$

Concerning the relation between the viscosity and the degree of polymerization, the following law has been identified (also in [4], Fig. 2):

$$ \eta(T, \alpha) = 2,14 \cdot \exp(-0,028 \cdot T + (12,08 + 0,016T) \cdot \alpha) \quad (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 the 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

<figure><img src="image_2.png" /><figcaption>Fig. 2 Evolution of the viscosity with the polymerization degree calculated with Eq. 8</figcaption></figure>
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:

$$ \mathbf{v} = -(\mathbf{K}/\eta) \cdot \nabla p \quad (9) $$

* the incompressibility of the fluid:

$$ \nabla \mathbf{v} = 0 \quad (10) $$

* the evolution of the resin fraction function I which is governed by the advection equation:

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

* the transport and the evolution of the polymerization degree which is governed by the advection equation:

$$ \frac{\partial \alpha}{\partial t} + \mathbf{v} \cdot \nabla \alpha = K \cdot \exp\left(\frac{-E_{\text{act}}}{RT}\right) \cdot \alpha^m \cdot (1-\alpha)^n \quad (12) $$

* the transport and the evolution of the temperature which is governed by the advection equation:

$$ \rho c_p(T) \frac{\partial T}{\partial t} + \mathbf{v} \cdot \nabla T = k(T) \Delta T + \rho w \Delta H \frac{\partial \alpha}{\partial t} \quad (13) $$

* the evolution of the viscosity which is governed by the chemo-rheological law defined previously:

$$ \eta(T, \alpha) = 2, 14 \cdot \exp(-0,028 \cdot T + (12,08 + 0,016T) \cdot \alpha) \quad (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) $\partial\Omega^-$;

- the flow is saturated at the nozzle::

$$ \forall t, I(\mathbf{x} \in \partial\Omega^-, t) = 1 \quad (15) $$

- at time $t=0$, the mould is empty:

$$ \forall \mathbf{x} \in \Omega, I(\mathbf{x}, t=0) = 0 \quad (16) $$

* for the temperature:

- the mould, and all the domain boundaries (excluding the injection nozzle) is at the imposed temperature $T_{\text{mould}}$;

- the fluid is injected at the temperature $T_{\text{inj}}$;

* for the polymerization degree:

- no polymerization occurs before the resin enter the mould:

$$ \forall t, \alpha(\mathbf{x} \in \partial\Omega^-, t) = 0 \quad (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 $\Omega$ [5]:

$$ \int_{\Omega} [\nabla p^* (\mathbf{K}/\mu)^{-1} \nabla p \cdot I + p^* p \cdot (1-I)] d\Omega = 0 \quad (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 [5]. 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:

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

The initial and boundary conditions are given respectively by Eqs. 15 and 16.

We now discretize the whole domain $\Omega$ 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:

$$ \int_{\Omega_e} \frac{\partial I}{\partial t} d\Omega + \int_{\partial \Omega_{e+}} I(\mathbf{v} \cdot \mathbf{n}) ds + \int_{\partial \Omega_{e-}} I(\mathbf{v} \cdot \mathbf{n}) ds = 0 \quad (20) $$

where $\partial \Omega_{e-}$ denotes the inflow cell boundary and $\partial \Omega_{e+}$ the outflow one, $\mathbf{n}$ is the unit outward vector defined on the cell boundary $\partial \Omega_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:

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

where the subscript refers to the considered element and the superscript to the time step, $|\Omega_e|$ denotes the surface of the element $e$, $I_{e-}$ defines the value of $I$ in the upstream neighbour, $\Omega_{e+}$ and $\Omega_{e-}$ are the inflow and outflow fluid volumes defined by:

$$ \left\{ \begin{array}{l} \Omega_{e+} = \left| \int_{\partial \Omega_{e+}} (\mathbf{v}_e \cdot \mathbf{n}) ds \right| \Delta t \\ \\ \Omega_{e-} = \left| \int_{\partial \Omega_{e-}} (\mathbf{v}_{e-} \cdot \mathbf{n}) ds \right| \Delta t \end{array} \right. \quad (22) $$

However, this scheme can introduce too diffuse flow fronts and certain instabilities related to small values of $I_e^n$. To overcome these effects, we introduce the parameters $\delta_e$ and $\delta_{e-}$ defined as:

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

that leads to:

$$ I_e^{n+1} = I_e^n - \delta_e I_e^n \frac{\Omega_{e+}}{|\Omega_e|} + \delta_{e-} I_{e-}^n \frac{\Omega_{e-}}{|\Omega_e|} \quad (24) $$

From a physical point of view these parameters ensure that an element must be almost fully filled ($I_e \simeq 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:

$$ \frac{\partial J}{\partial t} + \mathbf{v} \cdot \nabla J = S \quad (25) $$

where $J$ stands for the transported quantities ($T$ or $\alpha$) 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:

$$ J_e^{n+1} = J_e^n - \delta_e J_e^n \frac{\Omega_{e+}}{|\Omega_e|} + \delta_{e-} J_{e-}^n \frac{\Omega_{e-}}{|\Omega_e|} + S^{n+1} \Delta t \quad (26) $$

However, it has been demonstrated in [5] 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_e^{n+1} = \delta_{e-} J_{e-}^n \frac{\Omega_{e-}}{|\Omega_e|} + S^{n+1} \Delta t \quad (27) $$

This equation shows that, for $\Omega_{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 [3, 5] consists in considering the quantity $I \cdot J$ instead of $J$, that allows rewrite Eq. 26 as:

$$ J_e^{n+1} = \frac{I_e^n J_e^n}{I_e^{n+1}} - \delta_e \frac{I_e^n J_e^n}{I_e^{n+1}} \frac{\Omega_{e+}}{|\Omega_e|} + \delta_{e-} \frac{I_{e-}^n J_{e-}^n}{I_{e-}^{n+1}} \frac{\Omega_{e-}}{|\Omega_e|} + \frac{S^{n+1}}{I_e^{n+1}} \cdot \Delta t \quad (28) $$

Now, if we consider again the element that just begins its filling process at time $t_{n+1}$ and that consequently has $J_e^n = 0$, then by combining the Eqs. 24 and 28, we obtain:

$$ J_e^{n+1} = J_{e-}^n + S^{n+1} \Delta t \quad (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_e^{n+1} = J_{e-}^n \quad (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 $|\Omega|/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 noticed 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 [4] 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.

<table><caption>**Table 1** Parameters values used to performed the simulations</caption><thead><tr><th colspan="5">Chemical parameters (with respect to Eq. 4)</th></tr><tr><th>E<sub>act</sub> (kJ/mol)</th><th>m</th><th>n</th><th>R</th><th>ΔH (J/g)</th></tr></thead><tbody><tr><td>80.2</td><td>0.3</td><td>1.05</td><td>8.314</td><td>-53</td></tr></tbody><thead><tr><th colspan="6">Resin properties (with respect to Eq. 8)</th></tr><tr><th>ρ (kg/m<sup>3</sup>)</th><th>w</th><th>a<sub>k</sub></th><th>b<sub>k</sub></th><th>a<sub>cp</sub></th><th>b<sub>cp</sub></th></tr></thead><tbody><tr><td>1560</td><td>0,4</td><td>0,5836</td><td>-0,0012</td><td>1205,3</td><td>1,6</td></tr></tbody><thead><tr><th colspan="5">Process parameters</th></tr><tr><th>Preform permeability ()</th><th>V<sub>inj</sub> (m.s<sup>-1</sup>)</th><th>T<sub>inj</sub> (K)</th><th>T<sub>mould</sub> (K)</th><th></th></tr></thead><tbody><tr><td>10<sup>-7</sup></td><td>0,05</td><td>493</td><td>493</td><td></td></tr></tbody></table>
<figure><img src="image_3.png" /><figcaption>**Fig. 3** Fluid fraction function and polymerization state of the resin at $t = 8$ s for a simulation without chemo-rheological coupling</figcaption></figure>

Each quantity transported in our problem is governed by an advection equation like:

$$ \frac{\partial T}{\partial t} + \mathbf{v} \cdot \nabla T = S \qquad (31) $$

By using the just described strategy, this equation is discretized and integrated in each element as follows:

$$ J_e^{n+1} = \frac{I_e^n J_e^n}{I_e^{n+1}} - \delta_e \frac{I_e^n J_e^n}{I_e^{n+1}} \frac{\Omega_{e+}}{|\Omega_e|} + \delta_{e-} \frac{I_{e-}^n J_{e-}^n}{I_e^{n+1}} \frac{\Omega_{e-}}{|\Omega_e|} + \frac{S^{n+1}}{I_e^{n+1}} \Delta t \qquad (32) $$

In the following, we will focus on pure convection so we consider that $S^{n+1} = 0$ and then:

$$ J_e^{n+1} = \frac{I_e^n J_e^n}{I_e^{n+1}} - \delta_e \frac{I_e^n J_e^n}{I_e^{n+1}} \frac{\Omega_{e+}}{|\Omega_e|} + \delta_{e-} \frac{I_{e-}^n J_{e-}^n}{I_e^{n+1}} \frac{\Omega_{e-}}{|\Omega_e|} \qquad (33) $$

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 $\Delta I$ depending on the flow rate and the time step.

Because we introduce the parameters $\delta_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:

$$ \begin{cases} I_e^n < \delta_e \\ I_e^{n+1} = I_e^n + \Delta I_e^n > 1 \end{cases} \qquad (34) $$

at 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_e^{n+1}$. 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_e^n < \delta_e$, the fluid in element $e$ cannot be transferred to the neighbour elements located downstream.

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 non-negligible 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 $\delta_e$. Then, whatever the elementary fluid fraction $\Delta 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 $\delta_e$ should be proposed.

<figure><img src="image_4.png" /><figcaption>**Fig. 4** Fluid fraction function and polymerization state of the resin at $t = 8$ s for a simulation with chemo-rheological coupling</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>**Fig. 5** Void formation due to a too high injection velocity</figcaption></figure>

The excess of fluid is directly linked to the flow rate and to the time step. The flow rate is dictated by the injection flow rate $V_{inj}$ imposed by the process, and then we cannot use it to avoid or reduce the extra-filling.

Thus, the unique parameter that we can use to reduce it is the time step $\Delta t$.

The first idea is to choose a time step in order to ensure:

$$ \forall n, \forall e, I_e^{n+1} = I_e^n + \Delta I_e^n \le 1 $$

By this way, we ensure no excess. However, as pointed out in the last section, changing the time step depending on 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 process-ability window that ensure a correct filling of the mould.

<figure><img src="image_6.png" /><figcaption>**Fig. 6** SVD method applied to the separation time/space and $T_{inj}$ of the calculated polymerization degree $\alpha$</figcaption></figure>
Such work has already been performed for forming process optimization in [13]. The idea is to use the PGD method [14] 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(\mathbf{x}, t, T_{\text{mould}}) = \sum F_i(\mathbf{x}, t) \cdot G_i(T_{\text{mould}}) \quad (35)$$

The result of the SVD applied on the evolution of the polymerization degree is presented on Fig. 6. It can be noticed 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(\mathbf{x}, t, V_{\text{inj}}) = \sum F_i(\mathbf{x}, t) \cdot G_i(V_{\text{inj}}) \quad (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-chemical mechanical 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.

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

## References

1. Trochu F, Gauvin R, Gao DM (1993) Numerical analysis of the resin transfer molding process by the finite element method. Adv Polym Technol 12:329–342

2. Advani SG (1994) Flow and rheology in polymer composites manufacturing. Elsevier

3. Chinesta F, Mabrouki T, Ramón A (2002) Some difficulties in the flow front treatment in fixed mesh simulations of composites forming processes. In: 5th ESAFORM conference on material forming, Kracow, pp 295–298

4. Luisier A, Bourban P-E, Manson J-AE (2001) Time-temperature-transformation diagram for reactive processing of polyamide 12. J Appl Polym Sci 81:963–972

5. Sanchez F, Garcia JA, Gascón Ll, Chinesta F (2007) Towards an efficient numerical treatment of the transport problems in the resin transfer molding simulation. Comput Methods Appl Mech Eng 196:2300–2312
6. Scardovelli R, Zaleski S (1999) Direct numerical simulation of free-surface and interfacial flow. Annu Rev Fluid Mech 31: 567-603

7. Sethian JA, Smereka P (2003) Level set methods for fluid interfaces. Annu Rev Fluid Mech 35:341-372

8. Hirt CW, Nichols BD (1981) Volume-of-fluid (VOF) method for the dynamics of free boundaries. J Comput Phys 39:201-225

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

10. Cruchaga MA, Celentano DJ, Tezduyar TE (2005) Moving-interface computations with the edge-tracked interface locator technique (ETILT). Int J Numer Methods Fluids 47:451-469

11. Pilliot JE, Puckett EG (2004) Second-order accurate volume-of-fluid algorithms for tracking material interfaces. J Comput Phys 199:465-502

12. Sweby PK (1984) High resolution schemes using flux limiters for hyperbolic conservation laws. SIAM J Numer Anal 21:995-1011

13. Ghnatios C, Chinesta F, Cueto E, Leygue A, Poitou A, Breitkopf P, Villon P (2011) Methodological approach to efficient modeling and optimization of thermal processes taking place in a die: application to pultrusion. Compos A: Appl Sci Manuf 42: 1169-1178

14. Chinesta F, Ladevze P, Cueto E (2011) A short review on model order reduction based on proper generalized decomposition. Arch Comput Method Eng 18:395-404