Printer Friendly

Effect of extensional viscosity and wall quenching on modeling of mold fillings.


The use of numerical methods to simulate the filling of molding dies is important to the selection of die, runner, and gate geometries, process conditions for successful moldability, and to the production of high-quality molded parts. The results of modeling the energy and mechanical transport within the die often come in the form of temperature, hydrostatic and deviatoric stress, velocity, and phase distributions over the course of filling. These predictions of filling behavior are used to anticipate such problems as jetting, premature freeze-off, warpage, weld-line formation, variations in filler concentration, and cooling fracture. Results from numerical modeling are also used to anticipate process conditions such as filling pressure/stress requirements, die wall temperatures, and filling and packing schedules.

Numerical simulations often apply simplifying assumptions regarding the governing equations and molding material properties in order to make the analysis tractable. The realism of such simulations breaks down, then, when the analysis is applied to an inappropriate geometry or an inappropriate law or fit was applied to the material data. Molding material constitutive behavior is inseparably linked to its response to the prescribed process conditions and, thus, to the proper prediction of filling behavior.

Extensional deformation occurs when velocity gradients exist in the streamwise direction. It is the only deformation that can occur in irrotational fluid flow. During the filling of molding dies, extensional deformation occurs as the melt is convected through the cavity and must accelerate in the direction of its motion. This occurs, for example, when there are changes of flow cross section (e.g., at a gate, tapered wall, or thick/thin juncture). Because momentum is dissipated by wall friction as melt moves parallel to the die wall, there is shearing deformation as well. For those in engineering practice who are responsible for making realistic die filling simulations, it has been recognized that convergent and divergent passages in the die require careful treatment because they have distinctly different flows (1).

It has been common for filling codes that are capable of simulating the filling of complicated cavities to have viscosity models based entirely on the molding material's response to viscometric flow. There are many filling simulations that have shear-rate-dependent constitutive relations (2-6). The exclusive use of such constitutive models has the potential to limit a given code's worthy application for at least two reasons. First, the analysis of most injection molding needs to treat both shear and extensional deformation at the same time. Virtually any combination of die, material, and process conditions will result in material deformations that are not viscometric but combinations of shear and extension in significant portions of the die. The only common exception to this would be the molding of thin-walled plates and shells by materials that remain entirely fluid during filling. Second, there are many candidate molding materials whose constitutive behavior under extensional deformation cannot be predicted from their responses to shearing deformation. While the measurement of elongational viscosity often poses a significant experimental challenge (7-9) its presence bears upon filling behavior.

In terms of materials that are significant candidates for injection molding, a preponderance of evidence demonstrates that some of them possess differing constitutive behaviour under shear than under extension. Powder injection molding feedstocks have been shown to be far more resistive to flow through convergent capillary dies than would have been predicted from their shear viscosities (10). Dilute suspensions of fibers in fluid matrices have been predicted to have extensional viscosities that are an order of magnitude greater than the shear viscosity of the matrix (11). This has been confirmed experimentally (12-14). Some of the polymers and polymer solutions showing such dualistic behavior are polyacrylamide solutions (15), polystyrene melts (16), and other solutions (17-19). The difference between shear and elongational viscosity has been discussed in terms of internal fracture in highly loaded powder injection molding melts (10), unwinding and aligning of molecular chains in fluid polymers flows (20, 21), and in terms of interaction between stress and velocity fields in fiber-loaded flows (11).

Die filling simulations that utilize only shear-dependent viscosities, then, are not necessarily appropriate for predicting die filling behavior by materials having such dualistic behavior. However, there is a dearth of computational models that allow for the incorporation of distinct sensitivities of the viscosity to shear and extensional deformation rates. Barrie (22) included a pressure term for tensile viscosity effects in spreading disk flow. Kwon, et al. (23), simulated normal stress effects in conical die flows.

The inclusion of an extensional viscosity into an analysis often requires the commitment of large computational resources. The models of Schunk and Scriven (20) and of Evans and Heyes (24) were both run in supercomputer environments. Schunk and Scriven modeled that intake flow between opposed nozzles for a fluid characterized by two Carreau-like constitutive equations. This model has been applied to the modeling of the film coating process. Evans and Heyes have used a statistical mechanics model to predict the shear and elongational viscosity of monoatomic fluids. The intent of this was to prescribe constitutive laws for material continuum behavior. Dinh and Armstrong (25) analyzed the startup of steady shear and elongational flows of fiber suspensions individually to predict fiber orientation and bulk viscosity.

The approach described here is the inclusion of a dual viscosity model into a finite element formulation for quasi two-dimensional die filling. The viscosity model becomes the vehicle for distinguishing between elements that are parallel-walled flow passages (and therefore contain only shearing flow) and those that are convergent or divergent (and therefore contain a combination of shearing and extensional flow). The numerical analysis described here was run on a workstation in several hours for each modeling case.

The viscosity of the molding material is treated on a bulk level. It reduces to commonly used constituttve laws for viscometric and shear-free flows. Information on molecular or fiber orientation is not part of the results of this model. The inclusion of a configuration distribution function into a model greatly increases mathematical complexity (26). Thus far, such predictions have been done only for simple flow fields or materials (24, 26, 27) or for a specific flow and material at great computational effort (20).

