SUBROUTINE TYPE89 (TIME,XIN,OUT,T,DTDT,PAR,INFO,ICNTRL,*) C*********************************************************************** C* Copyright ASHRAE A Toolkit for Primary HVAC System Energy C* Calculation C*********************************************************************** C* PROGRAM: TYPE89 (PISCHIL2) C* C* LANGAGE: FORTRAN 77 C* C* PURPOSE: Parameter identification based on C* reciprocating chiller performances C* in steady-state regime.The pressure C* drop at the compressor exhaust is C* taken into account. C*********************************************************************** C* INPUT VARIABLES C* Ifluid Selection of the refrigerant (-) C* If Ifluid C* =1: Refrigerant 12 C* =2: Refrigerant 134a C* =3: Refrigerant 114 C* =4: Refrigerant 22 C* =5: Refrigerant 502 C* =6: Refrigerant 717 (Ammonia) C* xin(1) (-) C* N Number of available working points (-) C* xin(2) (-) C* Choice If Choice (-) C* =1: the supply water temperature is known for both C* evaporator and condenser C* =2: the water temperature is known at the evaporator C* supply and at the condenser exhaust C* =3: the exhaust water temperature is known for both C* evaporator and condenser C* =4: the water temperature is known at the evaporator C* exhaust and at the condenser supply C* xin(3) (-) C* PevG Guess of the cooling capacity (W) C* xin(4) (kJ/hr) C* PcdG Guess of the heat rejected in the condenser (W) C* xin(5) (kJ/hr) C* AUevi Lower bound of the evaporator heat transfer (W/K) C* coefficient C* xin(6) (kJ/hr/øC) C* AUevf Upper bound of the evaporator heat transfer (W/K) C* coefficient C* xin(7) (kJ/hr/øC) C* AUcdi Lower bound of the condenser heat transfer (W/K) C* coefficient C* xin(8) (kJ/hr/øC) C* AUcdf Upper bound of the condenser heat transfer (W/K) C* coefficient C* xin(9) (kJ/hr/øC) C* dAUev Increment (W/K) C* xin(10) (kJ/hr/øC) C* dAUcd Increment (W/K) C* xin(11) (kJ/hr/øC) C* Acyl Cross section of a cylinder (m**2) C* xin(12) (m**2) C* Csti Lower bound of the parameter Cst which is (-) C* the ratio of the nozzle throat area of a cylinder C* to the section of a cylinder C* xin(13) (-) C* Cstf Upper bound of the parameter Cst (-) C* xin(14) (-) C* dCst Increment applied to the parameter Cst (-) C* xin(15) (-) C* NcFL Number of cylinders used in full load (-) C* xin(16) (-) C* C* C* For the working point considered : C* Pev Cooling capacity (W) C* xin(17) (kJ/hr) C* Pcomp Power consumed by the compressor (W) C* xin(18) (kJ/hr) C* Mfrwev Evaporator water mass flow rate (kg/s) C* xin(19) (kg/hr) C* Mfrwcd Condenser water mass flow rate (kg/s) C* xin(20) (kg/hr) C* Twev Evaporator supply or exhaust water temperature (K) C* according to the value of Choice C* xin(21) (øC) C* Twcd condenser supply or exhaust water temperature (K) C* according to the value of Choice C* xin(22) (øC) C* C* OUTPUT VARIABLES C* AUev Value of the evaporator heat transfer (W/K) C* coefficient C* out(1) (kJ/hr/øC) C* AUcd Value of the condenser heat transfer (W/K) C* coefficient C* out(2) (kJ/hr/øC) C* Losses Constant part of the electromechanical (W) C* losses C* out(3) (kJ/hr) C* Alpha Loss factor allowing to define another (-) C* electromechanical loss which is assumed to C* be proportional to the isentropic power C* out(4) (-) C* Cf Clearance factor of the compressor (-) C* out(5) (-) C* VsFL Geometric displacement of the compressor (m**3/s) C* out(6) (m**3/hr) C* Cst Ratio of the nozzle throat area of a cylinder (-) C* to the section of a cylinder C* out(7) (-) C* Aex Nozzle throat area of a cylinder (m**2) C* out(8) (-) C* SEw Dimensionless standard deviation of the linear (-) C* regression over the power consumed by C* the compressor C* out(9) (-) C* SEv Dimensionless standard deviation of the linear (-) C* regression over the refrigerant volume flow rate C* entering the compressor C* out(10) (-) C* Wis Isentropic compression power calculated for (W) C* the working point considered C* out(11) (kJ/hr) C* Pcomp Power consumed by the compressor for (W) C* the working point considered C* out(12) (kJ/hr) C* pfactor Value of pfactor calculated for (-) C* the working point considered C* out(13) (-) C* V Refrigerant volume flow rate entering the (m**3/s) C* compressor for the working point considered C* out(14) (m**3/hr) C* Wreg Calculated values of the power consumed by the (W) C* compressor obtained by means of the linear regression C* , for the working point considered C* out(15) (kJ/hr) C* Vreg Calculated values of the volume flow rate (m**3/s) C* entering the compressor obtained by means of C* the linear regression, for the working point C* considered C* out(16) (m**3/hr) C* C* WATER PROPERTY C* CpWat Specific heat of liquid water (J/kg/K) C*********************************************************************** C MAJOR RESTRICTIONS: The surrouding heat exchanges are C neglected. The refrigerant leaves the C evaporator and the condenser as saturated C vapor and saturated liquid respectively. C The maximum number of working point is C equal to 200.The compression is assumed C to be isentropic. C Perfect gas properties are used. C C DEVELOPER: Jean Lebrun C Jean-Pascal Bourdouxhe C Marc Grodent C University of LiŠge, Belgium C C DATE: March 1, 1995 C C SUBROUTINES CALLED: PISID2 C TYPE93 (PISSIM2) C ERROR C LINKCK C*********************************************************************** C INTERNAL VARIABLES: C Effev Evaporator effectiveness C Effcd Condenser effectiveness C Tev Array containing the evaporating temperature (K) C calculated for each working point C Tcd Array containing the condensing temperature (K) C calculated for each working point C Pev93 Array containing the cooling capacity (W) C calculated by the routine TYPE93 for each C working point C Pcomp93 Array containing the power consumed by the (W) C compressor calculated by the routine TYPE93 C for each working point C F Function to minimize (W**2) C Fmin Lowest value of the function F (W**2) C C Sum1 and Sum2 are variables used in the calculation of the C function to minimize C Id,pfactorg,Vg and Wisg are storage arrays C*********************************************************************** INTEGER*4 INFO,INFO93,Choice DOUBLE PRECISION XIN,OUT,XIN93,OUT93 REAL Losses,Ifluid,NcFL REAL Mfrwev(200),Mfrwcd(200),Id(9) DIMENSION Pev(200),Pcomp(200),Pev93(200),Pcomp93(200),Twev(200), & Twcd(200),Tev(200),Tcd(200),pfactor(200),V(200), & Wis(200),pfactorg(200),Vg(200),Wisg(200),Wreg(200), & Vreg(200) DIMENSION XIN(22),OUT(16),INFO(15),PAR93(7),XIN93(9),OUT93(7), & INFO93(15) COMMON /LUNITS/ LUR,LUW,IFORM,LUK COMMON /SIM/ TIME0,TFINAL,DELT,IWARN COMMON /STORE/ NSTORE,IAV,S(5000) COMMON /CONFIG/ TRNEDT,PERCOM,HEADER,PRTLAB,LNKCHK,PRUNIT,IOCHEK, & PRWARN INFO(6)=16 INFO93(6)=7 DATA CpWat/4187./ C1*** Number of working points to be treated, refrigerant, Choice, C1*** guesses, boundaries, cylinder area and number of cylinders in C1*** full-load regime ITIME=INT(TIME) Ifluid=SNGL(xin(1)) N=IDINT(xin(2)) N2=2*N Choice=IDINT(xin(3)) PevG=SNGL(xin(4)/3.6) PcdG=SNGL(xin(5)/3.6) AUevi=SNGL(xin(6)/3.6) AUevf=SNGL(xin(7)/3.6) AUcdi=SNGL(xin(8)/3.6) AUcdf=SNGL(xin(9)/3.6) dAUev=SNGL(xin(10)/3.6) dAUcd=SNGL(xin(11)/3.6) Acyl=SNGL(xin(12)) Csti=SNGL(xin(13)) Cstf=SNGL(xin(14)) dCst=SNGL(xin(15)) NcFL=SNGL(xin(16)) IF (ITIME.LE.N) THEN C1*** FIRST PART : reception of the remaining inputs associated with C1*** the working point considered. These inputs are stored in arrays. C2*** REMAINING INPUTS 6 (converted in SI units) C2*** ****************** Pev(ITIME)=SNGL(xin(17)/3.6) Pcomp(ITIME)=SNGL(xin(18)/3.6) Mfrwev(ITIME)=SNGL(xin(19)/3600.) Mfrwcd(ITIME)=SNGL(xin(20)/3600.) Twev(ITIME)=SNGL(xin(21)+273.15) Twcd(ITIME)=SNGL(xin(22)+273.15) ENDIF IF (ITIME.EQ.N) THEN C1*** SECOND PART : identification and printing of the identified C1*** parameters and the dimensionless standard deviations Fmin=999E25 Nev=INT((AUevf-AUevi)/dAUev+1) Ncd=INT((AUcdf-AUcdi)/dAUcd+1) C1*** For each value of the evaporator heat transfer coefficient C1*** we calculate: AUev=AUevi DO 40 J=1,Nev C1*** For each value of the condenser heat transfer coefficient C1*** we calculate: AUcd=AUcdi DO 50 K=1,Ncd C1*** For each working point we calculate: DO 60 L=1,N C1*** Calculate the evaporator and condenser effectiveness Effev=1-EXP(-AUev/(CpWat*Mfrwev(L))) Effcd=1-EXP(-AUcd/(CpWat*Mfrwcd(L))) C1*** Calculate the evaporating and condensing temperature according C1*** to the value of Choice GOTO (63,64,65,66),Choice C2*** Twsuev and Twsucd are known 63 Tev(L)=Twev(L)-Pev(L)/(Effev*CpWat*Mfrwev(L)) Tcd(L)=Twcd(L)+(Pev(L)+Pcomp(L))/(Effcd*CpWat* & Mfrwcd(L)) GOTO 67 C2*** Twsuev and Twexcd are known 64 Tev(L)=Twev(L)-Pev(L)/(Effev*CpWat*Mfrwev(L)) Tcd(L)=Twcd(L)+(Pev(L)+Pcomp(L))/(1/Effcd-1)/ & (CpWat*Mfrwcd(L)) GOTO 67 C2*** Twexev and Twexcd are known 65 Tev(L)=Twev(L)+Pev(L)*(1-1/Effev)/(CpWat*Mfrwev(L)) Tcd(L)=Twcd(L)+(Pev(L)+Pcomp(L))*(1/Effcd-1)/ & (CpWat*Mfrwcd(L)) GOTO 67 C2*** Twexev and Twsucd are known 66 Tev(L)=Twev(L)+Pev(L)*(1-1/Effev)/(CpWat*Mfrwev(L)) Tcd(L)=Twcd(L)+(Pev(L)+Pcomp(L))/(Effcd*CpWat* & Mfrwcd(L)) 67 CONTINUE 60 CONTINUE C1*** Call the subroutine PISID2 in order to calculate the compressor C1*** parameters associated with the values of the heat transfer C1*** coefficients CALL PISID2 (Ifluid,N,Tev,Tcd,Pev,Pcomp,Acyl,Csti,Cstf, & dCst,NcFL,Losses,Alpha,VsFL,Cf,pfactor,V, & Wis,Aex,Cst,SEw) C1*** For each working point, call the subroutine PISSIM2 in order C1*** to calculate the refrigeratind capacity and the power consumed C1*** by the compressor associated with these parameters of the chiller DO 70 L=1,N par93(1)=AUev*3.6 par93(2)=AUcd*3.6 par93(3)=Losses*3.6 par93(4)=Alpha par93(5)=Aex par93(6)=Cf par93(7)=VsFL*3600 xin93(1)=DBLE(Ifluid) xin93(2)=DBLE(Mfrwev(L)*3600.) xin93(3)=DBLE(Mfrwcd(L)*3600.) xin93(4)=Choice xin93(5)=DBLE(Twev(L)-273.15) xin93(6)=DBLE(Twcd(L)-273.15) xin93(7)=DBLE(PevG*3.6) xin93(8)=DBLE(PcdG*3.6) xin93(9)=DBLE(NcFL) CALL TYPE93 (TIME,XIN93,OUT93,T,DTDT,PAR93,INFO93,ICNTRL, & *76) CALL LINKCK('TYPE89','TYPE93',1,99) 76 CONTINUE Pev93(L)=SNGL(out93(2)/3.6) Pcomp93(L)=SNGL(out93(3)/3.6) 70 CONTINUE C2*** Initialization of the variables used in the calculation of C2*** the function to minimize Sum1=0 Sum2=0 DO 80 L=1,N C2*** Up-to-date of the variables used in the calculation of C2*** the function to minimize Sum1=Sum1+((Pev(L)-Pev93(L))/Pev(L))**2 Sum2=Sum2+((Pcomp(L)-Pcomp93(L))/Pcomp(L))**2 80 CONTINUE C1*** Calculate the value of the function to minimize F=Sum1+Sum2 C1*** IF F is lower than the smallest value found so far THEN store C1*** the value of the variables associated with the two heat C1*** transfer coefficients considered IF (F.LT.Fmin) THEN C2*** Storage of the value of the function to minimize Fmin=F C2*** Storage of the chiller parameters Id(1)=AUev Id(2)=AUcd Id(3)=Losses Id(4)=Alpha Id(5)=VsFL Id(6)=Cf Id(7)=Cst Id(8)=Aex Id(9)=SEw C2*** For each working point, storage of the value of pfactor, of the C2*** refrigerant volume flow rate entering the compressor and of the C2*** isentropic compression power DO 90 L=1,N pfactorg(L)=pfactor(L) Vg(L)=V(L) Wisg(L)=Wis(L) 90 CONTINUE ENDIF C1*** Caculate a new value of the condenser heat transfer coefficient AUcd=AUcd+dAUcd 50 CONTINUE C1*** Calculate a new value of the evaporator heat transfer coefficient AUev=AUev+dAUev 40 CONTINUE C1*** Determine the calculated values of the power consumed by the compressor C1*** and the volume flow rate entering the compressor by means of the C1*** linear regressions DO 100 K=1,N Wreg(K)=Id(3)+(1+Id(4))*Wisg(K) Vreg(K)=Id(5)-Id(6)*Id(5)*pfactorg(K) 100 CONTINUE C2*** Calculate the dimensionless standard deviation of the C2*** linear regression over the refrigerant volume flow rate C2*** entering the compressor Slope=-Id(5)*Id(6) dN=SNGL(N) CALL ERROR (dN,pfactorg,Vg,Slope,SEv,*106) CALL LINKCK('TYPE89','ERROR',1,99) 106 CONTINUE C1*** "Printing" of the identified parameters and the dimensionless C1*** standard deviations out(1)=DBLE(Id(1)*3.6) out(2)=DBLE(Id(2)*3.6) out(3)=DBLE(Id(3)*3.6) out(4)=DBLE(Id(4)) out(5)=DBLE(Id(6)) out(6)=DBLE(Id(5)*3600.) out(7)=DBLE(Id(7)) out(8)=DBLE(Id(8)) out(9)=DBLE(Id(9)) out(10)=DBLE(SEv) ENDIF IF ((ITIME.GT.N).AND.(ITIME.LE.N2)) THEN C1*** THIRD PART : "printing" of the isentropic compression power, the C1*** actual and calculated power consumed by the compressor, the C1*** pfactor, the "actual" and calculated refrigerant volume flow rate out(11)=DBLE(Wisg(ITIME-N)*3.6) out(12)=DBLE(Pcomp(ITIME-N)*3.6) out(13)=DBLE(pfactorg(ITIME-N)) out(14)=DBLE(Vg(ITIME-N)*3600.) out(15)=DBLE(Wreg(ITIME-N)*3.6) out(16)=DBLE(Vreg(ITIME-N)*3600.) ENDIF C1*** FOURTH PART: For TIME= 2*N+1 to last simulation time, C1*** everything has been done; so, just wait for the last C1*** simulation time. RETURN 1 END SUBROUTINE PISID2 (Ifluid,N,Tev,Tcd,Pev,Pcomp,Acyl,Csti,Cstf, & dCst,NcFL,Losses,Alpha,VsFL,Cf,pfactor,V,Wis, & Aex,Cst,SEw) C*********************************************************************** C* Copyright ASHRAE A Toolkit for Primary HVAC System Energy C* Calculation C*********************************************************************** C* PROGRAM: PISID2 C* C* LANGUAGE: FORTRAN 77 C* C* PURPOSE: Parameter identification of a reciprocating C* compressor based on reciprocating chiller C* performances in steady-state regime. C* The pressure drop at the compressor C* exhaust is taken into account. C*********************************************************************** C* INPUT VARIABLES C* Ifluid Selection of the refrigerant (-) C* If Ifluid C* =1: Refrigerant 12 C* =2: Refrigerant 134a C* =3: Refrigerant 114 C* =4: Refrigerant 22 C* =5: Refrigerant 502 C* =6: Refrigerant 717 (Ammonia) C* N Number of available working points (-) C* Tev Array containing the evaporating temperature (K) C* for each working point C* Tcd Array containing the condensing temperature (K) C* for each working point C* Pev Array containing the cooling capacity (W) C* for each working point C* Pcomp Array containing the power consumed by the (W) C* compressor for each working point C* Acyl Section of a cylinder (m**2) C* Csti Lower bound of the parameter Cst which is (-) C* the ratio of the nozzle throat area of a cylinder C* to the section of a cylinder C* Cstf Upper bound of the parameter Cst (-) C* dCst Increment applied to the parameter Cst (-) C* NcFL Number of cylinders used in full load (-) C* C* OUTPUT VARIABLES C* Losses Constant part of the electromechanical (W) C* losses C* Alpha Loss factor allowing to define another (-) C* electromechanical loss which is assumed to C* be proportional to the isentropic power C* Cf Clearance factor of the compressor (-) C* VsFL Geometric displacement of the compressor (m**3/s) C* Cst Ratio of the nozzle throat area of a cylinder (-) C* to the section of a cylinder C* Aex Nozzle throat area of a cylinder (m**2) C* SEw Dimensionless standard deviation of the linear (-) C* regression over the power consumed by C* the compressor C* Wis Array containing the isentropic compression (W) C* power calculated for each working point C* pfactor Array containing the value of pfactor (-) C* calculated for each working point C* V Array containing the refrigerant volume (m**3/s) C* flow rate entering the compressor calculated C* for each working point C* C* REFRIGERANT PROPERTIES C* To Reference temperature (K) C* cpliq Mean specific heat in saturated liquid state (J/kg/K) C* hfo Enthalpy of the saturated liquid at the (J/kg) C* reference temperature C* cpvap Mean specific heat at constant pressure (J/kg/K) C* in superheated vapor state for saturation C* temperatures ranging from 253 K to 283 K C* cpvapcd Mean specific heat at constant pressure (J/kg/K) C* in superheated vapor state for saturation C* temperatures ranging from 303 K to 333 K C* hfgb Vaporization enthalpy at standard boiling (J/kg) C* point (101325 Pa) C* Tb Standard boiling temperature (K) C* Tc Critical temperature (K) C* b Coefficient used in the calculation of the (-) C* vaporization enthalpy C* r Gas constant (J/kg/K) C* Zeta Mean compressibility factor for saturation (-) C* temperatures ranging from 253 K to 283 K C* Zetacd Mean compressibility factor for saturation (-) C* temperatures ranging from 303 K to 333 K C* Gamma Mean isentropic coefficient (-) C* Acl First coefficient in the Clausius-Clapeyron (-) C* equation C* Bcl Second coefficient in the Clausius-Clapeyron (K) C* equation C*********************************************************************** C MAJOR RESTRICTIONS: The surrounding heat exchanges are C neglected.The maximum number of working C points is equal to 200. C The compression is assumed to be isentropic. C Perfect gas properties are used. C C DEVELOPER: Jean Lebrun C Jean-Pascal Bourdouxhe C Marc Grodent C University of LiŠge, Belgium C C DATE: March 1, 1995 C C SUBROUTINES CALLED: ERROR C PROPERTY C LINKCK C*********************************************************************** C INTERNAL VARIABLES: C pevap Array containing the evaporating pressure (Pa) C for each working point C pcond Array containing the condensing pressure (Pa) C for each working point C T1 Array containing the temperature at the (K) C evaporator exhaust for each working point C T1p Array containing the temperature (K) C after heating-up for each working point C dhfg Vaporization enthalpy (J/kg) C h1 Enthalpy at the evaporator exhaust (J/kg) C T3 Array containing the temperature at the (K) C condenser exhaust for each working point C h3 Enthalpy at the condenser exhaust (J/kg) C p2 Pressure at point 2 (Pa) C v2p Specific volume at point 2 prime (m**3/kg) C MfrRef Array containing the refrigerant mass (kg/s) C flow rate for each working point C Toler Relative error tolerance (-) C ErrRel Relative error (-) C C T1pp is a variable used in the iterative scheme C Sum1 to Sum8 and Y are variables used in the linear regressions C SEw1,T1p1,AlphaId,LossesId,CstId and AexId are storage variables C*********************************************************************** REAL Losses,LossesId,Ifluid,NcFL REAL MfrRef(200) DIMENSION Pev(200),Pcomp(200),Tev(200),Tcd(200),pfactor(200), & V(200),Wis(200),Y(200),pevap(200),pcond(200),T1(200), & T3(200),T1p(200),T1p1(200) DATA Toler/1E-5/ C2*** Selection of the refrigerant CALL PROPERTY (Ifluid,To,cpliq,hfo,cpvap,cpvapcd,hfgb,Tb,Tc, & b,r,Zeta,Zetacd,Gamma,Acl,Bcl,*1) CALL LINKCK('TYPE89','PROPERTY',1,99) 1 CONTINUE Gm1G=(Gamma-1)/Gamma C1*** For each working points we calculate: DO 20 I=1,N C1*** Calculate the evaporating and the condensing pressures pevap(I)=1000*EXP(Acl+Bcl/Tev(I)) pcond(I)=1000*EXP(Acl+Bcl/Tcd(I)) C1*** Calculate the temperature at the evaporator and C1*** condenser exhaust T1(I)=Tev(I) T3(I)=Tcd(I) C1*** Calculate the enthalpy at the evaporator and C1*** condenser exhaust dhfg=hfgb*((Tc-T1(I))/(Tc-Tb))**b h1=hfo+cpliq*(T1(I)-To)+dhfg h3=hfo+cpliq*(T3(I)-To) C1*** Calculate the refrigerant mass flow rate MfrRef(I)=Pev(I)/(h1-h3) Y(I)=Pcomp(I) 20 CONTINUE SEw1=1E12 J=INT((Cstf-Csti)/dCst+0.5)+1 C1*** For each possible value of the parameter Cst we calculate: DO 30 L=1,J Cst=Csti+FLOAT(L-1)*dCst C2*** Calculate the nozzle throat area of a cylinder Aex=Cst*Acyl C2*** Initialization of the variables used in the linear regression C2*** over the power consumed by the compressor Sum5=0 Sum6=0 Sum7=0 Sum8=0 DO 40 I=1,N C1*** Beginning of the loop C1*** First guess of the temperature after heating-up T1p(I)=T1(I) C2*** Calculate the specific volume at point 2 prime 50 v2p=Zeta*r*T1p(I)/pcond(I)*(pcond(I)/pevap(I))**Gm1G C2*** Calculate the pressure at point 2 p2=pcond(I)+MfrRef(I)**2*v2p/(2*(NcFL*Aex)**2) C2*** Calculate the isentropic compression power Wis(I)=MfrRef(I)*Zeta*r*T1p(I)*((p2/pevap(I))**Gm1G-1)/Gm1G T1pp=T1p(I) C1*** Recalculate the temperature after heating-up T1p(I)=T1(I)+(Pcomp(I)-Wis(I))/(MfrRef(I)*cpvap) ErrRel=ABS((T1p(I)-T1pp)/T1pp) C2*** If converged ,leave loop IF (ErrRel.GT.Toler) GOTO 50 C2*** Up-to-date the variables used in the linear regression C2*** over the power consumed by the compressor Sum5=Sum5+Pcomp(I) Sum6=Sum6+Wis(I) Sum7=Sum7+Wis(I)**2 Sum8=Sum8+Pcomp(I)*Wis(I) 40 CONTINUE C1*** Calculate the parameters Alpha and Losses of the compressor Alpha=(Sum5*Sum6-N*Sum8)/(Sum6**2-N*Sum7)-1 Losses=(Sum5-(1+Alpha)*Sum6)/N C1*** Calculate the dimensionless standard deviation of the linear C1*** regression over the power consumed by the compressor Slope=1+Alpha dN=SNGL(N) CALL ERROR(dN,Wis,Y,Slope,SEw,*56) CALL LINKCK('TYPE89','ERROR',1,99) 56 CONTINUE C1*** Compare this value of the dimensionless standard deviation C1*** with the one associated with the previous of Cst IF (SEw.LT.SEw1) THEN C2*** Storage of the dimensionless standard deviation SEw1=SEw C2*** Storage of the temperature after heating-up for each C2*** working point DO 60 K=1,N T1p1(K)=T1p(K) 60 CONTINUE C2*** Storage of the parameters Alpha,Losses,Cst and Aex C2*** of the compressor AlphaId=Alpha LossesId=Losses CstId=Cst AexId=Aex ENDIF 30 CONTINUE C2*** The parameters characterizing the power consumed by the compressor C2*** are now identified Alpha=AlphaId Losses=LossesId Cst=CstId Aex=AexId SEw=SEw1 C2*** Initialization of the variables used in the linear regression C2*** over the refrigerant volume flow rate entering the compressor Sum1=0 Sum2=0 Sum3=0 Sum4=0 C1*** For each working point we calculate: DO 70 I=1,N C1*** Calculate the refrigerant volume flow rate entering the C1*** compressor associated with the optimal value of Cst V(I)=MfrRef(I)*Zeta*r*T1p1(I)/pevap(I) v2p=Zeta*r*T1p1(I)/pcond(I)*(pcond(I)/pevap(I))**Gm1G p2=pcond(I)+MfrRef(I)**2*v2p/(2*(NcFL*Aex)**2) C1*** Calculate the value of pfactor associated with the optimal C1*** value of Cst pfactor(I)=(p2/pevap(I))**(1/Gamma)-1 C1*** Calculate the isentropic compression power associated with C1*** the optimal value of Cst Wis(I)=MfrRef(I)*Zeta*r*T1p1(I)*((p2/pevap(I))**Gm1G-1)/Gm1G C2*** Up-to-date the variables used in the linear regression over C2*** the refrigerant volume flow rate entering the compressor Sum1=Sum1+pfactor(I) Sum2=Sum2+pfactor(I)**2 Sum3=Sum3+V(I)*pfactor(I) Sum4=Sum4+V(I) 70 CONTINUE C1*** Calculate the two last parameters CfId and VsFLId of the compressor VsFL=(Sum1*Sum3-Sum2*Sum4)/(Sum1**2-N*Sum2) Cf=(N*Sum3-Sum1*Sum4)/(Sum1*Sum3-Sum2*Sum4) RETURN END & b,r,Zeta,Zetacd,Gamma,Acl,Bcl,*) C************************************************************************ C* Copyright ASHRAE A Toolkit for Primary HVAC System Energy C* Calculation C*********************************************************************** C* SUBROUTINE: PROPERTY C* C* LANGUAGE: FORTRAN 77 C* C* PURPOSE: Selection of the thermodynamic properties C* of a given refrigerant. C*********************************************************************** C* INPUT VARIABLES: C* Ifluid Selection of the refrigerant (-) C* If Ifluid C* =1: Refrigerant 12 C* =2: Refrigerant 134a C* =3: Refrigerant 114 C* =4: Refrigerant 22 C* =5: Refrigerant 502 C* =6: Refrigerant 717 (Ammonia) C*********************************************************************** C* OUTPUT VARIABLES: C* To Reference temperature (K) C* cpliq Mean specific heat in saturated liquid state (J/kg/K) C* hfo Enthalpy of the saturated liquid at the (J/kg) C* reference temperature C* cpvap Mean specific heat at constant pressure (J/kg/K) C* in superheated vapor state for saturation C* temperatures ranging from 253 K to 283 K C* cpvapcd Mean specific heat at constant pressure (J/kg/K) C* in superheated vapor state for saturation C* temperatures ranging from 303 K to 333 K C* hfgb Vaporization enthalpy at standard boiling (J/kg) C* point (101325 Pa) C* Tb Standard boiling temperature (K) C* Tc Critical temperature (K) C* b Coefficient used in the calculation of the (-) C* vaporization enthalpy C* r Gas constant (J/kg/K) C* Zeta Mean compressibility factor for saturation (-) C* temperatures ranging from 253 K to 283 K C* Zetacd Mean compressibility factor for saturation (-) C* temperatures ranging from 303 K to 333 K C* Gamma Mean isentropic coefficient (-) C* Acl First coefficient in the Clausius-Clapeyron (-) C* equation C* Bcl Second coefficient in the Clausius-Clapeyron (K) C* equation C*********************************************************************** C MAJOR RESTRICTION: Perfect gas approximation is used C C DEVELOPER: Claudio Saavedra C University of Concepcion, Chile C Marc Grodent, Jean-Pascal Bourdouxhe C University of LiŠge, Belgium C C DATE: March 1, 1995 C*********************************************************************** REAL Ifluid To=233.15 IF (Ifluid.EQ.1) THEN cpliq=917 hfo=0 cpvap=641.6 cpvapcd=779 hfgb=165300 Tb=243.4 Tc=385.2 b=0.37 r=68.7539 Zeta=0.9403 Zetacd=0.8670 Gamma=1.086 Acl=14.669 Bcl=-2443.13 ENDIF IF (Ifluid.EQ.2) THEN cpliq=1265 hfo=0 cpvap=892.5 cpvapcd=1144 hfgb=215100 Tb=246.9 Tc=374.3 b=0.376 r=81.4899 Zeta=0.9411 Zetacd=0.8610 Gamma=1.072 Acl=15.489 Bcl=-2681.99 ENDIF IF (Ifluid.EQ.3) THEN cpliq=925 hfo=0 cpvap=693.6 cpvapcd=784 hfgb=136100 Tb=276.9 Tc=418.9 b=0.359 r=48.6393 Zeta=0.9757 Zetacd=0.9260 Gamma=1.056 Acl=15.107 Bcl=-2908.73 ENDIF IF (Ifluid.EQ.4) THEN cpliq=1144 hfo=0 cpvap=710.4 cpvapcd=936 hfgb=233700 Tb=232.4 Tc=369.2 b=0.369 r=96.1426 Zeta=0.9300 Zetacd=0.8440 Gamma=1.114 Acl=15.070 Bcl=-2421.94 ENDIF IF (Ifluid.EQ.5) THEN cpliq=1090 hfo=0 cpvap=732 cpvapcd=965 hfgb=172500 Tb=227.8 Tc=355.4 b=0.374 r=74.4752 Zeta=0.9130 Zetacd=0.8150 Gamma=1.065 Acl=14.809 Bcl=-2312.21 ENDIF IF (Ifluid.EQ.6) THEN cpliq=4575 hfo=0 cpvap=2447.1 cpvapcd=3159 hfgb=1372900 Tb=239.8 Tc=405.6 b=0.396 r=488.2214 Zeta=0.9570 Zetacd=0.8960 Gamma=1.230 Acl=16.204 Bcl=-2772.39 ENDIF RETURN 1 END SUBROUTINE LINKCK(ENAME1,ENAME2,ILINK,LNKTYP) C*************************************************************************** C THIS SUBROUTINE WAS WRITTEN FOR TRNSYS 14.0 LINK CHECKING - THIS ROUTINE C IS CALLED BY OTHER SUBROUTINES WHEN AN UNLINKED SUBROUTINE HAS BEEN C FOUND. LINKCK IS NEEDED IN ORDER TO AVOID PUTTING COMMON BLOCKS LUNITS C AND CONFIG IN THE TRNSYS TYPES - JWT -- 3/93 C*************************************************************************** COMMON /LUNITS/ LUR,LUW,IFORM,LUK COMMON /CONFIG/ TRNEDT,PERCOM,HEADER,PRTLAB,LNKCHK,PRUNIT,IOCHEK, 1 PRWARN COMMON /SIM/TIME0,TFINAL,DELT,IWARN CHARACTER*1 TRNEDT,PERCOM,HEADER,PRTLAB,LNKCHK,PRUNIT,IOCHEK, 1 PRWARN CHARACTER*6 ENAME1,ENAME2 INTEGER ILINK,LNKTYP C ILINK = 1 --> GENERATE AN ERROR MESSAGE AND STOP TRNSYS C ILINK = 2 --> GENERATE A WARNING BUT DON'T STOP TRNSYS C ILINK = 3 --> TRNSYS HAS FOUND AN UNLINKED TYPE - GENERATE AN ERROR AND C STOP THE PROGRAM C ILINK = 4 --> WARN THE USER THAT A ROUTINE REQUIRES AN EXTERNAL FUNCTION C ENAME1 --> CALLING PROGRAM THAT NEEDED THE UNLINKED FILE C ENAME2 --> FILE THAT WAS NOT FOUND BY ENAME1 SUBROUTINE C LNKTYP --> TYPE NUMBER THAT IS UNLINKED IF((LNKCHK.EQ.'Y').OR.(LNKCHK.EQ.'y')) THEN IF(ILINK.EQ.1) THEN WRITE(LUW,20) 104,ENAME1,ENAME2 WRITE(LUW,15) CALL MYSTOP(104) ELSE IF(ILINK.EQ.2) THEN WRITE(LUW,20) 104,ENAME1,ENAME2 IWARN=IWARN+1 ELSE IF(ILINK.EQ.3) THEN WRITE(LUW,25) 105,LNKTYP,LNKTYP WRITE(LUW,15) CALL MYSTOP(105) ELSE IF(ILINK.EQ.4) THEN WRITE(LUW,35) LNKTYP,ENAME1,ENAME2 IWARN=IWARN+1 ELSE IF(ILINK.EQ.5) THEN WRITE(LUW,40) 105,LNKTYP,LNKTYP WRITE(LUW,15) CALL MYSTOP(105) ELSE WRITE(LUW,30) ENAME1 IWARN=IWARN+1 ENDIF ENDIF 15 FORMAT(//2X,47H*** SIMULATION TERMINATED WITH ERROR STATUS ***/) 20 FORMAT(//,1X,'***** ERROR *****',8X,'TRNSYS ERROR # ',I3,/1X,A6, 1' REQUIRES THE FILE "',A6,'" WHICH WAS CALLED BUT NOT LINKED.',/1X 1,'PLEASE LINK IN THE REQUIRED FILE AND RERUN THE SIMULATION.') 25 FORMAT(//,1X,'***** ERROR *****',8X,'TRNSYS ERROR # ',I3,/1X, 1'TYPE ',I3,' WAS CALLED IN THE TRNSYS INPUT FILE BUT NOT LINKED.', 1/1X,'LINK TYPE ',I3,' BEFORE RUNNING THIS SIMULATION.') 30 FORMAT(/1X,'*****WARNING*****',/1X,'THE LINKCK SUBROUTINE WAS CALL 1ED WITH AN INVALID OPERAND.',/1X,'THE PROGRAM WHICH CALLED LINKCK 1WITH THE IMPROPER OPERAND WAS ',A6,'.',/1X,'PLEASE MAKE SURE THAT 1THE CALLING PROGRAM IS FIXED OR UNLINKED SUBROUTINES MAY ',/1X,'GO 1 UNNOTICED.') 35 FORMAT(/1X,'*****WARNING*****',/1X,'UNIT ',I2,' ',A6,' REQUIRES TH 1E SUBROUTINE ',A6,/1X,'MAKE SURE THAT THIS SUBROUTINE IS LINKED IN 1 TO AVOID PROBLEMS. IT MAY ALREADY BE LINKED IN.',/) 40 FORMAT(//,1X,'***** ERROR *****',8X,'TRNSYS ERROR # ',I3,/1X, 1'TYPE',I3,' WAS CALLED IN THE TRNSYS INPUT FILE BUT NOT LINKED.', 1/1X,'A DUMMY TYPE SUBROUTINE WAS CALLED IN ITS PLACE. PLEASE LINK' 1,/1X,'TYPE',I3,' BEFORE RUNNING THIS SIMULATION OR TURN OFF THE CH 1ECK'/1X,'FOR UNLINKED SUBROUTINES OPTION IN THE CONFIGURATION FILE 1.') RETURN END SUBROUTINE TYPE93 (TIME,XIN,OUT,T,DTDT,PAR,INFO,ICNTRL,*) C************************************************************************ C* Copyright ASHRAE A Toolkit for Primary HVAC System Energy C* Calculation C*********************************************************************** C* SUBROUTINE: TYPE93 (PISSIM2) C* C* LANGUAGE: FORTRAN 77 C* C* PURPOSE: Numerical simulation of a chiller in C* steady-state regime. The routine mainly C* calculates the cooling capacity and C* the compressor consumption for specified C* working conditions. C* The pressure drop at the compressor exhaust C* is taken into account. C*********************************************************************** C* INPUT VARIABLES C* Ifluid Selection of the refrigerant (-) C* If Ifluid C* =1: Refrigerant 12 C* =2: Refrigerant 134a C* =3: Refrigerant 114 C* =4: Refrigerant 22 C* =5: Refrigerant 502 C* =6: Refrigerant 717 (Ammonia) C* xin(1) (-) C* Mfrwev Water mass flow rate in the evaporator (kg/s) C* xin(2) (kg/hr) C* Mfrwcd Water mass flow rate in the condenser (kg/s) C* xin(3) (kg/hr) C* Choice If Choice (-) C* =1: the supply water temperature is known for both C* evaporator and condenser C* =2: the water temperature is known at the evaporator C* supply and at the condenser exhaust C* =3: the exhaust water temperature is known for both C* evaporator and condenser C* =4: the water temperature is known at the evaporator C* exhaust and at the condenser supply C* xin(4) (-) C* Twev1 This value is equal to the evaporator supply or (K) C* exhaust water temperature according to the value C* of Choice C* xin(5) (øC) C* Twcd1 This value is equal to the condenser supply or (K) C* exhaust water temperature according to the value C* of Choice C* xin(6) (øC) C* PevG Guess of the cooling capacity (W) C* xin(7) (kJ/hr) C* PcdG Guess of the heat rejected in the condenser (W) C* xin(8) (kJ/hr) C* NcFL Number of cylinders used in full load (-) C* xin(9) (-) C* C* OUTPUT VARIABLES C* MfrRef Refrigerant mass flow rate (kg/s) C* out(1) (kg/hr) C* Pev Cooling capacity (W) C* out(2) (kJ/hr) C* Pcomp Power consumed by the compressor (W) C* out(3) (kJ/hr) C* Pcd Heat rejected in the condenser (W) C* out(4) (kJ/hr) C* COP Coefficient of performance (-) C* out(5) (-) C* Twev2 This value is equal to the evaporator exhaust (K) C* or supply water temperature according to the value C* of Choice C* out(6) (øC) C* Twcd2 This value is equal to the condenser exhaust (K) C* or supply water temperature according to the value C* of Choice C* out(7) (øC) C* C* PARAMETERS C* AUev Evaporator heat transfer coefficient (W/K) C* par(1) (kJ/øC/hr) C* AUcd Condenser heat transfer coefficient (W/K) C* par(2) (kJ/øC/hr) C* Losses Constant part of the electromechanical losses (W) C* par(3) (kJ/hr) C* Alpha Loss factor allowing to define another (-) C* electromechanical loss which is assumed to be C* proportional to the isentropic power C* par(4) (-) C* Aex Equivalent nozzle throat area of a cylinder (m**2) C* par(5) (m**2) C* Cf Clearance factor of the compressor (-) C* par(6) (-) C* VsFL Geometric displacement of the compressor (m**3/s) C* par(7) (m**3/hr) C* C* WATER PROPERTY C* CpWat Specific heat of liquid water (J/kg/K) C* C* REFRIGERANT PROPERTIES C* To Reference temperature (K) C* cpliq Mean specific heat in saturated liquid state (J/kg/K) C* hfo Enthalpy of the saturated liquid at the (J/kg) C* reference temperature C* cpvap Mean specific heat at constant pressure (J/kg/K) C* in superheated vapor state for saturation C* temperatures ranging from 253 K to 283 K C* cpvapcd Mean specific heat at constant pressure (J/kg/K) C* in superheated vapor state for saturation C* temperatures ranging from 303 K to 333 K C* hfgb Vaporization enthalpy at standard boiling (J/kg) C* point (101325 Pa) C* Tb Standard boiling temperature (K) C* Tc Critical temperature (K) C* b Coefficient used in the calculation of the (-) C* vaporization enthalpy C* r Gas constant (J/kg/K) C* Zeta Mean compressibility factor for saturation (-) C* temperatures ranging from 253 K to 283 K C* Zetacd Mean compressibility factor for saturation (-) C* temperatures ranging from 303 K to 333 K C* Gamma Mean isentropic coefficient (-) C* Acl First coefficient in the Clausius-Clapeyron (-) C* equation C* Bcl Second coefficient in the Clausius-Clapeyron (K) C* equation C*********************************************************************** C MAJOR RESTRICTIONS: The surrounding heat exchanges are C neglected. The refrigerant leaves the C evaporator and the condenser as saturated C vapor and saturated liquid respectively. C The compression is assumed to be isentropic. C Perfect gas properties are used. C C DEVELOPER: Jean Lebrun C Jean-Pascal Bourdouxhe C Marc Grodent C University of LiŠge, Belgium C C DATE: March 1, 1995 C C SUBROUTINE CALLED: PROPERTY C LINKCK C*********************************************************************** C INTERNAL VARIABLES: C Twsuev Evaporator supply water temperature (K) C Twexev Evaporator exhaust water temperature (K) C Twsucd Condenser supply water temperature (K) C Twexcd Condenser exhaust water temperature (K) C T1 Evaporating temperature (K) C T1p Temperature after the heating-up (K) C v1p Specific volume after the heating-up (m**3/kg) C p1 Evaporating pressure (Pa) C dhfg Vaporization enthalpy (J/kg) C h1 Enthalpy at the evaporator exhaust (J/kg) C T3 Condensing temperature (K) C h3 Enthalpy at the condenser exhaust (J/kg) C p3 Condensing pressure (Pa) C v2p Specific volume at point 2' (m**3/kg) C p2 Pressure at point 2 (Pa) C Effvol Volumetric effectiveness of the compressor (-) C Wis Isentropic compression power (W) C Effev Evaporator effectiveness (-) C Effcd Condenser effectiveness (-) C Iter1,Iter2 Loop counters (-) C TolRel Relative error tolerance (-) C ErrRel Relative error (-) C IterMax Iteration maximum (-) C dpex Pressure drop at the compressor exhaust (Pa) C C Pevp,Pcdp,T1pp and dpexp are variables used in the iterative C scheme C*********************************************************************** INTEGER*4 INFO DOUBLE PRECISION XIN,OUT REAL Mfrwev,Mfrwcd,MfrRef,Ifluid,NcFL,Losses DIMENSION PAR(7),XIN(9),OUT(7),INFO(15) COMMON /LUNITS/ LUR,LUW,IFORM,LUK COMMON /SIM/ TIME0,TFINAL,DELT,IWARN COMMON /STORE/ NSTORE,IAV,S(5000) COMMON /CONFIG/ TRNEDT,PERCOM,HEADER,PRTLAB,LNKCHK,PRUNIT,IOCHEK, & PRWARN INFO(6)=7 DATA TolRel,IterMax,CpWat/1E-05,200,4187/ C*** INPUTS 9 (converted in SI units) C************ Ifluid=SNGL(xin(1)) Mfrwev=SNGL(xin(2)/3600.) Mfrwcd=SNGL(xin(3)/3600.) Choice=SNGL(xin(4)) Twev1=SNGL(xin(5)+273.15) Twcd1=SNGL(xin(6)+273.15) PevG=SNGL(xin(7)/3.6) PcdG=SNGL(xin(8)/3.6) NcFL=SNGL(xin(9)) C*** PARAMETERS 7 (converted in SI units) C**************** AUev=par(1)/3.6 AUcd=par(2)/3.6 Losses=par(3)/3.6 Alpha=par(4) Aex=par(5) Cf=par(6) VsFL=par(7)/3600. C2*** Selection of the refrigerant CALL PROPERTY (Ifluid,To,cpliq,hfo,cpvap,cpvapcd,hfgb,Tb,Tc, & b,r,Zeta,Zetacd,Gamma,Acl,Bcl,*1) CALL LINKCK('TYPE93','PROPERTY',1,99) 1 CONTINUE Pev=PevG Pcd=PcdG Gm1G=(Gamma-1)/Gamma NChoice=INT(Choice) GOTO (10,20,30,40),NChoice C2*** The supply water temperature is known for both evaporator C2*** and condenser 10 Twsuev=Twev1 Twsucd=Twcd1 GOTO 50 C2*** The water temperature is known at the evaporator supply and C2*** at the condenser exhaust 20 Twsuev=Twev1 Twexcd=Twcd1 GOTO 50 C2*** The exhaust water temperature is known for both evaporator C2*** and condenser 30 Twexev=Twev1 Twexcd=Twcd1 GOTO 50 C2*** The water temperature is known at the evaporator exhaust and C2*** at the condenser supply 40 Twexev=Twev1 Twsucd=Twcd1 50 CONTINUE C1*** Calculate the evaporator and condenser effectivenesses Effev=1-EXP(-AUev/(CpWat*Mfrwev)) Effcd=1-EXP(-AUcd/(CpWat*Mfrwcd)) C1*** Beginning of the first loop Iter1=0 60 Iter1=Iter1+1 C1*** Calculate the evaporating temperature according to the C1*** information available IF ((NChoice.EQ.1).OR.(NChoice.EQ.2)) THEN T1=Twsuev-Pev/(Effev*CpWat*Mfrwev) ELSE T1=Twexev+(Pev/(CpWat*Mfrwev))*(1-1/Effev) ENDIF C1*** Calculate the evaporating pressure p1=1000*EXP(Acl+Bcl/T1) C1*** Calculate the enthalpy at the evaporator exhaust dhfg=hfgb*((Tc-T1)/(Tc-Tb))**b h1=hfo+cpliq*(T1-To)+dhfg C1*** Beginning of the second loop Iter2=0 70 Iter2=Iter2+1 C1*** Calculate the condensing temperature according to the information C1*** available IF ((NChoice.EQ.1).OR.(NChoice.EQ.4)) THEN T3=Twsucd+Pcd/(Effcd*CpWat*Mfrwcd) ELSE T3=Twexcd+(Pcd/(CpWat*Mfrwcd))*(1/Effcd-1) ENDIF C1*** Calculate the condensing pressure p3=1000*EXP(Acl+Bcl/T3) C2*** Calculate the enthalpy at the condenser exhaust h3=hfo+cpliq*(T3-To) C1*** Beginning of the third loop C2*** First guess of the temperature after the heating-up T1p=T1 C2*** Calculate the specific volumes after the heating-up and at point 2' 80 v1p=Zeta*r*T1p/p1 v2p=Zeta*r*T1p/p3*(p3/p1)**Gm1G C2*** First guess of the pressure drop dpex=1 90 p2=p3+dpex C2*** Calculate the volumetric efficiency of the compressor Effvol=1+Cf-Cf*(p2/p1)**(1/Gamma) C1*** Calculate the refrigerant mass flow rate MfrRef=Effvol*VsFL/v1p dpexp=dpex C2*** Recalculate the pressure drop dpex=MfrRef**2*v2p/(2*(NcFL*Aex)**2) ErrRel=ABS((dpex-dpexp)/dpexp) C2*** If converged, leave loop IF (ErrRel.GT.TolRel) GOTO 90 C2*** Calculate the isentropic compression power Wis=MfrRef*Zeta*r*T1p*((p2/p1)**Gm1G-1)/Gm1G T1pp=T1p C1*** Recalculate the temperature after the heating-up T1p=T1+(Losses+Alpha*Wis)/(MfrRef*cpvap) ErrRel=ABS((T1p-T1pp)/T1pp) C2*** If converged, leave the third loop IF (ErrRel.GT.TolRel) GOTO 80 C1*** Calculate the power consumed by the compressor Pcomp=Losses+Alpha*Wis+Wis Pcdp=Pcd C1*** Calculate the heat rejected in the condenser Pcd=Pev+Pcomp ErrRel=ABS((Pcd-Pcdp)/Pcdp) C2*** If converged, leave the second loop IF ((ErrRel.GT.TolRel).AND.(Iter2.LE.IterMax)) GO TO 70 C1*** Calculate the cooling capacity Pcd=Pcdp Pevp=Pev Pev=MfrRef*(h1-h3) ErrRel=ABS((Pev-Pevp)/Pevp) C2*** If converged, leave the first loop IF ((ErrRel.GT.TolRel).AND.(Iter1.LE.IterMax)) GO TO 60 Pev=Pevp C1*** Calculate the coefficient of performance COP=Pev/Pcomp GOTO (100,110,120,130),NChoice C1*** Calculate the evaporator and condenser exhaust water temperatures 100 Twev2=Twsuev-Pev/(CpWat*Mfrwev) Twcd2=Twsucd+Pcd/(CpWat*Mfrwcd) GOTO 140 C1*** Calculate the evaporator exhaust water temperature and the condenser C1*** supply water temperature 110 Twev2=Twsuev-Pev/(CpWat*Mfrwev) Twcd2=Twexcd-Pcd/(CpWat*Mfrwcd) GOTO 140 C1*** Calculate the evaporator and condenser supply water temperatures 120 Twev2=Twexev+Pev/(CpWat*Mfrwev) Twcd2=Twexcd-Pcd/(CpWat*Mfrwcd) GOTO 140 C1*** Calculate the evaporator supply water temperature and the condenser C1*** exhaust water temperature 130 Twev2=Twexev+Pev/(CpWat*Mfrwev) Twcd2=Twsucd+Pcd/(CpWat*Mfrwcd) 140 CONTINUE C*** OUTPUTS 7 (converted in TRNSYS units) C************* out(1)=DBLE(MfrRef*3600.) out(2)=DBLE(Pev*3.6) out(3)=DBLE(Pcomp*3.6) out(4)=DBLE(Pcd*3.6) out(5)=DBLE(COP) out(6)=DBLE(Twev2-273.15) out(7)=DBLE(Twcd2-273.15) RETURN 1 END SUBROUTINE ERROR (N,X,Y,Slope,ASyx,*) C************************************************************************ C* Copyright ASHRAE A Toolkit for Primary HVAC System Energy C* Calculation C*********************************************************************** C* SUBROUTINE: ERROR C* C* LANGAGE: FORTRAN 77 C* C* PURPOSE: Calculation of the dimensionless standard C* deviation of the linear regressions C*********************************************************************** C* INPUT VARIABLES C* N Number of available working points (-) C* X Array containing either the isentropic (W) or (-) C* compression power or the value of pfactor C* for each working point C* Y Array containing either the power (W) or (m**3/s) C* consumed by the compressor or the C* refrigerant volume flow rate entering the C* compressor for each working point C* Slope Slope of the straight line (-) or (m**3/s) C* C* OUTPUT VARIABLE C* ASyx Dimensionless standard deviation of the (-) C* linear regression C*********************************************************************** C MAJOR RESTRICTION: The maximum number of working points C is equal to 200. C C DEVELOPER: Jean-Pascal Bourdouxhe C Marc Grodent C University of LiŠge, Belgium C C DATE: March 1, 1995 C*********************************************************************** C INTERNAL VARIABLES C MeanX Mean value of either the isentropic (W) or (-) C compression power or the pfactor C MeanY Mean value of either the power (W) or (m**3/s) C consumed by the compressor or the C refrigerant volume flow rate entering C the compressor C C SumX and SumY are variables used to calculate these mean values C Xsq,Ysq and rsq are secondary variables used in the routine C*********************************************************************** REAL N,MeanX,MeanY DIMENSION X(200),Y(200) Nb=INT(N) SumX=0 SumY=0 DO 10 I=1,Nb SumX=SumX+X(I) SumY=SumY+Y(I) 10 CONTINUE C1*** Calculate the mean value of either the isentropic compression C1*** power or the pfactor MeanX=SumX/N C1*** Calculate the mean value of either the power consumed by the C1*** compressor or the refrigerant volume flow rate entering the C1*** compressor MeanY=SumY/N Xsq=0 Ysq=0 DO 20 I=1,Nb Xsq=Xsq+(X(I)-MeanX)**2 Ysq=Ysq+(Y(I)-MeanY)**2 20 CONTINUE rsq=Slope**2*Xsq/Ysq C1*** Calculate the standard deviation of the linear regression Syx=SQRT((1-rsq)*Ysq/(N-2)) C1*** Calculate the dimensionless standard deviation of the C1*** linear regression ASyx=Syx/MeanY RETURN 1 END