/********************************************/ /* 二次元キャビテイー流れ解析 */ /* */ /* */ /* */ /* */ /* */ /* ver. 1.0 */ /* */ /********************************************/ #include #include #include #include #define ID 102 /* ｃは０から配列を取る */ #define JD 102 /****************/ /** 変数宣言 **/ /****************/ double x[ID][JD],y[ID][JD],pr_x[ID][JD],pr_y[ID][JD]; double u[ID][JD],v[ID][JD],p[ID][JD]; double xx[ID][JD],yx[ID][JD],xe[ID][JD],ye[ID][JD],yaci[ID][JD]; double c1[ID][JD],c2[ID][JD],c3[ID][JD],c4[ID][JD],c5[ID][JD]; double uo[ID][JD],vo[ID][JD],un[ID][JD],vn[ID][JD]; int nx,ny,nxm1,nym1,nxp1,nyp1; int iter1,iter2,nstep,iprog,icon; double eps1,omega1,eps2,omega2; double re,rei,dt,time; char fn[25],res[15],tim[15],text[30],string[10]; FILE *fp1,*fp2,*fopen(); /**********************/ /** メインルーチン **/ /**********************/ main() { int istep; init (); f_read (); v_trans (); time=0.0 ; for(istep=1 ; istep <= nstep ; istep++){ time = time+dt ; printf( " \n"); printf( " 時間 = %8.4f \n",time); printf( " \n"); boundp (); poiss (); boundv (); navi () ; } f_write (); } /********************************/ /** ファイルデータの読み込み **/ /********************************/ int f_read () { int i,j; double xx,yy; printf (" 座標ファイル名を入力して下さい \n"); scanf ("%s",fn); puts (fn); printf("座標デ−タを読み込み中 \n"); fp1 = fopen(fn,"r"); fscanf(fp1,"%d %d",&nx,&ny); printf( "ｘ方向格子点= %d ｙ方向格子点= %d \n",nx,ny); nxm1 = nx-1 ; nym1 = ny-1 ; nxp1 = nx+1 ; nyp1 = ny+1 ; for (j=1;j<=ny;++j) { for (i=1;i<=nx;++i) { fscanf(fp1,"%lf %lf",&xx,&yy); x[i][j]=xx; y[i][j]=yy; } } return (0); } /********************************/ /** ファイルデータの書き込み **/ /********************************/ int f_write () { int i,j; printf (" 出力ファイル名を入力して下さい \n"); scanf ("%s",fn); fp1 = fopen(fn,"w"); fprintf(fp1,"%10d %10d\n",nx,ny); for (j=1;j<=ny;++j) { for (i=1;i<=nx;++i) { fprintf(fp1,"%10.7f %10.7f %15.7e\n",u[i][j],v[i][j],p[i][j]); } } } /***************************** ** subroutine initial data ** *****************************/ init () { /*** iteration parameter ***/ iter1=150; iter2= 30; /*** relaxation parameter ***/ omega1=0.9; omega2=1.0; /*** reynolds number ***/ re =3200.0; rei=1.0/re; /*** convergence parameter ***/ eps1=1.0e-3; eps2=1.0e-7; /*** pre-result reading yes=1 no=0 **/ icon=0; /*** step for history ***/ iprog=20000; /*** time step ***/ dt=0.002; nstep=75000; } /*******************************/ /*** compute stream function ***/ /*******************************/ v_trans () { int i,j; double yac,xxx,xxe,xee,yxx,yxe,yee,dx,dy ; double alpha,beta,gamma,yacs,yacc; /*** compute metric ***/ for(j=2 ; j<=ny-1 ; j++) { for(i=2 ; i<=nx-1 ; i++) { xx[i][j]=0.5*(x[i+1][j]-x[i-1][j]); yx[i][j]=0.5*(y[i+1][j]-y[i-1][j]); xe[i][j]=0.5*(x[i][j+1]-x[i][j-1]); ye[i][j]=0.5*(y[i][j+1]-y[i][j-1]); yac=xx[i][j]*ye[i][j]-xe[i][j]*yx[i][j]; yaci[i][j]=1.0/yac ; } } for(j=2 ; j<=ny-1 ; j++) { for(i=2 ; i<=nx-1 ; i++) { yacs=yaci[i][j]*yaci[i][j]; yacc=yacs*yaci[i][j]; xxx=x[i+1][j]-2.0*x[i][j]+x[i-1][j]; xee=x[i][j+1]-2.0*x[i][j]+x[i][j-1]; yxx=y[i+1][j]-2.0*y[i][j]+y[i-1][j]; yee=y[i][j+1]-2.0*y[i][j]+y[i][j-1]; xxe=0.25*(x[i+1][j+1]-x[i-1][j+1]-x[i+1][j-1]+x[i-1][j-1]); yxe=0.25*(y[i+1][j+1]-y[i-1][j+1]-y[i+1][j-1]+y[i-1][j-1]); alpha=xe[i][j]*xe[i][j]+ye[i][j]*ye[i][j]; beta =xx[i][j]*xe[i][j]+yx[i][j]*ye[i][j]; gamma=xx[i][j]*xx[i][j]+yx[i][j]*yx[i][j]; dx=alpha*xxx-2.0*beta*xxe+gamma*xee; dy=alpha*yxx-2.0*beta*yxe+gamma*yee; c1[i][j] = alpha*yacs; c2[i][j] =-2.0*beta*yacs; c3[i][j] = gamma*yacs; c4[i][j] = (xe[i][j]*dy-ye[i][j]*dx)*yacc; c5[i][j] = (yx[i][j]*dx-xx[i][j]*dy)*yacc; } } } /************************************* ** substituting boundary condition ** *************************************/ boundv () { int i,j; for(i=2 ; i<=nxm1 ; i++) { u[i][nyp1]= u[i][nym1]; v[i][nyp1]=-v[i][nym1]; u[i][ny] =1.0; v[i][ny] =0.0; } for(i=1 ; i<=nx ; i++) { u[i][1] = 0.0; v[i][1] = 0.0; u[i][0] =-u[i][2]; v[i][0] =-v[i][2]; } for(j=1 ; j<=ny ; j++) { u[1][j] = 0.0; v[1][j] = 0.0; u[0][j] =-u[2][j]; v[0][j] =-v[2][j]; u[nx ][j] = 0.0; v[nx ][j] = 0.0; u[nxp1][j] =-u[nxm1][j]; v[nxp1][j] =-v[nxm1][j]; } } /************************************* ** substituting boundary condition ** *************************************/ boundp () { int i,j; for(i=1 ; i<=nx ; i++) { p[i][ny] = (4.0*p[i][nym1]-p[i][ny-2])/3.0 ; } /*** on the body surface ***/ for(i=1 ; i<=nx ; i++) { p[i][1] = (4.0*p[i][2]-p[i][3])/3.0 ; } for(j=1 ; j<=ny ; j++){ p[1 ][j] = (4.0*p[2 ][j]-p[3 ][j])/3.0 ; p[nx][j] = (4.0*p[nxm1][j]-p[nx-2][j])/3.0 ; } } /********************************** ** poisson equation of pressure ** **********************************/ poiss () { int i,j,ip1,jp1,im1,jm1,loop; double dudx,dudy,dvdx,dvdy,po,dp,dcon,ddd,div; double rp,err; double ux,ue,vx,ve,rhs[ID][JD],dmyp[ID]; /***** right hand side of poisson equation *****/ for( j=2 ; j<=nym1 ; j++){ jp1=j+1; jm1=j-1; for( i=2 ; i<=nxm1 ; i++){ ux=.5*(u[i+1][j]-u[i-1][j]); ue=.5*(u[i][jp1]-u[i][jm1]); vx=.5*(v[i+1][j]-v[i-1][j]); ve=.5*(v[i][jp1]-v[i][jm1]); dudx=ux*ye[i][j]-ue*yx[i][j]; dvdy=ve*xx[i][j]-vx*xe[i][j]; dudy=ue*xx[i][j]-ux*xe[i][j]; dvdx=vx*ye[i][j]-ve*yx[i][j]; ddd=dudx*dudx+dvdy*dvdy+2.0*dudy*dvdx; dcon=ddd*yaci[i][j]*yaci[i][j]; /*** cal. divergence ***/ div=(dudx+dvdy)*yaci[i][j]; rhs[i][j]=dcon-div/dt; } } /***** slor iteration loop *****/ for( loop=1 ; loop<=iter1 ; loop++){ err=0.0; /**** cal. internal points ****/ for( j=2 ; j<=nym1 ; j++){ jp1=j+1; jm1=j-1; for( i=2 ; i<=nxm1 ; i++){ ip1=i+1; im1=i-1; po=p[i][j]; /***** finite difference equation for poisson eq. *****/ dp=c1[i][j]*(p[ip1][j]+p[im1][j]) +c2[i][j]*(p[ip1][jp1]-p[im1][jp1]-p[ip1][jm1]+p[im1][jm1])*.25 +c3[i][j]*(p[i][jp1]+p[i][jm1]) +c4[i][j]*(p[ip1][j]-p[im1][j])*.5 +c5[i][j]*(p[i][jp1]-p[i][jm1])*.5+rhs[i][j]; /***** residual *****/ rp=.5*dp/(c1[i][j]+c3[i][j])-po; dmyp[i]=po+omega1*rp; err=err+rp*rp; } for( i=2 ; i<=nxm1 ; i++){ p[i][j]=dmyp[i]; } } boundp () ; err=sqrt(err/(float)((nxm1-1)*(nym1-1))); if(err <= eps1) goto pass_sorp; } pass_sorp:; printf( "圧力反復回数 = %5d 残差 = %10.5e \n",loop,err); } /*********************************** ** solved navier-stokes equation ** ***********************************/ navi () { int i,j,im1,jm1,ip1,jp1,im2,jm2,ip2,jp2,loop; double conv,vis,cx,cy,cpx,cpy,vx,vy,du,dv,drx,dry,d12,err; double u1,u2,u3,u4,v1,v2,v3,v4; double coef[ID][JD],px[ID][JD],py[ID][JD],cu[ID][JD],cv[ID][JD], dmyu[ID],dmyv[ID]; for( j=2 ; j<=nym1 ; j++){ for( i=2 ; i<=nxm1 ; i++){ un[i][j]=uo[i][j]; vn[i][j]=vo[i][j]; uo[i][j]= u[i][j]; vo[i][j]= v[i][j]; /***** contravariant velocity component *****/ cu[i][j]=(uo[i][j]*ye[i][j]-vo[i][j]*xe[i][j])*yaci[i][j]; cv[i][j]=(vo[i][j]*xx[i][j]-uo[i][j]*yx[i][j])*yaci[i][j]; conv=0.5*(fabs(cu[i][j])+fabs(cv[i][j])); vis =2.0*rei*(c1[i][j]+c3[i][j]); coef[i][j]=1.0/(3.0+2.0*dt*(conv+vis)); /***** pressure term *****/ px[i][j]=( ye[i][j]*(p[i+1][j]-p[i-1][j]) -yx[i][j]*(p[i][j+1]-p[i][j-1]))*.5*yaci[i][j]; py[i][j]=(-xe[i][j]*(p[i+1][j]-p[i-1][j]) +xx[i][j]*(p[i][j+1]-p[i][j-1]))*.5*yaci[i][j]; } } d12=1.0/12.0; /*** slor iteration loop ***/ for(loop=1 ; loop <=iter2 ; loop ++ ){ err=0.0; for( j=2 ; j<=nym1 ; j++){ jp1=j+1; jp2=j+2; jm1=j-1; jm2=j-2; for( i=2 ; i<=nxm1 ; i++){ ip1=i+1; ip2=i+2; im1=i-1; im2=i-2; /***** convective term (kawamura scheme) *****/ u1=-u[ip2][j]+8.0*(u[ip1][j]-u[im1][j])+u[im2][j]; u2= u[ip2][j]-4.0*(u[ip1][j]+u[im1][j])+u[im2][j]; u3=-u[i][jp2]+8.0*(u[i][jp1]-u[i][jm1])+u[i][jm2]; u4= u[i][jp2]-4.0*(u[i][jp1]+u[i][jm1])+u[i][jm2]; v1=-v[ip2][j]+8.0*(v[ip1][j]-v[im1][j])+v[im2][j]; v2= v[ip2][j]-4.0*(v[ip1][j]+v[im1][j])+v[im2][j]; v3=-v[i][jp2]+8.0*(v[i][jp1]-v[i][jm1])+v[i][jm2]; v4= v[i][jp2]-4.0*(v[i][jp1]+v[i][jm1])+v[i][jm2]; cx=cu[i][j]*u1+fabs(cu[i][j])*u2+cv[i][j]*u3+fabs(cv[i][j])*u4; cy=cu[i][j]*v1+fabs(cu[i][j])*v2+cv[i][j]*v3+fabs(cv[i][j])*v4; cpx=-cx*d12-px[i][j]; cpy=-cy*d12-py[i][j]; /***** viscous term *****/ vx=c1[i][j]*(u[ip1][j]+u[im1][j]) +c2[i][j]*(u[ip1][jp1]-u[im1][jp1]-u[ip1][jm1]+u[im1][jm1])*.25 +c3[i][j]*(u[i][jp1]+u[i][jm1]) +c4[i][j]*(u[ip1][j]-u[im1][j])*.5 +c5[i][j]*(u[i][jp1]-u[i][jm1])*.5; vy=c1[i][j]*(v[ip1][j]+v[im1][j]) +c2[i][j]*(v[ip1][jp1]-v[im1][jp1]-v[ip1][jm1]+v[im1][jm1])*.25 +c3[i][j]*(v[i][jp1]+v[i][jm1]) +c4[i][j]*(v[ip1][j]-v[im1][j])*.5 +c5[i][j]*(v[i][jp1]-v[i][jm1])*.5; drx=cpx+vx*rei; dry=cpy+vy*rei; du=coef[i][j]*(4.0*uo[i][j]-un[i][j]+2.0*dt*drx)-u[i][j]; dv=coef[i][j]*(4.0*vo[i][j]-vn[i][j]+2.0*dt*dry)-v[i][j]; err=err+du*du+dv*dv; dmyu[i]=u[i][j]+omega2*du; dmyv[i]=v[i][j]+omega2*dv; } for( i=2 ; i<=nxm1 ; i++){ u[i][j]=dmyu[i]; v[i][j]=dmyv[i]; } } boundv (); err=sqrt(err/(float)((nxm1-1)*(nym1-1))); if(err <= eps2) goto pass_soruv; } pass_soruv:; printf( "流速反復回数 = %5d 残差 = %10.5e \n",loop,err); }