The analysis presented here involves estimation of all of the components of the rate-of-deformation tensor in each of the elements over the course of filling and computation of the corresponding dual viscosities and flow conductivities. Freezing of melt at the die walls and the corresponding narrowing of the flow passage is predicted. Although phase segregation in a multiphase molding material is not predicted directly, the stress gradients predicted can be indicators of likely separation in such molding material systems. The quantities calculated (especially flow fronts, hydrostatic pressure, temperature, and pressure and thermal gradients) are often related to molded part quality and can be applied as initial conditions to cooling and packing analyses, as well.


The analysis will be considered first in terms of the elemental flow field and the corresponding rate-of-deformation tensor, then a viscosity model will be applied and examined, and finally the momentum equation will be reduced to an elliptic form, which can be handled by a finite element formulation. The effect of wall freezing appears as a temperature-dependent specific heat in the energy equation and as a reduced flow passage width in the momentum equation.

Flow Field

The elemental flow field is treated as a quasi two-dimensional flow between nearly parallel walls and is examined over a triangular region. The flow field is considered here to be a superposition of Hele-Shaw and biaxial extensional flows. A generic example of the flow passage is shown in Fig. 1. The velocity field is symmetric about the x - y plane. The width of the flow passage is defined by a half-gap height, b(x, y), which varies linearly over the triangular region. Let the distance from the x - y plane to the walls be given by

b = b(x, y) = [Alpha]x + [Beta]y + c (1)

for x and y within the triangular region.

Further, if hydrostatic pressure is defined at each vertex of the triangular region, then pressure gradient may be approximated from the relation

[nabla]P [approximately equal to] ([P.sub.b] - [P.sub.a]/[x.sub.b], 1/[y.sub.c] ([P.sub.c] - [P.sub.a] + [x.sub.c]/[x.sub.b] ([P.sub.a] - [P.sub.b]))) (2)

It is assumed (for the moment) that hydrostatic pressure gradient is the predominant motivator of the flow. The elemental flux vector is then linearly related to the hydrostatic pressure gradient by the flow conductivity and half-gap height.

q = - S[nabla]P (3)

where mean velocity is related to the flux vector by

[Mathematical Expression Omitted]

and where flow conductivity is the same as defined by Hieber and Shen (6)

[Mathematical Expression Omitted]

Rate-of-Deformation Tensor

The velocity gradients within the passage are approximated from the quantities available, namely, the elemental mean velocities and the elemental flow passage geometry. The passage is here defined as mildly convergent when

[[[nabla]b]] [much less than] 1 (6)

By referring to the flow field as quasi two-dimensional, it is meant that u [much greater than] w and v [much greater than] w. By combining this velocity scaling with the assumption that the characteristic dimension in the z direction, b, is much less than the base and height of the element, it is concluded

du/dz [much greater than] dw/dx (7)

dv/dz [much greater than] dw/dy (8)

The gradient in the half-gap height ([nabla]b(x, y)) and the continuity equation are used to estimate the components of the velocity gradient tensor that arise from extensional deformation. The component of the gap gradient vector that is parallel to the velocity vector results in extensional deformation. The component that is normal to the velocity vector results in shearing deformation. At any x and y in the element, it is assumed that the integral of the velocity profile in the z direction is a constant. This is equivalent to the flux vector being constant over the element. The extensional portion of the velocity gradient tensor becomes the dyad product of the velocity and the thickness gradient.

[nabla][V.sub.ext] = V [cross product] [nabla][b.sup.T] (9)

Since the predominant components of the vorticity vector are in the x and y directions, the shear portion of velocity gradient tensor contains two terms, du/dz and dv/dz. The rate-of-deformation tensor then becomes

[Mathematical Expression Omitted]

Viscosity Model

A practical fluid model to apply to the filling of complicated die shapes by dualistic fluids would have several characteristics. It would be a function of in-variants of the rate-of-deformation tensor. It would, within the assumptions made for the flow field, have different sensitivities to extension than to shear rate, and bear resemblance to existing constitutive laws for viscometric and shearfree flows. Material viscosity models that are functions of the invariants of the rate-of-deformation tensor are among the simplest to implement because information of local-to-global coordinate orientation and material anisotropy need not be included in the analysis. The model used here is given by

[Tau] = 2[[Chi].sup.D] (11)


[Chi] = m(T)[(-4[II.sub.D]).sup.(n - 1)/2][(1 + q[absolute value of [[Lambda].sub.ext]]).sup.p] (12)


[II.sub.D] = 1/2 ([(trD).sup.2] - tr[(D).sup.2]) (13)

The quantity [[Lambda].sub.ext] is one of the eigenvalues of the rate-of-deformation tensor at the plane of the element (i.e., for z = 0). It is entirely a function of extension rate and is invariant under coordinate transformation.

[Mathematical Expression Omitted]

For the type of flow field modeled, the selected model becomes functionally different between the parallel-walled and convergent-walled cases. When the fluid is between parallel walls, [nabla]b = 0 and the constitutive law reduces to the often used simple power law relation (28, 29). When [nabla]b is nonzero, the flow passage is convergent or divergent, extensional or compressional deformation exists. The flow in the passage is considered convergent when

[Mathematical Expression Omitted]

in which case

p = c where c [greater than] 0 (16)

The flow is considered divergent when

[Mathematical Expression Omitted]

in which case

