\documentclass[12pt,a4paper]{report}%{book}
\usepackage{graphicx,amsmath,amssymb,xcolor}
%\usepackage{epstopdf}
\usepackage{makeidx}

\newcommand{\bfx}{{\bf x}}
\newcommand{\bfu}{{\bf u}}
\newcommand{\bfU}{{\bf U}}
\newcommand{\bfn}{{\bf n}}
\newcommand{\bfF}{{\bf F}}
\newcommand{\bfom}{{\boldsymbol\omega}}
\newcommand{\qword}[1]{\quad\hbox{\rm #1}\quad}
\newcommand{\pd}[2]{\frac{\partial #1}{\partial #2}}
\DeclareMathOperator\erfc{{\rm erfc}}
\DeclareMathOperator{\sech}{sech}
\newcommand{\cpicture}[2]{\begin{center}\includegraphics[width=#1\linewidth]{Figures/#2}\end{center}}
\newtheorem{ex}{Exercise}[chapter]

%\makeindex

\title{Fluid Dynamics II (D)}
\author{John Hinch}
\date{Draft 1, \quad \today}

\begin{document}
%\frontmatter

\maketitle

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\tableofcontents


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\chapter*{Preface}

These lecture notes are for the third-year course on fluid dynamics
in the Mathematical Tripos at Cambridge University.
They were developed while lecturing the course over the years on ten occasions.
The current lecturer may well present material differently and in a different order, which could be better.

The third-year course follows a second-year course on fluid dynamics.
Key results from that course are summarised in chapter \ref{chpt:rev}.
That second-year course included friction but only for unidirectional flows.
This course will generalise viscous friction to full three-dimensional flows, in chapter \ref{chpt:NS}.
Unidirectional flows are then revisited in the following chapter \ref{chpt:uni}, examining finer details than in the previous course.

The equations governing the motion of fluids, the Navier-Stokes equations, are non-linear partial differential equations, with few exact analytic solutions.
There is a major industry in solving the equations numerically, but that is not for this course.
To progress, this course makes approximations which are appropriate to various limiting conditions, and produce accurate answers in those conditions.
Chapter \ref{chpt:stokes} considers flows in which viscous friction dominates inertia, three-dimensional flows rather than the unidirectional flows of the previous course.
An important subclass of flows with negligible inertia occurs in thin layers, where locally the flow looks nearly unidirectional but globally the flow is three-dimensional, see chapter \ref{chpt:thinlayer}.
The opposite limit of large inertia and small viscosity is more subtle:
in most of the flow viscosity can be neglected, but there are thin layers, normally next to boundaries, where viscosity must be considered.
The idea of the effects of viscosity being confined to thin regions is first developed in some special flows in chapter \ref{chpt:vort}.
This is followed in chapter \ref{chpt:bl} by a study of viscous boundary layers.
Finally, many possible flows are unstable.
The course ends by examining the stability of a steady laminar flow in chapter \ref{chpt:stab}.


%\mainmatter
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\chapter{Revision of ideas and results from Fluids in IB}
\label{chpt:rev}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Continuum hypothesis}

We do not deal with the dynamics of individual molecules.
Molecular details are smeared out by averaging over small volumes to define properties like
\begin{itemize}
 \item density $\rho(\bfx,t)$, 
 \item velocity $\bfu(\bfx,t)$, and 
 \item pressure $p(\bfx,t)$.
\end{itemize}
The averaging volume needs to be small compared with the laboratory scale of interest, but large compared with molecules.
The processing of smearing out leads to derived forces, volume forces and surface forces, see later.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Time derivatives}

A {\em fluid particle\/}, sometime called a {\em material element\/} or a {\em Lagrangian point\/}, moves with the fluid, 
so its position $\bfx(t)$ satisfies
\[
\dot\bfx = \bfu(\bfx,t) .
\]

\noindent The rate of change of a quantity as seen by a fluid particle is written
$D/Dt$, given by the chain rule as
\[
 \frac{D}{Dt} \equiv \frac{\partial}{\partial t} + \bfu\cdot\nabla .
\]

\noindent In particular, the {\em acceleration\/} of a fluid particle is
\[
 \frac{D\bfu}{Dt} = \frac{\partial\bfu}{\partial t} + \bfu\cdot\nabla\bfu .
\]


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Mass conservation}

\cpicture{0.6}{mass}
Note that one must constrain $\bfu(\bfx,t)$, otherwise $u=1$ into a pipe and $u=2$ out of the pipe creates mass.

\noindent Because mass in neither created or destroyed, the mass density $\rho$ satisfies
\[
 \frac{\partial\rho}{\partial t} + \nabla\cdot(\rho\bfu) = 0 ,
 \quad\hbox{\rm or equivalently}\quad
 \frac{D\rho}{Dt} + \rho\nabla\cdot\bfu = 0.
\]
Note that the first form of this equation is typical conservation laws.
The quantity $\rho\bfu$ is called the {\em mass flux\/}.
For an {\em incompressible\/} fluid, the density of each material element is constant and so $D\rho/Dt = 0$. Hence
\[
 \nabla\cdot\bfu = 0.
\]

For this course, we restrict attention to incompressible fluids  with uniform density, i.e.~$\rho$ independent of $\bfx$ and $t$.

For two-dimensional flows, the condition $\nabla\cdot\bfu=0$ is automatically
satisfied by $\bfu = (\psi_y,-\psi_x,0)$, with a {\em streamfunction}
$\psi(x,y)$.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Kinematic boundary condition}

\begin{minipage}{0.55\linewidth}
Applying mass conservation to a region close
 to a boundary $S$, we have
\[
 \rho\bfu^-\cdot\bfn = \rho\bfu^+\cdot\bfn ,
\]
\end{minipage}
\begin{minipage}{0.4\linewidth}
  \cpicture{0.9}{kbc}
\end{minipage}

\medskip
\noindent i.e.~the normal component of velocity must be continuous across $S$.
{\color{darkgray}[Not so at an evaporating interface where $\rho$ would jump.]}
In particular, at a fixed boundary, $\bfu\cdot\bfn=0$.
And if the moving boundary of a fluid is given by $F(\bfx,t)=0$, then since the surface consists of material points:
\[
DF/Dt=0. 
\]
This form of the boundary condition is sometime more convenient for free-surface problems.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Momentum conservation}

\begin{minipage}{0.5\linewidth}
On the  {\em assumption\/} that the  only surface force across a material surface $\bfn dS$
 is given by a {\em pressure\/} as $p(\bfx,t)$ as $-p \bfn dS$,
\end{minipage}
\begin{minipage}{0.4\linewidth}
    \cpicture{0.7}{ds}
\end{minipage}

\medskip
\noindent then Newton's equation of motion is
\[
 \rho \frac{D\bfu}{Dt} = -\nabla p + \bfF(\bfx,t) ,
\]
where $\bfF(\bfx,t)$ is the force per unit volume 
(e.g.~gravity $\rho{\bf g}$).
This is {\em Euler's equation\/}.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Dynamic boundary condition}

On the same assumption, applying momentum conservation to a region close to
the boundary $S$ gives (in the absence of surface tension) gives
\[
 - p^-\bfn = -p^+\bfn.
\]
Thus the pressure is continuous across $S$. {\color{darkgray}[Again not so at an evaporating interface.]}

\medskip
{\em In this course, we abandon the {\sl inviscid\/} assumption of
\S1.5 \& \S1.6,
and include tangential frictional forces across material surfaces.}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{An example: steady flow past a circular cylinder}

The steady Euler equation with $\bfF=0 $ is satisfied by a {\em potential flow\/} $\bfu = \nabla\phi$  with  $\nabla\cdot\bfu = \nabla^2\phi = 0$,  and pressure from {\em Bernoulli\/} $p + \frac12\rho u^2 = \hbox{\rm const}$.

\medskip
\noindent\begin{minipage}{0.55\linewidth}
The solution with $\phi\sim Ux$ as $r\to\infty$
and 
$\bfu\cdot\bfn = \partial\phi/\partial r = 0$ on $r=a$ is
\[
 \phi = U(r + a^2/r)\cos\theta,
\]
with associated {\em streamfunction\/}
\[
 \psi=U(r-a^2/r)\sin\theta, 
\]
and tangential velocity $2U\sin\theta$ on $r=a$.
\end{minipage}
\begin{minipage}{0.39\linewidth}
  \cpicture{0.9}{streamline} 
\end{minipage}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Books for Part II}

\begin{itemize}
 \item D.J. Acheson {\sl Elementary Fluid Dynamics.\/} Oxford University Press -- {\sl elementary}
  \item G.K. Batchelor {\sl An Introduction to Fluid Dynamics\/}. Cambridge University Press. -- {\sl heavy style, for lecturers} 
 \item E. Guyon, J-P. Hulin, L. Petit \& C.D. Mitescu {\sl Physical Hydrodynamics.\/} Oxford University Press -- {\sl physicists' insights}
 \item Homsy, G.M. et al. {\sl Multi-media Fluid Mechanics.\/} Cambridge University Press 2000 -- {\sl lots of video clips, some beyond Part II}
 \item M. Van Dyke {\sl An Album of Fluid Motion.\/} -- {\sl good B\&W photos, now collectors item}
\end{itemize}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\chapter{Equations of motion for a Newtonian viscous fluid}
\label{chpt:NS}


%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Viscosity}

\cpicture{0.7}{Shearing}

\medskip
\noindent In viscous fluids there is a tangential component of the surface forces, e.g.~when a plate slides past another, there is a drag force $D$, proportional to the surface area $A$, velocity difference $V$, and inversely proportional to separation $h$
\[
 D = \mu \frac{AV}{h},
\]
with coefficient $\mu$, the {\em viscosity\/}, with units $ML^{-1}T^{-1}$; for water it is $1.1\,10^{-3}$, air $1.8\,10^{-5}$, and honey $10^{-1}$ $\rm kg\,m^{-1}s^{-1}$.

\medskip
How can we generalise this for one plate sliding over another to three-dimensional flows?
The velocity difference between the plates will become the  strain-rate, while the tangential force will become the stress. Both the strain-rate and the stress are  2nd order tensors. 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Rate of strain tensor}
\label{sec:strainrate}

Consider the velocity near to a fixed point, wlog ${\bf 0}$,
\[
 u_i(\bfx) = u_i({\bf 0}) + x_j\left.\pd{u_i}{x_j}\right\vert_{\bf 0} +
  \tfrac12 x_jx_k\left.\pd{^2u_i}{x_jx_k}\right\vert_{\bf 0} + \cdots
\]
We keep just the first two terms on the right hand side. Define the {\em symmetric\/} and the {\em antisymmetric\/} parts of the velocity gradient
\begin{align*}
  e_{ij} = \tfrac12\left(\pd{u_i}{x_j} + \pd{u_j}{x_i}\right)
   & \quad \hbox{\rm -- strain-rate tensor}\\
  \Omega_{ij} = \tfrac12\left(\pd{u_i}{x_j} - \pd{u_j}{x_i}\right)
   & \quad \hbox{\rm -- vorticity tensor}
\end{align*}

\noindent In 3D, there are three independent components of an antisymmetric 2nd order tensor
\[
 \Omega_{ij} =
  \begin{pmatrix}
   0 & \Omega_3 & -\Omega_2 \\
   -\Omega_3 & 0 & \Omega_1 \\
   \Omega_2 & -\Omega_1 & 0
   \end{pmatrix}
  = \epsilon_{ijk}\Omega_k
\]
Note the IB {\em vorticity vector\/}
   $\boldsymbol\omega = \nabla\times\bfu = \epsilon_{ijk}\pd{u_j}{x_k} = -2\boldsymbol\Omega$

\noindent So
\[
 u_i(\bfx) = u_i({\bf 0}) + e_{ij}x_j + \tfrac12(\boldsymbol\omega\times\bfx)_i + \cdots
 \]
 The last term describes a local solid-body rotation.

   Now $e_{ij}$ is a symmetric 2nd order tensor, so has real and orthogonal eigenvectors. Choose coordinate axes parallel to these evectors, so that $e$ is diagonal
\[
 e = \begin{pmatrix} e_1 & 0 & 0 \\ 0 & e_2 & 0 \\ 0 & 0 & e_3
     \end{pmatrix}.
\]
But incompressible is
\[
 0 = \nabla\cdot\bfu = e_1 + e_2 + e_3
\]
Wlog we take $e_1>0$ and $e_2<0$, so that the straining part of the flow is

\medskip
\begin{minipage}{0.5\linewidth}
\cpicture{0.6}{straining}
\end{minipage}
\begin{minipage}{0.5\linewidth}
 \[ \bfu= (e_1x_1, e_2x_2, e_3x_3) \]
\end{minipage}
\bigbreak

\noindent Example simple shear $\bfu = (\gamma y,0,0)$ breaks down into
  
\[
 \nabla u = \begin{pmatrix} 0 & 0  & 0 \\ \gamma&0 &0\\0&0&0 \end{pmatrix}
\quad  = \quad\underset{\rm strain}{\begin{pmatrix} 0 & \frac\gamma2 & 0 \\ \frac\gamma2 &0&0\\0&0&0 \end{pmatrix}}
\quad   + \quad\underset{\rm vorticity}{\begin{pmatrix} 0 & - \frac\gamma2 & 0 \\ \frac\gamma2 &0&0\\0&0&0 \end{pmatrix}}
\]

\cpicture{0.7}{SimpleShear}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Stress tensor}


%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Tractions}

In continuum mechanics, deformable solids or fluids, we smear out molecular details, producing two types of forces (derived forces):
\begin{itemize}
 \item {\em volume/body forces\/}, like gravity, that act on a volume, $\bfF(\bfx,t)\delta V$ for a small volume $\delta V$, and
 \item {\em surface tractions\/}, like pressure and tangential viscous forces, that act across a surface, $\boldsymbol\tau(\bfx,t;\bfn)\delta S$ for a small surface $\bfn\delta S$. The force is exerted by side into which $\bfn$ point. The force comes from short range ($\sim$ mean-free-path) activities.
\end{itemize}


%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Tensor}

(We will find that vector $\boldsymbol\tau \delta S$ is linear in vector $\bfn \delta S$.)

\medskip
\noindent
\begin{minipage}{0.45\linewidth}
  Take $\bfn\delta S$ to be the inclined surface of a small tetrahedron, whose other sides coincide with the coordinate planes.
\end{minipage}
\begin{minipage}{0.55\linewidth}
  \cpicture{0.8}{tetra}
\end{minipage}

\bigskip

\noindent Consider the force balance.
The volume and acceleration forces are $O(\rho g L^3)$, where $L$ is the linear size of the tetrahedron.
The surface forces are $O(pL^2)$ with a typical pressure $p = \rho g H$ ($H$ height of the atmosphere).
Hence for a small tetrahedron $L\ll H$, the surface forces are $O(H/L)$ larger, so must balance amongst themselves
\[
 \boldsymbol\tau(\bfn)\delta S
  + \boldsymbol\tau({\bf -1})\delta S_1
  + \boldsymbol\tau({\bf -2})\delta S_2
  + \boldsymbol\tau({\bf -3})\delta S_3 = 0.
\]

\noindent But by the geometry of the tetrahedron
\[
 \delta S_i = n_i\delta S.
\]
Hence the surface traction $\boldsymbol\tau(\bfx,t;\bfn)$ is linear in the direction $\bfn$
\[
 \tau_i(\bfx,t;\bfn) = \sigma_{ij}n_j,
\]
where the 2nd order {\em stress tensor\/} $\sigma_{ij}= \tau_i({\bf j})$, the $i$th component of force exerted across unit area in the $j$th direction by the $+$ side.
  
\medskip
\noindent Note we have used Newton III: $\boldsymbol\tau({\bf -1})\delta S_1= -\boldsymbol\tau({\bf 1})\delta S_1$, etc.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Symmetry}

(We shall find $\sigma_{ij} = \sigma_{ji}$.)

\medskip
\noindent
Consider the angular momentum balance on a small volume $V$, taking moments about a fixed point in $V$.
The moment of the surface forces is
\[
 \int_S \epsilon_{ijk} x_j \sigma_{kl} n_l \, dS,
\]
which is $O(pL^3)$. But the moment of the volume and acceleration forces is $O(\rho gL^4)$. So again for small $L\ll H$, the surface forces must balance amongst themselves. Then using the divergence theorem,
\[
 0 = \hbox{\rm above}
   = \int_V\pd{}{x_l}\left(\epsilon_{ijk}x_j\sigma_{kl}\right)\,dV
   = \int_V\epsilon_{ijk}\left(\sigma_{kj} + x_j\pd{\sigma_{kl}}{x_l}\right)\,dV.
\]
The 2nd term in 2nd  integral is again $O(\rho g L^4)$, so can be ignored.
Hence
\[
 \int_V \epsilon_{ijk}\sigma_{kj}\,dV = 0
\]
But $V$ is arbitrary, so the integrand must vanish
\[
 \epsilon_{ijk}\sigma_{kj} = 0 \qword{or} \sigma_{ij} = \sigma_{ji}.
\]

\cpicture{0.5}{couple}

Pictorially, balacing the couples exerted finds the stresses equal.


%%%%%%%%%%%%%%%%%%%%%%%%
\section{Constitutive equation for a Newtonian viscous fluid}

(To relate the stress tensor to the strain-rate tensor.)

