静止したまま温まる流体 — 自然対流とレイリー数1708

下から火をつけても、流体はしばらく動きません。床は熱く天井は冷たいのに、液体は静止したまま熱だけを上へ送ります。そしてある瞬間、縞模様の渦が一斉にスイッチオンします。この記事では、その「閾値」がなぜレイリー数1708なのかを追い、温められた壁がどれだけよく冷えるかを決めるヌセルト相関式まで進みます。最後には臨界値を自分で越えてみて、鉛直平板の熱損失をPythonで計算します。

浮力と粘性がぶつかる場所#

自然対流(外部ポンプなしに密度差が生む流れ)は三つの力の争いです。床を温めるとその上の流体が膨張して軽くなります。軽い塊は浮き上がろうとします。これが浮力です。

浮力だけなら流れは即座に始まるはずです。しかし二つのものが足を引っ張ります。粘性は浮き上がる塊を周囲との摩擦で引き下ろします。熱拡散は熱い塊の温度を素早くぼかし、浮力の原動力である温度差そのものを消します。

だから結論は単純な競争です。二つのブレーキ(粘性+熱拡散)が浮力に勝てば流体は静止します。浮力が勝てば対流がスイッチオンします。どちらが勝つかを一つの数で答えるのが本記事の狙いです。

ブシネスク近似 — 密度変化を一項に#

密度が温度で変わると方程式は煩雑になります。ブシネスク近似(密度を重力項の中だけ変数として扱う簡略化)はこの荷を軽くします。それ以外の場所では密度は定数、浮力項だけが温度に反応します。

体積熱膨張係数(温度が1度上がるときに密度が減る割合)を定義します。

ここで $\beta$ は膨張係数、$\rho$ は密度、$T$ は温度です。理想気体なら $\beta = 1/T$ にきれいに落ちます。

この定義で密度差を温度差に置き換えます。

$\rho_\infty$、$T_\infty$ は遠方の静止流体の密度と温度です。この一行が鉛直平板境界層の運動量方程式に浮力源を植え付けます。

$u$、$v$ は速度成分、$g$ は重力加速度、$\nu$ は動粘性係数です。右辺第一項が浮力、第二項が粘性抵抗です。エネルギー方程式は温度を拡散で運びます。

$\alpha$ は熱拡散率(熱が広がる速さ)です。浮力で結ばれたこの一対(運動量とエネルギー)が自然対流の骨格です。

グラスホフとレイリー — 二つの無次元数の分業#

無次元化すると浮力項は一つの数にまとまります。それがグラスホフ数です。

$T_s$ は壁温度、$L$ は特性長さです。グラスホフは強制対流でレイノルズ数が担う役を引き受けます。浮力が粘性をどれだけ圧倒するかを測ります。

しかし浮力の敵は粘性だけではありません。熱拡散も温度差を消します。この二つ目のブレーキを勘定に入れるにはプラントル数(運動量拡散と熱拡散の比、$\mathrm{Pr} = \nu/\alpha$)を掛けます。その結果がレイリー数です。

分母に粘性 $\nu$ と熱拡散 $\alpha$ が一緒に入ります。これが核心です。対流が点火するかは Gr 単独ではなく Ra が決めます。浮力が二つのブレーキの積に勝ってはじめて流れが生きます。

1708という閾値#

下から温められる流体層を線形安定性解析に入れると、きれいな臨界値が現れます。上下とも固体壁なら臨界レイリー数は $\mathrm{Ra}_c = 1708$ です。これより小さいと、あらゆる擾乱が粘性と熱拡散に食われて消えます。流体は静止し、熱は純粋な伝導だけで移動します。

Ra が1708を越えた瞬間、パターンがスイッチオンします。規則正しく並んで逆回転するロール、すなわちベナールセルです。セルの幅は層の厚さの約2倍です。

下のシミュレーションで直接操作してみましょう。

Ra を1708より下げると粒子はほぼ凍りつき、温度場は上下にきれいに分かれます。1708を越えると逆回転のセルがスイッチオンし、粒子が循環し始めます。Ra をさらに上げると循環が激しくなり、境界層が薄くなります。

Python — 鉛直平板のヌセルト数を追う#

