層流では2.0000、乱流では1.0322 — 速度水頭の前につく係数

平均速度だけ持ってベルヌーイに行くと半分が消える

同じ管に同じ流量が流れています。流量を断面積で割れば平均速度 $V$ が出ます。この $V$ で $V^2/2g$ を作り、エネルギー方程式に入れる。ここまでは誰でもやります。

ところが層流管内流では、この値は実際の運動エネルギーのちょうど半分です。流量も断面積も合って いるのに、エネルギーだけが100%ずれます。今日はその差がどこで生まれるのか、乱流ではなぜ3%まで 縮むのかを、積分して確かめます。

速度水頭を定義したとき、何を積分したのか

断面を通る運動エネルギー流束は、局所の運動エネルギーに局所の質量流束を掛けて積分した量です。

$u$ は局所の軸方向速度、$\rho$ は密度、$A$ は断面積です。被積分関数に $u$ が 3回 入ります。 2回はエネルギーから、1回は質量を運ぶところから来ました。

一方、教科書のエネルギー方程式は $\dot{m} V^2/2$ と書きます。$\dot{m} = \rho V A$ なので、これは $\rho V^3 A / 2$ です。両者の比が運動エネルギー補正係数 $\alpha$ です。運動量方程式にも同じ事情が あり、そちらは $u$ が2回入った $\beta$ になります。

どちらの係数も、プロファイルが平坦なほど1に近づきます。プロファイルの切替を層流から欠損則へ動かすと、右側の帳簿の高さの差が2から1.03付近まで折れます。

プロファイルを層流・べき乗則・速度欠損則と切り替えながら、右側の帳簿2本の高さの差を見ます。 トレーサーは局所速度で動くので、放物線では中心の点だけが遠くへ逃げ、乱流プロファイルではほぼ 一塊で進みます。この「一塊で進む度合い」がそのまま $\alpha$ です。

層流管内流で係数がちょうど2になる理由#

Hagen–Poiseuille 解は放物線です。$u = U_c \left[1 - (r/R)^2\right]$、$V = U_c/2$ です。 $s = r/R$ と置くと、分子も分母も閉じた形に落ちます。

$V^3 A = (U_c/2)^3 \pi R^2 = \pi R^2 U_c^3 / 8$ なので $\alpha = 2$ です。同じ手順で $\beta = 4/3$ が出ます。近似でも経験式でもなく、整数比です。

分子をもう一度開くと、2が整数で落ちる場所が見えます。$(1-s^2)^3 s = s - 3s^3 + 3s^5 - s^7$ で、0から1まで積分すると $1/2 - 3/4 + 1/2 - 1/8 = 1/8$ です。前の $2\pi R^2 U_c^3$ を掛けると分子は $\pi R^2 U_c^3/4$ です。分母 $V^3 A$ は $\pi R^2 U_c^3/8$ です。4分の1を8分の1で割るとちょうど2です。$\beta$ は次数が一つ低く、$u^2$ の積分は $4/3$ で止まります。表の1.3333はその有理数を小数第4位で切った値です。$n=7$ なら $2n^2/((n+1)(2n+1)) = 98/120 = 0.81666\ldots$ で、出力の0.8167と一致します。

$\alpha=2$ は「エネルギーの半分を落とした」という意味ではありません。平均速度という ひとつの数 では3乗積分の代わりにならない、という意味です。同じ放物線の閉じた解を速度から追うなら、 Poiseuille流の解析解が先にその断面を解いています。 あちらは速度、こちらはその速度の3乗積分です。

乱流プロファイルを指数ひとつで書く方法

乱流管内流は壁の近くだけ急に変わり、中心部はほぼ平坦です。この形を指数ひとつで近似したのが べき乗則です。

$y$ は壁からの距離、$n$ はレイノルズ数が大きくなるほど一緒に大きくなる指数です。通常は $n=7$ を 使います。このプロファイルの平均速度は積分が閉じて解析解があります。

$n=7$ なら $0.8167$ です。層流の $0.5$ と並べれば、プロファイルがどれだけ鈍くなったかがすぐ 見えます。$\alpha$ も同じ積分なので、$n$ ひとつ決まれば値が決まります。

Python でプロファイル4種類の帳簿を並べて書いた#

3つの系列（層流・べき乗則・速度欠損則）を同じ台形積分器に入れます。外部パッケージなし、標準 ライブラリだけです。半径方向40,000分割なら層流の解析解と小数第4位まで一致します。

import math
 
R, M = 1.0, 40000                      # 半径、半径方向の分割数
dr = R / M
rs = [i * dr for i in range(M + 1)]
A = math.pi * R * R
 
def ring_integral(vals):
    """f(r) を断面上で積分: ∫ f · 2πr dr (台形則)"""
    s = 0.0
    for i in range(M):
        a = vals[i] * 2.0 * math.pi * rs[i]
        b = vals[i + 1] * 2.0 * math.pi * rs[i + 1]
        s += 0.5 * (a + b) * dr
    return s
 
def ledger(u):
    """平均速度、運動エネルギー補正係数 alpha、運動量補正係数 beta"""
    vbar = ring_integral(u) / A
    al = ring_integral([x ** 3 for x in u]) / (vbar ** 3 * A)
    be = ring_integral([x ** 2 for x in u]) / (vbar ** 2 * A)
    return vbar, al, be
 
def poiseuille():
    return [1.0 - (r / R) ** 2 for r in rs]
 
def power_law(n):
    return [max(1.0 - r / R, 0.0) ** (1.0 / n) for r in rs]
 
def defect_law(ustar):
    """u/Vbar = 1 + (u*/Vbar)(3.75 - 2.5 ln(R/y)),  y = R - r"""
    out = []
    for r in rs:
        y = max(R - r, 1e-9)
        out.append(1.0 + ustar * (3.75 - 2.5 * math.log(R / y)))
    return out
 
print(f"{'profile':<20}{'Vbar/Uc':>9}{'alpha':>9}{'beta':>8}")
vb, al_lam, be = ledger(poiseuille())
print(f"{'Poiseuille (laminar)':<20}{vb:9.4f}{al_lam:9.4f}{be:8.4f}")
for n in (6, 7, 8, 10):
    vb, al, be = ledger(power_law(n))
    exact = 2.0 * n * n / ((n + 1.0) * (2.0 * n + 1.0))
    print(f"{'power law n=' + str(n):<20}{vb:9.4f}{al:9.4f}{be:8.4f}   Vbar exact {exact:.4f}")
 
print()
cases = [("laminar", al_lam)]
print(f"{'Re':>8}{'f':>8}{'u*/Vbar':>9}{'Uc/Vbar':>9}{'alpha':>8}{'beta':>8}")
for Re in (1e4, 1e5, 1e6, 1e7):
    f = (0.790 * math.log(Re) - 1.64) ** -2.0        # Petukhov の摩擦係数
    us = math.sqrt(f / 8.0)                          # u*/Vbar
    u = defect_law(us)
    vb, al, be = ledger(u)
    print(f"{Re:8.0e}{f:8.4f}{us:9.4f}{max(u) / vb:9.4f}{al:8.4f}{be:8.4f}")
    cases.append((f"turbulent Re={Re:.0e}", al))
 
print()
V, g = 3.0, 9.81                                     # 3 m/s の管内流、速度水頭の比較
base = V * V / (2 * g)
for name, al in (cases[0], cases[1], cases[3]):
    print(f"{name:<21}alpha={al:6.4f}  alpha·V²/2g = {al * base:6.4f} m  "
          f"alpha=1 missing {al * base - base:6.4f} m ({(al - 1) * 100:6.2f}%)")

profile               Vbar/Uc    alpha    beta
Poiseuille (laminar)   0.5000   2.0000  1.3333
power law n=6          0.7912   1.0768  1.0269   Vbar exact 0.7912
power law n=7          0.8167   1.0584  1.0204   Vbar exact 0.8167
power law n=8          0.8366   1.0460  1.0160   Vbar exact 0.8366
power law n=10         0.8658   1.0307  1.0106   Vbar exact 0.8658
 
      Re       f  u*/Vbar  Uc/Vbar   alpha    beta
   1e+04  0.0315   0.0627   1.2353  1.0835  1.0309
   1e+05  0.0180   0.0474   1.1779  1.0490  1.0177
   1e+06  0.0116   0.0381   1.1430  1.0322  1.0114
   1e+07  0.0081   0.0319   1.1195  1.0227  1.0080
 
laminar              alpha=2.0000  alpha·V²/2g = 0.9174 m  alpha=1 missing 0.4587 m (100.00%)
turbulent Re=1e+04   alpha=1.0835  alpha·V²/2g = 0.4970 m  alpha=1 missing 0.0383 m (  8.35%)
turbulent Re=1e+06   alpha=1.0322  alpha·V²/2g = 0.4735 m  alpha=1 missing 0.0148 m (  3.22%)

欠損則のスライダーを0から100%まで動かすと、レイノルズ数ごとにばらけていた曲線が $u^*$ の1本に重なります。

