クッタ・ジューコフスキーの翼理論を基礎から徹底解説

飛行機の翼はなぜ揚力を生み出すのでしょうか。翼の断面形状はどのように設計すれば効率的に空気から力を受けられるのでしょうか。これらの問いに対する最も基本的かつ美しい回答を与えるのが、クッタ・ジューコフスキーの定理（Kutta-Joukowski theorem）です。

この定理は驚くほどシンプルな形をしています。「単位スパンあたりの揚力は、流体の密度・一様流速・翼周りの循環の積に等しい」というものです。つまり、翼がどんなに複雑な形状をしていても、揚力を決めるのは翼周りの循環（circulation）というたった一つの量なのです。

クッタ・ジューコフスキーの定理を理解すると、以下のような幅広い分野への応用が開けます。

• 航空工学: 翼型設計の出発点。NACA翼型をはじめとする近代的な翼型は、ジューコフスキー翼の発展形として位置づけられます
• 風力発電: 風車ブレードの空力設計において、各断面の揚力係数を見積もる基礎理論です
• 船舶工学: 船のスクリュープロペラや舵の断面設計にも同じ理論が適用されます
• 複素解析の応用: 等角写像（conformal mapping）という数学の手法が物理問題を鮮やかに解く好例です
• スポーツ科学: 野球の変化球やサッカーのフリーキックにおけるマグヌス効果も、循環と揚力の関係で統一的に理解できます

本記事では、循環と揚力の本質的な関係を明らかにし、クッタ・ジューコフスキーの定理を厳密に導出します。さらに、ジューコフスキー変換（等角写像）を用いて円柱から翼形状を生成する方法を解説し、Pythonによる可視化を通じて理論の内容を直感的に確認します。

本記事の内容

• クッタ・ジューコフスキーの定理の物理的意味と数学的導出
• クッタ条件（後縁条件）の物理的必然性
• ジューコフスキー変換による翼形状の生成原理
• ジューコフスキー翼の幾何学的特性とパラメータ依存性
• 薄翼理論の基礎と揚力係数の解析的結果
• 出発渦と循環の起源
• Pythonによるジューコフスキー変換の可視化・翼形状のパラメトリック変化・圧力分布・循環と揚力の関係

前提知識

この記事を読む前に、以下の記事を読んでおくと理解が深まります。

• 渦度と循環 — 循環の定義と物理的意味。本記事の中心概念です
• ポテンシャル流れ — 非粘性・非回転流れの基本理論。円柱まわりの流れを扱います
• 流れ関数と速度ポテンシャル — 複素ポテンシャルの定義と性質
• ベルヌーイの定理 — 圧力分布の計算に必要です
• 複素関数 — 等角写像の数学的基盤

揚力と循環の関係 — 直感的な理解

翼はなぜ揚力を受けるのか

翼の揚力を理解するために、まず最も単純な状況を考えましょう。一様な流れの中に翼型の物体を置いたとき、翼の上面と下面で流速が異なります。上面では流れが速く、下面では流れが遅くなります。ベルヌーイの定理によれば、流速が大きいところでは圧力が低くなるため、上面の圧力が下面の圧力よりも低くなり、正味の上向きの力、すなわち揚力が発生します。

しかし、「なぜ上面で速く、下面で遅いのか」という本質的な問いに答えるには、流速の非対称性の原因を特定する必要があります。ここで登場するのが循環（circulation）の概念です。

循環が流速の非対称性を生む

渦度と循環で学んだように、循環 $\Gamma$ は閉曲線に沿った速度の線積分として定義されます。

$$ \begin{equation} \Gamma = \oint_C \bm{v} \cdot d\bm{l} \end{equation} $$

ここで、翼の周りの閉じた経路に沿って速度を積分したときに、$\Gamma \neq 0$ であることが揚力の源です。循環がゼロでなければ、翼の上面と下面を通る流体粒子は同じ速さでは進みません。循環が正（反時計回り）であれば上面の流速は一様流よりも速くなり、下面では遅くなります。この流速差が圧力差を生み、圧力差が揚力を生むのです。

大雑把に言えば、翼の周りに「渦を巻くような速度成分」が重ね合わさっている、と考えることができます。一様流と循環を持つ流れを足し合わせると、上面では一様流と循環流が同じ方向を向くため加速し、下面では逆向きのため減速します。

この直感を数学的に定式化したものが、次に導出するクッタ・ジューコフスキーの定理です。

クッタ・ジューコフスキーの定理の導出

定理の主張

クッタ・ジューコフスキーの定理は、二次元の定常・非粘性・非圧縮流れにおいて、物体に働く単位スパンあたりの揚力 $L’$ が次の式で与えられることを主張します。

$$ \begin{equation} L’ = \rho U \Gamma \end{equation} $$

ここで $\rho$ は流体の密度、$U$ は一様流の速度、$\Gamma$ は物体周りの循環です。この定理の驚くべき点は、揚力が物体の形状の詳細に依存せず、循環 $\Gamma$ のみで決まるということです。円柱であっても翼型であっても、同じ循環を持つならば同じ揚力が発生します。

また、この定理は同時に抗力がゼロであること（ダランベールのパラドックス）も導きます。非粘性の仮定のもとでは、流れの方向の力は消えてしまうのです。現実の抗力を説明するには境界層理論による粘性効果の導入が必要になります。

導出の準備 — 円柱まわりの循環付きポテンシャル流れ

導出のために、ポテンシャル流れで扱う円柱まわりの流れを出発点にします。複素平面上で半径 $a$ の円柱まわりに、一様流速 $U$ の流れと循環 $\Gamma$ を重ね合わせた複素速度ポテンシャル $w(\zeta)$ は次のように表されます。

$$ \begin{equation} w(\zeta) = U\left(\zeta + \frac{a^2}{\zeta}\right) + \frac{i\Gamma}{2\pi}\ln\zeta \end{equation} $$

ここで $\zeta$ は複素平面上の座標です。右辺の第1項 $U\zeta$ は一様流、第2項 $Ua^2/\zeta$ は円柱による二重極（doublet）、第3項は循環（渦）を表しています。流れ関数と速度ポテンシャルで学んだように、複素ポテンシャルの微分が複素速度を与えます。

ブラジウスの定理

物体に作用する力を計算するために、ブラジウスの定理（Blasius theorem）を用います。これは、二次元の定常ポテンシャル流れにおいて、物体に作用する力の $x$ 成分 $F_x$ と $y$ 成分 $F_y$ が、物体を囲む任意の閉曲線 $C$ 上の積分で与えられるというものです。

$$ \begin{equation} F_x – iF_y = \frac{i\rho}{2}\oint_C \left(\frac{dw}{d\zeta}\right)^2 d\zeta \end{equation} $$

この定理は、ベルヌーイの定理を使って圧力を速度で表し、物体表面上の圧力を積分することで導かれます。ベルヌーイの定理により、定常流の圧力は $p = p_\infty + \frac{1}{2}\rho(U^2 – |v|^2)$ と書けるため、表面上の圧力積分を複素速度ポテンシャルの言葉で整理するとブラジウスの定理が得られます。

複素速度の計算

複素速度ポテンシャルを微分して複素速度を求めます。

$$ \begin{equation} \frac{dw}{d\zeta} = U\left(1 – \frac{a^2}{\zeta^2}\right) + \frac{i\Gamma}{2\pi\zeta} \end{equation} $$

ブラジウスの定理に代入するために、この式を2乗します。

$$ \left(\frac{dw}{d\zeta}\right)^2 = U^2\left(1 – \frac{a^2}{\zeta^2}\right)^2 + 2U\left(1 – \frac{a^2}{\zeta^2}\right)\frac{i\Gamma}{2\pi\zeta} + \left(\frac{i\Gamma}{2\pi\zeta}\right)^2 $$

各項を展開すると次のようになります。

$$ \left(\frac{dw}{d\zeta}\right)^2 = U^2 – \frac{2U^2 a^2}{\zeta^2} + \frac{U^2 a^4}{\zeta^4} + \frac{iU\Gamma}{\pi\zeta} – \frac{iU\Gamma a^2}{\pi\zeta^3} – \frac{\Gamma^2}{4\pi^2\zeta^2} $$

