A numerical model for solving polymer melt flow

 

 

L. A. le Grange^I; G. P. Greyvenstein^II; W. J. de Kock^II; J. P. Meyer^II

^ISenior lecturer, Department of Mechanical Engineering, Potchefstroom University for CHE, Potchefstroom 2520, South Africa
^IIProfessor, Department of Mechanical Engineering, Potchefstroom University for CHE, Potchefstroom 2520, South Africa

 

 

---

ABSTRACT

A finite volume numerical method is presented for the prediction of unsteady, incompressible polymer melt flow. Most other researchers use finite element methods in numerical models simulating polymer melt flow. The continuity, momentum and energy equations are considered in terms of general curvilinear coordinates. A generalised Newtonian model is used to model the non-Newtonian viscosity which, in the present case, is a function of the deformation of the flow field and of temperature. The method may be used to obtain valuable information about the flow properties during the filling stage of the injection-moulding process. The method is tested by comparing the pressure drop versus volumetric flow rate for a case of the filling of a thin mould cavity with available computational data.

---

 

 

Nomenclature

Latin letters

Acoefficient, area

Asurface vector

b, c, dmetric coefficients

C specific heat, remaining terms

einternal energy

Eenergy

f any function

gbody forces

Gcontravariant velocity component

I flux over control volume face

J jacobian

kconductivity

ppressure

Ppressure; grid point control functions

q heat loss by conduction

Ssource term

ttime

Ttruncation error; temperature

u, v Cartesian velocity components

V volume

Vvelocity vector

X, y, zCartesian coordinate directions

Greek letters

α under-relaxation parameter

Γdiffusion coefficient

δdifference; Kronecker's delta

∆increment

∂derivative operator

derivative operator vector

ηcurvilinear coordinate; non-Newtonian viscosity

Π stress tensor

ẏrate-of-deformation tensor

κdilatational viscosity

μ.dynamic viscosity

τ shear stress tensor

ρdensity

Φscalar quantity

φ,ψpolynominal functions

ξcurvilinear coordinate

ζcurvilinear coordinate

Subscripts

einternal energy; centre of cell face: e

e, w, n, s, t, b centre of cell face: e, w, n, s, t, b

mreference to momentum

nbthe neighbours: E, W, N, S, T, b

p constant pressure

Φreference to a scalar quantity

ttotal

v constant volume

X, yCartesian coordinates

ξ, η, ζ derivatives with respect to curvilinear coordinates

1,2,3reference to x, y, z coordinate directions

Supercripts

cconvection

ddiffusion

dccross-derivative diffusion term

dnnormal-derivative diffusion term

Ttranspose

uCartesian velocity component

Φreference to a scalar quantity

* preliminary value

 

Introduction

Injection moulding is one of the major basic converting processes to manufacture plastic products. In this process molten plastic is injected into a mould after which it cools down rapidly. When the melt has set adequately, the completed part is ejected from the mould. During the past two decades the use of polymer materials or plastics has increased dramatically compared to the more traditionally used materials such as wood and metal. It is, however, more difficult to control dimensions and strength properties of plastic products with injection moulding. During the process the polymer melt cools down, which causes shrinkage of the product. The shrinkage problem provides some difficulty in controlling the dimensions of the component whereas the strength of moulded articles depends on flow properties during the injection moulding process. Some of these properties are the melt temperature, the volumetric flow rate, the injection pressure and cooling time.

The design of parts and moulds if often a trial and error process, backed by experience. The iterative nature makes the design process, up to the production of the first successful part, very expensive. The huge costs associated with the design process can be cut if the number of iterations is reduced. This can be done by obtaining a better understanding of the process of polymer melt flow.

As computers have become more powerful, the possibility of numerical stimulation of the injection-moulding process has emerged. Simulation models will not only enable the design engineer to design better moulds more cost-effectively but also to reduce the number of mould changes with the implication that lead times can be reduced. The value of computer simulation for analysing the injection-moulding process has previously been demonstrated [1, 2] in the design of numerous injection moulds. The development of this technology has made a major contribution to the processing of plastics. Although the usefulness of simulation models is clear, problems still exist which limit their application. These problems can mainly be attributed to inadequate materials modelling, complex mathematical models and the geometrical complexity of flow domains.

