Preparation

trailer Cellules de Bénard (1min)

Rayleigh–Bénard convection: cooking oil and small aluminium particles (5 min)

Content for today

Material laws

The effect of stress represented by \(\nabla p\) and \(\nabla \cdot\mathbb{T}\) terms.

\(\nabla p\) is called the pressure gradient and arises from the isotropic part of the stress tensor.

The anisotropic part of the stress tensor gives rise to \(\nabla \cdot\mathbb{T},\) describing viscous forces.

The stress terms p and \(\mathbb{T}\) are yet unknown: a force model is needed relating the stresses to the fluid motion.

Assumptions on the specific behavior of a fluid are made (based on observations).

The Cauchy stress tensor can be also written in matrix form: \[ \mathbb{T} = \left({\begin{matrix} \mathbf{T}^{(\mathbf{e}_1)} \\ \mathbf{T}^{(\mathbf{e}_2)} \\ \mathbf{T}^{(\mathbf{e}_3)} \\ \end{matrix}}\right) = \left({\begin{matrix} \sigma _{11} & \sigma _{12} & \sigma _{13} \\ \sigma _{21} & \sigma _{22} & \sigma _{23} \\ \sigma _{31} & \sigma _{32} & \sigma _{33} \\ \end{matrix}}\right) \equiv \left({\begin{matrix} \sigma _{xx} & \sigma _{xy} & \sigma _{xz} \\ \sigma _{yx} & \sigma _{yy} & \sigma _{yz} \\ \sigma _{zx} & \sigma _{zy} & \sigma _{zz} \\ \end{matrix}}\right) \equiv \left({\begin{matrix} \sigma _x & \tau _{xy} & \tau _{xz} \\ \tau _{yx} & \sigma _y & \tau _{yz} \\ \tau _{zx} & \tau _{zy} & \sigma _z \\ \end{matrix}}\right) \label{eqn:T} \]

where \(\sigma\) are the normal stresses and \(\tau\) are the shear stresses. From the Newton’s third law (actio est reactio), the stress tensor is symmetrical: \[ \mathbb{T} = \mathbb{T}^T \]

In the Figure the stress vectors \(\mathbf{T}^{(\mathbf{e}_{i})}\) can be decomposed in one normal stress and two shear stress components.

Components of stress in three dimensions.

Conservation of mass

Regardless of the flow assumptions, a statement of the conservation of mass is generally necessary. This is achieved through the mass continuity equation, given in its most general form as: \[ \frac{\partial \rho}{\partial t} + \nabla \cdot (\rho \mathbf{u}) = 0 \] or, using the substantive derivative: \[ \frac{D\rho}{Dt} + \rho (\nabla \cdot \mathbf{u}) = 0. \]

Navier-Stokes equations

(for the details: script)

The shear stress term \(\nabla \mathbb{T}\) becomes the quantity \(\mu \nabla^2 \mathbf{u}\), where \(\mu\) is the dynamic viscosity.

\[ \rho \left(\frac{\partial \mathbf{u}}{\partial t} + \mathbf{u} \cdot \nabla \mathbf{u}\right) = -\nabla p + \mu \nabla^2 \mathbf{u} + \mathbf{F}. \label{reyn0} \]

The vector field \(\mathbf{F}\) represents “other” (body force) forces. Typically gravity, but may include other fields (such as electromagnetic).

Often, these forces may be represented as the gradient of some scalar quantity. E.g.: Gravity in the z direction is the gradient of \(- \rho gz.\)

Since pressure shows up only as a gradient, this implies that solving a problem without any such body force can be mended to include the body force by modifying pressure.

It’s well worth observing the meaning of each term: \[ \overbrace{\rho \Big( \underbrace{\frac{\partial \mathbf{u}}{\partial t}}_{ \begin{smallmatrix} \text{Unsteady}\\ \text{acceleration} \end{smallmatrix}} + \underbrace{\mathbf{u} \cdot \nabla \mathbf{u}}_{ \begin{smallmatrix} \text{Advective} \\ \text{acceleration} \end{smallmatrix}}\Big)}^{\text{Inertia (per volume)}} = \overbrace{\underbrace{-\nabla p}_{ \begin{smallmatrix} \text{Pressure} \\ \text{gradient} \end{smallmatrix}} + \underbrace{\mu \nabla^2 \mathbf{u}}_{\text{Viscosity}}}^{\text{Divergence of stress}} + \underbrace{\mathbf{F}}_{ \begin{smallmatrix} \text{Other} \\ \text{body} \\ \text{forces} \end{smallmatrix}}. \label{gleichmitKlammer} \]

If temperature effects are also neglected, the only “other” equation: mass continuity equation. Under the incompressible assumption, density is a constant and it follows that the equation will simplify to: \[ \nabla \cdot \mathbf{u} = 0 \qquad . \label{consmass} \] This is more specifically a statement of the conservation of volume !

Exercise advection

The temperature at a point 50 km north of a station is 3\(^\circ\)C cooler than at the station. If the wind is blowing from the northeast at 20m/s and the air is being heated by radiation at a rate of 1\(^\circ\)C/h, what is the local temperature change at the station?

Elimination of the pressure term

in the 2D-Navier-Stokes equation.

2D: assume \(w = 0\) and no dependence of anything on z

