Spurt phenomena of the Johnson-Segalman fluid and related models

303

JournaI of Non-Newtonian Fluid Mechanics, 29 (1988) 303-335 Elsevier Science Publishers B.V., Amsterdam - Printed in The Netherlands

SPURT PHENOMENA OF THE JOHNSON-SEGALMAN AND RELATED MODELS

R.W. KOLKKA’, D.S. MALKUS and R.A. WORTHING i (Received

August

*, M.G. HANSEN

1, 1987; in revised form February

3, G.R. IERLEY

FLUID

’

25, 1988)

summary We examine the behavior of viscoelastic fluid models which exhibit local extrema of the steady shear stress. For the Johnson-Segalman and Giesekus models, a variety of steady singular solutions with jumps in shear rate are constructed and their stability to one dimensional disturbances analyzed. It is found that flow-rate versus imposed stress curves in slit-die flow fit experimental observation of the “spurt” phenomenon with some precision. The flow curves involve linearly stable singular solutions, but some assumptions on the dynamics of the spurt process are required. These assumptions are tested by a semi-implicit finite element solution technique which allows solutions to be efficiently integrated over the very long time-scale involved. The Johnson-Segalman model with added Newtonian viscosity is used in the calculations. It is found that the assumptions required to model spurt are satisfied by the dynamic model. The dynamic model also displays a characteristic “latency time” before the spurt ensues and a characteristic “shape memory” hysteresis in load/unload cycles. These as well as other features of the computed solutions should be observable experimentally. We conclude that constitutive equations with shear stress extrema are not necessarily flawed, that their predicted behavior may appear to be arrested “wall slip”, and that such behaviors may actually have been observed already.

’ Fluids Research Oriented Group and Department of Mathematics, Michigan Technological University, Houghton, MI 49931-1295. ’ Center for the Mathematical Sciences, Department of Engineering Mechanics, and Rheology Research Center, -University of Wisconsin, Madison, WI 53705. 3 F.R.O.G. and Department of Chemical Engineering, M.T.U. 0377-0257/88/$03.50

0 1988 Elsevier

Science Publishers

B.V.

304 1. Introduction Over the past two decades, the question of stability of viscoelastic parallel shear flows has come into prominence. It has often been suggested that analysis of the stability of viscoelastic parallel shear flows can lead to a better understanding of the complex mechanisms of various types of melt fracture, and turbulent drag reduction by polymer additives. In this article we focus on essential constitutive instabilities, as opposed to inertial instabilities. The latter are characterized by relatively high critical Reynolds number, Re, = 2000 + 6200, and small Weissenberg Number, We = 0(10p2). Constant shear viscosity fluids, such as the upper-convected Maxwell (UCM), Oldroyd-B, second-order-fluid (SOF), and of course Newtonian fluid [l], have been found to exhibit such instabilities. Monotonic shear-thinning models, such as the elastic Giesekus model [1,2], and inelastic models considered by Rimmer [3], can also exhibit interial instabilities. Spurt instabilities, which are referred to as “inherent constitutive” instabilities by Petrie and Denn [4], were originally experimentally observed by Tordella [5] and later by Vinogradov et al. [6], and by Blyler and Hart [7]. Such instabilities occur at relatively low Re and high We. A theoretical explanation was first put forth by Huseby [8] and Lin [9]. The ideas of Huseby were realized by McLeish and Ball [lo] in a modified Doi-Edwards model. These authors present a complete and apparently internally consistent explanation of the spurt phenomenon and can match the data of Vinogradov et al. with precision. We present an alternative idea which is more closely related to the dynamic analysis of Hunter and Slemrod [ll], which uses the Johnson-Segalman model. Our contribution is two-fold: First, our model is simple enough for use to consider its linear stability analysis and to compute fully dynamic, one-dimensional simulations of the spurt process. Second, our model can match Vinogradov et al’s data with at least equal precision of that of McLeish and Ball, but the nature of the dynamic process-in particular the shape of the load/unload hysteresis loop-is different in our model. These differences should be experimentally distinguishable. The essential idea is that if a viscoelastic fluid possesses two widely separated relaxation times, then it can exhibit a non-monotonic (double-valued) steady shear response. The presence of shear stress maxima and minima results in jumps of “spurts” as the driving pressure gradient or wall velocity is increased beyond critical values. The Johnson-Segalman model [1,12] with added Newtonian viscosity can exhibit such jump phenomena and is the model we analyze in detail in Section 3. A key difference between our spurt model and that of McLeish and Ball is our very small Newtonian viscosity term. We remark that any model which monotonically shear-thins cannot exhibit the spurt instability we examine.

305 2. Constructing steady solutions and linear stability analysis We analyze the spurt instability of the Johnson-Segalman model for the plane Poiseuille or slit-die flow case (hereafter PPF), and plane Couette flow case (hereafter PCF) in detail and extend our results to other fluid models and to axisymmetric capillary or pipe flow (hereafter HPF). We apply the results of Yerushahni, Katz and Shinnar [13] and will be concerned, as they were, with establishing the instability of certain flows, rather than with proving stability. We extend their results where necessary, always with the same aim. The Johnson-Segalman constitutive equation (with added Newtonian viscosity) is DT

T+h,+

=217 *+A,%

)

i

DT a Dt

= f--LIT+

i

TQ-

a(DT+

0) TD),

where A, is the relaxation time, A, the retardation time [l], 11the viscosity at zero shear, and D and D are the symmetric and antisymmetric parts of the velocity gradient tensor, respectively. The above may be viewed as simple addition of Newtonian viscosity to a Maxwell-like model [14]. Equation (1) is decomposed as T= Tl + T2,

T,+h,%

= 2r),D,

DT

“& = 271,D,

T+X,&

11=771+%,

A2= r/2vv

D+--

(2)

where q1 and q12 are non-Newtonian and Newtonian contributions to the viscosity, respectively, with corresponding stress contributions, Tl and T2_ Thus, the retardation time constant, A,, has the simple interpretation of adding Newtonian viscosity to the fluid model. We employ the following appropriate dimensionless variables,

-P*

T=T,d VU,’

pu,z’

p=

t=-t*uc d

’

The resulting constitutive and momentum equation are DT T-k E+ Du -= Dt

=2

D+EE+ (

-vp+v.T,

DD

(4)

306 where & = qX,/pd2 is the elasticity number, and c = q2/q. The elasticity number may be interpreted as a Deborah number, in which the characteristic fluid process time is chosen as viscous momentum diffusion. Again the

Fixed

x=-h/z

wall

Moving

x=+

wall

h/z

Fig. 1. (a) PCF and (b) PPF, showing actual flow and reduction to 1-D. The system of the xk is a standard rheological system [l]; z is a dimensionless coordinate with associated velocity U, x and u are dimensional analogues of z assumed that h is very much smaller than the other two dimensions. The PPF shows the nodes of a crude, graded mesh of linear elements.

coordinate crossstream and U. It is domain (b)

307 Weissenberg number, Reynolds number, and elasticity number are related by We=Re E. We begin by discussing the PCF case. The geometry of the 3-D flow and its reduction to 1-D are illustrated in Fig. la. The characteristic length, d, is taken to be the channel height, h, and the Reynolds number is based on the velocity of the moving plate. Previous studies have examined the one-dimensional simplification of eqns. (4) which are referred to as the equations of rectilinear shearing motion. In the present case these equations are given by

&( 7

+

zu,) +

(5)

E&Uzt,

24, f

7 =

&(zt+ (cl’- 1)7U,)+Z=~&(a2-1)U,2, where r is the shear stress, Z is a particular linear combination of the normal stresses TX, and T,,, and the subscripts of eqns. (5) indicate partial derivatives. We note that Squire’s theorem [15] does not apply for the Johnson-Segalman model, since the second normal stress difference does not vanish; thus if the analysis of eqns. (5) were to reveal stability, the stability picture would remain incomplete. It is very possible that onedimensional instability could be the responsible mechanism, even in the three-dimensional setting, but our primary purpose is to show what solutions are unlikely, due to their lack of stability. The velocity boundary conditions for eqns. (5) are U(0, t) = 0,

~(1, 1) = Re,

(6)

where Re is based on the upper wall or plate velocity. Initially, we consider the stability of the main branch of steady solutions to eqns. (5) and (6) which are given by u = U(z) = Re z, We (1 + cWe'(1 7 = --g (1 + we2(1 z

=

- ~2’)) _

