A mathematical model for RTM process simulations in geometries with irregular shapes

Composites made by using RTM processes have attractive properties to aerospace and automotive industries. However, it is mandatory optimize all process parameters. Some mathematical models developed to predict fluid flow in porous media may be used as a support tool for laboratory studies. These models present complex sets of differential partial equations that may be solved numerically, via discretization methods. This work presents a mathematical model to predict resin and air flow in RTM processes. The finite volume method was used to discretize the equations written in boundary fitted coordinates, without the need for any front tracking algorithm. To validate the methodology, numerical results for flow front positions in rectilinear flows were compared to analytical data. The model was also employed to describe the fluid flow in geometries with arbitrary irregular boundaries.


INTRODUCTION
Composite materials are used by humans since ancient times.Such materials are obtained from the combination of two or more materials and share the characteristics of its constituents.The possibility of creating a material that meets a market demand is the main advantage of composites, especially for automotive, aeronautical and aerospace applications.
The resin transfer molding (RTM) is a production method characterized by filling the mold that contains the reinforcing material, by means of the injection of a fluid resin.The technique allows the production of large quantities of parts, with good quality finishing, at low cost.The optimization process requires the understanding of the parameters involved during the filling stage: mold geometry, material properties and process parameters, such as temperature, injection flow rate or injection pressure.
Computer simulations are effective tools that can be helpful to support RTM optimization [1,2].Some of the studies conducted to simulate fluid flows during the mold filling stage in RTM processes use one-dimensional approaches (linear or radial injection).In most cases, these investigations are applied to the injection of resin inside rectangular geometries.
Most of the mathematical models for the resin injection in RTM processes consider only one phase, the resin.Such models are quite useful to analyze the macroscopic behavior of the fluid flow and understand some of its important properties.Despite of the use, some limitations of these models can be noticed because air phase within the mold cavity is ignored.Models which are based on two phases (air-resin) take into account the interaction between them, producing better results [3].For more complex problems, including threedimensional models with multiple phases, commercial software are generally preferred, but will imply higher costs.Commercial computational packages also don't allow access to the source code, making difficult or impossible to implement deeper changes to the mathematical and numerical modeling by user.
A wide variety of practical applications lies in multiphase flow of immiscible fluids in porous media, providing an overview around the models used, highlighting the application of certain techniques [15].Some works present a review regarding approaches to model and simulate RTM processes [1,15].
In some works, researchers have used discretization methods that need some specific technique for location of the resin front [4][5][6][7][8][9][10][11][12][13].When using the finite element method (FEM) or finite volume based on elements (MEFVC) some authors usually consider the transient problem as a series of quasi static problems, which implies limitations to the magnitude of the time step used [5].The finite volume method (FVM) is used in fluid flow problems because it ensures the conservation of the physical quantities [14].
Several researchers have preferred the use of Cartesian coordinates to study RTM process.Although the use of Cartesian grids is relatively easy, these grids are limited for the modeling of problems with irregular geometries.In such cases, it is preferable to use boundary fitted coordinates.In this sense, resin injection simulations in RTM processes using the MVF and boundary fitted coordinates are the goal of the present work.

