EnergyPlus™ Version 8.9.

0 Documentation

Engineering Reference

U.S. Department of Energy

February 27, 2018

Build: 1c5ba897d1
COPYRIGHT (c) 1996-2018 THE BOARD OF TRUSTEES OF THE UNIVERSITY OF ILLINOIS,
THE REGENTS OF THE UNIVERSITY OF CALIFORNIA THROUGH THE ERNEST ORLANDO
LAWRENCE BERKELEY NATIONAL LABORATORY, OAK RIDGE NATIONAL LABORATORY,
MANAGED BY UT-BATTELLE, ALLIANCE FOR SUSTAINABLE ENERGY, LLC, AND OTHER
CONTRIBUTORS. ALL RIGHTS RESERVED. NO PART OF THIS MATERIAL MAY BE REPRO-
DUCED OR TRANSMITTED IN ANY FORM OR BY ANY MEANS WITHOUT THE PRIOR WRIT-
TEN PERMISSION OF THE UNIVERSITY OF ILLINOIS OR THE ERNEST ORLANDO LAWRENCE
BERKELEY NATIONAL LABORATORY. ENERGYPLUS IS A TRADEMARK OF THE US DEPART-
MENT OF ENERGY.
Contents

2
Chapter 1

Overview

1.1 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.

1.2 General Modeling Overview

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.

1.3 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:

3
4 CHAPTER 1. OVERVIEW

Figure 1.1: EnergyPlus Program Schematic

Table 1.1: Simulation Flags

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.
SUBROUTINE ManageSimulation ! Main driver routine for this module
BeginSimFlag = .TRUE.
EndSimFlag = .FALSE.
CALL OpenOutputFiles
CALL GetProjectData
CALL GetEnvironmentInfo ! Get the number and type of Environments
DO Envrn = 1, NumOfEnvrn ! Begin environment loop ...
BeginEnvrnFlag = .TRUE.
EndEnvrnFlag = .FALSE.
WarmupFlag = .TRUE.
DayOfSim = 0
DO WHILE (( [Link]).OR.( WarmupFlag)) ! Begin day loop ...
DayOfSim = DayOfSim + 1
BeginDayFlag = .TRUE.
EndDayFlag = .FALSE.
DO HourOfDay = 1, 24 ! Begin hour loop ...
BeginHourFlag = .TRUE.
EndHourFlag = .FALSE.
DO TimeStep = 1, NumOfTimeStepInHour ! Begin time step (TINC) loop ...
BeginTimeStepFlag = .TRUE.
EndTimeStepFlag = .FALSE.
! Set the End_Flag variables to true if necessary. Note that each flag builds on
! the previous level. EndDayFlag cannot be .true. unless EndHourFlag is also .true.,
! etc.
! Note that the EndEnvrnFlag and the EndSimFlag cannot be set during warmup.
! Note also that BeginTimeStepFlag , EndTimeStepFlag , and the
! SubTimeStepFlags can/will be set/reset in the HVAC Manager.
IF (( [Link] )) THEN
1.4. WARMUP CONVERGENCE 5

CALL CloseOutputFiles
RETURN
END SUBROUTINE ManageSimulation

1.4 Warmup Convergence

At the beginning of each new environment (design day or runperiod) zone and surface temperatures
are initialized to 23C, and zone humidity ratios are initialized to the outdoor humidity ratio. The
first day of the environment is repeated until the loads/temperature convergence tolerance values
specified in the Building object are satisfied or until it reaches “maximum number of warmup days”.
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 intro-
duced 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 simula-
tion 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.
EnergyPlus determines warmup convergence in the following manner as shown in the Figure 1.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:

Tmax,prev − Tmax < Ttol (1.1)

Tmin,prev − Tmin < Ttol (1.2)

qh − qh,prev
< qtol (1.3)
qh

