1.2.1 Particle paths

A particle path is the path, \(\mathbfit{x}(t)\), of a fluid particle over a given time interval. At any time \(t\), the velocity of this particle is given by the velocity field, \(\mathbfit{u}\), evaluated at the position of the particle, meaning we need to solve the vector ODE \[\begin{equation} \mathchoice{\frac{{\mathrm d}\mathbfit{x}}{{\mathrm d}t}}{{\mathrm d}\mathbfit{x}/{\mathrm d}t}{{\mathrm d}\mathbfit{x}/{\mathrm d}t}{{\mathrm d}\mathbfit{x}/{\mathrm d}t}(t) = \mathbfit{u}\big(\mathbfit{x}(t),t\big), \label{particle-path} \end{equation}\] subject to some initial condition (paths have beginnings, after all!). In 3D Cartesian components, with \(\mathbfit{x} = (x,y,z)\) and \(\mathbfit{u} = (u,v,w)\), this is solving \[\mathchoice{\frac{{\mathrm d}x}{{\mathrm d}t}}{{\mathrm d}x/{\mathrm d}t}{{\mathrm d}x/{\mathrm d}t}{{\mathrm d}x/{\mathrm d}t} = u, \quad \mathchoice{\frac{{\mathrm d}y}{{\mathrm d}t}}{{\mathrm d}y/{\mathrm d}t}{{\mathrm d}y/{\mathrm d}t}{{\mathrm d}y/{\mathrm d}t} = v, \quad \mathchoice{\frac{{\mathrm d}z}{{\mathrm d}t}}{{\mathrm d}z/{\mathrm d}t}{{\mathrm d}z/{\mathrm d}t}{{\mathrm d}z/{\mathrm d}t} = w.\]

1.2.2 Streaklines

A streakline is a curve made up of all fluid elements that have passed through a given point in the past. It’s a snapshot at a given time.

Consider smoke from a chimney where the wind initially blows west and later shifts to the east. All the smoke passes through the chimney at \(\mathbfit{a}\), so the streakline at any given time is exactly the trail of smoke. But think about smoke particles emitted at the beginning, middle and end of our story: their particle paths do not always coincide with the streakline:

Because they’re easy to create, streaklines are often used in experiments to visualise flow.

To determine a streakline at a particular time, we just need to find the particles that previously passed through, or were released from, some given point, \(\mathbfit{a}\). That’s the same problem as finding particle paths, but with a different initial condition: if we name the release time \(\tau\) (different for each particle), we have to solve \[\mathchoice{\frac{{\mathrm d}\mathbfit{x}}{{\mathrm d}t}}{{\mathrm d}\mathbfit{x}/{\mathrm d}t}{{\mathrm d}\mathbfit{x}/{\mathrm d}t}{{\mathrm d}\mathbfit{x}/{\mathrm d}t} = \mathbfit{u}(\mathbfit{x},t), \quad \mathbfit{x}(\tau) = \mathbfit{a}.\]

Key difference:

• A particle path describes one particle over lots of time.

• A streakline describes lots of particles at one time.

1.2.3 Streamlines

A streamline is a line parallel to the velocity vector \(\mathbfit{u}(\mathbfit{x}, t)\) at a fixed time, \(t=t_0\). It is often used as a way to visualise a snapshot of the velocity field everywhere, because you can draw a streamline through any point, \(\mathbfit{x}_0\). A streamline is a curve \(\mathbfit{x}(s)\) satisfying \[\mathchoice{\frac{{\mathrm d}\mathbfit{x}}{{\mathrm d}s}}{{\mathrm d}\mathbfit{x}/{\mathrm d}s}{{\mathrm d}\mathbfit{x}/{\mathrm d}s}{{\mathrm d}\mathbfit{x}/{\mathrm d}s} = \mathbfit{u}\bigl(\mathbfit{x}(s), t_0\bigr),\] where \(t_0\) is fixed and \(s\) is a parameter that varies along the streamline. To obtain the streamline passing through \(\mathbfit{x}_0\), we also specify the condition \(\mathbfit{x}(0) = \mathbfit{x}_0\).

1.5.1 Stream functions to measure flux

The stream function has an additional physical interpretation. Take two points on the same streamline and draw a curve, \(C\), between them. What is the total flux through this curve?

We’ve seen flux already: it’s the total amount of fluid passing through the curve, \[\int_C \mathbfit{u} \cdot \widehat{\mathbfit{n}} \, {\mathrm d}s.\] No flow can pass through the streamline, because as we’ve seen already, the streamline is parallel to the velocity. So there’s only one way in and out of the region enclosed by the streamline and \(C\): by mass conservation, the total flow in through \(C\) must equal the total flow out through \(C\), i.e. \[\int_C \mathbfit{u}\cdot\widehat{\mathbfit{n}}\, {\mathrm d}s = 0.\]

So that’s the case when \(C\) goes between two points on the same streamline. Now what about a curve between two points on two different streamlines: what’s the flux through \(C\) this time? We again will calculate \(\oint \mathbfit{u}\cdot\widehat{\mathbfit{n}} \, {\mathrm d}s\) but we should first observe that flux can go in either direction through \(C\): ‘forwards’ or ‘backwards’. Whichever way we pick \(\widehat{\mathbfit{n}}\) will represent positive flux.

Let \(C\) be any simple curve connecting \(\mathbfit{x}_1\) to \(\mathbfit{x}_2\), with tangent vector \(\mathbfit{t}(s)\), \[\mathbfit{t} = \begin{pmatrix} t_x \\ t_y \end{pmatrix} .\] The unit normal to \(C\), may therefore be written as \[\widehat{\mathbfit{n}} = \frac{1}{|\mathbfit{t}|} \begin{pmatrix} t_y \\ -t_x \end{pmatrix} .\] This points down in the diagram, but more generally, it points clockwise relative to \(\mathbfit{x}_1\). The flux through \(C\) in this direction is \[\begin{aligned} \int_C\mathbfit{u}\cdot\widehat{\mathbfit{n}}\,{\mathrm d}s &= \int_C \begin{pmatrix} \mathchoice{\frac{\partial\psi}{\partial y}}{\partial\psi/\partial y}{\partial\psi/\partial y}{\partial\psi/\partial y} \\ - \mathchoice{\frac{\partial\psi}{\partial x}}{\partial\psi/\partial x}{\partial\psi/\partial x}{\partial\psi/\partial x} \end{pmatrix} \cdot \frac{1}{|\mathbfit{t}|} \begin{pmatrix} t_y \\ -t_x \end{pmatrix} {\mathrm d}s\\ &= \int_C\frac{\bm{\nabla}\psi\cdot\mathbfit{t}}{|\mathbfit{t}|}\,{\mathrm d}s \\ &= \int_C\bm{\nabla}\psi\cdot{\mathrm d}\mathbfit{x} = [\psi]_C = \psi(\mathbfit{x}_2)-\psi(\mathbfit{x}_1), \end{aligned}\] which is path independent.

We can conclude: the flux (in a clockwise direction about \(\mathbfit{x}_1\)) across any line between two points \(\mathbfit{x}_1\) and \(\mathbfit{x}_2\) is given by \(\psi(\mathbfit{x}_2)-\psi(\mathbfit{x}_1)\).

This has the consequence that if streamlines are drawn at regular intervals of \(\psi\), then just like how close contours of a map represent steep slopes, closer lines represent faster flow. (Have a think why if this isn’t immediately obvious.)

1.7.1 Vorticity

We have just seen the introduction of the vorticity, \[\mathbfit{\omega} = \bm{\nabla}\times\mathbfit{u}.\] Since it is twice the angular velocity of the fluid element, it is a measure of local rotation: how the fluid is swirling at each point in space. Flows where \(\bm{\nabla}\times\mathbfit{u}=\mathbf{0}\) everywhere are called irrotational.

So we have a rotating flow, and we have vorticity. But these do not necessarily coincide! We have to be careful not to judge vorticity – which is local rotation at every point – based on ‘overall’ rotation vibes. Here are a couple of illustrative examples:

The way to picture vorticity is to think of little paper windmills thrown into the flow:

If the windmills rotate, there is vorticity. In this case, the flow above each windmill in the top half-plane is faster than the flow beneath it, so naturally the windmills rotate clockwise.

The line vortex’s property of ‘high spin near the centre, not much spin away from the centre’ makes it a simple model for tornadoes, whirlpools and slightly more mundane things like emptying a bath. We will demonstrate the destructive power of tornadoes shortly when we prove an important result about vorticity conservation.

1.7.2 Circulation

These curious examples highlight that vorticity is local. If we want to quantify ‘how much rotation’ there is inside a region, vorticity alone is the wrong measurement. We need circulation.

Let \(C\) be a simple, smooth, closed curve in the fluid. Then the circulation around \(C\) is \[\begin{equation} \mathit{\Gamma} = \oint_C \mathbfit{u}\cdot{\mathrm d}\mathbfit{x}, \label{circulation-line-integral} \end{equation}\] that is, it counts how much \(\mathbfit{u}\) is in the direction of the curve, as we go around the curve.

lots of circulation	no circulation

If \(S\) denotes a surface spanning \(C\) (with its normal pointing to match the right-hand rule on \(C\)) then Stokes’ theorem gives \[\begin{align} \mathit{\Gamma} = \oint_C \mathbfit{u}\cdot{\mathrm d}\mathbfit{x}& = \int_S\bm{\nabla}\times\mathbfit{u}\cdot{\mathrm d}\mathbfit{S}\\& = \int_S\mathbfit{\omega}\cdot{\mathrm d}\mathbfit{S}. \label{circulation-surface-integral} \end{align}\] So the circulation is a measure of the overall rotation of the fluid in a certain region.

Although we might think this singularity is nasty, actually it’s quite neat: we have seen that we can represent a circulation over a whole region by a single delta function – a line vortex. We will see shortly that there are other types of singularities – those representing sources and sinks – and we can build interesting flows from them.

1.7.3 Vortex tubes

The line vortex is a pretty good model of real-life vortices if we zoom out a bit, since it is capable of producing lots of circulation, and typically most vorticity is near the centre of a real-life vortex anyway. But up close, ‘near’ can still be quite large: real-life vortices like tornadoes or whirlpools really do have vorticity within a reasonable distance of the centre. In tornadoes, this high-vorticity area coincides (due to a pressure drop and subsequent condensation) with the cloud you can see:

And of course, it would ultimately be very strange if all real-world rotating flows happened to be described by the line vortex’s velocity, \(\mathbfit{u}=(k/r)\widehat{\mathbfit{e}}_\theta\). We’ve already seen that commonplace rigid body rotation is described by a different equation, \(\mathbfit{u} = kr\widehat{\mathbfit{e}}_\theta\).

We can learn more by instead considering our vortex to have thickness. We can draw a curve in the fluid, like we did for circulation earlier, and ask how the circulation inside it is conserved throughout the flow.

So let’s take that curve, and again consider the vorticity at each point. But this time, we draw a line parallel to the vorticity at each point on the curve: just like how we defined streamlines as being parallel to velocity, but this time the lines are parallel to vorticity instead. If we consider our rigid body planar flow from before, for example, where \(\mathbfit{\omega}(\mathbfit{x})=2k\widehat{\mathbfit{e}}_z\), we get the world’s worst birthday cake:

We can keep drawing these lines, making them longer so that as we follow them, they are always parallel to \(\mathbfit{\omega}(\mathbfit{x},t)\) at every \(\mathbfit{x}\) at a fixed time: again, exactly like streamlines but this time for vorticity. These lines are called vortex lines – not to be confused with the line vortex we saw earlier – or sometimes vortex filaments. Taking the bundle of lines that go through \(C\), we end up with a wall of lines called a vortex tube.

What might one of these look like for a 3D flow?

The circulation around our original \(C\) turns out to be the same for all closed curves \(C\) around the tube, allowing us to define the vortex tube’s strength by this circulation, \(\mathit{\Gamma}\).

To confirm this, let’s pick two closed curves, \(C_1\) and \(C_2\), around the same vortex tube:

Labelling the surfaces at the ends of the tube \(S_1\) and \(S_2\), and the tubular surface by \(S_3\), then we can consider the difference between the two circulations \[\begin{aligned} \mathit{\Gamma}_2 - \mathit{\Gamma}_1 &= \int_{S_2}\mathbfit{\omega}\cdot{\mathrm d}\mathbfit{S} - \int_{S_1}\mathbfit{\omega}\cdot{\mathrm d}\mathbfit{S} \\ &= \int_{S_2}\mathbfit{\omega}\cdot\widehat{\mathbfit{n}}\,{\mathrm d}{S} + \int_{S_1}\mathbfit{\omega}\cdot\widehat{\mathbfit{n}}\,{\mathrm d}S \quad \text{because we want $\widehat{\mathbfit{n}}$ to point outward}\\ &= \int_{S_2}\mathbfit{\omega}\cdot\widehat{\mathbfit{n}}\,{\mathrm d}{S} + \int_{S_1}\mathbfit{\omega}\cdot\widehat{\mathbfit{n}}\,{\mathrm d}S + \int_{S_3}\underbrace{\mathbfit{\omega}\cdot\widehat{\mathbfit{n}}}_{=0}{\mathrm d}S \quad \text{because $\widehat{\mathbfit{n}} \perp \mathbfit{\omega}$ on $S_3$}\\ &= \int_{S_1+S_2+S_3}\mathbfit{\omega}\cdot\widehat{\mathbfit{n}}\,{\mathrm d}{S}\\ &= \int_V \bm{\nabla}\cdot\mathbfit{\omega} \, {\mathrm d}V \qquad \text{by the divergence theorem}\\ &= 0 \qquad \text{because $\bm{\nabla}\cdot\mathbfit{\omega}\equiv0$ always.} \end{aligned}\] So the two arbitrarily-chosen circulations are the same, and hence the strength of the vortex tube is well defined.

We’ve used the fact that vorticity is divergence-free here: you can convince yourself this is true since \[\bm{\nabla}\cdot\mathbfit{\omega} = \bm{\nabla}\cdot(\bm{\nabla}\times\mathbfit{u}) \equiv 0,\] using the ‘div curl = 0’ identity at the end there, which Wikipedia reminds us is ‘a special case of the vanishing of the square of the exterior derivative in the De Rham chain complex’ (obviously).

1.7.4 Flow reconstruction

As vortex rings demonstrate, it can often be easier to think about fluid flow in terms of the vorticity structures present rather than the flow itself. This leads us to wonder: if we prescribe \(\mathbfit{\omega}\), can we always recover \(\mathbfit{u}\)?

If the fluid is incompressible, \(\bm{\nabla}\cdot\mathbfit{u}=0\), then there exists a vector function \(\mathbfit{A}\) such that \[\mathbfit{u} = \bm{\nabla}\times\mathbfit{A}.\] This is like what we did when we derived the stream function (there, \(\mathbfit{A}\) was \(\psi \widehat{\mathbfit{e}}_z\)). The vorticity is therefore \[\mathbfit{\omega} = \bm{\nabla}\times\mathbfit{u} = -{\nabla^2}\mathbfit{A},\] which is Poisson’s equation for each component of the vector equation. If you’re taking PDEs, you’ll see that we can solve this by using the Green’s function, which gives us \[\mathbfit{A} = \frac{1}{4\pi}\int_{V} \frac{\mathbfit{\omega}(\mathbfit{x}')}{|\mathbfit{x}-\mathbfit{x}'|} {\mathrm d}V',\] and so the velocity is \[\mathbfit{u} = \bm{\nabla}\times \mathbfit{A} = \frac{1}{4\pi}\int_V \frac{\mathbfit{\omega}(\mathbfit{x}')\times(\mathbfit{x}-\mathbfit{x}')}{|\mathbfit{x}-\mathbfit{x}'|^3}{\mathrm d}V' =:\mathbfit{u}_v,\] where we have been able to happily bring the curl inside the integral because \(\bm{\nabla}\) is \({\mathrm d}\mathbfit{x}\), whereas the integral is \({\mathrm d}\mathbfit{x}'\).

This looks good – problem solved? Ah, but wait! Notice that \[\mathbfit{u} \quad \text{and} \quad \mathbfit{u} + \bm{\nabla}\phi,\] for some \(\phi\), have the same curl, since \[\bm{\nabla}\times\bm{\nabla}\phi \equiv 0 \quad \text{for all }\phi.\] This additional component, \(\bm{\nabla}\phi\), is called a potential flow with its defining characteristic being that it has no vorticity. We will take a look at potential flows properly in the next section.

The general solution is therefore \[\mathbfit{u} = \bm{\nabla}\phi + \mathbfit{u}_v.\] This is the Biot–Savart law. Where \(\mathbfit{u}_v\) carries the vorticity information, \(\phi\)’s job is to sort the boundary conditions. Setting the boundary conditions in your problem will set the conditions on \(\phi\), which solve the problem uniquely.

1.8.1 Solutions to Laplace’s equation

