SUBROUTINE TYPE74 (TIME,XIN,OUT,T,DTDT,PAR,INFO,ICNTRL,*) C*********************************************************************** C* Copyright ASHRAE A Toolkit for Primary HVAC System Energy C* Calculation C*********************************************************************** C* PROGRAM: TYPE74 (CENTCHID) C* C* LANGAGE: FORTRAN 77 C* C* PURPOSE: Parameter identification based on C* centrifugal chiller performances in C* steady-state regime 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* C* C* For the working point considered : C* Pev Cooling capacity (W) C* xin(6) (kJ/hr) C* Pcomp Power consumed by the compressor (W) C* xin(7) (kJ/hr) C* Mfrwev Evaporator water mass flow rate (kg/s) C* xin(8) (kg/hr) C* Mfrwcd Condenser water mass flow rate (kg/s) C* xin(9) (kg/hr) C* Twev Evaporator supply or exhaust water temperature (K) C* according to the value of Choice C* xin(10) (øC) C* Twcd condenser supply or exhaust water temperature (K) C* according to the value of Choice C* xin(11) (ø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 internal power C* out(4) (-) C* A Impeller exhaust area (m**2) C* out(5) (m**2) C* U Peripheral speed of the impeller (m/s) C* out(6) (m/s) C* Beta Angle between the direction of the vanes at the (-) C* impeller exhaust and the plan tangent to the C* circumference of the impeller C* (Beta is set equal to 130ø or 2.269 rad) C* out(7) (-) C* SEw Dimensionless standard deviation of the linear (-) C* regression over the power consumed by C* the compressor C* out(8) (-) C* ErrDetec This variable is equal to 1 if the routine does (-) C* not converge. Otherwise ErrDetec is equal to 0. C* out(9) (-) C* Win Internal compression power calculated for (W) C* the working point considered C* out(10) (kJ/hr) C* Pcomp Power consumed by the compressor for (W) C* the working point considered C* out(11) (kJ/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(12) (kJ/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 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. C The impeller and the diffuser are 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: CCHIFLSI C CENCOMID C PROPERTY 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 pratVt0 Limit pressure ratio (-) C Vpratt1 Limit volume flow rate (m**3/s) C Sum1 and Sum2 are the variables used in the calculation C of the function to minimize. C Fmin is a storage variable C Id is a storage array C*********************************************************************** INTEGER*4 INFO,INFO77,Choice DOUBLE PRECISION XIN,OUT,XIN77,OUT77 REAL Losses,MfrRef,Ifluid REAL Mfrwev(200),Mfrwcd(200),Id(2) DIMENSION Pev(200),Pcomp(200),Twev(200),Twcd(200),Tev(200), & Tcd(200),Win(200),Wreg(200) DIMENSION XIN(11),OUT(12),INFO(15),PAR77(7),XIN77(8), & OUT77(8),INFO77(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)=12 INFO77(6)=8 DATA CpWat/4187./ C1*** Number of working points to be treated, refrigerant, Choice, C1*** and guesses. 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) 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(6)/3.6) Pcomp(ITIME)=SNGL(xin(7)/3.6) Mfrwev(ITIME)=SNGL(xin(8)/3600.) Mfrwcd(ITIME)=SNGL(xin(9)/3600.) Twev(ITIME)=SNGL(xin(10)+273.15) Twcd(ITIME)=SNGL(xin(11)+273.15) ENDIF IF (ITIME.EQ.N) THEN C1*** SECOND PART : identification and printing of the identified C1*** parameters and the dimensionless standard deviations C1*** 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('TYPE74','PROPERTY',1,99) 1 CONTINUE C1*** Calculate the value of the heat transfer coefficients SumQev=0 SumQcd=0 SumMwev=0 SumMwcd=0 DO 30 I=1,N SumQev=SumQev+Pev(I) SumQcd=SumQcd+(Pev(I)+Pcomp(I)) SumMwev=SumMwev+Mfrwev(I) SumMwcd=SumMwcd+Mfrwcd(I) 30 CONTINUE QevMean=SumQev/N QcdMean=SumQcd/N MfrwevMean=SumMwev/N MfrwcdMean=SumMwcd/N EffevMean=1/(1+CpWat*MfrwevMean*5/QevMean) EffcdMean=1/(1+CpWat*MfrwcdMean*5/QcdMean) AUev=1000*DFLOAT(INT(-CpWat*MfrwevMean*LOG(1-EffevMean) & /1000)) AUcd=1000*DFLOAT(INT(-CpWat*MfrwcdMean*LOG(1-EffcdMean) & /1000)) C1*** For each working point we calculate: DO 40 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 (43,44,45,46),Choice C2*** Twsuev and Twsucd are known 43 Tev(L)=Twev(L)-Pev(L)/(Effev*CpWat*Mfrwev(L)) Tcd(L)=Twcd(L)+(Pev(L)+Pcomp(L))/(Effcd*CpWat & *Mfrwcd(L)) GOTO 47 C2*** Twsuev and Twexcd are known 44 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 47 C2*** Twexev and Twexcd are known 45 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 47 C2*** Twexev and Twsucd are known 46 Tev(L)=Twev(L)+Pev(L)*(1-1/Effev)/(CpWat & *Mfrwev(L)) Tcd(L)=Twcd(L)+(Pev(L)+Pcomp(L))/(Effcd*CpWat & *Mfrwcd(L)) 47 CONTINUE 40 CONTINUE C1*** Call the subroutine CENCOMID in order to calculate the C1*** parameters Losses and Alpha associated with the values C1*** of the heat transfer coefficients CALL CENCOMID(Ifluid,N,Tev,Tcd,Pev,Pcomp,Losses,Alpha, & Win,SEw,Vmax,pratMax,T1pMax,ErrDetec) IF (ErrDetec.EQ.1) GOTO 999 C1*** For each working point, call the routine CCHIFLSI in order C1*** to calculate the cooling capacity and the power C1*** consumed by the compressor associated with each value C1*** of the parameters A and U. Iter=0 Fmin=999E25 Beta=2.269 C2*** For each value of the limit pressure ratio, we calculate: pratVt0=1.6*pratMax DO 50 I=1,4 C2*** For each value of the limit volume flow rate, we calculate: Vpratt1=1.4*Vmax DO 60 J=1,4 C1*** Calculate the value of the parameters A and U associated C1*** with these limit values Gm1G=(Gamma-1)/Gamma U=((pratVt0**Gm1G-1)/Gm1G*Zeta*r*T1pMax)**0.5 A=-Vpratt1/(U*TAN(Beta)) C2*** Initialisation of the variables used in the calculation C2*** of the function to minimize Sum1=0 Sum2=0 C2*** For each working point, we calculate: DO 70 K=1,N par77(1)=AUev*3.6 par77(2)=AUcd*3.6 par77(3)=Losses*3.6 par77(4)=Alpha par77(5)=U par77(6)=A par77(7)=Beta xin77(1)=DBLE(Ifluid) xin77(2)=DBLE(Mfrwev(K)*3600.) xin77(3)=DBLE(Mfrwcd(K)*3600.) xin77(4)=Choice xin77(5)=DBLE(Twev(K)-273.15) xin77(6)=DBLE(Twcd(K)-273.15) xin77(7)=DBLE(PevG*3.6) xin77(8)=DBLE(PcdG*3.6) CALL TYPE77 (TIME,XIN77,OUT77,T,DTDT,PAR77,INFO77,ICNTRL,*72) CALL LINKCK('TYPE74','TYPE77 ',1,99) 72 CONTINUE MfrRef=SNGL(out77(1)/3600) PevS=SNGL(out77(2)/3.6) PcompS=SNGL(out77(3)/3.6) Pcd=SNGL(out77(4)/3.6) COP=SNGL(out77(5)) Twevap=SNGL(out77(6)+273.15) Twcond=SNGL(out77(7)+273.15) ErrDetec=SNGL(out77(8)) IF (ErrDetec.EQ.1) GOTO 999 C2*** Up-to-date of the variables used in the calculation C2*** of the function to minimize Sum1=Sum1+((PevS-Pev(K))/Pev(K))**2 Sum2=Sum2+((PcompS-Pcomp(K))/Pcomp(K))**2 70 CONTINUE C2*** Calculate the value of the function to minimize F=Sum1+Sum2 C1*** IF F is lower than the smallest value found so far THEN C1*** store these values of the variables A and U IF (F.LT.Fmin) THEN C2*** Storage of the value of the function to minimize Fmin=F C2*** Storage of the considered values of the parameters A and U Id(1)=A Id(2)=U ENDIF C1*** Calculate a new value of the limit volume flow rate Vpratt1=Vpratt1+0.2*Vmax 60 CONTINUE C1*** Calculate a new value of the limit pressure ratio pratVt0=pratVt0+0.2*pratMax 50 CONTINUE C1*** Determine the calculated values of the power consumed by the compressor C1*** by means of the linear regression DO 75 K=1,N Wreg(K)=Losses+(1+Alpha)*Win(K) 75 CONTINUE C1*** "Printing" of the identified parameters and the dimensionless C1*** standard deviation out(1)=DBLE(AUev*3.6) out(2)=DBLE(AUcd*3.6) out(3)=DBLE(Losses*3.6) out(4)=DBLE(Alpha) out(5)=DBLE(Id(1)) out(6)=DBLE(Id(2)) out(7)=DBLE(Beta) out(8)=DBLE(SEw) ENDIF IF ((ITIME.GT.N).AND.(ITIME.LE.N2)) THEN C1*** THIRD PART : "printing" of the internal compression power, the C1*** actual and calculated power consumed by the compressor out(10)=DBLE(Win(ITIME-N)*3.6) out(11)=DBLE(Pcomp(ITIME-N)*3.6) out(12)=DBLE(Wreg(ITIME-N)*3.6) 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. C2*** Remaining output 1 (the error status indicator) C2*** ****************** 999 out(9)=DBLE(ErrDetec) RETURN 1 END SUBROUTINE CENCOMID(Ifluid,N,Tev,Tcd,Pev,Pcomp,Losses,Alpha, & Win,SEw,Vmax,pratMax,T1pMax,ErrDetec) C*********************************************************************** C* Copyright ASHRAE A Toolkit for Primary HVAC System Energy C* Calculation C*********************************************************************** C* PROGRAM: CENCOMID C* C* LANGUAGE: FORTRAN 77 C* C* PURPOSE: Part of the parameter identification of a C* centrifugal compressor based on centrifugal C* chiller performances in steady-state regime 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* 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 internal power C* Win Array containing the internal compression (W) C* power calculated for each working point C* Vmax Maximum value of the refrigerant volume (m**3/s) C* flow rate at the impeller supply C* pratMax Maximum value of the ratio of the condensing (-) C* pressure to the evaporating pressure C* T1pMax Value of the refrigerant temperature after the (K) C* heating-up associated with the value of pratMax C* ErrDetec This variable is equal to 1 if the routine does (-) C* not converge. Otherwise ErrDetec is equal to 0. 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 maximum number of working C points is equal to 200. C The impeller and the diffuser are 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: ERROR C PROPERTY C LINKCK C*********************************************************************** C INTERNAL VARIABLES: C pevap Evaporating pressure (Pa) C pcond Condensing pressure (Pa) C T1 Evaporating temperature (K) C T1p Temperature after heating-up (K) 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 MfrRef Refrigerant mass flow rate (kg/s) C pratio Ratio of the condensing pressure to the (-) C evaporating pressure C v1p Specific volume after the heating-up (m**3/kg) C Vfr1p Volume flow rate after the heating-up (m**3/s) C TolRel Relative error tolerance (-) C ErrRel Relative error (-) C C T1pp is a variable used in the iterative scheme C Sum1 to Sum4 are variables used in the linear regression C Z is an array used in the linear regression over the power C consumed by the compressor C*********************************************************************** REAL Losses,MfrRef,Ifluid DIMENSION Pev(200),Pcomp(200),Tev(200),Tcd(200),Win(200),Z(200) DATA TolRel,IterMax/1E-4,200/ C1*** 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('TYPE74','PROPERTY',1,99) 1 CONTINUE Gm1G=(Gamma-1)/Gamma ErrDetec=0 SEvmin=1E25 Vmin=1E25 Vmax=0 C2*** Initialization of the variables used in the C2*** linear regression over the power consumed by the compressor Sum1=0 Sum2=0 Sum3=0 Sum4=0 C1*** For each working point we calculate: DO 10 I=1,N C1*** Calculate the evaporating, the condensing pressures C1*** and the value of pratio pevap=1000*EXP(Acl+Bcl/Tev(I)) pcond=1000*EXP(Acl+Bcl/Tcd(I)) pratio=pcond/pevap C1*** Calculate the temperature at the evaporator and C1*** condenser exhaust T1=Tev(I) T3=Tcd(I) C1*** Calculate the enthalpy at the evaporator and C1*** condenser exhaust dhfg=hfgb*((Tc-T1)/(Tc-Tb))**b h1=hfo+cpliq*(T1-To)+dhfg h3=hfo+cpliq*(T3-To) C1*** Calculate the refrigerant mass flow rate MfrRef=Pev(I)/(h1-h3) C1*** Beginning of the loop C1*** First guess of the temperature after the heating-up T1p=T1 Iter=0 20 Iter=Iter+1 Win(I)=MfrRef*Zeta*r*T1p*(pratio**Gm1G-1)/Gm1G T1pp=T1p C1*** Recalculate the temperature after the heating-up T1p=T1+(Pcomp(I)-Win(I))/(MfrRef*cpvap) ErrRel=ABS((T1p-T1pp)/T1pp) C2*** If converged ,leave the second loop IF ((ErrRel.GT.TolRel).AND.(Iter.LE.IterMax)) GOTO 20 IF (Iter.GT.Itermax) THEN ErrDetec=1 GOTO 40 ENDIF C1*** Compare the value of V1p associated with this working point C1*** with the minimum and the maximum value of V1p found so far v1p=Zeta*r*T1p/pevap Vfr1p=MfrRef*v1p IF (Vfr1p.LT.Vmin) THEN Vmin=Vfr1p pratMax=pratio T1pMax=T1p ENDIF IF (Vfr1p.GT.Vmax) THEN Vmax=Vfr1p ENDIF C2*** Up-to-date the variables used in the linear regression over C2*** the power consumed by the compressor Sum1=Sum1+Pcomp(I) Sum2=Sum2+Win(I) Sum3=Sum3+Win(I)**2 Sum4=Sum4+Pcomp(I)*Win(I) 10 CONTINUE C1*** Calculate the parameters Losses and Alpha Alpha=(Sum1*Sum2-N*Sum4)/(Sum2**2-N*Sum3)-1 Losses=(Sum1-(1+Alpha)*Sum2)/N C1*** Calculate the dimensionless standard deviation of the linear C1*** regression over the power consumed by the compressor DO 30 K=1,N Z(K)=Pcomp(K) 30 CONTINUE Slope=1+Alpha dN=SNGL(N) CALL ERROR(dN,Win,Z,Slope,SEw,*31) CALL LINKCK('TYPE74','ERROR',1,99) 31 CONTINUE 40 CONTINUE RETURN 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 PROPERTY (Ifluid,To,cpliq,hfo,cpvap,cpvapcd,hfgb,Tb,Tc, & 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 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