qc − qc,prev
< qtol (1.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
6 CHAPTER 1. OVERVIEW

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.
1.4. WARMUP CONVERGENCE 7

Figure 1.2: Flows of Warmup Convergence Checks

Chapter 2

Integrated Solution Manager

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 2.1.
Here, the operating point is at the intersection of the supply and demand curves.

Figure 2.1: Sequential Simulation Supply/Demand Relationship.

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.

8
2.1. BASIS FOR THE ZONE AND AIR SYSTEM INTEGRATION 9

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 simulta-
neous 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.

Figure 2.2: Schematic of Simultaneous Solution Scheme

In the sections which follow, the various individual functions of the integrated solution will be
described.

2.1 Basis for the Zone and Air System Integration

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.

dTz ∑ ∑
slN Nsurf aces N∑
zones

Cz = Q̇i + hi Ai (Tsi − Tz ) + ṁi Cp (Tzi − Tz ) + ṁinf Cp (T∞ − Tz ) + Q̇sys (2.1)

dt i=1 i=1 i=1

where:
∑
Nsl
Q̇i = sum of the convective internal loads
i=1
∑aces
Nsurf
hi Ai (Tsi − Tz ) = convective heat transfer from the zone surfaces
i=1
10 CHAPTER 2. INTEGRATED SOLUTION MANAGER

ṁ Cp (T∞ − Tz ) = heat transfer due to infiltration of outside air

∑inf
Nzones
i=1 ṁi Cp (Tzi − Tz ) = heat transfer due to interzone air mixing
Q̇sys = air systems output
Cz dT
dt
z
= energy stored in zone air
Cz = ρair Cp CT
ρair = zone air density
Cp = zone air specific heat
CT = sensible heat capacity multiplier (Detailed description is provided below)
If the air capacitance is neglected, the steady-state system output must be:

∑
Nsl ∑
Nsurf aces N∑
zones

− Q̇sys = Q̇i + hi Ai (Tsi − Tz ) + ṁi Cp (Tzi − Tz ) + ṁinf Cp (T∞ − Tz ) (2.2)

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, Qsys, can thus be formulated from the difference between the supply
air enthalpy and the enthalpy of the air leaving the zone as in Equation :

Q̇sys = ṁsys Cp (Tsup − Tz ) (2.3)

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 for Qsys in the heat balance Equation is shown in Equation :

∑
Nsl ∑aces
Nsurf ∑
Nzones
Cz dT z
= Q̇i + hi Ai (Tsi − Tz ) + ṁi Cp (Tzi − Tz )
dt
i=1 i=1 i=1 (2.4)
+ṁinf Cp (T∞ − Tz ) + ṁsys Cp (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 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:

dT ( )
= (δt)−1 Tz t − Tz t−δt + O (δt) (2.5)
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
2.1. BASIS FOR THE ZONE AND AIR SYSTEM INTEGRATION 11

Euler formula, Equation , was employed in Equation 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 , the formula for
updating the zone mean air temperature:
( )
Tzt −Tzt−δt ∑aces
Nsurf ∑
Nzones
Cz dt
+ Tzt hi Ai + ṁi Cp + ṁinf Cp + ṁsys Cp =
i=1 i=1
( )t−δt (2.6)
∑
N sl ∑
Nsurf aces ∑
Nzones
Q̇ti + ṁsys Cp Tsupply
t
+ hi Ai Tsi + ṁi Cp Tzi + ṁinf Cp T∞
i=1 i=1 i=1

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
∑
Nsl ∑aces
Nsurf ∑
Nzones
t
Q̇ti + ṁsys Cp Tsupply + Cz Tδtz + hi Ai Tsi + ṁi Cp Tzi + ṁinf Cp T∞
i=1 i=1 i=1
Tzt = ( ) (2.7)
Cz
∑aces
Nsurf ∑
Nzones
δt
+ hi Ai + ṁi Cp + ṁinf Cp + ṁsys Cp
i=1 i=1

Equation 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. How-
ever, 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. Approx-
imations 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:
( )
dTz −1 11 t 3 1 ( )
≈ (δt) Tz − 3Tzt−δt + Tzt−2δt − Tzt−3δt + O δt3 (2.8)
dt t 6 2 3

When this form for the derivative is used, equation changes to:

( 11 ) N∑sl ∑aces
Nsurf
Cz (δt)−1 T
6 z
t
− 3Tz
t−δt
+ 3 t−2δt
T
2 z
− 1 t−3δt
T
3 z
= Q̇ i + hi Ai (Tsi − Tz )
i=1 i=1 (2.9)
∑
Nzones
+ ṁi Cp (Tzi − Tz ) + ṁinf Cp (T∞ − Tz ) + ṁsys Cp (Tsup − Tz )
i=1

and the zone temperature update equation becomes:

12 CHAPTER 2. INTEGRATED SOLUTION MANAGER

∑
Nsl ∑aces
Nsurf N∑
zones (Cz)( )
Q̇i+ hiAiTsi+ ṁiCpTzi+ṁinf CpT∞+ṁsysCpTsupply− δt
−3Tzt−δt+32Tzt−2δt−13Tzt−3δt
Tzt=i=1 i=1 i=1
(11)Cz ∑aces
Nsurf N∑
zones

6 δt
+ hiA+ ṁiCp+ṁinf Cp+ṁsysC
i=1 i=1
(2.10)
This is the form historically used in EnergyPlus and is the current default referred to as 3rdOrder-
BackwardDifference in the ZoneAirHeatBalanceAlgorithm object. This algorithm requires zone air
temperatures at three previous time steps and uses constant temperature coefficients. The assump-
tion is that three previous time steps lengths are the same.
The AnalyticalSolution algorithm is an integration approach. While the 3rd 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

∑
N sl ∑aces
surf ∑
Nzones

 Q̇i + hi Ai Tsi + ṁi Cp Tzi +ṁinf Cp T∞ +ṁsys Cp Tsup


Tzt = Tzt−δt − i=1 i=1
Nsurf aces
i=1

∑ ∑
Nzones
hi Ai + ṁi Cp +ṁinf Cp +ṁsys Cp
 Nsurf aces
i=1 i=1 
∑ ∑
Nzones

 hi Ai + ṁi Cp +ṁinf Cp +ṁsys Cp


∗ exp − i=1 i=1
Cz
δt (2.11)

Nsurf aces
∑
N sl ∑ ∑
Nzones
Q̇i + hi Ai Tsi + ṁi Cp Tzi +ṁinf Cp T∞ +ṁsys Cp Tsup
+ i=1 i=1
Nsurf aces
i=1
∑ ∑
Nzones
hi Ai + ṁi Cp +ṁinf Cp +ṁsys Cp
i=1 i=1

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.

2.1.1 Zone Sensible Heat Capacity Multiplier

If the Zone Sensible Heat Capacity Multiplier = 1.0, this represents just the sensible heat capacitance
of the air volume in the specified zone. If the value is not defined, it is set to 1.0. This multiplier
can be greater than 1.0 if the zone air sensible heat capacity needs to be increased for stability
of the simulation. This multiplier increases the capacitance of the air volume by increasing the
zone volume that is used in the zone predictor-corrector algorithm in the simulation. This can be
done for numerical reasons, such as to increase the stability by decreasing the zone air temperature
2.2. SUMMARY OF PREDICTOR-CORRECTOR PROCEDURE 13

deviations at the time step level. Or it can be increased to try and account for the additional
capacitance in the air loop not specified in the zone, i.e. dampers, diffusers, duct work, etc., to
see the effect on the dynamics of the simulation. See the Input/Output Reference for additional
information (Object: ZoneCapacitanceMultiplier:ResearchSpecial).
In the source code below we see how the ZoneVolCapMultpSens increases the zone volume used
for the air ratio at the time step in the air system. This multiplier is constant throughout the
simulation.
AIRRAT(ZoneNum) = Zone(ZoneNum)%Volume \*** ZoneVolCapMultpSens **\* &
PsyRhoAirFnPbTdbW (OutBaroPress ,MAT(ZoneNum),ZoneAirHumRat(ZoneNum))\* &
PsyCpAirFnWTdb(ZoneAirHumRat(ZoneNum),MAT(ZoneNum))/( TimeStepSys \* SecInHour)

2.2 Summary of Predictor-Corrector Procedure

The predictor-corrector scheme can be summarized as follows:

• Using equation , an estimate is made of the air system energy required to balance the equation
with the zone air temperature equal to the setpoint temperature.
• With that quantity as a demand, the air system is simulated to determine its actual supply
capability at the time of the simulation. This will include a plant simulation if necessary.
• The actual air system capability is used in equation to calculate the resulting zone tempera-
ture.

2.3 Air System Control

Previously, the formulation of a new heat balance equation with an unsteady zone capacitance
term was discussed Equation . 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,
14 CHAPTER 2. INTEGRATED SOLUTION MANAGER

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 :

∑
Nsl ∑
Nsurf aces N∑
zones

Q̇load = Q̇i + hi Ai (Tsi − Tz ) + ṁi Cp (Tzi − Tz ) + ṁinf Cp (T∞ − Tz ) (2.12)

i=1 i=1 i=1

This is Equation without the term due to the air system. In addition, Tz 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. Qsys = Qload) 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 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 . 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 2.3. Simplified Single Zone Draw Through Air System.
The amount of heating or cooling provided by the air system in relation to the desired zone air
temperature is given by:

Q̇sys = ṁsys Cp η (Tsup − Tz, desired ) (2.13)

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 2.4. Simplified Variable Volume Air System. In VAV systems, the supply air
2.3. AIR SYSTEM CONTROL 15

Figure 2.3: Simplified Single Zone Draw Through Air System

temperature, as well as the supply air volume, are continuous functions of zone air temperature. As
shown in Figure 2.5. Idealized Variable Volume System Operation., when the 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.

Figure 2.4: Simplified Variable Volume Air System.

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.
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
16 CHAPTER 2. INTEGRATED SOLUTION MANAGER

Figure 2.5: Idealized Variable Volume System Operation.

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:

∑
Nsl ∑
Nsurf aces N∑
zones

Q̇0 = Q̇i + hi Ai Tsi + ṁi Cp Tzi + ṁinf Cp T∞ (2.14)

i=1 i=1 i=1

∑
Nsurf aces N∑
zones

Q̇slope = hi Ai + ṁi Cp + ṁinf Cp (2.15)

i=1 i=1

Equations and are derived, respectively, from the numerator and denominator of Equation 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 (2.16)
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
2.3. AIR SYSTEM CONTROL 17

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 (2.17)
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 and :
( )
Q̇sys,h = ηh Q̇h/c, max + Cp ρV̇min Tc/c − Tz,pred,heat (2.18)
( )( )
Q̇sys,c = Cp ρ ηc V̇max Tc/c − Tz,pred,cool (2.19)
Equation is valid for zone air temperatures below Th,upper , while Equation is valid for all tem-
peratures above this value. Equating the system output to the zone load, as given by Equation , 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,max Th,upper Cp ρV̇min Tc/c

Tz,pred,heat = + Q̇0 + (2.20)
Th,upper − Th,lower Q̇h/c,max
+ Cp ρV̇min + Q̇slope
Th,upper −Th,lower
√
B1 + B12 + B2
Tz,pred,cool = (2.21)
2
where,
ηc,min − C2
B1 = Tc/c + Tc,lower − (2.22)
C1
( ( ))
C3 ηc,min
B2 = 4 + Tc/c − Tc,lower (2.23)
C1 C1
and,
1 − ηc,min
C1 = (2.24)
Tc,upper − Tc,lower

Q̇slope
C2 = (2.25)
Cp ρV̇max
Q̇0
C3 = (2.26)
Cp ρV̇max
Once the predicted zone air temperature has been calculated from Equations and , 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:
18 CHAPTER 2. INTEGRATED SOLUTION MANAGER

V̇sup ply = ηc V̇max (2.27)

where the supply air fraction η c is computed from Equation . 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 + (2.28)
Cp ρV̇min
The reheat coil capacity fraction η h is determined by using Equation . Once Equation or ,
has been used, the supply air flow rate and temperature are known. These values are then used
in Equation 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.

2.4 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:

∑
Nsl ∑aces
Nsurf ∑
Nzones
ρair Vz CW dW
dt
z
= kgmasssched load + Ai hmi ρairz (Wsurf si − Wzt ) + ṁi (Wzi − Wzt )
i=1 i=1 i=1
+ṁinf (W∞ − Wzt ) + ṁsys (Wsup − Wzt )
(2.29)
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)−1 (Wzt − Wzt−δt ) + O(δt) (2.30)
dt
The zone air humidity ratio update at the current time step using the EulerMethod may be
expressed as follows:

