# Models for evaluating the saturation in Liquid Composites Molding: a numerical journey

Francisco Chinesta, Ll Gascon, Juan Antonio Garcia

► To cite this version:

Francisco Chinesta, Ll Gascon, Juan Antonio Garcia. Models for evaluating the saturation in Liquid Composites Molding: a numerical journey. 9e Colloque national en calcul des structures, CSMA, May 2009, Giens, France. hal-01418159

HAL Id: hal-01418159

https://hal.science/hal-01418159

Submitted on 16 Dec 2016

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
# Models for evaluating the saturation in LCM: a numerical journey

F. Chinesta¹, Ll. Gascon², J.A. Garcia²

¹ EADS Corporate Foundation International Chair - GEM UMR CNRS ECN
1 rue de la Noe, BP 92101, F-44321 Nantes cedex 3, France
Francisco.Chinesta@ec-nantes.fr

² Universidad Politecnica de Valencia
Camino de vera s/n, E-46071 Valencia, Spain
jugarcia@mcm.upv.es

**Abstract** – The main aim of this work is the proposal of a numerical procedure for solving the saturation equation in RTM process simulation. In order to analyze in more detail the progressive impregnation of a fibrous preform by a fluid resin, the numerical model here proposed considers the flow through a partially saturated medium, including the dependence of permeability on the saturation degree. The main numerical issues will be discussed.

**Keywords** – LCM, RTM, Saturation, Numerical simulation

## 1 Introduction

Liquid Composites Molding (LCM) processes are based on a proper resin impregnation of the reinforcement. Modeling and simulation play an important role in the development and optimization of molds production and in devising appropriate resin injection strategies. Minimization of mold filling time without losing part quality is an important issue in Resin Transfer Molding (RTM). Inadequate injection strategies tend to create macro and micro-voids in the part, the formation of which depends on fluid flow velocity.

Over the last decades, different numerical techniques have been used to model the mold filling process in RTM. In general, RTM process simulation requires an accurate treatment of the advection equation which governs the evolution of different fluid properties: fluid presence function, incubation time, temperature, concentration of reacting components, degree of saturation, etc. In some of our former works [1][2] we proposed a new procedure to integrate accurately the transport of different fluid properties. In this work, that numerical procedure is extended for solving the transport equation governing the evolution of the saturation. In this model the permeability is assumed to be a function of saturation [3], and then the continuity equation that governs the pressure distribution, includes a source term which depends on the saturation. The presence of this source term explains the quadratic spatial pressure distribution in an unsaturated porous medium noticed in unidirectional flows. To derive a closed model the saturation equation is considered, whose source term depends on the micro and macro voids content.
## 2 Governing equations

The model which describes the mould filling process is given by Darcy's law:

$$v = -\frac{K_{rel,s}}{\mu} \nabla p \qquad (1)$$

where $v$ is the velocity field, $p$ is the pressure, $K$ is the preform permeability tensor, $\mu$ is the fluid viscosity, and $K_{rel,s}$ denotes the relative permeability with respect to the saturated one. This modified Darcy law should be combined with the continuity equation, which should account for the loss of resin due to the perform impregnation. This can be done by incorporating a source term in the mass balance equation

$$\nabla \cdot v = -\phi \frac{\partial S}{\partial t} \qquad (2)$$

where $\phi$ is the porosity and $S$ the degree of saturation. In order to evaluate the permeability, we make use of the model proposed by Breard et al. [3] for predicting void formation in LCM. Then, the permeability in Darcy's law should be considered as the product of the geometrical permeability, $K_{sat}$, and a relative permeability $K_{rel}$ ($S$),

$$K_{rel,s}(S) = K_{rel}(S) K_{sat} \qquad (3)$$

where the relative permeability depends on the saturation degree

$$K_{rel}(S) = \left[ (1 - R_S^{1/\beta})S + R_S^{1/\beta} \right]^\beta \\ R_S \le K_{rel}(S) \le 1 \qquad (4)$$

