TWO

EXACT SOLUTIONS
OF THE
NAVIER-STOKES EQUATIONS

21 Introduction
Because of the great complexity of the full compressible Navier-Stokes equations,
no known general analytical solution exists. Hence, it is necessary to simplify
the equations either by making assumptions about the uid, about the ow
or about the geometry of the problem in order to obtain analytical solutions.
Typical assumptions are that the ow is laminar, steady, two-dimensional, the
uid incompressible with constant properties and the ow is between parallel
plates. By so doing it is possible to obtain analytical, exact and approximate
solutions to the Navier-Stokes equations.
Before proceeding let us clearly dene what is meant by analytical, exact and
approximate solutions. An analytical solution is obtained when the governing
boundary value problem is integrated using the methods of classical dierential
equations. The result is an algebraic expression giving the dependent variable(s)
as a function(s) of the independent variable(s). An exact solution is obtained
by integrating the governing boundary value problem numerically. The result
is a tabulation of the dependent variable(s) as a function(s) of the independent
variables(s). An approximate solution results when methods such as series expansion and the von Karman-Pohlhausen technique are used to solve the governing
boundary value problem (see Schlichting [Schl60], p. 239).

22 Analytical Solutions
Finding analytical solutions of the Navier-Stokes equations, even in the uncoupled case (see Section 110), presents almost insurmountable mathematical diculties due to the nonlinear character of the equations. However, it is possible to
nd analytical solutions in certain particular cases, generally when the nonlinear convective terms vanish naturally. Parallel ows, in which only one velocity
component is dierent from zero, of a two-dimensional, incompressible uid have
this characteristic. Examples for which analytical solutions exist are parallel ow

24

Laminar Flow Analysis

through a straight channel, Couette ow and Hagen-Poiseuille ow, i.e., ow in

a cylindrical pipe. Here we discuss parallel ow through a straight channel and
Couette ow.

23 Parallel Flow Through A Straight Channel

A ow is considered parallel if only one component of the velocity is dierent from
zero. In order to illustrate this concept, consider two-dimensional steady ow in a
channel with straight parallel sides (see Figure 21). This ow is two-dimensional
since the velocity is in the x direction and, as we shall see, its variation is in
the y direction. We consider the uid to be incompressible and to have constant
properties. Under these circumstances the momentum and energy equations are
uncoupled (see Section 110). Thus, we consider only the continuity and the
momentum equations, i.e., the velocity eld, and reserve our discussion of the
energy equation, i.e., the temperature eld, until Chapter 4. The continuity
and momentum equations for two-dimensional, steady, incompressible, constant
property ow are
continuity
v
u
+
=0
x
y

(1 49)

x momentum
u

u
2u
u
1 P
2u
+
+v
=
+
2
x
y
x
x
y2

(1 50a)

v
2v
v
1 P
2v
+
+v
=
+
x
y
y
x2
y2

(1 50b)

y momentum

Since the ow is constrained by the at parallel walls of the channel, no component of the velocity in the y direction is possible, i.e., v = 0. This implies

u( y)

y
C
L

2h
x

umax

Figure 21. Parallel ow through a straight channel.

Parallel Flow Through A Straight Channel

25

that v = 0 everywhere. Hence, the gradients in v are also equal to zero, i.e.,
v/y = v/x = 2 v/y 2 = 2 v/x 2 = 0.
From the continuity equation we have
u
v
=
=0
x
y
which implies that u2 /x2 = 0. The momentum equations are thus reduced to
dP
d 2u
= 2
dx
dy

(2 1)

P
=0
y

and

Since u is not a function of x, the distance along the axis, there is no physical
mechanism to provide a change in the pressure gradient. Thus, the pressure
gradient is considered constant, i.e., dP/dx = constant. Note that under these
conditions all of the nonlinear convective terms in the momentum equations
are eliminated. Equation (21) is a linear second-order ordinary dierential
equation. It is easily integrated twice to yield
u(y) =

1 1 dP 2
y + Ay + B
2 dx

(2 2)

where A and B are integration constants.

In order to evaluate the integration constants, we apply the no-slip boundary
conditions at the channel walls. The boundary conditions at the walls are
y = h

(2 3)

u=0

Using these boundary conditions to evaluate the integration constants, we have

A=0
B=

1 h2 dP
2 dx

Substitution into Eq. (22) yields

u(y) =

1 h2 dP
2 dx

y
h

(2 4)

Here we see that the velocity distribution in the channel is parabolic and symmetrical about the axis. The maximum velocity, which occurs at the center of
the channel, is given by
1 h2 dP
um =
(2 5)
2
dx

26

Laminar Flow Analysis

Introducing nondimensional variables

u=

u
um

y=

y
h
(2 6)

u = 1 y2

yields

This nondimensional velocity distribution is shown in Figure 21. Using Newtons law of friction, which is obtained with the help of Eqs. (115a) and (125),
the shearing stress at the channel walls is given by
(y=h) =

du
dy

(2 7)

(y=h)

or using nondimensional variables

(=1) =
y

um d
u
h d
y

(2 8)

(=1)
y

Thus, using Eqs. (25) and (26), the shearing stress at the channel walls is
(=1) = h
y

dP
dx

(2 9)

