rgwshjmn.cpp

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

//日光温室热环境模拟程序
//主程序
//2007年10月


#include "stdio.h"
#include "math.h"
#include "rgwshjhs.cpp"
#include "rgwsbxhs.cpp"


//-------------------------------------------------------------------------------------------------------
//墙体传热模拟程序

//在主函数中出现的变量说明
//ee -- 程序预演的控制误差；
//OK -- 程序预演是否达到要求的标志，未达到要求OK=0，达到了要求，OK=1；
//Tn0 -- 模拟开始时墙体内表面温度(初始温度场是否充分发展的判断参数)，℃；
//h_total -- 模拟的总时间长度，h；
//i_tatal -- 模拟的总步数；
//dstep -- 24小时的模拟运行步数；
//nstep -- 控制结果输出的间隔步数(=间隔的秒数/dtao)；
//month0 -- 模拟开始月(继续预演运行的起点)；
//date0 -- 模拟开始日(继续预演运行的起点)；
//timeb0 -- 模拟开始时刻(继续预演运行的起点)；
//tao_g -- 温室屋面透光率；
//TFG -- 透明(固定)覆盖层的种类，1为PE，2为PVC，3为玻璃，4为中空PC板；
//TGX -- 透明(固定)覆盖层在实际使用状态下的透光性，1透光性差，2透光性一般，3透光性好；
//MBX --  表示温室密闭性好坏的参数，1密闭性较差，2密闭性一般，3密闭性较好，4密闭性很好；
//WFG -- (活动)外覆盖保温被或草帘的保温性参数，1保温性较差，2保温性一般，3保温性较好，4保温性很好；
//FGZT -- 表示保温被或草帘卷放的参数，1卷起，2盖下；
//t0,t1 -- 揭帘(保温被)的时间和盖帘(保温被)的时间，北京时间，时(分化为小时，如7时30分，为7.5时)；
//TF -- 表示通风强弱的变量，根据日光温室管理时通风口开启大小确定，
//				TF=0，完全密闭，不通风，TF=1，弱通风，TF=2，一般通风，TF=3，较强通风，TF=4，强通风；

//vo -- 室外风速，m/s。
//IDS -- 室外与墙面平行的垂直表面太阳总辐射照度，或墙体表面吸收的太阳辐射热量，W/m2；
//Ii -- 室内地面吸收的太阳辐射热量，W/m2；
//In -- 北室外地面吸收的太阳辐射热量，W/m2；
//Is -- 南室外地面吸收的太阳辐射热量，W/m2；
//ro_w -- 墙内表面对太阳辐射的反射率；
//ro_i -- 室内地面对太阳辐射的反射率；
//ro_n -- 北室外地面对太阳辐射的反射率；
//ro_s -- 南室外地面对太阳辐射的反射率；
//Kg -- 前(南)屋面传热系数，W/(m2℃)；
//Kg1 -- 后屋面传热系数，W/(m2℃)；
//Aw  (AAD[20]) -- 后墙(北墙)面积，m2；
//Ag1 (AAD[21]) -- 后屋面面积，m2；
//Ag2 (AAD[22]) -- 日光温室横截面面积(也是东墙、西墙面积)，m2；
//Ag3 (AAD[23]) -- 前(南)屋面面积，m2；
//As  (AAD[24]) -- 温室的室内地面面积，m2；
//Vgh (AAD[25]) -- 温室内部容积，m3；
//Lv -- 温室内外空气交换量，m3/s；
//Qg1 -- 后屋面传出温室的热量，W
//Qg -- 前屋面传出温室的热量，W
//Qw -- 墙体传入温室内的热量，W
//Qs -- 地面传入温室内的热量，W
//Qv -- 空气交换传出温室的显热量，W
//Qe -- 室内蒸腾蒸发潜热量，W
//QQ -- 室内热量总收支，得热为正，失热为负，W

