全圧が1,080 Pa下がった場所の散逸はゼロだった — ベルヌーイ定数が流線を横切るとき

全圧が流線を横切って変わるのは損失ではなく渦度です。損失は全圧が流線に沿って変わるときにだけ生じます。

1738年、父は息子の本より6年早い刊行年を刷らせた#

Daniel Bernoulli は1738年に『Hydrodynamica』を出版しました。表紙には「Johann の息子」と記しています。父と和解したかったからです。その父 Johann はほぼ同じ内容を『Hydraulica』として別に出し、出版業者に圧力をかけて刊行年を1732年と刷らせました。息子より6年早く見えるようにしたのです。

一世紀半のち、Horace Lamb が整理をつけました。オイラー方程式を積分すればベルヌーイの定理が出る、というわけです。父子が争った二冊は、実は同じ式の二つの顔でした。ただしこの「積分すれば」には条件が一つ付いていて、CFD の後処理で全圧コンターを読み違えさせる原因は、まさにその条件にあります。

この記事ではその条件をランキン渦ひとつで測ります。先に結論を書くと、回転するコアの内側で全圧は 1,080 Pa 下がりますが、粘性散逸はちょうど0です。そして全圧が完全に平坦な外側では散逸が0ではありません。二つが逆に出ます。

下のシミュレーションでプローブをコアの内外へ動かしてみてください。

core a0.50 mOMEGA60 1/sprobe r0.25 m

u = 0.00 m/sp = 0.0 Pap0 = 0.0 Pad(p0)/dr = 0.0rho u omega = 0.0phi = 0.0000 W/m^3

Drag probe r across the dashed core edge. Inside the core the two right-hand numbers stay equal and nonzero, so p0 climbs with r while the red element rotates without changing shape. Outside, both go to zero — p0 is perfectly flat — yet the purple element keeps shearing, which is where a viscous fluid would actually lose energy.

probe r を動かしながら、右のグラフの緑の $p_0$ 曲線がどこで平坦になるかを見ます。そしてキャンバス左側の四角形が、いつ形を保ち、いつ潰れるかを一緒に見ます。この二つの観察が食い違う地点が、この記事の主題です。

オイラー方程式を積分するとベルヌーイになる — どの方向へ?

密度が一定で体積力のない定常非粘性流れの運動方程式はこうです。

(\mathbf{u}\cdot\nabla)\mathbf{u} = -\frac{1}{\rho}\nabla p

$\mathbf{u}$ は速度、 $p$ は圧力、 $\rho$ は密度です。左辺にベクトル恒等式を使います。

(\mathbf{u}\cdot\nabla)\mathbf{u} = \nabla\!\left(\frac{|\mathbf{u}|^2}{2}\right) - \mathbf{u}\times\boldsymbol{\omega}

$\boldsymbol{\omega} = \nabla\times\mathbf{u}$ は渦度です。二式を合わせ、全圧 $p_0 = p + \tfrac{1}{2}\rho|\mathbf{u}|^2$ でまとめると、残るのは一行です。

\nabla p_0 = \rho\,\mathbf{u}\times\boldsymbol{\omega}

これが Crocco 形式のオイラー方程式です(非圧縮・等エントロピー条件)。右辺が0でなければ全圧は空間で変わります。しかし、どの方向に変わるのかが肝心です。

$\mathbf{u}\times\boldsymbol{\omega}$ は $\mathbf{u}$ に垂直です。ですから両辺に $\mathbf{u}$ を内積すると右辺が消えます。

\mathbf{u}\cdot\nabla p_0 = 0

流線に沿えば $p_0$ は常に一定です。渦度があってもなくても関係ありません。これがベルヌーイの定理の正確な射程です。逆に $p_0$ があらゆる場所で同じであるには $\mathbf{u}\times\boldsymbol{\omega} = 0$ でなければならず、実質的に非回転流れを意味します。父子が争った二冊の違いはここにはなく、条件は Lamb が付けた脚注の側にありました。

ランキン渦: 回るコアと回らない外側

この条件が一つの流れの中で同時にオンとオフになる例がランキン渦です。半径 $a$ の内側は剛体のように回り、外側は自由渦です。

u_\theta(r) = \begin{cases} \Omega r, & r < a \\ \Omega a^2 / r, & r \ge a \end{cases} \qquad \omega_z(r) = \begin{cases} 2\Omega, & r < a \\ 0, & r \ge a \end{cases}

$\Omega$ はコアの角速度、 $u_\theta$ は周方向速度です。圧力は半径方向の運動量 $\mathrm{d}p/\mathrm{d}r = \rho u_\theta^2 / r$ を積分して得ます。外側では圧力が下がる分と動圧が上がる分がちょうど打ち消し合い、 $p_0$ が $p_\infty$ に固定されます。コアの内側では相殺しません。