留数定理の適用

ブラジウスの定理の積分 $\oint_C \left(\frac{dw}{d\zeta}\right)^2 d\zeta$ を評価するには、複素関数で学ぶ留数定理（residue theorem）を適用します。留数定理によれば、閉曲線 $C$ に沿った複素積分は、$C$ の内部にある特異点の留数（$\zeta^{-1}$ の係数）の和に $2\pi i$ を掛けたものに等しくなります。

$$ \begin{equation} \oint_C f(\zeta)\,d\zeta = 2\pi i \sum_k \text{Res}(f, \zeta_k) \end{equation} $$

先ほどの展開を見ると、$\zeta^{-1}$ の項は $\frac{iU\Gamma}{\pi\zeta}$ の1項だけです。他の項は $\zeta^0$（定数）、$\zeta^{-2}$、$\zeta^{-3}$、$\zeta^{-4}$ のべきであり、留数定理による積分ではこれらはすべてゼロになります。

したがって、$\zeta = 0$ における留数は $\frac{iU\Gamma}{\pi}$ であり、積分は次のようになります。

$$ \oint_C \left(\frac{dw}{d\zeta}\right)^2 d\zeta = 2\pi i \cdot \frac{iU\Gamma}{\pi} = -2U\Gamma $$

揚力とダランベールのパラドックス

この結果をブラジウスの定理に代入します。

$$ F_x – iF_y = \frac{i\rho}{2}(-2U\Gamma) = -i\rho U\Gamma $$

$x$ 方向（流れ方向）が抗力、$y$ 方向（流れに垂直）が揚力ですので、実部と虚部を比較すると次の結果が得られます。

$$ \begin{equation} F_x = 0, \quad F_y = \rho U \Gamma \end{equation} $$

すなわち、抗力はゼロであり、揚力は $L’ = \rho U\Gamma$ です。これがクッタ・ジューコフスキーの定理です。

$F_x = 0$ という結果はダランベールのパラドックス（d’Alembert’s paradox）と呼ばれます。非粘性流体の理論では抗力が発生しないのです。実際の翼には摩擦抗力や圧力抗力（形状抗力）が存在しますが、これらは境界層理論で説明される粘性効果に起因します。抗力と揚力の記事でも、揚力と抗力のメカニズムの違いを詳しく解説しています。

この定理には注目すべき点がもう一つあります。$L’ = \rho U\Gamma$ という式の中に、物体の形状に関する情報が一切含まれていないことです。円柱であろうと翼型であろうと、同じ循環 $\Gamma$ を持つ流れは同じ揚力を生みます。形状の役割は「どれだけの循環が生じるかを決める」ことにあり、揚力そのものは循環だけで決まるのです。

では、翼周りの循環の大きさはどのように決まるのでしょうか。その答えを与えるのがクッタ条件です。

クッタ条件 — 循環はなぜ一意に決まるのか

ポテンシャル流れの非一意性

ここまでの議論で、循環 $\Gamma$ は自由なパラメータとして扱ってきました。数学的には、ポテンシャル流れの解は任意の循環に対して存在します。円柱のような対称物体の場合、$\Gamma = 0$ から $\Gamma = \pm\infty$ まで、どの値でも流れの支配方程式（ラプラス方程式）を満たす正当な解です。

しかし、現実の流れでは循環の値は一意に決まります。この物理的な条件がなければ、理論は揚力を予測できません。必要なのは、数学に追加の物理条件を課すことです。

クッタ条件の内容

クッタ条件（Kutta condition）は、翼の後縁（trailing edge）で流れが滑らかに離脱するという条件です。より具体的には、後縁が鋭い尖端である翼型に対して、後縁における流速が有限であることを要求します。

$$ \begin{equation} |\bm{v}|_{\text{trailing edge}} < \infty \end{equation} $$

後縁がゼロ角度の尖端（cusp型）の場合、クッタ条件は「上面と下面の流れが後縁で同じ速度で合流する」ことを意味します。有限角度の後縁の場合は、「後縁に停留点が来る」すなわち後縁での速度がゼロになることを要求します。

物理的な意味

なぜクッタ条件が物理的に正しいのでしょうか。その理由は粘性の存在にあります。

もしクッタ条件が満たされなければ、流れは後縁を回り込んで下面から上面へ急激に方向転換しなければなりません。尖った後縁の周りで流体が回り込むとき、速度は後縁で無限大に発散します。しかし、実際の流体には必ず粘性があり、このような無限速度は実現できません。粘性が働いて流れの剥離（separation）が起こり、結果として後縁から滑らかに流れが離脱する形に落ち着くのです。

つまり、クッタ条件は「非粘性理論の中に、粘性の効果を間接的に組み込むための物理的な追加条件」と理解できます。これは非常に巧妙な工夫であり、完全に非粘性の方程式を解きながら、現実的な揚力を予測できる理由です。

循環の一意な決定

クッタ条件により、翼型の形状と迎え角（angle of attack）が与えられれば、循環 $\Gamma$ の値が一意に決まります。後縁で流速が有限となるような循環はただ一つしか存在しないからです。

例えば、薄い翼型の場合、クッタ条件から決まる循環は迎え角 $\alpha$ に比例します。

$$ \begin{equation} \Gamma = \pi c U \sin\alpha \approx \pi c U \alpha \quad (\alpha \ll 1) \end{equation} $$

ここで $c$ は翼弦長（chord length）です。この結果とクッタ・ジューコフスキーの定理を組み合わせると、揚力係数 $C_L$ が $\alpha$ に対して線形に増加するという、実験とよく一致する結果が得られます。

この循環がどのようにして物理的に「生まれる」のかという問題は、出発渦の議論で明らかになります。これについてはのちほど詳しく扱います。

まず、クッタ・ジューコフスキーの定理を翼型設計に結びつける重要な道具である、ジューコフスキー変換について見ていきましょう。

ジューコフスキー変換 — 等角写像による翼型の生成

等角写像とは

ジューコフスキー変換の背景にある数学的な道具は等角写像（conformal mapping）です。等角写像とは、複素関数による写像のうち、局所的に角度を保存するものを指します。つまり、元の平面で直交する2本の曲線は、写像後の平面でも直交します。

等角写像が流体力学で威力を発揮する理由は、ラプラス方程式の解が等角写像で保存されるためです。もしある領域でラプラス方程式を解いてポテンシャル流れの解を得たならば、等角写像でその領域を別の形に変換しても、変換後の解は変換後の領域におけるラプラス方程式の解になっています。

これは実用的に極めて強力です。円柱まわりのポテンシャル流れは解析的に完全に解けていますから、円柱を翼型に写す等角写像さえ見つかれば、翼型まわりの流れが自動的に求まるのです。

ジューコフスキー変換の定義

ジューコフスキー変換（Joukowski transformation）は、$\zeta$ 平面（円柱が存在する平面）から $z$ 平面（翼型が存在する平面）への等角写像で、次の式で定義されます。

$$ \begin{equation} z = \zeta + \frac{c_0^2}{\zeta} \end{equation} $$

ここで $c_0$ は正の実数定数です。この変換は、$\zeta$ が大きい（物体から遠い）ところでは $z \approx \zeta$ となるため、無限遠方での流れの性質を保存します。

変換の特異点は $\frac{dz}{d\zeta} = 0$ となる点、すなわち次を満たす点です。

$$ \frac{dz}{d\zeta} = 1 – \frac{c_0^2}{\zeta^2} = 0 \quad \Longrightarrow \quad \zeta = \pm c_0 $$

$\zeta = \pm c_0$ の2点がこの変換の臨界点であり、等角性が破れる点です。翼の前縁と後縁はこの臨界点に対応します。

円柱から翼型を生成する仕組み

ジューコフスキー変換で翼型を生成するには、$\zeta$ 平面での円の中心と半径を適切に選ぶ必要があります。具体的には、以下のように設定します。

