📄 muscl euler two dimentions.cpp
字号:
// Muscl Euler Two dimensions.cpp : Defines the entry point for the console application.
#include "stdafx.h"
#include "stdlib.h"
#include "math.h"
const int Nx=120, Ny=100;
double pi=atan(1.0)*4.0, gamma=1.4, R=287.06;
class CMyPoint
{
public:
double x;
double y;
};
double GetLength(CMyPoint a, CMyPoint b)
{
return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y));
}
void GenerateMesh(CMyPoint mesh[Nx+1][Ny+1])
{
int i=0, j=0;
double Len=10.0, High=14, delta=12.5*pi/180.0;
double LeftHigh=High+Len*tan(delta);
double xT[Nx+1], yT[Nx+1], xR[Ny+1], yR[Ny+1];
double xB[Nx+1], yB[Nx+1], xL[Ny+1], yL[Ny+1];
for (j=0; j<=Ny; j++)
{
xL[j]=0, yL[j]=j*LeftHigh/Ny;
xR[j]=2*Len, yR[j]=(LeftHigh-High)+j*High/Ny;
}
for (i=0; i<=Nx/2; i++)
{
xB[i]=i*2*Len/Nx;
yB[i]=xB[i]*tan(delta);
}
for (; i<=Nx; i++)
{
xB[i]=i*2*Len/Nx;
yB[i]=Len*tan(delta);
}
for (i=0; i<=Nx; i++)
{
xT[i]=i*2*Len/Nx;
yT[i]=LeftHigh;
}
for (i=0; i<=Nx; i++)
{
for (j=0; j<=Ny; j++)
{
double zeta=i*1.0/Nx, eta=j*1.0/Ny;
mesh[i][j].x=(1-eta)*xB[i]+eta*xT[i]+(1-zeta)*xL[j]+zeta*xR[j]
-( zeta*eta*xT[Nx]+zeta*(1-eta)*xB[Nx]+eta*(1-zeta)*xT[0]+(1-eta)*(1-zeta)*xB[0] );
mesh[i][j].y=(1-eta)*yB[i]+eta*yT[i]+(1-zeta)*yL[j]+zeta*yR[j]
-( zeta*eta*yT[Nx]+zeta*(1-eta)*yB[Nx]+eta*(1-zeta)*yT[0]+(1-eta)*(1-zeta)*yB[0] );
}
}
}
class CEuler
{
public:
CEuler();
double GetRho();
double GetU();
double GetV();
double GetH();
double GetP();
double GetArea();
void SetRho(double initRho);
void SetRhoU(double initRhoU);
void SetRhoV(double initRhoV);
void SetRhoE(double initRhoE);
double GetRhoU();
double GetRhoV();
double GetRhoE();
double GetLeftLen();
double GetRightLen();
double GetTopLen();
double GetBottomLen();
void SetWall();
CMyPoint NE;
CMyPoint NW;
CMyPoint SE;
CMyPoint SW;
CMyPoint nRight;
CMyPoint nTop;
double leftLen;
double rightLen;
double topLen;
double bottomLen;
private:
double rho;
double rhoU;
double rhoV;
double rhoE;
};
CEuler::CEuler()
{
}
double CEuler::GetRho()
{
return rho;
}
double CEuler::GetP()
{
return (gamma-1)*(rhoE-0.5*(rhoU*rhoU+rhoV*rhoV)/rho);
}
double CEuler::GetU()
{
return rhoU/rho;
}
double CEuler::GetV()
{
return rhoV/rho;
}
double CEuler::GetH()
{
return ( gamma*rhoE-(gamma-1)/2.0*(rhoU*rhoU+rhoV*rhoV) )/rho;// H=gamma*E-(gamma-1)/2*(u^2+v^2)*rho
}
double CEuler::GetRhoU()
{
return rhoU;
}
double CEuler::GetRhoV()
{
return rhoV;
}
double CEuler::GetRhoE()
{
return rhoE;
}
double CEuler::GetArea()
{
return fabs((SE.x-NW.x)*(NE.y-SW.y)-(SE.y-NW.y)*(NE.x-SW.x))/2.0;
}
double CEuler::GetLeftLen()
{
return GetLength(NW, SW);
}
double CEuler::GetRightLen()
{
return GetLength(NE, SE);
}
double CEuler::GetTopLen()
{
return GetLength(NE, NW);
}
double CEuler::GetBottomLen()
{
return GetLength(SE, SW);
}
void CEuler::SetRho(double initRho)
{
rho=initRho;
}
void CEuler::SetRhoU(double initRhoU)
{
rhoU=initRhoU;
}
void CEuler::SetRhoV(double initRhoV)
{
rhoV=initRhoV;
}
void CEuler::SetRhoE(double initRhoE)
{
rhoE=initRhoE;
}
void CEuler::SetWall()
{
double cosA=(SE.x-SW.x)/bottomLen;
double sinA=(SE.y-SW.y)/bottomLen;
double tempRhoU=rhoU*cosA*cosA+rhoV*cosA*sinA;
double tempRhoV=rhoU*sinA*cosA+rhoV*sinA*sinA;
rhoE=GetP()/(gamma-1)+(tempRhoU*tempRhoU+tempRhoV*tempRhoV)/rho/2.0;
rhoU=tempRhoU;
rhoV=tempRhoV;
}
double VanLeer(double a, double b, double c)
{
double left=b-a;
double right=c-b;
return (left*left+left*right)/(right*right+left*left+(1e-6))*right;
}
double coeffSHK[4]={0.0833, 0.2069, 0.4265, 1.0};
double tempRhoLeft[2][Nx][Ny], tempULeft[2][Nx][Ny], tempVLeft[2][Nx][Ny], tempHLeft[2][Nx][Ny];
double tempRhoRight[2][Nx][Ny], tempURight[2][Nx][Ny], tempVRight[2][Nx][Ny], tempHRight[2][Nx][Ny];
double tempRoeRho[2][Nx][Ny], tempRoeRhoU[2][Nx][Ny], tempRoeRhoV[2][Nx][Ny], tempRoeRhoE[2][Nx][Ny];
CEuler Q[Nx+1][Ny+1];
int _tmain(int argc, _TCHAR* argv[])
{
int i=0, j=0, k=0;
FILE *fp=NULL;
fp=fopen("Euler Values.dat", "w");
fprintf(fp, "VARIABLES=x, y, Pressure, Mach_Number, Enthalpy\n");
fprintf(fp, "ZONE, F=POINT, I=%d J=%d\n", Ny, Nx);
CMyPoint mesh[Nx+1][Ny+1];
GenerateMesh(mesh);
double Pinf=101325, Tinf=300.0, Minf=2.0;
double rho=Pinf/R/Tinf;
double rhoU=rho*Minf*sqrt(gamma*R*Tinf);
double rhoV=0;
double rhoE=Pinf*(1.0/(gamma-1)+gamma*Minf*Minf/2.0);
for (i=0; i<Nx; i++)
{
for (j=0; j<Ny; j++)
{
Q[i][j].SW=mesh[i][j];
Q[i][j].NW=mesh[i][j+1];
Q[i][j].NE=mesh[i+1][j+1];
Q[i][j].SE=mesh[i+1][j];
Q[i][j].SetRho(rho);
Q[i][j].SetRhoU(rhoU);
Q[i][j].SetRhoV(rhoV);
Q[i][j].SetRhoE(rhoE);
Q[i][j].leftLen=Q[i][j].GetLeftLen();
Q[i][j].rightLen=Q[i][j].GetRightLen();
Q[i][j].topLen=Q[i][j].GetTopLen();
Q[i][j].bottomLen=Q[i][j].GetBottomLen();
Q[i][j].nRight.x=-(Q[i][j].SE.y-Q[i][j].NE.y)/Q[i][j].rightLen;
Q[i][j].nRight.y=(Q[i][j].SE.x-Q[i][j].NE.x)/Q[i][j].rightLen;
Q[i][j].nTop.x=-(Q[i][j].NE.y-Q[i][j].NW.y)/Q[i][j].topLen;
Q[i][j].nTop.y=(Q[i][j].NE.x-Q[i][j].NW.x)/Q[i][j].topLen;
}
}
for (i=0; i<Nx; i++)
{
Q[i][0].SetWall();
}
int Nt=4000, t=0;
double dt=1, CFL=0.1;
for (int iter=0; iter<Nt; iter++)
{
printf("Iterate=%d\t\trho=%.6lf\n", iter, Q[20][10].GetRho());
for (i=0; i<Nx-1; i++)
{
for (j=0; j<Ny-1; j++)
{
tempRhoLeft[0][i][j]=Q[i][j].GetRho();
tempRhoRight[0][i][j]=Q[i+1][j].GetRho();
tempULeft[0][i][j]=Q[i][j].GetU();
tempURight[0][i][j]=Q[i+1][j].GetU();
tempVLeft[0][i][j]=Q[i][j].GetV();
tempVRight[0][i][j]=Q[i+1][j].GetV();
tempHLeft[0][i][j]=Q[i][j].GetH();
tempHRight[0][i][j]=Q[i+1][j].GetH();
tempRhoLeft[1][i][j]=Q[i][j].GetRho();
tempRhoRight[1][i][j]=Q[i][j+1].GetRho();
tempULeft[1][i][j]=Q[i][j].GetU();
tempURight[1][i][j]=Q[i][j+1].GetU();
tempVLeft[1][i][j]=Q[i][j].GetV();
tempVRight[1][i][j]=Q[i][j+1].GetV();
tempHLeft[1][i][j]=Q[i][j].GetH();
tempHRight[1][i][j]=Q[i][j+1].GetH();
}
}
double D[2], aU[2], aV[2], ac[2], aH[2], aRho[2], deltaRho[2], deltaU[2], deltaV[2], deltaP[2];
double beta[2][7], lamda[2][3], nx[2], ny[2];
for (i=0; i<Nx-1; i++)
{
for (j=0; j<Ny-1; j++)
{
nx[0]=Q[i][j].nRight.x, ny[0]=Q[i][j].nRight.y;
nx[1]=Q[i][j].nTop.x, ny[1]=Q[i][j].nTop.y;
for (k=0; k<2; k++)
{
aRho[k]=sqrt(tempRhoRight[k][i][j]*tempRhoLeft[k][i][j]);
D[k]=sqrt(tempRhoRight[k][i][j]/tempRhoLeft[k][i][j]);
aU[k]=(tempULeft[k][i][j]+D[k]*tempURight[k][i][j])/(1+D[k]);
aV[k]=(tempVLeft[k][i][j]+D[k]*tempVRight[k][i][j])/(1+D[k]);
aH[k]=(tempHLeft[k][i][j]+D[k]*tempHRight[k][i][j])/(1+D[k]);
ac[k]=sqrt((gamma-1)*(aH[k]-0.5*(aU[k]*aU[k]+aV[k]*aV[k])));
deltaRho[k]=tempRhoRight[k][i][j]-tempRhoLeft[k][i][j];
deltaU[k]=tempURight[k][i][j]-tempULeft[k][i][j];
deltaV[k]=tempVRight[k][i][j]-tempVLeft[k][i][j];
deltaP[k]=(gamma-1)/gamma*tempRhoRight[k][i][j]*(tempHRight[k][i][j]-0.5*(tempURight[k][i][j]*tempURight[k][i][j]+tempVRight[k][i][j]*tempVRight[k][i][j]))\
-(gamma-1)/gamma*tempRhoLeft[k][i][j]*(tempHLeft[k][i][j]-0.5*(tempULeft[k][i][j]*tempULeft[k][i][j]+tempVLeft[k][i][j]*tempVLeft[k][i][j]));
lamda[k][0]=nx[k]*aU[k]+ny[k]*aV[k];
lamda[k][1]=lamda[k][0]-ac[k];
lamda[k][2]=lamda[k][0]+ac[k];
beta[k][0]=fabs(lamda[k][0])*(deltaRho[k]-deltaP[k]/ac[k]/ac[k]);
beta[k][1]=fabs(lamda[k][2])*(deltaP[k]+aRho[k]*ac[k]*(nx[k]*deltaU[k]+ny[k]*deltaV[k]))/2.0/ac[k]/ac[k];
beta[k][2]=fabs(lamda[k][1])*(deltaP[k]-aRho[k]*ac[k]*(nx[k]*deltaU[k]+ny[k]*deltaV[k]))/2.0/ac[k]/ac[k];
beta[k][3]=beta[k][0]+beta[k][1]+beta[k][2];
beta[k][4]=ac[k]*(beta[k][1]-beta[k][2]);
beta[k][5]=fabs(lamda[k][0])*aRho[k]*(deltaU[k]-nx[k]*(nx[k]*deltaU[k]+ny[k]*deltaV[k]));
beta[k][6]=fabs(lamda[k][0])*aRho[k]*(deltaV[k]-ny[k]*(nx[k]*deltaU[k]+ny[k]*deltaV[k]));
tempRoeRho[k][i][j]=0.5*( tempRhoLeft[k][i][j]*(nx[k]*tempULeft[k][i][j]+ny[k]*tempVLeft[k][i][j])\
+tempRhoRight[k][i][j]*(nx[k]*tempURight[k][i][j]+ny[k]*tempVRight[k][i][j])\
-beta[k][3]);
tempRoeRhoU[k][i][j]=0.5*((tempRhoLeft[k][i][j]*(nx[k]*tempULeft[k][i][j]+ny[k]*tempVLeft[k][i][j]))*tempULeft[k][i][j]\
+nx[k]*(gamma-1)/gamma*tempRhoLeft[k][i][j]*(tempHLeft[k][i][j]-0.5*(tempULeft[k][i][j]*tempULeft[k][i][j]+tempVLeft[k][i][j]*tempVLeft[k][i][j]))\
+(tempRhoRight[k][i][j]*(nx[k]*tempURight[k][i][j]+ny[k]*tempVRight[k][i][j]))*tempURight[k][i][j]\
+nx[k]*(gamma-1)/gamma*tempRhoRight[k][i][j]*(tempHRight[k][i][j]-0.5*(tempURight[k][i][j]*tempURight[k][i][j]+tempVRight[k][i][j]*tempVRight[k][i][j])))\
-(aU[k]*beta[k][3]+nx[k]*beta[k][4]+beta[k][5])*0.5;
tempRoeRhoV[k][i][j]=0.5*((tempRhoLeft[k][i][j]*(nx[k]*tempULeft[k][i][j]+ny[k]*tempVLeft[k][i][j]))*tempVLeft[k][i][j]\
+ny[k]*(gamma-1)/gamma*tempRhoLeft[k][i][j]*(tempHLeft[k][i][j]-0.5*(tempULeft[k][i][j]*tempULeft[k][i][j]+tempVLeft[k][i][j]*tempVLeft[k][i][j]))\
+(tempRhoRight[k][i][j]*(nx[k]*tempURight[k][i][j]+ny[k]*tempVRight[k][i][j]))*tempVRight[k][i][j]\
+ny[k]*(gamma-1)/gamma*tempRhoRight[k][i][j]*(tempHRight[k][i][j]-0.5*(tempURight[k][i][j]*tempURight[k][i][j]+tempVRight[k][i][j]*tempVRight[k][i][j])))\
-(aV[k]*beta[k][3]+ny[k]*beta[k][4]+beta[k][6])*0.5;
tempRoeRhoE[k][i][j]=0.5*( tempRhoLeft[k][i][j]*(nx[k]*tempULeft[k][i][j]+ny[k]*tempVLeft[k][i][j])*tempHLeft[k][i][j]\
+tempRhoRight[k][i][j]*(nx[k]*tempURight[k][i][j]+ny[k]*tempVRight[k][i][j])*tempHRight[k][i][j]\
-(aH[k]*beta[k][3]+(nx[k]*aU[k]+ny[k]*aV[k])*beta[k][4]+aU[k]*beta[k][5]+aV[k]*beta[k][6]-ac[k]*ac[k]*beta[k][0]/(gamma-1)) );
}
}
}
for (t=0; t<4; t++)
{
for (i=1; i<Nx-1; i++)
{
for (j=1; j<Ny-1; j++)
{
dt=CFL*Q[i][j].bottomLen*(6e-4);
double temp=coeffSHK[t]*dt/Q[i][j].GetArea();
double sourse=coeffSHK[t]*dt/(Q[i][j].NE.y+Q[i][j].SE.y)*2.0;
rho= Q[i][j].GetRho()-temp*(Q[i][j].rightLen*tempRoeRho[0][i][j]-Q[i][j].leftLen*tempRoeRho[0][i-1][j]\
+Q[i][j].topLen*tempRoeRho[1][i][j]-Q[i][j].bottomLen*tempRoeRho[1][i][j-1])-sourse*Q[i][j].GetRhoV();
rhoU= Q[i][j].GetRhoU()-temp*(Q[i][j].rightLen*tempRoeRhoU[0][i][j]-Q[i][j].leftLen*tempRoeRhoU[0][i-1][j]\
+Q[i][j].topLen*tempRoeRhoU[1][i][j]-Q[i][j].bottomLen*tempRoeRhoU[1][i][j-1])-sourse*Q[i][j].GetRhoV()*Q[i][j].GetU();
rhoV= Q[i][j].GetRhoV()-temp*(Q[i][j].rightLen*tempRoeRhoV[0][i][j]-Q[i][j].leftLen*tempRoeRhoV[0][i-1][j]\
+Q[i][j].topLen*tempRoeRhoV[1][i][j]-Q[i][j].bottomLen*tempRoeRhoV[1][i][j-1])-sourse*Q[i][j].GetRhoV()*Q[i][j].GetV();
rhoE= Q[i][j].GetRhoE()-temp*(Q[i][j].rightLen*tempRoeRhoE[0][i][j]-Q[i][j].leftLen*tempRoeRhoE[0][i-1][j]\
+Q[i][j].topLen*tempRoeRhoE[1][i][j]-Q[i][j].bottomLen*tempRoeRhoE[1][i][j-1])-sourse*Q[i][j].GetRhoV()*Q[i][j].GetH();
Q[i][j].SetRho(rho);
Q[i][j].SetRhoU(rhoU);
Q[i][j].SetRhoV(rhoV);
Q[i][j].SetRhoE(rhoE);
}
}
for (i=0; i<Nx; i++)
{
for (j=0; j<1; j++)
{
Q[i][j].SetRho(Q[i][1].GetRho());
Q[i][j].SetRhoU(Q[i][1].GetRhoU());
Q[i][j].SetRhoV(Q[i][1].GetRhoV());
Q[i][j].SetRhoE(Q[i][1].GetRhoE());
}
Q[i][0].SetWall();
}
for (j=0; j<Ny; j++)
{
for (i=Nx-1; i<Nx; i++)
{
Q[i][j].SetRho( Q[Nx-2][j].GetRho());
Q[i][j].SetRhoU(Q[Nx-2][j].GetRhoU());
Q[i][j].SetRhoV(Q[Nx-2][j].GetRhoV());
Q[i][j].SetRhoE(Q[Nx-2][j].GetRhoE());
}
}
for (i=0; i<Nx; i++)
{
for (j=Ny-1; j<Ny; j++)
{
Q[i][j].SetRho(Q[i][Ny-2].GetRho());
Q[i][j].SetRhoU(Q[i][Ny-2].GetRhoU());
Q[i][j].SetRhoV(Q[i][Ny-2].GetRhoV());
Q[i][j].SetRhoE(Q[i][Ny-2].GetRhoE());
}
}
}
}
for (i=0; i<Nx; i++)
{
for (j=0; j<Ny; j++)
{
double Ma=sqrt((Q[i][j].GetRhoU()*Q[i][j].GetRhoU()+Q[i][j].GetRhoV()*Q[i][j].GetRhoV())/gamma/Q[i][j].GetP()/Q[i][j].GetRho());
fprintf(fp, "%.6lf\t%.6lf\t%.6lf\t%.6lf\t%.6lf\n", \
Q[i][j].SW.x, Q[i][j].SW.y, Q[i][j].GetP(), Ma, Q[i][j].GetRhoE()/Q[i][j].GetRho() );
}
}
fclose(fp);
return 0;
}⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -