! This component has been assigned Type Number 205. If that number conflicts with ! another user Type number, you will need to change it and recompile the appropriate ! dll. SUBROUTINE TYPE205(TIME,XIN,OUT,T,DTDT,PAR,INFO,ICNTRL,*) !DEC$ATTRIBUTES DLLEXPORT :: TYPE205 C************************************************************************ C THIS SUBROUTINE MODELS A SERPENTINE COLLECTOR. C A SERPENTINE COLLECTOR CAN BE MODELED USING THE CONVENTIONAL C HEADER-RISER PARALLEL FLOW COLLECTOR MODEL, IF THE NUMBER OF C TURNS ARE GREATER THAN ABOUT FIFTEEN. THE ONLY DIFFERENCE IS C THE CALCULATION OF THE HEAT TRANSFER COEFFICIENT, WHICH IS C CALCULATED FOR A LONG TUBE. C LAST MODIFIED 25 AUGUST 1997 -- MYRNA DAYAN C************************************************************************ DOUBLE PRECISION XIN,OUT INTEGER*4 INFO DIMENSION XIN(14),OUT(8),PAR(13),INFO(15) CHARACTER*3 YCHECK(14),OCHECK(8) INTEGER N,N_G,ITER REAL k,delta,D,D_i,T_a,L,W,C_p,T_in,T_pm,T_out REAL TauAlpha, L_tube_serpentine,F_1, RefInd REAL m_flat,U_L,F,beta,Q_useful,A_c,mu,Pr,kw REAL F_R,m_dot, h_fi_serp,pi,wind,U_be,Alpha REAL Nusselt_serp, Re_Serp, T_pm_old, change, E_p REAL I_D, I_T, I_H, XKL,theta,P_loss, rho, rho_g ! Set the version information for TRNSYS IF (INFO(7).EQ.-2) THEN INFO(12) = 15 RETURN 1 ENDIF IF (INFO(7).GE.0) GO TO 100 C FIRST CALL OF SIMULATION INFO(6)=8 INFO(9)=0 CALL TYPECK(1,INFO,14,13,0) DATA YCHECK/'TE1','MF1','TE1','IR1','VE1','IR1','IR1' .,'DM1','DG1','DG1','VS1','CP1','KT1','DN1'/ DATA OCHECK/'TE1','MF1','PW1','TE1','DM1','DM1','PR1','HT1'/ CALL RCHECK(INFO,YCHECK,OCHECK) T_pm=70.0 C SET PARAMETER VARIABLES 100 N =PAR(1) D_i =PAR(2) D =PAR(3) delta =PAR(4) L =PAR(5) W =PAR(6) k =PAR(7)/3.6 U_be =PAR(8)/3.6 E_p =PAR(9) Alpha =PAR(10) N_G =PAR(11) RefInd =PAR(12) XKL =PAR(13) c SET INPUT VARIABLES T_in =XIN(1) m_dot =XIN(2)/3600 T_a =XIN(3) I_T =XIN(4)/3.6 wind =XIN(5) I_h =XIN(6)/3.6 I_d =XIN(7)/3.6 rho_g =XIN(8) theta =XIN(9) beta =XIN(10) mu =XIN(11) C_p =XIN(12)*1000 kw =XIN(13)/3.6 rho =XIN(14) Pr =mu*C_p/kw C TO AVOID FLOATING POINT PROBLEMS IF(ABS(m_dot).LE.0.000001) m_dot=0.0 IF (I_T.GT.0.0.AND.Theta.LT.90) GO TO 200 C NO RADIATION TauAlpha=0.0 GO TO 300 200 TauAlpha=TauAlf(beta,N_G,XKL,RefInd,Alpha,I_d,I_h,I_T .,theta,rho_g) 300 pi=4.0*atan(1.0) A_c=N*W*L IF (m_dot.LE.0.)GO TO 400 change=10.0 T_pm=out(4) c HEAT TRANSFER COEFFICIENT L_tube_serpentine=N*(L+W)-W Call HEATTRANSFER(m_dot,L_tube_serpentine,D_i,mu,Pr, kw,h_fi_serp .,Re_serp,Nusselt_serp); c INITIAL GUESS VALUE FOR PLATE TEMPERATURE ITER=0 DO WHILE (change>0.0001.AND.ITER<10000) T_pm_old=T_pm Call LOSSCOEFFICIENT(T_a,T_pm,N_G,beta,wind,E_p,U_be .,U_L); m_flat=sqrt(U_L/(k*delta)) F=tanh(m_flat*(W-D)/2.0)/(m_flat*(W-D)/2.0) F_1=(1.0/U_L)/(W*(1.0/(U_L*(D+(W-D)*F))+1.0/ .(pi*D*h_fi_serp))) F_R =m_dot*C_p/(A_c*U_L)*(1.0-exp(-A_c*U_L*F_1/ . (m_dot*C_p))) Q_useful=A_c*F_R*(I_T*TauAlpha-U_L*(T_in-T_a)) T_pm=T_in+(Q_useful/A_c)/(F_R*U_L)*(1.0-F_R) change=abs(T_pm-T_pm_old) ITER=ITER+1 T_out=Q_useful/(m_dot*C_p)+T_in END DO call PRESSUREDROP(m_dot,D_i,mu,rho,N,L,W,P_loss) GO TO 500 c NO FLOW 400 Q_useful=0.0 change=10.0 ITER=0 U_L=OUT(8) DO WHILE (change>0.0001.AND.ITER<10000) T_pm_old=T_pm Call LOSSCOEFFICIENT(T_a,T_pm,N_G,beta,wind,E_p,U_be .,U_L); T_out=T_a+I_T*TauAlpha/U_L T_pm=T_out change =abs(T_pm-T_pm_old) ITER=ITER+1 END DO P_loss=0.0 c PRINT OUTPUTS FOR DEBUGGING c print*,F_R,T_pm,T_out,eta,Q_useful 500 OUT(1)=T_out OUT(2)=m_dot*3600 OUT(3)=Q_useful*3.6 OUT(4)=T_pm OUT(5)=F_R OUT(6)=TauAlpha OUT(7)=P_loss OUT(8)=U_L*3.6 RETURN 1 END c************************************************************************* Subroutine LOSSCOEFFICIENT(T_a,T_pm,N_G,beta,wind,Epsilon_p, .U_be,U_L) c CALCULATES THE COLLECTOR LOSS COEFFICIENT BASED ON A FUNCTION BY KLEIN IMPLICIT NONE INTEGER N_G REAL h_w, beta,c,e,T_pm,T_a,Epsilon_g,Epsilon_p,f REAL U_t,U_be,U_L,sigma,Wind, T_pm_temp, stf1, stf2 T_pm_temp=T_pm IF (T_pm.LE.T_a) T_pm_temp=T_a+1.0 c APPROXIMATION OF WIND HEAT TRANSFER COEFFICIENT h_w=5.7+3.8*Wind Epsilon_g=0.88 sigma=5.67*10.0**(-8) C USE KLEIN'S TOP LOSS CORRELATION * C=520.0*(1.0-0.000051*beta**2) * e=0.430*(1.0-100.0/(T_pm_temp+273.15)) * f=(1.0+0.089*h_w-0.1166*h_W*Epsilon_p)*(1.0+0.07866*N_G) * U_t=(N_G/(C/(T_pm_temp+273.15)*((T_pm_temp-T_a)/ * .(N_G+f))**e)+1.0/h_w)**(-1)+sigma*((T_pm_temp+273.15) * .+(T_a+273.15))*((T_pm_temp+273.15) * .**2+(T_a+273.15)**2)/((Epsilon_p+0.00591*N_G * .*h_w)**(-1)+(2.0*N_G+f-1.0+0.133*Epsilon_p)/Epsilon_g- * .N_G) F=(1.0-0.04*H_W+5.0E-04*H_W*H_W)*(1.0+0.091*N_G) C=365.9*(1.0-0.00883*beta+0.0001298*beta*beta) STF1=C/(T_pm_temp+273.15)*((T_pm_temp-T_a)/(N_G+F))**0.33 STF1=N_G/STF1+1.0/H_W STF1=1.0/STF1 STF2=1.0/(Epsilon_P+0.05*N_G*(1.0-Epsilon_P))+(2.*N_G+F-1.)/ .Epsilon_G-N_G STF2=SIGMA*((T_PM_TEMP+273.15)**2+(T_A+273.15)**2)* .((T_PM_TEMP+273.15)+(T_A+273.15))/STF2 U_t=(STF1+STF2)+U_BE U_L=(U_t+U_be) RETURN END c************************************************************************ Subroutine HEATTRANSFER(m_dot,L_tube,D_i,mu,Pr,kw,h_fi,Re,Nusselt) IMPLICIT NONE INTEGER Laminar, Turbulent REAL a,b,m_ht,nw,kw,Pr,Re,friction,Nusselt,D_i REAL L_tube,pi,h_fi,mu,m_dot c CONSTANTS a=0.0534 b=0.0335 m_ht=1.15 nw=0.82 pi=4.0*atan(1.0) Re=4.0*m_dot/(pi*D_i*mu) If (Re.LE.0.)THEN Nusselt=3.7 ELSE If (Re.LT.2100.AND.RE.GT.0.) THEN c LAMINAR friction=64.0/Re laminar=1 turbulent=0 Nusselt=3.7+a*(Re*Pr*D_i/L_tube)**m_ht/(1.0+b*(Re*Pr*D_i/L_tube) .**nw) Else c TURBULENT friction=(0.79*log(Re)-1.64)**(-2) laminar=0 turbulent=1 Nusselt= (friction/8.0)*(Re-1000)*Pr/(1.0+12.7*sqrt(friction/8.0)* .(Pr**(2.0/3.0)-1.0)) EndIf EndIf h_fi=Nusselt*kw/D_i RETURN END c***************************************************************************************** FUNCTION TauAlf(beta,N_G,XKL,RefInd,Alpha,I_d,I_h,I_T .,theta,rho_g) IMPLICIT NONE INTEGER N_G REAL beta,XKL,RefInd,Alpha,Rho_d,radconvert,TALN,theta REAL theta_sky,theta_ground,XKATDS,XKATDG,XKATB,XKAT,TALF REAL F_sky, F_gnd, ID_sky, ID_gnd, I_d, I_h,I_T,TauAlf, rho_g radconvert=0.017453 C COVER TRANSMITTANCE AT NORMAL INCIDENCE rho_d=-1.0 TALN=TALF(N_G,0.0,XKL,RefInd,Alpha,Rho_d) C USE THE RELATIONS OF BRANDEMUEHL FOR EFFECTIVE INCIDENCE ANGLES C FOR DIFFUSE RADIATION theta_sky=59.68-0.1388*beta+0.001497*beta*beta theta_ground=90.0-0.5788*beta+0.002693*beta*beta C DIFFUSE SKY RADIATION TAUALPHA RATIO XKATDS=TALF(N_G,theta_sky,XKL,RefInd,Alpha,Rho_d)/TALN C GROUND REFLECTED RADIATION TAUALPHA RATIO XKATDG=TALF(N_G,theta_ground,XKL,RefInd,Alpha,Rho_d)/TALN C BEAM RADIATION TAUALPHA RATIO XKATB=TALF(N_G,theta,XKL,RefInd,Alpha,Rho_d)/TALN C VIEW FACTORS F_SKY=(1.0+cos(beta*radconvert))/2.0 F_GND=(1.0-cos(beta*radconvert))/2.0 C SKY DIFFUSE RADIATION ID_SKY=F_SKY*I_D C GROUND DIFFUSE RADIATION ID_GND=Rho_g*F_GND*I_H C OVERALL TAUALPHA RATIO XKAT=(XKATB*(I_T-ID_SKY-ID_GND)+XKATDS*ID_SKY+XKATDG*ID_GND)/I_T TAUALf=TALN*XKAT RETURN END C***************************************************************** Subroutine PRESSUREDROP(m_dot,D_i,mu,rho,N,L,W,P_loss) IMPLICIT NONE INTEGER N REAL m_dot, D_i, mu, rho, L, pi, f,Re, v_serp REAL PipeLossCoefficient, L_eq, L_total, P_loss, W C Pressure Drop for the serpentine collector C Find the moody friction factor pi=4.0*atan(1.0) Re=4.0*m_dot/(pi*D_i*mu) If (Re.LT.0.001) Then !NO FLOW f=0.0 Else If (Re.LT.2100) Then !LAMINAR f=64.0/Re Else !TURBULENT f=(0.79*log(Re)-1.64)**(-2) End If End If v_serp=m_dot/(rho*pi*(D_i/2.0)**2) PipeLossCoefficient=0.5 L_eq=PipeLossCoefficient*D_i/f L_total=N*(L+W)-W+(2*N-2)*L_eq P_loss=f*L_total/D_i*rho*v_serp**2/2.0*10.0**(-3) RETURN END