直径を1%読み違えると流量は4%ずれる — 壁の半セルとPoiseuilleの定数

粘性スキームが完璧でも、壁を半セル内側に置けば流量はその4倍ずれます。

検証ケースが4%ずれたのに、スキームは正常だった#

円管層流は解析解を持つ数少ない検証ケースです。新しく書いた粘性項の離散化をここに掛けたところ、流量が4%低く出ました。まずスキームを疑うのが順序です。ところが格子を2倍に細かくしても、誤差は半分にならずそのまま残りました。収束しない誤差は離散化誤差ではなく、幾何の誤差です。

犯人は壁が立っている位置でした。半径を1%内側に取ると、流量は4%落ちます。この記事では、その倍率4がどこから来るのか、そして同じ指数4を1838年のPoiseuilleがガラス管の実験で逆向きに使ったことを追いかけます。

流量は直径の4乗に乗っている#

完全発達した円管流れで残るのは、軸方向の運動量バランス1本だけです。

\frac{\mu}{r}\frac{d}{dr}\left(r\,\frac{du}{dr}\right) = \frac{dp}{dx}

$\mu$ は粘性係数、 $r$ は中心軸から測った半径座標、 $u$ は軸方向速度です。中心で $du/dr = 0$ 、壁で $u = 0$ を課すと、放物線が1つ出てきます。

u(r) = \frac{\Delta p}{4\mu L}\left(R^{2} - r^{2}\right)

$\Delta p$ は長さ $L$ の区間での圧力降下、 $R$ は管半径です。これを断面で積分すると流量になります。

Q = \int_{0}^{R} u(r)\,2\pi r\,dr = \frac{\pi\,\Delta p\,D^{4}}{128\,\mu L}

$D = 2R$ は管の直径です。積分で $r\,dr$ が一度余分に掛かるため、 $R^{2}$ が $R^{4}$ になります。4乗の出どころはここだけです。

この指数が誤差に掛かると倍率になります。

\frac{\delta Q}{Q} = 4\,\frac{\delta D}{D} + \mathcal{O}\!\left(\left(\frac{\delta D}{D}\right)^{2}\right)

$\delta D$ は直径の誤差、 $\delta Q$ はそれによる流量の誤差です。壁位置の誤差は、ほかの検証量では普通1次で現れます。流量だけが4倍にして返してきます。

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

dD/D-2.5%dp/dx1.0

Drag dD/D to 1 % and watch the red bar: the wall moved by 1 %, the flow rate moved by 4.06 %. The tracer counters at the right are an independent measurement — they know nothing about the formula, and they still settle on 0.904 (currently 1.000 after 0 particles). Switch to wall pinned at cell centre and lower N: a mesh of 10 cells loses 19 % of the flow with a perfectly correct viscous scheme.

dD/D のスライダーで下側の管の壁だけを動かし、右側の赤いバーとカウンターを見てください。カウンターは公式を知らないまま出口を通過した粒子を数えるだけの独立した測定ですが、それでも $(1+e)^{4}$ に収束します。wall pinned at cell centre に切り替えてNを10まで下げると、粘性スキームが完全なまま流量が19%消えます。

1822年にNavierが空けた席、1838年にPoiseuilleが埋めた指数#

Navierは橋梁の技術者でした。流体抵抗は実測では明らかなのにEuler方程式には表現がない、という点に彼はとらわれていました。1822年、まだ弾性論の数学的な基礎が固まる前に、彼は速度のラプラシアンの形で粘性項を導き、Euler方程式に付け加えます。ダランベール以降半世紀以上停滞していた流体力学が、再び動き出した地点です。

同じ時期のÉcole PolytechniqueにはCauchyが教授として、Coriolisが講師としていました。回転座標系のエネルギー方程式が水車の実験から出てきたのも、この学校の同じ空気の中です。1816年に学校が一時閉鎖されると、入学したばかりのPoiseuilleは進路を変えて医学部へ進みました。

彼は医師になり、血流を測ります。血圧が血管の直径でどう変わるかが彼の問いでした。細い血管を模すためにガラス管を自分で引き、最も細いものは内径0.015 mm。髪の毛の5分の1の太さです。1838年に彼が発表したのは、こういう形でした。

Q = K''\,\frac{\Delta p\,D^{4}}{L}, \qquad K'' = \frac{\pi}{128\,\mu}

右側の $K''$ は後から埋められた値です。Poiseuille自身は粘度という概念を持っていませんでした。彼は $K''$ を単なる定数として書きました。実験が確定させたのは定数ではなく、指数だったのです。

Pythonで測った「壁の半セル」の値段#

上の運動量バランスを半径方向の有限体積で解きました。壁のすべりなし条件を壁面に課す場合と、最後のセル中心にそのまま固定する場合を並べて回します。

import math
 
MU, GRAD_P, RADIUS = 1.0e-3, 100.0, 1.0e-3   # Pa*s, Pa/m, m
 
 
def poiseuille_q(radius, mu=MU, grad_p=GRAD_P):
    """解析流量  Q = pi*G*R^4/(8 mu)"""
    return math.pi * grad_p * radius ** 4 / (8.0 * mu)
 
 
