红色是我加入的
这时我上面说的另一个程序！
#include<stdio.h>
#include<math.h>
void main()
{
double x[102][62],y[102][62];
int i,j,iter;
int im=101,jm=61;
double xe[102][62],xf[102][62],ye[102][62],yf[102][62];
double ex[102][62],fx[102][62],ey[102][62],fy[102][62];
double delta1=1.0,delta2=1.0; /*步长*/
double x1,y1;
double error1,error2;
double a,b,c;
double jacobi[102][62];
int iternum,type;
double k,r,ma0,p0,d0,u0,v0;
double p[102][62][4],d[102][62][4],u[102][62][4],v[102][62][4],e[102][62][4];
printf("出口为超音输入1，出口为亚音输入2\n");
printf("请输入出口的类型：");
scanf("%d",&type);
FILE*xy;
xy=fopen("xy.txt","w");
FILE*xe1,*xf1,*ye1,*yf1,*ex1,*fx1,*ey1,*fy1;
xe1=fopen("xe.txt","w");
xf1=fopen("xf.txt","w");
ye1=fopen("ye.txt","w");
yf1=fopen("yf.txt","w");
ex1=fopen("ex.txt","w");
fx1=fopen("fx.txt","w");
ey1=fopen("ey.txt","w");
fy1=fopen("fy.txt","w");
/*生成网格*/
for(i=1;i<=im;i++)
{
for(j=(jm+1)/2;j<=jm;j++)
{
x[i][j]=10.0/(im-1)*(i-1);
y[i][j]=(j-(jm+1)/2)*(1.398+0.347*tanh(0.8*x[i][j]-4))/(j-1);
}
}
for(i=1;i<=im;i++)
{
for(j=1;j<=jm/2;j++)
{
x[i][j]=x[i][jm+1-j];
y[i][j]=-y[i][jm+1-j];
}
}
iter=1;
do
{
error1=0.0;
for(i=2;i<=im-1;i++)
{
for(j=2;j<=jm-1;j++)
{
xe[i][j]=(x[i+1][j]-x[i-1][j])/(2*delta1);
xf[i][j]=(x[i][j+1]-x[i][j-1])/(2*delta2);
ye[i][j]=(y[i+1][j]-y[i-1][j])/(2*delta1);
yf[i][j]=(y[i][j+1]-y[i][j-1])/(2*delta2);
a=xf[i][j]*xf[i][j]+yf[i][j]*yf[i][j];
b=xe[i][j]*xf[i][j]+ye[i][j]*yf[i][j];
c=xe[i][j]*xe[i][j]+ye[i][j]*ye[i][j];
x1=(2*a*delta2*delta2*(x[i+1][j]+x[i-1][j])
-b*delta1*delta2*(x[i+1][j+1]-x[i-1][j+1]-x[i+1][j-1]+x[i-1][j-1])
+2*c*delta1*delta1*(x[i][j+1]+x[i][j-1]))
/(4*(a*delta2*delta2+c*delta1*delta1));
y1=(2*a*delta2*delta2*(y[i+1][j]+y[i-1][j])
-b*delta1*delta2*(y[i+1][j+1]-y[i-1][j+1]-y[i+1][j-1]+y[i-1][j-1])
+2*c*delta1*delta1*(y[i][j+1]+y[i][j-1]))
/(4*(a*delta2*delta2+c*delta1*delta1));
error2=fabs(x[i][j]-x1)+fabs(y[i][j]-y1);
if(error2>error1)
error1=error2;
x[i][j]=x1;
y[i][j]=y1;
}
}
iter=iter+1;
if(iter>10000)break;
}
while(error1>1e-4);
for(i=1;i<=im;i++)
{
for(j=1;j<=jm;j++)
{
fprintf(xy,"%lf, %lf\n",x[i][j],y[i][j]);
}
}
/*网格生成完毕*/
/*网格参数的计算*/
/*中间网格的参数*/
for(i=2;i<=im-1;i++)
{
for(j=2;j<=jm-1;j++)
{
xe[i][j]=(x[i+1][j]-x[i-1][j])/(2*delta1);
xf[i][j]=(x[i][j+1]-x[i][j-1])/(2*delta2);
ye[i][j]=(y[i+1][j]-y[i-1][j])/(2*delta1);
yf[i][j]=(y[i][j+1]-y[i][j-1])/(2*delta2);
}
}
/*边界网格的参数*/
/***下边界***/
xe[1][1]=(x[2][1]-x[1][1])/delta1;
ye[1][1]=(y[2][1]-y[1][1])/delta1;
xf[1][1]=(x[1][2]-x[1][1])/delta2;
yf[1][1]=(y[1][2]-y[1][1])/delta2;
for(i=2;i<=im-1;i++)
{
xe[i][1]=(x[i+1][1]-x[i][1])/(2*delta1);
ye[i][1]=(y[i+1][1]-y[i][1])/(2*delta1);
xf[i][1]=(x[i][2]-x[i][1])/delta2;
yf[i][1]=(y[i][2]-y[i][1])/delta2;
}
xe[im][1]=(x[im][1]-x[im-1][1])/delta1;
ye[im][1]=(y[im][1]-y[im-1][1])/delta1;
xf[im][1]=(x[im][2]-x[im][1])/delta2;
yf[im][1]=(y[im][2]-y[im][1])/delta2;
/***上边界***/
xe[1][jm]=(x[2][jm]-x[1][jm])/delta1;
ye[1][jm]=(y[2][jm]-y[1][jm])/delta1;
xf[1][jm]=(x[1][jm]-x[1][jm-1])/delta2;
yf[1][jm]=(y[1][jm]-y[1][jm-1])/delta2;
for(i=2;i<=im-1;i++)
{
xe[i][jm]=(x[i+1][jm]-x[i][jm])/(2*delta1);
ye[i][jm]=(y[i+1][jm]-y[i][jm])/(2*delta1);
xf[i][jm]=(x[i][jm]-x[i][jm-1])/delta2;
yf[i][jm]=(y[i][jm]-y[i][jm-1])/delta2;
}
xe[im][jm]=(x[im][jm]-x[im-1][jm])/delta1;
ye[im][jm]=(y[im][jm]-y[im-1][jm])/delta1;
xf[im][jm]=(x[im][jm]-x[im][jm-1])/delta2;
yf[im][jm]=(y[im][jm]-y[im][jm-1])/delta2;
/***进口边界***/
for(j=2;j<=jm-1;j++)
{
xe[1][j]=(x[2][j]-x[1][j])/delta1;
ye[1][j]=(y[2][j]-y[1][j])/delta1;
xf[1][j]=(x[1][j+1]-x[1][j-1])/(2*delta2);
yf[1][j]=(y[1][j+1]-y[1][j-1])/(2*delta2);
}
/***出口边界***/
for(j=2;j<=jm-1;j++)
{
xe[im][j]=(x[im][j]-x[im-1][j])/delta1;
ye[im][j]=(y[im][j]-y[im][j])/delta1;
xf[im][j]=(x[im][j+1]-x[im][j-1])/(2*delta2);
yf[im][j]=(y[im][j+1]-y[im][j-1])/(2*delta2);
}
for(i=1;i<=im;i++)
{
for(j=1;j<=jm;j++)
{
jacobi[i][j]=1.0/(xe[i][j]*yf[i][j]-xf[i][j]*ye[i][j]);
ex[i][j]=yf[i][j]/jacobi[i][j];
ey[i][j]=-xf[i][j]/jacobi[i][j];
fx[i][j]=-ye[i][j]/jacobi[i][j];
fy[i][j]=xe[i][j]/jacobi[i][j];
}
}
for(i=1;i<=im;i++)
{
for(j=1;j<=jm;j++)
{
fprintf(xe1,"%lf\n",xe[i][j]);
fprintf(xf1,"%lf\n",xf[i][j]);
fprintf(ye1,"%lf\n",ye[i][j]);
fprintf(yf1,"%lf\n",yf[i][j]);
fprintf(ex1,"%lf\n",ex[i][j]);
fprintf(fx1,"%lf\n",fx[i][j]);
fprintf(ey1,"%lf\n",ey[i][j]);
fprintf(fy1,"%lf\n",fy[i][j]);
}
}
/*网格参数计算完毕*/
/*迭代计算*/
k=1.4;
r=1.4;
ma0=1.5;
p0=47892.40;
d0=1.2218;
u0=ma0*sqrt(k*p0/d0);
v0=0;
iternum=1;
/*设置流场参数的初值*/
for(i=1;i<=im;i++)
{
for(j=1;j<=jm;j++)
{
p[i][j][1]=p0;
d[i][j][1]=d0;
u[i][j][1]=u0;
v[i][j][1]=v0;
e[i][j][1]=p[i][j][1]/(r-1)+0.5*d[i][j][1]/(pow(u[i][j][1],2)+pow(v[i][j][1],2));
}
}
}
谢谢