[TRNSYS-users] linking subroutines in trnsys16
Jeroen Van der Veken
Jeroen.VanderVeken at bwk.kuleuven.ac.be
Fri Jun 17 06:39:22 PDT 2005
Dear all,
I tried to use a type that calls different subroutines. I get the
following error, indicating that it did not link the subroutines well.
I do not see the solution to this error.
I attached the fortran files as well.
*** Fatal Error at time : 0.025000
Generated by Unit : Not applicable or not available
Generated by Type : Not applicable or not available
TRNSYS Message 104 : The TRNSYS processor has reported that a
subroutine was called that has not been found in the available TRNSYS
libraries.
Reported information : Reported by LINKCK
***** ERROR ***** TRNSYS ERROR # 104
TYPE905 FUEL REQUIRES THE FILE "ENTHALP " WHICH WAS CALLED BUT NOT
LINKED.
PLEASE LINK IN THE REQUIRED FILE AND RERUN THE SIMULATION.
Kind regards,
Leen
--
Ir.Jeroen Van der Veken
Afdeling Bouwfysica
Katholieke Universiteit Leuven
Kasteelpark Arenberg 40
3001 Heverlee
T: +32 16 32 13 47
F: +32 16 32 19 80
@: jeroen.vanderveken at bwk.kuleuven.be
NEW MAILADRES WITHOUT .AC !!
-------------- next part --------------
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
-------------- next part --------------
SUBROUTINE ENTHALP (Temp,I,Enthalpy,*)
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: November 8, 1993
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***********************************************************************
!export this subroutine for its use in external DLLs.
!DEC$ATTRIBUTES DLLEXPORT :: ENTHALP
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
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
More information about the TRNSYS-users
mailing list