So long as we have compatible boundary conditions, we can solve Laplace’s equation using separation of variables. This is true for all dimensions and for any coordinate system, but we’re going to look at 2D polar coordinates. Then, if \(\phi(r,\theta)=R(r)\mathit{\Theta}(\theta)\), \[\begin{aligned} {\nabla^2}\phi = \frac{1}{r}\mathchoice{\frac{\partial}{\partial r}}{\partial/\partial r}{\partial/\partial r}{\partial/\partial r}\left(r \mathchoice{\frac{\partial\phi}{\partial r}}{\partial\phi/\partial r}{\partial\phi/\partial r}{\partial\phi/\partial r}\right) + \frac{1}{r^2}\mathchoice{\frac{\partial^2 \phi}{\partial\theta^2}}{\partial^2 \phi/\partial\theta^2}{\partial^2 \phi/\partial\theta^2}{\partial^2 \phi/\partial\theta^2} &= 0 \\ \implies \frac{\mathit{\Theta}}{r}\mathchoice{\frac{{\mathrm d}}{{\mathrm d}r}}{{\mathrm d}/{\mathrm d}r}{{\mathrm d}/{\mathrm d}r}{{\mathrm d}/{\mathrm d}r}\left(r \mathchoice{\frac{{\mathrm d}R}{{\mathrm d}r}}{{\mathrm d}R/{\mathrm d}r}{{\mathrm d}R/{\mathrm d}r}{{\mathrm d}R/{\mathrm d}r}\right) + \frac{R}{r^2}\mathchoice{\frac{{\mathrm d}^2 \mathit{\Theta}}{{\mathrm d}\theta^2}}{{\mathrm d}^2 \mathit{\Theta}/{\mathrm d}\theta^2}{{\mathrm d}^2 \mathit{\Theta}/{\mathrm d}\theta^2}{{\mathrm d}^2 \mathit{\Theta}/{\mathrm d}\theta^2} &= 0 \\ \implies \frac{r}{R}\left(r R'\right)' &= -\frac{\mathit{\Theta}''}{\mathit{\Theta}}. \end{aligned}\] The key idea is that the left-hand side is a function of \(r\) only, and the right-hand side is a function of \(\theta\) only. So if we fix some \(r\) and vary \(\theta\), the right-hand side can’t change (because \(r\) is fixed), and vice-versa. This implies the equation is equal to a constant, \(\lambda\): \[\frac{r}{R}\left(r R' \right)' = -\frac{\mathit{\Theta}''}{\mathit{\Theta}} = \lambda.\] You’ll remember from Calculus I that the solutions to this equation depend on whether \(\lambda\) is positive, zero or negative.

1.8.2 Main characters of potential flow

Line sources and sinks:

The only axisymmetric potential, \(\phi=\phi(r)\), is the solution to [laplace-solution] when \(n=0\), \[\phi(r) = m \log r + C,\] where \(m\) and \(C\) are constants. Just like the stream function, the final constant \(C\) is arbitrary so we can set it to zero. What is the velocity represented here? \[\mathbfit{u} = \bm{\nabla}\phi = \mathchoice{\frac{\partial\phi}{\partial r}}{\partial\phi/\partial r}{\partial\phi/\partial r}{\partial\phi/\partial r}\widehat{\mathbfit{e}}_r+ \frac{1}{r}\mathchoice{\frac{\partial\phi}{\partial\theta}}{\partial\phi/\partial\theta}{\partial\phi/\partial\theta}{\partial\phi/\partial\theta}\widehat{\mathbfit{e}}_\theta= \frac{m}{r}\widehat{\mathbfit{e}}_r.\] This is a line source (\(m>0\)) or a sink (\(m<0\)) with strength \(m\), located at the origin. In fact, you saw this if you did the exercise which included [line-source-exercise], where you worked out that \(2m\pi\) was the flux produced by such a source.

line source	line sink

1.8.3 Flow past a cylinder

Let’s consider a static cylinder of radius \(a\) in a uniform flow with velocity \((U,0)\) at infinity. No flow can get into the cylinder, so we can ask ourselves: what is the velocity everywhere in the flow?

The upstream condition implies that away from the cylinder, we expect \(\mathbfit{u}=(U,0)\) which means \(\phi=Ux\). In polar coordinates, \[\phi(r\to\infty)=Ux=Ur\cos\theta.\] On the surface of the cylinder, we have the no normal flow condition. The normal to a cylinder is the radial direction, so we have \[\bm{\nabla}\phi\cdot\widehat{\mathbfit{n}}\bigr|_{r=a}=\left.\mathchoice{\frac{\partial\phi}{\partial r}}{\partial\phi/\partial r}{\partial\phi/\partial r}{\partial\phi/\partial r}\right|_{r=a}=0.\] What parts of our general mould, [laplace-solution], fits these conditions?

Well, actually, not much! The only bit that fits the form of the \(r\to\infty\) condition is \[\phi(r,\theta) = (A_1r + B_1r^{-1})(C_1\cos \theta),\] and comparing the constants, we see it has to be \[\phi(r,\theta) = (Ur + B_1r^{-1})\cos \theta.\] Now for the condition on the surface of the cylinder, we have that at \(r=a\), \[\begin{aligned} \mathchoice{\frac{\partial\phi}{\partial r}}{\partial\phi/\partial r}{\partial\phi/\partial r}{\partial\phi/\partial r} = (U - B_1a^{-2})\cos\theta &= 0 \\ \implies B_1 &= Ua^2. \end{aligned}\] So our final solution for a cylinder in a uniform flow is \[\begin{equation} \label{cylinder-in-uniform-flow} \phi = U\left(1 + \frac{a^2}{r^2}\right)r\cos\theta. \end{equation}\]

If we work out the velocity, \(\mathbfit{u}=\bm{\nabla}\phi\), we can plot the flow:

(the yellow circles mark stagnation points). Curiously, this also gives a valid flow inside the cylinder: in the context of our problem, this flow can be ignored as it just represents a ‘virtual’ flow which is compatible with the boundary condition on the cylinder.

1.8.4 Flow past a rotating cylinder

Things get more interesting if we now let the cylinder rotate! Recalling again that we can just add solutions to Laplace’s equation, we can add [cylinder-in-uniform-flow] to the line vortex, [line-vortex-potential]: \[\begin{equation} %\label{cylinder-in-uniform-flow} \phi = U\left(1 + \frac{a^2}{r^2}\right)r\cos\theta + \frac{\mathit{\Gamma}}{2\pi}\theta. \end{equation}\] What does this do? Well note that since we added on a term which only involves \(\theta\), the no normal flow condition around the cylinder of \(\mathchoice{\frac{\partial\phi}{\partial r}}{\partial\phi/\partial r}{\partial\phi/\partial r}{\partial\phi/\partial r}=0\) is still satisfied. If you think about what the line vortex is, this shouldn’t be a surprise.

A good way to understand these flows is to look for the stagnation points. Recalling the unfortunate fluid mixing problem in Problems Class 1, stagnation points can mark where fluid separates, and streamlines on either side can head in different directions around an obstacle. We can work out the velocity, \[\begin{aligned} \mathbfit{u} = \bm{\nabla}\phi &= \mathchoice{\frac{\partial\phi}{\partial r}}{\partial\phi/\partial r}{\partial\phi/\partial r}{\partial\phi/\partial r}\widehat{\mathbfit{e}}_r+ \frac1r \mathchoice{\frac{\partial\phi}{\partial\theta}}{\partial\phi/\partial\theta}{\partial\phi/\partial\theta}{\partial\phi/\partial\theta}\widehat{\mathbfit{e}}_\theta\\ &= U\left(1-\frac{a^2}{r^2}\right)\cos\theta \, \widehat{\mathbfit{e}}_r+ \left[\frac{\mathit{\Gamma}}{2\pi r}-U\left(1+\frac{a^2}{r^2}\right)\sin\theta\right]\widehat{\mathbfit{e}}_\theta, \end{aligned}\] and then look for where \(u_r=u_\theta = 0\).

We can see that \(u_r=0\) if \(r=a\) (not surprising!). Then \(u_\theta=0\) if \[\begin{equation} \label{u-theta-zero-gamma} \mathit{\Gamma} = 4\pi U a \sin \theta, \end{equation}\] so our two stagnation points are at \[r=a, \quad \sin\theta = \frac{\mathit{\Gamma}}{4\pi U a}.\] We can mark these and plot the flow. In this picture, we’ve chosen \(\mathit{\Gamma}<0\) (for interesting application reasons coming shortly) and so the stagnation points have \(\theta<0\):

So the cylinder is rotating, bringing most of the flow over the top in the same direction, but some flow underneath still goes against the direction of the cylinder. Do we always expect this to be the case? Surely not if we spin the cylinder fast enough…

You might also smell a rat if you realise that because \(\lvert\sin\theta\rvert\leq 1\), [u-theta-zero-gamma] only has a solution if \(|\mathit{\Gamma}|<4\pi U a\) (assuming \(U>0\)). So while our solution is good for low circulation, if the circulation is too large, something else is happening! What could that be?

The other option for \(u_r=0\) is, of course, when \(\cos\theta=0\), i.e. \(\theta = \pm\pi/2\). Under this condition, setting \(u_\theta=0\) gives us \[r^2 \mp \frac{\mathit{\Gamma}}{2\pi U}r +a^2 = 0.\] The case \(\theta=\pi/2\) gives us the negative branch of this equation. Since \(\mathit{\Gamma}\) is negative, we have all positive coefficients and hence there are no positive real solutions for the radius, \(r\). It is typically good for radii to be real and positive so we can discard this case.

However, \(\theta=-\pi/2\) gives us the positive branch of the equation, leading to the solutions \[r = a\left[\frac{|\mathit{\Gamma}|}{4\pi U a} \pm \sqrt{\left(\frac{\mathit{\Gamma}}{4\pi U a}\right)^2 - 1}\right].\] Given that \(|\mathit{\Gamma}| > 4\pi U a\), we get two real, positive roots: one greater than \(a\) and one less than \(a\) (a general property of expressions of the form \(a[b\pm\sqrt{b^2-1}]\)). We are only interested in the solution for \(r>a\), since that is the one that represents flow outside the cylinder.

All this gives us a flow that looks like:

The stagnation point is now below the cylinder.

What’s going on here? If we are in the frame of the cylinder, then this setup is equivalent to a tennis, cricket or golf ball heading leftward (into the wind), rotating backwards as it moves. This is backspin. Sports fans among you will know that backspin generates lift! This is due to the Magnus effect, and is because of a difference in pressure.

But hey, we haven’t spoken about pressure yet. So let’s get onto that quickly, covering one last thing before we do.

1.8.5 Method of images

In the last examples, we saw how a ‘virtual’ flow inside the cylinder – and therefore outside of the fluid domain – meant that the solution for the fluid potential in the entire plane matched the appropriate boundary conditions.

This was neat because, in general, solving Laplace’s equation in \(\mathbb{R}^2\) is easy, but solving it in a more complicated domain is hard. If we can instead add fictitious sources/sinks/vortices into \(\mathbb{R}^2\) to replicate the effect of a more complicated domain, we make life much easier for ourselves. This is called the method of images.

The classic example is of a line source near a wall.

In fact, the method of images works for more than just a straight wall. If the motion of the fluid in the plane is due to a distribution of singularities – our main characters – and there exists a curve \(C\) in the plane with no flow across it, then the system of singularities on one side of \(C\) is the image of the system on the other side.

The following two diagrams show the image system of a line source inside (i) a quadrant, and (ii) a wedge with angle \(\pi/4\).

Quick question: does this work for a wedge with an arbitrary angle \(\alpha\)? No, because you have to be able to reflect enough times to fill the plane without adding reflected singularities back into your original wedge. Mathematically this means \(\alpha = \pi/n\) for some integer \(n\).

Let’s have a look now at a system involving a line vortex: the idea is the same but we just have to be careful that rotation direction is reflected:

Recall that one reason we are interested in line vortices is because – if you zoom out enough – they’re pretty good approximations to real-life vortices, like hurricanes. And this last example demonstrates that line vortices aren’t fixed in place – they are instead moved by the flow, which may include the effects of other vortices, or just their reflections in some physical boundary.

If you want to predict the movement of hurricanes, then, it might be enough just to track their centres, and model how those move. This motivates the idea of tracking vorticity, which we will get to in the next chapter, once we have decided which forces we want our fluids to experience.

2.1.1 The momentum equation

As we did for conservation of mass, consider again a fixed volume \(V\) of fluid with surface \(S\). The momentum is \[\int_V \rho \mathbfit{u} \, {\mathrm d}V.\] The momentum changes for two reasons. The first is the application of forces – the classical mechanics bit of Newton’s second law. The second reason harks back to our conservation of mass derivation in 1.4. Once again, the volume is fixed, but fluid particles can enter and leave the volume. In doing so, they can add and subtract momentum from \(V\).

Quickly looking back at the conservation of mass derivation, we can note that mass is a scalar, whereas here, momentum is a vector. The same procedure (divergence theorem, product rule) will work if we start with a vector but we’d have to use the tensor version of these tools. That’s OK, but to make life easier we will instead show that the \(i\)th component of momentum, \[\int_V \rho u_i \, {\mathrm d}V\] is conserved (a scalar!), and this of course will imply that the (vector) momentum is also concerned.

The momentum change in \(V\) per unit time, is \[\begin{aligned} \mathchoice{\frac{{\mathrm d}}{{\mathrm d}t}}{{\mathrm d}/{\mathrm d}t}{{\mathrm d}/{\mathrm d}t}{{\mathrm d}/{\mathrm d}t} \int_V \rho u_i \, {\mathrm d}V &= -\int_S \underbrace{\rho u_i \underbrace{\mathbfit{u} \cdot \widehat{\mathbfit{n}}}_{\substack{\text{speed of flow}\\\text{in $\widehat{\mathbfit{n}}$-direction}}} \underbrace{{\mathrm d}S}_{\text{area}}}_{\text{momentum flux through ${\mathrm d}S$}} + \int_V F_i \, {\mathrm d}V \\ &= - \int_S \rho u_i\mathbfit{u}\cdot \,{\mathrm d}\mathbfit{S}+ \int_V F_i \, {\mathrm d}V, \end{aligned}\] remembering that the minus sign on the surface integral is because particles leaving the volume represent momentum loss. We’ll deal with the force bit shortly: \(\mathbfit{F}\) for now is the total force on \(V\) per unit volume.

Because the volume is fixed, we can again differentiate under the integral sign on the left-hand side and use the divergence theorem on the surface integral to get \[\int_V \mathchoice{\frac{\partial}{\partial t}}{\partial/\partial t}{\partial/\partial t}{\partial/\partial t}[\rho u_i] \, {\mathrm d}V = - \int_V \bm{\nabla}\cdot[\rho u_i\mathbfit{u}] \, {\mathrm d}V + \int_V F_i \, {\mathrm d}V.\] We can use the product rule on the left-hand side, and additionally use the product rule for divergence on the right-hand side, after choosing how to bracket the three factors \(\rho u_i \mathbfit{u}\) into two bits: \[\bm{\nabla}\cdot[\rho u_i\mathbfit{u}] = \bm{\nabla}\cdot[u_i (\rho \mathbfit{u})] = \rho\mathbfit{u}\cdot\bm{\nabla}u_i + u_i \bm{\nabla}\cdot(\rho \mathbfit{u}).\] This leaves us with \[\int_V \left[ {\color{orange}\mathchoice{\frac{\partial\rho}{\partial t}}{\partial\rho/\partial t}{\partial\rho/\partial t}{\partial\rho/\partial t}u_i} + \rho \mathchoice{\frac{\partial u_i}{\partial t}}{\partial u_i/\partial t}{\partial u_i/\partial t}{\partial u_i/\partial t} \right] {\mathrm d}V = \int_V \bigl[ - \rho\mathbfit{u}\cdot\bm{\nabla}u_i - {\color{orange}u_i \bm{\nabla}\cdot(\rho \mathbfit{u})} + F_i \bigr] {\mathrm d}V.\] The orange terms on each side cancel from the continuity equation, [continuity-equation-integral], leaving us with \[\begin{aligned} \int_V \rho \mathchoice{\frac{\partial u_i}{\partial t}}{\partial u_i/\partial t}{\partial u_i/\partial t}{\partial u_i/\partial t} {\mathrm d}V &= \int_V \bigl[-(\rho\mathbfit{u}\cdot\bm{\nabla})u_i + F_i \bigr] {\mathrm d}V \end{aligned}\] which we can rewrite with the material derivative, [material-derivative], to find \[\int_V \rho\mathchoice{\frac{\mathrm{D} u_i}{\mathrm{D} t}}{\mathrm{D} u_i/\mathrm{D} t}{\mathrm{D} u_i/\mathrm{D} t}{\mathrm{D} u_i/\mathrm{D} t} {\mathrm d}V = \int_V F_i \, {\mathrm d}V,\] or back in vector form, \[\int_V \rho\mathchoice{\frac{\mathrm{D} \mathbfit{u}}{\mathrm{D} t}}{\mathrm{D} \mathbfit{u}/\mathrm{D} t}{\mathrm{D} \mathbfit{u}/\mathrm{D} t}{\mathrm{D} \mathbfit{u}/\mathrm{D} t} {\mathrm d}V = \int_V \mathbfit{F} \, {\mathrm d}V.\] Surprise! It’s \(\mathbfit{F}=m\mathbfit{a}\)!

Now’s it time to think about what the forces are that act on our fluid volume.

2.1.2 Body forces

As we said at the top of the chapter, we consider that there are two types of forces. Body forces act everywhere on a fluid element, and you will be familiar with them from classical mechanics. It turns out to be somewhat convenient to write body forces instead as accelerations, \[\int_V \mathbfit{F}_\text{body} \, {\mathrm d}V = \int_V \rho \mathbfit{f} \, {\mathrm d}V,\] where \(\mathbfit{f}\) contains these accelerations, but hopefully you can see this is just a notation choice.

A common body force is gravity, \(\mathbfit{f} = -g\widehat{\mathbfit{e}}_z\), but other examples might be magnetic forces if the fluid conducts electricity, or inertial ‘forces’ such as the Coriolis and centrifugal forces when the fluid is in a rotating reference frame (something you’ll see in Geophysical and Astrophysical Fluid Dynamics IV). If \(\mathbfit{f} = \mathbf{0}\) we say the system is unforced.

2.1.3 Pressure

If we consider two touching fluid elements, or fluid volumes, why do they not just merge into each other? At the microscale, fluid molecules are bouncing off each other, all the way along our fluid element’s surface.

This is pressure, and from a fluid element’s perspective, it is pushing inwards (i.e. in the normal direction) all over the surface, \[\int_V \mathbfit{F}_\text{pressure} \, {\mathrm d}V = \int_S -p \widehat{\mathbfit{n}} \, {\mathrm d}S = \int_S -p \, {\mathrm d}\mathbfit{S}.\]

Pressure is serious business! From a problem-solving perspective, \(p(\mathbfit{x},t)\) is initially unknown to us; we have to find it at the same time as finding \(\mathbfit{u}(\mathbfit{x},t)\).

2.1.4 The stress tensor

This term, we are only going to consider body forces and pressure, but we can extend the idea of pressure to set up some notation for next term.

We said that pressure only acts normal to the surface, \[\begin{equation} \int_V \mathbfit{F}_\text{pressure} \, {\mathrm d}V = \int_S -p \, {\mathrm d}\mathbfit{S}, \label{f-pressure} \end{equation}\] but in general, we can have forces acting in any direction on the surface of a fluid element. Friction, for example, acts tangentially to surfaces and we will need that next term to introduce viscosity. We can encode forces acting on any surface by introducing the stress tensor or Cauchy stress tensor, \(\mathsfbfit{\normalsize{𝞼}}\), which is a \(3\times3\) rank two tensor, \[\normalsize{𝞼} = (\sigma_{ij}) = \begin{pmatrix} \sigma_{11} & \sigma_{12} & \sigma_{13}\\ \sigma_{21} & \sigma_{22} & \sigma_{23}\\ \sigma_{31} & \sigma_{32} & \sigma_{33} \end{pmatrix}.\]

Physically, the index \(i\) of \(\sigma_{ij}\) gives the coordinate direction in which the force acts, and the index \(j\) gives the orientation of the surface it is acting on. To find the stress on a surface with normal \(\widehat{\mathbfit{n}}\), take \(\normalsize{𝞼}\cdot\widehat{\mathbfit{n}}\), which gives a vector with components \(\sigma_{ij}n_j\), i.e. \[\int_V \mathbfit{F}_\text{stress} \, {\mathrm d}V = \int_S \mathsfbfit{\normalsize{𝞼}} \cdot \widehat{\mathbfit{n}} \, {\mathrm d}S = \int_S \mathsfbfit{\normalsize{𝞼}} \cdot {\mathrm d}\mathbfit{S}.\]

To specify \(\normalsize{𝞼}\), we need a model for how our particular fluid behaves. If the only internal force we want is pressure, which acts normal to the surfaces, then we can see from the diagram that we just want \(\sigma_{11}\), \(\sigma_{22}\) and \(\sigma_{33}\) to be nonzero. And specifically, we want them to be \(-p(\mathbfit{x},t)\) (negative because pressure points inwards), so in an inviscid fluid, the stress is \[\mathsfbfit{\normalsize{𝞼}} = -p\mathsfbfit{I} \quad \text{ or } \quad \sigma_{ij} = -p\delta_{ij}.\]

Remember that \(p(\mathbfit{x},t)\) is still a function! You may be forgiven for thinking ‘doesn’t \(\sigma_{11}\) on the blue side cancel out the \(-\sigma_{11}\) on the opposite side when you integrate over all the sides?’ The answer is not necessarily, because (in this case) the \(x\)-coordinate is different, and \(p\) could depend on \(x\).

2.1.5 The momentum equation revisited

Putting it all together, then, we have \[\begin{aligned} \int_V \rho\mathchoice{\frac{\mathrm{D} \mathbfit{u}}{\mathrm{D} t}}{\mathrm{D} \mathbfit{u}/\mathrm{D} t}{\mathrm{D} \mathbfit{u}/\mathrm{D} t}{\mathrm{D} \mathbfit{u}/\mathrm{D} t} {\mathrm d}V & = \int_V \mathbfit{F} \, {\mathrm d}V \\ &= \int_S \mathsfbfit{\normalsize{𝞼}} \cdot {\mathrm d}\mathbfit{S} + \int_V \rho \mathbfit{f} \, {\mathrm d}V \\ &= \int_{V}\left[\bm{\nabla}\cdot\normalsize{𝞼}^{\mathsf{T}} + \rho \mathbfit{f} \right]{\mathrm d}V \quad \textrm{(by the divergence theorem)},\\ &= \int_{V}\left[ - \bm{\nabla}p + \rho \mathbfit{f}\right]{\mathrm d}V. \end{aligned}\] In the second and third lines there we have used that \(\bm{\nabla}\cdot \normalsize{𝞼}^{\mathsf{T}}\) is a vector with components \(\displaystyle \bigl(\bm{\nabla}\cdot \normalsize{𝞼}^{\mathsf{T}}\bigr)_i = \nabla_j\sigma^{\mathsf{T}}_{ji} = \mathchoice{\frac{\partial\sigma_{ij}}{\partial x_j}}{\partial\sigma_{ij}/\partial x_j}{\partial\sigma_{ij}/\partial x_j}{\partial\sigma_{ij}/\partial x_j}\).

Of course, we didn’t need to go via the stress tensor to get this result where we just have pressure. [f-pressure] with the divergence theorem gives us \(-\bm{\nabla}p\) directly, but for next term it will be useful to have seen this.

Now, assuming smoothness of the integrand, the only way this result is possible for an arbitrary domain \(V\) is if \(\rho\), \(\mathbfit{u}\) and \(p\) obey the PDE \[\begin{equation} \rho\mathchoice{\frac{\mathrm{D} \mathbfit{u}}{\mathrm{D} t}}{\mathrm{D} \mathbfit{u}/\mathrm{D} t}{\mathrm{D} \mathbfit{u}/\mathrm{D} t}{\mathrm{D} \mathbfit{u}/\mathrm{D} t} = - \bm{\nabla}p + \rho\mathbfit{f} \qquad \iff \qquad \boxed{\mathchoice{\frac{\partial\mathbfit{u}}{\partial t}}{\partial\mathbfit{u}/\partial t}{\partial\mathbfit{u}/\partial t}{\partial\mathbfit{u}/\partial t} + (\mathbfit{u}\cdot\bm{\nabla})\mathbfit{u} = - \frac{1}{\rho}\bm{\nabla}p + \mathbfit{f}}. \label{eqn:mom} \end{equation}\] The form with a general stress is often known as the (Cauchy) momentum equation. The particular form we have here with pressure alone is known as the Euler equations (the plural is a relic from when it was thought of as three separate scalar equations).

2.3.1 Kinetic energy

Just as how kinetic energy of a particle in classical mechanics is \(\frac12 m v^2\), for a fluid occupying a fixed region \(V\), it is \[E = \frac12\int_V\rho_0|\mathbfit{u}|^2 \,{\mathrm d}V.\] Is this conserved? Well, \[\begin{aligned} \mathchoice{\frac{{\mathrm d}E}{{\mathrm d}t}}{{\mathrm d}E/{\mathrm d}t}{{\mathrm d}E/{\mathrm d}t}{{\mathrm d}E/{\mathrm d}t} &=\rho_0 \mathchoice{\frac{{\mathrm d}}{{\mathrm d}t}}{{\mathrm d}/{\mathrm d}t}{{\mathrm d}/{\mathrm d}t}{{\mathrm d}/{\mathrm d}t} \int_V \frac12 |\mathbfit{u}|^2 \, {\mathrm d}V \\ &= \rho_0 \int_V \frac12 \mathchoice{\frac{\partial}{\partial t}}{\partial/\partial t}{\partial/\partial t}{\partial/\partial t} (\mathbfit{u}\cdot\mathbfit{u}) \, {\mathrm d}V \quad \text{(because $V$ is fixed)}\\ &= \rho_0 \int_V \mathbfit{u}\cdot\mathchoice{\frac{\partial\mathbfit{u}}{\partial t}}{\partial\mathbfit{u}/\partial t}{\partial\mathbfit{u}/\partial t}{\partial\mathbfit{u}/\partial t} \, {\mathrm d}V. \end{aligned}\] Substituting in the expression for \(\mathchoice{\frac{\partial\mathbfit{u}}{\partial t}}{\partial\mathbfit{u}/\partial t}{\partial\mathbfit{u}/\partial t}{\partial\mathbfit{u}/\partial t}\) from the Euler equation, [eqn:mom], we have \[\begin{aligned} \mathchoice{\frac{{\mathrm d}E}{{\mathrm d}t}}{{\mathrm d}E/{\mathrm d}t}{{\mathrm d}E/{\mathrm d}t}{{\mathrm d}E/{\mathrm d}t} = \rho_0 \int_V \mathbfit{u}\cdot\left[-(\mathbfit{u}\cdot\bm{\nabla})\mathbfit{u} -\frac{1}{\rho_0}\bm{\nabla}p+\mathbfit{f} \right] {\mathrm d}V. \end{aligned}\] We’re now going to make the questionable decision to substitute in the vector identity \[\begin{equation} (\mathbfit{u}\cdot\bm{\nabla})\mathbfit{u} = (\bm{\nabla}\times\mathbfit{u})\times\mathbfit{u} + \frac12\bm{\nabla}\left(|\mathbfit{u}|^2\right), \label{u-dot-grad-u} \end{equation}\] but notice that when we dot it with \(\mathbfit{u}\), the cross product part disappears. Therefore \[\begin{aligned} \mathchoice{\frac{{\mathrm d}E}{{\mathrm d}t}}{{\mathrm d}E/{\mathrm d}t}{{\mathrm d}E/{\mathrm d}t}{{\mathrm d}E/{\mathrm d}t} = \rho_0 \int_V \mathbfit{u}\cdot\left[-\frac12\bm{\nabla}\left(|\mathbfit{u}|^2\right) -\frac{1}{\rho_0}\bm{\nabla}p+\mathbfit{f} \right] {\mathrm d}V. \end{aligned}\]

We can make progress if we take the body force \(\mathbfit{f}\) to be conservative, i.e. we can write it as \(\mathbfit{f} = -\bm{\nabla}\mathit{\Phi}\) (conventionally we pick the minus sign here). This is a restriction, but it’s an OK restriction because many forces, including gravity, can be written this way. Then \[\mathchoice{\frac{{\mathrm d}E}{{\mathrm d}t}}{{\mathrm d}E/{\mathrm d}t}{{\mathrm d}E/{\mathrm d}t}{{\mathrm d}E/{\mathrm d}t} = -\rho_0 \int_V \mathbfit{u}\cdot \bm{\nabla}\underbrace{\left[ \frac{1}{2}|\mathbfit{u}|^2 + \frac{p}{\rho_0} + \mathit{\Phi} \right]}_{:=H(\mathbfit{x},t)} {\mathrm d}V,\] where we’re introducing the scalar function \(H(\mathbfit{x},t)\).

We can bring the \(\mathbfit{u}\) inside the gradient by using the (rearranged) product rule and the fact that the fluid is incompressible, \[\mathbfit{u}\cdot\bm{\nabla}H = \bm{\nabla}\cdot(H\mathbfit{u}) - H\underbrace{\bm{\nabla}\cdot\mathbfit{u}}_{=0},\] to get \[\begin{equation} \mathchoice{\frac{{\mathrm d}E}{{\mathrm d}t}}{{\mathrm d}E/{\mathrm d}t}{{\mathrm d}E/{\mathrm d}t}{{\mathrm d}E/{\mathrm d}t} = -\rho_0 \int_V \bm{\nabla}\cdot(H\mathbfit{u}) {\mathrm d}V = -\rho_0 \int_S H\mathbfit{u}\cdot \,{\mathrm d}\mathbfit{S}, \label{almost-bern} \end{equation}\] and if we specify that there is no flow through the walls of the volume (i.e. \(\mathbfit{u}\cdot\widehat{\mathbfit{n}}=0\)), then \[\mathchoice{\frac{{\mathrm d}E}{{\mathrm d}t}}{{\mathrm d}E/{\mathrm d}t}{{\mathrm d}E/{\mathrm d}t}{{\mathrm d}E/{\mathrm d}t} = 0.\] Great! So kinetic energy is conserved when an ideal fluid sits in a conservative body force and has no normal flow. Conservation of kinetic energy makes physical sense, since an ideal incompressible fluid has no internal friction or heating to dissipate the energy.

If the flow is steady, we have a stronger result…

2.3.2 Bernoulli’s principle

Given the body force is conservative, we can use our new favourite vector identity to write the momentum equation, [eqn:mom], as \[\begin{align} \mathchoice{\frac{\partial\mathbfit{u}}{\partial t}}{\partial\mathbfit{u}/\partial t}{\partial\mathbfit{u}/\partial t}{\partial\mathbfit{u}/\partial t}+(\bm{\nabla}\times\mathbfit{u})\times\mathbfit{u}&=-\bm{\nabla}\left[\frac{p}{\rho_0}+\frac12 |\mathbfit{u}|^2 + \mathit{\Phi}\right] \label{eqn:b1} \\ &= -\bm{\nabla}H. \end{align}\] If the flow is steady, \(\mathchoice{\frac{\partial\mathbfit{u}}{\partial t}}{\partial\mathbfit{u}/\partial t}{\partial\mathbfit{u}/\partial t}{\partial\mathbfit{u}/\partial t} = 0\) and additionally we can dot both sides with \(\mathbfit{u}\) to get \[0 = \mathbfit{u}\cdot\bm{\nabla}H.\] \(H\) is known as the energy head, and this equation shows contours of \(H\) are parallel to \(\mathbfit{u}\): \(H\) is constant along streamlines.

This is Bernoulli’s principle: in an incompressible steady flow of an ideal fluid, with conservative body force \(\mathbfit{f}=-\bm{\nabla}\mathit{\Phi}\), then \[\begin{aligned} \boxed{H :=\frac{1}{2}|\mathbfit{u}|^2 + \frac{p}{\rho_0} + \mathit{\Phi}} \end{aligned}\] is constant along streamlines.

2.4.1 Circulation dynamics

We saw in 1.7.2 that if we want to quantify the amount of rotation in a region, we don’t use vorticity, which is local, but instead we measure circulation, \[\begin{aligned} \mathit{\Gamma} = \int_C \mathbfit{u} \cdot {\mathrm d}\mathbfit{x}. \end{aligned}\] Now consider a curve \(C(t)\) which evolves with the flow of the underlying fluid. How does the circulation evolve in time?

Firstly, since \(C\) evolves with the flow, this means we want to use the time derivative that follows the flow: \(\mathchoice{\frac{\mathrm{D} }{\mathrm{D} t}}{\mathrm{D} /\mathrm{D} t}{\mathrm{D} /\mathrm{D} t}{\mathrm{D} /\mathrm{D} t}\). This also allows us to bring it under the integral sign. It then acts on \(\mathbfit{u}\cdot{\mathrm d}\mathbfit{x}\) like the normal product rule: \[\begin{align} \mathchoice{\frac{\mathrm{D} \mathit{\Gamma}}{\mathrm{D} t}}{\mathrm{D} \mathit{\Gamma}/\mathrm{D} t}{\mathrm{D} \mathit{\Gamma}/\mathrm{D} t}{\mathrm{D} \mathit{\Gamma}/\mathrm{D} t} &= \oint_{C(t)} \mathchoice{\frac{\mathrm{D} }{\mathrm{D} t}}{\mathrm{D} /\mathrm{D} t}{\mathrm{D} /\mathrm{D} t}{\mathrm{D} /\mathrm{D} t} (\mathbfit{u} \cdot {\mathrm d}\mathbfit{x})\\ &= \oint_{C(t)}\left[ \mathchoice{\frac{\mathrm{D} \mathbfit{u}}{\mathrm{D} t}}{\mathrm{D} \mathbfit{u}/\mathrm{D} t}{\mathrm{D} \mathbfit{u}/\mathrm{D} t}{\mathrm{D} \mathbfit{u}/\mathrm{D} t}\cdot {\mathrm d}\mathbfit{x} + \mathbfit{u}\cdot\mathchoice{\frac{\mathrm{D} ({\mathrm d}\mathbfit{x})}{\mathrm{D} t}}{\mathrm{D} ({\mathrm d}\mathbfit{x})/\mathrm{D} t}{\mathrm{D} ({\mathrm d}\mathbfit{x})/\mathrm{D} t}{\mathrm{D} ({\mathrm d}\mathbfit{x})/\mathrm{D} t}\right].\label{DDtGamma} \end{align}\] We can bring in the Euler equation with a conservative body force, \(\mathbfit{f} = -\bm{\nabla}\mathit{\Phi}\), for the first term, \[\mathchoice{\frac{\mathrm{D} \mathbfit{u}}{\mathrm{D} t}}{\mathrm{D} \mathbfit{u}/\mathrm{D} t}{\mathrm{D} \mathbfit{u}/\mathrm{D} t}{\mathrm{D} \mathbfit{u}/\mathrm{D} t} = -\bm{\nabla}\left[\frac{p}{\rho_0} + \mathit{\Phi}\right].\] The first term in [DDtGamma] is therefore \[\begin{aligned} \oint_{C(t)}\mathchoice{\frac{\mathrm{D} \mathbfit{u}}{\mathrm{D} t}}{\mathrm{D} \mathbfit{u}/\mathrm{D} t}{\mathrm{D} \mathbfit{u}/\mathrm{D} t}{\mathrm{D} \mathbfit{u}/\mathrm{D} t}\cdot {\mathrm d}\mathbfit{x} &= \oint_{C(t)}-\bm{\nabla}\left[\frac{p}{\rho_0} + \mathit{\Phi}\right] \cdot {\mathrm d}\mathbfit{x} \\ & = 0, \end{aligned}\] since this is an integral of a continuous gradient around a closed loop (a property of conservative vector fields).

The first term was zero because of the Euler equations, but the second term, \[\begin{aligned} \oint_{C(t)}\mathbfit{u}\cdot\mathchoice{\frac{\mathrm{D} ({\mathrm d}\mathbfit{x})}{\mathrm{D} t}}{\mathrm{D} ({\mathrm d}\mathbfit{x})/\mathrm{D} t}{\mathrm{D} ({\mathrm d}\mathbfit{x})/\mathrm{D} t}{\mathrm{D} ({\mathrm d}\mathbfit{x})/\mathrm{D} t}, \end{aligned}\] will, in fact, always be zero. Let’s see why.

What is \(\mathchoice{\frac{\mathrm{D} ({\mathrm d}\mathbfit{x})}{\mathrm{D} t}}{\mathrm{D} ({\mathrm d}\mathbfit{x})/\mathrm{D} t}{\mathrm{D} ({\mathrm d}\mathbfit{x})/\mathrm{D} t}{\mathrm{D} ({\mathrm d}\mathbfit{x})/\mathrm{D} t}\)? It’s the rate of change of the length of a line segment, which is the difference in the speeds of the ends: so it’s actually just \({\mathrm d}\mathbfit{u}\)! So \[\begin{aligned} \oint_{C(t)}\mathbfit{u}\cdot\mathchoice{\frac{\mathrm{D} ({\mathrm d}\mathbfit{x})}{\mathrm{D} t}}{\mathrm{D} ({\mathrm d}\mathbfit{x})/\mathrm{D} t}{\mathrm{D} ({\mathrm d}\mathbfit{x})/\mathrm{D} t}{\mathrm{D} ({\mathrm d}\mathbfit{x})/\mathrm{D} t} &= \oint_{C(t)}\mathbfit{u}\cdot{\mathrm d}\mathbfit{u} \\ &= \frac 12 \oint_{C(t)} {\mathrm d}(|\mathbfit{u}|^2) = 0, \end{aligned}\] since this is integrating an exact derivative around a curve.

Putting both sides together, we have Kelvin’s circulation theorem: in an ideal, incompressible fluid with conservative body force, the circulation around a closed material curve \(C(t)\) is constant in time: \[\mathchoice{\frac{\mathrm{D} \mathit{\Gamma}}{\mathrm{D} t}}{\mathrm{D} \mathit{\Gamma}/\mathrm{D} t}{\mathrm{D} \mathit{\Gamma}/\mathrm{D} t}{\mathrm{D} \mathit{\Gamma}/\mathrm{D} t} = 0.\]

You are most likely to have seen the consequences of this watching aeroplanes. Here’s how, but first, flight school 101.

3.3.1 Solution for a specific wavenumber

So we have our simplified, linearised system. Let us first look for separable travelling plane wave solutions of the form \[\begin{equation} \phi(x,z,t) = X(x-ct)Z(z). \label{sep-of-vars} \end{equation}\]

In other words, we look for waves moving at a constant phase speed \(c\) in the \(x\)-direction.

Substituting this into [eqn:wws1] gives \[\begin{align} 0 = {\nabla^2}\phi &= \mathchoice{\frac{\partial^2 \phi}{\partial x^2}}{\partial^2 \phi/\partial x^2}{\partial^2 \phi/\partial x^2}{\partial^2 \phi/\partial x^2} + \mathchoice{\frac{\partial^2 \phi}{\partial z^2}}{\partial^2 \phi/\partial z^2}{\partial^2 \phi/\partial z^2}{\partial^2 \phi/\partial z^2} \\ &= X''Z + X Z'' \qquad \implies - \frac{X''}{X} = \frac{Z''}{Z} = \lambda,\label{vars-sepd} \end{align}\] for some constant \(\lambda\). You will recall this technique – separation of variables – from 1.8. We then have to go through the cases where \(\lambda\) is negative, zero and positive…unless we can decide now which case we want from the boundary conditions.

We can see straight away that when we solve [vars-sepd], one of \(X\) and \(Z\) will be periodic (sin and cos); and the other will be non-periodic (sinh and cosh). Which do we want?

Looking at the original diagram on , we would be happy with a solution which is periodic in \(x-ct\) (in other words, if you freeze time, you’d get a solution which is periodic in \(x\)); and looking at the floor and kinematic boundary conditions, [eqn:wws4] and [eqn:wws3], \[\begin{align} XZ'=0 &\quad\text{on }z=-h, \label{bc2}\\ XZ'=\mathchoice{\frac{\partial\eta}{\partial t}}{\partial\eta/\partial t}{\partial\eta/\partial t}{\partial\eta/\partial t} & \quad\text{on }z=0, \end{align}\] we don’t want a periodic solution in \(Z\). We therefore set \(\lambda=k^2\), and start by solving \[\frac{Z''}{Z} = k^2, \quad Z'(-h) = 0.\] The general solution to the equation is of course \[Z(z) = A\cosh kz + B\sinh kz,\] for constants \(A\) and \(B\), but it could equally be \[Z(z) = A\cosh[k(z+h)] + B\sinh[k(z+h)].\] Choosing this form makes implementing the first boundary condition, \(Z'(-h)=0\), easy: we end up with \(B=0\) and so \[Z(z) = A\cosh\big[k(z+h)\big].\] Meanwhile, the \(X\) equation gives us \[\begin{align} \qquad X(x-ct) &= \widetilde{C}\cos[k(x-ct)] + \widetilde{D}\sin[k(x-ct)] \label{minus-can-be-absorbed} \end{align}\] or equivalently \[\begin{aligned} X(x-ct) &= C \cos[k(x-ct) + \alpha]. \end{aligned}\] Since we can set \(t=0\) to be any reference point we want, we can choose it to be such that the phase, \(\alpha\), is zero, and hence \[\begin{aligned} X(x-ct) &= C \cos[k(x-ct)] = C \operatorname{Re}\bigl\{ \mathrm{e}^{\mathrm{i}k (x-ct)} \bigr\}. \end{aligned}\] Using the exponential form might see like overkill here, but it will prove a bit more convenient in the long run, especially next term when we look at instability. Practically we can treat it like a normal exponential throughout our calculations; we just need to remember to take the real part at the end.

Our solution, for a given \(k\), is therefore \(\phi=XZ\): \[\phi(x,z,t) = \operatorname{Re}\,\Bigl\{ A(k) \underbrace{\cosh\bigl[k(z+h)\bigr]}_{Z(z)}\underbrace{\mathrm{e}^{\mathrm{i}k (x-ct)}}_{X(x-ct)} \Bigr\},\] for a combined, real constant \(A(k)\). What does this mean? Read on…

3.3.2 The dispersion relation

Think about the form of our solution for \(\phi\): this is a periodic wave (\(X\)) with an amplitude (\(Z\)) that decreases exponentially with depth. By taking \(k>0\) without loss of generality (since the sign can be absorbed into the constants in [minus-can-be-absorbed]), we can interpret \(k\) as the wavenumber: the spatial frequency, equal to \(2\pi/\lambda\), where \(\lambda\) is the wavelength. (This interpretation requires \(k>0\) since wavelengths are inherently positive, although of course mathematically, negative \(k\) is fine.) To illustrate the role of \(k\), here’s \(\phi(x,0,0)\) for two different \(k\) values where \(A(k)\equiv1\):

However, we will see that \(c\) and \(k\) are not independent, since we also need to satisfy the kinematic and dynamic boundary conditions on the surface, [eqn:wws2] and [eqn:wws3]. Differentiating the dynamic condition by \(t\) and then combining the two conditions, we have \[\begin{aligned} 0 &= \mathchoice{\frac{\partial^2 \phi}{\partial t^2}}{\partial^2 \phi/\partial t^2}{\partial^2 \phi/\partial t^2}{\partial^2 \phi/\partial t^2} + g\mathchoice{\frac{\partial\eta}{\partial t}}{\partial\eta/\partial t}{\partial\eta/\partial t}{\partial\eta/\partial t} \\ &= \mathchoice{\frac{\partial^2 \phi}{\partial t^2}}{\partial^2 \phi/\partial t^2}{\partial^2 \phi/\partial t^2}{\partial^2 \phi/\partial t^2} + g\mathchoice{\frac{\partial\phi}{\partial z}}{\partial\phi/\partial z}{\partial\phi/\partial z}{\partial\phi/\partial z} \qquad \textrm{on $z=0$}. \end{aligned}\] Substituting our solution \(\phi(x,z,t)\) into this condition gives \[\begin{equation} -c^2k^2\cosh(kh)X + gk\sinh(kh)X = 0, \label{eqn:dispcond} \end{equation}\] and introducing the frequency \(\omega= ck\), we see that \(\omega\) (or \(c\)) and \(k\) must satisfy the dispersion relation, \[\begin{equation} \boxed{\omega^2 = gk\tanh(kh).} \label{eqn:disp} \end{equation}\] This says that the phase speed of wave crests, \(c\), depends on their wavenumber, \(k\). This equivalently says that the speed of a wave depends on its wavelength, \[c^2 = \frac{g\lambda}{2\pi} \tanh\left(\frac{2\pi h}{\lambda}\right).\] Waves of different wavelength travel at different speeds. For this reason, water waves are called dispersive.

We can then find the solution for the free surface \(\eta(x,t)\) by rearranging the dynamic boundary condition, [eqn:wws2], \[\begin{align} \eta(x,t) &= -\frac{1}{g}\left.\mathchoice{\frac{\partial\phi}{\partial t}}{\partial\phi/\partial t}{\partial\phi/\partial t}{\partial\phi/\partial t}\right|_{z=0} \\ & = \operatorname{Re}\left[\frac{\mathrm{i}\omega}{g}\phi(x,0,t)\right] \\ & = -\frac{A(k)\omega}{g}\cosh(kh)\,\operatorname{Im}\big[\mathrm{e}^{\mathrm{i}(kx-\omega t)}\big]\\ & = -\frac{A(k)\omega}{g}\cosh(kh)\sin(kx - \omega t). \label{eqn:eta} \end{align}\] Since the frequency, \(\omega\), is fixed by the dispersion relation, we can see that the solution has two ‘unknowns’: the wavenumber, \(k\), and the amplitude constant, \(A(k)\). Shortly we will use the principle of superposition to say that the full, general solution is the sum of these solutions over all possible values of \(k\), and the constants \(A(k)\) will be set by whatever initial condition is imposed (effectively how much energy is put into the system!). You will have come across this idea when you first studied separation of variables. But for now let’s look at some consequences of the dispersion relation.

What does the dispersion relation, \[\omega^2 = gk\tanh(kh) \iff c^2 = \frac{g\lambda}{2\pi} \tanh\left(\frac{2\pi h}{\lambda}\right),\] mean in practice? The best thing to do here is to plot \(c\) against \(\lambda\). As we do so, there are two important limiting cases. Recall we can take \(k>0\) without loss of generality.

1. Deep water, \(kh\gg 1\). As \(k\to\infty\), the dispersion relation becomes \[\begin{equation} \omega^2 = gk\tanh(kh) \to gk. \label{disp-rel-deep} \end{equation}\] So the wave speed \(c\to\sqrt{g/k}\sim\sqrt{\lambda}\) becomes independent of \(h\). This makes intuitive sense: the waves are so far away from the bottom that they don’t feel the bed. This also says long waves travel faster than short waves.

2. Shallow water/long wavelength, \(kh\to 0\). In this limit, \[\omega^2 = gk\tanh(kh) \to gk^2h,\] so these waves travel at a constant speed \(c=\sqrt{gh}\).

The following plot shows how the dispersion relation [eqn:disp] matches to the deep (dashed blue) and shallow (dotted orange) approximations:

Where are the shallow and deep water approximations useful in the real world?

3.3.3 Group velocity

As we’ve seen, our solution has a wavenumber \(k\), a constant which essentially we can choose. Because our governing equation, \({\nabla^2}\phi=0\), is linear, the general solution, then, is the sum of all possible solutions for all possible \(k\)s (this is the principle of superposition), \[\begin{equation} \phi(x,z,t) = \int_{-\infty}^\infty A(k)\underbrace{\cosh\big[k(z+h)\big]}_{Z(z)}\,\underbrace{\mathrm{e}^{\mathrm{i}(kx-\omega t)}}_{X(x-ct)}\,{\mathrm d}k. \label{eqn:wwgen} \end{equation}\] The integral is here, doing our summation over \(k\) for us, since there is no restriction that \(k\) is an integer.

What might \(\phi\) look like if we plotted it? If you have a specific function in mind, the Fourier integral theorem tells us that the function can be expressed as a sum of sine and cosine waves, each with its own wavenumber and coefficient. This is precisely the approach we’re using to construct our general solution, meaning \(\phi\) could take on virtually any form. The coefficients \(A(k)\), referred to as Fourier coefficients, determine the contribution of each wave to the overall function.

Suppose, then, that we have some localised bump in the fluid, \(\eta\), coming from a potential, \(\phi\):

This bump is a wave packet, formed by summing wave modes with a range of wavenumbers. But the dispersion relation has told us that waves with different wavenumbers travel at different speeds. So, just like at the start of a marathon, as the packet moves forwards, it spreads out:

The question is, at what speed is the whole packet moving?

We can make progress by knowing that in a localised wave packet, \(A(k)\) is typically clustered around one value, \(k_0\). Here’s \(A(k)\) for the example above (a so-called Gaussian wave packet):

Now, since \[\eta(x,t) = -\frac{1}{g}\left.\mathchoice{\frac{\partial\phi}{\partial t}}{\partial\phi/\partial t}{\partial\phi/\partial t}{\partial\phi/\partial t}\right|_{z=0},\] we can see that \(A(k)\) is shared between \(\eta\) and \(\phi\): the clustering around \(k_0\) is therefore also shared.

This gives us licence to Taylor expand \(\omega(k)\) around \(k_0\), \[\omega(k) \approx \omega(k_0) + (k-k_0)c_g, \qquad \text{where} \quad c_g=\left.\frac{{\mathrm d}\omega}{{\mathrm d}k}\right|_{k=k_0},\] and to see what happens when we put this approximation into [eqn:wwgen]: \[\begin{aligned} \phi(x,z,t) &\approx \int_{-\infty}^\infty A(k)\cosh\big[k(z+h)\big]\mathrm{e}^{\mathrm{i}kx - \mathrm{i}\omega(k_0)t - \mathrm{i}(k-k_0)c_g t}\,{\mathrm d}k\\ &= \underbrace{\mathrm{e}^{\mathrm{i}(k_0x - \omega(k_0)t)}}_{\hspace{0cm}\color{C0}(\rm i)}\underbrace{\int_{-\infty}^\infty A(k)\cosh\big[k(z+h)\big]\mathrm{e}^{\mathrm{i}(k-k_0)(x-c_g t)}\,{\mathrm d}k}_{\hspace{0cm}\color{C1}(\rm ii)}. \end{aligned}\] This is:

• A sinusoidal wave with wavenumber \(k_0\) and frequency \(\omega(k_0)\),

• whose envelope travels with speed \(c_g\) in the \(x\)-direction (the speed of the \(k_0\)th mode).

It is the speed of the moving envelope which perhaps best characterises the behaviour of the dispersive wave:

and so \(c_g\) gets the special name of the group velocity. See that the wave packet disperses as the wave modes around \(k_0\) travel faster and slower.

5.1.1 Conservation of energy

Having made this assumption, we now need to check whether this still leads to physical flows. In particular, energy of some form should still be conserved. For incompressible flows, we saw that kinetic energy was conserved, but this time that won’t be the case since kinetic energy can be converted to internal energy (effectively heat) by compressing the same mass into a smaller volume.

So let’s say make exactly that statement: in the absence of a body force, energy, \(E\), is conserved so long as we allow kinetic energy to convert into internal energy, \(e\). It turns out the exact form – we’ll prove it now – is \[\begin{align} \quad E = \int_V\rho\left(\frac{1}{2}|\mathbfit{u}|^2 + e(\rho)\right){\mathrm d}V, \quad \text{where} \quad e(\rho) = \int_0^\rho\frac{P(\sigma)}{\sigma^2}\,{\mathrm d}\sigma \label{E-is-conserved} \end{align}\] (assuming again that no fluid enters our control volume, i.e. \(\mathbfit{u}\cdot\widehat{\mathbfit{n}}=0\) on the boundary). We call \(e(\rho)\) the internal energy density or specific internal energy, meaning internal energy per unit mass. We will see shortly that we indeed experience it as heat.

Let’s now show we have the right form of \(e\) by proving \(E\) is conserved (take a deep breath). Consider each term of [E-is-conserved] separately, starting with the first term (kinetic energy). We’ll start in the same way we did for incompressible fluids in 2.3.1. Since the volume is fixed, we can take the derivative inside the integral to give \[\begin{aligned} \mathchoice{\frac{{\mathrm d}}{{\mathrm d}t}}{{\mathrm d}/{\mathrm d}t}{{\mathrm d}/{\mathrm d}t}{{\mathrm d}/{\mathrm d}t}\int_V\frac12\rho|\mathbfit{u}|^2\;\mathrm{d}V&=\int_V\frac12\mathchoice{\frac{\partial}{\partial t}}{\partial/\partial t}{\partial/\partial t}{\partial/\partial t}\left(\rho|\mathbfit{u}|^2\right){\mathrm d}V\\ &= \int_V\left(\frac12|\mathbfit{u}|^2\mathchoice{\frac{\partial\rho}{\partial t}}{\partial\rho/\partial t}{\partial\rho/\partial t}{\partial\rho/\partial t} + \rho\mathbfit{u}\cdot\mathchoice{\frac{\partial\mathbfit{u}}{\partial t}}{\partial\mathbfit{u}/\partial t}{\partial\mathbfit{u}/\partial t}{\partial\mathbfit{u}/\partial t}\right){\mathrm d}V\\ &= \int_V\left[\frac12|\mathbfit{u}|^2\mathchoice{\frac{\partial\rho}{\partial t}}{\partial\rho/\partial t}{\partial\rho/\partial t}{\partial\rho/\partial t} + \rho\mathbfit{u}\cdot\left(-(\mathbfit{u}\cdot\bm{\nabla})\mathbfit{u} - \frac{1}{\rho}\bm{\nabla}P \right)\right]{\mathrm d}V. \end{aligned}\] We use the same trick as for incompressible fluids and use the goated identity \[\begin{aligned} (\mathbfit{u}\cdot\bm{\nabla})\mathbfit{u} = (\bm{\nabla}\times\mathbfit{u})\times\mathbfit{u} + \frac12\bm{\nabla}(|\mathbfit{u}|^2), \end{aligned}\] along with the observation that the cross product term disappears when dotted with \(\mathbfit{u}\), to get \[\begin{aligned} \mathchoice{\frac{{\mathrm d}}{{\mathrm d}t}}{{\mathrm d}/{\mathrm d}t}{{\mathrm d}/{\mathrm d}t}{{\mathrm d}/{\mathrm d}t}\int_V\frac12\rho|\mathbfit{u}|^2\;\mathrm{d}V&=\int_V\left[\frac12|\mathbfit{u}|^2\mathchoice{\frac{\partial\rho}{\partial t}}{\partial\rho/\partial t}{\partial\rho/\partial t}{\partial\rho/\partial t} - \rho\mathbfit{u}\cdot\bm{\nabla}\left(\frac12|\mathbfit{u}|^2\right) - \mathbfit{u}\cdot\bm{\nabla}P \right]{\mathrm d}V. \end{aligned}\] The middle term on the right-hand side can be manipulated using the rearranged product rule for divergence, \[\begin{align} \mathbfit{A}\cdot\bm{\nabla}f = \bm{\nabla}\cdot(f \mathbfit{A})-f\bm{\nabla}\cdot\mathbfit{A}, \label{prod-rule-for-div} \end{align}\] and then the divergence theorem on the final term in the identity, to get \[\begin{aligned} \mathchoice{\frac{{\mathrm d}}{{\mathrm d}t}}{{\mathrm d}/{\mathrm d}t}{{\mathrm d}/{\mathrm d}t}{{\mathrm d}/{\mathrm d}t}\int_V\frac12\rho|\mathbfit{u}|^2\;\mathrm{d}V&= \int_V\left[\frac12|\mathbfit{u}|^2\mathchoice{\frac{\partial\rho}{\partial t}}{\partial\rho/\partial t}{\partial\rho/\partial t}{\partial\rho/\partial t} + \frac12|\mathbfit{u}|^2\bm{\nabla}\cdot\big(\rho\mathbfit{u}\big) - \mathbfit{u}\cdot\bm{\nabla}P \right]{\mathrm d}V - \int_{S}\frac12|\mathbfit{u}|^2\rho\mathbfit{u}\cdot\,{\mathrm d}\mathbfit{S}. \end{aligned}\] The first two terms on the right form a multiple of the continuity equation, so disappear, and the final term is integrating a multiple of \(\mathbfit{u}\cdot\widehat{\mathbfit{n}}\) around the boundary. But \(\mathbfit{u}\cdot\widehat{\mathbfit{n}}=0\) on the boundary by our original setup of the problem, which leaves us with \[\begin{align} \mathchoice{\frac{{\mathrm d}}{{\mathrm d}t}}{{\mathrm d}/{\mathrm d}t}{{\mathrm d}/{\mathrm d}t}{{\mathrm d}/{\mathrm d}t}\int_V\frac12\rho|\mathbfit{u}|^2\;\mathrm{d}V &= -\int_V\mathbfit{u}\cdot\bm{\nabla}P\;\mathrm{d}V.\label{eqn:prop1} \end{align}\] We can interpret this as saying kinetic energy is converted away when the pressure forces (\(-\bm{\nabla}P\)) act directly against the velocity. For what it’s worth, notice that this equation holds whether or not the fluid is barotropic.

Now consider the internal energy term in [E-is-conserved]. From the chain rule and the definition of \(e\), we have \[\begin{aligned} \mathchoice{\frac{\partial e}{\partial t}}{\partial e/\partial t}{\partial e/\partial t}{\partial e/\partial t} = \mathchoice{\frac{{\mathrm d}e}{{\mathrm d}\rho}}{{\mathrm d}e/{\mathrm d}\rho}{{\mathrm d}e/{\mathrm d}\rho}{{\mathrm d}e/{\mathrm d}\rho}\mathchoice{\frac{\partial\rho}{\partial t}}{\partial\rho/\partial t}{\partial\rho/\partial t}{\partial\rho/\partial t} = \frac{P}{\rho^2}\mathchoice{\frac{\partial\rho}{\partial t}}{\partial\rho/\partial t}{\partial\rho/\partial t}{\partial\rho/\partial t}, \end{aligned}\] so \[\begin{aligned} \mathchoice{\frac{{\mathrm d}}{{\mathrm d}t}}{{\mathrm d}/{\mathrm d}t}{{\mathrm d}/{\mathrm d}t}{{\mathrm d}/{\mathrm d}t}\int_V\rho e\;\mathrm{d}V&= \int_V\left(e\mathchoice{\frac{\partial\rho}{\partial t}}{\partial\rho/\partial t}{\partial\rho/\partial t}{\partial\rho/\partial t} + \rho\mathchoice{\frac{\partial e}{\partial t}}{\partial e/\partial t}{\partial e/\partial t}{\partial e/\partial t} \right)\;\mathrm{d}V\\ &= \int_V\left(e + \frac{P}{\rho} \right)\mathchoice{\frac{\partial\rho}{\partial t}}{\partial\rho/\partial t}{\partial\rho/\partial t}{\partial\rho/\partial t}\;\mathrm{d}V\\ &= -\int_V\left(e + \frac{P}{\rho} \right)\bm{\nabla}\cdot\big(\rho\mathbfit{u}\big)\;\mathrm{d}V\\ &= \int_V\rho\mathbfit{u}\cdot\bm{\nabla}\underbrace{\left(e + \frac{P}{\rho} \right)}_{\text{function of $\rho$}}\;\mathrm{d}V- \int_{S}\left(e + \frac{P}{\rho} \right)\underbrace{\rho\mathbfit{u}\cdot\,{\mathrm d}\mathbfit{S}}_{\mathbfit{u}\cdot\widehat{\mathbfit{n}}=0 \text{ on $S$}}\\ &= \int_V\rho\mathbfit{u}\cdot\left(\mathchoice{\frac{{\mathrm d}e}{{\mathrm d}\rho}}{{\mathrm d}e/{\mathrm d}\rho}{{\mathrm d}e/{\mathrm d}\rho}{{\mathrm d}e/{\mathrm d}\rho} + \frac{P'(\rho)}{\rho} - \frac{P}{\rho^2} \right)\bm{\nabla}\rho\;\mathrm{d}V\\ &= \int_V\rho\mathbfit{u}\cdot\left(\frac{P}{\rho^2} + \frac{P'(\rho)}{\rho} - \frac{P}{\rho^2} \right)\bm{\nabla}\rho\;\mathrm{d}V\\ &= \int_V\rho\mathbfit{u}\cdot\frac{P'(\rho)}{\rho}\bm{\nabla}\rho \;\mathrm{d}V\\ &= \int_V\mathbfit{u}\cdot\bm{\nabla}P\;\mathrm{d}V. \end{aligned}\] This is equal and opposite to [eqn:prop1], hence the total energy \(E\) is conserved! Phew…

5.1.2 A new Bernoulli

As with incompressible fluids, we can also find a conserved energy-related quantity along streamlines, generalising Bernoulli’s principle.

Recall that we derived Bernoulli’s principle in 2.3.2 by writing down the momentum equation for a conservative body force, \(\mathbfit{f}=-\bm{\nabla}\mathit{\Phi}\), \[\mathchoice{\frac{\partial\mathbfit{u}}{\partial t}}{\partial\mathbfit{u}/\partial t}{\partial\mathbfit{u}/\partial t}{\partial\mathbfit{u}/\partial t} + (\mathbfit{u}\cdot\bm{\nabla})\mathbfit{u} = - \frac{1}{\rho} \bm{\nabla}p - \bm{\nabla}\mathit{\Phi},\] and then using our favourite vector identity on the second term to turn the whole right-hand side into \(\bm{\nabla}[\text{stuff}]\). Dotting the equation with \(\mathbfit{u}\) and assuming steady flow gave us the result. But this time, it’s a bit harder to write the pressure term in the form \(\bm{\nabla}[\text{stuff}]\) because \(\rho\) isn’t a constant.

That is, unless we just decide we can and see what happens. Let’s introduce \(h(\rho)\) so that the momentum equation becomes \[\mathchoice{\frac{\partial\mathbfit{u}}{\partial t}}{\partial\mathbfit{u}/\partial t}{\partial\mathbfit{u}/\partial t}{\partial\mathbfit{u}/\partial t} + (\mathbfit{u}\cdot\bm{\nabla})\mathbfit{u} = - \bm{\nabla}h - \bm{\nabla}\mathit{\Phi}.\] What does \(h\) have to be to make this work? Well, \[\begin{aligned} \bm{\nabla}h(\rho) &= \frac{1}{\rho}\bm{\nabla}P(\rho) \\ \mathchoice{\frac{{\mathrm d}h}{{\mathrm d}\rho}}{{\mathrm d}h/{\mathrm d}\rho}{{\mathrm d}h/{\mathrm d}\rho}{{\mathrm d}h/{\mathrm d}\rho}\bm{\nabla}\rho &= \frac{1}{\rho}\mathchoice{\frac{{\mathrm d}P}{{\mathrm d}\rho}}{{\mathrm d}P/{\mathrm d}\rho}{{\mathrm d}P/{\mathrm d}\rho}{{\mathrm d}P/{\mathrm d}\rho}\bm{\nabla}\rho \\ \mathchoice{\frac{{\mathrm d}h}{{\mathrm d}\rho}}{{\mathrm d}h/{\mathrm d}\rho}{{\mathrm d}h/{\mathrm d}\rho}{{\mathrm d}h/{\mathrm d}\rho} &= \frac{1}{\rho}\mathchoice{\frac{{\mathrm d}P}{{\mathrm d}\rho}}{{\mathrm d}P/{\mathrm d}\rho}{{\mathrm d}P/{\mathrm d}\rho}{{\mathrm d}P/{\mathrm d}\rho} \implies h(\rho) = \int_0^\rho\frac{P'(\sigma)}{\sigma}\,{\mathrm d}\sigma. \end{aligned}\] And the method in 2.3.2 now works to show that the energy head \(H\) generalises to \[H = h(\rho) + \frac{1}{2}|\mathbfit{u}|^2+\mathit{\Phi}.\]

5.1.3 Some thermodynamics

A very common model for pressure is \[P(\rho) = k\rho^\gamma,\] which represents an ideal gas undergoing an adiabatic process (no additional heat added). The constant \(\gamma\) is called the adiabatic index and depends on the number of degrees of freedom of the molecules making up the gas. Specifically, \(\gamma=(n_{\rm dof} + 2)/n_{\rm dof}\). Diatomic molecules such as \({\rm N}_2\) or \({\rm O}_2\) have \(n_{\rm dof}=5\) (3 translational + 2 rotational), so \(\gamma=\tfrac75=1.4\). Thus air has \(\gamma\approx 1.4\).

Meanwhile, the internal energy density, \(e\), is related to temperature by \[e = \frac{c_VT}{\mu},\] where \(c_V\) is the molar heat capacity at constant volume (how much energy you have to add to one mole of a substance to increase the temperature by one unit), and \(\mu\) is the mass per mole of a substance. For an ideal gas, \(c_V=\tfrac12n_\text{dof}R\), where \(R\approx \qty{8.3}{\J\per\K\per\mole}\) is the universal gas constant.

If you combine this relation with the definition of \(e\) in [E-is-conserved], you might be able to see this is a generalisation of ‘\(pV=nRT\)’ from A-level physics.

5.2.1 Small-amplitude sound waves

We start from the unforced (\(\mathbfit{f}=\mathbf{0}\)) compressible barotropic Euler equations: \[\begin{align} \mathchoice{\frac{\partial\rho}{\partial t}}{\partial\rho/\partial t}{\partial\rho/\partial t}{\partial\rho/\partial t} + \bm{\nabla}\cdot(\rho\mathbfit{u}) &= 0,\label{eqn:ceul1}\\ \mathchoice{\frac{\partial\mathbfit{u}}{\partial t}}{\partial\mathbfit{u}/\partial t}{\partial\mathbfit{u}/\partial t}{\partial\mathbfit{u}/\partial t} + (\mathbfit{u}\cdot\bm{\nabla})\mathbfit{u} &= -\frac{1}{\rho}\bm{\nabla}p,\label{eqn:ceul2}\\ p &= P(\rho). \end{align}\] We will consider disturbances around a base state with uniform density \(\rho_0\), uniform pressure \(p_0=P(\rho_0)\), and zero velocity. Replacing \(\rho\to\rho_0+\varepsilon\rho_1\), \(p\to p_0+\varepsilon p_1\) and \(\mathbfit{u} = \varepsilon\mathbfit{u}_1\) in [eqn:ceul1] and [eqn:ceul2] gives \[\begin{align} \varepsilon\mathchoice{\frac{\partial\rho_1}{\partial t}}{\partial\rho_1/\partial t}{\partial\rho_1/\partial t}{\partial\rho_1/\partial t} + \varepsilon\rho_0\bm{\nabla}\cdot\mathbfit{u}_1 + \varepsilon^2\bm{\nabla}\cdot(\rho_1\mathbfit{u}_1) &= 0,\label{eqn:ceul11}\\ \varepsilon\mathchoice{\frac{\partial\mathbfit{u}_1}{\partial t}}{\partial\mathbfit{u}_1/\partial t}{\partial\mathbfit{u}_1/\partial t}{\partial\mathbfit{u}_1/\partial t} + \varepsilon^2(\mathbfit{u}_1\cdot\bm{\nabla})\mathbfit{u}_1 &= -\frac{\varepsilon}{\rho_0+\varepsilon\rho_1}\bm{\nabla}p_1,\label{eqn:ceul12}\\ p_0+\varepsilon p_1 &= P(\rho_0+\varepsilon\rho_1).\label{eqn:ceul13} \end{align}\]

To linearise the awkward terms on the right-hand sides, we use Taylor expansion to write \[\frac{1}{\rho_0+\varepsilon\rho_1} = \frac{1}{\rho_0}\left( \frac{1}{1+\varepsilon\frac{\rho_1}{\rho_0}} \right) = \frac{1}{\rho_0}\left(1 - \varepsilon\frac{\rho_1}{\rho_0} + \varepsilon^2\frac{\rho_1^2}{\rho_0^2} + \cdots \right)\] and also \[\begin{aligned} P(\rho_0+\varepsilon\rho_1) &= P(\rho_0) + \varepsilon\rho_1 P'(\rho_0) + \frac{\varepsilon^2\rho_1^2}{2!}P''(\rho_0) + \cdots\\ &= p_0 + \varepsilon\rho_1 P'(\rho_0) + \frac{\varepsilon^2\rho_1^2}{2!}P''(\rho_0) + \cdots, \end{aligned}\] so that the linearised equations (at order \(\varepsilon\) and dropping the subscript \(1\)) are \[\begin{align} \mathchoice{\frac{\partial\rho}{\partial t}}{\partial\rho/\partial t}{\partial\rho/\partial t}{\partial\rho/\partial t} + \rho_0\bm{\nabla}\cdot\mathbfit{u} &= 0,\label{eqn:lin1}\\ \mathchoice{\frac{\partial\mathbfit{u}}{\partial t}}{\partial\mathbfit{u}/\partial t}{\partial\mathbfit{u}/\partial t}{\partial\mathbfit{u}/\partial t} &= -\frac{1}{\rho_0}\bm{\nabla}p,\label{eqn:lin2}\\ p &= P'(\rho_0)\rho.\label{eqn:lin3} \end{align}\]

5.2.2 1D travelling wave solutions

We’ll start with 1D and build up. If we assume that \(\phi=\phi(x,t)\), then [eqn:soundPhi] reduces to \[\begin{equation} \mathchoice{\frac{\partial^2 \phi}{\partial t^2}}{\partial^2 \phi/\partial t^2}{\partial^2 \phi/\partial t^2}{\partial^2 \phi/\partial t^2} = c_0^2\mathchoice{\frac{\partial^2 \phi}{\partial x^2}}{\partial^2 \phi/\partial x^2}{\partial^2 \phi/\partial x^2}{\partial^2 \phi/\partial x^2}. \label{eqn:sound1d} \end{equation}\]

We’ve already seen last term that the wave equation admits travelling wave solutions of the form \[\begin{equation} \phi(x,t) = F(x-c_0t) + G(x+c_0t) \label{eqn:FG} \end{equation}\] and when given \(\phi(x,0)\) their particular forms can be found from a Fourier expansion. However, here we have two constraints on \(\phi\), given by specifying pressure, \(p\), as well as \(\mathbfit{u}\) at \(t=0\). Thus, we are trying to solve the initial value problem given by \[\begin{align} \mathchoice{\frac{\partial^2 \phi}{\partial t^2}}{\partial^2 \phi/\partial t^2}{\partial^2 \phi/\partial t^2}{\partial^2 \phi/\partial t^2} &= c_0^2\mathchoice{\frac{\partial^2 \phi}{\partial x^2}}{\partial^2 \phi/\partial x^2}{\partial^2 \phi/\partial x^2}{\partial^2 \phi/\partial x^2},\\ \phi(x,0)&=\phi_0(x), & \quad \text{(a condition on $\mathbfit{u}$)} \label{sound-1d-c1} \\ \mathchoice{\frac{\partial\phi}{\partial t}}{\partial\phi/\partial t}{\partial\phi/\partial t}{\partial\phi/\partial t}(x,0) &= \dot\phi_0(x). & \quad \text{(a condition on $p$)} \label{sound-1d-c2} \end{align}\] Let’s see what the solution looks like.

We first compute the two initial conditions. The first, [sound-1d-c1], is \[\begin{align} \phi_0(x) &= F(x) + G(x),\label{eqn:dal1} \end{align}\] The second condition, [sound-1d-c2], is \[\begin{align} \dot\phi_0(x) &= -c_0F'(x) + c_0G'(x), \label{eqn:dal2}\\ \implies \int_{a}^x\dot\phi_0(s)\,{\mathrm d}s &= -c_0F(x) + c_0G(x),\label{eqn:dal3} \end{align}\] where \(a\) is an arbitrary point. Solving [eqn:dal1] and [eqn:dal3] simultaneously gives \[\begin{align} F(x) = \frac{1}{2}\phi_0(x) - \frac{1}{2c_0}\int_{a}^x\dot\phi_0(s)\,{\mathrm d}s,\label{eqn:dalF}\\ G(x) = \frac{1}{2}\phi_0(x) + \frac{1}{2c_0}\int_{a}^x\dot\phi_0(s)\,{\mathrm d}s.\label{eqn:dalG} \end{align}\] This gives us the functional form for \(F\) and \(G\). Plugging this back into [eqn:FG] and assuming that \(a\) is between the diverging points defined by \(x-c_0t\) and \(x+c_0t\) we get \[\begin{aligned} \phi(x,t) &= \frac{1}{2}\phi_0(x-c_0t) - \frac{1}{2c_0}\int_{a}^{x-c_0t}\dot\phi_0(s)\,{\mathrm d}s + \frac{1}{2}\phi_0(x+c_0t) + \frac{1}{2c_0}\int_{a}^{x+c_0t}\dot\phi_0(s)\,{\mathrm d}s \\ &= \frac{1}{2}\Big[\phi_0(x-c_0t) + \phi_0(x+c_0t)\Big] + \frac{1}{2c_0}\int_{x-c_0t}^{x+c_0t}\dot\phi_0(s)\,{\mathrm d}s. \end{aligned}\] This elegant solution is known as D’Alembert’s solution. Let’s try it out!

How well does our linear approximation do when the full equations are solved?

5.2.3 1D standing wave solutions

Now we have an idea of the general properties of the waves, let start introducing more effects, starting with walls. Let’s assume the domain is now \(\{x : 0<x<L\}\). As for water waves, we can then look for standing wave solutions of the form \[\phi(x,t) = X(x)\sin(\omega t).\] Substituting this into our linear wave equation for \(\phi\) gives the ODE \[\begin{equation} -\omega^2X = c_0^2X'', \label{eqn:soundode} \end{equation}\] which of course has solution \[\begin{equation} X(x) = A \sin kx + B \cos kx, \quad k = \omega/c_0. \label{eqn:soundode-soln} \end{equation}\] The solution \(\phi\) must satisfy suitable boundary conditions at \(x=0\) and \(x=L\), which depend on the specific problem. These will either be in the form of a constraint on pressure or velocity.

5.2.4 3D travelling wave solutions

Let’s now generalise these concepts to 3D. In 1D, we saw that the wave equation, [eqn:soundPhi], admits a solution which is a superposition of a left-moving wave and a right-moving wave. In 3D, the solution will be a superposition of plane waves travelling in many different directions. Each component will have the form \[\phi(\mathbfit{x},t) = F(\mathbfit{k}\cdot\mathbfit{x} - \omega t)\] for some fixed vector \(\mathbfit{k}\) whose direction represents the direction of motion. Differentiating this expression, we have \[\begin{aligned} &\mathchoice{\frac{\partial^2 \phi}{\partial t^2}}{\partial^2 \phi/\partial t^2}{\partial^2 \phi/\partial t^2}{\partial^2 \phi/\partial t^2} = \omega^2 F''(\mathbfit{k}\cdot\mathbfit{x} - \omega t),\\ &{\nabla^2}\phi = \mathchoice{\frac{\partial^2 \phi}{\partial x^2}}{\partial^2 \phi/\partial x^2}{\partial^2 \phi/\partial x^2}{\partial^2 \phi/\partial x^2} + \mathchoice{\frac{\partial^2 \phi}{\partial y^2}}{\partial^2 \phi/\partial y^2}{\partial^2 \phi/\partial y^2}{\partial^2 \phi/\partial y^2} + \mathchoice{\frac{\partial^2 \phi}{\partial z^2}}{\partial^2 \phi/\partial z^2}{\partial^2 \phi/\partial z^2}{\partial^2 \phi/\partial z^2} = (k_x^2 + k_y^2 + k_z^2)F''(\mathbfit{k}\cdot\mathbfit{x} - \omega t), \end{aligned}\] so this will be a solution to [eqn:soundPhi] providing that \(\omega\) satisfies the dispersion relation \[\begin{equation} \omega^2 = c_0^2|\mathbfit{k}|^2. \label{eqn:soundis3} \end{equation}\]

The corresponding velocity perturbation is \[\mathbfit{u}(\mathbfit{x},t) = \bm{\nabla}\phi = \mathbfit{k} F'(\mathbfit{k}\cdot\mathbfit{x} - \omega t),\] so we see that the oscillations are in the direction of propagation, \(\mathbfit{k}\). This confirms our intuition from 1D that sound waves are longitudinal.

5.2.5 3D waveguide solutions

We saw in 1D that a domain of finite length supports a discrete spectrum of standing wave modes. Similarly, a domain of finite width, called a waveguide, supports a discrete spectrum of propagating modes. Think of the sounds of a didgeridoo or someone shouting down a corridor.

We can generalise this approach to 3D wave guides such as tubes with a square or circular cross sections. The former is a problem sheet question. Let’s try the latter as it’ll give us more practice with cylindrical coordinates!

5.3.1 The method of characteristics

In short, the method of characteristics says that solving PDEs is hard, so instead we should try to find curves in the parameter space – in this case, \((x,t)\) – along which the problem is easier to solve. Maybe along that curve the problem reduces to a nice ODE, or, even better, the solution is constant. Either way, to find the solution for given parameters, you just have to find what curve it’s on, and follow that curve along until you hit some known piece of information (an initial or boundary condition). Such curves are called characteristics.

The good news is that you’ve actually seen an example of this already! Recall that the general solution to the wave equation \[\mathchoice{\frac{\partial^2 y}{\partial t^2}}{\partial^2 y/\partial t^2}{\partial^2 y/\partial t^2}{\partial^2 y/\partial t^2} = c^2\mathchoice{\frac{\partial^2 y}{\partial x^2}}{\partial^2 y/\partial x^2}{\partial^2 y/\partial x^2}{\partial^2 y/\partial x^2},\] is \(y = F(x-ct) + G(x+ct)\). We could also have arrived at this result by making a change of variables of \(u=x-ct\) and \(v=x+ct\). Applying the chain rule a couple of times (check that you agree) we find \[\mathchoice{\frac{\partial^2 y}{\partial u \partial v}}{\partial^2 y /\partial u \partial v}{\partial^2 y /\partial u \partial v}{\partial^2 y /\partial u \partial v}=0 \implies y = F(u) + G(v) \quad \text{after integrating.}\] How is this any different? Well \(u=\text{constant}\) is a curve in the \((x,t)\)-plane along which \(F\) is constant. If we set \(x=x_0\) at \(t=0\), then \(u = x-ct=x_0\) is our curve. Similarly, \(v=\text{constant}\) is a curve along which \(G\) is constant, and we can use our initial condition to find our curve to be \(v=x+ct = x_0\). These curves are characteristics of the wave equation and, in this case, they represent the paths along which the waves propagate:

We can see that our right and left moving solutions are represented by red and blue lines in the plot. Every \(x_0\) will have its own lines, so we have curves passing through every value of \(x\) and \(t\). Importantly, as \(c\) is constant, each \(x\) and \(t\) has a unique pair of starting points on the line \(t=0\) which can be found by following the characteristics in reverse:

This means the wave equation is well-behaved and solutions are unique. As we’ll see, this isn’t true for the Euler equations!

5.3.2 Burgers’ equation and blow-up

Returning to the Euler equations, let’s start by looking at the inertial terms in the momentum equation (the left-hand side) in 1D: \[\begin{equation} \mathchoice{\frac{\partial u}{\partial t}}{\partial u/\partial t}{\partial u/\partial t}{\partial u/\partial t} + u\mathchoice{\frac{\partial u}{\partial x}}{\partial u/\partial x}{\partial u/\partial x}{\partial u/\partial x} = 0. \label{eqn:burgers} \end{equation}\] This 1D equation has its own name: the inviscid Burgers’ equation.

Suppose we have a curve \(x(t)\) in the \((x,t)\)-plane where \[\begin{equation} \mathchoice{\frac{{\mathrm d}x}{{\mathrm d}t}}{{\mathrm d}x/{\mathrm d}t}{{\mathrm d}x/{\mathrm d}t}{{\mathrm d}x/{\mathrm d}t} = u\big(x(t), t\big), \qquad x(0)=x_0. \label{eqn:char} \end{equation}\] Then along such a curve [eqn:burgers] may be written as \[\begin{equation} \mathchoice{\frac{\partial u}{\partial t}}{\partial u/\partial t}{\partial u/\partial t}{\partial u/\partial t} + \mathchoice{\frac{{\mathrm d}x}{{\mathrm d}t}}{{\mathrm d}x/{\mathrm d}t}{{\mathrm d}x/{\mathrm d}t}{{\mathrm d}x/{\mathrm d}t}\mathchoice{\frac{\partial u}{\partial x}}{\partial u/\partial x}{\partial u/\partial x}{\partial u/\partial x} = 0 \qquad \iff \qquad \mathchoice{\frac{{\mathrm d}u}{{\mathrm d}t}}{{\mathrm d}u/{\mathrm d}t}{{\mathrm d}u/{\mathrm d}t}{{\mathrm d}u/{\mathrm d}t} = 0. \label{eqn:charu} \end{equation}\] Along these curves, the PDE reduces to an ODE: these curves are therefore characteristics of Burgers’ equation! This time, what they represent differs from our wave equation example: [eqn:char] shows that these curves represent particle paths.

So far so good. But what about if the particles behind were moving faster? Then surely they would catch up?

5.3.3 Riemann invariants

Now let us put back in the other terms in the 1D Euler equations. For simplicity, we will consider an ideal gas, \(p=P(\rho)=k\rho^\gamma\). So the continuity and unforced momentum equations are \[\begin{align} \mathchoice{\frac{\partial\rho}{\partial t}}{\partial\rho/\partial t}{\partial\rho/\partial t}{\partial\rho/\partial t}+\bm{\nabla}\cdot(\rho\mathbfit{u})&=0 &\implies&& \mathchoice{\frac{\partial\rho}{\partial t}}{\partial\rho/\partial t}{\partial\rho/\partial t}{\partial\rho/\partial t} + \rho\mathchoice{\frac{\partial u}{\partial x}}{\partial u/\partial x}{\partial u/\partial x}{\partial u/\partial x} + u\mathchoice{\frac{\partial\rho}{\partial x}}{\partial\rho/\partial x}{\partial\rho/\partial x}{\partial\rho/\partial x} &= 0,\label{eqn:eul1a}\\ \mathchoice{\frac{\partial\mathbfit{u}}{\partial t}}{\partial\mathbfit{u}/\partial t}{\partial\mathbfit{u}/\partial t}{\partial\mathbfit{u}/\partial t}+(\mathbfit{u}\cdot\bm{\nabla})\mathbfit{u}&=-\frac{1}{\rho}\bm{\nabla}p &\implies&& \mathchoice{\frac{\partial u}{\partial t}}{\partial u/\partial t}{\partial u/\partial t}{\partial u/\partial t} + u\mathchoice{\frac{\partial u}{\partial x}}{\partial u/\partial x}{\partial u/\partial x}{\partial u/\partial x} + \frac{k\gamma\rho^{\gamma-1}}{\rho}\mathchoice{\frac{\partial\rho}{\partial x}}{\partial\rho/\partial x}{\partial\rho/\partial x}{\partial\rho/\partial x} &= 0\label{eqn:eul1b}. \end{align}\] We will see that these may be reduced to two sets of ODEs if we first eliminate \(\rho\) in favour of the variable \[c(\rho) = \sqrt{P'(\rho)} = \sqrt{k\gamma\rho^{\gamma-1}}.\] We will see later that it still makes sense to call this the sound speed, although it now varies in space. To eliminate \(\rho\) from [eqn:eul1a] and [eqn:eul1b], we use the fact that \[\begin{aligned} \mathchoice{\frac{\partial c}{\partial t}}{\partial c/\partial t}{\partial c/\partial t}{\partial c/\partial t} & = c'(\rho)\mathchoice{\frac{\partial\rho}{\partial t}}{\partial\rho/\partial t}{\partial\rho/\partial t}{\partial\rho/\partial t} \\ & = \sqrt{k\gamma}\frac{\gamma-1}{2}\rho^{(\gamma-3)/2}\mathchoice{\frac{\partial\rho}{\partial t}}{\partial\rho/\partial t}{\partial\rho/\partial t}{\partial\rho/\partial t} \\ & = \frac{\gamma -1}{2}\frac{c}{\rho}\mathchoice{\frac{\partial\rho}{\partial t}}{\partial\rho/\partial t}{\partial\rho/\partial t}{\partial\rho/\partial t} \\ \implies \mathchoice{\frac{\partial\rho}{\partial t}}{\partial\rho/\partial t}{\partial\rho/\partial t}{\partial\rho/\partial t} &= \frac{2\rho}{(\gamma-1)c}\mathchoice{\frac{\partial c}{\partial t}}{\partial c/\partial t}{\partial c/\partial t}{\partial c/\partial t}. \end{aligned}\] In the same way, \[\mathchoice{\frac{\partial c}{\partial x}}{\partial c/\partial x}{\partial c/\partial x}{\partial c/\partial x} = c'(\rho)\mathchoice{\frac{\partial\rho}{\partial x}}{\partial\rho/\partial x}{\partial\rho/\partial x}{\partial\rho/\partial x} \quad \implies \quad \mathchoice{\frac{\partial\rho}{\partial x}}{\partial\rho/\partial x}{\partial\rho/\partial x}{\partial\rho/\partial x} = \frac{2\rho}{(\gamma-1)c}\mathchoice{\frac{\partial c}{\partial x}}{\partial c/\partial x}{\partial c/\partial x}{\partial c/\partial x}.\] Substituting these into [eqn:eul1a] and [eqn:eul1b] (and multiplying the first equation by \(c/\rho\)) we find \[\begin{aligned} \frac{2}{\gamma-1}\mathchoice{\frac{\partial c}{\partial t}}{\partial c/\partial t}{\partial c/\partial t}{\partial c/\partial t} + c\mathchoice{\frac{\partial u}{\partial x}}{\partial u/\partial x}{\partial u/\partial x}{\partial u/\partial x} + \frac{2u}{\gamma-1}\mathchoice{\frac{\partial c}{\partial x}}{\partial c/\partial x}{\partial c/\partial x}{\partial c/\partial x} &= 0,\\ \mathchoice{\frac{\partial u}{\partial t}}{\partial u/\partial t}{\partial u/\partial t}{\partial u/\partial t} + u\mathchoice{\frac{\partial u}{\partial x}}{\partial u/\partial x}{\partial u/\partial x}{\partial u/\partial x} + \frac{2c}{\gamma-1}\mathchoice{\frac{\partial c}{\partial x}}{\partial c/\partial x}{\partial c/\partial x}{\partial c/\partial x} &= 0. \end{aligned}\] The trick is then to add or subtract these two equations, to find \[\begin{align} &\left[\mathchoice{\frac{\partial}{\partial t}}{\partial/\partial t}{\partial/\partial t}{\partial/\partial t} + (u + c)\mathchoice{\frac{\partial}{\partial x}}{\partial/\partial x}{\partial/\partial x}{\partial/\partial x} \right]\left[u + \frac{2c}{\gamma-1}\right] = 0,\label{eqn:rie1}\\ &\left[\mathchoice{\frac{\partial}{\partial t}}{\partial/\partial t}{\partial/\partial t}{\partial/\partial t} + (u - c)\mathchoice{\frac{\partial}{\partial x}}{\partial/\partial x}{\partial/\partial x}{\partial/\partial x} \right]\left[u - \frac{2c}{\gamma-1}\right] = 0.\label{eqn:rie2} \end{align}\] Thus, by exactly the same logic as at the start of 2.3.2, \[\begin{align} &F_+ = u + \frac{2c}{\gamma-1} \text{ is constant on curves $x_+(t)$ satisfying } \mathchoice{\frac{{\mathrm d}x_+}{{\mathrm d}t}}{{\mathrm d}x_+/{\mathrm d}t}{{\mathrm d}x_+/{\mathrm d}t}{{\mathrm d}x_+/{\mathrm d}t} = u+c,\label{eqn:rie3}\\ &F_- = u - \frac{2c}{\gamma-1} \text{ is constant on curves $x_-(t)$ satisfying } \mathchoice{\frac{{\mathrm d}x_-}{{\mathrm d}t}}{{\mathrm d}x_-/{\mathrm d}t}{{\mathrm d}x_-/{\mathrm d}t}{{\mathrm d}x_-/{\mathrm d}t} = u-c.\label{eqn:rie4} \end{align}\] The functions \(F_+\) and \(F_-\) are called Riemann invariants. If you know \(F_+\) and \(F_-\), then you know \(u\) and \(c\). And if you know \(c\), you know \(\rho\). And if you know \(\rho\), you know \(p\). So these two invariants tell us everything we normally ask for when solving a fluids problem: velocity and pressure.

The two curves are known as the \(+\) and \(-\) characteristics. They need no longer be straight lines since \(u\) and \(c\) need no longer be constant along them. Maybe something like this:

5.3.4 Simple waves

In general it is hard to find closed-form solutions for the characteristics. However, it is possible in some cases, called simple waves, in which one of the Riemann invariants is constant across a whole region.

6.1.1 The stress tensor is symmetric

Conservation of angular momentum for a fluid element means that \[\begin{align} \underbrace{\mathchoice{\frac{{\mathrm d}}{{\mathrm d}t}}{{\mathrm d}/{\mathrm d}t}{{\mathrm d}/{\mathrm d}t}{{\mathrm d}/{\mathrm d}t}\int_{V(t)}\mathbfit{x}\times\rho\mathbfit{u}\;\mathrm{d}V}_{\text{change in ang. mom.}} = \underbrace{\int_{S(t)}\mathbfit{x}\times(\normalsize{𝞼}\cdot\,{\mathrm d}\mathbfit{S})}_{\text{torque from boundary stress}} + \underbrace{\int_{V(t)}\mathbfit{x}\times\rho\mathbfit{f}\;\mathrm{d}V}_{\text{torque from body forces}},\label{eqn:angm} \end{align}\] where \(\mathbfit{x}\) is the position vector from the origin and \(V(t), S(t)\) are the volume and surface of the fluid element at a given time.

As we know, conservation of momentum leads to statements like ‘\(\mathbfit{F}=m\mathbfit{a}\)’, and so conservation of angular momentum leads to statements like ‘\(\mathbfit{x}\times\mathbfit{F}=\mathbfit{x}\times m\mathbfit{a}\)’: the form of this equation should not be surprising.

There is, however, something new here: our volume and surface now depend on time because this is a statement about a moving fluid element. So far, all of our integral equations have involved fixed volumes, where the fluid can pass through the volume. Not this time! The challenge that arises is how to bring the derivative inside the integral in the first term. For a fixed \(V\), \(\mathchoice{\frac{{\mathrm d}}{{\mathrm d}t}}{{\mathrm d}/{\mathrm d}t}{{\mathrm d}/{\mathrm d}t}{{\mathrm d}/{\mathrm d}t} \to \mathchoice{\frac{\partial}{\partial t}}{\partial/\partial t}{\partial/\partial t}{\partial/\partial t}\) as it moves inside; but if \(V\) moves, we have to account for that too.

The result we need is known as the Reynolds transport theorem, and no fluids course is complete without seeing it! It says, for a function \(f(\mathbfit{x},t)\), \[\mathchoice{\frac{{\mathrm d}}{{\mathrm d}t}}{{\mathrm d}/{\mathrm d}t}{{\mathrm d}/{\mathrm d}t}{{\mathrm d}/{\mathrm d}t} \int_{V(t)} f \rho \, {\mathrm d}V = \int_{V(t)} \mathchoice{\frac{\mathrm{D} f}{\mathrm{D} t}}{\mathrm{D} f/\mathrm{D} t}{\mathrm{D} f/\mathrm{D} t}{\mathrm{D} f/\mathrm{D} t} \rho \, {\mathrm d}V.\] That is to say, the rate of change of \(f\rho \, \delta V\) following the element is just \(\mathchoice{\frac{\mathrm{D} f}{\mathrm{D} t}}{\mathrm{D} f/\mathrm{D} t}{\mathrm{D} f/\mathrm{D} t}{\mathrm{D} f/\mathrm{D} t}\) multiplied by the mass \(\rho \, \delta V\), because the mass of the element is conserved. It can be generalised to vectors by applying it to each vector element in turn.

Applying the transport theorem to the left hand side of [eqn:angm] gives \[\begin{align} \mathchoice{\frac{{\mathrm d}}{{\mathrm d}t}}{{\mathrm d}/{\mathrm d}t}{{\mathrm d}/{\mathrm d}t}{{\mathrm d}/{\mathrm d}t}\int_{V(t)}\mathbfit{x}\times\rho\mathbfit{u}\;\mathrm{d}V&= \int_{V(t)}\mathchoice{\frac{\mathrm{D} }{\mathrm{D} t}}{\mathrm{D} /\mathrm{D} t}{\mathrm{D} /\mathrm{D} t}{\mathrm{D} /\mathrm{D} t}\big[\mathbfit{x}\times\mathbfit{u}\big]\rho \;\mathrm{d}V\\ &= \int_{V(t)}\bigg[\underbrace{\mathchoice{\frac{\mathrm{D} \mathbfit{x}}{\mathrm{D} t}}{\mathrm{D} \mathbfit{x}/\mathrm{D} t}{\mathrm{D} \mathbfit{x}/\mathrm{D} t}{\mathrm{D} \mathbfit{x}/\mathrm{D} t}\times\mathbfit{u}}_{\mathbfit{u}\times\mathbfit{u}=\mathbf{0}}{} + \mathbfit{x}\times\mathchoice{\frac{\mathrm{D} \mathbfit{u}}{\mathrm{D} t}}{\mathrm{D} \mathbfit{u}/\mathrm{D} t}{\mathrm{D} \mathbfit{u}/\mathrm{D} t}{\mathrm{D} \mathbfit{u}/\mathrm{D} t}\bigg]\rho \;\mathrm{d}V\\ &= \int_{V(t)}\mathbfit{x}\times\rho\mathchoice{\frac{\mathrm{D} \mathbfit{u}}{\mathrm{D} t}}{\mathrm{D} \mathbfit{u}/\mathrm{D} t}{\mathrm{D} \mathbfit{u}/\mathrm{D} t}{\mathrm{D} \mathbfit{u}/\mathrm{D} t}\;\mathrm{d}V. \label{eqn:stress1} \end{align}\]

Now consider the stress term in [eqn:angm]. Applying the divergence theorem ‘backwards’ gives \[\begin{align} \int_{S(t)}\big[\mathbfit{x}\times(\normalsize{𝞼}\cdot\,{\mathrm d}\mathbfit{S})\big]_i &= \int_{S(t)}\varepsilon_{ijk}x_j\sigma_{kl}\widehat{n}_l\,{\mathrm d}S\\ &= \int_{V(t)}\mathchoice{\frac{\partial}{\partial x_l}}{\partial/\partial x_l}{\partial/\partial x_l}{\partial/\partial x_l}\big[\varepsilon_{ijk}x_j\sigma_{kl}\big]\;\mathrm{d}V\\ &= \int_{V(t)}\varepsilon_{ijk}\left[\mathchoice{\frac{\partial x_j}{\partial x_l}}{\partial x_j/\partial x_l}{\partial x_j/\partial x_l}{\partial x_j/\partial x_l}\sigma_{kl}+ x_j\mathchoice{\frac{\partial\sigma_{kl}}{\partial x_l}}{\partial\sigma_{kl}/\partial x_l}{\partial\sigma_{kl}/\partial x_l}{\partial\sigma_{kl}/\partial x_l}\right]\;\mathrm{d}V\\ &= \int_{V(t)}\varepsilon_{ijk}\left[\delta_{jl}\sigma_{kl}+ x_j(\bm{\nabla}\cdot\normalsize{𝞼}^\mathsf{T})_k\right]\;\mathrm{d}V\\ &= \int_{V(t)}\left[\varepsilon_{ijk}\sigma_{kj} + \varepsilon_{ijk}x_j(\bm{\nabla}\cdot\normalsize{𝞼}^\mathsf{T})_k\right]\;\mathrm{d}V. \label{eqn:stress2} \end{align}\] Putting [eqn:stress1] and [eqn:stress2] into [eqn:angm] gives \[\int_{V(t)}\varepsilon_{ijk}x_j\left[\rho\mathchoice{\frac{\mathrm{D} u_k}{\mathrm{D} t}}{\mathrm{D} u_k/\mathrm{D} t}{\mathrm{D} u_k/\mathrm{D} t}{\mathrm{D} u_k/\mathrm{D} t} - (\bm{\nabla}\cdot\normalsize{𝞼}^\mathsf{T})_k - \rho f_k \right]\;\mathrm{d}V= \int_{V(t)}\varepsilon_{ijk}\sigma_{kj}\;\mathrm{d}V.\] The left-hand side vanishes by the Cauchy momentum equation, so angular momentum is conserved if and only if \[\int_{V(t)}\varepsilon_{ijk}\sigma_{kj}\;\mathrm{d}V= 0.\] This needs to hold for arbitrary \(V(t)\), so \[\begin{aligned} \varepsilon_{ijk}\sigma_{kj} = 0 \quad &\iff \quad \varepsilon_{imn}\varepsilon_{ijk}\sigma_{kj} = 0\\ &\iff \quad \big(\delta_{mj}\delta_{nk} - \delta_{mk}\delta_{nj}\big)\sigma_{kj} = 0\\ &\iff \quad \sigma_{nm} - \sigma_{mn} = 0 \quad \iff \sigma_{ij} = \sigma_{ji}. \end{aligned}\]

Phew! We made it! Now we also want that \(\sigma_{ij}\) reduces to the inviscid expression in the absence of viscosity, so we write \[\sigma_{ij} = -p\delta_{ij} + \tau_{ij},\] where \(\mathsfbfit{\tau} = (\tau_{ij})\) is called the deviatoric stress tensor. Since the diagonal components of \(\normalsize{𝞼}\) need no longer be equal, we will simply define the pressure \(p\) as the mean normal stress, \[p = -\frac{1}{3}\sigma_{ii} = -\frac{1}{3}(\sigma_{11}+\sigma_{22}+\sigma_{33}).\] It follows that the deviatoric part is traceless, \[\begin{equation} \tau_{ii} = 0. \label{eqn:dii} \end{equation}\]

We now have our general form for the elements of \(\normalsize{𝞼}\) so that each element depends on the other in an allowed physical way and in a way that reduces to the ideal case when the shear stresses are zero. Now we need a form for the shear stresses themselves.

6.2.1 Simple shear

To get an idea of what this means, let’s consider a shear flow, \[\mathbfit{u} = U(y)\widehat{\mathbfit{e}}_x,\] with \(U'(y) > 0\).

Consider the stress on the surface \(S\) given by \(y=y_0\). In an ideal fluid, the forces would all be normal to \(S\), so there would be no transfer of momentum across \(S\). But in reality, kinetic theory dictates that some molecules will diffuse across \(S\), transferring momentum. If a molecule moves from \(y<y_0\) to \(y>y_0\), it must be accelerated, so there must be a force in the \(x\)-direction. In a Newtonian fluid, we approximate the shear stress by the linear relation \[\sigma_{12} = \mu\mathchoice{\frac{{\mathrm d}U}{{\mathrm d}y}}{{\mathrm d}U/{\mathrm d}y}{{\mathrm d}U/{\mathrm d}y}{{\mathrm d}U/{\mathrm d}y},\] where the constant of proportionality \(\mu\) is the viscosity.

Since \(\widehat{\mathbfit{n}}=\widehat{\mathbfit{e}}_y\) for the volume \(y<y_0\), the force \(\normalsize{𝞼}\cdot\widehat{\mathbfit{n}}\) on the lower fluid is to the right (speeding it up). For the volume \(y>y_0\), the force is to the left (slowing it down) since \(\widehat{\mathbfit{n}}=-\widehat{\mathbfit{e}}_y\). So viscosity opposes relative motion between different parts of the fluid. Effectively it adds a drag force.

6.2.2 Shear in all directions

In general, we need to consider the viscous force exerted by all three components of \(\mathbfit{u}\) on surfaces of arbitrary orientation. So a Newtonian fluid is defined by the linear relation \[\begin{equation} \tau_{ij} = A_{ijkl}\mathchoice{\frac{\partial u_l}{\partial x_k}}{\partial u_l/\partial x_k}{\partial u_l/\partial x_k}{\partial u_l/\partial x_k} \label{eqn:newt} \end{equation}\] for some rank 4 tensor \(A_{ijkl}\). You can think of a rank 4 tensor as a matrix of matrices, or simply just that for every combination of \(ijkl\) there will be a specific number. We make the further physical assumption that \(A_{ijkl}\) must be isotropic.

The assumption that \(A_{ijkl}\) is isotropic means that it must take the form \[A_{ijkl} = \alpha\delta_{ij}\delta_{kl} + \beta\delta_{ik}\delta_{jl} + \gamma\delta_{il}\delta_{jk}.\] This is a general result for rank 4 tensors and can be proved using rotation matrices or geometrical arguments. For our purposes we’ll just use it directly. This means then that \[\begin{aligned} \tau_{ij} &= \alpha\delta_{ij}\delta_{kl}\mathchoice{\frac{\partial u_l}{\partial x_k}}{\partial u_l/\partial x_k}{\partial u_l/\partial x_k}{\partial u_l/\partial x_k}+ \beta\delta_{ik}\delta_{jl} \mathchoice{\frac{\partial u_l}{\partial x_k}}{\partial u_l/\partial x_k}{\partial u_l/\partial x_k}{\partial u_l/\partial x_k}+\gamma\delta_{il}\delta_{jk}\mathchoice{\frac{\partial u_l}{\partial x_k}}{\partial u_l/\partial x_k}{\partial u_l/\partial x_k}{\partial u_l/\partial x_k}, \\ &= \alpha\mathchoice{\frac{\partial u_k}{\partial x_k}}{\partial u_k/\partial x_k}{\partial u_k/\partial x_k}{\partial u_k/\partial x_k}\delta_{ij} + \beta\mathchoice{\frac{\partial u_j}{\partial x_i}}{\partial u_j/\partial x_i}{\partial u_j/\partial x_i}{\partial u_j/\partial x_i} + \gamma \mathchoice{\frac{\partial u_i}{\partial x_j}}{\partial u_i/\partial x_j}{\partial u_i/\partial x_j}{\partial u_i/\partial x_j}. \end{aligned}\] Since \(\normalsize{𝞼}\) is symmetric we must also have that \(\tau_{ji}=\tau_{ij}\). This implies that \(\gamma=\beta\) so \[\tau_{ij} = \alpha\mathchoice{\frac{\partial u_k}{\partial x_k}}{\partial u_k/\partial x_k}{\partial u_k/\partial x_k}{\partial u_k/\partial x_k}\delta_{ij} + \beta\left(\mathchoice{\frac{\partial u_j}{\partial x_i}}{\partial u_j/\partial x_i}{\partial u_j/\partial x_i}{\partial u_j/\partial x_i} + \mathchoice{\frac{\partial u_i}{\partial x_j}}{\partial u_i/\partial x_j}{\partial u_i/\partial x_j}{\partial u_i/\partial x_j}\right).\] Finally, applying \(\tau_{ii}=0\) gives \[0 = 3\alpha\mathchoice{\frac{\partial u_k}{\partial x_k}}{\partial u_k/\partial x_k}{\partial u_k/\partial x_k}{\partial u_k/\partial x_k} + 2\beta\mathchoice{\frac{\partial u_i}{\partial x_i}}{\partial u_i/\partial x_i}{\partial u_i/\partial x_i}{\partial u_i/\partial x_i},\] which implies \(3\alpha=-2\beta\) since \(\mathchoice{\frac{\partial u_k}{\partial x_k}}{\partial u_k/\partial x_k}{\partial u_k/\partial x_k}{\partial u_k/\partial x_k}\) is just a relabelling of \(\mathchoice{\frac{\partial u_i}{\partial x_i}}{\partial u_i/\partial x_i}{\partial u_i/\partial x_i}{\partial u_i/\partial x_i}\). Writing \(\mu=\beta\) gives \[\tau_{ij} = \mu\left(\mathchoice{\frac{\partial u_j}{\partial x_i}}{\partial u_j/\partial x_i}{\partial u_j/\partial x_i}{\partial u_j/\partial x_i} + \mathchoice{\frac{\partial u_i}{\partial x_j}}{\partial u_i/\partial x_j}{\partial u_i/\partial x_j}{\partial u_i/\partial x_j} - \frac{2}{3}\mathchoice{\frac{\partial u_k}{\partial x_k}}{\partial u_k/\partial x_k}{\partial u_k/\partial x_k}{\partial u_k/\partial x_k}\delta_{ij}\right).\] Thus, we finally have that the stress tensor for a Newtonian fluid takes the form \[\boxed{\sigma_{ij} = -p\delta_{ij} + \mu\left[\mathchoice{\frac{\partial u_j}{\partial x_i}}{\partial u_j/\partial x_i}{\partial u_j/\partial x_i}{\partial u_j/\partial x_i} + \mathchoice{\frac{\partial u_i}{\partial x_j}}{\partial u_i/\partial x_j}{\partial u_i/\partial x_j}{\partial u_i/\partial x_j} - \frac{2}{3}(\bm{\nabla}\cdot\mathbfit{u})\delta_{ij} \right],}\] where the constant \(\mu\) is called the viscosity. You might recognise the first two terms in the bracket as being equal to twice \(\mathsfbfit{E}\), the rate of strain tensor from 1.6.

6.4.1 Diffusion of vorticity

To see this, let’s take the curl of the momentum equation, [eqn:nscons]. Since the curl commutes with the Laplacian, and using our favourite vector identity to re-express the nonlinear velocity term, we get \[\begin{equation} \boxed{ \frac{\partial\mathbfit{\omega}}{\partial t} + \bm{\nabla}\times\big(\mathbfit{\omega}\times\mathbfit{u}\big) = \underbrace{\nu{\nabla^2}\mathbfit{\omega}}_{\text{diffusion}}}. \label{eqn:nsvort} \end{equation}\] This is the vorticity equation from 2.4, but now accounting for viscosity. Compare it with the inviscid version, [eqn:vort]: there, the left-hand side was unchanged, but the right-hand side was zero. We interpreted that to mean fluid elements with no vorticity do not gain any as they move. So, where does vorticity come from?

The answer lies in the new right-hand side. The Laplacian represents diffusion, scaled by the kinematic viscosity \(\nu\). This means vorticity spreads because of viscosity! Recall from Michaelmas that vorticity is local – it describes the rotation of a fluid element, moving as the element moves. Over time, the diffusion term spreads vorticity between elements and, eventually, throughout the fluid.

To illustrate this, we consider an example where vorticity is suddenly generated at one location and the motion of vorticity is purely through diffusion.

6.4.2 Viscous dissipation

Another way to characterise the effect of viscosity on a flow is to look at its effect on kinetic energy. Again, in a general sense, viscosity acts to decay and spread out the kinetic energy. The incompressible Navier–Stokes equations with conservative body force take the form \[\begin{align} \bm{\nabla}\cdot\mathbfit{u} &= 0,\label{eqn:nscons0}\\ \mathchoice{\frac{\partial\mathbfit{u}}{\partial t}}{\partial\mathbfit{u}/\partial t}{\partial\mathbfit{u}/\partial t}{\partial\mathbfit{u}/\partial t} + \big(\mathbfit{u}\cdot\bm{\nabla}\big)\mathbfit{u} &= - \frac{1}{\rho_0}\bm{\nabla}p + \nu{\nabla^2}\mathbfit{u} -\bm{\nabla}\mathit{\Phi}. \label{eqn:nscons-conservative} \end{align}\] Recall from last term, in 2.3.1, that kinetic energy is given by \[E=\frac{1}{2}\int_V \rho_0|\mathbfit{u}|^2\;\mathrm{d}V,\] and the corresponding Euler equations (\(\nu=0\)) conserve kinetic energy in a fixed, closed container, \[\mathchoice{\frac{{\mathrm d}E}{{\mathrm d}t}}{{\mathrm d}E/{\mathrm d}t}{{\mathrm d}E/{\mathrm d}t}{{\mathrm d}E/{\mathrm d}t} = 0.\] The presence of viscosity \(\nu\neq 0\), however, leads to the dissipation of this energy. Let’s see what that looks like.

OK, time for one more round of vector calculus. Proceeding as we did for the Euler equations in 2.3.1, we have \[\begin{aligned} \mathchoice{\frac{{\mathrm d}E}{{\mathrm d}t}}{{\mathrm d}E/{\mathrm d}t}{{\mathrm d}E/{\mathrm d}t}{{\mathrm d}E/{\mathrm d}t} &=\rho_0 \mathchoice{\frac{{\mathrm d}}{{\mathrm d}t}}{{\mathrm d}/{\mathrm d}t}{{\mathrm d}/{\mathrm d}t}{{\mathrm d}/{\mathrm d}t} \int_V \frac12 |\mathbfit{u}|^2 \, {\mathrm d}V \\ &= \rho_0\int_V\mathbfit{u}\cdot\mathchoice{\frac{\partial\mathbfit{u}}{\partial t}}{\partial\mathbfit{u}/\partial t}{\partial\mathbfit{u}/\partial t}{\partial\mathbfit{u}/\partial t}\;\mathrm{d}V\quad \text{(because $V$ is fixed)}\\ &= \rho_0\int_V\mathbfit{u}\cdot\left[- \frac{1}{\rho_0}\bm{\nabla}p - \big(\mathbfit{u}\cdot\bm{\nabla}\big)\mathbfit{u} + \nu{\nabla^2}\mathbfit{u} -\bm{\nabla}\mathit{\Phi} \right]\;\mathrm{d}V,\\ \end{aligned}\]then by the famous vector identity [u-dot-grad-u],\[\begin{aligned} &= - \rho_0\int_V\mathbfit{u}\cdot\bm{\nabla}\underbrace{\left[\frac{p}{\rho_0} + \frac12|\mathbfit{u}|^2 + \mathit{\Phi}\right]}_H\;\mathrm{d}V+ \mu\int_V\mathbfit{u}\cdot{\nabla^2}\mathbfit{u}\;\mathrm{d}V, \end{aligned}\]which by a reversed product rule and divergence theorem as in [almost-bern],\[\begin{aligned} &= -\rho_0\int_{S}H\mathbfit{u}\cdot\,{\mathrm d}\mathbfit{S}+ \mu\int_V\mathbfit{u}\cdot{\nabla^2}\mathbfit{u}\;\mathrm{d}V, \end{aligned}\]then since the volume is closed, we can once again say \(\mathbfit{u}\cdot\,{\mathrm d}\mathbfit{S}=\mathbfit{u}\cdot\widehat{\mathbfit{n}}\,{\mathrm d}S = 0\) on \(S\),\[\begin{aligned} &= \mu\int_V\mathbfit{u}\cdot{\nabla^2}\mathbfit{u}\;\mathrm{d}V. \end{aligned}\] In viscous fluids, though, we can say something stronger than the no-flux condition at the boundary: we have the no-slip condition. This will be handy shortly. For now, using the rearranged product rule (aka integration by parts), \[\begin{aligned} \int_V\mathbfit{u}\cdot{\nabla^2}\mathbfit{u}\;\mathrm{d}V&= \int_V u_k\mathchoice{\frac{\partial^2 u_k}{\partial x_j^2}}{\partial^2 u_k/\partial x_j^2}{\partial^2 u_k/\partial x_j^2}{\partial^2 u_k/\partial x_j^2}\;\mathrm{d}V,\\ &= \int_V\left[\mathchoice{\frac{\partial}{\partial x_j}}{\partial/\partial x_j}{\partial/\partial x_j}{\partial/\partial x_j}\left(u_k\mathchoice{\frac{\partial u_k}{\partial x_j}}{\partial u_k/\partial x_j}{\partial u_k/\partial x_j}{\partial u_k/\partial x_j}\right) - \mathchoice{\frac{\partial u_k}{\partial x_j}}{\partial u_k/\partial x_j}{\partial u_k/\partial x_j}{\partial u_k/\partial x_j}\mathchoice{\frac{\partial u_k}{\partial x_j}}{\partial u_k/\partial x_j}{\partial u_k/\partial x_j}{\partial u_k/\partial x_j} \right]\;\mathrm{d}V,\\ &= \int_{S}u_k\mathchoice{\frac{\partial u_k}{\partial x_j}}{\partial u_k/\partial x_j}{\partial u_k/\partial x_j}{\partial u_k/\partial x_j}\widehat{n}_j\,{\mathrm d}S - \int_V\mathchoice{\frac{\partial u_k}{\partial x_j}}{\partial u_k/\partial x_j}{\partial u_k/\partial x_j}{\partial u_k/\partial x_j}\mathchoice{\frac{\partial u_k}{\partial x_j}}{\partial u_k/\partial x_j}{\partial u_k/\partial x_j}{\partial u_k/\partial x_j}\;\mathrm{d}V, \quad\text{(divergence theorem)}\\ &= -\int_V\mathchoice{\frac{\partial u_k}{\partial x_j}}{\partial u_k/\partial x_j}{\partial u_k/\partial x_j}{\partial u_k/\partial x_j}\mathchoice{\frac{\partial u_k}{\partial x_j}}{\partial u_k/\partial x_j}{\partial u_k/\partial x_j}{\partial u_k/\partial x_j}\;\mathrm{d}V, \end{aligned}\] using the no-slip boundary condition \(\mathbfit{u}=\mathbf{0}\) on \(S\). What is this thing we’re left with? Well, it’s actually equivalent to \[\begin{aligned} \int_V\big|\bm{\nabla}\times\mathbfit{u}\big|^2\;\mathrm{d}V&= \int_V\varepsilon_{ijk}\varepsilon_{ilm}\mathchoice{\frac{\partial u_k}{\partial x_j}}{\partial u_k/\partial x_j}{\partial u_k/\partial x_j}{\partial u_k/\partial x_j}\mathchoice{\frac{\partial u_m}{\partial x_l}}{\partial u_m/\partial x_l}{\partial u_m/\partial x_l}{\partial u_m/\partial x_l}\;\mathrm{d}V\\ &= \int_V\Big(\delta_{jl}\delta_{km} - \delta_{jm}\delta_{kl} \Big)\mathchoice{\frac{\partial u_k}{\partial x_j}}{\partial u_k/\partial x_j}{\partial u_k/\partial x_j}{\partial u_k/\partial x_j}\mathchoice{\frac{\partial u_m}{\partial x_l}}{\partial u_m/\partial x_l}{\partial u_m/\partial x_l}{\partial u_m/\partial x_l}\;\mathrm{d}V\\ &= \int_V\left(\mathchoice{\frac{\partial u_k}{\partial x_j}}{\partial u_k/\partial x_j}{\partial u_k/\partial x_j}{\partial u_k/\partial x_j}\mathchoice{\frac{\partial u_k}{\partial x_j}}{\partial u_k/\partial x_j}{\partial u_k/\partial x_j}{\partial u_k/\partial x_j} - \mathchoice{\frac{\partial u_k}{\partial x_j}}{\partial u_k/\partial x_j}{\partial u_k/\partial x_j}{\partial u_k/\partial x_j}\mathchoice{\frac{\partial u_j}{\partial x_k}}{\partial u_j/\partial x_k}{\partial u_j/\partial x_k}{\partial u_j/\partial x_k} \right)\;\mathrm{d}V\\ &= \int_V\mathchoice{\frac{\partial u_k}{\partial x_j}}{\partial u_k/\partial x_j}{\partial u_k/\partial x_j}{\partial u_k/\partial x_j}\mathchoice{\frac{\partial u_k}{\partial x_j}}{\partial u_k/\partial x_j}{\partial u_k/\partial x_j}{\partial u_k/\partial x_j}\;\mathrm{d}V- \int_V\left[\mathchoice{\frac{\partial}{\partial x_j}}{\partial/\partial x_j}{\partial/\partial x_j}{\partial/\partial x_j}\left(u_k\mathchoice{\frac{\partial u_j}{\partial x_k}}{\partial u_j/\partial x_k}{\partial u_j/\partial x_k}{\partial u_j/\partial x_k} \right) - u_k\mathchoice{\frac{\partial^2 u_j}{\partial x_j \partial x_k}}{\partial^2 u_j /\partial x_j \partial x_k}{\partial^2 u_j /\partial x_j \partial x_k}{\partial^2 u_j /\partial x_j \partial x_k} \right]\;\mathrm{d}V,\\ &= \int_V\mathchoice{\frac{\partial u_k}{\partial x_j}}{\partial u_k/\partial x_j}{\partial u_k/\partial x_j}{\partial u_k/\partial x_j}\mathchoice{\frac{\partial u_k}{\partial x_j}}{\partial u_k/\partial x_j}{\partial u_k/\partial x_j}{\partial u_k/\partial x_j}\;\mathrm{d}V- \int_{S}u_k\mathchoice{\frac{\partial u_j}{\partial x_k}}{\partial u_j/\partial x_k}{\partial u_j/\partial x_k}{\partial u_j/\partial x_k}\widehat{n}_j\,{\mathrm d}S + \int_V u_k\mathchoice{\frac{\partial}{\partial x_k}}{\partial/\partial x_k}{\partial/\partial x_k}{\partial/\partial x_k}\Big(\bm{\nabla}\cdot\mathbfit{u}\Big)\;\mathrm{d}V,\\ &= \int_V\mathchoice{\frac{\partial u_k}{\partial x_j}}{\partial u_k/\partial x_j}{\partial u_k/\partial x_j}{\partial u_k/\partial x_j}\mathchoice{\frac{\partial u_k}{\partial x_j}}{\partial u_k/\partial x_j}{\partial u_k/\partial x_j}{\partial u_k/\partial x_j}\;\mathrm{d}V. \end{aligned}\] In the last step the boundary term vanished again since \(\mathbfit{u}=\mathbf{0}\) on \(S\), and the last term vanished by incompressibility (\(\bm{\nabla}\cdot\mathbfit{u}=0\)).

So in conclusion, in an incompressible viscous fluid with \(\mathbfit{f}=-\bm{\nabla}\mathit{\Phi}\) and a stationary boundary, \[\mathchoice{\frac{{\mathrm d}E}{{\mathrm d}t}}{{\mathrm d}E/{\mathrm d}t}{{\mathrm d}E/{\mathrm d}t}{{\mathrm d}E/{\mathrm d}t} = -\mu\int_V\left|\mathbfit{\omega}\right|^2\,{\mathrm d}V.\]

6.5.1 High Reynolds number limit

As \(\textit{Re}\to\infty\) the viscous term in [eqn:rey] becomes smaller and smaller. However, this is a singular limit: in other words, the solution of the Navier–Stokes equation does not (in general) tend to the solution of the Euler equation. This is because the viscous term has the highest derivative.

In the fluid case, there is an analogous layer above a solid boundary, where \(\mathbfit{u}\) has to depart from an inviscid flow in order to satisfy the no-slip condition. This is called a boundary layer. We already saw an example with the impulsively moved boundary. In that example, we saw that the boundary layer thickness grew like \(\sqrt{\nu t}\) with the fluid in the volume ‘feeling’ the effect of the moved boundary once the boundary layer propagated past it. More often we have the reverse situation as shown below.

Here we have a flow past a stationary object that slows down near the surface and stops at the surface itself (no slip). In that case it is the reduction in speed that propagates up from the surface within the boundary layer (of thickness \(\sim \sqrt{\nu t}\)).

6.5.2 Low Reynolds number limit

At the other extreme, where \(\textit{Re}\ll 1\), viscosity dominates. This applies to flows with extremely small length scales (e.g. biological applications like bacteria) or very viscous liquids like honey.

When \(\textit{Re}\ll 1\), we neglect \(\displaystyle\mathchoice{\frac{\mathrm{D} \mathbfit{u}}{\mathrm{D} t}}{\mathrm{D} \mathbfit{u}/\mathrm{D} t}{\mathrm{D} \mathbfit{u}/\mathrm{D} t}{\mathrm{D} \mathbfit{u}/\mathrm{D} t}\) and obtain the (unforced) Stokes equations, \[\boxed{ \begin{gathered} %\begin{align} \bm{\nabla}\cdot\mathbfit{u} = 0,\\ \mu{\nabla^2}\mathbfit{u} = \bm{\nabla}p. %\end{align} \end{gathered}}\] We still have the no-slip boundary condition.

The Stokes equations are linear, so are much easier to deal with (I recommend them). They also only admit one unique solution, \(\mathbfit{u}(\mathbfit{x})\), \(p(\mathbfit{x})\) in a domain \(V\) matching a given boundary condition \(\mathbfit{u}=\mathbfit{u}_S(\mathbfit{x})\) on \(S\). The boundary therefore drives the flow. The pressure \(p\) is unique up to an additive constant.

To see this, suppose there exist two solutions \(\mathbfit{u}_1\), \(\mathbfit{u}_2\) and consider the difference \(\mathbfit{v}=\mathbfit{u}_1-\mathbfit{u}_2\) with corresponding pressure difference \(q=p_1-p_2\). By linearity, we have \[\mu{\nabla^2}\mathbfit{v} = \bm{\nabla}q, \qquad \bm{\nabla}\cdot\mathbfit{v} = 0,\] and since \(\mathbfit{u}_1\) and \(\mathbfit{u}_2\) satisfy the same boundary condition, we have \(\mathbfit{v}=\mathbf{0}\) on \(S\). It follows that \[\begin{aligned} \int_V\mathbfit{v}\cdot\Big[\mu{\nabla^2}\mathbfit{v} - \bm{\nabla}q\Big]\;\mathrm{d}V&=0\\ \implies \quad \int_Vv_i\left[\mu\mathchoice{\frac{\partial^2 v_i}{\partial x_j \partial x_j}}{\partial^2 v_i /\partial x_j \partial x_j}{\partial^2 v_i /\partial x_j \partial x_j}{\partial^2 v_i /\partial x_j \partial x_j} - \mathchoice{\frac{\partial q}{\partial x_i}}{\partial q/\partial x_i}{\partial q/\partial x_i}{\partial q/\partial x_i}\right]\;\mathrm{d}V&= 0. \end{aligned}\] Writing both terms as results of the product rule, we get \[\begin{aligned} \mu\int_V\mathchoice{\frac{\partial}{\partial x_j}}{\partial/\partial x_j}{\partial/\partial x_j}{\partial/\partial x_j}\left[v_i\mathchoice{\frac{\partial v_i}{\partial x_j}}{\partial v_i/\partial x_j}{\partial v_i/\partial x_j}{\partial v_i/\partial x_j}\right]\;\mathrm{d}V- \mu\int_V\mathchoice{\frac{\partial v_i}{\partial x_j}}{\partial v_i/\partial x_j}{\partial v_i/\partial x_j}{\partial v_i/\partial x_j}\mathchoice{\frac{\partial v_i}{\partial x_j}}{\partial v_i/\partial x_j}{\partial v_i/\partial x_j}{\partial v_i/\partial x_j}\;\mathrm{d}V- \int_V\mathchoice{\frac{\partial}{\partial x_i}}{\partial/\partial x_i}{\partial/\partial x_i}{\partial/\partial x_i}\Big[v_i q\Big]\;\mathrm{d}V&= 0, \end{aligned}\] noting that we haven’t forgotten one of the contributions from the final term’s product rule expansion: it just equals zero, \(\mathchoice{\frac{\partial v_i}{\partial x_i}}{\partial v_i/\partial x_i}{\partial v_i/\partial x_i}{\partial v_i/\partial x_i}=\bm{\nabla}\cdot\mathbfit{v}=0\).

Applying the divergence theorem to the first and third terms gives \[\begin{aligned} \mu\int_{S}v_i\mathchoice{\frac{\partial v_i}{\partial x_j}}{\partial v_i/\partial x_j}{\partial v_i/\partial x_j}{\partial v_i/\partial x_j}\widehat{n}_j\,{\mathrm d}S - \mu\int_V\mathchoice{\frac{\partial v_i}{\partial x_j}}{\partial v_i/\partial x_j}{\partial v_i/\partial x_j}{\partial v_i/\partial x_j}\mathchoice{\frac{\partial v_i}{\partial x_j}}{\partial v_i/\partial x_j}{\partial v_i/\partial x_j}{\partial v_i/\partial x_j}\;\mathrm{d}V- \int_{S}qv_i\widehat{n}_i\,{\mathrm d}S &= 0\\ \implies\quad \mu\int_V\mathchoice{\frac{\partial v_i}{\partial x_j}}{\partial v_i/\partial x_j}{\partial v_i/\partial x_j}{\partial v_i/\partial x_j}\mathchoice{\frac{\partial v_i}{\partial x_j}}{\partial v_i/\partial x_j}{\partial v_i/\partial x_j}{\partial v_i/\partial x_j}\;\mathrm{d}V&= 0, \end{aligned}\] where we used the boundary condition \(\mathbfit{v}=\mathbf{0}\). The remaining integrand is \[\begin{aligned} \mathchoice{\frac{\partial v_i}{\partial x_j}}{\partial v_i/\partial x_j}{\partial v_i/\partial x_j}{\partial v_i/\partial x_j}\mathchoice{\frac{\partial v_i}{\partial x_j}}{\partial v_i/\partial x_j}{\partial v_i/\partial x_j}{\partial v_i/\partial x_j} &= \left(\mathchoice{\frac{\partial v_x}{\partial x}}{\partial v_x/\partial x}{\partial v_x/\partial x}{\partial v_x/\partial x}\right)^2 + \left(\mathchoice{\frac{\partial v_x}{\partial y}}{\partial v_x/\partial y}{\partial v_x/\partial y}{\partial v_x/\partial y}\right)^2+ \left(\mathchoice{\frac{\partial v_x}{\partial z}}{\partial v_x/\partial z}{\partial v_x/\partial z}{\partial v_x/\partial z}\right)^2 + \left(\mathchoice{\frac{\partial v_y}{\partial x}}{\partial v_y/\partial x}{\partial v_y/\partial x}{\partial v_y/\partial x}\right)^2 + \left(\mathchoice{\frac{\partial v_y}{\partial y}}{\partial v_y/\partial y}{\partial v_y/\partial y}{\partial v_y/\partial y}\right)^2+ \left(\mathchoice{\frac{\partial v_y}{\partial z}}{\partial v_y/\partial z}{\partial v_y/\partial z}{\partial v_y/\partial z}\right)^2 \nonumber\\ &\quad+ \left(\mathchoice{\frac{\partial v_z}{\partial x}}{\partial v_z/\partial x}{\partial v_z/\partial x}{\partial v_z/\partial x}\right)^2 + \left(\mathchoice{\frac{\partial v_z}{\partial y}}{\partial v_z/\partial y}{\partial v_z/\partial y}{\partial v_z/\partial y}\right)^2+ \left(\mathchoice{\frac{\partial v_z}{\partial z}}{\partial v_z/\partial z}{\partial v_z/\partial z}{\partial v_z/\partial z}\right)^2. \end{aligned}\] So the only way for the integral to vanish is if \(\mathbfit{v}\) is constant. And since \(\mathbfit{v}=\mathbf{0}\) on the boundary, \(\mathbfit{v} =\mathbf{0}\) everywhere. It follows that \(\mathbfit{u}_1=\mathbfit{u}_2\), hence \(\bm{\nabla}q=0\) and \(p_1=p_2+c\).

The key observation goes like this: if your boundary condition is \(\mathbfit{u}_S = \mathbf{0}\), clearly \(\mathbfit{u}=\mathbf{0}\) is a solution to the unforced Stokes equations. But by what we’ve just shown, it’s the only solution! So if we want non-zero flow, we therefore need \(\mathbfit{u}_S\neq \mathbf{0}\), which implies that the boundary itself is moving.