1. $\zeta$ 平面に中心が $\zeta_0 = -\epsilon + i\delta$ にある半径 $a$ の円を考えます。ここで $\epsilon > 0$ は小さな実数、$\delta$ も小さな実数です
2. 円が $\zeta = c_0$（ジューコフスキー変換の右側の臨界点）を必ず通るように半径を設定します。すなわち $a = |\zeta_0 – c_0| = \sqrt{(c_0 + \epsilon)^2 + \delta^2}$ とします
3. この円をジューコフスキー変換で $z$ 平面に写すと、翼型が得られます

各パラメータの役割を整理しましょう。

• $\epsilon$（実軸方向のずれ）: 翼の厚さを制御します。$\epsilon = 0$ なら厚さゼロの平板翼になります。$\epsilon$ が大きいほど翼は厚くなります
• $\delta$（虚軸方向のずれ）: 翼のキャンバー（反り）を制御します。$\delta = 0$ なら対称翼になります。$\delta > 0$ で上に反った翼型が得られます
• $c_0$: 翼弦長のスケールを決めます。ジューコフスキー翼の翼弦長はおよそ $4c_0$ になります

臨界点 $\zeta = c_0$ を円が通ることにより、この点でジューコフスキー変換の微分がゼロになり、$z$ 平面では後縁が尖った形状（尖点、cusp）になります。これがクッタ条件を適用するための幾何学的な基礎です。

もう一方の臨界点 $\zeta = -c_0$ は円の内部にあるため、$z$ 平面上の前縁は丸みを帯びた形状になります。この前縁の丸みは、翼の空力特性にとって重要な意味を持ちます。

では、実際にジューコフスキー変換がどのように円柱を翼型に写すのかを、Pythonで可視化してみましょう。

Pythonで可視化するジューコフスキー変換

変換の基本的な可視化

まず、$\zeta$ 平面上の円がジューコフスキー変換によって $z$ 平面上でどのような形状に写されるかを可視化します。同時に、変換前後の格子線の変形も描くことで、等角写像の性質を視覚的に確認します。

import numpy as np
import matplotlib.pyplot as plt

# --- パラメータ設定 ---
c0 = 1.0          # ジューコフスキー変換の定数
epsilon = 0.1      # 厚さパラメータ（実軸方向のずれ）
delta = 0.15       # キャンバーパラメータ（虚軸方向のずれ）

# 円の中心と半径
zeta_0 = -epsilon + 1j * delta
a = abs(zeta_0 - c0)  # 臨界点 c0 を通る半径

# --- ζ平面上の円を描く ---
theta = np.linspace(0, 2 * np.pi, 400)
zeta_circle = zeta_0 + a * np.exp(1j * theta)

# --- ジューコフスキー変換 ---
def joukowski(zeta, c0):
    return zeta + c0**2 / zeta

z_airfoil = joukowski(zeta_circle, c0)

# --- 格子線の生成（ζ平面の同心円と放射線） ---
r_vals = np.linspace(a * 1.05, 3.5, 8)
theta_grid = np.linspace(0, 2 * np.pi, 300)
theta_vals = np.linspace(0, 2 * np.pi, 24, endpoint=False)
r_grid = np.linspace(a * 1.02, 3.5, 200)

# --- 可視化 ---
fig, axes = plt.subplots(1, 2, figsize=(14, 6))

# ζ平面（変換前）
ax = axes[0]
ax.set_title(r"$\zeta$-plane (circle)", fontsize=14)
ax.plot(zeta_circle.real, zeta_circle.imag, "b-", linewidth=2, label="Circle")
ax.plot(zeta_0.real, zeta_0.imag, "r+", markersize=10, markeredgewidth=2)
ax.plot(c0, 0, "ko", markersize=6, label=r"$\zeta = c_0$")
ax.plot(-c0, 0, "ks", markersize=6, label=r"$\zeta = -c_0$")
# 同心円
for r in r_vals:
    zeta_g = zeta_0 + r * np.exp(1j * theta_grid)
    ax.plot(zeta_g.real, zeta_g.imag, "gray", linewidth=0.4, alpha=0.6)
# 放射線
for th in theta_vals:
    zeta_g = zeta_0 + r_grid * np.exp(1j * th)
    ax.plot(zeta_g.real, zeta_g.imag, "gray", linewidth=0.4, alpha=0.6)
ax.set_xlim(-4, 4)
ax.set_ylim(-3.5, 3.5)
ax.set_aspect("equal")
ax.legend(fontsize=9)
ax.grid(True, alpha=0.3)

# z平面（変換後）
ax = axes[1]
ax.set_title(r"$z$-plane (Joukowski airfoil)", fontsize=14)
ax.plot(z_airfoil.real, z_airfoil.imag, "b-", linewidth=2, label="Airfoil")
ax.plot(joukowski(c0 + 0j, c0).real, 0, "ko", markersize=6, label="Trailing edge")
# 同心円 → 変換後
for r in r_vals:
    zeta_g = zeta_0 + r * np.exp(1j * theta_grid)
    z_g = joukowski(zeta_g, c0)
    ax.plot(z_g.real, z_g.imag, "gray", linewidth=0.4, alpha=0.6)
# 放射線 → 変換後
for th in theta_vals:
    zeta_g = zeta_0 + r_grid * np.exp(1j * th)
    z_g = joukowski(zeta_g, c0)
    ax.plot(z_g.real, z_g.imag, "gray", linewidth=0.4, alpha=0.6)
ax.set_xlim(-4, 4)
ax.set_ylim(-3.5, 3.5)
ax.set_aspect("equal")
ax.legend(fontsize=9)
ax.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig("joukowski_transform.png", dpi=150, bbox_inches="tight")
plt.show()

上のグラフから、ジューコフスキー変換の2つの重要な特徴が読み取れます。

1. $\zeta$ 平面の円が $z$ 平面の翼型に写されている。前縁（左側）は丸みを帯びており、後縁（右側）は鋭い尖端になっています。この尖った後縁こそ、臨界点 $\zeta = c_0$ を円が通ることの帰結であり、クッタ条件を適用する幾何学的な根拠です

2. 格子線の直交性が保存されている。$\zeta$ 平面で直交していた同心円と放射線は、$z$ 平面に写された後も直交しています。これは等角写像の定義そのものであり、ラプラス方程式の解が保存されることの視覚的な確認になっています。翼型の近くでは格子が強く歪みますが、無限遠方に向かうにつれて歪みは小さくなり、一様流に漸近していきます

次に、パラメータ $\epsilon$ と $\delta$ を変化させると翼型がどのように変わるかを系統的に見てみましょう。

ジューコフスキー翼の形状とパラメトリック変化

パラメータと翼形状の関係

ジューコフスキー翼の形状は、$\epsilon$（厚さパラメータ）と $\delta$（キャンバーパラメータ）の2つの値で完全に決まります。この2つのパラメータを独立に変化させることで、薄い対称翼から厚いキャンバー翼まで、さまざまな翼型を生成できます。

直感的には次のように理解できます。$\epsilon$ を大きくすると、$\zeta$ 平面の円の中心が原点から左にずれるため、円が臨界点 $\zeta = -c_0$（前縁に対応）からより遠ざかり、前縁周辺での変形が大きくなって翼が厚くなります。$\delta$ を大きくすると、円の中心が虚軸方向にずれるため、変換後の翼型に上下の非対称性（キャンバー）が生まれます。

以下のPythonコードで、パラメータの変化による翼形状の変化を可視化します。

import numpy as np
import matplotlib.pyplot as plt

c0 = 1.0

def joukowski(zeta, c0):
    return zeta + c0**2 / zeta

def make_airfoil(epsilon, delta, c0, n=500):
    """ジューコフスキー翼の座標を返す"""
    zeta_0 = -epsilon + 1j * delta
    a = abs(zeta_0 - c0)
    theta = np.linspace(0, 2 * np.pi, n)
    zeta_circle = zeta_0 + a * np.exp(1j * theta)
    z = joukowski(zeta_circle, c0)
    return z

fig, axes = plt.subplots(2, 2, figsize=(12, 10))

# (a) 厚さパラメータ epsilon の変化（delta=0: 対称翼）
ax = axes[0, 0]
ax.set_title(r"(a) Varying thickness $\epsilon$ ($\delta=0$)", fontsize=12)
for eps in [0.0, 0.05, 0.10, 0.15, 0.20]:
    z = make_airfoil(eps, 0.0, c0)
    ax.plot(z.real, z.imag, label=rf"$\epsilon={eps:.2f}$")
