SUBROUTINE TYPE88(TIME,XIN,OUT,T,DTDT,PAR,INFO,ICNTRL,*) C*********************************************************************** C* PROGRAM: Merger of PISSIM1, PISSIM2 & PSIMPL C* C* LANGUAGE: FORTRAN 77 C* C* PURPOSE: Numerical simulation of a chiller in C* steady-state regime. The routine mainly C* calculates the refrigerating capacity and C* the compressor consumption for specified C* working conditions. C* The pressure drop at the compressor exhaust C* can be taken in account. C*********************************************************************** C* INPUT VARIABLES C* XIN(1)=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(2)=Mfrwev Water mass flow rate in the evaporator (Kg/hr) C* XIN(3)=Mfrwcd Water mass flow rate in the condenser (Kg/hr) C* XIN(4)=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(5)=Twev1 This value is equal to the evaporator supply or (C) C* exhaust water temperature according to the value C* of Choice C* XIN(6)=Twcd1 This value is equal to the condenser supply or (C) C* exhaust water temperature according to the value C* of Choice C* XIN(7)=PevG Guess of the refrigerating capacity (KJ/hr) C* XIN(8)=PcdG Guess of the heat rejected in the condenser (KJ/hr) C* XIN(9)=NcFL Number of cylinders used in full load. (-) C* If the pressure drop at the compressor exhaust C* is not taken into account and if there is no C* unloaded cylinder, NcFL MUST BE A C* NEGATIVE NUMBER (not 0) C* XIN(10)=Nc Number of loaded cylinders. (-) C* If there is no unloaded cylinders, Nc C* should be equal to NcFL. Therefore, C* if the pressure drop at the compressor exhaust C* is not taken into account and if there is no C* unloaded cylinder, Nc MUST BE THE SAME C* NEGATIVE NUMBER THAN THE ONE USED FOR NcFL C* C* OUTPUT VARIABLES C* OUT(1)=MfrRef Refrigerant mass flow rate (Kg/hr) C* OUT(2)=Pev Refrigerating capacity (KJ/hr) C* OUT(3)=Pcomp Power consumed by the compressor (KJ/hr) C* OUT(4)=Pcd Heat rejected in the condenser (KJ/hr) C* OUT(5)=COP Coefficient of performance (-) C* OUT(6)=Twev2 This value is equal to the evaporator exhaust (C) C* or supply water temperature according to the value C* of Choice C* OUT(7)=Twcd2 This value is equal to the condenser exhaust (C) C* or supply water temperature according to the value C* of Choice C* C* PARAMETERS C* PAR(1)=AUev Evaporator heat transfer coefficient (KJ/hr-C) C* PAR(2)=AUcd Condenser heat transfer coefficient (KJ/hr-C) C* PAR(3)=Losses Constant part of the electromechanical losses (KJ/hr) C* PAR(4)=Alpha Loss factor allowing to define another (-) C* electromechanical loss which is supposed to be C* proportional to the isentropic/internal(if C* Nc and NcFL are different) power ???? C* PAR(5)=Cf Clearance factor of the compressor (-) C* PAR(6)=VsFL Geometric displacement of the compressor (m**3/s) C* in full-load regime C* PAR(7)=Aex Equivalent nozzle throat area of a cylinder. (m**2) C* Aex is optional but needs to be provided if C* the pressure drop at the compressor exhaust C* is taken in account. C* or C* PAR(7)=Wpumping Internal power of the compressor when all C* the cylinders are unloaded. C* Wpumping only needs to be provided C* when Nc and NcFL are different. C* C* WATER PROPERTY C* CpWat Specific heat of liquid water (KJ/Kg-C) C* C* REFRIGERANT PROPERTIES C* To 0 C -> 273.15=To (K) C* cpliq Mean specific heat in saturated liquid state (KJ/Kg-C) C* hfo Enthalpy of the saturated liquid at the (KJ/Kg) C* reference temperature C* cpvap Mean specific heat at constant pressure (KJ/Kg-C) C* in superheated vapor state C* hfgo Enthalpy of vaporization at the reference (KJ/Kg) C* temperature C* r Gas constant (KJ/Kg-K) C* Zeta Mean compressibility factor (-) 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 When Nc and NcFL are different, the chiller C is supposed to work with only one compressor. C C DEVELOPER: Jean Lebrun C (of the TOOLKIT subroutine) Jean-Pascal Bourdouxhe C Marc Grodent C University of Liege, Belgium C C DATE: October 7, 1993 C Modified for TRNSYS: July 1994-Madison C Mark Nott C C SUBROUTINE CALLED: PROPERTY C*********************************************************************** C INTERNAL VARIABLES: C Twsuev Evaporator supply water temperature (C) C Twexev Evaporator exhaust water temperature (C) C Twsucd Condenser supply water temperature (C) C Twexcd Condenser exhaust water temperature (C) C T1 Evaporating temperature (C) C T1p Temperature after the heating-up (C) C v1p Specific volume after the heating-up (m**3/kg) C p1 Evaporating pressure (Pa) C h1 Enthalpy at the evaporator exhaust (KJ/kg) C T3 Condensing temperature (C) C h3 Enthalpy at the condenser exhaust (KJ/Kg) C p2 Condensing pressure if NcFL is negative (Pa) C or if Nc and NcFL are different C Pressure at point 2 otherwise C p3 Condensing pressure (if NcFL positive (Pa) C and NcFL=Nc) C v2p Specific volume at point 2' (if NcFl positive (m**3/kg) C and NcFL=Nc) C Effvol Volumetric effectiveness of the compressor (-) C Wis Isentropic compression power (KJ/hr) 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 (if NcFl positive and Nc=NcFL) C Vs Geometric displacement of the compressor C for specified working conditions C Vs is only used when Nc and NcFL are different C C Pevp,Pcdp,T1pp and dpexp are variables used in the iterative scheme C*********************************************************************** REAL*8 Mfrwev,Mfrwcd,MfrRef,Losses,Nc,NcFL INTEGER Ifluid INTEGER*4 INFO DIMENSION XIN(10),OUT(7),PAR(7),INFO(15) DOUBLE PRECISION XIN,OUT COMMON /SIM/ TIMEO,TFINAL,DELT,IWARN COMMON /LUNITS/ LUR,LUW,IFORM,LUK COMMON /STORE/ NSTORE,IAV,S(5000) COMMON /CONFIG/ TRNEDT,PERCOM,HEADER,PRLAB,LNKCHK,PRUNIT,IOCHECK . PRWARN DATA TolRel,IterMax/1E-05,100/ C INPUTS Ifluid = INT(XIN(1)+0.1) Mfrwev = XIN(2) Mfrwcd = XIN(3) Choice = XIN(4) Twev1 = XIN(5) Twcd1 = XIN(6) PevG = XIN(7) PcdG = XIN(8) NcFL = XIN(9) Nc = XIN(10) C PARAMETERS AUev = PAR(1) AUcd = PAR(2) Losses = PAR(3) Alpha = PAR(4) Cf = PAR(5) Vsfl = PAR(6) IF ((Nc.EQ.NcFL).AND.(NcFL.GT.0)) THEN Aex = PAR(7) Wpumping= 0 ELSEIF (Nc.NE.NcFL) THEN Wpumping= PAR(7) ELSE Wpumping=0 ENDIF INFO(6)=7 C Selection of the refrigerant CALL PROPERTY(Ifluid,To,cpliq,hfo,cpvap,hfgo,r,Zeta, & Gamma,Acl,Bcl,*5) 5 CONTINUE CpWat=4.187 Gm1G=(Gamma-1)/Gamma Pev=PevG Pcd=PcdG NChoice=INT(Choice) GOTO (10,20,30,40),NChoice C The supply water temperature is known for both evaporator C and condenser 10 Twsuev=Twev1 Twsucd=Twcd1 GOTO 50 C The water temperature is known at the evaporator supply and C at the condenser exhaust 20 Twsuev=Twev1 Twexcd=Twcd1 GOTO 50 C The exhaust water temperature is known for both C evaporator and condenser 30 Twexev=Twev1 Twexcd=Twcd1 GOTO 50 C The water temperature is known at the evaporator exhaust and C at the condenser supply 40 Twexev=Twev1 Twsucd=Twcd1 50 CONTINUE C Calculate the evaporator and condenser effectivenesses Effev=1-EXP(-AUev/(CpWat*Mfrwev)) Effcd=1-EXP(-AUcd/(CpWat*Mfrwcd)) C Calculate the geometric displacement of the compressor C (if Nc and NcFL are different) IF (Nc.NE.NcFL) Vs=Nc/NcFL*VsFL C Beginning of the first loop Iter1=0 60 Iter1=Iter1+1 C Calculate the evaporating temperature according to the C 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 C Calculate the evaporating pressure p1=1000*EXP(Acl+Bcl/(T1+To)) C Calculate the enthalpy at the evaporator exhaust h1=hfo+hfgo+cpvap*T1 C Beginning of the second loop Iter2=0 70 Iter2=Iter2+1 C Calculate the condensing temperature according to the information C 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 C Calculate the condensing pressure IF ((NcFL.GT.0).AND.(Nc.EQ.NcFL)) THEN p3=1000*EXP(Acl+Bcl/(T3+To)) ELSE p2=1000*EXP(Acl+Bcl/(T3+To)) ENDIF C Calculate the enthalpy at the condenser exhaust h3=hfo+cpliq*T3 C Calculate the volumetric efficiency of the compressor C either if NcFL is negative or if Nc and NcFL are different IF ((NcFL.LE.0).OR.(Nc.NE.NcFL)) THEN Effvol=1+Cf-Cf*(p2/p1)**(1/Gamma) ENDIF C Beginning of the third loop C First guess of the temperature after the heating-up T1p=T1 C Calculate the specific volume after the heating-up C and at point 2'if NcFl is greater than 0 and equal to Nc 80 v1p=Zeta*r*1000*(T1p+To)/p1 IF ((NcFL.GT.0).AND.(Nc.EQ.NcFL)) THEN v2p=Zeta*r*1000*(T1p+To)/p3*(p3/p1)**Gm1G C First guess of the pressure drop dpex=1 90 p2=p3+dpex C Calculate the volumetric efficiency of the compressor Effvol=1+Cf-Cf*(p2/p1)**(1/Gamma) C Calculate the refrigerant mass flow rate MfrRef=Effvol*VsFL/v1p*3600 dpexp=dpex C Recalculate the pressure drop dpex=(MfrRef/3600)**2*(v2p/1000)*1/(2*(NcFL*Aex)**2) ErrRel=ABS((dpex-dpexp)/dpexp) C If converged, leave loop IF (ErrRel.GT.TolRel) THEN GOTO 90 ENDIF ELSEIF (NcFL.LE.0) THEN C Calculate the refrigerant mass flow rate C (if NcFL is negative) MfrRef=Effvol*VsFL/v1p*3600 ELSE C Calculate the refrigerant mass flow rate C (if Nc and NcFL are different) MfrRef=Effvol*Vs/v1p*3600 ENDIF C Calculate the isentropic compression power Wis=MfrRef*Zeta*r*(T1p+To)*((p2/p1)**Gm1G-1)/Gm1G T1pp=T1p C Recalculate the temperature after the heating-up T1p=T1+(Losses+Alpha*Wis+(1+Alpha)*(1-Nc/NcFL)* & Wpumping)/(MfrRef*cpvap) ErrRel=ABS((T1p-T1pp)/T1pp) C If converged, leave the third loop IF (ErrRel.GT.TolRel) GOTO 80 C Calculate the power consumed by the compressor Pcomp=Losses+(1+Alpha)*(Wis+(1-Nc/NcFL)*Wpumping) Pcdp=Pcd C Calculate the heat rejected in the condenser Pcd=Pev+Pcomp ErrRel=ABS((Pcd-Pcdp)/Pcdp) C If converged, leave the second loop IF ((ErrRel.GT.TolRel).AND.(Iter2.LE.IterMax)) GO TO 70 C Calculate the refrigerating capacity Pcd=Pcdp Pevp=Pev Pev=MfrRef*(h1-h3) ErrRel=ABS((Pev-Pevp)/Pevp) C If converged, leave the first loop IF ((ErrRel.GT.TolRel).AND.(Iter1.LE.IterMax)) GO TO 60 Pev=Pevp C Calculate the coefficient of performance COP=Pev/Pcomp GOTO (100,110,120,130),NChoice C Calculate the evaporator and condenser exhaust water temperatures 100 Twev2=Twsuev-Pev/(CpWat*Mfrwev) Twcd2=Twsucd+Pcd/(CpWat*Mfrwcd) GOTO 140 C Calculate the evaporator exhaust water temperature and the condenser C supply water temperature 110 Twev2=Twsuev-Pev/(CpWat*Mfrwev) Twcd2=Twexcd-Pcd/(CpWat*Mfrwcd) GOTO 140 C Calculate the evaporator and condenser supply water temperatures 120 Twev2=Twexev+Pev/(CpWat*Mfrwev) Twcd2=Twexcd-Pcd/(CpWat*Mfrwcd) GOTO 140 C Calculate the evaporator supply water temperature and the condenser C exhaust water temperature 130 Twev2=Twexev+Pev/(CpWat*Mfrwev) Twcd2=Twsucd+Pcd/(CpWat*Mfrwcd) 140 CONTINUE C OUTPUT OUT(1) = MfrRef OUT(2) = Pev OUT(3) = Pcomp OUT(4) = Pcd OUT(5) = COP OUT(6) = Twev2 OUT(7) = Twcd2 RETURN 1 END