Approximations to simplify the complex mathematical models have been reported in the literature. The generalization of classical lubricatoin theory of non-isothermal, non-Newtonian flow, called the generalized Hele-Shaw formulation, has, according to Güceri, [3] become the standard way to formulate injection mould filling problems. Simulation models based on the Hele-Shaw formulation solve for one variable, i.e. pressure (p), only. The two velocity components (u, v) are then derived from the pressure.

Richardson [4] was the first to suggest a solution based on Hele-Shaw type of flow. Since 1975 various other detailed theoretical studies of two-dimensional flow in filling thin cavities, based on Hele-Shaw type of flow, have been reported as reviewed briefly by Wang, Hieber and Wang. [5]

The generalized Hele-Shaw model introduced by Hieber and Shen [6] provides simplified governing equations for non-isothermal, non-Newtonian and inelastic creeping flows in thin cavities. Besides the usual lubrication approximations, the velocity component in the thickness or gap direction is also neglected. Therefore the pressure is a function of only the other two spatial variables. In Hele-Shaw motion the ratio of inertia to viscous forces, just as in the case of the lubrication approximation, is less than one. If it exceeds unity the inertia forces become considerable and the motion deviates from the simple solution. They employ a hybrid numerical scheme in which the planar coordinates are described in terms of finite elements (triangular elements and quadratic shape functions are employed for p) and the gapwise and time derivatives are expressed in terms of finite differences. The numerical solution of the field equation for pressure is carried out by successive under-relaxation. If we compare this hybrid formulation with preceding work it is much more suitable for handling cavities of more general planar geometry.

Wang and Hieber [7] report on advances to simulate the cavity-filling stage in terms of the melt viscosity of an incompressible fluid using thin two-dimensional finite elements. They describe the filling stage as creeping, shear-dominated flow which might also be treated in terms of classical Hele-Shaw flow as generalised to a non-Newtonian, non-isothermal, non-steady situation. They developed the method of Hieber and Shen [6] further by employing a fixed finited element grid using a control volume approach which makes the melt-front advancement amenable to an automated algorithm. In addition, this development has been extended to handle thin, three-dimensional parts, treated as a union of the two-dimensional components.

In this paper a finite volume numerical method is presented for the prediction of unsteady, incompressible polymer melt flow as compared to finite element methods used by most other researchers simulating polymer melt flow. The governing equations are considered in terms of general curvilinear coordinates in order to handle complex geometries. With regard to the problems of inadequate materials modelling and complex mathematical equations, very elementary non-Newtonian flows are considered. In the subsequent section the equations governing non-Newtonian fluids are given.

 

Theoretical formulation

The motion of any fluid is described by the equations for the conservation of mass, momentum and energy. These equations can be expressed in the following vector form, which is independent of coordinate systems:

 

Newtonian fluid dynamics

Newtonian fluids behave according to Newton's law of viscosity. This law states that the shear force per unit area for simple Newtonian shear flow is proportional to the negative of the local velocity gradient, or

The appropriate generalisation for arbitrary, time-dependent flows is [8]

 

Non-Newtonian fluid dynamics

Non-Newtonian fluids are structurally complex. Because of their complicated nature, they behave quite differently from Newtonian fluids like water and gases. For Newtonian fluids the relation between the stress and the deformation of the velocity field is described by the constant Newtonian viscosity μ. In the case of non-Newtonian fluids, this relation is described by a more complicated expression. In general, the expression for non-Newtonian fluids is non-linear in velocity gradients and contains information on the time-dependent mechanical response of the fluid.

A very restricted type of non-Newtonian fluid is considered in this study. The flow is restricted to incompressible, steady shear flow and any time-dependent behaviour is omitted. This fluid type is described by the generalised Newtonian model. The generalised Newtonian model uses a non-Newtonian viscosity η, analogous to the Newtonian viscosity μ.

The non-Newtonian viscosity is a function of the deformation of the flow field and may also be a function of temperature and pressure.

Many empiricisms for the viscosity η (ẏ) are available to describe the behaviour of different non-Newtonian fluids. A few simple and popular models are the Power Law, the Eyring Model, the Ellis Model, and the Bingham Fluid Model. [9] These models are limited to steady-state shearing flows, and are thus in general inappropriate for the description of unsteady flows where the elastic response of polymeric fluids becomes important.

For non-Newtonian fluids, using the generalised Newtonian model, the momentum equation (2), can thus be expressed as

and the energy equation is expressed as

where S_m and S_e cointain source and residual terms.

 

Transformation of the governing equations

The dependent variables in the momentum and energy equations (equations (7) and (8)), all obey a generalised conservational principle. This enables the formulation of a generalised differential equation

where represents any one of the dependent variables, u, v, w or T; Γis the diffusion coefficient, and S^Φ a source term. The quantities Γand S^Φ are specific to a particular meaning of y. For example, in the momentum equation and energy equation, may include pressure terms and viscous heating terms respectively.

A transformation from Cartesian coordinates to general curvilinear coordinates is accomplished by using certain concepts from differential geometry and tensor analysis. The transformed equation is

The resulting equation (10) is of the same type as the original general differential equation (9), but is more complicated in that it contains more terms and variable coefficients. The domain, on the other hand, is greatly simplified since it is transformed to a rectangular region regardless of its shape in the physical space. This permits standard discrete representations of derivatives along the transformed boundary, without the need for interpolation.

 

Discretisation

In a finite volume approach equation (10) is integrated over a finite volume or cell. A typical cell, its centre point and the centre points of its neighbours are shown in figure 1. By using Gauss' theorem, [10] a volume integral can be converted into surface integrals over the six faces of a cell, thus

 

 

where, for example, f_e is the component of f normal to the "e" face of the cell and A_e is the area of that cell face. Equation (11), integrated over a control volume, can now be written in general form as

where each term comprises two distinct parts: - convection I_c and diffusion I_d, eg.

with the normal derivative term (diffusive) to the "e" face and the cross-derivative terms. The derivatives of Φ are approximated by

The values of Φ at cell faces are evaluated by linear interpolation in the physical domain.

After substituting expressions like (13), (14) and (15) into equation (11), a relation between Φ and the neighbouring values is obtained, i.e.

where A_nb is the neighbouring coefficients which account for the area, and the combined convection-diffusion influences at the control-volume faces, and Φ_nb is the neighbouring Φ-values.

 

Numerical solution method

An equation like (10) exists for every interior cell in the grid. The numerical solution procedure therefore entails solving a system of M equations in M unknowns, where M is the total number of cells. Such a system of equations can be written in matrix form as

where [A] is the coefficient matrix consisting of M × M coefficients is a column matrix having M entries, and {C} is also a column matrix with M entries.

Equation (19) can be solved in several ways. [11] Firstly, direct solution methods or matrix inversion offers exact solutions for the given coefficients. These methods are not well suited because they require storage of a large number of coefficients which can easily exceed the memory capacity of available computers. The second and more appropriate way is using iterative solution methods. The use of these methods requires little additional storage capacity and offers faster computing times. The Thomas algorithm, [12] a line iteration procedure, has been adopted as the iterative solver. This method sweeps through the computational grid line-by-line.

Iterative methods start the solution process from a set of guessed values for the dependent variables. The method then updates each variable and coefficients in the grid successively until a converged solution is obtained. When highly nonlinear equations, as in the present case, are being solved, iterative methods tend to become unstable and the solution can diverge. It is then necessary to slow down changes in Φ by under-relaxation. This can be done by writing equation (19) as,

where is the value of Φat node P, obtained in the previous iteration. The term in brackets represents the change in Φ_pfrom one iteration to the next. The change in α_Φ is reduced by an under-relaxation parameter α_Φ, in the range from zero to unity. The solution method converges when the term in brackets tends to zero.

 

Pressure calculation procedure

A difficulty exists in calculating the velocity field using the momentum equation. This is due to the fact that the pressure field is unknown. Since no obvious equation for obtaining the pressure exists, a pressure calculating procedure, based on the SIMPLE algorithm, [11] is used. The method produces the solution in primitive variable form, and is fully capable of solving recirculating flow problems. The main difference between the current method and the SIMPLE method is the use of a non-staggered grid, as opposed to the staggered grid arrangement used in the SIMPLE method. It is well-known that velocity-pressure decoupling arises with the use of a non-staggered variable arrangement. The technique of Rhie and Chow [13] is used to prevent the decoupling problems.

 

Boundary conditions

The boundary conditions which are used in the present study are those of inlet, outlet, symmetry plane and fixed wall boundaries. Velocities and temperatures are normally specified on the inlet and fixed wall boundaries. At the outlet and symmetry boundaries, zero-gradients for velocities and temperatures are assumed. As a result of the non-staggered variable arrangement on the grid, the pressures on the boundaries are needed in solving the momentum equations. At the walls, pressure values are calculated by assuming that the normal gradients to the walls are zero. The pressure values at inlet and outlet boundaries are obtained by extrapolating from the interior.

 

Overall solution algorithm

The sequence of operations for a complete iteration may be as follows: The field values for all dependent variables are initialised using a realistic guess; the momentum equations are solved to obtain preliminary u, v and w velocities; the pressures and velocities are updated using the SIMPLE algorithm; the energy equation is solved to obtain temperatures; the non-Newtonian viscosity is solved from the prevailing velocity field; and return to the first step for a new iteration cycle using the prevailing velocities, pressures, temperatures, and viscosities. The sequence is repeated until a final converged solution is reached.

 

Verification

To verify the numerical solution method, the filling of a thin cavity with polypropylene during the injection moulding process is considered. Sanou et al. [2] predicted the filling process with a simulation model and then performed it experimentally. The main difference between the numerical model used by Sanou et al. [2] and the current model, is that the current model can handle three-dimensional elliptic flows in complex geometries. The model used by Sanou et al. [2] is only suited to handle two-dimensional flows in thin rectangular cavities. Sanou et al. [2] published their results in the form of pressure drop versus volumetric flow rate, figure 2 shows a diagram of the cavity and the positions of the pressure transducers.

 

 

Sanou et al. [2] reports a nearly symmetrical flow at the position of the pressure transducers. Therefore a two-dimensional analysis was performed, using a 15 × 21 mesh. The grid lines increased in an exponential manner from the centre-line to the walls.

A series of simulations were done for different volumetric flow rates. Uniform inlet velocities in the range of 0,13 to 1,0 m/s were specified at the inlet. At the outflow boundary zero-gradients were imposed on the velocity components. A wall temperature of 30 °C and a melt temperature of 200 °C were specified.

To describe the viscosity in their numerical model Sanou et al.² used an asymptotic "power-law" function

where b = 0,0096 kg/m.s, n = 0.35,m T_b = 5900 K and C = 480 m/N. The material properties used for polypropylene were density ρ = 700 kg/m³, specific heat C_p = 3462 J/kgK and the conductivity

k = 1,507 × 10^-2 W/m°C.

The dependence of the viscosity on temperature is exponential, as can be seen from the current viscosity model in equation (21). This requires an unsteady analysis, specifically for the development of the temperature field. The numerical results of the pressure drop between the two flush-mounted transducers at points p1 and p2 (shown in figure 2), versus volumetric flow rate, are presented in figure 3.

 

 

It can be seen from figure 3 that the highest pressure drop occurs at low volumetric flow rates. At low volumetric flow rates the flowing melt is exposed to the cold mould walls for longer periods, which causes a higher heat loss to the walls. The low temperatures are responsible for high viscosities η, as can be seen from equation 21. The high viscosities are responsible for high shear stresses which cause high pressure drops. With higher flow rates there is less time for heat loss to the walls, which results in lower pressure drops. At the highest flow rate the pressure drop increases due to the high shear stresses caused by the high volumetric flow rate.

The pressure drops, according to the current simulation model, are close to the numerical results of Sanou et al. The present results are, as were the case with the reported pressure drops, approximately 30 percent higher than the experimental results of Sanou et al. These authors report that the discrepency between their numerical and experimental results is not easy to interpret. A few possible causes were, however, identified during the simulation process, using the current model. The first of these is the manner in which the filling process was simulated. In the current simulations the advancement of the melt front was not modelled. Instead, fully-developed flow through an open-ended, thin cavity was considered.

Unsteady heat flow from the melt to the walls was modelled for a time equal to that of the filling stage. The filling stage takes 0,29 seconds in the case of the lowest flow rate and 0,07 seconds for the highest flow rate. Variations in the velocity and temperature fields were calculated during these time increments, using hundred time-steps for each analysis. This method, however, allows for more heat transfer, which results in higher pressure drops, than when the advancing melt front would be taken into consideration.

The second cause of the high pressure drops is the inaccurate materials modelling. This includes the use of an incorrect value for the thermal conductivity k. A constant value for the conductivity was used in the present simulations. It was also found that the pressure drop is extremely sensitive to changes in conductivity. During the filling phase the warm melt near the walls solidifies due to high heat loss. The thermal conductivity in these layers near the walls is smaller than nearer to the core of the flow. This has the effect that less heat transer occurs from the hot flow core which will cause lower pressure drops. On the other hand, the higher shear stresses near the walls cause higher pressure drops. The modelling of the heat transfer during the filling stage is therefore an area which needs further research. Another very important inaccuracy in the current materials modelling, is the constant density assumed for polymer melts, which are in actual fact compressible. Due to this compressibility the time required to fill a closed cavity and the pressure drop along the cavity are affected in the experiments. Furthermore the semi-empirial power-law does not take the important pressure dependence of viscosity into account. In order to compare the results of the current model with those reported by Sanou et al.,² the same materials modelling had to be employed. However, to compare the results of the numerical model with experimental results more accurate, materials modelling will be required.

Although the current results were in satifactory agreement with the numerical results, the solution time of the simulation was unsatisfactory. This was due to severe under-relaxation used in updating the viscosity after each iteration. The viscosity decreases exponentially with temperature as can be seen from equation (21), which results in huge variations in viscosity during the cooling of the melt. The large discrepency between the numerical and experimental results, due mainly to inaccurate materials modelling, emphasizes the need for further work on the materials side before the results of the numerical model can be expected to correlate better with practice.

 

Conclusions

Knowledge of polymeric melt flow in injection mould cavities is of importance to both the designers of both the plastic component and the injection mould. In this study a finite volume method using a boundary-fitted coordinate system, to solve polymeric melt flow is presented. The method is based on the SIMPLE algorithm and makes use of the power-law formulation of shear stress and strain rate.

The method is verified by comparing predictions of pressure drop through a cavity with numerical results reported in the literature. The results illustrated the usefulness of the current simulation model as a basis for further work.

 

References

1. M. R. Kamal and P. G. Lafleur, "Computer simulation of injection moulding", Polymer engineering and science, 22 (27), 1066-1074, 1982.         [ Links ]

2. M. Sanou, B. Chung and C. Cohen, "Glass fiber filled thermoplastics II. Cavity filling and fiber orientation in injection moulding", Polymer engineering and science, 25, 1008-1016, 1985.         [ Links ]

3. S. O. Güceri, "Finite difference solution of Field Problems", Chapter 5 (In"Fundamentals of computer modelling for polymer processing", Ed:, C. L. Tucker III, Hanser Publishers, 1989.)

4. S. Richardson, "Hele-Shaw flows with a free boundary produced by the injection of fluid into a narrow channel", Journal of Fluid Mechanics, 56 (4), 609-618, 1972.         [ Links ]

5. V. W. Wang, C. A. Hieber and K. K. Wang, "Dynamic Simulation and Graphics for the Injection Moulding of Three-Dimensional Thin Parts", Polymer Engineering, 7 (1), 21-45, 1986.         [ Links ]

6. C. A. Hieber and S. F. Shen, "A Finite-element/Finite-difference Simulation of the Injection Moulding Process", Journal of Non-Newtonian Fluid Mechanics, 7 (1), 1-32, 1980.         [ Links ]

7. K. K. Wang and C. A. Hieber, "Injection Moulding Simulation", Proceedings of Manufacturing International, 87-95, 1988.

8. R. B. Bird, W. E. Stewart and E. N. Lightfoot, Transport phenomena, John Wiley, New York, 1960.

9. R. B. Bird, R. C. Armstrong and O. Hassager, Dynamics of polymeric liquids, Vol. 1, Wiley, New York, 1977.

10. E. Kreyszig, Advanced Engineering Mathematics, 4th Edn, John Wiley, New York, 1979.

11. S. V. Patankar, Numerical heat transfer and fluid flow, Hemisphere, Washington, 1980.

12. D. A. Anderson, J. C. Tannehill and R. H. Pletcher, Computational fluid mechanics and heat transfer, McGraw-Hill, Washington, 1984.

13. CM. Rhie and W. L. Chow, "Numerical study of the turbulent flow past an airfoil with leading edge separation", AIAA Journal, 21, 1525-1532, 1983.         [ Links ]

 

 

First received January 1992
Final version December 1992