Numerical solution of the laminar boundary layer equations for power-law fluids

Numerical solution of the laminar boundary layer equations for power-law fluids

175 Journal of Non-Newtonian Fluid Mechanics, 32 (1989) 175-195 Elsevier Science Publishers B.V., Amsterdam - Printed in The Netherlands NUMERICAL E...

175

Journal of Non-Newtonian Fluid Mechanics, 32 (1989) 175-195 Elsevier Science Publishers B.V., Amsterdam - Printed in The Netherlands

NUMERICAL EQUATIONS

SOLUTION OF THE LAMINAR FOR POWER-LAW FLUIDS

and TERJE H. TOFTEN

BOUNDARY

LAYER

*

Division of Applied Mechanics, Norwegian Institute of Technology N-7034 Trondheim-NTH (Norway) (Received September 29,1988;

in revised form January 12, 1989)

Summary

The power-law fluid boundary layer problem is formulated analogous to a turbulent flow problem. The deviation from Newtonian rheology is accounted for by an effective viscosity replacing the so called turbulent eddy viscosity. A numerical solution procedure for turbulent boundary layer problems is adapted to laminar boundary layer flows of power-law fluids. The transformed momentum equation is written as a set of first order equations, which is solved by the implicit Keller Box finite-difference scheme for parabolic problems. This method has second-order accuracy in all variables, and permits arbitrary net spacings. Newton’s method is used to solve the resulting system of nonlinear difference equations. Special attention is devoted to the extra nonlinearity introduced via the effective viscosity. The application of an incomplete linearization scheme, typical for turbulent flow computations, is compared with an alternative full linearization. Numerical results are presented for the similarity boundary layer flow along a flat plate, and for the non-similar boundary layer along a cylinder surface. The results are compared with other available numerical and series expansion solutions. The rate of convergence for the usual incomplete application of Newton’s method is compared with the convergence of the fully linearized scheme. It was observed that the latter led to substantial savings in computational effort, except for some similarity boundary layers of pseudoplastic fluids.

* Present address: Techno Consult A/S, KjDrboveien 23, N-1300 Sandvika, Norway. 0377-0257/89/\$03.50

0 1989 Elsevier Science Publishers B.V.

176 1. Introduction In spite of their relevance within modern manufacturing and processing industries, relatively modest attention has been devoted to non-Newtonian boundary layer flows as compared with their Newtonian counterparts. Nevertheless, it was demonstrated by Schowalter [l] and Acrivos et al.  that exact similarity solutions exist for boundary layer flow of power-law fluids along semi-infinite wedge-shaped bodies; i.e. for the same class of flows which was investigated by Falkner and Skan  in the Newtonian case. Schowalter [l] considered the non-Newtonian boundary layer equations in the absence of body forces, and demonstrated that similarity solutions are possible provided the external velocity U(x) is proportional to xm, where x is the coordinate along the surface of the body. The particular case m = 0, which corresponds to flow past a flat plate, was solved by Acrivos et al. , while the exact similarity solution for stagnation point flow (m = 1) was obtained by Shah . The similarity solutions for power-law fluid flow past a right angle wedge (m = l/3) and along a wedge with opening angle 2~/3 ( m = l/2) were provided by Lee and Ames  and Andersson and Irgens , respectively. Most of these similarity solutions were obtained by numerical techniques. Hsu and Cothern , however, used Meksyn’s method (method of steepest descent) and obtained similarity solutions as series of incomplete gamma functions for wedge angles from 0 to ?r radians. External flows which do not belong to the Falkner-Skan class U - xM do not permit similarity solutions of the momentum equation. Some different solution techniques for non-similar boundary layers of power-law fluids have therefore been proposed over the years. Bizzell and Slattery  used an extended Von &u-man-Pohlhausen integral method approach, which originally was carried over to non-Newtonian flat plate flow by Acrivos et al. . In a later paper Acrivos et al.  proposed an approximate calculation method, based on asymptotic solutions, for boundary layers with external flow of general form, while Serth and Kiser [lo] extended the Giirtler series method to include power-law fluids. Recently, the Merk-Chao series solution technique for Newtonian boundary layer flow was generalized to power-law fluids by Lin and Chem [ll] and Kim et al. . More recently, Nakayama and Koyarna  presented a general integral solution procedure for laminar boundary layer flow of power-law fluids past bodies of arbitrary geometrical configuration. This calculation procedure is an extension of the integral method for Newtonian  and non-Newtonian fluids [15,16] to arbitrary geometries. The Nakayama-Koyama approach is essentially a generalization of the Von Karman-Pohlhausen integral method for Newtonian fluids, which was first extended to include power-law fluids

