## 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.

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.

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.

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.

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).

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.

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:

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.

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.

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.

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=0}^{n=\infty} \sum\limits_{m=-n}^{n} \left[ \frac{(n+3)r^2\nabla\left[A_{n m}r^n\, Y^m_{n}(\theta,\varphi)\right]}{2\mu(n+1)(2n+3)} -\frac{n\vec{r} A_{n m}r^n\, Y^m_{n}(\theta,\varphi)}{\mu(n+1)(2n+3)}\right] + \sum_{n=1}^{n=\infty} \sum\limits_{m=-n}^{n} \left[(n-2)r^2\nabla\frac{B_{n m}\, Y^m_{n}(\theta,\varphi)}{2r^{n+1}\mu n (1-2n)} +\frac{\left(n+1\right)\vec{r} B_{n m}\, Y^m_{n}(\theta,\varphi)}{r^{n+1}\mu n(2n-1)}\right]+ \frac{B_{0m}\vec{r}}{2r} + \sum_{n=-\infty}^{n=\infty} \sum\limits_{m=-n}^{n} \nabla \left[\left(C_{n m}r^n + \frac{D_{n m}}{r^{n+1}}\right) Y^m_{n}(\theta,\varphi)\right]+ \sum_{n=-\infty}^{n=\infty} \sum\limits_{m=-n}^{n} \nabla \times \left[\vec{r} \left(E_{n m}r^n + \frac{F_{n m}}{r^{n+1}}\right)\, Y^m_{n}(\theta,\varphi) \right]$

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.

*Technically, the $B_{0m}$ term always has to be 0 since the corresponding flow fails to satisfy conservation of mass due to a mathematical quirk in vector calculus.

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:

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.

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.

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.

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:

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:

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.

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.

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.

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):

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.

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:

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.

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 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 }$.

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:

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.

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.

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:

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:

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:

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.

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$:

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.

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.

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.

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.

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.

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.

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.

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.

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.

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.

## Part VI: Transport Theory

The control volume analysis technique we saw and used previously is both astoundingly general and eminently pragmatic. But it has a shortcoming; namely, that we have to define a specific volume before we can make statements about the fluid inside of it. This means that all of the statements we can make about fluids using control volume analysis depend on the specific system we’re analyzing. If we want to make universal statements about the way fluids behave, i.e. describe the physical laws of fluid mechanics, we need to overcome this obstacle. As it turns out, this is easier than it seems.

Let’s briefly summarize what control volume analysis dictates about some quantity in a fluid volume; the rate of change of the total amount of that quantity inside a volume is equal to the total amount of that quantity being generated or destroyed within it, plus the amount of that quantity which is moving in/out of the volume through its boundary surface. If we refer to the quantity we are interested in as $\Xi$, and its associated density as $\xi$, we can generate the following equation for the rate of change of $\Xi$:

$\displaystyle \frac{\partial }{\partial t}\int_V \xi\ dV = \int_V \frac{d \xi}{d t}\ dV - \int_A \left(\vec{j}_\xi \cdot \hat{n}\right)\ dA$

There’s a neat mathematical trick we can do to turn that last surface integral into a volume integral called the divergence theorem. The details don’t matter too much as long as you trust me; the only important thing is that the term still represents what it always did, flow in/out of the control volume. Applying this trick, we find:

$\displaystyle \frac{\partial}{\partial t}\int_V \xi\ dV = \int_V \frac{d\xi}{d t}\ dV - \int_V \nabla \cdot \vec{j}_{\xi}\ dV$

Take a look at what we’ve managed to get; our equation, which represents the change in the total amount of some arbitrary fluid quantity, is exclusively in terms of adding up tiny contributions of different things within the control volume. So if we shrink the control volume down to a point, we don’t have to bother with the integrations—we can look at an equation that relates those contributions directly! Remember, the concept of the equation is the same, we’re just now using it to look at the way quantities accumulate within points rather than volumes.

To save ourselves considerable confusion from mathematical notation, I’m going to do what nearly every scientist does and relabel the $\frac{d \xi}{d t}$ term to something simpler, like $r$; it still represents the same thing, which is the generation/destruction of $\Xi$ through some physical process going on in that point. Performing this relabeling, and shrinking things down to a point, we obtain:

$\displaystyle \frac{\partial\xi}{\partial t}\ = r - \nabla \cdot \vec{j}_{\xi}$

