SUBROUTINE TYPE56 (TIME,XIN,OUT,T,DTDT,PAR,INFO,ICNTRL,*) C*********************************************************************** C* Copyright ASHRAE. Toolkit for HVAC System Energy Calculations C*********************************************************************** C* SUBROUTINE: MIXIAIR C* C* LANGUAGE: FORTRAN 77 C* C* PURPOSE: Calculate the mass flow rate of two C* entering air streams of a mixing box with C* a known leaving mass flow rate and the C* temperatures of all the streams. C*********************************************************************** C* INPUT VARIABLES DISCRIPTION(UNITS) SAMPLE VALUES C* XIN(1) MLvg Dry air mass flow rate of mixed stream (kg/s) 3.65 C* XIN(2) TLvg Temperature of mixed stream (C) 11.7877 C* XIN(3) T1En Entering temperature of stream 1 (C) 1.67 C* XIN(4) W1En Entering humidity ratio of stream 1 (-) .0017 C* XIN(5) T2En Entering temperature of stream 2 (C) 23.89 C* XIN(6) W2Ent Entering humidity ratio of stream 2 (-) .0092 C* C* OUTPUT VARIABLES C* OUT(1) M1Ent Dry air mass flow rate of stream 1 (kg/s) 2.00044 C* OUT(2) M2Ent Dry air mass flow rate of stream 2 (kg/s) 1.64956 C* OUT(3) WLvg Humidity ratio of mix air stream (-) .00508950 C* OUT(4) ErrStat Error flag (0=ok, 1=error) (-) 0.0 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 REQUIRED: MIXOAIR C FUNCTIONS REQUIRED: XITERATE C C REVISION HISTORY: None C C REFERENCE: None C*********************************************************************** C INTERNAL VARIABLES C deltaT Temperature difference of entering streams (C) C mEst Estimate of leaving flow (kg/s) C small Small number, in place of zero C target Target for mixed air temperatyure (C) C error Deviation of dependent variable in iteration C iter Iteration counter C icvg Iteration convergence flag C F1,F2 Previous values of dependent variable in XITERATE C X1,X2 Previous values of independent variable in XITERATE C*********************************************************************** DOUBLE PRECISION XIN, OUT INTEGER INFO,IOPT,NI,NP,ND REAL mEst,CPAIR,CPVAP,HFG DIMENSION XIN(6),OUT(4),INFO(15) CHARACTER*3 YCHECK(6),OCHECK(4) COMMON /LUNITS/LUR,LUW,IFORM,LUK DATA small/1.E-9/, itmax/20/ DATA CPAIR/1006.0/,CPVAP/1805.0/,HFG/250100.0/ DATA YCHECK/'MF2','TE1','TE1','DM1','TE1','DM1'/ DATA OCHECK/'MF2','MF2','DM1','DM1'/ ErrStat = 0 IOPT = -1 NI = 6 !CORRECT NUMBER OF INPUTS NP = 0 !CORRECT NUMBER OF PARAMETERS ND = 0 !CORRECT NUMBER OF DERIVATIVES MLVG = XIN(1) TLvg = XIN(2) T1Ent = XIN(3) W1Ent = XIN(4) T2Ent = XIN(5) W2Ent = XIN(6) IF (INFO(7).EQ.-1) THEN CALL TYPECK(IOPT,INFO,NI,NP,ND) C CHECKS TO SEE IF USER'S INFO MATCHES CORRECT NUMBERS CALL RCHECK(INFO,YCHECK,OCHECK) C CHECKS TO SEE IF INPUT AND OUTPUT UNITS MATCH ENDIF deltaT = T2Ent-T1Ent IF(ABS(deltaT).LT.small) deltaT=small C1*** Estimate the mass flow rate of stream 1 from a temperature balance M1Ent = (T2Ent-TLvg)/deltaT * MLvg C1*** Set iteration loop parameters target = TLvg C1*** BEGIN LOOP DO 100 iter = 1 ,itmax C1*** Calculate leaving air temperature and humidity for estimated flows M2Ent = MLvg-M1Ent CALL MIXOAIR (CPAIR,CPVAP,HFG,M1Ent,T1Ent,W1Ent,M2Ent, & T2Ent,W2Ent,mEst,TLvg,WLvg,ErrStat) C1*** Compare given leaving air temperature with estimated temperature C1*** and determine new estimate of flow error = TLvg-target M1Ent = XITERATE(M1Ent,error,X1,F1,X2,F2,iter,icvg) C1*** If converged, leave loop and RETURN IF (icvg .EQ. 1) GO TO 999 100 CONTINUE C1*** If not converged after itmax iterations, return error code WRITE(LUW,1005) itmax 1005 FORMAT(/1X,'*** ERROR IN SUBROUTINE MIXIAIR ***'/ & 1X,' Temperature has not converged after, 'I2, & ' iterations'/) ErrStat = 1 999 continue OUT(1)=M1Ent OUT(2)=M2Ent OUT(3)=WLVG OUT(4)=ERRSTAT RETURN 1 END SUBROUTINE MIXOAIR (CPAIR,CPVAP,HFG,M1Ent,T1Ent,W1Ent,M2Ent, & T2Ent,W2Ent,MLvg,TLvg,WLvg,ErrStat) C*********************************************************************** C* Copyright ASHRAE. Toolkit for HVAC System Energy Calculations C*********************************************************************** C* SUBROUTINE: MIXOAIR C* C* LANGUAGE: FORTRAN 77 C* C* PURPOSE: Calculate the leaving temperature, C* humidity ratio and mass flow rate of two C* mixed air streams by simple conservation. C*********************************************************************** C* INPUT VARIABLES C* M1Ent Dry air mass flow rate of stream 1 (kg/s) C* T1Ent Entering temperature of stream 1 (C) C* W1Ent Entering humidity ratio of stream 1 (-) C* M2Ent Dry air mass flow rate of stream 2 (kg/s) C* T2Ent Entering temperature of stream 2 (C) C* W2Ent Entering humidity ratio of stream 2 (-) C* C* OUTPUT VARIABLES C* MLvg Dry air mass flow rate of mixed stream (kg/s) C* TLvg Temperature of mixed stream (C) C* WLvg Humidity ratio of mixed stream (C) C* ErrStat Error flag (0=ok, 1=error) (-) C*********************************************************************** C MAJOR RESTRICTION: None C C DEVELOPER: Shauna Gabel, MS 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: ENTHALPY C DRYBULB C FUNCTIONS REQUIRED: None C C REVISION HISTORY: None C C REFERENCE: None C*********************************************************************** C* INTERNAL VARIABLES C* small Small number used in place of zero C*********************************************************************** INTEGER ErrStat REAL M1Ent,M2EnT,MLvg DATA small/1.E-9/ ErrStat = 0 C1*** Calculate the mass flow rate of the mixed stream. MLvg = M1Ent+M2Ent C1*** If leaving flow is zero, set leaving conditions to those of C1 stream 1 and RETURN. IF (ABS(MLvg) .LE. small) THEN WLvg = W1Ent TLvg = T1Ent ELSE C1*** Leaving flow is not zero. Proceed with calculations. C1*** Calculate the humidity ratio of the mixed stream WLvg = (M1Ent*W1Ent+M2Ent*W2Ent)/MLvg C1*** Calculate the mixed stream temperature from enthalpy and humidity h1Ent = ENTHALPY(CPAIR,CPVAP,HFG,T1Ent,W1Ent) h2Ent = ENTHALPY(CPAIR,CPVAP,HFG,T2Ent,W2Ent) hLvg = (M1Ent*h1Ent+M2Ent*h2Ent)/MLvg TLvg = DRYBULB(CPAIR,CPVAP,HFG,hLvg,WLvg) ENDIF RETURN END REAL FUNCTION DRYBULB (CPAIR,CPVAP,HFG,H,W) C*********************************************************************** C* Copyright ASHRAE. Toolkit for HVAC System Energy Calculations C*********************************************************************** C* FUNCTION: DRYBULB C* C* LANGUAGE: FORTRAN 77 C* C* PURPOSE: Calculate the dry bulb temperature of C* moist air from enthalpy and humidity. C*********************************************************************** C* INPUT VARIABLES: C* H Enthalpy (J/kg) C* W Humidity ratio (-) C* C* OUTPUT VARIABLES: C* Drybulb Dry bulb temperature (C) C* C* PROPERTIES: C* CpAir Specific heat of air (J/kg C) C* CpVap Specific heat of water vapor (J/kg C) C* Hfg Reference heat of vaporization of water (J/kg) C*********************************************************************** C MAJOR RESTRICTIONS: Uses perfect gas relationships C Fit for enthalpy of saturated water vapor 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: PROP.INC C SUBROUTINES CALLED: None C FUNCTIONS CALLED: None C C REVISION HISTORY: None C C REFERENCE: 1989 ASHRAE Handbook - Fundamentals C*********************************************************************** C1*** Calculate the dry bulb temperature as a function of enthalpy and C1*** humidity ratio. C2*** hDryAir = Prop(CpAir)*TDB C2*** hSatVap = Prop(Hfg) + Prop(CpVap)*TDB C2*** Enthalpy = hDryAir + W*hSatVap Drybulb = (H-Hfg*W)/(CpAir+CpVap*W) RETURN END REAL FUNCTION ENTHALPY (CPAIR,CPVAP,HFG,TDB,W) C*********************************************************************** C* Copyright ASHRAE. Toolkit for HVAC System Energy Calculations C*********************************************************************** C* FUNCTION: ENTHALPY C* C* LANGUAGE: FORTRAN 77 C* C* PURPOSE: Calculate the enthalpy of moist air. C*********************************************************************** C* INPUT VARIABLES: C* TDB Dry bulb temperature (C) C* W Humidity ratio (-) C* C* OUTPUT VARIABLES: C* Enthalpy Enthalpy of moist air (J/kg) C* C* PROPERTIES: C* CpAir Specific heat of air (J/kg C) C* CpVap Specific heat of water vapor (J/kg C) C* Hfg Reference heat of vaporization of water (J/kg) C*********************************************************************** C MAJOR RESTRICTIONS Uses perfect gas relationships C Fit for enthalpy of saturated water vapor 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: PROP.INC C SUBROUTINES CALLED: None C FUNCTIONS CALLED: None C C REVISION HISTORY: None C C REFERENCE: 1989 ASHRAE Handbook - Fundamentals C*********************************************************************** C1*** Calculate the enthalpy as a function of dry bulb temperature and C1*** humidity ratio. hDryAir = CpAir*TDB hSatVap = Hfg + CpVap*TDB Enthalpy = hDryAir + W*hSatVap 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