24 Couette Flow
Continuing our discussion of the analytical solutions, consider the ow between
two parallel innite at surfaces, one of which is moving in its plane with a
velocity U (see Figure 22). The ow is considered steady, two-dimensional and
incompressible with constant properties. Using the same physical and geometric
arguments presented in the previous discussion of channel ow, the governing
equations again reduce to
d 2u
dP
(2 1)
= 2
dx
dy

Moving Surface
U

y
x

Figure 22. Couette ow.

Couette Flow

27

No-slip boundary conditions are assumed to apply at both the moving and the
stationary surfaces, i.e.,
y=0 u=0

y=h u=U

(2 10a, b)

For simple Couette ow the pressure gradient is assumed zero. Thus, the governing equation becomes
d 2u
=0
dy2
Integrating twice and evaluating the constants of integration from the boundary conditions yields the linear velocity distribution shown in Figure 22.
u
y
=
U
h

(2 11)

or u = y

Turning now to nonsimple Couette ow, i.e., when the pressure gradient is
nonzero, the governing dierential equation is now Eq. (21) with the boundary
conditions given in Eqs. (210a, b). Integrating Eq. (21) twice again yields
u(y) =

1 1 dP 2
y + Ay + B
2 dx

(2 2)

where A and B are constants of integration. Again, the pressure gradient is

constant. Using the boundary conditions to evaluate A and B yields
B =0
A=

1
1 1 dP 2
U
h
h
2 dx

Substituting into Eq. (22), introducing the nondimensional variables

u=

u
U

y=

y
h

h2 dP

P =
2U dx

(2 12a, b, c)

and rearranging yields

u = y [ 1 + P (1 y) ]

(2 13)

From Eq. (213) we see that the shape of the nondimensional velocity distribu
tion is determined by the nondimensional pressure gradient, P . Nondimensional
are shown in Figure 23.
velocity proles for several values of P
The results shown in Figure 23 indicate that the slope of the velocity prole

28

Laminar Flow Analysis

Moving Surface
y

10

1 0 1 2

10 P

Backflow
3

u
U

Figure 23. Nondimensional velocity proles for nonsimple Couette ow.

Hence, the shearing stress at the stationary surface is zero for P = 1 and
= +1. Using Eqs. (213) and (214) we see that
at the moving surface for P

for P < 1, regions of backow, i.e., u 0, exist near the stationary surface.

For P > 1, the velocity in the ow exceeds the velocity of the moving plate.
Physically, a region of backow exists when the force due to the momentum of
the uid in the ow direction is overcome by the adverse pressure gradient in the
ow direction. Similarly, velocities greater than that of the moving plate occur
when a favorable pressure gradient in the ow direction adds to the momentum
of the uid in that direction.
One further analytical solution, the suddenly accelerated plane wall, is presented below in order to illustrate several solution techniques.

25 The Suddenly Accelerated Plane Wall

Consider the two-dimensional parallel ow of an incompressible uid near a at
plate which is suddenly accelerated from rest and moves in its own plane with
An adverse pressure gradient means the pressure increases in the ow direction.

The Suddenly Accelerated Plane Wall

29

a constant velocity U (see Figure 24). Recall that the Navier-Stokes equations
for two-dimensional incompressible, constant property ow are
u
u
u
1 P
2u 2u
+
+u
+v
=
+
t
x
y
x
x2

30

Laminar Flow Analysis

which requires f(x) = constant. The mathematical model is then based upon a
uniform pressure eld for all x, y and t equal to or greater than zero.
From the continuity equation we have
u
v
=
=0
x
y
and the governing dierential equations reduce to
2u
u
= 2
t
y

t>0

(2 16)

The initial and boundary conditions are

u=0

for all y

u=U y=0
u0 y

t0
t>0

(2 17a, b)

where for t > 0 the rst boundary condition expresses the fact that there is
no slip at the surface. The second boundary condition is read: In the limit as
y approaches innity the velocity approaches zero. It represents the physical
condition that the inuence of the wall, i.e., the eect of viscosity, decreases
asymptotically to zero as the distance above the plate increases.
The governing dierential equation, Eq. (216), has the same form as the equation describing the diusion of heat by conduction in the space y > 0 when at
t = 0 the wall temperature is suddenly changed. The kinematic viscosity, ,
which appears in the governing equation is sometimes called the momentum diffusivity. It plays the same role in the momentum equation that the thermal
diusivity, , does in the energy equation. Equation (216) is generally called
the diusion equation. The diusion equation is rst encountered in a classical
course in dierential equations. There, an analytical solution is obtained by the
classical separation of variables technique. Thus, it seems reasonable to look for a
solution in the same manner. In an attempt to separate the variables we assume
u = F (t)G(y)

(2 18)

Substituting into the dierential equation yields

F G = F G
or

F
G
=
= i 2
F
G

(2 19)

where here the prime denotes dierentiation with respect to the appropriate argument. Since each side of Eq. (219) is a function of only one variable, they can

Similarity Analysis

31

individually be satised only if they are each equal to a constant, in this case
i 2. Hence, we have
F
(2 20a)
= i 2
F
G
i 2
=
G
The solution of Eq. (220a) is

F = C1 ei

(2 20b)
(2 21a)

and to Eq. (220b)

G = C2 e

