Lagniappe: Compressible Flow

Throughout the entirety of my commentary before, I dealt solely with fluids whose density was always constant. Although this is usually an excellent approximation, there are in fact situations where the change of a fluid’s density is not only important but useful—in these situations, we unsurprisingly label the flow as compressible.

In Part II, I had mentioned that for fluids in general, a fluid’s density was a function of pressure: \rho = \rho(p). The reason we’re able to make the assumption that most fluids are incompressible is by guessing (and hoping) that this function is relatively “flat”; that is, that large changes in the pressure generate small changes in the density, small enough to ignore. If we wanted to fully understand compressible flow, we would need a full understanding of the function \rho(p), which is more in the realm of thermodynamics than fluid mechanics. But for many purposes, we only need to understand how the density of a fluid changes when its pressure changes—and if the range of fluid densities in a given scenario we’re studying is sufficiently small, then that density/pressure change ratio \frac{\partial\rho(p)}{\partial p} is essentially constant and only a property of the fluid in question.

It may seem like we need to overhaul the entire conceptual infrastructure I built up before to analyze compressible flow, but it’s actually surprisingly simple to handle! Consider the construction we were handling in Part V, where we used a control volume analysis inside of a pipe with a changing cross-section to study the key equations that governed steady flow within such a pipe; the fundamental equations of hydraulics. If I make the same assumptions I made in Part IV—the flow is steady, velocity always points directly into the pipe, etc.— but now consider that the densities at the inlet and outlet may be uniform, different values, control volume analysis dictates that the amount of mass flowing into the pipe needs to be the same as the amount of mass flowing out of the pipe. Mathematically, this translates into:

\displaystyle \int_{A_1}\rho_1 v_1\ dA\ = \int_{A_2}\rho_2 v_2\ dA

If I integrate this and only focus on average quantities as I did in Part IV, we find something suspiciously similar to the first fundamental law of hydraulics:

\displaystyle \rho_1 \langle v_{1} \rangle A_{1} = \rho_2 \langle v_{2} \rangle A_{2}

This naturally simplifies into the first fundamental law when the densities are the same.

A steady compressible flow through a pipe section with non-constant cross-sectional area. For the mass of fluid within the pipe to remain the same, the amount of mass flowing in needs to equal the amount of mass flowing out. For equal densities, the equation becomes identical to its incompressible counterpart.

Similarly, a control volume analysis on the momentum within such a pipe returns what is essentially the second fundamental law of hydraulics we derived in Part V, except with the inlet and outlet densities made distinct:

\displaystyle |\vec{H}| = \frac{\rho_1 + \rho_2}{2}|\vec{g}|\frac{A_1 + A_2}{2} \left(z_{1} - z_{2}\right) + \langle p_{1} \rangle A_1 - \langle p_{2} \rangle A_2 + \rho_1 \langle v_{1} \rangle^2 A_1 - \rho_2 \langle v_{2} \rangle^2 A_2

The density in the gravitational term is now the average of the densities at the inlet/outlet, and the momentum flux contributions into/out of the pipe now depend on their respective densities. So far, so good! But I’ll quickly mention that there’s another convenient but essentially equivalent formulation of these equations that will prove very useful—consider what happens if I take the compressible conservation of mass equation, move everything to one side, and divide both sides by the length L of the pipe in question:

\displaystyle \frac{\rho_2 \langle v_{2} \rangle A_{2} - \rho_1 \langle v_{1} \rangle A_{1}}{L}

In the limit where the pipe length becomes very, very small, the term in the above equation behaves like a derivative! And so by defining a “pipe direction” x, we can alternatively describe conservation of mass/the first fundamental equation through such a pipe using the following differential equation:

\displaystyle \frac{d\left(\rho \langle v \rangle A\right)}{dx} = 0

The same process generates the following differential description of the second fundamental equation of hydraulics:

\displaystyle -\frac{dH}{dx} = \frac{d\left(pA + \rho \langle v \rangle^2 A + \rho g Az\right)}{dx}

So now that we’ve gotten all that sorted out, what the heck can we actually do with these formulations of compressible hydraulics? Well, here’s a question we can tackle—what flow (if any) does a small change in fluid density generate?

To start, let’s consider steady flow within a horizontal pipe with a constant cross-sectional area, where the pressure and densities of the fluid in the inlet/outlet are uniform, unchanging and only slightly different. We can also make the pipe section very short—short enough to make the momentum losses within it negligible. This rapidly simplifies a lot of the content of the second fundamental equation:

\displaystyle 0 = p_{1} -  p_{2} + \rho_1 \langle v_{1} \rangle^2 - \rho_2 \langle v_{2} \rangle^2

We can then ask ourselves what the flow speed through such a pipe must be in order for the small differences in density to be sustained steadily. To do so, we can hypothesize that the speed throughout the pipe section is a uniform constant value a—this is a good approximation, since compressible conservation of mass predicts small speed changes for small density changes—and see how such a speed would depend on the changes in pressure/density.

A pipe of constant cross-section with slight differences in density & pressure at either end. Hypothetically, these density/pressure differences will cause a flow through the pipe of constant speed, which is related to the relationship between the change of pressure in the pipe and the change of density in the pipe.

Inserting this hypothesis into the second fundamental equation gets us:

\displaystyle 0 = p_{1} - p_{2} + \rho_1 a^2 - \rho_2 a^2

\displaystyle a^2 = -\frac{p_{2} - p_{1}}{\rho_2 - \rho_1}

Because we had said that the differences in pressure and density were very small, we find a familiar term appear when we make those differences very small:

\displaystyle a^2 = \lim_{\Delta p, \Delta \rho \to\ 0} -\frac{p_{2} - p_{1}}{\rho_2 - \rho_1} = -\frac{\partial p}{\partial \rho} = -\left(\frac{\partial \rho}{\partial p}\right)^{-1}

This is an interesting, nonintuitive result; if the density of a fluid changes somewhere a little bit, a flow towards the lower density direction will be generated with speed a = \sqrt{\frac{\partial \rho}{\partial p}^{-1}} that only depends on the way a specific fluid generates pressure differences as a result of internal density differences. This particular type of flow is commonly called sound, and the flow speed a is a property of a given fluid called its speed of sound. And because we had stated most fluids generate small density changes for large pressure differences \left(\frac{\partial \rho}{\partial p} << 1\right), we should expect speeds of sound a = \sqrt{\frac{\partial \rho}{\partial p}^{-1}} in general to be pretty big; 343 meters per second in air, and 1481 meters per second in water. That’s pretty fast!

There’s another question we can ask about compressible flow—how does it behave inside of horizontal pipes where the cross-sectional area is not constant? We can start by thinking about what happens for a fully incompressible flow inside such a pipe; as the first fundamental equation of hydraulics indicates, a decrease in cross-sectional area is directly connected to an increase in speed, and the second fundamental equation of hydraulics indicates that that speed increase must lead to a drop in pressure. One can imagine that if one takes this to the extreme, the drop in pressure becomes large enough that the density of the fluid begins to change. So what happens when that occurs? (Be warned: what follows will be an absolute mess of algebra.)

Let’s start by looking at the differential form of the first fundamental equation for compressible flow:

\displaystyle \frac{d\left(\rho \langle v \rangle A\right)}{dx} = 0

Expanding it out, we can spot the term introduced by compressibility immediately—the density change term.

\displaystyle \frac{d\rho}{dx}\langle v \rangle A + \rho\frac{d\langle v \rangle}{dx} A + \rho\langle v \rangle \frac{dA}{dx} = 0

Because we know that changes in the density are introduced by changes in pressure, we can introduce the speed of sound into the mix:

\displaystyle \frac{d\rho}{dx}= \frac{d\rho}{dp}\frac{dp}{dx} =\frac{1}{a^2}\frac{dp}{dx}

\displaystyle \frac{1}{a^2}\frac{dp}{dx}\langle v \rangle A + \rho\frac{d\langle v \rangle}{dx} A + \rho\langle v \rangle \frac{dA}{dx} = 0

Now we need to handle the pressure change term. Luckily, we have the differential form of the second fundamental equation to help us out! Applying it to this system, getting rid of gravitational terms since the pipe is horizontal (reasonable), and neglecting momentum loss (somewhat unreasonable, but oh well), we obtain the following:

\displaystyle 0 = \frac{d\left(pA + \rho \langle v \rangle^2 A\right)}{dx}

Expanding out, and noting that \displaystyle \frac{d\left(\rho \langle v \rangle A\right)}{dx} = 0, we get:

\displaystyle 0 = \frac{dp}{dx}A + p\frac{dA}{dx} + \rho \langle v \rangle A \frac{d\langle v \rangle}{dx}

Isolating the pressure derivative term:

\displaystyle \frac{dp}{dx} = - \frac{p}{A}\frac{dA}{dx} - \rho \langle v \rangle \frac{d\langle v \rangle}{dx}

Inserting this big mess into that first equation with the speed of sound in it that we were dealing with earlier leads to:

\displaystyle -\frac{1}{a^2}\left(\frac{p}{A}\frac{dA}{dx} + \rho \langle v \rangle \frac{d\langle v \rangle}{dx}\right)\langle v \rangle A + \rho\frac{d\langle v \rangle}{dx} A + \rho\langle v \rangle \frac{dA}{dx} = 0

Dividing everywhere by \displaystyle \rho \langle v \rangle A and grouping everything that’s multiplied by either the speed derivative or area derivative results in:

\displaystyle \frac{1}{\langle v \rangle} \left(1 - \frac{\langle v \rangle^2}{a^2}\right)\frac{d\langle v \rangle}{dx} + A \left(1 - \frac{p}{\rho a^2}\right) = 0

Look at that—we managed to group everything that involves the pressure and the density into one term! And luckily for us, that term is usually small enough to neglect—we had already claimed that the speed of sound is usually very large, so it stands to reason that pressures need to be very high in order for that term to matter. Neglecting that term, and solving for the area derivative, we find:

\displaystyle \frac{dA}{ds} = \frac{A}{U}\left(\frac{\langle v \rangle^2}{a^2} - 1\right)\frac{d\langle v \rangle}{dx} = 0

Alright, that’s enough algebra for today—now let’s try to pull some useful facts out of all this baloney.

Consider what would happen if I was a rocketry engineer and was trying to design a rocket nozzle—because I know that ejecting fluid at high speeds gives me more thrust, I want to try and design a nozzle that speeds up the fluid as much as possible within it. For a completely incompressible fluid, the obvious answer is to make the nozzle decrease in cross-sectional area, as this leads to a proportional speed increase by conservation of mass. You can also see this in the above equation when the flow speed is sufficiently low; area and speed are both positive quantities, and \left(\frac{\langle v \rangle^2}{a^2} - 1\right) is negative for speeds below the speed of sound, so the area derivative needs to be negative for the speed derivative to be positive. But as the flow keeps speeding up, something curious happens; once the flow gets faster than the speed of sound, the sign of the \left(\frac{\langle v \rangle^2}{a^2} - 1\right) term flips and becomes positive! Consequently, the area derivative needs to be positive for the speed derivative to be positive, which is the exact opposite of what was going on below the speed of sound. For speeds faster than the speed of sound, the nozzle’s cross-sectional area needs to start getting bigger in order to accelerate the fluid more.

As a result of all of this, your rocket nozzle design would look like a pipe that gets smaller in radius until the point you predict the flow hits the speed of sound, and then the pipe would start getting bigger again up to some radius based on your operating design. This is precisely why basically all rocket nozzles have that quasi-hourglass shape!

A schematic of a rocket nozzle. In the first section, the flow speed is below the speed of sound, and incompressible physics dominates—the cross-sectional area decreases to increase the speed. After the flow speed becomes larger than the speed of sound, compressible physics dominates—the cross-sectional area increases to further increase the flow speed. The total increase in flow speed is very large, causing a large imbalance in inlet/outlet momentum flux, which leads to a large thrust.

Physically, what’s happening is that two effects compete to speed up/slow down a compressible fluid—the effect you see in incompressible flow that causes flow to speed up/slow down when the cross-sectional area decreases/increases, and the slow down/speed up caused by the fluid density increasing/decreasing as a result of “squeezing”/”spreading out” when the cross-sectional area of the pipe decreases/increases. When the flow speed is below the speed of sound—commonly called subsonic flow—this first effect wins out. When the flow speed is above the speed of sound—commonly called supersonic flow—the second effect dominates. And the threshold between them can be quantified by the dimensionless number that is the ratio of the flow speed to the speed of sound of that fluid, commonly called the Mach number; when the Mach number is above 1, the flow is supersonic.

Part XIII: Boundary Layers, Lift & Drag

Now that we’ve said everything we can about flows far from immersed objects and flows near immersed objects, it makes sense to try and understand the flow everywhere in between. Even though we won’t be able to make too many quantitative statements about it—the math is way too complicated in that region—we will be able to make a couple of general statements about it.

The term boundary layer is often used as a catch-all term for any region between the surface of the object and the background flow at infinity, but here we use it specifically to refer to the regions in the fluid that are neither described by near or far field flow assumptions. Most of the time, this region is illustrated or described as a velocity field that monotonically increases in magnitude from 0 at the surface of the object to the background speed. Such a description is both shockingly common and always incorrect.

Although this can be proved rigorously, it suffices to look at the form of far field flow we derived two Parts ago; the velocity in the far field increases as we get closer to the object. As a result, the flow speed can’t just increase slowly until it gets to the background flow speed as we move away from the object. It needs to increase in the near field, and then keep increasing until it gets to a speed that’s larger than the background speed somewhere within the boundary layer, and then decay into the background speed in a way consistent with our expressions for far field flow. There isn’t a consistent name for this, but I like to call it the speed bump.

Two sketches of velocity profiles in the boundary layer. The common representation on the left assumes that the flow speed is always something between 0 (at the object surface) and the background flow speed U (at infinity). The more accurate representation on the right demonstrates the flow speed “overshooting” the background flow speed and then decaying into the background flow speed U at infinity, consistent with our findings for far field flow. This is the so-called “speed bump”.

Another qualitative feature of boundary layer flow that plays an important role in fluid dynamics is that it can reverse directions. In such a scenario, the near field flow around the object is pointing in the opposite direction as the far field flow associated with the object/flow, and so the flow needs to “turn around” somewhere in the boundary layer—this was the case for the sphere with a vortex pair behind it that we saw in the last Part, for example. Because the velocity changes continuously, as the flow “turns around”, there has to be some point at which the flow velocity becomes zero when moving away from the object. Around these points, the fluid effectively behaves as if it were near the surface of an imaginary object, and so any fluid bounded by the object and these imaginary surfaces effectively becomes locked in within it. This phenomenon is called flow separation, as the fluid flow appears to “separate” off of the object, and the flow bounded within the zero-velocity surfaces is said to be recirculating flow. Unsurprisingly, the region of fluid bounded by these surfaces is called the recirculation zone.

A sketch of the process of flow separation. Flow near the object is pointing in opposite directions, indicating that there is a region (the dotted line) where the flow velocity is zero within the fluid. Flow bounded from above by the dotted line lies within the recirculation zone, and fluid within this zone cannot leave it.

This is by-and-large all we can generically say about steady flow in the boundary layer. That being said, we’re going to find that—just as in pipe flow—flow in the boundary layer always becomes turbulent and unsteady once certain conditions are met, and we can find those conditions using the same technique we did for pipe flow; dimensional analysis.

To make things a little easier, let’s take a look at our tried-and-true scenario of a fixed rigid object immersed in a fluid with density \rho and viscosity \mu, with a background speed U and background pressure p_\infty. But this time, we’ll specifically make the immersed object a sphere with radius R. These parameters define everything about the flow around the sphere.

Steady flow past a sphere of radius R. The fluid has a uniform density \rho and viscosity \mu. At infinity far from the object, the flow is uniform with speed U and pressure p_\infty. This fully defines the flow everywhere.