METHODOLOGY 2.1 MATHEMATICAL MODELING
In simulations of RTM processes, the impregnation of the preform is usually modeled as a fluid resin flow through a homogeneous porous media in a macroscopic scale [1].In many cases, the resin injection process may be considered isothermal [17].Here, we use a model for a two-phase (resin-air) immiscible flow with no gravitational and capillarity effects in a homogeneous porous media.The mass conservation equation for any phase p (air / resin) of the system may be given by: (1) where superscript p indicates the phase p, φ is the porosity of the media, ρ p , Z p and are, respectively, density, mass fraction and velocity of phase p, and ρ m is the average density of the mixture.In Equation (1), m p represents the mass flow rate per unit of volume, defined as follows: where q p represents the volumetric flow rate of the phase.
The nature of the resin flow in RTM processes, allows the use of Darcy's law as the momentum equation [1].For horizontal multiphase systems, the mathematical formulation of the Darcy's law for phase p is written as [18]: (3) where K represents the permeability tensor, an intrinsic property of the porous media, λ p is defined as the phase mobility, given by: (4) In Equation ( 4), μ p is the phase viscosity and k rel, p is an important parameter that will appear when multiple phases are present in the system.It represents the relative permeability of the phase, which depends directly on its saturation.In the present work, it was assumed: k rel, r = S r and k rel, a = 1k rel, r .
Darcy's equation may be combined with mass conservation resulting in the following equation: ( where λ ~p is mobility rewritten to include phase density, as follows: Primary variables in system of equations represented by Equation ( 5) are the mass fractions (Z a , Z r ) and pressures (P a , P r ) of the phases air and resin.Neglecting capillary effects, we have P a = P r = P.The closing equation needed to solve this system ensures that the entire pore volume is completely filled by phases (Z a + Z r = 1).

NUMERICAL SOLUTION
Some terms in the governing equations introduce nonlinearities, for instance the phase mobility.So, to find an analytical solution for these equations would be quite impossible.In these cases, numerical solution may be a very interesting alternative.Among the different techniques cited, the finite volume method may be used, once it ensures the conservation of the physical quantities.Boundary-fitted coordinates make mathematical model more flexible to be applied to problems with irregular geometries.

Transformation of the domain
The use of boundary fitted coordinates requires that the nonrectangular mesh in the physical space be transformed to a rectangular hypothetical computational representation for it, as pictured in Figure 1.

Transformation of the governing equations
Governing equations also must be transformed and solved in the computational domain.
After transformation, equations can be written as follows: (7) where: Equations ( 10) have some terms defined as the metrics of the transformation which bring all the grid information to the transformed equations and relate the physical domain to its computational representation, given by: ξ x = Jy η , ξ y = -Jx η , η x = -Jy ξ and η y = Jx ξ [14].For 2D grids, J represents the relation between areas in computational and physical domains and can be calculated as follows: (11)

Integration of the governing equations
The use of the finite volume method requires the integration of the Equation ( 7) in space and time for an elementary volume shown on Figure 2.After integration, equation becomes: (12) where DV = DξDηDg is the volume dimensions on boundary-fitted coordinates system.
) Some terms with partial derivatives are present in the integrated equation.These terms may be approximated using finite differences, as shows the following samples for the west face: (13) (14) Phase mobility in Equation ( 12) must be evaluated in volume faces (w, e, n and s).In the present work it was used the Upwind Differencing Scheme.Calculation of the phase mobility in face e, as an example, is given by the phase velocity in that face calculated by the following Darcy's equation, rewritten in boundary-fitted coordinates, as follows: (15) where: flows in the opposite direction, i. e., , then .In the other faces, a similar idea may be used.More details about the numerical treatment can be found in [19].

Initial and boundary conditions
The boundary condition for the elementary volumes adjacent to the boundaries is no flow, i.e., it is assumed that all the boundaries are impermeable everywhere inside domain.Thus, the interaction with the surrounding is given by specifying the source / sink terms q p used in the governing equation.This condition will occur only for volumes with resin inlets or air vents.At these points, it is assumed that the volumetric flow rates of the phases are proportional to their mobility: (17) where superscript T is the total mobility of phases, resin and air.
For the resin inlets, we consider that only resin phase flows.There are two conditions possible for resin injection: constant flow rate or constant pressure.For constant injection flow rate, we have: q a = 0 (19) For injections under constant pressure, the boundary condition at the resin inlets is: (20) where subscript def indicates a predefined value.At the air vents, we can define total flow rate (resin + air) as the boundary condition as follows: q T = q def (21) According to Upwind scheme, for , mobility in face e will be .If resin Usually, in RTM processes, the condition for the air vent involves the value of the pressure.

