rgwsbxhs.cpp

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

//日光温室热环境模拟程序的函数
//本文件为备选函数，有些函数根据具体问题需改写
//2007年9月


int month,date;
double bw,dj,timeb,bmfw,bmqx,IDS;


//-------------------------------------------------------------------------------------------------------
//日期递进计算函数datecount()----------------------------------------------------------2007.9.4调试
//本函数根据起始日期和时间，计算经过一个步长的时间(s,3600s以内)后的月、日、时
//时间累计尚未进入下一日时，返回0；进入下一日但尚未跨月时，返回1；跨月时，返回2.

//全局变量说明
//AAI[30] -- 用于将整型局部变量的计算结果从函数内传送到外部的数组，各函数通用，
//                但不同函数分配使用不同的区段，日期递进计算函数使用AAI[0]～AAI[4]的区段
//AAD[80] -- 用于将双精度型局部变量的计算结果从函数内传送到外部的数组，各函数通用，
//                但不同函数分配使用不同的区段，日期递进计算函数使用AAD[15]～AAD[19]的区段

//调用本函数的实参依次为           函数调用后，下一步的日期与时间放在传送数组中对应元素中↓
//month -- 月  -----------------------------------------------------------------------AAI[1]
//date -- 日   -----------------------------------------------------------------------AAI[2]
//timeb -- 时(北京时间)，形式为xx.xx时(24时制)，分以下单位需化成时   -----------------AAD[16]
//dtao -- 时间步长，s；
//---------------------------------------------------------------------------------------------
int datecount(int month,int date,double timeb,double dtao)
{
	AAI[1]=month;	AAI[2]=date;	AAD[16]=timeb;

	AAD[16]=timeb+dtao/3600.0;
	if(AAD[16]<24.0)return(0);								//时间累计尚未进入下一日

	AAI[2]=AAI[2]+1;	AAD[16]=AAD[16]-24.0;

	if(AAI[2]>28 && AAI[1]==2)								//2月跨到3月，2月按平年28日计算
	{
		AAI[1]=AAI[1]+1;	AAI[2]=1;	return(2);
	}		

	if(AAI[2]>30)											//小月的月末跨月的情况
		if(AAI[1]==4 || AAI[1]==6 || AAI[1]==9 || AAI[1]==11)
		{
			AAI[1]=AAI[1]+1;	AAI[2]=1;	return(2);
		}

	if(AAI[2]>31)											//大月的月末跨月的情况
		if(AAI[1]==1 || AAI[1]==3 || AAI[1]==5 || AAI[1]==7 || AAI[1]==8 || AAI[1]==10 || AAI[1]==12)
		{
			AAI[1]=AAI[1]+1;	if(AAI[1]>12)AAI[1]=AAI[1]-12;	AAI[2]=1;	return(2);
		}
	return(1);												//跨日不能跨月的情况
}
//-------------------------------------------------------------------------------------------------------


//-------------------------------------------------------------------------------------------------------
//大气透明度计算函数dqtmd()
//执行本函数的结果返回大气透明度值

//调用本函数的实参依次为
//bw -- 北纬度数，度
//month -- 月份
//---------------------------------------------------------------------------------------------
double dqtmd(double bw,int month)
{	double x;
	x=month;
	if(bw<27.5)				return (0.0021*x*x - 0.0254*x + 0.7005);
	if(bw>27.5 && bw<=32.5)	return (0.0038*x*x - 0.0471*x + 0.7707);
	if(bw>32.5 && bw<=37.5)	return (0.0038*x*x - 0.0459*x + 0.7691);
	if(bw>37.5 && bw<=42.5)	return (0.0033*x*x - 0.0409*x + 0.772);
	if(bw>42.5)				return (0.0039*x*x - 0.0481*x + 0.787);
}	  
//-------------------------------------------------------------------------------------------------------


//-------------------------------------------------------------------------------------------------------
//计算太阳辐射照度的函数TYFZ()---------------------------------------------------------2007.9.6调试
//执行本函数的结果返回指定倾角表面的太阳总辐射照度值IDS，计算出水平面的直接辐射IDH和散射辐射ISH，W/m2，
//并将其以及一些中间计算结果存入数组AAD[]，以便调试检查和使用

//本函数执行将调用大气透明度计算函数dqtmd()

//调用本函数的实参依次为
//bw -- 北纬度数，度
//dj -- 东经度数，度
//month,date -- 月、日
//timeb -- 时(北京时间)，形式为xx.xx时(24时制)，分以下单位需化成时
//bmfw -- 表面方位角，度
//bmqx -- 表面倾斜角，度

