Microsoft PowerPoint - 先端GPGPUシミュレーション工学特論(web).pptx

数値流体力学への応用

（支配方程式，

CPUプログラム）

今回の内容



支配方程式



_{Taylor‐Green渦}

数値流体力学



数値計算を利用して，流体の挙動を計算



_{Computational Fluid Dynamics（略してCFD）}



計算機の性能向上に伴い，必要不可欠な設計ツー

ルとなっている



流体を取り扱う機器の性能評価



流体中を移動する物体が受ける抵抗の評価など



提案された当初は

_C

_O

_L

_O

_R

_F

_U

_L

_{Fluid Dynamicsと}

揶揄されていた

支配方程式



_{2次元直交(x–y)座標系}



x方向速度をu, y方向速度をvと記述



非圧縮性流れ



流体の密度



が変化しない



水は圧縮性



一般にマッハ0.3以下では非圧縮と見なす



粘性流れ



流体の粘り気により運動が妨げられる

支配方程式



連続の式およびNavier‐Stokes方程式



質量と運動量の保存式に対応



運動エネルギの保存は運動量保存式が兼ねる



x, y方向速度u, v，圧力pの3変数を連立

0













y

v

x

u























































































































2 2 2 2 2 2 2 2

1

1

y

v

x

v

y

p

y

v

v

x

v

u

t

v

y

u

x

u

x

p

y

u

v

x

u

u

t

u









連続の式

（質量保存式）

Navier‐Stokes方程式

（運動量保存式）

支配方程式



_{Navier‐Stokes方程式の回転を計算}



vに関する式をxで偏微分した式から，uに関する式をyで偏微

分した式を引く



時間，空間的に連続であることを考慮して微分の順序を交換



密度一定のため衝撃波のような不連続は発生しない



















































































































































2 2 2 2 2 2 2 2

1

1

y

u

x

u

x

p

y

u

v

x

u

u

t

u

y

y

v

x

v

y

p

y

v

v

x

v

u

t

v

x









支配方程式



渦度（Vorticity）



流体粒子の回転



流れ関数（Streamfunction）



二つの流れ関数の等高線（流線）間を通過する流量

y

u

x

v















v

x

u

y



















x

v





y

u







1



2



v

x











y

u









支配方程式



渦度方程式



渦度の移流・拡散方程式



流れ関数のPoisson方程式



流れ関数の定義式を渦度の定義式に代入



変数の数が流れ関数，渦度の2個に低減













































2 2 2 2

y

x

y

v

x

u

t

































2 2 2 2

y

x





























































y

y

x

x

y

u

x

v







渦度の移流

渦度の拡散

計算手順

1.

渦度方程式を計算して渦度を求める

2.

渦度を基に流れ関数のPoisson方程式を解いて

流れ関数を求める

3.

流れ関数の定義式から速度を求める

4.

1に戻って計算を繰り返す

渦度方程式の離散化



時間に前進差分，空間に中心差分を適用







































         2 1 , , 1 , 2 , 1 , , 1 1 , 1 , , , 1 , 1 , , 1 ,

2

2

2

2

Δy

Δx

Δy

v

Δx

u

Δt

n j i n j i n j i n j i n j i n j i n j i n j i n j i n j i n j i n j i n j i n j i

































































































         2 1 , , 1 , 2 , 1 , , 1 1 , 1 , , , 1 , 1 , , 1 ,

2

2

2

2

Δy

Δx

Δt

Δy

v

Δx

u

Δt

n j i n j i n j i n j i n j i n j i n j i n j i n j i n j i n j i n j i n j i n j i



























流れ関数の

Poisson方程式の離散化



中心差分による離散化



時刻

n+1における流れ関数と渦度の関係



流れ関数の定義式の離散化

1 , 2 1 1 , 1 , 1 1 , 2 1 , 1 1 , 1 , 1

2

  

2

  _  





_





_

_

n j i n j i n j i n j i n j i n j i n j i

Δy

Δx















Δx

v

Δy