\[ \rho \left(\frac{\partial u}{\partial t} + u \frac{\partial u}{\partial x} + v \frac{\partial u}{\partial y}\right) = -\frac{\partial p}{\partial x} + \mu \left(\frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2}\right) \] \[ \rho \left(\frac{\partial v}{\partial t} + u \frac{\partial v}{\partial x} + v \frac{\partial v}{\partial y}\right) = -\frac{\partial p}{\partial y} + \mu \left(\frac{\partial^2 v}{\partial x^2} + \frac{\partial^2 v}{\partial y^2}\right) \]

or using the total derivative

Non-dimensional parameters: The Reynolds number

For the case of an incompressible flow in the Navier-Stokes equations, assuming the temperature effects are negligible and external forces are neglected.

conservation of mass \[ \nabla \cdot \mathbf{u} = 0 \] conservation of momentum \[ \partial_t \mathbf{u} + ( \mathbf{u} \cdot \nabla) \mathbf{u} = - \frac{1}{\rho_0} \nabla p + \nu \nabla^2 \mathbf{u} \]

The equations can be made dimensionless by a length-scale L, determined by the geometry of the flow, and by a characteristic velocity U.

For analytical solutions, numerical results, and experimental measurements, it is useful to report the results in a dimensionless system (concept of dynamic similarity).

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)

cartoon

clouds

sketch

Bifurcations

Bifurcation youtube, Bifurcation Khan academy, Bifurcation K

Max and Moritz, not lazy, sawing secretly a gap in the bridge.
When now this act is over, you suddenly hear a scream:
“Hey, out! You billy goat!
Tailor, tailor, bitch, bitch, bitch!”

link

Bridge of Schneider Böck: Stability - Instability

And there he is on the bridge,
Cracks! The bridge is falling apart;

Bild

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.

Stability - Instability: Consumer-producer problem

A typical example could be the consumer-producer problem where the consumption is proportional to the (quantity of) resource.

For example:

\[ \frac{dx}{dt}=rx(1-x)-px \]

where

\[ rx(1-x) \]

is the logistic equation of resource growth:

Rate of reproduction proportional to the population, available resources

\(px\): Consumption, proportional to the resource x.

\[ x_{E 1} = 0, \mbox{ and } x_{E 2 } = 1 - \frac{p}{r} \]

Logistic equation of population growth

Verhulst: describe the self-limiting growth of a biological population with size N:

\[ \frac{dN}{dt}=r N \cdot \left( 1- \frac{N}{K} \right) \]

r growth rate and K carrying capacity.

the early, unimpeded growth rate is modeled by the first term \(r N.\)

Human population: logistic growth model

Bild

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

Saddle-node bifurcation: two fixed points collide

4) Graphical method: slope at equilibrium points

\[ \frac{dx}{dt}=b+x^2 \]

Bild

Convection in the Rayleigh-Benard system

Bild

Rayleigh (1916) temperature difference between the upper- and lower-surfaces \[ T(x, y, z=H) = \, T_0 \] \[ T(x, y, z=0) \, = \, T_0 + \Delta T \]

Furthermore \[ \rho = \rho_0 = const. \] except in the buoyancy term, where:

\[ \varrho = \varrho_0 (1 - \alpha(T-T_0)) \mbox{ with } \alpha > 0 \quad . \]

common feature of geophysical flows

No Convection Equilibrium: Diffusion

Bild

Diffusion: Temperature varies linearly with depth:

\[ T_{eq} = T_0 + \left(1 - \frac{z}{H}\right) \Delta T \]

No movement of particles:

\[ u = w= 0 \]

When this solution becomes unstable, convection should develop.

No Convection Equilibrium: Diffusion

Rayleigh-Bénard convection and bifurcation

Vorticity in the Rayleigh-Benard system

\[ 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) \]

Temperature in the Rayleigh-Benard system

\[ T_{eq} = T_0 + \left(1 - \frac{z}{H}\right) \Delta T \]

\[ \mbox{with } \quad T = T_{eq} + \Theta \quad , \quad \mbox{where } \quad \Theta \quad \mbox{is the anomaly} \]

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} \]

gives:

\[ 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 . \]

Galerkin approximation: Get a low-order model

Bild

\[ \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) \]

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 } \quad R_c = \pi^4 \frac{(1+a^2)^3}{a^2} \] \[=657.51 \mbox{ occurs when } 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.} \]

Bild

Lorenz system (also next lecture)

Bifurcation at \[ r = R_a/R_c = 1\]

Geometry constant \[b = 4(1+a^2)^{-1}\]

 
Famous low-order model:

\[ \dot X = -\sigma X + \sigma Y \]

\[ \dot Y = r X - Y - X Z \]

\[ \dot Z = -b Z + X Y \]

\[\mbox{dimensionless time } \quad t_d = \pi^2 H^{-2} (1+a^2) \kappa t,\]

\[ \mbox{ Prandtl number } \quad \sigma = \nu \kappa^{-1}, \]

Lorenz system r=24

left: 52%

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

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=3.5

r=3.5
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:1000){
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")

as an outlook: Conveyor Belt

Conveyor

outlook: Thermohaline ocean circulation

Overturning

Modelled meridional overturning streamfunction in Sv 10^6 = m^3 /s in the Atlantic Ocean. Grey areas represent zonally integrated smoothed bathymetry