rgwshjhs.cpp

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

//n -- 计算域内节点数量，本程序按最大100个内节点(加外节点共102个节点)编制，如需更多节点需修改数组定义；
//alf0,alfn -- 为二边界处的对流换热系数，对于第一类边界条件，相应对流换热系数取一个很大值即可，W/m2℃；
//tf0,tfn -- 依次为二边界处的流体温度，对于第一类边界条件，相应流体温度取为等于边界温度，℃；
//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；
//T[i] -- i节点温度(i=0～n+1)，本函数调用前为前一时刻的温度，调用后为计算时刻的温度，℃；
//S[i] -- i节点代表的控制区域内热源(i=0～n+1)，W/m3；
//A[i] -- 离散方程中的系数(i=0～n)；
//K[i] -- 离散方程中的系数(i=0～n+1)；
//P[i] -- 离散方程中的系数(i=1～n)；

//局部变量说明

//b[] -- 顺序存放三对角矩阵中三条对角线上的元素，逐行顺序存放，
//                        长度为：3*方程阶数-2，即：3*节点总数-2，或：3*(内节点数+2)-2，
//                                                              本程序内节点数最大100，所以最大长度304；
//d[] -- 常数向量中的各元素，长度等于方程组的阶数，即总节点数，三对角方程计算结束时存放方程组的解。
//本函数还可改进，将b[]中不变的部分放在静态参数计算中去，还可减少重复计算，提高运行速度。
//---------------------------------------------------------------------------------------------
void T0_T(int n,double alf0,double alfn,double tf0,double tfn,double Dx[],double T[],double S[],double A[],double K[],double P[])
{	int i;
	double b[306],d[102];

	b[0]=A[0]+K[0]+alf0;				b[1]=-A[0];
	for (i=1;i<=n;i++)
		{
			b[3*i-1]=-A[i-1];
			b[3*i]=P[i];
			b[3*i+1]=-A[i];
		}
	b[3*n+2]=-A[n];		b[3*n+3]=A[n]+K[n+1]+alfn;
	
	d[0]=K[0]*T[0]+alf0*tf0+S[0]*Dx[0];
	d[n+1]=K[n+1]*T[n+1]+alfn*tfn+S[n+1]*Dx[0];
	for(i=1;i<=n;i++) d[i]=K[i]*T[i]+S[i]*Dx[i];
	atrde(b,n+2,3*n+4,d);
	for(i=0;i<=n+1;i++) T[i]=d[i];
}
//-------------------------------------------------------------------------------------------------------



//以下为地面传热相关函数
//-------------------------------------------------------------------------------------------------------


//-------------------------------------------------------------------------------------------------------
//地下网格划分函数wghf()---------------------------------------------------------2007.9.7调试
//执行本函数的结果，将水平方向各节点间距放入dtx[i](i=0～31)，编号顺序从北到南；
//                  将垂直方向各节点间距放入dty[j](j=0～12)，编号顺序从地面往下；
//在计算域的边界处及各层材料的交界面处设置节点；
//水平方向从北到南分为北室外、墙下、室内、南室外共4段，每段相交处节间距较小，距离相交处越远节间距越大；
//垂直方向接近地面处节间距较小，距离地面越远节间距越大；

//调用本函数的实参依次为
//kd -- 日光温室跨度(净跨度，即从墙内表面至屋面南部根基处)，m；
//nw -- 墙体不同材料分层数量(nw=1～10)，本程序按最大10层编制；
//dw0[iw] -- 各分层墙体厚度(iw=0～nw-1,从外至内)，m；

//其他外部变量说明
//nn -- 计算域水平方向内节点数量，本程序按31个内节点(加上外节点共33个节点)编制，nn符号并未在本函数中出现；
//mm -- 计算域垂直方向内节点数量，本程序按12个内节点(加上外节点共14个节点)编制，mm符号并未在本函数中出现；
//dtx[i] -- 水平方向节点i与节点i+1之间的节间距离(i=0～n)，m；
//dty[j] -- 垂直方向节点j与节点j+1之间的节间距离(j=0～m)，m；

//局部变量说明
//dw -- 墙体总厚度，m；

