円筒座標系での勾配・発散・回転・ラプラシアン

物理問題には円筒対称性を持つものが多く存在します。同軸ケーブル、円柱導体周りの磁場、管内の流れなどでは、直交座標系よりも円筒座標系で記述するほうが自然で計算も簡潔になります。

本記事では、円筒座標系での微分演算子(勾配・発散・回転・ラプラシアン)を計量係数の概念から導出します。

本記事の内容

  • 円筒座標系の定義と座標変換
  • 計量係数(スケール因子)
  • 勾配 $\nabla f$ の導出
  • 発散 $\nabla\cdot\bm{A}$ の導出
  • 回転 $\nabla\times\bm{A}$ の導出
  • ラプラシアン $\nabla^2 f$ の導出
  • 応用例: 同軸ケーブルの電場
  • Pythonでの可視化

前提知識

円筒座標系の定義

円筒座標系 $(r, \theta, z)$ と直交座標系 $(x, y, z)$ の関係は:

$$ \boxed{x = r\cos\theta, \quad y = r\sin\theta, \quad z = z} $$

ここで $r \geq 0$, $0 \leq \theta < 2\pi$, $-\infty < z < \infty$ です。

逆変換は:

$$ r = \sqrt{x^2 + y^2}, \quad \theta = \arctan\frac{y}{x}, \quad z = z $$

単位ベクトル

各座標方向の単位ベクトルは位置に依存します:

$$ \bm{e}_r = \cos\theta\,\bm{e}_x + \sin\theta\,\bm{e}_y $$

$$ \bm{e}_\theta = -\sin\theta\,\bm{e}_x + \cos\theta\,\bm{e}_y $$

$$ \bm{e}_z = \bm{e}_z $$

重要な性質として、$\bm{e}_r$ と $\bm{e}_\theta$ は位置 $\theta$ に依存して向きが変わります。したがって:

$$ \frac{\partial\bm{e}_r}{\partial\theta} = \bm{e}_\theta, \quad \frac{\partial\bm{e}_\theta}{\partial\theta} = -\bm{e}_r $$

計量係数(スケール因子)

微小変位ベクトルは:

$$ d\bm{r} = \frac{\partial\bm{r}}{\partial r}dr + \frac{\partial\bm{r}}{\partial\theta}d\theta + \frac{\partial\bm{r}}{\partial z}dz $$

位置ベクトル $\bm{r} = r\cos\theta\,\bm{e}_x + r\sin\theta\,\bm{e}_y + z\,\bm{e}_z$ を微分すると:

$$ \frac{\partial\bm{r}}{\partial r} = \cos\theta\,\bm{e}_x + \sin\theta\,\bm{e}_y = \bm{e}_r $$

$$ \frac{\partial\bm{r}}{\partial\theta} = -r\sin\theta\,\bm{e}_x + r\cos\theta\,\bm{e}_y = r\,\bm{e}_\theta $$

$$ \frac{\partial\bm{r}}{\partial z} = \bm{e}_z $$

したがって微小変位は:

$$ d\bm{r} = dr\,\bm{e}_r + r\,d\theta\,\bm{e}_\theta + dz\,\bm{e}_z $$

計量係数(スケール因子)は各座標方向の微小変位の大きさです:

$$ \boxed{h_r = 1, \quad h_\theta = r, \quad h_z = 1} $$

微小体積要素は:

$$ dV = h_r h_\theta h_z \, dr\,d\theta\,dz = r\,dr\,d\theta\,dz $$

勾配(grad)

直交曲線座標系での勾配の一般公式は:

$$ \nabla f = \frac{1}{h_r}\frac{\partial f}{\partial r}\bm{e}_r + \frac{1}{h_\theta}\frac{\partial f}{\partial\theta}\bm{e}_\theta + \frac{1}{h_z}\frac{\partial f}{\partial z}\bm{e}_z $$

計量係数を代入すると:

