Blasiusの第一公式

循環をもつ一様流が円柱を過ぎる2次元ポテンシャル流れに対して Blasius の第一公式を使って円柱に働く力を計算してみます。

複素速度ポテンシャルで流れを表す

複素変数を $z = x + \mathrm{i}y$ 、一様流速 $U$ 、円柱半径 $a$ 、循環 $\Gamma$ 、迎角 $\alpha$ とおくと、完全流体中で循環をもつ一様流が円柱にあたるときの流れを表す複素速度ポテンシャルは以下の式で表される。

f(z) = U\left( z e^{-\mathrm{i}\alpha} + \frac{a^2 e^{\mathrm{i}\alpha}}{z} \right) + \frac{\mathrm{i}\Gamma}{2\pi}\log\left(\frac{z}{a}\right)

複素速度は複素速度ポテンシャルを微分することにより

\frac{df}{dz} = U\left( e^{-\mathrm{i}\alpha} - \frac{a^2 e^{\mathrm{i}\alpha}}{z^2} \right) + \frac{\mathrm{i}\Gamma}{2\pi z}

で与えられる。

Blasiusの第一公式で力を計算する

$\rho$ を流体密度、複素速度を $\frac{df}{dz}=u-iv$ 、積分路 $C$ を反時計回りとすると、Blasius の第一公式は以下の式で定義される。

F_x - \mathrm{i}F_y = \frac{\mathrm{i}\rho}{2} \oint_C \left(\frac{df}{dz}\right)^2 dz

計算の際には以下のパラメータを設定した。

一様流速 $U = 1.0$
円柱半径 $a = 1.0$
循環 $\Gamma =4\pi Ua$ ( $>0$ のとき円柱まわりの流れは時計回りで揚力は正の $y$ 方向に働く)
迎角 $\alpha = 0^\circ$
密度 $\rho = 1.0$

計算結果は以下のようになった。

$F_x: 0$
$F_y: 12.5663$

なお、ここでの $F_x,F_y$ は固定座標系での力の成分である。抗力・揚力を一様流方向に対して定義すると、抗力は常に 0 (ダランベールのパラドックス)、揚力の大きさは $\rho U\Gamma$ (Kutta–Joukowskiの定理) になる。

数値的には、円柱表面または円柱を囲む任意の円周を離散化し、 $\oint W^2dz$ を台形公式などで近似できる。被積分関数は円柱外部で解析的なので、円柱境界を横切らず、かつ円柱を一回囲むという条件を保つ限り、積分路を変形しても周回積分の値は変わらない。例えば、円柱表面ではなく、より大きな半径を持つ円周上で積分しても同じ結果が得られる。(円柱内部にある特異点の留数により積分結果が決まる。)

一方、今回のように複素速度 $\frac{df}{dz}$ が明示的に分かっている場合は、 $\left(\frac{df}{dz}\right)^2$ の Laurent 展開から $1/z$ の係数を特定することで解析的に計算できる。

W = \frac{df}{dz} = U\left( e^{-\mathrm{i}\alpha} - \frac{a^2 e^{\mathrm{i}\alpha}}{z^2} \right) + \frac{\mathrm{i}\Gamma}{2\pi z}

とおくと、

W^2 = U^2 e^{-2\mathrm{i}\alpha} + Ue^{-\mathrm{i}\alpha}\frac{\mathrm{i}\Gamma}{\pi z} + O\left( \frac{1}{z^2} \right)

留数は $Ue^{-\mathrm{i}\alpha}\frac{\mathrm{i}\Gamma}{\pi}$ なのでBlasiusの公式に代入すると

F_x - \mathrm{i}F_y = \frac{\mathrm{i}\rho}{2} \cdot2\pi \mathrm{i}\cdot Ue^{-\mathrm{i}\alpha}\frac{\mathrm{i}\Gamma}{\pi} = -\mathrm{i}\rho U \Gamma e^{-\mathrm{i}\alpha}

= -\mathrm{i}\rho U \Gamma \left( \cos\alpha - \mathrm{i}\sin\alpha \right)

= -\rho U \Gamma \sin\alpha-\mathrm{i}\rho U \Gamma \cos\alpha

\therefore F_x = -\rho U \Gamma \sin\alpha, F_y = \rho U \Gamma \cos\alpha

迎角 $\alpha = 0^\circ$ の場合は $F_x = 0$ , $F_y = \rho U \Gamma$ となり、働く力を解析的に計算することができた。

また、循環の量を変えながら流線と速度場を描いてみた。

迎角0度、循環3πUaの円柱まわりの流線 — **$\Gamma = 3\pi Ua$**

