こんにちは（@t_kun_kamakiri）(‘◇’)ゞ
↓前回の記事で2次元熱伝導方程式をFortranで解いて、アニメーションをgnuplotで作成するというのをしました。
前回の記事はこちら
アニメーションはこんな感じ
さらに、↓vtk形式としてparaviewでアニメーション作成するというのもしました。
前回の記事はこちら
↓paraviewによる基本的な操作は動画にしています。
カマキリ
「gnuplotのスクリプトを書く」のも「vtk形式でparaviewで読む」のも結構調べるの大変だった
上記のように色々試しましたが、結局アニメーション作成はPythonでさくっとやってしまおうと思いましたね・・・・(‘ω’)
本記事の内容
Pythonで2次元の熱伝導方程式を解いて、Pythonでアニメーションまで作成。
※Pythonの動作にはJupyter Notebookを使用しています。
Pythonの使用環境
Pythonの使用環境は色々あると思います。
以下に代表的なものを紹介しておきます。
②AnacondaをインストールしてJupyter notebookを使用
④visual studioをインストールしてPythonを使用
カマキリ
僕はAnacondaをインストールしてJupyter notebookを使用したいと思います。
2次元熱伝導方程式を離散化して解く
解きたい方程式は、下記の記事でも詳しく解説しています。
前回の記事はこちら
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*}
\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のコードを載せておきましょう。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 | %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() |
アニメーション作成の結果
作成されたアニメーション作成はこんな感じです。
2次元配列だとPythonは結構時間かかるんですね・・・(‘_’)
Follow @t_kun_kamakiri