def solve_pipe_fv(ncell, wall_at_face=True, radius=RADIUS, mu=MU, grad_p=GRAD_P):
    """軸対称の完全発達速度分布、環状セルncell個の有限体積"""
    dr = radius / ncell
    rf = [i * dr for i in range(ncell + 1)]          # 面の半径
    lo, dg, up, rhs = [0.0] * ncell, [0.0] * ncell, [0.0] * ncell, [0.0] * ncell
    for i in range(ncell):
        rhs[i] = -grad_p * (rf[i + 1] ** 2 - rf[i] ** 2) / 2.0
        if i > 0:
            w = mu * rf[i] / dr
            lo[i], dg[i] = w, dg[i] - w
        if i < ncell - 1:
            e = mu * rf[i + 1] / dr
            up[i], dg[i] = e, dg[i] - e
        else:
            if wall_at_face:                          # 壁面にすべりなし条件
                dg[i] -= mu * rf[i + 1] / (dr / 2.0)
            else:                                     # セル中心に固定したすべりなし条件
                lo[i], dg[i], up[i], rhs[i] = 0.0, 1.0, 0.0, 0.0
    for i in range(1, ncell):                         # Thomasアルゴリズム
        m = lo[i] / dg[i - 1]
        dg[i] -= m * up[i - 1]
        rhs[i] -= m * rhs[i - 1]
    u = [0.0] * ncell
    u[-1] = rhs[-1] / dg[-1]
    for i in range(ncell - 2, -1, -1):
        u[i] = (rhs[i] - up[i] * u[i + 1]) / dg[i]
    q = sum(u[i] * math.pi * (rf[i + 1] ** 2 - rf[i] ** 2) for i in range(ncell))
    return q
 
 
def fit_slope(xs, ys):
    """log y 対 log x の最小二乗の傾き"""
    lx = [math.log(x) for x in xs]
    ly = [math.log(y) for y in ys]
    n = len(lx)
    mx, my = sum(lx) / n, sum(ly) / n
    num = sum((lx[i] - mx) * (ly[i] - my) for i in range(n))
    den = sum((lx[i] - mx) ** 2 for i in range(n))
    return num / den
 
 
_seed = 20260901
 
 
def unit_normal():
    """どこで回しても同じ数が出るよう、LCGの上に載せたBox-Muller"""
    global _seed
    out = []
    for _ in range(2):
        _seed = (1103515245 * _seed + 12345) % (2 ** 31)
        out.append((_seed + 0.5) / 2 ** 31)
    return math.sqrt(-2.0 * math.log(out[0])) * math.cos(2 * math.pi * out[1])
 
 
def measured_exponent(dmin, dmax, ntube, noise, repeat=200):
    """Poiseuilleの実験: 流量は真値、直径は相対誤差を抱えて読む"""
    slopes = []
    for _ in range(repeat):
        ds, qs = [], []
        for k in range(ntube):
            f = k / (ntube - 1)
            d_true = dmin * (dmax / dmin) ** f
            qs.append(poiseuille_q(d_true / 2.0))
            ds.append(d_true * (1.0 + noise * unit_normal()))
        slopes.append(fit_slope(ds, qs))
    mean = sum(slopes) / len(slopes)
    sd = math.sqrt(sum((s - mean) ** 2 for s in slopes) / len(slopes))
    return mean, sd
 
 
print("[1] wall half a cell off  (R = 1.000 mm)")
print(" N   Q_face/Q_exact   Q_centre/Q_exact   (1-1/2N)^4")
qex = poiseuille_q(RADIUS)
for n in (10, 20, 40, 80):
    a = solve_pipe_fv(n, True) / qex
    b = solve_pipe_fv(n, False) / qex
    print(f"{n:3d}      {a:.4f}            {b:.4f}          {(1-1/(2*n))**4:.4f}")
 
print()
print("[2] 1% error on D vs the resulting error on Q")
for e in (0.005, 0.01, 0.02):
    print(f"  dD/D = {e*100:4.1f}%  ->  dQ/Q = {((1+e)**4-1)*100:5.2f}%")
 
print()
print("[3] exponent fitted from 12 tubes, 200 repeats")
print(" range of D            noise on D    exponent (mean +- sd)")
for (lo_d, hi_d, tag) in ((0.10e-3, 0.30e-3, "0.10 - 0.30 mm"),
                          (0.015e-3, 0.60e-3, "0.015- 0.60 mm")):
    for nz in (0.0, 0.01, 0.03):
        m, s = measured_exponent(lo_d, hi_d, 12, nz)
        print(f" {tag}        {nz*100:4.1f}%       {m:.3f} +- {s:.3f}")

[1] wall half a cell off  (R = 1.000 mm)
 N   Q_face/Q_exact   Q_centre/Q_exact   (1-1/2N)^4
 10      1.0100            0.8100          0.8145
 20      1.0025            0.9025          0.9037
 40      1.0006            0.9506          0.9509
 80      1.0002            0.9752          0.9752
 
