Numerical simulation of three-dimensional mold filling process in resin transfer molding using quasi-steady state and partial saturation formulations

Composites Science and Technology 62 (2002) 861–879 www.elsevier.com/locate/compscitech

Numerical simulation of three-dimensional mold filling process in resin transfer molding using quasi-steady state and partial saturation formulations A. Shojaeia,*, S.R. Ghaffariana,*, S.M.H. Karimianb a

Polymer Engineering Department, Amirkabir University of Technology, Tehran 15875-4413, Iran Aerospace Engineering Department, Amirkabir University of Technology, Tehran 15875-4413, Iran

b

Received 25 July 2001; received in revised form 3 December 2001; accepted 17 December 2001

Abstract This paper presents numerical simulations of a three-dimensional isothermal mold filling process based on the concept of control volume finite element method. Two numerical schemes are employed to track the flow front: one is based on the quasi-steady state formulation and the other is based on the formulation of partial saturation at flow front. The resulting codes named RTMS and RAPFIL, respectively, allow the prediction of flow front and pressure field in the three-dimensional molds with complex geometries. Since the preform permeability may be a function of fluid velocity, the proposed numerical schemes have accounted for velocitydependency of permeability. This introduces complexity in the numerical schemes due to the correlation between the preform permeability and the fluid velocity during the mold filling. The validity of the two schemes is evaluated by comparison with analytical solutions for simple geometries, and excellent agreements are observed. To illustrate the applicability of the developed codes, the mold filling process of some mold geometries is simulated. A close agreement is observed between the results obtained by both codes in all cases. The CPU time required for the RAPFIL is significantly less than the RTMS as a conventional method on a personal computer and hence, the RAPFIL shows to be quite CPU time-effective for the simulation of complicated three-dimensional geometries. # 2002 Published by Elsevier Science Ltd. Keywords: E. Resin transfer molding; Numerical simulation

1. Introduction Resin transfer molding (RTM) is a process which has found wide applications in polymer composite manufacturing due to its capability to produce parts at low cost, high performance and high volume. Moreover, RTM provides the advantages of low tooling costs, fabrication of large composite parts, highly loaded fiber content with a different form of fiber preform, closed mold and low emission of volatile materials, good potential for automation and so on [1,2]. For these reasons, this process has gained the attention of composite manufacturers especially in automotive industries [3,4]. All operations of the RTM process cycle can be divided into two categories: open mold operations such as mold cleaning, release agent applying, gel coat spraying, fiber * Corresponding authors. Tel./fax: +98-21-6468243. E-mail addresses: [email protected] (A. Shojaei), [email protected] (S.R. Ghaffarian).

loading and part finishing, and the other is closed mold operations involving mold filling and curing stages. In the closed mold operations, the thermally activated resin is injected through a gate or gates until the mold is fed up with resin, which is subsequently solidified by crosslinking reaction. Although the open mold operations may occupy a significant time, the closed mold operations have high potential for controlling the cycle time, performance of the produced part and cost. To achieve the full mechanical properties of fiber reinforced composite, every individual fiber in the preform must be completely wetted by the liquid resin and the void content must be reduced to its lowest possible level [5–7]. The fiber wetting depends on material properties, surface energy of the resin/fiber system, and their contact time [8–10]. In addition to the wetting problem, the macro voids or dry spots may occur due to flow front coalescence or race tracking. Pressure field and flow tracking during mold filling can also have a major role on mold design such as material selection, clamping force level, injection gates

0266-3538/02/$ - see front matter # 2002 Published by Elsevier Science Ltd. PII: S0266-3538(02)00020-9

862

A. Shojaei et al. / Composites Science and Technology 62 (2002) 861–879

Nomenclature A, A0 b, b0 a, b A, B, C f f g h [J1] [K] kij L N n p pinj pvac Q q qinj S [S] r, s, t tf tff
t

stiffness matrices force vectors constants in Eq. (A6) subsurfaces of a subcontrol volume fill fraction of a control volume fill fraction field gravity vector cavity thickness inverse of Jacobian transformation matrix permeability tensor components of permeability tensor mold length shape function outward unit normal surface vector gage pressure injection pressure vacuum level resin volume flow rate injection flow rate surface matrix of derivatives of the shape function local coordinates mold filling time time required to reach the flow front to a specified location time step size

and vent locations. Inappropriate flow can lead to creating a dry spot, take a long filling time, create high molding pressure, deform the mold, wash out the fiber preform and deform the fiber preform near the gate. These manufacturing variables are significantly influenced by fluid dynamics of the resin within the fiber preform. The fluid dynamics are largely affected by mold geometry, material properties and processing parameters. The first two are essentially dominated by part design while the last one is considered as a manufacturing window. The material properties include resin rheology and fiber architecture. The processing parameters involve injection strategies such as constant injecting pressure or flow rate. Furthermore, in nonisothermal processes, temperatures of inlet resin, fiber and mold walls can be considered as processing parameters. In such cases, viscosity of the resin can be affected by the temperature and reaction kinetic. It is, therefore, essential to optimize the processing conditions before fabrication of the mold. Optimization of the conditions is often expensive and time-consuming using the experimental analysis. Alternatively, computer simulation offers important information for design of tooling and prediction of processing conditions and

!

u u, v, w V Vc W [W] X, Y, Z x,y, z

velocity vector components of velocity vector volume volume of cavity width of cavity gravity vector flow front locations global Cartesian coordinates

Greek letters
,  parameters of Eq. (31)  viscosity of resin  porosity  resin density ! saturation  convergence criteria  flux term Subscripts and superscripts c current value within the iteration cycle i node or control volume j iteration level o old value or the value at previous time step Symbols ip Re SCV SS

integration point Reynolds number subcontrol volume subsurface of subcontrol volume

provides an efficient tool for slashing design and production costs. The objective of this paper is to simulate the threedimensional mold filling process including velocitydependent permeability. To this end, two computer codes are developed based on control volume finite element method (CV/FEM). At one, the flow tracking is implemented based on quasi-steady state formulation and in the other one, it is performed using the formulation on the basis of nodal partial saturation. The validation of both codes is evaluated by comparison with analytical solutions, and their applications are demonstrated by simulating some mold filling examples. Also, the results of both codes are compared in terms of closeness and CPU time.

2. Previous work Over the past decade, an extensive effort has been made to simulate the RTM process [11–21]. A review on the modeling and simulation approaches in RTM process was presented by Shojaei et al. [22]. In the most simulations, resin flow has been modeled as a two-dimensional problem. In such a case, the flow is ignored in thickness

A. Shojaei et al. / Composites Science and Technology 62 (2002) 861–879