As we saw for pipe flow, the Navier-Stokes equation (and by extension the physics of this problem) only involves pressure gradients, not the absolute pressure. And because the total pressure in the system is the sum of an absolute background pressure that is uniform everywhere and a perturbation pressure that isn’t, none of the pressure gradients are affected by the value of the background pressure! Consequently, the system is independent of the background pressure as long as the pressure isn’t so low that the liquid boils.

That leaves us with 4 dimensional parameters that define the flow based on 3 units of measurement—length, mass, and time. If we tried to construct a turbulence index \zeta as a function of these four parameters, we’d find that (thanks to the Buckingham Pi theorem) the turbulence index can only be a function of a single dimensionless number—the Reynolds number Re.

\displaystyle \zeta = f\left(Re\right) = f\left(\frac{\rho U R}{\mu}\right)

Just like we saw in pipes, whether or not the flow past an immersed sphere is steady or turbulent is entirely decided by the Reynolds number. But the transition to turbulence in flow past spheres (and immersed objects in general, for that matter) has a very rich structure that doesn’t just go from “nice and clean” to “chaotic and unpredictable”. This structure is largely present in all flow transitions for flows past immersed objects.

For steady flow past a sphere, the flow takes one of two forms—flow with no recirculation zone at the lowest Reynolds numbers, or flow with a small recirculation zone in the back consisting of two counter-rotating vortices at the small-but-not-smallest Reynolds numbers.

Two forms of steady flow past an immersed sphere. For the smallest Reynolds numbers, the flow does not recirculate, and the flow lines are symmetric in the “front” and “back” of the sphere. As the Reynolds number increases, a recirculation zone (or wake) forms behind the sphere, composed of two counter-rotating vortices.

When the Reynolds number increases enough, we find that the flow fails to be steady, but in a peculiar way—the vortices in the recirculation zone begin to wiggle perpendicular to the flow, and “detach” from the object, forming a streak of vortices called a von Karman vortex street. These emitted vortices travel downstream many distances longer than the radius of the sphere, and are the beginning of what we might recognize in common parlance as a wake. Although this flow is turbulent, it is quite simple and structured, and can be studied theoretically (although we will not do so here).

Flow past an immersed sphere at Reynolds numbers where flow begins to become unsteady. This unsteadiness begins to manifest as a von Karman vortex street, which is an unsteady, periodic downstream “emission” of counter-rotating vortices off of the back of the sphere. This forms a much larger recirculation zone than in the steady flow case.

As the Reynolds number increases further, the structure of the vortices being emitted from the back of the object decompose into a turbulent chaotic mess, and the flow appears to look well-behaved everywhere except for in a long messy “tail” immediately behind the object—this is a proper wake. Finally, increasing the Reynolds number even more results in the flow becoming turbulent everywhere in the boundary layer, not just in the wake. In this final stage, the flow has become fully turbulent and behaves very much as it does in turbulent pipe flow, where it possesses very little structure or order anywhere. These last two stages of turbulence are commonly referred to as subcritical and supercritical, respectively.

Subcritical (left) and supercritical (right) turbulent flow past an immersed sphere. In subcritical flow, the flow behind the sphere possesses a large recirculation zone (or wake) filled with chaotic unsteady flow, but the flow outside of this zone is steady. In supercritical turbulent flow, the flow everywhere in the boundary layer is turbulent.

Given that dimensional analysis gave us a way to at least classify observed flows past an immersed sphere, perhaps it will be useful to answer another extremely important question about flow past immersed objects; the forces that fluids exert on them.

Usually, the forces on an immersed object are split up into components that are parallel and perpendicular to the background flow direction—the parallel component is usually called the drag force, and the perpendicular component is called lift force. Drag and lift are both just different components of the same force, caused by both pressure-induced tractions and flow-induced molecular interaction tractions on the surface of the object. The general equation for the force on an immersed object is then:

\displaystyle \vec{F} = \int_S \vec{t}_s\ dA= \int_S -p \hat{n} + \boldsymbol \tau \cdot \hat{n}\ dA

\displaystyle \vec{F} = \int_S -p_\epsilon \hat{n} + \mu \left(\nabla \vec{v}_\epsilon + \nabla \vec{v}_\epsilon^T\right) \cdot \hat{n}\ dA

Notice the similarity to the expression we found in Part II for buoyant forces on a submerged object; this formula is a generalized version of it accounting for fluid flow.

However, we can’t solve for the flow in the boundary layer, and so we don’t really have a chance to derive anything for the net force acting on the object and its lift/drag components using this expression. But we can use dimensional analysis, coupled with some situational assumptions, to get surprisingly general expressions.

For example, we can try to find how many independent things with units of force we can construct out of the set of 4 dimensional parameters we have available for flow past a sphere. It turns out we can only make two, \mu U R and \rho U^2 R^2. This means that either the lift or the drag force on a sphere can be represented in the following form:

\displaystyle F_d\; \text{or}\; F_l = \sum\limits_{n \in \mathbb{R}} \alpha_n \left(\rho U^2 R^2 \right)^n   \left(\mu U R \right)^{1-n}

This is correct, but essentially useless; we need some other assumptions to get a workable, specific expression. One assumption we could make is that the surface tractions coming from the viscous effects are much smaller than those coming from the pressure perturbations; mathematically, that \mu |\left(\nabla \vec{v} + \nabla \vec{v}^T\right)| \cdot \hat{n} << |p| \hat{n} on the surface of the sphere. This is true if the viscosity is very small, or if the average velocity gradient on the surface of the sphere is on average zero because the flow in the boundary layer is chaotically flipping direction. In either scenario, the Reynolds number would be very large.

In this scenario, the drag/lift force is independent of the viscosity, and so anything dependent of it must have no bearing on the drag/lift force—namely, the \mu U R term. As a result, for high Reynolds numbers, we expect the drag and lift to have the following form:

\displaystyle F_d\; \text{or}\; F_l = \alpha_1\rho U^2 R^2

where \alpha_1 is a numerical constant independent of any of the other physical parameters in the system, a consequence of the dimensional crisis caused by removing the viscosity. This kind of drag, which is quadratic in the velocity, is usually called Newtonian drag.

Alternatively, we can consider what happens when the viscosity is very large, or when the Reynolds number is very low. In such a scenario, the terms dependent on the viscosity in the above sum should dominate, and anything independent of the viscosity shouldn’t influence the drag force; particularly, the \rho U^2 R^2 term. Removing the only parameter that solely appears in this expression, the density, we find the following expression for the drag/lift force at low Reynolds numbers:

\displaystyle F_d\; \text{or}\; F_l = \alpha_0 \mu U R

where \alpha_0 is also just a number. This kind of drag is called Stokes drag, and it is linear in the velocity of the object/background flow. If one attempts to solve for Stokes flow over a sphere assuming the near field flow equations apply everywhere (which is incorrect), one would find that \alpha_0 = 6 \pi for the drag force relationship.

Usually, fluid mechanicians plot the drag (or lift) force in terms of a nondimensional number called a drag (or lift) coefficient C_d\; \text{or}\; C_l, which is the drag force divided by \rho U^2 R^2 or some multiple of it:

C_d\; \text{or}\; C_l =  \frac{F_d\; \text{or}\; F_l}{\rho U^2 R^2}

This nondimensionalization of the drag/lift force is convenient because this nondimensional number can only be a function of nondimensional parameters constructed from the 4 physical variables that fully describe the physical set-up of flow past a sphere. And there’s only one nondimensional parameter we can make with those four variables; the Reynolds number.

If you were to plot the drag coefficient versus the Reynolds number Re for a sphere based on experimental results, this is what you’d find:

The drag coefficient as a function of the Reynolds number for a sphere, as determined by experiment. At low Reynolds numbers, the relationship between the drag force and the speed is linear, resulting in an inversely proportional relationship between the drag coefficient and the Reynolds number. At high Reynolds number, the drag force is quadratic in the speed, leading to the drag coefficient being independent of the Reynolds number at high Reynolds numbers.

This is all consistent with what we expected from dimensional analysis! Now we might ask, how do any of these results change when the object isn’t a sphere, but something totally different?

A simple way to illustrate this might be by considering what would happen if I put a spherical “nose” on the sphere with a radius R_2, directly in the front of the sphere. This parameter R_2 defines everything about the nose, and R along with R_2 fully define the geometry of this new class of immersed object.

Flow past an immersed object consisting of a sphere of radius R with a “nose” of radius R_2. The flow in this system is uniquely defined by the same parameters as for the sphere, with the addition of the “nose” radius.

If I did the whole dimensional analysis rigmarole I did for the sphere for this new object, I’d find the same things except for the presence of another, new dimensionless parameter: \frac{R_2}{R}. More importantly, if I tried to find expressions for the drag/lift force in the limits of low and high Reynolds numbers, I’d get stuck with a bunch of horrid sums again:

\displaystyle F_d\; \text{or}\; F_l\; \text{for low Re}\; = \sum\limits_{n \in  \mathbb{R}} \alpha_n \left(\mu U R \right)^n \left(\mu U R_2  \right)^{1-n}

\displaystyle F_d\; \text{or}\; F_l\; \text{for high Re}\; = \sum\limits_{n \in  \mathbb{R}} \alpha_n \left(\rho U^2 R^2 \right)^n   \left(\rho U^2 R_2^2   \right)^{1-n}

However, in either case, the \alpha‘s can only be a function of the only dimensionless objects we can construct; the radii ratio \frac{R_2}{R}. Because of that, we can actually factor out a couple of things out of the sum! Simplfying a bit, we’ll find that everything that isn’t a function of either radius winds up on the outside:

\displaystyle F_d\; \text{or}\; F_l\; \text{for low Re}\; =  \mu U \sum\limits_{n \in  \mathbb{R}} \alpha_n R^n R_2^{1-n}

\displaystyle F_d\; \text{or}\; F_l\; \text{for high Re}\; =  \rho U^2 \sum\limits_{n \in  \mathbb{R}} \alpha_n R^{2^n} R_2^{2^{1-n}}

Everything that’s within the sums is only a function of the geometry of the object, and so can be grouped into a geometric factor \Gamma which we choose to have units of length. With this, we find that:

\displaystyle F_d\; \text{or}\; F_l\; \text{for low Re}\; =  \mu U \Gamma_{\text{low}}

\displaystyle F_d\; \text{or}\; F_l\; \text{for high Re}\; =  \rho U^2  \Gamma_{\text{high}}^2

The great thing about this is that, if we were to add another geometric feature to the immersed object (say, “ears”), we’d wind up finding the exact same thing, albeit with a presumably more complicated geometric factor involving more geometric parameters. As a result, one can construct whatever kind of immersed object one wants using features, each defined by a single length, and the above results will still hold! For example, this indicates that at the very low and very high Reynolds number limits, the drag/lift force dependence on the viscosity, density, and speed are independent of the shape of the object. Anywhere in between, the geometry of the object does affect those relationships by affecting the values of the \alpha_n‘s.

Things to Think About

1. What do you think affects the magnitude of the speed bump, and why? Do you expect to be bigger at lower Reynolds numbers or higher ones?

2. What physical phenomena do you think affect flow separation? How could you encourage or ameliorate it?

3. Do you think that altering the geometry of an object can remove some of the “stages” of flow between non-recirculating and fully turbulent? How or why not?

4. Why do you think the drag coefficient for a sphere sharply drops after the flow becomes supercritically turbulent?

5. Approximating a swimmer as a funny-looking sphere, does a swimmer experience Stokes drag, Newtonian drag, or something in between?

Part XII: Near Field & Creeping Flow

Using the mathematical tools of asymptotics and perturbation theory, we’ve managed to find what the flow past rigid stationary objects—or the flow caused by rigid moving objects—looks like far away from those objects, in the region called the far field. Now we’d like to try and see if we can pull the same tricks to get the shape of flow in the region immediately next to an object; a region perhaps unsurprisingly referred to as the near field.

A rigid unmoving object is immersed in a uniform flow of constant pressure. We are currently interested in obtaining expressions for the flow in the region near the object, here surrounded by a transparent shell, called the near field. A spherical coordinate system has been illustrated for convenience.

One of the key reasons fluid mechanicists are interested in solving flows in the near field is to calculate forces on immersed objects. Just as we found in Part II, we’ll find that the force acting on an immersed object results from integrating the tractions \vec{t}_s on the surface of the object. The difference is that now we have additional contributions to the surface traction stemming from the flow-induced molecular forces in the fluid:

\displaystyle \vec{F} = \int_S \vec{t}_s\ dA = \int_S \left(-p \hat{n} + \boldsymbol \tau \cdot \hat{n}\right)\ dA

Because these tractions are being evaluated at the surface of the object, only the characteristics of the flow in the region immediately next to the object’s surface are going to influence the force on the object. Consequently, forces on objects are entirely determined by near field flow, and so we have good reason to try and characterize it.

Since we’re dealing with the same physical scenario as before, that of a steady uniform flow over a still rigid object (and its moving object counterpart, by extension), we can try to use the same perturbation flow approach we used before and expand the velocities and pressures into background and perturbation parts. This time, however, the background velocity we’re perturbing off of isn’t the velocity at infinity—we don’t care about what’s happening at infinity—it’s the velocity on the surface of the object.

Luckily, we have a very robust requirement for the velocity field near the object already. Our definition of fluid velocity we came up with in the very first Part, as the weighted average of the molecules in a tiny “molecular “net”, demands that the velocity transition continuously from the velocity of the object as we move away from it. This is called the no-slip condition, and for a still object, it means that \vec{v} = 0 on the surface of the immersed object. As a result, we can take our background velocity to be no velocity at all; \vec{v}_b = 0.

For the background pressure p_b, we don’t really have an obvious choice. The pressure assuredly changes over the surface of the object, so there isn’t a single value of the pressure at the surface that we can perturb off of. Consequently, we can simply take the background pressure to be the pressure of the flow at infinity, applied everywhere, since it’s the only pressure the physical set-up we’re considering is “giving” us. This has a considerable advantage in the sense that the net force across the surface of the object caused by the background pressure is identically zero, since the background pressure is the same everywhere on the surface of the object and integrating a surface normal over a closed surface is identically zero:

\displaystyle \vec{F}_{p_b} =  \int_S -p_b \hat{n}\ dA =  -p_b \int_S \hat{n}\ dA = 0

This also has the nice perk that \nabla p_b = 0, so it falls out of the perturbed Navier-Stokes equations.

Unfortunately, using the exact same physical set-up we used before for far field flow and these background quantities, we just wind up getting the original Navier-Stokes + conservation of mass system of equations for the perturbation flow:

\displaystyle \rho \nabla \vec{v}_\epsilon \cdot \left(\vec{v_b} + \vec{v_\epsilon}\right) = - \nabla p_\epsilon\ +\    \mu \nabla^2 \vec{v_\epsilon}

\displaystyle \nabla \cdot \vec{v}_\epsilon = 0

That certainly didn’t simplify things, but hopefully the following steps will.

We’d like to introduce a kinematic constraint just as we did for the far field case, but this time our constraint should be based on the behavior of the flow immediately next to the object, rather than the behavior far away from it. Because of the no-slip condition, the entire velocity field necessarily drops to zero as we move closer and closer to the object, which in this case is just the perturbation velocity field. This means that we can always define a region in which \rho\nabla\vec{v}\cdot\vec{v} is much smaller than \mu\nabla^2\vec{v}, no matter what the values of \nabla\vec{v} ,  \nabla^2\vec{v}, \rho and \mu are. This is equivalent to saying that the Reynolds tensor is effectively zero in this region, and it is precisely this region that we’ll consider as the near field.

With this in mind, it seems pretty obvious to take as our kinematic constraint that \nabla\vec{v}\cdot\vec{v} = 0. If we do that, we’ll obtain the following:

\displaystyle 0 = - \nabla p_\epsilon\ +\    \mu \nabla^2 \vec{v_\epsilon}

\displaystyle \nabla \cdot \vec{v}_\epsilon = 0

\nabla \vec{v}_\epsilon \cdot \vec{v_\epsilon} = 0

This form of the Navier-Stokes equation might look familiar—it’s exactly the form we had when we were solving for Hagen-Poiseuille flow. The difference is that now we don’t have the unidirectionality in the flow that we had when solving for pipe flow, so we can’t take this partial differential equation and turn it into an ordinary differential equation as we did before. Flows that we construct using this system of equations are usually called creeping flows or Stokes flows.

Doing a cute mathematical trick makes it really easy to solve for the pressure in a creeping flow. If you take the divergence of both sides of the Navier-Stokes equation for creeping flow, you get:

\displaystyle \nabla \cdot 0 = \nabla \cdot - \nabla p_\epsilon\ +\    \nabla \cdot \mu \nabla^2 \vec{v_\epsilon}

\displaystyle \nabla \cdot 0 =  \nabla \cdot - \nabla p_\epsilon\ +\    \mu \nabla^2 \left(\nabla \cdot \vec{v_\epsilon}\right)

\displaystyle  \nabla^2 p_\epsilon = 0

That last equation for the pressure is exactly the same equation that we had for the velocity potential in far field flow, Laplace’s equation! As a result, we know that the perturbation pressure has the same general expression as the velocity potential from far field flow:

\displaystyle p_\epsilon =  \sum \limits_{\ell=0}^{\infty}   \sum\limits_{m=-\ell}^{\ell}  \left(A_{\ell m}r^\ell + \frac{B_{\ell  m}}{r^{\ell+1}}\right)  Y^m_{\ell}(\theta,\varphi)

Using the general solution to Laplace’s equation in the form of spherical harmonics, you could then plug this solution for the perturbative pressure into the (simplified) Navier-Stokes equation, and you’d get your answer for the velocity field. Easier said than done, but I digress.

However, there’s a catch. The problem is that unlike irrotational flow, some Stokes flows don’t necessarily cause a perturbation in the pressure! To prove it, take a look at what happens if we a priori assume that the perturbative flow has a constant perturbative pressure (including the possibility that it’s zero), and so \nabla p_\epsilon = 0. With this assumption, we wind up finding:

\displaystyle 0 = - \nabla p_\epsilon\ +\    \mu \nabla^2 \vec{v_\epsilon}

\displaystyle 0 = \nabla^2 \vec{v_\epsilon}

\displaystyle \nabla^2 v_{\epsilon_r} = \nabla^2 v_{\epsilon_\theta} = \nabla^2 v_{\epsilon_\varphi} = 0

That means that, in the case of Stokes flows with uniform perturbative pressure, each component of the perturbation velocity satisfies Laplace’s equation, whose solutions have the exact same general forms as the solutions for the pressure—an infinite sum of solid harmonics. Most of the time, fluid dynamicists refer to each of these velocity fields as the homogeneous \vec{v}_h part for the uniform-pressure contribution and inhomogeneous \vec{v}_{\not h} part for the varying-pressure contribution of the near field flow:

\displaystyle \vec{v}_\epsilon = \vec{v}_h + \vec{v}_{\not h}

\displaystyle \nabla^2 \vec{v}_{h} = 0,\quad  \nabla^2 \vec{v}_{\not h} = \nabla p_\epsilon

p_\epsilon (r,\theta,\varphi) = p_h + p_{\not h}(r,\theta,\varphi)

Solving both equations for each of the components of the perturbation velocity nets you the following general form for creeping flow velocity fields, which is often called Lamb’s general solution. To skip you the trouble, I’ve written it below, where \vec{r} represents the position vector. Be warned; this is really ugly.

\displaystyle \vec{v}_\epsilon = \sum_{n=-\infty, n\neq1}^{n=\infty}  \sum\limits_{m=-n}^{n} \left[ \frac{(n+3)r^2\nabla\left[\left(A_{n m}r^n + \frac{B_{n  m}}{r^{n+1}}\right)  Y^m_{n}(\theta,\varphi)\right]}{2\mu(n+1)(2n+3)} - \frac{n\vec{r}  \left(A_{n m}r^n + \frac{B_{n  m}}{r^{n+1}}\right)  Y^m_{n}(\theta,\varphi)}{\mu(n+1)(2n+3)}\right] + \sum_{n=-\infty}^{n=\infty} \sum\limits_{m=-n}^{n} [\nabla \left(C_{n m}r^n + \frac{D_{n  m}}{r^{n+1}}\right) Y^m_{n}(\theta,\varphi)  + \nabla \times (\vec{r} \left(E_{n m}r^n + \frac{F_{n  m}}{r^{n+1}}\right) Y^m_{n}(\theta,\varphi) )]

Not only does this look like mathematical verborrhea, it also contains six distinct infinite sequences of coefficients. However, take a look at that \nabla \left(C_{n m}r^n + \frac{D_{n  m}}{r^{n+1}}\right) Y^m_{n}(\theta,\varphi) term inside of the second sum. This term, which is part of the homogeneous velocity field, has the same exact form as the general solution of far field flow in the previous Part; it’s the gradient of a velocity potential! As a result, the general solution for the flow velocity in the near field includes the general solution for the velocity in the far field. That doesn’t mean that we can accurately describe far field flows using the creeping flow equations, though; flows of this form always induce pressure variations in the far field thanks to Bernoulli’s equation, while the near field flows of this type never induce anything but uniform pressure changes everywhere by virtue of being part of the homogeneous velocity field.

In any case, we could try to take the general solution we got above and use the boundary conditions in the problem to restrict it into some relatively specific form we can derive conclusions about, like we did for far field flow. Unfortunately, we’ll find that’s not possible in general for two reasons.

The first is that the boundary condition associated with the object in the near field, that the velocity drop to zero on the surface of the object, doesn’t yield a massive restriction on the coefficients in the general solution of creeping flows. This is because our solution for the flow can still go to infinity at the origin if the origin is within the object, since the flow isn’t defined there anyways. That means that we can have terms that both grow and decay as we move closer to the origin, unlike what we found for far field flow (which just decays when moving closer to infinity). We can still use the condition that \vec{v} = 0 on the object’s surface to restrict the form of the flow, but it can only be done on a case-by-case basis as prescribed by the geometry of the object.

The second reason is associated with what near field flow should look like as we move away from the near field—intuition might tell you that the near field flow needs to eventually look like the imposed uniform flow as we move away from the object. But this isn’t necessarily true! Remember that our expressions & assumptions for near field flow are only valid in the region immediately next to the object, so expecting the near field flow to turn into the flow we prescribe at infinity is foolish in general. The near field flow should transition into something in a region somewhere between the near field and the far field, which will then turn into a dipole flow in the far field, which will then turn into a uniform flow at infinity.

That being said, sometimes constraining the near field flow by what it does at infinity isn’t a terrible approximation—fluid mechanicists do it all the time, and often come up with neat mathematical tricks to make the solutions constrained by the flow at infinity more accurate. However, doing so often leads to apparent paradoxes and nonsense results, so one needs to tread lightly in case one of these paradoxes pops up in your studies.

To give you an example, take a look at two different steady flows over a rigid unmoving cylinder, when one of these flows has a larger uniform velocity at infinity:

Sketches of two experimentally observed flows past a sphere; the scenario on the right has a larger uniform velocity at infinity. The flow on the left matches the Stokes flow prediction everywhere in the near field when using the boundary condition that the perturbative flow decay at infinity. The flow on the right does not, even in the near field region where the Stokes equations are always valid.

The solution one finds from the general form of creeping flows, combined with applying boundary conditions based on the uniform flow at infinity, matches the flow with no vortices almost identically everywhere in the flow. However, it doesn’t really seem to match the flow with the vortices; which would be fine, if the flows were dissimilar only after some distance away from the object. However, the flows don’t match even in the region we’d think they should, right next to the object where we said that flows are always creeping flows. How could this be happening?

The reason is that the flow in that region that isn’t quite the near or the far field is different, which means that the boundary conditions for the creeping flow in the near field change! So the flow immediately next to either sphere is still a creeping flow, they’re just different creeping flows because the boundary conditions for the creeping flow induced by the flow in the not-quite-near/not-quite-far region are different. I haven’t seen this mentioned in other fluid mechanics books, but I like to call it the matching paradox.

All of this should be motivation enough to dig into what flows look like in the not-quite-near, not-quite-far region. Formally, this region is often called the boundary layer, and it contains perhaps the most notable feature of flows past rigid immersed objects; wakes.

Things To Think About

1. How would you estimate some characteristic length representing how far the near field extends from the object?

2. Given the form of \boldsymbol \tau we found before, what does \boldsymbol \tau \cdot \hat{n} look like as a function of the velocity? What does it convey physically about the flow-induced forces on the object?

3. What other choices for the background pressure could you make, and how would they be better or worse than the choice we made?

4. How are the assumptions that we made to get the creeping flow equations here different from the assumptions we made to get the same equations in pipe flow?

5. What parts of the perturbation pressure induce a net force on an immersed objects?

6. When solving for creeping flow over a sphere, one can find a unique solution that both matches the boundary conditions both at infinity and at the surface of the object. Why is this solution only valid in the near field, given what we know about the flow behavior in the far field?

Part XI: Far Field & Irrotational Flow

Having said essentially everything we could about pipe flow, we now seek to characterize the last lingering question of basic fluid dynamics; flow past immersed objects.

In most descriptions of flow past immersed objects, and the only scenario we’ll cover here, we have some solid, rigid, stationary object surrounded by what can be approximated as an infinite amount of fluid. These types of flows tend to be called unbounded, or are said to be in an infinite domain, as they are defined over an infinite amount of space. That doesn’t necessarily mean we seek to characterize flow over all of space—just an infinite amount of it. All the way out at infinity, we’ll assume that the flow is a uniform flow, with a constant velocity in a single direction and a constant pressure; the presence of an embedded object then “bends” the flow close to it in a way consistent with the Navier-Stokes equation and with conservation of mass.

This physical set-up is convenient because it is equivalent to another ubiquitous scenario in fluid mechanics, that of a rigid object moving at a constant velocity through an infinitely large medium of otherwise stationary fluid. This scenario closely approximates the dynamics of a plane or bird in the sky, or of a submarine or fish in the ocean. The correspondence between this scenario and the first one I first mentioned comes from uniformly adding or subtracting the background flow; in this latter case, the flow will be a perturbation of a completely stationary flow by a rigid object moving with a velocity equal to the negative of the background velocity in the former case.

Two equivalent unbounded flows. On the left-hand side, a stationary object is placed in an infinite fluid with a uniform background flow and constant pressure at infinity. On the right-hand side, a rigid object is moving with a fixed speed and constant direction in a stationary background fluid with constant pressure at infinity. The velocity of the object in the right-hand scenario is equal and opposite to the velocity of the flow at infinity in the left-hand scenario. A spherical coordinate system at the origin has been illustrated for convenience.

With this in mind, one can consider the description of the velocity field/pressure field/flow to be a combination of the background flow and of some other velocity/pressure field which decays to nothing as you get further and further away from the object. Formally, this is called a flow perturbation, in the sense that the object “perturbs” the background flow only in a region close to the object. The velocity and pressure field would then be of the form:

\vec{v} = \vec{v}_b + \vec{v}_\epsilon

p = p_\infty + p_\epsilon

where the perturbation velocity field \vec{v}_\epsilon and perturbation pressure field p_\epsilon need to become negligibly small relative to the background velocity \vec{v}_b and pressure p_\infty respectively as we move further and further away from the immersed object.

With either set-up in mind, we could try to perform a mathematical analysis like we did for pipe flow and solve directly for the pressure and velocity everywhere. However, in this case, we don’t have nearly any of the simplifying assumptions we had with pipe flow; even though we can assume the flow is steady in the stationary immersed object case, the flow can potentially be in all three coordinate directions, and it can potentially change in all three directions. In addition, we know that turbulence is baked into the Navier-Stokes equations somewhere, meaning that even the steady assumption is bound to fail at some point. What can we do?

Well, here’s an idea; we could try analyzing the behavior of the flow only in a region “almost at infinity”. We can’t do it all the way at infinity, since all we’d get is the background flow by definition—but we can try to get a sense of what the Navier-Stokes equations and conservation of mass tell us about the perturbation flow as we move out further and further away from the object. Formally, we are attempting to understand the asymptotic behavior of the perturbation flow near infinity, in a region commonly referred to as the far field.

A still object is submerged in an unbounded background flow of uniform speed and constant pressure—a spherical shell around it has been constructed “almost at infinity”. We want to identify the characteristics of the flow in the question mark regions outside the shell, representing the region “almost at infinity” called the far field, using asymptotic methods. A spherical coordinate system has been drawn for convenience.

Solving the full Navier-Stokes equation plus conservation of mass in this scenario isn’t an option—not even the best mathematicians can do it—so what we’ll do is come up with an additional assumption about the perturbation velocity field far from the object that, when plugged in, makes the Navier-Stokes equation easier to solve. Such an assumption is commonly referred to as a kinematic constraint by dynamicists.

We want this kinematic constraint to have a couple of key features:

  1. It’s consistent with the perturbation velocity decaying to zero as we move further and further away from the object.
  2. It reduces the Navier-Stokes equation to something tractable.
  3. It doesn’t violate conservation of mass.

The idea is that now we’ll have three equations we can work with; and if we make the right kinematic constraint, we’ll be able to solve for the velocity or pressure straightforwardly and use the Navier-Stokes equation to calculate the fluid quantity we didn’t calculate before. Because we only care about the flow outside of the boundary of the spherical shell we drew above, the flows we obtain only need to be correct in the far field; the flows we obtain from this process are far field flows.

Let’s begin. By inserting the expression \vec{v} = \vec{v}_b + \vec{v}_\epsilon and p = p_\infty + p_\epsilon into the Navier-Stokes and incompressible conservation of mass equations for the scenario where we have a uniform background flow of constant speed and pressure, and under the assumption of no external forces and steady flow, we get:

\displaystyle \rho \nabla \vec{v}_\epsilon \cdot \left(\vec{v_b} + \vec{v_\epsilon}\right) = - \nabla p_\epsilon\ +\    \mu \nabla^2 \vec{v_\epsilon}

\displaystyle \nabla \cdot \vec{v}_\epsilon = 0

The criteria for the kinematic constraint on the velocity we listed above, namely the first one regarding the flow decay at infinity, induce a very specific behavior in the perturbation flow at the far field. Informally, the idea is that the flow far from the object should be unaffected by viscous effects—only the presence of the object is triggering viscous effects in the flow, and far enough away from it, those viscous effects should be trivial. This suggests that the ratio of convective to viscous effects (i.e. the Reynolds tensor) will be large in the far field, and so \nabla^2 \vec{v}_\epsilon \sim 0*. The \sim symbol indicates “asymptotic” equivalence, i.e. that the objects are identical as one approaches the asymptotic limit, which is at r \rightarrow \infty.

*I really wish I could prove this rigorously, but every strategy to do so that I could come up with is either pointlessly contrived or way too elaborate to show here. Oh well.

This creates a problem for our strategy of calculating far field flows; we know that our biggest obstacle to a straightforward solution is that nonlinear term, \displaystyle \rho \nabla \vec{v}_\epsilon \cdot \left(\vec{v_b} + \vec{v_\epsilon}\right), so it would make sense to try to get rid of it first. However, we know that in the far field, \nabla^2 \vec{v}_\epsilon \sim 0, so we can’t just get rid of it outright; it would make the Navier-Stokes equation just an equation for the pressure and not the velocity in the far field, which we don’t want. Somehow, we have to come up with a kinematic constraint that both simplifies the nonlinear convective term and completely gets rid of the viscous term. How are we going to do that?

The secret is in that alternative form of the Navier-Stokes equations I kept showing before—the Lamb form. If I plug in my expression for the velocity, I wind up with the following representation of the Navier-Stokes equation:

\displaystyle \rho\left(\frac{\nabla \left( \vec{v}_\epsilon \cdot \vec{v}_\epsilon \right)}{2} -  \vec{v}_\epsilon \times\left(\nabla\times\vec{v}_\epsilon \right) \right) =  -\nabla  p_\epsilon\ -\   \mu \nabla \times \left(\nabla \times \vec{v}_\epsilon \right)

This form makes a possible kinematic constraint obvious, one which simplifies the nonlinear terms and kills off the viscous term; the constraint that \nabla \times \vec{v}_\epsilon = 0. In general, the quantity \nabla \times \vec{v} is referred to as the vorticity \vec{\omega}, as it is proportional to the local rotation rate of the fluid. Small objects embedded in flows that satisfy this kinematic constraint do not rotate, and the flow itself is automatically unaffected by friction and rotation-related convective effects, as our Reynolds tensor calculation required for the far field.

If we take this as our constraint, we find the following system of equations:

\displaystyle \rho\left( \frac{\nabla \left( \vec{v} \cdot  \vec{v}\right)}{2}\right) =  -\nabla  p

\displaystyle \nabla \cdot \vec{v}_\epsilon = 0

\nabla \times \vec{v}_\epsilon = 0

Flows derived from this system of equations are alternatively referred to as either irrotational flows or potential flows. The latter naming convention comes from the fact that, thanks to the kinematic constraint, the velocity can be described as the gradient of a scalar quantity known as the velocity potential \Phi that satisfies the following equation:

\displaystyle \nabla^2 \Phi = 0,\quad \left(u = \nabla \Phi\right)

This equation is known as Laplace’s equation, and it is perhaps the easiest partial differential equation to solve (although solving it is still relatively complicated). In spherical coordinates, the general solution can be written as the sum of a bunch of different functions—to save you the trouble, here is the general solution to the Laplace equation in sum form:

\displaystyle \Phi = \sum \limits_{\ell=0}^{\infty}  \sum\limits_{m=-\ell}^{\ell} \left(A_{\ell m}r^\ell + \frac{B_{\ell m}}{r^{\ell+1}}\right) Y^m_{\ell}(\theta,\varphi)

where the A_{\ell m} and B_{\ell m} terms represent constant coefficients, and the Y^m_{\ell}(\theta,\varphi) terms are functions of only the angles called real spherical harmonics. As either of the indices in the sum (\ell or m) increases, these spherical harmonics become “wavier”. This type of sum expression is called a multipole expansion, as the waviness causes crests, or poles, in the resulting functions.

Taking the gradient of the velocity potential to find the velocity field, we find the following general expression for the components of all irrotational/potential flows:

\displaystyle v_{r_{\text{irrot.}}} =  \sum \limits_{\ell=0}^{\infty}  \sum\limits_{m=-\ell}^{\ell} \left(A_{\ell m} \ell r^{\ell-1} - \frac{B_{\ell m}\left(\ell + 1\right)}{r^{\ell+2}}\right) Y^m_{\ell}(\theta,\varphi)

\displaystyle v_{\theta_{\text{irrot.}}} =  \sum \limits_{\ell=0}^{\infty}  \sum\limits_{m=-\ell}^{\ell} \left(A_{\ell m}r^{\ell-1} + \frac{B_{\ell m}}{r^{\ell+2}}\right) \frac{\partial Y^m_{\ell}(\theta,\varphi)}{\partial \theta}

\displaystyle v_{\varphi_{\text{irrot.}}} =  \sum  \limits_{\ell=0}^{\infty}  \sum\limits_{m=-\ell}^{\ell} \left(A_{\ell  m}r^{\ell-1} + \frac{B_{\ell m}}{r^{\ell+2}}\right)\frac{1}{\sin{\theta}}\frac{\partial  Y^m_{\ell}(\theta,\varphi)}{\partial \varphi}

This all just seems like useless mathematical gobbledygook. But that’s because we haven’t restricted our form of the potential flow by a key principle—that the velocity decay in the far field. On inspection, this immediately indicates that every A_{\ell\geq1  m} term has to be zero, as the flow would fail to decay if they weren’t thanks to the r^{\ell-1} term. The flow coming from the A_{00} term is also identically zero, as Y^0_{0} is just a constant that is independent of either angle and the radial part of the flow is nullified by multiplication of the \ell = 0 factor. Finally, for a rigid object, B_{0m} = 0, as the flow generated by this term fails to satisfy conservation of mass in this scenario. This leaves us with:

\displaystyle v_{r_{\epsilon}} =  \sum \limits_{\ell=1}^{\infty}  \sum\limits_{m=-\ell}^{\ell} -\frac{B_{\ell m}\left(\ell + 1\right)}{r^{\ell+2}}Y^m_{\ell}(\theta,\varphi)

\displaystyle v_{\theta_{\epsilon}} =  \sum \limits_{\ell=1}^{\infty}  \sum\limits_{m=-\ell}^{\ell}\frac{B_{\ell m}}{r^{\ell+2}} \frac{\partial Y^m_{\ell}(\theta,\varphi)}{\partial \theta}

\displaystyle v_{\varphi_{\epsilon}} =  \sum  \limits_{\ell=1}^{\infty}  \sum\limits_{m=-\ell}^{\ell} \frac{B_{\ell m}}{r^{\ell+2} \sin{\theta}}\frac{\partial  Y^m_{\ell}(\theta,\varphi)}{\partial \varphi}

This still looks like a lot of terms. However, remember we only care about the behavior in the far field; and in the far field, no matter what nonzero values the B_{\ell m} coefficients have, the terms that decay the slowest will always eventually become much larger than any of the other terms. As a result, asymptotically, every far field flow always looks like the flow generated by the surviving terms with the lowest radial order—which for a steadily moving rigid object, are the B_{1m} terms. Even better, the spherical harmonics corresponding to Y^m_{1} are really just the same function rotated 90 degrees in either angle, so we can always rotate our coordinate system to get rid of the Y^{-1}_{1} and Y^{1}_{1} terms. This means we get the following universal form of far field perturbation flows for rigid immersed objects moving at a constant speed:

\displaystyle v_{r_{\epsilon}} \sim -\frac{2B_{1 0}}{r^{3}}\cos{\theta}

\displaystyle v_{\theta_{\epsilon}} \sim -\frac{B_{1 0}}{r^{3}}\sin{\theta}

\displaystyle v_{\varphi_{\epsilon}} \sim 0

Notice the velocity doesn’t depend on \varphi; this is because of the coordinate rotation trick I mentioned above.

I like to call this flow a dipole flow, in reference to the fact that the mathematical form of the flow field is identical to that of an electric dipole. It is important to mention that because of the analysis above, every flow perturbation caused by a steadily-moving rigid object, or by a rigid object embedded in a background flow, looks like a dipole flow in the far field. Here’s what that looks like in the latter case:

A far field flow corresponding to an object embedded in a uniform background flow. The flow is axisymmetric around the background flow axis, and is a sum of a uniform background flow and a dipole flow. The opaque grey sphere represents the near field, where the true flow behavior would be inconsistent with this flow. The background flow is deflected by the object about the background velocity axis, with the magnitude of the deflection dropping rapidly as one moves further from the object.

You might notice that the velocity vectors almost appear to be flowing over a sphere, as if the embedded object were a sphere itself. This is because every rigid object deflects flow in the far field as if it were a sphere—a mathematical consequence of the dipole flow dominating in the far field. As a result, every arbitrarily-shaped object has an effective radius r_{\text{eff}} such that a sphere with that radius deflects flow in the far field identically to the original object.

Switching over to the rigid object moving with constant velocity in a quiescent fluid case, here’s what that would look like:

A far field flow generated by a rigid object moving with a constant speed through a quiescent fluid. The flow is an axisymmetric dipole flow with an origin moving at a uniform speed. The fluid is pushed in the direction of the object’s motion in front of the object, pulled backwards in the opposite direction when near the object, and pulled forward again towards the object when behind it.

In the region immediately behind the object, the fluid is being pulled towards it; this phenomenon is called slipstreaming, and is exploited by nature and by humans to increase aerodynamic performance when traveling in groups. For example, this is why geese fly behind each other in V formations, or why cyclists and race car drivers like to drive immediately behind another car when they need a speed boost. The flow patterns of all objects moving steadily while immersed in quiescent fluids behave the same way far from the objects, so it is hopefully now unsurprising that the phenomenon is ubiquitous. They can only ever be different by the scaling factor B_{10}, which is called a dipole moment in electromagnetism contexts, and a quick fly-by dimensional analysis indicates that it must be proportional to the background/object speed U times the effective radius of the object cubed r_{\text{eff}}^3.

Now that we have our expression for the perturbation velocity field, we can use the (simplified) Navier-Stokes equation to get the expression for the pressure field in the static object case:

\displaystyle \rho\left( \frac{\nabla \left( \vec{v} \cdot  \vec{v}\right)}{2}\right) =  -\nabla  p

It might look like we’ll have to solve a partial differential equation, but notice that we’re taking the gradient of a scalar in all of the terms in this equation. As a result, we can collect all the terms inside of a single gradient operation and then “integrate it out”:

\displaystyle \nabla\left(p + \rho\frac{\vec{v} \cdot  \vec{v}}{2}\right) = 0

\displaystyle p + \rho\frac{\vec{v} \cdot  \vec{v}}{2} = C

where C is a constant.

This last equation is called Bernoulli’s equation, and provides a direct algebraic relationship between the pressure and the velocity of a fluid in the case where the viscous effects are negligible, as is the case (we claim) in the far field.

Because the equation holds at every point in the far field, including at infinity, we find that Bernoulli’s equation is satisfied even when the background flow is the only non-trivial flow:

\displaystyle p + \rho\frac{\vec{v} \cdot  \vec{v}}{2} =  p_\infty + \frac{\rho U^2}{2} = C

Noting that the total pressure is the sum of the background and perturbation pressures, and that the total velocity is the sum of the background and perturbation velocities, we can finally get an expression for the perturbation pressure:

\displaystyle p_\epsilon = \frac{B_{10} \rho  \left(2 r^3 U (3 \cos (2 \theta )+1)-B_{10} (3 \cos (2 \theta )+5)\right)}{4 r^6}

Although it may not be immediately obvious from the expression above, this perturbation pressure is always negative, which appears to be nonsense. However, notice that the total pressure, which is the thing we need to keep positive, isn’t necessarily always negative; what this result for the perturbation pressure indicates is that the effect of the object’s presence on the background pressure is to decrease the pressure, more and more as one gets closer to the object. Eventually, if the other assumptions we made for the flow don’t break down as we get closer to the object, the total pressure drops so much that it must become negative—indicating that our solution for the flow can’t be valid. This is yet another example of an existence failure for solutions of the Navier-Stokes equations, which we had seen before in pipe flow.

Practically, however, the liquid would boil due to the pressure reduction long before reaching that point. This occurs when the total pressure in the flow is equal to the vapor pressure p_{\text{vap}}. Using the equation for the perturbation pressure above, one can solve for the distance at which the liquid would vaporize as a function of the other system parameters:

r_{\text{vap}} = \sqrt[3]{\frac{\sqrt{B_{10}^2 \rho  \left(4 (3 \cos (2 \theta )+5) p_{\text{vap}}+\rho  (3 U \cos (2 \theta )+U)^2\right)}-B_{10} \rho  U (3 \cos (2 \theta )+1)}{4p_{\text{vap}}}}

For objects/background flows moving sufficiently fast, the pressure drop due to the presence of the object generates a vapor shell (also called a vapor cone) encasing the object, whose shape and location is described by the equation above. The equation for the vapor shell above gives us a lower estimate on what the “boundary” of the far field is, as the assumption we made when solving for far field flow clearly break down once we get inside the vapor shell due to the density change of the fluid. That being said, vapor shells rarely manifest except in the case of exceedingly rapid objects, in which case one needs to incorporate thermodynamics into the analysis above. In most cases, viscosity begins to become relevant long before hitting the vapor shell region.

The shape of the vapor shell predicted by the far field dipole pressure perturbation against a uniform background flow for one set of system parameters, while ignoring pressure-temperature coupling. Changing the parameters of the system causes the shape to change, but it is always roughly spheroidal. The vapor shell serves as an estimate of the boundary of the far field, as our flow assumptions fail within it.

This asymptotic analysis, using the irrotational kinematic constraint on the flow field, has netted us some very nice, universal results for flow far from an object embedded in a fluid. However, as the existence failure for the pressure near the object demonstrates, it’s not enough to get us some of the most important physical results we need to design objects in fluids for engineering purposes; namely, what the force on the object from the fluid flow is, and the nature of turbulence in flows over immersed objects. For this, we’ll need every trick in the book so far.

Things To Think About

1. To derive our kinematic constraint, we took the assumption that the viscous effects drop to zero in the far field. As a result, the Reynolds tensor should increase as one moves out into the far field. Is this correct, given the form of dipole flow?

2. Is the sum of two dipole flows also a potential flow? How would this be useful for calculating the flow caused by the motion of a group of distant objects, like airplanes?

3. What do you think would be the dominant behavior in the far field if we allowed the immersed object to expand/contract? You can make your life easier if you assume the expansion is really slow, so that the flow is still approximately steady.

4. Why can’t we use potential flow to describe behavior in the near field, i.e. very close to the object?

5. Instead of picking the irrotational kinematic constraint, we could have also simply stated \nabla^2 \vec{v} = 0s=3. Why didn’t I do that?

5. Instead of all this math, we could have tried doing a control volume analysis over some volume enclosing the object to get all the relevant information we needed. Why wouldn’t that be a good idea?

6. For the moving object case, consider a weather vane stuck on some fixed point in the far field while the object moves past it. What do you intuitively think happens to the weather vane as the objects moves past it, and how does the animation above help validate your intuition?

7. How do some of the other irrotational flows look like? Use the equations above to get a sense of what happens when you increase \ell.

Part X: Dimensional Analysis

We are, at the present, in a bit of a pickle. We have spent a lot of time deriving the laws of fluid mechanics to understand flow through a pipe, but we’ve found that even in something as simple as pipe flow, the inherent pathology of the laws of fluid mechanics causes pipe flow to spontaneously transition from something we can understand (Hagen-Poiseuille flow) to something we can’t (turbulence). We’ve also commented on the fact that experimentalists have noted that the transition from one to the other occurs when the Reynolds number, which we’ve usually taken to represent an approximation of the ratio of convective to viscous effects in a flow, passes a critical value—even though the actual ratio of convective to viscous effects in Hagen-Poiseuille flow is identically zero. How could this be?

To solve this mystery, we’re going to need apply a little bit of what may appear to be abstract nonsense. Let’s again take a look at a horizontal section of pipe of radius R, with an arbitrarily long length, and a fluid for whom we prescribe a uniform pressure p and a purely longitudinal, uniform velocity with speed |\vec{v}| at the inlet. The fluid itself has a constant density \rho and constant viscosity \mu, and is under the effects of gravity with an acceleration magnitude of |\vec{g}|. These quantities should fully define the flow throughout the pipe, whatever that flow may be.

A horizontal section of pipe with a radius R and an arbitrarily long length, containing a fluid with constant density and viscosity and under the effects of gravity. The pressure and velocity are prescribed uniformly at the inlet; as the flow passes through the pipe, it will either evolve into Hagen-Poiseuille flow or turbulent flow, something fluid-mechanical theory can’t tell us.

What was said above is worth repeating; every single characteristic of the flow downstream should be determined by a unique combination of these physical quantities. As a result, any physical quantity or property that we could possibly “read out” from this flow has to be a function of these 6 parameters and of these 6 parameters only.

One such readout could be some kind of turbulence index \zeta, a simple number that is 0 when the flow is Hagen-Poiseuille and 1 if the flow is turbulent. In accordance with the statement above, we can say this turbulence index \zeta is only a function of the above listed parameters:

