本記事の内容はこちら本記事では1次元の移流方程式を数値計算で解いた時、空間微分を以下の3パターンで解いてみたいと思います。x}=0\tag{1}\end{align*}数値計算にはFortranを使います。グラフのアニメーションはgnuplotを使います。カマキリ一刻も早く結果を確認したい人は「まとめ」へGo!!※プログラムのzipファイルがあります。x}=0\tag{1}\end{align*}※\(u\):一定速度1次元の移流方程式は、下記の絵のように何かの物理量を形を変えずに一定速度\(u\)で運ぶ方程式を意味しています。横軸は空間\(x\)ですが、縦軸は「密度」や「温度」だったり色々考えることができます。本記事では、温度\(T\)として扱うことに利ます。x}\tag{2}\end{align*}とすることで、物質微分を表していることがわかります。流れ(流速\(u\))に沿ってどのように物理量が運ばれるかを表している偏微分方程式です。t}=0\tag{3}\end{align*}と表現することもできるので、「なるほど、温度分布が流れに沿って形を変えずに運ばれるのかー」っていうのが見えてくるわけです。t}=0\end{align*}となります。これは、「流れに沿って密度が変化しない」・・・・つまり、流体の非圧縮性仮定というなります。1次元の移流方程式の一般解1次元の移流方程式は解析解が存在します。解析解がどのような形になっているかを知っていることで、この後数値計算したときの結果の正しさを評価することができます。1次元の移流方程式の一般解はとっても簡単です。以下が(1)式の一般解となります。今回は、速度\(u>0\)としているので(4)式の形が1次元の移流方程式の解であります。カマキリ本当かな?と思ったら(4)式を(1)式に代入してみよう(4)式を(1)式に代入して左辺が0になることを確認してみます。(x-ut)}=0\tag{6}\end{align*}やっぱり0になりましたね(^^)/1次元の移流方程式をFortranで解く1次元の移流方程式の厳密解がどのような振る舞いをするのかを理解しました。1次元の移流方程式:「物理量の分布(形)を変えずに速度\(u\)で運ぶ」これをFortranを使って数値計算をしてみたいと思います(^^)/実施する数値計算の空間差分は以下の3パターンです。前進差分で解く後退差分で解く中心差分で解くプログラムの全体は以下です。thermal_explicitもし、後退差分や中心差分で試したい場合は、以下の行を変更するだけです。グラフのアニメーションはgnuplotを使っています。実行方法は以下の記事に書いてありますし、本記事の最後にも書いてあります。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)式を数値計算で解くことになります。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}\)を求めるので、より大きく負の値に振り切れたり、場合によっては正の値をとったり・・・・・振動するという変な数値解が得られることが予想されます。安定に解けていないということです。実際シミュレーションしてみるとこのような結果になります。20190501_1ぐっちゃぐちゃになってしまいました。後退差分つぎは後退差分でやってみましょう。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)式を数値計算で解くことになります。T^{n}_{i}\tag{10}\end{align*}となります。x}<1\)としておけば、(10)式の右辺は必ず正の値を取っているので前進差分のような変な振動は起こらなさそうであることがわかります。実際シミュレーションしてみるとこのような結果になります。https://takun-physics.net/wp-content/uploads/2019/05/20190501_2.mp4前進差分のような変な振動はなくなりましたが、形を変えずに速度\(u\)で移動するという数値計算結果にはなっていないようです。T^{n}_{i}\)の部分が、\(i+1\)と\(i\)の割合で\(T^{n+1}_{i}\)を求めているため、角の部分がならされてしまうために起こっていると考えられます。空間の離散化をもっと細かくすれば改善されるかもしれませんが、計算時間がかかることになりますし、空間の刻みを細かくしただけでは根本的な解決には至らないと思われるので別の数値解法が必要であるかもしれません。次は中心差分で解いてみましょう。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*}となります。右辺に負の項が存在することになります。よって、数値解は振動するような振る舞いをすることが予想されます。実際シミュレーションしてみるとこのような結果になります。20190501_3前進差分はぐっちゃぐちゃになりましたが、中心差分での結果は振動しながら右へ進んでいっているのがわかります。結果は以下のような感じです。後退差分での結果が一番良い結果ということになりました。これは速度の空間微分を流れの風上側を使って解くため、風上差分と呼ばれる差分法です。速度の値が負の値を取っている場合は「前進差分の結果がまとも」ということになります。流れの向きを考えて差分法を考えないといけないという良い教訓になりました。流れが一方向だけという現象はめったにないので、どの方向から流れあったとしても安定的に解いてくれる数値解法を考える方が賢明かと思います。空間の2次精度である中心差分よりも、1次精度の後退差分の方が比較的良い結果を得られているのは面白いです。空間の2次精度とはテーラー展開の何次の項まで残すかに依るものですが、とにかく高次の項まで残せば良いというものではなく、解くべき方程式の特徴を掴んで数値計算することが大事であるということを教えてくれる結果となりました。本記事で紹介した計算結果は、下記よりzipファイルを入手して解凍することで再現できます。前進差分:20190430_iryu1.zip後退差分:20190430_iryu2.zip中心差分:20190430_iryu3.zipそして、以下のようにコマンドを打てばFortranが実行され、gnuplotでのアニメーションを確認することができます。【1階微分!中心差分の4次精度】テーラー展開から真面目に導出してみた。【1階微分!中心差分の4次精度】テーラー展開から真面目に導出してみた。-【CAE解析の目的とは?】CAEは実機の代替えなのか?精度を無視してパラメータスタディすることは意味がないのか?【Udemyがお勧め理由5つ】プログラミングを勉強するなら動画学習おすすめ【まとめ】比熱の温度依存性について。古典論と量子論から理解する。【理想気体の比熱の温度依存性(7)】シュレディンガー方程式の極座標表記(2原子分子)から動径方向だけを抽出しようコメントを残すコメントをキャンセルプロフィールカマキリ(^^)大学以上の物理を勉強してブログに書いていきます。大学の専攻は物性理論で、Fortranを使って数値計算をしていました。CAE解析(流体解析)・毎日note更新しています。内容は、「現在勉強中の内容」「日々思ったこと」を書いています。1DCAEの操作方法などを動画にしています。youtube・Qiitaではプログラミング言語の基本的な内容をまとめています。(10)ブログをメールで購読メールアドレスを記入して購読すれば、更新をメールで受信できます。