AirflowNetwork Model[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.
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:
˙mi=CiρΔPiμ
where:
˙mi

= Air mass flow rate at ith linkage [kg/s]

Ci

= Air mass flow coefficient [m3]

ΔPi

= Pressure difference across the ith linkage [Pa]

μ

= Air viscosity [Pas]

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=∑˙mi∑˙mi
Absolute airflow tolerance=∣∣∑˙mi∣∣
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:
ΔP=(Pn+ρV2n2)−(Pm+ρV2m2)+ρg(zn−zm)
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/m3]

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:
ΔP=Pn−Pm+PS+PW
where:
Pn,Pm

= 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 1. The equations used below are extracted from the COMIS Fundamentals manual (1990).
The air density is assumed to be a linear function of height:
ρi(z)=ρ0i+biz
The pressure difference is assumed to be linear and simulate the effect of turbulence:
ΔPt=Pt0+btz
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:
P1(z)−P2(z)=(P01−P02)−g[(ρ01z+b1z2/2)−(ρ02z+b2z2/2)]+(Pt0+btz)
The velocity at any level z is given by
v(z)=√2P1(z)−P2(z)ρ
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(b1−b2)z2/2+[g(ρ01−ρ02)−bt]z+(−P01+P02−Pt0)=0
This equation above can have two, one, or zero real solutions. The zero solution represents a oneway flow through the opening and may be expressed in the following equation:
˙m=Cdθz=H∫z=0ρν(z)Wdz
The one real solution represents twoway (bidirectional) flow, which may be written in the following equations.
˙m0,z1=Cdθz=z1∫z=0ρν(z)Wdz
˙mz1,H=Cdθz=H∫z=z1ρν(z)Wdz
The two real solutions represent threeway flow, which may be written in the following equations.
˙mz2,H=Cdθz=H∫z=z2ρν(z)Wdz
˙mz1,z2=Cdθz=z2∫z=z1ρν(z)Wdz
˙mz2,H=Cdθz=H∫z=z2ρν(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 horizontallypivoted window is shown in Figure 2.
The opening angle α (090°) is linearly proportional to the window opening factor (01.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(α))
h4=AxisHeight+(WindowHeight−AxisHeight)cos(α)
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:
Wpivot=
⎷11W2+1(2(AxisHeight−z)tan(α))2
The mass flow rate in the pivoted area becomes:
˙mpivot=Cdθz=h4∫z=h2ρν(z)Wpivotdz
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 horizontallypivoted window. The actual window width and height are used to calculate airflows for a horizontallypivoted window.
The schematic drawing of air flow patterns through a simple vertical opening (AirflowNetwork:Multizone:Component:SimpleOpening) is shown in Figure 3. 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(y)=P0n−ρngy
Pm(y)=P0m−ρmgy
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)=Cd√2Pn(y)−Pm(y)ρ
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=Pon−P0mg(ρn−ρm)orPom−P0ng(ρm−ρn)
When the neutral plane is within the opening (first pattern in Figure 3), twoway (bidirectional) flows occur. The total flow through a large opening is the sum of both flows.
˙m0,Y=Cdθy=Y∫y=0ρν(y)Wdy
˙mY,H=Cdθy=H∫y=Yρν(y)Wdy
When the neutral plane is below or above the large opening (second and third pattern in Figure 3), oneway flow occurs.
˙m=Cdθy=H∫y=0ρν(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 twoway 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 threeway 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 oneway flow is possible using the detailed and simple opening components which are meant for vertical or nearvertical openings. In reality, there are twoway 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 twoway 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 4. Forced and buoyancy airflows are described separately below.
Forced airflows[LINK]
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:
Scenario 1 PL=PU˙mU=˙mL=0
where:
PL

= Air pressure in the lower zone [Pa]

PU

= Air pressure in the upper zone [Pa]

˙mU

= Air mass flow rate from the lower zone to the upper zone driven by forced airflow pressure difference [kg/s]

˙mL

= Air mass flow rate from the upper zone to the lower zone driven by forced airflow pressure difference [kg/s]

Scenario 2
PL>PU
˙mU=0
˙mL=ρLACd(2ΔPρave)0.5
∂˙mL∂PL=ρLACd(12ΔPρave)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]

Scenario 3
PL<PU
˙mL=0
˙mU=−ρUACd(2ΔPρave)0.5
∂˙mU∂PL=ρUACd(12ΔPρave)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 bidirectional 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:
˙mbuo,max=ρave⋅0.055(gΔρD5Hρave)0.5
where:
˙mbuo,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 coexist, it is possible to have either unidirectional or bidirectional flows. For example, when the upward force due to pressure difference is greater than the buoyancy force (downward), unidirectional flow occurs. Bidirectional flow only occurs when the upward imposed force is less than the buoyancy force. The critical pressure between unidirectional and bidirectional 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):
ΔPFlood=∣∣
∣∣C2ShapegΔρD5H2A2∣∣
∣∣
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]

CShape={0.754foracircleopening0.942(w/D)forarectangleopening}
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 bidirectional flow due to the buoyancy force. When the pressure difference increases from 0 and is less than ΔPFlood, there is some bidirectional flow across the opening, but less than the maximum flow. If the pressure difference keeps increasing and exceeds ΔPFlood, there is no bidirectional flow. Cooper’s model assumes the buoyancy flow varies linearly with pressure difference.
˙mbuo=⎧⎨⎩˙mbuo,max⋅(1−ΔPΔPFlood)IfΔρ>0andΔPΔPFlood<10Otherwise⎫⎬⎭
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. PL=PU
˙mU=˙mbuo
˙mL=˙mbuo
∂˙mL∂PL=0
b. PL>PU
˙mU=˙mbuo
˙mL=ρLACd(2ΔPρave)0.5+˙mbuo
∂˙mL∂PL=ρLACd(12ΔPρave)0.5−˙mbuo,maxΔPFlood
where:
ρL

= Air density in the lower zone [kg/m3]

A

= Opening area [m2]

Cd

= Discharge coefficient [dimensionless]

Pave

= Average air density between the lower and upper zones [Pa]

ΔP

= Pressure difference PL  PU [Pa]