p = -d where d [greater than] 0 (18)

Many molding materials exhibit shear thinning and extension thickening (30). When viscometric flow is applied to this model, that is,

[Mathematical Expression Omitted]

the viscosity reduces to the power law

[Mathematical Expression Omitted]

When biaxial stretching flow is applied to this model, that is,

[Mathematical Expression Omitted]

the Trouton viscosity is

[Mathematical Expression Omitted]

and when [Mathematical Expression Omitted],

[Mathematical Expression Omitted]

For the flows considered here, the material is shear thinning when n [less than] 1 and extension thickening when

c [greater than] n - 1 (24)

Solidification at the die walls during filling is also part of the numerical analysis. Viscosity is then defined only for temperatures above the no-flow temperature. Viscosity was selected to vary inversely with the difference between the temperature and the no-flow temperature so that viscosity would rise rapidly near the temperature for which the material would be considered a solid. This is a two-constant empirical relation.

m(T) = A/T(z) - [] + B (25)

Comparison With Constitutive Laws

The model and its material constants can be selected such that it begins to replicate existing constitutive laws for purely extensional flows. Batchelor (11) has derived a constitutive law for the convergent flow of nondilute fiber suspensions in a Newtonian matrix in axisymmetric coordinates. The law is given by

[Mathematical Expression Omitted]

If biaxial stretching flow is prescribed in the elemental coordinate system and [Mathematical Expression Omitted], the models have the same extensional rate dependence when the exponent c is selected such that

c = -(n - 1) (27)

and the coefficient q is then selected such that

q = [([[Mu].sub.o](1 + 4/9 C[(AR).sup.2]/ln([Tau]/C)) / m [3/2.sup.(n - 1)/2]).sup.1/c] (28)

Application to Rheometric Data

The material constants q and c given in Eq 25 may be extracted in an approximate manner from rheometric data taken from conical dies. Both Cogswell (7) and Dealy (8) have pointed out the ambiguities existing in the interpretation of data from convergent flows in which the types of deformation are mixed and nonuniform. While measurement of viscosity in viscometric flow is relatively common, measurement of viscosity in a purely extensional flow field presents an extremely challenging experimental problem for many materials. It has been done relatively successfully for some materials in carefully controlled experiments (31, 32). In the absence of data from such experiments, however, the modeler may have to resort to more approximate results gotten by simpler methods in order to have the inputs necessary to predict the filling of a die.

In the work of Rosner, et at. (10), a highly loaded powder injection molding feedstock at molding temperature was passed through a set of convergent dies of differing cone angle. Apparent viscosity as a function of apparent shear rate is shown in Fig. 2a. In their presentation of the data, there is no estimation or accounting for the presence of extensional deformation rate in the conical dies, however.

If apparent shear rate is a slightly tapered conical die where averaged over its length, the result is

[Mathematical Expression Omitted]

The Rabinowitsch correction (33) converts this to a wall apparent shear rate.

[Mathematical Expression Omitted]

where lengthwise-average wall shear stress is taken to be

[[Tau].sub.w] = [Delta]P/2 [R.sub.L] - [R.sub.O] / L 1/ln ([R.sub.L]/[R.sub.O]) (31)

[Mathematical Expression Omitted]

For a power law material, viscosity data from straight and slightly tapered capillaries would fall along curves within a few percent of each other. However, in the results of Rosner et al., apparent viscosities for increasing cone angle differed by more than this. Generally, apparent viscosities decreased with increasing apparent shear rate. For the treatment of these data as presented here, additional reductions were applied.

The determination of the material constants q and c is done as follows. The die geometry and volume flow rate are used to compute a length-averaged apparent extension rate.

[Mathematical Expression Omitted]

By comparing the rate-of-deformation tensors between the elemental flow field and the conical die flow field, it is found that

[Mathematical Expression Omitted]

The data from the straight capillary are fit to a power law relation. Viscosity ratios are computed by dividing the apparent viscosity at the given cone angle by the power law fit to the straight capillary data for the same apparent shear rate. If it is assumed (for the moment) the [II.sub.D] is dominated by shear rate, the ratio is then related to q and c by

[Mathematical Expression Omitted]

Viscosity ratio is plotted as a function of apparent extension rate and a curve is fitted through the data. The data and curve fit are shown in Fig. 2b. q was found to be 3.52 x [10.sup.-3] s and c was 4.0. The portion of curve fit shown in Fig. 2b begins at 50 [s.sup.-1]. This is the lowest extension rate occurring in the simulations of the filling of the ring die cavity that are described here. Agreement among the curve fit and the data is better above this value than below it. The data points in Fig. 2b appear to approach a value of two as apparent extension rate approaches zero, whereas the ratio approaches one in Eq 35. The source of the discrepancy among curve fit at data at a small extension rate is likely to arise from the behavior of the test material in the capillary entry region. It is known that this type of highly loaded mixture has a yield value at very small deformation rates (34).

Determination of Rate-of-Deformation Tensor

In order to compute the invariants of the rate-of-deformation tensor, estimates need to be made of its components. For laminar flow of a power law fluid between parallel plates, the components of [nabla]V may be described completely by the flux components, [q.sub.x] and [q.sub.y], the power-law exponent, n, and the mean half-gap height [Mathematical Expression Omitted]. For the slightly convergent flows considered here, the velocity profile will be taken to be described by the same two parameters. The rate-of-deformation tensor is then

[Mathematical Expression Omitted]

Shear and Extension Rate Interdependence

Dependencies upon extensional and shear deformation rates are not distinctly separated from each other in this model. The second invariant of the rate-of-deformation tensor is a function of both extension and shear rates. It would perhaps seem best to have used a model in which shear rate sensitivity could be adjusted completely independently of extension rate sensitivity and vice versa. Such a model would not have been invariant under coordinate transformation, however, and would be more computationally tedious to implement in a finite element or finite difference model.

The simulation described here lends itself best to cavity meshes for which shear rate predominates in all the elemental flows. This holds when gap heights are specified to give mildly convergent passages (e.g., [[[nabla]b]] [less than or equal to] 0.25). Such a requirement does not necessarily preclude the modeler from simulation of the filling of, for example, cavities having sudden contractions. Depending upon the molding material and geometry, the flow field through a contraction may in fact be a more gradual convergent flow surrounded by a recirculating flow or dead zone. For the mildly convergent flows treated here, the change in elemental flow conductivity is on the order of 10% or less when extension rate is included in the second invariant. An example illustrating the dependency of conductivity upon elemental thickness gradient is given by Moller (35).

Momentum Transport

The important task in handling the momentum equation is to reduce it to a form that can be operated upon by the finite element technique. Such a form has been described by Hieber and Shen (6). For some cases in which the elemental flow field is quasi two-dimensional and laminar, this form may be retained. The momentum equation for steady, laminar flow with negligible body forces is

0 = [nabla] [multiplied by] [Pi] = [nabla]P + [nabla] [multiplied by] [Tau] (37)

The momentum equation reduces to

[nabla]P [approximately equal to] ([Chi] [d.sup.2]V/d[z.sup.2] + d[Chi]/dz dV/dz) = d/dz ([Chi] dV/dz) (38)

The details of this derivation are given by Moller (35).

The momentum equation is in the same form as at the initiation of Hieber and Shen's derivation (6) of a finite element model of die filling. The steps following from here in the derivation of the governing equation for application of the finite element method will not be repeated here. Their governing equation becomes an elliptic equation in pressure.

[nabla] [multiplied by] (S[nabla]P) = 0 (39)

where the flow conductivity for the no-slip wall conditions becomes

S [approximately equal to] [integral of] [z.sup.2]/[Chi] dz between limits b and 0 (40)

The primary difference between this finite element approach and that of Hieber and Shen has to do with the fact that flow conductivity cannot be written entirely as a function of pressure gradient. The pressure solution for a given time step is instead determined as follows. A solution is found using the power law portion of the viscosity. Using these nodal pressures, elemental mean velocities are computed. The velocities are then used to compute elemental conductivities, which include extensional deformation rate sensitivity. Nodal pressures are then computed for these new (and typically lower) flow conductivities. These new nodal pressures and their gradients are used to compute the position for the flow front for the next filling step. This sort of approach can be applied only when volume flow rates, and not pressures, are specified as boundary conditions at the gate nodes.

Energy Transport

During filling, enthalpy is convected into each active element from upstream, the melt will conduct heat into the die walls, some melt may freeze at the die walls, and the energy to deform the melt will be dissipated in the form of heat.

The temperatures at discrete locations in the gapwise direction are computed for each filling step using a finite difference method operating on the energy equation. The encoding of this equation for the case of constant specific heat with dissipation due only to shearing deformation is described by Najmi and Lee (2). For constant material density and thermal conductivity and temperature-dependent specific heat, the equation is as follows.

[Rho] Dh/Dt = k[[nabla].sup.2]T + [Sigma] (41)

Convection is taken to be only in the streamwise direction and thermal conduction only in the gapwise direction.

The combined momentum/energy solution, then, involves the solutions at two sets of nodes. The first are those used for the finite element solution for pressure. They are referred to as primary nodes. The second set are nodes lying along lines orthogonal to the plane of the element and located at equally spaced intervals from the element plane to the half-gap height. Each line intersects the centroid of its respective element. These are used for the solution of the energy equation and are referred to as the gapwise nodes.

Dissipation due to both extension and shear is included in this analysis. The value of the dissipation function is computed from the components of the rate-of-deformation tensor for each element. Enthalpy of solidification is now included in this analysis, as well through the use of specific heat with a solidification peak. At each time increment, the enthalpy and temperature profiles across the element gap are computed recursively until sufficient convergence is reached.

Wall Quenching (Freezing)

The thermal and mechanical consequences of quenching the melt at the die walls are included in this filling model. Thermally, freezing is accommodated through the use of a specific heat, which is a function of temperature. Enthalpy is treated as a continuous piecewise linear function of temperature. A schedule of material enthalpy as a function of temperature is computed from a schedule of specific heat given as an input file. Specific heat will contain a peak in the vicinity of the no-flow temperature, the integral of this peak being the enthalpy of solidification.

The mechanical consequence of quenching is that the half-gap thickness at any of the elements is decreased in proportion to the number of gapwise nodes that have become frozen since the previous time step. At the end of each time step, all gapwise nodes are examined. If any are below the no-flow temperature, they are noted. Gap-wise nodes are then relocated such that they lie in the still-molten portion of the flow. Temperatures for the newly relocated nodes are determined by interpolation between the old gapwise node temperatures. In this interpolation, the wall temperature is increased as well. The gapwise node locations for two successive time steps is shown in Fig. 3.3