//p -- 当地大气压力，Pa，默认值为101325.0Pa；
//tfn1 --室内空气温度(上一时刻)，℃；
//tfn -- 室内空气温度，℃；
//tf0 -- 室外空气温度，℃；
//fin1 --室内空气相对湿度(上一时刻)，(小数)；
//fin -- 室内空气相对湿度，(小数)；
//fi0 -- 室外空气相对湿度，(小数)；
//pw1  --室内空气水蒸汽压力(上一时刻)，Pa；
//pw  -- 室内空气水蒸汽压力，Pa；
//pw0 -- 室外空气水蒸汽压力，Pa；
//pws1 -- 室内空气饱和水蒸汽压力(上一时刻)，Pa；
//pws  -- 室内空气饱和水蒸汽压力，Pa；
//pws0 -- 室外空气饱和水蒸汽压力，Pa；
//dn1 --室内空气含湿量(上一时刻)，g/kg；
//dn -- 室内空气含湿量，g/kg；
//d0 -- 室外空气含湿量，g/kg；
//Gw -- 室内地面和植物蒸发蒸腾量，g/m2；
//Gws -- 室内地面蒸发量，g/m2；
//Gd -- 室内水凝结量，g/m2；
//Gv -- 逸出温室的水汽量，g/m2；
//tgn -- 前屋面内表面温度，℃；
//pgs -- 与tgn对应的饱和水蒸汽压力，Pa；
//k0 -- 计算水蒸腾蒸发的温度系数；
//k1 -- 地面状况系数，1.05土壤地面，无覆盖，1.0土壤地面，覆盖地膜或砖铺地面，0.95混凝土地面；
//k2 -- 地面潮湿系数，1.1极潮湿，积水多，1.05很潮湿，少量积水，1.0较潮湿，0.95较干燥，0.9很干燥；
//k3 -- 植物繁茂系数，0.8稀疏，0.9较稀疏，1.0一般，1.1较茂密，1.2很茂密；
//Kv -- 蒸腾蒸发系数，g/(m2sPa)；
//Kd -- 凝结系数，g/(m2sPa)；

