Printer Friendly

A model of advancing flow front in RIM.


Reaction injection molding (RIM) is a widely used process to manufacture exterior fascias in the automobile industry. In this process, a prepolymerized diisocyanate, and a mixture of prepolymerized polyol and amines are intimately mixed by impingement mixing. This low viscosity reactive mixture is then injected into a mold. An accurate simulation of this process poses a challenging problem. Some of the key features that need to be accounted for in the simulation, are: (1) simultaneous heat, mass, and momentum transfer, and resulting changes in the physical properties of the system; (2) evolving flow domain with contact lines at the walls; (3) dominating convective terms in energy transport; and (4) a necessity to know the space-time history of fluid elements. Although the flow front region is small, it controls the final residence time of the fluid elements. When fibers are added in the RIM system, the flow front region also determines the fiber orientation in the final product. To the best of our knowledge, none of the RIM models that were reported in the literature or are commercially available, account for all the above features. Most of the studies either made simplifying assumptions for computation in the flow front region (1-3), or did not compare their results with experiments (4, 5).

In this work, a two-dimensional model is developed for mold filling in RIM using a Petrov-Galerkin finite-element method with free surface parameterization. Dependence of viscosity on the conversion and temperature is represented by the Castro-Macosko model. To facilitate accurate predictions, no a priori assumptions are made regarding the shape, or the velocity distribution in the flow front region. The results are compared with the available experimental data (1).


The assumptions in the model are: constant thermal conductivity k, heat capacity [C.sub.p], and density [Rho]; laminar, two-dimensional flow at low Reynolds number, and at constant flow rate Q; negligible molecular diffusion and body forces; and uniform mixing at ideal stochiometric ratio of the reactants before entering the mold. We neglect microcellular foaming in RIM. Therefore, the assumption of incompressibility of reactive mixture is reasonable.

A kinetic rate expression for RIM reactions has been given by:

[Mathematical Expression Omitted],

where, [C.sub.i] is the diisocyanate concentration, R the gas-law constant, m the order of the reaction, [E.sub.r] the activation energy of the reaction, and [A.sub.r] the rate constant. Kinetics of all the complex reactions in RIM are lumped together in the above expression for process modeling. Although second-order rate expression is applicable to polyurethane chemistry, it may not be valid for polyurethane/polyurea and polyurea chemistries (6). The reaction between diisocyanate and aliphatic amines is assumed to be instantaneous in polyurea RIM. Thus, the reactive mixture enters the mold with conversion [X.sub.l] and temperature [T.sub.o], which have to be determined independently by mass and energy balance.

In general, inelastic generalized Newtonian model of viscosity, i.e., [Mathematical Expression Omitted] is appropriate, where [Tau] is the extra stress tensor, [Mathematical Expression Omitted] the rate-of-strain tensor, and [Eta] the viscosity. Although non-Newtonian effects are qualitatively observed in the previous work (7), these effects have not been incorporated in RIM simulations due to lack of quantitative data. The shear-rate dependence of viscosity is expected to be more pronounced in polyurea systems, in which reactive mixture entering the mold is already partially reacted. Since experimental characterization of non-Newtonian effects is unavailable, the results presented here are for Newtonian liquids. The Castro-Macosko viscosity function (1)

[Eta](X, T) = [Eta](X) [center dot] [Eta](T) = [A.sub.[Eta]] exp ([E.sub.[Eta]]/RT) [([X.sub.g]/[X.sub.g] - X).sup.(A + BX)], (1)

is used in all the results, where [X.sub.g] is the gel conversion, and [A.sub.[Eta]], [E.sub.[Eta]], A, and B are the constants.

The dimensionless equations of change are continuity equation:

[Nabla] [center dot] v = 0;(2)

conservation of momentum equation:

[Mathematical Expression Omitted];

mole balance equation;

Gz[T.sub.adb] [[Delta]X/[Delta]t + v [center dot] [Nabla]X] = [Dak.sub.r] [(1 - X).sup.m]; (4)

conservation of energy equation:

[Mathematical Expression Omitted];