//-------------------------------------------------------------------------------------------------------
void wghf(double kd, int nw, double dw0[])
{
	double dw;		int iw;

	dw=0.0;
	for(iw=0;iw<nw;iw++)dw=dw+dw0[iw];						//计算墙体总厚度

	dtx[0]=1.800;		dtx[1]=1.000;		dtx[2]=0.600;		dtx[3]=0.320;		//北面室外
	dtx[4]=0.180;		dtx[5]=0.100;

	dtx[6]=0.18*dw;	dtx[7]=0.32*dw;	dtx[8]=0.32*dw;	dtx[9]=0.18*dw;					//墙下

	dtx[10]=0.010*kd;	dtx[11]=0.015*kd;	dtx[12]=0.025*kd;	dtx[13]=0.040*kd;	//室内北半部
	dtx[14]=0.060*kd;	dtx[15]=0.090*kd;	dtx[16]=0.120*kd;	dtx[17]=0.140*kd;

	dtx[18]=0.140*kd;	dtx[19]=0.120*kd;	dtx[20]=0.090*kd;	dtx[21]=0.060*kd;	//室内南半部
	dtx[22]=0.040*kd;	dtx[23]=0.025*kd;	dtx[24]=0.015*kd;	dtx[25]=0.010*kd;

	dtx[26]=0.100;		dtx[27]=0.180;												//南面室外
	dtx[28]=0.320;		dtx[29]=0.600;		dtx[30]=1.000;		dtx[31]=1.800;

	dty[0]=0.002;	dty[1]=0.004;	dty[2]=0.008;	dty[3]=0.016;	dty[4]=0.025;	//深度方向
	dty[5]=0.045;	dty[6]=0.080;	dty[7]=0.120;	dty[8]=0.200;	dty[9]=0.300;
	dty[10]=0.500;	dty[11]=0.800;	dty[12]=1.200;
}
//-------------------------------------------------------------------------------------------------------


//-------------------------------------------------------------------------------------------------------
//二维导热差分方程组静态参数计算函数ewdrcfjc()------------------------------------2007.9.28调试
//本函数确定各节点控制区长度Dtx[]、Dty[]，材料参数lmds[][]，ros[][]，cs[][]，
//将存放在节点上的材料热导率值折算为控制区边界处的值lmdip[][]，lmdjp[][]，
//计算各节点材料密度值与比热容之乘积roc[][]，
//计算二维导热差分方程组中各静态参数AA[][]，BB[][]，KK[][]，PP[][]，
//节点温度TS[][]赋给初值(15℃)，控制区域内热源SS[][]按多数取值先全部赋值0.0，
//确定左、右、墙下和地下深处边界条件，tfl[]，tfr[]，tfu[]，tfd[]，alfl[]，alfr[]，alfu[]，alfd[]，

//调用本函数的实参依次为
//nn -- 计算域水平方向内节点数量，本程序按31个内节点(加上外节点共33个节点)编制，nn=31;
//mm -- 计算域垂直方向内节点数量，本程序按12个内节点(加上外节点共14个节点)编制，mm=12;
//dtx[i] -- 水平方向节点i与节点i+1之间的节间距离(i=0～nn)，m；
//dty[j] -- 垂直方向节点j与节点j+1之间的节间距离(j=0～mm)，m；
//lmdw[j] -- 室外第j层(j=0～mm)土壤的导热系数	(室外单元下标i=0～5，26～31)，	W/m℃；
//lmdx[j] -- 墙下第j层(j=0～mm)材料的导热系数	(墙下单元下标i=6～9)，			W/m℃；
//lmdn[j] -- 室内第j层(j=0～mm)土壤的导热系数	(室内单元下标i=10～25)，		W/m℃；
//row[j] --  室外第j层(j=0～mm)土壤的密度		(室外单元下标i=0～5，26～31)，	kg/m3；
//rox[j] --  墙下第j层(j=0～mm)材料的密度		(墙下单元下标i=6～9)，			kg/m3；
//ron[j] --  室内第j层(j=0～mm)土壤的密度		(室内单元下标i=10～25)，		kg/m3；
//cw[j] --   室外第j层(j=0～mm)土壤的比热容		(室外单元下标i=0～5，26～31)，	J/kg℃；
//cx[j] --   墙下第j层(j=0～mm)材料的比热容		(墙下单元下标i=6～9)，			J/kg℃；
//cn[j] --   室内第j层(j=0～mm)土壤的比热容		(室内单元下标i=10～25)，		J/kg℃；
//dtao -- 时间步长，s；

