CalcVerticalGroundHeatExchanger

Source Code

All Procedures

private subroutine CalcVerticalGroundHeatExchanger(GlheNum)

Help
Arguments

Type	Intent	Optional	Attributes		Name
integer				::	GlheNum
Source Code

SUBROUTINE CalcVerticalGroundHeatExchanger(GlheNum)
          ! SUBROUTINE INFORMATION:
          !       AUTHOR:          Dan Fisher
          !       DATE WRITTEN:    August, 2000
          !       MODIFIED         Arun Murugappan
          !       RE-ENGINEERED    na

          ! PURPOSE OF THIS SUBROUTINE:
          ! This is the main routine to simulate the operation of vertical
          ! closed-loop ground heat exchangers (GLHE).

          ! METHODOLOGY EMPLOYED:
          ! The borehole and fluid temperatures are calculated from the response to
          ! the current heat transfer rate and the response to the history of past
          ! applied heat pulses. The response to each pulse is calculated from a non-
          ! dimensionalized response function, or G-function, that is specific to the
          ! given borehole field arrangement, depth and spacing. The data defining
          ! this function is read from input.
          !
          ! The heat pulse histories need to be recorded over an extended period (months).
          ! To aid computational efficiency past pulses are continuously agregated into
          ! equivalent heat pulses of longer duration, as each pulse becomes less recent.


          ! REFERENCES:
          ! Eskilson, P. 'Thermal Analysis of Heat Extraction Boreholes' Ph.D. Thesis:
          !   Dept. of Mathematical Physics, University of Lund, Sweden, June 1987.
          ! Yavuzturk, C., J.D. Spitler. 1999. 'A Short Time Step Response Factor Model
          !   for Vertical Ground Loop Heat Exchangers. ASHRAE Transactions. 105(2): 475-485.

          ! USE STATEMENTS:
    USE DataPlant, ONLY : PlantLoop
    USE FluidProperties, ONLY: GetSpecificHeatGlycol, GetDensityGlycol

    IMPLICIT NONE    ! Enforce explicit typing of all variables in this routine

    ! SUBROUTINE ARGUMENT DEFINITIONS
    INTEGER                :: GlheNum


    !LOCAL BHORE HOLE PARAMETERS
    INTEGER                :: NumBholes
    REAL(r64)              :: FluidDensity
    REAL(r64)              :: BholeLength
    REAL(r64)              :: K_Ground
    REAL(r64)              :: K_Ground_Factor
    REAL(r64)              :: Cp_Fluid
    REAL(r64)              :: Tground
    REAL(r64)              :: ResistanceBhole           ! The thermal resistance of the borehole, (K per W/m]
    REAL(r64)              :: Gfuncval                  ! Interpolated G function value at a sub-hour
    REAL(r64)              :: ToutNew = 19.375d0
    REAL(r64)              :: FluidAveTemp
    REAL(r64)              :: GroundDiffusivity
    REAL(r64)              :: TimeSS                    !Steady state time
    REAL(r64)              :: TimeSS_Factor             !Steady state time factor for calculation
    REAL(r64)              :: XI
    REAL(r64)              :: C_1
    INTEGER                :: NumOfMonths               ! the number of months of simulation elapsed
    INTEGER                :: CurrentMonth              ! The Month upto which the Montly blocks are superposed
    REAL(r64)              :: SumQnMonthly              ! tmp variable which holds the sum of the Temperature diffrence
                                                        ! due to Aggregated heat extraction/rejection step
    REAL(r64)              :: SumQnHourly               ! same as above for hourly
    REAL(r64)              :: SumQnSubHourly            ! same as above for subhourly( with no aggreation]
    REAL(r64)              :: RQMonth
    REAL(r64)              :: RQHour
    REAL(r64)              :: RQSubHr
    INTEGER                :: I
    REAL(r64)              :: tmpQnSubHourly            ! current Qn subhourly value
    INTEGER                :: HourlyLimit               ! number of hours to be taken into account in superposition
    INTEGER                :: SubHourlyLimit            ! number of subhourlys to be taken into account in subhourly superposition
    REAL(r64)              :: SumTotal                  ! sum of all the Qn (load) blocks
    REAL(r64)              :: C0                        ! **Intermediate constants used
    REAL(r64)              :: C1                        ! **Intermediate constants used
    REAL(r64)              :: C2                        ! **in explicit  calcualtion of the
    REAL(r64)              :: C3                        ! **temperature at the U tube outlet.
    INTEGER, SAVE          :: PrevN =1                  ! The saved value of N at previous time step
    INTEGER                :: IndexN                    ! Used to index the LastHourN array
    LOGICAL,SAVE           :: UpdateCurSimTime = .TRUE. ! Used to reset the CurSimtime to reset after Warmupflag
    LOGICAL, SAVE          :: TriggerDesignDayReset = .FALSE.
    INTEGER                :: GlheInletNode             ! Inlet node number of the Glhe
    INTEGER                :: GlheOutletNode            ! Outlet node number of the Glhe