Such a reduction is gap thickness has two effects that can profoundly change the character of the filling. First, flow conductivity is strongly dependent upon half-gap thickness. For a Newtonian fluid, it varies as [b.sup.3] and for a power law fluid, as [b.sup.(2n + 1)/n]. With the use of freezing, the gap thickness decreases rapidly when thermal conduction terms are large compared with convection terms. Such would be the case as input volume flux to the die is decreased.

The second consequence is that freezing often alters the convergent/divergent character of the filling geometry. For example, filling between parallel plates becomes a combination of shear and extensional flow fields when freezing is occurring. The frozen layer just behind the flow front is comparatively thin because of the short time the melt has resided near the wall, while further upstream, the frozen layer is thicker. Still further upstream in the vicinity of the gate, the wall layer will be thin because of the hot melt's being convected in at the gate. If the material has the dual sensitivity character described previously, pressure distribution can be altered by inclusion of freezing.

The shedding of a layer of melt at the wall and effectively removing it from the flow is more realistic than declaring viscosity to become extremely large below a certain temperature for at least four reasons. First, the very large viscosity approach does not correspondingly alter the geometry of the flow passage. Elements having a uniform gap thickness cannot narrow or become convergent or divergent. Second, if the flow conductivity is evaluated at the bulk mean temperature for the element, the velocity and velocity gradient profiles inferred are unrealistic in the vicinity of the wall. Streamwise convection of energy in the vicinity of the wall becomes incorrect. Third, if the velocity profile is computed based upon a viscosity that varies from being very large at the wall to comparatively small in the warmer core flow, it is possible for the resulting profile to have an inflection point. Inflected velocity profiles can become unstable further downstream when momentum terms are included in the governing equation. When a strictly laminar analysis (as are most filling codes) generates an inflected velocity profile, concern over the appropriateness of the omission of momentum terms should arise. When it is realistic to declare the molding material to be solid rather than extremely viscous below some temperature, inflected profiles and their attendant concerns are far less common. Fourth, when viscous dissipation is included in the energy equation, there is an unrealistic heat source in the "frozen" region. In addition, the inclusion of very large numbers in a numerical analysis can lead to data overflows and underflows during computation.


Verification of the use of the dual viscosity model was done by applying the filling simulation to a set of material, geometry, and process conditions for which results already exist in the literature, Hieber et al. (36) have presented predictions and experimental results for the filling of a slightly tapered plaque cavity by polypropylene (Hercules Pro-fax 6523). The finite element mesh used here is the same as that used by Hieber et al. (36), except that elements representing a plenum that existed in the actual die cavity are also included. The mesh is shown in Fig. 4.

The junction between the plenum and plaque portions of the cavity presents the flow with a sudden contraction. The modeling of this flow field was handled as follows. Owing to its much smaller resistance to volume flux than the plaque, the plenum was assumed to fill entirely prior to any filling of the plaque. Once the plenum had filled, the flow from the plenum into the plaque was assumed to be predominantly two-dimensional. A means for estimating the angle between the converging and recirculating flows in a sudden contraction is provided by Cogswell (37). It was here estimated to be 23 [degrees]. The flow passage height for the elements in the plenum was specified based upon this angle rather than upon the full height of the plenum.

For these simulations, the function m(T) had the form

m(T) = [m.sub.o][e.sup.([T.sub.a]/T)] (42)

rather than as shown in Eq 25. A power law fit was applied to the shear viscosity data. The material constants were (6)

n = 0.323

[T.sub.a] = 1470 [degrees] K

[m.sub.o] = 540 kg/ms [(1/s).sup.1 - n] (43)

The process conditions were

[T.sub.w] = 300 [degrees] K

[T.sub.o] = 470 [degrees] K

Q = 7.1 [cm.sup.3]/s (44)

In terms of material constants to complete the dual viscosity model, pressure as a function of volume flow rate in capillaries having 10 [degrees] and 45 [degrees] half angles have been published by Kwon, et al. (23), for a similar material at similar temperatures. The data from the 10 [degrees] capillary were here reduced in the same manner as described previously, with the result that q = 1.407 sand c = 1/6.

The inclusion of the dual viscosity model improved the agreement between experimental results and predictions. Shown in Fig. 5 are the measured wall pressure histories and predictions of Hieber et al., along with predictions made here with both simple power law and dual viscosity models. The simple power law model agreed closely with the predictions of Hieber et al. Transducer 1 was located on the wall of the plenum. Its history was consistently underpredicted by the simple power law simulation, apparently because of the model's inability to include the material's extensional character in the flow transition from plenum into the plaque. The dual viscosity model was better able to account for this behavior.

The plaque portion of the cavity is tapered in the direction parallel to the line where the plenum and plaque are joined. The thickness gradient and velocity field, then, are nearly perpendicular during much of the filling, and extensional effects are comparatively small. Nonetheless, agreement was improved at Transducer 2 as well with the dual viscosity model.


The filling model was also applied to a die geometry and molding material properties typical of the powder injection molding process. The results from these filling cases have been assembled to more clearly illustrate the effects of dual viscosity and wall freezing and the process regimes in which they become most noticeable.