//其他外部变量说明
//Dtx[i] -- i节点代表的控制区域横向长度(i=0～nn+1)，
//                    内节点Dtx[i]=(dtx[i-1]+dtx[i])/2，边节点Dtx[0]=dtx[0]/2，Dtx[nn+1]=dtx[nn]/2，m；
//Dty[j] -- j节点代表的控制区域纵向长度(j=0～mm+1)，
//                    内节点Dty[j]=(dty[j-1]+dty[j])/2，边节点Dty[0]=dty[0]/2，Dty[mm+1]=dty[mm]/2，m；
//AA[i][j] -- 离散方程中的系数(i=0～nn，j=0～mm+1)；
//BB[i][j] -- 离散方程中的系数(i=0～nn+1，j=0～mm)；
//KK[i][j] -- 离散方程中的系数(i=0～nn+1，j=0～mm+1)；
//PP[i][j] -- 离散方程中间节点对应的对角线系数(i=0～nn，j=0～mm)；
//TS[i][j] -- i,j节点温度(i=0～nn+1,j=0～mm+1)，℃；
//SS[i][j] -- i,j节点代表的控制区域内热源(i=0～nn+1,j=0～mm+1)，W/m3；
//tfl[j]，tfr[j] -- 左，右边界处(j=0～mm+1)流体温度，取等于边界处土壤初始温度(以后保持不变)，℃；
//alfl[j]，alfr[j] -- 左，右边界对流换热系数(j=0～mm+1)，取等于0(第二类边界条件(绝热))，W/m2℃；
//tfu[i] -- 上边界处(i=0～nn+1)流体温度，
//                          i=0～6及i=26～32为室外气温，i=7～9为对应节点温度，i=10～25为室内气温，℃；
//tfd[i] -- 下边界处(i=0～nn+1)流体温度，取等于15℃；
//alfu[i]，alfd[i] -- 上，下边界对流换热系数(i=0～nn+1)，W/m2℃；

//Pbi0[i]，Pbim[i]，Pbj0[j]，Pbjn[j] -- 边界条件式中各系数(四边，i=1～nn，j=1～mm)；(可考虑改为局部变量)
//Pc00，Pcn0，Pcm0，Pcnm -- 边界条件式中各系数(四角)；(可考虑改为局部变量)

//TT[j] -- 离散方程矩阵中的温度向量，i,j节点数据逐行依次存放(元素个数=33*14=462)；
//DD[j] -- 离散方程矩阵中的常数向量(元素个数=33*14=462)；
//RR[i][j] -- 离散方程的系数矩阵中各元素(元素个数=462*462=213444)
//                     (i=0～[(nn+2)(mm+2)*(nn+2)(mm+2)-1],j=0～[(nn+2)(mm+2)*(nn+2)(mm+2)-1])。

