exam9.c

CSIMPLE2d CFD C源程序

/*
* exam9.c : example 9 for SIMPLE
*/
/*
* 注意C与FORTRAN的差别：
* 1. C中的数组下标是从0开始的,而FORTRAN则是从1开始的
*   因此本程序中的一些量(如L1,L2,L3,M1,M2,M3,NFMAX,nGam,nP,nRho,IPref,JPref)比FORTRAN程序中的小1
*   所有数组的下标从0开始,程序中数组下标中出现的数字比FORTRAN中的小1
* 2. C的数组是按行优先读的,而FORTRAN的则是按列读的
*   因此本程序中的F的定义不同,是F[NFMAX][I][J]
* 3. C用指针来实现FORTRAN的EQUIVALENCE功能。
*   如定义double (*U)[NSIZE],U=F[0],则U[i][j]等同于F[0][i][j];
*/

#include "staincl.h"
#include "extern.h"

/*global variable for user */
	double (*Flox)[NSIZE];
	double (*Fuel)[NSIZE];
	double (*Ox)[NSIZE];
	double (*H)[NSIZE];
	double (*T)[NSIZE];
	double (*Rsw)[NSIZE];
	double Tin;
	double Tw;
	double Ai;
	double Cp;
	double TRef;
	double Vin;
	double Vout;
	double Pr;
	double Amu;
	double Amup;
	double Wfu;
	double WOx;
	double WPr;
	double HFu;
	double ArFu;
	double Rswm;
	double Flowin;
	double Fl;

/*******************************************************************************
* functon name:	fnGrid
* description:	in this file you should do:
*			--initialize variable you defined
*			--set the filename of the output file, the step can be ignored 
*			--set the title of the variable that will be calculated
*			--set which variable will be calc and print
*			--set axis mode-nMode: 1-cartesian, 2-axisymmetric, 3-polar
*			--set the number of steps
*			--set the number of grids
*			--devide grids
********************************************************************************/
void fnGrid()
{
	int i;

	Flox=F[3];
	Fuel=F[4];
	Ox=F[5];
	H=F[6];
	T=F[7];
	Rsw=F[8];

	szTitle[0]="U";
	szTitle[1]="V";
	szTitle[2]="Str Fun";
	szTitle[4]="Fuel";
	szTitle[5]="Oxygen";
	szTitle[7]="T";
	szTitle[8]="R*Th";
	szTitle[10]="pressure";
	szTitle[11]="Density";
	Relax[0]=0.5;
	Relax[1]=0.5;
	Relax[10]=0.8;
	Relax[11]=0.5;
	for( i=0;i<=8;i++)
	{
		bSolve[i]=TRUE;
		bPrint[i]=TRUE;
	}
	bSolve[4]=FALSE;
	bSolve[5]=FALSE;
	bSolve[7]=FALSE;
	bPrint[3]=FALSE;
	bPrint[6]=FALSE;
	bPrint[10]=TRUE;
	bPrint[11]=TRUE;
	iLast=30;
	nMode=2;
	XL=1;
	YL=0.2;
	L1=6;//
	M1=9;//
	R[0]=0;

	fnUGrid();
}

/*******************************************************************************
* functon name:	fnStart
* description:	in this file you should do:
*			--initialize the 
********************************************************************************/
void fnStart()
{
	int i,j;
	double RATIO;
	Ai=4;
	Cp=1000;
    Tin=300;
    RhoCon=10;
    Wfu=16;
    WOx=32;
    WPr=28;
    HFu=3e7;
    ArFu=20;
	for( i=0;i<=L1;i++)
	{
		for( j=0;j<=M1;j++)
		{
			U[i][j]=20;
			U[1][j]=0;
			U[i][M1]=0;
			V[i][j]=0;
			Rsw[i][j]=0;
			P[i][j]=0;
			Flox[i][j]=-0.232/Ai;
			H[i][j]=Cp*Tin;
		}
	}
    Amu=0.1;
    Pr=0.7;
    Amup=Amu/Pr;
    U[1][M2]=20;
    U[1][2]=50;
    RATIO=1/(Wfu*(0.232/WOx+0.768/WPr));
    U[1][1]=(U[1][2]*YCVR[2]+U[1][M2]*YCVR[M2])*RATIO/(YCVR[1]*ArFu);
    H[0][1]=Cp*Tin+HFu;
    Flox[0][1]=1;
    Rsw[0][2]=R[2]*U[1][2]*0.2;
}

