# An Efficient Solver of the Saturation Equation in Liquid Composite Molding Processes

J. A. Garcia, L. I. Gascon, Francisco Chinesta, Edu Ruiz, François Trochu

► To cite this version:

J. A. Garcia, L. I. Gascon, Francisco Chinesta, Edu Ruiz, François Trochu. An Efficient Solver of the Saturation Equation in Liquid Composite Molding Processes. International Journal of Material Forming, 2010, 3 (2), pp.1295-1302. 10.1007/s12289-010-0681-8. hal-01004985

HAL Id: hal-01004985

https://hal.science/hal-01004985

Submitted on 11 Feb 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.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# An efficient solver of the saturation equation in liquid composite molding processes

J. A. García • Ll. Gascón • F. Chinesta • E. Ruiz • F. Trochu

J. A. García · L. Gascón
Universidad Politécnica de Valencia,
Camino Vera s/n,
46022 Valencia, Spain
e-mail: jugarcia@mcm.upv.es

F. Chinesta
EADS Corporate Foundation International Chair,
GEM—Centrale Nantes,
1 rue de la Noe, BP 92101, 44321 Nantes cedex 3, France

E. Ruiz · F. Trochu
Centre de recherche en plasturgie et composites (CREPEC),
Département de Génie Mécanique,
École Polytechnique de Montréal,
Montréal, Canada H3C 3A7

**Abstract** A major issue in Liquid Composite Molding Process (LCM) concerns the reduction of voids formed during the resin filling process. Reducing the void content increases the quality of the composite and improves its mechanical properties. Most of modeling efforts on process simulation of mold filling has been focused on the single phase Darcy's law, with resin as the only phase, ignoring the formation and transport of voids. The resin flow in a partially saturated region can be characterized as two phase flow through a porous medium. The mathematical formulation of saturation in LCM takes into account the interaction between resin and air as it occurs in a two phase flow. This model leads to the introduction of relative permeabilities as a function of saturation. The modified saturation equation is obtained as a result, which is a nonlinear advection-diffusion equation with viscous and capillary phenomena. In this work, a flux limiter technique has been used to solve a modified saturation equation for the LCM process. The implemented algorithm allows a numerical optimization of the injected flow rate which minimizes the micro/macroscopic void formation during

mold filling. Some preliminary numerical results are presented here in order to validate the proposed mathematical model and the numerical scheme. This formulation opens up new opportunities to improve LCM flow simulations and optimize injection molds.

**Keywords** Resin transfer molding · Transport problems · Saturation · Second order schemes · TVD · Flux limiter

## Introduction

Liquid Composites Molding (LCM) processes are based on the impregnation of the fibrous reinforcement by a liquid thermoset resin. Modeling and simulation play an important role in the development and optimization of production molds and in devising appropriate resin injection strategies. Minimization of mold filling time while improving quality of the part is an important issue. Inadequate injection strategies tend to create macro and microscopic voids in the part, the formation of which depends on the local resin flow velocity.

For single phase flow in a homogeneous porous medium, Darcy's law is an expression of momentum conservation at the macroscopic scale. When two or more fluid phases are present, the permeability in Darcy's original equation is replaced by an effective value to accommodate the presence of other phases. This parameter, called relative permeability is expressed as a function of the saturation level in the fibrous reinforcement. On the other hand, if capillary effects are neglected, the saturation equation reduces to a purely advection transport equation but, in general the saturation equation is an advection-diffusion equation which includes the capillary pressure effect. The equations that describe the RTM filling process with void formation are based on a two
phase flow model and lead to a coupled system of a nonlinear advection-diffusion equation for saturation and an elliptic equation for pressure and velocity. The hyperbolic nature of the saturation equation and its strong coupling through relative permeability represent a challenging numerical issue. In this paper, a quite simple but very accurate numerical approach is proposed to solve this complex flow behaviour.