u

n j i n j i n j i n j i n j i n j i

2

2

1 , 1 1 , 1 1 , 1 1 , 1 1 , 1 ,          



















計算条件



_{Taylor‐Green渦}



セル状の渦が流体の粘性によ

り減衰



移流と圧力勾配が釣り合うた

め，粘性の影響しか現れない



計算条件



拡散方程式の安定条件を満た

すように諸条件を設定



境界条件



流れ関数，渦度とも0

2



2



計算条件



_{Taylor‐Green渦}



セル状の渦が流体の粘性によ

り減衰



移流と圧力勾配が釣り合うた

め，粘性の影響しか現れない



計算条件



拡散方程式の安定条件を満た

すように諸条件を設定



境界条件



流れ関数，渦度とも0

2



2





2

−2

y

xsin

sin

2





#include<stdio.h> #include<stdlib.h> #include<math.h> #define Lx  (2.0*M_PI) #define Ly  (2.0*M_PI) #define Nx  101 #define Ny  Nx #define Nbytes (Nx*Ny*sizeof(double)) #define dx  (Lx/(Nx‐1)) #define dy  (Ly/(Ny‐1)) #define Kvisc (0.1) #define dx2  (dx*2.0) #define dy2  (dy*2.0) #define dxdx (dx*dx) #define dydy (dy*dy) #define dxdxdydy (dxdx*dydy) #define dxdy2 (2.0*(dxdx+dydy)) #define dt (0.01) #define Nt (int)(1.0/dt) #define ERR_TOL 1e‐12 #define Accel 1.7

パラメータ設定

taylor_green.c