direction and assumed to be a two-dimensional planar one. Since fiber preform may comprise a number of reinforcement layers stacked at different orientation, an averaged thickness permeability must be find to be included in the flow model. An expression for averaging of permeability in the thickness direction has been given by Bruschke et al. [23]. They suggested the applicability of their expression for behind the flow front where transverse flow between the layers does not occur. In nonisothermal processes, the heat flux from the mold wall cannot be eliminated due to strong presence of diffusive heat transfer. Thus, even in a two-dimensional flow modeling, the temperature distribution must be obtained using the three-dimensional heat transfer model involving convective term in planar direction and diffusive term in thickness direction. This is often called a 2.5-D simulation. In such a case, the temperature dependent variable in flow model, viscosity, must be averaged in the thickness direction of the cavity. Although two-dimensional model gives computationally fast solution, this assumption has reasonable results only for thin shell cavity, where the out of plane flows can be ignored. In practice, the parts may have thick sections and complicated geometries. In these cases or even in thin shell structures which have extremely low permeability in the thickness direction, there will be a significant out of plane flow. For these situations, twodimensional assumption is not able to accurately model the flow field and hence, three-dimensional simulation must be performed to consider the flow progression in all three directions. Young et al. [24], Young [25], Trochu et al. [26] and Lim and Lee [27] developed threedimensional simulations to study the filling process under isothermal and nonisothermal conditions. A major barrier to further development of three-dimensional simulation is its high computation time relative to twodimensional one. Thus, every method which speeds up the computation time can help to develop the applicability and practicality of three-dimensional simulations. Friedrichs and Guceri [28] used a hybrid numerical technique to model three-dimensional flow fields in RTM process. In order to achieve computational economy, they used a hybrid two/three-dimensional solution technique to simulate three-dimensional flow fields. According to them, three-dimensional flow field is only considered wherever all three velocity components are important. Maier et al. [29] presented a fast numerical method for isothermal flow in RTM. Their method was on the basis of adding a single row and column to the cholesky factorization of stiffness matrix derived from a finite element formulation for the pressure field. The numerical treatment of a moving boundary problem is an important issue in the simulation of mold filling. Although a variety of standard numerical techniques is available to simulate the process, but they need some additional computational efforts to solve the mold

863

filling stage due to the moving boundary nature of the problem. In a broad sense, the numerical procedures fall into two categories: (1) moving grid scheme, and (2) fixed grid scheme. A brief review of these schemes can be found elsewhere [14,22]. The moving grid scheme often gives accurate representation of flow front position. But, since the calculation domain should be remeshed at each time step, it is computationally very time consuming and very difficult for the mold with inserts or multiple injection gates. Coulter and Guceri [30] and Gauvin and Trochu [31] have proposed numerical simulations based on boundary fitted finite differences. Um and Lee [32] and Yoo and Lee [18] developed the boundary element method for mold filling process. A numerical simulation based on finite element method was proposed by Chan and Hwang [13,33] using quadrilateral elements to study isothermal and nonisothermal filling but the applications were limited to rectangular molds. These numerical methods were based on a moving grid scheme. Trochu et al. [14] developed a computer program using nonconforming finite elements to simulate the filling process. The flow progression in their method was on fixed grid scheme. Among the numerical techniques the method based on finite element (FE) and control volume (CV) formulations has become the most versatile and computationally efficient way to solve the filling process. This is an attractive method to simulate mold filling of a cavity with complex geometry including inserts. In this method one does not need to regenerate the mesh, since a fixed grid scheme is followed. Bruschke and Advani [11] simulated mold filling process using FE/CV formulations. They employed Galerkin finite elements to solve the pressure in the solution domain and used a control volume approach to find the flow front location at each time step. Young et al. [24] and Lin et al. [34] employed control volume finite element method (CV/ FEM) initially developed by Baliga and Patankar [35]. In these methods, a procedure based on quasi-steady concept is used to track the flow front. Therefore, time step size must be chosen so that the steady state assumption is satisfied at each step. Recently, some researchers have attempted to include the transient term in the formulation enabling one to remove the restriction in selection of time step size [20,36,37]. Voller and Peng proposed an iterative algorithm on the basis of control volume finite element formulation to search the flow front locations at isothermal condition [37] and validated their algorithm by comparison with analytical solution for 1-dimensional flow and with the previously-presented results especially for 2-D flow. However, their algorithm cannot be used to compute the pressure distribution during mold filling, therefore, velocity or pressure dependent parameters such as velocity-dependent permeability or shear-ratedependent viscosity cannot be included in the algorithm. Lin et al. [20] developed a finite element formulation based on the concept of partial saturation at flow front and

864

A. Shojaei et al. / Composites Science and Technology 62 (2002) 861–879

implicit time integration to simulate the isothermal filling process.

Eq. (3) can be rewritten in tensor notation as: !

u ¼

1 ½K ðrp  gÞ 

ð4Þ

3. The governing equations During mold filling, resin flows through a fiber bed. The filling process can be studied based on a microscopic or a macroscopic scale. Since the microscopic approach deals with the complicated local flow field between the fiber bundles, the macroscopic approach has become a common approach for simulation of the RTM process. In macroscopic approach, the mold filling process is regarded as a process of fluid flow through a porous medium.

The permeability tensor of the preform, [K], can be isotropic or anisotopic and is a function of fiber type and porosity. Also, it may be a function of fluid velocity [12]. The Darcy’s law applies to Newtonian fluids, ignores inertia effect and gives plug flow inside the cavity. Substituting the Darcy’s law into continuity equation results in an equation which can be used to solve the pressure field. This equation has been successfully used to simulate the flow of resin in a porous fiber bed, and its validity has been proven by many researchers [2].

3.1. The continuity equation On the basis of partial saturation concept, the mass balance at a point within the domain for an isothermal incompressible fluid flow inside a fiber preform expressed as [20]: 

@! ! ¼ r u @t

ð1Þ

in which the saturation level, !; is 1 for a fully saturated point 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, the following equation obtains for quasi-steady state situation: !

r u ¼ 0

ð2Þ

This is the governing equation for quasi-steady state formulation. 3.2. Darcy’s law It has been observed that liquid flow through porous media can be expressed by Darcy’s law, which has reasonable agreement with experimental results for low Reynolds number [38] defined as Re=u k1/2/, where k is the characteristic permeability of a porous medium and the other terms have been defined in the Nomenclature. Because of having low Reynolds number in the RTM process (Re < 1) Darcy’s law has also become the most commonly used equation for describing the flow within the mold, and is used as the momentum equation. This equation, when gravity is taken into account, is expressed as: 2

3 2 kxx u 1 4 v 5 ¼  4 kyx  k w zx

kxy kyy kzy

2 3 @p  g x7 3 kxz 6 6 @x 7 @p 7 kyz 56 6  gy 7 6 @y 7 kzz 4 @p 5  gz @z

ð3Þ

4. Numerical solution method In the present paper, CV/FEM is used for the numerical solution of the problem. This method is the most popular to solve the filling stage because of their simplicity in handling the moving boundary problems. In CV/FEM, a fixed grid approach is used in which there is no need to regenerate the mesh during flow progression. This makes the simulation to be rapid and effective for complicated geometries compared to moving grid approach. Solution domain is first divided into a finite number of elements enabling a well discretization of the complex geometries and then, the governing equation is discretized using the control volume method. The combination of geometric flexibility of finite element method with conservation property of the control volume technique makes the CV/FEM as a convenient numerical method for solving the moving boundary problem as well as the flow front progression during filling process. To begin the numerical formulation, the mold threedimensional cavity is divided into linear, hexahedral elements defined by eight nodes. A single hexahedral element with its local coordinate system is shown in Fig. 1(a). Control volumes are constructed by joining the centers of each face of the elements to the midpoint of edges and to the center of the elements creating polygonal control volumes. Therefore, each element will contain eight control volume octants from eight different control volumes. Each octant is called subcontrol volume (SCV) and its surfaces are called subsurfaces (SS). The control volumes which completely surround the interior nodes, are comprised of eight SCVs, however, the boundary nodes which are not completely surrounded by the control volumes, are comprised less than eight SCVs, and lie on a control volume face, edge or corner. A typical element divided into its SCVs with one removed is shown in Fig. 1(b). The six faces of each SCV are divided into two groups, those that are coincident with

A. Shojaei et al. / Composites Science and Technology 62 (2002) 861–879

zðt; s; rÞ ¼

8 X

Ni zi

865

ð7Þ

i¼1

p¼

8 X N i pi

ð8Þ

i¼1

where the Ni are shape functions, xi, yi, zi, and pi are the nodal values of x, y, z, and p at nodes 1–8 for i ranging 1–8 in an element [39]. The linear shape function Ni is equal to 1 at node i and vanishes on the other nodes in the element. Nodal fill fraction, f, is defined here to represent the status of each control volume within the entire mold. The fill fraction for each control volume represents the ratio of occupied volume by the resin to its total pore volume. For an empty control volume f=0, and f is 1 when the control volume is completely filled with resin. The completely filled control volumes are considered as a solution domain and the pressure is calculated at those, while in all partially filled or empty control volumes the pressure is atmospheric (see Fig. 2). According to this approach, the flow front lies over the control volumes where they are adjacent to filled control volumes and are not completely filled. The flow front boundary condition is imposed at these control volumes. 4.1. Quasi-steady state approach

Fig. 1. (a) A typical eight-node hexahedral element in global and local coordinate systems and (b) a flux element divided into subcontrol volumes with one removed.

the element faces and those that are in the interior of the element. The latter group will form the surfaces of the control volume and here, are labeled A, B and C for an SCV, as shown in Fig. 1(b). The x, y and z locations of any point within the element, as well as the value of the pressure, can be evaluated through the usual use of shape functions. Thus we have xðt; s; rÞ ¼

8 X Ni x i

ð5Þ

i¼1

yðt; s; rÞ ¼

8 X N i yi i¼1

4.1.1. Numerical formulation To discretize the governing equation for fluid flow, the continuity equation over an arbitrary control volume must be integrated. After, integrating the continuity equation over a fixed control volume, divergence theorem is used to obtain: ðð 1 sð ½K :ðrp   gÞdS ¼ 0 ð9Þ  in which Darcy’s law is used to substitute velocity components. Differential surface vectors, dS, has three components dSx, dSy and dSz. The integral of Eq. (9) is exact for a control volume and it can be applied over any of them. Since element by element assembly makes best use of the fact that all geometric information is defined on an elemental basis, the integrals are divided into components evaluated over the elemental SCVs that form the control volume. Therefore, if Eq. (9) is applied for an SCV with three discrete SSs, shown in Fig. 1(b), the following equation is obtained: ðð

ð6Þ

ðð 1 1 ð ½K :ðrp  gÞ:dB A ð ½K :ðrp  gÞ:dA þ B   ðð 1 þ ð10Þ C ð ½K :ðrp  gÞ:d þ  ¼ 0 

866

A. Shojaei et al. / Composites Science and Technology 62 (2002) 861–879

Fig. 2. Schematic view of the filling process.

in which  denotes the flow rates on the other surfaces of an SCV. The pressure gradient term, rp; in the above equation is simply derived in local coordinates by taking derivatives of shape functions and substitution of Eqs. (5)–(7). The resultant equation in closed form is expressed as follows [39]: rp ¼ ½J 1 ½S fpg

2 3 gx 

vector ½W ¼ 4 gy 5 and ipA illustrates the value of gz the matrices in the ip. Using this procedure for subsurfaces B and C, and summing all terms in an appropriate manner, discretized equation for a SCV is expressed as follows:

ð11Þ 

The discretization of the governing equation leads to create several integrals which must be performed over the SCV surfaces. Here, an approximation is used to convert the continuous surface integrations to a discrete form by evaluating them at the defined point on SCV subsurfaces. This point which is called integration point (ip), is located on midpoint of any SCV subsurface, as shown in Fig. 1(b). In this approximation, the integrals are replaced by flux term evaluated at the ip times its subsurface area. On the other hand, to perform these integrations, the argument of the integrals will be required at the ip. This approximation which was widely used by Schneider et al. [40,41], can help to further simplify the discretized equation. On the basis of this approximation and using Eq. (11), we have the following relation for subsurface A: ðð 1 1 < Ai > ½K ipA A ð ½K :ðrp  gÞ:dA ¼  ipA ð12Þ 

1 J ipA ½S ipA f pg  ½W ipA where Ai denotes the components of the vector of the subsurface A in x, y, and z directions (< Ak > = < Ax Ay Az > ), [W] represents the components of the gravity

1 1 < Ak > ½K ipA ½J1 ipA ½S ipA þ < Bk > ipA ipB  1 ½K ipB ½J1 ipB ½S ipB þ < Ck > ½K ipC ½J1 ipC ½S ipc fpg ipC  1 1 þ¼ < Ak > ½K ipA ½W ipA þ < Bk > ipA ipB  1 ½K ipB ½W ipB þ < Ck > ½K ipC ½W ipC ð13Þ ipC

Since elemental assembly is taken place to make the global matrices, the following closed form relation can be expressed for an arbitrary element in solution domain: 8  X 1 1 < Ak > ½K ipA ½J1 ipA ½S ipA þ < Bk > ipA ipB i¼1

½K ipB ½J1 ipB ½S ipB þ

 1 < Ck > ½K ipC ½J1 ipC ½S ipc fpg ipC i

8 8  X X 1 1 i ¼ < Ak > ½K ipA ½W ipA þ < Bk >   ip ipb A i¼1 i¼1  1 ½K ipB ½W ipB þ < Ck > ½K ipC ½W ipC ð14Þ ipC i

A. Shojaei et al. / Composites Science and Technology 62 (2002) 861–879

P8

