【第15回Python流体の数値計算】2次元ポアソン方程式をPythonで実装する。

こんにちは(@t_kun_kamakiri)。

前回の記事では、「2次元のラプラス方程式」をGoogle Colab上でPythonで実装しました。

前回の記事はこちら

【第14回Python流体の数値計算】2次元ラプラス方程式をPythonで実装する。

今回は、以下のような２次元のポアソン方程式をPythonで実装します。

なぜ、ポアソン方程式の解法を知る必要があるのかというと、ナビエストークス方程式を解く際に必要になるからです。（前回の記事で詳しく書きました。）

本件の基本的な内容はこちらのサイトにそってやっていきます。

この記事ではこんな人を対象にしています。

こんな人が対象

Pythonを使い始めたけどどう使うかわからない
流体の数値計算をはじめて勉強する人

本記事シリーズの最終目標は、ナビエストークス方程式をPythonで実装することですが、いきなりナビエストークス方程式を実装するのは難しいので各項の意味を確認しながら進めていきたいと思います。

今回の内容はこちら

２次元のポアソン方程式の数値計算をPythonで実装

カマキリ

今日もゆる～く学びます(^^)/

ポアソン方程式
ポアソン方程式を離散化
境界条件を設定
ポアソン方程式をPythonで実装
コンター図も作成してみる
まとめ
おすすめの書籍紹介

ポアソン方程式

まずは、

ポアソン方程式がどういう式なのか
ナビエストークス方程式のを解く際のポアソン方程式の使われ方

を復習しておきましょう。

流体現象を扱う際に「ナビエストークス方程式」と「質量保存則」を連立して解く必要があります。

ナビエストークス方程式

\begin{align*}
\frac{\partial u}{\partial t}+u\frac{\partial u}{\partial x}+v\frac{\partial u}{\partial y} = -\frac{1}{\rho}\frac{\partial p}{\partial x}+\nu \left(\frac{\partial^2 u}{\partial x^2}+\frac{\partial^2 u}{\partial y^2} \right)\tag{1}
\end{align*}

\begin{align*}
\frac{\partial v}{\partial t}+u\frac{\partial v}{\partial x}+v\frac{\partial v}{\partial y} = -\frac{1}{\rho}\frac{\partial p}{\partial y}+\nu\left(\frac{\partial^2 v}{\partial x^2}+\frac{\partial^2 v}{\partial y^2}\right) \tag{2}
\end{align*}

質量保存則（非圧縮条件$\nabla\cdot\boldsymbol{v}=0$から導出）

\begin{align*}
\frac{\partial^2 p}{\partial x^2}+\frac{\partial^2 p}{\partial y^2} = b(x,y)\tag{3.1}
\end{align*}

\begin{align*}
b(x,y)=-\rho\left(\frac{\partial u}{\partial x}\frac{\partial u}{\partial x}+2\frac{\partial u}{\partial y}\frac{\partial v}{\partial x}+\frac{\partial v}{\partial y}\frac{\partial v}{\partial y} \right)\tag{3.2}
\end{align*}

(3.1)式の形がポアソン方程式

このように、ナビエストークス方程式とポアソン方程式を連立させて解く必要があります。

右辺が0である場合のラプラス方程式、

\begin{align*}
\frac{\partial ^2 p}{\partial x^2} + \frac{\partial ^2 p}{\partial y^2} = 0
\end{align*}

は、前回の記事で扱いましたが、今回は以下のポアソン方程式を扱います。

ポアソン方程式

\begin{align*}
\frac{\partial ^2 p}{\partial x^2} + \frac{\partial ^2 p}{\partial y^2} = b(x,y)
\end{align*}

実際の流体の方程式としては、一般的には右辺が0ではないです。

なので、今回の趣旨はポアソン方程式をPythonを使って実装しようということになります。

ポアソン方程式として有名なのは、以下の式でしょうかね。

\begin{align*}
\frac{\partial ^2 \phi}{\partial x^2} + \frac{\partial ^2 \phi}{\partial y^2} = -\frac{\rho(\boldsymbol{r})}{\epsilon}
\end{align*}

これは、電磁気で習う点電荷まわりのスカラー場を求める偏微分方程式と全く同じですね。

圧力$p$を電場$E=-\nabla \phi$のスカラー場$\phi$と置き換えると良いでしょう。

ポアソン方程式を離散化