\medskip
\noindent `Constitutive' says what material is made of.
A fluid, as opposed to a deformable solid, moves continuously under deforming/stretching forces, but can be a rest under pressure.

In an incompressible ($\nabla\cdot\bfu=0$) viscous fluid, the stress tensor has two parts
\begin{itemize}
 \item an isotropic {\em pressure\/} part, as in IB, which produces a traction normal to the surface, with magnitude independent of the orientation of the surface, and
 \item a frictional term due to neighbouring fluid sliding past, called the {\em deviatoric\/} part,
\end{itemize}

\[
 \sigma_{ij} = -p\delta_{ij} + \sigma_{ij}^{\rm dev}\left(\pd{u_k}{x_l} \right).
\]
Here we measure the sliding of neighbouring fluid by $\nabla\bfu$, the first term in a Taylor series. This needs
$\ell_{\rm microstructure} (\sim 10^{-9}\,{\rm m})\ll\ell_{\rm laboratory}$ {\color{darkgray}-- may not be true for processing paper plup.}

\medskip
\noindent For common viscous fluids, we make two further assumptions
 \begin{itemize}
  \item $\sigma^{\rm dev}$ is {\em linear} and {\em instantaneous} 
    in $\partial u_i/\partial x_j$, which needs
    $t_{\rm micro} (\sim 10^{-12}\,{\rm s})\ll t_{\rm lab}$ {\color{darkgray}-- may not be true for processing flood.} This is expressed by
    \[
     \sigma^{\rm dev} = A_{ijkl}\pd{u_k}{x_l},
    \]
 \item the fluid is {\em isotropic} {\color{darkgray}-- may not be true for paper fibre suspension.} The general 4th order isotropic tensor is
    \[
    A_{ijkl} = \mu'\delta_{ij}\delta_{kl} + \mu''\delta_{ik}\delta_{jl}
       + \mu'''\delta_{il}\delta_{jk}.
    \]
 \end{itemize}
Thus
\[
  \sigma^{\rm dev} = \mu'\delta_{ij}\pd{u_k}{x_k} + \mu''\pd{u_i}{x_j}
     + \mu'''\pd{u_j}{x_i}
\]

\noindent But $\nabla\cdot\bfu=0$ (and $\sigma^{\rm dev}_{ii}=0$), so we may take the $\mu'$ term to be zero.
And we have symmetry $\sigma_{ij}=\sigma_{ji}$, so $\mu'''=\mu''$, say just $\mu$.
Thus
 \[
  \sigma_{ij} = -p\delta_{ij} + \mu\left(\pd{u_i}{x_j} + \pd{u_j}{x_i} \right),
 \]
  or very important
\[
 \boxed{ \sigma_{ij} = -p\delta_{ij} + 2\mu e_{ij} },
\]
with strain-rate tensor $e_{ij}$ from \S2.2.

\medskip
\noindent Note that the stress does not depend on the vorticity, as expected: solid-body rotation does not stress the fluid. One needs stretching for stress.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Momentum equation}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Cauchy equation}

(Same considerations for conservation of mass and of energy)
\bigskip
  
\noindent Use an inertial frame.
The momentum inside a fixed arbitrary (smooth) volume $V$ changes in time due to
\begin{itemize}
  \item[(i)] momentum carried across boundary $S$,
  \item[(ii)] volume forces, and 
  \item[(iii)] surface tractions.
\end{itemize}

\cpicture{0.4}{influx2}
 
\noindent The volume of fluid {\em entering\/} $V$ across $\delta S$ ($\bfn$ {\em outward} normal) in $\delta t$  is \hbox{$- \bfu\delta t\!\cdot\! \bfn \delta S$.}
Hence the 
momentum {\em entering\/} $V$ is $\rho\bfu\,(-\bfu\delta t\!\cdot\! \bfn \delta S)$.
Hence
\[
 \frac{d}{dt}\int_V \rho u_i\,dV
  = \int_S \boldsymbol-\rho u_i u_j n_j\,dS
  + \int_V F_i\,dV + \int_S \sigma_{ij}n_j \,dS .
 \]

 \noindent Now $V$ is fixed in time, so the left hand side is
\[
 \int_V\pd{}{t}\left(\rho u_i \right)\,dV.
 \]
 For the two surface integrals on the right hand side, we use the generalise divergence theorem
 \[
  \int_S T_{ij}n_j\,dS = \int_V \pd{T_{ij}}{x_j} \,dV,
 \]
 for
 \[
  \int_V\pd{}{t}(\rho u_i) \,dV =
  \int_V\left( -\pd{}{x_j}(\rho u_i u_j) + F_i + \pd{\sigma_{ij}}{x_j}\right)\,dV.
\]
But V is arbitrary, so the two integrands must be equal.
Now we can expand the derivatives of two products
\[
 \pd{}{t}(\rho u_i) = \pd{\rho}{t}u_i + \rho\pd{u_i}{t}, \quad
  -\pd{}{x_j}(\rho u_i u_j) = -\rho u_j\pd{u_i}{x_j} - u_i\pd{}{x_j}(\rho u_j).
\]
The first term of the first expression and last term of the second expression cancel by mass conservation
\[
 \pd{\rho}{t} + \pd{}{x_j}(\rho u_j) = 0,
\]
so we finally have Cauchy's momentum equation

\[ 
 \boxed{ \rho\left(\pd{u_i}{t} + u_j\pd{u_i}{x_j}\right)
  = F_i + \pd{\sigma_{ij}}{x_j}.   }
\]
This is for any continuum, a fluid or a deformable solid.
The left hand side is often written with the material derivative $Du_i/Dt$, so
\[
 \rho \frac{Du_i}{Dt} = F_i + \pd{\sigma_{ij}}{x_j}
\]

\noindent Sometimes the two terms from the surface integrals are combined into the
\[
 \hbox{\rm\em momentum flux} = \rho u_i u_j - \sigma_{ij}.
\]

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Navier-Stokes equation}
\label{sec:NSE}

  Substituting the expression for the stress in a Newtonian fluid from \S2.4, and assuming $\mu$ is constant
\[
 \pd{\sigma_{ij}}{x_j} = -\delta_{ij}\pd{p}{x_j}
   + \mu\left(\pd{^2u_i}{x_jx_j} + \pd{^2u_j}{x_ix_j}\right).
\]
The last term is $\partial(\nabla\cdot\bfu)/\partial x_i = 0$. Hence Cauchy's momentum equation becomes the Navier-Stokes equation

\[
 \boxed{ \rho\frac{D\bfu}{Dt} = \bfF + \mu\nabla^2\bfu }.
\]

\noindent The left hand side is from Navier and the right hand side  from Stokes.

\medskip
\noindent{\color{darkgray} In a non-inertial rotating frame,
$\frac{D\bfu}{Dt} \to \frac{D\bfu}{Dt} +2\bfom\times\bfu + \bfom\times(\bfom\times\bfx)$. }

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Dynamic boundary condition}

Mass conservation requires the normal component of the velocity, $\bfu\!\cdot\!\bfn$, to be continuous.
The momentum equation requires $\nabla\bfu$ to be finite onto the boundary, for else $\sigma$ would be infinite. Hence additionally we need the tangential component of the velocity to be continuous
\[
 \bfu - (\bfu\cdot\bfn)\bfn \quad\hbox{\rm continuous}.
\]
This is called the {\em no-slip\/} boundary condition.
Hence the full $\bfu$ is continuous {\color{darkgray}-- may not be true on super-hydrophobic surfaces like lotus leaves.}

If one applies momentum conservation to a small volume near to the surface, one concludes that momentum flux $\rho u_i u_j -\sigma_{ij}n_j$ must be continuous. But $\bfu$ is continuous, so
\[
 \boxed{ \sigma_{ij}n_j \quad{\rm continuous}}.
\]
Note that the not all the components of the stress tensor $\sigma_{ij}$ need be continuous.

\medskip{\color{darkgray}
At small scales, surface tension can be important, which modifies the boundary condition to
\[
 \sigma^+_{ij}n_j - \sigma^-_{ij}n_j
  = n_i\gamma\left(\frac{1}{R_1} + \frac{1}{R_2}\right),
\]  
where $\gamma$ is the surface tension and $R_i$ are the two principal radii of curvature of the surface. At the free surface between water and air, one can often ignore wind, so $\sigma^{\rm water}_{ij}n_j=0$.}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Energy equation}
\label{sec:energy}

We form the {\em mechanical energy equation} for a continuum by contracting $u_i$ with the Cauchy momentum equation
\[
 u_i\rho\left(\pd{u_i}{t} + u_j\pd{u_i}{x_j}  \right)
  = u_iF_i + u_i\pd{\sigma_{ij}}{x_j}.
\]
\begin{align*}
  {\rm LHS \;=\;}&
   \rho\pd{}{t}\left(\tfrac12u_iu_i\right)
    + \rho u_j\pd{}{x_j}\left(\tfrac12u_iu_i\right)
    + \overset{\hbox{{\rm by mass} = 0}}{\left(\tfrac12u_iu_i\right)\left(\pd{\rho}{t}
    + \pd{}{x_j}(\rho u_j)\right)} \\
    =\;& \pd{}{t}\left(\tfrac12\rho u_iu_i\right)
    + \pd{}{x_j}\left(u_j\tfrac12\rho u_iu_i\right), \\
  {\rm RHS \;=\;}&
    u_iF_i + \pd{}{x_j}(u_i\sigma_{ij}) - \pd{u_i}{x_j}\sigma_{ij}.
\end{align*}
Now integrate over a fixed volume $V$ and use the Divergence Theorem
\begin{multline*}
  \underset{(1)}{\frac{d}{dt}\int_V \tfrac12\rho u^2\,dV}
    + \underset{(2)}{\int_S \tfrac12\rho u^2(\bfu\!\cdot\!\bfn)\,dS} \\
  = \underset{(3)}{\int_V \bfu\!\cdot\!{\rm F}\,dV}
    + \underset{(4)}{\int_Su_i\sigma_{ij}n_j\,dS}
    - \underset{(5)}{\int_V\pd{u_i}{x_j}\sigma_{ij}\,dV}.
\end{multline*}
\begin{itemize}
\item[(1)] rate of change of kinetic energy in $V$ is due to
\item[(2)] flux of kinetic energy over boundary $S$,
\item[(3)] working against the volume forces inside $V$,
 \item[(4)] working against the surface tractions on $S$, and
 \item[(5)] {\em stress working} in $V$.      
 \end{itemize}

So far this is for a general continuum.
We now specialise to a Newtonian viscous fluid with
\[
 \sigma_{ij} = -p\delta_{ij} + 2\mu e_{ij}.
\]
The stress working term becomes
\[
 \sigma_{ij}\pd{u_i}{x_j} = -p\delta_{ij}\pd{u_i}{x_j} + 2\mu e_{ij}\pd{u_i}{x_j}.
\]
Now
 $-p\delta_{ij}\partial u_i/\partial x_j = - p \partial{u_j}/\partial{x_j} = 0$
 by mass conservation,
and $\partial u_i/\partial x_j = e_{ij} + \Omega_{ij}$, with symmetric $e_{ij}$, antisymmetric $\Omega_{ij}$, and the double contraction $e_{ij}\Omega_{ij}=0$ vanishes, as it does any symmetric 2nd order tensor contracted with any antisymmetric 2nd order tensor. Hence
\[
 \sigma_{ij}\pd{u_i}{x_j} = 2\mu e_{ij}e_{ij} = \Phi,
\]
Thus we find the rate of loss of mechanical energy due to friction is the {\em dissipation\/} $\Phi$ per unit volume of fluid.

In this course we do not discuss other forms of energy, e.g.~potential, heat, electrical,\ldots, which are needed for the full engergy equation, which with some thermodynamics yields an equation governing temperature.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Scalings}
\label{sec:scalings}

The governing Navier-Stokes equation is
\[
 \rho\left(\pd{\bfu}{t} + \bfu\cdot\nabla\bfu\right)
  = \bfF -\nabla p  +\mu\nabla^2\bfu
\]
We set $\bfF=0$ for this section.

\medskip
\begin{minipage}{0.47\linewidth}
 Consider flow past an obstacle
\end{minipage}
\begin{minipage}{0.47\linewidth}
  \cpicture{0.6}{scaling}
\end{minipage}

\medskip
\noindent Let the flow have a typical velocity $U$, vary over a typical length $L$, and vary over a typical time $T$. Then the flow depends
\[
 \bfu = \bfu(\bfx,t; \rho, \mu, U, L,T).
\]

Now to nondimensionalise.
  Define the nondimensional velocity $\bfu^* = \bfu/U$, position $\bfx^*=\bfx/L$, time $t^* = t/T$, and pressure (using a viscous scaling) $p^*=p\, L/\mu U$.
Alternatively one could have used an inertial scaling $p^{**}= p/\rho U^2$.
Substitute into the Navier-Stokes equation
\[
 \rho\left(\frac{U}{T}\pd{\bfu^*}{t^*}
   + \frac{U^2}{L}\bfu^*\cdot\nabla^*\bfu^*\right)
 = \frac{\mu U}{L^2}\left(-\nabla^* p^*  +\mu\nabla^{*2}\bfu^* \right),
\]
where $\nabla^* = \partial/\partial x^*$. Dividing by $\mu U/L^2$ we have
\[
 Re\left(St \pd{\bfu^*}{t^*} + \bfu^*\cdot\nabla^*\bfu^*\right)
  = -\nabla^* p^*  +\mu\nabla^{*2}\bfu^*  ,
\]
with two non-dimensional groups
\begin{align*}
  \hbox{\rm Reynolds number } Re\quad=\quad& {\rho UL}/{\mu}, \\
  \hbox{\rm Strouhal number } Sr\quad=\quad& {L}/{UT}.
\end{align*}
\begin{itemize}
 \item If $Sr\ll1$, then the flow is quasi-steady, and one can often ignore $\partial \bfu^*/\partial t^*$.
 \item If $Sr\gg1$, then the flow is rapidly oscillating, and one can often ignore $\bfu^*\cdot\nabla^*\bfu*$.
 \item If $Re\ll1$, then the flow is viscously dominated, and one can often ignore the left hand side, see chapter \ref{chpt:stokes}.
 \item If $Re\gg1$, then the flow is inertially dominated, and one can often ignore $\nabla^{*2}\bfu^*$ in most of the flow, but not everywhere, a dangerous approxamation, see chapter \ref{chpt:bl}.
\end{itemize}
Sometimes we use a Reynolds number for oscillating flows
\[
 Re_T = Re\,Sr = \frac{\rho L^2}{\mu T}.
\]
The ratio $\mu/\rho$ occurs often, so given its own name and symbol: the {\em kinematic viscosity\/}
\[
  \nu = \mu/\rho,
\]
with units $L^2T^{-1}$. We now call the original $\mu$ the {\em dynamic viscosity\/}.
$\mu$ is the viscosity per unit mass, while $\nu$ is the viscosity per unit volume.
Typical values are $\nu = 1.5\,10^{-5}\,{\rm m^2s^{-1}}$ for dry air, and
  $1.1\,10^{-6}\,{\rm m^2s^{-1}}$ water, equivalently $\tfrac34$ acre/year. Note that  air is more viscous than water in this measure.

\medskip
We now write the Reynolds number
\[
 Re = \frac{UL}{\nu}
\]
Values vary enormously:
\begin{itemize}
 \item E.g.~a low speed aircraft: $U=100\, {\rm m/s}$, $L=2\,{\rm m}$,
$\nu=1.5\,10^{-5}\,{\rm m^2/s}$, gives $Re=1.3\,10^7$.
  \item E.g.~a water droplet in a cloud  $U=1\,\mu{\rm m/s}$, $L= \frac12\,\mu{\rm m}$, $\nu =$ same, gives $Re=3\,10^{-8}$.
\end{itemize}

Two flows $\bfu(\bfx,t; \rho,\mu, U, L, T)$ with the same Reynolds number(s) $Re$ (and $Re_T$ if unsteady), with the same shape geometry, e.g.~sphere, (but possibly different $L$)
will have the same nondimensional flow $\bfu^*(\bfx^*,t^*; Re, Re_T)$ and are called {\em dynamically similar\/}.
This is basis of scale model testing; e.g.~wind tunnel testing of scaled aircraft and cars.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Governing equations for this course}

Mass conservation for incompressible fluid
\[
  \nabla\cdot\bfu = 0.
\]

\noindent Momentum conservation
\[
 \rho\left(\pd{\bfu}{t} + \bfu\cdot\nabla\bfu  \right) =
  \bfF - \nabla p + \mu\nabla^2\bfu.
  \]

\noindent Boundary conditions
\[
 \bfu \quad\hbox{{\rm continuous}},
  \qquad \sigma_{ij}n_j \quad\hbox{{\rm continuous}},
\]
with $\sigma_{ij} = -p\delta_{ij} + 2\mu e_{ij}$,
\quad $\displaystyle 2e_{ij} = \pd{u_i}{x_j} + \pd{u_j}{x_i}$.


\begin{ex}[Simple shear]\normalfont
Show that steady simple shear flow $\bfu = (\gamma y,0,0)$ is the
sum of a planar extensional flow (whose principal axes should be
determined) and a solid-body rotation.
Show that the Navier-Stokes equations are satisfied if the pressure is
constant and the body force vanishes.
If the flow is maintained between two plates at $y=0$ and $y=h$, find
the forces on the plates.
\end{ex}

\begin{ex}[Elliptic and hyberbolic flows]\normalfont
Consider the two-dimen\-sional linear flow
\[
 \bfu = (\alpha x - \tfrac12\omega y, -\alpha y + \tfrac12\omega x) .
\]
Confirm that this flow is incompressible and find its streamfunction.
Show that the streamlines are elliptic or hyperbolic according to
whether $\vert\alpha\vert \lessgtr \tfrac12\vert \omega \vert$.

Evaluate $\rho\bfu\cdot\nabla\bfu$ and find a pressure field to balance
it. 
Discuss the minimal or maximal nature of the pressure at the origin in
terms of the streamline pattern.

\end{ex}

\begin{ex}[Dissipation]\normalfont
Show for  a volume $V$ with a {\sl stationary\/} rigid
boundary that the total rate of dissipation of energy can be written
alternatively as 
\[
 2\mu\int_V e_{ij}e_{ij} \,dV = \mu \int_V \omega^2 \,dV,  
 \qquad{\rm where}\quad \omega=\vert\nabla\times\bfu\vert.
\]
It follows that if the flow is irrotational, there is no dissipation. Why?
\end{ex}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\chapter{Unidirectional flows}
\label{chpt:uni}

This chapter examines in finer detail flows which were presented in the earlier IB course.
We consider flows in one direction, say the $x$-direction, with no variation in the direction of flow $\partial/\partial x \equiv 0$, with possible variations in time and in directions perpendicular to the flow. Thus
\[
 \bfu(\bfx,t) = (u(y,z,t),0,0).
\] 
Then automatically satisfy mass conservation
\[
 \nabla\cdot\bfu = \pd{u}{x} \equiv 0.
\]
The left hand side of conservation of momentum reduces to
\[
 \rho\left(\pd{\bfu}{t} + \bfu\cdot\nabla \bfu  \right)
  = \left(\rho\pd{u}{t},0,0  \right),
 \]
 because there is no variation in the in the direction of the flow
 ($\bfu\cdot\nabla \equiv 0$),
 while right hand side becomes
 \[
  -\nabla p + \bfF + \mu\nabla^2\bfu =
   \left(-\pd{p}{x} + F_x + \mu\nabla^2 u,-\pd{p}{y} + F_y,-\pd{p}{z} + F_z \right).
\]
Then if $\bfF = 0$, $p$ is independent of $y$ and $z$, so $\partial p/\partial x$ is independent of $y$ and $z$, so can be written $dp/dx$. Thus we need to solve
\[
 \rho\pd{u}{t} = -\frac{dp}{dx} + \mu\left(\pd{^2u}{y^2} + \pd{^2u}{z^2}\right). 
\]
And because $u$ is independent of $x$, then $dp/dx$ is independent of $x$, so it is a constant.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Plane Couette flow}

(Maurice Couette c1888).

\cpicture{0.7}{Couette}

\noindent Take (i) independent of $z$, (ii) steady $\partial/\partial t \equiv 0$, (iii) $\bfF=0$, and (iv) $dp/dx=0$ because open to atmospheric pressure at both ends. Hence governing equation becomes
\[
 0 = 0 + \mu \frac{d^2u}{dy^2} \qword{in} 0\le y \le h.
\]
The boundary conditions are
\[
 u=V \qword{on} y=h, \quad u=0\qword{on} y=0.
 \]
 The solution is
 \[
 u = \gamma y \qword{with} \gamma = \frac{V}{h},
 \]
with strain-rate
\[
e_{ij} = \tfrac12\left(\pd{u_i}{x_j}+\pd{u_j}{x_i} \right) =
 \begin{pmatrix}0 & \gamma/2, & 0\\ \gamma/2 & 0 & 0\\ 0&0&0 \end{pmatrix},
\]
and stress
\[
\sigma_{ij} = -p\delta_{ij} + 2\mu e_{ij} =
  \begin{pmatrix}-p & \mu\gamma & 0\\ \mu\gamma & -p & 0\\ 0&0&-p \end{pmatrix}.
\]
Thus the force per unit area exerted by the top plate on the fluid, $y=h$, $\bfn=(0,1,0)$
\[
 \sigma_{ij}n_j = (\mu\gamma,-p,0),
\]
and the force per unit area exerted by the bottom plate on the fluid, $y=0$, $\bfn=(0,-1,0)$
\[
 \sigma_{ij}n_j = (-\mu\gamma,p,0).
\]
Note that the forces are equal and opposite because there is no source of momentum between the plates.

\noindent The rate of working by plates per unit area is
\[
 \hbox{\rm force $\times$ velocity} =(\mu\gamma,-p,0)\cdot(V,0,0) + )-\mu\gamma,p,0)\cdot(0,0,0) = \mu\gamma V.
\]
Now the dissipation of kinetic energy due to friction (\S\ref{sec:energy}) is
\[
 \Phi = 2\mu e_{ij}e_{ij} = 2\mu\left(\left(\frac{\gamma}{2}\right)^2+\left(\frac{\gamma}{2}\right)^2   \right) = \mu \gamma^2.
\]
This dissipation per unit area integrated through depth $h$ is
\[
 \int_{y=0}^{y=h} 2\mu e_{ij}e_{ij}\,dy = \mu \gamma^2 h.
\]
Note that it is equal to the rate of working of plates $\mu\gamma V$, because there is no other source of energy.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Plane Poiseuille flow}

Poiseuille 1838, Hagen 1839.

\cpicture{0.8}{Poiseuille}

\noindent Take (i) independent of $z$, (ii) steady $\partial/\partial t \equiv 0$, and (iii) $\bfF=0$. Flow driven by pressure gradient
\[
 \frac{dp}{dx} = \hbox{\rm constant, independent $x,y,z,t$}
  = - \frac{\Delta p}{L}.
\]
Hence the governing equation is
\[
 0 = -\frac{dp}{dx} + \mu\pd{^2u}{y^2} \qword{in} -h\le y\le h.
\]
The boundary conditions are
\[
 u = 0 \qword{on} y = \pm h.
 \]
 The solution is
 \[
  u = -\frac{1}{2\mu} \frac{dp}{dx}(h^2 - y^2).
 \]
The volume flux (per unit length in $z$-direction) is
 \[
 Q = \int_{-h}^h u\,dy = -\frac{2h^3}{3\mu}\frac{dp}{dx}.
 \]
The strain-rate is
\[
e_{ij} = \tfrac12\left(\pd{u_i}{x_j}+\pd{u_j}{x_i} \right) =
 \tfrac12 \frac{du}{dy}
 \begin{pmatrix}0 & 1 & 0\\ \ 1 & 0 & 0\\ 0&0&0 \end{pmatrix} =
  \frac{y}{2\mu} \frac{dp}{dx}\begin{pmatrix}0 & 1 & 0\\ \ 1 & 0 & 0\\ 0&0&0 \end{pmatrix},
\]
and the stress
\[
\sigma_{ij} = -p\delta_{ij} + 2\mu e_{ij} =
\begin{pmatrix}-p & y\tfrac{dp}{dx} & 0\\ y\tfrac{dp}{dx} & -p & 0\\ 0&0&-p \end{pmatrix}.
\]
The traction exerted by the top and bottom boundaries $y=\pm h$, \hfill \linebreak$\bfn = (0,\pm 1,0)$ is
\[
 \sigma_{ij}n_j = (h\frac{dp}{dx}, \mp p,0).
\]
Then total force (per unit length in  $z$-direction) over a $x$ length of $L$ is
\[
 (2h\frac{dp}{dx} L, 0, 0) = (-2h\Delta p,0,0),
\]
whereas net pressure force on the ends is $(2h\Delta p,0,0)$, i.e.~equal and opposite because there is no source of momentum in interior.

\noindent The integrated dissipation (per unit length in $z$-direction) is
\[
 \int_{x=0}^{x=L}\int_{y=-h}^{y=h}
   \left(2\mu e_{ij}e_{ij} = \mu\left(\frac{du}{dy}\right)^2 \right) \,dydx
 \quad= \frac{2h^3}{3\mu}\left(\frac{dp}{dx}\right)^2 L
 = Q\Delta p
\]
which is the rate of working by the flux $Q$ against the pressure drop $\Delta p$, because there is no other source of energy.

\paragraph{\bf Extensions}

\begin{itemize}
 \item Circular cross-section pipe, and elliptical cross-section pipe.
 \item Removing the coating of liquid inside a tube. Take steady, $F=0$, $dp/dx$ given constant (same in both air and water!).
  \cpicture{0.6}{Coating}
  On interface of water and air $u$ and $\mu\partial u/\partial r$ continuous.
\item Flow of two layers of liquid down an inclined plane. 
  \cpicture{0.5}{incline}  
  Take $\bfF= \rho (g\cos\alpha, g\sin\alpha, 0)$, $\rho$ and $\mu$ different in two liquids, $u=0$ on plane, $u$ and $\mu du/dy$ continuous on interface, and $\mu\, du/dy = 0$ and $p=0$ on top surface.
\end{itemize}

\medskip{\color{darkgray}
There are interesting questions of the stability of these steady laminar flows. Reynolds found by experiments in 1883 that pipe flow was unstable if $Re \gtrsim 1000$ when the flow became time-dependent and three-dimensional.}
    

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Impulsively started flat plate}
\label{sec:imp}

Rayleigh (1911) problem, but Stokes (1851)

\cpicture{0.5}{Rayleigh}
Take (i) $\bfF=0$, (ii) $dp/dx=0$, and (iii) initially at rest
\[
 u = 0, \qword{at} t=0, \qword{in} y>0.
\]
The plate $y=0$ starts at $t=0$ to move at a constant velocity $U$ in its own plane
\[
 u = U \qword{at} y=0, \qword{in} t>0.
\]
The momentum equation in the $x$-direction is
\[
 \rho\pd{u}{t} = \mu\pd{^2u}{y^2} \qword{for} t>0,\quad y>0.
\]

Before solving, note that we have a {\em diffusion equation} for $u$ -- information that the plate has started to move diffuses into the fluid with a diffusivity
\[
 \nu = \mu/\rho,
\]
the {\em kinematic viscosity\/} of \S\ref{sec:scalings}.
The solution can be found readily by Laplace transforms
\[
u(y,t) = U \erfc\left(\frac{u}{\sqrt{4\nu t}} \right),
\]
with error function complement
\[
 \erfc(z) = \frac{2}{\sqrt{\pi}}\int_z^\infty e^{-s^2}\,ds.
\]
An alternative method to derive the solution is to seek a ``similarity'' solution, which is possible because there is no specified length scale. First note that $u(y,t)$ must be linear in $U$. Second, on dimensional grounds $\rho, t, \mu, y$ can only occur in the non-dimensional group $\eta = y/\sqrt{4\nu t}$. Hence seek a solution of the form
\[
 u(y,t) = U f(\eta).
\]
Substitute into the diffusion equation
\[
 - \rho U f' \frac{\eta}{2t} = \mu U \frac{1}{4\nu t} f'',
\]
i.e.
\[
 f'' + 2\eta f' = 0, \qword{so} f'= - A e^{-\eta^2}, \qword{so the above solution.}
\]
Note that the region of the fluid effected $u \gtrsim \frac15 U$ is $y\lesssim \sqrt{\nu t}$. This distance $\sqrt{\nu t}$ through which momentum diffuses in time $t$ will recur many time in the course.
Note that the stress exerted by the fluid on the lower boundary is
\[
\mu\left
.\pd{u}{y}\right\vert_{y=0}
   = - \frac1{\sqrt\pi}\frac{\mu U}{\sqrt{\nu t}},
\]
where the second fraction easily comes from from dimensional analysis, while the first fraction is all that the mathematics contributes.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Oscillating flat plate}
\label{sec:osc}

Stokes

\noindent Look for steady oscillations after transients have decayed. As before
\begin{gather*}
  \pd{u}{t} = \nu \pd{^2u}{y^2} \qword{in} y>0 \\
  u = U\cos\omega t \qword{on} y=0 \\
  u \to 0 \qword{as} y \to \infty
\end{gather*}
Solution
\[
 u = Re\left\{U e^{i\omega t - \sqrt{\frac{i\omega}{\nu}}y}   \right\}
\]
Take $\omega>0$ ($\nu>0$), need $\sqrt{i} = (1+i)/\sqrt2$ for decay at $\infty$. So
\[
u(y,t) = U e^{-\kappa y} \cos(\omega t - \kappa y)
 \qword{where} \kappa = \sqrt{\frac{\omega}{2\nu}}
\]

\cpicture{0.7}{oscillating} 

The figure shows a 'propagation': the plate drags the fluid, the fluid above starts to move and through interia continues to move after the plate has stopped. Note activity in the second wavelength is  very very small,  $e^{-2\pi} \approx 10^{-3}$.

The depth of influence is $\delta = \sqrt{2\nu/\omega}$, i.e.~$\sqrt{\nu t}$ of \S\ref{sec:imp} with $t=\pi/\omega$. In water $\nu= 1.1\,10^{-6}\,{\rm m^2s^{-1}}$, hence $1\,{\rm Hz}$ gives $\delta=0.6\,{\rm mm}$, $1\,{\rm c/day}$ gives $0.17\,{\rm m}$ and $1\,{\rm c/yr}$ gives $3\,{\rm m}$.

The total rate of dissipation within the fluid (per unit area of the plate) is
\[
 \int_0^\infty \mu\left(\pd{u}{y}\right)^2\,dy
   \qword{by} \Phi = 2\mu e_{ij}e_{ij} = \mu (\partial u/\partial y)^2.
\]
This has a time-avergaed
\[
 \tfrac12 \mu U^2 \sqrt{\frac{\omega}{2\nu}}
\]
which equals the time-averaged rate of working against the surface stress
\[
 \mu U \sqrt{\frac{\omega}{\nu}} \cos(\omega t + \tfrac\pi4)
\]
{\color{darkgray} Need this for the slowing of pendulum - Stokes.}
Note the rate of working against the surface stress equals is not equal instantaneously to the rate of dissipation, because of an oscillating amount of kinetic energy.

\section{Flow between two rotation cylnders}

Taylor-Couette flow or circular Couette flow, Couette 1890.

\cpicture{0.5}{CircCouette}
 
\noindent The flow is not quite unidirection:
\[
 \bfu\cdot\nabla\bfu = - \frac{V^2}{r} \hat{\bf r},
 \]
 so centrifugal acceleration creates a pressure $p(r)$ which is low at middle.

The solution can be found by various methods (be warned of problems in curvilinear coordinates where $(\nabla^2\bfu)_\theta \neq \nabla^2u_\theta$) is
\[
 V = Ar + B\frac1r
\]
The first term is solid body rotation. The second term is the line vortex of IB, so a potential flow.
Fixing the constants so that $V=\Omega_1r_1$ at $r=r_1$ and $V=\Omega_2r_2$ at $r=r_2$
\[
 V = \frac{\Omega_1r_1^2-\Omega_2r_2^2}{r_1^2-r_2^2}\,r
      - \frac{(\Omega_1-\Omega_2)r_1^2r_2^2}{r_1^2-r_2^2}\,\frac1r
\]

{\color{darkgray}
To obtain the solution one can use the 
Navier-Stokes equations in cylindrical coodinates, which for the special flow
$\bfu = (0,V(r),0$ with pressure $p(r)$ are 
\begin{align*}
 r:  -\rho\frac{V^2}{r} =& - \pd{p}{r}, \\
 \theta:\qquad 0 =&
    \mu\left(\frac1r\pd{}{r}\left(r\pd{V}{r}\right) -\frac{V}{r^2}\right), \\
 \sigma_{r\theta} =& \mu r\pd{}{r}\left(\frac{V}{r} \right).
\end{align*}
}%color

\bigskip
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{ex}[Elliptical pipe]\normalfont
Fluid flows steadily through a cylindrical tube parallel to the
$z$-axis with velocity $  \bfu = (0,0,w(x,y)) $,
under a uniform pressure gradient $G=-dp/dz$.
Show that the Navier-Stokes equations with no body force are satisfied
provided
\[
 \nabla^2 w = -G/\mu,
\]
and state the appropriate boundary conditions.

For a tube with an elliptical cross-section with semi-axes
$a$ and $b$, show that
\[
 w = w_0\left(1-\frac{x^2}{a^2} - \frac{y^2}{b^2} \right) ,
\]
finding $w_0$.
Show that the volume flux (i.e.~the volume of fluid
passing any section of the tube per unit time) is given by
\[
 Q = \frac{\pi a^3 b^3 G}{4(a^2 + b^2) \mu} .
\]

Now specialise to circular cross-section, $b=a$.
Show that the viscous stress on the boundary, which you may take to
be $\sigma_{rz} = \mu \partial w/\partial r$, produces an axial force 
$4\pi\mu w_0 L$ on a length $L$ of the tube, and that this balances
the pressure difference exerted across the ends $LG\pi a^2$.
Further show that the dissipation within the tube is 
$2\pi\mu w_0^2 L$ and this is equal to the rate of working against the
pressure difference across the ends $LGQ$.
\end{ex}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{ex}[Inclined two-layer]\normalfont
Two incompressible fluids of the same density $\rho$ and
viscosities $\mu_B$ and $\mu_T$ flow steadily, one on top of the
other, down a plane inclined at an angle $\alpha$ to the horizontal.
The depths of the layers (normal to the plane) are uniform and equal
to $h_B$ and $h_T$ respectively.

Using coordinates $x$ down the plane and $y$ perpendicular to it, 
write down the boundary conditions on the plane, on the interface
between the two layers and on the top free surface.
Find the pressure field and velocity field in each fluid
on the assumption that they depend only on $y$.
Observe that the velocity profile in the bottom layer depends on
$h_T$ but not $\mu_T$. Why?
\end{ex}

%%%%%%%%%%%%%%%%%%%%%%%%
\begin{ex}[Dissipation for oscillating plate]\normalfont
A plane rigid boundary of a semi-infinite domain of fluid
oscillates in its own plane with velocity $U\cos\omega t$, and the
fluid is at rest at infinity.
Find the velocity field.
[Hint: use $e^{-\kappa(1+i)z}$ with $\kappa^2=\omega/2\nu$.]
Show that the time-averaged rate of dissipation of energy in the
fluid is
\[
 \tfrac12\rho U^2\left(\tfrac12\nu\omega\right)^{1/2}
\]
per unit area of the boundary.
Verify that this is equal to the time average of the rate of work of
the boundary on the fluid (per unit area).
\end{ex}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{ex}[Couette flow]\normalfont
Viscous fluid is contained in the space between two coaxial
cylinders $r=a$ and $b$~($>a$), which may be consider to be of
infinite length.
The inner cylinder rotates with steady angular velocity $\Omega$
about its axis and the outer cylinder is at rest.
The velocity field in the fluid is steady and of the form
$\bfu = (0,v(r),0)$ in cylindrical polar coordinates,
and the pressure varies only in the radial direction.
Look up the components of the Navier-Stokes equations in these
coordinates, say in Appendix 2 of Batchelor or Wikipedia.
[Alternatively work in Cartesians, with 
$\bfu = (yf(r), -xf(r),0)$ with $r^2 = x^2 + y^2$, using 
$\partial_x f = xf'/r$.]
Show that 
\[
 v(r) = Ar + B/r,
\]
where $A$ and $B$ are to be determined.
Calculate the torque per unit length that must be applied to the
inner cylinder to maintain the motion.
[Use the component $e_{\theta r}$ of the strain-rate tensor in
cylindrical polars, given by 
$2e_{\theta r} = r\, d(v/r)/dr$
in this flow.]
\end{ex}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{ex}[Oscillating channel flow]\normalfont
Fluid having kinematic viscosity $\nu$ and density $\rho$ is
confined between a fixed plate at $y=h$ and a plate at $y=0$ whose
velocity is $(U\cos\omega t,0,0)$, where $U$ is a constant.
There is no body force and the pressure is independent of $x$.
Explain the physical significance of the dimensionless number
$S=\omega h^2/\nu$.

Assuming that the flow remains time-periodic and unidirectional,
find expressions for the flow profile and the time-average rate of
working $\Phi$ per unit area by the plates on the fluid.

Sketch the velocity profile and evaluate $\Phi$ in the limits
$S\to0$ and $S\to\infty$, and explain why in these limits $\Phi$ 
becomes independent of $\omega$ and $h$ respectively.
\end{ex}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\chapter{Stokes flows}
\label{chpt:stokes}

This chapter is concerned with flows with a small Reynolds numbers $Re = UL/\nu \ll1$, sometimes called `creeping flows'.
The Reynolds number can be small for many reasons: 
the velocity $U$ can be small,  e.g.~$1\,{\rm cm/day}$ in oil reservoirs;
the lengthscale $L$ can be small, e.g.~$10\,\mu{\rm m}$ bacteria; 
or the viscosity $\nu$ can be large,  e.g.~$10^8\,{\rm m^2/s}$ molten glass.
An example more relevant to this course is a $1\,\mu{\rm m}$ water droplet falling under gravity in air at $0.1\,{\rm mm/s}$, so $Re=10^{-5}$. 

