Python

【第7回Python流体の数値計算】1次元拡散方程式を差分法で実装する。

こんにちは(@t_kun_kamakiri)。

今回は、1次元の拡散方程式」をPythonで実装するというのをやります。

前回の記事では、クーラン条件より「1次元の移流方程式」についての数値的安定性の内容をまとめましたが、「1次元の拡散方程式」も数値的安定性の条件が存在します。
本記事では、わざと数値的に不安定な状態での数値計算もやりますので、解の振る舞いがどのようになるのかを確認しながら見ていただければと思います。

☟前回の記事がまだの人はこちらから

前回の記事はこちら
 

本件の基本的な内容はこちらのサイトにそってやっていきます。

この記事ではこんな人を対象にしています。

  • Pythonを使い始めたけどどう使うかわからない
  • 流体の数値計算をはじめて勉強する人
 
本記事シリーズの最終目標は、ナビエストークス方程式をPythonで実装することですが、いきなりナビエストークス方程式を実装するのは難しいので各項の意味を確認しながら進めていきたいと思います。
 
 
 
今回の内容はこちら

  • 一次元拡散方程式の数値計算
  • 拡散方程式の安定性
 
 
 
拡散方程式というのはナビエストークス方程式の「左辺の第一項と右辺第二項」のみを使った偏微分方程式です。
 
 
 
では、まずは「1次元拡散方程式とは何か」という話から始めたいと思います。
 
スポンサーリンク

1次元の拡散方程式とは

 

拡散方程式について簡単に解説をしておこうと思います。

1次元の拡散方程式

\begin{align*}
\frac{\partial u}{\partial t}= \nu \frac{\partial^2 u}{\partial x^2} \tag{1}
\end{align*}

これは名前の通り、初期状態の形(速度、温度など)が系の中で均等になるように拡散する偏微分方程式です。

後ほど実際にPythonのコードを示しながらアニメーション作成を行いますので、その解の振る舞いを確認して理解を深めてもらえればと思います_(._.)_

拡散方程式を離散化する

 

拡散方程式を、

  • 時間の一階微分の項
  • 空間の二階微分の項

となっています。

時間の項に対しては、

時間微分の項(第一項)

\begin{align*}
\frac{\partial u}{\partial t} \approx \frac{u(t+\Delta t)-u(t)}{\Delta t}\tag{2}
\end{align*}
とします。
空間微分の項は、2次精度中心差分で行います。
 
ある位置\(x_0+\Delta x_0\)における\(x_0\)のまわりでテーラー展開を考えます。
 
\begin{align*}
u(x_0+\Delta x) = u(x_0) +\Delta x \frac{\partial u(x)}{\partial x}\bigg|_{x_0}  + \frac{\Delta x^2}{2} \frac{\partial ^2 u(x)}{\partial x^2}\bigg|_{x_0} + \frac{\Delta x^3}{3!} \frac{\partial ^3 u(x)}{\partial x^3}\bigg|_{x_0} + O(\Delta x^4)\tag{3}
\end{align*}
 

さらに、

ある位置\(x_0-\Delta x_0\)における\(x_0\)のまわりでテーラー展開を考えます。
\begin{align*}
u(x_0-\Delta x) = u(x_0) -\Delta x \frac{\partial u(x)}{\partial x}\bigg|_{x_0}  + \frac{\Delta x^2}{2} \frac{\partial ^2 u(x)}{\partial x^2}\bigg|_{x_0} – \frac{\Delta x^3}{3!} \frac{\partial ^3 u(x)}{\partial x^3}\bigg|_{x_0} + O(\Delta x^4)\tag{4}
\end{align*}
今ほしいのは、\(x_0\)における2階微分\(\frac{\partial ^2 u(x)}{\partial x^2}\bigg|_{x_0}\)なので(3)と(4)を足すことで、
 
\begin{align*}
\frac{\partial ^2 u}{\partial x^2}=\frac{u(x+\Delta x)-2u(x)+u(x-\Delta x)}{\Delta x^2} + O(\Delta x^2)
\end{align*}
 
得ることができます。
 
数値計算では離散化をしてある位置という任意の位置を指定することはできません。
つまり、空間は連続ではなく格子点でしか表現されていません。