177 by Acrivos et al.  and Bizzell and Slattery . In the Nakayama-Koyama approach, however, the forced convection heat transfer problem is solved along with the momentum boundary layer problem, in order to obtain the important heat transfer characteristics. Integral methods, like the Nakayama-Koyama approach, is attractive for rapid estimation of wall-friction and heat transfer rates. The heat transer rate between the wall and the fluid, for instance, is of particular interest in many industrial applications. However, to obtain accurate results for the heat flux, accurate knowledge of the wall friction is essential. Unfortunately, in using integral methods one should not expect the same approximation to reality in non-Newtonian boundary layer problems as is obtained in the corresponding Newtonian cases. The pioneering investigation on a momentum boundary layer along a flat plate by Acrivos et al. , for instance, revealed that integral-method predictions become inaccurate for highly dilatant or highly pseudoplastic fluids. This inaccuracy should obviously be attributed to the approximations for the streamwise velocity profile; a third-order polynomial was chosen by Acrivos et al. , while the widely used Pohlhausen polynomial of fourth degree was assumed in Refs. (8) and (13)-(16). Moreover, for free streams of the particular form U - x”‘, Andersson  recently demonstrated that the Nakayama-Koyama approach cannot be employed for highly pseudoplastic substances. Reliable numerical solutions of the complete conservation equations for mass and momentum can be obtained by modern finite element or finite difference techniques. By invoking the boundary layer approximations, however, the elliptic set of governing equations reduces to a parabolic system of partial differential equations. The parabolic nature of the boundary layer equations implies that the flow is not influenced by downstream events and that no boundary conditions are required on the downstream side of the calculation domain. Moreover, the numerical solution techniques for parabolic equations are of the “marching” type, in contrast to the more expensive “iterative” methods necessary for elliptic problems. The objective of the present investigation is therefore to present a general numerical calculation method, which accurately solves the two-dimensional boundary layer momentum equation for power-law fluids, without other assumptions involved than those inherent in the finite-difference approximations to the spatial derivatives. The emphasis will thus be on the numerical solution technique, while the validation of the rheological model and the boundary layer approximations is outside the scope of this paper. The method under consideration should be applicable to rather arbitrary external flow conditions, and it will be tested and compared with previous results for similarity as well as non-similarity boundary layer flows.

178 2. Physical model and governing equations 2.1 Basic equations and boundary conditions We consider two-dimensional steady flow of a viscous fluid. It is assumed that the laminar flow in the thin boundary layer between an inviscid free stream and a solid surface is governed by the boundary layer equations

au

au

uax+U;js.=u~+---.

dU

1 a7xr

(2)

The x-direction is along the surface and y normal to it. Here, eqn. (1) is the continuity equation, and eqn. (2) represents the simplified momentum equation in the streamwise direction. While u and u denote boundary layer velocity components in the streamwise and cross-stream directions, respectively, U(x) denotes the streamwise velocity component found from the inviscid equation of motion at y = 0, i.e. dU = ----1 dp TiY P dx

dz gYG*

(3)

The body force per unit mass is denoted by g and acts in the opposite direction of the z-axis. At the solid surface the usual impermeability and no-slip conditions u(x, 0) = 0,

(4a)

u(x, 0) = 0

(4b)

are applied. Outside the viscous boundary layer the streamwise velocity component u within the boundary layer should approach the external velocity U(x), while the tangential stress should vanish, i.e. u(x, Y) + U(x)

as

Y-,=)7

