rgwshjhs.cpp

对日光温室内的热湿小气候环境进行模拟预测。在确定的温室材料、构造、室外气候条件（气温、湿度、太阳辐射、风速等）以及温室的管理方式条件下

//日光温室热环境模拟程序的函数
//本文件为必用函数
//2007年9月


//全局变量定义-------------------------------------------------------------------------------------------

int AAI[30];
double AAD[80];

int n,nn,mm,nw;
double kd,dtao,tf0,tfn,alf0,alfn;
double dw0[10],lmd0[10],ro0[10],c0[10];
double dx[102],Dx[102],A[102],K[102],P[102],S[102],T[102];
double Dtx[33],Dty[14],dtx[33],dty[14];
double lmdw[14],lmdx[14],lmdn[14],row[14],rox[14],ron[14],cw[14],cx[14],cn[14];
double AA[33][14],BB[33][14],KK[33][14],PP[33][14],SS[33][14],TS[33][14];
double tfl[14],tfr[14],tfu[33],tfd[33],alfl[14],alfr[14],alfu[33],alfd[33];
double Pbi0[33],Pbim[33],Pbj0[14],Pbjn[14],Pc00,Pcn0,Pcm0,Pcnm;
double TT[462],DD[462],RR[462][462];

double Bgh[18],Hgh[18];



//以下为墙体传热相关函数
//-------------------------------------------------------------------------------------------------------


//-------------------------------------------------------------------------------------------------------
//墙体节点划分函数jdhf()---------------------------------------------------------2007.9.2调试
//执行本函数的结果返回墙体计算域的内节点数量n，各节点间距放入dx[i]，编号顺序从室外表面到室内表面；
//在计算域的边界处及各层材料的交界面处设置节点；
//在内外表面处的节间距最小(为0.001m)，离表面越远节间距越大(按0.001m递增)；
//但在各层材料的交界面两侧的节间距有可能小于相邻节间距；

//调用本函数的实参依次为
//nw -- 墙体不同材料分层数量(nw=1～10)，本程序按最大10层编制；
//dw0[iw] -- 各分层墙体厚度(iw=0～nw-1,从外至内)，总厚度一般限制在2.0m以下，要求为0.001m的整倍数，m；
//															位于内外表面的材料最小厚度限制为0.002m；											

//其他外部变量说明
//n -- 计算域内节点数量，本程序按最大100个内节点(加外节点共102个节点)编制，如需更多节点需修改数组定义；
//                                                          n符号并未在本函数中出现，是本函数返回的值
//dx[i] -- 节点i与节点i+1之间的节间距离(i=0～n)，m；