温められた壁はどれだけよく冷えるか? 答えはヌセルト数(伝導に対する対流熱伝達の比、$\mathrm{Nu} = hL/k$)です。鉛直等温平板には全域で通用するチャーチル–チュー相関式があります。60度の壁を20度の空気中に立てた場合を計算しましょう。

import numpy as np
 
def rayleigh_number(Ts, Tinf, L, beta, nu, alpha, g=9.81):
    # 浮力と (粘性 x 熱拡散) の比
    return g * beta * (Ts - Tinf) * L**3 / (nu * alpha)
 
def churchill_chu_vertical(Ra, Pr):
    # 鉛直等温平板 -- 全域で有効
    denom = (1.0 + (0.492 / Pr)**(9 / 16))**(8 / 27)
    return (0.825 + 0.387 * Ra**(1 / 6) / denom)**2
 
# 膜温度基準の空気物性
nu, alpha, Pr, k = 1.85e-5, 2.60e-5, 0.71, 0.027
Ts, Tinf, L = 60.0, 20.0, 0.30          # 60C 壁, 0.3 m の高さ
beta = 1.0 / (0.5 * (Ts + Tinf) + 273.15)  # 理想気体 1/Tf
 
Ra = rayleigh_number(Ts, Tinf, L, beta, nu, alpha)
Nu = churchill_chu_vertical(Ra, Pr)
h = Nu * k / L
q = h * (Ts - Tinf)                     # 単位面積あたりの熱損失
 
print(f"Ra = {Ra:.2e}")
print(f"Nu = {Nu:.1f},  h = {h:.2f} W/m^2K")
print(f"q  = {q:.1f} W/m^2")
 
# Ra = 7.03e+07
# Nu = 54.9,  h = 4.94 W/m^2K
# q  = 197.6 W/m^2

$\mathrm{Ra} \approx 7\times10^7$ なのでまだ層流域です。平板1平方メートルが約198ワットを失います。温度差や高さを変えると Ra はすぐに動きます。

下のグラフは同じ相関式を Ra 軸上に広げたものです。

Pr を下げると(液体金属側)曲線は下がり、上げると(油側)上がります。Ra マーカーを $10^9$ の先へ引くと、傾きが層流の $\mathrm{Ra}^{1/4}$ から乱流の $\mathrm{Ra}^{1/3}$ に切り替わります。

ヌセルト相関式のひとまとめ#

形状ごとに相関式は異なります。よく使うものをまとめました。

• 鉛直平板: $\mathrm{Nu} = 0.59\,\mathrm{Ra}^{1/4}$ (層流, $10^4$–$10^9$)、$\mathrm{Nu} = 0.10\,\mathrm{Ra}^{1/3}$ (乱流, $10^9$–$10^{13}$)。全域はチャーチル–チュー。
• 水平加熱板、上向き: $\mathrm{Nu} = 0.54\,\mathrm{Ra}^{1/4}$ ($10^4$–$10^7$)、$0.15\,\mathrm{Ra}^{1/3}$ ($10^7$–$10^{11}$)。
• 水平加熱板、下向き: $\mathrm{Nu} = 0.27\,\mathrm{Ra}^{1/4}$。
• 下から温める密閉空間: $\mathrm{Ra}_c = 1708$ を越えるとベナールセル。

水平板の特性長さは面積÷周長 $L = A_s/P$ で取ります。乱流域の指数 $1/3$ には隠れた意味があります。$\mathrm{Nu} \propto \mathrm{Ra}^{1/3}$ かつ $\mathrm{Ra} \propto L^3$ なので、$h = \mathrm{Nu}\,k/L$ で $L$ が相殺されます。乱流自然対流の熱伝達係数は大きさに依存しません。

覚えておくこと#

• 浮力の敵は二つ、粘性と熱拡散です。だから対流の開始は Gr 単独ではなく、二つを掛け合わせた Ra が決めます。
• $\mathrm{Ra} < 1708$ なら流体は静止(純粋伝導)、越えるとベナールセルがスイッチオンします。
• $\mathrm{Nu}$ は層流で $\mathrm{Ra}^{1/4}$、乱流で $\mathrm{Ra}^{1/3}$。後者では熱伝達係数が平板の大きさに依存しなくなります。