Poiseuille Flow

Jean Louis Marie Poiseuille, a French physicist and physiologist, was interested in human blood
flow and around 1840 he experimentally derived a “law” for flow through cylindrical pipes. It’s
extremely useful for all kinds of hydrodynamics such as plumbing, flow through hyperdermic
needles, flow through a drinking straw, flow in a volcanic conduit, etc. For this reason, it is
generally known as “pipe flow”. Actually, the cgs unit of viscosity, the Poise (P), was named
after Poiseuille and is still used in many engineering texts. A single Poise is equivalent to 10
Pa-s (the SI unit), thus making the Pa-s measure of viscosity the “peso” of fluid dynamics and
Poise equivalent to the dollar.
This is the first of many special cases of Navier-Stokes equation in which very simplified
situations can be solved analytically. Pipe flow is defined to be unidirectional, i.e. there is only
a single non-zero component of velocity and that component is both independent of distance
in the flow direction and has the same direction everywhere. The geometry is that of a long
cylindrical pipe with length l and radius a so the appropriate coordinate system is cylindrical
polar (r, θ, z). The pressures at each end of the pipe are P1 and P0 so the pressure gradient,
dP/dz, is constant everywhere in the pipe. The unidirectional nature of the problem means
ur = 0 and uθ = 0, thus the continuity equation is reduced to ∂u ∂z
z
= 0. This means that because
of the incompressibility constraint, at any value of z the velocity must both be a constant value
as well as have an identical velocity profile. Furthermore, any change in the flow will occur
everywhere in the pipe instantaneously. Of course, you already know this is true because you
have taken a shower without a pressure regulator so that when somebody else flushes a toilet
and cold water is diverted to refill the toilet tank, the pressure gradient in the pipes for cold
water drops, decreasing the flow of cold water and exposing you to the hot water alone — ouch!!
However, even
in the more general case of the Navier-Stokes equation that has an inertial term,
∂uz ∂uz ∂uz
ρ ∂t + ∂z , one can see that for steady flow ( ∂t = 0) the geometry of the problem and the
incompressibility of the fluid specify that the inertial term is exactly zero. So Poiseuille Flow is
not limited to the Stokes regime, but also occurs at higher Re and we’ll see that this is important.
This 1-D version of the momentum equation in cylindrical coordinates is then
 
dP 2 dP 1 ∂ ∂uz
− + η∇ uz = 0 or − +η r =0 (1)
dz dz r ∂r ∂r
We will try a solution of the form
1 dP 2

uz = r + c1 ln r + c2 (2)
4η dz
subject to the boundary conditions of no-slip side walls and finite force over the fluid length
uz 6= ∞ at r = 0
(3)
uz = 0 at r = a

Solving for the constants we now have

1 dP 2
r − a2


uz (r) = (4)
4η dz
This means the velocity profile of the flow has a parabolic shape with a maximum in the center
(r = 0) and is zero at the pipe walls. Also note that the flow is independent of the fluid density.

1
 As mentioned, the velocity is maximum at the center, which we can calculate

−dP a2
umax = at r = 0 (5)
dz 4η
Pressure gradients are normally defined to be negative, such that water flows from high pressure
to low pressure, so when P1 > P0 , umax is a positive quantity. It is also useful to calculate the
total flow rate through the pipe, so we integrate the velocity over the a cross-section of the pipe
Z a Z a
1 dP 2 −dP πa4 πa4 (P0 − P1 )
r − a2 dr =


Q= 2πr u(r)dr = 2πr = (6)
0 0 4η dz dz 8η 8ηl

The volumetric flow rate (units of volume/time or m3 /s) shows that for a given pressure gradient
and viscosity, the flow through the pipe is proportional to the radius of the pipe to the fourth
power. This is what Poiseuille demonstrated experimentally. The mean velocity is simply the
total flow normalized by the cross-sectional area of the pipe

−dP πa4 −dP a2 1

ū = 2
= = umax (7)
dz 8πa η dz 8η 2
The mean velocity is the result of the net force exerted on the fluid by the pressure gradient acting
to overcome the viscous drag from the pipe walls. The force (per unit length) from pressure is

(P0 − P1 ) dP
FP = πa2 = −πa2 (8)
l dz
This shows that the mean flow, ū, is related to the pressure force by ū = FP /(8πη) and so it is
linearly inversely proportional to η. Similarly, for a Newtonian fluid, viscous drag is proportional
to the shear (tangential) stress, σzr , which we can evaluate near the wall of the pipe, r = a
 
