Control Volume Finite Element Methods for Flow in Porous Media: Resin Transfer Molding

1. Introduction

The method of solution for the own problem is based on the control volume finite element approach which has been shown to be particularly well suited for a fast and efficient implementation of the Newton-Raphson linearization technique.

The basic idea of the control volume finite element approach is to obtain a discretized equation that mimics the governing mass conservation equation locally. A volume of influence, referred to as a control volume, is assigned to each node. The discretized equation for a given node then consists of a term describing the change in fluid mass storage for that volume which is balanced by the term representing the divergence of the fluid mass flux in the volume. The fluid mass flux will depend on the physical properties associated with the volume and the difference in the value of the primary variable between the node in question and its neighbors.

Discretization of the subsurface and the surface flow equations is identical except for the difference in dimensionality. For the sake of clarity, we present in this chapter, a detailed description of the control volume finite element method applied to discretize a simplified prototype continuity equation in Liquid composite molding (LCM). The final discretized equations for all subsurface domains and for surface flow are then presented without providing the details of the derivation.

Liquid composite molding (LCM) processes are routinely considered as a viable option to manufacture composite parts. In this process, a fibrous preform is placed in a mold. The mold is sealed and a liquid thermoset resin is injected to impregnate the fibrous preform. All LCM processes involve impregnation of the resin into a bed of fibrous network. The goal is to saturate all the empty spaces between the fibers with the resin before the resin gels and then solidifies. In RTM, resin is injected slowly and little or no heat transfer and chemical reaction takes place until the mold is filled. Thus, the entire cycle can be viewed as two separate events, fill and subsequent cure.

The mould filling is considered as one of the most critical and complicated stages throughout the entire RTM process. It has a great influence on the performance and quality of the final parts. However, it is hard to understand effects of the filling parameters on the flow front pattern during mold filling. Therefore, it is necessary to understand interrelationship among filling parameters, flow behavior during RTM, and physical properties of the final parts.

The present study concerns the numerical modeling of resin transfer molding techniques (RTM). From a mechanical point of view, these processes can be treated in the same way as the problems of fluids in porous media. Some of authors use methods based on the systems of curvilinear coordinates adapted to a border. However, this approach becomes limited during divisions or fusion of the flow fronts [1-4]. The modeling methods currently elaborated are based on a control volume finite element method (CVFEM). This type of approach was first presented by Wang and others and was adopted in the case of thin shell injection molding [5]. Fracchia was the first to apply the CVFEM to simulate the RTM process [6] and other researchers also followed this approach [7-16]. The application of these methods generate several commercial software: RTM-FLOT ( no longer supported ), PAM-RTM, MyRTM and LIMS.

In the last years, despite the introduction of alternative methods for simulation of flow in porous media BEM (Trochu [17]) and SPH (Krawczak [18]), the CVFEM method has been, usually, used to simulate resin flow in the RTM process. However, this numerical approach has some inherent drawbacks. First, the flow front is difficult to define with the exact location because of using fixed mesh system. This problem in the resin flow front location limits the accuracy of CVFEM method [19]. Mass conservation problems have also been reported with the use of this numerical approach [20,21]. Researchers have addressed these numerical problems and put forward methods to improve the conventional CVFEM method [19, 22].

In this study, the simulation of the resin flow in the RTM process is developed by the control volume finite element method (CVFEM) coupled with the equation of the free surface location. The equation is solved at each time step using nonconforming linear finite elements on triangles, which allow the conservation of the resin flow rate along inter-element boundaries [23]. At each time step, the velocity and pressure in the saturated domain is calculated. The effective velocity is used to update the front position. The filling algorithm determines the time increment needed to fill up completely at least one new element, then the boundary condition is updated and the flow front is advanced for the next iteration. The flow front is refined in an adaptive manner at each time step by using our mesh generator to add new nodes, to get a smoother flow front and reduce the error in the pressure at the flow front of a CVFEM simulation of resin flow in a porous medium. Thus, the position of the flow front, the time-lapse and the rate of the unsaturated zone are calculated at every step. Our results will be compared with the experimental and analytical models in the literature.