where the i¼1 i denotes the additional flux terms coming from each SCV in an element. Since this term is vanished for a CV, it will be ignored during making the global matrix. Eq. (14) is the working equation to solve the nodal pressure and can be applied to all elements within the solution domain yielding a set of linear algebraic equations by applying an appropriate assembly procedure. In matrix notation, it can be expressed as: Ap ¼ b

ð15Þ

A has the order of n, number of nodes in cavity, and each row denotes the mass balance for the corresponding control volume. All components of the empty and partially filled control volumes are zero. Therefore, A has several zero rows at the early stage of filling process. During filling algorithm, for a control volume that the fill fraction reaches to the value 1, a single row should be added to the stiffness matrix corresponding to the node number. When the material properties, fiber permeability and resin viscosity assumed to be constant, the components of the stiffness matrix will be constants during mold filling. In this case, the components of the stiffness matrix can be calculated once for all of the nodes in the cavity, rather than calculating and assembling at each time step. Therefore, stiffness matrix A is updated at each time step by adding those rows of pre-calculated stiffness matrix which correspond to filled control volumes. 4.1.2. Numerical scheme After the assembling of the stiffness matrix and imposing the boundary conditions, the system of equations is solved for unknown nodal pressures using an appropriate solver. Here, conjugate gradient method [42] is used to solve the nodal pressures in solution domain. It speeds up the calculation by one order of magnitude as compared with traditional solvers. After the pressure distribution is calculated for a given step, the velocity field including flow rate between CV surfaces can be computed using the Darcy’s law. When the permeability is a function of fluid velocity, an iterative solution algorithm is needed. The steps involved in the present computational scheme for simulation of three-dimensional mold filling process in isothermal condition are discussed below. Step 1: at the beginning of the simulation, the control volumes are constructed from the generated finite elements around each node, and then pore volumes of the control volumes are calculated. For the first time step, it is assumed that control volumes associated with the inlet nodes are fully filled with resin. It means that the fill fraction for these nodes are set to 1 and are considered as solution domain.

867

Step 2: the fiber permeability is calculated using a set of initial guesses for the velocity components when it is a function of velocity. Initial guess for the first time step is zero velocity. During the solution procedure, permeability values at the previous time step are used as initial guesses for the next time steps. With these values, the stiffness matrix is assembled and the nodal pressures are calculated. The permeability values are then re-calculated at ips using the newly computed velocity field and again the nodal pressures are calculated. This iterative procedure is continued until the nodal pressure values converge. The convergence criterion used in the present work is given by: 0

N P

11=2

jþ1 j 2 B jpi  pi j C Bi¼1 C B N C @ P jþ1 2 A jpi j i¼1

4

ð16Þ

where the convergence criterion, ; used in this study is 104 which is applied for total nodes in solution domain (N). In some cases, applying under-relaxation on the pressure can speed up the convergence. Step 3: once the nodal pressures are determined, the velocity field is calculated for completely filled control volumes connected to the flow front control volumes to compute the net flow rate entering into those control volumes. In order to do this, fluid velocity at ip of each connecting area is initially computed by Eq. (3) and then multiplied with that area to compute the volumetric flow rate between control volumes. Step 4: the fill fractions are updated by calculating the resin volume entering into control volumes using the following equation: Qitþ
t ¼ Qit þ

M X

qmi
t

ð17Þ

m¼1

where qmi is the flow rate between adjacent control volumes from m to i and M is the number of control volumes connected to control volume i. After computing the time step, the position of the flow front is obtained by interpolating the fill fraction field for the contour corresponding to a specified value, namely f=0.5. Step 5: the steps 2–4 are repeated for the new solution domain until the mold is completely filled. Darcy’s law is a steady state equation, while the filling process is a transient problem. This difficulty in solution algorithm can be removed by considering a quasi-steady state process in which one assumes a series of sequential steady conditions during filling process. For this purpose,

868

A. Shojaei et al. / Composites Science and Technology 62 (2002) 861–879

the time increment must be chosen in such a way that one control volume is completely filled. In practice, more than one control volume may be filled in a single time step. This usually occurs because of regular mesh and tolerance to define the fill fraction making a filled control volume, namely 0.99%. This restriction on time increment ensures the stability of the quasi-steady approximation.

an appropriate solution procedure. A set of nonlinear equations defined by Eq. (21) are solved having satisfied the following condintions:

4.2. Nodal partial saturation approach

4.2.2. Numerical scheme The present algorithm is proposed not only to compute the flow field and flow front but also to determine the pressure field during mold filling. This enables us to include the velocity dependent parameters, here velocitydependent permeability, in the algorithm and makes the algorithm to be convenient to incorporate the heat and mass transfer models in non-isothermal mold filling simulation. The steps involved in the present scheme are presented below.

4.2.1. Numerical formulation To develop a working equation for this scheme, Eq. (1) is integrated over an arbitrary control volume so that ð ððð @  !dV ¼  r:u dV ð18Þ @t V v in which subscript V denotes integration over a CV. Substituting Darcy’s law into the above equation and applying Divergence theorem on the right hand side of Eq. (18), leads to following expression: ð ðð @ 1  ð19Þ !dV ¼  ð ½K :ðrp   gÞ:dS @t V  s This is similar to Eq. (9) with an additional transient term on the left. On the basis of fill fraction concept, the integral at the left hand side of above equation for a control volume V can be replaced by product Vf, hence, the resulting equation is ðð @f 1 V  ð20Þ ð ½K :ðrp   gÞ:dS @t  s

0

b0i ¼


t Ai Vi 


t bi Vi 

For all nodes; 0 4 f 4 1

Step 1: repeat step 1 of the previous scheme. Step 2: if the permeability is a function of fluid velocity, it is calculated by a set of initial guesses as mentioned in step 2 of the previous scheme. Step 3: the flow field and the location of flow front are computed using an iterative algorithm [36,37]. To this end, after selecting a time step size, fill fraction field, f, is determined by an iterative procedure using Eq. (21) with zero gravity head: A0 p jþ1 ¼ fo  fjc

ð24Þ

ð22aÞ ð22bÞ

fcijþ1 ¼ max½0; minð1; fcij¼1 Þ

ð21Þ

The components of A0 and b0 are related to matrices A and b in Eq. (15) through the following expressions. Ai ¼

ð23Þ

To do this, fo is set to fill fraction field at the end of previous time step and the current value of f at each iteration, fc, is initialized with this value of fo at the beginning of iteration. After imposing the boundary conditions, the resultant set of equations are solved by appropriate solver ( currently conjugate gradient technique is employed as mentioned earlier) to find the nodal pressure. Since the current iterative fill fraction may have been incorrect, the fill fraction field is updated using unmodified Eq. (21).The updated value of fc is applied at every node, followed by the under/over shoot correction:

Further numerical treatment of right hand side of Eq. (20) similar to the procedure used for quasi-steady state formulation, Eqs. (9)–(14), and using an implicit time integration, gives A0 p ¼ b 0 þ fc  fo

