InitialInitHeatBalFiniteDiff – EnergyPlusFortran

Source Code

All Procedures

private subroutine InitialInitHeatBalFiniteDiff()

Help
Arguments

None
Source Code

SUBROUTINE InitialInitHeatBalFiniteDiff

          ! SUBROUTINE INFORMATION:
          !       AUTHOR         Linda Lawrie
          !       DATE WRITTEN   March 2012
          !       MODIFIED       na
          !       RE-ENGINEERED  na

          ! PURPOSE OF THIS SUBROUTINE:
          ! This routine performs the original allocate, inits and setup output variables for the
          ! module.

          ! METHODOLOGY EMPLOYED:
          ! na

          ! REFERENCES:
          ! na

          ! USE STATEMENTS:
  USE General,      ONLY : TrimSigDigits, RoundSigDigits
  USE DataSurfaces, ONLY : HeatTransferModel_CondFD
  USE DataHeatBalance, Only: HighDiffusivityThreshold, ThinMaterialLayerThreshold

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

          ! SUBROUTINE ARGUMENT DEFINITIONS:
          ! na

          ! SUBROUTINE PARAMETER DEFINITIONS:
          ! na

          ! INTERFACE BLOCK SPECIFICATIONS:
          ! na

          ! DERIVED TYPE DEFINITIONS:
          ! na

          ! SUBROUTINE LOCAL VARIABLE DECLARATIONS:
  INTEGER :: Lay
  INTEGER :: SurfNum
  CHARACTER(len=MaxNameLength) :: LayChar

  REAL(r64)             :: dxn        ! Intermediate calculation of nodal spacing. This is the full dx. There is
                                                  ! a half dxn thick node at each surface. dxn is the "capacitor" spacing.
  INTEGER                           :: ipts1      ! Intermediate calculation for number of full thickness nodes per layer. There
                                                  ! are always two half nodes at the layer faces.
  INTEGER :: Layer      ! Loop counter
  INTEGER :: OutwardMatLayerNum ! layer index, layer outward of the current layer
  INTEGER :: layerNode
  INTEGER :: Delt
  INTEGER :: ConstrNum    ! Loop counter
  INTEGER :: TotNodes     ! Loop counter
  INTEGER :: CurrentLayer ! Loop counter
  INTEGER :: Surf         ! Loop counter
  INTEGER :: index        ! Loop Counters

  REAL(r64) :: Alpha
  REAL(r64) :: Malpha
  REAL(r64) :: stabilitytemp
  REAL(r64) :: stabilitymoist
  REAL(r64) :: a
  REAL(r64) :: b
  REAL(r64) :: c
  REAL(r64) :: d
  REAL(r64) :: kt
  REAL(r64) :: RhoS
  REAL(r64) :: Por
  REAL(r64) :: Cp
  REAL(r64) :: Dv
  LOGICAL :: errorsFound
  REAL(r64) :: DeltaTimestep      ! zone timestep in seconds, for local check of properties
  REAL(r64) :: ThicknessThreshold ! min thickness consistent with other thermal properties, for local check

  ALLOCATE(ConstructFD(TotConstructs))
  ALLOCATE(SigmaR(TotConstructs))
  ALLOCATE(SigmaC(TotConstructs))

  ALLOCATE(SurfaceFD(TotSurfaces))
  ALLOCATE(QHeatInFlux(TotSurfaces))
  ALLOCATE(QHeatOutFlux(TotSurfaces))
!  ALLOCATE(QFluxInArrivSurfCond(TotSurfaces))
!  ALLOCATE(QFluxOutArrivSurfCond(TotSurfaces))
!  ALLOCATE(OpaqSurfInsFaceConductionFlux(TotSurfaces))
!  ALLOCATE(OpaqSurfOutsideFaceConductionFlux(TotSurfaces))
!  ALLOCATE(QFluxZoneToInSurf(TotSurfaces))
!  ALLOCATE(QFluxOutsideToOutSurf(TotSurfaces))

 ! And then initialize
  QHeatInFlux  = 0.d0
  QHeatOutFlux  = 0.d0
  !QFluxZoneToInSurf  = 0.d0
  !QFluxOutsideToOutSurf  = 0.d0
  !QFluxInArrivSurfCond  = 0.d0
  !QFluxOutArrivSurfCond  = 0.d0
  OpaqSurfInsFaceConductionFlux = 0.d0
  OpaqSurfOutsideFaceConductionFlux = 0.d0

  ! Setup Output Variables

  !  set a Delt that fits the zone time step and keeps it below 200s.

  fracTimeStepZone_Hour=1.0d0/REAL(NumOfTimeStepInHour,r64)

  DO index = 1, 20
    Delt = (fracTimeStepZone_Hour*SecInHour)/index   ! TimeStepZone = Zone time step in fractional hours
    IF ( Delt <= 200 )  EXIT
  END DO

  DO ConstrNum = 1, TotConstructs
    !Need to skip window constructions and eventually window materials
    IF(Construct(ConstrNum)%TypeIsWindow) CYCLE

    ALLOCATE(ConstructFD(ConstrNum)%Name(Construct(ConstrNum)%TotLayers))
    ALLOCATE(ConstructFD(ConstrNum)%Thickness(Construct(ConstrNum)%TotLayers))
    ALLOCATE(ConstructFD(ConstrNum)%NodeNumPoint(Construct(ConstrNum)%TotLayers))
    ALLOCATE(ConstructFD(ConstrNum)%DelX(Construct(ConstrNum)%TotLayers))
    ALLOCATE(ConstructFD(ConstrNum)%TempStability(Construct(ConstrNum)%TotLayers))
    ALLOCATE(ConstructFD(ConstrNum)%MoistStability(Construct(ConstrNum)%TotLayers))
      ! Node Number of the Interface node following each layer
