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

同じ Ra に 1/3 の指数を載せると熱流束は約 149 まで落ち、1e9 は高さ約 0.73 m で初めて届きます。

床を温めても、流体はしばらく止まったままです。動くのは熱だけです。縞の渦は一斉に点きます。その戸の数字がレイリー数 1708 です。

Ra のスライダーを 1708 より下げると粒子は止まり、温度は上下に分かれます。その上に上げると、互いに逆回転するロールが点きます。

Ra = 1.0e+4

critical Ra_c = 1708 · state: steady convection rolls (Benard cells)

閾値のすぐ上では循環は弱いです。スライダーをさらに押すと境界層は薄くなり、セルは強くなります。

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

自然対流は、ポンプなしに密度差がつくる流れです。床が熱くなると、その上の流体は膨張して軽くなります。軽い塊は浮かび上がろうとします。これが浮力です。

浮力だけなら流れはすぐ始まります。粘性は浮かぶ塊を摩擦で捕まえます。熱拡散は熱い塊の温度を消し、浮力が頼る温度差まで消します。

二つのブレーキが浮力に勝てば流体は静止します。浮力が勝てば対流が点きます。

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

ブシネスク近似は、密度を浮力項の中だけで温度の関数にします。

体膨張係数は、温度が 1 度上がったとき密度が減る割合です。

\beta \equiv -\frac{1}{\rho}\left(\frac{\partial \rho}{\partial T}\right)_P

β は膨張係数、ρ は密度、T は温度です。理想気体なら β は 1/T です。

密度差は温度差に変わります。

\frac{\rho_\infty - \rho}{\rho} \approx \beta\,(T - T_\infty)

ρ∞ と T∞ は、遠くの静止流体の密度と温度です。この一行が、鉛直平板の境界層の運動量式に浮力を埋めます。

u\frac{\partial u}{\partial x} + v\frac{\partial u}{\partial y} = g\beta\,(T - T_\infty) + \nu\frac{\partial^2 u}{\partial y^2}

u と v は速度、g は重力加速度、ν は動粘性係数です。右辺の第 1 項が浮力、第 2 項が粘性です。エネルギー式は温度を拡散で運びます。

u\frac{\partial T}{\partial x} + v\frac{\partial T}{\partial y} = \alpha\frac{\partial^2 T}{\partial y^2}

α は熱拡散係数です。

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

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

\mathrm{Gr} = \frac{g\beta\,(T_s - T_\infty)\,L^3}{\nu^2} \sim \frac{\text{浮力}}{\text{粘性力}}

Ts は壁温、L は特性長さです。グラスホフ数は、強制対流で長さで流れの居場所を分けるレイノルズ数がしていた仕事を引き受けます。

相手は粘性だけではありません。熱拡散も温度差を消します。その第二のブレーキはプラントル数、運動量拡散と熱拡散の比 Pr = ν/α を掛けて入ります。積がレイリー数です。

\mathrm{Ra} = \mathrm{Gr}\cdot\mathrm{Pr} = \frac{g\beta\,(T_s - T_\infty)\,L^3}{\nu\,\alpha}

分母に ν と α が一緒にいます。対流が点くかは Gr 一つではなく Ra が決めます。

1708という閾値#

下から温められる層を線形安定性解析に入れると、臨界値が出ます。上下が固体壁なら Ra_c は 1708 です。それより小さいと、擾乱は粘性と熱拡散に食われて消えます。流体は静止し、熱は伝導だけで動きます。

1708 を越えると、並んで逆に回るロール、ベナールセルが点きます。セルの幅は層の厚さの約 2 倍です。前の Ra スライダーがその戸をまたぎます。

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

壁がどれだけ冷えるかはヌセルト数 Nu = hL/k、伝導に対する対流です。鉛直等温平板にはチャーチル–チューの式があります。60 度の壁を 20 度の空気の中に、高さ 0.30 m で立てます。

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

Pr = 0.69Ra marker = 1e7.0

cyan: Churchill–Chu · green: laminar 0.59 Ra^1/4 · pink: turbulent 0.10 Ra^1/3

Pr のスライダーを液体金属側へ下げると曲線は下がり、油側へ上げると上がります。Ra マーカーを 1e9 の向こうへ押すと、傾きが 1/4 から 1/3 に替わります。

出力に出た数は Ra 7.03e7、Nu 54.9、h 4.94、熱流束 197.6 です。膜温度は 40 に 273.15 を足した 313.15 K で、膨張係数はその逆数、約 0.003193 です。粘性と熱拡散の積は 1.85e-5 掛ける 2.60e-5、すなわち 4.81e-10 です。高さ 0.30 m の 3 乗は 0.027 です。この断片がコードの重力 9.81 と温度差 40 に会って、印刷された 7.03e7 になります。Nu 54.9 に 0.027 を掛けて 0.30 で割ると 4.941 になり、出力の h 4.94 と揃います。4.94 掛ける 40 は 197.6 なので、熱流束の行も閉じます。この Ra を Pr 0.71 で割ると Gr は約 9.90e7 です。1708 と並べる数ではありません。

同じ Ra で層流の別名、0.59 掛ける Ra の 1/4 乗は約 54.0 です。出力の 54.9 はそれより約 0.9 上にあります。乱流の別名、0.10 掛ける Ra の 1/3 乗は約 41.3 です。その Nu で h を組み直すと約 3.72、熱流束は約 149 です。197.6 の約 4 分の 3 です。

Ra が 1e9 になるには 7.03e7 の約 14.2 倍が要ります。高さは 3 乗で入るので、0.30 m は約 0.73 m になって初めてその戸に届きます。それより手前のこの平板は 1/4 の指数側です。同じ空気、同じ 40 度で Ra が 1708 になる長さは約 8.7 mm です。その長さは密閉層の厚さです。0.30 m の平板の戸として 1708 を読むと、幾何が違います。

出力の 197.6 は平方メートルあたりの熱流束です。幅を高さと同じ 0.30 m にすると面積は 0.09 m^2 で、板 1 枚が失う熱は 197.6 掛ける 0.09、約 17.8 W です。

層流の別名だけをたどると、Nu は Ra の 1/4 乗、Ra は高さの 3 乗なので、h は高さの -1/4 乗で動きます。高さを 2 倍にすると h は 2 の -1/4 乗、約 0.841 倍になります。4.94 は約 4.15 になり、熱流束は約 166 まで減ります。1/3 の指数ではその減少が消え、h は高さに残りません。この出力の Ra は、まだその指数ではありません。

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

形ごとに式は違います。

鉛直平板: $\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}$ が $\mathrm{Ra}^{1/3}$ に比例し $\mathrm{Ra}$ が $L^3$ に比例するので、 $h = \mathrm{Nu}\,k/L$ から $L$ が落ちます。乱流自然対流の熱伝達率は大きさに縛られません。上の 0.30 m の平板は、まだその相殺を使う高さではありません。

7.03e7に乱流指数を載せると熱流束が削れます#

1/3 をこの出力に載せると Nu 54.9 は約 41.3 まで下がり、197.6 は約 149 まで削れます。高さを式から消す指数なので、0.30 m も一緒に消えます。層流の枝では高さを変えると h が追って動きます。指数を先に選ぶと、その動きが間違った理由で死にます。

1e9 は約 0.73 m で初めて来ます。セルが点いた密閉層でモードを三つだけ残すと、この熱流束はどこまで残るのか。その問いをロールを三つの式に折った Lorenz 模型が持っています。