Rayleigh–Bénard onset in 3D: the nonlinear Boussinesq equations, linearized by the IR

Family: convection · Script: examples/convection/rayleigh_benard_3d.py

Obtained

Ra = 1707.76178 at

Reference

k

 
Layer -1/2 < z < 1/2 between rigid isothermal walls, periodic in x and y (Fourier x Fourier x
Chebyshev), units of depth and thermal diffusion time:

    div u = 0
    dt T - lap T = -u . grad T
    dt u - Pr lap u + grad p - Pr Ra T e_z = -u . grad u

around the conduction state u = 0, T0 = -z. linearize() gives the eigenvalue problem; each
Fourier group (kx, ky) is one pencil. Onset: Ra_c = 1707.762 at |k| = 3.1163 (Reid & Harris 1958;
Chandrasekhar 1961, Ch. II, Table III). Isotropy: the growth rate depends on |k| only, so
(kx, ky) = (K, 0), (0, K) and (K cos 30°, K sin 30°) are all marginal at Ra_c.

Run:  python examples/convection/rayleigh_benard_3d.py

Run it: python examples/convection/rayleigh_benard_3d.py (in the conatus environment).

 1"""Rayleigh–Bénard onset in 3D: the nonlinear Boussinesq equations, linearized by the IR.
 2
 3Layer -1/2 < z < 1/2 between rigid isothermal walls, periodic in x and y (Fourier x Fourier x
 4Chebyshev), units of depth and thermal diffusion time:
 5
 6    div u = 0
 7    dt T - lap T = -u . grad T
 8    dt u - Pr lap u + grad p - Pr Ra T e_z = -u . grad u
 9
10around the conduction state u = 0, T0 = -z. linearize() gives the eigenvalue problem; each
11Fourier group (kx, ky) is one pencil. Onset: Ra_c = 1707.762 at |k| = 3.1163 (Reid & Harris 1958;
12Chandrasekhar 1961, Ch. II, Table III). Isotropy: the growth rate depends on |k| only, so
13(kx, ky) = (K, 0), (0, K) and (K cos 30°, K sin 30°) are all marginal at Ra_c.
14
15Run:  python examples/convection/rayleigh_benard_3d.py
16"""
17import numpy as np
18import scipy.optimize as so
19
20from conatus.expert import (IVP, CartesianCoordinates, Chebyshev, ComplexFourier, Domain, Field,
21                            Parameter, VectorField, div, dot, dt, eigenvalues, grad, integ, interp,
22                            lap, lift, discretize, linearize, unit_vector)
23from conatus.ir import Zero
24
25RA_C, K_C = 1707.762, 3.1163           # rigid-rigid onset (Reid & Harris 1958)
26
27
28def boussinesq_3d(kx: float, ky: float, Nz: int):
29    """Nonlinear Boussinesq IVP on a box holding exactly the wavevector (kx, ky) (4 modes a side)."""
30    cs = CartesianCoordinates("x", "y", "z")
31    x, y, z = cs.coords
32    # box lengths 2 pi / k hold the wavevector exactly (a zero component: any length, here 2 pi)
33    bx = ComplexFourier(x, 4, (0, 2 * np.pi / (kx or 1.0)), dealias=1)
34    by = ComplexFourier(y, 4, (0, 2 * np.pi / (ky or 1.0)), dealias=1)
35    bz = Chebyshev(z, Nz, (-0.5, 0.5), dealias=1)
36    D = Domain(cs, [bx, by, bz])
37    Dh = D.without(z)
38    u, T, p = VectorField("u", D), Field("T", D), Field("p", D)
39    tp = Field("tp", Domain(cs, []))
40    tT1, tT2 = Field("tT1", Dh), Field("tT2", Dh)
41    tu1, tu2 = VectorField("tu1", Dh), VectorField("tu2", Dh)
42    Ra, Pr = Parameter("Ra", RA_C), Parameter("Pr", 1.0)
43    ez = unit_vector(z)
44    ivp = IVP([u, T, p, tp, tT1, tT2, tu1, tu2], name="Boussinesq RB 3D")
45    ivp.add_equation((div(u) + tp, 0))
46    ivp.add_equation((dt(T) - lap(T) + lift(tT1, bz, -1) + lift(tT2, bz, -2), -dot(u, grad(T))))
47    ivp.add_equation((dt(u) - Pr * lap(u) + grad(p) - Pr * Ra * T * ez
48                      + lift(tu1, bz, -1) + lift(tu2, bz, -2), -dot(u, grad(u))))
49    ivp.add_equation((interp(T, z, "left"), 0.5))
50    ivp.add_equation((interp(T, z, "right"), -0.5))
51    ivp.add_equation((interp(u, z, "left"), 0))
52    ivp.add_equation((interp(u, z, "right"), 0))
53    ivp.add_equation((integ(p), 0))
54    T0 = Field("T0", Domain(cs, [bz]))
55    evp = linearize(ivp, around={u: Zero((cs,)), T: T0}, eigenvalue=Parameter("s", dtype=complex))
56    return evp, T0
57
58
59def growth(kx: float, ky: float, Ra: float, Nz: int = 24):
60    """Leading growth rate of the (kx, ky) group at Ra, with its residual."""
61    evp, T0 = boussinesq_3d(kx, ky, Nz)
62    plan = discretize(evp, ncc={T0: lambda X, Y, Z: -Z}, params={"Ra": Ra})
63    g = next(g for g in plan.groups if np.isclose(g[0].k, kx) and np.isclose(g[1].k, ky))
64    r = eigenvalues(plan, g)
65    i = int(np.argmax(r.values.real))
66    return r.values[i], r.residuals[i]
67
68
69def neutral_ra(kx: float, ky: float, Nz: int = 24) -> float:
70    """Ra where the leading (real) growth rate of the (kx, ky) group vanishes."""
71    return so.brentq(lambda Ra: growth(kx, ky, Ra, Nz)[0].real, 1600.0, 1800.0, xtol=1e-9)
72
73
74if __name__ == "__main__":
75    K = K_C
76    ras = []
77    for kx, ky in [(K, 0.0), (0.0, K), (K * np.cos(np.pi / 6), K * np.sin(np.pi / 6))]:
78        Ra = neutral_ra(kx, ky)
79        s, res = growth(kx, ky, Ra)
80        print(f"  k = ({kx:.4f}, {ky:.4f}), |k| = {np.hypot(kx, ky):.4f}: neutral Ra = {Ra:.8f} "
81              f"(2D reference {RA_C} at the minimum), sigma = {s:.1e}, residual {res:.1e}")
82        assert abs(s) < 1e-9 and res < 1e-14, (s, res)                # neutral, resolved
83        ras.append(Ra)
84    # Reid & Harris (1958) to 1e-3; isotropic: the same Ra for every direction of k to 1e-8
85    assert abs(ras[0] - RA_C) < 1e-3 and max(ras) - min(ras) < 1e-8 * ras[0], ras
86    evp, _ = boussinesq_3d(K, 0.0, 24)
87    print(evp.explain())