全亚声速等熵流动 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=5000;
double A[N];
double r=1.4;		//比热比
double C=0.5;		//柯朗数
double rho[N],T[N],u[N];
double rho_t[N],T_t[N],u_t[N];
double rho1[N],T1[N],u1[N];
double rho_t1[N],T_t1[N],u_t1[N];
double rho_t_av[N],T_t_av[N],u_t_av[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"<<"T\t"<<"u\n";
	for(i=0;i<N;i++)
	{
		x[i]=i*dx;
		if( x[i]<=1.5 )
		{
			A[i]=1.0 + 2.2*(x[i]-1.5)*(x[i]-1.5);
		}
		else
		{
			A[i]=1.0 + 0.2223*(x[i]-1.5)*(x[i]-1.5);
		}
		rho[i]=1.0 - 0.023*x[i];
		T[i]=1.0 - 0.009333*x[i];
		u[i]=0.05 + 0.11*x[i];
		outfile<<x[i]<<"\t"<<A[i]<<"\t"<<rho[i]<<"\t"<<T[i]<<"\t"<<u[i]<<endl;
	}
	outfile.close();

	for(j=0;j<timestep;j++)
	{
	//----------------预估步计算
		outfile.open("1.txt");
		outfile<<"x\t"<<"rho_t\t"<<"u_t\t"<<"T_t\t"<<"\n";
		for(i=0;i<N-1;i++)
		{
			rho_t[i]=( -u[i]*(rho[i+1]-rho[i]) -rho[i]*(u[i+1]-u[i]) -rho[i]*u[i]*(log(A[i+1])-log(A[i])) )/dx;
			u_t[i]=( -u[i]*(u[i+1]-u[i]) -((T[i+1]-T[i]) + T[i]/rho[i]*(rho[i+1]-rho[i]))/r )/dx;
			T_t[i]=( -u[i]*(T[i+1]-T[i]) -(r-1)*T[i]*((u[i+1]-u[i]) + u[i]*(log(A[i+1])-log(A[i]))) )/dx;
			outfile<<x[i]<<"\t"<<rho_t[i]<<"\t"<<u_t[i]<<"\t"<<T_t[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("2.txt");
		outfile<<"dt=\t"<<dt_min<<endl;
		outfile<<"x\t"<<"rho1\t"<<"u1\t"<<"T1\t"<<"\n";
		for(i=0;i<N-1;i++)
		{
			rho1[i]=rho[i]+rho_t[i]*dt_min;
			u1[i]=u[i]+u_t[i]*dt_min;
			T1[i]=T[i]+T_t[i]*dt_min;
			outfile<<x[i]<<"\t"<<rho1[i]<<"\t"<<u1[i]<<"\t"<<T1[i]<<endl;
		}	
		outfile.close();
	//----------------校正步计算
		outfile.open("3.txt");
		outfile<<"x\t"<<"rho_t1\t"<<"u_t1\t"<<"T_t1\t"<<"\n";
		for(i=1;i<N-1;i++)
		{
			rho_t1[i]=( -u1[i]*(rho1[i]-rho1[i-1]) -rho1[i]*(u1[i]-u1[i-1]) -rho1[i]*u1[i]*(log(A[i])-log(A[i-1])) )/dx;
			u_t1[i]=( -u1[i]*(u1[i]-u1[i-1]) -((T1[i]-T1[i-1]) + T1[i]/rho1[i]*(rho1[i]-rho1[i-1]))/r )/dx;
			T_t1[i]=( -u1[i]*(T1[i]-T1[i-1]) -(r-1)*T1[i]*((u1[i]-u1[i-1]) + u1[i]*(log(A[i])-log(A[i-1]))) )/dx;
			outfile<<x[i]<<"\t"<<rho_t1[i]<<"\t"<<u_t1[i]<<"\t"<<T_t1[i]<<endl;
		}	
		outfile.close();
	//----------------平均时间导数计算
		outfile.open("4.txt");
		outfile<<"x\t"<<"rho_t_av\t"<<"u_t_av\t"<<"T_t_av\t"<<"\n";
		for(i=1;i<N-1;i++)
		{
			rho_t_av[i]=(rho_t[i]+rho_t1[i])/2;
			u_t_av[i]=(u_t[i]+u_t1[i])/2;
			T_t_av[i]=(T_t[i]+T_t1[i])/2;
			outfile<<x[i]<<"\t"<<rho_t_av[i]<<"\t"<<u_t_av[i]<<"\t"<<T_t_av[i]<<endl;
		}	
		outfile.close();
	//----------------计算校正值
		outfile.open("5.txt");
		outfile<<"x\t"<<"rho2\t"<<"u2\t"<<"T2\t"<<"\n";
		for(i=1;i<N-1;i++)
		{
			rho2[i]=rho[i]+rho_t_av[i]*dt_min;
			u2[i]=u[i]+u_t_av[i]*dt_min;
			T2[i]=T[i]+T_t_av[i]*dt_min;

			outfile<<x[i]<<"\t"<<rho2[i]<<"\t"<<u2[i]<<"\t"<<T2[i]<<endl;
		}		
		outfile.close();
	//----------------边界点值计算-线性外插值
		outfile.open("7.txt");
		outfile<<"x\t"<<"rho2\t"<<"u2\t"<<"T2\t"<<"\n";
		u2[0]=2*u2[1]-u2[2];
		rho2[0]=1.0;
		T2[0]=1.0;

		u2[N-1]=2*u2[N-2]-u2[N-3];
		rho2[N-1]=2*rho2[N-2]-rho2[N-3];
		T2[N-1]=p_e/rho2[N-1];
	
		i=0;
		outfile<<x[i]<<"\t"<<rho2[i]<<"\t"<<u2[i]<<"\t"<<T2[i]<<endl;
		i=N-1;
		outfile<<x[i]<<"\t"<<rho2[i]<<"\t"<<u2[i]<<"\t"<<T2[i]<<endl;
		outfile.close();
	//----------------计算压力值
		outfile.open("6.txt");
		outfile<<"x\t"<<"p\n";
		for(i=0;i<N-1;i++)
		{
			p[i]=rho2[i]*T2[i];
			outfile<<x[i]<<"\t"<<p[i]<<endl;
		}	
		p[N-1] = p_e;
		outfile.close();

		for(i=0;i<N;i++)
		{
			u[i]=u2[i];
			rho[i]=rho2[i];
			T[i]=T2[i];
		}	
	//-------------------输出结果
		outfile.open("8.txt");
		outfile<<"x\t"<<"rho2\t"<<"u2\t"<<"T2\t"<<"Ma\t"<<"p_e\t"<<"\n";
		for(i=0;i<N;i++)
		{
			outfile<<x[i]<<"\t"<<rho2[i]<<"\t"<<u2[i]<<"\t"<<T2[i]<<"\t"<<u2[i]/sqrt(T2[i])<<"\t"<<p[i]<<endl;
		}		
		outfile.close();

	//-------------------输出中点结果
		mycout<<j<<"\t"<<rho2[15]<<"\t"<<u2[15]<<"\t"<<p[15]<<"\t"<<T2[15]<<"\t"<<u2[15]/sqrt(T2[15])<<endl;
	}
	mycout.close();
}