高校物理、大学受験、大学物理、プログラミングの記事を書きます。

【第15回Python流体の数値計算】2次元ポアソン方程式をPythonで実装する。

2020/06/26
 




この記事を書いている人 - WRITER -

こんにちは(@t_kun_kamakiri)。

前回の記事では、「2次元のラプラス方程式」をGoogle Colab上でPythonで実装しました。

前回の記事はこちら

 

今回は、以下のような2次元のポアソン方程式をPythonで実装します。

なぜ、ポアソン方程式の解法を知る必要があるのかというと、ナビエストークス方程式を解く際に必要になるからです。(前回の記事で詳しく書きました。)
 
 
 

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

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

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

2次元のポアソン方程式の数値計算をPythonで実装

 

カマキリ
今日もゆる~く学びます(^^)/
スポンサーリンク

ポアソン方程式

 

まずは、

  • ポアソン方程式がどういう式なのか
  • ナビエストークス方程式のを解く際のポアソン方程式の使われ方

を復習しておきましょう。

流体現象を扱う際に「ナビエストークス方程式」と「質量保存則」を連立して解く必要があります。

ナビエストークス方程式

\begin{align*}
\frac{\partial u}{\partial t}+u\frac{\partial u}{\partial x}+v\frac{\partial u}{\partial y} = -\frac{1}{\rho}\frac{\partial p}{\partial x}+\nu \left(\frac{\partial^2 u}{\partial x^2}+\frac{\partial^2 u}{\partial y^2} \right)\tag{1}
\end{align*}
 
\begin{align*}
\frac{\partial v}{\partial t}+u\frac{\partial v}{\partial x}+v\frac{\partial v}{\partial y} = -\frac{1}{\rho}\frac{\partial p}{\partial y}+\nu\left(\frac{\partial^2 v}{\partial x^2}+\frac{\partial^2 v}{\partial y^2}\right) \tag{2}
\end{align*}

 
質量保存則(非圧縮条件\(\nabla\cdot\boldsymbol{v}=0\)から導出)
\begin{align*}
\frac{\partial^2 p}{\partial x^2}+\frac{\partial^2 p}{\partial y^2} = b(x,y)\tag{3.1}
\end{align*}

 
\begin{align*}
b(x,y)=-\rho\left(\frac{\partial u}{\partial x}\frac{\partial u}{\partial x}+2\frac{\partial u}{\partial y}\frac{\partial v}{\partial x}+\frac{\partial v}{\partial y}\frac{\partial v}{\partial y} \right)\tag{3.2}
\end{align*}
(3.1)式の形がポアソン方程式
 
このように、ナビエストークス方程式とポアソン方程式を連立させて解く必要があります。
 
右辺が0である場合のラプラス方程式、
\begin{align*}
\frac{\partial ^2 p}{\partial x^2} + \frac{\partial ^2 p}{\partial y^2} = 0
\end{align*}
は、前回の記事で扱いましたが、今回は以下のポアソン方程式を扱います。
 
ポアソン方程式

\begin{align*}
\frac{\partial ^2 p}{\partial x^2} + \frac{\partial ^2 p}{\partial y^2} = b(x,y)
\end{align*}

実際の流体の方程式としては、一般的には右辺が0ではないです。

なので、今回の趣旨はポアソン方程式をPythonを使って実装しようということになります。
 
 

 
 
ポアソン方程式として有名なのは、以下の式でしょうかね。
 
\begin{align*}
\frac{\partial ^2 \phi}{\partial x^2} + \frac{\partial ^2 \phi}{\partial y^2} = -\frac{\rho(\boldsymbol{r})}{\epsilon}
\end{align*}
 

これは、電磁気で習う点電荷まわりのスカラー場を求める偏微分方程式と全く同じですね。
圧力\(p\)を電場\(E=-\nabla \phi\)のスカラー場\(\phi\)と置き換えると良いでしょう。
 
 
 

 
 

ポアソン方程式を離散化

 