( ) N∑sl ∑aces
Nsurf
ρair V CW (δt)−1 Wzt − Wzt−δt = kgmasssched load + Ai hmi ρairz (Wsurf si − Wzt )
i=1 i=1 (2.31)
∑
Nzones
+ ṁi (Wzi − Wzt ) + ṁinf (W∞ − Wzt ) + ṁsys (Wsup − Wzt )
i=1
2.4. MOISTURE PREDICTOR-CORRECTOR 19

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:
( 11 )
dWz W t
z − 3Wz
t−δt
+ 3
Wz
t−2δt
− f rac13Wz
t−3δt ( )
≈ 6 2
+ O δt3 . (2.32)
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.

ρair Vz CW
( 11 ) ∑aces
Nsurf ∑
Nzones ∑
N sl

δt 6
Wzt + Ai hmi ρairz Wzt + ṁi Wzt + ṁinf Wzt + ṁsys Wzt = kgmasssched load
i=1 i=1 i=1
∑aces
Nsurf ∑
Nzones
+ Ai hmi ρairz Wsurf si + ṁi Wzi + ṁinf W∞ + ṁsys Wsup
i=1 ( i=1 )
− ρair Vδtz CW −3Wzt−δt + 32 Wzt−2δt − 12 Wzt−3δt
(2.33)
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 Analyt-
icalSolution 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.

2.4.1 Moisture Prediction

For the moisture prediction case the equation is solved for the anticipated system response as shown
below.

