Nodes of different colours represent the following:
Solid arrows point from a parent (sub)module to the submodule which is descended from it. Dashed arrows point from a module being used to the module or program unit using it. Where possible, edges connecting nodes are given different colours to make them easier to distinguish in large graphs.
BEGIN SEQUENTIAL SUBSTITUTION to handle a lot of inter-related calcs
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| integer, | intent(in) | :: | GeneratorNum | |||
| logical, | intent(in) | :: | RunFlag | |||
| real(kind=r64), | intent(in) | :: | myload | |||
| logical, | intent(in) | :: | FirstHVACIteration |
Nodes of different colours represent the following:
Solid arrows point from a procedure to one which it calls. Dashed arrows point from an interface to procedures which implement that interface. This could include the module procedures in a generic interface or the implementation in a submodule of an interface in a parent module. Where possible, edges connecting nodes are given different colours to make them easier to distinguish in large graphs.
Nodes of different colours represent the following:
Solid arrows point from a procedure to one which it calls. Dashed arrows point from an interface to procedures which implement that interface. This could include the module procedures in a generic interface or the implementation in a submodule of an interface in a parent module. Where possible, edges connecting nodes are given different colours to make them easier to distinguish in large graphs.
SUBROUTINE CalcFuelCellGeneratorModel(GeneratorNum,RunFlag,MyLoad,FirstHVACIteration)
! SUBROUTINE INFORMATION:
! AUTHOR Brent Griffith
! DATE WRITTEN Aug 2005
! MODIFIED na
! RE-ENGINEERED na
! PURPOSE OF THIS SUBROUTINE:
! simulate a FuelCell generator using the Annex 42 model
! METHODOLOGY EMPLOYED:
! curve fit of performance data:
! many subdomains such as fuel and air compressors, wa
! REFERENCES: IEA/ECBCS Annex 42....
! USE STATEMENTS:
USE DataHVACGlobals, ONLY : FirstTimeStepSysFlag,TimeStepSys, SysTimeElapsed
USE CurveManager, ONLY : CurveValue
USE ScheduleManager, ONLY : GetCurrentScheduleValue
USE DataHeatBalFanSys, ONLY : ZT
USE DataEnvironment , ONLY: WaterMainsTemp
USE General , ONLY: SolveRegulaFalsi, RoundSigDigits
IMPLICIT NONE
! SUBROUTINE ARGUMENT DEFINITIONS:
INTEGER, INTENT(IN) :: GeneratorNum ! Generator number
LOGICAL, INTENT(IN) :: RunFlag ! TRUE when Generator operating
REAL(r64) , INTENT(IN) :: myload ! Generator demand
LOGICAL, INTENT(IN) :: FirstHVACIteration
! SUBROUTINE PARAMETER DEFINITIONS:
REAL(r64), PARAMETER :: KJtoJ = 1000.d0 !convert Kjoules to joules
! DERIVED TYPE DEFINITIONS
! na
! SUBROUTINE LOCAL VARIABLE DECLARATIONS:
REAL(r64) , save :: PpcuLosses ! losses in inverter [W]
REAL(r64) :: Pel ! DC power generated in Fuel Cell Power Module
REAL(r64) :: Pdemand
REAL(r64) :: Eel
REAL(r64) :: Tavg ! working average temperature
LOGICAL :: constrainedFCPM = .false. ! true if power prod is constrained for some reason
LOGICAL :: ConstrainedFCPMTrans = .false.
REAL(r64) :: PelDiff !
INTEGER :: iter ! loop index over repeating set of inter dependent calculaitons
REAL(r64) :: NdotO2 ! molar rate coeff working varible
REAL(r64) :: CpWater ! heat capacity of water in molar units
REAL(r64) :: WaterEnthOfForm ! Standard molar enthalpy of formation
REAL(r64) :: NdotFuel !fuel flow rate
REAL(r64) :: NdotStoicAir ! Air to match fuel molar rate coeff, working variable
REAL(r64) :: NdotExcessAir ! Air in excess of match for fuel
REAL(r64) :: NdotCO2ProdGas ! CO2 from reaction
REAL(r64) :: NdotH20ProdGas ! Water from reaction
REAL(r64) :: NdotCO2 !temp CO2 molar rate coef product gas stream
REAL(r64) :: NdotN2 !temp Nitrogen rate coef product gas stream
REAL(r64) :: Ndot02 !temp Oxygen rate coef product gas stream
REAL(r64) :: NdotH20 !temp Water rate coef product gas stream
REAL(r64) :: NdotAr !tmep Argon rate coef product gas stream
REAL(r64) :: Cp !temp Heat Capacity, used in thermochemistry units of (J/mol K)
REAL(r64) :: Hmolfuel !temp enthalpy of fuel mixture in KJ/mol
REAL(r64) :: Hmolair !temp enthalpy of air mixture in KJ/mol
REAL(r64) :: HmolProdGases ! enthalpy of product gas mixture in KJ/mol
REAL(r64) :: HLiqWater ! temp enthalpy of liquid water in KJ/mol No Formation
REAL(r64) :: HGasWater ! temp enthalpy of gaseous water in KJ/mol No Formation
INTEGER :: thisGas !loop index
REAL(r64) :: MagofImbalance ! error signal to control exiting loop and targeting product enthalpy
REAL(r64) :: tmpTotProdGasEnthalphy !
REAL(r64) :: Acc ! accuracy control for SolveRegulaFalsi
INTEGER :: MaxIter !iteration control for SolveRegulaFalsi
INTEGER :: SolverFlag !feed back flag from SolveRegulaFalsi
REAL(r64), Dimension(3) :: Par ! parameters passed in to SolveRegulaFalsi
! Par(1) = generator number index in structure
! Par(2) = targeted enthalpy (W)
! Par(3) = molar flow rate of product gases (kmol/s)
REAL(r64) :: tmpTprodGas
!unused REAL(r64) :: LHV !Lower Heating Value
LOGICAL :: ConstrainedStorage ! contrained overall elect because of storage
REAL(r64) :: PgridExtra !extra electric power that should go into storage but can't
REAL(r64) :: Pstorage ! power into storage (+), power from storage (-)
REAL(r64) :: PintoInverter ! power into inverter after storage interactions
REAL(r64) :: PoutofInverter ! power out of inverter after losses and including storage
REAL(r64) :: PacAncillariesTotal ! total AC ancillaries
!! begin controls block to be moved out to GeneratorDynamics module
!If no loop demand or Generator OFF, return
IF (.NOT. Runflag) THEN
! TODO zero out terms as appropriate
If (FuelCell(GeneratorNum)%FCPM%HasBeenOn) Then
!FuelCell just now beginning to shut down,
! set Day and Time of Last Shut Down
FuelCell(GeneratorNum)%FCPM%FractionalDayofLastShutDown = REAL(DayOfSim,r64) &
+ (INT(CurrentTime)+(SysTimeElapsed+(CurrentTime - INT(CurrentTime))))/HoursInDay
FuelCell(GeneratorNum)%FCPM%HasBeenOn = .false.
If (FuelCell(GeneratorNum)%FCPM%ShutDownTime > 0.0d0) FuelCell(GeneratorNum)%FCPM%DuringShutDown = .true.
endif
!TODO check to see if still in shut down mode and using fuel.
If (FuelCell(GeneratorNum)%FCPM%DuringShutDown) Then
endif
RETURN
END IF
IF ( .not. FuelCell(GeneratorNum)%FCPM%HasBeenOn ) Then
!fuel cell just turned on
! set Day and Time of Last STart Up
FuelCell(GeneratorNum)%FCPM%FractionalDayofLastStartUp = REAL(DayOfSim,r64) &
+ (INT(CurrentTime)+(SysTimeElapsed+(CurrentTime - INT(CurrentTime))))/HoursInDay
FuelCell(GeneratorNum)%FCPM%HasBeenOn = .true.
FuelCell(GeneratorNum)%FCPM%NumCycles = FuelCell(GeneratorNum)%FCPM%NumCycles + 1 ! increment cycling counter
If (FuelCell(GeneratorNum)%FCPM%StartUpTime > 0.0d0) FuelCell(GeneratorNum)%FCPM%DuringStartUp = .true.
ENDIF
!TODO deal with things when jump out if not running?
If ( .not. RunFlag) Then
Return
ENDIF
! Note: Myload (input) is Pdemand (electical Power requested)
Pdemand = myLoad
PacAncillariesTotal = 0.0d0
PpcuLosses = 0.0d0
Pstorage = 0.0d0
PgridExtra = 0.0d0
PoutofInverter = 0.0d0
ConstrainedFCPM = .false.
!!BEGIN SEQUENTIAL SUBSTITUTION to handle a lot of inter-related calcs
DO iter=1, 20
IF (iter > 1) then
Call FigurePowerConditioningLosses(GeneratorNum, PoutofInverter , PpcuLosses)
CALL FigureACAncillaries(GeneratorNum, PacAncillariesTotal)
Pdemand = myLoad + PacAncillariesTotal + PpcuLosses
ELSE
! control Step 1a: Figure ancillary AC power draws
CALL FigureACAncillaries(GeneratorNum, PacAncillariesTotal)
Pdemand = myLoad + PacAncillariesTotal
! Control Step 1b: Calculate losses associated with Power conditioning
Call FigurePowerConditioningLosses(GeneratorNum, Pdemand, PpcuLosses)
Pdemand = Pdemand + PpcuLosses
Pel = Pdemand
ENDIF
FuelCell(GeneratorNum)%Inverter%PCUlosses = PpcuLosses
! Control step 2: adjust for transient and startup/shut down constraints
Call FigureTransientConstraints(GeneratorNum, Pel, ConstrainedFCPMTrans, PelDiff)
! Control step 3: adjust for max and min limits on Pel
IF ( Pel < FuelCell(GeneratorNum)%FCPM%PelMin) THEN
PelDiff = PelDiff +(FuelCell(GeneratorNum)%FCPM%PelMin - Pel)
Pel = FuelCell(GeneratorNum)%FCPM%PelMin
ConstrainedFCPM = .true.
ENDIF
IF ( Pel > FuelCell(GeneratorNum)%FCPM%PelMax) THEN
PelDiff = PelDiff + (FuelCell(GeneratorNum)%FCPM%PelMax - Pel)
Pel = FuelCell(GeneratorNum)%FCPM%PelMax
ConstrainedFCPM = .true.
ENDIF
If (ConstrainedFCPM) then
ENDIF
FuelCell(GeneratorNum)%FCPM%Pel = Pel
!Now calculate FC models. return to controls and batter after
!Calculation Step 1. Determine electrical Efficiency Eel
IF (FuelCell(GeneratorNum)%FCPM%EffMode == NormalizedCurveMode) THEN
!Equation (8) in FuelCell Spec modified for normalized curve
Eel = CurveValue(FuelCell(GeneratorNum)%FCPM%EffCurveID, Pel/FuelCell(GeneratorNum)%FCPM%NomPel) &
*FuelCell(GeneratorNum)%FCPM%NomEff &
* (1.0d0 - FuelCell(GeneratorNum)%FCPM%NumCycles * FuelCell(GeneratorNum)%FCPM%CyclingDegradRat) &
* (1.0d0 - MAX((FuelCell(GeneratorNum)%FCPM%NumRunHours-FuelCell(GeneratorNum)%FCPM%ThreshRunHours), 0.0d0) &
*FuelCell(GeneratorNum)%FCPM%OperateDegradRat)
ELSEIF (FuelCell(GeneratorNum)%FCPM%EffMode == DirectCurveMode ) THEN
!Equation (8) in FuelCell Spec
Eel = CurveValue(FuelCell(GeneratorNum)%FCPM%EffCurveID, Pel) &
* (1.0d0 - FuelCell(GeneratorNum)%FCPM%NumCycles * FuelCell(GeneratorNum)%FCPM%CyclingDegradRat) &
* (1.0d0 - MAX((FuelCell(GeneratorNum)%FCPM%NumRunHours-FuelCell(GeneratorNum)%FCPM%ThreshRunHours), 0.0d0) &
*FuelCell(GeneratorNum)%FCPM%OperateDegradRat )
ENDIF
FuelCell(GeneratorNum)%FCPM%Eel = Eel
! Calculation Step 2. Determine fuel rate
NdotFuel = Pel / (Eel * FuelSupply(FuelCell(GeneratorNum)%FuelSupNum)%LHV * 1000000.0d0) !Eq. 10 solved for Ndot
FuelCell(GeneratorNum)%FCPM%NdotFuel = NdotFuel
If (Pel <= 0.0d0) Then
!TODO zero stuff before leaving
Pel = 0.0d0
FuelCell(GeneratorNum)%FCPM%Pel = 0.0d0
RETURN
Else
FuelCell(GeneratorNum)%FCPM%Pel = Pel
endif
! Calculation Step 3. Determine Air rate
IF (FuelCell(GeneratorNum)%AirSup%AirSupRateMode == ConstantStoicsAirRat) THEN !MEthod 1
NdotO2 = FuelSupply(FuelCell(GeneratorNum)%FuelSupNum)%StoicOxygenRate &
* FuelCell(GeneratorNum)%FCPM%NdotFuel &
* FuelCell(GeneratorNum)%AirSup%Stoics
FuelCell(GeneratorNum)%FCPM%NdotAir = NdotO2 / FuelCell(GeneratorNum)%AirSup%O2fraction
ELSEIF (FuelCell(GeneratorNum)%AirSup%AirSupRateMode == QuadraticFuncofPel) THEN !MEthod 2
FuelCell(GeneratorNum)%FCPM%NdotAir = CurveValue(FuelCell(GeneratorNum)%AirSup%AirFuncPelCurveID, Pel) &
* (1 + FuelCell(GeneratorNum)%AirSup%AirTempCoeff &
* FuelCell(GeneratorNum)%AirSup%TairIntoFCPM)
ELSEIF (FuelCell(GeneratorNum)%AirSup%AirSupRateMode == QuadraticFuncofNdot) THEN ! method 3
FuelCell(GeneratorNum)%FCPM%NdotAir = CurveValue(FuelCell(GeneratorNum)%AirSup%AirFuncNdotCurveID, &
FuelCell(GeneratorNum)%FCPM%NdotFuel) &
* (1 + FuelCell(GeneratorNum)%AirSup%AirTempCoeff &
* FuelCell(GeneratorNum)%AirSup%TairIntoFCPM)
ENDIF
! Calculation Step 4. fuel compressor power
FuelSupply(FuelCell(GeneratorNum)%FuelSupNum)%PfuelCompEl &
= CurveValue( FuelSupply(FuelCell(GeneratorNum)%FuelSupNum)%CompPowerCurveID, FuelCell(GeneratorNum)%FCPM%NdotFuel)
! calculation Step 5, Fuel Compressor (need outlet temperature)
If (FuelSupply(FuelCell(GeneratorNum)%FuelSupNum)%FuelTempMode == FuelInTempFromNode) then
FuelSupply(FuelCell(GeneratorNum)%FuelSupNum)%TfuelIntoCompress &
= Node(FuelSupply(FuelCell(GeneratorNum)%FuelSupNum)%NodeNum)%Temp
ELSEIF (FuelSupply(FuelCell(GeneratorNum)%FuelSupNum)%FuelTempMode == FuelInTempSchedule) then
FuelSupply(FuelCell(GeneratorNum)%FuelSupNum)%TfuelIntoCompress &
= GetCurrentScheduleValue(FuelSupply(FuelCell(GeneratorNum)%FuelSupNum)%SchedNum)
ENDIF
! evaluate heat capacity at average temperature usign shomate
Tavg = (FuelSupply(FuelCell(GeneratorNum)%FuelSupNum)%TfuelIntoCompress &
+ FuelSupply(FuelCell(GeneratorNum)%FuelSupNum)%TfuelIntoFCPM)/ 2.0d0
Call FigureFuelHeatCap(GeneratorNum, Tavg, Cp) ! Cp in (J/mol K)
! calculate a Temp of fuel out of compressor and into power module
IF (FuelCell(GeneratorNum)%FCPM%NdotFuel <= 0.0d0) THEN !just pass through, domain probably collapased in modeling
FuelSupply(FuelCell(GeneratorNum)%FuelSupNum)%TfuelIntoFCPM = FuelSupply(FuelCell(GeneratorNum)%FuelSupNum)%TfuelIntoCompress
ELSE
FuelSupply(FuelCell(GeneratorNum)%FuelSupNum)%TfuelIntoFCPM = ( &
(1.0d0 - FuelSupply(FuelCell(GeneratorNum)%FuelSupNum)%CompPowerLossFactor) &
*FuelSupply(FuelCell(GeneratorNum)%FuelSupNum)%PfuelCompEl &
/ (FuelCell(GeneratorNum)%FCPM%NdotFuel * Cp * 1000.0d0)) & !1000 Cp units mol-> kmol
+ FuelSupply(FuelCell(GeneratorNum)%FuelSupNum)%TfuelIntoCompress
ENDIF
! calc skin losses from fuel compressor
FuelSupply(FuelCell(GeneratorNum)%FuelSupNum)%QskinLoss = FuelSupply(FuelCell(GeneratorNum)%FuelSupNum)%CompPowerLossFactor &
*FuelSupply(FuelCell(GeneratorNum)%FuelSupNum)%PfuelCompEl
If (FuelSupply(FuelCell(GeneratorNum)%FuelSupNum)%QskinLoss < 0.0d0) then
! write(*,*) 'problem in FuelSupply%QskinLoss ', FuelSupply(FuelCell(GeneratorNum)%FuelSupNum)%QskinLoss
CALL ShowWarningError('problem in FuelSupply%QskinLoss '// &
TRIM(RoundSigDigits(FuelSupply(FuelCell(GeneratorNum)%FuelSupNum)%QskinLoss,3)))
FuelSupply(FuelCell(GeneratorNum)%FuelSupNum)%QskinLoss = 0.0d0
endif
! calculate tatal fuel enthalpy coming into power module
! (Hmolfuel in KJ/mol)
CAll FigureFuelEnthalpy(GeneratorNum, FuelSupply(FuelCell(GeneratorNum)%FuelSupNum)%TfuelIntoFCPM, Hmolfuel)
! units, NdotFuel in kmol/sec. Hmolfule in KJ/mol ,
! factor of 1000's to get to J/s or watts
FuelCell(GeneratorNum)%FCPM%TotFuelInEnthalphy = Hmolfuel * 1000.0d0 * FuelCell(GeneratorNum)%FCPM%NdotFuel * 1000.0d0
!Calculation Step 6, water compressor calculations
! calculate water consumption
FuelCell(GeneratorNum)%FCPM%NdotLiqwater = CurveValue(FuelCell(GeneratorNum)%WaterSup%WaterSupRateCurveID , &
FuelCell(GeneratorNum)%FCPM%NdotFuel )
! set inlet temp. (could move to init)
SELECT CASE (FuelCell(GeneratorNum)%WaterSup%WaterTempMode)
CASE (WaterInReformMains)
FuelCell(GeneratorNum)%WaterSup%TwaterIntoCompress = WaterMainsTemp
CASE (WaterInReformAirNode,WaterInReformWaterNode)
FuelCell(GeneratorNum)%WaterSup%TwaterIntoCompress = Node(FuelCell(GeneratorNum)%WaterSup%NodeNum)%Temp
CASE (WaterInReformSchedule)
FuelCell(GeneratorNum)%WaterSup%TwaterIntoCompress = GetCurrentScheduleValue(FuelCell(GeneratorNum)%WaterSup%SchedNum)
END SELECT
FuelCell(GeneratorNum)%WaterSup%PwaterCompEl = CurveValue(FuelCell(GeneratorNum)%WaterSup%PmpPowerCurveID, &
FuelCell(GeneratorNum)%FCPM%NdotLiqwater )
! 75.325 J/mol K Water at 0.1 MPa and 298 K, reference NIST WEBBOOK
call FigureLiquidWaterHeatCap(FuelCell(GeneratorNum)%WaterSup%TwaterIntoCompress, CpWater)
WaterEnthOfForm = -241.8264d0 !KJ/mol
IF (FuelCell(GeneratorNum)%FCPM%NdotLiqwater <= 0.0d0) THEN !just pass through, domain probably collapased in modeling
FuelCell(GeneratorNum)%WaterSup%TwaterIntoFCPM = FuelCell(GeneratorNum)%WaterSup%TwaterIntoCompress !
ELSE
FuelCell(GeneratorNum)%WaterSup%TwaterIntoFCPM = ((1 - FuelCell(GeneratorNum)%WaterSup%PmpPowerLossFactor) &
*FuelCell(GeneratorNum)%WaterSup%PwaterCompEl &
/ (FuelCell(GeneratorNum)%FCPM%NdotLiqwater * CpWater * 1000.0d0)) &
+ FuelCell(GeneratorNum)%WaterSup%TwaterIntoCompress
ENDIF
FuelCell(GeneratorNum)%WaterSup%QskinLoss = FuelCell(GeneratorNum)%WaterSup%PmpPowerLossFactor &
* FuelCell(GeneratorNum)%WaterSup%PwaterCompEl
If (FuelCell(GeneratorNum)%WaterSup%QskinLoss < 0.0d0) then
! write(*,*) 'problem in WaterSup%QskinLoss ',FuelCell(GeneratorNum)%WaterSup%QskinLoss
FuelCell(GeneratorNum)%WaterSup%QskinLoss = 0.0d0
endif
Call FigureLiquidWaterEnthalpy(FuelCell(GeneratorNum)%WaterSup%TwaterIntoFCPM, HLiqWater) ! HLiqWater in KJ/mol
FuelCell(GeneratorNum)%FCPM%WaterInEnthalpy = FuelCell(GeneratorNum)%FCPM%NdotLiqwater*HLiqWater*1000.0d0*1000.0d0
!Calculation Step 7, Air compressor
FuelCell(GeneratorNum)%AirSup%TairIntoBlower = Node(FuelCell(GeneratorNum)%AirSup%SupNodeNum)%Temp
FuelCell(GeneratorNum)%AirSup%PairCompEl = CurveValue(FuelCell(GeneratorNum)%AirSup%BlowerPowerCurveID, &
FuelCell(GeneratorNum)%FCPM%NdotAir)
Tavg = ( FuelCell(GeneratorNum)%AirSup%TairIntoBlower + FuelCell(GeneratorNum)%AirSup%TairIntoFCPM ) / 2.0d0
CALL FigureAirHeatCap(GeneratorNum, Tavg, Cp) ! Cp in (J/mol K)
! if PEMFC with stack cooler, then calculate stack cooler impacts
If (FuelCell(GeneratorNum)%StackCooler%StackCoolerPresent) then
FuelCell(GeneratorNum)%StackCooler%qs_cool = ( FuelCell(GeneratorNum)%StackCooler%r0 &
+ FuelCell(GeneratorNum)%StackCooler%r1 &
*(FuelCell(GeneratorNum)%StackCooler%TstackActual &
- FuelCell(GeneratorNum)%StackCooler%TstackNom) ) &
*( 1 + FuelCell(GeneratorNum)%StackCooler%r2*Pel &
+ FuelCell(GeneratorNum)%StackCooler%r3*Pel*Pel) * Pel
FuelCell(GeneratorNum)%FCPM%QdotStackCool = FuelCell(GeneratorNum)%StackCooler%qs_cool
ENDIF
!Figure heat recovery from Electrical Storage, power conditioning, and auxiliary burner
SELECT CASE (FuelCell(GeneratorNum)%AirSup%IntakeRecoveryMode)
CASE(RecoverBurnInvertBatt)
FuelCell(GeneratorNum)%AirSup%QintakeRecovery = FuelCell(GeneratorNum)%AuxilHeat%QairIntake &
+ FuelCell(GeneratorNum)%ElecStorage%QairIntake &
+ FuelCell(GeneratorNum)%Inverter%QairIntake
CASE(RecoverAuxiliaryBurner)
FuelCell(GeneratorNum)%AirSup%QintakeRecovery = FuelCell(GeneratorNum)%AuxilHeat%QairIntake
CASE(RecoverInverterBatt)
FuelCell(GeneratorNum)%AirSup%QintakeRecovery = FuelCell(GeneratorNum)%ElecStorage%QairIntake &
+ FuelCell(GeneratorNum)%Inverter%QairIntake
CASE(RecoverInverter)
FuelCell(GeneratorNum)%AirSup%QintakeRecovery = FuelCell(GeneratorNum)%Inverter%QairIntake
CASE(RecoverBattery)
FuelCell(GeneratorNum)%AirSup%QintakeRecovery = FuelCell(GeneratorNum)%ElecStorage%QairIntake
CASE(NoRecoveryOnAirIntake)
FuelCell(GeneratorNum)%AirSup%QintakeRecovery = 0.0d0
END SELECT
IF ( FuelCell(GeneratorNum)%FCPM%NdotAir <= 0.0d0 ) THEN !just pass through, domain probably collapased in modeling
FuelCell(GeneratorNum)%AirSup%TairIntoFCPM = FuelCell(GeneratorNum)%AirSup%TairIntoBlower
ELSE
FuelCell(GeneratorNum)%AirSup%TairIntoFCPM = (((1 - FuelCell(GeneratorNum)%AirSup%BlowerHeatLossFactor) &
* FuelCell(GeneratorNum)%AirSup%PairCompEl &
+ FuelCell(GeneratorNum)%AirSup%QintakeRecovery) &
/ (FuelCell(GeneratorNum)%FCPM%NdotAir * Cp * 1000.0d0)) & !1000 Cp units mol-> kmol
+ FuelCell(GeneratorNum)%AirSup%TairIntoBlower
ENDIF
FuelCell(GeneratorNum)%AirSup%QskinLoss = FuelCell(GeneratorNum)%AirSup%BlowerHeatLossFactor &
* FuelCell(GeneratorNum)%AirSup%PairCompEl
If (FuelCell(GeneratorNum)%AirSup%QskinLoss < 0.0d0) then
! write(*,*) 'problem in AirSup%QskinLoss ', FuelCell(GeneratorNum)%AirSup%QskinLoss
CALL ShowWarningError('problem in AirSup%QskinLoss '//TRIM(RoundSigDigits(FuelCell(GeneratorNum)%AirSup%QskinLoss,3)))
FuelCell(GeneratorNum)%AirSup%QskinLoss = 0.0d0
endif
CAll FigureAirEnthalpy(GeneratorNum, FuelCell(GeneratorNum)%AirSup%TairIntoFCPM, Hmolair) ! (Hmolair in KJ/mol)
! units, NdotAir in kmol/sec.; Hmolfuel in KJ/mol ,
! factor of 1000's to get to J/s or watts
FuelCell(GeneratorNum)%FCPM%TotAirInEnthalphy = Hmolair * 1000.0d0 * FuelCell(GeneratorNum)%FCPM%NdotAir * 1000.0d0
! calculation Step 8, Figure Product Gases
! figure stoic N dot for air
NdotO2 = FuelSupply(FuelCell(GeneratorNum)%FuelSupNum)%StoicOxygenRate &
* FuelCell(GeneratorNum)%FCPM%NdotFuel
NdotStoicAir = NdotO2 / FuelCell(GeneratorNum)%AirSup%O2fraction
! figure excess air rate
NdotExcessAir = FuelCell(GeneratorNum)%FCPM%NdotAir - NdotStoicAir
IF (NdotExcessAir < 0) THEN !can't meet stoichiometric fuel reaction
CALL ShowWarningError('Air flow rate into fuel cell is too low for stoichiometric fuel reaction')
CALL ShowContinueError('Increase air flow in GENERATOR:FC:AIR SUPPLY object:' &
//Trim(FuelCell(GeneratorNum)%AirSup%Name))
ENDIF
! figure CO2 and Water rate from products (coefs setup during one-time processing in gas phase library )
NdotCO2ProdGas = FuelCell(GeneratorNum)%FCPM%NdotFuel * FuelSupply(FuelCell(GeneratorNum)%FuelSupNum)%CO2ProductGasCoef
NdotH20ProdGas = FuelCell(GeneratorNum)%FCPM%NdotFuel * FuelSupply(FuelCell(GeneratorNum)%FuelSupNum)%H20ProductGasCoef
! set product gas constituent fractions (assume five usual components)
NdotCO2 = 0.0d0
NdotN2 = 0.0d0
Ndot02 = 0.0d0
NdotH20 = 0.0d0
NdotAr = 0.0d0
! Product gas constiuents are fixed (not a user defined thing)
!
Do thisGas = 1, FuelCell(GeneratorNum)%AirSup%NumConstituents
SELECT CASE (FuelCell(GeneratorNum)%AirSup%GasLibID(thisGas))
CASE (1)
! all the CO2 coming in plus the new CO2 from reactions
NdotCO2 = NdotCO2ProdGas + FuelCell(GeneratorNum)%AirSup%ConstitMolalFract(thisGas) * FuelCell(GeneratorNum)%FCPM%NdotAir
CASE (2)
! all the nitrogen comming in
NdotN2 = FuelCell(GeneratorNum)%FCPM%NdotAir * FuelCell(GeneratorNum)%AirSup%ConstitMolalFract(thisGas)
CASE (3)
! all the oxygen in the excess air stream
Ndot02 = NdotExcessAir * FuelCell(GeneratorNum)%AirSup%ConstitMolalFract(thisGas)
CASE (4)
! all the H20 comming in plus the new H20 from reactions and the H20 from water used in reforming
NdotH20 = NdotH20ProdGas + FuelCell(GeneratorNum)%AirSup%ConstitMolalFract(thisGas) * FuelCell(GeneratorNum)%FCPM%NdotAir
!+ FuelCell(GeneratorNum)%FCPM%NdotLiqwater
CASE (5)
! all the argon coming in.
NdotAr = FuelCell(GeneratorNum)%FCPM%NdotAir * FuelCell(GeneratorNum)%AirSup%ConstitMolalFract(thisGas)
CASE DEFAULT
END SELECT
ENDDO
FuelCell(GeneratorNum)%FCPM%NdotProdGas = NdotCO2 + NdotN2 + Ndot02 + NdotH20 + NdotAr
! now that we have the total, figure molar fractions
FuelCell(GeneratorNum)%FCPM%ConstitMolalFract(1) = NdotCO2 / FuelCell(GeneratorNum)%FCPM%NdotProdGas
! all the nitrogen comming in
FuelCell(GeneratorNum)%FCPM%ConstitMolalFract(2) = NdotN2 / FuelCell(GeneratorNum)%FCPM%NdotProdGas
! all the oxygen in the excess air stream
FuelCell(GeneratorNum)%FCPM%ConstitMolalFract(3) = Ndot02 / FuelCell(GeneratorNum)%FCPM%NdotProdGas
! all the H20 comming in plus the new H20 from reactions and the H20 from water used in reforming
FuelCell(GeneratorNum)%FCPM%ConstitMolalFract(4) = NdotH20 / FuelCell(GeneratorNum)%FCPM%NdotProdGas
! all the argon coming in.
FuelCell(GeneratorNum)%FCPM%ConstitMolalFract(5) = NdotAr / FuelCell(GeneratorNum)%FCPM%NdotProdGas
!HmolProdGases KJ/mol)
CALL FigureProductGasesEnthalpy(GeneratorNum, FuelCell(GeneratorNum)%FCPM%TprodGasLeavingFCPM, HmolProdGases)
! units, NdotProdGas in kmol/sec.; HmolProdGases in KJ/mol ,
! factor of 1000's to get to J/s or watts
FuelCell(GeneratorNum)%FCPM%TotProdGasEnthalphy = HmolProdGases * 1000.0d0 * FuelCell(GeneratorNum)%FCPM%NdotProdGas * 1000.0d0
! calculation Step 9, Figure Skin lossess
IF ( FuelCell(GeneratorNum)%FCPM%SkinLossMode == ConstantRateSkinLoss ) THEN
! do nothing just use QdotSkin
ELSEIF ( FuelCell(GeneratorNum)%FCPM%SkinLossMode == UADTSkinLoss ) then
! get zone air temp
FuelCell(GeneratorNum)%FCPM%QdotSkin = FuelCell(GeneratorNum)%FCPM%UAskin * &
( FuelCell(GeneratorNum)%FCPM%TprodGasLeavingFCPM &
- ZT(FuelCell(GeneratorNum)%FCPM%ZoneID ))
ELSEIF ( FuelCell(GeneratorNum)%FCPM%SkinLossMode == QuadraticFuelNdotSkin) then
FuelCell(GeneratorNum)%FCPM%QdotSkin = CurveValue( FuelCell(GeneratorNum)%FCPM%SkinLossCurveID, &
FuelCell(GeneratorNum)%FCPM%NdotFuel )
ENDIF
! calculation Step 10, AC FCPM power ancillaries
FuelCell(GeneratorNum)%FCPM%PelancillariesAC = FuelCell(GeneratorNum)%FCPM%ANC0 &
+ FuelCell(GeneratorNum)%FCPM%ANC1 * FuelCell(GeneratorNum)%FCPM%NdotFuel
! calculation Step 11, Dilution air
CAll FigureAirEnthalpy(GeneratorNum, FuelCell(GeneratorNum)%AirSup%TairIntoBlower, Hmolair) ! (Hmolair in KJ/mol)
! units, NdotDilutionAir in kmol/sec.; Hmolair in KJ/mol ,
! factor of 1000's to get to J/s or watts
FuelCell(GeneratorNum)%FCPM%DilutionAirInEnthalpy = Hmolair * 1000.0d0 * FuelCell(GeneratorNum)%FCPM%NdotDilutionAir * 1000.0d0
FuelCell(GeneratorNum)%FCPM%DilutionAirOutEnthalpy = FuelCell(GeneratorNum)%FCPM%DilutionAirInEnthalpy &
+ FuelCell(GeneratorNum)%FCPM%StackHeatLossToDilution
! calculation Step 12, Calculate Reforming water out enthalpy
Call FigureGaseousWaterEnthalpy(FuelCell(GeneratorNum)%FCPM%TprodGasLeavingFCPM, HGasWater)
FuelCell(GeneratorNum)%FCPM%WaterOutEnthalpy = HGasWater * 1000.0d0 * FuelCell(GeneratorNum)%FCPM%NdotLiqwater * 1000.0d0
! calculation Step 13, Calculate Heat balance
! move all terms in Equation 7 to RHS and calculate imbalance
MagofImbalance = - FuelCell(GeneratorNum)%FCPM%TotFuelInEnthalphy &
- FuelCell(GeneratorNum)%FCPM%TotAirInEnthalphy &
- FuelCell(GeneratorNum)%FCPM%WaterInEnthalpy &
- FuelCell(GeneratorNum)%FCPM%DilutionAirInEnthalpy &
- FuelCell(GeneratorNum)%FCPM%NdotFuel * FuelSupply(FuelCell(GeneratorNum)%FuelSupNum)%LHV * 1000000.0d0 &
- FuelCell(GeneratorNum)%FCPM%PelancillariesAC &
+ FuelCell(GeneratorNum)%FCPM%Pel &
+ FuelCell(GeneratorNum)%FCPM%TotProdGasEnthalphy &
+ FuelCell(GeneratorNum)%FCPM%WaterOutEnthalpy &
+ FuelCell(GeneratorNum)%FCPM%QdotStackCool &
+ FuelCell(GeneratorNum)%FCPM%QdotSkin &
+ FuelCell(GeneratorNum)%FCPM%DilutionAirOutEnthalpy
! Now find a new total prod Gas Enthalphy that would result in an energy balance
! TODO check signs...
tmpTotProdGasEnthalphy = FuelCell(GeneratorNum)%FCPM%TotProdGasEnthalphy - MagofImbalance
! solve for a new TprodGasLeavingFCPM using regula falsi method
Acc = 0.01d0 ! guessing need to refine
MaxIter = 150 ! guessing need to refine
SolverFlag = 0 !init
Par(1) = real(GeneratorNum,r64)
Par(2) = tmpTotProdGasEnthalphy
Par(3) = FuelCell(GeneratorNum)%FCPM%NdotProdGas
tmpTprodGas = FuelCell(GeneratorNum)%FCPM%TprodGasLeavingFCPM
CALL SolveRegulaFalsi(Acc, MaxIter, SolverFlag, tmpTprodGas, FuelCellProductGasEnthResidual, &
MinProductGasTemp,MaxProductGasTemp, Par)
IF (SolverFlag == -2) THEN
CALL ShowWarningError('CalcFuelCellGeneratorModel: Regula falsi problem, flag = -2, check signs, all positive')
ENDIF
IF (SolverFlag == -1) THEN
CALL ShowWarningError('CalcFuelCellGeneratorModel: '// &
'Regula falsi problem, flag = -1, check accuracy and iterations, did not converge')
ENDIF
IF (solverFlag > 0) THEN
FuelCell(GeneratorNum)%FCPM%TprodGasLeavingFCPM = tmpTprodGas
! write(*,*) 'Number of regula falsi iterations: ', solverFlag
ENDIF
! moved call to HeatBalanceInternalGains. Call FigureFuelCellZoneGains(GeneratorNum)
! Control Step 3 determine interaction with electrical storage
! How much power is really going into inverter?
PintoInverter = Pel + Pstorage ! Back out so we can reapply
CALL ManageElectStorInteractions(GeneratorNum, Pdemand, PpcuLosses, &
ConstrainedStorage, Pstorage, PgridExtra)
PintoInverter = Pel - Pstorage !
! refine power conditioning losses with more current power production
IF (FuelCell(GeneratorNum)%Inverter%EffMode == InverterEffConstant) THEN
PpcuLosses = (1.0d0 - FuelCell(GeneratorNum)%Inverter%ConstEff )* PintoInverter
ENDIF
IF (FuelCell(GeneratorNum)%Inverter%EffMode == InverterEffQuadratic) THEN
PpcuLosses = (1.0d0 - CurveValue(FuelCell(GeneratorNum)%Inverter%EffQuadraticCurveID, PintoInverter))* PintoInverter
ENDIF
PoutofInverter = PintoInverter - PpcuLosses
FuelCell(GeneratorNum)%ACPowerGen = PoutofInverter - FuelCell(GeneratorNum)%FCPM%PelancillariesAC &
- FuelCell(GeneratorNum)%AirSup%PairCompEl &
- FuelSupply(FuelCell(GeneratorNum)%FuelSupNum)%PfuelCompEl &
- FuelCell(GeneratorNum)%WaterSup%PwaterCompEl
FuelCell(GeneratorNum)%Inverter%PCUlosses = PpcuLosses
! model assumes air intake is drawn over power conditioner to recovery heat
FuelCell(GeneratorNum)%Inverter%QairIntake = FuelCell(GeneratorNum)%Inverter%PCUlosses
Call CalcFuelCellAuxHeater(GeneratorNum)
Call CalcFuelCellGenHeatRecovery(GeneratorNum)
! calculation Step 11, If imbalance below threshold, then exit out of do loop.
IF ( (ABS(MagofImbalance) < ABS(ImBalanceTol * FuelCell(GeneratorNum)%FCPM%Pel)) &
.AND. (Iter > 2) ) THEN
EXIT
ENDIF
ENDDO !sequential substitution loop
FuelCell(GeneratorNum)%FCPM%SeqSubstitIter = iter
FuelCell(GeneratorNum)%FCPM%RegulaFalsiIter = solverFlag
RETURN
END SUBROUTINE CalcFuelCellGeneratorModel