(54

7xy + 0

as

y+oO.

(5b)

In practical boundary layer calculations, however, the two conditions (4a) and (5a) are imposed on u and eqn. (4b) on u. Normally, the resulting numerical solution automatically satisfies the additional requirement (5b).

2.2 Rheological model In the present paper we consider only time independent inelastic nonNewtonian fluids which obey the Ostwald-de Waele rheological model. In

179 thin boundary layers the model for the shear stress component of the stress tensor reduces to

“-l au

au

ay=~“Bay,

(64

where B _ K

au n-l

- CLI aY I by order-of-magnitude arguments. This two-parameter equation is also known as the power-law model, the parameters K and n being the consistency coefficient and the power-law index, respectively. The function B, which depends on the shear rate, represents an effective dimensionless viscosity. 2.3 Transformed problem Boundary layer flows have a predominant flow direction and the boundary layer thickness 6 is small compared to a typical length in the main flow direction. Since the boundary layer thickness usually increases with increasing downstream distance, the basic equations are transformed in order to make the transformed boundary layer thickness a slowly varying function of X.

With this objective the partial differential equation (2) is transformed by means of the dimensionless variables \$d= Ux(Re,)-l’(R+l)f(x,

r)

17),

= z (Re,)l’(“+l),

where \Elis the physical stream function, defined such that u = a+/ay,

u=

-a+/ax,

(9)

and thereby automatically satisfying the continuity equation (1). The local Reynolds number Re, is defined as

where a is a dimensionless positive constant.

180

The mathematical problem defined in eqns. (2), (4) and (5) then transforms into ( 2~y--+l)ff~~+m[l-(f~y]

[\$f”]‘+

(11)

=x(f?&fq,

f ‘(x, 0) = 0,

(124

f (x, 0) = 0,

(12b)

f’(x,

d+l

[f”(x,

17)ln+0

\$+a,

as as

(134 9+

00,

(13’4

where f ’ = d f/dq has the particular physical significance of the normalized velocity (14)

f ‘(17) = u/u as is readily derived from eqns. (7)-(9). equation (10) the dimensionless variables

In the transformed

momentum

x dU

“=Edx

05)

and b= 1f”

y-1

06)

represent a driving parameter and a new effective viscosity, respectively. The dimensionless driving parameter m(x) depends only on the external flow velocity U(x) in eqn. (3). Of particular importance is the class of external flows which makes m independent of x, i.e. the Falkner-Skan flows U - xm. Then, with m being constant, the partial differential equation (11) in x and 17 reduces to an ordinary differential equation in n. 2.4 Resemblance with turbulent boundary layers The flow in a turbulent boundary layer is governed by a Reynolds averaged momentum equation of the same form as eqn. (2). In the turbulent case, however, the velocity components u and u denote mean velocities, and the shear stress includes a Reynolds stress term, i.e.

where B=l++. The overbar denotes averaging over fluctuating velocity components V’.

(17b) u’ and

181 The unknown Reynolds shear stress in eqn. (17a) is normally modelled in terms of an eddy or turbulent viscosity c, in accordance with Boussinesq’s eddy viscosity hypothesis. The introduction of the dimensionless effective viscosity B demonstrates the close resemblance between turbulent boundary layers and lam&r boundary layer flow of power-law fluids. It seems therefore likely that numerical solution techniques designated for turbulent boundary layer flow can be carried over to the numerical analysis of laminar boundary layers of power-law fluids. 3. Numerical formulation and solution procedure It has been demonstrated that the power-law boundary layer and the turbulent boundary layer are governed by the same streamwise momentum equation (2). The only distinction is the effective viscosity B, which accounts either for deviation from Newtonian behaviour or for turbulence. However, both effects increase the mathematical complexity of the governing equation by introducing an extra nonlinearity in the viscous term. Since numerous numerical analyses of turbulent boundary layer flows have been published over the years, we take advantage of the accumulated experiences in that field and use the well-established engineering calculation method presented by Keller and Cebeci [18,19], and described in detail for instance by Cebeci and Bradshaw  and Bradshaw et al. . Only an outline of the method is presented below, together with the necessary modifications required for handling non-Newtonian boundary layers. 3.1 The Keller Box scheme In order to solve the nonlinear partial differential equation (11) subject to the boundary conditions (12) and (13a), we use the implicit two-point finite difference method devised by Keller . This method, frequently referred to as the Box scheme, is unconditionally stable and permits arbitrary net spacings. Furthermore, it has second-order accuracy in all variables, and Newton’s method is used to solve the resulting nonlinear difference equations. First, we write the third-order differential equation (11) in the form of a first-order system, introducing new dependent variables g and q:

f’=g,

(184

g’ = 4,

(18b) 2mn-m+l n+l

i”) ( (/I

‘+

)fq+m(l--g’)=x(gg--qg),

084

182 with boundary conditions f(x,

0) = g(x, 0) = 0,

@a>

g(x,

77,) = 1.

(19b)

Next, finite difference approximations are introduced on discrete net points x0=0. 9 Xk=Xk-r+hk no=O;

k=l,2 j=l,2

nj=ni-l+dj

,..., ,...,

K-l, J-

1,

(2Oa) (2Ob)

and the variables f, g, 4 at points (xk, ~j) are approximated by fjk, gf, 4;. Following Cebeci and Bradshaw , eqns. (Isa) and (18b) are drscretized using centered-difference derivatives and averages about the midpoint (x“, nj_r&, whereas in eqn. (18~) the derivatives are centered about the midpoint (~~-l’~, ~j_l/2). The finite difference approximations to the differential momentum equations (18) then become

r;“-hk1=fdj(gT+gf_l)y

(214

gf - gf-1~ idj( 4; + qik_l)p

(2W

[(b)F-

+ a( q;Z;\$;

+ aI(fq

-&-;q;_;)

= R;_&

- a2( g2>;-; (214

where

-WZ[l-

(g’)j-;]

-WZk

(22)

and a=Xk-f

/

hk

3

(234

q = (2mn - m + l)/( n + 1) + (Y,

(23b)

(Y2= Wzk+ 1y.

(234

Note that Rk-’ defined in eqn. (22) involves only known quantities at x-level xkP1, except the driving parameter mk which is assumed known from the solution of the inviscid flow problem (3). The nonlinear system of algebraic equations (21) can now be linearized using Newton’s method, i.e. considering iterates of the type fci+‘) =fci) + Af(‘) for the unknowns hk, g; and qT_ The resulting linearized system of equations exhibits a block tridiagonal structure, the blocks consisting of 3 x 3 matrices, and the efficient block tridiagonal factorization scheme due to Keller  is used to solve this system.

183 3.2 Linearization of the diffusion term While the linearization of the laminar-flow equations for a Newtonian fluid (b = 1) is straightforward, the nonlinearity introduced via b for turbulent or non-Newtonian flows presents some difficulties. For turbulent boundary layers [20,21], Newton’s method is normally applied to all terms except for the turbulent diffusion terms where the values of b are assumed known from the previous iteration, i.e.

=b(i) ( qw + Aq”’ @4Yi+1)

),

With this incomplete linearization the difference equations are solved with prior knowledge of the important b-term, which represents the deviation from the linear Newtonian diffusion. Recently, Cebeci et al.  demonstrated that the consequences of this incomplete application of Newton’s method are that the solutions may oscillate and the convergence rate is slower than it need be. Accordingly, they applied Newton’s method to all terms in the momentum equation, and observed that for some typical turbulent boundary layer calculations the full application of Newton’s method provided a quadratic rate of convergence. In the context of the present investigation, the fulZ linearization of the discretized momentum equation (21) requires that also the nonlinear diffusion terms are linearized. With the effective viscosity b defined as in eqn. (16), we obtain (

bq)ci+l)

=

(

q”)(i+l)

~

(q”)(ojl

+ n\$J)

= b(O( q(i) + nL\q(i)),

(25)

where n is the power-law index. 3.3 Miscellaneous A computer program based on the Keller Box scheme, as adapted to non-Newtonian boundary layer problems, was written in Turbo-Pascal for IBM-compatible Personal Computers. The numerical scheme, which takes advantage of the parabolic nature of the problem and marches the solution in the streamwise direction, requires an initial guess for the independent variables f, g and q at the stagnation point x = 0. Here, these initial profiles were derived from the famous Pohlhausen fourth-order polynomial; see, for instance, Cebeci and Bradshaw  and White . The program may use either incomplete or full linearization of the discretized momentum equation (21). In order to improve the convergence of the Newton process, the program also allows for some over or under

184 relaxation of the Newton iterations. In the present investigation we used a constant grid spacing in the streamwise direction, while the spacing in the cross-stream direction was either constant or prescribed as a geometric progression. The program was tested on the classical Falkner-Skan problems for Newtonian boundary layers; i.e. n = 1 and m = constant. In order to compare the present results with the similarity solutions given by White , the arbitrary constant a was taken as 2/( m + 1). The predicted results for the dimensionless wall shear stress f”(0) agreed to within 0.02% with the values reported by White  for the free stream parameter p = 2m/(m + 1) between -0.18 and 10.0. For the particular p-value -0.19884 we obtained the slightly negative value f”(0) = -0.00264, while White got f”(0) = 0 which corresponds to flow separation. 4. Numerical results and discussions The numerical solution technique described above is capable of solving the momentum boundary layer problem for rather arbitrary free stream distributions. However, if the external flow varies according to U - xm the problem defined in eqns. (ll)-(13) reduces to a similarity problem and its solution f( 11) applies at any streamwise location. This important property of the mathematical problem is also reflected by the numerical scheme, which exactly reproduces the numerical solution at the stagnation point x = 0 at any downstream position. If the external flow does not belong to the Falkner-Skan class U - xm, on the other hand, the solution f(x, 17) is nonsimilar and must be calculated at every streamwise position x. Due to this important distinction numerical results for similarity and nonsimilarity boundary layer problems are presented separately below. 4.1 Similarity boundary layers The very simplest external flow, which belongs to the Falkner-Skan wedge flows U - x*, is the flow past a semi-infinite flat plate. For flat plate flow the free stream velocity U is independent of the streamwise position x and m = 0. In this particular case the third term on the left hand side of the momentum equation (11) vanishes. In Newtonian boundary layer theory this case is known as the Blasius problem, see, e.g. White . The corresponding power-law boundary layer problem was solved by Acrivos et al.  for n-values between 0.05 and 5.0, while a set of modified numerical results obtained by Shah  were presented in Table 1 of Ref. .

185 TABLE 1 Dimensionless wall velocity gradient f”(0) n

Ref. (2)

Present

for flat plate boundary layer (m = 0)

Deviation

Ref. (4)

Present

(W 0.05 0.1 0.2 0.3 0.5 1.0 1.5 2.0 2.5 3.0 4.0 5.0

1.40096 0.72986 0.50562 0.35429 0.33120 0.33206 0.36322 0.40150 0.43191 0.45956 0.51022 0.55170

1.49768 0.80554 0.48224 0.38808 0.33124 0.33206 0.36477 0.39979 0.43265 0.46248 0.51365 0.55550

6.90 10.37 - 4.62 9.54 0.01 0.00 0.42 - 0.43 0.17 0.64 0.67 0.69

Deviation (a)

1.5691 0.8787 0.5725 0.4749 0.4341 0.4696 0.5258 0.5766 0.6188 0.6540 0.7036 0.7741

1.5685 0.8786 0.5616 0.4748 0.4340 0.4696 0.5263 0.5766 0.6188 0.6540 0.7086 0.7487

-

0.04 0.01 1.90 0.02 0.02 0.00 0.10 0.00 0.00 0.00 0.71 - 3.28

In order to make comparisons with the results from Refs. (2) and (4), the arbitrary constant a, which appears in the transformation (7), (8) via Re, defined in eqn. (lo), was set equal to 1 and n + 1, respectively. Table 1 shows that the results obtained for the dimensionless wall velocity gradient f”(0) with the present numerical technique are closer to the solution of Shah  than to the results due to Acrivos et al. . The systematic deviation between the two sets of results are due to the different choice of the parameter a. The way of handling the nonlinearity in the viscous term is of particular importance in non-Newtonian fluid mechanics. It was observed that the numerical calculations converged to the same numerical solutions, irrespective of the linearization scheme employed. When the incomplete linearization was used, the calculations converged for all values of n considered. Converged solutions for the pseudoplastics could be obtained without any relaxation, but the convergence rate was significantly improved by using some over-relaxation in the Newton process. For the dilatant fluids, on the other hand, a significant amount of under-relaxation was necessary to obtain convergence for the higher n-values. The opposite tendency was observed when the problem was fu& Zinearized. Convergence could be achieved for any of the dilatant fluids considered without relaxation, but the convergence rate was significantly improved by using a small amount of under-relaxation. The calculation of the pseudoplastic cases required some under-relaxation to converge, and did not converge at all for the highly pseudoplastic fluids (n < 0.3).

186 1.0 q)

0.9 0.0 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0.0 0.0

1.0

2.0

3.0

4.0

.O

5.0 1

Fig. 1. Similarity velocity profiles for flat plate boundary layers (m = 0).

Some of the numerical results for the non-Newtonian flat plate boundary layer flow are displayed in Figs. 1-3. These results have been calculated with a = n + 1. Figure 1 shows the dimensionless velocity profiles f’(q) for some different values of the power law index n. It is observed that the form of the velocity profiles changes dramatically as n is varied. The slope of the profiles at the wall is strongly dependent on n, and this dependency is further demonstrated by the f “(0) data plotted in Fig. 2. The variation of f “(0) with n is, however, partly due to the choice made for the constant a. The cross-stream variation of the dimensionless velocity gradient f “( 77)is shown in Fig. 3 for some different n-values. Most of the solutions are observed to decrease monotonically from f “(0) at the wall to zero outside the viscous boundary layer. For the two highly dilatant cases n = 3 and

1.6 1.5 1.4 1.3 1.2 1.1 1.0 0.9 0.8 0.7 0.6 0.5 -1

0.4 0.0

1.0

2.0

3.0

4.0

5.0

6.0 n

Fig. 2. Predicted variation of dimensionless wall gradient f “(0) with power-law index n for m = 0.

187

6.

Fig. 3. Predicted cross-stream variation of the velocity gradient f”(q). layer m = 0.

Flat plate boundary

n = 5, however, oscillating solutions are obtained in the outer part of the boundary layer. This is obviously due to the observation by Acrivos et al. , that the boundary condition (13a) can only be imposed on the flat plate (i.e. m = 0) momentum equation if n < 2. For more dilatant substances they alternatively imposed both conditions (13a) and (13b) at a finite distance from the wall rather than at infinity. They therefore circumvented this problem by changing the boundary conditions, while the oscillating solutions displayed in Fig. 3 were obtained by imposing the standard outer condition (13a) only. The oscillations thus reflect that the mathematical problem defined in eqns. (ll)-(13a) is ill-posed for n > 2. However, the close correspondence between the present wall gradient data and the other results for n = 3 and n = 5 in Table 1, indicates that the flow in the near wall region is more or less unaffected by the flow in the outer part of the boundary layer. In this respect, it should be mentioned that the power-law model (6a) is valid only in regions with a certain minimum shear rate. If the velocity gradient tends to zero, for instance, the effective viscosity B in eqn. (6b) becomes infinite for n < 1, which is obviously unrealistic. The applicability of the power-law model in the outer part of a viscous boundary layer may therefore be questioned. The present numerical approach has also been used to calculate the wall velocity gradient f “(0) for stagnation point flow (m = 1) and for boundary layer flow past a wedge with opening angle 2~/3 (m = l/2). It was observed that the predictions for m = 1 agree to within 0.5% with those of Shah  and to within 1.2% as compared with the solutions obtained by Hsu

188 and Cothern  using Meksyn’s technique. For m = l/2 the results differed from those in Ref. (6) by not more than 0.25%. 4.2 Nonsimilar boundary layers An interesting and important example of nonsimilar boundary layers is the flow past a circular cylinder. From potential flow theory, see, for instance, White , the inviscid velocity distribution becomes U(X)/&

= 2 - sin(x/R),

(26)

where U, denotes the uniform stream velocity approaching the cylinder, x denotes the streamwise coordinate measured along the cylinder surface from the forward stagnation point, and R is the cylinder radius. The experimental investigation by Shah et al.  on a cylinder immersed in a non-Newtonian crossflow, revealed that the theoretical velocity distribution (26) represents a reasonable approximation to reality in the region from 0 to 60 degrees measured from the forward stagnation point, i.e. up to x/R = 7r/3 = 1.05. The alternative empirical external flow distribution, U(x)/&

= 2 - [0.92(x/R)

- O.l31(~,‘R)~],

(27)

was obtained by a least-square fit to the experimental pressure distribution along the cylinder surface. Shah et al.  also observed that the pressure distributions for water and the CMC 70 (carboxymethylcellulose) solution (0.63 < n < 0.79) had negligible difference. The development of a power-law viscous boundary layer along the cylinder surface has been predicted by the Keller Box scheme. The numerical results are presented in terms of the dimensionless friction parameter

(28) where

(29) Cf =

27,/m&

Re = Uz-“R”/(

(30) K/p),

(31)

are the wall shear stress, the dimensionless skin friction coefficient and the characteristic Reynolds number, respectively. The predicted variation of the wall friction along the cylinder surface is shown in Fig. 4. It is observed that maximum friction occurs at about x/R = 1 for all fluids considered, and the

189

&pe&l 2.4 2.2 2.0 1.8 1.6 1.4 1.2 1.0 0.8 0.6 0.4 0.2

-I

0.0 0.4

0.2

0.0

1.0

0.8

0.6

1.2

1.4

1.6

1.8 x

R

Fig. 4. Predicted variation of wall friction along cylinder surface. Theoretical free stream.

friction decreases continuously farther downstream until the separation point is reached at about 1.8 < x/R < 1.9. The wall friction in the stagnation point region is higher and the maximum friction is lower for the pseudoplastic fluid (n = 0.52) than in the Newtonian case. The opposite tendency is observed in the dilatant case (n = 1.6). The results in Fig. 4 were calculated using the theoretical free stream distribution (26). The data in Fig. 5 for a pseudoplastic boundary layer demonstrate that the main effect of using the alternative empirical external flow (27) is to reduce the wall friction along the cylinder surface and thereby promote the tendency of the flow to separate from the surface. The results

q

0.0,

0.0

Kim,

I 0.2

0.4

Jeng

81 Dewitt

I

I I

0.6

0.8

(1983)

I 1.0

1.2

I 1.4

1.6

I.0 li.

R

Fig. 5. Predicted variation of wall friction along cylinder surface for a pseudoplastic fluid (n = 0.52). Theoretical (-) and experimental (- - - - - -) free stream.

Fig. 6. Predicted variation of wall friction along cylinder surface for a dilatant fluid (n = 1.6). Theoretical (-) and experimental (- - - - - -) free stream.

are also compared with the Merk-Chao series solutions presented recently by Rim et al. . Unfortunately, the factor of 2 in the empirical free stream distribution (27) has been overlooked by some investigators [lo]-, . The tabulated results in Ref. (12) have therefore been resealed before they could be used for comparisons in Figs. 5 and 6. The effect on the boundary layer development induced by the two different free stream distributions, as displayed in Fig. 5, was also observed in the dilatant flow case (n = 1.6) shown in Fig. 6. In contrast with the similarity solutions presented in Section 4.1, the nonsimilar boundary layer problem considered in this Section required the iterative Newton process to be carried out at each streamwise location, the typical streamwise grid spacing being h = 0.05 R. At the upstream position x = 0 the flow starts as a stagnation point flow, and the transformed momentum equation (11) reduces to a similarity problem. The solution is then marched in the streamwise direction, and the initial data for the Newton process at the following x-locations are taken from the solution at nearest upstream station. While the convergence rate at x = 0, see Fig. 7, is as for the similarity problem discussed in Section 4.1, the solutions at the following downstream locations (Fig. 8) converge more rapidly. The convergence rates at x/R = 0.05 shown in Fig. 8 demonstrate that the convergence of the incomplete linearization scheme is crucially dependent upon the fluid properties, while the more efficient full scheme is less dependent upon n. In the latter case the error term in the wall gradient f”(0) is reduced from about 10e3 to below 10e6 with 2 iterations for all cases considered. Using the standard incomplete linearization, on the other

191 log/ Af”( O)l o-1 -2. -3-4-5-6(0) -7,

I

I

!

z

II

I

r

24

-l-2

-

-3. -4-5-6-

0

4

8

12

4

8

I 12

16

20

4

log/ Af”( O)\ O-l-2.

-4. -5-3-6-j (Cl -7,

0

I

NUMBER

I 16

II

20

1

24

OF ITERATIONS

Fig. 7. Convergence rate at x = 0. q full linearization, + incomplete linearization: pseudoplastic fluid n = 0.52, (b) dilatant fluid n = 1.6, and (c) Newtonian fluid n = 1.0.

(a)

hand, 8 iterations are required to reduce the error term by the same factor in the pseudoplastic case, and even more iterations are necessary for the dilatant boundary layer. 5. Concluding remarks The Keller Box scheme, which is frequently employed for turbulent boundary layer calculations, has successfully been adapted to laminar

192

logI

-8

(b)

0

I

I

I

2

4

(C) -8, I.8 2 0

I 4

I

I

I

6

I

8

I

I

I

10

12

-3 log/ Af”(

0)I -4-

-5-

-6.

-7-

I

Ia 6 NUMBER

I I I 8 ITERATION:’ 10 OF

Fig. 8. Convergence rate at x = 0.05. 0 full linearization, + incomplete linearization: (a) pseudoplastic fluid n = 0.52, (b) dilatant fluid n = 1.6, and (c) Newtonian fluid n = 1.0.

boundary layer flow of power-law fluids. By taking advantage of the resemblance between the turbulent eddy viscosity and the effective viscosity which accounts for non-Newtonian behaviour, it has been demonstrated that the power-law rheological model can be introduced via the effective viscosity concept. The extra nonlinearity introduced in the momentum equation has been treated by two different linearization schemes: the standard incomplete

193 linearization and an alternative full scheme. It has been observed that the latter led to substantial savings in computational effort, except for some similarity boundary layers of pseudoplastic fluids. It is conjectured that the presented method is applicable not only to power-law fluids, but to any inelastic “generalized Newton fluid” obeying, for instance, the truncated power-law of Spriggs, the Eyring model or the Carreau model (Bird et al. ). While the application of the incomplete linearization scheme to other fluids is straightforward, the full linearization certainly becomes somewhat more involved. The adaptation of a well-known calculation technique for standard wall boundary layers to include non-Newtonian behaviour, as demonstrated in the present paper, should also be feasible for other kinds of thin-shear-layer flows, like, for instance, jets, wakes and thin-film flows. Finally, it should be emphasized that the aim of the investigation was to verify the applicability of the Keller Box scheme to non-Newtonian boundary layer problems rather than to demonstrate the validity of the power-law model and the boundary layer approximations. The lack of experimental data for comparisons does therefore not weaken the conclusions arrived at. Nomenclature

riz B

cf

r” g h J K m n

P i

dimensionless constant, effective dimensionless viscosity, effective dimensionless viscosity (6b) and (17b), skin friction coefficient (30), cross-stream grid spacing (20b), dimensionless stream function (7), dimensionless velocity f ‘, also: gravitational acceleration, streamwise grid spacing (2Oa), number of net points across boundary layer, number of net points in the streamwise direction, also: coefficient of consistency, dimensionless driving parameter (15), power-law index, pressure, dimensionless velocity gradient g’ = f “, short-hand notation for right-hand side of eqn. (21c), also: cylinder radius,

194 Re

Rex U U'

u v V’ X

Y Z

characteristic Reynolds number (31), local Reynolds number (lo), streamwise velocity component, streamwise velocity fluctuation, external velocity, cross-stream velocity component, cross-stream velocity fluctuation, streamwise coordinate, cross-stream coordinate, vertical coordinate.

Greek symbols dimensionless parameter (23a), dimensionless parameter (23b), dimensionless parameter (23c), free stream parameter 2m/(m + l), boundary layer thickness, kinematic eddy viscosity, dimensionless variable (8), dynamic molecular viscosity, kinematic molecular viscosity, density, viscous stress, physical stream function (9). Subscripts

f

j W X

Y Cc

friction sequence number (2Ob), wall condition, streamwise direction, cross-stream direction, asymptotic value outside boundary layer, also: uniform stream velocity.

Superscripts

(i> k

iteration number in Newton process, sequence number (20a).

195 References 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27

W.R. Schowaher, AIChE J., 6 (1960) 24-28. A. Acrivos, M.J. Shah and E.E. Petersen, AIChE J., 6 (1960) 312-317. V.M. FaIkner and SW. Skan, PhiIos. Mag. VII, 12 (1931) 865-896. M.J. Shah, Ph.D. Thesis, Univ. CaIif., Berkeley, CA, 1961. S.Y. Lee and W.F. Ames, AIChE J., 12 (1966) 700-708. H.I. Andersson and F. Irgens, J. Non-Newtonian Fluid Mech., 27 (1988) 153-172. C.C. Hsu and J.H. Cothem, ASME Paper No. 71-FE-35, 1971. G.D. BizzeIl and J.C. Slattery, Chem. Eng. Sci., 17 (1962) 777-782. A. Acrivos, M.J. Shah and E.E. Petersen, Chem. Eng. Sci., 20 (1965) 101-105. R.W. Serth and K. M. Kiser, Chem. Eng. Sci., 22 (1967) 945-956. F.N. Lin and S.Y. Chem, Int. J. Heat Mass Transfer, 22 (1979) 1323-1329. H.W. Kim, D.R. Jeng and K.J. Dewitt, Int. J. Heat Mass Transfer, 26 (1983) 245-259. A. Nakayama and H. Koyama, W&me-Stofftibertrag., 22 (1988) 29-36. A. Nakayama, H. Koyama and S. Ohsawa, Int. J. Heat Mass Transfer, 26 (1983) 1721-1726. A. Nakayama, A.V. Shenoy and H. Koyama, Warme-Stofftibertrag., 20 (1986) 219-227. A. Nakayama and H. Koyama, Int. J. Heat Fhtid Flow, 7 (1986) 99-101. H.I. Andersson, Int. J. Heat Fluid Flow, 9 (1988) 343-346. H.B. Keller and T. Cebeci, Lecture Notes in Physics, 8 (1971) 92-100. H.B. Keller and T. Cebeci, Am. Inst. Aeronaut. Astronaut. J., 10 (1972) 1193-1199. T. Cebeci and P. Bradshaw, Momentum Transfer in Boundary Layers. McGraw-Hill/ Hemisphere, Washington, DC, 1977. P. Bradshaw, T. Cebeci and J.H. Whitelaw, Engineering Calculation Methods for Turbulent Flow. Academic Press, London, 1981. H.B. Keller, in: J. Bramble (Ed.), Numerical Solutions of Partial Differential Equations, Vol. II. Academic Press, New York, NY, 1970. H.B. Keller, SIAM J. Num. Anal., 11 (1974) 305-320. T. Cebeci, K.C. Chang and D.P. Mack, Am. Inst. Aeronaut. Astronaut. J., 22 (1984) 1819-1821. F.M. White, Viscous Fluid Flow, McGraw-Hill, New York, NY, 1974. M.J. Shah, E.E. Petersen and A. Acrivos, AIChE J., 8 (1962) 542-549. R.B. Bird, R.C. Armstrong and 0. Hassager, Dynamics of Polymeric Liquids, Vol. 1, 2nd ed. John Wiley & Sons, New York, NY, 1987.