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

こんにちは（@t_kun_kamakiri）(‘◇’)ゞ

↓前回の記事で2次元熱伝導方程式をFortranで解いて、アニメーションをgnuplotで作成するというのをしました。

前回の記事はこちら

Fortranからgnuplotを使って3次元アニメーション作成

アニメーションはこんな感じ

動画プレーヤー

Media error: Format(s) not supported or source(s) not found

ファイルをダウンロード: https://takun-physics.net/wp-content/uploads/2019/04/20190429_192249.mp4?_=1

00:00

ボリューム調節には上下矢印キーを使ってください。

さらに、↓vtk形式としてparaviewでアニメーション作成するというのもしました。

前回の記事はこちら

Fortranで解いた問題をParaviewを使ってvtkでアニメーション作成

↓paraviewによる基本的な操作は動画にしています。

動画プレーヤー

Media error: Format(s) not supported or source(s) not found

ファイルをダウンロード: https://takun-physics.net/wp-content/uploads/2019/04/20190430_143929.mp4?_=2

00:00

ボリューム調節には上下矢印キーを使ってください。

カマキリ

「gnuplotのスクリプトを書く」のも「vtk形式でparaviewで読む」のも結構調べるの大変だった

上記のように色々試しましたが、結局アニメーション作成はPythonでさくっとやってしまおうと思いましたね・・・・(‘ω’)

本記事の内容

Pythonで2次元の熱伝導方程式を解いて、Pythonでアニメーションまで作成。

※Pythonの動作にはJupyter Notebookを使用しています。

Pythonの使用環境
2次元熱伝導方程式を離散化して解く
Pythonで解いて2次元のアニメーション
アニメーション作成の結果
参考書

Pythonの使用環境

Pythonの使用環境は色々あると思います。

以下に代表的なものを紹介しておきます。

①Pythonの本家をインストール

②AnacondaをインストールしてJupyter notebookを使用

③PycharmをインストールしてPythonを使用

④visual studioをインストールしてPythonを使用

カマキリ

僕はAnacondaをインストールしてJupyter notebookを使用したいと思います。

2次元熱伝導方程式を離散化して解く

解きたい方程式は、下記の記事でも詳しく解説しています。

前回の記事はこちら

Fortranからgnuplotを使って3次元アニメーション作成

Pythonで解きたい方程式は↓こちらです。

2次元熱伝導方程式

$\begin{align*}\frac{\partial T}{\partial t}=\kappa\bigg(\frac{\partial^2T}{\partial x^2}+\frac{\partial^2T}{\partial y^2}\bigg)\tag{1}\end{align*}$

熱伝導方程式は、初期状態を与えれば(1)式に従って時間発展をする方程式です。

今回は初期状態を以下のようにしておきます。

中央だけ中央だけ温度が高く、

$x$ 方向の両端は温度０(K)の条件となっています。

また、

$x$ ,

$y$ の端は初期状態から変化させないような境界条件にしているので、時間発展させると中央の高い温度の部分から低い温度へ熱移動が生じて、下記のような変化を示します。

立体的に見るとこんな感じですね(^^)/

(1)式の偏微分方程式を離散化するにあたって、

時間微分の項（左辺）は陽解法
空間微分の項（右辺）は2次精度

で離散化で進めていきましょう。

離散かされた式が↓こちらです。

$\begin{align*}T^{n+1}(x_{i},y_{j})=T^{n}(x_{i},y_{j})+\kappa dt\frac{T^{n}(x_{i}+\Delta x,y_{j})-2T^{n}(x_{i},y_{j})-T^{n}(x_{i}-\Delta x,y_{j})}{\Delta x^2}\\ \kappa dt\frac{T^{n}(x_{i},y_{j}+\Delta y)-2T^{n}(x_{i},y_{j})-T^{n}(x_{i},y_{j}-\Delta y)}{\Delta y^2}\tag{2}\end{align*}$

今回は、これをPythonで解いてアニメーション作成までしようということになります。

Pythonで解いて2次元のアニメーション

Pythonのコードを載せておきましょう。

%matplotlib nbagg
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation


imax=100
jmax=imax
ntmax=2500
dt=0.004
kp=1
xmax=10
ymax=10
dx=xmax/imax
dy=ymax/jmax

    
X = np.linspace(0, xmax, imax)
Y = np.linspace(0, ymax, jmax)
temp=np.zeros([imax,jmax])

x2,y2=np.meshgrid(X,Y)

#initial state
for j in range(jmax):
    for i in range(imax):
        x=dx*i
        y=dy*j
        if (i<imax/2):
            temp[j,i]=x
        else:
            temp[j,i]=xmax-x


plt.pcolor(x2,y2,temp,cmap ='hsv')
plt.colorbar()

fig = plt.figure()
ims=[]
for t in range(0,ntmax):
    for j in range(1,jmax-1):
        for i in range(1,imax-1):
            temp[j,i]=temp[j,i]+kp*dt/(dx*dx)*(temp[j,i+1]-2*temp[j,i]+temp[j,i-1])+kp*dt/(dy*dy)*(temp[j+1,i]-2*temp[j,i]+temp[j-1,i])
    if (t%100==0):
        im = plt.imshow(temp,animated=True,cmap ='seismic')
        ims.append([im])


plt.colorbar()
plt.xlabel("x")
plt.ylabel("y")

ani = animation.ArtistAnimation(fig, ims, interval=50, blit=True,repeat_delay=10)

ani.save('anim.gif', writer="imagemagick")
ani.save('anim.mp4', writer="ffmpeg")
plt.show()