Subroutine Type30 !*********************************************************************************************** ! Collector array shading - this component will determine the incident radiation upon an array ! of collectors that shade one another. ! ! Revision History: ! 1993.05.00 - JWT ! 2004.02.00 - TPM: Updated to TRNSYS 16 coding standards. ! 2005.08.29 - DEB: Started adding PAR checks to avoid a divide by zero. ! 2005.09.28 - DAA Major changes to Mode2 - Parabollic collector array: ! 1. Changed number of possible outputs from 20 to 2. ! 2. Changed number of parameters to 5 ! 3. 'ModeAx' parameter is not needed anymore. ! 4. RCHECK was not being called. ! 2005.06.05 - JWT/DAA : Added new outputs to mode 1 ! 2009.02.18 - MJD/DEB : fixed the INFO(6) setting for mode 1. ! 2009.06.00 - TPM : conversion to version 17 coding standards ! !----------------------------------------------------------------------------------------------- ! Parabollic collector array !----------------------------------------------------------------------------------------------- ! ! PARAMETERS: ! Mode = 1 for flat plate collectors, 2 for parabollic collectors ! ! If Mode = 2 - Parabollic collectors: ! ModeAx - KEPT FOR BACKWARDS COMPATIBILITY - NOT USED SINCE VERSION 16.1 ! D_a - Distance between collector axes [m] ! W_a - Collector amplitude [m] ! N_r, - Number of rows ! Beta_ap - Slope of axes plane [deg] ! ! INPUTS ! I_bT - Incident beam on tilted surface [kJ/h-m^2] ! beta - Collector slope [deg] ! ! INTERMEDIATES ! theta_p - Intermediate angle [deg] ! P - Distance between axes in a plane normal to beam radiation [m] ! f_bs_1 - Shading fraction for one collector row [-] ! ! OUTPUTS ! I_bT_s - Shaded incident beam radiation [kJ/h-m^2] ! f_bs - Shading fraction for entire collector field [-] ! ! !*********************************************************************************************** ! Copyright © 2011 Solar Energy Laboratory, University of Wisconsin-Madison. All rights reserved. !export this subroutine for its use in external DLLs. !DEC$ATTRIBUTES DLLEXPORT :: TYPE30 !----------------------------------------------------------------------------------------------------------------------- Use TrnsysConstants Use TrnsysFunctions !----------------------------------------------------------------------------------------------------------------------- !Variable Declarations Implicit None !force explicit declaration of local variables Double Precision Time,Timestep Double Precision rad,hcl,widcl,slcl,sepcl,azcl,slfld,ht1,hb1,hd1,zen,azsol,hthor,hdhor,hbslfl,alb,usvfgr, & usvfsk,sep,d,theta1,theta6,alf,xl5,xl7,xl6,bb,xl8,svfgr1,svfgr2,svgr12,svfsk,svfgr,arvfsk,arvfgr,hh,ss,tsl, & delaz,sin1,zenaz,alt,theta2,hb,hs,xl2,xl1,xl3,a,as,st1,st2,fclb,fhb,hb1s,na15,da15,angl15,angl14,angl1,angl2, & angl3,grs,grus,cc,vfb1,angl10,angl11,angl12,angl13,ee,s,dfbbc,vfd1,vfd2,dfdbc,hgrbcl,hd1s,ht1s,cos1 Double Precision D_a !Distance between axes Double Precision W_a !Collector amplitude Double Precision N_r !Number of rows Double Precision Beta_ap !Slope of axes plane Double Precision I_bT !Incident beam radiation on tilted surface Double Precision beta !Collector slope Double Precision theta_p !Intermediate angle to calculate shading Double Precision P !Distance between axes in a plane normal to beam radiation Double Precision f_bs_1 !Shading fraction for one collector row Double Precision f_bs !Shading fraction of entire collector field Double Precision I_bT_s !Shaded incident radiation Integer mode,j,xnrow,modeax,CurrentUnit,CurrentType Data rad/0.01745329/ !----------------------------------------------------------------------------------------------------------------------- !Get the Global Trnsys Simulation Variables Time=getSimulationTime() Timestep=getSimulationTimeStep() CurrentUnit = getCurrentUnit() CurrentType = getCurrentType() !----------------------------------------------------------------------------------------------------------------------- !----------------------------------------------------------------------------------------------------------------------- !Set the Version Number for This Type If(getIsVersionSigningTime()) Then Call SetTypeVersion(17) Return Endif !----------------------------------------------------------------------------------------------------------------------- !----------------------------------------------------------------------------------------------------------------------- !Do All of the Last Call Manipulations Here If(getIsLastCallofSimulation()) Then Return Endif !----------------------------------------------------------------------------------------------------------------------- !----------------------------------------------------------------------------------------------------------------------- !Perform Any "End of Timestep" Manipulations That May Be Required If(getIsEndOfTimestep()) Then Return Endif !----------------------------------------------------------------------------------------------------------------------- !----------------------------------------------------------------------------------------------------------------------- !Do All of the "Very First Call of the Simulation Manipulations" Here If(getIsFirstCallofSimulation()) Then mode = JFIX(getParameterValue(1)+0.1) If (mode == 1) Then !Flat plate collector Call SetNumberofParameters(8) Call SetNumberofInputs(8) Call SetNumberofDerivatives(0) Call SetNumberofOutputs(5) Call SetInputUnits(1,'IR1') Call SetInputUnits(2,'IR1') Call SetInputUnits(3,'DG1') Call SetInputUnits(4,'DG1') Call SetInputUnits(5,'IR1') Call SetInputUnits(6,'IR1') Call SetInputUnits(7,'IR1') Call SetInputUnits(8,'DM1') Do j = 1,5 Call SetOutputUnits(j,'IR1') EndDo Else !Parabolic collector array Call SetNumberofParameters(6) Call SetNumberofInputs(2) Call SetNumberofDerivatives(0) Call SetNumberofOutputs(5) Call SetInputUnits(1,'IR1') Call SetInputUnits(2,'DG1') Call SetOutputUnits(1,'IR1') Call SetOutputUnits(2,'DM1') Call SetOutputUnits(3,'NAV') Call SetOutputUnits(4,'NAV') Call SetOutputUnits(5,'NAV') EndIf Return EndIf !----------------------------------------------------------------------------------------------------------------------- !----------------------------------------------------------------------------------------------------------------------- !Do All of the "Start Time" Manipulations Here - There Are No Iterations at the Intial Time If (getIsStartTime()) Then ! Read in the Values of the Parameters from the Input File mode = JFIX(getParameterValue(1)+0.1) If (mode == 1) Then !Flat plate collector hcl = getParameterValue(2) widcl = getParameterValue(3) slcl = getParameterValue(4)*rad sepcl = getParameterValue(5) xnrow = getParameterValue(6) azcl = getParameterValue(7) * rad slfld = getParameterValue(8) * rad ! Check the Parameters for Problems If (hcl <= 0.d0) Call FoundBadParameter(2,'Fatal','The heat exchanger mode must be between 1 and 7.') If ((slfld < 0.0) .or. (slfld > slcl) .or. (slfld >= 1.50)) Call FoundBadParameter(8,'Fatal','The field slope is out of range.') If (slcl < 0.0 .or. slcl > (90.*rad-0.00001)) Call FoundBadParameter(4,'Fatal','The collector slope is out of range.') If (ErrorFound()) Return ! Calculate collector view factors of sky and ground ! Unshaded collectors usvfgr=(1.-DCOS(slcl))/2.0 usvfsk=(1.+DCOS(slcl))/2.0 ! Shaded collectors sep=sepcl/DCOS(slfld) d=(sep**2.+hcl**2.-2.*sep*hcl*DCOS(slcl-slfld))**0.5 theta1=DCOS((sep**2.+d**2.-hcl**2.)/(2.*sep*d)) If (sep < 0.001) theta1 = 180.0*rad-slcl theta6 = theta1 - slfld alf = slcl-slfld xl5 = sepcl xl7 = sep xl6 = sepcl - (hcl * DCOS(slcl)) bb = sepcl * DTAN(slfld) - hcl * DSIN(slcl) xl8 = (bb**2. + xl6**2.)**0.5 If (bb <= 0.0) xl6 = xl8 svfgr1 = ((xl6-xl5) + (xl7-xl8)) / (2.*hcl) bb = (hcl**2.+sep**2.-(2.*hcl*sep*DCOS(180.*rad-alf)))**.5 svfgr2 = (hcl+sep-bb)/(2.*hcl) svgr12 = svfgr1 + svfgr2 svfsk = (((hcl + sep) - d) / (2.*hcl)) - svfgr1 svfgr = svfgr1 ! Average for array arvfsk=(svfsk*(xnrow-1.0)+usvfsk)/xnrow arvfgr=(svfgr*(xnrow-1.0)+usvfgr)/xnrow ! Calculate collector array angles hh=hcl*DCOS(slcl) ss=sepcl-hh tsl=DABS(slfld-slcl) ! Error insufficient collector separation If (ss < 0.) Then Call Messages(-1,'Insufficient Collector Row Separation','Fatal',CurrentUnit,CurrentType) Return EndIf ! Initialize outputs Do j = 1,5 Call SetOutputValue(j,0.d0) EndDo Else !Parabolic collector array ! Read in the Values of the Parameters from the Input File modeax = getParameterValue(2) !NOT USED - KEPT FOR BACKWARDS COMPATIBILITY D_a = getParameterValue(3) !Distance between axes W_a = getParameterValue(4) !Collector amplitude N_r = getParameterValue(5) !Number of rows Beta_ap = getParameterValue(6) * rad !Slope of axes plane ! Initialize outputs Call SetOutputValue(1,getInputValue(1)) !Unshaded radiation Call SetOutputValue(2,0.d0) !Shaded fraction EndIf Return Endif !----------------------------------------------------------------------------------------------------------------------- !----------------------------------------------------------------------------------------------------------------------- !ReRead the Parameters if Another Unit of This Type Has Been Called Last If(getIsReReadParameters()) Then mode = JFIX(getParameterValue(1)+0.1) If (mode == 1) Then !Flat plate collector hcl = getParameterValue(2) widcl = getParameterValue(3) slcl = getParameterValue(4)*rad sepcl = getParameterValue(5) xnrow = getParameterValue(6) azcl = getParameterValue(7) * rad slfld = getParameterValue(8) * rad ! Calculate collector view factors of sky and ground ! Unshaded collectors usvfgr=(1.-DCOS(slcl))/2.0 usvfsk=(1.+DCOS(slcl))/2.0 ! Shaded collectors sep=sepcl/DCOS(slfld) d=(sep**2.+hcl**2.-2.*sep*hcl*DCOS(slcl-slfld))**0.5 theta1=DCOS((sep**2.+d**2.-hcl**2.)/(2.*sep*d)) If (sep < 0.001) theta1 = 180.0*rad-slcl theta6 = theta1 - slfld alf = slcl-slfld xl5 = sepcl xl7 = sep xl6 = sepcl - (hcl * DCOS(slcl)) bb = sepcl * DTAN(slfld) - hcl * DSIN(slcl) xl8 = (bb**2. + xl6**2.)**0.5 If (bb <= 0.0) xl6 = xl8 svfgr1 = ((xl6-xl5) + (xl7-xl8)) / (2.*hcl) bb = (hcl**2.+sep**2.-(2.*hcl*sep*DCOS(180.*rad-alf)))**.5 svfgr2 = (hcl+sep-bb)/(2.*hcl) svgr12 = svfgr1 + svfgr2 svfsk = (((hcl + sep) - d) / (2.*hcl)) - svfgr1 svfgr = svfgr1 ! Average for array arvfsk=(svfsk*(xnrow-1.0)+usvfsk)/xnrow arvfgr=(svfgr*(xnrow-1.0)+usvfgr)/xnrow ! Calculate collector array angles hh=hcl*DCOS(slcl) ss=sepcl-hh tsl=DABS(slfld-slcl) Else !Parabolic collector array ! Read in the Values of the Parameters from the Input File modeax = getParameterValue(2) !NOT USED - KEPT FOR BACKWARDS COMPATIBILITY D_a = getParameterValue(3) !Distance between axes W_a = getParameterValue(4) !Collector amplitude N_r = getParameterValue(5) !Number of rows Beta_ap = getParameterValue(6) * rad !Slope of axes plane EndIf Endif !----------------------------------------------------------------------------------------------------------------------- !----------------------------------------------------------------------------------------- ! Flat plate collector array !----------------------------------------------------------------------------------------- If (mode == 1) Then ! Retrieve the current values of the inputs ht1 = getInputValue(1) hb1 = getInputValue(2) hd1 = ht1 - hb1 zen = getInputValue(3) * rad azsol = getInputValue(4) * rad hthor = getInputValue(5) hdhor = getInputValue(6) hbslfl = getInputValue(7) alb = getInputValue(8) ! Calculate incident angle in plane of collector azimuth If (ht1 <= 0.0) Then ! No solar Call SetOutputValue(1,getInputValue(1)) Call SetOutputValue(2,getInputValue(2)) Call SetOutputValue(3,getInputValue(1)-getInputValue(2)) Call SetOutputValue(4,0.d0) Call SetOutputValue(5,0.d0) Else If (zen >= 1.570) zen = 1.570 delaz = DABS(azsol - azcl) sin1 = DCOS(delaz) * DSIN(zen) cos1 = DCOS(zen) If (cos1 < 0.01 .and. sin1 >= 0.0) Then zenaz = 90.0*rad ElseIf (cos1 < 0.01 .and. sin1 < 0.0) Then zenaz = -90.0*rad Else zenaz = DATAN(sin1/cos1) EndIf alt = 90.0*rad - zenaz theta6 = theta1 - slfld If (alt <= 180.0*rad-slcl) Then If (alt >= theta6) goto 210 ! Calculate array beam irradiated fraction If (slfld+alt <= 0.0) Then hb=0.0 Else theta2 = rad*180.0 - alf - (slfld + alt) hb = sep * ((DSIN(alt+slfld)) / (DSIN(theta2))) If (hb > hcl) goto 210 EndIf hs = hcl - hb xl2 = DABS(ss*DTAN(delaz)) If (xl2 >= widcl) goto 210 ! Check for slcl = 90 degrees If (slcl >= 89.99*rad) Then a = hh*widcl as = (widcl-xl2)*xl1 fclb=1.0-as/a Else xl1 = hs * DCOS(slcl) xl3 = DABS((ss+xl1)*DTAN(delaz)) ! Calculate shaded collector fraction a = hh*widcl If (xl3 >= widcl) Then st1 = widcl-xl2 st2 = DABS(st1*DTAN(90.0*rad-delaz)) AS = (st1*st2)/2.0 Else AS=(widcl-((xl3+xl2)/2.0))*xl1 EndIf fclb=1.0-as/a EndIf goto 220 ! No beam shading 210 fclb = 1.0 220 fhb = (1.0 + fclb * (xnrow - 1.0)) / xnrow Else fhb = 0.0 EndIf hb1s=hb1*fhb ! Calculate beam radiation reflected to collectors from ground between collector rows If (alt < theta6) Then grus = 0.0 Else na15 = (DSIN(slcl))*hcl + (DTAN(slfld))*sepcl da15 = (DCOS(slcl))*hcl + sepcl angl15 = 180.*rad-DATAN(na15/da15) If (alt > angl15) goto 274 angl14 = 180.0*rad - slcl If (alt > angl14) goto 273 ! Shading of ground between collector rows from beam radiation when alt < 180 - slcl angl1 = 180.0*rad - 90.0*rad - slcl If (alt > 90.0*rad + angl1) Then grus = sep Else angl2 = angl1 + zenaz angl3 = 180.0*rad - angl2 - alf grs = hcl*(DSIN(angl2))/(DSIN(angl3)) grus = sep - grs EndIf EndIf ! Calculate view factor hcl to grus cc = (hcl**2.0+grus**2.0-2.0*hcl*grus*DCOS(180.*rad-alf))**.5 vfb1 = (hcl+grus-cc)/(2.0*hcl) goto 275 ! Shading of ground between collector rows from beam radiation when alt > 180-slcl 273 angl10 = 180.*rad-90.*rad-slcl angl11 = DABS(zenaz)-angl10 angl12 = 180.*rad-slcl+slfld angl13 = 180.*rad-angl11-angl12 ee = hcl*DSIN(angl12)/DSIN(angl13) s = hcl*DSIN(angl11)/DSIN(angl13) VFB1 = ((sep+ee)-(bb+s))/(2.0*hcl) If (s > sep) vfb1=0.0 If (vfb1 < 0.0) vfb1=0.0 goto 275 274 vfb1 = 0.0 275 dfbbc = (vfb1*hbslfl*alb*(xnrow-1.0))/xnrow ! Calculate diffuse radiation reflected to collector surface from ground between collectors ! Calculate view factor sep to sky vfd1 = (bb+d-(2.0*hcl))/(2.0*sep) ! Calculate view factor hcl to sep vfd2 = (hcl+sep-bb)/(2.0*hcl) dfdbc = (vfd1*vfd2*hdhor*alb*(xnrow-1.0))/xnrow hgrbcl = dfbbc + dfdbc ! Calculate array diffuse irradiated fraction hd1s = hdhor*arvfsk + hthor*arvfgr*alb + hgrbcl ! Calculate ht on shaded array ht1s=hd1s+hb1s ! Set the outputs Call SetOutputValue(1,ht1s) Call SetOutputValue(2,hb1s) Call SetOutputValue(3,hd1s) Call SetOutputValue(4,hdhor*arvfsk) !sky diffuse Call SetOutputValue(5,hthor*arvfgr*alb + hgrbcl) !grounf reflected diffuse EndIf Else !---------------------------------------------------------------- ! Parabollic collector array !---------------------------------------------------------------- ! Set inputs for parabollic trough mode I_bT = getInputValue(1) !Incident beam radiation beta = getInputValue(2) * rad !Collector slope ! Calculate shading ! Geometry for two collector rows theta_p = beta - beta_ap !Intermediate angle P = D_a * DCOS(theta_p) !Distance between axes in a plane normal to beam radiation ! Shading fraction for one collector row f_bs_1 = DMAX1(0.d0,1.0-P/W_a) ! Shading fraction of entire collector field f_bs = ((N_r - 1.d0 )*f_bs_1 )/N_r ! Effectice incident radiation on surface ! Parabollic collectors: only beam radiation I_bT_s = (1.d0 - f_bs)* I_bT !Shaded incident radiation ! Set outputs Call SetOutputValue(1,I_bT_s) !Shaded incident radiation Call SetOutputValue(2,f_bs) !Shading fraction EndIf Return End