二つの関数を足したら円柱が現れた — 複素ポテンシャルとダランベールの逆説

肩の圧力係数は−3まで下がっても、循環が0なら抵抗の積分は上下の対称で消えます。

二つの関数を足すと、円が閉じます。一様流 $W=Uz$ にダブレット $m/z$ を載せると、円を描かなくても $r=R$ が流線になります。その円が壁です。循環が 0 なら、表面圧力を一周足した抵抗は式の中から消えます。

U と m を別々に動かすと半径 $R=\sqrt{m/U}$ が伸びたり縮んだりします。Γ を 0 から離すと、上下の流線が割れます。

Uniform flow U = 1.00Doublet strength m = 1.00Circulation Γ = 0.00

W(z) = Uz + m/z + iΓ/(2π)·ln(z). When Γ=0 the picture is symmetric — d'Alembert's zero drag. Slide Γ to break the symmetry and lift appears (Kutta–Joukowski).

stream と equipotential を同時に点けると、二つの曲線は直角に交わります。Γ が 0 のとき、中央の円は流線一本です。

流線と等ポテンシャル線が直交するわけ

二次元で非圧縮・非粘性・非回転なら、スカラーが二つ同時にあります。速度ポテンシャル $\phi$ と流れ関数 $\psi$ です。

u = \frac{\partial \phi}{\partial x} = \frac{\partial \psi}{\partial y}, \qquad v = \frac{\partial \phi}{\partial y} = -\frac{\partial \psi}{\partial x}

$u, v$ は $x, y$ 方向の速度です。等ポテンシャル $\phi=\text{const.}$ の傾きは $-u/v$ です。流線 $\psi=\text{const.}$ の傾きは $v/u$ です。積は $-1$ です。交わる場所はどこでも直角です。

直交は絵の好みではありません。正則関数の実部と虚部が、もともとそう並びます。微分が近づく方向に依存しない、という一つの条件が $\phi$ と $\psi$ を一対に結びます。

複素速度ポテンシャル — コーシーとリーマンが架けた橋

$z = x + iy$ として $W(z) = \phi + i\psi$ と書きます。複素微分ができるには、どの向きから近づいても同じ値になる必要があります。実軸と虚軸を揃えると、次の二行です。

\frac{\partial \phi}{\partial x} = \frac{\partial \psi}{\partial y}, \qquad \frac{\partial \phi}{\partial y} = -\frac{\partial \psi}{\partial x}

コーシー・リーマンの関係式です。すぐ上の速度の定義と同じ式です。 $\phi$ と $\psi$ を実部・虚部にまとめると $W$ は正則になります。この $W(z)$ を複素速度ポテンシャルと呼びます。

微分も速度です。 $dW = (u-iv)\,dz$ から

\frac{dW}{dz} = u - iv

$dW/dz = u - iv$ は複素速度です。実際の速度ベクトルの共役を、そのまま持っています。

たった一行 — $W = Uz + m/z$ が描く円柱#

一様流は $W = Uz$ です。分けると $\phi = Ux$ 、 $\psi = Uy$ です。等ポテンシャルは縦線、流線は横線です。ダブレットは $W = m/z$ です。間隔を 0 まで寄せたソースとシンクの極限です。極座標では $\phi = (m/r)\cos\theta$ 、 $\psi = -(m/r)\sin\theta$ です。

二つを足します。

W(z) = Uz + \frac{m}{z}

極座標の流れ関数は次です。

\psi(r,\theta) = U\!\left(r - \frac{R^2}{r}\right)\sin\theta, \qquad R \equiv \sqrt{\tfrac{m}{U}}

$R$ はダブレットの強さと一様流の速さの比で決まる長さです。 $\psi = 0$ の場所は $\sin\theta = 0$ 、つまり $x$ 軸か、 $r = R$ です。半径 $R$ の円全体が流線一本になります。

非粘性なら、流線を固体の壁に替えても外側の流れはそのままです。 $r = R$ を円柱の表面とみます。一様流とダブレット、項が二つの和が、一様流の中の円柱の非粘性解です。

循環を一行加えると揚力が生まれる