In previous works [1, 2] different experiments were carried out to investigate the process of void formation. These experimental works shown that macrovoids tend to form during injection at low flow rates, due to capillary dominant effects, whereas high injection rates lead to microvoids formation. Therefore, an optimal impregnation velocity exists that minimizes the formation of such macro and microscopic voids. This is well explained by the two types of forces that induce the motion of the fluid trough the dual-scale porous media: the viscous and the capillary forces. These forces are included in the diffusivity coefficient in the saturation equation. Based on these observations, a simplified mathematical model has been proposed in order to model the RTM filling process with void formation. In this model, the diffusivity coefficient in the saturation equation has been replaced by a term which depends on the velocity.

This paper proposes a numerical solution to optimize the filling stage in RTM. This scheme enables process engineers to determine the optimum injection flow rate so that void formation is minimized, thus resulting in an improved quality of the part. Essential to this optimum process design is the numerical simulation of the modified saturation equation. Many numerical methods to solve this type of equations suffer from serious nonphysical oscillations, excessive numerical dispersion or a combination of both. The technique here used is based on a flux limiter fixed mesh strategy for solving the transport equation which governs the evolution of the degree of saturation of porous media [3]. A technique based on the flux limiter strategy was also used to simulate the volume fraction evolution in [4].

## Governing equations

The LCM process has been conventionally treated as a single-phase flow, with resin as the only phase, and then the formation and transport of voids cannot be modelled using this approach. The resin flow in a partially saturated region can be characterized as two phase flow through a porous medium. In this case, the permeability experienced by one phase depends on the degree of saturation of the reinforcement. Thus, the relative permeability depends on the saturation and describes how one phase flows in the presence of the second one. So, the saturation equation and the relative permeability are based on the two phase flow description. In

order to explain the equations proposed in this paper to simulate mold filling with void formation, a description of the two phase flow through a porous medium is needed.

In this case, the two phases will be referred as resin and air and designated by the subscript R and A respectively. Darcy's law can then be written for each phase j, as follows

$$v_j = -\lambda_j(S)\nabla p_j \quad (1)$$

where $v_j$ is the phase velocity, S is the degree of saturation of the reinforcement by the liquid resin, $p_j$ is the phase pressure, $\lambda_j(S) = K_{rj}(S)K/\mu_j$ is the phase mobility, with $K_{rj}(S)$ the relative permeability to phase j, $\mu_j$ the viscosity of phase j and K the permeability tensor, which is here taken to be diagonal.

The equations that describe mass conservation for the resin and air phases, are respectively given by

$$\nabla \cdot v_R = -\phi \frac{\partial S}{\partial t} \quad \text{and} \quad \nabla \cdot v_A = -\phi \frac{\partial (1-S)}{\partial t} \quad (2)$$

where $\phi$ denotes the porosity. Then, Eq. 2 leads to following assumption

$$\nabla \cdot (v_R + v_A) = 0 \quad \text{being } v = v_R + v_A \quad (3)$$

Combining Eqs. 1 and 2, the resulting equation for the saturation in its most general form gives

$$\phi \frac{\partial S}{\partial t} + \nabla \cdot (v_F(S)) = -\nabla \cdot (D_{cF}(S) \nabla S) \quad (4)$$

where

$$D_{cF}(S) = F(S)\lambda_A(S) \frac{\partial P_c}{\partial S} \quad (5)$$

$$\text{with } F(S) = \frac{\lambda_R(S)}{\lambda_R(S) + \lambda_A(S)}$$

is the nonlinear diffusivity coefficient due to capillary pressure $P_c$, defined as $P_c = P_A - P_R$.

For the particular case of $F(S)=S$, the saturation equation gives

$$\phi \frac{\partial S}{\partial t} + \nabla \cdot (v_S) = -\nabla \cdot (D_c(S) \nabla S) \quad (6)$$

where

$$D_c(S) = (1-S)\lambda_R(S) \frac{\partial P_c}{\partial S} = S\lambda_A(S) \frac{\partial P_c}{\partial S} \quad (7)$$

Replacing the total velocity $v$ by $v_R+v_A$ and simplifying, Eq. 7 can be rewritten as follows:

$$\phi \frac{\partial S}{\partial t} + \nabla \cdot (v_R S) = \nabla \cdot \left( (1-S) \lambda_R(S) \frac{\partial P_R}{\partial S} \nabla S \right) \quad (8)$$

