Python

【第14回Python流体の数値計算】2次元ラプラス方程式をPythonで実装する。

こんにちは(@t_kun_kamakiri)。

前回の記事では、「2次元のバーガース方程式」をGoogle Colabでアニメーション作成を行いました。

前回の記事はこちら

 

今回は、以下のような2次元のラプラス方程式をPythonで実装します。

なぜ、ラプラス方程式の解法を知る必要があるのかというと、ナビエストークス方程式を解く際に必要になるからです。(正確には、必要なのはポアソン方程式です。)
本記事では、そちらも詳しく解説します。
 
 

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

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

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

  • 2次元のラプラス方程式の数値的解法をなぜ知る必要があるのか
  • 2次元のラプラス方程式の数値計算をPythonで実装

 

カマキリ

今日もゆる~く学びます(^^)/

スポンサーリンク

2次元のラプラス方程式

 

今回は以下のラプラス方程式ついてPythonを用いて実装して数値計算を行いたいと思います。
 
\begin{align*}
\frac{\partial ^2 p}{\partial x^2} + \frac{\partial ^2 p}{\partial y^2} = 0
\tag{1}
\end{align*}
 
いきなりラプラス方程式の話をされても、なぜここでラプラス方程式の話をするのかと、疑問に思う人もいることでしょう!
 
 

  • ラプラス方程式とポアソン方程式の違い
  • ナビエストークス方程式の数値計算にポアソン方程式を解く場面がある

 

ラプラス方程式とポアソン方程式の違い

 

ラプラス方程式の形をしていても、一般的には右辺が0ではない場合、

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

ポアソン方程式と呼ばれています。

ポアソン方程式として最も有名な式としては、「電荷の分布\(\rho(\boldsymbol{r})\)があります。
周りの静電場\(\boldsymbol{E}\)」マクスウェル方程式\(\nabla\cdot\boldsymbol{E}=\frac{\rho(\boldsymbol{r})}{\epsilon}\)を満たし、スカラー場\(\phi\)としたとき、スカラー場の勾配が電場\(\boldsymbol{E}=-\nabla\phi\)なので、
\begin{align*}
\nabla^2 \phi =-\frac{\rho(\boldsymbol{r})}{\epsilon}
\end{align*}
 

電荷が0の場合は、ラプラス方程式になるというわけですね。

明らかにポアソン方程式の難しそうなので、まずは理解のために簡単そうなラプラス方程式について考えてみましょう。
というのが本記事の狙いでもあります。
 
 

ナビエストークス方程式の数値計算にポアソン方程式を解く場面がある

 
ところで、なぜラプラス方程式をここで持ち出したかというと、

ナビエストークス方程式の数値計算にポアソン方程式を解く場面があるからです。

非圧縮性流体として現象を取り扱う場合には、基本的には以下の2つの式を連立して解きます。
 
 
ナビエストークス方程式
 
\begin{align*}
\frac{\partial \boldsymbol{v}}{\partial t}+(\boldsymbol{v}\cdot\nabla)\boldsymbol{v}=-\frac{1}{\rho}\nabla p + \nu \nabla^2\boldsymbol{v}\tag{3}
\end{align*}

※ナビエストークス方程式は成分の数の式があります。

連続の式(非圧縮条件)
\begin{align*}
\nabla\cdot\boldsymbol{v}=0\tag{4}
\end{align*}
 
実際の流れは、すべての現象が(3)(4)式というわけではないで。
 
「ニュートン流体」「非圧縮の流れ(近似)」としているとても限定的な流れであることを理解しておきましょう。
もっと詳しく勉強したい方は☟こちらの記事を参考にしてください。
 
 
今回は、2次元で計算を行うため少々くどいですが、成分ごとの式を書き下すことにします。
 
 
【2次元での計算】
 
ナビエストークス方程式
 
\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{5}
\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{6}
\end{align*}
 
連続の式(非圧縮条件)
\begin{align*}
\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}=0\tag{7}
\end{align*}
 
  • ナビエストークスからは2本の式
  • 連続の式から1つの式
が出てくるので、合計で3本の式があります。

未知数は、\(\boldsymbol{v}\)の2つと圧力\(p\)の合計3つです。
なので、原理的に未知数である物理量は全て求まることになります。
それを各時刻で求めていけば良いという事になります。

しかし、言葉で簡単に言っても、特殊な境界条件や近似を用いないとナビエストークス方程式は解析的に取り扱うことができないのです。
つまり、ナビエストークス方程式には一般解が見つかっていないので、実際には数値計算を使って偏微分方程式を取り扱うしかないのが現状です。
 
しかし、連立方程式を解けばいいと言われてもちょっと困った・・・
なぜなら、(5)(6)式で求めた流速は(7)式を満たす必要があります。
 
そのように流速を選びながら圧力を同時に決めなければなりません。
考えただけでもちょっと工夫が必要ですよね。。
 
 
そこで、以下のようにして考えを改めます。
 
 
圧力に関しては非圧縮条件を満たすようにポアソン方程式を解きつつ、求めた圧力に関してナビエストークス方程式を解き流速を求めるという手法をとります。
 
 
言葉で説明してもよくわからないかもしれないので、実際にやります(^^)/
 