where $\beta$ is a fitting factor whose usual values range in the interval [0.4,0.8] [3].

Combining Eqs. (1) and (2) yields

$$\nabla \cdot (K_{unsat}(S)\nabla p) = \left(\frac{\phi\mu}{K_{sat}}\right) \frac{\partial S}{\partial t} \qquad (5)$$

Finally, to make a closed description of the flow, we assume the saturation equation proposed in [3][4], which governs the transport of the variable $S$:

$$\frac{\partial S}{\partial t} + \frac{v}{\phi} \cdot \nabla S = C \qquad (6)$$

where different expressions have been proposed for the source term $C$. Thus, Breard et al. [3] proposed:

$$C = \alpha v \nabla^2 S - Q_s \qquad (7)$$
Another widely used expression writes [4]:

$$C = \left[ \left( \alpha_M \nu + \frac{\alpha_m}{\nu} \right) (\nabla^2 S) - Q_s \right] \quad (8)$$

that incorporates two dispersive coefficients $\alpha_M$ and $\alpha_m$ of the macro and micro voids, respectively, and $Q_s$ is a function of the modified capillary number [3]. This expression reflects that high velocities tends to reduce mainly the macroscopic porosity originating microscopic filling defects whereas low values of the velocity tends to reduce the microscopic porosity because the dominant capillarity creating macroscopic filling defects. This mechanism is illustrated in Fig. 1.

<figure><img src="image_2.png" /><figcaption>Fig. 1. Void formation as a function of the injection flow speed.</figcaption></figure>

### 3 Discrete form of the advection equation

Eq. (6) can be integrated by applying a second-order upwind scheme preserving the TVD property [5], whose discrete form writes:

$$S_i^{n+1} = S_i^n - \lambda \frac{v_i^n}{\phi} (\hat{S}_{i+1/2}^{SW} - \hat{S}_{i-1/2}^{SW}) + (\Delta t) C_i^n \quad (9)$$

where $\lambda = \frac{\Delta t}{h}$, h represents the mesh size, $\Delta t$ the time step,

$$\hat{S}_{i+1/2}^{SW} = \hat{S}_{i+1/2}^{UP} + \frac{1}{2} \chi(r_{i+1/2}) \left( \mathrm{sign}(v_{i+1/2}) - \frac{\Delta t}{h} \frac{v_{i+1/2}}{\phi} \right) (S_{i+1} - S_i) \quad (10)$$

$$\hat{S}_{i+1/2}^{UP} = \frac{1}{2} \left\{ (S_i + S_{i+1}) - \mathrm{sign}(v_{i+1/2}) (S_{i+1} - S_i) \right\}$$

and $\chi(r)$ is the flux limiter function. The Superbee flux limiter reads [5]:
$$ \chi_{SB}(r) = \max\{0, \min\{2r, 1\}, \min\{r, 2\}\} \qquad (11) $$

where r is defined by comparing consecutive variations of the approximate numerical solution with respect to the flow direction [5]:

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

For the approximation of the source term, and considering the model (8):

$$ C_i = \left[ \left( \alpha_M v_i + \frac{\alpha_m}{v_i} \right) \left. \frac{\partial^2 S}{\partial x^2} \right|_i - Q_S \right]_i \approx \left[ \left( \alpha_M v_i + \frac{\alpha_m}{v_i} \right) \left( \frac{S_{i+1} - 2S_i + S_{i-1}}{h^2} \right) - Q_S \right]_i \qquad (13) $$

In general, $\chi(r)=0$ if $r \le 0$ that guarantees that the scheme will be of first-order in the neighbourhood of a discontinuity, since $r \le 0$ implies that the slopes of the solution have opposite signs. On the other hand, a second-order scheme results if $\chi(r)=1$. The flux limiter should be chosen verifying both conditions and maximizing the anti-diffusive flux. Obviously, if $\chi(r)=0, \forall r$, the proposed scheme reduces to the usual first-order upwind method, whereas the choice $\chi(r)=1, \forall r$, defines a fully second-order scheme. Hence, the proposed hybrid scheme is of second-order in smooth regions of the solution but it reduces to first-order method in the neighbourhood of discontinuities, avoiding the characteristic spurious numerical oscillations related to usual second-order methods in the presence of discontinuities, and at the same time reducing the excessive numerical diffusion induced by first-order upwind schemes.