/*******************************************************************************
* functon name:	fnDense
* description:	in this file you should do:
*			--calc the density if needed
********************************************************************************/
void fnDense()
{
	int i,j;
	double PROD,WAV,DEN;
	for( i=0;i<=L1;i++)
	{
		for( j=0;j<=M1;j++)
		{
			Fuel[i][j]=__max(0,Flox[i][j]);
			Ox[i][j]=__max(0,-Ai*Flox[i][j]);
			T[i][j]=(H[i][j]-Fuel[i][j]*HFu)/Cp;
			PROD=1-Fuel[i][j]-Ox[i][j];
			WAV=1/(Fuel[i][j]/Wfu+Ox[i][j]/WOx+PROD/WPr);
			DEN=RhoCon*WAV/T[i][j];
			if(iIter==0)
				Rho[i][j]=DEN;
			Rho[i][j]=DEN*Relax[nRho]+Rho[i][j]*(1-Relax[nRho]);
		}
	}

}

/*******************************************************************************
* functon name:	fnbound
* description:	in this file you should do:
*			--calc the boundary
********************************************************************************/
void fnBound()
{
	int j;
	double  FACTOR;
	if(iIter==0)
	{
		Flowin=0;
		for( j=1;j<=M2;j++)
			Flowin+=Rho[0][j]*U[1][j]*YCVR[j];
	}
	Fl=0;
    for( j=1;j<=M2;j++)
	{
		Fl+=Rho[L1][j]*U[L2][j]*YCVR[j];
		H[L1][j]=H[L2][j];
		Flox[L1][j]=Flox[L2][j];
	}
     FACTOR=Flowin/Fl;
    for(j=1;j<=M2;j++)
		U[L1][j]=U[L2][j]*FACTOR;
}

/*******************************************************************************
* functon name:	fnOutput
* description:	in this file you should do:
*			--print calc info to terminal or write to file
********************************************************************************/
void fnOutPut()
{
	if(iIter==0)
		fprintf(pFileOut,"\n    Iter    SMax    SSum      U[4][2]      T[4][2]   Ox[4][6]");
	fprintf(pFileOut,"\n   %4d  %10.3e  %10.3e  %10.3e  %10.3e   %10.3e",iIter,SMax,SSum,U[3][1],T[3][1],Ox[3][5]);

	if(iIter==iLast)
		fnPrint();
}

/*******************************************************************************
* functon name:	fnGamsor
* description:	in this file you should do:
*			--set the gam 
*			--set the sc,sp
********************************************************************************/
void fnGamSor()
{
	int i,j;
	double RhoM,AR;
	for( i=0;i<=L1;i++)
	{
		for( j=0;j<=M1;j++)
		{
			Gam[i][j]=Amu;
			if(nF==3||nF==6)
			{
				Gam[i][j]=Amup;
				Gam[i][M1]=0;
				Gam[0][j]=0;
			}
			Gam[L1][j]=0;
		}
	}
    if(nF==1)
	{
		for( i=1;i<=L2;i++)
		{
			for( j=2;j<=M2;j++)
			{
				Rswm=FY[j]*Rsw[i][j]+FYM[j]*Rsw[i][j-1];
				RhoM=FY[j]*Rho[i][j]+FYM[j]*Rho[i][j-1];
				Con[i][j]=RhoM*Rswm*Rswm/(RMN[j]*RMN[j]*RMN[j]);
				AP[i][j]=-Amu/(RMN[j]*RMN[j]);
			}
		}
	}
	if(nF==8)
	{
		for( i=1;i<=L2;i++)
		{
			for( j=1;j<=M2;j++)
			{
				AR=2*Amu/YCVR[j];
				Con[i][j]=AR*Rsw[i][j-1];
				AP[i][j]=-AR;
			}
		}
	}
}