The visco-elastic Ekman layer

Journal of Non-Newtonian Fluid Mechanics, 18 (1985) 87-100 Elsevier Science Publishers B.V., Amsterdam - Printed in The Netherlands

THE VISCO-ELASTIC

J.A. HORWITZ,

EKMAN LAYER

M.S. ENGELMAN

Department of Mathematics,

87

and S. ROSENBLAT

IIIinois Institute of Technology, Chicago, IL 40616 (U.S.A.)

(Received January 14, 1985)

Solutions are obtained for the flow in an Ekman layer in a visco-elastic fluid using various constitutive relations. For some models (those which exhibit no shear thinning in simple shear flow) it is found that the solution is identical with the Newtonian solution for all Weissenberg numbers. For other models it is shown that the solution differs from its Newtonian counterpart and may cease to exist when a critical Weissenberg number is exceeded. In these cases, when the solution exists it is found that the visco-elastic Ekman layer is thinner than the Newtonian Ekman layer.

1. Intruduction

The Ekman layer in fluid mechanics is a boundary layer formed by a balance between coriolis force, viscous stress and pressure. The driving mechanism of the flow in the layer is a combination of uniform rotation and uniform shear. Although the original investigations of Ekman [l] were concerned with the effects of wind stress on oceanic currents, it was subsequently realized that there is a wide variety of flow situations, and not merely in a geophysical context, in which the Ekman layer plays a crucial role. One especially important application is to the flow of fluids confined within rotating cylindrical containers. This aspect has been discussed in full by Greenspan [2], where the dominant part played by the Ekman layers at the cylinder ends in determining the interior flow, particularly during spin-up and spin-down, has been analyzed. Although the structure and stability of the Ekman layer have been thoroughly investigated for Newtonian fluids, there do not appear to have been analogous studies for the case of visco-elastic fluids. There are, however, practical situations where the flow field of a v&o-elastic fluid confined 0377-Q257/85/$03.30

0 1985 Elsevier Science Publishers B.V.

88

within a rotating cylinder is of interest as, for example, in the flight of a spinning projectile filled with a liquid containing polymer additives. Configurations such as these make the study of the v&o-elastic Ekman layer worthwhile. The simplest model (the so-called non-divergent Ekman layer) comprises a semi-infinite volume of fluid bounded by a rigid plane boundary, with the whole system rotating about an axis normal to the plane and with the boundary translating uniformly in its own plane. For a Newtonian fluid there is an exact solution to this problem [3]: the velocity vector lies in planes parallel to the boundary and varies only with distance from the boundary. This solution is displayed in eqn. (20) below. In this paper we obtain solutions to the problem of a non-divergent Ekman layer for a visco-elastic fluid using several different constitutive equations: upper convected Maxwell and Jeffreys, generalized Maxwell and Jeffreys, and Phan Thien-Tanner equations. The solutions we seek are the analogues of the Newtonian solution for the non-divergent Ekman layer: a time-independent flow in which the velocity vector lies in planes parallel to the boundary and varies only with distance from the boundary. For the two upper convected models we find that the solution is identical with the Newtonian solution for all values of the Weissenberg number (relaxation time), a result that can be attributed to the fact that upper convected models exhibit no shear thinning in simple flow configurations. In the case of the generalized Maxwell model and the Phan Thien-Tanner model we show that not only does the solution vary quantitatively with Weissenberg number but also that it ceases to exist beyond a critical value of this number. A similar phenomenon appears in the generalized Jeffreys model when the retardation time is sufficiently small, although for larger values of the retardation time the solution exists for all values of the Weissenberg number. The non-existence of solutions in certain parameter ranges relates only to the type of solution stipulated above, and does not preclude the possibility of other types of solutions, non-steady or non-planar. On the other hand it could indicate an inadequacy of the particular constitutive model for flows of this type. This issue remains an open question. 2. Formulation Consider a visco-elastic fluid of constant density p occupying the half-space adjacent to an infinite rigid plane. The whole system rotates with constant angular velocity St about an axis normal to the plane. A Cartesian coordinate system (x*, y*, z*) rotating with angular velocity Q is chosen such that the z*-axis coincides with the axis of rotation and the origin of coordinates lies

89

on the rigid boundary. With respect to this rotating coordinate equations of momentum and mass conservation are De* + 2pSkCv* = - v*p* ’ Dt

+ v*

system the

* 7*

(1)

