4A10 flow instability

Defining pertubations for flow instability

What is stability

• Every single frequecny and mode of pertubation that we input needs to be stable for stability.

• can characterise this with a laplace transform and teh value of s

Marginal stability

WE can also have marginal stability situations

• stationary basically standing wave solutions

• oscillatory rhyrthms

Methods of analysing stability.

We have two ways of analysing stability:

Linerarised equations

• looking at small pertubations

• solve for the growth rate

Energy analysis

• see if energy is released or needed for a given pertubation

Break up of liquid jet: Setup and background

Here we’re trying to find the critical reynolds number for a liquid jet:

• the physical mechanism is surface tension resulting in varying energy differences.

  • There is a squeeze from a smaller curvature increasing surface tension, balanced out by length wise curvature to some extent

• we will work out an energy argument based on surface tension

• Using conservation of volume to link how our curvatures in both directions

To work out growth rate:

• We will use a different method of using dimensional analysis to figure out the form of our growth rate equation
  

Break up of liquid jet: (note on surface tension)

Our surface tension can be defined in two ways:

• Tension (N/m) basically a force per unit length, a pull on the fluid

• Energy (J/m²) as a surface energy per unit area

To work out the pressure different due to surface tension

• This is identical to a hoop stress derivation

  • So ΔP = γ/R for a single curvature (γ = ΔP/R)

  • and ΔP = 2γ/R for spherical (γ = ΔP/2R)

  • And for compound curvature. $\Delta p = \gamma \left( \frac{1}{R_1} + \frac{1}{R_2} \right)$

We are just balancing and integrating our pressure force (can take as a sector) with our surface tension at an angle. The trig terms will cancel out.

Break up of liquid jet: applying pertubation and considering volume

We will apply a small amplitude harmonic disturbance

• $$r = \alpha + \beta cos kx$$

Working out how α and β vary, given that we know the overall volume is constant

• new volume is $\int_0^{\lambda}\pi r^2\,dx=\int_0^{\lambda}\pi(\alpha+\beta\cos kx)^2\,dx=\pi\lambda\left(\alpha^2+\frac{1}{2}\beta^2\right)$

so:

$$\underbrace{\pi \lambda \left( \alpha^{2} + \frac{1}{2} \beta^{2} \right)}_{\text{Volume of perturbed jet}} = \underbrace{\pi a^{2} \times \lambda}_{\text{Undisturbed volume}}$$

$$\Rightarrow \alpha = a \left( 1 - \frac{1}{2} \left( \frac{\beta}{a} \right)^{2} \right)^{1/2}$$

small amplitude:

β is much smaller the α, so can rewrite this as a binomial expansion

$$\alpha = a \left( 1 - \frac{1}{4} \left( \frac{\beta}{a} \right)^{2} \right)$$

Break up of liquid jet: considering surface area

Our surface tension is an energy per unit area

• To work out stability, we basically want to see if our pertubation increase or decreases the area

• this will tell us if the energy increase or decreases

Using the constraint from before, we can use this to work out the surface area. The result is shown in the picture

Detailed maths

Transform coordinate system along length to dx

$$= \int_{x=0}^{\lambda} 2\pi r \, ds \quad , \quad ds = \sqrt{1 + \left( \frac{dr}{dx} \right)^2} \cdot dx$$

now subbing in our expression for r both in 2πr and dr/dx

$ds=dx\left(1+(\beta k\sin kx)^2\right)^{1/2}\quad$ taking an expansion as β is really small:

$$= dx \left( 1 + \frac{1}{2} \beta^2 k^2 \sin^2 kx + O(\beta k)^4 \right)$$

Now subbing in the integral

$$\text{S.A.} = 2\pi \int_{0}^{\lambda} (\alpha + \beta \cos kx) \left( 1 + \frac{1}{2} \beta^2 k^2 \sin^2 kx \right) dx$$

quite easy to solve, basically in terms of chain rule if we get into the form sin(cos)^n or cos(sin)^n

this results in:

$$= 2\pi \lambda \left( \alpha + \frac{1}{4} \alpha \beta^2 k^2 \right)$$

now subbing in our expression for α

$\alpha = a \left( 1 - \frac{1}{4} \left( \frac{\beta}{a} \right)^2 + o \left( \frac{\beta}{a} \right)^4 \right)$

$$= 2\pi \lambda a \left( 1 + \frac{\beta^2}{4} \left( k^2 - \frac{1}{a^2} \right) \right)$$

$$= 2\pi \lambda a \left( 1 + \frac{\beta^2}{4a^2} (k^2 a^2 - 1) \right) \quad \text{if } \beta \text{ small}$$

Break up of liquid jet: stability result

Stability criteria can be seen in the image.

• With small β we can basically consider this as a balance of curvature

  • At long wavelengths we get an increase in curvature from the radius change (aximuthal) resulting in instability

  • But this is counteracted by some negative curvature from the axial curvature

  • at the balance these basically cancel each other out.

Break up of liquid jet: growth rate.

Here we are considering dimensional analysis: If we consider a perturbation in the form:
$r = a + e^{st} \beta \cos kx$

Our growth rate s has dimensions 1/T, this depends on

• wavenumber (1/L)

• original radius (a)

• density (M/L³)

• surface tension (M/T²)

we get $s \left( \frac{\rho a^3}{\gamma} \right)^{1/2} = f(ka)$

real variation

our actual variation in this non dimensionalised form can be seen below.

Break up of jet : limitations in analysis

effect of gravity

we have ignored gravity so far, we can examine this effect by comparing the hydrostatic pressure difference to the surface tension pressure difference

• Hydrostatic pressure difference: $\Delta p \sim \rho g a$

• Surface tension pressure diffference: $\Delta p \sim \gamma / a$

So our metric for whether this is significant is by dividing these 2

$$\rho g a \ll \frac{\gamma}{a} \implies a \ll \sqrt{\frac{\gamma}{\rho g}}$$

cannot predict drop size bc this is large disturbance

Basic recipe for linear analysis.

we have 6 basic steps for linear stability analysis.

1: Identify steady flow

• We want a solution for our unperturbed flow, generally want:

• Velocity $U(x)$, Pressure $P(x)$, and Interface position $\eta$

2: Identify governing equations/boundary counditions

• Governing equations like navier stokes and continuity, boundary condition like no slip, pressure at boundary etc.

3: Introduce small pertubations about the flow

• Add small pertubations to every variable:

• $$u(x, t) = U(x) + u'(x, t)$$

  $$p(x, t) = P(x) + p'(x, t)$$

  $$\eta(x, t) = a(x) + \eta'(x, t)$$

4: linearise governing equations

• Express our governing equations in terms of these pertubations

• Remember our pertubations are small so we can neglect their products (similar to boundary layer equations)

• allows us to remove many terms

5: solve linearised equations of motion

• Generally we want normal mode exponential solutions

• eg: $u'(x, z, t) \propto \hat{u}(z) e^{ikx + st}$ (shape term and then wave evolution term)

• This is a separation of variables approach

6: use boundary conditions and find S

• If s is on RHP then it is unstable

Why this is useful

• any disturbance can be expressed as a fourier series

• can look at the frequency component sensitivity

• we need every frequency to be stable for stability

Typical form of linear stability modal solutions

Linear stability analysis for plane jet: setup and introduction

Here we are considering a plane (not circular jet) so degrees of freedom (x,z,t)

We will following our linear analysis recipe to solve this:

• We will look at the interface position η, pressure p and velocity U

• Our flow is governed by continuity and the euler equation (inviscid navier stokes)

• We have a pressure boundary condition on the surface, dependent on surface tension, as well as a kinematic boundary condition (fluid particles on surface stay on surface)

• Goal here is to try eliminate variables

  • We will combine euler and continuity together to a single governing euqation

  • We get a second equation by subbing into the boundary conditions and combining togther

  • combine equations together to get a final PDE

• Now solve PDE with separation of variables with a Modal solution

• Use our boundary condition that we know this is solve for Z=a

• work out growth rate.

Linear stability analysis for plane jet: Governing equations and boundary conditions

We have two governing equations for this scenario:

• $\nabla \cdot \mathbf{u} = 0$ (mass conservation)

• $\frac{\partial \mathbf{u}}{\partial t} + (\mathbf{u} \cdot \nabla)\mathbf{u} = -\frac{1}{\rho}\nabla p$ (inviscid momentum equation)

We also have two boundary conditions:

Pressure boundary condition

• Our pressure on interface equals atmospheric pressure + a surface tension term

This is $p - p_{\infty} = \gamma \nabla \cdot \mathbf{\hat{n}} \quad \text{on } z = \eta(x,t)$

$\mathbf{\hat{n}} = \frac{\nabla F}{|\nabla F|}$, where $F = 0$ is the equation defining the surface.

Kinematic boundary condition

this is basically stating particles on the surface must remain on the surface, basically our material derivative on the surface must be zero.

Our surface function is $F(x, z, t) = z - \eta(x, t) = 0$ .

So: $\frac{DF}{Dt} = \left[ \frac{\partial F}{\partial t} + (\mathbf{u} \cdot \nabla)F \right] = 0$

• $$\frac{\partial F}{\partial t} = -\frac{\partial \eta}{\partial t}$$

• $$(\mathbf{u} \cdot \nabla)F = u\frac{\partial F}{\partial x} + v\frac{\partial F}{\partial y} + w\frac{\partial F}{\partial z} = -u\frac{\partial \eta}{\partial x} + 0 + w(1)$$

so we get $\frac{\partial\eta}{\partial t}+u\frac{\partial\eta}{\partial x}=w\text{ on }z=\eta(x,t)$ (w is the vertical velocity)

Linear stability analysis for plane jet: Adding pertubations and finding governing PDE

Basically our goal is to sub in all our values and then simplify so we get a single PDE

Introducing perturbations:

We are adding perturbations in all three variables (base velocity of zero from reference frame change)

• Velocity: $\mathbf{u} = 0 + \mathbf{u'}(x, z, t) \implies \mathbf{u} = \mathbf{u'}$

• Pressure: $p = P + p'(x, z, t) \implies p = p_{\infty} + p'$

• Interface: $\eta = a + \eta'(x, t)$

Overview of equations

Our governing equations are:

• $\nabla \cdot \mathbf{u} = 0$ (mass conservation)

• $\frac{\partial \mathbf{u}}{\partial t} + (\mathbf{u} \cdot \nabla)\mathbf{u} = -\frac{1}{\rho}\nabla p$ (inviscid momentum equation)

And boundary conditions are:

• $$p - p_{\infty} = \gamma \nabla \cdot \mathbf{\hat{n}} \quad \text{on } z = \eta(x,t)$$

• $$\frac{\partial\eta}{\partial t}+u\frac{\partial\eta}{\partial x}=\frac{\partial z}{\partial t}$$

subbing in to remove governing equations

