setup1.c

数值传热学NHT C语言源代码及解释 陶文铨院士 的经典例题中3个题目的解答

#include "global_var.h" 

// ugrid grid 以及setup1的铅版部分设定了网格
//这一部分可以参考图3.


void generate_grid()
{   int i,j;
    user_generate_velocity_grid();
	l2=l1-1;
	l3=l2-1;
	m2=m1-1;
	m3=m2-1;
    x[1]=xu[1];            //设定主控制体节点的x坐标    	
    for (i=2;i<=l2;i++)    
	{
		x[i]=0.5*(xu[i+1]+xu[i]);  
    }
	x[l1]=xu[l1];

	y[1]=yv[2];           ////设定主控制体节点的y坐标    	
    for (j=2;j<=m2;j++)    
	{
		y[j]=0.5*(yv[j+1]+yv[j]);  
    }
    y[m1]=yv[m1];

//xdif xcvs xcvi xcvip设定。参考图4	
	for (i=2;i<=l1;i++)
	{
		xdif[i]=x[i]-x[i-1];        
	}

	for (i=2;i<=l2;i++)
	{
		xcv[i]=xu[i+1]-xu[i];      
	}

	xcvs[3]=xdif[2]+xdif[3];  
    for (i=4;i<=(l2-1);i++)
    {
	   xcvs[i]=xdif[i];
	}
    xcvs[l2]=xdif[l2]+xdif[l1];

	for (i=3;i<=l3;i++)
   {
	   xcvi[i]=0.5*xcv[i];
	   xcvip[i]=xcvi[i];
   }
   xcvip[2]=xcv[2];
   xcvi[l2]=xcv[l2];

//ycv ydif ycvs ycvr ycvrs arx arxj arxjp设定。参考图5与图6
//原来FORTRAN的代码对于拄坐标与极坐标“似乎”是在直角坐标系加以校正
//校正因子包括sx sxmn, sx,sxmn又与r rmn相关
//不过在直角坐标系中，这四个参数都是1，为了代码的简洁。将它们都剔除了。


   for (j=2;j<=m1;j++)
   {
	   ydif[j]=y[j]-y[j-1];
   }

   for (j=2;j<=m2;j++)
   {
	   ycv[j]=yv[j+1]-yv[j];
   }
   
   ycvs[3]=ydif[2]+ydif[3];
   for (j=4;j<m2;j++)
   {
	   ycvs[j]=ydif[j];
   }
   ycvs[m2]=ydif[m2]+ydif[m1];

 
for (j=2;j<=m2;j++)
	{
	ycvr[j]=ycv[j];
	arx[j]=ycvr[j];
	}

   ycvrs[3]=ycvs[3];
	for (j=3;j<=m3;j++)
	{
		ycvrs[j]=ydif[j];
	}
   ycvrs[m2]=ycvs[m2];

   arxjp[2]=arx[2];
	for (j=3;j<m2;j++)
	{
		arxj[j]=0.5*arx[j];
		arxjp[j]=arxj[j];
	}
	arxj[m2]=arx[m2];

	if (mode!=1) 
	{
		printf("mode <>1. currently only x-y axis system supported\n");
		exit(0);
	}
//**********y方向上的变量设定结束***************

//插职的比例系数fx fxm fy fym设定
	for (j=3;j<=m3;j++)
	{
		fv[j]=arxjp[j]/arx[j];
		fvp[j]=1.0-fv[j];
	}

	for (i=3;i<=l2;i++)
	{
		fx[i]=0.5*xcv[i-1]/xdif[i];
		fxm[i]=1.0-fx[i];
	}
	
	fx[2]=0.0;
	fxm[2]=1;
	fx[l1]=1.0;
	fxm[l1]=0.0;


   for (j=3;j<=m2;j++)
   {
	   fy[j]=0.5*ycv[j-1]/ydif[j];
	   fym[j]=1.0-fy[j];
   }

   fy[2]=0.0;
   fym[2]=1.0;
   fy[m1]=1.0;
   fym[m1]=0.0;

//相同缺省将校正压力、速度、压力这些量设置为零

  for (j=1;j<=m1;j++)
  	  for (i=1;i<=l1;i++)
	  {
		  pc(i,j)=0;   //#define pc(i,j) (f[i][j][3])
          u(i,j)=0;    //#define u(i,j) (f[i][j][1])
 		  v(i,j)=0;    //#define v(i,j) (f[i][j][2])
		  p(i,j)=0;    //#define p(i,j)  (f[i][j][11])
		  con[i][j]=0;
		  ap[i][j]=0;
		  rho(i,j)=rhocon;
        }
  
}