v の項と v² の項を足したら両端は合い、真ん中で 42% 外れた

線形抵抗と二乗抵抗の和は Re≪60 と Re≫60 では合いますが、Re 10〜300 では後流が生む追加抵抗が抜け、終端速度が 30〜40% 大きく出ます。

同じ水、同じガラス玉、違う抵抗則

水の中にガラス玉を二つ落とします。一つは直径 0.05 mm、もう一つは 10 mm です。小さい玉が受ける抵抗は速度に比例します。大きい玉が受ける抵抗は速度の二乗に比例します。水も材質も同じです。違うのは大きさだけです。

18 世紀末のフランスの学者たちは、この二つの法則を別々に見つけました。そして両者を足した式 $F = a v^2 + b v$ にまとめました。この記事では、その和がどこで合い、どこで外れるかを数字で確かめます。読み終えると、二つの項の主役が入れ替わる Reynolds 数(慣性力と粘性力の比)と、単純な和がいちばん大きく外れる区間がわかります。

五つの問いに分けて進めます。

Q1. 1775 年から 1800 年まで、抵抗は v だったのか v² だったのか#

1752 年、ダランベールは理想流体では物体の抵抗がゼロになるというパラドックスに突き当たりました。理論が行き詰まると、実験が道を開きました。1775 年、フランス科学アカデミーは流体抵抗の実験チームを組みました。実務は Bossut が担当しました。

Bossut は模型船を水槽で引きました。精密に測ると、抵抗は速度の二乗に比例していました。一方、Du Buat はパイプ内の流れを測りました。彼が見た摩擦抵抗は、速度にそのまま比例していました。二つの実験は別々の答えを出したのです。

Prony はさらに精密な実験で両者を一つにまとめました。抵抗は $a v^2 + b v$ と書けます。Coulomb はこれを二種類の抵抗として整理しました。速度の二乗に比例するのが形状抵抗(form drag、物体の形がつくる圧力差)、速度に比例するのが摩擦抵抗(friction drag)です。この区別は後に Navier と Stokes の粘性流れの研究へつながります。

どちらが正しかったのでしょうか。どちらも正しかったのです。模型船は大きく速く、パイプ内の流れは遅く狭いものでした。二つの実験は同じ物理の両端を見ていました。その「端」を分ける数字が次の問いです。

Q2. 二つの項はどこから来るのか#

球に働く抵抗を二項で書くと次のようになります。

F = \underbrace{3\pi\mu d\, v}_{b\,v} + \underbrace{C_\infty \frac{\pi}{8}\rho d^2 v^2}_{a\,v^2}

$\mu$ は粘性係数、 $\rho$ は流体密度、 $d$ は球の直径、 $C_\infty$ は高速域の抗力係数(球ではおよそ 0.4)です。

第一項は Stokes 抵抗です。慣性がまったくない遅い流れ(クリーピング流れ)の厳密解です。第二項は Newton 抵抗です。前面にぶつかった流体の運動量がそのまま圧力に変わると考えると出てきます。

一つ注意があります。Stokes 抵抗 $3\pi\mu d v$ は純粋な表面摩擦ではありません。その 2/3 が摩擦で、1/3 は圧力抵抗です。ですから Coulomb の「摩擦対形状」より、「粘性支配対慣性支配」のほうが正確な呼び名です。

下のシミュレーションで実際に操作してみましょう。玉を落とすと、速度が上がるにつれて二つの抵抗が育ちます。その和が重さと等しくなる点が終端速度(terminal velocity、加速が止まる速度)です。

d = 0.30 mm fluid:

v = 0.00 mm/s · v_t = 1000 mm/s · Re = 0.00 · linear share 100.0% · quadratic 0.0%

直径スライダーを 0.02 mm から 20 mm まで動かしてみてください。小さい玉では青い直線(b·v)が重さの線と交わり、大きい玉ではオレンジの放物線(a·v²)が交わります。流体を glycerol に替えると、10 mm の玉でも Du Buat 側へ移ります。

Q3. 主役はいつ入れ替わるのか#

二つの項が等しくなる速度を求めればよいのです。

3\pi\mu d\, v = C_\infty \frac{\pi}{8}\rho d^2 v^2 \quad\Longrightarrow\quad Re \equiv \frac{\rho v d}{\mu} = \frac{24}{C_\infty}

$Re$ は球の直径を基準にした Reynolds 数です。 $C_\infty = 0.4$ なら $Re = 60$ です。

この結果が Bossut と Du Buat を和解させます。抵抗則を決めるのは流体の種類でも、速度そのものでもありません。 $Re$ ひとつです。パイプ内の遅い流れは $Re$ が小さく、模型船は $Re$ が数十万でした。先ほどのシミュレーションで水と glycerol を入れ替えたのも、結局は $Re$ を動かしたことになります。

