SUBROUTINE TYPE63 (TIME,XIN,OUT,T,DTDT,PAR,INFO,ICNTRL,*) C************************************************************************ C* Copyright ASHRAE A Toolkit for Primary HVAC System Energy C* Calculation C************************************************************************ C* SUBROUTINE: TYPE63 (TURBIFL) C* C* LANGUAGE: FORTRAN 77 C* C* PURPOSE: Calculates the shaft power when the gas C* turbine is running in steady-state regime C*********************************************************************** C* INPUT VARIABLES: C* Ifuel Selection of the fuel (-) C* If Ifuel C* =1: Light fuel oil (liquid fuel) C* =2: Heavy fuel oil (liquid fuel) C* =3: Domestic gas oil (liquid fuel) C* =4: Methane (gaseous fuel) C* xin(1) (-) C* Choice If Choice (-) C* = 1: the gas generator has to satisfy the C* strongest constraint between its maximum flue C* gas temperature at the exhaust of the C* adiabatic combustion chamber and its maximum C* shaft speed C* = 2: the gas generator is running at its maximum C* flue gas temperature at the exhaust of the C* adiabatic combustion chamber (no constraint C* on the value of the shaft speed) C* = 3: the gas generator is running at its maximum C* shaft speed (no constraint on the value of C* the flue gas temperature at the exhaust of the C* adiabatic combustion chamber) C* xin(2) (-) C* Ta Ambient air temperature (K) C* xin(3) (øC) C* pa Ambient air pressure (Pa) C* xin(4) (atm) C* T3max Maximum flue gas temperature at the exhaust (K) C* of the adiabatic combustion chamber C* xin(5) (øC) C* N1max Maximum gas generator shaft speed (1/s) C* xin(6) (rpm) C* N2 Shaft speed of the second expander (1/s) C* xin(7) (rpm) C* C* OUTPUT VARIABLES C* T3 Flue gas temperature at the exhaust of the (K) C* adiabatic combustion chamber C* out(1) (øC) C* N1 Gas generator shaft speed (1/s) C* out(2) (rpm) C* Wsh Available shaft power (W) C* out(3) (kJ/hr) C* Effic Gas turbine efficiency (-) C* out(4) (-) C* sfc Specific fuel consumption (kg/kW/h) C* out(5) (kg/kW/h) C* MfrFuel Fuel mass flow rate (kg/s) C* out(6) (kg/hr) C* MfrGas Flue gas mass flow rate (kg/s) C* out(7) (kg/hr) C* pratioComp Compressor pressure ratio (-) C* out(8) (-) C* T5 Flue gas temperature at the second expander (K) C* exhaust C* out(9) (øC) C* cpMeanGas Mean specific heat for the second expander (J/kg/K) C* out(10) (J/kg/hr) C* ErrDetec This variable is equal to 1 if the routine (-) C* does not converge.Otherwise, ErrDetec is C* equal to 0 . C* out(11) (-) C* C* PARAMETERS C* EffCE1 Compressor and first expander polytropic (-) C* effectiveness C* par(1) (-) C* Alaval Nozzle throat area (m**2) C* par(2) (m**2) C* Dlaval Second expander diameter (m) C* par(3) (m) C* Alpha1 Angle characterizing the orientation of the (-) C* Laval blade in the flow C* par(4) (-) C* RedMfr Reduced mass flow rate in the first (m*s*SQRT(K)) C* expander C* par(5) (m*s*SQRT(K)) C* MTPN1 Parameter characterizing the mass flow (m*s**2*K) C* rate in the compressor C* par(6) (m*s**2*K) C* C* FUEL PROPERTIES C* Cweight Weight of carbon in 1 kg of fuel (kg) C* FLHV The fuel lower heating value (J/kg) C* Tr Reference temperature at which the FLHV is (K) C* evaluated C* Cfuel Fuel specific heat (J/kg/K) C* C* AIR PROPERTIES C* RAir Air constant (J/kg/K) C*********************************************************************** C MAJOR RESTRICTIONS: Same constant polytropic effectiveness for C the compressor and the first expander. C Chocked regime in the compressor with C reduced flow rate proportional to reduced C shaft power. C Chocked regime in the first expander. C Negligeable pressure drops, combustion C and mechanical losses. C Second expander assimilated to an ideal C Laval stage with isentropic nozzle and C without any diffuser. 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: FAR C ENTHALP C FUEL C LINKCK C HAIR C TEMPFLUE C*********************************************************************** C INTERNAL VARIABLES: C Kmolp Array containing the (kmol constituent/kg fuel) C composition of the combustion products C Sum Kmol of flue gas per kg of fuel (kmol/kg) C cpMeanAir Mean specific heat of air (J/kg/K) C MfrAir Air mass flow rate (kg/s) C RGas Flue gas constant (J/kg flue gas/K) C p2 Pressure at the compressor exhaust (Pa) C p3 Pressure at the combustion chamber exhaust (Pa) C T2 Air temperature at the compressor exhaust (K) C h2 Air enthalpy at the compressor exhaust (J/kg) C h1 Air enthalpy at the compressor supply (J/kg) C GamComp Mean isentropic exponent for the compressor (-) C GamEx1 Mean isentropic exponent for the first expander (-) C GamEx2 Mean isentropic exponent for the second expander (-) C Tfuel Fuel temperature (K) C Fratio Fuel/air ratio (-) C h3s Flue gas enthalpy at the exhaust of (J/kg fuel) C the adiabatic combustion chamber C h3 Flue gas enthalpy at the exhaust of (J/kg flue gas) C the adiabatic combustion chamber C Wcomp Power consumed by the compressor (W) C Wexp1 Power given by the first expander (W) C p4 Pressure at the first expander exhaust (Pa) C h4 Gas enthalpy at the exhaust of the (J/kg flue gas) C first expander C T4 Flue gas temperature at the exhaust of the (K) C first expander C h5 Gas enthalpy at the exhaust of the (J/kg flue gas) C second expander C pCrit Critical pressure at the nozzle throat (Pa) C pThr Pressure at the nozzle throat (Pa) C U Second expander speed (m/s) C C41 Flow speed at the Laval blade supply (m/s) C MW Array containing the molecular weights (kg/kmol) C of the species (H2,O2,N2,CO2,H2O) C MGas Flue gas molecular weight (kg/kmol) C TolRel Relative error tolerance (-) C ErrRel Relative error (-) C C GamCompp,GamEx2p,p2p,N1p,T3p,Iter1,Iter2 and Itermax are variables C used in the iterative schema. C*********************************************************************** INTEGER*4 INFO DOUBLE PRECISION XIN,OUT REAL Kmolp(5),MW(5) REAL Ifuel,MfrAir,MfrFuel,MfrGas,N1,N1max,N1p,N2,MGas,MTPN1 DIMENSION PAR(6),XIN(7),OUT(11),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 COMMON/COMCP/PFCP(5,10) INFO(6)=11 DATA TolRel,IterMax,Pi/1E-04,300,3.14159265359/ DATA RAir/287.06/ DATA MW/2,32,28,44,18/ C*** INPUTS 7 (converted in SI units) C************ Ifuel=SNGL(xin(1)) Choice=SNGL(xin(2)) Ta=SNGL(xin(3)+273.15) pa=SNGL(xin(4)*101325) T3max=SNGL(xin(5)+273.15) N1max=SNGL(xin(6)/60) N2=SNGL(xin(7)/60) C*** PARAMETERS 6 (converted in SI units) C**************** EffCE1=par(1) Alaval=par(2) Dlaval=par(3) Alpha1=par(4) RedMfr=par(5) MTPN1=par(6) ErrDetec=0 CALL FUEL (Ifuel,Cweight,FLHV,Tr,Cfuel,*1) CALL LINKCK('TYPE63','FUEL',1,99) 1 CONTINUE C1*** Calculate the second expander speed U=Pi*Dlaval*N2 C1*** Determine the constraint to verify according to the value C1*** of Choice J=1 IF (Choice.EQ.3) THEN J=2 ENDIF 10 IF (J.EQ.1) THEN C2*** The gas generator is running at its maximum C2*** flue gas temperature at the exhaust of the C2*** adiabatic combustion chamber (no constraint C2*** on the value of the shaft speed) T3=T3max C2*** Guess of the value of the gas generator shaft speed N1=0.9*N1max ELSE C2*** The gas generator is running at its maximum C2*** shaft speed (no constraint on the value of C2*** the flue gas temperature at the exhaust of the C2*** adiabatic combustion chamber) N1=N1max C2*** Guess of the value of the flue gas temperature at the C2*** exhaust of the combustion chamber T3=0.9*T3max ENDIF C2*** Guess of the pressure at the compressor exhaust p2=8*pa C1*** Begin of the first loop Iter1=0 20 Iter1=Iter1+1 C1*** Calculate the air mass flow rate MfrAir=MTPN1*pa*N1/Ta C1*** Begin of the second loop Iter2=0 30 Iter2=Iter2+1 C2*** Calculate the mean isentropic coefficient for the compressor C2*** as well as the air temperature and enthalpy after compression Iter=0 GamComp=1.4 31 Iter=Iter+1 Gm1G=(GamComp-1)/GamComp T2=Ta*(p2/pa)**(Gm1G/EffCE1) CALL HAIR(T2,h2,*33) CALL LINKCK('TYPE63','HAIR',1,99) 33 CONTINUE CALL HAIR(Ta,h1,*34) CALL LINKCK('TYPE63','HAIR',1,99) 34 CONTINUE cpMeanAir=(h2-h1)/(T2-Ta) GamCompp=GamComp GamComp=cpMeanAir/(cpMeanAir-RAir) ErrRel=ABS((Gamcomp-Gamcompp)/Gamcompp) C2*** If converged, then leave loop IF ((ErrRel.GT.TolRel).AND.(Iter.LE.IterMax)) GOTO 31 C2*** Calculate the compression power Wcomp=MfrAir*(h2-h1) C1*** Calculate the fuel/air ratio Tfuel=333 CALL FAR(Ifuel,T2,Tfuel,T3,Fratio,Kmolp,h3s,*35) CALL LINKCK('TYPE63','FAR',1,99) 35 CONTINUE h3=h3s/(1+1/Fratio) C2*** Calculate the flue gas constant Sum=Kmolp(2)+Kmolp(3)+Kmolp(4)+Kmolp(5) MGas=(MW(2)*Kmolp(2)+MW(3)*Kmolp(3)+MW(4)*Kmolp(4)+MW(5)* & Kmolp(5))/Sum RGas=8314/MGas C1*** Calculate the fuel mass flow rate MfrFuel=MfrAir*Fratio C1*** Calculate the flue gas mass flow rate MfrGas=MfrAir+MfrFuel C1*** Calculate a new value of the pressure at the compressor exhaust p2p=p2 p2=MfrGas*SQRT(T3)/RedMfr ErrRel=ABS((p2-p2p)/p2p) C2*** If converged, then leave the second loop IF ((ErrRel.GT.TolRel).AND.(Iter2.LE.IterMax)) GOTO 30 IF (Iter2.GT.IterMax) THEN ErrDetec=1 GOTO 40 ELSE p2=p2p ENDIF C1*** Calculate the compressor pressure ratio pratioComp=p2/pa C1*** Calculate the pressure in the adiabatic combustion chamber p3=p2 C1*** Calculate the temperature at the first expander exhaust C1*** (Newton-Raphson method) Wexp1=Wcomp h4=h3-Wexp1/MfrGas CALL TEMPFLUE(Kmolp,Fratio,T3,h4,T4,*37) CALL LINKCK('TYPE63','TEMPFLUE',1,99) 37 CONTINUE C2*** Calculate the mean isentropic coefficient for the first C2*** expander cpMean=(h3-h4)/(T3-T4) GamEx1=cpMean/(cpMean-RGas) C1*** Calculate the pressure at the first expander exhaust Gm1G=(GamEx1-1)/GamEx1 p4=p3*(T4/T3)**(1/(Gm1G*EffCE1)) C2*** Calculate the mean specific heat and the mean isentropic C2*** coefficient for the second expander Iter=0 cpMeanGas=cpMean GamEx2=cpMeanGas/(cpMeanGas-RGas) 32 Iter=Iter+1 Gm1G=(GamEx2-1)/GamEx2 C41=SQRT(2*cpMeanGas*T4*(1-(pa/p4)**Gm1G)) Wsh=MfrGas*2*U*(C41*COS(Alpha1)-U) h5=h4-Wsh/MfrGas CALL TEMPFLUE(Kmolp,Fratio,T4,h5,T5,*38) CALL LINKCK('TYPE63','TEMPFLUE',1,99) 38 CONTINUE cpMeanGasp=cpMeanGas cpMeanGas=(h4-h5)/(T4-T5) GamEx2p=GamEx2 GamEx2=cpMeanGas/(cpMeanGas-RGas) ErrRel=ABS((GamEx2-GamEx2p)/GamEx2p) C2*** If converged, then leave loop IF ((ErrRel.GT.TolRel).AND.(Iter.LE.IterMax)) GOTO 32 GamEx2=GamEx2p cpMeanGas=cpMeanGasp C2*** Calculate the critical pressure pCrit=p4*(2/(GamEx2+1))**(1/Gm1G) C2*** Calculate the pressure at the nozzle throat pThr=MAX(pCrit,pa) C1*** Calculate a new value of the flue gas temperature at the C1*** adiabatic combustion chamber exhaust or of the gas C1*** generator speed according to the value of Choice prate=pThr/p4 Var=Alaval/SQRT(RGas)*SQRT(2/Gm1G)*prate**(1/GamEx2) & *SQRT(1-prate**Gm1G) MfrGas=Var*p4/SQRT(T4) IF (J.EQ.1) THEN MfrAir=MfrGas-MfrFuel N1p=N1 N1=MfrAir*Ta/(pa*MTPN1) ErrRel=ABS((N1-N1p)/N1p) ELSE T3p=T3 h3=h4+Wexp1/MfrGas CALL TEMPFLUE(Kmolp,Fratio,T3p,h3,T3,*39) CALL LINKCK('TYPE63','TEMPFLUE',1,99) 39 CONTINUE ErrRel=ABS((T3-T3p)/T3p) ENDIF C2*** If converged, then leave the first loop IF ((ErrRel.GT.TolRel).AND.(Iter1.LE.IterMax)) GOTO 20 IF (Iter1.GT.IterMax) THEN ErrDetec=1 GOTO 40 ELSE IF (J.EQ.1) THEN N1=N1p ELSE T3=T3p ENDIF ENDIF C1*** Calculate the gas turbine efficiency and the specific fuel C1*** consumption Effic=Wsh/(MfrFuel*FLHV) sfc=3600E3*MfrFuel/Wsh C1*** If Choice is equal to 1, each constraint has to be satisfied IF (Choice.EQ.1) THEN IF ((J.EQ.1).AND.(N1.GT.N1max)) THEN J=2 GOTO 10 ENDIF IF ((J.EQ.2).AND.(T3.GT.T3max)) THEN ErrDetec=1 GOTO 40 ENDIF ENDIF 40 CONTINUE C*** OUTPUTS 11 (converted in TRNSYS units) C************** out(1)=DBLE(T3-273.15) out(2)=DBLE(N1*60) out(3)=DBLE(Wsh*3.6) out(4)=DBLE(Effic) out(5)=DBLE(sfc) out(6)=DBLE(MfrFuel*3600.) out(7)=DBLE(MfrGas*3600.) out(8)=DBLE(pratioComp) out(9)=DBLE(T5-273.15) out(10)=DBLE(cpMeanGas) out(11)=DBLE(ErrDetec) RETURN 1 END SUBROUTINE FAR(Ifuel,Tair,Tfuel,Tadiab,Fratio,Kmolp,hprod,*) C************************************************************************ C* Copyright ASHRAE A Toolkit for Primary HVAC System Energy C* Calculation C*********************************************************************** C* SUBROUTINE: FAR C* C* LANGUAGE: FORTRAN 77 C* C* PURPOSE: Calculation of the fuel/air ratio, the C* composition and the enthalpy of the flue C* gas in an adiabatic combustion chamber C*********************************************************************** C* INPUT VARIABLES C* Ifuel Selection of the fuel (-) C* If Ifuel C* =1: Light fuel oil (liquid fuel) C* =2: Heavy fuel oil (liquid fuel) C* =3: Domestic gas oil (liquid fuel) C* =4: Methane (gaseous fuel) C* Tair Air temperature (K) C* Tfuel Fuel temperature (K) C* Tadiab Adiabatic temperature of the combustion products (K) C* C* OUTPUT VARIABLES C* Fratio Actual fuel/air ratio (-) C* Kmolp Array containing the (kmol constituent/kg fuel) C* composition of the combustion products C* This output is not given if the routine runs alone C* hprod Enthalpy of the combustion products (J/kg fuel) C* C* FUEL PROPERTIES C* Cweight Weight of carbon in 1kg of fuel (kg) C* FLHV Fuel lower heating value (J/kg) C* Tr Reference temperature at which the FLHV is C* evaluated (K) C* Cfuel Fuel specific heat (J/kg/K) C*********************************************************************** C MAJOR RESTRICTIONS: The combustion process is assumed to be C adiabatic. The combustion reaction is C assumed to be complete and dissociation is C not taken into account. C C DEVELOPER: Philippe Ngendakumana C Marc Grodent C Jean-Pascal Bourdouxhe C University of LiŠge, Belgium C C DATE: March 1, 1995 C C SUBROUTINE CALLED: ENTHALP C*********************************************************************** C INTERNAL VARIABLES C Excess Excess of air (-) C Hweight Weight of hydrogen in 1 kg of fuel (kg) C ParC Number of carbon atoms in equivalent hydrocarbon (-) C fuel C ParH Number of hydrogen atoms in equivalent (-) C hydrocarbon fuel C O2st Number of kmol of oxygen reacting with 1 kmol (kmol) C of fuel in a stoichiometric combustion process C Fratiost Stoichiometric fuel/air ratio (-) C Qas Fraction of the heat transfer in the (W) C air preheater C Qf Heat transfer in the fuel preheater (W) C Qr Heat transfer in the isothermal reactor (W) C MW Array containing the molecular weights (kg/kmol) C of the species (H2,O2,N2,CO2,H2O) C ZC Atomic number of carbon (-) C ParAir Number of kmol of nitrogen per kmol of oxygen (-) C (air composition) C hO2 Enthalpy of oxygen per kmol of oxygen (J/kmol) C hrO2 Enthalpy of oxygen per kmol of oxygen (J/kmol) C at the temperature Tr C hN2 Enthalpy of nitrogen per kmol of nitrogen (J/kmol) C hrN2 Enthalpy of nitrogen per kmol of nitrogen (J/kmol) C at the temperature Tr C hCO2 Enthalpy of carbon dioxide per kmol of (J/kmol) C carbon dioxide C hrCO2 Enthalpy of carbon dioxide per kmol of (J/kmol) C carbon dioxide at the temperature Tr C hH2O Enthalpy of water per kmol of water (J/kmol) C hrH2O Enthalpy of water per kmol of water (J/kmol) C at the temperature Tr C*********************************************************************** REAL Kmolp(5),MW(5) DATA MW/2,32,28,44,18/ DATA ZC/12/ CALL FUEL(Ifuel,Cweight,FLHV,Tr,Cfuel,*2) CALL LINKCK('FAR','FUEL',1,99) 2 CONTINUE C2*** The only fuel constituents are carbon and hydrogen Hweight=1-Cweight C2*** Calculate the number of carbon atoms in C2*** equivalent hydrocarbon fuel ParC=Cweight/ZC C2*** Calculate the number of hydrogen atoms in C2*** equivalent hydrocarbon fuel ParH=Hweight C2*** Calculate the number of kmol of nitrogen per kmol of oxygen C2*** (air composition) ParAir=0.79/0.21 C2*** Calculate the number of kmol of oxygen reacting with 1 kmol C2*** of fuel in a stoichiometric combustion process O2st=ParC+ParH/4 C2*** Calculate the stoichiometric fuel/air ratio Fratiost=1/(O2st*(MW(2)+ParAir*MW(3))) C1*** Calculate the amount of carbon dioxide and water in C1*** the combustion products Kmolp(4)=ParC Kmolp(5)=ParH/2 C1*** Calculate the heat transfer in the fuel preheater Qf=Cfuel*(Tr-Tfuel) C1*** Calculate a fraction of the heat transfer in the air preheater CALL ENTHALP(Tair,2,hO2,*3) CALL LINKCK('FAR','ENTHALP',1,99) 3 CONTINUE CALL ENTHALP(Tair,3,hN2,*4) CALL LINKCK('FAR','ENTHALP',1,99) 4 CONTINUE CALL ENTHALP(Tr,2,hrO2,*5) CALL LINKCK('FAR','ENTHALP',1,99) 5 CONTINUE CALL ENTHALP(Tr,3,hrN2,*6) CALL LINKCK('FAR','ENTHALP',1,99) 6 CONTINUE Qas=(hrO2-hO2)+ParAir*(hrN2-hN2) C1*** Calculate the heat transfer in the isothermal reactor Qr=FLHV C1*** Calculate the excess of air CALL ENTHALP(Tadiab,2,hO2,*10) CALL LINKCK('FAR','ENTHALP',1,99) 10 CONTINUE DhO2=hO2-hrO2 CALL ENTHALP(Tadiab,3,hN2,*11) CALL LINKCK('FAR','ENTHALP',1,99) 11 CONTINUE DhN2=hN2-hrN2 CALL ENTHALP(Tadiab,4,hCO2,*12) CALL LINKCK('FAR','ENTHALP',1,99) 12 CONTINUE CALL ENTHALP(Tr,4,hrCO2,*13) CALL LINKCK('FAR','ENTHALP',1,99) 13 CONTINUE DhCO2=hCO2-hrCO2 CALL ENTHALP(Tadiab,5,hH2O,*14) CALL LINKCK('FAR','ENTHALP',1,99) 14 CONTINUE CALL ENTHALP(Tr,5,hrH2O,*15) CALL LINKCK('FAR','ENTHALP',1,99) 15 CONTINUE DhH2O=hH2O-hrH2O Excess=(Qr-Qf-Kmolp(4)*DhCO2-Kmolp(5)*DhH2O-O2st*(Qas+ & ParAir*DhN2))/(O2st*(DhO2+ParAir*DhN2+Qas)) C1*** Calculate the actual fuel/air ratio Fratio=Fratiost/(1+Excess) C1*** Calculate the amount of oxygen and nitrogen in C1*** the combustion products Kmolp(2)=O2st*Excess Kmolp(3)=O2st*(1+Excess)*ParAir C1*** Calculate the enthalpy of the combustion products C1*** at the adiabatic temperature hprod=0 DO 20 KP=2,5 CALL ENTHALP(Tadiab,KP,hKP,*16) CALL LINKCK('FAR','ENTHALP',1,99) 16 CONTINUE hprod=hprod+Kmolp(KP)*hKP 20 CONTINUE RETURN 1 END SUBROUTINE ENTHALP (Temp,I,Enthalpy,*) C************************************************************************ C* Copyright ASHRAE A Toolkit for Primary HVAC System Energy C* Calculation C*********************************************************************** C* SUBROUTINE: ENTHALP C* C* LANGUAGE: FORTRAN 77 C* C* PURPOSE: Calculate the enthalpy (J/kmol) of each C* species (H2,O2,N2,CO2,H2O) at a given C* temperature C*********************************************************************** C* INPUT VARIABLES C* Temp Temperature at which enthalpy must be calculated (K) C* I Selection of the species to be considered (-) C* I=1: H2 C* I=2: O2 C* I=3: N2 C* I=4: CO2 C* I=5: H2O C* C* OUTPUT VARIABLES C* Enthalpy Enthalpy of the species (J/kmol) C*********************************************************************** C DEVELOPER: Philippe Ngendakumana C Marc Grodent C University of LiŠge, Belgium C C DATE: March 1, 1995 C C REFERENCE: A. Brohmer and P. Kreuter C FEV Motorentechnik GmbH & Co KG C Aachen, Germany C*********************************************************************** C INTERNAL VARIABLES C PFCP Array containing the coefficients used (J/kmol/K) C in the polynomial expressions C Tref Array containing the temperatures at which (K) C the reference enthalpies are calculated C href Array containing the reference enthalpies (J/kmol) C h Enthalpy of species I (J/kmol) C J Loop counter C*********************************************************************** COMMON/COMCP/PFCP(5,10) COMMON/THREF/Tref(5),href(5) h=href(I) Enthalpy=0 DO 10 J=1,10 h=h+((PFCP(I,J)*Temp**J)-(PFCP(I,J)*Tref(I)**J))/J 10 CONTINUE Enthalpy=h RETURN 1 END BLOCK DATA COMMON/COMCP/PFCP(5,10) COMMON/THREF/Tref(5),href(5) C1*** Coefficients are given for H2 DATA PFCP(1,1),PFCP(1,2),PFCP(1,3), $PFCP(1,4),PFCP(1,5),PFCP(1,6),PFCP(1,7), $PFCP(1,8),PFCP(1,9),PFCP(1,10)/ $ 2.12183E+04, 4.90483E+01,-1.18908E-01, 1.50167E-04, $-1.07285E-07, 4.66644E-11,-1.26418E-14, 2.08562E-18, $-1.91864E-22, 7.54661E-27/ C1*** Coefficients are given for O2 DATA PFCP(2,1),PFCP(2,2),PFCP(2,3), $PFCP(2,4),PFCP(2,5),PFCP(2,6),PFCP(2,7), $PFCP(2,8),PFCP(2,9),PFCP(2,10)/ $ 3.12398E+04,-2.51025E+01, 9.50643E-02,-1.29283E-04, $ 9.56020E-08,-4.25012E-11, 1.16866E-14,-1.94778E-18, $ 1.80410E-22,-7.12717E-27/ C1*** Coefficients are given for N2 DATA PFCP(3,1),PFCP(3,2),PFCP(3,3), $PFCP(3,4),PFCP(3,5),PFCP(3,6),PFCP(3,7), $PFCP(3,8),PFCP(3,9),PFCP(3,10)/ $ 3.10052E+04,-1.65866E+01, 4.37297E-02,-4.10720E-05, $ 2.08732E-08,-6.27548E-12, 1.11654E-15,-1.08777E-19, $ 4.47487E-24, 0.E0 / C1*** Coefficients are given for CO2 DATA PFCP(4,1),PFCP(4,2),PFCP(4,3), $PFCP(4,4),PFCP(4,5),PFCP(4,6),PFCP(4,7), $PFCP(4,8),PFCP(4,9),PFCP(4,10)/ $ 1.89318E+04, 8.20742E+01,-8.47204E-02, 5.92177E-05, $-2.92546E-08, 1.01523E-11,-2.39525E-15, 3.62658E-19, $-3.15882E-23, 1.19863E-27/ C1*** Coefficients are given for H2O DATA PFCP(5,1),PFCP(5,2),PFCP(5,3), $PFCP(5,4),PFCP(5,5),PFCP(5,6),PFCP(5,7), $PFCP(5,8),PFCP(5,9),PFCP(5,10)/ $ 3.42084E+04,-1.04650E+01, 3.61342E-02,-2.73709E-05, $ 1.12406E-08,-2.93883E-12, 5.25323E-16,-6.54907E-20, $ 5.27765E-24,-2.04468E-28/ C1*** Reference values are given for H2 DATA Tref(1),href(1)/2.E3,6.144129E7/ C1*** Reference values are given for O2 DATA Tref(2),Href(2)/2.E3,6.7926643E7/ C1*** Reference values are given for N2 DATA Tref(3),Href(3)/2.E3,6.485353E7/ C1*** Reference values are given for CO2 DATA Tref(4),Href(4)/2.E3,-2.9253172E8/ C1*** Reference values are given for H2O DATA Tref(5),Href(5)/2.E3,-1.5643141E8/ END SUBROUTINE FUEL (Ifuel,Cweight,FLHV,Tr,Cfuel,*) C************************************************************************ C* Copyright ASHRAE A Toolkit for Primary HVAC System Energy C* Calculation C*********************************************************************** C* SUBROUTINE: FUEL C* C* LANGUAGE: FORTRAN 77 C* C* PURPOSE: Selection of the properties of a given fuel C*********************************************************************** C* INPUT VARIABLE: C* Ifuel Selection of the fuel (-) C* If Ifuel C* =1: Light fuel oil (liquid fuel) C* =2: Heavy fuel oil (liquid fuel) C* =3: Domestic gas oil (liquid fuel) C* =4: Methane (gaseous fuel) C*********************************************************************** C* OUTPUT VARIABLES: C* Cweight Weight of carbon in 1kg of fuel (kg) C* FLHV The fuel lower heating value (J/kg) C* Tr Reference temperature at which the FLHV is (K) C* evaluated C* Cfuel Fuel specific heat (J/kg/K) C*********************************************************************** C DEVELOPER: Jean-Pascal Bourdouxhe C Marc Grodent C University of LiŠge, Belgium C C DATE: March 1, 1995 C*********************************************************************** REAL Ifuel IF (Ifuel.EQ.1) THEN Cweight=0.88 FLHV=40910E3 Tr=298.15 Cfuel=1840 ENDIF IF (Ifuel.EQ.2) THEN Cweight=0.89 FLHV=40430E3 Tr=298.15 Cfuel=1840 ENDIF IF (Ifuel.EQ.3) THEN Cweight=0.86 FLHV=42770E3 Tr=298.15 Cfuel=1880 ENDIF IF (Ifuel.EQ.4) THEN Cweight=0.75 FLHV=49997E3 Tr=289.15 Cfuel=2183 ENDIF RETURN 1 END SUBROUTINE HAIR(T,HA,*) C************************************************************************ C* Copyright ASHRAE A Toolkit for Primary HVAC System Energy C* Calculation C*********************************************************************** C* SUBROUTINE: HAIR C* C* LANGUAGE: FORTRAN 77 C* C* PURPOSE: The purpose of this function is to C* calculate the air enthalpy (expressed in C* J/kg air) at a given temperature C*********************************************************************** C* INPUT VARIABLES C* T Air temperature (K) C* OUTPUT VARIABLES C* HA Air enthalpy (J/kg air) C*********************************************************************** C MAJOR RESTRICTIONS: None. C C DEVELOPER: Jean-Pascal Bourdouxhe C Marc Grodent C University of LiŠge, Belgium C C DATE: March 1, 1995 C C SUBROUTINE CALLED: ENTHALP C LINKCK C*********************************************************************** C INTERNAL VARIABLES: C hO2 O2 enthalpy (J/kmol) C hN2 N2 enthalpy (J/kmol) C*********************************************************************** CALL ENTHALP(T,2,hO2,*1) CALL LINKCK('HAIR','ENTHALP',1,99) 1 CONTINUE CALL ENTHALP(T,3,hN2,*2) CALL LINKCK('HAIR','ENTHALP',1,99) 2 CONTINUE HA=(0.21/28.84)*hO2+(0.79/28.84)*hN2 RETURN 1 END SUBROUTINE TEMPFLUE(Kmolp,Fratio,Tguess,h,T,*) C************************************************************************ C* Copyright ASHRAE A Toolkit for Primary HVAC System Energy C* Calculation C*********************************************************************** C* SUBROUTINE: TEMPFLUE C* C* LANGUAGE: FORTRAN 77 C* C* PURPOSE: The purpose of this function is to C* calculate the flue gas temperature C* corresponding to a given flue gas C* enthalpy (expressed in J/kg flue gas) C* (Newton-Raphson method) C*********************************************************************** C* INPUT VARIABLES: C* Kmolp Array containing the (kmol constituent/kg fuel) C* composition of the combustion products C* Fratio Fuel/air ratio (-) C* Tguess Value used for the guess of the flue gas (K) C* temperature C* h Flue gas enthalpy (J/kg flue gas) C* C* OUTPUT VARIABLES: C* T Flue gas temperature (K) C*********************************************************************** C MAJOR RESTRICTIONS: None. C C DEVELOPER: Jean-Pascal Bourdouxhe C Marc Grodent C University of LiŠge, Belgium C C DATE: March 1, 1995 C C SUBROUTINES CALLED: ENTHALP C LINKCK C*********************************************************************** C INTERNAL VARIABLES C TolRel Relative error tolerance (-) C Fct Value of the function to be nullified (J/kg flue gas) C Dfct Value of the first derivative (J/(kg flue gas*K)) C ErrRel Relative error (-) C C hpi, hgcal, hgcal1, Sum1, Sum2, Tp are variables used in the C Newton-Raphson method. C*********************************************************************** REAL Kmolp(5) COMMON/COMCP/PFCP(5,10) DATA TolRel/1E-06/ C2*** First guess of the temperature T=0.9*Tguess 10 hgcal1=0 DO 20 J=2,5 CALL ENTHALP(T,J,hpi,*12) CALL LINKCK('TEMPFLUE','ENTHALP',1,99) 12 CONTINUE hgcal1=hgcal1+Kmolp(J)*hpi 20 CONTINUE hgcal=hgcal1/(1+1/Fratio) C2*** Calculate the function to nullify Fct=hgcal-h C2*** Calculate the value of the first derivative Sum1=0 DO 30 K=1,5 Sum2=0 DO 40 J=1,10 Sum2=Sum2+PFCP(K,J)*T**(J-1) 40 CONTINUE Sum1=Sum1+Kmolp(K)*Sum2 30 CONTINUE DFct=Sum1/(1+1/Fratio) C2*** A new estimated value is calculated Tp=T T=T-Fct/DFct ErrRel=ABS((T-Tp)/Tp) C2*** If converged, then leave the loop IF (ErrRel.GT.TolRel) GOTO 10 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