int main(void){ double *vrtx[2],*stmf,*u,*v; int curr=0,next=1; vrtx[curr] = (double *)malloc(Nbytes); stmf = (double *)malloc(Nbytes); u      = (double *)malloc(Nbytes); v      = (double *)malloc(Nbytes); vrtx[next] = (double *)malloc(Nbytes); //Poisson方程式を満たす初期状態を計算 init(vrtx[curr],stmf,u,v,vrtx[next]); computeStreamfunction_rbsor (stmf,vrtx[curr]);//流れ関数 computeVelocity(stmf,u,v);      //速度 //時間積分 for(int n=1;n<=Nt;n++){ //渦度の計算 computeVorticity (vrtx[curr],stmf,u,v,vrtx[next]); //流れ関数の計算 computeStreamfunction_rbsor (stmf,vrtx[next]); //速度の計算 computeVelocity(stmf,u,v); curr = next; next = 1‐curr; } free(vrtx[curr]); free(stmf); free(u); free(v); free(vrtx[next]); return 0; }

メイン

taylor_green.c

void init(double *vrtx, double *stmf, double *u, double *v, double *vrtx_next){ int i,j; for(j=0;j<Ny;j++){ for(i=0;i<Nx;i++){ double x = (double)i*dx; double y = (double)j*dy; vrtx[i+Nx*j] = 2.0*sin(2.0*M_PI*x/Lx)*sin(2.0*M_PI*y/Ly); stmf[i+Nx*j] = 0.0; v[i+Nx*j] = 0.0; u[i+Nx*j] = 0.0; vrtx_next[i+Nx*j] = 0.0; } } }

初期化

taylor_green.c

void computeVorticity (double *vrtx, double *stmf, double *u, double *v,double *vrtx_next){ int i,j; int ij,im1j,ip1j,ijm1,ijp1; double conv,visc; for(j=1;j<Ny‐1;j++) for(i=1;i<Nx‐1;i++){ ij = i +Nx*j; im1j = i‐1+Nx*j; ip1j = i+1+Nx*j; ijm1 = i +Nx*(j‐1); ijp1 = i +Nx*(j+1); //移流項の計算 conv = u[ij]*(vrtx[ip1j]‐vrtx[im1j])/dx2 + v[ij]*(vrtx[ijp1]‐vrtx[ijm1])/dy2; //粘性項の計算 visc = Kvisc*( (vrtx[im1j]‐2.0*vrtx[ij]+vrtx[ip1j])/dxdx +(vrtx[ijm1]‐2.0*vrtx[ij]+vrtx[ijp1])/dydy); //時間積分 vrtx_next[ij] = vrtx[ij] + dt*(‐conv+visc); } }

渦度方程式の計算（時間積分）

taylor_green.c

void computeStreamfunction_rbsor (double *stmf, double *vrtx){ int ite_SOR=0; double err_n,err_d,err_relative; double d_s; int ij,ip1j,im1j,ijm1,ijp1; do{ err_n = 0.0; err_d = 0.0; for(int j=1;j<Ny‐1;j++){ for(int i=1+j%2;i<Nx‐1;i+=2){ ij = i+Nx*j; ip1j = (i+1)+Nx*(j  ); im1j = (i‐1)+Nx*(j  ); ijp1 = (i )+Nx*(j+1); ijm1 = (i )+Nx*(j‐1); d_s =( (stmf[im1j]+stmf[ip1j])/dxdx +(stmf[ijm1]+stmf[ijp1])/dydy ‐(‐vrtx[ij]) //右辺が‐ωなのでマイナス付き )*dxdxdydy/dxdy2 ‐stmf[ij]; stmf[ij] += Accel*d_s; err_n += d_s*d_s; err_d += stmf[ij]*stmf[ij]; } }

流れ関数の計算（

RB-SOR法）

taylor_green.c

for(int j=1;j<Ny‐1;j++){ for(int i=2‐j%2;i<Nx‐1;i+=2){ ij = i+Nx*j; ip1j = (i+1)+Nx*(j  ); im1j = (i‐1)+Nx*(j  ); ijp1 = (i )+Nx*(j+1); ijm1 = (i )+Nx*(j‐1); d_s =( (stmf[im1j]+stmf[ip1j])/dxdx +(stmf[ijm1]+stmf[ijp1])/dydy ‐(‐vrtx[ij]) //右辺が‐ωなのでマイナス付き )*dxdxdydy/dxdy2 ‐stmf[ij]; stmf[ij] += Accel*d_s; err_n += d_s*d_s; err_d += stmf[ij]*stmf[ij]; } } if(err_d<1e‐20)err_d=1.0; err_relative = sqrt(err_n/err_d); ite_SOR++; }while(err_relative > ERR_TOL); }

流れ関数の計算（

RB-SOR法）

taylor_green.c

void computeStreamfunction_cg (double *stmf, double *vrtx){ double err; int ij,ip1j,im1j,ijm1,ijp1; double *p,*r,*Ap; double alph,beta,rr,pAp,bb; p  = (double *)malloc(Nbytes); r  = (double *)malloc(Nbytes); Ap = (double *)malloc(Nbytes); for(int j=0;j<Ny;j++) for(int i=0;i<Nx;i++){ ij = i+Nx*j; p [ij] = 0.0; r [ij] = 0.0; Ap[ij] = 0.0; } alph=0.0; beta=0.0; bb = 0.0; rr = 0.0; for(int j=1;j<Ny‐1;j++) for(int i=1;i<Nx‐1;i++){ ij = i+Nx*j; ip1j = (i+1)+Nx*(j  ); im1j = (i‐1)+Nx*(j  ); ijp1 = (i )+Nx*(j+1); ijm1 = (i )+Nx*(j‐1); r[ij] = ‐vrtx[ij]//右辺が‐ωなのでマイナス付き ‐( (stmf[im1j]‐2*stmf[ij]+stmf[ip1j])/dxdx +(stmf[ijm1]‐2*stmf[ij]+stmf[ijp1])/dydy); rr += r[ij]*r[ij]; bb += vrtx[ij]*vrtx[ij]; }

流れ関数の計算（共役勾配法）

taylor_green.c

do{ for(int j=0;j<Ny;j++) for(int i=0;i<Nx;i++){ ij = i+Nx*j; p [ij] = r[ij] + beta*p[ij]; } pAp = 0.0; for(int j=1;j<Ny‐1;j++) for(int i=1;i<Nx‐1;i++){ ij = i+Nx*j; ip1j = (i+1)+Nx*(j  ); im1j = (i‐1)+Nx*(j  ); ijp1 = (i )+Nx*(j+1); ijm1 = (i )+Nx*(j‐1); Ap[ij] = (p[im1j]‐2*p[ij]+p[ip1j])/dxdx +(p[ijm1]‐2*p[ij]+p[ijp1])/dydy; pAp += p[ij]*Ap[ij]; } alph = rr/pAp; rr = 0.0; for(int j=0;j<Ny;j++) for(int i=0;i<Nx;i++){ ij = i+Nx*j; stmf[ij] = stmf[ij] + alph* p[ij]; r   [ij] = r   [ij] ‐ alph*Ap[ij]; rr += r[ij]*r[ij]; } err = sqrt(rr/bb); beta = rr/(alph*pAp); }while(err > ERR_TOL); free(r); free(p); free(Ap); }

流れ関数の計算（共役勾配法）

taylor_green.c

void computeVelocity (double *stmf, double *u, double *v){ int i,j; int ij,im1j,ip1j,ijm1,ijp1; int im2j,ip2j,ijm2,ijp2; //計算領域内部 for(j=1;j<Ny‐1;j++) for(i=1;i<Nx‐1;i++){ ij = i +Nx*j; im1j = i‐1+Nx*j; ip1j = i+1+Nx*j; ijm1 = i +Nx*(j‐1); ijp1 = i +Nx*(j+1); u[ij] = (stmf[ijp1]‐stmf[ijm1])/dy2; v[ij] =‐(stmf[ip1j]‐stmf[im1j])/dx2; } //下境界 j=0; for(i=1;i<Nx‐1;i++){ ij = i +Nx*j; im1j = i‐1+Nx*j; ip1j = i+1+Nx*j; ijp1 = i +Nx*(j+1); ijp2 = i +Nx*(j+2); u[ij] = (‐3.0*stmf[ij]+4.0*stmf[ijp1] ‐stmf[ijp2])/dy2; v[ij] =‐(stmf[ip1j]‐stmf[im1j])/dx2; } //上境界 j=Ny‐1; for(i=1;i<Nx‐1;i++){ ij = i +Nx*j; im1j = i‐1+Nx*j; ip1j = i+1+Nx*j; ijm1 = i +Nx*(j‐1); ijm2 = i +Nx*(j‐2); u[ij] = ( 3.0*stmf[ij]‐4.0*stmf[ijm1] +stmf[ijm2])/dy2; v[ij] =‐(stmf[ip1j]‐stmf[im1j])/dx2; }

プログラム（速度の計算）

taylor_green.c

for(j=1;j<Ny‐1;j++){ //左境界 i=0; ij = i +Nx*j; im1j = i‐1+Nx*j; ip1j = i+1+Nx*j; ip2j = i+2+Nx*j; ijp1 = i +Nx*(j+1); ijm1 = i +Nx*(j‐1); u[ij] = (stmf[ijp1]‐stmf[ijm1])/dy2; v[ij] =‐(‐3.0*stmf[ij]+4.0*stmf[ip1j] ‐stmf[ip2j])/dx2; //右境界 i=Nx‐1; ij = i +Nx*j; im2j = i‐2+Nx*j; im1j = i‐1+Nx*j; ip1j = i+1+Nx*j; ijp1 = i +Nx*(j+1); ijm1 = i +Nx*(j‐1); u[ij] = (stmf[ijp1]‐stmf[ijm1])/dy2; v[ij] =‐( 3.0*stmf[ij]‐4.0*stmf[im1j] +stmf[im2j])/dx2; } //左下角 i=0   ;j=0; ij = i +Nx*j; ip1j = i+1+Nx*j; ip2j = i+2+Nx*j; ijp1 = i +Nx*(j+1); ijp2 = i +Nx*(j+2); u[ij] = (‐3.0*stmf[ij]+4.0*stmf[ijp1] ‐stmf[ijp2])/dy2; v[ij] =‐(‐3.0*stmf[ij]+4.0*stmf[ip1j] ‐stmf[ip2j])/dx2;

プログラム（速度の計算）

taylor_green.c

//右下角 i=Nx‐1;j=0; ij = i +Nx*j; im2j = i‐2+Nx*j; im1j = i‐1+Nx*j; ijp1 = i +Nx*(j+1); ijp2 = i +Nx*(j+2); u[ij] = (‐3.0*stmf[ij]+4.0*stmf[ijp1] ‐stmf[ijp2])/dy2; v[ij] =‐( 3.0*stmf[ij]‐4.0*stmf[im1j] +stmf[im2j])/dx2; //左上角 i=0   ;j=Ny‐1; ij = i +Nx*j; ip1j = i+1+Nx*j; ip2j = i+2+Nx*j; ijm1 = i +Nx*(j‐1); ijm2 = i +Nx*(j‐2); u[ij] = ( 3.0*stmf[ij]‐4.0*stmf[ijm1] +stmf[ijm2])/dy2; v[ij] =‐(‐3.0*stmf[ij]+4.0*stmf[ip1j] ‐stmf[ip2j])/dx2; //右上角 i=Nx‐1;j=Ny‐1; ij = i +Nx*j; im2j = i‐2+Nx*j; im1j = i‐1+Nx*j; ijm1 = i +Nx*(j‐1); ijm2 = i +Nx*(j‐2); u[ij] = ( 3.0*stmf[ij]‐4.0*stmf[ijm1] +stmf[ijm2])/dy2; v[ij] =‐( 3.0*stmf[ij]‐4.0*stmf[im1j] +stmf[im2j])/dx2; }

プログラム（速度の計算）

taylor_green.c

速度の計算



計算する位置に応じて差分に使う点が変化



計算領域内は

x, y方向とも中心差分

for(j=1;j<Ny‐1;j++) for(i=1;i<Nx‐1;i++){ ij = i +Nx*j; im1j = i‐1+Nx*j; ip1j = i+1+Nx*j; ijm1 = i +Nx*(j‐1); ijp1 = i +Nx*(j+1); u[ij] = (stmf[ijp1]‐stmf[ijm1])/dy2; v[ij] =‐(stmf[ip1j]‐stmf[im1j])/dx2; }

速度の計算



計算する位置に応じて差分に使う点が変化



下方向の境界



x方向は中心差分



y方向は片側差分（j+1,j+2を参照）

j=0; for(i=1;i<Nx‐1;i++){ ij = i +Nx*j; im1j = i‐1+Nx*j; ip1j = i+1+Nx*j; ijp1 = i +Nx*(j+1); ijp2 = i +Nx*(j+2); u[ij] = (‐3.0*stmf[ij]+4.0*stmf[ijp1] ‐stmf[ijp2])/dy2; v[ij] =‐(stmf[ip1j]‐stmf[im1j])/dx2; }

速度の計算



計算する位置に応じて差分に使う点が変化



上方向の境界



x方向は中心差分



y方向は片側差分（j‐1,j‐2を参照）

j=Ny‐1; for(i=1;i<Nx‐1;i++){ ij = i +Nx*j; im1j = i‐1+Nx*j; ip1j = i+1+Nx*j; ijm1 = i +Nx*(j‐1); ijm2 = i +Nx*(j‐2); u[ij] = ( 3.0*stmf[ij]‐4.0*stmf[ijm1] +stmf[ijm2])/dy2; v[ij] =‐(stmf[ip1j]‐stmf[im1j])/dx2; }

速度の計算



計算する位置に応じて差分に使う点が変化



左方向の境界



x方向は片側差分（i+1,i+2を参照）



y方向は中心差分

for(j=1;j<Ny‐1;j++){ i=0; ij = i +Nx*j; im1j = i‐1+Nx*j; ip1j = i+1+Nx*j; ip2j = i+2+Nx*j; ijp1 = i +Nx*(j+1); ijm1 = i +Nx*(j‐1); u[ij] = (stmf[ijp1]‐stmf[ijm1])/dy2; v[ij] =‐(‐3.0*stmf[ij]+4.0*stmf[ip1j] ‐stmf[ip2j])/dx2; }

速度の計算



計算する位置に応じて差分に使う点が変化



右方向の境界



x方向は片側差分（i‐1,i‐2を参照）



y方向は中心差分

for(j=1;j<Ny‐1;j++){ i=Nx‐1; ij = i +Nx*j; im2j = i‐2+Nx*j; im1j = i‐1+Nx*j; ip1j = i+1+Nx*j; ijp1 = i +Nx*(j+1); ijm1 = i +Nx*(j‐1); u[ij] = (stmf[ijp1]‐stmf[ijm1])/dy2; v[ij] =‐( 3.0*stmf[ij]‐4.0*stmf[im1j] +stmf[im2j])/dx2; }

速度の計算



計算する位置に応じて差分に使う点が変化



左下の境界



x方向は片側差分（i+1,i+2を参照）



y方向は片側差分（j+1,j+2を参照）

i=0   ;j=0; ij = i +Nx*j; ip1j = i+1+Nx*j; ip2j = i+2+Nx*j; ijp1 = i +Nx*(j+1); ijp2 = i +Nx*(j+2); u[ij] = (‐3.0*stmf[ij]+4.0*stmf[ijp1] ‐stmf[ijp2])/dy2; v[ij] =‐(‐3.0*stmf[ij]+4.0*stmf[ip1j] ‐stmf[ip2j])/dx2;

速度の計算



計算する位置に応じて差分に使う点が変化



右下の境界



x方向は片側差分（i‐1,i‐2を参照）



y方向は片側差分（j+1,j+2を参照）

i=Nx‐1;j=0; ij = i +Nx*j; im2j = i‐2+Nx*j; im1j = i‐1+Nx*j; ijp1 = i +Nx*(j+1); ijp2 = i +Nx*(j+2); u[ij] = (‐3.0*stmf[ij]+4.0*stmf[ijp1] ‐stmf[ijp2])/dy2; v[ij] =‐( 3.0*stmf[ij]‐4.0*stmf[im1j] +stmf[im2j])/dx2;

速度の計算



計算する位置に応じて差分に使う点が変化



左上の境界



x方向は片側差分（i+1,i+2を参照）



y方向は片側差分（j‐1,j‐2を参照）

i=0   ;j=Ny‐1; ij = i +Nx*j; ip1j = i+1+Nx*j; ip2j = i+2+Nx*j; ijm1 = i +Nx*(j‐1); ijm2 = i +Nx*(j‐2); u[ij] = ( 3.0*stmf[ij]‐4.0*stmf[ijm1] +stmf[ijm2])/dy2; v[ij] =‐(‐3.0*stmf[ij]+4.0*stmf[ip1j] ‐stmf[ip2j])/dx2;

速度の計算



計算する位置に応じて差分に使う点が変化



右上の境界



x方向は片側差分（i‐1,i‐2を参照）



y方向は片側差分（j‐1,j‐2を参照）

i=Nx‐1;j=Ny‐1; ij = i +Nx*j; im2j = i‐2+Nx*j; im1j = i‐1+Nx*j; ijm1 = i +Nx*(j‐1); ijm2 = i +Nx*(j‐2); u[ij] = ( 3.0*stmf[ij]‐4.0*stmf[ijm1] +stmf[ijm2])/dy2; v[ij] =‐( 3.0*stmf[ij]‐4.0*stmf[im1j] +stmf[im2j])/dx2;

計算結果

計算条件



_{2次元正方キャビティ流れ}



正方形断面のくぼみに水が満

たされている



くぼみにフタが置かれ，フタが

一定速度

_Uで移動



計算条件



レイノルズ数：

Re=UL/=1000



計算時間間隔：

U



t/



x=0.1



計算時間：

Ut/L=200まで

U L

境界条件



流れ関数



流れ関数の定義式から計算



全ての壁上で



=const.（一定値として0を採用）

0    



v 0      v x



0     u y



0     u y



より_=const. より_=const. より_=const. より=const.

境界条件



渦度



流れ関数の定義式と流れ関数のPoisson方程式から計算





















2 2 2 2

y

x

2 2

x













, 2 2 1 , 1 1 , 1 , 1 1 , Δx n j i n j i n j i n j i       _{ }













0 2 1 , 1 1 , 1 1 ,          Δx v n j i n j i n j i





1 , 1 1 , 1      in j n j i





側壁面の境界条件

壁面上で=0より 離散化されたPoisson方程式と流れ関数の定義式を連立すると 0 Right Wall on Left Wall on 2 2 2 1 , 1 2 1 , 1 1 , Δx Δx n j i n j i n j i         







1 ,  n j i



1 , 1   n j i



1 , 1   n j i



1 ,  n j i



Δx

境界条件



渦度



流れ関数の定義式と流れ関数のPoisson方程式から計算





















2 2 2 2

y

x

2 2

y













, 2 2 1 1 , 1 , 1 1 , 1 , Δy n j i n j i n j i n j i       _{ }













0 2 1 1 , 1 1 , 1 ,         Δy u n j i n j i n j i





1 1 , 1 1 ,      inj n j i





下壁面の境界条件

壁面上で=0より 離散化されたPoisson方程式と流れ関数の定義式を連立すると 0 Wall Bottom on 2 1 2 1 , 1 , Δy n j i n j i    _{ }

_



1 ,  n j i



1 1 ,   n j i



Δy

境界条件



渦度



流れ関数の定義式と流れ関数のPoisson方程式から計算





















2 2 2 2

y

x

2 2

y













, 2 2 1 1 , 1 , 1 1 , 1 , Δy n j i n j i n j i n j i       _{ }













U Δy u n j i n j i n j i         2 1 1 , 1 1 , 1 ,





1 1 , 1 1 , 2       inj n j i ΔyU





上壁面の境界条件

壁面上で=0より 離散化されたPoisson方程式と流れ関数の定義式を連立すると 0





on Top Wall 2 1 2 1 , 1 , ΔyU Δy n j i n j i      

_



1 ,  n j i



1 1 ,   n j i



Δy U

境界条件

2 1 1 , 1 , 1 , 2 0 Δy n j i n j i n j i       







2 1 , 1 1 , 1 , 2 0 Δx n j i n j i n j i       







2 1 , 1 1 , 1 , 2 0 Δx n j i n j i n j i       









1



2 1 , 1 , 1 , 2 0 Δy ΔyU n j i n j i n j i        







#include<stdio.h> #include<stdlib.h> #include<math.h> #define Re 1000 //レイノルズ数 #define Lx  0.01  //溝の幅（=高さ） #define Ly  Lx #define Uwall 0.1 //移動壁の速度 #define Kvisc (Lx*Uwall/Re) #define Nx  81 #define Ny  Nx #define Nbytes (Nx*Ny*sizeof(double)) #define dx  (Lx/(Nx‐1)) #define dy  (Ly/(Ny‐1)) #define dx2  (dx*2.0) #define dy2  (dy*2.0) #define dxdx (dx*dx) #define dydy (dy*dy) #define dxdxdydy (dxdx*dydy) #define dxdy2 (2.0*(dxdx+dydy)) #define ndt (0.1)         //CFL条件 #define dt (ndt*dx/Uwall) //時間刻み #define endnT (100.0)     //終了時間（無次元） #define Nt (int)(endnT/dt*Lx/Uwall) #define ERR_TOL 1e‐12 #define Accel 1.7

パラメータ設定

cavity.c

void init(double *vrtx, double *stmf, double *u, double *v, double *vrtx_next){ int i,j; for(j=0;j<Ny;j++) for(i=0;i<Nx;i++){ double x = (double)i*dx; double y = (double)j*dy; vrtx[i+Nx*j] = 0.0; stmf[i+Nx*j] = 0.0; v[i+Nx*j] = 0.0; u[i+Nx*j] = 0.0; vrtx_next[i+Nx*j] = 0.0; } //速度の境界条件．上壁面に速度Uwallを与える j=Ny‐1; for(i=0;i<Nx;i++){ u[i+Nx*j] = Uwall; } }

初期化

cavity.c

void computeVorticity (double *vrtx, double *stmf, double *u, double *v,double *vrtx_next){ int i,j; int ij,im1j,ip1j,ijm1,ijp1; double conv,visc; for(j=1;j<Ny‐1;j++){ //左壁面 i=0; ij = i +Nx*j; ip1j = i+1+Nx*j; vrtx[ij] = ‐2.0*stmf[ip1j]/dxdx; vrtx_next[ij] = vrtx[ij]; //右壁面 i=Nx‐1; ij = i +Nx*j; im1j = i‐1+Nx*j; vrtx[ij] = ‐2.0*stmf[im1j]/dxdx; vrtx_next[ij] = vrtx[ij]; }

渦度方程式の計算（時間積分）

2 1 , 1 1 , 2 Δx n j i n j i    _{ }

_



2 1 , 1 1 , 2 Δx n j i n j i    _{ }

_



cavity.c

//下壁面 j=0; for(i=0;i<Nx;i++){ ij = i+Nx* j; ijp1 = i+Nx*(j+1); vrtx[ij] = ‐2.0*stmf[ijp1]/dydy; vrtx_next[ij] = vrtx[ij]; } //上壁面（移動壁） j=Ny‐1; for(i=0;i<Nx;i++){ ij = i+Nx* j; ijm1 = i+Nx*(j‐1);

vrtx[ij] = ‐2.0*(stmf[ijm1]+Uwall*dy)/dydy; vrtx_next[ij] = vrtx[ij]; }

渦度方程式の計算（時間積分）

2 1 1 , 1 , 2 Δy n j i n j i    _{ }

_





1



2 1 , 1 , 2 ΔyU Δy n j i n j i      

_



U

cavity.c

//計算領域内部 for(j=1;j<Ny‐1;j++) for(i=1;i<Nx‐1;i++){ ij = i +Nx*j; im1j = i‐1+Nx*j; ip1j = i+1+Nx*j; ijm1 = i +Nx*(j‐1); ijp1 = i +Nx*(j+1); //移流項を計算 conv = u[ij]*(vrtx[ip1j]‐vrtx[im1j])/dx2 + v[ij]*(vrtx[ijp1]‐vrtx[ijm1])/dy2; //粘性項を計算 visc = Kvisc*( (vrtx[im1j]‐2.0*vrtx[ij]+vrtx[ip1j])/dxdx +(vrtx[ijm1]‐2.0*vrtx[ij]+vrtx[ijp1])/dydy); //時間積分 vrtx_next[ij] = vrtx[ij] + dt*(‐conv+visc); } }

渦度方程式の計算（時間積分）

cavity.c

void computeVelocity(double *stmf, double *u, double *v){ int i,j; int ij,im1j,ip1j,ijm1,ijp1; int im2j,ip2j,ijm2,ijp2; //境界条件は初期設定の際に設定済み //壁の移動速度は時間で変化しないので何もしなくてよい for(j=1;j<Ny‐1;j++) for(i=1;i<Nx‐1;i++){ ij = i +Nx*j; im1j = i‐1+Nx*j; ip1j = i+1+Nx*j; ijm1 = i +Nx*(j‐1); ijp1 = i +Nx*(j+1); u[ij] = (stmf[ijp1]‐stmf[ijm1])/dy2; v[ij] =‐(stmf[ip1j]‐stmf[im1j])/dx2; } }

プログラム（速度の計算）

cavity.c

計算結果