$$ \boxed{\nabla f = \frac{\partial f}{\partial r}\bm{e}_r + \frac{1}{r}\frac{\partial f}{\partial\theta}\bm{e}_\theta + \frac{\partial f}{\partial z}\bm{e}_z} $$

導出

スカラー場 $f$ の全微分は:

$$ df = \frac{\partial f}{\partial r}dr + \frac{\partial f}{\partial\theta}d\theta + \frac{\partial f}{\partial z}dz $$

一方、$df = \nabla f \cdot d\bm{r}$ であり、$d\bm{r} = dr\,\bm{e}_r + r\,d\theta\,\bm{e}_\theta + dz\,\bm{e}_z$ なので:

$$ df = (\nabla f)_r \, dr + (\nabla f)_\theta \cdot r\,d\theta + (\nabla f)_z \, dz $$

比較すると $(\nabla f)_r = \frac{\partial f}{\partial r}$, $(\nabla f)_\theta = \frac{1}{r}\frac{\partial f}{\partial\theta}$, $(\nabla f)_z = \frac{\partial f}{\partial z}$ が得られます。

発散(div)

直交曲線座標系での発散の一般公式は:

$$ \nabla\cdot\bm{A} = \frac{1}{h_r h_\theta h_z}\left[\frac{\partial}{\partial r}(h_\theta h_z A_r) + \frac{\partial}{\partial\theta}(h_r h_z A_\theta) + \frac{\partial}{\partial z}(h_r h_\theta A_z)\right] $$

計量係数を代入すると:

$$ \boxed{\nabla\cdot\bm{A} = \frac{1}{r}\frac{\partial}{\partial r}(rA_r) + \frac{1}{r}\frac{\partial A_\theta}{\partial\theta} + \frac{\partial A_z}{\partial z}} $$

導出

ガウスの発散定理の微分形から導出します。微小体積要素 $\Delta V = r\,\Delta r\,\Delta\theta\,\Delta z$ を考えます。

$r$ 方向のフラックス差は:

$$ [A_r(r+\Delta r)(r+\Delta r)\Delta\theta\Delta z – A_r(r) \cdot r\Delta\theta\Delta z] \approx \frac{\partial}{\partial r}(rA_r)\Delta r\Delta\theta\Delta z $$

$\theta$ 方向のフラックス差は:

$$ [A_\theta(\theta+\Delta\theta)\Delta r\Delta z – A_\theta(\theta)\Delta r\Delta z] \approx \frac{\partial A_\theta}{\partial\theta}\Delta r\Delta\theta\Delta z $$

$z$ 方向のフラックス差は:

$$ [A_z(z+\Delta z)r\Delta r\Delta\theta – A_z(z)r\Delta r\Delta\theta] \approx \frac{\partial A_z}{\partial z}r\Delta r\Delta\theta\Delta z $$

合計を $\Delta V = r\Delta r\Delta\theta\Delta z$ で割ると上記の結果が得られます。

回転(curl)

直交曲線座標系での回転の一般公式は:

$$ \nabla\times\bm{A} = \frac{1}{h_r h_\theta h_z}\begin{vmatrix} h_r\bm{e}_r & h_\theta\bm{e}_\theta & h_z\bm{e}_z \\ \frac{\partial}{\partial r} & \frac{\partial}{\partial\theta} & \frac{\partial}{\partial z} \\ h_r A_r & h_\theta A_\theta & h_z A_z \end{vmatrix} $$

計量係数を代入すると:

$$ \boxed{\nabla\times\bm{A} = \left(\frac{1}{r}\frac{\partial A_z}{\partial\theta} – \frac{\partial A_\theta}{\partial z}\right)\bm{e}_r + \left(\frac{\partial A_r}{\partial z} – \frac{\partial A_z}{\partial r}\right)\bm{e}_\theta + \frac{1}{r}\left(\frac{\partial(rA_\theta)}{\partial r} – \frac{\partial A_r}{\partial\theta}\right)\bm{e}_z} $$

