● 流体力学

流体力学は水や空気のような連続体の運動を解析する学問であるが、 その問題の多くは方程式が複雑なため、代数的な厳密解を得ることが できない。 そのため以前では、飛行機や自動車を設計する際にまわりの空気の流れの 様子を調べるには模型による風洞実験に頼るしかなかった。 しかし近年では、計算機の性能の飛躍的向上が大型計算機によるそれらの 流れのシミュレーションを可能にし、工学や気象の分野で著しい成果をあげる ようになった。 流体の計算手法にはさまざまな特殊な職人芸があり、それらはあまりに 専門的である。しかし流体の計算の基本的な部分には、すべての物理の 分野の計算に必要な大事な計算技術として、場の偏微分方程式の差分解法と 座標変換がある。 この章では、流体力学の典型的な問題である「円柱周りの2次元流れ」を 計算機で解くことにより、これらの技術を学ぼう。

■ 流体力学の基礎方程式

▲ 非圧縮 Navier-Stokes方程式

流体の流れが十分遅い場合には、流体の密度の変化を無視することが できる。この場合、速さと長さの基準を適当にとることにより 流体の運動を支配する方程式に現れるパラメータは Reynolds number と 呼ばれるただ一つの定数 Re のみとなる。方程式は連続の方程式 (1)式とNavier-Stokes方程式(以降NS方程式と略す) (2)式である。

(1)
(2)

(3)

(4)
(5)

(6)

(7)

(8)

■ 円柱周りの2次元流れ

2次元流れの典型的な例として円柱周りの流れを考えよう。 つまり、x軸正の方向から負の方向への一様な流れが、座標原点に 置かれた円形の障害物によってどう乱されるかを考える。

▲ Reynolds numberの決定

(9)

方程式(7)(8)式はデカルト座標で記述されているが 障害物が円形なので極座標に書き改めると便利である。 右図のメッシュのように円の表面近くを細かく、表面から 遠く離れた所は荒くメッシュする次の極座標系 (ξ,θ)系で考えることにする。

(11)
(12)

(13)

(14)

(15)
(16)

▲ 初期条件と境界条件

この方程式を解くためには初期条件と境界条件が必要である。 初期条件としては、v_x=-1,v_y=0の一様流となる ψ,ωを使う。すなわち、(12)式と (15)式より

(17)

円の表面(ξ=0)での境界条件は、流体が粘性流体であることから 円の表面においては流れの速度は0である。 従って(12)式より ψは表面で一定値をとることがわかる。ψは定数項の自由度がまだあり、 そこで表面での値を0と定義しよう。 さらに (15)式式よりωの表面境界条件

(18)

(19)

方程式 (15)式は未知の場 ψ のLaplacianが 与えられた場ωを 満たすことを要請する方程式であり、このような形式の方程式を解くには 次の手段をとる。

(ξ,θ)平面上の格子点(i,j)でのψ(ξ,θ),ω(ξ,θ)の 値をψ_{i,j},ω_{i,j}と表して(15)式を差分方程式で 表してみよう。2階微分を2次精度で差分化すると次式となる。

(20)

(21)

このような手段を緩和法(relaxation method)と呼び、 ほどほどに方程式を満たすようになった 解をself consistent(自己無撞着)な解と呼ぶ。

■ プログラムの実装

必要な知識が揃ったので、プログラムを作り始めよう。ソースファイルは fluid.cc である。

#include "cip.h"
#include "nxgraph.h"

