📄 rgwshjhs.cpp
字号:
}
for(i=1;i<=nn;i++)for(j=1;j<=mm;j++)
{ //内部
RR[j*(nn+2)+i][(j-1)*(nn+2)+i]=-BB[i][j-1]; RR[j*(nn+2)+i][j*(nn+2)+i-1]=-AA[i-1][j];
RR[j*(nn+2)+i][j*(nn+2)+i]=PP[i][j];
RR[j*(nn+2)+i][j*(nn+2)+i+1]=-AA[i][j]; RR[j*(nn+2)+i][(j+1)*(nn+2)+i]=-BB[i][j];
}
for(j=1;j<=mm;j++)
{ //右中
RR[j*(nn+2)+nn+1][(j-1)*(nn+2)+nn+1]=-BB[nn+1][j-1];
RR[j*(nn+2)+nn+1][j*(nn+2)+nn]=-AA[nn][j]; RR[j*(nn+2)+nn+1][j*(nn+2)+nn+1]=Pbjn[j];
RR[j*(nn+2)+nn+1][(j+1)*(nn+2)+nn+1]=-BB[nn+1][j];
}
RR[(mm+1)*(nn+2)][mm*(nn+2)]=-BB[0][mm]; //左下
RR[(mm+1)*(nn+2)][(mm+1)*(nn+2)]=Pcm0; RR[(mm+1)*(nn+2)][(mm+1)*(nn+2)+1]=-AA[0][mm+1];
for(i=1;i<=nn;i++)
{ //下中
RR[(mm+1)*(nn+2)+i][mm*(nn+2)+i]=-BB[i][mm];
RR[(mm+1)*(nn+2)+i][(mm+1)*(nn+2)+i-1]=-AA[i-1][mm+1];
RR[(mm+1)*(nn+2)+i][(mm+1)*(nn+2)+i]=Pbim[i];
RR[(mm+1)*(nn+2)+i][(mm+1)*(nn+2)+i+1]=-AA[i][mm+1];
}
RR[(mm+1)*(nn+2)+nn+1][mm*(nn+2)+nn+1]=-BB[nn+1][mm]; //右下
RR[(mm+1)*(nn+2)+nn+1][(mm+1)*(nn+2)+nn]=-AA[nn][mm+1];
RR[(mm+1)*(nn+2)+nn+1][(mm+1)*(nn+2)+nn+1]=Pcnm;
}
//-------------------------------------------------------------------------------------------------------
//求解方程组的函数agsdl()(高斯-赛德尔(Gauss-Seidel)迭代法《C常用算法程序集》)
//本函数求得方程组的解，放在数组x[]中；
//调用本函数的实参依次为
//a[][] -- 方程组的系数矩阵，元素总数为n×n(n为虚参)，本程序n对应的实参为(nn+2)*(mm+2)，
// nn(=31),mm(=12)分别为水平方向与垂直方向内节点数，所以该数组的元素总数为462×462=213444；
//b[] -- 方程组右端的常数向量；
//n -- 方程组的阶数，即变量总数，本程序为所有节点的数量(nn+2)*(mm+2)=462；
//x[] -- 方程组的解向量，原程序的算法是设定所有元素初始值为0.0，然后经迭代求解，
// 本程序为加快运算速度，采用上一时刻的温度值作为迭代的初始值，可大大减少迭代次数；
//eps -- 迭代求解运算的精度(误差)要求。
//-------------------------------------------------------------------------------------------------------
int agsdl(double a[][462],double b[],int n,double x[],double eps)
{ int i,j;
double p,t,s,q;
for (i=0; i<=n-1; i++)
{ p=0.0;
// x[i]=0.0; //原程序执行x[i]=0.0是取初值全为0，改为不执行此句，x[i]中存放上一时刻的值作为初值
for (j=0; j<=n-1; j++)
if (i!=j)
{p=p+fabs(a[i][j]);}
if (p>=fabs(a[i][i]))
{ printf("fail\n"); return(-1);}
}
p=eps+1.0;
while (p>=eps)
{ p=0.0;
for (i=0; i<=n-1; i++)
{ t=x[i]; s=0.0;
for (j=0; j<=n-1; j++)
if (j!=i) s=s+a[i][j]*x[j];
x[i]=(b[i]-s)/a[i][i];
q=fabs(x[i]-t)/(1.0+fabs(x[i]));
if (q>p) p=q;
}
}
return(1);
}
//-------------------------------------------------------------------------------------------------------
//-------------------------------------------------------------------------------------------------------
//二维温度传递函数TT0_TT()
//由上一时刻温度分布计算下一时刻温度的函数，TT[]在调用函数之初为上一时刻温度，调用结束为下一时刻温度
//本函数执行将调用求解方程组的函数agsdl()
//调用本函数的实参依次为
//nn -- 计算域水平方向内节点数量，本程序按31个内节点(加上外节点共33个节点)编制；
//mm -- 计算域垂直方向内节点数量，本程序按12个内节点(加上外节点共14个节点)编制；
//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；
//TT[j] -- 离散方程矩阵中的温度向量，i,j节点数据逐行依次存放(元素个数=33*14=462)；
//SS[i][j] -- i,j节点代表的控制区域内热源(i=0～nn+1,j=0～mm+1)，W/m3；
//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)；
//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])。
//alfl[j]，alfr[j] -- 左，右边界对流换热系数(j=0～mm+1)，取等于0(第二类边界条件(绝热))，W/m2℃；
//alfu[i]，alfd[i] -- 上，下边界对流换热系数(i=0～nn+1)，W/m2℃；
//tfl[j]，tfr[j] -- 左，右边界处(j=0～mm+1)流体温度，取等于边界处土壤初始温度(以后保持不变)，℃；
//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℃；
//其他外部变量说明
//Pbi0[i] -- 边界条件式中各系数(上边，i=1～nn)；(可考虑改为局部变量)
//Pc00，Pcn0 -- 边界条件式中各系数(左上，右上)；(可考虑改为局部变量)
//局部变量说明
//DD -- 离散方程矩阵右端的常数向量；
//eps -- 迭代精度。
//-------------------------------------------------------------------------------------------------------
void TT0_TT(int nn,int mm,double Dtx[],double Dty[],double TT[],double SS[][14],
double AA[][14],double BB[][14],double KK[][14],double RR[][462],
double alfl[],double alfr[],double alfu[],double alfd[],
double tfl[],double tfr[],double tfu[],double tfd[])
{ int i,j; double eps;
static double DD[462];
eps=0.01; //迭代精度
//重新计算上边界条件式中各参数
Pc00=alfu[0]*Dtx[0]+alfl[0]*Dty[0]+KK[0][0]+AA[0][0]+BB[0][0]; //左上
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]; //上中
Pcn0=alfu[nn+1]*Dtx[nn+1]+AA[nn][0]+KK[nn+1][0]+alfr[0]*Dty[0]+BB[nn+1][0]; //右上
//重新计算差分方程组矩阵中与上边界有关的变化的各系数
RR[0][0]=Pc00; //左上
for(i=1;i<=nn;i++) RR[i][i]=Pbi0[i]; //上中
RR[nn+1][nn+1]=Pcn0; //右上
//构造常数向量矩阵
DD[0]=alfu[0]*Dtx[0]*tfu[0]+alfl[0]*Dty[0]*tfl[0]+KK[0][0]*TT[0]+SS[0][0]*Dtx[0]*Dty[0]; //左上
for(i=1;i<=nn;i++) DD[i]=KK[i][0]*TT[i]+alfu[i]*Dtx[i]*tfu[i]+SS[i][0]*Dtx[i]*Dty[0]; //上中
DD[nn+1]=alfu[nn+1]*Dtx[nn+1]*tfu[nn+1]+alfr[0]*Dty[0]*tfr[0]+KK[nn+1][0]*TT[nn+1]+SS[nn+1][0]*Dtx[nn+1]*Dty[0];
//右上
for(j=1;j<=mm;j++)
{
DD[j*(nn+2)]=KK[0][j]*TT[j*(nn+2)]+alfl[j]*Dty[j]*tfl[j]+SS[0][j]*Dtx[0]*Dty[j]; //左中
for(i=1;i<=nn;i++) DD[j*(nn+2)+i]=KK[i][j]*TT[j*(nn+2)+i]+SS[i][j]*Dtx[i]*Dty[j]; //内部
DD[j*(nn+2)+nn+1]=KK[nn+1][j]*TT[j*(nn+2)+nn+1]+alfr[j]*Dty[j]*tfr[j]+SS[nn+1][j]*Dtx[nn+1]*Dty[j];
//右中
}
DD[(mm+1)*(nn+2)]=alfl[mm+1]*Dty[mm+1]*tfl[mm+1]+alfd[0]*Dtx[0]*tfd[0];
DD[(mm+1)*(nn+2)]=DD[(mm+1)*(nn+2)]+KK[0][mm+1]*TT[(mm+1)*(nn+2)]+SS[0][mm+1]*Dtx[0]*Dty[mm+1]; //左下
for(i=1;i<=nn;i++)
DD[(mm+1)*(nn+2)+i]=KK[i][mm+1]*TT[(mm+1)*(nn+2)+i]+alfd[i]*Dtx[i]*tfd[i]+SS[i][mm+1]*Dtx[i]*Dty[mm+1];
//下中
DD[(mm+2)*(nn+2)-1]=alfd[nn+1]*Dtx[nn+1]*tfd[nn+1]+alfr[mm+1]*Dty[mm+1]*tfr[mm+1];
DD[(mm+2)*(nn+2)-1]=DD[(mm+2)*(nn+2)-1]+KK[nn+1][mm+1]*TT[(mm+2)*(nn+2)-1]+SS[nn+1][mm+1]*Dtx[nn+1]*Dty[mm+1];
//右下
agsdl(RR,DD,(nn+2)*(mm+2),TT,eps);
}
//-------------------------------------------------------------------------------------------------------
//以下为日光温室热量平衡分析部分的函数
//-------------------------------------------------------------------------------------------------------
//-------------------------------------------------------------------------------------------------------
//日光温室各部分面积与容积的计算函数AandV()-----------------------------------2007.10.1调试
//执行本函数的结果，计算得出温室各部分面积与容积，存入AAD[]数组(使用AAD[20]～AAD[29]段)；
//调用本函数的实参依次为
//kd -- 日光温室跨度(净跨度，即从墙内表面至屋面南部根基处)，m；
//Lgh -- 日光温室长度(净长度，即东西端墙内表面间距离)，m；
//Bgh[i]，Hgh[i] -- 日光温室横截面中屋面i(=0～17)点的横坐标与竖坐标(高度)，以内墙面根部为原点，m；
// 0点为北墙内墙面，1点为屋脊，2点及以后为前(南)屋面各点，
// 屋面南部根基处Bgh[i]=kd，Hgh[i]=0.00，跨度以外各点均为Bgh[i]=0.00，Hgh[i]=0.00；
//其他外部变量说明
//AAD[20] -- 后墙(北墙)面积，m2；
//AAD[21] -- 后屋面面积，m2；
//AAD[22] -- 日光温室横截面面积(也是东墙、西墙面积)，m2；
//AAD[23] -- 前(南)屋面面积，m2；
//AAD[24] -- 温室的室内地面面积，m2；
//AAD[25] -- 温室内部容积，m3；
//-------------------------------------------------------------------------------------------------------
void AandV(double kd,double Lgh,double Bgh[18],double Hgh[18])
{
int i;
AAD[20]=Lgh*Hgh[0]; //后墙(北墙)面积
AAD[21]=Lgh*sqrt((Bgh[1]-Bgh[0])*(Bgh[1]-Bgh[0])+(Hgh[1]-Hgh[0])*(Hgh[1]-Hgh[0])); //后屋面面积
AAD[22]=0.0; //横截面面积(东墙、西墙面积)
for(i=1;i<18;i++)AAD[22]=AAD[22]+(Bgh[i]-Bgh[i-1])*(Hgh[i]+Hgh[i-1])/2.0;
AAD[23]=0.0; //前(南)屋面面积
for(i=2;i<18;i++)
{
AAD[23]=AAD[23]+Lgh*sqrt((Bgh[i]-Bgh[i-1])*(Bgh[i]-Bgh[i-1])+(Hgh[i-1]-Hgh[i])*(Hgh[i-1]-Hgh[i]));
if(Hgh[i]<0.000001)break; //计算到前屋面根基处为止
}
AAD[24]=Lgh*kd; //温室的室内地面面积
AAD[25]=Lgh*AAD[22]; //温室内部容积
}
//-------------------------------------------------------------------------------------------------------
//-------------------------------------------------------------------------------------------------------
//日光温室前屋面透光率计算函数Tao_g()------------------------------------------2007.10.5调试
//本函数根据日光温室外保温是否覆盖以及透明覆盖材料透光性的情况，根据经验计算返回前屋面透光率
//
//调用本函数的实参依次为
//TGX -- 透明(固定)覆盖层在实际使用状态下的透光性，1透光性差，2透光性一般，3透光性好；
//FGZT -- 表示保温被或草帘卷放的参数，1卷起，2放下；
double Tao_g(int TGX,int FGZT)
{
if(FGZT==2) return 0.0; //保温被或草帘放下的情况
if(FGZT==1) //保温被或草帘卷起的情况
{
if(TGX==1) return 0.55;
if(TGX==2) return 0.65;
if(TGX==3) return 0.75;
}
}
//-------------------------------------------------------------------------------------------------------
//-------------------------------------------------------------------------------------------------------
//日光温室前屋面传热系数计算函数KG()------------------------------------------2007.10.5调试
//本函数根据日光温室覆盖材料和外保温是否覆盖的情况以及室外风速，根据经验计算返回前屋面传热系数W/(m2℃)
//调用本函数的实参依次为
//TFG -- 透明(固定)覆盖层的种类，1为PE，2为PVC，3为玻璃，4为中空PC板；
//WFG -- (活动)外覆盖保温被或草帘的保温性参数，1保温性较差，2保温性一般，3保温性较好，4保温性很好；
//FGZT -- 表示保温被或草帘是否覆盖的参数，1未覆盖，2覆盖；
//vo -- 室外风速，m/s。
//局部变量说明
//Kg -- 前(南)屋面传热系数，W/(m2℃)；
//-------------------------------------------------------------------------------------------------------
double KG(int TFG,int WFG,int FGZT,double vo)
{
double Kg;
if(FGZT==1) //保温被或草帘卷起的情况
{
if(TFG==1)Kg=6.6; if(TFG==2)Kg=6.4; if(TFG==3)Kg=6.2; if(TFG==4)Kg=3.5;
Kg=Kg*sqrt(0.6+0.1*vo); return Kg;
}
if(FGZT==2) //保温被或草帘放下的情况
{
if(WFG==1)Kg=2.5; if(WFG==2)Kg=2.0; if(WFG==3)Kg=1.8; if(WFG==4)Kg=1.6;
Kg=Kg*sqrt(0.8+0.05*vo); return Kg;
}
}
//-------------------------------------------------------------------------------------------------------
//-------------------------------------------------------------------------------------------------------
//日光温室后屋面传热系数计算函数KG1()----------------------------------------2007.10.5调试
//本函数返回后屋面传热系数W/(m2℃)
//-------------------------------------------------------------------------------------------------------
double KG1()
{
return 0.50;
}
//-------------------------------------------------------------------------------------------------------
//-------------------------------------------------------------------------------------------------------
//温室内外空气交换量计算函数LV()----------------------------------------2007.10.5调试
//本函数根据日光温室管理情况和室外风速，确定室内外空气的交换量m3/s
//调用本函数的实参依次为
//TF -- 表示通风强弱的变量，根据日光温室管理时通风口开启大小确定，
// TF=0，完全密闭，不通风，TF=1，弱通风，TF=2，一般通风，TF=3，较强通风，TF=4，强通风；
//MBX -- 表示温室密闭性好坏的参数，1密闭性较差，2密闭性一般，3密闭性较好，4密闭性很好；
//FGZT -- 表示保温被或草帘是否覆盖的参数，1未覆盖，2覆盖；
//Vgh -- 温室内部容积，m3；
//vo -- 室外风速，m/s。
//局部变量说明
//nv -- 换气次数，次/h。
//-------------------------------------------------------------------------------------------------------
double LV(int TF,int MBX,int FGZT,double Vgh,double vo)
{
double nv;
if(TF==0)
{
if(FGZT==1) //保温被或草帘卷起的情况
{
if(MBX==1)nv=1.0; if(MBX==2)nv=0.7; if(MBX==3)nv=0.5; if(MBX==4)nv=0.3;
}
if(FGZT==2) //保温被或草帘放下的情况
{
if(MBX==1)nv=0.5; if(MBX==2)nv=0.3; if(MBX==3)nv=0.2; if(MBX==4)nv=0.15;
}
}
if(TF==1)nv=2.0;
if(TF==2)nv=4.0;
if(TF==3)nv=8.0;
if(TF==4)nv=16.0;
return nv*Vgh*(0.6+0.1*vo)/3600.0;
}
//-------------------------------------------------------------------------------------------------------
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -