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

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

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

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

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

下のシミュレーションを実際に操作してみましょう。

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

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

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

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

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

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

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

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

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

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

方程式が双子になる瞬間#

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

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

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

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

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

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

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

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

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

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

これが Chilton-Colburn の相似則(修正レイノルズ相似則)です。$Pr = 1$ なら元の式に戻ります。下のグラフを操作してみましょう。

オレンジ色の曲線(生の $St$)は $Pr$ に応じて上下に動きます。しかし水色の曲線($St \cdot Pr^{2/3}$)は $Pr$ が何であっても破線($C_f/2$)の上にぴたりと張り付いています。補正指数 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$ となります。

覚えておきたいこと#

• 速度境界層と温度境界層は同じ形の無次元方程式に従います。違いは拡散係数 $Re$ と $Re \cdot Pr$ だけです。
• $Pr = 1$、圧力勾配ゼロなら速度場と温度場が重なり、$C_f/2 = St$ となります。これがレイノルズの相似則です。
• 一般の流体では Chilton-Colburn の相似則 $C_f/2 = St \cdot Pr^{2/3}$ で摩擦から熱伝達を推定します。