原稿.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=IvReg +v.VvWfdV (12) + y-PI+(K+)VWidV- s s n.-pI+(lcf)WdS=O, (13) s -Brrc(j:Vv)-Dak,(l-X)” WdV 1 + s VT.VWdV- s (n.VT)WdS=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 the 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, Dont suppress wiggles-theyre 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.