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