SUBROUTINE TYPE67(TIME,XIN,OUT,T,DTDT,PAR,INFO,ICNTRL,*) C********************************************************************** 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 IT WAS MODIFIED BY J.K. CROSS ON 7/29/94 SO THAT IT WOULD MODEL * C COILS HAVING LESS THAN FOUR ROWS. * C********************************************************************** C DOUBLE PRECISION XIN,OUT COMMON /LUNITS/ LUR,LUW,IFORM,LUK COMMON /SIM/ TIME0,TFINAL,DELT,IWARN COMMON /CONFIG/ TRNEDT,PERCOM,HEADER,PRTLAB,LNKCHK,PRUNIT,IOCHEK, . PRWARN LOGICAL DRY, WET, MAXDRY REAL MA, MW, MHUA, KW, MHUW, LT, KT, KF, JF, MCONV REAL NTUD, NTUW, NTUO, NTUI, MSTAR, WPCV, NTUDP, NTUWP REAL LCOIL, LIMITW, LIMITD INTEGER MODE, FMODE , EMODE CHARACTER*1 TRNEDT,PERCOM,HEADER,PRTLAB,LNKCHK,PRUNIT,IOCHEK, . PRWARN INTEGER*4 INFO DIMENSION PAR(16), XIN(5), OUT(9), INFO(15) 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) = INT(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********************** FIRST CALL OF SIMULATION ********************* C C-- CALL TYPECK TO SET UP UNIT AND, IF NECESSARY, PRINT ERROR MESSAGES C IF(INFO(7) .EQ. -1) THEN NP=15 INFO(6)=11 CALL TYPECK(1, INFO, 5, NP, 0) MODE=ROUND(PAR(1)) IF(MODE.NE.1 .AND. MODE.NE.2) CALL TYPECK(-4,INFO,0,0,0) IU=1 FMODE=ROUND(PAR(13)) IF(FMODE.NE.1 .AND. FMODE.NE.2) CALL TYPECK(-4,INFO,0,0,0) 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 RETURN 1 ENDIF C CHECK FOR MORE THAN ONE COIL UNIT IF (IUNIT .EQ. INFO(1)) GOTO 100 IUNIT=INFO(1) C 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 MODE = ROUND(PAR(1)) 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)) ADIST = PAR(14) IF (FMODE .EQ. 2) DE = ADIST CDIST = PAR(15) PSYDAT(1) = PATM C 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 = SQRT(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/ALOG(RO/RI) C 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 WPCV = WATER PUMP CONTROL VARIABLE C 100 TA1 = XIN(1) WA1 = XIN(2) MA = XIN(3) TW1 = XIN(4) MW = XIN(5) WPCV = XIN(6) C C C************************* COIL ANALYSIS ************************* C C TA2 = DRY-BULB TEMPERATURE OF AIR LEAVING COIL C WA2 = ABSOLUTE 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 C C *** NO FLOW *** C 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 C FOLLOWING LINES ADDED TO TYPE 52!!! IF (WPCV .LT. 1.0E-6) THEN MW = 0.0 TA2 = TA1 ENDIF WA2 = WA1 TW2 = TW1 QCOIL = 0.0 QSENS = 0.0 QLAT = 0.0 QDAIR = 0.0 FDRY = 1.0 GOTO 500 END IF C C IF THE WATER PUMP CONTROL VARIABLE IS SET EQUAL TO 0, THEN C THE OUTLET WATER MASS FLOW RATE IS SET EQUAL TO O. ALL OTHER C OUTPUTS ARE SET AS ABOVE. THIS PROGRAM SEGMENT IS NECESSARY C IN ORDER TO USE THIS TYPE WITH "SOLVER 1". C IF (WPCV .LE. 1.0E-6) THEN MW = 0.0 TA2 = TA1 WA2 = WA1 TW2 = TW1 QCOIL = 0.0 QSENS = 0.0 QLAT = 0.0 QDAIR = 0.0 FDRY = 1.0 GOTO 500 END IF C C *** FLOW *** C C ** INSIDE FLUID COEFFICIENT ** C C FIND THE REYNOLDS NUMBER C GW = MW/NTUBES/AX TPROP = CTOF(TW1) IF ((32.0.GT.TPROP .OR. TPROP.GT.200.) .AND. . (TIME-OUT(11)) .GT. DELT/2.) THEN IF((PRWARN.EQ.'Y').OR.(PRWARN.EQ.'y')) THEN WRITE(LUW,2000) INFO(1),INFO(2),TIME,TPROP ENDIF IWARN=IWARN+1 OUT(11) = TIME END IF REI = GW*DI/MHUW(TPROP) C 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. C IF (REI .GE. 3000.) THEN HI = 0.023*KW(IU)/DI*REI**0.8*PRW(TPROP)**0.4 ELSE IF (REI .LE. 2000.) THEN HI = 4.36*KW(IU)/DI ELSE HIT = 0.023*KW(IU)/DI*REI**0.8*PRW(TPROP)**0.4 HIL = 4.36*KW(IU)/DI HI = HIL + (REI - 2000.0) * (HIT - HIL) / 1000.0 END IF C 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. C 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 C FIN EFFICIENCY CALCULATIONS C ALPHA = RO/RE BETA = RE*SQRT(2.*HDO/KF/FT) EFF = FIN(ALPHA,BETA) EFFD = 1. - FRATIO*(1.-EFF) C C ** DETERMINE NTU'S ** C THE AIR SPECIFIC HEAT IS FOR MOIST AIR C 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 C ** ENTHALPIES FOR ENTERING CONDITIONS ** C PSYDAT(2) = TA1 PSYDAT(6) = WA1 CALL PSYCH(TIME,INFO,IU,4,0,PSYDAT,EMODE,ISTAT,*101) CALL LINKCK('TYPE67','PSYCH ',1,99) 101 TDP = PSYDAT(5) HA1 = PSYDAT(7) PSYDAT(2) = TW1 PSYDAT(3) = TW1 CALL PSYCH(TIME,INFO,IU,1,0,PSYDAT,EMODE,ISTAT,*102) CALL LINKCK('TYPE67','PSYCH ',1,99) 102 HW1 = PSYDAT(7) C C 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 C THE CORRELATION FOR THE EFFECTIVENESS OF EACH TUBE PASS AND THE C OVERALL EFFECTIVENESS ARE BASED ON MATERIAL FOUND IN CHAPTER 2 OF C "COMPACT HEAT EXCHANGERS" BY KAYS AND LONDON. C CSTAR = CPM/RA/CPW(IU) NTUD = NTUO/(1. + NTUO/NTUI*CSTAR) NTUDP = NTUD/PAR(2) TEMP = -NTUDP GAMMA = 1.0 - EXP(TEMP) IF (CSTAR .LE. 1.0E-06) THEN EPSDP = GAMMA ELSE TEMP2 = -GAMMA*CSTAR EPSDP = (1.0/CSTAR)*(1.0 - EXP(TEMP2)) ENDIF IF (CSTAR .GE. 9.999E-01) THEN EPSD = (NROWS*EPSDP)/(1.0 + (NROWS - 1)*EPSDP) ELSE EPSD = (((1.0 - EPSDP*CSTAR)/(1.0 - EPSDP))**NROWS - 1.0)/ @ (((1.0 - EPSDP*CSTAR)/(1.0 - EPSDP))**NROWS - 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 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. C TW2W = TW2D ITER = 0 10 ITER = ITER + 1 PSYDAT(2) = TW2W PSYDAT(3) = TW2W CALL PSYCH(TIME,INFO,IU,1,0,PSYDAT,EMODE,ISTAT,*104) CALL LINKCK('TYPE67','PSYCH ',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*SQRT(2.*HWO/KF/FT) EFF = FIN(ALPHA,BETA) EFFW = 1. - FRATIO*(1.-EFF) ERAT = EFFW/EFFD NTUW = ERAT*NTUO/(1. + ERAT*NTUO/NTUI*MSTAR) NTUWP = NTUW/PAR(2) TEMP = -NTUWP GAMMA = 1.0 - EXP(TEMP) IF (MSTAR .LE. 1.0E-06) THEN EPSWP = GAMMA ELSE TEMP2 = -GAMMA*MSTAR EPSWP = (1.0/MSTAR)*(1.0 - EXP(TEMP2)) ENDIF IF (MSTAR .GE. 9.999E-01) THEN EPSW = (NROWS*EPSWP)/(1.0 + (NROWS - 1)*EPSWP) ELSE EPSW = (((1.0 - EPSWP*MSTAR)/(1.0 - EPSWP))**NROWS - 1)/ @ (((1.0 - EPSWP*MSTAR)/(1.0 - EPSWP))**NROWS - MSTAR) ENDIF C USE A SECANT METHOD TO CONVERGE ON TW2 G = TW2W - . (TW1 + EPSW*(HA1 - HW1)/RA/CPW(IU)) IF (ABS(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((AMAX1(1.0E-06,ABS(DGDT))),DGDT) TOLD = TW2W TW2W = TOLD - G2/DGDT IF (ABS(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 C CHECK TO SEE IF COIL COMPLETELY WET C 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 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.) C CK = NTUD*(1. - CSTAR) LIMITW = (TA1*(CSTAR-1.+CK/NTUO) + TDP*(1-CSTAR))/ . (CK/NTUO) LIMITD = (TA1*(CSTAR-(1.-CK/NTUO)*EXP(-CK)) + . TDP*(1.-CSTAR))/(1.-(1.-CK/NTUO)*EXP(-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,PAR) 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,PAR) IF ((GMAX .GT. 0.0 .AND. GMIN .GT. 0.0) .OR. . (GMAX .LT. 0.0 .AND. GMIN .LT. 0.0)) THEN IF (ABS(GMAX) .LT. ABS(GMIN) ) THEN DRY = MAXDRY WET = .NOT. MAXDRY ELSE WET = MAXDRY DRY = .NOT. MAXDRY ENDIF ELSE TW2 = TMIN ITER = -1 C C ITERATION LOOP TO DETERMINE TW2. C 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,PAR) IF ((ABS(G) .GT. TOL .AND. ITER .LT. IMAX) .AND. . (ABS(TMIN - TMAX) .GT. TOL)) THEN C IF (ITER .EQ. 0) THEN SGNMIN = SIGN(1.,G) C C USE THE MAXIMUM TW2 FROM THE WET AND DRY ANALYSIS AS AN INITIAL GUESS C TW2 = AMAX1(TMIN,AMIN1(AMAX1(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((AMAX1(1.0E-06,ABS(DGDT))),DGDT) TOLD = TW2 TW2 = TOLD - G2/DGDT G1 = G2 END IF C C NARROW THE BOUNDS TO THE SOLUTION C 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 C USE A BI-SECTION IF THE NEW VALUE OF TW2 IS OUTSIDE THE BOUNDS C IF (TW2.GE.TMAX .OR. TW2.LE.TMIN) . TW2 = (TMAX + TMIN)/2. GOTO 15 END IF END IF C C * END OF COMBINED COIL ANALYSIS C END IF ENDIF C C ** END OF COIL ANALYSIS C 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 C ** DRY ** IF(DRY) THEN TW2 = TW2D TA2 = TA2D HA2 = HA1 - CPM*(TA1 - TA2) WA2 = WA1 ELSE IF ((WET) .OR. MODE .EQ. 1) THEN C ** WET ** TW2 = TW2W HA2 = HA1 - EPSW*(HA1 - HW1) HSAVG = HA1 - (HA1 - HA2)/(1. - EXP(-ERAT*NTUO)) IF ((9.473.GT.HCONV(HSAVG) .OR. HCONV(HSAVG).GT.355.137) . .AND. (TIME-OUT(10)) .GT. DELT/2.) THEN IF((PRWARN.EQ.'Y').OR.(PRWARN.EQ.'y')) THEN WRITE(LUW,1000) INFO(1),INFO(2),TIME,HCONV(HSAVG) ENDIF IWARN=IWARN+1 OUT(10) = TIME END IF TSAVG = TSAT(HSAVG) TA2 = TSAVG + (TA1 - TSAVG)*EXP(-ERAT*NTUO) PSYDAT(2) = TA2 PSYDAT(7) = HA2 CALL PSYCH(TIME,INFO,IU,5,0,PSYDAT,EMODE,ISTAT,*105) CALL LINKCK('TYPE67','PSYCH ',1,99) 105 WA2 = PSYDAT(6) 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. - EXP(-(1. - FDRY)*ERAT*NTUO)) IF ((9.473.GT.HCONV(HSW) .OR. HCONV(HSW).GT.355.137) . .AND. (TIME-OUT(10)) .GT. DELT/2.) THEN IF((PRWARN.EQ.'Y').OR.(PRWARN.EQ.'y')) THEN WRITE(LUW,1000) INFO(1),INFO(2),TIME,HCONV(HSW) ENDIF IWARN=IWARN+1 OUT(10) = TIME END IF TSW = TSAT(HSW) TA2 = TSW + (TAX - TSW)*EXP(-(1.-FDRY)*ERAT*NTUO) PSYDAT(2) = TA2 PSYDAT(7) = HA2 CALL PSYCH(TIME,INFO,IU,5,0,PSYDAT,EMODE,ISTAT,*106) CALL LINKCK('TYPE67','PSYCH ',1,99) 106 WA2 = PSYDAT(6) ENDIF C 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 C RETURN 1 1000 FORMAT(/10X,'*** WARNING IN UNIT ',I2,' TYPE ',I2,' TIME ', .F10.4,' ***' ./5X,'THE CORRELATION FOR THE SATURATION TEMPERATURE USED IN THE', ./5X,'COOLING COIL WAS USED WITH AN ENTHALPY OF ',F9.3, . ' KJ/KG,',/5X,'WHICH IS OUT OF ITS INTENDED RANGE.') 2000 FORMAT(/10X,'*** WARNING IN UNIT ',I2,' TYPE ',I2,' TIME ', .F10.4,' ***' ./5X,'THE CORRELATIONS FOR THE VISCOSITY AND PRANDTL # OF WATER . USED',/5X,'IN THE COOLING COIL WERE USED WITH A VALUE . OF ',F9.3,' F,',/5X,'WHICH IS OUT OF THEIR INTENDED RANGE.') END C****************************************************************************** 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,PAR) INTEGER*4 INFO DIMENSION PSYDAT(9),INFO(15) REAL KF,NTUD,NTUW,NTUO,NTUI,MSTAR,PAR(16),NTUDP,NTUWP INTEGER EMODE DATA PATM/1.0/ EXPK = ((TDP-TW2)+CSTAR*(TA1-TDP))/ . (1.-NTUD*(1.-CSTAR)/NTUO)/(TA1-TW2) EXPK = AMAX1(1.0E-6,EXPK) FDRY = -1./NTUD/(1.-CSTAR)*ALOG(EXPK) FDRY = AMAX1(0.,AMIN1(1.,FDRY)) FNTUD = FDRY*NTUD FNTUDP= FNTUD/PAR(2) TEMP = -FNTUDP GAMMA = 1.0 - EXP(TEMP) IF (CSTAR .LE. 1.0E-06) THEN EPSDP = GAMMA ELSE TEMP2 = -GAMMA*CSTAR EPSDP = (1.0/CSTAR)*(1.0 - EXP(TEMP2)) ENDIF IF (CSTAR .GE. 9.999E-01) THEN EPSD = (NROWS*EPSDP)/(1.0 + (NROWS - 1)*EPSDP) ELSE EPSD = (((1.0 - EPSDP*CSTAR)/(1.0 - EPSDP))**NROWS - 1)/ @ (((1.0 - EPSDP*CSTAR)/(1.0 - EPSDP))**NROWS - CSTAR) ENDIF TWX = (TW2 - CSTAR*EPSD*TA1)/(1 - CSTAR*EPSD) TWX = AMAX1(TWX,TW1+1.E-04) PSYDAT(1) = PATM PSYDAT(2) = TWX PSYDAT(3) = TWX CALL PSYCH(TIME,INFO,IU,1,0,PSYDAT,EMODE,ISTAT,*107) CALL LINKCK('TYPE67','PSYCH ',1,99) 107 HWX = PSYDAT(7) CS = (HWX - HW1)/(TWX - TW1) MSTAR= CS/RA/CPW HWO = CS*HDO/CPM ALPHA = RO/RE BETA = RE*SQRT(2.*HWO/KF/FT) EFF = FIN(ALPHA,BETA) EFFW = 1. - FRATIO*(1.-EFF) ERAT = EFFW/EFFD NTUW = ERAT*NTUO/(1. + ERAT*NTUO/NTUI*MSTAR) FNTUW = (1 - FDRY)*NTUW FNTUWP= FNTUW/PAR(2) TEMP = -FNTUWP GAMMA = 1.0 - EXP(TEMP) IF (MSTAR .LE. 1.0E-06) THEN EPSWP = GAMMA ELSE TEMP2 = -GAMMA*MSTAR EPSWP = (1.0/MSTAR)*(1.0 - EXP(TEMP2)) ENDIF IF (MSTAR .GE. 9.999E-01) THEN EPSW = (NROWS*EPSWP)/(1.0 + (NROWS - 1)*EPSWP) ELSE EPSW = (((1.0 - EPSWP*MSTAR)/(1.0 - EPSWP))**NROWS - 1)/ @ (((1.0 - EPSWP*MSTAR)/(1.0 - EPSWP))**NROWS - MSTAR) ENDIF G = TWX*(1. - CSTAR*EPSW*EPSD) - TW1 - . CSTAR/CPM*EPSW*(HA1 - EPSD*CPM*TA1 - HW1) RETURN END C C****************************************************************************** C****************************************************************************** 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)) C FUNCTION FIN(ALPHA,BETA) REAL I0,I1,K0,K1 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**************************************************************************** C**************************************************************************** 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) COMMON /LUNITS/ LUR,LUW,IFORM,LUK REAL X,I0,I1,K0,K1,IT C C THE FOLLOWING DATA STATEMENTS CONTAIN THE COEFFICIENTS TO C THE POLYNOMIALS. C 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 WRITE(LUW,100) 164,67,67,X CALL MYSTOP(1001) RETURN END IF T=X/3.75 TT=T*T C C I0 C 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))))))))/(SQRT(X)*EXP(-X)) END IF C C I1 C 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))))))))/(SQRT(X)*EXP(-X)) END IF C C K0 C IF (X .LE. 0.0) THEN WRITE(LUW,100) X CALL MYSTOP(1001) RETURN END IF X1 = (X/2.)**2 X2 = 2./X IF (X .LE. 2.0) THEN K0=-ALOG(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)))))) . /(SQRT(X)*EXP(X)) END IF C C K1 C IF (X .LE. 2.0) THEN K1=(X*ALOG(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)))))) . /(SQRT(X)*EXP(X)) END IF RETURN C C FORMATS C 100 FORMAT(//,1X,'***** ERROR *****',8X,'TRNSYS ERROR # ',I3,/1X, .'UNIT ',I3,' TYPE ',I3,' COOLING COIL',/1X, .'THE BESSEL FUNCTION CALLED FROM THE COOLING COIL SUBROUTINE COULD . NOT BE'/1X,'EVALUATED AT THE GIVEN VALUE OF ',F5.2,'.') END