抗力を測れば熱伝達がわかる — レイノルズの相似則

Pr=7の水に補正なしのCf/2=Stを使うと、熱伝達率が約3.7倍に膨らみます。

風洞で模型の抗力だけを測り、その表面がどれほど速く冷えるかを計算できます。とんでもない話に聞こえますが、航空機の翼やタービン翼の熱設計に実際に使われている方法です。摩擦と熱伝達が同じ形の方程式に従うからです。相似がどこから来るかは、平板のコードで合わせてみました。残るのは、表面摩擦係数ひとつから熱伝達係数を推定する式です。

二つの境界層が並んで成長する

流体が壁の上を流れると、二種類の薄い層ができます。一つは速度が 0 から主流値へ回復する速度境界層(momentum boundary layer)です。もう一つは温度が壁温度から主流温度へ変わる温度境界層(thermal boundary layer)です。

二つの層は同じ場所から始まり、一緒に厚くなります。しかし厚さがいつも等しいわけではありません。どちらが速く広がるかを決めるのが Prandtl 数(運動量拡散率と熱拡散率の比)です。

Pr のスライダーを 0.1 から 10 へ動かすと、温度層と速度層の厚さの比が一緒に変わります。log10 Re_L のスライダーは、二つの層の絶対厚さだけを縮めます。

log₁₀ Re_L: 5.0  (Re_L = 1.0e+5)Prandtl number Pr: 1.00  (δ_T/δ = 1.00)C_f/2 = 1.050e-3  |  St = 1.050e-3
St · Pr^(2/3) = 1.050e-3  ← Pr≈1: C_f/2 = St

Pr を 1 に合わせると、二つの曲線がぴたりと重なります。Pr を 10 まで上げると、オレンジ色の温度境界層が青い速度境界層の内側へと薄くなります。この重なりが相似則の出発点です。

無次元化すると Re と Pr だけが残る#

平板上の定常二次元層流境界層では、運動量方程式とエネルギー方程式は次のように簡約されます。

u \frac{\partial u}{\partial x} + v \frac{\partial u}{\partial y} = -\frac{1}{\rho}\frac{dP}{dx} + \nu \frac{\partial^2 u}{\partial y^2}

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

ここで $u, v$ は速度成分、 $\nu$ は動粘性係数(運動量拡散率)、 $\alpha$ は熱拡散率です。

座標と変数を無次元化します。 $x^\* = x/L$ 、 $u^\* = u/U_\infty$ とし、温度は $T^\* = (T-T_s)/(T_\infty-T_s)$ と置きます。 $T_s$ は壁温度、 $T_\infty$ は主流温度です。代入すると、両方の方程式がわずか二つの無次元数で表されます。

u^\* \frac{\partial u^\*}{\partial x^\*} + v^\* \frac{\partial u^\*}{\partial y^\*} = -\frac{dP^\*}{dx^\*} + \frac{1}{Re_L}\frac{\partial^2 u^\*}{\partial y^{\*2}}

u^\* \frac{\partial T^\*}{\partial x^\*} + v^\* \frac{\partial T^\*}{\partial y^\*} = \frac{1}{Re_L \, Pr}\frac{\partial^2 T^\*}{\partial y^{\*2}}

$Re_L = U_\infty L / \nu$ は Reynolds 数(慣性力と粘性力の比)、 $Pr = \nu/\alpha$ は Prandtl 数です。二つの式の形はほぼ同じです。違いは拡散項の係数一つだけです。

Prandtl 数 — 速度と温度、どちらが厚いか#

Prandtl 数は流体の性質です。運動量が広がる速さと熱が広がる速さの比を測ります。

Pr = \frac{\nu}{\alpha} = \frac{\text{運動量拡散率}}{\text{熱拡散率}}

層流平板では、二つの境界層厚さは次の関係に従います。

\frac{\delta}{\delta_T} \approx Pr^{1/3}

$\delta$ は速度境界層厚さ、 $\delta_T$ は温度境界層厚さです。空気は $Pr \approx 0.7$ なので二つの層は近い厚さです。水は $Pr \approx 7$ なので温度層が薄くなります。エンジンオイルは $Pr$ が数百で温度層は非常に薄くなります。水銀のような液体金属は $Pr \ll 1$ で温度層がはるかに厚くなります。

空気の $Pr \approx 0.7$ では $Pr^{1/3}$ が約 0.89 です。温度層が速度層より少し厚くなります。水の $Pr \approx 7$ では $7^{1/3} \approx 1.91$ で、温度層は半分近くまで入ります。出力の Nu がその比です。Pr 1 で 104.99、Pr 0.7 で 92.95（0.885倍）、Pr 7 で 200.50（1.91倍）です。 $C_f/2 = 1.050 \times 10^{-3}$ は Pr を知りません。摩擦は速度場の分で、熱はその上に Pr を載せるだけです。運動論で μ と k を同じ式から取ると Pr が 1.98 倍ずれる係数を見ると、Pr=1 という前提がどこで割れるか分かります。

方程式が双子になる瞬間