//局部变量说明
//xd[i] -- 从墙内表面开始划分节点时的节点i与节点i+1之间的节间距离，m；
//imo -- 墙体外半部分节间数；
//imi -- 墙体内半部分节间数；
//di -- 已划分节点的墙体厚度，m；
//dd -- 按正常节间距离，预计下一节点划分后达到的墙体厚度，m；
//dm -- 已划分节点和正在划分节点的材料层的累计厚度，m；
//dw -- 墙体总厚度，m；
//dxmax -- 从墙表面往墙中依次划分节点中出现过的最大节间距，m；
//e -- 一个很小的数(=1e-10)，加入条件判断式中，避免因计算机计算位数舍入误差造成数值比较条件判断的错误；
//---------------------------------------------------------------------------------------------
int jdhf(int nw, double dw0[])
{
	int i,iw,imo,imi;		double di,dm,dw,dd,dxmax,e;
	double xd[51];

	e=1e-10;

	dw=0.0;
	for(iw=0;iw<nw;iw++)dw=dw+dw0[iw];						//根据各层材料厚度计算墙体总厚度

//外半部墙体节点划分--------------------
	dx[0]=0.001;	dxmax=dx[0];	di=0.001;	iw=0;	dm=dw0[iw];			//设置起点的数值
	for(i=1;  ;i++)
	{
		dd=di+dxmax+0.001;		//按正常节间距离(dxmax+0.001)，预计下一节点划分后达到的墙体厚度
		if(dd+e<dw/2.0)
		{
			if(dd+e<dm)						//正常划分节点的情况(节点仍在本材料层)
			{
				dx[i]=dxmax+0.001;
				dxmax=dx[i];
				di=dd;
				continue;
			}
			if(dd+e>=dm)					//如按正常节间距离划分，节点将到达或进入下一材料层的情况
			{
				dx[i]=dm-di;				//该节点划到材料边界处
				di=dm;
				iw=iw+1;
				dm=dm+dw0[iw];
				continue;
			}
		}
		if(dd+e>=dw/2.0)					//如按正常节间距离划分，节点将到达或进入内半墙的情况
		{
			if(dd-e<=dm)					//节点不会进入下一材料层
			{
				dx[i]=dw/2.0-di;		//该节点划到内外半墙的交界处(墙体总厚度中点)
				imo=i+1;
				break;
			}
			if(dd-e>dm)					//如按正常节间距离划分，节点将进入下一材料层中的情况
			{
				if(dm+e<dw/2.0)			//墙体总厚度中点进入下一材料层中的情况
				{
					dx[i]=dm-di;
					di=dm;
					iw=iw+1;
					dm=dm+dw0[iw];
					continue;
				}
				if(dm+e>=dw/2.0)			//墙体总厚度中点仍在本材料层中的情况
				{
					dx[i]=dw/2.0-di;	//该节点划到内外半墙的交界处(墙体总厚度中点)
					imo=i+1;
					break;
				}
			}
		}
	}

//内半部墙体节点划分--------------------
	xd[0]=0.001;	dxmax=xd[0];	di=0.001;		iw=nw-1;	dm=dw0[iw];
	for(i=1;  ;i++)
	{
		dd=di+dxmax+0.001;
		if(dd+e<dw/2.0)
		{
			if(dd+e<dm)
			{
				xd[i]=dxmax+0.001;
				dxmax=xd[i];
				di=dd;
				continue;
			}
			if(dd+e>=dm)
			{
				xd[i]=dm-di;
				di=dm;
				iw=iw-1;
				dm=dm+dw0[iw];
				continue;
			}
		}
		if(dd+e>=dw/2.0)
		{
			if(dd-e<=dm)
			{
				xd[i]=dw/2.0-di;
				imi=i+1;
				break;
			}
			if(dd-e>dm)
			{
				if(dm+e<dw/2.0)
				{
					xd[i]=dm-di;
					di=dm;
					iw=iw-1;
					dm=dm+dw0[iw];
					continue;
				}
				if(dm+e>=dw/2.0)
				{
					xd[i]=dw/2.0-di;
					imi=i+1;
					break;
				}
			}
		}
	}

	for(i=imo;i<=imo+imi;i++) dx[i]=xd[imo+imi-1-i];	//将xd[]颠倒顺序放入dx[]

	return imo+imi-1;									//返回内节点的数量n
}
//-------------------------------------------------------------------------------------------------------


//-------------------------------------------------------------------------------------------------------
//一维导热差分方程组静态参数计算函数ywdrcfjc()------------------------------------2007.9.2调试
//本函数确定各节点控制区长度Dx[]、材料参数lmd[]，ro[]，c[]，
//				并计算一维导热差分方程组中各静态参数A[]，K[]，P[]，设置初始温度，控制区域内热源全预置0

//调用本函数的实参依次为
//n -- 计算域内节点数量，本程序按最大100个内节点(加外节点共102个节点)编制，如需更多节点需修改数组定义；
//dx[i] -- 节点i与节点i+1之间的距离(i=0～n)，m；
//dw0[j] -- 各分层墙体厚度(j=0～nw-1,从外至内)，总厚度一般限制在2.0m以下，要求为0.001m的整倍数，m；
//                                                             位于内外表面的材料最小厚度限制为0.002m；											
//lmd0[j] -- j层材料的导热系数(j=0～nw-1)，W/m℃；
//ro0[j] -- j层材料的密度(j=0～nw-1)，kg/m3；
//c0[j] -- j层材料的比热容(j=0～nw-1)，J/kg℃；
//dtao -- 时间步长，s；

//其他外部变量说明
//Dx[i] -- i节点代表的控制区域长度(i=0～n+1)，
//                           内节点Dx[i]=(dx[i-1]+dx[i])/2，边节点Dx[0]=dx[0]/2，Dx[n+1]=dx[n]/2，m；
//A[i] -- 离散方程中的系数(i=0～n)；
//K[i] -- 离散方程中的系数(i=0～n+1)；
//P[i] -- 离散方程中的系数(i=1～n)；
//T[i] -- i节点温度(i=0～n+1)，℃；
//S[i] -- i节点代表的控制区域内热源(i=0～n+1)，W/m3；

//局部变量说明(以下只是对变量、下标等符号意义的说明，具体处理已在函数中作出，不需在程序使用中修改语句)
//lmd[i] -- i节点处导热系数(i=0～n)，(n+1节点的导热系数lmd[n+1]程序中不会用到)
//          控制区域内有二种材料时，节点取在二种材料交界点上，lmd[i]为i与i+1间的材料导热系数，W/m℃；
//ro[i] -- i节点代表的控制区域材料密度(i=0～n)，
//             n+1节点的密度ro[n+1]取作=ro[n]，所以程序中不出现ro[n+1]，直接使用ro[n]；
//                控制区域内有二种材料时，节点取在二种材料交界点上, ro[i]为i与i+1间的材料密度，kg/m3；
//c[i] -- i节点代表的控制区域材料比热容(i=0～n+1)，
//             n+1节点的比热容c[n+1]取作=c[n]，所以程序中不出现c[n+1]，直接使用c[n]；
//                控制区域内有二种材料时，节点取在二种材料交界点上, c[i]为i与i+1间的材料比热容，J/kg℃；
//        i节点为二种材料交界点时，c与ro的乘积取为二侧节距间材料c与乘积的加权平均值，
//dd -- i节点(包括i节点)之前节间距之和，m；
//dm -- j层材料之前(包括j层)之前材料厚度之和，m；
//e -- 一个很小的数(=1e-10)，加入条件判断式中，避免因计算机计算位数舍入误差造成数值比较条件判断的错误；
//---------------------------------------------------------------------------------------------
void ywdrcfjc(int n,double dx[],double dw0[],double lmd0[],double ro0[],double c0[],double dtao)
{
	double dd,dm,e;		int i,j;
	double lmd[102],ro[102],c[102];

	e=1e-10;
//存入各节点材料导热系数、密度、比热等参数
	dd=0.0;		j=0;	dm=dw0[j];
	for(i=0;i<=n;i++)
	{
		dd=dd+dx[i];
		if(dd-e>dm){j=j+1;	dm=dm+dw0[j];}
		lmd[i]=lmd0[j];	ro[i]=ro0[j];	c[i]=c0[j];
	}

//计算各节点控制区域长度
	Dx[0]=dx[0]/2.0;			Dx[n+1]=dx[n]/2.0;
	for(i=1;i<=n;i++) Dx[i]=(dx[i-1]+dx[i])/2.0;

//计算一维导热差分方程组中各静态参数
	for(i=0;i<=n;i++) A[i]=lmd[i]/dx[i];
	K[0]=ro[0]*c[0]*Dx[0]/dtao;
	K[n+1]=ro[n]*c[n]*Dx[n+1]/dtao;
	for(i=1;i<=n;i++) K[i]=(ro[i-1]*c[i-1]*dx[i-1]+ro[i]*c[i]*dx[i])/2.0/dtao;
	for(i=1;i<=n;i++) P[i]=K[i]+A[i]+A[i-1];

//设置初始温度，控制区域内热源全预置0
	for(i=0;i<=n+1;i++)T[i]=15.0;
	for(i=0;i<=n+1;i++)S[i]=0.0;
}
//-------------------------------------------------------------------------------------------------------


//-------------------------------------------------------------------------------------------------------
//求解三对角线方程组的函数atrde()(《C常用算法程序集》)
//本函数求得方程组的解，放在数组d[]中；

//调用本函数的实参依次为
//b[] -- 顺序存放三对角矩阵中三条对角线上的元素，逐行顺序存放，
//                        长度为：3*方程阶数-2，即：3*节点总数-2，或：3*(内节点数+2)-2，
//                                                             本程序内节点数最大100，所以最大长度304；
//n -- 方程组的阶数，即节点总数(注意本程序其他地方n是表示内节点数，加上2才是本处的节点总数n)；
//m -- 三对角矩阵中三条对角线上元素的个数，即b[]的长度；
//d[] -- 常数向量中的各元素，长度等于方程组的阶数，即总节点数；函数返回时存放方程组的解。
//---------------------------------------------------------------------------------------------
int atrde(double b[], int n, int m, double d[])
{ int k,j;
  double s;
  if (m!=(3*n-2))
    { printf("err\n"); return(-2);}
  for (k=0;k<=n-2;k++)
    { j=3*k; s=b[j];
      if (fabs(s)+1.0==1.0)
        { printf("fail\n"); return(0);}
      b[j+1]=b[j+1]/s;
      d[k]=d[k]/s;
      b[j+3]=b[j+3]-b[j+2]*b[j+1];
      d[k+1]=d[k+1]-b[j+2]*d[k];
    }
  s=b[3*n-3];
  if (fabs(s)+1.0==1.0)
    { printf("fail\n"); return(0);}
  d[n-1]=d[n-1]/s;
  for (k=n-2;k>=0;k--)
    d[k]=d[k]-b[3*k+1]*d[k+1];
  return(2);
}  
//-------------------------------------------------------------------------------------------------------


//-------------------------------------------------------------------------------------------------------
//一维温度传递函数T0_T()
//由上一时刻温度计算下一时刻温度，开始时温度存放于T[]中，函数执行结束时返回后一时刻温度于T[]中；

//本函数执行将调用求解三对角线方程组的函数atrde()

//调用本函数的实参依次为