p_0(r) - p_\infty = \rho\Omega^2 (r^2 - a^2), \qquad r < a

$r = 0$ での欠損は $\rho\Omega^2 a^2$ 、つまりコア縁の速度の二乗に $\rho$ を掛けた値です。 $\rho = 1.2$ 、 $\Omega = 60\,\mathrm{s^{-1}}$ 、 $a = 0.5\,\mathrm{m}$ なら $u_{\max} = 30$ m/s で、欠損は1,080 Pa です。

Python で両辺を同じ半径で測ってみた#

$\nabla p_0$ を中心差分で求め、 $\rho\,u_\theta\omega_z$ と突き合わせます。散逸関数 $\phi = \mu(2S_{r\theta})^2$ も一緒に測ります。

import math
 
RHO, MU = 1.2, 1.8e-5      # kg/m^3, Pa*s
OMEGA, A = 60.0, 0.5       # 1/s, m   -> u_max = 30 m/s
P_INF = 101325.0           # Pa
 
def rankine_velocity(r):
    return OMEGA * r if r < A else OMEGA * A * A / r
 
def rankine_spin(r):
    return 2.0 * OMEGA if r < A else 0.0
 
def rankine_pressure(r):
    if r >= A:
        return P_INF - 0.5 * RHO * (OMEGA * A * A / r) ** 2
    return P_INF - RHO * OMEGA**2 * A**2 + 0.5 * RHO * OMEGA**2 * r**2
 
def total_head(r):
    return rankine_pressure(r) + 0.5 * RHO * rankine_velocity(r) ** 2
 
def shear_dissipation(r):
    s = 0.0 if r < A else -OMEGA * A * A / r**2   # 2*S_rtheta
    return MU * s * s
 
def crocco_residual(r, h=1e-6):
    lhs = (total_head(r + h) - total_head(r - h)) / (2 * h)   # d(p0)/dr
    rhs = RHO * rankine_velocity(r) * rankine_spin(r)         # rho*(u x omega)_r
    return lhs, rhs
 
print("  r[m]   u[m/s]   p[Pa]      p0[Pa]     dp0/dr    rho*u*w    phi[W/m^3]")
for r in (0.10, 0.25, 0.40, 0.60, 1.00):
    lhs, rhs = crocco_residual(r)
    print("%6.2f %8.2f %10.1f %10.1f %9.1f %10.1f %11.4f"
          % (r, rankine_velocity(r), rankine_pressure(r), total_head(r),
             lhs, rhs, shear_dissipation(r)))
 
def bernoulli_gap(r1, r2):
    return total_head(r2) - total_head(r1)
 
print()
print("core  p0(0.00) - p0(%.2f) = %8.1f Pa" % (A, total_head(0.0) - total_head(A)))
print("outer p0(%.2f) - p0(1.00) = %8.1f Pa" % (A, total_head(A) - total_head(1.0)))
print("cross-streamline gap r=0.10 -> 0.40 : %8.1f Pa" % bernoulli_gap(0.10, 0.40))
print("along-streamline  gap r=0.40 -> 0.40 : %8.1f Pa" % bernoulli_gap(0.40, 0.40))
print("dissipation at r=0.25 (core)  : %.4f W/m^3" % shear_dissipation(0.25))
print("dissipation at r=0.60 (outer) : %.4f W/m^3" % shear_dissipation(0.60))

  r[m]   u[m/s]   p[Pa]      p0[Pa]     dp0/dr    rho*u*w    phi[W/m^3]
  0.10     6.00   100266.6   100288.2     864.0      864.0      0.0000
  0.25    15.00   100380.0   100515.0    2160.0     2160.0      0.0000
  0.40    24.00   100590.6   100936.2    3456.0     3456.0      0.0000
  0.60    25.00   100950.0   101325.0       0.0        0.0      0.0313
  1.00    15.00   101190.0   101325.0       0.0        0.0      0.0040
 
core  p0(0.00) - p0(0.50) =  -1080.0 Pa
outer p0(0.50) - p0(1.00) =      0.0 Pa
cross-streamline gap r=0.10 -> 0.40 :    648.0 Pa
along-streamline  gap r=0.40 -> 0.40 :      0.0 Pa
dissipation at r=0.25 (core)  : 0.0000 W/m^3
dissipation at r=0.60 (outer) : 0.0313 W/m^3

4列目と5列目がすべての半径で一致しています。Crocco の関係が小数点以下まで合っています。そしてコアの内側では $r = 0.10$ と $r = 0.40$ の間の全圧差が648 Pa です。同じ半径の二つの角度の間では0です。ベルヌーイを流線に沿って使えば合い、横切って使えば648 Pa を取りこぼします。

全圧が平坦な場所で散逸が0ではなかった#