i 2

1
2

+ C3 e

i 2

1
2

(2 21b)

The initial and boundary conditions are given by Eq. (217). There are three
possible cases, i = 0, 1. The case of i = 0 yields the trivial solution u = 0
everywhere and hence is discarded. For i =

32

Laminar Flow Analysis

We now seek a similarity solution. In particular, we seek a transformation of

variables which reduces the governing partial dierential equation to an ordinary
dierential equation. Since a partial dierential equation involves more than
one independent variable and an ordinary dierential equation only one, it is
reasonable to assume an independent variable transformation which attempts to
combine the two independent variables. Thus, we assume
(2 22)

= Bym tn

where B, m, n are, as yet, undetermined constants and is a transformed

independent variable. In addition, to nondimensionalize the equations we assume
a dependent variable transformation of the form
(2 23)

u = Af()

where A is again an as yet undetermined constant. Using these transformations

yields
u
f
(2 24)
=A
= ABnym tn1 f
t
t
where here the prime denotes dierentiation with respect to . Further
u

= A f = ABmym1 tn f
y
y
and

2u
= ABm(m 1)ym2 tn f + AB 2 m2 y2(m1) t2n f
y2

(2 25a)
(2 25b)

Substituting Eqs. (224) and (225) into the dierential equation, (Eq. 2-16),
yields
ABnym tn1 f = ABm(m 1)ym2 tn f + AB 2 m2 y2(m1) t2n f

(2 26)

Since A appears in each term it can be eliminated. Hence, its value is arbitrary
with respect to the dierential equation.
We now seek to select values of m, n, A, B such that reduction to a nondimensional ordinary dierential equation is achieved. Inspection of Eq. (226) shows
that for m = 1 the rst term on the right and the y dependence of the second
term are eliminated. Thus, we have
Bnyt(n1) f = B 2 t2n f
and

= Bytn

(2 27)
(2 28)

Using Eq. (228) allows Eq. (227) to be rewritten as

n t1 f = B 2 t2n f

(2 29)

Similarity Analysis

33

If n = 1/2 the time dependence is eliminated. Finally, choosing B 2 = 1/4 ,

where the 4 is introduced for later convenience, yields a nondimensional ordinary
dierential equation
f + 2f = 0
(2 30)
1 y

2 t

(2 31)

u = Af()

(2 32)

where

Now looking at the boundary conditions we have for t > 0, = 0 when y = 0

and when y . Hence
U
=0
f(0) =
(2 33)
A

f() 0
(2 34)
The rst of these boundary conditions, Eq. (233), is inconvenient in its present
form. Thus, we take A = U and have
=0

f(0) = 1

(2 35)

Note that the arbitrariness of A as revealed by the dierential equation is used

to achieve a simplied boundary condition.
A closed-form analytical solution to the boundary value problem given by
Eqs. (230), (234) and (235) is obtained by letting
=

df
=f
d

(2 36)

Upon substitution into Eq. (230) we have

+ 2 = 0

(2 37)

The solution of this ordinary dierential equation is (see [Murp60], Eq. 173,
p. 327)
2
= f = C1 e
(2 38)
Hence, after integrating

f = C1

e d + C2
2

(2 39)

The boundary condition at = 0 yields C2 = 1. The boundary condition as

yields
1
(2 40)
C1 =
2
e d
0

34

Laminar Flow Analysis

The solution is then

f() = 1

e d
2

(2 41)
d

Evaluation of the integral from zero to innity yields

2
f() = 1

e d
2

(2 42)

The second term on the right is the error function and 1 erf () is the complementary error function. Thus,
f() =

u
= erfc()
U

(2 43)

Here, a similarity analysis has been used to obtain a closed-form analytical solution to a problem which previously did not yield to any solution using separation
of variables.
Figure 25a gives y vs u/U for various values of time t. Note that the extent
of the viscous zone increases with increasing time. Due to the action of uid
viscosity, after an innite time the entire ow eld above the plate is moving with
the velocity of the plate. Figure 25b shows vs u/U . This gure illustrates
that the similarity variables collapse the solutions given in Figure 25a into a
single solution.
An exact solution to this boundary value problem can be obtained by numerical integration and compared with the analytical solution to illustrate the

t = 100

10
3
1
1
(a)

u
U

1
(b)

Figure 25. Velocity proles for the suddenly accelerated plate.

u
U

Two-Dimensional Stagnation Point Flow

35

accuracy of numerical integration techniques. Having demonstrated the utility

of similarity analysis, we proceed to apply this technique to problems which
generally do not yield closed-form analytical solutions (see Hansen [Hans64]).

27 Two-dimensional Stagnation Point Flow

Consider the two-dimensional steady ow of an incompressible viscous nonheatconducting uid impinging on a plate perpendicular to the ow direction (see
Figure 26). Assume that the ow at a large distance above the plate is given
by the corresponding inviscid (potential) ow and that no-slip conditions prevail
at the plate surface.
First, consider the inviscid ow solution. We take the plate to be at y = 0 and
the stagnation point at x = 0, y = 0 (see Figure 26). The ow is impinging on
the plate from the positive y direction. Under these circumstances the potential
ow solution yields the following expressions for the stream and velocity potential
functions (see [Rose63], p. 155)
(2 44)

