Simulation accuracy of a contraction flow of viscoelastic polymer liquid in liquid molding processes

Journal of Materials Processing Technology, 29 (1992) 351-371

351

Elsevier

Simulation accuracy of a contraction flow of viscoelastic polymer liquid in liquid molding processes Haiqing Gong School of Mechanical and Production Engineering, Nanyang Technological University, Nanyang Avenue 2263, Singapore

Selquk I. Gii~eri Thermal Engineering and Advanced Manufacturing Group, Department of Mechanical Engineering, University of Delaware, Newark, DE 19716, USA (Received May 15, 1991; accepted July 3, 1991 )

Industrial Summary To ensure the validity of the simulation results from computer software used in polymer processing, it is necessary to understand and be able to predict the simulation accuracy or suitability of constitutive models. This work is concerned with various issues relating to the accuracy of numerical simulations of viscoelastic polymer flows in liquid molding processes (RIM and RTM). A typical mold geometry, the abrupt contraction, is adopted for the investigation of the simulation accuracy of flows of Maxwell and 01droyd-B fluids through the contraction. The accuracy of the numerical solution is found to be low in the contraction region if upwind schemes are applied to stabilize the simulation. The loss of accuracy is due mainly to large gradients of the stress field in the contraction region where the artificial diffusivity is of the same order of magnitude as other terms in the constitutive model. The stress singularity at a sharp corner and a large stress gradient inherent in the Cartesian coordinate system are investigated analytically to determine the characteristics of the inaccuracy.

1. Introduction

Numerical simulation of viscoelastic fluid flows has become an increasingly important subject in resin injection molding (RIM) and resin transfer molding (RTM) processes in composite material manufacturing and other polymer processing technologies. Polymer liquid is composed of macro-molecular chains. During flow, the molecular chains experience various motions, such as entanglement, change of orientation, etc. Experiments show that the theological behavior of polymer liquid is quite different from that of the Newtonian fluid. One of the most prominent rheological characteristics of polymer fluid is its relaxation properties under stress or strain [ 1 ], whereas the Newtonian fluid 0924-0136/92/$05.00 © 1992 Elsevier Science Publishers B.V. All rights reserved.

352

and other viscometric non-Newtonian fluids do not exhibit any stress or strain relaxation. Another important characteristics of polymer liquid is its memory effect; i.e., the stress or strain in polymer liquids is a function of its history. Due to this property, flows of some polymer liquids can exhibit a recoil phenomenon, die swell, etc. [1]. Traditional efforts in the simulation technology of mold-filling processes are based on viscometric non-Newtonian constitutive models, such as the powerlaw fluids [ 2,3 ]. However, the viscoelasticity inherent in the polymer fluid flow was not incorporated into the simulations, and many important flow properties such as stress/strain relaxation, memory effect, large normal stress, etc. were not reflected in the numerical solutions. The importance of the viscoelasticity in simulations of the liquid molding process depends mainly on a dimensionless parameter, the Weissenberg number, defined by

we where 2 is the relaxation time of the polymer liquid, U is the characteristic velocity, and L is the characteristic length. Similar to the Weissenberg number, another commonly used dimensionless parameter, the Deborah number, is also given as

De=2/T

or

De.~.~,

where ~)is the characteristic shear rate, and T is the characteristic time of the flow process. Low values of De or We correspond to fluid-like behavior and high values to solid-like behavior. Note that even a liquid system with a small relaxation time, 2, can appear solid-like in a fast-deformation process, such as a polymer flow in complex geometries containing an abrupt contraction. Simulation of viscoelastic fluid flow is extremely difficult at large Deborah numbers because of high complexity and non-linearity of the governing equation, as well as stability problems of numerical schemes [4-8]. To overcome the numerical instability and eventually obtain numerical solutions, various forms of upwind schemes are used to stabilize the simulations. However, the damping of the numerical solution by the upwind schemes can cause severe accuracy problems, especially in the region where the mold boundary changes drastically. In such a region, the large velocity gradient induces a large stress gradient which contains a significant portion of solution components of large wave-numbers; therefore, the large stress gradient is subject to a severe damping by upwind schemes [9 ]. This work investigates various forms of large stress gradients in the contraction region and their influence on the accuracy of the numerical solutions. It shows that the artificial diffusivity term appearing in the first-order upwind scheme is of the same order of magnitude as other terms in the Maxwell constitutive model. Such an extra term can obviously create significant inaccuracy

353 in the contraction section, although the magnitude of the inaccuracy is unknown. Through a spectral accuracy analysis, the magnitude of the inaccuracy induced by the artificial diffusivity is characterized qualitatively at various stress gradient levels. The linear filtering of higher-order is found to be a more accurate stabilization method than the upwind scheme in the sense that it can preserve the solution components of lower wave-numbers. In addition, the stress singularity at a sharp corner and the large stress gradient inherent in the Cartesian coordinate system are investigated to exhibit the influence of various process parameters, such as We and contraction angle, on the magnitude of the stress gradient. 2. A simulation procedure

