SUBROUTINE TYPE905 (TIME,XIN,OUT,T,DTDT,PAR,INFO,ICNTRL,*) !DEC$ATTRIBUTES DLLEXPORT :: TYPE 905 C************************************************************************ C* Copyright ASHRAE A Toolkit for Primary HVAC System Energy C* Calculation C*********************************************************************** C* SUBROUTINE: TYPE905 (ENGIFLSI) C* C* LANGUAGE: FORTRAN 77 C* C* PURPOSE: Calculates the shaft power when the gas C* engine is running in steady-state regime. C*********************************************************************** C* INPUT VARIABLES: C* Fratio Fuel/air ratio (-) C* xin(1) (-) C* Ta Air temperature (K) C* xin(2) (øC) C* MfrW Water mass flow rate (kg/s) C* xin(3) (kg/hr) C* Twsu Supply water temperature (K) C* xin(4) (øC) C* pa Air pressure (Pa) C* xin(5) (atm) C* N Rotation speed (1/s) C* xin(6) (rpm) C* C* OUTPUT VARIABLES C* CPgas Mean specific heat of the combustion (J/kg/K) C* products C* out(1) (J/kg/øC) C* Twex Exhaust water temperature (K) C* out(2) (øC) C* Tgex Flue gas temperature at the exhaust of the (K) C* gas-water heat exchanger C* out(3) (øC) C* Wsh Shaft power (W) C* out(4) (kJ/hr) C* MfrFuel Gas mass flow rate (kg/s) C* out(5) (kg/hr) C* MfrGas Flue gas mass flow rate (kg/s) C* out(6) (kg/hr) C* Effic Gas engine efficiency (-) C* out(7) (-) C* ErrDetec = 1: the ratio of water capacity flow rate to (-) C* flue gas capacity flow rate is too small ( <1 ) C* = 2: the rotation speed for specified working C* conditions is lower than the minimum rotation C* speed (Nmin) C* = 3: the rotation speed for specified working C* conditions is greater than the maximum C* rotation speed (Nmax) C* = 4: the routine does not converge C* In these cases, the routine stops running; C* otherwise this variable is equal to 0. C* out(8) (-) C* C* PARAMETERS C* i Intermittency factor (-) C* par(1) (-) C* Vs Swept volume corresponding to all the (m**3) C* cylinders C* par(2) (m**3) C* Athroat Nozzle throat area (m**2) C* par(3) (m**2) C* EffiInt Internal efficiency (-) C* par(4) (-) C* Tlo Torque associated with the mechanical losses (N*m) C* and the auxiliary consumptions C* par(5) (N*m) C* AUgwNoEng Gas-water heat transfer coefficient in nominal (W/K) C* conditions C* par(6) (kJ/hr/øC) C* MfrGasNom Flue gas mass flow rate in nominal conditions (kg/s) C* par(7) (kg/hr) C* AUwenvEng Water-environment heat transfer coefficient (W/K) C* par(8) (kJ/hr/øC) C* C* AIR PROPERTIES C* CpAir Air specific heat (J/kg/K) C* RAir Air constant (J/kg/K) C* GammaAir Air isentropic coefficient (-) C* C* WATER PROPERTIES C* CpWat Specific heat of liquid water (J/kg/K) 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: It is assumed that the water-environment C heat transfer coefficient as well as the C nozzle throat area, the internal efficiency C ,the fuel/air ratio and the torque C associated with the mechanical losses and C the auxiliary consumptions are constant. C Air-fuel mixing properties are the same as C for pure air. C The gas-water heat transfer coefficient is C function of the flue gas mass flow rate. C C DEVELOPER: Jean Lebrun C Marc Grodent C Jean-Pascal Bourdouxhe C Mark Nott C University of LiŠge, Belgium C C DATE: March 1, 1995 C C SUBROUTINES CALLED: TYPE99 (COMBCH) C ENTHALP C FUEL C LINKCK C*********************************************************************** C INTERNAL VARIABLES C Nmin Minimum rotation speed (1/s) C Nmax Maximum rotation speed (1/s) C VfrCyl Volume flow rate corresponding to all the (m**3/s) C cylinders C vCyl Specific volume at the cylinder supply (m**3/kg) C p3 Pressure at the cylinder supply (Pa) C pcritic Critical pressure (Pa) C v1 Air specific volume (m**3/kg) C Wpumping Pumping losses (W) C Win Internal power (W) C Wlo Power associated with the mechanical losses (W) C and the auxiliary consumptions C hg0f Flue gas enthalpy at the exhaust of the (J/kg gas) C adiabatic combustion chamber C hg0 Flue gas enthalpy at the exhaust of (J/kg flue gas) C the adiabatic combustion chamber C hgsu Flue gas enthalpy at the supply of (J/kg flue gas) C the gas-water heat exchanger C Tg0 Flue gas temperature at the exhaust of the (K) C adiabatic combustion chamber C Tgsu Flue gas temperature at the supply of the (K) C gas-water heat exchanger C Twexs Water temperature at the flue gas-water heat (K) C exchanger exhaust C TolRel Relative error tolerance (-) C Crgas Capacity flow rate of the combustion products (W/K) C Crw Water capacity flow rate (W/K) C Fct Value of the function to be nullified (K) C Dfct Value of the first derivative (-) C Effgw Effectiveness of the gas-water heat exchanger (-) C ErrRel Relative error (-) C hgex Gas enthalpy at the exhaust of the (J/kg flue gas) C flue gas-water heat exchanger C Qgw Flue gas-water heat transfer (W) C Qwenv Water-environment heat transfer (W) C AUgwEng Gas-water heat transfer coefficient (W/K) C C Sum1,Sum2,Jm1,Dhgex,DCPgas,Dcrgas,Deffgw,hgcal1,hgcal,Tgsup,hgex1 C and Tgexp are variables used in the Newton-Raphson method. C*********************************************************************** INCLUDE 'c:/Trnsys15/Include/param.inc' INTEGER*4 INFO,INFO99 DOUBLE PRECISION XIN,OUT,XIN99,OUT99 REAL Kmolp(5) REAL Ifuel,MfrW,MfrFuel,MfrGas,N,Nmin,Nmax,i,MfrGasNom DIMENSION PAR(8),XIN(6),OUT(8),INFO(15), & XIN99(5),OUT99(7),INFO99(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) ! Set the version information for TRNSYS if (INFO(7) == -2) then INFO(12) = 15 return 1 endif INFO(6)=8 INFO99(6)=7 DATA TolRel,Nmin,Nmax,CpWat,Pi/1E-05,8,85,4187,3.14159265359/ DATA CpAir,RAir,GammaAir/1005,287.06,1.4/ C*** INPUTS 6 (converted in SI units) C************ Fratio=SNGL(xin(1)) Ta=SNGL(xin(2)+273.15) Mfrw=SNGL(xin(3)/3600.) Twsu=SNGL(xin(4)+273.15) pa=SNGL(xin(5)*101325) N=SNGL(xin(6)/60) C*** PARAMETERS 8 (converted in SI units) C**************** i=par(1) Vs=par(2) Athroat=par(3) EffiInt=par(4) Tlo=par(5) AUgwNoEng=par(6)/3.6 MfrGasNom=par(7)/3600. AUwenvEng=par(8)/3.6 C2*** The gaseous fuel used is methane Ifuel=4 CALL FUEL (Ifuel,Cweight,FLHV,Tr,Cfuel,*1) CALL LINKCK('TYPE905','FUEL',1,99) 1 CONTINUE C1*** Test on the value of the rotation speed IF (N.LT.Nmin) THEN ErrDetec=2 GOTO 90 ELSE IF (N.GT.Nmax) THEN ErrDetec=3 GOTO 90 ENDIF ENDIF C1*** Calculate the critic pressure at the nozzle throat Gm1G=(GammaAir-1)/GammaAir pcritic=pa*(2/(GammaAir+1))**(1/Gm1g) C2*** Calculate the volume flow rate corresponding to all C2*** the cylinders VfrCyl=i*N*Vs C2*** Calculate the pressure at the cylinder supply if we C2*** assumed to be in sonic regime at the nozzle throat v1=Rair*Ta/pa MfrGas=Athroat/v1*SQRT(2*CpAir*Ta)*SQRT((pcritic/pa)** & (2/GammaAir)*(1-(pcritic/pa)**Gm1G)) vCyl=VfrCyl/MfrGas p3=RAir*Ta/vCyl C2*** Compare the pressure at the cylinder supply with the C2*** critic pressure IF (p3.GT.pcritic) THEN C2*** No sonic regime at the nozzle throat; calculate the pressure C2*** at the cylinder supply by means of the Newton-Raphson method C2*** First guess of the value of the pressure at the cylinder supply p3=0.9*pa 5 CONTINUE C2*** Calculate the function to be nullified Fct=Athroat/v1*SQRT(2*CpAir*Ta)*SQRT((p3/pa)** & (2/GammaAir)*(1-(p3/pa)**Gm1G))-VfrCyl*p3/(RAir* & Ta) C2*** Calculate the value of the first derivative prate=p3/pa Den=SQRT(prate**(2/GammaAir)*(1-prate**Gm1G)) DFct=Athroat/v1*SQRT(2*CpAir*Ta)*(2/GammaAir*prate & **((2-GammaAir)/GammaAir)-(GammaAir+1)/ & GammaAir*prate**(1/GammaAir))/(2*pa*Den)- & VfrCyl/(RAir*Ta) C2*** A new estimated value is calculated p3p=p3 p3=p3-Fct/DFct ErrRel=ABS((p3-p3p)/p3p) C2*** If converged, leave the loop IF (ErrRel.GT.TolRel) GOTO 5 IF (p3.LT.pcritic) THEN ErrDetec=4 GOTO 90 ENDIF C2*** Calculate the flue gas mass flow rate vCyl=RAir*Ta/p3 MfrGas=VfrCyl/vCyl ENDIF C1*** Calculate the gas mass flow rate MfrFuel=MfrGas*(Fratio/(1+Fratio)) C2*** Calculate the internal power Win=EffiInt*MfrFuel*FLHV C2*** Calculate the pumping loss Wpumping=i*N*Vs*(pa-p3) C2*** Calculate the mechanical losses and the auxiliary consumptions Wlo=Tlo*2*Pi*N C1*** Calculate the shaft power Wsh=Win-Wpumping-Wlo C1*** Calculate the gas-water heat transfer coefficient AUgwEng=AUgwNoEng*(MfrGas/MfrGasNom)**0.65 C1*** Calculate the adiabatic temperature, the fuel/air ratio as well as C1*** the enthalpy (expressed in J/kg fuel) and composition of the C1*** combustion products xin99(1)=DBLE(Ifuel) xin99(2)=1 xin99(3)=DBLE(Fratio) xin99(4)=DBLE(Ta-273.15) xin99(5)=DBLE(Ta-273.15) CALL TYPE908 (TIME,XIN99,OUT99,T,DTDT,PAR99,INFO99,ICNTRL,*7) CALL LINKCK('TYPE905','TYPE908 ',1,99) 7 CONTINUE Fratio=SNGL(out99(1)) Tg0=SNGL(out99(2)+273.15) Kmolp(2)=SNGL(out99(3)) Kmolp(3)=SNGL(out99(4)) Kmolp(4)=SNGL(out99(5)) Kmolp(5)=SNGL(out99(6)) hg0f=SNGL(out99(7)) C2*** The flue gas enthalpy at the exhaust of the adiabatic C2*** combustion chamber is expressed in J/kg (flue gas) hg0=hg0f/(1+1/Fratio) C1*** Calculate the flue gas enthalpy at the supply of the gas-water C1*** heat exchanger hgsu=hg0-Wsh/MfrGas C1*** Calculate the flue gas temperature at the supply of the C1*** gas-water heat exchanger C2*** First guess of the flue gas temperature at the heat exchanger C2*** supply Tgsu=Tg0/2 10 hgcal1=0 DO 20 J=2,5 CALL ENTHALP (Tgsu,J,hpi,*11) CALL LINKCK('TYPE905','ENTHALP',1,99) 11 CONTINUE hgcal1=hgcal1+Kmolp(J)*hpi 20 CONTINUE hgcal=hgcal1/(1+1/Fratio) C2*** Calculate the function to nullify Fct=hgcal-hgsu 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)*Tgsu**(J-1) 40 CONTINUE Sum1=Sum1+Kmolp(K)*Sum2 30 CONTINUE DFct=Sum1/(1+1/Fratio) C2*** A new estimated value is calculated Tgsup=Tgsu Tgsu=Tgsu-Fct/DFct ErrRel=ABS((Tgsu-Tgsup)/Tgsup) C2*** If converged, then leave the loop IF (ErrRel.GT.TolRel) GOTO 10 C2*** First guess of the exhaust flue gas temperature Tgex=Tgsu/2 C1*** Calculate the exhaust flue gas enthalpy (expressed in J/kg fuel) 50 hgex1=0 DO 60 J=2,5 CALL ENTHALP (Tgex,J,hpi,*51) CALL LINKCK('TYPE905','ENTHALP',1,99) 51 CONTINUE hgex1=hgex1+Kmolp(J)*hpi 60 CONTINUE C2*** The exhaust flue gas enthalpy is expressed in J/kg gas hgex=hgex1/(1+1/Fratio) C1*** Calculate the flue gas mean specific heat CPgas=(hgsu-hgex)/(Tgsu-Tgex) C1*** Calculate a new estimated value of the exhaust flue gas C1*** temperature by using the Newton-Raphson method C2*** Calculate the value of the function to be nullified Crgas=MfrGas*CPgas Crw=MfrW*CpWat C1*** Determine the value of ErrDetec IF (Crgas.GT.Crw) THEN ErrDetec=1 GOTO 90 ELSE ErrDetec=0 ENDIF par1=EXP(-AUgwEng*(1/Crgas-1/Crw)) Effgw=(1-par1)/(1-Crgas*par1/Crw) Fct=Effgw*(Tgsu-Twsu)-Tgsu+Tgex C2*** Calculate the value of the first derivative Sum1=0 DO 70 K=2,5 Sum2=0 DO 80 J=1,10 Jm1=J-1 Sum2=Sum2+PFCP(K,J)*Tgex**Jm1 80 CONTINUE Sum1=Sum1+Sum2*Kmolp(K) 70 CONTINUE Dhgex=Sum1/(1+1/Fratio) DCPgas=(hgsu-hgex-Dhgex*(Tgsu-Tgex))/(Tgsu-Tgex)**2 DCrgas=MfrGas*DCPgas DEffgw=(AUgwEng*DCrgas*par1*(1/Crw-1/Crgas)/Crgas+DCrgas*par1* & (1-par1)/Crw)/(1-(Crgas/Crw)*par1)**2 Dfct=(Tgsu-Twsu)*DEffgw+1 Tgexp=Tgex C2*** The new estimated value is calculated Tgex=Tgex-Fct/Dfct ErrRel=ABS((Tgex-Tgexp)/Tgexp) C2*** If converged, leave loop IF (ErrRel.GT.TolRel) GO TO 50 C1*** Calculate the gas-water heat transfer Qgw=MfrGas*(hgsu-hgex) C1*** Calculate the exhaust water temperature Twexs=Twsu+Qgw/(MfrW*CpWat) Twex=Ta+(Twexs-Ta)/EXP(AUwenvEng/(MfrW*CpWat)) C1*** Calculate the water-environment heat transfer Qwenv=MfrW*CpWat*(Twexs-Twex) C1*** Calculate the gas engine efficiency Effic=Wsh/(MfrFuel*FLHV) 90 CONTINUE C*** OUTPUTS 8 (converted in TRNSYS units) C************* out(1)=DBLE(CPgas) out(2)=DBLE(Twex-273.15) out(3)=DBLE(Tgex-273.15) out(4)=DBLE(Wsh*3.6) out(5)=DBLE(MfrFuel*3600.) out(6)=DBLE(MfrGas*3600.) out(7)=DBLE(Effic) out(8)=DBLE(ErrDetec) RETURN 1 END