\zeta = f(\rho, |\vec{v}|, |\vec{g}|, p, \mu, R)

Alright, this helps us out a little but not much. That being said, we’ve actually constrained the problem far more than it appears we have—to demonstrate, take a look at what happens if I claimed the following:

\zeta = \rho

From a mathematical perspective, there’s no problem with this, but physically, this doesn’t make a lick of sense; a number like the turbulence index can’t be equal to a density! This is like saying “A dog is three oranges” or “The sky is looking pretty book today”—it just doesn’t possess meaning as a physical description, because we aren’t relating objects that are the same type. In our case, whether or not objects are of the same type is determined by their physical dimensions. In fact, notice that it’s impossible to have any expression of the form \zeta = f(\rho) make physical sense, because there’s no mathematical function depending on the density that can take a density and return a simple number! And the same goes for every other physical variable that describes the system.

As a result of this, we find that we need to get a number on the right-hand side of our equality, and to get that to happen, we need to have our right hand side to be a function of just numbers as well. Now we face another question; how many independent number quantities can we construct out of those six physical parameters, and what are they? To do this, let’s list out all the quantities in the system and their units (in SI):

Physical quantityDensity \rhoViscosity \muRadius RPressure pSpeed |\vec{v}|Gravitational acceleration |\vec{g}|
Units (SI)\frac{kg}{m^3} \frac{kg}{m\, s} m \frac{kg}{m\, s^2} \frac{m}{s} \frac{m}{s^2}

Notice that the units of each quantity are built from multiplying and dividing three fundamental units: meters representing distance, seconds representing time, and kilograms representing mass. To get numbers with no units, it follows that we have to multiply and divide these quantities by each other to “cancel out” any and all fundamental units they may possess.

Let’s start out by trying to cancel out the gravitational acceleration term |\vec{g}|. As expected, it possess units of acceleration \frac{m}{s^2}, and so we it seems natural to find quantities with units of meters, seconds, or some combination of them to try and cancel its units out. If we decided to divide |\vec{g}| by the speed squared |\vec{v}|, that would get rid of all of our time units; multiplying that ratio by the radius R nets us something that doesn’t have units! Note that we didn’t need to cancel gravity out in any specific way; I just picked out a process to try and make the unit removal as simple as possible. In any case, let’s go ahead and call this number \Pi_1:

\Pi_1 = \frac{|\vec{g}|R }{ |\vec{v}|^2}

Note that we could do basically whatever we wanted to this quantity and still get a dimensionless number; we could square it, cube it, multiply it by a constant, exponentiate it, you name it; it will still be a totally valid number. In fact, most fluid mechanicists like to use the square root of the inverse of this number out of habit, which they call the Froude number Fr:

Fr = \sqrt{\frac{1}{\Pi_1}} = \frac{|\vec{v}|}{\sqrt{|\vec{g}|R}}

Alright, we’ve got one number down—let’s try to build another one. To guarantee independence between our numbers, let’s try to start the process of canceling units by picking some physical quantity that didn’t show up in our numbers before, like the pressure p. It has units of force per unit area \frac{kg}{m\, s^2}, so we can try canceling out the time units the same way we did for \Pi_1, by dividing by the speed squared. We can then get rid of the mass units by dividing by the density, leaving us with a quantity possessing only units of distance squared, which we can readily get rid of by dividing again by the radius squared. This second number we’ve just found, \Pi_2, is also commonly called the Euler number Eu:

Eu = \Pi_2 = \frac{p}{\rho|\vec{v}|^2 R^2}

Neat! With these two dimensionless numbers in hand, we’ve utilized every physical quantity in the system expect for the viscosity; it makes sense to start our search for a third number there. Viscosity has units of \frac{kg}{m\, s}, so we can try to divide by the speed and by the density to get rid of the time and kilogram units respectively. That leaves behind a single distance unit, which we can get rid of by dividing by the radius. This third number might look a little familiar if we explicitly show it:

\Pi_3 = \frac{\mu}{ \rho |\vec{v}| R} = \frac{1}{Re}

This is just a sneaky inverted form of our good friend the Reynolds number.

Now we can ask, are there any other independent dimensionless numbers we can construct out of these quantities? As it turns out, there isn’t—although there’s an infinity of different numbers we can make out of these quantities, they’d all wind up being functions of the three we saw above. This is because of a famous theorem called the Buckingham pi theorem, which states that the number of dimensionless numbers one can build out of a set of physical quantities is equal to the number of physical quantities minus the number of fundamental units they’re built from. (Proving this involves some relatively simple linear algebra, but it isn’t really germane to what we’re doing here.) In our case, we have 6 quantities and 3 units, so we get 3 dimensionless numbers as I stated before.

Given these numbers, we now have much tighter relationship between the turbulence index \zeta and the variables in the system:

\zeta = f(\Pi_1, \Pi_2, \Pi_3) = g(Fr, Eu, Re)

Remember that changing the specific dimensionless numbers we use as arguments doesn’t really matter, whether they be \Pi_1 or Fr—all the change does is alter the structure of the function, which we don’t know anyways.

Nearly all of the time, fluid mechanicians discard the dependence of the turbulence index on the Froude and Euler numbers by making some assumptions on the role of pressure in the flow. For example, in Hagen-Poiseuille flow, the absolute value of the pressure at either the inlet or the outlet didn’t play a physical role, just their difference; if we assume this to be true for turbulent flow as well, then the flow is independent of the pressure, which means it must be independent of the Euler number. Similarly, gravity didn’t play a role in the flow because all it did was generate a pressure gradient perpendicular to the flow, which doesn’t factor into the equation we used to solve for the Hagen-Poiseuille flow. All it does is make the absolute pressure at the inlet non-uniform, which we’ve already considered irrelevant, making the turbulence index be independent on the Froude number as well.

This gives us the desired relationship that experimentalists observe:

\boxed{\zeta = g(Re)}

As a result, the Reynolds number—which we saw before as an approximation to the ratio of convective to viscous effects—is now playing a dual role as a dimensionless parameter solely responsible for determining whether pipe flow is turbulent or not. And given our predetermined form of the turbulence index \zeta, which outputs one if the flow is turbulent and 0 if the flow is Hagen-Poiseuille flow, we know everything about the function g(...) except for the location where the “switch” from 0 to 1 happens. Experimentalists measure the switch to happen at about Re \approx 2300.

This process that we have just described required understanding absolutely nothing about the nature of the turbulent solutions to the Navier-Stokes equation, or the nature of the process by which flow transitions from Hagen-Poiseuille to turbulent, other than some assumptions on the role of pressure. In fact, it required understanding nothing about fluid mechanics at all except for understanding the physical quantities required to uniquely define a fluid flow.* This process is called dimensional analysis, and is far more powerful than it may appear at first glance.

*Granted, the entire reason we had to resort to this was because we weren’t finding a unique flow from this mathematical setup, but we are hypothetically observing only one type of flow for a given set of conditions, which is the flow we care about.

Consider what would happen if I as an experimentalist wanted to understand the critical average flow speed |\vec{v_{\text{crit}}}| of some fluid through a pipe of radius R_0 at which turbulence occurs, but only had access to a pipe of radius R_1. Well, because the turbulence index is only a function of Reynolds number, all I have to do is match them:

Re_{\text{crit}} =  \frac{ \rho |\vec{v_{\text{crit}}}|_1 R_1} {\mu} =  \frac{ \rho |\vec{v_{\text{crit}}}|_0 R_0} {\mu}

By doing some experiments to find the critical speed |\vec{v_{\text{crit}}}|_1 in the pipe I have access to, I can then solve for the critical speed |\vec{v_{\text{crit}}}|_0 I’m interested in:

|\vec{v_{\text{crit}}}|_0  =  \frac{|\vec{v_{\text{crit}}}|_1 R_1}  {R_0 }

I can now determine the critical speed at which turbulence occurs for a given fluid flowing through a pipe of arbitrary radius by doing experiments on just one pipe. Such a set of systems is said to possess similitude, and it is a key design tool for fluid dynamicists; an engineer designing a humongous oil pipeline can rest assured that an appropriately scaled & sped-up model of the pipeline in his lab will show the same non-turbulent/turbulent behavior as his final gargantuan product. One could also change the densities and viscosities too—as long as the Reynolds number is the same, the non-turbulent/turbulent behavior will be the same.

Two pipes with different radii and inlet velocities. Their inlet velocities are calibrated such that they both possess the same Reynolds number, which implies they will show the same non-turbulent/turbulent behavior. Engineers would refer to such systems as similar.

Most engineers would call it a day here, but we’re not done yet. Understanding how to trigger turbulence is one thing, but it doesn’t really give us any information about the direct physical consequence of it; how it affects the behavior of the flow in the pipe. More specifically, we don’t know anything about the pressure drop in a pipe flow when the flow is turbulent. For this, we can use a slightly different flavor of dimensional analysis.

To make straightforward comparisons with the results we found before from Hagen-Poiseuille flow, we can define some variable \left< \frac{\partial p}{\partial z} \right> representing the average pressure drop per length in a pipe flow, Hagen-Poiseuille or turbulent. This quantity has units of pressure per length, or \frac{kg}{m^2 s^2}. If we ignore the roles of absolute pressure and gravity by the arguments we justified above, we’ll inevitably find an expression for the average pressure drop per unit length of the form:

\left< \frac{\partial p}{\partial z} \right> = f\left(\rho, \mu,  |\vec{v}| , R\right)

We know that dimensional analysis restricts this expression more. In fact, we know that whatever is on the right-hand side has to possess units of pressure per length, and so we should set about constructing physical quantities with units of pressure per length out of the four quantities we listed above.

To save you the trouble, there’s only two such independent quantities we could construct: \frac{\rho|\vec{v}|^2}{R} and \frac{\mu|\vec{v}|}{R^2}. Consequently, one could expect that the expression representing \left< \frac{\partial p}{\partial z} \right> could be something of the following form:

\displaystyle \left< \frac{\partial p}{\partial z} \right> = \sum\limits_{n \in \mathbb{R}} \alpha_n \left( \frac{\rho|\vec{v}|^2}{R} \right)^n  \left(\frac{\mu|\vec{v}|}{R^2}\right)^{1-n}

Those \alpha_n are dimensionless, and so must be a function only of the single dimensionless quantity we can construct out of this system—the Reynolds number.

Now, recall what we said about turbulent flows before; their Reynolds tensors are decidedly non-zero. In fact, it’s reasonable to expect that as we crank up the Reynolds number in a turbulent pipe flow, the Reynolds tensor is going to increase as well. Additionally, we could crank it up enough such that the convective effects completely dominate the viscous effects! In that scenario, we’d find that all of the flow characteristics are virtually independent of the viscosity.

If that’s the case, then the expression for the pressure drop per length we found above is tightly constrained by the fact that every term with the viscosity in it has to vanish, since the pressure drop per length can’t depend on it. This leaves us with just the n = 1 term:

\displaystyle \left< \frac{\partial p}{\partial z} \right> =  \alpha_1 \frac{\rho|\vec{v}|^2}{R}

We had originally stated that \alpha_1, like the other coefficients, had to be a function of just the Reynolds number. But alas—the Reynolds number is a function of the viscosity! The only way to reconcile the fact that the viscosity can’t influence the physics with this Reynolds-viscosity dependence is if \alpha_1 wasn’t a function of anything; if it was just some arbitrary constant number, totally independent of the Reynolds number. I like to call this phenomenon, where a reduction in the number of dependent physical quantities makes it impossible to construct a dimensionless number, a dimensional crisis. Sounds cool, doesn’t it?

Anyways, we’ve found the following result at high-Reynolds numbers for pipe flow:

\displaystyle \left< \frac{\partial p}{\partial z} \right>_{\text{high Re}} = \alpha_1 \frac{\rho|\vec{v}|^2}{R}

where \alpha_1 is just some constant number.

As a sanity check, we could try the same procedure when the Reynolds number is low and see if what we get matches the Hagen-Poiseuille result. In that scenario, the Reynolds tensor should be small (it is in fact identically zero), and so the viscosity-proportional term should dominate the behavior of the system. To that end, the average pressure drop per length should be independent of the only physical variable not present in the \frac{ \mu|\vec{v}| }{R^2} term, the density.

That means we kill off every term in the above sum except for the n=0 term, and again trigger a dimensional crisis, leading to the following expression at the low-Reynolds number limit:

\displaystyle \left< \frac{\partial p}{\partial z} \right>_{\text{low Re}} = \alpha_0 \frac{ \mu|\vec{v}| }{R^2}

This precisely matches the expression we got before for Hagen-Poiseuille flow; that expression lets us determine that \alpha_0 = 8. Up to that factor of 8, we could have avoided all the differential equation-solving we did in the previous Part by making the assumptions we did above and performing dimensional analysis; we would have gotten the same result knowing absolutely nothing about differential equations. That’s pretty awesome!

I should mention that fluid mechanicists do a very cute trick when studying this average pressure drop per length; they divide the pressure drop per unit length by one of those constructed quantities with units of pressure per length, almost always the \frac{\rho|\vec{v}|^2}{R} expression—this expression is called the Fanning friction factor f:

\displaystyle f = \frac{ \left< \frac{\partial p}{\partial z} \right>}{ \frac{\rho|\vec{v}|^2}{R}}

The reason they do this is because now you wind up getting a dimensionless number, and we already showed that this number can only be a function of the only other dimensionless number we can construct in this system, if it’s a function of anything at all—the Reynolds number. Hence, plotting the Fanning friction factor against the Reynolds number tells me the pressure drop behavior of every single possible pipe flow system of the type we describe here by virtue of similarity.

Dividing the previous results we got by this \frac{\rho|\vec{v}|^2}{R} expression, we find the following expressions for the friction factor at low and high Reynolds numbers:

\displaystyle f_{\text{low Re}} =  \frac{\alpha_0 \frac{ \mu|\vec{v}| }{R^2}}{ \frac{\rho|\vec{v}|^2}{R} } = \alpha_0 \frac{\mu}{\rho|\vec{v}|R} = \frac{\alpha_0}{Re}

\displaystyle f_{\text{high Re}} =  \frac{ \alpha_1 \frac{\rho|\vec{v}|^2}{R} }{ \frac{\rho|\vec{v}|^2}{R} } = \alpha_1

Having done no calculus, and relying almost exclusively on dimensional analysis, we have established reasonable hypotheses for the behavior of the Fanning friction factor as a function of the Reynolds number; in a logarithmic plot, it will show up as a straight line with slope of -1 at the start, followed by a vertical kink somewhere due to the discontinuous transition from Hagen-Poiseuille to turbulent flow, while eventually flattening out into a straight horizontal line as the Reynolds number gets higher.

Lo and behold, look what the experimentalists got:

A logarithmic plot of the Fanning friction factor versus the Reynolds number in a pipe as determined by experimental data. The friction factor is proportional to the inverse of the Reynolds number for low-Re, and independent of the Reynolds number at high-Re. This diagram is called a Moody diagram, and fluid mechanicists often also include an empirical value called surface roughness in it to justify pipe irregularity effects. All surface roughness values lead to the same qualitative behavior as shown above.

That’s one hell of a magic trick! Using the behavior at the “ends” to come close to describing the behavior in the “middle” is squarely in the ballpark of a field of math called asymptotics, and for my next trick, we’re going to use it to answer the last remaining question of introductory fluid mechanics—how to characterize flow over an immersed object.

Things to Think About

1. Before we described the Reynolds number as being an approximation of the ratio of convective effects to viscous effects. How would you describe the Froude and Euler numbers in a similar way? Is there a single way to do so?

2. Try to find a physical system you don’t understand and use the dimensional analysis techniques I demonstrated above to get some practical answers from them. When is it helpful? When isn’t it?