「二つの極限をつなぐ」同じ形の問題はほかにもあります。粒子熱伝達の Knudsen 遷移領域では、連続体と自由分子の二つの極限を調和平均でつなぎました。Wood の音速では、圧縮率を体積分率で足し合わせました。極限の混ぜ方が結果を変えます。Prony が選んだのは単純な足し算でした。

Q4. では足すだけでよいのか#

確かめてみましょう。水中のガラス玉( $\rho_p = 2500$ kg/m³)を六つ、静止状態から放します。二項抵抗で運動方程式を積分します。

m \frac{dv}{dt} = (\rho_p - \rho) V g - b v - a v^2

$m$ は球の質量、 $V$ は体積、 $g$ は重力加速度です。付加質量(added mass、球が引きずる周囲流体の慣性)は無視しました。

終端速度を基準値と比べます。基準は Schiller–Naumann 相関式 $C_D = \frac{24}{Re}(1 + 0.15 Re^{0.687})$ です。 $Re < 1000$ で実験と数 % 以内で合うことが知られています。

import math
 
# 水中を沈むガラス玉:Du Buat(線形)対 Bossut(二乗)の抵抗
RHO_P, RHO, MU, G = 2500.0, 1000.0, 1.0e-3, 9.81
C_INF = 0.4          # 二項則の二乗(形状)抵抗係数
 
def weight(d):
    return (RHO_P - RHO) * math.pi * d**3 / 6 * G
 
def prony_drag(d, v):
    lin = 3 * math.pi * MU * d * v                 # Du Buat / Stokes:  b v
    quad = C_INF * math.pi / 8 * RHO * d**2 * v**2  # Bossut / Newton:  a v^2
    return lin, quad
 
def vt_prony(d):
    a = C_INF * math.pi / 8 * RHO * d**2
    b = 3 * math.pi * MU * d
    return (-b + math.sqrt(b * b + 4 * a * weight(d))) / (2 * a)
 
def cd_schiller(re):
    return 24 / re * (1 + 0.15 * re**0.687) if re < 1000 else 0.44
 
def vt_reference(d):
    lo, hi = 1e-9, 10.0
    for _ in range(200):
        v = 0.5 * (lo + hi)
        re = RHO * v * d / MU
        f = cd_schiller(re) * 0.5 * RHO * v * v * math.pi * d * d / 4
        lo, hi = (v, hi) if f < weight(d) else (lo, v)
    return 0.5 * (lo + hi)
 
def settle(d, steps_per_tau=400):
    m = RHO_P * math.pi * d**3 / 6
    vt = vt_prony(d)
    tau = vt / ((1 - RHO / RHO_P) * G)       # 加速の時間スケール
    dt = tau / steps_per_tau
    v, t = 0.0, 0.0
    while v < 0.99 * vt:
        lin, quad = prony_drag(d, v)
        v += dt * (weight(d) - lin - quad) / m
        t += dt
    return t, tau, vt
 
print(" d[mm]   Re_t    vt[mm/s]  ref[mm/s]  err%   lin%   t99/tau  Stokes/vt  Newton/vt")
for dmm in [0.05, 0.1, 0.3, 1.0, 3.0, 10.0]:
    d = dmm * 1e-3
    t99, tau, vt = settle(d)
    lin, quad = prony_drag(d, vt)
    ref = vt_reference(d)
    v_st = weight(d) / (3 * math.pi * MU * d)
    v_nw = math.sqrt(weight(d) / (C_INF * math.pi / 8 * RHO * d * d))
    re = RHO * vt * d / MU
    print(f"{dmm:6.2f} {re:7.2f} {vt*1e3:10.2f} {ref*1e3:10.2f} {100*(vt-ref)/ref:6.1f} "
          f"{100*lin/(lin+quad):6.1f} {t99/tau:8.2f} {v_st/vt:10.2f} {v_nw/vt:10.2f}")
re_x = 24 / C_INF  # 3 pi mu d v = C_inf pi/8 rho d^2 v^2
print(f"crossover Re where lin = quad: 24/C_inf = {re_x:.1f}")

 d[mm]   Re_t    vt[mm/s]  ref[mm/s]  err%   lin%   t99/tau  Stokes/vt  Newton/vt
  0.05    0.10       2.04       1.98    2.9   99.8     4.59       1.00      24.27
  0.10    0.81       8.07       7.29   10.6   98.7     4.55       1.01       8.68
  0.30   17.16      57.21      40.21   42.3   77.8     3.92       1.29       2.12
  1.00  193.49     193.49     145.95   32.6   23.7     2.93       4.22       1.14
  3.00 1121.20     373.73     365.75    2.2    5.1     2.70      19.69       1.03
 10.00 6973.63     697.36     667.76    4.4    0.9     2.65     117.23       1.00
crossover Re where lin = quad: 24/C_inf = 60.0

