Onset of rotating Rayleigh–Bénard convection on the f-plane (Chandrasekhar 1961, Ch. III)¶
Family: convection · Script: examples/convection/rotating_rayleigh_benard.py
Obtained |
one pencil for all Ta (√Ta declared); free walls: Ra_c(Ta) = closed form to ≤ 2e-15 (Ta = 0…1e5); rigid: 1756.347 / 2151.341 / 4712.042 (Ta = 1e2, 1e3, 1e4) |
Reference |
Chandrasekhar 1961, ch. III |
A Boussinesq layer 0 < z < 1 heated from below rotates about the vertical with Taylor number
Ta = (2 Omega d^2 / nu)^2. Units: depth d, thermal diffusion time d^2 / kappa. A normal mode
exp(i k x + sigma t) independent of y (the horizontal plane is isotropic, so this loses no
generality) has all three velocity components; with D = d/dz and Pr = nu / kappa:
sigma u = Pr (D^2 - k^2) u - i k p + Pr sqrt(Ta) v (Coriolis: -2 Omega x u)
sigma v = Pr (D^2 - k^2) v - Pr sqrt(Ta) u
sigma w = Pr (D^2 - k^2) w - D p + Pr Ra T
sigma T = (D^2 - k^2) T + w (linear conduction profile)
0 = i k u + D w
Free (stress-free) walls: w = Du = Dv = T = 0. Onset is stationary for Pr > 0.6766 and then
(Chandrasekhar 1961, §26, closed form)
Ra(k) = [(pi^2 + k^2)^3 + pi^2 Ta] / k^2, Ra_c = min_k Ra(k).
Rigid walls (u = v = w = T = 0) have no closed form (Chandrasekhar's Table, Ch. III); we report
our values with their convergence in N.
Run: python examples/convection/rotating_rayleigh_benard.py
Run it: python examples/convection/rotating_rayleigh_benard.py (in the conatus environment).
1"""Onset of rotating Rayleigh–Bénard convection on the f-plane (Chandrasekhar 1961, Ch. III).
2
3A Boussinesq layer 0 < z < 1 heated from below rotates about the vertical with Taylor number
4Ta = (2 Omega d^2 / nu)^2. Units: depth d, thermal diffusion time d^2 / kappa. A normal mode
5exp(i k x + sigma t) independent of y (the horizontal plane is isotropic, so this loses no
6generality) has all three velocity components; with D = d/dz and Pr = nu / kappa:
7
8 sigma u = Pr (D^2 - k^2) u - i k p + Pr sqrt(Ta) v (Coriolis: -2 Omega x u)
9 sigma v = Pr (D^2 - k^2) v - Pr sqrt(Ta) u
10 sigma w = Pr (D^2 - k^2) w - D p + Pr Ra T
11 sigma T = (D^2 - k^2) T + w (linear conduction profile)
12 0 = i k u + D w
13
14Free (stress-free) walls: w = Du = Dv = T = 0. Onset is stationary for Pr > 0.6766 and then
15(Chandrasekhar 1961, §26, closed form)
16
17 Ra(k) = [(pi^2 + k^2)^3 + pi^2 Ta] / k^2, Ra_c = min_k Ra(k).
18
19Rigid walls (u = v = w = T = 0) have no closed form (Chandrasekhar's Table, Ch. III); we report
20our values with their convergence in N.
21
22Run: python examples/convection/rotating_rayleigh_benard.py
23"""
24import numpy as np
25import scipy.optimize as so
26
27from conatus.expert import EVP, Parameter, ParametricEVP, channel, critical, diff, interp, sqrt
28
29
30def rotating_layer(N: int, walls: str, Pr: float = 1.0):
31 """Linear stability pencil in (Ra, k, Ta) at fixed Pr; walls = 'free' or 'rigid'. The Coriolis
32 parameter sqrt(Ta) is not a Laurent polynomial in Ta: it is declared as a derived pencil
33 coordinate, so one pencil serves every Ta (it was rebuilt per Ta before)."""
34 ch = channel(N, base=None, interval=(0.0, 1.0)) # Chebyshev in z on [0, 1]
35 z = ch.y
36 u, v, w, p, T = (ch.field(n) for n in ("u", "v", "w", "p", "T"))
37 sigma = Parameter("sigma", dtype=complex)
38 Ra, k, Ta = Parameter("Ra", 1e3), Parameter("k", 3.0), Parameter("Ta", 1e3)
39 f = sqrt(Ta) # Coriolis parameter in these units
40 Dz = lambda q: diff(q, z) # noqa: E731
41 lap = lambda q: Dz(Dz(q)) - k ** 2 * q # noqa: E731
42 pb = EVP([u, v, w, p, T], eigenvalue=sigma, name=f"rotating RB ({walls} walls)")
43 pb.add_equation((sigma * u - Pr * lap(u) + 1j * k * p - Pr * f * v, 0))
44 pb.add_equation((sigma * v - Pr * lap(v) + Pr * f * u, 0))
45 pb.add_equation((sigma * w - Pr * lap(w) + Dz(p) - Pr * Ra * T, 0))
46 pb.add_equation((sigma * T - lap(T) - w, 0))
47 pb.add_equation((1j * k * u + Dz(w), 0))
48 # two conditions per second-order equation; add_bc lifts their taus into modes N-1, N-2
49 wall = (lambda q: Dz(q)) if walls == "free" else (lambda q: q) # stress-free or no-slip
50 for i, cond in enumerate((wall(u), wall(v), w, T)):
51 pb.add_bc(interp(cond, z, 0.0), interp(cond, z, 1.0), into=i)
52 return ParametricEVP(pb, [Ra, k, Ta], derived=[f]), pb
53
54
55def ra_free_exact(Ta: float) -> tuple[float, float]:
56 """Chandrasekhar's closed form for free walls: (Ra_c, k_c). dRa/dk = 0 gives, with x = k^2,
57 2 x^3 + 3 pi^2 x^2 - pi^6 - pi^2 Ta = 0 (one positive root)."""
58 x = so.brentq(lambda x: 2 * x ** 3 + 3 * np.pi ** 2 * x ** 2 - np.pi ** 6 - np.pi ** 2 * Ta,
59 1e-6, 1e4, xtol=1e-15)
60 return ((np.pi ** 2 + x) ** 3 + np.pi ** 2 * Ta) / x, float(np.sqrt(x))
61
62
63_PENCILS: dict = {}
64
65
66def onset(N: int, Ta: float, walls: str, k_bounds: tuple[float, float],
67 ra_bracket: tuple[float, float]):
68 """Critical Ra minimized over k: neutral curve by tracked Newton on Re(sigma), Brent in k.
69 One pencil per (N, walls), reused for every Ta."""
70 if (N, walls) not in _PENCILS:
71 _PENCILS[N, walls] = rotating_layer(N, walls)[0]
72 P = _PENCILS[N, walls]
73 # the stationary onset mode is the least stable one; critical() re-checks it by a full QZ
74 return critical(P, "Ra", bracket=ra_bracket, growth="real", over={"k": k_bounds},
75 values={"Ta": Ta}, xtol=1e-11), P
76
77
78if __name__ == "__main__":
79 print("Free walls: Ra_c(Ta) against Chandrasekhar's closed form (Pr = 1, stationary onset)")
80 for Ta in (0.0, 1e2, 1e3, 1e4, 1e5):
81 ex, kx = ra_free_exact(Ta)
82 for N in (16, 24):
83 cp, P = onset(N, Ta, "free", (0.8 * kx, 1.25 * kx), (0.8 * ex, 1.25 * ex))
84 print(f" Ta={Ta:8.0f} N={N}: Ra_c = {cp.value:.10f} (exact {ex:.10f}, rel. {cp.value / ex - 1:+.1e}), "
85 f"k_c = {cp.over['k']:.6f} (exact {kx:.6f}), sigma = {cp.mode.sigma:.1e}, residual {cp.mode.residual:.1e}")
86 assert abs(cp.value / ex - 1) < 1e-10 and abs(cp.over["k"] - kx) < 1e-5, (Ta, N, cp.value, ex)
87 print(cp.explain())
88 print("Rigid walls (no closed form; to compare with Chandrasekhar's table), convergence in N:")
89 for Ta in (0.0, 1e2, 1e3, 1e4, 1e5):
90 ex, kx = ra_free_exact(Ta) # free-wall values only bracket the search
91 kb = max(kx, 3.1)
92 for N in (24, 32, 48):
93 cp, _ = onset(N, Ta, "rigid", (0.7 * kb, 1.3 * kb), (0.5 * ex, 3.0 * ex))
94 print(f" Ta={Ta:8.0f} N={N}: Ra_c = {cp.value:.8f}, k_c = {cp.over['k']:.6f}, residual {cp.mode.residual:.1e}")
95 if Ta == 0.0: # Reid & Harris (1958), rigid walls
96 assert abs(cp.value - 1707.762) < 1e-3, cp.value
97 if Ta == 1e3: # Chandrasekhar (1961), Table VIII: 2151.3
98 assert abs(cp.value - 2151.34119768) < 1e-6, cp.value