On the whole, our study is concerned with the simulation of isothermal filling of moulds in RTM process while adopting the CVFEM and VOF method, taking into account the presence of obstacles and the thickness variation of the reinforcement. The elaborated code allows calculating the position of the injection points and vents, and injection pressure in order to optimize the process parameters. We present a mesh generator for 1D, 2D and 2.5 D geometries for one or several points of injection and vents, and a numerical method for the simulation of the RTM process. Numerical examples are used to validate and assess the applicability of the developed code in the case of anisotropic reinforcements, multilayered, several injection points and the existence of inserts.

2. Mathematical formalism

2.1. Continuity equations and Darcy law

In RTM process and when filling the mould, the resin flow passes through a bed of fibers, The process of injection in the mould is treated as part of the flows of fluids inside a porous medium.

On the basis of partial saturation concept, the mass balance at a point within the domain of an isothermal incompressible fluid flow inside a fiber preform can be expressed as [24] :

$$\varphi \frac{\partial s}{\partial t}=-\nabla .q$$E1

Where q is the volumetric flow rate per unit area, $\varphi$ is the porosity and the saturation level s is 1 for a fully saturated node and its value ranges between 0 and 1 for a partially saturated point. If the transient term on the left hand side of the above equation is removed (saturate case), the following equation for quasi-steady state situation is obtained:

$$\frac{\partial q}{\partial x_i}=0$$E2

As the fluid flows through the pores of the preform, the interstitial velocity of the resin is given by :

$$v_i=\frac{q_i}{\varphi }$$E3

Where $$v_i$$ the intrinsic phase average resin velocity within the pores and $\varphi$ is the porosity of the solid.

Using the assumptions that the preform is a porous medium and that the flow is quasi-steady state, the momentum equation can be replaced by Darcy's law:

$$q_i=-\frac{K_{ij}}{\mu }\frac{\partial P}{\partial x_j}$$E4

where $$\mu$$ is the fluid viscosity, $$K_{ij}$$ is the permeability tensor of the preform, and $$P$$ is the fluid pressure.

Assuming that the resin is incompressible and substituting (4) into (2) gives the governing differential equation of the flow :

$$\frac{\partial }{\partial x_i}\left(\frac{K_{ij}}{\mu }\frac{\partial P}{\partial x_j}\right)=0$$E5

This second order partial differential equation can be solved when the boundary conditions are prescribed. Two common boundary conditions for the inlet to the mould are either a prescribed pressure condition:

$$P_{inlet}=P_{inlet}(t)$$E6

Or a prescribed flow rate condition:

$$Q_n=n_i\frac{K_{ij}}{\mu }\frac{\partial P}{\partial x_j}$$E7

Where $$Q_n$$ is the volumetric flow rate and $$n_i$$ is the normal vector to the inlet.

The boundary conditions along the flow front are as follows:

$$P_{front}=0$$E8

Since the resin cannot pass through the mould wall, the final boundary condition necessary to solve equation (5) is that the velocity normal to the wall at the boundary of the mould must be zero:

$$\overline{v}.\overline{n}=0$$E9

Where n is the vector normal to the mould wall.

3. Discretization of the domain by CV /FEM -VOF method

3.1. Delaunay triangulations

In mathematics and computational geometry, a Delaunay triangulation for a set P of points in the plane is a triangulation DT (P) such that no point in P is inside the circumcircle (figure 1) of any triangle in DT (P). Delaunay triangulations maximize the minimum angle of all the angles of the triangles in the triangulation; they tend to avoid triangles with high aspect ratio [25].

Suppose P = {p₁,….., p_n} is a point set in d dimensions. The convex hull of d+1 affinely independent points from P forms a Delaunay simplex if the circumscribed ball of the simplex contains no point from P in its interior. The union of all Delaunay simplices forms the Delaunay diagram, DT (P). If the set P is not degenerate then the DT (P) is a simplex decomposition of the convex hull of P.

The geometric dual of Delaunay Diagram is the Voronoi Diagram, which consists of a set of polyhedra V₁,..., V_n, one for each point in P, called the Voronoi Polyhedra. Geometrically, V_i is the set of points $$p\in {\Re }^d$$whose Euclidean distance to p_i is less than or equal to that of any other point in P. We call p_i the center of polyhedra V_i. For more discussion, see [26, 27].

The DT has some very desired properties for mesh generation. For example, among all triangulations of a point set in 2D, the DT maximizes the smallest angle, it contains the nearest-neighbors graph, and the minimal spanning tree. Thus Delaunay triangulation is very useful for computer graphics and mesh generation in two dimensions. Moreover, discrete maximum principles will only exist for Delaunay triangulations. Chew [28] and Ruppert [29] have developed Delaunay refinement algorithms that generate provably good meshes for 2D domains.