//---------------------------------------------------------------------------------------------
//以下程序中，后面带有“//$$”符号的行，为根据不同计算问题需修改数据的语句
main()
{
	int i,j,i0,j0,i_tatal,month0,date0,dstep,nstep,OK;
	int TFG,TGX,MBX,WFG,FGZT,TF;
	double vo,h_total,IDS,Ii,In,Is,tao_g,ro_w,ro_i,ro_n,ro_s,timeb0,t0,t1,Tn0,ee;
	double Lgh,Vgh,Aw,Ag1,Ag2,Ag3,As,Lv,Kg,Kg1,Qg1,Qg,Qw,Qs,Qv,Qe,QQ;

	double p,fin,fi0,pw,pw0,pws,pws0,dn,d0,Gw,Gws,Gd,Gv,tgn,pgs,k0,k1,k2,k3,Kv,Kd;
	double tfn1,dn1,pw1,pws1,fin1;

//静态数据输入和参数计算------------------------------------------------------

//总体静态数据输入-----------------------------------
	dtao=60.0;				h_total=360.0;			//$$输入时间步长，模拟的总时数
	ee=0.005;				nstep=3600.0/dtao;		//$$设定程序预演误差和结果输出间隔(=间隔秒数/dtao)

	month=12;		date=21;		timeb=0.0;		//$$输入模拟初始月、日、(北京)时间

	bw=40.0;		dj=118.4;						//$$输入当地经纬度
	kd=8.0;			Lgh=50.0;						//$$输入温室跨度、长度
	TFG=1;				//$$确定透明(固定)覆盖层的种类，1--PE，2--PVC，3--玻璃，4--中空PC板
	TGX=2;				//$$确定透明(固定)覆盖层在实际使用状态下的透光性，1透光差，2透光一般，3透光好
	MBX=2;				//$$确定温室密闭性好坏，1密闭性较差，2密闭性一般，3密闭性较好，4密闭性很好
	WFG=2;				//$$保温被或草帘的保温性，1保温性较差，2保温性一般，3保温性较好，4保温性很好

	p=101325.0;			//$$输入当地大气压力

//以下输入日光温室横截面中屋面的横坐标Bgh[](m)与竖坐标(高度)Hgh[](m)，以内墙面根部为原点，
//0点的Bgh[0]，Hgh[0]为内墙面横坐标与高度，1点的Bgh[1]，Hgh[1]为屋脊横坐标与高度，
//2点及以后为前(南)屋面各点(最多给出16点)，依次给出横坐标与高度，
//横坐标等于跨度的点，Bgh[]=跨度，Hgh[]=0.00，以后各点均取为Bgh[]=0.00，Hgh[]=0.00
//注意：折线形屋面需给出屋面各组成平面相交处坐标，
//拱形屋面至少给出3点(不包括屋脊与前屋面根基处)，且其中至少包括前屋面根基以内0.3～1.5m范围内的一点。
	Bgh[0]=0.00;	Bgh[1]=1.50;	Bgh[2]=2.50;	Bgh[3]=3.50;	Bgh[4]=4.50;	Bgh[5]=5.50;	//$$
	Hgh[0]=2.50;	Hgh[1]=3.30;	Hgh[2]=3.08;	Hgh[3]=2.82;	Hgh[4]=2.48;	Hgh[5]=2.06;	//$$

	Bgh[6]=6.50;	Bgh[7]=7.50;	Bgh[8]=8.00;	Bgh[9]=0.00;	Bgh[10]=0.00;	Bgh[11]=0.00;	//$$
	Hgh[6]=1.49;	Hgh[7]=0.67;	Hgh[8]=0.00;	Hgh[9]=0.00;	Hgh[10]=0.00;	Hgh[11]=0.00;	//$$

	Bgh[12]=0.00;	Bgh[13]=0.00;	Bgh[14]=0.00;	Bgh[15]=0.00;	Bgh[16]=0.00;	Bgh[17]=0.00;	//$$
	Hgh[12]=0.00;	Hgh[13]=0.00;	Hgh[14]=0.00;	Hgh[15]=0.00;	Hgh[16]=0.00;	Hgh[17]=0.00;	//$$

	AandV(kd,Lgh,Bgh,Hgh);												//求算温室各部分面积与容积
	Aw=AAD[20];		Ag1=AAD[21];	Ag2=AAD[22];	Ag3=AAD[23];	As=AAD[24];		Vgh=AAD[25];

//墙体静态数据输入和参数计算-------------------------
	bmfw=0.0;	bmqx=90.0;		ro_w=0.10;				//$$墙体方位，墙面反射率

	nw=3;												//$$输入墙体材料层数
	dw0[0]=0.240;	dw0[1]=0.060;	dw0[2]=0.240;		//$$输入墙体各层材料厚度(从外至内)
	lmd0[0]=0.81;	lmd0[1]=0.035;	lmd0[2]=0.81;		//$$输入墙体各层材料导热系数(从外至内)
	ro0[0]=1800.0;	ro0[1]=20.0;	ro0[2]=1800.0;		//$$输入墙体各层材料密度(从外至内)
	c0[0]=880.0;	c0[1]=2000.0;	c0[2]=880.0;		//$$输入墙体各层材料比热容(从外至内)

	n=jdhf(nw,dw0);										//划分墙体节点，计算内节点数
	ywdrcfjc(n,dx,dw0,lmd0,ro0,c0,dtao);	//设定初始温度，计算墙体一维导热差分方程组静态参数

//地下静态数据输入和参数计算-------------------------
	nn=31;			mm=12;								//$$输入地下网格划分横向内节点和竖向内节点数量
	ro_i=0.05;		ro_n=0.10;		ro_s=0.10;			//$$室内、外地面对太阳辐射的反射率
//输入地下各层材料导热系数
//室内外土壤分四层次，分界处深度0.1m、0.3m、0.8，从上往下，导热系数、密度增大，比热容不变
//在上部同等深度，室外土壤导热系数、密度大于或等于室内，但比热容不变
//墙下0.8m以上为砖砌体，各热工参数与墙体相同，1.3m以下各热工参数与室内外土壤相同
//		室外			墙下			室内						分层
	lmdw[0] =1.000;	lmdx[0] =0.81;	lmdn[0] =0.600;		//$$深度0.000m～0.002m处
	lmdw[1] =1.000;	lmdx[1] =0.81;	lmdn[1] =0.600;		//$$深度0.002m～0.006m处
	lmdw[2] =1.000;	lmdx[2] =0.81;	lmdn[2] =0.600;		//$$深度0.006m～0.014m处
	lmdw[3] =1.000;	lmdx[3] =0.81;	lmdn[3] =0.600;		//$$深度0.014m～0.030m处	上部
	lmdw[4] =1.000;	lmdx[4] =0.81;	lmdn[4] =0.600;		//$$深度0.030m～0.055m处
	lmdw[5] =1.000;	lmdx[5] =0.81;	lmdn[5] =0.600;		//$$深度0.055m～0.100m处____________
	lmdw[6] =1.200;	lmdx[6] =0.81;	lmdn[6] =1.000;		//$$深度0.100m～0.180m处	中上部
	lmdw[7] =1.200;	lmdx[7] =0.81;	lmdn[7] =1.000;		//$$深度0.180m～0.300m处____________
	lmdw[8] =1.200;	lmdx[8] =0.81;	lmdn[8] =1.200;		//$$深度0.300m～0.500m处	中下部
	lmdw[9] =1.200;	lmdx[9] =0.81;	lmdn[9] =1.200;		//$$深度0.500m～0.800m处____________
	lmdw[10]=1.400;	lmdx[10]=1.400;	lmdn[10]=1.400;		//$$深度0.800m～1.300m处
	lmdw[11]=1.400;	lmdx[11]=1.400;	lmdn[11]=1.400;		//$$深度1.300m～2.100m处	下部
	lmdw[12]=1.400;	lmdx[12]=1.400;	lmdn[12]=1.400;		//$$深度2.100m～3.300m处
//输入地下各层材料密度
	row[0]=1200.0;	rox[0]=1800.0;	ron[0]=1000.0;		//$$深度0.000m～0.002m处
	row[1]=1200.0;	rox[1]=1800.0;	ron[1]=1000.0;		//$$深度0.002m～0.006m处
	row[2]=1200.0;	rox[2]=1800.0;	ron[2]=1000.0;		//$$深度0.006m～0.014m处
	row[3]=1200.0;	rox[3]=1800.0;	ron[3]=1000.0;		//$$深度0.014m～0.030m处	上部
	row[4]=1200.0;	rox[4]=1800.0;	ron[4]=1000.0;		//$$深度0.030m～0.055m处
	row[5]=1200.0;	rox[5]=1800.0;	ron[5]=1000.0;		//$$深度0.055m～0.100m处____________
	row[6]=1500.0;	rox[6]=1800.0;	ron[6]=1200.0;		//$$深度0.100m～0.180m处	中上部
	row[7]=1500.0;	rox[7]=1800.0;	ron[7]=1200.0;		//$$深度0.180m～0.300m处____________
	row[8]=1500.0;	rox[8]=1800.0;	ron[8]=1500.0;		//$$深度0.300m～0.500m处	中下部
	row[9]=1500.0;	rox[9]=1800.0;	ron[9]=1500.0;		//$$深度0.500m～0.800m处____________
	row[10]=1800.0;	rox[10]=1800.0;	ron[10]=1800.0;		//$$深度0.800m～1.300m处
	row[11]=1800.0;	rox[11]=1800.0;	ron[11]=1800.0;		//$$深度1.300m～2.100m处	下部
	row[12]=1800.0;	rox[12]=1800.0;	ron[12]=1800.0;		//$$深度2.100m～3.300m处
//输入地下各层材料比热容
	cw[0]=1400.0;	cx[0]=880.0;	cn[0]=1400.0;		//$$深度0.000m～0.002m处
	cw[1]=1400.0;	cx[1]=880.0;	cn[1]=1400.0;		//$$深度0.002m～0.006m处
	cw[2]=1400.0;	cx[2]=880.0;	cn[2]=1400.0;		//$$深度0.006m～0.014m处
	cw[3]=1400.0;	cx[3]=880.0;	cn[3]=1400.0;		//$$深度0.014m～0.030m处	上部
	cw[4]=1400.0;	cx[4]=880.0;	cn[4]=1400.0;		//$$深度0.030m～0.055m处
	cw[5]=1400.0;	cx[5]=880.0;	cn[5]=1400.0;		//$$深度0.055m～0.100m处____________
	cw[6]=1400.0;	cx[6]=880.0;	cn[6]=1400.0;		//$$深度0.100m～0.180m处	中上部
	cw[7]=1400.0;	cx[7]=880.0;	cn[7]=1400.0;		//$$深度0.180m～0.300m处____________
	cw[8]=1400.0;	cx[8]=880.0;	cn[8]=1400.0;		//$$深度0.300m～0.500m处	中下部
	cw[9]=1400.0;	cx[9]=880.0;	cn[9]=1400.0;		//$$深度0.500m～0.800m处____________