Mold flow of polymer liquid involves a change of fluid boundary, i.e., the fluid boundary expends as the flow proceeds. Such a moving-boundary flow problem is a transient process. So far, the numerical solution of a transient moving-boundary flow of viscoelastic liquid in two dimensions is impossible to obtain due to the limitation of existing computer capability. However, in the flow regions where the transients of flow properties become small, the governing equation can be treated as a steady-state one. Such regions include: (1) flow regions far away from the front of the expanding domain, where the relaxation of flow properties almost finishes; (2) regions where the mold boundary changes drastically, such as an abrupt contraction, so that the transient variations of flow properties is insignificant compared to the corresponding spatial variations. As far as the simulation accuracy is concerned, abruptly changed mold-boundary configurations, such as an abrupt contraction, are the most important ones to investigate, in the sense that large stress gradients in these regions cause severe discretization inaccuracy. Therefore, in this work, the investigation of simulation accuracy in the contraction region is conducted based on a steady-state flow process. To study the solution characteristics of viscoelastic flows in the region where the mold boundary changes drastically, a 4:1 contraction geometry is employed in this work. In the contraction region, the velocity gradient is large, which increases greatly the stiffness of the system of the governing equations and the corresponding numerical scheme. To resolve the large stiffness, iterative methods have to be strongly implicit. For this reason, an ADI consistent upwind scheme is proposed to simulate the flow of the Oldroyd-B fluid through a planar 4:1 contraction. In addition, a staggered mesh is used in the computation to avoid the two-grid wavy solution pattern [8] and achieve a higher computational efficiency. This simulation provides solution characteristics for further studies of the simulation accuracy in later sections.

354

2.1. Numerical formulation The current study is based on a primitive variable (u, v,p) formulation. This approach is advantageous computationally since it is equally applicable to twoand three-dimensional viscoelastic fluid flows. The governing equations of the steady planar flow of viscoelastic fluid with the Oldroyd-B model [ 1,4 ] can be written as

Ux-~-vy=O p( uux + vuy) +Px =ax + zy -~ ~ll (Uxx "~-Uyy) p(UVx + v v y ) +py =z~ + Ty + ~/l(Vx~ +vyy) ~ ( Uax + vay - 2au~ - 2 zuy ) + a = 2tl2 Ux

A ( u ~ + w y - a V x --Tuy) + z = t/2 (uy +v~) ,~(uT~ +VTy - 2 z v x - 27vy) + 7=2~72vy

where u and v are velocity components, p is pressure, a and 7 are extra normal stresses, z is extra shear stress, and ~ is a characteristic relaxation time for the fluid. The first equation is the continuity equation, the second and the third equations are the m o m e n t u m equations, and the last three equations represent the constitutive model. The above equation can be represented in a block form as

(I)

Aqx + Bqy + C q = G ( q x x +qyy )

where q = (p, u, v, a, p, ?)T

A=

B=

tO

1

0

0

1

pu

0

- 1

0 0 0 - 2 ()~a+~) o o o o [0 0 0 pv 1 0 0 - 2A~ 0 - (,~7+~) 0 0

0 0 -1 0 Au Au o 0 0 2u/

pu 0 - ( ~ a + tl) ( - 2~) 1 0 pv 0 0 -2(A7+~)

0 0 0 Av 0 0

0 0 - 1 0 0 -1 0 0 0 Av 2v 0

355

¢0 0 0

C=

1 1 1

fo }71 }71

G=

0 0 0

In t h i s work, }71/(}71 "Jff}72) = 0.111.

The artificial compressibility method, proposed by Chorin [10] for Newtonian flows and analysed by the authors [11] for viscoelastic flows, is used as a solution procedure for eqn. (1). Using this method the divergence-free constraint on the velocity field has been replaced by a time-evolution equation for the pressure, i.e.,

