【第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の完全初心者は書籍で学ぶとよい

流体の勉強をしたい方

数値流体の書籍

COMMENT コメントをキャンセル

2次元の拡散方程式

初期状態の設定

境界条件の設定

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

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

まとめ

おすすめの書籍紹介

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

流体の勉強をしたい方

数値流体の書籍

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

【第1回Python流体の数値計算】目標と環境構築。まずはここからはじめよう！

【数値流体(1次元移流方程式)】ExcelからPythonコードを実行する。

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

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

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

【Python初心者】クラスの継承の基本定義とコードの書き方を解説する。

【第2回Python流体の数値計算】Numpyの使い方を確認する。

【Python初心者】クラスの基本定義とコードの書き方を解説する。

Pythonの環境構築のためのインストール先まとめ

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

【Python初心者】数値型の3つの種類を覚えておこう

【Python初心者】if文を使って条件分岐をしてみよう！

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

【Python初心者】リスト型を学ぼう！

【第8.2回Python流体の数値計算】バーガース方程式の非線形項は異なる波数間での相互作用が生じるという意味。

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

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

Python(gspred)でスプレッドシートを操作する。

【Python初心者】Numpyを使って行列計算の使い方をまとめる。

COMMENT コメントをキャンセル