c. PL<PU
˙mL=˙mbuo
˙mU=−ρUACd(2ΔPρave)0.5+˙mbuo
∂˙mU∂PL=ρUACd(12ΔPρave)0.5+˙mbuo,maxΔPFlood
Sloping Plane[LINK]
When a staircase is introduced as shown in Figure 5, 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):
Aeff=Asinα(1+cosα)
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 6 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 6.
Legend Description
Forced downward 
Forced flow rate from upper to lower at P P < 0 
Forced upward 
Forced flow rate from lower to upper at P P > 0 
Buoyancy upward 
Total upward flow rate due to buoyancy only at P P < 0 
Buoyancy downward 
Total downward flow rate due to buoyancy only at P P > 0 
Combined downward 
Total downward flow at P P < 0 (Forced downward + buoyancy upward) 
Combined upward 
Total upward flow at P P > 0 (Forced upward + buoyancy downward) 
Wind pressure calculations[LINK]
The wind pressure is determined by Bernoulli’s equation, assuming no height change or pressure losses:
pw=CpρV2ref2
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)
Vref=Vmet(δmetzmet)αmet(zδ)α
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 = “AVERAGESURFACE 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=0.6⋅ln[1.248−0.703sin(α/2)−1.175sin2(α)+0.131sin3(2αG)+0.769cos(α/2)+0.07G2sin2(α/2)+0.717cos2(α/2)]
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 30degree increments

For walls of a high rise building, a twodimensional array of surfaceaveraged 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 twodimensional array of surfaceaveraged 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:
pw,n=Cp,nρV2ref2
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:
{P}∗={P}−{C}
where the correction vector, {C}, is computed by the matrix relationship:
[J]{C}={B}
{B} is a column vector with each component given by:
Bn=∑i˙mi
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:
Jn,m=∑i∂˙m∂Pm
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:
Pn∗=Pn−Cn/(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 fixedpoint 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:
˙mCpdTdx=UP(T∞−T)
where:
Cp

= Specific heat of airflow [J/kgK]

˙m

= Airflow rate [kg/s]

P

= Perimeter of a duct element [m]

T

= Temperature as a field variable [°C]

T∞

= Temperature of air surrounding the duct element [°C]

U

= Overall heat transfer coefficient [W/m2K]

U=11hi+1ho+∑tjkj
hi

= Inside heat transfer coefficient [W/m2K]

ho

= Outside heat transfer coefficient [W/m2K]

tj

= Thickness at jth layer [m]

kj

= Thermal conductivity at jth layer [W/mK]

The outlet air temperature at the end of the duct (x = L) is:
To=T∞+(Ti−T∞)⋅exp(−UA˙mCp)
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=˙mCp(T∞−Ti)[1−exp(−UA˙mCp)]
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.
Duct Radiation[LINK]
The ductsurface userdefined view factors object defines the view factor from the duct to the zone surfaces. Because the view factors, duct surface area, and temperatures are known, the heat gain to the duct can be directly calculated and applied to the duct air heat balance to determine the outlet temperature. The duct heat gain is then treated as a direct heat loss from each participating surface. Once the heat loss to the duct for each surface is calculated, the surface temperatures are adjusted to maintain the surface energy balance.
Radiation is treated as point loads heat source/sink on the zone surfaces and the duct. These point loads will be added or subtracted using the same implementation strategy as the convective load for the ducts to the zone air. The radiative rate of energy transfer between the duct (d) and a given surface (s) for two grey surfaces can be written as shown below.
˙Qd→s=σ(T4d−T4s)1−ϵdfeAdϵd+1feAdFd→s+1−ϵsAsϵs
where:
σ

= SteffanBoltzmann constant [W/m2K4]

Td

= Duct surface temperature [C]

Ts

= Zone surface temperature [C]

ϵd

= Duct grey body emissivity []

ϵs

= Zone surface grey body emissivity []

Ad

= Duct surface area [m2]

As

= Zone surface area [m2]

fe

= Duct exposure fraction []

Fd→s

= Ducttosurface view factor

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:
˙mdWdx=UmP(W∞−W)
where:
˙m

= Airflow rate [kg/s]

P

= Perimeter of a duct element [m]

W

= Humidity ratio [kg/kg]

W∞

= Humidity ratio of air surrounding the duct element [kg/kg]

Um

= Overall moisture transfer coefficient [kg/m2s]

Um=11hm,i+1hm,o+∑tjDj
hm,i

= Inside moisture transfer coefficient [kg/m2s]

hm,o

= Outside moisture transfer coefficient [kg/m2s]

tj

= Thickness at jth layer [m]

Dj

= Moisture diffusivity at jth layer [kg/ms]

The outlet air humidity ratio at the end of the duct (x = L) is:
Wo=W∞+(Wi−W∞)⋅exp(−UmA˙m)
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
Qm=˙m⋅(W∞−Wi)[1−exp(−UmA˙m)]
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:
{Mm}[W]=[Bm]
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 forcedair 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:
MCPairflow=˙minfCp+∑(˙mmixCp)
MCPTairflow=˙minfCpTamb+∑(˙mmixCpTzone)
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]

˙minf

= Incoming air mass flow rate from outdoors [kg/s]

˙mmix

= Incoming air mass flow rate from adjacent zones [kg/s]

Tamb

= Outdoor air drybulb temperature [°C]

Tzone

= Adjacent zone air temperature [°C]

The latent load items from multizone load calculations may be written as follows:
Mairflow=˙minf+∑˙mmix
MWairflow=˙minfWamb+∑˙mmixWzone
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]

˙minf

= Incoming air mass flow rate from outdoors [kg/s]

˙mmix

= 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=∑jQcond(i,j)+∑jQleak(i,j)
QADS,m,i=∑jQcond,m(i,j)+∑jQleak,m(i,j)
where:
QADS,i

= Total sensible load in the ith zone due to ADS losses [W]

Qcond(ij)

= Duct wall conduction loss at the jth duct located in the ith zone [W]

Qleak(ij)

= Sensible supply leak loss at the jth linkage located in the ith zone [W]

QADS,m,i

= Total latent load in the ith zone due to ADS losses [kg/s]

Qcond,m(ij)

= Duct wall vapor diffusion loss at the jth duct located in the ith zone [kg/s]

Qleak,m(ij)

