プログラミング言語

1次元の移流方程式を前進差分、後退差分、中心差分で解く。【Fortranコードあり】

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

 

本記事の内容はこちら

本記事では1次元の移流方程式を数値計算で解いた時、空間微分を以下の3パターンで解いてみたいと思います。

  • 前進差分(1次精度)で解く
  • 後退差分(1次精度)で解く
  • 中心差分(2次精度)で解く

果たしてどれがまともな数値解を得られるのか・・・

1次元の移流方程式
 
\begin{align*}\frac{\partial T(x,t)}{\partial t}+u\frac{\partial T(x,t)}{\partial x}=0\tag{1}\end{align*}
数値計算にはFortranを使います。
グラフのアニメーションはgnuplotを使います。
 

 

カマキリ

一刻も早く結果を確認したい人は「まとめ」へGo!!

※プログラムのzipファイルがあります。

 

1次元の移流方程式とは何か

 

1次元の移流方程式
 
\begin{align*}\frac{\partial T(x,t)}{\partial t}+u\frac{\partial T(x,t)}{\partial x}=0\tag{1}\end{align*}
 
※\(u\):一定速度
 
1次元の移流方程式は、下記の絵のように何かの物理量を形を変えずに一定速度\(u\)で運ぶ方程式を意味しています。
 
横軸は空間\(x\)ですが、縦軸は「密度」や「温度」だったり色々考えることができます。
本記事では、温度\(T\)として扱うことに利ます。
 
この移流方程式は、
\begin{align*}\frac{D }{D t}=\frac{\partial }{\partial t}+u\frac{\partial }{\partial x}\tag{2}\end{align*}
とすることで、物質微分を表していることがわかります。
 
つまり、流れ(流速\(u\))に沿ってどのように物理量が運ばれるかを表している偏微分方程式です。
 
(1)式を、
\begin{align*}\frac{D T(x,t)}{D t}=0\tag{3}\end{align*}
と表現することもできるので、「なるほど、温度分布が流れに沿って形を変えずに運ばれるのかー」っていうのが見えてくるわけです。
 
 
もっと実用的な偏微分方程式の中でこのような移流方程式が出てくるのは、流体の基礎方程式ですかね(‘ω’)
 
 
「流体の基礎方程式」の記事はこちら
 
 
例えば、1次元移流方程式を密度の方程式と考えるならば、
\begin{align*}\frac{D \rho(x,t)}{D t}=0\end{align*}
となります。
 
これは、「流れに沿って密度が変化しない」・・・・つまり、流体の非圧縮性仮定というなります。
 
 

1次元の移流方程式の一般解

 
 
1次元の移流方程式は解析解が存在します。
 
解析解がどのような形になっているかを知っていることで、この後数値計算したときの結果の正しさを評価することができます。
 
 
1次元の移流方程式の一般解はとっても簡単です。
 

 
以下が(1)式の一般解となります。
 
もし、速度\(u\)が正なら
\begin{align*}T(x,t)=f(x-ut)\tag{4}\end{align*}
速度\(u\)で右へ進む解
 
もし、速度\(u\)が負なら
\begin{align*}T(x,t)=g(x+ut)\tag{5}\end{align*}
速度\(u\)で左へ進む解
 

 
今回は、速度\(u>0\)としているので(4)式の形が1次元の移流方程式の解であります。
 
 
カマキリ

本当かな?と思ったら(4)式を(1)式に代入してみよう

 
(4)式を(1)式に代入して左辺が0になることを確認してみます。
 
 
\begin{align*}\frac{\partial T}{\partial (x-ut)}(-u)+u\frac{\partial T}{\partial (x-ut)}=0\tag{6}\end{align*}
やっぱり0になりましたね(^^)/
 
 

1次元の移流方程式をFortranで解く

 

1次元の移流方程式の厳密解がどのような振る舞いをするのかを理解しました。

1次元の移流方程式:「物理量の分布(形)を変えずに速度\(u\)で運ぶ」

これをFortranを使って数値計算をしてみたいと思います(^^)/