It can be observed that the diffusivity coefficient vanishes when $S=0$ or $S=1$, but for the unsaturated region
this term represents a diffusive flux component, which depends on the velocity and degree of saturation.

Assuming that the diffusivity coefficient depends on the resin flow velocity leads to the following equation

$$ \phi \frac{\partial S}{\partial t} + \nabla \cdot (v_R S) = \nabla \cdot (D_V \nabla S) \quad (9) $$

where the diffusivity coefficient $D_V$ includes the velocity effect

$$ D_V = \left( \alpha_m v_R^2 + \frac{\alpha_M}{v_R} \right) \quad (10) $$

In this model, $\alpha_M$ and $\alpha_m$ represent the dispersive coefficients of the macro and microscopic voids respectively, as identified in experimental observations, so that the source term describes the hydrodynamic dispersion due to the presence of pores [1, 5].

Since the main interest is to simulate the flow of the resin phase, Darcy's law and conservation mass for the resin phase have been considered, i.e., Eqs. 1 and 2. Combining both equations for the resin phase yields

$$ \nabla \cdot (K_R(S) \nabla p) = \left( \frac{\phi \mu}{K_{sat}} \right) \frac{\partial S}{\partial t} \quad (11) $$

where the relative permeability depends on the saturation degree. In this study, the permeability model proposed by Breard et al. [5] has been considered. Permeability as a function of saturation is calculated as follows:

$$ R_S \le K_R(S) \le 1 \\ K_R(S) = \left[ (1 - R_S^{1/\beta}) S + R_S^{1/\beta} \right]^\beta \quad (12) $$

where $R_S$ is a fitting factor whose usual values range in the interval [0.4,0.8]. Finally, to obtain a closed description of the flow, we assume the modified saturation equation which governs the transport of the variable S:

$$ \phi \frac{\partial S}{\partial t} + \nabla \cdot (vS) = \nabla \cdot \left( \left( \alpha_m v^2 + \frac{\alpha_M}{v} \right) \nabla S \right) \quad (13) $$

<figure><img src="image_2.png" /><figcaption>Fig. 1 Void content definition</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>Fig. 2 Void content logarithmic functions versus impregnation velocity (m/s)</figcaption></figure>

The saturation S takes a unit value in the saturated domain, a zero value in the empty region and varies between 0 and 1 in the partially saturated region. Note that Eq. 11 for pressure is an elliptic equation coupled to the saturation equation Eq. 13 through the saturation term. On the other hand, the saturation equation is a nonlinear advection equation coupled to the pressure equation through Darcy velocity.

The simulation of the filling process involves the following operations at each time step:

1) calculate the saturation dependent permeability and the source term in Eq. 11;

2) calculate the pressure distribution by applying a standard finite element discretization of Eq. 11;

3) calculate the velocity field from Darcy's law for the resin (1);

4) update the saturation distribution by integrating Eq. 13 using a flux limiter technique.

The boundary conditions are given by:

* The pressure gradient in the normal direction to the mold walls vanishes.

* The pressure or the flow rate is specified on the inflow boundary.

* The pressure is zero in the empty part of mold.

<figure><img src="image_4.png" /><figcaption>Fig. 3 Saturation curve at different impregnation velocity</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>**Fig. 4** Analysis of the saturation calculation and convergence for different mesh grids</figcaption></figure>

The modified saturation equation

For the sake of simplicity from now on we only consider one-dimensional models. In the one-dimensional case Eq. 13 writes:

$$ \frac{\partial S}{\partial t} + \frac{\partial}{\partial x} \left( \frac{v}{\phi} S \right) = \frac{\partial}{\partial x} \left( \frac{1}{\phi} \left( \alpha_m v^2 + \frac{\alpha_M}{v} \right) \frac{\partial S}{\partial x} \right) \quad (14) $$

We use a flux-limiter strategy to approximate numerically the saturation equation. Defining the flux as

$$ F = \frac{v}{\phi} S - \left( \frac{1}{\phi} \left( \alpha_m v^2 + \frac{\alpha_M}{v} \right) \frac{\partial S}{\partial x} \right) \quad (15) $$

