【第12回Python流体の数値計算】2次元拡散方程式をGoogle Colabでアニメーション作成する。

こんにちは(@t_kun_kamakiri)。

前回の記事では、「2次元の移流方程式」をGoogle Colabでアニメーション作成を行いました。

前回の記事はこちら

【第11回Python流体の数値計算】2次元移流方程式をGoogle Colabでアニメーション作成する。

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

今日作成する動画は以下のような感じになります。

１次元の拡散方程式については、以前の記事でも書きましたが今回はその２次元バージョンというわけです。

カマキリ

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

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

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

こんな方を対象

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

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

今回の内容はこちら

２次元の拡散方程式の数値計算
２次元のアニメーション作成

2次元の拡散方程式
1. 初期状態の設定
2. 境界条件の設定
2次元の拡散方程式をPythonで実装する
どうせならアニメーションにしたい
まとめ
おすすめの書籍紹介

2次元の拡散方程式

まずは、解くべき方程式の確認を行いましょう。

\begin{align*}
\frac{\partial u}{\partial t} = \nu \frac{\partial ^2 u}{\partial x^2} + \nu \frac{\partial ^2 u}{\partial y^2}\\
\frac{\partial v}{\partial t} = \nu \frac{\partial ^2 v}{\partial x^2} + \nu \frac{\partial ^2 v}{\partial y^2}\tag{1}
\end{align*}

これを数値計算できるように空間微分に対して2次精度の中心差分を施して以下のように離散化します。

\begin{align*}
\frac{u_{i,j}^{n+1} – u_{i,j}^n}{\Delta t} = \nu \frac{u_{i+1,j}^n – 2 u_{i,j}^n + u_{i-1,j}^n}{\Delta x^2} + \nu \frac{u_{i,j+1}^n-2 u_{i,j}^n + u_{i,j-1}^n}{\Delta y^2}\\
\frac{v_{i,j}^{n+1} – v_{i,j}^n}{\Delta t} = \nu \frac{v_{i+1,j}^n – 2 v_{i,j}^n + v_{i-1,j}^n}{\Delta x^2} + \nu \frac{v_{i,j+1}^n-2 v_{i,j}^n + v_{i,j-1}^n}{\Delta y^2}\tag{2}
\end{align*}

空間微分に対しては、数値計算の安定性のために後退差分で行います。

今、ほしい情報はn+1ステップでの$u^{n+1}$や$v^{n+1}$の値ですので、式変形して、

\begin{align*}
u_{i,j}^{n+1} = u_{i,j}^n &+ \frac{\nu \Delta t}{\Delta x^2}(u_{i+1,j}^n – 2 u_{i,j}^n + u_{i-1,j}^n) \\
&+ \frac{\nu \Delta t}{\Delta y^2}(u_{i,j+1}^n-2 u_{i,j}^n + u_{i,j-1}^n)
\\
v_{i,j}^{n+1} = v_{i,j}^n &+ \frac{\nu \Delta t}{\Delta x^2}(v_{i+1,j}^n – 2 v_{i,j}^n + v_{i-1,j}^n) \\
&+ \frac{\nu \Delta t}{\Delta y^2}(v_{i,j+1}^n-2 v_{i,j}^n + v_{i,j-1}^n)\tag{3}
\end{align*}

上記のような偏微分方程式の場合は、時々刻々と変化する物理量（今の場合は流速）に対して解くので初期状態を設定しないと解くための種がありません。
というわけで、初期状態を設定します。

初期状態の設定

初期状態を以下のようにします。

\begin{align*}
u,\ v\ = \begin{cases}\begin{matrix}
2 & \text{for } x,y \in (0.5, 1)\times(0.5,1) \cr
1 & \text{それ以外}
\end{matrix}\end{cases}\tag{4}
\end{align*}

さらに、数値計算は有限の大きさに対して解くので計算領域の端に対して数値を設けないと問題を解くことができません。

なので、以下のように境界条件を設定します。

境界条件の設定

境界条件を以下のようにします。