∂uz −a dP −a dP
σzr |r=a = η |r=a = η = (9)
∂r 2η dz 2 dz

Similar to the non-dimensional drag coefficient of the Stokes sphere, cD , we can determine a
friction factor, f , which describes the effect of drag. We use the shear stress evaluated at the
wall as a characteristic stress and normalize that value by a characteristic pressure ( 12 ρf ū2 ) in
which we use the mean velocity
σzr |r=a −4a dP
f= 1 = (10)
ρ ū
2 f
2 ρf ū2 dz
If we substitute for just one of the ū, then f looks like
−4a 1 dP −4a 8η dP 32η
f= = −dP

= (11)
ρf ū ū dz ρf ū a2 dz dz ρf ūa

If we choose a characteristic length scale as the diameter of the pipe, D = 2a, then we have
64η 64
f= = (12)
ρf ūD Re

This relationship holds until the transition into the turbulent flow regime at Re ∼ 2000 − 3000.

2
Channel Flow
Another unidirectional flow is the flow between two rigid plates driven by a pressure gradient.
This is actually just Poiseuille flow in Cartesian geometry (with ẑ the same direction as in
cylindrical polar) so the pressure gradient and resultant flow are both only in the x direction
(uy = uz = 0) and the velocity profile varies with z. The geometry has the x-axis along the
mid-plane of the channel, and since the channel has height h, the channel walls are at ±h/2.
The governing equations are
∂ux
∂x
=0
(13)
2
− dP
dx
2
+ η∇ ux = 0 or − dP
dx
+ η ∂∂zu2x =0
Since dP/dx is constant, this is a second order O.D.E. and the integration is straightforward
1 dP 2
ux = z + c2 z + c1 (14)
2η dx
The boundary conditions are from the mirror symmetry along the mid-plane (ux (z) = ux (−z))
and no-slip at the walls (ux (z = ±h2
) = 0) We can now solve for the constants of integration and
get the velocity profile. All the same insights from Poiseuille flow in a pipe are applicable here.
1 dP  2
z − (h/2)2

ux = (15)
2η dx
The velocity profile is again parabolic in shape and constant everywhere.

Couette Flow
Couette flow is similar to channel flow and has the same geometry but with an important modi-
fication. Instead of the pressure gradient driving the flow, it is driven by the motion of one of the
boundaries and that motion is parallel to the direction of the channel ( dP dx
is in fact absent from
this problem). The assumption is that some external force is applied to move the wall and that
applied force simply scales with the viscosity of the fluid. Depending on the reference frame you
choose to do the problem in, the top or bottom plate can be moving at some velocity (U0 ) or they
can both move in opposite directions at (U0 /2). The most convenient choice for the coordinate
system is to have a stationary plate at z = 0 and a moving plate at z = h so again the channel
has height h. The governing equations for a shear driven flow are even simpler than for channel
flow since now dP
dx
=0
∂ux
∂x
=0
(16)
2 ∂ 2 ux
0 = η∇ ux or 0 = η ∂z2
Twice integrating this second order O.D.E. gives the solution ux = c1 z + c2 . The boundary
conditions are again no-slip velocity boundary conditions at the stationary and moving walls, so
ux (z = 0) = 0 and ux (z = h) = U0 . The solution for the ux is again a constant
z
ux = U0 (17)
h
The velocity profile in a shear driven flow is again identical for all values of x, varies linearly
with distance from the moving wall, and is independent of both density and viscosity. Also note
that the shear stress is also constant everywhere
∂ux U0
σxz = η =η (18)
∂z h