For all nodes; p 5 pvac For all nodes where 0 4 f<1; p ¼ pvac

in which Ai and bi represent the components of ith row 0 in matrices A and b, and Ai and b0i are their correspondence in the present formulation. Eq. (21) is a working equation to determine the fill fraction field and position of the flow front by employing

ð25Þ

where subscript i denotes ith node. This correction is used to account for control volumes that complete filling in the given time step. This procedure continues until changes in the sum of the fill fractions between iterations falls below a small value; namely 106.

A. Shojaei et al. / Composites Science and Technology 62 (2002) 861–879

Step 4: the step 3 is only used to determine the fill fraction field for a selected time step size. After calculation of the fill fraction field, the pressure field is computed by imposing the appropriate boundary condition and solving this relation: 0

Ap¼b

0

ð26Þ

869

For the case of constant pressure, an specified pressure is set to inlet nodes: p ¼ pinj

ð29Þ

The pressure is set to zero gage pressure or a specified vacuum level at the flow front: p ¼ pvac

ð30Þ

This gives the pressure field for filled region at the end of each time step. 5. Validation of numerical simulation Step 5: the fiber permeability is updated by newly computed pressure field as described in step 2 of quasi-steady state scheme.

The simulation results of the developed computer codes, RTMS and RAPFIL, are compared with analytical solutions of some simple geometries for which the exact

Step 6: steps 2–5 continue until the pressure fields converge. Then flow front is obtained by interpolating the fill fraction field for f=0.5. Step 7: Steps 2–6 are repeated until the mold is completely filled. Both numerical schemes, quasi-steady state and nodal partial saturation methods, are implemented in FORTRAN 77 on a personal computer, Pentium 650 MHz. The computer code based on quasi-steady state approach is given the name RTMS and the other one on the basis of nodal partial saturation is called RAPFIL. The flow charts of the numerical schemes developed for both codes are illustrated in Figs. 3 and 4. 4.3. Boundary conditions In order to compute the mathematical model, the governing equation must be solved in conjunction with an appropriate set of boundary conditions. Fig. 5 illustrates the boundaries used in the present study. The boundary conditions are divided into three regions: mold wall, inlet gate and flow front. At the mold walls there is no flow in the normal direction and the condition is expressed as: @p !

¼ 0

ð27Þ

@n

There are two types of resin injection systems at the inlet gate: constant flow rate and constant pressure. At the constant flow rate, a prescribed flow rate is set to the control volumes enclosing the inlet nodes and the condition in those control volumes is as: ðð Sð

1 ½K :ðrp  gÞ : dS ¼ qinj 

ð28Þ Fig. 3. Flow chart of the numerical scheme for RTMS.

870

A. Shojaei et al. / Composites Science and Technology 62 (2002) 861–879

solution are known, and by checking convergence with mesh refinement to evaluate the accuracy of the numerical simulation in terms of mold filling time, pressure and flow front location. Also, the effect of time step size on the simulations resulted by the RAPFIL is studied. In some cases that there are no closed form solution for

pressure distribution, only the mass balance evaluation is performed for constant flow rate boundary condition. The material properties and process parameters used here for analytical solutions are listed in Tables 1 and 2. Table 1 shows three types of fiber preform, where kx, ky, kz denote their permeabilities in principal directions. 5.1. One-dimensional flow The Geometric details of this test case is shown in Fig. 6. The resin is injected at one surface of a rectangular mold and air is removed from the other surface. This can be regarded as a one-dimensional flow. The analytical solutions for the location of flow front, pressure distribution at time tff and mold filling time for both cases of constant flow rate and constant pressure at the inlet gate are given in Table 3 and the details of derivations have been given by Cai [43]. The expression of mold filling time, tf ¼  Vc =qinj is a general form for mold filling time of any geometry of mold cavity for resin injection under constant flow rate and is applicable to one-, twoand three-dimensional cases. Table 1 Permeability of three types of fiber preform (=0.5) used in the simulation Fiber

kx (m2)

ky (m2)

kz (m2)

1 2 3

109 109 109+107 v

109 2.51010 109+107 v

109 1010 109

Table 2 Process parameters used in the flow simulation studies Parameter

Value

3

5107 150 0.4 1150

qinj (m /s) pinj (kPa)  (Pa. s)  (kg/m3)

Fig. 4. Flow chart of the numerical scheme for RAPFIL. Table 3 Analytical expressions for one-dimensional flow

Fig. 5. Schematic view of the boundaries.

Parameter

Constant flow rate

Mold filling time

tf ¼

Pressure

pinj ðtff Þ ¼

Flow front position

Xðtff Þ ¼

Vc hWL ¼ qinj qinj qinj Xðtff Þ hWk

qinj tff hW

Constant pressure tf ¼

L 2 2kpinj

pðxÞ ¼ pinj ð1 

Xðtff Þ ¼

x Þ Xðtff Þ

rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2kpinj tff 

A. Shojaei et al. / Composites Science and Technology 62 (2002) 861–879

871

Fig. 6. Geometric details of one-dimensional test problem: (a) injection at horizontal position and (b) injection at vertical position.

In the case of vertical injection [see Fig. 6(b)], inlet pressure under constant injection flow rate is modified to the following relation (see Appendix): qinj  þ gÞZðtff Þ pinj ðtff Þ ¼ ð hWk

pinj ¼ ð31Þ

The location of flow front for filling of a vertical cavity under a constant injection pressure is as follows: Zðtff Þ 
Zðtff Þ  2 lnð1  Þ tff ¼



qinj Xðtff Þ þ pvac kS

ð35Þ

Flow front location under constant injection pressure 1 2 b a X ðtff Þ  pinj Xðtff Þ ¼ pinj tff 2  

ð36Þ

5.2. Two-dimensional flow in a square mold

ð33Þ

and for flow front location under constant injection pressure, we have 1 2 k X ðtff Þ ¼ ðpinj  pvac Þtff 2 

qinj Xðtff Þ aS þ bqinj

ð32Þ

where
¼ gk=  ¼ kpinj =: The detailed derivation of Eq. (32) can be found elsewhere [37]. When a vacuum is applied at the flow front, the following expressions are obtained (see Appendix): For inlet pressure under constant flow rate pinj ðtff Þ ¼

Inlet pressure under constant injection flow rate

ð34Þ

For the case of linear-velocity dependency for preform permeability, the following relations derived in the Appendix can be used for:

In this test case, a square mold of size 0.2 m0.2 m0.005 m with an injection gate at the center is used. The fluid is injected at a constant flow rate so that the analytical mold filling time can be simply calculated by dividing pore volume of the mold cavity to the injection flow rate. The fluid flows in two directions in Cartesian coordinates, while, the flow front has a circular shape for an isotropic preform and an elliptical shape for an orthotropic preform before the flow front reaches to the mold walls. The ratio of major and minor axes in an elliptical flow front is equal to the square root of the ratio of the principal permeabilities. 5.3. Results For a one-dimensional flow problem, one element in the thickness and two elements in the width of the mold are used while the number of elements with equal sizes

