Fortran

【2次元熱伝導方程式】Fortranからgnuplotでグラフの3次元アニメーション作成

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

前回下記のような記事を出しまして、

前回の記事で示したアニメーションはこんな感じ

 

カマキリ

1次元の熱伝導方程式のアニメーションを作ったのです

1次元グラフのアニメーションができるのならば、gnuplotの「splot」を使って3次元アニメーションも作れるなと思ったので、記事にします(^^)/

本記事の内容

2次元の熱伝導方程式をFortranで解いて、gnuplotで3次元アニメーションを作る(^^)/

※Fortranのコードとスクリプトは本記事にアップされています。
※FortranのコードはFortran90/95の書き方にしているつもりですが正確に書けているかは不明です。
※FortranはWSLを使ってインストールしています。
※gnuplotはこのサイトを参考にすれば良いかと思います。
僕の使用環境はWindows10でWindows Subsystem for Linux(WSL)を使って、UbuntuからFortranやgnuplotをインストールして使用しています。
カマキリ

説明がくどいな~って思った方は「まとめ」へGo!!

スポンサーリンク

2次元熱伝導方程式をFortranで解く

まずは解くべき方程式を明らかにしておきましょう。

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\)方向の両端は温度0(K)の条件となっています。
また、\(x\),\(y\)の端は初期状態から変化させないような境界条件にしているので、時間発展させると中央の高い温度の部分から低い温度へ熱移動が生じて、下記のような変化を示します。
今回は、これをgnuplotを使ってアニメーションを作ろうということになります。

全体の流れ

アニメーション作成までの全体の流れを押さえておくと少し気持ちのハードルは下がるかもしれないので、先に全体の流れを示しておきます。

  1. 2次元熱伝達方程式は偏微分方程式であるため、それを数値的に解くためには離散化してあげる必要があります。
    今回は、「時間に対しては陽解法」「空間に対しては2次精度」で解くことを考えています。
  2. 次にForranを使って離散化した方程式を数値的に解くことを考えます。
    (※コードは後ほど)
  3. Fortranで出力したデータをgnuplotで順番に読み込んでいくためのスクリプトをFortranを使って作成します。

これらのことができれば、あとはgnuplotで「3.」で作成されたスクリプトを実行するだけアニメーションを作成することができます。

※1次元熱伝導方程式を題材にしたのは、簡単な方程式の方がアニメーション作成の方法を理解しやすいと思ったからです。
ゆえに、離散化も特にこだわってはいません。

2次元熱伝導方程式の離散化

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*}
この偏微分方程式を離散化するにあたって、
  • 時間微分の項(左辺)は陽解法
  • 空間微分の項(右辺)は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}\end{align*}
では、これをFortranで解きましょう(^^)/

Fortranで解く

離散化された方程式をFortranで解くためのプログラムを書いていきましょう。

以下がプログラムです。

↓前回の記事の次元を1次元から2次元に拡張しただけです。

↑今回はgnuplotの3次元表示のsplotでの表示の仕方をpm3d(カラーマップ表示)でしたかったのですが、上記のようにプロットするデータは同じ数字が続いている場合には改行をいれないといけないようでして、改行を入れるようにプログラムを書いています。
↑ちなみに「data」という名前のディレクトリを自分で作っておく必要があります。
以上がプログラムの中身です。
ここまでできたら、以下のコマンドで計算を実行させます。

実行ができれば「data」という名前のディレクトリの中に、以下のように大量のデータが保存されていることが確認できます。

これらのデータが指定したタイムステップで出力された温度の空間分布のデータです。

あとは、これらのデータを順番にgnuplotで読み込めばアニメーションが完成ということになります。

そのためには、gnuplotのスクリプトが必要になります。

スクリプトは書くのが少々面倒なので、スクリプト自体をFortranで書くことにします。

アニメーション用のgnuplotのスクリプトをFortranで作成

gnuplotのスクリプトを作成するためのFortranのプログラムがこちらです。

↑簡単に説明すると「do文」のところで先ほどの「data」という名前のディレクトリ内の大量のデータを順番に読み込む作業を行っています。

では、以下のコマンドを実行してみます。

「./a.out >movie」とするとmovieというファイルが出来上がり、その中にFortranの「print文」の内容を全て書き出してくれます。

この「movie」というファイルがgnuplotのスクリプトとなります。

以上のことが終えれば、あとはgnuplotで「movie」という名前のスクリプトを走らせるだけです。

アニメーションを作成してみよう(^^)/

「movie」というスクリプトは2通りの方法で実行できます。

方法1

↑このようにコマンドを打っても良いですし、

方法2

↑このように「gnuplot」と打って、「load “movie”」とスクリプトファイルをloadしても良いです。

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

 

カマキリ

できた

まとめ

ごちゃごちゃ説明してしまいましたが、下記のzipファイルを解凍して下記のコマンドを打ってお試しください。

20190429_thermal2D.zip

数値計算のためのFortran90/95プログラミング入門(第2版)

数値計算のためのFortran90/95プログラミング入門(第2版)

牛島 省
3,465円(12/07 01:49時点)
Amazonの情報を掲載しています
数値計算のためのFortran90/95プログラミング入門(第2版)・アンサーブック: 演習問題の解答と解説

数値計算のためのFortran90/95プログラミング入門(第2版)・アンサーブック: 演習問題の解答と解説

牛島 省
2,514円(12/06 17:15時点)
発売日: 2022/11/07
Amazonの情報を掲載しています
みんなのFortran―基礎から発展まで―

みんなのFortran―基礎から発展まで―

松本 敏郎, 野老山 貴行
3,520円(12/06 15:53時点)
Amazonの情報を掲載しています
Fortran ハンドブック

Fortran ハンドブック

田口 俊弘
3,520円(12/06 05:49時点)
Amazonの情報を掲載しています

COMMENT

メールアドレスが公開されることはありません。 が付いている欄は必須項目です