i = K xy
i =

K 2
(x y2 )
2

(2 45)

where the i subscript indicates the inviscid ow solution. Dierentiation yields

Ui and Vi , the x and y components of the inviscid ow velocity, i.e.,
Ui =

i
=Kx
y

Vi =

(2 46)

i
= K y
x

(2 47)

0,0
Figure 26. Two-dimensional stagnation point ow.

36

Laminar Flow Analysis

In a potential ow the incompressible Bernoulli equation is applicable and is

used to obtain a relationship between the pressure at any point in the ow and
the stagnation point. Thus, we have
P0 Pi =

1
1
Ui2 + Vi2 = K 2 (x2 + y2 )
2
2

(2 48)

where the zero subscript indicates the stagnation point.

Returning to the viscous ow case, we have that the governing equations for
the dynamics of the two-dimensional steady ow of a viscous incompressible
constant property uid are
continuity
v
u
+
=0
x
y

(1 49)

u
2u
u
1 P
2u
+
+v
=
+
x
y
x
x2
y2

(1 50a)

v
2v
v
1 P
2v
+
+v
=
+
x
y
y
x2
y2

(1 50b)

momentum

The continuity equation is integrated by introducing a stream function such that

u=

v=

(2 49)

Hence, the momentum equations become

and

1 P
+ (xxy + yyy )
y xy x yy =
x
1 P
y xx x xy =
+ (xxx + xyy )
y

(2 50)
(2 51)

The appropriate boundary conditions at the plate surface, i.e., at y = 0, are

u = v = 0 or y = x = 0. At the stagnation point, i.e., at x = 0, y = 0, we
have P = P0 . Recall that at a large distance above the plate it was assumed
that the inviscid ow is recovered. Hence, as y , u Ui , v Vi . Thus,
u Kx and v Ky, and the pressure eld approaches the inviscid pressure
eld given by Eq. (248).
Notwithstanding our experience with the suddenly accelerated at plate, we
seek a solution by separation of variables. Thus, we assume
(x, y) = F (x)G(y) = KxG(y)

(2 52)

where K is a constant. Assuming a solution of this form is essentially equivalent

to predetermining the solution as a function of x. This form of the solution is

Two-Dimensional Stagnation Point Flow

37

justied on the basis of the required behavior at innity, i.e., as y , u =

y = KxG (y) Kx. Here the prime denotes dierentiation with respect to the
appropriate argument. Our result also implies that as y , G (y) 1. Since
Eqs. (250) and (251) represent two equations in the two dependent variables
and P , we must also make an assumption about the form of the viscous pressure
eld. Assuming that the functional form of the viscous pressure eld is closely
approximated by that for the inviscid pressure eld, we take
P0 P =

1
K 2 x2 + F (y)
2

(2 53)

Substitution of Eqs. (252) and (253) into the governing dierential equations
Eqs. (250) and (251) yields
K 2 xG G K 2 xGG = K 2 x + KxG
and

K 2 GG =

K2
F + KG
2

where here the prime denotes dierentiation with respect to the argument y. Dividing the rst of these equations through by K 2 x and the second by K 2 yields
G G GG = 1 +
and

GG =

G
K

F G
2
K

(2 54)
(2 55)

Immediately we note that the governing partial dierential equations have been
reduced to ordinary dierential equations. Since no equation which is a function
of x results, we know that we have correctly chosen the functional form with
respect to x in Eqs. (252) and (253).
From Eq. (252) the boundary conditions at y = 0 are
G(0) = G (0) = 0

(2 56)

The rationale for this assumption is illustrated by assuming

P0 P =

1
(u2 + v 2 )
2

Substituting for u and v as obtained from Eq. (252), we have

P0 P =

1
2
K 2 [(xG ) + G2 ]
2

Realizing that as y G 1, we nally assume for simplicity that the pressure

eld has the form given by Eq. (253). However, note that F (y) is taken as an arbitrary
function of y for greater generality.

38

Laminar Flow Analysis

At x = 0, y = 0, i.e., at the stagnation point, from Eq. (253) we have

P0 P =

1
K F (0) = 0
2

Since at x = 0, y = 0, P = P0
F (0) = 0

(2 57)

Recalling that at a large distance above the plate it is assumed that the inviscid
ow solution is recovered, we write
u = y = KxG U Kx
Hence as y
G (y) 1

(2 58)

Although an analytical solution of the two-point asymptotic boundary value

problem given by Eqs. (254) to (258) does not presently exist, exact numerical
solutions can be obtained. However, to obtain any insight into the results, it is
necessary to parameterize the solutions with respect to /K. This is inconvenient and also expensive in terms of computation time. Hence, before proceeding
further we look for ane stretching transformations for the dependent and independent variables which remove the factor /K, i.e., we nondimensionalize the
equations. Thus, we let
G(y) = f()

(2 59a)

F (y) = g()

(2 59b)

= y

(2 59c)

where , , and are nonzero constants.

Substituting into Eqs. (254) and (255) yields
f

ff =

ff =

1

+
f
2
()
K

1

g
f
2
2
K

(2 60)
(2 61)

Before choosing , and we look at the boundary conditions. The reason