ある位置\(x_0\)と書かずに、ある位置(格子点)\(i\)と書くのが一般的です。
そうすると上で行った計算は以下のように書くことができます。

\begin{align*}
u_{i+1} = u_i + \Delta x \frac{\partial u}{\partial x}\bigg|_i + \frac{\Delta x^2}{2} \frac{\partial ^2 u}{\partial x^2}\bigg|_i + \frac{\Delta x^3}{3!} \frac{\partial ^3 u}{\partial x^3}\bigg|_i + O(\Delta x^4)\tag{5}
\end{align*}

\begin{align*}
u_{i-1} = u_i – \Delta x \frac{\partial u}{\partial x}\bigg|_i + \frac{\Delta x^2}{2} \frac{\partial ^2 u}{\partial x^2}\bigg|_i – \frac{\Delta x^3}{3!} \frac{\partial ^3 u}{\partial x^3}\bigg|_i + O(\Delta x^4)\tag{6}
\end{align*}

(5)と(6)を足して、
\begin{align*}
u_{i+1} + u_{i-1} = 2u_i+\Delta x^2 \frac{\partial ^2 u}{\partial x^2}\bigg|_i + O(\Delta x^4)
\end{align*}

 \(\frac{\partial ^2 u}{\partial x^2}\bigg|_i\)について少し整理をすると、
\begin{align*}\frac{\partial ^2 u}{\partial x^2}\bigg|_i=\frac{u_{i+1}-2u_{i}+u_{i-1}}{\Delta x^2} + O(\Delta x^2)\end{align*}
となります。

2次以降は無視をして計算をすることになるので、2次精度の微分ということになります。

以上より、(1)式の拡散方程式の離散化は以下のように書けることがわかりました。
 
\begin{align*}
\frac{u_{i}^{n+1}-u_{i}^{n}}{\Delta t}=\nu\frac{u_{i+1}^{n}-2u_{i}^{n}+u_{i-1}^{n}}{\Delta x^2}\tag{7}
\end{align*}
 
今、ほしい情報は「ある時刻置ける\(u(x)\)」です。
 
未来を知りたいのですから、今のステップ数を\(n\)とすると、ちょっと先の未来は\(n+1\)ステップです。
\begin{align*}
u_{i}^{n+1}=u_{i}^{n}+\frac{\nu\Delta t}{\Delta x^2}(u_{i+1}^{n}-2u_{i}^{n}+u_{i-1}^{n})\tag{8}
\end{align*}
 
 
では、さっそくPythonを使った拡散方程式の数値計算をと・・・・行きたいところですが、前回の記事で確認したような数値的な安定性条件を見ておきます
 

フォン・ノイマンのフーリエ解析法

 
今回使う方法は、「フォン・ノイマンのフーリエ解析法」と呼ばれるものです。
名前だけで圧倒されそうな解法ですが、ポイントは以下となります。

【ポイント】

  1. 関数を波に分解→フーリエ変換(or フーリエ級数展開)
  2. 時間が経つと波の振幅が大きくなる(大きくなると不安定)

 
フーリエ変換について学びたい方は☟こちらからどうぞ!

「フーリエ変換」の記事はこちら
 
では、本題に行きます(^^)/
 

時間が経つと波の振幅が大きくなる(大きくなると不安定)

 
(1)式を離散化して、時々刻々と変動する数値解を安定的求めるには、計算を進めるにあたって数値的な誤差が増大しないことが条件になります。
移流方程式には厳密解が存在するので、厳密解をとして\(\tilde{u}^{n}_{i}\)、数値解を\({u}^{n}_{i}\)とします。
そうすると、厳密解と数値解には誤差があるので誤差を\(\tilde{\epsilon}^{n}_{i}\)と置いておきしょう。

\begin{align*}
u^{n}_{i}=\tilde{u}^{n}_{i}+\tilde{\epsilon}^{n}_{i}
\end{align*}

これを(8)式に代入すると、誤差に関する方程式を得ることができます。
\begin{align*}
\tilde{\epsilon}^{n+1}_{i}=\tilde{\epsilon}^{n}_{i}+\frac{\nu\Delta t}{\Delta x^2}(\tilde{\epsilon}^{n}_{i+1}-2\tilde{\epsilon}^{n}_{i}+\tilde{\epsilon}^{n}_{i-1})\tag{8-1}
\end{align*}