[2] 1% error on D vs the resulting error on Q
  dD/D =  0.5%  ->  dQ/Q =  2.02%
  dD/D =  1.0%  ->  dQ/Q =  4.06%
  dD/D =  2.0%  ->  dQ/Q =  8.24%
 
[3] exponent fitted from 12 tubes, 200 repeats
 range of D            noise on D    exponent (mean +- sd)
 0.10 - 0.30 mm         0.0%       4.000 +- 0.000
 0.10 - 0.30 mm         1.0%       3.996 +- 0.033
 0.10 - 0.30 mm         3.0%       3.993 +- 0.098
 0.015- 0.60 mm         0.0%       4.000 +- 0.000
 0.015- 0.60 mm         1.0%       4.000 +- 0.009
 0.015- 0.60 mm         3.0%       3.995 +- 0.028

1つ目の表の2列目は、格子を2倍にするたびに誤差が4分の1になります。2次収束です。3列目はそうなりません。誤差は半分ずつしか減らず、4列目の $(1-1/2N)^{4}$ と小数第3位まで一致します。

この一致が診断です。セル中心にすべりなし条件を固定した瞬間、計算は半径 $R - \Delta r/2$ の管を解いています。離散化誤差ではなく、別の管です。そしてその半セルが流量では4倍になります。N = 20なら半径は2.5%小さく、流量は9.75%小さくなります。

壁をノードの中間に置く問題は、有限体積だけのものではありません。格子ボルツマンのbounce-backも、壁を2つのノードの中間に立てます。どちらにせよ、「壁がどこにあるとコードが信じているか」を先に確定させる必要があります。

指数4はノイズの中でも生き残るか#

Poiseuilleの状況は私たちと逆でした。彼は $D$ を知って $Q$ を予測したのではなく、測った $Q$ から指数を読み取らねばなりませんでした。しかも直径は最も測りにくい量です。内径0.015 mmのガラス管を1%の精度で読む作業を想像してみてください。

3つ目の表がその実験です。直径にだけ相対誤差を入れ、 $\log Q$ – $\log D$ の傾きを200回引き直しました。管の直径範囲が0.10–0.30 mmのとき、1%のノイズは傾きに $\pm 0.033$ の揺れを作ります。範囲を0.015–0.60 mmまで広げると、同じノイズ、同じ本数で $\pm 0.009$ に減ります。

理由は回帰の見積もり1本で見えます。

\sigma_{n} \approx \frac{n\,\sigma_{\delta D/D}}{\sqrt{N}\;\sigma_{\ln D}}

$n = 4$ は真の指数、 $N$ は管の本数、 $\sigma_{\ln D}$ は $\log D$ 軸上で管が散らばっている幅です。狭い範囲では $4 \times 0.01 / (\sqrt{12} \times 0.331) = 0.035$ 、広い範囲では $0.010$ 。表の0.033と0.009がそれです。

指数の精度は、定規を精密にして得るものではありません。 $\log D$ 軸のてこを長く取って得るものです。Poiseuilleが髪の毛の5分の1の太さまで下りた理由がここにあります。

noise on D3.0%Dmax/Dminx3tubes12

Leave the bench running with Dmax/Dmin at 3: the histogram spreads out and the fitted exponent wanders (mean 4.000, sd 0.000 over 0 fits). Now push the range to 40 without touching the noise — the same clumsy ruler, the same number of tubes, and the histogram collapses onto 4. The exponent is bought with the lever arm in log D, not with a better ruler.

Dmax/Dmin を3に置いてヒストグラムが広がるのを見た後、ノイズはそのままで範囲だけ40まで押してみてください。定規も本数も同じなのに、分布が4に集まります。

Stokesが $K''$ を $\pi/128$ に置き換えた後#

Poiseuilleの $K''$ が $\pi/128\mu$ に変わったのは、単位の整理ではありませんでした。実験が残した空の定数の席に、Navier方程式の $\mu$ が入った出来事です。Stokesがこの法則をNavier方程式から数学的に導いて見せたことで、管で測った値と方程式に入れる物性が同じものだと確認されました。

その後、この式は向きが反転します。指数を測る実験ではなく、 $\mu$ を測る粘度計になったのです。CGSの粘度単位poiseは彼の名前から来ています。方程式はPoiseuilleの実験とStokesの補強を経てNavier–Stokes方程式として定着し、2000年にはクレイ研究所の7大難問の1つになりました。

私たちのコードでも同じ構造が繰り返されます。壁関数の定数、有効直径、接触角 — 実験が形を与え、理論が係数を埋めます。形が $D^{4}$ のように急なら、係数より先に幾何が誤差を支配します。

円管の検証で私が最初に見るもの

流量がずれたときにスキームから開けると時間を捨てます。順序はこうです。

まず格子を2倍に細かくして、誤差が4分の1になるかを見ます。減らなければ離散化の問題ではありません。次に誤差比を $(1-1/2N)^{4}$ のような幾何因子と照合します。合えば壁位置です。最後に流量の誤差を直径の誤差で割ります。その商が4に近ければ、半径1つで全部が説明できるという意味です。

Poiseuilleは直径を知り得なかったので指数を読み、私たちは指数を知っているので直径を逆算できます。同じ $D^{4}$ を両側から使っているわけです。