ポアソン方程式は偏微分の方程式なので、離散化して代数的に取り扱えるようにしないと数値計算ができません。

離散化の方法としては、空間に対しては2次精度中心差分にして$i,j$での$p_{i,j}$として以下のようにして求めます。

\begin{align*}
\frac{p_{i+1,j}^{n}-2p_{i,j}^{n}+p_{i-1,j}^{n}}{\Delta x^2}+\frac{p_{i,j+1}^{n}-2 p_{i,j}^{n}+p_{i,j-1}^{n}}{\Delta y^2}=b_{i,j}^{n}\tag{4}
\end{align*}

よって、

\begin{align*}
p_{i,j}^{n}=\frac{(p_{i+1,j}^{n}+p_{i-1,j}^{n})\Delta y^2+(p_{i,j+1}^{n}+p_{i,j-1}^{n})\Delta x^2-b_{i,j}^{n}\Delta x^2\Delta y^2}{2(\Delta x^2+\Delta y^2)}\tag{5}
\end{align*}

特に、$\Delta x= \Delta y$の場合は、

\begin{align*}
p_{i,j}^{n}=\frac{p_{i+1,j}^{n}+p_{i-1,j}^{n}+p_{i,j+1}^{n}+p_{i,j-1}^{n}-b_{i,j}^{n}\Delta x^2}{4}\tag{6}
\end{align*}

例えば(1,1)の格子点上の圧力は上記のように以下のようになります。

青文字は境界条件で値が与えられているので、既知の値です。
これを見ればわかるように、(i,j)の格子点での圧力はその周りの圧力の平均値として計算しています。

境界条件を設定

ポアソン方程式は、時間発展の方程式ではなく、境界条件によって空間分布の解が与えられる方程式なので、境界条件を設定する必要があります。

今回与える境界条件は以下のようにします。

$x=0, 2$：$p=0$
$y=0, 1$：$p=0$

次に$b(\boldsymbol{r})$の条件を設定しましょう。

$i=\frac{1}{4}nx, j=\frac{1}{4}ny$：$b_{i,j}=100$
$i=\frac{3}{4}nx, j=\frac{3}{4}ny$：$b_{i,j}=-100$

ちょっとわかりにくいかもしれないので、絵にしてみました。

では、以下で実際にポアソン方程式をGoogle Colaboratoryを使ってPythonのコードを書きながら理解を深めていきたいと思います。

ポアソン方程式をPythonで実装

ひとつひとつ解説をしながらコードを書いていきます。

※ポアソン方程式は、前回の☟こちらの記事で書いているので是非参考にしてみてください。

前回の記事はこちら

【第14回Python流体の数値計算】2次元ラプラス方程式をPythonで実装する。

【ポイント】

必要なライブラリをインポートする。

コードを以下のように書きます。

import numpy as np
from matplotlib import pyplot as plt
from matplotlib import cm
from mpl_toolkits.mplot3d import Axes3D

import numpy as np

from matplotlib import pyplot as plt

from matplotlib import cm

from mpl_toolkits.mplot3d import Axes3D

NumpyとMatplotlibをインポートします。
記述ミスがなければエラーなく進みます。

では、次にポアソン方程式における「初期の状態」と「境界条件」を課します。

【ポイント】

初期状態を設定
境界条件を設定

コードを以下のように書きます。

# 条件設定
nx = 50
ny = 50
nt  = 100
xmin = 0
xmax = 2
ymin = 0
ymax = 1

dx = (xmax - xmin) / (nx - 1)
dy = (ymax - ymin) / (ny - 1)

# 初期状態
p  = np.zeros((ny, nx))
pd = np.zeros((ny, nx))
b  = np.zeros((ny, nx))
x  = np.linspace(xmin, xmax, nx)
y  = np.linspace(xmin, xmax, ny)

# 点源
b[int(ny / 4), int(nx / 4)]  = 100
b[int(3 * ny / 4), int(3 * nx / 4)] = -100

# 条件設定

nx = 50

ny = 50

nt = 100

xmin = 0

xmax = 2

ymin = 0

ymax = 1

dx = (xmax - xmin) / (nx - 1)

dy = (ymax - ymin) / (ny - 1)

# 初期状態

p = np.zeros((ny, nx))

pd = np.zeros((ny, nx))

b = np.zeros((ny, nx))

x = np.linspace(xmin, xmax, nx)