(5)式を\(y\)で偏微分
 
\begin{align*}
&\frac{\partial }{\partial x}\frac{\partial u}{\partial t}+\frac{\partial u}{\partial x}\frac{\partial u}{\partial x}+u\frac{\partial^2 u}{\partial x^2}+ \frac{\partial v}{\partial x}\frac{\partial u}{\partial y}+v\frac{\partial^2 u}{\partial x\partial y} \\
&= -\frac{1}{\rho}\frac{\partial^2 p}{\partial x^2}+\nu \frac{\partial }{\partial x}\left(\frac{\partial^2 u}{\partial x^2}+\frac{\partial^2 u}{\partial y^2} \right)\tag{8}
\end{align*}
 
(6)式を\(x\)で偏微分
 
\begin{align*}
&\frac{\partial }{\partial y}\frac{\partial v}{\partial t}+\frac{\partial u}{\partial y}\frac{\partial v}{\partial x}+u\frac{\partial^2 v}{\partial y\partial x}+ \frac{\partial v}{\partial y}\frac{\partial v}{\partial y}+v\frac{\partial^2 v}{\partial x\partial y^2} \\
&=-\frac{1}{\rho}\frac{\partial^2 p}{\partial y^2}+\nu \frac{\partial }{\partial y}\left(\frac{\partial^2 u}{\partial x^2}+\frac{\partial^2 u}{\partial y^2} \right)\tag{9}
\end{align*}
 
このようにして辺々足して、非圧縮条件を使えば、
 
 
\begin{align*}
\frac{\partial^2 p}{\partial x^2}+\frac{\partial^2 p}{\partial y^2} = -\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{10}
\end{align*}
 
計算はだるいですが、一度やってみてください。
 
つまり、(5)(6)(7)を解く代わりに以下の式を解けばよいという事になります。
 
 
 
【2次元での計算】※実際に解く方程式
 
ナビエストークス方程式
\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{5}
\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{6}
\end{align*}
 
ポアソン方程式※連続の式(非圧縮条件)
 
\begin{align*}
\frac{\partial^2 p}{\partial x^2}+\frac{\partial^2 p}{\partial y^2} = -\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{11}
\end{align*}
 

この3つの式を解けばよいという事になります。
で、3つ目の式がポアソン方程式です。

ポアソン方程式は数値的な解法としては有名なので、このような形に式変形しました。

速度にいては、ナビエストークス方程式を解きつつ、圧力に対してポアソン方程式を連立させながら解くという事になります。

だから、ここでポアソン方程式を学んでおこうというわけです。

ですが、
ポアソン方程式よりも簡単な右辺が0であるラプラス方程式を学びましょうというのが本記事の趣旨です。
 

ラプラス方程式の離散化

 

くどいですが、今回の記事では☟こちらの方程式を扱う・・・・

 

ポアソン方程式※連続の式(非圧縮条件)
 
\begin{align*}
\frac{\partial^2 p}{\partial x^2}+\frac{\partial^2 p}{\partial y^2} = b(x,y)\tag{12}
\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{13}
\end{align*}
 
 
わけではなく、\(b(x,y)=0\)とした場合に、
 
ラプラス方程式
 
 
\begin{align*}
\frac{\partial^2 p}{\partial x^2}+\frac{\partial^2 p}{\partial y^2} =0\tag{14}
\end{align*}
 
これを空間微分に対して、2次精度で離散化すると、
 
\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 – 2p_{i,j}^n + p_{i, j-1}^n}{\Delta y^2} = 0
\end{align*}
 
となります。
 
 
 
つまり、☟これを解けばよいという事です。
\begin{align*}
p_{i,j}^n = \frac{\Delta y^2(p_{i+1,j}^n+p_{i-1,j}^n)+\Delta x^2(p_{i,j+1}^n + p_{i,j-1}^n)}{2(\Delta x^2 + \Delta y^2)}\tag{15}
\end{align*}
 
 
もし、\(\Delta x = \Delta y\)とした場合は、
 
\begin{align*}
p_{i,j}^n = \frac{1}{4}\big(p_{i+1,j}^n+p_{i-1,j}^n+p_{i,j+1}^n + p_{i,j-1}^n\big)\tag{16}
\end{align*}
 
(16)式からわかるように、\((i,j)\)の格子点での圧力\(p_{i,j}\)はまわり4点の圧力の平均であることがわかります。
 
式を見ればわかる通り、ラプラス方程式は時間微分の項がありません。離散化した場合では、\(p^{n+1}\)がないってことです。

つまり、ラプラス方程式は平衡状態に達した場合における偏微分方程式です。

有限の空間で行われる数値計算においては、ラプラス方程式は境界条件によって空間分布の状態が決定するため「境界値問題」とも言われています。