is that a particular choice of , and which eliminates /K from the dierential equation may result in boundary conditions involving /K. Under these
circumstances, the solution requires parameterization with respect to the boundary conditions. Transformation of Eqs. (256) to (258) at = 0 yields
f(0) = f (0) = g(0) = 0

(2 62)

Two-Dimensional Stagnation Point Flow

and as
f () =

39

(2 63)

From an examination of Eqs. (260), (261) and (263) we see that an appropriate choice is = 1, / = K/ and = 22 . After solving for , , we have
=

1/2

1/2

=2

(2 64)

and the governing two-point asymptotic boundary value problem is

f +f + (1 f 2 ) = 0
g = f + ff

(2 65)
(2 66)

with boundary conditions

= 0 f(0) = f (0) = g(0) = 0

(2 67)

f () 1

(2 68)

Here, we see that not only are the momentum equations uncoupled from the energy equation but in addition the x and y momentum equations, Eqs. (265) and
(266), respectively, are also uncoupled. Thus, we can solve the x momentum
equation, Eq. (265), independently of the y momentum equation, Eq. (266),
and subsequently use the solution of Eq. (265) in obtaining that of Eq. (266).
That is, we use the values of f () and f () obtained from a solution of
Eq. (265) in Eq. (266) to solve for g().
Equation (266) can be directly integrated to yield
g=

f2
+ f + constant
2

From the boundary conditions given in Eq. (256) we see that the constant of
integration is zero, hence
f2
(2 69)
+f
g=
2
No known closed-form analytical solution of the remaining two-point asymptotic
boundary value problem is available. Thus, we look for an exact numerical solution. The procedure is to seek an exact solution of Eq. (2-65) and to subsequently
use that solution to obtain g().
Equation (265) is a third-order nonlinear ordinary dierential equation. Numerical integration of this equation requires a knowledge of f(0), f (0) and f (0)
A dierential equation is nonlinear if powers and/or products of the dependent variable
and/or its derivatives occur; e.g., y(dy/dx) + y = 0 is nonlinear, while x(dy/dx) + y = 0
is linear.

40

Laminar Flow Analysis

to start the integration. However, from the given boundary conditions only f(0)
and f (0) are known. The third required boundary condition is specied at innity. The procedure is to estimate the unknown value of f (0) and then perform
the numerical integration out to some large value of which we call max . max
is taken to be equivalent to innity. When the integration has proceeded to max
the value of f (max ) is compared with the required asymptotic value of one.
If f (max ) is within some specied small value of one then the outer boundary
condition is said to be satised and we have a solution of the governing two-point
asymptotic boundary value problem. If not, we estimate a new value of f (0)
and repeat the procedure. Since max = , f (max ) cannot equal one precisely.
Thus, we consider the outer boundary condition to be satised if
f (max ) = 1

(2 70)

.
where 1 is some small quantity, say 5 107 This is the so-called shooting method.

28 Iteration Scheme
In general the rst estimate for f (0) does not yield a solution. Arbitrary guessing of subsequent estimates of f (0), of course, proves to be quite inecient.
Hence, a logical method of determining the new estimates for f (0) must be
used. The Newton-Raphson method (see Appendix B) is frequently used for
estimating the unknown gradients needed to obtain numerical solutions to linear
and nonlinear dierential equations. Although the Newton-Raphson iteration
scheme assures convergence to the required outer boundary condition at max ,
it does not insure asymptotic convergence to the specied outer boundary condition required for this two-point asymptotic boundary value problem. Before
continuing the discussion, recall that asymptotic convergence implies that as
f () 1 its rst derivative approaches zero, i.e., f () 0 as max .
Considering the Newton-Raphson iteration scheme in this context reveals that
it does not insure asymptotic convergence; in fact, the Newton-Raphson iteration scheme might be satised by f (max ) = 1 + and (say) f (max ) = 0.5.
Hence, at max + , f () would not satisfy the outer boundary condition. In
addition, if the initial estimate of f (0) is very far from the correct value the
solutions tend to diverge. Under these circumstances the derivatives required in
the Newton-Raphson method cannot be calculated in any meaningful manner.
In order to eliminate the problems associated with the Newton-Raphson iteration scheme, the Nachtsheim-Swigert [Nach65] iteration scheme is used. This
technique is fully discussed in Appendix B. In summary the Nachtsheim-Swigert
iteration scheme is so structured that asymptotic convergence to the correct
outer boundary conditions is assured. This is accomplished by imposing the
additional condition that
(2 71)
f (max ) = 2

Numerical Solution

41

where 2 is some small quantity. The new estimates of f (0) are obtained such
that the sum of the squares of the errors, i.e., 1 2 + 2 2 , in satisfying the asymptotic boundary condition is a minimum. Thus, convergence is achieved in a least
squares sense.
A numerical algorithm to solve boundary value problems is discussed in detail
in Appendix D.

29 Numerical Solution
The results of a typical run of the stag2d program described in Section D12
are shown below.
f (0) = 1.25

0
6.

1.25
0.222783

0
1.57637

0
6.42524

0
1.61137

0
6.48878

0
1.03421

0
5.41737

0
1.00402

0
5.35978

0
1.00049

0
5.35303

0
1.00006

0
5.35221

0
1.00001

0
5.35211

0
1.

0
5.3521