When the Reynolds number is small, one can ignore the inertial terms of the left hand side of the Navier-Stokes equations.
This produces the Stokes equations, the governing equations for this chapter
\begin{gather*}
 0  = \bfF -\nabla p + \mu\nabla^2\bfu, \\
 \text{with}\quad
    \nabla\cdot\bfu = 0,
\end{gather*}
with boundary conditions of $\bfu(\bfx,t)$ given. One then has to calculate the full flow field $\bfu(\bfx,t)$ and surface stress $\sigma_{ij}n_j$.
Note: Stokes theory for $Re\ll1$ usually works for $Re<2$.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Some simple properties}

Solving the Stokes equations rapidly brings some messy details of algebra,
so first we examine some simple yet very useful propoerties before facing up to the nasty calculations.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Instantaneity}

With no $\partial\bfu/\partial t$ term (no inertia to build up to terminal velocity) the response is instantaneous. As the boundary conditions might change in time, the forces will change in time corresponding to the instantaneous boundary conditions, a so called 'quasi-steady' response.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Linearity}

With no nonlinear $\bfu\cdot\nabla\bfu$ term, the velocity $\bfu(\bfx)$ will be linear in the boundary conditions, and so will be the surface stresses and resultant forces.

For example the drag force $\bfF$ on a particle is linear in its instantaneous velocity $\bfU(t)$,
\[
  \bfF(t) = {\bf A} \bfU(t),
\]
where the second order tensor ${\bf A}$ will depend on the size, shape and orientation of the particle.

This simple result can be applied to an axisymmetric particle to break up the contributions to the force into the linearly additive contributions from the motion parallel to the axis $\bfU_\parallel$ and perpendiculat $\bfU_\perp$,

\begin{minipage}{0.6\linewidth}

\[
 \bfF(t) = \alpha \bfU_\parallel(t) + \beta\bfU_\perp(t),
\]

with two scalars $\alpha$ and $\beta$.
\end{minipage}
\begin{minipage}{0.38\linewidth}
\cpicture{0.45}{linearity}
\end{minipage}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Reversible in time}

Consider a force $\bfF(t)$ applied for a certain time $0\le t\le t_1$.
Now reverse force and its history, i.e. 
$\bfF(t)= -\bfF(2t_1-t)$ in $t_1\le t\le 2t_1$
The flow $\bfu(\bfx,t)$ will reverse and its history reverses.
Hence all fluid particles return to their starting position. 
There as some nice videos demonstarting this.

The reversibility in time means that one cannot swim at $Re\ll1$ by reversible flapping, i.e.~a recipricating motion. Again there are some nice videos of G.I.Taylor demonstration this. To swim at low Reynolds numbers one needs something like a propagating wave or a helical motion.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Reversible in space}

Linearity along with some symmetry of geometry requires certain components of the velocity to vanish.

\begin{itemize}
\item A sphere sedimenting next to a vertical wall does not migrate towards or away from the wall at $Re\ll1$. For suppose not:
  \cpicture{0.9}{wallreverse}
  First and last pictures should be identical, so the velocity component perpendicular to the wall must vanish.
 \item Two equal spheres fall without separating.
 \item An ellipsoid, or any other particle with three perpendicular planes of symmetry, falls under gravity without rotating in an unbounded flow.
 \item Two rigid spheres in a shear flow (possibly unequal, possibly next to a rigid wall) resume their original undisturbed streamlines after a collision.
   \cpicture{0.8}{twospheres}
\end{itemize}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Harmonic functions}

If $\bfF= 0$, then
\begin{align*}
 \nabla\cdot{\rm momentum} \qquad\to&\qquad \nabla^2p=0 \\
 \nabla\times{\rm momentum} \qquad\to&\qquad \nabla^2\omega=0, \quad \omega=\nabla\times\bfu \\
 \nabla^2\,{\rm momentum} \qquad\to&\qquad \nabla^2\nabla^2\bfu=0
   \qword{using} -\nabla \nabla^2p=0\\
\end{align*}
Hence $p$ and $\omega$ are harmonic ($\nabla^2=0$) and $\bfu$ is biharmonic.


%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Flow past a sphere}

We now confront the nasty algebra for the solution of uniform flow $\bfU$ past a fixed rigid sphere of radius $a$.

\subsection{The solution}

The solution is more important than its  derivation
\[
 \bfu  =  \bfU \left(1 - \frac{3a}{4r} - \frac{a^3}{4r^3} \right)
    + \bfx (\bfU\cdot\bfx)\left(-\frac{3a}{4r^3} + \frac{3a^3}{4r^5}\right), 
\]
\[
 p  = - \frac{3a\mu\bfU\cdot\bfx}{2r^3} \qword{and}
  \boldsymbol\sigma\cdot\bfn\big\vert_{r=a}  = \frac{3\mu}{2a}\bfU .
\]
Note that the stress on the surface does not vary with position, so that it is easily integrated to evaluated the Stokes drag on the sphere
\[
\int_{r=a} \boldsymbol\sigma \cdot\bfn \,dS =
 4\pi a^2 \, \frac{3\mu}{2a}\bfU = 6\pi\mu a\bfU .
\]
Note that of the 6, 4 is from pressure and 2 is from the tangential viscous friction.
Note that the drag is $4\pi\mu aU$ for a bubble with boundary conditions $\bfu\cdot\bfn=0$ and $(\sigma_{ij}n_j)_{\rm tangential} = 0$.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Method 1}

The simpest derivation that I know starts by noting that the linearity of the Stokes equations means that $\bfu(\bfx)$ must be linear in $\bfU$.
Further, the problem has spherical symmetry about the centre of the
sphere, which  we take as the origin. 
The velocity and pressure fields must therefore take the forms
\begin{eqnarray*}
 \bfu(\bfx) & = & \bfU f(r) + \bfx (\bfU\cdot\bfx) g(r), \\
 p(\bfx) & = & \mu (\bfU\cdot\bfx) h(r),
\end{eqnarray*}
where $r = |\bfx|$, and $f$, $g$ and $h$ are functions of scalar $r$ to
be determined.