ax.set_xlim(-2.5, 2.5)
ax.set_ylim(-1.0, 1.0)
ax.set_aspect("equal")
ax.legend(fontsize=9)
ax.grid(True, alpha=0.3)

# (b) キャンバーパラメータ delta の変化（epsilon=0.1）
ax = axes[0, 1]
ax.set_title(r"(b) Varying camber $\delta$ ($\epsilon=0.1$)", fontsize=12)
for d in [0.0, 0.05, 0.10, 0.15, 0.20]:
    z = make_airfoil(0.1, d, c0)
    ax.plot(z.real, z.imag, label=rf"$\delta={d:.2f}$")
ax.set_xlim(-2.5, 2.5)
ax.set_ylim(-1.0, 1.0)
ax.set_aspect("equal")
ax.legend(fontsize=9)
ax.grid(True, alpha=0.3)

# (c) 厚い翼型
ax = axes[1, 0]
ax.set_title(r"(c) Thick airfoil: $\epsilon=0.20$, $\delta=0.10$", fontsize=12)
z = make_airfoil(0.20, 0.10, c0)
ax.fill(z.real, z.imag, alpha=0.3, color="skyblue")
ax.plot(z.real, z.imag, "b-", linewidth=2)
ax.set_xlim(-2.5, 2.5)
ax.set_ylim(-1.0, 1.0)
ax.set_aspect("equal")
ax.grid(True, alpha=0.3)

# (d) 薄い翼型
ax = axes[1, 1]
ax.set_title(r"(d) Thin airfoil: $\epsilon=0.02$, $\delta=0.05$", fontsize=12)
z = make_airfoil(0.02, 0.05, c0)
ax.fill(z.real, z.imag, alpha=0.3, color="lightyellow")
ax.plot(z.real, z.imag, "b-", linewidth=2)
ax.set_xlim(-2.5, 2.5)
ax.set_ylim(-1.0, 1.0)
ax.set_aspect("equal")
ax.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig("joukowski_parametric.png", dpi=150, bbox_inches="tight")
plt.show()

このパラメトリック可視化から、以下の特徴が確認できます。

1. 図(a): $\epsilon = 0$ のとき翼は厚さゼロの平板（直線）になり、$\epsilon$ を増やすにつれて翼の厚さが系統的に増加します。$\delta = 0$ のため、すべて上下対称な翼型です。対称翼はゼロ迎え角でゼロ揚力を生じるため、正の迎え角で飛行しなければなりません

2. 図(b): $\delta$ を増やすと翼に上向きの反り（キャンバー）が加わります。キャンバーがある翼はゼロ迎え角でも正の揚力を生じるため、実用的な翼型のほとんどはキャンバーを持っています。キャンバーの増加は失速迎え角を下げる傾向があるため、設計上のトレードオフが存在します

3. 図(c)(d): 厚い翼型は構造的に強く燃料タンクの容積も確保できる一方、抗力が増加します。薄い翼型は抗力が小さいが構造的に弱くなります。航空機の翼設計では、この相反する要求のバランスを取る必要があります

ジューコフスキー翼の幾何学的特性

ジューコフスキー翼には、パラメータから読み取れるいくつかの定量的な特性があります。

翼弦長: 前縁と後縁の距離です。$\epsilon$ と $\delta$ がともに小さい場合、翼弦長は近似的に次のように表されます。

$$ \begin{equation} c \approx 4c_0\left(1 + \epsilon\right) \end{equation} $$

$\epsilon = 0$ のとき翼弦長はちょうど $4c_0$ になります。

最大厚さ: 翼の上面と下面の間の最大距離です。近似的に次の式で表されます。

$$ \begin{equation} t_{\max} \approx 3\sqrt{3}\,c_0\,\epsilon \end{equation} $$

最大厚さは翼弦の約 $25\%$ の位置に現れます。

最大キャンバー: 翼弦線（前縁と後縁を結ぶ直線）と翼の中心線（上面と下面の中点を結ぶ線）の間の最大距離です。$\delta$ に比例します。

これらの定量的な関係は、航空工学の初期において翼型設計の指針として重要な役割を果たしました。現代のNACA翼型やスーパークリティカル翼型はジューコフスキー翼の発展形ですが、基本的な設計パラメータの考え方はここに起源があります。

次に、ジューコフスキー翼まわりの流れ場を計算し、圧力分布を求めてみましょう。

翼周りの流れと圧力分布

翼周りの複素速度

$\zeta$ 平面での循環付き円柱まわりの複素速度ポテンシャルは、先ほど示したように次の式です。

$$ w(\zeta) = U\left(\zeta + \frac{a^2}{\zeta}\right) + \frac{i\Gamma}{2\pi}\ln\zeta $$

ここで、一様流速 $U$ は迎え角 $\alpha$ を持つ場合、$\zeta$ 平面では $U e^{-i\alpha}$ の方向に流れます。円の中心が $\zeta_0$ にある場合の正確な複素速度ポテンシャルは次のようになります。

$$ \begin{equation} w(\zeta) = U e^{-i\alpha}(\zeta – \zeta_0) + \frac{U a^2 e^{i\alpha}}{\zeta – \zeta_0} + \frac{i\Gamma}{2\pi}\ln(\zeta – \zeta_0) \end{equation} $$

$\zeta$ 平面での複素速度は、$w$ を $\zeta$ で微分して得られます。

$$ \begin{equation} \frac{dw}{d\zeta} = U e^{-i\alpha} – \frac{U a^2 e^{i\alpha}}{(\zeta – \zeta_0)^2} + \frac{i\Gamma}{2\pi(\zeta – \zeta_0)} \end{equation} $$

クッタ条件による循環の決定

$z$ 平面での複素速度は、連鎖律により次のように計算されます。

$$ \begin{equation} \frac{dw}{dz} = \frac{dw/d\zeta}{dz/d\zeta} \end{equation} $$

ここで $\frac{dz}{d\zeta} = 1 – \frac{c_0^2}{\zeta^2}$ です。

後縁は $\zeta = c_0$ に対応し、この点でジューコフスキー変換の微分 $\frac{dz}{d\zeta}$ がゼロになります。したがって、$z$ 平面での速度 $\frac{dw}{dz}$ が有限であるためには、分子 $\frac{dw}{d\zeta}$ も $\zeta = c_0$ でゼロでなければなりません。これがクッタ条件の数学的な表現です。

$\frac{dw}{d\zeta}\bigg|_{\zeta = c_0} = 0$ を要求すると、循環 $\Gamma$ は次のように決定されます。

$$ Ue^{-i\alpha} – \frac{Ua^2 e^{i\alpha}}{(c_0 – \zeta_0)^2} + \frac{i\Gamma}{2\pi(c_0 – \zeta_0)} = 0 $$

$c_0 – \zeta_0 = c_0 + \epsilon – i\delta = a e^{-i\beta}$ と書くと（$\beta$ は $\zeta_0$ から $c_0$ への方向角）、整理すると次の式が得られます。

$$ \begin{equation} \Gamma = 4\pi a U \sin(\alpha + \beta) \end{equation} $$

ここで $\beta = \arctan\left(\frac{\delta}{c_0 + \epsilon}\right)$ です。$\beta$ は翼のキャンバーに対応する角度であり、$\delta = 0$（対称翼）のときは $\beta = 0$ になります。

この結果から、揚力係数 $C_L$ は次のように計算されます。クッタ・ジューコフスキーの定理 $L’ = \rho U \Gamma$ と揚力係数の定義 $L’ = \frac{1}{2}\rho U^2 c \cdot C_L$（$c$ は翼弦長）を組み合わせると次のようになります。

$$ \begin{equation} C_L = \frac{2\Gamma}{Uc} = \frac{8\pi a}{c}\sin(\alpha + \beta) \approx 2\pi(\alpha + \beta) \end{equation} $$

最後の近似は $\alpha + \beta$ が小さい場合に成り立ち、$a \approx c/4$ の近似を使っています。この $C_L \approx 2\pi\alpha$ という結果は、薄翼理論の古典的な結果と一致しています。