= Latent supply leak loss at the jth linkage located in the ith 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 partload 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 partload 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 Partload ratio + Airflow during ADS off cycle (1.0 – ADS Partload 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
˙mfan=∑ni˙mi,terminal1−∑kjFj
where:
˙mfan

= The supply fan flow rate

˙mi,terminal

= The flow rate at the ith terminal, which is determined in the subroutine SimVAV in the HVACSingleDuctSystem module

n

= Number of terminals

Fj

= 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=˙mfan,max˙mfan,cal
˙mi,terminal,final=˙mi,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

˙mfan,max

= The maximum supply fan flow rate by input

˙mfan,cal

= The calculated supply fan flow rate

˙mi,terminal,final

= The final flow rate at each terminal adjusted by the ratio

The inputs of both AirflowNetwork:IntraZone:Node and AirflowNetwork:IntraZone:Linkage objects are treated as normal model nodes and linkages. Therefore, the calculation procedures are the same as above.
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:
Ttz=∑Nsli=1˙Qi+∑Nsurfacesi−1hiAiTsi+MCPTairflow+˙msysCpTsupply+QADSz−(Czδt)(−3Tt−δtz+32Tt−2δtz−13Tt−3δtz)116Czδt+∑Nsurfacesi=1hiA+MCPairflow+˙msysC
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=∑kgmass,sched−loads+MWairflow+˙msys,inWsys+surfs∑i=1AihmiρairzWsurfsi+QADSm,z
where MWairflow is the sum of mass flow rate multiplied by humidity ratio for infiltration and mixing and QADS,m,zis the added 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 PredictorCorrector” found elsewhere in this document.
The available outputs from the AirflowNetwork model are described in the EnergyPlus Input Output Reference manual.
Single Sided Wind Pressure Coefficient Algorithm[LINK]
The approach in the single sided (SS) model is based on correlating the flow rate into the room through one or more openings with local flow conditions at the openings. These are in one façade of a rectangular building of width WB, depth DB and height HB. There is a vertical wind profile U(z) approaching the building, which generates a local flow as=0 parallel to the façade (whose value depends on location, building geometry and wind direction).
Consider first the 2aperture case, where air flows in through one opening of area Q′ and exits through another of equal area. Let us denote the flow rate of ambient air into the room by Q, and define a dimensionless flow rate Δp(t) by Q′=[apΔcp+aσσΔcp]1/2+asU′L where ¯¯¯¯¯¯¯Δp is a reference velocity taken as the undisturbed wind speed σΔp at a reference height ¯¯¯¯¯¯¯Δp. There is a pressure difference ¯¯¯¯¯¯¯Δp1/2 between the two opening locations, with timeaveraged mean value as and standard deviation aσ.
There are three principal mechanisms determining the flow rate:
Mean pressure difference: a nonzero mean value for the difference in pressure at the two opening locations, as, drives a flow that scales as ap.[ssmeanmech]
Unsteady pressure difference: turbulence drives a fluctuating pressure difference between the two openings, resulting in a net addition of air from the outside. This is characterized by the standard deviation of the pressure difference signal, aσ and the associated flow rate adds to that arising from as. It includes periodic Strouhal forcing (‘pumping’) at relevant wind angles. Note that this contribution is nonzero even when the mean pressure difference is zero.[ssunsteadymech]
Shear layer mechanism: when there is small or zero mean pressure difference between the two opening locations, the local flow parallel to the façade causes a mixing layer to grow over the length of the opening resulting in the incorporation of ambient air. This scales with as=0.[ssshearmech]
Our model therefore views the flow rate as comprising contributions from each of these as appropriate, i.e. θ=φw−φn. where s′=0.32, A1 and Δcp are constants, the mean and fluctuating pressure difference terms have been written as pressure coefficients Q′=(apΔcp+aσσΔcp)1/2 Δcp=f(s′,θ) ρ being the ambient air density and σΔcp=g(s′,θ)
The pressuredifference mechanisms [ssmeanmech] and [ssunsteadymech] are dominant for the 2opening case and [ssshearmech] is negligible, while, in contrast, in the case of a single opening, shear is dominant. Hence we take σΔcp for 2opening cases and ap for one opening. Henceforth we focus on the 2opening case, although a procedure is given at the end of the section for how to approximate a greater number by just two.
Figure 7 shows the geometry of this case in plan view. The two openings are a distance s apart in a façade of width WB whose normal has azimuthal angle aσ. The wind azimuthal direction is φw, so that its direction relative to the façade is f(s′,θ)=s′⋅Π(θ), In these definitions, Δcp, while σΔcp. Note that if θ initially falls outside the range indicated it should be replaced by Q.
According to the above discussion, the flow rate is modeled as comprising two terms, arising, respectively, from the mean and fluctuating components of the pressure difference between the two aperture locations: g(s′,θ)=Σ0+s′⋅Σ(θ).
There are two parts to deriving the detailed form of the model:
Express the mean and fluctuating pressures in terms of other physical quantities that are more readily accessible. EnergyPlus does not have the information on variations in pressure that is needed to compute Δp and Δcp>0. Analysis of wind tunnel pressure data (Linden et al., 2013) allows for the derivation of simple expressions for these quantities in terms of the wind direction and opening separation s.[ssitemthefirst]
Find the coefficients s′=0.32 and A1. This is achieved here by analysis of ventilation rate data (Linden et al., 2013), and using linear regression on the data.[ssitemthesecond]
Estimation of pressure difference coefficients[LINK]
The analysis of pressure data carried out showed that the contributions to the righthand side of Equation [eqn:ssflowtwoterms] can be written to first order as Σ(θ)=0.423−1.63×10−3θ and Σ0=0.24 where the dimensionless separation aσ, QAinUref=(0.01+[0.173Π(θ)+0.042Σ(θ)]s′)1/2 and Δp1212ρV2=1C2d(Q/AeffV)2 That is, to a first approximation the pressure difference mean and standard deviation vary linearly with opening separation, each modulated by its own function of relative wind angle θ.
The formulae for the components are as follows. Aeff=A1A2(A21+A22)1/2 ΔcE+p,12=1C2d(AinAeff)2(UrefV)2Q′2 pw=ρa[Vref(zwin]2Cp/2 Vref(z)=Umet(δmetz)αmet(zδ)α
Figures 8 and 9 show an example of Δp and Δcp>0 generated by Equations [eqn:ssfirstformula] through [eqn:sslastformula], respectively, compared with the data.
This completes item [ssitemthefirst] above, i.e. the approximation of the two pressure coefficient terms in Equation [eqn:ssflowtwoterms] in terms of the dimensionless opening separation and the relative wind angle. There remains item [ssitemthesecond], i.e. fitting the flow rate data to the formulagenerated pressure parameters to determine the constants s′=0.32 and A1. When this was carried out, the best fit was found to occur with A2 and Cd. Figure 10 shows the comparison between the flow rate prediction and the experimental data.
Combining with the previous results we have ΔcE+p,12=1C2d[Vref(10)Vref(zwin)]2(0.02+[0.346Π(θ)+0.084Σ(θ)]⋅s′)
The flow Q is defined to be positive when opening #1 is the inlet and opening #2 is the outlet. Since Δp12 is defined as Aeff, this can occur when ΔcE+p,12, which in turn corresponds to certain relative wind angles, i.e. Vref(10)Vref(zwin)=(10zwin)α
Note that the sense of Q is not necessarily welldefined for all such angles, e.g. when the unsteady contribution in Equation [eqn:ssflowtwoterms] is significant: nevertheless, we retain this definition for all angles.
Inversion of flow rate to equivalent pressure difference[LINK]
Finally, we invert the flow rate to obtain the equivalent pressure difference. The purpose of this procedure is to allow incorporation of the model into the network airflow part of EnergyPlus, which requires pressure coefficients to calculate the flow rate: in this way the model will become essentially transparent to the user. A nodal representation of the 2aperture singlesided case is shown in Figure 11.
Here, 0 is the ambient, 1 is at the upstream outer surface, 2 is at the downstream outer surface, and R is the room.
Let the openings have areas Uref and (Ain/Aeff)2=2, not necessarily equal, and a discharge coefficient A1=A2=Ain, which is assumed equal, and let Q be the flow rate through the system. Then it can be shown that the pressure difference between the two openings, Aeff=Ain/√2, is related to the flow rate by θ=−0.0028(100−PPD)2+0.3419(100−PPD)−6.6275 where the effective area (Uref/V) is defined by Aeff=A1A2(A21+A22)1/2 and V is a reference velocity. The lefthand side of Equation [eqn:ssdeltap12] is a pressure coefficient, pw, say, and our task is to express this in terms of the flow rate in Equation [eqn:sscombinedeqn] but using the definition of a pressure coefficient appropriate to EnergyPlus.
The first step is to rewrite Equation [eqn:ssdeltap12] as follows: ΔcE+p,12=1C2d(AinAeff)2(UrefV)2Q′2 where Vref(zwin) is the nondimensional flow rate of Equation [eqn:sscombinedeqn], which is expressed in terms of Q, the window area V=Vref(zwin) and the wind speed at 10 m height, Umet.
The factor zmet=10 m, since δmet=270 m, and therefore αmet=0.14. The factor Vref(10) relating the velocity scale used in the correlation with that used in the EnergyPlus pressure coefficient may be simplified by noting the following.
EnergyPlus assumes the surface pressure due to wind, +0.5ΔcE+p,12, is derived from a pressure coefficient Cp using the local wind speed at window height, −0.5ΔcE+p,12, i.e. pw=ρa[Vref(zwin]2Cp/2 Hence the reference velocity ΔcE+p,12.
EnergyPlus characterizes the vertical wind profile as a power law in z, defined in terms of the local atmospheric boundary layer depth δ and a power law exponent α. EnergyPlus also uses an “interpolation” formula to derive the above local wind speed from that measured at the meteorological site, s′=0.75, according to Vref(z)=Umet(δmetz)αmet(zδ)α where the subscript “met” denotes conditions at the meteorological site, and the default values are s′=0.32, δmet=270 m and αmet=0.14. The user can change the site parameters from their default values.
The correlation in Equation [eqn:sscombinedeqn] assumes the reference velocity is the wind speed at 10 m, which will be taken as the local velocity at the building, Vref(10). Combining these results, the pressure coefficient representing the difference in pressure between the two openings is given by ΔcE+p,12=1C2d[Vref(10)Vref(zwin)]2(0.02+[0.346Π(θ)+0.084Σ(θ)]⋅s′) where the velocity ratio Vref(10)Vref(zwin)=(10zwin)α
Thus, if we assign a pressure coefficient +0.5ΔcE+p,12 to opening 1 and −0.5ΔcE+p,12 to opening 2, then this will provide the necessary pressure difference to give the flow in Equation [eqn:sscombinedeqn] in a network context.
Note: the pressure coefficient pw refers to the difference in pressure between the two openings and is defined in terms of the wind speed at window height. If this is to be combined with a background pressure coefficient for the façade it is important to ensure the two are defined using the same reference velocity.
Testing/Validation/Data Source(s)[LINK]
The SS model was developed using wind tunnel pressure and flow rate data acquired by CPP Wind Engineering & Air Quality Consultants, who carried out two series of tests: (a) Closed Box tests, with a sealed building, to obtain pressure and velocity data and (b) Ventilation tests, with one or two apertures open, to obtain flow rate data. See Linden et al. (2013) for full details.
Each experimental configuration consisted of a cuboidal building with or without a set of blocks representing nearby buildings.
Building: either 2story, with WB = 47.5 cm, DB = 19.3 cm, HB = 9.9 cm, or 4story, with HB = 19.8 cm. The model has a nominal scale of 1:70 (based on earlier tests).
Environment: building either isolated or at center of a 3 x 3 array of equalsized blocks, which were of same aspect ratio and either low (2story height) or high (4story height), and with either narrow spacing (approximately 0.5WB between blocks) or wide spacing (approximately WB), giving a total of 5 different environments.
Approach flow: suburban boundary layer, approximately 10 m/s at 1 m height. The model was mounted on a turntable allowing any desired wind direction to be used. Increments of s′=0.75 and s′=0.32 were used for the Closed Box and Ventilation runs, respectively.
Room and openings: the building envelope contained a secondfloor room of dimensions WRM = 45.9 cm, DRM = 17.7 cm and HRM = 5.0 cm. A total of 15 opening positions were available, distributed around the room. These were sealed for the Closed Box tests but opened in various combinations for the Ventilation runs.
Sensors: 56 pressure transducers, arranged around the opening positions, including 32 at midopening height, recorded static pressure, while hot film sensors at 6 different positions recorded flow speed near the surface. For Ventilation runs the decay in concentration of a tracer was measured at two locations inside the room.
The data used to develop the model consisted of 6 scenarios: 3 building/environment cases,
2story isolated building
4story isolated building
2story building with low, widelyspaced blocks
combined with 2 choices of opening position on the South (long) façade, namely wide separation (s′=0.75) and medium separation (s′=0.32).
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.
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 drybulb temperature Tout. The comfort temperature calculation may be based on two curves and a boundary temperature point.
Tcomf={LowTempCurve,minTout<x≤BoundaryPoi
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.
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:
˙mi=CiρΔPiμ
where:
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=∑˙mi∑˙mi
Absolute airflow tolerance=∣∣∑˙mi∣∣
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:
ΔP=(Pn+ρV2n2)−(Pm+ρV2m2)+ρg(zn−zm)
where:
By rearranging terms and adding wind pressure impacts, the above equation may be rewritten in the format used by the airflow network model:
ΔP=Pn−Pm+PS+PW
where:
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 [fig:thegeneralproblemofgravitationalflow]
The schematic drawing of a possible air flow pattern through a detailed vertical opening (AirflowNetwork:Multizone:Component:DetailedOpening) is shown in Figure 1. The equations used below are extracted from the COMIS Fundamentals manual (1990).
The air density is assumed to be a linear function of height:
ρi(z)=ρ0i+biz
The pressure difference is assumed to be linear and simulate the effect of turbulence:
ΔPt=Pt0+btz
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:
P1(z)−P2(z)=(P01−P02)−g[(ρ01z+b1z2/2)−(ρ02z+b2z2/2)]+(Pt0+btz)
The velocity at any level z is given by
v(z)=√2P1(z)−P2(z)ρ
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(b1−b2)z2/2+[g(ρ01−ρ02)−bt]z+(−P01+P02−Pt0)=0
This equation above can have two, one, or zero real solutions. The zero solution represents a oneway flow through the opening and may be expressed in the following equation:
˙m=Cdθz=H∫z=0ρν(z)Wdz
The one real solution represents twoway (bidirectional) flow, which may be written in the following equations.
˙m0,z1=Cdθz=z1∫z=0ρν(z)Wdz
˙mz1,H=Cdθz=H∫z=z1ρν(z)Wdz
The two real solutions represent threeway flow, which may be written in the following equations.
˙mz2,H=Cdθz=H∫z=z2ρν(z)Wdz
˙mz1,z2=Cdθz=z2∫z=z1ρν(z)Wdz
˙mz2,H=Cdθz=H∫z=z2ρν(z)Wdz
where:
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 horizontallypivoted window is shown in Figure 2.
Schematic drawing of a horizontallypivoted window [fig:schematicdrawingofahorizontallypivoted]
The opening angle α (090°) is linearly proportional to the window opening factor (01.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(α))
h4=AxisHeight+(WindowHeight−AxisHeight)cos(α)
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:
Wpivot= ⎷11W2+1(2(AxisHeight−z)tan(α))2
The mass flow rate in the pivoted area becomes:
˙mpivot=Cdθz=h4∫z=h2ρν(z)Wpivotdz
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 horizontallypivoted window. The actual window width and height are used to calculate airflows for a horizontallypivoted window.
The schematic drawing of air flow patterns through a simple vertical opening (AirflowNetwork:Multizone:Component:SimpleOpening) is shown in Figure 3. The equations used below are available from Walton (1989).
Schematic of large opening and associated three flow patterns [fig:schematicoflargeopeningandassociated]
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(y)=P0n−ρngy
Pm(y)=P0m−ρmgy
where:
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)=Cd√2Pn(y)−Pm(y)ρ
where:
The neutral height, Y, where the velocity of the air is zero, may be calculated in the following equation:
Y=Pon−P0mg(ρn−ρm)orPom−P0ng(ρm−ρn)
When the neutral plane is within the opening (first pattern in Figure 3), twoway (bidirectional) flows occur. The total flow through a large opening is the sum of both flows.
˙m0,Y=Cdθy=Y∫y=0ρν(y)Wdy
˙mY,H=Cdθy=H∫y=Yρν(y)Wdy
When the neutral plane is below or above the large opening (second and third pattern in Figure 3), oneway flow occurs.
˙m=Cdθy=H∫y=0ρν(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 twoway 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 threeway 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 oneway flow is possible using the detailed and simple opening components which are meant for vertical or nearvertical openings. In reality, there are twoway 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 twoway 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 4. Forced and buoyancy airflows are described separately below.
Air movements across a horizontal opening [fig:airmovementsacrossahorizontalopening]
Forced airflows[LINK]
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:
Scenario 1 PL=PU˙mU=˙mL=0
where:
Scenario 2
PL>PU
˙mU=0
˙mL=ρLACd(2ΔPρave)0.5
∂˙mL∂PL=ρLACd(12ΔPρave)0.5
where:
Scenario 3
PL<PU
˙mL=0
˙mU=−ρUACd(2ΔPρave)0.5
∂˙mU∂PL=ρUACd(12ΔPρave)0.5
where:
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 bidirectional 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:
˙mbuo,max=ρave⋅0.055(gΔρD5Hρave)0.5
where:
Combined airflows[LINK]
When forced and buoyancy flows coexist, it is possible to have either unidirectional or bidirectional flows. For example, when the upward force due to pressure difference is greater than the buoyancy force (downward), unidirectional flow occurs. Bidirectional flow only occurs when the upward imposed force is less than the buoyancy force. The critical pressure between unidirectional and bidirectional 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):
ΔPFlood=∣∣ ∣∣C2ShapegΔρD5H2A2∣∣ ∣∣
where:
CShape={0.754foracircleopening0.942(w/D)forarectangleopening}
where:
As mentioned above, when the air pressure difference between two zones is zero there is the maximum bidirectional flow due to the buoyancy force. When the pressure difference increases from 0 and is less than ΔPFlood, there is some bidirectional flow across the opening, but less than the maximum flow. If the pressure difference keeps increasing and exceeds ΔPFlood, there is no bidirectional flow. Cooper’s model assumes the buoyancy flow varies linearly with pressure difference.
˙mbuo=⎧⎨⎩˙mbuo,max⋅(1−ΔPΔPFlood)IfΔρ>0andΔPΔPFlood<10Otherwise⎫⎬⎭
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. PL=PU
˙mU=˙mbuo
˙mL=˙mbuo
∂˙mL∂PL=0
b. PL>PU
˙mU=˙mbuo
˙mL=ρLACd(2ΔPρave)0.5+˙mbuo
∂˙mL∂PL=ρLACd(12ΔPρave)0.5−˙mbuo,maxΔPFlood
where:
c. PL<PU
˙mL=˙mbuo
˙mU=−ρUACd(2ΔPρave)0.5+˙mbuo
∂˙mU∂PL=ρUACd(12ΔPρave)0.5+˙mbuo,maxΔPFlood
Sloping Plane[LINK]
A Staircase is attached to the horizontal opening. [fig:astaircaseisattachedtothehorizontal]
When a staircase is introduced as shown in Figure 5, 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):
Aeff=Asinα(1+cosα)
where:
Note: the hydraulic diameter calculation is based on the effective opening area, while the opening depth remains the same.
Figure 6 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 6.
Flow rates at different pressure differences [fig:flowratesatdifferentpressuredifferences]
Wind pressure calculations[LINK]
The wind pressure is determined by Bernoulli’s equation, assuming no height change or pressure losses:
pw=CpρV2ref2
where:
Vref may be expressed as (Ref, Local Wind Speed Calculations)
Vref=Vmet(δmetzmet)αmet(zδ)α
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 = “AVERAGESURFACE 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=0.6⋅ln[1.248−0.703sin(α/2)−1.175sin2(α)+0.131sin3(2αG)+0.769cos(α/2)+0.07G2sin2(α/2)+0.717cos2(α/2)]
where:
For walls of a high rise building, a twodimensional array of surfaceaveraged 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 twodimensional array of surfaceaveraged 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:
pw,n=Cp,nρV2ref2
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:
{P}∗={P}−{C}
where the correction vector, {C}, is computed by the matrix relationship:
[J]{C}={B}
{B} is a column vector with each component given by:
Bn=∑i˙mi
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:
Jn,m=∑i∂˙m∂Pm
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:
Pn∗=Pn−Cn/(1−r)
where:
This method is similar to a Steffensen iteration (Conte and de Boor 1972) which is used as a fixedpoint 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]
General[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:
˙mCpdTdx=UP(T∞−T)
where:
U=11hi+1ho+∑tjkj
The outlet air temperature at the end of the duct (x = L) is:
To=T∞+(Ti−T∞)⋅exp(−UA˙mCp)
where:
The heat transfer by convection to ambient, Q, is:
Q=˙mCp(T∞−Ti)[1−exp(−UA˙mCp)]
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:
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.
Duct Radiation[LINK]
The ductsurface userdefined view factors object defines the view factor from the duct to the zone surfaces. Because the view factors, duct surface area, and temperatures are known, the heat gain to the duct can be directly calculated and applied to the duct air heat balance to determine the outlet temperature. The duct heat gain is then treated as a direct heat loss from each participating surface. Once the heat loss to the duct for each surface is calculated, the surface temperatures are adjusted to maintain the surface energy balance.
Radiation is treated as point loads heat source/sink on the zone surfaces and the duct. These point loads will be added or subtracted using the same implementation strategy as the convective load for the ducts to the zone air. The radiative rate of energy transfer between the duct (d) and a given surface (s) for two grey surfaces can be written as shown below.
˙Qd→s=σ(T4d−T4s)1−ϵdfeAdϵd+1feAdFd→s+1−ϵsAsϵs
where:
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:
˙mdWdx=UmP(W∞−W)
where:
Um=11hm,i+1hm,o+∑tjDj
The outlet air humidity ratio at the end of the duct (x = L) is:
Wo=W∞+(Wi−W∞)⋅exp(−UmA˙m)
where:
The moisture transfer by convection to ambient, Qm, is
Qm=˙m⋅(W∞−Wi)[1−exp(−UmA˙m)]
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:
{Mm}[W]=[Bm]
where:
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 forcedair 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:
MCPairflow=˙minfCp+∑(˙mmixCp)
MCPTairflow=˙minfCpTamb+∑(˙mmixCpTzone)
where:
The latent load items from multizone load calculations may be written as follows:
Mairflow=˙minf+∑˙mmix
MWairflow=˙minfWamb+∑˙mmixWzone
where:
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=∑jQcond(i,j)+∑jQleak(i,j)
QADS,m,i=∑jQcond,m(i,j)+∑jQleak,m(i,j)
where:
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 partload 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 partload 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
where:
Zone load
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
˙mfan=∑ni˙mi,terminal1−∑kjFj
where:
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=˙mfan,max˙mfan,cal
˙mi,terminal,final=˙mi,terminal⋅R
where:
Airflow Calculation Procedure using inputs of Intrazone nodes and linkages[LINK]
The inputs of both AirflowNetwork:IntraZone:Node and AirflowNetwork:IntraZone:Linkage objects are treated as normal model nodes and linkages. Therefore, the calculation procedures are the same as above.
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:
Ttz=∑Nsli=1˙Qi+∑Nsurfacesi−1hiAiTsi+MCPTairflow+˙msysCpTsupply+QADSz−(Czδt)(−3Tt−δtz+32Tt−2δtz−13Tt−3δtz)116Czδt+∑Nsurfacesi=1hiA+MCPairflow+˙msysC
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=∑kgmass,sched−loads+MWairflow+˙msys,inWsys+surfs∑i=1AihmiρairzWsurfsi+QADSm,z
where MWairflow is the sum of mass flow rate multiplied by humidity ratio for infiltration and mixing and QADS,m,zis the added 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 PredictorCorrector” found elsewhere in this document.
Model Output[LINK]
The available outputs from the AirflowNetwork model are described in the EnergyPlus Input Output Reference manual.
Single Sided Wind Pressure Coefficient Algorithm[LINK]
The approach in the single sided (SS) model is based on correlating the flow rate into the room through one or more openings with local flow conditions at the openings. These are in one façade of a rectangular building of width WB, depth DB and height HB. There is a vertical wind profile U(z) approaching the building, which generates a local flow as=0 parallel to the façade (whose value depends on location, building geometry and wind direction).
Consider first the 2aperture case, where air flows in through one opening of area Q′ and exits through another of equal area. Let us denote the flow rate of ambient air into the room by Q, and define a dimensionless flow rate Δp(t) by Q′=[apΔcp+aσσΔcp]1/2+asU′L where ¯¯¯¯¯¯¯Δp is a reference velocity taken as the undisturbed wind speed σΔp at a reference height ¯¯¯¯¯¯¯Δp. There is a pressure difference ¯¯¯¯¯¯¯Δp1/2 between the two opening locations, with timeaveraged mean value as and standard deviation aσ.
There are three principal mechanisms determining the flow rate:
Mean pressure difference: a nonzero mean value for the difference in pressure at the two opening locations, as, drives a flow that scales as ap.[ssmeanmech]
Unsteady pressure difference: turbulence drives a fluctuating pressure difference between the two openings, resulting in a net addition of air from the outside. This is characterized by the standard deviation of the pressure difference signal, aσ and the associated flow rate adds to that arising from as. It includes periodic Strouhal forcing (‘pumping’) at relevant wind angles. Note that this contribution is nonzero even when the mean pressure difference is zero.[ssunsteadymech]
Shear layer mechanism: when there is small or zero mean pressure difference between the two opening locations, the local flow parallel to the façade causes a mixing layer to grow over the length of the opening resulting in the incorporation of ambient air. This scales with as=0.[ssshearmech]
Our model therefore views the flow rate as comprising contributions from each of these as appropriate, i.e. θ=φw−φn. where s′=0.32, A1 and Δcp are constants, the mean and fluctuating pressure difference terms have been written as pressure coefficients Q′=(apΔcp+aσσΔcp)1/2 Δcp=f(s′,θ) ρ being the ambient air density and σΔcp=g(s′,θ)
The pressuredifference mechanisms [ssmeanmech] and [ssunsteadymech] are dominant for the 2opening case and [ssshearmech] is negligible, while, in contrast, in the case of a single opening, shear is dominant. Hence we take σΔcp for 2opening cases and ap for one opening. Henceforth we focus on the 2opening case, although a procedure is given at the end of the section for how to approximate a greater number by just two.
Figure 7 shows the geometry of this case in plan view. The two openings are a distance s apart in a façade of width WB whose normal has azimuthal angle aσ. The wind azimuthal direction is φw, so that its direction relative to the façade is f(s′,θ)=s′⋅Π(θ), In these definitions, Δcp, while σΔcp. Note that if θ initially falls outside the range indicated it should be replaced by Q.
According to the above discussion, the flow rate is modeled as comprising two terms, arising, respectively, from the mean and fluctuating components of the pressure difference between the two aperture locations: g(s′,θ)=Σ0+s′⋅Σ(θ).
Plan view of building with 2opening façade. [fig:ssplanviewtwoopenings]
There are two parts to deriving the detailed form of the model:
Express the mean and fluctuating pressures in terms of other physical quantities that are more readily accessible. EnergyPlus does not have the information on variations in pressure that is needed to compute Δp and Δcp>0. Analysis of wind tunnel pressure data (Linden et al., 2013) allows for the derivation of simple expressions for these quantities in terms of the wind direction and opening separation s.[ssitemthefirst]
Find the coefficients s′=0.32 and A1. This is achieved here by analysis of ventilation rate data (Linden et al., 2013), and using linear regression on the data.[ssitemthesecond]
Estimation of pressure difference coefficients[LINK]
The analysis of pressure data carried out showed that the contributions to the righthand side of Equation [eqn:ssflowtwoterms] can be written to first order as Σ(θ)=0.423−1.63×10−3θ and Σ0=0.24 where the dimensionless separation aσ, QAinUref=(0.01+[0.173Π(θ)+0.042Σ(θ)]s′)1/2 and Δp1212ρV2=1C2d(Q/AeffV)2 That is, to a first approximation the pressure difference mean and standard deviation vary linearly with opening separation, each modulated by its own function of relative wind angle θ.
The formulae for the components are as follows. Aeff=A1A2(A21+A22)1/2 ΔcE+p,12=1C2d(AinAeff)2(UrefV)2Q′2 pw=ρa[Vref(zwin]2Cp/2 Vref(z)=Umet(δmetz)αmet(zδ)α
Figures 8 and 9 show an example of Δp and Δcp>0 generated by Equations [eqn:ssfirstformula] through [eqn:sslastformula], respectively, compared with the data.
Mean pressure difference, Δp predicted by Equations [eqn:ssfirstformula] and [eqn:sssecondformula], compared with experimental data. The façade containing the openings is facing South and the separation is large (s′=0.75). [fig:ssfigure2upper]
Standard deviation of pressure difference, Δcp>0, predicted by Equations [eqn:ssthirdformula] and [eqn:sslastformula], compared with experimental data. The façade containing the openings is facing South and the separation is large (s′=0.75). [fig:ssfigure2lower]
This completes item [ssitemthefirst] above, i.e. the approximation of the two pressure coefficient terms in Equation [eqn:ssflowtwoterms] in terms of the dimensionless opening separation and the relative wind angle. There remains item [ssitemthesecond], i.e. fitting the flow rate data to the formulagenerated pressure parameters to determine the constants s′=0.32 and A1. When this was carried out, the best fit was found to occur with A2 and Cd. Figure 10 shows the comparison between the flow rate prediction and the experimental data.
Combining with the previous results we have ΔcE+p,12=1C2d[Vref(10)Vref(zwin)]2(0.02+[0.346Π(θ)+0.084Σ(θ)]⋅s′)
The flow Q is defined to be positive when opening #1 is the inlet and opening #2 is the outlet. Since Δp12 is defined as Aeff, this can occur when ΔcE+p,12, which in turn corresponds to certain relative wind angles, i.e. Vref(10)Vref(zwin)=(10zwin)α
Note that the sense of Q is not necessarily welldefined for all such angles, e.g. when the unsteady contribution in Equation [eqn:ssflowtwoterms] is significant: nevertheless, we retain this definition for all angles.
Comparison of predicted flow rates with experimental data. Each graph shows experimental data (red diamonds) plotted against model prediction (gray squares). Cases in the lefthand column, (a)(c), are for Q′=Q/(AinUref), those in righthand column, (d)(f), Ain. The three rows are for 2story isolated, 4story isolated and 2story with low density surroundings, respectively. [fig:ssfigure3inoriginal]
Inversion of flow rate to equivalent pressure difference[LINK]
Finally, we invert the flow rate to obtain the equivalent pressure difference. The purpose of this procedure is to allow incorporation of the model into the network airflow part of EnergyPlus, which requires pressure coefficients to calculate the flow rate: in this way the model will become essentially transparent to the user. A nodal representation of the 2aperture singlesided case is shown in Figure 11.
Pressure network for 2opening singlesided room. [fig:sspressurenetworktwoopenings]
Here, 0 is the ambient, 1 is at the upstream outer surface, 2 is at the downstream outer surface, and R is the room.
Let the openings have areas Uref and (Ain/Aeff)2=2, not necessarily equal, and a discharge coefficient A1=A2=Ain, which is assumed equal, and let Q be the flow rate through the system. Then it can be shown that the pressure difference between the two openings, Aeff=Ain/√2, is related to the flow rate by θ=−0.0028(100−PPD)2+0.3419(100−PPD)−6.6275 where the effective area (Uref/V) is defined by Aeff=A1A2(A21+A22)1/2 and V is a reference velocity. The lefthand side of Equation [eqn:ssdeltap12] is a pressure coefficient, pw, say, and our task is to express this in terms of the flow rate in Equation [eqn:sscombinedeqn] but using the definition of a pressure coefficient appropriate to EnergyPlus.
The first step is to rewrite Equation [eqn:ssdeltap12] as follows: ΔcE+p,12=1C2d(AinAeff)2(UrefV)2Q′2 where Vref(zwin) is the nondimensional flow rate of Equation [eqn:sscombinedeqn], which is expressed in terms of Q, the window area V=Vref(zwin) and the wind speed at 10 m height, Umet.
The factor zmet=10 m, since δmet=270 m, and therefore αmet=0.14. The factor Vref(10) relating the velocity scale used in the correlation with that used in the EnergyPlus pressure coefficient may be simplified by noting the following.
EnergyPlus assumes the surface pressure due to wind, +0.5ΔcE+p,12, is derived from a pressure coefficient Cp using the local wind speed at window height, −0.5ΔcE+p,12, i.e. pw=ρa[Vref(zwin]2Cp/2 Hence the reference velocity ΔcE+p,12.
EnergyPlus characterizes the vertical wind profile as a power law in z, defined in terms of the local atmospheric boundary layer depth δ and a power law exponent α. EnergyPlus also uses an “interpolation” formula to derive the above local wind speed from that measured at the meteorological site, s′=0.75, according to Vref(z)=Umet(δmetz)αmet(zδ)α where the subscript “met” denotes conditions at the meteorological site, and the default values are s′=0.32, δmet=270 m and αmet=0.14. The user can change the site parameters from their default values.
The correlation in Equation [eqn:sscombinedeqn] assumes the reference velocity is the wind speed at 10 m, which will be taken as the local velocity at the building, Vref(10). Combining these results, the pressure coefficient representing the difference in pressure between the two openings is given by ΔcE+p,12=1C2d[Vref(10)Vref(zwin)]2(0.02+[0.346Π(θ)+0.084Σ(θ)]⋅s′) where the velocity ratio Vref(10)Vref(zwin)=(10zwin)α
Thus, if we assign a pressure coefficient +0.5ΔcE+p,12 to opening 1 and −0.5ΔcE+p,12 to opening 2, then this will provide the necessary pressure difference to give the flow in Equation [eqn:sscombinedeqn] in a network context.
Note: the pressure coefficient pw refers to the difference in pressure between the two openings and is defined in terms of the wind speed at window height. If this is to be combined with a background pressure coefficient for the façade it is important to ensure the two are defined using the same reference velocity.
Local velocity as function of wind angle for wind tunnel data analyzed. On the 2story graph ‘WP85’ refers to Warren and Parkins (1985), whose measurement location is most closely matched by the “Center” curve. [fig:sslocalvelocity]
Wind tunnel data, averaged over all cases, plotted with curve to fit data. [fig:sswindtunneldata]
Testing/Validation/Data Source(s)[LINK]
The SS model was developed using wind tunnel pressure and flow rate data acquired by CPP Wind Engineering & Air Quality Consultants, who carried out two series of tests: (a) Closed Box tests, with a sealed building, to obtain pressure and velocity data and (b) Ventilation tests, with one or two apertures open, to obtain flow rate data. See Linden et al. (2013) for full details.
Each experimental configuration consisted of a cuboidal building with or without a set of blocks representing nearby buildings.
Building: either 2story, with WB = 47.5 cm, DB = 19.3 cm, HB = 9.9 cm, or 4story, with HB = 19.8 cm. The model has a nominal scale of 1:70 (based on earlier tests).
Environment: building either isolated or at center of a 3 x 3 array of equalsized blocks, which were of same aspect ratio and either low (2story height) or high (4story height), and with either narrow spacing (approximately 0.5WB between blocks) or wide spacing (approximately WB), giving a total of 5 different environments.
Approach flow: suburban boundary layer, approximately 10 m/s at 1 m height. The model was mounted on a turntable allowing any desired wind direction to be used. Increments of s′=0.75 and s′=0.32 were used for the Closed Box and Ventilation runs, respectively.
Room and openings: the building envelope contained a secondfloor room of dimensions WRM = 45.9 cm, DRM = 17.7 cm and HRM = 5.0 cm. A total of 15 opening positions were available, distributed around the room. These were sealed for the Closed Box tests but opened in various combinations for the Ventilation runs.
Sensors: 56 pressure transducers, arranged around the opening positions, including 32 at midopening height, recorded static pressure, while hot film sensors at 6 different positions recorded flow speed near the surface. For Ventilation runs the decay in concentration of a tracer was measured at two locations inside the room.
The data used to develop the model consisted of 6 scenarios: 3 building/environment cases,
2story isolated building
4story isolated building
2story building with low, widelyspaced blocks
combined with 2 choices of opening position on the South (long) façade, namely wide separation (s′=0.75) and medium separation (s′=0.32).
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.
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 drybulb temperature Tout. The comfort temperature calculation may be based on two curves and a boundary temperature point.
Tcomf={LowTempCurve,minTout<x≤BoundaryPoi