## Stokes's and Poiseuille's Laws

Here are two interesting examples of the flow of viscous fluids

### The Results

Let's begin by stating the formulas that are the goal of this analysis. They deal with the movement of viscous fluids, commonly air and water, and are some of the rare analytical results of fluid mechanics of practical value. We will study the dynamics of viscous fluids in this article, and show how these results can be obtained by analysis.

Poiseuille was interested in the flow of a fluid, actually blood, in narrow tubes, and found that the discharge Q in cc/s was related to the pressure difference Δp in dynes/cm2 by Q = (πa4/8Lμ) Δp. Here, a is the radius of the tube, and L its length, both in cm, and μ is the absolute viscosity of the liquid in dyne-s/cm2 or poise. The viscosity of water is 0.0101 poise at 20°C. Blood, with which Poiseuille was concerned, is not a simple Newtonian fluid, and its flow is a complicated problem, but water, oil and such obey the equation very well. The catch is that the flow must be rather slow, or the flow will become turbulent, obeying different rules and the resistance to flow will increase. The flow is guaranteed to be laminar (not turbulent) provided that the Reynolds number R = DVρ/μ is below 2000, as found from experiment. Here, D is the diameter of the tube, 2a (cm), V is the average velocity of flow Q/πa2 (cm/s), μ the absolute viscosity, and ρ the density in gm/cc. This number is dimensionless, and is the same in any consistent system of units. For a 1/2" pipe, R is 2000 for an average velocity of 15.7 cm/s of water, or Q = 20 cc/s, which is not a very rapid flow. We should also be aware that the laminar flow is not established immediately at the entrance to the tube, but takes a little length for this to happen, in which the liquid at the center speeds up, and that at the wall of the pipe is retarded.Stokes solved the problem of the steady motion of a sphere with velocity V through a still viscous fluid (or the flow of the fluid past a stationary sphere, which is the same thing). He found that the force exerted on the sphere by the fluid was F = 6πμaV, where F is in dyne, a in cm, V in cm/s, and μ is the absolute viscosity. Note that the force is proportional to the velocity, the rule we often assume for a frictional force. Again, the motion must not be too fast, and a Reynolds number R = DVμ/ρ can be used to decide what is too fast. This has the same form as in Poiseuille's equation, but the flow is different, and the critical value may be different. If V is too great, again the flow becomes turbulent, with the formation of a turbulent wake behind the sphere, and now the force is proportional to the square of the velocity. We usually write F = C A(ρV2/2), where C is the dimensionless drag coefficient. Where Stokes's formula applies, C = 24/R. The critcal Reynolds number is about 1 for air, and the drag coefficient for high R is about 0.4, dropping to about 0.2 when the boundary layer changes from laminar to turbulent. This is a complex and interesting subject that you can read more about in the references.

The streamlines in Stokes flow dip in toward the sphere, so to avoid disturbing the flow pattern, there must be no nearby boundaries. A ball dropping in a tube about twice its diameter would be significantly affected (the drag would increase). Its area of application is really to rather small spheres in air, such as cloud particles, or to slow motion in a viscous fluid, such as glycerin or oil. It was famously used by R. A. Millikan in the determination of the electronic charge, to measure the weight of the drops by how fast they fell. The drops were so small that a correction for the mean free path of air molecules was necessary for an accurate result.

The problem can also be solved neglecting viscosity, and the results are entirely different. Some are applicable to reality, but others are not. When a sphere moves in an ideal fluid, it suffers no drag whatsoever, but once started in motion, would continue forever. The effective mass of the sphere is increased by half the mass of the fluid displaced, a result that is quite valid and observable. It can be applied to the motion of bubbles in beer, for example. Most real-life problems of fluid flow depend greatly on viscosity, even when it would seem to have a small effect. The drag on a body moving in a fluid is mostly caused by the lower pressure in the turbulent wake behind it, and the drag is reduced by making this wake as small as possible. The pressure on the front of the body has only a small effect. Trucks with the "streamlining" in front of the trailer would do better to put it behind!

