原稿.pdf

Pergamon Computers in revised form 4 May 1994 Abstract A numerical scheme for computing the advancement of a flow front and related velocity pressure confersion and temperature distributions during mold filling in reaction injection molding RIM is described in this work In the RIM process the convective term in the energy equation is dominant Therefore the numerical scheme has incorporated a Petrov Galerkin finite element method to suppress spurious oscillations and to improve accuracy of the calculations The other feature of the numerical scheme is that the flow front locations are computed simultaneously with primary variables by using a surface parameterization technique The numerical results compare well with the reported experimental data Improved accuracy obtained by this numerical scheme in the flow front region is expected to assist in the predictions of the fiber orientations and the bubble growth in RIM which are determined primarily by the flow front region I INTRODUCTION Reaction injection molding RIM is a widely used process to manufacture exterior fascias in the automobile industry In this process a prepolymerized isocyanate and a polyol amine mixture are mixed together and injected into a mold where polymerization occurs A fountain flow effect in the advancing flow front region during the mold filling stage plays an important role in determining the residence time of the fluid elements and in controlling the fiber orientations in the final product 11 An accurate s imulation of this flow front however poses a challenging problem Evolving flow domain with advancing flow front requires updating of the numerical grids and prediction of the moving boundary at every time step Low thermal conductivity of the material high flow rates in the RIM process and highly exothermic rapid reactions result in convection dominated energy transport equation which needs a special numerical treatment Besides moving contact lines near the walls need suitable boundary conditions that do not introduce numerical instability A numerical scheme that incorporates all these complex features of the RIM process is required for accurate predictions near the flow front region Previous studies either have made simplifying assumptions regarding the flow front region 2 3 conservation of momentum equation Re v Vv pV I v rcj Gz7 v VX Dak l X mole balance equation 4 5 conservation of energy equation Gz g v VT V T Brrc j Vv Darc l X m L 1 6 where v is the velocity vector q the rate of strain tensor t the time p the pressure and k is the dimensionless rate constant defined as exp E R l T l T The equations are made dimensionless using the average velocity V half of the thickness of the mold H and the temperature T and the viscosity qO r X 0 T T at the inlet of the mold All the dimensionless groups and their definitions are listed in Table 1 The boundary conditions in terms of dimensionless variables are 1 at the walls v 0 no slip T T 2 at the mid plane aTjay 0 3 at the inlet v fully developed flow T 1 X X 4 at the contact line n PI 2 0 full slip 5 at the flow front n PI 2 0 force balance n v ah at 0 kinematic condition Table 1 Dimensionless groups in governing equations where AH is the heat of reaction AT the adiabatic temperature rise and C the initial concentration of isocyanate GZ Graetz number VHpC lk Re Reynolds number HVlrlo K viscosity ratio 41 Br Brinkmann number to V lkT Da Damkohler number AH H C kT A exp E RT T adb adiabatic temperature rise AT IT Flow front advancement in reaction injection molding 57 where a and vY are the components of the velocity vector v II the unit normal vector r the extra stress tensor h the location vector of the flow front and Twal the dimensionless temperature at the mold wall The details of incorporating the boundary conditions in the numerical analysis are explained in the next section 3 NUMERICAL ANALYSIS In the finite element formulation the unknown velocities temperature and conversion are expanded in terms of the biquadratic basis functions 4 the pressure in terms of the bilinear basis functions ll i and the flow front shape h in terms of the quadratic basis functions 7 where l and q are the coordinates in isoparametric transformation defined as i 1 i l in the isoparametric domain 1 4 1 1 q L 5 where Pe is the local element Peclet number VA D A the element size and c s is the cubic polynomial 5 8 5 l t 1 The index i 1 corresponds to the vertex nodes and i 2 58 NITIN R ANTIJRKAR corresponds to the centroid nodes in the element The standard one dimensional convec tiondiffusion problem has exact solution at the nodes if 25 9 c Pe 2 tanh Pe 2 l 3 Pe coth Pe 4 X Pe coth Pe 4 1 la c2 16 Pe 4 coth Pe 4 1 lb In a two dimensional problem the tensorial product of equations 10 and 11 provides the function c in the weighting functions described in equation 9 The local Peclet number is computed for each three node group based on the average velocities at the relevant boundaries in the two dimensional element 9 There are six such groups three in the x direction and three in the y direction and thus there are 12 upwinding parameters E The calculations of the Peclet number involve linear distances which essentially neglect the curvilinear sides of the elements However it is a good approximation since flow front is not severely deformed in our problem The diffusivities are l Gz for the energy equation and is K R for the momentum equation The Petrov Galerkin weighted residual equations are R V v dl O s RL IvRe g v Vv WfdV 12 y PI K VWidV s s n pI lcf W dS O 13 s Brrc j Vv Dak l X W dV 1 s VT VW dV s n VT W dS O 15 V S RI s n v ah for v T and X at the inlet of the mold and for vY at the mid plane axis of symmetry are applied by substituting the boundary conditions for the equations The natural boundary conditions namely the symmetry conditions at the mid plane the full slip zero friction condition at the contact point and the zero force at the free surface are implemented by substituting the boundary terms in the residual equations The kinematic boundary condition at the flow front is incorporated as the governing equation for predicting the flow front locations The weak form of energy equation is extended to the flow front boundary by evaluating the boundary terms instead of by imposing any unknown essential or natural boundary conditions 27 Such free boundary condition as denoted by Papanastasiou et al 27 minimizes the energy functional among all possible choices at least for various types of creeping flows and has been successfully used in several applications including those with high Reynolds numbers The spatial discretization reduces the time dependent equations 12 16 to a system of ordinary differential equations M 2 R q 0 17 Flow front advancement in reaction injection molding 59 where q is the vector of all nodal unknowns such as pressure velocities along x and y directions temperature conversion and locations of the flow front The matrix M is the mass matrix and R the residual vector The time derivatives are discretized by a standard first order implicit scheme Note that the temporal derivatives need to be adjusted for moving tessellation according to aq dq aq dx 3 at dt 8X dt 18 19 where left hand side represents the local change of variable with time the first term on tilt right hand side rlepresents the total change of variable with time while the second term on the right hand side represents convective changes due to moving finite element grids When equation 19 is substituted in equation 18 a nonlinear algebraic system of equations is obtained These equations are solved by the Newton iteration scheme Since the flow front locations themselves are unknowns in these equations one must be careful in obtaining the derivatives of the residual equations with respect to the flow front locations as the Jacobian of the isoparametric mapping also depends on these locations 28 The linearized equations are solved using the frontal numerical technique 29 The tessellation is updated at each iteration by the newly found flow front location values which are determined simultaneously with other variables The solution at previous time step is used as initial guess in the Newton iteration for the next time step The computations were carried out on an IBM RISC 6000 machine In all the calculations a dimensionless time step of 0 05 and seven elements along the transverse direction with finer tessellation near the wall were sufficient to give solutions that are independent of mesh refinement and time step size Initial length of an element along x direction was 0 1 The mesh was regenerated when the length of the element exceeded 0 2 4 RESULTS The axial velocity at the tip of the flow front is expected to be equal to the average velocity of the flow front In all our calculations the dimensionless axial velocity u was equal to 1 000 with a tolerance of 0 1 In the numerical computations the reaction rate can be forced to be zero If the mold temperature is the same as the material temperature then the results are expected to simulate isothermal injection molding process The shape of the flow front in this simulation was identical to that reported in the earlier study on injection molding of Newtonian liquids l Numerical results were further validated by comparing them with the experimental pressure rise data 2 at the inlet of the rectangular mold Castro and Macosko characterized rheokinetic and thermal properties of two polyurethane RIM systems in their paper We performed three numerical runs corresponding to three of their experiments involving these RIM systems at various filling times refer to Table 2 for details Although the filling times in the numerical runs were identical to those in the experiments the aspect ratio 6 L 2H of the rectangular mold was kept at 25 to save computational time When the aspect ratio was increased in proportion to the average velocity while keeping the filling time constant the differences in the numerical results of velocities temperatures and conversions were less than 5 The predicted pressure rise data was scaled by Table 2 Summary of numerical experiments where X is the maximum conversion r the maximum temperature and r the minimum temperature in the mold at the end of filling Experiment 1 9 4 l 2 912913 3 6 738 Material exp exp 2200 2H cm 0 32 0 32 0 64 t 25 0 25 0 25 0 Filling time set I 11 2 65 2 46 T C 74 I 21 3 65 0 To 49 3 50 8 60 0 X 0 85 0 85 0 65 X mtix 0 05 0 14 0 60 T I 08 I 03 I 18 T nil 1 00 0 93 1 00 60 50 r Run NITIN R ANTURKAR P Time s Fig 1 Comparison of the experimental data represented by symbols with the numerical predictions represented by lines the appropriate aspect ratio of the experimental investigation The comparisons of the measured and predicted pressure rise at the inlet of the rectangular mold are shown in Fig 1 Note that the increase in viscosity is marginal for first two experiments because the filling times are much smaller than the gelling times Therefore the pressure rise at the entrance of the mold is linear and represents the length of the mold filled However the injection rate is much slower in the third experiment Therefore viscosity increases substantially during mold filling and the pressure rise curve is not linear In all these cases the predictions compared well with the experimental data In fact the predictions were better using our model than those in 2 for experiment 3 in which substantial extent of reaction occurs during filling These better agreements are probably due to the accurate simulation of the flow front and that of the heat transfer in thicker mold without any simplifying assumptions The accuracy of the Petrov Galerkin method is compared with conventional Galerkin method in Fig 2 In this figure temperature profiles along the flow direction at various transverse locations are plotted for the two formulations when the flow front reaches at the end of the mold Severe artificial oscillations are observed near the flow front at all transverse locations when the temperature profiles are obtained by conventional Galerkin formulation These oscillations disappear when the Petrov Galerkin formulation is used with identical mesh refinement Thus the Petrov Galerkin formulation is clearly superior for accurate calculations in the advancing flow front in RIM h 1 16 I 08 Num Method Petrov Galerkin I oo I I I I I I I I 0 60 0 65 0 70 0 75 0 80 0 85 0 90 0 95 1 oo x e flow direction Fig 2 Comparison of temperature profiles along the flow directon at various transverse locations obtained by the Galerkin and Petrov Galerkin formulations for experiment 3 Flow front advancement in reaction injection molding 61 0 8 0 6 I I I I I 0 0 77 0 78 0 79 0 80 0 81 0 82 0 83 X E flow direction Fig 3 Prediction of the shape of the flow fronts just before the end of mold filling for experiments 1 and 3 Finally representative results of the front profiles corresponding to experiments 1 and 3 at the end of mold filling are illustrated in Fig 3 As material advances in the mold the maximum temperature moves from the wall to the center due to the heat evolution from the reaction which can not be removed by conduction The conversion follows the temperature profile However the viscosity is strongly dependent on the conversion and the rise in temperature near the wall is not large enough to offset the viscosity rise due to higher conversion Therefore there is more resistance to the flow near the wall This surging effect is evident in the front profile of experiment 3 In contrast to the previous assumptions 2 4 those flow fronts are not flat which is expected from the experimental observations 11 5 CONCLUSIONS Accurate predictions near the flow front region play a major role in predicting fiber orientations and bubble growth in reaction injection molding RIM The accuracy is improved in this model by incorporating three important features 1 No assumptions were made regarding the shape of the flow front or the distribution of primary variables in the flow front region 2 The Petrov Galerkin formulation was incorporated instead of conventional Galerkin formu lation to avoid spurious oscillations due to dominant convection terms in the energy equation 3 The kinematic boundary conditions on the free surface is incorporated as a governing equation and thereby the shape of the flow front is calculated simultaneously along with other variables using t he free surface parameterization technique The numerical results are validated by comparing them with the experimental pressure rise data The results are in excellent agreement with the experiments even close to the gel point The spurious oscillations observed near the flow front region in the Galerkin method are suppressed by using the Petrov Galerkin method which is superior for accurate predictions of the advancing flow front in RIM 1 2 3 4 5 6 1 REFERENCES D J Coyle J W Blake and C W Macosko The kinematics of fountain flow in mold filling AZChE J 33 1168 1987 J M Castro and C W Macosko Studies of mold filling and curing in the reaction injection molding AIChE J 28 250 1982 N P Vespoli and C C Marken Heat transfer and reaction effects during mold filling of fast reacting polyurea RIM systems Presented at the Annual AIChE Meefing New York 1987 M A Garcia C W Macosko S Subbiah and S 1 Guceri Modeling of reactive filling in complex cavities Int Polym Processing 6 73 1991 C D Lack and C A Silebi Numerical simulation of reactive injection molding in a radial flow geometry Polym Eng Sci 28 434 1988 R E Hayes H IH Dannelongue and P A Tanguy Numerical simulation of mold filling in reaction injection molding Polym Eng Sci 31 842 1991 N R Anturkar An advancing flow front in RIM Polym Eng Sci in press 1995 62 NITIN R ANTURKAR 8 P M Gresho and R L Lee Don t suppress wiggles they re telling you something In Finite Element Methods for Convection Dominated Flows Edited by T J R Hughes AMD Vol 34 ASME New York 1979 9 P M Adornato and R A Brown Petrov Galerkin methods for natural convection in directional solidification of binary alloys Int J Numer Meth Fluids 7 761 1987 10 A N Brooks and T J R Hughes Streamline upwind PetrovGalerkin formulation for convection dominated flows with particular emphasis on the incompressible Navier Stokes equations Comput Meth Appl Mech Eng 32 199 1982 11 T J R Hughes M Mallet and A Mizukami A new finite element formulation for computational fluid dynamics II Beyond SUPG Comput Meth Appl Mech Eng 54 341 1986 12 A C Galeao and E G Dutra do Carmo A consistent approximate upwind PetrovGalerkin method for convection dominated problems Comput Meth Appl Mech Eng 68 83 1988 13 P A Tanguy and J M Grygiel A slightly compressible transient finite element model of the packing phase in injection molding Polym Eng Sci 33 1229 1993 14 J M Grygiel and P A Tanguy Finite element solution for advection dominated thermal flows Comput Meth Appl Mech Eng 93 277 1991 15 T J R Hughes L P Franca and G M Hulbert A new finite element formulation for computational fluid dynamics VIII The Galerkin least squares method for advectivediffusive equations Comput Meth Appl Mech Eng 73 173 1989 16 C W Macosko Fundamentals of Reaction Injection Ma ding 1st edition Hanser New York 1989 17 D F Mielewski N R Anturkar and D R Bauer The rheokinetics of RIM systems in solution ANTEC 92 Proceedings of the Socie