where v* is the velocity vector, t* is the time, p* is the pressure and I* is the extra-stress tensor. To begin with, we take as the constitutive relation between stress and deformation rate the generalized Jeffreys model [4] 7* +

xy[zB*7*

- +a( 7* . y* + y* * 7*)]

= po [ y* + x; (LB***

- a** * y*)]

(3)

where 9* is the Jaumarm derivative D =-+;(a*. 9* Dt*

- *o*)

and where the deformation-rate are defined by o* = v*l)* + v*v*T,

tensor and the vorticity tensor, respectively,

a* = v*u* - v*&.

(9 Equation (3) contains four parameters: the zero-shear rate viscosity p,,, a relaxation time XT, a retardation time A:, and the dimensionless parameter a, 0
(6)

and v*+Oasz*+cc.

(7)

The preceding system of equations and boundary mensionalized in the following way. We set (x*, y*, z*) = (po/Poy2(x, 0*=1/v

0 9

p* = P@,

y, 4, T* = /bo(

conditions

is non-di-

t* = (/h_dPnQly2tT

pti&2/~o)1’2T

(8)

and also write D = Qi, where i denotes the unit vector along the positive z-axis. Then eqns. (1) and (2) become, in dimensionless notation

90

where Re is the Reynolds number defined by the formula Re =

( pVt/poQ)1’2.

The constitutive relation (3) becomes r+X[97-$2(7*~+~~7)]

=t+rx(gy-a**0)

where

x = hf(pi2v;/poy2,

E =

X:/q.

03)

The parameter ? can be regarded as a Weissenberg number for this problem, and is in fact a dimensionless relaxation time; e is a retardation parameter. The tensors y and w are the dimensionless forms of eqns. (5). The boundary conditions (6) and (7) have the non-dimensional forms v= (0,1,O) on z = 0

(14)

and v--,Oasz-,oo.

(15) We look for a time-independent solution of the governing equations having no velocity component normal to the boundary and where the field quantities depend only on the z-coordinate. In other words, we seek a solution of the form v= (U(Z),U(Z),O),

7=1(z)

(16) which is the generalization of Barcilon’s solution for a Newtonian fluid [3]. In this case it is easy to show that the momentum equations reduce to - 2v = r13,p, 2u = r23,Z

(17) where 7i3 and r23 are the appropriate components of the extra-stress tensor. These equations have to be solved subject to the boundary conditions (14) and (15), with the extra stress being determined from eqn. (12). The solutions are determined separately for the respective constitutive models in the following sections. 3. Upper convected models The upper convected models result from setting a = 1 in eqn. (12). Having done this, we then substitute the representation (16) into eqn. (12), and it is then a matter of straightforward algebra to calculate the stress components that are required in eqn. (17). We find that 713 =

U,,

‘23 =

‘z

(18)

whereupon eqns. (17) become U z,? =

-2v,

v,, = 2u.

(19)

91

These equations are identical with those obtained for a Newtonian fluid, and the solutions subject to the boundary conditions (14) and (15) are 2.4= e-” sin z ,

u=e-‘cosz

(20) Equations (20) form the classical Newtonian solution for the Newtonian Ekman layer, and one sees that the upper convected models have no effect on the velocity field. This result can be understood if one recalls that upper convected models exhibit no shear thinning in u&directional shear flow. The Ekman layer, albeit not uni-directional, is essentially a shear flow driven by boundary motions and not by pressure gradients, and therefore it might be relatively unsurprising that the velocity field in upper convected models is identical with that of a Newtonian fluid. The other extra-stress components determined by eqn. (12) are found to be 711=2x(1-+&

722=2h(l-C)U,2,

r33=o

(21)

and 7*2= 2A(l - ‘)UzU*.

(22) One sees from eqn. (21) that the normaI stress components sir and 722 are proportional to the Weissenberg number, as is the case in uni-directional simple shear. At the boundary z = 0 the velocity is in the y-direction, and therefore one can define primary and secondary normal stress differences respectively. Using eqns. (20) and (21) we find that 733 - ql, 722 - 733 and 722-

2X(1 -c>,

(23) which coincide with well-known formulas for uni-directional simple shear in upper convected models. In the Ekman layer, of course, these quantities vary with distance from the boundary, which is not the case with simple shear. 7331z4

=

lj,

-

7*1jr_l)

=

-2x(1

-c)

4. Generalized Maxwell model We take the constitutive relation (12) with c = 0 and substitute the representation (16). The component equations for the stress tensor are found to be 711-h(l+a)u,7,3=0 712

- gh(1 + u)( u,T23+ Uz713)= 0

T13 -3X(1

+ a)uzT3, +Q(l

- U)(#,T,, + U&) = u,

722_X(l+LI)U,T23=0 T23- $h(l + a)u,7,, + )A(1 - u)( u,7i2 + r&722)= u,

733+ A(1 - U)(Z4,7~3 + U_.T23)=0.

(24)

92

We solve these to obtain 713 =

u,/(l

+ A2MZ),

7z3= Q/(1 + A2M2)

(25)

where M2 = u,z + v,z

(26)

and A2 =X2(1 - Cz’).

(27)

Hence eqns. (17) become $

[ z&,(1 + PM’)]

= -2u,

$

[ u,/(l

+ AW)]

=

22.4

and the boundary conditions are U=o,u=1onz=o

(29a)

and ~+O,u+Oasz-,cc.

(29b) Equations (28) are non-linear: therefore a numerical solution is required. We have used Runge-Kutta procedures coupled with a shooting technique to handle the fact that we have a two-point boundary-value problem. Some care needs be taken to ensure that the far end-point is sufficiently distant that the solution behaves in .a proper asymptotic way. Actually, it is easy to see from eqns. (28) and (29b) that the solution approaches the Newtonian solution for large z; this observation is both helpful from the computational viewpoint, and also indicates that v&o-elastic effects are confined to some region close to the rigid boundary. Before the results of the .computations are presented it is instructive to rewrite eqns. (28) in a somewhat different form. Expanding these equations by differentiating out the left-hand sides we obtain [l-n’(z+U,‘)]~,,-2A2u,uZu,,

2A2u,~,vZ, = -241

+ [l + A’( u,2- u;)] II,, = 2u(l+

The determinant found to be

+ A2M2)* A2M2)2.

(30)

of the coefficients of the second derivatives u, and u, is

A=l--A4M4

(31) from which it can be inferred that bounded solutions may not exist if A crosses through zero for some value of z. In fact, since M2 + 0 as z + cc by virtue of the boundary conditions, a necessary condition for a bounded solution to exist is that A2M2 < 1

(32)

93 for all z. In the Newtonian case we find from eqn. (20) that M2 decreases monotonically with increasing z, and that M’(O), that is, the value of M2 at z = 0, is 2. Assuming tentatively that M2 is a monotonically decreasing function of z in the viscoelastic case, we may replace eqn. (32) by the requirement A2M2(0) K 1.

(33)

Computations show that M’(O) increases as A increases froom the value 2 at A = 0. It follows that a value A* of A is reached, at and beyond which the inequality (32) is violated and no solution having the structure of eqns.

M(O) E-

3.0 & -0.2 E = 0.09

1.0 0

0.1

0.2

0.3

0.4

0.5

0.6

A

Fig. 1. The variation of M(O), the value of (u: + u:)~/’ at the rigid boundary, with generalized Weissenberg number A, for the generalized Jeffreys model. The curve terminates for the upper convected model (E = 0), for generalized Jeffreys models with retardation parameter c < 4, and also for the Phan Thien-Tanner model.

94

(16) exists. It turns out that A.* = 0.40 (34) and consequently a solution of the stipulated form exists only in the range O A* there is no such solution. The numerical integrations show that M2 is indeed a decreasing function of z in the range (35) that the solution exists, so that M2(0) is the maximum value of M2. Figure 1 shows the variation of M(0) with A in the range (35). It should be noted that the slope of the curve increases with increasing A, and becomes infinite when A = A* where the curve terminates. At this point M*(O) = 2.5, obtained from A*M*(O) = 1 (36)

Fig. 2. The velocity profiles for u and u in the upper convected model. The Newtonian the viscoelastic profile at the terminal Weissenberg number profile is indicated by denoted by - - - - - -.

95 The velocity components U(Z) and u(z) do not change much from their Newtonian forms as A increases from 0 to A*. Figure 2 shows the velocity profiles for A = 0 and A = A*; the difference between them, which is slight, is seen to be a displacement of the profiles towards the boundary in the viscoelastic case. In other words, the layer is a little thinner than in the Newtonian case. From eqns (24) we determine the normal stress components to be ~~~=~(1+a)u,Z/(1+A*M*),~~~=X(1+a)u,2/(1+A*M*), 733 =

-X(1 - u)M*/(l

+ A*M*).

(37)

From these we infer that the primary and secondary ferences at the wall are, respectively 722

-~~~=h[(l+a)u;+(l

733 -

-X[(l

711 =

-a)M*]/(l

-u)M*+(l

normal stress dif-

+A*M*)1,_,

+a)&(1

+A2M2)(z_0.

(38)

Note that these quantities are not in a fixed ratio as h varies. 5. Generalized Jeffreys model

We again take the constitutive law (12) with 0 < c c 1. The component equations for the stress tensor are now found to be T*r-X(1 + U)#,T,, = -CX(1 + U)UZ + u)(u,r*, + uz7*3)= -Ch(1 + u)u;u,

712 -3X(1 713 -

3X(1 + U)z4,T3,+ $A(1 - a)( 24,Tii+ UzT**)= u,

T**-X(1 + u)u,7*, = -CX(l + a)$ ‘23 -

$h(l + U)UzT3,+ ih(l

733 +

x(1

-

a)(

u,713 +

- U)(U,T*, + z&T**) = 0,

uzT23) =

Ch(1 - a)( u; + u;).

(39)

From these we obtain ,ri3 = u,(l + rA*A4*)/(1 + AM*),

723= ~(1 + cA*M*)/(l

+ A2M2)

(40)

whereupon eqns. ‘(17) become $[

u,(l + cA*M*)/(l

-& [ u,(l + rA*M*)/(l

+ A*M*)] = -2u, + A*M*)] = 2u

which have to be solved numerically, (29).

(41)

subject to the boundary

conditions

96

Before proceeding with the computations we note that, as in the previous section, there may be parameter ranges for which no solution (of the given type) exists. On performing the differentiations on the left-hand sides of eqns. (41) we find that the determinant of the coefficients of the second derivatives u,, and u,, is A = (1 i tA2M2)[ cA4M4 + (3~ - l)A2M2 + l] .

(42)

By considering the quantity cA4M4 + (3~ - 1)A2M2 + 1 as a quadratic for A2M2 we see that A2M2 _ 1 - 3C+ /(l - 9e)(l - E) 2E

(43)

is a necessary and sufficient condition for A to have a zero at some point. Recalling the condition 0 -Cc < 1 we can infer the following. (i) When E > 6 there is no value of A2M2 for which A vanishes. Therefore a solution exists for all A when E > 4. (ii) When 0 < c < $ the solution ceases to exist if and when A2M2 attains one of the values given by eqn. (43). For small c these values are approximately A2M2 = 1 + 4c,

l/3

(44)

and since the smaller of these will be encountered first, we conclude that a necessary condition for a bounded solution when 0 < z -C$ is approximately A2M2 < 1 + 4~. Assuming again the monotonicity A2M2(0) -C1 + 4~.

(45) of M2 we can replace eqn. (45) by (46)

These inferences are supported by the computations. We find that when 0 c c -C6 the solution ceases to exist beyond a value A* = A*(E), where A* increases with c from 0.400 at E= 0 to 0.479 at e = 4. The numerical integrations confirm that M is a monotonically decreasing function of z, and the variation of M(0) with A in the range 0 G A < A*(c) is illustrated in Fig. 1. The behavior when 0 < E < 6 is very similar to that when c = 0: the solution terminates at A = A*(e) with infinite slope. When c > $ the solution exists for all A. The variation of M(0) with A in this range is shown in Fig. 1 for one particular value (C = 0.2) and on a different scale in Fig. 3. Note that M(0) is not always a monotone function of A: it may attain a maximum value at a relatively small A and then decrease, possibly increasing again at larger A. However, as c approaches 1 these curves of M(0) versus A tend to become monotone. As A -+ 00 we would expect the solutions to approach the solutions of the

97

-I M(O)

45

Fig. 3. The variation of M(O), the value of (u,” + u~)‘/~ at the rigid boundary, with generalized Weissenberg number A for Jeffreys models with retardation parameter E > $.

limiting forms of eqns. (41), namely CU,, = -2V,

CVZZ = 2u

and the computations u = e-I sin 5,

(47)

do confirm this. These solutions are

v = e+ cos 5

(48)

where 5 = Z/C ‘I2 from which the asymptotes for M(0) are easily calculated. What eqns. (48) do indicate, moreover, is that the layer is thinner than in the Newtonian case by a factor &* when A is large. Some typical profiles of velocity are shown in Fig. 4 in the case c = 0.2 and for several values of A. The thinning of the layer with increasing A is apparent. It is worth remarking that the limit A + cc leading to eqns. (47) and (48) can formally be taken for any E > 0. We find, however, that when c < $ the solution (48) is an isolated solution in the sense that it cannot be reached by any continuous path through increasing A. Finally, a comment regarding the cut-off value c = 6. In m&directional shear flow with a generalized Jeffreys model the condition z > $ is necessary in order that the shear stress be a strictly monotone increasing function of the shear rate. This critical value of z appears here in much the same way, except that now it separates regions of differing solution characteristics.

98

0

I

2

3

I

I

I

c

4

5

6

L

Fig. 4. The velocity profiles for u and u in the corotational model with c = 0.2. The ; the case A =1 is denoted by .-.--.; the Newtonian profile (A = 0) is indicated by case A = 50 is denoted by - - - - - -.

6. Phan Thien-Tanner

model

In order to test whether the results obtained above hold more generally we have determined the solution when the constitutive relation is that proposed by Phan Thien and Tanner [5] on the basis of network theory. With the same notation as used previously, this constitutive equation, in dimensionless form, is [l+~(tr7)]7+X[~~-_)a(7.~+~.7)]

=y

where K> 0 is an additional parameter. On substituting

(49) eqn. (16) into eqn.

99 (49) and performing the necessary algebraic manipulations,

we arrive eventually at the following expressions for the relevant stress components u,Q/( Q2 + A2M2),

713 =

(1 + KF)F=

2Xa(u,r,,

+

~23 =

u,Q/( Q2 + A2M2)

UzT23)

(50) 01)

where F = tr7 = 7ri + 722+ 733

(52)

and Q=l+uF

(53)

The other quantities appearing here have the same meanings as before. To determine the solution, we need to integrate eqns. (17) with the stress components 713 and ~23 given by eqns. (50). The latter, however, are only given implicitly because of the coupling with the trace F through eqn. (51); this added complication has to be resolved first. We note two special limiting cases: (i) When a = 0 (the analogue of the co-rotational model), eqn. (51) gives F = 0, and therefore Q = 1. In this case
uz/Q,

T23 =

G/Q

(54)

while eqn. (51) becomes (1+ ~cF)~F= 2XM2

(55)

This equation is a cubic for F, whose solution will lead to explicit forms for ri3 and r23. In particular, when K is small, eqn. (55) gives F = 2XM2,

Q=1+~-2hM’

(56)

whereupon eqns (54) become 713 =

~,/(l

+

ZKXM’),

~23 = u,/(I

+

2~hM’)

(57)

These have the same forms as eqns. (25) with A2 = ~KX. Thus when a = 1 and K is small, the solution ceases to exist beyond .some value of A which depends on K. Beyond these two special limiting cases we can also consider the more general case with 0 -Ca -z 1 and K small. To order K, we find from eqns. (51) and (53) that Q= 1+

K - 2Aa(

upI

•t ~~~23)

Now using eqns. (50) for the stress components

w

we find that eqn. (58)

100

becomes Q = 1+

K -

[2XaQM2/(

Q2

+

n’M’)]

(59)

which can be solved for small K to give Q= 1+

K -

[2A&4”/(1

+ A2M2)]

(60)

Finally, on substituting eqn. (60) into eqns. (50), we obtain explicit forms for the stress components r13 and 723. The solutions of eqns. (17) with these stress components have been computed numerically for several (small) values of K. We have found that in all cases the solution terminates at a finite value of X (or A); thus the Phan Thien-Tanner model behaves just like the generalized Maxwell model and the generalized Jeffreys model with sufficiently small retardation parameter. All these three classes of models have in common two features; in simple shear flow they exhibit shear thinning and the stress is not a strictly monotone function of the rate of shear. It seems reasonable to infer that the presence of these characteristics causes the breakdown of solution. Acknowledgements

The work described in this paper was supported by the US Army Research Office under Contract No. DAAG29-82-K-0061. The authors are greatly indebted to Dr. J.R.A. Pearson; the present paper is an expanded version of an earlier draft regarding which he made a number of valuable suggestions. References 1 W.W. Pkman, On the influence of the earth’s rotation on ocean currents, Ark. Mat. Astron. Fys, 2 (1905) 1. 2 H.P. Greenspan, The Theory of Rotating Fluids, Cambridge Univ. Press, 1968. 3 V. Barcilon, Stability of a non-divergent Ekman layer, Tellus, 17 (1965) 53. 4 C.J.S. Petrie, Elongational Flows, Pitman, London, 1979. 5 N. Phan Thien and R.I. Tanner, A new constitutive equation derived from network theory, J. Non-Newtonian Fluid Mech., 2 (1977) 353.