PredictedSystemLoad = ṁsys ∗ (Wzt − Wsup )

(2.34)
MassFlow ∗ HumRat = kgsecair kgkgwater
air
= kgsec
water

Since the program provides three solution algorithms, the moisture prediction from each solution
algorithm is given below.

[Link] EulerMethod
For this solution algorithm, the air mass balance for the predicted air system load or response is:
20 CHAPTER 2. INTEGRATED SOLUTION MANAGER

( )
PredictedSystemLoad[kgW ater/sec]=ρairVzCW (δt)−1 Wsetpoint
t −Wzt−δt
[ ]
∑
N sl ∑
Nsurf aces ( ) Nzones
∑ ( ) ( )
− kgmassschedload+ Aihmiρairz Wsurf si−Wsetpoint +
t ṁi Wzi−Wsetpoint +ṁinf W∞−Wsetpoint
t t (2.35)
i=1 i=1 i=1

[Link] ThirdOrderBackwardDifference
For this solution algorithm, the air mass balance for the predicted system load or response is given
below:
( )
( ) Nsurf ∑aces ∑
Nzones
PredictedSystemLoad[kgW ater/sec]= ρairV z CW 11
δt 6
+ A h ρ
i mi airz + ṁ i+ ṁinf Wz −
t

[ i=1 i=1
] (2.36)
∑
N sl ∑
Nsurf aces
∑
Nzones ( )
kgmassschedload+ AihmiρairzWsurf si+ ṁiWzi+ṁinf W∞+ρairV z CW
δt
3W t−δt 3
z − W
2 z
t−2δt 1
+ W
3 z
t−3δt
i=1 i=1 i=1

Then, using the following substitutions, the air mass balance equation becomes:

∑
Nsurf aces N∑
zones

A= Ai hmi ρairz + ṁi + ṁinf (2.37)

i=1 i=1

∑
Nsl ∑
Nsurf aces N∑
zones

B= kgmasssched load + Ai hmi ρairz Wsurf si + ṁi Wzi + ṁinf W∞ (2.38)

i=1 i=1 i=1

ρair Vz CW
C= (2.39)
δt
[ ]
PredictedSystemLoad [kgW ater / sec] = 11 ∗ C + A ∗ WSetP oint −
[ ( 6 )] (2.40)
B + C ∗ 3Wzt−δt − 32 Wzt−2δt + 31 Wzt−3δt

[Link] AnalyticalSolution
For this solution algorithm, the air mass balance for the predicted air system load or response is
given below:
[ ]
∑aces
Nsurf ∑
Nzones
PredictedSystemLoad [kgW ater / sec] = Ai hmi ρairz + ṁi Wzi + ṁinf ∗
  N i=1 i=1

∑aces
surf ∑
Nzones

 t  Ai hmi ρairz + ṁi Wzi +ṁinf


Wsetpoint − Wzt−δt ∗ exp − i=1
ρair Vz CW
i=1
δt ∗

  Nsurf aces
−1 (2.41)
∑ ∑
Nzones

  Ai hmi ρairz + ṁi +ṁinf


1 − exp − i=1
ρair Vz CW
i=1
δt −
( )
∑
N sl ∑aces
Nsurf ∑
Nzones
kgmasssched load Ai hmi ρairz Wsurf si + ṁi Wzi + ṁinf W∞
i=1 i=1 i=1
2.4. MOISTURE PREDICTOR-CORRECTOR 21

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.

2.4.2 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 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.

[Link] EulerMethod
∑
Nsl ∑aces
Nsurf N∑
zones t−δt
kgmassschedload+ AihmiρairzWsurf si+ ṁiWzi+ṁinf W∞+ṁsysWsup+ρairVzCW Wzδt
Wzt=i=1 i=1 i=1

ρairVzCW ∑aces
Nsurf N∑
zones

δt
+ Aihmiρairz+ ṁi+ṁinf +ṁsys
i=1 i=1
(2.42)

[Link] ThirdOrderBackwardDifference
∑
N sl ∑
Nsurf aces
∑
Nzones
kgmassschedload+ AihmiρairzWsurf si+ ṁiWzi+ṁinf W∞+ṁsysWsup+ρairVz CW
δt
(3Wzt−δt−23Wzt−2δt+13Wzt−3δt)
i=1 i=1 i=1
Wzt=
ρairVz CW (11) ∑ ∑
Nsurf aces Nzones
δt 6
+ Aihmiρairz+ ṁi+ṁinf +ṁsys
i=1 i=1
(2.43)

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:
If (ZoneSupplyAirMassFlowRate > 0.0) Then

∑
Nsurf aces N∑
zones

A= Ai hmi ρairz + ṁi + ṁinf + ṁsys (2.44)

i=1 i=1

∑
Nsl ∑
Nsurf aces N∑
zones

B= kgmasssched load + Ai hmi ρairz Wsurf si + ṁi Wzi + ṁinf W∞ + ṁsys Wsup (2.45)
i=1 i=1 i=1
22 CHAPTER 2. INTEGRATED SOLUTION MANAGER

ρair Vz CW
C= (2.46)
δt
Else If (ZoneSupplyAirMassFlowRate < = 0.0) Then

∑
Nsurf aces N∑
zones

A= Ai hmi ρairz + ṁi + ṁinf + ṁExhaust (2.47)

i=1 i=1

∑
Nsl ∑
Nsurf aces N∑
zones

B= kgmasssched load + Ai hmi ρairz Wsurf si + ṁi Wzi + ṁinf W∞ + ṁExhaust W∞ (2.48)
i=1 i=1 i=1

ρair Vz CW
C= (2.49)
δt
End If
Inserting in the parameters A, B and C above in the air mass balance equation, it simplifies to:
[ ( )]
B + C ∗ 3W t−δt
− 3
W t−2δt
+ 1
W t−3δt
Wzt = z
( 11 ) 2 z 3 z
(2.50)
6
∗ C + A

[Link] AnalyticalSolution
 Nsurf aces

