« ThinkVantage 指紋認証ユーティリティー | トップページ | NTEmacs + MinGW + MSYSでWindows上にC言語開発環境を構築してみる - その1 »

2008/10/05

2次元キャビティ流れ

2次元の正方形キャビティ内の流れを「流れ関数-渦度法」で変形した支配方程式を差分を使って解く数値計算プログラムを作っている。まだ未完成でいろいろ修正点もあるけど一応それらしい解が得られている。

Psi2_2 

この図は流れ関数の等高線をプロットしたもの。

大雑把に言えば、このプログラムは移流拡散方程式(Navieor-Stoks方程式のrotをとったもの)とポアソン方程式(渦度と流れ関数の関係式)を解いているだけだけど、(キャビティ流れの)理論的な部分を完全に理解していないので他人に説明できるようなレベルじゃない。(この状態だと、"このようにすると何故か動いてくれる"といった説明しかできない)

あと未完成だけどソースコード:「cavity.c」

/*
  [2008]
   2次元 Cavity Flow Simulation
*/

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

/* Cavity Parameters */
#define L  1.0    // 1辺の長さ
#define N  100    // メッシュの数
#define DX (L/N)  // メッシュ間隔(x)
#define DY (L/N)  // メッシュ間隔(y)

/* Fluid Parameters */
#define U  1.0       // 流速
#define RE 100       // レイノルズ数
#define NU (U*L/RE)  // 動粘性率

/* CFL Condition */
#define CFL     0.3                                    // 安全係数alpha
#define DT_ADV  (1.0/(U/DX))                           // 移流時間
#define DT_DIFF (1.0/(2*NU*(1.0/DX/DX + 1.0/DY/DY)))   // 拡散時間
#define DT      (CFL*(1.0/(1.0/DT_ADV + 1.0/DT_DIFF))) // 時間刻み

/* Calc Step */
#define TMAX 10000   // 計算回数

/* Poisson : Gauss-Sidel Parameters */
#define LOOP_MAX (N*N/4)   // ガウスザイデルの計算回数上限
#define EPS      (1.0e-4)  // 精度


/* Define New Data Type */
typedef double Array[N+1][N+1];


/* Initialize Array Type */
void initArray(Array var)
{
  int i, j;
  for(i = 0; i <= N; i++){
    for(j = 0; j <= N; j++)
      var[i][j] = 0.0;
  }
}


/* Omega's Boundary Conditions */
void setBoundary_omega(Array omega, Array psi)
{
  int ix, iy;
  for(ix = 0; ix <= N; ix++){
    omega[ix][0] = -2.0*(psi[ix][0+1]     )/DY/DY;   // bottom
    omega[ix][N] = -2.0*(psi[ix][N-1]+U*DY)/DY/DY;   // top
  }
  for(iy = 0; iy <= N; iy++){
    omega[0][iy] = -2.0*psi[0+1][iy]/DX/DX;          // left
    omega[N][iy] = -2.0*psi[N-1][iy]/DX/DX;          // right
  }
}

/* Psi's Boundary Conditions */
void setBoundary_psi(Array psi)
{
  int ix, iy;
  for(ix = 0; ix <= N; ix++){
    psi[ix][0] = 0.0;   // bottom
    psi[ix][N] = 0.0;   // top
  }
  for(iy = 0; iy <= N; iy++){
    psi[0][iy] = 0.0;   // left
    psi[N][iy] = 0.0;   // right
  }
}
 

/* Operator */
double advection(Array omega, Array psi, int x, int y)
{
  return ( (psi[x+1][y]-psi[x-1][y])*(omega[x][y+1]-omega[x][y-1])
	   -(psi[x][y+1]-psi[x][y-1])*(omega[x+1][y]-omega[x-1][y])
	   )/4.0/DX/DY;
}

double laplacian(Array omega, int x, int y)
{
  return ( (omega[x+1][y]-2.0*omega[x][y]+omega[x-1][y])/DX/DX
	   +(omega[x][y+1]-2.0*omega[x][y]+omega[x][y-1])/DY/DY );
}

