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. The impact of thermal capacity will be addressed in future upgrades to this model.
Pressure and Airflow Calculations[LINK]
The EnergyPlus airflow network consists of a set of nodes linked by airflow components. Therefore, it is a simplified airflow model, compared to detailed models such as those used in computational fluid dynamics (CFD) models. The node variable is pressure and the linkage variable is airflow rate. A brief description is presented below. A detailed description of the airflow network model may be found in the work of Walton (1989), Dols and Walton (2002), and Walton and Dols (2003).
Initialization[LINK]
Newton’s method is used to solve for node air pressures and it requires an initial set of values for the node pressures. There are two initialization methods available. The first is linear initialization and equivalent to Initialization flag = 0. These initial values may be obtained by including in each airflow component a linear approximation relating airflow to pressure drop:
where
= Air mass flow rate at i-th linkage [kg/s]
Ci = Air mass flow coefficient [m3]
i= Pressure difference across the i-th linkage [Pa]
*µ= Air viscosity [Pa-*s]
This initialization handles stack effects very well and tends to establish the proper direction for the airflows. The linear approximation is provided by the laminar regime.
The second initialization method assumes the initial pressures are zero and uses Newton’s method directly.
Convergence criteria[LINK]
Conservation of air mass flow rate at each linkage provides the convergence criterion. When the sum of mass flow rates in all the linkages approaches zero within the convergence tolerance, the solution has converged. The solution is assumed to have converged when the sum is less than the convergence value, in order to reduce the number of iterations and obtain sufficient accuracy. There are two convergence criteria used in the AirflowNetwork model: Relative airflow convergence tolerance and Absolute airflow convergence tolerance.
Relative airflow tolerance = 
Absolute airflow tolerance = 
The relative airflow tolerance is equivalent to the ratio of the absolute value of the sum of all network airflows to the sum of the network airflow magnitudes. The absolute airflow tolerance is the summation of the absolute value of all network airflows. The solution has converged when both of these convergence criteria have been met.
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:
where
ΔP= Total pressure difference between nodes n and m [Pa]
Pn, Pm= Entry and exit static pressures [Pa]
Vn, Vm= Entry and exit airflow velocities [m/s]
*ρ= Air density [kg/m*3]
g= Acceleration due to gravity [9.81 m/s2]
zn, zm= Entry and exit elevations [m]
By rearranging terms and adding wind pressure impacts, the above equation may be rewritten in the format used by the airflow network model:
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 135. The equations used below are extracted from the COMIS Fundamentals manual (1990).
The air density is assumed to be a linear function of height:
The pressure difference is assumed to be linear and simulate the effect of turbulence:
The reference pressures on each side are given at the bottom of the opening. By assuming the Bernoulli hypothesis on both sides, the pressure difference can be defined at any level of z as:
The velocity at any level z is given by
The locations of the two possible neutral planes are given by an equilibrium in pressure which leads to a zero velocity point. By assuming the left terms in the equation above to be zero, one may have:
This equation above can have two, one, or zero real solutions. The zero solution represents a one-way flow through the opening and may be expressed in the following equation:
The one real solution represents two-way (bi-directional) flow, which may be written in the following equations.
The two real solutions represent three-way flow, which may be written in the following equations.
where
Cd= discharge coefficient [dimensionless]
θ= Area reduction factor [dimensionless]
W= Opening width [m]
The discharge coefficient, opening width, opening height, and start height factor are modulated based on opening factors. A detailed description of opening factor calculations may be found in the Input Output Reference (Ref. AirflowNetwork:Multizone:Zone, AirflowNetwork:Multizone:Surface, and
AirflowNetwork:Multizone:Component:DetailedOpening).
The above calculation procedure is used for a normal rectangular window. For a horizontally pivoted rectangular window, the calculation procedure is slightly different. A schematic drawing of a horizontally-pivoted window is shown in Figure 136.
The opening angle α (0-90°) is linearly proportional to the window opening factor (0-1.0). An opening factor of 1.0 is equal to an opening angle of 90°. The heights in the pivoted area are expressed as:
When z < h2 or z > h4, where z is the distance from the bottom of the window, the integration procedure is the same as the procedure for a normal rectangular window. When h2 PU
= 0
where:
ρL= Air density in the lower zone [kg/m3]
A= Opening area [m2]
Cd= Discharge coefficient [Dimensionless]
ρave = Average air density between the lower and upper zones [Pa]
ΔP= Pressure difference PL -PU [Pa]
PL < PU
= 0
where:
ρU= Air density in the upper zone [kg/m3]
A= Opening area [m2]
Cd= Discharge coefficient [Dimensionless]
ρave = Average air density between the lower and upper zones [Pa]
ΔP= Pressure difference PL -PU [Pa]
Buoyancy airflows[LINK]
Buoyancy flow only occurs when the air density in the upper zone is greater than the air density in the lower zone. The flow is bi-directional and the amount of upper flow is equal to the lower flow across the opening. The following discussion assumes the air density in the upper zone is greater than the air density in the lower zone. Otherwise, the buoyancy flow rate is equal to zero. It is also assumed that the maximum buoyancy flow occurs when the pressure difference across the opening due to forced airflows is zero. The maximum buoyancy flow may be expressed as a part of Cooper’s model:
where:
= Buoyancy mass flow rate at zero forced airflow pressure difference [kg/s]
g= Gravity acceleration [m/s2]
DH= Hydraulic diameter [m]
ρave = Average air density between the lower and upper zones [kg/m3]
Δρ= Density difference between the lower and upper zones [kg/m3]
Combined airflows[LINK]
When forced and buoyancy flows co-exist, it is possible to have either unidirectional or bi-directional flows. For example, when the upward force due to pressure difference is greater than the buoyancy force (downward), unidirectional flow occurs. Bi-directional flow only occurs when the upward imposed force is less than the buoyancy force. The critical pressure between unidirectional and bi-directional flows is called the purge pressure (Tan and Jaluria 1992). The purge pressure is a function of opening geometry and the buoyancy force (ΔP/(gΔρDH)) and may be expressed as (Cooper 1998):
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]
where:
w= Opening width [m]
D= Opening depth [m]
As mentioned above, when the air pressure difference between two zones is zero there is the maximum bi-directional flow due to the buoyancy force. When the pressure difference increases from 0 and is less than |ΔPFlood|, there is some bi-directional flow across the opening, but less than the maximum flow. If the pressure difference keeps increasing and exceeds |ΔPFlood|, there is no bi-directional flow. Cooper’s model assumes the buoyancy flow varies linearly with pressure difference.
The total air flow across the opening is based on superposition of the forced and buoyancy flows, and may be expressed for three different pressure difference scenarios as follows:
PL = PU
PL > PU
where:
ρL= Air density in the lower zone [kg/m3]
A= Opening area [m2]
Cd= Discharge coefficient [dimensionless]
Ρave = Average air density between the lower and upper zones [Pa]
ΔP= Pressure difference PL -PU [Pa]
PL < PU
Sloping plane
When a staircase is introduced as shown in Figure 139, 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):
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 140 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 140.
Legend Description
|
Forced downward
|
Forced flow rate from upper to lower at PL-PU < 0
|
|
Forced upward
|
Forced flow rate from lower to upper at PL-PU > 0
|
|
Buoyancy upward
|
Total upward flow rate due to buoyancy only at PL-PU < 0
|
|
Buoyancy downward
|
Total downward flow rate due to buoyancy only at PL-PU > 0
|
|
Combined downward
|
Total downward flow at PL-PU < 0 (Forced downward + buoyancy upward)
|
|
Combined upward
|
Total upward flow at PL-PU > 0 (Forced upward + buoyancy downward)
|
: Legend Description
Wind pressure calculations[LINK]
The wind pressure is determined by Bernoulli’s equation, assuming no height change or pressure losses:
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):
Cp is a function of location on the building envelope and wind direction. When Wind Pressure Coefficient Type = “INPUT”, the Cp values are explicitly defined in the input for AirflowNetwork:Multizone:Wind Pressure Coefficient Values. When Wind Pressure Coefficient Type = “AVERAGE-SURFACE CALCULATION” and the building shape is rectangular, the program uses the following equations to calculate wind pressure coefficient (Cp) values for different wind directions. For a low rise building, the normalized surface pressure coefficient may be written as (Swami and Chandra 1988):
where
Cp,n= Cp value at a given angle between wind direction and the outward normal of the surface under consideration [dimensionless]
α= Angle between wind direction and outward normal of wall under consideration [deg]
G= Natural log of the ratio of the width of the wall under consideration to the width of the adjacent wall [dimensionless]
n = Index of incident angle at 30-degree increments
For walls of a high rise building, a two-dimensional array of surface-averaged wind pressure coefficients is generated based on wind incident angle and side ratio. The wind pressure coefficients are provided in 2001 ASHRAE Fundamentals Handbook, p. 16.5, Fig. 7, “Surface Averaged Wall Pressure Coefficients for Tall Buildings”. The original work was performed by Atkins et al. (1979). The incident angle has an increment of 30 degrees. The side ratio values are 0.25, 1.0, and 4.0. For a given incident angle and building aspect ratio, the program uses linear interpolation to calculate the corresponding wind pressure coefficient Cp,n.
For the roof of a high rise building, a two-dimensional array of surface-averaged wind pressure coefficients is also generated based on wind incident angle and side ratio. The wind pressure coefficients are provided in 2001 ASHRAE Fundamentals Handbook, p. 16.6, Fig. 9, “Surface Averaged Roof Pressure Coefficients for Tall Buildings”. The original work was performed by Holmes (1986). The incident angle has an increment of 30 degrees. The side ratio values are 0.25, 1.0, and 4.0. At a given wind incident angle and building aspect ratio, the program uses linear interpolation to calculate the corresponding wind pressure coefficient Cp,n.
The wind surface pressure at the given incident angle can be calculated by combining the above two equations:
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:
where the correction vector, {C}, is computed by the matrix relationship:
{B} is a column vector with each component given by:
where n is the node number and i indicates all flow paths connecting node n to other nodes, and [J] is the square Jacobian matrix whose elements are given by:
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:
where
r= the ratio of Cn for the current iteration to its value for the previous iteration [dimensionless]
Cn= Correction value at the nth node [Pa]
Pn~~= Estimated pressure at the nth node [Pa]
Pn*~~= Corrected pressure at the nth node used in the next iteration [Pa]
This method is similar to a Steffensen iteration (Conte and de Boor 1972) which is used as a fixed-point iteration method for individual nonlinear equations.
The iteration correction method presented in the above equation gives a variable factor. When the solution is close to convergence, the solution method converges quadratically. By limiting cases where the value of r is less than some value, such as -0.5, the solution will not interfere with the rapid convergence. It has not been proven that the convergence acceleration equation will always lead to convergence, but it can be shown that it will not prevent convergence. Newton’s method converges when the estimated solution values are within some distance, called the radius of convergence, or the correct solution. Applying the convergence acceleration equation when -1 < r <0, will cause a smaller correction than Newton’s method, which therefore, can not force the iterations outside the radius of convergence. When r<-1, the solution diverges in an oscillatory fashion. When r>1, the solution also diverges, but in a nonoscillatory manner. For 0<r<1, the solution is approached from one direction. In all three cases, the convergence acceleration equation applies as long as r is truly constant over several iterations. However, for the last case, this involves a true extrapolation of correction factor which is very sensitive to the accuracy of r. This is most extreme for the case of r=1, which would cause an infinite correction.
Node Temperature Calculations[LINK]
A brief description of the air node temperature calculation is given below. A detailed description can be found in the work of Swami et al. (1992). The following equation is used to calculate temperature distribution across a duct element at the given airflow rate and inlet air temperature:
where
Cp= Specific heat of airflow [J/kg•K]
= Airflow rate [kg/s]
P= Perimeter of a duct element [m]
T= Temperature as a field variable [°C]
= Temperature of air surrounding the duct element [°C]
U= Overall heat transfer coefficient [W/m2•K]
hi= Inside heat transfer coefficient [W/m2•K]
ho= Outside heat transfer coefficient [W/m2•K]
tj= Thickness at j-th layer [m]
kj= Thermal conductivity at j-th layer [W/m•K]
The outlet air temperature at the end of the duct (x=L) is:
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:
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:
where
{M}= Airflow matrix
[T]= Temperature vector
[B]= Given boundary conditions
The zone air temperatures and primary air loop component (fan and coils) outlet conditions are used as prescribed conditions in the AirflowNetwork model. In addition, the temperature difference across zone loop components (terminal units) is held constant during the calculations. For example, thermal zone temperatures calculated during the previous system time step are used as prescribed temperatures when calculating all other node temperatures. The zone air temperature is assumed constant (prescribed) throughout the AirflowNetwork iterative solution. The fan and coil outlet air temperatures, and terminal unit temperature differences are assumed constant within an AirflowNetwork iteration. The sensible heat gains calculated during the AirflowNetwork solution are then used to predict a new zone air temperature.
Node Humidity Ratio Calculations[LINK]
A brief description of the air node humidity ratio calculation is given below. A detailed description can found in the work of Swami et al. (1992). The following equation is used to calculate humidity ratio distribution across a duct element at the given airflow rate and inlet air humidity ratio:
where
= Airflow rate [kg/s]
P= Perimeter of a duct element [m]
W= Humidity ratio [kg/kg]
= Humidity ratio of air surrounding the duct element [kg/kg]
Um= Overall moisture transfer coefficient [kg/m2•s]
hm,i~~= Inside moisture transfer coefficient [kg/m2•s]
*hm,o= Outside moisture transfer coefficient [kg/m*2•s]
tj= Thickness at j-th layer [m]
Dj= Moisture diffusivity at j-th layer [kg/m•s]
The outlet air humidity ratio at the end of the duct (x=L) is:
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
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:
where
{Mm}= Airflow matrix
[W]= Humidity ratio vector
[Bm]= Given boundary conditions
The zone air humidity ratios and primary air loop component (fan and coils) outlet conditions are used as prescribed conditions in the AirflowNetwork model. For example, thermal zone humidity ratios calculated during the previous system time step are used as prescribed humidity ratios when calculating all other node humidity ratios. The zone air humidity ratio is assumed constant (prescribed) throughout the AirflowNetwork iterative solution. The coil outlet air humidity ratio is assumed constant within an AirflowNetwork iteration. The latent heat gains calculated during the AirflowNetwork solution are then used to predict a new zone air humidity ratio.
Sensible and Latent Load Calculations[LINK]
The zone sensible and latent loads calculated in the AirflowNetwork model consist of multizone, duct conduction and leakage. In addition, the impact of infiltration and mixing is accounted for in this calculation. The multizone load only includes incoming airflows from outside (infiltration) and other adjacent zones (mixing) with and without forced-air fan operation. It is divided into two terms: variable and constant. The constant term is the sum of the mass flow rate multiplied by the specific heat for both infiltration and mixing. The variable term includes the impact of zone and outdoor air temperature. Each of these terms is used in the zone energy balance equation. The sensible load items from the multizone load calculations may be written as follows:
where
MCPairflow = Sum of air mass flow rate multiplied by specific heat for infiltration and mixing [W/K]
MCPTairflow = Sum of air mass flow rate multiplied by specific heat and temperature for infiltration and mixing [W]
= Incoming air mass flow rate from outdoors [kg/s]
= Incoming air mass flow rate from adjacent zones [kg/s]
Tamb= Outdoor air dry-bulb temperature [°C]
Tzone= Adjacent zone air temperature [°C]
The latent load items from multizone load calculations may be written as follows:
where
Mairflow = Sum of air mass flow rates for infiltration and mixing [kg/s]
MWairflow = Sum of air mass flow rate multiplied by humidity ratio for infiltration and mixing [kg/s]
= Incoming air mass flow rate from outdoors [kg/s]
= Incoming air mass flow rate from adjacent zones [kg/s]
Wamb= Outdoor air humidity ratio [kg/kg]
Wzone= Adjacent zone air humidity ratio [kg/kg]
The air distribution system (ADS) loads due to duct conduction and leakage depend on the air distribution system component (e.g., duct) location. The air distribution system sensible and latent loads are calculated for each zone as follows:
where
QADS,i= Total sensible load in the i-th zone due to ADS losses [W]
Qcond(ij)= Duct wall conduction loss at the j-th duct located in the i-th zone [W]
Qleak(ij)= Sensible supply leak loss at the j-th linkage located in the i-th zone [W]
QADS,m,i= Total latent load in the i-th zone due to ADS losses [kg/s]
Qcond,m(ij)~~= Duct wall vapor diffusion loss at the j-th duct located in the i-th zone [kg/s]
Qleak,m(ij) = Latent supply leak loss at the j-th linkage located in the i-th zone [kg/s]
Impacts of Supply Air Constant Volume Fan Control on Load: Cycling vs. Continuous[LINK]
The AirflowNetwork model currently allows two types of constant volume fans: Fan:ConstantVolume and Fan:OnOff. The Fan:ConstantVolume object has only one type of supply air fan operation mode: continuous fan, cycling compressor (ContinuousFanWithCyclingCompressor). However, the Fan:OnOff has two types of supply air fan operation modes: cycling fan, cycling compressor (CyclingFanAndCompressor) or continuous fan, cycling compressor (ContinuousFanWithCyclingCompressor). The CyclingFanAndCompressor operation mode is frequently referred to as “AUTO fan”, where the compressor(s) and supply air fan operate in unison to meet the zone heating or cooling load, and cycle off together when the heating or cooling load has been met. The ContinuousFanWithCyclingCompressor operation mode is often referred to as “fan ON”, where the compressor(s) cycle on and off to meet the zone heating or cooling load but the supply air fan operates continuously regardless of compressor operation. The supply air fan operation mode is specified in an HVAC system object based on a given fan operation mode schedule (e.g., AirLoopHVAC:UnitaryHeatCool object).
The determination of the zone sensible and latent loads caused by multizone airflows and forced air distribution system operation is dependent on the supply air fan operation mode (see Sensible and Latent Load Calculations section above). The zone loads calculated by the AirflowNetwork model are added to zone sensible and latent balances in the ZonePredictorCorrector module to calculate zone air temperatures and humidity ratios (see Integration of the AirflowNetwork Model section below). For the case of continuous fan/cycling compressor, the zone loads during forced air distribution system operation are calculated with the system design air mass flow rate without modification, since the system air node conditions (temperature and humidity) reflect the average values for the time step considering the coil/fan on and off periods during the time step.
For the case of cycling fan/cycling compressor, where the forced air distribution system can operate for a portion of the simulation time step, the airflows are determined based on the air distribution system part-load ratio (ratio of the average air mass flow rate for the time step divided by the design air mass flow rate). The airflows for the AirflowNetwork:Distribution:Linkage objects are reported during the air distribution system on cycle, since no airflow is assumed during the system off cycle. The airflows for the AirflowNetwork:Multizone:Surface objects are weighted by the air distribution system part-load ratio. The zone loads are the sum of energy losses during both the air distribution system on and off cycle at each system time step. The energy losses when the air distribution system is on are calculated using the system “on” air flow rate multiplied by the air distribution system run time fraction. The energy losses when the air distribution system is off are obtained from the multizone airflow calculations (without forced air through the air distribution system) and these losses are multiplied by (1.0 - system run time fraction), assuming no airflows through the air distribution system when the fan is off. The formulas used to calculate airflows and zone loads are given below:
Airflow
Airflow = Airflow during ADS on cycle * ADS Part-load ratio + Airflow during ADS off cycle * (1.0 – ADS Part-load ratio)
where ADS = Air Distribution System
Zone load
System run time fraction = Maximum run time fraction of coils and fans in the air distribution system
Zone energy losses = Zone energy loss during ADS on cycle * System run time fraction + Zone energy loss during ADS off cycle * (1.0 - System run time fraction)
The calculation of loads due to multizone airflow, without forced air distribution system operation, is performed when the HVAC system is first simulated during a simulation time step (FirstHVACIteration = True). The calculation of loads due to air distribution system operation is performed on subsequent iterations within the same simulation time step when the mass flow rate at the supply air fan inlet is greater than 0.0 (and FirstHVACIteration = False).
Airflow Calculation Procedure using A Supply Variable Air Volume Fan[LINK]
The AirflowNetwork model currently also allows a variable air volume fan type as Fan:VariableVolume. The allowed terminal type is AirTerminal:SingleDuct:VAV:Reheat only. Other types of terminals will be added later.
In general, the supply fan air flow rate in a VAV central system is determined by a sum of terminal flow rates when the AirflowNetwork model is not applied. When the AirflowNetwork model is applied and the supply air fan flow rate is given, each terminal flow is determined by pressure resistance of each supply air pathway. It is possible that the delivered air flow rate from the pressure resistance at each terminal may be totally different from the desired flow rate determined by terminal units. Therefore, it is not easy to meet both requirements. The following compromised approach, including possible supply and return leaks in an air distribution system, is implemented.
Set up terminal airflows in the AirflowNetwork module based on the SimVAV subroutine in the HVACSingleDuctSystem module.
Require AirflowNetwork:Distribution:Component:LeakageRatio objects to define supply leaks, so that the values of the Effective Leakage Ratio field are used to decide the supply fan flow rates. The base of the ratio will be actual supply fan flow rate, instead of the maximum fan flow rate used in the constant volume fan systems.
Assign the supply fan airflow rate based on the sum of all terminal flow rates and all supply leak ratios until it reaches the maximum fan flow rate
where
= The supply fan flow rate
= The flow rate at the ith terminal, which is determined in the subroutine SimVAV in the HVACSingleDuctSystem module
n= Number of terminals
= 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. ****
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
= The maximum supply fan flow rate by input
= The calculated supply fan flow rate
= The final flow rate at each terminal adjusted by the ratio
Integration of the AirflowNetwork Model[LINK]
The loads calculated by the AirflowNetwork model are integrated into the EnergyPlus heat balance equation in a similar manner as described elsewhere in this document in the section “Basis for the Zone and System Integration”. The mass flow rate summations and sensible and latent loads described in the previous section are included in the calculation of zone temperature and humidity ratio.
The revised zone temperature update equation becomes:
Where MCPTairflow is the sum of mass flow rate multiplied by specific heat and temperature for infiltration and mixing, QADS,z is the added total sensible load in the zone due to Air Distribution System losses, and MCPairflow is the sum of mass flow rate multiplied by specific heat for infiltration and mixing as calculated in the AirflowNetwork model described above.
The revised coefficient (B) used in the zone humidity ratio calculation is shown below:
Where MWairflow is the sum of mass flow rate multiplied by humidity ratio for infiltration and mixing and QADS,m,z~~is 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 Predictor-Corrector” found elsewhere in this document.
The available outputs from the AirflowNetwork model are described in the EnergyPlus Input Output Reference manual.
Atkins, R. E., J. A. Peterka, and J. E. Cermak. 1979. “Averaged pressure coefficients for rectangular buildings,” Wind Engineering, Proceedings of the Fifth International Conference 7:369-80, Fort Collins, CO. Pergamon Press, NY.
Bolmqvist, C. and M. Sandberg, 2004, “Air Movements through Horizontal Openings in Buildings – A Model Study,” International Journal of Ventilation, Vol. 3, No. 1, pp. 1-9
COMIS Fundamentals. 1990. Edited by Helmut E. Feustel and Alison Rayner-Hooson, LBL-28560, Lawrence Berkeley Laboratory, Berkeley, CA
Conte, S. D. and C de Boor. 1972. Elementary Numerical Analysis: an Algorithmic Approach, McGraw-Hill.
Cooper, L., 1989, “Calculation of the Flow Through a Horizontal Ceiling/Floor Vent,” NISTIR 89-4052, National Institute of Standards and Technology, Gaithersburg, MD
Dols, W. S. & G. N. Walton. 2002. “CONTAMW 2.0 User Manual,” NISTIR 6921, National Institute of Standards and Technology, Gaithersburg, Maryland
Holmes, J. D. 1986. Wind Loads on low-rise buildings: The structural and environmental effects of wind on buildings and structures, Chapter 12, Faculty of Engineering, Monash University, Melbourne, Australia
Swami, M. V. and S. Chandra. 1988. Correlations for pressure distribution on buildings and calculation of natural-ventilation airflow, ASHRAE Transactions 94(1988) (Pt 1), pp. 243-266.
Swami, M. V., L. Gu, & V. Vasanth. 1992. “Integration of Radon and Energy Models for Building,” FSEC-CR-553-92, Florida Solar Energy Center, Cocoa, Florida
Tan, Q. and Y. Jaluria, 1992, “Flow through Horizontal Vents as Related to Compartment Fire Environments,” NIST-GCR-92-607, National Institute of Standards and Technology, Gaithersburg, Maryland
Walton, G. N. 1989. “AIRNET – A Computer Program for Building Airflow Network Modeling,” NISTIR 89-4072, National Institute of Standards and Technology, Gaithersburg, Maryland
Walton, G. N. & W. S. Dols. 2003. “CONTAM 2.1 Supplemental User Guide and Program Documentation,” NISTIR 7049, National Institute of Standards and Technology, Gaithersburg, Maryland
AirflowNetwork Model[LINK]
Overview[LINK]
The AirflowNetwork model provides the ability to simulate the performance of an air distribution system, including supply and return leaks, and calculate multizone airflows driven by outdoor wind and forced air during HVAC system operation. The pressure and airflow model described here was developed based on AIRNET (Walton 1989). This detailed model is used to simulate thermal conduction and air leakage losses for constant volume air distribution systems (e.g., in residential or light commercial buildings). The multizone airflow calculations are performed at the HVAC system time step which, among other benefits,.allows for modeling hybrid ventilation systems.
Model Description[LINK]
The input object AirflowNetwork:SimulationControl provides access to the airflow network method, which consists of a set of nodes connected by airflow components through linkages. The objects AirflowNetwork:Multizone:Zone, AirflowNetwork:Multizone:ExternalNode, and AirflowNetwork:Distribution:Node represent airflow nodes. The objects AirflowNetwork:Multizone:Surface and AirflowNetwork:Distribution:Linkage represent airflow linkages. The other objects with a relationship between pressure and airflow represent airflow components.
The AirflowNetwork model consists of three sequential steps:
The pressure and airflow calculations determine pressure at each node and airflow through each linkage given wind pressures and forced airflows. Based on the airflow calculated for each linkage, the model then calculates node temperatures and humidity ratios given zone air temperatures and zone humidity ratios. Using these node temperatures and humidity ratios, the sensible and latent loads from duct system conduction and leakage are summed for each zone. The sensible and latent loads obtained in this step are then used in the zone energy balance equations to predict HVAC system loads and to calculate the final zone air temperatures, humidity ratios, and pressures.
The present AirflowNetwork model may only be applied to a single heating and cooling system that uses a single air distribution system (a single AirLoopHVAC object). The model excludes the impact of the air and duct system thermal capacitance. The impact of thermal capacity will be addressed in future upgrades to this model.
Pressure and Airflow Calculations[LINK]
The EnergyPlus airflow network consists of a set of nodes linked by airflow components. Therefore, it is a simplified airflow model, compared to detailed models such as those used in computational fluid dynamics (CFD) models. The node variable is pressure and the linkage variable is airflow rate. A brief description is presented below. A detailed description of the airflow network model may be found in the work of Walton (1989), Dols and Walton (2002), and Walton and Dols (2003).
Initialization[LINK]
Newton’s method is used to solve for node air pressures and it requires an initial set of values for the node pressures. There are two initialization methods available. The first is linear initialization and equivalent to Initialization flag = 0. These initial values may be obtained by including in each airflow component a linear approximation relating airflow to pressure drop:
where
Ci = Air mass flow coefficient [m3]
*µ= Air viscosity [Pa-*s]
This initialization handles stack effects very well and tends to establish the proper direction for the airflows. The linear approximation is provided by the laminar regime.
The second initialization method assumes the initial pressures are zero and uses Newton’s method directly.
Convergence criteria[LINK]
Conservation of air mass flow rate at each linkage provides the convergence criterion. When the sum of mass flow rates in all the linkages approaches zero within the convergence tolerance, the solution has converged. The solution is assumed to have converged when the sum is less than the convergence value, in order to reduce the number of iterations and obtain sufficient accuracy. There are two convergence criteria used in the AirflowNetwork model: Relative airflow convergence tolerance and Absolute airflow convergence tolerance.
Relative airflow tolerance =
Absolute airflow tolerance =
The relative airflow tolerance is equivalent to the ratio of the absolute value of the sum of all network airflows to the sum of the network airflow magnitudes. The absolute airflow tolerance is the summation of the absolute value of all network airflows. The solution has converged when both of these convergence criteria have been met.
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:
where
ΔP= Total pressure difference between nodes n and m [Pa]
Pn, Pm= Entry and exit static pressures [Pa]
Vn, Vm= Entry and exit airflow velocities [m/s]
*ρ= Air density [kg/m*3]
g= Acceleration due to gravity [9.81 m/s2]
zn, zm= Entry and exit elevations [m]
By rearranging terms and adding wind pressure impacts, the above equation may be rewritten in the format used by the airflow network model:
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 135. The equations used below are extracted from the COMIS Fundamentals manual (1990).
The air density is assumed to be a linear function of height:
The pressure difference is assumed to be linear and simulate the effect of turbulence:
The reference pressures on each side are given at the bottom of the opening. By assuming the Bernoulli hypothesis on both sides, the pressure difference can be defined at any level of z as:
The velocity at any level z is given by
The locations of the two possible neutral planes are given by an equilibrium in pressure which leads to a zero velocity point. By assuming the left terms in the equation above to be zero, one may have:
This equation above can have two, one, or zero real solutions. The zero solution represents a one-way flow through the opening and may be expressed in the following equation:
The one real solution represents two-way (bi-directional) flow, which may be written in the following equations.
The two real solutions represent three-way flow, which may be written in the following equations.
where
Cd= discharge coefficient [dimensionless]
θ= Area reduction factor [dimensionless]
W= Opening width [m]
The discharge coefficient, opening width, opening height, and start height factor are modulated based on opening factors. A detailed description of opening factor calculations may be found in the Input Output Reference (Ref. AirflowNetwork:Multizone:Zone, AirflowNetwork:Multizone:Surface, and
AirflowNetwork:Multizone:Component:DetailedOpening).
The above calculation procedure is used for a normal rectangular window. For a horizontally pivoted rectangular window, the calculation procedure is slightly different. A schematic drawing of a horizontally-pivoted window is shown in Figure 136.
The opening angle α (0-90°) is linearly proportional to the window opening factor (0-1.0). An opening factor of 1.0 is equal to an opening angle of 90°. The heights in the pivoted area are expressed as:
When z < h2 or z > h4, where z is the distance from the bottom of the window, the integration procedure is the same as the procedure for a normal rectangular window. When h2 PU
where:
ρL= Air density in the lower zone [kg/m3]
A= Opening area [m2]
Cd= Discharge coefficient [Dimensionless]
ρave = Average air density between the lower and upper zones [Pa]
ΔP= Pressure difference PL -PU [Pa]
PL < PU
where:
ρU= Air density in the upper zone [kg/m3]
A= Opening area [m2]
Cd= Discharge coefficient [Dimensionless]
ρave = Average air density between the lower and upper zones [Pa]
ΔP= Pressure difference PL -PU [Pa]
Buoyancy airflows[LINK]
Buoyancy flow only occurs when the air density in the upper zone is greater than the air density in the lower zone. The flow is bi-directional and the amount of upper flow is equal to the lower flow across the opening. The following discussion assumes the air density in the upper zone is greater than the air density in the lower zone. Otherwise, the buoyancy flow rate is equal to zero. It is also assumed that the maximum buoyancy flow occurs when the pressure difference across the opening due to forced airflows is zero. The maximum buoyancy flow may be expressed as a part of Cooper’s model:
where:
g= Gravity acceleration [m/s2]
DH= Hydraulic diameter [m]
ρave = Average air density between the lower and upper zones [kg/m3]
Δρ= Density difference between the lower and upper zones [kg/m3]
Combined airflows[LINK]
When forced and buoyancy flows co-exist, it is possible to have either unidirectional or bi-directional flows. For example, when the upward force due to pressure difference is greater than the buoyancy force (downward), unidirectional flow occurs. Bi-directional flow only occurs when the upward imposed force is less than the buoyancy force. The critical pressure between unidirectional and bi-directional flows is called the purge pressure (Tan and Jaluria 1992). The purge pressure is a function of opening geometry and the buoyancy force (ΔP/(gΔρDH)) and may be expressed as (Cooper 1998):
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]
where:
w= Opening width [m]
D= Opening depth [m]
As mentioned above, when the air pressure difference between two zones is zero there is the maximum bi-directional flow due to the buoyancy force. When the pressure difference increases from 0 and is less than |ΔPFlood|, there is some bi-directional flow across the opening, but less than the maximum flow. If the pressure difference keeps increasing and exceeds |ΔPFlood|, there is no bi-directional flow. Cooper’s model assumes the buoyancy flow varies linearly with pressure difference.
The total air flow across the opening is based on superposition of the forced and buoyancy flows, and may be expressed for three different pressure difference scenarios as follows:
PL = PU
PL > PU
where:
ρL= Air density in the lower zone [kg/m3]
A= Opening area [m2]
Cd= Discharge coefficient [dimensionless]
Ρave = Average air density between the lower and upper zones [Pa]
ΔP= Pressure difference PL -
PU [Pa]PL < PU
Sloping plane
When a staircase is introduced as shown in Figure 139, 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):
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 140 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 140.
: Legend Description
Wind pressure calculations[LINK]
The wind pressure is determined by Bernoulli’s equation, assuming no height change or pressure losses:
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):
Cp is a function of location on the building envelope and wind direction. When Wind Pressure Coefficient Type = “INPUT”, the Cp values are explicitly defined in the input for AirflowNetwork:Multizone:Wind Pressure Coefficient Values. When Wind Pressure Coefficient Type = “AVERAGE-SURFACE CALCULATION” and the building shape is rectangular, the program uses the following equations to calculate wind pressure coefficient (Cp) values for different wind directions. For a low rise building, the normalized surface pressure coefficient may be written as (Swami and Chandra 1988):
where
Cp,n= Cp value at a given angle between wind direction and the outward normal of the surface under consideration [dimensionless]
α= Angle between wind direction and outward normal of wall under consideration [deg]
G= Natural log of the ratio of the width of the wall under consideration to the width of the adjacent wall [dimensionless]
n = Index of incident angle at 30-degree increments
For walls of a high rise building, a two-dimensional array of surface-averaged wind pressure coefficients is generated based on wind incident angle and side ratio. The wind pressure coefficients are provided in 2001 ASHRAE Fundamentals Handbook, p. 16.5, Fig. 7, “Surface Averaged Wall Pressure Coefficients for Tall Buildings”. The original work was performed by Atkins et al. (1979). The incident angle has an increment of 30 degrees. The side ratio values are 0.25, 1.0, and 4.0. For a given incident angle and building aspect ratio, the program uses linear interpolation to calculate the corresponding wind pressure coefficient Cp,n.
For the roof of a high rise building, a two-dimensional array of surface-averaged wind pressure coefficients is also generated based on wind incident angle and side ratio. The wind pressure coefficients are provided in 2001 ASHRAE Fundamentals Handbook, p. 16.6, Fig. 9, “Surface Averaged Roof Pressure Coefficients for Tall Buildings”. The original work was performed by Holmes (1986). The incident angle has an increment of 30 degrees. The side ratio values are 0.25, 1.0, and 4.0. At a given wind incident angle and building aspect ratio, the program uses linear interpolation to calculate the corresponding wind pressure coefficient Cp,n.
The wind surface pressure at the given incident angle can be calculated by combining the above two equations:
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:
where the correction vector, {C}, is computed by the matrix relationship:
{B} is a column vector with each component given by:
where n is the node number and i indicates all flow paths connecting node n to other nodes, and [J] is the square Jacobian matrix whose elements are given by:
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:
where
r= the ratio of Cn for the current iteration to its value for the previous iteration [dimensionless]
Cn= Correction value at the nth node [Pa]
Pn~~= Estimated pressure at the nth node [Pa]
Pn*~~= Corrected pressure at the nth node used in the next iteration [Pa]
This method is similar to a Steffensen iteration (Conte and de Boor 1972) which is used as a fixed-point iteration method for individual nonlinear equations.
The iteration correction method presented in the above equation gives a variable factor. When the solution is close to convergence, the solution method converges quadratically. By limiting cases where the value of r is less than some value, such as -0.5, the solution will not interfere with the rapid convergence. It has not been proven that the convergence acceleration equation will always lead to convergence, but it can be shown that it will not prevent convergence. Newton’s method converges when the estimated solution values are within some distance, called the radius of convergence, or the correct solution. Applying the convergence acceleration equation when -1 < r <0, will cause a smaller correction than Newton’s method, which therefore, can not force the iterations outside the radius of convergence. When r<-1, the solution diverges in an oscillatory fashion. When r>1, the solution also diverges, but in a nonoscillatory manner. For 0<r<1, the solution is approached from one direction. In all three cases, the convergence acceleration equation applies as long as r is truly constant over several iterations. However, for the last case, this involves a true extrapolation of correction factor which is very sensitive to the accuracy of r. This is most extreme for the case of r=1, which would cause an infinite correction.
Node Temperature Calculations[LINK]
A brief description of the air node temperature calculation is given below. A detailed description can be found in the work of Swami et al. (1992). The following equation is used to calculate temperature distribution across a duct element at the given airflow rate and inlet air temperature:
where
Cp= Specific heat of airflow [J/kg•K]
P= Perimeter of a duct element [m]
T= Temperature as a field variable [°C]
U= Overall heat transfer coefficient [W/m2•K]
hi= Inside heat transfer coefficient [W/m2•K]
ho= Outside heat transfer coefficient [W/m2•K]
tj= Thickness at j-th layer [m]
kj= Thermal conductivity at j-th layer [W/m•K]
The outlet air temperature at the end of the duct (x=L) is:
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:
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:
where
{M}= Airflow matrix
[T]= Temperature vector
[B]= Given boundary conditions
The zone air temperatures and primary air loop component (fan and coils) outlet conditions are used as prescribed conditions in the AirflowNetwork model. In addition, the temperature difference across zone loop components (terminal units) is held constant during the calculations. For example, thermal zone temperatures calculated during the previous system time step are used as prescribed temperatures when calculating all other node temperatures. The zone air temperature is assumed constant (prescribed) throughout the AirflowNetwork iterative solution. The fan and coil outlet air temperatures, and terminal unit temperature differences are assumed constant within an AirflowNetwork iteration. The sensible heat gains calculated during the AirflowNetwork solution are then used to predict a new zone air temperature.
Node Humidity Ratio Calculations[LINK]
A brief description of the air node humidity ratio calculation is given below. A detailed description can found in the work of Swami et al. (1992). The following equation is used to calculate humidity ratio distribution across a duct element at the given airflow rate and inlet air humidity ratio:
where
P= Perimeter of a duct element [m]
W= Humidity ratio [kg/kg]
Um= Overall moisture transfer coefficient [kg/m2•s]
hm,i~~= Inside moisture transfer coefficient [kg/m2•s]
*hm,o= Outside moisture transfer coefficient [kg/m*2•s]
tj= Thickness at j-th layer [m]
Dj= Moisture diffusivity at j-th layer [kg/m•s]
The outlet air humidity ratio at the end of the duct (x=L) is:
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
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:
where
{Mm}= Airflow matrix
[W]= Humidity ratio vector
[Bm]= Given boundary conditions
The zone air humidity ratios and primary air loop component (fan and coils) outlet conditions are used as prescribed conditions in the AirflowNetwork model. For example, thermal zone humidity ratios calculated during the previous system time step are used as prescribed humidity ratios when calculating all other node humidity ratios. The zone air humidity ratio is assumed constant (prescribed) throughout the AirflowNetwork iterative solution. The coil outlet air humidity ratio is assumed constant within an AirflowNetwork iteration. The latent heat gains calculated during the AirflowNetwork solution are then used to predict a new zone air humidity ratio.
Sensible and Latent Load Calculations[LINK]
The zone sensible and latent loads calculated in the AirflowNetwork model consist of multizone, duct conduction and leakage. In addition, the impact of infiltration and mixing is accounted for in this calculation. The multizone load only includes incoming airflows from outside (infiltration) and other adjacent zones (mixing) with and without forced-air fan operation. It is divided into two terms: variable and constant. The constant term is the sum of the mass flow rate multiplied by the specific heat for both infiltration and mixing. The variable term includes the impact of zone and outdoor air temperature. Each of these terms is used in the zone energy balance equation. The sensible load items from the multizone load calculations may be written as follows:
where
MCPairflow = Sum of air mass flow rate multiplied by specific heat for infiltration and mixing [W/K]
MCPTairflow = Sum of air mass flow rate multiplied by specific heat and temperature for infiltration and mixing [W]
Tamb= Outdoor air dry-bulb temperature [°C]
Tzone= Adjacent zone air temperature [°C]
The latent load items from multizone load calculations may be written as follows:
where
Mairflow = Sum of air mass flow rates for infiltration and mixing [kg/s]
MWairflow = Sum of air mass flow rate multiplied by humidity ratio for infiltration and mixing [kg/s]
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:
where
QADS,i= Total sensible load in the i-th zone due to ADS losses [W]
Qcond(ij)= Duct wall conduction loss at the j-th duct located in the i-th zone [W]
Qleak(ij)= Sensible supply leak loss at the j-th linkage located in the i-th zone [W]
QADS,m,i= Total latent load in the i-th zone due to ADS losses [kg/s]
Qcond,m(ij)~~= Duct wall vapor diffusion loss at the j-th duct located in the i-th zone [kg/s]
Qleak,m(ij) = Latent supply leak loss at the j-th linkage located in the i-th zone [kg/s]
Impacts of Supply Air Constant Volume Fan Control on Load: Cycling vs. Continuous[LINK]
The AirflowNetwork model currently allows two types of constant volume fans: Fan:ConstantVolume and Fan:OnOff. The Fan:ConstantVolume object has only one type of supply air fan operation mode: continuous fan, cycling compressor (ContinuousFanWithCyclingCompressor). However, the Fan:OnOff has two types of supply air fan operation modes: cycling fan, cycling compressor (CyclingFanAndCompressor) or continuous fan, cycling compressor (ContinuousFanWithCyclingCompressor). The CyclingFanAndCompressor operation mode is frequently referred to as “AUTO fan”, where the compressor(s) and supply air fan operate in unison to meet the zone heating or cooling load, and cycle off together when the heating or cooling load has been met. The ContinuousFanWithCyclingCompressor operation mode is often referred to as “fan ON”, where the compressor(s) cycle on and off to meet the zone heating or cooling load but the supply air fan operates continuously regardless of compressor operation. The supply air fan operation mode is specified in an HVAC system object based on a given fan operation mode schedule (e.g., AirLoopHVAC:UnitaryHeatCool object).
The determination of the zone sensible and latent loads caused by multizone airflows and forced air distribution system operation is dependent on the supply air fan operation mode (see Sensible and Latent Load Calculations section above). The zone loads calculated by the AirflowNetwork model are added to zone sensible and latent balances in the ZonePredictorCorrector module to calculate zone air temperatures and humidity ratios (see Integration of the AirflowNetwork Model section below). For the case of continuous fan/cycling compressor, the zone loads during forced air distribution system operation are calculated with the system design air mass flow rate without modification, since the system air node conditions (temperature and humidity) reflect the average values for the time step considering the coil/fan on and off periods during the time step.
For the case of cycling fan/cycling compressor, where the forced air distribution system can operate for a portion of the simulation time step, the airflows are determined based on the air distribution system part-load ratio (ratio of the average air mass flow rate for the time step divided by the design air mass flow rate). The airflows for the AirflowNetwork:Distribution:Linkage objects are reported during the air distribution system on cycle, since no airflow is assumed during the system off cycle. The airflows for the AirflowNetwork:Multizone:Surface objects are weighted by the air distribution system part-load ratio. The zone loads are the sum of energy losses during both the air distribution system on and off cycle at each system time step. The energy losses when the air distribution system is on are calculated using the system “on” air flow rate multiplied by the air distribution system run time fraction. The energy losses when the air distribution system is off are obtained from the multizone airflow calculations (without forced air through the air distribution system) and these losses are multiplied by (1.0 - system run time fraction), assuming no airflows through the air distribution system when the fan is off. The formulas used to calculate airflows and zone loads are given below:
Airflow
Airflow = Airflow during ADS on cycle * ADS Part-load ratio + Airflow during ADS off cycle * (1.0 – ADS Part-load ratio)
where ADS = Air Distribution System
Zone load
System run time fraction = Maximum run time fraction of coils and fans in the air distribution system
Zone energy losses = Zone energy loss during ADS on cycle * System run time fraction + Zone energy loss during ADS off cycle * (1.0 - System run time fraction)
The calculation of loads due to multizone airflow, without forced air distribution system operation, is performed when the HVAC system is first simulated during a simulation time step (FirstHVACIteration = True). The calculation of loads due to air distribution system operation is performed on subsequent iterations within the same simulation time step when the mass flow rate at the supply air fan inlet is greater than 0.0 (and FirstHVACIteration = False).
Airflow Calculation Procedure using A Supply Variable Air Volume Fan[LINK]
The AirflowNetwork model currently also allows a variable air volume fan type as Fan:VariableVolume. The allowed terminal type is AirTerminal:SingleDuct:VAV:Reheat only. Other types of terminals will be added later.
In general, the supply fan air flow rate in a VAV central system is determined by a sum of terminal flow rates when the AirflowNetwork model is not applied. When the AirflowNetwork model is applied and the supply air fan flow rate is given, each terminal flow is determined by pressure resistance of each supply air pathway. It is possible that the delivered air flow rate from the pressure resistance at each terminal may be totally different from the desired flow rate determined by terminal units. Therefore, it is not easy to meet both requirements. The following compromised approach, including possible supply and return leaks in an air distribution system, is implemented.
Set up terminal airflows in the AirflowNetwork module based on the SimVAV subroutine in the HVACSingleDuctSystem module.
Require AirflowNetwork:Distribution:Component:LeakageRatio objects to define supply leaks, so that the values of the Effective Leakage Ratio field are used to decide the supply fan flow rates. The base of the ratio will be actual supply fan flow rate, instead of the maximum fan flow rate used in the constant volume fan systems.
Assign the supply fan airflow rate based on the sum of all terminal flow rates and all supply leak ratios until it reaches the maximum fan flow rate
where
n= Number of terminals
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. ****
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
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:
Where MCPTairflow is the sum of mass flow rate multiplied by specific heat and temperature for infiltration and mixing, QADS,z is the added total sensible load in the zone due to Air Distribution System losses, and MCPairflow is the sum of mass flow rate multiplied by specific heat for infiltration and mixing as calculated in the AirflowNetwork model described above.
The revised coefficient (B) used in the zone humidity ratio calculation is shown below:
Where MWairflow is the sum of mass flow rate multiplied by humidity ratio for infiltration and mixing and QADS,m,z~~is 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 Predictor-Corrector” found elsewhere in this document.
Model Output[LINK]
The available outputs from the AirflowNetwork model are described in the EnergyPlus Input Output Reference manual.
References[LINK]
Atkins, R. E., J. A. Peterka, and J. E. Cermak. 1979. “Averaged pressure coefficients for rectangular buildings,” Wind Engineering, Proceedings of the Fifth International Conference 7:369-80, Fort Collins, CO. Pergamon Press, NY.
Bolmqvist, C. and M. Sandberg, 2004, “Air Movements through Horizontal Openings in Buildings – A Model Study,” International Journal of Ventilation, Vol. 3, No. 1, pp. 1-9
COMIS Fundamentals. 1990. Edited by Helmut E. Feustel and Alison Rayner-Hooson, LBL-28560, Lawrence Berkeley Laboratory, Berkeley, CA
Conte, S. D. and C de Boor. 1972. Elementary Numerical Analysis: an Algorithmic Approach, McGraw-Hill.
Cooper, L., 1989, “Calculation of the Flow Through a Horizontal Ceiling/Floor Vent,” NISTIR 89-4052, National Institute of Standards and Technology, Gaithersburg, MD
Dols, W. S. & G. N. Walton. 2002. “CONTAMW 2.0 User Manual,” NISTIR 6921, National Institute of Standards and Technology, Gaithersburg, Maryland
Holmes, J. D. 1986. Wind Loads on low-rise buildings: The structural and environmental effects of wind on buildings and structures, Chapter 12, Faculty of Engineering, Monash University, Melbourne, Australia
Swami, M. V. and S. Chandra. 1988. Correlations for pressure distribution on buildings and calculation of natural-ventilation airflow, ASHRAE Transactions 94(1988) (Pt 1), pp. 243-266.
Swami, M. V., L. Gu, & V. Vasanth. 1992. “Integration of Radon and Energy Models for Building,” FSEC-CR-553-92, Florida Solar Energy Center, Cocoa, Florida
Tan, Q. and Y. Jaluria, 1992, “Flow through Horizontal Vents as Related to Compartment Fire Environments,” NIST-GCR-92-607, National Institute of Standards and Technology, Gaithersburg, Maryland
Walton, G. N. 1989. “AIRNET – A Computer Program for Building Airflow Network Modeling,” NISTIR 89-4072, National Institute of Standards and Technology, Gaithersburg, Maryland
Walton, G. N. & W. S. Dols. 2003. “CONTAM 2.1 Supplemental User Guide and Program Documentation,” NISTIR 7049, National Institute of Standards and Technology, Gaithersburg, Maryland
Documentation content copyright © 1996-2014 The Board of Trustees of the University of Illinois and the Regents of the University of California through the Ernest Orlando Lawrence Berkeley National Laboratory. All rights reserved. EnergyPlus is a trademark of the US Department of Energy.
This documentation is made available under the EnergyPlus Open Source License v1.0.