$Re=10^4$ の $\alpha=1.0835$ と $Re=10^6$ の $\alpha=1.0322$ は、重なった線の上で $f$ だけを変えて積分し直した値です。曲線の本数を増やしても、重なった線は太くなりません。

数値積分が層流で $2.0000$ と $1.3333$ をそのまま再現しました。べき乗則の平均速度も解析解と4桁 まで同じです。積分器が正しいということです。

Re が消える座標 — 速度欠損則#

表の2番目のブロックは、べき乗則ではなく速度欠損則で作りました。対数法則を中心線で1回、任意の 点で1回書いて引くと、粘性係数がまるごと消えます。

$u^* = \sqrt{\tau_w/\rho}$ は摩擦速度です。右辺には $\nu$ も $Re$ もありません。壁からどれだけ 離れているかだけが残ります。ここに断面平均を取ると $U_c/u^* - V/u^* = 3.75$ という定数がひとつ 落ち、$u^*/V = \sqrt{f/8}$ で摩擦係数とつながります。粘性係数 $\nu$ は引く過程でまるごと抜け、 壁からの距離だけが残ります。上で重ねた線が、この定数3.75です。

$\alpha$ と $\beta$ を同じ表に載せる#

乱流側の4行を $f$ で割ると規則が見えます。$\alpha - 1$ は $f$ の2.7〜2.8倍、$\beta - 1$ は $f$ の0.98倍です。つまり2つの係数はプロファイルの形ではなく 摩擦係数ひとつ に縛られて います。粗い管ほど $f$ が大きく、$\alpha$ も一緒に大きくなります。滑らかな管を高レイノルズで 回せば $\alpha$ は1.02まで下がります。

$Re=10^6$ の行を直接割ります。$\alpha-1 = 0.0322$、$f = 0.0116$ なので $0.0322/0.0116 \approx 2.78$ です。$Re=10^4$ は $0.0835/0.0315 \approx 2.65$ です。2.7〜2.8倍という話は、この2行から出ています。$\beta$ はもっと近いです。$Re=10^6$ で $0.0114/0.0116 \approx 0.98$、$Re=10^4$ でも $0.0309/0.0315 \approx 0.98$ です。運動量の補正は摩擦係数とほぼ同じ大きさです。

3 m/s、$g = 9.81$ なら $\alpha=1$ の速度水頭は $V^2/2g = 9/19.62 = 0.4587$ m です。層流は2を掛けて0.9174 m になり、落ちた0.4587 m がその半分です。$Re=10^6$ の1.0322を同じ底に掛けると0.4735 m、差は0.0148 m です。$Re=10^4$ では落ちる量が0.0383 m、8.35%です。損失水頭が数メートルの配管では0.0148 m は埋もれ、層流の0.4587 m は帳簿に残ります。

べき乗則と欠損則が同じレイノルズ数で少し違う値を出すことも表に見えます。べき乗則は壁で勾配が 無限大、中心線で折れます。2つの近似は、近似である分だけ食い違います。

この係数が実際に数字を動かす場所

最後のブロックが3 m/s の管内流の速度水頭です。乱流で $\alpha$ を1にすると、$Re=10^6$ で 0.0148 m を落とします。配管1区間の損失水頭が数メートルなら無視できます。層流では0.4587 m を 落とします。同じ計算で100%です。

数字が動く場所は分かれています。粘度の大きい流体や細い管は層流なので $\alpha=2$ です。入口領域は プロファイルがまだ育っている途中で、$\alpha$ は1と2の間を通ります。出口面の全圧を面積平均で取るか 質量流束重みで取るかは、$\int u^2$ と $\int u^3$ のどちらを使うかという問題です。2つの値が違って 出たらソルバーが間違ったのではなく、平均の定義が違うのです。

$\alpha$ は定数のように書かれていますが、決めるのはプロファイルです。手で1を打ち込むと、層流の帳簿が半分ずれます。

教科書が $\alpha=1$ と書いて先へ進む条件#

$\alpha=1$ はプロファイルが完全に平坦なときだけ成り立ちます。乱流でこの近似が通るのは、乱流が 平坦だからではなく $f$ が小さいからです。$f$ が0.01あたりなら $\alpha$ は1.03で、3%は配管計算の 他の不確かさに埋もれます。

ですから判断基準は「乱流かどうか」ではなく「$f$ がいくつか」です。遷移が管内の層を何枚にも 分ける過程はレイノルズ臨界と管内遷移にあります。 粗い管、低いレイノルズ数、層流 — どれかひとつでも当てはまれば $\alpha$ を積分し直します。計算は 上のコードのように1回です。