【第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

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

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

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')

#オブジェクトを生成

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')

【結果】

このように圧力の分布もコンターで可視化できます。

まとめ

今回は、2次元のポアソン方程式をPythonで実装してみました。

本記事のシリーズを読むことでPythonを使ったナビエストークス方程式の実装まではいけると思いますので、是非最後までお付き合いくださいm(__)m

次回は、いよいよナビエストークス方程式をPythonで実装してみたいと思います(^^)/

次回の記事はこちら

【第16回Python流体の数値計算】2次元キャビティ流れをPythonで実装する。

ご質問に対する回答

ブログのコメントに「b[int(ny / 2), int(nx / 2)] = -10.e-8」の電荷として解けないかという質問がありましたので、全体のコードを載せておきます。
下記のコードをmain.pyに書き、python main.pyで実行すれば動作するかと思います。

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

# 条件設定
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 / 2), int(nx / 2)]  =  -10.e-8

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

	
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)
plt.show()

from matplotlib import pyplot as plt

from matplotlib import cm

from mpl_toolkits.mplot3d import Axes3D

import numpy as np

# 条件設定

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 / 2), int(nx / 2)] = -10.e-8

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

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)

plt.show()

POSTED COMMENT

牧田より:

2020年6月26日 11:35

(5)式から(6)式への式変形のところで質問があります。
Δx=Δyとした場合、(6)式のb_(i,j)^nの部分は、b_(i,j)^nではなく、(b_(i,j)^n)(Δy)^2 or (b_(i,j)^n)(Δx)^2ではないのですか？

私の計算ミスでしたらごめんなさい。

よろしくお願いします。

返信
- korokoro より:
  
  2020年6月26日 12:01
  
  ご指摘ありがとうございます。
  おっしゃるとおりです。失礼致しました。
  修正したいと思います。
  数値計算の方は、修正の必要は無く正しいです。
  
  ありがとうございました。
  
  返信
ねこより:

2022年9月26日 16:44

最近Pythonを使い始めとても参考になる内容をありがとうございます。
今回使ったコードを電位分布、電場分布を求める場合に真ん中に電荷があるときは点源のところをb[int(ny / 2), int(nx / 2)] = -10.e-8
で使えますか？

返信
ねこまるより:

2022年9月26日 16:51

最近Pythonを触り始めとても参考になる内容をありがとうございます。
初心的な質問で申し訳ありませんが
このコードを真ん中に電荷がある領域の電位分布を求める時の点源のコードはb[int(ny / 2), int(nx / 2)] = -10.e-8だとうまくいきませんでしたがどう変えたらいいでしょうか

返信
- korokoro より:
  
  2022年9月27日 15:45
  
  コメントいただきありがとうございます。
  うまくいかないというのはエラーで計算ができないということでしょうか。
  私の手元で以下のようなコードをmain.pyに書いて、python main.pyと実行するとポアソン方程式を解いてくれているように思えます。
  一度お確かめください。なお、Pythonはインデントをつけないとエラーになるのでコピペする際には十分注意してください。
  コードはブログの一番最後にも書きましたのでご確認お願い致します。
  from matplotlib import pyplot as plt
  from matplotlib import cm
  from mpl_toolkits.mplot3d import Axes3D
  import numpy as np
  
  # 条件設定
  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 / 2), int(nx / 2)] = -10.e-8
  
  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
  
  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)
  plt.show()
  
  返信
ねこまるより:

2022年9月27日 22:37

丁寧な回答ありがとうございます。とりあえず一つの点で電荷を置く場合は理解できましたが1ｍの領域の真ん中に10ｃｍの電荷を置く、などここからここまでの範囲に電荷を置く場合のコードがどうなるか分かりません。お時間があればご回答頂けないでしょうか？

返信