Figure 1.

A Delaunay triangulation in the plane with circumcircles.

3.2. Discretization domain

In processes such as Resin Transfer Moulding (RTM), numerical simulations are usually performed on a fixed mesh, on which the numerical algorithm predict the displacement of the flow front. Error estimations can be used in the numerical algorithm to optimize the mesh for the finite element analysis. The mesh can also be adapted during mould filling to follow the shape of the moving boundary. In CVFEM, the calculation domain is first discretized using finite elements, and then each element is further divided into sub-volumes. For the discretization of the calculation domain in FEM, we developed a mesh generator (figure 2) allowing to generate 2 and 2.5 dimensional, unstructured Delaunay and constrained Delaunay triangulations in general domains.

Figure 2.

Discretization of calculation domain

3.3. Domain discretization CV/FEM

To use the method CV/FEM coupled with VOF, the mould is first divided into finite elements. Around each nodal location, a control volume is constructed by subdividing the elements into smaller volumes. These control volumes are used to track the location of the flow front.

The calculation domain is in a finite number of triangular elements. After connecting the centroides of the elements with the middles of the elements borders, the calculation domain another time being divided in a number of polygonal control volume, as indicated in figure 3. The borders of any element of the control volume constitute the control surface.

Figure 3.

Discretization of the calculation domain during CV/FEM

3.4. Resin front tracking

The control volumes can be empty, partially filled, or completely filled. The amount of fluid in each control volume is monitored by a quantity called the fill factor. It is the ratio of the volume of fluid to the total volume of the control volume. The fill factor takes values from 0 to 1 where 0 represents totally empty and 1 represents totally full. The control volume method tracks the flow front by determining which control volumes are partially filled and connecting them to form the flow front. The numerical flow front is composed from the nodes with partially filled control volumes as shown in figure 4. The location of the fluid in the control volume cannot be identified, therefore the exact shape of the flow front is not known. Thus, the mesh density can affect seriously the accuracy of numerical solution of the flow front.

For any control volume and after integration equation 5, we obtained the following relationship:

$${\displaystyle {\int }_{cv}\frac{\partial }{\partial x_i}}\left(\frac{K_{ij}}{\mu }\frac{\partial P}{\partial x_j}\right)d\Omega ={\displaystyle {\int }_{cs}\frac{K_{ij}}{\mu }\frac{\partial P}{\partial x_i}}\overrightarrow{n}ds=0$$E10

Where s,Ω, $\overrightarrow{n}$ are the control surface, the control volume and the normal vector of the control surface, respectively. CV and CS represent the control volume and the control surface domains respectively.

Figure 4.

Treatment of the flow front during the fixed meshing method.

Moreover, in every iteration, the calculation matrix contains only elements that have at least one node with a filling ratio unity f=1. This approach requires a rigor during the development of code. However, the time of treatment of the problem is ameliorated.

4. Numerical simulation

4.1. Pressure

During the resolution of the pressure field, we adopted Galerkin’s approximation to represent the distribution of the field of pressure.

Using the procedure outlined by Reddy [30], the finite element formulation of Equation (5) was found to be:

$$\left[K_{ij}^e\right]\left\{P_j^e\right\}=\left\{F_i^e\right\}$$E11

Where

$$K_{ij}^e={\displaystyle \underset{{\Omega }_e}{\int }\frac{K_{\alpha \beta }}{\mu }}\frac{\partial {\psi }_i}{\partial x_{\alpha }}\frac{\partial {\psi }_j}{\partial x_{\beta }}d\Omega$$E12

And

$$F_i^e={\displaystyle {\int }_{{\Omega }_e}f{\Psi }_{i}d\Omega +{\displaystyle {\int }_{{\Gamma }_e}{Q}_n}}{\Psi }_{i}d\Gamma$$E13

Here $${\Omega }_e$$ is the domain of an element. $${\Gamma }_e$$ is the surface of an element, $$P_j^e$$ is the pressure at each node, f is a volumetric source term, $$Q_n$$ is a specified fluid flux through the face of the element, and $${\Psi }_i$$ is a linear interpolation function.

4.2. Velocity