3
Classification of PDEs and types of Boundary Conditions
Any PDE can be classified using the method of characteristics which determines if the PDE
is either hyperbolic, elliptic, or parabolic. Both Laplace’s equation and Poisson’s equation are
classified as elliptical, and is a common class of equation one encounters in fluid dynamics. Other
examples include of the wave equation (hyperbolic) and the diffusion equation (parabolic). It is
important to understand which class of equation you are attempting to solve, in particular if you
are using numerical methods, because the stability or success of the numerical method applied
to one class of equation may be a completely unstable or be an unsuccessful approach if applied
to a different class of PDE.
The primary variable is the variable in the governing equation (either PDE or ODE) and every
primary variable always has an associated secondary variable. The secondary variable is usually
the derivative of the primary variable and always has a physical meaning that is often a quantity
of interest. In fluid dynamics the primary variable is velocity and the secondary variable is stress.
Another example is heat transfer in which the primary variable is temperature and the secondary
variable is heat flux.
In order to obtain a solution to any PDE, boundary conditions must be specified. There are
two types of boundary conditions that can be applied: those that specify the primary dependent
variable on boundary and those that specify a secondary variable on the boundary, and usually
the derivative is taken normal to the boundary. The first type of boundary condition is called
an essential boundary condition and when solving an elliptic class of equation it is known as
a Dirichlet boundary condition. The second type of boundary condition is called a natural
boundary condition and when solving an elliptic class of equation it is known as a Neumann
boundary condition. It is quite ok, and even somewhat common, to have mixed types of boundary
conditions along different parts of the boundary. For example, one portion of the boundary will
specify a Dirichlet boundary condition and another portion will specify a Neumann boundary
condition. However, it is impossible to specify both types of conditions at the same point of any
portion of the boundary. Thus, if the temperature is specified, the heat flux will be determined
(or vice-versa) but it can never happen that both are specified at the same place. Similarly, if
the stress is specified at a given point, the velocity will be solved for on the boundary at the
same point.
This is actually quite a powerful, and useful, thing to know, especially in situations like Couette
flow and channel flow, which have the same geometry. It is actually possible to combine the simple
solutions from both problems because 1) they are both linear ODEs we can use the principle
of superposition and 2) the solutions were arrived upon by applying the same type of boundary
condition. Both problems specified the velocity on the walls and therefore both applied Dirichlet
boundary conditions. We can then write the solution of Couette flow that now includes a pressure
gradient by simply transforming the channel flow solution to a coordinate system with the bottom
wall at z = 0 so
z 1 dP  2 
ux = U0 + z − hz (19)
h 2η dx
A simple model of asthenospheric counterflow is motivated by a shear flow driven by plate motions
on the surface. The shear flow sets up a pressure gradient in the the opposite direction which
drives an associated channel flow underneath the shear flow (a return flow). This is the same as
the above problem, except the direction of the pressure gradient is reversed
z 1 dP  2 
ux = U0 − z − hz (20)
h 2η dx

4
Steam Function
The stream function, ψ, is both an illustrative and useful approach to apply to fluid dynamics
as it can provide relatively quick solutions to 2-D incompressible flow problems. The major
drawback of the stream function is that it is basically limited entirely to 2-D incompressible flow
problems. The stream function is like a potential field in that only the difference in ψ between
two points has any physical meaning (the absolute value of ψ is arbitrary). Lines of constant ψ
are called stream lines and give an excellent visual representation of the flow, however, only in a
2-D geometry. In 2-D, the incompressibility constraint is from the continuity equation
∂ux ∂uy
∇ · u = 0 or + =0 (21)
∂x ∂y
The definition of the stream function is then
−∂ψ
ux = ∂y
(22)
∂ψ
uy = ∂x

The stream function satisfies the continuity equation

−∂ 2 ψ ∂ 2ψ
+ =0 (23)
∂x∂y ∂y∂x
The stream function can also be substituted into the Stokes equation
 3 
∂ ψ ∂3ψ
0 = −dP
dx
− η ∂ 2 x∂y
+ ∂3y

  (24)
−dP ∂3ψ ∂3ψ
0= dy
+η ∂3x
+ ∂ 2 y∂x

Now eliminate the pressure term using the same technique that was applied earlier when solving
for the flow around a Stokes sphere, i.e. take partial derivatives w.r.t. the other dimension
h  3 i
∂ −dP ∂ ψ ∂3ψ
0 = ∂y dx
− η ∂ 2 x∂y
+ ∂3y

h  i (25)
∂ −dP ∂3ψ ∂3ψ
0= ∂x dy
+η ∂3x
+ ∂ 2 y∂x

and then subtracting the resulting equations we get

∂ 4ψ ∂ 4ψ ∂ 4ψ
0= + + (26)
∂ 4 x ∂ 2 x∂ 2 y ∂ 4y
Rearranging the derivatives we now have
 2
∂2
 2
∂2

∂ ∂
0= + + ψ (27)
∂ 2x ∂ 2y ∂ 2x ∂ 2y
This equation can be recognized as the Laplacian operator (∇2 ) being applied twice to ψ
0 = ∇2 ∇2 ψ