872

A. Shojaei et al. / Composites Science and Technology 62 (2002) 861–879

in the length of the mold is varied. The fiber used in this study is type 1 in Table 1. For both cases of the inlet boundary conditions, very excellent results were observed for this filling pattern even in a very coarse mesh when the permeability was a constant. The comparison between the numerical and analytical results for

both developed numerical codes are shown in Figs 7–9. In these examples, different time step sizes are used for the RAPFIL. The simulated flow front and inlet pressure are compared with the analytical solutions in Fig. 7 for a constant injection flow rate. Fig. 8 illustrates the analytical and numerical results for the flow front location

Fig. 7. Simulations resulted by the RTMS and the RAPFIL for onedimensional flow with a constant injection flow rate; (a) flow front location; ( b) inlet pressure and (c) inlet pressure with gravity.

Fig. 8. Simulations resulted by the RTMS and the RAPFIL for onedimensional flow with a constant injection pressure; (a) flow front location; (b) flow front location with gravity and (c) pressure distribution when the flow front reaches to X=0.24 m.

A. Shojaei et al. / Composites Science and Technology 62 (2002) 861–879

and pressure distribution during mold filling under the constant injection pressure. An excellent agreement is observed between the numerical and analytical results for this test problem. Also, the results are excellent for large time step size for the RAPFIL with more iterations within the convergence tolerance. As it is seen in the figures, the effect of gravity are not significant on the results. In fact its effect is important when the inlet pressure or inlet flow rate is low. In the next case, a linear function of fluid velocity for fiber preform permeability is employed to verify the accuracy of the iterative procedure used for simulation of the mold filling (see type 3 in Table 1). It was observed that the numerical filling time tended to the analytical value as the mesh was refined so that for a 30 elements in the length direction the error fell below 1 percent for the both computer codes. Fig. 9 compares

Fig. 9. Simulations resulted by the RTMS and the RAPFIL for onedimensional flow with a velocity-dependent permeability, Number of elements in length direction=30; (a) inlet pressure with a constant injection flow rate and (b) flow front location with a constant injection pressure.

873

the numerical results and theoretical solutions for onedimensional mold filling under the constant flow rate and pressure inlet boundary conditions. From the figure it is clear that the numerical predictions agree well with the analytical data. A small impregnation length which is observed in the figures at the onset of the filling process, shows the inlet control volumes that are assumed as fully filled at the beginning of the simulation. For the central injection in a square mold, the simulated flow front has the expected elliptical shape (see Fig. 10). The ratio of the minor to major axes of the flow front was compared with the analytical value for different mesh numbers in planer direction when the flow front along the major axis reaches near the mold

Fig. 10. Simulated flow front locations at two filling stage resulted by; (a) the RTMS and (b) the RAPFIL, solid line for
t=20 s and dotted line for
t=2 s.

874

A. Shojaei et al. / Composites Science and Technology 62 (2002) 861–879

walls. It was observed that error in flow front location is reduced by increasing planer mesh number. The simulated flow front location by the both codes for a mesh number of 3030 in planer direction are shown in Fig. 10. Fiber preform type 2 in Table 1 is used in this study. The error in the flow front location by the RTMS falls below 0.5% for the mesh number described above. As shown in Fig. 10, the results by the RAPFIL for two time step sizes, 2 and 20 s, are exactly similar and in this case the error falls below 0.5% too. The lager time step size takes slightly shorter CPU time, but the small one gives more information regarding the flow progression during mold filling. Therefore, time step sizes in the RAPFIL do not affect the prediction of flow front

Fig. 11. Geometry of the cubic mold with a line gate.

Fig. 12. Simulated flow front locations with constant permeability at equal time intervals of 8 s, by; (a) the RTMS and (b) the RAPFIL.

Fig. 13. Simulated flow front locations with velocity-dependent permeability at equal time intervals of 8 s, by; (a) the RTMS and (b) the RAPFIL.

A. Shojaei et al. / Composites Science and Technology 62 (2002) 861–879

875

locations. The analytical mold filling time for this test problem is 200 s. The simulated mold filling time by the RTMS is 193.5 s while the mold is filled by the RAPFIL after 10 steps for
t=20 s and 100 steps for
t=2 s which the results agree well with analytical solution.

Fig. 14. Simulated flow front locations with anisotropic preform at equal time intervals of 100 sec, by; (a) the RTMS and (b) the RAPFIL,
t=5 s.

Fig. 16. Simulated flow front locations at equal time intervals of 6 s, resulted by; (a) the RTMS and (b) the RAPFIL,
t=1 s.

Fig. 15. Geometric details of complex shaped mold.

876

A. Shojaei et al. / Composites Science and Technology 62 (2002) 861–879

6. Applicability of the developed computer codes Three test mold geometries are used here to demonstrate the applications of the developed computer codes, RTMS and RAPFIL, for simulation of filling process in

RTM under isothermal condition. The material properties and injection parameters are given in Table 2. The first case is a simple cubic mold with dimensions of 0.15 m0.15 m0.15 m. The fluid is injected through a line gate located along a side of the bottom surface, as

Fig. 17. Simulated pressure distribution near the end of mold filling resulted by; (a) the RTMS and (b) the RAPFIL,
t=1 s.

A. Shojaei et al. / Composites Science and Technology 62 (2002) 861–879

shown in Fig. 11. In such a case, the fluid progresses in x and z directions. The mold is filled under a constant injection pressure with two types of permeability: the first is an isotropic preform (type 1 in Table 1) and the second is a fiber preform that the permeabilities in x and y directions are functions of the fluid velocity (type 3 in Table 1). The total number of nodes and elements are 1331 and 1000, respectively, and the time step size used for the RAPFIL is 2 s. Fig. 12 shows the flow front locations obtained from the both codes throughout the mold filling with 8 sec time intervals for the isotropic preform. In this case, the fluid progresses in an equal position at x and z directions due to the same permeabilities in these directions, kx=kz. As expected, the flow front merges at front side of top surface of the mold as illustrated in Fig. 12. The simulated filling time is 76.66 s for RTMS and 39 time step with
t=2 s for the RAPFIL which shows close agreement between the results obtained by both codes. The flow fronts predicted by the RTMS and the RAPFIL for a velocity-dependent permeability are shown in Fig. 13 and again a close result is observed from both codes. The flow front progression in x direction is faster than the one in z direction due to the low permeability of the fiber preform in the z direction. This is because the permeability in x direction is a function of fluid velocity. Thus, when the fluid reaches to the end location of the bottom surface of the mold, the flow front has not reached to the top surface of the mold yet. In this case, predicted mold filling times by both codes are close ( 56.4 s for RTMS and 29 time steps with
t=2 s for RAPFIL). A hollow rectangular shaped mold, shown in Fig. 14, is the next considered geometry. The outer dimensions of the mold are 0.16 m0.12 m0.04 m and the inner ones are 0.08 m0.08 m0.04 m. Number of elements and nodes are 512 and 860, respectively. A point gate located at a corner of the bottom surface of the mold is used in the simulation. The fluid is injected at the constant injection pressure with a value illustrated in Table 2. The fiber preform used in this study is an anisotropic one and its permeabilities are kx: ky: kz=10:2.5:1. (see type 2 in Table 1 ). The predicted flow fronts are shown in Fig. 14. As shown in the figure, the fluid is divided into two sections, one into the front arms and the other into the back arms of the mold and then, the flow fronts merge at the other corner side of the mold. Also, at early stage of the mold filling, there is a lead-lag in the flow front between the top and bottom surfaces of the mold due to the low permeability of the preform in z direction. As it is seen in the figure, the flow front locations and mold filling time predicted by both codes are very close. The last case is the simulation of the filling process for a complicated three-dimensional geometry illustrated in Fig. 15. It has a square hole at the top. The geometry is discretized into 2076 nodes and 994 elements. The fluid