After the pressures are calculated, the velocities are calculated at the centroid of each element using the volumetric flux equation 4 :

$$v_i=-\frac{K_{ij}}{\mu \varphi }\frac{\partial P}{\partial x_j}$$E14

4.3. Calculation of the parameters of filling

The control volume method tracks the flow front by identifying the controls volume partially filled, and connecting them to form the flow front. The numerical flow front is made from the nodes with the partially filled control volumes.

4.3.1. Flow rate calculation

It is assumed that the velocity of the fluid is constant throughout each element (figure.5).

$$Q_{en}={\overline{v}}_e{\overline{a}}_{en}=-{\displaystyle {\int }_o^A\frac{K_{ij}}{\mu }}\frac{\partial P}{\partial x_j}\cdot \overrightarrow{n}ds$$E15

Where $Q_{en}$is the volumetric flow rate in the control volume (n) from element (e), $${\overline{v}}_{en}$$ is the fluid velocity in the element, and $${\overline{a}}_{en}$$ is the area vector for the sub-volume.

Figure 5.

Calculation of the filling velocity (CV/FEM)

4.3.2. Fill factor calculations

After the flow rates in each control volume have been calculated, the fill factors can be updated. Given the current time step, the fill factors from the previous step, the calculated flow rate, and the volume of each CV, the new fill factors can be calculated with:

$$f_n^{i+1}=f_n^i+\frac{\Delta t{\displaystyle \sum {}_e{Q}_{en}}}{V_n}$$E16

where $$f_n$$ is the fill factor, $$\Delta t$$ is the time step, $$V_n$$ is the volume of the control volume, and the superscripts indicate time level.

4.3.3. Time step calculation

The time step for the next iteration must be calculated before the solution can proceed. The optimal time step would be where the fluid just fills one control volume. If a larger step were chosen, the flow front would over-run the control volume and a loss of mass from the system would result. The time to fill the partially filled control volume “n” is calculated with the following relation:

$$\Delta t_n=\frac{(1-f_n)V_n}{{\displaystyle \sum {}_e{Q}_{en}}}$$E17

Once $$\Delta t_n$$ has been calculated for all the partially filled control volumes, the smallest $$\Delta t$$ is chosen as the time step for the next iteration.

5. Adaptive mesh

The numerical schemes used in mould filling simulations are usually based on a time dependent resolution of an unsteady (free surface) boundary value problem. The boundary of the filled area in the mould cavity is constantly evolving, and it is difficult to generate a mesh suitable for all the successive calculation steps of a filling simulation. The fluid front cannot be approximated with a fine precision with an isotropic mesh. Such a mesh would have to be very fine everywhere in the geometrical domain in order to provide an accurate approximation. This would lead to time consuming calculations, although a fine mesh is required only in the vicinity of the flow front and near the inlet gates. For this reason, several researchers have proposed to construct a new mesh of the fluid saturated domain at each time step (Bechet et al. [31] for Eulerian scheme, Muttin et al. [32] for Lagrangian schemes). This approach is long in terms of computer time and rather complex, especially in the case of obstacles, merging flow fronts and for 3D problems Kang and Lee [19] proposed an algorithm, referred to as the Floating Imaginary Nodes and Elements (FINE) method, to get a smoother flow front and reduce the error in the pressure at the flow front of a CVFEM simulation of resin flow in a porous medium. With the FINE method, imaginary new nodes were added at the estimated flow front and the flow front elements were divided into two separate regions: the area of resin and the area of air. Thus, the flow front element was refined in an adaptive manner at each time step. In this study, the generation of the mesh is realized by a code developed by our team, included like a module of the code of numerical simulation. The mesh generator allows to discretize the field in unstructured triangular elements with possibility of local refinement (static and dynamic) and inclusion of inserts. The development of the mesh generator code and its use in our solution allows the refinement of the critical zone (0<f<1) in each iteration.

In the CVFEM process, the numerical flow front is composed from the nodes with partially filled control volumes. So, since the location of the fluid in the control volume is not known, the exact shape of the flow front cannot be identified. In our numerical code, for the zone positioned in the flow front, we have developed a technique of local refinement. This technique uses an iterative method to refine repeatedly an initial triangulation of Delaunay, by inserting new nodes in the triangulation until satisfying size criteria and the shape of the elements (figure.6). The number of calculation points necessary to characterize accurately the deformations of the flow front decreases, and thus, the numerical computing time is reduced.

