| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| integer, | intent(in) | :: | Delt | |||
| integer, | intent(in) | :: | I | |||
| integer, | intent(in) | :: | Lay | |||
| integer, | intent(in) | :: | Surf | |||
| real(kind=r64), | intent(in), | DIMENSION(:) | :: | T | ||
| real(kind=r64), | intent(inout), | DIMENSION(:) | :: | TT | ||
| real(kind=r64), | intent(in), | DIMENSION(:) | :: | Rhov | ||
| real(kind=r64), | intent(inout), | DIMENSION(:) | :: | RhoT | ||
| real(kind=r64), | intent(inout), | DIMENSION(:) | :: | RH | ||
| real(kind=r64), | intent(in), | DIMENSION(:) | :: | TD | ||
| real(kind=r64), | intent(inout), | DIMENSION(:) | :: | TDT | ||
| real(kind=r64), | intent(in), | DIMENSION(:) | :: | EnthOld | ||
| real(kind=r64), | intent(inout), | DIMENSION(:) | :: | EnthNew | ||
| integer, | intent(in) | :: | GSiter |
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 IntInterfaceNodeEqns(Delt,I,Lay,Surf,T,TT,Rhov,RhoT,RH,TD,TDT,EnthOld,EnthNew,GSiter)
! SUBROUTINE INFORMATION:
! AUTHOR Richard Liesen
! DATE WRITTEN November, 2003
! MODIFIED May 2011, B. Griffith, P. Tabares, add first order fully implicit, bug fixes, cleanup
! RE-ENGINEERED Curtis Pedersen, Changed to Implit mode and included enthalpy. FY2006
! PURPOSE OF THIS SUBROUTINE:
! calculate finite difference heat transfer for nodes that interface two different material layers inside construction
! METHODOLOGY EMPLOYED:
! na
! REFERENCES:
! na
! USE STATEMENTS:
! na
IMPLICIT NONE ! Enforce explicit typing of all variables in this routine
! SUBROUTINE ARGUMENT DEFINITIONS:
INTEGER, INTENT(IN) :: Delt ! Time Increment
INTEGER, INTENT(IN) :: I ! Node Index
INTEGER, INTENT(IN) :: Lay ! Layer Number for Construction
INTEGER, INTENT(IN) :: Surf ! Surface number
INTEGER, INTENT(IN) :: GSiter ! Iteration number of Gauss Seidell iteration
REAL(r64),DIMENSION(:), INTENT(In) :: T !INSIDE SURFACE TEMPERATURE OF EACH HEAT TRANSFER SURF.
REAL(r64),DIMENSION(:), INTENT(InOut) :: TT !INSIDE SURFACE TEMPERATURE OF EACH HEAT TRANSFER SURF.
REAL(r64),DIMENSION(:), INTENT(In) :: Rhov !INSIDE SURFACE TEMPERATURE OF EACH HEAT TRANSFER SURF.
REAL(r64),DIMENSION(:), INTENT(InOut) :: RhoT !INSIDE SURFACE TEMPERATURE OF EACH HEAT TRANSFER SURF.
REAL(r64),DIMENSION(:), INTENT(In) :: TD !OLD NODE TEMPERATURES OF EACH HEAT TRANSFER SURF IN CONDFD.
REAL(r64),DIMENSION(:), INTENT(InOut) :: TDT !NEW NODE TEMPERATURES OF EACH HEAT TRANSFER SURF IN CONDFD.
REAL(r64),DIMENSION(:), INTENT(InOut) :: RH !RELATIVE HUMIDITY.
REAL(r64),DIMENSION(:), INTENT(In) :: EnthOld ! Old Nodal enthalpy
REAL(r64),DIMENSION(:), INTENT(InOut) :: EnthNew ! New Nodal enthalpy
! SUBROUTINE PARAMETER DEFINITIONS:
! REAL(r64), PARAMETER :: NinetyNine=99.0d0
! INTERFACE BLOCK SPECIFICATIONS:
! na
! DERIVED TYPE DEFINITIONS:
! na
! SUBROUTINE LOCAL VARIABLE DECLARATIONS:
INTEGER :: ConstrNum
INTEGER :: MatLay
INTEGER :: MatLay2
INTEGER :: IndVarCol
INTEGER :: DepVarCol
REAL(r64) :: kt1
REAL(r64) :: kt2
REAL(r64) :: kt1o
REAL(r64) :: kt2o
REAL(r64) :: kt11
REAL(r64) :: kt21
REAL(r64) :: Cp1
REAL(r64) :: Cp2
REAL(r64) :: Cpo1
REAL(r64) :: Cpo2
REAL(r64) :: RhoS1
REAL(r64) :: RhoS2
REAL(r64) :: Delx1
REAL(r64) :: Delx2
REAL(r64) :: Enth1New
REAL(r64) :: Enth2New
REAL(r64) :: Enth1Old
REAL(r64) :: Enth2Old
REAL(r64) :: Rlayer ! resistance value of R Layer
REAL(r64) :: Rlayer2 ! resistance value of next layer to inside
REAL(r64) :: QSSFlux ! Source/Sink flux value at a layer interface
LOGICAL :: RlayerPresent = .FALSE.
LOGICAL :: RLayer2Present = .FALSE.
ConstrNum=Surface(surf)%Construction
MatLay = Construct(ConstrNum)%LayerPoint(Lay)
MatLay2 = Construct(ConstrNum)%LayerPoint(Lay+1)
! Set Thermal Conductivity. Can be constant, simple linear temp dep or multiple linear segment temp function dep.
kt1o = Material(MatLay)%Conductivity
kt11= MaterialFD(MatLay)%tk1
kt1= kt1o + kt11*((TDT(I)+TDT(I-1))/2.d0 -20.d0)
IF( SUM(MaterialFD(MatLay)%TempCond(1:3,2)) >= 0.0d0) THEN ! Multiple Linear Segment Function
IndVarCol = 1 ! temperature
DepVarCol = 2 ! thermal conductivity
kt1 = terpld(MaterialFD(MatLay)%numTempCond,MaterialFD(MatLay)%TempCond,(TDT(I)+TDT(I-1))/2.d0,IndVarCol,DepVarCol)
ENDIF
RhoS1 = Material(MatLay)%Density
Cpo1 = Material(MatLay)%SpecHeat ! constant Cp from input file
Delx1 = ConstructFD(ConstrNum)%Delx(Lay)
Rlayer = Material(MatLay)%Resistance
kt2o = Material(MatLay2)%Conductivity
kt21= MaterialFD(MatLay2)%tk1
kt2= kt2o + kt21*((TDT(I)+TDT(I+1))/2.d0 -20.d0)
IF( SUM(MaterialFD(MatLay2)%TempCond(1:3,2)) >= 0.0d0) THEN ! Multiple Linear Segment Function
IndVarCol = 1 ! temperature
DepVarCol = 2 ! thermal conductivity
kt2 =terpld(MaterialFD(MatLay2)%numTempCond,MaterialFD(MatLay2)%TempCond,(TDT(I)+TDT(I+1))/2.d0,IndVarCol,DepVarCol)
ENDIF
RhoS2 = Material(MatLay2)%Density
Cpo2 = Material(MatLay2)%SpecHeat
Delx2 = ConstructFD(ConstrNum)%Delx(Lay+1)
Rlayer2 = Material(MatLay2)%Resistance
Cp1=Cpo1 ! Will be reset if PCM
Cp2=Cpo2 ! will be reset if PCM
IF ( Surface(Surf)%HeatTransferAlgorithm == HeatTransferModel_CondFD )THEN !HT Algo issue
!Calculate the Dry Heat Conduction Equation
RLayerPresent = .FALSE.
RLayer2Present = .false.
! Source/Sink Flux Capability ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
QSSFlux = 0.0d0
IF (Surface(Surf)%Area >0.0d0 .and. &
Construct(ConstrNum)%SourceSinkPresent .and. Lay == Construct(ConstrNum)%SourceAfterLayer) then
QSSFlux = QRadSysSource(Surf)/Surface(Surf)%Area &
+ QPVSysSource(Surf)/Surface(Surf)%Area ! Includes QPV Source
END IF
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
IF (Material(MatLay)%ROnly .or. Material(MatLay)%Group == 1) RLayerPresent = .TRUE.
IF (Material(MatLay2)%ROnly .or. Material(MatLay2)%Group ==1) Rlayer2Present =.TRUE.
IF ( RLayerPresent .and. RLayer2Present) THEN
TDT(I)=(Rlayer2*TDT(I-1) + Rlayer*TDT(I+1))/(Rlayer + Rlayer2) ! two adjacent R layers
ELSE IF ( RLayerPresent .and. .not. RLayer2Present ) THEN ! R-layer first
! Check for PCM second layer
IndVarCol = 1 ! temperature
DepVarCol = 2 ! thermal conductivity
IF( Sum(MaterialFD(MatLay)%TempEnth(1:3,2)) < 0.0d0 &
.and. Sum(MaterialFD(MatLay2)%TempEnth(1:3,2)) > 0.0d0) THEN ! phase change material Layer2, Use TempEnth Data
Enth2Old = terpld(MaterialFD(MatLay2)%numTempEnth,MaterialFD(MatLay2)%TempEnth,TD(I),IndVarCol,DepVarCol)
Enth2New = terpld(MaterialFD(MatLay2)%numTempEnth,MaterialFD(MatLay2)%TempEnth,TDT(I),IndVarCol,DepVarCol)
EnthNew(I) = Enth2New ! This node really doesn't have an enthalpy, this gives it a value
IF (ABS(Enth2New-Enth2Old) <= smalldiff .or. ABS(TDT(I)-TD(I)) <= smalldiff) THEN
Cp2 = Cpo2
ELSE
Cp2 = MAX(Cpo2,(Enth2New -Enth2Old)/(TDT(I)-TD(I)))
END IF
ENDIF
! R layer first, then PCM or regular layer.
SELECT CASE (CondFDSchemeType)
CASE (CrankNicholsonSecondOrder)
TDT(I) =( 2.d0*delt*Delx2*QSSFlux*Rlayer - delt*Delx2*TD(I) - delt*kt2*Rlayer*TD(I) + &
Cp2*Delx2**2*RhoS2*Rlayer*TD(I) + delt*Delx2*TD(I-1) + delt*kt2*Rlayer*TD(I+1) + delt*Delx2*TDT(I-1) + &
delt*kt2*Rlayer*TDT(I+1))/(delt*Delx2 + delt*kt2*Rlayer + Cp2*Delx2**2.d0*RhoS2*Rlayer)
CASE (FullyImplicitFirstOrder)
TDT(I) =( 2.d0*delt*Delx2*QSSFlux*Rlayer + Cp2*Delx2**2*RhoS2*Rlayer*TD(I) + 2.d0*delt*Delx2*TDT(I-1) + &
2.d0*delt*kt2*Rlayer*TDT(I+1))/(2.d0*delt*Delx2 + 2.d0*delt*kt2*Rlayer + Cp2*Delx2**2.d0*RhoS2*Rlayer)
END SELECT
IF ((TDT(I) > MaxSurfaceTempLimit) .OR. &
(TDT(I) < MinSurfaceTempLimit) ) THEN
TDT(I) = Max(MinSurfaceTempLimit,Min(MaxSurfaceTempLimit,TDT(I))) ! +++++ Limit Check
! CALL CheckFDSurfaceTempLimits(I,TDT(I))
ENDIF
ELSE IF (.not. RLayerPresent .and. RLayer2Present) THEN ! R-layer second
! check for PCM layer before R layer
IndVarCol = 1 ! temperature
DepVarCol = 2 ! thermal conductivity
IF( SUM(MaterialFD(MatLay)%TempEnth(1:3,2)) > 0.0d0 &
.and. SUM(MaterialFD(MatLay2)%TempEnth(1:3,2)) < 0.0d0) THEN ! phase change material Layer1, Use TempEnth Data
Enth1Old = terpld(MaterialFD(MatLay)%numTempEnth,MaterialFD(MatLay)%TempEnth,TD(I),IndVarCol,DepVarCol)
Enth1New = terpld(MaterialFD(MatLay)%numTempEnth,MaterialFD(MatLay)%TempEnth,TDT(I),IndVarCol,DepVarCol)
EnthNew(I) = Enth1New ! This node really doesn't have an enthalpy, this gives it a value
IF (ABS(Enth1New-Enth1Old) <= smalldiff .or. ABS(TDT(I)-TD(I)) <= smalldiff) THEN
Cp1=Cpo1
ELSE
Cp1=Max(Cpo1,(Enth1New -Enth1Old)/(TDT(I)-TD(I)))
END IF
ENDIF
SELECT CASE (CondFDSchemeType)
CASE (CrankNicholsonSecondOrder)
TDT(I)=(2.d0*delt*Delx1*QSSFlux*Rlayer2 - delt*Delx1*TD(I) - &
delt*kt1*Rlayer2*TD(I) + Cp1*Delx1**2*RhoS1*Rlayer2*TD(I) + delt*kt1*Rlayer2*TD(I-1) + &
delt*Delx1*TD(I+1) + delt*kt1*Rlayer2*TDT(I-1) + delt*Delx1*TDT(I+1))/ &
(delt*Delx1 + delt*kt1*Rlayer2 + Cp1*Delx1**2*RhoS1*Rlayer2)
CASE (FullyImplicitFirstOrder)
TDT(I)=(2.d0*delt*Delx1*QSSFlux*Rlayer2 + Cp1*Delx1**2*RhoS1*Rlayer2*TD(I) + 2.d0*delt*kt1*Rlayer2*TDT(I-1) + &
2.d0*delt*Delx1*TDT(I+1))/ (2.d0*delt*Delx1 + 2.d0*delt*kt1*Rlayer2 + Cp1*Delx1**2*RhoS1*Rlayer2)
END SELECT
IF ((TDT(I) > MaxSurfaceTempLimit) .OR. &
(TDT(I) < MinSurfaceTempLimit) ) THEN
TDT(I) = Max(MinSurfaceTempLimit,Min(MaxSurfaceTempLimit,TDT(I))) ! +++++ Limit Check
! CALL CheckFDSurfaceTempLimits(I,TDT(I))
ENDIF
ELSE ! Regular or Phase Change on both sides of interface
! Consider the various PCM material location cases
Cp1 = Cpo1 ! Will be changed if PCM
Cp2 = Cpo2 ! Will be changed if PCM
IndVarCol = 1 ! temperature
DepVarCol = 2 ! thermal conductivity
IF( Sum(MaterialFD(MatLay)%TempEnth(1:3,2)) > 0.0d0 &
.and. Sum(MaterialFD(MatLay2)%TempEnth(1:3,2)) > 0.0d0) THEN ! phase change material both layers, Use TempEnth Data
Enth1Old = terpld(MaterialFD(MatLay)%numTempEnth,MaterialFD(MatLay)%TempEnth,TD(I),IndVarCol,DepVarCol)
Enth2Old = terpld(MaterialFD(MatLay2)%numTempEnth,MaterialFD(MatLay2)%TempEnth,TD(I),IndVarCol,DepVarCol)
Enth1New = terpld(MaterialFD(MatLay)%numTempEnth,MaterialFD(MatLay)%TempEnth,TDT(I),IndVarCol,DepVarCol)
Enth2New = terpld(MaterialFD(MatLay2)%numTempEnth,MaterialFD(MatLay2)%TempEnth,TDT(I),IndVarCol,DepVarCol)
EnthNew(I) = Enth1New ! This node really doesn't have an enthalpy, this gives it a value
IF (ABS(Enth1New-Enth1Old) <= smalldiff .or. ABS(TDT(I)-TD(I)) <= smalldiff) THEN
Cp1 = Cpo1
ELSE
Cp1 = MAX(Cpo1,(Enth1New -Enth1Old)/(TDT(I)-TD(I)))
END IF
IF (ABS(Enth2New-Enth2Old) <= smalldiff .or. ABS(TDT(I)-TD(I)) <= smalldiff) THEN
Cp2 = Cpo2
ELSE
Cp2 = Max(Cpo2,(Enth2New -Enth2Old)/(TDT(I)-TD(I)))
END IF
ELSE IF( SUM(MaterialFD(MatLay)%TempEnth(1:3,2)) > 0.0d0 &
.and. SUM(MaterialFD(MatLay2)%TempEnth(1:3,2)) < 0.0d0) THEN ! phase change material Layer1, Use TempEnth Data
Enth1Old = terpld(MaterialFD(MatLay)%numTempEnth,MaterialFD(MatLay)%TempEnth,TD(I),IndVarCol,DepVarCol)
Enth1New = terpld(MaterialFD(MatLay)%numTempEnth,MaterialFD(MatLay)%TempEnth,TDT(I),IndVarCol,DepVarCol)
EnthNew(I) = Enth1New ! This node really doesn't have an enthalpy, this gives it a value
IF (ABS(Enth1New-Enth1Old) <= smalldiff .or. ABS(TDT(I)-TD(I)) <= smalldiff) THEN
Cp1 = Cpo1
ELSE
Cp1 = Max(Cpo1,(Enth1New -Enth1Old)/(TDT(I)-TD(I)))
END IF
Cp2 = Cpo2
ELSE IF( SUM(MaterialFD(MatLay)%TempEnth(1:3,2)) < 0.0d0 &
.and. SUM(MaterialFD(MatLay2)%TempEnth(1:3,2)) > 0.0d0) THEN ! phase change material Layer2, Use TempEnth Data
Enth2Old = terpld(MaterialFD(MatLay2)%numTempEnth,MaterialFD(MatLay2)%TempEnth,TD(I),IndVarCol,DepVarCol)
Enth2New = terpld(MaterialFD(MatLay2)%numTempEnth,MaterialFD(MatLay2)%TempEnth,TDT(I),IndVarCol,DepVarCol)
EnthNew(I) = Enth2New ! This node really doesn't have an enthalpy, this gives it a value
IF (ABS(Enth2New-Enth2Old) <= smalldiff .or. ABS(TDT(I)-TD(I)) <= smalldiff) THEN
Cp2 = Cpo2
ELSE
Cp2 = MAX(Cpo2,(Enth2New -Enth2Old)/(TDT(I)-TD(I)))
END IF
Cp1 = Cpo1
END If ! Phase change material check
SELECT CASE (CondFDSchemeType)
CASE (CrankNicholsonSecondOrder)
! Regular Internal Interface Node with Source/sink using Adams Moulton second order
TDT(I) = (2.d0*delt*Delx1*Delx2*QSSFlux - delt*Delx2*kt1*TD(I) - delt*Delx1*kt2*TD(I) + &
Cp1*Delx1**2.d0*Delx2*RhoS1*TD(I) + Cp2*Delx1*Delx2**2.d0*RhoS2*TD(I) + delt*Delx2*kt1*TD(I-1) + &
delt*Delx1*kt2*TD(I+1) + delt*Delx2*kt1*TDT(I-1) + delt*Delx1*kt2*TDT(I+1))/ &
(delt*Delx2*kt1 + delt*Delx1*kt2 + Cp1*Delx1**2.d0*Delx2*RhoS1 + Cp2*Delx1*Delx2**2.d0*RhoS2)
CASE (FullyImplicitFirstOrder)
! first order adams moulton
TDT(I) = (2.d0*delt*Delx1*Delx2*QSSFlux + &
Cp1*Delx1**2.d0*Delx2*RhoS1*TD(I) + Cp2*Delx1*Delx2**2.d0*RhoS2*TD(I) + &
2.d0*delt*Delx2*kt1*TDT(I-1) + 2.d0*delt*Delx1*kt2*TDT(I+1))/ &
(2.d0*delt*Delx2*kt1 + 2.d0*delt*Delx1*kt2 + Cp1*Delx1**2.d0*Delx2*RhoS1 + Cp2*Delx1*Delx2**2.d0*RhoS2)
END SELECT
IF ((TDT(I) > MaxSurfaceTempLimit) .OR. &
(TDT(I) < MinSurfaceTempLimit) ) THEN
TDT(I) = Max(MinSurfaceTempLimit,Min(MaxSurfaceTempLimit,TDT(I))) ! +++++ Limit Check
! CALL CheckFDSurfaceTempLimits(I,TDT(I))
ENDIF
IF (Construct(ConstrNum)%SourceSinkPresent .and. Lay == Construct(ConstrNum)%SourceAfterLayer) Then
TcondFDSourceNode(Surf)= TDT(I) ! transfer node temp to Radiant System
TempSource(Surf) = TDT(I) ! Transfer node temp to DataHeatBalSurface module.
ENDIF
!+++++++++++++++++++++++++++++++
END IF ! end of R-layer and Regular check
END IF !End of the CondFD if block
RETURN
END SUBROUTINE IntInterfaceNodeEqns