鍋底から始まったカオス — レイリー・ベナール対流とローレンツ方程式

層の厚さを2倍にするとレイリー数は8倍です。その8倍はρ=8で、カオスの閾値24.74にはまだ届きません。

1900年、アンリ・ベナールは薄い油の層を底から温めて、液体が自ら六角形のハニカムに割れるのを見ました。かき混ぜてはいません。16年後にレイリー卿がその閾値を一つの数で書き、さらに47年後、ローレンツは同じ系からカオスを取り出しました。

閾値はレイリー数1708です。下は伝導、上はロールです。ロールを三つの式へ切ると加熱パラメータρが残り、24.74を超えると固定点は軌道を捕まえられません。

Raのスライダーを1708より下へ下げると粒子はほとんど止まり、その上へ押すと互いに逆回りのロールが点火します。

Ra = 1,000

Ra_c = 1708전도 (still) — conduction

1708より下では温度は水平の層として積もるだけです。閾値を越えてスライダーをさらに押すと、同じロールがより速く回ります。

静止した液体はなぜ動き出すのか

下が熱く、上が冷たい層です。温められた下側は膨張して軽くなります。軽いものが下、重いものが上。浮力はこの積み重ねをひっくり返そうとします。

すぐには噴き上がりません。粘性は流れを押さえ、熱拡散は上昇する塊の熱を周囲へ逃がして浮力を殺します。どちらもブレーキです。

浮力が勝てば循環、つまり対流ロールが生まれます。ブレーキが勝てば液体はそのままで、熱は伝導だけを通ります。浮力が経路の途中で使い切られると、ロールまで届きません。煙道を1 m下げてドラフト7 Paが消えた場合が、その予算です。

レイリー数 — ひっくり返す力と引き止める力の比

綱引きを一つの無次元数へ押さえるとレイリー数です。

\mathrm{Ra} = \frac{g \beta \, \Delta T \, L^3}{\nu \, \alpha}

$g$ は重力加速度、 $\beta$ は熱膨張係数(温度が上がると体積が増える度合い)、 $\Delta T$ は上下の温度差、 $L$ は液体層の厚さ、 $\nu$ は動粘性係数、 $\alpha$ は熱拡散係数です。分子はひっくり返そうとする浮力、分母は引き止める粘性と拡散です。

大きい項は $L^3$ です。層が厚くなるとレイリー数は三乗で増えます。薄い油膜はなかなか対流しません。深い海やマントルはごく小さな温度差でも動きます。鍋を少し深くするだけで閾値が一気に近づく理由です。

臨界レイリー数1708 — 対流が点火する閾値#

レイリーの線形安定性解析は、境界ごとに数が変わります。上下が自由な理想境界では閾値は $\mathrm{Ra}_c = 27\pi^4/4 \approx 657.5$ 。実験で一般的な剛体壁では $\mathrm{Ra}_c \approx 1708$ です。

$\mathrm{Ra} < \mathrm{Ra}_c$ なら、小さな揺らぎは粘性と拡散に食われて消えます。液体は静止へ戻ります。 $\mathrm{Ra} > \mathrm{Ra}_c$ を超えた瞬間、特定の波長の揺らぎが自ら育ちます。対流ロールの種です。

1708を657.5で割ると約2.6です。同じ流体、同じ厚さなら、剛体壁は温度差を約2.6倍要求します。厚さで払うなら三乗根なので、 $2.6^{1/3} \approx 1.37$ 倍で足ります。自由滑りの蓋を剛体へ替える代金がこれです。

対流ロールの形とヌッセルト数

閾値直後のロールは特定の波長を選びます。臨界波数は $k_c \approx 3.117/L$ です。一波長の長さは $2\pi/k$ で、 $3.117$ で割ると約 $2.02L$ です。隣り合うロール二つがその波長を埋め、ロール一つの幅は半分の $L$ 付近です。正方形だという言い方は、この割り算です。最も低い Ra で点火するモードが、その波長を持ちます。

熱がどれだけ余分に出るかは**ヌッセルト数(全熱伝達と純伝導熱伝達の比)**です。

\mathrm{Nu} = \frac{q_\text{total}}{q_\text{conduction}}

$q_\text{total}$ は実際の熱流束、 $q_\text{conduction}$ は対流がないときの伝導熱流束です。伝導だけなら $\mathrm{Nu}=1$ 。対流が強まると実験では $\mathrm{Nu} \sim \mathrm{Ra}^{1/3}$ 程度です。この指数だけを読むと、Ra が8倍のとき Nu は2倍です。閾値の直後には当てはまりません。そこでは Nu がまだ1の近くです。

ローレンツが三つの方程式へ縮めた対流

1963年、エドワード・ローレンツは温度場と流れ場の最も粗いモードだけを残しました。ガラーキン切断です。偏微分方程式が常微分方程式三つへ縮みます。