The refinement of position and shape of the front flow consists in adding the new nodes to the initials meshes (triangulation of Delaunay) in the zone of the flow front. The integration of these refined nodes, in the computer code, is conditioned by the value of the filling rate. The algorithm adopted for the mesh generator makes possible to generate first the standard Delaunay elements and initial nodes for the calculation domain. Then, in a second time, these elements are re-meshed by a technique of addition of nodes (figure.7). The criteria to be respected during refinement are:

The meshing coincides at the interfaces: the meshing technique used implies the coincidence of the grid at the interfaces between the neighbouring elements (of the first standard grid). The triangulation of Delaunay is based directly on the contour nodes discretization which are discretized only once.

An automatic identification of the nodes and the elements of refinement: During the refinement of the meshes created, the initial nodes and elements are identified by a traditional technique of programming called “the coloring”. This technique consists in affecting a particular code to define the elements of the sub-domains. During the resolution of the equations of the linear system, the integration of the nodes and elements resulting from refinement, in the numerical code, is conditioned by the value of the filling rate. Only, the elements with partially filled nodes are taken into account.

Association of the sub-domains to the principal element: During the refinement of the initial elements, the numerical algorithm affects a code to each sub-domain realized by the mesh generator. This code corresponds to the principal element generating the sub-domains. Thus, the mesh code generates a structure of data that permits to associate the sub-domains, resulting from the refinement process, to their principal element.

Figure 6.

Refinement of the triangular element.

Figure 7.

Refinement by inserting new nodes in the triangulation.

In the present approach, a local mesh generator module has been developed and integrated to the principal code of numerical simulation of the RTM process. This approach is useful to carry out the instruction of the meshing at each step time during the execution of the program, which improves calculation time. The meshing technique adopted is based on two concepts. The first one is the restricted triangulations of Delaunay and the second is the Delaunay refinement applied to elements of the flow front (figure.8).

Starting from the reference level of meshing G₀, we define for each element E_1,ps, the condition of refinement as :