特に $z$ 成分の $\frac{1}{r}\frac{\partial(rA_\theta)}{\partial r}$ に注意してください。これは直交座標系の場合($\frac{\partial A_\theta}{\partial r}$)とは異なります。

ラプラシアン

スカラーラプラシアン $\nabla^2 f = \nabla\cdot(\nabla f)$ に発散と勾配の公式を代入します。

$\nabla f$ の $r$ 成分は $\frac{\partial f}{\partial r}$ なので、$\nabla\cdot(\nabla f)$ の第1項は:

$$ \frac{1}{r}\frac{\partial}{\partial r}\left(r\frac{\partial f}{\partial r}\right) $$

同様に他の項も計算すると:

$$ \boxed{\nabla^2 f = \frac{1}{r}\frac{\partial}{\partial r}\left(r\frac{\partial f}{\partial r}\right) + \frac{1}{r^2}\frac{\partial^2 f}{\partial\theta^2} + \frac{\partial^2 f}{\partial z^2}} $$

第1項を展開すると:

$$ \frac{1}{r}\frac{\partial}{\partial r}\left(r\frac{\partial f}{\partial r}\right) = \frac{\partial^2 f}{\partial r^2} + \frac{1}{r}\frac{\partial f}{\partial r} $$

直交座標系のラプラシアン $\frac{\partial^2}{\partial x^2} + \frac{\partial^2}{\partial y^2} + \frac{\partial^2}{\partial z^2}$ と比べると、$\frac{1}{r}\frac{\partial f}{\partial r}$ の項が余分に現れます。

応用例: 同軸ケーブルの電場

内半径 $a$、外半径 $b$ の同軸ケーブルで、内導体に線電荷密度 $\lambda$ [C/m] が与えられた場合を考えます。

対称性から電場は $\bm{E} = E_r(r)\,\bm{e}_r$ の形です。ガウスの法則の微分形 $\nabla\cdot\bm{E} = \rho/\epsilon_0$ を使います。

導体間($a < r < b$)では $\rho = 0$ なので:

$$ \nabla\cdot\bm{E} = \frac{1}{r}\frac{\partial}{\partial r}(rE_r) = 0 $$

$$ \frac{\partial}{\partial r}(rE_r) = 0 \implies rE_r = C \implies E_r = \frac{C}{r} $$

境界条件(ガウスの法則の積分形)から $C = \lambda/(2\pi\epsilon_0)$ と決まり:

$$ \boxed{E_r = \frac{\lambda}{2\pi\epsilon_0 r} \quad (a < r < b)} $$

Pythonでの可視化

import numpy as np
import matplotlib.pyplot as plt

# --- (1) 円筒座標系の基底ベクトルの可視化 ---
fig, axes = plt.subplots(1, 3, figsize=(18, 5))

theta_pts = np.linspace(0, 2*np.pi, 8, endpoint=False)
for th in theta_pts:
    x0, y0 = 1.5*np.cos(th), 1.5*np.sin(th)
    # e_r
    axes[0].arrow(x0, y0, 0.4*np.cos(th), 0.4*np.sin(th),
                  head_width=0.08, color='blue', lw=1.5)
    # e_theta
    axes[0].arrow(x0, y0, -0.4*np.sin(th), 0.4*np.cos(th),
                  head_width=0.08, color='red', lw=1.5)

circle = plt.Circle((0, 0), 1.5, fill=False, ls='--', color='gray')
axes[0].add_patch(circle)
axes[0].set_xlim(-2.5, 2.5); axes[0].set_ylim(-2.5, 2.5)
axes[0].set_aspect('equal')
axes[0].set_title('Unit vectors $\\hat{e}_r$ (blue), $\\hat{e}_\\theta$ (red)', fontsize=12)
axes[0].set_xlabel('$x$'); axes[0].set_ylabel('$y$')
axes[0].grid(True, alpha=0.3)