∑
N sl ∑ ∑
Nzones

 kgmasssched load + Ai hmi ρairz Wsurf si + ṁi Wzi +ṁinf W∞ +ṁsys Wsup

Wzt =  Wzt−δt − i=1 i=1
Nsurf aces
i=1
∗
∑ ∑
Nzones
Ai hmi ρairz + ṁi +ṁinf +ṁsys
 Nsurf aces
i=1  i=1

∑ ∑
Nzones

 Ai hmi ρairz + ṁi +ṁinf +ṁsys


exp − i=1 i=1
ρair Vz CW
δt +

Nsurf aces
∑
N sl ∑ ∑
Nzones
kgmasssched load + Ai hmi ρairz Wsurf si + ṁi Wzi +ṁinf W∞ +ṁsys Wsup
i=1 i=1 i=1
Nsurf aces
∑ ∑
Nzones
Ai hmi ρairz + ṁi +ṁinf +ṁsys
i=1 i=1
(2.51)
The above solutions are implemented in the Correct Zone Air Humidity Ratio step in Ener-
gyPlus. 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 HeatBalanceAl-
gorithm 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.
2.5. CARBON DIOXIDE PREDICTOR-CORRECTOR 23

2.4.3 Which moisture buffering model is best?

The ’correct’ moisture buffering model depends on the questions being answered by the building
energy simulation. Previous research (Woods et al., 2013a) has shown that using the effective
capacitance model to account for moisture buffering of materials will provide a good estimate of
energy use when humidity is not being actively controlled. See the InputOutput Reference for
additional information on the object ZoneCapacitanceMultiplier:ResearchSpecial. This model has
some limitations (Woods et al., 2013b):

• it will not accurately predict indoor humidity (or thermal comfort),

• it will not accurately predict energy use when humidity is being actively controlled, and

• it will not provide insight into the moisture content and potential moisture problems associated
with a specific wall construction.

The effective moisture penetration depth (EMPD) model will address the first two concerns
above: it can accurately predict indoor humidity, and can accurately predict energy use associated
with controlling humidity. The EMPD model requires more user input than the effective capacitance
model, specifically some of the moisture properties of the materials in the building. For more
information, see the Effective Moisture Penetration Depth Model section in this document.
Like the EMPD model, the combined heat, air, and moisture transfer (HAMT) model addresses
the first two issues discussed above for the effective capacitance model. It also addresses the third,
by providing temperature and moisture profiles through composite building walls, and helping to
identify surfaces with high surface humidity. The HAMT model requires a few more user inputs on
moisture properties of materials than the EMPD model, and this model also increases the required
simulation time by an order of magnitude. For more information on this model, see the Combined
Heat and Moisture Transfer (HAMT) Model section in this document.
Note that the EMPD and HAMT models above ensure accurate calculations of the effect of
moisture buffering, but it will only be accurate relative to reality when given appropriate inputs for
the material properties.
Woods, J., J. Winkler, D. Christensen, Moisture modeling: Effective moisture penetration depth
versus effective capacitance, in Thermal Performance of the Exterior Envelopes of Whole Buildings
XII International Conference. 2013a: Clearwater, FL.
Woods, J., Winkler, J, and Christensen, D. Evaluation of the Effective Moisture Penetration
Depth Model for Estimating Moisture Buffering in Buildings, NREL/TP-5500-57441, 2013b.

2.5 Carbon Dioxide Predictor-Corrector

The transient air mass balance equation for the change in zone air carbon dioxide concentration
may be expressed as follows:

t ∑
Nsl ∑
Nzones
ρair Vz CCO2 dC z
= kgmasssched load ∗ 1.06 + ṁi (Czi − Czt ) + ṁinf (C∞ − Czt )
dt
i=1 i=1 (2.52)
+ṁsys (Csup − Czt )
24 CHAPTER 2. INTEGRATED SOLUTION MANAGER

where:
∑
Nsl
kgmasssched load = sum of scheduled internal carbon dioxide loads. The zone air density is used
i=1
to convert the volumetric rate of carbon dioxide generation from user input into mass generation
rate∑ [kg/s].The coefficient of 106 is used to make the units of carbon dioxide as ppm.
Nzones
i=1 ṁi (Czi − Czt ) = carbon dioxide transfer due to interzone air mixing [ppm-kg/s]
Czi = carbon dioxide concentration in the zone air being transferred into this zone [ppm]
ṁinf (C∞ − Czt ) = carbon dioxide transfer due to infiltration and ventilation of outdoor air [ppm-
kg/s]
C∞ = carbon dioxide concentration in outdoor air [ppm]
ṁsys (Csup − Czt ) = carbon dioxide transfer due to system supply [ppm-kg/s]
Csup = carbon dioxide concentration in the system supply airstream [ppm]
ṁsys = air system supply mass flow rate [kg/s]
t
ρair Vz dCdt
z
= carbon dioxide storage term in zone air [kg/s]
t
Cz = 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)−1 (Czt − Czt−δt ) + O(δt) (2.53)
dt
The zone air carbon dioxide concentration update at the current time step using the Euler-
Method may be expressed as follows:

( ) N∑sl ∑
Nzones
ρair VZ CCO2 (δt)−1 Czt − Czt−δt = kgmasssched load ∗ 106 + ṁi (Czi − Czt )
i=1 i=1 (2.54)
+ṁinf (C∞ − Czt ) + ṁsys (Csup − Czt )

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 calculations. This
algorithm is the default choice and is defined as ThirdOrderBackwardDifference in the ZoneAirHeat-
BalanceAlgorithm object.
The third order derivative derived from a Taylor Series expansion is defined as:

dCzt C − 3Czt−δt + 32 Czt−2δt − 13 Czt−3δt

11 t
≈ 6 z
+ O(δt3 ) (2.55)
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
2.5. CARBON DIOXIDE PREDICTOR-CORRECTOR 25

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.

