理系暇人たちのホワイトボード

　暇なので、gnuplotを使って遊んでました。

　今回は、複素速度ポテンシャルに基づいた、半無限物体を通る流れを描画しようと思い、やってみたのでブログに挙げてみます。

前提事項
- 1. 複素速度ポテンシャル
- 2. 流れ場の描画に用いる事項
半無限物体を通る流れについて考える

前提事項

1. 複素速度ポテンシャル

　本題に入る前に、流体力学における複素速度ポテンシャルとは何たるかを軽く解説したいと思います。

　そもそも、流れ場が複素速度ポテンシャルで表すことができるのは特定の場合に限られており、それは「非圧縮かつ渦なし」の場合である、と認識しています。そのような場合について、どういう形で流れ場が複素速度ポテンシャルで表すことができるのかを以下で実際に数式を用いて見ていきましょう…

ここからは、2次元流れを考えていきます。

まず、非圧縮であることを考えると流体力学の重要な基本式「連続の式」は

$\frac{\partial u}{\partial x} + \frac{\partial v}{\partial y} = 0　\tag{1-a}$

つまり

$div \, \boldsymbol{v} = 0 \tag{1-b}$

と、表されます。簡単な形ですね。ちなみに

$\boldsymbol{v} = \left( u, v \right) \tag{2}$

は流速のことです。

更に渦なしであることを考慮していきます。渦なしの場合、渦度 $\boldsymbol{\omega}$ とすると

$\boldsymbol{\omega} = rot\, \boldsymbol{v} = \boldsymbol{0} \tag{3}$

となりますね。ここでベクトル解析の恒等式：

$rot\, grad\, \phi = \boldsymbol{0} \tag{4}$

を思い出して、(3)に適用すれば、流速 $\boldsymbol{v}$ は

$\boldsymbol{v} = grad\, \phi \tag{5}$

と書き表すことができます。つまり $\phi$ は速度ポテンシャルということになりますね。

さて、再び連続の式(1-a)式に戻りましょう。この式を満たせばよいのですから、 $u, v$ を

$u = \frac{\partial \psi}{\partial y}, \,\, v = - \frac{\partial \psi}{\partial x} \tag{6}$

とすれば(1-a)は常に満足することがわかると思います。ここで導入した $\psi$ は流体力学では「流れの関数」と呼ばれるものであり、簡単に物理的意味を説明すると

「流れの関数 $\psi$ をある方向に微分すると、その方向から直角に右に回った方向の速度が得られるようなもの」

といったところでしょうかね。詳しい説明は流体の教科書に載っているはずなのでここでは割愛させてくださいｗ（後のほうでもう少し説明を加えますが。）

以上のことを一旦まとめておくと、流速は

$u = \frac{\partial \phi}{\partial x} = \frac{\partial \psi}{\partial y} \tag{7}$

$y = \frac{\partial \phi}{\partial y} = - \frac{\partial \psi}{\partial x} \tag{8}$

のように、スカラー値 $\phi, \psi$ を用いて表されます。

ところで、いきなり $\phi, \psi$ などといった訳のわからん変数（ $\phi$ は速度ポテンシャル、 $\psi$ は流れの関数なんですけれどもねｗ)を導入して何がしたいんや…と突っ込みたくなるところですが

ここでようやく複素速度ポテンシャルの登場です…！

その意味を述べる前に、複素速度ポテンシャル $f$ はこれまでに出てきた変数を用いて

$f = \phi + i \psi \tag{9}$

と表される関数です。大学で数学を勉強している人なら気づいたかもしれませんが、(9)式のように複素速度ポテンシャルに登場する $\phi, \psi$ は(7), (8)式からわかるように「コーシー・リーマンの関係式」を満たしています。(7), (8)式は常に成り立つようにしてあるので、このとき $f$ は複素平面全体で正則となり、微分可能となっています。

色々書きましたが、結局のところ複素速度ポテンシャルを使うと何が嬉しいかというと

「二次元の流れ場を複素関数論を用いて簡単に調べることができる」

ことだと思っています。（実際はどうなんでしょうかw）

2. 流れ場の描画に用いる事項

　これまでで複素速度ポテンシャルを用いて流れを表せるということを説明しておいて言うのもなんですが、実際に私がgnuplotを用いて流れを描画する際には複素速度ポテンシャル自体は用いませんでした…ｗ