//全局变量说明
//AAD[80] -- 用于将双精度型局部变量的计算结果从函数内传送到外部的数组，各函数通用，
//                但不同函数分配使用不同的区段，太阳辐射照度相关函数使用AAD[0]～AAD[14]的区段
//IDS -- 指定倾角表面的太阳总辐射照度，W/m2  (是函数返回值，但未在本函数中出现)

//局部变量说明                                                      传送数组中对应元素↓
//time -- 当地平均太阳时，形式为xx.xx时(24时制)，分以下单位需化成时    ------------AAD[0]
//sj -- 时间角，度    -------------------------------------------------------------AAD[1]
//tycw -- 太阳赤纬角，度    -------------------------------------------------------AAD[2]
//tygd -- 太阳高度角，弧度    -----------------------------------------------------AAD[3]
//tyfw -- 太阳方位角，弧度    -----------------------------------------------------AAD[4]
//p -- 大气透明度    --------------------------------------------------------------AAD[5]
//I0 -- 大气层外边界处法向太阳辐射照度，W/m2    -----------------------------------AAD[6]
//IDN -- 法向太阳直接辐射，W/m2    ------------------------------------------------AAD[7]
//ID -- 指定位置表面上的太阳直接辐射，W/m2    -------------------------------------AAD[8]
//IDH -- 水平面太阳直接辐射    ----------------------------------------------------AAD[9]
//IS -- 指定位置表面上的太阳散射辐射，W/m2    -------------------------------------AAD[10]
//ISH -- 水平面太阳散射辐射(实际上为天空散射，仅计算天空散射辐射)，W/m2    --------AAD[11]
//d -- 1月1日至所计算日之前天数之和
//m -- 大气质量
//sinh -- 太阳高度角的正弦
//cosbt -- 太阳光线与表面法线夹角的余弦
//temp -- 中间变量

//---------------------------------------------------------------------------------------------
double TYFZ(double bw,double dj,int month,int date,double timeb,double bmfw,double bmqx)
{
	double time,sj,tycw,tygd,tyfw,p,I0,IDN,ID,IDH,IS,ISH;
	double d,dp,m,sinh,cosbt,temp;

	if	(month==1) {I0=1405.0;	d=date;}
	else if	(month==2) {I0=1394.0;	d=31+date;}
	else if	(month==3) {I0=1378.0;	d=59+date;}
	else if	(month==4) {I0=1353.0;	d=90+date;}
	else if	(month==5) {I0=1334.0;	d=120+date;}
	else if	(month==6) {I0=1316.0;	d=151+date;}
	else if	(month==7) {I0=1308.0;	d=181+date;}
	else if	(month==8) {I0=1315.0;	d=212+date;}
	else if	(month==9) {I0=1330.0;	d=243+date;}
	else if	(month==10){I0=1350.0;	d=273+date;}
	else if	(month==11){I0=1372.0;	d=304+date;}
	else if	(month==12){I0=1392.0;	d=334+date;}

	dp=3.1416/180.0;

	tycw=23.45*cos((360.0*(d-172.0)/365.0)*dp);
	time=timeb-(120.0-dj)/15.0;		sj=15.0*(time-12.0);
	sinh=cos(bw*dp)*cos(tycw*dp)*cos(sj*dp)+sin(bw*dp)*sin(tycw*dp);
	tygd=asin(sinh);

	//太阳方位角的计算，多数书籍上说法不准确，见《建筑热过程》86年版上的红字批注
	temp=(sin(tygd)*sin(bw*dp)-sin(tycw*dp))/cos(tygd)/cos(bw*dp);
	if(temp<=-1.0)tyfw=180.0;	//由于计算机的数字存储方式产生的误差，可能出现temp绝对值略大于1的情况
	if(temp>=1.0)tyfw=0.0;
	if(temp>-1.0 && temp<1.0)tyfw=acos(temp);
	if(sj<0)tyfw=(-1.0)*tyfw;	//太阳方位角正负与时角相同，acos()只能计算出正角，负角只能根据时角判断

	cosbt=sin(tygd)*cos(bmqx*dp)+cos(tygd)*sin(bmqx*dp)*cos(tyfw-bmfw*dp);

	p=dqtmd(bw,month);
	m=1.0/sin(tygd);

	if(sinh<=0.0) {IDN=0.0;	ID=0.0; IDH=0.0; IS=0.0; ISH=0.0;}

	if(sinh>0.0)
	{
		IDN=I0*pow(p,m);
		IDH=IDN*sin(tygd);
		ID=IDN*cosbt;	if(cosbt<0.0)ID=0.0;		//cosbt<0.0，有bt>90度，即光线照到该表面背面去了

		ISH=I0/2.0*(1.0-pow(p,m))/(1.0-1.4*log(p))*sin(tygd);
		IS=ISH*cos(bmqx*dp/2.0)*cos(bmqx*dp/2.0);
	}

//以下将局部变量的计算结果放入传送数组中，以便以后在函数外部使用或调试中输出检查
	AAD[0]=time;
	AAD[1]=sj;	AAD[2]=tycw;	AAD[3]=tygd;	AAD[4]=tyfw;	AAD[5]=p;
	AAD[6]=I0;	AAD[7]=IDN;		AAD[8]=ID;		AAD[9]=IDH;		AAD[10]=IS;		AAD[11]=ISH;

	return (ID+IS);
}
//-------------------------------------------------------------------------------------------------------