Die Geometry

The die geometry selected for testing the code was a sleeve-shaped part having a runner of square cross section, a fan gate, and a convergent transition section between the runner and gate. The geometry and its finite element mesh are shown in Fig. 6. This type of entry geometry is of particular interest because the convergent transition section is not uncommon in die design practice and it lends itself to the illustration of the effect of the implementation of the viscosity model in the finite element/finite difference formulation. The die total volume was 1.835 [cm.sup.3].

Test Matrix

Several runs were made to model the filling of the ring-shaped die geometry. The molding machine process parameters - die wall temperature, input volume flow rate, and melt temperature - were selected to match settings used in actual molding practice. They are listed in Table 1. The runs were selected such that the effects of wall freezing and dual viscosity [TABULAR DATA FOR TABLE 1 OMITTED] could be examined for their effects upon pressure and temperature when acting both alone and in consort.

Specific heat and thermal conductivity were taken from Fox (38). A plot of feedstock specific heat as a function of temperature is shown in Fig. 7. When wall freezing was not enabled, specific heat and the temperature-dependency of viscosity were treated differently than for the freezing-enabled cases. Specific heat was specified as a constant value having a magnitude equal to the specific heat above the melting peak. This has the consequence of making the amount of heat that must be removed from the melt to bring it to the no-flow temperature greater in the freezing-enabled case than in the freezing disabled case. Below the no-flow temperature in the freezing-disabled case, the value of m was buffered to be no larger than [10.sup.4] x A/[] + B.

All cases were run on a SPARCstation (Sun Microsystems, Mountain View, California) IPC running at 17.4 MIPS and having 8 megabytes of RAM. Each case ran to completion in batch mode in approximately three hours.


Dual Viscosity

The effect of the inclusion of a sensitivity to extension rate shows in the gate pressure histories and in pressure distribution contours. As soon as the flow front enters the convergent section of the gate, entrance pressure rises rapidly. When the dual viscosity model was enabled, the magnitude of the rise was larger. For this die geometry, this occurs at [similar to] 20% of the total time to fill. Once the convergent section of the gate has filled, the pressure traces for dual and power law models rise at approximately the same rate as each other as the remaining elements (all of which have parallel walls) are filled. The effect of the inclusion of extensional sensitivity becomes more pronounced with increasing fill rate. The magnitude of the difference in gate pressure histories between the power law and dual viscosity cases increases with increasing filling volume flux. Gate pressure histories for three filling rates and both viscosity cases are shown in Fig. 8.

Pressure across the die does not rise by the same factor or amount when the dual viscosity was enabled. Inclusion of a dual viscosity has the effect of increasing hydrostatic pressure gradients in the convergent section. Pressure contours on half-views of the plane of gate nodes are shown in Fig. 9. Shown are isobars just prior to the end of filling. At the high filling rate, pressure along the convergent section of the gate falls by 0.3 MPa in the power law case, whereas it declines by 0.9 MPa in the dual viscosity case.

Wall Quenching (Freezing)

The effect of the inclusion of quenching of melt at the die walls is more complicated than the effect of the inclusion of the dual viscosity. When compared to the case in which wall freezing is not enabled, it is possible for the resulting gate pressure histories to be lower. Such results are shown in Fig. 10.

The decrease is due in part to the way in which the temperature-dependent portion of the viscosity is evaluated in each case. In the freezing-disabled element, temperature is always evaluated from the mid-plane of the passage to the original die wall. In the freezing-enabled element, it is evaluated only from the plane of the element to the edge of the liquid zone. Since frozen regions at the wall have been shed from the flow path, the temperature distribution across the gap is more uniform and higher than for the freezing-disabled case. Bulk mean temperature contours for various freezing-disabled and freezing-enabled cases are shown in Fig. 11.

The second reason for the lower resultant gate pressure histories for the freezing-enabled cases is that, by the inclusion of an enthalpy of freezing, the melt carries more heat into the die than in the freezing-disabled cases. This has the effect of keeping temperatures higher further into the filling. The third reason is that the wall temperature becomes the temperature of recently frozen material and not the original wall temperature, and as a consequence, less heat is conducted into the wall.

The comparatively higher temperatures and narrower flow passages in freezing-enabled cases have opposing effects on flow conductivity. For the material thermal and constitutive properties used, the elevated temperatures were more than enough to offset the narrowing of flow passages. When freezing and the simple power law were enabled, the code would compute lower flow conductivities due only to the reduction in half-gap height (and not due to passage convergence) as freezing occurred. Contours of gap height ratio across the gate portion of the mesh at the end of filling are shown in Fig. 12. Gap height ratio is defined as mean elemental gap height at the end of filling over initial mean gap height.

The results here show the effect of wall quenching only; a simple power law viscosity model was used. The effects of dual viscosity and wall quenching occurred in different process regimes for the thermal and mechanical properties selected. The effect of the inclusion of wall freezing upon gate pressure was greatest at the low fill rate. As fill rate decreases, the characteristic time for freezing becomes smaller with respect to the time to fill, and more narrowing of the passage occurs. The effect of dual viscosity, however, was far less noticeable than the effect of wall freezing at the low fill rate because of the correspondingly low extension rates.