877

is injected simultaneously through two injection points located at both sides of the square hole under the constant injection flow rate. An isotropic fiber preform (type 1 in Table 1) is used in this simulation. The flow front positions calculated using both codes versus time are shown in Fig. 16, and a very close agreement is observed between the both results. At early stages of the mold filling, fluid progresses radially from the injection ports until it reaches to the square hole, then the flow fronts merge together under the lower sides of the hole. From the simulated flow front positions, the locations of weld lines and vents can be estimated which help one to prevent air entrapment inside the molded part. Fig. 17 illustrates the pressure fields resulted by both codes near the completion of mold filling and shows a close result from both codes. As shown in the figure, maximum pressure is near the inlet gates. The resulting approximation is 58.02 sec for RTMS and 60 sec for RAPFIL (60 time steps with
t=1 s) to fill the mold versus the analytical value obtained by dividing the total pore volume of the cavity to total flow rate, 59.3 s. It is worth noting that the CPU time for RTMS is more than 100 h on our computer system, Pentium 650 MHz, while it is less than 7 h for the RAPFIL representing the CPUeffectiveness of the RAPFIL. Although the CPU-effectiveness of the RAPFIL was observed in our previous test cases, but this result can be valuable for simulation of three-dimensional complicated geometries in which the number of nodes are usually high.

7. Conclusions In this paper, two computer codes, RTMS and RAPFIL, have been developed to simulate three-dimensional filling process as well as flow front locations and pressure field during mold filling. This information can help one to design of the tooling and processing conditions before a mold is built. An important feature of the developed codes is their flexibility to consider some filling problems such as velocity-dependent permeability, gravity contribution and vacuum effect. An excellent agreement is observed for both codes on comparison with analytical solutions. By using the codes, mold filling process of some geometries has been simulated to investigate their applicability, and a close agreement has been observed between the results of both codes. Also, the RAPFIL shows a considerable CPU time-effectiveness for simulation of a complicated three-dimensional cavity.

Appendix. Analytical solution for 1-D isothermal flow Let us observe the flow of fluid into a 1-D cavity as shown in Fig. 6. For this problem, Darcy0 s law and continuity equation are expressed as:

878

A. Shojaei et al. / Composites Science and Technology 62 (2002) 861–879

k dp u ¼  ð þ gÞ  dx

ðA1Þ

du ¼0 dx

ðA2Þ

Substituting Darcy0 s law into Eq. (A2) yields: d k dp ð ð þ gÞÞ ¼ 0 dx  dx

ðA10Þ

For the case of velocity dependent permeability with no vacuum at the flow front, pvac=0, Eq. (A10) becomes: ðA3Þ

This is the governing equation for 1-D flow, and the boundary conditions are  p ¼ pinj for constant injection pressure at x ¼ 0 q ¼ qinj for constant injection flow rate ðA4aÞ p ¼ pvac at x ¼ Xðtff Þ

  1 2 b X ðtff Þ  ðpinj  pvac ÞXðtff Þ 2  a ðpinj  pvac Þtff ¼ 

1 2 b a X ðtff Þ  pinj Xðtff Þ ¼ pinj tff 2  

ðA11Þ

If one assumes a constant permeability with a vacuum at the flow front, Eq. (A10) is simplified into: 1 2 k X ðtff Þ ¼ ðpinj  pvac Þtff 2 

ðA12Þ

ðA4bÞ

in which moving front, X(tff) is expressed as:

u¼

dX dt

Analytical expressions for 1-D flow under constant injection flow rate. ðA5Þ The solution of Eq. (A3) leads to following expression for pressure distribution:

We assume a linear function of fluid velocity for permeability as: pðx; tff Þ ¼ pinj ðtff Þ þ k ¼ a þ bu

ðA6Þ

The solution of Eq. (A3) is ðA7Þ

Substitution of this equation into Eq. (A1) without gravity contribution, using Eqs. (A5) and (A6), and further simplification leads to the following expression: ðXðtÞ 

b dXðtÞ a ðpinj  pvac ÞÞ ¼ ðpinj  pvac Þ  dt 

ðA13Þ

By substituting Eq. (A13) into Eq. (A1), we get the following relation for inlet pressure with no vaccum at the flow front

Analytical expressions for 1-D flow under constant injection pressure

pvac  pinj pðx; tff Þ ¼ pinj þ x Xðtff Þ

pvac  pinj ðtff Þ x Xðtff Þ

ðA8Þ

The boundary conditions for this differential equation are  XðtÞ ¼ 0 at t ¼ 0 ðA9Þ XðtÞ ¼ Xðtff Þ at t ¼ tff The solution of Eq. (A8) leads to the following relation for the flow front location at time tff:

qinj  þ gÞZðtff Þ pinj ðtff Þ ¼ ð kS

ðA14Þ

in which Z(tff) denotes the position of the flow front according to notations in Fig. 6. For horizontal mold filling with a vaccum at the flow front, inlet pressure is obtained by: pinj ðtff Þ ¼

qinj Xðtff Þ þ pvac kS

ðA15Þ

When the permeability is a function of fluid velocity, the following equation is obtained by substituiting Eq. (A13) into Darcy0 s law and considering Eq. (A6): pinj ¼

qinj Xðtff Þ aS þ bqinj

in which pvac has been set to zero.

ðA16Þ

A. Shojaei et al. / Composites Science and Technology 62 (2002) 861–879

