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.
Type | Intent | Optional | Attributes | Name | ||
---|---|---|---|---|---|---|
character(len=*), | intent(in) | :: | GshpType | |||
character(len=*), | intent(in) | :: | GshpName | |||
integer, | intent(in) | :: | GSHPNum | |||
real(kind=r64) | :: | 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 CalcGshpModel( GSHPType, GSHPName, GSHPNum, MyLoad,FirstHVACIteration )
! SUBROUTINE INFORMATION:
! AUTHOR
! DATE WRITTEN Sept. 1998
! MODIFIED April 1999
! September 2002, SJR
! RE-ENGINEERED Mar2000
! PURPOSE OF THIS SUBROUTINE: This routine performs
! METHODOLOGY EMPLOYED: under development
! REFERENCES:
! USE STATEMENTS:
USE DataHVACGlobals, ONLY : TimeStepSys ,SysTimeElapsed, FirstTimeStepSysFlag
USE FluidProperties
USE General, ONLY: TrimSigDigits
USE DataPlant, ONLY: PlantLoop
USE DataBranchAirLoopPlant, ONLY: MassFlowTolerance
USE PlantUtilities, ONLY: SetComponentFlowRate
IMPLICIT NONE
! SUBROUTINE ARGUMENT DEFINITIONS:
CHARACTER(len=*), INTENT(IN) :: GshpType ! type ofGSHP
CHARACTER(len=*), INTENT(IN) :: GshpName ! user specified name ofGSHP
INTEGER , INTENT(IN) :: GSHPNum ! GSHP Number
REAL(r64) :: MyLoad ! Operating Load
LOGICAL, INTENT(IN) :: FirstHVACIteration
! SUBROUTINE PARAMETER DEFINITIONS:
REAL(r64), PARAMETER :: gamma = 1.114d0 ! Expnasion Coefficient
REAL(r64), PARAMETER :: HeatBalTol = 0.0005d0
REAL(r64), PARAMETER :: RelaxParam = 0.6d0
REAL(r64), PARAMETER :: SmallNum = 1.0d-20
INTEGER, PARAMETER :: IterationLimit = 500
! INTERFACE BLOCK SPECIFICATIONS
! na
! DERIVED TYPE DEFINITIONS
! na
! SUBROUTINE LOCAL VARIABLE DECLARATIONS:
REAL(r64) :: SourceSideEffect ! Source Side effectiveness
REAL(r64) :: LoadSideEffect ! Load Side effectiveness
REAL(r64) :: SourceSideTemp ! Source Side temperature °C
REAL(r64) :: LoadSideTemp ! Load Side temperature °C
REAL(r64) :: SourceSideUA ! Source Side heat transfer coefficient w/k
REAL(r64) :: LoadSideUA ! Load Side heat transfer coefficient W/k
REAL(r64) :: SourceSidePressure ! Source Side pressure Pascals
REAL(r64) :: LoadSidePressure ! Load Side pressure Pascals
REAL(r64) :: SuctionPr ! Suction Pressure pascals
REAL(r64) :: DischargePr ! Discharge Pressure pascals
REAL(r64) :: CompressInletTemp ! Compressor inlet temperature °C
REAL(r64) :: PressureDrop ! Suction Pressure drop °C
REAL(r64) :: ClearanceFactor ! Clearance factor
REAL(r64) :: PistonDisp ! Compressor piston displacement m3
REAL(r64) :: ShTemp ! Superheat temperature °C
REAL(r64) :: LosFac ! Loss factor used to define the electromechanical loss for compressor
REAL(r64) :: MassRef ! mass flow rate of refrigerant Kg/s
REAL(r64) :: SourceSideOutletEnth ! Enthalpy at Source Side pressure Joules
REAL(r64) :: LoadSideOutletEnth ! Enthalpy at Condensor Pressure Joules
REAL(r64) :: initialQSource ! Guess Source Side Heat rate Joules
REAL(r64) :: initialQLoad ! Guess Load Side Heat rate Joules
REAL(r64) :: qual ! quality
REAL(r64) :: SuperHeatEnth
REAL(r64) :: T110
REAL(r64) :: T111
REAL(r64) :: CompSuctionTemp
REAL(r64) :: CompSuctionEnth
REAL(r64) :: CompSuctionDensity
REAL(r64) :: PowerLosses
REAL(r64) :: CompSuctionSatTemp
REAL(r64) :: HighPressCutOff
REAL(r64) :: LowPressCutOff
CHARACTER(len=25) :: ErrString
REAL(r64) :: DutyFactor
INTEGER :: IterationCount
REAL(r64), SAVE :: CurrentSimTime = 0.0d0
REAL(r64), SAVE :: PrevSimTime = 0.0d0
LOGICAL, SAVE :: OneTimeFlag = .TRUE.
! Nodes
INTEGER :: SourceSideInletNode ! Source Side inlet node number, water side
INTEGER :: SourceSideOutletNode ! Source Side outlet node number, water side
INTEGER :: LoadSideInletNode ! Load Side inlet node number, water side
INTEGER :: LoadSideOutletNode ! Load Side outlet node number, water side
INTEGER :: LoopNum
INTEGER :: LoopSideNum
REAL(r64) :: CpSourceSide ! local temporary for fluid specific heat
REAL(r64) :: CpLoadSide ! local temporary for fluid specific heat
! LOAD LOCAL VARIABLES FROM DATA STRUCTURE (for code readability)
PressureDrop = GSHP(GSHPNum)%CompSucPressDrop
ClearanceFactor = GSHP(GSHPNum)%CompClearanceFactor
PistonDisp = GSHP(GSHPNum)%CompPistonDisp
ShTemp = GSHP(GSHPNum)%SuperheatTemp
LosFac = GSHP(GSHPNum)%LossFactor
SourceSideUA = GSHP(GSHPNum)%SourceSideUACoeff
LoadSideUA = GSHP(GSHPNum)%LoadSideUACoeff
PowerLosses = GSHP(GSHPNum)%PowerLosses
HighPressCutOff = GSHP(GSHPNum)%HighPressCutOff
LowPressCutOff = GSHP(GSHPNum)%LowPressCutOff
! REPORT VAR
GSHPReport(GSHPNum)%Running = 0
! Init Module level Variables
GSHP(GSHPNum)%MustRun = .TRUE. ! Reset MustRun Flag to TRUE
LoadSideWaterMassFlowRate = 0.0d0 ! Load Side mass flow rate, water side
SourceSideWaterMassFlowRate = 0.0d0 ! Source Side mass flow rate, water side
Power = 0.0d0 ! power consumption
QLoad = 0.0d0 ! heat rejection from Load Side coil
QSource = 0.0d0
LoadSideInletNode = GSHP(GSHPNum)%LoadSideInletNodeNum
LoadSideOutletNode = GSHP(GSHPNum)%LoadSideOutletNodeNum
SourceSideInletNode = GSHP(GSHPNum)%SourceSideInletNodeNum
SourceSideOutletNode = GSHP(GSHPNum)%SourceSideOutletNodeNum
LoopNum = GSHP(GSHPNum)%LoadLoopNum
LoopSideNum = GSHP(GSHPNum)%LoadLoopSideNum
IF(PrevSimTime .NE. CurrentSimTime)THEN
PrevSimTime = CurrentSimTime
END IF
! CALCULATE THE SIMULATION TIME
CurrentSimTime = (dayofSim-1)*24 + hourofday-1 + (timestep-1)*timestepZone + SysTimeElapsed
! initialize event time array when the environment simulation begins
IF(CurrentSimTime == 0.0d0 .AND. OneTimeFlag)THEN
OneTimeFlag = .FALSE.
END IF
IF(CurrentSimTime > 0.0d0 )OneTimeFlag = .TRUE.
IF(MyLoad > 0.0d0)THEN
GSHP(GSHPNum)%MustRun = .TRUE.
GSHP(GSHPNum)%IsOn = .TRUE.
ELSE
GSHP(GSHPNum)%MustRun = .FALSE.
GSHP(GSHPNum)%IsOn = .FALSE.
END IF
!*******Set flow based on "run" flags**********
! Set flows if the heat pump is not running
IF( .NOT. GSHP(GSHPNum)%MustRun )THEN
LoadSideWaterMassFlowRate = 0.d0
Call SetComponentFlowRate(LoadSideWaterMassFlowRate, &
LoadSideInletNode, LoadSideOutletNode, &
GSHP(GSHPNum)%LoadLoopNum, GSHP(GSHPNum)%LoadLoopSideNum, &
GSHP(GSHPNum)%LoadBranchNum, GSHP(GSHPNum)%LoadCompNum)
SourceSideWaterMassFlowRate = 0.d0
Call SetComponentFlowRate(SourceSideWaterMassFlowRate, &
SourceSideInletNode, SourceSideOutletNode, &
GSHP(GSHPNum)%SourceLoopNum, GSHP(GSHPNum)%SourceLoopSideNum, &
GSHP(GSHPNum)%SourceBranchNum, GSHP(GSHPNum)%SourceCompNum)
!now initialize simulation variables for "heat pump off"
QLoad = 0.0d0
QSource = 0.0d0
Power = 0.0d0
LoadSideWaterInletTemp = Node(LoadSideInletNode)%Temp
LoadSideWaterOutletTemp = LoadSideWaterInletTemp
SourceSideWaterInletTemp = Node(SourceSideInletNode)%Temp
SourceSideWaterOutletTemp = SourceSideWaterInletTemp
RETURN !if heat pump is not running return without simulation
! Set flows if the heat pump is running
ELSE ! the heat pump must run, request design flow
LoadSideWaterMassFlowRate = GSHP(GSHPNum)%LoadSideDesignMassFlow
CALL SetComponentFlowRate(LoadSideWaterMassFlowRate, &
LoadSideInletNode, LoadSideOutletNode, &
GSHP(GSHPNum)%LoadLoopNum, GSHP(GSHPNum)%LoadLoopSideNum, &
GSHP(GSHPNum)%LoadBranchNum, GSHP(GSHPNum)%LoadCompNum)
SourceSideWaterMassFlowRate = GSHP(GSHPNum)%SourceSideDesignMassFlow
CALL SetComponentFlowRate(SourceSideWaterMassFlowRate, &
SourceSideInletNode, SourceSideOutletNode, &
GSHP(GSHPNum)%SourceLoopNum, GSHP(GSHPNum)%SourceLoopSideNum, &
GSHP(GSHPNum)%SourceBranchNum, GSHP(GSHPNum)%SourceCompNum)
! get inlet temps
LoadSideWaterInletTemp = Node(LoadSideInletNode)%Temp
SourceSideWaterInletTemp = Node(SourceSideInletNode)%Temp
!if there's no flow, turn the "heat pump off"
IF(LoadSideWaterMassFlowRate < MassFlowTolerance .OR. &
SourceSideWaterMassFlowRate < MassFlowTolerance)THEN
LoadSideWaterMassFlowRate = 0.d0
Call SetComponentFlowRate(LoadSideWaterMassFlowRate, &
LoadSideInletNode, LoadSideOutletNode, &
GSHP(GSHPNum)%LoadLoopNum, GSHP(GSHPNum)%LoadLoopSideNum, &
GSHP(GSHPNum)%LoadBranchNum, GSHP(GSHPNum)%LoadCompNum)
SourceSideWaterMassFlowRate = 0.d0
Call SetComponentFlowRate(SourceSideWaterMassFlowRate, &
SourceSideInletNode, SourceSideOutletNode, &
GSHP(GSHPNum)%SourceLoopNum, GSHP(GSHPNum)%SourceLoopSideNum, &
GSHP(GSHPNum)%SourceBranchNum, GSHP(GSHPNum)%SourceCompNum)
QLoad = 0.0d0
QSource = 0.0d0
Power = 0.0d0
LoadSideWaterInletTemp = Node(LoadSideInletNode)%Temp
LoadSideWaterOutletTemp = LoadSideWaterInletTemp
SourceSideWaterInletTemp = Node(SourceSideInletNode)%Temp
SourceSideWaterOutletTemp = SourceSideWaterInletTemp
RETURN
END IF
END IF
!***********BEGIN CALCULATION****************
! initialize the source and load side heat transfer rates for the simulation
initialQSource = 0.0d0
initialQLoad = 0.0d0
IterationCount = 0
CpSourceSide = GetSpecificHeatGlycol(PlantLoop(GSHP(GSHPNum)%SourceLoopNum)%FluidName, &
SourceSideWaterInletTemp, &
PlantLoop(GSHP(GSHPNum)%SourceLoopNum)%FluidIndex, &
'CalcGshpModel')
CpLoadSide = GetSpecificHeatGlycol(PlantLoop(GSHP(GSHPNum)%LoadLoopNum)%FluidName, &
LoadSideWaterInletTemp, &
PlantLoop(GSHP(GSHPNum)%LoadLoopNum)%FluidIndex, &
'CalcGshpModel')
! Determine effectiveness of Source Side (the Evaporator in heating mode)
SourceSideEffect = 1.0d0- EXP( -SourceSideUA / &
(CpSourceSide * SourceSideWaterMassFlowRate))
!Determine effectiveness of Load Side the condenser in heating mode
LoadSideEffect = 1.0d0- EXP ( -LoadSideUA / &
(CpLoadSide * LoadSideWaterMassFlowRate))
LOOPLoadEnth: DO ! main loop to solve model equations
IterationCount = IterationCount+1
! Determine Source Side tempertaure
SourceSideTemp = SourceSideWaterInletTemp - initialQSource/ &
(SourceSideEffect * CpSourceSide * SourceSideWaterMassFlowRate)
! To determine Load Side temperature condenser
LoadSideTemp = LoadSideWaterInletTemp + initialQLoad/ &
(LoadSideEffect * CpLoadSide * LoadSideWaterMassFlowRate)
! Determine the evaporating and condensing pressures
SourceSidePressure = GetSatPressureRefrig(GSHPRefrigerant,SourceSideTemp,GSHPRefrigIndex,'CalcGSHPModel:SourceSideTemp')
LoadSidePressure = GetSatPressureRefrig(GSHPRefrigerant,LoadSideTemp,GSHPRefrigIndex,'CalcGSHPModel:LoadSideTemp')
! check cutoff pressures
IF (SourceSidePressure < LowPressCutOff) THEN
CALL ShowSevereError(ModuleCompName//'="'//trim(GSHPName)//'" Heating Source Side Pressure Less than the Design Minimum')
CALL ShowContinueError('Source Side Pressure='//TRIM(TrimSigDigits(SourceSidePressure,2))// &
' and user specified Design Minimum Pressure='//TRIM(TrimSigDigits(LowPressCutoff,2)))
CALL ShowFatalError('Preceding Conditions cause termination.')
END IF
IF (LoadSidePressure > HighPressCutOff)THEN
CALL ShowSevereError(ModuleCompName//'="'//trim(GSHPName)//'" Heating Load Side Pressure greater than the Design Maximum')
CALL ShowContinueError('Load Side Pressure='//TRIM(TrimSigDigits(LoadSidePressure,2))// &
' and user specified Design Maximum Pressure='//TRIM(TrimSigDigits(HighPressCutOff,2)))
CALL ShowFatalError('Preceding Conditions cause termination.')
END IF
! Determine Suction Pressure at compressor inlet
SuctionPr = SourceSidePressure - PressureDrop
! Determine Discharge Pressure at compressor exit
DischargePr = LoadSidePressure + PressureDrop
! check cutoff pressures
IF (SuctionPr < LowPressCutOff) THEN
CALL ShowSevereError(ModuleCompName//'="'//trim(GSHPName)//'" Heating Suction Pressure Less than the Design Minimum')
CALL ShowContinueError('Heating Suction Pressure='//TRIM(TrimSigDigits(SuctionPr,2))// &
' and user specified Design Minimum Pressure='//TRIM(TrimSigDigits(LowPressCutoff,2)))
CALL ShowFatalError('Preceding Conditions cause termination.')
END IF
IF (DischargePr > HighPressCutOff)THEN
CALL ShowSevereError(ModuleCompName//'="'//trim(GSHPName)//'" Heating Discharge Pressure greater than the Design Maximum')
CALL ShowContinueError('Heating Discharge Pressure='//TRIM(TrimSigDigits(DischargePr,2))// &
' and user specified Design Maximum Pressure='//TRIM(TrimSigDigits(HighPressCutOff,2)))
CALL ShowFatalError('Preceding Conditions cause termination.')
END IF
! Determine the Source Side Outlet Enthalpy
qual=1.0d0
SourceSideOutletEnth = GetSatEnthalpyRefrig(GSHPRefrigerant, SourceSideTemp, qual, &
GSHPRefrigIndex,'CalcGSHPModel:SourceSideTemp')
! Determine Load Side Outlet Enthalpy
qual= 0.0d0
LoadSideOutletEnth = GetSatEnthalpyRefrig(GSHPRefrigerant,LoadSideTemp,qual, &
GSHPRefrigIndex,'CalcGSHPModel:LoadSideTemp')
! Determine superheated temperature of the Source Side outlet/compressor inlet
CompressInletTemp = SourceSideTemp + ShTemp
! Determine the enathalpy of the super heated fluid at Source Side outlet
SuperHeatEnth = GetSupHeatEnthalpyRefrig(GSHPRefrigerant,CompressInletTemp,SourceSidePressure, &
GSHPRefrigIndex,'CalcGSHPModel:CompressInletTemp')
! Determining the suction state of the fluid from inlet state involves interation
! Method employed...
! Determine the saturated temp at suction pressure, shoot out into the superheated region find the enthalpy
! check that with the inlet enthalpy ( as suction loss is isenthalpic). Iterate till desired accuracy is reached
CompSuctionSatTemp = GetSatTemperatureRefrig(GSHPRefrigerant,SuctionPr,GSHPRefrigIndex,'CalcGSHPModel:SuctionPr')
T110 = CompSuctionSatTemp
!Shoot into the super heated region
T111 = CompSuctionSatTemp + 80
! Iterate to find the Suction State - given suction pressure and superheat enthalpy
LOOP: DO
CompSuctionTemp = 0.5d0 * ( T110 + T111 )
CompSuctionEnth = GetSupHeatEnthalpyRefrig(GSHPRefrigerant,CompSuctionTemp,SuctionPr, &
GSHPRefrigIndex,'CalcGSHPModel:CompSuctionTemp')
IF (ABS(CompsuctionEnth-SuperHeatEnth)/SuperHeatEnth < .0001d0) THEN
EXIT LOOP
END IF
IF ( CompsuctionEnth < SuperHeatEnth ) THEN
T110 = CompSuctionTemp
ELSE
T111 = CompSuctionTemp
END IF
END DO LOOP
! Determine the Mass flow rate of refrigerant
CompSuctionDensity = GetSupHeatDensityRefrig(GSHPRefrigerant, CompSuctionTemp, SuctionPr, &
GSHPRefrigIndex,'CalcGSHPModel:CompSuctionTemp')
MassRef = PistonDisp * CompSuctionDensity * (1+ClearanceFactor-ClearanceFactor* &
((DischargePr/SuctionPr)**(1.d0/gamma)))
! Find the Source Side Heat Transfer
QSource = MassRef * ( SourceSideOutletEnth - LoadSideOutletEnth )
! Determine the theoretical power
Power = PowerLosses+(MassRef*gamma/(gamma-1) * SuctionPr /CompSuctionDensity/LosFac * &
((DischargePr/SuctionPr)**((gamma-1)/gamma) - 1))
! Determine the Loadside HeatRate (QLoad)
QLoad = Power + QSource
! convergence and iteration limit check
IF(ABS((QLoad - initialQLoad)/(initialQLoad+SmallNum)) < HeatBalTol .OR. IterationCount>IterationLimit) THEN
IF(IterationCount>IterationLimit)then
CALL ShowWarningError(ModuleCompName//' did not converge')
CALL ShowContinueErrorTimeStamp(' ')
CALL ShowContinueError('Heatpump Name = '//TRIM(GSHP(GSHPNum)%Name))
WRITE(ErrString,*) ABS(100.0d0*(QLoad - initialQLoad)/(initialQLoad+SmallNum))
CALL ShowContinueError('Heat Inbalance (%) = '//TRIM(ADJUSTL(ErrString)))
WRITE(ErrString,*) QLoad
CALL ShowContinueError('Load-side heat transfer rate = '//TRIM(ADJUSTL(ErrString)))
WRITE(ErrString,*) Qsource
CALL ShowContinueError('Source-side heat transfer rate = '//TRIM(ADJUSTL(ErrString)))
WRITE(ErrString,*) SourceSideWaterMassFlowRate
CALL ShowContinueError('Source-side mass flow rate = '//TRIM(ADJUSTL(ErrString)))
WRITE(ErrString,*) LoadSideWaterMassFlowRate
CALL ShowContinueError('Load-side mass flow rate = '//TRIM(ADJUSTL(ErrString)))
WRITE(ErrString,*) SourceSideWaterInletTemp
CALL ShowContinueError('Source-side inlet temperature = '//TRIM(ADJUSTL(ErrString)))
WRITE(ErrString,*) LoadSideWaterInletTemp
CALL ShowContinueError('Load-side inlet temperature = '//TRIM(ADJUSTL(ErrString)))
END IF
EXIT LOOPLoadEnth
ELSE ! update load
initialQLoad= initialQLoad+ RelaxParam*(QLoad-initialQLoad)
initialQSource = initialQSource + RelaxParam*(QSource - initialQSource)
END IF
END DO LOOPLoadEnth
!Control Strategy
IF(ABS(MyLoad) < QLoad) THEN
DutyFactor = ABS(MyLoad)/QLoad
QLoad = ABS(MyLoad)
Power = DutyFactor * Power
QSource = QSource * DutyFactor
! Determine the Exterior fluid temperature at the Load Side oulet and eveporator outlet...
! Refrigerant = "Steam"
LoadSideWaterOutletTemp = LoadSideWaterInletTemp + QLoad/(LoadSideWaterMassFlowRate * &
CpLoadSide)
SourceSideWaterOutletTemp = SourceSideWaterInletTemp - QSource/(SourceSideWaterMassFlowRate * &
CpSourceSide)
RETURN
END IF
LoadSideWaterOutletTemp = LoadSideWaterInletTemp + QLoad/(LoadSideWaterMassFlowRate * &
CpLoadSide)
SourceSideWaterOutletTemp = SourceSideWaterInletTemp - QSource/(SourceSideWaterMassFlowRate * &
CpSourceSide )
! REPORT VAR
GSHPReport(GSHPNum)%Running = 1
RETURN
END SUBROUTINE CalcGshpModel