C语言求解欧拉方程组
以下是我的需要写代码的公式一下是代码
程序代码：//改进的Lax-Wendroff的MacCormack格式.cpp
//利用差分格式求解二维欧拉方程问题
//date 20160725-0801
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#define GAMA 1.4//气体常数
#define MIN(x,y) (((x)<(y))?(x):(y))
#define L 0.05//涡核尺寸
#define Lx (20.0*L)//计算区域
#define Ly (20.0*L)
#define Tmax 3//总时间
#define Sf 0.8//时间步长因子
#define nx 40//网格数
#define ny 40
//全局变量
double rho[nx+2][ny+2], ux[nx+2][ny+2], uy[nx+2][ny+2], e[nx+2][ny+2];//rho(i,j,t)
double par_rho[nx+2][ny+2], par_ux[nx+2][ny+2], par_uy[nx+2][ny+2],par_e[nx+2][ny+2],par_p[nx+2][ny+2];//partial_rho(i,j,t)
double rho_ave[nx+2][ny+2], ux_ave[nx+2][ny+2], uy_ave[nx+2][ny+2], e_ave[nx+2][ny+2],p_ave[nx+2][ny+2];//rho_ave(i,j,t+dt)
double par_rho_ave[nx+2][ny+2], par_ux_ave[nx+2][ny+2],par_uy_ave[nx+2][ny+2], par_e_ave[nx+2][ny+2];//partial_rho_ave(i,j,t+dt)
double rho_t[nx+2][ny+2],ux_t[nx+2][ny+2],uy_t[nx+2][ny+2], e_t[nx+2][ny+2];//rho(i,j,t+dt)
double dx=Lx/nx,dy=Ly/ny;
//函数声明
void Initial(double rho[nx+2][ny+2],double ux[nx+2][ny+2],double uy[nx+2][ny+2],double e[nx+2][ny+2]);
double Cal_time_step(double rho[nx+2][ny+2],double ux[nx+2][ny+2],double uy[nx+2][ny+2],double e[nx+2][ny+2]);
void Mac_Cormack_2DSolver(double rho[nx+2][ny+2],double ux[nx+2][ny+2],double uy[nx+2][ny+2],double e[nx+2][ny+2],double dt);
void boundary(double rho_t[nx+2][ny+2],double ux_t[nx+2][ny+2],double uy_t[nx+2][ny+2],double e_t[nx+2][ny+2]);
void output(double rho[nx+2][ny+2],double ux[nx+2][ny+2],double uy[nx+2][ny+2],double e[nx+2][ny+2]);
//主函数
void main()
{
double t,dt;
double rho1,u1,v1,p1;
int i,j,w=0;
Initial(rho, ux, uy, e); //初始化
FILE *fp1;
fp1=fopen("initial_result.txt","w");
for(i=0;i<=nx+1;i++)
for(j=0;j<=ny+1;j++)
{
rho1=rho[i][j];
u1=ux[i][j];
v1=uy[i][j];
p1=(GAMA-1)*(e[i][j]-0.5*rho1*(u1*u1+v1*v1));
fprintf(fp1,"%e %e %e %e %e %e %e \n",i*dx,j*dy,rho1,u1,v1,p1,e[i][j]);
}
fclose(fp1);
t=0;
while(t<Tmax)
{
//dt=0.01;
dt=Cal_time_step(rho,ux, uy, e);//计算时间步长dt
t+=dt;
printf("t=%10g ,dt=%10g\n",t,dt);
Mac_Cormack_2DSolver( rho, ux,uy,e,dt);//二维差分格式
//exit(0);
w++;
}
output(rho,ux,uy, e);//数据保存
}
/*-----------------------------------------------------------------
初始化;
--------------------------------------------------------------------*/
void Initial(double rho[nx+2][ny+2],double ux[nx+2][ny+2],double uy[nx+2][ny+2],double e[nx+2][ny+2] )//初始化
{
int i,j;
double pinf=1.0, rouinf=1.0, epsllion=0.3, alpha=0.204;//初始条件
double x0=Lx/2.0, y0=Ly/2.0, r0, tau,theta,eps=0.0;
for(j=0;j<=ny+1;j++)
for(i=0;i<=nx+1;i++)
{
theta=atan((j*dy-y0)/(i*dx-x0));
r0=sqrt((i*dx-x0)*(i*dx-x0)+(j*dy-y0)*(j*dy-y0));
tau=r0/L;
rho[i][j]=rouinf*(pow((1-(GAMA-1)/(4*alpha*GAMA)*epsllion*epsllion*(exp(2*alpha*(1-tau*tau)))),(1/(GAMA-1)))-1)+eps;
ux[i][j]= epsllion*tau*exp(alpha*(1-tau*tau))*sin(theta)+eps;
uy[i][j]=-epsllion*tau*exp(alpha*(1-tau*tau))*cos(theta)+eps;
e[i][j]=rho[i][j]/(GAMA-1)+1/2*rho[i][j]*(ux[i][j]*ux[i][j]+uy[i][j]*uy[i][j])+eps;
//printf("%e\t%e\t%e\t%e\t%e\t%e\t%e\n\n",theta,r0,tau,rho[i][j],ux[i][j],uy[i][j],e[i][j]);
}
printf("初始化已完成\n\n");
}
/*--------------------------------------------------------------------
计算时间步长返回: 时间步长。
--------------------------------------------------------------------*/
double Cal_time_step(double rho[nx+2][ny+2],double ux[nx+2][ny+2],double uy[nx+2][ny+2],double e[nx+2][ny+2])
{
int i,j;
double maxvel,p0,u0,v0,vel,rho0,e0;
maxvel=1e-50;
for(i=1;i<=nx;i++)
for(j=1;j<=ny;j++)
{
//printf("rho=%e %e %e %e\n",rho[i][j],ux[i][j],uy[i][j],e[i][j]);
rho0=rho[i][j];
u0=ux[i][j];
v0=uy[i][j];
e0=e[i][j];
p0=(GAMA-1)*(e0-0.5*rho0*(u0*u0+v0*v0));
vel=sqrt(fabs(GAMA*p0/rho0))+sqrt(u0*u0+v0*v0); //CFL准则
// printf("%e %e\n",p0,vel);
if(vel>maxvel) maxvel=vel;
//printf("%e\t%e\t%e\t%e\t%e\n",rho0,u0,v0,p0,vel);
}
return Sf*MIN(dx,dy)/maxvel;
}
/*--------------------------------------------------------------------
差分格式
--------------------------------------------------------------------*/
void Mac_Cormack_2DSolver(double rho[nx+2][ny+2],double ux[nx+2][ny+2],double uy[nx+2][ny+2],double e[nx+2][ny+2],double dt)
{
int i,j;
double A1,A2,A3,A4,B1,B2,B3,C1,C2,C3,D1,D2,D3,D4,pi0j0,pi0j1,pi1j0,pi1j1;
//预测值
for(i=1;i<nx+1;i++)
{
for(j=1;j<ny+1;j++)
{
//printf("rho=%e %e %e %e\n",rho[i][j],ux[i][j],uy[i][j],e[i][j]);}}
pi0j0=(GAMA-1)*(e[i][j]-0.5*rho[i][j]*ux[i][j]*ux[i][j]);
pi0j1=(GAMA-1)*(e[i][j+1]-0.5*rho[i][j+1]*ux[i][j+1]*ux[i][j+1]);
pi1j0=(GAMA-1)*(e[i+1][j]-0.5*rho[i+1][j]*ux[i+1][j]*ux[i+1][j]);
pi1j1=(GAMA-1)*(e[i+1][j+1]-0.5*rho[i+1][j+1]*ux[i+1][j+1]*ux[i+1][j+1]);
A1=rho[i][j]*(ux[i+1][j]-ux[i][j])/dx;
A2=ux[i][j]*(rho[i+1][j]-rho[i][j])/dx;
A3=rho[i][j]*(uy[i][j+1]-uy[i][j])/dy;
A4=uy[i][j]*(rho[i][j+1]-rho[i][j])/dy;
par_rho[i][j]=-(A1+A2+A3+A4);
B1=ux[i][j]*(ux[i+1][j]-ux[i][j])/dx;
B2=uy[i][j]*(ux[i][j+1]-ux[i][j])/dy;
B3=1.0/rho[i][j]*(pi1j0-pi0j0)/dx;
par_ux[i][j]=-(B1+B2+B3);
C1=ux[i][j]*(uy[i+1][j]-uy[i][j])/dx;
C2=uy[i][j]*(uy[i][j+1]-uy[i][j])/dy;
C3=1.0/rho[i][j]*(pi0j1-pi0j0)/dy;
par_uy[i][j]=-(C1+C2+C3);
D1=ux[i][j]*(e[i+1][j]-e[i][j])/dx;
D2=uy[i][j]*(e[i][j+1]-e[i][j])/dy;
D3=pi0j0/rho[i][j]*(ux[i+1][j]-ux[i][j])/dx;
D4=pi0j0/rho[i][j]*(uy[i][j+1]-uy[i][j])/dy;
par_e[i][j]=-(D1+D2+D3+D4);
//printf("par_rho=%e,ux=%e,uy=%e,e=%e\n",par_rho[i][j],par_ux[i][j],par_uy[i][j],par_e[i][j]);
}
}
/*校正值，暂时未考虑人工黏滞项*/
for(i=1;i<nx+1;i++)
for(j=1;j<ny+1;j++)
{
rho_ave[i][j]=rho[i][j]+par_rho[i][j]*dt;//粘滞性加在这里
ux_ave[i][j]= ux[i][j]+ par_ux[i][j]*dt;//粘滞性加在这里
uy_ave[i][j]= uy[i][j]+ par_uy[i][j]*dt;//粘滞性加在这里
e_ave[i][j]= e[i][j]+ par_e[i][j]*dt;//粘滞性加在这里
p_ave[i][j]=(GAMA-1)*(e_ave[i][j]-0.5*rho_ave[i][j]*ux_ave[i][j]*ux_ave[i][j]);
}
double F1,F2,F3,F4,G1,G2,G3,H1,H2,H3,J1,J2,J3,J4;
for(i=1;i<nx+1;i++)
for(j=1;j<ny+1;j++)
{
F1=rho_ave[i][j]*(ux_ave[i][j]-ux_ave[i-1][j])/dx;
F2=ux_ave[i][j]*(rho_ave[i][j]-rho_ave[i-1][j])/dx;
F3=rho_ave[i][j]*(uy_ave[i][j]-uy_ave[i][j-1])/dy;
F4=uy_ave[i][j]*(rho_ave[i][j]-rho_ave[i][j-1])/dy;
par_rho_ave[i][j]=-(F1+F2+F3+F4);
G1=ux_ave[i][j]*(ux_ave[i][j]-ux_ave[i-1][j])/dx;
G2=uy_ave[i][j]*(ux_ave[i][j]-ux_ave[i][j-1])/dy;
G3=1.0/rho_ave[i][j]*(p_ave[i][j]-p_ave[i-1][j])/dx;
par_ux_ave[i][j]=-(G1+G2+G3);
H1=ux_ave[i][j]*(uy_ave[i][j]-uy_ave[i-1][j])/dx;
H2=uy_ave[i][j]*(uy_ave[i][j]-uy_ave[i][j-1])/dy;
H3=1.0/rho_ave[i][j]*(p_ave[i][j]-p_ave[i][j-1])/dy;
par_uy_ave[i][j]=-(H1+H2+H3);
J1=ux_ave[i][j]*(e_ave[i][j]-e_ave[i-1][j])/dx;
J2=uy_ave[i][j]*(e_ave[i][j]-e_ave[i][j-1])/dy;
J3=p_ave[i][j]/rho_ave[i][j]*(ux_ave[i][j]-ux_ave[i-1][j])/dx;
J4=p_ave[i][j]/rho_ave[i][j]*(uy_ave[i][j]-uy_ave[i][j-1])/dy;
par_e_ave[i][j]=-(J1+J2+J3+J4);
}
//预测加校正
for(i=1;i<nx+1;i++)
for(j=1;j<ny+1;j++)
{
rho_t[i][j]=rho[i][j]+0.5*(par_rho[i][j]+par_rho_ave[i][j]);
ux_t[i][j]= ux[i][j]+0.5*( par_ux[i][j]+ par_ux_ave[i][j]);
uy_t[i][j]= uy[i][j]+0.5*( par_uy[i][j]+ par_uy_ave[i][j]);
e_t[i][j]= e[i][j]+0.5*( par_e[i][j]+ par_e_ave[i][j]);
}
//加入边界条件
boundary( rho_t, ux_t, uy_t, e_t);
//把t+dt时刻值赋给t，做为下一时刻的起始值
for(i=0;i<=nx+1;i++)
for(j=0;j<=ny+1;j++)
{
rho[i][j]=rho_t[i][j];
ux[i][j]= ux_t[i][j];
uy[i][j]= uy_t[i][j];
e[i][j]= e_t[i][j];
}
}
/*--------------------------------------------------------------------
边界条件
--------------------------------------------------------------------*/
void boundary(double rho_t[nx+2][ny+2],double ux_t[nx+2][ny+2],double uy_t[nx+2][ny+2],double e_t[nx+2][ny+2])
{
int i,j;
//左边界
for(j=0;j<=nx+1;j++)
{
i=0;
rho_t[i][j]=2.0*rho_t[1][j]-rho_t[2][j];
ux_t[i][j]= 2.0*ux_t[1][j]- ux_t[2][j];
uy_t[i][j]= 2.0*uy_t[1][j]- uy_t[2][j];
e_t[i][j]= 2.0*e_t[1][j]- e_t[2][j];
}
//右边界
for(j=0;j<=nx+1;j++)
{
i=nx+1;
rho_t[i][j]=2.0*rho_t[nx][j]-rho_t[nx-1][j];
ux_t[i][j]= 2.0*ux_t[nx][j]- ux_t[nx-1][j];
uy_t[i][j]= 2.0*uy_t[nx][j]- uy_t[nx-1][j];
e_t[i][j]= 2.0*e_t[nx][j]- e_t[nx-1][j];
}
//下边界
for(i=0;i<=ny+1;i++)
{
j=0;
rho_t[i][j]=2.0*rho_t[i][1]-rho_t[i][2];
ux_t[i][j]= 2.0*ux_t[i][1]- ux_t[i][2];
uy_t[i][j]= 2.0*uy_t[i][1]- uy_t[i][2];
e_t[i][j]= 2.0*e_t[i][1]- e_t[i][2];
}
//上边界
for(i=0;i<=ny+1;i++)
{
j=ny+1;
rho_t[i][j]=2.0*rho_t[i][ny]-rho_t[i][ny-1];
ux_t[i][j]= 2.0*ux_t[i][ny]- ux_t[i][ny-1];
uy_t[i][j]= 2.0*uy_t[i][ny]- uy_t[i][ny-1];
e_t[i][j]= 2.0*e_t[i][ny]- e_t[i][ny-1];
}
}
/*-------------------------------------------------------
输出文件格式数据
---------------------------------------------------------*/
void output(double rho[nx+2][ny+2],double ux[nx+2][ny+2],double uy[nx+2][ny+2],double e[nx+2][ny+2])
{
int i,j;
FILE *fp;
double rho1,u1,v1,p1;
fp=fopen("result0.txt","w");
for(i=0;i<=nx+1;i++)
for(j=0;j<=ny+1;j++)
{
rho1=rho[i][j];
u1=ux[i][j];
v1=uy[i][j];
p1=(GAMA-1)*(e[i][j]-0.5*rho1*(u1*u1+v1*v1));
fprintf(fp,"%20f %20f %20.10e %20.10e %20.10e %20.10e %20.10e \n",i*dx,j*dy,rho1,u1,v1,p1,e[i][j]);
}
fclose(fp);
}
不知道那里有误？帮忙看一下，可以运行结果不对。