//局部变量说明(以下只是对变量、下标等符号意义的说明，具体处理已在函数中作出，不需在程序使用中修改语句)
//lmdip[i][j] -- i,j节点(i=0～nn,j=0～mm+1)控制区右边界(i+)的热导率值，W/m℃；
//lmdjp[i][j] -- i,j节点(i=0～nn+1,j=0～mm)控制区下边界(j+)的热导率值，W/m℃；
//lmds[i][j] -- i,j节点代表的控制区域材料热导率(i=0～nn,j=0～mm)，控制区域内有二种材料时，
//                                              i与i+1、j与j+1间材料的热导率存入lmds[i][j]，W/m℃；
//ros[i][j] --  i,j节点代表的控制区域材料密度值(i=0～nn,j=0～mm)，控制区域内有二种材料时，
//                                              i与i+1、j与j+1间材料的密度值存入ros[i][j]，kg/m3；
//cs[i][j] --   i,j节点代表的控制区域材料比热容(i=0～nn,j=0～mm)，控制区内有二种材料时，
//                                              i与i+1、j与j+1间材料的比热容存入cs[i][j]，J/kg℃；
//roc[i][j] -- 为ros[i][j]*cs[i][j]，i,j点位于分界处时，按不同材料的控制区面积取该乘积的加权平均值；
//
//---------------------------------------------------------------------------------------------
void ewdrcfjc(int nn,int mm,double dtx[],double dty[],double lmdw[],double lmdx[],double lmdn[],
			  double row[],double rox[],double ron[],double cw[],double cx[],double cn[],double dtao)
{
	int i,j;
	double lmdip[33][14],lmdjp[33][14],lmds[32][13],ros[32][13],cs[32][13],roc[33][14];	

//存入各节点材料导热系数、密度、比热等参数
	for(i=0;i<6;i++)for(j=0;j<=mm;j++)								//北面室外
		{lmds[i][j]=lmdw[j];	ros[i][j]=row[j];	cs[i][j]=cw[j];}

	for(i=6;i<10;i++)for(j=0;j<=mm;j++)								//墙下
		{lmds[i][j]=lmdx[j];	ros[i][j]=rox[j];	cs[i][j]=cx[j];}

	for(i=10;i<26;i++)for(j=0;j<=mm;j++)							//室内
		{lmds[i][j]=lmdn[j];	ros[i][j]=ron[j];	cs[i][j]=cn[j];}

	for(i=26;i<32;i++)for(j=0;j<=mm;j++)							//南面室外
		{lmds[i][j]=lmdw[j];	ros[i][j]=row[j];	cs[i][j]=cw[j];}

//计算各节点控制区域长度
	Dtx[0]=dtx[0]/2.0;			Dtx[nn+1]=dtx[nn]/2.0;
	Dty[0]=dty[0]/2.0;			Dty[mm+1]=dty[mm]/2.0;

	for(i=1;i<=nn;i++) Dtx[i]=(dtx[i-1]+dtx[i])/2.0;
	for(j=1;j<=mm;j++) Dty[j]=(dty[j-1]+dty[j])/2.0;

//将存放在节点上的材料热导率值折算为控制区边界处的值lmdip,lmdjp
	for(i=0;i<=nn;i++)
	{								//i节点控制区i+边界的热导率值
		lmdip[i][0]=lmds[i][0];		lmdip[i][mm+1]=lmds[i][mm];
		for(j=1;j<=mm;j++)
		lmdip[i][j]=2.0*Dty[j]/(dty[j-1]/lmds[i][j-1]+dty[j]/lmds[i][j]);
	}
		
	for(j=0;j<=mm;j++)
	{								//j节点控制区j+边界的热导率值
		lmdjp[0][j]=lmds[0][j];		lmdjp[nn+1][j]=lmds[nn][j];
		for(i=1;i<=nn;i++)
		lmdjp[i][j]=2.0*Dtx[i]/(dtx[i-1]/lmds[i-1][j]+dtx[i]/lmds[i][j]);
	}

//计算密度与比热容乘积，并取相邻区域平均值
	roc[0][0]=ros[0][0]*cs[0][0];			roc[nn+1][0]=ros[nn][0]*cs[nn][0];		//左上，右上
	roc[0][mm+1]=ros[0][mm]*cs[0][mm];		roc[nn+1][mm+1]=ros[nn][mm]*cs[nn][mm];	//左下，右下
	
	for(i=1;i<=nn;i++)																//上边
		roc[i][0]=(ros[i-1][0]*cs[i-1][0]*dtx[i-1]+ros[i][0]*cs[i][0]*dtx[i])/2.0/Dtx[i];
	for(i=1;i<=nn;i++)																//下边
		roc[i][mm+1]=(ros[i-1][mm]*cs[i-1][mm]*dtx[i-1]+ros[i][mm]*cs[i][mm]*dtx[i])/2.0/Dtx[i];
	for(j=1;j<=mm;j++)																//左边
		roc[0][j]=(ros[0][j-1]*cs[0][j-1]*dty[j-1]+ros[0][j]*cs[0][j]*dty[j])/2.0/Dty[j];
	for(j=1;j<=mm;j++)																//右边
		roc[nn+1][j]=(ros[nn][j-1]*cs[nn][j-1]*dty[j-1]+ros[nn][j]*cs[nn][j]*dty[j])/2.0/Dty[j];

	for(i=1;i<=nn;i++)
	{																				//中部
		for(j=1;j<=mm;j++)
		{
			roc[i][j]=ros[i-1][j-1]*cs[i-1][j-1]*dtx[i-1]*dty[j-1]+ros[i-1][j]*cs[i-1][j]*dtx[i-1]*dty[j];
			roc[i][j]=roc[i][j]+ros[i][j-1]*cs[i][j-1]*dtx[i]*dty[j-1]+ros[i][j]*cs[i][j]*dtx[i]*dty[j];
			roc[i][j]=roc[i][j]/4.0/Dtx[i]/Dty[j];
		}
	}

//一些变量赋给初值，或按多数取值先全部赋给初值，少数不同取值或随时间变化的变量在以后修改
	for(i=0;i<=nn+1;i++)for(j=0;j<=mm+1;j++){SS[i][j]=0.0;	TS[i][j]=15.0;}
													//内热源按多数取值先全部赋值0.0，节点温度赋给初值
	for(j=0;j<=mm+1;j++)								//左右边界条件
		{
			tfl[j]=TS[0][j];			alfl[j]=0.0;	//第二类边界条件(绝热)
			tfr[j]=TS[n+1][j];			alfr[j]=0.0;	//第二类边界条件(绝热)
		}
	for(i=0;i<=6;i++)									//上边界(北室外)，暂取值，在动态参数部分修改
		{
			tfu[i]=5.0;					alfu[i]=23.0;
		}
	for(i=7;i<=9;i++)									//上边界(墙下)，第二类边界条件(绝热)
		{
			tfu[i]=TS[i][0];			alfu[i]=0.0;
		}
	for(i=10;i<=25;i++)									//上边界(室内)，暂取值，在动态参数部分修改
		{
			tfu[i]=15.0;				alfu[i]=8.7;
		}
	for(i=26;i<=32;i++)									//上边界(南室外)，暂取值，在动态参数部分修改
		{
			tfu[i]=5.0;					alfu[i]=23.0;
		}
	for(i=0;i<=nn+1;i++)								//地下深处边界条件
		{
			tfd[i]=15.0;				alfd[i]=1000000.0;
		}

//计算二维导热差分方程组中各静态参数
	for(i=0;i<=nn;i++)for(j=0;j<=mm+1;j++)AA[i][j]=lmdip[i][j]*Dty[j]/dtx[i];
	for(i=0;i<=nn+1;i++)for(j=0;j<=mm;j++)BB[i][j]=lmdjp[i][j]*Dtx[i]/dty[j];
	for(i=0;i<=nn+1;i++)for(j=0;j<=mm+1;j++)KK[i][j]=roc[i][j]*Dtx[i]*Dty[j]/dtao;
	for(i=1;i<=nn;i++)for(j=1;j<=mm;j++)PP[i][j]=BB[i][j-1]+AA[i-1][j]+KK[i][j]+AA[i][j]+BB[i][j];

//计算边界条件式中各参数(四边及四角)
	for(i=1;i<=nn;i++)										
	{
		Pbi0[i]=alfu[i]*Dtx[i]+AA[i-1][0]+KK[i][0]+AA[i][0]+BB[i][0];				//上中
		Pbim[i]=BB[i][mm]+AA[i-1][mm+1]+KK[i][mm+1]+AA[i][mm+1]+alfd[i]*Dtx[i];		//下中
	}

	for(j=1;j<=mm;j++)
	{
		Pbj0[j]=BB[0][j-1]+alfl[j]*Dty[j]+KK[0][j]+AA[0][j]+BB[0][j];				//左中
		Pbjn[j]=BB[nn+1][j-1]+AA[nn][j]+KK[nn+1][j]+alfr[j]*Dty[j]+BB[nn+1][j];		//右中
	}

		Pc00=alfu[0]*Dtx[0]+alfl[0]*Dty[0]+KK[0][0]+AA[0][0]+BB[0][0];					//四角
		Pcn0=alfu[nn+1]*Dtx[nn+1]+AA[nn][0]+KK[nn+1][0]+alfr[0]*Dty[0]+BB[nn+1][0];
		Pcm0=BB[0][mm]+alfl[mm+1]*Dty[mm+1]+KK[0][mm+1]+AA[0][mm+1]+alfd[0]*Dtx[0];
		Pcnm=BB[nn+1][mm]+AA[nn][mm+1]+KK[nn+1][mm+1]+alfr[mm+1]*Dty[mm+1]+alfd[nn+1]*Dtx[nn+1];

//计算差分方程组矩阵中的各系数
	for(i=0;i<=nn+1;i++)for(j=0;j<=mm+1;j++){TT[j*(nn+2)+i]=TS[i][j];DD[j*(nn+2)+i]=0.0;}	//温度与常数向量

	for(i=0;i<(nn+2)*(mm+2);i++)	for(j=0;j<(nn+2)*(mm+2);j++)	RR[i][j]=0.0;			//先多数赋予0.0

		RR[0][0]=Pc00;			RR[0][1]=-AA[0][0];		RR[0][nn+2]=-BB[0][0];						//左上
	for(i=1;i<=nn;i++)
		{																							//上中
			RR[i][i-1]=-AA[i-1][0];	RR[i][i]=Pbi0[i];	RR[i][i+1]=-AA[i][0];	RR[i][nn+2+i]=-BB[i][0];
		}
			RR[nn+1][nn]=-AA[nn][0];RR[nn+1][nn+1]=Pcn0;	RR[nn+1][2*nn+3]=-BB[nn+1][0];			//右上
	for(j=1;j<=mm;j++)
		{																							//左中
			RR[j*(nn+2)][(j-1)*(nn+2)]=-BB[0][j-1];	RR[j*(nn+2)][j*(nn+2)]=Pbj0[j];
			RR[j*(nn+2)][j*(nn+2)+1]=-AA[0][j];		RR[j*(nn+2)][(j+1)*(nn+2)]=-BB[0][j];