Dear TRNSYS users
I tried to modify a trnsys component (type 160) by introducing a new electrochemical model. I have followed the procedure and I recived this message:
“Fatal Error at time : 4001.000000
Generated by Unit : 8
Generated by Type : 208
TRNSYS Message 6 : The TRNSYS Program attempted to evaluate a logarithm (base 10) on a value that is less than or equal to zero. Please check the EQUATION formulation and re-run the simulation
Reported information: Not available”
Then I tried to create a new component that have the same program as Type 160 (the original one) and I received the same message.
Attached is the program of type206 the modified version of type160.
SUBROUTINE TYPE206(TIME,XIN,OUT,T,DTDT,PAR,INFO,ICNTRL,*)
C**********************************************************************************************
C ØYSTEIN ULLEBERG
C IFE
C
C TYPE160: ADVANCED ALKALINE ELECTROLYZER - Version 4.5
C
C Versions:
C 1.0 1998 - Originally developed for ØU's PhD-thesis [2]
C 2.0 2000-08-07 - Modified for TRNSYS 14.2 W/IISiBat, ØU
C 3.0 2001-03-14 - Modified for TRNSYS 15.0 w/IISiBat3, ØU
C 4.0 2001-10-10 - Modifications for TYPE160 TRNSYS15.0 w/IISiBat3,
C modularization of TYPE160, new ELYFARAD equation,
C improved thermal models in ELYOFF and ELYTHERM,
C improved source code readability, ØU
C 4.1 2002-01-23 - Bugs in FORMAT statements fixed
C 4.2 2002-04-14 - NUMSTR removed from declaration
C 4.3 2003-10-30 - Checking voltage & current limits are skipped
C 4.4 2004-08-17 - Updated to TRNSYS 16.0 coding standard, DB
C 4.5 2005-03-21 - fixed a bug. This component resets values in its own XIN() array. I
C added an internal array called XIN_INT(), which is filled from the
C normal XIN() array at each timestep. Internal modifications are now
C made to the XIN_INT() array and not to XIN() itself.
C 4.6 2006-07-27 - DAA: Added a check so Type160 does not issue a warning if it is turned
C off and the voltage is low.
C
C *
C MAIN EQUATIONS: *
C --------------- *
C *
C THE MAIN REACTION - DECOMPOSITION OF WATER IS : *
C *
C H2O -> H2 + 1/2 O2 *
C *
C THE MODEL CONSISTS OF THREE MAIN PARTS: *
C *
C 1. THE I-U CHARACTERISTICS USED IS OF THE FORM: *
C *
C U = Urev(T) + r(T)*I/A + s*log10[(t(T)*I/A+1] *
C *
C where: *
C *
C DG(T) *
C Urev(T) = ----- reversible voltage *
C n*F *
C *
C r(T) = r1 + r2*T resistance in electrolyte *
C s = s1 overvoltage on electrodes *
C t(T) = t1 + t2/T + t3/T^2 overvoltage on electrodes *
C *
C P = constant *
C *
C *
C 2. THERMAL ENERGY BALANCE: *
C *
C Qstore = Qgen - Qloss - Qcool *
C *
C where: *
C *
C Qstore = Ct*dT/dt stored heat *
C Qgen = ncells*U*I*[1-eta_e] generated heat *
C Qloss = 1/Rt*(T-Troom) heat losses *
C Qcool = n_H2O*cp_hH2O*(TH2Oin-TH2Oout) cooling *
C *
C *
C 3. OXYGEN & HYDROGEN PRODUCTION: *
C *
C (I/A)^2 *
C Eta_f = ------------ * a2 Faraday efficiency *
C a1 + (I/A)^2 *
C *
C I *
C n_H2 = Eta_f * ncells * --- molar H2 production rate *
C n*F *
C *
C n_O2 = 1/2 * n_H2 molar O2 production rate *
C *
C *
C SUBROUTINES *
C *
C TYPE160 calls the following 6 subroutines: *
C *
C 1. ELYOFF: *
C - Simplified calculations when the electrolyzer is switched OFF *
C *
C 2. ELYGIBBS: *
C - Thermodynamics in a water splitting reaction (per cell basis) *
C *
C 3. ELYELEC: *
C - Electrochemical model calculations (per cell basis) *
C *
C 4. ELYFARAD: *
C - Faraday efficiency calculations *
C *
C 5. ELYSTACK *
C - Overall stack calculations *
C *
C 6. ELYTHERM *
C - Thermal energy equation calculaltions (per stack basis) *
C *
C *
C ASSUMPTIONS: *
C See reference [1]. *
C *
C REFERENCES: *
C 1. Ulleberg Ø (2001) Modeling of advanced alkaline electrolyzers: *
C A systems Simulation Approach. Paper to be published in 2001. *
C 2. Ulleberg Ø (1998) Stand-Alone Power Systems for the Future: *
C Optimal Design, Operation & Control of Solar-Hydrogen Energy *
C Systems, PhD thesis. Norwegian University of Science and *
C Technology, Trondheim. ISBN 82-471-0344-3. *
C *
C *
C PARAMETERS - TRNSYS DECK: *
C 1. TMODE [1-3] TEMPERATURE MODE *
C MODE 1 = T IS GIVEN AS INPUT XIN(7) *
C MODE 2 = T FROM SIMPLE THERMAL MODEL *
C MODE 3 = T FROM COMPLEX THERMAL MODEL *
C 2. AREA [m^2] AREA OF ELECTRODES *
C 3. NCELLS [#] NUMBER OF CELLS IN SERIES *
C 4. NSTACKS [#] NUMBER OF STACKS IN PARALLEL PER UNIT *
C 5. IDMAX [mA/cm^2] MAXIMUM ALLOWABLE CURRENT DENSITY *
C 6. TMAX [°C] MAXIMUM ALLOWABLE OPERATING TEMPERATURE *
C 7. UCMIN [V/cell] MINIMUM ALLOWABLE OPERATING CELL VOLTAGE *
C 8. R_T [K/W] THERMAL RESISTANCE *
C 9. TAU_T [hr] THERMAL TIME CONSTANT *
C 10. ELYTYPE [#] TYPE OF ELECTOLYZER *
C 11. LUD [#] LOGICAL UNIT OF PARAMETER FILE *
C *
C PARAMETERS - EXTERNAL FILE: *
C 1. R1 [ohm m^2] OHMIC RESISTANCE *
C 2. R2 [ohm m^2/C] OHMIC RESISTANCE *
C 3. S1 [V] OVERVOLTAGE ON ELECTRODES *
C 4. T1 [m^2/A] OVERVOLTAGE ON ELECTRODES *
C 5. T2 [m^2 °C/A] OVERVOLTAGE ON ELECTRODES *
C 6. T3 [m^2 °C^2/A] OVERVOLTAGE ON ELECTRODES *
C 7. A1 [mA/cm2] FARADAY EFFICIENCY *
C 8. A2 [0...1] FARADAY EFFICIENCY *
C 9. HX1 [W/°C] HEAT EXCHANGER COEFFICIENT *
C 10. HX2 [W/°C per A] HEAT EXCHANGER COEFFICIENT *
C *
C INPUTS: *
C 1. SWITCH [0 or 1] ELECTROLYZER OPERATING SWITCH *
C 2. IELY [A] CURRENT THROUGH ELECTROLYZER *
C 3. PELY [bar] PRESSURE *
C 4. TROOM [°C] AMBIENT (ROOM) TEMPERATURE *
C 5. TCOOLIN [°C] TEMPERATURE - COOLING WATER INLET *
C 6. VCOOL [Nm3/h] COOLING WATER FLOW RATE *
C 7. TELY/TINI [°C] TEMERATURE - ELECTROLYZER *
C (TELY->TCMODE=1, TINI->TCMODE=2 OR 3) *
C *
C VARIABLES: *
C T [°C] TEMPERATURE *
C P [bar] PRESSURE *
C *
C SUBSCRIPTS: *
C COOL COOLING WATER *
C ELY ELECTROLYZER *
C 0,STD STANDARD CONDITIONS *
C REF REFERENCE *
C TOT TOTAL/OVERALL *
C H2 HYDROGEN *
C O2 OXYGEN *
C H2O WATER *
C GAS IDEAL GAS *
C FLAG FOR ITERATION PURPOSES *
C DIFF FOR ITERATION PURPOSES *
C OLD FOR ITERATION PURPOSES *
C INI VALUE AT BEGINNING OF TIMESTEP *
C *
C OUTPUTS: *
C 1. IELY [A] CURRENT DRAWN BY ELECTROLYZER *
C 2. UELY [V] TOTAL VOLTAGE ACROSS ELECTROLYZER *
C 3. PTOT [W] TOTAL POWER DRAWN BY ELECTROLYZER *
C 4. VDOT_H2 [Nm3/hr] HYDROGEN PRODUCTION RATE *
C 5. VDOT_O2 [Nm3/hr] OXYGEN PRODUCTION RATE *
C 6. ETA_TOT [0...1] OVERALL EFFICIENCY *
C 7. ETA_E [0...1] ENERGY EFFICIENCY *
C 8. ETA_F [0...1] FARADAY EFFICIENCY *
C 9. QGEN [W] HEAT GENERATED BY ELECTROLYZER *
C 10. QLOSS [W] THERMAL HEAT LOSS *
C 11. QCOOL [W] AUXILARY COOLING *
C 12. QSTORE [W] ENERGY STORAGE RATE *
C 13. TELY [°C] TEMERATURE - ELECTROLYZER *
C 14. TCOOLOUT [°C] TEMPERATURE - COOLING WATER OUTLET *
C 15. IDENSITY [mA/cm2] CURRENT DENSITY *
C 16. UCELL [V/CELL] VOLTAGE ACROSS A SINGLE CELL *
C 17. UREV [V/CELL] REVERSIBLE VOLTAGE *
C 18. UTN [V/CELL] THERMONEUTRAL VOLTAGE *
C 19. void for testing purposes *
C 20. void for testing purposes *
C *
C**********************************************************************************************
!export this subroutine for its use in external DLLs.
!DEC$ATTRIBUTES DLLEXPORT :: TYPE206
C USE STATEMENTS
USE TrnsysFunctions
USE TrnsysConstants
IMPLICIT NONE !force explicit declaration of local variables
C TRNSYS DECLARATIONS
DOUBLE PRECISION XIN,OUT,TIME,PAR,T,DTDT,STORED
INTEGER*4 INFO(15),NPMAX,NI,NO,ND,NS,IUNIT,ITYPE,ICNTRL
CHARACTER*3 OCHECK,YCHECK
C SET THE MAXIMUM NUMBER OF PARAMETERS(NP),INPUTS(NIMAX),OUTPUTS(NO),AND DERIVATIVES(ND)
C THAT MAY BE SUPPLIED FOR THIS TYPE
PARAMETER (NPMAX=20,NI=7,NO=18,ND=0,NS=2)
C REQUIRED TRNSYS DIMENSIONS
DIMENSION XIN(NI),OUT(NO),PAR(NPMAX),YCHECK(NI),OCHECK(NO),
& STORED(NS)
C COMMON BLOCKS
!none required by this Type
C LOCAL VARIABLE DECLARATIONS
DOUBLE PRECISION ELYPAR(10)
INTEGER TMODE,NSTACKS,NELY,ELYTYPE,LUD,I,NP
CHARACTER (len=60) ELYTEXT
CHARACTER (len=30) IDENSITYStr,IDMAXStr,UCELLStr,UCMINStr,
1TELYStr,TMAXStr
CHARACTER (len=maxMessageLength) T160Msg
DOUBLE PRECISION IDMAX,TMAX,UCMIN
DOUBLE PRECISION R1,R2,S1,T1,T2,T3,A1,A2,HX1,HX2
DOUBLE PRECISION IDENSITY,UCELL,IELY,TELY,TINI,TNEW
DOUBLE PRECISION BETA,DIF,TIME0,DELT,DIFF,XIN_INT(NI)
INTEGER*4 IS,TDEBUG1,TDEBUG2,LUW
LOGICAL DeckPars,LUPars
C DATA STATEMENTS
DATA IUNIT/0/
DATA TDEBUG1/-2/,TDEBUG2/-1/
DATA YCHECK/'DM1','EL1','PR1','TE1','TE1','SF1','TE1'/
DATA OCHECK/'EL1','EL2','PW2','SF1','SF1','DM1','DM1','DM1','PW2',
& 'PW2','PW2','PW2','TE1','TE1','EL3','EL4','EL4','EL4'/
C--------------------------------------------------------------------------------------------------
C GET GLOBAL TRNSYS SIMULATION VARIABLES
TIME0 = getSimulationStartTime()
DELT = getSimulationTimeStep()
LUW = getListingFileLogicalUnit()
C--------------------------------------------------------------------------------------------------
C SET THE VERSION INFORMATION FOR TRNSYS
IF(INFO(7).EQ.-2) THEN
INFO(12)=16
RETURN 1
ENDIF
C--------------------------------------------------------------------------------------------------
C PERFORM LAST CALL MANIPULATIONS
IF (INFO(8).EQ.-1) THEN
RETURN 1
ENDIF
C--------------------------------------------------------------------------------------------------
C PERFORM POST CONVERGENCE MANIPULATIONS
IF(INFO(13).GT.0) THEN
CALL getStorageVars(STORED,NS,INFO)
STORED(1)=STORED(2)
CALL setStorageVars(STORED,NS,INFO)
RETURN 1
ENDIF
C--------------------------------------------------------------------------------------------------
C PERFORM FIRST CALL MANIPULATIONS
IF (INFO(7).EQ.-1) THEN
!retrieve unit and type number for this component from the INFO array
IUNIT=INFO(1)
ITYPE=INFO(2)
!reserve space in the OUT array using INFO(6)
INFO(6) = NO
!set the way in which this component is to be called
INFO(9) = 1
!set required number of spots in the single precision storage structure
INFO(10) = 0
!reserve required number of spots in the double precision storage structure
CALL setStorageSize(NS,INFO)
!check for the correct number of PARAMETERS.
IF (INFO(4).EQ.11) THEN !some parameters are contained in the input file, others in an external file.
DeckPars = .FALSE.
LUPars = .TRUE.
NP=11
ELSEIF (INFO(4).EQ.19) THEN !all parameters are contained in the input file.
DeckPars = .TRUE.
LUPars = .FALSE.
NP=19
ELSE
CALL TYPECK(4,INFO,0,0,0)
ENDIF
!call the type check routine to compare what this type requires to what is in the input file
CALL TYPECK(1,INFO,NI,NP,ND)
!call the input-output check subroutine to set the correct input and output units.
CALL RCHECK(INFO,YCHECK,OCHECK)
!RETURN TO THE CALLING PROGRAM
RETURN 1
ENDIF
C--------------------------------------------------------------------------------------------------
C PERFORM INITIAL TIMESTEP MANIPULATIONS
IF (TIME.LT.(TIME0+DELT/2.)) THEN
!set the UNIT number for future calls
IUNIT=INFO(1)
ITYPE=INFO(2)
!determine whether parameters have been provided in the deck or in the external file.
IF (INFO(4).EQ.11) THEN
DeckPars = .FALSE.
LUPars = .TRUE.
ELSEIF (INFO(4).EQ.19) THEN
DeckPars = .TRUE.
LUPars = .FALSE.
ENDIF
!read the parameter values that are always contained in the deck.
!NOTE: not all parameters are required by the main routine.
TMODE = JFIX(PAR(1)+0.1d0)
NSTACKS = JFIX(PAR(4)+0.1d0)
IDMAX = PAR(5)
TMAX = PAR(6)
UCMIN = PAR(7)
!read the parameters that may be in the deck and may be in the external file.
IF (DeckPars.EQ..FALSE.) THEN !the remaining parameters are in an external file
ELYTYPE = JFIX(PAR(10)+0.1)
LUD = JFIX(PAR(11)+0.1)
NELY = 10
!get the remaining parameters from the external file.
CALL PARREAD(NELY,LUD,ELYTYPE,ELYPAR,ELYTEXT,INFO)
R1 = ELYPAR(1)
R2 = ELYPAR(2)
S1 = ELYPAR(3)
T1 = ELYPAR(4)
T2 = ELYPAR(5)
T3 = ELYPAR(6)
A1 = ELYPAR(7)
A2 = ELYPAR(8)
HX1 = ELYPAR(9)
HX2 = ELYPAR(10)
ELSE !the remaining parameters are in the deck
R1 = PAR(10) ; ELYPAR(1) = R1
R2 = PAR(11) ; ELYPAR(2) = R2
S1 = PAR(12) ; ELYPAR(3) = S1
T1 = PAR(13) ; ELYPAR(4) = T1
T2 = PAR(14) ; ELYPAR(5) = T2
T3 = PAR(15) ; ELYPAR(6) = T3
A1 = PAR(16) ; ELYPAR(7) = A1
A2 = PAR(17) ; ELYPAR(8) = A2
HX1 = PAR(18) ; ELYPAR(9) = HX1
HX2 = PAR(19) ; ELYPAR(10) = HX2
ELYTEXT = 'Electrolyzer description unavailable'
ENDIF
!write electrolyzer information to the list file
WRITE(LUW,110) INFO(1),ELYTEXT,R1,R2,S1,T1,T2,T3,A1,A2,HX1,HX2
!check the parameters for problems and RETURN if any are found
IF ((TMODE.LT.1).OR.(TMODE.GT.3)) CALL TYPECK(-4,INFO,0,1,0)
!set initial values of the global storage array variables
STORED(1:NS) = XIN(7)
CALL setStorageVars(STORED,NS,INFO)
!set initial output values
OUT(1:NO)=0.d0
RETURN 1
ENDIF
C--------------------------------------------------------------------
C THIS IS AN ITERATIVE CALL TO THIS COMPONENT ***
C--------------------------------------------------------------------
C RE-READ THE PARAMETERS IF ANOTHER UNIT OF THIS TYPE HAS BEEN CALLED SINCE THE LAST
C TIME THEY WERE READ IN
IF(INFO(1).NE.IUNIT) THEN
IUNIT = INFO(1)
ITYPE = INFO(2)
!reread parameter values (those contained in TRNSYS input file)
!determine whether parameters have been provided in the deck or in the external file.
IF (INFO(4).EQ.11) THEN
DeckPars = .FALSE.
LUPars = .TRUE.
ELSEIF (INFO(4).EQ.19) THEN
DeckPars = .TRUE.
LUPars = .FALSE.
ENDIF
!read the parameter values that are always contained in the deck.
!NOTE: not all parameters are required by the main routine.
TMODE = JFIX(PAR(1)+0.1d0)
NSTACKS = JFIX(PAR(4)+0.1d0)
IDMAX = PAR(5)
TMAX = PAR(6)
UCMIN = 0
!read the parameters that may be in the deck and may be in the external file.
IF (DeckPars.EQ..FALSE.) THEN !the remaining parameters are in an external file
ELYTYPE = JFIX(PAR(10)+0.1)
LUD = JFIX(PAR(11)+0.1)
NELY = 10
!get the remaining parameters from the external file.
CALL PARREAD(NELY,LUD,ELYTYPE,ELYPAR,ELYTEXT,INFO)
R1 = ELYPAR(1)
R2 = ELYPAR(2)
S1 = ELYPAR(3)
T1 = ELYPAR(4)
T2 = ELYPAR(5)
T3 = ELYPAR(6)
A1 = ELYPAR(7)
A2 = ELYPAR(8)
HX1 = ELYPAR(9)
HX2 = ELYPAR(10)
ELSE !the remaining parameters are in the deck
R1 = PAR(10) ; ELYPAR(1) = R1
R2 = PAR(11) ; ELYPAR(2) = R2
S1 = PAR(12) ; ELYPAR(3) = S1
T1 = PAR(13) ; ELYPAR(4) = T1
T2 = PAR(14) ; ELYPAR(5) = T2
T3 = PAR(15) ; ELYPAR(6) = T3
A1 = PAR(16) ; ELYPAR(7) = A1
A2 = PAR(17) ; ELYPAR(8) = A2
HX1 = PAR(18) ; ELYPAR(9) = HX1
HX2 = PAR(19) ; ELYPAR(10) = HX2
ELYTEXT = 'Electrolyzer description unavailable'
ENDIF
ENDIF
C---START OF MAIN PROGRAM----------------------------------------------
C This component resets values in its own XIN() array for internal calculations. So that the
C global values in the XIN() array are not modified, the contents of the XIN() array need to
C be read into an internal array for those internal calulations.
XIN_INT(1:NI) = XIN(1:NI)
C Retrieve stored values for this component
CALL getStorageVars(STORED,NS,INFO)
C Get the Values of the Inputs
IELY = XIN(2)/NSTACKS !Current per stack
XIN_INT(2) = IELY !Resetting current for internal calculations
TINI = STORED(1) !Initial electrolyzer temperature
TELY = TINI !Temperature (first guess)
C---ELECTROLYZER IS SWITCHED 'OFF'---
IF ((JFIX(XIN_INT(1)+0.01) .EQ. 0).OR.(IELY.LT.1E-6)) THEN
CALL ELYOFF(ELYPAR,PAR,XIN_INT,OUT,DELT,TELY,INFO)
C---ELECTROLYZER IS SWITCHED 'ON'---
ELSE
C Thermal steady state operation (TMODE=1)
IF (TMODE.EQ.1) THEN
TINI = XIN_INT(7) !T is given
TELY = TINI !constant T
CALL ELYGIBBS(ELYPAR,PAR,XIN_INT,OUT,TELY,INFO) !Thermodynamics
IF (ErrorFound () ) RETURN 1
CALL ELYELEC(ELYPAR,PAR,XIN_INT,OUT,TELY,INFO) !Electrochemistry
IF (ErrorFound () ) RETURN 1
CALL ELYFARAD(ELYPAR,PAR,XIN_INT,OUT) !Efficiency calcs.
CALL ELYSTACK(ELYPAR,PAR,XIN_INT,OUT) !Stack calculations
CALL ELYTHERM(ELYPAR,PAR,XIN_INT,OUT,DELT,TINI,TELY) !Thermal energy calcs
C Thermal non-steady state operation (TMODE=2 or 3)
ELSEIF ((TMODE.EQ.2).or.(TMODE.EQ.3)) THEN
BETA = 0.6
DIFF = 10.
I = 0
DO WHILE ((DIFF.gt.0.0001).and.(I.lt.100))
CALL ELYGIBBS(ELYPAR,PAR,XIN_INT,OUT,TELY,INFO) !Thermodynamics
IF (ErrorFound () ) RETURN 1
CALL ELYELEC(ELYPAR,PAR,XIN_INT,OUT,TELY,INFO) !Electrochemistry
IF (ErrorFound () ) RETURN 1
CALL ELYFARAD(ELYPAR,PAR,XIN_INT,OUT) !Efficiency
CALL ELYSTACK(ELYPAR,PAR,XIN_INT,OUT) !Stack calcs
CALL ELYTHERM(ELYPAR,PAR,XIN_INT,OUT,DELT,TINI,TELY) !Thermal energy
TNEW = OUT(13) !T from ELYTHERM
DIFF = DABS(TELY-TNEW) !T difference
TELY = BETA*TNEW+(1-BETA)*TELY !Resetting T
I = I+1
END DO
ENDIF
ENDIF
C Warnings
IDENSITY = OUT(15)
UCELL = OUT(16)
IF (IDENSITY.GT.IDMAX) THEN
WRITE(IDENSITYStr,*) JFIX(IDENSITY+0.1d0)
WRITE(IDMAXStr,*) JFIX(IDMAX+0.1d0)
T160Msg = 'The electrolyzer current ('//TRIM(ADJUSTL(IDENSITYStr
1))//' [mA/cm^2]) is too high. The maximum allowable electrolyzer c
1urrent is '//TRIM(ADJUSTL(IDMAXStr))//' [mA/cm^2].'
CALL MESSAGES(-1,T160Msg,'warning',IUNIT,ITYPE)
ENDIF
IF ((UCELL.LT.UCMIN).and.
& (JFIX(XIN_INT(1)+0.01).NE.0)) THEN
WRITE(UCELLStr,*) UCELL
WRITE(UCMINStr,*) UCMIN
T160Msg = 'The electrolyzer voltage ('//TRIM(ADJUSTL(UCELLStr))
1//' [V]) is too low. The minimum allowable electrolyzer voltage i
1s '//TRIM(ADJUSTL(UCMINStr)) //' [V].'
CALL MESSAGES(-1,T160Msg,'warning',IUNIT,ITYPE)
ENDIF
IF (TELY.GT.TMAX) THEN
WRITE(TELYStr,*) TELY
WRITE(TMAXStr,*) TMAX
T160Msg = 'The electrolyzer temperature ('//TRIM(ADJUSTL(TELYStr
1))//' [C]) is too high. The maximum allowable electrolyzer tempera
1ture is '//TRIM(ADJUSTL(TMAXStr)) //' [C].'
CALL MESSAGES(-1,T160Msg,'warning',IUNIT,ITYPE)
ENDIF
C WRITE DEBUGGING INFORMATION
IF((TIME.GE.TDEBUG1).AND.(TIME.LE.TDEBUG2)) THEN
WRITE(LUW,*) 'UNIT = ',INFO(1)
WRITE(LUW,*) 'TYPE = ',INFO(2)
WRITE(LUW,*) 'TIME = ',TIME
WRITE(LUW,*) 'ITER = ',INFO(7)
WRITE(LUW,*) 'PELYTOT = ',OUT(3)*NSTACKS
WRITE(LUW,*) ' '
ENDIF
C STORE THE FINAL ELECTROLYZER TEMPERATURE AS THE NEXT TIMESTEP'S INITIAL TEMPERATURE.
STORED(2) = OUT(13) !TELY from ELYTHERM
CALL setStorageVars(STORED,NS,INFO)
C---OUTPUTS---
OUT(1) = IELY* NSTACKS ! IELY
OUT(2) = OUT(2) ! USTACK
OUT(3) = OUT(3)* NSTACKS ! PTOT
OUT(4) = OUT(4)* NSTACKS ! VH2
OUT(5) = OUT(5)* NSTACKS ! VO2
OUT(6) = OUT(6) ! ETA_TOT
OUT(7) = OUT(7) ! ETA_E
OUT(8) = OUT(8) ! ETA_F
OUT(9) = OUT(9)* NSTACKS ! QGEN
OUT(10) = OUT(10)* NSTACKS ! QLOSS
OUT(11) = OUT(11)* NSTACKS ! QCOOL
OUT(12) = OUT(12)* NSTACKS ! QSTORE
OUT(13) = OUT(13) ! TELY; final temperature
OUT(14) = OUT(14) ! TCOOLOUT
OUT(15) = OUT(15) ! IDENSITY
OUT(16) = OUT(16) ! UCELL
OUT(17) = OUT(17) ! UREV
OUT(18) = OUT(18) ! UTN
110 FORMAT(//2X,72H***************************************************
.*********************,//30X,13HS U M M A R Y,//24X,23HELECTROLYZER
. PARAMETERS,7H - UNIT,I3,//4X,A50,//4X,5HR1 = ,E10.4,3X,9H[ohm m^2
.],/4X,5HR2 = ,E10.4,3X,11H[ohm m^2/C],/4X,5HS1 = ,E10.4,3X,3H[V],
./4X,5HT1 = ,E10.4,3X,7H[m^2/A],/4X,5HT2 = ,E10.4,3X,10H[m^2 °C/A],
./4X,5HT3 = ,E10.4,3X,12H[m^2 °C^2/A],/4X,5HA1 = ,F12.4,1X,3H[-],/4
.X,5HA2 = ,F12.4,1X,7H[0...1],/4X,6HHX1 = ,F11.4,1X,6H[W/°C],/4X,6H
.HX2 = ,F11.4,1X,12H[W/°C per A],//2X,72H**************************
.**********************************************,/)
150 FORMAT(/16X,25H*********WARNING*********,/16X,4HUNIT,I3,2X,6HTIME
.=, F9.3,/16X,25HINCOMPLETE CONVERGENCE!!!,/16X,25H****************
.*********,/)
RETURN 1
END
C---END OF ELECTROLYZER SUBROUTINE--------------------------------------
C***********************************************************************
C EVERYTHING BELOW THIS LINE BELONGS TO THE ELECTROLYZER MODEL *
C***********************************************************************
SUBROUTINE ELYOFF(ELYPAR,PAR,XIN_INT,OUT,DELT,TINI,INFO)
C***********************************************************************
C *
C ELECTROLYZER SWITCHED OFF *
C *
C NOTE! THERMAL CALCULATIONS ON A PER STACK BASIS *
C *
C PARAMETERS: *
C TMODE TEMPERATURE MODE (1=T GIVEN, 2 & 3 = T CALCULATED) *
C NCELLS NUMBER OF CELLS IN SERIES [#] *
C R_T THERMAL RESISTANCE [K/W] *
C TAU_T THERMAL TIME CONSTANT [hr] *
C *
C INPUTS: *
C TROOM AMBIENT (ROOM) TEMPERATURE [°C] *
C DELT TIME STEP [hr] *
C TINI INITIAL ELECTROLYZER TEMPERATURE [°C] *
C *
C VARIABLES: *
C C_T THERMAL CAPACITY OF ELECTROLYZER [J/K] *
C UREST RESTING VOLTAGE [V] *
C TINI INITIAL ELECTROLYZER TEMPERATURE [°C] *
C *
C OUTPUTS: *
C 1. IELY CURRENT DRAWN BY ELECTROLYZER [A] *
C 2. UELY TOTAL VOLTAGE ACROSS ELECTROLYZER [V] *
C 3. PTOT TOTAL POWER DRAWN BY ELECTROLYZER [W] *
C 4. VDOT_H2 HYDROGEN PRODUCTION RATE [Nm3/hr] *
C 5. VDOT_O2 OXYGEN PRODUCTION RATE [Nm3/hr] *
C 6. ETA_TOT OVERALL EFFICIENCY [-] *
C 7. ETA_E ENERGY EFFICIENCY [-] *
C 8. ETA_F FARADAY EFFICIENCY [-] *
C 9. QGEN HEAT GENERATED BY ELECTROLYZER [W] *
C 10. QLOSS THERMAL HEAT LOSS [W] *
C 11. QCOOL AUXILARY COOLING [W] *
C 12. QSTORE ENERGY STORAGE RATE [W] *
C 13. TELY ELECTROLYZER TEMPERATURE [°C] *
C 14. TCOOLOUT TEMPERATURE - COOLING WATER OUTLET [°C] *
C 15. IDENSITY CURRENT DENSITY [mA/cm2] *
C 16. UCELL VOLTAGE ACROSS A SINGLE CELL [V/CELL] *
C 17. UREV REVERSIBLE VOLTAGE [V/CELL] *
C 18. UTN THERMONEUTRAL VOLTAGE [V/CELL] *
C *
C***********************************************************************
IMPLICIT NONE
DOUBLE PRECISION ELYPAR(10),PAR(11),XIN_INT(7),OUT(20)
DOUBLE PRECISION NCELLS,R_T,TAU_T,C_T,AA,BB
DOUBLE PRECISION QGEN,QCOOL,QLOSS,QSTORE,UREV,UTN,UREST
DOUBLE PRECISION TROOM,TINI,DELT,TELY
DOUBLE PRECISION BETA,DIFF,TOLD,TNEW
INTEGER TMODE,I,INFO(15)
C PARAMETERS
TMODE = JFIX(PAR(1)+0.1)
NCELLS = PAR(3)
R_T = PAR(8)
TAU_T = PAR(9)
C INPUTS
TROOM = XIN_INT(4)
C THERMAL MODEL
C (On a per stack basis)
C_T = TAU_T*3600/R_T !3600 s = 1 hr
QGEN = 0.
QCOOL = 0.
C TEMPERATURE (FINAL)
IF (TMODE.EQ.1) THEN
TELY = TINI
QLOSS = 1/R_T*(TELY-TROOM)
ELSEIF (TMODE.EQ.2) THEN
TELY = TINI
BETA = 0.6
DIFF = 10.
I = 0
DO WHILE ((DIFF.gt.0.0001).and.(I.lt.100))
TOLD = TELY !intermediate T
QLOSS = 1/R_T*(TELY-TROOM)
TELY = TINI+(DELT*3600)/C_T*(QGEN-QCOOL-QLOSS) !1hr=3600s
TNEW = TELY
DIFF = DABS(TNEW-TOLD) !T difference
TELY = BETA*TNEW+(1-BETA)*TELY !Resetting T
I = I+1
END DO
ELSEIF (TMODE.EQ.3) THEN
AA = 1/C_T*(1/R_T)
BB = 1/C_T*(QGEN+TROOM/R_T)
TELY = (TINI-BB/AA)*DEXP(-AA*DELT*3600)+BB/AA !1hr=3600s
QLOSS = 1/R_T*(TELY-TROOM)
ENDIF
QSTORE = QGEN-QCOOL-QLOSS
C RESTING VOLTAGE
CALL ELYGIBBS(ELYPAR,PAR,XIN_INT,OUT,TELY,INFO)
UREV = OUT(17)
UTN = OUT(18)
UREST = UREV
C OUTPUTS
OUT(1) = 0. !IELY
OUT(2) = UREST*NCELLS !UELY
OUT(3) = 0. !PTOT
OUT(4) = 0. !VDOT_H2
OUT(5) = 0. !VDOT_O2
OUT(6) = 0. !ETA_TOT
OUT(7) = 0. !ETA_E
OUT(8) = 0. !ETA_F
OUT(9) = QGEN
OUT(10) = QLOSS
OUT(11) = QCOOL
OUT(12) = QSTORE
OUT(13) = TELY !Final temperature
OUT(14) = XIN_INT(5) !TCOOLOUT=TCOOLIN
OUT(15) = 0. !IDENSITY
OUT(16) = UREST !UCELL
OUT(17) = UREV
OUT(18) = UTN
RETURN
END
SUBROUTINE ELYGIBBS(ELYPAR,PAR,XIN_INT,OUT,TELY,INFO)
C***********************************************************************
C *
C THIS SUBROUTINE CALCULATES THE THERMODYNAMICS IN A *
C WATER SPLITTING REACTON (WATER ELECTROLYSIS): *
C *
C H2O --> 1/2 O2 + H2 *
C *
C *
C DATA: *
C TREF REFERENCE TEMPERATURE [C] *
C PREF REFERENCE PRESSURE [bar] *
C FARAD FARADAY CONSTANT [C/mol] or [As/mol] *
C RGAS UNIVERSAL GAS CONSTANT [J/K-mol] *
C NELEC NUMBER OF ELECTRONS PER MOLECULE OF H2 [#] *
C *
C INPUTS: *
C TELY TEMPERATURE [°C] *
C PELY PRESSURE [bar] *
C *
C VARIABLES: *
C DH ENTHALPY CHANGE FOR DECOMPOSITION OF WATER [J/mol] *
C DS ENTROPY CHANGE FOR DECOMPOSITION OF WATER [J/mol-K] *
C DG CHANGE IN GIBBS ENERGY - DECOMPOSITION OF WATER [J/mol] *
C *
C OUTPUTS: *
C UTN THERMONEUTRAL VOLTAGE [V/cell ] *
C UREV THERMODYNAMIC REVERSIBLE VOLTAGE [V/cell] *
C *
C***********************************************************************
IMPLICIT NONE
DOUBLE PRECISION ELYPAR(10),PAR(11),XIN_INT(7),OUT(20)
DOUBLE PRECISION TELY,PELY,UTN,UREV
DOUBLE PRECISION TREF,PREF,FARAD,RGAS,T_H2,T_O2,T_H2O,DH_H2,
1DH_O2,DH_H2O
DOUBLE PRECISION DH0F_H2O,DH0F_H2,DH0F_O2,S0F_H2O,S0F_H2,S0F_O2
DOUBLE PRECISION S_H2,S_O2,S_H2O,CP0_H2O,CP0_H2,CP0_O2
DOUBLE PRECISION DH,DS,DG
INTEGER NELEC,INFO(15)
DATA TREF/25.0/ !C
DATA PREF/1.0/ !bar
DATA FARAD/96485/ !As/mol
DATA RGAS/8.3145/ !J/K-mol
DATA NELEC/2/ !mol electrons/mol H2
DATA DH0F_H2O/-286E3/, DH0F_H2/0/, DH0F_O2/0/ !J/mol
DATA S0F_H2O/70/, S0F_H2/131/, S0F_O2/205/ !J/K-mol
DATA CP0_H2O/75/, CP0_H2/29/, CP0_O2/29/ !J/K-mol
C INPUTS
PELY = XIN_INT(3)
C TEMPERATURE
T_H2 = TELY
T_O2 = TELY
T_H2O = TELY
C ENTHALPY
DH_H2 = CP0_H2*(T_H2-TREF)+DH0F_H2
DH_O2 = CP0_O2*(T_O2-TREF)+DH0F_O2
DH_H2O = CP0_H2O*(T_H2O-TREF)+DH0F_H2O
DH = DH_H2+0.5*DH_O2-DH_H2O
C ENTROPY - CHECK LOGARITHM ARGUMENTS BEFORE TRYING THEM.
IF (((T_H2+273.15)/(TREF+273.15)).LE.0) THEN
CALL MESSAGES(8,'','fatal',INFO(1),INFO(2))
RETURN
ENDIF
IF ((PELY/PREF).LE.0) THEN
CALL MESSAGES(8,'','fatal',INFO(1),INFO(2))
RETURN
ENDIF
S_H2 = CP0_H2*DLOG((T_H2+273.15)/(TREF+273.15))
& -RGAS*DLOG(PELY/PREF)+S0F_H2
IF (((T_O2+273.15)/(TREF+273.15)).LE.0) THEN
CALL MESSAGES(8,'','fatal',INFO(1),INFO(2))
RETURN
ENDIF
IF ((PELY/PREF).LE.0) THEN
CALL MESSAGES(8,'','fatal',INFO(1),INFO(2))
RETURN
ENDIF
S_O2 = CP0_O2*DLOG((T_O2+273.15)/(TREF+273.15))
& -RGAS*DLOG(PELY/PREF)+S0F_O2
IF (((T_H2O+273.15)/(TREF+273.15)).LE.0) THEN
CALL MESSAGES(8,'','fatal',INFO(1),INFO(2))
RETURN
ENDIF
S_H2O = CP0_H2O*DLOG((T_H2O+273.15)/(TREF+273.15))+S0F_H2O
DS = S_H2+0.5*S_O2-S_H2O
C GIBBS FREE ENERGY
DG = DH-(TELY+273.15)*DS
C THERMONEUTRAL VOLTAGE (PER CELL)
UTN = DH/(NELEC*FARAD)
C REVERSIBLE VOLTAGE (PER CELL)
UREV = DG/(NELEC*FARAD)
C OUTPUTS
OUT(17) = UREV
OUT(18) = UTN
RETURN
END
SUBROUTINE ELYELEC(ELYPAR,PAR,XIN_INT,OUT,TELY,INFO)
C***********************************************************************
C *
C ELECTROCHEMICAL MODEL *
C *
C PARAMETERS: *
C AREA AREA OF ELECTRODES [m^2] *
C R1 OHMIC RESISTANCE [ohm-m^2] *
C R2 " [ohm-m^2/C] *
C R3 " [ohm-m^2/C^2] *
C S1 PARAMETER FOR OVERVOLTAGE ON ELECTRODES [V] *
C S2 " " [V/C] *
C S3 " " [V/C^2] *
C T1 PARAMETER FOR OVERVOLTAGE ON ELECTRODES [m^2/A] *
C T2 " " [m^2/A-C] *
C T3 " " [m^2/A-C^2] *
C *
C INPUTS: *
C TELY TEMPERATURE [C] *
C IELY CURRENT DRAWN BY ELECTROLYZER [A] *
C UREV THERMODYNAMIC REVERSIBLE VOLTAGE [V/cell] *
C *
C OUTPUTS: *
C IDENSITY CURRENT DENSITY [mA/cm2] *
C UCELL VOLTAGE ACROSS A SINGLE CELL [V/CELL] *
C *
C***********************************************************************
IMPLICIT NONE
DOUBLE PRECISION ELYPAR(10),PAR(11),XIN_INT(7),OUT(20),TELY
DOUBLE PRECISION IELY,UREV
DOUBLE PRECISION R1,R2,R3,S1,S2,S3,T1,T2,T3,AREA
DOUBLE PRECISION UOHMIC,UOVERPOT,UCELL,IDENSITY
INTEGER INFO(15)
C PARAMETERS
R1 = ELYPAR(1)
R2 = ELYPAR(2)
S1 = ELYPAR(3)
T1 = ELYPAR(4)
T2 = ELYPAR(5)
T3 = ELYPAR(6)
AREA = PAR(2)
C INPUTS
IELY = XIN_INT(2)
UREV = OUT(17)
C CURRENT-VOLTAGE CHARACTERISTIC (PER CELL)
IF (((T1+T2/TELY+T3/(TELY**2))*IELY/AREA+1).LE.0) THEN
CALL MESSAGES(6,'','fatal',INFO(1),INFO(2))
RETURN
ENDIF
IF (IELY.EQ.0) THEN
UOHMIC = (R1+R2*TELY)*IELY/AREA
UOVERPOT = S1*DLOG10((T1+T2/TELY+T3/(TELY**2))*IELY/AREA+1)
UCELL = UREV+UOHMIC+UOVERPOT
ELSE
UOHMIC = 8.3145*TELY*LOG((IELY+0.1)/0.1)/(2*0.9*96484)
T1= exp(1267*(1/303-1/TELY))
UOHMIC = UOHMIC + IELY*0.2/(AREA*10000*(0.005139*14+0.00326)*T1)
UOVERPOT = (8.3145*TELY/(2*96484))*LOG(1+IELY)
UCELL = UREV+UOHMIC+UOVERPOT
ENDIF
C CURRENT DENSITY
IDENSITY = IELY/AREA/10 !10 mA/cm2 = 1 A/m2
C OUTPUTS
OUT(15) = IDENSITY
OUT(16) = UCELL
RETURN
END
SUBROUTINE ELYFARAD(ELYPAR,PAR,XIN_INT,OUT)
C***********************************************************************
C *
C THIS SUBROUTINE CALCULATES THE FARADAY EFFICIENCY *
C *
C PARAMETERS: *
C A1 PARAMETER FOR FARADAY EFFICIENCY [mA/cm2] *
C *
C INPUTS: *
C IDENSITY CURRENT DENSITY [mA/cm2] *
C *
C OUTPUTS: *
C ETA_F FARADAY EFFICIENCY [0-1] *
C *
C***********************************************************************
IMPLICIT NONE
DOUBLE PRECISION ELYPAR(10),PAR(11),XIN_INT(7),OUT(20)
DOUBLE PRECISION A1,A2,IDENSITY,ETA_F
C PARAMETERS
A1 = ELYPAR(7)
A2 = ELYPAR(8)
C INPUTS
IDENSITY = OUT(15)
C FARADAY EFFICIENCY
ETA_F = IDENSITY**2/(A1+IDENSITY**2)*A2
C OUTPUTS
OUT(8) = ETA_F
RETURN
END
SUBROUTINE ELYSTACK(ELYPAR,PAR,XIN_INT,OUT)
C***********************************************************************
C *
C STACK EFFICIENCY, FLOW RATES, VOLTAGE, AND POWER *
C *
C DATA: *
C TREF REFERENCE TEMPERATURE [C] *
C PREF REFERENCE PRESSURE [bar] *
C FARAD FARADAY CONSTANT [C/mol] or [As/mol] *
C RGAS UNIVERSAL GAS CONSTANT [J/K-mol] *
C NELEC NUMBER OF ELECTRONS PER MOLECULE OF H2 [#] *
C *
C PARAMETERS: *
C NCELLS NUMBER OF CELLS IN SERIES [#] *
C *
C INPUTS: *
C IELY CURRENT THROUGH ELECTROLYZER [A] *
C ETA_F FARADAY EFFICIENCY [0-1] *
C UCELL VOLTAGE ACROSS A SINGLE CELL [V/CELL] *
C UTN THERMONEUTRAL VOLTAGE [V/cell ] *
C *
C VARIABLES: *
C RHO_STD DENSITY OF GAS AT STANDARD CONDITIONS [mol/m3] *
C NDOT_H2 HYDROGEN PRODUCTION RATE [mol/s] *
C NDOT_O2 OXYGEN PRODUCTION RATE [mol/s] *
C *
C OUTPUTS: *
C UELY TOTAL VOLTAGE ACROSS ELECTROLYZER [V] *
C PTOT TOTAL POWER DRAWN BY ELECTROLYZER [W] *
C VDOT_H2 HYDROGEN PRODUCTION RATE [Nm3/hr] *
C VDOT_O2 OXYGEN PRODUCTION RATE [Nm3/hr] *
C ETA_TOT OVERALL EFFICIENCY [-] *
C ETA_E ENERGY EFFICIENCY [-] *
C *
C***********************************************************************
IMPLICIT NONE
DOUBLE PRECISION ELYPAR(10),PAR(11),XIN_INT(7),OUT(20)
DOUBLE PRECISION IELY,UTN,UCELL,ETA_F
DOUBLE PRECISION FARAD,TSTD,PSTD,RGAS
DOUBLE PRECISION RHO_STD,NDOT_H2,NDOT_O2
DOUBLE PRECISION ETA_E,ETA_TOT,VDOT_H2,VDOT_O2,UELY,PTOT
INTEGER NELEC,NCELLS
DATA NELEC/2/ !mol electrons/mol H2
DATA FARAD/96485/ !As/mol
DATA RGAS/8.3145/ !J/K-mol
DATA TSTD/0/ !C
DATA PSTD/1.0/ !bar
C PARAMETERS
NCELLS = JFIX(PAR(3)+0.1d0)
C INPUTS
IELY = XIN_INT(2)
ETA_F = OUT(8)
UCELL = OUT(16)
UTN = OUT(18)
C ENERGY EFFICIENCY
ETA_E = UTN/UCELL
C OVERALL EFFICIENCY
ETA_TOT = ETA_E*ETA_F
C HYDROGEN & OXYGEN PRODUCTION---
NDOT_H2 = ETA_F*NCELLS*IELY/(NELEC*FARAD)
NDOT_O2 = 0.5*NDOT_H2
C PRODUCTION FLOW RATES (STANDARD CUBIC METERS)
RHO_STD = PSTD*1E5/(RGAS*(TSTD+273.15)) !1 bar = 1E5 Pa
VDOT_H2 = NDOT_H2/RHO_STD*3600 !Nm3/hr
VDOT_O2 = NDOT_O2/RHO_STD*3600 !Nm3/hr
C STACK VOLTAGE
UELY = NCELLS*UCELL
C STACK POWER
PTOT = UELY*IELY
C OUTPUTS
OUT(2) = UELY
OUT(3) = PTOT
OUT(4) = VDOT_H2
OUT(5) = VDOT_O2
OUT(6) = ETA_TOT
OUT(7) = ETA_E
RETURN
END
SUBROUTINE ELYTHERM(ELYPAR,PAR,XIN_INT,OUT,DELT,TINI,TELYINI)
C***********************************************************************
C *
C THERMAL ENERGY MODEL (ON A PER STACK BASIS) *
C *
C DATA: *
C TREF REFERENCE TEMPERATURE [C] *
C PREF REFERENCE PRESSURE [bar] *
C FARAD FARADAY CONSTANT [C/mol] or [As/mol] *
C RGAS UNIVERSAL GAS CONSTANT [J/K·mol] *
C NELEC NUMBER OF ELECTRONS PER MOLECULE OF H2 [#] *
C *
C PARAMETERS: *
C TMODE TEMPERATURE MODE (1=T GIVEN, 2 & 3 = T CALCULATED) *
C NCELLS NUMBER OF CELLS IN SERIES [#] *
C R_T THERMAL RESISTANCE [K/W] *
C TAU_T THERMAL TIME CONSTANT [hr] *
C HX1 PARAMETER FOR HEAT EXCHANGER COEFFICIENT [W/°C] *
C HX2 PARAMETER FOR HEAT EXCHANGER COEFFICIENT [W/°C per A] *
C *
C INPUTS: *
C IELY CURRENT THROUGH ELECTROLYZER [A] *
C TROOM AMBIENT (ROOM) TEMPERATURE [C] *
C TCOOLIN TEMPERATURE - COOLING WATER INLET [C] *
C VCOOL COOLING WATER FLOW RATE [Nm3/h] *
C UCELL VOLTAGE ACROSS A SINGLE CELL [V/CELL] *
C UTN THERMONEUTRAL VOLTAGE [V/cell ] *
C DELT TIME STEP [hr] *
C TINI INITIAL ELECTROLYZER TEMPERATURE [°C] *
C TELYINI ELECTROLYZER TEMERATURE - AT ITERATION N [°C] *
C *
C VARIABLES: *
C UA_HX OVERALL HEAT TRANSFER COEFFICIENT f(IELY) [W/°C] *
C CP_H2O HEAT CAPACITY OF COOLING WATER [W/K] *
C C_T THERMAL CAPACITY OF ELECTROLYZER [J/K] *
C *
C OUTPUTS: *
C QGEN HEAT GENERATED BY ELECTROLYZER [W] *
C QLOSS THERMAL HEAT LOSS [W] *
C QCOOL AUXILARY COOLING [W] *
C QSTORE ENERGY STORAGE RATE [W] *
C TELYFIN ELECTROLYZER TEMERATURE - AT ITERATION N+1 [°C] *
C TCOOLOUT TEMPERATURE - COOLING WATER OUTLET [°C] *
C *
C***********************************************************************
IMPLICIT NONE
DOUBLE PRECISION ELYPAR(10),PAR(11),XIN_INT(7),OUT(20)
DOUBLE PRECISION R_T,TAU_T,HX1,HX2,RHO_H2O,M_H2O,CP0_H2O
DOUBLE PRECISION IELY,TROOM,VCOOL,UCELL,UTN
DOUBLE PRECISION QGEN,QLOSS,QCOOL,QSTORE,UA_HX,CP_H2O,C_T,AA,BB
DOUBLE PRECISION TCOOLIN,TCOOLOUT,DELT,TINI,TELYINI,TELYFIN,TELY
INTEGER TMODE,NCELLS
DATA M_H2O/18.016/ !g/mol
DATA RHO_H2O/1000/ !kg/m3
DATA CP0_H2O/75/ !J/K-mol
C PARAMETERS
TMODE = JFIX(PAR(1)+0.1)
NCELLS = PAR(3)
R_T = PAR(8)
TAU_T = PAR(9)
HX1 = ELYPAR(9)
HX2 = ELYPAR(10)
C INPUTS
IELY = XIN_INT(2)
TROOM = XIN_INT(4)
TCOOLIN = XIN_INT(5)
VCOOL = XIN_INT(6)
UCELL = OUT(16)
UTN = OUT(18)
C TEMPERATURE (INITIAL)
TELY = TELYINI
C GENERATED THERMAL ENERGY
QGEN = NCELLS*IELY*(UCELL-UTN)
C HEAT LOSSES TO AMBIENT
QLOSS = 1/R_T*(TELY-TROOM)
C AUXILIARY COOLING DEMAND
CP_H2O = VCOOL*RHO_H2O/M_H2O*CP0_H2O*1000/3600 !1kg=1000g, 1hr=3600s
IF (TMODE.EQ.1) THEN
QCOOL = QGEN-QLOSS
TCOOLOUT = TCOOLIN+QCOOL/CP_H2O
ELSE
UA_HX = HX1+HX2*IELY
TCOOLOUT = TCOOLIN+(TELY-TCOOLIN)*(1-DEXP(-UA_HX/CP_H2O))
QCOOL = CP_H2O*(TCOOLOUT-TCOOLIN)
ENDIF
C THERMAL CAPACITANCE
C_T = TAU_T*3600/R_T !3600 s = 1 hr
C TEMPERATURE (FINAL)
IF (TMODE.EQ.1) THEN
TELYFIN = TINI
ELSEIF (TMODE.EQ.2) THEN
TELYFIN = TINI+(DELT*3600)/C_T*(QGEN-QCOOL-QLOSS) !1hr=3600s
ELSEIF (TMODE.EQ.3) THEN
AA = 1/C_T*(1/R_T+CP_H2O*(1-DEXP(-UA_HX/CP_H2O)))
BB = 1/C_T*(QGEN+TROOM/R_T
& +CP_H2O*(1-DEXP(-UA_HX/CP_H2O))*TCOOLIN)
TELYFIN = (TINI-BB/AA)*DEXP(-AA*DELT*3600)+BB/AA !1hr=3600s
ENDIF
C TOTAL THERMAL ENERGY BALANCE
QSTORE = QGEN-QCOOL-QLOSS
C OUTPUTS
OUT(9) = QGEN
OUT(10) = QLOSS
OUT(11) = QCOOL
OUT(12) = QSTORE
OUT(13) = TELYFIN
OUT(14) = TCOOLOUT
RETURN
END