175
Journal of NonNewtonian Fluid Mechanics, 32 (1989) 175195 Elsevier Science Publishers B.V., Amsterdam  Printed in The Netherlands
NUMERICAL EQUATIONS
SOLUTION OF THE LAMINAR FOR POWERLAW FLUIDS
HELGE I. ANDERSSON
and TERJE H. TOFTEN
BOUNDARY
LAYER
*
Division of Applied Mechanics, Norwegian Institute of Technology N7034 TrondheimNTH (Norway) (Received September 29,1988;
in revised form January 12, 1989)
Summary
The powerlaw 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 powerlaw fluids. The transformed momentum equation is written as a set of first order equations, which is solved by the implicit Keller Box finitedifference scheme for parabolic problems. This method has secondorder 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 nonsimilar 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, N1300 Sandvika, Norway. 03770257/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 nonNewtonian boundary layer flows as compared with their Newtonian counterparts. Nevertheless, it was demonstrated by Schowalter [l] and Acrivos et al. [2] that exact similarity solutions exist for boundary layer flow of powerlaw fluids along semiinfinite wedgeshaped bodies; i.e. for the same class of flows which was investigated by Falkner and Skan [3] in the Newtonian case. Schowalter [l] considered the nonNewtonian 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. [2], while the exact similarity solution for stagnation point flow (m = 1) was obtained by Shah [4]. The similarity solutions for powerlaw 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 [5] and Andersson and Irgens [6], respectively. Most of these similarity solutions were obtained by numerical techniques. Hsu and Cothern [7], 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 FalknerSkan class U  xM do not permit similarity solutions of the momentum equation. Some different solution techniques for nonsimilar boundary layers of powerlaw fluids have therefore been proposed over the years. Bizzell and Slattery [8] used an extended Von &umanPohlhausen integral method approach, which originally was carried over to nonNewtonian flat plate flow by Acrivos et al. [2]. In a later paper Acrivos et al. [9] 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 powerlaw fluids. Recently, the MerkChao series solution technique for Newtonian boundary layer flow was generalized to powerlaw fluids by Lin and Chem [ll] and Kim et al. [12]. More recently, Nakayama and Koyarna [13] presented a general integral solution procedure for laminar boundary layer flow of powerlaw fluids past bodies of arbitrary geometrical configuration. This calculation procedure is an extension of the integral method for Newtonian [14] and nonNewtonian fluids [15,16] to arbitrary geometries. The NakayamaKoyama approach is essentially a generalization of the Von KarmanPohlhausen integral method for Newtonian fluids, which was first extended to include powerlaw fluids
177 by Acrivos et al. [2] and Bizzell and Slattery [8]. In the NakayamaKoyama 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 NakayamaKoyama approach, is attractive for rapid estimation of wallfriction 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 nonNewtonian 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. [2], for instance, revealed that integralmethod 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 thirdorder polynomial was chosen by Acrivos et al. [2], 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 [17] recently demonstrated that the NakayamaKoyama 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 twodimensional boundary layer momentum equation for powerlaw fluids, without other assumptions involved than those inherent in the finitedifference 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 nonsimilarity boundary layer flows.
178 2. Physical model and governing equations 2.1 Basic equations and boundary conditions We consider twodimensional 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 xdirection 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 crossstream 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 zaxis. At the solid surface the usual impermeability and noslip 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 Ostwaldde 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 nl
 CLI aY I by orderofmagnitude arguments. This twoparameter equation is also known as the powerlaw model, the parameters K and n being the consistency coefficient and the powerlaw 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”