Equation 14 can be integrated by applying a second-order upwind scheme preserving the TVD property [6]. The discrete form writes

$$ S_{i}^{n+1} = S_{i}^{n} - \lambda \left( \hat{F}_{i+1/2}^{SW} - \hat{F}_{i-1/2}^{SW} \right) \quad (16) $$

where $\lambda = \frac{\Delta t}{h}$ and

$$ \begin{aligned} \hat{F}_{i+1/2}^{SW} &= \hat{F}_{i+1/2}^{UP} + \frac{1}{2} \chi(r_{i+1/2}) (\mathrm{sign}(v_{i+1/2}) - \frac{\Delta t}{h} v_{i+1/2}) (F_{i+1} - F_i) \\ \hat{F}_{i+1/2}^{UP} &= \frac{1}{2} \{(F_i + F_{i+1}) - \mathrm{sign}(v_{i+1/2}) (F_{i+1} - F_i)\} \end{aligned} \quad (17) $$

Here $h$ represents the mesh size, $\Delta t$ the time step and $\chi(r)$ is the flux limiter function. The superscript UP is

<figure><img src="image_6.png" /><figcaption>**Fig. 5** Comparison of the first and second-order schemes used to solve the saturation problem</figcaption></figure>

associated with the first-order upwind scheme and the superscript SW with the second-order scheme using a modified flux limiter technique (in our case the superbee flux limiter). It is well known that upwinding is an essential part of any numerical scheme for hyperbolic equations. In our construction of numerical fluxes, the upwind direction is determined by the sign of $v$ in the element. The numerical scheme respects the balance that occurs between the flux gradient associated to the advective part and the source term related to the diffusive part when $\partial S/\partial t$ is small. This choice is consistent with implementations employed to simulate numerically hyperbolic conservation laws with source terms, where the upwind contribution of the source term is necessary to the convergence of schemes.

The superbee flux limiter is defined as follows [6]:

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

The coefficient $r_{i+1/2}$ has been defined in Eq. 17 by comparing consecutive variations of the approximate numerical solution for the saturation with respect to the flow direction

$$ r_{i+1/2} = \begin{cases} \frac{S_i - S_{i-1}}{S_{i+1} - S_i} & \text{if } v_{i+1/2} \ge 0 \\ \frac{S_{i+2} - S_{i+1}}{S_{i+1} - S_i} & \text{if } v_{i+1/2} \le 0 \end{cases} \quad (19) $$

Flux limiters are used in numerical schemes to solve problems in science and engineering, particularly in fluid dynamics, described by hyperbolic conservation laws. They are used in high resolution schemes to avoid the spurious oscillations that would otherwise occur with high order spatial discretisation schemes due to shocks, discontinuities or sharp changes in the solution domain. The use of flux limiters, together with an appropriate high resolution scheme, makes the solutions total variation diminishing (TVD). The limiter function is constrained to be equal to zero, if $r \le 0$. Therefore, when the limiter is equal to zero (sharp gradient, opposite slopes or zero gradient), the flux is represented by a low resolution scheme. Alternatively, when the limiter is equal to 1 (smooth solution), it is represented by a high resolution scheme. This prevents the spurious numerical oscillations associated with the conven-
<figure><img src="image_7.png" /><figcaption>Fig. 6 Saturation and pressure profiles for a 1D RTM filling at constant flow rate</figcaption></figure>

tional second-order methods in the presence of discontinuities and reduces the excessive numerical diffusion introduced by the first-order upwind schemes. There are various limiters with different switching characteristics. In this work, the superbee limiter has been used, which combines the stability and accuracy corresponding to first and second order schemes.

## Numerical simulations

