SUBROUTINE ManageElectCenterStorageInteractions(LoadCenterNum,StorageDrawnPower,StorageStoredPower)
! SUBROUTINE INFORMATION:
! AUTHOR B. Griffith
! DATE WRITTEN June-August 2008
! MODIFIED BG May 2009, added EMS
! BN (FSEC) Feb 2010 (pass out two storage values)
! Y. KyungTae & W. Wang July-August, 2011 Added a battery model
! RE-ENGINEERED na
! PURPOSE OF THIS SUBROUTINE:
! manage controls and calculations related to electrical storage in Electric load center
! METHODOLOGY EMPLOYED:
! started from fuel cell storage routine and modified for general use in load center
! needs lots of work
! REFERENCES:
! na
! USE STATEMENTS:
USE General, ONLY: ReallocateRealArray
USE CurveManager, ONLY: CurveValue
USE DataHVACGlobals, ONLY: TimeStepSys,SysTimeElapsed
USE ScheduleManager, ONLY: GetCurrentScheduleValue
USE DataGlobals , ONLY: TimeStep, TimeStepZone, SecInHour, BeginEnvrnFlag, WarmupFlag, HourOfDay
IMPLICIT NONE ! Enforce explicit typing of all variables in this routine
! SUBROUTINE ARGUMENT DEFINITIONS:
INTEGER, INTENT(IN) :: LoadCenterNum ! load center number, index for structure
REAL(r64), INTENT(OUT) :: StorageDrawnPower ! Electric Power Draw Rate from storage units
REAL(r64), INTENT(OUT) :: StorageStoredPower ! Electric Power Store Rate from storage units
! SUBROUTINE PARAMETER DEFINITIONS:
! na
! INTERFACE BLOCK SPECIFICATIONS:
! na
! DERIVED TYPE DEFINITIONS:
! na
! SUBROUTINE LOCAL VARIABLE DECLARATIONS:
REAL(r64) :: tmpPdraw = 0.0D0 ! power draw from storage, working var
REAL(r64) :: tmpPcharge = 0.0D0 ! power charge to storage, working var
LOGICAL :: drawing = .false. ! true if drawing power
LOGICAL :: charging = .false. ! true if charging
INTEGER :: ElecStorNum = 0
REAL(r64) :: Pgensupply = 0.0D0
REAL(r64) :: Pdemand = 0.0D0 !
REAL(r64) :: PpcuLosses = 0.0D0 !
REAL(r64) :: Pstorage = 0.0D0
LOGICAL, ALLOCATABLE,Save, DIMENSION(:) :: MyEnvrnFlag
LOGICAL,SAVE :: MyOneTimeFlag = .true.
LOGICAL, ALLOCATABLE, SAVE, DIMENSION(:) :: MyWarmupFlag ! flag for init after warmup complete
REAL(r64) :: TimeElapsed ! Fraction of the current hour that has elapsed (h)
INTEGER :: BinNum = 0
REAL(r64) :: Input0 = 0.0D0
REAL(r64) I0 ! initial guess of current
REAL(r64) T0 ! initial guess of T(I)
REAL(r64) q0 ! initial charge
REAL(r64) qmax ! maximum capacity
REAL(r64) qmaxf ! maximum capacity, a function of current(I)
REAL(r64) k ! change rate from chemically bound charge to available charge
REAL(r64) c ! fraction available charge capacity to total capacity
REAL(r64) TotalSOC ! total charge (ThisTimeStepAvailable+ThisTimeStepBound)
REAL(r64) Ef ! effective internal voltage source, V
REAL(r64) E0c ! fully charged internal battery voltage
REAL(r64) Volt ! terminal voltage
REAL(r64) Pw ! power required
REAL(r64) E0d ! fully discharged internal battery voltage
REAL(r64) InternalR ! internal resistance
REAL(r64) XF ! normalized maximum capacity at the given current
REAL(r64) X ! normalized maximum capacity at the given current
REAL(r64) Inew ! converged current
REAL(r64) Tnew ! charge of discharge time, defined by T=qmaxf/I
REAL(r64) Imax ! maximum current
REAL(r64) Numpar ! number of battery in parallel
REAL(r64) Numser ! number of battery in series
REAL(r64) Numbattery ! number of battery (Numpar*Numser)
REAL(r64) initialCharge ! initial charge in Ah
REAL(r64) dividend !
REAL(r64) divisor !
REAL(r64) newAvailable ! new available charge in Ah
REAL(r64) newBound ! new bound charge in Ah
REAL(r64) ::error=0.0D0 ! error in iterative process
REAL(r64) ::Pactual=0.0D0 ! actual Power output
REAL(r64) ::RHS=0.0D0 ! right hand side of a equation
REAL(r64) ::I=0.0D0 ! current
REAL(r64) ::DeltaSOC1 ! difference of fractional SOC between this time step and last time step
REAL(r64) ::DeltaSOC2 ! difference of fractional SOC between last time step and last two time step
Integer ::SaveArrayBounds !Maximum size for the arrays used for rainflow counting
If ( .NOT. (ElecLoadCenter(LoadCenterNum)%StoragePresent)) RETURN
ElecStorNum = ElecLoadCenter(LoadCenterNum)%StorageModelNum
!_____________________________________
! begin initalizations
! perform the one time initializations
IF (MyOneTimeFlag) THEN
! initialize the environment and sizing flags
ALLOCATE(MyEnvrnFlag(NumElecStorageDevices))
ALLOCATE(MyWarmupFlag(NumElecStorageDevices))
MyEnvrnFlag = .true.
MyOneTimeFlag = .false.
MyWarmupFlag = .FALSE.
END IF
IF (BeginEnvrnFlag .AND. MyEnvrnFlag(ElecStorNum)) THEN
IF(ElecStorage(ElecStorNum)%StorageModelMode == SimpleBucketStorage) THEN
ElecStorage(ElecStorNum)%LastTimeStepStateOfCharge = ElecStorage(ElecStorNum)%StartingEnergyStored
ElecStorage(ElecStorNum)%ThisTimeStepStateOfCharge = ElecStorage(ElecStorNum)%StartingEnergyStored
ElecStorage(ElecStorNum)%PelNeedFromStorage = 0.d0
ElecStorage(ElecStorNum)%PelFromStorage = 0.d0
ElecStorage(ElecStorNum)%PelIntoStorage = 0.d0
ElecStorage(ElecStorNum)%QdotConvZone = 0.d0
ElecStorage(ElecStorNum)%QdotRadZone = 0.d0
ElecStorage(ElecStorNum)%TimeElapsed = 0.d0
ElecStorage(ElecStorNum)%ElectEnergyinStorage = 0.d0
ElecStorage(ElecStorNum)%StoredPower = 0.d0
ElecStorage(ElecStorNum)%StoredEnergy = 0.d0
ElecStorage(ElecStorNum)%DecrementedEnergyStored = 0.d0
ElecStorage(ElecStorNum)%DrawnPower = 0.d0
ElecStorage(ElecStorNum)%DrawnEnergy = 0.d0
ElecStorage(ElecStorNum)%ThermLossRate = 0.d0
ElecStorage(ElecStorNum)%ThermLossEnergy = 0.d0
ELSEIF(ElecStorage(ElecStorNum)%StorageModelMode == KiBaMBattery) THEN
initialCharge = Elecstorage(ElecStorNum)%MaxAhCapacity * Elecstorage(ElecStorNum)%StartingSOC
Elecstorage(ElecStorNum)%LastTwoTimeStepAvailable = initialCharge * Elecstorage(ElecStorNum)%AvailableFrac
Elecstorage(ElecStorNum)%LastTwoTimeStepBound = initialCharge * (1.0d0-Elecstorage(ElecStorNum)%AvailableFrac)
Elecstorage(ElecStorNum)%LastTimeStepAvailable = initialCharge * Elecstorage(ElecStorNum)%AvailableFrac
Elecstorage(ElecStorNum)%LastTimeStepBound = initialCharge * (1.0d0-Elecstorage(ElecStorNum)%AvailableFrac)
Elecstorage(ElecStorNum)%ThisTimeStepAvailable = initialCharge * Elecstorage(ElecStorNum)%AvailableFrac
Elecstorage(ElecStorNum)%ThisTimeStepBound = initialCharge * (1.0d0-Elecstorage(ElecStorNum)%AvailableFrac)
IF (ElecStorage(ElecStorNum)%LifeCalculation .eq. Battery_LifeCalculation_Yes) THEN
ElecStorage(ElecStorNum)%count0 = 2 ! Index 1 is for initial SOC, so new input starts from index 2.
ElecStorage(ElecStorNum)%B10(1) = Elecstorage(ElecStorNum)%StartingSOC ! the initial fractional SOC is stored as the reference
ElecStorage(ElecStorNum)%X0 = 0.0D0
ElecStorage(ElecStorNum)%OneNmb0 = 0.0D0
ElecStorage(ElecStorNum)%Nmb0 = 0.0D0
ElecStorage(ElecStorNum)%BatteryDamage = 0.0D0
ENDIF
ENDIF
MyEnvrnFlag(ElecStorNum) = .FALSE.
MyWarmupFlag(ElecStorNum) = .TRUE.
ENDIF
IF (.NOT. BeginEnvrnFlag) MyEnvrnFlag(ElecStorNum) = .TRUE.
IF (MyWarmupFlag(ElecStorNum) .AND. (.NOT. WarmUpFlag) ) THEN
! need to reset initial state of charge at beginning of environment but after warm up is complete
IF(ElecStorage(ElecStorNum)%StorageModelMode == SimpleBucketStorage) THEN
ElecStorage(ElecStorNum)%LastTimeStepStateOfCharge = ElecStorage(ElecStorNum)%StartingEnergyStored
ElecStorage(ElecStorNum)%ThisTimeStepStateOfCharge = ElecStorage(ElecStorNum)%StartingEnergyStored
ELSEIF(ElecStorage(ElecStorNum)%StorageModelMode == KiBaMBattery) THEN
initialCharge = Elecstorage(ElecStorNum)%MaxAhCapacity * Elecstorage(ElecStorNum)%StartingSOC
Elecstorage(ElecStorNum)%LastTwoTimeStepAvailable = initialCharge * Elecstorage(ElecStorNum)%AvailableFrac
Elecstorage(ElecStorNum)%LastTwoTimeStepBound = initialCharge * (1.0d0-Elecstorage(ElecStorNum)%AvailableFrac)
Elecstorage(ElecStorNum)%LastTimeStepAvailable = initialCharge * Elecstorage(ElecStorNum)%AvailableFrac
Elecstorage(ElecStorNum)%LastTimeStepBound = initialCharge * (1.0d0-Elecstorage(ElecStorNum)%AvailableFrac)
Elecstorage(ElecStorNum)%ThisTimeStepAvailable = initialCharge * Elecstorage(ElecStorNum)%AvailableFrac
Elecstorage(ElecStorNum)%ThisTimeStepBound = initialCharge * (1.0d0-Elecstorage(ElecStorNum)%AvailableFrac)
IF (ElecStorage(ElecStorNum)%LifeCalculation .eq. Battery_LifeCalculation_Yes) THEN
ElecStorage(ElecStorNum)%count0 = 2 ! Index 1 is for initial SOC, so new input starts from index 2.
ElecStorage(ElecStorNum)%B10(1) = Elecstorage(ElecStorNum)%StartingSOC ! the initial fractional SOC is stored as the reference
ElecStorage(ElecStorNum)%X0 = 0.0D0
ElecStorage(ElecStorNum)%OneNmb0 = 0.0D0
ElecStorage(ElecStorNum)%Nmb0 = 0.0D0
ElecStorage(ElecStorNum)%BatteryDamage = 0.0D0
ENDIF
ENDIF
MyWarmupFlag(ElecStorNum) = .FALSE.
ENDIF
TimeElapsed = HourOfDay + TimeStep * TimeStepZone + SysTimeElapsed
If (ElecStorage(ElecStorNum)%TimeElapsed /= TimeElapsed) Then !time changed, update last with "current" result from previous time
! When program comes here, this means the first iteration of a new system time step. Battery life calculation needs to be considered
! for the Kinetic battery model.
IF (ElecStorage(ElecStorNum)%StorageModelMode == KiBaMBattery .AND. &
ElecStorage(ElecStorNum)%LifeCalculation.eq.Battery_LifeCalculation_Yes) THEN
! At this point, the current values, last time step values and last two time step values have not been updated, hence:
! "ThisTimeStep*" actually points to the previous one time step
! "LastTimeStep*" actually points to the previous two time steps
! "LastTwoTimeStep" actually points to the previous three time steps
! Calculate the fractional SOC change between the "current" time step and the "previous one" time step
DeltaSOC1 = ElecStorage(ElecStorNum)%ThisTimeStepAvailable + ElecStorage(ElecStorNum)%ThisTimeStepBound &
- ElecStorage(ElecStorNum)%LastTimeStepAvailable - ElecStorage(ElecStorNum)%LastTimeStepBound
DeltaSOC1 = DeltaSOC1/Elecstorage(ElecStorNum)%MaxAhCapacity
! Calculate the fractional SOC change between the "previous one" time step and the "previous two" time steps
DeltaSOC2 = ElecStorage(ElecStorNum)%LastTimeStepAvailable + ElecStorage(ElecStorNum)%LastTimeStepBound &
- ElecStorage(ElecStorNum)%LastTwoTimeStepAvailable - ElecStorage(ElecStorNum)%LastTwoTimeStepBound
DeltaSOC2 = DeltaSOC2/Elecstorage(ElecStorNum)%MaxAhCapacity
! DeltaSOC2 = 0 may occur at the begining of each simulation environment.
! DeltaSOC1 * DeltaSOC2 means that the SOC from "LastTimeStep" is a peak or valley. Only peak or valley needs
! to call the rain flow algorithm
IF ( (DeltaSOC2 == 0) .OR. ((DeltaSOC1 * DeltaSOC2) .lt. 0) ) THEN
! Because we cannot determine whehter "ThisTimeStep" is a peak or valley (next time step is unknown yet), we
! use the "LastTimeStep" value for battery life calculation.
Input0 = (ElecStorage(ElecStorNum)%LastTimeStepAvailable + ElecStorage(ElecStorNum)%LastTimeStepBound)/ &
Elecstorage(ElecStorNum)%MaxAhCapacity
ElecStorage(ElecStorNum)%B10(ElecStorage(ElecStorNum)%count0) = input0
! The arrary size needs to be increased when count = MAXRainflowArrayBounds. Please note that (MAXRainflowArrayBounds +1)
! is the index used in the subroutine RainFlow. So we cannot reallocate array size until count = MAXRainflowArrayBounds +1.
IF (ElecStorage(ElecStorNum)%count0 == MAXRainflowArrayBounds) THEN
SaveArrayBounds=MAXRainflowArrayBounds+1
CALL ReallocateRealArray(ElecStorage(ElecStorNum)%B10,SaveArrayBounds,MAXRainFlowArrayInc)
SaveArrayBounds=MAXRainflowArrayBounds+1
CALL ReallocateRealArray(ElecStorage(ElecStorNum)%X0,SaveArrayBounds,MAXRainFlowArrayInc)
! Decrement by 1 is needed because the last input parameter of RainFlow is MAXRainflowArrayBounds+1
MAXRainflowArrayBounds=SaveArrayBounds-1
ENDIF
Call Rainflow(ElecStorage(ElecStorNum)%CycleBinNum,input0,ElecStorage(ElecStorNum)%B10, &
ElecStorage(ElecStorNum)%X0,ElecStorage(ElecStorNum)%count0, &
ElecStorage(ElecStorNum)%Nmb0,ElecStorage(ElecStorNum)%OneNmb0,MAXRainflowArrayBounds+1)
ElecStorage(ElecStorNum)%BatteryDamage=0.0d0
DO BinNum = 1, ElecStorage(ElecStorNum)%CycleBinNum
! Battery damage is calculated by accumulating the impact from each cycle.
ElecStorage(ElecStorNum)%BatteryDamage = ElecStorage(ElecStorNum)%OneNmb0(BinNum)/ &
Curvevalue(ElecStorage(ElecStorNum)%LifeCurveNum, &
(Real(BinNum,r64)/Real(ElecStorage(ElecStorNum)%CycleBinNum,r64))) &
+ ElecStorage(ElecStorNum)%BatteryDamage
ENDDO
ENDIF
ENDIF
! Updating "LastTimeStep" and "LastTwoTimeStep" should be done after calling the rain flow algorithm.
! Otherwise, it is not possible to determine whehter "LastTimeStep" is a peak or valley value.
ElecStorage(ElecStorNum)%LastTimeStepStateOfCharge = ElecStorage(ElecStorNum)%ThisTimeStepStateOfCharge
Elecstorage(ElecStorNum)%LastTwoTimeStepAvailable = Elecstorage(ElecStorNum)%LastTimeStepAvailable
Elecstorage(ElecStorNum)%LastTwoTimeStepBound = Elecstorage(ElecStorNum)%LastTimeStepBound
Elecstorage(ElecStorNum)%LastTimeStepAvailable=Elecstorage(ElecStorNum)%ThisTimeStepAvailable
Elecstorage(ElecStorNum)%LastTimeStepBound=Elecstorage(ElecStorNum)%ThisTimeStepBound
ElecStorage(ElecStorNum)%TimeElapsed = TimeElapsed
ENDIF
! end Inits
!Begin determine Power demand request
If (ElecLoadCenter(LoadCenterNum)%InverterPresent) Then
! this will be lagged from last calc.
PpcuLosses = Inverter(ElecLoadCenter(LoadCenterNum)%InverterModelNum)%ThermLossRate
ELSE
PpcuLosses = 0.0D0
ENDIF
! Pdemand = Sum (ElecLoadCenter(LoadCenterNum)%ElecGen%PowerRequestThisTimestep) + PpcuLosses
! Modified when improved electric load center disptach procedure was implemented
Pdemand = ElecLoadCenter(LoadCenterNum)%TotalPowerRequest + PpcuLosses
! END determine Power demand request.
! BEGIN determine available generation supply
SELECT CASE (ElecLoadCenter(LoadCenterNum)%BussType)
CASE (ACBussStorage)
Pgensupply = Sum(ElecLoadCenter(LoadCenterNum)%ElecGen%ElectProdRate)
CASE (DCBussInverterDCStorage)
ElecLoadCenter(LoadCenterNum)%DCElectProdRate = Sum(ElecLoadCenter(LoadCenterNum)%ElecGen%DCElectProdRate)
Pgensupply = ElecLoadCenter(LoadCenterNum)%DCElectProdRate
CASE (DCBussInverterACStorage)
Pgensupply = Inverter(ElecLoadCenter(LoadCenterNum)%InverterModelNum)%ACPowerOut
END SELECT
! End determine available generation
!initialize locals
tmpPdraw = 0.d0
tmpPcharge = 0.d0
drawing = .false.
charging = .false.
Pstorage = 0.d0
IF (ElecStorage(ElecStorNum)%StorageModelMode == KiBaMBattery) THEN
Numpar = ElecStorage(ElecStorNum)%ParallelNum
Numser = ElecStorage(ElecStorNum)%SeriesNum
Numbattery = Numpar*Numser
InternalR = Elecstorage(ElecStorNum)%InternalR
qmax = Elecstorage(ElecStorNum)%MaxAhCapacity
E0c = ElecStorage(ElecStorNum)%ChargedOCV
E0d = ElecStorage(ElecStorNum)%DischargedOCV
k = ElecStorage(ElecStorNum)%ChargeConversionRate
c = ElecStorage(ElecStorNum)%AvailableFrac
ENDIF
!*****************************************************************************************************************
! step 1 figure out what is desired of electrical storage system
If (Pgensupply < (Pdemand )) THEN
!draw from storage
tmpPdraw = Pdemand - Pgensupply
tmpPcharge = 0.d0
drawing = .true.
charging = .false.
ELSEIF (Pgensupply > (Pdemand )) THEN
!add to storage
tmpPcharge = Pgensupply - Pdemand
tmpPdraw = 0.d0
charging = .true.
drawing = .false.
ELSEIF (Pgensupply == (Pdemand )) THEN
!do nothing
tmpPcharge = 0.d0
tmpPdraw = 0.d0
charging = .false.
drawing = .false.
ENDIF
! EMS override -- intent to draw, charge or hold?
IF ((ElecStorage(ElecStorNum)%EMSOverridePelFromStorage) &
.AND. (.NOT. ElecStorage(ElecStorNum)%EMSOverridePelIntoStorage) ) THEN
! EMS is calling for specific discharge rate
tmpPdraw = MAX(ElecStorage(ElecStorNum)%EMSValuePelFromStorage, 0.0D0)
tmpPcharge = 0.0D0
drawing = .TRUE.
charging = .FALSE.
ELSEIF((.NOT. ElecStorage(ElecStorNum)%EMSOverridePelFromStorage) &
.AND. (ElecStorage(ElecStorNum)%EMSOverridePelIntoStorage)) THEN
! EMS is calling for specific charge rate
tmpPcharge= MAX(ElecStorage(ElecStorNum)%EMSValuePelIntoStorage, 0.0D0)
tmpPdraw = 0.0D0
drawing = .FALSE.
charging = .TRUE.
ELSEIF ((ElecStorage(ElecStorNum)%EMSOverridePelFromStorage) &
.AND. (ElecStorage(ElecStorNum)%EMSOverridePelIntoStorage) ) THEN
! EMS is asking to override both
IF ( ElecStorage(ElecStorNum)%EMSValuePelIntoStorage > ElecStorage(ElecStorNum)%EMSValuePelFromStorage) THEN
tmpPcharge= ElecStorage(ElecStorNum)%EMSValuePelIntoStorage - ElecStorage(ElecStorNum)%EMSValuePelFromStorage
tmpPdraw = 0.0D0
drawing = .FALSE.
charging = .TRUE.
ELSEIF ( ElecStorage(ElecStorNum)%EMSValuePelIntoStorage < ElecStorage(ElecStorNum)%EMSValuePelFromStorage) THEN
tmpPdraw = ElecStorage(ElecStorNum)%EMSValuePelFromStorage - ElecStorage(ElecStorNum)%EMSValuePelIntoStorage
tmpPcharge = 0.0D0
drawing = .TRUE.
charging = .FALSE.
ELSE !they equal just hold
tmpPdraw = 0.0D0
tmpPcharge = 0.0D0
drawing = .FALSE.
charging = .FALSE.
ENDIF
ENDIF
!*****************************************************************************************************************
! step 2, figure out what is possible for electrical storage draws/charges
If (charging) then
IF (ElecStorage(ElecStorNum)%StorageModelMode == SimpleBucketStorage) THEN
IF (ElecStorage(ElecStorNum)%LastTimeStepStateOfCharge >= ElecStorage(ElecStorNum)%MaxEnergyCapacity) THEN
! storage full! no more allowed!
tmpPcharge = 0.d0
! Constrained = .true.
charging = .FALSE.
ENDIF
IF (tmpPcharge > ElecStorage(ElecStorNum)%MaxPowerStore) THEN
tmpPcharge = ElecStorage(ElecStorNum)%MaxPowerStore
! Constrained = .true.
ENDIF
!now add energy to storage from charging
IF ((ElecStorage(ElecStorNum)%LastTimeStepStateOfCharge &
+ tmpPcharge *TimeStepSys*SecInHour*ElecStorage(ElecStorNum)%EnergeticEfficCharge) &
< ElecStorage(ElecStorNum)%MaxEnergyCapacity) THEN
ElecStorage(ElecStorNum)%ThisTimeStepStateOfCharge = ElecStorage(ElecStorNum)%LastTimeStepStateOfCharge &
+ tmpPcharge *TimeStepSys*SecInHour*ElecStorage(ElecStorNum)%EnergeticEfficCharge
ELSE ! would over charge this time step
tmpPcharge = (ElecStorage(ElecStorNum)%MaxEnergyCapacity - ElecStorage(ElecStorNum)%LastTimeStepStateOfCharge) &
/(TimeStepSys*SecInHour*ElecStorage(ElecStorNum)%EnergeticEfficCharge)
! Constrained = .true.
ElecStorage(ElecStorNum)%ThisTimeStepStateOfCharge = ElecStorage(ElecStorNum)%LastTimeStepStateOfCharge &
+ tmpPcharge *TimeStepSys*SecInHour*ElecStorage(ElecStorNum)%EnergeticEfficCharge
ENDIF
Pstorage = tmpPcharge
ENDIF
IF (ElecStorage(ElecStorNum)%StorageModelMode == KiBaMBattery) THEN
!*************************************************
!The sign of power and current is negative in charging
!*************************************************
Pw = -tmpPcharge/Numbattery
q0 = ElecStorage(ElecStorNum)%LastTimeStepAvailable + ElecStorage(ElecStorNum)%LastTimeStepBound
I0 = 1.0d0 ! Initial assumption
T0 = abs(qmax/I0) ! Initial Assumption
qmaxf = qmax*k*c*T0/(1.0d0-exp(-k*T0)+c*(k*T0-1.0d0+exp(-k*T0))) !Initial calculation of a function qmax(I)
Xf = q0/qmaxf
Ef = E0d + CurveValue(ElecStorage(ElecStorNum)%ChargeCurveNum,Xf) !E0d+Ac*Xf+Cc*Xf/(Dc-Xf) (use curve)
Volt = Ef-I0*InternalR
Inew = Pw/Volt
Tnew = qmaxf/abs(Inew)
error = 1.0d0
DO WHILE (error.gt.0.0001d0) !Iteration process to get converged current(I)
I0 = Inew
T0 = Tnew
qmaxf = qmax*k*c*T0/(1.0d0-exp(-k*T0)+c*(k*T0-1.0d0+exp(-k*T0)))
Xf = q0/qmaxf
Ef = E0d+CurveValue(ElecStorage(ElecStorNum)%ChargeCurveNum,Xf) !E0d+Ac*Xf+Cc*Xf/(Dc-Xf) (use curve)
Volt = Ef-I0*InternalR
Inew = Pw/Volt
Tnew = abs(qmaxf/Inew) ! ***Always positive here
error = abs(Inew-I0)
ENDDO
dividend = -k*c*qmax + k*ElecStorage(ElecStorNum)%LastTimeStepAvailable*exp(-k*TimeStepSys) + &
q0*k*c*(1.0d0-exp(-k*TimeStepSys))
divisor = 1.0d0 - exp(-k*TimeStepSys) + c*(k*TimeStepSys-1+exp(-k*TimeStepSys))
Imax = dividend/divisor
! Below: This is the limit of charging current from Charge Rate Limit (input)
Imax = Max(Imax,-(qmax-q0)*ElecStorage(ElecStorNum)%MaxChargeRate)
IF (abs(I0).le.abs(Imax)) THEN
I0 = Pw/Volt
Pactual = I0*Volt
ELSE
I0 = Imax
qmaxf = 80.0d0 !Initial assumption to solve the equation using iterative method
error = 10.0d0 !Initial assumption ...
DO WHILE (error.gt.0.001d0)
! *** I0(current) should be positive for this calculation
RHS=(qmax*k*c*qmaxf/abs(I0))/(1.0d0-exp(-k*qmaxf/abs(I0))+c*(k*qmaxf/abs(I0)-1.0d0+exp(-k*qmaxf/abs(I0))))
error=abs(qmaxf-RHS)
qmaxf=RHS
ENDDO
ENDIF
ENDIF ! end KiBaM charging
ENDIF !charging
IF (drawing) THEN
IF (ElecStorage(ElecStorNum)%StorageModelMode == SimpleBucketStorage) THEN
IF (ElecStorage(ElecStorNum)%LastTimeStepStateOfCharge <= 0.d0) THEN
! storage empty no more allowed!
tmpPdraw = 0.0d0
drawing = .FALSE.
ENDIF
IF (tmpPdraw > ElecStorage(ElecStorNum)%MaxPowerDraw) THEN
tmpPdraw = ElecStorage(ElecStorNum)%MaxPowerDraw
ENDIF
!now take energy from storage by drawing (amplified by energetic effic)
IF ((ElecStorage(ElecStorNum)%LastTimeStepStateOfCharge &
- tmpPdraw *TimeStepSys*SecInHour/ElecStorage(ElecStorNum)%EnergeticEfficDischarge) > 0.0d0 ) Then
ElecStorage(ElecStorNum)%ThisTimeStepStateOfCharge = ElecStorage(ElecStorNum)%LastTimeStepStateOfCharge &
- tmpPdraw *TimeStepSys*SecInHour/ElecStorage(ElecStorNum)%EnergeticEfficDischarge
ELSE !would over drain storage this timestep so reduce tmpPdraw
tmpPdraw = ElecStorage(ElecStorNum)%LastTimeStepStateOfCharge * ElecStorage(ElecStorNum)%EnergeticEfficDischarge &
/(TimeStepSys*SecInHour)
ElecStorage(ElecStorNum)%ThisTimeStepStateOfCharge = ElecStorage(ElecStorNum)%LastTimeStepStateOfCharge &
- tmpPdraw *TimeStepSys*SecInHour/ElecStorage(ElecStorNum)%EnergeticEfficDischarge
ElecStorage(ElecStorNum)%ThisTimeStepStateOfCharge = MAX(ElecStorage(ElecStorNum)%ThisTimeStepStateOfCharge, 0.0D0)
ENDIF
ElecStorage(ElecStorNum)%ThermLossRate = tmpPdraw * (1.0d0/ElecStorage(ElecStorNum)%EnergeticEfficDischarge - 1.0d0)
Pstorage = - tmpPdraw
ENDIF ! simple discharging
IF (ElecStorage(ElecStorNum)%StorageModelMode == KiBaMBattery) THEN
!**********************************************
!The sign of power and current is positive in discharging
!**********************************************
Pw = tmpPdraw/Numbattery
q0 = ElecStorage(ElecStorNum)%LastTimeStepAvailable+ElecStorage(ElecStorNum)%LastTimeStepBound
I0 = 10.0d0 ! Initial assumption
T0 = qmax/I0 ! Initial Assumption
qmaxf = qmax*k*c*T0/(1.0d0-exp(-k*T0)+c*(k*T0-1.0d0+exp(-k*T0))) !Initial calculation of a function qmax(I)
Xf = (qmax-q0)/qmaxf
Ef = E0c + CurveValue(ElecStorage(ElecStorNum)%DischargeCurveNum,Xf) !E0d+Ac*Xf+Cc*X/(Dc-Xf)
Volt = Ef-I0*InternalR
Inew = Pw/Volt
Tnew = qmaxf/Inew
error = 1.0d0
DO WHILE (error.gt.0.0001d0) !Iteration process to get converged current(I)
I0 = Inew
T0 = Tnew
qmaxf= qmax*k*c*T0/(1.0d0-exp(-k*T0)+c*(k*T0-1.0d0+exp(-k*T0)))
Xf = (qmax-q0)/qmaxf
Ef = E0c + CurveValue(ElecStorage(ElecStorNum)%DischargeCurveNum,Xf) !E0c+Ad*Xf+Cd*X/(Dd-Xf)
Volt = Ef - I0*InternalR
Inew = Pw/Volt
Tnew = qmaxf/Inew
error= abs(Inew-I0)
ENDDO
dividend = k*ElecStorage(ElecStorNum)%LastTimeStepAvailable*exp(-k*TimeStepSys) + &
q0*k*c*(1.0d0-exp(-k*TimeStepSys))
divisor = 1.0d0 - exp(-k*TimeStepSys) + c*(k*TimeStepSys-1.0d0+exp(-k*TimeStepSys))
Imax = dividend/divisor
Imax = Min(Imax,ElecStorage(ElecStorNum)%MaxDischargeI)
IF (abs(I0).le.Imax) THEN
I0 = Pw/Volt
Pactual = I0*Volt
ELSE
I0 = Imax
qmaxf = 10.0d0 !Initial assumption to solve the equation using iterative method
error = 10.0d0 !Initial assumption ...
DO WHILE (error.gt.0.001d0)
RHS = (qmax*k*c*qmaxf/I0)/(1.0d0-exp(-k*qmaxf/I0)+c*(k*qmaxf/I0-1+exp(-k*qmaxf/I0)))
error = abs(qmaxf-RHS)
qmaxf = RHS
ENDDO
Xf = (qmax-q0)/qmaxf
Ef = E0c+CurveValue(ElecStorage(ElecStorNum)%DischargeCurveNum,Xf)
Volt = Ef-I0*InternalR
ENDIF
IF (Volt.lt.ElecStorage(ElecStorNum)%CutoffV) THEN
I0 = 0.d0
ENDIF
ENDIF ! end KiBaM discharging
ENDIF !drawing
! correct if not available.
IF (GetCurrentScheduleValue(ElecStorage(ElecStorNum)%AvailSchedPtr) == 0.d0) THEN
IF ((.NOT. ElecStorage(ElecStorNum)%EMSOverridePelFromStorage) &
.AND. (.NOT. ElecStorage(ElecStorNum)%EMSOverridePelIntoStorage)) THEN
charging = .FALSE.
drawing = .FALSE.
ENDIF
ENDIF
IF (ElecStorage(ElecStorNum)%StorageModelMode == SimpleBucketStorage) THEN
IF ((.NOT. charging) .AND. ( .NOT. drawing)) THEN
ElecStorage(ElecStorNum)%ThisTimeStepStateOfCharge = ElecStorage(ElecStorNum)%LastTimeStepStateOfCharge
ElecStorage(ElecStorNum)%PelIntoStorage = 0.d0
ElecStorage(ElecStorNum)%PelFromStorage = 0.d0
Pstorage = 0.d0
ENDIF
IF (Pstorage >= 0.d0) THEN
ElecStorage(ElecStorNum)%PelIntoStorage = Pstorage
ElecStorage(ElecStorNum)%PelFromStorage = 0.0d0
ENDIF
IF (Pstorage < 0.d0) THEN
ElecStorage(ElecStorNum)%PelIntoStorage = 0.d0
ElecStorage(ElecStorNum)%PelFromStorage = - Pstorage
ENDIF
ElecStorage(ElecStorNum)%ElectEnergyinStorage = ElecStorage(ElecStorNum)%ThisTimeStepStateOfCharge
ElecStorage(ElecStorNum)%StoredPower = ElecStorage(ElecStorNum)%PelIntoStorage
ElecStorage(ElecStorNum)%StoredEnergy = ElecStorage(ElecStorNum)%PelIntoStorage * TimeStepSys * SecInHour
ElecStorage(ElecStorNum)%DecrementedEnergyStored = -1.0D0 * ElecStorage(ElecStorNum)%StoredEnergy
ElecStorage(ElecStorNum)%DrawnPower = ElecStorage(ElecStorNum)%PelFromStorage
ElecStorage(ElecStorNum)%DrawnEnergy = ElecStorage(ElecStorNum)%PelFromStorage * TimeStepSys * SecInHour
ElecStorage(ElecStorNum)%ThermLossRate = MAX(ElecStorage(ElecStorNum)%StoredPower &
* (1.0D0 - ElecStorage(ElecStorNum)%EnergeticEfficCharge ) , &
ElecStorage(ElecStorNum)%DrawnPower &
* (1.0D0 - ElecStorage(ElecStorNum)%EnergeticEfficDischarge ) )
ElecStorage(ElecStorNum)%ThermLossEnergy = ElecStorage(ElecStorNum)%ThermLossRate * TimeStepSys * SecInHour
IF (ElecStorage(ElecStorNum)%ZoneNum > 0) THEN ! set values for zone heat gains
ElecStorage(ElecStorNum)%QdotconvZone = (1.0D0 - ElecStorage(ElecStorNum)%ZoneRadFract)* ElecStorage(ElecStorNum)%ThermLossRate
ElecStorage(ElecStorNum)%QdotRadZone = (ElecStorage(ElecStorNum)%ZoneRadFract) * ElecStorage(ElecStorNum)%ThermLossRate
ENDIF
StorageStoredPower=ElecStorage(ElecStorNum)%StoredPower
StorageDrawnPower=ElecStorage(ElecStorNum)%DrawnPower
ENDIF ! Outputs for simple battery model
IF (ElecStorage(ElecStorNum)%StorageModelMode == KiBaMBattery) THEN
IF ((.NOT. charging) .AND. ( .NOT. drawing)) THEN
ElecStorage(ElecStorNum)%ThisTimeStepAvailable = ElecStorage(ElecStorNum)%LastTimeStepAvailable
ElecStorage(ElecStorNum)%ThisTimeStepBound = ElecStorage(ElecStorNum)%LastTimeStepBound
I0 = 0.d0
Volt = 0.d0
q0 = ElecStorage(ElecStorNum)%LastTimeStepAvailable+ElecStorage(ElecStorNum)%LastTimeStepBound
ELSE
newAvailable = ElecStorage(ElecStorNum)%LastTimeStepAvailable*exp(-k*TimeStepSys) + &
(q0*k*c-I0)*(1.0d0-exp(-k*TimeStepSys))/k - I0*c*(k*TimeStepSys-1.0d0+exp(-k*TimeStepSys))/k
newBound = ElecStorage(ElecStorNum)%LastTimeStepBound*exp(-k*TimeStepSys)+q0*(1.0d0-c)* &
(1.0d0-exp(-k*TimeStepSys))-I0*(1.0d0-c)*(k*TimeStepSys-1.0d0+exp(-k*TimeStepSys))/k
ElecStorage(ElecStorNum)%ThisTimeStepAvailable = Max(0.0d0 , newAvailable)
ElecStorage(ElecStorNum)%ThisTimeStepBound = Max(0.0d0 , newBound)
ENDIF
Pactual = I0*Volt
TotalSOC = ElecStorage(ElecStorNum)%ThisTimeStepAvailable + ElecStorage(ElecStorNum)%ThisTimeStepBound
!output1
IF (TotalSOC.gt.q0) THEN
ElecStorage(ElecStorNum)%StorageMode = 2
ElecStorage(ElecStorNum)%StoredPower = Volt*I0*Numbattery
ElecStorage(ElecStorNum)%StoredEnergy = Volt*I0*Numbattery*TimeStepSys*SecInHour
ElecStorage(ElecStorNum)%DecrementedEnergyStored = -1.0D0 * ElecStorage(ElecStorNum)%StoredEnergy
ElecStorage(ElecStorNum)%DrawnPower = 0.0d0
ElecStorage(ElecStorNum)%DrawnEnergy = 0.0d0
ELSEIF (TotalSOC.lt.q0) THEN
ElecStorage(ElecStorNum)%StorageMode = 1
ElecStorage(ElecStorNum)%StoredPower = 0.0d0
ElecStorage(ElecStorNum)%StoredEnergy = 0.0d0
ElecStorage(ElecStorNum)%DecrementedEnergyStored = 0.0d0
ElecStorage(ElecStorNum)%Drawnpower = Volt*I0*Numbattery
ElecStorage(ElecStorNum)%Drawnenergy = Volt*I0*Numbattery*TimeStepSys*SecInHour
ELSE
ElecStorage(ElecStorNum)%StorageMode = 0
ElecStorage(ElecStorNum)%StoredPower = 0.0d0
ElecStorage(ElecStorNum)%StoredEnergy = 0.0d0
ElecStorage(ElecStorNum)%DecrementedEnergyStored = 0.0d0
ElecStorage(ElecStorNum)%DrawnPower = 0.0d0
ElecStorage(ElecStorNum)%DrawnEnergy = 0.0d0
ENDIF
ElecStorage(ElecStorNum)%AbsoluteSOC = TotalSOC*Numbattery
ElecStorage(ElecStorNum)%FractionSOC = TotalSOC/qmax
ElecStorage(ElecStorNum)%BatteryCurrent = I0*Numpar
ElecStorage(ElecStorNum)%BatteryVoltage = Volt*Numser
ElecStorage(ElecStorNum)%ThermLossRate = InternalR*I0**2*Numbattery
ElecStorage(ElecStorNum)%ThermLossEnergy = InternalR*I0**2*TimeStepSys*SecInHour*Numbattery
IF (ElecStorage(ElecStorNum)%ZoneNum > 0) THEN ! set values for zone heat gains
ElecStorage(ElecStorNum)%QdotconvZone = ( (1.0D0 - ElecStorage(ElecStorNum)%ZoneRadFract) * &
ElecStorage(ElecStorNum)%ThermLossRate) * Numbattery
ElecStorage(ElecStorNum)%QdotRadZone = ((ElecStorage(ElecStorNum)%ZoneRadFract) * &
ElecStorage(ElecStorNum)%ThermLossRate)*Numbattery
ENDIF
StorageStoredPower=ElecStorage(ElecStorNum)%StoredPower
StorageDrawnPower=ElecStorage(ElecStorNum)%DrawnPower
ENDIF ! Outputs for kibam model
RETURN
END SUBROUTINE ManageElectCenterStorageInteractions