/* Solve Omega */
void solveOmega(Array omega, Array psi)
{
  int ix, iy;
  static Array omega_new;

  for(ix = 1; ix < N; ix++){
    for(iy = 1; iy < N; iy++){
      omega_new[ix][iy] = omega[ix][iy]
	+DT*(advection(omega, psi, ix, iy) + NU*laplacian(omega, ix, iy));
    }
  }

  for(ix = 1; ix < N; ix++){
    for(iy = 1; iy < N; iy++){
      omega[ix][iy] = omega_new[ix][iy];
    }
  }
  
}


/* Solve Psi using Gauss-Sidel Method */
void solvePsi(Array psi, Array omega)
{
  int    ix, iy, iter;
  double psiMax = 1.0e-7;
  
  for(iter = 1; iter <= LOOP_MAX; iter++){
    double errMax = 0.0;
    for(ix = 1; ix < N; ix++){
      for(iy = 1; iy < N; iy++){
	double err     = 0.0;
	double psi_old = psi[ix][iy];
	double psi_new = (DX*DX*DY*DY*omega[ix][iy]
			  +DY*DY*psi[ix+1][iy  ]
			  +DY*DY*psi[ix-1][iy  ]
			  +DX*DX*psi[ix  ][iy+1]
			  +DX*DX*psi[ix  ][iy-1])/(2.0*(DX*DX+DY*DY));
	psi[ix][iy] = psi_new;
	if(psiMax < fabs(psi_new))
	  psiMax = psi_new;
	err = (fabs(psi_new - psi_old))/psiMax;
	if(errMax < err)
	  errMax = err;
      }
    }
    if(errMax < EPS) break;
  }
}


/* Output Data */
void output(FILE *fp, Array psi)
{
  int ix, iy;
  for(ix = 0; ix <= N; ix++){
    for(iy = 0; iy <= N; iy++){
      fprintf(fp, "%g %g %g\n", ix*DX, iy*DY, psi[ix][iy]);
    }
    fprintf(fp, "\n");
  }
}

/* Plot Data */
void plot_psi(FILE *gp, Array psi)
{
  int ix, iy;
  fprintf(gp, "set contour\n");
  fprintf(gp, "unset surface\n");
  fprintf(gp, "splot '-' w l\n");
  for(ix = 0; ix <= N; ix++){
    for(iy = 0; iy <= N; iy++){
      fprintf(gp, "%g %g %g\n", ix*DX, iy*DY, psi[ix][iy]);
    }
    fprintf(gp, "\n");
  }
  fprintf(gp, "end\n");
  fflush(gp);
}

/* MAIN */
int main(void)
{
  static Array psi;     // Stream Function
  static Array omega;   // Vorticity
  int  it;
  FILE *fp = fopen("output.dat", "w");
  FILE *gp = popen("gnuplot -persist", "w");

  /* ===Initialize=== */
  initArray(psi);
  initArray(omega);

  /* ===Set Boundary Condition=== */
  setBoundary_psi(psi);
  setBoundary_omega(omega, psi);

  /* ===Solve Omega & Psi=== */
  for(it = 1; it <= TMAX; it++){
    solveOmega(omega, psi);
    solvePsi(psi, omega);
    setBoundary_psi(psi);
    setBoundary_omega(omega, psi);
    printf("Iterator:%d\r", it);
  }
  printf("\n");

  /* ===Output Data=== */
  output(fp, psi);
  fclose(fp);

  /* ===Plot Data=== */
  plot_psi(gp, psi);
  pclose(gp);
  
  return 0;
}

|

« ThinkVantage 指紋認証ユーティリティー | トップページ | NTEmacs + MinGW + MSYSでWindows上にC言語開発環境を構築してみる - その1 »

C」カテゴリの記事

Numerical Computing」カテゴリの記事

コメント

コメントを書く



(ウェブ上には掲載しません)


コメントは記事投稿者が公開するまで表示されません。



トラックバック

この記事のトラックバックURL:
http://app.cocolog-nifty.com/t/trackback/222903/42707209

この記事へのトラックバック一覧です: 2次元キャビティ流れ:

« ThinkVantage 指紋認証ユーティリティー | トップページ | NTEmacs + MinGW + MSYSでWindows上にC言語開発環境を構築してみる - その1 »