Kelvin–Helmholtz instability of the tanh mixing layer (Michalke 1964), inviscid and viscous

Family: shear · Script: examples/shear/kelvin_helmholtz.py

Obtained

max ω_i = 0.189702 at α = 0.44492 (L = 20, N = 384); viscous → inviscid in O(1/Re)

Reference

Michalke 1964, JFM 19, Table 1: 0.1897 at α = 0.4446

 
Base flow U(y) = tanh(y) (lengths in units of the half-thickness, velocities in units of the
half velocity difference). A normal mode v(y) exp(i alpha (x - c t)) of the cross-stream velocity
obeys the Rayleigh equation (inviscid) or the Orr–Sommerfeld equation (viscous):

    (U - c)(v'' - alpha^2 v) - U'' v = (1 / (i alpha Re)) (v'''' - 2 alpha^2 v'' + alpha^4 v)

with the temporal growth rate omega_i = alpha Im(c). Michalke (1964, J. Fluid Mech. 19, 543-556,
Table 1 / Fig. 3) finds the maximum omega_i = 0.1897 at alpha = 0.4446 (c_r = 0 by symmetry).

The unbounded layer is truncated to y in [-L, L] with v = 0 (and v' = 0 when viscous) at y = +-L.
The eigenfunction decays like exp(-alpha |y|), so the truncation error is ~exp(-2 alpha L); both
L and N are explicit and their convergence is printed.

Run:  python examples/shear/kelvin_helmholtz.py

Run it: python examples/shear/kelvin_helmholtz.py (in the conatus environment).

 1"""Kelvin–Helmholtz instability of the tanh mixing layer (Michalke 1964), inviscid and viscous.
 2
 3Base flow U(y) = tanh(y) (lengths in units of the half-thickness, velocities in units of the
 4half velocity difference). A normal mode v(y) exp(i alpha (x - c t)) of the cross-stream velocity
 5obeys the Rayleigh equation (inviscid) or the Orr–Sommerfeld equation (viscous):
 6
 7    (U - c)(v'' - alpha^2 v) - U'' v = (1 / (i alpha Re)) (v'''' - 2 alpha^2 v'' + alpha^4 v)
 8
 9with the temporal growth rate omega_i = alpha Im(c). Michalke (1964, J. Fluid Mech. 19, 543-556,
10Table 1 / Fig. 3) finds the maximum omega_i = 0.1897 at alpha = 0.4446 (c_r = 0 by symmetry).
11
12The unbounded layer is truncated to y in [-L, L] with v = 0 (and v' = 0 when viscous) at y = +-L.
13The eigenfunction decays like exp(-alpha |y|), so the truncation error is ~exp(-2 alpha L); both
14L and N are explicit and their convergence is printed.
15
16Run:  python examples/shear/kelvin_helmholtz.py
17"""
18import numpy as np
19import scipy.optimize as so
20
21from conatus.expert import EVP, Parameter, ParametricEVP, channel, diff, interp, verify
22
23MICHALKE = {"alpha": 0.4446, "omega_i": 0.1897}      # Michalke (1964), Table 1
24
25
26def mixing_layer(N: int, L: float, viscous: bool):
27    """Rayleigh (viscous=False) or Orr–Sommerfeld (viscous=True) pencil in alpha (and Re)."""
28    ch = channel(N, base=None, interval=(-L, L))                # Chebyshev in y on [-L, L]
29    y, v, U = ch.y, ch.field("v"), ch.field("U")
30    c, alpha, Re = Parameter("c", dtype=complex), Parameter("alpha", 0.44), Parameter("Re", 1e3)
31
32    def L2(f):                                                  # d^2/dy^2 - alpha^2
33        return diff(diff(f, y), y) - alpha ** 2 * f
34
35    Upp = diff(diff(U, y), y)                                   # U'' from the NCC, not by hand
36    # (U - c) L2 v - U'' v [- L2 L2 v / (i alpha Re)] = 0; the eigenvalue c multiplies L2 v
37    rayleigh = U * L2(v) - c * L2(v) - Upp * v
38    pb = EVP([v], eigenvalue=c, name="Orr-Sommerfeld (tanh)" if viscous else "Rayleigh (tanh)")
39    if viscous:
40        pb.add_equation((rayleigh - L2(L2(v)) / (1j * alpha * Re), 0))
41        pb.add_bc(interp(v, y, -L), interp(v, y, L), interp(diff(v, y), y, -L), interp(diff(v, y), y, L))
42        params = [alpha, Re]
43    else:
44        pb.add_equation((rayleigh, 0))
45        pb.add_bc(interp(v, y, -L), interp(v, y, L))
46        params = [alpha]
47    # U'' is computed from the sampled tanh: on [-15, 15] at N = 192 its Chebyshev tail is 1.1e-5
48    # (the interpolation error of tanh, 1e-10, amplified by d^2/dy^2; largest at |y| = L where
49    # U'' ~ 0 and the mode ~ exp(-alpha L)); accepted explicitly up to 2e-5. The gallery
50    # test checks the growth rate against Michalke, and __main__ its convergence in N and L.
51    return ParametricEVP(pb, params, ncc={U: np.tanh}, ncc_tail_tol=2e-5), pb
52
53
54def most_unstable(P: ParametricEVP, **fixed):
55    """max over alpha of omega_i = alpha Im(c): full QZ once, then the KH mode is tracked in alpha."""
56    # KH mode = largest Im(c) among |c| < 2 (|U| <= 1 bounds physical phase speeds, Howard's
57    # semicircle); the viscous tau-spurious modes have |c| ~ 1e3 and are excluded (count logged)
58    seed = P.least_stable(growth="imag", select=lambda s: abs(s) < 2, alpha=0.44, **fixed)
59
60    def minus_omega_i(a):
61        m = P.track(seed.sigma, seed.x, left=False, alpha=a, **fixed)
62        return -a * m.sigma.imag
63
64    opt = so.minimize_scalar(minus_omega_i, bounds=(0.40, 0.49), method="bounded",
65                             options={"xatol": 1e-7})
66    mode = P.track(seed.sigma, seed.x, alpha=float(opt.x), **fixed)
67    return float(opt.x), mode
68
69
70if __name__ == "__main__":
71    print("Inviscid (Rayleigh), convergence in N and L:")
72    for L, N in [(10, 192), (15, 192), (15, 256), (20, 384)]:
73        P, pb = mixing_layer(N, L, viscous=False)
74        a, m = most_unstable(P)
75        print(f"  L={L:2d} N={N:3d}: alpha_max = {a:.6f}, omega_i = {a * m.sigma.imag:.8f}, "
76              f"c = {m.sigma:.2e}, residual {m.residual:.1e}")
77    print(f"  Michalke (1964): alpha_max = {MICHALKE['alpha']}, omega_i = {MICHALKE['omega_i']}")
78    # Michalke's 4 printed digits, at the largest L and N (L = 20, N = 384)
79    assert abs(a * m.sigma.imag - MICHALKE["omega_i"]) < 5e-5 and abs(a - MICHALKE["alpha"]) < 5e-4, (a, m.sigma)
80    print(pb.explain())
81    print(verify(P, m).explain())
82    print("Viscous (Orr–Sommerfeld), L=15, N=256: omega_i -> inviscid value as Re -> infinity")
83    for Re in (1e2, 1e3, 1e4):
84        P, _ = mixing_layer(256, 15, viscous=True)
85        a, m = most_unstable(P, Re=Re)
86        print(f"  Re = {Re:7.0f}: alpha_max = {a:.6f}, omega_i = {a * m.sigma.imag:.8f}, residual {m.residual:.1e}")