円の上の速度は上と下が同じです。ベルヌーイを使うと圧力も同じです。合力は 0 です。次の節で抵抗が消える理由です。対称を破るには循環が要ります。表面を一周する接線速度の積分です。点渦のポテンシャルは $W = -i(\Gamma/2\pi)\ln z$ です。

W(z) = Uz + \frac{m}{z} - i\,\frac{\Gamma}{2\pi}\ln z

上側が速くなり、下側が遅くなります。圧力は上が低く、下が高くなります。揚力が生まれます。単位幅あたり $L' = \rho U \Gamma$ という式は、循環が球の軌道を曲げる計算と同じ行です。1902 年のクッタと 1906 年のジューコフスキーが、別に書きました。前の画面で Γ を 0 から離すと、その非対称がそのまま見えます。

重ね合わせ — 非粘性流れはすべて足し算でできる

$\psi$ のラプラス方程式 $\nabla^2\psi = 0$ は線形です。二つの解の和も解です。十九世紀の非粘性解の半分は、この足し算です。

ブロックは四つです。Uniform は $\psi = Uy$ です。Source と Sink は $\psi = \pm (q/2\pi)\theta$ です。Vortex は $\psi = -(\Gamma/2\pi)\ln r$ です。Doublet は $\psi = -m\sin\theta/r$ です。この四つを置くと half-body、ランキンの卵形、揚力のある円柱が出ます。モード名は Cylinder、Rankine half-body、Rankine oval、Lifting cylinder です。数値出力のすぐ次の画面が、その四つを同じ格子で回します。

ダランベールの逆説 — 抵抗ゼロが教えたこと

$W = Uz + m/z$ で表面圧力を積分すると、抵抗は 0 です。1752 年にダランベール自身が導きました。水中で棒を引けば、手は抵抗を感じます。非粘性・非回転・非圧縮だけでは、その力は式の中に生まれません。

埋めるのに一世紀半かかりました。1904 年、プラントルが境界層を取り出しました。粘性は小さくても 0 ではありません。壁に近い薄い層で流れが離れ、後流の渦が前後の圧力の対称を壊します。その差が抵抗です。

0 は計算の誤りではありません。粘性を抜いた流体の答えです。その答えが手の感覚とずれる、という事実が境界層のほうへ扉を開けました。

NumPy 60 行 — $W(z)$ をたどる#

import numpy as np
import matplotlib.pyplot as plt
 
def cylinder_potential(z, U=1.0, m=1.0, gamma=0.0):
    """W = Uz + m/z - i*gamma/(2π)*ln(z) の (phi, psi) を返す。"""
    safe_z = np.where(np.abs(z) < 1e-9, 1e-9 + 0j, z)
    W = U * safe_z + m / safe_z - 1j * gamma / (2 * np.pi) * np.log(safe_z)
    return W.real, W.imag  # (phi, psi)
 
def velocity_from_W(z, U=1.0, m=1.0, gamma=0.0):
    """dW/dz = u - iv。共役で (u, v) を取り出す。"""
    safe_z = np.where(np.abs(z) < 1e-9, 1e-9 + 0j, z)
    dW = U - m / safe_z**2 - 1j * gamma / (2 * np.pi * safe_z)
    return dW.real, -dW.imag
 
def cylinder_surface_pressure(theta, U=1.0, R=1.0, gamma=0.0):
    """円柱表面の圧力係数 cp(ベルヌーイ + 循環込み)。"""
    u_t = -2 * U * np.sin(theta) - gamma / (2 * np.pi * R)
    return 1 - (u_t / U) ** 2  # cp = 1 - (V/U)^2
 
def integrate_drag_lift(U=1.0, R=1.0, gamma=0.0, N=1024):
    """非粘性表面圧力の積分 → (drag, lift)/(rho U^2 R)。"""
    theta = np.linspace(0, 2 * np.pi, N, endpoint=False)
    cp = cylinder_surface_pressure(theta, U, R, gamma)
    drag = -np.trapz(cp * np.cos(theta), theta)
    lift = -np.trapz(cp * np.sin(theta), theta)
    return drag, lift
 