　というのも、流れを視覚化するのに流線を考えるのですが、流線の様子は前述の「流れの関数」のみを考えることで調べることができるからです。それについて以下で少し説明したいと思います。

　流れの関数 $\psi$ の物理的意味は先ほど説明したように、ある方向に微分するとその方向から直角に右に回った方向の速度が得られるような関数でしたね。したがって、流れと同じ向き、すなわち流速と同じ向きの $\psi$ の方向微分は0ということになるので

$\frac{d \psi}{d s} = 0 \tag{10}$

ここで $d s$ は流線に沿った微小な線素です。すなわち

$\psi (x, y) = const. (流線に沿って）\tag{11}$

が成り立ちます。したがって、(11)式をx-y平面上に描けば、流れの様子が視覚的にわかるといった寸法です。

　以上、やや長い前置きでしたが、流れの様子を描く準備が整ったのでこの後で実際に半無限物体を通る流れを描画していきます。

半無限物体を通る流れについて考える

1. 流れの関数

　半無限物体を通る流れの複素速度ポテンシャルは

$f = U z + m\,{\rm ln}z \ \ \ (U, m : const.) \tag{12}$

のように一様流のポテンシャルと湧き出しのポテンシャルの組み合わせによって表すことができます（以下の結果から実際にそうなることが見て分かると思います）。 $U$ は一様流の流速、 $m$ は湧き出しの強さを表します。

　ここで $z = re^{i\theta}$ と極座標で表すと、流れの関数 $\psi$ は

$\psi (r, \theta)= U r \,{\rm sin}\, \theta + m \theta \tag{13}$

と表すことができます。(13)式を(11)式の左辺に持っていき、(11)式の右辺の定数を小刻みに変えて描画することにより、半無限物体周辺の流れ場の様子を描画していきます。

2. 結果

　以下に、(13)式の $U, m$ をそれぞれ変えて流れの様子を描画した結果を示します。

(i) U = 0.1, m = 1

f:id:sKrut:20200919114932p:plain — 図1. U = 0.1, m = 1 の場合

分かりやすくするため、半無限物体の表面に相当する流線（よどみ点を通る流線）を太い黒線で表し、内部を青線、実際に物体の周りを流れると考えられる部分の流れを赤線で表しました（以下同様）。凝視すると目がバグりそうｗ

　ちなみに、物体の先端付近を拡大してみると以下のようになりました：

f:id:sKrut:20200919115522p:plain — 図2. U = 0.1, m = 1 の場合（拡大図）

前に述べたように、図の右側に向かう一様流の中に湧き出し点があり、それらの組み合わせにより半無限物体の周りの流れができていることが確認できると思います。

(ii) U = 0.1, m = 5

f:id:sKrut:20200919184500p:plain — 図3. U = 0.1, m = 5 の場合

　m = 5に変えた場合、図1と比較して、物体のy方向の幅が大きくなり、周辺の流れが大きく迂回するようになっていることがわかります。図1の場合と比べてmの値が大きく湧き出しの強さが強くなっているので、この結果はイメージと合致していますね。

(iii) U = 0.5, m = 1

f:id:sKrut:20200919185006p:plain — 図4. U = 0.5, m = 1 の場合

　今度はU = 5と変更しており、図1の場合と比べて流速が大きいような状況です。物体の形は棒のように細長くなっています。

3. 感想など

　青い線は、(11)式が $r \gt 0$ に対して $\theta$ が1つしか解を持たないので $r \gt 0$ の範囲で $r$ を小刻みに変えていき、それぞれの $r$ に対してNewton法を用いて解 $\theta$ を求めましたが、外側の赤い線は同じ $r \gt 0$ に対して $\theta$ が2つの解を持ちうるので、 $\theta$ を $0 \rightarrow \pi$ で変化させてそれぞれについて $r$ を求めました。このような感じで何となくソースコードを書いて数値計算しただけなので、数値計算の知識などはほとんど用いていません。 $\theta$ の範囲を $0 \rightarrow 2\pi$ として計算してみたところ描画結果が支離滅裂になったので、数値計算を勉強を進めてその原因などを調べていきたいと思いました。