要するに境界条件を決めたら物理量の空間分布は決まりますよってことです。

 
境界条件は、基本的には計算領域の一番端をどのようにするかで決まります。
 
今回行う境界条件は以下のように設定します。
 
【ディリクレ境界条件】
  • \(p=0\) at \(x=0\)
  • \(p=y\) at \(x=2\)

【ノイマン境界条件】
  • \(\frac{\partial p}{\partial y}=0\) at \(y=0, 1\)

境界条件を可視化すると☟このようになります。

 

この境界条件をもとにラプラス方程式を解くと、
 
\begin{align*}
p(x,y)=\frac{x}{4}-4\sum_{n=1,odd}^{\infty}\frac{1}{(n\pi)^2\sinh2n\pi}\sinh n\pi x\cos n\pi y\tag{17}
\end{align*}
となるとこちらのサイトでは言っていますが、
 
カマキリ

どうやったら(17)式になるんだろう・・・・って考え中です。

 
 

 
 
Twitterでも質問してみましたが、情報を得ることができず・・・・誰かわかったら教えてください。
(ラプラス方程式は変数分離をして、境界条件を課すことでできるというのが僕の中の定石ですが、上手くいきませんでした_(._.)_)
 
また別途、ラプラス方程式の解析解についてと数値計算の解法についての解説記事を買いたいと思います。
 
では、以下で実際にラプラス方程式をGoogle Colaboratoryを使ってPythonのコードを書きながら理解を深めていきたいと思います。
 

ラプラス方程式をPythonで実装

 
ひとつひとつ解説をしながらコードを書いていきます。
※関数にまとめておいた方がいい部分はあとでまとめたいと思います。
 
 
 
【ポイント】

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

 

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

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

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

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

 

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

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

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

そうは言っても数値計算では、初期の形状は確認しておいた方が良いですよね。

【ポイント】

Matplotlibを使って結果を可視化する。
 
コードを以下のように書きます。
 
【結果】
👆自分が想定していた境界条件になっていればOKです。
 
では、いよいよ上記で設定した境界条件をもとにラプラス方程式を解いてみましょう!
 
 
【ポイント】

\(p_{i,j}^n = \frac{1}{4}\big(p_{i+1,j}^n+p_{i-1,j}^n+p_{i,j+1}^n + p_{i,j-1}^n\big)\)について解くコードを実装します。
 
 

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

これで最終状態の圧力のプロファイルが完成しています。

何を行っているのかを「4×4」の格子状の点で簡単にして、回折します。

\begin{align*}
p_{i,j}^n = \frac{1}{4}\big(p_{i+1,j}^n+p_{i-1,j}^n+p_{i,j+1}^n + p_{i,j-1}^n\big)\tag{15}
\end{align*}

ですから、例えば(1,1)の格子点上の圧力は上記のように以下のようになります。

 

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

これを全ての格子点上で行うことを考えます。

今求めたいのは黒文字で書かれた圧力ですので、行列計算で求めることを考えます。

これで、未知の量について解くことができるのですが、

格子点数が多いとこのように行列計算で求めることは結構厳しいですよね。

この行列計算を行ってくれる方法としていろいろあるのですが、有名なのをいくつか挙げておきます。

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

などがあります。

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

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

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

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

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

数値計算の常識

数値計算の常識

伊理 正夫, 藤野 和建
2,750円(10/01 10:31時点)
発売日: 1985/06/03
Amazonの情報を掲載しています

その中でも今回は、反復法の中でも簡単な「ヤコビ法」で解くことにします。

 

 

収束判定は、

としています。

全ての格子点上の値の合計と前のステップでの値の合計とを比較して、「1e-4(0.01%)」以下の誤差になれば収束したと判定しています。

 

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

【ポイント】

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

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

【結果】

このような結果になりました(^^)/

 

関数にまとめる

 

今のままで良いのですが、関数にまとめてコードをきれいにまとめておきます。

☟初期状態を設定します。

次にラプラス方程式を解く部分を関数にまとめます。

ラプラス方程式を解いた後に、「return p」でpを返す関数を作成しました。

次に、結果を可視化する関数を用意します。

これで準備完成です(^^)/

初期の分布を見てみましょう。

【結果】

これは境界条件だけを可視化したもので、特に意味があるものではありませんね。

このあと、ラプラス方程式を解くことによって意味のある分布になります。

これでラプラス方程式を解いたことになります。
いちおう収束のために何回ループしたのかも計算しておきました。
「1132回」収束のためにループ計算をしています。

では、最終状態の結果を見てみます。

【結果】

 

まとめ

 

今回は、2次元のラプラス方程式をPythonで実装してみました。

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

次回の記事はこちら

おすすめの書籍紹介

 

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

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

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

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

 

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

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

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

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

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

流体の勉強をしたい方

 

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

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

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

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

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

数値流体の書籍

 

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

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

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

数値流体力学 第2版

数値流体力学 第2版

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

 

【プロフィール】

カマキリ
(^^)

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

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

 

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

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

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