ポテンシャル流れを描く（python, matplotlib, 複素関数）

概要

複素ポテンシャル流れを描くスクリプトを作成しました。
与えられた複素関数に基づいて複素平面上にポテンシャル流線を描画し、複素積分を利用して流れ場の特性を可視化します。
以下のサイトを参考にしました。
https://qiita.com/Helmet/items/50f92cbd6fd212e43a2c

以下は概要の数式および結果の一例です。

円柱周り流れ

流速をU, 円柱半径をRとして

複素ポテンシャル：W(z)=U(Z+\frac{R^2}{Z})

共役複素速度：w(z)=\frac{dW}{dz}=U(Z-(\frac{R}{Z})^2)

pythonの場合、42行あたりを以下のように書きます

コードの41行～

# 被積分複素関数（複素ポテンシャル）の定義
def W_potential(Z):
    return U_velocity*(Z+radius**2/Z) # 二重湧きだし

# 共役複素速度の定義
def w_velocity(Z):
    return U_velocity*(np.ones_like(Z)-(radius/Z)**2) # 二重湧きだし

出力はこのようになります。

元文献と同じく、クリックすれば積分が出来ます。
赤いベクトルはその点の流速です。

出力

sum w(z)*dz = (Γ:-0.09) + (Q:-0.22)*i
sum w(z)*dz = (Γ:-0.14) + (Q:-0.29)*i
sum w(z)*dz = (Γ:-0.16) + (Q:-0.38)*i
sum w(z)*dz = (Γ:-0.17) + (Q:-0.48)*i
sum w(z)*dz = (Γ:-0.16) + (Q:-0.58)*i
sum w(z)*dz = (Γ:-0.07) + (Q:-0.77)*i
sum w(z)*dz = (Γ:0.02) + (Q:-0.88)*i
sum w(z)*dz = (Γ:0.09) + (Q:-1.06)*i

流量Qと循環Γ

流量をQ, 循環をGとして

湧き出し複素ポテンシャル：W(z)=\frac{Q}{2\pi} \ln(Z)

湧き出し共役複素速度：w(z)=\frac{dW}{dz}=\frac{Q}{2\pi} \ln(Z)

自由渦複素ポテンシャル：W(z)=\frac{Q}{2\pi} \ln(Z)

自由渦共役複素速度：w(z)=\frac{dW}{dz}=\frac{Q}{2\pi} \ln(Z)

pythonの場合、42行あたりを以下のように書きます

コードの41行～

# 被積分複素関数（複素ポテンシャル）の定義
def W_potential(Z):
    return Q/(2*np.pi) * np.log(Z) # 吸い込みQ<0・湧きだしQ>0
    # return 1j*(-G)/(2*np.pi) * np.log(Z) # 自由渦

# 共役複素速度の定義
def w_velocity(Z):
    return Q/(2*np.pi)*1/Z # 吸い込みQ<0・湧きだしQ>0```
    # return 1j*-G/(2*np.pi)*1/Z # 自由渦

Q=1, G=1としました。
出力はこのようになります。
湧き出し：

自由渦：

積分結果

湧き出し：sum w(z)*dz = (Γ:-0.12) + (Q:0.97)*i
自由渦　：sum w(z)*dz = (Γ:0.99) + (Q:0.15)*i

以上より、
複素共役速度の積分=複素ポテンシャルの実部と虚部が
循環G=1と流量Q=1
に大体一致したので合ってそうかなと思います。

使用したモジュール

pip

pip install matplotlib
pip install japanize-matplotlib
pip install numpy

書いたスクリプト

complex_fluid.py

import japanize_matplotlib
import numpy as np
import matplotlib.pyplot as plt

# TODO: リセット・道の保存・循環とかのテキスト表示

###################################
"準備"
LX, LY=[-2,2],[-2,2] #描画範囲

gridwidth=0.01 #矢印を描く密度

# 画面準備(1行*2列個)
fig, ax = plt.subplots(1, 2)
fig.set_size_inches(15, 7) #画像サイズ

# タイトルを付ける
ax[0].set_title('複素ポテンシャル:W = f(z)')
ax[1].set_title('共役複素速度:w = df(z)/dz = u-vj')
# 格子を加える
ax[0].grid()
ax[1].grid()
# 範囲指定
ax[0].set_xlim(LX)
ax[0].set_ylim(LY)
ax[1].set_xlim(LX)
ax[1].set_ylim(LY)

