ASTR/ATOC-5410: Fluid Instabilities, Waves, and Turbulence September 12, 2016, Axel Brandenburg

Handout 3: Rayleigh–Bénard problem

Convection plays an important role in geo- and astrophysics. Without it, the would be no turbulence and no magnetic fields in the Sun, for example. This problem has been studied for a long time and goes back to early work of Count Rumford. In a historical account, Brown (1957) mentions in this connection the year 1797,

The decisive statement about ascending currents is made here:

But the word convection cas coined only much later.

Figure 1: Convec- tion experiment. A detailed account of the subsequent discoveries that went into our theoretical understanding of this phenomenon is given by Chandrasekhar (1961). After the original experimental work of Bénard (1900) and theoretical work of Rayleigh (1916), the mathematical formulation that is used nowadays in many textbooks goes back to work by Jeffreys (1926). It was also Jeffreys (1930), who generalized the theory to the compressible case, where the temperature gradient is to be replaced by the superadiabatic gradient, which is the famous Schwarzschild (1906) criterion.

Discrepancies between Bénard’s original experiments and the theoretical calculations are explained by the importance of surface tension, which was ignored in the early work. This leads to the so-called Marangoni effect which describes the mass transfer along an interface between two fluids due to surface tension gradient.

1 Governing equations

The best studied case is that of very weak stratification, so one can adopt the Boussinesq approximation. In the absence of rotation and magnetic field, we have

∇ · u = 0,  (1)
\frac{Du}{Dt} = −\frac{1}{ρ_{0}} ∇P + \frac{ρ}{ρ_{0}} g + ν∇^{2}u,  (2)
\frac{DT}{Dt} = κ∇^{2}T,  (3)

where g = (0, 0, −g) is gravity. The essential point here is that the fluid is incompressible and the density is almost constant except for a weak response to changes in the temperature, which leads to a temperature-dependent buoyancy of the form

ρ = ρ_{0} + (\frac{∂ρ}{∂T} )_{P} (T − T_{0}).  (4)

This coefficient is negative and it is proportional to what is usually called α, i.e., −ρ_{0}^{−1}(∂ρ/∂T)_{P} = α. For a perfect gas, α = 1/T, but one also wants to consider water and other liquids, where αT is of the order of 10^{6}; see Chandrasekhar (1961) for details.

The temperature is assumed fixed on the two boundaries, so the hydrostatic solution must obey ∇^{2}T = 0, which leads to a linear temperature profile, T_{0}(z) = T_{00} − βz, where β is the negative temperature gradient. The set of linearized equations can then be written in the form

∇ · u_{1} = 0  (5)
\frac{∂u_{1}}{∂t} = −\frac{1}{ρ_{0}} ∇P_{1} + αT_{1}gẑ + ν∇^{2}u_{1}  (6)
\frac{∂T_{1}}{∂t} − βu_{z1} = κ∇^{2}T_{1}.  (7)

2 Eliminating the pressure

A common trick is to apply the double curl. With the first curl, we eliminate the pressure, and with the second curl, the buoyancy term has just a z component, i.e.,

∇ × ∇ × ( \begin{matrix}0 \\ 0 \\ αgT_{1} \end{matrix} ) = ∇ × ( \begin{matrix}∂_{x} \\ ∂_{y} \\ ∂_{z} \end{matrix} ) × ( \begin{matrix}0 \\ 0 \\ αgT_{1} \end{matrix} ) = ( \begin{matrix}∂_{x} \\ ∂_{y} \\ ∂_{z} \end{matrix} ) × ( \begin{matrix}+αgT_{1,y} \\ −αgT_{1,x} \\ 0 \end{matrix} ) = αg ( \begin{matrix}T_{1,xz} \\ T_{1,yz} \\ −∇_{⊥}^{2}T_{1} \end{matrix} ) , (8)

where ∇_{⊥}^{2} = ∂_{x}^{2} + ∂_{y}^{2} is the horizontal Laplacian and commas denote partial derivatives. Furthermore, making use of the relation ∇ × ∇× = −∇^{2} + ∇∇· and the fact that ∇ · u_{1} = 0, we have ∇ × ∇ × u_{1} = −∇^{2}u_{1}, and so we arrive directly at

(\frac{∂}{∂t} − ν∇^{2}) ∇^{2}u_{1} = αgẑ∇_{⊥}^{2}T_{1},  (9)

or, more specifically for the z component of the velocity, at

(\frac{∂}{∂t} − ν∇^{2}) ∇^{2}u_{z1} = αg∇_{⊥}^{2}T_{1},  (10)

and for the temperature we have analogously

(\frac{∂}{∂t} − κ∇^{2}) T_{1} = βu_{z1}.  (11)

This is a system of two partial differential equations that can easily be combined into one:

(\frac{∂}{∂t} − ν∇^{2}) (\frac{∂}{∂t} − κ∇^{2}) ∇^{2}u_{z1} = αβg∇_{⊥}^{2}u_{z1}.  (12)

(This equation applies equally to T_{1}.) Furthermore, it is advantageous to nondimensionalize the equations by measuring length in units of d, which is the vertical thickness of the layer, so we write [x] = d, and time in units of the viscous time, [t] = d^{2}/ν.

(\frac{∂}{∂t} − ∇^{2}) (\frac{ν}{κ} \frac{∂}{∂t} − ∇^{2}) ∇^{2}u_{z1} = \frac{αβgd^{4}}{νκ} ∇_{⊥}^{2}u_{z1}.  (13)