y = np.linspace(xmin, xmax, ny)

# 点源

b[int(ny / 4), int(nx / 4)] = 100

b[int(3 * ny / 4), int(3 * nx / 4)] = -100

これで、初期状態の圧力のプロファイルができています。
初期状態と言っても、境界条件を課しただけですので初期の状態には何の意味もありません。

あくまで境界条件をもとにラプラス方程式を解いて得られた圧力分布に意味があるので、初期状態はその最終状態を得るための種を与えているにすぎません。

では、いよいよ上記で設定した境界条件をもとにラプラス方程式を解いてみましょう！

【ポイント】

\begin{align*}
p_{i,j}^{n}=\frac{p_{i+1,j}^{n}+p_{i-1,j}^{n}+p_{i,j+1}^{n}+p_{i,j-1}^{n}-b_{i,j}^{n}\Delta x^2}{4}\tag{6}
\end{align*}

コードを以下のように書きます。

for it in range(nt):

    pd = p.copy()

    p[1:-1,1:-1] = (((pd[1:-1, 2:] + pd[1:-1, :-2]) * dy**2 +
                    (pd[2:, 1:-1] + pd[:-2, 1:-1]) * dx**2 -
                    b[1:-1, 1:-1] * dx**2 * dy**2) / 
                    (2 * (dx**2 + dy**2)))

    p[0, :] = 0
    p[ny-1, :] = 0
    p[:, 0] = 0
    p[:, nx-1] = 0

for it in range(nt):

pd = p.copy()

p[1:-1,1:-1] = (((pd[1:-1, 2:] + pd[1:-1, :-2]) * dy**2 +

(pd[2:, 1:-1] + pd[:-2, 1:-1]) * dx**2 -

b[1:-1, 1:-1] * dx**2 * dy**2) /

(2 * (dx**2 + dy**2)))

p[0, :] = 0

p[ny-1, :] = 0

p[:, 0] = 0

p[:, nx-1] = 0

これで最終状態の圧力の分布が計算されました。

前回のラプラス方程式を解いた際には、全ての格子点上の値の合計と前のステップでの値の合計とを比較して、「1e-4（0.01%)」以下の誤差になれば収束したと判定していましたが、今回は収束条件を反復回数（100回）にしています。

※反復計算には「ヤコビ法」と呼ばれる方法で計算しています。

では、次に収束した圧力分布を可視化してみてみましょう。

【ポイント】

Matplotlibを使って結果を可視化する。

コードを以下のように書きます。

def plot2D(x, y, p):
    fig = plt.figure(figsize=(11, 7), dpi=100)
    ax = fig.gca(projection='3d')
    X, Y = np.meshgrid(x, y)
    surf = ax.plot_surface(X, Y, p[:], rstride=1, cstride=1, cmap=cm.viridis,linewidth=0, antialiased=False)
    ax.view_init(30, 225)
    ax.set_xlabel('$x$')
    ax.set_ylabel('$y$')
    ax.set_zlabel('$p$')

def plot2D(x, y, p):

fig = plt.figure(figsize=(11, 7), dpi=100)

ax = fig.gca(projection='3d')

X, Y = np.meshgrid(x, y)

surf = ax.plot_surface(X, Y, p[:], rstride=1, cstride=1, cmap=cm.viridis,linewidth=0, antialiased=False)

ax.view_init(30, 225)

ax.set_xlabel('$x$')

ax.set_ylabel('$y$')

ax.set_zlabel('$p$')

これで、可視化するための関数が完成したので、以下のように関数を呼び出して圧力分布を可視化してみましょう。

plot2D(x, y, p)

1	plot2D(x, y, p)

【結果】

このようになればOKです(^^)/

コンター図も作成してみる

Matplotlibの公式ドキュメントは非常にわかりやすい（たぶんPythonがわかりやすいから）ので、コンター図もサクッとさくせいできます。

#オブジェクトを生成
fig = plt.figure(figsize=(11,7), dpi=100)
fig

xmin = 0
xmax = 2
ymin = 0
ymax = 2

x  = np.linspace(xmin, xmax, nx)
y  = np.linspace(xmin, xmax, ny)

X, Y = np.meshgrid(x, y)
plt.contourf(X, Y, p, alpha=0.5, cmap=cm.viridis) 
plt.colorbar()
plt.xlabel('X')
plt.ylabel('Y')