ポアソン方程式は偏微分の方程式なので、離散化して代数的に取り扱えるようにしないと数値計算ができません。

離散化の方法としては、空間に対しては2次精度中心差分にして\(i,j\)での\(p_{i,j}\)として以下のようにして求めます。
 
\begin{align*}
\frac{p_{i+1,j}^{n}-2p_{i,j}^{n}+p_{i-1,j}^{n}}{\Delta x^2}+\frac{p_{i,j+1}^{n}-2 p_{i,j}^{n}+p_{i,j-1}^{n}}{\Delta y^2}=b_{i,j}^{n}\tag{4}
\end{align*}
 
よって、
 
\begin{align*}
p_{i,j}^{n}=\frac{(p_{i+1,j}^{n}+p_{i-1,j}^{n})\Delta y^2+(p_{i,j+1}^{n}+p_{i,j-1}^{n})\Delta x^2-b_{i,j}^{n}\Delta x^2\Delta y^2}{2(\Delta x^2+\Delta y^2)}\tag{5}
\end{align*}
 
特に、\(\Delta x= \Delta y\)の場合は、
 
\begin{align*}
p_{i,j}^{n}=\frac{p_{i+1,j}^{n}+p_{i-1,j}^{n}+p_{i,j+1}^{n}+p_{i,j-1}^{n}-b_{i,j}^{n}\Delta x^2}{4}\tag{6}
\end{align*}
 
例えば(1,1)の格子点上の圧力は上記のように以下のようになります。

 
青文字は境界条件で値が与えられているので、既知の値です。
これを見ればわかるように、(i,j)の格子点での圧力はその周りの圧力の平均値として計算しています。

境界条件を設定

 

ポアソン方程式は、時間発展の方程式ではなく、境界条件によって空間分布の解が与えられる方程式なので、境界条件を設定する必要があります。

今回与える境界条件は以下のようにします。

  • \(x=0,  2\):\(p=0\)
  • \(y=0,  1\):\(p=0\)
次に\(b(\boldsymbol{r})\)の条件を設定しましょう。
 
  • \(i=\frac{1}{4}nx, j=\frac{1}{4}ny\):\(b_{i,j}=100\)
  • \(i=\frac{3}{4}nx, j=\frac{3}{4}ny\):\(b_{i,j}=-100\)

ちょっとわかりにくいかもしれないので、絵にしてみました。

 

では、以下で実際にポアソン方程式をGoogle Colaboratoryを使ってPythonのコードを書きながら理解を深めていきたいと思います。
 

ポアソン方程式をPythonで実装

 
ひとつひとつ解説をしながらコードを書いていきます。
※ポアソン方程式は、前回の☟こちらの記事で書いているので是非参考にしてみてください。
 
前回の記事はこちら
 
 
 
【ポイント】

  1. 必要なライブラリをインポートする。

 
コードを以下のように書きます。

NumpyとMatplotlibをインポートします。
記述ミスがなければエラーなく進みます。

 
では、次にポアソン方程式における「初期の状態」と「境界条件」を課します。
 
【ポイント】

  1. 初期状態を設定
  2. 境界条件を設定

 
コードを以下のように書きます。

これで、初期状態の圧力のプロファイルができています。
初期状態と言っても、境界条件を課しただけですので初期の状態には何の意味もありません。

あくまで境界条件をもとにラプラス方程式を解いて得られた圧力分布に意味があるので、初期状態はその最終状態を得るための種を与えているにすぎません

 

では、いよいよ上記で設定した境界条件をもとにラプラス方程式を解いてみましょう!
 
 
【ポイント】

\begin{align*}
p_{i,j}^{n}=\frac{p_{i+1,j}^{n}+p_{i-1,j}^{n}+p_{i,j+1}^{n}+p_{i,j-1}^{n}-b_{i,j}^{n}\Delta x^2}{4}\tag{6}
\end{align*}
 
 

コードを以下のように書きます。

これで最終状態の圧力の分布が計算されました。