# 4 Numerical results

In order to analyze the accuracy and efficiency of the proposed numerical scheme just described for the discretization of the transport equation governing the evolution of the saturation, where the source term is given by Eq.(8), we analyze some one-dimensional numerical experiments. A mold of 0.5m length is considered. A constant injection pressure is prescribed ($10^{-5}$ Pa) being the saturated permeability $K_{sat}$ and the resin viscosity of $5.10^{-11}$m$^2$ and 0.1Pa.s, respectively. For the numerical simulation, we consider a time step of 0.1 seconds, $\alpha_M = \alpha_m=10^{-10}$, $\beta=0.8$ and a constant value $R_S=0.8$. The domain is assumed initially empty, except the first element that represents the injection nozzle that is assumed full-filled.

To analyze the influence of mesh size on the simulated results, three uniform meshes are considered, consisting of 30, 60 and 150 nodes, respectively. The associated numerical solutions are depicted in Fig. 2. for a filling time of 700 seconds, using the first-order scheme (i.e. $\chi(r)=0$) and the Sweby flux limiter scheme. We can notice that the convergence is faster when the flux limiter second-order scheme is considered.
<figure><img src="image_3.png" /><figcaption>Fig. 2. Numerical results for the saturation with the first-order scheme (left) and with the Sweby flux limiter (right).</figcaption></figure>

Fig. 2. proves also that the diffusivity of the scheme decreases with the mesh refinement. It can be noticed that the use of a first-order scheme introduces a significant numerical over-diffusion. Obviously, there are two terms that contribute to the front smoothing, the one related to the source term, and the one purely numerical introduced by the discretization scheme. The last one can be reduced by using higher-order schemes and fine enough meshes.

Obviously, if the flux limiter is activated during the simulation the convergence rate cannot exceed a unit value, however, because the higher accuracy in the treatment of regions in which the solution is evolving smoothly (full filled domain) the expected accuracy of second order technique should be better than the one obtained with a first order upwinding. Figure 3 proves this expected behaviour.

<figure><img src="image_4.png" /><figcaption>Fig. 3. Convergence analysis</figcaption></figure>

The pressure profiles exhibit a linear part in the full-filled regions and a non-linear one in the unsaturated domain located in the flow front neighbourhood. Figure 4 illustrates this behaviour.
<figure><img src="image_5.png" /><figcaption>Fig. 4. Pressure profiles</figcaption></figure>

**5 Conclusions**

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

**References**

[1] J.A. García, Ll. Gascón, F. Chinesta, A Fixed Mesh Numerical Method for Modelling the Flow in Liquid Composites Moulding Processes Using a Volume of Fluid Technique, *Computer Methods in Applied Mechanics and Engineering*, 192/7-8, 877-893, 2003.

[2] F. Sánchez, J.A. García, Ll. Gascón, F. Chinesta, Towards an efficient numerical treatment of the transport problems in the Resin Transfer Molding Simulation, *Computer Methods in Applied Mechanics and Engineering*, 196, 2300-2312, 2007.

[3] J. Breard, A. Saouab, G. Bouquet, Numerical simulation of void formation in LCM, *Composites Part A: Applied Science and Manufacturing*, 34, 517-523, 2003.

[4] F. Trochu, E. Ruiz, V. Achim, S. Soukane, New Approaches to Accelerate Calculations and Improve Accuracy of Numerical Simulations in Liquid Composite Molding, *International Conference on Flow Processes in Composite Materials*, Delaware, USA, July 7-9, 2004.

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