The numerical analysis described here is a more advanced predictive and diagnostic tool for molding die and process design. It serves to predict the pressure requirement for the filling of a part that is an assemblage of quasi two-dimensional flows. It includes a dual viscosity, freezing at the die walls, and a variable specific heat. The resulting pressure, temperature, and filling behavior are distinctly different than if these effects had not been included. While the viscosity model given in Eq 12 was incorporated into the numerical analysis, any model that is a function of material constants and invariants of the rate-of-deformation tensor could have been encoded.

The inclusion of the dual viscosity resulted in increased pressure gradients. In molding practice, large pressure gradients are known to have the potential to lead to molding difficulties and reduced-quality parts. In powder injection molding, power/binder separation can occur in such regions. In terms of molding, this creates cycle-to-cycle variation in injection pressure due to variation in solid volume fraction within the flow passage. In terms of part quality, a part having nonuniform density will typically not retain the desired finished dimensions after subsequent de-binding and sintering. In polymer molding, locations of severe pressure gradients are points where molecular uncoiling, orientation, and scission are likely to Occur.

The inclusion of wall freezing changes the deformation fields within the die, and the resulting effect on pressure and temperature is ultimately related to the various characteristic times relevant to the filling process. There is the time to fill the die, characteristic times associated with shear and extension rates, and the time for the melt to freeze in the presence of a characteristic thermal gradient. The reduction of these to dimensionless groups and their application to a rational approach to die and process design will be the topic of future work.


1. A viscosity model has been incorporated into a die filling model having differing sensitivities to shear and extensional deformation rates. The results are distinctly different than those derived from strictly shear-rate-dependent models.

2. When the candidate molding material has a dual-viscosity behavior, the results from the model described here may be more accurate predictors of difficult molding and poor quality molded parts than would the results from a model based upon a strictly shear-rate-dependent viscosity.

3. For the geometry and material properties selected, the effects of die wall quenching and dual viscosity became most pronounced at opposite extremes of the filling rate process envelope. Significant pressure increases due to an interaction between the two effects were not noted. However, this is not anticipated to hold true for all combinations of material and process characteristic time scales.

4. The approach to die wall quenching presented here results in more realistic temperature and enthalpy convection profiles than an approach involving the use of very large viscosities.

5. The conditions predicted in the die at the end of filling can be used for subsequent cooling and residual stress analysis. The inclusion of a dual viscosity and wall freezing alters the specification of these conditions.

6. The model generates results quickly in a modest computing environment and increases the number of die cavities whose filling can be realistically simulated.


The results reported in this paper are based upon work supported by the General Motors Corporation.


AR = Aspect ratio.

b = Element half-gap thickness.

C = Concentration.

[C.sub.p] = Specific heat.