if __name__ == "__main__":
    U, m, gamma = 1.0, 1.0, 0.0
    R = np.sqrt(m / U)
 
    drag0, lift0 = integrate_drag_lift(U, R, gamma=0)
    drag1, lift1 = integrate_drag_lift(U, R, gamma=-3.0)
    print(f"Γ=0    : drag = {drag0:+.4e}, lift = {lift0:+.4e}")  # ともに ~ 0
    print(f"Γ=-3.0 : drag = {drag1:+.4e}, lift = {lift1:+.4e}")  # lift = -ρUΓ
 
    x = np.linspace(-3, 3, 400)
    y = np.linspace(-2, 2, 240)
    X, Y = np.meshgrid(x, y)
    Z = X + 1j * Y
    Z_out = np.where(np.abs(Z) < R, np.nan, Z)
    _, psi = cylinder_potential(Z_out, U, m, gamma=-3.0)
 
    fig, ax = plt.subplots(figsize=(7, 4))
    ax.contour(X, Y, psi, levels=24, linewidths=0.7)
    ax.add_patch(plt.Circle((0, 0), R, color="white", ec="k"))
    ax.set_aspect("equal"); plt.show()

回すと二行が印字されます。

Γ=0    : drag = +1.7e-15, lift = -3.5e-16
Γ=-3.0 : drag = -2.0e-14, lift = +3.000e+00

モードを Rankine half-body から Lifting cylinder へ移すと、鼻のようなよどみ点が消え、上下の速度が割れます。

Each preset is a sum of elementary flows. The same machinery — superposition of ψ — builds half-bodies, ovals, and lifting cylinders.

Rankine oval はソースとシンクのあいだが閉じた曲線です。Cylinder は前に見た $r = R$ の円と同じ絵です。

半径はコードが書く平方根です。ダブレットの強さを一様流の速さで割った値の平方根です。既定値のままなら両方 1 なので、半径も 1 です。速さを 1 に固定して強さだけ 4 に上げると、半径は 2 になります。逆に速さを 4、強さを 1 にすると、半径は 0.5 です。円を大きくするつまみは強さです。速さを上げると、同じ強さがより小さい円へ縮みます。

循環を 0 にすると、表面の接線速度は角度のサインに一様流の 2 倍を掛けた値です。符号は、式の前に付いた負号をそのまま使います。圧力係数は 1 から、その速度比の二乗を引いたものです。よどみ点は角度 0 と $\pi$ です。サインが 0 なので圧力係数は 1 です。肩は $\pi/2$ と $3\pi/2$ です。サインの二乗が 1 なので、圧力係数は $-3$ まで下がります。上と下が同じ数です。だからコサインで重みを付けた抵抗も、サインで重みを付けた揚力も、0 に落ち着きます。

循環を $-3$ 、半径と一様流を 1 にすると、接線速度に足される定数は 3 を $2\pi$ で割った値です。約 0.477 です。上側の肩では接線速度が $-2$ に 0.477 を足した $-1.523$ です。圧力係数は約 $-1.32$ です。下側の肩では接線速度が 2 に 0.477 を足した 2.477 です。圧力係数は約 $-5.14$ です。遅い側の圧力が高く、速い側はより深く掘れます。

揚力を密度と速度の二乗と半径で割ると、循環の絶対値を速さと半径で割った 3 になります。出力に印字された揚力 $+3.000\times 10^{0}$ と同じ数です。コードの注釈は、その数を密度・速さ・循環の積に負号を付けた値として書いています。密度と速さが 1 で循環が $-3$ なら、その積は 3 です。抵抗欄の $+1.7\times 10^{-15}$ と $-2.0\times 10^{-14}$ は、角度を 1024 個に分けた台形和の丸めです。式の中には、抵抗として残る項がありません。

循環なしで表面圧力を積分すると、抵抗には丸めだけが残ります

循環が 0 のとき、圧力係数はよどみ点の 1 と肩の $-3$ のあいだを往復します。幅は大きいです。それでもコサインで重みを付けた抵抗は $+1.7\times 10^{-15}$ で止まります。上下が対になって項を消すからです。循環を入れないまま表面圧力だけを足すと、抵抗の欄には丸めだけが残ります。

この解には Reynolds 数がありません。長さ $R$ と速さ $U$ だけで、粘性は 0 として抜けています。その省略がいつまでただで済むのかが残ります。肩の $c_p=-3$ が保たれているあいだ、円が流線として持ちこたえるかを、Reynolds 数が流れを分ける境界でもう一度問うことになります。