(28)
2
And this is known as the Biharmonic operator ((∇2 ) = ∇4 ) which we can use to write
0 = ∇4 ψ (29)
There are well-known solutions to this equation and it is also valid for non-Cartesian geometries.

5
Corner Flow
The situation of a subduction zone is in some ways analogous to one variation of the classic corner
flow problem in fluid dynamics. In this version, two rigid plates (infinite in extent) converge at
a point where the advancing plate (plate A) dips at an angle below the back-arc plate (plate
B). We will use the point of convergence as the origin of a 2-D cylindrical coordinate system
with plates on the surface (the line at θ = 0). The angle that plate A makes between itself on
the surface and the dipping portion is defined as θa and the “’dip angle” between plates A and
B is defined as θb (and assumed to be acute). Plate B is assumed to remain stationary while
plate A is moving on the surface at velocity ur = −U0 (towards the origin) and along the dip
angle at ur = U0 (away from the origin). For both plates, the velocities in the θ direction are
assumed to be zero (uθ = 0). Notice that there are no body forces in this problem, and that the
Stokes flow is driven entirely by the velocity boundary conditions (which themselves are driving
by some applied force but since it is not a body force it is irrelevant). The governing equations
for Stokes flow are simply ∇ · τ = 0 and ∇ · u = 0. Expanding the momentum equation out
into the components of total stress
∂τrr 1 ∂τθr
∂r
+ r ∂θ
=0
(30)
∂τrθ 1 ∂τθθ
∂r
+ r ∂θ
=0
We can use the constitutive relationship between total stress and strain rate, τ = −P I + 2ηD

τrr = −P + σrr = −P + 2η ε̇rr

τθθ = −P + σθθ = −P + 2η ε̇θθ (31)

τrθ (= τθr ) = σrθ = 2η ε̇rθ

And rewrite the total stress with the strain rate having terms of the velocity gradients
τrr = −P + 2η ∂u
∂r
r

1 ∂uθ ur


τθθ = −P + 2η r ∂θ
+ r
(32)

1 ∂ur ∂uθ uθ


τrθ = η r ∂θ
+ ∂r
− r

Notice that if we add the normal components of stress together we get

 
∂ur 1 ∂uθ ur
τrr + τθθ = −2P + 2η + + (33)
∂r r ∂θ r
The 2nd term on the RHS vanishes since ∇ · u = 0, and because the fluid is isotropic the
expressions for normal stresses become
τrr = −P
(34)
τθθ = −P
Using these allows us to express the momentum equation entirely in terms of P and τrθ
− ∂P
∂r
+ 1 ∂
τ
r ∂θ rθ
=0
(35)
− 1r ∂P
∂θ
+ ∂
τ
∂r rθ
=0

6
and this can be rewritten as
∂uθ uθ
− ∂P + η 1r ∂θ
∂ 1 ∂ur


∂r r ∂θ
+ ∂r
− r
=0
(36)
∂uθ uθ
− 1r ∂P ∂ 1 ∂ur


∂θ
+ η ∂r r ∂θ
+ ∂r
− r
=0

Seems all that manipulation didn’t help as the momentum equation still looks a little ugly.
Luckily, it was shown that solving ∇4 ψ = 0 will also give the solution for velocity, and it turns
out the stream function is a more convenient way to approach the problem. In 2-D cylindrical,
the Laplacian is
1 ∂ 2ψ
 
2 1 ∂ ∂ψ
∇ψ= r + 2 2 (37)
r ∂r ∂r r ∂θ
Considering the geometry of the problem has plates of infinite extent with constant relative
velocity, the solution for velocity everywhere is expected to be independent of r. This means the
equation is separable and we will use a solution of the form

ψ = R(r)T (θ) and (38)

1 ∂ψ
ur = r ∂θ
(39)
−∂ψ
uθ = ∂r
Simple substitution for ψ gives
1 ∂T (θ)
ur = R(r) (40)
r ∂θ
which means R(r) = r and then ψ = rT (θ) which upon substituting back into the Biharmonic
equation gives
∂ 4T ∂ 2T
+ 2 +T =0 (41)
∂θ4 ∂θ2
The 4th order PDE has now been reduced to a 4th order ODE which has a general solution of
the form
T (θ) = A sin θ + B cos θ + Cθ sin θ + Dθ cos θ (42)
and there are also 4 boundary conditions but these are given as velocities so we need ur and uθ
∂T (θ)
ur = ∂θ
= A cos θ − B sin θ + C(sin θ + θ sin θ) + D(cos θ − θ sin θ)
(43)
uθ = −T (θ)

At this point its a good idea to break the problem into two portions and solve for the stream
function in each domain. The obvious choice for the two domains is the “back-arc region” formed
by the (acute) dip angle between the subducting plate and overriding plate and the “fore-arc
region” underneath the subducting plate. The flows are identical along the boundary of the
subducting plate, and this line is known as the separatrix. The boundary conditions are then

ur (θ = 0) = −U0 in the fore-arc region

ur (θ = 0) = 0 in the back-arc region
uθ (θ = 0) = 0 in both regions (44)
uθ (θ = θb ) = 0 along the separatrix
uθ (θ = θb ) = U0 along the separatrix

7
Each region has 4 boundary conditions to solve for the 4 unknowns constants, and after a lot of
algebra one arrives at the solution
−rU0 [(θa −θ) sin θ−θ sin(θa −θ)]
ψa = θa +sin θa
or more simply, ψa = −rU0 fa (θ) in the fore-arc region
(45)
rU0 [(θb −θ) sin θb sin θ−θb θ sin(θb −θ)]
ψb = θb2 −sin2 θb
or more simply, ψb = rU0 fb (θ) in the back-arc region
The velocities in each region are readily obtained through differentiation of ψ: ur = −U0 fa0 (θ)
and uθ = U0 fa (θ) in the fore arc and ur = U0 fb0 (θ) and uθ = −U0 fb (θ) in the back arc. In order
to obtain the pressure, we need to go back to the momentum equation and use the fact that
τrr = τθθ = −P which itself is related to the shear (tangential) stress
τrθ = η Ur0 [−fa00 (θ) − fa (θ)] in the fore-arc region
(46)
τrθ = η Ur0 [fa00 (θ) + fa (θ)] in the back-arc region
This helps obtain the pressure solution
2U0 η [sin θ−sin(θ−θa )]
Pa (r, θ) = r θa +sin θa
in the fore-arc region
(47)
−2U0 η [θb sin(θb −θ)−sin θb sin θ]
Pb (r, θ) = r θb2 −sin2 θb
in the back-arc region
Inspection of these solutions reveals that Pa is always a positive quantity and Pb is always
a negative quantity. A positive pressure below the subducting plate implies compression or
upward force on the surface. A negative pressure in the mantle wedge indicates that there is a
suction between the subducting plate and overriding plate. This corner flow suction acts as a
hydrodynamic lift that is proportional to the pressure difference above and below the slab. The
lift is found by integrating P (r, θ) along the dip angle, θa , over a length l. The torque exerted
by lift is balanced by gravity through the weight of the slab with thickness h and density ∆ρ.
h i
sin2 θb
Rl
Tf low = 0 [Pa (r, θa ) − Pb (r, θb )] r dr = 2U0 ηl (π−θsin θb
b )+sin θb
+ 2 2
θb −sin θb
(48)
Tgravity = 12 ∆ρghl2 cos θb
Both torques can be normalized by a characteristic torque, 2U0 ηl, which then allows one to find
the critical dip angle, θc that determines when the torque derived from gravity is balanced by
the lift generated by circulation in the mantle wedge. For any angle smaller than θc , the torque
exerted on the slab by mantle flow will exceed the weight of the slab, and assuming the velocities
remain constant, a positive feedback will occur such that θ decreases to zero. This critical angle
was determined by Stevenson and Turner, (1977), to be 63◦ for which they found the net torque
was about 2 times the characteristic torque. Assuming a 100 km thick slab that is 600 km in
length subducts at 6 cm/yr and has ∆ρ =80 kg/m3 gives
∆ρghl
2∼ which can be used to estimate the upper mantle viscosity, η = π × 1021 Pa s (49)
4ηU0
Clearly, θc = 63◦ is too large as many slabs are observed to have dip angles shallower than this
estimate, so obviously there must be many other important factors. One of the more important
factors is the non-Newtonian rheology of the mantle wedge as studied by Tovish et al. (1978)
who found this reduced θc = 54◦ for a power law fluid with n=3. There are also reasons for θc to
be larger, as slabs with finite lateral extent allow for a 3-D component of the mantle flow (i.e the
toroidal flow) around slab edges which reduces the pressure differential Dvorkin et al. (1993).