Vol.16 No. 4

CHIN ESE J OL RN AL OF AERON A L T ICS

November 2003

Non-isothermal Mold Filling and Curing Simulation for Resin Transfer Molding

1 1 2 1 .3

CHEN Ren-iang , GUI Bing ' , LI Ming-cheng , LIANG Zhi-yong (1. College of Aeronautics & Astronautics, Nanjing University of Aeronautics & Astronautics,

Nanjing 210016, China) (2. Department of Mathmatics, Nanjing Forest University, Nanjing 210018, China) (3. College of Material & Engineering, Bey ing University of Aeronautics & A stronautics,

Bej ing 100083, China)

Abstract: A numerical model of 2. 5D non-isothermal resin transfer molding simulation is developed for thin part based on the control volume/ finite element method. The non-uniform temperature distribution and the heat generation during the filling stage are modeled with the lump«! temperature system and the species balance. Numerical algorithm of the simulation are studied. The molding simulation for a part is performed to show the ef fectiveness of simulating filling time, temperature distribu tion and cu ring degree.

Key words: resin transfer molding; finite-element analysis; mathematical model; numerical simulat ion

^№),2003, 16(4): 247- 252.

««miia» iiae * w ae, m%r ssewii

a ~ fisffts mm i w mmrn , * m b w ©*«a i^a*« »fit «ai^ag^ttiia^iasfftfittiiies»

1000-9361(2003)04-0247-06 A

Resin Transfer Molding (RTM) has become an attractive method in which liquid react ant cures in a mold with fabric reinforcement to produce a composite part- Interest in these technologies has grown dramatically, spurred by the military's need for affordability and high performance and the auto industry's emphasis on higher productivity and low er-volume m arket niches.

There are so many technical issues concerned with RTM processes'13]. Because of the anisotropic nature of the reinforcement fibers, the flow pattern in the mold can be very complicated. Cure of resins utilizes heat to cause a crosslinking react ion .

Received date: 2003-01-09; Revision received date: 2003-07-09 Artrle L RL: http: //www. lik^Kne^ cn/cja/2°°3/04/0247/

The rate of the reaction is dependent on the temperature and extent of resin conversion. Heat generation and transfer can result in a non-uniform temperature distribution in the mold.

Owing to the interaction of different variables and the complexity of the process, mold and process design are currently conducted using trial and error. Experimental methods for optimization of mold and process design may be too time consuming and economically prohibitive. A compute numerical simulation model capable of predicting resin flow behavior, pressure gradient inside the mold, and curing process is needed to optimize injection

Wishing House. All rights reserved, http://www.cnki.net

CHEN Ren-liang, GLIBing, Li Ming-cheng, LIANG Zhi-yong

port positions, injection pressure, filling time, temperature, and curing degree-

A number of studies have developed non-isothermal flow models with a heat transfer model'4-9]. These models are coupled by the fluid viscosity, which is the function of both temperature and cure conversion- However, most of them are concentrated on two-dimensional by simplifying the heat transfer form. In the actual case, the thermal effects are three-dimensional and require three-dimensional analysis. There are some studies on the non-isothermal RTM mold filling and curing processing with true three-dimensional analysis' 10 13] . Unfortunately, the huge computation time limits its application- Moreover, three-dimensional model makes the numerical algorithm more complicated.

In order to improve the two-dimensional model and avoid the huge computation time and numerical algorithm difficulty in the three-dimensional model, a 2. 5D simulation RTM non-isothermal model is proposed for thin part in this paper. The meaning of 2. 5D is that the flow is 2D while the heat transfer is 3D. In this 2. 5D model, the flow is determined by Darcy's law and the heat transfer and resin reaction are developed based on the lumped temperature system ( i. e. local thermal equilibrium between resin and fiber). The control volume/finite element method (CV/FEM) is used to solve the distributions of pressure, temperature and curing degree during the RTM processing.

1 Simulation Scheme

1. 1 Governing equations

In order to simulate the mold filling process of RTM, several assumptions must be made to simplify the problem. In general, the fiber reinforcement in the mold cavity is assumed to be rigid during the molding process. There is no exchange of mass between the solid and the fluid. The fluid is considered as incompressible and its viscosity is independent of shear rate ( i. e. , the fluid is Newtonian) . Inertia effect is neglected because of the low resin flow. Surface tension is negligible as com-

pared with the dominant viscous force.

Under these assumptions, Darcy's law for flow through porous media can be used as the momentum equation.