圧力分布の計算

ベルヌーイの定理により、翼面上の圧力係数 $C_p$ は次の式で与えられます。

$$ \begin{equation} C_p = 1 – \left|\frac{v}{U}\right|^2 \end{equation} $$

ここで $v$ は翼面上での流速です。$z$ 平面の速度は $\frac{dw}{dz} = \frac{dw/d\zeta}{dz/d\zeta}$ から計算できるため、$\zeta$ 平面の円上の各点で速度を求め、ジューコフスキー変換の微分で割ることで翼面上の速度が得られます。

以下のPythonコードで、ジューコフスキー翼まわりの圧力分布を可視化します。

import numpy as np
import matplotlib.pyplot as plt

# --- パラメータ ---
c0 = 1.0
epsilon = 0.10
delta = 0.08
alpha_deg = 5.0
alpha = np.radians(alpha_deg)
U = 1.0

# 円の中心と半径
zeta_0 = -epsilon + 1j * delta
a = abs(zeta_0 - c0)

# --- ζ平面の円上の点 ---
n = 1000
theta = np.linspace(0, 2 * np.pi, n, endpoint=False)
zeta = zeta_0 + a * np.exp(1j * theta)

# --- クッタ条件から循環を計算 ---
beta = np.arctan2(delta, c0 + epsilon)
Gamma = 4 * np.pi * a * U * np.sin(alpha + beta)

# --- ζ平面の複素速度 ---
dw_dzeta = (U * np.exp(-1j * alpha)
            - U * a**2 * np.exp(1j * alpha) / (zeta - zeta_0)**2
            + 1j * Gamma / (2 * np.pi * (zeta - zeta_0)))

# --- ジューコフスキー変換の微分 ---
dz_dzeta = 1.0 - c0**2 / zeta**2

# --- z平面の速度 ---
dw_dz = dw_dzeta / dz_dzeta

# --- 圧力係数 ---
speed = np.abs(dw_dz)
Cp = 1.0 - (speed / U)**2

# --- z平面の翼型座標 ---
z = zeta + c0**2 / zeta
x = z.real
y = z.imag

# --- 上面と下面の分離 ---
# 後縁に最も近い点を基準にして分割
trailing_idx = np.argmin(np.abs(zeta - c0))
# θ方向で上面・下面を判定
upper = (theta > theta[trailing_idx]) & (theta < theta[trailing_idx] + np.pi)
lower = ~upper

fig, axes = plt.subplots(2, 1, figsize=(10, 9))

# 翼型の形状
ax = axes[0]
ax.fill(x, y, alpha=0.2, color="skyblue")
ax.plot(x, y, "b-", linewidth=1.5)
ax.set_title(rf"Joukowski airfoil ($\epsilon={epsilon}$, $\delta={delta}$, "
             rf"$\alpha={alpha_deg}°$)", fontsize=13)
ax.set_xlabel("x", fontsize=12)
ax.set_ylabel("y", fontsize=12)
ax.set_aspect("equal")
ax.set_xlim(-2.5, 2.5)
ax.set_ylim(-1.2, 1.2)
ax.grid(True, alpha=0.3)

# 圧力係数分布
ax = axes[1]
ax.plot(x[upper], Cp[upper], "r-", linewidth=1.5, label="Upper surface")
ax.plot(x[lower], Cp[lower], "b-", linewidth=1.5, label="Lower surface")
ax.invert_yaxis()  # 慣例: Cpは下向きが正
ax.set_title(r"Pressure coefficient $C_p$ distribution", fontsize=13)
ax.set_xlabel("x", fontsize=12)
ax.set_ylabel(r"$C_p$", fontsize=12)
ax.legend(fontsize=11)
ax.grid(True, alpha=0.3)
ax.axhline(0, color="k", linewidth=0.5)

plt.tight_layout()
plt.savefig("joukowski_pressure.png", dpi=150, bbox_inches="tight")
plt.show()

print(f"循環 Γ = {Gamma:.4f}")
print(f"揚力係数 C_L ≈ {2 * Gamma / (U * 4 * c0):.4f}")

圧力分布のグラフから、次の重要な特徴が読み取れます。

1. 上面の $C_p$ は下面よりも低い（グラフ上では上側）。これは上面の流速が下面よりも大きいことを意味し、この圧力差が揚力を生み出しています。$C_p$ 分布で上面の曲線が下面の曲線を「囲む」面積が揚力に比例します

2. 前縁付近で $C_p$ の値が大きく変動している。前縁の停留点付近では $C_p = 1$（速度ゼロ）に達し、そこから上面では急激に加速して $C_p$ が大きく負になります。この急加速は実際の翼でも前縁近くの吸い込みピーク（suction peak）として観測されます

3. 後縁に向かって上面と下面の $C_p$ が収束している。これはクッタ条件の帰結です。後縁で流れが滑らかに合流するため、上面と下面の圧力が等しくなります

この圧力分布の計算は、等角写像とクッタ・ジューコフスキーの定理が「円柱まわりの解ける問題」から「翼型まわりの実用的な圧力分布」を導く手順を如実に示しています。

ここまでは翼型の形状と圧力分布を個別に見てきました。次に、循環 $\Gamma$ と揚力の定量的な関係をさまざまな迎え角で確認しましょう。

循環と揚力の関係 — 定量的な検証

揚力係数の迎え角依存性

クッタ・ジューコフスキーの定理と、クッタ条件から決まる循環の式を組み合わせると、揚力係数 $C_L$ は迎え角 $\alpha$ の線形関数になることが予測されます。

$$ \begin{equation} C_L = 2\pi(\alpha + \beta) + C_{L,\text{correction}} \end{equation} $$

ここで $\beta$ はキャンバーに由来するゼロ揚力迎え角のずれ、$C_{L,\text{correction}}$ は翼の有限厚さによる小さな補正です。薄翼近似では揚力傾斜 $\frac{dC_L}{d\alpha} = 2\pi$ [rad$^{-1}$] $\approx 0.1097$ [deg$^{-1}$] という普遍的な値が得られます。

以下のPythonコードで、迎え角を変化させたときの循環と揚力係数を計算し、薄翼理論の予測と比較します。

import numpy as np
import matplotlib.pyplot as plt

# --- パラメータ ---
c0 = 1.0
epsilon = 0.10
delta = 0.08
U = 1.0
rho = 1.225  # 空気の密度 [kg/m^3]（参考値）

zeta_0 = -epsilon + 1j * delta
a = abs(zeta_0 - c0)
beta = np.arctan2(delta, c0 + epsilon)
chord = 4 * c0 * (1 + epsilon)  # 近似的な翼弦長

# --- 迎え角の範囲 ---
alpha_deg = np.linspace(-10, 15, 100)
alpha_rad = np.radians(alpha_deg)

# --- クッタ条件から循環を計算 ---
Gamma = 4 * np.pi * a * U * np.sin(alpha_rad + beta)

# --- 揚力（クッタ・ジューコフスキーの定理） ---
L_prime = rho * U * Gamma  # 単位スパンあたり [N/m]

# --- 揚力係数 ---
CL_exact = 2 * Gamma / (U * chord)

# --- 薄翼理論の予測 ---
CL_thin = 2 * np.pi * (alpha_rad + beta)

# --- ゼロ揚力迎え角 ---
alpha_L0_deg = -np.degrees(beta)

fig, axes = plt.subplots(1, 2, figsize=(13, 5.5))

# 揚力係数 vs 迎え角
ax = axes[0]
ax.plot(alpha_deg, CL_exact, "b-", linewidth=2, label="Joukowski (exact)")
ax.plot(alpha_deg, CL_thin, "r--", linewidth=1.5, label=r"Thin airfoil: $2\pi(\alpha+\beta)$")
ax.axhline(0, color="k", linewidth=0.5)
ax.axvline(alpha_L0_deg, color="gray", linewidth=0.8, linestyle=":",
           label=rf"$\alpha_{{L=0}} = {alpha_L0_deg:.2f}°$")