In order to analyze the accuracy and efficiency of the proposed numerical scheme, a set of numerical examples are presented. The numerical simulations were carried out in the same conditions and with different constant flow rates. A mold of one-meter length is considered. A constant injection pressure is prescribed ($10^5 Pa$), and the saturated permeability $K_{sat}$ and resin viscosity are set to $10^{-8} m^2$ and $0.1 Pa \cdot s$, respectively. In the numerical simulation, we consider a time step small enough, and a constant value $R_s = 0.4$. The domain is assumed initially empty, except the first element that represents the injection nozzle that is assumed to be fully filled. The dispersive coefficients of the macro and microscopic voids, $\alpha_m$ and $\alpha_M$ are 1 and $10^{-7}$, respectively. The impregnation velocity is usually related to the void content using the dimensionless capillary number ($Ca$), that is

$$Ca = \frac{\mu v}{\gamma \cos(\theta)} \qquad (20)$$

<figure><img src="image_8.png" /><figcaption>Fig. 7 Estimation of the void content by Least Square Method for different impregnation velocities</figcaption></figure>

where $\mu$ is the viscosity of the fluid, $\gamma$ the surface tension at the interface air/resin, $v$ is the fluid velocity and $\theta$ is the contact angle between the resin and the fibers.

The proposed model employs the degree of diffusion in the saturation as a variable to measure the void content. The quantitative relationship between these two parameters is not known. However, the qualitative relationship is clear. In order to evaluate the void content, we define the Void Index as follows:

$$ \text{Void Index}(\%) = \frac{\text{Unsaturated Void}}{\text{Saturated Volume} + \text{Unsaturated Filled Volume}} * 100 $$

where the Unsaturated Void Volume is the volume occupied by bubbles in the unsaturated part or the mold, and the Unsaturated Filled Volume is the volume of resin in the unsaturated volume, as is shown in Fig. 1.

As illustrated in Fig. 2, due to the impregnation of a double scale porous medium, an optimum impregnation resin velocity (or $Ca$) exists, so that the injection velocity controls the void content and a proper optimal filling strategy must be developed. The left side of the curve in Fig. 2 represents (negative slope) the formation of macroscopic voids, whereas the right side (positive slope) represents the microscopic void formation. For this set of experiments, we can conclude that the optimal velocity is about 0.005 m/s.

Results in Fig. 2 have been obtained at constant flow rate in order to demonstrate the influence of the velocity in the Void Content. The experiments were developed at constant injection pressure and for 1D case, then the velocity must be the same in the whole domain in the case of saturated flow. If micro/macro voids appear the velocity will increase in the partiality filled region. The proposed method to optimize the macro/micro content is controlled by the injection velocity and not by the flow front velocity. In the case of constant injection pressure control, the optimized pressure can be obtained by the Darcy's law.

Figure 3 represents, at the same flow front position, the degree of saturation at three different velocities, 0.001, 0.005 and 0.1 m/s using the second order scheme described by Eq. 17. This theoretical flow front position
<figure><img src="image_9.png" /><figcaption>**Fig. 8** Geometry of the part</figcaption></figure>

is defined by the saturated flow, that is $v_i*t$. We can observe that a velocity of 0.005 m/s produces the minimum void content.

To analyze the influence of the mesh size on the simulated results, two meshes with different nodal distributions are considered, consisting of 30 and 150 nodes respectively. The associated numerical solutions are depicted in Fig. 4 (at the left) for a filling time of 50 s, using the first-order scheme (Eqs. 16–17 with $\chi(r)=0$) and the using the superbee flux limiter scheme (Eq. 17 with $\chi(r)$ defined by Eq. 18). We can notice that the convergence is faster when the flux limiter second-order scheme is considered. The errors quantified by the L2-norm for different mesh sizes have been plotted in Fig. 4 (at the right). We can notice that the order of convergence aproximately twice that of the first-order upwind scheme [3].

A comparison of numerical results for the saturation obtained using the first and the flux limiter second-order scheme in the same conditions is illustrated in Fig. 5. Because the diffusion of saturation has important physical

implications, the control of the numerical diffusion is critical. Obviously, two factors contribute to smooth the flow front: one is related to the source term, and the other is a purely numerical effect introduced by the discretization scheme. The classical first order scheme produces spurious numerical diffusion, and then inaccurate results. It can be also appreciated that the numerical diffusion is lower, as expected, when the flux limiter technique is employed. Then the second order scheme becomes a useful tool to obtain accurate results.