誤差\(\tilde{\epsilon}^{n}_{i}\)、時刻\(n\)ステップと格子点\(i\)で以下のようにフーリエ級数展開を行います。
\begin{align*}
\tilde{\epsilon}^{n}_{i}=\sum A^{n}_{k}e^{Ik(i\Delta x)}\tag{9}
\end{align*}

これは、波数\(k\)の波\(e^{Ik(i\Delta x)}\)が\(A^{n}_{k}\)で重み付けされた形になっています。
※虚数は、格子点の\(i\)と紛らわしいので「虚数\(I\)」と表記しました。
 
これを(8-1)式に代入します。
(8-1)式は線形な方程式なので(9)式の和のほとんどは書かなくてもよくなります。
\begin{align*}
A^{n+1}_{k}=A^{n}_{k}+\alpha\bigg(A^{n}_{k}e^{Ik(i\Delta x)}-2A^{n}_{k}+A^{n}_{k}e^{-Ik(i\Delta x)}\bigg)\tag{10}
\end{align*}

※\(\alpha=\frac{\nu\Delta t}{\Delta x^2}\)と置いています。
これを今、以下のようにまとめます。
\begin{align*}
g& = \frac{A^{n+1}_{k}}{A^{n}_{k}}\\
& =1+\alpha\big(e^{Ik(i\Delta x)}-2+e^{-Ik(i\Delta x)}\big)\\
& =1-2\alpha\big(1-\cos\phi\big)
\end{align*}

※\(\phi=Ik\Delta x\)と置きました。
ちょっとわかりにくいので絵でイメージすると以下のようになります。
 
