# A flux limiter strategy for solving the saturation equation in resin transfer molding process simulation

J. A. Garcia, L. I. Gascon, Francisco Chinesta

► To cite this version:

J. A. Garcia, L. I. Gascon, Francisco Chinesta. A flux limiter strategy for solving the saturation equation in resin transfer molding process simulation. Composites Part A: Applied Science and Manufacturing, 2010, Flow Processes in Composite Materials, 41 (1), pp.78-82. 10.1016/j.compositesa.2009.09.008. hal-01004832

HAL Id: hal-01004832

https://hal.science/hal-01004832

Submitted on 15 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
# A FLUX LIMITER STRATEGY FOR SOLVING THE SATURATION EQUATION IN RESIN TRANSFER MOLDING PROCESS SIMULATION

J.A. García<fnref n="1" />,<fnref n="3" />, Ll. Gascón<fnref n="1" />, F. Chinesca<fnref n="2" />

**SUMMARY:** The main aim of this work is to propose 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 proposed here considers the flow through a partially saturated medium, including the dependence of permeability on the saturation degree. The model consists of an elliptic equation governing the pressure distribution and a transport hyperbolic equation describing the evolution of the saturation in RTM. A global flux limiter fixed mesh strategy is proposed for solving the transport equation with a source term. The flux limiter method has the ability to limit the extra numerical diffusion introduced by standard first-order schemes. This formulation can lead to improvements of existing RTM flow simulation codes and optimize the injection process. Preliminary numerical results are presented to validate this new approach.

**KEYWORDS:** Resin Transfer Molding (RTM), fixed mesh simulation, transport problems, flux limiter

## 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], a new procedure was proposed

<footnote n="1">¹Universidad Politécnica de Valencia. Camino Vera s/n, 46022, Valencia (Spain) and</footnote>

<footnote n="2">²Laboratoire de Mécanique des Systèmes et des Procédés – ENSAM Paris, 151 Boulevard de l’Hôpital, F-75013 Paris (France)</footnote>

<footnote n="3">³Corresponding author’s Email: jugarcia@mcm.upv.es</footnote>
to integrate accurately the transport of the different above mentioned fluid properties. In this work, the numerical procedure is extended for solving the transport equation governing the evolution of saturation in RTM during mold filling. 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 that depends on saturation. The presence of this source term explains the quadratic spatial pressure distribution in an unsaturated porous medium noticed in unidirectional flows. To derive a closed model the saturation equation is considered whose source term depends on the micro and macro voids content.

The technique here used is based on a flux limiter fixed mesh strategy for solving the transport equation which governs the evolution of degree of saturation. It is well-known that spurious oscillations in the computed solution can be avoided by using first-order upwind discretization schemes, but unfortunately these first order schemes introduce an excessive numerical diffusion. One alternative for avoiding non-physical oscillations, reducing the extra numerical diffusion, lies in the use of second-order numerical fluxes in the regions where the solution is smooth enough that reduce to first-order in presence of discontinuities. This technique was successfully used for computing the volume fraction evolution in [1].

## GOVERNING EQUATIONS

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

$$v = -\frac{K}{\mu} \nabla p \qquad (1)$$

combined with the continuity equation, which should account for the loss of resin due to the perform impregnation. This can be done by incorporating a source term in the mass balance equation

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

where $v$ is Darcy velocity, $K$ is the preform permeability tensor, $\mu$ is the fluid viscosity, $\phi$ is the porosity of the fibrous reinforcement, $p$ is the resin pressure and $S$ the degree of saturation. In order to evaluate the permeability, we make use of the model proposed by Breard et al. [3] to predict 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_{unsat}(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)$$
and $\beta$ is a fitting factor whose usual values range in the interval [0.4,0.8] [3]. Combining Eqn. 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 close the problem, we assume that saturation is governed by the following transport equation [3, 5]:

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

where the source term can take different forms. According to the model proposed by Breard et al. [3], the source terms writes:

$$ C = \alpha v \nabla^2 S - Q_s \qquad (7) $$

or

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

as in the model of Ruiz et al. [5]. In this latter approach, $\alpha_M$ and $\alpha_m$ represent the dispersion coefficients of the macro and micro voids, respectively, and $Q_s$ is a function of the modified capillary number [3].

The saturation $S$ takes a unit value in the saturated domain, a zero value in the empty region and varies between 0 and 1 in the partially saturated region. The whole domain has been designated by $\Omega$ and its boundary by $\partial\Omega$. Initially, we assume the condition

$$ S(x, t=0) = \begin{cases} 0 & x \in \Omega \\ 1 & x \in \partial\Omega^{-} \end{cases} \qquad (9) $$

where the saturation function on the inflow boundary $\partial\Omega^{-}$ is assumed unchanged during the filling process, that is:

$$ S(x \in \partial\Omega^{-}, t) = 1 \qquad (10) $$

The filling process simulation involves at each time step:

1. The calculation of saturation dependent permeability as well as the calculation of the source term of Eqn. (5).

2. The calculation of the pressure distribution by applying a standard finite element discretization of Eqn. (5).

3. The calculation of the velocity field from Darcy's law (1).
4. Updating of saturation is done by integrating Eqn.(6) using a flux limiter technique.

The boundary conditions are given by:

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

* The pressure or the flow rate is specified on the inflow boundary (injection nozzle).

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

## THE ADVECTION EQUATION