\begin{aligned} \dot{x} &= \sigma (y - x) \\ \dot{y} &= x(\rho - z) - y \\ \dot{z} &= xy - \beta z \end{aligned}

$x$ は対流ロールの循環の強さ、 $y$ は上昇・下降する柱の温度差、 $z$ は鉛直温度分布の歪みです。 $\sigma$ はプラントル数(運動量拡散と熱拡散の比)、 $\beta$ は幾何定数、 $\rho$ は $\mathrm{Ra}/\mathrm{Ra}_c$ です。

$\rho$ を上げることは、鍋をより強く温めることです。 $\rho<1$ なら軌道は原点へ向かいます。伝導です。 $1<\rho<24.74$ なら原点ではない固定点へ集まります。定常ロールです。 $\rho>24.74$ を超えると、どの固定点にも着地せず二つの翼の間を回ります。カオスです。

定常点は三つの導関数を0に置いた場所です。 $\sigma \neq 0$ なら $y=x$ です。 $x \neq 0$ なら $z=\rho-1$ です。残る式は $x^2 = \beta(\rho-1)$ です。コードの $\beta = 8/3$ 、 $\rho = 15$ なら $z=14$ 、 $x^2 = 112/3 \approx 37.3$ です。真ん中の出力 steady roll は、尾の4000点で $x$ の幅が0.1より狭いという意味です。 $\rho=28$ なら $z=27$ 、 $x^2=72$ で二つの翼がその位置にあり、出力 chaos はその尾が0.1の門を通れなかったという意味です。

Python — 対流の三つの顔を積分する#

4次ルンゲ・クッタです。 $\sigma=10$ 、 $\beta=8/3$ 、時間刻み0.01、20000歩。尾は最後の4000点です。

import numpy as np
 
SIGMA, BETA = 10.0, 8.0 / 3.0
 
def lorenz_deriv(state, rho):
    x, y, z = state
    return np.array([
        SIGMA * (y - x),        # ロール循環の強さの変化
        x * (rho - z) - y,      # 上昇・下降の温度差の変化
        x * y - BETA * z,       # 鉛直温度歪みの変化
    ])
 
def integrate_rk4(rho, dt=0.01, n=20000, s0=(0.6, 0.4, 0.2)):
    s = np.array(s0, dtype=float)
    traj = np.empty((n, 3))
    for i in range(n):
        k1 = lorenz_deriv(s, rho)
        k2 = lorenz_deriv(s + 0.5 * dt * k1, rho)
        k3 = lorenz_deriv(s + 0.5 * dt * k2, rho)
        k4 = lorenz_deriv(s + dt * k3, rho)
        s = s + (dt / 6.0) * (k1 + 2 * k2 + 2 * k3 + k4)
        traj[i] = s
    return traj
 
def convection_regime(rho):
    tail = integrate_rk4(rho)[-4000:]          # 過渡を捨てた尾の部分
    spread = tail[:, 0].max() - tail[:, 0].min()
    if abs(tail[-1, 0]) < 1e-3 and rho < 1.0:
        return "conduction"                     # 原点へ収束
    if spread < 0.1:
        return "steady roll"                    # 固定点へ収束
    return "chaos"                              # 二つの翼の間で振動
 
for rho in (0.8, 15.0, 28.0):
    print(f"rho={rho:5.1f} -> {convection_regime(rho)}")

出力は3行です。

rho=  0.8 -> conduction
rho= 15.0 -> steady roll
rho= 28.0 -> chaos

ρ (Ra/Ra_c) = 28.0

카오스 — 나비 날개 (chaotic convection)

가로축 x = 롤 순환 세기, 세로축 z = 수직 온도 왜곡. 노란 점 = 고정점 C±.

ρのスライダーを24.74より下に置くと、軌道は黄色い固定点の一つへ入ります。超えると蝶の翼が開きます。reset は初期点を少しずらします。最初は重なっていた二つの軌道が、やがて分かれます。

1708を超えただけで終わったと書くと失敗する#

厚さを2倍にすると $L^3$ のため Ra は8倍です。ちょうど1708に乗っていた層をその深さにすると $\rho = \mathrm{Ra}/\mathrm{Ra}_c = 8$ です。ローレンツの定常ロール区間は1から24.74の間です。8はその中にあります。カオスまで行くにはρを24.74まで上げる必要があり、同じ温度差なら厚さは $24.74^{1/3} \approx 2.9$ 倍です。閾値直後のロールをカオスと書くと、この間隔の分だけ外れます。

その間隔は、 $\beta$ を定数のまま $\Delta T$ だけを上げて作っています。温度差が30 Kを超えた瞬間、この定数のほうが先に嘘になりはしないでしょうか。その温度でブシネスクが崩れる側を開くと、1708と24.74のどちらが先に数字を失うかが見えます。