###################################
"ベクトル場作成"
# 変数定義
U_velocity = 1 # 一様流流速 
a = np.pi/3 # 一様流角度
Q = 1 # 湧きだし流量
G = 1 # 循環
radius = 1 # 円柱周り流れ半径

# 被積分複素関数（複素ポテンシャル）の定義
def W_potential(Z):
    # return U_velocity*np.exp(1j*-a)*Z # 一様流
    # return Q/(2*np.pi) * np.log(Z) # 吸い込みQ<0・湧きだしQ>0
    return 1j*(-G)/(2*np.pi) * np.log(Z) # 自由渦
    # return U_velocity*(Z+radius**2/Z) # 円柱周り流れ
    # return U_velocity*(Z+radius**2/Z) + 1j*(-G)/(2*np.pi) * np.log(Z)# 回転あり円柱周り流れ

# 共役複素速度の定義
def w_velocity(Z):
    # return U_velocity*np.exp(1j*-a)*np.ones_like(Z) # 一様流
    # return Q/(2*np.pi)*1/Z # 吸い込みQ<0・湧きだしQ>0
    return 1j*-G/(2*np.pi)*1/Z # 自由渦
    # return U_velocity*(np.ones_like(Z)-(radius/Z)**2) # 円柱周り流れ
    # return U_velocity*(np.ones_like(Z)+radius**2/Z) + 1j*-G/(2*np.pi)*1/Z # 回転あり円柱周り流れ

#零点(特異点)は手打ちする
Xzero, Yzero = 0, 0

X, Y= np.meshgrid(np.arange(LX[0], LX[1], gridwidth), np.arange(LY[0], LY[1], gridwidth))
Z = (X + Y*1j)
W = W_potential(Z) # 複素ポテンシャル
U = W.real
V = W.imag
w = w_velocity(Z) # 共役複素速度
u = w.real
v = -1*w.imag

#正規化したかったら以下もつける
R = np.sqrt(U**2 + V**2)
U = U / R
V = V / R

r = np.sqrt(u**2 + v**2)
u = u / r
v = v / r


###################################
"複素平面描画"
# 
ax[0].plot(Xzero,Yzero,'o',color='black')
strm = ax[0].streamplot(X, Y, U, V, color=R, linewidth=2, cmap='summer')
# plt.colorbar(strm.lines) # 複素ポテンシャルのグラフはそんなに重要じゃないので。

ax[1].plot(Xzero,Yzero,'o',color='black')
strm = ax[1].streamplot(X, Y, u, v, color=r, linewidth=2, cmap='summer')
plt.colorbar(strm.lines)

"クリックイベント（積分）"
C_start = []
C_end = []
dz = []
sum_fz_dz = 0+0j
def draw_C_vec(event):
    global C_start
    global C_end
    global dz
    global sum_fz_dz
    
    if event.button == 1:
    # if event.button == 2: # ホイールクリックなら
        # 始点取得
        if C_start == []: 
            C_start = [event.xdata, event.ydata]
            ax[1].plot(C_start[0], C_start[1],
                       "o", color="blue")

            # 更新
            fig.canvas.draw()
            print("get_start")

        # 終点取得＋処理
        else:
            C_end = [event.xdata, event.ydata]
            # dzベクトル取得
            dz = (C_end[0]-C_start[0]) + (C_end[1]-C_start[1]) * 1j

            # 要素の計算 (z_i, f(z), f(z)*dz)
            z = C_start[0]+C_start[1]*1j
            f_z = w_velocity(z)

            fz_dz = f_z * dz


            " 以下グラフ描画 "
            # 経路C(dz_i)をプロット
            ax[1].arrow(C_start[0],C_start[1],
                        dz.real, dz.imag,
                        width=0.01,color="blue")
            # クリック点プロット
            ax[1].plot(C_end[0], C_end[1],
                        "o", color="blue")
            # f(z)ベクトル = 流速ベクトル
            ax[1].arrow(C_start[0],C_start[1],
                        f_z.real,  -f_z.imag,
                        width=0.01,color="red")

            # 積分結果(循環と流量)計算
            sum_fz_dz += fz_dz

            # 3桁までconsoleに表示
            print(f"sum w(z)*dz = (Γ:{sum_fz_dz.real:.2f}) + (Q:{sum_fz_dz.imag:.2f})*i")
            
            # グラフ更新
            fig.canvas.draw()
            # 値の更新
            C_start = C_end.copy()
            
# クリックに反応するよう設定
cid = fig.canvas.mpl_connect('button_press_event', draw_C_vec)

# 描写開始
plt.draw()
plt.show()