- Email: [email protected]
A purely hyperbolic model for unsteady viscoelastic flow
Journal of Non-Newtonian Fluid Mechanics, 32 (1989) 197-224 Elsevier Science Publishers B.V., Amsterdam - Printed in The Netherlands
A PURELY HYPERBOLIC VISCOELASTIC FLOW
197
MODEL FOR UNSTEADY
F.R. PHELAN, Jr., M.F. MALONE and H.H. WINTER Department of Chemical Engineering, Goessmann Luboratoly, University of Massachusetts Amherst, MA 01003 (U.S.A.) (Received June 27,1988;
in revised form December 22,1988)
Summary A purely hyperbolic model for unsteady, two-dimensional viscoelastic flow is developed. The model uses a pressure-density relationship in the continuity equation to obtain a dynamic equation for the pressure, and the Giesekus constitutive equation with zero retardation time to describe the extra stress. The small degree of compressibility allowed by the model is negligible for practical purposes but important mathematically because it changes the unsteady governing equations from parabolic-hyperbolic to a purely hyperbolic system. This allows the application of a numerical method designed specifically for solving hyperbolic problems, and thus, to better treat the hyperbolic nature of the flow which is dominant at high Deborah numbers. The method is used to predict the velocity and stress fields in a square cavity. The results for the creeping, Newtonian limit compare well with previously published calculations. Viscoelastic flow computations are presented for Deborah numbers up to 4. Comparison of the results obtained on several grids show that the calculations increase in smoothness with mesh refinement, and that the scheme is numerically stable. The calculations predict flows exhibiting significant elastic effects. A transition from single vortex flow to flow containing multiple vortices, occurs under creeping flow conditions for Deborah numbers greater than 1.
1. Introduction The numerical simulation of viscoelastic flow has been an active area of research in computational fluid mechanics for over a decade. The persistent 0377-0257/89/$03.50
0 1989 Elsevier Science Publishers B.V.
198 failure of many numerical schemes when elastic effects are significant has been an outstanding problem and is well documented in the literature [1,2]. Various investigators have laid blame for this problem on the use of fundamentally unsound constitutive equations, failure to properly treat stress singularities, improper specification of stress boundary conditions, excessive approximation error due to an inadequate ability to resolve the large gradients which may occur in these flows, or bifurcation from two-dimensional flow to three-dimensional or time-dependent flow. While all of these are probably problematic to some degree, recent studies indicate that the root of the numerical difficulties is associated with the use of numerical schemes that are inappropriate for the mathematical problem. Such schemes give rise to fictitious limit points that are not an actual property of the mathematical system but are due solely to the numerical approximation procedure [3,4]. Joseph et al. [5] were the first to suggest that the complicated mathematical nature of the equations that are commonly used to simulate viscoelastic flows may be related to the numerical difficulties. They classified the equations governing the steady, two-dimensional flow of the upper convected Maxwell fhtid, as well as those governing a class of Oldroyd models of which the Maxwell fluid is a special case. They showed that a change in type from elliptic to hyperbolic occurs in the equation for the vorticity when certain conditions depending on the local velocity and stresses fields and on the material parameters of the fluid are satisfied. Joseph and Saut [6] have classified the equations governing the unsteady flow of viscoelastic fluids with constitutive equations of the form
(1) where S/St represents the upper-convected derivative. They have shown that a loss of evolution (Hadamard instability) occurs when certain critical values of the normal stress are reached. Dupret and Marchal [7] have found that a Johnson-Segahnan Fluid exhibits the same behavior. In some problems, a change of type in a model is reflected in the actual flow and is a stable physical phenomenon, e.g., the onset of shocks in transonic gas flow. However, the analysis of Joseph and Saut [6] shows that a change of type in viscoelastic flow problems is not always stable. In problems without inertia, the point at which a change of type occurs in the steady problem is exactly the same as that at which the corresponding unsteady problem becomes non-evolutionary. Thus, a steady-state numerical simulation in which a change of type occurs is sometimes susceptible to Hadamard instabilities which grow in severity on more refined computational meshes. Several investigators [8,9] have linked the occurence of a
199 change of type in their numerical schemes to a subsequent loss of convergence. While the use of constitutive equations which change type and become non-evolutionary is certainly problematic, it does not entirely explain the problems that have limited past computations. Renardy [lo] has shown that a change of type cannot occur for the Upper Convected Maxwell Fluid at zero Reynolds number. This contradicts the results of various numerical studies which have shown that a change of type in the flow an Upper Convected Maxwell Fluid at zero Reynolds number led to a loss of convergence [ll]. This contradiction suggests that the loss of convergence in traditional numerical schemes is due to errors which lead to an artificial change of type (and thus, a non-evolutionary problem) and that numerical simulations can be improved by developing schemes which properly treat the mathematical type of the system. Such schemes could not be expected to obtain solutions when the mathematical problem really loses evolution, for in such cases the problem is ill-posed. However, they would provide a means for preventing numerically induced instabilities in problems that are actually well-posed. The studies of Joseph and co-workers have emphasized the importance of recognizing the type of the system in models for viscoelastic flow. One of the more important results that has emerged is that a wide range of viscoelastic models are mixed systems, being hyperbolic-elliptic in the steady-state and hyperbolic-parabolic in the unsteady-state. This situation is difficult to treat numerically because most algorithms are designed for systems of equations of a single type and methods for mixed systems are not well developed. However, numerical methods for purely hyperbolic systems are available, primarily from studies of compressible flow. LeBlanc [12] formulated a one-dimensional model problem for viscoelastic flows by ignoring the conservation of mass constraint. Numerical treatments of this simple problem displayed many of the essential features of more realistic models and the method of characteristics provided the only satisfactory solutions. This simplified approach provided the starting point for this study. In what follows, a purely hyperbolic model for unsteady two-dimensional flows is developed. Purely hyperbolic behavior is obtained by allowing a small degree of compressibility in the flow. With the incompressible constraint relaxed, the continuity of mass condition provides an evolution equation for the pressure. This equation together with the momentum balance and the Giesekus equation with zero retardation time as the constitutive model (or any expression that can be written in the form of eqn. (1)) constitutes a completely hyperbolic system in time and space provided certain critical normal stress limits which depend upon the Deborah number
200 are not violated. Outside of these limits the model becomes non-evolutionary and the problem is not well-posed. Although this does reflect deficiencies in the constitutive model, these are not particularly relevant to the development of a numerical scheme. Numerical solutions are developed using an algorithm designed for solving systems of hyperbolic equations that arise in gas dynamics. As an example, the flow in a driven cavity is examined. 2. Model formulation 2. I. Basic equations The model employed her consists of the Cauchy momentum balance, the Giesekus [13] constitutive equation and the equation of continuity of mass. In the absence of body forces and for a zero retardation time, these are written as Do pa--V.7+(1/3)V(trr)+VP=O, h
DT i~-VVT*T-T*vv
(2)
1-qq=
-T*(I+(ah/~)T),
where P is formulated here as the mean normal stress which accounts for the extra term involving the gradient of the trace of the stress in the momentum balance (tr 7). A model based on the assumption of incompressible fluid behavior inevitably results in a system of equations of mixed type. In order to obtain a purely hyperbolic system, nearly incompressible fluids are considered as a perturbation of the incompressible state. A first-order Taylor expansion for the density (p) at constant temperature together with the definition of the isothermal compressibility gives P=Po[I
+ (P-P,)/819
where /3 is the inverse of the isothermal compressibility. into eqn. (4) yields the expression
(5) Inserting eqn. (5)
+g+[p+(P-P,)I(v*v)=O. One further simplification is made in this last expression. Because a pertubation of the incompressible state is desired, only cases for which ((P PO)//?1 -=K1 are of interest. Thus, eqn. (6) is simplified to DP ~+/3(v~v)=o.
201 Equations (2), (3) and (7) make up the model to be studied here. These can be written in dimensionless form as E-v*7+(1/3)v(trT)+VP=O, RDU D
(8) -+=
-7=(1+&T),
(9)
(10) where R = (pVL/p) is the Reynolds number, D = (XV/L) is the Deborah number and u, T, P and /3 are the dimensionless velocity, extra stress, pressure and compressibility. For a two-dimensional flow where ?J= [u, VI,
(11)
and
S
‘=Q
Q
I I
(12)
T’
the components nates as
of eqns. (8), (9) and (10) are written in Cartesian coordi-
tLJ+ UU, + Vu, + R-‘(P,
- (2/3)S,
c: + Uv, + WY + R-I( Py + (l/3)$ St+ Us,+
vs,-Au,-2Qu,=
uq+
Pt+ up,+
(13a)
- Q, - (2/3)T,)
= 0,
03b)
Q2)), -(Q/D+(I~Q(S+T)),
V5-2Q<-4= vPy+p(ux+
= 0,
-(S/D++‘+
Q,+UQ,+VQ,-(&2)U,-(&2)V,=
T,+
- Q,, + (1/3)T,)
-(T/D+cx(T~+Q~)), v,)=o,
034
(13d) (13e) W)
where a = 2( S + l/D), h==(T+l/D), and the subscripts t, x and y denote partial-differentiation those variables.
(14) (15) with respect to
2.2. Classification The system of equations (13), can be written in the form wt+A*w,+B*w,=f,
(16)
202 where w is a vector of length six with components, w = [U, V, S, Q, T, PIT, the matrices A and B contain the coefficients of W, and wv, and f is a vector of forcing terms. Equation (16) represents a quasi-linear system of first-order partial differential equations and in order to treat such a system numerically, we must first classify it mathematically [14]. In general, the system of equations (16) is evolutionary in (x, y, t) if for every unit vector v = [ vl, v2] the eigenvalue problem (-XI++4
+Y,B)e=O
07)
has real eigenvalues only, and hyperbolic in (x, y, t) if in addition, there are six linearly independent eigenvectors, e, [6]. However, we present here the classification information that is relevant with regard to the finite difference algorithm to be employed for solving the model equations, and that is the conditions for an evolutionary and hyperbolic problem for the cases v = [l, 0] and v = [0, 11. The information obtained for these particular cases determines separately the conditions for an evolutionary and hyperbolic in the (x-t) and (r-t) planes, respectively [14,15]. The system is evolutionary in (x, t) if the eigenvalues of A are all real. In addition, the system is also hyperbolic if A can be written in the form A =X*A,*P,
(1%
where the diagonal matrix A, contains the eigenvalues of A and X-’ is the matrix of left eigenvectors. The same can be said of the behavior in (y, t) with respect to the eigenvalues and eigenvectors of B. The eigenvalues of the matrix A are L,,
= u,
@a>
L,,
= u,
Wb)
A,,, = u+ /272x,
(1W
A,,,=
U-
Wd)
A,,,=
tiJ+/@GjpY,
(194
A,,,=
u-J@zjq?x,
ww
\ia/zR,
and the eigenvalues of the matrix B are
v, A,,, = v, A,,, =
(204 (2Ob)
203
A,,, =
v-t \lqiiF,
A,,, =
v-
@x,
(2Od)
A,,, =
v-
/@ETjjp.
@Of)
(W
The compressibility parameter, 8, takes on positive values only, and thus, an inspection of eqns. (19) and (20) shows that A and B have real eigenvalues if * S>
-I/D
(21)
-I/D.
(22)
and T>
When the inequalities (21) and (22) are satisfied, the coefficient matrices A and B each have six linearly independent eigenvectors and they can be decomposed by the similarity transforms A =X4,*X-f, B=
Y*A;
(23) Y-r.
(24)
Under these conditions, the equations of the model are both evolutionary and hyperbolic. 2.3. Parameter limits It is not possible to establish general apriori parameter limits on (Yand D for which one of the conditions (21) or (22) will be violated. However, it is useful to examine two simple flows for which analytical solutions can be obtained to determine under what conditions a loss of evolution may occur. First, a pressure-driven shear flow in a slit with a gap width of 2 is considered. The analytical solutions for the shear and normal stresses for a Giesekus fluid with zero retardation time and non-zero (Y are given by Q=
-Y,
(25)
S = (1,‘2aD)(
- 1+ I/g(Y) + &xD’YU,),
(26)
T= (1/‘2aD)(
- I+ /g(y)),
(27)
* Note added in proof: M. Renardy [22] has pointed out that conditions (21) and (22) are necessary but not sufficient to guarantee real eigenvalues in (17). These could be replaced by more general conditions involving the eigenvalues of the extra stress as discussed in [23].
204 where u,=
Y[l+(2a
g(Y)
= 1 - 4cy2D2Y2.
- l)HV]/[(2~
- 1) + \lg(y)12,
(28) (2%
The shear stress has been scaled to a value of unity at the walls of the flow and the pressure-gradient has been absorbed in the scaling. Since the function g(Y) appears as a square-root term in several places, the first condition that must be satisfied is g(1) 2 0. This requirement yields the condition aD I l/2.
(30)
The above condition tells us nothing about hyperbolicity; it is simply a requirement that needs to be satisfied in order that the dimensionless second normal stress be real-valued. One of the conditions for hyperbolicity can be violated if one of the normal stress components, S or T, becomes too large and negative. In this flow, S is always positive and T is always negative. Thus, we require that (Y and D satisfy (1/2cwD)( -1 + m)
> -l/D.
(31)
This yields the condition a( D2 + 1) I 1.
(32)
A physical interpretation of this condition can be obtained by solving eqn. (13d) for the viscosity function in steady shear flow. By doing this, it is found that q(y) = (1+ DT)/(l
+ aD(S
+ T)).
(33)
A violation of the inequality (32) corresponds to a negative viscosity which is a violation of thermodynamics. In other flow situations, condition (32) might not be conservative enough since it has been derived under the stipulation that 1Q I_ = 1. Consider the more general expression for T in a pressure-driven shear flow
(34) If the right hand side of eqn. (34) is substituted into the inequality (22), it is found that the inequality will hold provided that (Y and D satisy a(4D2Q2
+ 1) < 1.
(35)
This more conservative condition becomes necessary in the simulation of flow in a geometry where the outlet is smaller than the inlet, for example, a contraction flow. In such a problem, if the upstream shear stress is scaled
205 such that ] Q ] max= 1, then the downstream value will be much larger than one. The values of (Y and D must then be chosen such that (35) is satisfied where the value of Q in this expression is taken as ] Q ]zet. Finally, planar extensional flow is considered. The velocity field for this flow is given by u= [ix,
-iY],
(36)
and the normal stress components are given by s = (l/2&)(
- (1- 2iD)
+ {(l - 2iD)2 + 8icrD
T = (l/2&)(
- (1+ 2iD)
+ {(l + 2iD)’
),
- 8&D ).
(37) (38)
In this flow, T is negative. In order that the condition (22) not be violated, we require that (l/2&)(
- (1+ 2iD)
+ J(1+
2iD)2 - l&D)
2 -l/D,
(39
which yields the condition (Y
(40)
This condition is not restrictive at all since the Giesekus equation is only valid up to (Y= 1. Thus, somewhat surprisingly, it is independent of D. Based on these results, it is expected there may be parameter limits that must be observed in flow simulations with this model. However, these parameter limits reflect particular deficiencies in the constitutive equation that has been chosen for this study - not the specific procedure that has been developed here. 3. Numerical algorithms 3.1. The Split Coefficient Matrix (SCM) method The Split Coefficient Matrix (SCM) method developed by Chakravarthy et al. [16], is a finite difference method which uses information on signal propagation provided by the theory of characteristics to solve systems of first-order hyperbolic equations. This method is similar in philosophy to the method of characteristics, but much easier to apply to two-dimensional problems. The major advantage of the SCM method over the method of characteristics is that a conventional, rectangular finite difference grid can be employed for the calculations. In the latter method, the characteristic mesh in not known a priori and must be calculated as part of the solution.
206 Detailed discussions of the SCM method are provided by Chakravarthy et al. [16] and Anderson et al. [14]. The application of the SCM method to the system represented by eqn. (16) is relatively straightforward. First, the coefficient matrices A and B are each split into two parts. To accomplish this splitting, the similarity transforms (23) and (24) are used to rewrite (16) as W,+x*A;x-‘*W,+
Y*A,= Y%$=f.
(41)
The diagonal matrices A, and A, are then partitioned negative parts according to the formulas
into positive and
A,=Ax’+A,,
(42)
Ay=A;
(43)
+A,,
where the matrices with the ‘ + ’ superscripts contain the positive eigenvalues of the original matrices and those with the ‘ - ’ superscripts contain the negative eigenvalues. Using these last two equations, we define the splitting of the the matrices A and B as A=A++A-,
w
B=B++B-,
(45)
where the split coefficient matrices are given by A+=X.A,+.X-l,
(46)
A-=X.A,.X-l,
(47)
B+=Y*A;Y-I,
(4%
B-=Y*AJ*Y-?
(49)
Using eqns. (44-49) we can then rewrite (41) as w,+A+
lw,+A-•w,+B+*w,+B-•w,=f.
(50)
Equation (50) is still entirely equivalent to the original system of equations (16). In the second part of the SCM method procedure, finite differences are used to approximate the spatial derivatives in eqn. (50). To correctly treat the hyperbolicity of the system, backward differences are used to approximate those derivatives which are multiplied by A + and B+ and forward differences are used to approximate those derivatives which are multiplied by A - and B-. The SCM methods procedure can then be represented symbolically as w,=f-(A+*w,+
+A-•w,-
+B+vv;
where the ‘ + ’ and ‘ - ’ superscripts
+B-*I+$-),
on W, and wy denote backward
(51)
and
207 forward difference approximations, respectively. Equation (51) represents a set of six first-order, ordinary differential equations to be integrated at every internal point on the finite difference grid. 3.2. Boundary conditions At the boundaries of the computational grid, some of the characteristics point outward from the domain, some are perpendicular to the boundary and some point inward. This is illustrated for the characteristics in the x - t and y - t planes along the boundaries of the geometry shown in Fig. 1. The information propagated along the characteristics that point inward cannot be used in the integration of the equations at the boundary. Thus, the number of equations along the boundary must be reduced by the number of inward pointing characteristics and replaced by the same number of boundary conditions. To illustrate this, the procedure for the approximation of eqn. (41) along the boundary whose characteristics are shown in part a) of Fig. 1 is outlined. The normal SCM method splitting procedure is used to handle the x-direction derivative terms. Backward differences are used to approximate the y-direction derivative terms. Applying this procedure to these respective terms and rearranging, eqn. (41) is rewritten as wt+ Y*h,*Y-l*~$+
=f-a,
(52)
where a=A+
l
w,’ + A -
l
w, .
(53)
Equation (52) is then multiplied by Y-r to yield Y1
l
W,= c,
(54)
where c=
-A,
l
rl*wj+
+ r’*(
f -a).
(55)
Equations (54) are the “compatibility equations” for the y-t plane. Note that each term of the vector c is a function of only one of the eigenvalues of the matrix B. To correctly reduce the number of equations and introduce the boundary conditions, the compatibility equations containing the eigenvalues which point into the computational domain are discarded. For the case of no-slip (and in general, whenever the Reynolds number is small and n u = 0), the eigenvalues X,,, and X,,, are negative. Since that is the case being considered here, the fourth and sixth compatibility are replaced with the no-slip boundary conditions for the two components of the velocity vector. In order to introduce this into the numerical scheme in a straightforl
208 ward manner, the compatibility replaced with the equations
equations that have been thrown out are
Q=O
(W
and q=o.
(57)
The system of equations (54) are thus modified to the set New*=&
(58)
where the vector c^ is derived from c by setting the fourth and sixth components of c to zero and the matrix N is given by,
NE
0
0
Y;;' Y;-,'Y;;' Y;-,’
0
0
Y*
Y$
0
0
Y$
YG1 Ygl
YG1
r;;r
0 Y&l
0 Y;;’
0 I&’
0
0
0
010 0 0 100
Y,3'
YG1
(59)
The matrix N is then inverted to give t:he the following system of equations for the rate-of-change of the dependent variables at the boundary w,=N-‘*c^.
(60)
Following this procedure automatically yields the desired conditions U, = 0 and v = 0 without any further changes in the numerical scheme. While eqn. (60) can be used to find the pressure changes along the boundary, in the course of testing the scheme it was found that a more accurate alternative was to instead, solve for the pressure by enforcing the normal component of the momentum balance. Along a no-slip boundary this is written as ii-vP=ri~v*r.
(61)
For the boundary shown in Fig. 1, eqn. (61) reduces to P,,= -+S,,+Q,+fr,.
(62)
By applying backward differences to the derivatives in the y-direction and central differences to the derivatives in the x-direction we can derive the following expression for the value of P at time step n + 1 at each point (i, j) along this boundary Peal = P,~~‘l-
:( S~+l-
+ ;( T;;ti+l - q;pr),
S,?~~l) + (k/2h)( Q~~~j - Q~~~j) (63)
209
d’
X
Fig. 1. General structure of the characteristic y-f and b) x-t planes.
where
h = Ax and
directions along no-slip boundaries
in
the a)
k = Ay. The treatment of the other boundaries is com-
pletely analogous. There is another situation that warrants mentioning and that is the treatment of comer points where the flow boundaries meet. At these points there are inward pointing characteristics in both the x and y directions. In this case, the boundary treatment shown above cannot be used because the problem requires that the information propagated along both these sets of characteristics be removed from the equations. To accomplish this a timesplitting scheme is employed in conjunction with the normal SCM method techniques used to split the spatial terms. This procedure is illustrated for the approximation of eqn. (41) at the comer labeled c) in Fig. 1. The time-splitting technique consists of partitioning eqn. (41) into two parts as follows ;wt’~‘+x.A,.x-i.w,-
= if
(64)
=
(65)
and $l+@‘+ Y*h;
Y-‘*w,
if.
Each of these split parts are then solved using the same type of boundary procedure illustrated for non-corner boundary points. First, the split compatibility equations are written as +x-f,,,,(i)
= b9
3 pp+p
= c,
(66) (67)
where
b=+X-lf-&.X-‘.w,-
(68)
and
c=:r’f-A;r’v$-.
(69)
210 Then, the appropriate compatibility equations are discarded and replaced with the conditions U, and v = 0. For the case in question, the third and fifth compatibility equations (i.e., the third and fifth components of the vectors b and c) are discarded in both the x and y directions. This procedure yields the equations n/f. &)
= 5
(70)
jjl.
= 2,’
(71)
,;,
where 6 and 2 contain the components of the vectors b and c except for the third and fifth components which are zero and M and N are given by 0
0
x,-,l
x,-,*
x,-,’
x,-,’
0
0
x,-,l
xi1
x,-,l
x,-,l
0
0
0 x,-,’
0 xi1
0 x,-,l
0 x$
M=1O 010
N=
0
0
0
0
0
x,-,’
x,-,’
x,-,’
x,-,l
0
0
Y;;’
YG’
YE1
Y,-,l
0
0
Y2i1
Y,-,’
YG1
Y&l
0
0
0
Yi’
Y&l
Y&l
100 0
0
YG1
’
(72)
(73)
The matrices M and N are then inverted and eqns. (70) and (71) are combined to yield the following system of equations for the rate-of-change of the dependent variables at the corner ,,,,=M-l.&+N-l.c^.
(74) To solve for the pressure at the comer points the momentum balance was used. Since there is no unique normal vector at a comer the two equations used along each of the respective sides were averaged and an equation for the new value of P was derived by applying appropriate finite-difference expressions. 4. Simulation of cavity flow 4.1. Description of geometry and procedure The algorithm described above was used to calculate the velocity, stress and pressure fields for the flow in a square cavity with dimensions: - 1 I x
t l.._m_!
211
Y
x
Fig. 2. Driven cavity flow geometry.
I 1, - 1 I y I 1. This geometry is depicted in Fig. 2. The x-component of the velocity along the upper wall of the cavity was assumed to be a constant of value 1 except at the comer points where it was zero. The velocity was zero on all the other walls and comer points of the flow as well. The initial conditions used for the calculations were II = 0, T = 0 and P = 1 over the entire cavity. Roache [17] discussed the problem of setting a pressure datum in a cavity flow for dynamic calculations. Because the specification of the velocity components completely satisfies the boundary condition requirements for this problem, there is no physically acceptable point at which to set a pressure datum. The procedure used here was to pick a point in the flow at random, and subtract the value of P, at that point from all the other values of Pt in the flow. This has the effect of keeping the value of the pressure constant at this point but still allowing the correct relative changes between points to occur. Roache also suggested that for this problem the final values of pressure be corrected by the addition or subtraction of a numerical constant such that the integrals of the initial and final pressure distributions over the cavity are equal. This procedure was used here but is not absolutely necessary for fluids that behave incompressibly or nearly incompressibly because only pressure differences are important in such a problem. Several methods were tested for the integration of the discretized equation set including a fourth-order Runge-Kutta, an Adams-Bashford-Moulton predictor-corrector and a simple Euler scheme. Because the time-step involved in the calculations was usually of the order of 10e5, it was found that the additional computational expenses incurred by the use of the higher order methods was not worth the very small gains in accuracy, especially, since all steady-state results are independent of the integrator employed. Thus, the majority of the calculations were performed using Euler’s method. Irrespective of the integrator used, the following time step limitations were observed to provide stability
At < MWmx>@xAY/@x+ AY ))
(75)
212 for r < 1, and At < 1/(2r.
A,,)(
AX Ay/( AZ + Ay))
(76)
for r > 1, where
(77) and r=/i@.
(78)
The maximum value of the divergence in the flow is controlled by the value of the parameter j3. Fairly large values of B are needed to insure a nearly incompressible flow. For this particular problem, values between 100 and 10,000 were satisfactory. The differences that one observes between a calculation at a low value of /3 and one at a high value of p are generally very small at a macroscopic level but not necessarily at a microscopic level. For example, two sets of calculations with constant values of R, D and (Y, but different values of B typically predict the same vortex locations and overall flow pattern, but small differences with regard to the exact shape of the streamlines may be observed. 4.2. Newtonian limit Initial calculations were performed with R = D = (Y= 0.001, and fl= 10,000, on a 17 X 17 grid. This represents a creeping, Newtonian flow limit. These calculations were undertaken to test the accuracy of the scheme by comparing with other Newtonian calculations of the same flow that are available in the literature. It should be noted that the scheme is not suited to handle the absolute limits of R and D = 0. Changing the scaling from a Deborah number to an Elasticity number [18] would allow one to accomodate the limit R = 0; however, under no condition is a zero value of either of these elasticity scaling choices acceptable since this would change the nature of the problem. The streamlines for D = 0.001 are shown in Fig. 3. The streamline values were normalized to a value of 1 along all four boundaries of the flow. The calculations show a vortex located at the point x = 0 and y = 1.5, and a difference between the streamline value at the vortex and the wall of -0.203. This is in good agreement with the results of Pan and Acrivos [19] and other investigators who calculated a vortex center of x = 0.0 and y = 1.52 and a streamline difference of -0.2 for the creeping flow of a Newtonian fluid. The small differences that exist between the two calculations are probably attributable to a number of factors such as the non-zero Reynolds and Deborah number, and finite compressibility used in our
213
L
Fig. 3. Steady-state streamlines for creeping, Newtonian flow.
Fig. 4. Steady-state pressure field for creeping, Newtonian flow.
214
Fig. 5. Steady-state first principal stress difference distribution
for creeping, Newtonian
flow.
calculations, and also possibly to the fact that the grid used by Pan and Acrivos was much more refined than the grid used here. The pressure field and first principal stress difference calculated for this flow are shown in Figs. 4 and 5. 4.3. Viscoehstic j7ow Calculating the start-up flow for a viscoelastic fluid from rest as was done for the Newtonian simulations turned out to be very time-consuming. Therefore, to save on CPU-time it was decided to use instead the steady-state results from the simulation for one Deborah number as input for a simulation at a higher Deborah number, but with the numerical value of the plate speed increased. Thus, the plate speed was used as a parameter to increase the Deborah number. Viscoelastic flow computations were performed with R -C 0.001, a = 0.25, /3= 500 for Deborah numbers from 1 to 4 on 9 X 9 and 17 X 17 grids. Calculations were considered to be steady when the changes in the values of the flow variables were of order 10v6 per 100 time steps and the integral drag and normal forces along the upper wall of the flow had assumed steady values. The steady-state streamlines for D = 1, 2, 3 and 4 for both grids are
a)
(cl
Cd)
Fig. 6. Steady-state streamlines for viscoelastic flow calculated on a 9X 9 grid: (a) D =l, R=O.2~10-~; (b) D=2, R=0.4~10-~; (c) D=3, R=0.6X10-4; and(d) D=4, R= 0.8 x 10-4.
shown in Figs. 6(a-d) and 7(a-d), respectively. The computations on both grids predict very similar flow behavior and the results can be considered convergent. At D = 1, a small secondary vortex which grows with increasing D has formed in the upper left-hand corner of the flow. At D = 2, this secondary vortex has grown much larger and fills the entire bottom left-hand portion of the cavity. The calculation done on the 17 x 17 grid for D = 2 also predicts a third vortex along the right-hand side of the cavity. At D = 3, the secondary vortex occupies over 75% of the bottom portion of the cavity
216 a)
Cd)
Fig. 7. Steady-state streamlines for viscoelastic flow calculated on a 17 X 17 grid: (a) D = 1, R = 0.2~10-~; (b) D = 2, R = O.4x1O-3; (c) D = 3, R = 0.6~10-~; and (d) D = 4, R = 0.8x 10-4.
and the primary vortex which drives the flow is confined to a thin layer next to the upper plate. The calculation for D = 4 is very similar to this except that at D = 4 the intensity of the vortices has increased. The pressure distributions predicted by both sets of calculations are shown in Figs. 8(a-d) and 9(a-d). The increased smoothness of the calculations performed in the 17 x 17 grid indicate the convergent behavior of the scheme. The pressure contours for the viscoelastic calculations show significant deviation from the Newtonian case. In the Newtonian calculations, the pressure distribution is divided into a low pressure side on the left and a
217 La)
(b)
(dj
Fig. 8. Steady-state pressure field for viscoelastic flow calculated on a 9X9 grid: (a) D =l, R = 0.2x 10-3; (b) D = 2, R = 0.4x 10-3; (c) D = 3, R = 0.6x 10-4; and (d) D = 4, R = 0.8 x 10-4.
high pressure side on the right. In the viscoelastic calculations, the contours become somewhat circular and the high pressure region is confined to the inner portion of the flow near the center of the cavity while the area along all the walls of the flow are low pressure regions. Information on the extra stress in the flow is presented in terms of quantities that might be readily measurable in experiments with viscoelastic fluids in a cavity flow. The contours for the first principal stress difference (defined as A@ = [(S - T)’ + 4Q2]‘12) predicted by the calculations on both grids are shown in Figs. lO(a-d) and ll(a-d). The calculations are in
218 (a)
(b)
(d)
Fig. 9. Steady-state pressure field for viscoelastic flow calculated on a 17 X 17 grid: (a) D = 1, R = 0.2x 10-3; (b) D = 2, R = 0.4x 10-3; (c) D = 3, R = 0.6 x 10-4; and (d) D = 4, R = 0.8 x 10-4.
good agreement as far as the contour patterns go, but differ somewhat in terms of the contour values. This is due to the fact that at the singular upper left-hand comer of the flow the velocity gradient is twice as high on the more refined mesh. This higher gradient generates a greater stress level at that point and this effect persists throughout the flow. The contours for the viscoelastic calculations show evidence of sharp gradients in stress along all the walls of the flow in contrast to the Newtonian calculations in which the large stress gradients are confined to the region along the upper wall of the flow.
219 (a)
(b)
Cd)
Fig. 10. Steady-state first principal stress difference distribution lated on a 9x9 grid: (a) D=l, R=0.2X10m3; (b) D=2, R=0.6x 10-4; and(d) D=4, R=0.8x 10-4.
for viscoelastic flow cakuR=0.4X10m3; (c) D=3,
The integral drag and normal forces along the upper wall of the flow as a function of time for the calculations on the 17 X 17 grid are shown in Figs. 12(a)-(d). These forces are given by the equations Fd =
J’
Q
dx,
(79)
-1
F,=
l(Po-P+T)
/-1
dx.
(b)
Id)
Fig. 11. Steady-state first principal stress difference distribution lated on a 17x17 grid: (a) D=l, R=0.2X10-3; (b) D=2, R = 0.6 x 10-4; and (d) D = 4, R = 0.8 X 10-4. (a)
(b)
(cl
for viscoelastic flow cakuR=0.4X10-3; (c) D= 3,
Cd)
1
0.5
0
0
20 0
DRAG
40 TIME +
60
80
NORMAL
Fig. 12. Integral drag force and normal force as a function of time for viscoelastic flow calculated on a 17x17 grid: (a) D ‘1, R = 0.2~10~~; (b) D = 2, R = O.4X1O-3; (c) D = 3, R=0.6x10-4; and(d) D=4, R=0.8x10-4.
221 5. Discussion A purely hyperbolic model and numerical method for unsteady, viscoelastic flow has been developed and applied to the flow in a driven cavity. The model allows a slight degree of compressibility in order to obtain a completely hyperbolic system of equations. The equations are solved using a finite difference method that was developed for the treatment of hyperbolic systems. The model has been formulated in this way specifically to try to work around the fact that many of the traditional systems of equations for viscoelastic flow (i.e., those presuming incompressibility and no redardation time) are mixed systems, for which no general numerical algorithms are available. The computations for a fluid with parameters corresponding to a creeping, Newtonian flow show good agreement with calculations of the same flow that are present in the literature. Based on this, it seems safe to presume that allowing compressibility in the flow is not a major impediment to obtaining solutions that are comparable with completely incompressible calculations. The viscoelastic flow computations illustrate two positive points about the numerical scheme developed here. First, the calculations on the more refined mesh agree with those obtained on the coarser mesh and are better in quality. This is a strong indication that the scheme is convergent with mesh refinement which is especially encouraging in light of the fact that the velocity field for this flow contains two singularities. Secondly, the scheme is stable. Flows exhibiting significant elastic effects at Deborah numbers an order of magnitude higher than those obtained in previous cavity flow calculations [20] have been obtained in this study. While no claim is made for being able to integrate to unlimited Deborah numbers, we note that the calculations presented here have been limited to D = 4 due to the heavy computational expenses associated with this method-not numerical instability. These results indicate that the hyperbolic differencing scheme employed to solve our model is successful because the mathematical problem is correctly treated numerically. Based on the recent success of other numerical schemes [3,4], correct numerical treatment of the mathematical problem seems to be the key to obtaining stable and convergent viscoelastic flow calculations. An outstanding feature of the calculations is the prediction of multiple vortices in the cavity flow under creeping flow conditions at Deborah numbers greater than unity. While this has not yet been witnessed experimentally, it seems to be similar in nature to the phenomenon observed by Parlato [21] in experiments with viscoelastic fluids in a four-roll mill. At high Deborah numbers the fluid could no longer be dragged through the mill, but instead, the flow reversed itself and emptied the apparatus of material.
222 Acknowledgements The authors are grateful to General Electric, Plastics Division, and Exxon for support of this research, and the National Science Foundation and the Pittsburgh Supercomputing Center for computational resources. Notation a a A A+ Ab i B B B+ Bc ,x Fl D/Dt ei E f h k M ,Y ; P PO
Q
R S t
Vector used in compatibility equations containing x-direction derivatives terms Expression which is a function of S, D and R Coefficient matrix of x-direction derivatives Split coefficient matrix of x-direction derivatives Split coefficient matrix of x-direction derivatives Vector containing the right-hand side of the compatibility equations Modified version of b with certain components set to zero Expression which is a function of T, D and R Coefficient matrix of y-direction derivatives Split coefficient matrix of y-direction derivatives Split coefficient matrix ,of y-direction derivatives Vector containing the right-hand side of the compatibility equations Modified version of c with certain components set to zero, Deborah number Material derivative Eigenvector Elasticity number Vector of forcing terms in viscoelastic model Denotes Ax Denotes Ay Matrix used in compatibility expressions; modified version of x-’ Unit normal at wall Matrix used in compatibility expressions; modified version of Y-l Pressure Reference pressure Dimensionless shear stress Reynolds number Dimensionless first normal stress Time
223 At T u V U
Ix
X X-l AY Y Y-l
Change in time Dimensionless second normal stress Dimensionless velocity, x-direction Dimensionless velocity, y-direction Velocity vector Vector of dependent variables Grid-spacing in x-direction Similarity transform matrix for A Similarity transform matrix for A; matrix of left eigenvectors Grid-spacing in y-direction Similarity transform matrix for B Similarity transform matrix for B
Greek
Giesekus parameter Compressibility parameter Rate of strain tensor Viscosity parameter in Giesekus equation Shear viscosity function Relaxation time i* characteristic direction in the x-t plane 1‘th characteristic direction in the y-t plane Estimate of largest eigenvalue in x-t and y-t planes Matrix containing eigenvalues of the matrix A Matrix containing positive eigenvalues of the matrix A Matrix containing negative eigenvalues of the matrix A Matrix containing eigenvalues of the matrix B Matrix containing positive eigenvalues of the matrix B Matrix containing negative eigenvalues of the matrix B Density Reference density Unit vector Component of v Component of v Extra stress tensor Stream function Upper-convected derivative References 1 M.J. Crochet, A.R. Davis and K. Walters, Numerical Simulation of Non-Newtonian Flow, Elsevier, Amsterdam, 1984.
224 2 R. Keunings, in: C.L. Tucker III (Ed.), Fundamentals of Computer Modeling for Polymer Processing. Carl Hanser Verlag, Miichen, 1988, Ch. 9. 3 J.M. Marchal and M.J. Crochet, J. Non-Newtonian Fluid Mech., 26 (1986) 77-114. 4 R.C. King, M.R. Apelian, R.C. Armstrong and R.A. Brown, Proceedings of 5th Workshop on Numerical Methods in Non-Newtonian Fluid Mechanics, Lake Arrowhead, CA, 1987. J. Non-Newtonian Fluid Mechanics, 29 (1988) 147-216. 5 D.D. Joseph, M.M. Renardy and J.C. Saut, Arch. Rat. Mech. Anal., 87(3) (1985) 213. 6 D.D. Joseph and J.C. Saut, J. Non-Newtonian Fluid Mech., 20 (1986) 117-141. 7 F. Dupret and J.M. Marchal, J. Non-Newtonian Fluid Mech., 20 (1986) 143-171. 8 J.H. Song and J.Y. Yoo, J. Non-Newtonian Fluid Mech., 24 (1987) 221-243. 9 R.A. Brown, R.C. Armstrong, A.N. Beris and P.-W. Yeh, Comput. Meth. Appl. Mech. Eng., 58 (1986) 201-226. 10 M.M. Renardy, Z. Angew. Math. Mech., 65 (1985) 449-451. 11 F. Dupret, J.M. Marchal and M.J. Crochet, J. Non-Newtonian Fluid Mech., 18 (1985) 173-185. 12 J.V. LeBlanc, Ph.D. Dissertation, Univ. of Massachusetts, Amherst, MA, 1985. 13 H. Giesekus, J. Non-Newtonian Fluid Mech., 11 (1982) 69-109. 14 D.A. Anderson, S.R. Chakravarthy and M.D. Salas, Computational Fluid Mechanics and Heat Transfer. McGraw-Hill, New York, NY, 1983. 15 E.C. Zachmanoglou and D.W. Thoe, Introduction to Partial Differential Equations with Applications, Williams and Wilkins, 1976. 16 S.R. Chakravarthy, Ph.D. Thesis, Iowa State University, IA, 1979. S.R. Chakravarthy, D.A. Anderson and M.D. Salas, AIAA Paper No. 80-0268, Pasadena, CA, 1980. 17 P.J. Roache, Computational Fluid Dynamics. Hermosa, Revised ed., Albuquerque, NM, 1982, pp. 184-185. 18 J.Y. Yoo and D.D. Joseph, J. Non-Newtonian Fluid Mech., 19 (1985) 15-41. 19 F. Pan and A. Acrivos, J. Fluid Mech., 28 (1967) 634-655. 20 M.A. Mendelson, P.-W. Yeh, R.C. Armstrong and R.A. Brown, J. Non-Newtonian Fluid Mech., 10 (1982) 31-54. 21 P. Parlato, M. Chem. Eng. Thesis, University of Delaware, DE, 1969. 22 M. Renardy, personal communication, 1988. 23 M. Renardy, W.J. Hrusa and J.A. Nobel, Mathematical Problems in Viscoelasticity. Longman Scientific & Technical, New York, NY, 1987, pp. 33-34.
A PURELY HYPERBOLIC VISCOELASTIC FLOW
197
MODEL FOR UNSTEADY
F.R. PHELAN, Jr., M.F. MALONE and H.H. WINTER Department of Chemical Engineering, Goessmann Luboratoly, University of Massachusetts Amherst, MA 01003 (U.S.A.) (Received June 27,1988;
in revised form December 22,1988)
Summary A purely hyperbolic model for unsteady, two-dimensional viscoelastic flow is developed. The model uses a pressure-density relationship in the continuity equation to obtain a dynamic equation for the pressure, and the Giesekus constitutive equation with zero retardation time to describe the extra stress. The small degree of compressibility allowed by the model is negligible for practical purposes but important mathematically because it changes the unsteady governing equations from parabolic-hyperbolic to a purely hyperbolic system. This allows the application of a numerical method designed specifically for solving hyperbolic problems, and thus, to better treat the hyperbolic nature of the flow which is dominant at high Deborah numbers. The method is used to predict the velocity and stress fields in a square cavity. The results for the creeping, Newtonian limit compare well with previously published calculations. Viscoelastic flow computations are presented for Deborah numbers up to 4. Comparison of the results obtained on several grids show that the calculations increase in smoothness with mesh refinement, and that the scheme is numerically stable. The calculations predict flows exhibiting significant elastic effects. A transition from single vortex flow to flow containing multiple vortices, occurs under creeping flow conditions for Deborah numbers greater than 1.
1. Introduction The numerical simulation of viscoelastic flow has been an active area of research in computational fluid mechanics for over a decade. The persistent 0377-0257/89/$03.50
0 1989 Elsevier Science Publishers B.V.
198 failure of many numerical schemes when elastic effects are significant has been an outstanding problem and is well documented in the literature [1,2]. Various investigators have laid blame for this problem on the use of fundamentally unsound constitutive equations, failure to properly treat stress singularities, improper specification of stress boundary conditions, excessive approximation error due to an inadequate ability to resolve the large gradients which may occur in these flows, or bifurcation from two-dimensional flow to three-dimensional or time-dependent flow. While all of these are probably problematic to some degree, recent studies indicate that the root of the numerical difficulties is associated with the use of numerical schemes that are inappropriate for the mathematical problem. Such schemes give rise to fictitious limit points that are not an actual property of the mathematical system but are due solely to the numerical approximation procedure [3,4]. Joseph et al. [5] were the first to suggest that the complicated mathematical nature of the equations that are commonly used to simulate viscoelastic flows may be related to the numerical difficulties. They classified the equations governing the steady, two-dimensional flow of the upper convected Maxwell fhtid, as well as those governing a class of Oldroyd models of which the Maxwell fluid is a special case. They showed that a change in type from elliptic to hyperbolic occurs in the equation for the vorticity when certain conditions depending on the local velocity and stresses fields and on the material parameters of the fluid are satisfied. Joseph and Saut [6] have classified the equations governing the unsteady flow of viscoelastic fluids with constitutive equations of the form
(1) where S/St represents the upper-convected derivative. They have shown that a loss of evolution (Hadamard instability) occurs when certain critical values of the normal stress are reached. Dupret and Marchal [7] have found that a Johnson-Segahnan Fluid exhibits the same behavior. In some problems, a change of type in a model is reflected in the actual flow and is a stable physical phenomenon, e.g., the onset of shocks in transonic gas flow. However, the analysis of Joseph and Saut [6] shows that a change of type in viscoelastic flow problems is not always stable. In problems without inertia, the point at which a change of type occurs in the steady problem is exactly the same as that at which the corresponding unsteady problem becomes non-evolutionary. Thus, a steady-state numerical simulation in which a change of type occurs is sometimes susceptible to Hadamard instabilities which grow in severity on more refined computational meshes. Several investigators [8,9] have linked the occurence of a
199 change of type in their numerical schemes to a subsequent loss of convergence. While the use of constitutive equations which change type and become non-evolutionary is certainly problematic, it does not entirely explain the problems that have limited past computations. Renardy [lo] has shown that a change of type cannot occur for the Upper Convected Maxwell Fluid at zero Reynolds number. This contradicts the results of various numerical studies which have shown that a change of type in the flow an Upper Convected Maxwell Fluid at zero Reynolds number led to a loss of convergence [ll]. This contradiction suggests that the loss of convergence in traditional numerical schemes is due to errors which lead to an artificial change of type (and thus, a non-evolutionary problem) and that numerical simulations can be improved by developing schemes which properly treat the mathematical type of the system. Such schemes could not be expected to obtain solutions when the mathematical problem really loses evolution, for in such cases the problem is ill-posed. However, they would provide a means for preventing numerically induced instabilities in problems that are actually well-posed. The studies of Joseph and co-workers have emphasized the importance of recognizing the type of the system in models for viscoelastic flow. One of the more important results that has emerged is that a wide range of viscoelastic models are mixed systems, being hyperbolic-elliptic in the steady-state and hyperbolic-parabolic in the unsteady-state. This situation is difficult to treat numerically because most algorithms are designed for systems of equations of a single type and methods for mixed systems are not well developed. However, numerical methods for purely hyperbolic systems are available, primarily from studies of compressible flow. LeBlanc [12] formulated a one-dimensional model problem for viscoelastic flows by ignoring the conservation of mass constraint. Numerical treatments of this simple problem displayed many of the essential features of more realistic models and the method of characteristics provided the only satisfactory solutions. This simplified approach provided the starting point for this study. In what follows, a purely hyperbolic model for unsteady two-dimensional flows is developed. Purely hyperbolic behavior is obtained by allowing a small degree of compressibility in the flow. With the incompressible constraint relaxed, the continuity of mass condition provides an evolution equation for the pressure. This equation together with the momentum balance and the Giesekus equation with zero retardation time as the constitutive model (or any expression that can be written in the form of eqn. (1)) constitutes a completely hyperbolic system in time and space provided certain critical normal stress limits which depend upon the Deborah number
200 are not violated. Outside of these limits the model becomes non-evolutionary and the problem is not well-posed. Although this does reflect deficiencies in the constitutive model, these are not particularly relevant to the development of a numerical scheme. Numerical solutions are developed using an algorithm designed for solving systems of hyperbolic equations that arise in gas dynamics. As an example, the flow in a driven cavity is examined. 2. Model formulation 2. I. Basic equations The model employed her consists of the Cauchy momentum balance, the Giesekus [13] constitutive equation and the equation of continuity of mass. In the absence of body forces and for a zero retardation time, these are written as Do pa--V.7+(1/3)V(trr)+VP=O, h
DT i~-VVT*T-T*vv
(2)
1-qq=
-T*(I+(ah/~)T),
where P is formulated here as the mean normal stress which accounts for the extra term involving the gradient of the trace of the stress in the momentum balance (tr 7). A model based on the assumption of incompressible fluid behavior inevitably results in a system of equations of mixed type. In order to obtain a purely hyperbolic system, nearly incompressible fluids are considered as a perturbation of the incompressible state. A first-order Taylor expansion for the density (p) at constant temperature together with the definition of the isothermal compressibility gives P=Po[I
+ (P-P,)/819
where /3 is the inverse of the isothermal compressibility. into eqn. (4) yields the expression
(5) Inserting eqn. (5)
+g+[p+(P-P,)I(v*v)=O. One further simplification is made in this last expression. Because a pertubation of the incompressible state is desired, only cases for which ((P PO)//?1 -=K1 are of interest. Thus, eqn. (6) is simplified to DP ~+/3(v~v)=o.
201 Equations (2), (3) and (7) make up the model to be studied here. These can be written in dimensionless form as E-v*7+(1/3)v(trT)+VP=O, RDU D
(8) -+=
-7=(1+&T),
(9)
(10) where R = (pVL/p) is the Reynolds number, D = (XV/L) is the Deborah number and u, T, P and /3 are the dimensionless velocity, extra stress, pressure and compressibility. For a two-dimensional flow where ?J= [u, VI,
(11)
and
S
‘=Q
Q
I I
(12)
T’
the components nates as
of eqns. (8), (9) and (10) are written in Cartesian coordi-
tLJ+ UU, + Vu, + R-‘(P,
- (2/3)S,
c: + Uv, + WY + R-I( Py + (l/3)$ St+ Us,+
vs,-Au,-2Qu,=
uq+
Pt+ up,+
(13a)
- Q, - (2/3)T,)
= 0,
03b)
Q2)), -(Q/D+(I~Q(S+T)),
V5-2Q<-4= vPy+p(ux+
= 0,
-(S/D++‘+
Q,+UQ,+VQ,-(&2)U,-(&2)V,=
T,+
- Q,, + (1/3)T,)
-(T/D+cx(T~+Q~)), v,)=o,
034
(13d) (13e) W)
where a = 2( S + l/D), h==(T+l/D), and the subscripts t, x and y denote partial-differentiation those variables.
(14) (15) with respect to
2.2. Classification The system of equations (13), can be written in the form wt+A*w,+B*w,=f,
(16)
202 where w is a vector of length six with components, w = [U, V, S, Q, T, PIT, the matrices A and B contain the coefficients of W, and wv, and f is a vector of forcing terms. Equation (16) represents a quasi-linear system of first-order partial differential equations and in order to treat such a system numerically, we must first classify it mathematically [14]. In general, the system of equations (16) is evolutionary in (x, y, t) if for every unit vector v = [ vl, v2] the eigenvalue problem (-XI++4
+Y,B)e=O
07)
has real eigenvalues only, and hyperbolic in (x, y, t) if in addition, there are six linearly independent eigenvectors, e, [6]. However, we present here the classification information that is relevant with regard to the finite difference algorithm to be employed for solving the model equations, and that is the conditions for an evolutionary and hyperbolic problem for the cases v = [l, 0] and v = [0, 11. The information obtained for these particular cases determines separately the conditions for an evolutionary and hyperbolic in the (x-t) and (r-t) planes, respectively [14,15]. The system is evolutionary in (x, t) if the eigenvalues of A are all real. In addition, the system is also hyperbolic if A can be written in the form A =X*A,*P,
(1%
where the diagonal matrix A, contains the eigenvalues of A and X-’ is the matrix of left eigenvectors. The same can be said of the behavior in (y, t) with respect to the eigenvalues and eigenvectors of B. The eigenvalues of the matrix A are L,,
= u,
@a>
L,,
= u,
Wb)
A,,, = u+ /272x,
(1W
A,,,=
U-
Wd)
A,,,=
tiJ+/@GjpY,
(194
A,,,=
u-J@zjq?x,
ww
\ia/zR,
and the eigenvalues of the matrix B are
v, A,,, = v, A,,, =
(204 (2Ob)
203
A,,, =
v-t \lqiiF,
A,,, =
v-
@x,
(2Od)
A,,, =
v-
/@ETjjp.
@Of)
(W
The compressibility parameter, 8, takes on positive values only, and thus, an inspection of eqns. (19) and (20) shows that A and B have real eigenvalues if * S>
-I/D
(21)
-I/D.
(22)
and T>
When the inequalities (21) and (22) are satisfied, the coefficient matrices A and B each have six linearly independent eigenvectors and they can be decomposed by the similarity transforms A =X4,*X-f, B=
Y*A;
(23) Y-r.
(24)
Under these conditions, the equations of the model are both evolutionary and hyperbolic. 2.3. Parameter limits It is not possible to establish general apriori parameter limits on (Yand D for which one of the conditions (21) or (22) will be violated. However, it is useful to examine two simple flows for which analytical solutions can be obtained to determine under what conditions a loss of evolution may occur. First, a pressure-driven shear flow in a slit with a gap width of 2 is considered. The analytical solutions for the shear and normal stresses for a Giesekus fluid with zero retardation time and non-zero (Y are given by Q=
-Y,
(25)
S = (1,‘2aD)(
- 1+ I/g(Y) + &xD’YU,),
(26)
T= (1/‘2aD)(
- I+ /g(y)),
(27)
* Note added in proof: M. Renardy [22] has pointed out that conditions (21) and (22) are necessary but not sufficient to guarantee real eigenvalues in (17). These could be replaced by more general conditions involving the eigenvalues of the extra stress as discussed in [23].
204 where u,=
Y[l+(2a
g(Y)
= 1 - 4cy2D2Y2.
- l)HV]/[(2~
- 1) + \lg(y)12,
(28) (2%
The shear stress has been scaled to a value of unity at the walls of the flow and the pressure-gradient has been absorbed in the scaling. Since the function g(Y) appears as a square-root term in several places, the first condition that must be satisfied is g(1) 2 0. This requirement yields the condition aD I l/2.
(30)
The above condition tells us nothing about hyperbolicity; it is simply a requirement that needs to be satisfied in order that the dimensionless second normal stress be real-valued. One of the conditions for hyperbolicity can be violated if one of the normal stress components, S or T, becomes too large and negative. In this flow, S is always positive and T is always negative. Thus, we require that (Y and D satisfy (1/2cwD)( -1 + m)
> -l/D.
(31)
This yields the condition a( D2 + 1) I 1.
(32)
A physical interpretation of this condition can be obtained by solving eqn. (13d) for the viscosity function in steady shear flow. By doing this, it is found that q(y) = (1+ DT)/(l
+ aD(S
+ T)).
(33)
A violation of the inequality (32) corresponds to a negative viscosity which is a violation of thermodynamics. In other flow situations, condition (32) might not be conservative enough since it has been derived under the stipulation that 1Q I_ = 1. Consider the more general expression for T in a pressure-driven shear flow
(34) If the right hand side of eqn. (34) is substituted into the inequality (22), it is found that the inequality will hold provided that (Y and D satisy a(4D2Q2
+ 1) < 1.
(35)
This more conservative condition becomes necessary in the simulation of flow in a geometry where the outlet is smaller than the inlet, for example, a contraction flow. In such a problem, if the upstream shear stress is scaled
205 such that ] Q ] max= 1, then the downstream value will be much larger than one. The values of (Y and D must then be chosen such that (35) is satisfied where the value of Q in this expression is taken as ] Q ]zet. Finally, planar extensional flow is considered. The velocity field for this flow is given by u= [ix,
-iY],
(36)
and the normal stress components are given by s = (l/2&)(
- (1- 2iD)
+ {(l - 2iD)2 + 8icrD
T = (l/2&)(
- (1+ 2iD)
+ {(l + 2iD)’
),
- 8&D ).
(37) (38)
In this flow, T is negative. In order that the condition (22) not be violated, we require that (l/2&)(
- (1+ 2iD)
+ J(1+
2iD)2 - l&D)
2 -l/D,
(39
which yields the condition (Y
(40)
This condition is not restrictive at all since the Giesekus equation is only valid up to (Y= 1. Thus, somewhat surprisingly, it is independent of D. Based on these results, it is expected there may be parameter limits that must be observed in flow simulations with this model. However, these parameter limits reflect particular deficiencies in the constitutive equation that has been chosen for this study - not the specific procedure that has been developed here. 3. Numerical algorithms 3.1. The Split Coefficient Matrix (SCM) method The Split Coefficient Matrix (SCM) method developed by Chakravarthy et al. [16], is a finite difference method which uses information on signal propagation provided by the theory of characteristics to solve systems of first-order hyperbolic equations. This method is similar in philosophy to the method of characteristics, but much easier to apply to two-dimensional problems. The major advantage of the SCM method over the method of characteristics is that a conventional, rectangular finite difference grid can be employed for the calculations. In the latter method, the characteristic mesh in not known a priori and must be calculated as part of the solution.
206 Detailed discussions of the SCM method are provided by Chakravarthy et al. [16] and Anderson et al. [14]. The application of the SCM method to the system represented by eqn. (16) is relatively straightforward. First, the coefficient matrices A and B are each split into two parts. To accomplish this splitting, the similarity transforms (23) and (24) are used to rewrite (16) as W,+x*A;x-‘*W,+
Y*A,= Y%$=f.
(41)
The diagonal matrices A, and A, are then partitioned negative parts according to the formulas
into positive and
A,=Ax’+A,,
(42)
Ay=A;
(43)
+A,,
where the matrices with the ‘ + ’ superscripts contain the positive eigenvalues of the original matrices and those with the ‘ - ’ superscripts contain the negative eigenvalues. Using these last two equations, we define the splitting of the the matrices A and B as A=A++A-,
w
B=B++B-,
(45)
where the split coefficient matrices are given by A+=X.A,+.X-l,
(46)
A-=X.A,.X-l,
(47)
B+=Y*A;Y-I,
(4%
B-=Y*AJ*Y-?
(49)
Using eqns. (44-49) we can then rewrite (41) as w,+A+
lw,+A-•w,+B+*w,+B-•w,=f.
(50)
Equation (50) is still entirely equivalent to the original system of equations (16). In the second part of the SCM method procedure, finite differences are used to approximate the spatial derivatives in eqn. (50). To correctly treat the hyperbolicity of the system, backward differences are used to approximate those derivatives which are multiplied by A + and B+ and forward differences are used to approximate those derivatives which are multiplied by A - and B-. The SCM methods procedure can then be represented symbolically as w,=f-(A+*w,+
+A-•w,-
+B+vv;
where the ‘ + ’ and ‘ - ’ superscripts
+B-*I+$-),
on W, and wy denote backward
(51)
and
207 forward difference approximations, respectively. Equation (51) represents a set of six first-order, ordinary differential equations to be integrated at every internal point on the finite difference grid. 3.2. Boundary conditions At the boundaries of the computational grid, some of the characteristics point outward from the domain, some are perpendicular to the boundary and some point inward. This is illustrated for the characteristics in the x - t and y - t planes along the boundaries of the geometry shown in Fig. 1. The information propagated along the characteristics that point inward cannot be used in the integration of the equations at the boundary. Thus, the number of equations along the boundary must be reduced by the number of inward pointing characteristics and replaced by the same number of boundary conditions. To illustrate this, the procedure for the approximation of eqn. (41) along the boundary whose characteristics are shown in part a) of Fig. 1 is outlined. The normal SCM method splitting procedure is used to handle the x-direction derivative terms. Backward differences are used to approximate the y-direction derivative terms. Applying this procedure to these respective terms and rearranging, eqn. (41) is rewritten as wt+ Y*h,*Y-l*~$+
=f-a,
(52)
where a=A+
l
w,’ + A -
l
w, .
(53)
Equation (52) is then multiplied by Y-r to yield Y1
l
W,= c,
(54)
where c=
-A,
l
rl*wj+
+ r’*(
f -a).
(55)
Equations (54) are the “compatibility equations” for the y-t plane. Note that each term of the vector c is a function of only one of the eigenvalues of the matrix B. To correctly reduce the number of equations and introduce the boundary conditions, the compatibility equations containing the eigenvalues which point into the computational domain are discarded. For the case of no-slip (and in general, whenever the Reynolds number is small and n u = 0), the eigenvalues X,,, and X,,, are negative. Since that is the case being considered here, the fourth and sixth compatibility are replaced with the no-slip boundary conditions for the two components of the velocity vector. In order to introduce this into the numerical scheme in a straightforl
208 ward manner, the compatibility replaced with the equations
equations that have been thrown out are
Q=O
(W
and q=o.
(57)
The system of equations (54) are thus modified to the set New*=&
(58)
where the vector c^ is derived from c by setting the fourth and sixth components of c to zero and the matrix N is given by,
NE
0
0
Y;;' Y;-,'Y;;' Y;-,’
0
0
Y*
Y$
0
0
Y$
YG1 Ygl
YG1
r;;r
0 Y&l
0 Y;;’
0 I&’
0
0
0
010 0 0 100
Y,3'
YG1
(59)
The matrix N is then inverted to give t:he the following system of equations for the rate-of-change of the dependent variables at the boundary w,=N-‘*c^.
(60)
Following this procedure automatically yields the desired conditions U, = 0 and v = 0 without any further changes in the numerical scheme. While eqn. (60) can be used to find the pressure changes along the boundary, in the course of testing the scheme it was found that a more accurate alternative was to instead, solve for the pressure by enforcing the normal component of the momentum balance. Along a no-slip boundary this is written as ii-vP=ri~v*r.
(61)
For the boundary shown in Fig. 1, eqn. (61) reduces to P,,= -+S,,+Q,+fr,.
(62)
By applying backward differences to the derivatives in the y-direction and central differences to the derivatives in the x-direction we can derive the following expression for the value of P at time step n + 1 at each point (i, j) along this boundary Peal = P,~~‘l-
:( S~+l-
+ ;( T;;ti+l - q;pr),
S,?~~l) + (k/2h)( Q~~~j - Q~~~j) (63)
209
d’
X
Fig. 1. General structure of the characteristic y-f and b) x-t planes.
where
h = Ax and
directions along no-slip boundaries
in
the a)
k = Ay. The treatment of the other boundaries is com-
pletely analogous. There is another situation that warrants mentioning and that is the treatment of comer points where the flow boundaries meet. At these points there are inward pointing characteristics in both the x and y directions. In this case, the boundary treatment shown above cannot be used because the problem requires that the information propagated along both these sets of characteristics be removed from the equations. To accomplish this a timesplitting scheme is employed in conjunction with the normal SCM method techniques used to split the spatial terms. This procedure is illustrated for the approximation of eqn. (41) at the comer labeled c) in Fig. 1. The time-splitting technique consists of partitioning eqn. (41) into two parts as follows ;wt’~‘+x.A,.x-i.w,-
= if
(64)
=
(65)
and $l+@‘+ Y*h;
Y-‘*w,
if.
Each of these split parts are then solved using the same type of boundary procedure illustrated for non-corner boundary points. First, the split compatibility equations are written as +x-f,,,,(i)
= b9
3 pp+p
= c,
(66) (67)
where
b=+X-lf-&.X-‘.w,-
(68)
and
c=:r’f-A;r’v$-.
(69)
210 Then, the appropriate compatibility equations are discarded and replaced with the conditions U, and v = 0. For the case in question, the third and fifth compatibility equations (i.e., the third and fifth components of the vectors b and c) are discarded in both the x and y directions. This procedure yields the equations n/f. &)
= 5
(70)
jjl.
= 2,’
(71)
,;,
where 6 and 2 contain the components of the vectors b and c except for the third and fifth components which are zero and M and N are given by 0
0
x,-,l
x,-,*
x,-,’
x,-,’
0
0
x,-,l
xi1
x,-,l
x,-,l
0
0
0 x,-,’
0 xi1
0 x,-,l
0 x$
M=1O 010
N=
0
0
0
0
0
x,-,’
x,-,’
x,-,’
x,-,l
0
0
Y;;’
YG’
YE1
Y,-,l
0
0
Y2i1
Y,-,’
YG1
Y&l
0
0
0
Yi’
Y&l
Y&l
100 0
0
YG1
’
(72)
(73)
The matrices M and N are then inverted and eqns. (70) and (71) are combined to yield the following system of equations for the rate-of-change of the dependent variables at the corner ,,,,=M-l.&+N-l.c^.
(74) To solve for the pressure at the comer points the momentum balance was used. Since there is no unique normal vector at a comer the two equations used along each of the respective sides were averaged and an equation for the new value of P was derived by applying appropriate finite-difference expressions. 4. Simulation of cavity flow 4.1. Description of geometry and procedure The algorithm described above was used to calculate the velocity, stress and pressure fields for the flow in a square cavity with dimensions: - 1 I x
t l.._m_!
211
Y
x
Fig. 2. Driven cavity flow geometry.
I 1, - 1 I y I 1. This geometry is depicted in Fig. 2. The x-component of the velocity along the upper wall of the cavity was assumed to be a constant of value 1 except at the comer points where it was zero. The velocity was zero on all the other walls and comer points of the flow as well. The initial conditions used for the calculations were II = 0, T = 0 and P = 1 over the entire cavity. Roache [17] discussed the problem of setting a pressure datum in a cavity flow for dynamic calculations. Because the specification of the velocity components completely satisfies the boundary condition requirements for this problem, there is no physically acceptable point at which to set a pressure datum. The procedure used here was to pick a point in the flow at random, and subtract the value of P, at that point from all the other values of Pt in the flow. This has the effect of keeping the value of the pressure constant at this point but still allowing the correct relative changes between points to occur. Roache also suggested that for this problem the final values of pressure be corrected by the addition or subtraction of a numerical constant such that the integrals of the initial and final pressure distributions over the cavity are equal. This procedure was used here but is not absolutely necessary for fluids that behave incompressibly or nearly incompressibly because only pressure differences are important in such a problem. Several methods were tested for the integration of the discretized equation set including a fourth-order Runge-Kutta, an Adams-Bashford-Moulton predictor-corrector and a simple Euler scheme. Because the time-step involved in the calculations was usually of the order of 10e5, it was found that the additional computational expenses incurred by the use of the higher order methods was not worth the very small gains in accuracy, especially, since all steady-state results are independent of the integrator employed. Thus, the majority of the calculations were performed using Euler’s method. Irrespective of the integrator used, the following time step limitations were observed to provide stability
At < MWmx>@xAY/@x+ AY ))
(75)
212 for r < 1, and At < 1/(2r.
A,,)(
AX Ay/( AZ + Ay))
(76)
for r > 1, where
(77) and r=/i@.
(78)
The maximum value of the divergence in the flow is controlled by the value of the parameter j3. Fairly large values of B are needed to insure a nearly incompressible flow. For this particular problem, values between 100 and 10,000 were satisfactory. The differences that one observes between a calculation at a low value of /3 and one at a high value of p are generally very small at a macroscopic level but not necessarily at a microscopic level. For example, two sets of calculations with constant values of R, D and (Y, but different values of B typically predict the same vortex locations and overall flow pattern, but small differences with regard to the exact shape of the streamlines may be observed. 4.2. Newtonian limit Initial calculations were performed with R = D = (Y= 0.001, and fl= 10,000, on a 17 X 17 grid. This represents a creeping, Newtonian flow limit. These calculations were undertaken to test the accuracy of the scheme by comparing with other Newtonian calculations of the same flow that are available in the literature. It should be noted that the scheme is not suited to handle the absolute limits of R and D = 0. Changing the scaling from a Deborah number to an Elasticity number [18] would allow one to accomodate the limit R = 0; however, under no condition is a zero value of either of these elasticity scaling choices acceptable since this would change the nature of the problem. The streamlines for D = 0.001 are shown in Fig. 3. The streamline values were normalized to a value of 1 along all four boundaries of the flow. The calculations show a vortex located at the point x = 0 and y = 1.5, and a difference between the streamline value at the vortex and the wall of -0.203. This is in good agreement with the results of Pan and Acrivos [19] and other investigators who calculated a vortex center of x = 0.0 and y = 1.52 and a streamline difference of -0.2 for the creeping flow of a Newtonian fluid. The small differences that exist between the two calculations are probably attributable to a number of factors such as the non-zero Reynolds and Deborah number, and finite compressibility used in our
213
L
Fig. 3. Steady-state streamlines for creeping, Newtonian flow.
Fig. 4. Steady-state pressure field for creeping, Newtonian flow.
214
Fig. 5. Steady-state first principal stress difference distribution
for creeping, Newtonian
flow.
calculations, and also possibly to the fact that the grid used by Pan and Acrivos was much more refined than the grid used here. The pressure field and first principal stress difference calculated for this flow are shown in Figs. 4 and 5. 4.3. Viscoehstic j7ow Calculating the start-up flow for a viscoelastic fluid from rest as was done for the Newtonian simulations turned out to be very time-consuming. Therefore, to save on CPU-time it was decided to use instead the steady-state results from the simulation for one Deborah number as input for a simulation at a higher Deborah number, but with the numerical value of the plate speed increased. Thus, the plate speed was used as a parameter to increase the Deborah number. Viscoelastic flow computations were performed with R -C 0.001, a = 0.25, /3= 500 for Deborah numbers from 1 to 4 on 9 X 9 and 17 X 17 grids. Calculations were considered to be steady when the changes in the values of the flow variables were of order 10v6 per 100 time steps and the integral drag and normal forces along the upper wall of the flow had assumed steady values. The steady-state streamlines for D = 1, 2, 3 and 4 for both grids are
a)
(cl
Cd)
Fig. 6. Steady-state streamlines for viscoelastic flow calculated on a 9X 9 grid: (a) D =l, R=O.2~10-~; (b) D=2, R=0.4~10-~; (c) D=3, R=0.6X10-4; and(d) D=4, R= 0.8 x 10-4.
shown in Figs. 6(a-d) and 7(a-d), respectively. The computations on both grids predict very similar flow behavior and the results can be considered convergent. At D = 1, a small secondary vortex which grows with increasing D has formed in the upper left-hand corner of the flow. At D = 2, this secondary vortex has grown much larger and fills the entire bottom left-hand portion of the cavity. The calculation done on the 17 x 17 grid for D = 2 also predicts a third vortex along the right-hand side of the cavity. At D = 3, the secondary vortex occupies over 75% of the bottom portion of the cavity
216 a)
Cd)
Fig. 7. Steady-state streamlines for viscoelastic flow calculated on a 17 X 17 grid: (a) D = 1, R = 0.2~10-~; (b) D = 2, R = O.4x1O-3; (c) D = 3, R = 0.6~10-~; and (d) D = 4, R = 0.8x 10-4.
and the primary vortex which drives the flow is confined to a thin layer next to the upper plate. The calculation for D = 4 is very similar to this except that at D = 4 the intensity of the vortices has increased. The pressure distributions predicted by both sets of calculations are shown in Figs. 8(a-d) and 9(a-d). The increased smoothness of the calculations performed in the 17 x 17 grid indicate the convergent behavior of the scheme. The pressure contours for the viscoelastic calculations show significant deviation from the Newtonian case. In the Newtonian calculations, the pressure distribution is divided into a low pressure side on the left and a
217 La)
(b)
(dj
Fig. 8. Steady-state pressure field for viscoelastic flow calculated on a 9X9 grid: (a) D =l, R = 0.2x 10-3; (b) D = 2, R = 0.4x 10-3; (c) D = 3, R = 0.6x 10-4; and (d) D = 4, R = 0.8 x 10-4.
high pressure side on the right. In the viscoelastic calculations, the contours become somewhat circular and the high pressure region is confined to the inner portion of the flow near the center of the cavity while the area along all the walls of the flow are low pressure regions. Information on the extra stress in the flow is presented in terms of quantities that might be readily measurable in experiments with viscoelastic fluids in a cavity flow. The contours for the first principal stress difference (defined as A@ = [(S - T)’ + 4Q2]‘12) predicted by the calculations on both grids are shown in Figs. lO(a-d) and ll(a-d). The calculations are in
218 (a)
(b)
(d)
Fig. 9. Steady-state pressure field for viscoelastic flow calculated on a 17 X 17 grid: (a) D = 1, R = 0.2x 10-3; (b) D = 2, R = 0.4x 10-3; (c) D = 3, R = 0.6 x 10-4; and (d) D = 4, R = 0.8 x 10-4.
good agreement as far as the contour patterns go, but differ somewhat in terms of the contour values. This is due to the fact that at the singular upper left-hand comer of the flow the velocity gradient is twice as high on the more refined mesh. This higher gradient generates a greater stress level at that point and this effect persists throughout the flow. The contours for the viscoelastic calculations show evidence of sharp gradients in stress along all the walls of the flow in contrast to the Newtonian calculations in which the large stress gradients are confined to the region along the upper wall of the flow.
219 (a)
(b)
Cd)
Fig. 10. Steady-state first principal stress difference distribution lated on a 9x9 grid: (a) D=l, R=0.2X10m3; (b) D=2, R=0.6x 10-4; and(d) D=4, R=0.8x 10-4.
for viscoelastic flow cakuR=0.4X10m3; (c) D=3,
The integral drag and normal forces along the upper wall of the flow as a function of time for the calculations on the 17 X 17 grid are shown in Figs. 12(a)-(d). These forces are given by the equations Fd =
J’
Q
dx,
(79)
-1
F,=
l(Po-P+T)
/-1
dx.
(b)
Id)
Fig. 11. Steady-state first principal stress difference distribution lated on a 17x17 grid: (a) D=l, R=0.2X10-3; (b) D=2, R = 0.6 x 10-4; and (d) D = 4, R = 0.8 X 10-4. (a)
(b)
(cl
for viscoelastic flow cakuR=0.4X10-3; (c) D= 3,
Cd)
1
0.5
0
0
20 0
DRAG
40 TIME +
60
80
NORMAL
Fig. 12. Integral drag force and normal force as a function of time for viscoelastic flow calculated on a 17x17 grid: (a) D ‘1, R = 0.2~10~~; (b) D = 2, R = O.4X1O-3; (c) D = 3, R=0.6x10-4; and(d) D=4, R=0.8x10-4.
221 5. Discussion A purely hyperbolic model and numerical method for unsteady, viscoelastic flow has been developed and applied to the flow in a driven cavity. The model allows a slight degree of compressibility in order to obtain a completely hyperbolic system of equations. The equations are solved using a finite difference method that was developed for the treatment of hyperbolic systems. The model has been formulated in this way specifically to try to work around the fact that many of the traditional systems of equations for viscoelastic flow (i.e., those presuming incompressibility and no redardation time) are mixed systems, for which no general numerical algorithms are available. The computations for a fluid with parameters corresponding to a creeping, Newtonian flow show good agreement with calculations of the same flow that are present in the literature. Based on this, it seems safe to presume that allowing compressibility in the flow is not a major impediment to obtaining solutions that are comparable with completely incompressible calculations. The viscoelastic flow computations illustrate two positive points about the numerical scheme developed here. First, the calculations on the more refined mesh agree with those obtained on the coarser mesh and are better in quality. This is a strong indication that the scheme is convergent with mesh refinement which is especially encouraging in light of the fact that the velocity field for this flow contains two singularities. Secondly, the scheme is stable. Flows exhibiting significant elastic effects at Deborah numbers an order of magnitude higher than those obtained in previous cavity flow calculations [20] have been obtained in this study. While no claim is made for being able to integrate to unlimited Deborah numbers, we note that the calculations presented here have been limited to D = 4 due to the heavy computational expenses associated with this method-not numerical instability. These results indicate that the hyperbolic differencing scheme employed to solve our model is successful because the mathematical problem is correctly treated numerically. Based on the recent success of other numerical schemes [3,4], correct numerical treatment of the mathematical problem seems to be the key to obtaining stable and convergent viscoelastic flow calculations. An outstanding feature of the calculations is the prediction of multiple vortices in the cavity flow under creeping flow conditions at Deborah numbers greater than unity. While this has not yet been witnessed experimentally, it seems to be similar in nature to the phenomenon observed by Parlato [21] in experiments with viscoelastic fluids in a four-roll mill. At high Deborah numbers the fluid could no longer be dragged through the mill, but instead, the flow reversed itself and emptied the apparatus of material.
222 Acknowledgements The authors are grateful to General Electric, Plastics Division, and Exxon for support of this research, and the National Science Foundation and the Pittsburgh Supercomputing Center for computational resources. Notation a a A A+ Ab i B B B+ Bc ,x Fl D/Dt ei E f h k M ,Y ; P PO
Q
R S t
Vector used in compatibility equations containing x-direction derivatives terms Expression which is a function of S, D and R Coefficient matrix of x-direction derivatives Split coefficient matrix of x-direction derivatives Split coefficient matrix of x-direction derivatives Vector containing the right-hand side of the compatibility equations Modified version of b with certain components set to zero Expression which is a function of T, D and R Coefficient matrix of y-direction derivatives Split coefficient matrix of y-direction derivatives Split coefficient matrix ,of y-direction derivatives Vector containing the right-hand side of the compatibility equations Modified version of c with certain components set to zero, Deborah number Material derivative Eigenvector Elasticity number Vector of forcing terms in viscoelastic model Denotes Ax Denotes Ay Matrix used in compatibility expressions; modified version of x-’ Unit normal at wall Matrix used in compatibility expressions; modified version of Y-l Pressure Reference pressure Dimensionless shear stress Reynolds number Dimensionless first normal stress Time
223 At T u V U
Ix
X X-l AY Y Y-l
Change in time Dimensionless second normal stress Dimensionless velocity, x-direction Dimensionless velocity, y-direction Velocity vector Vector of dependent variables Grid-spacing in x-direction Similarity transform matrix for A Similarity transform matrix for A; matrix of left eigenvectors Grid-spacing in y-direction Similarity transform matrix for B Similarity transform matrix for B
Greek
Giesekus parameter Compressibility parameter Rate of strain tensor Viscosity parameter in Giesekus equation Shear viscosity function Relaxation time i* characteristic direction in the x-t plane 1‘th characteristic direction in the y-t plane Estimate of largest eigenvalue in x-t and y-t planes Matrix containing eigenvalues of the matrix A Matrix containing positive eigenvalues of the matrix A Matrix containing negative eigenvalues of the matrix A Matrix containing eigenvalues of the matrix B Matrix containing positive eigenvalues of the matrix B Matrix containing negative eigenvalues of the matrix B Density Reference density Unit vector Component of v Component of v Extra stress tensor Stream function Upper-convected derivative References 1 M.J. Crochet, A.R. Davis and K. Walters, Numerical Simulation of Non-Newtonian Flow, Elsevier, Amsterdam, 1984.
224 2 R. Keunings, in: C.L. Tucker III (Ed.), Fundamentals of Computer Modeling for Polymer Processing. Carl Hanser Verlag, Miichen, 1988, Ch. 9. 3 J.M. Marchal and M.J. Crochet, J. Non-Newtonian Fluid Mech., 26 (1986) 77-114. 4 R.C. King, M.R. Apelian, R.C. Armstrong and R.A. Brown, Proceedings of 5th Workshop on Numerical Methods in Non-Newtonian Fluid Mechanics, Lake Arrowhead, CA, 1987. J. Non-Newtonian Fluid Mechanics, 29 (1988) 147-216. 5 D.D. Joseph, M.M. Renardy and J.C. Saut, Arch. Rat. Mech. Anal., 87(3) (1985) 213. 6 D.D. Joseph and J.C. Saut, J. Non-Newtonian Fluid Mech., 20 (1986) 117-141. 7 F. Dupret and J.M. Marchal, J. Non-Newtonian Fluid Mech., 20 (1986) 143-171. 8 J.H. Song and J.Y. Yoo, J. Non-Newtonian Fluid Mech., 24 (1987) 221-243. 9 R.A. Brown, R.C. Armstrong, A.N. Beris and P.-W. Yeh, Comput. Meth. Appl. Mech. Eng., 58 (1986) 201-226. 10 M.M. Renardy, Z. Angew. Math. Mech., 65 (1985) 449-451. 11 F. Dupret, J.M. Marchal and M.J. Crochet, J. Non-Newtonian Fluid Mech., 18 (1985) 173-185. 12 J.V. LeBlanc, Ph.D. Dissertation, Univ. of Massachusetts, Amherst, MA, 1985. 13 H. Giesekus, J. Non-Newtonian Fluid Mech., 11 (1982) 69-109. 14 D.A. Anderson, S.R. Chakravarthy and M.D. Salas, Computational Fluid Mechanics and Heat Transfer. McGraw-Hill, New York, NY, 1983. 15 E.C. Zachmanoglou and D.W. Thoe, Introduction to Partial Differential Equations with Applications, Williams and Wilkins, 1976. 16 S.R. Chakravarthy, Ph.D. Thesis, Iowa State University, IA, 1979. S.R. Chakravarthy, D.A. Anderson and M.D. Salas, AIAA Paper No. 80-0268, Pasadena, CA, 1980. 17 P.J. Roache, Computational Fluid Dynamics. Hermosa, Revised ed., Albuquerque, NM, 1982, pp. 184-185. 18 J.Y. Yoo and D.D. Joseph, J. Non-Newtonian Fluid Mech., 19 (1985) 15-41. 19 F. Pan and A. Acrivos, J. Fluid Mech., 28 (1967) 634-655. 20 M.A. Mendelson, P.-W. Yeh, R.C. Armstrong and R.A. Brown, J. Non-Newtonian Fluid Mech., 10 (1982) 31-54. 21 P. Parlato, M. Chem. Eng. Thesis, University of Delaware, DE, 1969. 22 M. Renardy, personal communication, 1988. 23 M. Renardy, W.J. Hrusa and J.A. Nobel, Mathematical Problems in Viscoelasticity. Longman Scientific & Technical, New York, NY, 1987, pp. 33-34.