3. Sometimes dimensional analysis lets us construct differential equations we can then solve. Consider dipping a hot sphere of radius R\  (m) with a specific heat C\ (\frac{kg\,m^2}{s^2 K}) into a vat of cold liquid with a fixed ambient temperature far from the sphere. If you assume the rate of change of the temperature with time \frac{dT}{dt} is only a function of these two parameters, of the current temperature of the sphere T (K) , and of the heat transfer coefficient h\  (\frac{kg}{s^3 K}) , find an expression for \frac{dT}{dt} using dimensional analysis. Compare what you get to what engineers use. Try solving it, if you can! No knowledge of thermal physics is needed.

4. Physicists in the early 20th century noticed that particles at the quantum scale behave like waves. Knowing nothing about quantum mechanics other than that this wave is a function of the particle’s mass m (kg), the particle’s speed u (\frac{m}{s}), and the Planck constant h (\frac{kg\, m^2}{s}), determine the wavelength (m) and frequency (\frac{1}{s}) of this matter wave up to a proportionality factor. Compare it to what quantum physicists theorized.

5. Do you agree with my arguments for why the Euler and Froude numbers don’t matter for turbulent pipe flow? Can you think of other fluid-mechanical systems where you should have to consider them?

6. Think about how you would exploit similarity to analyze the flow of a substance like oil through a pipe versus a substance like water. What would you do?

Part IX: Pipe Flow & Turbulence

With the analytical tool of the Navier-Stokes equation in hand, we can finally begin to make statements about the precise shapes and speeds of flows. Perhaps the best place where we can make those statements without running afoul of the pathology of the Navier-Stokes equation—and trust me, we will run afoul of it—is in studying flow in pipes or channels, which we described using hydraulic equations in Part V.

Let’s again consider steady incompressible flow through a horizontal section of pipe with length L and with a uniform cross-section. This time, we’re going to try and explicitly define what the momentum loss \vec{H} is using the Navier-Stokes equation, and get a formula for the shape and speed of the flow inside of the pipe as a function of the other free parameters in the system. To do this, we’ll need to do a couple of tricks to bring everything down from the integral formulation we described in Part V—the fundamental equations of hydraulics—to a differential one.

Now, the first fundamental equation of hydraulics tells us that, because the fluid is incompressible, conservation of mass requires that the speed at the inlet of the pipe be the same as the speed at the outlet, meaning that all the momentum coming into the pipe from flux at the inlet must leave at the outlet. In addition, all the momentum being added into the pipe due to gravity is acting perpendicular to the momentum loss; this immediately indicates that the momentum loss in steady flow within such a pipe must necessarily be compensated by a drop in the pressure (and only a pressure drop) along the flow direction.

A horizontal pipe of length L with constant cross-sectional area and the associated tractions/force densities acting on it. Conservation of mass indicates that the speeds on either end are equal; this means the loss coming from \vec{H} can only be balanced by a drop in the pressure between the ends due to conservation of momentum.

A nice thing about this result is that it is totally independent of the length of the pipe section! We can make the pipe as long or as short as we want and the above statement still holds. Now I’m going to do a classic math trick; I’m going to pick some horizontal line across the pipe, and represent the pressure through that arbitrary line as a function of the distance from the inlet z with a Taylor series off of the value at the inlet portion of the line:

p = p_1 + G_1 z + G_2 z^2 ...

If I wanted to represent the pressure drop, I could just subtract by the inlet value:

\Delta p = G_1 z + G_2 z^2 ...

If I take the length of the pipe and shrink it down more and more, the squared term is going to get smaller than the linear term as long as both terms are nonzero. This means that, for an infinitesimally short length of pipe, the pressure drop is always strictly linear in the length, and so the pressure drop per length is a constant, G_1.

Now think about this; if I assume the momentum loss is independent of the pressure, I can take the section of pipe and split it up into a bunch of infinitesimally tiny segments one behind the other, which should all be essentially identical to each other; sure, the pressure at each of the inlets might be different, but the pressure drop should be the same, as well as everything else! This implies the pressure drop is linear in all of the tiny segments, each with the same pressure drop per length G_1, so the pressure drop per length in a finite section of such a pipe must be constant (under the given assumptions, of course).

This gives us what we need to analyze pipe flow through a horizontal section of pipe with constant cross-sectional area using Navier-Stokes. In this scenario, based on the argument above, we are going to impose some arbitrary pressure drop per length through the pipe G_1, and see what the resulting flow field looks like. The argument above guarantees (or at least, suggests) that imposing such a pressure drop is physically sensible, and will lead to a sensible solution of the Navier-Stokes equation. Before we do, however, notice a very curious quirk in our logic—if we assume the pressure drop per length is constant, that means that the pressure in the pipe through some arbitrary horizontal line is going to be of the form:

p = p_1 + G_1 z

where p_1 can be the pressure in an arbitrary location of the pipe. No matter how big that arbitrary reference pressure is, there is always some length of pipe in one direction that causes the pressure to go negative, which makes no physical sense! In that sense, we find that any solution we get for such a problem can’t apply to arbitrarily long pipes, as we will invariably run into a region of pipe where the solution we find simply cannot exist—an existence failure, as we saw in Part VII, and a harbinger of things to come.

Alright, let’s get on with the solution already. We want to solve the incompressible Navier-Stokes equation for steady flow inside of a cylindrical pipe with a uniform cross-sectional pressure and a constant pressure drop per length in the “pipe” direction:

\displaystyle \rho\left(\frac{\partial \vec{v}}{\partial   t}\ + \nabla \vec{v} \cdot \vec{v} \right) = \vec{f} - \nabla p\ +\    \mu \nabla^2 \vec{v}

Let’s list out a couple of observations and “common-sense” assumptions that get us to a point where we can handle this equation:

  1. We assumed the flow is steady already, so the time rate of change term is zero: \displaystyle \rho\frac{\partial \vec{v}}{\partial   t} = 0.
  2. There’s no flow in any direction other than in the “pipe direction”/longitudinal direction; this means there’s no flow in the radial direction or the angular direction (swirling flow): \vec{v} = [v_r\ v_\theta\ v_z] = [0\ 0\ v_z]
  3. Everything about the flow is totally independent of the pipe angle, so we can kill off any terms involving changes in the angular direction: \frac{\partial ...}{\partial \theta} = 0
  4. Gravity is the only external force and acts downwards, so it only causes pressure gradients perpendicular to the flow direction.

If we take a look at the incompressible conservation of mass equation and take into consideration these assumptions, we see something that is fairly obvious:

\nabla \cdot \vec{v} = \frac{1}{r}\frac{\partial(r v_r)}{\partial r} +  \frac{1}{r}\frac{\partial(v_\theta)}{\partial \theta} + \frac{\partial v_z}{\partial z} = 0

\frac{\partial  v_z}{\partial z} = 0

The flow doesn’t change as we move along the pipe/longitudinal direction. This is going to do something drastically nice for our problem; it’s going to kill off the convective flux term we spent so much time whining about!

\nabla \vec{v} \cdot \vec{v} = v_z \frac{\partial v_z}{\partial z} = 0

That means that the Reynolds tensor for this problem is exactly 0, and we shouldn’t have to worry about shocks and all of those other horrible things we saw in the previous Part because the Navier-Stokes equation becomes linear. In fact, it becomes something almost too simple to recognize:

\displaystyle 0 = - \nabla p\ +\    \mu  \nabla^2 \vec{v} +\rho \vec{g}

The Navier-Stokes equation is telling us that, as we saw before, the pressure drop down the pipe length needs to be balanced out by a momentum loss or viceversa; but our knowledge of rheology in a fluid described by Navier-Stokes is telling us specifically how that momentum loss is related to the velocity of fluid in the pipe, giving us the connection to the velocity we needed!

Recalling that the equation we’re dealing with relates vectors, we can now begin to look at what this equation tells us for the individual components of these vectors. If I looked at the component of the vectors pointing in the direction of gravity in this equation, which I’ll label as the \hat{y} direction, I’d get a simple relationship that looks very familiar:

0 = -\frac{\partial p}{\partial y} + \rho g

This is just the hydrostatic equation we saw before! That means that the pressure of a fluid in a horizontal pipe increases in the direction of gravity consistent with what you’d observe in a horizontal pipe full of standing fluid. This also has the convenient effect of ensuring that gravity doesn’t affect the flow; it only affects the absolute value of the pressure, which we assumed previously as not affecting the momentum loss.

Now we can take a look at the components of the vector that are perpendicular to gravity, pointing down the length of the pipe. Applying the assumptions we listed above and doing some simplifications, we find that:

\displaystyle 0 =  -\frac{G_1}{\mu} + \frac{1}{r}\frac{\partial(r \frac{\partial v_z}{\partial r})}{\partial r} = \frac{1}{r}\frac{\partial v_z}{\partial r} + \frac{\partial^2 v_z}{\partial r^2}

Not only is this equation an ordinary differential equation—we’re only taking derivatives in the radial direction r—but it’s also linear in the velocity! Doing some calculus tricks to solve it, we find the general solution:

v_z(r)=c_1 \log (r)+c_2+\frac{G_1 r^2}{4 \mu}

We need c_1 = 0 or the velocity will go to negative infinity at the center of the pipe, and we need the velocity at the edge of the pipe to be zero to ensure that the velocity doesn’t discontinuously change at the pipe/fluid interface—giving us a value of c_2 such that we finally have our velocity for pipe flow:

v_z(r)=\frac{G_1 (r^2 - R^2)}{4 \mu }

\boxed{\vec{v} = \left[0 \enspace 0\enspace \frac{G_1 (r^2 - R^2)}{4 \mu }\right]}

This type of flow is called Hagen-Poiseuille flow, and is perhaps the most famous and well-known flow in fluid mechanics. It’s parabolic in shape, is 0 at the pipe edges by definition, peaks at a value of v_{max}=\frac{G_1 R^2}{4 \mu } at the center of the pipe, and has a cross-sectional average of \langle v \rangle = \frac{G_1 R^2}{8 \mu }.

A view of Hagen-Poiseuille flow within a pipe for specific parameter values. The velocity drops to 0 at the pipe edges, denoted in black, while the velocity peaks in the center. We have assumed there is no flow in any direction other than the longitudinal one.

This finally gives us an expression for the momentum loss in a horizontal, constant cross-section pipe as a function of the velocity, by solving for the pressure drop and noting the drop in pressure is equivalent to the momentum loss from the arguments we made at the start:

|\vec{H}| = \frac{8 \mu \langle v \rangle L}{R^2}

In a kinder, more just world, this would be all there is to it. Unfortunately, when researchers started directly looking at pipe flows, they noticed that everything went according to theory—until they cranked the speed up high enough:

A simulation of pipe flow. Notice how the observed flow profile is totally inconsistent with the Hagen-Poiseuille flow solution we obtained above.

What the hell is going on? Well, as the existence failure for the pressure had been warning us, it appears as if there’s another fundamental problem we didn’t spot when we began to set up pipe flow; non-uniqueness! As it turns out, there might not be just one possible flow for a given condition on a pipe inlet, but a potentially infinite number of them—and from the experiments, it seems like those other flows possess all the things we had assumed away when we derived Hagen-Poiseuille flow; they’re unsteady, they have nonzero convective momentum flux, they have nonzero flow in all three pipe directions, etcetra! These flows, which undoubtedly have nonzero Reynolds tensors, are called turbulent flows, and are omnipresent in all of fluid mechanics.

Comparison sketches between Hagen-Poiseuille and turbulent pipe flows. Hagen-Poiseuille flow is orderly and symmetric; turbulent pipe flow possesses eddies, whorls, and few if any symmetries. When the parameters of the fluid/flow are changed in specific ways, the flow transitions from one to the other.

One way fluid mechanicists explain what is happening is that when the average speed is low through the pipe, nature “picks” out the Hagen-Poiseuille flow from all the other possible valid pipe flows because the litany of assumptions we made to get to Hagen-Poiseuille flow made sense; but when the average speed is sufficiently large enough, nature changes its mind and picks out one of the extremely complicated chaotic flows instead. In mathematics, this is commonly referred to as a bifurcation. Nobody has ever been able to come up with an analytical expression for any of these alternate chaotic flows, precisely determine a theoretical criteria for this transition, or even determine if this non-uniqueness hypothesis is true for certain—and so we are left with experiments and numerical simulations to fill this knowledge gap for now.

But what experimentalists could conclude was curious; they observed that once a specific dimensionless number crossed a threshold, the Hagen-Poiseuille flow would switch into turbulent flow and back. That dimensionless number is the Reynolds number, which we saw in the previous Part as an approximation to the ratio of convective to viscous effects. But that interpretation doesn’t make any sense here; the convective effects in Hagen-Poiseuille flow are exactly zero, so saying that the flow transitions from laminar to turbulent because the ratio of convective to viscous effects in pipe flows becomes too large is gobbledygook. Clearly, the Reynolds number here represents something else, and to discover what it represents, we’ll need to take stroll down a very curious field of mathematical physics called dimensional analysis.

Things To Think About

1. List out all the assumptions I made to justify the Hagen-Poiseuille flow. When do these assumptions make sense? When don’t they?

2. In a real pipe system, what happens when the pressure drops “too much”? Do the Navier-Stokes equations apply then, and if not, what would you need to change?

3. How would you try to understand turbulence? Do you agree or disagree with the explanation above?

4. What happens if you try to include gravity? Are you able to get a solution to the Navier-Stokes equations then?

5. What happens if the pipe has a constant cross-section that isn’t circular? Could you get a solution similar to Hagen-Poiseuille, and if you do, what would drive the differences?

6. Based on the results we got for pressure in the downwards and longitudinal direction, can you write a complete expression for the pressure in a pipe section? If you want to practice changing coordinate frames, use cylindrical coordinates.

Part VIII: Navier-Stokes, Existence & Uniqueness

At long last, we’ve finally managed to derive the famous Navier-Stokes equation; the cornerstone of most of fluid mechanics. Notice all the assumptions we had to make to get to it! To refresh your memory, this is what it looks like, along with descriptions of each term in it:

The Navier-Stokes equation, which describes momentum density transport in a simple, Newtonian, isotropic, incompressible fluid. The terms on the left-hand side are mathematical consequences of the concept of transport, while the terms on the right-hand side (save for the external force density) are defined by intermolecular fluid interactions.

There are many different, equivalent ways to write the Navier-Stokes equation, but the most common is this one, which incorporates a couple of simplifications thanks to conservation of mass:

\boxed{\displaystyle \rho\left(\frac{\partial \vec{v}}{\partial  t}\ + \nabla \vec{v} \cdot \vec{v} \right) = \vec{f} - \nabla p\ +\   \mu \nabla^2 \vec{v}}

(Note that some people write the dot product in the flux term backwards based on their notation for the gradient, but it’s all just preference.) Another equivalent form I briefly mentioned when discussing transport theory, which is particularly useful when discussing the curvature of flows, is commonly known as the Lamb form:

\displaystyle \rho\left(\frac{\partial \vec{v}}{\partial  t}\ + \frac{\nabla \left( \vec{v} \cdot \vec{v}\right)}{2} - \vec{v}\times\left(\nabla\times\vec{v}\right) \right) = \vec{f} - \nabla p\ -\   \mu \nabla \times \left(\nabla \times \vec{v}\right)

In whatever form it takes, this equation (along conservation of mass and varying simplifying assumptions) will be the chief mathematical tool we’ll use to determine the shape and speeds of flows. That being said, it turns out that the Navier-Stokes equation is notoriously difficult to deal with except in the simplest of situations, and it’s all thanks to the flux term:

\rho \nabla \vec{v} \cdot \vec{v}

This little mathematical pest has been the bane of fluid mechanicists for a solid hundred years, and it induces such a monumental headache on any mathematician attempting to solve the Navier-Stokes equation that there is actually a million-dollar bounty out for anyone who is even able to prove that the Navier-Stokes equation always has “sensible” solutions for a physical input—this is one of the most important mathematical problems ever conceived!