\(g = \frac{A^{n+1}_{k}}{A^{n}_{k}}\)が1より大きいという事は、ある波数\(k\)の重みがどんどん増大してやがて数値的には無視できないくらいの誤差になるということを意味しています。
つまり、「n+1ステップでの波数\(k\)の重み」と「nステップでの波数\(k\)の重み」\(g = \frac{A^{n+1}_{k}}{A^{n}_{k}}\)を考えたときに、
\begin{align*}
\left\{\begin{matrix}
|g|\leq 1\,\,\,\,(安定)\\
|g|> 1\,\,\,\,(不安定)
\end{matrix}\right.
\end{align*}

という判別ができます。
どんな\(\phi\)に対しても\(|g|\leq 1\)であれば安定解を得るということですね。
つまり・・・
\begin{align*}
-1\leq \alpha \leq 1
\end{align*}

すなわち、
\begin{align*}
-1\leq 1-2\alpha\big(1-\cos\phi\big)\leq 1\tag{11}
\end{align*}

となります。
今の場合、\(g\)はどんな\(\phi\)に対しても1より小さいので問題になってくるのは、\(-1\leq 1-2\alpha\bigg(1-\cos\phi\bigg)\)です。
一番厳しい条件は\(\cos\phi=-1\)なので、どんな\(\phi\)に対しても(11)が成立するためには、
\begin{align*}
\alpha\leq \frac{1}{2}\tag{12}
\end{align*}

となります。
ここで、数値計算においては\(\Delta x\)を決めた後に時間幅\(\Delta t\)をいくつにしようか考えることが多いので、\(\Delta t\)についての条件式に書き換えておきます。
\begin{align*}
\Delta t\leq \frac{\nu\Delta t}{2\Delta x^2}\tag{13}
\end{align*}

このように、空間幅\(\Delta x\)に対して時間幅\(\Delta t\)は(13)式を満たしていなければ安定ではない解となってしまうということを意味しています。
では、いよいよPythonを使って「1次元の拡散方程式」を解いてみましょう。
 
(13)の条件を「満たすもの」と「満たさないもの」に分けててシミュレーションをしてみます。

カマキリ

ちなみに(13)の条件って名前があるんですかね?(/・ω・)/

「\(\alpha=\)0.1、0.5、1.0」と変えてシミュレーション

 

では、Google Colaboratoryを使ってコードを書いて、解の振る舞いをアニメーションで見て理解を深めていきたいと思います。

(13)式より、\(\alpha\)は0.5より小さくないと安定解を得ることができないとわかったので、色々値を変えて検証していきます。

【ポイント】
以下の3パターンで解析をします。

\(\alpha\)解の振る舞い
0.1
0.5
1.0
 
 
基本となるPythonコードは以下となります。
前回の記事と変わった点は、

  • 偏微分方程式を解く部分が「u[i] = un[i] + nu * dt / dx**2 * (un[i+1] – 2 * un[i] + un[i-1])」となっています。
  • un[i+1]があるため空間領域を「for i in range(1, nx – 1):」としています。
カマキリ

では、検証スタート!

\(\alpha =0.1\)

 

この条件は、安定的に解ける条件のはずですがどうなるでしょうか。
☟動画を見てください。

安定していますね。

\(\alpha =0.5\)

 

この条件は、ギリギリ安定的に解ける条件のはずですがどうなるでしょうか。
☟動画を見てください。

じゃっかんギザギザしていますね。

\(\alpha =1.0\)

 

この条件は、不安定な解の振る舞いをする条件のはずですがどうなるでしょうか。
☟動画を見てください。

見てられないくらいめちゃくちゃな解の振る舞いをしてしまいました!!

結果まとめ

 

\(\alpha\)解の振る舞い
0.5安定
拡散する
1.0ちょっと不安定
じゃっかんガタガタ
1.5不安定
見てられない

\begin{align*}
\left\{\begin{matrix}
|g|\leq 1\,\,\,\,(安定)\\
|g|> 1\,\,\,\,(不安定)
\end{matrix}\right.
\end{align*}

の条件と一致している数値計算の結果ですね。

まとめ

 
本記事は、
  • 一次元拡散方程式の数値計算
  • 拡散方程式の安定性
についてまとめました。
拡散方程式の空間微分を中心差分で離散化した際は以下の条件を満たすように時間幅を設定しないと不安定になるということを確認しました。
 
\begin{align*}
\Delta t\leq \frac{\nu\Delta t}{2\Delta x^2}\tag{13}
\end{align*}
 
このように流体の数値計算においては、条件設定次第では解が不安定になったりする可能性があるためとても注意を払わないといけないということがわかりますね。
 
 
 

Pythonの完全初心者は書籍で学ぶとよい

 
Python自体が不安だという方は、こちらの書籍から勉強しても良いかも知れません。
☟こちらの本がPython初心者が挫折することなく勉強できる本です。
(本記事のようのPython使用環境と異なりますが、とてもわかりやすいので全く問題ありません)

確かな力が身につくPython「超」入門 (確かな力が身につく「超」入門シリーズ)

確かな力が身につくPython「超」入門 (確かな力が身につく「超」入門シリーズ)

鎌田 正浩
2,728円(10/31 20:12時点)
発売日: 2016/03/16
Amazonの情報を掲載しています
僕自身もこの書籍ではじめは勉強しました。
小難しい余計なことが書いていないし、年末の一週間くらいで読むことができたので、初心者の方にはとてもお勧めです。
 

流体の勉強をしたい方

 
流体の数値計算以前に、流体力学を勉強しないと!!って思っている方は以下の参考書からはじめてみるのが良いと思います。

流体力学 (JSMEテキストシリーズ)

流体力学 (JSMEテキストシリーズ)

日本機械学会
2,075円(10/31 11:49時点)
発売日: 2005/04/01
Amazonの情報を掲載しています

工学よりに書かれていて実用的な内容も触れながら、流体の理論的な内容も易しく書かれていますので大学1年生で偏微分を学んだ方にとってはちょうどよい参考書です。

数値計算を勉強したい方

 
数値計算は微分を離散化して代数的に解くため、どうしても誤差が発生します。
それをいかに誤差が少なく、安定的に解くかがカギになっており、数値的解法は多く提案されています。
例えば以下の参考書のような、数値計算に特化した参考書を一冊手元に置きながら勉強を進めることで「流体の数値計算」の理解が深まると思います。
 
数値計算 (理工系の基礎数学 8)

数値計算 (理工系の基礎数学 8)

高橋 大輔
3,190円(11/01 08:23時点)
発売日: 1996/02/16
Amazonの情報を掲載しています
 
 
 
 

【プロフィール】

カマキリ
(^^)

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

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

 

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

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

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