実施する数値計算の空間差分は以下の3パターンです。

  • 前進差分で解く
  • 後退差分で解く
  • 中心差分で解く
 
プログラムの全体は以下です。(例:前進差分)
 
 
もし、後退差分や中心差分で試したい場合は、以下の行を変更するだけです。
 
 
グラフのアニメーションはgnuplotを使っています。
実行方法は以下の記事に書いてありますし、本記事の最後にも書いてあります。
 
 
 
 

前進差分

 

 
\begin{align*}\frac{T^{n+1}_{i}-T^{n}_{i}}{\Delta t}+u\frac{T^{n}_{i+1}-T^{n}_{i}}{\Delta x}=0\end{align*}
なので、\(\alpha=\frac{u\Delta t}{\Delta x}>0\)と置くと、
 
 
 
\begin{align*}T^{n+1}_{i}=T^{n}_{i}-\alpha\bigg(T^{n}_{i+1}-T^{n}_{i}\bigg)\tag{7}\end{align*}
となり、(7)式を数値計算で解くことになります。
 
これをもう少し式変形変形すると、
 
\begin{align*}T^{n+1}_{i}=(1+\alpha)T^{n}_{i}-\alpha T^{n}_{i+1}\tag{8}\end{align*}
 
これを数値計算で解くと、以下のようにある位置\(i\)では\(n\)ステップ目の値\(T^{n}_{i}\)が正の値であったとしても、\(n+1\)ステップ目の値\(T^{n+1}_{i}\)は負の値を取ることになります。
 
 
 
そして、\(T^{n+1}_{i}\)(負の値)を使って次のステップ\(T^{n+2}_{i}\)を求めるので、より大きく負の値に振り切れたり、場合によっては正の値をとったり・・・・・振動するという変な数値解が得られることが予想されます。
 
つまり、全然安定に解けていないということです。
 
実際シミュレーションしてみるとこのような結果になります。
 

  1. 20190501_1


ぐっちゃぐちゃになってしまいました。
 
 

後退差分

 

つぎは後退差分でやってみましょう。

 

\begin{align*}\frac{T^{n+1}_{i}-T^{n}_{i}}{\Delta t}+u\frac{T^{n}_{i}-T^{n}_{i-1}}{\Delta x}=0\end{align*}
 
なので、\(\alpha=\frac{u\Delta t}{\Delta x}>0\)と置くと、
 
 
 
\begin{align*}T^{n+1}_{i}=T^{n}_{i}-\alpha\bigg(T^{n}_{i}-T^{n}_{i-1}\bigg)\tag{9}\end{align*}
となり、(9)式を数値計算で解くことになります。
 
これをもう少し式変形変形すると、
 
\begin{align*}T^{n+1}_{i}=(1-\alpha)T^{n}_{i}+\alpha T^{n}_{i}\tag{10}\end{align*}
 
となります。
 
これは、\(\alpha=\frac{u\Delta t}{\Delta x}<1\)としておけば、(10)式の右辺は必ず正の値を取っているので前進差分のような変な振動は起こらなさそうであることがわかります。
 
実際シミュレーションしてみるとこのような結果になります。
 
前進差分のような変な振動はなくなりましたが、形を変えずに速度\(u\)で移動するという数値計算結果にはなっていないようです。
 
 
これは、(11)式のように
 
\begin{align*}T^{n+1}_{i}=(1-\alpha)T^{n}_{i}+\alpha T^{n}_{i}\tag{11}\end{align*}
 
\((1-\alpha)T^{n}_{i}+\alpha T^{n}_{i}\)の部分が、\(i+1\)と\(i\)の割合で\(T^{n+1}_{i}\)を求めているため、角の部分がならされてしまうために起こっていると考えられます。
 
空間の離散化をもっと細かくすれば改善されるかもしれませんが、計算時間がかかることになりますし、空間の刻みを細かくしただけでは根本的な解決には至らないと思われるので別の数値解法が必要であるかもしれません。
 
 

中心差分

 

次は中心差分で解いてみましょう。

 