Now differentiating gives
\[
 \pd{u_i}{x_j} = U_ix_j f'/r +  \delta_{ij}U_nx_n g +
   x_iU_j g + x_ix_jU_nx_n g'/r.
\]
Then contracting $i$ with $j$ in this expression, we have the incompressibility condition
\[
 0 = \nabla\cdot\bfu = U_nx_n(f'/r + 4g + rg'),
\]
our first relation between the unknown functions.
Differentiating again for momentum equation
\begin{eqnarray*}
 \mu\nabla^2 u_i & = & \mu U_i\left(f''+ 2f'/r + 2g\right) +
         \mu x_iU_nx_n\left(g'' + 6g'/r\right) \\
 \nabla_ip & = & \mu U_i h + \mu x_iU_nx_n h'/r
\end{eqnarray*}

Hence the governing equations give three relations between the unknown functions
\[
 f'/r + 4g + rg' = 0,\quad
 f''+ 2f'/r + 2g = h \qword{and}
 g'' + 6g'/r = h'/r .
\]
Eliminating $h$ and then $f$ yields
\[
 r^2g''' + 11rg'' + 24 g' = 0.
\]
This differential equation is homogeneous in $r$ had so has solutions of the form $g=r^\alpha$. Substituting one finds $\alpha = 0$, $-3$ and $-5$. 
The associated other two functions are 
$f = -(\alpha+4)r^{\alpha+2}/(\alpha+2)$ and
$h = -(\alpha+5)(\alpha +2) r^\alpha$.


Hence the general solution of the assumed form linear in $\bfU$ is
\begin{align*}
 \bfu(\bfx) & =  
   \bfU \left(-2Ar^2 + B + Cr^{-1} -\tfrac13Dr^{-3} \right) + 
 \bfx (\bfU\cdot\bfx)\left( A + Cr^{-3} + Dr^{-5} \right), \\
 p(\bfx) & =  \mu (\bfU\cdot\bfx) \left(-10A + 2Cr^{-3}  \right).
\end{align*}
We shall need the stress exerted across a spherical surface with
unit normal $\bfn = \bfx/r$
\[
 \boldsymbol\sigma \cdot\bfn = 
  \bfU \left(-3Ar + 2Dr^{-4} \right) + 
 \bfx (\bfU\cdot\bfx)\left( 9Ar^{-1} - 6Cr^{-4} - 6Dr^{-6} \right).
\]

Applying the boundary conditions of no flow on the rigid sphere and uniform flow far from the sphere, we find the coefficients
\[
 A = 0, \qquad B=1, \qquad C = -\tfrac34a \qword{and} D = \tfrac34a^3,
\]
And hence the solution given earlier.

{\color{darkgray}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Method 2}

An alternative derivation uses a Stokes streamfunction for the axisymmetric flow
\[
 u_r = \frac{1}{r^2\sin\theta}\pd\Psi\theta \qword{and}
 u_\theta = - \frac{1}{r\sin\theta}\pd\Psi{r} .
\]
The vorticity equation (curl of the momentum equation, to eliminate
the pressure) is then at low Reynolds numbers (after a lot of algebra)
\[
 {\cal D}^2 {\cal D}^2 \Psi = 0 \qword{where}
 {\cal D}^2 = \pd{^2}{r^2} + \frac{\sin\theta}{ r^2}\pd{}{\theta}
   \left(\frac{1}{\sin\theta}\pd{}{\theta} \right) .
\]
The uniform flow at infinity has $\Psi = \frac12Ur^2\sin^2\theta$, 
so one tries $\Psi = F(r)\sin^2\theta$, 
and finds (after a lot more work) $F = Ar^4 + Br^2 + Cr + D/r$.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Method 3}

One can show (Papkovich-Neuber) that the general solution of the 
Stokes equation can be
expressed in terms of a vector harmonic function $\boldsymbol\phi(\bfx)$
(i.e.~$\nabla^2\boldsymbol\phi = 0$) 
\[
\bfu = 2\boldsymbol\phi - \nabla(\bfx\cdot\boldsymbol\phi) \qword{with} p = - 2\mu\nabla\cdot\boldsymbol\phi,
\]
\[
 \text{and}\quad
 \sigma_{ij} = 2\mu 
  \left(\delta_{ij}\pd{\phi_n}{x_n} - x_k\pd{^2\phi_k}{x_i \partial x_j}\right)
\]
both requiring some detailed calculations.
Linearity and spherical symmetry then give
\[
 \boldsymbol\phi = \alpha\bfU\frac{a}{r} + \beta\bfU\cdot\nabla \nabla \frac{a^3}{r} ,
\]
with coefficients $\alpha$ and $\beta$ to be determined by applying the
boundary conditions, again  a lot of work.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Method 4}

In Stokes flow, the pressure and vorticity are harmonic functions.
Using linearity and spherical symmetry, they must take the form
\[
 p = \mu A \bfU\cdot\bfx/r^3 \qword{and}
 \nabla\wedge\bfu = B \bfU\wedge\bfx/r^3,
\]
withe coefficients $A$ and $B$ to be determined.
The final step to find $\bfu(\bfx)$ is tedious.

}%end color gray

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Sedimentation of a rigid sphere}

Having found the Stokes flow past a sphere, we can use it to find the Stokes settling velocity under gravity of an isolated sphere.
Balancing the forces on the falling sphere, with densities $\rho_s$ of sphere and $\rho_f$ of fluid
\begin{align*}
   0 &\quad =& \rho_s\tfrac{4\pi}3 a^3{\bf g} &\quad- &\rho_f\tfrac{4\pi}3 a^3{\bf g} &\quad -& 6\pi\mu a\bfU \\
 \text{no inertia} & & \text{weight} && \text{buoyancy} & & \text{Stokes drag}
\end{align*}
So we find the Stokes settling velocity
\[
 \bfU_{S} = \frac{2\Delta\rho a^2{\bf g}}{9\mu}
\]
E.g.~$1\,\mu{\rm m}$ sphere, $\Delta\rho=10^3\,{\rm kg\,m^{-3}}$, 
water $\mu = 10^{-3}\,{\rm Pa\,s}$ 
gives $U=2\mu{\rm m/s}$, i.e.~falls through diameter in a second.
Check $Re=10^{-6}$, which is small.

Using the general solution for Stokes flow inside and outside a sphere, one can find (tough student exercise!) that the drag on a fluid sphere is
\[
 \bfF = - 2\pi\frac{2\mu_f +3\mu_s}{\mu_f + \mu_s}\mu_f a \bfU
\]

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Far field}

At large distances, $r\gg a$, the sphere appears just to apply a point force. The disturbance flow is
\[
 \bfu' \sim \bfU\left(-\frac{3a}{4r}\right) 
 + \bfx(\bfU\cdot\bfx)\left(-\frac{3a}{4r^3}\right).
 \]
 Now $\bfF=6\pi\mu a\bfU$, so
 \[
 \bfu' = \frac1{8\pi\mu}\left(\bfF\frac1r
  + \bfx(\bfF\cdot\bfx)\frac1{r^3}\right).
  \]
  The far field depends only on the net force and is independent of the shape of the particle.
 It is the Greens function for Stokes flows, called a {\em Stokeslet\/}.
 The $F/4\pi\mu r$ is the familiar Greens function for a harmonic function, with remainder to ensure an incompressible flow. The form can alternatively be found by Fourier Transforms.
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{A rotating sphere}

The Stokes flow outside a sphere rotating at angular velocity $\boldsymbol{\Omega}$ is simply
\[
 \bfu(\bfx) = \boldsymbol{\Omega}\times\bfx \frac{a^3}{r^3}.
\]
A potential flow, so satisfies Stokes equations.
Hence the couple on the sphere  is (student exercise)
\[
 {\bf G} = - 8\pi\mu a^3 \boldsymbol{\Omega}.
\]

{\color{darkgray}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Stokes flow past an ellipsoid}

For an ellispoid with principle semi-diameters $a_1, a_2, a_3$, Oberbeck (1876) found
\[
 \text{Force}\quad F_1 = - \frac{16\pi\mu U_1}{L + a_1^2K_2},
 \qword{and Couple}
   G_1 = - \frac{16\pi\mu(a_2^2 + a_3^2)}{3(a_2^2K_2 + a_3^2K_3)},
\]
where
\[{
 L = \int_0^\infty \frac{d\lambda}{\Delta(\lambda)} 
  \qword{and}
 K_i =  \int_0^\infty \frac{d\lambda}{(a_i^2+\lambda)\Delta(\lambda)},
}\]
with $\Delta^2 =(a_1^2+\lambda)(a_3^2+\lambda)(a_3^2+\lambda)$.

For a disk $a_1 \ll a_2=a_3$
\[
 F_1 \sim 16\pi\mu a_2U_1, \quad F_2 \sim \tfrac{32}{3}\mu a_2 U_2,
 \quad G_i \sim \tfrac83\mu a_2^3\Omega_i,
\]
while for a rod $a_1 \gg a_2=a_3$, 
\[
 F_1 \sim \frac{4\pi\mu a_1 U_1}{\ln - \tfrac12}, \quad
 F_2 \sim \frac{8\pi\mu a_1 U_2}{\ln + \tfrac12}, \quad
 G_1 \sim \tfrac{16}{3}\pi\mu a_1a_2^2\Omega_1, \quad
 G_2 \sim \frac{\tfrac83\pi\mu a_1^3\Omega_2}{\ln - \tfrac12},
\]
where $\ln = \ln\tfrac{2a_1}{a_2}$.

The important conclusion is that in Stokes flow the drag depends on the largest linear dimension and is otherwise rather insensitive to shape of the particle. 
In translation the drag is approximately (within a factor of 4) $6\pi\mu aU$ with $a$ the largest diameter.

}%color darkgray


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{More properties of Stokes flows}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{A useful result}
\label{sec:useful}

Let $\bfu^S(\bfx)$ and $p^S(\bfx)$ be a Stokes flow with force $\bfF=0$ in $V$, with its $e^S_{ij}$ and $\sigma^S_{ij}$. Let $\bfu(\bfx)$ be any other incompressible flow. Then claim
\[
 \int_V 2\mu e^S_{ij}e_{ij}\,dV = \int_S \sigma^S_{ij}n_ju_i\, dA.
\]

\noindent Now
\[
 2\mu e^S_{ij} = \sigma^S_{ij} + p^S\delta_{ij} \qword{and} 
 p^S\delta_{ij}e_{ij} = p^S\nabla\cdot\bfu = 0,
 \qword{so} 2\mu e^S_{ij}e_{ij} = \sigma^S_{ij}e_{ij}.
\]
And
\[
 \sigma^S_{ij} = \sigma^S_{ji}, \qword{so}
  \sigma^S_{ij}e_{ij} = \sigma^S_{ij}\pd{u_i}{x_j}
   = \pd{}{x_j}\left(\sigma^S_{ij} u_i\right) - \pd{\sigma^S_{ij}}{x_j}u_i.
\]
The last term vanishes because $\bfF=0$. Then by the divergence theorem
\[
 \int_V 2\mu e^S_{ij}e_{ij}\,dV = \int_V \sigma^S_{ij}e_{ij}\,dV =
  \int_V \pd{}{x_j}\left(\sigma^S_{ij} u_i\right)\,dV =
  \int_S \sigma^S_{ij}n_ju_i\, dA.
\]

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Minimum dissipation}

Let $\bfu(\bfx)$ and $\bfu^S(\bfx)$ be two incompressible flows in $V$, both satisfying the same boundary condition $\bfu=\bfu^S=\bfU(\bfx)$ given on $S$.
Let $\bfu^S$ also satisfy the Stokes equation with $\bfF=0$ in $V$.


\noindent Then
\begin{gather*}
 \int_V2\mu e_{ij}e_{ij}\,dV = \int_V2\mu e^S_{ij}e^S_{ij}\,dV \\
  + \int_V2\mu(e_{ij}-e^S_{ij})(e_{ij}-e^S_{ij})\,dV
  + \int_V4\mu e^S_{ij}(e_{ij}-e^S_{ij})\,dV.
\end{gather*}
The middle integral on the right hand side is positive.
The last integral is of the form of the useful result \S\ref{sec:useful}
\[
 \int_V4\mu e^S_{ij}(e_{ij}-e^S_{ij})\,dV = 
 \int_S2\sigma^S_{ij}n_j(u_i-u^S_i)\,dA
\]
But on $S$, $u_i-u^S_i = U_i-U_i=0$, so the last integral vanishes. 
Hence
\[
 \int_V2\mu e_{ij}e_{ij}\,dV \ge \int_V2\mu e^S_{ij}e^S_{ij}\,dV,
\]
i.e.~the Stokes flow $\bfu^S(\bfx)$ has the minimum dissipation out of all incompressible flows satisfying the boundary condition

{\bf Warning this is for the same geometry.} One is not comparing two flows in different geometry, and so one cannot use the result to select the geometry which has the  minimum dissipation.

\medskip
As a flow at non-zero Reynolds number is an incompressible flow which is not a Stokes flow, 
the drag will larger at non-zero Reynolds numbers.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Uniqueness}

If $\bfu^1(\bfx)$ and $\bfu^2(\bfx)$ are two Stokes flows in $V$ satisfying the same boundary conditions, then minimum dissipation gives
\[
 \int_V2\mu(e^1_{ij}-e^2_{ij})(e^1_{ij}-e^2_{ij})\,dV = 0
\]
Hence
\[
 e^1_{ij}-e^2_{ij} = 0 \qword{in $V$,}
\]
i.e.~$\bfu^1-\bfu^2$ is strainless, 
i.e.~a solid body translation $+$ rotation, 
i.e.~zero by the boundary conditions. 
Hence
\[
 \bfu^1(\bfx) = \bfu^2(\bfx) \qword{in $V$.}
\]
Hence Stokes flows are unique.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Geometric bounding}

An application of minumum dissipation.
Consider a rigid cube with sides of length $2L$ moving at $\bfU$ through a viscous fluid, resulting in a drag force $\bfF^{\rm cube}$.

Let $\bfu^S(\bfx)$ be Stokes flow outside the cube. Then dissipation
\[
 \int_{\hbox{\rm outside cube}}2\mu e^S_{ij}e^S_{ij}\,dV
  = \int_{\hbox{\rm cube}}\sigma^S_{ij}n_ju^S_i\,dS
 = -\bfU\cdot\bfF^{\rm cube},
\]
the rate of working by the surface forces on the cube.

Now consider a sphere which just contains the cube, and hence has a radius $a=\sqrt3L$, also moving at $\bfU$.
Define a second flow
\[
 \bfu(\bfx) = \begin{cases}
                \text{the Stokes flow for sphere} & \text{outside sphere,} \\
                \bfU & \text{in gap.}
              \end{cases}
\]

\noindent Then for this second flow
\[
  \int_{\hbox{outside cube}}2\mu e_{ij}e_{ij}\,dV
  = \int_{\hbox{outside sphere}}2\mu e_{ij}e_{ij}\,dV,
\]
because $e_{ij}=0$ in gap,
\[
 = \text{rate of working by sphere}
 = 6\pi\mu\sqrt3L \bfU\cdot\bfU
\]

\noindent Hence minimum dissipation bounds drag $\bfF$ on cube
\[
 -\bfF^{\rm cube}\cdot\bfU \le 6\pi\mu\sqrt3L\bfU\cdot\bfU
\]

\noindent Similarly for sphere just contained inside cube
\[
 6\pi\mu L\bfU\cdot\bfU \le -\bfF^{\rm cube}\cdot\bfU
\]

\noindent Student exercise: bound for tetrahedron (not so tight).

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Reciprocal theorem}

Let $\bfu_1(\bfx)$ and $\bfu_2(\bfx)$ be two Stokes flows inside the same volume $V$ with different boundary conditions on $S$. Then by the useful result \S\ref{sec:useful}, extended to include volume forces ${\bf f}_i(\bfx)$
\begin{gather*}
 \int_V \bfu_1\cdot{\bf f}_2\,dV + \int_S \bfu_1\cdot\boldsymbol\sigma_2\cdot\bfn\,dA
 = \int_V2\mu {\bf e}_1:{\bf e}_2\,dV \\
 = \int_V \bfu_2\cdot{\bf f}_1\,dV + \int_S \bfU_2\cdot\boldsymbol\sigma_1\cdot\bfn\,dA
\end{gather*}
i.e. work done by one velocity on the forces of the other is vice versa.
This is the  Greens theorem in any other subject.

\subsection{Symmetry of resistance matrix}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

An application of the reciprocal theorem.
For a rigid particle moving at $\bfU$ with no volume forces, the drag force $\bfF$ is linear in $\bfU$
\[
 \bfF = {\bf A} \bfU
\]
with second order tensor ${\bf A}$. 
Now consider a particle moving at two velocites, first $\bfU^1$ and then $\bfU^2$. The reciprocal theorem gives
\[
 \int u^1_i\sigma^2_{ij}n_j\,dS = \int u^2_i\sigma^1_{ij}n_j\,dS
\]
But by boundary condition
\[
 \int u^1_i\sigma^2_{ij}n_j\,dS = U^1_i\int \sigma^2_{ij}n_j\,dS = U^1_iF^2_i,
\]
and similarly for the right hand side, so
\[
 \bfU^1\cdot \bfF^2 = \bfU^2 \cdot \bfF^1
 \qword{i.e.} \bfU^{T1}{\bf A}\bfU^2 = \bfU^{2T} {\bf A}\bfU^1
\] 
But $\bfU^1$ and $\bfU^2$ are arbitrary. Hence the second order tensor ${\bf A}$ is symmetric.

\medskip
Now consider  a particle moving in a general rigid body motion in fluid at rest a infinity, translating at $\bfU(t)$ and rotating {(about a selected point)} at $\boldsymbol{\Omega}(t)$, with no volume forces. The force $\bfF(t)$ and couple ${\bf G}(t)$ (about the same selected point) are given by linearity
\[
 \begin{pmatrix} \bfF \\ {\bf G} \end{pmatrix} =
\begin{pmatrix} {\bf A} \; {\bf B} \\ {\bf C}\; {\bf D} \end{pmatrix} 
\begin{pmatrix} \bfU \\ \boldsymbol{\Omega} \end{pmatrix}
\]
with second order tensors ${\bf A}$, ${\bf B}$, ${\bf C}$ and ${\bf D}$ depending on the size, shape and orientation of the particle.
The reciprocal theorem then gives for the two rigid body motions
\[
 \bfU_1\cdot\bfF_2 + \boldsymbol{\Omega}_1\cdot{\bf G}_2 =
 \bfU_2\cdot\bfF_1 + \boldsymbol{\Omega}_2\cdot{\bf G}_1
\]
This is true for all $\bfU_1$ etc, so
\[
 {\bf A} = {\bf A}^T,\quad {\bf B} = {\bf C}^T \qword{and} {\bf D} = {\bf D}^T.
\]
The symmetry $ {\bf B} = {\bf C}^T $ means
\[
 \text{force due to rotating} = \text{couple due to translating.}
\]

\medskip
\noindent These symmetries plus the geometric symmetry of a cube give
\[
 {\bf A} \;\&\; {\bf D} \qword{diagonal, and} {\bf B} = {\bf C}^T = 0 
\]
Thus drag on a cube is parallel to velocity, also true for symmetric tetrahedron.
One needs a ``corkscrew'' feature for ${\bf B} \neq 0$.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Flows in a corner}
\label{sec:corner}

There are many flows which are locally two-dimensionally near a corner
\begin{itemize}
 \item Scraping fluid off a rigid surface
  \cpicture{0.7}{Scraping}
 \item Moving three-phase contact of a spreading drop
  \cpicture{0.7}{MCL}
 \item Hinged plates
  \cpicture{0.7}{Hinge}
 \item Source flow
  \cpicture{0.7}{Source}
 \item Flow past a corner
  \cpicture{0.5}{Corner}
 \item Flow past a slot in a wall
  \cpicture{0.35}{Slot}
  This is equivalent to the corner flow with an angle of zero.
\end{itemize}

It is best to tackle all these 2D flows with a streamfunction $\psi$
\[
 u_r = \frac1r\pd\psi\theta, \quad u_\theta = -\pd\psi{r},
\]
so automatically satisfying the incompressibility, $\nabla\cdot\bfu=0$.

The vorticity is $\boldsymbol\omega = \nabla\times\bfu = (0,0,-\nabla^2\psi)$.
Eliminate the pressure in the Stokes equations by taking its curl forms the {\em vorticity equation\/}
\[
 -\nabla^2\omega = \nabla^4\psi = 0,
\]
i.e.~$\psi$ is a biharmonic function.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Source flow}
\label{subsec:source}

(Jeffrey \& Hamel)

If there is a volume flux $Q$ (per unit $z$-direction), then $u_r \propto Q/r$, so  we try
\[
 \psi = Q f(\theta)
\]
so that $u_r=Qf'/r$ and $u_\theta\equiv0$. The biharmonic equations gives
\[
 f^{iv} + 4f'' = 0,
\]
with general solution
\[
 f = A\cos2\theta + B\sin2\theta + C + D\theta.
\]
Looking for solutions in which $u_r$ is symmetric about $\theta=0$ gives $A=C=0$. Applying a no-slip boundary condition $f'=0$ at $\theta=\pm\alpha$  and a given volume flux
\[
 Q = \int_{-\alpha}^\alpha u_r \,rd\theta = Q(f(\alpha)-f(-\alpha)),
\]
we have
\[
 f = \tfrac12 \frac{\sin2\theta - 2\theta\cos2\alpha}{\sin2\alpha - 2\alpha\cos2\alpha}.
\]
Check $f$ becomes parabolic for $\alpha\ll1$. This solution has a problem at $2\alpha=257^\circ$ when the denominator vanishes -- see \S\ref{sec:backtosource}.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Flow past a corner}
\label{subsec:corner}

(Moffatt 1964)

This flow is not forced by the boundary condition on $\theta=\pm\alpha$ where $\psi_\theta$ and $\psi_r$ vansih, but is forced by disturbances at large $r$ outside the local solution.

We look for eigensolutions (solutions without scale)
\[
 \psi(r,\theta) = r^\lambda(\theta).
\]
Then the biharmonic equation gives
\[
 \nabla^4\psi = \nabla^2\left(r^{\lambda-2}(f''+\lambda^2 f)   \right)
   = r^{\lambda-4}\left((D^2+\lambda^2)(D^2 + (\lambda-2)^2)f  \right) = 0.
\]
This has a general solution
\[
 f = A\cos\lambda\theta + B\sin\lambda\theta
   + C\cos(\lambda-2)\theta + D\sin(\lambda-2)\theta.
\]

There are symmetric flows with $A=C=0$ and antisymmetric flows $B=D=0$. Consider the latter. The no-flux boundary condition $f=0$ at $\theta=\pm\alpha$ reduces the solutuon to
\[
 f = \cos(\lambda-2)\alpha\cos\lambda\theta
      - \cos\lambda\alpha\cos(\lambda-2)\theta.
\]
The eigenvalue $\lambda$ is determined by applying the no-slip boundary condition $f'=0$ and $\theta=\alpha$
\[
0 = -\lambda\cos(\lambda-2)\alpha\sin\lambda\alpha +
   (\lambda-2)\cos\lambda\alpha\sin(\lambda-2)\alpha.
\]
For vertex angles $2\alpha<146^\circ$, all roots for $\lambda$ are complex
\[
 \qword{e.g.~for $2\alpha=90^\circ$} \lambda = 2.37\pm i0.56, \quad
  \lambda = 3.44\pm i0.74, \ldots
\]
For $2\alpha>146^\circ$, the smallest root is real. There is a similar soltion for the symmetric flows, with the critical angle of $159^\circ$.

To understand the complex roots, write
\begin{gather*}
 \psi = r^{\lambda_r+i\lambda_i} f_\lambda(\theta) +\hbox{\rm complex conjugate} \\
    r^{\lambda_r}\left[\cos(\lambda_i \ln r)\hbox{\rm Re}(f)
                       -\sin(\lambda_i \ln r)\hbox{\rm Im}(f)\right]
\end{gather*}
Thus at all fixed $\theta$, $\psi$ oscillates indefinitely as $r\to0$, i.e.~we have a sequence of eddies of ever decreasing size.

\cpicture{0.7}{Moffatt}

\noindent The centres satisfy
\[
 \lambda_i\ln r_{n+1} = \lambda_i\ln r_n -\pi
 \qword{i.e.} r_{n+1}/r_n = e^{-\pi/\lambda_i}.
\]
The magnitude of the circulation (magnitude of $\psi$) decreases by
\[
 \left(\frac{r_{n+1}}{r_n}\right)^{\lambda_r} = e^{-\pi\lambda_r/\lambda_i}
 \quad\sim \frac1{2000}\qword{for $90^\circ$}.
\]
Note when $2\alpha<146^\circ$ the fluid avoids stretching by the sequence of roughly solid body rotations.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Back to source flow}
\label{sec:backtosource}

Beyond $2\alpha=257^\circ$ the source flow of \S\ref{subsec:source} with $\psi\propto r^0$ is dominated by one of \S\ref{subsec:corner}'s eigensoltions with $\psi\propto r^\lambda$ ($\lambda=0$ at $257^\circ$ rising to $\lambda=\tfrac12$ at $360^\circ$).

In the source flow with inlet radius $r_1$ and outlet radius $r_2$, it was assumed that the flow in $r_1\ll r\ll r_2$ would be a scaleless {\em similarity solution\/} not depending on either $r_1$ or $r_2$. While such a flow exists, it is only realised (in the real world) if disturbances from $r=r_1$ and from $r=r_2$ decay into the interior $r_1\ll r\ll r_2$. They do not decay from $r=r_1$ when $2\alpha>257^\circ$.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\bigskip
\begin{ex}[Drag on a cube]\normalfont
A force is applied to a cube at its centre in a direction normal to
one flat surface. 
Using reversibility in space, show that the cube moves in the
direction of the applied force, also without rotating.
Now using linearity, deduce that in all orientations a cube of uniform
density sediments vertically without rotating.
[{\sl Hints: resolve force into components, and isotropy.}]

[** What of a tetrahedron, an ellipsoid and a helix? **]
\end{ex}

%%%%%%%%%%%%%%%%
\begin{ex}[Two equal spheres]\normalfont
Show that in Stokes flow two equal spheres arbitrarily aligned fall under gravity at constant separation, i.e.~neither separating nor coming closer together. 

\end{ex}

%%%%%%%%%%%%%%%%
\begin{ex}[Strainless flow]\normalfont
If the strain-rate tensor ${\bf e}(\bfx)$ vanishes throughout a
connected region, show that the flow is rigid body motion.

[{\sl Hint: first show $\partial^2u_1/\partial x_2 \partial x_3 \equiv
0$.}]

Show that if the surface traction is specified on a bounding
surface, then the Stokes flow in the interior is unique to within the
addition of a rigid body motion.
\end{ex}

%%%%%%%%%%%%%%%%
\begin{ex}[Rotating sphere]\normalfont
Derive the Stokes flow outside a rotating rigid sphere 
\[
 \bfu(\bfx) = {\boldsymbol\Omega}\times\bfx\frac{a^3}{r^3} \qword{and} p = 0.
\]
Show that the couple exerted on the sphere is $-8\pi\mu a^3{\boldsymbol\Omega}$.
\end{ex}

%%%%%%%%%%%%%%%%
\begin{ex}[Papkovich-Neuber]\normalfont
If ${\bf A}(\bfx)$ is a vector harmonic function,
i.e.~$\nabla^2{\bf A} = 0$, show that
$$
 \bfu = 2{\bf A} - \nabla({\bf A}\cdot\bfx) \qword{and}
 p = -2\mu\nabla\cdot{\bf A}
$$
satisfy the Stokes equation. Calculate the stress tensor.

For a sphere of radius $a$ translating at velocity ${\bf V}$ through a fluid which is otherwise at rest, the harmonic function takes the form
\[
 {\bf A} = \alpha a {\bf V} \frac{1}{r} + 
   \beta a^3 \left({\bf V}\cdot\nabla\right)\nabla \frac{1}{r},
\]
(Why?) Find the constants $\alpha$ and $\beta$.
\end{ex}

%%%%%%%%%%%%%%%%
\begin{ex}[Spherical bubble]\normalfont
Consider a spherical bubble of radius $a$ in a uniform flow $\bfU$.
Recall the expression obtained in lectures for the Stokes flow outside a sphere of the form
\[
 \bfu(\bfx) = \bfU f(r) + \bfx(\bfU\cdot\bfx)g(r).
\]
Applying boundary conditions on $r=a$ of no normal component of velocity and no tangential component of surface traction, find the flow $\bfu(\bfx)$.
Find the drag force $4\pi \mu a \bfU$.
\end{ex}

%%%%%%%%%%%%%%%%
\begin{ex}[Faxen]\normalfont
When a rigid sphere of radius $a$ translates with velocity ${\bf U}$ through unbounded fluid at rest at infinity, it may be shown that the traction per unit area, $ \boldsymbol\sigma \cdot {\bf n}$, exerted by the sphere on the fluid has the uniform value  $3\mu{\bf U}/a$ over the sphere surface. 
Find the drag on the sphere.

Suppose that the same sphere is free of external forces and is placed with its centre at the origin in an unbounded Stokes flow given in the absence of the sphere as ${\bf u}^*({\bf x})$.
By applying the reciprocal theorem to the perturbation to the flow generated by the presence of the sphere, and assuming this tends to zero sufficiently rapidly at infinity, show that the instantaneous velocity of the centre of the sphere is
\[
 {\bf V} = \frac1{4\pi a^2}\int_{r=a}{\bf u}^*({\bf x})\, dS.
 \]
\end{ex}

%%%%%%%%%%%%%%%%
\begin{ex}[Bounds]\normalfont
Find upper and lower bounds for the couple on a tetrahedron
rotating about its centre in a viscous fluid.
\end{ex}

%%%%%%%%%%%%%%%%
\begin{ex}[Sperical annulus]\normalfont
A spherical annulus of incompressible viscous liquid occupies
the region $R_1(t) < r < R_2(t)$ between two free surfaces on which
pressures (normal traction) $P_1(t)$ and $P_2(t)$ are applied.
The resulting flow is spherically symmetric.
Show (neglecting inertia and surface tension)
\[
 \frac{d}{dt}\left(R_1^3\right) \;=\; \frac{\pi(P_1 - P_2)}{\mu V}
 R_1^3\left(R_1^3+ 3V/4\pi\right) ,
\]
where $V$ is the constant volume of the liquid.
[{\sl Hints: $u_r=A/r^2$ (why?) and $\sigma_{rr} = -p + 2\mu\partial u/\partial r$ in this
flow.}]

Show that if $P_1-P_2$ is maintained positive and constant, then $R_1$
becomes infinite in a finite time. What happens if 
 $P_1-P_2$ is maintained negative and constant.
\end{ex}

%%%%%%%%%%%%%%%%
\begin{ex}[Hinged plates]\normalfont
Fluid is contained in the region $-\alpha < \theta < \alpha$
between two rigid hinged plates. 
Thus the velocity components in plane
polar coordinates satisfy
\[
u_r = 0, \qquad u_\theta = \mp \omega r \qword{on} \theta = \pm \alpha .
\]
Neglecting inertia forces, show that {\sl a} solution to the Stokes problem
may be found in the form
\[
 \psi = \tfrac12\omega r^2 g(\theta) 
\]
and find the function $g(\theta)$.
Deduce the pressure field $p(r,\theta)$. Discuss the limitations of
the model. [{\sl Does denominator of $g$ vanish?}]
\end{ex}

%%%%%%%%%%%%%%%%
\begin{ex}[Channel entry]\normalfont
Viscous fluid is contained between two planes $y=\pm a$ and a
two-dimensional flow with streamfunction $\psi(x,y)$ is generated by
some agency (e.g.~a rotating cylinder) near $x=y=0$. 
It is required to find the form of the flow field for large positive
$x$. 
Find the general solution of $\nabla^4\psi=0$ of the form
\[
 \psi = f(y)e^{-kx} \qquad {\rm Re}\, k > 0 ,
\]
for which $f(y)$ is an even function of $y$, and hence show that $k$ is
determined by the equation
\[
 2ka + \sin 2ka = 0.
\] 
Show that this equation as no real roots. The equation has complex
roots, that with the smallest real part being 
$2ka = 4.2 \pm 2.3i$. Sketch the streamlines of the flow.
\end{ex}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\chapter{Flows in thin layers}
\label{chpt:thinlayer}

The is variously known as {\em lubrication theory, slowly varying, and the long wavelength approximation\/}. All the applications are essentially Poiseuille flow at low Reynolds numbers. We will follow Reynolds' approach. An alternative approach based on asymptotic expansions works to produce the same answer, but with many more pages of algebra. Above all, these flows are very useful.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Thrust bearing}
\label{sec:thrust}

This problem has simple algebra. It has applications of gear teeth and knee cartilage.
\cpicture{0.8}{Thrust}
The problem is easiest tackled in a fixed geometry, i.e.~move with lower boundary.


There are 4 steps.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Geometry}

The gap is
   \[
     h(x) = d_1 + \alpha x \qword{with} \alpha = \frac{d_2-d_1}{L}.
   \]
To be a thin film, we require $h\ll L$ and $h'\ll1$.  

%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Unidirectional flow}   

At leading order the flow is $\bfu\sim(u(y),0,0)$ so that $\bfu\cdot\nabla\bfu\sim0$. Then the momentum equations reduces to
\begin{alignat*}{2}
  \hbox{\rm $x$-mtm:}\quad 0 &= \pd{p}{x} + \mu\pd{^2u}{y^2},
      &\quad\qword{because no} \pd{^2}{x^2} \qword{by} h\ll L, \\
  \hbox{\rm $y$-mtm:}\quad 0 &= \pd{p}{y},
      &\quad\qword{because flow}\sim (u(y),0,0).
\end{alignat*}
Hence to this approximation $p$ is independent of $y$, so $\partial p/\partial x$ is independent of $y$, which helps the integration of the $x$-momentum equation. We need boundary conditions
\[
 u = -U \qword{on} y=0, \quad u=0 \qword{on} h=h(x).
\]
So integrating the $x$-momentum equation
\[
 u = -\frac{1}{2\mu} \frac{dp}{dx} y(h-y) - U\frac{h-y}{h}.
\]
In this, we still need to find $dp(x)/dx$.  

%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Mass conservation}

This is the clever trick of Reynolds to avoid a tedious asymptotic analysis. The total volume flux (per unit $z$-direction)
\[
 Q = \int_0^hu\,dy = -\frac{h^3}{12\mu}\frac{dp}{dx} - \tfrac12 Uh.
\]
Global mass conservation requires $Q$ to be constant, independent of $x$. We do not yet know value of $Q$, but that it is constant gives the $x$-variation of $dp/dx$
\[
 \frac{dp}{dx} = -\frac{12\mu Q}{h^3(x)} - \frac{6\mu U}{h^2(x)}.
\]
Integrating from $p=p_0$ at $x=0$, with $h(x) = d_1 + \alpha x$
\[
 p(x) = p_0
     + \frac{6\mu Q}{\alpha}\left(\frac{1}{h^2(x)} - \frac{1}{d_1^2}\right)
     + \frac{6\mu U}{\alpha}\left(\frac{1}{h(x)} - \frac{1}{d_1}\right).
\]
This gives the variation of $p$ with $x$ with one undetermined constant, $Q$. But $p=p_0$ also at $x=L$, so
\[
  Q = - \frac{Ud_1d_2}{d_1 + d_2}.
\]
Note with this $Q$, $p$ is a maximum, $dp/dx=0$, when $h(x) = 2d_1d_2/(d_1+d_2)$. A parabolic component of the flow is driven by the pressure gradient away from the pressure maximum, in order to keep the total flux constant, see figure.
\cpicture{0.8}{Thrust2}
When $d_2\gg d_1$, $Q=-Ud_1$, the maximum pressure is at $h(x)=2d_1$ and
\begin{align*}
 \mu\left.\pd{u}{y}\right\vert_{y=0} &\qword{reverses sign at} h=\tfrac32d_1, \\
 \mu\left.\pd{u}{y}\right\vert_{y=h} &\qword{reverses sign at} h= 3d_1.
\end{align*}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Forces}

The forces are calculated from the surface tractions $\sigma_{ij}n_j$ with
\[
\sigma_{ij}= -p\delta_{ij} + \mu\left(\pd{u_i}{x_j} + \pd{u_j}{x_i}\right).
\]
On the lower boundary $\bfn = (0,-1,0)$.
But be warned there is {\em danger\/} on the upper boundary, where it is necessary to take account of the small slope, so $\bfn = (-h',1,0)$, where a denominator of $\sqrt{1 + h'^2}$ has been dropped because $h'^2\ll 1$.
The normal force on the lower plate (per unit $z$-width, and in excess of that from the atmospheric pressure $p_0$) is
\[
 \int_0^L (p-p_0)\,dx = \frac{6\mu U}{\alpha^2}
 \left(\ln\frac{d_2}{d_1} - 2\frac{d_2-d_1}{d_2+d_1}  \right).
\]
The tangential force on the lower plate (per unit $z$-width) is
\begin{gather*}
 \int_0^L \mu\left.\pd{u}{y}\right\vert_{y=0}\,dx 
     = \int_0^L \left(-\frac{h}{2}\frac{dp}{dx} + \frac{\mu U}{h} \right)\,dx \\
  = \frac{4\mu U}{\alpha}\left(\ln\frac{d_2}{d_1} - \frac32\frac{d_2-d_1}{d_2+d_1} \right).
\end{gather*}
The equal and opposite normal force on the top plate is straight forward to calculate. For the tangential force on the upper plate, one must include a contribution from the large pressure multiplied by the small slope $h'$.

Note that the $1/\alpha^2$ normal force is must larger than the $1/\alpha$ tangential force ($\alpha\ll 1$). Thus one can make low friction bearings -- all due to the large pressure built up over the long (thin) layer.
{\color{darkgray}
In real gears the pressure is so high the the solid surfaces deform elastically, typically doubling the gap, and the viscosity changes with pressure, typically exponentially.
}%end color

\bigskip
The condition to neglect the inertial terms can be refined to
\[
 1 \gg \frac{\vert \rho\bfu\cdot\nabla\bfu\vert}{\vert\mu\nabla^2\bfu\vert}
     = \frac{\rho U^2/L}{\mu U/h^2}  = \frac{Uh}{\nu}\frac{h}{L},
\]
so it is possible to have a moderate Reynolds number ${\rm Re} = Uh/\nu$ if the layer is compensatingly thin $h\ll L$. The combination $Re\, h/L$ is called the ``reduced Reynolds number''.  

The geometry can be wrapped round a circle to make a journal bearing
 \cpicture{0.4}{Journal}
where the pressure variation ensures a constant flux around the bearing.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Cylinder approaching a wall}

\medskip
\cpicture{0.7}{WallCylinder}
Instantaneous separation $d \ll a$.   

Repeating the four key steps of \S\ref{sec:thrust}, now more briefly.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Geometry of thin gap}

\begin{align*}
  h(x) &= d + a\left(1-\cos\theta \sim \tfrac12\theta^2
    \sim \tfrac12\left(\frac{x}{a}\right)^2  \right), \\
  &\sim d\left(1 + \frac{x^2}{2ad} \right).
\end{align*}
We shall find that all the action is where $h \lesssim 3d$, i.e.~where
\[
 x \sim 2\sqrt{ad} \begin{cases}
                     \ll a, & \hbox{so $\theta$ small, so above approx OK,} \\
                     \gg d, & \hbox{so approx unidirectional flow.}
                   \end{cases}
\]

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Unidirectional flow}

Now the boundary condition is $u=0$ on both $y=0$ and $h(x)$. So as before
\[
  u = - \frac{1}{2\mu} \frac{dp}{dx} y(h-y).
\]

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Mass conservation}

The total volume flux (per unit $z$-width), as before
\[
 Q = \int_0^h u\,dy = - \frac{h^3}{12\mu} \frac{dp}{dx}.
\]
Now the top boundary is moving down at $V$, so that the volume flux out of $[0,x]$ is
\[
   Q(x) = Vx.
\]
Hence we have the pressure gradient
\[
 \frac{dp}{dx} = - \frac{12\mu Vx}{h^3(x)}
     = - \frac{12\mu Vx}{d^3\left(1 + \frac{x^2}{2ad} \right)^3}.
\]
Integrating
\[
    p = p_0 + \frac{6\mu Va}{d^2\left(1 + \frac{x^2}{2ad} \right)^2}.
\]
Here we can see that most of the pressure is where $h\lesssim 3d$, as claimed earlier.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Forces}

The normal force (per unit $z$-width) on lower flat plate in excess of that from atmospheric pressure $p_0$ is
\begin{align*}
  \int_{-\infty}^\infty (p&-p_0)\,dx, \qword{where ``$\infty$'' is where $h\gtrsim 3d$} \\
   &= \frac{6\mu Va}{d^2}\sqrt{2ad}\left(\int_{-\infty}^\infty \frac{d\xi}{(1 + \xi^2)^2} \quad = \pi\right) \\
   &= 6\sqrt2\pi\left(\frac{a}{d}\right)^{3/2}\mu V. 
\end{align*}
If the cylinder is falling under own weight, $F=Mg$, then $\dot{d} = -V$, so $d\propto t^{-2}$, and it takes an infinite time to touch
{\color{darkgray}(but surface roughness)}.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Other thin layers}

\begin{itemize}
 \item Journal bearing,
 \item Sphere approaching a wall,
 \item Cylinder falling parallel to a vertical wall (does not rotate!),
 \item Two cylinders counter rotating, for milling steel,
 \item Thin layers on inclined plane, for lava flow,
 \item Thin layer dripping off a rotating horizontal cylinder (Moffatt).
\end{itemize}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Hele-Shaw cell}

This device was invented to visualise two-dimensional potential flow for the design of ships.

\cpicture{0.7}{HeleShaw}

Flow is in the thin gap of thickness $h$ between two fixed plates, with obstacles of size $L$. If $h \ll L$, the flow is mostly in $xy$-plane.
The no-slip boundary condition is $\bfu=0$ on $z=0$ and $z=h$. Then
\begin{align*}
  \hbox{$\bfu$ varies on $h$ } &\hbox{in $z$-direction ($\perp$ to flow)}, \\
  \hbox{$\bfu$ varies on $L$ } &\hbox{in $x$ and $y$ directions ($\parallel$ to flow).}   
\end{align*}
Then we can ignore inertia if
\[
  \frac{\vert \rho\bfu\cdot\nabla\bfu\vert}{\vert\mu\nabla^2\bfu\vert}
     = \frac{\rho U^2/L}{\mu U/h^2}  = \frac{Uh}{\nu}\frac{h}{L} \ll 1.
\]
The momentum equation then becomes, by $\partial_x, \partial_y \sim 1/L \ll 1/h \sim \partial_z$,
\[
 0 = -\nabla p + \mu \pd{^2\bfu}{z^2}.
\]
But $\bfu \sim (u(x,y,z),v(x,y,z),0)$ do $\partial p/\partial z =0$. So $\nabla p$ is in the $xy$-plane and independent of $z$. Hence
\[
 \bfu = -\frac{1}{2\mu}\nabla_2 p\; z(h-z),
\]
where $\nabla_2 = (\partial_x,\partial_y,0)$. Integrating through the depth, the volume flux is
\[
 {\bf q}_2 = \int_0^h \bfu\,dz = - \frac{h^3}{12\mu}\nabla_2p.
\]
Integrating
\[
 0 = \nabla\cdot\bfu = \pd{u}{x} + \pd{v}{y} + \pd{w}{z}
\]
across the depth
\[
0 = \pd{}{x}\left(\int_0^hu\,dz\right) + \pd{}{y}\left(\int_0^hv\,dz\right) + w(h) - w(0).
\]
This requires $h$ to be independent of $x$ and $y$, although the generalisations work. Now the vertical velocity $w$ vanishes on the boundaries $z=0$ and $z=h$ so that the last two terms are zero. Hence for the two fixed flat plates
\[
 \nabla_2\cdot{\bf q}_2 = 0.
\]{}
Substituting in the expression for ${\bf q}_2$, and taking $h$ and $\mu$ to be constant,
\[
 \nabla_2^2p = 0.
\]
Hence although the flow is dominating by viscosity, the flow is a (IB) potential flow $\bfu \propto \nabla p$ with $\nabla^2 p = 0$. Hence the Hele-Shaw device is used to visualise two-dimensional potential flows, e.g.~aerofoils.
Note although $\bfu$ varies in magnitude with $z$, it has the same direction for all $z$, which is essential when using dye to mark the flow pattern.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Flow past a circular cylinder}

Now
\[
 {\bf q}_2 = -\frac{h^3}{12\mu}\nabla_2\,p \qword{and} \nabla_2\cdot{\bf q}_2 = 0, \qword{so} \nabla_2^2 p = 0 \qword{in} r \ge a.
\]
There is uniform flow at large distances
\[
 {\bf q}_2 \to (Q,0) \qword{in Cartesian coordinates as $r\to\infty$}.
\]
Hence
\[
 p \to - \frac{12\mu Q}{h^3} r\cos\theta.
\]
The condition of no flow through the cylinder is
\[
  {\bf q}_2\cdot\bfn = 0 \qword{on} r=a, \qword{so} q_r=0,
  \qword{so} \pd{p}{r} = 0.
\]
Hence
\[
 p = - \frac{12\mu Q}{h^3}\cos\theta\left( r + \frac{a^2}{r}\right).
\]

The same equation governs for similar reasons seepage flows through porous rock strata, as in oil reservoirs, aquifers and filters.
  

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Gravitational spreading on a horizontal surface}

\cpicture{0.7}{SpreadingDrop}
How does a liquid drop spread under gravity? -- due to $\partial h/\partial r$.


We assume (i) axisymmetry, (ii) a thin layer $\vert\partial h/\partial r\vert \ll1$, (iii) a viscous flow so that we can ignore inertia (Reduced Reynolds number small), (iv)~negligble surface tension and no problems of the moving contact line, and (v) finally the table is horizontal.

The motion is mainly horizontal, so that the $z$ momentum equation is
\[
 0 = - \pd{p}{z} - \rho g + \mu 0,
\]
with the free surface boundary condition $p=p_0$ at $z=h(r,t)$. Hence the pressure distribution is hydrostatic
\[
 p(r,z,t) = p_0 + \rho g(h-z).
\]
The higher pressure on the the table ($z=0$) in the centre leads to the spreading. The $r$-momentum equation is
\[
 0 = - \pd{p}{r} + \mu\pd{^2u}{z^2},
\]
with no slip boundary condition $u=0$ at $z=0$ and the free surface boundary condition of no tangential stress $\mu\partial u/\partial z = 0$ on $z=h$. Now our hydrostatic pressure gives
\[
 \pd{^2u}{z^2} = \frac{1}{\mu}\pd{p}{r} = \frac{\rho g}{\mu}\pd{h}{r}.
\]
But $\partial h/\partial r$ is independent of $z$ so
\[
 u = - \frac{g}{2\nu} \pd{h}{r}\, z(2h-z).
\]
The volume flux out of a cylinder of radius $r$ is
\[
 Q(r) = 2\pi r \int_0^h u\,dz = -2\pi r \frac{gh^3}{3\nu}\pd{h}{r}.
 \]
 Mass conservation between cylinders $r$ and $r+\delta r$ is

\begin{minipage}{0.45\linewidth}
\[
 2\pi\delta r \pd{h}{t} + Q(r+\delta r) - Q(r) = 0,
 \]
 so
 \[
 2\pi r \pd{h}{t} + \pd{Q}{r} = 0.
 \]
\end{minipage}
\begin{minipage}{0.45\linewidth}

 \cpicture{0.8}{SpreadingDrop2}

\end{minipage}

\medskip
\noindent Hence subsituting in the expression for $Q$, we have a nonlinear  diffusion equation governing the spreading of the drop
\[
 \pd{h}{t} = \frac{g}{3\nu}\frac{1}{r}\pd{}{r}\left(r h^3\pd{h}{r}\right).
\]
 
\noindent Note the condition for the edge of the spreading drop is
\[
 h=0 \qword{at} r=a(t), \qword{so there} Q=0.
\]

\noindent Multiplying the above diffusion equation by $r$ and integrating from $r=0$ to $r=a$ gives
\[
 \int_0^a r\pd{h}{t}\,dr = \left. \frac{g}{3\nu}rh^3\pd{h}{r}\right\vert_0^a = 0.
\]
Hence
\[
 \int_0^a h2\pi r\,dr = V, \qword{a constant volume of the drop.}
\]

\medskip
\noindent The nonlinear diffusion equation (with interesting behaviour where $h\to0$) can be solved easily numerically as a initial value problem. One discovers that all solutions tend in time to a {\em similarity solution\/} (of self-similar form)
\[
 h(r,t) = H(t) f(\eta) \qword{where} \eta = \frac{r}{R(t)}.
\]
First we find the scaling functions for the time-varying height $H(t)$ and radius $R(t)$, and then the ordinary differential equation for the shape $f(\eta)$. Global mass conservations gives (dropping all numerical factors)
\[
 HR^2 = V,
\]
while the governing diffusion eqaution gives
\[
 \frac{H}{t} = \frac{g}{\nu}\frac{H^4}{R^2}.
\]
Hence
\[
 H = \left(\frac{\nu V}{gt} \right)^{1/4} \qquad
 R = \left(\frac{gtV^3}{\nu} \right)^{1/8}. 
\]
Substituting into the governing diffusion equation
\[
 - \frac{H}{4t} + H\frac{df}{d\eta}\frac{-\eta}{8t}
  = \frac{g}{3\nu} \frac{H^4}{R^2} \frac{1}{\eta}\frac{d}{d\eta}
   \left(\eta f^3 \frac{df}{d\eta}  \right).
\]
i.e.
\[
 - \tfrac{1}{4}f - \tfrac{1}{8}\eta f' = \tfrac{1}{3}\frac{1}{\eta}
  \left(\eta f^3 f'  \right)'.
\]
To be solved subject to the volume constraint
\[
 2\pi\int_0^{\eta_0} f\eta\,d\eta = 1, \qword{where} f(\eta_0) = 0.
\]
At this stage it is common to have to solve the ordinary equation numerically, but in this example it can be solved analytically. If one first multiplies by $\eta$, the equation has a first integral
\[
 -\tfrac{1}{8} \eta^2 f = \tfrac{1}{3}\eta f^3 f' + \hbox{const}.
\]
The constant vanishes because $f'$ is finite at $f=0$.
Hence
\[
 \tfrac{1}{3} f^2 f' = - \tfrac{1}{8}\eta, \qword{so}
 \tfrac{1}{9} f^3 = \tfrac{1}{16}(\eta_0^2-\eta^2), \qword{so}
 f = \left(\tfrac{9}{16}\right)^{1/3}
 (\eta_0^2-\eta^2)^{1/3}.
\]
Finally the volume constraint is
\[
 1 = 2\pi\int_0^{\eta_0}f\eta\,d\eta = 2\pi\left(\tfrac{9}{16}\right)^{1/3}
   \left[\int_0^{\eta_0} (\eta_0^2-\eta^2)^{1/3} \eta\,d\eta = \tfrac{3}{8}\eta_0^{8/3}\right]
\]
Hence
\[
 \eta_0 = \left(\frac{2^{10}}{3^5\pi^3} \right){1/8} = 0.7792,
\]
and the drop spreads according to
\[
 r = \eta_0 \left(\frac{gtV^3}{\nu}\right)^{1/8}.
\]
The $1/8$ power-law for the time dependence means increasing slow spreadins, e.g.~if it takes 1 second for the first 1 cm, it will take 4 minutes to go to 2 cm and 18 hours to go to 4 cm.


\bigskip
Note that the assumption of a similarity shape has turned the partial differential equation into an ordinary differential equation. Similarity solutions may or may not exist, and may or may not be stable.

\bigskip
%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{ex}[Sphere approaching a wall]\normalfont
A rigid sphere of radius $a$ falls through a fluid of viscosity
$\mu$ under gravity towards a horizontal rigid plane.
Use lubrication theory to show that, when the minimum gap $h_0$ is
very small, the speed of approach of the sphere is
\[
 h_0W/6\pi\mu a^2, 
\]
where $W$ is the weight of the sphere corrected for buoyancy.
\end{ex}

%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{ex}[Flow through a gap]\normalfont
Oil is forced by a pressure difference $\Delta p$ through the
narrow gap between two parallel circular cylinders of radius $a$ with
axes $2a+b$ apart.
Show that, provided $b\ll a$ and $\rho b^3\Delta p \ll \mu^2 a$, the
volume flux is approximately
\[
 \frac{2b^{5/2}\Delta p}{9\pi a^{1/2} \mu}
\]
when the cylinders are fixed.

Show also that when the two cylinders rotate with angular velocities
$\Omega_1$ and $\Omega_2$ in opposite directions, the change in the
volume flux is
\[
 \tfrac23 ab(\Omega_1+ \Omega_2) .
\]
\end{ex}

%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{ex}[Hovering on an air-table]\normalfont
A disk hovers on a cushion of air above an air-table -- a fine porous plate through which a constant flux of air is pumped.
Let the disk have a radius $R$ and a weight $Mg$ and hover at a low height $h$ ($h \ll R$) above the air-table.
Let the volume flux of air, which has density $\rho$ and viscosity $\mu$, be $w_0$ across unit surface area.
The conditions are such that $\rho w_oh^2/\mu R\ll 1$.
Explain the significance of this restriction.

Find the pressure distribution in the air under disk.
Show that this pressure can balance the weight of the disk if
\[
 h = R \left(\frac{3\pi\mu R w_o }{2Mg}\right)^{1/3}.
\] 
\end{ex}

%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{ex}[Dripping from rotating horizontal cylinder]\normalfont
A viscous fluid coats the outer surface of a cylinder of radius
$a$ which rotates with angular velocity $\Omega$ about its axis which
is horizontal.
The angle $\theta$ is measured from the horizontal on the rising side.
Show that the volume flux per unit length $Q(\theta,t)$ is related to
the thickness $h(\theta,t)$ of the fluid layer by
\[
 Q = \Omega a h - \frac{g}{3\nu}h^3\cos\theta ,
\] 
and deduce an evolution equation for $h(\theta,t)$.

Consider now the possibility of a steady state with $Q=\rm const$,
$h=h(\theta)$. 
Show that a steady solution with $h(\theta)$ continuous and 
2$\pi$-periodic
exists only if
\[
 \Omega a > (9Q^2g/4\nu)^{1/3} .
\]
\end{ex}

%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{ex}[Spreading of 2D drop]\normalfont
A two-dimensional drop $h(x,t)$ spreads on a horizontal table.
Assuming that the drop has become a thin layer, find how the drops spreads. 
[It is not possible to integrate the volume in closed form.]
\end{ex}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\chapter{Vorticity generation and confinement}
\label{chpt:vort}


We now switch to high Reynolds numbers. While inertia dominates, viscosity has a role to play, particularly in special thin regions.
This chapter is a preparation for the following chapter on boundary layers.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Vorticity equation}

The vorticity $\bfom = \nabla\times\bfu$ is the local angular velocity, see \S\ref{sec:strainrate}. Hence the vorticity equation is the angular momentum equation for fluid motion.

The momentum equation (Navier-Stokes), see \S\ref{sec:NSE}, is
\[
 \rho\frac{D\bfu}{Dt} = \bfF - \nabla p + \mu\nabla^2\bfu.
\]
The fluid acceleration may be written as
\[ 
 \frac{D\bfu}{Dt} = \pd{\bfu}{t} + (\bfu\cdot\nabla)\bfu
  = \pd{\bfu}{t} + \nabla\left(\tfrac12 u^2\right) - \bfu\times(\nabla\times\bfu) .
\]
We now assume that the force is conservation
\[  
  \bfF = - \nabla\Phi,
\]
and that $\rho$ and $\mu$ are constants.\hfil
\break\noindent Taking the curl of the momentum equation, we have
\[ 
  \rho\left(\pd{\bfom}{t} - \nabla\times(\bfu\times\bfom)\right)
  = \mu\nabla^2\bfom
\]
Now
\[
 \nabla\times(\bfu\times\bfom)
  = (\bfom\cdot\nabla)\bfu + \bfu(\nabla\cdot\bfom)
  - (\bfu\cdot\nabla)\bfom - \bfom(\nabla\cdot\bfu)
\]
The two divergences vanish, because $\bfom$ is a curl and because the flow is incompressible . Hence the vorticity equation
\[
 \boxed{
  \pd{\bfom}{t} + \underset{(1)}{(\bfu\cdot\nabla)\bfom}
   = \underset{(2)}{(\bfom\cdot\nabla)\bfu} + \underset{(3)}{\mu\nabla^2\bfom}.
       }
\]
Term (1) is the advection of vorticity with the flow. Term (2) is stretching of vorticity. Term (3) is the diffusion of vorticity. To build up an understanding of the vorticity equation, we will examine in the following sections first stretching, then diffusion, then balancing advection and diffusion, and finally a combination of advection, stretching and diffusion.

Note with the assumption of conservative forces and $\rho$ constant, there is no source of vorticity in the interior of the flow.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Stretching vorticity}

Also called the Ballerina effect, and important for tornadoes and hurricanes.

Recall the motion of a material line element $\boldsymbol{\delta\ell}$ moving with the fluid. The two ends at $\bfx$ and $\bfx+\boldsymbol{\delta\ell}$ move in $\delta t$ with the flow according to
%\cpicture{0.8}{LineElement}
\begin{align*}
  \bfx  \quad&\to\quad \bfx + \delta t\,\bfu(\bfx,t), \\
  \bfx+\boldsymbol{\delta\ell}  \quad&\to\quad
  \bfx +\boldsymbol{\delta\ell} + \delta t\,\bfu(\bfx+\boldsymbol{\delta\ell},t)
\end{align*}
Hence for small $\boldsymbol{\delta\ell}$
\[
  \boldsymbol{\delta\ell} \quad\to\quad 
  \boldsymbol{\delta\ell} + \delta t\, (\boldsymbol{\delta\ell}\cdot\nabla)\bfu.
\]
or
\[
 \frac{d}{dt}\boldsymbol{\delta\ell} = (\boldsymbol{\delta\ell}\cdot\nabla)\bfu.
\]
The component of the right hand side perpendicular to $\boldsymbol{\delta\ell}$ rotates $\boldsymbol{\delta\ell}$, while the component parallel stretches it.
Comparing with the vorticity equation, we see that vorticity is rotated and stretched just like a material line element.

The stretching is just the conservation of angular momentum.
Consider a cylinder of radius $a_1$ and length $\ell_1$ spinning at $\omega_1$ which becomes $a_2$, $\ell_2$ and $\omega_2$
\cpicture{0.6}{Ballerina}
Mass conservation gives
\[
  \rho a_1^2\ell_1 = \rho a_2^2\ell_2,
\]
while angular momentum conservation gives
\[
  \rho a_1^2\ell_1 a_1^2\omega_1 = \rho a_2^2\ell_2 a_2^2\omega_2.
\]
Hence 
\[
 \omega_2 = \omega_1\frac{\ell_2}{\ell_1}.
\]

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Diffusion of vorticity}

In \S\ref{sec:imp} we studied a plate impulsively started, by solving the diffusion equation
\[
 \pd{u}{t} = \nu \pd{^2u}{y^2},
\]
and found that momentum diffused from the plate, influencing ($u\gtrsim \tfrac15 U$) a region next to the plate $y \lesssim \sqrt{\nu t}$.

Alternatively, we can view the initial condition as a discontinuity in $u(y)$ on the boundary, i.e.~a delta function of vorticity $\omega = - \partial u/\partial y$. So we could say that vorticity was generate at the boundary at $t=0+$ and then diffused into the interior
\[
 \pd{\omega}{t} = \nu\pd{^2\omega}{y^2}.
\]

In \S\ref{sec:osc} we studied an oscillating plate, and found that alternating signed momentum diffused into the interior to a depth $y=\sqrt{2\nu/\omega} \approx \sqrt{\nu t}$ with $t = \pi/\omega$. Alternatively, we can view the flow as alternating signed vorticity diffusing from the boundary into the interior.

In both cases, vorticity is generated at the boundary, not in the interior.


%%%%%%%%%%%%%%%%%%%
\section{Wall with suction}

This is an artificial, but simple, example of the confinement of vorticity. It does have some relevance to cross-flow filtration and boundary-layer control.

The boundary is a porous plate with a low pressure on its other side which produces a suction of a volume flux $V$ per unit area.

\cpicture{0.6}{Suction}

\noindent Far from the wall, there is a cross-flow $U$ and the suction $V$.
We take $\bfF=0$ and $\nabla p=0$.
We try a flow field $\bfu=(u(y),-V,0)$ with $u\to U$ as $y\to\infty$, and no-slip $u=0$ on $y=0$.

\noindent With this assumed form of the flow,  mass conservation $\nabla\cdot\bfu=0$ is automatically satisfied.
The $x$-momentum equation is
\[
 -\rho V \pd{u}{y} = \mu \pd{^2u}{y^2},
\]
with solution
\[
  u = U\left(1 - e^{-Vy/\nu} \right).
\]
Thus vorticity $-\partial u/\partial y$ is confined to a region next to the wall, of thickness $\delta = \nu/V$.

Vorticity diffuses away from the wall, in time $t$ to a distance $\sqrt{\nu t}$, but is advected back by the suction flow, moving a distance $Vt$. Balancing diffusion away with advection back occurs at
\[
 t = \nu/V^2 \qword{i.e.~at thickness} \delta = \nu/V.
\]
Alternatively one can balance terms in the momentum equation
\[
 -\rho V\pd{u}{y} = O\left(\frac{\rho VU}{\delta}\right), \quad
  \mu\pd{^2u}{y^2} = O(\left(\frac{\mu U}{\delta^2} \right),
\]
balancing gives $\delta=\nu/V$.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Stagnation point flow on rigid boundary}
\label{sec:stag}

This is a second example of the confinement of vorticity, in which there is a balance between vorticity diffusing away from the boundary and being advected back to it. For flow $U$ past an obstacle size $L$, some flow will pass one side and some the other. The dividing streamline reaches the obstacle at a stagnation point where the velocity vanishes, and there the velocity gradient (strain-rate) is $E = O(U/L)$. In this example, we examine the flow in the neighbourhood of the stagnation point, where we may take the boundary to be flat. We also consider the two-dimensional case.

\cpicture{0.8}{Stagnation}

\noindent We take a pure straining field far from the plate
\[
 \bfu \to (Ex,-Ey,0) \qword{as} y \to\infty,
\]
with strain-rate $E$, and no slip $\bfu=0$ at $y=0$. With a pure straining flow at large distances, there is no vorticity there.

We estimate the thickness of the region of confined vorticity by balancing terms in the vorticity equation
\[
 \rho v\pd{u}{y} = O\left(\rho E\delta \frac{Ex}{\delta}\right), \quad
  \mu\pd{^2u}{y^2} = O\left(\mu\frac{Ex}{\delta^2} \right),
\]
Hence
\[
 \delta = \sqrt{\frac{\nu}{E}}.
\]
This is just $\sqrt{\nu t}$ again with $t=1/E$. Because $\delta$ does not depend on $x$, we will obtain an exact solution of the Navier-Stokes equation, which is quite rare.

We non-dimensionalise the problem using $\delta$ to scale $y$, $Ex$ to scale $u$, $E\delta$ to scale $v$, i.e.~we seek a solution of the form
\[
 u = Exf(\eta), \quad v=E\delta g(\eta), \qword{with} \eta = y/\delta.
\]
Mass conservation gives
\[
 0 = \pd{u}{x} + \pd{v}{y} = Ef + E\delta g' \frac1\delta,
\]
i.e.~is satsified for all $x$ and $y$ if
\[
 f + g' = 0.
\]
The $x$-momentum equation is
\begin{gather*}
  \rho\left(u\pd{u}{x} + v\pd{u}{y}\right) = -\pd{p}{x}
    + \mu\left(\pd{^2u}{x^2} + \pd{^2u}{y^2}\right), \\
  \rho\left(Exf\,Ef + E\delta g\,Exf'\frac1{\delta}\right) = - \pd{p}{x},
    + \mu\left(0 + Ex f''\frac1{\delta^2} \right), \\
  \rho E^2 x\left(f^2 + gf' - f'' \right) = - \pd{p}{x},  
\end{gather*}
where we have used $\delta^2=\nu/E$.
The $y$-momentum is
\begin{gather*}
  \rho\left(u\pd{v}{x} + v\pd{v}{y}\right) = -\pd{p}{y}
    + \mu\left(\pd{^2v}{x^2} + \pd{^2v}{y^2}\right), \\
  \rho\left(Exf\,0 + E\delta g\,E\delta g'\frac1\delta\right) = -\pd{p}{y}
    + \mu\left(0 + E\delta g''\frac{1}{\delta^2} \right), \\
  \rho E^2\delta\left(gg' -g'' \right) = - \pd{p}{y}.
\end{gather*}
Eliminating the pressure by cross-differentiating (forming the vorticity equation), one finds that both momentum equations can be satisfied for all $x$ and $y$ if
\[
 (f^2 + gf' - f'')' = 0,
\]
as $\partial_x(p_y) = 0$.
The boundary and far conditions are
\[
 \qword{as} \eta\to\infty, \quad f\to1, \quad g \to -\eta,
  \qword{and at} \eta=0, \quad f=0, \quad g=0.
\]
Then integrating once, using $f\to1$ as $\eta\to\infty$,
\[
  f^2 + gf' - f'' = 1.
\]
Substituting in $f=-g'$, we have
\[
 g'^2 - gg'' + g''' = 1.
 \]
 This has to be integrated numerically. One finds
 \cpicture{0.8}{Stagnation2}
 The horizontal velocity $u$ is within 1\% of the far field value by $y=2.4\delta$. The behaviour $g\sim-(\eta-0.65)$ as $\eta\to\infty$ means that the flow looks like an inviscid pure straining with slip at a wall at $y=0.65\delta$. This is called the displacement thickness, corresponding to a deficit in the volume flux because of the fluid slowed down next to the wall.

 The pressure is
 \[
 p = \tfrac12\rho E^2x^2 + \tfrac12 E^2\delta^2\int_0^\eta (g''-gg') \,d\eta.
 \]

 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 \section{Burgers vortex}

This is another exact solution of the Navier-Stokes equation. It is steady, and has a balance between advection, stretching and diffusion of vorticity.
\cpicture{0.6}{Burgers}
The flow is a combination of an axisymmetric swirling flow around an axis and an axisymmetric straining motion, in polar coordinates
\[
 \bfu = (0,v(r),0) + \alpha(-r,0,2z),
\]
or in Cartesian coordinates
\[
 \bfu = (-\tfrac{y}{r}v,\tfrac{x}{r}v,0) + \alpha(-x,-y,2z).
\]
This combination satisfies mass conservation, in polar coordinates
\[
 \nabla\cdot\bfu = \frac{1}{r}\pd{}{r}(ru_r) + \frac{1}{r}\pd{u_\theta}{\theta} + \pd{u_z}{z} = -2\alpha + 0 + 2\alpha = 0,
\]
or in Cartesian coordinates
\[
 \pd{u_x}{x} + \pd{u_y}{y} + \pd{u_z}{z} =
  -y\left(\frac{v}{r}\right)'\frac{x}r - \alpha +
  x\left(\frac{v}{r}\right)'\frac{y}r - \alpha + 2\alpha = 0.
\]
The vorticity comes only from the swirling motion, in polars
\[
 \bfom = \nabla\times\bfu = (0,0,\frac{1}{r}\pd{}{r}(rv)) = (0,0,\omega(r)),
\]
or in Cartesian the $z$-component is
\[
 \omega = \pd{}{x}\left(\frac{x}{r}v\right)- \pd{}{y}\left(-\frac{y}{r}v\right)
 = \frac{2v}{r} + \frac{x^2+y^2}{r}\pd{}{r}\left(\frac{v}{r}\right)
  = \frac{1}{r}\pd{}{r}(rv).
\]
The vorticity equation
\[
 \pd{\bfom}{t} + (\bfu\cdot\nabla)\bfom = (\bfom\cdot\nabla)\bfu + \nu\nabla^2\bfom
\]
only has a non-zero component in the $z$-direction and $\omega$ only varies with $r$, so
\[
 \underset{0\, {\rm steady}}{\pd{\omega}{t}}
 -\underset{{\rm advection}}{\alpha r\pd{\omega}{r}}
 =\underset{{\rm stretching}}{\omega\pd{}{z}(2\alpha z)}
 +\underset{{\rm diffusion}}{\nu\frac{1}{r}\pd{}{r}\left(r\pd{\omega}{r}\right)},
\]
i.e.~multiplying by $r$
\[
 0 = \alpha r^2\pd{\omega}{r} + 2\alpha r\omega + \nu\pd{}{r}\left(r\pd{\omega}{r} \right).
\]
Integrating
\[
 0 = \alpha r^2 \omega + \nu r\pd{\omega}{r} + \hbox{\rm const}.
\]
The constant vanishes by requiring the vortex to decay at large distance, $\omega\to0$ as $r\to\infty$, so
\[
 \omega = \omega_0 e^{-\alpha r^2/2\nu}.
\]
This distribution of vorticity represents a steady balance between advection, stretching and diffusion.

Now we have above an expression relating the vorticity to the swirling velocity
\[
 \omega = \frac{1}{r}\pd{}{r}(rv), 
\]
so integrating
\[
 v = \frac{1}{r}\frac{\nu\omega_0}{\alpha}\left(1 - e^{-\alpha r^2/2nu}\right).
\]
Here a constant of integration has been chosen so that $v$ is finite as $r\to0$.
Near to the axis $r\ll\delta=\sqrt{\alpha/\nu}$,
\[
 v \sim \tfrac12\omega_0r,
\]
a solid body rotation, and far from the axis $r\gg\delta$
\[
 v \sim \frac{1}{r}\frac{\nu\omega_0}{\alpha},
\]
the potential flow of a line vortex.
Thus viscosity regularises the singular core of a line vortex.

In the Burgers vortex, again vorticity is confined by advection in against diffusion out
\[
 v_r\pd{\omega}{r} = O\left(\alpha\delta\frac{\omega}{\delta}\right), \quad
 \nu\pd{^2\omega}{r^2} = O\left(\nu\frac{\omega}{\delta^2}\right),
\]
balancing gives $\delta = \sqrt{\alpha/\nu}$. 

\bigskip
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{ex}[Channel with cross-flow]\normalfont
The walls of a channel are porous and separated by a distance $d$.
Fluid is driven through the channel by a pressure gradient 
$G = -\partial p/\partial x$, and at the same time suction is applied
to one wall of the channel providing a cross flow with uniform
transverse component of velocity $V$, fluid being supplied at this
rate at the other wall. 
Find and sketch the steady velocity and vorticity distributions in the
fluid (i) when $Vd/\nu\ll1$ and (ii) when $Vd/\nu\gg1$. 
\end{ex}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{ex}[Cylindrical annulus with cross-flow]\normalfont
Viscous fluid fills an annulus $a<r<b$ between a long stationary
cylinder $r=b$ and a long cylinder $r=a$ rotating at angular velocity
$\Omega$.
Find the axisymmetric velocity field, ignoring end effects.

Suppose now that the two cylinders are porous, and a pressure
difference is applied so that there is a radial flow
$-Va/r$.
Find the new steady flow around the cylinder when $Va/\nu<2$ 
and $Va/\nu>2$. 
Comment on the flow structure when \hbox{$Va/\nu\gg1$}.

Find the torque that must be applied to maintain the motion.
\end{ex}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{ex}[Diffusion of vorticity in time]\normalfont
Starting from the Navier-Stokes equations for incompressible
viscous flow with conservative forces, obtain the vorticity equation
\[
 \frac{D\bfom}{Dt} = \bfom\cdot\nabla\bfu + \nu \nabla^2\bfom .
\]
Interpret the terms in the equation.

At time $t=0$ a concentration of vorticity is created along the
$z$-axis, with the same circulation $\Gamma$ around the axis at each $z$.
The fluid is viscous and incompressible, and for $t>0$ has only an
azimuthal velocity $v$, say.
Show that there is a similarity solution of the form 
$vr/\Gamma = f(\eta)$, where $r=(x^2+y^2)^{1/2}$ and $\eta$ is a
suitable similarity variable.
Further show that all conditions are satisfied by
\[
 f(\eta) = \frac{1}{2\pi}(1-e^{-\eta}), \quad \eta = r^2/4\nu t .
\]
Show also that the total vorticity in the flow remains constant at
$\Gamma$ for all $t>0$. Sketch $v$ as a function of $r$.
\end{ex}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{ex}[Diffusion and stretching of vorticity in time]\normalfont
Calculate the vorticity $\bfom$ associated with the velocity
field 
\[
 \bfu = \left(-\alpha x - yf(r,t),\, -\alpha y + xf(r,t), \,2\alpha
 z\right), 
\]
where $\alpha$ is a positive constant, and $f(r,t)$ depends on 
$r=(x^2+y^2)^{1/2}$ and time $t$.
Hence show that the velocity field represents a dynamically possible
motion if $f(r,t)$ satisfies
\[
 2f + r\pd{f}{r} = A \gamma(t)e^{-\gamma(t)r^2},
\]
where
\[
 \gamma(t) = 
  \frac{\alpha}{2\nu}\left(1\pm e^{-2\alpha(t-t_0)}\right)^{-1},
\]
and $A$ and $t_0$ are constants.

Show that in the case where the minus sign is taken $\gamma$ is
approximately $1/(4\nu(t-t_0))$ when $t$ only just exceeds
$t_0$. Which terms in the vorticity equation dominate when this
approximation holds?
\end{ex}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\chapter{Boundary layers at high Reynolds numbers}
\label{chpt:bl}

In this chapter we shall see vorticity confined to layers very close to the wall, because in a fast flow it will not have time to diffuse far.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Euler and Prandtl limits}

The full incompressible Navier-Stokes equations are
\begin{gather*}
  \nabla\cdot\bfu = 0 \\
  \rho \frac{D\bfu}{Dt} = \bfF - \nabla p + \mu\nabla^2\bfu.
\end{gather*}
At high Reynolds numbers, $Re = UL/\nu \gg1$ for flow $U$ past a body of size $L$, {\em one might argue\/} that the viscous terms should be neglected, which would lead to
\[
   \rho \frac{D\bfu}{Dt} = \bfF - \nabla p,
\]
the Euler equation. This {\em argument is dangerous\/}. Mathematically, one has to give up one boundary condition, the no-slip condition $\bfu_{\rm tang}=0$ associated with viscous stress being finite.

The Euler limit is OK sometimes, for water waves, implusive motion, and lift on an aerofoil.

But not always. A major divergence had arisen in the 1900s between experiments and ``theoretical'', i.e.~Euler-equation, studies. The paradox was resolved by Prandtl in 1902.

{\em The error in the argument\/} is to assume that the flow adopts the lengthscale of the geometry, $L$. Instead one sees thin regions with large spatial gradients where $\mu\nabla^2\bfu$ (small $\mu$, large $\nabla^2$) is comparable with $\rho(\bfu\cdot\nabla)\bfu$. These thin layers permit one to apply the no-slip boundary condition.

The time to pass a body of length $L$ at velocity $U$ is $t=L/U$. During this time, vorticity generated at the boundary diffuses by viscosity (the influence of the no-slip boundary) a distance
\[
 \delta = \sqrt{\nu t} = \sqrt{\nu L/U} = L \left(\frac{UL}{\nu}\right)^{-1/2}.
\]
Thus at high Reynolds numbers $Re = UL/\nu \gg 1$, the vorticity is confined to a thin region (thin compared with the body). This chapter studies these thin regions, called {\em boundary layers\/}, in what is called the {\em Prandtl limit.\/} Outside these thin layers the Euler equation is OK. Occasionally these thin layers are not on against the boundary, as in the example of a jet, see \S\ref{sec:2djet}.

{\color{darkgray}
%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Matched Asymptotic Expansions}

This section gives the mathematical background for the Euler/Prandtl limits, providing the machinery to go beyond the leading order which we concentrated on in subsquent sections.

%%%%%%%%%%%%%%%
\subsection{Model problem}

To find $y(x;\epsilon)$ with small positive parameter $\epsilon$
\[
 \epsilon y' + y = e^{-x}, \quad x\ge 0, \qword{and} y = 0 \qword{at} x=0.
\]
The key feature is a small parameter multiplying the highest derivative. The simple model is chosen to have a simple exact solution
\[
 y = \frac{1}{1-\epsilon}\left(e^{-x} - e^{-x/\epsilon}\right).
\]
\cpicture{0.6}{MAS1}

%%%%%%%%%%%%%%%%
\subsection{OUTER}
 
The asymptotic expansion for fixed $x\ne0$ as $\epsilon\to0$ is called the ``Outer'' expansion
\[
 y \sim e^{-x} + \epsilon e^{-x} + \epsilon^2 e^{-x}.
\]
The $e^{-x/\epsilon}$ term is exponentially small, i.e.~after all powers $\epsilon^n$ in this expansion for $x$ fixed.
\cpicture{0.6}{MAS2}  

%%%%%%%%%%%%%%%%
\subsection{INNER}

The exponentially small term can be made visible by making an alternative asymptotic expansion with $\xi = x/\epsilon$ fixed as $\epsilon\to0$, in the so-called ``Inner'' expansion
\[
  y = \frac{1}{1-\epsilon}\left(e^{-\epsilon\xi} - e^{-\xi}\right) \quad
  \sim\quad (1 - e^{-\xi}) + \epsilon(1-\xi - e^{-\xi}).
\]
\cpicture{0.6}{MAS3}


\bigskip
We have obtained these two asymptotic expansions from the exact solution, which was available because the model problem was simple. Typically in real problems, like the nonlinear partial differential Navier-Stokes equation, an exact solution is not available. For such problems, can one first approximate the question in order to obtain an exact solution of an approximate question instead of the above approximation of an exact solution to the original unapproximated question.

%%%%%%%%%%%%%%%%%
\subsection{Seek the outer}

To find the outer expansion, one seeks a solution for $x\ne0$ by posing
\[
 y(x;\epsilon) \sim y_0(x) + \epsilon y_1(x).
\]
Substituting into the governing equation and comparing terms with the same powers of $\epsilon$ on the two sides of the equation, one has
\begin{alignat*}{3}
 \qword{At $O(\epsilon^0)$:}& y_0 = e^{-x}, &\qword{hence} y_0 = e^{-x},\\
 \qword{At $O(\epsilon^1)$:}& y'_0 +y_1=0, &\qword{hence} y_1 = e^{-x}.\\
\end{alignat*}
Note that the boundary condition $y=0$ at $x=0$ cannot be applied because $x=0$ is not in the outer domain.

%%%%%%%%%%%%%%%%%
\subsection{Seek the inner}

To find the inner expansion where the boundary condition is applied, one first introduces the {\em stretched variable\/} $x = \epsilon\xi$. In terms of this variable the original governing equation becomes
\[
 \frac{dy}{d\xi} + y = e^{\epsilon\xi} \quad\sim\quad 1 - \epsilon\xi
\]
One then poses the expansion
\[
 y(x;\epsilon) \sim Y_0(\xi) + \epsilon Y_1(\xi).
\]
Substituting into the new form of the governing equation and comparing terms with the same powers of $\epsilon$ on the two sides of the equation, one has
\begin{alignat*}{4}
  \qword{At $O(\epsilon^0)$:}& Y'_0+ Y_0 = 1, &\qword{with} Y_0(0)=0,
    &\qword{hence} Y_0 = 1- e^{-\xi},\\
  \qword{At $O(\epsilon^1)$:}& Y'_1 +Y_1=-\xi, &\qword{with} Y_1(0)=0,
      &\qword{hence} Y_1 = 1-\xi -e^{-\xi}.\\
\end{alignat*}

In more complicated problems, there are sometimes undetermined constants of integration which have to be found by ``matching'', i.e.~looking at 
\[
 \underset{\displaystyle x\to0}{\rm OUTER} \quad\equiv\quad
    \underset{\displaystyle \xi\to\infty}{{\rm INNER}}.
\]
}%color


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Boundary layer eqautions}
\label{sec:bleqn}