This type of equation is called a transport equation, and the process of using this type of equation to understand something is called a transport theory. Sometimes physicists call this kind of equation a continuity equation, which is in my opinion silly and confusing. Moving the flux term to the other side of the equation, we get its final version:

$\displaystyle \boxed{\frac{\partial \xi}{\partial t}\ + \nabla \cdot \vec{j}_{\xi} = r}$

This equation is the most important equation in science. I am biased, sure, but you can model nearly everything with some version of it: fluid and solid mechanics, chemical kinetics, heat transfer, most of electromagnetism, general relativity, and a decent slice of quantum mechanics. Even statistical mechanics, if you play your cards right! This equation lets you mathematically understand the movement of anything through points in space, as long as you’re able to understand the ways that your quantity of interest can move through space (its fluxes) and the ways your quantity is generated or destroyed at certain points (its sources and sinks).

Here’s a simple and very useful example. If we wanted to look at setting up an equation for mass, you would first identify the corresponding quantity representing mass per unit volume (density $\rho$) and write a transport equation for it down:

$\displaystyle \frac{\partial \rho}{\partial t}\ + \nabla \cdot \left(\rho \vec{v}\right) = r_{\rho}$

Remember that all fluid mass movement is described by fluid velocity $\vec{v}$, so no need to include a diffusion term. In addition, in the absence of something goofy like nuclear reactions or relativistic phenomena, we know mass can’t be destroyed or created. That means that the $r_{\rho}$ term representing the generation and destruction of mass has to be zero! Therefore, we get:

$\displaystyle \frac{\partial \rho}{\partial t}\ + \nabla \cdot \left(\rho \vec{v}\right) = 0$

This equation describes the conservation of mass, and is the second-most important equation in fluid mechanics. If we additionally assumed that the liquid was incompressible, something interesting happens; the density can’t change with respect to time or space, so any term proportional to that needs to vanish in the above equation. Expanding it out and getting rid of the density change terms, we get:

$\displaystyle \frac{\partial \rho}{\partial t}\ + \nabla \rho \cdot \vec{v} + \rho \left(\nabla \cdot \vec{v}\right) = 0$

$\displaystyle \nabla \cdot \vec{v} = 0$

This equation refers to conservation of mass when a liquid is incompressible, and is sometimes called the incompressibility equation.

Even though mass in general is conserved, there might be some chemical reactions going on in our system of interest that change mass from one type of chemical to another. In that situation, you would set up transport equations for every chemical of interest in your system, and wind up with a system of equations for the chemical densities:

$\displaystyle \frac{\partial \rho_i}{\partial t}\ + \nabla \cdot \left(\rho_i \vec{v}\right) - D_i \nabla \rho_i = r_{i}$

where the subscripts indicate a specific chemical in the system, and each source/sink term representing chemical reactions can depend on every other chemical density in the system. Such a system of equations is usually referred to as a reaction network. In chemical systems where there’s no fluid flow, only diffusion and reactions occur, which are studied under the (apt) name of reaction-diffusion systems. Most chemical engineers ignore both the movement of fluid within the system and gradients in the densities of the chemicals, leaving only the reactions to drive the physical changes; this approximation is commonly called the CSTR approximation.

We can construct transport equations for energy as well; one very useful type of energy transport equation involves looking at the movement of thermal energy through a solid. The generic version of the equation should look something like this, according to transport theory:

$\displaystyle \frac{\partial e}{\partial t}\ + \nabla \cdot \vec{j}_{e} = r_e$

As it turns out, the heat energy density of a point is equal to the temperature $T$ of the point, times the mass density $\rho$ times the heat capacity $c$. In addition, because the solid doesn’t flow, the only way energy moves through the solid is diffusively. As a result, we update the equation above to find:

$\displaystyle \frac{\partial e}{\partial t}\ + \nabla \cdot \left[-D_e \nabla e\right] = r_e$

$\displaystyle \frac{\partial \left(c \rho T\right)}{\partial t}\ + \nabla \cdot \left[-D_e \nabla \left(c \rho T \right)\right] = r_e$

Most thermal engineers like to assume the material properties of a solid that’s transferring heat stay constant and are the same everywhere, so we can pull those properties (the heat capacity, density, and heat diffusivity) out of the gradients and simplify:

$\displaystyle c \rho \frac{\partial T}{\partial t}\ - c \rho D_e \nabla \cdot \nabla T = r_e$

