exam4.c

CSIMPLE2d CFD C源程序

/*
* exam4.c : example 4 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"

#define PI 3.14159
/*global variable for user */
	double (*T)[NSIZE];
	double (*W)[NSIZE];
	double Amu;
	double Cond;
	double Cp;
	double RhoCp;
	double dPdZ;
	double dTdZ;
	double Tw;
	double Tb;
	double WBAR;
	double Re;
	double FRe;
	double ANu;

/*******************************************************************************
* 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()
{
	W=F[3];
	T=F[4];
	Amu=1;
	Cond=1;
	RhoCp=1;
	dPdZ=-2000;
	dTdZ=100;
	Tw=1;

	szTitle[3]="W/WBAR";
	szTitle[4]="Theta";
	bSolve[3]=TRUE;
	bPrint[3]=TRUE;
	bPrint[4]=TRUE;
	iLast=6;
	nMode=3;
	nTimes[3]=4;
	nTimes[4]=4;
	XL=15*PI/180;
	YL=1;
	R[0]=1;
	L1=6;//
	M1=6;//

	fnUGrid();
}

/*******************************************************************************
* functon name:	fnStart
* description:	in this file you should do:
*			--initialize the 
********************************************************************************/
void fnStart()
{
	int i,j;
	for(i=0;i<=L1;i++)
	{
		for(j=0;j<=M1;j++)
		{
			W[i][j]=0;
			T[i][j]=Tw;
		}
	}
}

/*******************************************************************************
* functon name:	fnDense
* description:	in this file you should do:
*			--calc the density if needed
********************************************************************************/
void fnDense()
{

}

/*******************************************************************************
* functon name:	fnbound
* description:	in this file you should do:
*			--calc the boundary
********************************************************************************/
void fnBound()
{
	int i,j;
	double Asum,Wsum,Tsum,AR,WR,DH,Wp,Htp,Htc;

    Asum=0;
    Wsum=0;
    Tsum=0;
    for( i=1;i<=L2;i++)
	{
		for( j=1;j<=M2;j++)
		{
			if(i!=1||j>3)
			{
				AR=XCV[i]*YCVR[j];
				WR=W[i][j]*AR;
				Wsum+=WR;
				Asum+=AR;
				Tsum+=WR*T[i][j];
			}
		}
	}
    WBAR=Wsum/Asum;
    Wp=(R[0]+R[M1])*XL+(1+XCV[1])*(RMN[4]-R[0]);
    DH=4*Asum/Wp;
    Re=RhoCon*WBAR*DH/Amu;
    FRe=-2*dPdZ*DH/(RhoCon*WBAR*WBAR+1e-30)*Re;
    Tb=Tsum/(Wsum+1e-30);
	Htp=Wp-R[M1]*XL;
	Htc=RhoCp*Wsum*dTdZ/((Tw-Tb+1e-30)*Htp);
    ANu=Htc*DH/Cond;
    if(iIter==3)
	{
      bSolve[4]=TRUE;
      bSolve[3]=FALSE;
	}
}

/*******************************************************************************
* functon name:	fnOutput
* description:	in this file you should do:
*			--print calc info to terminal or write to file
********************************************************************************/
void fnOutPut()
{
	int i,j;
	if(iIter==0)
	{
		fprintf(pFileOut,"\n   Iter     F*Re    Nu");
	}
	fprintf(pFileOut,"\n   %4d   %11.3e   %11.3e ",iIter,FRe,ANu);

	if(iIter==iLast)
	{
	    for( i=1;i<=L1;i++)
		{
			for( j=1;j<=M1;j++)
			{
				W[i][j]/=(WBAR+1e-30);
				T[i][j]=(T[i][j]-Tw)/(Tb-Tw+1e-30);
			}
		}
		fnPrint();
	}
}

/*******************************************************************************
* functon name:	fnGamsor
* description:	in this file you should do:
*			--set the gam 
*			--set the sc,sp
********************************************************************************/
void fnGamSor()
{
	int i,j;

	for( i=0;i<=L1;i++)
	{
		for( j=0;j<=M1;j++)
		{
			Gam[i][j]=Amu;
			if(nF==4)//4
				Gam[i][j]=Cond;
			Gam[L1][j]=0;
			Gam[0][j]=0;
			if(nF==4)
				Gam[i][M1]=0;
			if(j<4)
				Gam[1][j]=1e10;
			if(nF==4&&j<=3)
				Gam[1][j]=10*Cond;
		}
	}
	for(i=1;i<=L2;i++)
	{
		for( j=1;j<=M2;j++)
		{
			Con[i][j]=-dPdZ;
			if(nF==4)
				Con[i][j]=-dTdZ*W[i][j]*RhoCp;
		}
	}

}