where v is the velocity vector, VP is the pressure gradient, V is the viscosity, and K is the permeability tensor. For the part with small thickness, the flow can be modeled as two-dimensional flow

dx dp l^vj

where vx, vy are the flow velocity components and Kij (i, j = x, y) is the component of the permeability tensor defined in the two-dimensional Cartesian coordinate system.

For an incompressible fluid, the continuity equation is

dvx cV.v

2x + cv = 0

Eq. (2) and ( 3) are combined and integrated over a control volume to have

[ nx nv ]

Kxv Kvv

►ds = 0 (4)

where nx, nv are the components of the vector normal to the surface of the integrated volume.

In RTM process, if the temperature is assumed to be in local equilibrium at any specific time due to its slow filling speed (i. e. , the temperature is t he same for both resin fluid and fiber mat at each local position), the equations of energy and species balance can be expressed as

PCp ^ + Pr Cpr-

V T = k V2 T + PAHm°(5)

Ox - „ ° ft + x= m

where T is the temperature, Wis the resin conversion, m is the mass generation rate, 'Pis the fiber mat porosity and AH is the volumetric heat generation per unit volume. The lumped thermal param-

et ers are defined as 3ublishing He p = (pf /_) / (ptw r + pM f) tp://www.cnki.net

November 2003

Non-isothermal Mold Filling and Curing Simulation _for Resin T ransfer Molding_

Cp = cprw r + cp w f k = ( kr kf) / ( kr w f+ kfw r)

w r= A

w f = 1 -

where P is the density, cpis the specific heat coefficient and k is the thermal conductivity- The subscripts r andf denote resin and fiber respectively. 1- 2 Boundary conditions

The boundary conditions for equation 4 through 6 are as follows At the flow front