迎角0度、循環4πUaの円柱まわりの流線 — **$\Gamma = 4\pi Ua$**

迎角0度、循環5πUaの円柱まわりの流線 — **$\Gamma = 5\pi Ua$**

迎角0度、循環3πUaの円柱まわりの速度場 — **$\Gamma = 3\pi Ua$**

迎角0度、循環4πUaの円柱まわりの速度場 — **$\Gamma = 4\pi Ua$**

迎角0度、循環5πUaの円柱まわりの速度場 — **$\Gamma = 5\pi Ua$**

迎角45°の場合。

迎角45度、循環3πUaの円柱まわりの流線 — **$\Gamma = 3\pi Ua$**

迎角45度、循環4πUaの円柱まわりの流線 — **$\Gamma = 4\pi Ua$**

迎角45度、循環5πUaの円柱まわりの流線 — **$\Gamma = 5\pi Ua$**

迎角45度、循環3πUaの円柱まわりの速度場 — **$\Gamma = 3\pi Ua$**

迎角45度、循環4πUaの円柱まわりの速度場 — **$\Gamma = 4\pi Ua$**

迎角45度、循環5πUaの円柱まわりの速度場 — **$\Gamma = 5\pi Ua$**

まとめ

実際に計算してみて初めてその式を自分が自由に扱える感覚が身につくような気がする。
“Allez en avant, et la foi vous viendra” (前進すれば信仰が訪れる) はダランベールの言葉らしい。
圧力分布も計算したので時間があれば今後追記します。

今回のコード

今回の記事で使用したコードを載せておきます。

1
# -*- coding: utf-8 -*-
2

3
import numpy as np
4
import matplotlib.pyplot as plt
5
from scipy.integrate import trapezoid
6

7
# フォント指定(文字化け対策)
8
plt.rcParams['font.family'] = 'Meiryo'
9

10
# パラメータ設定
11
U = 1.0                 # 一様流速
12
a = 1.0                 # 円柱の半径
13
Gamma = 4*np.pi*U*a     # 循環
14
alpha = np.radians(0)   # 迎角
15
#alpha = np.radians(45) # 迎角
16
rho = 1.0               # 密度
17
N = 100                 # 格子数
18

19
# 複素速度ポテンシャル関数
20
def complex_potential(z, U, a, Gamma, alpha):
21
    f = U * (z * np.exp(-1j * alpha) + a**2 * np.exp(1j * alpha) / z) + (1j * Gamma / (2 * np.pi)) * np.log(z / a)
22
    return f
23

24
# 複素速度関数 (df/dz)
25
def complex_velocity(z, U, a, Gamma, alpha):
26
    df_dz = U * (np.exp(-1j * alpha) - a**2 * np.exp(1j * alpha) / z**2) + (1j * Gamma / (2 * np.pi * z))
27
    return df_dz
28

29
x = np.linspace(-4, 4, N)
30
y = np.linspace(-4, 4, N)
31
X, Y = np.meshgrid(x, y)
32
Z = X + 1j * Y
33

34
# 円柱内部へのマスク処理
35
mask = np.abs(Z) > a
36
Z_mask = np.where(mask, Z, np.nan)
37

38
# 複素速度ポテンシャル関数
39
f = complex_potential(Z_mask, U, a, Gamma, alpha)
40
phi = np.real(f) # 速度ポテンシャル
41
psi = np.imag(f) # 流れ関数
42

43
plt.figure(figsize=(8, 8))
44
plt.contour(X, Y, psi, levels=50, cmap='viridis')
45
circle = plt.Circle((0, 0), a, color='black', fill=False)
46
plt.gca().add_artist(circle)
47
plt.xlim(-4, 4)
48
plt.ylim(-4, 4)
49
plt.gca().set_aspect('equal')
50
plt.title(f'循環のある円柱を過ぎる流れ(alpha=0, $\Gamma$=4$\pi$Ua)')
51
plt.xlabel('x')
52
plt.ylabel('y')
53

54
plt.show()
55

56
# ブラジウスの第一公式の適用
57
theta = np.linspace(0, 2 * np.pi, 1000)
58
z_surf = a * np.exp(1j * theta)
59
df_dz = complex_velocity(z_surf, U, a, Gamma, alpha)
60
force_complex = (1j * rho / 2) * trapezoid((df_dz**2), z_surf)
61

62
Fx = force_complex.real
63
Fy = -force_complex.imag
64

65
print(f"Fx: {Fx:.4f}")
66
print(f"Fy: {Fy:.4f}")