0p 8~-[+ux +vy=O Correspondingly, transient terms are added to eqn. (1), which forms a transient system, i.e.,

Dqt +Aqx +Bqy +Cq=G(qxx +q~ )

(2)

where the transient term Dqt is added to the left-hand side of eqn. (1), with

~u D--

~v ~

where ~p, ~,, ~v, ~o, ~, ~y are iteration parameters. Transient system (2) is used equivalently as an iterative procedure in solving system (1) [11]. The solution of system (1) is an asymptotic solution, when t-~ ~ , of transient system (2). The major advantage of using system (2) is that no additional iteration is needed to satisfy Vu = 0 at each time step. As

356

analysed by Gong and Giiqeri [11], system (2) is a well-posed mixed hyperbolic-parabolic system and it can be solved by existing numerical methods for compressible flow problems. To facilitate the simulations, a numerical grid generation technique is used. With this approach, a set of partial differential equations describes the mapping between the physical domain (x, y) and the computational domain (~, y). Generation of a solution using this technique begins by mapping the physical domain having complex geometry onto a computational domain having regular geometry, such as the rectangle used in the current study. This technique continues by posing an equivalent problem in the computational coordinate system. Then, the transformed governing equations are solved in the computational domain using traditional finite-difference methods. The mapping process generates a curvilinear coordinate system having coordinate lines, or coordinate surfaces, that coincide with the boundaries of the physical domain. This makes representing the domain shape and applying the boundary conditions easy. In general, numerical grid generation combines the geometric flexibility of finite-element methods with the computational simplicity of finite-difference methods. The following Poisson-type elliptic grid generation pair [12,13] is used in the simulation of eqn. (2). The domain transformation from (x, y) to (~, ~/) is given as,

,Txx

and the transformation from (~, ~/) to (x, y) is, fll X~ - 2,62 x~, + fla x,, 1+ j 2 ( Px~ + Qx.) = 0

(3)

flly~ - 2flays,7 + flay.. + j 2 ( py~ + Qy.) = 0

(4)

where the geometric coefficients are given as fll = x.2 + y .2, fla = xex. +y~y., fla = x ~ + y 3, and the Jacobian of the transformation is given as J = x e y , - x . y e . P and Q, the grid control functions used to attract mesh points or lines to areas where high resolution is desired [ 12 ], are given as, M

P(~, r/)= - ~ a ~ i ( ~ - ~ i ) e x p ( - Q i l ~ + c j 4 ) i=1 M

N

-- ~, ~., b¢,,~(~-~i)exp[-dc~,j[ ( ~ - ~ i ) 2+ (r/-rb) 211/21 i=1 j=l M

Q(~, r/) = - ~ a..(t/-t/i)exp(-c.Jrl+c,v/rlj ) i=1 M

- •

N

E b,7,.(~l-~b)exp[-d.,.J[(~-~i) 2+ (~/-~b) 21 1 / 2 ]

i=lj=l

357

The effect of the coefficient a¢i is to attract the constant-~lines to a specified ~i line. Likewise, a~ also governs the constant4/line attraction to a specified ~/jline. The effects of the coefficients b~iJ and b~ij are to attract constant-~ and 41 lines to a specifiedpoint (~, ~/). The coefficients, c~i, c~, d¢,.j and d,,j function as decay factors which determine the range of the concentration effect of the and ~/lines, respectively. The elliptic grid generator described above is applied to the 4:1 planar contraction with a sharp re-entrant corner to generate a 55 by 13 mesh, shown in Fig. 1. Transformations (3) and (4) are solved using the SOR method [ 12 ] subject to the Dirichlet-type boundary conditions that consist of the coordinate values of the boundary nodes describing the flow domain. The lengths of the upstream and downstream sections are taken as 10 and 20, respectively. The curvilinear coordinate lines shown in this mesh are attracted to ~i= 29, ~/j= 13, respectively, through the grid control functions with a¢~= 70, c~ =0.4, b~iJ =0, a~ = 15, c~ =0.4, and b,i.j =0. The next step of the numerical mapping involves the transformations of the various derivatives of the dependent variables in the physical domain to those in the computational domain. This operation, involving application of the chain rule, yields [12,8] qx = [Y~q¢ --Ycq~ ] / J qy = [ - - x . q ¢ + x c q . ] / J qxx = axxq¢ + bx~q. + Cxxq¢¢+ dx~q~. + e~xq~. qyy = anq¢ + byy q. + cyyq¢¢ + d:~q.~ + eyyq¢. qxy = a~yq¢ + bxyqn + Cxyq~¢ -1-dxyq.. "b exyq¢~ Substitution of these expressions into eqn. (2) gives, D q t -[-t~q¢ -{-B q , + Cq = (~q¢¢ -{-A q , , + .#q¢,,

(5)

where D=D

I0 Fig. 1. M e s h on a 4:1 abrupt contraction.

30

358

B = ~xA + tl#B - ( bxx + byy ) G

C=C G=(Cxx+c~)G I~= (dxx +d:~)G I~= - (exx +ery)G and the transformation coefficients J, a~, a~, axy, b~x, b~, b~y, Cx~,cry, Cxy, dxx, dry, dxy, exx, eyy,e~y are obtained after mathematical re-arrangements [ 12 ]. In this study the staggered mesh proposed by Harlow and Welch [14] is adopted. This mesh was used first by Gatski and Lumley [15], and later by Zick [ 16 ] in simulations of viscoelastic flows with the Oldroyd model using a streamline-vorticity formulation. In the staggered mesh Pij, aij and 7ij-denote the p, a and ~ at the centre of the cell, respectively, uij and v/j denote the u and v at the right boundary and the upper boundary of the cell, respectively, and zij denotes the ~ at the upper-right corner of the cell. When the elliptic mesh generation is applied to the staggered mesh, four sets of transformation coefficients (J, a~, a~y, axy, b~x, bye, bxy, Cxx, c:~, c~y, dx~, d:~, d~y, e~x, eyy, exy) are obtained for the centre, right boundary, upper boundary and upper-right corner of the cell, respectively.

2.2. A D I upwind scheme The factored implicit scheme (ADI) proposed by Douglas and Gunn [17] has been used successfully by researchers such as Beam and Warming [18] in computations of compressible flows. Using this scheme system (5) is factored as the following

Dqt + (A¢ +A~ )q+ Cq=I~q¢, 7

(6)

where the differential operations A¢ and A, are defined as

A.=[~2-A 02] Introducing an intermediate step, denoted by *, the time discretization of the system (6) can be written in the form of a two-step scheme as follows: q*-qn+D

-

q~ +I --q* + D

l~t(A ~qn + A ~ q * +

Cq*) = D

- ~at (Pq~.)~

- l~t(A~q~+ I+ A ~q* + (~qn + ~) = D

- ~At (Pq~.)n

359

where the first and the second expressions represent a prediction and a correction step, respectively. The above scheme can be further expressed in a simpler A-form formulation [18]as [D+ At(A; + $!)]Aq*=

- At(A;” +A:

+C)q” -At(Fq&’

[D+ At(Ay + $C)]Aqn+ l=DAq*

(7) (8)

where Aq* =q*_qnandAqn+l=qn+l_qn For most numerical schemes used in viscoelastic flow simulations, upwind schemes are used to stabilize the numerical simulations of viscoelastic flows [ 7,8,15,16,19-241. The upwind difference scheme as proposed by Peyret [ 251 is applied to the convective terms A, and A,, of (7,8) with the notations ’ and ‘. Using Peyret’s upwind scheme the derivative terms denoted as (a/a,) ’ in the term AT, and (a/a,) ” in the term At represent the following operations [8,25 19

(-g=[

wJ(;)++

(f>“=[

~1+%&$‘+

(Ifs,)(;)-] w(~)-]

Etl=sign V, where the superscript + denotes the forward differencing and the superscript - denotes the backward differencing. In this work, V is the projection of the velocity vector along the q-axis, that is, V= ux,,+ uyv The same procedure is applied to the term A ;” with cC= sign U and U= UX,+ uye Scheme (7,8) is almost fully implicit except for the cross derivative term pq&. This can be seen by eliminating the intermediate term Aq* in the scheme (7,8), which yields q n+l-qn+D-lAt(A~+A~+e)qn+l=o-lAt(~q~~))” Scheme (7,8) is also coupled for the unknowns (Ap*, Au*, Au*, AcJ*, AT*, Ay*) and (Apn+‘, Aun+l, Aun+l, ArY'+l,A7"+l, Ar”“), so that the iterative procedure requires an inversion of a 6 x 6 block tridiagonal matrix of a size of MxN at each iteration [ 26,8]. The simulation result shown in Fig. 2 indicates that the axial velocity increases drastically in the contraction region, which induces large stress gradients in this region. The profile of the normal stress Qalong the center line of the flow domain is shown in Fig. 3. The stress gradient increases with the Deborah number, as does the stress peak level; therefore, the simulation accuracy is more difficult to maintain at large D,, as explained in later sections. A very large normal stress gradient is found along the horizontal line y= 1, as

360 6! i

De = 0.00 De = 0.70

/f ~ \

De = 2.45

f

De = 6.60

4

\ ~.

I - ......

I . . . . . .

~ --

r.,

]:' 3

2~-

1

0

0

I

~

I

I

I

I

5

10

15

20

25

30

35

X Fig. 2. A x i a l v e l o c i t y a l o n g t h e c e n t e r line.

De = 0.00

i \ f t

De = 0.70 De = 2.45

P: f:, / /::

4

t/

De = 6.60

\ \ \ \\ \

•

-2

0

NN

I

t

I

I

I

I

5

10

15

20

25

30

X Fig. 3. N o r m a l s t r e s s a a l o n g t h e c e n t e r line.

35

361 140 De =0.00

120

De = 0.70 1

100

........

i

De = 2.45 De =6.60

8O U

6O 4O

........"

20 0 -20

I

0

5

10

i

I

I

]

15

20

25

30

X Fig. 4. Normal stress a along the line y = I.

25OO

2O00

1500

1000

5O0

0

I

I

I

I

I

I

I

1

2

3

4

5

6

7

De Fig. 5. Normal stress gradient at the grid near the sharp corner on the downstream boundary. s h o w n in Fig. 4. T h e s m a l l n o r m a l s t r e s s i n c r e a s e s s u d d e n l y t o a m u c h h i g h e r level a t t h e c o n t r a c t i o n a n d a s h a r p p e a k exists a t t h e r e - e n t r a n t c o r n e r . T h e level of s u c h a n o r m a l s t r e s s g r a d i e n t is p l o t t e d in Fig. 5 as t h e f u n c t i o n o f De

35

362

which is verified analytically in later sections. 2.3. Filtering stabilization method

If a solution is represented by a Fourier series, it can be shown that the damping of solutions by the upwind scheme depends on the wave-number of each component of the Fourier series [9]. It is further shown [9] that the upwind scheme severely damages the solution components of the lower wavenumbers which are the major part of the solution. To avoid the excessive damping of the lower wave-number components and hence increase the accuracy of the numerical solution, a filtering stabilization method of 8th order is proposed, as follows,

qin,+l__.n±l/2 (1)16 -- -~ ,j ~

(~

16

16

n+l/2

+ 8,1 ) (qi,j

)+ (½

)32

[~16~16"1[__n+1/2 ~ v ~ url ] ~ t l i , j

)

where qinJ- 1/2 is the undamped solution obtained from schemes without using a stabilization method, and 5~n (qij) is the central-difference operator corresponding to the 2nth-order derivative qi,j" 2n For example, 6(qj) = q j + l - q j - 1 . It is found that the above filtering method preserves the solution components of lower wave-numbers so that the accuracy of the solution is enhanced. The detailed damping characteristics of the filter at various wave-numbers can be found in reference [9 ]. 3. Large stress gradient in the contraction region

Several solution characteristics contribute to the large stress gradient in the contraction region. Firstly, the stresses are singular at the sharp re-entrant corner, which cause significant simulation inaccuracy and discretization error which induce the breakdown of the simulations. Lipscomb et al. [27] obtained the explicit form of singularity at the sharp corner for the second-order fluid (SOF), which shows that the singularity of the shear stress for the SOF is one order higher than that for Newtonian fluid, and the difference of the two stresses is proportional to We, that is, ~= r - °'4~55/1(0) + Wer -0"911°/2(0) where r denotes the distance from the sharp corner, and 0 is the angular coordinate in the polar coordinate system with its origin at the sharp corner. However, for a Maxwell-type fluid, no explicit form of singularity has been derived. A perturbation analysis on the Giesekus fluid can only show that, for small Weissenberg numbers, there exists a critical We beyond which the mesh refinement computation fails in the vicinity of the sharp corner of a size 0 ~
363

j.

,

x

Fig. 6. Sharp comerof angle a=3n/2. Secondly, the mold geometry of the abrupt contraction can also induce a large stress gradient. Suppose that point 0 denotes the sharp corner, and I denotes a grid {e.g. due to a finite difference) on the downstream boundary from the corner, with x i - X o = Ax, as shown in Fig. 6. On this section of the boundary, u = v = ux = vx= 0. Therefore, the constitutive model, the upper convected Maxwell model (UCM), in eqn. (1) at I in the Cartesian coordinate system (x, y) becomes, 2 ai-2WeUy,l,

TI=Uy,I ,

~I=0

(9)

Similarly, in the vicinity of the corner, the grid J is on the upstream boundary from the corner with y j - y o = A y , and u=v=uy=vy=O. Hence, the constitutive model in eqn. (1) at J gives, aj-~0,

Tj-'Vx,j,

~j=2Wev2,j

(10)

At small We, (Vu)I=O(VUNewtonian) I and (VV)j=O(VVNewtonian)J; thus (9) and (10) become,

al=WeO(1/(Ax)2n),

aj.~.O

~, = 0 ( 1 / ( A x ) n ) ,

T j = O ( 1 / ( Ay) n )

71--~-0,

~ j = W e O ( 1 / ( ~ ) 2n)

where n=0.455 for Newtonian flow [27,28]. At small We, the above relations show that the singularity of the UCM model, in the directions of two boundaries, are of the same of order as that of the SOF. It can be easily proved that the lower convected Maxwell model has the same properties. The stress gradients between I and J (in eqns. (9) and (10)) are of the following orders of magnitude,

364

v a = ( ai - a j ) /,~=O ( Wo ( Vu )2/,~)

Vr( T,I- rj) /~=O(Au/,~)

(11)

VT= (TI--TJ)/(~=O( We(VU)2/{~) Relation (11 ) shows that as a mesh refinement proceeds for the UCM model, an extra difficulty is introduced to numerical simulations due to the large gradients of normal stresses, a and 7, between grid I and J, and such large gradients are proportional to the Weissenberg number. It is expected that the large stress gradients can extend from I and J on the boundary to the rest of region around the corner. Such large stress gradients in the numerical simulations using the ~UCM model may cause uncontrollable discretization error [27,2932 ], as well as intensify the non-linear interaction between velocity field and stress field [8,9]. 4. Solution a c c u r a c y The solution of the above contraction flow shows that the normal stress gradient in the contraction region is large, and its level is approximately proportional to the Deborah number. Such a large stress gradient makes the upwind differencing different from the central differencing in their magnitudes. Consider a Maxwell/Oldroyd constitutive equation at the center line of a symmetric flow domain (such as the 4:1 contraction domain) [8,32], where v=O:

da

~+

1-

2Deux

2Ux

Dett a-De u

(12)

which has the following modified equation if the first-order upwind scheme is used [8,26],

da 1 - 2Deux

~-~

Deu

2ux . d2a ff=~eU"b AX dx2

(13)

The magnitude of the artificial diffusivity term, Axaxx, can be evaluated based on the original stress field (the one without using upwind scheme). The gradient of original stress field is given by the original constitutive equation (12), i.e.,

da 2Deux-1 dx Deu

2ux Deu

Thus, the second-order derivative in the artificialdiffusivityterm is written as,

365 d2a

u(2Deuxxa+(2Deux-1)a~ +2uxx) Deu 2

dx 2 -

Ux((2D~u,,-1)a+ 2u~) DeU 2

The relation between the artificial diffusivity term and other terms in the modified constitutive equation (12) is, d2a_] Ax-~J=Ax(2D, ux-l)[d-d-d-~]

(2uDeuxx-ux( 2D~u~- 1 ) )Deu [1- 2D~u~al -2~u~- I [ D~u +Ax2(UUxx-U2x)[2Ux] hX

2ux

[_D-~uJ

(14)

where "[ ]" on the right-hand side of relation (14) denotes the terms appearing in the modified constitutive equation (12). Relation (14) shows that the magnitude of each term in eqn. (12) depends on the velocity field, i.e., Uxand u~. The flow in an abrupt contraction is taken as an example for the above magnitude analysis. In such a contraction flow, u~ and u~x experience drastic changes from the entrance, Xl, to the exit of the contraction region, x2, as shown in Fig. 7. Since the length of the contraction region is small, i.e., x2- Xl = O (Ax), and the difference between the velocities inside and outside the contraction region is finite (assume a step change in velocity field in the contraction region [8] ), the following relations hold [8]:

(du)

=0 (~X--1)'

XCXl,X2

=O(Ax-2)

kd~2]X~Xl,X2

Thus, the magnitude of the artificial diffusivity becomes

I d2a-]

a])+ O(F2Ux]~ Deu

\LDeuJ./

The above relation shows that the artificial diffusivity of the first-order upwind scheme has the same order of magnitude as the rest of the terms in the constitutive equation in the contraction region, and it does not vanish as Ax-* 0. Therefore, the first-order upwind scheme can cause a significant loss of simulation accuracy in the region of the abrupt contraction. A direct analysis of the magnitude of the simulation inaccuracy of the upwind scheme is difficult. Such an analysis requires the solution of the modified equation (13) in order to evaluate its deviation from the solution of the original equation (12). However, a qualitative characterization of the inaccuracy can

366

/

U1 I

u2

I

X 1

X2

Fig. 7. Velocitycharacteristic alongthe center line of the abrupt contraction. be achieved by a spectral analysis on various levels of stress gradients. Consider a stress solution of a large gradient, x (n), where n denotes a grid point in the physical space. The discrete Fourier transform of x (n), X ( m ) , is given by the following 1-D discrete Fourier transform pair: N--1

X ( m + l ) = ~ x(n+l)e j(2~mn/N)

(15.1)

n~O

1N--1

x(n+ l ) = ~ ~__oX(m+l )e-J(2~mn/N)

(15.2)

where m denotes wave-number, and N denotes the number of mesh points in domain discretization. The solution of a large stress gradient (shown as the original signal in Fig. 8) is represented by Whittaker's cardinal function [33],

x(n) =

sin(n~Ax/h) n~Ax/h

and its Fourier transform in wave-number space, obtained using relation ( 15.1 ), is shown in Fig. 9, i.e., X(m)=

h,

-~/h~m<~/h

0,

elsewhere

where h is a parameter controlling the gradient and peak level of the stress signal. If the upwind scheme (13) is used, the damping of the Fourier components of the solution of a simple finite difference scheme [9] is given as IG[ = [ (1-De/2+Ddx c o s f l - (1-2Deux)At)2+ (De~2sinfl)2] 1/2 which is plotted in Fig. 9 as the damped signal, where G=an+l/an, /2= uAt/Ax, the wave-number fl= 2nrn/N, and At is an iteration parameter. This curve shows the Fourier transform of the solution after the damping

367

original si'gnal ~

damped signal

x(n)

0

h n

Fig. 8. The original physical signal and the damped physical signal due to the use of the upwind scheme.

X(m) original signal

i

0

g/h

m

Fig. 9. The Fourier transform of the original signal and the damping curve of the upwind scheme.

368 due to the upwind scheme. The corresponding damped solution in the physical domain (obtained using relation (15.2)) is shown in Fig. 8 as the damped signal. It is noted that the original solution in the form of a peak is greatly damaged especially in its magnitude and gradient. This analysis may explain some existing simulation results using upwind schemes [20,32] where the magnitude of the stress peak at the contraction section becomes small at large

De. The advantage of the filtering method in achieving a higher accuracy can be seen in the following. As the order of the filter increases (from 1 to 4), the width of the undamped region increases (as shown in Fig. 11 ). The damping effect ofa 1-D filter [9] is IG I = 1 -

sin2nkA~/2

----1 -- sin2nfl/2

which is plotted as the darted line in Fig. 11 for n = 1 and 4. The corresponding filtered solutions in the physical domain are plotted in Fig. 10. It shows that as the order of the filter, n, increases, the peak solution is damaged less; thus, the numerical solution using the higher-order filter is more accurate than that using the first-order upwind scheme. In general, the loss of the simulation accuracy decreases in the region of a smaller gradient. That is, as the parameter h in Whittaker's cardinal function increases, the gradients of the physical solution in Fig. 8 and Fig. 10 decrease; thus, the Fourier transform of the original solution shrinks towards the region

x(n) ~figinal signal

I

2~ order = 4

o

h

n

Fig. 10. The original physical signal and the damped physical signal due to the use of high-order filtering.

369 X(m) original signal

o

~/h

ITl

Fig. 11. The Fourier transform of the original signal and the damping curve for high-order filtering.

of lower wave-numbers, as shown in Figs. 9 and 11, and few Fourier components are damped by the upwind scheme and the filtering method. 5. C o n c l u s i o n s

To ensure the validity of simulation results of a liquid molding process involving a drastic mold boundary change, the solution accuracy that a software yields must be understood and characterized. This work showed that loss of accuracy occurs mainly in the region where the mold boundary changes abruptly. The test case, the Maxwell/Oldroyd-B fluid flow through a 4:1 contraction domain, shows that the large stress gradients formed at the contraction section and the sharp re-entrant corner are the major sources of inaccuracy in numerical solutions. A qualitative estimate of the inaccuracy of the numerical simulations using upwind schemes is achieved. Firstly, an extra term, the artificial diffusivity, is found to be introduced to the constitutive model by the first-order upwind scheme, which is of the same of order of magnitude as other terms in the constitutive model. Secondly, the damping of the upwind scheme on the stress gradients of various levels is characterized through a spectral analysis, which shows that large gradients at the contraction section and the re-entrant corner are damaged severely by the upwind scheme. It is found also that a new filtering technique can improve the accuracy of the simulation results by preserving the solution components of lower wave-numbers. In addi-

370 tion, the large stress gradients at the re-entrant corner are characterized to reveal the influence of the Deborah number on the simulation accuracy.

References 1 R.B. Bird, R.C. Armstrong and O. Hassager, Dynamics of Polymeric Liquids, Vol. 1, Wiley, New York, 1987. 2 S. Subbiah, D.I. Trafford and S.I. Gii~eri, Non-isothermal flow of polymers into two-dimensional, thin cavity molds: A numerical grid generation approach, Int. Heat Mass Transfer, 32 (1989) 415-434. 3 J. Coulter and S.I. Gtiqeri, Resin impregnation during the manufacturing of composite materials subject to prescribed injection rate, J. Reinforced Plast. Compos., 7 (1988) 200-219. 4 M.J. Crochet, A.R. Davies and.K. Walters, Numerical Simulation of Non-Newtonian Flow, Elsevier, Amsterdam, 1984. 5 R. Keunings, Simulation of viscoelastic fluid flows in: C.L. Tucker III (Ed.), Fundamentals of Computer Modeling for Polymer Processing, Carl Hanser, New York, 1989. 6 S.A. White, A.D. Gotsis and D.G. Baird, Review of the entry flow problem: Experimental and numerical, J. Non-Newton. Fluid Mech., 24 (1987) 121-160. 7 R.A. Brown, R.C. Armstrong, A.N. Beris and P-W. Yeh, Galerkin finite element analysis of complex viscoelastic flows, Comput. Methods Appl. Mech. Eng., 58 (1986) 201-226. 8 H. Gong, An analysis of viscoelastic flow simulations, Ph.D. Thesis, University of Delaware, 1990. 9 H. Gong and S.I. Giiceri, Non-linear computational instability, damping and filtering in viscoelastic flow simulations, J. Non-Newton. Fluid Mech., 38 (1991) 247-271. 10 A.J. Chorin, A numerical method for solving incompressible viscous flow problems, J. Comput. Phys., 2 (1967) 12-26. 11 H. Gong and S.I. Giiqeri, An investigation of an inherent linear instability of an iterative scheme in viscoelastic flow simulation, J. Non-Newton. Fluid Mech., submitted. 12 S.I. G0qeri, Finite difference solution of fluid problems, in C.L. Tucker III (Ed.), Fundamentals of Computer Modeling for Polymer Processing, Carl Hanser, New York, 1989. 13 D.L. Trafford, S.I. Gfiqeri and M.J. Crochet, A numerical grid generation approach to the injection molding problem, in: NUMIFORM--Int. Conf. Numerical Methods in Industrial Forming, Gotenburg, Sweden, 1988. 14 F.H. Harlow and J.E. Welch, Numerical calculation of time-dependent viscous incompressible flow with free surface, Phys. Fluids, 8(12) {1965). 15 T.B. Gatski and J.L. Lumley, Steady flow of a non-Newtonian fluid through a contraction, J. Comput. Phys., 27 (1978) 42-70. 16 A.A. Zick, Numerical simulation of non-Newtonian fluid flow, Ph.D. Thesis, Stanford University, 1983. 17 J. Douglas and J.E. Gunn, A general formulation of alternating direction methods--Part I. Parabolic and hyperbolic problems, Numer. Math., 6 (1964) 428-453. 18 R.M. Beam and R.F. Warming, An implicit factored scheme for the compressible NavierStokes equations, AIAA J., 16 (1978) 393-401. 19 J.M. Marchal and M.J. Crochet, A new mixed finite element for calculating viscoelastic flow, J. Non-Newton. Fluid Mech., 26 {1987) 77-114. 20 B. Debbaut, J.M. Marchal and M.J. Crochet, Numerical simulation of highly viscoelastic flows through an abrupt contraction, J. Non-Newton. Fluid Mech., 29 (1988) 119-146.

371 21

22 23 24

25 26 27 28 29 30

31 32 33

J.H. Song and J.Y. Yoo, Numerical simulation of viscoelastic flow through a sudden contraction using a type dependent difference method, J. Non-New ton. Fluid Mech., 24 ( 1987 ) 221243. M. Fortin and A. Fortin, A new approach for the FEM simulation of viscoelastic flows, J. Non-Newton. Fluid Mech., 32 (1989) 295-310. S. Pilitsis, Calculations of steady-state viscoelastic flow in a periodically constricted tube, Ph.D. Thesis, University of Delaware, 1990. R.C. King, M.R. Apelian, R.C. Armstrong and R.A. Brown, Numerically stable finite element techniques for viscoelastic calculations in smooth and singular geometries, J. Non-Newton. Fluid Mech., 29 (1988) 147-216. R. Peyret and T.D. Taylor, Computational Methods for Fluid Flow, Springer, New York, 1983. D.A. Anderson, J.C. TannehiU and R.H. Pletcher, Computational Fluid Mechanics and Heat Transfer, Hemisphere, Washington, DC, 1984. G.G. Lipscomb, R. Keunings and M.M. Denn, Implications of boundary singularities in complex geometries, J. Non-Newton. Fluid Mech., 24 (1987) 85-96. H.K. Moffatt, Viscous and resistive eddies near a sharp comer, J. Fluid Mech., 18 (1964) 118. R. Keunings, On the high Weissenberg number problem, J. Non-Newton. Fluid Mech., 20 (1986) 209-226. M.A. Mendelson, P.W. Yeh, R.A. Brown and R.C. Armstrong, Approximation error in finite element calculation of viscoelastic fluid flows, J. Non-Newton. Fluid Mech., 10 (1982) 3154. J.R. Rosenberg and R. Keunings, Further results on the flow of a Maxwell fluid through an abrupt contraction, J. Non-Newton. Fluid Mech., 29 (1988) 295-302. S.L. Josse and B.A. Finlayson, Reflections on the numerical viscoelastic flow problem, J. Non-Newton. Fluid Mech., 16 (1984) 13-36. R. Vichnevesky and J.B. Bowles, Fourier Analysis of Numerical Approximation to Hyperbolic Equations, SIAM, Philadelphia, 1982.