拟一维超音速流动 激波捕捉法 C++代码

拟一维超音速流动 激波捕捉法 C++代码

#include <iostream>
#include <fstream>
#include <math.h>
using namespace std;const int N=61;		//节点数
double dx=3.0/(N-1);		//网格间距
double x[N];
double dt[N],dt_min;
int timestep=1400;
double A[N];
double r=1.4;		//比热比
double C=0.5;		//柯朗数
double C_x=0.2;		//人工粘性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.6784;double a,b;
double u_1[N],p_1[N];
double S1[N],S2[N],S3[N],S1_1[N],S2_1[N],S3_1[N];int main()
{int i,j;ofstream outfile;ofstream mycout;mycout.open(&#");mycout<<"xt"<<1.5<<"n";mycout<<"timet"<<"rho2t"<<"u2t"<<"pt"<<"T2t"<<"Man";//----------------初始化outfile.open(&#");outfile<<"xt"<<"At"<<"rhot"<<"ut"<<"Vt"<<"Tt"<<"U1t"<<"U2t"<<"U3n";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 if(x[i]<=2.1){rho[i]=0.634 - 0.702*(x[i]-1.5);T[i]=0.833 - 0.4908*(x[i]-1.5);}else{rho[i]=0.5892 - 0.10228*(x[i]-2.1);T[i]=0.93968 - 0.0622*(x[i]-2.1);}u[i] = 0.59/(rho[i]*A[i]);p[i] = rho[i]*T[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、Joutfile.open(&#");outfile<<"xt"<<"F1t"<<"F2t"<<"F3t"<<"J2n";for(i=0;i<N-1;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(&#");outfile<<"xt"<<"dU1_dtt"<<"dU2_dtt"<<"dU3_dtt"<<"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(&#");outfile<<"dt=t"<<dt_min<<endl;outfile<<"xt"<<"U1_1t"<<"U2_1t"<<"U3_1t"<<"n";for(i=0;i<N;i++){a = C_x * fabs(p[i+1]-2*p[i]+p[i-1])/(p[i+1]+2*p[i]+p[i-1]);S1[i] = a * (U1[i+1]-2*U1[i]+U1[i-1]);S2[i] = a * (U2[i+1]-2*U2[i]+U2[i-1]);S3[i] = a * (U3[i+1]-2*U3[i]+U3[i-1]);U1_1[i] = U1[i] + dU1_dt[i]*dt_min + S1[i];U2_1[i] = U2[i] + dU2_dt[i]*dt_min + S2[i];U3_1[i] = U3[i] + dU3_dt[i]*dt_min + S3[i];}	i=15;outfile<<x[i]<<"t"<<U1_1[i]<<"t"<<U2_1[i]<<"t"<<U3_1[i]<<endl;outfile.close();//----------------校正步计算outfile.open(&#");outfile<<"xt"<<"rho_1t"<<"T_1t"<<"F1_1t"<<"F2_1t"<<"F2_1t"<<"dU1_dt_1t"<<"dU2_dt_1t"<<"dU3_dt_1t"<<"n";for(i=0;i<N;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];u_1[i] = U2_1[i]/U1_1[i];T_1[i] = (r-1)*(U3_1[i]/U1_1[i] - r/2*pow(U2_1[i]/U1_1[i],2));p_1[i] = rho_1[i]*T_1[i];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(&#");outfile<<"xt"<<"dU1_dt_avt"<<"dU2_dt_avt"<<"dU3_dt_avt"<<"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(&#");outfile<<"xt"<<"U1_2t"<<"U2_2t"<<"U3_2t"<<"n";for(i=1;i<N-1;i++){b = C_x * fabs(p_1[i+1]-2*p_1[i]+p_1[i-1])/(p_1[i+1]+2*p_1[i]+p_1[i-1]);S1_1[i] = b * (U1_1[i+1]-2*U1_1[i]+U1_1[i-1]);S2_1[i] = b * (U2_1[i+1]-2*U2_1[i]+U2_1[i-1]);S3_1[i] = b * (U3_1[i+1]-2*U3_1[i]+U3_1[i-1]);U1_2[i] = U1[i] + dU1_dt_av[i]*dt_min +S1_1[i];U2_2[i] = U2[i] + dU2_dt_av[i]*dt_min +S2_1[i];U3_2[i] = U3[i] + dU3_dt_av[i]*dt_min +S3_1[i];}	i=15;outfile<<x[i]<<"t"<<U1_2[i]<<"t"<<U2_2[i]<<"t"<<U3_2[i]<<endl;outfile.close();	//----------------边界点值计算-线性外插值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];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];u2[N-1] = U2_2[N-1] / U1_2[N-3];U3_2[N-1] = p_e*A[N-1]/(r-1) + r/2*U2_2[N-1]*u2[N-1];//----------------计算原始变量outfile.open(&#");outfile<<"xt"<<"rho2t"<<"u2t"<<"T2t"<<"pt"<<"Mat"<<"U1_2t"<<"U2_2t"<<"U3_2t"<<"n";for(i=0;i<N;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];}	rho2[0]=1.0;T2[0]=1.0;p[N-1] = p_e;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();
}