gooブログはじめました!

写真付きで日記や趣味を書くならgooブログ

流体力学の計算を楽しむ

2016-02-07 07:57:42 | ブログ
 流体力学を少しかじったのが縁で、一様流中に静止する円柱について、その流線を計算し、図示してみたくなった。流線とは、時間の経過とともに流れが変化しない定常流の場合には、流れのみちすじを示す曲線のことである。

 一般に流体について数学的モデルをつくるときの手順は、次のようなものであろう。まず流体の流線についてのデータを実験的に採取し、そのデータを表すような関数を求める。このとき、実変数の関数よりも、複素数を変数とする複素関数を求めるのがよいということになる。

 ここで、複素関数として、複素速度ポテンシャルという物理量を考えるのが正規の手順となっている。

 複素速度ポテンシャルfは、iを虚数単位として、
     f=速度ポテンシャル+i×流れの関数
で与えられる。速度ポテンシャルは、流体に流速を生ぜしめるようなポテンシャル関数である。流れの関数は、流線を示す関数の集まりと考えられる。2次元の流れの場合、流れの関数(x,y)=定数が流線を表す関数である。

 速度ポテンシャルと流れの関数とは、コーシー・リーマンの関係式で結ばれているので、流れの関数が求められれば、速度ポテンシャルを計算によって求めることができる。あるいは、複素関数f(z)をローラン展開したものの虚数部分が流れの関数を表現するように近似できれば、その実数部分はこの関係式を満たすことになる。

 以下の説明は、縮まない完全流体(粘性のない流体)の2次元の渦なし運動に関するものである。

 一様な流れについて、x軸方向の速度をUとすると、複素速度ポテンシャルは、zを複素変数として
   f=Uz, U>0
で与えられる。z=x+iyであるから、
   速度ポテンシャル=Ux
   流れの関数=Uy
となる。

 次に、この一様流中の原点をその中心として静止する円柱を置いたときの複素速度ポテンシャルを考える。

 ここで、x-y平面上の点(a,0)に「わき出し」、点(-a,0)に「すいこみ」のある「二重わき出し」を考える。わき出した流体と同量の流体がすいこまれるとする。aを無限に小さくとると、その複素速度ポテンシャルfは、係数/zとなる。

 このような「二重わき出し」の流線は、y>0面内にあり原点に接する半径の異なる複数の円で構成される流線群と、y<0面にあり、この流線群と対称な流線群とより構成される。

 静止する円柱によってひきおこされる流れは、円柱の中心に「二重わき出し」をおいたときの流れと同等であるという。境界条件として、円柱の半径をaとし、一様流の速度をUとすると、その複素速度ポテンシャルは、Ua^2/zとなる(^はべき乗を表す)。

 速度ポテンシャルはラプラス方程式の解であるから、方程式は線形、流れの関数が解となるような方程式も線形であるから、上記2つの複素速度ポテンシャルは重ね合わせることができるので、結局、一様流中に静止する円柱に関する複素速度ポテンシャルは、
   f=Uz+Ua^2/z     (1)
となる。

 ここで、一様な流れは現実の流れに基づいているが、「二重わき出し」による流れは仮想的な流れであるから、両者を重ね合わせるということは、常識に反するような気がする。

 しかし、理論というものは、所詮は仮説であり、その結果がいかに現実をよく説明できるかだけが問題になるのだから、上記のような考えでよいのであろう。

 円柱から遠く離れると(Ua^2よりzがはるかに大きくなると)、(1)式のfはUzとなり、一様流を表す。

 簡単のためにa=1とする。z=x+iyとおいて、fの虚数部分をとり出すと、
   流れの関数=Uy-Uy/(x^2+y^2)
が得られる。充分遠方の流線の1つは、yの一定値をcとすると、Ucとしてよいから、この流れと同一のみちすじとなる流線を表す関数は、
   y-y/(x^2+y^2)=c
すなわち、
   y^3-cy^2+(x^2-1)y=cx^2
で与えられる。

 そうすると、c、すなわちyの一定値を設定し、そのcに関してxのキザミを順次変えてyについての3次方程式を解けばよい。これで1つの流線の軌跡が求められるから、cを変えて別の流線について3次方程式を作成し、同様の計算を繰り返せばよい。あるいは、xの一定値を設定し、そのxに関してcのキザミを順次変えてyについての3次方程式をつくってもよい。

 この関数の形をみると、絶対値が同一のxが正でも負でも同じ結果となるから、流線はy軸に関して対称となる。また、c>0,y>0とするこの関数の両辺に-1を掛けても成り立つから、c<0に変えると同一のyが負に変わるだけであり、流線はx軸に関しても対称である。

 c=0に設定した流れのみちすじは、x<-1ではx軸上を流れ、x=-1で上下に分離し、x=1まで円柱の円周にまとわりつき、x>1からは再びx軸上を流れることがわかる。

 (1)式の実数部分は速度ポテンシャルであるから、これをxで偏微分すると流速のx成分、yで偏微分すると流速のy成分が求められる。

 特に、円柱の表面の流速を求めるために、(1)式でz=aexp(it)とおくと(tはx軸からの偏角)、fは実数部分のみとなる。従って、速度ポテンシャルは、
   f=Uaexp(it)+Uaexp(-it)=2Uacost
により求められる。表面の流速は、
   v(t)=(1/a)(fのtに関する微分値)=-2Usint
で求められる。-がつくのは、tの進む方向が流速の方向と逆になるためである。t=+パイ/2(角度90度)または-パイ/2では、一様流の速度Uの2倍になる。

 円柱表面での流体の圧力は、ベルヌーイの定理からおよその見当がつけられる。すなわち、流速が0のとき圧力は最大であり、流速の増加とともに圧力が減少していく。

 そうすると、円柱表面の各点にかかる圧力の大きさはx軸に関して対称、y軸に関しても対称となる。圧力ベクトルの方向に関しても同様である。その結果、円柱表面全体が受ける圧力の総和は0である。言い換えれば、一様な流れに垂直におかれた円柱は、流体から力をうけないことになる。このことは、あきらかに日常の経験と反しているので、ダランベールの背理と呼ばれる。その原因は、流体が円柱表面と接するときに働く粘性力を無視したためである。

 Excelのグラフ機能を使って、流線の計算結果をグラフ表示することにした。

 流線は、y>0領域に6本、y<0領域に6本、合計12系列つくることにし、-3<x<3を設定する列にxのキザミ値を入力し、各系列に計算結果のデータを入力した。円柱に沿った流線のデータを入力するために、ExcelのSQRT関数を利用した。



 参考文献
 今井功著「物理テキストシリーズ 流体力学」(岩波書店)
 金丸隆志著「理系のためのExcelグラフ入門」(ブルーバックス)