原稿.pdf

1、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

2、 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

3、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 ori

4、entations 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 mixtu

5、re 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

6、 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 i

7、n 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

8、 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=Da

9、k.(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

10、 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 o

11、f 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 (kine

12、matic 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 Damk

13、ohler 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

14、, 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

15、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

16、), 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

17、- 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 functio

18、ns 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

19、 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

20、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;

21、 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 surfac

22、e 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

23、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 successful

24、ly 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

25、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 deri

26、vatives 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 rig

27、ht-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

28、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

29、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 R

30、ISC/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-di

31、rection 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

32、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

33、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 po

34、lyurethane 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 a

35、spect 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

36、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 (

37、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 symb

38、ols 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

39、 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 substant

40、ially 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. T

41、hese 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 pro

42、files 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 convention

43、al 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

44、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 inject

45、ion 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

46、 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

47、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 contra

48、st 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 ac

49、curacy 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-

50、 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

51、 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 Gal

52、erkin 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.

53、 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 (198

54、7). 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

55、, 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-the

56、yre 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.

57、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. Mizukam

58、i, 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,

59、 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 Gale