前回のラプラス方程式を解いた際には、全ての格子点上の値の合計と前のステップでの値の合計とを比較して、「1e-4(0.01%)」以下の誤差になれば収束したと判定していましたが、今回は収束条件を反復回数(100回)にしています。

※反復計算には「ヤコビ法」と呼ばれる方法で計算しています。

では、次に収束した圧力分布を可視化してみてみましょう。

【ポイント】

Matplotlibを使って結果を可視化する。
 
 

コードを以下のように書きます。

これで、可視化するための関数が完成したので、以下のように関数を呼び出して圧力分布を可視化してみましょう。

【結果】

このようになればOKです(^^)/

コンター図も作成してみる

 

Matplotlibの公式ドキュメントは非常にわかりやすい(たぶんPythonがわかりやすいから)ので、コンター図もサクッとさくせいできます。

【結果】

このように圧力の分布もコンターで可視化できます。

まとめ

 

今回は、2次元のポアソン方程式をPythonで実装してみました。

本記事のシリーズを読むことでPythonを使ったナビエストークス方程式の実装まではいけると思いますので、是非最後までお付き合いくださいm(__)m

次回は、いよいよナビエストークス方程式をPythonで実装してみたいと思います(^^)/

次回の記事はこちら

おすすめの書籍紹介

 

ネット上には無料で多くのことが学べるのですが、ネット情報だけだと情報が偏ったりします。

やっぱり書籍を手に取って体系だって学んだ方が良い場合も多いです。

初心者でもので体系だって学ぶことができる、且つあまり難しすぎない書籍をここで紹介しておきます。

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

 

Python自体が不安だという方は、こちらの書籍から勉強しても良いかも知れません。

☟こちらの本がPython初心者が挫折することなく勉強できる本です。
(本記事のようのPython使用環境と異なりますが、とてもわかりやすいので全く問題ありません)

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

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

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

流体の勉強をしたい方

 

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

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

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

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

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

数値流体の書籍

 

本記事で紹介している数値計算手法は、多くの数値解法のごく一部です。

特に、非線形方程式で記述できる流体現象のようなものは注意深く数値解法を選択しなくてはいけません。

以下の書籍は、値段が高いですが有限体積法をメインに数値解法がとても詳しく解説されています。

数値流体力学 第2版

数値流体力学 第2版

H.K.Versteeg, W.Malalasekera
10,450円(07/10 20:14時点)
発売日: 2011/05/31
Amazonの情報を掲載しています

 

数値計算の書籍

 

今回紹介したポアソン方程式は「ヤコビ法」と呼ばれる反復計算を行っていますが、その他の有名な数値的解法のをいくつか挙げておきます。

  • ガウス消去法(直接法)
  • ガウス・ジョルダン法(直接法)
  • ヤコビ法(反復法)
  • ガウス・ザイデル法(反復法)
  • SOR法(反復法)

などがあります。

詳しく勉強したい方は、以下の参考書がとてもわかりやすいです。

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

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

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

数値計算の闇について詳しく知りたい方は、☟こちらの参考書も合わせてお読みください。

数値計算の常識

数値計算の常識

伊理 正夫, 藤野 和建
8,250円(07/11 05:46時点)
発売日: 1985/06/03
Amazonの情報を掲載しています

 

この記事を書いている人 - WRITER -

Comment

  1. 牧田 より:

    (5)式から(6)式への式変形のところで質問があります。
    Δx=Δyとした場合、(6)式のb_(i,j)^nの部分は、b_(i,j)^nではなく、(b_(i,j)^n)(Δy)^2 or (b_(i,j)^n)(Δx)^2ではないのですか?

    私の計算ミスでしたらごめんなさい。

    よろしくお願いします。

    • korokoro より:

      ご指摘ありがとうございます。
      おっしゃるとおりです。失礼致しました。
      修正したいと思います。
      数値計算の方は、修正の必要は無く正しいです。

      ありがとうございました。

コメントを残す

Copyright© 宇宙に入ったカマキリ , 2020 All Rights Reserved.