D = Rate of deformation tensor = 1/2 [[nabla]V + [([nabla]V.sup.T]].

II = Second tensor invariant = 1/2[((tr D).sup.2] - tr([D.sup.2])).

h = Enthalpy.

k = Thermal conductivity.

I = Characteristic dimension of two-dimensional element.

L = Capillary length.

m = Material constant.

n = Material shear-thinning exponent or step count index.

P = Pressure.

q = Flux vector.

q = Material constant.

Q = Volume flow rate.

R = Radius.

S = Flow conductivity.

T = Temperature.

u = Velocity in x direction.

v = Velocity in y direction.

w = Velocity in z direction.

V = Velocity vector (u, v, w).

x = Cartesian coordinate.

y = Cartesian coordinate.

z = Cartesian coordinate.

[Alpha] = Half-gap height gradient in x direction.

[Beta] = Half-gap height gradient in y direction.

[Mathematical Expression Omitted] = Shear rate.

[Mathematical Expression Omitted] = Extension rate.

[Lambda] = Tensor eigenvalue.

[Mu] = Shear viscosity.

[Pi] = Total stress tensor.

[Chi] = Dual viscosity.

[Phi] = Dissipation function. ([[Tau].sub.ij]([Delta] [u.sub.i]/[Delta][x.sub.i])).

[Rho] = Density.

[Tau] = Deviatoric stress tensor.

tr = Trace.

D/Dt = Substantial derivative.

[nabla] = Gradient.

[absolute value of] = Absolute value.

[[ ]] = Magnitude.

[Mathematical Expression Omitted] = Time derivative.

[multiplied by] = Dot product.

[Mathematical Expression Omitted] = Average.


a = First vertex node or activation.

b = Second vertex node.

c = Third vertex node.

app = Apparent.

bs = Biaxial stretching.

ext = Extensional.

L = Capillary exit.

nf = No flow.

o = Interstitial fluid or at die inlet.

O = Capillary entrance.

s = Solid.

st = Straight.

tp = Tapered.

z = z direction.

1 = In direction of flow.

2 = Orthogonal to direction of flow.


1. G. Engelstein, Plast. Eng., April 1993, p. 41.

2. L. A. Najmi and D. Lee, Polym. Eng. Sci., 31, 1137 (1991).

3. J. R. Gaspervich, Int. J. Powder Metall., 27, 133 (1991).

4. C. M. Wang, R. L. Leonard, R. A. Posteraro, and T. J. McCabe, in Powder Injection Molding Symposium - 1992, P. H. Booker, J. Gaspervich, and R. M. German, eds., Metal Powder Industries Federation, Princeton, N.J. (1992).

5. B. Friedrichs, S. I. Guceri, S. Subbiah, and M. C. Altan, in Processing of Polymers & Polymeric Composites, Vol. 19, A. A. Tseng and S. K. Soh, eds., ASME, New York (1990).

6. C. A. Hieber and S. F. Shen, J. Non-Newtonian Fluid Mech., 7, 1 (1980).

7. F. N. Cogswell, J. Non-Newtonian Fluid Mech., 4, 23 (1978).

8. J. M. Dealy, J. Non-Newtonian Fluid Mech., 4, 9 (1978).

9. C. W. Macosko, in Fluid Mechanics Measurements, R. J. Goldstein, ed., Hemisphere, Washington, D.C. (1983).

10. M. J. Rosner, X. Zheng, M. Kojima, R. A. Posterato, and J. T. Lindt, in Powder Injection Molding Symposium - 1992, P. H. Booker, J. Gaspervich, and R. M. German, eds., Metal Powder Industries Federation, Princeton, N.J. (1992).

11. G. K. Batchelor, J. Fluid Mech., 46, 813 (1971).

12. J. Mewis and A. B. Metzner, J. Fluid Mech., 62, 593 (1974).

13. T. E. Kizior and F. A. Seyer, Trans. Soc. Rheol., 18, 271 (1974).

14. C. B. Weinberger and J. D. Goddard, Int. J. Multiphase Flow, 1, 465 (1974).

15. K. J. Mikkelson, C. W. Macosko, and G. G. Fuller, in Proc. Xth Intern. Cong. Rheol., Sydney, Australia (1988).

16. H. Munstedt, Rheol. Acta, 14, 1077 (1975).

17. K. Bald and A. B. Metzner, Trans. Soc. Rheol., 21, 237 (1977).

18. G. G. Fuller, C. A. Cathey, B. Hubbard, and B. E. Zebrowski, J. Rheol., 31, 235 (1987).

19. H. M. Laun and H. Munstedt, Rheol. Acta 18, 492 (1979).

20. P. R. Schunk and L. E. Scriven, J. Rheol., 34, 1085 (1990).

21. E. J. Hinch, Phys. Fluids, 20, (10), S22 (1977).

22. I. T. Barrie, SPE J., 27, 64 (1971).

23. T. H. Kwon, S. F. Shen, and K. K. Wang, Polym. Eng. Sci., 26, 214 (1986).

24. M. W. Evans and D. M. Heyes, Comput. Phys. Commun., 62, 249 (1991).

25. S. M. Dinh and R. C. Armstrong, J. Rheol., 28, 207 (1984).

26. R. B. Bird, R. C. Armstrong, O. Hassager, and C. F. Curtiss, Dynamics of Polymeric Liquids, Vol. 2, Wiley, New York (1977).

27. S. G. Advani and C. L. Tucker III, J. Rheol., 31, 751 (1987).

28. W. Ostwald, Kolloid-Z., 36, 99 (1925).

29. A. de Waele, Oil Color Chem. Assoc. J., 6, 33 (1923).

30. R. B. Bird, R. C. Armstrong, and O. Hassager, Dynamics of Polymeric Liquids, Vol. 1, Wiley, New York (1987).

31. A. E. Everage, Jr. and R. L. Ballman, Nature, 273, 213 (18 May 1978).

32. H. H. Winter, C. W. Macosko, and K. E. Bennett, Rheol. Acta, 18, 323 (1979).

33. B. Rabinowitsch, Z. phisik Chem., Abt. A, 145, 1 (1929).

34. J. C. Moller and D. Lee, Intl. J. Powder Metall., 30, 103 (1994).

35. J. C. Moller, PhD thesis, Rensselaer Polytechnic Institute, Troy, N.Y. (1994).

36. C. A. Hieber, L. S. Socha, S. F. Shen, K. K. Wang, and A. I. Isayev, Polym. Eng. Sci., 23, 20 (1983).

37. F. N. Cogswell, Polym. Eng. Sci., 12, 64 (1972).

38. R. T. Fox, PhD thesis, Rensselaer Polytechnic Institute, Troy, N.Y. (1992)
COPYRIGHT 1995 Society of Plastics Engineers, Inc.
No portion of this article can be reproduced without the express written permission from the copyright holder.
Copyright 1995 Gale, Cengage Learning. All rights reserved.

Article Details
Printer friendly Cite/link Email Feedback
Author:Moller, James C.; Lee, Daeyong; Kibbel, Bradley W.; Mangapora, Leslie
Publication:Polymer Engineering and Science
Date:Sep 1, 1995
Previous Article:Fracture toughness of acrylonitrile-butadiene-styrene by J-integral methods.
Next Article:Effect of mean strain on the cyclic deformation and stress relaxation in polypropylene.

Related Articles
Shear and elongational behavior of linear low-density and low-density polyethylene blends from capillary rheometry.
Elongational, hysteresis, and oscillatory flows of complex fluids.
Finite element analysis modeling of powder injection molding filling process including yield stress and slip phenomena.
Computer simulation of weld lines in injection molded poly(methyl methacrylate).
Control of polymer properties by melt vibration technology: a review.
Transient shear and extensional properties of biodegradable polycaprolactone.
Extensional flow of polymeric dispersions.
Verification of extensional viscosity effects in injection mold filling simulation.
Rheological characterization and process analysis.

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