2次元キャビティ流れ
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; }
| 固定リンク
「C」カテゴリの記事
- CG法(2008.12.09)
- NTEmacs + MinGW + MSYSでWindows上にC言語開発環境を構築してみる - その5(2008.10.10)
- NTEmacs + MinGW + MSYSでWindows上にC言語開発環境を構築してみる - その4(2008.10.09)
- NTEmacs + MinGW + MSYSでWindows上にC言語開発環境を構築してみる - その3(2008.10.08)
- NTEmacs + MinGW + MSYSでWindows上にC言語開発環境を構築してみる - その2(2008.10.07)
「Numerical Computing」カテゴリの記事
- BiCGSTAB法(2008.12.23)
- CG法(2008.12.09)
- 2次元キャビティ流れ(2008.10.05)
- Box-Muller法で発生させた乱数で正規分布を描く(2008.09.06)
- グラフの表示を遅延(2008.07.14)
この記事へのコメントは終了しました。
コメント