The saturation equation has a hyperbolic character, needing for appropriate stabilizations. In this section we propose a flux-limiter technique for its solution. For the sake of simplicity from now on we only consider one dimensional models.

In the one-dimensional case Eqn. (6) reads:

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

and defining the flux as $F = \frac{v}{\phi} S$, Eqn. (11) can be integrated by applying a second-order upwind scheme preserving the TVD property [6], whose discrete form writes:

$$ S_i^{n+1} = S_i^n - \lambda (\hat{F}_{i+1/2}^{SW} - \hat{F}_{i-1/2}^{SW}) + (\Delta t) C_i^n \qquad (12) $$

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}) \left( \operatorname{sign}(v_{i+1/2}) - \frac{\Delta t}{h} v_{i+1/2} \right) (F_{i+1} - F_i) \\ \hat{F}_{i+1/2}^{UP} &= \frac{1}{2} \left\{ (F_i + F_{i+1}) - \operatorname{sign}(v_{i+1/2}) (F_{i+1} - F_i) \right\} \end{aligned} \qquad (13) $$

Here $h$ represents the mesh size, $\Delta t$ the time step and $\chi(r)$ is the flux limiter function. For the approximation of the source term, and considering the model (7)

$$ C_i = \left( \alpha v_i \left. \frac{\partial^2 S}{\partial x^2} \right|_i - Q_s \right) \approx \left( \alpha v_i \left( \frac{S_{i+1} - 2S_i + S_{i-1}}{h^2} \right) - Q_s \right) \qquad (14) $$

or using the expression of the source term given by (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] \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] \qquad (15) $$

The superscript UP is associated with the first-order upwind scheme and the superscript SW with the second-order scheme using a modified flux limiter technique (in our case the Sweby flux limiter). To define the coefficient $r_{i+1/2}$ in Eqn. (13), we propose to define its value
comparing consecutive variations of the approximate numerical solution with respect to the flow direction [6]:

$$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 (16)$$

The Sweby flux limiter reads [6]:

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

In general, $\chi(r) = 0$ if $r \le 0$ what 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, it is necessary that $\chi(1) = 1$ to obtain a second-order scheme in smooth regions of the solution. The limiter function should be chosen verifying both conditions and maximizing the antidiffusive flux. We remark in Eqn. (13) that the choice $\chi(r) = 0, \forall r$ gives rise to the original first-order upwind method, whereas the choice $\chi(r) = 1, \forall r$ defines a fully second-order scheme. Hence, the proposed hybrid scheme, with $0 \le \chi(r) \le 1$ in Eqn. (13), is of second-order in smooth regions of the solution but it is a first-order method when a discontinuity is present. This allows avoid the spurious numerical oscillations associated with the conventional second-order methods in the presence of discontinuities and reduces the excessive numerical diffusion introduced by the first-order upwind schemes.

## NUMERICAL SIMULATIONS

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

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

Fig. 1 proves also that the diffusivity of the scheme decreases with the mesh refinement. It can be noticed that the use of a first-order scheme introduces a significant numerical 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.

Fig. 2 depicts the saturation profile each 100 seconds by applying both the first-order and the Sweby flux limiter scheme when 150 nodes are used in the space discretization. We can appreciate the significant impact of the discretization scheme on the computed solution. It can be also appreciated that the diffusion is lower, as expected, when flux limiter techniques are employed.

<figure><img src="image_3.png" /><figcaption>Fig. 2 Saturation profiles each 100 seconds, using the first-order scheme (left) and Sweby flux limiter (right).</figcaption></figure>

To quantify the convergence different computations were performed using different mesh sizes (30, 60, 90, 120 and 150 nodes, respectively). The computed results illustrated in Fig. 3 represent the partially filled domain to the total length of the mold for both discretization schemes. The pressure distribution varies linearly within each saturated element and is continuous between elements. The saturation and pressure distributions for a filling time of 700 seconds are represented in Fig. 4. It can be noticed the presence of three different behaviors: fully impregnation, unsaturated and empty domains. In the fully saturated regions the model reduces to the usual Darcy model, leading to a linear pressure distribution. In the partially saturated regions, the pressure distribution becomes parabolic, as announced in [3]. Finally, in the empty domain the pressure vanishes.
<figure><img src="image_4.png" /><figcaption>Fig. 3 Convergence analysis.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 4 Pressure and saturation at $t = 700$ seconds (left); pressure distributions (right).</figcaption></figure>

## CONCLUSIONS

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

## ACKNOWLEDGEMENTS

This research work is supported by a grant from the Ministerio de Educación y Ciencia (MEC), project DPI2007-66723-C02-01.
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 (2003) 877-893.

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", *Comput. Methods Appl. Mech. Eng.*, 196 (2007), 2300-2312.

3. J. Breard, A. Saouab, G. Bouquet, "Numerical Simulation of Void Formation in LCM", *Composites Part A: Applied Science and Manufacturing*, 34 (2003), 517-523.

4. E. Ruiz, V. Achim, S. Soukane, F. Trochu, J. Bréard, "Optimization of Injection Flow Rate to Minimize Micro/Macro-Voids Formation in Resin Transfer Molded Composites", *Composites Science and Technology*, 66 (2006), 475-486.

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

6. Sweby P.K., "High Resolution Schemes Using Flux Limiters for Hyperbolic Conservation Laws", *SIAM J. Numer. Anal.* (1984), 995-1011.