a2))

=

conSt-,

(7)

(1 - 6) We2(1 - u’)

E(1 + We2(1 - a”))

= constm

Yerushalmi, Shinnar and Katz [13] examined the linear stability of PCF for nonmonotone shear-thinning viscoelastic fluids whose disturbance equations can be expressed in the form .cP(a/at)7”=%(a/at)ii,,

(8)

for polynominals Johnson-Segalman

9 and 22. The disturbance equations for the model can be expressed in this form. The first example

308 of Section 5 of Ref. 13 deals with the special case of our problem in which a = 0. It is an easy matter to modify these results to account for a f 0, from which it follows that solutions for which a?/a( We) > 0 ( < 0) are linearly stable (linearly unstable). Further examination of the roots reveals some interesting facets. The characteristic equation is cubic in So, and there are three roots for a given n denoted by (s,)i, (s,)~, and (s,)+ In the limit as n -+ 00, two of the roots satisfy

(9) therefore,

where B and C are the coefficients in the quadratic equation, es: + Bs,+ C = 0. Simple singular perturbation analysis of the characteristic equation reveals that the third root behaves as

+ o(1). (%A= --c&(n7r)2 Thus, in the presence of Newtonian viscosity (e # 0), none of the roots becomes positive and unbounded as the mode number, n -+ cc, and therefore the instability is not a Hadamard instability of the type considered by Joseph, Renardy and Saut [16]. However, the n = 00 mode grows faster than any other, as is the case in the Hadamard instability. The fact that the Johnson-Segalman (hereafter JS) and upper-convected Giesekus model with added Newtonian viscosity will not exhibit Hadamard instabilities (loss of evolution to infinitesimally short wavelength disturbances) for the PCF case, can also be demonstrated from the calculations of Joseph and Saut [17]. The steady 2-D vorticity equation may change type, but the PCF will not become Hadamard unstable. Numerical examination of the roots also reveals the presence of stable complex eigenvalues with O(1) decay rates. This fact could account for the experimentally observed periodicity occurring at the onset of instability that Yerushalmi et al. cite. We now consider results for the upper-convected Giesekus model [1,2]. The dimensionless form of the upper-convected Giesekus model with added Newtonian viscosity is given by

where a is the mobility parameter,

v the viscosity ratio, and E is the

309 previously defined elasticity number. The one dimensional rectilinear shearing equations corresponding to eqns. (5) are U,=P+vu,,+7r*, u + cx&(02 + T:) + &( a, - 27rU,) = 0, T1+(I1&(U+~)T1+&(Tlt-~:U,)=U,,’ z

+

a&(

7;

+

2:“)

+

EZ,

=

(13)

0

For the PCF case, steady solutions, a, ?r, 3, are determined from the set of nonlinear algebraic equations, a+aE(Z2+?~)-2&7rRe=0,

where the total shear stress, 7 is given by (15)

?=Tl+vRe.

Due to the complexity of the above equations, general analytic solutions for arbitrary (Y are not available *. However, a closed form solution can be found in the special case (Y= 1 and is given by

Linearization of eqn. (14) about eqn. (16) results in the following set of linear stability equations: ii, = viizz + ?I, )

We note that the above set of equations cannot be reduced to the form of eqn. (8), but solutions of the form ii(z, t) = c(z) e”‘, etc. are appropriate, and a boundary value problem of the form, L.

u

II

-

(

h(We;

s

s) + v i

ii = 0,

(18)

G(O)= G(1) = 0,

* In this context, the pertinent values of cx are l/2 I (YI 1 for the upper-convected Giesekus model to exhibit shear stress extrema.

310 where h(We;

3) =

[(1+&~)2-4&2We2~~][(1-We?,&)(1+&s)+2~e&~~-2&2~~(l+We2)] [(1+&s)3+[-4&2~~(l+We2)+4&~~We](l+&s)][1+&s+2We~~&]

’

(19) is obtained. The eigenvalues, s,, are determined from (na)*(h(We;

s)+v)

+s=o.

(20)

Neutral stability requires h(We;O)+v=O=l-We*+v(l+We*)*=j$-,

(21)

which is again in direct accordance with the JS results, though the details of the analysis differ due to the normality in eqns. (12). We note that eqns. (5) or (13) do not exhibit bifurcation in the classical sense [18]. There are not steady finite amplitude solutions with a sin( nrz) structure, bifurcating from the loss of stability point. The PCF solutions simply jump from one branch to the other, very similar to the “snap-through” phenomena in the buckling of an elastic arch (see Fig. 2). This can be seen rather dramatically in numerical simulations of PCF, using the JS model and the techniques which are described in Section 3. We start with a linear velocity profile and consistent r and 2 at t = 0 implying that We, < We < We,. Because of numerical errors in representation, their amplification, and introduction of new rounding errors during temporal integration, in effect we have a small perturbation of eqns. (7). The steady profile is linearly unstable to the perturbation. In the numerical simulation, a double snapthrough ensues, as is illustrated by Fig. 2; in part of the interval z E (0, 1) the fluid jumps to a state with We, and part to a state with We,, resulting in a layered flow with two shear rates and jumps at the interfaces between them. The velocity is continuous at the interfaces, and the We based on plate velocity is achieved. But locally in the flow, this We is so unstable as to be practically unattainable. The general observation is that, while the unstable flow may look quite steady, the dynamic process and numerical errors supply enough “kick”, so that unstable states are never achieved in long-time simulations. We note that if the velocity boundary condition on the plate in eqn. (6) is replaced by a stress condition (corresponding to “problem C” of Ref. 13), the snap-through is from one homogeneous shear flow to another at constant stress. The double snap-through to a layered solution is more relevant to our analysis of PPF. We also note that there may be bifurcation in the classical sense if the two or more complete three-dimensional equations are considered.

311

We

Fig. 2. Double-valued inverse of a non-monotonic shear stress. Arrows converge on stable steady states and diverge from unsteady ones. Solutions with local values of We = We, snap through to We, or We, locally.

The numerical simulations show that nonclassical, weak solutions are possible and are achieved as an alternative to linearly unstable flows with the same imposed plate velocity. The complete analysis of weak solutions and their stability takes us beyond the scope of this paper; however, there is an obvious extension of our results to weak solutions which have a finite number of jumps in shear rate and constant u,, 7, Z, etc. in between; T is given by eqns. (7) between jumps with We based on ~localshear rate, and has the same value in each layer. The stability analysis holds piecewise between jumps (using an appropriately shifted velocity profile for the base flow), and we reach the same conclusion: Solutions with local shear rates corresponding to the descending portion of the curve of Fig. 2 are so unstable as to be practically unattainable. A model, which is very similar in spirit to the Johnson-Segalman model, is one considered by Hunter and Slemrod [ll] * . The one-dimensional rectilinear shearing equations are given by &7,+7=&a(u,), u,=P+7

P)

(22)

where P is the dimensionless pressure gradient (P = 0 for PCF) and a( uz) is a nonmonotonic function with the same shape as Fig. 3. Thus, our models have qualitatively identical steady shear behavior. However, this model

* Hunter and Slemrod actually considered the PPF case, but the results for the PCF case can be readily inferred from their analysis.

312

Fig. 3. Schematic of a non-monotonic steady stress with two local extrema, such as occurs with the JS model with Newtonian viscosity. Shear stress, T, is plotted versus a dimensionless measure of the shear rate, We. The critical values are for the two local extrema.

exhibits Hadamard instability, in the absence of Newtonian viscosity, between the shear stress extrema, in contrast to our JS model. We also remark that the Hunter-Slemrod model exhibits the jump phenomenon in absence of Newtonian viscosity, in contrast to our simple JS model. More complicated JS models can be constructed with multiple nonzero relaxation times which are similar in this regard to Hunter-Slemrod. Reference 16 applies directly to the linearized equations, and we deduce that (s&

= - &[l

+ {l - 4E20’(IVf?)(n~)2],

(23)

and for large n, n-,cQ. (24) (%l)* = n?r{_, If a’( We) < 0, the flow is unstable, and one of the growth rates becomes unbounded as n + co. This loss of evolution is avoided in Ref. 11 by the addition of an arbitrarily small amount of Newtonian viscosity. We now consider the linear stability of PPF. The geometry of the 3-D flow and its reduction to 1-D are illustrated in Fig. lb. The characteristic length, d, is taken to be the channel height, h, and the Reynolds number is based on the centerline velocity. The dimensionless one-dimensional rectilinear shearing equations for the PPF case, of the Johnson-Segalman model (with added Newtonian viscosity), are u,=P+T,, &(7+ 224,) + 7 = u, + C&U,,, &(z~+(az-l)rU,)+Z=e&(a2-l)u,2,

(25)

313

with boundary conditions, u(*l,

7(0, t) = 0.

t) =o,

(26)

The only steady solutions of eqns. (25), subject to boundary conditions (26), are such that 2.4= u(z),

7 =.Y(u’),

Z=Y(U’),

(27)

where Y( U’) and Y( U’) at any point, z, are as in eqns. (7), with We replaced by EU’. Note that 7 and Z are thus functions of z in PPF. U(z) satisfies G(U’,

P;&, c, u) = e&'(l

- LZ’)U’~ + ~~(l-

a2)PzU’2

+ U’+

Pz = 0. (28)

Note that this result holds as written for the class of weak solutons with a finite number of jumps in U’. We have employed a Chebyshev polynomial expansion to obtain solutions of eqns. (25) subject to boundary conditions (26), both below and above criticality. The program was checked by corroborating the results of van Schaftingen and Crochet [19], who considered the HPF case of the above problem, where E 2 l/9, i.e., the model simply shear-thins with no shear stress maxima or minima. The steady shear stress function 7 is as in Fig. 3, where the critical values of the shear rate, at the maximum and minimum (corresponding to We, and We,, resp.), are given by

m

&‘=_

u;=_

m

JZ&(l - a2)1’2 ’ f*(E)

= f

-

O&(1 - a2)1’2 ’

(29)

3 f (( + - 3)2 - ,)I”,

and the corresponding

values of the pressure gradient are

1 f-(4 pc= r d 2(1 -u’)

(1+ ef-(4/2) (1 +f_(4/2)

1 f+(c) pm= e J--2(1 -u’)

(1+ ef+(4/2) (l+f+(e)/2)

’

(30) *

We first briefly examine the Maxwell limit, (e = 0), for which there is only a

Fig. 4. Schematic of a non-monotonic shear stress with an absolute maximum, such as occurs in the JS model without Newtonian viscosity. The value of We, indicated is for the JS model.

shear stress maximum; this is illustrated in Fig. 4. In this case an exact solution to eqn. (28) can be obtained and is given by

Jl /l -

1-

2E2(1 l- a2)P

u(z)=

+

2

1+

[ L In i

1 - 4E2(1 - a2)P2z2

-i

In i

4E2(1 4E2(1 - a2)P2z2 a2)P2z2

- In 1z 1

1-

1 - 4E2(1 - a2)P2

1+

1-4&*(1-a*)P*

:

-

-/l -

II’

1

4e2(1 - a2)P2 (31)

The functional dependence of We on P can be readily obtained from eqn. (31) as We=

1 2&(1-

-$

a2)P

In

1 - + ln( E2P2(1 - a’)) [

1 - Jl - 4E2(1 - a2)P2

-

Jl -

I 1+ Jl - 4e2(1 - a2)P2

4E2(1 - a2)P2

1 (32 .

We observe that classical steady solutions do not even exist for values of the pressure gradient and Weissenberg number beyond PC=

1

2&(1 - .*)l’*

we

’

=

c

l-h-l2

(1 _ u*)i/* *

(33)

Numerical solutions without Newtonian viscosity, using the random choice [20] method, show the development of a ‘slug-like’ flow with a flat profile

315 and unboundedly growing shear rate in a very thin layer near the wall. A physical interpretation-if there is one-is not clear at this time. Addition of Newtonian viscosity results in existence of steady solutions for all values of the pressure gradient. There are unique C” solutions for P -c P,,,, and for P > PC the solutions will have a jump discontinuity in U’. For values of P such that Pm I P s PC,jumps may or may not occur. Linear stability analysis of the solutions to eqn. (28) results in the following second order non-constant coefficient eigenvalue problem for the spatial part of the disturbance velocity, G, (g(7-7,

s)a’)’

- sii = 0,

G(1) = z?‘(O) = 0,

(34)

(1+&s)2+(-l+&s(l+r)+&2s~+3~)&Z(1-~~)U~~+~(1-.~)~&4U~4 g= [(l+ &s)*+&*(1-2)uq1+

&2(1-a2)U'2]

.

Multiplying the first of eqns. (34) by ii and subsequent integration by parts results in the nonlinear Rayleigh quotient s = - J(olg2’ d#0”

dz.

(35)

Taking s = 0 in eqn. (35) results in the following condition on U’, 1 + (3e - 1)E2(1 - a2)U’2 + eE4(1 - a2)2u’4

= 0,

(36)

which simply requires that U’ = UC’or U,,!,.Linear stability theory thereby states that only steady velocity profiles which are composed of positive sloping portions of Fig. 3 are stable to infinitesimal one-dimensional disturbances. For values of P such that P >, Pm, there is an uncountable infinity of possible solutions not ruled out be our analysis; see Fig. 5. If the values of P is say P*, then the possible profiles are those which jump at any value of P in between Pm and P*. When a jump does occur while the driving gradient is being increased, we will associate this with a phenomenon which we will call “Vinogradov’s spurt” [6]. The actual fluids Vinogradov et al. [6] experimented with were monodisperse polymer solutions of polybutadiene and polyisoprene; we will discuss this phenomenon and give data in Section 3. We shall define spurt as the nature of experimental observations or mathematical solutions which is characterized by a period of time during which the throughtput increases dramatically (the dynamic spurt process). Or alternatively spurt is characterized by a dramatic increase in steady throughput as a function of driving gradient (the asymptotic steady consequences of dynamic spurt). These flows or solutions have two dimensional surfaces across which the shear rate jumps discontinuously, which, after Ref.

316 (a)

(b)

z

Z

Fig. 5. Profiles for two of the infinity of PPF solutions with a fixed driving gradient, P*. At the wall, 7 = hP/2. Profile under column (a) does not jump; profile under (b) is a “bottom jumper”. The dashed profiles are Newtonian parabolas for the purposes of comparison.

11, we shall call singuhr surfaces. The 1-D solutions have the discontinuities in U’ mentioned earlier. We shall refer to such flows or solutions as singular when jumps occur and nonsingular otherwise.

Fig. 6. A profile with P/P, =1.00014, nian parabola for comparison.

showing “apparent slip.” Dashed profile is a Newto-

317

Fig. 7. Top jumping solutions under loading conditions. 6 increases, and x, moves to the right as loading progresses. Computational coordinate has origin shifted to the left-hand die wall, and only the left half-channel is plotted. Figure is a schematic because for Vinogradov’s data, the critical solution would not be on scale with solutions which show significant layer thickness.

Figure 6 shows a case where the jump has occurred at the top (shear stress maximum). The velocity profiles below and above the critical pressure gradient, PC,are shown. We see that the post-critical solution is composed of an extremely thin high-shear wall layer together with a parabolic-like core. In Fig. 6, the pressure gradient is only very slightly over critical, and the

Fig. 8. Shape memory in unloading (from curve A in Fig. 7); x, stays fixed until only bottom jumping is possible at curve B. Curves below B are bottom jumpers; curve A’ is an intermediate jumper. x, moves to the left during bottom jumping. This Figure is a schematic for reasons cited under Fig. 7. Computational coordinate is the same as for Fig. 7.

318 Unloading 474

I

P = 3.0 P = 2.7

P = 2.4 P = 2.0 U 237

P=l.O 118.E &= 1.0, E = 0.001, a = 0.9793

0 -1

I 0

Fig. 9. Actual computed curves for representative data showing fixed layer position in unloading, representing “shape memory.” Profiles in which layer becomes smaller after further unloading would not be easily visible to this scale (see schematic of Fig. 8).

resulting layer is very thin. How layer thickness depends on pressure gradient will be detailed in Section 3, if it is not already obvious. Two hundred Chebyshev polynomials were used to resolve the wall layer. Such a phenomenon is generally referred to as “apparent slip”, even though there is no actual slip, due to extreme thinness of the wall layer, it may appear as slip to the experimenter. However, as the pressure gradient is increased, the wall layer thickens and eventually envelops the entire region. This aspect is depicted in Figs. 7-9 for solutions jumping at various different points on the curve of Fig. 3, using representative data described below. Additionally, we note that Joseph and Saut [16] have shown that the PPF case of a Johnson-Segalman model with e > 0 will not exhibit the Hadamard instability, associated with change of type. We also note that the conclusions from our present analysis are valid for the HPF case. The axisymmetric form of the one-dimensional shearing equations for the JS model are given by

T+&(~+ZW,)=W,+e&w,,,

(37)

Z+&(z~+(a2-l)7W,)=~&W~, w(l,

t) = 7(0, t) = 0,

where w( r, t) is the axial velocity (in the z-direction). The above system is virtually identical to the planar case. The only difference is in the form of the r-derivative on the shear stress, 7, in the first equation. The stability

319 analysis proceeds as previously described, and a nonlinear Rayleigh quotient analogous to eqn. (35) is obtained and is given by -

.S=

J

o1rg(W’, s)( G’)2dr 1 J0

+*r

dr

(38)

3

where g( IV’, s) is the same function that appears in eqn. (35). Thus we have exactly the same stability results as in the PPF case. Thus the spurt instability is not geometry dependent, as is the case for the inertial instabilities, which has been detailed for the Newtonian case by Patera and Orszag [2Il. To summarize the conclusions of our stability analysis, the most important observation is that in either PCF with an imposed plate velocity, or in PPF, solutions which would have shear rates on the descending portion of a non-monotonic constitutive curve cannot ever be expected to occur in the real world-at least if the real world is reflected in any of the several models we have analyzed. It is therefore irrelevant to dismiss such a model because it implies that stress decreases with increasing positive strain rate: it does not so imply. Descending portions of the curve are avoided in these flows by the formation of layers in which the shear rate “jumps across” the inadmissable range, while the total physical stress remains continuous, as it must to balance forces. Furthermore, though stability analysis does not identify which stable solutions are obtained for such a model, it shows that the dynamic equations retain evolutionarity, and thus the dynamic problem appears to remain mathematically well-posed “beyond the shear stress maximum.” Thus there are mathematical implications of such constitutive equations which bear investigation before their physical validity can be assessed. Finally, the instabilities occur in a similar fashion in slit-die or capillary geometry and appear to depend on material properties, not geometry. 3. Dynamic behavior of the Johnson-Segalman

model-Numerical

results

The construction of steady solutions and stability analysis of the previous sections give us no indication of which of the infinitude of solutions is selected when the dynamic system is given rheologically reasonable initial data but only of which steady solutions will not be selected. This is a critical issue to resolve, as is illustrated by Fig. 10, where We is plotted versus P using the techniques of Section 2 to obtain steady solutions for data representative of Ref. 6. Top jumping and bottom jumping solutions can lead to significantly different curves. To explore this question further we

320 92

We 46

0

Fig. 10. “Bottom-jumping” branch, emanating from Pm, and “top-jumping” branch emanating from PC. Flows which have driving gradients larger than P,,, could be on either branch or in between and consequently differ markedly. Dynamic information is required to determine where in between these extremes actual solutions fall.

turn to the Johnson-Segalman model and numerical simulations of PPF. In a forthcoming publication, Malkus et al. [20] analyze the dynamic behavior of this model in general shearing flows; they focus on the mathematical nature of the boundary-initial value problem, the analogy to problems in gas dynamics, and the ramifications for numerical methods. We present some numerical results here, calculated by one of the methods of Ref. 20, with the aim of establishing what an experimenter might actually see if the phenomena described in the previous sections represent physical phenomena rather than mathematical artifact. To facilitate comparison with Ref. 20 and with experimental results, we use the dimensional form of eqns. (25). We also will find it useful to take the Newtonian viscosity contribution out of the stress equations and put it into the momentum equation, leaving the stress equations in terms of the non-Newtonian contribution only, as in the second of eqns. (2) and in eqns. (12) and (13) for the Giesekus model. This is straight-forward, and the result is pu, = 0, + vu,, + P, q=(Z+p)u,-ACT, 2, = (a’-

l)au,

(3% - AZ,

where u(x, t) is the velocity field, x E (-h/2, + h/2) is the dimensional analogue of z, the cross-stream coordinate; u is the non-Newtonian part of the shear stress (which corresponds to the 1-3 component of Tl in eqn. (2)),

321

P the applied pressure gradient, h = l/h, the reciprocal of the (single) relaxation time, p/X the non-Newtonian viscosity at zero shear, v the Newtonian viscosity (dimensional analogue to v of Section 2), and 2 has dimensions of stress and is similar to the quantity introduced in eqns. (5) (but has the opposite sign in steady flows). The parameter - 1 5 a I 1 is nondimensional as before. The v-term is the Newtonian viscosity term, moved to the momentum equation. When v = 0 the system is hyperbolic with wave speeds 0, f {(z ; with v # 0, the system is not easy to classify, but experience shows that wave propagation and many aspects (but not all) of hyperbolic character are retained if v is small enough. We will always take v small or zero, and not insist that it be taken large enough to suppress the ‘shear stress maximum’ phenomenon. The numerical method we employ finds its origins in a Newmark-P [22] method with implicitly treated Newtonian viscosity, an implicit/explicit treatment of the constitutive equation, and no strain-stiffness term. The diplacements can thus be dispensed with, and with them the p parameter; this method is then very much in the spirit of Hughes-Liu-Brooks [23]. We will use linear 1-D finite elements [22], ‘line elements’, which are illustrated in a crude graded mesh in Fig. lb. As in the figure, only the left half of the channel is discretized, or ( -h/2,0). In the computations reported here, 40 elements were used in a graded mesh with the smallest elements nearest the wall. In Ref. 20, evidence is given to support our contention that the resulting level of accuracy is sufficient to warrant the conclusions we draw here. We define the usual finite element ‘B-matrix’ [22], for linear line elements. For element. e this is

[&?I= ~H,ll,

(40)

e

where the subscript e indicates that a nonuniform mesh may be employed. The discretization, in terms of nodal values of u(x, t), of eqns. (39) is

[M](ic”+‘}

+v[K,]{u”+‘}

+ {Fg’}

= {F”+l},

(41)

where [M] is the usual mass matrix, and {F”+‘} is a consistent load vector based P [22]. The matrix [K,, ] is the Newtonian dissipation matrix evaluated with constant, unit viscosity. { FvTz’} is the dissipation due to the viscoelastic extra stress; in the dynamic case it cannot be computed directly from the shear rate. Instead, we define auxiliary stress variables, u,” and 22, for each element which are updated at each time step by forward differencing the incremental form of the constitutive equations. We proceed by splitting the viscoelastic term into two pieces in the momentum equation to result in

[M]{ti”+l } +(v+Aty)[&]{~“+~}

= {..+l}

- {&?}.

(42)

322 The Atp-term is the portion of the incremental stress in eqns. (39) due to pu,, which will be treated implicitly in algorithm (42). { fi$,z’ } is the remaining part of the viscoelastic dissipation, predicted from the previous time step, so that it is assembled from element vectors

J[B,I ‘((l-Ath)u,“+Z,“[B,]{u,“}) e

dx,

(43)

where { ue} is the element nodal-value vector. The velocity approximation of the Newmark method is used: {U ‘+l}

= {u”}

+ At@

- y){ P}

+ h{ P+l}),

where y 2 l/2 is a Newmark parameter. This equation is solved and substituted into eqn. (42)-which is linearly implicit but explicit-to give { u’+l} in terms of known quantities, starting data. The stress is then updated using the latest values of according to -~A~)~J~+A~(Z,“+/L)[B,]{~“+‘}, 0,‘+I=(1 Z;+l= (1 - XAt)Z,n + (a* - l)u,“+‘[ B,] { zP+l}.

(44) for { tin+‘) nonlinearly from initial strain rate,

(45)

Simple matrix stability analysis in the linear case (i.e., a = + 1) of this algorithm shows that there is no elastic wave speed restriction on the time step but there is a restriction relating At and l/X which proves to be of no practical consequence, since accuracy seems to require that At be significantly smaller than l/A. Numerical experience shows that there is a nonlinear stability limit when a # + 1, but the practical time step can be as large as 0(106) element wave traversal times (based on wave speed at zero shear). The method can compute solutions above and below the critical stress for the shear stress maximum. For example, in a computation with p/A = 9.8 X 107, l/X = 159.5, a = 0.98, Y = 0.00142 p/h, and h = 96 mm, with a pressure gradient of 95% of critical, this method is stable with At = 5 X 10’ element wave traversals, using 40 finite elements in the half channel. A steady solution is obtained in fewer than 32 time steps, which are of negligible computational cost. The method is every bit as robust in computing solutions with post-critical driving gradients. In such problems, the physical process of the spurt phenomenon takes several relaxation times to evolve, and more time steps are required to obtain post-critical solutions with singular surfaces. 3.1. Vinogradov’s polyisoprene data Vinogradov et al. [6] report spurt phenomena in capillary flow for certain polybutadienes and polyisoprenes at room temperature. We chose to model

323 the latter for several reasons not relevant to the discussion at hand. The polyisoprene samples are labelled PI-l through PI-8, and their molecular weights are given in Table 1. The samples are fairly monodisperse, with polydispersities [l] ranging from 2.03 at the lower molecular weights to typical values of 1.1 at the higher molecular weights. Vinogradov and co-workers found that the zero shear viscosities and dominant relaxation times obeyed the following power laws: q=BMa, l/X = CMB. Values of (Y= 3.2 and p = 3.3 can be determined from a logarithmic regression line. Simple rubbery liquid theory suggests that (Y and p should be equal, and experiment suggests that the value should be about 3.4 [l]. As a compromise, we took (Y= /3 = 3.3 for simplicity and adjusted B and C to minimize the mean square residuals. To two digits we find that B = 9.6 X lo-=, and C = 1.6 X 10-17, when-following Vinogradov et al.-cgs units are employed. Thus we have a simple model of Vinogradov et al’s data which actually fits it rather well. From this point on, we define “Vinogradov’s data” to be the power laws of eqns. (46), with the given values of the constants evaluated at the molecular weights given in Table 1. The fit has a maximum error on the order of 50% in relaxation times for PI-2 and 3 but has less than 12% in error for fluids PI-4 through 8. The viscosities fit the given power laws quite well indeed, as may easily be verified by the reader. So our modelled Vinogradov data is simple to describe, and we believe it matches the data to a level commensurate with the precision of the other modelling assumptions. The p of eqns. (39) is an elastic modulus, and by eqns. (46) is a constant, p = 6 x lo5 dyne/cm* (60 kPa) for all samples. Next we determine v by observing that sample PI-2 appears not to spurt, but PI-3 does. We presume

TABLE 1 Vingradov et al.‘s raw data on which power laws are based (Molecular weight, M = aW, dominant relaxation time, l/h, and polydispersity, aW/E” are taken from polyisoprene data of Ref. 6 Fluid M x10-’

1)x1o-6 A*= l/h M,/M,

PI-1

PI-2

PI-3

PI-4

PI-5

PI-6

PI-7

PI-8

1.06 0.480 0.8 2.03

1.48 1.26 1.3 1.61

2.82 10.00 10.0 1.14

3.55 19.95 1.14

3.80 28.18 46.8 1.10

4.22 39.81 66.1 1.05

5.75 100.00 1.02

6.02 109.60 182.0 1.10

324

I 4

’

I 5

’

1 6

’

LOGt (dyne /cm2 ) Fig. 11. Flow curves from Vinogradov et ah’s spurt data for capillary flow. Nominal shear rate, based on throughput, is plotted against nominal shear stress, based on driving gradient. Fluid is high-molecular weight polyisoprene at 22O C. The figure represents a tracing of an enlargement of the figure of Ref. 6; where the points are tightly clustered, it was not possible to identify the fluid symbols with complete certainty.

that we have e = l/9 for PI-2. We should point out that we do not believe that there is any solvent or solvent-like action involved here. The Newtonian viscous term in the current study represents the effect of a relaxation time much shorter than the dominant relaxation time, modelled here as a zero relaxation time. Modelling the short relaxation time by a Newtonian viscous term can only be approximately correct in a restricted range; at higher stresses the modelled curve and Vinogradov’s data in Fig. (11) would diverge. As is borne out by Fig. (12), we restrict our attention to a range of stresses in which the model and data for PI-2 are in close agreement. Under the assumption that PI-2 (raw data) is just subcritical, we deduce that Y = 1.4 x lo4 Pas and presume that this does not vary with M. Note that this is a pure presumption which can only be justified by our resulting ability to reproduce the spurt effect with the molecular weight dependence observed by Vinogradov et al. Note also that this is a reasonably large value of v, consistent with the hypothesis that the Newtonian viscosity actually represents some macromolecular process. Nevertheless, the viscosities of fluids PI-3 through 8 are so high that e ranges between 0.01484 and 0.001433 for these fluids (modelled by eqns. (46)). Since Vinogradov et al. do not specify the density of the test fluids, we assume that p = 1 g/cm3

325

Fig. 12. Simulated flow curves for slit-die flow, using Johnson-Segalman rectilinear shear flow model. Nominal shear rate (planar) is plotted against nominal shear stress in analogy to Fig. 11. Coalescent points are ones at which one or more data points fell in the same position, to graphical accuracy.

independent of M. This seems to be a reasonable materials ([l], p_ 207).

assumption

for these

3.2. Fitting Vinogradov’s spurt data Under the above assumptions, E turns out to be so small that its effect on

PCof eqns. (33) is negligible. This simplifies the modelling of Vinogradov et al.‘s spurt data, as we shall now explain. Vinogradov et al.‘s key observation is that the critical PC (dimensional) for spurt to ensue is independent of M. This is illustrated in Fig. (ll), where an nominal shear rate, S, based on volumetric flow-rate is plotted versus nominal shear stress. Here we use slit-die flow (PPF) instead of capillary flow (HPF), and the corresponding definition of S is u(x)

dx.

We correspond capillary diameter of Ref. 6 to die height and use h = 96 mm. Note that a correspondence between dimensional pressure gradient, P, and nominal shear stress is given by

rwall= u + vu, = hP/2,

(48)

in steady PPF. Henceforth, when we refer to nominal shear stresses, we shall

326

by implication, mean stress at the die wall unless specifically indicated to the contrary. Using eqns. (3) in eqns. (29) results in dimensional analogues u; and u,’ to UC’ and UA, respectively, which leads to two critical nominal stresses, rtop= T( u;),

and 7bot= r( u,‘).

(49)

Asymptotic expansion of the critical shear rates (taking positive values, without loss of generality) in the indicated powers of c leads to u, -X(1 - a2)-1’2 + O(E), u,’ - e-112 u, + 0(r”2), P Top

-

rbot

-2

2(1

_

J

(l_a2)

a2)1/2

+

o(4

CLVX +w’2)-

Thus “ top-jumping” critical stress, rtop, is essentially independent of M; to leading order, it is the critical stress in what is referred to as the “Maxwell limit” in conjunction with Fig. 3 and eqns. (33). For purposes of comparison, it should be noted that in Ref. 20, E = VA/~; this definition leads to a critical value of l/8 for shear stress extrema but is asymptotically equivalent to our definition here for small values. Equations (50) are the same with either definition. The fact that the observed critical stress is independent of molecular weight suggests it may depend on the elastic modulus, p. This observation is also consistent with rubbery liquid theory [l]; there the elastic modulus is due to temporary entanglements. The near constancy of the modulus over a wide range of molecular weights in Vinogradov et al.‘s raw data suggests that the density of these entanglements depends more on the nature of isoprene than the lengths of the chains involved in the polymerization of PI-l through 8. The fact that the critical stress, rtop, in our model is also independent of M (to O(C)) when a is independent of M, suggests using rtop to determine a. By ignoring terms of order c and setting 7t0p of equations (50), to the observed value of 1.5 X lo6 dyne/cm2 (0.15 MPa) [6], we obtain a value of a = 0.98. While this seems to be a high value of a in terms of what Johnson and Segalman believe to be typical [12], we shall see shortly that this value gives very little shear-thinning below u, , a feature which can be observed in Fig. (ll), and is pointedly remarked upon in Ref. 6. Note that the “bottom-jumping” critical stress, Tbot, is dependent on M, and this is a crucial feature of our model which plays an important role in the predicted hysteresis. The error involved in the above expansions decreases

327 with c and thus M; for fluid PI-3 the expansion is accurate to about 4% for 7t0p and better than 1% for T,,~~and is lower for the remaining critical fluids, PI-4 through 8. We will employ these asymptotic expansions for the critical numbers of eqns. (SO) hereafter. We have presumed that Vingogradov’s observed critical stress is given by rtop. In the light of the constructions of Section 2, this assumption seems unwarranted, since solutions with layers which lead to spurt phenomena can be constructed with any critical stress between 7t0p and rbot. In fact, the model of McLeish and Ball is based on the presumption that for the Doi-Edwards fluid employed, the stress at the singular surface falls from 7t0p to rbot during the dynamic spurt process; this comprises a key difference between the two models. Our own presumption can be argued to be valid for the Hunter-Slemrod model of eqns. (22) [ll]. We shall use numerical dynamic simulation to argue that our more complicated model acts in precise qualitative analogy to the Hunter-Slemrod model. First some terminology. We define loading to be occurring when P dP/dt 2 0, which includes start-up. We define unloading to be occurring when P dP/dt < 0, which includes buck-off. Start-up and back-off have, for constants PO and PI, P=P,+P,H(t-

to),

(51)

where for start-up, POPI 2 0 or P,PI < 0 and 1PI 12 2 1PO (and for back-off, P,PI < 0 and 1PI 1< 2 1PO );in either case, the initial data at t = t, for eqns. (39) is consistent with a steady solution with P = PO,and H(t) is the unit step function. Note that the zero solution is the unique steady solution with P,, = 0, and our definitions imply that starting with zero initial data and nonzero P, is the start-up, regardless of the sign of P,. Our key numerical finding is that 7t0p is the dynamically chosen critical stress in start-up experiments (and loading in general). This conclusion holds in simulations which start from Pa = 0, in which P, determines the final nominal stress, as well as simulations which proceed “stepwise”, in which PI is the step increment, and the flow is allowed to steady out between steps. In Fig. 12 we plot the results of the simulation of fluids PI-2,3,4, and 7 in start-up experiments covering a range of nominal stress from zero to slightly more than two times rtor,. The data for the remaining PI fluids is omitted in order not to clutter the plot, but the agreement is every bit as good for those fluids. Note that there is little evidence of shear-thinning until the spurt occurs. One difference between our model’s predictions and Fig. 11 is that our model predicts that all flow curves will tend to coalesce into a single, M-independent curve at nominal stresses significantly higher than critical. There is more scatter in Vinogradov et al.‘s curves than there is in our model; however, Vinogradov et al.‘s curves to tend to cluster in this region,

328

with much less molecular-weight-dependence than at lower nominal stresses. One explanation for the scatter may be that-as is borne out by our simulations-the spurt process for the modelled fluids involves surprisingly long time-scales. Reference 6 refers to the spurt process as unsteady; if our model is a good one, we predict that in principle, steady, spurted flow can be achieved, though this might pose significant experimental difficulties. If such steady flows can be achieved, our model predicts M-independent throughput at sufficiently high nominal stress. How this matching of the spurt data fits into a broader picture of the dynamic process which includes our predictions of unloading behavior will now be explained. 3.3. Major conclusions 3.3.1. Latency, spurt and arrest There is a “latency period” after start-up from zero during which the velocity profile looks like a subcritical one; the pseudo-steady profile for PI-7, started from zero at 1.2 7t0,, and sampled at t = 300 s is shown in Fig. 13a. It should be emphasized that, though the flow appears steady during latency, it is not. (Though a latent flow could deceive a numerical termination criterion which is not tight enough.) We can be sure of this: No steady, nonsingular solution of eqns. (39) exists at nominal stress above 7_,, and the solution of Fig. 13a is surely the finite element approximation to a nonsingular solution [20]. The only alternative is that the solution of Fig. 13a is not steady. Close inspection of latent solutions confirms this-there are small changes with time which are numerically significant. The latency period depends on M and the fraction of critical nominal stress which is imposed; in this simple model, latency evidently could be made to last arbitrarily long, but in the real world or more complicated models finite amplitude disturbances or infinitesimal two-dimensional disturbances would, before too long, kick the system to the next phase: spurt. As illustrated in Fig. 14, spurt is a period during which the throughput grows rapidly with time; nevertheless, since the throughput can increase by orders of magnitude, it can take many hundreds of seconds to become essentially steady. Figure 13b depicts the same experiments as Fig. 13a, but samples during the spurt at time t = 500 s. The singular surface is not fully developed, as is indicated by the roundness near the corner at the developing singularity. The finite element solution of Fig. 13b does show a jump in the shear rate because the finite element method described earlier approximates shear rates with piecewise constants; these jumps can be large when the gradients of strain rate are large. Mesh refinement [20] shows a definite roundness in similar flows. This observation is in accord with the analysis of Ref. 20 which argues that the surface becomes truly singular only in the

329 700x10*

I

I

I

600 500 Ll (A)

0.00

0.01

0.020

0.02 )( 0.03

I

I

I 0.01

I 0.02

I

I

0.01

0.02

0.04

I

I

I 0.03

I 0.04

0.05

0.015 w

u

0.010 -. 0.005 -

o.oM) 0.00 0.10

CC)

x

I

0.05

I

U

0.00

Fig. 13. data for of scale same as 10 -

X 0.03

0.04

0.05

Velocity profiles during the three phases of the simulated spurt phenomenon, using fluid PI-7: (A) latency, (B) dynamic spurt, and (C) arrest. Note that there is a change over several orders of magnitude from (A) to (C). Computational coordinate is the for Figs. 7 and 8. I

I

I

86S t

4-

2-

0.

0

1000

2000

I 3000

4000

t

Fig. 14. Nominal shear rate as a function of time in seconds, during the simulated spurt of fluid PI-7. The plot is scaled by a short burst of high flux in which Newtonian viscosity dominates. This cannot be seen on the time-scale of the plot. The flat portion of the curve is latency, follwed by dynamic spurt, characterized by high frequency oscillations. As the spurt arrests, the funnel-shaped envelopes which appear are numerical in origin; they are generated by progress of layer from one element to the next.

10

S

a 6 4

1: 0

100

200

300

400

500

Fig. 15. Detail of simulated spurt for PI-7 at early time, using l/S th the timestep of Fig. 14. The initial burst is better resolved but still not easily visible. Oscillations at the initiation of dynamic spurt are resolved with many points per period, showing that they are of lower frequency than numerical sampling and thus may be real.

infinite time limit. But as shown by Fig. 13c, it is essentially fully developed by time t = 4000 s. This regime of fully developed singular flow with two distinct layers represents the final phase of the dynamic process: arrest of the spurt. Since it would take so long to develop, and might be subject to a variety of secondary instabilities (including, perhaps, true wall slip), this arrested spurt may be very hard or impossible to observe. A few words about the character of Fig. 14 are in order. The dark irregular shape of the curve is due to high frequency oscillations of the calculated S. They appear to be at a frequency close to that of the numerical sampling frequency [22]. Closer inspection shows they are of lower frequency, as is illustrated by using a timestep of l/8 th that of Fig. 14 to generate Fig. 15. We conclude that, as the singular surface moves and develops, there may well be real oscillations in S. On the other hand, the funnel-shaped oscillation envelopes which begin to appear at about t = 1500 s in Fig. 14 are surely numerical in origin. This matter will be dealt with more fully in Ref. 20, but in summary, the funnels result from the fact that the motion of the singular surface is ‘quantitized’ by the spatial discretization and the fact that a singular surface must coincide with an element boundary. A funnel is generated each time the surface ‘pops’ past a finite element. We believe, at the very least, that Fig. 14 is reflective of the mean trend. 3.3.2. Hysteresis Here we consider a loading sequence followed by an unloading sequence in which we perform a series of startups from one stress level to the next,

331 followed by a sequence of back-offs. The following observations apply to essentially steady solutions, in which the singular surface is sharply defined. We start with P,, = 2r0/h and let PI = 0.2 X P,,, letting I’,, be replaced by P,, + PI in each succeeding cycle. At each load we solve the dynamic equations and let them steady out, using the previous solution as initial data. The solutions behave like those previously obtained during loading. That is, subcritical parabolic-like solutions are obtained until P > PtO,, = 2rtop/h. From then on top-jumping spurt solutions are found. For the asymptotic steady solutions r = Px, so during loading only, the thickness of the high shear rate layer is

(52) where 7imp is the nominal shear stress associated with P. The numerical solutions typically put xc, the location of the singular surface, at the nearest element boundary to the exact x,. A schematic of such a sequence is shown in Fig. 7; a schematic is required as it is hard to show motion of the singular surface and keep the velocities on the same scales. During unloading, i.e. when P is being decreased by factors of 1.2, from are two phases in the some value, rreV-- that is we ‘back-off’ from T,,,-there numerical results: a) Unloading with fixed 6, where

Thus, if we denote

the stress at the singularity

Ph rtop rjmp=pxc=

2

7,,,

=

‘imp’top 7,,

by 7jmp,

.

b) Bottom-jumping unloading. When the value reaches rbot it can go no lower, so that when 7bot ’

TimPTtop ~

7rev

(54) of rjm,, in unloading

(55)

’

or when 7imp satisfies 7bot<-

7rev

TimP
7rev

7top

’

(56)

then 7jmp= rbot, and 6 + 0 with 6 = 0 when rimp = rbot. For lower 7imp, subcritical solutions are obtained. A schematic of such a sequence is pictured in Fig. 8. A portion of the back-off in which the layer does not move can be plotted from actual computations and is given in Fig. 9.

332 161

I

I

I

I

I

I

1412IO-

S

B6-

Fig. 16. Simulated flow curve for PI-7 through a loading sequence (start-ups) and unloading sequence (back-offs). Spurt starts at 7top and the loop opens at the top during the first back-off. It closes at TV,. This figure represents the signature of shape memory hysteresis.

Note that for a bottom-jumping

h --

rbot

Xc = 2 Timp.

solution we have

(57)

The numerical solutions show that for the first value of 7impcorresponding to eqn. (56), eqn. (57) is satisfied to computational accuracy. The phenomenon of (b) is called “shape memory” by Hunter and Slemrod [ll]. Figure 16 illustrates a hysteresis loop for fluid PI-3, loaded to 7,eV= 2.629 x lo6 dyne/cm2 and unloaded until the loop closes. A salient feature of shape memory hysteresis is that the loop never retraces itself until the main branch is rejoined, but always opens from T,~“.It moves from the top-jumping branch of Fig. 10 to rejoin the bottom-jumping branch at some point. Note that Fig. 10 plots We and Figure 16 S, but the pictures are directly comparable. The shape of the loop is a consequence of the fact that x, always moves toward the centerline of the die during loading but always stays fixed for some fraction of an unloading sequence. This feature of our model is identical to what is observed by Hunter and Slemrod [ll], but in marked contrast to the model described by McLeish and Ball [lo], in which the supercritical part of the flow curve is retraced, with the opening of a loop only when 7tii, becomes less than Q_,. If the flows we simulate can actually be realized in experiment, this difference could serve to validate one or the other of our models. A prediction which may be easier to verify, but would not differentiate the models, is the molecular weight-dependence of the width of the hysteresis loop when it rejoins the main branch. According

333

to the adySiS above, this iS at 7imp = Tbot. If we substitute of eqns. (56) into the last of eqns. (53), we find that 7bot= 2C-2M-“‘2&(1

- a2)-1’2.

the power laws (58)

For PI-3, this gives Tbot= 7.4 X lo5 dyne/cm2, which is in accord with Fig. 16 to graphical accuracy. For PI-7, the loop would close at 2.293 x lo5 dyne/cm2. We thus have a critical stress in loading which is independent of molecular weight and a closing of the loop which is strongly dependent on molecular weight. 3.3.3 Experimental verification The experimental verification of our model then rests on the following predictions: (i) No hysteresis or spurt for 7,, I 7t0p, or fluids with M low enough so that E 2 l/9. (ii) Ttclpis virtuaZZy independent of M, as is S in flows at high enough driving gradient, achieved in a loading sequence. @) rbot is associated with the subcritical closing of the hysteresis loop, is strongly dependent on molecular weight and obeys the scaling law rbot

-

0(M-“‘2),

where, for Vinogradov et al.‘s polyisoprene data, (Y= 3.3. (iv) Finite layer thickness. If flows of the type we describe can be sustained without secondary instability for long enough to be observed, we predict that layers having a 6 which is significant with respect to h/2 will be observed, in distinct contrast to true wall slip. But we also predict that 6 is a continuous function of 7imp/7top in supercritical loading. Our spurt model is a continuous one with very Zarge gradients. This is in distinct contrast to the predictions of McLeish and Ball [lo]. (v) Shape memory hysteresis. If spurted flows can be sustained long enough to back the driving gradient off, we predict that no part of the loading curve will be retraced until subcritical closing. This is in distinct contrast to the model of McLeish and Ball. 4. Discussion

We are not proposing that our results constitute a new explanation of wall slip. There surely is ample evidence that breakdown of surface bonding between fluid and die wall give rise to behavior similar to what we describe here [24-261. Such a phenomenon could indeed be responsible for Vinogradov et al.‘s spurt observations. We are suggesting that there may be phenomena which give the appearance of slip but which are governed solely

334 by material properties of fluid. If these phenomena can be observed, they could well be confused with slip, and because of the possibility-or even likelihood-of two and three-dimensional primary or secondary instabilities, the phenomena we describe could be an initiating mechanism in true slip. At least, we believe that our results make a strong argument to the effect that bulk polymer properties cannot be dismissed as major contributing factors in slip situations. Polymer constitutive models with shear stress maxima cannot be dismissed out-of-hand as somehow “unphysical” for at least two reasons: (i) “At the shear stress maximum” shear stress cannot be said to “go down” with steady shear strain rate, since classical steady solutions cease to exist in PPF and cease to be linearly stable even to one-dimensional disturbances in PCF; (ii) The equations retain evolutionarity “beyond the shear stress maximum” so that asymptotic, nonclassical steady solutions exist, and these seem to describe observed physical behavior. Total physical stress never decreases with increasing strain rate in such a solution. Observation (i) points out the value of our stability analysis: While it is not complete enough to establish what solutions will occur, it has a crucial message in establishing what solutions will not occur. As is corroborated by numerical simulation, the decreasing portion of the steady stress-strain rate curve is essentially unattainable. Further study of more complicated classes of instabilities is essential to determine how stable the spurt solutions are and thus how likely it is that the mechanism we have identified here is solely reponsible for observed spurt; this may well be the case, since our model can reproduce the spurt effect with some precision. In summary, we believe that the validity of the shear stress maximum feature of a constitutive model cannot be definitively assessed without consideration of its dynamic and asymptotic implications. Acknowledgements The research of D.S. Malkus on the development of numerical methods for viscoelastic flow is sponsored by AFOSR Grant No. 85-0141. His participation in interdisciplinary research in viscoelasticity and rheology is sponsored by the Center for the Mathematical Sciences University of Wisconsin, Madison, with receives support for research in these areas from AR0 (Grant No. D7AAL03-87-K-0036), NSF (Grant No. DMS-871205Q and AFOSR (Grant No. 87-0191). The authors would like to thank J.A. Nohel, M. Slemrod, and B. Plohr of the Center for Mathematical Sciences for many valuable discussions concerning the spurt instability of the Johnson-Segalman model. They also express their gratitude to T.C.B. McLeish for his cogent explanation of his own ideas about spurt (which predate our own).

335 References 1 R.B. Bird, 0. Hassager and R. Armstrong (Vols. 1 and 2) and CF. Curtiss (Vol. 2) Dynamics of Polymeric Liquids, 2nd ed. Wiley, New York, 1987. 2 H. Giesekus, J. Non-Newtonian Fluid Mech., 11 (1982) 69. 3 P.L. Rimmer, Rheol. Acta, 10 (1971) 601. 4 C.J.S. Petrie and M.M. Denn, AIChE J., 22 (1976) 209. 5 J.P. Tordella, in: Rheology Theory and Applications 5, F. Eirich (Ed.), Academic Press, New York, NY, 1969, p. 57. 6 G.V. Vinogradov, A. Ya. Malkin, Yu. G. Yanovskii, E.K. Borisenkova, B.V. Yarlykov and G.V. Berezhnaya, J. Polym. Sci., A-2, 10 (1972) 1061. 7 L.L. Blyler, jr. and A.C. Hart, jr., Polym. Eng. Sci., 10 (1970) 193. 8 T.U. Huseby, Trans. Sot. Rheol., 10 (1966) 181. 9 Y.-H. Lin, J. Rheol., 29 (1985) 605. 10 T.C.B. McLeish and R.C. Ball, J. Polym. Sci., B, 24 (1986) 1735. 11 J.K. Hunter and M. Slemrod, Phys. Fluids, 26 (1983) 2345. 12 M.W. Johnson and D. Segalman, J. Non-Newtonian Fluid Mech., 2 (1977) 278. 13 J. Yerushalmi, S. Katz, and R. Shinnar, Chem. Eng. Sci., 25 (1970) 1891. 14 M. Renardy, W.J. Hrusa and J.A. Nohel, Mathematical Problems in Viscoelasticity, Surveys in Pure and Applied Mathematics, Vol. 35. Longman Group Limited (formerly Pitman II monographic series), copublished with John Wiley and Sons, New York, 1987. 15 H.B. Squire, Proc. R. Sot. London, A, 142 (1933) 621. 16 D.D. Joseph, M. Renardy and J.C. Saut, Arch. Ration. Mech. Anal., 87 (1985) 213. 17 D.D. Joseph and J.C. Saut, J. Non-Newtonian Fluid Mech., 20 (1986) 117. 18 G. Iooss and D.D. Joseph. Elementary Stability Theory and Bifurcation, Springer, N.Y., 1980. 19 J.J. van Schaftingen and M.J. Crochet, J. Non-Newtonian Fluid Mech., 18 (1985) 335. 20 D.S. Malkus, J.A. Nohel and B. Plohr, Technical Report No. xx, Ctr. for Math. Sci., Univ. of Wis., Madison, in preparation. 21 A.T. Patera and S.A. Orszag, J. Fluid Mech., 112 (1981) 467. 22 T.J.R. Hughes, The Finite Element Method. Prentice-Hall, Englewood Cliffs, NJ, 1987. 23 T.J.R. Hughes, W.K. Liu and A.N. Brooks, J. Comput. Phys., 30 (1979) 1. 24 D.S. Kalika and M.M. Denn, J. Rheol. 31 (1987) 815. 25 W.R. Schowalter, J. Non-Newtonian Fluid Mech., 29 (1988) 25. 26 A.V. Ramamurthy, J. Rheol., 30 (1986) 337.