//-------------------------------------------------------------------------------------------------------
//以下为内外侧流体温度、相对湿度、热源以及对流换热系数
//---------------------------------------------------------------------------------------------
double tf_0(double h)		//h为24时制的时数
{
	double pi,tp,a1,a2,f1,f2,tmax,tmin,t;
	pi=3.1416;	tp=-0.016;	a1=0.48;	a2=0.08;	f1=3.8;		f2=0.005;
	tmax=5.0;		tmin=-10.0;
	t=(tmax+tmin)/2.0+(tmax-tmin)*(tp+a1*cos(2.0*h*pi/24.0-f1)+a2*cos(4.0*h*pi/24.0-f2));
	return (t);
}
//---------------------------------------------------------------------------------------------
double fi_0(double h)		//h为24时制的时数(fi为小数)
{
	double pi,fip,a1,a2,f1,f2,fimax,fimin,fi;
	pi=3.1416;	fip=-0.016;	a1=0.48;	a2=0.08;	f1=3.8;		f2=0.005;
	fimax=0.90;		fimin=0.40;
	fi=(fimax+fimin)/2.0-(fimax-fimin)*(fip+a1*cos(2.0*h*pi/24.0-f1)+a2*cos(4.0*h*pi/24.0-f2));
	return (fi);
}
//---------------------------------------------------------------------------------------------
double tf_n(double h)		//h为24时制的时数
{
	double pi,tp,a1,a2,f1,f2,tmax,tmin,t;
	pi=3.1416;	tp=-0.016;	a1=0.48;	a2=0.08;	f1=3.8;		f2=0.005;
	tmax=25.0;		tmin=10.0;
	t=(tmax+tmin)/2.0+(tmax-tmin)*(tp+a1*cos(2.0*h*pi/24.0-f1)+a2*cos(4.0*h*pi/24.0-f2));
	return (t);
}
//---------------------------------------------------------------------------------------------
double alf_0()
{
return(23.0);
}
//---------------------------------------------------------------------------------------------
double alf_n()
{
return(8.7);
}
//---------------------------------------------------------------------------------------------


//---------------------------------------------------------------------------------------------
//确定风速的函数v_o()

//---------------------------------------------------------------------------------------------
double v_o()
{
return(4.0);
}
//---------------------------------------------------------------------------------------------


//---------------------------------------------------------------------------------------------
//确定揭帘(保温被)时间的函数t_0()
//本函数根据日光温室的管理制度，给出某月份和日期揭帘(保温被)的时间
//调用本函数的实参依次为
//month,date -- 月、日

//---------------------------------------------------------------------------------------------
double t_0(int month,int date)
{
	if(month==10)	return 6.5;
	if(month==11)	return 7.0;
	if(month==12)	return 7.5;
	if(month==1)	return 7.5;
	if(month==2)	return 7.0;
	if(month==3)	return 6.5;
					return 6.0;
}
//---------------------------------------------------------------------------------------------

//---------------------------------------------------------------------------------------------
//确定盖帘(保温被)时间的函数t_1()
//本函数根据日光温室的管理制度，给出某月份和日期盖帘(保温被)的时间
//调用本函数的实参依次为
//month,date -- 月、日

//---------------------------------------------------------------------------------------------
double t_1(int month,int date)
{
	if(month==10)	return 18.0;
	if(month==11)	return 17.5;
	if(month==12)	return 17.0;
	if(month==1)	return 17.0;
	if(month==2)	return 17.5;
	if(month==3)	return 18.0;
					return 19.0;
}
//---------------------------------------------------------------------------------------------