Salt fingers: stationary double-diffusive onset in a layer with free walls (Stern 1960; Baines & Gill 1969, J. Fluid Mech. 37, 289-306)¶
Family: convection · Script: examples/convection/salt_fingers.py
Obtained |
R_S/τ − R_T = 27π⁴/4 to ≤ 4e-16 (4 pairs) |
Reference |
Stern 1960; Baines & Gill 1969, JFM 37 |
Warm salty water above cold fresh water: temperature is stabilizing (Rayleigh number R_T > 0)
and salt, which diffuses tau = kappa_S / kappa_T < 1 times slower, is destabilizing (R_S > 0).
Units: depth d, thermal diffusion time d^2 / kappa_T; normal mode exp(i k x + sigma t), D = d/dz:
sigma u = Pr (D^2 - k^2) u - i k p
sigma w = Pr (D^2 - k^2) w - D p + Pr (R_T T - R_S S) (heat and salt both make it lighter)
sigma T = (D^2 - k^2) T - w (stable temperature gradient)
sigma S = tau (D^2 - k^2) S - w (unstable salinity gradient)
0 = i k u + D w
Free walls: w = Du = T = S = 0. With sin(pi z) modes and sigma = 0, the stationary (finger)
threshold is R_S / tau - R_T = min_k (pi^2 + k^2)^3 / k^2 = 27 pi^4 / 4 at k^2 = pi^2 / 2,
independent of Pr (Baines & Gill 1969, §2).
Run: python examples/convection/salt_fingers.py
Run it: python examples/convection/salt_fingers.py (in the conatus environment).
1"""Salt fingers: stationary double-diffusive onset in a layer with free walls (Stern 1960;
2Baines & Gill 1969, J. Fluid Mech. 37, 289-306).
3
4Warm salty water above cold fresh water: temperature is stabilizing (Rayleigh number R_T > 0)
5and salt, which diffuses tau = kappa_S / kappa_T < 1 times slower, is destabilizing (R_S > 0).
6Units: depth d, thermal diffusion time d^2 / kappa_T; normal mode exp(i k x + sigma t), D = d/dz:
7
8 sigma u = Pr (D^2 - k^2) u - i k p
9 sigma w = Pr (D^2 - k^2) w - D p + Pr (R_T T - R_S S) (heat and salt both make it lighter)
10 sigma T = (D^2 - k^2) T - w (stable temperature gradient)
11 sigma S = tau (D^2 - k^2) S - w (unstable salinity gradient)
12 0 = i k u + D w
13
14Free walls: w = Du = T = S = 0. With sin(pi z) modes and sigma = 0, the stationary (finger)
15threshold is R_S / tau - R_T = min_k (pi^2 + k^2)^3 / k^2 = 27 pi^4 / 4 at k^2 = pi^2 / 2,
16independent of Pr (Baines & Gill 1969, §2).
17
18Run: python examples/convection/salt_fingers.py
19"""
20import numpy as np
21
22from conatus.expert import EVP, Parameter, ParametricEVP, channel, critical, diff, interp
23
24RA_FREE = 27 * np.pi ** 4 / 4 # 657.511..., the free-wall Rayleigh–Bénard number
25
26
27def finger_layer(N: int, R_T: float, tau: float, Pr: float = 7.0):
28 """Linear stability pencil in (R_S, k) at fixed R_T, tau, Pr (free walls)."""
29 ch = channel(N, base=None, interval=(0.0, 1.0)) # Chebyshev in z on [0, 1]
30 z = ch.y
31 u, w, p, T, S = (ch.field(n) for n in ("u", "w", "p", "T", "S"))
32 sigma = Parameter("sigma", dtype=complex)
33 R_S, k = Parameter("R_S", 1e3), Parameter("k", 2.2)
34 Dz = lambda q: diff(q, z) # noqa: E731
35 lap = lambda q: Dz(Dz(q)) - k ** 2 * q # noqa: E731
36 pb = EVP([u, w, p, T, S], eigenvalue=sigma, name=f"salt fingers (R_T={R_T:g}, tau={tau:g}, Pr={Pr:g})")
37 pb.add_equation((sigma * u - Pr * lap(u) + 1j * k * p, 0))
38 pb.add_equation((sigma * w - Pr * lap(w) + Dz(p) - Pr * (R_T * T - R_S * S), 0))
39 pb.add_equation((sigma * T - lap(T) + w, 0))
40 pb.add_equation((sigma * S - tau * lap(S) + w, 0))
41 pb.add_equation((1j * k * u + Dz(w), 0))
42 # free walls; add_bc lifts the two taus of each second-order equation into modes N-1, N-2
43 for i, cond in enumerate((Dz(u), w, T, S)):
44 pb.add_bc(interp(cond, z, 0.0), interp(cond, z, 1.0), into=i)
45 return ParametricEVP(pb, [R_S, k]), pb
46
47
48def finger_onset(N: int, R_T: float, tau: float, Pr: float = 7.0):
49 """Critical R_S (growth Re sigma = 0) minimized over k; exact: tau (R_T + 27 pi^4 / 4)."""
50 P, _ = finger_layer(N, R_T, tau, Pr)
51 ex = tau * (R_T + RA_FREE)
52 return critical(P, "R_S", bracket=(0.8 * ex, 1.25 * ex), growth="real",
53 over={"k": (1.8, 2.7)}, xtol=1e-11), ex
54
55
56if __name__ == "__main__":
57 print("Stationary finger threshold R_S,c = tau (R_T + 27 pi^4 / 4), k_c = pi / sqrt(2)")
58 for R_T, tau in [(0.0, 0.01), (1e3, 0.01), (1e4, 0.1), (1e5, 1 / 3)]:
59 for N in (16, 24):
60 cp, ex = finger_onset(N, R_T, tau)
61 print(f" R_T={R_T:8.0f} tau={tau:.3f} N={N}: R_S,c = {cp.value:.8f} (exact {ex:.8f}, "
62 f"rel. {cp.value / ex - 1:+.1e}), k_c = {cp.over['k']:.6f} "
63 f"(exact {np.pi / np.sqrt(2):.6f}), residual {cp.mode.residual:.1e}")
64 assert abs(cp.value / ex - 1) < 1e-10 and abs(cp.over["k"] - np.pi / np.sqrt(2)) < 1e-5, (R_T, tau, N)
65 print(cp.explain())