The AirflowNetwork model provides the ability to simulate
the performance of an air distribution system, including
supply and return leaks, and calculate multizone airflows
driven by outdoor wind and forced air during HVAC system
operation. The pressure and airflow model described here was
developed based on AIRNET (Walton 1989). This detailed model
is used to simulate thermal conduction and air leakage losses
for constant volume air distribution systems (e.g., in
residential or light commercial buildings). The multizone
airflow calculations are performed at the HVAC system time
step which, among other benefits,.allows for modeling hybrid
ventilation systems.
The input object AirflowNetwork:SimulationControl
provides access to the airflow network method, which consists
of a set of nodes connected by airflow components through
linkages. The objects AirflowNetwork:Multizone:Zone,
AirflowNetwork:Multizone:ExternalNode, and AirflowNetwork:Distribution:Node
represent airflow nodes. The objects
AirflowNetwork:Multizone:Surface and AirflowNetwork:Distribution:Linkage
represent airflow linkages. The other objects with a
relationship between pressure and airflow represent airflow
components.
The AirflowNetwork model consists of three sequential
steps:
Pressure and airflow calculations
Node temperature and humidity calculations
Sensible and latent load calculations
The pressure and airflow calculations determine pressure at
each node and airflow through each linkage given wind
pressures and forced airflows. Based on the airflow calculated
for each linkage, the model then calculates node temperatures
and humidity ratios given zone air temperatures and zone
humidity ratios. Using these node temperatures and humidity
ratios, the sensible and latent loads from duct system
conduction and leakage are summed for each zone. The sensible
and latent loads obtained in this step are then used in the
zone energy balance equations to predict HVAC system loads and
to calculate the final zone air temperatures, humidity ratios,
and pressures.
The present AirflowNetwork model may only be applied to a
single heating and cooling system that uses a single air
distribution system (a single AirLoopHVAC
object). The model excludes the impact of the air and duct
system thermal capacitance. The impact of thermal capacity
will be addressed in future upgrades to this model.
The EnergyPlus airflow network consists of a set of nodes
linked by airflow components. Therefore, it is a simplified
airflow model, compared to detailed models such as those used
in computational fluid dynamics (CFD) models. The node
variable is pressure and the linkage variable is airflow rate.
A brief description is presented below. A detailed description
of the airflow network model may be found in the work of
Walton (1989), Dols and Walton (2002), and Walton and Dols
(2003).
Newton’s method is used to solve for node air pressures and
it requires an initial set of values for the node pressures.
There are two initialization methods available. The first is
linear initialization and equivalent to Initialization flag =
0. These initial values may be obtained by including in each
airflow component a linear approximation relating airflow to
pressure drop:
\(\dot{m}_i\) = Air mass
flow rate at i-th linkage [kg/s]
Ci = Air mass flow coefficient
[m3]
\(\Delta P_i\) = Pressure
difference across the i-th linkage [Pa]
\(\mu\) = Air viscosity
[Pa-s]
This initialization handles stack effects very well and
tends to establish the proper direction for the airflows. The
linear approximation is provided by the laminar regime.
The second initialization method assumes the initial
pressures are zero and uses Newton’s method directly.
Conservation of air mass flow rate at each linkage provides
the convergence criterion. When the sum of mass flow rates in
all the linkages approaches zero within the convergence
tolerance, the solution has converged. The solution is assumed
to have converged when the sum is less than the convergence
value, in order to reduce the number of iterations and obtain
sufficient accuracy. There are two convergence criteria used
in the AirflowNetwork model: Relative airflow convergence
tolerance and Absolute airflow convergence tolerance.
The relative airflow tolerance is equivalent to the ratio
of the absolute value of the sum of all network airflows to
the sum of the network airflow magnitudes. The absolute
airflow tolerance is the summation of the absolute value of
all network airflows. The solution has converged when both of
these convergence criteria have been met.
A linkage used in the AirflowNetwork model has two nodes,
inlet and outlet, and is linked by a component which has a
relationship between airflow and pressure. The pressure
difference across each component in a linkage is assumed to be
governed by Bernoulli’s equation:
ΔP= Total pressure difference between nodes n and m
[Pa]
Pn, Pm= Entry and exit
static pressures [Pa]
Vn, Vm= Entry and exit
airflow velocities [m/s]
*ρ= Air density [kg/m*3]
g= Acceleration due to gravity [9.81
m/s2]
zn, zm= Entry and exit
elevations [m]
By rearranging terms and adding wind pressure impacts, the
above equation may be rewritten in the format used by the
airflow network model:
\[\Delta P = {P_n} - {P_m} +
{P_S} + {P_W}\]
where
\(P_n, P_m\) = Total
pressures at nodes n and m [Pa]
PS = Pressure difference due to density
and height differences [Pa]
PW = Pressure difference due to wind
[Pa]
The Input Output Reference provides the relationship
between airflow and pressure for the most of the components
(Ref. AirflowNetwork Model). The relationship between airflow
and pressure for the
AirflowNetwork:Multizone:Component:DetailedOpening,
AirflowNetwork:Multizone:Component:SimpleOpening, and
AirflowNetwork:Multizone:Component:HorizontalOpening objects
are provided in detail in this reference.
The general problem of
gravitational flow through a vertical opening
The schematic drawing of a possible air flow pattern
through a detailed vertical opening
(AirflowNetwork:Multizone:Component:DetailedOpening) is shown
in Figure 147. The equations used below are extracted from the
COMIS Fundamentals manual (1990).
The air density is assumed to be a linear function of
height:
\[{\rho _i}(z) = {\rho _{0i}} +
{b_i}z\]
The pressure difference is assumed to be linear and
simulate the effect of turbulence:
\[\Delta {P_t} = {P_{t0}} +
{b_t}z\]
The reference pressures on each side are given at the
bottom of the opening. By assuming the Bernoulli hypothesis on
both sides, the pressure difference can be defined at any
level of z as:
The locations of the two possible neutral planes are given
by an equilibrium in pressure which leads to a zero velocity
point. By assuming the left terms in the equation above to be
zero, one may have:
This equation above can have two, one, or zero real
solutions. The zero solution represents a one-way flow through
the opening and may be expressed in the following
equation:
The discharge coefficient, opening width, opening height,
and start height factor are modulated based on opening
factors. A detailed description of opening factor calculations
may be found in the Input Output Reference (Ref.
AirflowNetwork:Multizone:Zone,
AirflowNetwork:Multizone:Surface, and
The above calculation procedure is used for a normal
rectangular window. For a horizontally pivoted rectangular
window, the calculation procedure is slightly different. A
schematic drawing of a horizontally-pivoted window is shown in
Figure 148.
Schematic drawing of a
horizontally-pivoted window
The opening angle α (0-90°) is linearly proportional to the
window opening factor (0-1.0). An opening factor of 1.0 is
equal to an opening angle of 90°. The heights in the pivoted
area are expressed as:
When z < h2 or z > h4, where z is the distance from
the bottom of the window, the integration procedure is the
same as the procedure for a normal rectangular window. When
h2<z<h4, the window width W in the above equations is
modified as:
It should be pointed out that the discharge coefficient is
modulated based on opening factors, while opening width,
opening height, and start height factor do not apply for a
horizontally-pivoted window. The actual window width and
height are used to calculate airflows for a
horizontally-pivoted window.
The schematic drawing of air flow patterns through a simple
vertical opening
(AirflowNetwork:Multizone:Component:SimpleOpening) is shown in
Figure 149. The equations used below are available from Walton
(1989).
Schematic of large opening and
associated three flow patterns
The air density for each node is assumed to be constant.
The hydrostatic equation is used to relate pressures at
various heights for each node:
\[Pn\left( y \right) = {P_{0n}}
- {\rho _n}gy\]
\[{P_m}(y) = {P_{0m}} - {\rho
_m}gy\]
where
P0n, P0m = pressure at nodes
(zones) n and m at y = 0, the reference elevation of the
opening [Pa]
ρn, ρm = air densities of
zones n and m [kg/m3]
Pn, Pm = reference pressures
of zones n and m [Pa]
It is assumed that the velocity of the airflow as a
function of height is given by the orifice equation (Brown and
Solvason 1962):
When the neutral plane is within the opening (first pattern
in Figure 149), two-way (bi-directional) flows occur. The
total flow through a large opening is the sum of both
flows.
The opening width is modulated based on opening factors. A
detailed description of opening factor calculations may be
found in the Input Output Reference (Ref.
AirflowNetwork:Multizone:Zone,
AirflowNetwork:Multizone:Surface, and
AirflowNetwork:Multizone:Component:DetailedOpening).
The above two objects are used to simulate airflows across
large vertical openings. The simple opening component
(AirflowNetwork:Multizone:Component:SimpleOpening) assumes the
pressure difference across the opening is a function of height
varied from opening bottom to top, so that two-way flow may be
obtained if appropriate (Walton 1989). The Detailed Opening
component (AirflowNetwork:Multizone:Component:DetailedOpening)
assumes both the pressure difference across the opening and
air density are a function of height, so that three-way flow
may be obtained (COMIS 1990). If these opening models would be
used for horizontal openings, the pressure difference across
the opening and air density remain constant, so that only
one-way flow is possible using the detailed and simple opening
components which are meant for vertical or near-vertical
openings. In reality, there are two-way flows (air
recirculation) across a large horizontal opening caused by
buoyancy due to temperature and pressure difference and forced
flow driven by air pressure difference only. Therefore, a
horizontal opening component
(AirflowNetwork:Multizone:Component:HorizontalOpening) is
available to simulate airflows across large horizontal
openings with the possibility of two-way flow by combining
forced and buoyancy airflows together.
The model for horizontal openings consists of forced
airflow, buoyancy airflow, purge pressure and sloping plane.
The model is mainly from a NIST report presented by Cooper
(1989). The sloping plane (Bolmqvist and Sandberg 2004)
portion of the model was added to allow for staircase
simulations.
For simplicity, a two zone building (upper and lower zones)
connected by a large horizontal opening is used to describe
the model, as shown in Figure 150. Forced and buoyancy
airflows are described separately below.
Air movements across a
horizontal opening
Forced airflows
The air mass flow rate is determined by the pressure
difference across the opening. The relationship between
pressure and airflow is the same as AIRNET for a component
(see AirflowNetwork:Multizone:Component:SimpleOpening
description above). Since the height of the opening is
constant, the forced airflow is unidirectional. A positive
value for pressure difference indicates flow direction is from
the lower zone to the upper zone across the opening, while a
negative value represents flow in the opposite direction. The
following description addresses forced air mass flow rates and
partial derivatives for three possible scenarios of pressure
difference:
Buoyancy flow only occurs when the air density in the upper
zone is greater than the air density in the lower zone. The
flow is bi-directional and the amount of upper flow is equal
to the lower flow across the opening. The following discussion
assumes the air density in the upper zone is greater than the
air density in the lower zone. Otherwise, the buoyancy flow
rate is equal to zero. It is also assumed that the maximum
buoyancy flow occurs when the pressure difference across the
opening due to forced airflows is zero. The maximum buoyancy
flow may be expressed as a part of Cooper’s model:
When forced and buoyancy flows co-exist, it is possible to
have either unidirectional or bi-directional flows. For
example, when the upward force due to pressure difference is
greater than the buoyancy force (downward), unidirectional
flow occurs. Bi-directional flow only occurs when the upward
imposed force is less than the buoyancy force. The critical
pressure between unidirectional and bi-directional flows is
called the purge pressure (Tan and Jaluria 1992). The purge
pressure is a function of opening geometry and the buoyancy
force (ΔP/(gΔρDH)) and may be expressed as (Cooper
1998):
As mentioned above, when the air pressure difference
between two zones is zero there is the maximum bi-directional
flow due to the buoyancy force. When the pressure difference
increases from 0 and is less than |ΔPFlood|, there
is some bi-directional flow across the opening, but less than
the maximum flow. If the pressure difference keeps increasing
and exceeds |ΔPFlood|, there is no bi-directional
flow. Cooper’s model assumes the buoyancy flow varies linearly
with pressure difference.
The total air flow across the opening is based on
superposition of the forced and buoyancy flows, and may be
expressed for three different pressure difference scenarios as
follows:
A Staircase is attached to the
horizontal opening.
When a staircase is introduced as shown in Figure 151, the
effective opening area will be used to replace A (opening
area) in the above equations. The effective area may be
estimated as (Bolmqvist and Sandberg 2004):
\[{A_{eff}} = A\sin \alpha (1 +
\cos \alpha )\]
where:
Aeff= Effective area of horizontal opening
[m2]
A= Area of horizontal opening [m2]
α= Angle between the stair plane and horizontal opening
plane [degrees]
Note: the hydraulic diameter calculation is based on the
effective opening area, while the opening depth remains the
same.
Figure 152 demonstrates possible forced and buoyancy flow
rates at different ratios of pressure difference to purging
pressure across a horizontal opening when the upper zone air
density is greater than the lower zone air density. The
pressure difference is the lower zone pressure minus the upper
zone pressure. Otherwise, the buoyancy flow is zero. In
addition, when the absolute ratio is above 1, the buoyancy
flow is also zero. The following table provides a brief
description for the legend listed in Figure 152.
Legend Description
Legend
Description
Forced downward
Forced flow rate from upper to
lower at PL-PU < 0
Forced upward
Forced flow rate from lower to
upper at PL-PU > 0
Buoyancy upward
Total upward flow rate due to
buoyancy only at PL-PU < 0
Buoyancy downward
Total downward flow rate due to
buoyancy only at PL-PU > 0
Combined downward
Total downward flow at
PL-PU < 0 (Forced downward + buoyancy
upward)
Combined upward
Total upward flow at
PL-PU > 0 (Forced upward + buoyancy
downward)
Cp is a function of location on the
building envelope and wind direction. When Wind Pressure
Coefficient Type = “INPUT”, the Cp values
are explicitly defined in the input for
AirflowNetwork:Multizone:Wind Pressure Coefficient Values.
When Wind Pressure Coefficient Type = “AVERAGE-SURFACE
CALCULATION” and the building shape is rectangular, the
program uses the following equations to calculate wind
pressure coefficient (Cp) values for
different wind directions. For a low rise building, the
normalized surface pressure coefficient may be written as
(Swami and Chandra 1988):
Cp,n = Cp value at
a given angle between wind direction and the outward normal of
the surface under consideration [dimensionless]
α = Angle between wind direction and outward normal of wall
under consideration [deg]
G = Natural log of the ratio of the width of the
wall under consideration to the width of the adjacent wall
[dimensionless]
n = Index of incident angle at 30-degree
increments
For walls of a high rise building, a two-dimensional array
of surface-averaged wind pressure coefficients is generated
based on wind incident angle and side ratio. The wind pressure
coefficients are provided in 2001 ASHRAE Fundamentals
Handbook, p. 16.5, Fig. 7, “Surface Averaged Wall Pressure
Coefficients for Tall Buildings”. The original work was
performed by Atkins et al. (1979). The incident angle has an
increment of 30 degrees. The side ratio values are 0.25, 1.0,
and 4.0. For a given incident angle and building aspect ratio,
the program uses linear interpolation to calculate the
corresponding wind pressure coefficient
Cp,n.
For the roof of a high rise building, a two-dimensional
array of surface-averaged wind pressure coefficients is also
generated based on wind incident angle and side ratio. The
wind pressure coefficients are provided in 2001 ASHRAE
Fundamentals Handbook, p. 16.6, Fig. 9, “Surface Averaged Roof
Pressure Coefficients for Tall Buildings”. The original work
was performed by Holmes (1986). The incident angle has an
increment of 30 degrees. The side ratio values are 0.25, 1.0,
and 4.0. At a given wind incident angle and building aspect
ratio, the program uses linear interpolation to calculate the
corresponding wind pressure coefficient
Cp,n.
The wind surface pressure at the given incident angle can
be calculated by combining the above two equations:
Based on the relationship between airflow rate and pressure
drop for each component, a system of equations for all
components can be assembled together in an n x n square
matrix, where n is the number of nodes. Newton’s method is
used to iteratively solve for the air pressure at each node. A
new estimated vector for all node pressures, {P}*, is computed
from the current estimated vector of node pressures, {P},
by:
\[{\left\{ P \right\}^*} =
\left\{ P \right\} - \left\{ C \right\}\]
where the correction vector, {C}, is computed by the matrix
relationship:
\[\left[ J \right]\left\{ C
\right\} = \left\{ B \right\}\]
{B} is a column vector with each component given by:
where n is the node number and i indicates all flow paths
connecting node n to other nodes, and [J] is the square
Jacobian matrix whose elements are given by:
The convergence tolerance is used to check the sum of the
mass flow rates by applying mass conservation. The convergence
acceleration equation shown below is used to correct the node
pressures to more rapidly obtain a solution. By assuming a
constant ratio of correction values from one iteration to the
next, the following method is applied:
\[{P_n}^* = {P_n} - {C_n}/(1 -
r)\]
where
r= the ratio of Cn for the current
iteration to its value for the previous iteration
[dimensionless]
Cn= Correction value at the
nth node [Pa]
Pn= Estimated pressure at the
nth node [Pa]
Pn*= Corrected pressure at
the nth node used in the next iteration [Pa]
This method is similar to a Steffensen iteration (Conte and
de Boor 1972) which is used as a fixed-point iteration method
for individual nonlinear equations.
The iteration correction method presented in the above
equation gives a variable factor. When the solution is close
to convergence, the solution method converges quadratically.
By limiting cases where the value of r is less than some
value, such as -0.5, the solution will not interfere with the
rapid convergence. It has not been proven that the convergence
acceleration equation will always lead to convergence, but it
can be shown that it will not prevent convergence. Newton’s
method converges when the estimated solution values are within
some distance, called the radius of convergence, or the
correct solution. Applying the convergence acceleration
equation when -1 < r < 0, will cause a smaller
correction than Newton’s method, which therefore, can not
force the iterations outside the radius of convergence. When r
< -1, the solution diverges in an oscillatory fashion. When
r > 1, the solution also diverges, but in a nonoscillatory
manner. For 0 < r < 1, the solution is approached from
one direction. In all three cases, the convergence
acceleration equation applies as long as r is truly constant
over several iterations. However, for the last case, this
involves a true extrapolation of correction factor which is
very sensitive to the accuracy of r. This is most extreme for
the case of r=1, which would cause an infinite correction.
A brief description of the air node temperature calculation
is given below. A detailed description can be found in the
work of Swami et al. (1992). The following equation is used to
calculate temperature distribution across a duct element at
the given airflow rate and inlet air temperature:
The outlet air temperature can be calculated using the
above equation at the given inlet air temperature. Since the
inlet temperature at one linkage is the outlet temperature for
the connected linkage, the outlet air temperatures at all
nodes are solved simultaneously. A square linear system
assembled by the AirflowNetwork model is expressed below:
\[\{ M\} [T] = [B]\]
where
{M}= Airflow matrix
[T]= Temperature vector
[B]= Given boundary conditions
The zone air temperatures and primary air loop component
(fan and coils) outlet conditions are used as prescribed
conditions in the AirflowNetwork model. In addition, the
temperature difference across zone loop components (terminal
units) is held constant during the calculations. For example,
thermal zone temperatures calculated during the previous
system time step are used as prescribed temperatures when
calculating all other node temperatures. The zone air
temperature is assumed constant (prescribed) throughout the
AirflowNetwork iterative solution. The fan and coil outlet air
temperatures, and terminal unit temperature differences are
assumed constant within an AirflowNetwork iteration. The
sensible heat gains calculated during the AirflowNetwork
solution are then used to predict a new zone air
temperature.
A brief description of the air node humidity ratio
calculation is given below. A detailed description can found
in the work of Swami et al. (1992). The following equation is
used to calculate humidity ratio distribution across a duct
element at the given airflow rate and inlet air humidity
ratio:
The outlet air humidity ratio can be calculated using the
above equation at the given inlet air humidity ratio. Since
the inlet humidity ratio at one linkage is the outlet humidity
ratio for the connected linkage, the outlet air humidity ratio
at all nodes are solved simultaneously. A square linear system
assembled by the AirflowNetwork model is expressed below:
\[\{ {M_m}\} [W] =
[{B_m}]\]
where
{Mm}= Airflow matrix
[W]= Humidity ratio vector
[Bm]= Given boundary conditions
The zone air humidity ratios and primary air loop component
(fan and coils) outlet conditions are used as prescribed
conditions in the AirflowNetwork model. For example, thermal
zone humidity ratios calculated during the previous system
time step are used as prescribed humidity ratios when
calculating all other node humidity ratios. The zone air
humidity ratio is assumed constant (prescribed) throughout the
AirflowNetwork iterative solution. The coil outlet air
humidity ratio is assumed constant within an AirflowNetwork
iteration. The latent heat gains calculated during the
AirflowNetwork solution are then used to predict a new zone
air humidity ratio.
The zone sensible and latent loads calculated in the
AirflowNetwork model consist of multizone, duct conduction and
leakage. In addition, the impact of infiltration and mixing is
accounted for in this calculation. The multizone load only
includes incoming airflows from outside (infiltration) and
other adjacent zones (mixing) with and without forced-air fan
operation. It is divided into two terms: variable and
constant. The constant term is the sum of the mass flow rate
multiplied by the specific heat for both infiltration and
mixing. The variable term includes the impact of zone and
outdoor air temperature. Each of these terms is used in the
zone energy balance equation. The sensible load items from the
multizone load calculations may be written as follows:
where
MCPairflow = Sum of air mass flow rate
multiplied by specific heat for infiltration and mixing
[W/K]
MCPTairflow = Sum of air mass flow rate
multiplied by specific heat and temperature for infiltration
and mixing [W]
= Incoming air mass flow rate from outdoors [kg/s]
= Incoming air mass flow rate from adjacent zones
[kg/s]
Tamb= Outdoor air dry-bulb temperature
[°C]
Tzone= Adjacent zone air temperature
[°C]
The latent load items from multizone load calculations may
be written as follows:
where
Mairflow = Sum of air mass flow rates
for infiltration and mixing [kg/s]
MWairflow = Sum of air mass flow rate
multiplied by humidity ratio for infiltration and mixing
[kg/s]
= Incoming air mass flow rate from outdoors [kg/s]
= Incoming air mass flow rate from adjacent zones
[kg/s]
Wamb= Outdoor air humidity ratio
[kg/kg]
Wzone= Adjacent zone air humidity ratio
[kg/kg]
The air distribution system (ADS) loads due to duct
conduction and leakage depend on the air distribution system
component (e.g., duct) location. The air distribution system
sensible and latent loads are calculated for each zone as
follows:
QADS,i= Total sensible load in the i-th
zone due to ADS losses [W]
Qcond(ij)= Duct
wall conduction loss at the j-th duct located in the i-th zone
[W]
Qleak(ij)= Sensible supply leak loss at
the j-th linkage located in the i-th zone [W]
QADS,m,i= Total latent load in the i-th
zone due to ADS losses [kg/s]
Qcond,m(ij)= Duct
wall vapor diffusion loss at the j-th duct located in the i-th
zone [kg/s]
Qleak,m(ij) = Latent supply leak loss
at the j-th linkage located in the i-th zone [kg/s]
Impacts
of Supply Air Constant Volume Fan Control on Load: Cycling
vs. Continuous[LINK]
The AirflowNetwork model currently allows two types of
constant volume fans: Fan:ConstantVolume
and Fan:OnOff.
The Fan:ConstantVolume
object has only one type of supply air fan operation mode:
continuous fan, cycling compressor
(ContinuousFanWithCyclingCompressor). However, the Fan:OnOff
has two types of supply air fan operation modes: cycling fan,
cycling compressor (CyclingFanAndCompressor) or continuous
fan, cycling compressor (ContinuousFanWithCyclingCompressor).
The CyclingFanAndCompressor operation mode is frequently
referred to as “AUTO fan”, where the compressor(s) and supply
air fan operate in unison to meet the zone heating or cooling
load, and cycle off together when the heating or cooling load
has been met. The ContinuousFanWithCyclingCompressor operation
mode is often referred to as “fan ON”, where the compressor(s)
cycle on and off to meet the zone heating or cooling load but
the supply air fan operates continuously regardless of
compressor operation. The supply air fan operation mode is
specified in an HVAC system object based on a given fan
operation mode schedule (e.g., AirLoopHVAC:UnitaryHeatCool
object).
The determination of the zone sensible and latent loads
caused by multizone airflows and forced air distribution
system operation is dependent on the supply air fan operation
mode (see Sensible and Latent Load Calculations section
above). The zone loads calculated by the AirflowNetwork model
are added to zone sensible and latent balances in the
ZonePredictorCorrector module to calculate zone air
temperatures and humidity ratios (see Integration of the
AirflowNetwork Model section below). For the case of
continuous fan/cycling compressor, the zone loads during
forced air distribution system operation are calculated with
the system design air mass flow rate without modification,
since the system air node conditions (temperature and
humidity) reflect the average values for the time step
considering the coil/fan on and off periods during the time
step.
For the case of cycling fan/cycling compressor, where the
forced air distribution system can operate for a portion of
the simulation time step, the airflows are determined based on
the air distribution system part-load ratio (ratio of the
average air mass flow rate for the time step divided by the
design air mass flow rate). The airflows for
the AirflowNetwork:Distribution:Linkage
objects are reported during the air distribution system on
cycle, since no airflow is assumed during the system off
cycle. The airflows for the
AirflowNetwork:Multizone:Surface objects are weighted by the
air distribution system part-load ratio. The zone
loads are the sum of energy losses during both the
air distribution system on and off cycle at each system time
step. The energy losses when the air distribution system is on
are calculated using the system “on” air flow rate multiplied
by the air distribution system run time fraction. The energy
losses when the air distribution system is off are obtained
from the multizone airflow calculations (without forced air
through the air distribution system) and these losses are
multiplied by (1.0 - system run time fraction), assuming no
airflows through the air distribution system when the fan is
off. The formulas used to calculate airflows and zone loads
are given below:
Airflow
Airflow = Airflow during ADS on cycle * ADS Part-load ratio
+ Airflow during ADS off cycle * (1.0 - ADS Part-load
ratio)
System run time fraction = Maximum run time fraction of
coils and fans in the air distribution system
Zone
energy losses = Zone
energy loss during ADS on cycle * System run time fraction +
Zone
energy loss during ADS off cycle * (1.0 - System run time
fraction)
The calculation of loads due to multizone airflow, without
forced air distribution system operation, is performed when
the HVAC system is first simulated during a simulation time
step (FirstHVACIteration = True). The calculation of loads due
to air distribution system operation is performed on
subsequent iterations within the same simulation time step
when the mass flow rate at the supply air fan inlet is greater
than 0.0 (and FirstHVACIteration = False).
Airflow
Calculation Procedure using A Supply Variable Air Volume
Fan[LINK]
The AirflowNetwork model currently also allows a variable
air volume fan type as Fan:VariableVolume.
The allowed terminal type is AirTerminal:SingleDuct:VAV:Reheat
only. Other types of terminals will be added later.
In general, the supply fan air flow rate in a VAV central
system is determined by a sum of terminal flow rates when the
AirflowNetwork model is not applied. When the AirflowNetwork
model is applied and the supply air fan flow rate is given,
each terminal flow is determined by pressure resistance of
each supply air pathway. It is possible that the delivered air
flow rate from the pressure resistance at each terminal may be
totally different from the desired flow rate determined by
terminal units. Therefore, it is not easy to meet both
requirements. The following compromised approach, including
possible supply and return leaks in an air distribution
system, is implemented.
Set up terminal airflows in the AirflowNetwork module based
on the SimVAV subroutine in the HVACSingleDuctSystem
module.
Require AirflowNetwork:Distribution:Component:LeakageRatio
objects to define supply leaks, so that the values of the
Effective Leakage Ratio field are used to decide the supply
fan flow rates. The base of the ratio will be actual supply
fan flow rate, instead of the maximum fan flow rate used in
the constant volume fan systems.
Assign the supply fan airflow rate based on the sum of all
terminal flow rates and all supply leak ratios until it
reaches the maximum fan flow rate
If the calculated supply fan flow rate is above the maximum
limit of the supply fan flow rate, and the supply leak ratios
remain the same, the supply fan flow rate is set to the
maximum limit, and the terminal flow rates are reduced
proportionally weighted by a ratio of the maximum supply fan
flow rate by input to the calculated supply fan flow rate.
Therefore, a sum of all terminal rates and all supply leak
rates is equal to the maximum supply fan rate.
The loads calculated by the AirflowNetwork model are
integrated into the EnergyPlus heat balance equation in a
similar manner as described elsewhere in this document in the
section “Basis for the Zone
and System Integration”. The mass flow rate summations and
sensible and latent loads described in the previous section
are included in the calculation of zone temperature and
humidity ratio.
The revised zone temperature update equation becomes:
Where MCPTairflow is the sum of mass flow rate
multiplied by specific heat and temperature for infiltration
and mixing, QADS,z is the added total sensible load
in the zone due to Air Distribution System losses, and
MCPairflow is the sum of mass flow rate multiplied
by specific heat for infiltration and mixing as calculated in
the AirflowNetwork model described above.
The revised coefficient (B) used in the zone humidity ratio
calculation is shown below:
Where MWairflow is the sum of mass flow rate
multiplied by humidity ratio for infiltration and mixing and
QADS,m,zis theadded total latent
(moisture) load in the zone due to Air Distribution System
losses from the AirflowNetwork model described above. This
coefficient is used in the prediction of moisture as described
in the section “Moisture Predictor-Corrector” found elsewhere
in this document.
The AirflowNetwork:OccupantVentilationControl
object enhances the AirflowNetwork model and provides more
practical and advanced controls for window opening and closing
operations, based on Marais & Teichmann (2014). This
control includes minimum open and closed time control, indoor
thermal comfort control, and opening and closing probability
controls.
It should be pointed that the open elapsed time and closed
elapsed time are not independent. In other words, when one of
the elapsed times value is greater than 0, then other elapsed
time value must be equal to 0. The model is either tracking a
window as open, in which case the open elapsed time grows, or
as closed, in which case the closed elapsed time grows.
The calculation procedures are presented in the following
steps:
Step 1: Open elapsed time check
The model checks the open elapsed time first. When the open
elapsed time is greater than 0 and less than the minimum
opening time, a window will remain open. Otherwise, the model
goes to Step 2.
Open elapsed time > minimum open time
Step 2: Closed elapsed time check
This step checks the closed elapsed time. When the closed
elapsed time is greater than 0 and less than the minimum
closed time, a window will remain closed. Otherwise, the model
goes to Step 3.
Closed elapsed time > minimum closed time
Step 3: Elapsed time only?
When either the open elapsed time or the closed elapsed
time is long enough (greater than the minimum time) and no
other checks are needed, the model returns to the other
ventilation control defined in the Ventilation Control Mode
field in AirflowNetwork:MultiZone:Zone or
AirflowNetwork:MultiZone:Surface. If other checks are needed,
the following steps will be performed.
Step 4: Thermal comfort temperature calculation
The thermal comfort check requires the thermal comfort
temperature and the comfort band. The comfort temperature of
Tcomf is calculated as a function of the outdoor dry-bulb
temperature Tout. The comfort temperature calculation may be
based on two curves and a boundary temperature point.
After calculation of the comfort temperature and band, the
thermal comfort check will be performed using the zone air
operative temperature Tg. The check consists of upper and
lower boundary checks.
The upper boundary check checks window opening status and
uses the following logic:
\[T_g > \left(T_{comf} +
\theta\right)\]
If the above logic check is true, the opening probability
check will be performed. The detailed description is given in
the opening probability section.
If the above logic check is false, no action is needed.
If the opening probability check is true and upper boundary
check is satisfied, a window will be opened, regardless of
open/closed status at the previous time step. If the opening
probability check is false and upper boundary check is
satisfied, a window will remain at the status from the
previous time step.
The lower boundary check will follow the upper boundary
check to check the window closing status, using the following
logic.
\[T_g < \left(T_{comf} -
\theta\right)\]
If the above logic check is true, the closing probability
check will be performed. The detailed description is given in
the closing probability section.
If the above logic check is false, no action is needed.
If the closing probability check is true and lower boundary
check is satisfied, a window will be closed, regardless of
open/closed status at the previous time step. If the closing
probability check is false and lower boundary check is
satisfied, a window will remain at the status from the
previous time step.
The output variables from the model are open status,
opening probability status, and closing probability status.
The detailed description of opening status is given in the
Airflow Network Outputs section in the Input Output
Reference.
Note: The upper and lower boundary checks are not
independent. In other words, when one of boundary check is
satisfied, the other check will be dissatisfied.
Opening probability control provides an optional random
number check. The control logic of opening probability is
described in the following steps.
Step 1: Closed elapsed time check
This check requires that closed elapsed time is longer than
the minimum closing time.
Closed time > minimum closed time
If the time duration is not long enough, the output is
false so that a window remains closed.
If the time duration is long enough, an occupancy check is
performed.
Step 2: Occupancy check
If a zone is not occupied and the occupancy check is
requested, the output is false. If the zone is occupied, the
next step is to check the setpoints using zone air temperature
at the previous time step as a reference.
Step 3: Setpoint check
There are 5 temperature control types. The following types
are available:
No control: Bypass
Single heating setpoint: If Tzon > setpoint, go to
next step. Otherwise, return false.
Single cooling setpoint: If Tzon < setpoint, go to
next step. Otherwise, return false.
Single heating and cooling setpoint: no action by
returning false
Dual heating and setpoints: If heating setpoint <
Tzon < Cooling setpoint, go to next step. Otherwise, return
false.
Step 4: Select bypass or opening probability check
A choice is provided at this stage so that the opening
probability check may be performed or bypassed. If bypassed,
the output will be true to open a window. If performed,
probability will be determined from a schedule.
Step 5: Perform opening probability check
The opening probability (OP) value is determined from a
schedule.
OP = Schedule value or specific function
If the probability value is greater than a random number,
the output is true.
OP > random number (random number is between 0 and
1)
The control logic of closing probability is described as
follows.
Step 1: Open elapsed time check
This check requires that open elapsed time is longer than
the minimum opening time.
Open time > minimum open time
If the time duration is not long enough, the output is
false and the window remains open. Otherwise, Step 2 will be
performed.
Step 2: Select by-pass or closing probability check
A choice is provided at this step so that the closing
probability check may be performed or bypassed. If bypassed,
the output will be true to close a window. If performed, the
closing probability will be calculated. The closing
probability (CP) value is given from a schedule.
CP = Schedule value or specific function
If the closing probability check is performed and the
closing probability is greater than a random number, the
output is true.
CP > random number (random number is between 0 and
1)
Atkins, R. E., J. A. Peterka, and J. E. Cermak. 1979.
“Averaged pressure coefficients for rectangular buildings,”
Wind Engineering, Proceedings of the Fifth International
Conference 7:369-80, Fort Collins, CO. Pergamon Press, NY.
Bolmqvist, C. and M. Sandberg, 2004, “Air Movements through
Horizontal Openings in Buildings - A Model Study,”
International Journal of Ventilation, Vol. 3, No. 1,
pp. 1-9
COMIS Fundamentals. 1990. Edited by Helmut E. Feustel and
Alison Rayner-Hooson, LBL-28560, Lawrence Berkeley Laboratory,
Berkeley, CA
Conte, S. D. and C de Boor. 1972. Elementary Numerical
Analysis: an Algorithmic Approach, McGraw-Hill.
Cooper, L., 1989, “Calculation of the Flow Through a
Horizontal Ceiling/Floor Vent,” NISTIR 89-4052, National
Institute of Standards and Technology, Gaithersburg, MD
Dols, W. S. & G. N. Walton. 2002. “CONTAMW 2.0 User
Manual,” NISTIR 6921, National Institute of Standards and
Technology, Gaithersburg, Maryland
Holmes, J. D. 1986. Wind Loads on low-rise buildings: The
structural and environmental effects of wind on buildings and
structures, Chapter 12, Faculty of Engineering, Monash
University, Melbourne, Australia
Marais, J. M. & C. Teichmann, “Window Simulation
Methods Required for Manual Window
Ventilated Buildings,” Fifth German-Austrian IBPSA Conference,
September 22-24, 2014, RWTH Aachen University, Germany
Swami, M. V. and S. Chandra. 1988. Correlations for
pressure distribution on buildings and calculation of
natural-ventilation airflow, ASHRAE Transactions 94(1988) (Pt
1), pp. 243-266.
Swami, M. V., L. Gu, & V. Vasanth. 1992. “Integration
of Radon and Energy Models for Building,”
FSEC-CR-553-92, Florida Solar Energy Center, Cocoa,
Florida
Tan, Q. and Y. Jaluria, 1992, “Flow through Horizontal
Vents as Related to Compartment Fire Environments,”
NIST-GCR-92-607, National Institute of Standards and
Technology, Gaithersburg, Maryland
Walton, G. N. 1989. “AIRNET - A Computer Program for Building
Airflow Network Modeling,” NISTIR 89-4072, National Institute
of Standards and Technology, Gaithersburg, Maryland
Walton, G. N. & W. S. Dols. 2003. “CONTAM 2.1
Supplemental User Guide and Program Documentation,” NISTIR
7049, National Institute of Standards and Technology,
Gaithersburg, Maryland
AirflowNetwork Model[LINK]
Overview[LINK]
The AirflowNetwork model provides the ability to simulate the performance of an air distribution system, including supply and return leaks, and calculate multizone airflows driven by outdoor wind and forced air during HVAC system operation. The pressure and airflow model described here was developed based on AIRNET (Walton 1989). This detailed model is used to simulate thermal conduction and air leakage losses for constant volume air distribution systems (e.g., in residential or light commercial buildings). The multizone airflow calculations are performed at the HVAC system time step which, among other benefits,.allows for modeling hybrid ventilation systems.
Model Description[LINK]
The input object AirflowNetwork:SimulationControl provides access to the airflow network method, which consists of a set of nodes connected by airflow components through linkages. The objects AirflowNetwork:Multizone:Zone, AirflowNetwork:Multizone:ExternalNode, and AirflowNetwork:Distribution:Node represent airflow nodes. The objects AirflowNetwork:Multizone:Surface and AirflowNetwork:Distribution:Linkage represent airflow linkages. The other objects with a relationship between pressure and airflow represent airflow components.
The AirflowNetwork model consists of three sequential steps:
Pressure and airflow calculations
Node temperature and humidity calculations
Sensible and latent load calculations
The pressure and airflow calculations determine pressure at each node and airflow through each linkage given wind pressures and forced airflows. Based on the airflow calculated for each linkage, the model then calculates node temperatures and humidity ratios given zone air temperatures and zone humidity ratios. Using these node temperatures and humidity ratios, the sensible and latent loads from duct system conduction and leakage are summed for each zone. The sensible and latent loads obtained in this step are then used in the zone energy balance equations to predict HVAC system loads and to calculate the final zone air temperatures, humidity ratios, and pressures.
The present AirflowNetwork model may only be applied to a single heating and cooling system that uses a single air distribution system (a single AirLoopHVAC object). The model excludes the impact of the air and duct system thermal capacitance. The impact of thermal capacity will be addressed in future upgrades to this model.
Pressure and Airflow Calculations[LINK]
The EnergyPlus airflow network consists of a set of nodes linked by airflow components. Therefore, it is a simplified airflow model, compared to detailed models such as those used in computational fluid dynamics (CFD) models. The node variable is pressure and the linkage variable is airflow rate. A brief description is presented below. A detailed description of the airflow network model may be found in the work of Walton (1989), Dols and Walton (2002), and Walton and Dols (2003).
Initialization[LINK]
Newton’s method is used to solve for node air pressures and it requires an initial set of values for the node pressures. There are two initialization methods available. The first is linear initialization and equivalent to Initialization flag = 0. These initial values may be obtained by including in each airflow component a linear approximation relating airflow to pressure drop:
\[\dot{m}_i = C_i \rho \left( \frac{\delta P_i}{\mu} \right)\]
where
\(\dot{m}_i\) = Air mass flow rate at i-th linkage [kg/s]
Ci = Air mass flow coefficient [m3]
\(\Delta P_i\) = Pressure difference across the i-th linkage [Pa]
\(\mu\) = Air viscosity [Pa-s]
This initialization handles stack effects very well and tends to establish the proper direction for the airflows. The linear approximation is provided by the laminar regime.
The second initialization method assumes the initial pressures are zero and uses Newton’s method directly.
Convergence criteria[LINK]
Conservation of air mass flow rate at each linkage provides the convergence criterion. When the sum of mass flow rates in all the linkages approaches zero within the convergence tolerance, the solution has converged. The solution is assumed to have converged when the sum is less than the convergence value, in order to reduce the number of iterations and obtain sufficient accuracy. There are two convergence criteria used in the AirflowNetwork model: Relative airflow convergence tolerance and Absolute airflow convergence tolerance.
Relative airflow tolerance = \(\frac{{\,\,\left| {\,\sum\limits_{} {{{\mathop m\limits^ \bullet }_{_i}}} } \right|\,\,}}{{\sum\limits_{} {\left| {{{\mathop m\limits^ \bullet }_{_i}}} \right|} }}\)
Absolute airflow tolerance = \(\left| {\sum {{{\mathop m\limits^ \bullet }_{_i}}} } \right|\)
The relative airflow tolerance is equivalent to the ratio of the absolute value of the sum of all network airflows to the sum of the network airflow magnitudes. The absolute airflow tolerance is the summation of the absolute value of all network airflows. The solution has converged when both of these convergence criteria have been met.
Linkage models[LINK]
A linkage used in the AirflowNetwork model has two nodes, inlet and outlet, and is linked by a component which has a relationship between airflow and pressure. The pressure difference across each component in a linkage is assumed to be governed by Bernoulli’s equation:
\[\Delta P = \left( {Pn + \frac{{\rho V_n^2}}{2}} \right) - \left( {Pm + \frac{{\rho V_m^2}}{2}} \right) + \rho g\left( {zn - zm} \right)\]
where
ΔP= Total pressure difference between nodes n and m [Pa]
Pn, Pm= Entry and exit static pressures [Pa]
Vn, Vm= Entry and exit airflow velocities [m/s]
*ρ= Air density [kg/m*3]
g= Acceleration due to gravity [9.81 m/s2]
zn, zm= Entry and exit elevations [m]
By rearranging terms and adding wind pressure impacts, the above equation may be rewritten in the format used by the airflow network model:
\[\Delta P = {P_n} - {P_m} + {P_S} + {P_W}\]
where
\(P_n, P_m\) = Total pressures at nodes n and m [Pa]
PS = Pressure difference due to density and height differences [Pa]
PW = Pressure difference due to wind [Pa]
The Input Output Reference provides the relationship between airflow and pressure for the most of the components (Ref. AirflowNetwork Model). The relationship between airflow and pressure for the AirflowNetwork:Multizone:Component:DetailedOpening, AirflowNetwork:Multizone:Component:SimpleOpening, and AirflowNetwork:Multizone:Component:HorizontalOpening objects are provided in detail in this reference.
The schematic drawing of a possible air flow pattern through a detailed vertical opening (AirflowNetwork:Multizone:Component:DetailedOpening) is shown in Figure 147. The equations used below are extracted from the COMIS Fundamentals manual (1990).
The air density is assumed to be a linear function of height:
\[{\rho _i}(z) = {\rho _{0i}} + {b_i}z\]
The pressure difference is assumed to be linear and simulate the effect of turbulence:
\[\Delta {P_t} = {P_{t0}} + {b_t}z\]
The reference pressures on each side are given at the bottom of the opening. By assuming the Bernoulli hypothesis on both sides, the pressure difference can be defined at any level of z as:
\[{P_1}(z) - {P_2}(z) = ({P_{01}} - {P_{02}}) - g\left[ {({\rho _{01}}z + {b_1}{z^2}/2) - ({\rho _{02}}z + {b_2}{z^2}/2)} \right] + ({P_{t0}} + {b_t}z)\]
The velocity at any level z is given by
\[v(z) = \sqrt {2\frac{{{P_1}(z) - {P_2}(z)}}{\rho }} \]
The locations of the two possible neutral planes are given by an equilibrium in pressure which leads to a zero velocity point. By assuming the left terms in the equation above to be zero, one may have:
\[g({b_1} - {b_2}){z^2}/2 + \left[ {g({\rho _{01}} - {\rho _{02}}) - {b_t}} \right]z + ( - {P_{01}} + {P_{02}} - {P_{t0}}) = 0\]
This equation above can have two, one, or zero real solutions. The zero solution represents a one-way flow through the opening and may be expressed in the following equation:
\[\dot m = {C_d}\theta \int\limits_{z = 0}^{z = H} {\rho \nu (z)Wdz} \]
The one real solution represents two-way (bi-directional) flow, which may be written in the following equations.
\[{\dot m_{0,z1}} = {C_d}\theta \int\limits_{z = 0}^{z = z1} {\rho \nu (z)Wdz} \]
\[{\dot m_{z1,H}} = {C_d}\theta \int\limits_{z = z1}^{z = H} {\rho \nu (z)Wdz} \]
The two real solutions represent three-way flow, which may be written in the following equations.
\[{\dot m_{z2,H}} = {C_d}\theta \int\limits_{z = z2}^{z = H} {\rho \nu (z)Wdz} \]
\[{\dot m_{z1,z2}} = {C_d}\theta \int\limits_{z = z1}^{z = z2} {\rho \nu (z)Wdz} \]
\[{\dot m_{z2,H}} = {C_d}\theta \int\limits_{z = z2}^{z = H} {\rho \nu (z)Wdz} \]
where
Cd= discharge coefficient [dimensionless]
θ= Area reduction factor [dimensionless]
W= Opening width [m]
The discharge coefficient, opening width, opening height, and start height factor are modulated based on opening factors. A detailed description of opening factor calculations may be found in the Input Output Reference (Ref. AirflowNetwork:Multizone:Zone, AirflowNetwork:Multizone:Surface, and
AirflowNetwork:Multizone:Component:DetailedOpening).
The above calculation procedure is used for a normal rectangular window. For a horizontally pivoted rectangular window, the calculation procedure is slightly different. A schematic drawing of a horizontally-pivoted window is shown in Figure 148.
The opening angle α (0-90°) is linearly proportional to the window opening factor (0-1.0). An opening factor of 1.0 is equal to an opening angle of 90°. The heights in the pivoted area are expressed as:
\[h2 = AxisHeight(1 - \cos (\alpha ))\]
\[h4 = AxisHeight + (WindowHeight - AxisHeight)\cos (\alpha )\]
When z < h2 or z > h4, where z is the distance from the bottom of the window, the integration procedure is the same as the procedure for a normal rectangular window. When h2<z<h4, the window width W in the above equations is modified as:
\[{W_{pivot}} = \sqrt {\frac{1}{{\frac{1}{{{W^2}}} + \frac{1}{{{{(2(AxisHeight - z)\tan (\alpha ))}^2}}}}}} \]
The mass flow rate in the pivoted area becomes:
\[{\dot m_{pivot}} = {C_d}\theta \int\limits_{z = h2}^{z = h4} {\rho \nu (z){W_{pivot}}dz} \]
It should be pointed out that the discharge coefficient is modulated based on opening factors, while opening width, opening height, and start height factor do not apply for a horizontally-pivoted window. The actual window width and height are used to calculate airflows for a horizontally-pivoted window.
The schematic drawing of air flow patterns through a simple vertical opening (AirflowNetwork:Multizone:Component:SimpleOpening) is shown in Figure 149. The equations used below are available from Walton (1989).
The air density for each node is assumed to be constant. The hydrostatic equation is used to relate pressures at various heights for each node:
\[Pn\left( y \right) = {P_{0n}} - {\rho _n}gy\]
\[{P_m}(y) = {P_{0m}} - {\rho _m}gy\]
where
P0n, P0m = pressure at nodes (zones) n and m at y = 0, the reference elevation of the opening [Pa]
ρn, ρm = air densities of zones n and m [kg/m3]
Pn, Pm = reference pressures of zones n and m [Pa]
It is assumed that the velocity of the airflow as a function of height is given by the orifice equation (Brown and Solvason 1962):
\[v(y) = {C_d}\sqrt {2\frac{{{P_n}(y) - {P_m}(y)}}{\rho }} \]
where
Cd= discharge coefficient [dimensionless]
ρ = density of the air going through the opening [kg/m3]
The neutral height, Y, where the velocity of the air is zero, may be calculated in the following equation:
\[Y = \frac{{{P_{on}} - {P_{0m}}}}{{g({\rho _n} - {\rho _m})}}\,\,\,\,\,\,\,\,or\,\,\,\,\,\,\,\frac{{{P_{om}} - {P_{0n}}}}{{g({\rho _m} - {\rho _n})}}\]
When the neutral plane is within the opening (first pattern in Figure 149), two-way (bi-directional) flows occur. The total flow through a large opening is the sum of both flows.
\[{\dot m_{0,Y}} = {C_d}\theta \int\limits_{y = 0}^{y = Y} {\rho \nu (y)Wdy} \]
\[{\dot m_{Y,H}} = {C_d}\theta \int\limits_{y = Y}^{y = H} {\rho \nu (y)Wdy} \]
When the neutral plane is below or above the large opening (second and third pattern in Figure 149), one-way flow occurs.
\[{\dot m_{}} = {C_d}\theta \int\limits_{y = 0}^{y = H} {\rho \nu (y)Wdy} \]
The opening width is modulated based on opening factors. A detailed description of opening factor calculations may be found in the Input Output Reference (Ref. AirflowNetwork:Multizone:Zone, AirflowNetwork:Multizone:Surface, and AirflowNetwork:Multizone:Component:DetailedOpening).
The above two objects are used to simulate airflows across large vertical openings. The simple opening component (AirflowNetwork:Multizone:Component:SimpleOpening) assumes the pressure difference across the opening is a function of height varied from opening bottom to top, so that two-way flow may be obtained if appropriate (Walton 1989). The Detailed Opening component (AirflowNetwork:Multizone:Component:DetailedOpening) assumes both the pressure difference across the opening and air density are a function of height, so that three-way flow may be obtained (COMIS 1990). If these opening models would be used for horizontal openings, the pressure difference across the opening and air density remain constant, so that only one-way flow is possible using the detailed and simple opening components which are meant for vertical or near-vertical openings. In reality, there are two-way flows (air recirculation) across a large horizontal opening caused by buoyancy due to temperature and pressure difference and forced flow driven by air pressure difference only. Therefore, a horizontal opening component (AirflowNetwork:Multizone:Component:HorizontalOpening) is available to simulate airflows across large horizontal openings with the possibility of two-way flow by combining forced and buoyancy airflows together.
The model for horizontal openings consists of forced airflow, buoyancy airflow, purge pressure and sloping plane. The model is mainly from a NIST report presented by Cooper (1989). The sloping plane (Bolmqvist and Sandberg 2004) portion of the model was added to allow for staircase simulations.
For simplicity, a two zone building (upper and lower zones) connected by a large horizontal opening is used to describe the model, as shown in Figure 150. Forced and buoyancy airflows are described separately below.
Forced airflows
The air mass flow rate is determined by the pressure difference across the opening. The relationship between pressure and airflow is the same as AIRNET for a component (see AirflowNetwork:Multizone:Component:SimpleOpening description above). Since the height of the opening is constant, the forced airflow is unidirectional. A positive value for pressure difference indicates flow direction is from the lower zone to the upper zone across the opening, while a negative value represents flow in the opposite direction. The following description addresses forced air mass flow rates and partial derivatives for three possible scenarios of pressure difference:
PL = PU
\({\mathop m\limits^ \bullet _{_U}}\)= \({\mathop m\limits^ \bullet _L}\)= 0
where:
PL= Air pressure in the lower zone [Pa]
PU= Air pressure in the upper zone [Pa]
\({\mathop m\limits^ \bullet _{_U}}\)= Air mass flow rate from the lower zone to the upper zone driven by forced airflow pressure difference [kg/s]
\({\mathop m\limits^ \bullet _L}\)= Air mass flow rate from the upper zone to the lower zone driven by forced airflow pressure difference [kg/s]
PL > PU
\({\mathop m\limits^ \bullet _{_U}}\)= 0
\[{\mathop m\limits^ \bullet _L} = {\rho _L}A{C_d}{\left( {\frac{{2\Delta P}}{{{\rho _{ave}}}}} \right)^{0.5}}\]
\[\frac{{{{\mathop {\partial m}\limits^ \bullet }_L}}}{{\partial {P_L}}} = {\rho _L}A{C_d}{\left( {\frac{1}{{2\Delta P{\rho _{ave}}}}} \right)^{0.5}}\]
where:
ρL= Air density in the lower zone [kg/m3]
A= Opening area [m2]
Cd= Discharge coefficient [Dimensionless]
ρave = Average air density between the lower and upper zones [Pa]
ΔP= Pressure difference PL -PU [Pa]
PL < PU
\({\mathop m\limits^ \bullet _{_L}}\)= 0
\[{\mathop m\limits^ \bullet _U} = - {\rho _U}A{C_d}{\left( {\frac{{2\left| {\Delta P} \right|}}{{{\rho _{ave}}}}} \right)^{0.5}}\]
\[\frac{{\partial \mathop {{m_U}}\limits^ \bullet }}{{\partial {P_L}}} = {\rho _U}A{C_d}{\left( {\frac{1}{{2\left| {\Delta P} \right|{\rho _{ave}}}}} \right)^{0.5}}\]
where:
ρU= Air density in the upper zone [kg/m3]
A= Opening area [m2]
Cd= Discharge coefficient [Dimensionless]
ρave = Average air density between the lower and upper zones [Pa]
ΔP= Pressure difference PL -PU [Pa]
Buoyancy airflows[LINK]
Buoyancy flow only occurs when the air density in the upper zone is greater than the air density in the lower zone. The flow is bi-directional and the amount of upper flow is equal to the lower flow across the opening. The following discussion assumes the air density in the upper zone is greater than the air density in the lower zone. Otherwise, the buoyancy flow rate is equal to zero. It is also assumed that the maximum buoyancy flow occurs when the pressure difference across the opening due to forced airflows is zero. The maximum buoyancy flow may be expressed as a part of Cooper’s model:
\[{\mathop m\limits^ \bullet _{buo,\max }} = {\rho _{ave}}*0.055{\left( {\frac{{g\left| {\Delta \rho } \right|D_{_H}^5}}{{{\rho _{ave}}}}} \right)^{0.5}}\]
where:
\({\mathop m\limits^ \bullet _{_{buo,\max }}}\) = Buoyancy mass flow rate at zero forced airflow pressure difference [kg/s]
g = Gravity acceleration [m/s2]
DH = Hydraulic diameter [m]
ρave = Average air density between the lower and upper zones [kg/m3]
Δρ = Density difference between the lower and upper zones [kg/m3]
Combined airflows[LINK]
When forced and buoyancy flows co-exist, it is possible to have either unidirectional or bi-directional flows. For example, when the upward force due to pressure difference is greater than the buoyancy force (downward), unidirectional flow occurs. Bi-directional flow only occurs when the upward imposed force is less than the buoyancy force. The critical pressure between unidirectional and bi-directional flows is called the purge pressure (Tan and Jaluria 1992). The purge pressure is a function of opening geometry and the buoyancy force (ΔP/(gΔρDH)) and may be expressed as (Cooper 1998):
\[\Delta {P_{Flood}} = \left| {C_{Shape}^2\frac{{g\Delta \rho D_H^5}}{{2{A^2}}}} \right|\]
where:
ΔPFlood = Purging pressure [Pa]
g= Gravity acceleration [m/s2]
DH= Hydraulic diameter [m]
A= Opening area [m2]
Δρ= Density difference between the lower and upper zones [kg/m3]
CShape = Shape factor [dimensionless]
\[C_{Shape}^{} = \left\{ \begin{array}{l}0.754\;\;\;{\rm{for}}\;{\rm{a}}\;{\rm{circle}}\;{\rm{opening}}\\0.942(w/D)\;{\rm{for}}\;{\rm{a}}\;{\rm{rectangle}}\;{\rm{opening}}\end{array} \right\}\]
where:
w= Opening width [m]
D= Opening depth [m]
As mentioned above, when the air pressure difference between two zones is zero there is the maximum bi-directional flow due to the buoyancy force. When the pressure difference increases from 0 and is less than |ΔPFlood|, there is some bi-directional flow across the opening, but less than the maximum flow. If the pressure difference keeps increasing and exceeds |ΔPFlood|, there is no bi-directional flow. Cooper’s model assumes the buoyancy flow varies linearly with pressure difference.
\[{\mathop m\limits^ \bullet _{buo}} = \left\{ \begin{array}{l}{\mathop m\limits^ \bullet _{buo,\max }} * \left( {1 - \frac{{\left| {\Delta P} \right|}}{{\left| {\Delta {P_{Flood}}} \right|}}} \right)\;\,\,If\;\Delta \rho > 0\;and\;\frac{{\left| {\Delta P} \right|}}{{\left| {\Delta {P_{Flood}}} \right|}} < 1\\0\;\,\,\,Otherwise\end{array} \right\}\]
The total air flow across the opening is based on superposition of the forced and buoyancy flows, and may be expressed for three different pressure difference scenarios as follows:
PL = PU
\[{\mathop m\limits^ \bullet _{_U}} = {\mathop m\limits^ \bullet _{buo}}\]
\[{\mathop m\limits^ \bullet _{_L}} = {\mathop m\limits^ \bullet _{buo}}\]
\[\frac{{{{\mathop {\partial m}\limits^ \bullet }_L}}}{{\partial {P_L}}} = 0\]
PL > PU
\[{\mathop m\limits^ \bullet _{_U}} = {\mathop m\limits^ \bullet _{buo}}\]
\[{\mathop m\limits^ \bullet _L} = {\rho _L}A{C_d}{\left( {\frac{{2\Delta P}}{{{\rho _{ave}}}}} \right)^{0.5}} + {\mathop m\limits^ \bullet _{buo}}\]
\[\frac{{{{\mathop {\partial m}\limits^ \bullet }_L}}}{{\partial {P_L}}} = {\rho _L}A{C_d}{\left( {\frac{1}{{2\Delta P{\rho _{ave}}}}} \right)^{0.5}} - \frac{{{{\mathop m\limits^ \bullet }_{buo,\max }}}}{{\left| {\Delta {P_{Flood}}} \right|}}\]
where:
ρL= Air density in the lower zone [kg/m3]
A= Opening area [m2]
Cd= Discharge coefficient [dimensionless]
Ρave = Average air density between the lower and upper zones [Pa]
ΔP= Pressure difference PL -
PU [Pa]PL < PU
\[{\mathop m\limits^ \bullet _{_L}} = {\mathop m\limits^ \bullet _{buo}}\]
\[{\mathop m\limits^ \bullet _U} = - {\rho _U}A{C_d}{\left( {\frac{{2\left| {\Delta P} \right|}}{{{\rho _{ave}}}}} \right)^{0.5}} + {\mathop m\limits^ \bullet _{buo}}\]
\[\frac{{\partial \mathop {{m_U}}\limits^ \bullet }}{{\partial {P_L}}} = {\rho _U}A{C_d}{\left( {\frac{1}{{2\left| {\Delta P} \right|{\rho _{ave}}}}} \right)^{0.5}} + \frac{{{{\mathop m\limits^ \bullet }_{buo,\max }}}}{{\left| {\Delta {P_{Flood}}} \right|}}\]
Sloping plane
When a staircase is introduced as shown in Figure 151, the effective opening area will be used to replace A (opening area) in the above equations. The effective area may be estimated as (Bolmqvist and Sandberg 2004):
\[{A_{eff}} = A\sin \alpha (1 + \cos \alpha )\]
where:
Aeff= Effective area of horizontal opening [m2]
A= Area of horizontal opening [m2]
α= Angle between the stair plane and horizontal opening plane [degrees]
Note: the hydraulic diameter calculation is based on the effective opening area, while the opening depth remains the same.
Figure 152 demonstrates possible forced and buoyancy flow rates at different ratios of pressure difference to purging pressure across a horizontal opening when the upper zone air density is greater than the lower zone air density. The pressure difference is the lower zone pressure minus the upper zone pressure. Otherwise, the buoyancy flow is zero. In addition, when the absolute ratio is above 1, the buoyancy flow is also zero. The following table provides a brief description for the legend listed in Figure 152.
Wind pressure calculations[LINK]
The wind pressure is determined by Bernoulli’s equation, assuming no height change or pressure losses:
\[{p_w} = {C_p}\rho \frac{{{V_{ref}}^2}}{2}\]
where:
pw= Wind surface pressure relative to static pressure in undisturbed flow [Pa]
ρ= Air density [kg/m3]
Vref= Reference wind speed at local height [m/s]
Cp= Wind surface pressure coefficient [dimensionless]
Vref may be expressed as (Ref, Local Wind Speed Calculations):
\[{V_{ref}} = {V_{met}}{\left( {\frac{{{\delta _{met}}}}{{{z_{met}}}}} \right)^{{\alpha _{met}}}}{\left( {\frac{z}{\delta }} \right)^\alpha }\]
Cp is a function of location on the building envelope and wind direction. When Wind Pressure Coefficient Type = “INPUT”, the Cp values are explicitly defined in the input for AirflowNetwork:Multizone:Wind Pressure Coefficient Values. When Wind Pressure Coefficient Type = “AVERAGE-SURFACE CALCULATION” and the building shape is rectangular, the program uses the following equations to calculate wind pressure coefficient (Cp) values for different wind directions. For a low rise building, the normalized surface pressure coefficient may be written as (Swami and Chandra 1988):
\[{C_{p,n}} = 0.6*\ln \left[ \begin{array}{l}1.248 - 0.703\sin (\alpha /2) - 1.175{\sin ^2}(\alpha ) + 0.131{\sin ^3}(2\alpha G)\\ + \,0.769\cos (\alpha /2) + 0.07{G^2}{\sin ^2}(\alpha /2) + 0.717{\cos ^2}(\alpha /2)\end{array} \right]\]
where
Cp,n = Cp value at a given angle between wind direction and the outward normal of the surface under consideration [dimensionless]
α = Angle between wind direction and outward normal of wall under consideration [deg]
G = Natural log of the ratio of the width of the wall under consideration to the width of the adjacent wall [dimensionless]
n = Index of incident angle at 30-degree increments
For walls of a high rise building, a two-dimensional array of surface-averaged wind pressure coefficients is generated based on wind incident angle and side ratio. The wind pressure coefficients are provided in 2001 ASHRAE Fundamentals Handbook, p. 16.5, Fig. 7, “Surface Averaged Wall Pressure Coefficients for Tall Buildings”. The original work was performed by Atkins et al. (1979). The incident angle has an increment of 30 degrees. The side ratio values are 0.25, 1.0, and 4.0. For a given incident angle and building aspect ratio, the program uses linear interpolation to calculate the corresponding wind pressure coefficient Cp,n.
For the roof of a high rise building, a two-dimensional array of surface-averaged wind pressure coefficients is also generated based on wind incident angle and side ratio. The wind pressure coefficients are provided in 2001 ASHRAE Fundamentals Handbook, p. 16.6, Fig. 9, “Surface Averaged Roof Pressure Coefficients for Tall Buildings”. The original work was performed by Holmes (1986). The incident angle has an increment of 30 degrees. The side ratio values are 0.25, 1.0, and 4.0. At a given wind incident angle and building aspect ratio, the program uses linear interpolation to calculate the corresponding wind pressure coefficient Cp,n.
The wind surface pressure at the given incident angle can be calculated by combining the above two equations:
\[{p_{w,n}} = {C_{p,n}}\rho \frac{{{V_{ref}}^2}}{2}\]
Solution method[LINK]
Based on the relationship between airflow rate and pressure drop for each component, a system of equations for all components can be assembled together in an n x n square matrix, where n is the number of nodes. Newton’s method is used to iteratively solve for the air pressure at each node. A new estimated vector for all node pressures, {P}*, is computed from the current estimated vector of node pressures, {P}, by:
\[{\left\{ P \right\}^*} = \left\{ P \right\} - \left\{ C \right\}\]
where the correction vector, {C}, is computed by the matrix relationship:
\[\left[ J \right]\left\{ C \right\} = \left\{ B \right\}\]
{B} is a column vector with each component given by:
\[{B_n} = \sum\limits_i {\mathop {{m_i}}\limits^ \bullet }\]
where n is the node number and i indicates all flow paths connecting node n to other nodes, and [J] is the square Jacobian matrix whose elements are given by:
\[{J_{n,m}} = \sum\limits_i {\frac{{\partial \mathop m\limits^ \bullet }}{{\partial {P_m}}}}\]
Convergence acceleration[LINK]
The convergence tolerance is used to check the sum of the mass flow rates by applying mass conservation. The convergence acceleration equation shown below is used to correct the node pressures to more rapidly obtain a solution. By assuming a constant ratio of correction values from one iteration to the next, the following method is applied:
\[{P_n}^* = {P_n} - {C_n}/(1 - r)\]
where
r= the ratio of Cn for the current iteration to its value for the previous iteration [dimensionless]
Cn= Correction value at the nth node [Pa]
Pn= Estimated pressure at the nth node [Pa]
Pn*= Corrected pressure at the nth node used in the next iteration [Pa]
This method is similar to a Steffensen iteration (Conte and de Boor 1972) which is used as a fixed-point iteration method for individual nonlinear equations.
The iteration correction method presented in the above equation gives a variable factor. When the solution is close to convergence, the solution method converges quadratically. By limiting cases where the value of r is less than some value, such as -0.5, the solution will not interfere with the rapid convergence. It has not been proven that the convergence acceleration equation will always lead to convergence, but it can be shown that it will not prevent convergence. Newton’s method converges when the estimated solution values are within some distance, called the radius of convergence, or the correct solution. Applying the convergence acceleration equation when -1 < r < 0, will cause a smaller correction than Newton’s method, which therefore, can not force the iterations outside the radius of convergence. When r < -1, the solution diverges in an oscillatory fashion. When r > 1, the solution also diverges, but in a nonoscillatory manner. For 0 < r < 1, the solution is approached from one direction. In all three cases, the convergence acceleration equation applies as long as r is truly constant over several iterations. However, for the last case, this involves a true extrapolation of correction factor which is very sensitive to the accuracy of r. This is most extreme for the case of r=1, which would cause an infinite correction.
Node Temperature Calculations[LINK]
A brief description of the air node temperature calculation is given below. A detailed description can be found in the work of Swami et al. (1992). The following equation is used to calculate temperature distribution across a duct element at the given airflow rate and inlet air temperature:
\[\mathop m\limits^ \bullet {C_p}\frac{{dT}}{{dx}} = UP(T{}_\infty - T)\]
where
Cp = Specific heat of airflow [J/kg•K]
\(\mathop m\limits^ \bullet\) = Airflow rate [kg/s]
P = Perimeter of a duct element [m]
T = Temperature as a field variable [°C]
\(T{}_\infty\) = Temperature of air surrounding the duct element [°C]
U = Overall heat transfer coefficient [W/m2•K]
\[U = \frac{1}{{\frac{1}{{{h_i}}} + \frac{1}{{{h_o}}} + \sum {\frac{{{t_j}}}{{{k_j}}}} }}\]
hi = Inside heat transfer coefficient [W/m2•K]
ho = Outside heat transfer coefficient [W/m2•K]
tj = Thickness at j-th layer [m]
kj = Thermal conductivity at j-th layer [W/m•K]
The outlet air temperature at the end of the duct (x=L) is:
\[{T_o} = {T_\infty } + ({T_i} - {T_\infty })*\exp \left( { - \frac{{UA}}{{\mathop {m{C_p}}\limits^ \bullet }}} \right)\]
where
Ti = Inlet air temperature [°C]
To = Outlet air temperature [°C]
T∞ = Temperature of air surrounding the duct element [°C]
A = Surface area (Perimeter * Length) [m2]
The heat transfer by convection to ambient, Q, is:
\[Q = \mathop m\limits^ \bullet {C_p}({T_\infty } - {T_i})\left[ {1 - \exp \left( { - \frac{{UA}}{{\mathop {m{C_p}}\limits^ \bullet }}} \right)} \right]\]
The outlet air temperature can be calculated using the above equation at the given inlet air temperature. Since the inlet temperature at one linkage is the outlet temperature for the connected linkage, the outlet air temperatures at all nodes are solved simultaneously. A square linear system assembled by the AirflowNetwork model is expressed below:
\[\{ M\} [T] = [B]\]
where
{M}= Airflow matrix
[T]= Temperature vector
[B]= Given boundary conditions
The zone air temperatures and primary air loop component (fan and coils) outlet conditions are used as prescribed conditions in the AirflowNetwork model. In addition, the temperature difference across zone loop components (terminal units) is held constant during the calculations. For example, thermal zone temperatures calculated during the previous system time step are used as prescribed temperatures when calculating all other node temperatures. The zone air temperature is assumed constant (prescribed) throughout the AirflowNetwork iterative solution. The fan and coil outlet air temperatures, and terminal unit temperature differences are assumed constant within an AirflowNetwork iteration. The sensible heat gains calculated during the AirflowNetwork solution are then used to predict a new zone air temperature.
Node Humidity Ratio Calculations[LINK]
A brief description of the air node humidity ratio calculation is given below. A detailed description can found in the work of Swami et al. (1992). The following equation is used to calculate humidity ratio distribution across a duct element at the given airflow rate and inlet air humidity ratio:
\[\mathop m\limits^ \bullet \frac{{dW}}{{dx}} = {U_m}P(W{}_\infty - W)\]
where
\(\mathop m\limits^ \bullet\) = Airflow rate [kg/s]
P= Perimeter of a duct element [m]
W = Humidity ratio [kg/kg]
\(W{}_\infty\) = Humidity ratio of air surrounding the duct element [kg/kg]
Um = Overall moisture transfer coefficient [kg/m2•s]
\[{U_m} = \frac{1}{{\frac{1}{{{h_{m,i}}}} + \frac{1}{{{h_{m,o}}}} + \sum {\frac{{{t_j}}}{{{D_j}}}} }}\]
hm,i = Inside moisture transfer coefficient [kg/m2•s]
hm,o = Outside moisture transfer coefficient [kg/m2•s]
tj = Thickness at j-th layer [m]
Dj = Moisture diffusivity at j-th layer [kg/m•s]
The outlet air humidity ratio at the end of the duct (x=L) is:
\[{W_o} = {W_\infty } + ({W_i} - {W_\infty })*\exp \left( { - \frac{{{U_m}A}}{{\mathop m\limits^ \bullet }}} \right)\]
where
Wi= Inlet air humidity ratio [kg/kg]
Wo= Outlet air humidity ratio [kg/kg]
A= Surface area (Perimeter * Length) [m2]
The moisture transfer by convection to ambient, Qm, is
\[{Q_m} = \mathop m\limits^ \bullet ({W_\infty } - {W_i})\left[ {1 - \exp \left( { - \frac{{{U_m}A}}{{\mathop m\limits^ \bullet }}} \right)} \right]\]
The outlet air humidity ratio can be calculated using the above equation at the given inlet air humidity ratio. Since the inlet humidity ratio at one linkage is the outlet humidity ratio for the connected linkage, the outlet air humidity ratio at all nodes are solved simultaneously. A square linear system assembled by the AirflowNetwork model is expressed below:
\[\{ {M_m}\} [W] = [{B_m}]\]
where
{Mm}= Airflow matrix
[W]= Humidity ratio vector
[Bm]= Given boundary conditions
The zone air humidity ratios and primary air loop component (fan and coils) outlet conditions are used as prescribed conditions in the AirflowNetwork model. For example, thermal zone humidity ratios calculated during the previous system time step are used as prescribed humidity ratios when calculating all other node humidity ratios. The zone air humidity ratio is assumed constant (prescribed) throughout the AirflowNetwork iterative solution. The coil outlet air humidity ratio is assumed constant within an AirflowNetwork iteration. The latent heat gains calculated during the AirflowNetwork solution are then used to predict a new zone air humidity ratio.
Sensible and Latent Load Calculations[LINK]
The zone sensible and latent loads calculated in the AirflowNetwork model consist of multizone, duct conduction and leakage. In addition, the impact of infiltration and mixing is accounted for in this calculation. The multizone load only includes incoming airflows from outside (infiltration) and other adjacent zones (mixing) with and without forced-air fan operation. It is divided into two terms: variable and constant. The constant term is the sum of the mass flow rate multiplied by the specific heat for both infiltration and mixing. The variable term includes the impact of zone and outdoor air temperature. Each of these terms is used in the zone energy balance equation. The sensible load items from the multizone load calculations may be written as follows:
where
MCPairflow = Sum of air mass flow rate multiplied by specific heat for infiltration and mixing [W/K]
MCPTairflow = Sum of air mass flow rate multiplied by specific heat and temperature for infiltration and mixing [W]
= Incoming air mass flow rate from outdoors [kg/s]
= Incoming air mass flow rate from adjacent zones [kg/s]
Tamb= Outdoor air dry-bulb temperature [°C]
Tzone= Adjacent zone air temperature [°C]
The latent load items from multizone load calculations may be written as follows:
where
Mairflow = Sum of air mass flow rates for infiltration and mixing [kg/s]
MWairflow = Sum of air mass flow rate multiplied by humidity ratio for infiltration and mixing [kg/s]
= Incoming air mass flow rate from outdoors [kg/s]
= Incoming air mass flow rate from adjacent zones [kg/s]
Wamb= Outdoor air humidity ratio [kg/kg]
Wzone= Adjacent zone air humidity ratio [kg/kg]
The air distribution system (ADS) loads due to duct conduction and leakage depend on the air distribution system component (e.g., duct) location. The air distribution system sensible and latent loads are calculated for each zone as follows:
\[{Q_{ADS,i}} = \sum\limits_j {{Q_{cond(i,j)}} + } \sum\limits_j {{Q_{leak(i,j)}}} \]
\[{Q_{ADS,m,i}} = \sum\limits_j {{Q_{cond,m(i,j)}} + } \sum\limits_j {{Q_{leak,m(i,j)}}} \]
where
QADS,i= Total sensible load in the i-th zone due to ADS losses [W]
Qcond(ij)= Duct wall conduction loss at the j-th duct located in the i-th zone [W]
Qleak(ij)= Sensible supply leak loss at the j-th linkage located in the i-th zone [W]
QADS,m,i= Total latent load in the i-th zone due to ADS losses [kg/s]
Qcond,m(ij)= Duct wall vapor diffusion loss at the j-th duct located in the i-th zone [kg/s]
Qleak,m(ij) = Latent supply leak loss at the j-th linkage located in the i-th zone [kg/s]
Impacts of Supply Air Constant Volume Fan Control on Load: Cycling vs. Continuous[LINK]
The AirflowNetwork model currently allows two types of constant volume fans: Fan:ConstantVolume and Fan:OnOff. The Fan:ConstantVolume object has only one type of supply air fan operation mode: continuous fan, cycling compressor (ContinuousFanWithCyclingCompressor). However, the Fan:OnOff has two types of supply air fan operation modes: cycling fan, cycling compressor (CyclingFanAndCompressor) or continuous fan, cycling compressor (ContinuousFanWithCyclingCompressor). The CyclingFanAndCompressor operation mode is frequently referred to as “AUTO fan”, where the compressor(s) and supply air fan operate in unison to meet the zone heating or cooling load, and cycle off together when the heating or cooling load has been met. The ContinuousFanWithCyclingCompressor operation mode is often referred to as “fan ON”, where the compressor(s) cycle on and off to meet the zone heating or cooling load but the supply air fan operates continuously regardless of compressor operation. The supply air fan operation mode is specified in an HVAC system object based on a given fan operation mode schedule (e.g., AirLoopHVAC:UnitaryHeatCool object).
The determination of the zone sensible and latent loads caused by multizone airflows and forced air distribution system operation is dependent on the supply air fan operation mode (see Sensible and Latent Load Calculations section above). The zone loads calculated by the AirflowNetwork model are added to zone sensible and latent balances in the ZonePredictorCorrector module to calculate zone air temperatures and humidity ratios (see Integration of the AirflowNetwork Model section below). For the case of continuous fan/cycling compressor, the zone loads during forced air distribution system operation are calculated with the system design air mass flow rate without modification, since the system air node conditions (temperature and humidity) reflect the average values for the time step considering the coil/fan on and off periods during the time step.
For the case of cycling fan/cycling compressor, where the forced air distribution system can operate for a portion of the simulation time step, the airflows are determined based on the air distribution system part-load ratio (ratio of the average air mass flow rate for the time step divided by the design air mass flow rate). The airflows for the AirflowNetwork:Distribution:Linkage objects are reported during the air distribution system on cycle, since no airflow is assumed during the system off cycle. The airflows for the AirflowNetwork:Multizone:Surface objects are weighted by the air distribution system part-load ratio. The zone loads are the sum of energy losses during both the air distribution system on and off cycle at each system time step. The energy losses when the air distribution system is on are calculated using the system “on” air flow rate multiplied by the air distribution system run time fraction. The energy losses when the air distribution system is off are obtained from the multizone airflow calculations (without forced air through the air distribution system) and these losses are multiplied by (1.0 - system run time fraction), assuming no airflows through the air distribution system when the fan is off. The formulas used to calculate airflows and zone loads are given below:
Airflow
Airflow = Airflow during ADS on cycle * ADS Part-load ratio + Airflow during ADS off cycle * (1.0 - ADS Part-load ratio)
where ADS = Air Distribution System
Zone load
System run time fraction = Maximum run time fraction of coils and fans in the air distribution system
Zone energy losses = Zone energy loss during ADS on cycle * System run time fraction + Zone energy loss during ADS off cycle * (1.0 - System run time fraction)
The calculation of loads due to multizone airflow, without forced air distribution system operation, is performed when the HVAC system is first simulated during a simulation time step (FirstHVACIteration = True). The calculation of loads due to air distribution system operation is performed on subsequent iterations within the same simulation time step when the mass flow rate at the supply air fan inlet is greater than 0.0 (and FirstHVACIteration = False).
Airflow Calculation Procedure using A Supply Variable Air Volume Fan[LINK]
The AirflowNetwork model currently also allows a variable air volume fan type as Fan:VariableVolume. The allowed terminal type is AirTerminal:SingleDuct:VAV:Reheat only. Other types of terminals will be added later.
In general, the supply fan air flow rate in a VAV central system is determined by a sum of terminal flow rates when the AirflowNetwork model is not applied. When the AirflowNetwork model is applied and the supply air fan flow rate is given, each terminal flow is determined by pressure resistance of each supply air pathway. It is possible that the delivered air flow rate from the pressure resistance at each terminal may be totally different from the desired flow rate determined by terminal units. Therefore, it is not easy to meet both requirements. The following compromised approach, including possible supply and return leaks in an air distribution system, is implemented.
Set up terminal airflows in the AirflowNetwork module based on the SimVAV subroutine in the HVACSingleDuctSystem module.
Require AirflowNetwork:Distribution:Component:LeakageRatio objects to define supply leaks, so that the values of the Effective Leakage Ratio field are used to decide the supply fan flow rates. The base of the ratio will be actual supply fan flow rate, instead of the maximum fan flow rate used in the constant volume fan systems.
Assign the supply fan airflow rate based on the sum of all terminal flow rates and all supply leak ratios until it reaches the maximum fan flow rate
\[{\dot m_{fan}} = \frac{{\mathop \sum \nolimits_i^n {{\dot m}_{i,terminal}}}}{{1 - \mathop \sum \nolimits_j^k {F_j}}}\]
where
\({\dot m_{fan}}\)= The supply fan flow rate
\({\dot m_{i,terminal}}\)= The flow rate at the ith terminal, which is determined in the subroutine SimVAV in the HVACSingleDuctSystem module
n= Number of terminals
\({F_j}\)= The fraction of the supply fan flow rate at the jth supply leak, given in the AirflowNetwork:Distribution:Component:LeakageRatio objects.
k= Number of supply leaks
If the calculated supply fan flow rate is above the maximum limit of the supply fan flow rate, and the supply leak ratios remain the same, the supply fan flow rate is set to the maximum limit, and the terminal flow rates are reduced proportionally weighted by a ratio of the maximum supply fan flow rate by input to the calculated supply fan flow rate. Therefore, a sum of all terminal rates and all supply leak rates is equal to the maximum supply fan rate.
\[R = \frac{{{{\dot m}_{fan,max}}}}{{{{\dot m}_{fan,cal}}}}\]
\[{\dot m_{i,terminal,final}} = {\dot m_{i,terminal}}*R\]
where
R = The ratio of the maximum fan flow rate given in the inputs to the requested fan flow rate based on the above equation
\({\dot m_{fan,max}}\) = The maximum supply fan flow rate by input
\({\dot m_{fan,cal}}\) = The calculated supply fan flow rate
\({\dot m_{i,terminal,final}}\) = The final flow rate at each terminal adjusted by the ratio
Integration of the AirflowNetwork Model[LINK]
The loads calculated by the AirflowNetwork model are integrated into the EnergyPlus heat balance equation in a similar manner as described elsewhere in this document in the section “Basis for the Zone and System Integration”. The mass flow rate summations and sensible and latent loads described in the previous section are included in the calculation of zone temperature and humidity ratio.
The revised zone temperature update equation becomes:
\[T_z^t = \frac{{\sum\limits_{i = 1}^{{N_{sl}}} {\mathop {{Q_i}}\limits^ \bullet } + \sum\limits_{i = 1}^{{N_{surfaces}}} {{h_i}} {A_i}{T_{si}} + MCPTairflow + {{\mathop m\limits^ \bullet }_{sys}}{C_p}{T_{{\rm{supply}}}} + QADS,z - \left( {\frac{{{C_z}}}{{\delta t}}} \right)\left( { - 3T_z^{t - \delta t} + \frac{3}{2}T_z^{t - 2\delta t} - \frac{1}{3}T_z^{t - 3\delta t}} \right)}}{{\left( {\frac{{11}}{6}} \right)\frac{{{C_z}}}{{\delta t}} + \sum\limits_{i = 1}^{{N_{surfaces}}} {{h_i}} A + MCPairflow + {{\mathop m\limits^ \bullet }_{sys}}C}}\]
Where MCPTairflow is the sum of mass flow rate multiplied by specific heat and temperature for infiltration and mixing, QADS,z is the added total sensible load in the zone due to Air Distribution System losses, and MCPairflow is the sum of mass flow rate multiplied by specific heat for infiltration and mixing as calculated in the AirflowNetwork model described above.
The revised coefficient (B) used in the zone humidity ratio calculation is shown below:
\[B = \sum {k{g_{mas{s_{{\rm{Sched Loads}}}}}}} + M{W_{airflow}} + {\dot m_{sy{s_{in}}}}{W_{sys}} + \sum\limits_{i = 1}^{surfs} {{A_i}{h_{mi}}{\rho _{ai{r_z}}}{W_{surf{s_i}}}} + QADS,m,z\]
Where MWairflow is the sum of mass flow rate multiplied by humidity ratio for infiltration and mixing and QADS,m,z
is theadded total latent (moisture) load in the zone due to Air Distribution System losses from the AirflowNetwork model described above. This coefficient is used in the prediction of moisture as described in the section “Moisture Predictor-Corrector” found elsewhere in this document.Model Output[LINK]
The available outputs from the AirflowNetwork model are described in the EnergyPlus Input Output Reference manual.
Occupant Ventilation Control[LINK]
The AirflowNetwork:OccupantVentilationControl object enhances the AirflowNetwork model and provides more practical and advanced controls for window opening and closing operations, based on Marais & Teichmann (2014). This control includes minimum open and closed time control, indoor thermal comfort control, and opening and closing probability controls.
Procedures of occupant ventilation control[LINK]
It should be pointed that the open elapsed time and closed elapsed time are not independent. In other words, when one of the elapsed times value is greater than 0, then other elapsed time value must be equal to 0. The model is either tracking a window as open, in which case the open elapsed time grows, or as closed, in which case the closed elapsed time grows.
The calculation procedures are presented in the following steps:
Step 1: Open elapsed time check
The model checks the open elapsed time first. When the open elapsed time is greater than 0 and less than the minimum opening time, a window will remain open. Otherwise, the model goes to Step 2.
Step 2: Closed elapsed time check
This step checks the closed elapsed time. When the closed elapsed time is greater than 0 and less than the minimum closed time, a window will remain closed. Otherwise, the model goes to Step 3.
Step 3: Elapsed time only?
When either the open elapsed time or the closed elapsed time is long enough (greater than the minimum time) and no other checks are needed, the model returns to the other ventilation control defined in the Ventilation Control Mode field in AirflowNetwork:MultiZone:Zone or AirflowNetwork:MultiZone:Surface. If other checks are needed, the following steps will be performed.
Step 4: Thermal comfort temperature calculation
The thermal comfort check requires the thermal comfort temperature and the comfort band. The comfort temperature of Tcomf is calculated as a function of the outdoor dry-bulb temperature Tout. The comfort temperature calculation may be based on two curves and a boundary temperature point.
\[T_{comf} = \left\{ \begin{array}{ll} \rm{Low Temp Curve,} & \min T_{out} < x \leq \rm{Boundary Point} \\ \rm{High Temp Curve,} & \rm{Boundary Point} \leq x < \max T_{out} \end{array}\right.\]
Step 5: Thermal band calculation
The comfort band is calculated as a function of the design personal dissatisfaction PPD.
\[\theta = -0.0028 \left(100-\rm{PPD}\right)^2 + 0.3419 \left(100-\rm{PPD}\right) – 6.6275\]
This equation is valid for PPD 0% to 35%
Step 6: Upper boundary check of thermal comfort
After calculation of the comfort temperature and band, the thermal comfort check will be performed using the zone air operative temperature Tg. The check consists of upper and lower boundary checks.
The upper boundary check checks window opening status and uses the following logic:
\[T_g > \left(T_{comf} + \theta\right)\]
If the above logic check is true, the opening probability check will be performed. The detailed description is given in the opening probability section.
If the above logic check is false, no action is needed.
If the opening probability check is true and upper boundary check is satisfied, a window will be opened, regardless of open/closed status at the previous time step. If the opening probability check is false and upper boundary check is satisfied, a window will remain at the status from the previous time step.
\[T_g > \left(T_{comf} + \theta\right) \&\& \rm{OpeningProbability}\]
Step 7: Lower boundary check of thermal comfort
The lower boundary check will follow the upper boundary check to check the window closing status, using the following logic.
\[T_g < \left(T_{comf} - \theta\right)\]
If the above logic check is true, the closing probability check will be performed. The detailed description is given in the closing probability section.
If the above logic check is false, no action is needed.
If the closing probability check is true and lower boundary check is satisfied, a window will be closed, regardless of open/closed status at the previous time step. If the closing probability check is false and lower boundary check is satisfied, a window will remain at the status from the previous time step.
\[T_g < \left(T_{comf} - \theta\right) \&\& \rm{ClosingProbability}\]
The output variables from the model are open status, opening probability status, and closing probability status. The detailed description of opening status is given in the Airflow Network Outputs section in the Input Output Reference.
Note: The upper and lower boundary checks are not independent. In other words, when one of boundary check is satisfied, the other check will be dissatisfied.
Procedures of opening probability control[LINK]
Opening probability control provides an optional random number check. The control logic of opening probability is described in the following steps.
Step 1: Closed elapsed time check
This check requires that closed elapsed time is longer than the minimum closing time.
Closed time > minimum closed time
If the time duration is not long enough, the output is false so that a window remains closed.
If the time duration is long enough, an occupancy check is performed.
Step 2: Occupancy check
If a zone is not occupied and the occupancy check is requested, the output is false. If the zone is occupied, the next step is to check the setpoints using zone air temperature at the previous time step as a reference.
Step 3: Setpoint check
There are 5 temperature control types. The following types are available:
No control: Bypass
Single heating setpoint: If Tzon > setpoint, go to next step. Otherwise, return false.
Single cooling setpoint: If Tzon < setpoint, go to next step. Otherwise, return false.
Single heating and cooling setpoint: no action by returning false
Dual heating and setpoints: If heating setpoint < Tzon < Cooling setpoint, go to next step. Otherwise, return false.
Step 4: Select bypass or opening probability check
A choice is provided at this stage so that the opening probability check may be performed or bypassed. If bypassed, the output will be true to open a window. If performed, probability will be determined from a schedule.
Step 5: Perform opening probability check
The opening probability (OP) value is determined from a schedule.
OP = Schedule value or specific function
If the probability value is greater than a random number, the output is true.
OP > random number (random number is between 0 and 1)
Otherwise, the result will be false.
Procedures of closing probability control[LINK]
The control logic of closing probability is described as follows.
Step 1: Open elapsed time check
This check requires that open elapsed time is longer than the minimum opening time.
Open time > minimum open time
If the time duration is not long enough, the output is false and the window remains open. Otherwise, Step 2 will be performed.
Step 2: Select by-pass or closing probability check
A choice is provided at this step so that the closing probability check may be performed or bypassed. If bypassed, the output will be true to close a window. If performed, the closing probability will be calculated. The closing probability (CP) value is given from a schedule.
CP = Schedule value or specific function
If the closing probability check is performed and the closing probability is greater than a random number, the output is true.
CP > random number (random number is between 0 and 1)
Otherwise, the result will be false.
References[LINK]
Atkins, R. E., J. A. Peterka, and J. E. Cermak. 1979. “Averaged pressure coefficients for rectangular buildings,” Wind Engineering, Proceedings of the Fifth International Conference 7:369-80, Fort Collins, CO. Pergamon Press, NY.
Bolmqvist, C. and M. Sandberg, 2004, “Air Movements through Horizontal Openings in Buildings - A Model Study,” International Journal of Ventilation, Vol. 3, No. 1, pp. 1-9
COMIS Fundamentals. 1990. Edited by Helmut E. Feustel and Alison Rayner-Hooson, LBL-28560, Lawrence Berkeley Laboratory, Berkeley, CA
Conte, S. D. and C de Boor. 1972. Elementary Numerical Analysis: an Algorithmic Approach, McGraw-Hill.
Cooper, L., 1989, “Calculation of the Flow Through a Horizontal Ceiling/Floor Vent,” NISTIR 89-4052, National Institute of Standards and Technology, Gaithersburg, MD
Dols, W. S. & G. N. Walton. 2002. “CONTAMW 2.0 User Manual,” NISTIR 6921, National Institute of Standards and Technology, Gaithersburg, Maryland
Holmes, J. D. 1986. Wind Loads on low-rise buildings: The structural and environmental effects of wind on buildings and structures, Chapter 12, Faculty of Engineering, Monash University, Melbourne, Australia
Marais, J. M. & C. Teichmann, “Window Simulation Methods Required for Manual Window Ventilated Buildings,” Fifth German-Austrian IBPSA Conference, September 22-24, 2014, RWTH Aachen University, Germany
Swami, M. V. and S. Chandra. 1988. Correlations for pressure distribution on buildings and calculation of natural-ventilation airflow, ASHRAE Transactions 94(1988) (Pt 1), pp. 243-266.
Swami, M. V., L. Gu, & V. Vasanth. 1992. “Integration of Radon and Energy Models for Building,” FSEC-CR-553-92, Florida Solar Energy Center, Cocoa, Florida
Tan, Q. and Y. Jaluria, 1992, “Flow through Horizontal Vents as Related to Compartment Fire Environments,” NIST-GCR-92-607, National Institute of Standards and Technology, Gaithersburg, Maryland
Walton, G. N. 1989. “AIRNET - A Computer Program for Building Airflow Network Modeling,” NISTIR 89-4072, National Institute of Standards and Technology, Gaithersburg, Maryland
Walton, G. N. & W. S. Dols. 2003. “CONTAM 2.1 Supplemental User Guide and Program Documentation,” NISTIR 7049, National Institute of Standards and Technology, Gaithersburg, Maryland
Documentation content copyright © 1996-2026 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. EnergyPlus is a trademark of the US Department of Energy.
This documentation is made available under the EnergyPlus Open Source License v1.0.