haymodeling.cpp

本程序的功能在于实现用Hay模型计算太阳能集热棚最佳倾斜角度

//////////////////////////////////////////////////////////////////////////////////////
//                                                                                  //
//          本程序的功能在于实现用Hay模型计算太阳能集热棚最佳倾斜角度               //
//                                                                                  //
//                         以昆明地区为例进行计算                                   //
//////////////////////////////////////////////////////////////////////////////////////
#include<stdio.h>
#include<math.h>
#define PI 3.1415926
#define CON 0.0174533  //PI/180的值；

void main()
{   
   	int i,k;
    int optmise;
	double qj=0;  //qj为集热棚倾斜角度;
	double p=0.2;  //定义p为地面反射率；
	double wd=25.02; //定义wd为当地纬度；
	double r=0; //倾斜面的方位角；
	double ws[12];  //定义ws为水平面上的日落时角；	
	double cw[12]={-20.84,-13.32,-2.4,9.46,18.78,23.04,21.11,13.28,1.97,-9.84,-19.02,-23.12};  //cw为赤纬角;	
	double N[12]={17,47,75,105,135,162,198,228,258,288,318,344};  //日子数
 	double H[12]={42622,46582,59242,64544,58053,47216,43077,45998,42850,38400,38225,37785};//昆明地区总辐射量(十年平均值)       单位为0.01兆焦耳/平方米
    double Hd[12]={11838,14251,20396,23342,27761,29862,30658,29360,25534,19409,14248,12713};//昆明地区散射辐射量(十年平均值)    单位为0.01兆焦耳/平方米
	double Ho[12];  //Ho为大气层外水平面上的太阳辐射量
	double Rb[12];  //Rb为太阳直射辐射修正因子；
	double Rd[12];  //Rd为太阳散射辐射修正因子；
	double HT[12];  //HT为倾斜集热棚上月平均一天的辐射量；
	double a[12],b[12],c[12],D[12];
    double wss[12],wsr[12];//倾斜集热棚面上的日出时角wsr和日落时角wss； 
    double Tot[901]; 
	double max;

    for(i=0;i<=900;i++)  //对集热棚倾斜角度从0~90度进行循环，步进为0.1度。
  
	{  
		double sum=0;

	    for(k=0;k<=11;k++)  //对一年12个月进行循环计算

		{	
			//计算水平面上的日落时角ws；
	        ws[k]=acos(-tan(wd*CON)*tan(cw[k]*CON))*180/PI;  

		//*************************************Bushell法计算集热棚倾斜面上的日出和日落时角****************************************************//
		//***************************该方法可以避免集热棚倾斜角度qj为0和集热棚方位角r为0时出现奇点********************************************//
			//计算中间变量a,b,c,D
	     	a[k]=sin(cw[k]*CON)*(sin(wd*CON)*cos(qj*CON)-cos(wd*CON)*sin(qj*CON)*cos(r*CON));                                         
        
		    b[k]=cos(cw[k]*CON)*(cos(wd*CON)*cos(qj*CON)+sin(wd*CON)*sin(qj*CON)*cos(r*CON));

    		c[k]=cos(cw[k]*CON)*sin(qj*CON)*sin(r*CON);

    		D[k]=sqrt(b[k]*b[k]+c[k]*c[k]);

            //计算倾斜集热棚面上的日出时角wsr和日落时角wss； 
    		wsr[k]=-(ws[k]<(fabs(-acos(-a[k]/D[k])+asin(c[k]/D[k]))*180/PI)? ws[k]:(fabs(-acos(-a[k]/D[k])+asin(c[k]/D[k]))*180/PI)); 

			wss[k]=ws[k]<((acos(-a[k]/D[k])+asin(c[k]/D[k]))*180/PI)? ws[k]:((acos(-a[k]/D[k])+asin(c[k]/D[k]))*180/PI);  		
    	//***********************************************************************************************************************************//

    		//计算太阳直射辐射修正因子Rb；       
         
    		Rb[k]=((cos(qj*CON)*sin(cw[k]*CON)*sin(wd*CON))*(wss[k]-wsr[k])*CON-(sin(cw[k]*CON)*cos(wd*CON)*sin(qj*CON)*cos(r*CON))*(wss[k]-wsr[k])*CON+(cos(wd*CON)*cos(cw[k]*CON)*cos(qj*CON))*(sin(wss[k]*CON)-sin(wsr[k]*CON))+(cos(cw[k]*CON)*cos(r*CON)*sin(wd*CON)*sin(qj*CON))*(sin(wss[k]*CON)-sin(wsr[k]*CON))-(cos(cw[k]*CON)*sin(qj*CON)*sin(r*CON))*(cos(wss[k]*CON)-cos(wsr[k]*CON)))/(2*(cos(wd*CON)*cos(cw[k]*CON)*sin(ws[k]*CON)+(CON)*ws[k]*sin(wd*CON)*sin(cw[k]*CON)));
		
			//计算大气层外水平面上的太阳辐射量Ho

			Ho[k]=24*3600*1367e-4*(1+0.033*cos(360*N[k]*CON/365)*(cos(wd*CON)*cos(cw[k]*CON)*sin(ws[k]*CON)+2*PI*ws[k]*sin(wd*CON)*sin(cw[k]*CON)/360))/PI;
	    
			//计算太阳散射辐射修正因子Rd；
		
    		Rd[k]=(H[k]-Hd[k])*Rb[k]/Ho[k]+(1-(H[k]-Hd[k])/Ho[k])*(1+cos(qj*CON))/2;

            //计算倾斜集热棚上月平均一天的辐射量HT；

			HT[k]=(H[k]-Hd[k])*Rb[k]+Hd[k]*Rd[k]+H[k]*p*(1-cos(qj*CON))/2;
   	
    		sum+=HT[k]; 

		}

 	Tot[i]=sum;

	qj+=0.1;

	}
/////////////////甄别出全年接收到最多的辐射量的集热棚倾斜角度/////////////////////////
   max=Tot[0];
    for(k=0;k<=901;k++)		
		if(Tot[k]>=max)
		{   
			max=Tot[k];	  
	        optmise=k;
		}   
////////////////////////////////////////////////////////////////////////////////////////
		
    printf("最佳的集热棚倾斜角度： %f deg\n与之对应的全年辐射量：%f 单位为0.01兆焦耳/平方米\n",(optmise*0.1),max);  
	
}