( 11 ) ∑
Nzones ∑
N sl
ρair Vz CCO2
δt 6
Czt + ṁi Czt + ṁinf Czt + ṁsys Czt = kgmasssched load ∗ 106
i=1 i=1
(2.56)
∑
Nzones
+ ṁi Czi + ṁinf C∞ + ṁsys Csup − ρair Vz CCO2
δt
(−3Czt−δt + 3 t−2δt
C
2 z
− 1 t−3δt
C
3 z
)
i=1

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 Analyt-
icalSolution 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.

2.5.1 Carbon Dioxide Prediction

For the carbon dioxide concentration prediction case, the equation is solved for the anticipated
system response as shown below.

( )
PredictedSystemLoad = ṁsys Csup − Czt (2.57)

Since the program provides three solution algorithms, the carbon dioxide prediction from each
solution algorithm is given below.

[Link] EulerMethod

For this solution algorithm, the air mass balance for the predicted air system load or response is:

( t )
PredictedSystemLoad [kg/ sec] = ρair VZ CCO2 (δt)−1 Csetpoint − Czt−δt
[N ]
∑sl ∑
Nzones ( ) ( )
− kgmasssched load ∗ 10 +
6
ṁi Czi − Csetpoint + ṁinf C∞ − Csetpoint
t t (2.58)
i=1 i=1

[Link] ThirdOrderBackwardDifference

For this solution algorithm, the air mass balance for the predicted system load or response is given
below:
26 CHAPTER 2. INTEGRATED SOLUTION MANAGER

[ ]
( 11 ) ∑
Nzones
PredictedSystemLoad [kg/ sec] = ρair Vz CCO2
6
+δt
ṁi + ṁinf ∗ Csetpoint
t

[N i=1 ]
∑sl ∑
Nzones ( t−δt 3 t−2δt 1 t−3δt )
− kgmasssched load ∗ 10 +
6
ṁi Czi + ṁinf C∞ + ρair Vz CCO2
δt
3Cz − 2 Cz + 3 Cz
i=1 i=1

(2.59)

[Link] AnalyticalSolution

For this solution algorithm, the air mass balance for the predicted air system load or response is
given below:

  Nzones 
[N ] ∑
∑
zones ṁi +ṁ inf
PredictedSystemLoad[kg/sec]= ṁi+ṁinf ∗Csetpoint
t
−Czt−δt∗exp− ρi=1
air VZ CCO2
δt∗
i=1
  −1
∑
Nzones (N )
ṁi+ṁinf ∑sl N∑
zones
1−exp− i=1
δt − kgmassschedload∗10 + 6
ṁiCzi+ṁinf C∞
ρairVZ CCO2
i=1 i=1

(2.60)
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 = ṁsys Csup − Czt ≈ ṁOA,z C∞ − Csetpoint
t
(2.61)

where:
ṁ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.

2.5.2 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 pre-
dicting 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.
2.6. GENERIC CONTAMINANT PREDICTOR-CORRECTOR 27

[Link] EulerMethod
∑
Nsl ∑
Nzones t−δt
kgmasssched load ∗ 106 + ṁi Czi + ṁinf C∞ + ṁsys Csup + ρair VZ CCO2 Czδt
i=1 i=1
Czt = (2.62)
ρair VZ CCO2 ∑
Nzones
δt
+ ṁi + ṁinf + ṁsys
i=1

[Link] ThirdOrderBackwardDifference
∑
Nsl N∑
zones
kgmassschedload∗106+ ṁiCzi+ṁinf C∞+ṁsysCsup+ρairVδt
Z CCO2
(3Czt−δt−23Czt−2δt+13Czt−3δt)
Czt =i=1 i=1
( ) N∑
ρairVZ CCO2 11
zones

δt 6
+ ṁi+ṁinf +ṁsys
i=1
(2.63)

[Link] AnalyticalSolution
   
∑
N sl ∑
Nzones ∑
Nzones
kgmasssched load ∗106 + ṁi Czi +ṁinf C∞ +ṁsys Csys ṁi +ṁinf +ṁsys
Czt = Czt−δt − i=1
∑
Nzones
i=1  ∗ exp − i=1
ρair VZ CCO2
δt
ṁi +ṁinf +ṁsys
i=1
∑
N sl ∑
Nzones
kgmasssched load ∗106 + ṁi Czi +ṁinf C∞ +ṁsys Csys
+ i=1 ∑
Nzones
i=1

ṁi +ṁinf +ṁsys

i=1
(2.64)
The above solutions are implemented in the Correct Zone Air Carbon Dioxide step in the Zone
Contaminant Predictor Corrector module of EnergyPlus.

2.6 Generic Contaminant Predictor-Corrector

The transient air mass balance equation for the change in zone air generic contaminant concentration
may be expressed as follows:

t
dCf,z ∑
Nsource N∑
sink ∑
Nzones ( )
ρair Vz Mf or dt
= ρair Gf,i ∗ 1.06 −ρair Rf,i Cf,z + ṁi Cf,z,i − Cf,z
t

( i=1 ) ( i )
∑ i=1
C
(2.65)
+ṁinf Cf,∞ − Cf,z + ṁsys Cf,sup − Cf,z + hj ρAj ( ks,j
t t
j
− Cf,z ) + Sf (Cf,z
t−δt
)
j

where:
∑
Nsource
ρair Gf,i = Sum of internal generic contaminant loads from sources in a zone or interior
i=1
surfaces.
The zone air density is used to convert the volumetric rate of generic contaminant generation
from user input into mass generation rate [kg/s].The coefficient of 106 is used to make the units of
generic contaminant as ppm.
N∑
sink
ρair Rf,i Cf,z = Sum of removal rate from sinks in a zone or interior surfaces [ppm-kg/s]
i
28 CHAPTER 2. INTEGRATED SOLUTION MANAGER

∑Nzones ( )
ṁi Cf,z,i − Cf,z
i=1
t
= Generic contaminant transfer due to interzone air mixing [ppm-kg/s]
Cf,z,i(= Generic contaminant
) concentration in the zone air being transferred into this zone [ppm]
ṁinf Cf,∞ − Cf,z = Generic contaminant transfer due to infiltration and ventilation of outdoor
t