We consider the two-dimensional case.

\cpicture{0.6}{BL}

\noindent  Let the tangential velocity $u$ vary in the flow $x$-direction on a lengthscale $L$, and vary across the flow on the smaller lengthscale of the boundary layer thickenss $\delta(x)$, i.e.~we assume $\delta\ll L$. As one comes out of the boundary layer, $y/\delta\to\infty$, the tangentail velocity tends to a given free-stream flow $u(x,y,t)\to U(x,t)$, while the transverse $y$-velocity ``decays'' (to be interpretted later).

\noindent Within the boundary layer, the tangentail velocity will be $u=O(U)$, but the $y$-velocity will be much smaller. We can estimate its size from the incompressibility condition
\[
 0 = \pd{u}{x}+\pd{v}{y} \quad\to\quad \frac{U}{L} = \frac{v}{\delta}
   \quad\to\quad v = O\left(U\frac{\delta}{L}\right).
\]
In other words, $v$ is $U$ deflected through a small angle $\delta/L$.

%%%%%%%%%%%
\paragraph{First boundary layer approximation}
We negelect $\displaystyle\pd{^2}{x^2}$ compared with $\displaystyle\pd{^2}{y^2}$.  Thus in the $x$-momentum equation becomes
\[
  \rho\left(\pd{u}{t} + u\pd{u}{x} + v\pd{u}{y}\right) =
    -\pd{p}{x} + \mu\pd{^2u}{y^2}.
\]
\begin{itemize}
 \item Note we keep $\displaystyle v\pd{u}{y}$ because although $v$ is small $\displaystyle \pd{}{y}$ is large:
  \[
   v\pd{u}{y} = O\left(U\frac{\delta}{L}\frac{U}{\delta}\right) =
    O\left(u\pd{u}{x}  \right).
  \]
 \item Note the need to balance viscous and inertial terms gives the boundary layer thickness:
   \[
    \rho\frac{U^2}{L} = \mu\frac{U}{\delta^2} \quad\to\quad
     \delta = \left(\frac{\nu L}{U}\right)^{1/2} =
     L \left(\frac{UL}{\nu}\right)^{-1/2},
   \]
   so that the boundary is thin when the Reynolds number is large.
 \item Outside the boundary layer, $y/\delta\to\infty$, $u\to U(x,t)$ and the equation reduces to
   \[
      \rho\left(\pd{U}{t} + U\pd{U}{x}\right) = -\pd{p}{x}.
   \]
   This gives an estimate for the pressure variations {\em along\/} the flow
   \[
    \underset{x}{\Delta} p = O(\rho U^2)
   \]