$\displaystyle \frac{\partial T}{\partial t}\ - \frac{D_e}{c \rho}\left(\nabla \cdot \nabla T\right) = \frac{r_e}{c \rho}$

Cleaning up the vector calculus term, expressing $\frac{D_e}{c \rho}$ as a single term $\alpha$ representing the thermal diffusivity, and expressing $\frac{r_e}{c \rho}$ as a single thermal generation term $h$, we find:

$\displaystyle \frac{\partial T}{\partial t}\ - \alpha \nabla^2 T = h$

Congrats! You just derived the heat equation using transport theory. Notice how many assumptions we had to make! And funnily enough, that heat diffusivity term $D_e$ we saw above is almost always never called that; it’s confusingly referred to as the thermal conductivity.

Hopefully this has convinced you of the generality and usefulness of transport theory. And with transport theory, we can now construct the most important equation in fluid mechanics; the transport equation for momentum.

If transport theory is to be believed, all I have to do is identify the momentum density, which is $\rho \vec{v}$, and write a transport equation for it that looks like this:

$\displaystyle \frac{\partial\left(\rho \vec{v}\right)}{\partial t}\ + \nabla \cdot \vec{j}_{\rho \vec{v}} = r_{\rho \vec{v}}$

So far, so good. But if I go and write the advective flux term explicitly, what we find seems a bit confusing:

$\displaystyle \frac{\partial (\rho \vec{v})}{\partial t}\ + \nabla \cdot \rho \textrm{?`}\vec{v} \vec{v}\,\textrm{?} + \nabla \cdot \vec{j}^{\left[\rho \vec{v}\right]}_{\text{diff}} = r_{\rho \vec{v}}$

How should we multiply these two vectors? It shouldn’t be a dot product, because then we’d get a scalar, and we can’t take the divergence of that. As it turns out, the right way to multiply these two vectors is with something called an outer product, which looks like this: $\otimes$. And the outer product of two vectors is a weird mathematical object called a tensor, which we will see a lot of in the next section.

Equipped with this knowledge, we rework the momentum transport equation to:

$\displaystyle \frac{\partial (\rho \vec{v})}{\partial t}\ + \nabla \cdot \left(\rho \vec{v} \otimes\vec{v}\right) + \nabla \cdot \vec{j}^{\left[\rho \vec{v}\right]}_{\text{diff}} = r_{\rho \vec{v}}$

It might make you uncomfortable to have that tensor expression up there; most students don’t learn how to do calculus with tensors until after fluid mechanics (if they ever do). Luckily, we can break the term up into terms that can be handled entirely with vector calculus, leading us to the Lamb form of the momentum transport equation:

$\displaystyle \frac{\partial (\rho \vec{v})}{\partial t}\ + \frac{1}{2} \nabla \left(\rho \vec{v} \cdot \vec{v}\right) + \left(\nabla \times \rho \vec{v}\right)\times \vec{v} + \nabla \cdot \vec{j}^{\left[\rho \vec{v}\right]}_{\text{diff}} = r_{\rho \vec{v}}$

This is helpful to understand what the tensor term in the equation means, but I personally don’t find it helpful for keeping track of what the equation means; which is that it’s a transport equation for momentum in a fluid. We’ll revisit the Lamb form later, but right now I’d like to stick with the $\nabla \cdot \left(\rho \vec{v} \otimes\vec{v}\right)$ notation for conceptual clarity.

At this point, it seems tempting to simply say that the only source or sink of momentum in a fluid point comes from body forces that are external to the fluid. But we know that’s not true! We know pressure gradients cause changes in momentum within a fluid, and the momentum for that is coming entirely from the “storage” of molecular energy within the fluid. In addition, we’re missing another huge factor; friction. With friction, the opposite happens—energy in the movement of a macroscale fluid gets dissipated into microscale molecular energy of random motion. The common factor between pressure-induced forces and friction is that both of them are causing macroscale, fluid effects as a result of microscale molecular phenomena.

Noting that momentum density diffusion is also a molecular effect, I find it helpful to lump in all the molecular effects together into a simple term $\vec{r}_{\textrm{molecular}}$ and rewrite the momentum transport equation in the following way:

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

Although it may not look like it yet, this is the most general form of the mathematical basis for essentially all of fluid mechanics: the Cauchy momentum equation.