Engineeringreference Energy Plus
Engineeringreference Energy Plus
EnergyPlus
Engineering Reference
COPYRIGHT © 1996-2011 The Board of Trustees of the University of Illinois and the Regents of the University of California through
the Ernest Orlando Lawrence Berkeley National Laboratory.
All Rights Reserved. No part of this material may be reproduced or transmitted in any form or by any means without the prior written
permission of the University of Illinois or the Ernest Orlando Lawrence Berkeley National Laboratory.
EnergyPlus is a Trademark of the US Department of Energy.
TABLE OF CONTENTS
Overview ..........................................................................................................................................1
Moisture Prediction.................................................................................................. 17
References.............................................................................................................. 32
4/1/13 i
TABLE OF CONTENTS
Conduction Through The Walls ..................................................................................... 34
References.............................................................................................................. 42
References.............................................................................................................. 47
Overview ................................................................................................................. 47
References.............................................................................................................. 52
Overview ................................................................................................................. 52
EMPD Nomenclature............................................................................................... 55
References.............................................................................................................. 56
4/1/13 ii
TABLE OF CONTENTS
References.............................................................................................................. 59
References.............................................................................................................. 72
Introduction ............................................................................................................. 99
References............................................................................................................ 105
4/1/13 iii
TABLE OF CONTENTS
Fortran Variable Descriptions ................................................................................ 111
References............................................................................................................ 113
References............................................................................................................ 123
References............................................................................................................ 132
4/1/13 iv
TABLE OF CONTENTS
Perez Direct/Diffuse Splitting Model ...................................................................... 139
References............................................................................................................ 139
References............................................................................................................ 164
Diffuse Reflection of Beam Solar and Sky Solar Radiation .......................................... 167
Rays...................................................................................................................... 168
4/1/13 v
TABLE OF CONTENTS
Sky Solar Radiation Diffusely Reflected from Obstructions.................................... 169
4/1/13 vi
TABLE OF CONTENTS
References............................................................................................................ 196
References............................................................................................................ 200
4/1/13 vii
TABLE OF CONTENTS
Edge-Of-Glass Effects........................................................................................... 293
Heat Balance Equations for Shading Device and Adjacent Glass.......................... 312
References............................................................................................................ 341
Infiltration............................................................................................................... 343
Reference.............................................................................................................. 349
4/1/13 viii
TABLE OF CONTENTS
Temperature Difference Controlled Air Exchange.................................................. 349
References............................................................................................................ 353
4/1/13 ix
TABLE OF CONTENTS
Zone Supply Air Path............................................................................................. 363
In addition, to avoid the need for overly complex solver routines, there are
some restrictions on the placement of pumps within a particular half-loop.
There are two general types of pumps, loop pumps and branch pumps. A
pump that is the first component on the first branch (between A and B) is
termed a “loop pump” while any pump in the parallel section (between Ci
and Di) is termed a “branch pump”. The simplest and most common
arrangement is to have one loop pump on the supply side inlet. In plant
demand half-loops pumps can be placed only in the inlet branch. This will
allow simulation of primary-secondary systems. For more information on
pumps and pump placement rules, see the section on
PipingSystem:Underground Simulation ................................................................. 368
The pump is quite simply the component that drives the flow (also see
PipingSystem:Underground Simulation ................................................................. 374
4/1/13 x
TABLE OF CONTENTS
Heat Recovery Loop Systems ............................................................................... 390
System Design Flow Rate and Load Summation and Adjustment ......................... 417
Coil:Cooling:Water................................................................................................. 425
4/1/13 xi
TABLE OF CONTENTS
Coil:Cooling:Water:DetailedGeometry Sizing......................................................... 428
4/1/13 xii
TABLE OF CONTENTS
Zone Outdoor Air Design Data .................................................................................... 456
References............................................................................................................ 458
Demand Limiting..........................................................................................................................459
Impacts of Supply Air Constant Volume Fan Control on Load: Cycling vs.
Continuous ............................................................................................................ 529
Airflow Calculation Procedure using A Supply Variable Air Volume Fan ................ 530
References............................................................................................................ 531
4/1/13 xiii
TABLE OF CONTENTS
Air System Distribution Terminals................................................................................ 533
Variable Air Volume Single Duct Reheat and No Reheat Air Terminals ................. 534
Variable Air Volume Heating and Cooling Single Duct Reheat and
NoReheat Air Terminal .......................................................................................... 535
Constant Volume Single Duct Four Pipe Induction Air Terminal ............................ 538
Fan Powered Induction Series and Parallel Single Duct Reheat Air
Terminal ................................................................................................................ 540
Variable Air Volume Fan Powered Single Duct Air Terminal.................................. 543
Dual Duct Dedicated Outside Air Terminal with VAV Cooling ................................ 551
4/1/13 xiv
TABLE OF CONTENTS
Electric Chiller Model Based on Condenser Leaving Temperature ........................ 597
Coils............................................................................................................................ 610
4/1/13 xv
TABLE OF CONTENTS
Steam-Based Air Heating Coil ............................................................................... 733
Variable Speed Water to Air Heat Pump (Heating & Cooling)................................ 740
Outdoor Air Damper Controller for Zone Energy Recovery Ventilator .................... 765
One and Two Speed Cooling Towers and Evaporative Fluid Coolers.................... 776
4/1/13 xvi
TABLE OF CONTENTS
References............................................................................................................ 806
References............................................................................................................ 839
4/1/13 xvii
TABLE OF CONTENTS
Variable Refrigerant Flow Heat Pump Model ......................................................... 898
Air System Air-To-Air Sensible and Latent Effectiveness Heat Exchanger ............ 924
Air System Air-To-Air Balanced Flow Desiccant Heat Exchanger .......................... 937
References............................................................................................................ 967
References............................................................................................................ 975
Pipes........................................................................................................................... 975
Pumps......................................................................................................................... 984
4/1/13 xviii
TABLE OF CONTENTS
Summary of Pump Rules....................................................................................... 984
References............................................................................................................ 992
ChillerHeaterPerformance:Electric:EIR........................................................................ 996
References.......................................................................................................... 1003
4/1/13 xix
TABLE OF CONTENTS
High Temperature Radiant Heater Model ............................................................ 1028
References.......................................................................................................... 1088
4/1/13 xx
TABLE OF CONTENTS
Multizone Heating Average.................................................................................. 1097
4/1/13 xxi
TABLE OF CONTENTS
High Temperature Turn On.................................................................................. 1139
References.......................................................................................................... 1166
4/1/13 xxii
TABLE OF CONTENTS
Connected Water Use Equipment ....................................................................... 1196
Earthtube............................................................................................................. 1242
4/1/13 xxiii
TABLE OF CONTENTS
Thermal Chimney Model ..................................................................................... 1248
Controls............................................................................................................... 1253
On-Site Generation....................................................................................................................1256
Micro-Cogenerator............................................................................................... 1291
4/1/13 xxiv
TABLE OF CONTENTS
Curves based on a single independent variable .................................................. 1307
Adjustments......................................................................................................... 1323
ExpressAsCashFlows.......................................................................................... 1329
ComputeTaxAndDepreciation.............................................................................. 1329
Special Modules/Reporting........................................................................................................1331
4/1/13 xxv
TABLE OF CONTENTS
Other Energy-Related Pollutants and Sources of Other Information .................... 1342
References.......................................................................................................... 1342
4/1/13 xxvi
Overview Document Overview
Overview
Document Overview
This document is organized to give you the best possible look into the EnergyPlus
calculations. First, the concepts of modeling in EnergyPlus are presented. These include
descriptions of the zone heat balance process, air loop/plant loop processes as well as other
important processes for the building simulation.
Discussions during the modeling process may reference specific “object names” as found in
the Input/Output Reference document.
The remainder of the document focuses on individual models.
The EnergyPlus program is a collection of many program modules that work together to
calculate the energy required for heating and cooling a building using a variety of systems
and energy sources. It does this by simulating the building and associated energy systems
when they are exposed to different environmental and operating conditions. The core of the
simulation is a model of the building that is based on fundamental heat balance principles.
Since it is relatively meaningless to state: “based on fundamental heat balance principles”,
the model will be described in greater detail in later sections of this document in concert with
the FORTRAN code which is used to describe the model. It turns out that the model itself is
relatively simple compared with the data organization and control that is needed to simulate
the great many combinations of system types, primary energy plant arrangements,
schedules, and environments. The next section shows this overall organization in schematic
form. Later sections will expand on the details within the blocks of the schematic.
Simulation Manager
The simulation manager of EnergyPlus is contained in a single module. The main subroutine
is shown below. Flow within the entire program is managed using a series of flags. These
paired flags, in order (from the highest to the lowest) are:
4/1/13 1
Overview Simulation Manager
BeginSimulationFlag EndSimulationFlag
BeginEnvironmentFlag EndEnvironmentFlag(one to many days)
BeginDayFlag EndDayFlag
BeginHourFlag EndHourFlag
BeginTimeStepFlag EndTimeStepFlag
There is also a WarmupFlag to signal that the program is in warmup state. The operation of
these flags can be seen in the following subroutine. The advantage of using the flag system
is that any subroutine throughout the code can determine the exact state of the simulation by
checking the status of the flags.
4/1/13 2
Overview Warmup Convergence
Warmup Convergence
Since everything in EnergyPlus is based on the foundation of the loads simulation, it stands
to reason that any inaccuracies in the loads calculation will result in inaccuracies of similar or
larger magnitude in the HVAC calculations. In the presumably limited cases where
convergence was not truly achieved before the actual simulation began, it is unknown how
much error would be introduced into the results. While simulations that last longer (annual vs.
design day) will hopefully have any initial condition problems balanced by the shear number
of days in the simulation, shorter simulations—particularly those used for sizing—could result
in relatively large errors. The simulation results could be unreliable and inaccurate when
steady periodic conditions are not achieved. Therefore, it is important to properly determine
when there is enough temperature and flux history terms to start an EnergyPlus simulation
since this has a potential economic and energy impact on buildings that use EnergyPlus in
design.
4/1/13 3
Overview Warmup Convergence
EnergyPlus determines warmup convergence in the following manner as shown in the Figure
2 below. The process of the convergence checks begins by tracking four parameters such
including the maximum zone air temperature, the minimum zone air temperature, the
maximum heating load, and the maximum cooling load for individual zone. It is note that
these convergence checks are only in effective in simulations with at least one zone since the
criteria is solely based on the maximum and minimum values obtained from an individual
zone. Differences in these parameters between two consecutive days are then compared
with the convergence tolerance values at the end of the day during the warmup period. For
example, the maximum and minimum air temperature and the percentage difference of zone
load for each zone at 9:00AM during the second to last warmup is compared to the values at
9:00AM last warmup day as follows:
qh qh , prev
qtol (3)
qh
qc qc , prev
qtol (4)
qc
where Tmax,prev is the maximum zone temperature of previous day, Tmax is the maximum
zone temperature of current day, Ttol is the value of temperature tolerance, qh,prev is the
maximum heating load of previous day, qh, is the maximum heating load of current day, qtol
is the value of load tolerance, qc,prev is the maximum cooling load of previous day, and qc, is
the maximum cooling load of current day.
Note that a minimum load of 100W is used to establish a fraction for the maximum loads
when they are less than the minimum. This is done to avoid a false negative indication for the
percentage load difference that may appear when zonal loads are very small. The
convergence checks are repeated until passed for all zones. EnergyPlus assumes that the
warmup period has been reached steady-periodic when these four parameters are within
tolerance. Finally, temperature and load differences between the last two warmup days for
individual zone at each time step in the last warmup day are reported so that users can easily
track whether or not the warmup period has converged. The input parameters and output
related to the warmup period are discussed in the Input-Output Reference.
4/1/13 4
Overview Warmup Convergence
4/1/13 5
Integrated Solution Manager Warmup Convergence
EnergyPlus is an integrated simulation. This means that all three of the major parts, building,
system, and plant, must be solved simultaneously. In programs with sequential simulation,
such as BLAST or DOE-2, the building zones, air handling systems, and central plant
equipment are simulated sequentially with no feedback from one to the other. The sequential
solution begins with a zone heat balance that updates the zone conditions and determines
the heating/cooling loads at all time steps. This information is fed to the air handling
simulation to determine the system response; but that response does not affect zone
conditions. Similarly, the system information is passed to the plant simulation without
feedback. This simulation technique works well when the system response is a well-defined
function of the air temperature of the conditioned space. For a cooling situation, a typical
supply and demand situation is shown schematically in the Figure 3. Here, the operating
point is at the intersection of the supply and demand curves.
However, in most situations the system capacity is dependent on outside conditions and/or
other parameters of the conditioned space. The simple supply and demand situation above
becomes a more complex relationship and the system curve is not fixed. The solution should
move up and down the demand curve. This doesn’t happen in sequential simulation methods
and the lack of feedback from the system to the building can lead to nonphysical results. For
example, if the system provides too much cooling to a conditioned space the excess is
reported by the program as "overcooling". Other categories of unmatched loads exist and are
similarly reported by the program. While this kind of reporting enables the affected system or
plant components to be properly sized, the system designer would, in most cases, prefer to
see the actual change in zone temperature. The same mismatches can occur between the
system and plant simulations when they are simulated sequentially.
To obtain a simulation that is physically realistic, the elements have to be linked in a
simultaneous solution scheme. The entire integrated program can be represented as a series
of functional elements connected by fluid loops as shown in Figure “Schematic of
Simultaneous Solution Scheme”. In EnergyPlus all the elements are integrated and controlled
by the Integrated Solution Manager. The loops are divided into supply and demand sides,
and the solution scheme generally relies on successive substitution iteration to reconcile
supply and demand using the Gauss-Seidell philosophy of continuous updating.
4/1/13 6
Integrated Solution Manager Basis for the Zone and Air System Integration
In the sections which follow, the various individual functions of the integrated solution will be
described.
The basis for the zone and air system integration is to formulate energy and moisture
balances for the zone air and solve the resulting ordinary differential equations using a
predictor-corrector approach. The formulation of the solution scheme starts with a heat
balance on the zone air.
N surfaces
dTz N sl N zones
Cz Qi hi Ai Tsi Tz m C T i p zi Tz m inf C p T Tz Q sys (5)
dt i 1 i 1 i 1
where:
N sl
Q
i 1
i = sum of the convective internal loads
N surfaces
i 1
hi Ai Tsi Tz = convective heat transfer from the zone surfaces
m C T
i 1
i p zi Tz = heat transfer due to interzone air mixing
4/1/13 7
Integrated Solution Manager Basis for the Zone and Air System Integration
N sl N surfaces N zones
Q sys Q i hi Ai Tsi Tz m C T i p zi Tz m inf C p T Tz (6)
i 1 i 1 i 1
Air systems provide hot or cold air to the zones to meet heating or cooling loads. The system
energy provided to the zone, Q sys , can thus be formulated from the difference between the
supply air enthalpy and the enthalpy of the air leaving the zone as in Equation (7):
This equation assumes that the zone supply air mass flow rate is exactly equal to the sum of
the air flow rates leaving the zone through the system return air plenum and being exhausted
directly from the zone. Both air streams exit the zone at the zone mean air temperature. The
result of substituting Equation (7) for Q sys in the heat balance Equation (5) is shown in
Equation (8):
N sl N surfaces N zones
Cz z Q i
dT
dt i 1
i 1
hi Ai Tsi Tz m C T
i 1
i p zi Tz
(8)
m inf C p T Tz m sys C p Tsup Tz
The sum of zone loads and air system output now equals the change in energy stored in the
zone. Typically, the capacitance Cz would be that of the zone air only. However, thermal
masses assumed to be in equilibrium with the zone air could be included in this term.
EnergyPlus provides three different solution algorithms to solve the zone air energy and
moisture balance equations. These are defined in the Algorithm field in the
ZoneAirHeatBalanceAlgorithm object: 3rdOrderBackwardDifference, EulerMethod and
AnalyticalSolution. The first two methods to solve Equation (8) use the finite difference
approximation while the third uses an analytical solution. A short description is given below.
In order to calculate the derivative term with respect to time, a finite difference approximation
may be used, such as:
t Tz t Tz t t t
dT 1
(9)
dt
The use of numerical integration in a long time simulation is a cause for some concern due to
the potential build-up of truncation error over many time steps. In this case, the finite
difference approximation is of low order that further aggravates the problem. However, the
cyclic nature of building energy simulations should cause truncation errors to cancel over
each daily cycle so that no net accumulation of error occurs, even over many days of
simulation (Walton, 1990). The Euler formula, Equation (9), was employed in Equation (8) to
replace the derivative term. All the terms containing the zone mean air temperature were then
grouped on the left hand side of the equation. Since the remaining terms are not known at
the current time, they were lagged by one time step and collected on the right hand side. This
manipulation resulted in Equation (10), the formula for updating the zone mean air
temperature:
4/1/13 8
Integrated Solution Manager Basis for the Zone and Air System Integration
One final rearrangement was to move the lagged temperature in the derivative approximation
to the right side of the equation. The explicit appearance of the zone air temperature was
thus eliminated from one side of the equation. An energy balance equation that includes the
effects of zone capacitance was then obtained by dividing both sides by the coefficient of Tz:
t t
N sl
m C T Tz N surfaces N zones
Q t
z
C
t
h AT m iC pTzi m inf C pT
t
i sys p supply i i si
Tzt
i 1 i 1 i 1
(11)
Cz
N surfaces N zones
hi Ai m i C p m inf C p m sysC p
t i 1 i 1
Equation (11) could be used to estimate zone air temperatures, and is defined as the
EulerMethod, one of the three solution algorithms provided in the
ZoneAirHeatBalanceAlgorithm object. However, it can severely limit the time step size under
some conditions. To improve on this, higher order expressions for the first derivative, with
corresponding higher-order truncation errors, were developed. The goal of this approach was
to allow for the use of larger time steps in the simulation than would be possible using the first
order Euler form, without experiencing instabilities. Approximations from second through fifth
order were tried as reported by Taylor, et al. (1990) with the conclusion that the third order
finite difference approximation, shown below, gave the best results:
1 11
t Tzt 3Tzt t Tzt 2t Tzt 3t O t 3
dTz 3 1
(12)
dt t 6 2 3
When this form for the derivative is used, equation (10) changes to:
N surfaces
1 11
N sl N zones
C z t Tzt 3Tzt t Tzt 2t Tzt 3 t Qi
3 1
6 2 3 i 1
i 1
hi Ai Tsi Tz m C T
i 1
i p zi Tz
(13)
m inf C p T Tz m sys C p Tsup Tz
and the zone temperature update equation becomes:
Cz 3T t t 3 T t 2 t 1 T t 3 t
N sl N surfaces N zones
Q i hi AT
i si
m C T i p zi
m inf C pT m sys C pTsupply z
t 2
z
3
z
Tzt i 1 i 1 i 1
(14)
11 Cz
N surfaces N zones
6 t
h A m C
i i p
m inf C p m sys C
i 1 i 1
This is the form historically used in EnergyPlus and is the current default referred to as
3rdOrderBackwardDifference in the ZoneAirHeatBalanceAlgorithm object. This algorithm
requires zone air temperatures at three previous time steps and uses constant temperature
coefficients. The assumption is that three previous time steps lengths are the same.
4/1/13 9
Integrated Solution Manager Basis for the Zone and Air System Integration
rd
The AnalyticalSolution algorithm is an integration approach. While the 3 order finite
difference approximation provides stability without requiring a prohibitively small time step,
the method still has truncation errors and requires a fixed time step length for the previous
three simulation time steps. Therefore, different time step lengths for the previous three
simulation time steps may make the temperature coefficients invalid.
The AnalyticalSolution algorithm provides a possible way to obtain solutions without
truncation errors and independent of time step length. In addition, the algorithm only requires
the zone air temperature for one previous time step, instead of three previous time steps as
required by the 3rdOrderBackwardDifference algorithm. The integrated (analytical) solution
for Eq. (4) may be expressed as follows:
N sl
N surfaces N zones
Q i h AT
i i si m iC pTzi m inf C pT m sys C pTsup
Tzt Tzt t i 1 i 1 i 1
N surfaces N zones
i 1
hi Ai m i C p m inf C p m sys C p
i 1
N surfaces N zones
hi Ai m i C p m inf C p m sys C p
*exp i 1 i 1
t
Cz
N sl N surfaces N zones
Qi i si
hi AT m C T i p zi m inf C pT m sys C pTsup
i 1 i 1
N surfaces
i 1
N zones
i 1
hi Ai m C
i 1
i p m inf C p m sys C p
(15)
Since the load on the zone drives the entire process, that load is used as a starting point to
give a demand to the air system. Then a simulation of the air system provides the actual
supply capability and the zone temperature is adjusted if necessary. This process in
EnergyPlus is referred to as a Predictor/Corrector process. It is summarized below.
Code Reference: the ZoneTempPredictorCorrector module performs the
calculations.
4/1/13 10
Integrated Solution Manager Summary of Predictor-Corrector Procedure
Previously, the formulation of a new heat balance equation with an unsteady zone
capacitance term was discussed Equation (7). In this equation the updated zone temperature
was calculated by removing its explicit dependence from the right hand side and lagging, by
one time step, the unknown terms on that side. However, the right hand side still contains
implicit dependencies on the zone temperature through the air system control logic; the need
for heating or cooling in the zones, is based on zone temperature. In real buildings the control
system consists of one or more sensing units in the zone, such as a wall thermostat that
samples the air temperature and sends signals to a control unit. The controller looks at the
difference between the actual zone temperature and the desired temperature to ascertain if
heating or cooling is required and then sends appropriate signals to the air system
components to drive the zone temperature closer to the desired value.
Although many control systems use only the zone air temperature to control the air system,
most modern energy management systems consider many other variables, such as outside
environment conditions. Simulating such controllers would seem to be relatively
straightforward in a simulation especially since some of the more complex control problems,
such as managing duct pressures and flow rates, are not always modeled. However, real
controllers have an advantage because they can sample zone conditions, and thus update air
system response, on a time scale much shorter than any characteristic time of the air system
or zone. Thus the feedback between zone and air system usually results in steady or, at
worst, slowly oscillating zone conditions and air system operation unless the air system is
grossly oversized. On the other hand, the numerical model is only able to sample zone
conditions at discrete time intervals. In the interest of minimizing computation time, these
intervals need to be as long as possible. Frequently, they are of the order of, or longer than,
the characteristic times of the air system and zones, except in the case of small air system
capacity in relation to zone capacitance. This situation has the potential for unstable feedback
between the zone and air system, resulting in an oscillatory or diverging solution.
Prior to implementing the new heat balance method (3rdOrderBackwardDifference) in
IBLAST, several air system control strategies were considered. The primary objective was
selection of a control method that would be numerically stable over a reasonable range of
conditions, realistic from the standpoint of looking and operating like an actual air system
controller, and flexible enough to be applied to all current and projected systems. The method
actually implemented in IBLAST, and later EnergyPlus, took advantage of the computational
model's "knowledge" of how much energy enters or leaves the zone as a function of zone air
temperature i.e., the zone load. The real controller, on the other hand, does not have this
information. The net zone load is given by Equation (16):
4/1/13 11
Integrated Solution Manager Air System Control
N sl N surfaces N zones
Q load Q i hi Ai Tsi Tz m C T i p zi Tz m inf C p T Tz (16)
i 1 i 1 i 1
This is Equation (6) without the term due to the air system. In addition, T z is now the desired
zone temperature as defined by the control system setpoints that must be specified for each
zone. An assumption was made that if the air system has sufficient capacity (based on the
desired zone air temperature) to meet the zone conditioning requirements (i.e. Q = Q load )
sys
at the desired zone air temperature then those requirements will be met. On the other hand, if
the air system cannot provide enough conditioning to the zone to maintain the desired
temperature, then the air system provides its maximum output to the zone and the zone air
temperature is allowed to "float." Equation (16) was used to calculate the air system output
required to maintain the desired zone air temperature; the actual zone temperature update
was accomplished using Equation (11). This method was called predictive system energy
balance. It has many characteristics of a predictor-corrector method since the air system
response is first approximated based on a predicted zone temperature and then the actual
change in zone temperature is determined from that air system response. The predictive air
system energy balance method required that the system controls on air mass flow rate,
supply air temperature, etc., be formulated as a function of the zone air temperature.
However, this was not a serious drawback. The first example considered was a single zone
draw through air system. Typically, such systems have a cooling coil and heating coil in
series, and constant air volume flow rate. Single zone draw through systems run at maximum
capacity when turned on; so the only way to regulate net air system output and keep the zone
air temperature within the desired range is to turn the air system on and off. A simplified
schematic of this system type is shown in Figure 5. Simplified Single Zone Draw Through Air
System.
OUTSIDE
AIR
MIXING
C/C H/C
BOX
RELIEF
AIR CONSTANT
VOLUME FAN
RETURN AIR
ZONE
The amount of heating or cooling provided by the air system in relation to the desired zone air
temperature is given by:
where is the fraction of the time step that the air system is turned on and varies between 0
and 1. The supply air temperature is also implicitly limited by the effectiveness of the coils
and the operating parameters of the central plant components. These interactions are
discussed later.
A far more complex, though again simplified, air system is the variable air volume (VAV)
system, shown in Figure 6. Simplified Variable Volume Air System. In VAV systems, the
supply air temperature, as well as the supply air volume, are continuous functions of zone air
temperature. As shown in Figure 7. Idealized Variable Volume System Operation., when the
4/1/13 12
Integrated Solution Manager Air System Control
zone air temperature is between Tcl and Tcu, cooling is required and the air system varies
the supply air flow rate while maintaining a constant supply air temperature. When the zone
air temperature is between Thl and Thu, heating is required and air is supplied at a constant
minimum flow rate while the supply air temperature is varied.
OUTSIDE VARIABLE
AIR VOLUME FAN
MIXING
C/C
BOX
RELIEF
AIR
DAMPER
RETURN AIR
ZONE H/C
The next figure (Idealized variable volume system operation) shows idealized behavior of a
VAV system; in practice, the air flow rate and temperature are not exact linear functions of
zone air temperature.
FRACTION
COLD
TEMPERATURE
DECK
MIN.
T hl T hu T cl T cu
ZONE TEMPERATURE
Figure 7. Idealized Variable Volume System Operation.
As long as a VAV system has sufficient capacity, the zone air temperatures can be expected
to vary within the limits defining the range of operation of the air damper, when cooling, or the
throttling range of the reheat coil, when the air system is heating. This means that the desired
zone air temperature, used to predict the air system response, is variable and must be
calculated in order to determine the air system output. For the purposes of this calculation,
the following definitions were found useful:
N sl N surfaces N zones
Q 0 Q i i si
hi AT m C T i p zi m inf C p T (18)
i 1 i 1 i 1
N surfaces N zones
Q slope hi Ai m C i p m inf C p (19)
i 1 i 1
4/1/13 13
Integrated Solution Manager Air System Control
Equations (18) and (19) are derived, respectively, from the numerator and denominator of
Equation (14) but with the system related terms omitted. Also excluded from these
expressions are the effects of zone capacitance.
When a zone requires cooling, the VAV system is designed to provide air to that zone at a
constant supply air temperature. The amount of cooling is matched to the load by dampers in
the supply air duct that vary the air volume flow rate of being supplied to the zone. Assuming
that the volume flow rate varies linearly with zone air temperature, the volume flow rate of
supply air normalized to the maximum flow rate, or supply air fraction, is given by:
Tz Tc , lower
c c , min 1 c , min ;c , min c 1.0 (20)
Tc , upper Tc , lower
Normally, the minimum supply air fraction c,min must be greater than zero to ensure a
supply of fresh air sufficient to eliminate contaminants from the zone.
Conversely, when heating is required in a zone, the VAV system becomes a constant volume
flow rate system with a variable supply air temperature. The dampers are set to provide air to
the zone at the minimum supply air fraction. Throttling the hot water supply to the reheat coil,
which effectively alters the coil’s heating capacity, modulates the supply air temperature.
Again, assuming the heat energy output varies linearly with zone air temperature and
normalizing with respect to the maximum coil output gives the following result:
Th , upper Tz
h ;0 h 1.0 (21)
Th , upper Th , lower
Observe that when h is equal to zero, the zone is supplied with air at the cooling coil outlet
temperature at the minimum air fraction. Because the control strategies of the VAV system
are different whether the air system is heating or cooling, two equations are necessary to
describe the air system output in terms of h and c. These expressions are as shown in
Equations (22) and (23):
Equation (22) is valid for zone air temperatures below Th,upper, while Equation (23) is valid
for all temperatures above this value. Equating the system output to the zone load, as given
by Equation (16), the definitions of c and h were then used to develop expressions for the
predicted zone air temperature in the cases of heating and cooling:
Q h / c ,maxTh,upper C p VminTc / c
Tz , pred ,heat Q 0
Q h / c ,max
(24)
Th,upper Th,lower
C p Vmin Q slope
Th,upper Th,lower
B1 B12 B2
Tz , pred ,cool (25)
2
4/1/13 14
Integrated Solution Manager Air System Control
where,
c ,min C2
B1 Tc / c Tc ,lower (26)
C1
C
B2 4 3 Tc / c c ,min Tc ,lower (27)
C1 C1
and,
1 c ,min
C1 (28)
Tc ,upper Tc ,lower
Q slope
C2 (29)
C V
p max
Q0
C3
Cp Vmax
(30)
Once the predicted zone air temperature has been calculated from Equations (24) and (25),
the air system response may be determined. When a zone requires cooling the system
supply air temperature is constant at the cooling coil outlet temperature and the volume flow
rate is given by:
where the supply air fraction c is computed from Equation (20). When heating is required by
the zone, the air system provides air at the minimum volume flow rate and at a temperature
given by:
h Q h / c ,max
Tsup ply Tc / c
C p Vmin
(32)
The reheat coil capacity fraction h is determined by using Equation (21). Once Equation (31)
or (32), has been used, the supply air flow rate and temperature are known. These values are
then used in Equation (11) to calculate the updated zone air temperature. The equations
describing VAV system operation may be solved without iteration if the cooling coil outlet
temperature is constant, i.e. if the coil has infinite capacity, and if the reheat coil capacity
varies linearly with zone air temperature. This is not the case, either in practice or in
simulations, when realistic coil models are used. Therefore, an iteration scheme was
developed that solved these equations simultaneously with the coil performance models.
4/1/13 15
Integrated Solution Manager Moisture Predictor-Corrector
Moisture Predictor-Corrector
The transient air mass balance equation for the change in the zone humidity ratio = sum of
internal scheduled latent loads + infiltration + system + multizone airflows + convection to the
zone surfaces may be expressed as follows:
N surfaces
dWz Nsl
m W
N zones
airVz CW kg masssched load Ai hmi airz Wsurfsi W
z
t
i zi Wzt
dt i 1 i 1 i 1
where
CW = humidity capacity multiplier (See the InputOutput Reference for additional
information on the object ZoneCapacitanceMultiplier:ResearchSpecial)
In the same manner as described above for zone air temperature (ref. Basis for the Zone and
Air System Integration), the solution algorithms provided in the ZoneAirHeatBalanceAlgorithm
object are also applied to zone air moisture calculations.
In order to calculate the derivative term with respect to time, the first order backward finite
difference method, defined as the EulerMethod in the ZoneAirHeatBalanceAlgorithm object,
may be used:
dW
t (Wzt Wzt t ) O ( t )
1
dt
The zone air humidity ratio update at the current time step using the EulerMethod may be
expressed as follows:
N surfaces
m W
N sl N zones
airVCW t
1
Wzt Wzt t kgmasssched load Ai hmi airz Wsurfsi Wzt i zi Wzt
i 1 i 1 i 1
To preserve the stability of the calculation of the zone humidity ratio, the third order
differential approximation, derived by a Taylor Series and used in the calculation of the next
time step’s zone air temperature, is also applied to the zone air humidity ratio calculations.
This algorithm is the default choice and is defined as 3rdOrderBackwardDifference in the
ZoneAirHeatBalanceAlgorithm object.
The third order derivative derived from a Taylor Series expansion is defined as:
dWz
6 Wz 3Wz
11 t t t
23 Wzt 2 t 13 Wzt 3 t
O t 3 .
dt t t
The coefficients of the approximated derivative are very close to the coefficients of the
analogous Adams-Bashforth algorithm. Then the approximated derivative is substituted into
the mass balance and the terms with the humidity ratio at past time steps are all put on the
right hand side of the equation. This third order derivative zone humidity ratio update
increases the number of previous time steps that are used in calculating the new zone
humidity ratio, and decreases the dependence on the most recent. The higher order
derivative approximations have the potential to allow the use of larger time steps by
smoothing transitions through sudden changes in zone operating conditions.
4/1/13 16
Integrated Solution Manager Moisture Predictor-Corrector
airVz CW 11 t
N surfaces N zones N sl N surfaces
t 6
Wz i 1
Ai hmi airz W
z
t
mW
i 1
i z
t
m inf W m sysW kgmasssched load
z
t
z
t
i 1
i 1
Ai hmi airz Wsurfsi
N zones
airVz CW 3 t 2 t 1 t 3 t
mW
i 1
i zi m inf W m sysWsup
t
t t
3Wz Wz
2
Wz
2
This gives us the basic air mass balance equation that will be solved two different ways, one
way for the predict step and one way for the correct step.
Since the third choice of solution algorithms uses an integration approach, defined as
AnalyticalSolution in the ZoneAirHeatBalanceAlgorithm object, it does not require any
approximations and has no truncation errors. The solutions in both prediction and correction
are provided below in detail.
Moisture Prediction
For the moisture prediction case the equation is solved for the anticipated system response
as shown below.
Since the program provides three solution algorithms, the moisture prediction from each
solution algorithm is given below.
EulerMethod
For this solution algorithm, the air mass balance for the predicted air system load or response
is:
ThirdOrderBackwardDifference
For this solution algorithm, the air mass balance for the predicted system load or response is
given below:
V C 11 Nsurfaces N zones
PredictedSystemLoad [kgWater / sec] air z W Ai hmi airz m i m inf Wzt
t 6 i 1 i 1
N sl N surfaces N zones
VC 3 1
kg masssched load Ai hmi airz Wsurfsi mW i zi m inf W air z W 3Wzt t Wzt 2 t Wzt 3 t
i 1 i 1 i 1 t 2 3
Then, using the following substitutions, the air mass balance equation becomes:
N surfaces N zones
A i 1
Ai hmi airz m m
i 1
i inf
4/1/13 17
Integrated Solution Manager Moisture Predictor-Corrector
N sl N surfaces N zones
B kg masssched load Ai hmi airz Wsurfsi mW
i zi m inf W
i 1 i 1 i 1
airVz CW
C
t
11
PredictedSystemLoad [kgWater / sec] * C A *WSetPo int
6
t t 3 t 2 t 1 t 3 t
B C * 3Wz 2 Wz Wz
3
AnalyticalSolution
For this solution algorithm, the air mass balance for the predicted air system load or response
is given below:
N surfaces N zones
PredictedSystemLoad [kgWater / sec] Ai hmi airz mW i zi m inf *
i 1 i 1
N surfaces N zones
i mi airz mW
A h i zi m inf
W t Wzt t *exp i 1 i 1
t *
setpoint airVz CW
1
N surfaces N zones
i mi airz m i m inf
A h
1 exp i 1 i 1
t
airVz CW
N sl N surfaces N zones
masssched load i mi airz surfsi mW
kg A h W i zi m inf W
i 1 i 1 i 1
At the prediction point in the simulation, the system air mass flows are not known; therefore,
the system response is approximated. The predicted air system moisture load is then used in
the system simulation to achieve the best results possible. The system simulation
components that have moisture control will try to meet this predicted moisture load. For
example, humidifiers will look for positive moisture loads and add moisture at the specified
rate to achieve the relative humidity setpoint. Likewise, dehumidification processes will try to
remove moisture at the specified negative predicted moisture load to meet the relative
humidity setpoint.
After the system simulation is completed the actual response from the air system is used in
the moisture correction of step, which is shown next.
Moisture Correction
For the correct step the expanded air mass balance equation is solved for the final zone
humidity ratio at the current time step. When the air system is operating, the mass flow for the
system outlet includes the infiltration mass flow rate, therefore the infiltration mass flow rate is
4/1/13 18
Integrated Solution Manager Moisture Predictor-Corrector
not included as a separate term in the air mass balance equation. But when the air system is
off, the infiltration mass flow in is then exhausted out of the zone directly.
In the same manner as described above for predicting the moisture load to be met by the air
system, the zone air moisture correction calculation will be described individually for the three
solution algorithms.
EulerMethod
Wzt t
N sl N surfaces N zones
airVz CW
N surfaces N zones
t
i 1
Ai hmi airz m m
i 1
i inf m sys
ThirdOrderBackwardDifference
airVz CW
N sl N surfaces N zones
3 1
kg masssched load Ai hmi airz Wsurfsi m W i zi m inf W m sysWsup
t
(3Wzt t Wzt 2 t Wzt 3 t )
2 3
Wzt i 1 i 1 i 1
airVz CW 11
N surfaces N zones
t
6
i 1
Ai hmi airz m m
i 1
i inf m sys
Using the same A, B, and C parameters from the prediction step modified with actual zone
mass flows with the air system ON and OFF result in:
N surfaces N zones
A i 1
Ai hmi airz m m
i 1
i inf m sys
N sl N surfaces N zones
B kg masssched load Ai hmi airz Wsurfsi mW
i zi m inf W m sysWsup
i 1 i 1 i 1
airVz CW
C
t
Else If (ZoneSupplyAirMassFlowRate <= 0.0) Then
N surfaces N zones
A i 1
Ai hmi airz m m
i 1
i inf m Exhaust
N sl N surfaces N zones
B kg masssched load Ai hmi airz Wsurfsi mW
i zi m inf W m ExhaustW
i 1 i 1 i 1
airVz CW
C
t
End If
4/1/13 19
Integrated Solution Manager Carbon Dioxide Predictor-Corrector
Inserting in the parameters A, B and C above in the air mass balance equation, it simplifies
to:
AnalyticalSolution
N sl N surfaces N zones
kg masssched load Ai hmi airzWsurfsi mW i zi m infW m sysW sup
Wzt t i 1 i 1 i 1
Wz
t N surfaces N zones
*
Ai hmi airz m i m inf m sys
i 1 i 1
N surfaces N zones
i mi airz m i m inf m sys
A h
exp i 1 i 1
t
airVz CW
N sl N surfaces N zones
kg masssched load
i 1
i 1
Ai hmi airzWsurfsi m W
i 1
i zi m infW m sysW sup
N surfaces N zones
i 1
Ai hmi airz m
i 1
i m inf m sys
The above solutions are implemented in the Correct Zone Air Humidity Ratio step in
EnergyPlus. This moisture update equation is used for the Conduction Transfer Function
(CTF) heat balance algorithm, in addition to the effective moisture penetration depth (EMPD)
with conduction transfer function heat balance algorithm. The equations are identical except
that the convection to the zone surfaces is non-zero for the moisture penetration depth case.
This moisture update allows both methods to be updated in the same way, with the only
difference being the additional moisture capacitance of the zone surfaces for the Effective
Moisture Penetration Depth (EMPD) solution approach.
When the HAMT (Combined Heat And Moisture Finite Element) defined in the
HeatBalanceAlgorithm object is applied, the moisture update equations are also the same as
the equations used in the effective moisture penetration depth (EMPD) with conduction
transfer function solution algorithm.
The transient air mass balance equation for the change in zone air carbon dioxide
concentration may be expressed as follows:
dCzt N sl N zones
airVz CCO 2 kgmass sched load *1.0 mi Czi Czt minf C Czt msys Csup Czt
6
dt i 1 i 1
where:
4/1/13 20
Integrated Solution Manager Carbon Dioxide Predictor-Corrector
N sl
kg
i 1
masssched load = sum of scheduled internal carbon dioxide loads. The zone air density is
used to convert the volumetric rate of carbon dioxide generation from user input into mass
6
generation rate [kg/s].The coefficient of 10 is used to make the units of carbon dioxide as
ppm.
N zones
m C
i 1
i zi Czt = carbon dioxide transfer due to interzone air mixing [ppm-kg/s]
C zi = carbon dioxide concentration in the zone air being transferred into this zone [ppm]
m inf C C zt = carbon dioxide transfer due to infiltration and ventilation of outdoor air
[ppm-kg/s]
C = carbon dioxide concentration in outdoor air [ppm]
m sys Csup Czt = carbon dioxide transfer due to system supply [ppm-kg/s]
Csup = carbon dioxide concentration in the system supply airstream [ppm]
m sys = air system supply mass flow rate [kg/s]
dC zt
airVz = carbon dioxide storage term in zone air [kg/s]
dt
C zt = zone air carbon dioxide concentration at the current time step [ppm]
air = zone air density [kg/m3]
Vz = zone volume [m3]
CCO2 = carbon dioxide capacity multiplier [dimensionless] (See the InputOutput Reference for
additional information on the object ZoneCapacitanceMultiplier:ResearchSpecial)
In the same manner as described above for zone air temperature (ref. Basis for the Zone and
Air System Integration), the solution algorithms provided in the ZoneAirHeatBalanceAlgorithm
object are also applied to the zone air carbon dioxide calculations.
In order to calculate the derivative term with respect to time, the first order backward finite
difference method, defined as the EulerMethod in the ZoneAirHeatBalanceAlgorithm object,
may be used:
dCzt
t (Czt C zt t ) O( t )
1
dt
The zone air carbon dioxide concentration update at the current time step using the
EulerMethod may be expressed as follows:
N sl N zones
airVZ CCO 2 t
1
C t
z C t t
z kg mass sched load *10
6
m C i zi Czt m inf C Czt m sys Csup Czt
i 1 i 1
To preserve the stability of the calculation of the zone carbon dioxide concentration, the third
order differential approximation, derived by a Taylor Series and used in the calculation of the
next time step’s zone air temperature, is also applied to the zone air carbon dioxide
4/1/13 21
Integrated Solution Manager Carbon Dioxide Predictor-Corrector
dCzt
11 t
Cz 3Czt t 32 C zt 2 t 13 Czt 3 t
6
O( t 3 )
dt t
The coefficients of the approximated derivative are very close to the coefficients of the
analogous Adams-Bashforth algorithm. Then the approximated derivative is substituted into
the mass balance and the terms with the carbon dioxide concentration at past time steps are
all put on the right-hand side of the equation. This third order derivative zone carbon dioxide
update increases the number of previous time steps that are used in calculating the new zone
carbon dioxide concentration, and decreases the dependence on the most recent. The higher
order derivative approximations have the potential to allow the use of larger time steps by
smoothing transitions through sudden changes in zone operating conditions.
This gives us the basic air mass balance equation that will be solved two different ways, one
way for the predict step and one way for the correct step.
Since the third choice of solution algorithms uses an integration approach, defined as
AnalyticalSolution in the ZoneAirHeatBalanceAlgorithm object, it does not require any
approximations and has no truncation errors. The solutions in both prediction and correction
are provided below in detail.
Carbon Dioxide Prediction
For the carbon dioxide concentration prediction case, the equation is solved for the
anticipated system response as shown below.
Since the program provides three solution algorithms, the carbon dioxide prediction from
each solution algorithm is given below.
EulerMethod
For this solution algorithm, the air mass balance for the predicted air system load or response
is:
4/1/13 22
Integrated Solution Manager Carbon Dioxide Predictor-Corrector
ThirdOrderBackwardDifference
For this solution algorithm, the air mass balance for the predicted system load or response is
given below:
V C N zones
t
PredictedSystemLoad [kg / sec] air z CO 2
t
11
6
m m i inf * Csetpoint
i 1
N sl N zones
airVz CCO 2
kg masssched load *106 m C i zi m inf C
t
3Czt t 32 Czt 2 t 13 Czt 3 t
i 1 i 1
AnalyticalSolution
For this solution algorithm, the air mass balance for the predicted air system load or response
is given below:
N zones
N zones
m i m inf
PredictedSystemLoad [kg / sec] m i m inf * Csetpoint
t
C zt t *exp i 1 t *
i 1 airVZ CCO 2
1
N zones
m i m inf N sl N zones
1 exp i 1 t kg masssched load *106 m i Czi m inf C
airVZ CCO 2 i 1 i 1
At the prediction point in the simulation, the system air mass flows are not known; therefore,
the system response is approximated. The predicted air system carbon dioxide load is then
used in the system simulation to achieve the best results possible. If a central HVAC system
provides the outdoor flow rate from a Controller:MechanicalVentilation object, the outdoor
airflow rate may be approximated as:
PredictedSystemLoad m sys Csup C zt m OA, z C Csetpoint
t
where:
m OA, z = supply outdoor airflow rate into the controlled zone [kg/s]
The above approximation is based on the assumption that the carbon dioxide concentration
at the outdoor air (OA) mixer inlet is equal to the zone air outlet concentration level, and the
carbon dioxide level at the zone supply air inlet is equal to the level at the outlet node of the
OA mixer.
After the system simulation is completed the actual response from the air system is used in
the carbon dioxide correction step, which is shown next.
Carbon Dioxide Correction
For the correct step the expanded air mass balance equation is solved for the final zone
carbon dioxide concentration at the current time step. In the same manner as described
above for predicting the carbon dioxide load to be met by the air system, the zone air carbon
dioxide correction calculation will be described individually for the three solution algorithms.
4/1/13 23
Integrated Solution Manager Generic Contaminant Predictor-Corrector
EulerMethod
N sl N zones
C zt t
kg masssched load *106 m i Czi m inf C m sys Csup airVZ CCO 2
t
C zt i 1 i 1
ThirdOrderBackwardDifference
N sl N zones
airVZ CCO 2 3 1
kgmasssched load *106 m C i zi m inf C m sys Csup
t
(3Czt t Czt 2 t Czt 3 t )
2 3
C zt i 1 i 1
N sl N zones
N zones
kg masssched load *10 mi Czi minf C msys Csys
6
m i m inf m sys
C zt C zt t i 1 i 1 *exp i 1 t
N zones
airVZ CCO 2
i 1
m i m inf m sys
N sl N zones
m
i 1
i m inf m sys
The above solutions are implemented in the Correct Zone Air Carbon Dioxide step in the
Zone Contaminant Predictor Corrector module of EnergyPlus.
Generic Contaminant Predictor-Corrector
The transient air mass balance equation for the change in zone air generic contaminant
concentration may be expressed as follows:
dC tf , z N source N zones
m C C tf , z m inf C f , C tf , z
Nsink
airVz M for
dt
i 1
air G f ,i *1.0 air 6
R i
f ,i C f ,z
i 1
i f , z ,i
m sys C f ,sup C tf , z h j Aj (
Cs , j
C f , z ) S f (C tf , z t )
j kj
where:
N source
i 1
air G f ,i = Sum of internal generic contaminant loads from sources in a zone or interior
surfaces.
The zone air density is used to convert the volumetric rate of generic contaminant generation
6
from user input into mass generation rate [kg/s].The coefficient of 10 is used to make the
units of generic contaminant as ppm.
Nsink
air R i
f ,i C f , z = Sum of removal rate from sinks in a zone or interior surfaces [ppm-kg/s]
4/1/13 24
Integrated Solution Manager Generic Contaminant Predictor-Corrector
N zones
m C
i 1
i f , z ,i C tf , z = Generic contaminant transfer due to interzone air mixing [ppm-kg/s]
C f , z ,i = Generic contaminant concentration in the zone air being transferred into this zone
[ppm]
m inf C f , C tf , z = Generic contaminant transfer due to infiltration and ventilation of
outdoor air [ppm-kg/s]
C f , = Generic contaminant concentration in outdoor air [ppm]
m sys C f ,sup C tf , z = Generic contaminant transfer due to system supply [ppm-kg/s]
C f ,sup = Generic contaminant concentration in the system supply airstream [ppm]
m sys = Air system supply mass flow rate [kg/s]
dC tf , z
airVz = Generic contaminant storage term in zone air [ppm-kg/s]
dt
C tf , z = Zone air generic contaminant concentration at the current time step [ppm]
air = Zone air density [kg/m3]
Vz = Zone volume [m3]
Cs , j
h A ( k
j
j j C f ,z ) = Generic contaminant transport through diffusion between
j
interior surfaces and zone air
S f (C tf,z t ) = Generic contaminant generation or removal rate as a function of zone air generic
contaminant level at the previous time step
Mfor = Generic contaminant capacity multiplier [dimensionless] (See the InputOutput
Reference for additional information on the object
ZoneCapacitanceMultiplier:ResearchSpecial)
In the same manner as described above for zone air temperature (ref. Basis for the Zone and
Air System Integration), the solution algorithms provided in the ZoneAirHeatBalanceAlgorithm
object are also applied to the zone air carbon dioxide calculations.
In order to calculate the derivative term with respect to time, the first order backward finite
difference method, defined as the EulerMethod in the ZoneAirHeatBalanceAlgorithm object,
may be used:
dC tf , z
t (C tf , z C tf , z t ) O ( t )
1
dt
The zone air generic contaminant concentration update at the current time step using the
EulerMethod may be expressed as follows:
4/1/13 25
Integrated Solution Manager Generic Contaminant Predictor-Corrector
N source N zones
C C tf , z t m C C tf , z
Nsink
airVz M for t
1 t
f ,z air G f ,i *1.06 air R f ,i C f , z i f , z ,i
i 1 i i 1
dC tf , z 11 t
C f ,z 3C tf, z t 32 C tf , z2 t 13 C tf , z3 t
6
O( t 3 )
dt t
The coefficients of the approximated derivative are very close to the coefficients of the
analogous Adams-Bashforth algorithm. Then the approximated derivative is substituted into
the mass balance, and the terms with the carbon dioxide concentration at past time steps are
all put on the right-hand side of the equation. This third order derivative zone carbon dioxide
update increases the number of previous time steps that are used in calculating the new zone
generic contaminant concentration and decreases the dependence on the most recent. The
higher order derivative approximations have the potential to allow the use of larger time steps
by smoothing transitions through sudden changes in zone operating conditions.
airVz M for 11 t Nsink N zones
air R f ,i C f , z m C m inf C tf , z m sys C tf , z h j Aj C tf , z
t t
f ,z
C
t
i f ,z
6 i i 1 j
N source N zones
Cs , j
i 1
air G f ,i *1.06 m C
i 1
i f , z ,i m inf C f , m sys C f ,sup h j Aj
j kj
Sf
airVz M for
t
3C tf,z t 23 C tf , z2 t 13 C tf, z3 t
This gives us the basic air mass balance equation that will be solved in two different ways,
one way for the predict step and one way for the correct step.
Since the third choice of solution algorithms uses an integration approach, defined as
AnalyticalSolution in the ZoneAirHeatBalanceAlgorithm object, it does not require any
approximations and has no truncation errors. The solutions in both prediction and correction
are provided below in detail.
Generic Contaminant Prediction
For the generic contaminant concentration prediction case, the equation is solved for the
anticipated system response as shown below.
Since the program provides three solution algorithms, the generic contaminant prediction
from each solution algorithm is given below.
EulerMethod
For this solution algorithm, the air mass balance for the predicted air system load or response
is:
4/1/13 26
Integrated Solution Manager Generic Contaminant Predictor-Corrector
i 1 i i 1
m inf C f , Csetpoint h j Aj ( s , j Csetpoint )
C
j kj
ThirdOrderBackwardDifference
For this solution algorithm, the air mass balance for the predicted system load or response is
given below:
N zones N sink
PredictedSystemLoad [kg / sec] m i m inf air R f ,i h j Aj *
i 1 i j
N zones N sink
i m
m inf air R f ,i h j Aj
C t C t t
*exp i 1 i j
t *
setpoint z
airVZ M FOR
1
N zones Nsink
i
m
m inf air R f ,i h j Aj
1 exp i 1 i j
t
airVZ M FOR
N source N zones
C
air f ,i
m i C f , z ,i m inf C f , h j Aj s , j S f
6
G *1.0
i 1 i 1 j kj
At the prediction point in the simulation, the system air mass flows are not known; therefore,
the system response is approximated. The predicted air system generic contaminant load is
then used in the system simulation to achieve the best results possible. If a central HVAC
system provides the outdoor flow rate from a Controller:MechanicalVentilation object, the
outdoor airflow rate may be approximated as:
4/1/13 27
Integrated Solution Manager Generic Contaminant Predictor-Corrector
N source N zones
Cs , j C tf , z t
i 1
air G f ,i *1.0 6
m C
i 1
i f , z ,i m inf C f , m sys C f ,sup h j Aj
kj
airVZ M FOR
t
Sf
t j
C f ,z Nsink N zones
airVz M for t air R m m m sys h j Aj
1
f ,i i inf
i i 1 j
ThirdOrderBackwardDifference
N source N zones
i 1
air G f ,i *1.06 m C
i 1
i f , z ,i m inf C f , m sys C f ,sup
Cs , j airVZ M FOR 3 1
h j Aj (3C tf, z t C tf , z2 t C tf , z3 t ) S f
kj t 2 3
C tf , z
j
1 11
Nsink N zones
airVz M for t air R f ,i m m i inf m sys h j Aj
6 i i 1 j
4/1/13 28
Integrated Solution Manager Summary of Time Marching Solution
AnalyticalSolution
N source N zones
C
air G f ,i *1.0 6
m i C f , z ,i m inf C f , m sys C f , sys h j Aj s , j
kj
C tf, z t *
i 1 i 1 j
C tf , z
N zones N sink
m i m inf m sys air R f ,i h j Aj
i 1 i j
N zones Nsink
mi minf msys air R f ,i h j A j
exp t
i 1 i j
airVZ M FOR
N source N zones
Cs , j
i 1
air G f ,i *1.0
6
m C
i 1
i f , z ,i m inf C f , m sys C f , sys h j Aj
j kj
Sf
N zones N sink
i 1
m i m inf m sys air R f ,i h j Aj
i j
The above solutions are implemented in the Correct Zone Air Generic Contaminant step in
the Zone Contaminant Predictor Corrector module of EnergyPlus.
EnergyPlus models building performance over time spans of days to years using a time
marching method based on timesteps. This section provides more information on issues
related to timestep formulation.
Summary of Timestep Model Formulation
An EnergyPlus simulation covers a certain period of time, such as a day or a year, that is
broken down into a series of discrete bins of time that are referred to as timesteps. The
program marches through time by recalculating model equations at each timestep. The figure
below diagrams some of these basic concepts.
4/1/13 29
Integrated Solution Manager Summary of Time Marching Solution
Most models in EnergyPlus are quasi-steady energy balance equations used to predict the
conditions present during each timestep. Various input data and boundary conditions for the
models are time-varying and a “staircase” approach is used where such values are
determined for a particular timestep and then essentially held constant over the entire
timestep. Predictions for state variables, such as temperature, are averages over the
timestep. Predictions for summed variables, such as energy use, are simple totals over the
timestep.
EnergyPlus produces time-series results for selected output variables at selected
frequencies. The time values associated with the time-series data, or timestamps, are output
for the end of the timestep, but the values of the variables are for the entire bin of time prior to
the timestamp. When data are reported at coarser frequencies such as hourly, then the
results are averages or simple totals for all the timesteps that are contained within the larger
bin of time.
To simplify solutions that would otherwise need to be simultaneous, models sometimes use
data that are “lagged” which means that the values used in the calculations for the current
timestep are actually results from the previous timestep. Many models just use the most
current results available and so may use lagged data for a first iteration, but then use current
data that are not lagged for subsequent iterations.
Zone Update Method
A zone is not necessarily a single room but is usually defined as a region of the building or a
collection of rooms subject to the same type of thermal control and having similar internal
load profiles that, subsequently, can be grouped together. Zones can interact with each other
thermally through adjacent surfaces and by intermixing of zone air. In EnergyPlus, the
conditions in each zone are updated by Equation (11), which uses previously calculated
values of the zone conditions. This means that EnergyPlus does not have to iterate to find a
self-consistent solution of the updated zone conditions. However, because heat transfer
through each zone's surfaces and interzone mixing of air still occur, the new space
temperatures must be computed at the same simulation time and on the same time step in all
zones, even though conditions in one zone may be changing much more rapidly than
conditions in the other zones. We have previously documented the method used to update
4/1/13 30
Integrated Solution Manager Summary of Time Marching Solution
the zone temperature at each time step. This method, called the predictor corrector method,
has proved to be satisfactory over a wide range of conditions.
Variable Timestep
The need for a variable timestep in EnergyPlus was identified during development of its
predecessor IBLAST. Prior to the integration of the central plant simulation in IBLAST, a time
step t for the zone temperature update of 0.25 hours (15 minutes) was found to give stable
results without a large increase in computation time. The first step in integrating plant was to
implement the detailed coil models and coil control strategies without actually adding the
plant models themselves. This meant that the user had to specify the coil water inlet
temperature and the maximum coil inlet water flow rate to run the simulation. The real life
analogy would be a chiller of very large, though not infinite, capacity. The coil capacity was
controlled by adjusting the water flow rate, but the effect of the plant on the chilled water
temperature was eliminated. After implementation of this step, experience with the program
showed that updating the zone temperatures on a fixed time step frequently resulted in
instabilities unless a very short time step was used. However, as the time step got shorter the
time required to execute the program got prohibitively high.
Clearly, an adaptive time step was required. This would shorten the time step to maintain
stability of the zone air energy balance calculation when zone conditions were changing
rapidly and expand it to speed computation when zone conditions were relatively unchanging.
But, the adaptive time step could not be applied easily to the surface heat transfer
calculations, even using interpolation methods to determine new temperature and flux
histories. The problem of updating the zone temperature was resolved by using a two-time-
step approach in which the zone air temperature is updated using an adaptive time step that
ensures stability. In this two time level scheme, the contributions to the zone loads from the
surfaces, and user specified internal loads are updated at the default or user specified time
step that is constant. This is often referred to as the “zone” time step. The contributions to the
zone loads from HVAC system response, infiltration, and air mixing are updated at a second
variable time step, referred to as the “system” time step. The system time step is limited to
between one minute and the zone time step. The lower limit of the system time step can be
increased to larger than one minute by the user with a System Convergence Limits object
(which may be useful to decrease simulation run times at the expense of some accuracy).
The program’s decision to adapt the time step is made by first using the usual zone time step
and executing the full predictor-corrector calculations to find resulting zone temperatures. The
maximum temperature change experienced by any zone in the model is determined. If this
maximum zone temperature change is more than a preset limit of 0.3°C, then the simulation
switches to using the shorter system time step. The number of system time steps (within a
particular zone time step) is modeled from the results for the maximum zone temperature
change (just obtained from the corrector) by assuming that temperature change is linear.
The number of system time steps indicated by the temperatures is:
EnergyPlus takes the smallest of these two, truncates them to a whole number and
calculates a system time step as:
4/1/13 31
Integrated Solution Manager Summary of Time Marching Solution
4/1/13 32
Integrated Solution Manager Summary of Time Marching Solution
Taylor, R.D., C.O. Pedersen, D.E. Fisher, R. J. Liesen, L.K. Lawrie. 1991. Impact of
Simultaneous Simulation of Buildings and Mechanical Systems in Heat Balance Based
Energy Analysis Programs on System Response and Control, Conference Proceedings
IBPSA Building Simulation '91, Nice, France, August 20-22, 1991.
4/1/13 33
Surface Heat Balance Manager / Processes Conduction Through The Walls
(t ) X jTo ,t j Y jTi ,t j
qko (33)
j 0 j 0
where q” is heat flux, T is temperature, i signifies the inside of the building element, o signifies
the outside of the building element, t represents the current time step, and X and Y are the
response factors.
While in most cases the terms in the series decay fairly rapidly, the infinite number of terms
needed for an exact response factor solution makes it less than desirable. Fortunately, the
similarity of higher order terms can be used to replace them with flux history terms. The new
solution contains elements that are called conduction transfer functions (CTFs). The basic
form of a conduction transfer function solution is shown by the following equation:
nz nz nq
qki (t ) Z oTi ,t Z jTi ,t j YoTo ,t Y jTo ,t j j qki ,t j (34)
j 1 j 1 j 1
nz nz nq
(t ) YoTi ,t Y jTi ,t j X oTo ,t X jTo ,t j j qko ,t j
qko (35)
j 1 j 1 j 1
4/1/13 34
Surface Heat Balance Manager / Processes Conduction Through The Walls
both the interior and exterior surface as well as some of the previous flux values at the
interior surface.
The final CTF solution form reveals why it is so elegant and powerful. With a single, relatively
simple, linear equation with constant coefficients, the conduction heat transfer through an
element can be calculated. The coefficients (CTFs) in the equation are constants that only
need to be determined once for each construction type. The only storage of data required
are the CTFs themselves and a limited number of temperature and flux terms. The
formulation is valid for any surface type and does not require the calculation or storage of
element interior temperatures.
Calculation of Conduction Transfer Functions
The basic method used in EnergyPlus for CTF calculations is known as the state space
method (Ceylan and Myers 1980; Seem 1987; Ouyang and Haghighat 1991). Another
common, older method used Laplace transformations to reach the solution; the Laplace
method was used in BLAST (Hittle, 1979; Hittle & Bishop, 1983). The basic state space
system is defined by the following linear matrix equations:
d x
A x B u
dt
y C x Du
where x is a vector of state variables, u is a vector of inputs, y is the output vector, t is time,
and A, B, C, and D are coefficient matrices. Through the use of matrix algebra, the vector of
state variables (x) can be eliminated from the system of equations, and the output vector (y)
can be related directly to the input vector (u) and time histories of the input and output
vectors.
This formulation can be used to solve the transient heat conduction equation by enforcing a
finite difference grid over the various layers in the building element being analyzed. In this
case, the state variables are the nodal temperatures, the environmental temperatures
(interior and exterior) are the inputs, and the resulting heat fluxes at both surfaces are the
outputs. Thus, the state space representation with finite difference variables would take the
following form:
T1
d
T1
Tn T
A B i
dt
Tn To
T1
q "i Ti
q " C D T
o Tn o
where T1, T2, ..., Tn-1, Tn are the finite difference nodal temperatures, n is the number of
nodes, Ti and To are the interior and exterior environmental temperatures, and qi and qo
are the heat fluxes (desired output).
4/1/13 35
Surface Heat Balance Manager / Processes Conduction Through The Walls
Seem (1987) shows that for a simple one layer slab with two interior nodes as in Figure 7 and
convection at both sides the resulting finite difference equations are given by:
dT1 T T
C hA To T1 2 1
dt R
dT2 T T
C hA Ti T2 1 2
dt R
q "i h Ti T2
q "o h T1 To
where:
R ,
kA
c p A
C , and
2
A is the area of the surface exposed to the environmental temperatures.
In matrix format:
dT1 1 hA 1 hA
dt RC C T1 C 0
RC To
dT2 1 1 hA T2 hA Ti
0
dt RC RC C C
q "o 0 h T1 0 h To
q " h 0 T h 0 T
i 2 i
4/1/13 36
Surface Heat Balance Manager / Processes Conduction Through The Walls
To T1 T2 Ti
1 R 1
hA hA
C C
The important aspect of the state space technique is that through the use of matrix algebra
the state space variables (nodal temperatures) can be eliminated to arrive at a matrix
equation that gives the outputs (heat fluxes) as a function of the inputs (environmental
temperatures) only. This eliminates the need to solve for roots in the Laplace domain. In
addition, the resulting matrix form has more physical meaning than complex functions
required by the Laplace transform method.
The accuracy of the state space method of calculating CTFs has been addressed in the
literature. Ceylan and Myers (1980) compared the response predicted by the state space
method to various other solution techniques including an analytical solution. Their results
showed that for an adequate number of nodes the state space method computed a heat flux
at the surface of a simple one layer slab within 1% of the analytical solution. Ouyang and
Haghighat (1991) made a direct comparison between the Laplace and state space methods.
For a wall composed of insulation between two layers of concrete, they found almost no
difference in the response factors calculated by each method.
Seem (1987) summarizes the steps required to obtain the CTF coefficients from the A, B, C,
and D matrices. While more time consuming than calculating CTFs using the Laplace
Transform method, the matrix algebra (including the calculation of an inverse and exponential
matrix for A) is easier to follow than root find algorithms. Another difference between the
Laplace and State Space methods is the number of coefficients required for a solution. In
general, the State Space method requires more coefficients. In addition, the number of
temperature and flux history terms is identical (nz=nq). Note that as with the Laplace method
that the actual number of terms will vary from construction to construction.
Two distinct advantages of the State Space method over the Laplace method that are of
interest when applying a CTF solution for conduction through a building element are the
ability to obtain CTFs for much shorter time steps and the ability to obtain 2- and 3-D
conduction transfer functions. While not implemented in the Toolkit, both Seem (1987) and
Strand (1995) have demonstrated the effectiveness of the State Space method in handling
these situations that can have important applications in buildings.
Conduction Transfer Function (CTF) Calculations in EnergyPlus
Conduction transfer functions are an efficient method to compute surface heat fluxes
because they eliminate the need to know temperatures and fluxes within the surface.
However, conduction transfer function series become progressively more unstable as the
4/1/13 37
Surface Heat Balance Manager / Processes Conduction Through The Walls
time step decreases. This became a problem as investigations into short time step
computational methods for the zone/system interactions progressed because, eventually, this
instability caused the entire simulation to diverge. This phenomenon was most apparent for
thermally massive constructions with long characteristic times and, correspondingly, requiring
a large number of terms in the CTF series. This indicates that the problem is related to round-
off and truncation error and is in no way an indictment of the CTF method itself. Methods that
develop CTF series from finite difference approximations to the heat conduction equation
(Meyers, 1980; Seem, 1987) were considered to address this problem. Seem's method did
give better accuracy and stability at short time steps than the current BLAST technique but,
the method still had difficulty computing stable CTF series for time steps of less than 1/4 hour
for the heaviest constructions in the BLAST library.
The zone heat gains consist of specified internal heat gains, air exchange between zones, air
exchange with the outside environment, and convective heat transfer from the zone surfaces.
Of these, the surface convection load requires the most complicated calculations because a
detailed energy balance is required at the inside and outside surface of each wall, floor, and
roof. In addition, the transient heat conduction in the material between the surfaces must be
solved. This solution gives the inside and outside temperatures and heat fluxes that must be
known in order to calculate the convection component to the zone load for each zone
surface. BLAST uses a conduction transfer function CTF method attributed to Hittle (1980) to
solve the transient conduction problem for each surface. The method results in a time series
of weighting factors that, when multiplied by previous values of the surface temperatures and
fluxes and the current inside and outside surface temperatures, gives the current inside and
outside heat flux. The method is easily applied to multilayered constructions for which
analytical solutions are unavailable. In addition, determining the series of CTF coefficients is
a one-time calculation, making the method much faster than finite difference calculations.
A problem with CTF methods is that the series time step is fixed; that is, a CTF series
computed for a one hour time step takes information at t-1 hours, t-2 hours, etc. and
computes conditions at the current time t. As time advances the oldest term in the input
series is dropped and the data moved back one time step to allow the newest value to be
added to the series. For convenience, the time step used to determine the CTF series should
be the same as the time step used to update the zone mean air temperature in the zone
energy balance. But, as the time step used to calculate the CTF series gets shorter, the
number of terms in the series grows. Eventually, with enough terms, the series becomes
unstable due to truncation and round-off error. Heavy constructions, such as slab-on-grade
floors (12" heavyweight concrete over 18" dirt), have accuracy and stability problems at time
steps as large as 0.5 hours when modeled by Hittle's CTF method. In an attempt to
overcome this problem, Hittle's method was replaced by Seem's method (1987) in IBLAST.
This resulted in some improvement in stability at shorter time steps, but not enough to allow
IBLAST to run at a 0.1 hour time step without restricting the types of surfaces that could be
used.
Even though CTF methods require that values of the surface temperatures and fluxes be
stored for only a few specific times before the current time, the temperature and flux histories
are, actually, continuous functions between those discrete points. However, there is no way
to calculate information at these intermediate times once a series has been initialized. The
terms in the temperature and flux histories are out of phase with these points. However, they
can be calculated by shifting the phase of the temperature and flux histories by only a fraction
of a time step. This procedure would allow a CTF series computed for a time step t, to be
used to compute information at times t+t/2, t+t/3, t+t/4, or any other arbitrary fraction of
the time step, so long as the surface temperatures and flux values were still t apart. Several
ways of doing this are described below.
The method shown in the Figure 10 maintains two sets of histories out of phase with each
other. The figure shows how this would work for two sets of histories out of phase by one
half of a time step. More sets of temperature and flux histories could be used, allowing the
simulation time step to take on values: 1/3, 1/4, 1/5, etc., of the minimum time step allowed
for the CTF calculations. The time step between inputs to the CTF series would be the
4/1/13 38
Surface Heat Balance Manager / Processes Conduction Through The Walls
smallest convenient interval at which the CTF series is stable. This scenario is illustrated in
this figure for two separate sets of temperature and flux histories. Cycling through each
history, in order, allowed calculations of the zone energy balance to be performed with
updated surface information at a shorter time step than one CTF history series would
otherwise allow. This method required no interpolation between the series once each set of
histories was initialized. However, if the smallest time step for a stable CTF series was large
compared to the zone temperature update time step, significant memory was required to
store all the sets of histories.
x x x x x x x history 1
o o o o o o o history 2
x x x x x
o o o o o
dt x x x x x
time
Figure 10. Multiple, staggered time history scheme
Another method is shown in Figure 11. Sequential interpolation of new histories that uses
successive interpolations to determine the next set of temperature and flux histories. The
current history is interpolated directly from the previous history set using the required time
phase shift between the two. This method required permanent storage for only one set of
temperature and flux histories at a time, but smoothed out temperature and flux data as more
interpolations were performed. As a result, at concurrent simulation times current values of
history terms were different form previous "in phase" history terms. This was unacceptable
from, a physical point of view, because it allowed current information to change data from a
previous time.
x x x x x x x history 1
o o o o o o o history 2
A final method, shown in Figure 12. Master history with interpolation, was something of a
hybrid of the previous two methods. One "master" history set was maintained and updated
for all time; this solved the problem of current events propagating information backwards in
time. When surface fluxes needed to be calculated at times out of phase with this master
history a new, temporary history was interpolated from the master values. This method
proved to be the best of the three options described because it eliminated propagation of
information backwards in time and only required concurrent storage of two sets of
temperature and flux histories. This method was subsequently incorporated into the IBLAST
program in conjunction with Seem's procedure for calculating the coefficients of the CTF
series.
4/1/13 39
Surface Heat Balance Manager / Processes Conduction Through The Walls
x x x x x
x x history 1
o o o o o o o history 2
x x x x x
o o o o o
dt x x x x x
time
Figure 12. Master history with interpolation
4/1/13 40
Surface Heat Balance Manager / Processes Conduction Through The Walls
The case where a resistance-only layer is defined anywhere except the inner or outer layer of
a construction is handled by treating the “no mass” layer as a single node layer. This will
result in a node at each interface as in the standard material layer cases. When a “no mass”
material is present, the R-Value only layer will not add any thermal capacitance to the nodes
at the interfaces at either side of the material. It will simply add resistance between the two
nodes.
From the EnergyPlus code, the A matrix (AMat) is assigned with values at the interface using
the following equations (taken from the actual code):
cap = ( rho(Layer)*cp(Layer)*dx(Layer) &
+rho(Layer+1)*cp(Layer+1)*dx(Layer+1) ) * 0.5D0
4/1/13 41
Surface Heat Balance Manager / Processes Conduction Finite Difference Solution Algorithm
References
Ceylan, H. T., and G. E. Myers. 1980. Long-time Solutions to Heat Conduction Transients
with Time-Dependent Inputs. ASME Journal of Heat Transfer, Volume 102, No. 1, pp. 115-
120.
Hittle, D. C. 1979. Calculating Building Heating and Cooling Loads Using the Frequency
Response of Multilayered Slabs, Ph.D. Thesis, University of Illinois, Urbana, IL.
Hittle, D. C., and R. Bishop. 1983. An Improved Root-Finding Procedure for Use in
Calculating Transient Heat Flow Through Multilayered Slabs. International Journal of Heat
and Mass Transfer, Vol. 26, No. 11, pp. 1685-1693.
Ouyang, K., and F. Haghighat. 1991. A Procedure for Calculating Thermal Response Factors
of Multi-layered Walls--State Space Method. Building and Environment, Vol. 26, No. 2, pp.
173-177.
Seem, J. E. 1987. Modeling of Heat Transfer in Buildings, Ph.D. Thesis, University of
Wisconsin, Madison, WI.
Strand, R. K. 1995. Heat Source Transfer Functions and Their Application to Low
Temperature Radiant Heating Systems, Ph.D. Thesis, University of Illinois, Urbana, IL.
Taylor, R. D., C.O. Pedersen, D.E. Fisher, R. J. Liesen, L.K. Lawrie. 1990. Simultaneous
Simulation of Buildings and Mechanical Systems in Heat Balance Based Energy Analysis
Programs, Proceedings of the 3rd International Conference on System Simulation in
Buildings, Liege, Belgium, December 3-5, 1990.
Taylor, R.D., C.O. Pedersen, D.E. Fisher, R. J. Liesen, L.K. Lawrie. 1991. Impact of
Simultaneous Simulation of Buildings and Mechanical Systems in Heat Balance Based
Energy Analysis Programs on System Response and Control, Conference Proceedings
IBPSA Building Simulation '91, Nice, France, August 20-22, 1991.
4/1/13 42
Surface Heat Balance Manager / Processes Conduction Finite Difference Solution Algorithm
change energy accurately. The implicit formulation for an internal node is shown in the
equation below.
kW
Ti j11 Ti j 1
kE
Ti j11 Ti j 1
j 1 Δx Δx
Ti Ti j
1 (36)
C p Δx
Δt 2
kW
Ti j1 Ti j kE
Ti j1 Ti j
Δx Δx
Where:
T = node temperature
Subscripts:
i = node being modeled
i+1 = adjacent node to interior of construction
i-1 = adjacent node to exterior of construction
j+1 = new time step
j = previous time step
t = calculation time step
x = finite difference layer thickness (always less than construction layer thickness)
Cp = specific heat of material
kw = thermal conductivity for interface between i node and i+1 node
kE = thermal conductivity for interface between i node and i-1 node
= density of material
Then, this equation is accompanied by a second equation that relates enthalpy and
temperature.
C p Δx
Ti j 1 Ti j
kW
Ti j11 Ti j 1
kE
Ti j11 Ti j 1
Δt Δx Δx
For both schemes, EnergyPlus uses the following four types of nodes, as shown in the figure
below (1) interior surface nodes, (2) interior nodes, (3) material interface nodes and (4)
external surface nodes. The grid for each material is established by specifying a half node for
each edge of a material and equal size nodes for the rest of the material. Equations such as
(36) are formed for all nodes in a construction. The formulation of all node types is basically
the same.
4/1/13 43
Surface Heat Balance Manager / Processes Conduction Finite Difference Solution Algorithm
In the CondFD model, surface discretization depends on the thermal diffusivity of the material
(α) and time step (Δt) selected, as shown in the equation below. The default value of 3 for the
space discretization constant, C, is basically the inverse of the Fourier Number
Fo t / x 2
and is based on the stability requirement for the explicit mode that
requires values higher than 2, or a Fourier number lower than 0.5. However, CondFD uses
implicit schemes that do not have the same stability requirements as the explicit mode. Thus,
the default 3 was originally set rather arbitrary. As of version 7, the value of this constant can
be controlled by the user with the input field called Space Discretization Constant in the
HeatBalanceSettings:ConductionFiniteDifference input object. The discretization method
allows CondFD to assign different node spacing or grid size to different material layers in a
wall or roof, as building walls and roofs typically consist of several layers of different materials
having different thermal properties.
Δx C α Δt
The actual integer number of nodes for each layer is then calculated by rounding off the
result from dividing the length of the material layer by the result of the equation above. After
this, Δx is recalculated by dividing the length of the material by the number of nodes. A full
node is equal to two half nodes. Lower values for the Space Discretization Constant yield
more nodes, with higher values yield fewer nodes.
Because the solution is implicit, a Gauss-Seidell iteration scheme is used to update to the
new node temperatures in the construction and under-relaxation is used for increased
stability. The Gauss-Seidell iteration loop is the inner-most solver and is called for each
surface. It is limited to 30 iterations but will exit early when the sum of all the node
temperatures changes between the last call and the current call, normalized by the sum of
the temperature values, is below 0.000001C. This convergence criteria is typically met after
3 iterations, except when PCMs are simulated as it takes an average of 2-3 more iterations
when PCM are changing phase. If the number if iterations needed to met convergence
criteria start to increase, an automatic internal relaxation factor stabilities the solution and in
most cases keep the number of iterations less than 10.
EnergyPlus also uses a separate, outer iteration loop across all the different inside surface
heat balances so that internal long-wave radiation exchange can be properly solved. For
CTF formulations, this iteration is controlled by a maximum allowable temperature difference
of 0.002C for inside face surface temperatures from one iteration to the next (or a limit of 100
iterations). CondFD uses the same default value for allowable temperature difference as
CTF. However, this parameter was found to often need to be smaller for stability and so the
inside surface heat balance manager uses a separate allowable maximum temperature
4/1/13 44
Surface Heat Balance Manager / Processes Conduction Finite Difference Solution Algorithm
difference when modeling CondFD. The user can control the value of the relaxation factor by
using the input field called Inside Face Surface Temperature Convergence Criteria in the
HeatBalanceSettings:ConductionFiniteDifference input object. In addition, if the program
detects that there is instability by watching for excessive numbers of iterations in this outer
loop and may decrease the relaxation factor. Users can also output the number of iterations
inside of CondFD loop for each surface and the outer internal heat balance loop for each
zone with “CondFD Inner Solver Loop Iterations” and “Heat Balance Inside Surfaces
Calculation Iterations” respectively.
Because of the iteration scheme used for CondFD, the node enthalpies get updated each
iteration, and then they are used to develop a variable Cp if a phase change material is being
simulated. This is done by including a third equation for Cp.
hi,new hi,old
Cp
Ti,new Ti,old (38)
The iteration scheme assures that the correct enthalpy, and therefore the correct Cp is used
in each time step, and the enthalpy of the material is accounted for accurately. Of course, if
the material is regular, the user input constant Cp is used.
The algorithm also has a provision for including a temperature coefficient to modify the
thermal conductivity. The thermal conductivity is obtained from:
where:
ko is the 20C value of thermal conductivity (normal idf input)
k1 is the change in conductivity per degree temperature difference from 20C
As of Version 7, the CondFD implementation was changed to evaluate the thermal
conductivity at the interface between nodes, as shown below. In this case, EnergyPlus uses a
linear interpolation between nodal points.
C p Δx
Ti j 1 Ti j 1
kW
Ti j11 Ti j 1
kE
Ti j11 Ti j 1 k T j
i 1 Ti j k T j
i 1
Ti j
2
W E
Δt Δx Δx Δx Δx
Where,
kW
k j 1
i 1 ki j 1
2
kE
k j1
i 1 k ij1
2
These additional property information values are put into the input file as explained in the
Input/Output Reference Document, but it consists simply of a value for k1 and set of enthalpy
temperature pairs that describe the enthalpy of the phase change material in straight line
segments with respect to temperature.
4/1/13 45
Surface Heat Balance Manager / Processes Conduction Finite Difference Solution Algorithm
A graph showing the effect of a large PCM on the outside surface of a zone is shown below.
The phase change temperature was 30C, and the flat temperature response during the
phase change is obvious. This example was run with a zone time step of one minute to show
that such small time steps can be done with the finite difference solution technique. It is more
efficient to set the zone time step shorter than those used for the CTF solution algorithm. It
should be set to 20 time steps per hour or greater, and can range up to 60. The finite
difference algorithm actually works better with shorter zone time steps. The computation time
has a minimum at a zone time step around two minutes (30 time steps/hr), and increases for
shorter or longer zone time steps.
4/1/13 46
Surface Heat Balance Manager / Processes Combined Heat and Moisture Transfer (HAMT) Model
References
Pedersen C.O., Enthalpy Formulation of conduction heat transfer problems involving latent
heat, Simulation, Vol 18, No. 2, February 1972
Versteeg, H. and Malalasekra, W. 1996. An introduction to computational fluid dynamic: the
finite volume method approach. Prentice Hall.
Tabares-Velasco, P.C. and Griffith, B. 2012. Diagnostic Test Cases for Verifying Surface
Heat Transfer Algorithms and Boundary Conditions in Building Energy Simulation Programs,
Journal of Building Performance Simulation, doi:10.1080/19401493.2011.595501
Tabares-Velasco, P.C., Christensen, C. and Bianchi, M. 2012. Verification and Validation of
EnergyPlus Phase Change Material Model for Opaque Wall Assemblies, Building and
Environment 54: 186-196.
Overview
The combined heat and moisture transfer finite (HAMT) solution algorithm is a completely
coupled, one-dimensional, finite element, heat and moisture transfer model simulating the
movement and storage of heat and moisture in surfaces simultaneously from and to both the
internal and external environments. As well as simulating the effects of moisture buffering,
HAMT is also be able to provide temperature and moisture profiles through composite
building walls, and help to identify surfaces with high surface humidity.
HAMT Nomenclature
w
Dependencies on moisture content are indicated by a superscript , on heat by a superscript
h v
and vapor pressure by a superscript .
Table 2. Combined Heat and Moisture Transfer Model Nomenclature
4/1/13 47
Surface Heat Balance Manager / Processes Combined Heat and Moisture Transfer (HAMT) Model
2
Dw m /s Liquid Transport Coefficient
2
A m Contact Surface area
3
Vi m Cell Volume
t s Time
s Time step between calculations
x m Distance between cell centres
C h J/C Heat Capacitance of cell i
i
H T w T T
k hv
x x
(40)
T x x
The three terms in equation (40) describe the storage, transport and generation of heat
respectively.
w w w T
x x x x
D (41)
The three terms in equation (41) describe the storage of moisture, the transport of liquid
moisture and the transport of vapor respectively. The equation to calculate the vapor diffusion
coefficient in air ( ) used in the third term of both equations, is also taken from Künzel,
2 10 7
T 273.15
0.81
(42)
Pambient
4/1/13 48
Surface Heat Balance Manager / Processes Combined Heat and Moisture Transfer (HAMT) Model
The heat storage capacity ( H ) depends on the moisture content w of the material by
T
the following equation.
H
c cw w (43)
T
The moisture content of the material w and the vapor diffusion resistance factor μ depend on
the relative humidity inside the material. The parameters w , k
w
and D w are also
moisture dependent parameters.
The following sections describe how the above equations are used within the HAMT model.
Surfaces, Material Layers and Cells
“Surfaces” are made of a number of layers of potentially any combination of materials. Each
surface is split into its constituent materials and is then split up further into cells through its
depth. HAMT will generate no more than 10 cells per material with widths that are thinner
near the boundaries of each material where most changes are expected and detail is needed.
Heat Transfer
Equation 1 can be re-written and used to describe the heat storage and transfer through the
th
i cell in a surface.
p 1 p 1
Ti p 1 ij p jp 1 pip 1
ci i c wi Vi Ti Ti kijw Aij j
p
T
w
hv Aij (44)
j xij j ij xij
In the one dimensional case there are only two adjacent cells each labelled j. The heat
v
generated due to vaporisation qi can be calculated separately.
ij p jp 1 pip 1
qiv hv Aij (45)
j ij xij
adds
Rearranging equation (40) and including other sources of heat ( qi ) such as radiation from
other surfaces in the calculation gives the temperature in a cell in the next time step as,
T jp 1 Ti p
j
Rijh
qiv qiadds Cih
Ti p 1 (46)
Cih 1
j h
Rij
xij
where Ci ci i c wi Vi is thermal heat capacitance of cell i and Rijh
h w
kij Aij is the
thermal resistance between cells i and j.
th
This equation can be solved using the Gauss-Seidel iteration technique. The i cell
th
temperature is calculated whilst the j cell temperatures are kept as up to date as possible.
4/1/13 49
Surface Heat Balance Manager / Processes Combined Heat and Moisture Transfer (HAMT) Model
The iteration is stopped when the maximum difference between two consecutive calculations
in all cells is less than a threshold of 0.002°C.
Moisture Content w
The moisture content (w) of a cell is needed for the calculation of the heat transfer through
the cell as it affects the thermal resistance and heat capacitance. The moisture content of
cells is calculated from the relative humidity (RH) of the material. The relationship between w
and the RH for each material is known as the sorption isotherm and measured data points
are entered into EnergyPlus as a series of coordinates. HAMT interpolates between the
measurements to obtain the moisture content of a material for any RH value. The sorption
isotherm input is via the MaterialProperty:HeatAndMoistureTransfer:SorptionIsotherm object
and is described in the Input Output Reference document.
Porosity P
The porosity of a material (P) is an input variable and defined as the maximum fraction, by
volume, of a material that can be taken up with moisture. It is used to calculate the maximum
point on the sorption isotherm curve. The porosity is entered for each material via the
MaterialProperty:HeatAndMoistureTransfer:Settings object, as described in the Input Output
Reference document.
w
Moisture Dependant Thermal Conductivity k
w
The thermal conductivity (k ) of the cell is determined by interpolating between data points of
thermal conductivity versus the moisture content of the material, entered into EnergyPlus via
the MaterialProperty:HeatAndMoistureTransfer:ThermalConductivity object. The moisture
content is determined via the sorption isotherm which gives the moisture content as a
function of Relative Humidity.
Moisture Dependant Moisture Diffusion Coefficient μ
This is used in the third term of equation (40) to describe the heat transfer due to vapor
movement. It is determined by interpolating between data points of moisture diffusion
coefficient versus the moisture content of the material, entered into EnergyPlus via the
MaterialProperty:HeatAndMoistureTransfer:Diffusion object. A simple linear interpolation is
used to obtain the conductivity between measured points.
Moisture Transfer
Moisture, as well as heat, is transported through materials as either liquid (w) or vapor (p).
There are two different potentials that control the movement though the material. Liquid
transfer is driven by differences in relative humidity whereas vapor transfer is driven by
differences in vapor pressure. Materials also have a capacity to store moisture. Equation (41)
can be re-written for a discrete cell in a continuous material.
dw p 1 i p jp 1 i p 1 ij p jp 1 pip 1
Vi i kij Aij Aij (47)
di j xij j ij xij
th
Equation (47) can be rearranged to provide the relative humidity of the i cell in the next time
step.
jp 1 pip 1 w i
p
j R w j R v Ci
i
p 1
w
ij ij
(48)
Ci 1 p sat
j w j iv
Rij Rij
4/1/13 50
Surface Heat Balance Manager / Processes Combined Heat and Moisture Transfer (HAMT) Model
xij
Rijw (49)
dw
Aij Dijw
d
ij xij
is the moisture resistance between cells i and j and Rijv
Aij ij is the vapor resistance
between cells i and j.
Equation (48) can be used together with the heat equation (46) in an alternate step by step
fashion to calculate the new temperature and relative humidity profiles for each cell for the
next time step.
w
Liquid Transport Coefficient D
The Moisture Dependant Liquid Transport Coefficient is entered as a series of moisture
density and liquid transport coefficient data points. There are two different coefficients, one
for suction, where the surface is wet due to rain, and one for redistribution where the surface
is no longer wet. If the weather file has a rain flag it is used to switch between these two
types of coefficient. HAMT-SUCTION and HAMT-REDISTRIBUTION.
4/1/13 51
Surface Heat Balance Manager / Processes Effective Moisture Penetration Depth (EMPD) Model
Reference DataSets (MoistureMaterials.idf). The properties were synthesised from the Annex
24 database [Kumar Kumaran, M. (1996)], supplemented, when required, by data from the
database of the WUFI model [WUFI (1999)] and are therefore not related to any unique,
measured material. Users should consult material property catalogues and other primary
sources when the properties of a specific material are required.
Moisture and heat from the surfaces are used by EnergyPlus to calculate the room air
temperature and moisture content. EnergyPlus with HAMT works best with as short a time
step as possible. However the optimum time step which gives a good prediction for a short
computing time will very much depend on the nature of the weather and type of building.
Buildings with frequent and large changes in internal and external temperature will need a
small time step, maybe even 60 steps per hour. Slowly evolving temperatures and relative
humidity’s will not require such a short time step and 20, or even 6, steps per hour may be
sufficient.
References
Künzel, H.M. (1995) Simultaneous Heat and Moisture Transport in Building Components.
One- and two-dimensional calculation using simple parameters. IRB Verlag 1995
Holman, J.P. (2002) Heat Transfer, Ninth Edition. McGraw-Hill
Winterton, R.H.S. (1997) Heat Transfer. (Oxford Chemistry Primers; 50) Oxford University
Press
Kumar Kumaran, M. (1996) IEA ANNEX 24, Final Report, Volume 3
WUFI (1999) version 2.2 Simultaneous Heat and Moisture Transport in Building components.
Fraunhofer IBP, Holzkirchen, Germany
Overview
Moisture has little effect on heating system performance, but a profound effect on the
performance of air conditioning systems. In order to accurately describe building
performance during periods when cooling is needed, it is very important to know the moisture
conditions of the building. If one assumes that all building moisture is contained in the room
air, then one ignores the fact that the materials that bound the room (e.g. wall surfaces,
furnishings, linens, etc.) store and release moisture. Thus, to assume that the only moisture
that effects cooling system performance is contained in the room air is a false, and it can lead
to significant error in the prediction of room moisture conditions and cooling system loads.
The EMPD (Effective Moisture Penetration Depth) model is a simplified, lumped approach to
simulate surface moisture adsorption and desorption.
EMPD Model Description
The EMPD concept assumes that a thin layer (δM) close to the wall surface behaves
dynamically and exchanges moisture with the air domain when exposed to cyclic air moisture
pulses. For short periods where the cyclic integral of the total moisture adsorption and
desorption is near zero (i.e. there is no net moisture storage), the EMPD concept has been
shown to be a reasonable approximation of reality (Kerestecioglu et al, 1989). In other
words, the following constraint must be met:
2 dU
1 d
d 0 (50)
where, τ2-τ1 denotes the finite time interval over which the equation holds. The EMPD model
assumes no spatial distribution of moisture content across the thickness (L) of the solid;
4/1/13 52
Surface Heat Balance Manager / Processes Effective Moisture Penetration Depth (EMPD) Model
rather, a thin layer (δM) of uniform moisture content (U) is assumed to represent the total
moisture content of the solid. This may be mathematically stated as:
L
0
U ( x)dx U M (51)
For most building materials, the equilibrium moisture sorption isotherm can be defined by the
following general equation (Kerestecioglu et al. 1988):
U a b c d (52)
where
W*
(53)
Wsat *
and
1 4111
Wsat * exp 23.7093 * (54)
Rv aT *
T 35.45
* *
Given that U=U(W ,T ), the moisture content may be differentiated with respect to time in the
following manner:
du U dW * U dT * dW * dT *
A B (55)
d W * d T * d d d
T
where AT and Bρ are the isothermal moisture capacity and thermo-gradient coefficient,
respectively. From Eqs. (53), (53) and (54), they can be expressed as:
ab b cd d
AT (56)
W*
and
1 4111
B * * 2
*(ab b cd d ) (57)
T (T 35.45)
The lumped mass transfer equation for the i-th solid domain may be written as
dU i
( A b M ) i hM ,i Ai (Wr Wi * ) (58)
d
Using Eqs. (55), (56), (57) and (58), one obtains the final equation needed for closure
moisture transfer at internal surface.
4/1/13 53
Surface Heat Balance Manager / Processes Effective Moisture Penetration Depth (EMPD) Model
dWi* dTi*
( Ai b M AT )i hM ,i Ai (Wr Wi ) ( Ab M B )i
*
(59)
d d
The energy equation for the envelope contains the surface temperature and is given by the
conduction equation
dT
C p ( k T ) (60)
d
k T qT " hT (T * Tr ) hM (W * Wr ) (61)
A more detailed account of the numerical solution procedure can be found in Kerestecioglu et
al. (1988).
EMPD Value Determination
An effective moisture penetration depth may be determined from either experimental or
detailed simulation data by using actual surface areas and moisture vapor diffusivity. An
empirical function derived from the detailed simulation may be used to determine the EMPD
value (Kerestecioglu et al, 1989):
where
(63)
Figure 16 gives the EMPD values to be used for various vapor diffusivities evaluated at
different ambient excitations.
4/1/13 54
Surface Heat Balance Manager / Processes Effective Moisture Penetration Depth (EMPD) Model
12
8
ξ=0.02
ξ=0.04
6
ξ=0.08
ξ=0.16
4
0
0 5 10 15 20
2 6
Vapor diffusivity, Dv (m /s x 10 )
Figure 16. Limit of Effective Penetration Depth Values for Various Vapor Diffusivities at Different Ambient
Excitations.
EMPD Nomenclature
2
A = Area [m ]
3
AT = Isothermal moisture capacity [m /kg]
Bρ = Thermo-gradient coefficient [kg/kg-K]
Cp = Specific heat [J/kg.K]
2
hM = Convective mass transfer coeff. [kg/m -s]
2
hT = Convective heat transfer coeff. [W/m -K]
k = Thermal conductivity [W/m-K]
L = Length [m]
2
q"T = Imposed heat flux [W/m ]
Rv = Ideal gas constant [461.52 J/kg-K]
T = Temperature [K]
U = Moisture content [kg/kg]
W = Humidity ratio [kg/kg]
Greek letters
δM = Effective penetration depth for moisture equation [m]
λ = Heat of vaporization [J/kg]
3
ρ = Density [kg/m ]
τ = Time [s]
φ = Relative humidity [0 to 1]
ξ = Ambient moisture excitation rate [1/h]
4/1/13 55
Surface Heat Balance Manager / Processes Outside Surface Heat Balance
Shortwave radiation,
including direct,
reflected, and diffuse Conduction
sunlight into wall,
qko
Longwave radiation
from the
environment
Convective exchange
with outside air Wall
Outside
Face
qsol qLWR
qconv
qko
0 (64)
where:
qsol = Absorbed direct and diffuse solar (short wavelength) radiation heat flux.
= Net long wavelength (thermal) radiation flux exchange with the air and surroundings.
qLWR
= Convective flux exchange with outside air.
qconv
= Conduction heat flux (q/A) into the wall.
qko
All terms are positive for net flux to the face except the conduction term, which is traditionally
taken to be positive in the direction from outside to inside of the wall. Simplified procedures
generally combine the first three terms by using the concept of a sol-air temperature. Each of
these heat balance components is introduced briefly below.
4/1/13 56
Surface Heat Balance Manager / Processes Outside Surface Heat Balance
qsol is calculated using procedures presented later in this manual and includes both direct
and diffuse incident solar radiation absorbed by the surface face. This is influenced by
location, surface facing angle and tilt, surface face material properties, weather conditions,
etc.
External Longwave Radiation
is a standard radiation exchange formulation between the surface, the sky, and the
qLWR
ground. The radiation heat flux is calculated from the surface absorptivity, surface
temperature, sky and ground temperatures, and sky and ground view factors.
The longwave radiation heat exchange between surfaces is dependent on surface
temperatures, spatial relationships between surfaces and surroundings, and material
properties of the surfaces. The relevant material properties of the surface, emissivity and
absorptivity , are complex functions of temperature, angle, and wavelength for each
participating surface. However, it is generally agreed that reasonable assumptions for
building loads calculations are (Chapman 1984; Lienhard 1981):
each surface emits or reflects diffusely and is gray and opaque ( = , = 0, = 1- )
each surface is at a uniform temperature
energy flux leaving a surface is evenly distributed across the surface,
the medium within the enclosure is non-participating.
These assumptions are frequently used in all but the most critical engineering applications.
4/1/13 57
Surface Heat Balance Manager / Processes Outside Surface Heat Balance
"
qLWR qgnd qsky
qair
(65)
where
=Stefan-Boltzmann constant
Linearized radiative heat transfer coefficients are introduced to render the above equation
more compatible with the heat balance formulation,
q "LWR hr , gnd (Tgnd Tsurf ) hr , sky (Tsky Tsurf ) hr ,air (Tair Tsurf ) (67)
where
Fgnd (Tsurf
4
Tgnd
4
)
hr , gnd (68)
Tsurf Tgnd
4/1/13 58
Surface Heat Balance Manager / Processes Outside Surface Heat Balance
Fsky (Tsurf
4
Tsky
4
)
hr , sky (69)
Tsurf Tsky
Fair (Tsurf
4
Tair4 )
hr ,air (70)
Tsurf Tair
The longwave view factors to ground and sky are calculated with the following expressions
(Walton 1983):
where is the tilt angle of the surface. The view factor to the sky is further split between sky
and air radiation by:
The ground surface temperature is assumed to be the same as the air temperature. The final
forms of the radiative heat transfer coefficients are shown here.
Fgnd (Tsurf
4
Tair4 )
hr , gnd (74)
Tsurf Tair
Fsky (Tsurf
4
Tsky
4
)
hr , sky (75)
Tsurf Tsky
Fsky 1 (Tsurf
4
Tair4 )
hr ,air (76)
Tsurf Tair
References
ASHRAE. 1993. 1993 ASHRAE Handbook – Fundamentals. Atlanta: American Society of
Heating, Refrigerating, and Air-Conditioning Engineers, Inc.
th
Chapman, A. J. 1984. Heat Transfer, 4 Edition, New York: Macmillan Publishing Company.
Lienhard, J. H. 1981. A Heat Transfer Textbook, Englewood Cliffs, N.J.: Prentice-Hall, Inc.
McClellan, T. M., and C. O. Pedersen. 1997. Investigation of Outside Heat Balance Models
for Use in a Heat Balance Cooling Load Calculation. ASHRAE Transactions, Vol. 103, Part 2,
pp. 469-484.
Walton, G. N. 1983. Thermal Analysis Research Program Reference Manual. NBSSIR 83-
2655. National Bureau of Standards.
Atmospheric Variation
All buildings are located in the troposphere, the lowest layer of the atmosphere. The
troposphere extends from sea level to an altitude of 11 km. Throughout the troposphere, air
4/1/13 59
Surface Heat Balance Manager / Processes Outside Surface Heat Balance
temperature decreases almost linearly with altitude at a rate of approximately 1°C per 150 m.
Barometric pressure decreases more slowly. Wind speed, on the other hand, increases with
altitude.
Because the atmosphere changes with altitude (defined as height above ground in this case),
tall buildings can experience significant differences in local atmospheric properties between
the ground floor and the top floor. Buildings interact with the atmosphere through convective
heat transfer between the outdoor air and the exterior surfaces of the building envelope, and
through the exchange of air between the outside and inside of the building via infiltration and
ventilation.
Impetus for using this modeling is illustrated in the next table. Using a 70 story (284 meters)
building as an example, the atmospheric variables are significant.
Table 4. Atmospheric Variables at Two Different Altitudes above Ground Level
Tz Tb L H z H b (77)
where
Tz = air temperature at altitude z
Tb = air temperature at the base of the layer, i.e., ground level for the troposphere
L = air temperature gradient, equal to –0.0065 K/m in the troposphere
Hb = offset equal to zero for the troposphere
Hz = geopotential altitude.
The variable Hz is defined by:
Ez
Hz (78)
E z
4/1/13 60
Surface Heat Balance Manager / Processes Outside Surface Heat Balance
where
E = 6,356 km, the radius of the Earth
z = altitude.
For the purpose of modeling buildings in the troposphere, altitude z refers to the height above
ground level, not the height above sea level. The height above ground is calculated as the
height of the centroid, or area-weighted center point, for each zone and surface.
The air temperature at ground level, Tb, is derived from the weather file air temperature by
inverting the equation above:
Ezmet
Tb Tz , met L Hb (79)
E zmet
where
Tz,met = weather file air temperature (measured at the meteorological station)
zmet = height above ground of the air temperature sensor at the meteorological station.
The default value for zmet for air temperature measurement is 1.5 m above ground. This
value can be overridden by using the Site:WeatherStation object.
Local Wind Speed Calculation
Chapter 16 of the Handbook of Fundamentals (ASHRAE 2005). The wind speed measured
at a meteorological station is extrapolated to other altitudes with the equation:
met
z
Vz Vmet met (80)
zmet
where
z = altitude, height above ground
Vz = wind speed at altitude z
= wind speed profile exponent at the site
= wind speed profile boundary layer thickness at the site
zmet = height above ground of the wind speed sensor at the meteorological station
Vmet = wind speed measured at the meteorological station
met = wind speed profile exponent at the meteorological station
met = wind speed profile boundary layer thickness at the meteorological station.
The wind speed profile coefficients , , met, and met, are variables that depend on the
roughness characteristics of the surrounding terrain. Typical values for and are shown in
the following table:
4/1/13 61
Surface Heat Balance Manager / Processes Outside Surface Heat Balance
The terrain types above map to the options in the Terrain field of the Building object. The
Terrain field can be overridden with specific values for and by using the
Site:HeightVariation object.
The default value for zmet for wind speed measurement is 10 m above ground. The default
values for met and met are 0.14 and 270 m, respectively, because most meteorological
stations are located in an open field. These values can be overridden by using the
Site:WeatherStation object.
Outdoor/Exterior Convection
Heat transfer from surface convection is modeled using the classical formulation:
where
Qc = rate of exterior convective heat transfer
hc,ext = exterior convection coefficient
A = surface area
Tsurf = surface temperature
Tair = outdoor air temperature
Substantial research has gone into the formulation of models for estimating the exterior
convection coefficient. Since the 1930's there have been many different methods published
for calculating this coefficient, with much disparity between them (Cole and Sturrock 1977;
Yazdanian and Klems 1994). More recently Palyvos (2008) surveyed correlations cataloging
some 91 different correlations into four categories based on functional form of the model
equation. EnergyPlus therefore offers a wide selection of different methods for determining
values for hc,ext. The selection of model equations for hc,ext can be made at two different levels.
The first is the set of options available in the input object
SurfaceConvectionAlgorithm:Outside that provides a way of broadly selecting which model
equations are applied throughout the model. The input objects
SurfaceProperty:ConvectionCoefficients and
SurfaceProperty:ConvectionCoefficients:MultipleSurface also provide ways of selecting which
model equations or values are applied for specific surfaces. These basic options are
identified by the key used for input and include:
SimpleCombined
TARP
MoWiTT
DOE-2
AdaptiveConvectionAlgorithm
4/1/13 62
Surface Heat Balance Manager / Processes Outside Surface Heat Balance
Note that when the outside environment indicates that it is raining, the exterior surfaces (exposed to wind)
are assumed to be wet. The convection coefficient is set to a very high number (1000) and the outside
temperature used for the surface will be the wet-bulb temperature. (If you choose to report this variable,
you will see 1000 as its value.)
where
h = heat transfer coefficient
Vz = local wind speed calculated at the height above ground of the surface centroid
D, E, F = material roughness coefficients
The roughness correlation is taken from Figure 1, Page 22.4, ASHRAE Handbook of
Fundamentals (ASHRAE 1989). The roughness coefficients are shown in the following table:
Table 6. Roughness Coefficients D, E, and F.
Note that the simple correlation yields a combined convection and radiation heat transfer coefficient.
Radiation to sky, ground, and air is included in the exterior convection coefficient for this algorithm.
All other algorithms yield a convection only heat transfer coefficient. Radiation to sky, ground, and air is
calculated automatically by the program.
4/1/13 63
Surface Heat Balance Manager / Processes Outside Surface Heat Balance
TARP ALGORITHM
TARP, or Thermal Analysis Research Program, is an important predecessor of EnergyPlus
(Walton 1983). Walton developed a comprehensive model for exterior convection by
blending correlations from ASHRAE and flat plate experiments by Sparrow et. al. In older
versions of EnergyPlus, prior to version 6, the “TARP” model was called “Detailed.” The
model was reimplemented in version 6 to use Area and Perimeter values for the group of
surfaces that make up a facade or roof, rather than the single surface being modeled.
Table 7. Nomenclature List of Variables.
4/1/13 64
Surface Heat Balance Manager / Processes Outside Surface Heat Balance
The Detailed, BLAST, and TARP convection models are very similar. In all three models,
convection is split into forced and natural components (Walton 1981). The total convection
coefficient is the sum of these components.
hc h f hn (83)
The forced convection component is based on a correlation by Sparrow, Ramsey, and Mass
(1979):
1/ 2
PV
h f 2.537W f R f z (84)
A
where
or
Leeward is defined as greater than 100 degrees from normal incidence (Walton 1981).
The surface roughness multiplier Rf is based on the ASHRAE graph of surface conductance
(ASHRAE 1981) and may be obtained from the following table:
Table 8. Surface Roughness Multipliers (Walton 1981).
1
h 1.31 T 3 (86)
4/1/13 65
Surface Heat Balance Manager / Processes Outside Surface Heat Balance
For (T < 0.0 AND an upward facing surface) OR (T > 0.0 AND an downward facing
surface) an enhanced convection correlation is used:
1
9.482 T 3
h (87)
7.283 cos
1
1.810 T 3
h (88)
1.382 cos
2 b
A Constant W/(m K(m/s) -
B Constant - -
2 4/3
Ct Turbulent natural W/(m K ) -
convection constant
2
hc Surface exterior W/(m K) -
convective heat transfer
coefficient
Tso Outside surface C/K -
temperature
T Temperature difference C/K -
between the surface
and air
The MoWiTT model is based on measurements taken at the Mobile Window Thermal Test
(MoWiTT) facility (Yazdanian and Klems 1994). The correlation applies to very smooth,
vertical surfaces (e.g. window glass) in low-rise buildings and has the form:
2
1
hc Ct T 3 aV zb
2
(89)
Constants a, b and turbulent natural convection constant Ct are given in Table 10. The
original MoWiTT model has been modified for use in EnergyPlus so that it is sensitive to the
local suface’s wind speed which varies with the height above ground. The original MoWiTT
model was formulated for use with the air velocity at the location of the weather station. As of
Version 7.2, EnergyPlus uses the “a” model coefficients derived by Booten et al. (2012)
rather than the original values from Yazdanian and Klems (1994).
4/1/13 66
Surface Heat Balance Manager / Processes Outside Surface Heat Balance
NOTE: The MoWiTT algorithm may not be appropriate for rough surfaces, high-rise surfaces, or surfaces
that employ movable insulation.
Table 10. MoWiTT Coefficients (Yazdanian and Klems 1994, Booten et al. 2012)
Wind Ct a b
Direction
2 4/3 2 b
(Units) W/m K W/m K(m/s) -
Windward 0.84 3.26 0.89
Leeward 0.84 3.55 0.617
DOE-2 Model
Table 11. Nomenclature List of Variables.
2 b
a Constant W/(m K(m/s) -
b Constant - -
2
hc Surface exterior W/(m K) -
convective heat
transfer coefficient
2
hc,glass Convective heat W/(m K) -
transfer coefficient for
very smooth surfaces
(glass)
2
hn Natural convective W/(m K) -
heat transfer
coefficient
Rf Surface roughness - -
multiplier
Tso Outside surface C/K -
temperature
T Temperature C/K -
difference between
the surface and air,
Angle between the radian -
ground outward
normal and the
surface outward
normal
The DOE-2 convection model is a combination of the MoWiTT and BLAST Detailed
convection models (LBL 1994). The convection coefficient for very smooth surfaces (e.g.
glass) is calculated as:
4/1/13 67
Surface Heat Balance Manager / Processes Outside Surface Heat Balance
2
hc , glass hn2 aVzb (90)
hn is calculated using Equation (87) or Equation (88) . Constants a and b are given in Table
10.
For less smooth surfaces, the convection coefficient is modified according to the equation
4/1/13 68
Surface Heat Balance Manager / Processes Outside Surface Heat Balance
1
PV 2
h f 2.53R f z
A
Sparrow Leeward
Sparrow et al. (1979) conducted flat plate measurements and develop the following
correlation for finite-size flat plates oriented to leeward.
1
2.53 PVz 2
hf Rf
2 A
MoWITT Windward
As discussed above, Yazdanian and Klems (1994) used outdoor laboratory measurements to
develop the following correlation for smooth surfaces oriented to windward. Booten et al.
(2012) developed revised coefficients for use with local surface wind speeds.
2
hc 0.84 T 2.38V 0.89 2
1
3
z
This model equation is for the total film coefficient and includes the natural convection
portion. Therefore it should not be used in conjunction with a second natural convection
model equation.
MoWITT Leeward
Yazdanian and Klems (1994) used outdoor laboratory measurements to develop the following
correlation for smooth surfaces oriented to leeward. Booten et al. (2012) developed revised
coefficients for use with local surface wind speeds.
4/1/13 69
Surface Heat Balance Manager / Processes Outside Surface Heat Balance
2
hc 0.84 T 2.86V 0.617 2
1
3
z
This model equation is for the total film coefficient and includes the natural convection
portion. Therefore it should not be used in conjunction with a second natural convection
model equation.
Blocken
Blocken et al. (2009) developed a set of correlations for windward facing outdoor surfaces
using numerical methods (key: BlockenWindward).
h f 4.6 V100.89
m ; 11.25
h f 5.0 V100.80
m ;11.25 33.75
h f 4.6 V100.84
m ; 33.75 56.25
h f 4.5 V100.81
m ; 56.25 100.0
Where V10m is the air velocity at the location of the weather station and θ is the angle of
incidence between the wind and the surface in degrees. This model is only applicable to
windward surfaces and lacks a natural convection component and therefore cannot be used
on its own but only within the adaptive convection algorithm for the outside face.
Clear
Clear et al. (2003) developed correlations from measurements for horizontal roofs on two
commercial buildings. In EnergyPlus the implementation uses the model for natural
convection plus turbulent forced convection (eq. 8A in the reference) and applies it to the
center point of each surface section that makes up the roof.
k 1 k 4 1
hc 0.15 RaLn3 R f 0.0296 Re x 5 Pr 3
Ln x
Where
x is the distance to the surface centroid from where the wind begins to intersect the roof.
In EnergyPlus this is currently simplified to half the square root of the roof surface.
Area
Ln of overall roof
Perimeter
Gr
ln 1 Lx 2
Re x
is the weighting factor for natural convection (suppressed at
GrLx
1 ln 1
Re 2x
high forced convection rates)
4/1/13 70
Surface Heat Balance Manager / Processes Outside Surface Heat Balance
g 2 Ln 3 T
GrLn is the Grashof number
Tf 2
Vz x
Re x is the Reynolds number at x
Pr is the Prandtl number
This model only claims to be applicable to horizontal roof surfaces so it may not be applicable
to tilted roofs. It combines natural and forced convection and therefore should not be used in
conjunction with yet another natural convection model.
Emmel
Emmel et al. (2007) developed a set of correlations for outdoor surfaces using numerical
methods. The following equations are for vertical surfaces (key: EmmelVertical):
h f 5.15 V100.81
m ; 22.5
h f 3.34 V100.84
m ; 22.5 67.5
h f 4.78 V100.71
m ; 67.5 112.5
h f 4.05 V100.77
m ;112.5 157.5
h f 3.54 V100.76
m ;157.5 180.0
Where V10m is the air velocity at the location of the weather station and θ is the angle of
incidence between the wind and the surface in degrees. The following equations are used for
horizontal (roof) surfaces (key: EmmelRoof):
h f 5.11 V100.78
m ; 22.5
h f 4.60 V100.79
m ; 22.5 67.5
h f 3.67 V100.85
m ; 67.5 90.
Where θ is the angle of incidence between the wind and the longest edge of the roof surface
in degrees.
This model is for all wind directions but lacks a natural convection component. The model
was developed for simple, rectangular low-rise buildings. It is available only within the
adaptive convection algorithm for the outside face
Nusselt Jurges
Perhaps the oldest equation for wind-driven convection was developed by Nusselt and
Jurges (1922). Palyvos (2008) casts their model in simplified form in SI units as:
hc 5.8 3.94 Vz
Where Vz is the wind velocity in m/s, in EnergyPlus that velocity is adjusted for height above
ground using the z axis coordinate of the surface’s centroid and the site wind model. This
model can be applied to all surfaces and the relatively large constant is assumed to represent
the natural convection portion of a total convection coefficient. The model is not sensitive to
wind direction nor surface roughness.
4/1/13 71
Surface Heat Balance Manager / Processes Outside Surface Heat Balance
McAdams
A venerable equation for wind-driven convection was developed by McAdams (1954) which
Palyvos (2008) casts in SI units as:
hc 5.7 3.8 Vz
Where Vz is the wind velocity in m/s that has been adjusted for height above ground using the
z axis coordinate of the surface’s centroid. This model can be applied to all surfaces and the
relatively large constant is assumed to represent the natural convection portion of a total
convection coefficient. The model is not sensitive to wind direction nor surface roughness.
Mitchell
A useful geometric scale based on building volume is used in an equation developed by
Mitchell (1976). The wind-driven convection equation is cast by Palyvos as:
8.6 Vz0.6
hf
L0.4
Where Vz is the wind velocity in m/s that has been adjusted for height above ground using the
z axis coordinate of the surface’s centroid and L is the cube root of the building’s total
volume. EnergyPlus interprets this as the sum of the volume of all the zones in the input file.
Exterior/External Conduction
The conduction term, , can in theory be calculated using a wide variety of heat conduction
qko
formulations. Typically in EnergyPlus, the Conduction Transfer Function (CTF) method is
used. The available models are described in this section: Conduction Through The Walls.
References
ASHRAE. 1981. 1981 ASHRAE Handbook – Fundamentals, Atlanta: American Society of
Heating, Refrigerating, and Air-Conditioning Engineers, Inc.
ASHRAE. 1989. 1989 ASHRAE Handbook – Fundamentals, Atlanta: American Society of
Heating, Refrigerating, and Air-Conditioning Engineers, Inc.
ASHRAE. 1993. 1993 ASHRAE Handbook – Fundamentals, Chapter 3, Heat Transfer, I-P &
S-I Editions, Atlanta: American Society of Heating, Refrigerating, and Air-Conditioning
Engineers, Inc.
ASHRAE. 2001. 2001 ASHRAE Handbook – Fundamentals, Atlanta: American Society of
Heating, Refrigerating, and Air-Conditioning Engineers, Inc.
ASHRAE. 2005. 2005 ASHRAE Handbook – Fundamentals, Chapter 16, Air Flow Around
Buildings, Atlanta: American Society of Heating, Refrigerating, and Air-Conditioning
Engineers, Inc.
Booten, C., N. Kruis, and C. Christensen. 2012. Identifying and Resolving Issues in
EnergyPlus and DOE-2 Window Heat Transfer Calculations. National Renewable Energy
Laboratory. NREL/TP-5500-55787. Golden, CO.
Cole, R. J., and N. S. Sturrock. 1977. The Convective Heat Exchange at the External Surface
of Buildings. Building and Environment, Vol. 12, p. 207.
Ellis, P.G., and P.A. Torcellini. 2005. "Simulating Tall Buildings Using EnergyPlus",
Proceedings of the Ninth International IBPSA Conference, Building Simulation 2005,
Montreal, Canada, August 15-18, 2005.
Lawrence Berkeley Laboratory (LBL). 1994. DOE2.1E-053 source code.
4/1/13 72
Surface Heat Balance Manager / Processes Inside Heat Balance
Sparrow, E. M., J. W. Ramsey, and E. A. Mass. 1979. Effect of Finite Width on Heat Transfer
and Fluid Flow about an Inclined Rectangular Plate. Journal of Heat Transfer, Vol. 101, p.
204.
U.S. Standard Atmosphere. 1976. U.S. Government Printing Office, Washington, D.C.
Walton, G. N. 1981. Passive Solar Extension of the Building Loads Analysis and System
Thermodynamics (BLAST) Program, Technical Report, United States Army Construction
Engineering Research Laboratory, Champaign, IL.
Walton, G. N. 1983. Thermal Analysis Research Program Reference Manual. NBSSIR 83-
2655. National Bureau of Standards.
Yazdanian, M. and J. H. Klems. 1994. Measurement of the Exterior Convective Film
Coefficient for Windows in Low-Rise Buildings. ASHRAE Transactions, Vol. 100, Part 1, p.
1087.
The heart of the heat balance method is the internal heat balance involving the inside faces
of the zone surfaces. This heat balance is generally modeled with four coupled heat transfer
components: 1) conduction through the building element, 2) convection to the air, 3) short
wave radiation absorption and reflectance and 4) longwave radiant interchange. The incident
short wave radiation is from the solar radiation entering the zone through windows and
emittance from internal sources such as lights. The longwave radiation interchange includes
the absorption and emittance of low temperature radiation sources, such as all other zone
surfaces, equipment, and people.
The heat balance on the inside face can be written as follows:
qLWX qSW
qLWS
qki qsol
qconv
0 (92)
where:
qLWX = Net longwave radiant exchange flux between zone surfaces.
= Net short wave radiation flux to surface from lights.
qSW
= Longwave radiation flux from equipment in zone.
qLWS
qki = Conduction flux through the wall.
= Transmitted solar radiation flux absorbed at surface.
qsol
= Convective heat flux to zone air.
qconv
Each of these heat balance components is introduced briefly below.
4/1/13 73
Surface Heat Balance Manager / Processes Inside Heat Balance
Longwave radiation
exchange with other
surfaces in zone
4/1/13 74
Surface Heat Balance Manager / Processes Inside Heat Balance
qi , j Ai Fi , j Ti 4 T j4
4/1/13 75
Surface Heat Balance Manager / Processes Inside Heat Balance
Transmitted Solar
Transmitted solar radiation is also distributed over the surfaces in the zone in a prescribed
manner. It would be possible to calculate the actual position of beam solar radiation, but that
would involve partial surface irradiation, which is inconsistent with the rest of the zone model
that assumes uniform conditions over an entire surface. The current procedures incorporate
a set of prescribed distributions. Since the heat balance approach can deal with any
distribution function, it is possible to change the distribution function if it seems appropriate.
Convection to Zone Air
The convection flux is calculated using the heat transfer coefficients as follows:
hc (Ts Ta )
qconv (89)
The inside convection coefficients (hc) can be calculated using one of many different models.
Currently the implementation uses coefficients based on correlations for natural, mixed, and
forced convection.
Interior Conduction
This contribution to the inside surface heat balance is the wall conduction term, qki shown in
Equation (30). This represents the heat transfer to the inside face of the building element.
Again, a CTF formulation is used to determine this heat flux.
Interior Convection
There are many different modeling options available in EnergyPlus for inside convection
coefficients, hc. There are four different settings to direct how EnergyPlus managers select hc
models during a simulation. There are numerous individual model equations for hc in
EnergyPlus to cover different situations that arise from surface orientations, room airflow
conditions, and heat flow direction. Additionally, in many cases multiple researchers have
developed competing models for the same situations that differ and there is no way to
declare one is better than another. An overall default for the simulation is selected in the
“SurfaceConvectionAlgorithm:Inside” object and can be overridden by selecting a different
option in a zone description. These models are explained in the following sections. In
addition to the correlation choices described below, it is also possible to override the
convection coefficients on the inside of any surface by using the
“SurfaceProperty:ConvectionCoefficients” object in the input file to set the convection
coefficient value on the inside of any surface. The values can be specified directly or with
schedules. Specific details are given in the Input Output Reference document.
Adaptive Convection Algorithm
Beausoleil-Morrison (2000, 2002) developed a methodology for dynamically managing the
selection of hc equations called adaptive convection algorithm. The algorithm is used to
select among the available hc equations for the one that is most appropriate for a given
surface at a given time. As Beausoleil-Morrison notes, the adaptive convection algorithm is
intended to be expanded and altered to reflect different classification schemes and/or new hc
equations. The implementation in EnergyPlus has been modified from the original in the
following ways:
An input mechanism is provided (see the
SurfaceConvectionAlgorithm:Inside:AdapativeModelSelections object) so that the
user can customize the specific selections of hc equations that are applied for
different flow regimes and surface orientations. The changes apply in a general way
to the entire model (but can be overridden by setting surface properties).
To avoid requiring additional user input on the position of ZoneHVAC-type equipment
within a zone, there is no distinction between zones that have convective zone heater
4/1/13 76
Surface Heat Balance Manager / Processes Inside Heat Balance
equipment located underneath the windows and those that have convective heaters
located away from the windows. This applies to the air flow regime associated with
convective zone heaters. Using Beausoleil-Morrison’s terminology, regimes B1 and
B2 are combined into just one B regime.
To avoid requiring additional user input on the position of ZoneHVAC-type equipment
within a zone, there is no distinction between surfaces that are directly blown on the
fan and those that are away from the fan for the air flow regime associated with
mechanical circulation from a zone fan (ZoneHVAC type equipment).
The correlation for horizontal free jet developed by Fisher (1995) is not used. Ceiling
diffuser models are used for all mechanical circulation from central air system. This
decision was made for two reasons: (1) to avoid requiring additional user input on the
position of, and momentum generated by, air terminal units, and (2) because Fisher
(1995) found that the Coanda effect is so significant that in practice a free horizontal
jet is difficult to maintain and mechanical-driven room airflows generally attach to
surfaces and tend to match the flow regime of a ceiling diffuser much more often than
a free jet.
EnergyPlus supports arbitrary geometry so surfaces can be tilted with respect to
vertical or horizontal. Beausoleil-Morrison’s adaptive convection algorithm was
originally structured to use hc equations that have no functional dependence on
surface tilt angle. However, tilted surfaces do perform differently than vertical or
horizontal surface when buoyancy forces are significant. Therefore, the EnergyPlus
implementation expands the structure of the algorithm to include additional
categories for tilted surfaces. The hc equations developed by Walton (1983) are
selected as the defaults for tilted surfaces because they have a functional
dependence on tilt angle.
Fohanno and Polidari (2006) produced a new hc equation for vertical walls inside
buildings with a simple buoyancy flow regime. They used a theoretical approach
based on integral formalism and uniform heat flux (rather than uniform temperature)
that covers both laminar and turbulent flow situations. In EnergyPlus, this model is
selected as the default in place of the model by Alamdari and Hammond (1983) for
vertical walls.
Karadag (2009) produced a new hc equation for ceiling surfaces that are actively
chilled. He used computation fluid dynamics and various sized rooms and
temperature conditions. In EnergyPlus, this model is selected as the default for
surfaces that have active, in-ceiling cooling (in place of the model by Alamdari and
Hammond (1983) for unstable ceilings).
International Standard Organization (ISO) completed Standard 15099-2003 which
includes hc equations for the inside face of windows. EnergyPlus strives to adhere to
formal modeling Standards where possible. Therefore the implementation includes a
larger structure for the adaptive algorithm that includes additional categories for
windows in all flow regimes and ISO 15099-2003 models are used as the default for
windows in natural convection flow regimes. The ISO 15099 model applies to
various tilt angles.
Goldstein and Novosalec (2010) produced new hc equations for forced air situations
with ceiling slot diffusers along perimeters with significant glazing fractions. They
used experiments with full-sized test room. These new equations are selected as the
default for windows, ceilings and floors when there is an active central air system.
Interior mass surfaces are assigned the hc equation that would apply (stable or
unstable) to a horizontal, upward facing surface for each flow regime.
The algorithm switches between forced, mixed, and natural flow regimes by
calculating the Richardson number, Ri = Gr/Re^2, for the zone. Large values of Ri
indicate buoyancy dominates, while small values indicate forced flows dominate. To
distinguish between opposing Zone unit type equipment (with fans) are assumed to
4/1/13 77
Surface Heat Balance Manager / Processes Inside Heat Balance
force air up walls, and central air type equipment (with diffusers) are assumed to
force air down walls.
The adaptive convection algorithm implemented in EnergyPlus for the inside face has a total
of 45 different categories for surfaces and 29 different options for hc equation selections. The
following table summarizes the categories and the default assignments for hc equations. The
individual hc equations are documented below.
Table 1. Inside Convection Categories and Assignments
4/1/13 78
Surface Heat Balance Manager / Processes Inside Heat Balance
4/1/13 79
Surface Heat Balance Manager / Processes Inside Heat Balance
1
h 1.31 T 3
4/1/13 80
Surface Heat Balance Manager / Processes Inside Heat Balance
1
9.482 T 3
h
7.283 cos
Unstable refers to the direction of heat flow and the associated buoyancy relative to the
surfaces. Unstable is when the natural tendency is to enhance flow in the sense that rising
warmer air, or falling cooler air, is free to move away from the surface. This is usually bound
at a minimum of .1 in EnergyPlus. This is a component of the TARP overall algorithm
described below.
Walton Stable Horizontal Or Tilt
Walton (11983) developed the following equation by fitting curves from various sources.
1
1.810 T 3
h
1.382 cos
Stable refers to the direction of heat flow and the associated buoyancy relative to the
surfaces. Stable is when the natural tendency is to retard flow in the sense that rising warmer
air, or falling cooler air, is driven against the surface. This is usually bound at a minimum of .1
in EnergyPlus. This is a component of the TARP overall algorithm described below.
Fisher Pedersen Ceiling Diffuser Walls
Fisher and Pedersen 1997) developed the following equation from laboratory chamber
measurements.
4/1/13 81
Surface Heat Balance Manager / Processes Inside Heat Balance
1
T 5
h 0.6 2
Dh
where,
4A
Dh 2
, hydraulic diameter of horizontal surface, A is area (m ) and P is the perimeter
P
(m) of the entire zone.
1
1
4
6
6
T 6
h 1.4 3
1
1.63 T
Dh
1
1
4
6
6
T 6
h 1.5 3
1
1.23 T
H
where,
H is the characteristic height for the surface. In EnergyPlus this is the zone’s ceiling
height (which could be larger than the height of an individual surface when wall are
subdivided into more than one surface).
Khalifa Eq3 Wall Away From Heat
Khalifa (1989) conducted experiments with test chambers and developed correlations for
certain types of surfaces. One of them, identified as “Equation 3” in original reference, is for
convectively heated zones and applies to the inside surfaces of walls away from the heat
source:
h 2.07 T
0.23
4/1/13 82
Surface Heat Balance Manager / Processes Inside Heat Balance
convectively heated zones and applies to the inside surfaces of ceilings away from the heat
source:
h 2.72 T
0.13
h 1.98 T
0.32
h 2.30 T
0.24
h 3.10 T
0.17
0.308
2.175 T
h
Dh0.076
where,
4A
Dh 2
, hydraulic diameter of horizontal surface, A is area (m ) and P is the perimeter
P
(m) of the entire zone (all of the adjacent floor surfaces if more than one in the zone).
0.293
1.823 T
h
Dh0.076
4/1/13 83
Surface Heat Balance Manager / Processes Inside Heat Balance
where,
4A
Dh 2
, hydraulic diameter of wall surface, A is area (m ) and P is the perimeter (m)
P
of the entire wall (all of the adjacent wall surfaces if more than one along the wall).
1
6
3
6 3
T 4
1 3
3
6
Tsurf TSAT
h 1.5 1.23 T
1
0.199 0.190 ACH 0.8
H T
where,
TSAT is the supply air temperature at the diffuser.
Here the reference temperature is the zone air temperature rather than the diffuser supply air
temperature.
Beausoleil Morrison Mixed Opposing Wall
Beausoleil-Morrison (2000) used blending techniques to combine correlations originally
developed by Alamdari and Hammond (1983) and Fisher and Pedersen (1997) to create the
following correlation is for walls where the flow driving forces from mechanical forces are
opposed by the driving forces from buoyancy.
3
1
3
1
4
6
6 3
1.5 T 1.23 T 13 Tsurf TSAT
0.8
6
T
6
h max 0.8 1.5 1
1.23 T 3
H
0.8 Tsurf TSAT 0.199 0.190ACH 0.8
T
4/1/13 84
Surface Heat Balance Manager / Processes Inside Heat Balance
1
1
5
3 3 3
1
6
6
3
3
T 4
1
3
6
Tsurf TSAT 0.8
h 1.4 1.63 T 3
1
0.159 0.116 ACH
Dh T
Beausoleil Morrison Mixed Stable Ceiling
Beausoleil-Morrison (2000) used blending techniques to combine correlations originally
developed by Alamdari and Hammond (1983) and Fisher and Pedersen (1997) to create the
following correlation is for ceilings where the flow driving forces include both mechanical
forces and thermally Stable buoyancy.
1
1
5
3 3 3
T Tsurf TSAT
h 0.6 0.166 0.484 ACH
0.8
DH T
T 4
1
3
6
Tsurf TSAT 0.8
h 1.4 3
1
1.63 T 0.166 0.484 ACH
Dh T
Fohanno Polidori Vertical Wall
Fohanno and Polidori (2006) developed the following equation for hc for vertical walls under
simple buoyancy flow conditions.
4/1/13 85
Surface Heat Balance Manager / Processes Inside Heat Balance
T 4
1
where,
g f qcH 4
Ra*H Pr f
k f 2f
h 3.1 T
0.22
hi Nu
H
where,
is the thermal conductivity of air, and
H is the height of the window.
The Rayleigh number based on height, RaH , is calculated using,
2 H 3 g c p Tsurf ,i Tair
RaH
Tm , f
where,
is the density of air
g is the acceleration due to gravity,
c p is the specific heat of air,
is the dynamic viscosity of air, and
Tm , f is the mean film temperature in Kelvin given by,
Tm , f Tair
1
4
Tsurf ,i Tair
4/1/13 86
Surface Heat Balance Manager / Processes Inside Heat Balance
There are four cases for the Nusselt correlation that vary by the tilt angle in degrees, , and
are based on heating conditions. For cooling conditions (where Tsurf ,i Tair ) the tilt angle is
complemented so that 180
Case A. 0 15
1
Nu 0.13RaH3
1
e0.72 5
Racv 2.5 105
sin
1 1
Nu 0.13 RaH3 RaCV3 0.56 RaCV sin 4 ; RaH RaCV
1
1
Nu 0.58RaH5 ; RaH 1011
The material properties are evaluated at the mean film temperature. Standard EnergyPlus
psychrometric functions are used for and c p . Thermal conductivity is calculated using,
This correlation depends on the surface temperature of the room-side glazing surface and is
therefore included inside the window heat balance iteration loop.
Goldstein Novoselac Ceiling Diffuser Window
Goldstein and Novoselac (2010) used laboratory chamber measurements to develop
convection correlations for perimeter zones with highly glazed spaces served by overhead
slot-diffuser-based air systems. The following are for bare windows in such spaces.
For WWR<50% with window in upper part of wall:
4/1/13 87
Surface Heat Balance Manager / Processes Inside Heat Balance
V
0.8
h 0.117
L
V
0.8
h 0.093
L
V
0.8
h 0.103
L
Where,
WWR is the window to wall ratio.
L is the length of exterior wall with glazing in the zone.
V is the air system flow rate in m3/s.
Goldstein Novoselac Ceiling Diffuser Walls
Goldstein and Novoselac (2010) used laboratory chamber measurements to develop
convection correlations for perimeter zones with highly glazed spaces served by overhead
slot-diffuser-based air systems. The following are for exterior walls in such spaces.
For walls located below a window
V
0.8
h 0.063
L
V
0.8
h 0.093
L
V
0.8
h 0.048
L
Separate from the above model structure, there are also other comprehensive algorithm
structures which are described below.
4/1/13 88
Surface Heat Balance Manager / Processes Inside Heat Balance
TARP Algorithm
The comprehensive natural convection model, accessed using the keyword “TARP,”
correlates the convective heat transfer coefficient to the surface orientation and the difference
between the surface and zone air temperatures (where T = Surface Temperature - Air
Temperature). The algorithm is taken directly from Walton (1983). Walton derived his
algorithm from ASHRAE literature which can now be found for example in the ASHRAE
Handbook (HoF 2001), Table 5 on p. 3.12, which gives equations for natural convection heat
transfer coefficients in the turbulent range for large, vertical plates and for large, horizontal
plates facing upward when heated (or downward when cooled). A note in the text also gives
an approximation for large, horizontal plates facing downward when heated (or upward when
cooled) recommending that it should be half of the facing upward value. Walton adds a curve
fit as a function of the cosine of the tilt angle to provide intermediate values between vertical
and horizontal. The curve fit values at the extremes match the ASHRAE values very well.
For no temperature difference OR a vertical surface the following correlation is used:
1
h 1.31 T 3 (90)
For (T < 0.0 AND an upward facing surface) OR (T > 0.0 AND an downward facing
surface) an enhanced convection correlation is used:
1
9.482 T 3
h (91)
7.283 cos
1
1.810 T 3
h (92)
1.382 cos
h = 3.076
h = 0.948
4/1/13 89
Surface Heat Balance Manager / Processes Inside Heat Balance
h = 4.040
h = 2.281
h = 3.870
Floor Correlation
exp. data
14
12
correlation
10
h (W/m**2-K)
8
6
4
2
0
0 50 100 150
ACH
For ceilings:
4/1/13 90
Surface Heat Balance Manager / Processes Inside Heat Balance
Ceiling Correlation
exp. data
50
40 correlation
h (W/m**2-K)
30
20
10
0
0 50 100 150
ACH
For Walls:
Wall Correlation
exp. data
25
20 correlation
h (W/m**2-K)
15
10
0
0 50 100 150
ACH
4/1/13 91
Surface Heat Balance Manager / Processes Inside Heat Balance
and a Trombe zone. The rest of the zone heat balance is the same, e.g., transmitted solar,
long-wave radiation between surfaces, etc.
For a vertical cavity, the correlation from ISO 15099 is:
1
NU1=0.0673838 Ra 3 for 5E4 < Ra < 1E6
A
0.272
NU2 0.242 Ra
NU=MAX(NU1,NU2)
where
Nu = Nusselt number
Ra = Rayleigh number
A = aspect ratio of cavity
This is then used in EnergyPlus as follows:
Net convection coefficient from glazing to wall is:
hnet k NU
L
where
k = conductivity of air
L = air gap thickness
Convection coefficient applied to each wall separately and actually used in the zone heat
balance is:
hc 2hnet
References
Alamdari, F. and G.P. Hammond. 1983. Improved data correlations for buoyancy-driven
convection in rooms. Building Services Engineering Research & Technology. Vol. 4, No. 3.
ASHRAE. 1985. 1985 ASHRAE Handbook – Fundamentals, Atlanta: American Society of
Heating, Refrigerating, and Air-Conditioning Engineers, Inc.
ASHRAE. 2001. 2001 ASHRAE Handbook – Fundamentals, Atlanta: American Society of
Heating, Refrigerating, and Air-Conditioning Engineers, Inc.
Awbi, H.B. and A. Hatton. 1999. Natural convection from heated room surfaces. Energy and
Buildings 30 (1999) 233-244.
Beausoleil-Morrison, I. 2000. The adaptive coupling of heat and air flow modeling within
dynamic whole-building simulations. PhD. Thesis. University of Strathclyde, Glasgow, UK.
Ellis, Peter G. 2003. Development and Validation of the Unvented Trombe Wall Model in
EnergyPlus. Master's Thesis, University of Illinois at Urbana-Champaign.
Fisher, D.E. and C.O. Pedersen. 1997. "Convective Heat Transfer in Building Energy and
Thermal Load Calculations", ASHRAE Transactions, Vol. 103, Pt. 2.
4/1/13 92
Surface Heat Balance Manager / Processes Adiabatic Boundary Conditions
Fohanno, S., and G. Polidori. 2006. Modelling of natural convective heat transfer at an
internal surface. Energy and Buildings 38 (2006) 548 - 553
Goldstein, K. and A. Novoselac. 2010. Convective Heat Transfer in Rooms With Ceiling Slot
Diffusers (RP-1416). HVAC&R Research Journal TBD
Karadag, R. 2009. New approach relevant to total heat transfer coefficient including the effect
of radiation and convection at the ceiling in a cooled ceiling room. Applied Thermal
Engineering 29 (2009) 1561-1565
Khalifa AJN. 1989. Heat transfer processes in buildings. Ph.D. Thesis, University of Wales
College of Cardiff, Cardiff, UK.
ISO. 2003. ISO 15099:2003. Thermal performance of windows, doors, and shading devices –
Detailed calculations. International Organization for Standardization.
Walton, G. N. 1983. Thermal Analysis Research Program Reference Manual. NBSSIR 83-
2655. National Bureau of Standards (now NIST). This is documentation for “TARP.”
One of the important technical features of EnergyPlus is that the fundamental zone model
includes infrared (IR) radiation exchange among all surfaces within the zone. A zone can
consist of a single room or multiple rooms having similar thermal characteristics. The other
fundamental characteristic of a zone in EnergyPlus is that the air within the zone is modeled
with the assumption that it is well-stirred, so it is characterized by a single uniform air
temperature. Normally, this corresponds to the desired physical state of the air within a zone.
Practically, this is accomplished with the mixing caused by the air distribution system or even
with the addition of ceiling fans. However, there are situations where the well-stirred model
condition does not fit the physical conditions within a space (e.g., large atria spaces,
displacement ventilation systems, and under floor air distribution systems). In these cases,
the air temperature varies spatially within the zone, but the surfaces in the regions that have
different air temperatures still exchange IR radiation with the other surfaces in the zone, even
those in regions having a different air temperature.
Historically, several approaches have been taken to model such situations. They all involve
unrealistically modifying the convective heat transfer characteristics within a single zone to
produce reasonable total heat transfer for the zone while maintaining the zone IR radiation in
its basic form. The IRT model takes a different, more fundamental approach. The space is
divided into subzones, each having the basic well-stirred air model, but surfaces in these
subzones are able to exchange IR radiation with other surfaces throughout the original
space. Any convective air exchange between subzones is handled using the existing flexible
capabilities within EnergyPlus. In other words, the subzones are standard EnergyPlus zones
but they have been given the capability of allowing IR radiation to be exchanged with
surfaces in adjacent zones.
4/1/13 93
Surface Heat Balance Manager / Processes Infrared Radiation Transfer Material
The same arrangement can be used for simulating a underfloor air distribution system
(UFAD). The two zones represent the lower occupied (mixed) zone and the upper stratified
zone. Since the upper zone is being modeled as a mixed EnergyPlus zone, it is not precisely
the stratified zone concept. However, if a user has concern about having the entire upper part
of the space at a single average temperature, the space could be modeled with two stacked
upper zones. In that case the stratified temperature profile would be established by the
relative size of the mixing from the convective plumes. The user would have to supply those
estimates from external knowledge of the behavior of UFAD systems and plumes.
4/1/13 94
Surface Heat Balance Manager / Processes Infrared Radiation Transfer Material
4/1/13 95
Surface Heat Balance Manager / Processes Infrared Radiation Transfer Material
q13
2
A T14 T34 (95)
where:
2
q is the heat flux in W/m
σ is the Stephan Boltzman constant
2
A is the plate area in m , and
T is the temperature in K.
Equation (95) shows that the presence of a black body surface between a source and a sink
reduces the heat flux by a factor of two. The same result occurs when the IRT surface is
between two zones in EnergyPlus. In that case the adjacent zones behave as black body
cavities at some equivalent temperature. In order to account for this reduction, the IRT area
must be doubled. This can be done without any difficulty in the EnergyPlus radiant exchange
routine because the radiation view factors are determined by an approximate procedure that
is based on the areas of the surfaces. Thus, doubling the surface area of the IRT surface
results in the correct transfer of radiation through the IRT surface. The doubling will occur
automatically in the program as described in the Input Output Reference document.
It should be noted that, because of the black body behavior of the IRT surface, any visible or
solar short wavelength radiation incident on the surface will be absorbed and included with
the long wavelength (IR) exchange with the adjacent zone. No energy will be lost, but zones
with IRT surfaces should not be used in any lighting analyses.
Radiation Transfer Surface Details
Specifying an Infrared Transparent (IRT) surface
The Infrared Transparent (IRT) surface is similar to a resistance-only surface. The idd object
for this type of surface is shown below. The fields indicate that the surface will actually
participate in the transfer of visible and solar radiation by doing a wavelength transformation
and making all short wavelength radiation that is incident on the surface into long wavelength
radiation and having it participate in the long wavelength radiant exchange. The
Material:InfraredTransparent object requires only a name. All other parameters are set
internally.
The Infrared Transparent surface should not participate in a convective/conductive exchange
between the zones it separates. In order to minimize this effect, the
SurfaceProperty:ConvectionCoefficients object must be used. Outside and Inside values for
the surface’s convection coefficients should be on the order of .1. Further examples are
given in the Input Output Reference document.
Behavior Checks
The behavior of multi zones separated with infrared transparent surfaces can be checked
with a simple comparison. Begin with a single zone model as shown below. This model has a
south-facing window, and four walls exposed to wind and sun, and a roof exposed to wind
and sun.
4/1/13 96
Surface Heat Balance Manager / Processes Infrared Radiation Transfer Material
The single zone model will be compared with a stacked three zone model that has zones
separated by interzone infrared transparent surfaces. This model is shown below.
4/1/13 97
Surface Heat Balance Manager / Processes Infrared Radiation Transfer Material
The two upper zones have south facing windows whose total area is the same as the area of
the single window in the single zone model. The top and the sides are again exposed to sun
and wind. The separating surfaces are modeled as IRT surfaces. All zones in both models
are controlled at the same setpoint temperature using purchased air.
The sensible heating results are shown below. The results show the sum of the sensible
cooling load for the three stacked zones and the single zone. It is clear that the IRT surfaces
are very effective in transmitting infrared radiation between the zones. Some small
differences, such as those shown, will occur because of the conversion from short
wavelength solar to infrared through the special IRT dividing surfaces.
4/1/13 98
Surface Heat Balance Manager / Processes Transparent Insulation Material (TIM)
Introduction
The input object “SurfaceControl:MovableInsulation” allows modeling Transparent Insulation
Materials (TIM) that were originally designed for use in solar collector systems, where there
was a need to increase the insulation in the solar collector without dramatically reducing solar
energy transmittance. Transparent Insulation provides both these properties, insulation from
heat loss and transmittance of solar energy. The combination of these properties is
achieved, because Transparent Insulation is a transmitter of short wave radiation but a
barrier to longwave radiation. Therefore short wave solar radiation passes through the
Transparent Insulation and longwave heat radiation is insulated by the transparent insulation.
Incident solar energy falling on the transparent insulation is reflected and re-reflected within
the material and eventually falls on the absorber. In addition, transparent insulation materials
also have increase thermal resistance due to conduction in comparison to standard glass.
Transparent Insulation is now used in the housing industry as a passive solar feature. It is
attached to the walls of houses for insulation and solar energy gains are transmitted to the
house during the right ambient conditions. The walls of the house act as a thermal mass,
absorbing the sunlight at the surface and converting it into heat which is slowly transmitted to
the inside of the house.
4/1/13 99
Surface Heat Balance Manager / Processes Transparent Insulation Material (TIM)
Figure 25. Energy Flows of Opaquely and Transparently Insulated Walls (Wood and Jesch 1993).
While both types of insulation reduce energy losses from the building via conduction through
the building surfaces, transparent insulation allows solar radiation to penetrate deeper into
the surface construction. This increases the construction internal temperature and can result
in heat being conducted into the building under the proper weather conditions. This can be
seen in the lower half of the above figure during a sunny day. The temperature plot shows a
maximum between the transparent insulation and the rest of the surface construction. As a
result, the temperature gradient results in heat transfer from this point into the interior space,
causing a heating effect on the zone. Thus, the advantage of transparent insulation is that,
like opaque insulation, it reduces winter heat transfer losses during low or no solar conditions
and has the possibility of providing heating during sunny winter days. It should be noted that
this same effect in summer could be detrimental to the cooling loads of a building since the
introduction of solar radiation closer to the space will increase the solar heating within the
zone. Most systems counteract this with a shading device or with sophisticated transparent
insulation systems.
Types of Transparent Insulation Materials
Transparent insulation can be classified into four general categories:
Absorber Parallel Covers
Cavity Structures
Absorber Vertical Covers
Quasi-Homogeneous Structures
Cross-sections of each of these types is shown in the figure below. The arrows in these
diagrams indicate solar rays and the path these rays trace as they are transmitted through
the transparent insulation layer. The most advantageous set-up (see absorber-parallel
below) would send most of the rays downward towards the interior of the building while
minimizing the rays that are reflected back to the exterior environment.
4/1/13 100
Surface Heat Balance Manager / Processes Transparent Insulation Material (TIM)
Figure 26. Geometrical Categories of Classification for Transparent Insulation Material (Wood and Jesch
1993).
TIM
TIM
TIM
Incident Solar Gain
Wall
QSM
WALL
WALL
OUTSIDE INSIDE
QSO
Figure 27. Cross Section of TIM and wall, showing energy flow
The total solar gain on any exterior surface is a combination of the absorption of direct and
diffuse solar radiation given by
4/1/13 101
Surface Heat Balance Manager / Processes Transparent Insulation Material (TIM)
Ss
Incident Solar ( I b cos I s Fss I g Fsg ) (97)
S
Where,
= solar absorptance of the surface
= angle of incidence of the sun’s rays
S = area of the surface
Ss = sunlit area of the surface
Ib = intensity of the beam (direct) radiation
Is = intensity of the sky diffuse radiation
Ig = intensity of the beam (direct) radiation
Fss = angle factor between the surface and the sky
Fsg = angle factor between the surface and the ground
Now,
The model for TIM is simplified in that it assumes that absorption of solar radiation takes
place at the inside and outside of the TIM only, not throughout the material. In addition, the
model assumes that the solar radiation absorbed during the first pass through the TIM affects
the outside surface of the TIM while the solar radiation reflected at the outer wall surface that
gets absorbed during the back reflection will affect the inside TIM surface (which is also the
outside surface of the wall). Thus, the heat absorbed at the outside of the TIM is as shown in
Equation (96).
The heat absorbed at the inside of the TIM/outside of the wall includes two components. The
first component is the amount of solar that is transmitted through the TIM and absorbed at the
inside of the wall. This is characterized by the following equation:
First pass solar absorbed by wall TIM Incident Solar wall
The amount of solar absorbed by the TIM and aggregated at the inside surface of the TIM
(outside wall surface) is:
Amount of back reflection absorbed byTIM TIM Incident Solar 1 wall TIM (99)
The heat absorbed at the interface between the wall and the TIM includes both of these
components. Thus, QSO is equal to:
QSO TIM Incident Solar wall 1 wall TIM (100)
Substituting the definition for QSM into this equation and rearranging results in:
QSM
QSO TIM wall 1 wall TIM
TIM
4/1/13 102
Surface Heat Balance Manager / Processes Transparent Insulation Material (TIM)
QSO TIM QSM wall 1 wall (101)
TIM
Where,
QSM = Short wave radiant flux absorbed at surface of Movable Insulation
QSO = Short wave radiant flux absorbed at surface of Wall.
TIM = Absorptance of TIM
TIM = Transmittance of TIM.
WALL = Absorptance of Wall.
WALL = Reflectance of Wall surface
Following is the FORTRAN Code used in the HeatBalanceSurfaceManager module, to
determine the short wave radiation absorbed on outside of movable insulation and the short
wave radiation absorbed on outside of opaque surface of the wall.
IF (Surface(SurfNum)%MaterialMovInsulExt.GT.0) &
CALL EvalOutsideMovableInsulation(SurfNum,HMovInsul,RoughIndexMovInsul,AbsExt)
IF (HMovInsul > 0) THEN ! Movable outside insulation in place
QRadSWOutMvIns(SurfNum) = QRadSWOutAbs(SurfNum)*AbsExt &
/Material(Construct(ConstrNum)%LayerPoint(1))%AbsorpSolar
! For Transparent Insulation
QRadSWOutAbs(SurfNum) = Material(Surface(SurfNum)%MaterialMovInsulExt)%Trans &
*QRadSWOutMvIns(SurfNum)* &
( (Material(Construct(ConstrNum)%LayerPoint(1))%AbsorpSolar/AbsExt) &
+(1-Material(Construct(ConstrNum)%LayerPoint(1))%AbsorpSolar) )
4/1/13 103
Surface Heat Balance Manager / Processes Transparent Insulation Material (TIM)
The following two tables shows data for two series of runs. The first “summer table”
illustrates the execution of a summer design day. The second “winter table” shows winter
conditions with clearness=0 (the typical default for a winter design day) and clearness=1 (to
illustrate solar radiation with other winter conditions). Test cases included no movable
insulation, moveable opaque insulation, and TIM on the exterior (south wall unless otherwise
noted). Savings reported are heating and cooling loads for the design days only. The results
showed that the TIM model was performing reasonably well and was producing results that
were within expectations.
4/1/13 104
Surface Heat Balance Manager / Processes Surface Heat Balance With Moveable Insulation
References
P.O. Braun, A. Goetzberger, J. Schmid, and W.Stahl. Transparent Insulation of Building
Facades- Steps from Research to Commercial applications, Fraunhofer Institute for Solar
Energy Systems, Oltmannsstrasse 22, D-7800 Freiburg, Germany.
Thermotropic materials and Systems for Overheating Protection.
http://www.ise.fhg.de/Projects/Solbuild/materials.html
Robert Hausner. Arbeitsgemeinschaft Erneuerbare energie, Transparent Insulation- Areas of
Application, Society for Renewable Energy. http://www..aee.at/verz/english/tin.html
Werner J.Platzer. Transparent Insulation materials: a review, Fraunhofer Institute for Solar
Energy Systems, Oltmannsstr. 5, D-79100 Freiburg, Germany.
Volker Wittwer. The use of Transparent Insulation Materials and Optical Switching Layers in
Window Systems, Fraunhofer Institute for Solar Energy Systems, Oltmannsstr. 5, D-79100
Freiburg, Germany.
M. Wood and L.F. Jesch. 1993. Transparent insulation technology: a technical note, Ambient
Press Limited.
Façade Modules with back-ventilated Transparent Insulation- Research and Development
toward Series Application. http://www.ise.fhg.de/Projects/development99/art4.html
Two 0-Energy Houses, http://www.smartarch.nl/smartgrid/items/oo5_chur.htm
Advanced Building Technologies – Transparent Insulation Materials ( TIM ).
http://www.enermodal.com/advancedtech/transp.html
Transparent Insulation, http://www.esv.or.at/service/info-material/diverse/twd/index_e.htm
G. Verbeeck, H. Hens. Transparent Insulation: an alternative solution for summer discomfort.
Die neue Transparenz: Warmedamm-Verbund- system StoTherm Solar.
E.Lindauer, H.Leonhardt. Brauchwasservorerwarmmung mit transparent gedammten
Bauteilen ( Hybridsystem ), Fraunhofer- Institut fur Bauphysik.
In contrast to the internal surface heat balance that treats all surfaces simultaneously, the
external thermal balance for each surface is performed independent of all other surfaces.
This implies that there is no direct interaction between the individual surfaces.
4/1/13 105
Surface Heat Balance Manager / Processes Surface Heat Balance With Moveable Insulation
TARP includes four possible representations for the basic outside surface heat balance. The
first two depend on which of the optimal surface conductance algorithms the user selects.
The simple outside surface conductance that includes both the convective and thermal
interchange between the surface and the environment in a single coefficient, is represented
by the thermal network in Figure 28. Equation (102) can also be expressed as:
The detailed outside surface conductance model considers convection and radiant
interchange with the sky and with the ground as separate factors. Its use in the outside
thermal balance is shown in Figure 29. In this case, equation (102) can be expanded to give
KOPt +Y0 TI t -X 0 TOt + HA Ta -TOt +HS Ts -TOt +HG Tg -TOt +QSO = 0 (105)
The third and fourth representations occur when the outside surface has been covered with
movable insulation. The insulation has a conductance of UM. The thermal network in Figure
30 represents this case. The insulation must be mass-less because it is not generally
possible to perform a correct thermal balance at the juncture of two surfaces each modeled
by CTF.
The equation for the thermal balance between the surface and the insulation is
Depending on whether or not the detailed or simple algorithm for surface conductance is
being used, there are two expressions for TM, the outside temperature of the insulation. For
the simple conductance:
4/1/13 106
Surface Heat Balance Manager / Processes Surface Heat Balance With Moveable Insulation
In this case the values of HA, HS and HG must be found by using an estimated value of TM
in place of TO.
QSO
KOP
TA TO TI
1/Ho 1/Yo
TS QSO
1/HS
KOP
TA TO TI
1/HA 1/Yo
1/HG
TG
QSO
KOP
TM TO TI
1/UM 1/Yo
4/1/13 107
Surface Heat Balance Manager / Processes Surface Heat Balance With Moveable Insulation
effecting the value of TOt . When Y0 is large, TO and TI can so strongly be coupled that
separate outside and inside heat balances do not work because the environment and zone
temperatures have negligible influence on the heat balances. The TARP uses the inside
surface heat balance to couple TOt with TZ and TR. These two temperatures are less
strongly influenced by TO and allow a reasonable heat balance. On the first heat balance
iteration, TZ and TR are the values at time t-1. The user may optionally require that TOt be
recomputed with every iteration of TIt . In this case TZ and TR have values from the previous
iteration and a true simultaneous solution is achieved. In most conventional constructions,
recomputing TOt does not significantly change the computed zone loads and temperatures.
The inside surface heat balance is given by
The surface heat balances can be combined in eight ways according to conditions for
calculations of the outside surface temperature
Y0
F1 (112)
Z 0 HI HR
UM
F2
UM HO
(113)
UM
F3
UM HA HS HG
(114)
4/1/13 108
Surface Heat Balance Manager / Processes Surface Heat Balance With Moveable Insulation
KOPt QSO F1 KIPt +QSI+HITZ+HR TR +F3 QSM+HATa +HSTs +HG Tg
TO t = (122)
X 0 +UM-F3 UM-F1 Y0
4/1/13 109
Surface Heat Balance Manager / Processes Surface Heat Balance With Moveable Insulation
4/1/13 110
Surface Heat Balance Manager / Processes Surface Heat Balance With Moveable Insulation
4/1/13 111
Surface Heat Balance Manager / Processes Surface Heat Balance With Moveable Insulation
coefficient
2
HExtSurf(SurfNum) Outside Convection HO, HA W/m K Overall
Coefficient outside
surface
conductance
2
HGround Radiant exchange HG W/m K Radiative
(linearized) conductance
coefficient (outside
surface to
ground
temperature
2
HmovInsul Conductance or "h" UM W/m K Conductance
value of movable of Movable
insulation insulation
2
HSky Radiant exchange HS W/m K Radiative
(linearized) conductance
coefficient (outside
surface to
sky radiant
temperature
MAT(ZoneNum) Zone temperature TZ C Temperature
of zone air
2
NetLWRadToSurf(SurfNum) Net interior HR*TR W/m Net surface
longwave radiation to surface
to a surface from radiant
other surfaces exchange
2
QRadSWInAbs(SurfNum) Short-wave radiation QSI W/m Short wave
absorbed on inside radiant flux
of opaque surface absorbed at
inside of
surface
2
QRadSWOutAbs(SurfNum) Short wave radiation QSO W/m Short wave
absorbed on outside radiant flux
opaque surface absorbed at
outside of
surface
2
QRadSWOutMvIns(SurfNum) Short wave radiation QSM W/m Short wave
absorbed on outside radiant flux
of movable absorbed at
insulation surface of
movable
insulation
2
QRadThermInAbs(SurfNum) Thermal Radiation W/m Longwave
absorbed on inside radiant flux
surfaces from internal
gains
SkyTemp Sky temperature Ts C Sky temp
TempExt Exterior surface TM, Ta C Temperature
temperature or of external
4/1/13 112
Surface Heat Balance Manager / Processes Ground Heat Transfer Calculations using C and F Factor Constructions
References
Walton, G.N. 1983. “The Thermal Analysis Research Program Reference Manual Program
(TARP)”, National Bureau of Standards (now National Institute of Standards and
Technology).
Building energy code and standards like ASHRAE 90.1, 90.2 and California Title 24 require
the underground wall constructions and slabs-on-grade or underground floors not to exceed
certain maximum values of C-factor and F-factor, which do not specify detailed layer-by-layer
materials for the constructions. If using the normal approach (layer by layer) of ground
constructions with EnergyPlus, users would need to create a pseudo wall or floor construction
to match the thermal performance such as thermal mass effect and U-factor, and rely on the
EnergyPlus Basement and Slabs tools to generate the monthly ground temperatures.
A simplified approach is introduced to create equivalent constructions and model the ground
heat transfer through underground walls and ground floors for the building energy code
compliance calculations. The approach is to create constructions based on the user defined
C or F factor with two layers: one concrete layer (0.15 m thick) with thermal mass, and one
fictitious insulation layer with no thermal mass.
Three objects are used in the C and F factor calculations:
Construction:CfactorUndergroundWall
Construction:FfactorGroundFloor
Site:GroundTemperature:FCfactorMethod
Site:GroundTemperature:FCfactorMethod is used only by the underground walls or slabs-on-
grade or underground floors defined with C-factor (Construction:CfactorUndergroundWall)
and F-factor (Construction:FfactorGroundFloor) method for code compliance calculations
where detailed construction layers are unknown. Only one such ground temperature object
can be included. The monthly ground temperatures for this object are close to the monthly
outside air temperatures delayed by three months. If user does not input this object in the IDF
file, it will be defaulted to the 0.5m set of monthly ground temperatures from the weather file if
they are available.
Detailed description of the three objects is in the Input Output Reference document. The
following section describes how the equivalent material layers are created based on the C-
factor and F-factor.
Slab-on-grade and Underground Floors Defined with F-factors
The steady state heat transfer through the floor is calculated as,
4/1/13 113
Surface Heat Balance Manager / Processes Ground Heat Transfer Calculations using C and F Factor Constructions
Where,
Q is the steady state heat transfer through the floor in Watt
2
Area is the area of the floor in m
Ueff is the effective heat transfer coefficient including the floor construction, the soil, and the
thermal resistance of the interior and exterior air films.
Tair,in is the indoor air temperature in °C
Tair,out is the outside air temperature in °C
Pexp is the exposed perimeter of the floor in m
F-Factor is the heat transfer through the floor, induced by a unit temperature difference
between the outside and inside air temperature, on the per linear length of the exposed
perimeter of the floor. The unit is W/m·K.
Therefore,
Where,
2
Reff is the effective thermal resistance in m ·K/W, including the soil and the floor construction
Rfilm,in and Rfilm, out are the air film resistance of the inside and outside surfaces, respectively.
2
The outside air film resistance Rfilm, out = 0.03 m ·K/W. The inside air film resistance Rfilm,in =
2 2 2
0.125 m ·K/W, which is the average of the 0.14 m ·K/W for heat flow up and 0.11 m ·K/W for
heat flow down.
Approximate the thermal mass of the floor construction with a 6-inch (0.15 m) heavy
concrete, and use a fictitious insulation layer with no thermal mass to match the thermal
resistance of the construction.
We have,
Where,
2
Rfic is the thermal resistance of the fictitious insulation layer in m ·K/W.
2
Rcon is the thermal resistance of the concrete layer in m ·K/W.
Properties of the concrete layer are:
Thickness = 0.15 m
Finally,
4/1/13 114
Surface Heat Balance Manager / Processes Ground Heat Transfer Calculations using C and F Factor Constructions
For slabs-on-grade as shown in the figure above, the exposed perimeter is (2A + C) for the
Dining area, and (2B + C) for the Kitchen area. For underground floors with no exposed
2 2
perimeter, the Reff can be assumed a big value such as 1000 hr·ft ·°F/Btu (177 m ·K/W).
Underground Walls Defined with C-Factors
The steady state heat transfer through the underground wall is calculated as,
Where,
Q is the heat transfer through the wall in Watt
Area is the area of the wall in m2
Ueff is the effective heat transfer coefficient including the floor construction, the soil, and the
thermal resistance of the interior and exterior air films.
2
Reff is the effective thermal resistance in m ·K/W, including the soil and the wall construction
Rfilm,in and Rfilm, out are the air film resistance of the inside and outside surfaces, respectively.
C-factor is the time rate of steady-state heat flow through unit area of the construction,
induced by a unit temperature difference between the body surfaces. The C-Factor unit is
2
W/m ·K. The C-factor does not include soil or air films.
Rsoil is the effective R-value of the soil. Reference values from Table C6.10.1 of the SI version
of the ASHRAE Standard 90.1-2010 are as follows:
4/1/13 115
Surface Heat Balance Manager / Processes Ground Heat Transfer Calculations using C and F Factor Constructions
2
A fairly good linear regression (R = 0.9967) for the above data is,
Rsoil = 0.0607 + 0.3479 * Depth
Approximate the thermal mass of the wall construction with a 6-inch (0.15 m) heavy concrete,
and use a fictitious insulation layer with no thermal mass to match the thermal resistance of
the construction. Then we have the thermal resistance of the insulation layer,
Where,
Rfic is the thermal resistance of the fictitious insulation layer
2
Rcon is the thermal resistance of the concrete layer in m ·K/W.
Properties of the concrete layer are:
Thickness = 0.15 m
4/1/13 116
Advanced Surface Concepts Exterior Naturally Vented Cavity
qsol qLWR
, Env qconv
, Env qLWR
,cav qconv
,cav qsource
0 (123)
4/1/13 117
Advanced Surface Concepts Exterior Naturally Vented Cavity
where:
qsol is absorbed direct and diffuse solar (short wavelength) radiation heat flux.
qLWR , Env is net long wavelength (thermal) radiation flux exchange with the air and
surroundings.
, Env = surface convection flux exchange with outside air.
qconv
,cav is net long wavelength (thermal) radiation flux exchange with the outside face of the
qLWR
underlying surface(s).
,cav = surface convection flux exchange with cavity air.
qconv
qsource is a source/sink term that accounts for energy exported out of the control volume
when the baffle is a hybrid device such as a photovoltaic panel.
All terms are positive for net flux to the baffle. Each of these heat balance components is
introduced briefly below.
External SW Radiation
qsol is calculated using procedures presented elsewhere in this manual and includes both
direct and diffuse incident solar radiation absorbed by the surface face. This is influenced by
location, surface facing angle and tilt, shading surfaces, surface face material properties,
weather conditions, etc. The baffle blocks all shortwave radiation from reaching the
underlying surface.
External LW Radiation
qLWR , Env is a standard radiation exchange formulation between the surface, the sky, the
ground, and the atmosphere. The radiation heat flux is calculated from the surface
absorptivity, surface temperature, sky, air, and ground temperatures, and sky and ground
view factors. Radiation is modeled using linearized coefficients. The baffle blocks all
longwave radiation.
External Convection
, Env is modeled using the classical formulation: qconv
qconv = hco(Tair - To) where hco, is the
convection coefficient. The hco is treated in the same way as an outside face with
ExteriorEnvironment conditions. In addition, when it is raining outside, we assume the baffle
gets wet and model the enhanced surface heat transfer using a large value for hco .
Cavity LW Radiation
,cav is a standard radiation exchange formulation between the baffle surface and the
qLWR
underlying heat transfer surface located across the cavity. Radiation is modeled using
linearized coefficients.
Cavity Convection
,cav is modeled using the classical formulation: qconv
qconv = hcp(Tair - To) where hcp, is the
convection coefficient. The value for hcp is obtained from correlations used for window gaps
from ISO (2003) standard 15099.
Substituting models into (113) and solving for Ts ,baff yields the following equation:
4/1/13 118
Advanced Surface Concepts Exterior Naturally Vented Cavity
Ts ,baff
I h
s T
co amb
hr ,atmTamb hr ,skyTsky hr , gnd Tamb hr ,cavTso hc ,cav Ta ,cav qsource
h hr ,air hr , sky hr , gnd hr ,cav hc ,cav
(124)
co
where,
I s is the incident solar radiation of all types [W/m2],
is the solar absorptivity of the baffle [dimensionless],
hr ,atm is the linearized radiation coefficient for the surrounding atmosphere [W/m2·K],
Tamb is the outdoor drybulb from the weather data, also assumed for ground surface [ºC],
hr , sky is the linearized radiation coefficient for the sky [W/m2·K],
Tsky is the effective sky temperature [ºC],
hr , gnd is the linearized radiation coefficient for the ground [W/m2·K],
hr ,cav is the linearized radiation coefficient for the underlying surface [W/m2·K],
Tso is the temperature of the outside face of the underlying heat transfer surface [ºC],
hco is the convection coefficient for the outdoor environment [W/m2·K],
hc ,cav is the convection coefficient for the surfaces facing the plenum [W/m2·K], and
Ta ,cav is the drybulb temperature for air in the cavity [ºC].
4/1/13 119
Advanced Surface Concepts Exterior Naturally Vented Cavity
where,
Q vent is the net rate of energy added from natural ventilation – where outdoor ambient air
exchanges with the cavity air.
Q co is the net rate of energy added by surface convection heat transfer with the underlying
surface.
Q c,baff is the net rate of energy added by surface convection heat transfer with the collector.
And substituting into (125) yields the following equation:
Ta ,cav
hc ,cav ATso m vent c pTamb hc ,cav ATs ,baff
h c ,cav A m vent c p hc ,cav A
where,
m vent is the air mass flow from natural forces [kg/s]
Modeling natural ventilation air exchanges in a general way is challenging. Simplistic
engineering models are used to model m vent resulting from natural buoyancy and wind
forces. Reasoning that the configuration is similar to single-side natural ventilation, we elect
to use correlations for natural ventilation presented as equations (29) and (30) in Chapter 26.
of ASHRAE HOF (2001).
m vent V tot
where,
is the density of air [kg/m3], and
V
Vtot wind Vthermal is the total volumetric flow rate of air ventilating in and out of the
cavity.
V wind Cv AinU
V thermal C D Ain 2 g H NPL Ta ,cav Tamb / Ta ,cav (if Ta ,cav Tamb )
V thermal C D Ain 2 g H NPL Tamb Ta ,cav / Tamb (if Tamb Ta ,cav and baffle is vertical)
Cv is the effectiveness of the openings that depends on opening geometry and the
orientation with respect to the wind. ASHRAE HoF (2001) indicates values ranging from 0.25
to 0.6. This value is available for user input.
C D is the discharge coefficient for the opening and depends on opening geometry. This
value is available for user input.
Mass continuity arguments lead to modeling the area of the openings as one half of the total
area of the openings, so we have:
4/1/13 120
Advanced Surface Concepts Exterior Naturally Vented Cavity
A
Ain
2
situation.
Underlying Heat Transfer Surface
The exterior baffle and cavity are applied to the outside of a heat transfer surface. This
surface is modeled using the usual EnergyPlus methods for handling heat capacity and
transients – typically the CTF method. These native EnergyPlus heat balance routines are
used to calculate Tso . The exterior baffle and cavity system is coupled to the underlying
surface using the SurfaceProperty:OtherSideConditionsModel mechanism. The exterior
naturally vented cavity model provides values for hr ,cav , Ts ,baff , hc ,cav , and Ta ,cav for use with
the heat balance model calculations for the outside face of the underlying surface (described
elsewhere in this manual).
Solar and Shading Calculations
The exterior vented cavity model uses standard EnergyPlus surfaces in order to take
advantage of the detailed solar and shading calculations. Solar radiation incident on the
surface includes beam and diffuse radiation, as well as radiation reflected from the ground
and adjacent surfaces. Shading of the collector by other surfaces, such as nearby buildings
or trees, is also taken into account.
Local Wind Speed Calculations
The outdoor wind speed affects terms used in modeling. The wind speed in the weather file
is assumed to be measured at a meteorological station located in an open field at a height of
10 m. To adjust for different terrain at the building site and differences in the height of
building surfaces, the local wind speed is calculated for each surface.
The wind speed is modified from the measured meteorological wind speed by the equation
(ASHRAE 2001):
amet
a
z
U Vmet met (126)
zmet
where z is the height of the centroid of the system, zmet is the height of the standard
meteorological wind speed measurement, and a and are terrain-dependent coefficients.
is the boundary layer thickness for the given terrain type. The values of a and are shown in
the following tables:
4/1/13 121
Advanced Surface Concepts Exterior Naturally Vented Cavity
Tair Tsurf
hc
qconv
First, hco is the convection coefficient for the baffle surface facing the outdoors. It is modeled
in exactly the same way as elsewhere in EnergyPlus and will depend on the user setting for
Outside Convection Algorithm – Outside Surface Heat Balance entry elsewhere in this
document.
Second, hc ,cav is the convection coefficient for baffle surfaces facing the cavity. This
coefficient is applied to both the baffle and the underlying surface. The convection coefficient
is modeled in the same way used in EnergyPlus to model air gaps in windows. These
correlations vary by Rayleigh number and surface tilt and are based on the work of various
research including Hollands et. al., Elsherbiny et. al., Wright, and Arnold. The formulations
are documented in ISO (2003) standard 15099. The routines were adapted from Subroutine
NusseltNumber in WindowManager.f90 (by F. Winkelmann), which itself was derived from
Window5 subroutine “nusselt”.
Radiation Coefficients
Exterior vented cavity modeling requires calculating up to four different linearized coefficients
for radiation heat transfer. Whereas radiation calculations usually use temperature raised to
the fourth power, this greatly complicates solving heat balance equations for a single
temperature. Linearized radiation coefficients have the same units and are used in the same
manner as surface convection coefficients and introduce very little error for the temperature
levels involved.
The radiation coefficient, hr ,cav , is used to model thermal radiation between the collector
surface and the outside face of the underlying heat transfer surface. We assume a view
factor of unity. It is calculated using:
where,
4/1/13 122
Advanced Surface Concepts Green Roof Model (EcoRoof)
Overview
The input object Material:RoofVegetation provides a model for green roofs (aka ecoroofs or
vegetated roofs) that are becoming increasingly common for both new and retrofit buildings.
There is widespread recognition and a growing literature of measured data that suggest
green roofs can reduce building energy consumption. Currently, however, there are few
design tools available to assist developers and architects in assessing the likely magnitude of
energy savings associated with various implementation options (e.g., soil type/depth,
irrigation options, plant type). As a result there is a significant need for a quantitative and
physically-based building energy simulation tool that represents the effects of green roof
constructions. Such a tool would facilitate more rapid spread of green roof technologies and
make it possible to account for green roof benefits in state energy codes and related energy
efficiency standards such as LEED.
In response to the need for green roof design tools a computational model of the heat
transfer processes involved on a vegetated roof has been developed. This model accounts
for:
long wave and short wave radiative exchange within the plant canopy,
plant canopy effects on convective heat transfer,
evapotranspiration from the soil and plants, and
heat conduction (and storage) in the soil layer
The ability to track moisture-dependent thermal properties is not implemented yet due to
stability issues in the CTF scheme, but is under development for use with the finite difference
solution scheme made available in EnergyPlus starting in version 2.
As implemented in EnergyPlus the green roof module allows the user to specify “ecoroof” as
the outer layer of a rooftop construction using a “Material:RoofVegetation” object. The user
can then specify various aspects of the green roof construction including growing media
depth, thermal properties, plant canopy density, plant height, stomatal conductance (ability to
transpire moisture), and soil moisture conditions (including irrigation).
The model formulation includes the following:
simplified moisture balance that allows precipitation, irrigation, and moisture transport
between two soil layers (top and root zone).
4/1/13 123
Advanced Surface Concepts Green Roof Model (EcoRoof)
soil and plant canopy energy balance based on the Army Corps of Engineers’ FASST
vegetation models (Frankenstein and Koenig), drawing heavily from BATS
(Dickenson et al.) and SiB (Sellers et al.).
soil surface (Tg) and foliage (Tf) temperature equations are solved simultaneously
each time step, inverting the CTF to extract heat flux information for the energy
balance calculation.
The detailed energy balance analysis and resulting equations, being rather complicated, are
summarized here. The interested reader is referred to the FASST documentation cited herein
for the complete development. The end result is a set of two simultaneous equations for
temperature—one for the soil surface and the other for the foliage.
Green Roof Model Description
As with a traditional roof, the energy balance of an green roof is dominated by radiative
forcing from the sun. This solar radiation is balanced by sensible (convection) and latent
(evaporative) heat flux from soil and plant surfaces combined with conduction of heat into the
soil substrate. This energy balance is illustrated in Figure 34. The variables introduced in this
figure are defined in the equations that follow.
The energy budget analysis follows the Fast All Season Soil Strength (FASST) model
developed by Frankenstein and Koenig for the US Army Corps of Engineers. FASST was
developed, in part, to determine the ability of soils to support manned and unmanned
vehicles and personnel movement. In order to accomplish this, however, FASST tracks the
energy and moisture balance (including ice and snow) within a vegetated soil. It is a one-
dimensional model that draws heavily from other plant canopy models including BATS
(Dickinson et al.) and SiB (Sellers et al.). We have implemented FASST here with only a few
modifications to adapt it for use with a relatively thin soil layer. The sign convention used
assumes all heat fluxes are positive when energy is absorbed into the layer.
In the following discussion this energy budget is divided into a budget for the foliage layer (Ff)
and a budget for the ground surface (Fg). The various parameterizations for latent and
sensible heat flux are described in some detail and then the equation set is reduced to the
4/1/13 124
Advanced Surface Concepts Green Roof Model (EcoRoof)
simultaneous solution of two equations involving the temperatures of the foliage and ground
surface.
Energy budget in the foliage layer
The foliage energy balance is given by:
f g f
Ff f I S (1 f ) f Iir f T f 4
1
Tg 4 Tf 4 H f L f
In addition to convective and sensible heat transfer this equation accounts for both the short
and longwave radiation absorbed by the vegetation, including the effects of multiple
reflections. The sensible and latent heat flux terms (Hf and Lf) are somewhat complicated and
therefore discussed in some detail below.
Sensible heat flux in the foliage layer
The sensible heat transfer between the leaf surface and near-canopy air (Hf) is influenced by
the temperature difference between them, wind speed, and Leaf Area Index (LAI). The Leaf
Area Index is the dimensionless ratio of the projected leaf area for a unit ground area (Oke).
In contrast fractional vegetative cover (f) is the ratio of shaded ground surface to total
ground surface area. The sensible heat flux is given by:
In this equation the constant 1.1 accounts for heat transfer from the stems, twigs and limbs
(Deardorff). The properties of air near the foliage are modeled using the average from the
foliage and instrument conditions:
af 0.5( a f )
where a is the density of air at the instrument height and f is the density of air at the leaf
temperature. The air temperature within the foliage is estimated by:
where, Ta is the air temperature at the instrument height in Kelvin Tf, is leaf temperature in
Kelvin and Tg, is the ground surface temperature in Kelvin. The foliage wind speed is
estimated as:
Here W is the larger of 2.0 m/s or the actual wind speed above the canopy (Hughes et al.)
f
and C hn is the transfer coefficient at near-neutral atmospheric stability conditions:
2
Z Z
C K v ln a f d
f 2
hn
Zo
where Kv, is von Karmen’s constant (0.4), Za is the instrument height, Zd is the zero
displacement height in meters (height above soil within which the wind speed is effectively
f
zero), and Z o is the foliage roughness length scale (m). The formulations for zero
displacement height, roughness length are based on Balick et al.:
4/1/13 125
Advanced Surface Concepts Green Roof Model (EcoRoof)
Z d 0.701Z 0.979
f
Z o 0.131Z 0.997
f
0.3(m / s )
C f 0.01* 1
Waf (m / s )
rs ,min
rs f1 f 2 f 3
LAI
Here, rs,min is the minimum stomatal resistance. The actual stomatal resistance at any time is
proportional to this minimum resistance and inversely proportional to LAI. The stomatal
resistance is further modified by fractional multiplying factors that relate to incoming solar
radiation and atmospheric moisture. As found in Frankenstein and Koenig the inverses of the
multiplying factors f1, f2, and f3 are given by:
1 0.004* I s 0.005
min 1,
f1 0.81*(0.004* I s 1)
0 when r
1
r
f2 when r max
max r
1
exp g d (e f , sat ea
f3
Here, r, is the residual moisture content (defined as the amount of moisture in soil when
plants begin to wilt), max is the maximum moisture content (defined as the maximum amount
of moisture a particular type of soil can hold and above which run off occurs), and is the
average soil moisture in the root zone. The residual moisture content is typically around 0.01
3 3
m /m (Frankenstein and Koenig). The maximum moisture content depends upon the soil, but
3 3
generally varies from 0.3 to 0.6 m /m (Guymon et al.). In the expression for f3, gd is a plant
specific characteristic that is only non-zero for trees, ef,sat is the saturated vapor pressure at
the leaf temperature, and ea is the air vapor pressure.
Resistance to moisture exchange offered by the boundary layer formed on the leaf surface is
known as aerodynamic resistance. It is measured in units of (s/m) and is influenced by wind
speed, surface roughness and stability of the atmosphere (Oke). It is formulated as:
4/1/13 126
Advanced Surface Concepts Green Roof Model (EcoRoof)
1
ra
c f Waf
The combined effect of aerodynamic and stomatal resistances to vapor diffusion is integrated
into a foliage surface wetness factor:
ra
r
ra rs
This surface wetness factor is simply a ratio of the aerodynamic resistance to the total
resistance. When the aerodynamic resistance is small the wetness factor approaches zero
(leaf surfaces remain dry as surface moisture is readily evaporated). As the aerodynamic
resistance increases in importance relative to stomatal resistance the wetness factor
approaches 1.0 (moisture readily travels to the leaf surfaces, but is not easily evaporated).
The latent heat flux is then given by:
Here lf , is the latent heat of vaporization (J/kg), qf,sat is the saturation mixing ratio at the leaf
surface temperature, and qaf is the mixing ratio of the air within the canopy. As developed in
Frankenstein and Koenig the mixing ratio within the canopy can be determined from:
1 f 0.6 1 r " 0.11 M g
where the factor Mg (ranging from 0 to 1) is the ratio of volumetric moisture content to the
porosity of the soil (Koenig). The latent heat of vaporization (lf) is the amount of energy
required to convert a unit mass of water to vapor. It is measured in units of J/kg and is
inversely proportional to the temperature. From Henderson-Sellers it is estimated as:
2
Tf
l f 1.91846*106
T f 33.91
f g f 4 T
Fg (1 f ) I s (1 g ) g Iir g Tg4
1
Tg Tf4 H g Lg K * g
z
4/1/13 127
Advanced Surface Concepts Green Roof Model (EcoRoof)
As with the energy equation for the foliage this equation represents sensible heat flux (Hg),
latent heat flux (Lg) and the multiple reflections associated with long and short wave radiation.
The final term on the right side gives the conduction of heat into the soil substrate.
Sensible heat flux in the soil layer
Sensible heat flux between the soil surface and air in its vicinity is dependent on the
temperature difference between them and the wind speed within the canopy. It is given as
H g ag C p ,a ChgWaf (Taf Tg )
g
where C h is the bulk transfer coefficient and ag is the density of air near the soil surface
3
(kg/m ) given by:
pa pg
pag
2
The ground and foliage bulk transfer coefficients, in turn, are given by:
2
Kv
Chng rch1
Za
ln Z og
And
2
Kv
Chnf
Za Zd
ln Z of
g f
where Z o and Z o are the ground and foliage roughness lengths, rch is turbulent Schmidt
number (0.63), and Kv is the von Karman constant (0.4).
The condition of the atmosphere (h) is determined as stable or unstable based on the sign of
the bulk Richardson number:
2 gZ a Taf Tg
Rib
T af Tg Waf2
4/1/13 128
Advanced Surface Concepts Green Roof Model (EcoRoof)
The atmospheric stability factor is then given by Businger and Lumley and Panofsky as:
1.0
for Rib 0
1.0 16.0 Rib
0.5
h
1.0
for Rib 0
1.0 5.0 Rib
g
Here Ce is the bulk transfer coefficient, lg is the latent heat of vaporization at the ground
surface temperature, qaf is the mixing ratio at the foliage-atmosphere interface, and qf is the
mixing ratio at the ground surface, given by:
qg M g qg , sat 1 M g qaf
The bulk transfer coefficient for latent heat exchange is analogous to that for sensible heat
exchange and is given by:
g
where Cen is the near ground bulk transfer coefficient for Latent heat flux and e is the latent
heat exchange stability correction factor (assumed to be the same as h).
Linearization
th 4 4
In order to solve the foliage and soil heat budget equations, the 4 order terms Tf and Tg
and mixing ratio terms qg,sat and qf,sat are linearized as given by Deardorff:
4
n1 n
4 3
n n 1
T f T f 4 T f T f T f
n
4
n1 n
4 3
n n 1
Tg Tg 4 Tg Tg Tg
n
n+1 n+1
Here Tf and Tg are the current time step leaf and ground surface temperatures in Kelvin.
n n
Tf and Tg are the corresponding temperatures at the previous time step.
The saturation mixing ratio at the ground and leaf surface temperatures are given as:
q
qg , sat Tgn 1 qsat Tgn sat Tgn1 Tgn
T Tgn
4/1/13 129
Advanced Surface Concepts Green Roof Model (EcoRoof)
q
q f , sat T fn1 qsat T fn sat T fn 1 T fn
T T fn
n
where qsat(Tg ) is the saturation mixing ratio at the previous time step and is formulated as
given in Garratt:
0.622e Tgn
qsat T
n
g
P e T
g
n
Here the saturation vapor pressure e* (Pa) is evaluated at the ground temperature from the
n
previous time step (Tg ) as:
Tgn 273.15
e 611.2 exp 17.67 n
T 29.65
*
g
The derivative of saturation mixing ratio at the previous time step is given by:
dq 0.622 P de
2
dTg P 0.378 e
n
dTgn
Here, the derivative of the saturation vapor pressure can be calculated from the Clausius-
Clapeyron equation:
de
lg e Tgn
dTgn
2
Rv Tgn
Where Rv is the gas constant for water vapor and lg is the latent heat of vaporization at the
soil surface temperature.
The corresponding saturation mixing ratio relations for the leaf surfaces can be obtained by
replacing Tg with Tf in the above relations.
Final Equations
After linearization the final equations are of the form:
The coefficients in these equations result from the direct combination of the equations from
the above development. The interested reader is directed to the papers by Frankenstein and
Koenig for the complete and somewhat complicated expressions.
This final set of equations is then solved simultaneously to obtain Tg and Tf . One key
difference in our implementation of the FASST algorithm is that the conduction terms in the
4/1/13 130
Advanced Surface Concepts Green Roof Model (EcoRoof)
g g
equations for C1 and C2 are solved by inverting the Conduction Transfer Functions (CTF)
within the EnergyPlus solution scheme.
Green Roof Nomenclature
C1, C2, C3 = coefficients in linearized temperature equations
g
Ce = latent heat flux bulk transfer coefficient at ground layer
Cf = bulk heat transfer coefficient
g
Ch = sensible heat flux bulk transfer coefficient at ground layer
f
Chn = near-neutral transfer coefficient at foliage layer
g
Chn = near-neutral transfer coefficient at ground layer
Cp,a = specific heat of air at constant pressure (1005.6 J/kg k)
*
e = saturation vapor pressure (Pa)
f1 = multiplying factor for radiation effect on stomatal resistance
f2 = multiplying factor for moisture effect on stomatal resistance
f3 = additional multiplying factor for stomatal resistance
2
Ff = net heat flux to foliage layer (W/m )
2
Fg = net heat flux to ground surface (W/m )
gd = plant specific characteristic related to stomatal resistance
2
Hf = foliage sensible heat flux (W/m )
2
Hg = ground sensible heat flux (W/m )
I s = total incoming short wave radiation (W/m )
2
4/1/13 131
Advanced Surface Concepts Green Roof Model (EcoRoof)
4/1/13 132
Advanced Surface Concepts Green Roof Model (EcoRoof)
Dickinson, R.E., A. Henderson-Sellers, P.J. Kennedy, and M.F. Wilson. 1986. Biosphere-
Atmosphere Transfer Scheme (BATS) for the NCAR community climate model. NCAR
Technical Note, TN-275+STR.
ECMWF. 2002. European Centre for Medium-Range Weather Forecasts, Integrated Forecast
System. Documentation, CY25R1 (Operational implementation 9 April 2002).
http://www.ecmwf.int/research/ifsdocs/CY25r1/Physics/Physics-08-03.html.
Frankenstein, S., and G. Koenig. 2004. FASST Vegetation Models. U. S. Army Engineer
Research and Development Center, Cold regions Research and Engineering Laboratory,
ERDC/CRREL Technical Report TR-04-25.
Frankenstein, S., and G. Koenig. 2004. Fast All-season Soil Strength (FASST). U.S. Army
Engineer Research and Development Center, Cold regions Research and Engineering
Laboratory, ERDC/CRREL Special Report SR-04-1.
Garratt, J.R. 1992. The Atmospheric Boundary Layer, Cambridge university press.
Gates, D.M. 1980. Biophysical Ecology. New York: Springer-Verlag
Guymon, G.L., R.L. Berg, and T.V. Hromadka. 1993. Mathematical Model of Frost Heave and
Thaw Settlement in Pavements. U.S. Army Cold Regions Research and Engineering
Laboratory, CRREL Report 93-2.
Henderson-Sellers, B. 1984. “A New Formula for Latent Heat of Vaporization of water as
function of temperature”, Quarterly Journal Royal Meteorological Society, 10 pp. 1186-1190.
Hughes, P.A., T.J.L. McComb, A.B. Rimmer, and K.E. Turver. 1993. “A mathematical model
for the prediction of temperature of man-made and natural surfaces”, International Journal of
Remote Sensing 14 (7), pp. 1383-1412.
Koenig, G.G. 1994. Smart Weapons Operability Enhancement (SWOE) Joint Test and
Evaluation (JT and E) Program: Final Report. Dr. James P. Welch, Joint Test Director,
SWOE JT and E, SWOE Report 94-10, Annex D.
Lumley, J. L. and Panofsky, H. A. 1964. ‘The structure of Atmospheric Turbulence’.
Interscience Monographs and Texts in Physics and Astronomy, Vol. XII. Wiley, New York.
Oke, T.R. 1987. Boundary Layer Climates, University Press, Cambridge
Sellers, P.J., Y. Mintz, Y.C. Sud, and A. Dalcher. 1986. A simple biosphere model (SiB) for
use within general circulation models. Journal of Atmospheric Science, 43 (6), pp. 505-532.
4/1/13 133
Climate, Sky and Solar/Shading Calculations Climate Calculations
Climate Calculations
The location of the facility under analysis is critical for the determination of energy
consumption, heating/cooling loads, daylighting potential, and a host of other calculations. In
EnergyPlus, both external (i.e, weather files supplied from others) and internal (i.e., solar
position, design day temperature/humidity/solar profiles) data is used during simulations.
The “Site:Location” input object includes parameters (Latitude, Longitude, Elevation,
Timezone) that allow EnergyPlus to calculate the solar position (using Latitude, Longitude
and Timezone) for any day of the year as well as supply the standard barometric pressure
(using elevation). Solar position modeling is discussed in more detail in both the Sky
Radiance and Shading Calculation sections that directly follow this section.
Weather files have hourly or sub-hourly data for each of the critical elements needed during
the calculations (i.e., Dry-Bulb Temperature, Dew-Point Temperature, Relative Humidity,
Barometric Pressure, Direct Normal Radiation, Diffuse Horizontal Radiation, Total & Opaque
Sky Cover, Wind Direction, Wind Speed) as well as some auxiliary data such as Rain or
Snow that assist in certain calculational aspects. Weather file excerpts such as might be used
in sizing calculations also have this breadth of data. The input object
“SizingPeriod:DesignDay” describes design days (meant to mimic ASHRAE design
conditions but in a whole day profile) using certain characteristics for the day and then
EnergyPlus supplies the remaining portions to complete outdoor conditions needed for
program execution. SizingPeriod:DesignDay are perhaps the best objects for sizing
equipment as the ASHRAE specified design conditions can be input AND weather files may
or may not have the conditions necessary to size equipment properly.
Two other objects, however, can be used at times: SizingPeriod:WeatherFileDays and
SizingPeriod:WeatherFileConditionType. With the first of these, one specifies a set of
weather file days similar to a RunPeriod (but will be used for the Sizing calculations). In the
second, the extreme or typical conditions that are calculated for the weather file during
processing can be used by name. Of course, either of these can be used as a measure of
usage over small periods.
The ASHRAE Handbook of Fundamentals describes their criteria in creating design condition
synopses. “Design data based on dry-bulb temperature represent peak occurrences of the
sensible component of ambient outdoor conditions. Design values based on wet-bulb
temperature are related to the enthalpy of the outdoor air. Conditions based on dewpoint
relate to the peaks of the humidity ratio. The designer, engineer, or other user must decide
which set(s) of conditions and probability of occurrence apply to the design situation under
consideration.”
EnergyPlus Design Day Temperature Calculations
In EnergyPlus, the typical design day input includes a “high” and a “low” dry-bulb temperature
for the day. As these two temperatures are insufficient to represent a full 24 hour period, the
program uses a “range multiplier” profile to represent the full day’s temperatures:
4/1/13 134
Climate, Sky and Solar/Shading Calculations Climate Calculations
1.2
1
Multiplier Value
0.8
0.6
0.4
0.2
0
AM
AM
AM
AM
AM
AM
AM
AM
AM
AM
AM
AM
PM
PM
PM
PM
PM
PM
PM
PM
PM
PM
PM
PM
0
0
00
00
00
00
00
00
00
00
00
00
00
00
00
00
00
00
00
00
:0
:0
:0
:0
:0
:0
1:
2:
3:
4:
5:
6:
7:
8:
9:
1:
2:
3:
4:
5:
6:
7:
8:
9:
10
11
12
10
11
12
Figure 35. Default Daily Temperature Range Profile
The multipliers are taken from the ASHRAE 2009 HOF, Table 6, p. 14.11.. More explicitly,
EnergyPlus creates an air temperature for each timestep by using the entered maximum dry-
bulb temperature in conjunction with the entered daily range and the above multiplier values.
The actual equation used is shown below:
where
Tcurrent= Air temperature of current Hour of Day
TMax= User supplied Max Dry-bulb Temperature
Trange= User supplied Daily Temperature Range
TMultiplier= Range multiplier as shown on the above graph
The range multiplier values represent typical conditions of diurnal temperatures (i.e. the low
temperature for the day occurring about 5:00 AM and the maximum temperature for the day
occurring about 3:00 PM. Note that EnergyPlus does not shift the profile based on the time
of solar noon as is optionally allowed by ASHRAE procedures.
ASHRAE research indicates that dry-bulb and wet-bulb temperatures typically follow the
same profile, so EnergyPlus can use the default daily temperature profile to generate
humidity conditions based on maximum and range of wet-bulb temperature.
Since this default temperature profile may not be applicable to all locations, the user can give
a different profile as part of the design day definition.
Sky Radiation Modeling
2,
EnergyPlus calculates the Horizontal Infrared Radiation Intensity in Wh/m if it is missing on
the weather file or for design days, from the Opaque Sky Cover field as shown in the
following explanation.
where
2
Horizontal_IR = horizontal IR intensity {W/m }
Skyemissivity = sky emissivity
2 4
Sigma = Stefan-Boltzmann constant = 5.6697e-8 {W/m -K }
4/1/13 135
Climate, Sky and Solar/Shading Calculations Climate Calculations
Temperaturedewpoint
Skyemissivity .787 .764ln 1. .0224 N .0035 N .00028 N
2 3
273.
where
Temperaturedewpoint = dewpoint temperature {K}
N = opaque sky cover {tenths}
Example: Clear sky (N=0), Temperaturedrybulb = 273+20=293 K, Temperaturedewpoint =
273+10=283 K:
Skyemissivity = 0.787 + 0.764*0.036 = 0.815
2
Horizontal_IR = 0.815*5.6697e-8*(293**4) = 340.6 W/m
References for these calculations are contained in the references section at the end of this
list of fields. (Walton, 1983) (Clark and Allen, 1978).
EnergyPlus Sky Temperature Calculation
The default calculation for sky temperature is:
.25
Horizontal _ IR
SkyTemperature TemperatureKelvin
Sigma
Where
SkyTemperature = Sky radiative temperature {C}
Horiizontal_IR = Horizontal Infrared Radiation Intensity as described in the previous section
2
{Wh/m }
TemperatureKelvin = Temperature conversion from Kelvin to C, i.e. 273.15
The Sky Temperature can also be set by the user from several options using the
WeatherProperty:SkyTemperature object.
EnergyPlus Design Day Solar Radiation Calculations
Similarly, calculating solar irradiance (solar modeling) is one of the important effects to be
accomplished. Several solar models exist with varying complexity.
ASHRAE Clear Sky Solar Model
The default model used is the ASHRAE Clear Sky model. The ASHRAE clear sky model
described in ASHRAE HOF 2005 Chapter 31 can be used to estimate hourly clear-day solar
radiation for any month of the year in U.S. or similar temperate climates in the northern
hemisphere. EnergyPlus calculations extend the clear sky application to both northern and
southern hemispheres. Note that the Clear Sky model has been updated in the ASHRAE
HOF 2009 (see ASHRAE Revised Clear Sky Model below).
At the earth’s surface on a clear day, direct normal irradiation is represented by
A
Direct Normal Irradiation
(128)
exp B
sin
where
A = apparent solar irradiation at air mass m = 0 (Table 17)
4/1/13 136
Climate, Sky and Solar/Shading Calculations Climate Calculations
Io Equation Declination A B C
{W/m2} of Time {degrees} {W/m2} {} {}
{minutes}
Jan 1416 -11.2 -20.0 1202 0.141 0.103
Feb 1401 -13.9 -10.8 1187 0.142 0.104
Mar 1381 -7.5 0.0 1164 0.149 0.109
Apr 1356 1.1 11.6 1130 0.164 0.120
May 1336 3.3 20.0 1106 0.177 0.130
Jun 1336 -1.4 23.45 1092 0.185 0.137
Jul 1336 -6.2 20.6 1093 0.186 0.138
Aug 1338 -2.4 12.3 1107 0.182 0.134
Sep 1359 7.5 0.0 1136 0.165 0.121
Oct 1380 15.4 -10.5 1166 0.152 0.111
Nov 1405 13.8 -19.8 1190 0.144 0.106
Dec 1417 1.6 -23.45 1204 0.141 0.103
For locations where clear, dry skies predominate (e.g., at high elevations) or, conversely,
where hazy and humid conditions are frequent, values found by using Equation (128) and
Table 17 should be multiplied by the clearness numbers in Threlkeld and Jordan (1958),
reproduced as Figure 5 in Chapter 33 of the 2007 ASHRAE Handbook—HVAC Applications.
The Clear Sky model usually over estimates the amount of solar radiation available to the
building.
ASHRAE Revised Clear Sky Model (“Tau Model”)
The ASHRAE 2009 HOF introduced a revised clear sky model based on location-specific
optical depths for direct and diffuse radiation. These values are tabulated by month for all
5564 locations in the ASHRAE design data that accompanies the 2009 HOF.
The model requires air mass, m, calculated as follows --
Eb Eo exp b m ab
4/1/13 137
Climate, Sky and Solar/Shading Calculations Climate Calculations
Ed Eo exp d m ad
where:
Eb = beam normal irradiance, W/m2
Ed = diffuse horizontal irradiance,
Eo = extraterrestrial normal irradiance,
m = air mass
b and d = beam and diffuse optical depths (from ASHRAE climatic design data)
ab and ad = beam and diffuse air mass exponents (see below)
Values of b and d are location-specific and vary during the year. They embody the
dependence of clear sky solar radiation upon local conditions, such as elevation, precipitable
water content, and aerosols.
The air mass exponents ab and ad were correlated to b and d through the following
empirical relationships:
Studies done as part of ASHRAE research projects show that the revised tau model
produces more physically plausible irradiance values than does the traditional clear sky
model. In particular, diffuse irradiance values are more realistic.
Zhang-Huang Solar Model
The Zhang-Huang solar model was developed for initial use in modeling typical
meteorological years for China weather sites. This model seems to be good for other
locations as well. Using total cloud cover, dry-bulb temperature, relative humidity, and wind
speed as the independent variables, the total (global horizontal) solar radiation is estimated
by:
I 0 sin(h) c0 c1 CC c2 CC 2 c3 Tn Tn 3 c4 c5 Vw d
I
k
Where
2
I = estimated hourly solar radiation, W/m
2
I0= global solar constant, 1355 W/m
h=solar altitude angle, i.e, the angle between the horizontal and the line to the sun
CC = cloud cover
= relative humidity, %
Tn , Tn 3 = dry-bulb temperature at hours n (current) and n-3, respectively
Vw = Wind speed, m/s
c0 , c1 , c2 , c3 , c4 , c5 , d , k = regression coefficients
The constants were determined by analysis from measured data and are as follows:
4/1/13 138
Climate, Sky and Solar/Shading Calculations Design Week Specification
A design week (or longer) may be specified in EnergyPlus by selecting a RunPeriod of one
week. One-week custom weather files may be constructed if extreme design week
conditions are desired or one can use the specific dates from the EPW weather files that
designate the extreme conditions on the weather file. The design week simulation requires
that the simulation reach a ‘steady-periodic’ state for the one week time period. This is
achieved manually by increasing the ‘number of years of simulation’ until the weekly profile
does not change from one year to the next. The ‘steady-periodic’ output for the last week of
the simulation represents the design week simulation results.
Using the SizingPeriod:WeatherFileDays object, one could size equipment from a longer
period.
4/1/13 139
Climate, Sky and Solar/Shading Calculations Sky Radiance Model
In EnergyPlus the calculation of diffuse solar radiation from the sky incident on an exterior
surface takes into account the anisotropic radiance distribution of the sky. For this
distribution, the diffuse sky irradiance on a surface is given by
AnisoSkyMultipliersurface DiffuseSolarIrradiance
Where
Diffuse Solar Irradiance is the diffuse solar irradiance from the sky on the ground.
surface is the surface being analyzed.
AnisoSkyMultiplier is determined by surface orientation and sky radiance distribution, and
accounts for the effects of shading of sky diffuse radiation by shadowing surfaces such as
overhangs. It does not account for reflection of sky diffuse radiation from shadowing surfaces.
The sky radiance distribution is based on an empirical model based on radiance
measurements of real skies, as described in Perez et al., 1990. In this model the radiance of
the sky is determined by three distributions that are superimposed (see Figure 36)
(1) An isotropic distribution that covers the entire sky dome;
(2) A circumsolar brightening centered at the position of the sun;
(3) A horizon brightening.
Circumsolar brightening
Isotropic dome (concentrated at center of sun)
Horizon brightening
(concentrated at horizon)
Figure 36. Schematic view of sky showing solar radiance distribution as a superposition of three
components: dome with isotropic radiance, circumsolar brightening represented as a point source at the
sun, and horizon brightening represented as a line source at the horizon.
The proportions of these distributions depend on the sky condition, which is characterized by
two quantities, clearness factor and brightness factor, defined below, which are determined
from sun position and solar quantities from the weather file.
The circumsolar brightening is assumed to be concentrated at a point source at the center of
the sun although this region actually begins at the periphery of the solar disk and falls off in
intensity with increasing angular distance from the periphery.
The horizon brightening is assumed to be a linear source at the horizon and to be
independent of azimuth. In actuality, for clear skies, the horizon brightening is highest at the
horizon and decreases in intensity away from the horizon. For overcast skies the horizon
brightening has a negative value since for such skies the sky radiance increases rather than
decreases away from the horizon.
4/1/13 140
Climate, Sky and Solar/Shading Calculations Sky Radiance Model
Table 18. Variables in Anisotropic Sky Model and Shadowing of Sky Diffuse Radiation
4/1/13 141
Climate, Sky and Solar/Shading Calculations Sky Radiance Model
where
Ihm / Io
where
m = relative optical air mass
2
Io = extraterrestrial irradiance (taken to have an average annual value of 1353 W/m );
and the sky clearness factor is
4/1/13 142
Climate, Sky and Solar/Shading Calculations Sky Radiance Model
(Ih I ) / Ih Z 3
1 Z 3
where
I = direct normal solar irradiance
κ = 1.041 for Z in radians
The factors Fij are shown in the following table. The Fij values in this table were provided by
R. Perez, private communication, 5/21/99. These values have higher precision than those
listed in Table 6 of Perez et al., 1990.
Table 19. Fij Factors as a Function of Sky Clearness Range.
ε Range 1.000-1.065 1.065-1.230 1.230-1.500 1.500-1.950 1.950-2.800 2.800-4.500 4.500-6.200 > 6.200
F11 -0.0083117 0.1299457 0.3296958 0.5682053 0.8730280 1.1326077 1.0601591 0.6777470
F12 0.5877285 0.6825954 0.4868735 0.1874525 -0.3920403 -1.2367284 -1.5999137 -0.3272588
F13 -0.0620636 -0.1513752 -0.2210958 -0.2951290 -0.3616149 -0.4118494 -0.3589221 -0.2504286
F21 -0.0596012 -0.0189325 0.0554140 0.1088631 0.2255647 0.2877813 0.2642124 0.1561313
F22 0.0721249 0.0659650 -0.0639588 -0.1519229 -0.4620442 -0.8230357 -1.1272340 -1.3765031
F23 -0.0220216 -0.0288748 -0.0260542 -0.0139754 0.0012448 0.0558651 0.1310694 0.2506212
24
I
Irradiance from horizon without obstructions
i
i 1
th
where Ii is the unobstructed irradiance on the surface from the i interval, SFi is the sunlit
th
fraction from radiation coming from the i interval, and the sums are over intervals whose
center lies in front of the surface. SFi is calculated using the beam solar shadowing method
th
as though the sun were located at the i horizon point. Here
where
E (θi) = radiance of horizon band (independent of θ)
dθ = 2π/24 = azimuthal extent of horizon interval (radians)
O O O
θi = 0 , 15 , … , 345
αi = incidence angle on surface of radiation from θi
The corresponding ratio for the isotropic sky dome is given by
4/1/13 143
Climate, Sky and Solar/Shading Calculations Sky Radiance Model
24 6
I
Irradiance from dome without obstructions
ij
i 1 j 1
where (i,j) is a grid of 144 points (6 in altitude by 24 in azimuth) covering the sky dome, Iij is
th
the unobstructed irradiance on the surface from the sky element at the ij point, SFij is the
th
sunlit fraction for radiation coming from the ij element, and the sum is over points lying in
front of the surface. Here
I ij E ( i , j ) cos j d d cos ij
where
E (θi,φj) = sky radiance (independent of θ and φ for isotropic dome)
dθ = 2π/24 = azimuthal extent of sky element (radians)
dφ = (π/2)/6 = altitude extent of sky element (radians)
O O O
θi = 0 , 15 , … , 345
O O O
φj = 7.5 , 22.5 , … , 82.5
αij = incidence angle on surface of radiation from (θi,φj)
Because the circumsolar region is assumed to be concentrated at the solar disk, the
circumsolar ratio is
where SFsun is the beam sunlit fraction. The total sky diffuse irradiance on the surface with
shadowing is then
Rhorizon and Rdome are calculated once for each surface since they are independent of sun
position.
With shadowing we then have:
AnisoSkyMult = I’sky /DifSolarRad.
Shadowing of Sky Long-Wave Radiation
EnergyPlus calculates the sky long-wave radiation incident on exterior surfaces assuming
that the sky long-wave radiance distribution is isotropic. If obstructions such as overhangs are
present the sky long-wave incident on a surface is multiplied by the isotropic shading factor,
Rdome, described above. The long-wave radiation from these obstructions is added to the
long-wave radiation from the ground; in this calculation both obstructions and ground are
assumed to be at the outside air temperature and to have an emissivity of 0.9.
4/1/13 144
Climate, Sky and Solar/Shading Calculations Shading Module
Shading Module
Solar Position
Current solar position is described in terms of three direction cosines that are convenient for
determining the angle of incidence of the sun’s rays on a building surface. The following
procedure is used to determine the direction cosines. The values of the solar declination
4/1/13 145
Climate, Sky and Solar/Shading Calculations Shading Module
angle, , and the equation of time, , are based on Astronomical Algorithms, Meeus. Solar
declination is a function of local/site latitude.
The fractional year is calculated, in radians:
2
( day _ of _ year )
366
From this fractional year, the equation of time and solar declination angle are calculated. For
each time step (time value = fractional hour), the hour angle is calculated from:
HourAngle 1512 TimeValue EquationOfTime TimeZoneMeridian Longitude
TimeZoneMeridian is the standard meridian for the location’s time zone {GMT +/-}.
Solar HourAngle (H) gives the apparent solar time for the current time period (degrees);
HourAngle is positive before noon, negative after noon. It is common astronomical practice
to express the hour angle in hours, minutes and seconds of time rather than in degrees. You
can convert the hour angle displayed from EnergyPlus to time by dividing by 15. (Note that 1
hour is equivalent to 15 degrees; 360 of the Earth’s rotation takes place every 24 hours.)
The relationship of angles in degrees to time is shown in the following table:
Table 20. Relationship of Angles (degrees) to Time
4/1/13 146
Climate, Sky and Solar/Shading Calculations Shading Module
Surface Geometry
Shadow calculations first require that the building surfaces be described geometrically.
Surfaces are described by the coordinates of their vertices in a three dimensional Cartesian
coordinate system. This Right-hand coordinate system has the X-axis pointing east, the Y-
axis pointing north, and the Z-axis pointing up (see figure below). The azimuth angle () of a
surface is the angle from the north axis to the projection onto the X-Y plane of a normal to the
surface (clockwise positive). The surface tilt angle () is the angle between the Z-axis and the
normal to the surface. The vertices are recorded in counter-clockwise sequence (as the
surface is viewed from outside its zone).
During surface entry, surfaces are checked for convex or non-convex shape. If non-convex
and inappropriate (used as a receiving surface) then a severe error is produced telling the
user that shadowing calculations may be inaccurate.
Similarly collinear points (or as noted below, points within 1 mm distance) are removed
unless removing would make an illegal surface (less than 3 points). But degenerate collinear
surfaces should be removed – they make the shadowing routines do extra work which takes
extra time.
Collinear – points that essentially form a “line” rather than a surface shape.
Resolution of 1mm or less – near collinear points.
Note that the resolution on surfaces/shadowing is 1 mm – using resolution beyond that will result in
truncation of the shadowing.
4/1/13 147
Climate, Sky and Solar/Shading Calculations Shading Module
The GlobalGeometryRules object specifies to EnergyPlus how the surface vertices will be
presented in the input file. Of pertinent interest here is that the user may specify the vertices
in either “relative” or “world” coordinates. Regardless of input specifications, when vertices
are reported, they are reported in world coordinates, starting at the upper-left-corner (4-sided
surface) and are listed counter-clockwise.
Relative Coordinate Transformation
When vertices are specified in “relative” coordinates, there can be a “building” north axis as
well as a “zone” north axis. The building north axis/coordinate system is a rotation of b
degrees from the global/world coordinate system. The global coordinates of zone origins are
related to the building relative coordinates by:
Z zo Z br (131)
Where
zo – represents Zone Origin
br – represents the Zone Origin as input (relative to building origin)
The zone may also be rotated z degrees relative to the building coordinates. Origins of zone
surfaces are then given relative to the zone coordinate system. The global coordinates of the
surface origins are calculated by:
4/1/13 148
Climate, Sky and Solar/Shading Calculations Shading Module
A surface azimuth angle relative to the zone coordinate system (s) is converted to a global
azimuth by:
s z b (135)
The surface tilt angle () is not changed by these rotations about the Z-axis.
The coordinates of the surface vertices are given in a coordinate system in the plane of the
surface relative to the second vertex as shown for surfaces in Figure 39. The X-axis of the
surface coordinate system is a horizontal line through the second vertex. The global
coordinates of the surface vertices are given by:
X ' X X so (139)
Z ' Z Z so (141)
Ysr X ' sin cos Y ' cos cos Z ' sin (143)
Shadow Projection
All architectural forms are represented by plane polygons. This can give good accuracy even
for curved surfaces: a sphere can be approximated by the 20 nodes of an icosahedron with
only 3 percent error in the shadow area cast by the sphere. Consider how a solid object,
which is composed of a set of enclosing plane polygons, casts a shadow. Figure 40 shows a
box shaped structure on a horizontal surface. The structure consists of a top (surface 1) and
four vertical surfaces (2 and 3 visible to the observer and 4 and 5 not visible). The sun is
positioned behind and to the right of the structure and a shadow is cast onto the horizontal
surface (the ground).
Surfaces 1, 4, and 5 are in sunlight; 2 and 3 are in shade. It is possible to think of the
structure's shadow as the combination of shadows cast by surfaces 1, 2, 3, 4 and 5 or by 1, 4
and 5, or by surfaces 2 and 3. This last combination of shadow casting surfaces is the
simplest. In the EnergyPlus shadow algorithm every surface is considered to be one of the
4/1/13 149
Climate, Sky and Solar/Shading Calculations Shading Module
surfaces that enclose a solid, and only those surfaces that are not sunlit at a given hour are
considered shadowing surfaces.
The expressions in equation (144) are the direction cosines of the surface:
The cosine of the angle of incidence of the sun's rays on the surface are given by the dot
product of surface and sun direction cosines.
4/1/13 150
Climate, Sky and Solar/Shading Calculations Shading Module
This is done by finding, through linear interpolation, the points on the perimeter of the SP,
which intersect the plane of the RP. These points become new vertices of the SP, which
together with the other positive vertices define a clipped SP that casts only a real shadow.
A vertex located at (x, y, z) relative to the RP coordinate system casts a shadow to a point in
the plane of the RP given by
z a
x' x (149)
cos
z b
y' y (150)
cos
where
and
More explicitly, a casting surface – a shadow casting surface or general casting surface – is
one that casts a shadow on other surfaces. A receiving surface – a shadow receiving surface
– is one that receives shadows from other surfaces (i.e. casting surfaces). A back surface –
an inside surface – is one that may be partially sunlit/receive solar transmission for interior
solar distribution.
4/1/13 151
Climate, Sky and Solar/Shading Calculations Shading Module
Homogeneous Coordinates
Two-dimensional homogeneous coordinate techniques are used to determine the vertices of
shadow overlaps. In homogeneous coordinates, points and lines are represented by a single
form that allows simple vector operations between those forms [Newman-Sproul]. A point (X,
Y) is represented by a three element vector (x, y, w) where x = w*X, y = w*Y, and w is any
real number except zero. A line is also represented by a three element vector (a, b, c). The
directed line (a, b, c) from point (x1, y1, w1) to point (x2, y2, w2) is given by:
(a, b, c) ( x1 , y1 , z1 ) ( x2 , y2 , z2 ) (151)
The sequence in the cross product is a convention to determine sign. The condition that a
point (x, y, w) lie on a line (a, b, c) is that
( a , b , c ) ( x , y , w) 0 (152)
( a, b, c ) ( x / w, y / w,1) 0 (153)
the point is to the left of the line. If it is less than zero, the point is to the right of the line. The
intercept (x, y, w) of line (a1, b1, c1) and line (a2, b2, c2) is given by:
Note that the use of homogeneous coordinates as outlined above provides a consistent
method and notation for defining points and lines, for determining intercepts, and for
determining whether a point lies to the left, to the right, or on a line. Normalization provides
the means for transforming to and from homogeneous notation and Cartesian coordinates.
Thus, if (X, Y) is a Cartesian coordinate pair, its homogeneous coordinates are (X, Y, 1).
Similarly, the homogeneous coordinates (x, y, w) can be transformed to the Cartesian point
with coordinates (x/w, y/w).
Polygon Clipping Algorithms
Two methods for polygon clipping (treating of overlapping shadows) are currently in use in
EnergyPlus.
Convex Weiler - Atherton
Sutherland – Hodgman
The original EnergyPlus method for polygon clipping is a special version of the Weiler-
Atherton model (Weiler, Atherton, 1977). It was developed to be sufficiently general to clip
concave polygons with holes. The implementation in the current version of EnergyPlus,
however, does not support concave shadowing surfaces or holes. The relative computational
complexity is preserved – the algorithm is carried out in four steps. For example, if A and B
are polygons (see Figure 42).
1) A call to INCLOS determines which vertices of X lie within Y.
2) A second call determines which vertices of Y lie within X.
3) If neither polygon is contained completely within the other, INTCPT is called to collect
points of intersection between X and Y.
4) Since the points are usually gathered out of order, they must then be oriented.
The Sutherland-Hodgman algorithm (Sutherland, Hodgman, 1974) is less complex compared
to the Weiler-Atherton method and is well-suited to clipping convex polygons. In actuality,
4/1/13 152
Climate, Sky and Solar/Shading Calculations Shading Module
only convex shading surfaces are currently supported by EnergyPlus. Let X be a polygon
called the “subject polygon” (SP) and Y be a polygon called the “clipping polygon” (CP). The
method performs the computation by iterating over the edges of the CP and removing points
from the SP that fall in the clipping plane, i.e. points that fall to the left of the edge of the CP.
Intersections between the clip edge and the edges of the SP are added appropriately, and
points falling outside of the clipping plane, i.e. to the right of the edge of the CP, are added
the output polygon as well. This resultant polygon is stored and the process is repeated for
the rest of the clip edges in CP. The process is analogous to cutting off pieces of the SP one-
by-one with respect to each edge of the CP. The result is ordered and identical to the polygon
produced by the Weiler-Atherton method.
Overlapping Shadows
After transforming the shadows onto the plane of the receiving surface, the basic job of the
shadow algorithm is to determine the area of the overlap between the polygons representing
the shadows and the polygon representing the receiving surface. Concave surfaces are
supported only for exterior wall heat transfer surfaces, when using SutherlandHodgman
option. Concave shading devices are not supported by the this option. Neither concave
shading devices nor concave exterior wall heat transfer surfaces are supported by the
ConvexWeilerAtherton clipping routine.
When only convex shading devices are considered, this provides a great simplification. The
overlap between two convex polygons (i.e. projections of shading devices via the direction of
the sun) is another convex polygon. Coordinate and projection transformations of a convex
polygon produce another convex polygon. Any non-convex polygon can be constructed as
the union of convex ones.
For ConvexWeilerAtherton, there is considerable simplification if only convex (no interior
angle > 180 ) polygons are considered. The overlap between two convex polygons is another
convex polygon. Coordinate and projection transformations of a convex polygon produce
another convex polygon. Any non-convex polygon can be constructed as a sum of convex
ones.
The vertices that define the overlap between two convex polygons, A and B, consist of:
the vertices of A enclosed by B
the vertices of B enclosed by A
and the intercepts of the sides of A with the sides of B
In Figure 42, point a is the result of rule 1, point c is the result of rule 2, and points b and d
result from rule 3. The overlap of A and B is the polygon a-b-c-d. Figure 43 shows an overlap
where all of the vertices of B are enclosed by A. Figure 44 shows an overlap defined only by
the intercepts of A and B. Figure 45 shows a more complex overlap.
4/1/13 153
Climate, Sky and Solar/Shading Calculations Shading Module
Coordinate transformation retains the order of the vertices of a polygon, while a projection
reverses the order. The sequence of vertices of the receiving polygons should be reversed so
it and all shadow polygons will have the same sequence.
A point is enclosed by a clockwise, convex polygon if the point lies to the right of all sides (or
does not lie to the left of any side) of the polygon. The intercept of two sides may not lie
beyond the ends of either side. These are "line segments" rather than "lines". It is possible to
tell if line segments A and B intercept within their end points by noting that the ends of A must
lie on both sides of B, and the ends of B must lie on both sides of A. This should be done
before the intercept is calculated.
4/1/13 154
Climate, Sky and Solar/Shading Calculations Shading Module
Once the vertices are determined, they must be sorted into clockwise order for the area to be
computed. Given a closed, planar polygon of n sequential vertices (x1, y1), (x2, y2) …, (xn, yn),
its area is given:
n
Area 1
2 (x y
i 1
i i 1 xi 1 yi ) (155)
If two shadows overlap the receiving surface, they may also overlap each other as in Figure
46. The vertices of this overlap can be computed. The areas of all overlaps can be
computed. The total sunlit area can be expressed as the sum of all polygon areas given a
proper sign on each of the areas.
4/1/13 155
Climate, Sky and Solar/Shading Calculations Shading Module
Partially transparent shadowing surfaces can also be modeled by giving a transparency () to
every shadowing polygon. Let of the receiving polygon be one. Then the of every overlap
of polygons i and j is the product of i and j The shaded area is then computed by summing
Ai(1 - i) for all overlap polygons.
It is easy to determine the sunlit area of a window once all the shadow and overlap vertices
on the wall have been computed. Consider wall 2 of Figure 37. First, the wall is considered a
simple rectangle and the window on it is ignored. The shadow overlapping is performed and
the sunlit portion of the gross wall area is computed. Then the window rectangle is
overlapped with the shadow to determine its sunlit area. The sunlit area of the window is
subtracted from the gross wall sunlit area to determine the net wall sunlit area. During this
calculation it is not necessary to recompute the shadows, because they were precisely
determined on the wall.
When the SutherlandHodgman option is selected, the overlap is computed using the
Sutherland-Hodgman algorithm for polygon clipping when. Let X be a polygon called the
“subject polygon” (SP) and Y be a polygon called the “clipping polygon” (CP). The method
performs the computation by iterating over the edges of the CP and removing points from the
SP that fall in the clipping plane, i.e. points that fall to the left of the edge of the CP. If it is to
the left of any edge, it the point does not overlap with the CP. Intersections between the clip
edge and the edges of the SP are added appropriately, and points falling outside of the
4/1/13 156
Climate, Sky and Solar/Shading Calculations Shading Module
clipping plane, i.e. to the right of the edge of the CP, are added the output polygon as well.
This resultant polygon is stored and the process is repeated for the rest of the clip edges in
CP. The process is analogous to cutting off pieces of the SP one-by-one with respect to each
edge of the CP. Note that the SP may be concave, but the CP may not. This means that the
exterior wall surfaces may be concave, while shading devices may not be concave.
Solar Gains
The total solar gain on any exterior surface is a combination of the absorption of direct and
diffuse solar radiation given by
S
Qso I b cos s I s Fss I g Fsg (156)
S
where
a =solar absorptance of the surface
A =angle of incidence of the sun's rays
S =area of the surface
Ss = sunlit area
Ib =intensity of beam (direct) radiation
Is =intensity of sky diffuse radiation
Ig =intensity of ground reflected diffuse radiation
Fss = angle factor between the surface and the sky
Fsg = angle factor between the surface and the ground
For the surface of a building located on a featureless plain
1 cos
Fss (157)
2
and
1 cos
Fsg (158)
2
If the surface is shaded the program modifies Fss by a correction factor that takes into
account the radiance distribution of the sky (see “Shadowing of Sky Diffuse Solar Radiation”).
Shading of ground diffuse solar radiation is not calculated by the program. It is up to the user
to estimate the effect of this shading and modify the input value of Fsg accordingly.
Solar Distribution
As discussed in the Input Output Reference (Object: Building), the field Solar Distribution, in
the “Building” input object, determines how EnergyPlus will treat beam solar radiation
entering a zone through exterior windows. There are five choices: MinimalShadowing,
FullExterior, FullInteriorAndExterior, FullExteriorWithReflections, and
FullInteriorAndExteriorWithReflections.
MinimalShadowing
In this case, there is no exterior shadowing except from window and door reveals. All beam
solar radiation entering the zone is assumed to fall on the floor, where it is absorbed
according to the floor's solar absorptance. Any reflected by the floor is added to the
transmitted diffuse radiation, which is assumed to be uniformly distributed on all interior
4/1/13 157
Climate, Sky and Solar/Shading Calculations Shading Module
surfaces. If no floor is present in the zone, the incident beam solar radiation is absorbed on
all interior surfaces according to their absorptances. The zone heat balance is then applied at
each surface and on the zone's air with the absorbed radiation being treated as a flux on the
surface.
FullExterior
In this case, shadow patterns on exterior surfaces caused by detached shading, wings,
overhangs, and exterior surfaces of all zones are computed. As for MinimalShadowing,
shadowing by window and door reveals is also calculated. Beam solar radiation entering the
zone is treated as for MinimalShadowing.
FullExteriorWithReflections
This case is the same interior distribution as the preceding option but uses exterior reflections
as well (see the section Solar Radiation Reflected from Exterior Surfaces for further
explanation).
FullInteriorAndExterior
This is the same as FullExterior except that instead of assuming all transmitted beam solar
falls on the floor the program calculates the amount of beam radiation falling on each surface
in the zone, including floor, walls and windows, by projecting the sun's rays through the
exterior windows, taking into account the effect of exterior shadowing surfaces and window
shading devices.
If this option is used, you should be sure that the surfaces of the zone totally enclose a
space. This can be determined by viewing the eplusout.dxf file with a program like
AutoDesk’s Volo View Express. You should also be sure that the zone is convex. Examples
of convex and non-convex zones are shown in Figure 47. The most common non-convex
zone is an L-shaped zone. (A formal definition of convex is that any straight line passing
through the zone intercepts at most two surfaces.) If the zone’s surfaces do not enclose a
space or if the zone is not convex you should use Solar Distribution = FullExterior instead of
FullInteriorAndExterior.
If you use FullInteriorAndExterior the program will calculate how much beam radiation
falling on an interior window is absorbed by the window, how much is reflected back into the
zone, and how much is transmitted into the adjacent zone. (Interior windows are assumed to
have no shading device).
If you use FullInteriorAndExterior the program will also calculate how much beam radiation
falling on the inside of an exterior window (from other windows in the zone) is absorbed by
the window, how much is reflected back into the zone, and how much is transmitted to the
outside. In this calculation the effect of an interior or exterior shading device, if present, is
accounted for.
FulInteriorAndlExteriorWithReflections
This case is the same interior distribution as the preceding option but uses exterior reflections
as well (see Solar Radiation Reflected from Exterior Surfaces for further explanation).
4/1/13 158
Climate, Sky and Solar/Shading Calculations Shading Module
4/1/13 159
Climate, Sky and Solar/Shading Calculations Shading Module
N surf
QS (ZoneNum)* A Q
i 1
i i SW ( ZoneNum)
where
i = zone surface number counter
Nsurf = number of heat transfer surfaces in zone
Ai = surface area [m2]
i = inside solar absorptance for an opaque surface, or, for a window, = back diffuse
transmittance plus back diffuse system absorptance of glass layers and shading device
layers (if present)
Solving this equation for QS gives:
Q ( ZoneNum)
QS ( ZoneNum) N surf
SW
QSW ( ZoneNum) *VMULT ( ZoneNum) (160)
AbsInsSurf A
i 1
i i
where
1
VMULT ( ZoneNum) N surf
[m -2 ]
AbsIntSurf * A
i 1
i i
Qsw is given by
QS SW QD( ZoneNum)
ZoneIntGain( ZoneNum)%QLTSW
ZoneIntGain( ZoneNum)%T _ QLTSW [W]
where
4/1/13 160
Climate, Sky and Solar/Shading Calculations Shading Module
where
BTOTZone = total beam solar incident on the zone’s exterior windows that is transmitted as
1
beam or diffuse.
BABSZone = total beam solar absorbed inside the zone.
BTOTZone is given by:
N extwin
BTOTZone TBmAll * SunlitFract * CosInc * Area * InOutprojSLFracMult
i 1
i i i i i
2
+ Diffuse entering zone from beam reflected by window inside reveal surfaces
+ Diffuse transmitted by windows from beam reflected by outside reveal surfaces
– Beam absorbed by window inside reveal surfaces
Here,
TBmAll = beam-to-beam plus beam-to-diffuse transmittance of window
SunlitFract = fraction of window irradiated by sun
CosInc = cosine of solar incidence angle on window
Area = glazed area of window [m2]
InOutProjSLFracMult = shadowing factor due to inside and outside projections of window
frame and/or divider (= 1.0 if there is no frame or divider).
BABSZone is given by the following sum (see Figure 48):
BABSZone = Beam absorbed by opaque inside surfaces3
+ Beam transmitted through the zone’s interior windows +
+ Beam transmitted back out of the zone’s exterior windows +
+ Beam absorbed by the zone’s exterior and interior windows +
+ Beam absorbed by inside daylighting shelves
1
For beam incident on an exterior window we have the following: For transparent glass with no shade or blind there is only
beam-to-beam transmission. For diffusing glass, or if a window shade is in place, there is only beam-to-diffuse transmission. If
a window blind is in place there is beam-to-diffuse transmission, and, depending on slat angle, solar profile angle, etc., there
can also be beam-to-beam transmission.
2
See “Beam Solar Reflection from Window Reveal Surfaces.”
3
If Solar Distribution = FullInteriorAndExterior in the Building object, the program calculates where beam solar from exterior
windows falls inside the zone. Otherwise, all beam solar is assumed to fall on the floor.
4/1/13 161
Climate, Sky and Solar/Shading Calculations Shading Module
Aoverlap(IW) Z1 B Z2
Aoverlap(B) IW
EW
Aoverlap(D)
Figure 48. Vertical section through a two-zone building showing where transmitted beam solar falls. Some
of the beam solar from exterior window EW is absorbed by the floor, D, interior wall, B, and interior
window, IW. Some is transmitted by IW to the adjacent zone, Z2. Aoverlap is the irradiated area of a
surface projected back onto the plane of EW. Beam reflected by D, B and IW contributes to the interior
short-wave radiation flux in Z1.
If zone ZoneNum shares interior windows with other zones, QS(ZoneNum) is modified to take
into account short-wave radiation received from the other zones through these windows:
QS ( ZoneNum) QS ( ZoneNum)
FractDifShortZtoZ (OtherZoneNum, ZoneNum)*
other
zones
where
FractDifShortZtoZ(OtherZoneNum,ZoneNum) = “diffuse solar exchange factor” = fraction of
short-wave radiation in OtherZoneNum that is transmitted to ZoneNum. This factor is
calculated in subroutine ComputeDifSolExcZonesWIZWindows taking into account multiple
reflection between zones. For example, for two zones means that some of the radiation
transmitted from Zone1 to Zone2 is reflected back to Zone1, and some of this is in turn
reflected back to Zone2, etc.
Interior Beam Radiation
4
The inside beam solar irradiance factor in (159) is given by:
4
For the purposes of the surface heat balance calculation, any beam solar radiation absorbed by a surface is assumed to be
uniformly distributed over the surface even though in reality it is likely to be concentrated in one or more discrete patches on
the surface.
4/1/13 162
Climate, Sky and Solar/Shading Calculations Shading Module
where
ldif,back = the system diffuse solar absorptance of layer l for irradiance from the back side
lbeam
,back = the system beam solar absorptance of layer l for irradiance from the back side
The interior diffuse short-wave radiation transmitted by an interior window to the adjacent
zone is given by
where
dif ( SurfNum) = diffuse transmittance of the interior window
The interior beam solar radiation transmitted by an interior window to the adjacent zone is
5
TBmi is zero if the window has diffusing glass or a shade. TBmi can be > 0 if a blind is present and the slat angle, solar
profile angle, etc., are such that some beam passes between the slats.
4/1/13 163
Climate, Sky and Solar/Shading Calculations Shading Module
N extwin
BeamSolarRad * beam
(SurfNum) TBmi * Aoverlapi (SurfNum) * CosInci [W]
i 1
where
beam
( SurfNum) is the beam-to-beam transmittance of the interior window at the
angle of incidence of beam solar from the exterior window on the interior window. The
program does not track where this radiation falls in the adjacent zone: it is counted as diffuse
radiation in that zone. Therefore,
Daylighting Ground Reflected Solar Modifier is used to modified the basic monthly ground
reflectance when snow is on the ground (from design day input or weather data values).
Values can range from 0.0 to 1.0.
References
ASHRAE. 2005. Handbook of Fundamentals, Chapter 31, Atlanta: ASHRAE.
ASHRAE. 2007. HVAC Applications, Chapter 33, Atlanta, ASHRAE.
Zhang, Qingyuan, Joe Huang, and Siwei Lang. 2002. "Development of Typical Year Weather
Data for Chinese Locations", American Society of Heating Refrigeration and Air-Conditioning
Engineers, ASHRAE Transactions, Vol 108, Part 2.
Threlkeld, J.L. and R.C. Jordan. 1958. Direct solar radiation available on clear days.
ASHRAE Transactions 64:45.
Groth, C. C., and Lokmanhekim, M. 1969. "Shadow - A New Technique for the Calculation of
Shadow Shapes and Areas by Digital Computer," Second Hawaii International Conference on
System Sciences, Honolulu, HI, January 22-24, 1969.
Walton, G.N. 1983. “The Thermal Analysis Research Program Reference Manual Program
(TARP)”, National Bureau of Standards (now National Institute of Standards and
Technology).
4/1/13 164
Climate, Sky and Solar/Shading Calculations Shading Module
4/1/13 165
Solar Radiation Reflected from Exterior Surfaces Shading Module
EnergyPlus has an option to calculate beam and sky solar radiation that is reflected from
exterior surfaces and then strikes the building. This calculation occurs if “withReflections” is
used on the SolarDistribution option in the Building object. For zones with detailed
6
daylighting, these reflections are also considered in the daylight illuminance calculations.
The reflecting surfaces fall into three categories:
Shadowing surfaces. These are surfaces like overhangs or neighboring buildings
entered with objects Shading:Site:Detailed, Shading:Building:Detailed, or
Shading:Zone:Detailed Examples are shown in Figure 49.
These surfaces can have diffuse and/or specular (beam-to-beam) reflectance values that
are specified with the ShadingProperty:Reflectance object.
Exterior building surfaces. In this case one section of the building reflects solar
radiation onto another section (and vice-versa). See Figure 50.
Opaque building surfaces (walls, for example) are assumed to be diffusely reflecting.
Windows and glass doors are assumed to be specularly reflecting. The reflectance
values for opaque surfaces are calculated by the program from the Solar Absorptance
and Visible Absorptance values of the outer material layer of the surface’s construction.
The reflectance values for windows and glass doors are calculated by the program from
the reflectance properties of the individual glass layers that make up surface’s
construction assuming no shading device is present and taking into account inter-
reflections among the layers.
The ground surface. Beam solar and sky solar reflection from the ground is calculated
even if “withReflections” is not used (the default). But in this case the ground plane is
considered unobstructed, i.e., the shadowing of the ground by the building itself or by
obstructions such as neighboring buildings is ignored. This shadowing is taken into
account only if “WithReflections” is used in the Solar Distribution field (in “Building” input
object) (Figure 51). In this case the user-input value of ground view factor is not used.
Sky diffuse
Beam
Beam
Figure 49. Examples of solar reflection from shadowing surfaces in the Shading series of input objects.
Solid arrows are beam solar radiation; dashed arrows are diffuse solar radiation. (a) Diffuse reflection of
beam solar radiation from the top of an overhang. (b) Diffuse reflection of sky solar radiation from the top
of an overhang. (c) Beam-to-beam (specular) reflection from the façade of an adjacent highly-glazed
building represented by a vertical shadowing surface.
6
A different method from that described here is used for calculating reflections from daylighting shelves (see “Daylighting
Shelves”).
4/1/13 166
Solar Radiation Reflected from Exterior Surfaces Diffuse Reflection of Beam Solar and Sky Solar Radiation
Wall
Window
Figure 50. Solar reflection from building surfaces onto other building surfaces. In this example beam
solar reflects from a vertical section of the building onto a roof section. The reflection from the window is
specular. The reflection from the wall is diffuse.
Beam Beam
Ground Ground
A B C
Figure 51. Shadowing by the building itself affects beam solar reflection from the ground. Beam-to-
diffuse reflection from the ground onto the building occurs only for sunlit areas, A and C, not for shaded
area, B. Shadowing by the building also affects sky solar reflection from ground (not shown).
A ray-tracing method is used to calculate beam solar and sky solar radiation that is diffusely
reflected onto each of a building’s exterior surfaces (walls, roofs, windows and doors), called
here “receiving surfaces.” The calculation begins by generating a set of rays proceeding into
the outward hemisphere at each receiving point on a receiving surface. Then it determinines
whether each ray hits the sky, ground or an obstruction. The radiance at the hit point from
reflection of incident beam or sky solar is determined and the contribution of this radiance to
the receiving surface is calculated, added to the contribution from other hit points, and
averaged over the receiving points. Separate calculations are done for beam-to-diffuse and
sky solar reflection from all obstructions and beam-to-diffuse and sky solar reflection from the
ground. (For beam-to-beam reflection see “Beam Solar Radiation Specularly Reflected from
Obstructions,” below.)
4/1/13 167
Solar Radiation Reflected from Exterior Surfaces Diffuse Reflection of Beam Solar and Sky Solar Radiation
Receiving points
An N-sided surface is assigned N receiving points with the following coordinates, expressed
in terms of the surface vertex coordinates:
N
Pij aki vkj , i 1,3; j 1, 3
k 1
where
Pij = jth coordinate of the ith receiving point
vkj = jth coordinate of the kth surface vertex
If N = 3: akj = 3/5 if k = i; = 1/5 otherwise
N 1 1
If N > 3: akj = if k = i; = otherwise
2N 2N
For example, for a vertical 3m by 5m rectangle with vertices (0,0,3), (0,0,0), (5,0,0) and
(5,0,3), this expression gives receiving points at (1.25,0,2.25), (1.25,0,0.75), (3.75,0,0.75)
and (3.75,0,2.25), as shown in Figure 52.
5m
Figure 52. Vertical rectangular exterior heat transfer surface showing location of receiving points for
calculating incident solar radiation reflected from obstructions.
Rays
A total of 90 rays are sent out into the exterior hemisphere surrounding each receiving point.
An upgoing ray may hit an obstruction or the sky. A downgoing ray may hit an obstruction or
the ground. See Figure 53.
In subroutine InitSolReflRecSurf, the following is determined for each ray, i, :
1. Unit vector in direction of ray
2. Cosine of angle between ray and plane of receiving surface ( cos i )
4/1/13 168
Solar Radiation Reflected from Exterior Surfaces Diffuse Reflection of Beam Solar and Sky Solar Radiation
5. For the surface containing the hit point: the surface number, the solar reflectance (
obs ,i if an obstruction), and the surface unit vector at the hit point pointing into the
hemisphere containing the receiving point
6. If the ray is downgoing and hits the ground, the coordinates of the ground hit point
7. Distance from receiving point to hit point
14 15
13
12
Rays from the
receiving point
11
Obstruction
10
Receiving surface
9
(window)
8 One of the
window’s
7 receiving points
Ground plane
6 5 4 3 21
Figure 53. Two-dimensional schematic showing rays going outward from a point on a receiving surface.
Rays 1-6 hit the ground, rays 7-11 hit an obstruction, and rays 12-15 hit the sky.
ReflFacSkySolObs(RecSurfNum)
N rec N ray
1
N rec
Hit
1 i 1
ViewFacSkyobs ,i DifShdgRatioIsoSkyobs ,i obs ,i cos i /
obs ,i
where
RecSurfNum is the receiving surface number,
N rec is the number of receiving points,
N ray is the number of rays,
“obs,i” denotes the obstruction hit by ray i,
Hitobs,i = 1 if ray i hits an obstruction, = 0 otherwise,
ViewFacSkyobs,i = unobstructed sky view factor of the obstruction = (1 cos tiltobs ) / 2 ,
DifShdgRatioIsoSkyobs,i = (obstructed sky irradiance on obstruction)/(unobstructed sky
irradiance on obstruction)
In this equation the product ViewFacSky*DifShdgRatioIsoSky is the sky irradiance at the hit
point divided by the horizontal sky irradiance taking into account shadowing of sky diffuse
radiation on the obstruction by other obstructions, and assuming that the radiance of the sky
is uniform. Note that we ignore secondary reflections here and in the following sections. In
4/1/13 169
Solar Radiation Reflected from Exterior Surfaces Diffuse Reflection of Beam Solar and Sky Solar Radiation
the present case this means that the irradiance at the hit point due to reflection of sky
radiation from the ground or from other obstructions is not considered.
The above reflection factor is used in the timestep calculation to find the irradiance on a
receiving surface due to sky radiation reflected from obstructions:
QRadSWOutIncSkyDiffReflObs(RecSurfNum) =
DifSolarRad * ReflFacSkySolObs(RecSurfNum) (W/m2)
Obstruction 5
11
4 Receiving surface
12 (window)
3
13 One of the
2 window’s
14
1 15 receiving points
Ground plane
One of the
ground hit
points
Figure 54. Two-dimensional schematic showing rays going upward from a ground hit point.
The factor for reflection of sky radiation from the ground onto a receiving surface is calculated
in subroutine CalcSkySolDiffuseReflFactors. It is given by:
ReflFacSkySolGnd(RecSurfNum)
1
N rec N ray N ray
N rec
1 ( Hit gnd ,i d i cos i / ) Hitsky , j ( i ) cos j (i ) d j (i ) /
i 1 j (i )
where
j (i ) denotes an upgoing ray from the ground point hit by ray i from the receiving point,
4/1/13 170
Solar Radiation Reflected from Exterior Surfaces Diffuse Reflection of Beam Solar and Sky Solar Radiation
QRadSWOutIncSkyDiffReflGnd(RecSurfNum) =
DifSolarRad * gnd * ReflFacSkySolGnd(RecSurfNum) (W/m2)
where gnd is the solar reflectance of the ground, which is assumed to be uniform over the
ground plane but may vary monthly (because of snow cover, for example).
Beam Solar Radiation Diffusely Reflected from Obstructions
This calculation is similar to that for sky solar reflected from obstructions. However, to find the
radiance at a hit point on an obstruction a line is drawn from the hit point to center of the sun.
From this line it is determined (1) if there is an obstruction between the hit point and the sun,
in which case it is assumed that no beam solar reaches the hit point; and (2) if beam solar
does reach the hit point, what the the solar angle of incidence at that point is.
The calculation is done for the hourly sun positions on each of the design days. It is also
done for hourly sun positions on selected days during the weather file run period (the same
days for which the shadowing calculations are done).
The factor for diffuse reflection of beam solar radiation from obstructions onto a receiving
surface is calculated in subroutine CalcBeamSolDiffuseReflFactors. It is given by:
N rec N ray
1
ReflFacBmToDiffSolObs(RecSurfNum,IHr)
N rec
Hit
1 i 1
obs ,i Hitobs ,i , sun d i cos i obs ,i cos sun ,obs ,i
where
QRadSWOutIncBmToDiffReflObs(RecSurfNum) = BeamSolarRad *
(WeightNow * ReflFacBmToDiffSolObs(RecSurfNum,HourOfDay) +
WeightPreviousHour * ReflFacBmToDiffSolObs(RecSurfNum,PreviousHour))
where BeamSolarRad is the timestep value of beam normal solar intensity (W/m2), and
WeightNow and WeightPreviousHour are time-averaging factors.
Beam Solar Radiation Diffusely Reflected from the Ground
This calculation is the same as that for beam solar diffusely reflected from obstructions
except that only rays from a receiving point that hit the ground are considered. The factor for
4/1/13 171
Solar Radiation Reflected from Exterior Surfaces Diffuse Reflection of Beam Solar and Sky Solar Radiation
diffuse reflection of beam solar from the ground onto a receiving surface is calculated in
subroutine CalcBeamSolDiffuseReflFactors. It is given by:
N rec N ray
1
ReflFacBmToDiffSolGnd(RecSurfNum,IHr)
N rec
Hit
1 i 1
gnd ,i Hit gnd ,i ,sun d i cos gnd ,i cos sun , gnd
where
Hit gnd ,i = 1 if ray i hits the ground, = 0 otherwise,
Hit gnd ,i , sun = 1 if the line from ray i’s hit point ot the sun is unobstructed, = 0 otherwise,
sun , gnd = angle of incidence of sun on ground (= solar zenith angle).
This factor is used in the timestep calculation to find the diffuse irradiance on a receiving
surface due to beam solar diffusely reflected from the ground:
Mirror image of
E Sun
sun about EF
B A
P
Q
R
C
J Receiving
surface
Specularly D
reflecting Receiving
S K
obstruction point
Figu