!    LOGICAL, SAVE      :: Allocated = .FALSE.
    INTEGER          :: AGG
    INTEGER          :: SubAGG
    INTEGER          :: LoopNum
    INTEGER          :: LoopSideNum


          !set local glhe parameters

  NumBholes    = VerticalGlhe(GlheNum)%NumBoreholes
  BholeLength  = VerticalGlhe(GlheNum)%BoreholeLength
  GlheInletNode = VerticalGlhe(GlheNum)%GlheInletNodeNum
  GlheInletTemp = Node(GlheInletnode)%Temp
  Cp_fluid     = GetSpecificHeatGlycol(PlantLoop(VerticalGlhe(GlheNum)%LoopNum)%FluidName, &
                                       GlheInletTemp, &
                                       PlantLoop(VerticalGlhe(GlheNum)%LoopNum)%FluidIndex, &
                                       'CalcVerticalGroundHeatExchanger')

  Tground      = VerticalGlhe(GlheNum)%TempGround
  FluidDensity = GetDensityGlycol(PlantLoop(VerticalGlhe(GlheNum)%LoopNum)%FluidName, &
                                       GlheInletTemp, &
                                       PlantLoop(VerticalGlhe(GlheNum)%LoopNum)%FluidIndex, &
                                       'CalcVerticalGroundHeatExchanger')
  K_Ground     = VerticalGlhe(GlheNum)%KGround
  K_Ground_Factor = 2.0d0 * PI * K_Ground
  AGG      = VerticalGlhe(GlheNum)%AGG
  SubAGG     = VerticalGlhe(GlheNum)%SubAGG
  GroundDiffusivity = VerticalGlhe(GlheNum)%KGround/VerticalGlhe(GlheNum)%CpRhoGround

  ! calculate annual time constant for ground conduction
  TimeSS = (VerticalGlhe(GlheNum)%BoreholeLength**2/(9.d0*GroundDiffusivity)) /SecInHour/8760.0d0
  TimeSS_Factor = TimeSS*8760.d0

  GlheOutletNode      = VerticalGlhe(GlheNum)%GlheOutletNodeNum
  LoopNum = VerticalGlhe(GlheNum)%LoopNum
  LoopSideNum = VerticalGlhe(GlheNum)%LoopSideNum

  GlheMassFlowRate = MDotActual

  IF(TriggerDesignDayReset.AND. Warmupflag)UpdateCurSimTime=.TRUE.
  IF(DayofSim .EQ. 1 .AND. UpdateCurSimTime)THEN
    CurrentSimTime = 0.0d0
    PrevTimeSteps = 0.0d0
    DO I = 1,NumVerticalGlhes
      VerticalGlhe(I)%QnHr           = 0.0d0
      VerticalGlhe(I)%QnMonthlyAgg   = 0.0d0
      VerticalGlhe(I)%QnSubHr   = 0.0d0
      VerticalGlhe(I)%LastHourN = 1
    END DO
    N = 1
    UpdateCurSimTime = .FALSE.
    TriggerDesignDayReset = .FALSE.
  END IF

  CurrentSimTime = (dayofSim-1)*24 + hourofday-1 + (timestep-1)*timestepZone + SysTimeElapsed  !+ TimeStepsys
  LocHourOfDay = MOD(CurrentSimTime,HrsPerDay)+1
  LocDayOfSim = CurrentSimTime/24 + 1

  IF(DayofSim .GT. 1) THEN
    UpdateCurSimTime = .TRUE.
  END IF

  IF(.NOT. Warmupflag)THEN
    TriggerDesignDayReset = .TRUE.
  ENDIF

  IF (CurrentSimTime .LE. 0.0d0)THEN
    PrevTimeSteps = 0.0d0 ! this resets history when rounding 24:00 hours during warmup avoids hard crash later
    GlheOutletTemp = GlheInletTemp
    GlheMassFlowRate = MDotActual
    CALL CalcAggregateLoad(GlheNum)     !Just allocates and initializes PrevHour array
    RETURN
  END IF

  ! Store currentsimtime in PrevTimeSteps only if a time step occurs

  IF( PrevTimeSteps (1) /= CurrentSimTime) THEN
    PrevTimeSteps = EOSHIFT(PrevTimeSteps,-1, CurrentSimTime)
    N = N + 1
  END IF
  IF(N /= PrevN) THEN
    PrevN = N
    DO I = 1,NumVerticalGlhes
      VerticalGlhe(I)%QnSubHr = EOSHIFT(VerticalGlhe(I)%QnSubHr,-1,LastQnSubHr(I))
    END DO
  END IF

  CALL CalcAggregateLoad(GlheNum)

  ! Update the borehole resistance each time
  CALL BoreholeResistance(GlheNum, ResistanceBhole)

  SumTotal = 0.0d0 !Objexx:Uninit Line added to assure SumTotal initialized when used in: GlheBoreholeTemp = TGround - SumTotal
  IF(N .EQ. 1) THEN
    IF(MDotActual .LE. 0.0d0) THEN
      tmpQnSubHourly = 0.0d0
      FluidAveTemp = Tground
      ToutNew = GlheInletTemp
    ELSE
      XI = LOG( CurrentSimTime / (TimeSS_Factor)  )
      CALL INTERP(GlheNum,XI,GfuncVal)

      C_1 = (BholeLength * NumBholes )/( 2.d0 * MDotActual * Cp_Fluid)
      tmpQnSubHourly = (Tground - GlheInletTemp)/ &
                           (GfuncVal/(K_Ground_Factor ) + ResistanceBhole + C_1)
      FluidAveTemp = Tground - tmpQnSubHourly * ResistanceBhole
      ToutNew = Tground - tmpQnsubHourly * (GfuncVal /  &
                                          (K_Ground_Factor  ) + ResistanceBhole - C_1)
    END IF
  ELSE
      ! no monthly super position
    IF(CurrentSimTime .LT.(HrsPerMonth + AGG + SubAGG)) THEN

      ! Calculate the Sub Hourly Superposition
      SumQnSubHourly = 0.0d0
      IF(INT(CurrentSimTime).LT.SubAGG)THEN
        IndexN = INT(CurrentSimTime)+1
      ELSE
        IndexN = SubAGG+1
      END IF
      SubHourlyLimit = N - VerticalGlhe(GlheNum)%LastHourN(IndexN) !Check this when running simulation

      SUBHRLY_LOOP:DO I = 1, SubHourlyLimit
        IF(I.EQ.SubHourlyLimit) THEN
          IF( INT(CurrentSimTime)>=SubAGG)THEN
            XI = LOG( (CurrentSimTime - PrevTimeSteps(I+1)) / (TimeSS_Factor) )
            CALL INTERP(GlheNum,XI,GfuncVal)
            RQSubHr =  GfuncVal/(K_Ground_Factor )
            SumQnSubHourly = SumQnSubHourly + (VerticalGlhe(GlheNum)%QnSubHr(I)-&
                                      VerticalGlhe(GlheNum)%QnHr(IndexN)) * RQSubHr
          ELSE
            XI = LOG( (CurrentSimTime - PrevTimeSteps(I+1)) / (TimeSS_Factor) )
            CALL INTERP(GlheNum,XI,GfuncVal)
            RQSubHr =  GfuncVal/(K_Ground_Factor)
            SumQnSubHourly = SumQnSubHourly + VerticalGlhe(GlheNum)%QnSubHr(I)* RQSubHr
          END IF
          EXIT SUBHRLY_LOOP
        END IF
          !PrevTimeSteps(I+1) This is "I+1" because PrevTimeSteps(1) = CurrentTimestep
        XI = LOG( (CurrentSimTime - PrevTimeSteps(I+1)) / (TimeSS_Factor) )
        CALL INTERP(GlheNum,XI,GfuncVal)
        RQSubHr = GfuncVal / (K_Ground_Factor)
        SumQnSubHourly = SumQnSubHourly + &
                           (VerticalGlhe(GlheNum)%QnSubHr(I)-&
                           VerticalGlhe(GlheNum)%QnSubHr(I+1)) * RQSubHr
      END DO SUBHRLY_LOOP

          ! Calculate the Hourly Superposition

      HourlyLimit =  INT(CurrentSimTime)
      SumQnHourly = 0.0d0
      HOURLY_LOOP:DO I = SubAGG+1, HourlyLimit
        IF(I.EQ.HourlyLimit) THEN
          XI = LOG( CurrentSimTime / (TimeSS_Factor) )
          CALL INTERP(GlheNum,XI,GfuncVal)
          RQHour =  GfuncVal/(K_Ground_Factor)
          SumQnHourly = SumQnHourly + VerticalGlhe(GlheNum)%QnHr(I)*RQHour
          EXIT HOURLY_LOOP
        END IF
        XI = LOG( (CurrentSimTime - INT(CurrentSimTime) + I) / (TimeSS_Factor))
        CALL INTERP(GlheNum,XI,GfuncVal)
        RQHour  = GfuncVal/(K_Ground_Factor )
        SumQnHourly = SumQnHourly + &
                            (VerticalGlhe(GlheNum)%QnHr(I)-&
                            VerticalGlhe(GlheNum)%QnHr(I+1)) * RQHour
      END DO HOURLY_LOOP


      ! Find the total Sum of the Temperature difference due to all load blocks
      SumTotal = SumQnSubHourly + SumQnHourly

       !Calulate the subhourly temperature due the Last Time steps Load
      XI = LOG( (CurrentSimTime - PrevTimeSteps(2)) / (TimeSS_Factor) )
      CALL INTERP(GlheNum,XI,GfuncVal)
      RQSubHr = GfuncVal / (K_Ground_Factor )

      IF(MDotActual .LE. 0.0d0)THEN
        tmpQnsubHourly = 0.0d0
        FluidAveTemp = Tground - SumTotal ! Q(N)*RB = 0
        ToutNew = GlheInletTemp
      ELSE
         !Dr.Spitler's Explicit set of equations to calculate the New Outlet Temperature of the U-Tube
        C0 = RQSubHr
        C1 = Tground - (SumTotal - VerticalGlhe(GlheNum)%QnSubHr(1)*RQSubHr)
        C2 = BholeLength * NumBholes / ( 2.d0 * MDotActual *Cp_Fluid)
        C3 = MDotActual * Cp_Fluid / ( BholeLength * NumBholes )
        tmpQnsubHourly = (C1 - GlheInletTemp)/( ResistanceBhole + C0 - C2 + ( 1 / C3 ))
        FluidAveTemp = C1 - (C0 + ResistanceBhole) * tmpQnsubHourly
        ToutNew = C1 + (C2 - C0 - ResistanceBhole) * tmpQnsubHourly
      END IF

    ELSE ! Monthly Aggregation and super position

      NumOfMonths = (CurrentSimTime+1)/HrsPerMonth

      IF(CurrentSimTime .LT. ((NumOfMonths)*HrsPerMonth)+AGG+SubAGG)THEN
        CurrentMonth = NumOfMonths - 1
      ELSE
        CurrentMonth = NumOfMonths
      END IF

         !monthly superposition
      SumQnMonthly = 0.0d0
      SUMMONTHLY :DO I = 1,CurrentMonth
        IF(I.EQ.1) THEN
          XI = LOG(CurrentSimTime/ (TimeSS_Factor))
          CALL INTERP(GlheNum,XI,GfuncVal)
          RQMonth =  GfuncVal/(K_Ground_Factor )
          SumQnMonthly = SumQnMonthly + VerticalGlhe(GlheNum)%QnMonthlyAgg(I) * RQMonth
          CYCLE SUMMONTHLY
        END IF
        XI = LOG((CurrentSimTime - (I-1)*HrsPerMonth )/ (TimeSS_Factor))
        CALL INTERP(GlheNum,XI,GfuncVal)
        RQMonth =  GfuncVal/(K_Ground_Factor )
        SumQnMonthly = SumQnMonthly + &
                         (VerticalGlhe(GlheNum)%QnMonthlyAgg(I)- &
                         VerticalGlhe(GlheNum)%QnMonthlyAgg(I-1))* RQMonth
      END DO SUMMONTHLY

       ! Hourly Supr position
      HourlyLimit = INT( CurrentSimTime - CurrentMonth* HrsPerMonth)
      SumQnHourly = 0.0d0
      HOURLYLOOP: DO I = 1 + SubAGG, HourlyLimit
        IF(I.EQ.HourlyLimit) THEN
          XI = LOG( (CurrentSimTime - INT(CurrentSimTime) + I) / (TimeSS_Factor) )
          CALL INTERP(GlheNum,XI,GfuncVal)
          RQHour =  GfuncVal/(K_Ground_Factor )
          SumQnHourly = SumQnHourly + (VerticalGlhe(GlheNum)%QnHr(I) -&
                                 VerticalGlhe(GlheNum)%QnMonthlyAgg(CurrentMonth)) * RQHour
          EXIT HOURLYLOOP
        END IF
        XI = LOG( (CurrentSimTime - INT(CurrentSimTime) + I) / (TimeSS_Factor) )
        CALL INTERP(GlheNum,XI,GfuncVal)
        RQHour  = GfuncVal / (K_Ground_Factor )
        SumQnHourly = SumQnHourly + &
                          (VerticalGlhe(GlheNum)%QnHr(I)-&
                          VerticalGlhe(GlheNum)%QnHr(I+1)) * RQHour
      END DO HOURLYLOOP

          ! Subhourly Superposition
      SubHourlyLimit = N - VerticalGlhe(GlheNum)%LastHourN(SubAGG+1)
      SumQnSubHourly = 0.0d0
      SUBHRLOOP:  DO I = 1 ,SubHourlyLimit
        IF(I.EQ.SubHourlyLimit) THEN
          XI = LOG( (CurrentSimTime - PrevTimeSteps(I+1)) / (TimeSS_Factor) )
          CALL INTERP(GlheNum,XI,GfuncVal)
          RQSubHr =  GfuncVal/(K_Ground_Factor )
          SumQnSubHourly = SumQnSubHourly + (VerticalGlhe(GlheNum)%QnSubHr(I)-&
                                      VerticalGlhe(GlheNum)%QnHr(SubAgg+1)) * RQSubHr
          EXIT SUBHRLOOP
        END IF
        XI = LOG( (CurrentSimTime - PrevTimeSteps(I+1)) / (TimeSS_Factor) )
        CALL INTERP(GlheNum,XI,GfuncVal)
            RQSubHr = GfuncVal / (K_Ground_Factor )
            SumQnSubHourly = SumQnSubHourly + &
                                 (VerticalGlhe(GlheNum)%QnSubHr(I)-&
                                  VerticalGlhe(GlheNum)%QnSubHr(I+1)) * RQSubHr
      END DO SUBHRLOOP

      SumTotal = SumQnMonthly + SumQnHourly + SumQnSubHourly

      !Calulate the subhourly temperature due the Last Time steps Load

      XI = LOG( (CurrentSimTime - PrevTimeSteps(2)) / (TimeSS_Factor) )
      CALL INTERP(GlheNum,XI,GfuncVal)
      RQSubHr = GfuncVal / (K_Ground_Factor )

      IF(MDotActual .LE. 0.0d0)THEN
        tmpQnsubHourly = 0.0d0
        FluidAveTemp = Tground - SumTotal ! Q(N)*RB = 0
        ToutNew = GlheInletTemp
      ELSE
      ! Explicit set of equations to calculate the New Outlet Temperature of the U-Tube
        C0 = RQSubHr
        C1 = Tground - (SumTotal - VerticalGlhe(GlheNum)%QnSubHr(1)*RQSubHr)
        C2 = BholeLength * NumBholes / ( 2 * MDotActual *Cp_Fluid)
        C3 = MDotActual * Cp_Fluid / ( BholeLength * NumBholes )
        tmpQnsubHourly = (C1 - GlheInletTemp)/( ResistanceBhole + C0 - C2 + ( 1 / C3 ))
        FluidAveTemp = C1 - (C0 + ResistanceBhole) * tmpQnsubHourly
        ToutNew = C1 + (C2 - C0 - ResistanceBhole) * tmpQnsubHourly
      END IF
    END IF !  end of AGG OR NO AGG
  END IF ! end of N  = 1 branch
  GlheBoreholeTemp = TGround - SumTotal
  !Load the QnSubHourly Array with a new value at end of every timestep

  !Load the report vars
  LastQnSubHr(GlheNum)  = tmpQnsubHourly
  GlheOutletTemp = ToutNew
  QGlhe = tmpQnsubHourly
  GlheAveFluidTemp = FluidAveTemp
  GlheRB = ResistanceBhole
  GlheMassFlowRate = MDotActual

END SUBROUTINE CalcVerticalGroundHeatExchanger
CalcVerticalGroundHeatExchanger Subroutine

Source Code

All Procedures

private subroutine CalcVerticalGroundHeatExchanger(GlheNum)

Uses:

Arguments

Calls

Called By

Source Code

Source Code

All Procedures
