Python

【第6回Python流体の数値計算】1次元移流方程式のクーラン条件を確認!数値的安定性(差分法)



こんにちは(@t_kun_kamakiri)。

今回は、Pythonを使って1次元の移流方程式」のクーラン条件のお話をします。

クーラン条件は数値的安定性とも関係があるので、必ず学んでおいた方が良いですね。
本記事では、わざと数値的に不安定な状態での数値計算もやりますので、解の振る舞いがどのようになるのかを確認しながら見ていただければと思います。

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

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

  • Pythonを使い始めたけどどう使うかわからない
  • 流体の数値計算をはじめて勉強する人

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

前回の記事はこちら
 
 
前回は「1次元の移流方程式」のアニメーション作成をしましたが、今回はそのアニメーションを使って数値的安定性について考えていきたいと思います。
 
 
今回の内容はこちら
移流方程式における安定性条件(クーラン条件)を考える。

 
では、Google Colaboratoryを使ってコードを書きながら理解を深めていきたいと思います。
 
スポンサーリンク

前回作成した全体のコード:1次元移流方程式

 

前回は、1次元移流方程式を後退差分で解くプログラムをPythonで書きました。

解いている方程式は、

1次元移流方程式の離散化(後退差分)

\begin{align*}
u^{n+1}_{i}=u^{n}_{i}-\frac{c\Delta t}{\Delta x}\big(u^{n}_{i}-u^{n}_{i-1}\big)\tag{1}
\end{align*}

です。

数値的安定性

 

(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}\tag{2}
\end{align*}

今考えるべきポイントは、
いかに数値誤差\(\tilde{\epsilon}^{n}_{i}を小さくすることができるか\)
です。

(2)式を(1)式に代入します。

\begin{align*}
\tilde{u}^{n+1}_{i}+\tilde{\epsilon}^{n+1}_{i}&=\tilde{u}^{n}_{i}+\tilde{\epsilon}^{n}_{i}\\
&-\frac{c\Delta t}{\Delta x}\big(\tilde{u}^{n}_{i}-\tilde{u}^{n}_{i-1}\big)\\
&-\frac{c\Delta t}{\Delta x}\big(\tilde{\epsilon}^{n}_{i}-\tilde{\epsilon}^{n}_{i-1}\big)\tag{3}
\end{align*}

じゃっかんごちゃごちゃしてしまいましたが、\(\tilde{u}^{n}_{i}\)は厳密解なのですからそもそも、

\begin{align*}
\tilde{u}^{n+1}_{i}=\tilde{u}^{n}_{i}-\frac{c\Delta t}{\Delta x}\big(\tilde{u}^{n}_{i}-\tilde{u}^{n}_{i-1}\big)
\end{align*}

を満たしています。
移流方程式が線形な方程式であるため、線形解である\(u^{n}_{i}=\tilde{u}^{n}_{i}+\tilde{\epsilon}^{n}_{i}\)のうちの右辺のどちらも移流方程式を満たすという事です。

今、考えるべきは以下の式です。

\begin{align*}
\tilde{\epsilon}^{n+1}_{i}=\tilde{\epsilon}^{n}_{i}-\frac{c\Delta t}{\Delta x}\big(\tilde{\epsilon}^{n}_{i}-\tilde{\epsilon}^{n}_{i-1}\big)\tag{4}
\end{align*}

(4)式は厳密解というよりは、そもそも誤差なので\(\tilde{\epsilon}^{n}_{i}=0\)の予定ですが、(4)式がある以上はこの式に従ってもらいながら発散されるととても困ることになります・・・・

 

この誤差\(\tilde{\epsilon}^{n}_{i}\)が時間とともに増大するかどうかを考えて、増大するような数値解法だと解が不安定になるというわけです。

カマキリ

どうやって増大するかどうかを調べるのか?

例えば、以下のようにまとめると数値的安定性条件が出てきそうなのがわかります。

\begin{align*}
\tilde{\epsilon}^{n+1}_{i}=\bigg(1-\frac{c\Delta t}{\Delta x}\bigg)\tilde{\epsilon}^{n}_{i}+\tilde{\epsilon}^{n}_{i-1}\tag{4-1}
\end{align*}

この式を見たときに、

\begin{align*}
1-\frac{c\Delta t}{\Delta x}\geq 0
\end{align*}

じゃないと、解が安定しないのではないかというのが予想できます。
なぜなら、(4-1)式の左辺は「n+1ステップ」の「格子点i」での誤差に対して、右辺は「nステップ」の「格子点iとi-1」の誤差です。
あるステップのときのある場所ではマイナスの値をとったりプラスの値を取ったりして解が不連続になりやすいと予想できます。
そうすると、不連続点の位置での微分は値が発散するということになるので(今は離散化しているのでせいぜい大きな値をとるくらい)、解が不安定になりやすい・・・・とちょっぴり予想できます。

カマキリ

しかし、予想だけではいけないので別の方法で安定条件を導きます!

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

 

今回使う方法は、「フォン・ノイマンのフーリエ解析法」と呼ばれるものです。

名前だけで圧倒されそうな解法ですが、ポイントは以下となります。

【ポイント】

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

 

フーリエ変換の簡単なお話をしておきましょう(^^)/

関数を波に分解→フーリエ変換(or フーリエ級数展開)

 

フーリエ変換は下記のようにある関数\(f(x)\)を、\(\sin\)波や\(\cos\)波や正弦波\(e^{ikx}\)などの様々な異なる波数の波に分解して波の空間(波数空間)に変換する方法のことです。

どんな関数も異なる波数の\(\sin\)波の和で書いてやるのです。↓こんな感じです。
※以下の場合は周期関数\(f(x)\)に関して、例を挙げています。

こうやって考えると、「おや?\(\sin2\pi x\)(周期1)の波が1個か・・・・\(\sin10\pi x\)(周期1/10)の波が10個か・・・」と、波のくせに数えれるという恩恵を受けることができます。

  • 「1個、10個、3個」とか書いたものがフーリエ変換(or フーリエ級数展開)した後の関数(波数空間での関数)のこと
  • 「1個、10個、3個」は各波数の波の振幅

と考えることができます。

    結局複雑な波形でも、ある波数を持つ波の重み付けによる重ね合わせで表現できるのです。

つまり、数値的に不安定な解とはこういうことです。

「1個、10個、3個」などと書いた各波数の波の振幅時間が経つにつれて、どんどん大きくなってしまう場合が不安定な数値解である。

誤差\(\tilde{\epsilon}^{n}_{i}\)をフーリエ級数展開して、ある波数の波が増大する可能性があるのかどうかを調べることで、解の安定性条件が導くことができます。

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

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

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

 

誤差\(\tilde{\epsilon}^{n}_{i}\)、時刻\(n\)ステップと格子点\(i\)で以下のようにフーリエ級数展開を行います。

\begin{align*}
\tilde{\epsilon}^{n}_{i}=\sum A^{n}_{k}e^{Ik(i\Delta x)}\tag{5}
\end{align*}

これは、波数\(k\)の波\(e^{Ik(i\Delta x)}\)が\(A^{n}_{k}\)で重み付けされた形になっています。
※虚数は、格子点の\(i\)と紛らわしいので「虚数\(I\)」と表記しました。

これを(4)式に代入します。
(4)式は線形な方程式なので(5)式の和のほとんどは書かなくてもよくなります。

\begin{align*}
A^{n+1}_{k}=A^{n}_{k}-\alpha\bigg(A^{n}_{k}-A^{n}_{k}e^{-Ik(i\Delta x)}\bigg)\tag{6}
\end{align*}

※\(\alpha=\frac{c\Delta t}{\Delta x}\)と置いています。

これを今、以下のようにまとめます。

\begin{align*}
g& = \frac{A^{n+1}_{k}}{A^{n}_{k}}\\
& =1-\alpha\bigg(1-e^{-Ik(i\Delta x)}\bigg)\\
& = 1-2\alpha\sin^2\frac{\phi}{2}-I\alpha\sin\phi\\
& = x+yI
\end{align*}

※実部と虚部に分けました。
※\(\phi=Ik\Delta x\)と置きました。

だから、

\begin{align*}
|g|=\bigg(4\sin^2\frac{\phi}{2}\big(\alpha-\frac{1}{2}\big)-\sin^2\frac{\phi}{2}+1\bigg)\tag{7}
\end{align*}

となります。
ちょっとわかりにくいので絵でイメージすると以下のようになります。

 

\(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*}

という判別ができます。

じゃー、変数\(\alpha\)として\(g\)のグラフを描いて\(g\leq 1\)となる条件を導けばよいという事になります。2次関数のグラフを描いて安定領域を探せば、

\begin{align*}
0\leq \alpha \leq 1
\end{align*}

すなわち、

\begin{align*}
0\leq \frac{c\Delta t}{\Delta x}\leq 1\tag{8}
\end{align*}

となります。

(8)式を満たせば、\(|g|\leq 1\)を満たしますね。

先に予想していたような結果を得ることができました。

ちなみに(8)式をクーラン条件といいますので覚えておきましょう(^^)/

\(\alpha=\frac{c\Delta t}{\Delta x}\)をクーラン数といいます。

では、クーラン数の値を変えてシミュレーションをしてみましょう。

クーラン数を「\(\alpha=\)0.5、1、1.5」と変えてシミュレーション

 

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

\(\alpha\)解の振る舞い
0.5
1.0
1.5

 

カマキリ

さて、どうなるでしょう?

Pythonコードを書きます。

ここから\(\alpha\)の値を変えて、アニメーションを見てみます。

\(\alpha =0.5\)

 

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

解が安定はしていますが、「+x方向に進みながら波形の形が変わっている」という結果でした。

 

\(\alpha =1.0\)

 

続いて・・・

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

おや(/・ω・)/

安定的に解けてるし、「+x方向に形を変えずに進んでいる・・・」!

\(\alpha =1.5\)

 

これはどうでしょうか?

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

 

動画を見ると、みごとにめちゃくちゃな解の振る舞いをしてしまいました!!

結果まとめ

 

\(\alpha\)解の振る舞い
0.5安定
形が変わって進む
1.0安定
形が変わらず
1.5不安定
見てられない

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

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

まとめ

 

本記事の内容は「\(\alpha>0\)かつ後退差分をした際」の安定条件が、

\begin{align*}
0\leq \frac{c\Delta t}{\Delta x}\leq 1\tag{8}
\end{align*}

となることを示しました。

(8)式をクーラン条件といいます。

ここで注意しておきたいのは、本記事の内容は「\(\alpha>0\)かつ後退差分をした際」の安定条件についてですが、

\(\alpha<0\)である場合に関しては、後退差分ではいかなる条件も不安定な解の振る舞いをします。

もし、\(\alpha<0\)の条件で数値計算を行うのであれば前進差分で解く方が安定性が良いです。
それは以下の記事でもまとめていますので、良ければ参考にしてください。
コード自体はFortranを使用していますが、コードは無視して動画を見るだけでも勉強になると思います。

こちらの記事もどうぞ
 
上記の記事で以下の結果を得ています。
 
 
 

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ではプログラミング言語の基本的な内容をまとめています。
関連記事もどうぞ