二つの無次元方程式をもう一度見てみましょう。圧力勾配がなく( $dP^\*/dx^\* = 0$ )、 $Pr = 1$ ならどうなるでしょうか。運動量方程式とエネルギー方程式が文字通り同じ形になります。境界条件も同じです。壁で $u^\* = 0, T^\* = 0$ 、主流で $u^\* = 1, T^\* = 1$ です。

同じ方程式に同じ境界条件なら、解も同じです。つまり $u^\*(x^\*, y^\*) = T^\*(x^\*, y^\*)$ です。速度分布と温度分布が重なります。これが相似則の核心です。二つの過程が同じ無次元方程式に従うなら、互いに交換可能なのです。

レイノルズの相似則: Cf/2 = St#

壁での勾配が等しければ、摩擦と熱伝達も結びつきます。壁せん断応力から定まる摩擦係数と、壁熱流束から定まる Nusselt 数が、同じ微分で定義されるからです。

C_f = \frac{\tau_s}{\tfrac{1}{2}\rho U_\infty^2}, \qquad Nu_L = \frac{hL}{k_f} = \left.\frac{\partial T^\*}{\partial y^\*}\right|_{y^\*=0}

$\tau_s$ は壁せん断応力、 $h$ は対流熱伝達係数、 $k_f$ は流体の熱伝導率です。 $u^\* = T^\*$ なら壁勾配が等しいので、次が成り立ちます。

\frac{C_f}{2} Re_L = Nu_L

Stanton 数(修正 Nusselt 数)を導入するともっとすっきりします。

St \equiv \frac{h}{\rho U_\infty c_p} = \frac{Nu_L}{Re_L \, Pr}

$c_p$ は定圧比熱です。これを上式に入れるとレイノルズの相似則が得られます。

\frac{C_f}{2} = St \qquad (Pr \approx 1,\ dP/dx \approx 0)

摩擦係数の半分がそのまま Stanton 数です。抗力を測れば熱伝達がわかります。

Pr が 1 でないとき — Chilton-Colburn#

ほとんどの流体は $Pr \ne 1$ です。幸い小さな補正で救えます。温度境界層の厚さが $Pr^{1/3}$ だけ違う点を反映すると、次の式が得られます。

\frac{C_f}{2} = St \cdot Pr^{2/3} \qquad (0.6 < Pr < 60)

これが Chilton-Colburn の相似則(修正レイノルズ相似則)です。 $Pr = 1$ なら $Pr^{2/3}=1$ で元の式に戻ります。補正指数 2/3 は、厚さ比 $Pr^{1/3}$ と Stanton 数の分母にある Pr が残したものです。 $St = Nu/(Re\,Pr)$ で Nu が $Pr^{1/3}$ に比例するなら、指数は $1/3-1=-2/3$ です。両辺に $Pr^{2/3}$ を掛けると Pr は式から消えます。

Python — 平板で摩擦と熱を一緒に測る#

相似則を数値で確かめましょう。Blasius 方程式で速度場を解いて $C_f$ を得ます。同じ速度場でエネルギー方程式を解いて $Nu$ を得ます。そして $St \cdot Pr^{2/3}$ が $C_f/2$ と合うか見ます。

import numpy as np
 
def blasius_profile(eta_max=10.0, n=2000):
    """Blasius 方程式 f''' + 0.5 f f'' = 0 を RK4 + シューティングで解く。"""
    deta = eta_max / n
    def rhs(y):              # y = [f, f', f'']
        f, fp, fpp = y
        return np.array([fp, fpp, -0.5 * f * fpp])
    def integrate(s):        # s = f''(0)
        y = np.array([0.0, 0.0, s])
        Y = [y.copy()]
        for _ in range(n):
            k1 = rhs(y)
            k2 = rhs(y + 0.5 * deta * k1)
            k3 = rhs(y + 0.5 * deta * k2)
            k4 = rhs(y + deta * k3)
            y = y + deta / 6 * (k1 + 2 * k2 + 2 * k3 + k4)
            Y.append(y.copy())
        return np.array(Y)
    lo, hi = 0.1, 1.0        # f'(inf)=1 を満たすよう f''(0) を二分探索
    for _ in range(60):
        s = 0.5 * (lo + hi)
        if integrate(s)[-1, 1] > 1.0:
            hi = s
        else:
            lo = s
    eta = np.linspace(0.0, eta_max, n + 1)
    return eta, integrate(s), s   # s -> f''(0) ~ 0.332
 