\end{itemize}

%%%%%%%%%%%%%%
\paragraph{Second boundary layer approximation}
Consider the $y$-momentum equation
\[
  \rho\left(\pd{v}{t} + u\pd{v}{x} + v\pd{v}{y}\right) =
    -\pd{p}{y} + \mu\pd{^2v}{y^2}.
\]
This gives an estimate for the pressure variations {\em across\/} the flow
\[
 \underset{y}{\Delta} p = O\left(\rho U^2\frac{\delta^2}{L^2}\right).
\]
Hence the variation across the boundary layer is $\delta^2/L^2$ smaller than along the flow. Hence $\partial p/\partial x$ varies little across the boundary layer from its value outside. Thus we have the governing boundary layer equations
\begin{gather*}
  \boxed{\rho\left(\pd{u}{t} + u\pd{u}{x} + v\pd{u}{y}\right) =
    \rho\left(\pd{U}{t} + U\pd{U}{x}\right)  + \mu\pd{^2u}{y^2},}\\
 \boxed{\pd{u}{x} + \pd{v}{y} = 0.}
\end{gather*}

In three-dimensional flows, there are two momenta equations (quasi-Cartesian) for the two directions within the locally flat layer, plus a mass conservation equation (in proper curvilinear coordinates).

{\color{darkgray}
The missing $\displaystyle\pd{^2}{x^2}$ term makes the steady boundary layer equation parabolic in $x$ instead of elliptic, so can be solved numerically by stepping downstream.
}%color


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Boundary layer on a flat plate}
\label{sec:flatplate}

Blassius 1908. A simple application of the boundary layer equations, although an important flow.