# --- (2) 同軸ケーブルの電場 ---
a, b = 0.5, 2.0  # 内半径, 外半径
lam = 1e-9  # 線電荷密度 [C/m]
eps0 = 8.854e-12

r = np.linspace(0.1, 3.0, 500)
E_r = np.where((r >= a) & (r <= b), lam / (2*np.pi*eps0*r), 0)
# 内導体内部・外導体外部は 0

axes[1].plot(r, E_r, 'b-', lw=2)
axes[1].axvline(a, color='gray', ls='--', alpha=0.5, label=f'$a={a}$')
axes[1].axvline(b, color='gray', ls='-.', alpha=0.5, label=f'$b={b}$')
axes[1].fill_between([0, a], 0, max(E_r)*1.1, alpha=0.1, color='orange', label='Inner conductor')
axes[1].fill_between([b, 3], 0, max(E_r)*1.1, alpha=0.1, color='green', label='Outer conductor')
axes[1].set_xlabel('$r$ [m]', fontsize=12)
axes[1].set_ylabel('$E_r$ [V/m]', fontsize=12)
axes[1].set_title('Electric field in coaxial cable', fontsize=13)
axes[1].legend(fontsize=9)
axes[1].grid(True, alpha=0.3)
axes[1].set_ylim(0, max(E_r)*1.1)

# --- (3) 電場の2次元ベクトルプロット ---
r_2d = np.linspace(a, b, 10)
theta_2d = np.linspace(0, 2*np.pi, 16, endpoint=False)
R, Th = np.meshgrid(r_2d, theta_2d)
X2 = R * np.cos(Th)
Y2 = R * np.sin(Th)
Er = lam / (2*np.pi*eps0*R)
Ex = Er * np.cos(Th)
Ey = Er * np.sin(Th)

circle_a = plt.Circle((0, 0), a, fill=True, color='orange', alpha=0.3)
circle_b = plt.Circle((0, 0), b, fill=False, ls='--', color='green', lw=2)
axes[2].add_patch(circle_a)
axes[2].add_patch(circle_b)
axes[2].quiver(X2, Y2, Ex/np.max(Er), Ey/np.max(Er), color='blue', alpha=0.7)
axes[2].set_xlim(-2.5, 2.5); axes[2].set_ylim(-2.5, 2.5)
axes[2].set_aspect('equal')
axes[2].set_title('Electric field (2D cross section)', fontsize=13)
axes[2].set_xlabel('$x$'); axes[2].set_ylabel('$y$')
axes[2].grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('cylindrical_coordinates.png', dpi=150, bbox_inches='tight')
plt.show()

まとめ

本記事では、円筒座標系 $(r, \theta, z)$ での微分演算子を導出しました。

  • 座標変換: $x = r\cos\theta$, $y = r\sin\theta$, $z = z$
  • 計量係数: $h_r = 1$, $h_\theta = r$, $h_z = 1$
  • 勾配: $\nabla f = \frac{\partial f}{\partial r}\bm{e}_r + \frac{1}{r}\frac{\partial f}{\partial\theta}\bm{e}_\theta + \frac{\partial f}{\partial z}\bm{e}_z$
  • 発散: $\nabla\cdot\bm{A} = \frac{1}{r}\frac{\partial}{\partial r}(rA_r) + \frac{1}{r}\frac{\partial A_\theta}{\partial\theta} + \frac{\partial A_z}{\partial z}$
  • ラプラシアン: $\nabla^2 f = \frac{1}{r}\frac{\partial}{\partial r}\left(r\frac{\partial f}{\partial r}\right) + \frac{1}{r^2}\frac{\partial^2 f}{\partial\theta^2} + \frac{\partial^2 f}{\partial z^2}$
  • 同軸ケーブルの電場: $E_r = \frac{\lambda}{2\pi\epsilon_0 r}$($1/r$ に比例)

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