亚音速-超音速等熵流动 守恒型 C++代码

#include <iostream>
#include <fstream>
#include <math.h>
using namespace std;

const int N=31;		//节点数
double dx=0.1;		//网格间距
double x[N];
double dt[N],dt_min;
int timestep=1400;
double A[N];
double r=1.4;		//比热比
double C=0.5;		//柯朗数

double rho[N],T[N],u[N];
double U1[N],U2[N],U3[N],F1[N],F2[N],F3[N],J2[N];
double dU1_dt[N],dU2_dt[N],dU3_dt[N];
double U1_1[N],U2_1[N],U3_1[N];
double rho_1[N],T_1[N];
double F1_1[N],F2_1[N],F3_1[N],J2_1[N];
double dU1_dt_1[N],dU2_dt_1[N],dU3_dt_1[N];
double dU1_dt_av[N],dU2_dt_av[N],dU3_dt_av[N];
double U1_2[N],U2_2[N],U3_2[N];
double rho2[N],T2[N],u2[N];
double p[N];
double p_e=0.93;

int main()
{
	int i,j;
	ofstream outfile;
	ofstream mycout;
	mycout.open("midpoint.txt");
	mycout<<"x\t"<<1.5<<"\n";
	mycout<<"time\t"<<"rho2\t"<<"u2\t"<<"p\t"<<"T2\t"<<"Ma\n";
//----------------初始化
	outfile.open("0.txt");
	outfile<<"x\t"<<"A\t"<<"rho\t"<<"u\t"<<"V\t"<<"T\t"<<"U1\t"<<"U2\t"<<"U3\n";
	for(i=0;i<N;i++)
	{
		x[i]=i*dx;
		A[i]=1.0 + 2.2*(x[i]-1.5)*(x[i]-1.5);
		if(x[i]<=0.5)
		{
			rho[i]=1.0;
			T[i]=1.0;
		}
		else if(x[i]<=1.5)
		{
			rho[i]=1.0 - 0.366*(x[i]-0.5);
			T[i]=1.0 - 0.167*(x[i]-0.5);
		}
		else
		{
			rho[i]=0.634 - 0.3879*(x[i]-1.5);
			T[i]=0.833 - 0.3507*(x[i]-1.5);
		}
		u[i] = 0.59/(rho[i]*A[i]);

		U1[i] = rho[i]*A[i];
		U2[i] = u[i]*rho[i]*A[i];
		U3[i] = rho[i]*A[i]*( T[i]/(r-1) +r/2*u[i]*u[i] );
		outfile<<x[i]<<"\t"<<A[i]<<"\t"<<rho[i]<<"\t"<<u[i]<<"\t"<<T[i]<<"\t"<<U1[i]<<"\t"<<U2[i]<<"\t"<<U3[i]<<endl;
	}
	outfile.close();

	for(j=0;j<timestep;j++)
	{
	//----------------计算F、J
		outfile.open("1.txt");
		outfile<<"x\t"<<"F1\t"<<"F2\t"<<"F3\t"<<"J2\n";
		for(i=0;i<N;i++)
		{
			F1[i] = U2[i];
			F2[i] = pow(U2[i],2)/U1[i] + (r-1)/r * (U3[i] - r/2*pow(U2[i],2)/U1[i]);
			F3[i] = r*U2[i]*U3[i]/U1[i] - r*(r-1)/2 * pow(U2[i],3)/pow(U1[i],2);
			J2[i] = 1/r * rho[i]*T[i] *(A[i+1]-A[i])/dx; 
		}
		i=15;
		outfile<<x[i]<<"\t"<<F1[i]<<"\t"<<F2[i]<<"\t"<<F3[i]<<"\t"<<J2[i]<<endl;
		outfile.close();
	//----------------预估步计算
		outfile.open("2.txt");
		outfile<<"x\t"<<"dU1_dt\t"<<"dU2_dt\t"<<"dU3_dt\t"<<"\n";
		for(i=0;i<N-1;i++)
		{
			dU1_dt[i] = -(F1[i+1]-F1[i])/dx;
			dU2_dt[i] = -(F2[i+1]-F2[i])/dx + J2[i];
			dU3_dt[i] = -(F3[i+1]-F3[i])/dx;
		}
		i=15;
		outfile<<x[i]<<"\t"<<dU1_dt[i]<<"\t"<<dU2_dt[i]<<"\t"<<dU3_dt[i]<<endl;
		outfile.close();
	//----------------计算限制时间步长
		dt_min=1.0;
		for(i=0;i<N;i++)
		{
			dt[i]=C*dx/(u[i]+sqrt(T[i]));
			if(dt_min>dt[i])
			{dt_min=dt[i];}
		}
	//----------------计算预估值
		outfile.open("3.txt");
		outfile<<"dt=\t"<<dt_min<<endl;
		outfile<<"x\t"<<"U1_1\t"<<"U2_1\t"<<"U3_1\t"<<"\n";
		for(i=0;i<N-1;i++)
		{
			U1_1[i] = U1[i] + dU1_dt[i]*dt_min;
			U2_1[i] = U2[i] + dU2_dt[i]*dt_min;
			U3_1[i] = U3[i] + dU3_dt[i]*dt_min;
		}	
		i=15;
		outfile<<x[i]<<"\t"<<U1_1[i]<<"\t"<<U2_1[i]<<"\t"<<U3_1[i]<<endl;
		outfile.close();
	//----------------校正步计算
		outfile.open("4.txt");
		outfile<<"x\t"<<"rho_1\t"<<"T_1\t"<<"F1_1\t"<<"F2_1\t"<<"F2_1\t"<<"dU1_dt_1\t"<<"dU2_dt_1\t"<<"dU3_dt_1\t"<<"\n";
		for(i=0;i<N-1;i++)
		{
			F1_1[i] = U2_1[i];
			F2_1[i] = pow(U2_1[i],2)/U1_1[i] + (r-1)/r * (U3_1[i] - r/2*pow(U2_1[i],2)/U1_1[i]);
			F3_1[i] = r*U2_1[i]*U3_1[i]/U1_1[i] - r*(r-1)/2 * pow(U2_1[i],3)/pow(U1_1[i],2);
			rho_1[i] = U1_1[i]/A[i];
			T_1[i] = (r-1)*(U3_1[i]/U1_1[i] - r/2*pow(U2_1[i]/U1_1[i],2));
			J2_1[i] = 1/r * rho_1[i]*T_1[i] *(A[i]-A[i-1])/dx; 
						
			dU1_dt_1[i] = -(F1_1[i] - F1_1[i-1])/dx;
			dU2_dt_1[i] = -(F2_1[i] - F2_1[i-1])/dx + J2_1[i];
			dU3_dt_1[i] = -(F3_1[i] - F3_1[i-1])/dx;
			outfile<<x[i]<<"\t"<<rho_1[i]<<"\t"<<T_1[i]<<"\t"<<F1_1[i]<<"\t"<<F2_1[i]<<"\t"<<F3_1[i]<<"\t"<<dU1_dt_1[i]<<"\t"<<dU2_dt_1[i]<<"\t"<<dU3_dt_1[i]<<endl;
		}	
		outfile.close();
	//----------------平均时间导数计算
		outfile.open("5.txt");
		outfile<<"x\t"<<"dU1_dt_av\t"<<"dU2_dt_av\t"<<"dU3_dt_av\t"<<"\n";
		for(i=1;i<N-1;i++)
		{
			dU1_dt_av[i] = (dU1_dt[i] + dU1_dt_1[i])/2;
			dU2_dt_av[i] = (dU2_dt[i] + dU2_dt_1[i])/2;
			dU3_dt_av[i] = (dU3_dt[i] + dU3_dt_1[i])/2;
		}	
		i=15;
		outfile<<x[i]<<"\t"<<dU1_dt_av[i]<<"\t"<<dU2_dt_av[i]<<"\t"<<dU3_dt_av[i]<<endl;
		outfile.close();
	//----------------计算校正值
		outfile.open("6.txt");
		outfile<<"x\t"<<"U1_2\t"<<"U2_2\t"<<"U3_2\t"<<"\n";
		for(i=1;i<N-1;i++)
		{
			U1_2[i] = U1[i] + dU1_dt_av[i]*dt_min;
			U2_2[i] = U2[i] + dU2_dt_av[i]*dt_min;
			U3_2[i] = U3[i] + dU3_dt_av[i]*dt_min;
		}	
		i=15;
		outfile<<x[i]<<"\t"<<U1_2[i]<<"\t"<<U2_2[i]<<"\t"<<U3_2[i]<<endl;
		outfile.close();
	//----------------计算原始变量
		outfile.open("result.txt");
		outfile<<"x\t"<<"rho2\t"<<"u2\t"<<"T2\t"<<"p\t"<<"Ma\t"<<"U1_2\t"<<"U2_2\t"<<"U3_2\t"<<"\n";
		for(i=1;i<N-1;i++)
		{
			rho2[i] = U1_2[i]/A[i];
			u2[i] = U2_2[i]/U1_2[i];
			T2[i] = (r-1)*(U3_2[i]/U1_2[i] - r/2*pow(U2_2[i]/U1_2[i],2));
			p[i] = rho2[i]*T2[i];
			
		}		
	//----------------边界点值计算-线性外插值
		U1_2[0]=2*U1_2[1]-U1_2[2];
		U2_2[0]=2*U2_2[1]-U2_2[2];
		U3_2[0]=2*U3_2[1]-U3_2[2];
		u2[0]=2*u2[1]-u2[2];
		rho2[0]=1.0;
		T2[0]=1.0;
		p[0] = rho2[0]*T2[0];

		U1_2[N-1]=2*U1_2[N-2]-U1_2[N-3];
		U2_2[N-1]=2*U2_2[N-2]-U2_2[N-3];
		U3_2[N-1]=2*U3_2[N-2]-U3_2[N-3];
		u2[N-1]=2*u2[N-2]-u2[N-3];
		rho2[N-1]=2*rho2[N-2]-rho2[N-3];
		T2[N-1]=2*T2[N-2]-T2[N-3];
		p[N-1] = rho2[N-1]*T2[N-1];
	
		for(i=0;i<N;i++)
		{
			outfile<<x[i]<<"\t"<<rho2[i]<<"\t"<<u2[i]<<"\t"<<T2[i]<<"\t"<<p[i]<<"\t"<<u2[i]/sqrt(T2[i])<<"\t"<<U1_2[i]<<"\t"<<U2_2[i]<<"\t"<<U3_2[i]<<endl;
		
			u[i]=u2[i];
			rho[i]=rho2[i];
			T[i]=T2[i];
		}	
	//-------------------输出中点结果
		mycout<<j<<"\t"<<rho2[15]<<"\t"<<u2[15]<<"\t"<<p[15]<<"\t"<<T2[15]<<"\t"<<u2[15]/sqrt(T2[15])<<endl;
	}
	mycout.close();
}