!    ALLOCATE(ConstructFD(ConstrNum)%InterfaceNodeNums(Construct(ConstrNum)%TotLayers))


    TotNodes = 0
    SigmaR(ConstrNum) =0.0d0
    SigmaC(ConstrNum) =0.0d0

    DO Layer = 1, Construct(ConstrNum)%TotLayers    ! Begin layer loop ...

          ! Loop through all of the layers in the current construct. The purpose
          ! of this loop is to define the thermal properties  and to.
          ! determine the total number of full size nodes in each layer.
          ! The number of temperature points is one more than this
          ! because of the two half nodes at the layer faces.
          ! The calculation of dxn used here is based on a standard stability
          ! criteria for explicit finite difference solutions.  This criteria
          ! was chosen not because it is viewed to be correct, but rather for
          ! lack of any better criteria at this time.  The use of a Fourier
          ! number based criteria such as this is probably physically correct.
          !  Change to implicit formulation still uses explicit stability, but
          ! now there are special equations for R-only layers.

      CurrentLayer = Construct(ConstrNum)%LayerPoint(Layer)


      ConstructFD(ConstrNum)%Name(Layer) = Material(CurrentLayer)%Name
      ConstructFD(ConstrNum)%Thickness(Layer) = Material(CurrentLayer)%Thickness

     ! Do some quick error checks for this section.

      IF (Material(CurrentLayer)%ROnly  ) THEN  ! Rlayer

        !  These values are only needed temporarily and to calculate flux,
        !   Layer will be handled
        !  as a pure R in the temperature calc.
        ! assign other properties based on resistance

        Material(CurrentLayer)%SpecHeat  = 0.0001d0
        Material(CurrentLayer)%Density = 1.d0
        Material(currentLayer)%Thickness = 0.1d0  !  arbitrary thickness for R layer
        Material(CurrentLayer)%Conductivity= &
        Material(currentLayer)%Thickness/Material(CurrentLayer)%Resistance
        kt = Material(CurrentLayer)%Conductivity
        ConstructFD(ConstrNum)%Thickness(Layer) = Material(CurrentLayer)%Thickness

        SigmaR(ConstrNum) = SigmaR(ConstrNum) + Material(CurrentLayer)%Resistance  ! add resistance of R layer
        SigmaC(ConstrNum) = SigmaC(ConstrNum)  + 0.0d0                               !  no capacitance for R layer

        Alpha = kt/(Material(CurrentLayer)%Density*Material(CurrentLayer)%SpecHeat)

        mAlpha = 0.d0

      ELSE IF (Material(CurrentLayer)%Group == 1 ) THEN   !  Group 1 = Air

                !  Again, these values are only needed temporarily and to calculate flux,
                !   Air layer will be handled
                !  as a pure R in the temperature calc.
                !
                ! assign
                ! other properties based on resistance

        Material(CurrentLayer)%SpecHeat  = 0.0001d0
        Material(CurrentLayer)%Density = 1.d0
        Material(currentLayer)%Thickness = 0.1d0  !  arbitrary thickness for R layer
        Material(currentLayer)%Conductivity= &
        Material(currentLayer)%Thickness/Material(CurrentLayer)%Resistance
        kt = Material(CurrentLayer)%Conductivity
        ConstructFD(ConstrNum)%Thickness(Layer) = Material(CurrentLayer)%Thickness

        SigmaR(ConstrNum) = SigmaR(ConstrNum) + Material(CurrentLayer)%Resistance  ! add resistance of R layer
        SigmaC(ConstrNum) = SigmaC(ConstrNum)  + 0.0d0                               !  no capacitance for R layer

        Alpha = kt/(Material(CurrentLayer)%Density*Material(CurrentLayer)%SpecHeat)
        mAlpha = 0.d0
      ELSE IF (Construct(ConstrNum)%TypeIsIRT) THEN  ! make similar to air? (that didn't seem to work well)
        CALL ShowSevereError('InitHeatBalFiniteDiff: Construction ="'//trim(Construct(ConstrNum)%Name)//  &
             '" uses Material:InfraredTransparent. Cannot be used currently with finite difference calculations.')
        IF (Construct(ConstrNum)%IsUsed) THEN
          CALL ShowContinueError('...since this construction is used in a surface, the simulation is not allowed.')
          errorsFound=.true.
        ELSE
          CALL ShowContinueError('...if this construction were used in a surface, the simulation would be terminated.')
        ENDIF
        CYCLE
!            ! set properties to get past other initializations.
!          Material(CurrentLayer)%SpecHeat  = 0.0001d0
!          Material(CurrentLayer)%Density = 0.0001d0
!          Material(currentLayer)%Thickness = 0.1d0  !  arbitrary thickness for R layer
!          Material(currentLayer)%Conductivity= &
!    -      Material(currentLayer)%Thickness/Material(CurrentLayer)%Resistance
!          kt = Material(CurrentLayer)%Conductivity
!          ConstructFD(ConstrNum)%Thickness(Layer) = Material(CurrentLayer)%Thickness
!
!          SigmaR(ConstrNum) = SigmaR(ConstrNum) + Material(CurrentLayer)%Resistance  ! add resistance of R layer
!          SigmaC(ConstrNum) = SigmaC(ConstrNum)  + 0.0                               !  no capacitance for R layer
!
!          Alpha = kt/(Material(CurrentLayer)%Density*Material(CurrentLayer)%SpecHeat)
!          mAlpha = 0.d0
      ELSE
       !    Regular material Properties
        a = Material(CurrentLayer)%MoistACoeff
        b = Material(CurrentLayer)%MoistBCoeff
        c = Material(CurrentLayer)%MoistCCoeff
        d = Material(CurrentLayer)%MoistDCoeff
        kt = Material(CurrentLayer)%Conductivity
        RhoS = Material(CurrentLayer)%Density
        Por = Material(CurrentLayer)%Porosity
        Cp = Material(CurrentLayer)%SpecHeat
         ! Need Resistance for reg layer
        Material(CurrentLayer)%Resistance = Material(currentLayer)%Thickness/Material(CurrentLayer)%Conductivity
        Dv = Material(CurrentLayer)%VaporDiffus
        SigmaR(ConstrNum) = SigmaR(ConstrNum) + Material(CurrentLayer)%Resistance  ! add resistance
        SigmaC(ConstrNum) = SigmaC(ConstrNum)  + &
        Material(CurrentLayer)%Density*Material(CurrentLayer)%SpecHeat*Material(CurrentLayer)%Thickness
        Alpha = kt/(RhoS*Cp)
        !   Alpha = kt*(Por+At*RhoS)/(RhoS*(Bv*Por*Lambda+Cp*(Por+At*RhoS)))
        MAlpha = 0.d0

        !check for Material layers that are too thin and highly conductivity (not appropriate for surface models)
        IF (Alpha > HighDiffusivityThreshold .AND. .NOT. Material(CurrentLayer)%WarnedForHighDiffusivity) THEN
          DeltaTimestep      = TimeStepZone * SecInHour
          ThicknessThreshold = SQRT(Alpha * DeltaTimestep * 3.d0)
          IF (Material(CurrentLayer)%Thickness < ThicknessThreshold) THEN
            CALL ShowSevereError('InitialInitHeatBalFiniteDiff: Found Material that is too thin and/or too highly conductive,' &
                                    //' material name = '  // TRIM(Material(CurrentLayer)%Name))
            CALL ShowContinueError('High conductivity Material layers are not well supported by Conduction Finite Difference, ' &
                                    //' material conductivity = ' //TRIM(RoundSigDigits(Material(CurrentLayer)%Conductivity, 3)) &
                                    //' [W/m-K]')
            CALL ShowContinueError('Material thermal diffusivity = ' //TRIM(RoundSigDigits(Alpha, 3)) //' [m2/s]')
            CALL ShowContinueError('Material with this thermal diffusivity should have thickness > ' &
                                      //  TRIM(RoundSigDigits(ThicknessThreshold , 5)) // ' [m]')
            IF (Material(CurrentLayer)%Thickness < ThinMaterialLayerThreshold) THEN
              CALL ShowContinueError('Material may be too thin to be modeled well, thickness = ' &
                                      // TRIM(RoundSigDigits(Material(currentLayer)%Thickness, 5 ))//' [m]')
              CALL ShowContinueError('Material with this thermal diffusivity should have thickness > ' &
                                      //  TRIM(RoundSigDigits(ThinMaterialLayerThreshold , 5)) // ' [m]')
            ENDIF
            Material(CurrentLayer)%WarnedForHighDiffusivity = .TRUE.
          ENDIF
        ENDIF

      END IF  !  R, Air  or regular material properties and parameters

     ! Proceed with setting node sizes in layers

      DXN=SQRT(Alpha*Delt*SpaceDescritConstant)  ! The Fourier number is set using user constant

       ! number of nodes=thickness/spacing.  This is number of full size node spaces across layer.
      Ipts1=INT(Material(CurrentLayer)%Thickness/DXN)
       !  set high conductivity layers to a single full size node thickness. (two half nodes)
      IF(Ipts1 <= 1) Ipts1 = 1
      IF (Material(CurrentLayer)%ROnly .or. Material(CurrentLayer)%Group == 1 ) THEN

        Ipts1 = 1  !  single full node in R layers- surfaces of adjacent material or inside/outside layer
      END If

      Dxn=Material(CurrentLayer)%Thickness/REAL(Ipts1,r64)  ! full node thickness

      StabilityTemp = Alpha*Delt/dxn**2
      StabilityMoist = malpha*Delt/dxn**2
      ConstructFD(ConstrNum)%TempStability(Layer) = StabilityTemp
      ConstructFD(ConstrNum)%MoistStability(Layer) = StabilityMoist
      ConstructFD(ConstrNum)%Delx(Layer) = DXN

      TotNodes = TotNodes + Ipts1   !  number of full size nodes
      ConstructFD(ConstrNum)%NodeNumPoint(Layer) = Ipts1   !  number of full size nodes
    END DO   !  end of layer loop.


    ConstructFD(ConstrNum)%TotNodes = TotNodes
    ConstructFD(ConstrNum)%DeltaTime = Delt


  END DO            ! End of Construction Loop.  TotNodes in each construction now set

  ! now determine x location, or distance that nodes are from the outside face in meters
  DO ConstrNum = 1, TotConstructs
     IF (ConstructFD(ConstrNum)%TotNodes > 0) THEN
       ALLOCATE(ConstructFD(ConstrNum)%NodeXlocation(ConstructFD(ConstrNum)%TotNodes + 1 ))
       ConstructFD(ConstrNum)%NodeXlocation = 0.d0 ! init them all
       Ipts1 = 0 ! init counter
       DO Layer = 1, Construct(ConstrNum)%TotLayers
         OutwardMatLayerNum = Layer - 1
         DO layerNode = 1, ConstructFD(ConstrNum)%NodeNumPoint(Layer)
           Ipts1 = Ipts1 +1
           IF (Ipts1 == 1) THEN
             ConstructFD(ConstrNum)%NodeXlocation(Ipts1) = 0.d0 ! first node is on outside face

           ELSEIF (LayerNode == 1) THEN
             IF ( OutwardMatLayerNum > 0 .AND. OutwardMatLayerNum <= Construct(ConstrNum)%TotLayers ) THEN
             ! later nodes are Delx away from previous, but use Delx from previous layer
               ConstructFD(ConstrNum)%NodeXlocation(Ipts1) = ConstructFD(ConstrNum)%NodeXlocation(Ipts1 - 1) &
                                                           + ConstructFD(ConstrNum)%Delx(OutwardMatLayerNum)
             ENDIF
           ELSE
             ! later nodes are Delx away from previous
             ConstructFD(ConstrNum)%NodeXlocation(Ipts1) = ConstructFD(ConstrNum)%NodeXlocation(Ipts1 - 1) &
                                                           + ConstructFD(ConstrNum)%Delx(Layer)
           ENDIF

         ENDDO
       ENDDO
       Layer = Construct(ConstrNum)%TotLayers
       Ipts1 = Ipts1 +1
       ConstructFD(ConstrNum)%NodeXlocation(Ipts1) = ConstructFD(ConstrNum)%NodeXlocation(Ipts1 - 1) &
                                                           + ConstructFD(ConstrNum)%Delx(Layer)
     ENDIF
  ENDDO

  DO Surf = 1,TotSurfaces
    IF (.not. Surface(Surf)%HeatTransSurf) CYCLE
    IF (Surface(Surf)%Class == SurfaceClass_Window) CYCLE
    IF (Surface(Surf)%HeatTransferAlgorithm /= HeatTransferModel_CondFD) CYCLE
!    IF(Surface(Surf)%Construction <= 0) CYCLE  ! Shading surface, not really a heat transfer surface
    ConstrNum=Surface(surf)%Construction
!    IF(Construct(ConstrNum)%TypeIsWindow) CYCLE  !  Windows simulated in Window module
    TotNodes = ConstructFD(ConstrNum)%TotNodes

     !Allocate the Surface Arrays
    ALLOCATE(SurfaceFD(Surf)%T(TotNodes + 1))
    ALLOCATE(SurfaceFD(Surf)%TOld(TotNodes + 1))
    ALLOCATE(SurfaceFD(Surf)%TT(TotNodes + 1))
    ALLOCATE(SurfaceFD(Surf)%Rhov(TotNodes + 1))
    ALLOCATE(SurfaceFD(Surf)%RhovOld(TotNodes + 1))
    ALLOCATE(SurfaceFD(Surf)%RhoT(TotNodes + 1))
    ALLOCATE(SurfaceFD(Surf)%TD(TotNodes + 1))
    ALLOCATE(SurfaceFD(Surf)%TDT(TotNodes + 1))
    ALLOCATE(SurfaceFD(Surf)%TDTLast(TotNodes + 1))
    ALLOCATE(SurfaceFD(Surf)%TDOld(TotNodes + 1))
    ALLOCATE(SurfaceFD(Surf)%TDreport(TotNodes + 1))
    ALLOCATE(SurfaceFD(Surf)%RH(TotNodes + 1))
    ALLOCATE(SurfaceFD(Surf)%RHreport(TotNodes + 1))
    ALLOCATE(SurfaceFD(Surf)%EnthOld(TotNodes + 1))
    ALLOCATE(SurfaceFD(Surf)%EnthNew(TotNodes + 1))
    ALLOCATE(SurfaceFD(Surf)%EnthLast(TotNodes + 1))


    !Initialize the allocated arrays.
    SurfaceFD(Surf)%T              = TempInitValue
    SurfaceFD(Surf)%TOld           = TempInitValue
    SurfaceFD(Surf)%TT             = TempInitValue
    SurfaceFD(Surf)%Rhov           = RhovInitValue
    SurfaceFD(Surf)%RhovOld        = RhovInitValue
    SurfaceFD(Surf)%RhoT           = RhovInitValue
    SurfaceFD(Surf)%TD             = TempInitValue
    SurfaceFD(Surf)%TDT            = TempInitValue
    SurfaceFD(Surf)%TDTLast        = TempInitValue
    SurfaceFD(Surf)%TDOld          = TempInitValue
    SurfaceFD(Surf)%TDreport       = TempInitValue
    SurfaceFD(Surf)%RH             = 0.0d0
    SurfaceFD(Surf)%RHreport       = 0.0d0
    SurfaceFD(Surf)%EnthOld        = EnthInitValue
    SurfaceFD(Surf)%EnthNew        = EnthInitValue
    SurfaceFD(Surf)%EnthLast       = EnthInitValue
  END DO


  DO SurfNum=1,TotSurfaces
    IF (.not. Surface(SurfNum)%HeatTransSurf) CYCLE
    IF (Surface(SurfNum)%Class == SurfaceClass_Window) CYCLE
    IF (Surface(SurfNum)%HeatTransferAlgorithm /= HeatTransferModel_CondFD) CYCLE
 !   If(SolutionAlgo == UseCondFD .or. SolutionAlgo == UseCondFDSimple)Then
!
!    CALL SetupOutputVariable('CondFD Outside Surface Heat Flux [W/m2]',   QFluxOutArrivSurfCond(SurfNum), &
!                             'Zone','State',TRIM(Surface(SurfNum)%Name))
!    CALL SetupOutputVariable('CondFD Inside Surface Heat Flux [W/m2]',    QFluxInArrivSurfCond(SurfNum), &
!                             'Zone','State',TRIM(Surface(SurfNum)%Name))
!    CALL SetupOutputVariable('CondFD Outside Heat Flux to Surface [W/m2]',QFluxOutsideToOutSurf(SurfNum), &
!                             'Zone','State',TRIM(Surface(SurfNum)%Name))
!    CALL SetupOutputVariable('CondFD Inside Heat Flux to Surface [W/m2]', QFluxZoneToInSurf(SurfNum), &
!                             'Zone','State',TRIM(Surface(SurfNum)%Name))

    CALL SetupOutputVariable('CondFD Inner Solver Loop Iteration Count [ ]', SurfaceFD(SurfNum)%GSloopCounter, &
                             'Zone','Sum',Surface(SurfNum)%Name)


    TotNodes = ConstructFD(Surface(SurfNum)%Construction)%TotNodes  ! Full size nodes, start with outside face.
    DO Lay = 1,TotNodes +1  ! include inside face node
      CALL SetupOutputVariable('CondFD Surface Temperature Node '//TRIM(TrimSigDigits(Lay))//' [C]',&
           SurfaceFD(SurfNum)%TDreport(Lay),  &
          'Zone','State', Surface(SurfNum)%Name)
    END DO

  ENDDO  ! End of the Surface Loop for Report Variable Setup

  CALL ReportFiniteDiffInits  ! Report the results from the Finite Diff Inits

  RETURN

END SUBROUTINE InitialInitHeatBalFiniteDiff
InitialInitHeatBalFiniteDiff Subroutine

Source Code

All Procedures

private subroutine InitialInitHeatBalFiniteDiff()

Uses:

Arguments

Calls

Called By

Source Code

Source Code

All Procedures