y1
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 FalknerSkan 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 powerlaw 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 powerlaw fluids. 3. Numerical formulation and solution procedure It has been demonstrated that the powerlaw 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 wellestablished engineering calculation method presented by Keller and Cebeci [18,19], and described in detail for instance by Cebeci and Bradshaw [20] and Bradshaw et al. [21]. Only an outline of the method is presented below, together with the necessary modifications required for handling nonNewtonian 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 twopoint finite difference method devised by Keller [22]. This method, frequently referred to as the Box scheme, is unconditionally stable and permits arbitrary net spacings. Furthermore, it has secondorder accuracy in all variables, and Newton’s method is used to solve the resulting nonlinear difference equations. First, we write the thirdorder differential equation (11) in the form of a firstorder system, introducing new dependent variables g and q:
f’=g,
(184
g’ = 4,
(18b) 2mnm+l n+l
i”) ( (/I
‘+
)fq+m(lg’)=x(ggqg),
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=Xkr+hk no=O;
k=l,2 j=l,2
nj=nil+dj
,..., ,...,
Kl, J
1,
(2Oa) (2Ob)
and the variables f, g, 4 at points (xk, ~j) are approximated by fjk, gf, 4;. Following Cebeci and Bradshaw [20], eqns. (Isa) and (18b) are drscretized using centereddifference 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  gf1~ idj( 4; + qik_l)p
(2W
[(b)F
(b)?l]/adj
+ a( q;Z;$;
+ aI(fq
&;q;_;)
= R;_&
 a2( g2>;; (214
where
WZ[l
(g’)j;]
WZk
(22)
and a=Xkf
/
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 xlevel 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 [23] is used to solve this system.
183 3.2 Linearization of the diffusion term While the linearization of the laminarflow equations for a Newtonian fluid (b = 1) is straightforward, the nonlinearity introduced via b for turbulent or nonNewtonian 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 bterm, which represents the deviation from the linear Newtonian diffusion. Recently, Cebeci et al. [24] 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 powerlaw index. 3.3 Miscellaneous A computer program based on the Keller Box scheme, as adapted to nonNewtonian boundary layer problems, was written in TurboPascal for IBMcompatible 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 fourthorder polynomial; see, for instance, Cebeci and Bradshaw [20] and White [25]. 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 crossstream direction was either constant or prescribed as a geometric progression. The program was tested on the classical FalknerSkan 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 [25], 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 [25] for the free stream parameter p = 2m/(m + 1) between 0.18 and 10.0. For the particular pvalue 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 FalknerSkan 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 FalknerSkan wedge flows U  x*, is the flow past a semiinfinite 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 [25]. The corresponding powerlaw boundary layer problem was solved by Acrivos et al. [2] for nvalues between 0.05 and 5.0, while a set of modified numerical results obtained by Shah [4] were presented in Table 1 of Ref. [9].
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 [4] than to the results due to Acrivos et al. [2]. 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 nonNewtonian 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 overrelaxation in the Newton process. For the dilatant fluids, on the other hand, a significant amount of underrelaxation was necessary to obtain convergence for the higher nvalues. 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 underrelaxation. The calculation of the pseudoplastic cases required some underrelaxation 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 nonNewtonian flat plate boundary layer flow are displayed in Figs. 13. 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 crossstream variation of the dimensionless velocity gradient f “( 77)is shown in Fig. 3 for some different nvalues. 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 powerlaw index n for m = 0.
187
6.
Fig. 3. Predicted crossstream 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. [2], 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 illposed 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 powerlaw 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 powerlaw 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 [4] and to within 1.2% as compared with the solutions obtained by Hsu
188 and Cothern [7] 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 [25], 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. [26] on a cylinder immersed in a nonNewtonian 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 leastsquare fit to the experimental pressure distribution along the cylinder surface. Shah et al. [26] 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 powerlaw 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 MerkChao series solutions presented recently by Rim et al. [12]. Unfortunately, the factor of 2 in the empirical free stream distribution (27) has been overlooked by some investigators [lo][13], [16]. 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 xlocations 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 o1 2. 3456(0) 7,
I
I
!
z
II
I
r
24
l2

3. 456
0
4
8
12
4
8
I 12
16
20
4
log/ Af”( O)\ Ol2.
4. 536j (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 powerlaw fluids. By taking advantage of the resemblance between the turbulent eddy viscosity and the effective viscosity which accounts for nonNewtonian behaviour, it has been demonstrated that the powerlaw 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 powerlaw fluids, but to any inelastic “generalized Newton fluid” obeying, for instance, the truncated powerlaw of Spriggs, the Eyring model or the Carreau model (Bird et al. [27]). 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 wellknown calculation technique for standard wall boundary layers to include nonNewtonian behaviour, as demonstrated in the present paper, should also be feasible for other kinds of thinshearlayer flows, like, for instance, jets, wakes and thinfilm flows. Finally, it should be emphasized that the aim of the investigation was to verify the applicability of the Keller Box scheme to nonNewtonian boundary layer problems rather than to demonstrate the validity of the powerlaw 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), crossstream 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), powerlaw index, pressure, dimensionless velocity gradient g’ = f “, shorthand notation for righthand 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, crossstream velocity component, crossstream velocity fluctuation, streamwise coordinate, crossstream 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, crossstream 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) 2428. A. Acrivos, M.J. Shah and E.E. Petersen, AIChE J., 6 (1960) 312317. V.M. FaIkner and SW. Skan, PhiIos. Mag. VII, 12 (1931) 865896. M.J. Shah, Ph.D. Thesis, Univ. CaIif., Berkeley, CA, 1961. S.Y. Lee and W.F. Ames, AIChE J., 12 (1966) 700708. H.I. Andersson and F. Irgens, J. NonNewtonian Fluid Mech., 27 (1988) 153172. C.C. Hsu and J.H. Cothem, ASME Paper No. 71FE35, 1971. G.D. BizzeIl and J.C. Slattery, Chem. Eng. Sci., 17 (1962) 777782. A. Acrivos, M.J. Shah and E.E. Petersen, Chem. Eng. Sci., 20 (1965) 101105. R.W. Serth and K. M. Kiser, Chem. Eng. Sci., 22 (1967) 945956. F.N. Lin and S.Y. Chem, Int. J. Heat Mass Transfer, 22 (1979) 13231329. H.W. Kim, D.R. Jeng and K.J. Dewitt, Int. J. Heat Mass Transfer, 26 (1983) 245259. A. Nakayama and H. Koyama, W&meStofftibertrag., 22 (1988) 2936. A. Nakayama, H. Koyama and S. Ohsawa, Int. J. Heat Mass Transfer, 26 (1983) 17211726. A. Nakayama, A.V. Shenoy and H. Koyama, WarmeStofftibertrag., 20 (1986) 219227. A. Nakayama and H. Koyama, Int. J. Heat Fhtid Flow, 7 (1986) 99101. H.I. Andersson, Int. J. Heat Fluid Flow, 9 (1988) 343346. H.B. Keller and T. Cebeci, Lecture Notes in Physics, 8 (1971) 92100. H.B. Keller and T. Cebeci, Am. Inst. Aeronaut. Astronaut. J., 10 (1972) 11931199. T. Cebeci and P. Bradshaw, Momentum Transfer in Boundary Layers. McGrawHill/ 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) 305320. T. Cebeci, K.C. Chang and D.P. Mack, Am. Inst. Aeronaut. Astronaut. J., 22 (1984) 18191821. F.M. White, Viscous Fluid Flow, McGrawHill, New York, NY, 1974. M.J. Shah, E.E. Petersen and A. Acrivos, AIChE J., 8 (1962) 542549. R.B. Bird, R.C. Armstrong and 0. Hassager, Dynamics of Polymeric Liquids, Vol. 1, 2nd ed. John Wiley & Sons, New York, NY, 1987.