Dynamics II
Lecture: May 5, 2024 (Monday), 14:00 Prof. Dr. Gerrit Lohmann
Rayleigh-Benard convection (and Rayleigh number)
Experiments:
Rayleigh–Bénard convection: cooking oil and small aluminium particles (5 min),
Rayleigh Benard
Thermal Convection 3D Simulation (2 min)
Linear stability analysis
Consider the continuous dynamical system described by \[ \dot x=f(x,\lambda)\quad \] A bifurcation occurs at \[(x_E,\lambda_0)\] if the Jacobian matrix \[ \textrm{d}f/dx (x_E,\lambda_0)\] has an Eigenvalue with zero real part.
Example: transcritical bifurcation
a fixed point interchanges its stability with another fixed point as the control parameter is varied. Bifurcation at \(r=0\).
\[ \frac{dx}{dt}=rx (1-x) \, \]
The two fixed points are 0 and 1. When r is negative, the fixed point at 0 is stable and 1 is unstable. But for \(r>0\), 0 is unstable and 1 is stable.
Saddle-node bifurcation: two fixed points collide
The normal form: \[ \frac{dx}{dt}=r+x^2 \]
\[ r<0: \mbox{ a stable equilibrium point at } -\sqrt{-r} \mbox{ and an unstable one at } \sqrt{-r} \qquad \qquad \]
\(r=0\): exactly one equilibrium point, saddle-node fixed point.
\(r>0\): no equilibrium points. Saddle-node bifurcations may be associated with hysteresis loops.
Saddle-node bifurcation: two fixed points collide
1) Calculate the equilibrium points:
\[ f(x) = r + x^2 = 0 \]
\[ x_{E1,E2} = \pm \sqrt{-r} \mbox{ for } r \le 0 \]
2) Linear stability conditions for the equilibrium points
\[ f'(x) = 2x \]
\[ f'(x_{E1}) = 2\sqrt{-r} > 0 \quad \mbox{unstable} \]
\[ f'(x_{E2}) = -2\sqrt{-r} < 0 \quad \mbox{stable} \]
For r=0: equilibrium point \(x_E=0\) which is indifferent (not stable/unstable).
3) Potential or Lyapunov Method
\[ \frac{dx}{dt}=b+x^2 = - \frac{d}{dx} \left( - b x - \frac{x^3}{3} \right) = - \frac{d}{dx} V (x) \]
Global analysis including basins of attraction for \(x_{E2}: (-\infty,3)\)
\[ x_{E1} = 2\sqrt{9} = 6 \quad \mbox{unstable} \]
\[ x_{E1} = -2\sqrt{9} = -6 \quad \mbox{stable} \]
4) Graphical method: slope at equilibrium points
\[ \frac{dx}{dt}=b+x^2 \]
filled points: positive slope => unstable
open points: negative slope => stable
No Convection Equilibrium: Diffusion
\[ D_t u = - \frac{1}{\rho_0} \partial_x p + \nu \nabla^2 u \label{eqref:einse} \qquad \qquad \qquad \qquad \qquad \qquad (1) \]
\[ D_t w = - \frac{1}{\rho_0} \partial_z p + \nu \nabla^2 w + g (1- \alpha (T-T_0)) \qquad \qquad (2)\]
\[ \partial_x u + \partial_z w = 0 \qquad \qquad \qquad \qquad \qquad \qquad \qquad (3) \]
\[ D_t T = \kappa \nabla^2 T \qquad \qquad \qquad \qquad \qquad \qquad \qquad (4) \]
Diffusion: Temperature varies linearly with depth:
\[ T_{eq} = T_0 + \left(1 - \frac{z}{H}\right) \Delta T \] fullfills (4)
No movement of particles:
\[ u = w= 0 \] fullfills (1,2,3)
with the hydrostatic vertical pressure variation
\[ \partial_z p = g \, \rho = g \, \rho_0 \left(1- \alpha (1 - \frac{z}{H}) \cdot \Delta T \right) \qquad \qquad (2)\]
No Convection.
When this solution becomes unstable, convection should develop. Spontaneous symmetry breaking.
Vorticity in the Rayleigh-Benard system
Introduce vorticity to have (1,2,3) in one equation:
\[ D_t \left( \nabla^2 \Psi\right) = \nu \nabla^4 \Psi - g \alpha \frac{\partial \Theta}{\partial x} \]
\[ \mbox{with } \quad T = T_{eq} + \Theta \quad , \quad \mbox{where } \quad \Theta \quad \mbox{is the anomaly from the diffusion solution.} \]
In vorticity equation:
\[ \frac{\partial }{\partial x} g (1- \alpha (T_{eq} + \Theta -T_0)) = - g \alpha \frac{\partial }{\partial x} \Theta \]
Temperature equation:
\[ D_t T = D_t T_{eq} + D_t \Theta = w \cdot \frac{- \Delta T}{H} + D_t \Theta = - \frac{\Delta T}{H} \frac{\partial \Psi}{\partial x} + D_t \Theta \]
\[ D_t \Theta \quad = \frac{\Delta T}{H} \frac{\partial \Psi}{\partial x} + \kappa \nabla^2 \Theta \quad \]
Non-dimensional Rayleigh-Benard system
\[ \frac{1}{T} \frac{1}{L^2} \frac{L^2}{T} D_t \left( \nabla_d^2 \Psi_d\right) = \nu \frac{1}{L^4} \frac{L^2}{T} \nabla_d^4 \Psi_d - g \alpha \frac{ \Delta T}{L} \frac{\partial \Theta_d}{\partial x_d} \]
\[ \frac{\Delta T}{T} D_t \Theta_d \quad = \frac{\Delta T}{H} \frac{L^2}{ T L} \frac{\partial \Psi_d}{\partial x_d} + \kappa \frac{\Delta T}{L^2} \nabla_d^2 \Theta_d \quad \]
This yields (remember \(L=H\)) \[ D_t \left( \nabla_d^2 \Psi_d\right) = \nu \frac{T}{H^2} \nabla_d^4 \Psi_d - g \alpha \frac{T^2 \Delta T}{H} \frac{\partial \Theta_d}{\partial x_d} \] \[ D_t \Theta_d \quad = \frac{\partial \Psi_d}{\partial x_d} + \kappa \frac{T}{H^2} \nabla_d^2 \Theta_d \quad \]
\[ \mbox{Inserting} \quad T= H^2/\kappa \] \[ \mbox{Rayleigh number} \quad R_a = \frac{g \alpha H^3 \Delta T}{\nu \kappa}\] \[ \mbox{Prandtl number} \quad \sigma = \frac{ \nu}{ \kappa} \]
\[ D_t \left( \nabla_d^2 \Psi_d\right) = \frac{ \nu}{ \kappa} \nabla_d^4 \Psi_d - g \alpha \frac{H^3 \Delta T}{\kappa^2} \frac{\partial \Theta_d}{\partial x_d} \label{eqref:psieqnn3} \] \[ D_t \Theta_d \quad = \frac{\partial \Psi_d}{\partial x_d} + \nabla_d^2 \Theta_d \quad . \]
Finally, we get the following dynamics:
\[ \underline{D_t \left( \nabla_d^2 \Psi_d\right) = \sigma \nabla_d^4 \Psi_d - R_a \sigma \frac{\partial \Theta_d}{\partial x_d} \quad} \]
\[ \underline{D_t \Theta_d \quad = \frac{\partial \Psi_d}{\partial x_d} + \nabla_d^2 \Theta_d \quad} \]
Galerkin approximation: Get a low-order model
\[ \mbox{ Saltzman (1962): Expand } \Psi, \Theta \mbox{ in double Fourier series in x and z: } \]
\[ \Psi (x,z,t) \, = \, \sum_{k=1}^\infty \sum_{l=1}^\infty \Psi_{k,l} (t) \, \, \sin \left(\frac{k \pi a}{H} x \right) \, \times \, \sin \left(\frac{ l \pi}{H} z \right) \] \[ \Theta (x,z,t) \, = \, \sum_{k=1}^\infty \sum_{l=1}^\infty \Theta_{k,l} (t) \cos \left(\frac{k \pi a}{H} x \right) \, \times \, \sin \left( \frac{l \pi}{H} z \right) \]
From partial differential equations to ordinary differential equations for \(\Psi_{k,l}(t)\) and \(\Theta_{k,l}(t)\).
Approximation: Just 3 Modes X(t), Y(t), Z(t)
\[ \frac{a}{1+a^2} \, \kappa \, \Psi = X \sqrt{2} \sin\left(\frac{\pi a}{H} x \right) \sin\left(\frac{\pi}{H} z \right) \]
\[ \pi \frac{R_a}{R_c} \frac{1}{\Delta T} \, \Theta = Y \sqrt{2} \cos\left(\frac{\pi a}{H} x\right) \sin\left(\frac{\pi}{H} z \right) - Z \sin\left(2 \frac{\pi}{H} z \right) \]
Rayleigh number Ra: Buoyancy & Viscosity
\[ \mbox{Motion develops if } \quad R_a = \frac{g \alpha H^3 \Delta T}{\nu \kappa} \quad \mbox{exceeds a critical } \quad R_c \]
As the Rayleigh number increases, the gravitational force becomes more dominant. The critical Rayleigh number can be obtained analytically for a number of different boundary conditions by doing a perturbation analysis on the linearized equations in the stable state.
\(R_c = \pi^4 \frac{\left(1+a^2\right)^3}{a^2} = \pi^4 \frac{27}{4}= 657.51\) with \(a^2 = 1/2\)
\[\mbox{When } R_a < R_c,\mbox{ heat transfer is due to conduction} \]
\[\mbox{When } R_a > R_c, \mbox{ heat transfer is due to convection.} \]
Lorenz system
Bifurcation at \[ r = R_a/R_c = 1\]
Famous low-order model:
\[ \dot X = -\sigma X + \sigma Y \]
\[ \dot Y = r X - Y - X Z \]
\[ \dot Z = -b Z + X Y \]
Lorenz system r=0.9
r=0.9
s=10
b=8/3
dt=0.01
x=1.1
y=0.1
z=11.1
vx<-c(0)
vy<-c(0)
vz<-c(0)
for(i in 1:100){
x1=x+s*(y-x)*dt
y1=y+(r*x-y-x*z)*dt
z1=z+(x*y-b*z)*dt
vx[i]=x1
vy[i]=y1
vz[i]=z1
x=x1
y=y1
z=z1}
plot(vx,type="l",xlab="time",ylab="x")
plot(vy,type="l",xlab="time",ylab="y")
Lorenz system r=24
r=24
s=10
b=8/3
dt=0.01
x=0.1
y=0.1
z=0.1
vx<-c(0)
vy<-c(0)
vz<-c(0)
for(i in 1:10000){
x1=x+s*(y-x)*dt
y1=y+(r*x-y-x*z)*dt
z1=z+(x*y-b*z)*dt
vx[i]=x1
vy[i]=y1
vz[i]=z1
x=x1
y=y1
z=z1}
plot(vx,vy,type="l",xlab="x",ylab="y")
plot(vy,vz,type="l",xlab="y",ylab="z")
North Atlantic Current & Gulfstream
brings warm water northward where it cools.
returns southward as a cold, deep, western-boundary current.
Gulf Stream carries 40 Sv of 18°C water northward.
Of this, 15 Sv return southward in the deep western boundary current at a temperature of 2°C.
How much heat is transported northward ?
Calculation:
\[ \underbrace{ c_p}_{4.2 \cdot 10^3 Ws/(m^3 kg)} \, \cdot \, \underbrace{ \rho }_{10^3 kg/m^3 } \, \cdot \, \underbrace{\Phi}_{15 \cdot 10^6 m^3/s} \, \cdot \, \underbrace{\Delta T}_{(18-2) K } = 1 \cdot 10^{15} W \]
The flow carried by the conveyor belt loses 1 Petawatts (PW), close to estimates of Rintoul and Wunsch (1991)
The deep bottom water from the North Atlantic is mixed upward in other regions and ocean, and it makes its way back to the Gulf Stream and the North Atlantic. Thus most of the water that sinks in the North Atlantic must be replaced by water from the far South Atlantic and Pacific Ocean.
Ocean Conveyor Belt
Conveyor Belt: Industry
Conveyor belt circulation
The the conveyor is driven by deepwater formation in the northern North Atlantic.
The conveyor belt metaphor necessarily simplifies the ocean system, it is of course not a full description of the deep ocean circulation.
Broecker’s concept provides a successful approach for global ocean circulation, although several features can be wrong like the missing Antarctic bottom water, the upwelling areas etc..
metaphor inspired new ideas of halting or reversing the ocean circulation and put it into a global climate context.
interpretation of Greenland ice core records indicating different climate states with different ocean modes of operation (like on and off states of a mechanical maschine).
Thermohaline ocean circulation
Modelled meridional overturning streamfunction in Sv 10^6 = m^3 /s in the Atlantic Ocean. Grey areas represent zonally integrated smoothed bathymetry
Estimates of overturning ?
It is observed that water sinks in to the deep ocean in polar regions of the Atlantic basin at a rate of 15 Sv. (Atlantic basin: 80,000,000 km^2 area * 4 km depth.)
– How long would it take to ‘fill up’ the Atlantic basin?
– Supposing that the local sinking is balanced by large-scale upwelling, estimate the strength of this upwelling.
Hint: Upwelling = area * w
– Compare this number with that of the Ekman pumping!
Estimates of overturning: Solution
Timescale T to ‘fill up’ the Atlantic basin:
\[ T = \frac{ 80 \cdot 10^{12} \, m^2 \cdot 4000 \, m}{15 \cdot 10^6 \, m^3 s^{-1}} = 2.13 \cdot 10^{10} s = 676 \;years\]
Overturning is balanced by large-scale upwelling:
\[ area \cdot w = 15 \cdot 10^6 \, m^3 s^{-1}\]
\[ w = 0.1875 \cdot 10^{-6} m\;s^{-1} = 5.9 \cdot 10^{-15} m \; y^{-1}. \]
Ekman pumping \[ w_E \simeq 32 \, \, m \; y^{-1}. \]