The reason this flux term is such a problem is because of a property called nonlinearity, which is best shown rather than told. To demonstrate, let’s consider a silly non-physical problem in 1-D, where we consider two imaginary fluids whose velocity is solely determined by the following equations:

Fluid 1Fluid 2
1-D Equation\rho v\frac{dv}{dx} = f \mu \frac{d^2v}{dx^2} = f
3-D Equivalent \rho \nabla \vec{v} \cdot \vec{v} = \vec{f} \mu \nabla^2 \vec{v} = \vec{f}

Applying some calculus tricks, what you’ll wind up finding is that the expressions for the velocities in each of the fluids is markedly different:

Fluid 1Fluid 2
1-D Solutionv(x) =\pm \frac{\sqrt{2} \sqrt{c_1 \rho +f x}}{\sqrt{\rho }} v(x) = c_2 x+c_1+\frac{f x^2}{2 \mu }

Note that for Fluid 2, if I set the constants to 0, the velocity is linearly proportional to the force. This means that if I double the force, the velocity everywhere is going to be doubled, and so on—a property unsurprisingly called linearity. For Fluid 1, however, this isn’t true; setting the constant to 0 and doubling the force instead leads to an increase in the velocity by a factor of \sqrt{2}! As innocuous as it seems, nearly every single tool ever made by mathematicians and scientists to solve partial differential equations like the Navier-Stokes equation relies on linearity, meaning any attempts to solve the Navier-Stokes equation relies on conceptual machinery on the outside of most of mathematical history. But that’s not even the worst part; as it turns out, nonlinearity is often comorbid with a host of horrendous properties, and the velocity of Fluid 1 is no exception.

For example, consider what would happen if I attempt to find the flow velocity for Fluid 1 if I take the velocity at the origin to be zero, v(0) = 0; because of that \pm sign in front of the velocity expression for Fluid 1, I actually have two possible fluid velocities! This is incredibly strange from a physical perspective; one could conjecture a situation where I have two exact copies of this system, where I’m applying the same force under the same conditions to each of them, and somehow get different flows. This goes straight against any rational understanding of physics, and is called non-uniqueness. Even worse, no matter what value of the constant c_1 I choose based on some physical information, or which expression for the velocity of Fluid 1 I pick, there is always some value of x past which my velocity suddenly turns into an imaginary number! This means that in some regions, we don’t even have an expression for the flow velocity, which is even more mind-boggling from a physical perspective and is referred to by mathematicians as existence failure (this is a really cool name). And this is only in one dimensions, folks—just imagine how much worse it gets in 3.

Two solutions for the velocity of Fluid 1, each solving v(x)\frac{dv(x)}{dx} = 1 with zero velocity at the origin, corresponding to a fluid with f, \rho = 1. Neither solution is defined in the region x < 0.

When you include time into the equations, things get even worse; you can start off with a flow that makes sense, and find it evolving into something seemingly nonsensical! Take, for example, what happens if you now take the ODE describing Fluid 1 and add a time rate-of-change term, turning it into a PDE:

\frac{\partial v}{\partial t} + \rho v\frac{\partial v}{\partial x} = f

This equation represents momentum transport in a one-dimensional fluid with no intermolecular interactions, and is commonly called Burgers’ equation. Take a look at what happens if you take a perfectly nice initial flow field for such a fluid and let it evolve in time for a constant unit f and \rho:

The solution to the PDE \frac{\partial v}{\partial t} +  v\frac{\partial v}{\partial x} = 1 for a smooth initial condition. The function describing the velocity eventually develops a discontinuity in space, causing the solution to become ill-defined. The numerical solver fails to handle this phenomenon, and glitches out once it develops.

After time passes, the velocity becomes discontinuous, and jumps immediately from one value to another; a phenomenon commonly referred to as a shock or shock wave. Mathematically, this is a problem not just because \frac{dv}{dx} = \infty at the discontinuity, but because the velocity has effectively two values there, causing it to be ill-defined! What happens here, from a philosophical perspective, is astounding; we’ve gone and defined an object mathematically using the rules of calculus, and the object has then evolved over time into something that violates the very mathematical principles we used to define it. This is exactly the problem (or at least, one of them) people run into when trying to solve the Navier-Stokes equation; the objects one uses to define flows for fluids can very rapidly outgrow the conceptual framework we used to define them.

Consequently, nearly every analytical attempt at using the Navier-Stokes equations to describe a flow necessitates ensuring that the momentum flux term in the equations is either negligible or handled indirectly. One way to verify this negligiblity is by taking the ratio of the magnitudes of each component of the flux term with the magnitudes of the components of the only other term that involves spatial gradients of the velocity—the viscosity term:

\mathbf{Re} = \rho \nabla \vec{v} \cdot \vec{v} \oslash \mu \nabla^2 \vec{v}

The \oslash symbol represents that we are dividing every term of a tensor by every other term of another tensor, effectively acting like the division counterpart of the outer product \otimes. Since each of these quantities is a 3-D vector, this ratio actually contains 9 distinct ratios, and I refer to it as the Reynolds tensor. The Reynolds tensor is useful in the sense that it is telling us whether or not, for a given flow, we should expect to observe the same kind of mind-boggling behavior we saw above in the Burgers’ equation, where it would be located, and it what directions we expect it to manifest. Because calculations involving the Reynolds tensor are pretty unwieldy, fluid mechanicists almost always single out a specific component of the Reynolds tensor, guess “characteristic” values for the density/viscosity/velocity/velocity gradients/lengths that they assume are valid everywhere in the flow, and come up with a single scalar ratio called the Reynolds number which they use instead. Sometimes this leads to spurious assumptions and misdrawn conclusions, but alas.

The components of the Reynolds tensor in Cartesian coordinates, along with a typical Reynolds number-type approximation. Calculating the Reynolds tensor is often unwieldy, so most scientists simply estimate a single number based on quantities in the problem they expect to see and call it the Reynolds number, shown below.

Nearly all of the time, we use the Reynolds tensor/number to do sanity checks on any simplifying assumptions we make to the Navier-Stokes equation—either by using a conjectured Reynolds tensor/number to justify removing some term from the Navier-Stokes equations, or by calculating the Reynolds tensor/number for a solution to a simplified Navier-Stokes equation to validate whether or not the simplification was sensible in the first place. And one setting where we’ll be able to get rid of that flux term and, as a result, easily construct specific conclusions about flow shapes and speeds is a setting we’ve looked at before—pipe flow.

Things to Think About

1. Some people think adding a term representing viscosity to equations like the ones described above help get rid of the shocks & other strange behaviors we saw just now. Try solving the ODE \rho v(x)\frac{dv}{dx} - \frac{d^2 v}{dx^2} = f for a bunch of different \rho/f/\mu using some math computer programs to see whether or not that line of thinking is correct.

2. How would neglecting the flux term in a problem change the physics of the fluid you would observe, and would it make sense? Is there a fluid that can possess a zero flux term?

2. In what situations do you think simply using the Reynolds number instead of the Reynolds tensor would be sufficient? In what situations would it fail?

3. Try calculating the Reynolds tensor for simple flows; maybe a constant flow, or a flow that changes linearly in one direction in space. What components are zero, and which aren’t? What do they correspond to?

4. Shock waves appear to make no sense in a differential equations sense, but seem to make perfect sense physically; is there some mathematical tool we were using before that would be able to make mathematical sense of shock waves?

5. In Part II, we talked about how pressure ensures liquid doesn’t just accumulate into a hyperdense thin film at the bottom of a cup. How is that connected to the shock wave phenomenon we saw above?

6. If you want to get a glimpse of how linearity is useful for solving differential equations, solve \frac{dv}{dt} = f for f=1 and 2 for arbitrary constants. Can the sum of both of these solutions represent a solution for \frac{dv}{dt} = 3?

Part VII: Rheology

As we saw during our discussion on transport theory, the principal equation of fluid dynamics really boils down to a transport equation for the momentum in a fluid, or for the transport of momentum density through fluid points:

\displaystyle \frac{\partial\left(\rho\vec{v}\right)}{\partial t}\ + \nabla \cdot \left(\rho \vec{v} \otimes \vec{v}\right) = \vec{f} + \vec{r}_{mol}

where the \vec{r}_{mol} term represents forces stemming from molecular effects. It is this molecular term which exclusively “identifies” a fluid, as everything else in the governing equation above was generated through a very general mathematical framework describing the movement of “things” in space rather than through specific physical insight associated with fluids. Consequently, understanding this connection between continuum-scale momentum and molecular effects is a humongously important task for physicists and engineers alike, and is the principal endeavor of rheology.

It’s often extremely difficult to generate refined models for \vec{r}_{mol} from first principles, so often the process of determining an \vec{r}_{mol} for a given fluid consists of making some basic assumptions about the fluid and then testing the fluid experimentally to see if the behavior of the fluid is consistent with those assumptions. Rheology, as a result, is largely an experimental science.

This doesn’t mean we’re simply left to fiddle around with rheometers from now on, however—in fact, we’ll rapidly find that a couple of simple assumptions and some not-so-simple tensor calculus will lead us to an equation that applies to nearly all fluids.

Firstly—and as we saw in Part 1—when we talk about the density or velocity of a fluid at a point, we’re really talking about the average velocity or mass/volume ratio inside a small “molecular net” around a point. One thing we haven’t talked about yet is on the consequences of formulating a theory of fluid mechanics strictly on these average quantities, and how fluctuations around those averages affects the way fluid molecules exert forces on each other at the macroscale. In order to do so, we first need to make a relatively sensible hypothesis about the nature of these molecular forces, a hypothesis I call the “fluctuation hypothesis”; that all random fluctuations of fluid quantities occur as a result of intermolecular forces in the fluid. Based on this, we can try to formulate the Cauchy momentum equation we derived before that excludes an implicit intermolecular force term—but for a randomly fluctuating density and velocity—and equate any differences we see between either momentum equation to intermolecular forces.