上のコードをブラウザで直接動かしてみましょう。同じ玉が二つ並んで沈みます。一つは Prony の二項抵抗、もう一つは Schiller–Naumann の基準抵抗を受けます。

d = 0.30 mm C∞ = 0.40 ρp = 2500 kg/m³ steps/frame 8

step 0 · Prony v = 0.00 · ref v = 0.00 mm/s · t99/τ = …
Re_t 17.16 · vt 57.21 · ref 40.21 · err 42.3% · lin 77.8% · Stokes/vt 1.29 · Newton/vt 2.12

既定値 d = 0.30 mm では err 42.3%、lin 77.8%、t99/τ 3.92 となり、出力ブロックと一致します。右の曲線は直径全体にわたる誤差です。C∞ を 0.6 に上げると大きい玉側の誤差が負に転じ、ρp を下げると同じ直径でも Re_t が下がって山が右へずれます。

出力から三つのことが読み取れます。第一に、両端はよく合います。0.05 mm で 2.9%、3 mm で 2.2% です。第二に、片方の項だけだと反対の端で崩れます。Stokes だけなら 10 mm の玉が 117 倍速くなり、Newton だけなら 0.05 mm の玉が 24 倍速くなります。第三に、問題は真ん中です。0.3 mm( $Re_t = 17$ )で 42%、1 mm( $Re_t = 193$ )で 33% 外れます。

加速時間も法則に従って変わります。線形抵抗だけなら速度は指数関数的に近づくので、 $t_{99}/\tau = \ln 100 = 4.61$ です。二乗抵抗だけなら $\tanh$ で近づくので、 $t_{99}/\tau = \operatorname{artanh} 0.99 = 2.65$ です。出力の 4.59 と 2.65 がこの二つの極限で、中間の直径はその間を通ります。

Q5. なぜ真ん中で外れるのか#

二項則を抗力係数に書き直すと、原因が見えてきます。

C_D^{\text{Prony}} = \frac{24}{Re} + C_\infty, \qquad C_D^{\text{S–N}} = \frac{24}{Re}\left(1 + 0.15\, Re^{0.687}\right)

どちらも $Re \to 0$ で $24/Re$ に近づきます。違いは括弧の中の第二項です。

$Re = 17$ を入れてみます。Prony は $1.41 + 0.40 = 1.81$ 、Schiller–Naumann は $1.41 \times 2.05 = 2.90$ を与えます。抵抗が基準より 38% 小さいのです。重さは同じなので、Prony の玉はより速く落ちてやっと釣り合います。その結果が 42% 大きい終端速度です。

足りない抵抗はどこから来るのでしょうか。 $Re$ が 1 を超えると、前後対称だった流れが崩れます。後ろに遅い後流ができ、 $Re \approx 20$ 付近で表面の流れが剥がれて再循環領域が生まれます。この後流の圧力抵抗は、粘性と慣性が混ざって生じます。 $v$ にも $v^2$ にもきれいには属しません。Prony の足し算には、この混ざった項の居場所がないのです。

だから二つの項が同じくらいの大きさのとき、ちょうど $Re \approx 60$ の近くで、足し算がいちばん大きく外れます。歴史上の二つの実験がそれぞれ正しかった理由も同じです。どちらもこの区間の外にありました。

二項則がいまも生きている場所

外れる区間があっても、 $a v^2 + b v$ の形は消えませんでした。多孔質媒体の Ergun 式と Darcy–Forchheimer 式はまさにこの形です。パイプ損失の Prony 式も同じ根を持ちます。多孔質媒体では、空隙スケールの $Re$ が広い範囲にわたっても二項の和がかなりよく合います。係数を実験に合わせて決めるからです。

Euler–Lagrange 粒子コードでは事情が違います。粒子一つの $Re$ が、ちょうどこの問題の区間(1〜1000)に入ることが多いのです。噴霧液滴、流動層の粒子、沈降粒子がそうです。そこで $24/Re + 0.4$ を使うと、終端速度が 30〜40% ずれます。沈降時間や滞留時間も同じだけずれます。

コードの抗力モデルを選ぶときは、次の順で確認するとよいでしょう。

まず粒子 $Re$ の分布範囲を取り出します。1 未満なら Stokes、1000 超なら一定の $C_D$ で十分です。
1〜1000 の区間に質量が多くかかっているなら、Schiller–Naumann のように中間項を持つ相関式を使います。
単一粒子の沈降テストを一つ回し、上の表の 0.3 mm と 1 mm の行と比べます。誤差 30% が出たら、二項の足し算が隠れているサインです。

Bossut の船と Du Buat のパイプは、同じ水を違う $Re$ で見ていました。その間の空白は、足し算だけでは埋まりません。

v の項と v² の項を足したら両端は合い、真ん中で 42% 外れた — Prony の二項抵抗則と Re=60