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