• Continuity: $\nabla \cdot \mathbf{u'} = 0$

• Euler: $\frac{\partial \mathbf{u'}}{\partial t} + \underbrace{(\mathbf{u'} \cdot \nabla)\mathbf{u'}}_{\text{product of small terms}} = -\frac{1}{\rho}\nabla (P + p')$

  • $\rho \frac{\partial \mathbf{u'}}{\partial t} = -\nabla p'$ (remember this is a VECTOR equation)

Combining both equations by taking the divergence of both sides, we get a laplacian equation:

$$\nabla^2 p' = \frac{\partial^2 p'}{\partial x^2} + \frac{\partial^2 p'}{\partial z^2} = 0$$

subbing into the pressure boundary condition

$$p - p_{\infty} = \gamma \nabla \cdot \mathbf{\hat{n}}$$

finding unit normal

$\mathbf{\hat{n}} = \frac{\nabla F}{|\nabla F|} = \frac{-\frac{\partial \eta'}{\partial x} \mathbf{i} + 0\mathbf{j} + \mathbf{k}}{\sqrt{(\frac{\partial \eta'}{\partial x})^2 + (0)^2 + (1)^2}} \approx -\frac{\partial \eta'}{\partial x} \mathbf{i} + 0\mathbf{j} + \mathbf{k}$ (neglect squared term at bottom)

• This is basically just our gradient as maybe expected, as our vector does not need to be normalised due to small movements

our pressure is

$$\frac{\partial \eta'}{\partial t} = w' \quad \text{on } z = a$$

$$p' = -\gamma \frac{\partial^2 \eta'}{\partial x^2} \quad \text{on } z = a$$

subbing into kinematic boundary condition

$$\frac{\partial \eta}{\partial t} + u \frac{\partial \eta}{\partial x} = w \quad \text{on } z = \eta(x,t)$$

we get:

$\frac{\partial \eta'}{\partial t} = w' \quad \text{on } z = a$ (once we neglect small terms)

combining boundary conditions

We can combine to eliminate η

differentiating pressure twice to get $\frac{\partial\eta'}{\partial t}$

$\frac{\partial^2 p'}{\partial t^2} = -\gamma \frac{\partial^2}{\partial x^2} \left( \frac{\partial^2 \eta'}{\partial t^2} \right)$ now we sub in our z velocity

$\frac{\partial^2 p'}{\partial t^2} = -\gamma \frac{\partial^2}{\partial x^2} \left( \frac{\partial w'}{\partial t} \right)$ now we use the Z component of our combined euler equation

• $\rho \frac{\partial \mathbf{u'}}{\partial t} = -\nabla p'$ , Z component is $\rho\frac{\partial w'}{\partial t}=-\frac{\partial p'}{\partial z}$

combined we get

$$\frac{\partial^2 p'}{\partial t^2} = -\gamma \frac{\partial^2}{\partial x^2} \left( -\frac{1}{\rho} \frac{\partial p'}{\partial z} \right)$$

boundary equation PDE

rearranging this our final PDE is:

$$\rho \frac{\partial^2 p'}{\partial t^2} = \gamma \frac{\partial^2}{\partial x^2} \left( \frac{\partial p'}{\partial z} \right) \quad \text{on } z = a$$

Linear stability analysis for plane jet: Finding solutions for PDE

We’re looking for normal mode solutions, (basically doing separation of variables)

$$p'(x, z, t) = \hat{p}(z)e^{ikx + st}$$

Now plugging this into the governing equation first (laplace equation):

$\nabla^2 p' = 0$.

$\frac{d^2}{dz^2}\hat{p}(z) - k^2 \hat{p}(z) = 0$ this reduces it down to an ODE

this has a hyperbolic solution:

$$\hat{p}(z) = Ae^{-kz} + Be^{+kz} \quad \text{for constants } A \text{ and } B$$

• We clearly have symmetry about Z, so the constants must be the same

• So our perturbation is in the form: $p' = A(e^{-kz} + e^{+kz})e^{ikx + st}$

to find growth rate we sub into boundary condition

$$\rho \frac{\partial^2 p'}{\partial t^2} = \gamma \frac{\partial^2}{\partial x^2} \left( \frac{\partial p'}{\partial z} \right) \quad \text{on } z = a$$

• $$\rho s^2 A(e^{ka} + e^{-ka})e^{ikx + st} \Big|_{z=a} = \gamma \frac{\partial^2}{\partial x^2} \left[ kA(e^{kz} - e^{-kz})e^{ikx + st} \right] \Big|_{z=a}$$

• $\rho s^2 A(e^{ka} + e^{-ka}) = \gamma (-k^2) [kA(e^{ka} - e^{-ka})]$ can eliminate A and this is a tanh equation

• $$s^2 = -\frac{\gamma}{\rho a^3} (ka)^3 \tanh(ka)$$

verdict

Tanh is positive for all positive values of ka

• so s² = -ve

• so s must be purely imaginary neutrally stable

• makes sense if we just think about the energies for 2 seconds, clearly our areas stay the same

Result from dimensional analysis

$$s \left( \frac{\rho a^{3}}{\gamma} \right)^{1/2} = f(ka) \quad \text{or} \quad s^{2} \left( \frac{\rho a^{3}}{\gamma} \right) = F(ka)$$

matches form of $s^2 = -\frac{\gamma}{\rho a^3} (ka)^3 \tanh(ka)$

Linear stability analysis for circular jet: Growth rate

As we know for energy analysis this is unstable

Linear analysis tells us our growth rate (wiith very complicated maths)

$$s^2 = \frac{\gamma}{\rho a^3} \frac{I_n'(\alpha)}{I_n(\alpha)} \alpha (1 - \alpha^2 - n^2) \quad \text{where } \alpha = ka$$

Our $I_n$ are modified bessel functions of the first kind, basically like our hyperbolic functions but for radially varying systems.

Linear stability analysis for circular jet: Modes

• Only unstable for long wave axisymmetric modes

• Stable for azimuthal modes, strongly varying capillary force restore interface position.

Axisymmetric n = 0

$s^2 = \frac{\gamma}{\rho a^3} \frac{I_n'(\alpha)}{I_n(\alpha)} \alpha (1 - \alpha^2 - n^2) \quad \text{where } \alpha = ka$ $(1-\alpha^2-n^2)$

then (1-α²) is positive if α (ka) < 1

Azimuthal n ≠ 0

We have a boundary condition on n, so it must take integer values due to the rotation symmetry

immediately this results in $(1-\alpha^2-n^2)$ being negative

this results in a stable jet.

We can see these modes on the right

Linear stability analysis for circular jet: Development of pertubations

Linear stability analysis for circular jet: Droplet size predictions

Although linear analysis isn’t really meant for this, because it considers small pertubations, we can perform a surprisingly accurate estimation for the droplet size.

• Basically we take maximum growth wavelength (ka =0.697)

• see the volume contained within that wavelength

• this gives us d = 3.8a

Taylor couette flow: overview

Two rotating cylinders with a “narrow gap”

• This results in the formation of cells

• This is an exchange of stability, as we go from one stable mode to another stable mode.

steps for solving

• We will first look at the steady state

Taylor courette flow: Solving steady state base flow

Base flow

Here we are assuming a purely circular flow, so the flow is only in the Uθ direction, and this must only vary the radius.
Our flow structure is thus:

$\mathbf{u}(r) = \underbrace{0}_{\text{radial}} \mathbf{e}_r + \underbrace{U_{\theta}(r)}_{\text{circumferential}} \mathbf{e}_{\theta} + \underbrace{0}_{\text{vertical}} \mathbf{e}_z$

Solving for base flow structure

Our base flow is governed by two governing equations, navier stokes and continuity, writing these out in cylindrical coordinates we get:

Navier stokes

• Radial: $\frac{\partial u_r}{\partial t} + (\mathbf{u} \cdot \nabla)u_r - \frac{u_{\theta}^2}{r} = -\frac{1}{\rho}\frac{\partial p}{\partial r} + \nu \left[ \nabla^2 u_r - \frac{u_r}{r^2} - \frac{2}{r^2}\frac{\partial u_{\theta}}{\partial \theta} \right]$

• Circumferential: $\frac{\partial u_{\theta}}{\partial t} + (\mathbf{u} \cdot \nabla)u_{\theta} + \frac{u_r u_{\theta}}{r} = -\frac{1}{\rho r}\frac{\partial p}{\partial \theta} + \nu \left[ \nabla^2 u_{\theta} + \frac{2}{r^2}\frac{\partial u_r}{\partial \theta} - \frac{u_{\theta}}{r^2} \right]$

• Vertical: $\frac{\partial u_z}{\partial t}+(\mathbf{u}\cdot\nabla)u_{z}=-\frac{1}{\rho}\frac{\partial p}{\partial z}+\nu\left[\nabla^2u_{z}\right]$

Continuity this is the vector equation:

• $$\nabla \cdot \mathbf{u} = \frac{1}{r}\frac{\partial}{\partial r}(r u_r) + \frac{1}{r}\frac{\partial}{\partial \theta}(u_{\theta}) + \frac{\partial}{\partial z}(u_z) = 0$$

Plugging in our solutions to simplify

We can simplify given that U is only in the circumferential direction, and only varies with θ.

• Continuity: Automatically satisfied by our flow assumptions, every term is zero

• Radial: $-\frac{u_{\theta}^2}{r}=-\frac{1}{\rho}\frac{\partial p}{\partial r}$

• Circumferential: $\frac{\partial u_{\theta}}{\partial t} = \nu \left[ \frac{\partial^2 u_{\theta}}{\partial r^2} + \frac{1}{r} \frac{\partial u_{\theta}}{\partial r} - \frac{u_{\theta}}{r^2} \right]$

  • (splitting laplacian), pressure must be constant since its periodic

• Vertical: $0=-\frac{1}{\rho}\frac{\partial p}{\partial z}$ No vertical pressure gradient

Now for steady flow

$$\frac{d^2u_{\theta}}{d r^2}+\frac{1}{r}\frac{d u_{\theta}}{d r}-\frac{u_{\theta}}{r^2}=0$$

Solving ODE and adding boundary conditions

Use the substitution $r = e^x$ for Euler-Cauchy equation so $x=ln\left(r\right)$

For first derivative:

• $$\frac{dU_{\theta}}{dr} = \frac{dU_{\theta}}{dx} \frac{dx}{dr} = \frac{dU_{\theta}}{dx} \left( \frac{1}{r} \right)$$

• $$\implies r \frac{dU_{\theta}}{dr} = \frac{dU_{\theta}}{dx}$$

For second derivative:

• $\frac{d^2 U_{\theta}}{dr^2}=\frac{d}{dr}\left(\frac{dU_{\theta}}{dr}\right)=\frac{d}{dr}\left(\frac{1}{r}\frac{dU_{\theta}}{dx}\right)$ (now use product rule)

• $$\frac{d^2 U_{\theta}}{dr^2}=\frac{d}{dr}\left(\frac{1}{r}\frac{dU_{\theta}}{dx}\right)=\frac{1}{r^2}\left(\frac{d^2 U_{\theta}}{dx^2}-\frac{dU_{\theta}}{dx}\right)$$

Plugging back together we get:

• $$\frac{1}{r^2}\left(\left(\frac{d^2 U_{\theta}}{dx^2}-\frac{dU_{\theta}}{dx}\right)+\left(\frac{dU_{\theta}}{dx}\right)-U_{\theta}\right)=0$$

• $$\frac{d^2 U_{\theta}}{dx^2} - U_{\theta} = 0$$

result

This is solved with:

$$U_{\theta} = Ae^x + Be^{-x}$$

Giving us a solution of:

$$U_{\theta}(r) = Ar + \frac{B}{r}$$

• we just apply no slip condition to give A and B

• $$A = \frac{\Omega_2 r_2^2 - \Omega_1 r_1^2}{r_2^2 - r_1^2}$$

• $$B = \frac{(\Omega_1 - \Omega_2)r_1^2 r_2^2}{r_2^2 - r_1^2}$$

Taylor courette flow: Physical reason for instability:

Above a critical Taylor number (1708) our system becomes unstable. This is defined as: T = \frac{2(\Omega_1 r_1^2 - \Omega_2 r_2^2)\overline{\Omega}d^3}{\nu^2 r_1} > T_{\text{crit}} \approx 1708

• We can see how on the graph, our growth rate curve increases with increasing taylor number, and how we only get marginal growth at the critical taylor number.

• at even higher taylor numbers, our flow changes wave number

Physical explanation

The physical explanation here is pretty much identical to a ball on a spring which is spun around.

Angular momentum

• If we consider a toroid of fluid with radius r₁ and velocity U₁ , it’s angular momentum is r₁U₁

• Now if this is displaced slightly it must conserve angular momentum, with radius r₂ and U₂

The velocity decreases with: $U_2^{\prime}=\left(\frac{r_{1}}{r_{2}}\right)U_1$

• The velocity decrease but we don’t know if this is enough to stay within our centripetal force boudary

Centripetal force:

Now at equilibrium our force holding the fluid in place is centripetal force.

$-\frac{1}{\rho} \frac{\partial p}{\partial r} = -\frac{U_{2}^{2}}{r_{2}} \text{}$ at r₂

If we perturb our fluid from U₁ , and the resulting U₂’ is LOWER than U₂, then our system is stable.

Finding stability criteria

At our boundary, our perturbation will result in U₂’ = U₂

So setting this equal we get:

$$\left( \frac{r_{1}}{r_{2}} U_{1} \right)^{2} \le U_{2}^{2} \implies (r_{1}U_{1})^{2} \le (r_{2}U_{2})^{2} \text{}$$

$$(r_{2}^{2}\Omega_{2})^{2} - (r_{1}^{2}\Omega_{1})^{2} \ge 0 \text{}$$

• this means our angular momentum must INCREASE with RADIUS for stability

• this is an inviscid prediction which matches well at high speeds.

Taylor courette flow: Linear stability analysis

Here we’re doing our same 6 step process:

we want to find the marginal stability state for when the flow regime switches

Setting up system and plugging in pertubations

• Same with our base equations of navier stokes and continuity

Our perturbations are:

• $$p = P(r) + p'(r, z, t)$$

• $$u_{\theta} = U_{\theta} + u_{\theta}'(r, z, t)$$

• $$u_{z}=0+u_{z}^{\prime}(r,z,t)$$

• $$u_r = 0 + u_r'(r, z, t)$$

Plugging this into our cylindrical navier stokes and continuity, eliminating products of small terms, we get:

• Radial: $\frac{\partial u_r'}{\partial t} - 2\frac{U_{\theta} u_{\theta}'}{r} = -\frac{1}{\rho}\frac{\partial p'}{\partial r} + \nu\left(\nabla^2 u_r' - \frac{u_r'}{r^2}\right)$

• Circumferential: $\frac{\partial u_{\theta}'}{\partial t} + u_r' \frac{dU_{\theta}}{dr} + \frac{u_r' U_{\theta}}{r} = \nu\left(\nabla^2 u_{\theta}' - \frac{u_{\theta}'}{r^2}\right)$

• Vertical: $\frac{\partial u_z'}{\partial t} = -\frac{1}{\rho}\frac{\partial p'}{\partial z} + \nu\left(\nabla^2 u_z'\right)$

• Continuity: $\frac{u_r'}{r} + \frac{\partial u_r'}{\partial r} + \frac{\partial u_z'}{\partial z} = 0$

Narrow gap approximation

we can further simplify our equations with a narrow gap approximation.

• Basically curvature is identical on both faces

• any gradients which are order X/d so any derivatives with respect to r generally, much bigger than say a u/r term

Our equations are now: (also moving our operators together)

• Radial: $\left(\frac{\partial}{\partial t}-\nu\tilde{\nabla}^2\right)u_{r}^{\prime}-\left(\frac{2 U_\theta}{r}\right)u_{\theta}^{\prime}=-\frac{1}{\rho}\frac{\partial p'}{\partial r}\text{ }$

• Circumferential: $\left(\frac{\partial}{\partial t}-\nu\tilde{\nabla}^2\right)u_{\theta}^{\prime}+(2A)u_{r}^{\prime}=0\text{ }$

  • A is our base flow term, since we $u_{r}^{\prime}\left(\frac{dU_{\theta}}{dr}+\frac{U_{\theta}}{r}\right)$

  • So can plug in $\frac{dU_{\theta}}{dr} = A - \frac{B}{r^2}$ and $\frac{U_{\theta}}{r} = A + \frac{B}{r^2}$

• Vertical:

  • $$\left(\frac{\partial}{\partial t}-\nu\tilde{\nabla}^2\right)u_{z}^{\prime}=-\frac{1}{\rho}\frac{\partial p'}{\partial z}$$

Finding our PDE

now our goal from our narrow gap equations

• Remove pressure term in radial equation using vertical equation

  • $$\frac{\partial^2}{\partial z^2}(\text{Radial Eq. 1}) - \frac{\partial}{\partial z}\frac{\partial}{\partial r}(\text{Vertical Eq. 2})$$

• use continuity to remove the the Uz term

This gives us two PDEs in just $u_r$ and $u_\theta$

$$\left( \frac{\partial}{\partial t} - \nu \tilde{\nabla}^2 \right) \tilde{\nabla}^2 u_r' = \left( \frac{2U_\theta}{r} \right) \frac{\partial^2 u_\theta'}{\partial z^2}$$

$$\left( \frac{\partial}{\partial t} - \nu \tilde{\nabla}^2 \right) u_\theta' + (2A) u_r' = 0$$

Solving PDE

we’re using normal mode solutions (separation of variables) to find the solution

• $$u_{r}^{\prime}=\underbrace{\hat{u}_{r}(r)}_{\text{Unknown radial}}\cdot\underbrace{\cos(nz)}_{\text{From observation}}\cdot\underbrace{e^{st}}_{\text{Growth rate}}$$

• $$u_\theta' = \hat{u}_\theta(r) \cdot \cos(nz) \cdot e^{st}$$

This results in two ODEs

• $$[\nu (\frac{d^2}{dr^2} - n^2) - S] \hat{u}_\theta = 2A \cdot \hat{u}_r$$

• $$[\nu (\frac{d^2}{dr^2} - n^2) - S] [\frac{d^2}{dr^2} - n^2] \hat{u}_r = \left( \frac{2U_\theta n^2}{r} \right) \hat{u}_\theta$$

We can note that our constant for the second equation is not constant

similar rotation approximation

If our rotation is similar though, our base flow circular velocity doesn’t really vary. so can do

• $$\frac{U_{\theta}}{r}=\overline{\Omega}=\frac{1}{2}(\Omega_1+\Omega_2)$$

This gives us:

$$\left[ \nu \left( \frac{d^2}{dr^2} - n^2 \right) - S \right] \left[ \frac{d^2}{dr^2} - n^2 \right] \hat{u}_r = \left( \frac{2 U_\theta n^2}{r} \right) \hat{u}_\theta$$

and if we use equation 3 to get rid of $u_{\theta}$

$$\left[ \nu \left( \frac{d^2}{dr^2} - n^2 \right) - S \right]^2 \left[ \frac{d^2}{dr^2} - n^2 \right] \hat{u}_r = (4\overline{\Omega} n^2 A) \hat{u}_r$$

this is an ODE in just R. Our growth rate S drops out

Boundary conditions:

Now to find boundary conditions we have

• No flow through cylinder walls: $\hat{u}_r = 0 \quad \text{at } r = r_1, r_2$

• No slip on walls $\hat{u}_{z}=0\quad\text{at }r=r_1,r_2$ so also $\frac{\partial u_z'}{\partial z} = 0$

• Continuity: As $\frac{\partial u_r'}{\partial r} + \frac{\partial u_z'}{\partial z} = 0$ ∴ $\frac{d\hat{u}_r}{dr} = 0 \quad \text{at } r = r_1, r_2$

• Also $\hat{u}^{\prime}_{\theta}=0$

Now subbing this into our ODE equation, we can get an equation value at r = r₁ and r₂

$$\left[ \nu \left( \frac{d^2}{dr^2} - n^2 \right) - S \right] \left[ \frac{d^2}{dr^2} - n^2 \right] \hat{u}_r = \left( \frac{2U_\theta n^2}{r} \right) \hat{u}_\theta$$

• no Uθ so right hand size gets removed

we get:

$\frac{d^4 \hat{u}_r}{dr^4} - \left( 2n^2 + \frac{S}{\nu} \right) \frac{d^2 \hat{u}_r}{dr^2} = 0 \quad \text{on } r = (r_1, r_2)$ (this is our 6th boundary condition)

this is because our equation for Uᵣ is 6th order

so we need 6 boundary conditions for uᵣ, 3 on each wall.

Solving ODE

Our system is now:

ODE:

• $$\left[ \nu \left( \frac{d^2}{dr^2} - n^2 \right) - S \right]^2 \left[ \frac{d^2}{dr^2} - n^2 \right] \hat{u}_r = (4\overline{\Omega} n^2 A) \hat{u}_r$$

6 boundary conditions at r = r₁ and r₂

• $\hat{u}_r = 0$:

• $$\frac{d\hat{u}_r}{dr} = 0$$

• $$\frac{d^4\hat{u}_{r}}{dr^4}-\left(2n^2+\frac{S}{\nu}\right)\frac{d^2\hat{u}_{r}}{dr^2}=0$$

Non dimensionalise

Defining two variables to non dimensionlise our lengths

• $X = \frac{r - r_1}{d}$ so X = 0 at inner gap and X = 1 at outer gap

• Axial wavenumber: $a = nd$

This results in the ODE: below with the same boundary conditions:

$\left( \frac{d^2}{dX^2} - a^2 \right)^3 \hat{u}_r = -T \cdot (a^2) \hat{u}_r$ where the taylor number $T = \frac{2(\Omega_1 r_1^2 - \Omega_2 r_2^2)\overline{\Omega}d^3}{\nu^2 r_1} \approx \frac{\text{inertia}}{\text{viscosity}}$

This is extremely hard to solve

Boundary condition simplification

if we simplify the boundary condition $\frac{d^4 \hat{u}_r}{dr^4} - \left( 2n^2 + \frac{S}{\nu} \right) \frac{d^2 \hat{u}_r}{dr^2} = 0 \quad \text{on } r = (r_1, r_2)$

to $\frac{d^4 \hat{u}_r}{dr^4}=0$

Now we just end up with a sinusoidal solution:

$$\hat{u}_r = \sin(N\pi X) \quad \text{where } N = 1, 2, 3 \dots \text{ (Equation 7)}$$

this gives us a critical taylor number of 658

• this is compared to the actual value of 1708 without the simplified boundary conditions

• This is when we set the wave number to 1

Result structure

With our unsimplified structure, we find our critical wavelength by finding the minimum of the taylor number vs wave number curve, for N = 1

• we can note in our critical solution where a = 3.1, we can almost square cells as 3.1 is basically π

Taylor courette flow: Overview of analysis

Where would we find shear flow instabilities

Vorticies in shear flow

If we have a finite discontinuity, this results in an infinitely thin vortex sheet

Examples of shear flow instabilities

Solution for shear flow instabilities

If we apply a small pertubation to a boundary between fast moving and slow moving fluid

$$\eta'(x, t) = \hat{\eta} e^{ikx + st}$$

This gives us the solution where the growth mode is:

$s = -i\frac{1}{2}kU \pm \frac{1}{2}kU$

• We have a pole on the right hand plane, this is unstable and oscillatory

if we write our our wave equation:

$$\eta'(x, t) = \hat{\eta} e^{ik(x - \frac{1}{2}Ut)} e^{\frac{1}{2}kUt}$$

• always unstable

Dispersion relationship for shear flow instabilities

For our solution for the shear flow instability, we can write S in terms of the frequencies.

Our general solution is:

$s = -i\frac{k}{2}(U_1 + U_2) \pm \frac{k}{2}(U_1 - U_2)$ (now subbing in s = -iw) for an a wave solution

$$\omega = \left( \frac{1}{2}k(U_1 + U_2) \right) \pm i \left( \frac{1}{2}k(U_1 - U_2) \right)$$

We can note two things

• The propagation component of our wave is equal to the mean of our two velocities: $\frac12(U_1+U_2)$

• The growth component of our wave is equal to the difference of our two velocities: $\frac12(U_1-U_2)$

• Only ever stable with no shear

further characteristisics of dispersion relation

• Our dispersion relationship isn’t dispersive, our different wavelengths travel at the same speed

• But it is distortive, as high wave numbers will grow much more quickly than small wave numbers

uniform vs stratified shear flow instability

stratified case has mechanism like bouyancy that can supress instability

Shear flow linear analysis: base case

Our base case is shown below, where we have equal pressure at the interface but a discrete velocity jump:

Shear flow linear analysis: governing equations:

• We have an incompressible flow so have continuity expressed like this

• Irrotational flow (shear is contained very locally)

• define flow in terms of velocity potential, allowing for laplacian expression of velocity

Our flow is defined by our laplacian velocity potential above and below our interface

Shear flow linear analysis: boundary conditions

We have three boundary conditions

• Pressure boundary condition: pressure must be the same across the interface

• Kinematic boundary condition: Particles on interface must remain on interface

• Localised disturbance: Initial disturbance occurs in finite region

Pressure boundary condition:

We are defining this using the unsteady bernoulli equation defined in 3A1 databook

$$\frac{p}{\rho} + \frac{1}{2}|\mathbf{u}|^2 + gz + \frac{\partial \phi}{\partial t} = \text{constant throughout flow field}$$

so:

• Below interface: $\frac{\partial \phi_2}{\partial t} = -\left( \frac{p}{\rho} + \frac{\mathbf{u}_2^2}{2} + gz + G_2(t) \right)$

• Above interface: $\frac{\partial \phi_1}{\partial t} = -\left( \frac{p}{\rho} + \frac{\mathbf{u}_1^2}{2} + gz + G_1(t) \right)$

So setting pressures equal (as this is defined at the interface)

$$\frac{p}{\rho} = -\left( \frac{\partial \phi_1}{\partial t} + \frac{\mathbf{u}_1^2}{2} + gz + G_1 \right) = -\left( \frac{\partial \phi_2}{\partial t} + \frac{\mathbf{u}_2^2}{2} + gz + G_2 \right)$$

For steady state

$$\Rightarrow \left( \frac{1}{2}U_1^2 + G_1 \right) = \left( \frac{1}{2}U_2^2 + G_2 \right) \text{ on } z=0$$

Kinematic boundary condition:

Basically we are saying the material derivative across the interface is zero

Our equation of the surface can be defined as:

$$\therefore F(x, z, t) = z - \eta(x, t) = 0 \quad \text{where } F=0 \text{ on } z = \eta(x, t)$$

Now the material derivative is:

$$\frac{DF}{Dt} = \left[ \frac{\partial F}{\partial t} + (\mathbf{u} \cdot \nabla)F \right] = 0$$

Expanding this out we get:

$\frac{\partial \eta}{\partial t} + u \frac{\partial \eta}{\partial x} = w$ (basically change in boundary height + horizontal transport term must equal vertical velocity)

Localised disturbance

We assume our initial disturbance occurs within a finite extent, so our stream function reverts back to our original value when we move far away.

• $$\nabla\phi_2\rightarrow U_2\mathbf{i}\quad\text{as}\quad z\rightarrow-\infty\text{ }$$

• $$\nabla\phi_1\rightarrow U_1\mathbf{i}\quad\text{as}\quad z\rightarrow+\infty\text{ }$$

Shear flow linear analysis: inputting perturbations

Adding perturbations

We are now adding in our perturbations:

• Perturb velocity: $\mathbf{u} = \mathbf{U} + \mathbf{u}'(x, z, t)$

• Perturb pressure: $p = P + p'(x, z, t)$

• Perturb 'interface': $\eta = 0 + \eta'(x, t)$

Subbing into governing equations

subbing this into our laplacian governing equations, we get:

\nabla^2 \phi_1' = 0 \quad \text{for } z > 0

\nabla^2 \phi_2' = 0 \quad \text{for } z < 0

Changing z coordinate

We also replace our boundary position $z = \eta'$ to $z = 0$

• this is because η is very small, we can see this with a taylor expansion of our velocity potential

• $$\phi'(x, \eta', t) \approx \phi'(x, 0, t) + \underbrace{\eta' \cdot \frac{\partial \phi'}{\partial z} \bigg|_{z=0}}_{\text{The "Second-Order" Term}} + \dots$$

• almost no difference

Subbing into boundary conditions

Our pressure boundary condition is:

• $$\left(\frac{\partial\phi_1}{\partial t}+\frac{\mathbf{u}_1^2}{2}+gz+G_1\right)=\left(\frac{\partial\phi_2}{\partial t}+\frac{\mathbf{u}_2^2}{2}+gz+G_2\right)$$

• subbing in our perturbation terms, this results in binomial expansion about u₁²

• $$\left(G_1+\frac{U_1^2}{2}+U_1\frac{\partial\phi_1'}{\partial x}+\frac{(\nabla\phi_1')^2}{2}+g\eta^{\prime}+\frac{\partial\phi_1'}{\partial t}\right)=\left(G_2+\frac{U_2^2}{2}+U_2\frac{\partial\phi_2'}{\partial x}+\frac{(\nabla\phi_2')^2}{2}+g\eta^{\prime}+\frac{\partial\phi_2'}{\partial t}\right)\quad\text{on }z=\eta^{\prime}\text{ }$$

  • can simplify using base state since z=η’ and z=0 is basically the same, so $\Rightarrow \left( \frac{1}{2}U_1^2 + G_1 \right) = \left( \frac{1}{2}U_2^2 + G_2 \right) \text{ on } z=0$

  • our grad squared term is small and can be removed

This results in:

$$\left(U_1\frac{\partial\phi_1'}{\partial x}+\frac{\partial\phi_1'}{\partial t}\right)=\left(U_2\frac{\partial\phi_2'}{\partial x}+\frac{\partial\phi_2'}{\partial t}\right)\quad\text{on }z=0\text{ }$$

Kinematic boundary condition is:

$$\frac{\partial \eta}{\partial t} + u \frac{\partial \eta}{\partial x} = w$$

• removing small terms, we get:

• Upper layer: $\frac{\partial\eta'}{\partial t}+U_1\frac{\partial\eta'}{\partial x}=\frac{\partial\phi_1'}{\partial z}\quad\text{on }z=0\text{ }$

• Lower layer: $\frac{\partial\eta'}{\partial t}+U_2\frac{\partial\eta'}{\partial x}=\frac{\partial\phi_2'}{\partial z}\quad\text{on }z=0$

Localised disturbance boundary condition

This with our perturbations is just:

• $$\nabla\phi^{\prime}_2\rightarrow0\text{ as}\quad z\rightarrow-\infty\text{ }$$

• $$\nabla\phi^{\prime}_1\rightarrow0\quad\text{as}\quad z\rightarrow+\infty\text{ }$$

Shear flow linear analysis: finding solution

Finding ODE

We want normal solutions to our linearised small disturbance laplace equations. this is simply:

• \nabla^2 \phi_1' = 0 \quad \text{for} \quad z > 0

• \nabla^2 \phi_2' = 0 \quad \text{for} \quad z < 0

Our normal mode solutions are of the form:

• $$\phi_1^{\prime}=\underbrace{\hat{\phi}_1(z)}_{\substack{\text{vertical structure}}}\cdot e^{i(\overbrace{k}^{\substack{\\\text{wavenumber }k}})x+(\overbrace{s}^{\substack{\text{growth} \\ \text{rate } s}})t}$$

• $$\phi_2' = \hat{\phi}_2(z) \cdot e^{ikx+st}$$

• $$\eta' = \hat{\eta} \cdot e^{ikx+st}$$

Plugging this into our laplacian we get:

• Upper layer \frac{d^2}{dz^2}\hat{\phi}_1(z) - k^2\hat{\phi}_1(z) = 0 \implies \hat{\phi}_1(z) = Ae^{kz} + Be^{-kz}, \quad z > 0

  • We need this to be bounded so: \hat{\phi}_1(z) = Be^{-kz} \quad z > 0

• Lower layer \frac{d^2}{dz^2}\hat{\phi}_2(z)-k^2\hat{\phi}_2(z)=0\implies\hat{\phi}_2(z)=Ce^{kz}+De^{-kz},\quad z<0\text{ }

  • Again bounded so: \hat{\phi}_2(z) = Ce^{kz} \quad z < 0.

Applying boundary equations

Our solutions are in the form:

$$\phi_1^{\prime}=\underbrace{Be^{-kz}}_{=\hat{\phi}_1}\cdot e^{ikx+st}\quad\phi_2^{\prime}=\underbrace{Ce^{kz}}_{=\hat{\phi}_2}\cdot e^{ikx+st}\text{ }$$

Now we can apply our boundary conditions:

• Pressure:$\left(U_1\frac{\partial\phi_1'}{\partial x}+\frac{\partial\phi_1'}{\partial t}\right)=\left(U_2\frac{\partial\phi_2'}{\partial x}+\frac{\partial\phi_2'}{\partial t}\right)\quad\text{on }z=0\text{ }$

• Kinematic:

  • Upper layer: $\frac{\partial\eta'}{\partial t}+U_1\frac{\partial\eta'}{\partial x}=\frac{\partial\phi_1'}{\partial z}\quad\text{on }z=0\text{ }$

  • Lower layer: $\frac{\partial\eta'}{\partial t}+U_2\frac{\partial\eta'}{\partial x}=\frac{\partial\phi_2'}{\partial z}\quad\text{on }z=0$

The kinematic boundary conditions give:

$B = -(s + U_1 ik) \frac{\hat{\eta}}{k}$ and $C = (s + U_2 ik) \frac{\hat{\eta}}{k}$

The pressure boundary conditions give:

$$(s + U_1 ik) B = (s + U_2 ik) C$$

This gives me a quadratic that we can solve for S

$$2s^2 + s 2ik(U_1 + U_2) - k^2(U_2^2 + U_1^2) = 0$$

Solution:

This gives us the solution:
$s = -i\frac{k}{2}(U_1 + U_2) \pm \frac{k}{2}(U_1 - U_2)$

we cannot form a dimensionless group because there is no stabilising force against the shear

Shear flow linear analysis: Physical mechanism and rayleight’s inflexion point theorem

We can see the physical mechanism for the instability

The vorticity of the sheet results in our sine wave perturbation being amplified

Inflection point theorem

Instead if we assume a stream function: $\psi = f(z)e^{st + ikx}$

The momentum equation results in this result: $\int s_i U''(z) \frac{|f|^2}{|s + ikU|^2} dz = 0$

• If we have an unstable mode, so sᵢ > 0

• For this to work U’’(z) needs to change sign

• so ONE of the conditions (but not only) for instability is the existence of an inflection point

Consequence of this

• Explains boundary layer separation whenever we form an inversion

• Why jets are unstable

Summary of kelvin helmholtz instability

This is the result of shear instability being stabilised with buoyancy. The growth rate is given by this equation:

$$s = -ik \frac{\rho_1 U_1 + \rho_2 U_2}{\rho_1 + \rho_2} \pm \left[ \frac{k^2 \rho_1 \rho_2 (U_1 - U_2)^2}{(\rho_1 + \rho_2)^2} - gk \frac{(\rho_2 - \rho_1)}{\rho_1 + \rho_2} \right]^{1/2}$$

• The buoyancy acts to stabilise the growth
  

Kelvin helmholtz linear analysis: Base solution

the base solution can be seen in the image

now we have varying pressures

Kelvin helmholtz linear analysis: Boundary conditions

These are identical to our shear layer except for the pressure boundary condition which varies with density now.

So:

Pressure: $\rho_1\left(\frac{\partial\phi_1}{\partial t}+\frac{\mathbf{u}_1^2}{2}+gz+G_1\right)=\rho_2\left(\frac{\partial\phi_2}{\partial t}+\frac{\mathbf{u}_2^2}{2}+gz+G_2\right)$ densities are now shown

• Our steady state is now:$\rho_1\left(\frac{\mathbf{u}_1^2}{2}+G_1\right)=\rho_2\left(\frac{\mathbf{u}_2^2}{2}+G_2\right)$

Kinematic:-

• Upper layer: $\frac{\partial\eta'}{\partial t}+U_1\frac{\partial\eta'}{\partial x}=\frac{\partial\phi_1'}{\partial z}\quad\text{on }z=0\text{ }$

• Lower layer: $\frac{\partial\eta'}{\partial t}+U_2\frac{\partial\eta'}{\partial x}=\frac{\partial\phi_2'}{\partial z}\quad\text{on }z=0$

Localisation:

• $$\nabla\phi_2\rightarrow U_2\mathbf{i}\quad\text{as}\quad z\rightarrow-\infty\text{ }$$

• $$\nabla\phi_1\rightarrow U_1\mathbf{i}\quad\text{as}\quad z\rightarrow+\infty\text{ }$$

Kelvin helmholtz linear analysis: Introducing perturbations

Still adding the same pertubations

• Perturb velocity: $\mathbf{u} = \mathbf{U} + \mathbf{u}'(x, z, t)$

• Perturb pressure: $p = P + p'(x, z, t)$

• Perturb 'interface': $\eta = 0 + \eta'(x, t)$

subbing this into our laplacian governing equations, we get:

\nabla^2 \phi_1' = 0 \quad \text{for } z > 0

\nabla^2 \phi_2' = 0 \quad \text{for } z < 0

Same governing equations as before

Changing z coordinate

Same analysis in shear flow where we ignore the small distance caused by the perturbations and take our values at z = 0

Linearised boundary equations

Our kinematic and localisation boundary conditions are exactly the same, when we’re linearising.

Kinematic

• Upper layer: $\frac{\partial\eta'}{\partial t}+U_1\frac{\partial\eta'}{\partial x}=\frac{\partial\phi_1'}{\partial z}\quad\text{on }z=0\text{ }$

• Lower layer: $\frac{\partial\eta'}{\partial t}+U_2\frac{\partial\eta'}{\partial x}=\frac{\partial\phi_2'}{\partial z}\quad\text{on }z=0$

localisation:

• $$\nabla\phi^{\prime}_2\rightarrow0\text{ as}\quad z\rightarrow-\infty\text{ }$$

• $$\nabla\phi^{\prime}_1\rightarrow0\quad\text{as}\quad z\rightarrow+\infty\text{ }$$

pressure

we can no longer ignore the hydrostatic term as before:

so original equation is:

• $$\rho_1\left(G_1+\frac{U_1^2}{2}+U_1\frac{\partial\phi_1'}{\partial x}+\frac{(\nabla\phi_1')^2}{2}+g\eta^{\prime}+\frac{\partial\phi_1'}{\partial t}\right)=\rho_2\left(G_2+\frac{U_2^2}{2}+U_2\frac{\partial\phi_2'}{\partial x}+\frac{(\nabla\phi_2')^2}{2}+g\eta^{\prime}+\frac{\partial\phi_2'}{\partial t}\right)\quad\text{on }z=\eta^{\prime}\text{ }$$

• can still simplify with base state since z=η’ and z=0 is basically the same, so $\Rightarrow \left( \frac{1}{2}U_1^2 + G_1 \right) = \left( \frac{1}{2}U_2^2 + G_2 \right) \text{ on } z=0$

• our grad squared term is small and can be removed, but no longer can remove hydrostatic term ρgη’

Find normal mode solutions

Plugging in our normal mode solutions into laplace’s equation:$\phi_2' = \hat{\phi}_2(z) \cdot e^{ikx+st}$ etc..

We have the exact same form (since same governing equations)

so:

$$\phi_1^{\prime}=\underbrace{Be^{-kz}}_{=\hat{\phi}_1}\cdot e^{ikx+st}\quad\phi_2^{\prime}=\underbrace{Ce^{kz}}_{=\hat{\phi}_2}\cdot e^{ikx+st}\text{ }$$

exactly the same using kinematic to find B and C

$B = -(s + U_1 ik) \frac{\hat{\eta}}{k}$ and $C = (s + U_2 ik) \frac{\hat{\eta}}{k}$

different pressure boundary condition for different relationship between B and C

$$\rho_1 \left[ (s + ikU_1)B + g\hat{\eta} \right] = \rho_2 \left[ (s + ikU_2)C + g\hat{\eta} \right]$$

• Compare with $(s + U_1 ik) B = (s + U_2 ik) C$ for same density

We get our quadratic now:

$$\rho_1 [kg - (s + U_1 ik)^2] = \rho_2 [kg + (s + U_2 ik)^2]$$

Result:

this result in a growth rate: (once solved with quadratic formula)

$$s = - ik \frac{\rho_1 U_1 + \rho_2 U_2}{\rho_1 + \rho_2} \pm \left[ \frac{k^2 \rho_1 \rho_2 (U_1 - U_2)^2}{(\rho_1 + \rho_2)^2} - gk \frac{(\rho_2 - \rho_1)}{\rho_1 + \rho_2} \right]^{1/2}$$

Physical implication of kelvin helmholtz

Our poles for kelvin helmholtz lie on:

$$s = - ik \frac{\rho_1 U_1 + \rho_2 U_2}{\rho_1 + \rho_2} \pm \left[ \frac{k^2 \rho_1 \rho_2 (U_1 - U_2)^2}{(\rho_1 + \rho_2)^2} - gk \frac{(\rho_2 - \rho_1)}{\rho_1 + \rho_2} \right]^{1/2}$$

• This is unstable if the second term is real

• Thus our stable criteria is if the second term is imaginary: so

\left[\frac{k^2\rho_1\rho_2(U_1-U_2)^2}{(\rho_1+\rho_2)^2}-gk\frac{(\rho_2-\rho_1)}{\rho_1+\rho_2}\right]<0

This results in a stability criteria of:

\Rightarrow \text{instability for wavenumbers } k > g \frac{(\rho_2^2 - \rho_1^2) / \rho_1 \rho_2}{(U_1 - U_2)^2}

• Low wavelengths are stabilised by bouyancy

• This is because it involves moving large masses of fluids, while at smaller scale the shear forces dominate.

Differences between uniform shear and kelvin helmholtz

Energy analysis of general stratification

We can form a criterion for stratification by considering work.

• We are considering two particles swapping places

• Change in work from bouyancy

• Change in kinetic energy

This is unstable if ΔKE > work done against bouyancy

Looking at a continuous strata where our density and velocity varies continuously with height.

Work done against bouyancy:

$$F_{buoy} = \{ \underbrace{\rho(z)}_{\text{density of particle}} - \underbrace{[\rho(z) + a \frac{d\rho}{dz} + \dots]}_{\rho(z+a) \text{ density of environment}} \} gV \approx -\left( a \frac{d\rho}{dz} \right) gV$$

$$\text{work done} = \int_{a=0}^{a=\delta z} F_{buoy} \cdot d\mathbf{r} = \int_{0}^{\delta z} -\left( a \frac{d\rho}{dz} gV \right) da = -\frac{1}{2}(\delta z)^2 \frac{d\rho}{dz} gV$$

change in kinetic energy:

• Kinetic energy before switching: (density changes don’t affect KE, only hydrostatic)

  • $$KE_{before} = \frac{1}{2}m_1 u^2 + \frac{1}{2}m_2 (u + \delta u)^2 \quad \text{where } m_1 \approx m_2 = \rho_0 V$$

• Kinetic energy after swapping: (assuming they take mean velocity)

  • $$KE_{after} = \frac{1}{2}m_1 \left( \frac{u + (u + \delta u)}{2} \right)^2 + \frac{1}{2}m_2 \left( \frac{u + (u + \delta u)}{2} \right)^2$$

change is thus:

$$\Delta KE = KE_{before} - KE_{after} = \frac{\rho_0 V}{4}(\delta u)^2$$

Unstable if:

\frac{\rho_0 V}{4}(\delta u)^2 > -2 \times \frac{1}{2} (\delta z)^2 \frac{d\rho}{dz} gV (two added for two particles)

If we’re taking this as a derivative by δu/δz ≈ du/dz

\frac{1}{4} > \frac{-\frac{g}{\rho_0} \frac{d\rho}{dz}}{\left( \frac{du}{dz} \right)^2} = \frac{N^2}{\left( \frac{du}{dz} \right)^2}

This dimensionless ratio is the Richardson Number ($Ri$).

$$Ri = \frac{N^2}{(du/dz)^2}$$

Condition for instability: $Ri < \frac{1}{4}$

Rayleigh bernard convection overview

Here we have a similar flow structure to taylor courette, modelling the flow between two plates at different temperatures:

• The steady state solution is just diffusive flow, where we have heat diffusion from one side to another

• We have a convective mode solution, where bouyancy results in instability

Analogy to taylor courette

There is an exact analogy to thin gap taylor courette

• In taylor courette we have the taylor number $T = \frac{2(\Omega_1 r_1^2 - \Omega_2 r_2^2)\overline{\Omega}d^3}{\nu^2 r_1} \approx \frac{\text{inertia}}{\text{viscosity}}$

• In rayleigh bernard we have the rayleigh number: Ra = \frac{g \beta \Delta T d^3}{\kappa \nu} > Ra_{crit}

• Both of these result in stability when we exceed 1708

There is a perfect physical analogy here too

• In taylor courette, if we enter a rotating frame. We have an instability driven by a centrifugal force, with a fixed boundary at the two walls.

  • This centrifugal force is constant with r just like a constant g for a thin gap

• In rayleigh bernard. We have an instability driven by bouyancy (gravity field) (temperature), with a fixed boundary at two walls

• Base case both have a linear potential field (centrifugal driven, and gravity) (when small gap)

• This is destabilising against viscosity

This is why we result in the same maths

Rayleigh bernard linear analysis: Governing equations and boundary equations

Governing equations:

Our governing equations are:

• Navier stokes: $\frac{\partial \mathbf{u}}{\partial t} + (\mathbf{u} \cdot \nabla) \mathbf{u} = -\frac{1}{\rho} \nabla p + \nu \nabla^2 \mathbf{u} + \mathbf{g}$ (assume viscosity constant with temperature)

• continuity: $\frac{D\rho}{Dt} + \rho \nabla \cdot \mathbf{u} = 0 \quad \xrightarrow{\text{Density variations small}} \quad \frac{D\rho}{Dt} = 0 \implies \nabla \cdot \mathbf{u} = 0$ (assume constant density for continuity)

• Thermal diffusion: $\frac{\partial T}{\partial t} + (\mathbf{u} \cdot \nabla) T = \kappa \nabla^2 T$

• Equation of state: $\rho = \tilde{\rho} [1 - \beta(T - \tilde{T})]$

Compared to taylor-courette

• Same navier stokes and continuity to govern the fluid flow

• Our thermal diffusion is lumped within navier stokes, as it is momentum diffusion that results the change in fluid energy

• No equation of state, also lumped into navier stokes as a centrifugal force term.

Boundary conditions:

Temperature BC

• $T = T_{cool}$ on $z = d$ for all $t$.

• $T = T_{warm}$ on $z = 0$ for all $t$.

This is analogous to our no slip condition in taylor courette, as both of these boundary conditions define a fluid energy within our system

Velocity BC

1. No flow through boundaries: Vertical velocity $w = 0$ at the walls.

2. No slip: Horizontal velocities $u = v = 0$ at the walls.

Same as taylor-courette

We can see how for the boundary conditions too, taylor courette sort of lumps terms together

Rayleigh bernard linear analysis: Base flow solution

This base flow involves no motion, so our base flow is purely a heat diffusion equation.

Solving steady temperature diffusion equation:

$\frac{\partial T}{\partial t} + (\mathbf{U} \cdot \nabla) T = \kappa \nabla^2 T$

• No transport or unsteady term, so just a laplacian solution:

$$0 + 0 = \kappa \frac{d^2}{dz^2} T_0(z)$$

Plugging in our boundary conditions this has the following solution:

$$T_0(z) = T_{warm} - \frac{z}{d}(T_{warm} - T_{cool})$$

With analogous pressure and density solutions:

Density: $\rho_0(z) = \tilde{\rho} [1 - \beta(T_0(z) - \tilde{T})]$

Pressure: $\frac{dP_0}{dz} = -\rho_0(z)g$

Analogy to taylor couette

This is exactly analogous to the thin gap taylor couette solution

• In a thin gap we get a linear increase in azimuthal velocity

• In a rotating frame the flow is also steady

• But like here, we have an decrease in flow energy from one plate to the other, but driven by momentum diffusion rather than thermal diffusion

  • Our centrifugal potential energy decreases

  • Compared to our gravitational potential energy

Rayleigh bernard linear analysis: Introducing perturbations and finding PDE

Perturbations:

We are introducing the following perturbations:

• Velocity: $\mathbf{u} = \mathbf{0} + \mathbf{u}' \quad \text{[cite: 173]}$

• Temperature $T = T_0(z) + T' \quad \text{[cite: 173]}$

• Density: $\rho = \rho_0(z) + \rho' \quad \text{[cite: 174]}$

• pressure $p = P_0(z) + p' \quad \text{[cite: 174]}$

Plugging into governing equations

Remember our governing equations are:

• Navier stokes: $\frac{\partial \mathbf{u}}{\partial t} + (\mathbf{u} \cdot \nabla) \mathbf{u} = -\frac{1}{\rho} \nabla p + \nu \nabla^2 \mathbf{u} + \mathbf{g}$ (assume viscosity constant with temperature)

• continuity: $\frac{D\rho}{Dt} + \rho \nabla \cdot \mathbf{u} = 0 \quad \xrightarrow{\text{Density variations small}} \quad \frac{D\rho}{Dt} = 0 \implies \nabla \cdot \mathbf{u} = 0$ (assume constant density for continuity)

• Thermal diffusion: $\frac{\partial T}{\partial t} + (\mathbf{u} \cdot \nabla) T = \kappa \nabla^2 T$

• Equation of state: $\rho = \tilde{\rho} [1 - \beta(T - \tilde{T})]$

Plugging into navier stokes

When plugging into navier stokes we get:

$$(\rho_0+\rho^{\prime})\left[\frac{\partial\mathbf{u}'}{\partial t}+(\mathbf{u}^{\prime}\cdot\nabla)\mathbf{u}^{\prime}\right]=-\nabla(P_0+p^{\prime})+(\rho_0+\rho^{\prime})\nu\nabla^2\mathbf{u}^{\prime}+(\rho_0+\rho^{\prime})\mathbf{g}\quad$$

Now ignoring products of small temrs and introducing a hydrostatic balance $\frac{dP_0}{dz} = -\rho_0 g$

$$\rho_0 \frac{\partial \mathbf{u}'}{\partial t} = -\nabla p' + \rho_0 \nu \nabla^2 \mathbf{u}' + \rho'\mathbf{g} \quad$$

Because density changes are small, we are rewriting it in terms of a constant reference density

$$\tilde{\rho} \frac{\partial \mathbf{u}'}{\partial t} = -\nabla p' + \tilde{\rho} \nu \nabla^2 \mathbf{u}' + \rho'\mathbf{g} \quad \text{[cite: 194]}$$

Other governing equations

• Continuity: we get $\nabla \cdot \mathbf{u}' = 0$

• Thermal diffusion: we get: $\frac{\partial T'}{\partial t}+w^{\prime}\left(\frac{dT_0}{dz}\right)=\kappa\nabla^2T^{\prime}$

• Equation of state: $\rho' = -\beta \tilde{\rho} T'$

Combining linearised governing equations:

Our goal is an expression for vertical velocity This is because convection gives vertical velocity

Eliminate pressure in NS

We can take the curl to eliminate pressure (as it is a potential field)

$$\left( \frac{\partial}{\partial t} - \nu \nabla^2 \right) \nabla \times \mathbf{u}' = -\beta (\nabla T') \times \mathbf{g}$$

If we take the curl again, we can simplify further with the identities

• $\nabla \times \nabla \times \mathbf{u}' = \nabla(\nabla \cdot \mathbf{u}') - \nabla^2 \mathbf{u}'$ and $\nabla \cdot \mathbf{u}' = 0$, leaving just a laplacian term

$\Rightarrow -\left( \frac{\partial}{\partial t} - \nu \nabla^2 \right) \nabla^2 \mathbf{u}' = -\nabla \times (\beta (\nabla T') \times \mathbf{g})$ which simplifies with more vector identities to:

• $\nabla\times(\mathbf{F}\times\mathbf{G})=(\mathbf{G}\cdot\nabla)\mathbf{F}-(\mathbf{F}\cdot\nabla)\mathbf{G}+\mathbf{F}(\nabla\cdot\mathbf{G})-\mathbf{G}(\nabla\cdot\mathbf{F})$ so:

• $$\nabla \times [\beta (\nabla T') \times \mathbf{g}] = (\mathbf{g} \cdot \nabla) \beta (\nabla T') - \underbrace{(\beta (\nabla T') \cdot \nabla) \mathbf{g}}_{=0 \text{ as } \nabla \mathbf{g} = 0} + \underbrace{\beta (\nabla T') (\nabla \cdot \mathbf{g})}_{=0 \text{ as } \nabla \cdot \mathbf{g} = 0} - \mathbf{g} (\nabla \cdot \beta (\nabla T'))$$

$$\left( \frac{\partial}{\partial t} - \nu \nabla^2 \right) \nabla^2 \mathbf{u}' = \beta \left[ (\mathbf{g} \cdot \nabla) \nabla T' - \mathbf{g} \nabla^2 T' \right]$$

looking at vertical component

Using the fact gravity only acts in vertical direction g = (0,0,-g)

$$\left( \frac{\partial}{\partial t} - \nu \nabla^2 \right) \nabla^2 w' = \beta g \left( \frac{\partial^2}{\partial x^2} + \frac{\partial^2}{\partial y^2} \right) T'$$

eliminate temperature perturbation (so we just have our equation in w’)

We will use the energy equation to do this, rewriting our energy equation we get:

• $\frac{\partial T'}{\partial t}+w^{\prime}\left(\frac{dT_0}{dz}\right)=\kappa\nabla^2T^{\prime}$ rearranging we get:

• $$\left( \frac{\partial}{\partial t} - \kappa \nabla^2 \right) T' = -w' \frac{dT_0}{dz} \quad [cite: 267]$$

PDE solution

Subbing this into the vertical component equation we get:

$$\left( \frac{\partial}{\partial t} - \kappa \nabla^2 \right) \left( \frac{\partial}{\partial t} - \nu \nabla^2 \right) \nabla^2 w' = -\beta g \left( \frac{\partial^2}{\partial x^2} + \frac{\partial^2}{\partial y^2} \right) w' \frac{dT_0}{dz} \quad [cite: 270]$$

Rayleigh bernard linear analysis: Solving PDE

Forming ODE and comparing to taylor-courette

Our PDE is now: $\left(\frac{\partial}{\partial t}-\kappa\nabla^2\right)\left(\frac{\partial}{\partial t}-\nu\nabla^2\right)\nabla^2w^{\prime}=-\beta g^{\prime}\frac{dT_0}{dz}\left(\frac{\partial^2}{\partial x^2}+\frac{\partial^2}{\partial y^2}\right)w^{\prime}$

• compare with taylor-courette:

  • $$\left( \frac{\partial}{\partial t} - \nu \tilde{\nabla}^2 \right) \tilde{\nabla}^2 u_r' = \left( \frac{2U_\theta}{r} \right) \frac{\partial^2 u_\theta'}{\partial z^2}$$

  • $\left( \frac{\partial}{\partial t} - \nu \tilde{\nabla}^2 \right) u_\theta' + (2A) u_r' = 0$ if we sub in for uθ in the derivative

  • $$\left(\frac{\partial}{\partial t}-\nu\nabla^2\right)\left(\frac{\partial}{\partial t}-\nu\nabla^2\right)\nabla^2u_{r}^{\prime}=-(2A)\frac{2U_\theta}{r}\left[\frac{\partial^2}{\partial z^2}\right]u_{r}^{\prime}$$

• Exactly in the same form: $\underbrace{\left(\frac{\partial}{\partial t}-\text{Diff}_1\right)\left(\frac{\partial}{\partial t}-\text{Diff}_2\right)}_{\text{The Two Brakes}}\underbrace{\nabla^2\Psi^{\prime}}_{\text{Stabilizing Curvature}}=\underbrace{-\mathcal{C}}_{\text{Coupling}}\underbrace{\left(\frac{d\Phi_0}{d\eta}\right)}_{\text{Energy Gradient}}\underbrace{\nabla_{\perp}^2\Psi^{\prime}}_{\text{Symmetry-Breaking Laplacian}}\text{ }$

• we have different geometry though, since taylor-courette our cells must be stacked vertically, while in rayleigh bernard they can spread out in a plane

Applying modal solutions

We assume a solution of the form:

$$w' = \hat{w}(z) f(x,y) e^{st}$$

• $\hat{w}(z)$: The unknown vertical structure of the roll.

• $f(x,y)$: The horizontal pattern; we assume there is no preferred horizontal direction.

• $e^{st}$: Introduces the growth rate ($s$) to determine if the perturbation grows (s > 0) or decays.

This gives us the ODE

$$\left[ \nu \left( \frac{d^2}{dz^2} - a^2 \right) - s \right] \left[ \kappa \left( \frac{d^2}{dz^2} - a^2 \right) - s \right] \left( \frac{d^2}{dz^2} - a^2 \right) \hat{w} = \left( \beta g \frac{dT_0}{dz} a^2 \right) \hat{w}$$

Workings for the Substitution:

• Temporal derivative: $\frac{\partial}{\partial t} w' = s f(x,y) e^{st} \hat{w}$

• Spatial Laplacian: $\nabla^2 w' = \left[ \left( \frac{\partial^2 f}{\partial x^2} + \frac{\partial^2 f}{\partial y^2} \right) e^{st} + f e^{st} \frac{d^2}{dz^2} \right] \hat{w} = \left( \frac{d^2}{dz^2} - a^2 \right) f e^{st} \hat{w}$

Applying boundary conditions

We need 6 boundary conditions for w’

• No penetration boundary condition: $w^{\prime}=0\implies\hat{w}=0\space at \space z=(0,d)$

• no slip boundary condition $\frac{d\hat{w}}{dz} = 0 \quad \text{on } z = (0, d)$ (from no slip and continuity)

• energy anchor boundary condition $T' = 0$ on $z = (0, d)$

  • Applying this simplified navier stokes $\left( \frac{\partial}{\partial t} - \nu \nabla^2 \right) \nabla^2 w' = \beta g \left( \frac{\partial^2}{\partial x^2} + \frac{\partial^2}{\partial y^2} \right) T'$ so we get this in terms of w’

  • We get: $\left.\left(\frac{\partial}{\partial t}-\nu\nabla^2\right)\nabla^2w^{\prime}\right|_{z=0,d}=0$ which when differentiated and with our modal solution gets us

  • $$\left.\frac{d^4 \hat{w}}{dz^4}-\left(2a^2+s/\nu\right)\frac{d^2 \hat{w}}{dz^2}\right|_{z=0,d}=0$$

Again analogous to taylor-courette, we have the same two no penetration and no slip boundary conditions

• This condition: $\frac{d^4 \hat{u}_r}{dr^4} - \left( 2k^2 + \frac{s}{\nu} \right) \frac{d^2 \hat{u}_r}{dr^2} = 0 \quad \text{on } r = (r_1, r_2) \text{}$ is equivalent to a $u_\theta' = 0$ condition in uᵣ This is an energy anchor condition

Solution:

we get an identical solution to taylor courette as shown in picture.

if we solve with the simplified boundary condition set so:

• No-Slip: $\hat{w} = 0$ $\quad \longrightarrow \quad$ Stress-Free: $\hat{w} = 0$ (no change)

• No-Slip: $\frac{d\hat{w}}{dz} = 0$ $\quad \longrightarrow \quad$ Stress-Free: $\frac{d^2\hat{w}}{dz^2} = 0$ (from no slip to no shear stress)

• No-Slip: $\frac{d^4 \hat{w}}{dz^4} - (2a^2 + s/\nu)\frac{d^2 \hat{w}}{dz^2} = 0$ $\quad \longrightarrow \quad$ Stress-Free: $\frac{d^4 \hat{w}}{dz^4} = 0$ (follows from previous assumptions)
  

This enables our sinusoidal solution:

$$\hat{w}=\sin\left(N\pi\frac{z}{d}\right)$$

stability criteria

With our simplified boundary conditions, we’re looking when the wave number is 1

We find a definition for the rayleigh number in terms of our wavenumber by plugging the sin solution into the ODE

$$Ra = \frac{d^4}{a^2} \left[ \frac{N^2\pi^2}{d^2} + a^2 \right]^3$$

• Critical rayleigh number is the minimum value of this when N = 1, equals 658

Our real value is 1708 like with taylor courette

Flow structure of rayleigh-bernard

Bernard Marangoni

This looks similar to rayleigh bernard, but the driving force is NOT a density gradient due to temperature but rather a surface tension gradient due to temperature.

• This is defined by a different Maragoni number

• This is unstable if$Ma = -\frac{d\gamma}{dT} \frac{(T_0 - T_1) d}{\rho \nu \kappa}$ > 80

• Often both effects are combined, described by if \frac{Ra}{Ra_{crit}} + \frac{Ma}{Ma_{crit}} > 1

• Rayleigh dominates for thicker layers

Role of symmetry breaking in instability

Instability can often been seen in the perspective of breaking symmetries:

The 4 Core Examples of Symmetry Breaking

Type of Symmetry	The "Perfect" (Symmetric) State	The Broken (Instability) State	Example
Reflection	A flat plate aligned perfectly with a steady flow.	The plate spontaneously bends to one side.	A cantilevered postcard in the wind.
Time-Translation	A steady input (like blowing air at a constant rate).	The output becomes oscillatory (vibration).	Blowing a raspberry or vocal cords vibrating.
Spatial Translation	An infinitely long, uniform thread of liquid.	The thread pinches off into individual droplets.	Rayleigh-Plateau instability (dripping tap).
Axisymmetry	Wine coating a glass in a perfect, smooth circle.	The wine falls in distinct, repeating "tears" or legs.	Marangoni effect (tears of wine).

Role of positive feedback in breaking symmetry

Symmetries are often broken through positive feedback loops.

Example of aeroelastic divergence:

• If we consider a wing with a centre of lift infront of the torsional center

• Doing a torque balance: $T = \underbrace{\frac{1}{2} C_L \rho U^2 c w d \sin \theta \cos \theta}_{\text{Fluid Torque (Destabilizing)}} - \underbrace{\kappa \theta}_{\text{Spring Torque (Restoring)}}$

• We can see we get an increasing torque, if:

• \frac{C_L \rho U^2 c w d}{2\kappa} > \frac{\theta}{\sin \theta \cos \theta}

This was a big issue with early flight.

Formalising analysis of lienar stability theory

Comparing with the picture there are 8 steps to linear analysis.

1. State and state variables

We first need to find the quantities we are using the define the system, these are our state variables.

• For the aeroelastic plate example, the system state is only defined by θ(t)

2. Governing dynamics

We then need to work out the equations that link our state variables together and define our system.

• $\frac{d\mathbf{x}}{dt} = \mathcal{F}(\mathbf{x}; p)$ (our equation for how our variable changes over time)

• This depends on our state variables (x) like position, and system parameters (p) like density, freestream velocity etc.

Enslaved variables:

• Many variables respond so quickly that they are “enslaved” to others

• For instance our lift force is enslaved to θ, since we assume the flow field responds fast enough to reach a steady state

• As such we are describing the system with just θ

symmetries and approximations

• Many systems are invariant in our state variables, and have symmetries

• Ie our plate situation is symmetric about θ, and also time invariant (where a shift in time doesn’t change F)

3. Steady state

We first want to find a steady state solution, this should obey the symmetry which is being broken by the instability ie:

• Time-Translation Symmetry: We look for a state where $$\theta$</span> or velocity is constant. It looks the same whether you look now or ten minutes from now.</p></li><li><p><strong>Spatial Translation Symmetry:</strong> In a liquid thread (Rayleigh-Plateau), the steady state is a perfect cylinder. Every point along the axis is identical.</p></li><li><p><strong>Reflection Symmetry:</strong> For your hinged plate, the steady state is <span>$\theta = 0$</span>. This is the only state that is identical to its mirror image. </p></li></ul><p></p><p><strong>Floquet analysis</strong></p><p>Some systems don’t have a steady state in time, such a periodic systems like a swing</p><ul><li><p>We then consider if there is growth after a period.</p></li></ul><p></p><h4 id="8f04a236-725d-4f2c-a268-6362b7a323e9" data-toc-id="8f04a236-725d-4f2c-a268-6362b7a323e9" collapsed="false" seolevelmigrated="true">4. Introducing a perturbation</h4><p>We are now introducing a perturbation onto our steady state. Ie nudging our state,.</p><p>$\mathbf{x} = \mathbf{x}_0 + \mathbf{x}' \quad (1.7)$</p><p>As our steady state is by definition not changing, we can analyse the system by how this perturbation varies </p><p>$\frac{d\mathbf{x}'}{dt} = \mathcal{F}(\mathbf{x}_0 + \mathbf{x}'; p) \quad (1.8)$</p><p></p><h4 id="7e237f6a-0706-450a-b858-7c821cb310b4" data-toc-id="7e237f6a-0706-450a-b858-7c821cb310b4" collapsed="false" seolevelmigrated="true">5. Linearisation</h4><p>By considering only very small perturbations we are able to linearise our governing equations </p><p>$\mathcal{F}(x_0 + x'; p) = \underbrace{\mathcal{F}(x_0; p)}_{= 0} + \underbrace{\frac{\delta \mathcal{F}}{\delta x}(x_0; p) x'}_{\text{Linear Term}} + \text{Higher Order Terms (H.O.T.)}$</p><ul><li><p>Slight caveat of our governing equations must be smooth (differentiable) but this is true for most systems</p></li></ul><p></p><h4 id="d26c3f97-efb9-43b3-b718-22c704baefce" data-toc-id="d26c3f97-efb9-43b3-b718-22c704baefce" collapsed="false" seolevelmigrated="true">6. Non dimensionalising</h4><p>This is an optional step, but often helps understand the system</p><ul><li><p>Such as normalising with a length scale, or finding a time scale</p></li><li><p>Helps find non dimensional parameters that define the system like <strong>reynolds numbers, rayleigh number</strong> </p></li></ul><p></p><h4 id="69aa0805-8626-469e-bf73-889d73f9e03c" data-toc-id="69aa0805-8626-469e-bf73-889d73f9e03c" collapsed="false" seolevelmigrated="true">7. Normal modes and exponential growth</h4><ul><li><p>This is because we have linearised our system so superposition applies</p></li><li><p>Because of the time invariance an exponential solution must be apply,</p></li><li><p>As such we have an eigenvalue problem,$s\hat{\mathbf{x}} = \mathcal{L}\hat{\mathbf{x}}$</p></li></ul><p></p><h4 id="c052cc37-4153-4335-bbdf-d85c15351647" data-toc-id="c052cc37-4153-4335-bbdf-d85c15351647" collapsed="false" seolevelmigrated="true">8. Conditions for instability</h4><ul><li><p>This is basically if we have any poles s, where have positive real values</p></li></ul><table style="min-width: 75px;"><colgroup><col style="min-width: 25px;"><col style="min-width: 25px;"><col style="min-width: 25px;"></colgroup><tbody><tr><td colspan="1" rowspan="1" style="animation: auto ease 0s 1 normal none running none; appearance: none; background: none 0$\Re(s_1) < 0$</span></p></td></tr><tr><td colspan="1" rowspan="1" style="animation: auto ease 0s 1 normal none running none; appearance: none; background: none 0$\Re(s_1) = 0$</span></p></td></tr><tr><td colspan="1" rowspan="1" style="animation: auto ease 0s 1 normal none running none; appearance: none; background: none 0$\Re(s_1) > 0$$

Example of a structural system : Mass spring damper

Modelling equations

This is described by a second order differential equation

$m\ddot{x} + b\dot{x} + kx = f$ (where f is a driving function). This is already linear

We can also take an energy perspective, multiplying the force by the velocity to get the work.

• so for this is it $\frac{d}{dt}\left[\frac{1}{2}m\left(\frac{dx}{dt}\right)^2+\frac{1}{2}kx^2\right]=-b\left(\frac{dx}{dt}\right)^2+f\frac{dx}{dt}$

• We have our conservative KE and PE terms on LHS. Damper and energy addition via external work too.

Matrix form

For two degrees of freedom (with no damping) say if we have a system in the form: (for pitch and heave)

• $$m\ddot{h}+kh+kx_{s}\sin\theta=L$$

• $$I\ddot{\theta}+kx_{s}h\cos\theta+\kappa\theta=Lx_{a}\cos\theta$$

we can rewrite this coupled equation in terms of a mass and stiffness matrix. where the state vector $\mathbf{x} = [h, \theta]^T$

$\mathbf{M}\ddot{\mathbf{x}} + \mathbf{Kx} = 0$ where $\mathbf{M} = \begin{pmatrix} m &amp; 0 \\ 0 &amp; I \end{pmatrix}$ and $\mathbf{K} = \begin{pmatrix} k &amp; kx_s \\ kx_s &amp; \kappa \end{pmatrix}$

energy formulation

Our kinetic and potential energy of this system can be written as the below, where we have energy conservation.

$$\frac{d}{dt}\left[\frac{1}{2}\dot{\mathbf{x}}^{T}\mathbf{M}\dot{\mathbf{x}}+\frac{1}{2}\mathbf{x}^{T}\mathbf{Kx}\right]=0$$

Eigenvector equation:

We want to solve the for eigenvalues which are our natural frequencies. As $x\ddot = -\omega² x$

• This is the determinant of $(K - \omega^2 M)\hat{x} = 0$

• We would then plug in our eigenvalues to find the eigenvectors which are our mode shapes

• Remember to note symmetric systems have orthogonal eigenvectors

Non dimensionalisation:

We can non dimensionlise our variables.

$$t=\tilde{t}\sqrt{\frac{m}{k}},\quad h^{\prime}=R_{g}\tilde{h},\quad\theta^{\prime}=\tilde{\theta},\quad\delta=\frac{2x_s}{R_g},\quad\text{and}\quad\omega_{a}^2=\frac{\kappa}{kR_g^2}$$

Basically dividing by our natural frequency of heave, radius of gyration, our ratio between heave to CoM and radius of gyration and our ratio of torsional to translational stiffness

Very large systems and rayleigh’s quotient

We can represent increasingly large systems with just a larger state vector.

$$\mathbf{M}\frac{d^{2}\mathbf{x}'}{dt^{2}}+\mathbf{Kx}^{\prime}=0$$

We can represent any response as a linear combination of our eigenvector modes, as a modal summation

$$\sum_{i=1}^{n}(\ddot{a}_{i}+\omega_{i}^2a_{i})\mathbf{M}\hat{\mathbf{x}}_{i}=\mathbf{f}$$

As these modes are orthogonal, if we dot this with an eigenvector $\hat{\mathbf{x}}_j$, it only leaves the i = j term.

Hence: $(\ddot{a}_j + \omega_j^2 a_j) (\hat{\mathbf{x}}_j^T \mathbf{M} \hat{\mathbf{x}}_j) = \hat{\mathbf{x}}_j^T \mathbf{f}$

and: $\ddot{a}_{j}+\omega_{j}^2a_{j}=\frac{\hat{\mathbf{x}}_j^T \mathbf{f}}{\hat{\mathbf{x}}_j^T \mathbf{M} \hat{\mathbf{x}}_j}$

Example of a structural system : Continuous string system

This is now a continuous system, with a stretched string or membrane.

Derivation of response

if we consider a free small element, our vertical force balance

• $(\mu \Delta x) \frac{\partial^2 y}{\partial t^2} = T \left[ \left( \frac{\partial y}{\partial x} \right)_{x+\Delta x} - \left( \frac{\partial y}{\partial x} \right)_x \right] + f \Delta x$ plugging into f = ma, we get

• $\mu \frac{\partial^2 y}{\partial t^2} = T \left[ \frac{\left( \frac{\partial y}{\partial x} \right)_{x+\Delta x} - \left( \frac{\partial y}{\partial x} \right)_x}{\Delta x} \right] + f$ this is basically the definition of a derivative, as Δx → 0

• $$\mu\frac{\partial^2 y}{\partial t^2}=T\frac{\partial^2 y}{\partial x^2}+f$$

Result

This is a PDE with a separable MODAL solution $y^{\prime}(x,t)=\sum_{n=1}^{\infty}a_{n}(t)\hat{y}_{n}(\tilde{x})$

Our modes will take a harmonic solution:

$$\hat{y}(\tilde{x}) = A \sin(\omega \tilde{x}) + B \cos(\omega \tilde{x}) \text{}$$

• For this stretched string example, plugging in our boundary conditions we get $\hat{y}_{n}(\tilde{x})=\sin(n\pi\tilde{x})$

Frequency response function

we are then calculating our frequency response function.

$$\sum_{n=1}^{\infty}(\ddot{a}_{n}+\omega_{n}^2a_{n})\hat{y}_{n}(\tilde{x})=\tilde{f}(\tilde{x},\tilde{t})$$

Now taking the integral dot product with respect to modes

$$\int_0^1\hat{y}_{m}(\tilde{x})\left[\sum_{n=1}^{\infty}(\ddot{a}_{n}+\omega_{n}^2a_{n})\hat{y}_{n}(\tilde{x})\right]d\tilde{x}=\int_0^1\hat{y}_{m}(\tilde{x})\tilde{f}(\tilde{x},\tilde{t})d\tilde{x}$$

From our orthogonality condition, only the modes where n=m are non zero

$$\ddot{a}_{m}+\omega_{m}^2a_{m}=\frac{\int_{0}^{1} \hat{y}_m(\tilde{x}) \tilde{f}(\tilde{x}, \tilde{t}) d\tilde{x}}{\int_{0}^{1} \hat{y}_m(\tilde{x})^2 d\tilde{x}}$$

This probably makes more sense in a transfer function notation, where we will consider only our modal shapes.

$$a_{m}\left(\omega_{m}^2-\omega^2\right)_{}=\frac{\int_{0}^{1} \hat{y}_m(\tilde{x}) \tilde{f}(\tilde{x}, \tilde{t}) d\tilde{x}}{\int_{0}^{1} \hat{y}_m(\tilde{x})^2 d\tilde{x}}$$

now combining our mass normalisation in the denominator with our modal dot product above, we recover the 3C6 equation

$H = \sum \frac{u_j^{(n)} u_k^{(n)}}{\omega_n^2 - \omega^2}$

Example of a structural system : Euler beam

This is a 4th order continuous system,

we will first derive it for the case with no tension

Derivation of response

Again if we consider a small element, this is basically a vertical force balance but with shear forces

• From classical beam theory our shear force is $V=\frac{\partial M}{\partial x}=B\frac{\partial^3 y}{\partial x^3} \quad \text(\space {where} \space B=EI)$

• Looking at a small element on both sides, again we’re looking at the difference in shear over Δx

• This results in our derivative term:

$$\mu\frac{\partial^2 y}{\partial t^{2}}+B\frac{\partial^4 y}{\partial x^{4}}=f$$

Considering the free response assuming a separable modal and harmonic solution $\hat{y} = exp(\omega t)exp(kx)$

• $$\mu\frac{\partial^2 y}{\partial t^{2}}+B\frac{\partial^4 y}{\partial x^{4}}=0$$

• $$Bk^4 \hat{y}+ \mu\omega^2\hat{y}=0$$

• now defining k as wave number of our modes, given we need an exponential modal solution

$$\omega=k^2\sqrt{\frac{B}{\mu}}=k^2\sqrt{\frac{EI}{\mu}}$$

Our full equation if we include tension as well, is:

$$\mu\frac{\partial^{2}y}{\partial t^{2}}+B\frac{\partial^{4}y}{\partial x^{4}}-T\frac{\partial^{2}y}{\partial x^{2}}=f$$

Boundary conditions for beams

We have a characteristic equation for k in the form of k⁴ = β, as such we have four roots on the real and imaginary axes. This results in the solution:

$\hat{y}(x) = C_1 e^{kx} + C_2 e^{-kx} + C_3 e^{ikx} + C_4 e^{-ikx}$ which can be written much more neatly as:

$$\hat{y}(x) = A \sin(kx) + B \cos(kx) + C \sinh(kx) + D \cosh(kx)$$

• This is true with tension too as this just super imposes more harmonic solutions

To find our boundary conditions, we can basically differentiate this, we know:

• No end translation (pinned or fixed): y = 0

• Fixed end (no rotation): dy/dx = 0

• no shear (free end): d²y/dx² = 0

• no moment (free end/pinned end): d³y/dx³ = 0

From this we get expressions that can help define our mode shapes, such as the 3C6 databook expressions:
$\cos \alpha \cosh \alpha + 1 = 0, \quad \alpha_1 = 1.8751$

$$\cos \alpha \cosh \alpha - 1 = 0, \quad \alpha_1 = 4.7300$$

$$\tan \alpha - \tanh \alpha = 0, \quad \alpha_1 = 3.9266$$

Also we can take into account the symmetry of the beam to work out what functions exist or not

we can see the solutions for a fixed free beam.

Modelling fluids

Our fluids can generally be expressed with the navier stokes equal which models
$\rho\left[\frac{\partial\mathbf{u}}{\partial t}+(\mathbf{u}\cdot\nabla)\mathbf{u}\right]+\nabla p=\mathbf{f}+\mu\nabla^2\mathbf{u}$

• This is a momentum equation

• f is external volumetric force

Non dimensionalising NS

We can non dimensionalise the navier stokes equations, for turbulant flows relevant length scales are:

• $\mathbf{u} = U\tilde{\mathbf{u}}$ rescaling with a free stream velocity

• $\ell = \mu / (\rho U)$ this is basically assuming constant reynolds = 1, for length scaling. In turbulent regime flow isn’t very Re dependant

• $p = \rho U^2 \tilde{p}$ scale with dynamic pressure

• $\tilde{\mathbf{f}} = \mathbf{f} L / \rho U^2$ scale with dynamic pressure and length scale

thus we get:

$$\frac{L}{UT}\frac{\partial\tilde{\mathbf{u}}}{\partial\tilde{t}}+(\tilde{\mathbf{u}}\cdot\tilde{\nabla})\tilde{\mathbf{u}}+\tilde{\nabla}\tilde{p}=\tilde{\mathbf{f}}+Re^{-1}\tilde{\nabla}^2\tilde{\mathbf{u}}$$

Differences between low reynolds and high reynolds flows.

Low reynolds flow

Here viscosity completely dominates inertia, so our acceleration term and material derivative term can be removed, this leaves:

$-\nabla p+\mu\nabla^2\mathbf{u}=\mathbf{f}$ These are the stoke equations

• This is linear

• no du/dt term so it is instantaneous

• reversible (fluid can be unmixed)

High reynolds flow

this is generally the flow regime of our everyday life with big objects moving at high speed.

• Inertial forces dominate over viscous forces and flow can almost be considered inviscid

• Except in a boundary layer where viscous effects dominate

We can model many important results in terms of inviscid potential flow + a boundary layer effect

High reynolds example: Thin airfoil lift

An example of high reynolds flow is thin airfoil lift

• This is completely modelled with potential flow theory, where we are examining flow over a spinning cylinder

• Now remapping this flow where the cylinder turns into a flat plate using a joukowski transform

• The vortex we are using to model contains the vorticity that would be carried within a boundary layer in the real world.

This leads to the result:

$$F_{L}=\frac{1}{2}\rho U^2cC_{L}(\alpha,Re),\quad\text{where}\quad C_{L}(\alpha,Re)\approx2\pi\alpha,$$

This predicts no drag (D’alembert paradox) boundary layer theory is needed to model this, with drag being shown to be

$$F_{D}=\frac{1}{2}\rho U^2cC_{D}$$

High reynolds example: Bluff body

1. The Breakdown of Potential Flow

In potential flow theory, a fluid is assumed to be inviscid and "slip" perfectly around a body, meeting up again at the back to recover the pressure. This leads to D'Alembert's Paradox, where the math predicts zero drag because the front and back pressures cancel out.

• Separation: For bluff bodies, the high Reynolds number flow cannot handle the "climb" against the increasing pressure on the back half of the shape.

• Vorticity Shedding: The boundary layer detaches (separates), shedding vorticity into the region behind the body.

• The Wake: This separation creates a wake—a region of low-pressure, turbulent fluid that prevents the pressure recovery predicted by potential flow.

The wake is very unsteady, but we can average out the unsteadiness and define an average drag and lift:

$$\overline{F}_{D}=\frac{1}{2}\rho U^2A\overline{C}_{D}$$

$$\overline{F}_{L}=\frac{1}{2}\rho U^2A\overline{C}_{L}$$

Discussion of statistical steadiness

The parameter $UT/L$, defines how quickly our system changes. If it is large, it means our system moves slowly relative to the flow field

• However: if it is small this our system moves at a similar speed to the flow field, no more spectral separation. this is the realm of our vortex induced vibrations.

Simplified NS example: Unsteady fluid carrying conduit

This is an example of how we can simplify the complex 3d navier stokes equations to 1D.

Physical explanation

We can see a flexible pipe, carrying a fluid with constant density and a constant velocity across the crosssection

• The centreline of the conduit is defined by X(s,t)

Finding velocity

• Our velocity is at magnitude U and in the same direction as X, hence our velocity is the material derivative of X

• $$\mathbf{u}=\frac{\partial\mathbf{X}}{\partial t}+U\frac{\partial\mathbf{X}}{\partial s}=\left(\frac{\partial}{\partial t}+U\frac{\partial}{\partial s}\right)\mathbf{X}$$

Finding acceleration

• Our fluid acceleration is defined as the material derivative of velocity

• $$\mathbf{a} = \left( \frac{\partial}{\partial t} + U \frac{\partial}{\partial s} \right) \mathbf{u} = \left( \frac{\partial}{\partial t} + U \frac{\partial}{\partial s} \right)^2 \mathbf{X}$$

Force on pipe

Now because the crosssection is uniform, and the fluid is inviscid, the pressure is constant. As such the only force on the pipe is from the fluid acceleration.

$$\mathbf{F}=-\rho A\mathbf{a}=-\rho A\left(\frac{\partial}{\partial t}+U\frac{\partial}{\partial s}\right)^2\mathbf{X}$$

Examples

now considering two examples

Semi circular conduit

here we are considering a steady conduit with square side b, it forms a quarter circle with radius R

• There is a centrifugal force on the conduit at the 45 degree direction

• Can solve easily by considering inlet and outlet momentums

  • If you hate yourself you can do the integral of force along the length

Flexible conduit

Here we have a pipe which the centreline is varying sinusoidally given by $\mathbf{X}=s\hat{e}_{x}+a\Re(e^{i\omega t-iks})\hat{e}_{y}$

Using our force function, our force per unit length is $\mathbf{F}=\rho Aa(\omega+kU)^2\Re(e^{i\omega t-iks})\hat{e}_{y}$

• Can see that we have an amplifying force that can result in instability

Simplified NS example: Thin layer flow

Now we are considering a thin layer of fluid which is viscous but incompressible. Our goal is to find the fluid pressure to find the force on the wall.

Mass continuity

First our fluid is incompressible so we can add teh condition: $\frac{\partial h}{\partial t}+\frac{\partial q}{\partial x}=0,\quad\text{where}\quad q=\int_0^{h}u(x,y,t)\,dy$

Considering extreme viscosity cases

We can consider two limits for our navier stokes equations, our euler equations which are inviscid and our viscous dominated stoke’s equations

Viscous case

Now our fluid is dominated by viscosity, so our navier stokes reduces down.

• $\rho\left[\frac{\partial\mathbf{u}}{\partial t}+(\mathbf{u}\cdot\nabla)\mathbf{u}\right]+\nabla p=\mathbf{f}+\mu\nabla^2\mathbf{u}$ (no momentum transport terms and no external force), we are left with only pressure and momentum diffusion

• $$\nabla p=\mu\nabla^2\mathbf{u}$$

in 1D this is: $\frac{\partial p}{\partial x}=\mu\frac{\partial^2 u}{\partial y^2}\implies u=-\frac{1}{2\mu}\frac{\partial p}{\partial x}y(h-y)$

this is our poiseuille profile

Inviscid case

Now we neglect our viscous terms so navier stokes goes to the euler equations

• $$\rho\left[\frac{\partial\mathbf{u}}{\partial t}+(\mathbf{u}\cdot\nabla)\mathbf{u}\right]+\nabla p=0$$

We need to guess the velocity profile, a common guess is a top hat profile ie $q = uh$.

We can use euler equations directly to find our pressure (or unsteady bernoulli which is derived from euler)

Example, onset of roll waves.

We have a thin film of fluid rapidly flowing down at an incline θ, this forms an undulation at the surface. Determine the onset criteria for these undulations

Assumptions

• Since it says rapidly, we can approximate our fluid as inviscid

• We will take wall drag as the Chezy friction coefficient: $d(u,h)=C_{f}\rho\frac{u^2}{h}$

plugging into governing equations

• Continuity transforms to: $\frac{\partial h}{\partial t}+\frac{\partial(uh)}{\partial x}=0$

• Navier stokes goes to: $\frac{\partial u}{\partial t}+u\frac{\partial u}{\partial x}+g\cos\theta\frac{\partial h}{\partial x}=g\sin\theta-C_{f}\frac{u^2}{h}$

Linear stability analyisis

Now lets consider our perturbations:

• $u = U + u’$ and $h = H + h’$

plugging this into navier stokes and removing small terms, we get:

$$\frac{\partial h'}{\partial t}+U\frac{\partial h'}{\partial x}+H\frac{\partial u'}{\partial x}=0$$

$$\frac{\partial u'}{\partial t}+U\frac{\partial u'}{\partial x}+g\cos\theta\frac{\partial h'}{\partial x}=C_{f}\left(\frac{2U}{H}u^{\prime}-\frac{U^2}{H^2}h^{\prime}\right)$$

Now solving for normal mode solutions in form $(u', h') = (\hat{u}, \hat{h}) e^{st} e^{ikx}$

we get:

• $$s\hat{h} + ikU\hat{h} + ikH\hat{u} = 0$$

• $$s\hat{u}+ikU\hat{u}+g\cos\theta ik\hat{h}=g\sin\theta\left(2\frac{\hat{u}}{U}-\frac{\hat{h}}{H}\right)$$

At the threshold our real component is zero, so we are on the imaginary axis

• looking at the real and imaginary components of the equation

• we get the threshold:

$$\frac{U^2}{gH \cos\theta}=2\quad\text{or, equivalently,}\quad\tan\theta=2C_{f}$$

This is the square of the froude number, showing how this non-dimensional quantity is relevant for fluid waves.

our stability is when we have a breakdown of friction as shown by the tanθ relationship