= (1- Ac„fVn( Tfo - T ) D

( F 0) = Fm

where F is the resin volume fraction of a control volume. At the inlet gate

( a) For the constant pressure mold filling

P = Po, T = T„„in, and (X = 0 ( b) For the constant flow rate mold filling

■ 1 " n 1

JJ ^ [ n ]

■Kyx

Kxv KyvJ

ds = - Q

T = T res At mold boundaries

l^vj 0= 0

% = o T = T -

where Q is the volume injecting rate at the inlet gate, Tfo is the initial fiber temperature that is the same as the mold wall temperature Twall and Tasinis the initial resin temperature. At the flow front, the parameter F represents the status of each control volume in the flow domain. F is zero for an empty control volume, while it is 1 for an entirely filled control volume. For a partially filled control volume, F is equal to the volume fraction occupied by the resin. Pressure is assumed vanished at the flow front where the temperature and conversion follow the energy and mass balance of a control

mass. For the case of specified injecting flow rate,

© 1094-3010 China Academic Journal Elec__

a specified flow rate is assigned to the control vol-

umes enclosing the inlet nodes. For the case of the specified injecting pressure, a specified pressure is assigned to the inlet nodes. For the outlets, the pressure depends on the filling condition. If the vacuum is used during mold filling, the outlet pressure is equal to the vacuum pressure; otherwise, the outlet pressure is zero. 1. 3 Computation schemes

Eqs. ( 4) -( 6) are the governing formulations of non-isothermal mold processes. The solution of them requires the use of a numerical method to find the distribution of pressure, temperature and resin conversion within the filled domain. CV/ FEM has many advantages in solving the moving boundary problems'5'8'91.

In the CV/FEM, the entire calculation domain is first divided into a number of elements. The control volume is formed based on the nodes of the adjacent elements. Each element is divided into several sub-areas by the lines connecting the centroid of the element to the midpoint/center of each side/ surface. A control volume is composed of several sub-area/volumes, which have a common node at the center of the control volume. Fig. 1 is the typical control volume distribution.

(a) T rian gl e (b) T rian glean d q uar dran g le

Fig. 1 Typical control volume distribution For 2.5 D non-sothermal case, the heat transfer must be solved in three dimensions. Thus the cavity should be discretised through the thickness, which involves sub-dividing each control volume defined for the flow analysis into a number of layers. T he energy balance is then carried out within each of these sub-domains. The layers through the thickness can be divided with collocate method'14].

Fig. 2 shows the diagram of computation schemes for 2. 5D non-isothermal RT M simulation based on CV/FEM. At the beginning of mold fill-iirn^ahon pr°gram ass'j¡nes]-that thecon-

CHEN Ren-iang, GL 1 Bing, Li Ming-cheng, LIANG Zhi-yong

Fig. 2 The diagram of computation schemes trol volumes enclosing the inlet nodes are filled with fluid. Then, with the boundary conditions mentioned earlier and the assumed viscosity, the pressure calculation is conducted to determine the pressures of nodes that are currently filled with fluid. After the pressure field is determined, the velocity field can be evaluated according to Darcy's law. Then, the temperature and conversion can be solved by the same way as the pressure based on the current velocity field. The resulting temperatures and conversions are used to compute the update viscosity. New values of viscosity are used to calculate the pressure field again, which leads to the same calculation procedure.

The Courant Number in the above numerical scheme is used in the improvement of convergence and accuracy of the solution. The convergence problem comes from the convection terms in equations ( 5) and ( 6) . In order to reduce the divergence, the convection terms of the two equations are evaluated using the upwind scheme. The upwind scheme adds a numerical diffusion term to the system and reduces the instability. However, it increases the truncation errors. The truncation errors may result in a m eaningless solution even though the scheme is stable. They are sensitive to

the flow velocity, mesh size and time step if the

n • . • .mi . . . ■ c

flow is convection dominant. A sub-time step is

used in the numerical model to improve the computational accuracy of temperature and conversion with reasonable Courant Number'2]. This sub-time technique can ensure the solutions to be as accurate as possible, although the truncation error is still very difficult to control because of the varying mesh sizes in the model.

2 Case Study and Discussion

A part in reference [9] is taken to demonstrate the use of the method described above. Figure 3 shows the finite elements generated by PA-TRAN software at the mid-plane of the part. The geometry of the part includes a circular cut-out with resin injected through a 20 mm diameter inlet port. Table 1 and 2 list the material data and molding parameters for simulation.

Fig. 3 Finite element mesh at the mid-plane of the p art

Tabl e 1 Material data

Re sin Fiber

PI ( kg• m- 3) 1120 2540

Cp/(J'(kg -K)-1) 1400 840

k/( W '(m'K)- 1) 0. 13 0. 1/0. 01

A H/(J'm-3) 2. 5768X 108

<*> 0. 7

K/m2 1.0X10-8

Table 2 Molding parameters

Injection rate ( constant) / ( m3* s ') 6. 5X 10- 5

Cavity thickness/ mm 3

Resin inlet temperat ure/K 310

In itial perform temperature/ K 350

Mold surface temperature/K 370

The reaction rate is f («, T) = 300exp

(1- o):

The relationship among resin viscosity, temperature and degree of conversion used is

1. 45 X 104

V = 2. 56 x 10- 21 exp Publishing House. All rights reserved. http:^/www.cnki.net

November 2003

Non-isothermal Mold Filling and Curing Simulation fo r Resin T ransfer M olding

0. 695

1. 502

0.487- a

Fig. 4(a) shows the predicted flow pattern. In this figure, the initial radial flow quickly becomes rectilinear until the cut-out is reached. Once the flow passes the cut-out the tw o flow fronts merge, after which the flow becomes rectilinear again. The total fill time for this case is approximately 3.8 seconds. Fig. 5(b) is the flow patter in reference [9]. It can be found that the predicted flow pattern is in a good agreement with those in reference [ 9] . The difference in filling time comes from the fact that the model in reference [ 9] does not consider the in-plane heat conduction which affects the resin viscosity.

(b) Temperature field in reference [9] Fig. 5 Temperature field

(a) Predicted conversion field

(a) Predicted flow patta'n and filling time

(b) Flow patta'n and filling time from reference [9] Fig. 4 Flow pattern and filling time Fig. 5(a) and 6(a) are the distributions of predicted temperature and cure convection for the

(b) Conversion field in reference [ 9] Fig. 6 Flow pattern and filling time

layer located between the mold surface and the mid-plane, For this particular case, the maximum degree of conversion is 0. 0752, occurring at the location furthest from the injection gate (where the resin has the greatest residence time)- mid-plane-For this particular case, the maximum This value is significantly lower than the resin gel point of 0- 487, which indicates that the heat transfer during the filling stage is caused mainly by the convection- The temperature at each location decreases after the cold resin touches the hot fiber- After the resin moves forward, colder resin arrives which resultes in low er temperature, but still higher than © 1994-(a) Predict«! temperature field that at the injection port, at this location- As a re-

CHEN Ren-liang, GL 1 Bing, Li Ming-cheng, L1ANG Zhi-yong

sult, the temperature shows a gradual increase as the filling time increases in the whole filling domain.

Fig. 5(b) and 6(b) are the temperature and conversion field at the same layer made by reference [ 9]. There is a good agreement between them. The difference between them also comes from the effect of in-plane heat conduction.

3 Conclusion

A 2.5 D non-isothermal RT M simulation model is developed for the thin part based on the CV/FEM. The lumped temperature system and the species balance are used in modeling the non-uniform temperature distribution and the heat generation during filling stage. The convergence and accuracy of numerical schemes are studied and simulation for a part is performed. Compared with 2D model, 2. 5D model do affect the flow pattern, filling time, distributions of temperature and cure conversion in RT M process.

To verify the process model, experimental work is needed. Results from the experimental study and the comparison with the numerical simulation are the future study in this area.

References

| 1| Tucker C L. Heat transfer and reaction issues in liquid composite molding | J| . Polymer Composites J, 1996, 17 ( 1) : 60- 72.

1 2| GutowskiT G. Advanced composite manufacturing 1 M !.

Portland: John Wiley & Sons Inc, 1997. 395- 399. |31 AntonelliD, Farina A. Resin transfer molding: mathematical modeling and numerical simulations | J | . Composites Part

A J, 1999, 30( 12): 1367- 1385. |4| Lin, RJ, LeeLJ, Liou M L. N on-isothermal mold tilling and curing simulation in thin cavities with preplaced fiber mats | J | . international Polymer Process J, 1991, 6(4): 356- 369.

| 5 | Lin R J, Lee L J , Liou M L. Mold filling and curing analysis in liquid composite molding | J| . Polymer Composites J,

1993, 14( 1): 71- 81. | 6| Trochu F, Gauvin R, Gao D M. Numerical analysis of the

resin transfer molding process by the finite element method | J| . Advances in Polymer ' Technology J, 1993, 12(4): 329 - 342.

| 7| Bruschke M V, Advani S G. A numerical approach to model non-isothermal viscous flow through fibrous media with free surface | J| . international Journal N umerical methods in

Fluids J, 1994, 19(7): 575- 603. | 8| Lee L J, Young W B, Lin R J. Molding filling and cure modeling of RTM and SRIM processes | J|. Composite Structure J, 1994, 27(2): 109- 120. | 9| Rudd C D, Long A C, Kendall K N, et al. Liquid moulding technologies | M | . Cambridge E ngland: W oodhead Publishing Limited, 1997. 286- 288. | 101 Young W B. Three-dimensional non-isothermal mold filling simulations in resin transfer molding | J| . Polymer Composites J, 1994, 15(2): 118- 127. | 11 Chan A W, Morgan R J. Computer modeling of liquid composite molding for 3-dimensional complex shaped structures | A | . In: Proceedings of 19th ASM/ESD Advanced Composites Conference | C| . Dearborn, 1994. 341- 345. | 121 Kang M K, Lee W I. A flow front refinement technique for the numerical simulation of the resin transfer molding process | J| . Composites Science and Technology J, 1999, 59 ( 11) : 1663- 1674. | 131 Lim ST, Lee W I. An analysis of the three-dimensional resin transfer Mold filling process | J| . Composites Science and Technology J, 2000, 60(7): 961- 975. | 141 Dunn D. A chebyshev differentiation method | J| . Computer Physics Communications J, 1996, 96( 1) : 10- 16.

Biographies:

CHEN Ren -liang Born in 1963, he received B. S-, M.S. and Ph.D. in aircraft design from Nanjing Aeronautical

Institute in 1985, 1988 and 1998. From 2000 to 2002, he worked as post-doctoral in the field of RTM simulation at Florida Composite Center for Advanced Technologies ( FC2AT) of USA. Several scientific papers have been publish«! in various periodicals. Tel: ( 025 ) 4892141, 4892152, E-mail: crlae@nuaa.edu.cn GUI Bin Born in 1958, she received M.S. in «imputation mathematics from Nanjing Aeronautical Institute in 1988. She is the associate professor in Nanjing Forest University and the candidate Ph. D. student in computation mathematics at Nanjing University of Aeronautics & Astronautics. Tel: (025 ) 3735120

LI Ming-cheng Born in 1957, he receive! M . S. in mechanics from Nanjing Aeronautical Institute in 1988. He is the associate processor of mechanics in Nanjing University of Aa-onautics & Astronautics Tel: ( 025) 4892231, 4892152.

© 1994-2010 China Academic Journal Electronic Publishing House. All rights reserved, http://www.cnki.net