//---- definition of fundermental constants
//---- experimental settings
const double RE         = 100.0;           // Reynolds number
const double IRE        = 1.0/RE;          // inverse of RE
const double dT         = 0.0625;          // temporal step
const int    XI_NUM     = 24;              // radial cells
const int    TH_NUM     = 48;              // angular cells
const double dXI        = 0.125;           // radial step
const double dTH        = 2*M_PI/(TH_NUM-2(; // angular step
const double IdXI       = 1.0/dXI;         // inverse of dXI
const double IdTH       = 1.0/dTH;         // inverse of dTH

//---- index indicators
const int    XI_SUR     = 0;               // index of surface
const int    XI_FAR     = XI_NUM-1;        // index of farthest
const int    TH_STA     = 1;               // index of start angle
const int    TH_END     = TH_NUM-2;        // index of end angle

//---- graphics settings
const int    WIN_WIDTH  = 256;             // width of window
const int    WIN_HEIGHT = 256;             // height of window
const int    CYLINDER_R = 24;              // radius of cylinder
const double V_MAG      = 16.0;            // magnify factor of velocity

//---- prototype of functions
void InitFlow( void );
void SetBorder( void );
void ModifyPsi( void );
void UpdateOmg( void );
void DrawFlow( void );

//---- field variables
double Psi[XI_NUM][TH_NUM];  // Flow Function field
double Omg[XI_NUM][TH_NUM];  // Vorticity field 

class ExpTable
{
  double expo[4*XI_NUM];
public:
  ExpTable( void ){
    for( int i=-2*XI_NUM ; i<2*XI_NUM ; i++ ){
      expo[2*XI_NUM+i] = exp( i*dXI );
    }
  }
  inline double operator () ( int i ){
    return( expo[2*XI_NUM+i] );
  }
} Exp;

class SinTable
{
  double sine[TH_NUM];
public:
  SinTable( void ){
    for( int j=TH_STA ; j<=TH_END ; j++ ){
      sine[j] = sin( j*dTH );
    }
  }
  inline double operator () ( int j ){
    return( sine[j] );
  }
} Sin;

class CosTable
{
  double cosi[TH_NUM];
public:
  CosTable( void ){
    for( int j=TH_STA ; j<=TH_END ; j++ ){
      cosi[j] = cos( j*dTH );
    }
  }
  inline double operator () ( int j ){
    return( cosi[j] );
  }
} Cos;

class XTable
{
  int x[XI_NUM][TH_NUM];
public:
  XTable( void ){
    for( int i=XI_SUR ; i<=XI_FAR ; i++ ){
      for( int j=TH_STA ; j<=TH_END ; j++ ){
	x[i][j] = WIN_WIDTH/2  + int(CYLINDER_R*Exp(i)*Cos(j));
      }
    }
  }
  inline int operator () ( int i, int j ){
    return( x[i][j] );
  }
} X;

class YTable
{
  int y[XI_NUM][TH_NUM];
public:
  YTable( void ){
    for( int i=XI_SUR ; i<=XI_FAR ; i++ ){
      for( int j=TH_STA ; j<=TH_END ; j++ ){
	y[i][j] = WIN_HEIGHT/2 + int(CYLINDER_R*Exp(i)*Sin(j));
      }
    }
  }
  inline int operator () ( int i, int j ){
    return( y[i][j] );
  }
} Y;

//---- intialize flow as uniformal flow
void InitFlow( void )
{
  for( int i=XI_SUR ; i<=XI_FAR ; i++ ){
    for( int j=TH_STA ; j<=TH_END ; j++ ){
      Psi[i][j] = -Exp(i)*Sin(j);    Omg[i][j] = 0.0;
    }
  }
}

//---- main function
int main( void )
{
  XEvent ev;
  NXOpenWindow( "2D flow around a cylinder", WIN_WIDTH, WIN_HEIGHT );
  NXEnable( NX_DOUBLEBUFFER );

  InitFlow();           // inttialize flow as uniformal flow
  do{                   // evolution loop
    SetBorder();        // set border condtions at surface and far away
    ModifyPsi();        // solve phi by eq13
    UpdateOmg();        // calculate next omg by eq10
    DrawFlow();         // draw velocity filed

    NXCheckEvent( NX_NOWAIT, ev );
  }while( ev.type != KeyPress );

  NXCloseWindow();

  return(0);
}

//---- set border conditions on surface and far away
void SetBorder( void )
{
  for( int j=TH_STA ; j<=TH_END ; j++ ){    // anguler loop
    Psi[XI_SUR][j] = 0.0;                    Omg[XI_SUR][j] = (-2.0*IdXI*IdXI) * Psi[XI_SUR+1][j];    // surface
    Psi[XI_FAR][j] = -Exp(XI_FAR)*Sin(j);    Omg[XI_FAR][j] = 0.0;                                    // far away
  }
}

//---- derive flow function field (psi) from Laplace equation
void ModifyPsi( void )
{
  int i, j, loop=8;
  while( loop-- ){
    for( i=XI_SUR+1 ; i<=XI_FAR-1 ; i++ ){  // radial loop
      Psi[i][TH_STA-1] = Psi[i][TH_END];      Psi[i][TH_END+1] = Psi[i][TH_STA];    // cyclic condition
      for( j=TH_STA ; j<=TH_END ; j++ ){    // anguler loop
        Psi[i][j] = (0.5/(IdXI*IdXI+IdTH*IdTH))*(
          + (IdXI*IdXI)*(Psi[i+1][j]+Psi[i-1][j]) + (IdTH*IdTH)*(Psi[i][j+1]+Psi[i][j-1]) + Exp(2*i)*Omg[i][j] );
      }
    }
  }
}

//---- update vorticity field (omg)
void UpdateOmg( void )
{
  int i, j;
  for( i=XI_SUR+1 ; i<=XI_FAR-1 ; i++ ){  // radial loop
    Omg[i][TH_STA-1] = Omg[i][TH_END];    // cyclic condition
    Omg[i][TH_END+1] = Omg[i][TH_STA];
    for( j=TH_STA ; j<=TH_END ; j++ ){    // anguler loop
      Omg[i][j] +=
          dT*Exp(-2*i)*( + (IdXI*0.5) * (Psi[i+1][j]-Psi[i-1][j]) * (IdTH*0.5) * (Omg[i][j+1]-Omg[i][j-1])
                         - (IdTH*0.5) * (Psi[i][j+1]-Psi[i][j-1]) * (IdXI*0.5) * (Omg[i+1][j]-Omg[i-1][j]) )
        + dT*Exp(-2*i)*IRE*( + (IdXI*IdXI) * (Omg[i+1][j]-2*Omg[i][j]+Omg[i-1][j])
                             + (IdTH*IdTH) * (Omg[i][j+1]-2*Omg[i][j]+Omg[i][j-1]) );
    }
  }
}

//---- derive velocity field and draw them 
void DrawFlow( void )
{
  int i, j;
  double Vxi, Vth, Vx, Vy;
  NXClearWindow();
  for( i=XI_SUR+1 ; i<=XI_FAR-1 ; i++ ){  // radial loop
    for( j=TH_STA ; j<=TH_END ; j++ ){    // anguler loop
      Vxi = (+0.5*IdTH) * Exp(-i)*(Psi[i][j+1]-Psi[i][j-1]);
      Vth = (-0.5*IdXI) * Exp(-i)*(Psi[i+1][j]-Psi[i-1][j]);
      Vx  = Vxi*Cos(j) - Vth*Sin(j);
      Vy  = Vxi*Sin(j) + Vth*Cos(j);

      NXDrawMoveto( X(i,j), Y(i,j) );
      NXDrawLinerel( int(V_MAG*Vx), int(V_MAG*Vy) );
    }
  }
  NXFlush();
}

(22)
(23)

(24)

(25)

(26)

これに対して次のような差分化を陰的(implicit)な差分化と呼ぶ。

(27)

(28)

なお、(27)式で u_j^{n+1} を計算するためには、 すべての j の u_j^{n+1} に関する(27)式の 連立方程式を解く計算を行うか、緩和法を用いる。

偏微分方程式の差分解法は非常に奥が深いので、興味のある読者には 専門書を参照されることを勧める。 本書で強調しておくことは、計算機を安易に信じて単純に計算を行うことの 危険性と、その危険性を回避するための人間の英知と努力があることである。