SUBROUTINE TYPE67(TIME,XIN,OUT,T,DTDT,PAR,INFO,ICNTRL,*) C*********************************************************************** C* Copyright ASHRAE. Toolkit for HVAC System Energy Calculations C*********************************************************************** C* FUNCTION: UAHX C* C* LANGUAGE: FORTRAN 77 C* C* PURPOSE: Calculate the UA of a heat exchanger C* using the effectiveness-Ntu relationships C* given the entering capacity rate and C* temperature of each flow stream, the C* heat transfer rate under these conditions C* and the heat exchanger configuration. C*********************************************************************** C* INPUT VARIABLES Sample Value C* XIN(1) Cap1 Capacity rate of stream 1 (W/C) 5000.00 C* XIN(2) In1 Inlet state of stream 1 (C) 25.00 C* XIN(3) Cap2 Capacity rate of stream 2 (W/C) 4000.00 C* XIN(4) In2 Inlet state of stream 2 (C) 9.00 C* XIN(5) Q Heat transfer rate (W) 4000.00 C* XIN(6) ConfigHX Heat exchanger configuration (-) 1.0 C* 1 - Counterflow C* 2 - Parallel flow C* 3 - Cross flow, both streams unmixed C* 4 - Cross flow, both streams mixed C* 5 - Cross flow, stream 1 unmixed C* 6 - Cross flow, stream 2 unmixed C* C* OUTPUT VARIABLES C* OUT(1) UAHX Overall heat transfer coefficient (W/C) 5753.64 C* OUT(2) ErrStat Error status indicator,0=ok,1=error(-) 0.0 C*********************************************************************** C MAJOR RESTRICTIONS: Models coil using effectiveness Ntu model C C DEVELOPER: Michael J. Brandemuehl, PhD, PE C University of Colorado at Boulder C C DATE: January 1, 1992 C C INCLUDE FILES: None C SUBROUTINES CALLED: HEATEX C FUNCTIONS CALLED: XITERATE C C REVISION HISTORY: None C C REFERENCE: None C*********************************************************************** C INTERNAL VARIABLES C qEstimate Estimated heat transfer in iteration (W) C ua Estimated heat transfer coefficient (W/C) C error Deviation of dependent variable in iteration C icvg Iteration convergence flag C iter Iteration index C itmax Maximum number of iterations C F1,F2 Previous values of error in iteration C X1,X2 Previous values of independent variable in iteration C*********************************************************************** REAL OUT1,OUT2 DOUBLE PRECISION XIN, OUT DIMENSION XIN(6),OUT(2),INFO(15) INTEGER INFO,IOPT,NI,NP,ND CHARACTER*3 YCHECK(6),OCHECK(2) COMMON /LUNITS/LUR,LUW,IFORM,LUK DATA itmax/20/ DATA YCHECK/'NAV','TE1','NAV','TE1','PW2','DM1'/ DATA OCHECK/'NAV','DM1'/ ErrStat = 0 IOPT = -1 NI = 6 !CORRECT NUMBER OF INPUTS NP = 0 !CORRECT NUMBER OF PARAMETERS ND = 0 !CORRECT NUMBER OF DERIVATIVES CAP1 =XIN(1) IN1 =XIN(2) CAP2 =XIN(3) IN2 =XIN(4) Q =XIN(5) ConfigHX=XIN(6) IF (INFO(7).EQ.-1) THEN CALL TYPECK(IOPT,INFO,NI,NP,ND) C CHECKS TO SEE IF USER'S INPUT AND CORRECT NUMBERS MATCH CALL RCHECK(INFO,YCHECK,OCHECK) C CHECKS TO SEE IF INPUT AND OUTPUT UNITS MATCH ENDIF C1*** Check for Q out of range (effectiveness > 1) IF(ABS(Q) .GT. ABS(MIN(Cap1,Cap2)*(In1-In2))) THEN WRITE(LUW,1001) 1001 FORMAT(/1X,'*** ERROR IN SUBROUTINE UAHX ***'/ & 1X,' Given Q is impossible for given inlet states'/) ErrStat = 1 ENDIF C1*** Estimate UAHX ua = ABS(Q/(In1-In2)) C1*** BEGIN LOOP to iteratively calculate UAHX DO 100 iter = 1,itmax C1*** Calculate heat transfer rate for estimated UAHX CALL HEATEX (Cap1,In1,Cap2,In2,ua,ConfigHx,out1,out2) qEstimate = Cap1*(In1-out1) C1*** Calculate new estimate for UAHX error = ABS(qEstimate) - ABS(Q) ua = XITERATE(ua,error,X1,F1,X2,F2,iter,icvg) C1*** If converged, leave loop IF (icvg .EQ. 1) GO TO 110 100 CONTINUE C1*** If not converged after itmax iterations, return error code WRITE(LUW,1005) itmax 1005 FORMAT(/1X,'*** ERROR IN SUBROUTINE UAHX ***'/ & 1X,' UA has not converged after',I2, & ' iterations'/) ErrStat = 1 110 CONTINUE UAHX = ua OUT(1)=UA OUT(2)=ERRSTAT RETURN 1 END SUBROUTINE HEATEX (Cap1,In1,Cap2,In2,UA,ConfigHX,Out1,Out2) C*********************************************************************** C* Copyright ASHRAE. Toolkit for HVAC System Energy Calculations C*********************************************************************** C* SUBROUTINE: HEATEX C* C* LANGUAGE: FORTRAN 77 C* C* PURPOSE: Calculate the outlet states of a simple C* heat exchanger using the effectiveness-Ntu C* method of analysis. C*********************************************************************** C* INPUT VARIABLES C* Cap1 Capacity rate of stream 1 (W/C) C* In1 Inlet state of stream 1 (C) C* Cap2 Capacity rate of stream 2 (W/C) C* In2 Inlet state of stream 2 (C) C* UA Overall heat transfer coefficient (W/C) C* ConfigHX Heat exchanger configuration (-) C* 1 - Counterflow C* 2 - Parallel flow C* 3 - Cross flow, both streams unmixed C* 4 - Cross flow, both streams mixed C* 5 - Cross flow, stream 1 unmixed C* 6 - Cross flow, stream 2 unmixed C* C* OUTPUT VARIABLES C* Out1 Outlet state of stream 1 (C) C* Out2 Outlet state of stream 2 (C) C*********************************************************************** C MAJOR RESTRICTIONS: None C C DEVELOPER: Shauna Gabel C Michael J. Brandemuehl, PhD, PE C University of Colorado at Boulder C C DATE: January 1, 1992 C C INCLUDE FILES: None C SUBROUTINES CALLED: None C FUNCTIONS CALLED: None C C REVISION HISTORY: None C C REFERENCE: Kays, W.M. and A.L. London. 1964. C Compact Heat Exchangers, 2nd Ed., McGraw- C Hill: New York. C*********************************************************************** C* INTERNAL VARIABLES: C* cMin Minimum capacity rate of the streams (W/C) C* cMax Maximum capacity rate of the streams (W/C) C* cRatio Ratio of minimum to maximum capacity rate C* ntu Number of transfer units (-) C* effectiveness Heat exchanger effectiveness (-) C* qMax Maximum heat transfer possible (W) C*********************************************************************** REAL ntu,qMa,large,CAP1,IN1,CAP2,IN2,UA,CONFIGHX, &OUT1,OUT2,CMIN,CMAX,CRATIO,E,EFFECTIVENESS,ETA DATA small/1.E-15/, large/1.E15/ C1*** Ntu and Cmin/Cmax (cRatio) calculations cMin = MIN(Cap1,Cap2) cMax = MAX(Cap1,Cap2) IF( cMax .EQ. 0.) THEN cRatio = 1. ELSE cRatio = cMin/cMax ENDIF IF( cMin .EQ. 0.) THEN ntu = large ELSE ntu = ua/cMin ENDIF C1*** Calculate effectiveness for special limiting cases mode = NINT(ConfigHX) IF(ntu .LE. 0) THEN effectiveness = 0. ELSE IF(cRatio .LT. small) THEN C2*** Cmin/Cmax = 0 and effectiveness is independent of configuration effectiveness = 1 - EXP(-ntu) C1*** Calculate effectiveness depending on heat exchanger configuration ELSE IF (mode .EQ. 1) THEN C2*** Counterflow IF (ABS(cRatio-1.) .LT. small) THEN effectiveness = ntu/(ntu+1.) ELSE e=EXP(-ntu*(1-cRatio)) effectiveness = (1-e)/(1-cRatio*e) ENDIF ELSE IF (mode .EQ. 2) THEN C2*** Parallel flow effectiveness = (1-EXP(-ntu*(1+cRatio)))/(1+cRatio) ELSE IF (mode .EQ. 3) THEN C2*** Cross flow, both streams unmixed eta = ntu**(-0.22) effectiveness = 1 - EXP((EXP(-ntu*cRatio*eta)-1)/(cRatio*eta)) ELSE IF (mode .EQ. 4) THEN C2*** Cross flow, both streams mixed effectiveness = ((1/(1-EXP(-ntu)))+ & (cRatio/(1-EXP(-ntu*cRatio)))-(1/(-ntu)))**(-1) ELSE C2*** One stream is mixed and one is unmixed. Determine whether the C2*** minimum or maximum capacity rate stream is mixed. IF ( (ABS(Cap1-cMin).LT.small .AND. mode.EQ.5) .OR. & (ABS(Cap2-cMin).LT.small .AND. mode.EQ.6) ) THEN C2*** Cross flow, stream with minimum capacity rate unmixed effectiveness = (1-EXP(-cRatio*(1-EXP(-ntu))))/cRatio ELSE C2*** Cross flow, stream with maximum capacity rate unmixed effectiveness = 1-EXP(-(1-EXP(-ntu*cRatio))/cRatio) ENDIF ENDIF C1*** Determine leaving conditions for the two streams qMax = MAX(cMin,small)*(In1-In2) Out1 = In1 - effectiveness*qMax/MAX(Cap1,small) Out2 = In2 + effectiveness*qMax/MAX(Cap2,small) RETURN END REAL FUNCTION XITERATE (X0,F0,X1,F1,X2,F2,ICount,ICvg) C*********************************************************************** C* Copyright ASHRAE. Toolkit for HVAC System Energy Calculations C*********************************************************************** C* C* SUBROUTINE: XITERATE C* C* LANGUAGE: FORTRAN 77 C* C* PURPOSE: Iterately solves for the value of X which C* satisfies F(X)=0. Given Xi,F(Xi) pairs, C* the subroutine tests for convergence and C* provides a new guess for the value of the C* independent variable X. C*********************************************************************** C* INPUT VARIABLES C* F0 Current value of the function F(X) C* X0 Current value of X C* F1,F2 Two previous values of F(Xi) C* X1,X2 Two previous values of X C* C* NOTE: F1,X1,F2,X2 MUST BE STORED AND SAVED IN CALLING C* ROUTINE. THEY NEED NO INITIALIZATION C* C* ICount Number of iterations C* C* OUTPUT VARIABLES C* XIterate New estimate of X for F(X)=0 C* ICvg Convergence flag ICvg = 0: Not converged C* ICvg = 1: Converged C*********************************************************************** C DEVELOPER: Michael J. Brandemuehl, PhD, PE C University of Colorado at Boulder C C DATE: January 1, 1992 C C INCLUDE FILES: None C SUBROUTINES CALLED: None C FUNCTIONS CALLED: None C C REVISION HISTORY: None C C REFERENCE: None C*********************************************************************** C INTERNAL VARIABLES C small Small number used in place of zero C mode Number of points used in fit C mode = 1: Use XPerburb to get new X C mode = 2: Linear equation to get new X C mode > 2: Quadratic equation to get new X C coef(i) Coefficients for quadratic fit C F(X) = coef(1) + coef(2)*X + coef(3)*X*X C check Term under radical in quadratic solution C FiQ,XiQ Double precision values of Fi,Xi C slope Slope for linear fit C tolRel Relative error tolerance C xPerturb Perturbation applied to X to initialize iteration C*********************************************************************** DOUBLE PRECISION coef(3),check,F0Q,F1Q,F2Q,X0Q,X1Q,X2Q DATA tolRel/1.E-5/,xPerturb/0.1/,small/1.E-9/ C1*** Check for convergence by comparing change in X IF ((ABS(X0-X1) .LT. tolRel*MAX(ABS(X0),small) .AND. & ICount .NE. 1) .OR. F0 .EQ. 0.) THEN XIterate = X0 ICvg=1 RETURN ENDIF C1*** Not converged. C2*** If after the second iteration there are enough previous points to C2 fit a quadratic for the new X. If the quadratic fit is not C2 applicable, mode will be set to 1 or 2 and a new X will be C2 determined by incrementing X from xPerturb or from a linear fit. ICvg=0 mode=ICount 10 IF (mode .EQ. 1) THEN C1*** New guess is specified by xPerturb IF (ABS(X0) .GT. small) THEN XIterate = X0*(1.+xPerturb) ELSE XIterate = xPerturb ENDIF ELSEIF (mode .EQ. 2) THEN C1*** New guess calculated from LINEAR FIT of most recent two points SLOPE=(F1-F0)/(X1-X0) IF(slope.EQ.0) THEN mode=1 GO TO 10 ENDIF XIterate=X0-F0/SLOPE ELSE C1*** New guess calculated from QUADRATIC FIT C1*** If two Xi are equal, set mode for linear fit and return to top IF (X0 .EQ. X1) THEN X1=X2 F1=F2 mode=2 GO TO 10 ELSEIF (X0 .EQ. X2) THEN mode=2 GO TO 10 ENDIF C1*** Determine quadratic coefficients from the three data points C1*** using double precision. F2Q=F2 F1Q=F1 F0Q=F0 X2Q=X2 X1Q=X1 X0Q=X0 coef(3)=((F2Q-F0Q)/(X2Q-X0Q)-(F1Q-F0Q)/(X1Q-X0Q))/(X2Q-X1Q) coef(2)=(F1Q-F0Q)/(X1Q-X0Q)-(X1Q+X0Q)*coef(3) coef(1)=F0-(coef(2)+coef(3)*X0Q)*X0Q C1*** If points are colinear, set mode for linear fit and return to top IF (ABS(coef(3)) .LT. 1.D-10) THEN mode=2 GO TO 10 ENDIF C1*** Check for precision. If the coefficients do not accurately C1*** predict the given data points due to round-off errors, set C1*** mode for a linear fit and return to top. IF (ABS((coef(1)+(coef(2)+coef(3)*X1Q)*X1Q-F1Q)/F1Q) .GT. & 1.D-4) THEN mode=2 GO TO 10 ENDIF C1*** Check for imaginary roots. If no real roots, set mode to C1*** estimate new X by simply incrementing by xPerturb check=coef(2)**2-4*coef(1)*coef(3) IF (check .LT. 0) THEN C1*** Imaginary roots -- go back to linear fit mode=2 GO TO 10 ELSEIF (check .GT. 0) THEN C1*** Real unequal roots -- determine root nearest to most recent guess XIterate=(-coef(2)+SQRT(check))/coef(3)/2 xOther=-XIterate-coef(2)/coef(3) IF (ABS(XIterate-X0) .GT. ABS(xOther-X0)) XIterate=xOther ELSE C1*** Real Equal Roots -- one solution XIterate=-coef(2)/coef(3)/2 ENDIF ENDIF C1*** Set previous variable values for the next iteration IF (mode .LT. 3) THEN C1*** No valid previous points to eliminate. X2=X1 F2=F1 X1=X0 F1=F0 ELSE C1*** Eliminate one previous point based on sign and magnitude of F(X) C2*** Keep the current point and eliminate one of the previous ones. IF (F1*F0 .GT. 0 .AND. F2*F0 .GT. 0) THEN C2*** All previous points of same sign. Eliminate one with biggest F(X) IF (ABS(F2) .GT. ABS(F1)) THEN X2=X1 F2=F1 ENDIF ELSE C1*** Points of different sign. C1*** Eliminate the previous one with the same sign as current F(X). IF (F2*F0 .GT. 0) THEN X2=X1 F2=F1 ENDIF ENDIF X1=X0 F1=F0 ENDIF RETURN END