exam11.c

CSIMPLE2d CFD C源程序

/*
* exam11.c : example 11 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 (*T)[NSIZE];
	double Tin;
	double Tair;
	double Uin;
	double Vout;
	double Pr;
	double Amu;
	double Amup;
	double Depth;
	double Flow;
	double Htc;
	double Cp;
	double Fdz;

/*******************************************************************************
* 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;

	T=F[3];

	szTitle[0]="U";
	szTitle[1]="V";
	szTitle[2]="Str Fun";
	szTitle[3]="temprature";
	szTitle[10]="pressure";
	Relax[0]=0.5;
	Relax[1]=0.5;
	for( i=0;i<=2;i++)
	{
		bSolve[i]=TRUE;
		bPrint[i]=TRUE;
	}
	bPrint[2]=FALSE;
	bPrint[10]=TRUE;
	bPrint[3]=TRUE;
	iLast=40;
	XL=1;
	YL=2.0;
	L1=6;//
	M1=11;//

	fnUGrid();
}

/*******************************************************************************
* functon name:	fnStart
* description:	in this file you should do:
*			--initialize the 
********************************************************************************/
void fnStart()
{
	int i,j;
	Depth=0.2;
	Flow=4;
	Uin=100;
	Tin=500;
	Tair=100;
	Htc=2;
	for( i=0;i<=L1;i++)
	{
		for( j=0;j<=M1;j++)
		{
			U[i][j]=0;
			V[i][j]=Vout;
			T[i][j]=Tair;
		}
	}
	Pr=5;
	Cp=1.2;
	Amu=1;
	Amup=Amu/Pr;
	Fdz=Flow/Depth;
}

/*******************************************************************************
* 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 n;
	if(iIter==35)
	{
		for( n=0;n<=2;n++)
			bSolve[n]=FALSE;
		bSolve[3]=TRUE;
	}
}

/*******************************************************************************
* 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      V[4][7]      U[6][4]    T[5][9]");
	fprintf(pFileOut,"\n   %4d  %10.3e  %10.3e  %10.3e  %10.3e   %10.3e ",iIter,SMax,SSum,V[3][6],U[5][3],T[4][8]);

	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 SOR;
 	for( i=0;i<=L1;i++)
	{
		for( j=0;j<=M1;j++)
		{
			Gam[i][j]=Amu;
			if(nF==3)
			{
				Gam[i][j]=Amup;
				Gam[i][0]=0;
				Gam[i][M1]=0;
				Gam[0][j]=0;
				Gam[L1][j]=0;
			}
		}
	}
	if(nF==0)    
	{
		SOR=0.5*Fdz/(XCVS[3]*YCV[3]);
		Con[3][3]=SOR*Uin;
		AP[3][3]=-SOR;
		SOR=0.5*Fdz/(XCVS[4]*YCV[3]);
		Con[4][3]=SOR*Uin;
		AP[4][3]=-SOR;
	}
	else if(nF==1) 
	{
		SOR=0.5*Fdz/(XCV[3]*YCVS[3]);
		AP[3][3]=-SOR;
		SOR=0.5*Fdz/(XCV[3]*YCVS[4]);
		AP[3][4]=-SOR;
	}
	else if(nF==2)
	{
		Con[3][3]=Fdz/(XCV[3]*YCV[3]);
		Con[4][8]=-Fdz/(XCV[4]*YCV[8]);
	}
	else if(nF==3)
	{
		SOR=Htc/(Cp*Depth);
  		for( i=1;i<=L2;i++)
		{
			for( j=1;j<=M2;j++)
			{
				Con[i][j]=SOR*Tair;
				AP[i][j]=-SOR;
			}
		}
		SOR=Fdz/(XCV[3]*YCV[3]);
		Con[3][3]+=SOR*Tin;
		AP[3][3]+=-SOR;
	}
}