f (0) = 1.251
0
6.

1.251
0.237166

f (0) = 1.23368
0
6.

1.23368
1.24093e2

f (0) = 1.23272
0
6.

1.23272
1.45136e3

f (0) = 1.2326
0
6.

1.2326
1.75541e4

f (0) = 1.23259
0
6.

1.23259
2.1323e5

f (0) = 1.23259
0
6.

1.23259
2.60013e6

f (0) = 1.23259
0
6.

1.23259
3.25885e7

Convergence achieved

42

Laminar Flow Analysis

Table 21. Solutions for two-dimensional stagnation point ow.

0
0.2
0.4
0.6
0.8
1.
1.2
1.4
1.6
1.8
2.
2.2
2.4
2.6
2.8
3.
3.2
3.4
3.6
3.8
4.
4.2
4.4
4.6
4.8
5.
5.2
5.4
5.6
5.8
6.

1.232588
1.03445
0.846325
0.675171
0.525131
0.398013
0.293776
0.211003
0.147351
9.99638e2
6.58254e2
4.20396e2
2.60203e2
1.55973e2
9.04887e3
5.07797e3
2.7549e3
1.44421e3
7.31269e4
3.57497e4
1.6868e4
7.67926e5
3.37239e5
1.42847e5
5.83719e6
2.30344e6
8.80753e7
3.29701e7
1.24507e7
5.12017e8
2.623e8

0
0.226612
0.414456
0.566281
0.685937
0.777865
0.846671
0.896809
0.932348
0.956834
0.973217
0.983853
0.990549
0.994634
0.997046
0.998424
0.999186
0.999593
0.999803
0.999908
0.999958
0.999982
0.999992
0.999997
0.999999
1.
1.
1.
1.
1.
1.

f
0
2.33223e2
8.80566e2
0.186701
0.312423
0.459227
0.622028
0.796652
0.979779
1.16886
1.36197
1.55776
1.75525
1.95381
2.153
2.35256
2.55233
2.75221
2.95215
3.15212
3.35211
3.5521
3.7521
3.9521
4.1521
4.3521
4.5521
4.7521
4.9521
5.1521
5.3521

g
0
0.226884
0.418333
0.583709
0.734742
0.88331
1.04013
1.21414
1.41233
1.63995
1.9007
2.19716
2.53101
2.90331
3.31474
3.76569
4.25637
4.78692
5.3574
5.96784
6.61828
7.3087
8.03912
8.80954
9.61997
10.4704
11.3608
12.2912
13.2616
14.2721
15.3225

Here the initial guess for the unknown initial condition f (0) for the x momentum
equation (Eq. 265) is taken as 1.25. After perturbing the initial guess for f (0)
in order to calculate the Nachtsheim-Swigert iteration derivatives at the edge of
the boundary layer, convergence within the required error values of f (max ) =
1 106 and f (max ) 1 = 1 106 occurs after six iterations.
Solution of the y momentum equation is obtained by using these results in
Eq. (269). The results are tabulated in Table 21 and shown in Figure 27.
Since no prior knowledge of the solution of equations of this type was available, the
value of f (0) = 1.25 was arrived at by initially limiting the value of max as discussed
in Appendix B.

Numerical Solution

43

44

Laminar Flow Analysis

'
f ()

f ''
(0)
0.8

1.50
2

1.25
1.23

1.20

6
1.0

2
Figure 28. Eect of f (0) on f ().

is somewhat improved in this case. These modications are left as an exercise

for the reader (see Problem 25).

210 Axisymmetric Stagnation Point Flow

The previously obtained solution for two-dimensional stagnation point ow can
be extended to the case of an axisymmetric stream impinging on a plane wall.
The solution obtained is representative of that near the forward stagnation point
of an axisymmetric blunt body. Figure 29 illustrates the problem under discussion. Here x and y are the radial and axial directions, respectively, with
u and v the velocity components in the x and y directions, respectively. The
plane is assumed to be perpendicular to the ow direction, with the stagnation
point as the center of the coordinate system. The momentum equations for
steady axisymmetric incompressible ow with constant properties are obtained
by transforming Eqs. (149) and (150) into cylindrical polar coordinates. They
are then (see Goldstein [Gold38], p. 143)
u

1 u
2u
u
u
1 P
2u
u
+
+v
=
+
2+
2
x
y
x
x
x x
x
y2

(2 72)

v
1 v
v
1 P
2v
2v
+
+v
=
+
+
x
y
y
x2
x x
y2

(2 73)

The axisymmetric continuity equation is

(xu) +
(xv) = 0
x
y

(2 74)

Axisymmetric Stagnation Point Flow

45

0,0
Figure 29. Axisymmetric stagnation point ow.

The boundary conditions are again taken to be no-slip at the plate surface
and no-mass-transfer through the plate surface. Further, it is required that the
inviscid solution be recovered at a large distance from the body. For axisymmetric inviscid stagnation point ow, the velocity components are given by (see
Moore [Moor64], p. 79)
Ui = Kx

(2 75)

Vi = 2Ky

(2 76)

and the pressure distribution by

P0 Pi =

1
K 2 x2 + 4y2
2

(2 77)

where the zero subscript indicates the stagnation point at x = y = 0.

