圆柱绕流程序，结果为tecplot文件_.c程序结果用tecplot, fortran输出二进制数据到tecplot-Java文

#include"common.h" void main() { double tet,vrx,vtt,vrp,vrm,vtp,vtm,vrc,vtc,vrt,vtx,err,rhs,temp; double velo[2],exx,vrl,vrn,vtl,vtn; fstream id1; id1.open("end.dat",ios::out); na=121,nb=40,nx=na+1,ny=nb+1,re=100,dt=0.008,dy=0.08,nmax=10000,kk=40,coef=1; dx=2*pi/double(nx-2); dxi=1/dx; dyi=1/dy; rei=1/re; dti=1/dt; dx2=dxi*dxi; dy2=dyi*dyi; fct=1/(2*dx2+2*dy2); for(j=0;j<ny;j++) { for(i=0;i<nx;i++) { te1=dx*i; radi=exp(dy*j); vr[i][j]=0; vt[i][j]=0; p[i][j]=0; } } for(j=0;j<ny;j++) { for(i=0;i<nx;i++) { tet=dx*i; x[i][j]=exp(dy*(j-1))*cos(tet); y[i][j]=exp(dy*(j-1))*sin(tet); } } for(l=1;l<nmax;l++) { for(i=0;i<nx;i++) { vr[i][0]=0; vr[i][1]=0; vt[i][0]=0; } for(i=0;i<nx;i++) { te1=dx*i; vr[i][ny-1]=cos(te1); vt[i][ny-1]=-sin(te1); vr[i][ny-2]=cos(te1); } for(j=0;j<ny;j++) { vr[0][j]=vr[nx-2][j]; vr[nx-1][j]=vr[1][j]; vt[0][j]=vt[nx-2][j]; vt[nx-1][j]=vt[1][j]; vt[nx][j]=vt[2][j]; } for(j=1;j<ny-1;j++) { for(i=1;i<nx-1;i++) { vrx=(vr[i][j+1]-vr[i][j])*dyi; vtt=(vt[i+1][j]-vt[i][j])*dxi; vrp=0.25*(vr[i+1][j+1]+vr[i+1][j]+vr[i][j+1]+vr[i][j]); vrm=0.25*(vr[i][j+1]+vr[i][j]+vr[i-1][j+1]+vr[i-1][j]); vtp=0.25*(vt[i+1][j+1]+vt[i+1][j]+vt[i][j+1]+vt[i][j]); vtm=0.25*(vt[i+1][j]+vt[i+1][j-1]+vt[i][j]+vt[i][j-1]); vrc=0.5*(vr[i][j+1]+vr[i][j]); vtc=0.5*(vt[i+1][j]+vt[i][j]); vrt=(vrp-vrm)*dxi; vtx=(vtp-vtm)*dyi; q[i][j]=-vrx*vrx-2*vrt*vtx-vtt*vtt-vrc*vrc-2*(vrc*vtt-vtc*vtx)+exp(dy*(j-1))*(vrx+vrc+vtt)*dti; } } for(k=0;k<kk;k++) { err=0; for(i=1;i<nx-1;i++) { p[i][0]=p[i][1]; p[i][ny-1]=p[i][ny-2]; } for(j=0;j<ny;j++) { p[0][j]=p[nx-2][j]; p[nx-1][j]=p[1][j]; } for(j=1;j<ny-1;j++) { for(i=1;i<nx-1;i++) { rhs=((p[i+1][j]+p[i-1][j])*dx2+(p[i][j+1]+p[i][j-1])*dy2-q[i][j])*fct; err=err+pow(rhs-p[i][j],2); p[i][j]=p[i][j]*(1-coef)+rhs*coef; } } m=k; if(err<0.00001) { k=kk; } } temp=large; for(i=0;i<nx;i++) { for(j=1;j<ny;j++) { if(temp>p[i][j]) { temp=p[i][j]; } } } for(i=0;i<nx;i++) { for(j=1;j<ny;j++) { p[i][j]=p[i][j]-temp; } } if((l+1)%1==0) { cout.precision(5); cout.setf(ios::showpoint); double time=l*dt; cout<<"当前时间:"<<time<<setw(10)<<"内迭代:"<<m<<setw(10)<<"误差:"<<err<<endl; } for(j=1;j<ny-1;j++) { for(i=1;i<nx-1;i++) { vtc=0.25*(vt[i+1][j]+vt[i][j]+vt[i+1][j-1]+vt[i][j-1]); exx=exp(-j*dy); vrl=(vr[i+1][j]-2*vr[i][j]+vr[i-1][j])*dx2+(vr[i][j+1]-2*vr[i][j]+vr[i][j-1])*dy2-vr[i][j]- (vt[i+1][j]+vt[i+1][j-1]-vt[i][j]-vt[i][j-1])*dxi; vrn=vr[i][j]*(vr[i][j+1]-vr[i][j-1])*dyi*0.5+vtc*(vr[i+1][j]-vr[i-1][j])*dxi*0.5-vtc*vtc; temp1[i][j]=vr[i][j]+dt*(-exx*(vrn+(p[i][j]-p[i][j-1])*dyi)+exx*exx*rei*vrl); } } for(j=1;j<ny-1;j++) { for(i=1;i<nx-1;i++) { vrc=0.25*(vr[i][j+1]+vr[i][j]+vr[i-1][j+1]+vr[i-1][j]); exx=exp(-j*dy); vtl=(vt[i+1][j]-2*vt[i][j]+vt[i-1][j])*dx2+(vt[i][j+1]-2*vt[i][j]+vt[i][j-1])*dy2-vt[i][j]+ (vr[i][j+1]+vr[i][j]-vr[i-1][j+1]-vr[i-1][j])*dxi; vtn=vrc*(vt[i][j+1]-vt[i][j-1])*dyi*0.5+vt[i][j]*(vt[i+1][j]-vt[i-1][j])*dxi*0.5+vrc*vt[i][j]; temp2[i][j]=vt[i][j]+dt*(-exx*(vtn+(p[i][j]-p[i-1][j])*dxi)+exx*exx*rei*vtl); } } for(j=1;j<ny-1;j++) { for(i=1;i<nx-1;i++) { vr[i][j]=temp1[i][j]; } } for(j=1;j<ny-1;j++) { for(i=1;i<nx-1;i++) { vt[i][j]=temp2[i][j]; } } if(l%200==0) { id1<<"VARIABLES = \"X\",\"Y\",\"U\",\"V\",\"P\" "<<endl; id1<<"ZONE T=\"a\", I= "; id1<<nx-1; id1<<" J= "; id1<<ny; id1<<" f=point"<<endl; for(j=0;j<ny;j++) { for(i=0;i<nx-1;i++) { tet=dx*i; velo[0]=vr[i][j]*cos(tet)-vt[i][j]*sin(tet); velo[1]=vr[i][j]*sin(tet)+vt[i][j]*cos(tet); id1.precision(10); id1.setf(ios::showpoint); id1<<x[i][j]<<' '<<setw(20)<<y[i][j]<<' '<<setw(20)<<velo[0]<<' '<<setw(20)<<velo[1]<<' '<<setw(20)<<p[i][j]<<endl; } } } } }