ax.set_xlabel(r"Angle of attack $\alpha$ [deg]", fontsize=12)
ax.set_ylabel(r"Lift coefficient $C_L$", fontsize=12)
ax.set_title(r"$C_L$ vs $\alpha$", fontsize=13)
ax.legend(fontsize=10)
ax.grid(True, alpha=0.3)

# 循環 vs 迎え角
ax = axes[1]
ax.plot(alpha_deg, Gamma, "b-", linewidth=2)
ax.axhline(0, color="k", linewidth=0.5)
ax.axvline(alpha_L0_deg, color="gray", linewidth=0.8, linestyle=":")
ax.set_xlabel(r"Angle of attack $\alpha$ [deg]", fontsize=12)
ax.set_ylabel(r"Circulation $\Gamma$ [m$^2$/s]", fontsize=12)
ax.set_title(r"$\Gamma$ vs $\alpha$", fontsize=13)
ax.grid(True, alpha=0.3)

# 揚力傾斜の数値を表示
dCL_dalpha = np.gradient(CL_exact, alpha_rad)
slope_center = dCL_dalpha[len(dCL_dalpha) // 2]
axes[0].text(0.05, 0.92, rf"$dC_L/d\alpha = {slope_center:.3f}$ rad$^{{-1}}$"
             + f"\n(thin airfoil theory: {2*np.pi:.3f})",
             transform=axes[0].transAxes, fontsize=10,
             verticalalignment="top",
             bbox=dict(boxstyle="round", facecolor="wheat", alpha=0.5))

plt.tight_layout()
plt.savefig("joukowski_CL_vs_alpha.png", dpi=150, bbox_inches="tight")
plt.show()

print(f"キャンバー角 β = {np.degrees(beta):.3f}°")
print(f"ゼロ揚力迎え角 α_L0 = {alpha_L0_deg:.3f}°")
print(f"揚力傾斜 dC_L/dα = {slope_center:.4f} rad⁻¹ (理論値 2π = {2*np.pi:.4f})")

このグラフから、以下の重要な結果が確認できます。

1. $C_L$ は迎え角 $\alpha$ に対してほぼ線形に変化する。ジューコフスキー翼の厳密計算（青の実線）と薄翼理論の予測（赤の破線）はよく一致しています。わずかな差は翼の有限厚さに起因するもので、$\epsilon \to 0$ の極限では完全に一致します

2. 揚力傾斜 $dC_L/d\alpha$ は $2\pi$ rad$^{-1}$ に近い値です。これは薄翼理論の普遍的な予測であり、翼型の具体的な形状によらない理論的な上限値です。実際の翼では粘性効果と有限翼幅の効果により、この値よりもやや小さくなります

3. ゼロ揚力迎え角 $\alpha_{L=0}$ が負の値を取ります。これはキャンバーの効果です。キャンバーのある翼型は、迎え角ゼロでも正の揚力を発生させます。ゼロ揚力迎え角は $\alpha_{L=0} \approx -\beta$ で与えられ、これはキャンバーの大きさに比例します

4. 循環 $\Gamma$ も迎え角に対して線形に変化する。$\Gamma = 4\pi aU\sin(\alpha + \beta)$ という式の $\sin$ を小角近似すれば線形になりますが、$\pm 10°$ 程度の範囲では $\sin$ 関数との差はほとんど見えません

この線形的な $C_L$-$\alpha$ 関係は、実際の航空機翼で広く観測されており、失速角（stall angle）までの範囲で良い近似を与えます。失速はこのポテンシャル流れの理論では扱えず、境界層理論に基づく解析が必要です。

ここまでジューコフスキー翼の数値的な検証を行いました。次に、この結果の理論的な裏付けとなる薄翼理論について解説します。

薄翼理論の基礎

薄翼理論の考え方

薄翼理論（thin airfoil theory）は、翼の厚さが翼弦長に比べて十分に小さいときに成り立つ近似理論です。この理論は、翼面を渦面（vortex sheet）で置き換えるというアイデアに基づいています。

翼が薄ければ、翼面上の渦の分布を翼弦線（前縁から後縁を結ぶ直線）上に「折り畳む」ことができます。翼弦線上の各点 $x$ に渦強度 $\gamma(x)$ が分布しているとすると、翼弦線上の任意の点における法線方向の流速がゼロになるという境界条件（流体が翼面を貫通しないという条件）から、$\gamma(x)$ を決定する積分方程式が得られます。

$$ \begin{equation} \frac{1}{2\pi}\int_0^c \frac{\gamma(\xi)}{x – \xi}\,d\xi = U\left(\alpha – \frac{dy_c}{dx}\right) \end{equation} $$

ここで $c$ は翼弦長、$y_c(x)$ はキャンバー線の形状（翼弦線からの高さ）、$\alpha$ は迎え角です。右辺は「翼面に沿って流れが平行でなければならない」という条件を、翼面の傾きと迎え角で表したものです。

この積分方程式は、コーシーの主値積分を含むため一見難しそうに見えますが、変数変換 $x = \frac{c}{2}(1 – \cos\theta)$ を施すとフーリエ級数の形に帰着し、解析的に解くことができます。

対称翼の結果

キャンバーがない対称翼（$y_c = 0$）の場合、渦分布は次のように求まります。

$$ \begin{equation} \gamma(\theta) = 2U\alpha\frac{1 + \cos\theta}{\sin\theta} \end{equation} $$

この渦分布は前縁（$\theta = 0$）で無限大に発散し、後縁（$\theta = \pi$）でゼロになります。前縁での発散は、薄翼理論の前縁近似が破綻していることを示しており、実際の翼の前縁では丸みがあるため速度は有限に収まります。後縁で $\gamma = 0$ となることは、クッタ条件が自動的に満たされていることの表れです。

翼全体の循環は、渦分布の積分として得られます。

$$ \Gamma = \int_0^c \gamma(x)\,dx = \pi c U \alpha $$

したがって、クッタ・ジューコフスキーの定理により揚力係数は次のようになります。

$$ \begin{equation} C_L = \frac{2\Gamma}{Uc} = 2\pi\alpha \end{equation} $$

これが薄翼理論の最も有名な結果です。揚力傾斜 $dC_L/d\alpha = 2\pi$ は、翼型の形状によらない普遍的な値であり、先ほどのジューコフスキー翼の数値計算で確認した値と一致しています。

キャンバーのある翼の結果

キャンバー線 $y_c(x)$ がある一般的な翼型では、キャンバーの効果はゼロ揚力迎え角 $\alpha_{L=0}$ のシフトとして現れます。

$$ \begin{equation} C_L = 2\pi(\alpha – \alpha_{L=0}) \end{equation} $$

ゼロ揚力迎え角は、キャンバー線のフーリエ係数から次のように計算されます。

$$ \alpha_{L=0} = -\frac{1}{\pi}\int_0^\pi \frac{dy_c}{dx}(\cos\theta – 1)\,d\theta $$

揚力傾斜 $2\pi$ はキャンバーの有無によらず同じ値を取ります。つまり、キャンバーは $C_L$-$\alpha$ グラフを上下にシフトさせるだけで、傾きは変えません。

薄翼理論はさらに、ピッチングモーメント（空力中心まわりのモーメント）も予測します。空力中心（aerodynamic center）は翼弦の $1/4$ の位置にあり、この点まわりのモーメント係数 $C_{M,c/4}$ は迎え角に依存しない定数になるという重要な結果が得られます。この性質は翼の安定性設計において中心的な役割を果たします。

薄翼理論は数学的に美しいだけでなく、実用的にも驚くほど正確です。薄い翼型に対しては、失速角の手前まで実験結果とよく一致します。しかし、翼の厚さが大きい場合や、失速後の非線形領域では適用できません。

次に、翼の周りに循環がどのようにして「生まれる」のか、その物理的な起源である出発渦について解説します。

出発渦と循環の起源

静止状態から流れが始まるとき

ここまで、クッタ条件によって循環の値が一意に決まることを見てきました。しかし、渦度と循環で学んだように、非粘性・非圧縮流体の保存則（ケルビンの循環定理）によれば、渦なしの流れから循環は自発的に生まれません。

$$ \begin{equation} \frac{D\Gamma}{Dt} = 0 \quad \text{（ケルビンの循環定理）} \end{equation} $$

翼が静止している状態では、流れ場全体の循環はゼロです。翼が動き始めた後も、翼を含む十分大きな閉曲線に沿った循環はゼロのまま保たれなければなりません。では、翼周りの循環 $\Gamma$ はどこから来るのでしょうか。

出発渦の発生

答えは出発渦（starting vortex）にあります。翼が急に動き始めると、最初の瞬間には循環のないポテンシャル流れが形成されます。この流れではクッタ条件が満たされておらず、後縁で流れが急激に回り込もうとします。

しかし、実際の流体には粘性があるため、後縁での急激な回り込みは維持できません。後縁付近の境界層が剥離し、渦が後縁から放出されます。この渦が出発渦です。

出発渦は翼の循環と逆符号の循環を持っています。ケルビンの循環定理により、翼を含む大きな閉曲線に沿った全循環はゼロのまま保たれるため、翼周りに循環 $+\Gamma$ が生じるとき、出発渦は $-\Gamma$ の循環を持って後方に流されていきます。

$$ \begin{equation} \Gamma_{\text{wing}} + \Gamma_{\text{starting vortex}} = 0 \end{equation} $$

翼が十分な距離を移動した後、出発渦は遠くに流されてしまい、翼の近傍の流れにはほとんど影響を与えなくなります。結果として、翼のまわりにはクッタ条件を満たすちょうど正しい値の循環が残り、定常的な揚力が発生します。

出発渦の実験的観察

出発渦は理論的な概念にとどまらず、実験で直接観察できます。水槽中で翼型を急に動かすと、後縁から明瞭な渦が放出される様子が見えます。また、航空機が離陸する際にも出発渦が生成されています。この渦は後方乱気流の一成分であり、後続の航空機にとって危険な乱気流となり得ます。

出発渦の概念は、循環の発生が粘性という「非理想的な」物理を本質的に必要とすることを教えてくれます。非粘性のポテンシャル流れの理論だけでは循環は生まれず、粘性が後縁での流れの剥離を通じて循環を「作り出す」のです。クッタ条件は、この粘性効果の帰結を非粘性理論に後から組み込む、巧みな方法だったのです。

ここまでの理論的な内容を踏まえ、最後にPythonを使って出発渦と循環の生成過程、ならびに循環と揚力のクッタ・ジューコフスキーの関係を定量的に確認しましょう。

Pythonで確認する循環と揚力の定量的関係

翼周りの流線の可視化

まず、循環の有無による翼周りの流れの違いを流線で可視化します。循環がゼロの場合とクッタ条件で決まる循環がある場合を比較することで、循環の物理的な効果を視覚的に確認します。

import numpy as np
import matplotlib.pyplot as plt

# --- パラメータ ---
c0 = 1.0
epsilon = 0.10
delta = 0.08
alpha_deg = 5.0
alpha = np.radians(alpha_deg)
U = 1.0

zeta_0 = -epsilon + 1j * delta
a = abs(zeta_0 - c0)
beta = np.arctan2(delta, c0 + epsilon)

# --- 格子点（z平面）の生成 ---
nx, ny = 400, 300
x = np.linspace(-4, 4, nx)
y = np.linspace(-3, 3, ny)
X, Y = np.meshgrid(x, y)
Z = X + 1j * Y

def inverse_joukowski(z, c0):
    """ジューコフスキー変換の逆変換（外部領域の分岐を選択）"""
    return 0.5 * (z + np.sqrt(z**2 - 4 * c0**2))

def stream_function(Z, c0, zeta_0, a, U, alpha, Gamma):
    """z平面上の流れ関数を計算"""
    zeta = inverse_joukowski(Z, c0)
    # 円の内部をマスク
    inside = np.abs(zeta - zeta_0) < a * 0.98
    # 複素ポテンシャル
    w = (U * np.exp(-1j * alpha) * (zeta - zeta_0)
         + U * a**2 * np.exp(1j * alpha) / (zeta - zeta_0)
         + 1j * Gamma / (2 * np.pi) * np.log(zeta - zeta_0))
    psi = w.imag
    psi[inside] = np.nan
    return psi

# --- クッタ条件の循環 ---
Gamma_kutta = 4 * np.pi * a * U * np.sin(alpha + beta)

# --- 翼型の形状 ---
theta_af = np.linspace(0, 2 * np.pi, 500)
zeta_af = zeta_0 + a * np.exp(1j * theta_af)
z_af = zeta_af + c0**2 / zeta_af

fig, axes = plt.subplots(1, 2, figsize=(14, 5.5))

for idx, (Gamma, title) in enumerate([
    (0.0, r"No circulation ($\Gamma=0$)"),
    (Gamma_kutta, rf"Kutta condition ($\Gamma={Gamma_kutta:.2f}$)")
]):
    ax = axes[idx]
    psi = stream_function(Z, c0, zeta_0, a, U, alpha, Gamma)

    # 流線の描画
    levels = np.linspace(-3, 3, 50)
    ax.contour(X, Y, psi, levels=levels, colors="steelblue",
               linewidths=0.6, alpha=0.8)
    # 翼型
    ax.fill(z_af.real, z_af.imag, color="gray", alpha=0.8)
    ax.plot(z_af.real, z_af.imag, "k-", linewidth=1.5)

    ax.set_title(title, fontsize=13)
    ax.set_xlim(-3.5, 3.5)
    ax.set_ylim(-2.5, 2.5)
    ax.set_aspect("equal")
    ax.set_xlabel("x", fontsize=11)
    ax.set_ylabel("y", fontsize=11)
    ax.grid(True, alpha=0.2)

plt.tight_layout()
plt.savefig("joukowski_streamlines.png", dpi=150, bbox_inches="tight")
plt.show()

流線の比較から、以下の重要な違いが読み取れます。

1. 循環なし（左図）では、後縁を回り込む流れが見られる。下面から上面に向かって後縁を急激に回り込む流線が確認できます。この回り込みは後縁で無限大の速度を要求するため、実際の粘性流体では実現できません。停留点（流速ゼロの点）が上面側にあり、物理的に不自然な流れパターンです

2. クッタ条件を満たす循環（右図）では、後縁で流れが滑らかに離脱している。後縁に停留点が来ており、上面と下面の流れが後縁で合流しています。流線パターンは全体として翼の上側で密（高速）、下側で疎（低速）になっており、これが揚力を生む圧力差に対応しています

揚力の数値検証 — クッタ・ジューコフスキーの定理の確認

最後に、ジューコフスキー翼の表面上の圧力を数値積分して得られる揚力が、クッタ・ジューコフスキーの定理 $L’ = \rho U\Gamma$ と一致するかを直接確認します。これは、導出の正しさを独立に検証するテストです。

import numpy as np
import matplotlib.pyplot as plt

# --- パラメータ ---
c0 = 1.0
epsilon = 0.10
delta = 0.08
U = 1.0
rho = 1.0  # 無次元化

zeta_0 = -epsilon + 1j * delta
a = abs(zeta_0 - c0)
beta = np.arctan2(delta, c0 + epsilon)

# --- 複数の迎え角で検証 ---
alpha_list_deg = np.linspace(-8, 12, 30)
L_KJ = []      # クッタ・ジューコフスキーの定理による揚力
L_pressure = [] # 圧力分布の数値積分による揚力

n_surface = 5000  # 翼面上の分割数

for alpha_deg in alpha_list_deg:
    alpha = np.radians(alpha_deg)

    # クッタ条件から循環
    Gamma = 4 * np.pi * a * U * np.sin(alpha + beta)

    # クッタ・ジューコフスキーの定理
    L_KJ.append(rho * U * Gamma)

    # 翼面上の速度と圧力の計算
    theta = np.linspace(0, 2 * np.pi, n_surface, endpoint=False)
    zeta = zeta_0 + a * np.exp(1j * theta)

    # ζ平面の速度
    dw_dzeta = (U * np.exp(-1j * alpha)
                - U * a**2 * np.exp(1j * alpha) / (zeta - zeta_0)**2
                + 1j * Gamma / (2 * np.pi * (zeta - zeta_0)))

    # ジューコフスキー変換の微分
    dz_dzeta = 1.0 - c0**2 / zeta**2

    # z平面の速度
    dw_dz = dw_dzeta / dz_dzeta
    speed = np.abs(dw_dz)

    # ベルヌーイの定理から圧力
    p_inf = 0.0  # 無限遠方の基準圧力
    p = p_inf + 0.5 * rho * (U**2 - speed**2)

    # z平面の翼面座標
    z = zeta + c0**2 / zeta

    # 圧力の面積分（揚力 = -∮ p dy）
    dz = np.gradient(z)
    dy = dz.imag
    L_num = -np.sum(p * dy).real
    L_pressure.append(L_num)

L_KJ = np.array(L_KJ)
L_pressure = np.array(L_pressure)

fig, axes = plt.subplots(1, 2, figsize=(13, 5.5))

# 揚力の比較
ax = axes[0]
ax.plot(alpha_list_deg, L_KJ, "b-", linewidth=2.5, label=r"$L' = \rho U \Gamma$ (K-J theorem)")
ax.plot(alpha_list_deg, L_pressure, "ro", markersize=5, alpha=0.7,
        label="Pressure integration")
ax.set_xlabel(r"Angle of attack $\alpha$ [deg]", fontsize=12)
ax.set_ylabel(r"Lift per unit span $L'$", fontsize=12)
ax.set_title("Verification of Kutta-Joukowski theorem", fontsize=13)
ax.legend(fontsize=10)
ax.grid(True, alpha=0.3)
ax.axhline(0, color="k", linewidth=0.5)

# 誤差の確認
ax = axes[1]
relative_error = np.abs(L_pressure - L_KJ) / (np.abs(L_KJ) + 1e-15) * 100
ax.semilogy(alpha_list_deg, relative_error, "k-", linewidth=1.5)
ax.set_xlabel(r"Angle of attack $\alpha$ [deg]", fontsize=12)
ax.set_ylabel("Relative error [%]", fontsize=12)
ax.set_title("Relative error between two methods", fontsize=13)
ax.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig("joukowski_verification.png", dpi=150, bbox_inches="tight")
plt.show()

print(f"最大相対誤差: {np.max(relative_error):.4f}%")
print(f"平均相対誤差: {np.mean(relative_error):.4f}%")

この検証結果から、次のことが確認できます。

1. クッタ・ジューコフスキーの定理による揚力（青線）と、圧力分布の数値積分による揚力（赤丸）が正確に一致している。すべての迎え角で両者の値が重なっており、$L’ = \rho U\Gamma$ という単純な式が翼面上の圧力分布の積分と数値的に矛盾しないことが実証されています

2. 相対誤差は非常に小さい（典型的には $0.1\%$ 以下）。残っている誤差は翼面上の離散化に起因する数値誤差であり、分割数を増やせばさらに小さくなります。理論的には両者は厳密に一致するべきであり、数値計算がそれをよく裏付けています

3. 揚力は迎え角に対して線形に変化する。ゼロ揚力迎え角は負の値であり、キャンバーの効果が確認できます。この線形関係は、循環が迎え角の正弦関数で決まることの直接的な帰結です

この結果は、ブラジウスの定理を通じた導出が数値的にも正しいことを示しており、クッタ・ジューコフスキーの定理の妥当性を独立に確認するものです。

定理の適用範囲と限界

定理が成り立つ条件

クッタ・ジューコフスキーの定理は、以下の条件のもとで正確に成り立ちます。

1. 二次元流れ: 翼のスパン方向に無限に長い（端の効果がない）
2. 定常流れ: 時間変化しない
3. 非粘性流体: 粘性がない（オイラーの方程式が支配方程式）
4. 非圧縮性流体: 密度が一定
5. 翼の周りの流れがポテンシャル流れ: 渦度が翼面上と後流にのみ存在

これらの条件は実際の流れでは厳密には成り立ちません。しかし、レイノルズ数が十分に大きく、マッハ数が低い（非圧縮に近い）条件下では、定理は優れた近似を与えます。

三次元翼への拡張

実際の翼は有限の長さ（有限スパン）を持つため、翼端での流れの回り込みが生じます。翼の下面の高圧領域から上面の低圧領域へ空気が翼端を回り込むことで、翼端渦（tip vortex）が生成されます。この翼端渦は誘導抗力（induced drag）を引き起こし、揚力を減少させます。

有限翼への拡張はプラントルの揚力線理論（Prandtl’s lifting-line theory）で扱われ、揚力傾斜は次のように修正されます。

$$ \begin{equation} \frac{dC_L}{d\alpha} = \frac{2\pi}{1 + 2/AR} \end{equation} $$

ここで $AR$ はアスペクト比（翼幅$^2$/翼面積）です。アスペクト比が大きい（細長い翼の）場合は二次元の結果 $2\pi$ に近づき、小さい場合は揚力傾斜が減少します。グライダーが細長い翼を持つのは、この効果を最小化して高い揚力効率を得るためです。

粘性効果と失速

クッタ・ジューコフスキーの定理は、揚力が迎え角の増加とともに永遠に増え続けることを予測します。しかし、実際の翼では迎え角がある臨界値（失速角、典型的には $12°$ から $18°$）を超えると、翼上面の境界層が剥離し、揚力が急激に低下する失速（stall）が発生します。

失速の予測にはポテンシャル流れ理論だけでは不十分であり、粘性を含むナビエ・ストークス方程式の解析、または実験データに基づく経験的手法が必要です。

圧縮性効果

流速がマッハ数 0.3 を超えると、連続の式における密度変化が無視できなくなり、非圧縮の仮定が破綻します。亜音速域（マッハ 0.3 から 0.8 程度）では、プラントル・グラウエルト則（Prandtl-Glauert rule）による補正が適用できます。

$$ \begin{equation} C_{L,\text{compressible}} = \frac{C_{L,\text{incompressible}}}{\sqrt{1 – M_\infty^2}} \end{equation} $$

ここで $M_\infty$ は一様流のマッハ数です。超音速域ではさらに異なる理論（アッカレートの公式など）が必要になります。

これらの限界を理解しておくことで、クッタ・ジューコフスキーの定理を適切な範囲で信頼をもって使うことができます。

まとめ

本記事では、クッタ・ジューコフスキーの翼理論を基礎から解説しました。

• クッタ・ジューコフスキーの定理 $L’ = \rho U\Gamma$ は、二次元の非粘性・非圧縮流れにおいて、揚力が翼周りの循環のみで決まるという美しい結果です。この定理はブラジウスの定理と留数定理を用いて導出しました
• クッタ条件は、後縁で流速が有限であるという物理的要請であり、粘性の効果を非粘性理論に間接的に取り込む役割を果たします。この条件により循環の値が一意に決まり、揚力の予測が可能になります
• ジューコフスキー変換 $z = \zeta + c_0^2/\zeta$ は、円柱を翼型に写す等角写像です。厚さパラメータ $\epsilon$ とキャンバーパラメータ $\delta$ を調整することで、さまざまな翼形状を生成できます
• 薄翼理論は、揚力傾斜 $dC_L/d\alpha = 2\pi$ という普遍的な結果を与えます。キャンバーはゼロ揚力迎え角をシフトさせますが、揚力傾斜は変えません
• 出発渦は、翼が動き始めるときに後縁から放出される渦であり、ケルビンの循環定理を満たしつつ翼周りに循環を生成するメカニズムです
• Pythonによる可視化で、ジューコフスキー変換の幾何学的な動作、翼型の圧力分布、循環と揚力の線形関係、および定理の数値的検証を確認しました

クッタ・ジューコフスキーの定理は1世紀以上前に確立された理論ですが、その物理的な洞察は現代の航空工学においても設計の出発点であり続けています。循環という単一の量が揚力を決めるという深い結果は、流体力学と複素関数論の美しい交差点です。

次のステップとして、以下の記事も参考にしてください。

• ポテンシャル流れ — 円柱まわりの流れなどポテンシャル流れの基礎理論
• 境界層理論 — クッタ条件の物理的根拠となる粘性効果の理論
• 抗力と揚力 — 揚力と抗力のメカニズムの全体像
• ストークスの定理 — 循環と渦度の関係を結ぶ数学的基盤