where, v is the velocity vector, t the time, p the pressure, T the temperature, and [k.sub.r] is the dimensionless rate constant, defined as exp[(-[E.sub.r]/R)(1/T-1/[T.sub.o])]. The above equations are made dimensionless using the average velocity V, the half of the thickness of the mold H, and the temperature [T.sub.o] and the viscosity [[Eta].sub.o](= [Eta](X = [X.sub.l], T = [T.sub.o])) 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.sub.wall] = 0 (no-slip), T = [T.sub.wall];

2. at the mid-plane: [Delta]T/[Delta]y = O, [Delta][v.sub.x]/[Delta]y = O, [v.sub.y] = 0;

3. at the inlet: v = fully developed flow, T = 1, X = [X.sub.l]; 4. at the contact line: n [center dot] (-pI + [Tau]) = 0 (full-slip);

5. at the flow front: n [center dot] (-pI + [Tau]) = 0 (force balance), n [center dot] (v - [Delta]h/[Delta]t) = 0 (kinematic condition):


where [v.sub.x] and [v.sub.y] are the components of the velocity vector v, n the unit normal vector, and h the location vector of the flow front.

The essential boundary conditions for v and T at the walls, for v, T, and X at the inlet of the mold, and for [v.sub.y] 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 equations for predicting the flow front location. 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 (8). Such "free boundary condition" 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.


Since adequate characterization of polyurethane/polyurea or polyurea materials, and associated experimental data are not available, the computations were carried out for polyurethane systems, characterized by Castro and Macosko. Rheological, thermal, and kinetic parameters of two different RIM systems (experimental and RIM 2200) have been reported in their paper (1). We performed four numerical runs, which correspond to four of their experiments, identified in brackets next to Run #s in Table 2. The filling times in the numerical runs were identical to those in the experiments. However, the aspect ratio [Epsilon] (= L/2 H) of the rectangular mold was kept at 25 to save computational time. When the length of the mold was increased ([Epsilon] = 400) in Run #3 by increasing average velocity for the same filling time, the differences in conversion, velocities, and temperatures at any identical x/[Epsilon], were less than 5%. In other words, the results are similar for molds with various lengths as long as filling times are identical.

The Petrov-Galerkin isoparametric finite-element method with free surface parameterization is used to solve the governing equations (Eqs 1 to 5) with above-mentioned boundary conditions. All variable are expanded in terms of biquadratic basis functions, except pressure, which is expanded in terms of bilinear basis functions. The locations of the flow front are computed simultaneously along with other variables by satisfying the kinematic boundary condition at the flow front. The y-coordinates of all the nodes are fixed, and the x-coordinates are left as unknowns. The resulting ordinary differential equations with respect to time are solved by simple, first-order, implicit time-marching scheme with full Newton iteration. The frontal routine is used to solve the final set of linear equations.
Table 2. Summary of Numerical Runs.

Run# 1(9/4/1) 2(9/29/3) 3(6,7,8) 4(9/21/2)

Material exp. exp. 2200 exp.
2H, cm 0.32 0.32 0.64 0.16
[Epsilon] 25.0 25.0 25.0 25.0
filling time, s 1.11 2.65 2.46 5.72
[T.sub.wall], [degrees] C 74.1 27.3 65.0 70.0
[T.sub.o] 49.3 50.8 60.0 54.0
[X.sub.g] 0.85 0.85 0.65 0.85
[X.sub.max] 0.05 0.14 0.60 0.78
[T.sub.max] 1.08 1.03 1.18 1.24
[T.sub.min] 1.00 0.93 1.00 1.00

[X.sub.max] = maximum conversion.
[T.sub.max] = maximum temperature.
[T.sub.min] = minimum temperature in the mold at the end of filling.

The formulation of the Galerkin finite-element method is customary, and can be found elsewhere. However, a regular Galerkin technique is known to be numerically unstable for convection dominated equations. Since Graetz number in energy equation is usually high in RIM, an artificial oscillating temperature profile is computed along the flow direction. These oscillations can be resolved by a very fine tessellation. However, such time consuming approach can be avoided by using an upwinding numerical scheme. such as the Petrov-Galerkin method. The scheme involves introduction of artificial diffusivity in the upstream direction through weighting functions in the Galerkin formulation. A typical weighting function is,