We recall that the coefficient of viscosity μ is defined by Newton's equation F/A = μ(dv/dy), where F is the shearing force on an area A, and dv/dy is the velocity gradient normal to the surface. If μ is a constant independent of the velocity of shearing, the fluid is called Newtonian. Many, if not most, fluids met with in engineering are Newtonian. Some, like paint and blood, are not, and there is a variety of interesting behaviors in some other substances, especially those that are a mixture of phases. μ is measured in poise, dyne-s/cm2.

### Stress and Strain

To begin our study, we must review the nature of stress and strain. This will require the use of quantities with a number of components in three-dimensional space, such as vectors, as well as more complex quantities. These can be represented explicitly, and were indeed for many years, but this involves tedious repetitive writing. Vector analysis provides an excellent shorthand, and is widely understood and used. For more complex quantities, this was extended to what were called dyadics, but these proved rather cumbersome, and seem to have disappeared. What turns out to be very convenient is representation by Euclidean tensors using indices. This essentially writes all the components, as was originally done, but in a very compact and useful form. The vector a becomes the quantity ai, where we understand that the index i takes the values 1, 2 and 3, or x, y and z, or whatever. The scalar produt of a and b is aibi, where the repeated index implies summation over its range of values, called contraction. If you are not familiar with this, a short course can be found at Euclidean Tensors. I shall use the index notation in what follows.

Suppose that xi is a point in space occupied by a solid body, and that ξi is the displacement of that point when the body is strained. Then, εij = ∂iξj is the strain, or ratio of displacement to distance. The symbol ∂i is short for ∂/∂i. This is a new kind of quantity, depending on two directions. Suppose this region of the body isn't really strained, just rotated bodily by an angle θ. Considering the displacements that result from such a rotation, it is clear that ε ij = -εji = θ. We really don't want these if we are studying strains, and can get rid of them by forming Φij = (εij + εji)/2, which is called symmetrizing Φ. What's left is just the pure strain, described by the symmetric tensor Φij = Φji. There are six independent components of this tensor, three diagonal and three off- diagonal. Among other things, this means that the order of the indices is immaterial, so we do not have to worry about it. In a fluid, instead of the displacements, we use the rate of displacement, the time derivative of ξi instead (which is just the fluid velocity v), so Φ becomes the rate of strain tensor.

Strain causes stress, so we must find some way to specify stress. If we imagine a plane cutting the solid, then the part of the solid on one side of the plane exerts a force on the other side. The ratio of this force to the area considered is the stress. Given a normal unit vector to specify the orientation of the plane, then the stress is another vector that is a function of the direction of the normal vector, as well as of position. We expect, then, that stress can be described by a quantity Ψij. By considering the equilibrium of a small volume, the fact that it must not rotate is reflected in the symmetry of Ψ. Like Φ, then, &Psi is symmetric and has six independent components.

In a fluid at rest, the only stress that can exist is a pressure normal to any area, and this pressure is the same in any direction. This is, in fact, the definition of a fluid. Therefore, considering areas normal to the coordinate axes, we find that -pnj = Ψijnj, where ni is a normal vector, and a pressure is considered a negative stress. This means that Ψij = -pδij, where δij is the Kronecker delta, which is 1 when the two indices are equal, 0 otherwise. In solids or in viscous fluids, the stress may also have a component in the plane, called a shear stress. In this case, Ψ will acquire off-diagonal components.

