boundary.h

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

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

void Boundary(double U[4][N+3][M+3],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;
	double kx,ky,kx1,ky1,detak;
	double pb,roub,ub,vb,Eb;
	double pp,roup,up,vp,Ep;
	double pa,roua,ua,va,Ea,aa;
	double pr,rour,ur,vr,Er,ar;
	double pl,roul,ul,vl,El,al;
		 
	for(j=2;j<M+1;j++)
	{
		kx=L_x[2][j]/Area[2][j];
		ky=L_y[2][j]/Area[2][j];
		detak=sqrt(kx*kx+ky*ky);
		kx1=kx/detak;
		ky1=ky/detak;

		rour=U[0][2][j];
		ur=U[1][2][j]/U[0][2][j];
		vr=U[2][2][j]/U[0][2][j];
		Er=U[3][2][j]/U[0][2][j];
		pr=(Er-(ur*ur+vr*vr)/2)*rour*(k-1);
		ar=sqrt(k*pr/rour);

		pb=pr;
		roub=rour;
		ub=ur-kx1*(kx1*ur+ky1*vr);
		vb=vr-ky1*(kx1*ur+ky1*vr);
		Eb=pb/((k-1)*roub)+(ub*ub+vb*vb)/2;

		pp=2*pb-pr;
		roup=2*roub-rour;
		up=2*ub-ur;
		vp=2*vb-vr;
		Ep=pp/((k-1)*roup)+(up*up+vp*vp)/2;

		U[0][0][j]=U[0][1][j]=roup;
		U[1][0][j]=U[1][1][j]=roup*up;
		U[2][0][j]=U[2][1][j]=roup*vp;
		U[3][0][j]=U[3][1][j]=roup*Ep;
	}

	for(j=2;j<M+1;j++)
	{
		kx=L_x[N+1][j]/Area[N+1][j];
		ky=L_y[N+1][j]/Area[N+1][j];
		detak=sqrt(kx*kx+ky*ky);
		kx1=kx/detak;
		ky1=ky/detak;

		rour=U[0][N][j];
		ur=U[1][N][j]/U[0][N][j];
		vr=U[2][N][j]/U[0][N][j];
		Er=U[3][N][j]/U[0][N][j];
		pr=(Er-(ur*ur+vr*vr)/2)*rour*(k-1);
		ar=sqrt(k*pr/rour);

		pb=pr;
		roub=rour;
		ub=ur-kx1*(kx1*ur+ky1*vr);
		vb=vr-ky1*(kx1*ur+ky1*vr);
		Eb=pb/((k-1)*roub)+(ub*ub+vb*vb)/2;

		pp=2*pb-pr;
		roup=2*roub-rour;
		up=2*ub-ur;
		vp=2*vb-vr;
		Ep=pp/((k-1)*roup)+(up*up+vp*vp)/2;

		U[0][N+1][j]=U[0][N+2][j]=roup;
		U[1][N+1][j]=U[1][N+2][j]=roup*up;
		U[2][N+1][j]=U[2][N+2][j]=roup*vp;
		U[3][N+1][j]=U[3][N+2][j]=roup*Ep;
	}
//////////////////////////////////////////////////////////////////////////////////////////////////
//下为超声速进口
	for(i=0;i<N+3;i++)
		for(j=0;j<2;j++)
		{
			U[0][i][j]=rou0;
			U[1][i][j]=rou0*u0;
			U[2][i][j]=rou0*v0;
			U[3][i][j]=p0/(k-1)+0.5*rou0*(u0*u0+v0*v0);
			//cout<<U[3][i][j]<<endl;
		}
//////////////////////////////////////////////////////////////////////////////////////////////////
//下为亚声速进口
/*	for(i=0;i<N+3;i++)
		{
			kx=H_x[i][2]/Area[i][2];
			ky=H_y[i][2]/Area[i][2];
			detak=sqrt(kx*kx+ky*ky);
			kx1=kx/detak;
			ky1=ky/detak;

			roul=U[0][i][2];
			ul=U[1][i][2]/U[0][i][2];
			vl=U[2][i][2]/U[0][i][2];
			El=U[3][i][2]/U[0][i][2];
			pl=(El-(ul*ul+vl*vl)/2)*roul*(k-1);
			al=sqrt(k*pl/roul);

			roua=rou0;
			pa=p0;
			ua=u0;
			va=v0;
			
			pb=0.5*(pa+pl+roul*al*(kx1*(ua-ul)+ky1*(va-vl)));
			roub=roua+(pb-pa)/(al*al);
			ub=ua+kx1*(pa-pb)/(roul*al);
			vb=va+ky1*(pa-pb)/(roul*al);
			Eb=pb/((k-1)*roub)+(ub*ub+vb*vb)/2;

			pp=2*pb-pl;
	    	roup=2*roub-roul;
			up=2*ub-ul;
			vp=2*vb-vl;
			Ep=pp/((k-1)*roup)+(up*up+vp*vp)/2;

			U[0][i][0]=U[0][i][1]=roup;
			U[1][i][0]=U[1][i][1]=roup*up;
			U[2][i][0]=U[2][i][1]=roup*vp;
			U[3][i][0]=U[3][i][1]=roup*Ep;
		}*/
//////////////////////////////////////////////////////////////////////////////////////////////////
//下为超声速出口
/*	for(i=0;i<N+3;i++)
		for(j=M+1;j<M+3;j++)
		{
			U[0][i][j]=U[0][i][M];
			U[1][i][j]=U[1][i][M];
			U[2][i][j]=U[2][i][M];
			U[3][i][j]=U[0][i][M]*pout0/(k-1)+(pow(U[1][i][M],2)+pow(U[2][i][M],2))/(2.0*U[0][i][M]);
		}*/
////////////////////////////////////////////////////////////////////////////////////////////////
//下为亚声速出口
	for(i=0;i<N+3;i++)
	{
		kx=H_x[i][M+1]/Area[i][M+1];
		ky=H_y[i][M+1]/Area[i][M+1];
		detak=sqrt(kx*kx+ky*ky);
		kx1=kx/detak;
		ky1=ky/detak;

		roua=U[0][i][M];
		ua=U[1][i][M]/U[0][i][M];
		va=U[2][i][M]/U[0][i][M];
		Ea=U[3][i][M]/U[0][i][M];
		pa=(Ea-(ua*ua+va*va)/2)*roua*(k-1);
		aa=sqrt(k*pa/roua);

		
		pb=pout0;
		roub=roua+(pb-pa)/(aa*aa);
		ub=ua+kx1*(pa-pb)/(roua*aa);
		vb=va+ky1*(pa-pb)/(roua*aa);
		Eb=pb/((k-1)*roub)+(ub*ub+vb*vb)/2;

  		pp=2*pb-pa;
	    roup=2*roub-roua;
		up=2*ub-ua;
		vp=2*vb-va;
		Ep=pp/((k-1)*roup)+(up*up+vp*vp)/2;

		U[0][i][M+1]=U[0][i][M+2]=roup;
		U[1][i][M+1]=U[1][i][M+2]=roup*up;
		U[2][i][M+1]=U[2][i][M+2]=roup*vp;
		U[3][i][M+1]=U[3][i][M+2]=roup*Ep;
	}
}