$$Condition\Re :\{\begin{array}{l}\exists nod{e}_i\in E_{1,ps}/f_{nod{e}_i}=1\\ \exists nod{e}_j\in E_{1,ps}/f_{nod{e}_j}<1\end{array}$$

Figure 8.

Identification of the elements for refinement

Figure 9.

Example of numerical simulation in RTM process using refinement technique.

The operation for detecting which elements of G₀ must be refined, is repeated at each iteration of the resolution process. If the condition $$\Re$$ has value “true” for an element E_1,ps of the level G₀, the numerical algorithm applies the refinement process and create sub-domains. Thus, the element E_1,ps is divided into an number of sub-domains (7 in figure.9), in each direction to preserve good properties of connection between the domains, which creates a new meshing level G_l. The solutions at the refined G_l are initialized by a prolongation of the same numerical algorithm of the G₀ level. The system of equation is solved successively at each level of grid. Each refined element has a pressure boundary conditions of Dirichlet type. These operations are repeated at each iteration of the numerical resolution until the mould is completely filled.

The advantage of this technique is to increase the local precision of calculations while preserving the properties of a meshing. On the other hand, the resolution process at each meshing level can be carried out starting from the same system of equations, the same approximations and the same solver.

6. Numerical algorithm

Large complex structures need computational models that accurately capture both the geometric and physical phenomena. This may involve the use of flow simulations as warranted by geometry, thickness, and fiber preforms employed. There also exists a need for accuracy improvements by refining the discretization of the computational domain. All these have serious impact on the computational time and power requirements. Physical modeling and computational algorithms and methodologies play an important role in computational times. For large-scale computations, it becomes critical to have algorithms that are physically accurate and permit faster solution of the computational domain. It is not only essential to design efficient parallel algorithms, data structures, and communication strategies for highly scalable parallel computing, it is also very important to have improved computational algorithms and methodologies to further improve the computational performance of large-scale simulations.

In this study, the algorithm adopted uses techniques for the optimization of the execution time. For example, the management of the memory by the dynamic allocation of the tables and matrices allow the optimization of the resources machines. Also, the utilization of the pointers in the definition of the variables of the problems ensures the code the adaptation to the size of the data to be treated. In the same way, the adoption of algorithm of the sparse matrix for the inversion impacts the processing time seriously. Also in every iteration, the calculation matrix is dynamic and contains only elements that have a node with a filling ratio unity f=1. This approach imposes a rigor during the development of code, however, the time of treatment of the problem is ameliorated.

In this study, the numerical adopted is based on the various computational steps involved in adaptive meshing that accurately capture both the geometric and physical phenomena., CVFEM methodology for the update of the filled regions and refinement technique for flow front advancement are summarized below.

Figure 10.

Numerical Algorithm of simulation of filling of moulds in RTM process

7. Results and discussion

7.1. Adaptive meshing

In the first case (figure.11), we present examples prove the accuracy provided by our adaptive meshing technique. These examples treat an injection mold for radial rectangular 400mmx400mm with a uniform thickness of 4mm. The injection is located at the center of the mold. The meshing generates the following data (table.1). The running time is related to a machine CPU 2Ghz Core Duo 2GB RAM

	Using Technique; Adaptive Mesh	Number of principal elements	Number of secondary elements	Number of principal Nodes	Number of secondary Nodes		Time; Consuming
Example A	No	160	0	97	0		12s
Example B	Yes	160	800	97	506		27s
Example C	Yes	160	3520	97	1972		97s

Table 1.

Data of meshing

Figure 11.

Numerical simulation with different type of meshing

In the second case (figure.12), we present two examples for comparison. These examples are related to the gain of the execution time for our adaptive meshing technique. These two examples treat an injection mold for radial rectangular 400mmx400mm with a uniform thickness of 4 mm. The first example is related to a discretization without adaptive meshing. This model has the same level of precision as the second example related to the adaptive mesh (same total number of nodes and elements). The meshing generates the following data (table.2). The running time is related to a machine CPU 2Ghz Core Duo 2GB RAM

	Using Technique; Adaptive Mesh	Number of principal elements	Number of secondary elements	Number of principal Nodes	Number of secondary Nodes		Time; Consuming
Example A	No	3680	0	2069	0		452s
Example B	Yes	160	3520	97	1972		97s

Table 2.

Data of meshing

Figure 12.

Numerical simulation with different type of meshing

7.2. Validation of the results relating to RTM flow with uniform thickness

In the validation examples presented, the moulds used, have a rectangular cavity, with dimensions: (1000 × 200) mm² and (400 × 400) mm², respectively. The first example is relative to the unidirectional validation (1D) of our numerical results, whereas the second is used in radial injection (2D). The fluid viscosity μ = 0,109 Pa.s, the pressure injection P_inj = 2 × 10⁵ Pa, the permeability K = 2,65 × 10⁻¹⁰ m² and the porosity is φ= 0,696.

7.2.1. Unidirectional validation

The mould used has a rectangular cavity of dimensions (1000 × 200) mm², the thickness is uniform and equal to 4 mm. The resin is injected from left side of the mould and the vents are placed on the right-side of the cavity. Under these conditions, the kinetics of the flow can be obtained analytically by the equations of table.3. The comparison of the three kinetics of the front flow obtained is presented in figure 13. It shows a good concordance with the solution obtained.

Pressure	Front’s position	Filling time

Table 3.

Analytical solutions of unidirectional simulations in RTM

Figure 13.

Unidirectional Validation

7.2.2. Bidirectional validation

The mould used has rectangular cavity of dimensions 400x400 mm², and the thickness is uniform 4mm. The resin is injected from the central point of the mould, the reinforcement is isotropic and the permeability is the same one in the various directions. Under these conditions, the form and the position of the flow front can be obtained analytically by the equation 18. The comparison of the two front’s kinetics (analytical and numerical) is presented in figure.14. It shows a good concordance with the solution obtained.

$$\frac{r_f^2}{2}\ln \left(\frac{r_f}{r_0}\right)+\frac{1}{4}\left(r_0^2-r_f^2\right)=\frac{K}{\phi \mu }(P_f-P_{inlet})t$$E19

With r_f is the radius of the front flow in a time t.

7.3. Validation of the results relating to RTM flow with variable thickness

7.3.1. Analytical validation

RTM process can be used to produce pieces with complex geometry. In the industry of the composite, the plates employed often consist of reinforcements with a variable number of plies and stacking sequences. A correct simulation of this process requires taking into account all these parameters. Lonné makes a modeling according to a formalism derived from the Thomson-Haskell method for the prediction of these geometrical variations on the ultrasound transmission [33].

Figure 14.

Bidirectional Validation with K₁₁=K₂₂=K

The reinforcement variation generates different pressures when closing the mold. Under the impact of the compressibility or the relaxation of the mold plates, a variation occurs in the volume and the pores distribution through the fabric and influences permeability and porosity. The mechanical performance of resin transfer molding depend on the fiber volume fraction [34], microstructure of the preform [35–36], void content [37], and impregnation parameters [38]. In most cases, mechanical properties of composite parts can be improved by increasing fiber volume fraction. Higher fiber volume fractions, however, require increased injection pressure and longer time to fill up the mould, which may significantly affect the properties of the final part. Patel et al. [39] molded composite parts containing glass fibers at constant injection pressure.

The study of Chen and al [40] showed that the initial compressibility of reinforcements is essentially related to that of the pores. This compressibility or « relaxation » effect directly influences the global volume and the distribution of the pores. During the mold closing, the variation of the reinforcement thickness generates, under the compressibility or the relaxation effects, a variation of the pores volume and their distribution through the fabric. The works of Buntain and Bickerton [41] were oriented to the way that compressibility affects permeability. Their results clearly showed that permeability (a property required to be perfectly controlled for a correct simulation of the flow front and the distribution of the pressure) was closely related to the pore volumetric fraction. Several models have been proposed to estimate the value of the permeability for various porous media. Capillary models such as those proposed by Carman [42] and Gutowski et al. [43] use the fiber radius and porosity to predict the permeability, but several discrepancies with experimental data have been reported [43–47]. Theoretical models have also been developed for different idealized media structures [44–47]. Most models may not give accurate prediction of permeabilities since fibrous mats used in RTM are often more complex than the idealized unit cell patterns used in theoretical derivations. Thus due to the lack of adequate predictive models, permeability of RTM preforms are usually determined experimentally.

A number of permeability measurement methods have been developed; however, there is no standardized measurement technique for RTM applications. Trevino et al. [48] and Young et al. [49] determined unidirectional permeability using two pressure transducers at each of the inlet and the exit in conjunction with an equation based on Darcy’s law applied to their flow geometry. Calhoun et al. [50] presented a technique based on placing several pressure transducers at various locations inside the mold. Adams and Rebenfeld [51–53] developed a technique that quantifies the planar permeability using the position and shape of the advancing resin front as a function of time. A transparent mold was used to enable the monitoring of the advancing front. Other techniques based also on the observation of the moving resin front are common in the literature [43,51-59]. However, transparent mold walls may not have enough rigidity to avoid deflection, which has been shown to perturbate the measured permeability values [45,48,49,54]. In this context, the University of Plymouth radial flow permeability appartus was enhanced by the use of a laminate of two 25 mm toughened float glass sheets [60] as the upper mould tool to achieve a similar stiffness to the aluminium mould base [61].

In this study, the variation of the plies number (figure.15.a) and the stacking sequence are modeled by the variation of permeability and porosity. During the standard approach, these parameters are defined as an intrinsic property of the global discretized domain. In our approach, the permeability and porosity are defined at the level of the element. The comparison of the two kinetic fronts (analytical and numerical) is presented in figure.15.b ; it shows a good concordance with the solution obtained.

The analytical model for this type of reinforcement is indeed the prolongation of the linear model already adopted in the case of a medium with isotropic permeability. To ensure the accuracy of our numerical results, the elements of the initial meshing belong only to the one of the two zones. Also, in the static refinement of meshing, the creation of the new refined elements respects the condition of the uniformity of the permeability within the element (the element of refinement must belong only to the one of the two zones).

In addition, to illustrate the thickness effect, we present within the framework of the analytical validation, the case of a radial flow through a multi-thickness reinforcement.

In the setting of a bi-dimensional flow, we used the reinforcement with variation of the number of plies. The mould cavity had a uniform thickness.

1. Reinforcement with a constant thickness (10 plies) (Figure 16.a)

2. Reinforcement with half 10 plies and half 20 plies (Figure 16.b)

3. Reinforcement with ¼ 10 plies and ¾ 20 plies (Figure 16.c)

Figure 15.

a) Reinforcement with multiple thickness, b) Evolution of the pressure-position for multiple thickness

For a reinforcement with uniform thickness (figure16.a), the permeability in the principal directions is constant. The shape of the flow front for a radial injection is a circle whose center is the point of injection. For the case of two different reinforcements thicknesses (figure 16.b, and 16.c), the position and the shape of the flow are variable according to the resistance presented by the fabric (permeability).

Thus, for the figure 16.b where each half of the reinforcement is characterized by a fixed value of the permeability, it is quite clear that the solution of the Darcy law concerning the shape of flow front, in each zone, has a form of half-circles spaced according to the value of the permeability. The connection between the two shapes respects the continuity of the flow front. Finally, and in order to ensure the required precision, we proceeded to each step of time, with a refinement in the zone of connection of the two parts. The same approach is adopted for the figure16.c.

Figure 16.

a) Flow Front for a bi-dimensional injection with Reinforcement with a constant thickness, b) Flow Front for a bi-dimensional injection with Reinforcement half 10 plies and half 20 plies, c) Flow Front for a bi-dimensional injection with Reinforcement ¼ 10 plies and ¾ 20 plies

7.4. Simulation of the pressure distribution for a multiple – Thickness reinforcement

The pressure distribution directly influences the injection pressure, the time needed to fill up the mould, the position of the injection points and vents and especially the mechanical properties of the final piece. The approach adopted by our team, use the permeability at the element’s level. The thickness variation of the element directly influences its permeability and its porosity. The integration of these parameter’s variation inside the code of resolution, gives more precision.