The boundary conditions at the surface become u = v = 0 at y = 0. Further,
at the stagnation point, y = 0, x = 0, P = P0 . At a large distance from the
plate, i.e., as y , u Ui Kx, v Vi 2Ky.
The continuity equation is automatically satised by introducing a stream
function of the form
u=

1
x y

v=

1
x x

(2 78)

Alternatively, we may take U = Kx/2, V = Kx, in which case the pressure distrii
i
bution is
1
x 2
P0 Pi = K 2
+ y2
2
2

46

Laminar Flow Analysis

Introducing the stream function into the momentum equations yields

1
1
1 2
y xy 2 x yy 3 y =
x2
x
x
1 P
1
1
1

+
yyy + xxy 2 xy
x
x
x
x
and

1
1
1
x xy 2 xx y + 3 xy =
x2
x
x
1 P
1
1
1
1

xxx + xyy + 3 x 2 xx
y
x
x
x
x

(2 79)

(2 80)

By analogy with two-dimensional stagnation point ow, we seek solutions of

the viscous momentum equations by separation of variables. To this end we
assume that
(2 81)
(x, y) = Kxa G(y)
where a is an as yet undetermined constant. Choosing this form for (x, y)
presupposes that the x dependence is a power law. We again assume that
P0 P =

1
K 2 x2 + F (y)
2

(2 82)

Substituting into the momentum equations, the x momentum equation is

(a 1)K 2 x2a3 G 2 aK 2 x2a3GG =
K 2 x + Kxa1 G + a(a 2)Kxa3 G

(2 83)

and the y momentum equation becomes

2aK 2 x2a4 GG =

K2F
a(a 2)2 Kxa4 G + aKxa2 G
2

(2 84)

where here the prime denotes dierentiation with respect to the argument y.
Inspection of Eqs. (283) and (284) shows that if a = 2 these partial dierential
equations reduce to ordinary dierential equations. Thus, after dividing by K 2
and 4K 2 in the x and y momentum equations, respectively, we have

G + 2GG + 1 G 2 = 0
K
and

GG =

G
8
2K

(2 85)
(2 86)

Here, note that Eqs. (285) and (286) are very similar in form to Eqs. (254)
and (255). Again, the equations are nondimensionalized by seeking ane

Axisymmetric Stagnation Point Flow

47

stretching transformations of the form

G(y) = f()

(2 87a)

F (y) = g()

(2 87b)

= y

(2 87c)

where , , are nonzero constants. Substituting into Eqs. (285) and (286)
yields
f
and

2K
K 1
ff +
1 ()2 f

3

ff =

=0

1

g
f
8 2
2K

(2 88)
(2 89)

Before choosing particular values for , and the boundary conditions must
be investigated. The transformed surface boundary conditions are
=0

(2 90a)

f(0) = f (0) = g(0) = 0

The transformed boundary condition at innity is

f()

(2 90b)

Hence, in order to obtain a simplied boundary condition at innity we take

= 1. Finally, we take
2K
=1

= 82

and

Hence

2K

1/2

2K

1/2

4
K

(2 91)

The governing two-point asymptotic boundary value problem is then given as

f

+ ff +

1
(1 f 2 ) = 0
2

g =ff + f

(2 92)
(2 93)

with boundary conditions

=0

f(0) = f (0) = g(0) = 0

f() 1

(2 94)
(2 95)

48

Laminar Flow Analysis

As was the case for two-dimensional stagnation point ow, the transformed y
momentum equation is immediately integrable. Hence
g() =

f2
+f
2

(2 96)

Comparing Eqs. (292) and (265), we see that they dier only by the factor
of one-half multiplying the last term. The solution of Eq. (292) is obtained
by the modications to the eqmot2d routine used with stag2d. Specically,
eqmot2d becomes
equations of motion two-dimensional axisymmetric stagnation point ow
subroutine eqmot(eta, x(), param(), f())
f(1) = x(3)x(1) 0.5(1 x(2)x(2))
f(2) = x(1)
f(3) = x(2)
return
A typical run using the above modied eqmot2d routine is
f (0) = 1.23259

0
6.

1.232588
0.423879

0
3.29404

0
11.6197

0
3.30122

0
11.6401

0
0.881218

0
4.87142

0
1.01441

0
5.23484

0
0.998322

0
5.19087

0
1.0002

0
5.19599

f (0) = 1.233588
0
6.

1.233588
0.425169

f (0) = 0.912905
0
6.

0.912905
0.022935

f (0) = 0.929481
0
6.

0.929481
2.7727e3

f (0) = 0.92747
0
6.

0.92747
3.22936e4

f (0) = 0.927705
0
6.

0.927705
3.79112e5

f (0) = 1.23259

Axisymmetric Stagnation Point Flow

49

f (0) = 0.927677
0
6.

0.927677
4.31852e6

0
0.999977

0
5.19539

0
1.

0
5.19546

0
1.

0
5.19545

f (0) = 0.92768
0
6.

0.92768
6.21323e7

f (0) = 0.92768
0
6.

0.92768
4.34516e8

Convergence achieved
As shown for two-dimensional stagnation point ow, arbitrary estimation of
the value of f (0) can lead to results for the Nachtsheim-Swigert iteration derivatives, which do not yield convergence to the required outer boundary condition
f (max ) 1. Hence, it is desirable to have some knowledge of an approximate
range of values for f (0). If either an analytical or exact solution of a boundary
value problem similar to that under investigation is known, then it is reasonable
to use the value of f (0) obtained for the known solution as the rst estimate for
the problem under investigation. Considering the similarity of Eqs. (292) and
(265), the rst estimate of f (0) is taken to be that for two-dimensional stagnation point ow, i.e., 1.232588. Convergence to within f (max ) = 1 106 and
f (max ) 1 = 1 106 in seven iterations yields a value of f (0) = 0.927680.
Note that the correction value for the last iteration is less than 5 106. The
complete solution including the y momentum equation (see Eq. 296) is shown in
Table 22. The value of f (0) = 0.927680 agrees with that originally calculated
by Homann [Homa36].
The nondimensional velocity and pressure function results for axisymmetric
stagnation point ow are shown graphically in Figure 210. Also shown are the
results for two-dimensional stagnation point ow.
From an analysis of the previous assumptions and transformations we see that
for both two-dimensional and axisymmetric ow the local velocity components
and the stream function are related to their respective inviscid values by
u
= f ()
Ui

v
f ()
=
Vi

f()

=
i

(2 97a, b, c)

Further, the pressure eld for two-dimensional viscous stagnation point ow is

P0 P =

2 2 2g() 2
y
K x +
2
2

(2 98)

and that for axisymmetric viscous stagnation point ow

P0 P =

K 2
2g()
x2 + 4 2 y 2
2

(2 99)

50

Laminar Flow Analysis

Table 22. Solutions for axisymmetric stagnation point ow.

0
0.2
0.4
0.6
0.8
1.
1.2
1.4
1.6
1.8
2.
2.2
2.4
2.6
2.8
3.
3.2
3.4
3.6
3.8
4.
4.2
4.4
4.6
4.8
5.
5.2
5.4
5.6
5.8
6.

f
0.92768
0.82771
0.728152
0.630002
0.53477
0.444284
0.360449
0.284977
0.219151
0.163652
0.1185
8.30998e2
0.056379
3.69739e2
2.34217e2
1.43227e2
8.4508e3
4.80892e3
2.63824e3
1.39496e3
7.10674e4
3.48767e4
1.6484e4
7.50188e5
3.2869e5
1.38626e5
5.62723e6
2.19832e6
8.26459e7
2.99049e7
1.04212e7

f
0
0.175537
0.33111
0.466891
0.583305
0.681115
0.761462
0.825853
0.876098
0.914204
0.94225
0.962255
0.976069
0.985294
0.991248
0.994959
0.997192
0.998488
0.999213
0.999605
0.999808
0.99991
0.999959
0.999982
0.999993
0.999997
0.999999
1.
1.
1.
1.

f
0
0.017887
6.88836e2
0.149011
0.254348
0.381092
0.525629
0.684612
0.855027
1.03424
1.22004
1.41061
1.60453
1.80073
1.99843
2.19708
2.39631
2.59589
2.79567
2.99556
3.1955
3.39547
3.59546
3.79546
3.99545
4.19545
4.39545
4.59545
4.79545
4.99545
5.19545

g
0
0.175697
0.333483
0.477993
0.615651
0.753731
0.899605
1.0602
1.24163
1.44903
1.6865
1.95716
2.26332
2.60661
2.9881
3.40854
3.86835
4.36782
4.9071
5.48628
6.10542
6.76453
7.46363
8.20272
8.98181
9.80091
10.66
11.5591
12.4982
13.4773
14.4964

Referring to Eqs. (248) and (277) shows that the viscous pressure elds are
given by similar modications of the inviscid pressure elds. Looking at the
tabulated results, we see that close to the plate surface 2g()/ 2 > 1 and hence
the local pressure P is larger than in the inviscid case, whereas far from the
plate, i.e., for larger values of , 2g()/ 2 < 1. Hence, the local pressure P is
less than in the inviscid case. Further, the eect is less for axisymmetric than
for two-dimensional viscous stagnation point ow.
Returning to the stream function and velocity components, we see from the
numerical results that the viscous stream function and the velocity components

Axisymmetric Stagnation Point Flow

5
f'
()

0.8

f()

f ''
()

10

1.0

51

8
g()

0.6

f'
()

f()

0.4

0.2

Axisymmetric
Two Dimensional
0

g()

f ''
()

Figure 210. Nondimensional results for axisymmetric stagnation point ow.

are always less than the corresponding inviscid results. Further, the eect is
greater for axisymmetric than for two-dimensional stagnation point ow.
Finally, we investigate the shearing stress at the wall. From the results given
in Chapter 1, we see that the shearing stress at the wall for both axisymmetric
and two-dimensional stagnation point ow is given by
=

u
y

y=0

(2 100)

In transformed coordinates, the shearing stress at the wall for both two-dimensional and axisymmetric stagnation point ow is
= Kxf (0)

(2 101)

Hence, the shearing stress at the wall is greater for two-dimensional stagnation
point ow than for axisymmetric stagnation point ow.
The technique developed in this chapter for solving an asymptotic two-point
boundary value problem is useful in succeeding chapters, where solutions for the
velocity proles in laminar boundary layer ows are obtained.