📄 rgwshjmn.cpp
字号:
cw[10]=1400.0; cx[10]=1400.0; cn[10]=1400.0; //$$深度0.800m～1.300m处
cw[11]=1400.0; cx[11]=1400.0; cn[11]=1400.0; //$$深度1.300m～2.100m处 下部
cw[12]=1400.0; cx[12]=1400.0; cn[12]=1400.0; //$$深度2.100m～3.300m处
wghf(kd,nw,dw0); //地下部分网格划分
ewdrcfjc(nn,mm,dtx,dty,lmdw,lmdx,lmdn,row,rox,ron,cw,cx,cn,dtao);
//设定初始温度，计算地下二维导热差分方程组静态参数
//输入与蒸发蒸腾有关的静态参数(本部分为暂用程序，需进一步修改)
k1=1.0; //$$地面状况系数，1.05土壤地面，无覆盖，1.0土壤地面，覆盖地膜或砖铺地面，0.95混凝土地面
k2=1.0; //$$地面潮湿系数，1.1极潮湿，积水多，1.05很潮湿，少量积水，1.0较潮湿，0.95较干燥，0.9很干燥
k3=1.0; //$$植物繁茂系数，0.8稀疏，0.9较稀疏，1.0一般，1.1较茂密，1.2很茂密
//逐步递推流程----------------------------------------------------------------
i=0; OK=0;
month0=month; date0=date; timeb0=timeb; //记下模拟开始时间(继续预演运行的起点)
dstep=3600.0*24.0/dtao; i_tatal=h_total*3600.0/dtao;
tfn=15.0; fin=0.80; //室内空气状态初始假设值
pws=610.78*exp(17.2694*tfn/(tfn+237.3));
pw=fin*pws; dn=622.0*pw/(p-pw);
while(i<=i_tatal)
{
//总体动态数据---------------------------------------
if(i==0)Tn0=T[n+1]; //记下程序预演初始时刻墙体内表面温度
j=datecount(month,date,timeb,dtao); //月、日和时间递推
month=AAI[1]; date=AAI[2]; timeb=AAD[16];
tf0=tf_0(i*dtao/3600.0); //$$确定室外一侧的气温,i*dtao/3600.0为时数
fi0=fi_0(i*dtao/3600.0); //$$确定室外一侧的相对湿度(小数)
TF=0; //$$确定是否通风，0为密闭不通风，1～4通风弱～强
IDS=TYFZ(bw,dj,month,date,timeb,bmfw,bmqx); //太阳辐射照度计算
vo=v_o(); //室外风速
t0=t_0(month,date); t1=t_1(month,date); //揭帘(保温被)的时间t0，盖帘(保温被)的时间t1
FGZT=2; if(timeb>t0 && timeb<t1)FGZT=1; if(AAD[9]+AAD[11]<30.0)FGZT=2;
//$$根据预设的时间判断外保温揭开还是覆盖，但如室外辐射过弱(界限值可设定)，仍保持不揭开的状态
Kg=KG(TFG,WFG,FGZT,vo); Kg1=KG1(); //前、后屋面的传热系数
tao_g=Tao_g(TGX,FGZT); //确定透光率
Lv=LV(TF,MBX,FGZT,Vgh,vo); //室内外空气交换量
//墙体动态数据---------------------------------------
IDS=tao_g*(1.0-ro_w)*IDS; //室内墙面吸收的太阳辐射热量
alf0=alf_0(); alfn=alf_n(); //确定墙体外侧、内侧表面换热系数
S[0]=(1.0-ro_w)*0.5*AAD[11]/Dx[0];
//墙北面吸收的太阳辐射(水平面散射辐射的1/2)热流量转成边节点控制体积中的内热源
S[n+1]=IDS/Dx[n+1]; //墙南表面热流量转成边节点控制体积中的内热源
T0_T(n,alf0,alfn,tf0,tfn,Dx,T,S,A,K,P); //递推下一时刻墙体温度分布
//水汽收支量-------------------------------(本部分为暂用程序，需进一步修改)
//室内地面和植物蒸发蒸腾量
k0=15.0+0.05*tfn*tfn; Kv=1e-6*k0*k1*k2*k3;
Gw=Kv*(pws-pw); Gws=(1.0-0.8*k3)*Gw;
//凝结水量
Kd=3.0e-6; tgn=tfn-Kg*(tfn-tf0)/alfn;
pgs=610.78*exp(17.2694*tgn/(tgn+237.3));
Gd=Kd*(pw-pgs); if(Gd<0.0)Gd=0.0;
//逸出温室的水汽量
pws0=610.78*exp(17.2694*tf0/(tf0+237.3)); pw0=fi0*pws0;
d0=622.0*pw0/(p-pw0); Gv=Lv*1.2*(dn-d0)/As;
//地下动态数据---------------------------------------
Ii=tao_g*(1.0-ro_i)*(AAD[9]+AAD[11])-Gws*2442.0; //室内地面吸收的太阳辐射热量
In=0.9*(1.0-ro_n)*AAD[11]; //北室外地面吸收的太阳辐射热量，考虑地面蒸发带走热量，折减10%
Is=0.9*(1.0-ro_s)*(AAD[9]+AAD[11]); //南室外地面吸收的太阳辐射热量，考虑地面蒸发带走热量，折减10%
for(i0=0;i0<=6;i0++) //上边界(北室外)
{
tfu[i0]=tf0; alfu[i0]=alf_0(); SS[i0][0]=In/Dty[0];
}
for(i0=10;i0<=25;i0++) //上边界(室内)
{
tfu[i0]=tfn; alfu[i0]=alf_n(); SS[i0][0]=Ii/Dty[0];
}
for(i0=26;i0<=32;i0++) //上边界(南室外)
{
tfu[i0]=tf0; alfu[i0]=alf_0(); SS[i0][0]=Is/Dty[0];
}
TT0_TT(nn,mm,Dtx,Dty,TT,SS,AA,BB,KK,RR,alfl,alfr,alfu,alfd,tfl,tfr,tfu,tfd); //递推下一时刻地下温度分布
//计算各部分热量收支项-------------------------------
Qg1=Kg1*Ag1*(tfn-tf0); //后屋面传出的热量
Qg=Kg*Ag3*(tfn-tf0); //前屋面传出的热量
Qw=Aw*alfn*(T[n+1]-tfn); //墙体传入的热量
Qs=0.0;
for(i0=11;i0<=25;i0++)Qs=Qs+Dtx[i0]*alfu[i0]*(TT[i0]-tfn);Qs=Lgh*Qs;//地面传入的热量
Qv=Lv*1.2*1030*(tfn-tf0); //空气交换传出的显热量
Qe=2442.0*(Gw-Gd)*As; //室内蒸腾蒸发潜热量
QQ=Qw+Qs-Qg1-Qg-Qv-Qe; //室内热量总收支，得热正，失热负
//下一时刻室内空气状态(先保存上一时刻状态值，仅用于输出数据)
tfn1=tfn; tfn=tfn+QQ*dtao/1.2/1030/Vgh;
dn1=dn; dn=dn+As*dtao*(Gw-Gd-Gv)/1.2/Vgh;
pw1=pw; pw=p*dn/(622+dn);
pws1=pws; pws=610.78*exp(17.2694*tfn/(tfn+237.3));
fin1=fin; fin=pw/pws;
if(pw>pws) //出现该情况说明下一时刻气温有可能下降至露点温度以下
{ //则下一时刻含湿量及水蒸汽压维持不变，气温为露点温度
pws=pw; fin=1.0; tfn=237.3/(17.2694/(log(pw)-log(610.78))-1.0);
}
//程序进程控制---------------------------------------
if(i==dstep) //程序预演是否可以结束的判断：
{ //以24小时周期前后墙体内表面温度变化作为判断指标
if(fabs(Tn0-T[n+1])>ee) //如果变化大于设定误差ee，返回初始条件继续预演
{
i=0; month=month0; date=date0; timeb=timeb0; continue;
}
else OK=1;
}
i=i+1;
if(i>i_tatal) break;
if(fmod(i,2*nstep)==0.0) //调试的中间结果输出
printf("%3.0f时 室内气温%6.3f 相对湿度%6.4f 墙外表%6.3f 墙内表%6.3f 地表%6.3f\n",timeb,tfn1,fin1,T[0],T[n+1],TT[16]);
//数据输出
if(fmod(i,nstep)==0.0 && OK==1) //输出温室参数、状态与管理数据，每隔nstep步输出一次结果数据
{
printf("\n步长%5.1fs 误差%8.5f℃\n\n",dtao,ee);
printf("覆盖材料 ");if(TFG==1)printf(" PE ");if(TFG==2)printf("PVC ");if(TFG==3)printf("玻璃");if(TFG==4)printf(" PC ");
printf(" 透光性 ");if(TGX==1)printf("差");if(TGX==2)printf("一般");if(TGX==3)printf("好");
printf(" 密闭性 ");if(MBX==1)printf("较差");if(MBX==2)printf("一般");if(MBX==3)printf("较好");if(MBX==4)printf("很好");
printf(" 外覆盖保温性 ");if(WFG==1)printf("较差");if(WFG==2)printf("一般");if(WFG==3)printf("较好");if(WFG==4)printf("很好");
printf("\n");
printf("温室参数 跨度%5.1f 长度%6.1f 脊高%5.1f 后墙高%5.1f 容积%6.1f\n",kd,Lgh,Hgh[1],Hgh[0],Vgh);
printf("各部面积 后墙%6.1f 后屋面%6.1f 横截面%5.1f 前屋面%6.1f 地面%6.1f\n",Aw,Ag1,Ag2,Ag3,As);
printf("\n外覆盖卷放状态 ");if(FGZT==1)printf("卷起");if(FGZT==2)printf("盖下");
printf(" 通风强弱 ");if(TF==0)printf("密闭");if(TF==1)printf("弱");if(TF==2)printf("一般");if(TF==3)printf("较强");if(TF==4)printf("强");
printf(" 揭帘时间%5.2f 盖帘时间%5.2f\n\n",t0,t1);
}
if(fmod(i,nstep)==0.0 && OK==1) //输出太阳总辐射照度计算数据，每隔nstep步输出一次结果数据
{
printf("%2d月%2d日 北纬%6.2f° 东经%6.2f°表面方位角%6.2f°表面倾斜角%6.2f°\n",month,date,bw,dj,bmfw,bmqx);
printf("%5.2f时 时角%6.2f°太阳赤纬%6.2f°高度%6.2f°方位%6.2f°大气透明度%5.2f\n",timeb,AAD[1],AAD[2],AAD[3]/3.1416*180.0,AAD[4]/3.1416*180.0,AAD[5]);
printf("太阳I0=%5.0f IDN=%6.2f ID=%6.2f IDH=%6.2f IS=%6.2f ISH=%6.2f I=%6.2f\n\n",AAD[6],AAD[7],AAD[8],AAD[9],AAD[10],AAD[11],IDS);
}
if(fmod(i,nstep)==0.0 && OK==1) //输出温室内外热交换量，每隔nstep步输出一次结果数据
{
printf("传热系数 后屋面%6.3f 前屋面%6.3f 透光率%6.3f 室外风速%6.3f 换气量%6.3f\n",Kg1,Kg,tao_g,vo,Lv);
printf("传热 后屋面%8.1f 前屋面%8.1f 墙体%8.1f 地面%8.1f\n",Qg1,Qg,Qw,Qs);
printf(" 换气显热%8.1f 潜热%8.1f 总热量%8.1f\n\n",Qv,Qe,QQ);
}
if(fmod(i,nstep)==0.0 && OK==1) //输出温室内水汽变化量，每隔nstep步输出一次结果数据
{
printf("k0%7.3f k1%5.2f k2%5.2f k3%5.2f Kv%10.8f Kd%10.8f tgn%6.4f pgs%5.1f\n",k0,k1,k2,k3,Kv,Kd,tgn,pgs);
printf("蒸发蒸腾%8.6f 地面%8.6f 凝结%8.6f 逸出%8.6f 总增加%8.6f\n\n",Gw,Gws,Gd,Gv,Gw-Gd-Gv);
}
if(fmod(i,nstep)==0.0 && OK==1) //输出下一时刻空气状态，每隔nstep步输出一次结果数据
{
printf("室外 气温%6.3f 相对湿度%6.4f 水汽压力%6.1f 饱和压力%6.1f 含湿量%6.3f\n",tf0,fi0,pw0,pws0,d0);
printf("室内始 气温%6.3f 相对湿度%6.4f 水汽压力%6.1f 饱和压力%6.1f 含湿量%6.3f\n",tfn1,fin1,pw1,pws1,dn1);
printf("室内终 气温%6.3f 相对湿度%6.4f 水汽压力%6.1f 饱和压力%6.1f 含湿量%6.3f\n",tfn,fin,pw,pws,dn);
scanf("%c",&i0);
}
if(fmod(i,nstep)==0.0 && OK==1) //输出墙体温度分布数据，每隔nstep步输出一次结果数据
{
printf(" 外表面传出热量%8.3fW/m2 内表面传出热量%8.3fW/m2\n",alf0*(T[0]-tf0),alfn*(T[n+1]-tfn1));
printf("%4.2f时 tf0%6.2f S0%8.1f alf0%5.2f alfn%5.2f Sn%8.1f tfn%6.2f\n",timeb,tf0,S[0],alf0,alfn,S[n+1],tfn1);
for(j0=0;j0<=n/10;j0++)
{
printf("点编号"); for(i0=10*j0;i0<10*(j0+1);i0++)printf("%6d",i0); printf("\n");
printf("温度℃ "); for(i0=10*j0;i0<10*(j0+1);i0++)printf("%6.2f",T[i0]); printf("\n");
}
scanf("%c",&i0);
}
if(fmod(i,nstep)==0.0 && OK==1) //输出地下温度分布数据，每隔nstep步输出一次结果数据
{
printf(" %4.2f时 室外北 室内北 室内中 室内南 室外南\n",timeb);
printf("地面上气温℃ %8.3f %8.3f %8.3f %8.3f %8.3f\n",tfu[3],tfu[14],tfu[18],tfu[22],tfu[29]);
printf("换热系数W/m2℃ %8.3f %8.3f %8.3f %8.3f %8.3f\n",alfu[3],alfu[14],alfu[18],alfu[22],alfu[29]);
printf("地面热源 %8.2fW/m3 %8.2fW/m3 %8.2fW/m3 %8.2fW/m3 %8.2fW/m3\n",SS[3][0],SS[14][0],SS[18][0],SS[22][0],SS[29][0]);
printf("地面传热 %8.2fW/m2 %8.2fW/m2 %8.2fW/m2 %8.2fW/m2 %8.2fW/m2\n",
alfu[2]*(TT[3]-tfu[3]),alfu[14]*(TT[14]-tfu[14]),alfu[18]*(TT[18]-tfu[18]),alfu[22]*(TT[22]-tfu[22]),alfu[29]*(TT[29]-tfu[29]));
printf("室外北地下、墙下温度℃\n");
for(j0=0;j0<=mm+1;j0++) {for(i0=0;i0<=11;i0++)printf("%6.2f",TT[j0*(nn+2)+i0]);printf("\n");}
printf("室内地下温度℃\n");
for(j0=0;j0<=mm+1;j0++) {for(i0=12;i0<=24;i0++)printf("%6.2f",TT[j0*(nn+2)+i0]);printf("\n");}
printf("室外南地下温度℃\n");
for(j0=0;j0<=mm+1;j0++) {for(i0=25;i0<=32;i0++)printf("%6.2f",TT[j0*(nn+2)+i0]);printf("\n");}
scanf("%c",&i0);
}
}
}
//-------------------------------------------------------------------------------------------------------
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -