Three-dimensional finite element analysis of process-induced residual stress in resin transfer moldi

Joural of Harbini Instilute of Technology (New Series)， Vol. 15, No. 2, 2008Three-dimensional finite element analysis of process-inducedresidual stress in resin transfer molding processDAI Fu-hong, ZHANG Bo -ming, DU Shan-yi戴福洪，张博明，杜善义(Center for Compoite Materials, Habin Insitute of Technology, Harbin 150001 ,China, E maildaih@e hit. edu. cn)Abstract: A three-dimensional finite element analysis of process-induced residual stress in resin transfer mold-ing ( RTM) process is presented. The finite element method ( FEM) was employed to solve the coupled equa-tions involved in the transient heat transfer and the cure kinetics of the resin, and the distributions of intemnaltemperature and cure degree of the composite at any instant time were oblained. The slf-consistent field micero-mechanics model Was used to predict the cure-dependent mechanical properties of the composites. Thermal ex-pansion and cure shrinkage were included in the analysis. The thermo-elastic mechanical goveming equationswere solved using the incremental stress-strain relationship based FEM and the residual stress development waspredicted. The present resuls were validated by the comparisons with the pertinent lterature. The numericalexample of a half eylinder was presented. The results show that it is necessary to carry out the three-dimensionalanalysis due to the complex distributions of temperatures, cure degrees and process induced stress for thickparts, which can be predicted at any point within composite structures in the present analysis.Key words: cure ; residual stess ; finite element method; resin transfer molding ( RTM )CLC number: TB331Document code: AArticle ID: 1005-9113( 2008 )02-0271-06Determination of cure cycle and process-inducedcarbon/ epoxy laminates and it showed the initial size ofresidual stress are two of the most critical issues inlaminates and autoclave pressure can significantly afectcomposite manufacturing process. The poor design ofthe residual deformation.cure cycles may make the residual stress significantlyThe resin volume shrinkage is another importantdevelop tring the process, which can have a great im-mechanism contributions to process-induced stress be-pact on the performance of composite structures.sides thermal expansion mismatch. The effective me-Process-induced residual stress can reduce both ulti-chanical properties of resin can change drastically asso-mate strength and fatigue life, and is more important tociated with the cross-link polymerization reaction.composite parts with strict tolerance requirements. THTherefore，thermal and chemical strains and the cure-numerical simulation of process is a efective tool for re-dependent material properties should be included in thedueing process costs and improving part quality.analysis of residual stress. Bogetti' . .91 simulated the cu-The finite element method is widely used to simu-ring process and analyzed the curing process-inducedlate curing process for the sake of its adaplation to solvestress for autoclave thick section parts. It pointed outthe complex shape problems-4.。 These simulations of that the complex spatially varying thermal and degree ofcuring process can predict the distributions of tempera-cure gradients can lead to residual stress developmenttures and degree of cure within the composite materi-greatly. White[10.employed two dimensional simula-als. The residual stress analysis is necessary to carriedtion and visco-elastic model to study the curing processout if one wants to evaluate how these gradients of tem-and residual stess. Teplinsky [121 combined the one-di-perature and degree of cure influence on the quality ofmensional simulation withincremental laminate theoryparts. Chen'5 used a visco-elastic micro-mechanicalto predict the process-induced stress. ABQUS softwaremodel to study the effect of the cooling rate on the cu-was used to analyze the residual stress during autoclavering process induced residual stress. Andrew Johns-process for thick section woven fabric composite partstonto proposed a plane strain finite element model forby Huang'simulation of the process-induced deformation duringAlthough a quite amount of studies were made toautoclave processing. Tarsha Kurdi的investigated themodel compression molding process, little is aimed toinluence of autoclave pressure on curing stresses inmodel中国煤化工sS in resin transferReceived 2005 -04 -05.CNMHGSponsored by the National Natural Science Foundation of China( 'GranNo.105022016) a.. epnmn Lwvgau m munoranding Young Teachers inHarbin Insitute of Technology( Crant No. IITQNJS. 2006. 020).●271.Journal of Harbin Institute of Technology (New Series)，Vol 15, No.2, 2008molding ( RTM ) process. Golestanian!4J stated thated from Refs. [2, 11] and are listed as Eqs. (3,4):he first used the cure - dependent mechanical properties indc= kgexp(- AE./RT)a"(1 -a")，(3)predicting residual stresses in RTM process. However,his model is based on shell element , which neglect thede=(h, +ka")(1-a").(4)variation of temperature in the direction of thickness. Itdmay be successful in predicting residual stresses in thinwhere .laminates. The temperature, degree of cure and resid-h = A.exp(-△E /RT),ual stress can change from point to point within com-后= A2exp(-△E 2/RT).posite materials for thick section parts. It is necessarywhere R universal gas constant, T is temperature, ho,to develop three- dimensional model to achieve the spa-A1 ,A2,OE.,OE l,OE: are experimental constants.tially varying distributions of temperature and degree ofA variational principle is applied to Eq. (1) andcure and evaluate the residual stresses at any point.after the virtual temperature (8T) is introduced andThis paper is focused on developing three-dimen-volume integration is taken over the domain, Eq. (1)sional finite element program to simulate curing processis rewritten as Eq, (5):and analyze process induced stress within the complexshape thick parts. A three-dimensional finite element(.Cp )8TdV= [(V(k.V7)) +(1 -)q, )8TdV.formulation is proposed to avoid limitations of one-andtwo-dimensional analysis. The self-consistent field mi-The eight node three- dimensional elements is usedcro-mechanics model is used to predict cure dependentto discretized Eq. (5). Temperature (T) at any pointmechanical properties. The cure simulation and residu-within an element is interpolated as in Eq. (6):al analysis are performned and results are compared withthe relative literature.T= ENI.(6)where N; and T are shape functions and the nodal tem-1 Approachperature at ith node , respectively.The final discretized governing equation can be1.1 Cure SimulationThe thermal curing process of resin matrix com-expressed as follows:CT+KT= F.(7)posites is a themo-chemo coupled process. The resin istaken as not flowing in the stage of cuing and the con-where C is the thermal capacitance matrix, T is thevection thermal conduction is neglected. Thus the ther-nodal temperature vector, K is thermal conductivily ma-rix, F is thermal load vector. They are expressed asmal conduction equation with chemical reaction as writ-Eq. (8):ten in Eq. (1) can be used to describe a full three di-mensional cure processt4J :C= E fo.cenN.NdV,ρ.Cp i= V(k.VT) +(1-f)9, (1)K= 2(VNk。VN,)dV，wherep,, Cp are density, specific heat and k。is ther-(8)mal conductivity tensor, respectively. They can be de-F= 2」.N.(1-f) q,dV,termined according to mixture lawρ。=fpz+ (1 -f)p,T= aTat’: fpscp+(1 -f)p.Cp，Eq. (8) can be solved using the direct time inte-gration method. The thermal load F from chemical re-hjk,p.k。=fekhy+(1 -f)p,k,"action of resin is directly related to the degree of curewhere, f stands for volume fraction of fber, subscriptionin the resin. The degree of cure at the time step (n +1), C%+1 can be approximated as Eq. (9):r stands for resin, f stands for fiber, c stands for com-ar%+ =a.+(n) Ot.(9)posite. The intermal thermal resourceq, is the heat re-action of resin, described as Eq. (2):where an and (O) are degree of cure and the rate ofq, =ρ,H,la(2)cure at the time step (n).wherep, is density of resin, H, is total reaction heat of中国煤化工- use of the formula a-unit mass resin when resin finishes the curing reaction.boveresponding boundarycond;YHC N M H Gent must sficientlyis cure rate of resin, a is degree of cure.dismalr for accurate cure simulation, and is taken as 10 s inThe cure kinetics of resin in this paper are select-the present work , which gives a converged solution.●272●Journal of Harbin Institute of Technology (New Series)，Vol. 15, No.2, 20081.2 Material Model1.2. 1 Cure dependent resin modulus[Ca] = E U%[T]I[C][T].). (15)The modulus of resin is assumed to follow a simplewhere [Cen] is effective siffness matrix of unit cell,fmrule of mixture as shown in Eq. (10)[9)，is fiber volume fraction of mx sort composite, [T]m isE, = (1 -xmJ)E; +axmosE, ,(10)the coordinate transfer matrix of stress and strain be-a-Qe .tween the local coordinate system of of my sort compos-Qmol =(11)Xxend一age .ite and the global unit cell coordinate system.where E, and E, are the assumed fully uncured and1.2.4 Effective thermal expansion cofficients andfully cured temperature dependent resin modulus.shrinkage strain of compositesThe simple expression is used since it was shownThe thermal expansion and shrinkage stains areto offer an good representation of modulus for the resinalso dependent on the fiber and resin constituent prop-considered in this study 。The present model as-erties, and fiber volume fraction. A simple model onsumes the poisson ratio and the fiber properties arethe basis of a rule of mixture is used to predict the ef-constant during cure. The instantaneous resin shearfective thermal expansion coefficients and shrinkagemodulus during cure is determined based on the iso-strain of composites [9 .tropic material relation.syEf +8,E,(1 -f)E8: =EJ+E,(1-)(16)C, =2(2(1 +v,)(12)8r =e3 = (ey+Vnyey)f+ (er +128,)(1 -f) -I.2.2 Resin volumetric shrinkage modelResin shrinkage occurs during cure and provide[nf+na(1 -151.65.5.1-0A]. (17)Ef +E,(1 -J)an important source of internal loading. Assuming a u-where subscription 1 stands for longitudinal direction ,2niform strain contraction for all principal strain compo-stands for transverse direction. When the effectivenents, the incremental isotropic resin shrinkage strainthermal expansion cofficient are evaluated, the formu-of a unit volume element, Ae" , and the incrementallaan = ε,C2 = E2 is taken, and when the effectiveresin volumetric shrinkage have a relationship given byshrinkage strains are evaluated, the formulaei = 8，Eq. (13).&2 = E2 is taken.Oe"= (VT + OV,) -1.(13)The fiber is assumed to be zero shrinkage strain inA given incremental change in the degree of cure,this study.Sa, and the asociated incremental change in specific1.3 Process Induced Residual Stressvolume of resin, OV,, can be related to the total specif-The cure simulation is carried out to yield the tem-ic volume shrinkage of the completely cured resin, VTIperature and degree of cure distributions within the com-(as shown in Eq. (14)):posite parts in a single time increment before the stressOV, = AaVa.and deformation are calculated. The instantaneous efec-The cure shrinkage strain in the resin during curetive material properties and the resin shrinkage load areis the cumulative sum of all the incremental contribu-computed according to these distributions as above distions. The fiber itself is assumed not to undergo anycussed. These results are taken as the input parameterschemical contraction during cure.for residual stress calculations at the current time step.1.23 Effective elasic modulus of composies unit cellThe total incremental strain△e is given byThe composite mechanical properties strongly de-Oe=Oε'+Oe"+Oε'18)pend on the fiber and resin constituent properties, andwhere ε° is strain induced by mechanical load. Thefiber volume fraction. The self-consistent field micro-thermal strain Oe" can be expressed asmechanical model is widely used to compute the instan-Aε" = a●AT.(19)taneous mechanical properties for unidirectional fiberThe incremental form of stress and strain relation-reinforced composites'9. Huang[04) used TEXCADship is give bymodel to predict the effective unit cell modulus. We a-[Aσ] = [Ca]{[Ae*] +[Qe"] +[Qe]}. (20)dopt the similar method in this paper. At first the com-Eq. (20) is solved with the use of the finite elementposite unit cell is classified into N sorts of unidirection-method at each time step. The total strain can be ob-al reinforced composite (if necessary, the resin is alsotained taking the cumulative sum of incremental strain oftreated as one sort of such composite). Then the self-each time step. The details can be found elsewhere 19,13] .consistent model is employed to obtain effective me-chanical properties of Nu sort composite. Finally, the中国煤化工effective of unit cell is taken to be the superposition of2.1total N sorts of composite mechanical properties on the.CNM HGom Bgti 0basis of theirselves spatial directions.validate the present program. Material system is about.273.Joumal of Harbin Institute of Technology (New Series), Vol. 15, No.2, 2008glass/polyester. Tab: 1 gives the composite thermalTab.2 Cure kinetic parameters of Glass/Polyes-properties and Tab. 2 gives the cure kinetic parame-ter for two-dimensional exampleters. The characteristic values of polyester resin duringcure are shown in Tab. 3 and the mechanical proper-nko/(g-) SE。/(J. mol-')Hr/(J.kg-')ties of glass and polyester are listed as Tab.4. The ex-0.524 1.476 6. 1667x1021. 674 x 10s77500ample consists of a unidirectional laminate with the di-mensions of with 0. 1524 x0. 1524 x0. 0254 m. Themodel of cure kinetic is as shown in Eq. (3).Tab.3 Resin characteristic elastic modulus during cureTab.1 Thermal properties of Glass/Polyester for two-di-ProperticsE:/MPaEm /MPamensional examplePolyester2.7572. 757 xI0p%/(kg.m-3) cp/(J/(WC)) Kaz = Kssw/(mC) Kxu/Ka3Epoxy3. 4473.447 x10'189012600. 2163Tab.4 Fiber and resin constituent mechanical propertiesProperties E /MPaEr /MPa'121_ P423Gn /MPa Gr3 /MPa Gr3/MPa an/ 9a2/-1Glass7. 308x104 7.308x104 0.220.220.22 2.992x10* 2.992x10* 2.992x10* 5.04x10-6 5.04 x10-6Eq. (6)Eq.(6) 0.400.40 0. 407.20x10-5 7.20x10-9Eq.(6) 0.35 0.35 0. 35Eq.(6)5.76x10-55. 76x10-9The same boundary conditions are used in the ex-present solutions and the Bogetti's is %8. The reasonample , namely the mold temperature is applied at theis that Bogetti used the laminated theory, while thetop and bottom surface of the part and the isolate heatthree-dimensional finite elements are used in the presentcondition is applied at the lateral surfaces. The temper-analysis. The significant self-equilibrating stresses re-ature history at the centroid of the part is as shown inmain after complete cure. The magnitude of the resinFig.1. The results agree well with the Bogetti's. Thvolumetric shrinkage strongly affects the stress develop-present predicted temperature peak is 126.9 C, where-ment. These results are agreement with Bogetti" s.as the Bogetti's is 126 C or so. Since the exact curecycle curve cannot be obtained , the inflexion of the time0.5in the cure cycle curve may be a lttle different from itsoriginal value of the Bogti's. The eror leads to that旨0.3the time when the temperature peak occurs is 168 min(Bogetti's is 164 min) for the present simulations.6%(3% )2.0%14_ Presant 3D FEM12-0.3 |- Bogrti brinedplate theory10- 3D-0.5参80-10-8-6-4-2024681012Transverse stress/MPa6Fig.2 Comparisons between the present analysis with the4(Bogetti' s results for Example 12(2.2 Example Two100200300400Material system is about glass/ epoxy, and fibert/ minvolume fraction is 60%. The effective unit cell is asFig. 1 The history of temperature of example oneshown in Fig. 3. The simple unit cell used here is forThe in-plane transverse stress for three diferentits very small thickness of 0.18 mm. The third exam-volume shrinkage rates 6% ,3% ,1% and 0% are shownple is a half evlinder with a height0. 1 m, an inner ra-中国煤化工0.005. A tolal volu-in Fig. 2. A lttle difference between the present solu-ned to be 5%. Tab. 5tions and the Bogetti's is observed, especially near thefYHCN M H Gopeties and Tab. 6mid-plane of the part. The maximum error between thegives the cure kinetic parameters. The characteristic●274●Journal of Harbin Institute of Technology (New Series)， Vol. 15, No. 2, 2008values of epoxy resin during cure are seen in Tab. 3and the mechanical properties of glass and epoxy arelisted as Tab.4. The model of cure kinetic is as shownin Eq. (4).Tab. 5 Thermal properties of Glass/Epoxy for half cylin-der examplep/(kg.m -3)C,/ (J/(W . C))Kw/(m.C)Class2560712. 350. 137Fig.3 The effect unit cell of composite materials for halfEpoxy115082108. 673.cylinder exampleTab.6 Cure kinetic parameters of Glass/Epoxy for half cylinder exampleAA:OEErHrm/(s-1)/(J●mol-' )/(] . mol-1)/(J.kg")0.9580.841 .1.759 x10*3. 888x10%6.513 x10f5.408 x10f459224In residual stress analysis, perfect bonding be-rections on the global coordinate system ( as shown intween the laminate and the mould wall is assumed.Fig.6). The deformed shape is plotted ( as shown inThough the assumption still is a quite open issue, hereFig. 7). The maximum displacements along three dif-we follow the previous convention to take it. Theferent directions, U,,U,,U, are 1. 11 ,0. 545 and 0. 59process rinduced residual stresses are solved as discussmm respectively.above at each time step. The tolal stress is the cumula-tive sum of each step solutions.180The developments of the temperature are shown in160Fig.4. It can be seen that temperatures at central loca-40 t二. at centertions are higher than those at the positions close to themold120surface. It resulted in a temperature difference of about24.2 C in the thickness direction. The significantR 100temperature overshoots are observed from 36 min. The8Qtwo different distributions of the temperature at the time60of30 and 38 min are shown respectively in Fig. 5. It40can be found that the regions with highest temperaturetends to offset towards the center as the resin cures. Italso indicates that the fraction of residual stress before0 50~100150200 250 300the cooldown is about 50 percent of its maximum val-t/ minue. The stress σ is greatest among the three axial di-Fig.4 The history of temperature for half cylinder exampleA=80.363.A=80.825.182-81.089V9=82475C-81.815”C=84.125D=82.540D)=85.775E=83.266E=87.425F-83.992F-89.074G=84.718C=90.724H=85.444H=92.3741-94.024中国煤化工(a)t = 30 minYHCNMHGFig.5 The distributions of temperatures at tume 30 min and 38 minJoumal of Harbin Institule of Technology (New Series), Vol. 15, No.2, 2008200of temperatures, curedegrees and process-induced stress for thick parts,150---0.which can be predicted at any point within compositestructures in the present analysis.令100-References:[1] Choe M A, Lee M H, Chang J. Thre-dimensional simula-tion of the curing step in the resin transfer molding process.Polymer Composites, 1999 ,20(4) :543 -552.[2] Liu Xiaolin, CrouchI C, Lam Y C. Simulation of heattransfer and cure in pulrusion with a general-purpose finite-505C100 150 200 250 300element package. Composites Sciencend Technology，2000 ,60:857n/min[3] Park H C, Goo N s, Min K J, et al. Three-dimensionalFig.6 The history of reidual stresses for half eylinder ex-cure simulation of composite structures by the fnite elementamplemethod. Composite Structure , 2003 ,62:51 -57.[4] Park H C, Lee S W. Cure simulation of thick compositestructures using the finite element method. Journal of Com-posite Materials, 2001 ,35(3) :188 - 201.[5] Chen Yu, Xia Zihui, Elyin F. Evolution of residual stres-ses induced during curing processing using a viscoelasticmicromechanical model. Joumal of Composite Materials,2001 ,35(6) :522 -542.[6] Johnston A, Vazini R, Poursartip A. A plane strain modelfor process-induced defornation of laminated compositestructures. Joumal of Composite Materials , 2001 ,35(16):1435 - 1469.[7. ] Tarsha-Kurdi K E, Olivier P. Thermoviscoelastie analysisof residual curing stresses and the influence of autoclavepressure on these stresses in carbon/epoxy laminales. Com-posites Science and Technology , 2002 ,62 :559 - 565.Fig.7 The deformation of half cylinder example[8] Bogeti T A, Gillespie Jr J W. Two-dimesional cure sirmula-tion of thick thrmoseling composites. Joumal of Composite3 ConclusionMaterials, 1991 ,25(3) :239 -273. .9] BogttiT A, Gillespie JrJ W. Process-induced stress andA three- dimensional finite element analysis of cu-deformation in thick-section thermoset composite laminates.Joumnal of Composite Materials , 1992 ,26(5) :626 - 660.ring and process induced residual stress in Resin[ 10]While s R, Hahn H T. Proces modeling of composite ma-Transfer Molding ( RTM) process is presented. The fi-terials: residual stress development during cure( Part 1):nite element method ( FEM) was employed to solve theModel formuation. Joumal of Composite Materials, 1992,coupled equations involved in the transient heat trans-26( 16) :2402 -2422.fer and the cure kinetics of the resin and the distribu-[11]White s R, Hahn H T. Proces modeling of copoite ma-tions of intemnal temperature and cure degree within theterials: reidual stress development during cure(Part II):composite at any instant time are obtained. A compari-Experimental validation. Journal of Composite Materials ,1992 ,26( 16) :2423 -2453.son with the Bogetti ' s results validated the present[ 12 ]Teplinsky s, Gutman E M. Computer simulation of processsimulations. The self-consistent field micro-mechanicsinduced stress and strain development during cure of thick-model was used to predict the cure dependent compos-section thermosetting composites. Computational Materialsite mechanical properties. Thermal expansion and cureScience , 1996 ,6:71 - 76.shrinkage were included in the analysis. The thermo-e-[ 13 ]Huang Xiaogang, Gillespie JrJ w, Bogetti T. Process in-lastic mechanical goverming equations were solved withduced stress for woven fabrie thick section composite struc-the use of the incremental stress-strain relationshiptures. Composite Structures, 2000, 49 :303 - 312.based FEM and the residual stress development was[14] Golestania H, El-CGizawy A s. Modeling of process induced residual stresses in resin transfer modeled compositespredicted. The numerical examples including a halfwith woven fber mats. Journal of Composite Materials,cylinder were presented . The results show that it is2001 ,35(1):1513 - 1528.necessary to carry out the three dimensional analysis中国煤化工MYHCNMHG●276●