Figure 6 represents the saturation (left) and pressure fields (rigth) obtained with the flux limiter technique for different times for a constant injection rate.

Finally, this strategy was introduced in commercial software [7] and the Void Index behaviour adjusted by the least square method as shown in Fig. 7. Numerical results for the Void Index in Fig. 7 have been obtained using the technique described in previous Section. The curve with the negative slope represents the macroscopic void distribution, and the other one the microscopic void distribution. The

<figure><img src="image_10.png" /><figcaption>**Fig. 9** Optimized injection velocity</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>**Fig. 10** Filling times and Void Content for two simulations, one carried out at injection pressure (1bar) and one using the optimized injection strategy</figcaption></figure>

optimal impregnation velocity is 0.005 m/s that yields a Void Index of around 20%.

The following values of physical parameters have been considered:

* Micro Void Index function (Vᵢ_m): Vᵢ_m = 2.031v⁺ 0.1958

* Macro Void Index function (Vᵢ_M): Vᵢ_M = -36.65v⁺ 0.3621

* Optimal velocity: 0.005 m/s

* Capillary coefficient: 0.02

* Optimal capillary number: 0.025

* Preform thickness: 0.005 m

* Permeability: 10⁻⁸ m²

* Viscosity: 0.1 Pa*s

* Porosity: 0.4

A complex part has been simulated as example of void optimization. This example is part of an ecological vehicle for the tourism industry (see Fig. 8). The preform and the resin are, respectively, fibre glass and polyester.

The injection nozzle is a line located along the border of the part. The first simulation has been carried out at constant injection pressure of 1 bar; the impregnation velocity is changing during the mold filling and then the Void content. For the second simulation, as is represented in Fig. 9, an optimized injection velocity has been employed in order to minimize the Void content.

Figure 10 shows the pressure, and void content for a constant injection pressure (left) and optimized injection velocity (right). As observed in the bottom right image the

void content in the second simulation has been reduced to a minimum.

## Conclusions

In the present work, a new numerical procedure is proposed to predict an optimized injection flow rate to minimize micro/macroscopic voids in RTM filling. The macro/microscopic void content are given as a function of the impregnation velocity. The calculated distribution obtained by this new numerical method reproduces experimental observations [5]. For this purpose, a modified advection-diffusion equation describing the evolution of saturation has been discretized by using a flux limiter upwind scheme, which has the ability to limit the extra numerical diffusion introduced by standard first-order schemes. In the numerical simulation of the modified saturation equation two terms contribute to smooth the flow front: one is related to the source term, and the other purely numerical term is introduced by the discretization scheme. The last effect can be reduced by using higher-order numerical schemes. Numerical results confirm that first-order schemes exhibit an excessive and no realistic diffusion due to the numerical approximation of the advective term, while the flux-limiter scheme shows less extra-diffusive effects. Thus, the flux limiter proposed improves significantly the results with respect to the first-order solutions. This technique is able to predict the optimum injection flow rate so as to minimize the
distribution of the micro/macroscopic void content in the finished part.

**Acknowledgments** This research work is supported by a grant from the Ministerio de Educación y Ciencia (MEC), project DPI2007-66723-C02-01.

References

1. Trochu F, Ruiz E, Achim V, Soukane S (2004) New approaches to accelerate calculations and improve accuracy of numerical simulations in liquid composite molding, 7th International Conference on Flow Processes in Composite Materials (FPCM-7), Delaware, USA July 2004

2. Ruiz E, Achim V, Soukane S, Trochu F, Bréard J (2006) Optimization of injection flow rate to minimize micro/macro-voids formation in resin transfer molded composites. Compos Sci Technol 66:475-486

3. García JA, Gascón L, Chinesta F (2010) A flux limiter strategy for solving the saturation equation in RTM process simulation. Composites: Part A 41:78-82

4. Sánchez F, García JA, Gascón L, 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

5. Breard J, Saouab A, Bouquet G (2003) Numerical simulation of void formation in LCM. Composites: Part A 34:517-523

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

7. PAM-RTM 2008 © ESI Group