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

こんにちは(@t_kun_kamakiri)。

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

前回の記事はこちら

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

👆こちらの記事は移流方程式ですが、正確には移流項の速度が定数である場合を行いました。
なので、前回の記事では「線形の移流方程式」を扱ったという言い方の方が正しいかもしれませんね(‘ω’)

本日は、移流項が速度に依存している「2次元の移流方程式」をPythonで実装するというのをやります。

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

カマキリ

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

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

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

こんな方が対象

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

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

今回の内容はこちら

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

2次元の移流方程式
1. 初期状態の設定
2. 境界条件の設定
2次元の移流方程式をPythonで実装する
まとめ
おすすめの書籍紹介

2次元の移流方程式

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

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

これを離散化します。

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

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

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

\begin{align*}
u_{i,j}^{n+1} = u_{i,j}^n – u_{i,j} \frac{\Delta t}{\Delta x} (u_{i,j}^n-u_{i-1,j}^n) – v_{i,j}^n \frac{\Delta t}{\Delta y} (u_{i,j}^n-u_{i,j-1}^n)\\
v_{i,j}^{n+1} = v_{i,j}^n – u_{i,j} \frac{\Delta t}{\Delta x} (v_{i,j}^n-v_{i-1,j}^n) – v_{i,j}^n \frac{\Delta t}{\Delta y} (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で実装する

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

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

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

カマキリ

以下が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 = 80
c = 1
dx = 2 / (nx - 1)
dy = 2 / (ny - 1)
alpha = 0.2
dt = alpha * (dx/c)

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

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  

v = np.ones((ny, nx))  
vn = np.ones((ny, nx))
v[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 updata(nt):
  un = u.copy()
  vn = v.copy()
  
  u[1:, 1:] = (un[1:, 1:] - 
                 (un[1:, 1:] * c * dt / dx * (un[1:, 1:] - un[1:, :-1])) -
                  vn[1:, 1:] * c * dt / dy * (un[1:, 1:] - un[:-1, 1:]))
  # v[1:, 1:] = (vn[1:, 1:] -
  #                (un[1:, 1:] * c * dt / dx * (vn[1:, 1:] - vn[1:, :-1])) -
  #                vn[1:, 1:] * c * dt / dy * (vn[1:, 1:] - vn[:-1, 1:]))
    
  u[0, :] = 1
  u[-1, :] = 1
  u[:, 0] = 1
  u[:, -1] = 1
    
  # v[0, :] = 1
  # v[-1, :] = 1
  # v[:, 0] = 1
  # v[:, -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 = 80

c = 1

dx = 2 / (nx - 1)

dy = 2 / (ny - 1)

alpha = 0.2

dt = alpha * (dx/c)

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

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

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

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

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

v[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 updata(nt):

un = u.copy()

vn = v.copy()

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

(un[1:, 1:] * c * dt / dx * (un[1:, 1:] - un[1:, :-1])) -

vn[1:, 1:] * c * dt / dy * (un[1:, 1:] - un[:-1, 1:]))

# v[1:, 1:] = (vn[1:, 1:] -

# (un[1:, 1:] * c * dt / dx * (vn[1:, 1:] - vn[1:, :-1])) -

# vn[1:, 1:] * c * dt / dy * (vn[1:, 1:] - vn[:-1, 1:]))

u[0, :] = 1

u[-1, :] = 1

u[:, 0] = 1

u[:, -1] = 1

# v[0, :] = 1

# v[-1, :] = 1

# v[:, 0] = 1

# v[:, -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, updata, interval = 100,frames = 50)
HTML(ani.to_html5_video())

from IPython.display import HTML

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

HTML(ani.to_html5_video())

【結果】

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

線形移流方程式の場合は、位置によって速度の違いがあっても一定速度で進むという一般解があるのに対して、

今回の移流方程式では、☟この絵のように、速度が速い位置では速く進むというアニメーションになっていますよね。

まとめ

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

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

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

2次元の移流方程式

初期状態の設定

境界条件の設定

2次元の移流方程式をPythonで実装する

まとめ

おすすめの書籍紹介

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

流体の勉強をしたい方

数値流体の書籍

COMMENT コメントをキャンセル

2次元の移流方程式

初期状態の設定

境界条件の設定

2次元の移流方程式をPythonで実装する

まとめ

おすすめの書籍紹介

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

流体の勉強をしたい方

数値流体の書籍

【第8.1回Python流体の数値計算】バーガース方程式を差分法で実装する。

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

【第3回Python流体の数値計算】「参照渡し」と「値渡し」の違いを理解する。

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

【第6回Python流体の数値計算】１次元移流方程式のクーラン条件を確認！数値的安定性（差分法）

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

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

【Python初心者】 関数のlamda式、ジェネレーター、デコレーターを理解しよう。

【第5回Python流体の数値計算】移流方程式をGoogle Colaboratoryでアニメーション作成。

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

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

【第8.3回Python流体の数値計算】バーガース方程式の一般解（コール・ホップ変換によって拡散方程式になる）

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

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

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

【初心者歓迎！Python流体の数値計算】「2次元キャビティ流れ」と「ポアズイユ流れ」までをまとめました！

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

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

【Python入門講座】タイタニック号のデータ分析

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

COMMENT コメントをキャンセル

【Python初心者】関数のlamda式、ジェネレーター、デコレーターを理解しよう。