最後の列がひっくり返っています。全圧が1,080 Pa 下がるコア内側の散逸が0で、全圧がぴたりと平坦な外側の散逸が0ではありません。

理由は、散逸が回転ではなく変形に結びついているからです。剛体回転は流体要素を回すだけで、潰しません。ひずみ速度テンソルが0なので $\phi = 0$ です。自由渦は逆です。渦度は0ですが $u_\theta \propto 1/r$ なので、内側と外側が違う速度で通り過ぎます。要素はせん断され続けます。

上のシミュレーションの四角形がこれをそのまま見せてくれます。コアの内側(赤)では正方形がそのまま回り、外側(紫)では平行四辺形へと崩れます。回転系で何が保存され何が保存されないかはコリオリ力とロタルピーで扱いましたが、ここでも軸は同じです。回ることと仕事をすることは、別の帳簿に記されます。

五つの流れを同じ表に並べる

渦度と散逸は独立です。四つの組み合わせがすべて実在します。

流れ	$\boldsymbol{\omega}$	流線に沿う $p_0$	流線を横切る $p_0$	粘性散逸
一様流	0	一定	一定	0
自由渦 ( $r > a$ )	0	一定	一定	> 0
ランキンコア ( $r < a$ )	$2\Omega$	一定	1,080 Pa 変化	0
せん断流入、非粘性	$\ne 0$	一定	変化	0
粘性後流	$\ne 0$	減少	変化	> 0

表で太字にしたマスが一つしかない点が重要です。損失と呼べるのは、流線に沿って $p_0$ が減る場合だけです。残る四行で $p_0$ が空間的に変わるのは、すべて渦度の幾何であって、エネルギーが消えたわけではありません。

全圧コンターが二つの理由で同じに見えるとき

実務でこの区別が崩れる場所は、出口断面の全圧コンターです。大気境界層や発達した管内流れを流入条件として与えると、最初のセルから $p_0$ が $y$ 方向に変わります。損失はまだ0です。後流が作った全圧欠損も同じコンターとして描かれます。こちらは本物の損失です。

下では二つのチャネルを同じカラースケールで並べて回してみます。

inlet shear S0.15wake deficit D0.32speed1.0x

A spread 0 PaA loss 0 PaB spread 0 PaB loss 0 Pa

Press match outlet spread — the two outlet profiles now cover the same range of p0, so a contour plot of the outlet cannot tell them apart. The dots can: in A every parcel keeps the color it entered with, in B the parcels that pass the body change color on the way out. Only a color change along a streamline is a loss.

match outlet spread を押すと二つのチャネルの出口 $p_0$ の範囲が揃います。出口コンターだけでは区別が不可能になる、ということです。代わりに点の色が流れる間に変わるかどうかを見ます。上のチャネルの点は色をそのまま保ち、下のチャネルの点だけが物体を通り過ぎながら色を変えます。

ですから損失係数を計算するとき、基準値を断面平均 $p_0$ に取ると、せん断流入では損失がないのに 0でない値が出ます。基準はその流線の入口 $p_0$ でなければなりません。流量加重平均を使う場合でも、入口と出口の両方を同じやり方で平均して初めて、差に損失だけが残ります。

数値的にもう一つ付きまといます。格子が粗いと回転領域で $p_0$ が人為的に平坦になります。数値拡散が渦度をならすと $\nabla p_0 = \rho\,\mathbf{u}\times\boldsymbol{\omega}$ の右辺が小さくなり、コアの欠損が実際より浅く出ます。渦コアを貫く格子を見るとき、全圧欠損の深さは渦度解像度の指標として使えます。円管の検証で直径の誤差が流量に四乗で効いていたことと同じで、ここでも誤差は目につきにくい量を通って入ってきます。

ダランベールが1752年にぶつかったのも同じ場所だった#

Johann は最後までニュートンの粘性理論を受け入れず、息子の Daniel も弟子のオイラーも同じでした。非粘性理論だけで物体まわりの流れを解くと、抗力は0になります。1752年にダランベールがこれを発表したとき、学界がパニックに陥った理由です。ダランベールが同じ時期に波動方程式でぶつかった問題と根は一つです。方程式が許す解の範囲をどこまでと見るのか、という問いです。

いま私たちが使う表現で書き直すとこうなります。非回転・非粘性の流れでは $p_0$ が全領域で一つの定数になり、そうなると物体の前後の圧力が対称になって積分が0になります。抗力を作るには、どこかで $p_0$ が流線に沿って下がらなければなりません。その場所を用意するのが粘性と、それが壁で生産する渦度です。

全圧コンターを開いたときに問うべきは、だから「どれだけ下がったか」ではありません。「流線に沿って下がったのか、横切って下がったのか」です。前者なら損失で、後者なら渦度です。