In an isotropic solid, the relation between stress and strain, or between Ψ and Φ, takes the form Ψij = [kΦkk -(2μ/3)Φkkij + 2μΦij. The quantity Φkk is the sum of the diagonal elements of Φ, also called its trace. It is equal to ΔV/V, the relative volume expansion of the material, and called the dilatation. The trace does not depend on the orientation of the axes, which should be clear from its meaning. In general, each of the six components of Ψ should depend on the six components of Φ, giving 36 constants of proportionality if we assume Hooke's Law, that stress is proportional to strain. If the medium is isotropic, or the same in any direction, then these reduce to only two constants, the modulus of expansion k, p = -k(ΔV/V) (positive p causes a decrease in V), and the modulus of shear, μ. Young's modulus, for example, is Y = 3kμ/(k + μ/3) and Poisson's ratio is σ = (3k - 2μ)/(6k + 2μ). This is much simpler than it appears. Try contracting Ψ and some normal vectors and see what stresses you get, for simple Φ's.

In the case of a fluid, we find analogously that Ψij = [-p - (2μ/3)Φkkij + 2μΦij, in terms of the pressure p and the absolute viscosity μ, which relates a shear stress to the rate of shear strain. The analogy is clear. Φkk = Tr Φ = ∂kvk = div v is the dilatation. For an incompressible fluid, this quantity is zero, and the stress tensor become even simpler. Unlike solids, fluids are almost always isotropic, so all this is less of an approximation. We are now able to express the stresses in a viscous fluid in terms of the fluid velocity.

### Equations of Motion

The motion of a fluid in the Eulerian picture is defined by the velocity field v or vi, which is the velocity of the particle of fluid that happens to be at P(x,y,z) at time t. A time dt later, this particle has moved off a distance vi t and has been replaced by another. If we are interested in the time rate of change of some property of the particle of fluid, which we shall call F, then dF/dt must "follow the particle," or dF/dt = ∂F/∂t + viiF. This is called the substantial derivative, and is one of the first things explained in fluid mechanics.

The mass of fluid must be conserved in the motion, or dρ/dt + ρ∂ivi = 0, where the dilatation appears again. To find how ρ as a point function behaves, we use the substantial derivative to find ∂ρ/∂t + viiρ + ρ∂ivi = 0, or ∂ρ/∂t + ∂i(ρvi). This is the equation of continuity that must hold for any possible motion.

Newton's Second Law for a particle of fluid is ρ(dvi/dt) = ρFi + ∂jΨij, where we have used the divergence theorem to get the force per unit volume from the stress tensor. ρFi is an external force per unit mass, such as gravity. The total force on the particle is the sum of the external force and the force given by the stress tensor. Using the expression of Ψ in terms of Φ, and the substantial derivative, we find: ρ(∂vi/∂t) + ρvkkvi = ρFi - ∂ip + (μ/3)∂i(∂kvk) + μ∂kkvi. This is called the Navier-Stokes equation, derived by Navier (1822), Poisson (1829), Saint-Venant (1843) and Stokes (1845) in successively more modern form. The second term on the left-hand side makes it nonlinear in v, which is extremely unpleasant if we are looking for analytical solutions.

If we specialize to the slow motion of an incompressible fluid, the Navier-Stokes equation becomes a bit simpler: ρ(∂vi/∂t) = ρFi - ∂ip + μ∂kpart;kvi.

### Vorticity

As elements of a solid body can rotate, so can elements of a fluid. If we antisymmetrize the rate-of-strain tensor, that is form Ωij = (εij - εji)/2, we eliminate the pure strain and are left with the rotation. This is an antisymmetric tensor with three independent (off-diagonal) components. These three components form a vector (1/2)ω = (1/2) curl v. curl v is the vorticity of the motion, a field quantity like any other. From Stokes's Theorem, the line integral of the velocity around a closed curve is equal to the flux of vorticity through it. The diagram shows how the antisymmetric part of the strain is related to the angle of rotation. Taking the time derivative gives the relation between vorticity and curl v. The lines of unit length and the strains in the diagram are actually multiplied by some small number so that the change in the lengths of the axes is infinitesimal, and the strains perpendicular to both. This is hard to represent clearly in a diagram.

If we take the curl of the Navier-Stokes equation (in its slow, incompressible form), we find that ∂(curl v)/∂t = μ div grad curl v, since curl commutes with div grad (write it out, if you don't believe it). This is the same equation satisfied by the diffusion of heat for each component of the vorticity. We see that vorticity diffuses through the fluid with the diffusion coefficient equal to the viscosity. If the fluid is non-viscous, then ∂(curl v)/∂t = 0, and curl v = constant. This is a simple form of Kelvin's Vorticity Theorem, stating that for an ideal fluid, the vorticity is a constant. If the motion began from rest, then this says that the vorticity will be zero everywhere ever after, and the flow will be irrotational. Even the smallest amount of viscosity will allow vorticity to be created, and it may dominate the flow pattern. In the present case, with Poiseuille's and Stokes's Laws, we expect vorticity.

### Flow Around a Sphere

The problem is shown at the right. A rigid sphere of radius a is placed in a uniform flow of velocity V in the x-direction at large distances. The problem is symmetrical about the x-axis, so the velocity in the xy-plane has no z component, and conditions are the same in any plane through the x-axis. We assume the flow to be steady, slow and incompressible. The velocity must satisfy the equation of continuity, div v = 0, Newton's Second Law, div grad v = (1/μ) grad p, and must vanish on the surface of the sphere. At great distances, its x-component must be V.

This expression for v gives the required vorticity, satisfies the equation of continuity, and now if it can be made to satisfy the boundary condtions by proper choices of A, B and C, it will be the desired solution. The behavior at large distances requires B = V. The vanishing of the velocity on the surface of the sphere demands that A = -3aV/4 and C = -a3 V/4. For the components of the velocity, we have vr = V(1 - 3a/2r + a3/2r3) cos θ and vθ = -V(1 - 3a/4r - a3/4r3) sin θ. The z-component of curl v in the xy-plane is -(3aV/4r2) sin θ.

Now we can use the stress tensor to find the stresses on the surface of the sphere, and by integrating them over the sphere to find the drag. This is made much easier by choosing special x,y and z axes at a point on the sphere at angle θ with the x-axis, with the x-axis radially outward, the y axis in the direction of increasing θ, and the z-axis normal to the xy-plane, as shown in the diagram. Then, ∂vx/∂x = ∂vr/r∂θ, ∂vy/∂y = ∂vr/r∂θ + vr/r, ∂vy/∂x = ∂vθ/∂r, ∂vy/∂y = ∂vθ/r∂θ - vθ/r. Don't forget that in the second and fourth equation the unit vector also changes, and gives a contribution.

The only contributions that do not vanish on the surface of the sphere are the normal pressure, P = (3μV/2a) cos θ, and the shear stress U = μ(∂vy/∂x) = -(3μV/2a) sin θ. The net stress on unit area is P cos θ - U sin θ = 3μV/2a in the x-direction, so that the total drag becomes 4πa2(3μV/2a) = 6πμaV, the desired result. P and U can be separately integrated over the area, of course, and it is found that the pressure contributes 1/3, the shear 2/3 of the total drag. Note that there is a positive pressure over the face of the sphere facing the flow, and a negative pressure over the other face. The shear force is a maximum at the equator, and is zero on the axis.

Page, and other authors, copy Lamb's pictures of the stream functions for the irrotational flow around a sphere and for the viscous flow, and remark on how very different they are, repeating Lamb. Indeed they are, but they refer to two different situations, the first where the sphere is moving, and the second where the fluid is moving. The actual difference is not as great as may appear from these diagrams, but, of course, there is a significant difference nonetheless.

### Flow Through a Tube

Poiseuille's problem is quite a bit easier. We can solve it using the definition of viscosity and a little statics. Consider a cylinder of fluid of thickness dr at a radius r, and of length L, and let this cylinder move at velocity v = v(r). The fluid on the side towards the wall moves a little slower, that towards the center a little faster. The velocity gradient causes the outer fluid to retard our cylinder, the inner fluid to urge it forward. In equilibrium, the velocity gradient will arrange itself so that the net force is balanced by the pressure difference at the ends of the cylinder. This is expressed mathematically by (d/dr)[2πrLμ(dv/dr)]dr = -2πrΔp dr. Integrating once, we find 2πrLμ(dv/dr) = -2π(r2/2)Δp, or dv/dr = -(Δp/2Lμ) r + A, where A is a constant. By symmetry, when r = 0, dv/dr = 0, so A = 0. Integrating a second time, v(r) = -(Δp/4Lμ)r2 + B = 0, where B is a second constant of integration. Since v = 0 when r = a, the evaluation of B gives v(r) = (Δp/4Lμ)(a2 - r2), which is the parabolic velocity profile of Poiseuille flow. As with the sphere, the vorticity forms rings around the x-axis. It is easy to find the vorticity from the velocity profile.

The rate of discharge Q is found by integrating the velocity over the cross section, Q = 2π∫v(r)rdr. The integral is easily performed after making the substitution u = r2, with the result that Q = (πa4/8Lμ)Δp. The rate of discharge is proportional to the pressure difference Δp divided by the length L, or the pressure gradient, to the fourth power of the radius of the tube, and inversely proportional to the absolute viscosity. The mass discharge is proportional to the kinematic viscosity μ/ρ.

### Flow in a Pipe

Flow in cylindrical pipes is a subject of great importance in engineering, and a great deal of empirical data has been built up that allows the solution of most practical problems with relative ease. For details, see Daugherty in the References. We will only make the connection with Poiseuille flow here.

In general, practical flow in pipes is turbulent, not laminar. The engineer speaks of "head loss" instead of pressure loss. Pressure and head are directly related through p = ρgh = γh, where p is pressure in dyne/cm2, h is head in cm, ρ is the density in gm/cm3 and γ is specific weight in dynes/cm3. Of course, engineers use other units, such as lb/ft2, ft, slug/ft3 and lb/ft3. The head can be measured directly by a manometer at any point, though actual manometers are seldom used these days. A certain volume rate of flow Q through a length L of pipe is associated with a loss in head of h/L, a dimensionless ratio.

The engineer uses the equation h = f(L/D)(V2/2g), where D is the inside diameter of the pipe, V is the average velocity Q/A, where A = πD2/4, h the head loss in a length L of pipe, and f the Fanning friction factor. The friction factor is a function of the Reynolds number DVμ/ρ and, for turbulent flow, of the relative roughness of the pipe, e/D, where e is the absolute roughness. The usual connection is given in the famous chart by L. F. Moody, which can be found in Daugherty. Turbulent flow is usual when the Reynolds number exceeds 2000. This is much larger than for a sphere, because the pipe tends to keep things smooth. Below 2000, Poiseuille flow governs.

Poiseuille's equation, expressed in terms of head loss, is h = 32μLV/γD2. Comparing this with the Fanning equation, we find that f = 64/R, where R is the Reynolds number. This is represented by a line in the Moody chart showing that f decreases steadily to about 0.032, where the critical Reynolds number is reached, and then jumps abruptly to perhaps twice as much before it begins to decrease again as R increases.

### References

H. Lamb, Hydrodynamics, 6th ed. (New York: Dover, 1945). The classic treatise on hydrodynamics, first published in 1879, and still an excellent reference. It uses components exclusively. Viscosity is treated in Chapter XI.

R. L. Daugherty and J. B. Franzini, Fluid Mechanics with Engineering Applications 6th ed. (New York: McGraw-Hill, 1965). The first edition of Prof. Daugherty's text was in 1916. A good reference for engineering applications, especially pipe flow, which is treated in detail. Theoretical hydrodynamics is avoided in this undergraduate text, but Poiseuille's problem is solved.

L. Page, Introduction to Theoretical Physics, 3rd ed. (New York: D. Van Nostrand, 1952). Page has a concise treatment of elasticity and hydrodynamics that makes an excellent review. Dyadics are used.