椭圆型偏微分方程生成O型网格_流场_流场网格_网格；偏微分方程_翼型网格

#include<stdio.h> #include<math.h> #include<stdlib.h> int i=0,j=0,N,M; double PI=3.1415926,theta[1001],alpha,beta,gama,xi=1.0,eta=1.0,vel=5.0,xerror_max,yerror_max,J,ferror_max,e=1.E-7; double X[1001][1001],Y[1001][1001],XN[1001][1001],YN[1001][1001],Xerror[1001][1001],Yerror[1001][1001],F[1001][1001],FN[1001][1001],Ferror[1001][1001],U[1001][1001],V[1001][1001],P[1001][1001]; FILE *OUT; FILE *CP; FILE *XCP; void main() { printf("请输入X和Y方向的网格数（不超过1000）:N M\n"); scanf("%d%d",&N,&M); for(i=0;i<N+1;i++) { theta[i]=2*PI/N*i; X[i][0]=0.5*(1+cos(theta[i])); if(i<N/2) Y[i][0]=0.6*(-0.1015*pow((X[i][0]),4.0)+0.2843*pow((X[i][0]),3.0)-0.3576*pow((X[i][0]),2.0)-0.1221*X[i][0]+0.2969*pow((X[i][0]),0.5)); else Y[i][0]=-0.6*(-0.1015*pow((X[i][0]),4.0)+0.2843*pow((X[i][0]),3.0)-0.3576*pow((X[i][0]),2.0)-0.1221*X[i][0]+0.2969*pow((X[i][0]),0.5)); }//给定物面边界条件; for(i=0;i<N+1;i++) { theta[i]=2*PI/N*i; X[i][M]=30*cos(theta[i]); Y[i][M]=30*sin(theta[i]); }//给定远场边界条件; for(j=0;j<M+1;j++) for(i=0;i<N+1;i++) { X[i][j]=X[i][0]+j*(X[i][M]-X[i][0])/M; Y[i][j]=Y[i][0]+j*(Y[i][M]-Y[i][0])/M; }//插值初始化; for(;;) { for(j=1;j<M;j++) for(i=0;i<N+1;i++) { XN[i][j]=X[i][j]; YN[i][j]=Y[i][j]; if(i==0) { alpha=pow((0.5*(X[i][j+1]-X[i][j-1])/eta),2.0)+pow((0.5*(Y[i][j+1]-Y[i][j-1])/eta),2.0); beta=(0.5*(X[i+1][j]-X[N-1][j])/xi)*(0.5*(X[i][j+1]-X[i][j-1])/eta)+(0.5*(Y[i+1][j]-Y[N-1][j])/xi)*(0.5*(Y[i][j+1]-Y[i][j-1])/eta); gama=pow((0.5*(X[i+1][j]-X[N-1][j])/xi),2.0)+pow((0.5*(Y[i+1][j]-Y[N-1][j])/xi),2.0); X[i][j]=0.5*(alpha*(X[i+1][j]+X[N-1][j])+gama*(X[i][j+1]+X[i][j-1])-0.5*beta*(X[i+1][j+1]+X[N-1][j-1]-X[i+1][j-1]-X[N-1][j+1]))/(alpha+gama); Y[i][j]=0.5*(alpha*(Y[i+1][j]+Y[N-1][j])+gama*(Y[i][j+1]+Y[i][j-1])-0.5*beta*(Y[i+1][j+1]+Y[N-1][j-1]-Y[i+1][j-1]-Y[N-1][j+1]))/(alpha+gama); } else if(i==N) { alpha=pow((0.5*(X[i][j+1]-X[i][j-1])/eta),2.0)+pow((0.5*(Y[i][j+1]-Y[i][j-1])/eta),2.0); beta=(0.5*(X[1][j]-X[i-1][j])/xi)*(0.5*(X[i][j+1]-X[i][j-1])/eta)+(0.5*(Y[1][j]-Y[i-1][j])/xi)*(0.5*(Y[i][j+1]-Y[i][j-1])/eta); gama=pow((0.5*(X[1][j]-X[i-1][j])/xi),2.0)+pow((0.5*(Y[1][j]-Y[i-1][j])/xi),2.0); X[i][j]=0.5*(alpha*(X[1][j]+X[i-1][j])+gama*(X[i][j+1]+X[i][j-1])-0.5*beta*(X[1][j+1]+X[i-1][j-1]-X[1][j-1]-X[i-1][j+1]))/(alpha+gama); Y[i][j]=0.5*(alpha*(Y[1][j]+Y[i-1][j])+gama*(Y[i][j+1]+Y[i][j-1])-0.5*beta*(Y[1][j+1]+Y[i-1][j-1]-Y[1][j-1]-Y[i-1][j+1]))/(alpha+gama); } else { alpha=pow((0.5*(X[i][j+1]-X[i][j-1])/eta),2.0)+pow((0.5*(Y[i][j+1]-Y[i][j-1])/eta),2.0); beta=(0.5*(X[i+1][j]-X[i-1][j])/xi)*(0.5*(X[i][j+1]-X[i][j-1])/eta)+(0.5*(Y[i+1][j]-Y[i-1][j])/xi)*(0.5*(Y[i][j+1]-Y[i][j-1])/eta); gama=pow((0.5*(X[i+1][j]-X[i-1][j])/xi),2.0)+pow((0.5*(Y[i+1][j]-Y[i-1][j])/xi),2.0); X[i][j]=0.5*(alpha*(X[i+1][j]+X[i-1][j])+gama*(X[i][j+1]+X[i][j-1])-0.5*beta*(X[i+1][j+1]+X[i-1][j-1]-X[i+1][j-1]-X[i-1][j+1]))/(alpha+gama); Y[i][j]=0.5*(alpha*(Y[i+1][j]+Y[i-1][j])+gama*(Y[i][j+1]+Y[i][j-1])-0.5*beta*(Y[i+1][j+1]+Y[i-1][j-1]-Y[i+1][j-1]-Y[i-1][j+1]))/(alpha+gama); } } for(j=1;j<M;j++) for(i=0;i<N+1;i++) { Xerror[i][j]=fabs(X[i][j]-XN[i][j]); Yerror[i][j]=fabs(Y[i][j]-YN[i][j]); } xerror_max=Xerror[0][0]; yerror_max=Yerror[0][0]; for(j=1;j<M;j++) for(i=0;i<N+1;i++) { if(Xerror[i][j]>=xerror_max) xerror_max=Xerror[i][j]; if(Yerror[i][j]>=yerror_max) yerror_max=Yerror[i][j]; } if(xerror_max<=e && yerror_max<=e) break; }//网格生成; OUT=fopen("grid.dat","wb+"); fprintf(OUT,"TITLE=\"2D grid\"\r\nVARIABLES=\"X\"\,\"Y\"\r\nZONE I=%d J=%d\r\n",N+1,M+1); for(j=0;j<M+1;j++) for(i=0;i<N+1;i++) fprintf(OUT,"%lf %lf\r\n",X[i][j],Y[i][j]); printf("网格已生成，请查看输出文件grid.dat\n"); fclose(OUT);//网格文件输出; for(i=0;i<N+1;i++) { F[i][0]=0; F[i][M]=vel*X[i][M]; } for(j=1;j<M;j++) for(i=0;i<N+1;i++) F[i][j]=F[i][0]+j*(F[i][M]-F[i][0])/M; //速度势初始化; for(;;) { for(i=0;i<N+1;i++) for(j=0;j<M;j++) { FN[i][j]=F[i][j]; if(j==0) { if(i==0) { F[i-1][j]=F[N-1][j]; X[i-1][j]=X[N-1][j]; Y[i-1][j]=Y[N-1][j]; beta=((X[i+1][j]-X[i-1][j])/2)*(((-3)*X[i][j]+4*X[i][j+1]-X[i][j+2])/2)+((Y[i+1][j]-Y[i-1][j])/2)*(((-3)*Y[i][j]+4*Y[i][j+1]-Y[i][j+2])/2); gama=pow((X[i+1][j]-X[i-1][j])/2,2)+pow((Y[i+1][j]-Y[i-1][j])/2,2); F[i][j]=((F[i+1][j]-F[i-1][j])*beta/gama-(4*F[i][j+1]-F[i][j+2]))/(-3); } else if(i==N) { F[i+1][j]=F[1][j]; X[i+1][j]=X[1][j]; Y[i+1][j]=Y[1][j]; beta=((X[i+1][j]-X[i-1][j])/2)*(((-3)*X[i][j]+4*X[i][j+1]-X[i][j+2])/2)+((Y[i+1][j]-Y[i-1][j])/2)*(((-3)*Y[i][j]+4*Y[i][j+1]-Y[i][j+2])/2); gama=pow((X[i+1][j]-X[i-1][j])/2,2)+pow((Y[i+1][j]-Y[i-1][j])/2,2); F[i][j]=((F[i+1][j]-F[i-1][j])*beta/gama-(4*F[i][j+1]-F[i][j+2]))/(-3); } else { beta=((X[i+1][j]-X[i-1][j])/2)*(((-3)*X[i][j]+4*X[i][j+1]-X[i][j+2])/2)+((Y[i+1][j]-Y[i-1][j])/2)*(((-3)*Y[i][j]+4*Y[i][j+1]-Y[i][j+2])/2); gama=pow((X[i+1][j]-X[i-1][j])/2,2)+pow((Y[i+1][j]-Y[i-1][j])/2,2); F[i][j]=((F[i+1][j]-F[i-1][j])*beta/gama-(4*F[i][j+1]-F[i][j+2]))/(-3); } } else { if(i==0) { F[i-1][j]=F[N-1][j]; F[i-1][j-1]=F[N-1][j-1]; F[i-1][j+1]=F[N-1][j+1]; X[i-1][j]=X[N-1][j]; Y[i-1][j]=Y[N-1][j]; alpha=pow((X[i][j+1]-X[i][j-1])/2,2)+pow((Y[i][j+1]-Y[i][j-1])/2,2); beta=((X[i+1][j]-X[i-1][j])/2)*((X[i][j+1]-X[i][j-1])/2)+((Y[i+1][j]-Y[i-1][j])/2)*((Y[i][j+1]-Y[i][j-1])/2); gama=pow((X[i+1][j]-X[i-1][j])/2,2)+pow((Y[i+1][j]-Y[i-1][j])/2,2); F[i][j]=0.5*(alpha*(F[i+1][j]+F[i-1][j])+gama*(F[i][j+1]+F[i][j-1])-0.5*beta*(F[i+1][j+1]+F[i-1][j-1]-F[i+1][j-1]-F[i-1][j+1]))/(alpha+gama); } else if(i==N) { F[i+1][j]=F[1][j]; F[i+1][j-1]=F[1][j-1]; F[i+1][j+1]=F[1][j+1]; X[i+1][j]=X[1][j]; Y[i+1][j]=Y[1][j]; alpha=pow((X[i][j+1]-X[i][j-1])/2,2)+pow((Y[i][j+1]-Y[i][j-1])/2,2); beta=((X[i+1][j]-X[i-1][j])/2)*((X[i][j+1]-X[i][j-1])/2)+((Y[i+1][j]-Y[i-1][j])/2)*((Y[i][j+1]-Y[i][j-1])/2); gama=pow((X[i+1][j]-X[i-1][j])/2,2)+pow((Y[i+1][j]-Y[i-1][j])/2,2); F[i][j]=0.5*(alpha*(F[i+1][j]+F[i-1][j])+gama*(F[i][j+1]+F[i][j-1])-0.5*beta*(F[i+1][j+1]+F[i-1][j-1]-F[i+1][j-1]-F[i-1][j+1]))/(alpha+gama); } else { alpha=pow((X[i][j+1]-X[i][j-1])/2,2)+pow((Y[i][j+1]-Y[i][j-1])/2,2); beta=((X[i+1][j]-X[i-1][j])/2)*((X[i][j+1]-X[i][j-1])/2)+((Y[i+1][j]-Y[i-1][j])/2)*((Y[i][j+1]-Y[i][j-1])/2); gama=pow((X[i+1][j]-X[i-1][j])/2,2)+pow((Y[i+1][j]-Y[i-1][j])/2,2); F[i][j]=0.5*(alpha*(F[i+1][j]+F[i-1][j])+gama*(F[i][j+1]+F[i][j-1])-0.5*beta*(F[i+1][j+1]+F[i-1][j-1]-F[i+1][j-1]-F[i-1][j+1]))/(alpha+gama); } } } for(j=0;j<M;j++) for(i=0;i<N+1;i++) { Ferror[i][j]=fabs(F[i][j]-FN[i][j]); } ferror_max=Ferror[0][0]; for(j=0;j<M;j++) for(i=0;i<N+1;i++) { if(Ferror[i][j]>=ferror_max) ferror_max=Ferror[i][j]; } if(ferror_max<=e) break; }//速度势求解完毕; for(i=0;i<N+1;i++) { U[i][M]=vel; V[i][M]=0; } for