References [1] Potter K. Resin transfer moulding. London: Chapman & Hall, 1997. [2] Advani SG, Bruschke MV. Resin transfer moulding flow phenomena in polymeric composites. In: Advani SG, editor. Flow and rheology in polymer composites manufacturing. Amsterdam BV: Elsevier Science, 1994 [Chapter 12]. [3] Owen MJ, Middleton V, Hutcheon KF. The development of resin transfer molding (RTM) for volume manufacture. Proc I Mech E Design in Composite Materials 1989:107–14. [4] Owen MJ, Rudd CD, Middleton V, Kendall KN, Revill ID. Resin transfer molding (RTM) for automotive components. Composite Material Technology ASME 1991;PD-Vol. 37:177–83. [5] Hayward JS, Harris B. Processing factors affecting the quality of resin transfer molding composites. Plastics and Rubber Processing and Applications 1989;11:191–8. [6] Rudd CD, Kendall KN. Towards a manufacturing technology for high volume production of composite components. Proc I Mech E Part B 1992;206:77–91. [7] Patel N, Rohatgi V, Lee LJ. Influence of processing and material variables on resin-fiber interface in liquid composite molding. Polymer Composites 1993;14:161–72. [8] Rudd CD, Long AC, Kendall KN, Mangin CGE. Liquid moulding technologies. Cambridge: Woodhead Publishing Limited, 1997. [9] Patel N, Lee LJ. Modeling of void formation and removal in liquid composite molding. Part I: wettability analysis. Polymer Composites 1996;17:96–103. [10] Chen Y, Davis HT, Macosko CW. Wetting of fiber mats for composites manufacturing: I. visualization experiments. AIChE Journal 1995;41:2261–73. [11] Bruschke MV, Advani SG. A finite element/control volume approach to mold filling in anisotropic porous media. Polymer Composites 1990;11:398–405. [12] Young WB, Rupel K, Han K, Lee LJ, Liou M. Analysis of resin injection molding in molds with preplaced fiber mats. II: numerical simulation and experiments of mold filling. Polymer Composites 1991;12:30–8. [13] Chan AW, Hwang ST. Modeling of the impregnation process during resin transfer molding. Polymer Engineering and Science 1991;31:1149–56. [14] Trochu F, Gauvin R, Gao DM. Numerical analysis of resin transfer molding process by the finite element method. Advances in Polymer Technology 1993;12:329–42. [15] Lee LJ, Young WB, Lin RJ. Mold filling and cure modeling of RTM and SRIM processes. Composites Structures 1994;27:109–20. [16] Bruschke MV, Advani SG. A numerical approach to model nonisothermal viscous flow through fibrous media with free surfaces. International Journal For Numerical Methods in Fluids 1994;19: 575–603. [17] Kang MK, Lee WI, Yoo JY, Cho SM. Simulation of mold filling process during resin transfer molding. Journal of Materials Processing & Manufacturing Science 1995;3:297–313. [18] Yoo YE, Lee WI. Numerical simulation of the resin transfer mold filling process using the boundary element method. Polymer Composites 1996;17:368–74. [19] Mal O, Couniot A, Dupret F. Nonisothermal simulation of the resin transfer molding process. Composites Part A 1998;29:189– 98. [20] Lin M, Hahn HT, Huh H. A finite element simulation of resin transfer molding based on partial nodal saturation and implicit time integration. Composites Part A 1998;29:541–50. [21] Lam YC, Joshi SC, Liu XL. Numerical simulation of the moldfilling process in resin-transfer molding. Composite Science and Technology 2000;60:845–55.

879

[22] Shojaei A, Ghaffarian SR, Karimian SMH. Modeling and simulation approaches in resin transfer molding process—a review. Polymer Composites 2002:[in press]. [23] Bruschke MV, Luce TL, Advani SG. Effective in-plane permeability of multi-layered RTM preforms. Proc. 7th Tech. Conf. of American Society for Composites. University of Delaware 1992: 103–12. [24] Young WB, Han K, Fong LH, Lee LJ, Liou M. Flow simulation in molds with preplaced fiber mats. Polymer Composites 1991;12: 391–403. [25] Young WB. Three-dimensional nonisothermal mold filling simulation in resin transfer molding. Polymer Composites 1994;15: 118–27. [26] Trochu F, Boudreault JF, Gao DM, Gauvin R. Three-dimensional flow simulation for the resin transfer molding process. Materials and Manufacturing Processes 1995;10:21–6. [27] Lim ST, Lee WI. An analysis of the three-dimensional resintransfer mold filling process. Composite Science and Technology 2000;60:961–74. [28] Friedrichs B, Guceri SI. A hybrid numerical technique to model 3-D flow fields in resin transfer molding processes. Polymer Engineering and Science 1995;35:1834–51. [29] Maier RS, Rohaly TF, Advani SG, Fickie KD. A fast numerical method for isothermal resin transfer mold filling. International Journal For Numerical Methods in Engineering 1996;39:1405–17. [30] Coulter JP, Guceri S. Resin impregnation during the manufacturing of composite materials subject to prescribed injection rate. Journal of Reinforced Plastics and Composites 1988;7:200– 19. [31] Gauvin R, Trochu F. Comparison between numerical and experimental results for mold filling in resin transfer molding. Plastics, Rubber and Composites processing and Applications 1993;19:151–7. [32] Um MK, Lee WI. A study on the mold filling process in resin transfer molding. Polymer Engineering and Science 1991;31:765– 71. [33] Chan AW, Hwang ST. Modeling nonisothermal impregnation of fibrous media with reactive polymer resin. Polymer Engineering and Science 1992;32:310–8. [34] Lin R, Lee LJ, Liou M. Non-isothermal mold filling and curing simulation in thin cavities with preplaced fiber mats. International Polymer Processing 1991;VI 4:356–69. [35] Baliga BR, Patankar SV. A new finite element formulation for convection-diffusion problems. Numerical Heat Transfer 1980;3: 393–409. [36] Voller VR, Peng S. An algorithm for analysis of polymer filling of molds. Polymer Engineering and Science 1995;35:1758–65. [37] Voller VR, Peng S, Chen YF. Numerical solution of transient, free surface problems in porous media. International Journal For Numerical Methods in Engineering 1996;39:2889–906. [38] Das BM. Advanced soil mechanics. McGraw-Hill, 1983. [39] Reddy JN. An introduction to finite element method. Singapore: McGraw-Hill Book Company, 1985. [40] Schnieder GE, Zedan M. Control-volume based finite-element formulation of the heat conduction equation. Spacecraft Thermal Control, Design, and Operation, Progress in Astronatics and Aeronautics 1983;86:305–27. [41] Schneider GE, Raw MJ. Control volume finite-element method for heat transfer and fluid flow using collocated variables—1. Computational procedure. Numerical Heat Transfer 1987;11: 363–90. [42] Press WH, Flannery BP, Teukolsky SA, Vetterling WT. Numerical recipes-the art of scientific computing. Cambrige: Cambrige University Press, 1986. [43] Cai Z. Simplified mold filling simulation in resin transfer molding. Journal of Composites Materials 1992;26:2606–29.