\begin{align*}
u = 1,\ v = 1 \text{ for } \begin{cases} \begin{matrix}x=0,2\cr y=0,2 \end{matrix}\end{cases}
\end{align*}

では、Pythonを使って2次元の移流方程式を解いてみましょう。

Pythonの使用環境はGoogle Colabとします。

2次元の拡散方程式をPythonで実装する

まずは必要なライブラリをインポートします。

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

import numpy as np

from matplotlib import pyplot as plt

from matplotlib import cm

from mpl_toolkits.mplot3d import Axes3D

import matplotlib.animation as animation

次に、「条件設定」と「初期状態設定」を行います。

#条件設定
nx = 31
ny = 31
nt = 17
nu = .05
dx = 2 / (nx - 1)
dy = 2 / (ny - 1)
sigma = .25
dt = sigma * dx * dy / nu

x = np.linspace(0, 2, nx)
y = np.linspace(0, 2, ny)

#初期状態の設定
u = np.ones((ny, nx))  
un = np.ones((ny, nx))

u[int(.5 / dy):int(1 / dy + 1),int(.5 / dx):int(1 / dx + 1)] = 2

#条件設定

nx = 31

ny = 31

nt = 17

nu = .05

dx = 2 / (nx - 1)

dy = 2 / (ny - 1)

sigma = .25

dt = sigma * dx * dy / nu

x = np.linspace(0, 2, nx)

y = np.linspace(0, 2, ny)

#初期状態の設定

u = np.ones((ny, nx))

un = np.ones((ny, nx))

u[int(.5 / dy):int(1 / dy + 1),int(.5 / dx):int(1 / dx + 1)] = 2

では、初期状態のプロファイルを確認してみましょう！
Matplotlibを使って可視化してみます。

fig = plt.figure()
ax = fig.gca(projection='3d')
X, Y = numpy.meshgrid(x, y)
surf = ax.plot_surface(X, Y, u, rstride=1, cstride=1, cmap=cm.viridis,
        linewidth=0, antialiased=False)

ax.set_xlim(0, 2)
ax.set_ylim(0, 2)
ax.set_zlim(1, 2.5)

ax.set_xlabel('$x$')
ax.set_ylabel('$y$');

fig = plt.figure()

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

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

surf = ax.plot_surface(X, Y, u, rstride=1, cstride=1, cmap=cm.viridis,

linewidth=0, antialiased=False)

ax.set_xlim(0, 2)

ax.set_ylim(0, 2)

ax.set_zlim(1, 2.5)

ax.set_xlabel('$x$')

ax.set_ylabel('$y$');

【結果】

初期状態はこのようになっています。

では、次にスライスを使って(3)式の偏微分方程式の離散化を数値計算で求めましょう。

def diffuse(nt):
    u[int(.5 / dy):int(1 / dy + 1),int(.5 / dx):int(1 / dx + 1)] = 2  
    
    for n in range(nt + 1): 
        un = u.copy()
        u[1:-1, 1:-1] = (un[1:-1,1:-1] + 
                        nu * dt / dx**2 * 
                        (un[1:-1, 2:] - 2 * un[1:-1, 1:-1] + un[1:-1, 0:-2]) +
                        nu * dt / dy**2 * 
                        (un[2:,1: -1] - 2 * un[1:-1, 1:-1] + un[0:-2, 1:-1]))
        u[0, :] = 1
        u[-1, :] = 1
        u[:, 0] = 1
        u[:, -1] = 1

    
    fig = plt.figure()
    ax = fig.gca(projection='3d')
    surf = ax.plot_surface(X, Y, u[:], rstride=1, cstride=1, cmap=cm.viridis,
        linewidth=0, antialiased=True)
    ax.set_zlim(1, 2.5)
    ax.set_xlabel('$x$')
    ax.set_ylabel('$y$')

def diffuse(nt):

u[int(.5 / dy):int(1 / dy + 1),int(.5 / dx):int(1 / dx + 1)] = 2

for n in range(nt + 1):

un = u.copy()