VALIDATION WITH ANALYTICAL SOLUTION (REGULAR GEOMETRY)
The first case studied here is related to an unidirectional flow with resin injection under constant flow rate.This problem describes the resin flow through an isotropic fibrous medium [20].Figure 3 shows mold cavity dimensions and how resin flows inside the mold.Simulation parameters are shown in Table 1.
Figure 4 shows the distribution of resin injection points and air vents as well as the grid with 50x10 control volumes that discretizes the geometry.To simulate correctly the resin flow across the borders, points were distributed along all the volumes, as illustrated in Figure 4.The total flow rate in each border is divided by the quantity of volumes.

290
A mathematical model for RTM process simulations in geometries with irregular shapes Figure 5 shows the similarities between numerical and analytical results for the flow front advance.The analytical solution for the flow front position as a function of time can be calculated as [8]: The curve shows that the complete cavity filling occurs approximately at t = 155 s, what can be calculated analytically as follows [8]: (25) Figure 6 illustrates the comparison between the numerical and analytical values for the injection pressure as a function of the time.In this case, analytical results were obtained by the following model [8]:

APPLICATION (IRREGULAR BORDERS)
The geometry with irregular borders used in this case is 2 mm deep, as shown in Figure 7.
The values of the parameters used in the simulation are listed in Table 2. To discretize the domain a boundary fitted grid was generated, containing 44 x 35 elementary volumes, illustrated in Figure 8.This figure also shows schematically the distribution of resin inlets and air vents, located at the two opposite sides of the cavity.Here, as in the first case, total flow rate is divided by the number of volumes.

292
A mathematical model for RTM process simulations in geometries with irregular shapes The evolution of the flow front position may be seen in Figure 9, for three different times.The flow front shape changes according to the border edge.Afterwards, the flow front reaches the end of the cavity unevenly, hitting the corners faster.
The pressure field inside cavity is according to the front position, as can be seen in Figure 10 for an elapsed time of 35 s.From the analysis of those figures we can see that the pressure increases from atmospheric level only after resin reaches that area.Note that pressure level is higher at the location where resin is being injected, i. e., along the border, being even higher in the center of that area.
For investigate how injection pressure changes with time during the whole process, some "virtual sensors" were distributed as shown in Figure 11.Pressure profile along time for both sensors is seen in Figure 12.In the beginning of the simulation, pressure levels for both sensors are similar.As more resin is injected in the cavity, pressure in the sensor 2 (center) becomes higher and this difference persists until the end of the process.
Figure 11 also shows a longitudinal line used to evaluate pressure levels inside mold cavity.Pressure profiles along this line may be seen in Figure 13, for three elapsed times (5, 35 and 85s).These curves show that pressure decreases linearly from the injection points (x = 0 m), where pressure is maximum, until atmospheric value in regions where there is only air.For the others time values, curve is not completely straight, indicating that pressure profile may depend on the local geometry shape.

CONCLUSIONS
This work presented a mathematical model based on the Darcy's law and its numerical solution to simulate a two-phase flow (air-resin) through an isotropic porous media.The set of governing equations was discretized by the finite volume method using boundary fitted coordinates.The solution for the two unknowns of the problem, pressure and mass fraction, was obtained implicitly via Newton's method.Some numerical results for unidirectional flow were validated with analytical model.Methodology was also able to simulate resin flow within irregular shaped geometries, predicting estimation of the flow front shape, position evolution in time and pressure levels inside the mold cavity during the filling stage.We conclude that the model was able to predict correctly resin flow behavior inside the porous media.The use of boundary fitted coordinates allowed application of the model to geometries with irregular shapes.

1 :
Figure 1: Schematic diagram of the transformation from physical to computational domain

Figure 2 :
Figure 2: Elementary volume for integration

Figure 3 :Figure 4 :
Figure 3: Geometry for the unidirectional flow

Figure 5 :Figure 6 :
Figure 5: Numerical and analytical results of the resin front position along the time

Figure 7 :Figure 8 :
Figure 7: Geometry with irregular shape for the multidirectional flow