The plate starts at $x=0$. The external flow $U$ is constant.
\cpicture{0.7}{Blassius}
With no $x$-length scale imposed in the definition of the problem, we take $L=x$ itself, so heading for a similarity solution. The estimate of the boundary layer thickness of \S\ref{sec:bleqn} gives
\[
 \delta(x) = \sqrt{\frac{\nu x}{U}},
\]
which is just $\sqrt{\nu t}$ with $t=x/U$. To satisfy the mass conservation, we introduce a streamfunction $\psi(x,y)$ with
\[
 u = \pd{\psi}{y}, \quad v= - \pd{\psi}{x}.
\]
Now $u=O(U)$ and $y=O(\delta)$, so $\psi = O(U\delta)$. Hence we seek a solution  of the (similarity) form
\[
  \psi = U\delta(x)f(\eta), \qword{where} \eta= \frac{y}{\delta(x)}.
\]
We need
\[
\left(\pd{\delta}{x}\right)_y = \frac{\delta}{2x}, \qword{and}
 \left(\pd{\eta}{x}\right)_y = - \frac{\eta}{2x}.
\]
Then
\[
 u = Uf', \quad v = -U\frac{\delta}{2x}f - U\delta f'\frac{-\eta}{2x}
     = \frac{U\delta}{2x}(-f+\eta f').
\]
The steady boundary layer equation
\[
   \rho\left(u\pd{u}{x} + v\pd{u}{y}\right) =
    \rho U\underset{=0}{\pd{U}{x}}  + \mu\pd{^2u}{y^2},\\ 
\]
becomes
\[
 \rho\left(Uf'Uf''\frac{-\eta}{2x}
           + \frac{U\delta}{2x}(-f+\eta f')Uf''\frac1\delta \right) 
  = \mu Uf'''\frac{1}{\delta^2}.
\]
With $\delta^2 = \mu x/\rho U$, this equation is satisfied at all $x$ and $y$ if
\[
 -\tfrac12 ff'' = f'''.
\]
This is an ordinary differential equation to be solved subject to
\[
 f = f' = 0 \qword{at} \eta=0, \qword{and} f'\to1 \qword{as} \eta\to\infty.
\]
The equation is solved numerically, shooting from $\eta=0$ with $f''(0)=0.332$.
\cpicture{0.7}{Blassius2}
Plotted is $f'$, i.e.~$u$. Also
\[
 f \sim \eta-1.72 \qword{as} \eta\to\infty.
\]

\bigskip
One can apply the above local solution to a finite length plate. The drag on each side (top and bottom) is per unit $z$-width
\begin{gather*}
  \int_0^l \left[\left.\mu\pd{u}{y}\right\vert_{y=0} =
    \frac{\mu U f''(0)}{\delta(x) = \sqrt{\nu x/U}} \right]\,dx \\
  = 2f''(0) \mu U \left(\frac{Ul}{\nu}\right)^{1/2}
  = \tfrac12\rho U^2l \left(\frac{Ul}{\nu}\right)^{-1/2} 4f''(0).
\end{gather*}
This is the correct leading order drag at $Re\gg1$ on a finite plate. There are corrections from the leading edge where $x=\delta(x)$, and bigger corrections from the trailing edge where a local cusp in the velocity profile produces a high acceleration and so a large pressure change.

\bigskip
Note that
\[
 v = \frac{U\delta}{2x}(-f + \eta f' \to 1.72) \qword{as} y/\delta\to\infty.
\]
As the boundary layer grows downstream, less fluid is moving at $U$, so an excess must be pushed out of vertically. One says that the boundary layer is ``detraining''. There is an effective boundary to the inviscid outer flow at $y=1.72\delta(x)$. This is called the ``displacement thickness''
\[
 \delta^* = \int_0^\infty (U-u)\,dy/U.
\]

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Two-dimensional momentum jet}
\label{sec:2djet}

Schlichting 1932. At high Reynolds numbers, one observes a jet with flow confined to a thin region in which friction acts. Thus we shall use the boundary layer equations as an approximation to the Navier-Stokes equation, even though there is no boundary.
\cpicture{0.7}{Jet}
We need to find the scalings of the thickness of the jet $\delta(x)$ and the centre line velocity $U(x)$. Balancing the inertial and viscous terms in the boundary layer equations
\[
 \frac{\rho U^2}{x} = \frac{\mu U}{\delta^2},
\]
where $x$ is used for the length-scale along the flow, as no other length is given. For a second relation between $\delta$ and $U$ we use a global result that the momentum flux (per unit $z$-width) in the jet is the same at each $x$, because there is no friction from $\infty$ and we have dropped the $\displaystyle\pd{^2}{x^2}$ terms, i.e. it is claimed that
\[
 F = \int_{-\infty}^{\infty} \rho u^2\,dy \qword{is independent of $x$.}
\]
To check
\[
 \frac{dF}{dx} = 2\int_{-\infty}^{\infty} \rho u \pd{u}{x}\,dy.
\]
Using the steady boundary layer equations
\[
= 2\int_{-\infty}^{\infty}\left(-\rho v\pd{u}{y} - \frac{dp}{dx} + \mu\pd{^2u}{y^2}\right)\,dy.
\]
Now by $u\to0$ outside the jet as $y\to\infty$ we have $dp/dx=0$. Then integrating by parts
\[
 = -2\rho vu\bigg\vert_{-\infty}^{\infty} +2\int_{-\infty}^{\infty} \rho\pd{v}{y}u\,dy
  +\mu\pd{u}{y}\bigg\vert_{-\infty}^{\infty}.
\]
Now $u\to0$ and $\partial u/\partial y\to0$ outside the jet, and by incompressibility $\partial v/\partial y = - \partial u/\partial x$. Hence
\[
 = - \frac{dF}{dx},
\]
equal and opposite to the first expression, which must therefore vanish. Hence we have a second relation between $\delta$ and $U$
\[
 F = \rho U^2\delta.
\]
Solving
\[
 \delta = \left(\frac{\rho x^2\nu^2}{F}\right)^{1/3}, \quad
 U = \left(\frac{F^2}{\rho^2 x \nu}\right)^{1/3}.
\]
As in the previous section about a flat plate, we use a streamfunction $\psi=O(U\delta)$, and seek a similarity solution (self-similar profile of the jet)
\[
 \psi(x,y) = U(x)\delta(x)f(\eta) \qword{with} \eta = \frac{y}{\delta(x)}.
\]
We need
\[
 \left(\pd{(U\delta)}{x}\right)_y = \frac{U\delta}{3x}, \quad
 \left(\pd{\eta}{x}\right)_y = -\frac{2\eta}{3x}.
\]
Then
\[
 u = \pd{\psi}{y} = Uf', \quad
 v = -\pd{\psi}{x} = -\frac{U\delta}{3x} - U\delta f'\frac{-2\eta}{3x}
   = \frac{U\delta}{x}\left(-\tfrac13f+\tfrac23\eta f'\right).
\]
The steady boundary layer equation
\[
 \rho\left(u\pd{u}{x} + v\pd{u}{y}\right) = -\frac{dp}{dx} + \mu\pd{^2u}{y^2},
\]
becomes, with $dp/dx=0$ from outside the jet,
\[
 \rho\left(Uf'\left(\frac{-U}{3x}f' + Uf''\frac{-2\eta}{3x}\right)
     +\frac{U\delta}{x}\left(-\tfrac13f+\tfrac23\eta f'\right)uf''\frac{1}{\delta}\right) = \frac{\mu U}{\delta^2}f'''.
\]
With $\delta^2=x\mu/\rho U$ this equation is satisfied for all $x$ and $y$ if
\[
 -\tfrac13f'^2 - \tfrac13ff'' = f'''.
\]
This nonlinear ordinary differential equation is to be solved subject to
\begin{gather*}
  f' \to 0 \qword{as} \eta \to \infty, \qword{no flow outside jet,}\\
  f = 0 \qword{at} \eta=0, \qword{defines centre of jet,} \\
  \int_{-\infty}^{\infty} f'^2\,dy = 1, \qword{momentum flux normalisation.}
\end{gather*}
Integrating once
\[
 -\tfrac13 ff' = f'' + \hbox{\rm const}.
\]
The constant vanishes from conditions outside the jet.
Integrating again
\[
 \tfrac16\left(k^2 - f^2\right) = f',
\]
with constant $k$. Integrating again
\[
 f = k \tanh\frac{k}{6}(\eta-\eta_0),
\]
with $\eta_0=0$ from the centre line at $\eta=0$.
Now applying the momentum flux normalisation
\[
 1 = \int_{-\infty}^{\infty} f'^2\,dy
 = \int_{-\infty}^{\infty} f'\tfrac16(k^2-f^2)\,dy
 = \tfrac16(k^2f - \tfrac13 f^3)\Big\vert_{-\infty}^{\infty}
 = \tfrac29k^3,
 \]
so $k=(9/2)^{1/3}$.

\medskip
The mass flux in the jet (per unit $z$-width) is
\[
  \int_{-\infty}^{\infty} \rho u\,dy
  = \rho U\delta \int_{-\infty}^{\infty} f'\,d\eta
  = 2k\rho U(x)\delta(x) \propto x^{1/3}.
\]
Thus the jet drags with it more and more fluid in a process of ``entrainment'', as in shower curtains.
\cpicture{0.5}{Jet2}
The Reynolds number $Re=U\delta/\nu\propto x^{1/3}$ increases, so at large $x$ the jet  becomes unstable. One can also solve the axisymmetric ``round'' jet, momentum wakes behind towed objects, momentum-free wakes behind self-propelled objects, and shear layers.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Effect of acceleration and deceleration of the external stream}
\label{sec:FS}

Falkner-Scan 1930. Consider a boundary layer on a rigid wall, starting at $x=0$, with flow outside the boundary layer
\[
 u(x) = U_0(x/\ell)^m.
\]
This flow occurs in potential (irrotational) flow around a corner of angle $\pi/(m+1)$.

\begin{minipage}{0.3\linewidth}
  $m=1$, angle $=\frac{\pi}{2}$
\end{minipage}
\begin{minipage}{0.3\linewidth}
  \cpicture{0.7}{FS1}
\end{minipage}
\begin{minipage}{0.3\linewidth}
  Stagnation flow, \S\ref{sec:stag}.
\end{minipage}

\vspace{20pt}
\begin{minipage}{0.3\linewidth}
  $m=0$, angle $=\pi$
\end{minipage}
\begin{minipage}{0.3\linewidth}
  \cpicture{0.7}{FS2}
\end{minipage}
\begin{minipage}{0.3\linewidth}
  Flat plate, uniform

  flow, \S\ref{sec:flatplate}.
\end{minipage}

\vspace{20pt}
\begin{minipage}{0.3\linewidth}
  $m=-\tfrac13$, angle = $\frac{3\pi}{2}$
\end{minipage}
\begin{minipage}{0.3\linewidth}
  \cpicture{0.7}{FS3}
\end{minipage}
\begin{minipage}{0.3\linewidth}
  Flow decelerates from

  $x=0$ to $\infty$.
\end{minipage}

In the boundary layer equations, the pressure gradient is set by the external flow
\[
 -\frac{dp}{dx} = \rho U \frac{dU}{dx} = \rho m \frac{U^2}{x}.
\]
Then balancing these inertial terms with the viscous terms
\[
 \rho \frac{U^2}{x} = \mu \frac{U}{\delta^2},
\]
gives the boundary layer thickness
\[
 \delta(x) = \sqrt{\frac{\nu x}{U(x)}} \propto x^{(1-m)/2}.
\]
We seek a similarity solution using a streamfunction of the form
\[
 \psi(x,y) = U(x) \delta(x) f(\eta) \qword{with} \eta = \frac{y}{\delta(x)}.
\]
We need
\[
 U\delta \propto x^{(m+1)/2},\quad \eta \propto x^{(m-1)/2}.
\]
Then
\[
 u = \pd{\psi}{y} = U f', \quad
 v = -\pd{\psi}{x} = -\frac{m+1}{2} \frac{U\delta}{x} f
                     - U\delta f' \frac{m-1}{2}\frac{\eta}{x}.
\]
Substituting into the boundary layer equation yields, after some cancellations
\[
 m f'^2 - \tfrac12(m+1)ff'' = m + f'''.
\]
[Check $m=1$ with \S\ref{sec:stag}, and $m=0$ with \S\ref{sec:flatplate}.] 
This equation for $f$ is to be solved subject to
\[
 f=f'=0 \qword{at} \eta=0, \quad f'\to 1 \qword{as} \eta\to\infty.
\]
This is solved numerically, shooting from $\eta=0$ with $f(0)=f'(0)=0$ and a guess for $f''(0)$, adjusting the guess until $f'\to1$ as $\eta\to\infty$.

\medskip
\begin{minipage}{0.5\linewidth}
  Results:
  \begin{align*}
    m &\quad f''(0) \\
    1 &\quad 1.233 \\
    0.5 &\quad 0.900 \\
    0 &\quad 0.322 \\
    -0.05 &\quad 0.213
  \end{align*}
  All with similar velocity profile.
\end{minipage}
\begin{minipage}{0.4\linewidth}
  \cpicture{0.4}{FS4}
\end{minipage}

\medskip
\noindent But also for small $m<0$ there is a second solutions

\begin{minipage}{0.5\linewidth}
  e.g. $m=-0.05$ with $f''(0)=-0.098$.

  with reverse flow near to the wall.
\end{minipage}
\begin{minipage}{0.4\linewidth}
  \cpicture{0.57}{FS5}
\end{minipage}

\bigskip
The solutions with reverse flow are probably unstable and impossible to set up in an experiment.
The two branches come together at $m=-0.0904$, corresponding to a wedge angle of $180^\circ + 18^\circ$. 

\noindent In $-1<m<-0.0904$, the profiles oscillate about $f'\to1$, with a range of $f''(0)$.
In $m<-1$, the oscillations do not decay as $\eta\to\infty$.
All $m<-0.0904$ are unrealistic.

The solutions are acceptable for accelerating external flows, $m>0$, and for slightly decelerating flows $-0.0904<m<0$.
If the external flow is slowing down,

\cpicture{0.5}{FS6}

\noindent mass conservation gives an advection away from the wall. With advection and diffusion away from the wall, it is not possible to confine the vorticity near to the wall, so the boundary layer structure breaks down. Typically the boundary layer separates from the wall and head off in the interior of the flow.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Flow past a body at high Reynolds number}

This section is entirely qualatitive, based on observations of experiments.

\subsection{Streamline body at zero lift}

\bigskip
\cpicture{0.7}{aerofoil1}
\begin{itemize}
 \item Solve Euler problem for the potential flow, with $\bfu\cdot\bfn=0$ but allow slip, $\bfu_{\rm tang}$, to find the tangential velocity $U(x)$ just outside the boundary layer.
 \item Solve numerically the boundary layer equationns, with external pressure gradient
   \[
    - \frac{dp}{dx} = \rho U\frac{dU}{dx},
   \]
   to find the tangential viscous stress $\mu\partial u/\partial y$ evaluated on the boundary $y=0$. Integrate this stress over the boudary to obtain the ``skin friction'' drag force (per unti $z$-width)
   \[
     O\left(\rho U^2\ell \left(\frac{U\ell}{\nu}\right)^{-1/2}\right).
   \]
\end{itemize}

%%%%%%%%%%%%%%%
\subsection{Bluff body}

A bluff body is a non-streamline body, such as a sphere or cylinder.
\bigskip
\cpicture{0.7}{Bluff1}

The potential flow past the body has a maximum and afterwards decelerates.
If one uses this potential flow as the external flow $U(x)\propto \sin\theta$ in the boundary layer equations, one finds that the solution blows up at $\theta=104.5^\circ$. One observes in experiments

\cpicture{0.9}{Bluff2}

The separation point is roughly at $\theta=80^\circ$ when laminar at $Re=10^{5-}$, but moves to $120^\circ$ at $Re=10^{5+}$ when the flow becomes turbulent. The turbulence reduces the cross-sction of the wake, which thus reduces the drag, so trip-wires are used to trigger turbulence at lower $Re$.
\[
 \qword{Drag} \sim \tfrac12\rho U_0^2 A C_D,
\]
with $A$ cross-sectionaal area of the body and $C_D$ a drag coefficient which depends on the shape, $0.3$ for a sphere.

%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Streamline body with lift}

It is necessary to keep the ``angle of attack'' less than $18^\circ$, \S\ref{sec:FS}, so that there is only slight deceleration, so that the boundary layer remains attached.
%\bigskip
\cpicture{0.6}{aerofoil2}

\noindent A clockwise circulation makes the flow faster along the top surface than along the bottom. By Bernoilli (in the potential flow) the pressure is lower on the top surface than on the bottom surface. Hence there is an upward lift force.

\bigskip
If the angle of attack is large, the aerofoil ``stalls'', and the flow is like a bluff body.
\cpicture{0.7}{aerofoil3}
The drag jumps by a factor of $Re^{1/2}$ and lift is lost.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Flow in a wedge}


\cpicture{0.3}{JHamel1}

Jeffery-Hamel 1914.
In \S\ref{sec:corner} we looked at a source flow in a wedge at low Reynolds numbers. There is an exact (similarity) solution of the Navier-Stokes equation at all Reynolds numbers with a radial flow, which we shall now derive before considering the behaviour at high Reynolds numbers.

Consider the flow in two dimensions with a volume flux (per unit $z$-width). We use a streamfunction for the radial flow
\[
 \psi = Qf(\theta), \qword{so} u_r = \frac1r\pd{\psi}{\theta} = \frac{Q}{r}f', \quad u_\theta = - \pd{\psi}{r} = 0,
\]
The vorticity is $\bfom=(0,0,\omega)$ with
\[
 \omega = - \nabla^2\psi = \frac{Q}{r^2}f''.
\]
The vorticity equation is
\[
 \underset{=0,\;\rm steady}{\pd{\bfom}{t}} + (\bfu\cdot\nabla)\bfom
 = \underset{=0,\; 2D}{(\bfom\cdot\nabla)\bfu} + \nu\nabla^2\bfom,
\]
becomes
\[
 u_r\pd{\omega}{r} = \nu\nabla^2\omega, \qword{i.e.}
  -\frac{2Q}{\nu} f'f'' = f''''+4f''.
\]
The Reynolds number is $Q/\nu$. The vorticity equation is to be solved subject to boundary conditions
\begin{align*}
  f'(\pm\alpha) = 0& \qword{no slip,}\\
  f(\alpha)-f(-\alpha) = 1& \qword{volume flux normalisation,} \\
  f(-\alpha) = 0& \qword{additive constant.}
\end{align*}
The problem can be solved numerically by shooting from $\theta=-\alpha$ with $f=f'=0$ and guesses for $f''$ and $f'''$. These guesses are adjusted until $f'(\alpha)=0$ and $f(\alpha)=1$. From the solution one can evaluate the maximum velocity, $f'(\theta)$ maximising over $\theta$. Alternatively, one can integrate once, then multiply by $f''$ to integrate a second time, then invert to produce an elliptic integral of $f'$ for $\theta$.

At low Reynolds numbers, $\vert Q/\nu\vert\ll 1$, one recovers the solutions of \S\ref{sec:corner}, with vaguely parabolic velocity profiles whether the flow is out, $Q>0$, or in, $Q<0$.

\begin{minipage}{0.45\linewidth}
  \cpicture{0.5}{JHamel2}
\end{minipage}
\begin{minipage}{0.45\linewidth}
  \cpicture{0.5}{JHamel3}
\end{minipage}
\medskip

At high Reynolds numbers with in flow $Q<0$, there is away from the walls (Euler limit) a uniform potential flow $u_r = - \vert Q\vert/2\alpha r$. Near the walls there are boundary layers. These boundary layers are governed by
\[
 f''' + \frac{Q}{\nu}\left(f'^2 - \frac{1}{4\alpha^2}\right) = 0,
\]
where we have ignored a $4f''$ term compared with $f''''$ and then integrated once. The equation can be integrated further to find near the boundary $\theta=\alpha$
\[
 f' = \frac{1}{2\alpha}\left(1 - 3\sech^2\left(\sqrt{\frac{-Q}{2\alpha\nu}}(\alpha-\theta) + \sech^{-1}\frac{1}{\sqrt{3}}\right)\right).
\]
\cpicture{0.35}{JHamel4}

At high Reynolds numbers with out flow $Q>0$, a central jet develops at moderate $Re$, with reverse flow near the walls as $Re$ increases.

\medskip
\begin{minipage}{0.45\linewidth}
  \cpicture{0.5}{JHamel5}
  \centerline{$\alpha=0.5$,
  $Q/\nu=10$}
\end{minipage}
\begin{minipage}{0.45\linewidth}
  \cpicture{0.5}{JHamel6}
  \centerline{$\alpha=0.5$,
  $Q/\nu=20$}
\end{minipage}
\medskip

\noindent There are also asymmetric solutions, many.

The radial out flow is decelerating, so it cannot confine the vorticity to be near to the walls. The type of jet depends on the details of the nozzle. The flow often does not take the assumed similarity form.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Boundary layer on a free surface}

Consider the free surface between air and water.
Because $\mu_{\rm air}\ll\mu_{\rm water}$ and $\rho_{\rm air}\ll\rho_{\rm water}$, the air exerts only a constant pressure $p_0$ (unless there are strong winds) on the water. Hence the boundary conditions on the water flow are
\begin{align*}
  \qword{Mass:}&  u_{\rm normal} = \bfu\cdot\bfn = \hbox{\rm normal velocity of moving interface,} \\
  \qword{Momentum:}& \boldsymbol{\sigma}\cdot\bfn = - p_0\bfn \quad(+\hbox{\rm surface tension effects).}
\end{align*}
Thus the rigid boundary condition $\bfu_{\rm tang}=0$ is replaced by the free boundary condition $(\boldsymbol{\sigma}\cdot\bfn)_{\rm tang}=0.$

Potential flow next to a free surface will not satify this condition, so we need a boundary layer, a new simpler type of boundary layer. One can show in general that for steady potential flow with tangential velocity $U$
\[
 \sigma_{\rm tang} = 2\mu\frac{U}{R},
\]
where $R$ is the radius of curvature of the surface.

\begin{minipage}{0.45\linewidth}
  \[
   v \sim U\frac{x}{R}
  \]
\end{minipage}
\begin{minipage}{0.45\linewidth}
 \cpicture{0.8}{FreeS}
\end{minipage}

Irrotational flow gives
\[
 \pd{u}{y} = \pd{v}{x} \qword{so} = \frac{U}{R},
 \qword{hence result.}
\]


The failure of the potential flow to satisfy the stress-free boundary condition $\sigma_{\rm tang}=0$ is corrected by a thin boundary layer. Only a small change $\Delta U$ in the velocity $U$ is required to correct the value of the tangential stress.

\cpicture{0.5}{FreeS2}

The tangential stress in the potential flow $\mu U/R$ is corrected by $+\mu\Delta U/\delta$ to go to zero on the boundary. Hence
\[
 \Delta U = O\left(U\frac{\delta}{R}\right),
\]
with boundary layer thickness $\delta = \sqrt{\nu t} = \sqrt{\nu R/U}$. Because $\Delta U \ll U$, the boundary layer equations can be linearised to leading order, and additionally the boundary layer will not separate.

\medskip
The leading order boundary layer equations yield zero drag, so one has to look to the corrections. Fortunately there is a trick which bypasses the boundary layer calculation -- one uses dissipation.
For example, for a bubble rising under gravity, the total rate of dissipation is equal to the rate of working by gravity against the drag force
\[
 \int 2\mu e_{ij}e_{ij}\,dV = - \bfF_{\rm drag}\cdot{\bf V}_{\rm rise}.
 \]
Now the strain-rate in the boundary layer is the same size as in the potential flow, while the volume of the boundary layer is much smaller than the volume of potential flow. Hence one can ignore the boundary layer in the left hand side integral and just use the potential flow to calculate the dissipation.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Rise velocity of a spherical bubble at $Re\gg1$}

To keep the bubble spherical, one needs surface tension, which means the diameter needs to be less than $5\,{\rm mm}$ in water.

First we need to find the potential flow $\bfu=\nabla\phi$ for a spherical bubble of radius $a$ rising at ${\bf V}$ through a fluid at rest.
\begin{align*}
 &\nabla^2\phi = 0 \qword{in} r\ge a, \\
 &\pd{\phi}{r} = {\bf V}\cdot\bfn = {\bf V}\cdot\bfx/a  \qword{on} r=a, \\
 &\phi\to0 \qword{as} r\to\infty.
\end{align*}
The solution is
\[
 \phi = - \frac{({\bf V}\cdot\bfx) a^3}{2r^3},
\]
so
\[
 \bfu = - {\bf V}\frac{a^3}{2r^3} + ({\bf V}\cdot\bfx)\bfx \frac{3a^3}{2r^5}.
\]

\paragraph{Direct calculation of dissipation.} Messy, with long steps skipped.

\begin{gather*}
  u_i = -V_i\frac{a^3}{2r^3} + V_nx_nx_i \frac{3a^3}{2r^5},\\
  e_{ij} = \frac{a^3}{2r^5}\left(3V_ix_j + 3V_jx_i + 3V_nx_n\delta_{ij}
                 - 15V_nx_nx_ix_j/r^2\right),\\
  e_{ij}e_{ij} = \frac{9a^6}{4r^{10}}\left(2V^2r^2+4({\bf V}\cdot\bfx)^2\right).
\end{gather*}
Integrating, using $dV=d\Omega r^2dr$, $\int\bfx\bfx\,d\Omega=(4\pi r^2/3){\bf I}$,
\[
  \int 2\mu e_{ij}e_{ij}\,dV
   = 2\mu\frac{9a^6}{4}(2V^2 + 4V^2\tfrac13)\frac{4\pi}{5a^5} = 12\pi\mu a V^2.
\]
Hence the drag force is
\[
 -12\pi\mu a {\bf V}.
\]
This is for $Re\gg1$. Recall the results for $Re\ll1$, $6\pi\mu aV$ for a rigid sphere, and $4\pi\mu a V$ for a free sphere.

\paragraph{Extra theory for less algebra}

\begin{gather*}
  u_i = \pd{\phi}{x_i},\\
  e_{ij} = \tfrac12\left(\pd{i_i}{x_j}+\pd{u_j}{x_i}\right) = \pd{^2\phi}{x_i\partial x_j}, \\
  2\mu e_{ij}e_{ij}
   = \mu\left[\pd{^2}{x_i^2}\left(\pd{\phi}{x_j}\pd{\phi}{x_j} \right)
           -2\pd{\phi}{x_j}\pd{^2}{x_i^2}\left(\pd{\phi}{x_j}\right)\right].
\end{gather*}
The last term vanishes by $\nabla^2\phi=0$, so $\partial^2_i(\partial_j\phi)=0$. And for the first term $\partial_i\partial_i$ is a divergence of a gradient, so we can use the divergence theorem for the general result in potential flow
\[
 \int_V 2\mu e_{ij}e_{ij}\,dV = \int_S \mu\pd{}{n}\left(\vert\nabla\phi\vert^2\right)\,dS.
\]

For the potential flow for a rising sphere
\[
 u^2 = \frac{a^6}{4r^6}\left(V^2 + 3({\bf V}\cdot\bfx)^2/r^2\right).
\]
Then remembering the normal out of the fluid volume is the negative radial direction, 
\[
 \pd{u^2}{n} = - \pd{u^2}{r} = \frac{6u^2}{r}.
\]
The dissipation is then
\[
 \int_{r=a}\mu\frac{6a^6}{4a^7}\left(V^2 + 3({\bf V}\cdot\bfx)^2/r^2\right)
 = \mu\frac{6}{4a}V^2(1 + 3\tfrac13)4\pi a^2
 = 12\pi\mu a V^2,
\]
again!

As an example, $a=\tfrac13\,{\rm mm}$, water, $Re\gg1$ theory gives $V=10^{-1}\,{\rm m/s}$, and $Re=30$.
Gas bubbles exapnd as they rise, so go faster, so expand faster, and so champagne explodes out of a bottle.
In the potential flow past a sphere, there is a lower pressure at the equator by Bernoulli, which makes larger bubbles become oblates. Larger ones oscillate from side to side as they rise, spin round helices, become spherical cap bubbles.

\bigskip
%%%%%%%%%%%%%%%%%%%%%
\begin{ex}[Wind over a lake]\normalfont
Wind blowing over a reservoir exerts at the water surface a
uniform tangential stress~$S$ which is normal to, and away from, a
straight side of the reservoir.
Use dimensional analysis, based both on balancing the inertial and
viscous forces in a thin boundary layer and on the imposed boundary
condition, to find order-of-magnitude estimates for the
boundary-layer thickness $\delta(x)$ and the surface velocity $U(x)$
as functions of distance $x$ from the shore.
Using the boundary-layer equations, find the ordinary differential
equation governing the non-dimensional function $f$ defined by
\[
  \psi(x,y) = U(x)\delta(x)f(\eta),\qword{where} \eta = y/\delta(x) .
\]
What are the boundary conditions on $f$?
\end{ex}

%%%%%%%%%%%%%%%%%%%%%
\begin{ex}[Wall jet]\normalfont
A steady two-dimensional jet of fluid runs along a plane rigid
wall, the fluid being at rest far from the wall.
Use the boundary-layer equations to show that the quantity
\[
 P = \int_0^\infty u(y)\left(\int_y^\infty u(y')^2\,dy'\right)\,dy
\]
is independent of the distance $x$ along the wall.
Find order-of-magnitude estimates for the boundary-layer thickness and
velocity as functions of $x$.

Show that in the analogous axisymmetric wall jet spreading out radially
the velocity varies like $r^{-3/2}$. 
\end{ex}

%%%%%%%%%%%%%%%%%%%%%
\begin{ex}[Another stress driven flow]\normalfont
A tangential stress is applied to the boundary of an incompressible viscous fluid in $y>0$, which is otherwise at rest,
\[
 \mu\frac{\partial u}{\partial y} = - Kx^2 \quad\text{at $y=0$ for $x> 0$}
\]
with $K>0$.

Find order-of-magnitude estimates for the boundary layer thickness $\delta(x)$ and tangential surface velocity $U(x)$.

Using the boundary-layer equations, find the ordinary differential equation governing the non-dimensional function $f$ defined in the streamfunction
\[
 \psi(x,y) = U(x)\delta(x)f(\eta), \quad\text{where $\eta = y/\delta(x)$.}
\]
What are the boundary conditions on $f$?

\end{ex}

%%%%%%%%%%%%%%%%%%%%%
\begin{ex}[Corner flow]\normalfont
Show that the streamfunction $\psi(r,\theta)$ for a steady
two-dimensional flow satisfies
\[
 -\frac{1}{r} \pd{(\psi,\nabla^2\psi)}{(r,\theta)} =
  \nu\nabla^4\psi .
\]
Show further that this equation admits solutions of the form
\[
 \psi = Qf(\theta) ,
\]
if $f$ satisfies
\[
 f'''' + 4f'' + {2Q\over\nu}f'f'' = 0 .
\]
\end{ex}

%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{ex}[Decay of vibration of a drop]\normalfont
Show that the rate of dissipation of mechanical energy in an
incompressible fluid is $2\mu e_{ij}e_{ij}$ per unit volume, where
$e_{ij}$ is the rate-of-strain tensor and $\mu$ is the viscosity.

A finite mass of incompressible fluid, of viscosity $\mu$ and density
$\rho$ is held in the shape of a sphere $r<a$ by surface tension.
It is set into a mode of small oscillations in which the velocity
field may be taken to have Cartesian components
\[
 u = \beta x, \quad v = -\beta y, \quad w=0 .
\]
where $\beta \propto \exp(-\epsilon t)\sin\omega t$.
Assuming that $\epsilon\ll\omega$, calculate the dissipation rate
averaged over a cycle (ignoring the slowly varying factor 
$\exp(-\epsilon t)$) and hence show that
$\epsilon = 5\mu/\rho a^2$.
You may assume that the total energy of the oscillation is twice the
kinetic energy averaged over a cycle.
Why is is permissible to ignore the details of the boundary layer near
$r=a$?
\end{ex}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\chapter{Stability of a unidirectional inviscid flow}
\label{chpt:stab}

Many steady laminar flows at $Re\gg1$ are unstable, and breakdown to turbulence which is a time-dependent flow with many spatial scales. The subject of hydrodynamical stability is huge and complicated, so here we have one simple example plus a few extensions.

One starts with a steady base flow, and then adds small perturbations. Do any of the perturbations grow or do they all decay.

Here we shall make a linear stability calculation in which we ignore quantities quadratic or higher in the small perturbations. There are some important recent developments in which strictly nonlinear instabilities have been found in linearly stable flows.

We shall consider an inviscid two-dimensional flow of a unidirectional base flow.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Kelvin-Helmholtz instability}

This problem is very like the water waves in the earlier IB course.
At one stage it was thought that this instability might explain the generation of water waves on the sea. Now we understand that it does not.

%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{The problem}

\cpicture{0.4}{KH1}

\noindent Take the base flow to be a shear layer, ignoring a thin viscous boundary layer where the discontinuity is smoothed out.

\noindent Perturb the interface between the two streams to
\[
 y = \eta(x,t).
\]
The small perturbation requires that the slope is small,
$\displaystyle\left\vert\pd{\eta}{x}\right\vert\ll1 $

Now add velocity perturbations on both sides. The inviscid flow is irrotational, so the perturbations are potential flows
\[
 \bfu = \begin{cases} (U_1,0) + \nabla\phi_1, &\qword{in} y>\eta, \\
                      (U_2,0) + \nabla\phi_2, &\qword{in} y<\eta. \end{cases}
\]
Mass conservation $\nabla\cdot\bfu=0$ gives
\begin{align}
  \nabla^2\phi_1 = 0 &\qword{in} y>\eta, \label{eq:kh1}\\
  \nabla^2\phi_2 = 0 &\qword{in} y<\eta. \label{eq:kh2}
\end{align}
The perturbations decaying at infinity gives
\begin{align}
  \phi_1 \to 0 &\qword{as} y\to+\infty, \label{eq:kh3}\\
  \phi_2 \to 0 &\qword{as} y\to-\infty. \label{eq:kh4}
\end{align}

The kinematic boundary condition (mass conservation) is that $y=\eta(x,t)$ remains a material surface
\begin{gather*}
 \frac{D}{Dt}(\eta-y) = 0, \\ \qword{i.e.}
  \pd{\eta}{t} + \left(U+\pd{\phi}{x}\right)\pd{\eta}{x} = \pd{\phi}{y},
  \qword{on} y=\eta.
\end{gather*}
The second term in the bracket is a quadratic term, so is ignored.
We now move the evaluation of the term on the right hand side from $y=\eta$ to $y=0$ using a Taylor series
\[
 \left.\pd{\phi}{y}\right\vert_{\eta} = \left.\pd{\phi}{y}\right\vert_0
 + \eta \left.\pd{^2\phi}{y^2}\right\vert_0
 + \tfrac12 \eta^2 \left.\pd{^3\phi}{y^3}\right\vert_0 + \cdots .
\]
We can ignore the second quadratic and subsequent higher order terms. Hence
\begin{align}
 \pd{\eta}{t} + U_1\pd{\eta}{x} =& \left.\pd{\phi_1}{y}\right\vert_0, \label{eq:kh5}\\
 \pd{\eta}{t} + U_2\pd{\eta}{x} =& \left.\pd{\phi_2}{y}\right\vert_0. \label{eq:kh6}
\end{align}
\begin{minipage}{0.6\linewidth}
  Note different $v=\partial\phi/\partial y$ on two sides,

  because different $U$.

  This is just mass conservation.
\end{minipage}
\begin{minipage}{0.35\linewidth}
  \cpicture{0.7}{KH2}
\end{minipage}
 
\medskip
The dynamic boundary condition is that the pressure is continuous across the interface. We find the pressure from the IB expression for pressure in a potential flow. Euler's equation for the inviscid flow is
\[
 \rho\left(\pd{\bfu}{t} + \nabla^2(\tfrac12 u^2) - \bfu\times\bfom \right)
  = - \nabla p - \nabla\Phi.
\]
where for the conservative force we take gravity $\Phi=\rho g y$. In irrotational flow, the vorticity vanishes $\bfom\equiv0$. In the time derivative we put $\bfu=\nabla\phi$. Then all the terms are a gradient, so we can integrate to
\[
 \rho\left(\pd{\phi}{t} + \tfrac12 u^2 \right) + p + \Phi = f(t),
  \qword{independent of $x$.}
\]
Linearising
\[
 u^2 = U^2 + 2U\pd{\phi}{x} + \left(\pd{\phi}{x}\right)^2,
\]
where we can ignore the last quadratic term.
Moving the evaluation from $y=\eta$ to $y=0$ using a Taylor series and ignoring nonlinear terms, the condition that pressure is continuous on $y=\eta$ is
\begin{equation}
  \rho_1\left(\left.\pd{\phi_1}{t}\right\vert_0
              + \left.U_1\pd{\phi_1}{x}\right\vert_0 +g\eta \right)
  = \rho_2\left(\left.\pd{\phi_2}{t}\right\vert_0
              + \left.U_2\pd{\phi_2}{x}\right\vert_0 +g\eta \right) + f(t). \label{eq:kh7}
\end{equation}
Equations (\ref{eq:kh1}--\ref{eq:kh7}) constitute the problem to solve in order to find whether the perturbations grow or decay. {\em We first consider the case of $\rho_1 =\rho_2$\/}.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Growth rate}

Now the governing equations involve derivatives $\displaystyle\pd{}{t}$ and $\displaystyle\pd{}{x}$ and no special value of $t$ and of $x$, unlike $\displaystyle\pd{}{y}$ which has the privileged value $y=0$. Hence we can Fourier transform in $t$ and $x$, looking for solutions proportional to
\[
 e^{ikx+\sigma t},
\]
with {\em wavenumber\/} $k$ (the wavelength is $2\pi/k$) and {\em growth-rate\/} $\sigma$. The growth-rate may be complex, in which case $Re(\sigma)$ gives the growth and $Im(\sigma)$ gives the propagation of the perturbations as waves. Fourier transforming turns the differentials to algebra. The problem reduces to finding the eigenvalue relationship
\[
 \sigma(k),
\]
called the {\em dispersion relation\/}.

Let
\[
 \eta = Ae^{ikx+\sigma t}.
\]
Then the potentials will vary similarly in $x$ and $t$. Equations (\ref{eq:kh1}) and (\ref{eq:kh2}), then give $e^{\pm ky}$ variations. The decay conditions (\ref{eq:kh2}) and (\ref{eq:kh4}) then give
\begin{align*}
  \phi_1 =  B e^{ikx+\sigma t} e^{-ky} &\qword{in} y>\eta, \\
  \phi_2 =  C e^{ikx+\sigma t} e^{+ky} &\qword{in} y<\eta.
\end{align*}
The kinematic boundary conditions (\ref{eq:kh5}) and (\ref{eq:kh6}) then relate the coefficients $B$ and $C$ to $A$
\begin{align*}
  & (\sigma + ikU_1)A = -kB, \\
  & (\sigma + ikU_2)A = kC.
\end{align*}
Finally in the dynamic boundary condition (\ref{eq:kh7}) the function $f(t)$ (independent of $x$) has no variation like $e^{ikx+\sigma t}$, so
\[
 (\sigma + ikU_1)B = (\sigma + ikU_2)C.
\]
Substituting $B$ and $C$ in terms of $A$ yields
\[
 (\sigma + ikU_1)^2 + (\sigma + ikU_2)^2  = 0.
\]
Solving we have finally the dispersion relation
\[
 \sigma(k) = -ik\frac{U_1+U_2}{2} \pm k \frac{U_1-U_2}{2}
\]
This means disturbances vary as
\[
  e^{\displaystyle ik(x-\tfrac12(U_1+U_2)t)}\,e^{\displaystyle \pm\tfrac12(U_1-U_2)kt}.
\]
The first factor shows the disturbance propagating with the mean velocity $\tfrac12(U_1+U_2)$, while the second factor shows two modes, one growing and one decaying. That there is one growing mode means that the flow is unstable. Worse, there is a problem as $k\to\infty$, which we shall return to later.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{The physical mechanism}

We switch to a frame moving with the mean velocity, in which the interface is not propagating.

\cpicture{0.7}{KH3}

The crowding of the streamlines over the peaks, increases the flow, with a corresponding decrease in the troughs. The need to accelerate to the peak and decelerate after requires a pressure minimum on the peak -- Bernoulli. Similarly there is a pressure maximum in the troughs. 

\cpicture{0.7}{KH4}

These pressures generate an acceleration of the fluid, which grow the perturbation of the interface. Hence the flow is unstable.

Actually there are two modes, one growing and one decaying. The growing one is explained by the figures above. So is the decaying mode, but more subtly. In the decaying mode, the velocity of the (left) peak is downwards. An upward acceleration of a downward velocity is a decreasing velocity, i.e.~a decaying mode.

An alternative description of the mechanism uses vorticity. The jump in the base velocity across the interface makes the interface a sheet of vorticity.

\cpicture{0.7}{KH5}

\noindent Again switch to a frame moving with the mean velocity. In this frame, the vorticity in the interface at the (left) peak is advected to the right, while the vorticity at the (right) trough is advected to the left. This leads to an accumulation of vorticity at mid-point between the peak and the trough. This vorticity spins the interface  to the left upwards and the interface to the right downwards. Hence the disturbance grows. There is the same subtly issue of using this picture to explain the decaying mode.

\cpicture{0.7}{KH5b}

In a nonlinear development, the vorticity accumulates further at the mid-points to form strong vortices. These vortices later pair up, and the vorticity within them merging to from a new larger vortex.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Spatial growth}
\label{sec:spgrowth}

The analysis in \S\ref{sec:spgrowth} was for {\em temporal growth\/}, in which one starts with a perturbation {\em sinusoidal in space\/} and asks if {\em it grows in time\/}.
However the typical experiment for a shear layer is to introduce a perturbation {\em sinusoidal in time\/} at $x=0$ and look to see if {\em it grows is space\/}.

\cpicture{0.7}{KH6}

To investigate the spatial stability, one sets $\sigma=i\omega$ with $\omega$ real and $k$ becomes complex.

\noindent We have the same governing equations, and so they result in the same dispersion relation, i.e.~changing $\sigma$ to $i\omega$
\[
 (i\omega + ikU_1)^2 + (i\omega + ikU_2)^2 = 0.
\]
Solving for $k$
\[
 k = -\omega\frac{U_1+U_2}{U_1^2+U_2^2} \pm i\omega\frac{U_1-U_2}{U_1^2+U_2^2}.
\]
This means disturbances vary as
\[
 e^{\displaystyle i\omega\left(t - x\frac{U_1+U_2}{U_1^2+U_2^2}\right)}\,\,e^{\displaystyle \mp\omega x\frac{U_1-U_2}{U_1^2+U_2^2} }.
\]
Again the second factor gives one growing and one decaying mode. Hence the flow is unstable.

Note that the propagation speed $(U_1^2+U_2^2)/(U_1+U_2)$ has a different value to the temporal analysis $\tfrac12(U_1+U_2)$.
This ``propagation'' is for the wave crests. A later course about waves will show that the disturbance energy propagates with $\partial\omega/\partial k$ rather than $\omega/k$. Usually disturbances propagate as they grow, called a ``convective instability'', so decay in time at any fixed point. Those that do not convect are called ``absolute instabilities''.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Ultra-violet divergence}

Returning to the the temporal analysis, the growth rate was found to be
\[
 \sigma_r = Re(\sigma) = \tfrac12k(U_1-U_2).
\]
This grows without bound as $k\to\infty$ (short waves). This makes the original problem ill-posed.

However for short waves, one cannot ignore the viscous boundary layer that smooths the velocity discontinuity. Hence our analysis is limited to
\[
k < \frac{U}{\nu}.
\]
A full analysis requires the smoothed velocity profile plus the viscous equations for the perturbations. One finds that the increasing growth rate levels off at $O(U^2/\nu)$.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Stabilisation of long waves by gravity}

So far we have studied the case of equal density fluids. To bring in gravity, we must make the densities different. We consider the case of a light fluid above a heavy one, $\rho_1<\rho_2$.

The temporal stability analysis gives
\[
 \rho_1(\sigma + ikU_1)^2 +  \rho_(\sigma + ikU_2)^2 + (\rho_2+\rho_1)gk=0.
\]
Solving
\[
 \sigma = \frac{-ik(\rho_1U_1+\rho_2U_2)
   \pm\sqrt{\rho_1\rho_2(U_1-U_2)^2k^2 - (\rho_1+\rho_2)(\rho_2-\rho_1)gk}}
        {\rho_1+\rho_2}.
\]
Hence for sufficiently long waves
\[
 k < \frac{(\rho_1+\rho_2)(\rho_2-\rho_1)}{\rho_1\rho_2}\frac{g}{(U_1-U_2)^2},
\]
both roots are purely imaginary and there is no growth or decay.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Stabilisation by gravity (for long waves) and by surface tension (for short waves)}
For two immiscible fluids, there will be an interfacial tension (surface energy). which generates a jump in the pressure across the interface (dynamic boundary condition) equal to the value of the surface tension, $\gamma$, multiplied by the curvature of the surface, which in the small slope limit is $\partial^2\eta/\partial x^2$, i.e.
\[
  \rho_1\left(\left.\pd{\phi_1}{t}\right\vert_0
              + \left.U_1\pd{\phi_1}{x}\right\vert_0 +g\eta \right)
  = \rho_2\left(\left.\pd{\phi_2}{t}\right\vert_0
              + \left.U_2\pd{\phi_2}{x}\right\vert_0 +g\eta \right)
              + \gamma\pd{^2\eta}{x^2}.
\]
This leads to a dispersion relation
\[
 \sigma = \frac{-ik(\rho_1U_1+\rho_2U_2)
   \pm k\sqrt{\rho_1\rho_2(U_1-U_2)^2
              - (\rho_1+\rho_2)\left[\frac{(\rho_2-\rho_1)g}{k} + \gamma k\right]}}
        {\rho_1+\rho_2}.
\]
The square bracket has a minimum value of $\sqrt{(\rho_2-\rho_1)g\gamma}$ at $k^2=(\rho_2-\rho_1)g/\gamma$.
Hence $\sigma$ will be purely imaginary (no growing mode) for all wavenumbers $k$ if the velocity difference is smaller than a critical value
\[
 (U_1-U_2)^2 < \frac{(\rho_1+\rho_2)\sqrt{(\rho_2-\rho_1)g\gamma}}{\rho_1\rho_2},
\]
i.e.~the flow is stable if $\vert U_1 - U_2\vert$ is below this critical value.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Other instabilities}

This chapter has consider the stability of just one flow. There are many others which have been studied. A quick extension of our analysis of the Kelvin-Helmholtz instability is to the Raleigh-Taylor instability, in which there is no base flow, $U_1=0=U_2$, and there is heavy fluid above light, $\rho_1>\rho_2$.

One can study the Raleigh-Plateau instability of a jet in which surface tension tries to reduce the surface area/energy of a cylinder by producing a line of spherical drops.

In the Hele-Shaw geometry, if a low viscosity fluid pushes out a high viscosity fluid, perturbations to a plane interface develop a viscous fingering instability which enables the mobile low viscosity fluid to bypass the sluggish high viscosity fluid.

Flow in the Couette apparatus of an annular gap between two coaxial cylinders is
unstable through a centrifugal instability if the inner cylinder rotates faster than the outer one.

Poiseuille flow in a pipe is unstable at high Reynolds numbers. This instability needs viscosity, although small. The instability is strictly nonlinear.

There are many other instabilities. 
Adding new physics brings new instabilities. Thus allowing the density to vary with temperature and including gravity, one has the possibility of hot fluid rising and cold sinking. Thus if one heats sufficiently strongly the bottom flat horizontal plate below a layer of fluid, heat is transferred across the layer more efficiently by flow  compared with molecular conduction.

\bigskip
%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{ex}[Jet next to a wall]\normalfont
A vortex sheet of strength $U$ is located at a distance $h$
above a rigid wall $y=0$ and is parallel to it, so that the fluid
velocity $(u,0,0)$ is
\[
 u = \begin{cases} U & \qword{in} 0<y<h, \\
            0 & \qword{in} y>h .\end{cases}
\]
Suppose now that the sheet is perturbed slightly to the position 
$y=h + \eta_0e^{ik(x-ct)}$ where $k>0$ is real but $c$ may be complex.
Show that
\[
 c = U/(1\pm i\sqrt{\tanh kh}) .
\]
Deduce that
\begin{itemize}
 \item the sheet is unstable to disturbances of all wavelengths;
 \item for short waves ($kh\gg1$) the growth rate $k{\rm Im}(c)$
is  $\tfrac12Uk$ and the wave propagation speed ${\rm Re}(c)$ is
$\tfrac12U$, as if the wall were absent;
 \item for long waves ($kh\ll1$) the growth rate is 
$Uk\sqrt{kh}$ (so that the wall inhibits the growth of long waves) and
the propagation speed is $U$.
\end{itemize}
\end{ex}

%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{ex}[Free jet]\normalfont
A two-dimensional jet in the $x$-direction has velocity profile
\[  
    u =\begin{cases}0 & \qword{in} y>h, \\
              U & \qword{in} -h<y<h, \\
              0 & \qword{in} y<-h .\end{cases}
\]
The vortex sheets at $y=\pm h$ are perturbed to
\[
 y = \begin{cases}+h + \eta_1 e^{ik(x-ct)},& \\
               -h + \eta_2e^{ik(x-ct)} .&\end{cases}
\]
Show that the jet is unstable to a `varicose' instability for which
$\eta_1=-\eta_2$ (identical to that of question 5), and also to a
`sinuous' instability for which $\eta_1=\eta_2$ and
\[
 c = U/(1\pm i\sqrt{\coth kh}) .
\]
[The growth rates at small $kh$ are again $Uk\sqrt{kh}$. Hence thin
jets (e.g.~smoke filaments) can suffer rather slowly growing sinuous
instabilities.]
\end{ex}


%\backmatter
%\printindex

\end{document}
