A novel hybrid numerical technique to model 3-D fountain flow in injection molding processes

Journal of Non-Newtonian

Fluid Mechanics,

Elsevier

B.V., Amsterdam

Science Publishers

49 (1993) 141-113

141

A novel hybrid numerical technique to model 3-D fountain flow in injection molding processes ’ B. Friedrichs Department (Received

and S.I. Giigeri

of Mechanical October

*

Engineering,

University of Delaware,

8. 1992; in revised form February

Newark,

DE 19716 (USA)

15, 1993)

Abstract A novel hybrid two-dimensional (2-D)/three-dimensional (3-D) numerical technique is presented to model 3-D fountain flows in thin cavities as encountered in injection molding processes. At the fountain flow region, where all three velocity components are significant, the governing 3-D fluid flow equations are solved by using a pressure Poisson formulation. Behind the flow front, where out of plane flows are negligible, the 2-D Hele-Shaw formulation is employed, largely reducing the number of unknowns in comparison to a fully-three-dimensional formulation. Boundary fitted coordinate systems (BFCS) together with the finite difference method (FDM) are used to solve the governing equations on a non-staggered grid. The formulation is capable of handling non-linearities in the material behavior due to the shear-thinning characteristics of typical resin systems. Results are presented for isothermal flow of Newtonian and shear-thinning fluids through diverging and converging flow sections. Keywords: fountain flow; grid generation; molding; three-dimensional analysis

Hele-Shaw

flow; hybrid 2-D/3-D

method;

injection