def thermal_slope(eta, fvals, pr):
    """エネルギー方程式 theta'' + 0.5 Pr f theta' = 0 を解いて theta'(0) を返す。"""
    deta = eta[1] - eta[0]
    def integrate(g):        # g = theta'(0)
        th, thp = 0.0, g
        for i in range(len(eta) - 1):
            fi = fvals[i]
            for _ in range(1):   # 1 ステップの RK4 (f を固定)
                k1 = (thp, -0.5 * pr * fi * thp)
                k2 = (thp + 0.5 * deta * k1[1], -0.5 * pr * fi * (thp + 0.5 * deta * k1[1]))
                k3 = (thp + 0.5 * deta * k2[1], -0.5 * pr * fi * (thp + 0.5 * deta * k2[1]))
                k4 = (thp + deta * k3[1], -0.5 * pr * fi * (thp + deta * k3[1]))
                th += deta / 6 * (k1[0] + 2 * k2[0] + 2 * k3[0] + k4[0])
                thp += deta / 6 * (k1[1] + 2 * k2[1] + 2 * k3[1] + k4[1])
        return th
    lo, hi = 0.0, 2.0        # theta(inf)=1 を満たすよう二分探索
    for _ in range(60):
        g = 0.5 * (lo + hi)
        if integrate(g) > 1.0:
            hi = g
        else:
            lo = g
    return g                 # theta'(0) ~ 0.332 Pr^(1/3)
 
Re_x = 1.0e5
eta, Y, fpp0 = blasius_profile()
Cf_half = fpp0 / np.sqrt(Re_x)          # 0.332 / sqrt(Re)
print(f"f''(0) = {fpp0:.4f}   Cf/2 = {Cf_half:.3e}")
for pr in (0.7, 1.0, 7.0):
    thp0 = thermal_slope(eta, Y[:, 0], pr)
    Nu = thp0 * np.sqrt(Re_x)
    St = Nu / (Re_x * pr)
    print(f"Pr={pr:4.1f}  Nu={Nu:7.2f}  St*Pr^(2/3)={St * pr ** (2/3):.3e}")

出力は次のようになります。

f''(0) = 0.3321   Cf/2 = 1.050e-03
Pr= 0.7  Nu= 92.95  St*Pr^(2/3)=1.050e-03
Pr= 1.0  Nu=104.99  St*Pr^(2/3)=1.050e-03
Pr= 7.0  Nu=200.50  St*Pr^(2/3)=1.050e-03

三つの Prandtl 数すべてで、 $St \cdot Pr^{2/3}$ が $C_f/2 = 1.05 \times 10^{-3}$ と一致します。 $Pr = 1$ のとき $Nu = 0.332\,Re^{1/2}$ であり、ただちに $St = C_f/2$ となります。

marker Pr のスライダーを 0.7、1、7 に置くと、水色の補正 Stanton 数が破線 $C_f/2$ に付いているかどうかがそのまま読めます。log10 Re_x を動かすと二つの曲線は一緒に下がり、Pr を動かしても水色は破線を離れません。

log₁₀ Re_x: 5.0marker Pr: 0.70St(Pr=0.70) = 1.332e-3  |  C_f/2 = 1.050e-3
St · Pr^(2/3) = 1.050e-3 (= C_f/2)

オレンジ色の曲線(生の $St$ )は $Pr$ に応じて上下に動きます。水色の曲線( $St \cdot Pr^{2/3}$ )は $Pr$ が何であっても破線( $C_f/2$ )の上に付いています。補正指数 2/3 が $Pr$ の影響を打ち消します。

三行が同じ $1.050 \times 10^{-3}$ で終わるのは、割り算が指数を戻すからです。Pr 1 では $St = 104.99/10^5 = 1.050 \times 10^{-3}$ です。 $Pr^{2/3}=1$ なので補正はありません。Pr 0.7 の生の St は $92.95/(10^5 \times 0.7) = 1.328 \times 10^{-3}$ です。 $0.7^{2/3} \approx 0.788$ を掛けると約 $1.05 \times 10^{-3}$ に戻ります。Pr 7 は逆です。生の St は $200.50/(10^5 \times 7) = 2.864 \times 10^{-4}$ です。 $7^{2/3} \approx 3.66$ を掛けて初めて $C_f/2$ と合います。

$f''(0)=0.3321$ も同じ場所です。 $C_f/2 = 0.3321/\sqrt{10^5}$ で、 $\sqrt{10^5} \approx 316.2$ なので $0.3321/316.2 = 1.050 \times 10^{-3}$ です。Blasius の壁勾配ひとつが摩擦の全部を決めます。エネルギー式は速度場 $f$ を係数として借りるだけです。温度分布が Pr で変わっても $f''(0)$ はそのままです。

水で Pr の 2/3 乗を外すと Stanton 数が約 3.7 倍に膨らみます#

Pr 7 で生の St $2.864 \times 10^{-4}$ と $C_f/2 = 1.050 \times 10^{-3}$ の比は 3.67 です。 $7^{2/3}$ と同じです。水にレイノルズの相似則 $C_f/2 = St$ をそのまま使うと、熱伝達率を約 3.7 倍大きく書いてしまいます。空気は $0.7^{2/3} \approx 0.79$ なので、同じ誤りは二割の内側に収まります。式に書いてある窓は Pr 0.6 から 60 までです。

圧力勾配がゼロでなければ、圧力項は運動量式にだけ残り、エネルギー式にはありません。二つの式は同じ解を持ちません。はく離が起きると壁面せん断はほとんど消え、熱流束はその比につきません。粗い球で抗力危機が表面摩擦を急に落とすところを見ると、 $C_f$ だけを信じて熱を外挿すると数字がどこで割れるかが分かります。