u[1:-1, 1:-1] = (un[1:-1,1:-1] +

nu * dt / dx**2 *

(un[1:-1, 2:] - 2 * un[1:-1, 1:-1] + un[1:-1, 0:-2]) +

nu * dt / dy**2 *

(un[2:,1: -1] - 2 * un[1:-1, 1:-1] + un[0:-2, 1:-1]))

u[0, :] = 1

u[-1, :] = 1

u[:, 0] = 1

u[:, -1] = 1

fig = plt.figure()

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

surf = ax.plot_surface(X, Y, u[:], rstride=1, cstride=1, cmap=cm.viridis,

linewidth=0, antialiased=True)

ax.set_zlim(1, 2.5)

ax.set_xlabel('$x$')

ax.set_ylabel('$y$')

計算時間短縮のためvについては解いていません。
メイン部分は関数にまとめておくことにします。
そうすることでntを引数に入れることでほしい時間帯（正確にはステップ数）でのuのプロファイルを可視化することができます
こちらの記事で扱ったようにスライスを使っています。

では、「nt=10(ステップ数が10）」の場合のuのプロファイルを見てみましょう！

diffuse(10)

1	diffuse(10)

【結果】

では、「nt=30(ステップ数が30）」の場合のuのプロファイルは？

こんな感じで拡散しているのがわかりますね(^^♪

どうせならアニメーションにしたい

先ほどのように静止画でuのプロファイルを見るよりも、やっぱりアニメーションにしたいですよね。

アニメーションの作成には以下の２つあります。

matplotlib.animation.ArtistAnimation：グラフのオブジェクトのリストを格納することでアニメーションを作成する。
matplotlib.animation.FuncAnimation：グラフ行進用の関数を用意して逐一更新しながらアニメーションを作成する。

今回は、「matplotlib.animation.FuncAnimation」と使ってアニメーションを作成しています。

カマキリ

以下がPythonのコードです！

import numpy as np
from matplotlib import pyplot as plt
from matplotlib import cm
from mpl_toolkits.mplot3d import Axes3D ##library for 3d projection plots
import matplotlib.animation as animation


#条件設定
nx = 31
ny = 31
nt = 17
nu = .05
dx = 2 / (nx - 1)
dy = 2 / (ny - 1)
sigma = .25
dt = sigma * dx * dy / nu

x = np.linspace(0, 2, nx)
y = np.linspace(0, 2, ny)

#初期状態
u = np.ones((ny, nx))  
un = np.ones((ny, nx))
u[int(.5 / dy):int(1 / dy + 1),int(.5 / dx):int(1 / dx + 1)] = 2  



fig = plt.figure()
ax = fig.gca(projection='3d')
ax.set_zlim(1, 2.5)
ax.set_xlabel('$x$')
ax.set_ylabel('$y$')

ax = fig.gca(projection='3d')
X, Y = np.meshgrid(x, y)
surf = ax.plot_surface(X, Y, u, rstride=1, cstride=1, cmap=cm.viridis,linewidth=0, antialiased=False)


def diffuse(nt):
  un = u.copy()
  u[1:-1, 1:-1] = (un[1:-1,1:-1] + 
                          nu * dt / dx**2 * 
                          (un[1:-1, 2:] - 2 * un[1:-1, 1:-1] + un[1:-1, 0:-2]) +
                          nu * dt / dy**2 * 
                          (un[2:,1: -1] - 2 * un[1:-1, 1:-1] + un[0:-2, 1:-1]))
  u[0, :] = 1
  u[-1, :] = 1
  u[:, 0] = 1
  u[:, -1] = 1

    
    
  ax.set_title('time = '+ str(round(nt*dt,2))+ 'sec, '+ str(nt)+'step')
  ax.plot_surface(X, Y, u[:], rstride=1, cstride=1, cmap=cm.viridis,linewidth=0, antialiased=True)

import numpy as np

from matplotlib import pyplot as plt

from matplotlib import cm

from mpl_toolkits.mplot3d import Axes3D ##library for 3d projection plots

import matplotlib.animation as animation

#条件設定

nx = 31

ny = 31

nt = 17

nu = .05

dx = 2 / (nx - 1)

dy = 2 / (ny - 1)

sigma = .25

dt = sigma * dx * dy / nu

x = np.linspace(0, 2, nx)

y = np.linspace(0, 2, ny)

#初期状態

u = np.ones((ny, nx))

un = np.ones((ny, nx))

u[int(.5 / dy):int(1 / dy + 1),int(.5 / dx):int(1 / dx + 1)] = 2

fig = plt.figure()

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

ax.set_zlim(1, 2.5)

ax.set_xlabel('$x$')

ax.set_ylabel('$y$')

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

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

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

def diffuse(nt):

un = u.copy()

u[1:-1, 1:-1] = (un[1:-1,1:-1] +

nu * dt / dx**2 *

(un[1:-1, 2:] - 2 * un[1:-1, 1:-1] + un[1:-1, 0:-2]) +

nu * dt / dy**2 *

(un[2:,1: -1] - 2 * un[1:-1, 1:-1] + un[0:-2, 1:-1]))

u[0, :] = 1

u[-1, :] = 1

u[:, 0] = 1

u[:, -1] = 1

ax.set_title('time = '+ str(round(nt*dt,2))+ 'sec, '+ str(nt)+'step')

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

※計算時間短縮のために、$y$方向の速度$v$については解かないようにしています。

from IPython.display import HTML
ani = animation.FuncAnimation(fig, diffuse, interval = 100,frames = 50)
HTML(ani.to_html5_video())

from IPython.display import HTML

ani = animation.FuncAnimation(fig, diffuse, interval = 100,frames = 50)

HTML(ani.to_html5_video())

【結果】

じゃっかん残像を残しながらアニメーションが進みますが・・・・アニメーションは完成です。

まとめ

今回は、2次元の拡散方程式をPythonで実装してアニメーション作成を行いました。

前回の記事の内容とあまり変わらず復習程度のないようになりましたが、本記事のシリーズを読むことでPythonを使ったナビエストークス方程式の実装まではいけると思いますので、是非最後までお付き合いくださいm(__)m

【第12回Python流体の数値計算】2次元拡散方程式をGoogle Colabでアニメーション作成する。

2次元の拡散方程式

初期状態の設定

境界条件の設定

2次元の拡散方程式をPythonで実装する

どうせならアニメーションにしたい

まとめ

おすすめの書籍紹介

Pythonの完全初心者は書籍で学ぶとよい

流体の勉強をしたい方

数値流体の書籍

【第13回Python流体の数値計算】2次元バーガース方程式をGoogle Cola...

Pythonで2次元の熱伝導方程式をアニメーションが簡単にできた

【第9回Python流体の数値計算】Numpyでスライスを使って計算を速くする！

【Pythonプログラム】衝撃波菅の理論とOpenFOAMの結果比較

【電子書籍出します】OpenFOAM（自宅で極める流体解析）2次元円柱まわりの流れ

Python ファイルの読み込み方法

【Python初心者】別ファイルの関数やクラスの実行（モジュール化）

Pythonで1次元熱伝導方程式のグラフのアニメーションが簡単にできた

【Python】ひもでつながれた2質点の放物線運動

【第5回Python流体の数値計算】移流方程式をGoogle Colaborator...

【Python入門講座】Pandasの基礎を身に付けよう

【第8.3回Python流体の数値計算】バーガース方程式の一般解（コール・ホップ変換...

モンテカルロ法を用いて円周率をPythonで近似計算する

【Python初心者】繰り返し文（for文）を学ぼう！

【第17回Python流体の数値計算】2次元ナビエストークス方程式！ポアズイユ流れを...

【第4回Python流体の数値計算】１次元移流方程式を差分法で実装する。

【Python初心者】辞書型を学ぼう！

【第11回Python流体の数値計算】2次元移流方程式をGoogle Colabでア...

【デイトラPythonコース評判と受講後】SNSマーケティングの自動化のためのコース...

【Python初心者】タプル型を学ぼう！

COMMENT コメントをキャンセル

初心者向け

中級者向け