[Mathematical Expression Omitted]

where, [[Phi].sup.i] is the biquadratic basis function used to expand unknown velocities, temperature, and conversion, [Xi] and [Eta] the coordinates in the isoparametric transformation and the artificial diffusivity is introduced through the function [Zeta], which depends on the local velocity field, and the appropriate diffusivity D associated with each governing equation. The diffusivities are 1/Gz for the energy equation and [Kappa]/Re for the momentum equation. For other governing equations, the basis functions themselves are used as weighting functions. Additional details of the Petrov-Galerkin method are reported elsewhere (9). A typical example of the temperature profiles obtained by these two formulations are shown in Fig. 1. The Petrov-Galerkin formulation is clearly preferable.

When the flow front moves, all the nodes in the flow domain are proportionately moved. When the length of any element exceeds predetermined value, mesh is regenerated by dividing each element into two along the flow direction. Interpolation of variables to new nodes is simple, and can be rapidly computed for biquadratic elements. The initial flow domain is a single row of elements across the entrance of the mold. When the advancing flow front hits the wall at the end of the mold, the adaptive time-step procedure is used to minimize the artificial loss of mass through the walls.

The computations were carried out on an IBM RISC/6000 machine. In all the calculations, a time-step of 0.05 and seven elements along the transverse direction (with finer tessellation near the wall) were sufficient to give mesh- and time-step-independent solutions. Initial length of an element along the x-direction was 0.1. The mesh was regenerated as soon as this length exceeded 0.2. When the aspect ratio [Epsilon] of the mold was equal to 25, the required CPU time was approximately 25 min for any typical run listed in Table 2. As the length of the mold increases, the increased CPU time is not only due to increased volume to fill, but also due to more number of elements required in a larger flow domain. Since finer tessellations are required for Galerkin formulation, the CPU time required is at least 15 times higher than that required for Petrov-Galerkin formulation for computations. The CPU time required for each time-step at the initial stages of mold filling is small due to less number of elements in a smaller flow domain.


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 [v.sub.x] was equal to 1.00 with a tolerance of 0.01%. Besides, when the calculations were carried out in the absence of reaction and heat transfer, the shape of the flow front was identical to that reported in the earlier study on Newtonian liquids (10). Numerical results were further validated by comparing them with the experimental data on the pressure rise (1) at the inlet of the mold. These comparisons are shown in Fig. 2 for three different runs. The values of [X.sub.max] and [T.sub.max] in Table 2 indicate that viscosity essentially was constant for Run #1, whereas, it was close to the gel conversion for Run #3. The initial nonlinear pressure rise in Run #2 was due to increase in viscosity associated with decrease in temperature ([T.sub.wall]/[T.sub.o]). In all these cases, the predictions compared well with the experimental data. In fact, the predictions were better using our model than those in Ref. 1 for Run #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 heat transfer in thicker mold without any simplifying assumptions.

The temperature and conversion profiles for several axial positions are plotted in Figs. 3 and 4 respectively for Run #4 at the end of mold filling. Qualitatively, the temperature profile is similar to that predicted in Ref. 1. As material advances in the mold, the maximum temperature moves from the wall towards the center due to heat evolution from the reaction. The temperatures are much lower near the wall at x/[Epsilon] = 0.92, than those in Ref. 1 model. Note that the boundary condition of T = [T.sub.wall] used in our model, cannot be enforced in the simplified flow front region in Ref. 1. Therefore, the temperatures of fluid elements near the wall are closer to [T.sub.wall], particularly for thin molds in Run #4. The conversion profile shown in Fig. 4 is similar to the temperature profile.

The axial velocity along the transverse direction is plotted in Fig. 5 at various locations along the flow direction for Run #3. These velocity profiles are plotted for constant node numbers along x, rather than for constant x/[Epsilon]. Since the flow front itself is not flat, the constant node numbers along x do not correspond to constant x (refer to the section on Numerical Analysis). The velocity profile at the inlet is the standard, parabolic Newtonian velocity profile at constant viscosity. However, similar to Fig. 4, the conversion increases at a faster rate near the wall compared to the conversion at the center. The increase in the temperature near the wall is, however, much less to offset the viscosity rise due to higher conversion. Therefore, there is more resistance to the flow near the wall, and the material at the center has much higher axial velocity than the corresponding Newtonian flow of constant viscosity. However, at the flow front, the axial velocity at the center reduces to the average velocity of the flow front. Since the fluid element at the contact line is allowed to slip, the velocity at that point is not equal to 0. The corresponding flow front for Run #3 is compared in Fig. 6 with the flow front for Run #1, in which the change in viscosity is marginal. The above effect of "surging" of material near the center is also apparent from the shape of the flow front for Run #3.


We have developed a model of the advancing flow front in RIM. No assumptions were made to separate flow front region, regarding the shape of the front and velocities, conversion, and temperature distribution in the front. The accuracy of the flow front region is considerably improved by using the Petrov-Galerkin formulation, rather than the Galerkin formulation. The model is general, and can be used for various RIM systems. However, the results presented here are for well-characterized polyurethane systems with available experimental data.

The predictions of the model for the pressure rise are in excellent agreement with the experimental data, even close to the gel point. The predictions of the conversion and temperature profiles near the flow front are considerably affected by the type of assumptions made in the calculations of this region. The shape of the flow front itself also depends on the dynamics of all these variables and their effect on viscosity and kinetics, and vice versa. Although these computations are more time-consuming than the ones with simplified assumptions, they provide more refined predictions. These predictions are better suited in estimating fiber orientation and bubble growth in final RIM parts, in which the flow front region plays the most important role. The issues related to fiber orientation and bubble growth are being investigated in this ongoing model development with ultimate goals to improve mechanical properties and to reduce scrap.


1. J. M. Castro and C. W. Macosko, AIChE J., 28, 250 (1980).

2. 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 Meeting, New York (1987).

3. M. A. Garcia, C. W. Macosko, S. Subbiah, and S. I. Guceri, Inem. Polym. Process., 6, 73 (1991).

4. C. D. Lack and C. A. Silebi, Polym. Eng. Sci., 28, 434 (1988).

5. R. E. Hayes, H. H. Dannelongue, and P. A. Tanguy, Polym. Eng. Sci., 31, 842 (1991).

6. T. J. Hsu and L. J. Lee, Polym. Eng. Sci., 28, 955 (1988).

7. K. J. Wang, Y. J. Huang, and L. J. Lee, Polym. Eng. Sci., 30, 654 (1990).

8. T. C. Papanastasiou, N. Malamataris, and K. Ellwood, Int. J. Numer. Meth. Fluids, 14, 587 (1992).

9. N. R. Anturkar, submitted for publication in Compt. Meth. Appld. Mech. Eng. (1993).

10. D. J. Coyle, J. W. Blake, and C. W. Macosko, AIChE J., 33, 1168 (1987).

@ @
COPYRIGHT 1994 Society of Plastics Engineers, Inc.
No portion of this article can be reproduced without the express written permission from the copyright holder.
Copyright 1994 Gale, Cengage Learning. All rights reserved.

Article Details
Printer friendly Cite/link Email Feedback
Title Annotation:reaction injection molding
Author:Anturkar, Nitin R.
Publication:Polymer Engineering and Science
Date:Oct 15, 1994
Previous Article:The onset of wall slip and sharkskin melt fracture in capillary flow.
Next Article:Shear and elongational behavior of linear low-density and low-density polyethylene blends from capillary rheometry.

Related Articles
Three-dimensional modeling of reaction injection molding.
Reaction injection molding process.
Reaction injection molding: analyzing the filling stage of a complex product with a highly viscous thermoset.
Bubble growth in reaction injection molded parts foamed by ultrasonic excitation.
Reaction Injection Molding Process of Glass Fiber Reinforced Polyurethane Composites.
A 3-D Finite Element Model for Gas-Assisted Injection Molding: Simulations and Experiments.
Metal injection molding: 3D modeling of nonisothermal filling.
Simulation of injection-compression molding for optical media.
Three-dimensional simulation of primary and secondary penetration in a clip-shaped square tube during a gas-assisted injection molding process.

Terms of use | Copyright © 2018 Farlex, Inc. | Feedback | For webmasters