The model we treated is a wiring cover made of glass fibers and an Isophthalic Polyester resin. The model has the specificity of including a square insert inside the piece, with the possibility to vary the thickness around the insert (square zone, Figure 17).

Figure 17.

Piece with insert and reinforcement multiple thickness

During the moulding in RTM process, the resin injected infiltrates in empty spaces between fibers. However, a minor modification of the characteristics of the preform in specific places (around the insert for example), can cause significant deviations in the flow and the results on the final properties of the part can be disastrous generating the rejection of the whole process.

In order to simulate the filling process for 3D Dimension, logically a 3D model would be required. Since the thickness of composite parts is often much smaller than its length and width, thin part assumptions can be used for these simulation models [62]. For example, the resin flow in the thickness direction (here denoted as z) is neglected. Therefore, these models, although they describe 3D geometries, they are often called 2.5 D flow models [63].

The figures presented in 18.a, 18.b and 19 clearly show the impact of taking into account the variation thickness on the accuracy of the simulation of the flow front and the distribution of the pressure for 2.5 D models. These elements are paramount to control the parameters of the moulding and to obtain the required properties of the final part. the increase in thickness on a particular area by inserting new plies of reinforcement, leading to decreased permeability and increased porosity.

During numerical resolution, after the pressures are calculated, the velocities are calculated at the centroid of each element. Thus, given the reduced permeability and increased porosity, the velocity decreases while increased the thickness. In Figure 18.b we note the delaying of flow front in this area.

Figure 18.

a) 2D Simulation without taking into account the effect of thickness variation, b) 2D Simulation when taking into account the effect of thickness variation

Figure 19.

2.5D simulation of pressure distribution with thickness variation

The approach adopted during this study ensures flexibility during the simulation of the heterogeneities of the problem. The properties of the medium are calculated at the level of the element of the mesh. Having developed an adaptive mesh generator specific to the team, we are able to ensure a better adaptation of the elements of grid to describe physical specificities generated in each problem. The results obtained on the figures of this paragraph show the relevance of the present approach.

8. Conclusion

During this study, we developed a mesh generator and a numerical code to simulate the filling of an isothermal mould in RTM process, by adopting the CV/FEM and VOF method. An adaptive static meshing to model the variation of the thickness, and dynamic to refine the flow front are used. This approach is useful to carry out the meshing at each step of time during the resolution process, which improves precision and calculation time. The algorithm adopted treats permeability and porosity at the level of the element of the mesh. The effects of the variations of the plies number and the stacking sequence, around the inserts for example, are modeled by the variation of permeability and porosity. The flows around obstacles and through the reinforced area near inserts are simulated by the present formalism with a local refinement of the meshing. The results obtained during the numerical simulation show a good concordance with the results: analytical, experimental and numerical. We can note the effectiveness of the numerical model developed in the predicted flow front and the distribution of pressure. An excellent reproduction of the form of the front and a good precision of its position are obtained. In our next studies, we will be interested in applying a similar approach to simulate saturation effects.