air [ppm-kg/s]
Cf,∞ (= Generic contaminant
) concentration in outdoor air [ppm]
ṁsys Cf,sup − Cf,z = Generic contaminant transfer due to system supply [ppm-kg/s]
t

Cf,sup = Generic contaminant concentration in the system supply airstream [ppm]

ṁsys = Air system supply mass flow rate [kg/s]
dC t
ρair Vz dtf,z = Generic contaminant storage term in zone air [ppm-kg/s]
t
Cf,z = Zone air generic contaminant concentration at the current time step [ppm]
ρair = Zone air density [kg/m3 ]
Vz = Zone volume [m3 ]
∑ C
hj ρAj ( ks,j
j
− Cf,z ) = Generic contaminant transport through diffusion between interior sur-
j
faces and zone air
t−δt
Sf (Cf,z ) = 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:
t
dCf,z
= (δt)−1 (Cf,z
t
− Cf,z
t−δt
) + O(δt) (2.66)
dt
The zone air generic contaminant concentration update at the current time step using the
EulerMethod may be expressed as follows:

( t ) Nsource
∑ N∑
sink N∑
zones ( )
ρairVzMf or(δt)−1 Cf,z −Cf,z
t−δt
= ρairGf,i∗1.06−ρair Rf,iCf,z+ ṁi Cf,z,i−Cf,z
t

( ) ( i=1 )
∑ Cs,j
i i=1 (2.67)
+ṁinf Cf,∞−Cf,z +ṁsys Cf,sup−Cf,z + hj ρAj ( kj −Cf,z)+Sf
t t
j

To preserve the stability of the calculation of the zone generic contaminant 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 calcula-
tions. This algorithm is the default choice and is defined as ThirdOrderBackwardDifference in the
ZoneAirHeatBalanceAlgorithm object.
The third order derivative resulting from a Taylor Series expansion is defined as:

t
dCf,z 11 t
C
6 f,z
− 3Cf,z
t−δt t−2δt
+ 32 Cf,z − 13 Cf,z
t−3δt
≈ + O(δt3 ) (2.68)
dt δt
2.6. GENERIC CONTAMINANT PREDICTOR-CORRECTOR 29

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.

ρair Vz Mf or ( 11 ) t
N∑
sink
t
∑
Nzones
t t t
∑ t
δt 6
Cf,z + ρair Rf,i Cf,z + ṁi Cf,z + ṁinf Cf,z + ṁsys Cf,z + hj ρAj Cf,z
i i=1 j
∑
Nsource ∑
Nzones ∑ Cs,j
= ρair Gf,i ∗ 1.06 + ṁi Cf,z,i + ṁinf Cf,∞ + ṁsys Cf,sup + hj ρAj kj
+ Sf
i=1
ρ V M ( i=1
) j
− air δtz f or −3Cf,zt−δt t−2δt
+ 32 Cf,z − 31 Cf,z
t−3δt

(2.69)
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 Analyt-
icalSolution 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.

2.6.1 Generic Contaminant Prediction

For the generic contaminant concentration prediction case, the equation is solved for the anticipated
system response as shown below.
( )
PredictedSystemLoad = ṁsys Cf,sup − Cf,z
t
(2.70)

Since the program provides three solution algorithms, the generic contaminant prediction from
each solution algorithm is given below.

[Link] EulerMethod

For this solution algorithm, the air mass balance for the predicted air system load or response is:

( )
PredictedSystemLoad [kg/ sec] = ρair Vz Mf or (δt)−1 Csetpoint − Cf,z t−δt
[N ]
∑
source N∑
sink ∑
Nzones
− ρair Gf,i ∗ 1.0 −ρair
6
Rf,i Csetpoint + ṁi (Cf,z,i − Csetpoint ) + Sf
[ i=1 i i=1 ]
(2.71)
∑ C
− ṁinf (Cf,∞ − Csetpoint ) + hj ρAj ( ks,j j
− Csetpoint )
j
30 CHAPTER 2. INTEGRATED SOLUTION MANAGER

[Link] ThirdOrderBackwardDifference
For this solution algorithm, the air mass balance for the predicted system load or response is given
below:

[ ]
ρ Vz Mf or(11) N∑
sink ∑
Nzones
PredictedSystemLoad[kg/sec]= air δt 6
Csetpoint− ρair Rf,iCsetpoint+ ṁiCsetpoint+ṁinf Csetpoint
[ i i=1
]
∑ ∑
Nsource ∑
Nzones ∑ C (2.72)
− hj ρAj Csetpoint+ ρairGf,i∗1.06+ ṁiCf,z,i+ṁinf Cf,∞+ hj ρAj ks,j +Sf
j
ρ
j
Vz Mf or
( i=1
) i=1 j
t−δt 3 t−2δt 1 t−3δt
+ air δt −3Cf,z +2Cf,z −3Cf,z

[Link] AnalyticalSolution
For this solution algorithm, the air mass balance for the predicted air system load or response is
given below:

[ ]
∑
Nzones N∑
sink ∑
PredictedSystemLoad [kg/ sec] = ṁi + ṁinf + ρair Rf,i + hj ρAj ∗
  i=1 i

j
∑
Nzones N∑
sink ∑
ṁi +ṁinf +ρair Rf,i + hj ρAj
 t  
Csetpoint − Czt−δt ∗ exp − δt ∗
i=1 i j
ρair VZ MF OR

  −1 (2.73)
∑
Nzones N∑
sink ∑
ṁi +ṁinf +ρair Rf,i + hj ρAj
  
1 − exp − δt −
i=1 i j
ρair VZ MF OR

( )
∑
Nsource ∑
Nzones ∑ Cs,j
ρair Gf,i ∗ 1.06 + ṁi Cf,z,i + ṁinf Cf,∞ + hj ρAj kj
+ Sf
i=1 i=1 j

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:
( ) ( )
PredictedSystemLoad = ṁsys Cf,sup − Cf,z
t
≈ ṁOA,z Cf,∞ − Csetpoint
t
(2.74)
where:
ṁOA,z = Supply outdoor airflow rate into the controlled zone [kg/s]
The above approximation is based on the assumption that the generic contaminant concentration
at the outdoor air (OA) mixer inlet is equal to the zone air outlet concentration level, and the generic
contaminant 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
generic contaminant correction step, which is shown next.
2.7. ZONE AIR MASS FLOW CONSERVATION 31

2.6.2 Generic Contaminant Correction

For the correct step, the expanded air mass balance equation is solved for the final zone generic
contaminant 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.

[Link] EulerMethod
∑
Nsource ∑
Nzones ∑ C C t−δt
ρair Gf,i ∗1.06 + ṁi Cf,z,i +ṁinf Cf,∞ +ṁsys Cf,sup + hj ρAj ks,j +ρair VZ MF OR f,z
δt
+Sf
j
t
Cf,z = i=1 i=1
N∑
sink ∑
Nzones
j
∑
(2.75)
ρair Vz Mf or (δt)−1 +ρair Rf,i + ṁi +ṁinf +ṁsys + hj ρAj
i i=1 j

[Link] ThirdOrderBackwardDifference
∑
Nsource ∑
Nzones
ρair Gf,i ∗ 1.06 + ṁi Cf,z,i + ṁinf Cf,∞ + ṁsys Cf,sup
i=1
∑ Cs,j
i=1
+ hj ρAj kj
+ ρair VZ MF OR
δt
t−δt
(3Cf,z − 32 Cf,z
t−2δt t−3δt
+ 13 Cf,z ) + Sf (2.76)
t j
Cf,z = N∑
sink ∑
Nzones ∑
ρair Vz Mf or (δt)−1 ( 11
6 )
+ρair Rf,i + ṁi +ṁinf +ṁsys + hj ρAj
i i=1 j

[Link] AnalyticalSolution
 
∑
Nsource ∑
Nzones ∑ C
ρair Gf,i ∗1.06 + ṁi Cf,z,i +ṁinf Cf,∞ +ṁsys Cf,sys + hj ρAj ks,j
 t−δt j 
t
Cf,z = Cf,z − i=1
∑
Nzones
i=1
N∑
sink ∑
j
∗
ṁi +ṁinf +ṁsys +ρair Rf,i + hj ρAj
 i=1 i
 j

∑
Nzones N∑
sink ∑
ṁi +ṁinf +ṁsys +ρair Rf,i + hj ρAj
 
exp − i=1
ρair VZ MF OR
i j
δt + (2.77)

∑
Nsource ∑
Nzones ∑ C
ρair Gf,i ∗1.06 + ṁi Cf,z,i +ṁinf Cf,∞ +ṁsys Cf,sys + hj ρAj ks,j +Sf
i=1 i=1 j j

∑
Nzones N∑
sink ∑
ṁi +ṁinf +ṁsys +ρair Rf,i + hj ρAj
i=1 i j

The above solutions are implemented in the Correct Zone Air Generic Contaminant step in the
Zone Contaminant Predictor Corrector module of EnergyPlus.

2.7 Zone Air Mass Flow Conservation

2.7.1 Overiew
The zone air mass flow conservation object, ZoneAirMassFlowConservation, activates zone
air mass flow balance calculations. This feature is available only for controlled zones
(ZoneVAC:EquipmentConnections) which also have either a zone mixing or infiltration ob-
ject. The user may specify that zone mixing, infiltration, or both can be overridden to balance the
zone air mass flows. The following rules apply:
32 CHAPTER 2. INTEGRATED SOLUTION MANAGER

• If there are no zone mixing flows to adjacent zones, then the zone air mass flow is balanced
by setting the Zone Mixing objects mass flow rate to zero.
• If there are no zone exhaust fans defined and there are no zone mixing objects specified, then
a zone in an air loop is always balanced.
• Infiltration mass flow is included in the zone air mass flow balance depending upon one of
three options: none (all infiltration is assumed to be self-balanced), all zones, or only zones
that serve as a source zone for zone mixing objects.
• The base infiltration mass flow rate (calculated based on user inputs) may be controlled one
of two ways for zone air mass flow balance purposes: adjust the base infiltration up or down as
needed to balance the zone air mass flow, or assume the base infiltration rate is self-balanced
and add infiltration if needed to balance the zone air mass flow.
• Optional user inputs can override the default return air flow rate.

The zone air mass flow conservation equation always includes: supply air flow rates, return air
flow rates, and zone exhaust fan flow rates. Zone mixing and infiltration object flow rates may be
included depending upon the selected options. A particular zone can be a source zone, receiving
zone, or both depending on the number of ZoneMixing objects specified for that zone.

2.7.2 Return Air Flow Rate Calculations

The return air flow rate is calculated one of several ways, depending on whether ZoneAirMassFlow-
Conservation is active, if there is more than one return node, and if Zone Return Air Node 1 Flow
Rate Fraction Schedule Name or Zone Return Air Node 1 Flow Rate Basis Node or NodeList Name
has been specified.

[Link] Base Zone Total Return Flow

The total zone return flow calculation with no ZoneAirMassFlowConservation is the sum of the
zone inlets less the sum of the zone exhausts adjusted for any balanced exhaust fan flow:

ṁR = M AX (0.0, ṁS − [ṁEX,tot − ṁEXF,bal ]) (2.78)

If ZoneAirMassFlowConservation is active, then the total zone return air flow rate is

ṁR = M AX (0.0, ṁS − [ṁEX,tot − ṁEXF,bal ] + [ṁXR − ṁXS ]) (2.79)

[Link] Zone Return Air Node 1

The first return air node is treated differently to maintain backward compatibility. If
one or more Zone Return Air Node 1 Flow Rate Basis Nodes are specified in the Zone-
HVAC:EquipmentConnections object, then the mass flow rate for the first return air node
is:
∑
ṁR,1 = ReturnF lowSchedule ∗ ṁBasisN ode,j (2.80)
j

If there is only one return node and no Zone Return Air Node 1 Flow Rate Basis Nodes, then:

ṁR,1 = ReturnF lowSchedule ∗ ṁR (2.81)