1. Introduction Modeling of injection molding has been an area of considerable research over the last two decades. These studies resulted in several computational algorithms that assist in the design and fabrication of plastic parts using injection molding processes [ 11. Modeling of the flow in these processes represents several major challenges since the flow is inherently transient, * Corresponding author. Present address: Mechanical Engineering, Chicago, Chicago, Illinois 60607, USA. ’ This paper is dedicated to Professor Hermann Janeschitz-Kriegl 0377-0257/93/$06.00

0

1993 - Elsevier

Science Publishers

University

of Illinois at

on his 70th birthday.

B.V. All rights reserved

141

B. Friedrichs and S.I. Giigeri 1 J. Non-Newtonian

Fluid Me&.

49 (1993) 141-l 73

non-isothermal, non-Newtonian, and includes a free surface moving through cavities of highly irregular geometries. Due to limited computational resources, the three-dimensional (3-D) flow problem has customarily been simplified to a two-dimensional (2-D) problem, based on the observations of Hele-Shaw [2]. In this approach the cavity is assumed to be thin, and out of plane flows are neglected. It was first proposed for flow of polymer melts by Richardson [3] who looked at isothermal, Newtonian flow in simple geometric shapes. Over the years the Hele-Shaw formulation has become the basis for many computational codes that either use a Eulerian mesh together with a control volume (CV) approach (also referred to as flow analysis network (FAN)) to track the motion of the flow front [4- lo] or remesh either locally at the flow front [ 1 I] or completely [ 12- 141 such that the mesh coincides at all times with the domain covered by the fluid. The Hele-Shaw formulation, however, cannot capture the details of the flow field whenever out of plane flows are present such as in the vicinity of the flow front (see Fig. 1). Here, fluid typically spills outward towards the walls due to the no-slip situation leading to the “fountain flow” phenomenon first called as such by Rose [ 151 and experimentally observed in injection molding as reported in Refs. 16-18. Whenever the transient location of fluid particles is of importance, as for temperature, degree of cure, or short fiber orientation predictions, the incorporation of the fountain flow phenomenon becomes necessary to achieve satisfactory accuracy in process modeling, see e.g. Ref. 19. Although the fountain flow can be partially accounted for in Hele-Shaw flows by

Hele-Shaw (2-D flow)

Fig. 1. Flow in a thin cavity.

fountain flow (3-D flow)

B. Friedrichs and S.I. GiiGeri 1 J. Non-Newtonian

Fluid Mech. 49 (1993) 141 -I 73

143

following heuristic approaches such as that of Lord and Williams [20], much attention has been directed towards a more accurate two-dimensional analysis of the flow in the gapwise direction only. Early contributors to the field were Bhattarcharji and Savic [21] who developed an analytical model for isothermal, Newtonian flow in a pipe that is dragged over a stationary piston. Further analytical work that relates to injection molding processes was done by Tadmor [22] for an isothermal power-law fluid and by Castro and Macosko [ 191 and Manas-Zloczower et al. [23] for a reactive fluid. It should be noted that all these studies assume a flat flow front. The shape of the flow front was allowed to develop freely in simulations where a tube or parallel plates are sliding over the fluid with a stationary flow front. While Givler et al. [24] and Coyle et al. [25] solved this problem for an isothermal, Newtonian fluid, Mavridis et al. [26] carried the modeling to isothermal power-law fluids. The same authors later extended the work to non-isothermal Leonov fluids [ 271. The advance of a fluid into a stationary cavity represents a more complex problem since the domain is continuously deforming and enlarging. A remeshing procedure was first employed by Behrens et al. [28] who solved for the case of an isothermal, Newtonian fluid. A similar technique was also used by Mavridis et al. [29] who solved for two Newtonian fluids forming a weld line and by Fauchon et al. [30] for a Newtonian fluid moving through a corner. More recently Hayes et al. [31] solved for the advancement of a reacting fluid into the two-dimensional cross-section of a T-joint. A stationary mesh in conjunction with the marker-and-cell (MAC) technique was used by Gogos et al. [ 181 to model the advancement of a power-law fluid. Kamal and co-workers [32, 331 utilized the same approach for non-isothermal flow of a White-Metzner fluid. While some of these studies treat the problem in detail, it should be noted that this kind of analysis does not give any information about the nature of the planar flow through the mold cavity. Therefore, the predictions can only give qualitative understanding in a true three-dimensional flow situation. It is the purpose of this paper to bridge the gap between existing in-plane and gapwise analyses by solving for the fully three-dimensional flow field. To achieve computational economy a novel hybrid numerical technique is proposed that only solves for the 3-D flow field wherever all three-velocity components are significant; elsewhere, the 2-D Hele-Shaw formulation is retained. 2. Governing equations 2.1. 3-D Jlow region Typical Reynolds numbers encountered in injection molding processes using polymers are much less than unity [ 341. Therefore, the inertia terms

144

B. Friedrichs and S.I. Giigeri 1 J. Non-Newtonian

Fluid Mech. 49 (1993) 141-I 73

can be neglected, leading to a quasi-steady-state creeping flow type formulation. In a three-dimensional flow configuration the governing equations become continuity:

V *24= 0,

(1)

momentum:

V*o=O;a=z-pl

(2)

r = r(j)?.

(3)

rheological

behavior:

The fluid behavior is modeled as generalized Newtonian, i.e. z,, = 2q(j)&/ ax, z, = y(lj)(au/ay + au/ax), . . . . Amongst various possibilities, the highly shear-dependent viscosity can be represented by the Carreau model [35] yI =

qcc,+ (qo- &)[ 1 + (a#]““-

‘)‘2),

where y1 is the dimensionless power-law index which is less than unity for shear-thinning fluids. It ranges from 0.2 to 0.8 for most polymer melts. Furthermore, i is the square root of i times the second invariant of the rate of strain tensor, I a time constant and q. and q= are the Newtonian viscosities at zero and very high shear rates, respectively. For the Newtonian case y1 is equal to unity. It should be noted that the governing equations include no explicit expression for the pressure which acts as a source term in the momentum equation and adjusts itself in a way to satisfy continuity everywhere in the domain [36]. In order to overcome the difficulties associated with the missing explicit equation for pressure, the pressure Poisson equation formulation of Harlow and Welch [37] has been adopted. The first step in this approach is to consider the transient momentum equation

au

pz+vp

Then,

v*

=v-z. taking

the divergence,

p;+vp =v.z,

changing the order of differentiation this equation becomes 0 f

and considering

a finite time increment,

A

v.~n+‘_v.Un

v*p =V*(V*z)

-p

At

.

(7)

To ensure continuity at the new time step 12+ 1, the first term in the numerator on the right-hand side of eqn. (7) is set to zero. Equation (7) now gives the expression needed for pressure.

B. Fr~~drichs and XI. G&eri 1 J. ~o~-~e~~onian

Fluid Mech. 49 (1993) 141 I 173

14.5

This method is suitable for transient [37] as well as for steady-state [38-401 problems. It should be noted that for (quasi) steady-state problems, the transient term in the momentum equation (5) is only reintroduced as part of the pressure Poisson fo~ulation and does not alter the (quasr) steady-state character of the original problem. Equation (5) is then marched in time while updating the pressure (eqn. (7)) and the viscosity (eqn. (4)) until (quasi) steady-state is reached. The transient of this marching procedure, however, have no physical meaning. Once (quasi) steady-state is reached, the flow front is advanced and eqns. (S), (7) and (4) are solved again. It should be noted that the pressure Poisson equation (7) is not obtained from physical principles but derived mathematically; therefore, the boundary conditions on the pressure need to be consistent with the original governing equations which only allow for the specification of either the velocity or the traction on the boundaries [41]. These boundary conditions can be obtained by taking the normal component of the momentum equation as Vj?*n=(V*r)*n,

(8)

leading to Neumann type conditions straints to the system [42].

which introduce

no additional

con-

region

2.2. 2-Djow

Whenever the gapwidth is small in comparison with the planar dimensions, and the gapwise velocity component is negligible, an order of magnitude analysis shows that the momentum equations (2) can be reduced to

-.&-+& [,,,,(,,,,,lco,

_ aP(x> Y)

x momentum:

(9) aP(x9 Y) -~+~[“ai;(x~~,-)]=O,

y momentum:

\’

leading to the Hele-Shaw formulation. Integrating eqn. (9) with respect to z twice and using the no-slip boundary condition at the top and bottom walls, i.e. u = v = 0 at z = &h/2, gives hi2

8P 4x,

Y,

z)

=

-

ax

aP D(X,p, 2) = - -

ay

rw*

i

d,-*

~

sr

3 Y

s m

-

7

*

&

--y’

*

B. Friedrichs and XI. GiiGeri 1 J. Non-Newtonian Fluid Mech. 49 (1993) 141-I 73

146

Defining

fluidity

as

h/2$ dZ -

S(x, J) = s0

YI

the mean velocities

(11) ’

can then be found

by performing

integration

2s ap U=-Tax7

(12)

2s ap “=-we

Defining

the stream

function,

which satisfies continuity,

w u=Y’

V=-&,

as

(13)

W

and combining stream function

it with eqn. (12), a unified is obtained as

flow equation

in terms of the

(14) Although this equation can also be written in terms of the pressure, the stream function formulation is preferred in this study, because a prescribed flow rate boundary condition is applied at the inlet. 3. The hybrid solution technique Experimental observations have shown that the fountain flow region extends only of the order of one gapwidth backward from the free surface [25]. This region is, therefore, small in comparison to the planar length scales of the flow field (see Fig. 1). Thus, it is advantageous to only solve for three-dimensional flow in the fountain flow region using a three-dimensional mesh which is trailed by a 2-D mesh to solve for the Hele-Shaw flow. Figure 2 depicts a typical hybrid mesh generated for a planar how. This “hammerhead” approach therefore greatly reduces the number of nodes in the system, thus rendering a 3-D fountain flow analysis suitable for workstation based computations. A critical decision in constructing the hammerhead mesh is the determination of the pseudo-boundaries separating the two- and three-dimensional zones. Based on the experimental observations and the numerous numerical studies on 2-D fountain flow mentioned before, a length of 1.5-3 times the

B. Friedrichs and S.I. Giiqeri 1 J. Non-Newtonran

Fluid Mech. 49 (1993) 141-l 73

147

Flow

-(Hele-Shaw

2-D Mesh approach)-, 3-D Mesh (primitive variable approach)

Fig. 2. Hammerhead situation.

mesh resulting

from the hybrid

solution

technique

in a fountain

gapwidth was chosen for the extent of the 3-D flow region backward the contact line. The results later confirm this choice.

flow

from

4. Boundary conditions Figure 3 depicts various boundary conditions that need to be specified to solve the governing eqns. (5) and (14) for the fountain flow problem. At the inlet, the volume flow rate is specified leading to a Dirichlet type boundary condition for the stream function $. The boundary conditions at the side walls in the Hele-Shaw domain are also of Dirichlet type leading to full slip. This is consistent with the simplifications introduced in the momentum equations (9) which imply that the flow is unbounded in its plane. At the

interface: u = u(x,y,z) v = v(x,y,z) W=Q I

interface: g=

Cl(Y)

I

Top and bottom wall: u=v=w=o

side walls (Hele-Shaw): w = const.

Fig. 3. Boundary

conditions.

side walls (3-D): v_- Q= 0 (g._n)*t,=O (a- _n).t,=o

B. Friedrich

148

and S.I. Giiqeri 1 J. Non-Newtonian Fluid Mech. 49 (1993) 141~ I73

top and bottom walls, the no-slip condition holds which was used for the integration of the momentum equations (9). This interface between the 2-D Hele-Shaw and the 3-D flow domain behind the flow front is formed by an overlap of one grid layer in the direction of the flow. At the downstream end of the overlap, the normal derivative of the stream function is computed from the average velocities in the 3-D flow region using eqn. (13). At the upstream end of the overlap the planar velocity components u and u are computed from the stream function distribution using eqn. (10) written in terms of the stream function instead of pressure. Naturally, the gapwise velocity component w is zero. These conditions form the upstream boundary conditions for the 3-D flow region. At the top and bottom wall of that region, the no-slip condition is maintained. To be consistent with the Hele-Shaw formulation, a slip boundary condition is applied at the side walls which provides two boundary conditions. The third condition comes from the impermeability of the cavity boundaries. Since the pressure drops out from these equations, no additional boundary condition is needed. At the free surface a force balance gives [29,43] (15) where a and b denote positions on either side of the free surface (see Fig. 3). Ca is the capillary number defined as the ratio of viscous to surface tension forces Ca

=gU,

(16)

For creeping flow of highly viscous fluids this number is large (values for the surface tension y of various polymer melts can be found in Ref. 44) and eqn. ( 15) can finally be written as o-n

=o,

which needs to be four equations for the pressure at the determined from previous discussion equation.

(17) solved together with the continuity equation to provide the four unknowns U, 21,w and p. Although this specifies free surface, it is worthwhile to note that the pressure is the continuity constraint which is consistent with the on the boundary conditions for the pressure Poisson

5. The moving contact line Because of the no-slip boundary condition specified at the top and bottom walls of the cavity, there exists an infinite force singularity at the

B. Friedrichs and S.I. Giigeri 1 J. Non-Newtonian

contact point at time 1,

Fluid Mech. 49 (1993) 141-l 73

149

contact point at time tn+,

flow

Fig. 4. Flow front

advancement.

moving contact line. For a review of this problem see Dussan [45]. According to the same author, the problem is not a kinematic but a dynamic incompatibility, among the modeling assumptions such as no-slip at the wall and no stresses at the free surface. To the best of our knowledge there is no model that draws on the physics and resolves the issue. One way to remedy this problem is to introduce an artificial slip coefficient at the vicinity of the contact line [27,33,46-481. In this work, however, the approach of Behrens et al. [ 281 has been adopted which avoids the computation of the velocities at the contact line. Instead, the contact point is held fixed while the flow front is advanced using the computed velocities at time t, (Fig. 4). Then, a curve is fitted to the new location of the nodes along the flow front. The intersection of the curve with the cavity wall forms the new contact point at time tn+,. Thus, every position of the cavity wall is covered by a different fluid particle which resembles a rolling motion of the fluid as suggested by Dussan and Davis [49]. Although strict convergence of the solution cannot be established under strong spatial refinement, second-order convergence with respect to the flow rate at the entrance and the free surface was found for the practicable range of lo-30 nodes through the gap [50]. This suggests that the solution is sufficiently accurate with respect to the velocity field in this range of nodes. It should also be noted that the amount of fluid which appears to be lost through the walls computationally, is small in comparison with the discretization error for the overall problem and can therefore be neglected [28,50]. 6. Numerical implementation 6.1 Grid generation Boundary fitted coordinate systems (BFCS) in conjunction with the finite difference method (FDM) were chosen for the solution of the governing equations in which an irregularly shaped domain in the physical x, y, z space is mapped into a regularly shaped domain in the computational 5, r,~,v

B. Friedrichs and XI. Giigeri / J. Non-Newtonian

150

Physical

Domain

Fluid Mech. 49 (1993) 141-173

Computational

Domain

X

Fig. 5. 3-D mapping of an irregularly computational domain.

shaped

physical

domain

into

a regularly

shaped

space (see Fig. 5). Numerical grid generation as proposed by Thompson et al. [51,52] is used to generate the grid in the physical domain. The grid generation equations for the 3-D flow domain can be written as v’< =g”P, V2q =g”Q,

(18)

V’,t = g=R, where P, Q and R are grid control functions that enable the control of grid line concentration. The quantities g”, gz2 and g3’ are the diagonal components of the contravariant metric tensor and defined as g” = v;l

. V(‘,

(19)

where 5’ = [, ye, v for i = 1, 2, 3, respectively. independent variables, eqn. ( 18) becomes g”(x;:

+ Px,) +g”(x,,,

dependent

+ Pyt) +gz2(yqq +

= 0,

Q.Y,J +g33(~,,,. + Ryv)

(20)

+ 2(giZyc:, +g’3_v,;,,+ gZ3y,,.) = 0, g”(z<; + P,;) +g=(+

and

+ Qx,) +g3’(x ,,,,+ Rx,.)

-t 2(g ‘%~o + g?Ycv + g%,,) g”(y,,

Interchanging

+ Qz,, +gj3(zv,, + Rz,)

+ 2(g ‘~~~,]+ g ‘3z;,,+ g~3z,zr)= 0, which can be solved using second-order

finite difference

expressions.

B. Friedrichs and XI. Giiferi 1 J. Non-Newtonian Fluid Mech. 49 (1993) 141-l 73

151

The boundary conditions in the present study are of Dirichlet type, such that grids need to be specified on the bounding surfaces of the domain. These can be generated on a v = const surface, for example, by using [52]:

(21) which are similar to eqn. ( 18) except that the Laplacian is replaced with the second-order Beltramian. Interchanging dependent and independent variables renders gl’(xC; + Px;) + g2’(x,, + Qx,) +g12xrq = n’;(k, + k,,), +

Qv,> +g’2ysq = &(k + k,h

g”h

+ Pyt) +g”b,,

g”(+

+ Pq) + g”(z,, + Qz,) + g”z;, = n;(k, + k,,),

(22)

which, again, can be solved using centered, second-order finite difference expressions. In the above expressions n’;, n; and n; are the components of the outward unit normal vector to a v = const. surface and k, + k,, is the sum of the principal curvatures to be evaluated from an equation of the form F(x, y, z) = J?. But on the free surface such a function is not known. Although it is, in principle, possible to determine it locally from the location of surrounding nodes, this is too inaccurate for the present purpose. Instead, the nodes on the free surface are merely redistributed along the two varying computational coordinates on that surface according to the specified distribution along the bounding edges. This avoids undesirable node concentration due to the motion of the flow front. A particularly useful characteristic of the above generation equation is that the grid control functions P, Q and R can be determined from the distribution on the edges of the domain and then interpolated into the interior using linear (2-D) and transfinite (3-D) interpolation. The control function P, for example, can be found in l-D, i.e. along a varying < edge, as [52] p=_*

(23) x; which follows directly from eqn. (20). Node concentration along a varying 5 edge towards the i = m node, for example, can be achieved by using an expression of the form [53]

s

$+3(q,

(24)

where s, is the total length, s the incremental length and B the concentration factor. In this study eqn. (24) has been used to concentrate nodes towards

152

B. Friedrichs and S.I. GiiGeri / J. Non-Newtonian Fluid Mech. 49 (1993) 141-I 73

the contact lines in the direction of the flow. Additional information on three-dimensional grid generation and mapping in materials processing can be found in Refs. 50 and 54. The generation of the grid in the 2-D Hele-Shaw domain follows along the same lines and is therefore not discussed here any further. 6.2 Transformation

and discretization

of the governing equations

In order to achieve maximum accuracy for a given number of nodes and possible high viscosity gradients, the governing equations (5) (7) and (4) should be solved in their conservative form [55]. In other words, the coefficients of the derivative terms, such as viscosity, must either be constant or, if variable, their derivatives may not appear anywhere in the partial differential equation (PDE) [56]. Also to retain conservation in the discretized version on an irregularly shaped domain, the derivatives of a quantity fin 3-D, with respect to the physical coordinates x, y and z have to be stated in their geometrically conservative form which can be written as [ 571

JJ;- = (fJL>t + (fJ~ls>~+ (fJv.A, Jfi = (fJi”,>t+ CLJry),+ (fJy&

(25)

JL = (fJSz>t+ (fJ~z>~+ (fJv,h, where J is the Jacobian of the transformation and c,, q,, . . . are the components of the contravariant base vectors which can be written in terms of the covariant base vectors as V+z,

x a,+)

J = a, - (a2 x a3).

(i = 1, 2, 3)

(i, j, k) cyclic,

(26)

Here, a = (x~,, yrl, z;,. These expressions permit the cancellation of the fluxes across common cell boundaries, thus ensuring overall conservation. For second-order derivatives eqns. (25) need to be applied in a consecutive manner. For example, ar, /ax, is discretized as a a.y LX

(27)

B. Friedrichs and XI. Giigeri / J. Non-Newtonian

Fluid Mech. 49 (1993) 141-l 73

153

It is important that the metrics at the half-integer points may not be averaged but computed. Otherwise, a uniform flow field is not preserved [57]. Therefore, the nodal coordinates of the half-integer points need to be found by, for example, averaging. All other variables, such as viscosity, can be computed by averaging among the values at the full integer nodes. The discretization of Neumann type boundary conditions follows along the same lines. A good review on the topic can be found in Ref. 57. 6.3 Solution

technique for algebraic equations

The momentum equations are solved using the explicit Euler scheme for the time derivative. The maximum time step is limited by [58,59]

At<

(

-+ r.J,k

>min

if 4,,,k < 0,

(28)

is the multiplication factor of the velocity component of interest where Al.J,k in the discretized version of the right-hand side of eqn. (5). Although this time step is small, it is similar to the process time and other techniques such as fully implicit or approximate factorization show no savings in CPU time [50]. The convergence criterion was chosen to be Un+‘+Un

I

un

I-.
The pressure Poisson equation and the grid generation solved using successive over-relaxation (SOR).

(29) equations

are

6.4 Overall solution procedure The flow simulation is initialized with a small prespecified mesh at the inlet and an initial guess for the flow variables. The momentum equation (5) is then marched in time while updating the pressure using the Poisson equation (7). During this procedure the viscosity and the stream function in the Hele-Shaw domain are updated periodically after a specified number of velocity updates. Once convergence is achieved the nodes on the free surface are advanced using the first-order explicit scheme AX = u At, Ay = v At,

(30)

AZ = w At, while holding the nodes along the contact line fixed, as mentioned before. It should be noted that the time step in eqn. (30) is user-specified and only used to advance the flow front. The momentum equation (5), instead, is

154

B. Friedrichs and S.I. Giiqeri / J. Non-Xewtoniart

Fluid Mech. 49 (1993) 141-l 73

marched using a different time step computed from eqn. (28). Every coordinate line on the free surface that may intersect the cavity’s top and bottom wall is fitted using Akima’s method 1601. This is a one-dimensional interpolation from a set of given data points using piecewise third-order polynomials which has proven to produce smooth results. The new contact points, i.e. the new contact line, is found by intersecting the fitted curves with the cavity walls. To avoid undesirable node concentration, the nodes along each computational coordinate line on the free surface are redistributed algebraically. Finally, the mesh is regenerated on the boundaries and in the interior using the new position of the flow front. Whereas a 3-D mesh is generated downstream of the interface of the 2-D and 3-D flow domain, a planar 2-D mesh is generated upstream in the Hele-Shaw domain. Using the field variables from the previous time step as the initial guess for the next time step, the procedure is repeated until the cavity is filled. Since the domain is continuously deforming and enlarging, lines of nodes are added periodically in the downchannel direction of the Hele-Shaw domain to avoid either too fine or too coarse meshes. To achieve this, nodes are distributed along the cavity side walls which not only specify the cavity shape but also trigger the addition of node lines when passed by the flow front. It should be noted that the number of grid nodes in the 3-D flow domain is held constant due to only small changes in its size in the downchannel direction. 7. Results Isothermal flow of a Newtonian and a shear-thinning fluid through diverging and converging cavity sections of constant gapwidth is considered to demonstrate the proposed technique. Results are presented for the shear-thinning fluid unless stated otherwise. The viscosity of the Newtonian fluid was chosen to be 100 Pas. Material data for the Carreau model were taken for a polystyrene (PS) melt at 18OC with q. = 14 800 Pas, qL = 0 Pa.s, 3, = 1.04 s, IZ= 0.4 [ 121. The volume flow rate at the inlet was chosen to be 200 cm3 s-’ for the diverging flow configuration in order to obtain wall shear rates of 0 ( IO” s-‘) found in typical injection molding processes. Figure 6(a) shows the top view of the cavity shape and the initial mesh at the inlet. Due to symmetry, only the left half is considered. Figure 6(b) shows an intermediate mesh while Fig. 6(c) shows the final mesh after 0.3 s of mold filling using a maximum time step of 0.001 s for the flow front advancement. An enlarged side view of the downchannel symmetry line of the final mesh is depicted in Fig. 6(d). Here, the interface between the 2-D and the 3-D flow region is clearly seen. Although only a planar 2-D mesh

B. Friedrichs and S.I. Gibgeri /J. Non-Newtonian Fluid Mech. 49 (1993) 141-173

1.55

is needed in the Hele-Shaw domain, a grid is generated algebraically in the gapwise direction. The nodes of this grid form the collocation points for the fluidity integral given by eqn. (11) and the velocity integral given by eqn. (10). This grid, however, is not needed for the stream function computation,

Fig. 6 (a).

Fig. 6 (b).

156

B. Friedrichs and S.I. Giiceri 1 J. Non-Newtonian Fluid Mech. 49 (1993) 141-I 73

interface 2-D

r-1

3-D

Fig. 6. (a) Initial mesh and boundary of cavity; (b) intermediate mesh; (c) final mesh; (d) side view of final mesh along the downchannel symmetry line at the flow front.

thus reducing the mesh in Fig. 6(d)to the hammerhead mesh in Fig. 2. It should be noted that grid lines are concentrated strongly toward the flow front using eqn (24). In order to generate a high quality grid in the semicircular part of the fountain flow region, the last 4 grid lines/layers in the downchannel direction are collapsed into the contact point/line (see Fig. 6(d). The initial mesh was chosen to have a semicircular flow front with 1.5 nodes in the gapwise direction. Figure 7 shows the dimensionless flow front elongation which is defined as Elongation

=

Downchannel

distance

between tip the flow front and contact half gapwidth

line (31)

B. Friedrichs and XI. Giiqeri 1 J. Non-Newtonian

00

;

I

0 00

0 05

Fluid Mech. 49 (1993) 141 -I 73

I

010

I

time

Fig. 7. Dimensionless

flow front elongation

I

0 15

020

157

1

025

0 30

[s]

for diverging

flow configuration.

For a flat flow front this ratio is zero while it is unity for a semicircular front. Values for the dimensionless flow front elongation were taken at the center of the flow front. It is seen that the front shape remains close to semicircular in accordance with Hoffman [61] and experimental observations for flow in a tube [28]. Figure 7 also shows that in the straight section of the channel a constant value of 1.Ol for the Newtonian and 1.15 for the shear-thinning fluid are held rather steadily after some initial oscillations. They only drop slightly after the diverging section has been reached. These values correspond well to a grid refinement study performed for 2-D parallel plate flow using the same formulation [50]. It is seen in Fig. 8 that for a Newtonian fluid using between 15 and 41 nodes in the gapwise direction, the values lie well within the reported range of 0.84- 1.04 in other numerical studies [25,26,28,29,30,62,63]. Up to now, no experimental data have been reported for the elongation in parallel plate flow. Behrens et al. [28], however, were able to simulate the elongation to within the error of their experiments for pipe flow. Using the same numerical technique for parallel plate flow they predicted a value of 0.94. The results of this work are seen to converge towards this value. In addition, no experimental data have been reported for shear-thinning fluids at realistic injection molding process conditions. But, the asymptotic value of 1.04 reached for the shear-thinning

B. Friedrichs and S.I. Giigeri / J. Non-Newtonian Fluid Mech. 49 (1993) 141-I 73

158 14

: z

oa-

3 0 c

06-

: al -E

04-

0 .-

E ;

.-0

02-

00

#

of

I 40

I 30

I 20

10

nodes

in

Fig. 8. Convergence of flow front elongation

the

gapwise

direction

for 2-D parallel plate flow.

corresponds well to the work of Mavridis et al. [26] who predicted a value of 1.O for a power-law fluid with an exponent of 0.5 under similar flow conditions. Figure 9(a) shows the velocity u in the .X direction at the end of the simulation for the shear-thinning fluid. The drawing is enlarged in the gapwise direction for clarity. Thus, the flow front is stretched into an almost straight line. Only the bottom half of the domain is displayed due to symmetry; therefore, the top surface in Fig. 9(a) corresponds to the midplane of the cavity. Acceleration of the fluid is noticeable at the bend forming the entrance to the diverging section. The parabolic velocity profile and the slowdown of the fluid in the fountain flow region can be observed at the side. Figure 9(b) shows this more clearly in a side view of the fountain flow region. Here, both the top and bottom half are displayed and the figure is enlarged in every direction so that the real flow front shape is depicted. Figure 9(c) depicts the same situation for the Newtonian fluid. It should be noticed that the upstream, parabolic velocity profile is more protruding. Moreover, the deceleration occurs closer to the flow front. These differences will have a large impact on fluid particle relocation as demonstrated later. polystyrene

B. Friedrichs and XI. G&eri 1 J. Non-Newtonian Fluid Mech. 49 (1993) 141-l 73

159

Fig. 9. (a) Velocity component u at the bottom half of the cavity for the shear-thinning fluid (min. O.Om s-I max, 1.614ms~‘); (b) side view of (a) at the flow front; (c) side view for Newtonian fluid (min, 0.0 m s-l, max, 1.423m S-‘).

Fig. 10. (a) Velocity component u at the bottom fluid (min, 0.0 m s-‘; max, 0.557 m SC’).

half of the cavity

for the shear-thinning

160

B. Friedrichs and S.I. GiiGeri 1 J. Non-Newtonian

Fluid Mech. 49 (1993) 141-l 73

The other planar velocity component u in the y direction is depicted in Fig. 10. Looking at the mid-plane, it can be observed, again, how the fluid is slowing down when approaching the flow front. The slight deflection in the velocity contours is due to the transition from the Hele-Shaw to the 3-D flow domain. While the velocities in the 3-D domain are solved using second-order finite difference expressions, they are solved by Romberg integration [64] of eqn. (10) in the Hele-Shaw domain. Grid refinement was shown to improve the smoothness [50]. Figure 11 shows the gapwise velocity component u’ in the z direction in full view. Figure 11(b) shows the side of Fig. 11(a) in the fountain flow region. Here, the spilling motion can be observed where fluid is flowing towards the cavity’s top and bottom walls. It is seen that the out of plane flows are restricted to the immediate vicinity of the flow front, thus justifying the previously mentioned minimum size of the 3-D flow region of 1.5 times the gapwidth. Figure 1 l(c) shows the gapwise velocity component

Fig. 11. (a) Velocity component w for the shear-thinning fluid (min, -0.054 m SK’, max. 0.054 m SC’); (b) side view of (a) at the flow front; (c) side view for Newtonian fluid (min, -0.086 m SK’; max, 0.086 m SK’).

B. Friedrichs and XI. Giigeri 1 J. Non-Newtonian

Fluid Mech. 49 (1993) 141-f 73

161

0.03000 0.02500 0.02000 0.01500 0.01000 0.00500 0.00000 Fig. 12. Stream function distribution (kin, 0.0 s-‘; max, 0.040 s-l).

in the Hele -Shaw domain

Fig. 13. Viscosity distribution at the bottom (min. 144.2 Pas; max, 14 796.9 Pas).

for the shear-thinning

fluid

half of the cavity for the shear-thinning

fluid

for the Newtonian fluid. Due to the more protruding shape of the upstream velocity profile, the out of plane flows are more pronounced. The maximum value of the gapwise velocity is also seen to be located closer to the flow front. For completeness, the stream function distribution in the Hele-Shaw domain is shown in Fig. 12.

162

B. Friedrichs and S.I. Giigeri 1 J. Non-Newtonian Fluid Mech. 49 (1993) 141-l 73

The viscosity distribution is depicted in Fig. 13 in the bottom half of the domain. Due to extensional flow in the diverging section, the viscosity drops along the mid-plane from the zero shear rate viscosity at the inlet. The deformation of fluid particles in the fountain flow region manifests itself in a further drop close to the flow front. It is worthwhile to note that the steepest viscosity gradients are found in the gapwise direction. Here, the viscosity drops by two orders of magnitude within a short distance from the mid-plane. This, almost shock-like behavior, necessitated the use of a conservation formulation. As a consequence, satisfactory results were achievable using only 15 nodes in the gapwise direction. The pressure distribution (Fig. 14) is seen to be smooth despite the fact that pressure and velocities were computed at the same nodal locations, i.e. on a non-staggered grid. Due to the continuity constraint, the pressure is observed not to be zero at the flow front, although the deviation is small. In order to visualize the impact of the three-dimensional fountain flow on the fluid particle movement, numerical tracers were introduced at a constant rate at the inlet. While in each batch 20 particles were introduced along the width of the inlet, 40 were introduced through the half-gapwidth, i.e. the top half of the domain. These particles were relocated using eqn. (30). Since the tracer particle locations do not need to coincide with the grid nodes, a second-order Taylor series expansion around the closest node was used to interpolate for the velocities. Figure 15(a) shows the position of some tracer particles for the shear-thinning fluid when the simulation was stopped.

2500000.

Fig. 14. Pressure distribution at the bottom (min, 134 903 Pa; max, 22 399 292 Pa).

half of the cavity for the shear-thinning

fluid

B. Fried6ehs and S.L Gii~'eri l J, Non-Newtonicm Fluid Mech. 49 (1993) 141-173

163

15 14

"K

13 12 11

[a)

~0

7

15

3

13

12 (b)

11 10

"~,

9

6 Fig. 15. (a) Tracer particle location at the top half of the cavity for thc shear-thinning fluid: tracer particles were introduced slightly above the mid-plane of the cavity; (b) tracer particle location for the Newtonian fluid.

164

B. Friedrichs and XI. Giiceri /J.

Non-Newtonian Fluid Mech. 49 (1993) 141-173

These particles were introduced at the inlet just above the mid-plane. The numbering of the batches is done in chronological order. Therefore, batch 1 was introduced first and batch 15 last. It is seen that batch 1 has already experienced the fountain flow while batches 2 and 3 are in the process of doing so. The situation is much different for the Newtonian fluid. Here, more batches have experienced the fountain flow due to the more protruding upstream velocity profile. As a consequence, batch 1 is found further upstream, because it has reached the fountain flow region earlier. The impact of the diverging flow is well observed with batch 4. While the part close to the downchannel symmetry line is only starting to feel the fountain flow, the part at the other side has already gone through it and has, therefore, reached the top wall. The slight waviness in batch 2 is a consequence of the velocity interpolation. Some particles that are supposed to end up right at the top wall do so slightly below, thus still moving with the fluid. The situation changes when particles are traced that were introduced halfway between the mid-plane and the top wall (Fig. 16). Here, the shear-thinning fluid is moving faster than the Newtonian. Thus, batches 1 and 2 have experienced more of the fountain flow in Fig. 16(a) than in Fig. 16(b). This can also be observed by taking a side view along the downchannel symmetry line of the cavity (Fig. 17). The plots are enlarged in the gapwise direction for clarity. The characteristic V shapes formed by the individual batches in such fountain flow situations [65] are seen to occupy different positions and to take on different shapes due to the different constitutive behavior of the fluids. While the V shapes formed by the Newtonian fluid are further upstream and greater in number (Fig. 17(b)), they are stretched out over a wider distance for the shear-thinning fluid (Fig. 17(a)). Figure 18 shows the streaklines of tracer particles introduced at 3 different positions along the mid-plane of the cavity. Here, the three-dimensionality is perhaps best demonstrated since particles of subsequent batches take on different paths through the cavity in the direction perpendicular to the flow. Again, the different fluid behaviors are seen to influence the particle relocation dramatically. To further demonstrate the capabilities of the solution technique the same cavity shape was chosen at the inlet but followed by a 30” converging section. The flow rate was lowered to 100 cm3 s-’ with all other input parameters remaining unchanged. Figure 19 shows the location of the tracer particles introduced at the mid-plane of the cavity. Again, more batches have experienced the fountain flow in the Newtonian case than in the shear-thinning case. The fluid, however, is deformed differently as compared to the diverging case (Fig. 15). The computations were carried out on an IBM RS6000~550. The CPU time was around 30 h for all the presented cases. It should be mentioned

Non-Newtonian

Fluid Mech. 49 (1993) 141 --I73

165

Fig. 16. (a) Tracer particle location at the top half of the cavity for the shear-thinning fluid; tracer particles were introduced half way between the mid-plane and the top wall of the cavity; (b) tracer particle location for the Newtonian fluid.

B. Friedrichs and S.I. GiiGeri1 J. Non-Newtonian Fluid Mech. 49 (1993) 141 -I 73

166

15

14

13

12

11

10

9

8

7

6

5

4

3

(a)

15

14

13

11

12

10

98765

(b)

Fig. 17. (a) Side view of tracer particle location along the downchannel symmetry the shear-thinning fluid; (b) side view of tracer particle location for the Newtonian

line for fluid.

that the Newtonian case can be speeded up, because the term V * (V * t) is zero and can therefore be dropped from the computation of the pressure Poisson equation (7). 8. Conclusions A novel hybrid 2-D/3-D numerical technique was presented to model 3-D fountain flow in thin cavities as encountered in injection molding processes. At the flow front, where all three velocity components are significant, the 3-D fluid flow equations for creeping flow of a generalized Newtonian fluid are solved using a pressure Poisson formulation. Trailing the region, a 2-D Hele-Shaw type formulation is used, leading to a “hammerhead” approach. Due to the small size of the 3-D flow region, significant reduction in nodes, i.e. unknowns, was achieved as compared to a fully three-dimensional simulation. The governing equations were solved on a non-staggered grid using boundary fitted coordinate systems in conjunction with the finite difference technique. Material non-linearities due to the shear-thinning behavior of typical thermoplastic resin systems are solved for by periodically updating the viscosity during the iteration process. Results were

B. Friedrichs and S.I. Giigeri 1 J. Non-Newtonian

Fluid Mech. 49 (1993) 141 -I 73

Fig. 18. (a) Streaklines of tracer particles at the top half of the cavity for the shear-thinning fluid; (b) streaklines of tracer particles for the Newtonian fluid.

167

168

B. Friedrichs and S.L Gi~geri ] J. Non-Newtonian Fluid Mech. 49 (1993) 141-173

Fig. 19. (a) Tracer particle location at the top half of the cavity for the shear-thinning fluid; tracer particles were introduced slightly above the mid-plane of the cavity: (b) tracer particle location for the Newtonian fluid.

B. Friedrichs and S.I. Giiceri 1 J. Non-Newtonian

Fluid Mech. 49 (1993) 141-l 73

169

presented for the isothermal filling of diverging and converging cavity sections using a shear-thinning and a Newtonian resin. The impact of the 3-D fountain flow on the flow field was studied by tracking numerical tracer particles. The observed flow patterns compared well to existing two-dimensional numerical and experimental studies. Acknowledgments We are grateful to Dr. S.D. Gilmore and Mr. C.A. Moore for their help with computer graphics. This work was partially funded by the Graduate Studies Office at the University of Delaware. Nomenclature

2

w.k

B Ca f

g” h

iA k J

k k,, 2

m n n

P

P, Q, R s St

s t

u, 0, w U

f.4

v

u x, Y, z

covariant base vectors multiplication factor concentration factor capillary number field variable diagonal components of the contravariant metric gapwidth indices Jacobian principal curvatures maximum value of i dimensionless power law index outward unit normal pressure grid control functions incremental length total length fluidity time velocity components in the x, y, and z directions velocity vector mean velocities in the x and y directions characteristic velocity coordinates in the physical domain

tensor

Greek letters Y

i

surface tension square root of l/2 times the second invariant tensor

of the rate of strain

B. Friedrichs and S.I. Giigeri 1 J. Non-Newtonian

rate of strain tensor viscosity viscosity at zero shear rate viscosity at infinite shear rate time constant coordinates in the computational density total stress tensor stress stream function

Fluid Mech. 49 (1993) 141- 173

domain

References 1 L.T. Manzione (Ed.), Application of Computer Aided Engineering in Injection Molding, Hanser Publishers, Munich, 1987. 2 H.S. Hele-Shaw, The flow of water, Nature (London). 58 (1898) 34-36. 3 S. Richardson, Hele-Shaw flows with a free boundary produced by the injection of fluid into a narrow channel, J. Fluid Mech., 56(4) (1972) 609-618. 4 E. Broyer, C. Gutfinger and Z. Tadmor, A theoretical model for the cavity filling process in injection molding, Trans. Sot. Rheology, 19(3) (1975) 423-444. 5 C.A. Hieber and S.F. Shen, A finite-element/finite-difference simulation of the injection molding filling process, J. Non-Newtonian Fluid Mech., 7 ( 1980) l-32. 6 H.H. Chiang, C.A. Hieber and K.K. Wang, A unified simulation of the filling and post-filling stages in injection molding. Part I: Formulation. Polym. Eng. Sci., 31(2) (1991) 116-124. 7 J.-L. Willien, J.-F. Agassant and M. Vincent, Numerical simulation of mold filling for thermoplastic injection, The Sixth Annual Meeting of the Polymer Processing Society, PPS, Brookfield, CT. 1990. 8 P. Filz, Simulation of the injection moulding process for crosslinking materials-Current situation and future possibilities, Kunststoffe, 79( 10) ( 1989) 1057- 1061. 9 T. Matsuoka, J.-I. Takabatake, A. Koiwai, Y. Inoue, S. Yamamoto and H. Takahashi, Integrated simulation to predict warpage of injection molded parts., in ANTEC ‘90, Society of Plastics Engmeers, SPE, Brookfield, CT, 1990, pp. 369-374. 10 P.A. Tanguy and R. Lacroix, A 3D mold filling study with significant heat effects, Int. Polym. Proc., 6( 1) (1991) 19-2.5. 11 A. Couinot, L. Dheur, 0. Hansen and F. Dupret. A finite element method for simulating injection molding of thermoplastics, m J.F. Thompson (Ed.), Numiform 89, Balkema, Rotterdam, 1989. 12 S. Subbiah, D.L. Trafford and S.1. GiiGeri. Non-isothermal flow of polymers into two-dimensional, thin cavity molds: A numerical grid generation approach, Int. J. Heat Mass Transfer. 32( 3) (1989) 415-434. 13 B. Friedrichs, S. i. GiiGeri, S. Subbiah and M.C. Altan, Simulation and analysis of mold filling processes with polymer-fiber composites; TGMOLD. in A.A. Tseng and S.K. Soh (Eds.), Processing of Polymers and Polymeric Composites, The American Society of Mechanical Engineers, ASME. New York, 1990, pp. 73391. 14 M.A. Garcia, C.W. Macosko. S. Subbiah and S.i. GiiGeri, Modelling of reactive filling in

B. Friedrichs and S.I. GiiGeri 1 J. Non-Newtonian

15 16 17 18 19 20 21

22 23 24 25 26 27 28 29 30 31 32

33 34 35 36 37

Fluid Mech. 49 (1993) 141-l 73

171

complex cavities: Comparison of fountain flow approximations. Int. Polym. Proc., 6( 1) (1991) 73-82. W. Rose, Fluid-fluid interfaces in steady motion, Nature, 191 (4785) (1961) 242-243. R.S. Spencer and G.D. Gilmore, Some flow phenomena in the injection molding of polystyrene, J. Colloid Sci., 6 (1951) 1188132. L.R. Schmidt, A special mold and tracer technique for studying shear and extensional flows in a mold cavity during injection molding, Polym. Eng. Sci., 14( 11) (1974) 7977800. C.G. Gogos, C. Huang and L.R. Schmidt, The process of cavity filling including the fountain flow in injection molding, Polym. Eng. Sci., 26(20) (1986) 1457-1466. J.M. Castro and C.W. Macosko, Studies of mold filling and curing in the reaction-injection molding process, AIChE J., 28( 2) (1987) 250-260. H.A. Lord and G. Williams, Mold filling studies for the injection mold filling of thermoplastic materials, Polym. Eng. Sci., 15(8) (1975) 5699582. S. Bhattacharji and P. Savic, Real and apparent non-Newtonian behavior in viscous pipe flow of suspensions driven by a fluid piston, in Proc. 1965 Heat Transfer and Fluid Mechanics Institute, 1965, pp. 248-262. Z. Tadmor, Molecular orientation in injection molding. J. Appl. Polym. Sci., 18 (1974) 1753-1772. Manas-Zloczower, J.W. Blake and C.W. Macosko, Space-time distribution in filling a mold, Polym. Eng. Sci., 27( 16) (1987) 1229-1235. R.C. Givler, M.J. Crochet and R.B. Pipes, Numerical prediction of fiber orientation in dilute suspensions, J. Composite Mater., 17 (1983) 330-343. D.J. Coyle, J.W. Blake and C.W. Macosko, The kinematics of fountain flow in mold filling, AIChE J. 33(7) (1987) 1168-1177. H. Mavridis, A.N. Hrymak and J. Vlachopoulos, Finite element solution of fountain flow in injection molding, Polym. Eng. Sci., 26(7) (1986) 449-454. H. Mavridis. A.N. Hrymak and J. Vlachopoulos, The effect of fountain flow on molecular orientation in injection molding, J. Rhoelogy, 32(6) (1988) 639-661. R.A. Behrens, M.J. Crochet, C.D. Denson and A.B. Metzner, Transient free-surface flows: motion of a fluid advancing in a tube, AIChE J., 33(7) (1987) 1178-l 186. H. Mavridis, A.N. Hrymak and J. Vlachopoulos, Transient free-surface flows in injection mold filling, AIChE J., 34(3) (1988) 403-410. D. Fauchon, H.H. Dannelongue and P.A. Tanguy, Numerical simulation of the advancing front in injection molding, Int. Polym. Proc., 6 (1991) 13-18. R.E. Hayes, H.H. Dannelongue and P.A. Tanguy, Numerical simulation of mold filling in reaction injection molding, Polym. Eng. Sci., 31( 11) (1991) 842-848. P.G. Lafleur and M.R. Kamal, A structure-oriented computer simulation of the injection molding of viscoelastic crystalline polymers. Part I: Model with fountain flow, packing, solidification, Polym. Eng. Sci., 26( 1) (1986) 92- 102. M.R. Kamal, SK. Goyal and E. Chu, Simulation of injection mold filling of viscoelastic polymer with fountain flow, AIChE J. 34( 1) (1988) 94-106. AI. Isayev (Ed.), Injection and Compression Molding Fundamentals, Marcel Dekker, New York, 1987. P.J. Carreau, Rheological equations from molecular network theories, Ph.D. Thesis, University of Wisconsin, 1986. S.V. Patankar, Numerical Heat Transfer and Fluid Flow, Hemisphere, Washington, D.C., 1980. F.H. Harlow and J.E. Welch, Numerical calculation of time-dependent viscous incompressible flow of fluid with free surface, Phys. Fluids, 8( 12) (1965) 2182-2189.

172

B. Friedrichs and S.I. Giigeri / J. Non-Newtonian

Fluid Mech. 49 (1993) 141-l 73

38 S. Abdallah, Numerical solutions for the incompressible Navier-Stokes equations in primitive variables using a non-staggered grid, II. J. Comput. Phys.. 70 (1987) 193202. 39 M.L. Mansour and A. Hamed, Implicit solution of the incompressible Navier-Stokes equations on a non-staggered grid, J. Comput. Phys., 86 (1990) 147- 167. 40 K. N. Ghia. W.L. Hankey and J.K. Hodge, Use of primitive variables in the solution of imcompressible Navier-Stokes equations, AIAA J., 17(3) (1979) 298-301. 41 P.M. Gresho and R.L. Sani, On pressure boundary conditions for the incompressible Navier-Stokes equations, Int. J. Numer. Methods Fluids, 7 (1987) 1111-1145. 42 P.J. Roache, Short communication: A comment on the paper “Finite difference methods by J.C. Strikwerda, Int. J. Numer. for the Stokes and Navier-Stokes equations” Methods Fluids, 8 (1988) 1459-1463. 43 R.L. Panton, Incompressible Flow, Wiley, New York, 1984. 44 B.B. Sauer and N.V. DiPaolo, Surface tension and dynamic wetting of polymers using the Wilhelmy method: applications to high molecular weights and elevated temperatures, J. Colloid Interface Sci., 144(2) (1991) 527-537. 45 E.B. Dussan, On the spreading of liquids on solid surfaces: Static and dynamic contact lines, Annu. Rev. Fluid Mech., 11 (1979) 371-400. 46 E.B. Dussan. The moving contact line: The slip boundary condition, J. Fluid Mech., 77(4) (1976) 665-684. 47 C. Huh and L.E. Striven, Hydrodynamic model of steady movement of a solid/liquid/ fluid contact line, J. Colloid Interface Sci., 35( 1) (1971) 855101. 48 H.P. Greenspan, On the motion of a small viscous droplet that wets a surface, J. Fluid Mech., 84( 1) (1978) 125-143. 49 E.B. Dussan and S.H. Davis, On the motion of a fluid-fluid interface along a solid surface, J. Fluid Mech., 65( 1) (1974) 71-95. 50 B. Friedrichs, Modeling of three-dimensional flow fields in injection and resin transfer molding processes, Ph.D. Thesis, Department of Mechanical Engineering, University of Delaware, 1993. 51 J.F. Thompson (Ed.), Numerical Grid Generation, Elsevier, New York, 1982. 52 J.F. Thompson, Z.U.A. Warsi and C.W. Mastin, Numerical Grid Generation-Foundations and Applications, Elsevier, New York, 1985. 53 J.P. Coulter. Resin impregnation during the manufacturing of composite materials., Ph.D. Thesis, Department of Mechanical Engineering, University of Delaware, 1988. 54 S.D. Gilmore and S.i. Gtiqeri, Solidification in anisotropic thermoplastic composites. Polym. Composites, 1 l( 6), December 1990. 55 B. Friedrichs and S.I. Giiceri. A numerical technique for the internal creeping flow of highly shear-thinning fluids in general curvilinear coordinates, Int. J. Comput. Fluid Dynamics, 1992, submitted for publication. 56 D.A. Anderson, J.C. Tannehill and R.H. Pletcher, Computational Fluid Mechanics and Heat Transfer, Series in Computational Methods in Mechanics and Thermal Sciences, Hemisphere. New York, 1984. 57 J.F. Thompson, Z.U.A. Warsi and C.W. Mastin, Boundary-fitted coordinate systems for numerical solution of partial differential equations-a review, J. Comput. Phys., 47 (1982) l-108. 58 S.D. Gihnore, Thermal and residual stress analysis in processing of thermoplastic composites, Ph.D. Thesis, Department of Mechanical Engineering, University of Delaware. 1991.

B. Friedrichs and S.I. Giigeri / J. Non-Newtonian

Fluid Mech. 49 (1993) 141-l 73

173

59 M.N. ijzigik, Basic Heat Transfer, McGraw-Hill, New York, 1977. 60 H. Akima, A new method of interpolation and smooth curve fitting based on local procedures, J. Assoc. Comput. Machinery, 17(4) (1970) 589-602. 61 R.L. Hoffman, A study of the advancing interface, J. Colloid Interface Sci., 50(2) (1975) 228-241. 62 K.K. Wang, S.F. Shen, C. Cohen, C.A. Hieber and S. Jahanmir, Computer-aided injection molding system, Progress Report 5, Cornell University, College of Engineering, Ithaca, NY, 1978. 63 K.K. Wang, S.F. Shen, C. Cohen and C.A. Hieber, Computer-aided injection molding system, Progress Report 6, Cornell University, College of Engineering, Ithaca, NY, 1979. 64 W.H. Press, B.P. Flannery, S.A. Teukolsky and W.T. Vetterling, Numerical Recipies, Cambridge University Press 1986. 65 A.N. Beris, Fluid elements deformation behind an advancing flow front, J. Rheol., 31(2) (1987) 121-124.