csimple2d.c

CSIMPLE2d CFD C源程序

		fprintf(pFileOut,"               computation for axisymmetric situation\n");
		break;
	case 3:
		fprintf(pFileOut,"               computation in polar coordinates\n");
		break;
	default:
		break;
	}
}

/*********************************************************************************
* function name:	fnSetUp2
* description:		solve the eqs
* input param:		-none-
* out param:		-none-
*********************************************************************************/
void fnSetUp2()
{
	int i,j,n;
	double rel,GM,GMM,VOL,APT,AREA,SXT,SXB,ARHO;
	double FL,FLM,FLP;

	//COEFFICIENTS FOR THE U EQUATION
	nF=0;
    if(bSolve[nF])
	{
		IST=2;
		JST=1;

		fnGamSor();

		rel=1.0-Relax[nF];
		for(i=2;i<=L2;i++)
		{
			FL=XCVI[i]*V[i][1]*Rho[i][0];
			FLM=XCVIP[i-1]*V[i-1][1]*Rho[i-1][0];
			Flow=R[0]*(FL+FLM);
			Diff=R[0]*(XCVI[i]*Gam[i][0]+XCVIP[i-1]*Gam[i-1][0])/YDif[1];
			
			fnDiFlow();

			AJM[i][1]=ACof+max(0,Flow);
		}

		for(j=1;j<=M2;j++)
		{
			Flow=ARX[j]*U[1][j]*Rho[0][j];
			Diff=ARX[j]*Gam[0][j]/(XCV[1]*SX[j]);

			fnDiFlow();

			AIM[2][j]=ACof+max(0,Flow);

			for(i=2;i<=L2;i++)
			{
				if(i!=L2)
				{
					FL=U[i][j]*(FX[i]*Rho[i][j]+FXM[i]*Rho[i-1][j]);
					FLP=U[i+1][j]*(FX[i+1]*Rho[i+1][j]+FXM[i+1]*Rho[i][j]);
					Flow=ARX[j]*(FL+FLP)/2;
					Diff=ARX[j]*Gam[i][j]/(XCV[i]*SX[j]);
				}
				else
				{
					Flow=ARX[j]*U[L1][j]*Rho[L1][j];
					Diff=ARX[j]*Gam[L1][j]/(XCV[L2]*SX[j]);
				}

				fnDiFlow();

				AIM[i+1][j]=ACof+max(0,Flow);
				AIP[i][j]=AIM[i+1][j]-Flow;
				if(j!=M2)
				{
					FL=XCVI[i]*V[i][j+1]*(FY[j+1]*Rho[i][j+1]+FYM[j+1]*Rho[i][j]);
					FLM=XCVIP[i-1]*V[i-1][j+1]*(FY[j+1]*Rho[i-1][j+1]+FYM[j+1]*Rho[i-1][j]);
					GM=Gam[i][j]*Gam[i][j+1]/(YCV[j]*Gam[i][j+1]
						+YCV[j+1]*Gam[i][j]+1e-30)*XCVI[i];
					GMM=Gam[i-1][j]*Gam[i-1][j+1]/(YCV[j]*Gam[i-1][j+1]
						+YCV[j+1]*Gam[i-1][j]+1e-30)*XCVIP[i-1];
					Diff=RMN[j+1]*2*(GM+GMM);
				}
				else
				{
					FL=XCVI[i]*V[i][M1]*Rho[i][M1];
					FLM=XCVIP[i-1]*V[i-1][M1]*Rho[i-1][M1];
					Diff=R[M1]*(XCVI[i]*Gam[i][M1]+XCVIP[i-1]*Gam[i-1][M1])/YDif[M1];
				}
				Flow=RMN[j+1]*(FL+FLM);

				fnDiFlow();

				AJM[i][j+1]=ACof+max(0,Flow);
				AJP[i][j]=AJM[i][j+1]-Flow;

				VOL=YCVR[j]*XCVS[i];
				APT=(Rho[i][j]*XCVI[i]+Rho[i-1][j]*XCVIP[i-1])/(XCVS[i]*Dt);
				AP[i][j]-=APT;
				Con[i][j]+=APT*U[i][j];
				AP[i][j]=(-AP[i][j]*VOL+AIP[i][j]+AIM[i][j]+AJP[i][j]+AJM[i][j])/Relax[nF];
				Con[i][j]=Con[i][j]*VOL+rel*AP[i][j]*U[i][j];
				DU[i][j]=VOL/(XDif[i]*SX[j]);
				Con[i][j]+=DU[i][j]*(P[i-1][j]-P[i][j]);
				DU[i][j]/=AP[i][j];
			}
		}
		fnSolve();
	}
//COEFFICIENTS FOR THE V EQUATION
	nF=1;
    if(bSolve[nF])
	{
		IST=1;
		JST=2;

		fnGamSor();

		rel=1-Relax[nF];
		for(i=1;i<=L2;i++)
		{
			AREA=R[0]*XCV[i];
			Flow=AREA*V[i][1]*Rho[i][0];
			Diff=AREA*Gam[i][0]/YCV[1];

			fnDiFlow();

			AJM[i][2]=ACof+max(0,Flow);
		}

		for(j=2;j<=M2;j++)
		{
			FL=ARXJ[j]*U[1][j]*Rho[0][j];
			FLM=ARXJP[j-1]*U[1][j-1]*Rho[0][j-1];
			Flow=FL+FLM;
			Diff=(ARXJ[j]*Gam[0][j]+ARXJP[j-1]*Gam[0][j-1])/(XDif[1]*SXMN[j]);

			fnDiFlow();

			AIM[1][j]=ACof+max(0,Flow);

			for(i=1;i<=L2;i++)
			{
				if(i!=L2)
				{
					FL=ARXJ[j]*U[i+1][j]*(FX[i+1]*Rho[i+1][j]+FXM[i+1]*Rho[i][j]);
					FLM=ARXJP[j-1]*U[i+1][j-1]*(FX[i+1]*Rho[i+1][j-1]+FXM[i+1]*Rho[i][j-1]);
					GM=Gam[i][j]*Gam[i+1][j]/(XCV[i]*Gam[i+1][j]
						+XCV[i+1]*Gam[i][j]+1e-30)*ARXJ[j];
					GMM=Gam[i][j-1]*Gam[i+1][j-1]/(XCV[i]*Gam[i+1][j-1]
						+XCV[i+1]*Gam[i][j-1]+1.0e-30)*ARXJP[j-1];
					Diff=2*(GM+GMM)/SXMN[j];
				}
				else
				{
					FL=ARXJ[j]*U[L1][j]*Rho[L1][j];
					FLM=ARXJP[j-1]*U[L1][j-1]*Rho[L1][j-1];
					Diff=(ARXJ[j]*Gam[L1][j]+ARXJP[j-1]*Gam[L1][j-1])/(XDif[L1]*SXMN[j]);
				}
				Flow=FL+FLM;

				fnDiFlow();

				AIM[i+1][j]=ACof+max(0,Flow);
				AIP[i][j]=AIM[i+1][j]-Flow;
				if(j!=M2)
				{
					AREA=R[j]*XCV[i];
					FL=V[i][j]*(FY[j]*Rho[i][j]+FYM[j]*Rho[i][j-1])*RMN[j];
					FLP=V[i][j+1]*(FY[j+1]*Rho[i][j+1]+FYM[j+1]*Rho[i][j])*RMN[j+1];
					Flow=(FV[j]*FL+FVP[j]*FLP)*XCV[i];
					Diff=AREA*Gam[i][j]/YCV[j];
				}
				else
				{
					AREA=R[M1]*XCV[i];
					Flow=AREA*V[i][M1]*Rho[i][M1];
					Diff=AREA*Gam[i][M1]/YCV[M2];
				}

				fnDiFlow();

				AJM[i][j+1]=ACof+max(0,Flow);
				AJP[i][j]=AJM[i][j+1]-Flow;
				VOL=YCVRS[j]*XCV[i];
				SXT=SX[j];
				if(j==M2)
					SXT=SX[M1];
				SXB=SX[j-1];
				if(j==2)
					SXB=SX[0];
				APT=(ARXJ[j]*Rho[i][j]*(SXT+SXMN[j])/2
					+ARXJP[j-1]*Rho[i][j-1]*(SXB+SXMN[j])/2)/(YCVRS[j]*Dt);
				AP[i][j]-=APT;
				Con[i][j]+=APT*V[i][j];
				AP[i][j]=(-AP[i][j]*VOL+AIP[i][j]+AIM[i][j]+AJP[i][j]+AJM[i][j])/Relax[nF];
				Con[i][j]=Con[i][j]*VOL+rel*AP[i][j]*V[i][j];
				DV[i][j]=VOL/YDif[j];
				Con[i][j]+=DV[i][j]*(P[i][j-1]-P[i][j]);
				DV[i][j]/=AP[i][j];
			}
		}
		fnSolve();
	}
//COEFIICIENTS FOR THE PRESSURE CORRECTION EQUATION
	nF=2;
    if(bSolve[nF])
	{
		IST=1;
		JST=1;

		fnGamSor();

		SMax=0;
		SSum=0;

		for(i=1;i<=L2;i++)
			for(j=1;j<=M2;j++)
				Con[i][j]*=YCVR[j]*XCV[i];
		for(i=1;i<=L2;i++)
		{
			Con[i][1]+=R[0]*XCV[i]*Rho[i][0]*V[i][1];
			AJM[i][1]=0;
		}
		for(j=1;j<=M2;j++)
		{
			ARHO=ARX[j]*Rho[0][j];
			Con[1][j]+=ARHO*U[1][j];
			AIM[1][j]=0;
			for(i=1;i<=L2;i++)
			{
				if(i!=L2)
				{
					ARHO=ARX[j]*(FX[i+1]*Rho[i+1][j]+FXM[i+1]*Rho[i][j]);
					Flow=ARHO*U[i+1][j];
					Con[i][j]-=Flow;
					Con[i+1][j]+=Flow;
					AIP[i][j]=ARHO*DU[i+1][j];
					AIM[i+1][j]=AIP[i][j];
				}
				else
				{
					ARHO=ARX[j]*Rho[L1][j];
					Con[i][j]-=ARHO*U[L1][j];
					AIP[i][j]=0;
				}

				if(j!=M2)
				{
					ARHO=RMN[j+1]*XCV[i]*(FY[j+1]*Rho[i][j+1]+FYM[j+1]*Rho[i][j]);
					Flow=ARHO*V[i][j+1];
					Con[i][j]-=Flow;
					Con[i][j+1]+=Flow;
					AJP[i][j]=ARHO*DV[i][j+1];
					AJM[i][j+1]=AJP[i][j];
				}
				else 
				{
					ARHO=RMN[M1]*XCV[i]*Rho[i][M1];
					Con[i][j]-=ARHO*V[i][M1];
					AJP[i][j]=0;
				}
				AP[i][j]=AIP[i][j]+AIM[i][j]+AJP[i][j]+AJM[i][j];
				PC[i][j]=0;
				SMax=max(SMax,fabs(Con[i][j]));
				SSum+=Con[i][j];
			}
		}
		fnSolve();

	//COMEE HERE TO CORRECT THE PRESSURE AND VELOCITIES
		for(i=1;i<=L2;i++)
		{
			for(j=1;j<=M2;j++)
			{
				P[i][j]+=PC[i][j]*Relax[nP];
				if(i!=1)
					U[i][j]+=DU[i][j]*(PC[i-1][j]-PC[i][j]);
				if(j!=1)
					V[i][j]+=DV[i][j]*(PC[i][j-1]-PC[i][j]);
			}
		}
	}
//COEFFICIENTS FOR OTHER EQUATIONS
	IST=1;
    JST=1;
    
	for(n=3;n<=nFMax;n++)
	{
		nF=n;
		if(bSolve[nF])
		{
			fnGamSor();

			rel=1-Relax[nF];
			for(i=1;i<=L2;i++)
			{
				AREA=R[0]*XCV[i];
				Flow=AREA*V[i][1]*Rho[i][0];
				Diff=AREA*Gam[i][0]/YDif[1];

				fnDiFlow();

				AJM[i][1]=ACof+max(0,Flow);
			}
			for(j=1;j<=M2;j++)
			{
				Flow=ARX[j]*U[1][j]*Rho[0][j];
				Diff=ARX[j]*Gam[0][j]/(XDif[1]*SX[j]);

				fnDiFlow();

				AIM[1][j]=ACof+max(0,Flow);
				for(i=1;i<=L2;i++)
				{
					if(i!=L2)
					{
						Flow=ARX[j]*U[i+1][j]*(FX[i+1]*Rho[i+1][j]+FXM[i+1]*Rho[i][j]);
						Diff=ARX[j]*2*Gam[i][j]*Gam[i+1][j]
							/((XCV[i]*Gam[i+1][j]+XCV[i+1]*Gam[i][j]+1e-30)*SX[j]);
					}
					else
					{
						Flow=ARX[j]*U[L1][j]*Rho[L1][j];
						Diff=ARX[j]*Gam[L1][j]/(XDif[L1]*SX[j]);
					}

					fnDiFlow();

					AIM[i+1][j]=ACof+max(0,Flow);
					AIP[i][j]=AIM[i+1][j]-Flow;
					AREA=RMN[j+1]*XCV[i];
					if(j!=M2)
					{
						Flow=AREA*V[i][j+1]*(FY[j+1]*Rho[i][j+1]+FYM[j+1]*Rho[i][j]);
						Diff=AREA*2*Gam[i][j]*Gam[i][j+1]/(YCV[j]*Gam[i][j+1]+
							YCV[j+1]*Gam[i][j]+1e-30);
					}
					else 
					{
						Flow=AREA*V[i][M1]*Rho[i][M1];
						Diff=AREA*Gam[i][M1]/YDif[M1];
					}

					fnDiFlow();

					AJM[i][j+1]=ACof+max(0,Flow);
					AJP[i][j]=AJM[i][j+1]-Flow;
					VOL=YCVR[j]*XCV[i];
					APT=Rho[i][j]/Dt;
					AP[i][j]-=APT;
					Con[i][j]+=APT*F[nF][i][j];
					AP[i][j]=(-AP[i][j]*VOL+AIP[i][j]+AIM[i][j]
						+AJP[i][j]+AJM[i][j])/Relax[nF];
					Con[i][j]=Con[i][j]*VOL+rel*AP[i][j]*F[nF][i][j];
				}
			}
			fnSolve();
		}
	}

    Time+=Dt;
    iIter++;
    if(iIter>=iLast)
		bStop=TRUE;
}