As a consequence of nondimensionalization, two nondimensional numbers have appeared, the Prandtl and Rayleigh numbers,

Pr = \frac{ν}{κ}, Ra = \frac{αβgd^{4}}{νκ} .  (14)

respectively. Thus, we have

(\frac{∂}{∂t} − ∇^{2}) (Pr \frac{∂}{∂t} − ∇^{2}) ∇^{2}u_{z1} = Ra ∇_{⊥}^{2}u_{z1}.  (15)

3 Spectral ansatz for a specific boundary condition

Before considering the case of general boundary conditions, let us first consider a case where we can assume the solution to be of the form u_{z1} = û_{z1}(z) e^{σt+ik·x},

(σ + k^{2}) (Pr σ + k^{2}) k^{2} = Ra k_{⊥}^{2}.  (16)

thus, we have

Pr σ^{2} + (1 + Pr) σk^{2} + k^{4} = Ra \frac{k_{⊥}^{2}}{k^{2}} .  (17)

At this point, we can already say that there exists a nontrivial solution with σ = 0 when Ra = k^{6}/k_{⊥}^{2}. Such a solution corresponds to the marginally excited, nonoscillatory state.^{1}

If σ = 0 is the critical condition of interest, then the wavevector of interest is that corresponding to the longest wavelength that fits into the domain. Thus now depends on the boundary conditions. In fact, the only boundary condition that would be allowed with the adopted spectral ansatz is the stress-free condition, with zero temperature fluctuation, i.e.,

u_{1x,z} = u_{1y,z} = u_{1z} = T_{1} = 0 (on z = ±d/2).  (18)

In that case, we have

u_{1x} ∝ sin k_{z}z, u_{1y} ∝ sin k_{z}z, u_{1z} ∝ cos k_{z}z, T_{z} ∝ cos k_{z}z with k_{z} = π/d.  (19)

Figure 2: Ra(k_{⊥}^{2}) = (k_{⊥}^{2} + k_{z}^{2})^{3}/k_{⊥}^{2}. Now that k_{z} is fixed, we would still be free to vary k_{⊥}. The preferred value of k_{⊥} is the one that minimizes Ra(k_{⊥}^{2}) = (k_{⊥}^{2} + k_{z}^{2})^{3}/k_{⊥}^{2}. Thus we solve dRa/dk_{⊥}^{2} = 0, i.e.,

dRa/dk_{⊥}^{2} = 3(k_{⊥}^{2} + k_{z}^{2})^{2}/k_{⊥}^{2} − (k_{⊥}^{2} + k_{z}^{2})^{3}/k_{⊥}^{4} = 0. (20)

Thus, 3 k_{⊥}^{2} − (k_{⊥}^{2} + k_{z}^{2}) = 0, and therefore

k_{⊥}^{2} = k_{z}^{2}/2 = π^{2}/2,  (21)

1There could in principle also be another condition for a marginally excited solution where Re σ = 0, but Im σ ≠ 0, which would be oscillatory and thus with nonvanishing σ.

Figure 3: Dispersion relation for different values of Ra and Pr; increasing values are indicated by increasing thickness.

or k_{⊥} = π/\sqrt{2} ≈ 2.22; see Figure 2.

For k_{⊥} = π/\sqrt{2} and k_{z} = π, the value of Ra is

Ra = (k_{⊥}^{2} + k_{z}^{2})^{3}/k_{⊥}^{2} = (π^{2}/2 + π^{2})^{3}/(π^{2}/2) = \frac{27}{4}π^{4} ≈ 657.5 .  (22)

which is the critical value for the onset of convection in a horizontally infinitely extended layer. Obviously, if one does a simulation in a horizontally periodic domain, it won’t be infinitely extended, and therefore k_{⊥} can only be a multiple of the smallest horizontal domain size.

4 Dispersion relation

Note that have not yet solved the dispersion relation. Let us now solve Equation (17) for fixed k_{z} = π as a function of k_{⊥}^{2}. This, we have k^{2} = k_{⊥}^{2} + π^{2}. Let us rewrite Equation (17) in standard form as

σ^{2} + \frac{1 + Pr}{Pr} σk^{2} + \frac{k^{4}}{Pr} − \frac{Ra}{Pr} \frac{k_{⊥}^{2}}{k^{2}} = 0.  (23)

and write the solutions as

σ_{±} = −\frac{1 + Pr}{2 Pr}k^{2} ± \sqrt{\frac{(1 + Pr)^{2}}{4 Pr^{2}} k^{4} − \frac{k^{4}}{Pr} + \frac{Ra}{Pr} \frac{k_{⊥}^{2}}{k^{2}}} .  (24)

The result is plotted in Figure 3.

References

Brown, S. C., “Count Rumford Discovers Thermal Convection,” Daedalus 86, 340-343 (1957).

Chandrasekhar, S. Hydrodynamic and Hydromagnetic Stability. Dover Publications, New York (1961).

Schwarzschild, K., “On the equilibrium of the Sun’s atmosphere,” Nachr. Königl. Ges. Wissensch. Göttingen. Math.-phys. Klasse 195, 41-53 (1906).

$Header: /var/cvs/brandenb/tex/teach/ASTR_5410/3_Rayleigh_Benard/notes.tex,v 1.12 2016/09/12 21:32:52 brandenb Exp $