Let’s start by defining a “true” fluctuating density \rho^* and velocity \vec{v}^* at a point, both of which fluctuate independently around an average \rho and \vec{v}. This lets us split up the true quantities into an average part, which we’ve been dealing with from the beginning, and a fluctuating part (indicated with ') that is on average zero:

\rho^* = \rho + \rho'
\vec{v}^* = \vec{v} + \vec{v}'

Since these quantities are transported in a fluid just like their averaged counterparts, we can construct a momentum equation based on them that is just as valid as the one we saw before:

\displaystyle \frac{\partial\left(\rho^*\vec{v}^* \right)}{\partial t}\ + \nabla \cdot \left(\rho^* \vec{v}^* \otimes \vec{v}^* \right) = \vec{f}

Expanding out in terms of average and fluctuating components leads to a lengthy, messy expression:

\displaystyle \frac{\partial\left(\rho\vec{v}\right)}{\partial t}\ + \frac{\partial\left(\rho'\vec{v}'\right)}{\partial t}\ + \nabla \cdot \left(\rho \vec{v} \otimes \vec{v}\right) + \nabla \cdot \left(\rho' \vec{v} \otimes \vec{v}\right) + \nabla \cdot \left(\rho \vec{v}' \otimes \vec{v}\right)\ + \nabla \cdot \left(\rho \vec{v} \otimes \vec{v}'\right) + \nabla \cdot \left(\rho' \vec{v}' \otimes \vec{v}\right) + \nabla \cdot \left(\rho' \vec{v} \otimes \vec{v}'\right) + \nabla \cdot \left(\rho \vec{v}' \otimes \vec{v}'\right)\ + \nabla \cdot \left(\rho' \vec{v}' \otimes \vec{v}'\right) = \vec{f}

This just looks like a big jumble of symbols, so we have to get to work on simplifying this a little bit. Luckily, we can get rid of nearly everything in one fell swoop by averaging everything out! This will get rid of every term that contains a single fluctuating term (since those terms are normal numbers multiplied by something that’s on average zero) and everything that has a \rho'\vec{v}' in it (since both are independent and average to zero, they multiply to something that’s on average zero). This results in a far simpler expression:

\displaystyle \frac{\partial\left(\rho\vec{v}\right)}{\partial t}\ + \nabla \cdot \left(\rho \vec{v} \otimes \vec{v}\right) + \nabla \cdot \left(\rho \overline{\vec{v}' \otimes \vec{v}'}\right) = \vec{f}

where the \overline{\phantom{x}} indicates an average of something that wasn’t already defined as an average. Curiously, the \rho \overline{\vec{v}' \otimes \vec{v}'} doesn’t zero out! This is because even though each individual \vec{v}' is on average zero, the combination \vec{v}'\otimes\vec{v}' isn’t, since it’s “kind of” the square of an individual \vec{v}'. The presence of that extra term shows that considering fluctuations in fluid mechanical quantities actually does lead to a verifiable change in the momentum transport equations we’ve been setting up! And since we had attributed every possible extra term in the fluctuating form of the Cauchy momentum equation to intermolecular forces—that was our fluctuation hypothesis—then for both momentum equations to match, it has to be true that:

\vec{r}_{mol} = -\nabla \cdot \left(\rho \overline{\vec{v}' \otimes \vec{v}'}\right)

This gives us a look at how molecular-scale phenomena in a fluid leads to forces at the macroscale, but doesn’t go very far in actually giving us a relationship between those molecular quantities and the averaged-out properties of a fluid we use to formulate fluid mechanics. Remember, \vec{v}' could be basically anything! So what fluid mechanicists do is collect everything inside of the parenthesis, flip the sign, and label the grouped object as a new, unknown second-order tensor \boldsymbol\sigma commonly referred to as the stress. This leads to the form of the Cauchy momentum equation that most engineers \& scientists use:

\displaystyle \frac{\partial\left(\rho\vec{v}\right)}{\partial t}\ + \nabla \cdot \left(\rho \vec{v} \otimes \vec{v}\right) = \vec{f} + \nabla \cdot \boldsymbol\sigma

Equivalently, the stress tensor \boldsymbol\sigma is just a second-order tensor that represents, for a given fluid point, the traction/force per unit area \vec{t} acting on that element thanks to its neighboring fluid point in the \hat{n} direction. If you’ve studied solid mechanics, you’ve likely heard the term stress thrown around before; it means the same thing there as it does here. In fact, the only difference between a solid and a fluid is that the stress is a function of different continuum parameters in each type of substance. Notice we never did anything fluid-specific when we derived the momentum transport equation above!

Luckily, we already know one piece of the puzzle of \boldsymbol\sigma; it has to contain the forces coming from pressure gradients, since we spent two Parts discussing how pressure (and the forces it induces) is an entirely molecular effect. We can make our lives easier and pull that pressure gradient term out, making sure we make it negative so that everything is consistent with the hydrostatic equation we derived in Part II:

\displaystyle \frac{\partial\left(\rho\vec{v}\right)}{\partial   t}\ + \nabla \cdot \left(\rho \vec{v} \otimes \vec{v}\right) = \vec{f} - \nabla p +  \nabla \cdot \boldsymbol\tau

That new tensor \boldsymbol\tau is usually called the deviatoric stress tensor, which I think is a silly name.

The contents of the deviatoric stress tensor are in some sense what we’re really looking for, since the notion of pressure and pressure gradients is universal in solids and liquids. And here, we can classify fluids rheologically into two important & distinct branches based on a key characteristic of their behavior when standing still; fluids whose static behavior is fully characterized by pressure gradients (commonly called simple fluids), and fluids for which it is not (commonly called non-Newtonian fluids). For the latter, the molecules in the fluid will usually have some kind of extra molecular force that “causes” them to try to remain in some specific macroscopic shape or configuration—just like a solid. Some examples of this are ketchup (it doesn’t like to flow out unless you smack the bottle or pump it out), mucus, pancake mix, melted rubber, etcetra; there’s too many to count. A simple experimental way to figure out whether or not a fluid is simple or non-Newtonian is to pour it into onto a surface and then tilt the surface very, very slightly; if the fluid flows at even a slight angle, then it’s safe to say the fluid behaves like a simple fluid. (Professionals use rheometers). Modeling non-Newtonian fluids is very tricky, and usually done on a case-by-case basis; we won’t concern ourselves here with understanding how they’re modeled, other than that their additional static molecular forces usually depend on the same thing that they depend on for solids, which is some displacement off of an “equilibrium” configuration like in a spring.

The “tilt test”; a puddle of simple fluid in blue and a puddle of non-Newtonian fluid in yellow are placed on a surface. When the surface is tilted slightly, the simple fluid loses its form and flows while the non-Newtonian fluid deforms but doesn’t flow/preserves form.

Looking back at simple fluids, modeling them is far easier than modeling non-Newtonian fluids since the assumption on the behavior of simple fluids induces quite a big constraint on that deviatoric stress tensor \boldsymbol\tau; if we take the assumption that every molecular effect that induces forces in an unmoving fluid is by definition captured through pressure gradients, then we find that \boldsymbol\tau can’t possibly contain any terms involving the physics of the fluid when it’s standing still. It is entirely a function of the properties of the fluid’s motion, which is why I like to call \boldsymbol\tau the hydrodynamic stress tensor when modeling simple fluids.

With that in hand, let’s continue with some extra assumptions. The next assumption we can make is that \boldsymbol\tau, which we’ve said depends on the motion of the fluid, is independent of the absolute velocity of the fluid. This is pretty obvious when you think about it; if the stress of a fluid directly depended on the velocity, we would be able to observe a fluid deform differently depending on whether or not we were walking forwards or backwards relative to the fluid, since velocity is always relative to the observer. Instead, we can assume that \boldsymbol\tau is only a function of velocity gradients, which is almost always true unless you’re moving extremely fast. Mathematically, this means:

\displaystyle \frac{\partial\left(\rho\vec{v}\right)}{\partial t}\ + \nabla \cdot \left(\rho \vec{v} \otimes \vec{v}\right) = \vec{f} - \nabla p +  \nabla \cdot \left(\boldsymbol\tau\!\left[\nabla \vec{v}\right]\right)

Here the brackets just indicate that \boldsymbol\tau is a function of \nabla \vec{v}. But remember; when talking about velocity gradients, we have to consider the gradients in each direction of the velocity components in each direction. So instead of representing three pieces of independent information like \vec{v} does, velocity gradients represents nine, meaning that it’s—you guessed it—a tensor.

As a matter of fact, the assumption above puts a further restriction on \boldsymbol\tau; if I decided to somehow spin myself around a static fluid I’m observing in the lab, I would observe a nonzero velocity gradient even though the fluid is standing perfectly still. As a result, a simple fluid’s deviatoric stress tensor can’t depend on every part of the velocity gradient tensor, just some parts. Some behind-the-scenes tensor manipulations show that in a simple fluid, the deviatoric stress tensor can only depend on the symmetric part of the velocity gradient tensor. This property is commonly referred to as objectivity, and the symmetric part of the velocity gradient tensor is usually called the strain-rate tensor \mathbf{E}, calculated from the velocity gradients by:

\displaystyle \mathbf{E} = \frac{\nabla\vec{v} + \nabla\vec{v}^\intercal}{2}

This still doesn’t really get us anywhere; the nine components of \boldsymbol\tau are so far completely arbitrary functions of the nine components of \mathbf{E}. But we can section simple fluids out into two other categories; fluids for which \boldsymbol\tau is strictly proportional to the strain-rate tensor, which are commonly called Newtonian fluids, and fluids for which it is not, which are called generalized Newtonian fluids. (For what it’s worth, I think the naming conventions for Newtonian/non-Newtonian/generalized Newtonian fluids is confusing and outright ridiculous, but there’s not much I can do about it.) Examples of generalized Newtonian fluids include blood, nail polish, paint, syrup, etcetra. For generalized Newtonian fluids, the behavior of the fluid is properly “liquid”, but the magnitude of the molecular forces induced by the flow of the fluid are complicated functions of the components of the gradients in the velocity. Models for these functions are more of an art than a science, based largely on experimental results; here are some examples.

For a Newtonian liquid, however, we now find that there is a linear relationship between the second-order strain-rate tensor \mathbf{E}, with nine independent components, and the second-order deviatoric stress tensor \boldsymbol\tau, also with nine independent components. This means that we can represent the linear relationship by a fourth-order tensor with 81 different components (yikes!) which I’ll call the viscosity tensor \mathbf{M}:

\displaystyle \boldsymbol\tau = \mathbf{M}:\mathbf{E}

For better or worse, there’s an extra assumption that nearly always applies to fluids of this type; the molecular interactions are independent of direction, and as a result the stress tensor is independent of the orientation of the velocity gradients. Fluids for which this holds true are called isotropic, and fluids for which it’s not are called anisotropic. I’ve never heard of an anisotropic simple Newtonian fluid, but maybe you’ll be the first to find one!

For a simple Newtonian isotropic fluid, doing some behind-the-scenes tensor calculus reveals that our updated conservation of momentum equation has to take the following form once we all these assumptions:

\displaystyle \frac{\partial\left(\rho\vec{v}\right)}{\partial  t}\ + \nabla \cdot \left(\rho \vec{v} \otimes \vec{v}\right) = \vec{f} -  \nabla p\ +\  \mu \nabla^2 \vec{v}\ +\ \left(\lambda + \mu\right)\nabla\left(\nabla\cdot\vec{v}\right)

This absolute mess of an equation is called the compressible Navier-Stokes equation, and is the equation nearly everyone uses when dealing with flow of compressible fluids. The two symbols \mu and \lambda are undetermined proportionality factors that are the only surviving components of that viscosity tensor \mathbf{M}, and are usually determined through experiment; the first is universally referred to as the dynamic viscosity or just viscosity, and the second is so inconsistently defined in the literature that the only name we can give this thing that rheologists wouldn’t argue over is the profoundly unhelpful moniker of second Lamé parameter.

If we further assume that the fluid is incompressible and that the viscosity is the same everywhere, we can polish the equation up a lot more and pull the density out of a lot of expressions, get rid of that last term that depends on \lambda, and wind up with:

\boxed{\displaystyle  \frac{\partial\left(\rho\vec{v}\right)}{\partial  t}\ + \nabla \cdot \left(\rho \vec{v} \otimes \vec{v}\right)  = \vec{f} -   \nabla p\ +\  \mu \nabla^2 \vec{v}}

This is the incompressible Navier-Stokes equation, also known as just the Navier-Stokes equation, and it is the most famous equation in fluid dynamics; understanding it and analyzing its components will be the sole subject of the next Part. Given how involved this sequence of assumptions is for an introduction to fluid mechanics, I’ve put a table below describing how each assumption we made brought us from the Cauchy momentum equation at the beginning of this Part to the Navier-Stokes equation above.

AssumptionMathematical ConsequenceWhen It’s Wrong
Fluctuation hypothesis\vec{r}_{mol} = -\nabla \cdot \left(\rho \overline{\vec{v}' \otimes \vec{v}'}\right)“Exotic” molecular interactions
Simplicity \vec{r}_{mol} = -\nabla p + \nabla \cdot \boldsymbol\tau Material has solid-like properties, pressure defined differently
Objectivity\boldsymbol\tau =  \boldsymbol\tau[\mathbf{E}]Velocities/rotations are extremely fast
Linearity\boldsymbol\tau = \mathbf{M}:\mathbf{E}Molecular forces have a complicated relationship with velocity gradients
Isotropy\boldsymbol\tau = 2 \mu \mathbf{E} + \lambda \text{Tr}(\mathbf{E})\mathbf{I}I don’t know (molecules would likely have to be very asymmetric)
Incompressibility\boldsymbol\tau = \mu \nabla \vec{v}Density is very low, fluid speeds are comparable to molecular speed (see Part II)

Things to Think About

1. What physical effects do you think the viscosity and second Lamé parameter quantify? Can you describe flows where only one of the terms is nonzero?

2. Can you think of other ways to test if a fluid is simple or not? How about to test if a fluid is Newtonian or not?

3. If you’re brave enough to handle some tensors, can you find why the form of the deviatoric stress tensor goes from 2 \mu \mathbf{E} + \lambda Tr(\mathbf{E})\mathbf{I} to \mu \nabla \vec{v} when we take the incompressibility assumption? Remember the definition of the strain-rate tensor!

4. What statements can you make about the stress tensor based just on the fact that \boldsymbol{\sigma} = -\rho \overline{\vec{v}' \otimes \vec{v}'}?

4. What would the Navier-Stokes equation look like if the fluid was a simple Newtonian isotropic fluid but with a non-uniform viscosity?

5. Let’s say that you could describe the forces per unit area in the definition of the stress as coming from an energy density gradient. How would you define the stress in that case, and what would be the associated directions?

Intermezzo: Tensors

Before we continue on with our study of fluid mechanics, it’s best to make a quick pit stop and make sense of that strange concept I mentioned previously; that of a tensor, which we saw in the form of that strange multiplication \rho \vec{v} \otimes \vec{v}. The “proper” way of defining a tensor is hotly debated by mathematicians and physicists of different fields, so my idea here is not necessarily to rigorously define them, but to give you a brief sense of what a tensor is and how we use them in fluid mechanics. In my experience, that is best illustrated by revisiting the concept of a vector, like the velocity \vec{v}.

Knowing the velocity of an object tells us two things; the speed of the object, which is just some number like 20 miles per hour, and the direction in which it is traveling. But as anyone that’s said “No, I meant my left” can attest to, direction is relative—whenever we talk about velocity in its most general sense, we usually describe it in terms of three independent components representing the speeds in specific, perpendicular directions that we’ve arbitrarily defined in the 3-dimensional space we live in (usually denoted with x, y, z). Changing those arbitrarily defined reference directions must then necessarily change the expressions for the components, commonly referred to as projections, even though it doesn’t physically change the velocity of the object itself.

A velocity vector, in black, along with its three components (or projections) onto a set of three distinct directions represented by coordinate axes. Each component represents the speed of the object in the direction the component represents. If I spun the axes around, the components would change, but the vector itself wouldn’t.

That being said, velocity doesn’t intrinsically possess three directions—as we stated above, it possesses just one, the direction in which the object is instantly traveling in. With that in mind, a tensor is just a mathematical object that possesses any number of directions. People usually refer to the amount of directions a tensor possesses as its order, and so we can straightforwardly spot that the velocity is a first-order tensor because of its single “natural” direction.

If the velocity is a first-order tensor, then it might seem obvious to state that the mathematical object \rho \vec{v} \otimes \vec{v}, which we said represented the momentum density flux, possesses two directions, and is as such a second-order tensor. In this case, those directions are straightforward to spot; the direction of the momentum density \rho \vec{v}, and the direction of the flux, which is just the velocity \vec{v}. As a matter of principle, you can always spot the directions of any tensor defined through a sequence of outer products of vectors by looking at the direction of each vector. In fact, both of those directions are the same in this specific case, since \rho \vec{v} \otimes \vec{v} and \vec{v} point in the same direction!

So now that we know what that tensor means, we can think a bit more concretely about how to represent it mathematically. Well, just like the velocity can be represented using three components with reference to some arbitrary coordinate system, the momentum density flux is the product of every component of the 3-D vector \rho \vec{v} with every component of the 3-D vector \vec{v}, so that the tensor \rho \vec{v} \otimes \vec{v} is represented by 3\times3 = 9 independent components. To make sure we distinguish which direction is which, mathematicians usually visually represent second-order tensors using a visual “square” of numbers called a matrix.

These are the nine different components of the momentum density flux tensor, associating every component of the momentum density with every component of the flux. They’ve been placed in a geometric array called a matrix that highlights the different directions associated with each part of the momentum density flux tensor.

This isn’t the end of the story on tensors, though; remember that the Cauchy momentum equation we saw before includes the term \nabla \cdot \left(\rho \vec{v} \otimes \vec{v}\right), not just the momentum density flux. Hence, we need to be able to understand the effects of trying to take derivatives, like the divergence and gradient, on tensors. To do that, let’s repeat our formula on testing it out on the velocity—which we know is a first-order tensor—and see how to extrapolate from there.

The divergence of the velocity at a point, from a purely mathematical standpoint, tells you what the rate of change of the velocity is at that point as you move in each of those predefined spatial coordinates, and then adds all of those rates of changes up into a single direction-less number that tells you if the flow is “accumulating” or “diminishing” at that point. Make sure to notice that the divergence of the velocity is truly directionless; if you looked at a flow upside down, the direction of the flow would certainly change, but the amount by which the flow is “accumulating” at a point doesn’t. This is easily extrapolated to a general, tensorial case; the divergence of a tensor takes the spatial rate of change of the tensor in the direction of one of the tensor’s directional elements, usually defined to be the “last” element, and adds them up. This causes a tensor to “lose” a direction when acted upon by the divergence—meaning that the second-order tensor \rho \vec{v} \otimes \vec{v}, when acted on by the divergence, gives you a first-order tensor/vector, just as the momentum transport equation required.

The divergence of the momentum density flux tensor (order 2). Here, we take the spatial rate of change of each of the 9 elements in the flux direction, and add them up to obtain 3 elements representing a vector (order 1). Note the directions of the partial derivatives match the directions of the flux velocities.

Lastly, let’s take a look at the gradient of the velocity \nabla \vec{v}. Here we are still obtaining the spatial rates of change of the velocity, but this time we’re not adding anything up; in fact, we’re obtaining the spatial rate of change of every component of the velocity in each and every one of those arbitrarily defined directions of space. 3 different possible rates of change for each of the three components indicates that the velocity gradient possesses 9 distinct pieces of information, revealing that the gradient of the velocity is a second-order tensor. As such, we can extrapolate to the general case and state that the gradient of a tensor is another tensor of a single higher order.

The velocity gradient tensor. The two “directions” associated with this tensor are connected to the velocity and to the direction of the rate of change.

Explicitly defining the direction that gets added as a result of taking the gradient is trickier to spot than in the outer product case, where the directions of the resulting tensor are just the directions of each component of the tensor; but the direction exists, and the vector representing it can be found by solving the following equation for the vector \vec{m}:

\nabla \vec{v} \cdot \vec{m} = \vec{v}

(The dot product here is just representing standard matrix-column vector multiplication.) This vector doesn’t really have a name, nor a more intuitive description other than the above as far as I can tell; but I call it the gradial vector for book-keeping, and would love to find a nice application for it.

These explanations are by no means rigorous or comprehensive; there are many incredible and beautiful aspects to the study of tensors that I sadly have to exclude for the sake of coherence, but I invite you to read up on them if you’re interested in the subject. My current favorite book on the topic is Henry Block’s Introduction to Tensor Analysis, which is tragically out-of-print but has been preserved electronically by the thoughtfulness of faculty members in the Theoretical & Applied Mechanics field at Cornell. In any case, this is just about everything we need to know about tensors to get a handle on all the mathematical machinery abound in introductory fluid mechanics; and although the principal language of fluid mechanics is vector calculus, it’s good to get a handle on tensors to deal with an often-ignored part of the world of fluid mechanics—the process by which molecules exert forces on each other when a fluid moves, known as the study of rheology.

Things to Think About

1. What is the difference between a second-order tensor and a matrix? Are all matrices second-order tensors? Are all second-order tensors matrices?

2. If a fluid is incompressible, how does the expression for \nabla \cdot \left(\rho \vec{v} \otimes \vec{v}\right) simplify? Remember the incompressibility equation, \nabla \cdot \vec{v} = 0.

3. How would you describe the gradient of a second-order tensor? Extrapolate from the principles shown here for the velocity gradient.

4. Consider the gradient of any simple scalar (order 0 tensor) function. If you changed around the coordinate axes, how would the gradient change? Would it change in a different way than a velocity vector would?

5. What happens if you take the divergence of the gradient of the velocity? How is that different from taking the gradient of the divergence of the velocity? The first expression is called the vector Laplacian, and will feature prominently in the following sections.