\begin{align*}\frac{T^{n+1}_{i}-T^{n}_{i}}{\Delta t}+u\frac{T^{n}_{i+1}-T^{n}_{i-1}}{2\Delta x}=0\end{align*}
 
なので、\(\alpha=\frac{u\Delta t}{\Delta x}>0\)と置くと、
 
 
 
\begin{align*}T^{n+1}_{i}=T^{n}_{i}-\frac{1}{2}\alpha\bigg(T^{n}_{i+1}-T^{n}_{i-1}\bigg)\tag{12}\end{align*}
となり、(12)式を数値計算で解くことになります。
 
これをもう少し式変形変形すると、
 
\begin{align*}T^{n+1}_{i}=-\frac{\alpha}{2}T^{n}_{i+1}+T^{n}_{i}+\frac{\alpha}{2}T^{n}_{i-1}\tag{13}\end{align*}
 
となります。
 
また右辺に負の項が存在することになります。
 
よって、数値解は振動するような振る舞いをすることが予想されます。
 
実際シミュレーションしてみるとこのような結果になります。
 

  1. 20190501_3


 
前進差分はぐっちゃぐちゃになりましたが、中心差分での結果は振動しながら右へ進んでいっているのがわかります。
 
 

まとめ

 

1次元の移流方程式を、

  • 前進差分
  • 後退差分
  • 中心差分

で解いてみました。

 

結果は以下のような感じです。

 

後退差分での結果が一番良い結果ということになりました。

これは速度の空間微分を流れの風上側を使って解くため、風上差分と呼ばれる差分法です。

速度\(u\)の値が正の値を取っていたので、「後退差分の結果がまとも」という結果になりましたが、

速度の値が負の値を取っている場合は「前進差分の結果がまとも」ということになります。

流れの向きを考えて差分法を考えないといけないという良い教訓になりました。

しかし、必ず流れが一方向だけという現象はめったにないので、どの方向から流れあったとしても安定的に解いてくれる数値解法を考える方が賢明かと思います。

空間の2次精度である中心差分よりも、1次精度の後退差分の方が比較的良い結果を得られているのは面白いです。

空間の2次精度とはテーラー展開の何次の項まで残すかに依るものですが、とにかく高次の項まで残せば良いというものではなく、解くべき方程式の特徴を掴んで数値計算することが大事であるということを教えてくれる結果となりました。

 

Fortranソースコードのzipファイルと実行の方法

 

本記事で紹介した計算結果は、下記よりzipファイルを入手して解凍することで再現できます。

 
そして、以下のようにコマンドを打てばFortranが実行され、gnuplotでのアニメーションを確認することができます。
 
数値計算のためのFortran90/95プログラミング入門(第2版)

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

牛島 省
3,520円(09/20 16:22時点)
発売日: 2020/01/28
Amazonの情報を掲載しています
Fortran ハンドブック

Fortran ハンドブック

田口俊弘
3,344円(09/21 03:35時点)
発売日: 2015/07/22
Amazonの情報を掲載しています

【プロフィール】

カマキリ
(^^)

大学の専攻は物性理論で、Fortranを使って数値計算をしていました。
CAEを用いた流体解析は興味がありOpenFOAMを使って勉強しています。

プロフィール記事はこちら

 

大学学部レベルの物理の解説をします 大学初学者で物理にお困りの方にわかりやすく解説します。

このブログでは主に大学以上の物理を勉強して記事にわかりやすくまとめていきます。

  • ・解析力学
  • ・流体力学
  • ・熱力学
  • ・量子統計
  • ・CAE解析(流体解析)
  • note
    noteで内容は主に「プログラミング言語」の勉強の進捗を日々書いています。また、「現在勉強中の内容」「日々思ったこと」も日記代わりに書き記しています。
  • youtube
    youtubeではオープンソースの流体解析、構造解析、1DCAEの操作方法などを動画にしています。
    (音声はありません_(._.)_)
  • Qiita
    Qiitaではプログラミング言語の基本的な内容をまとめています。
関連記事もどうぞ