com_fg.h

C语言实现二维无粘欧拉方程的求解

#include "iostream.h"
#include "math.h"
#include "stdio.h"
#include "COM_FG_LRUD.h"

void COM_FG(double U[4][N+3][M+3],double F[4][N+3][M+4],double G[4][N+4][M+3],
			double dt[1],double Area[N+3][M+3],
		    double L_x[N+4][M+3],double L_y[N+4][M+3],double H_x[N+3][M+4],double H_y[N+3][M+4])
{
	int i,j,m,n;
	double D,u,v,E,a;
	double lemda_x_max=0,lemda_y_max=0;
	double LEMDA_x[4][4],Tx[4][4],INV_Tx[4][4],LEMDA_y[4][4],Ty[4][4],INV_Ty[4][4];
    double Ax[4][4],Ax1[4][4],Ay[4][4],Ay1[4][4],A_U[4];
    double F_LAR[4],U_LSR[4],G_DAU[4],U_DSU[4];

    double U_L[4][N+3][M+4],U_R[4][N+3][M+4],U_U[4][N+4][M+3],U_D[4][N+4][M+3];
	double F_L[4][N+3][M+4],F_R[4][N+3][M+4],G_U[4][N+4][M+3],G_D[4][N+4][M+3];
	double rfa,bata,cigma;

	FILE *fp;
	fp=fopen("check.dat","w");
	COM_FG_LRUD(U,U_L,U_R,U_U,U_D,F_L,F_R,G_U,G_D,L_x,L_y,H_x,H_y);

	for(i=2;i<N+1;i++)
		for(j=2;j<M+2;j++)
		{
			D=sqrt(U_R[0][i][j]/U_L[0][i][j]);
			u=(U_L[1][i][j]/U_L[0][i][j]+D*U_R[1][i][j]/U_R[0][i][j])/(1+D);
			v=(U_L[2][i][j]/U_L[0][i][j]+D*U_R[2][i][j]/U_R[0][i][j])/(1+D);
			E=(U_L[3][i][j]/U_L[0][i][j]+D*U_R[3][i][j]/U_R[0][i][j])/(1+D);
			a=sqrt(k*(k-1)*(E-(u*u+v*v)/2));

			rfa=H_x[i][j]/Area[i][j];
			bata=H_y[i][j]/Area[i][j];
			cigma=sqrt(rfa*rfa+bata*bata);

			LEMDA_x[0][0]=fabs(rfa*u+bata*v);
			LEMDA_x[0][1]=LEMDA_x[0][2]=LEMDA_x[0][3]=0.0;
			LEMDA_x[1][1]=fabs(rfa*u+bata*v);
			LEMDA_x[1][0]=LEMDA_x[1][2]=LEMDA_x[1][3]=0.0;
			LEMDA_x[2][2]=fabs(rfa*u+bata*v-a*cigma);
			LEMDA_x[2][0]=LEMDA_x[2][1]=LEMDA_x[2][3]=0.0;
			LEMDA_x[3][3]=fabs(rfa*u+bata*v+a*cigma);
			LEMDA_x[3][0]=LEMDA_x[3][1]=LEMDA_x[3][2]=0.0;
			
			if((fabs(rfa*u+bata*v)+a*cigma)>lemda_x_max)
				lemda_x_max=fabs(rfa*u+bata*v)+a*cigma;

			Tx[0][0]=(1-k)/(a*a);
			Tx[0][1]=0.0;
			Tx[0][2]=-1/(2*a*cigma);
			Tx[0][3]=1/(2*a*cigma);

			Tx[1][0]=(1-k)*u/(a*a);
			Tx[1][1]=-bata/(cigma*cigma);
			Tx[1][2]=rfa/(2*cigma*cigma)-u/(2*a*cigma);
			Tx[1][3]=rfa/(2*cigma*cigma)+u/(2*a*cigma);

			Tx[2][0]=(1-k)*v/(a*a);
			Tx[2][1]=rfa/(cigma*cigma);
			Tx[2][2]=bata/(2*cigma*cigma)-v/(2*a*cigma);
			Tx[2][3]=bata/(2*cigma*cigma)+v/(2*a*cigma);

			Tx[3][0]=0.5*(1-k)*(u*u+v*v)/(a*a);
			Tx[3][1]=(rfa*v-bata*u)/(cigma*cigma);
			Tx[3][2]=((a*rfa*u+a*bata*v)/cigma-(u*u+v*v)/2-(a*a)/(k-1))/(2*a*cigma);
			Tx[3][3]=((a*rfa*u+a*bata*v)/cigma+(u*u+v*v)/2+(a*a)/(k-1))/(2*a*cigma);

			
			INV_Tx[0][0]=(u*u+v*v)/2-(a*a)/(k-1);
			INV_Tx[0][1]=-u;
			INV_Tx[0][2]=-v;
			INV_Tx[0][3]=1.0;

			INV_Tx[1][0]=-rfa*v+bata*u;
			INV_Tx[1][1]=-bata;
			INV_Tx[1][2]=rfa;
			INV_Tx[1][3]=0.0;

			INV_Tx[2][0]=-rfa*u-bata*v+(u*u+v*v)*cigma*(1-k)/(2*a);
			INV_Tx[2][1]=rfa+(k-1)*u*cigma/a;
			INV_Tx[2][2]=bata+(k-1)*v*cigma/a;
			INV_Tx[2][3]=(1-k)*cigma/a;

			INV_Tx[3][0]=-rfa*u-bata*v-(u*u+v*v)*cigma*(1-k)/(2*a);
			INV_Tx[3][1]=rfa-(k-1)*u*cigma/a;
			INV_Tx[3][2]=bata-(k-1)*v*cigma/a;
			INV_Tx[3][3]=(k-1)*cigma/a;

			for(m=0;m<4;m++)
				for(n=0;n<4;n++)
				{
					Ax[m][n]=Tx[m][0]*LEMDA_x[0][n]+Tx[m][1]*LEMDA_x[1][n]
					        +Tx[m][2]*LEMDA_x[2][n]+Tx[m][3]*LEMDA_x[3][n];
				}
			for(m=0;m<4;m++)
				for(n=0;n<4;n++)
				{
					Ax1[m][n]=Ax[m][0]*INV_Tx[0][n]+Ax[m][1]*INV_Tx[1][n]
							 +Ax[m][2]*INV_Tx[2][n]+Ax[m][3]*INV_Tx[3][n];
				}

			for(m=0;m<4;m++)
				{
					F_LAR[m]=F_L[m][i][j]+F_R[m][i][j];
					U_LSR[m]=(U_L[m][i][j]-U_R[m][i][j])*Area[i][j];
				}
			for(m=0;m<4;m++)
				{
					A_U[m]=Ax1[m][0]*U_LSR[0]+Ax1[m][1]*U_LSR[1]
					      +Ax1[m][2]*U_LSR[2]+Ax1[m][3]*U_LSR[3];
					F[m][i][j]=0.5*(F_LAR[m]+A_U[m]);
				}
			
		}
	for(j=2;j<M+1;j++)
	    for(i=2;i<N+2;i++)
		{
			D=sqrt(U_U[0][i][j]/U_D[0][i][j]);
			u=(U_D[1][i][j]/U_D[0][i][j]+D*U_U[1][i][j]/U_U[0][i][j])/(1+D);
			v=(U_D[2][i][j]/U_D[0][i][j]+D*U_U[2][i][j]/U_U[0][i][j])/(1+D);
			E=(U_D[3][i][j]/U_D[0][i][j]+D*U_U[3][i][j]/U_U[0][i][j])/(1+D);

			a=sqrt(k*(k-1)*(E-(u*u+v*v)/2));
			rfa=L_x[i][j]/Area[i][j];
			bata=L_y[i][j]/Area[i][j];
			cigma=sqrt(rfa*rfa+bata*bata);

			LEMDA_y[0][0]=fabs(rfa*u+bata*v);
			LEMDA_y[0][1]=LEMDA_y[0][2]=LEMDA_y[0][3]=0.0;
			LEMDA_y[1][1]=fabs(rfa*u+bata*v);
			LEMDA_y[1][0]=LEMDA_y[1][2]=LEMDA_y[1][3]=0.0;
			LEMDA_y[2][2]=fabs(rfa*u+bata*v-a*cigma);
			LEMDA_y[2][0]=LEMDA_y[2][1]=LEMDA_y[2][3]=0.0;
			LEMDA_y[3][3]=fabs(rfa*u+bata*v+a*cigma);
			LEMDA_y[3][0]=LEMDA_y[3][1]=LEMDA_y[3][2]=0.0;

			if((fabs(rfa*u+bata*v)+a*cigma)>lemda_y_max)
				lemda_y_max=fabs(rfa*u+bata*v)+a*cigma;

			Ty[0][0]=(1-k)/(a*a);
			Ty[0][1]=0.0;
			Ty[0][2]=-1/(2*a*cigma);
			Ty[0][3]=1/(2*a*cigma);

			Ty[1][0]=(1-k)*u/(a*a);
			Ty[1][1]=-bata/(cigma*cigma);
			Ty[1][2]=rfa/(2*cigma*cigma)-u/(2*a*cigma);
			Ty[1][3]=rfa/(2*cigma*cigma)+u/(2*a*cigma);

			Ty[2][0]=(1-k)*v/(a*a);
			Ty[2][1]=rfa/(cigma*cigma);
			Ty[2][2]=bata/(2*cigma*cigma)-v/(2*a*cigma);
			Ty[2][3]=bata/(2*cigma*cigma)+v/(2*a*cigma);

			Ty[3][0]=0.5*(1-k)*(u*u+v*v)/(a*a);
			Ty[3][1]=(rfa*v-bata*u)/(cigma*cigma);
			Ty[3][2]=((a*rfa*u+a*bata*v)/cigma-(u*u+v*v)/2-(a*a)/(k-1))/(2*a*cigma);
			Ty[3][3]=((a*rfa*u+a*bata*v)/cigma+(u*u+v*v)/2+(a*a)/(k-1))/(2*a*cigma);

			
			INV_Ty[0][0]=(u*u+v*v)/2-(a*a)/(k-1);
			INV_Ty[0][1]=-u;
			INV_Ty[0][2]=-v;
			INV_Ty[0][3]=1.0;

			INV_Ty[1][0]=-rfa*v+bata*u;
			INV_Ty[1][1]=-bata;
			INV_Ty[1][2]=rfa;
			INV_Ty[1][3]=0.0;

			INV_Ty[2][0]=-rfa*u-bata*v+(u*u+v*v)*cigma*(1-k)/(2*a);
			INV_Ty[2][1]=rfa+(k-1)*u*cigma/a;
			INV_Ty[2][2]=bata+(k-1)*v*cigma/a;
			INV_Ty[2][3]=(1-k)*cigma/a;

			INV_Ty[3][0]=-rfa*u-bata*v-(u*u+v*v)*cigma*(1-k)/(2*a);
			INV_Ty[3][1]=rfa-(k-1)*u*cigma/a;
			INV_Ty[3][2]=bata-(k-1)*v*cigma/a;
			INV_Ty[3][3]=(k-1)*cigma/a;

			for(m=0;m<4;m++)
				for(n=0;n<4;n++)
				{
					Ay[m][n]=Ty[m][0]*LEMDA_y[0][n]+Ty[m][1]*LEMDA_y[1][n]
					        +Ty[m][2]*LEMDA_y[2][n]+Ty[m][3]*LEMDA_y[3][n];
				}

			for(m=0;m<4;m++)
				for(n=0;n<4;n++)
				{
					Ay1[m][n]=Ay[m][0]*INV_Ty[0][n]+Ay[m][1]*INV_Ty[1][n]
					         +Ay[m][2]*INV_Ty[2][n]+Ay[m][3]*INV_Ty[3][n];
				}

			for(m=0;m<4;m++)
				{
					G_DAU[m]=G_D[m][i][j]+G_U[m][i][j];
					U_DSU[m]=(U_D[m][i][j]-U_U[m][i][j])*Area[i][j];
				}
			for(m=0;m<4;m++)
				{
					A_U[m]=Ay1[m][0]*U_DSU[0]+Ay1[m][1]*U_DSU[1]
					      +Ay1[m][2]*U_DSU[2]+Ay1[m][3]*U_DSU[3];
					G[m][i][j]=0.5*(G_DAU[m]+A_U[m]);
				}
		}
		if(lemda_x_max>=lemda_y_max)
			dt[0]=cfl/lemda_x_max;
		else 
			dt[0]=cfl/lemda_y_max;

}