/*********************************************************************************
* function name:	fnUGrid
* description:		devide the grid uniformly
* input param:		-none-
* out param:		-none-
*********************************************************************************/
void fnUGrid()
{
	int i,j;
	double dx,dy;

	XU[1]=0;
    dx=XL/(L1-1);
    for(i=2;i<=L1;i++)
		XU[i]=XU[i-1]+dx;

    YV[1]=0;
    dy=YL/(M1-1);
    for(j=2;j<=M1;j++)
		YV[j]=YV[j-1]+dy;
}

/*********************************************************************************
* function name:	fnPrint
* description:		print the result
* input param:		-none-
* out param:		-none-
*********************************************************************************/
void fnPrint()
{
	int i,j,n,m;
	int IFST,JFST;
	double RHOM,PREF;
	if(bPrint[2])
	{
		//CALCULATE THE STREAM FUNTION
		F[2][1][1]=0;
		for(i=1;i<=L1;i++)
		{
			if(i!=1)
				F[2][i][1]=F[2][i-1][1]-Rho[i-1][0]*V[i-1][1]*R[0]*XCV[i-1];
			for(j=2;j<=M1;j++)
			{
				RHOM=FX[i]*Rho[i][j-1]+FXM[i]*Rho[i-1][j-1];
				F[2][i][j]=F[2][i][j-1]+RHOM*U[i][j-1]*ARX[j-1];
			}
		}
	} 
	
	if(bPrint[nP])
	{
		//CONSTRUCT BOUNDARY PRESSURES BY EXTRAPOLATION
		for(j=1;j<=M2;j++)
		{
			P[0][j]=(P[1][j]*XCVS[2]-P[2][j]*XDif[1])/XDif[2];
			P[L1][j]=(P[L2][j]*XCVS[L2]-P[L3][j]*XDif[L1])/XDif[L2];
		}
		for(i=1;i<=L2;i++)
		{
			P[i][0]=(P[i][1]*YCVS[2]-P[i][2]*YDif[1])/YDif[2];
			P[i][M1]=(P[i][M2]*YCVS[M2]-P[i][M3]*YDif[M1])/YDif[M2];
		}
		P[0][0]=P[1][0]+P[0][1]-P[1][1];
		P[L1][0]=P[L2][0]+P[L1][1]-P[L2][1];
		P[0][M1]=P[1][M1]+P[0][M2]-P[1][M2];
		P[L1][M1]=P[L2][M1]+P[L1][M2]-P[L2][M2];
		PREF=P[IPRef][JPRef];
		for(j=0;j<=M1;j++)
			for(i=0;i<=L1;i++)
				P[i][j]-=PREF;
	}
	//输出坐标系
	fprintf(pFileOut,"\n");
	for(n=0;n<=L1/7;n++)
	{
		fprintf(pFileOut,"\nI=  ");
		for(i=0;i<=min(6,L1-7*n);i++)
			fprintf(pFileOut,"   %4d  ",i+7*n);
		if(nMode==3)
			fprintf(pFileOut,"\nTH=  ");
		else
			fprintf(pFileOut,"\nX=  ");
		for(i=0;i<=min(6,L1-7*n);i++)
			fprintf(pFileOut,"%9.2e ",X[i+7*n]);
	}

	fprintf(pFileOut,"\n");
	for(n=0;n<=M1/7;n++)
	{
		fprintf(pFileOut,"\nJ=  ");
		for(j=0;j<=min(6,M1-7*n);j++)
			fprintf(pFileOut,"   %4d  ",j+7*n);
		fprintf(pFileOut,"\nY=  ");
		for(j=0;j<=min(6,M1-7*n);j++)
			fprintf(pFileOut,"%9.2e ",Y[j+7*n]);
	}

	for(n=0;n<=nGam;n++)
	{
		if(bPrint[n])
		{
			fprintf(pFileOut,"\n\n*************************%s*************************",szTitle[n]);
			IFST=0;
			JFST=0;
			if(n==0||n==2)
				IFST=1;
			if(n==1||n==2)
				JFST=1;
			for(m=0;m<=(L1-IFST)/7;m++)
			{
				fprintf(pFileOut,"\nI=");
				for(i=0;i<=min(6,(L1-IFST)-7*m);i++)
					fprintf(pFileOut,"   %7d ",i+7*m+IFST);
				for(j=M1;j>=JFST;j--)
				{
					if(j==M1)
						fprintf(pFileOut,"\nJ=%3d ",j);
					else
						fprintf(pFileOut,"\n  %3d ",j);
					for(i=0;i<=min(6,(L1-IFST)-7*m);i++)
						fprintf(pFileOut,"%10.2e ",F[n][i+7*m+IFST][j]);
				}
			}
		}
	}
}

/***************************************************************
function name:	power
description:	calc the value of a^b
intput param:	a,b 
output param:	a^b
****************************************************************/
double power(double a,int b)
{
	int i;
	double d;
	d=1;
	for(i=0;i<b;i++)
		d*=a;
	return d;
}