SUBROUTINE TYPE52(TIME,XIN,OUT,T,DTDT,PAR,INFO,ICNTRL,*) C----------------------------------------------------------------------------------------------------------------------- C THIS ROUTINE DETERMINES THE THERMAL PERFORMANCE OF A COOLING COIL * C USING AN EFFECTIVENESS MODEL AS DESCRIBED IN "METHODOLOGIES FOR * C THE DESIGN AND CONTROL OF CHILLED WATER SYSTEMS", PHD THESIS, * C SOLAR ENERGY LABORATORY, UNIVERSITY OF WISCONSIN, J.E. BRAUN, 1988 * C**************************** PARAMETERS *************************** C C MODE = SIMPLE OR DETAILED CALCULATIONS C IU = UNITS: 1-SI, 2-ENGLISH C NROWS = NUMBER OF ROWS OF TUBES C NTUBES = NUMBER TUBES PER ROW C LT = LENGTH OF TUBES IN EACH ROW C WD = WIDTH OF DUCT PERPENDICULAR TO THE TUBES C DO = TUBE OUTER DIAMETER C DI = TUBE INNER DIAMETER C KT = TUBE THERMAL CONDUCTIVITY C FT = FIN THICKNESS C FS = FIN SPACING C NFINS = NUMBER OF FINS ON ONE TUBE C KF = FIN THERMAL CONDUCTIVITY C FMODE = CONTINUOUS FLAT PLATE OR ANNULAR FINS C ADIST = DISTANCE BETWEEN TUBE CENTERS ACROSS ROW C DE = DIAMETER OF ANNULAR FIN C CDIST = DISTANCE BETWEEN TUBE ROW CENTERLINES (FLOW DIRECTION) C----------------------------------------------------------------------------------------------------------------------- ! Copyright © 2004 Solar Energy Laboratory, University of Wisconsin-Madison. All rights reserved. C----------------------------------------------------------------------------------------------------------------------- C USE STATEMENTS USE TrnsysFunctions C----------------------------------------------------------------------------------------------------------------------- C----------------------------------------------------------------------------------------------------------------------- C REQUIRED BY THE MULTI-DLL VERSION OF TRNSYS !DEC$ATTRIBUTES DLLEXPORT :: TYPE52 !SET THE CORRECT TYPE NUMBER HERE C----------------------------------------------------------------------------------------------------------------------- C----------------------------------------------------------------------------------------------------------------------- C TRNSYS DECLARATIONS IMPLICIT NONE DOUBLE PRECISION XIN,OUT,TIME,PAR,T,DTDT,TIME0,TFINAL,DELT INTEGER*4 INFO(15),NP,NI,NOUT,ND,NPAR,IUNIT,ITYPE,ICNTRL, 1 NSTORED CHARACTER*3 YCHECK,OCHECK C----------------------------------------------------------------------------------------------------------------------- C----------------------------------------------------------------------------------------------------------------------- C USER DECLARATIONS - SET THE MAXIMUM NUMBER OF PARAMETERS (NP), INPUTS (NI), C OUTPUTS (NOUT), AND DERIVATIVES (ND) THAT MAY BE SUPPLIED FOR THIS TYPE PARAMETER (NP=16,NI=5,NOUT=12,ND=0,NSTORED=0) C----------------------------------------------------------------------------------------------------------------------- C----------------------------------------------------------------------------------------------------------------------- C REQUIRED TRNSYS DIMENSIONS DIMENSION XIN(NI),OUT(NOUT),PAR(NP),YCHECK(NI),OCHECK(NOUT), 1 T(ND),DTDT(ND) C----------------------------------------------------------------------------------------------------------------------- C----------------------------------------------------------------------------------------------------------------------- C DECLARATIONS AND DEFINITIONS FOR THE USER-VARIABLES DOUBLE PRECISION MA,MW,MHUA,KW,MHUW,LT,KT,KF,JF,MCONV,NTUD,NTUW, & NTUO,NTUI,MSTAR,LCOIL,LIMITW,LIMITD,RNUM,TCONV,A1,A0,HCONV,H,B0, & B1,CTOF,C1,C0,TS,HS,TSAT,PRW,CPA,CPW,CPV,HFG,PRA,PATM,TOL,PI,HSAT & ,WD,DO,DI,FT,FS,ADIST,DE,CDIST,PSYDAT,RI,RO,RE,AC,FH,AX,AI,ATO,AF & ,AO,FRATIO,WA1,TW1,UM,TA1,TA2,WA2,TW2,QCOIL,QSENS,QLAT,RHA2,RHA1, & QDAIR,FDRY,GW,TPROP,REI,HI,HIT,HIL,DH,CJ1,CJ2,GA,REO,HDO,ALPHA, & BETA,EFF,FIN,EFFD,CPM,RA,TDP,HA1,HW1,CSTAR,TEMP,C,EPSD,TA2D,TS2D, & TW2W,HW2,CS,HWO,EFFW,ERAT,EPSW,G,G1,TOLD,G2,DGDT,TS1W,CK,TW2D, & TMAX,TMIN,GMAX,TWX,GMIN,SGNMIN,DMAX1,DMIN1,HA2,HSAVG,TSAVG,TAX, & HAX2,HAX,HSW,TSW, & GZ,NU_LAM,RE_TUR,FD,NU_TUR !Added by Avatar Lee on 17/07/2013 INTEGER MODE,FMODE,EMODE,ROUND,IU,IMAX,NROWS,NTUBES,NFINS,IWARN, & ISTAT,ITER LOGICAL DRY,WET,MAXDRY DIMENSION CPA(2),MHUA(2),CPW(2),KW(2),CPV(2),PSYDAT(9),HFG(2), &A0(2),A1(2),B0(2),B1(2),C0(2),C1(2),MCONV(2) DATA CPA/1.005, 0.240/, MHUA/0.0659, 0.0443/, CPW/4.186, 1.000/, . KW/2.043, 0.328/, CPV/1.884, 0.450/, . HFG/2452., 1054./, PRA/0.72/, PATM/1.0/ DATA TOL/1.E-04/, IMAX/20/, PI/3.14159/, IUNIT/0/ DATA A0/0., 32./, A1/1., 1.8/, B0/0., 7.687/, B1/1.,0.43/ . C0/32., 0./, C1/1.8, 1./, MCONV/1.488, 1.0/ C----------------------------------------------------------------------------------------------------------------------- C EMODE DETERMINES HOW ERRORS ARE HANDLED FROM CALLS TO THE PSYCH C SUBROUTINE. IF EMODE IS: 0 - NO ERROR MESSAGES WILL BE PRINTED, C 1 - ERROR MESSAGES WILL BE PRINTED ONLY ONCE PER SIMULATION, C 2 - ERROR MESSAGES WILL BE PRINTED EVERY TIMESTEMP THAT THEY OCCUR. C DATA EMODE/1/ C C*********************** STATEMENT FUNCTIONS ************************** C C-- ROUND-OFF REAL NUMBERS TO NEAREST INTEGER C ROUND(RNUM) = JFIX(RNUM + SIGN(0.5,RNUM)) C C-- UNIT CONVERSIONS FOR TEMPERATURE AND ENTHALPY C C TEMPERATURE: SI TO ENGLISH (IU=2), SI TO SI (IU=1) C ENTHALPY: ENGLISH TO SI (IU=2), SI TO SI (IU=1) C TEMPERATURE: ENGLISH TO ENGLISH (IU=2), SI TO ENGLISH (IU = 1) C TCONV(T) = T*A1(IU) + A0(IU) HCONV(H) = (H - B0(IU))/B1(IU) CTOF(T) = T*C1(IU) + C0(IU) C C-- CORRELATION FOR SATURATION TEMPERATURE IN TERMS OF SATURATION C ENTHALPY: CORRELATION IS IN SI UNITS AND IS GOOD FROM 9.473 KJ/KG C TO 355.137 KJ/KG (0 C TO 55 C). C (CORRELATION IS FOR A TOTAL SYSTEM PRESSURE OF 1 ATMOSPHERE) C TS(HS) = -5.79013 + 6.64030E-01*HS - 5.07802E-03*HS**2 + . 2.80381E-05*HS**3 - 9.47051E-08*HS**4 + . 1.72758E-10*HS**5 - 1.29547E-13*HS**6 C TSAT(HSAT) = TCONV(TS(HCONV(HSAT))) C C-- CORRELATION FOR KINEMATIC WATER VISCOSITY IN TERMS OF TEMPERATURE: C CORRELATION IS IN ENGLISH UNITS (LBM/FT-HR AND DEG. F) AND C IS GOOD FROM 32 F TO 200 F. C MHUW(T) = MCONV(IU)*(7.67201 - 1.39290E-01*T + 1.23322E-03*T**2 . - 5.32425E-06*T**3 + 8.87455E-09*T**4) C C-- CORRELATION FOR PRANTL NUMBER OF WATER IN TERMS OF TEMPERATURE: C TEMPERATURE IS IN ENGLISH UNITS (DEG. F) AND THE CORRELATION IS C GOOD FROM 32 F TO 200 F. C PRW(T) = 2.73327E01 - 6.30480E-01*T + 7.56704E-03*T**2 - . 5.04074E-05*T**3 + 1.74546E-07*T**4 - 2.43923E-10*T**5 C----------------------------------------------------------------------------------------------------------------------- C TRNSYS FUNCTIONS TIME0=getSimulationStartTime() TFINAL=getSimulationStopTime() DELT=getSimulationTimeStep() C SET THE VERSION INFORMATION FOR TRNSYS IF(INFO(7).EQ.-2) THEN INFO(12)=16 RETURN 1 ENDIF C----------------------------------------------------------------------------------------------------------------------- C----------------------------------------------------------------------------------------------------------------------- C DO ALL THE VERY LAST CALL OF THE SIMULATION MANIPULATIONS HERE IF (INFO(8).EQ.-1) THEN RETURN 1 ENDIF C----------------------------------------------------------------------------------------------------------------------- C----------------------------------------------------------------------------------------------------------------------- C PERFORM ANY "AFTER-ITERATION" MANIPULATIONS THAT ARE REQUIRED IF(INFO(13).GT.0) THEN RETURN 1 ENDIF C----------------------------------------------------------------------------------------------------------------------- C----------------------------------------------------------------------------------------------------------------------- C DO ALL THE VERY FIRST CALL OF THE SIMULATION MANIPULATIONS HERE IF (INFO(7).EQ.-1) THEN C RETRIEVE THE UNIT NUMBER AND TYPE NUMBER FOR THIS COMPONENT FROM THE INFO ARRAY IUNIT=INFO(1) ITYPE=INFO(2) C SET SOME INFO ARRAY VARIABLES TO TELL THE TRNSYS ENGINE HOW THIS TYPE IS TO WORK INFO(6)=12 INFO(9)=1 INFO(10)=0 C SET THE REQUIRED NUMBER OF INPUTS, PARAMETERS AND DERIVATIVES THAT THE USER SHOULD SUPPLY IN THE INPUT FILE NPAR=15 C CALL THE TYPE CHECK SUBROUTINE TO COMPARE WHAT THIS COMPONENT REQUIRES TO WHAT IS SUPPLIED CALL TYPECK(1,INFO,NI,NPAR,ND) C SET THE VARIABLE TYPES DATA YCHECK/'TE1','DM1','MF1','TE1','MF1'/ DATA OCHECK/'TE1','DM1','MF1','TE1','MF1','PW1','PW1','PW1', 1 'DM1','PC1','NAV','NAV'/ C CALL THE INPUT-OUTPUT CHECK SUBROUTINE TO SET THE CORRECT INPUT AND OUTPUT UNITS CALL RCHECK(INFO,YCHECK,OCHECK) C RETURN TO THE CALLING PROGRAM RETURN 1 ENDIF C----------------------------------------------------------------------------------------------------------------------- C----------------------------------------------------------------------------------------------------------------------- C DO ALL OF THE INITIAL TIMESTEP MANIPULATIONS HERE - THERE ARE NO ITERATIONS AT THE INTIAL TIME IF (TIME.LT.(TIME0+DELT/2.D0)) THEN C SET THE UNIT NUMBER FOR FUTURE CALLS IUNIT=INFO(1) ITYPE=INFO(2) C READ IN THE VALUES OF THE PARAMETERS IN SEQUENTIAL ORDER MODE=ROUND(PAR(1)) IF(MODE.NE.1 .AND. MODE.NE.2) CALL TYPECK(-4,INFO,0,0,0) IU=1 NROWS = ROUND(PAR(2)) NTUBES = ROUND(PAR(3)) LT = PAR(4) WD = PAR(5) DO = PAR(6) DI = PAR(7) KT = PAR(8) FT = PAR(9) FS = PAR(10) NFINS = ROUND(PAR(11)) KF = PAR(12) FMODE=ROUND(PAR(13)) IF(FMODE.NE.1 .AND. FMODE.NE.2) CALL TYPECK(-4,INFO,0,0,0) ADIST = PAR(14) IF (FMODE .EQ. 2) DE = ADIST CDIST = PAR(15) PSYDAT(1) = PATM C** COIL AREAS ** C C RI & RO - INSIDE AND OUTSIDE TUBE RADIUS C RE - ANNULAR FIN RADIUS OF EQUIVALENT RADIUS FOR C CONTINUOUS FLAT FINS C FH - FIN HEIGHT C AI & AO - INSIDE AND OUTSIDE COIL SURFACE AREAS C FRATIO - RATIO OF FIN AREA TO OUTSIDE COIL SURFACE AREA C LCOIL - LENGTH OF THE COIL SECTION IN THE AIR FLOW DIRECTION C (SEE KAYS & LONDON, "COMPACT HEAT EXCHANGERS") C AC - FLOW CROSS SECTIONAL AREA C AX - CROSS SECTIONAL AREA OF ONE TUBE C RI = DI/2. RO = DO/2. IF (FMODE .EQ. 1) THEN RE = DSQRT(ADIST*CDIST/PI) AC = LT*WD-NFINS*FT*(NTUBES*ADIST+ADIST/2.) . -(NTUBES*LT*DO-(NFINS*FT*DO)) ELSE RE = DE/2. AC = LT*WD-NTUBES*LT*DO-NTUBES*NFINS*(DE-DO)*FT END IF LCOIL = NROWS * CDIST FH = RE - RO AX = PI*RI**2 AI = NTUBES*NROWS*PI*DI*LT ATO = NTUBES*NROWS*PI*DO*(LT-NFINS*FT) AF = NFINS*NTUBES*NROWS*2.*PI*(RE**2-RO**2) AO = ATO + AF FRATIO = AF/AO C C** TUBE HEAT TRANSFER COEFFICIENT ** C BASED UPON INSIDE AREA C UM = KT/RI/DLOG(RO/RI) C PERFORM ANY REQUIRED CALCULATIONS TO SET THE INITIAL VALUES OF THE OUTPUTS HERE OUT(1) = XIN(1) OUT(2) = XIN(2) OUT(3) = XIN(3) OUT(4) = XIN(4) OUT(5) = XIN(5) OUT(6) = 0. OUT(7) = 0. OUT(8) = 0. OUT(9) = 1.0 OUT(10) = 0.0 C RETURN TO THE CALLING PROGRAM RETURN 1 ENDIF C----------------------------------------------------------------------------------------------------------------------- C----------------------------------------------------------------------------------------------------------------------- C *** ITS AN ITERATIVE CALL TO THIS COMPONENT *** C----------------------------------------------------------------------------------------------------------------------- C----------------------------------------------------------------------------------------------------------------------- C RE-READ THE PARAMETERS IF ANOTHER UNIT OF THIS TYPE HAS BEEN CALLED SINCE THE LAST TIME THE PARAMETERS C WERE READ IN IF(INFO(1).NE.IUNIT) THEN C RESET THE UNIT NUMBER IUNIT=INFO(1) ITYPE=INFO(2) C READ IN THE VALUES OF THE PARAMETERS IN SEQUENTIAL ORDER MODE=ROUND(PAR(1)) IF(MODE.NE.1 .AND. MODE.NE.2) CALL TYPECK(-4,INFO,0,0,0) IU=1 NROWS = ROUND(PAR(2)) NTUBES = ROUND(PAR(3)) LT = PAR(4) WD = PAR(5) DO = PAR(6) DI = PAR(7) KT = PAR(8) FT = PAR(9) FS = PAR(10) NFINS = ROUND(PAR(11)) KF = PAR(12) FMODE=ROUND(PAR(13)) IF(FMODE.NE.1 .AND. FMODE.NE.2) CALL TYPECK(-4,INFO,0,0,0) ADIST = PAR(14) IF (FMODE .EQ. 2) DE = ADIST CDIST = PAR(15) PSYDAT(1) = PATM C** COIL AREAS ** C C RI & RO - INSIDE AND OUTSIDE TUBE RADIUS C RE - ANNULAR FIN RADIUS OF EQUIVALENT RADIUS FOR C CONTINUOUS FLAT FINS C FH - FIN HEIGHT C AI & AO - INSIDE AND OUTSIDE COIL SURFACE AREAS C FRATIO - RATIO OF FIN AREA TO OUTSIDE COIL SURFACE AREA C LCOIL - LENGTH OF THE COIL SECTION IN THE AIR FLOW DIRECTION C (SEE KAYS & LONDON, "COMPACT HEAT EXCHANGERS") C AC - FLOW CROSS SECTIONAL AREA C AX - CROSS SECTIONAL AREA OF ONE TUBE C RI = DI/2. RO = DO/2. IF (FMODE .EQ. 1) THEN RE = DSQRT(ADIST*CDIST/PI) AC = LT*WD-NFINS*FT*(NTUBES*ADIST+ADIST/2.) . -(NTUBES*LT*DO-(NFINS*FT*DO)) ELSE RE = DE/2. AC = LT*WD-NTUBES*LT*DO-NTUBES*NFINS*(DE-DO)*FT END IF LCOIL = NROWS * CDIST FH = RE - RO AX = PI*RI**2 AI = NTUBES*NROWS*PI*DI*LT ATO = NTUBES*NROWS*PI*DO*(LT-NFINS*FT) AF = NFINS*NTUBES*NROWS*2.*PI*(RE**2-RO**2) AO = ATO + AF FRATIO = AF/AO C C** TUBE HEAT TRANSFER COEFFICIENT ** C BASED UPON INSIDE AREA C UM = KT/RI/DLOG(RO/RI) ENDIF C----------------------------------------------------------------------------------------------------------------------- C****************************** INPUTS **************************** C C TA1 = ENTERING AIR DRY BULB TEMPERATURE (C, F) C WA1 = ENTERING AIR HUMIDITY RATIO C MA = MASS FLOW RATE OF AIR (KG/HR, LBM/HR) C TW1 = ENTERING WATER TEMPERATURE (C, F) C MW = MASS FLOW RATE OF WATER (KG/HR LBM/HR) C TA1 = XIN(1) WA1 = XIN(2) MA = XIN(3) TW1 = XIN(4) MW = XIN(5) C CHECK THE INLET AIR STATE PSYDAT(1) = PATM PSYDAT(2) = TA1 PSYDAT(6) = WA1 CALL PSYCHROMETRICS(TIME,INFO,IU,4,0,PSYDAT,EMODE,ISTAT,*99) CALL LINKCK('TYPE52','PSYCHROMETRICS ',1,99) 99 TA1 = PSYDAT(2) WA1 = PSYDAT(6) RHA1 = PSYDAT(4) C----------------------------------------------------------------------------------------------------------------------- C************************* COIL ANALYSIS ************************* C C TA2 = DRY-BULB TEMPERATURE OF AIR LEAVING COIL C WA2 = ABSOLUTE HUMIDITY OF AIR LEAVING COIL C RHA2 = RELATIVE HUMIDITY OF AIR LEAVING COIL C TW2 = DRY-BULB TEMPERATURE OF WATER LEAVING COIL C QCOIL = TOTAL HEAT TRANSFER TO COIL C QSENS = SENSIBLE HEAT REMOVED FROM MOIST AIR C QLAT = LATENT HEAT REMOVED FROM MOIST AIR C QDAIR = SENSIBLE HEAT REMOVED FROM AIR ONLY (DRY AIR) C FDRY = FRACTION OF COIL WHICH IS DRY C GW = WATER MASS FLOW RATE DIVIDED BY CROSS-SECTIONAL AREA C REI = INSIDE REYNOLDS NUMBER C HI = INSIDE COEFFICIENT OF CONVECTION C GA = AIR MASS FLOW RATE DIVIDED BY CROSS-SECTIONAL AREA C REO = OUTSIDE REYNOLDS NUMBER C CPM = MOIST AIR HEAT CAPACITY C NTUO = NTU'S OUTSIDE C NTUI = NTU'S INSIDE C TDP = AIR DEW-POINT C HA1 = ENTHALPY OF ENTERING AIR C HA2 = ENTHALPY OF EXITING AIR C *** NO FLOW *** IF (MA .LE. 1.0E-06 .OR. MW .LE. 1.0E-06) THEN TA2 = TA1 IF (MA .LE. 1.0E-02 .AND. MW .GT. 1.0E-02) TA2 = TW1 WA2 = WA1 RHA2 = RHA1 TW2 = TW1 QCOIL = 0.0 QSENS = 0.0 QLAT = 0.0 QDAIR = 0.0 FDRY = 1.0 GOTO 500 END IF C *** FLOW *** C ** INSIDE FLUID COEFFICIENT ** C FIND THE REYNOLDS NUMBER GW = MW/NTUBES/AX TPROP = CTOF(TW1) IF ((32.0.GT.TPROP .OR. TPROP.GT.200.) .AND. . (TIME-OUT(12)) .GT. DELT/2.) THEN CALL MESSAGES(-1,'THE CORRELATIONS FOR THE VISCOSITY AND PRANDT &L # OF WATER USED IN THE COOLING COIL MODEL WERE USED WITH A VALUE & OUTSIDE OF THEIR INTENDED RANGE','WARNING',INFO(1),INFO(2)) IWARN=IWARN+1 OUT(12) = TIME END IF REI = GW*DI/MHUW(TPROP) C Added by Avatar Lee on 17/07/2013 C Calculation of laminar Nusselt number based on Nellis G, Klein S. Heat Transfer. New York: Cambridge University Press, 2009. GZ = DI*MIN(REI,2300.0D0)*PRW(TPROP)/LT NU_LAM = 3.66+(0.049+0.02/PRW(TPROP))*GZ**1.12/(1.0+0.065*GZ**0.7) C Calculation of turbulent Nusselt number based on Gnielinski V. New equations for heat and mass transfer in turbulent pipe and channel flow. C International Chemical Engineering 1976;16:359-68. RE_TUR = MAX(REI,3000.0D0) FD = 4.0*(-0.0015702/LOG(RE_TUR)+0.3942031/(LOG(RE_TUR))**2.0 &+2.5341533/(Log(RE_TUR))**3.0) NU_TUR = (FD/8.0)*(RE_TUR-1000.0)*PRW(TPROP)*(1.0+(DI/LT)**0.7) &/(1.0+12.7*SQRT(FD/8.0)*(PRW(TPROP)**(2.0/3.0)-1.0)) C Revision ended... C IF THE REYNOLDS NUMBER IS GREATER THEN 3000, USE THE CORRELATION C FOR TURBULENT WATER FLOW IN A TUBE. IF THE REYNOLDS NUMBER IS C LESS THAN 2000, USE THE CORRELATION FOR LAMINAR WATER FLOW IN A C TUBE. FOR A REYNOLDS NUMBER BETWEEN 2000 AND 3000, A LINEAR C RELATIONSHIP BETWEEN THE CONVECTION COEFFICIENTS FOR THE TURBULENT C AND LAMINAR CASES IS USED. THE INLET TEMPERATURE IS USED FOR C PROPERTY EVALUATIONS. IF (REI .GE. 3000.) THEN C HI = 0.023*KW(IU)/DI*REI**0.8*PRW(TPROP)**0.4 HI = NU_TUR*KW(IU)/DI !Added by Avatar Lee on 17/07/2013 ELSE IF (REI .LE. 2300.) THEN C HI = 4.36*KW(IU)/DI HI = NU_LAM*KW(IU)/DI !Added by Avatar Lee on 17/07/2013 ELSE C HIT = 0.023*KW(IU)/DI*REI**0.8*PRW(TPROP)**0.4 HIT = NU_TUR*KW(IU)/DI !Added by Avatar Lee on 17/07/2013 C HIL = 4.36*KW(IU)/DI HIL = NU_LAM*KW(IU)/DI !Added by Avatar Lee on 17/07/2013 HI = HIL + (REI - 2000.0) * (HIT - HIL) / 1000.0 END IF C ** OUTSIDE COEFFICIENT FOR DRY SURFACES ** C CORRELATION FOR OUTSIDE HEAT TRANSFER COEFFICIENT IS FROM C "FINNED TUBE HEAT EXCHANGER: CORRELATION OF DRY SURFACE HEAT C TRANSFER DATA," A.H. ELMAHDY, ASHRAE TRANSACTIONS. THE C HYDRAULIC DIAMETER IS: 2 * THE FLOW LENGTH OF THE HEAT EXCHANGER * C THE FLOW CROSS SECTIONAL AREA / TOTAL HEAT TRANSFER AREA. C THE FIN EFFICIENCY EQUATION IS FOR ANNULAR FINS. STRAIGHT FINS C ARE TREATED AS ANNULAR FINS BY FINDING AN EQUIVALENT ANNULAR FIN C WITH THE SAME SURFACE AREA. DH = 4.*LCOIL*AC/AO CJ1 = 0.159*(FT/FH)**0.141*(DH/FT)**0.065 CJ2 = -0.323*(FT/FH)**0.049*(FS/FT)**0.077 GA = MA/AC REO = GA*DH/MHUA(IU) JF = CJ1*REO**CJ2 HDO = GA*CPA(IU)*JF/PRA**0.6667 C FIN EFFICIENCY CALCULATIONS ALPHA = RO/RE BETA = RE*DSQRT(2.*HDO/KF/FT) EFF = FIN(ALPHA,BETA) EFFD = 1. - FRATIO*(1.-EFF) C ** DETERMINE NTU'S ** C THE AIR SPECIFIC HEAT IS FOR MOIST AIR CPM = CPA(IU) + WA1*CPV(IU) NTUO = EFFD*HDO*AO/(MA*CPM) NTUI = AI/(1./HI + 1./UM)/(MW*CPW(IU)) RA = MW/MA C ** ENTHALPIES FOR ENTERING CONDITIONS ** PSYDAT(2) = TA1 PSYDAT(6) = WA1 CALL PSYCHROMETRICS(TIME,INFO,IU,4,0,PSYDAT,EMODE,ISTAT,*101) CALL LINKCK('TYPE52','PSYCHROMETRICS ',1,99) 101 TDP = PSYDAT(5) HA1 = PSYDAT(7) PSYDAT(2) = TW1 PSYDAT(3) = TW1 CALL PSYCHROMETRICS(TIME,INFO,IU,1,0,PSYDAT,EMODE,ISTAT,*102) CALL LINKCK('TYPE52','PSYCHROMETRICS ',1,99) 102 HW1 = PSYDAT(7) C ** DRY ANALYSIS ** C DETERMINES THE AIR EFFECTIVENESS (EPSD) AND OUTLET WATER TEMPERATURE C (TW2D) ASSUMING THE COIL IS DRY C IF THE SURFACE TEMPERATURE AT THE OUTLET IS GREATER THAN THE AIR C ENTERING DEWPOINT THEN COIL IS COMPLETELY DRY, OTHERWISE DO THE WET C ANALYSIS CSTAR = CPM/RA/CPW(IU) NTUD = NTUO/(1. + NTUO/NTUI*CSTAR) TEMP = -NTUD*(1.-CSTAR) IF (TEMP .LT. 50.0) THEN C = DEXP(TEMP) EPSD = (1. - C)/(1. - CSTAR*C) ELSE EPSD = 1/CSTAR ENDIF TW2D = TW1 + EPSD*CSTAR*(TA1 - TW1) TA2D = TA1 - EPSD*(TA1 - TW1) TS2D = TW1 + CSTAR*NTUD/NTUI*(TA2D - TW1) C DRY = TS2D .GT. TDP FDRY = 1.0 IF(.NOT. DRY) THEN FDRY = 0.0 C ** WET COIL ** C DETERMINES THE AIR EFFECTIVENESS (EPSW) AND OUTLET WATER TEMPERATURE C (TW2W) ASSUMING THE COIL IS WET. C THE AVERAGE SATURATION SPECIFIC HEAT (CS) DEPENDS UPON THE OUTLET C CONDITIONS; THEREFORE AN ITERATIVE SOLUTION. C THE WET SURFACE HEAT TRANSFER COEFFICIENT (HWO) AND THEREFORE THE FIN C EFFICIENCY (EFFW) DEPEND UPON THE SATURATION SPECIFIC HEAT. C IF THE SURFACE TEMPERATURE AT WATER OUTLET IS LESS THAN THE DEWPOINT C OF THE ENTERING AIR THEN THE COIL IS COMPLETELY WET. TW2W = TW2D ITER = 0 10 ITER = ITER + 1 PSYDAT(2) = TW2W PSYDAT(3) = TW2W CALL PSYCHROMETRICS(TIME,INFO,IU,1,0,PSYDAT,EMODE,ISTAT,*104) CALL LINKCK('TYPE52','PSYCHROMETRICS ',1,99) 104 HW2 = PSYDAT(7) CS = (HW2 - HW1)/(TW2W - TW1) MSTAR = CS/RA/CPW(IU) HWO = CS*HDO/CPM ALPHA = RO/RE BETA = RE*DSQRT(2.*HWO/KF/FT) EFF = FIN(ALPHA,BETA) EFFW = 1. - FRATIO*(1.-EFF) ERAT = EFFW/EFFD NTUW = ERAT*NTUO/(1. + ERAT*NTUO/NTUI*MSTAR) TEMP = -NTUW*(1.-MSTAR) IF (TEMP .LT. 50.0) THEN C = DEXP(TEMP) EPSW = (1. - C)/(1. - MSTAR*C) ELSE EPSW = 1/MSTAR ENDIF C USE A SECANT METHOD TO CONVERGE ON TW2 G=TW2W-(TW1 + EPSW*(HA1 - HW1)/RA/CPW(IU)) IF (DABS(G) .GT. TOL .AND. ITER .LT. IMAX) THEN IF (ITER .EQ. 1) THEN G1 = G TOLD = TW2W TW2W = TOLD - G1 ELSE G2 = G DGDT = (G2 - G1)/(TW2W - TOLD) DGDT = SIGN((DMAX1(1.0E-06,DABS(DGDT))),DGDT) TOLD = TW2W TW2W = TOLD - G2/DGDT IF (DABS(TW2W-TOLD) .LT. TOL) TW2W = TOLD - G2 G1 = G2 END IF GOTO 10 END IF TS1W = TW2W + CSTAR/CPM*NTUW/NTUI*(HA1 - HW2) ENDIF C CHECK TO SEE IF COIL COMPLETELY WET WET = ((.NOT. DRY) .AND. TS1W .LT. TDP) IF(.NOT.(DRY) .AND. .NOT.(WET)) THEN IF (MODE .EQ. 1) THEN DRY = TW2D .GT. TW2W FDRY = -1.0 ELSE IF (MODE .EQ. 2) THEN C ** COMBINED WET AND DRY COIL PERFORMANCE ** C C AN ITERATIVE BY-SECTION METHOD IS USEDD TO FIND TW2. A FUNCTION C G, WHICH EQUALS ZERO AT THE SOLUTION, IS DEFINED FROM AN ENERGY C BALANCE. THE VALUE OF G IS FOUND FOR TW2 CORRESPONDING TO FDRY=0 C AND FDRY=1. IF BOTH VALUES OF G ARE POSITIVE OR BOTH ARE NEGATIVE, C THE LOGICAL VARIABLES "WET" AND "DRY" ARE APPROPRIATELY SET. (IE. C THE COIL WILL BE CONSIDERED WET IF THE WET LIMIT GIVES THE G VALUE C CLOSEST TO 0, OR THE COIL WILL BE CONSIDERED DRY IF THE DRY LIMIT C GIVES THE G VALUE CLOSEST TO 0.) CK = NTUD*(1. - CSTAR) LIMITW=(TA1*(CSTAR-1.+CK/NTUO)+TDP*(1-CSTAR))/(CK/NTUO) LIMITD=(TA1*(CSTAR-(1.-CK/NTUO)*DEXP(-CK))+TDP*(1.-CSTAR))/ & (1.-(1.-CK/NTUO)*DEXP(-CK)) IF(LIMITD .GT. LIMITW) THEN TMAX = 0.999999 * LIMITD TMIN = 1.000001 * LIMITW MAXDRY = .TRUE. ELSE TMAX = 0.999999 * LIMITW TMIN = 1.000001 * LIMITD MAXDRY = .FALSE. ENDIF CALL CALCG(TMAX,TDP,TA1,TW1,HA1,HW1,RO,RA,RE,KF,FT,NTUO, & NTUI,NTUD,HDO,CPM,CSTAR,CPW(IU),EFFD,FRATIO,INFO, & IU,GMAX,FDRY,EPSD,TWX,ERAT,EMODE) CALL CALCG(TMIN,TDP,TA1,TW1,HA1,HW1,RO,RA,RE,KF,FT,NTUO, & NTUI,NTUD,HDO,CPM,CSTAR,CPW(IU),EFFD,FRATIO,INFO, & IU,GMIN,FDRY,EPSD,TWX,ERAT,EMODE) IF ((GMAX .GT. 0.0 .AND. GMIN .GT. 0.0) .OR. . (GMAX .LT. 0.0 .AND. GMIN .LT. 0.0)) THEN IF (DABS(GMAX) .LT. DABS(GMIN) ) THEN DRY = MAXDRY WET = .NOT. MAXDRY ELSE WET = MAXDRY DRY = .NOT. MAXDRY ENDIF ELSE TW2 = TMIN ITER = -1 C ITERATION LOOP TO DETERMINE TW2. 15 ITER = ITER + 1 CALL CALCG(TW2,TDP,TA1,TW1,HA1,HW1,RO,RA,RE,KF,FT,NTUO, & NTUI,NTUD,HDO,CPM,CSTAR,CPW(IU),EFFD,FRATIO, & INFO,IU,G,FDRY,EPSD,TWX,ERAT,EMODE) IF ((DABS(G) .GT. TOL .AND. ITER .LT. IMAX) .AND. . (DABS(TMIN - TMAX) .GT. TOL)) THEN IF (ITER .EQ. 0) THEN SGNMIN = SIGN(1.,G) C USE THE MAXIMUM TW2 FROM THE WET AND DRY ANALYSIS AS AN INITIAL GUESS TW2 = DMAX1(TMIN,DMIN1(DMAX1(TW2W,TW2D),TMAX)) GOTO 15 ELSE IF (ITER .EQ. 1) THEN G1 = G TOLD = TW2 TW2 = TOLD - G1 ELSE G2 = G DGDT = (G2 - G1)/(TW2 - TOLD) DGDT = SIGN((DMAX1(1.0E-06,DABS(DGDT))),DGDT) TOLD = TW2 TW2 = TOLD - G2/DGDT G1 = G2 END IF C NARROW THE BOUNDS TO THE SOLUTION IF (G1 .LT. 0.) THEN IF (SGNMIN .LT. 0.) THEN TMIN = TOLD ELSE TMAX = TOLD END IF ELSE IF (SGNMIN .GT. 0.) THEN TMIN = TOLD ELSE TMAX = TOLD END IF END IF C USE A BI-SECTION IF THE NEW VALUE OF TW2 IS OUTSIDE THE BOUNDS IF (TW2.GE.TMAX.OR.TW2.LE.TMIN) TW2 = (TMAX + TMIN)/2. GOTO 15 END IF END IF C * END OF COMBINED COIL ANALYSIS END IF END IF C ** END OF COIL ANALYSIS C C *** NET COIL PERFORMANCE *** C FOR MODE 1, UTILIZE THE RESULTS OF THE ANALYSIS THAT GIVES THE C GREATEST HEAT TRANSFER (I.E. HIGHEST WATER RETURN TEMPERATURE) C IF COIL IS NOT COMPLETELY WET OR DRY C ** DRY ** IF(DRY) THEN TW2 = TW2D TA2 = TA2D HA2 = HA1 - CPM*(TA1 - TA2) WA2 = WA1 RHA2 = RHA1 ELSE IF ((WET) .OR. MODE .EQ. 1) THEN C ** WET ** TW2 = TW2W HA2 = HA1 - EPSW*(HA1 - HW1) HSAVG = HA1 - (HA1 - HA2)/(1. - DEXP(-ERAT*NTUO)) IF ((9.473.GT.HCONV(HSAVG) .OR. HCONV(HSAVG).GT.355.137) . .AND. (TIME-OUT(11)) .GT. DELT/2.) THEN CALL MESSAGES(-1,'THE CORRELATION FOR THE SATURATION TEMPERA &TURE USED IN THE COOLING COIL MODEL WAS USED WITH AN ENTHALPY OUTS &IDE OF ITS INTENDED RANGE','WARNING',INFO(1),INFO(2)) IWARN=IWARN+1 OUT(11) = TIME END IF TSAVG = TSAT(HSAVG) TA2 = TSAVG + (TA1 - TSAVG)*DEXP(-ERAT*NTUO) PSYDAT(2) = TA2 PSYDAT(7) = HA2 CALL PSYCHROMETRICS(TIME,INFO,IU,5,0,PSYDAT,EMODE,ISTAT,*105) CALL LINKCK('TYPE52','PSYCHROMETRICS ',1,99) 105 WA2 = PSYDAT(6) RHA2 = PSYDAT(4) ELSE C ** COMBINED ** HA2 = HA1 - RA*CPW(IU)*(TW2 - TW1) TAX = TA1 - EPSD*(TA1 - TWX) HAX2 = HA1 - RA*CPW(IU)*(TW2 - TWX) HAX = HA1 - CPM*(TA1 - TAX) HSW = HAX + (HA2 - HAX)/(1. - DEXP(-(1. - FDRY)*ERAT*NTUO)) IF ((9.473.GT.HCONV(HSW) .OR. HCONV(HSW).GT.355.137) . .AND. (TIME-OUT(11)) .GT. DELT/2.) THEN CALL MESSAGES(-1,'THE CORRELATION FOR THE SATURATION TEMPERA &TURE USED IN THE COOLING COIL MODEL WAS USED WITH AN ENTHALPY OUTS &IDE OF ITS INTENDED RANGE','WARNING',INFO(1),INFO(2)) IWARN=IWARN+1 OUT(11) = TIME END IF TSW = TSAT(HSW) TA2 = TSW + (TAX - TSW)*DEXP(-(1.-FDRY)*ERAT*NTUO) PSYDAT(2) = TA2 PSYDAT(7) = HA2 CALL PSYCHROMETRICS(TIME,INFO,IU,5,0,PSYDAT,EMODE,ISTAT,*106) CALL LINKCK('TYPE52','PSYCHROMETRICS ',1,99) 106 WA2 = PSYDAT(6) RHA2 = PSYDAT(4) ENDIF C----------------------------------------------------------------------------------------------------------------------- C SET THE OUTPUTS FROM THIS MODEL IN SEQUENTIAL ORDER AND GET OUT C****************************** OUTPUTS ***************************** C QCOIL = MW*CPW(IU)*(TW2 - TW1) QLAT = MA*(WA1 - WA2)*HFG(IU) QSENS = QCOIL - QLAT C 500 OUT(1) = TA2 OUT(2) = WA2 OUT(3) = MA OUT(4) = TW2 OUT(5) = MW OUT(6) = QCOIL OUT(7) = QSENS OUT(8) = QLAT OUT(9) = FDRY OUT(10) = 100.*DMIN1(1.,(DMAX1(0.,RHA2))) C----------------------------------------------------------------------------------------------------------------------- C EVERYTHING IS DONE - RETURN FROM THIS SUBROUTINE AND MOVE ON RETURN 1 END C----------------------------------------------------------------------------------------------------------------------- C THIS SUBROUTINE CALCULATES G (WHICH IS BASED ON AN ENERGY BALANCE ON THE C WET SIDE OF THE COIL), FDRY AND EPSD. C SEE VARIABLE DEFINITIONS IN THE MAIN PROGRAM. CPW IS ASSUMED TO BE C PASSED IN THE CORRECT UNITS. C SUBROUTINE CALCG(TW2,TDP,TA1,TW1,HA1,HW1,RO,RA,RE,KF,FT,NTUO, . NTUI,NTUD,HDO,CPM,CSTAR,CPW,EFFD,FRATIO, . INFO,IU,G,FDRY,EPSD,TWX,ERAT,EMODE) IMPLICIT NONE INTEGER*4 INFO DIMENSION PSYDAT(9),INFO(15) DOUBLE PRECISION KF,NTUD,NTUW,NTUO,NTUI,MSTAR,PATM,EXPK,TDP,TW2, & CSTAR,TA1,FDRY,TEMP,C,EPSD,TWX,TW1,PSYDAT,HWX,CS,HW1,RA,HWO,HDO, & CPM,ALPHA,RO,RE,BETA,FT,EFF,FIN,EFFW,FRATIO,ERAT,EFFD,CPW,EPSW,G, & HA1 INTEGER EMODE,IU,ISTAT DATA PATM/1.0/ EXPK = ((TDP-TW2)+CSTAR*(TA1-TDP))/ & (1.-NTUD*(1.-CSTAR)/NTUO)/(TA1-TW2) EXPK = DMAX1(1.0E-6,EXPK) FDRY = -1./NTUD/(1.-CSTAR)*DLOG(EXPK) FDRY = DMAX1(0.D0,DMIN1(1.D0,FDRY)) TEMP = -FDRY*NTUD*(1. - CSTAR) IF (TEMP .LT. 50.0) THEN C = DEXP(TEMP) EPSD = (1. - C)/(1. - CSTAR*C) ELSE EPSD = 1/CSTAR ENDIF TWX = (TW2 - CSTAR*EPSD*TA1)/(1 - CSTAR*EPSD) TWX = DMAX1(TWX,TW1+1.E-04) PSYDAT(1) = PATM PSYDAT(2) = TWX PSYDAT(3) = TWX CALL PSYCHROMETRICS(1.0D0,INFO,IU,1,0,PSYDAT,EMODE,ISTAT,*107) CALL LINKCK('TYPE52','PSYCHROMETRICS ',1,99) 107 HWX = PSYDAT(7) CS = (HWX - HW1)/(TWX - TW1) MSTAR= CS/RA/CPW HWO = CS*HDO/CPM ALPHA = RO/RE BETA = RE*DSQRT(2.*HWO/KF/FT) EFF = FIN(ALPHA,BETA) EFFW = 1. - FRATIO*(1.-EFF) ERAT = EFFW/EFFD NTUW = ERAT*NTUO/(1. + ERAT*NTUO/NTUI*MSTAR) TEMP = -(1. - FDRY)*NTUW*(1. - MSTAR) IF (TEMP .LT. 50.0) THEN C = DEXP(TEMP) EPSW = (1. - C)/(1. - MSTAR*C) ELSE EPSW = 1/MSTAR ENDIF G=TWX*(1.-CSTAR*EPSW*EPSD)-TW1-CSTAR/CPM*EPSW*(HA1-EPSD*CPM*TA1 & -HW1) RETURN END C****************************************************************************** C THIS FUNCTION CALCULATES THE FIN EFFICIENCY (EFFECTIVENESS) C OF AN ANNULAR FIN OF CONSTANT THICKNESS. C C ALPHA = RADIUS AT FIN BASE / RADIUS AT FIN TIP C BETA = RADIUS AT FIN TIP * C (SQRT (2 * CONVECTION COEFFICIENT / C FIN CONDUCTIVITY * FIN THICKNESS)) FUNCTION FIN(ALPHA,BETA) IMPLICIT NONE DOUBLE PRECISION FIN,I0,I1,K0,K1,ALPBET,ALPHA,BETA,XI0,XI1,XK0, & XK1,YI0,YI1,YK0,YK1 ALPBET = ALPHA * BETA CALL BESSEL(ALPBET,I0,I1,K0,K1) XI0 = I0 XI1 = I1 XK0 = K0 XK1 = K1 CALL BESSEL(BETA,I0,I1,K0,K1) YI0 = I0 YI1 = I1 YK0 = K0 YK1 = K1 FIN=2.*ALPHA/BETA/(1.- ALPHA**2)*(XK1*YI1-XI1*YK1)/(XK0*YI1 & +XI0*YK1) RETURN END C**************************************************************************** C THIS SUBROUTINE USES POLYNOMIAL APPROXIMATIONS TO EVALUATE C THE BESSEL FUNCTIONS. THE APPROXIMATIONS ARE FROM ABRAMOWITZ C AND STEGUN, HANDBOOD OF MATHEMATICAL FUNCTIONS, DOVER C PUBLICATIONS, INC., NEW YORK, NY. C SUBROUTINE BESSEL(X,I0,I1,K0,K1) IMPLICIT NONE DOUBLE PRECISION X,I0,I1,K0,K1,IT,A0,A1,A2,A3,A4,A5,A6,B0,B1,B2, & B3,B4,B5,B6,B7,B8,C0,C1,C2,C3,C4,C5,C6,D0,D1,D2,D3,D4,D5,D6,D7, & D8,E0,E1,E2,E3,E4,E5,E6,F0,F1,F2,F3,F4,F5,F6,G0,G1,G2,G3,G4,G5, & G6,H0,H1,H2,H3,H4,H5,H6,T,TT,X1,X2 C THE FOLLOWING DATA STATEMENTS CONTAIN THE COEFFICIENTS TO C THE POLYNOMIALS. C I0 DATA A0/1.0/,A1/3.5156229/,A2/3.0899424/,A3/1.2067492/ DATA A4/0.2659732/,A5/0.0360768/,A6/0.0045813/ C I0 DATA B0/0.39894228/,B1/0.01328592/,B2/0.00225319/ DATA B3/-0.00157565/,B4/0.00916281/,B5/-0.02057706/ DATA B6/0.02635537/,B7/-0.01647633/,B8/0.00392377/ C I1 DATA C0/0.5/,C1/0.87890594/,C2/0.51498869/,C3/0.15084934/ DATA C4/0.02658733/,C5/0.00301532/,C6/0.00032411/ C I1 DATA D0/0.39894228/,D1/-0.03988024/,D2/-0.00362018/ DATA D3/0.00163801/,D4/-0.01031555/,D5/0.02282967/ DATA D6/-0.02895312/,D7/0.01787654/,D8/-0.00420059/ C K0 DATA E0/-0.57721566/,E1/0.4227842/,E2/0.23069756/ DATA E3/0.0348859/,E4/0.00262698/,E5/0.0001075/,E6/0.0000074/ C K0 DATA F0/1.25331414/,F1/-0.07832358/,F2/0.02189568/ DATA F3/-0.01062446/,F4/0.00587872/,F5/-0.0025154/ DATA F6/0.00053208/ C K1 DATA G0/1.0/,G1/0.15443144/,G2/-0.67278579/,G3/-0.18156897/ DATA G4/-0.01919402/,G5/-0.00110404/,G6/-0.00004686/ C K1 DATA H0/1.25331414/,H1/0.23498619/,H2/-0.0365562/ DATA H3/0.01504268/,H4/-0.00780353/,H5/0.00325614/ DATA H6/-0.00068245/ C IF (X .LT. -3.75) THEN CALL MESSAGES(-1,'THE BESSEL FUNCTION CALLED FROM THE COOLING C &OIL SUBROUTINE COULD NOT BE EVALUATED AT THE GIVEN VALUE','FATAL', &-1,-1) RETURN END IF T=X/3.75 TT=T*T C I0 IF (X .LE. 3.75) THEN I0=A0+TT*(A1+TT*(A2+TT*(A3+TT*(A4+TT*(A5+TT*A6))))) ELSE IT=1/T I0=(B0+IT*(B1+IT*(B2+IT*(B3+IT*(B4+IT*(B5+IT*(B6+IT* . (B7+IT*B8))))))))/(DSQRT(X)*DEXP(-X)) END IF C I1 IF (X .LE. 3.75) THEN I1=(C0+TT*(C1+TT*(C2+TT*(C3+TT*(C4+TT*(C5+TT*C6))))))*X ELSE IT=1/T I1=(D0+IT*(D1+IT*(D2+IT*(D3+IT*(D4+IT*(D5+IT*(D6+IT* . (D7+IT*D8))))))))/(DSQRT(X)*DEXP(-X)) END IF C K0 IF (X .LE. 0.0) THEN CALL MESSAGES(-1,'THE BESSEL FUNCTION CALLED FROM THE COOLING C &OIL SUBROUTINE COULD NOT BE EVALUATED AT THE GIVEN VALUE','FATAL', &-1,-1) RETURN END IF X1 = (X/2.)**2 X2 = 2./X IF (X .LE. 2.0) THEN K0=-DLOG(X/2)*I0+E0+X1*(E1+X1*(E2+X1*(E3+X1*(E4+X1* . (E5+X1*E6))))) ELSE K0=(F0+X2*(F1+X2*(F2+X2*(F3+X2*(F4+X2*(F5+X2*F6)))))) . /(DSQRT(X)*DEXP(X)) END IF C K1 IF (X .LE. 2.0) THEN K1=(X*DLOG(X/2.)*I1+G0+X1*(G1+X1*(G2+X1*(G3+X1*(G4+X1* . (G5+X1*G6))))))/X ELSE K1=(H0+X2*(H1+X2*(H2+X2*(H3+X2*(H4+X2*(H5+X2*H6)))))) . /(DSQRT(X)*DEXP(X)) END IF RETURN END