NISTIR 89-4072

~~14'f 01' "0

~~ ~ ~ '+~

: ~

* ~ *

Co ~

~ ..$I

~~ ~q.'

or ...~ 1'-17'£5 Of

George N. Walton

u.s. DEPARTMENT OF COMMERCE

National Institute of Standards and Technology

(Formerly National Bureau of Standards)

National Engineering Laboratory

Center for Building Research

Building Environment Division

Gaithersburg, MD 20899

April 1989

Sponsored by:

U.S. DEPARTMENT OF COMMERCE

National Institute of Standards and Technology

Gaithersburg, MD 20899

u.s. Department of Energy

Washington, DC 20585
George N. Walton

u.s. DEPARTMENT OF COMMERCE

National Institute of Standards and Technology

(Formerly National Bureau of Standards)

National Engineering Laboratory

Center for Building Research

Building Environment Division

Gaithersburg, MD 20899

April 1989

~~~~1 0' ('0++ National Bureau of Standards became the

~~q. ~ ~ National Institute of Standards and Technology

: ~ on August 23, 1988, when the Omnibus Trade and

* ~

.* Competitiveness Act was signed. NIST retains

\ § all NBS functions. Its new programs will encourage

~4'4 ~$ improved use of technology by U.S. industry.

.r1'04rrs of to'

Sponsored by:

U.S. DEPARTMENT OF COMMERCE

National Institute of Standards and Technology

Gaithersburg, M D 20899

U.S. Department of Energy

Washington, DC 20585

u.s. DEPARTMENT OF COMMERCE

Robert Mosbacher, Secretary

NATIONAL INSTITUTE OF STANDARDS

AND TECHNOLOGY

Raymond G. Kammer, Acting Director
ABSTRACT

In spite of its importance, the analysis of airflows has significantly lagged

the modeling of other building features because of limited data, computational

difficulties, and incompatible methods for analyzing different flows. Methods

have been developed to analyze airflows in HVAC ducts and to estimate infiltration but the interaction between building HVAC systems and infiltration

airflows has seldom been studied. This report describes a computer program

for modeling networks of airflow elements, such as openings, ducts, and fans.

It emphasizes the numerical aspects of an airflow network method which would

provide a unified approach to building airflow calculations. It also

discusses the limitations of the method and poorly understood factors that

could profit from further research.

iii
ACKNOWLEDGEMENTS

This work was sponsored by the US Department of Energy, Office of Building

and Community Systems, contract DEAIO136CE21013 modification no. 003. The

author also wishes to thank Dr. James Axley for the long discussions on many

of the ideas used in the AIRNET program.

DISCLAIMER

The program described in this report is furnished by the government and is

accepted and used by any recipient with the express understanding that the

United States Government makes no warranty, expressed or implied, concerning

the accuracy, completeness, reliability, usability, or suitability for any

particular purpose of the information and data contained in this program or

furnished in connection therewith, and the United States shall be under no

liability whatsoever to any person by reason of any use made thereof. This

program belongs to the government. Therefore, the recipient further agrees

not to assert any proprietary rights therein or to represent this program to

anyone as other than a government program.

iv
CONTENTS

ABSTRACT ..iii

iv

ACKNOWLEDGEMENTS.

DISClAIMER iv

1 INTRODUCTION

2. AN AIRFLOW NETWORK METHOD. .

1. Modular Approach. ...

2. Newton's Method. ...

3. Solution of the Equations

4. Convergence Criteria. .

5. Convergence Acceleration

6. Linear Initialization.

2

2

2

3

4

4

5

..

~

3. 5

6

7

9

11

13

~

1. ELEMENT Powerlaw MobELS Flow Elements

2. Ducts .."'...

5. 3. 4 .Fans Quadratic Doorways Flow Elements

4 13

14

14

2. 1SAMPLE A Simple CALCULATIONS. Comparison Test Cases. of Methods. .

5 15

15

16

17

DIRECTIONS FOR FUTURE RESEARCH. .

1. Alternate Solution Methods. .

3. 2. Other General Element Limitations. Models. 

6 SUMMARY AND CONCLUSIONS. ... 17

7 REFERENCES 18

FIGURES 20
Appendix A: AIRNET User's Manual. .

3. 2. 1. The Intractive Introduction AIRNET data input. files. ..

4. Wind pressure calculations. ...

5. Equation ordering in the Jacobian

6. 7Sample Compiling output AIRNET file. 

31

31

31

32

36

36

37

40

.. ..

43

43

46

51

53

55

56

58

59

61

63

64

Appendix 7. 6. 910. 2. 1B: Quadratic Executiontimetests Powerlawelementtests Doorwayelementtest Constant Introduction AIRNET flow Validation flow element element Tests. tests. test. 

65

65

67

69

72

74

75

75

76

77

Appendix \8. 7. 5. 6. 1. C: Relationshipsforcracks Newton'smethod Fitting Constant Calculation powerlaw power flowelement Details. fan coefficients. element. 

LIST OF FIGURES

20

21

22

23

24

25

26

27

28

29

30

Figure 1. Plan of a Portion of a Building Served by a VAV System

Figure 3. 54. 2. Oscillating Typical Airflow Dooray Network Fan Coordinate Pressure Performance for System VAV Corrections. Curve. System and Three of Figure Flow Patterns. 1 ...

Figure 8. 6. 7910. Three One Fans Powerlaw Floor in Airflow Series Parallel Airfow of the Networks (from Elements 4 (from Floor Osborne) for Osborne) / in Doorways. 36 Series. Room Timing Test Case. ..

Figure 11. Proposed Structure of Airflow Network Program.

vi
1 INTRODUCTION

Air movement models have been developed for estimating airflows in buildings.

These airflows include infiltration, natural ventilation, interroom airflows

through various openings including doorways, and flows through the HVAC

system. The numerical estimation of average characteristics of such airflows

is useful for moisture and contaminant dispersal analysis, including the

design of smoke control systems, and heat transfer analysis including load

and energy calculations. In spite of its importance, the analysis of airflows

has significantly lagged the modeling of other building features because of

limited data, computational difficulties, and incompatible methods for

analyzing different flows. This is particularly true of the combined building

and HVAC system simulation. Methods have been developed to analyze airflows

in ducts (ASHRAE, 1985, ch 33) and to estimate infiltration (Liddament &

Thompson, 1982) and ventilation (ASHRAE, 1985, ch 22), but the intimate

relationship between these processes has seldom been studied. When it has,

the results have sometimes been surprising (Persily, 1985).

Relatively few methods that could be applied to both processes have been

developed within the building simulation community and described in detail.

Several computer models developed for smoke control analysis are reviewed by

Said (1988). Models for building energy analysis have been developed by

Clarke (1985) and Walton (1984). All of these methods are based upon the idea

that there is a simple nonlinear relationship between the flow through an

opening and the relative air pressure difference across it, and that a

building can be considered to be composed of a large number of rooms which are

connected by openings to each other and to the outside. This is a network of

rooms (nodes) and openings (connections) which is conceptually similar to the

air handling system network where the connections are the ductwork and the

nodes are the ductwork junctions. Conservation of mass for the flows into and

out of each node leads to a set of simultaneous nonlinear equations which are

solved iteratively for the airflows. This can be called an "airflow network"method. 

Its relationship to pipe network methods will be discussed. Such an

analysis is also sometimes referred to as a multi-chamber or multi-cell method

(ASHRAE, 1985, p 22.13). This report draws extensively on Axley's airflow

element (1987) and contaminant element (1988) methods which are, in turn,

based on numerical methods associated with finite element modeling techniques.

Modeling of airflows requires: (1) determination of the location and

mathematical characterization of the airflow paths, (2) determination of the

boundary conditions (primarily wind pressure), (3) calculation of the

resulting airflows, and (4) a user-friendly framework in which to do the

analysis. Progress has been made in such vital areas as wind pressure

estimates (Swami & Chandra, 1988) and interroom airflows (Barakat, 1987).

Unfortunately, it is often thought that a network model is so complex that it

requires a mainframe computer for its solution (ASHRAE, 1985, p 22.13) and is

therefore, impractical. This apparent impracticality discourages gathering

the data which is necessary for the use of network models.

This report will shown that a network model is practical. It will emphasize

the numerical aspects of the airflow network method which allow it to provide

a practical, unified approach to building airflow calculations. Details of

the program AIRNET, a microcomputer implementation of this airflow network

method, will be discussed. It will also discuss the limitations of the method

1
and poorly understood factors that could profit from further research.

2. AN AIRFLOW NETWORK METHOD

An airflow network consists basically of a set of nodes connected by airflow

elements. The nodes may represent rooms, connection points in ductwork, or

the ambient environment. The airflow elements correspond to discrete airflow

passages such as doorways, construction cracks, ducts, and fans. Figure 1 is

a sketch of a portion of a building consisting of two rooms, a hallway, and

air distribution equipment representing a VAV system. Figure 2 shows an

airflow network superimposed on the physical structure of figure 1. The large

dots are nodes and the connecting lines are the various airflow elements.

2.1 Modular Approach

The network approach makes the development of element models, excitation

models, and solution method somewhat independent. the computer program

modules will obviously mirror the theoretical modules with input and output

modules added to create a useful simulation tool. The various modules provide

a toolkit for the analyst to consider a practically infinite variety of systemmodels.

For this study an airflow network simulation computer program, AIRNET, was

developed from an earlier airflow analysis program (Walton, 1984). The new

program consists of:

(1) a process for establishing an initial set of values to start the

iterative solution process,

(2) a solution method for nonlinear equations consisti~g of a traditional

Newton's method combined with Steffensen iteration to accelerate convergence,

(3) airflow element subroutines which compute the flow rate and derivative of

the flow with respect to pressure difference needed to form the Jacobian

matrix used in Newton's method,

(4) a separate process for transferring the above data into the Jacobian

matrix (called the element assembly process), and

(5) solution of the simultaneous linear equations involving the Jacobianmatrix.

discussion will begin with the solution method

Newton's Method

Each airflow element, i, relates the mass flow rate, wi' through the element

due to the pressure drop, APi' across it. Conservation of mass at each node

is equivalent to the mathematical statement that the sum of the mass flows

equal zero (or the mass generated, as in the case of a fire) at each node.

The flows are related nonlinearly to the pressures at the nodes thus

requiring the iterative solution of a set of nonlinear equations. In Newton's

method (Conte & de Boor, 1972, p 86), a new estimate of the vector of all node

pressures, {P}*, is computed from the current estimate of pressures, {P}, by

{P}* = {P} -{C} (1

2
where the correction vector, {C}, is computed by the matrix relationship

[J] {C} = {B} (2)

{B} is a column vector with each element given by

Bn = ~ Wi

i

(3)

where n is the node number and i indicates all flow paths connecting node n to

other nodes, and [J] is the square (i.e. Nby N for a network of N nodes)

Jacobian matrix whose elements are given by

= ~ .<2..!i l. 8P

.m1.

In.m

In equations (3) and (4) wi and 8wi/8Pm are evaluated using the current

estimate of pressure {P}. The AIRNET program contains subroutines for each

airflow element which return the mass flow rates and the partial derivative

values for a given pressure difference input.

2.3 Solution of the Equations

Equation (2) represents a set of linear equations which must be set up and

solved for each iteration until a convergent solution of the set of nonlinear

equations is achieved. In its full form [J] requires computer memory for N2

values, and a standard Gauss elimination solution has execution time

proportional to N3. Sparse matrix methods can be used to reduce both the

storage and execution time requirements.

A skyline solution process following the method of Dhatt (1984, pp. 282-192)

was chosen. This method can be used to solve equations with symmetric or

nonsymmetric matrices. It stores no zero values above the highest nonzero

element in the columns above the diagonal and no zero values to the left of

the first nonzero value in each row below the diagonal. Analysis of the

element models will show that

Note that the degree of sparsity of Jacobian matrix is dependent on the

ordering of the nodes. Ordering can be improved by various algorithms or

rules-of-thumb. Also note that it is easy to define an airflQw network which

has no unique solution. One requirement for solution is that at least one of

the node pressures be known. This is usually the ambient node. All nodes

must be linked, through some path, to a known pressure. There may be several

known pressure nodes. The airflow network method allows two types of nodes:

those with known or unknown pressures. In AIRNET the constant pressure nodes

are included in the system of equations and equation (2) processed so as to

not change those node pressures. This gives an added flexibility in defining

the airflow network with special processing maintaining the symmetric set of

3
equations. The form of the equations for known pressure nodes combined with

the condition in equation (5) for unknown pressures nodes is a sufficient

condition for the Jacobian to be nonsingular (Axley, 1987). AIRNET ispresently. 

set up so that the ambient node pressure is zero causing the

computed node pressures to be values relative to the true ambient pressure.

This helps maintain numerical significance in calculating ~P.

Convergence Criteria

Conservation of mass at each node provides the convergence criterion. That

is, if ~ wi = 0 for all nodes for the current system pressure estimate, the

solution has converged. Many iterations can be saved and sufficient accuracy

attained by testing for relative convergence at each node

I~Wil /~IWil < f

(6)

with a test to prevent division by zero. The size of ! can be established by

considering the use of the calculated airflows, such as in an energy balance.

Convergence Acceleration

Numerical tests of Newton's method solution indicated occasional instances of

very slow convergence, always with oscillating corrections on successive

iterations. This is depicted graphically for the successive values of

pressure at a single node in figure 3. In the case shown each successive

pressure correction is a constant ratio of the previous correction. The

observed corrections come close to this pattern. By assuming a constant

ratio, it is simple to extrapolate the corrections to an assumed solution:

Pn* = Pn -Cn/(l~r) (7)

where r is the ratio of Cn for the current iteration to its value for the

previous iteration. This extrapolated value of node pressure is used in the

next Newton iteration. At every other iteration, there are two pressure

correction values which may be used for an extrapolation. This method is

similar to a Steffensen iteration (Conte & de Boor, 1972, p 54) which is used

with a fixed point iteration method for individual nonlinear equations.

The oscillating corrections have been observed by other investigators (Wood,

1981; Demuren, 1986). Demuren uses a constant relaxation factor of 0.5 to

prevent the oscillations. The iteration correction method presented in

equation (7) gives a variable factor. When the solution is close to

convergence, Newton's method iterations converge quadratically. By limiting

the application of equation (7) to cases where r is less than some value such

as -0.5, it will not interfere with the rapid convergence. Tests by the

author confirm that this is faster than the constant relaxation factor. It

has not been proven that equation (7) 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, of the correct solution. Applying equation (7) when

-1 < r < 0, will cause a smaller correction than Newton's method, which,

therefore, cannot force the iterations outside the radius of convergence.

4
The meanings of other values of r are also interesting. 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, equation (7) 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.

2.6 Linear Initialization

Newton's method requires an initial set of values for the node pressures.

These may be obtained by including in each airflow element model a linear

approximation relating the flow to the pressure drop:

+ bioAP (8) w. ;:: C 1

Conservation of mass at each node leads to a set of linear equations of the

form

[A] {P} = (B)

The coefficient matrix [A] in equation (9) has the same sparsity pattern as

[J] in equation (2) allowing use of the same sparse matrix solution process

for both equations. This initialization handles stack effects very well and

tends to establish the proper directions for the flows. The linear approximation is conveniently provided by the laminar regime of the element models

described below, but it also may be provided by a secant approximation to the

actual nonlinear behavior.

The initialization has been made optional in AIRNET. When solving a set of

similar problems, such as when the node temperatures or wind pressures are

changed by small amounts, it may be preferable to use the previous solution

for the node pressures as the initial values for the new problem.

3 ELEMENT MODELS

Flow within each airflow element is assumed to be governed by Bernoulli'sequation:

~p = (Pl + pV12/2) -(P2 + pV22/2) + pg(Zl -Zz)

where

~p = total pressure drop between points 1 and 2

Pl, Pz = entry and exit static pressuresVl' 

Vz = entry and exit velocities

p = fluid density

g = acceleration of gravity (9.81 m/sz)

Zl' Zz = entry and exit elevations.

The following parameters apply to the nodes: pressure, temperature (to

compute density and viscosity), and height. The node height values are used

to determine stack effect pressures. When the node represents a room, the

5
airflow elements may connect with the room at other than its reference height.

Appendix C.3 shows how to use the hydrostatic equation to relate the pressure

difference across a flow element to the heights of the element ends and the

node heights, assuming the air in the room is at constant temperature.

Pressure terms can be rearranged and a possible wind pressure for building

envelope openings added to give

AP = Pn -Pm + PS + PW

where

Pn' Pm = total pressures at nodes n and ill

PS = pressure difference due to density and height differences, and

PW = pressure difference due to wind.

Equation (11) establishes a sign convention for direction of flow: positive

is from node n to node m. Since the airflow elements will be described by a

relationship of the form w = f(AP), the partial derivative needed for [J] in

equation (4) are related by 8w/8Pm = -8w/8Pn which establishes the relation

in equation (5).

3.1 Powerlaw Flow Elements

Most infiltration models are based on the following empirical (powerlaw)

relationship between the flow and the pressure difference across a crack or

opening in the building shell:

(AP>O)

or

Wi : C JPn (6P)X

Wi : -C JPm (-6P)X (AP<O)

where

w. = mass flow rate of air through element i from node n to node m,

].C = flow coefficient,

p = air density of node n or m,

AP = total pressure loss across the element (Pn -Pm), and

x = flow exponent.

Theoretically, the value of the flow exponent should lie between 0.5 and 1.0.

Large openings are characterized by values very close to 0.5, while values

near 0.65 have been found for small crack-like openings. The form of

equations (12) is suggested by the orifice equation:

Q = Cd A JUPlP

where

Q = volumetric flow (Q = w/p) ,

Cd = discharge coefficient, and

A = orifice opening area.

Equation (12) should be considered a correlation rather than a physical law.

It can be used with the element leakage area formulation which has been used

to characterize openings for infiltrations calculations (ASHRAE, 1985, p

22.16). The author has used it to describe flows through ducts to an accuracy

of about 2% over a range of flow rates that vary by a factor of four. Such a

variation would be found in a VAV system.

6
The primary advantage of equations (12) for describing airflow components is

the simple calculation of the partial derivatives for the Newton's method

solution of the simultaneous equations:

(14a)

£!i -

8Pm --x Wi / A-P (14b)

However, there is also a problem with equations (14): the derivatives become

unbounded as the pressure drop (and the flow) go to zero. A simple way to

avoid this problem is suggested by what physically happens at low flow rates:

the physical character of the flow (and the form of the equation) changes. It

goes from turbulent to laminar. Equation (12) can be replaced by

w:: K p lJ.P /.u (15) where

K = flow coefficient,

p. = viscosity.

The partial derivatives are simple constants

~i

8Pn

= K p/..u (16a)

£.!i = -K p/,u (16b)

8Pm

The origin of this laminar relationship is shown by the duct equations in the

next section. This technique has been independently discovered and used by

several researchers (Axley, 1987; Isaacs, 1980). Although there is physical

reason for using equation (15) at low pressure drops, its real purpose is to

assure convergence of the equations when AP approaches zero for one of the

many flow paths in a complex network, instead of accurately representing

airflows which are too small to be of interest.

The AIRNET function for powerlaw elements calculates flows using both the

laminar and the turbulent mod~ls, equations (12) and (15), and selects the

method giving the smaller magnitude flow. There is a discontinuity in the

derivative of the w(AP) curve where the two equations intersect. This discontinuity is a violation of one of the sufficient conditions for convergence of

Newton's method (Conte & de Boor, 1972, p 86). However, numerical tests

conducted by the author for flows at that point using a small airflow network

have shown no convergence problem.

Ducts

The theory of flows in ducts (and pipes) is well established and summarized in

the ASHRAE Handbook of Fundamentals (ASHRAE, 1985, ch 33). More extensive

treatment is given by Blevins (1984) in a long chapter on pipe and duct flow.

Analysis is based on Bernoulli's equation and its assumptions.

The friction losses in a section o.f duct or pipe are given by
A-Pf = foL/DopV2/2

where

f = friction factor,

L = duct length, and

D = hydraulic diameter.

The dynamic losses due to fittings and so forth are given by

l1Pd = Co pV2/2

where

Co = dynamic loss coefficient

Total pressure losses are given by

~p = ~Pf + ~ ~Pd

Since w = pVA, where A is the cross section area,

w = 2pA2/(fL/D + ~ Co)]\ .Ap\ (20)

The friction factor can be computed using the nonlinear Colebrook equation

(ASHRAE, 1985, p 2.9)

where

f = roughness dimension, and

Re = Reynolds number = pVDjp = wDjpA.

This nonlinear equation may be readily solved using the following iterative

expression derived from equation (21) by Newton's method:

= g -[g -a + ~ In(l+g~) g' / [1 + -yfJ/(l+gfJ)

where

g::: 1 / if,

Q ::: 1.14 -1 In(E/D),

p::: 9.3/(ReoE/D), and

1::: 2o1og(e) ::: 0.868589

The convergent solution is achieved in 2 or 3 iterations of equation (22)

using g = Q as a starting value. If the value of g has been saved from the

previous time it was computed for a particular duct element, and the flow rate

has not changed greatly, only one iteration of equation (22) will be needed to

compute the friction factor.

The exact derivatives of equation (20) are difficult to compute. However,

reasonable convergence is achieved by assuming the the coefficients in

equation (20) are constant giving

_1 8w = ~ w. / l1.P 8P J. (

n

and

aWi 8P = -~ Wi / AP

m

(23b)

8
The above description of duct flow applies only in the fully turbulent flow

regime above a Reynolds number of about 4000. When the Reynolds number is

below about 2000, the flow is laminar with the laminar friction loss

described by

~Pf = k/ReoL/DopV2/2

= p./pokL/(2AD2) ow

where

k = laminar friction factor.

Laminar dynamic losses are given by

APd = Ko pV2/2

= Ko/(2pA2).w2

where

Ko = laminar dynamic loss coefficient

Expressions (24) and (25) lead to a quadratic equation for mass flow in terms

of total pressure drop:

a w2 + b w + c = 0

where

a = Ko/(2pA2)

b = jjkL/(2pAD2)

C = lAP!

giving

w = sign(AP) [(b2+4ac)% -b] / 2a

The partial derivatives are given by

8w 1

aFn: (b2+4ac)%

8w -1

aFro: (b2+4ac)%

The derivatives at AP: 0 are finite (i.e :tl/b)

3.3 Doorways

Flows through large openings (e.g. doorways) tend to be more complex with the

possibility of flows in opposite directions in different parts of the

opening. The temperature and resulting density differences between two rooms

may mean that the stack effect causes a positive pressure difference at the

top of the doorway and a negative pressure difference at the bottom (or viceversa

A summary of research on heat transfer through doorways is presented

by Barakat (1987). Most research has attempted to develop dimensionless

correlations (using Nusselt, Prandlt, and Grashoff numbers) of the form

NuD/Pr = CoGrDb

where

b is approximately 0.5 and

9
C lies between 0.22 and 0.33 depending on the temperature difference used

for the correlation. It has been shown that such a heat transfer is equivalent to an airflow which can be modeled by powerlaw elements (Walton, 1982)

by dividing the total opening into several smaller openings having the same

total area but configured to properly account for the magnitude and direction

of airflows at different heights in the opening.

An alternative approach is to create a single airflow element which accounts

for the flow over the entire opening. A simple theory which estimates the

stack induced air flow through a large opening in a vertical partition is

given by Brown and Solvason (1962). The derivation of the doorway element is

based on the model shown in Figure 4.

By assuming that the air density in each room is constant, the hydrostatic

equation is used to relate pressures at various heights in each room: I

Pan = Pn + Png(~ -ho) and POrn = Pm + Pmg(~ -ho)

Pn (Y) = POn -Pngy and Pm (Y) = POm -Pmgy

Following Brown & Solvason (1962) it is assumed that the velocity of t~e

airflow as a function of height is given by the orifice equation: i

V(y) = Cd[2(Pn(y)-Pm(y))/p]1/2

where

Cd = discharge coefficient, and

p = density of the air going through the area

At the neutral height, Y, the velocity of the air is zero.

this must occur when Pn(y) = Pm(y). From equations (30)

From equation (31

y = Pnn-Pnm

g(pn-pm)

-Lnm p -pLnn

-g(pm-pn)

If 0 < Y < H, there is a two-way airflow through the doorway. If Pn = Pm, the

neutral height cannot be computed, but, since there is no possibility ~f twoway flow, the doorway can be considered a simple orifice opening. I

Define ~p = Pn-Pm and a transformed height coordinate z = Y -y. Then the

pressure difference across the doorway is given by

A.P(z) = -gzA.p

The mass flow through the doorway above the neutral height is given by

z=H-Y

J (pV)i W dz

z=O

w a= i = n or m

and the mass flow below the neutral height by

10
z=o

I (pV)i W dz

z=-Y

Wb = i = n or m

Whether the subscript i should be nor m depends on the direction of flow.

Integration of equations (34) and (35) gives several different solutiops for

the airflow depending on the value of Y and the sign of ~p. IiDefining:

fa = /H-y/3/2 fb = ly13/2 (36a) and

G' 

= W Cd [2/(gIAPI)]1/2 fa' == IH-yll/2 fb' = IYll/2 (36b)

gives the following equations for flows and derivatives

Case 1 Y:s:O

Ap > 0 w = -G JP: Ifa-fb 8w/8Pn = G' JP:- Ifa'-fb'

Ap < 0: w = G JP:" I fa -fb 8w/8Pn = G' J~ Ifa'-fb'

Case 2 Y~H

/J.p > 0

~p < 0: (37d)

Case 3 O<Y<H

= -G JP:- fa

= G JP-: fa

= G JP-: fa

= -G JP:- fa

~p > 0 8w/8Pn = G' JP:- fa'

JP-: fb'

JP:- fa'

JP:- fb'

Wa (37e)

8w/8Pn = G' w, (37f)

Ap < 0: 8w/8Pn = G' Wa (37g)

8w/8Pn = G' w, (37h)

with aw/apm = -8wj8Pn (37i)

This model of a doorway tends to be faster that the multiple opening

approach. However, it also complicates the assembly process for the Jacobian

matrix because one or two flows may exist (i.e., case 3 above). More importantly, development of the doorway element model requires knowledge of the

vertical temperature profile used in the node model (here assumed to be

constant) in order to compute the pressure difference as a function of height

across the opening. This requirement compromises the independence of the

modularity of airflow network program.

3.4 Fans

The theory of flows induced by fans is summarized in the ASHRAE Equipment

Handbook (ASHRAE, 1983, ch 3). More extensive treatment is given by Osborne

(1977). Fan performance is normally characterized by a performance curve such

11
as shown in figure 5. This curve relates the total pressure ~ to the flow

rate for a given fan speed and air density. The fan performance curve is well

represented by one or more cubic polynomials:

P = ao + al W + aZ wZ + a3 W3

Multiple polynomials might conveniently be obtained by a cubic spline fit to

the performance data. There are two important factors to note on the shape of

the fan performance curve. First, it is described by a relationship of the

form P(w) instead of w(P) which would be more appropriate for the calculation

of flow and partial derivatives. The basic shape of. the performance curve

cannot be well represented by a simple polynomial with P as the independentvariable. 

Equation (38) requires an iterative solution to determine the flowrate. 

A modified false-position method (Conte & de Boor, 1972, p 31) works

quickly and reliably. Fortunately, the derivative dw/dP is simply the

inverse of dP/dw, which is a simple expression for a polynomial.

Second, it is common for the performance curve to contain points of contraflecture, with up to three different flow rates possible at certain values of

fan pressure. This causes difficulty in solving for the flow rate and, more

importantly, has points where the derivative goes to infinity. However, it is

usually not recommended that the fan operate in the region of the contraflecture points. Therefore, the fan can be modeled with a performance curve

that does not include the contraflecture so long as the user checks that the

air distribution system does not permit operation in that region.

It is easy to identify the points of contraflecture from the coefficients of

the polynomial by the evaluation of simple derivatives:

P = aO + al W + aZ W Z + a3 W3

(39a)

p' = al + 2a2w + 3a3w2 (39b)

P" = 2a2 + 6a3w

p' = 0 at the points of contraflecture Solving equation (39b) for w gives:

If (4a22 -12ala3) > 0, there are two real points of contraflecture, with

equation (39c) defining the highest root (a3 must be negative to give the

typical fan curve). If P"(w) > 0, the point is a maximum; if P" = 0, it is a

point of inflection.

The performance of a given fan at various speeds and air densities can be

related to a single fan performance curve through the "fan laws".

Wl / Wz = (N1Pl) / (Nzpz)

and

Pl / P2 = (N12Pl

/ 

(N22p2)

where

w = volume flow rate,

P = total pressure rise,

N = rotational speed,

12
p = density

These laws are valid if all flow conditions at the two speeds are similar.

particular, they will not apply at very low flows where fully turbulent

conditions have not been developed.

In

Numerical tests with AIRNET for flows at the laminar-turbulent transition

indicate some convergence difficulty: about twice as many iterations as usual

are needed for convergence. In one case the iterations showed potential

divergence with r < -1, but the convergence acceleration algorithm saved the

cases tested and produced a solution.

3.5 Quadratic Flow Elements

Baker, Sharples, and Ward (1987) indicate that infiltration openings can be

more accurately modeled by a quadratic relationship of the form

AP = A Q + B Q2 (Q,6P ~ 0) (43a) and

AP = A Q -B Q2 (Q,AP:S 0) (43b)

This form can be used as an airflow element by solving the quadratic equation

for w (= pQ). Letting a = Alp and b = B/p2 allows equations (43a,b) to be

rewritten as

~p = a w + b w2 (w,6P ~ 0)

and

AP = a w -b w2 (W,l:.P:S 0)

These quadratic equations solve as

w = (Ja2+4bAcP -a) I 2b (AP ~ 0)

and

w = (a -Ja2-4bAP) I 2b (AP:S 0)

The partial derivatives are given by

8w/8Pn = + 1 / (a+2blwl) (46a)

8w/8Pm = -1 / (a+2b Iwl) (46b)

Equations (45) and (46) require that a and b both be nonzero to prevent a

division by zero. There is no problem as ~p and w go to zero --the

derivatives are finite.

4 SAMPLE CALCULATIONS

Several simple airflow networks have been analyzed to demonstrate the

procedure described in the previous sections.

13
Simple Test Cases

Figure 6 shows two cases involving powerlaw elements in series. The first

case, consisting of three nodes and two flow elements, can be considered as

modeling a room with small 0.01 m2 openings on opposite sides with wind

pressure driving flow through the room. The second case divides the single

room into two room with a partition containing a large 2.00 m2 opening. This

case with a very low resistance (large opening) mixed with large resistances

(small openings) is difficult to solve with some methods (Walton, 1982;

Clarke, 1985, p 206). In both cases AIRNET required only two iterations and

computed the expected nearly identical flows.

Figure 7 shows three cases involving doorway elements. In the first case a

0.8 m by 2.0 m doorway connects two rooms with a 4°C temperature difference.

The computed two-way airflow is 0.259 kg/so In the second case ten 0.16 m2

airflow openings at different heights are used to represent a doorway. The

computed two-way airflow is 0.261 kg/so The third case represents six rooms

in series connected by doorway elements. The computed flows are identical to

the first case. All three cases were solved in two iterations.

Figure 8 shows a test involving two fans in parallel. This is a problem in

the textbook by Osborne (1977, p 75). The computed pressures agree with the

text to within 2 Pa (about 0.5%) and the flows to within 1.5%. These

differences are probably due to the inaccuracies in the polynomial fit for the

fan performance curve and in the graphical solution used in the text.

Figure 9 shows a test with two fans in series (Osborne, p 76). The computed

room pressure differ from the value reported by Osborne by 2.5 Pa and the

airflows differ by less than 1%.

Figure 10 shows one floor of a 36-room airflow network created to test

execution time for a larger network. This test case describes a four story

building with six rooms, a hallway, an elevator shaft, a stairwell, and a node

representing ambient on each floor. The nodes representing the elevator shaft

and stairwell on each floor are connected by very large (2.0m2) openings.

Similar openings connect each room to the hallways. Very small (0.01m2)

openings connect the building nodes to the outside. Intermediate size

openings (0.lm2) connect the large vertical shafts to the hallways. This case

was solved, with a 0.01% convergence criterion, in 5 iterations and 2.89

seconds on a PC compatible computer (4.77 MHz 8088 CPU with 8087 math

coprocessor). This is about 16 times faster than the predecessor to AIRNET,

the AIRMOV program (Walton, 1984), could solve this sample problem.

Appendix B discusses some of these and other test cases in greater detail

A Comparison of Methods

The ESP building thermal simulation program (ABACUS, 1986) includes a

separate program, ESPAIR, for calculating airflows. ESPAIR was compared to

the AIRNET program. Both programs were recompiled and run on a workstation

computer using the 36-room test case. ESPAIR solved this case using default

5% convergence in about 8800 iterations requiring a total of 150 seconds.

AIRNET solved it using the default 0.01% convergence in 5 iterations requiring

0.16 seconds, or about 1000 times faster.

14
This extreme difference in calculation times occurs partially because of the

difficulty which the ESPAIR algorithm has with large openings (Clarke, 1985, p

206). Limiting all the openings to an area of 0.01m2 allowed ESPAIR to reach

a solution in only 137 iteration and 2.10 seconds. AIRNET was also somewhat

faster for this case: 2 iterations and 0.06 seconds, or about 35 times faster

than ESPAIR. Greater accuracy in the ESPAIR solution (0.5% convergence)

required more iterations (22,000) and more time (400 seconds). This is a

particularly extreme comparison, since it was concerned with a relatively

large problem and large interroom openings where ESPAIR is weakest and does

not consider overhead calculations such as I/O, but it does demonstrate the

potential of the new method. The ESPAIR results may explain why airflow

network calculations have a reputation for being slow.

5. DIRECTIONS FOR FUTURE RESEARCH

5.1 Alternate Solution Methods

Although the simple tests of the AIRNET program and its comparison to ESPAIR

look very promising, some important questions remain. The most important

question concerns the reliability of the method for solving the airflow

network equations. Solution of the nonlinear equations has been demonstrated

in several tests but has not been mathematically proven. The literature for

the solution of similar equations may be helpful. The airflow network is very

similar to a pipe network with the flow resistance of openings and ducts

corresponding to the resistance of pipes and fans corresponding to pumps.

Much of the theory for computing fluid flows in pipe networks is described by

Jeppson (1976). The basic flow phenomena are nonlinear and must be described

by a set of nonlinear algebraic equations. These equations may be expressed

in terms of the unknown flows in the pipes (referred to as loop equations) or

the unknown heads at the junctions (node equations). The equations are

derived from a form of Kirchoff's circuit laws: (1) the sum of flows into a

junction equals the sum of outward flows, and (2) the total headloss around

any loop in the system must be zero. Wood and Rayes (1981) give an excellent

comparison of several algorithms. Five methods are described and tested;

three are based on the loop equations and two on the node equations. The

least reliable methods (those least likely to converge to the correct

solution) are the method that adjusts each loop flow individually, the method

that adjusts each node head individually, and the method that adjusts the node

heads simultaneously.

It is interesting to note that ESPAIR solves the airflow problem with a

version of the algorithm which adjusts node heads individually. Among the

airflow algorithms used in smoke control algorithms, Klote & Fothergill (1983)

use individual node head adjustments while Sander (1974) uses the simultaneous

node head adjustment algorithm, both of which are among the least reliable

methods, according to Wood and Rayes. The method used in AIRNET also does

simultaneous node head adjustment, but it is so different that it should be

evaluated separately. It addresses the two problems observed by Wood andRayes: 

(1) Large openings (low resistances) give inexact flows because small

differences in the computed pressures lead to large differences in the flows.

This is solved by stringent requirements on mass balance convergence at each

15
node. 

Such accuracy is not costly, because Newton's method is quadratically

convergent --near the solution each iteration greatly improves the accuracy

(2) Failure of the node adjustment method to converge because of oscillating

corrections is handled by the Steffensen iteration applied to the Newton's

method correction factors.

The two simultaneous loop methods have a good history of convergence for pipe

network problems. On the other hand, they are more difficult to set up than

the node methods since independent loops must be defined; they tend to

require the solution of more simultaneous equations; the equations do not have

the very desirable feature of diagonal dominance; they tend to be less sparse

than the node equations; and some airflow elements may be difficult to

implement. The doorway model may be difficult because it can have either one

or two flows which may make it especially difficult to define the loops.

Of particular interest to the idea of establishing a general modular program

is that the loop methods require the airflow elements to compute pressure drop

as a function of flow rate, which is opposite the requirement for the node

method. For some of the airflow elements, such as powerlaw elements, the

transformation is simple. Others are described more naturally in one form

than the other. For example, the duct and fan models are described more

naturally for the loop method. This indicates the need to consider the

solution technique in the development of element models.

The work of Wood and Rayes indicates that several apparently reasonable

solution methods for the simultaneous nonlinear equations are not very

reliable. The ideal solution to the question of reliability would be

mathematical proofs of the convergent nature of the solution algorithm and the

limitations on the element models. Such proofs may be difficult to achieve

for nonlinear systems. Alternatively, extensive tests of different methods

would give some confidence as to their reliability.

5.2 Other Element Models

The modular structure of AIRNET would allow many more airflow elements to be

developed. These elements could provide either new capabilities or more

accurate simulation. The necessary requirements for element models are w(AP)

uniquely defined for all AP, and bounded derivatives for all AP.

It appears possible to represent dampers as variable flow resistance elements

in an airflow network. The relationship between resistance and actuator

position could be represented by a polynomial. The flow characteristics of

some airflow elements may depend significantly on the direction of flow. In

pipe networks check valves act in such a manner. These could be represented

by elements with separate performance curves applied to different pressure

drop or flow regimes.

Much more work could be done on the development of the doorway models.

Complex flow patterns involving boundary layer flows can exist. These

patterns are related to the geometries and surface and air temperature

distributions in the adjoining rooms. For example, Hill (1986) uses a model

which incorporates nonuniform temperatures in the rooms which leads to

multiple neutral pressure levels in the doorway and compares the computed

16
airflows to measured flows. Here the intimate relation between the doorway

element model and the node models is important. The constant temperature node

model could be expanded to three more complex models: (1) temperature varies

uniformly with height, (2) two uniform temperature layers in the room, and (3)

two layers each having uniformly varying temperature. It may be necessary to

develop several doorway models to account for different type of airflow.

Detailed doorway calculations would then involve methods to identify which

model to use.

The experimental data base for two-way flows between nodes at different

heights (through stairs and elevator shafts) appears insufficient to develop

element models. It should be possible to extend the airflow network method to

include 2- and 3-dimensional fluid elements for the detailed modeling of

airflows within rooms. Of course, this would greatly increase computation

time.

5.3 General Limitations

The simple airflow network method outlined above has some inherent limitations. These include inability to quickly model airflow patterns within a

room or to model the transient airflows caused by short-term transients in

wind pressure distributions. These effects can possibly be approximated by

dividing rooms into several nodes and adding transient flows to the average

flows, but the direct modeling of such effects would greatly increase

calculation time and would probably be impractical for most engineering

analysis. The existence of such known limitations, not to mention unknown

factors, makes experimental validation of airflow network calculations

essential.

It must be expected that uncertainty in the input parameters will always

limit the absolute accuracy of airflow calculations. However, a network model

based on physical laws will be useful for evaluating design alternatives

because relative changes in flow values should be fairly accurate. Modularity

can be used further in the design of an airflow analysis program. Figure 11

shows a structure for such a program similar to that used in ESP (ABACUS,

1986). The program separates the evaluation of wind pressures from the

airflow calculation to allow alternate inputs: manual entries, measured

values, or simulated values. Input of airflow element data would involve a

data base of element data. The computed airflows go to an output file which

could be used in either indoor air quality or loads calculations. The entire

procedure could be incorporated into an energy analysis program.

6 SUMMARY AND CONCLUSIONS

This report has discussed how an airflow network method can be used to

provide a unified model of major building airflows. Of particular importance

is the idea of modularity. ESP's modularity made the comparison test

possible. It is often very difficult to isolate a single computational

feature of a monolithic program. AIRNET includes modularity of the airflow

elements, allowing elements with greatly different flow characteristics to be

connected to the core algorithm by a common interface. More air,flow elements

could be added. The sparse matrix solution of the simultaneous equations

17
involving the Jacobian matrix allows larger systems of equations to be handled

without the full execution time penalty of using the complete matrix. By

separating the solution and matrix assembly processes, faster solution

processes could be easily substituted.

The performance of AIRNET relative to ESPAIR indicates that it is practical to

solve the flow network in detail. Solution of complex airflow networks for

the steady-state case is practical on current small computers. Solution of

the dynamic case for many timesteps is now possible. The use of small

computers will make advanced user input features available which could

significantly aid in the airflow analysis process.

Research is still needed in several areas. These include determination of the

most reliable airflow network solution method, a mathematical analysis of the

network flow equations and the solution method, development of additional

airflow elements (especially improved large opening models), experimental

validation of the simplifying assumptions in the element models and network

method, expansion of the wind pressure and airflow element performance

database, and modeling of intraroom effects by simplified methods and by

integration with microscopic modeling methods.

7 REFERENCES

ABACUS, 1986. "ESP, A Building and Plant Energy Simulation System, Version 5

Release 3", University of Strathclyde, Scotland.

ASHRAE, 1983. ASHRAE Handbook -1983 Egui~ment. ASHRAE, Atlanta GA.

ASHRAE, 1985 ASHRAE Handbook -1985 Fundamentals, ASHRAE, Atlanta GA

Axley, J., 1987. Indoor Air Quality Modelin~ Phase II Re~ort, US National

Bureau of Standards report NBSIR 87-3661.

Axley, J., 1988. "Multi-Zone Dispersal Analysis by Element Assembly",

submitted to Building and Environment.

Baker, P.H., Sharples, S., & Ward, I.C., 1987. "Air FLow Through Cracks"

Building and Environment, Vol. 22, No 4, Pergamon, pp. 293-304.

Barakat, S.A., 1987. "Inter-Zone Convective Heat Transfer in Buildings: A

Review", ASME Journal of Solar Engineering, Vol. 109, May.

Blevins, R.D., 1984. Applied Fluid D~namics Handbook, Van Nostrand Reinhold.

Brown, W.G. & Solvason, K.R., 1962. "Natural Convection Through Rectangular

Openings in Partitions -1: Vertical Partitions", International Journal of

Heat and Mass Transfer, Vol. 5, pp. 859-867.

Clarke, J.A., 1985

Bristol and Boston

Energy Simulation in Building Design. Adam Hilger Ltd

Conte, S.D. & de Boor, C., 1972.

Algorithmic A~~roach, McGraw-Hill.

Elementarx 

Numerical Analxsis: an

18
Dhatt, G., Touzot, G., & Catin, G., 1984.

Dis~layed, John Wiley & Sons, New York.

The Finite Element Method

Demuren, A.a. & Ideriah, F.J.K., 1986. "Pipe Network Analysis by Partial

Pivoting Method", Journal of Hydraulic Engineering, ASCE, Vol. 112,

pp. 327-333.

Hill, D., Kirkpatrick, A. & Burns, P., 1986. "Analysis and Measurements of

Interzone Natural Convection Heat Transfer in Buildings", Journal of Solar

Energv Engineering, ASME, Vol. 108, August.

Isaacs, L.T. & Mills, K.G., 1980. "Linear Theory Methods for Pipe NetworkAnalysis", 

Journal of the H~draulics Division. Proceedings ASCE, July,

pp. 1191-1201. (See also author's closure, Jan 1982, p. 153)

Jeppson, R.W., 1976

Ann Arbor, MI, 1976

Analysis of Flow in PiDe Networks, Ann Arbor Science

Klote, J.H. & Fothergill, J.W.Jr., 1983.

Buildings", ASHRAE, Atlanta GA.

"Qesign of Smoke Control Systems for

Liddament, M. & Thompson, C., 1982. "Mathematical models of air infiltration

-a brief review and bibliography", Technical Note AIC 9, The Air Infiltration

Centre, Bracknell, England.

Persily, A.K. & Grot, R.A., 1985. "The Airtightness of Office Building

Envelopes", Proceedings of the ASHRAE-DOE-BTECC Conference on the Thermal

Performance of the Exterior Envelo~es of Buildings III, Clearwater Beach FL

Osborne, W.C., 1977 ~, 2nd edition, Pergamon Press

Said, M.N.A., 1988.

30, No 4, April.

"A Review of Smoke Control Models", ASHRAE Journal, Vol.

Sander, D.M., 1974. "Fortran IV Program to Compute Air Infiltration in

Buildings", NRC Canada, Div. Bldg. Res., Computer Program No. 37.

Swami, M.V. & Chandra, S., 1988. "Correlations for Pressure Distribution on

Buildings and Calculation of Natural-Ventilation Airflow", ASHRAE Transactions, Vol. 94, Part 1.

Walton, G.N., 1982. "Airflow and Multiroom Thermal Analysis", ASHRAE Transactions, Vol. 88, Part 2.

Walton, G.N., 1984. "A Computer Algorithm for Predicting Infiltration and

Interroom Airflows", ASHRAE Transactions, Vol. 90, Part 1.

Wood, D.J. & Rayes, A.G., 1981. "Reliability of Algorithms fo~ Pipe Network

Analysis", Journal of the Hydraulics Division, ASCE, Vol. 107, pp. 1145-1161

19
Figure 1. Plan of a Portion of a Building Served by a VAV System

20
Figure 2. Airflow Network for VAV System of Figure 1

21
c( 0 )

Q)

~

='

cn

cn

Q)

~

~

'0

Q)

+J

='

~0

u

c( 2 )

...,.

4 iteration

c( 3 )

c( 1 )

Figure 3. Oscillating Pressure Corrections

22
with

hd = height of the base of the doorway,

~, ~ = heights of nodes nand m,

Pn, Pm = reference pressures of rooms (nodes) nand m,

Pn, Pm = air densities of rooms (nodes) nand m,

Pan' Pam = reference pressures of rooms nand m at base of doorway,

H = height of doorway, and W = width of doorway.

neutral plane

above doorway

neutral plane

below doorway

neutral plane

within doorway

Figure 4 Doorway Coordinate System and Three Flow Patterns

23
Figure 5. Typical Fan Performance Curve

24

T=18°C T=22°C

1. .2

P = 0

Doorway simulated

by doorway element

Doorway simulated

by multiple

powerlawelements

Five doorway elements in series

Figure 7 Three Airflow Networks for Doorways

26
600 2.25

1.9 2

kW

400 Pa

F2

/

0-0

200

m3/s

2.

3.0

1 2 3 4 5 6 m3/s

Figure 8. Fans in Parallel (from Osborne)

27
r.~~M~-~~~~~~~:~:_F 1 ~ 3 F 2 R 2

Room 4 6

1 2

800 5

,

I

I

I

I

I

I

I

I

I

I

~5 m3/s

j

-'

F1 -R1 600

400

200

J /1\

Atmosoheric Pressure Pa 0 ~ I I

\ 4 m3/s

\

\

\

\

3.48 \

m3/s \

! 1

0.38 m3/s

2

3.1

-200

-400

~-R2

-600

-800

Figure 9. Fans in Series (from Osborne)

28
.10 (ambient)

a = 0.01 m2 orifice

b = 0.10 m2 orifice

c = 2.00 m2 orifice

Elements:

Nodes 8 and 9 connected to similar nodes

on adjacent floors by type c openings

Figure 10. One Floor of the 4 Floor / 36 Room Timing Test Case

29
Figure 11 Proposed Structure of Airflow Network Program

30
APPENDIX A: AIRNET User's Manual

Introduction

AIRNET is a program for testing airflow network calculations which might

be used for modeling of infiltration, natural ventilation, interroom, and

mechanical system airflows. The program has been written for IBM PC' and

compatible computers using the Turbo C v 2.0.. compiler. The C source code is

included on the distribution disk to allow modification of the program for

research. Execution of the program is controlled interactively, but most of

the problem data is stored in files which are created by the user.

Interactive in~ut

There are two primary forms of interactive input.

a question of the form

The first consists of

> Question? (y/n)

to which the user can press the Y and ENTER keys for a "yes" response or N

and ENTER for a "no". The second form allows the entry of numeric data:

> enter parameter [min = Pl, max = P2, default = P3]

where the information in the brackets in the minimum permitted value, the

maximum and the default. Pressing only enter is equivalent to entering the

default value. If a default is not given, a value must be entered.

Execution of AIRNET begins with several questions about the input and

output files. The user must first provide the name of the network data file,

which contains the description of the airflow network and elements. Second,

the name of the wind pressure coefficients file is entered. Third is an

option to name an output file which will hold the results of the calculations

The names of the input files are written to the output file to help document

the simulation.

The user may then enter several run control parameters. The following

parameters are requested:

(1) "output control flag". Higher values cause more writing to the output

file: 0 gives node pressures; 1 adds element flow rates; 2 adds data on the

iterations; 3 adds an echo of the input data which can be useful in locating

errors in the input files; higher numbers can be used to dump intermediate

calculations if the program is recompiled with appropriate debugging

parameters defined.

(2) "skip initialization flag". 0 means use th~ linear initialization

method; 1 means use the current values for node pressures to begin the

Newton's method.

(3) "maximum iterations". This is the number of Newton's method iterations

tried until it is assumed the solution is not converging.

* Internation Business Machines Personal Computer -copyright/trademark

** Turbo C version 2.0 by Borland International -trademark

31
(4) "maximum pressure change (Pa)". Limit the maximum change in node

pressure during a single Newton iteration.

(5) "relative airflow convergence". The solution at a node is assumed

convergent if I ~ flows I / ~ I flows I < this factor.

(6) "absolute airflow convergence (kg/s)". The solution at a node is assumed

convergent if I ~ flows I < this factor.

(7) "convergence acceleration limit". If the ratio of successive pressure

corrections is less than this limit, use Steffensen acceleration algorithm.

Processing of the run control parameters ceases when the user responds that

all values are correct.

The user may then enter several weather values related to wind pressure.

(1) "ambient temperature (C)".

(2) "barometric pressure (Pa)". These two values are used to determine the

air density. The barometric pressure should be the absolute pressure for the

site, i.e. not corrected to sea level.

(3) "wind speed (m/s)".

(4) "wind direction". This is the direction, in degrees clockwise from

north, that the wind is blowing from.

Again, processing of the weather data ceases when the user responds that all

values are correct. See section A.4 for details on the wind pressure

calculation.

Both data files are then processed including the allocation of memory to

store the simulation data. This allows the size of the problem to be limited

only by the amount of available memory. Warnings and error messages may be

given. The user can then print information relative to the structure of the

Jacobian matrix. This will indicate how efficiently the matrix is stored for

a sparse solution. See section A.5 for more details on this subject.

The user can interactively change some of the network description data

from the input file. The user can change the values of all control values

from their default settings of 1. See the information with the airflow

elements for uses of control values. The revised network description

parameters may be printed. The user must then explicitly allow the airflow

network calculations to proceed. This allows a chance to revise any of the

above data.

After computing the airflows, the user can continue to execute AIRNET.

The program returns to the point where the run control parameters are

reviewed. The user can change these parameters, the weather data, the control

values, or some of the airflow network parameters to create a modified

problem, which can then be solved. The following network parameters can be

altered: the height, temperature, or pressure of each node; the link heights,

element, wind profile, or wind pressure modifier of each link.

A.3 The AIRNET data files

Most of the data required to simulate an airflow network is contained in

two data files. Files are used, rather than interactive input, because of the

potential size and complexity of airflow network problems. ASCII files are

used for portability.

32
A.3.1 Wind coefficients file

The wind coefficients file contains data relating the wind pressure to

the wind direction. The first line of the file is a title which is echoed to

the output. The remaining lines each consist of a profile name (of up to 15

characters) followed by 16 wind pressure modifiers corresponding to wind

directions of 0°, 22.5°,45°, 67.5°, etc. See section A.4 for details on the

wind pressure calculation. Each wind profile must have a different name.

A.3.2 Network definition file

The AIRNET network data file contains a complete description of the airflow network. The file contains node data, element data, and linkage data.

It is intended that the network data file be self documenting. Node,

element, and linkage names may each be up to 15 character long (no imbeddedblanks). 

Names should be chosen for physical significance. In addition, an

optional comment can be placed on any line after the last data entry. Empty

lines may be inserted to improve readability. An asterisk in the first column

of a line may be used to signify end-of-data; any information after that is

not processed. The first line of data is a title which is echoed on the

output listings.

In the following description keywords are shown underlined.

A.3.2.1 Node data

The nodes in the airflow network can represent either rooms or ductwork

linkage points.. Each node description has the following form:

ht name type temp pres

!lQQg.

name

ht

temp

pres

-This word identifies the following data as node data.

-The name identifies the node for later reference and in the output

listings. Each node must have a different name.

-node type (single character) (y = variable or unknown pressure, £ =

constant or known pressure, ~ = ambient node) The ambient node is a

constant pressure node whose temperature is set to ambient

temperature by the program.

-reference height (m) at which the node pressure is computed.

~ node air temperature (C).

-node pressure (Pa) (this value is not required for type v nodes).

Notes:

There must be at least one constant pressure (or ambient) node in the

network in order to have a solution to the simultaneous equations. In

addition, every node must be connected by some path to a constant pressure

node. This condition is tested by the program.

The sequence of nodes determines the sparsity of the Jacobian matrix for

the Newton's method solution and can, therefore, significantly effect execution time. The nodes can be easily reordered with a text processor. See

section A.5 for more details.

33
A.3.2.2 Element data

Each airflow element begins with an identifier (element), a name, and an

element type followed by performance parameters on one or more lines. Each

element must have a different name.

Powerlawelement:

element ill name init lam turb expt

init

lam

turb

expt

-coefficient for linear initialization (this may be the same as lam)

-coefficient for laminar flow;

-coefficient for turbulent flow;

-pressure difference exponent;

Detailed duct element (Darcy-Weisback model / Colebrook friction correlation)

~ len

lflc ldlc

element dh area rgh

init

name

tdlc

len

dh

area

rgh

tdlc

lflc

ldlc

init

-length (m);

-hydraulic diameter (m);

-cross section area (m2);

-roughness dimension (m);

-turbulent dynamic loss coefficient;

-laminar friction loss coefficient;

-laminar dynamic loss coefficient;

-laminar initialization coefficient.

Doorway element:

name Q:Q.!:

dtmin ht

init lam

wd cd

element turb expt

init -expt: powerlaw coefficients for low temperature difference;

dtmin -minimum temperature difference for two-way flow (C);

ht -height of doorway (m);

wd -width of doorway (m);

cd -doorway discharge coefficient.

Constant flow element

element name ill flow

flow -rated mass flow (kg/s);

Control parameter = actual mass flow / rated mass flow

Note that since flow is not a function of pressure, constant flow elements

cannot be counted as part of the linkage of variable pressure nodes to

constant pressure nodes (A.3.2.1).

34
Detailed fan element

element name

rdens

all

a2l

init -expt: powerlaw coefficients for fan speed = 0;

rdens -reference fluid (air) density (kg/m3);

fdf -free delivery flow (kg/s);

sop -shutoff pressure (Pa);

ltt -laminar/turbulent transition (RPM/rated RPM);

nr -number of ranges for performance curve;

mfl -minimum mass flow of range 1 (kg/s);

all-a14: range 1 polynomial coefficients for f(p) at rated RPM;

mf2 -maximum mass flow rate of range 1 (kg/s) =

minimum mass flow of range 2;a2l-a24: 

range 2 polynomial coefficients for f(p) at rated RPM;

mf2 -maximum mass flow rate of range 2 (kg/s);

mfn -maximum mass flow rate of range nr (kg/s);

Control parameter = actual RPM / rated RPM

Constant power fan element:

element name £Ri upo

upo -useful power output of fan (W);

Control parameter = actual power / rated power

Quadratic flow element:

element name ill a b

-flow coefficient;

-flow2 coefficient

a

b

A.3.2.3 Linkage data

The linkage data defines the airflow network. Each linkage description

consists of the identifier (link) followed by six or seven parameters.

link. name node-l ht-l node-2 ht-2 element wind wpmod

name

node-l

ht-l

node-2

-This user assigned name identifies the link in the output listings.

-name of the first node. This refers to the names which were

assigned with the node data.

-height (m) of the linkage point relative to the first node height.

-name of the second node. Node-l and node-2 define the direction of

airflow: flow from node-l to node-2 is positive, the reverse flow

35
I is negative.

ht-2 -height (m) of the linkage point relative to the second node height.

element -airflow element name. This refers to the names which were assigned

with the element data.

wind I -wind pressure coefficients profile. This is the name of the profile

given in the wind pressure coefficients file. The name "null" is

reserved to mean that there is no wind pressure on a given link.

wpmod -This wind pressure modifier is applied for the particular link.

See section A.4 for details.

Note that the wind coefficient profile name must exist in the wind

coefficients file and that the node and element names must be defined in the

network file before they are referenced by a link.

Wind Pressure Calculations

The wind pressure, PWi' acting on link i is proportional to the kinetic

energy of the airstream: ~PaVa2. The ambient air density, Pa, is determined

from the weather data parameters: temperature and barometric pressure. The

wind velocity, Va, is also a weather data parameter. The wind pressure is

also a function of the wind direction, D. D is a weather data parameter; the

function is given in the wind coefficients file in terms of 16 wind pressure

modifiers for evenly spaced wind directions. The wind direction pressure

modifier, f(D), is determined from the wind coefficients by linear interpolation. Wind speed, and therefore pressure, also varies with height and local

shielding of the surface. These variations can be included in the link wind

pressure modifier, wpmod. Therefore, the wind pressure acting at a given link

is

PWi = wprnod .f(D) .~PaVa2 (A.l)

This procedure should give considerable flexibility in defining the wind

pressure, but it does not say what coefficients should be used. That is left

to the user.

One particular aspect of the input should also be noted. A positive wind

pressure tends to increase flow in the positive direction and vice versa for a

negative value of wind pressure. Flow direction is determined by the ordering

of the two nodes in the link command. Therefore, in order for positive wind

pressures to push air into the building and negative wind pressures to pull

air out, each link which represents an opening in the envelope of the building

should be defined with the ambient node first and the interior node second.

Equation Ordering in the Jacobian

The sequence in which the nodes are listed in the network data file

determines the ordering of the simultaneous equations and pattern in which the

Jacobian is filled. Different equation ordering can significantly effect the

sparsity of the Jacobian and, therefore, the memory requirements and execution

time of the program. Memory requirements increase linearly with the number of

nonzeros in the Jacobian; computation time increases more than linearly. See

chapter 2 of Sparse Matrix Technolog~ (Academic Press, 1984) by S. Pissanetzky

for a detailed discussion of this topic.

36
Consider the following example. Nodes B, C, D, and E are linked to node

A but not to each other. If only node A is constant pressure and are ordered

A ~ E, the diagram of the upper part of the symmetric Jacobian will look like

(using. to indicate initial nonzero values, 0 to indicate zeros, and + for

new nonzero values after factoring):

.0 0 0

.0 0

.0

.+ + +

.+ +

.+

.

before and after factoring, respectively, showing that this numbering scheme

results in non-zero elements, or "fill", below the diagonal.

If the nodes are ordered B ~ E, A, the diagram of the Jacobian will look

, .,

.000.

.0 0 .

.0 .

..

.

both before and after factoring, thus completely avoiding fill. See section

B.1O for another example of equation reordering.

To reduce the computations associated with fill, a rule of thumb for the

ordering of nodes is that the node which is linked to the greatest number of

other nodes should appear later in the node list. In buildings consisting of

many similar levels this is partially accomplished by grouping all nodes on

each level together.

Sample Out~ut File

The following output file for a simple test network illustrates its

contents by notes which are indicated in brackets: [..].

[1 powerlaw test input file

wind coefficients file -no profiles

Run control data:

output control flag:

skip initialization flag:

maximum iterations:

maximum pressure change:

relative airflow convergence:

absolute airflow convergence:

convergence acceleration limit:

[2]

3

0

20

500 Pa

0.0001

1e-06 kg/s

-0.5

Weather data:

ambient temperature:

barometric pressure:

wind speed:

wind direction:

20 C

101325 Pa

0 m/s

0 (N=O, E=90,

37
node

node

node

node

node-l

node-2

node-3

node-4

0.0

0.0

0.0

0.0

20.0

20,0

20.0

20.0

c 50.0

v

v

c -50.0

element

element

arf-O.OOO! 

plr 7.2e-9

arf-O.OI00 plr 7.2e-6

7.2e-9

7.2e-6

8.48528e-5

0.00848528

0.5

0.5

orf-O.OOOl m"2

orf-O.Ol m"2

link

link

link

link-l

link-2

link-3

node-l

node-2

node-3

0.0

0.0

0.0

node-2

node-3

node-4

0.0

0.0

0.0

orf-O.OOOl

orf-O.O1OO

orf-O.OOOl

null

null

null

*********

Time to read data files 0.44 [4

1234

1: *2: 

**3: 

*4: 

*

1234

Initial fill fraction: 0.500000

Final fill fraction 0.500000

Number of nonzero elements in AU: 2

Number of nonzero elements in AL: 2

Unallocated memory: 54212 bytes

dens

1.204742

1.204148

1.204148

1.203554

wp

0.000000

0.000000

0.000000

visc

0.000018

0.000018

0.000018

0.000018

dp

0.000000

0.000000

0.000000

-1.2647e-O2 -S.OOOOe+Ol [8

1.000000e+OO

1.000000e+OO

3.170580e-O3

2.151550e-O2

[9

5.502771e-Ol

5.501204e-Ol

1.852654e-O2

6. 146604e-O3

node temp

node-l 20.000000

node-2 20.000000

node-3 20.000000

node-4 20.000000

link sp

link-l 0.000000

link-2 0.000000

link-3 0.000000

InitializationPn: 

5.0000e+Ol 3.7316e-02

Begin iteration 1Rev: 

2 1.422954e-03 3.414537e-02Rev: 

3 -1.422954e-03 -3.416223e-02

Begin iteration 2

Rev: 2 -1.919686e-03 -2.790587e-02

Rev: 3 1.919686e-03 2.793733e-02

Begin iteration 3Rev: 

2 3.77570ge-04 4.51476ge-03Rev: 

3 -3.77570ge-04 -4.508156e-03

1.000000e+OO

1.000000e+OO

1.401177e-O2

1.065476e-O2

38
1.000000e+OO

1.000000e+OO

1.475415e-O2

9.916436e-O3

1.000000e+OO

1.000000e+OO

1.483198e-O2

9.832351e-O3

Begin iteration 4Rev: 

2 -1.189833e-O4 -7.423793e-O4Rev: 

3 1.189833e-O4 7.383248e-O4

Begin iteration 5

Rev: 2 -1.083831e-O5 -7.782953e-O5Rev: 

3 1.083835e-O5 8.408471e-O5

Begin iteration 6

Fin: 1 O.OOOOOOe+OO O.OOOOOOe+OO

Fin: 2 -8.919034e-O8 1.316844e-O3Fin: 

3 8.91478ge-O8 1.316844e-O3Fin: 

4 O.OOOOOOe+OO O.OOOOOOe+OO

1.204742e+OO 5.000000e+Ol

1.204148e+OO 1.483198e-O2

1.204148e+OO 9.832351e-O3

1.203554e+OO -5.000000e+Ol

Time to compute node pressures:

Number of iterations: 6

0.66

Link:link-!:

link-2:

link-3:

i

1

2

3

flow-l

6.584668e-O4

6.583776e-O4

6.584668e-O4

n

1

2

3

m

2

3

4

e

1

2

1

flow-2

O.OOOOOOe+OO

O.OOOOOOe+OO

O.OOOOOOe+OO

[13]

pdrop

4.998517e+Ol

4.99962ge-O3

5.000983e+Ol

no

no

no

no

n pres sUInf

1 5.000000e+Ol -6.584668e-O4

2 1.483198e-O2 8.917414e-O8

3 9.832351e-O3 -8.917414e-O8

4 -5.000000e+Ol 6.584668e-O4

Notes:

[1] This is the echo of the input file titles.

[2] These are the values of the run control parameters and weather data.

[3] This is the echo of the input files as set by the output control flag.

[4] The input processing time is accurate to ~0.02 second.

[5] This sketch shows the structure of the Jacobian. *'s indicate array

positions filled by the airflow element functions; +'s indicate positions

filled when the matrix is factored (See Appendix C.1). The two fill fractions

compare the positions in the sparse matrix to a full square matrix. AU and AL

refer to the upper and lower triangular matrices.

[6] Memory is allocated to store the problem parameters and variables and the

various arrays. The unallocated memory indicates that larger problems can be

accommodated.

Notes 7 through 11 are for output generated by an output control flag >= 2.

[7] The air temperature, density, and viscosity are reported for each node.

Nodes are identified by name and number, which is the sequence in which they

occur in the input file. The stack pressure, wind pressure, and total

pressure drop are reported for each link.

[8] These are the pressures computed by linear initialization and which are

used to start the Newton's method iterations.

[9] At each iteration 4 values are printed for each variable pressure node:

the sum of the flows (kg/s) into the node, the computed pressure correction

factor, a correction acceleration factor (1.0 = no correction), and the new

estimate of node pressure.

[10] On the second iteration oscillating pressure corrections have been

observed for both nodes and correction acceleration factors computed.

39

Nodededede de:-1:-2:-3:-4:
[11] When the convergence criteria have been met, 4 values are printed for

all nodes: ~ flows, ~ Iflowsl, air density, and node pressure. (Flows are

not computed for the constant pressure nodes.)

[12] The time to compute air properties and node pressures is accurate to

:to.O2 second.

[13] The airflows through each link are reported if the output control flag

>= 1. The following values are printed for each link: the link name, link

number, 1st node number, 2nd node number, element number, pressure drop across

the link, mass flow rate (kg/s), 2nd mass flow rate (possible for doorways).

[14] The final node values are reported if the output control flag >= O. The

following values are printed for each node: the node name, node number, the

node pressure (Pa) relative to ambient, and the sum of the flows into the

node and the sum of the flows ignoring direction of flow.

A.7 Compiling AIRNET

A.7.1 Distribution diskette

The distribution diskette is a standard 360K byte floppy disk containing

the following files:

size

file (k bytes) description

AIRNET.EXE 58 program executable on a PC compatible microcomputer

COMSEP.C 10 compilable code for the preprocessor program

AIRNET.ALL 126 all header files and functions for the AIRNET program

ELEMENT.EXE 56 program to help generate element models

ELEMENT.STR 6 file used by ELEMENT.EXE to limit program length

ELEMENT.ALL 59 all header files and functions for the ELEMENT program

NOWIND wind coefficients file containing no wind profiles

WIND wind coefficients file with some wind profiles

AFDATA.PL1 input data for powerlaw element test #1

AFDATA.PL2 input data for powerlaw element test #2

AFDATA.PL3 input data for powerlaw element test #3

AFDATA.ST1 input data for stack effect test #1

AFDATA.ST2 input data for stack effect test #2

AFDATA.WP1 input data for wind pressure test #1

AFDATA.WP2 input data for wind pressure test #2

AFDATA.DW1 input data for duct element test #1

AFDATA.DW2 input data for duct element test #2

AFDATA.DR1 input data for doorway element test #1

AFDATA.CF1 input data for constant flow element test #1

AFDATA.CF2 input data for constant flow element test #2

AFDATA.FN1 input data for fan element test #1

AFDATA.FN2 input data for fan element test #2

AFDATA.FN3 input data for fan element test #3

AFDATA.QF1 input data for quadratic flow element test #1

AFDATA.QF2 input data for quadratic flow element test #2

AFDATA.08F input data for 8-floor execution time test

40
A.7.2 The COMSEP program

The COMSEP program is used to assist in the compilation of the AIRNET

and ELEMENT programs. COMSEP combines and separates files involved in a

project and creates several useful project files. First, compile COMSEP.

COMSEP automatically creates files names based on the project name entered by

the user. COMSEP presents the user several possible actions to be performed.

(1) "combine subfiles":

them into a single file.

single line of the form:

This option takes a list of filenames and combines

Each of the individual files must be headed by a

/*subfile filename ******/

Note that this line is in the form of a C comment and does not effect

compilation of the file.

(2) "separate into subfiles"; This option takes a single file and separates

it into individual files based on header lines indicated above. AIRNET.ALL

must be separated for compilation. This process creates the list of subfiles

which is used with the other processes. Subfile names may not be duplicated.

(3) "create compile batch file": This option creates a batch file for use

with the command-line version of Turbo C, TCC.EXE. Be sure to include the

".bat" extension in the name of the batch file. The compile options used to

create AIRNET.EXE were "-ms -c -d -f87 -G -0 -ZIt. Other options may be

appropriate for a particular computer or problem. In particular, the -ml

large memory option will allow the program to use the entire memory of the

computer for solving large problems. Large problems almost require the use of

a math coprocessor.

(4) "create link response file": This option creates a response file for use

with TLINK.EXE to link the object files and libraries. A particular directory

structure is assumed (see below), but the link response file can be edited to

suit the user's individual needs. AIRNET.EXE was linked with the small memory

(s) and math emulation options. Emulation uses the coprocessor if one ispresent. 

The coprocessor greatly improves execution time. Compiling with the

math coprocessor option reduces the program size by about lOk bytes. It does

not significantly improve execution speed, but will cause the program to fail

if a math coprosessor is not present. Therefore, the emulation option was

chosen for the distribution version of AIRNET.EXE.

(5) "create library response file": This option creates a response file for

use with TLIB.EXE for the creation of an object library.

(0) "exit": Exit the preprocessor program

41
A.7.3 Compilation

The following directory structure is assumed:

...\TURBOC contains the Turbo C programs;

...\TURBOC\INC contains the include files;

...\TURBOC\LIB contains the libraries and COx.OBJ files;

...\TURBOC\AIRNET contains the AIRNET program files.

The file TURBOC. CFG must be appropriately defined and \TURBOC must be set in

the PATH command.

Use COMSEP to split AIRNET.ALL into subfiles, create a compile file, and

create a link response file. AIRNET.ALL contains three header files. Copy

these files to the \INC directory. Execute the compile batch program to

create object files. Use the linker and link response file to create an

executable version of AIRNET.

The header file AIRNET.H includes several definable parameters (TESTI,

TESTC, and TESTS) which were used to test the program during development.

When these parameters are not defined, the dump functions (dumps.c, dumpv.c,

dumpa.c, dumpf.c, and dumpn.c) are not called and can be removed from the link

response file to reduce the size of AIRNET.EXE. These functions increase the

size of the program by about 5k bytes. They are not included in AIRNET.EXE on

the distribution diskette.

A.7.4 The ELEMENT Program

The ELEMENT program is an aid for the creation of airflow elements for

the AIRNET program. All input is interactive; an output file will include

element descriptions that can be transferred to the AIRNET network data file

by a word processor. Both ELEMENT.EXE and ELEMENT.STR must be in the default

directory when ELEMENT is executed. The compilation process for ELEMENT.ALL

is the same as for AIRNET.ALL.

1 :

2'

3:

4:

5:

6:

7:

8:

9:10:

There are presently ten capabilities:

convert the physical characteristics of an orifice to a Rlx element;

convert the physical characteristics of a crack to a Rlx element;

convert one or two pressure/flow test points to a Rlx element; i

convert the physical characteristics of a duct to a Rlx element; !

convert the physical characteristics of a duct to a ~ element;

convert the physical characteristics of a doorway to a QQI element;

convert the performance data of a fan to a ~ element;

convert the physical characteristics of a crack to a ~ element;!

convert two pressure/flow test points to a ~ element;

convert the physical characteristics of a duct to a ~ element. !

42
APPENDIX B: AIRNET Validation Tests

B.l Introduction

Validation is an essential part of the development of a computer

simulation. Inaccurate results are not always due to program errors as shown

in the SERI report on validation of building energy analysis programs (1983)

which identified seven error sources classified into two groups. External

sources are those which are not under the control of the developer of the

computer code. These errors include:

1. differences between the actual weather around the building and the weather

used in the simulation;

2. differences between the actual effect of occupant behavior and those

effects assumed by the user;

3. user error, including inappropriate simplifying assumptions, in deriving

the input files; and

4. differences between the actual thermal and physical properties of the

building and those input by the user.

Internal error sources are those contained within the coding of theprogram. 

They include:

1. differences between the actual heat/mass transfer mechanisms and the

algorithmic representations of those mechanisms;

2. differences between the actual interactions of heat/mass transfer

mechanisms and those interactions between the algorithms; and

3. coding errors.

Three types of tests have been used to validate building energy analysisprograms. 

One is comparison to other simulation programs. Another is

comparison to analytically calculated results. The third is comparison to

experimental data. The following table from the SERI validation report

summarizes the advantages and disadvantages of each method.

VALIDATION TECHNIQUES

Technique Advantages Disadvantages 1 

Comparative No input uncertainty No truth standard

Relative test Any level of complexity

of different Inexpensive

programs Many comparisons possible f.

Analytical

Test of

numerical

solution

No input uncertainty

Exact truth standard

given the simplicity

of the model

Inexpensive

Does not test the model

Limited to cases for

which analytical

solutions can be

derived

Empirical Approximate truth Measurement involves some

Comparison to standard within accuracy input uncertainty

measured building of data acquisition High quality, detailed

performance Any level of complexity measurements are time

consuming & expensive

This report will concentrate on analytic validation tests which are best

for revealing coding errors.

43
B.l.l Simplified Analysis of Airflow Elements in Parallel and in Series

It will be necessary to develop analytic solutions for simple airflow

networks to test the AIRNET program. The analysis of DC electric networks is

simple because the elements are linear. The following discussion will show a

similar means of network analysis for nonlinear airflow elements.

(a) Laminar flow: w = (pip.) K AP

Parallel

Flows through elements 1 and 2 are given by:

Wl = (p/jJ.) Kl (Pl-PZ) W2 = (pljJ.) K2 (Pl-P2)

Wl + W2 = (p/J1.) (Kl+K2) (Pl-P2)

Wl + W2 = (pip.) Ke (Pl-P2) where Ke = Kl + K2

I Ke=IKi I

Therefore, for elements in parallel:

as long as p and ~ are constant

(B .1.1)

Series:

.Kl

Pl P3

Pressure drops through elements 1 and 2 are given by:

Pi-P2 = (~/p) Wi / Ki P2-PS = (jJ/p) W2 / K2

Pl-P3 = (Pl-P2)+(P2-P3) = (Jl./p) (wl/Kl + w2/K2)

Since Wi = Wz = W

Pl-P3 = (~/p) w (l/Kl + 1/K2

or

W = (pIJL) Ke (Pl-P3) where liKe = l/K! + l/K2

I liKe = I l/Ki I

(B.l.2) Therefore, for elements in series:

as long as p and ~ are constant

These relationships are, of course, identical to the relationships for

electric conductances.

w = jp C jAi' (b) Turbulent flow:

Parallel:

44
Flows through elements 1 and 2 are given by

Wl = Jp C1 !Pl-P2 Wz = Jp Cz J~

Wl + W2 = jp (C1+C2) !Pl-P2

w. + W2 = Jp Ce J~ where Ce = Cl + Cz

Therefore, for elements in parallel: I Ce = L Ci I

as long as p is constant (B.l.3)

Series:

.C1

Pl

Pressure drops through elements 1 and 2 are given by:

(1/P)(Wl/C1)2 P2-P3 = (1/p)(w2/C2)2

(Pl-PZ)+(PZ-P3) = (1/P)[(Wl/C1)Z + (WZ/CZ)Z

Since Wl = W2 = W

= (lip) w2(l/C12 + l/C22

or

W := Jp Ce J~ where 1/Ce2 = 1/C12 + 1/C22

Therefore, for elements in series: Il/Ce2 = L 1/Ci2 I

as long as p is constant

(B.l.4)

These relationships will allow the creation of relatively complex airflow

networks whose analytic solutions can be used to test AIRNET. The test cases

will be limited to constant p and~, that is, no temperature or height

differences in the networks. Slightly different versions of the turbu~ent

relationships (equations B.l.3 and B.l.4) appear in Klote (1983, pp 331-36).

B.l.2 Thermal Properties of Air

Two primary values are known for each node n: the temperature, Tn. and

the pressure, Pn, relative to the absolute ambient pressure. PB. From these

values the nodal density and viscosity are computed:

(B.l.S) Pn = 0.00348380 (PB+Pn)j(Tn+273.15)

and

(B .1. 6) ~n = 0.0000171432 + 0.00000004828oTn

where

Tn is in °C and PB and Pn are in Pa

It would be relatively easy to substitute the proper psychrometric relationships to compute these values more accurately and include humidity as a factor

in air density.

45
B.2 Powerlaw Element Tests

B.2.! Powerlaw Element Test #1

This is a test of the powerlaw element calculations in both the turbulent

and laminar regimes and of solution convergence at the laminar-turbulent

transition.

The air in all nodes is at 20°C; the ambient pressure is 101325 Pa

Pn = 0.00348380 (PB+Pn)/(Tn+273.15) = 1.20415 kg/m3 (Pn = 0)

~n = 0.0000171432 + 0.000000048280Tn = 0.0000181088 kg/ms (B. 2 .1)

(B.2.2)

The coefficients for the powerlaw elements used in this test are computed

from an orifice model. The critical properties of an orifice are its area, A,

hydraulic diameter, D, and discharge coefficient, c. The orifice flow

equation is

Q = c A JUPjp (B. 2 .3)

which can be rearranged to the powerlaw element equation

w (= pQ) = C J~ (B.2.4)

by letting

I c = c A J21 (B. 2.5)

Assuming that the flow is laminar below a certain Reynolds number (Re = pVD/~

= wD/~A) and that the laminar and turbulent flow equations give the same flow

at that point, gives

W = Jl,ARe/D = pKAP/Jl, and w = p.ARe /D = cAj2PAi>

K = 1J.2ARe/pDAP and A.P = J'2Re2/2pc2D2

I K = 2ADc2/Re I

(B.2.6)

The first powerlaw element assumes D = 0.1 m, A = 0.01 m2, c = 0.6, and

transition at a Reynolds number of 100. These values give Cl = 0.00848528,

and Kl = 0.0000072. Laminar flow occurs if wl (= ~ARe/D) S 0.000181088 kg/so

The second powerlaw element assumes D = 0.2 m, A = 0.04 m2, and the same c andRe. 

This gives C2 = 0.03394113, and K2 = 0.0000576. Laminar flow occurs if

W2 S 0.000362176 kg/so

In this test the two powerlaw airflow elements are arranged in series:

.C1

P,

.

P2 P3

46
For turbulent flow in both elements, the equivalent powerlaw coefficient is

Ce = [ 1/C12 + 1/C22 -\ = 0.00823193

(B.2.7)

For laminar flow in both elements, the equivalent coefficient is

K = e

ilK! 

+ i/K2 -1 = 0.0000064

(B. 2 .8)

Substituting these two coefficients and thermal properties into the appropriate expressions give

WI = o. 00903321 j~ and WL = 0.425570 (Pl-P3) (B.2.9)

Flow is turbulent in both elements if Pl-P3 ~ 0.00160752 Fa. Flow is laminar

in both elements if Pl-P3 ~ 0.000042552 Fa. AIRNET should give the following

flow values as a function of P1 (assuming P3 = 0):

Pl (Pa)

1.000

0.100

0.010

0.00161

0.0016076

0.0000426

0.000042

0.00004

w (kg/s)

0.00903321

0.00285655

0.000903321

0.000362456

0.000362541

0.000018108

0.000017874

0.000017023

Notes:

(1) AIRNET computed the expected flows in 4 or fewer iterations.

(2) If two equal elements were used in the airflow network, the AIRNET

initialization algorithm would estimate P2 = (P1+P3)/2, which is exactly

correct. The iterative solution process would not be tested.

(3) When flow is laminar in both elements, the initialization gives the exact

solution.

The followinng input data is include as file AFDATA.PLl on the AIRNET

distribution diskette.

powerlaw test #1 input file

node

node

node

node-l

node-2

node-3

0.0

0.0

0.0

20.0

20.0

20.0

c 1.0

v

c 0.0

orf-O.O!

orf-O.O4

plr

plr

orifice -0.01 m"2

orifice -0.04 m"2

element

element

7.2e-6 7.2e-6 0.00848528

5.76e-5 5.76e-5 0.03394113

0.5

0.5

link

link

link-l

link-2

node-l

node-2

0.0

0.0

node-2

node-3

0.0

0.0

orf-O.Ol

orf-O.O4

null

null

*********

47
B.2.2 Powerlaw Element Test #2

This is a test of the AIRNET convergence acceleration algorithm. It

involves a simple case which can be very difficult to solve. In this test

three powerlaw airflow elements are arranged in series:

.Ci

Pi

.

P4

When C2 is very large compared to C1 and C3, as in the case of a doorway

connecting two rooms which have only small cracks connecting to ambient,

convergence can be very slow with some algorithms.

This test uses a series of different size orifices, each with c = 0.6 and

D2 = A, varying in size from 0.0001 m2 to 100.0 m2. The computed flows can be

checked by the equivalent powerlaw coefficient (turbulent flow) which is

c = 1/C12 + 1/C22 + 1/C32 -~

e (B.2.10)

AIRNET should give the following flow values (Pl-P4 = 100 Pa) for each combination of three elements listed. The numbers of iterations are the results of

the test.

w (kg/s)

0.00053758

0.00065676

0.00065839

0.00065840

0.00065840

0.00065840*

0.00065840*

Iterations

2 (2)

5 (5)

8 (10)

8 (29)

7 (83)

5 (>100)

5 (>100)

Elements

1-1-1

1-2-1

1-3-1

1-4-11-5-1

1-6-1

1-7-1

w (kg/s)

0.00092195

0.00537585

0.00656764

0.00658388

0.00658404

0.00658404

0.00658404*

Elements

2-1-2

2-2-2

2-3-2

2-4-2

2-5-2

2-6-2

2-7-2

Iterations

5 (5)

2 (2)

5 (5)

8 (10)

9 (29)

11 (93)

9 (>100)

Notes:

(1) AIRNET computed the expected flows.

(2) * indicates cases with laminar flow in the center element.

(3) The iterations in parnetheses are without convergence acceleration.

powerlaw test #2 input file

orf -0.0001m"2

orf -0.001 m"2

orf -0.01 m"2

orf -0.1 m"2

orf -1.0 m"2

orf -10.0 m"2

orf -100.0 m"2

node

node

node

node

element

element

element

element

element

element

element

link

link

link I

8.48528e-5

0.000848528

0.00848528

0.0848528

0.848528

8.48528

84.8528

orf-O.OOOl

orf-O.OOOl

orf-O.OOOl

0.5

0.5

0.5

0.5

0.5

0.5

0.5

null

null

null

node-l 

c 0.0 20.0 50.0node-2 

v 0.0 20.0node-3 

v 0.0 20.0

node-4 c 0.0 20.0 -50.0

orf-O.OOOl plr 7.2e-9 7.2e-9

orf-O.OOl plr 2.276ge-7 2.276ge-7

orf-O.Ol plr 7.2e-6 7.2e-6

orf-O.l plr 2.276ge-4 2.276ge-4

orf-l.O plr 0.0072 0.0072

orf-10. plr 0.22769 0.22769

orf-100. plr 7.2 7.2

link-l node-l 0.0 node-2 0.0

link-2 node-2 0.0 node-3 0.0

link-3 node-3 0.0 node-4 0.0

*********

48
Powerlaw Element Test #3

This is a test of the convergence of the AIRNET iterative solution method

for a relatively complex network of airflow elements. It involves 12 nodes

and 20 powerlaw airflow elements arranged in series and parallel:

Using the conversions for parallel airflow elements gives:

.Cb:

P3

:. 

Cc

P4

-====Cg ====- ~vl 0 ===-===Cl1 ===-===C12 ===-===C13 ===-===C1 4

Pl P2 Ps P6 P7 Pa

.-=C15=8

Pll P12

-.Cd .C2 0

Pg P1O

"'16

where

Ca = Cl + C2 + C3'

Cc = C6 + C7 + Cs, and

Cb = C4 + Cs,

Cd = C17 + CiS + Cig.

the conversions for series airflow elements gives:

:::::::: 

] .- c -

Ir -15- .

IIP11 P12

where

2 + 1/C14 2 -\ and

,

1/Ca2 + 1/Cb2 + 1/Cc2 ]-~,

1/Cl02 + 1/Cl12 + 1/C122 + 1/C13

1/C162 + 1/Cd2 + 1/C2o2 ]-~

c = e

Cf =

C =

&

Again using the conversion for parallel airflow elements gives:

-==Cg

Pl

:.

p Ch

2

:8=C15

Pll P12 where

Ch = Ce + Cf + Cg

49
Again using the conversion for series airflow elements gives:

.P12 P1 .

J. where

Ci = 1/Cg2 + 1/Ch2 + 1/C15 2 -%

In order to test AIRNET convergence, powerlaw coefficient values are

created for orifices of greatly differning areas:

number c D(rn) A(rn2) C K

1 .6 0.1 0.01 0.00848528 7.2e-6

2 .6 1.0 1.0 0.848528 0.0072

3 r .6 2.0 4.0 3.394113 0.0576

4 .6 0.05 0.0025 0.00212132 9.0e-7

5 .6 0.06 0.0036 0.00305470 1.555e-6

6 .6 1.0 1.0 0.848528 0.0072

7 .6 1.0 1.0 0.848528 0.0072

8 .6 1.0 1.0 0.848528 0.0072

9 .6 1.0 1.0 0.848528 0.0072

10 .6 1.0 1.0 0.848528 0.0072

11 .6 0.01 0.0001 0.000084853 7.2e-9

12 .6 1.0 1.0 0.848528 0.0072

13 .6 0.02 0.0004 0.000339411 5.76e-8

14 .6 2.0 4.0 3.394113 0.0576

15 .6 1.0 1.0 0.848528 0.0072

16 .6 0.02 0.0004 0.000339411 5.76e-8

17 .6 2.0 4.0 3.394113 0.0576

18 .6 0.01 0.0001 0.000084853 7.2e-9

19 .6 1.0 1.0 0.848528 0.0072

20 .6 0.03 0.0009 0.000763675 1.944e-7

which gives w = 0.0611024 kgjs for Pl-P12 = 100 Pa and p = 1.20415 kgjm3.

AIRNET computed the expected flow in 12 iterations. This indicates that

the number of iterations increases with the complexity of the network.

However, this large number of iterations is for a combination of large and

small airflow resistances chosen to be difficult to solve. The solution

required 157 iterations without convergence acceleration. AIRNET solved the

same network containing 20 identical airflow coefficients (C = .00848528) in

only 6 iterations. Tests on different computers and with different compilers

indicate that the number of iterations for this network configuration is very

sensitive to round-off errors.

See the input file AFDATA.PL3 on the diskette.

50

The coefficients for the simplified networks are:

Ca = 4.251126

Cb = 0.00517602

Cc = 2.545584

Cd = 4.242726

Ce = b.00517601

Cf = 0.00008232

Cg = 0.00031016

~ = 0.00556848

C, = 0.00556824
B.3 Stack Effect Tests

Once the performance of the powerlaw element model has been verified, it

can be used to test the stack pressure and wind pressure calculations.

B.3.1 Stack Effect Test #1

The following figure shows two rooms represented by airflow nodes 1 and 2

connected by two simple openings at different heights.

-z = 10m

z = 0

node 1

.

Tl = OC

node 2

.

T2 = 20C

The standard equation for stack effect (Klote, 1983, p12) is

h.P = .gr. ( .!. -.!.) h

R Ti T2 (1)

where

g = acceleration of gravity (9.80 m/s2),

P = absolute pressure (101325. Fa),

R = gas constant (l/R = 0.0034838),

Ti = absolute temperature of node i (OK), and

h = distance from the neutral plane (m).

The values in parentheses are the numbers used in AIRNET. They give gP/R =

3459.36 (Klote gives 3460 for this value). Since the openings are identical,

the height of the neutral plane should be 5 m. Substituting into equation (1)

gives ~P = 4.3202 Fa.

AIRNET gives slightly different results for ~P at the two openings

because the relationship between flow rate and pressure drop depends on the

direction of flow. The computed flow rate is 0.019692 kg/s at pressure drops

of -4.472866 and 4.167554 Pa for the upper and lower openings, respectively.

The average pressure drop (neglecting sign) through the two openings is 4.3202

Fa. An extra iteration produces a very small change in these results.

stack test #1 input file

0.0

0.0

0.0

20.0

node-l 0.0

node-2

node

node

c

v

orf-O.Ol 7.2e-6 7.2e-6 0.00848528 0.5 orf -0.01 rn"2 element plr

node-2

node-2

0.0

10.0

orf-O.Ol

orf-O.Ol

null

null

link-l

link-2

node-l

node-l

0.0

10.0

link

link

*********

51
B.3.2 Stack Effect Test #2

The following figure shows two rooms represented by airflow nodes 1 and 2

connected by three airflow elements which connect to node 2 at differentheights.

-t z = 5m1 

z = 0

-z = -5

node 1

.

node 2

.

: 

node 1

: .

node 2

.

: node 1

: .

node 2

.

link 1 link 2 link 3

The pressure of node 1 is 5 Pa greater than node 2. The temperature of node 1

is 0 C; node 2 is 20 C. Because of this temperature difference, the flow

through element 2 will ge less than the flow through element 1, and the flow

through element 3 will be greater ( w2 < wl < W3). The standard stack effect

equation indicates a pressure difference of about I4.3202 Pa due to the height

changes. That is, ~P1 = 5.0 Fa, ~P2 ~ 0.68 Fa, and ~P3 ~ 9.32 Fa.

The first AIRNET solution is computed with unknown airflow directions.

This gives: ~Pl = 5.0 Fa, ~P2 ~ 2.84 Fa, and ~P3 ~ 7.16 Fa. Recalculation

with known flow directions gives the expected results. See the discussion

about stack effect calculation in Appendix C.

Also note that AIRNET "solved" a case where ill the node pressures are

specified. Computer algorithms often run into numerical problems on such

trivial cases --they provide good tests for algorithmic errors.

stack test #2 input file

node

node

node-l

node-2

0.0

0.0

0.0

20.0

5.0

0.0

c

c

element orf-O.Ol plr 7.2e-6 7.2e-6 0.00848528 0.5 orf -0.01 m"2

link

link

link

link-l

link-2

link-3

node-l

node-l

node-l

0.0

0.0

0.0

node-2

node-2

node-2

0.0

5.0

-5.0

orf-O.Ol

orf-O.Ol

orf-O.Ol

null

null

null

*********

52
Wind Pressure Tests

AIRNET does not include a model for converting environmental variables

such as wind speed and direction to wind pressure. It does include a flexible

way for the user to transfer data from a wind pressure model (or measurements)

to the appropriate airflow element. The following tests are designed to

insure that the AIRNET portion of that data transfer occurs correctly.

The wind pressure tests require use of the wind coefficients file named

The contents of this file are:

wind coefficients file

1.00 0.924 0.707 0.383 0.00 -.383 -.707 -.924 -1.0 -.924 -.707

0.00 0.383 0.707 0.924 1.00 0.924 0.707 0.383 0.00 -.383 -.707

-1.0 -.924 -.707 -.383 0.00 0.383 0.707 0.924 1.00 0.924 0.707

0.00 -.383 -.707 -.924 -1.0 -.924 -.707 -.383 0.00 0.383 0.707

north

east

south

west

1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.

-1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1". -1.

plus-one

minus-one

The four lines of data shown with ". .." are actually complete on the

file; they are truncated here because of margin requirements. The four

profiles given direction names have coefficients equal to the cosine of angle

between the wind direction and the direction indicated.

B.4.1 Wind Pressure Test #1

The wind pressure tests are based on a simple three-element airflow

network:

Pl P2 P3 P4

Powerlaw coefficients of 0.000848528, 0.00848528, and 0.000848528 are

assumed for the three elements. These combine to give an effective powerlaw

coefficient of 0.000598506 .

Assuming standard barometric pressure, an ambient temperature of 20°C,

and a wind speed of 5.0 m/s gives a potential wind pressure of (~pV2) 15.0519

Fa. Test #1 uses wind coefficient profiles that are independent of wind

direction. Both links use a wind pressure modifier of 2. These combine to

give a total wind pressure across the network of 60.2075 Fa. Therefore, the

expected airflow through the network is 0.00509605 kg/so

53
wind pressure test #1 input file

node

node

node

node

node-l

node-2

node-3

node-4

0.0

0.0

0.0

0.0

c 0.0

v

v

c 0.0

orf-O.OOl

orf-O.Ol0

element

element

plr 2.276ge-7 2.276ge-7

plr 7.2e-6 7.2e-6

0.000848528

0.00848528

0.5

0.5

orf -0.001 rnA2

orf -0.01 rnA2

link

link

link

link-l

link-2

link-3

node-l

node-2

node-4

0.0

0.0

0.0

node-2

node-3

node-3

0.0

0.0

0.0

orf-O.OOl

orf-O.Ol0

orf-O.OOl

plus-one 2JO

null

minus-one 2.10

*********

Note how the order of the nodes in link-l and link-3 is from the "outside" to

the "inside" which requires a reversal in the numbering.

B.4.2 Wind Pressure Test #2

This test is almost the same as the first test except the "north" and

"south" wind coefficient profiles are used, instead of "plus-one" and "minusone". Since these profiles vary as cosD, where D is the wind direction, the

total wind pressure across the network of 60.207SocosD Pa, and the expected

airflow through the network is IO.OOSO960SojlcosDI kg/so

wind

directions

0, 360

30, 330

60, 300

90, 270

120, 240

150, 210

180

expected

flow rate

0.00509605 kg/s

0.00474241

0.00360345

0.0

-0.00360345

-0.00474241

-0.00509605

computed

flow rate

0.00509604 kg/s

0.00470292

0.00357087

0.0

-0.00357087

-0.00470292

-0.00509604

The difference between the computed and the expected flow rates is apparently

caused by the linear interpolation of the wind coefficient profiles.

54

20.0

20.0

20.0

20.0
B.5 Duct Element Tests

This is a test of the AIRNET Darcy-Weisbach-Colebrook duct model. The

expected flow rates are based on the ASHRAE friction chart (ASHRAE, 1985, pp

33.5 & 33.26). It involves a simple network in which three serial duct

elements are in parallel with a duct element of equivalent length. The flows

in both sides of the network should be equal.

~e=C2 

=:=C3=== 11

P2 P3 IreP4

\

-vl

P1 ..

Two design points were selected from the chart. The first point is for a duct

250 rom in diameter (area = 0.04909 m2) with a pressure drop of 0.9 Pa/m of

duct, giving a flow velocity of 4.0 m/s and a flow rate of 195 l/s. The

second point is for a duct 630 rom in diameter (area = 0.31172 m2) with a

pressure drop of 4.0 Pa/m of duct, giving a flow velocity of 16.0 m/s and a

flow rate of 5000 l/s. The absolute roughness dimension of both ducts is 0.15

rom. Duct lengths are 2 m, 3 m, 5 m, and 10 m for elements 1 through 4,

respectively.

In the first case the expected flow rate is 0.24 kg/s while in the second

case it is about 6.0 kg/so The values computed by AIRNET are 0.245 and 6.19

kg/s, respectively. These values are within 3% of those on the ASHRAE chart.

duct (Darcy-Weisbach-Colebrook model) test #1 input file

20.0

20.0

20.0

20.0

node 9.0

node

node

node

node-l

node-2

node-3

node-4

0.0

0.0

0.0

0.0

c

v

v

c 0.0

0.00015

128.

0.00015

128.

0.00015

128.

0.00015

128.

duct-2x25 dwc 2.0 duct 2m x 25cm

0.0

duct-3x25 dwc 3.0

0.0

duct-5x25 dwc 5.0

0.0

duct-10x25 dwc 10.0

0.0

0.25

64.0

0.25

64.0

0.25

64.0

0.25

64.0

0.04909

0.0

0.04909

0.0

0.04909

0.0

0.04909

0.0

element

element duct 3m x 25cm

element duct Sm x 2Scm

element duct 10m x 25crn

null

null

null

null

node-2

node-3

node-4

node-4

0.0

0.0

0.0

0.0

duct-2x25

duct-3x25

duct-5x25

duct-1Ox25

0.0

0.0

0.0

0.0

link

link

link

link

link-l

link-2

link-3

link-4

node-l

node-2

node-3

node-l

*********

55
B.6 Doorway Element Test

This test compares two different ways to simulate the two-way airflows

through doorways. The test represents a 0.8 m by 2.0 m doorway connecting two

rooms with a 4°C temperature difference as shown below. The computed results

can be compared to the correlation presented by Weber and Kearney (1980) for

heat transfer through a doorway:

Nu = aoProJGr (B. 6.1) where

Nu = Nusselt number = hoH/k,

Pr = Prandtl number = c o~/k,

Gr = Grashof number = p~ogopoH3oAT/~2, and

a = experimentally determined coefficient.

In addition,

p = density,

~ = viscosity,

Cp = specific heat,

k = thermal conductivity,

g = acceleration of gravity,

P = coefficient of thermal expansion = -AP/(poAT) ,

H = doorway height,

h = convection coefficient = q/(WoHoAT) ,

q = heat flux rate = wocpoAT,

W = doorway width, and

w = mass flow rate.

Substituting into equation (B.6.!) and simplifying gives:

w = aoWoJpogoApoH3 (B.6.2)

Assuming the doorway behaves as an orifice leads to the relationship Q =

c/3, where c is the orifice coefficient. Substituting the values Q = 0.26, p

= 1.20415 kg/m3, ~p = 0.0164312 kg/m3, g = 9.80 m/s2, H = 2.0 m, and W = 0.8 m

gives w = 0.25906 kg/so

The doorway element test compares the doorway element model with an

approximation composed of multiple orifice elements.

Tl=18°C T2=22°C T1=18°C T2=22°C

-z=l-

-0.5-

1.

P1 = 0

1. .2

Pi = 0

.2 -z==o-

--.5-

~

-z = -1 -

Doorway simulated

by multiple

powerlawelements

Doorway simulated

by doorway element

56
In the first case a doorway element with a discharge coefficient of 0.78

is used giving computed two-way airflows of 0.25913 and 0.25899 kg/so These

are essentially exact. In the second case ten 0.16 m2 powerlaw openings at

equally spaced heights are used to represent a doorway. The computed two-way

airflow is 0.261 kg/s which is within 1% of the exact value. Both cases were

solved in two iterations.

doorway test #1 input file

node

node

node-l

node-2

0.0

0.0

18.0

22.0

c 0.0

v

orf-O.16

dor-l.60

element

element

plr 0.0015575 0.0015575 0.176494 .500

dor 0.015575 0.015575 1.76494 .500

0.0001 2.0 0.8 0.78

opening-.16m"2

doorway-l.6m"2

link

link

link

link

link

link

link

link

link

link

link

link-l

link-2

link-3

link-4

link-S

link-G

link-7

link-8

link-9

link-lO

link-ll

node-l

node-l

node-l

node-l

node-l

node-l

node-l

node-l

node-l

node-l

node-l

o.

o.

o.

o.

o.

-0.

-0.

-0.

-0.-0.

-1.

node-2

node-2

node-2

node-2

node-2

node-2

node-2

node-2

node-2

node-2

node-2

0.9

0.7

0.5

0.3

0.1

-0.1

-0.3

-0.5

-0.7

-0.9

-1.0

orf-O.16

orf-O.16

orf-O.16

orf-O.16

orf-O.16

orf-O.16

orf-O.16

orf-O.16

orf-O.16

orf-O.16

dor-l.60

null

null

null

null

null

null

null

null

null

null

null

*********

Note the heights used with link-!! for the doorway element. In this case the

reference heights of the two nodes is one meter above the floor, so the

doorway element must be located at -!.O relative to each node. !

57

9

7

5

3

1

1

3

5

7

9

0
Constant Flow Element Tests

Constant Flow Test #1

In this test a powerlaw airflow elements and a constant flow element are

arranged in series:

.Cl

Pl P3

The flow rate is set to 1.0 kg/so Since the flow through both elements must

be the same and the powerlaw coefficient is 0.0848528, the expected pressure

drop across the powerlaw element is 115.342 Fa.

constant flow test #1 input file

node

node

node

node-l

node-2

node-3

0.0

0.0

0.0

20.0

20.0

20.0

c 0.0

v

c 0.0

element

element

orf-0.1000 plr 2.276ge-4 2.276ge-4

flow-1.000 cfr 1.0

0.0848528 0.5 orf -0.1 m~2

constant flow of 1.000 kg/s

link

link

link-l

link-2

node-l

node-2

0.0

0.0

node-2

node-3

0.0

0.0

orf-0.1000

flow-l .000

null

null

*********

B.7.2 Constant Flow Test #2

This test uses the same airflow network as powerlaw test #3 except that a

constant flow element is added at node 13. Node 12 becomes an unknown

pressure and new node 13 must be a known pressure. If the constant flow is

set to 0.0611025 kgjs and the pressure of node 1 set to 50 Fa, this test

should give the same pressures and flows as powerlaw test #3.

It does give the same answers (to within the convergence limits) and does

this in only 6 iterations. See the input file AFDATA.CF2 on the diskette.

58
B.8 Fan Element Tests

The fan element tests are based on problems in the textbook by Osborne

(1977, pp 71-75). All tests use a fan whose performance curve was developed

from the following characteristics:

volume flow Total pressure

0.0 m3/s 750 Pa

1.0 755

2.0 730

3.0 590

4.0 275

The performance curve (assuming a density of 1.204 kg/m3) is:

w ; 764.429 -18.2922.P + 19.4633.p2 -7.63940.p3

This curve does not include the contraflexure evident in the table of

characteristics. This was done by increasing the pressure rise at zero flow

until the contraflexure was eliminated (765 Pa).

B.B.! Fan Element Test #1

This corresponds to problem 1 in Osborne. It consists of a fan in series

with a powerlaw element. The expected answer is 3.167 kg/s (2.63 m3/s) at a

total fan pressure of 660 Pa. The computed results are 3.158 kg/s at 660.3

Pa.

B.8.2 Fan Element Test #2

This corresponds to problem 5 in Osborne. It consists of two fans in

parallel and three powerlaw elements forming a complete circuit.

-I: 1 =.=C1 ~-.P ~L.r-:; 

1 ~F 2 =.~C2 ~ P 4

I P3

computed

291.1 Pa

2.909 kg/s

3.616 kg/s

6.524 kg/s

values

P4 -Pi

Wi

Wz

W3

expected

293. Pa

2 .91 kg/ s (2.42 rn3 / s)

3 .61 kg/ s (3.00 rn3 / s)

6.53 kg/s (5.42 rn3/s)

The differences can be attributed to the inaccuracies in the fit of the

fan performance coefficients and the inaccuracies in the graphic procedures

used to compute the textbook answers.

59
B.B.! Fan Element Test #3

This corresponds to problem 6 in Osborne. It consists of two fans in

series and three powerlaw elements forming an open network.

values

P3 -P5

Wl

Wz

W3

computed

72.5 Pa

4.171 kg/s

3.713 kg/s

0.459 kg/s

expected

70. Pa

4.19 kgjs (3.48 m3js)

3.73 kgjs (3.10 m3js)

0.46 kgjs (0.38 m3js)

fan test #3 input file: comparison to Osborne problem #6

node-l

node-2

node-3

node-4

node-S

0.0

0.0

0.0

0.0

0.0

no

no

no

no

no

c

v

v

v

c

7.2e-6 0.084853

0.10 3 -100.0

19.4633 -7.63940

19.4633 -7.63940

19.4633 -7.63940

1.S8Se-OS 0.19005

1.110e-OS 0.13306

4.092e-06 0.04907

fan 3.0e-5 .500 Osborne fan A

750.0 5.00

-18.2922

-18.2922

-18.2922

plr 1.0e-5

plr 1.0e-5

plr 1.0e-5

fan-A

1.204

764.429

764.429

764.429

rest1

rest2

rest3

element

1.05.0100 no contraflexure

0.5

0.5

0.5

resistance 1

resistance 2

resistance 3

element

element

element

node-2

node-3

node-4

node-S

node-S

0.0

0.0

0.0

0.0

0.0

null

null

null

null

null

link-l

link-2

link-3

link-4

link-S

node-l

node-2

node-3

node-3

node-4

0.0

0.0

0.0

0.0

0.0

fan-A

rest!

fan-A

rest3

rest2

link

link

link

link

link

*********

60

20.0

20.0

20.0

20.0

20.0

dede

de

de

de
B.9 Quadratic Flow Element Test

The test of the quadratic flow element is an extension of the duct

element test. Is also shows the possibilities of modeling ducts components by

powerlaw or quadratic flow elements. The airflow network is:

Element Do is a 10 meter long duct element; elements Dl' Dz' and D3 are 2, 3,

and 5 meter long duct elements. Elements Cl' Cz' and C3 are powerlaw elements

whose coefficients were generated by the ELEMENT program to correspond to the

duct elements. The elements were fitted at Reynolds numbers of 10000 and

40000. Elements Ql' Qz, and Q3 are quadratic flow elements with coefficients

generated to correspond to the duct elements. There should be equal flows

through each branch of the network. The results of the test are summarized

below:

AP

0.01

0.02

0.04

0.08

0.16

0.32

0.64

1.28

2.56

5.12

10.24

20.48

40.96

Fa

0.0049276

0.0073755

0.0109604

0.0161748

0.0237102

0.0345369

0.0500175

0.0720693

0.1033972

0.1478221

0.2107488

0.2998266

0.4258877

Fd

0.0049276

0.0073755

0.0109604

0.0161748

0.0237108

0.0345369

0.0500174

0.0720691

0.1033964

0.1478200

0.2107428

0.2998018

0.4258392

Fp diff

0.0045302 -8.06%

0.0080704 +9.39%

0.0116082 +5.91%

0.0166968 +3.23%

0.0240167 +1.29%

0.0345440 +0.02%

0.0496870 -0.66%

0.0714681 -0.83%

0.1027976 -0.58%

0.1478615 +0.03%

0.2126813 +0.92%

0.3059127 +2.03%

0.4400455 +3.32%

Fq diff Re

0.0040626 -17.55% 1311.

0.0066006 -10.53% 2G53.

0.0103411 -5.65% 3050.

0.0157490 -2.63% 4502.

0.0234849 -0.95% 6597.

0.0344891 -0.14% 9612.

0.0500971 +0.16% 13920.

0.0722028 +0.18% 20057.

0.1034881 +0.09% 28776.

0.1477486 -0.05% 41139.

0.2103540 -0.18% 58652.

0.2988997 -0.31% 83443.

0.4241279 -0.41% 118526.

where

p.P = Pi -PSI

Fo = mass flow (kg/s) through element Do,

Fd = mass flow through the 3 duct elements branch,

F = mass flow through the powerlaw elements branch,

FP = mass flow through the quadratic flow elements branch, q .

diff = percent d1fference compared to Fol and

Re = Reynolds number computed from Fo.

The agreement for both fitted elements is quite good in the turbulent region,

Re > 4000, with the quadratic flow element giving the best agreement. Good

agreement is not expected in the transition and laminar region, Re < 4000.

61
quadratic flow relationship test #1 input file

node

node

node

node

node

node

node

node

o.

o.

o.

o.

o.

o.

o.

o.

20.0

20.0

20.0

20.0

20.0

20.0

20.0

20.0

source

duct-l

duct-2

dplr-l

dplr-2

dqfr-l

dqfr-2

sink

c 0.32

v

v

v

v

v

v

c 0.0

element duct-2x25 dwc 2 0.25 0.0490874 0.0015 duct -2m long by 25cm dia

0 64 0 64

element duct-3x25 dwc 3 0.25 0.0490874 0.0015 duct -3m long by 25cd dia

0 64 0 64

element duct-5x25 dwc 5 0.25 0.0490874 0.0015 duct -5m long by 25cm dia

0 64 0 64

element duct-10x25 dwc 10 0.25 0.0490874 0.0015 duct -10m long by 25cm dia

0 64 0 64

element dplr-2x25 plr 3.40642e-05 3.40642e-05 0.133079 0.524431 plr-duct ..

element dplr-3x25 plr 2.27095e-05 2.27095e-05 0.107587 0.524431 plr-duct ..

element dplr-5x25 plr 1.36257e-05 1.36257e-05 0.0823032 0.524431 plr-duct .

element dqfr-2x25 qfr 0.310242 44.8089 qfr-duct -2m long by 25cm dia

element dqfr-3x25 qfr 0.465364 67.2134 qfr-duct -3m long by 25cm dia

element dqfr-5x25 qfr 0.775606 112.022 qfr-duct -5m long by 25cm dia

link

link

link

link

link

link

link

link

link

link

duct-Q

duct-l

duct-2

duct-3

dplr-l

dplr-2

dplr-3

dqfr-l

dqfr-2

dqfr-3

o.

o.

o.

o.

o.

o.

o.

o.

o.

o.

sink

duct-l

duct-2

sink

dplr-l

dplr-2

sink

dqfr-l

dqfr-2

sink

o.

o.

o.

o.

o.

o.

o.

o.

o.

o.

duct-1Ox25

duct-2x25

duct-3x25

duct-5x25

dplr-2x25

dplr-3x25

dplr-5x25

dqfr-2x25

dqfr-3x25

dqfr-5x25

null

null

null

null

null

null

null

null

null

null

source

source

duct-l

duct-2

source

dplr-l

dplr-2

source

dqfr-l

dqfr-2

*********

Note:

The a and b values for the quadratic flow elements are proportional to the

length of each element. This could be convenient for modeling duct networks.

62

,0

,0

,0,0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0
B.1O Execution Time Tests

These tests are designed to provide simple benchmarks for relatively

large airflow networks. Analytic solutions have not been developed for these

cases. The runs were made on a PC compatible computer (4.77 MHz 8088 CPU with

8087 math coprocessor). This is a minimal microcomputer by today's standards.

The tests are based on the 9-node floor module mentioned in the main

report and shown in figure 10. Test faces were developed for buildings of up

to 32 floors. Versions of AIRNET using both the small and large memory

modules are compared in the following table.

Small memory model

floors nodes

4 37

8 73

16 145

Large memory model

16 145

32 289

calculation time

input airflows

2.8 s 1.6 s

5.3 3.3

9.8 6.7

iterations

5

5

5

unallocated memory

z 44000 bytes

31000

7000

5

5

12.4

24.7

7.6

15.3

333000

265000

Four points to note are that (1) calculation time increases almost

linearly with the number of nodes because of the sparse equation solver,

(2) the input files process slowly because of their length, (3) the large

memory model requires slightly more time than the small model, and (4) this C

version of AIRNET is about 40% faster than the FORTRAN version used in the

comparisons in the main report.

These network data files can be easily modified to show the effect of

improper node order on sparsity and performance by altering the sequence of

nodes on each floor. In the 37 node case the number of nonzero elements in

[A] increases from 127 to 253 (an increase of 1008 bytes) and execution time

increases from 1.6 to 2.6 seconds. The change in the fill pattern of the

Jacobian is shown in the two figures below. !

.

.

.

..

..

...

.+.

~.

.+. ++. .++ .++

.++

.++

..++

..++

....+.

.~

.+.

++.

.++

.++

.++

.++

..++

..++

...

.+.

.~

.+.

++.

.~

.+. e+

e

.

1234567890123456789012345678901234567 1234567890123456789012345678901234567

63
B.ll References

ASHRAE, 1985. ASHRAE Handbook -1985 Fundamentals, ASHRAE, Atlanta GA.

Judkoff, R., Wortman, D., O'Doherty, R. and Burch, J., 1983. "A Methodology

for Validating Building Energy Analysis Simulations", Solar Energy Research

Institute report SERI/TR-254-1508, Golden, CO.

Klote, J.H. & Fothergill, J.W.Jr., 1983.

Buildings", ASHRAE, Atlanta GA.

"Design of Smoke Control Systems for

Osborne, W.C., 1977 ~, 2nd edition, Pergamon Press.

Weber, D.D., and Kearney, R.J., 1980. "Natural Convective Heat Transfer

Through an Aperture in Passive Solar Heated Buildings", 5th National Passive

Solar Conference, American Solar Energy Society, pp. 1037-1041.

64
APPENDIX C: Calculational Details

c.! Newton's Method

This section reviews the theory of Newton's method for simultaneous

nonlinear equations (Stoecker, 1971) as applied to the airflow equations.

Assume a case consisting of three nodes which are to be solved for three

unknown pressures Pi, Pz, and P3' Conservation of mass at each node is used

to set up three nonlinear equations in terms of the unknown pressures:

fi (Pi' P2' P3) = 0

(C.l.la)

f2 (Pl' P2' Pa) = 0

(C.l.lb)

f3 (Pl' P2' P3) = 0

(C .1.lc

If the correct solution is Plc, P2c' and P3c, a Taylor series expansion

retaining only the low-order terms gives:

f1(Pl' 

PZ' P3) z

(C.l.2)

8f1(Plc. 

P2c. P3C)

f1(P1C' P2c. P3c) +

(Pl -Plc)

8Pl

8f1(Plc, P2c' P3c) 8fl (Plc, P2c, P3c)

+ (Pz -Pzc + (P3i -P3c

8P2 aP3

and similarly for f2 and f3, where the variables Pl, P2, and P3 are only an

approximation to the correct solution.

bi ~ aii(Pi -Pic) + a12(P2 -P2c) + a13(P3 -P3c) (C.l.3a)

b2 ~ a21(Pl -P1c + azz (Pz -Pz c) + aZ3 (P3 -P3 c) (C.l.3b)

b3 ~ a31(Pl -Plc + a32 (P2 -P2 c) + a33 (P3 -P3 c) (C.l.3c)

where

bi = fi (Pl' P2, P3 (C.l.3d) and

8fi (Pic P2c. P3c aij =

(C.l.3e)

8P.J

However, since the exact pressure values are not yet known, compute the aij

terms at the approximate pressure values:

aft (Pl' P2, P3

(C.l.3f) aij ~

8Pj

Equations (C.l.3) are a set of simultaneous linear algebraic equations in

terms of the unknowns: (Pl -Plc), (Pz -Pzc)' and (P3 -P3c). Solution of

65
the equations for these values is used to improve the approximate pressurevalues:

Pl* = Pl -(Pl -P1c) (C.l.4a)

P2* = P2 -(P2 -P2c) (C.l.4b)

P3* = P3 -(P3 -P3c) (C.l.4c)

The new pressure values replace the old pressure values in equations (C.l.3)

and solution continues iteratively until a sufficiently accurate set of

pressure values has been computed.

These relationships are expressed in vector/matrix form in the main

paper: Each airflow element, i, relates the mass flow rate, Wi' through the

element to the pressure drop, APi, across it. A new estimate of the vector of

all node pressures, {P}*, is computed from the current estimate of pressures,

{P}, by

(P)* = (P) -(C) (C.l.4')

where the correction vector, {C}, is computed by the matrix relationship

:J] {C} = {B} (C.l.3')

{B} is a column vector with each element given by

Bn = ~ Wi

i

(C.l.3d')

where n is the node number and i indicates all flow paths connecting node n to

other nodes, and [J] is the square (i.e. N by N for a network of N nodes)

Jacobian matrix whose elements are given by

In.m = I --

a;i (C.l.3f')

.mJ.

Since the airflows for the various airflow elements are computed from

relationships of the form w = f(~P), where

~P = Pn -Pm + PS + PW,

Pn, Pm = total pressures at nodes nand m,

PS = pressure difference due to density and height differences, and

PW = pressure difference due to wind,

the partial derivative needed for [J] in equation (C.l.3f') are related by

8w/8Pm = -8w/8Pn (C.l.5)

These last two equations indicate that the Jacobian has a strong diagonal and

is symmetric, two properties which can be used to speed the factoring of theJacobian.

66
Solvin!!: Simultaneous Linear Eauations

The following discussion is designed to provide a brief introduction to

the method used to solve simultaneous equations in AIRNET. See Pissanetzky

(1984) and Dhatt & Touzot (1984) for details.

The general form of the problem is

[A] (X) = (B) (C.2.!) where

[A] is a square matrix of coefficients,

{X} is a vector of solution values,

{B} is a matrix of coefficients.

It is often computationally preferable to factor the [A]matrix into triangular

matrices which are simple to solve nUmerically.

Original problem: [A] {X} = {B}

Factored problem: [L] [U] {X} = {B} (C.2.2)

Forward substitution [L] (B') = (B) (C.2.3)

Back substitution: [U] {X} = {B'} (C.2.4)

The processing of the matrices is dependent on the way the matrices are

represented. In this case we are interested in the skyline storage method

[L] [U] = [A]

1

L21

L31

...

0

1

L32

...

0

0

1

..

Dl

0

0

...

UZ1

Dz

0

...

U31

U32

D3

...

E1

R21

R31

...

CZ1

Ez

R3Z

...

C31

C32

E3

...

=

n = 1:

n = 2:

1.Dl = El implies

L2 1 E1 = R2 1

1.U21 = C21

L21U21 + D2 = E2

L31El = R31

L3 1 U2 1 + L3 2 E2 = R32

l'U31 = C31

L21U31 + U32 = C32

L31U31 + ~2U32 + D3 = E3etc.

3:

D1 = El

Lzl = R21/Dl

U21 = C21

D2 = E2 -LzlU21

L31 = R31/Dl

L32 = (R32 -L31U21)/D2

U31 = C31

U32= C32 -L21U31

D3 = E3 -L31U31 -L32U32

n = 4:

67
Factoring algorithm:

for all values of n from 1 to N:

~i = (~i i-1

-I ~mUim) / Di

m=1

i = 1, 2, .

, n-l

i = 1, 2, ..., n-1

Dn = ~

n-l

-I ~mUnm

m=l

Forward substitution algorithm:

for all values of n from 1 to N:

Back substitution algorithm:

for all values of n from N to 1:

Bn = Bn /Dn

Bi = Bi -UniBn i = 1, 2, ..., n-1

(~ = Bn when finished)

The development of the skyline method is based on the simple observation that

array locations above the topmost nonzero element in any column above the

diagonal and to the left of the first nonzero element in any row below the

diagonal must remain zero during factoring. It is unnecessary to allocate

array space for those locations that are always zero. The [A] matrix is

reorganized into the main diagonal vector of N elements, the lower triangular

matrix of N variable length rows, and the upper triangular matrix of N

variable length columns. The factoring and solution algorithms are then

reorganized to accommodate this new data structure with its new numbering of

array locations. All operations can be performed in place, i.e., without

creating separate arrays for setting up the problem and for the factoring.

68

i-l

Uni = Cni -I LimUnm

m=l

n-l

Bn = Bn -L ~j Bj

j=l
Arrays in C

2-dimensional arrays

2-D arrays can be handled directly or by arrays of pointers. Reviewing: 

declarations (Kernighan & Ritchie, pp 103-110) for a 2-D array:

double

colwnns.

aa [3] [4] ; is a declaration for an array of 3 rows and 4

aa[O] [.0]

aa [0] [1]

aa[O] [2]

aa[0][3]

aa [1] [0]

aa [ 1] [1]

aa[l] [2]

aa[l] [3]

aa[2] [0]

aa [2] [1]

aa[2] [2]

aa[2] [3]

aa and aa[O] --> where aa[i][j] is both a location in

memory and the C access syntax for it

"-->" means "points to".

aa[l -->

aa[2] -->

Note that the array indices start at zero for both columns and rows.

When aa is passed to a function the number of columns must be given in

the declaration in the function:

double aa [ ] [4] ;

It is not necessary to state the number of rows (the first dimension in a

multidimensional array).

Now consider the use of an array of pointers:

double

*pa[4];

is a declaration for an array of 4 pointers

aa(O,O)

aa(O,l)

aa(O,2)

aa(O,3)

aa(l,O)

aa(l,l)

aa(1,2)

aa(1,3)

aa(2,O)

aa(2,1)

aa(2,2)

aa(2,3)

pa[O] --> where aa(i,j) is a location in memory

for the appropriate value.

par! -->

pa[2] -->

pa[3] -->

Element aa(O,O) is accessed by *(pa+O+O) (= *pa) or pa[O][O]. Element

aa(1,3) is accessed by *(pa[l]+3) or pa[l][3]. The pointer past the end of

the array is used to determine the end of the last row during dumps. The

pointer vector must be dimensioned one more than the number of rows in the

array.

69
Since there is a problem in writing functions which are passed 2-D arrays

and since it is desirable to always use the same method of representing such

arrays, we will always represent a 2-D array by means of an array of pointers.

For arrays with very few columns, the addition of an array of pointers

represents a significant increase in storage. In many such cases it may be

more desirable to use a l-D array of structures, where the structure element

names can help document the program.

We often want to use arrays with elements 1 through N instead of 0

through N-l. In this case, point to the location just before each row:

double *pa[S]; is a declaration for an array of 5 pointers

pa[l] -->

para] -->

---

aa(l,l)

aa(1,2)

aa(1,3)

aa(1,4)

aa(2,1)

aa(2,2)

aa(2,3)

aa(2,4)

aa(3,1)

aa(3,2)

aa(3,3)

aa(3,4)

->

->

pa[4] -->

The zero pointer is not used, and the final pointer is needed to

determine the end of the last row during dumps. This means that the pointer

array will have to be dimensioned two more than ,the number of rows (if we want

row N accessed by pa[N] instead of pa[N-l]). Element aa(1,3) is accessed by

*(pa[l]+3) or pa[l] [3]. ,

The user must know whether the array begins at zero or at one. :

70
C.3.2 Variable Length Rows

We may also want to use arrays with variable length rows.

the use of an array of pointers is essential: In this case

double

*pa[S];

is a declaration for an array of 5 pointers.

aa(O,O)

aa(O,l)

aa(O,2)

aa(l,O)

aa(l,l)

aa(2,O)

aa(2,1)

aa(3,O)

aa(3,1)

aa(3,2)

aa(3,3)

pa[O] 

-->

pa[l -->

pa[2] -->

pa[3] -->

pa[4 -->

Element aa(O,O) is accessed by *(pa+O+O) (= *pa) or pa[O][O]. Element

aa(2,1) is accessed by *(pa[2]+1) or pa[2][1]. The pointer to the end of the

array is used to determine the end of the last row during dumps. The pointer

vector must be dimensioned one more than the number of rows in the array.

In the case of arrays with elements 1 through N instead of 0 through N-l,

point to the location just before each row:

double

*pa[6];

is a declaration for an array of 6 pointers

pari] -->

para] -->

---

aa(l,l)

aa(1,2)

aa(1,3)

aa(2,1)

aa(2,2)

aa(3,1)

aa(3,2)

aa(4,1)

aa(4,2)

aa(4,3)

aa(4,4)

pa[2] -->

pa[3] -->

pa[4] -->

parS] -->

The zero pointer is not used, and the final pointer is needed to

determine the end of the last row during dumps. This means that the pointer

array will have to be dimensioned two more than the number of rows. Element

aa(1,3) is accessed by *(pa[1]+3) or pa[1][3].

The user must know whether the array begins at zero or at one.

71
C.4 Stack Effect

The following figure shows two rooms represented by airflow nodes nand

m, respectively. It is assumed that each node can be characterized by a

single temperature and a single static pressure at some height relative to a

common data plane. The two rooms are shown with an airflow element connecting

them. The inlet and outlet to the element are at different heights from each

other and from the nodes representing the rooms to show the entire calculation

for height differences.

node n-.

IZn node m

.-

Zm

Zl

reference height

Analysis of airflow through the element i is based on Bernoulli's

equation and its assumptions.

A.Pi = (P1 + pV12/2) -(P2 + pV22/2) + pg(Zl -Z2) (C.4.1) where

APi = sum of all friction and dynamics losses (Pa),

Pl' Pz = entry and exit static pressures (Pa),Vl' 

Vz = entry and exit velocities (m/s) ,

p = density of fluid flowing through the element (kg/m3

g = acceleration of gravity (9.81 m/sz), and

Zl' Zz = entry and exit elevations (m).

P may be either Pn or Pm depending on the direction of flow. One possibility

is to let P equal the average of the two densities and accept the inaccuracy

in the computed airflow. The other possibility is to select the pressure

based on the most recently computed flow direction and let the iterative

solution of the airflow equation operate.

Equation (C.4.1) defines a sign convention for the direction of airflow:

positive flow is from point 1 to point 2 --node n to node m. Equation

(C.4.1) can be simplified for use in the airflow computer algorithm by

defining several related terms. Dynamic pressures are the pv2/2 terms in

equation (C.4.1) and total pressure is defined to be P = p + pV2/2. If nodes

nand m represent rooms, the dynamic pressures are effectively zero. If the

nodes are part of a duct network, there will be a positive dynamic pressure.

The pressures at the inlet and outlet of the airflow element can be

related to the node pressures by the hydrostatic law:

Pl = Pn + png(Zn -Zl) = Pn -Pngh1 where (C.4.2a) h1 = Zl-Zn"

and

P2 = Pm + pmg(Zm -Z2) = Pn -Pmgh2 where (C.4.2b) hz = ZZ-Zmo

72
The relative heights, h1 and hz, are a convenient way of expressing the

element inlet and outlet heights. Similar flow elements in similar rooms tend

to have the same relative heights even though the rooms are at differentheights. 

If the element is part of a duct network, the relative heights will

be zero. Equation (C.4.1) is thus reduced to

APi = P n -Pm + P g (Zn +h1 -Zm -hZ) -Pn gh1 + Pm ghz

(C.4.3)

The terms [pg(Zn+hl-zm-hz) -Pngh1 + Pmghz

the stack pressure, PSi' acting on element i.

direction:

] 

can be collectively calledFor 

airflow in the positive

PSi = png(Zn-Zm) + hzg(Pm-Pn)

(C.4.4a)

For flow in the negative direction:

PSi = Pmg(Zn-Zm) + h1g(Pm-Pn

(C.4.4b)

If the direction of flow is unknown, use

PSi = (g/2)[(pn+pm)(zn-zm) + (Pm-Pn)(h1+hZ)]

(C.4.4c)

Equation (C.4.4c) could also be used if it is assumed that the air temperature

within the element is the average of the two nodal temperatures.

Elements in building envelope surfaces connect the rooms to ambient air

which may be moving at significant speed. This speed is reduced to zero as

the air meets the building surface; the dynamic pressure of the ambient air is

converted to an additional pressure on the element, PWi, which will be called

the wind pressure. Wind pressure may act in either a positive or negative

fashion on the element airflow. Equation (C.4.3) is thus reduced to

APi = Pn -Pm + PSi + PWi (C.4.5)

Calculation of the stack pressure interacts with the iterative solution

for the nodal pressures. It is possible to recompute PSi at every iteration,

but this can mean that as the estimated pressures change with each iteration

flow directions may also change leading to different values for PSi during the

iterations. These changing estimates of the stack pressure interfere with the

convergence of the nodal pressure solution. It appears to be better to

compute PSi for each component based on the most recent node pressures and to

leave these values unchanged during the solution process. The new pressure

values can be used to recompute PSi and then again solve for pressures in a

very few iterations because the pressures will probably change by only very

small amounts. Because of the relatively small effect of flow direction on

the overall solution, it is probably not necessary to recompute PSi during a

given timestep in a transient building performance simulation.

More complicated stack effects can be created by assuming the nodal

temperature varies with height. This could be used to model the effects of

various temperature stratifications. Three fairly simple temperature

variations come to mind. First, the temperature may vary linearly with

height. Second, each node could consist of two regions of constant

temperature, perhaps with the height of the interface varying. Third, the

first two models could be combined.

73
C.5 Fitting Dowerlaw coefficients

Data is seldom available for direct substitution into the powerlaw

equation:

w = c Jp (AP)X (c.s.!)

where

w = mass flow rate of air,

C = flow coefficient,

~p = pressure difference

p = air density, and

x = flow exponent.

It is usually necessary to convert the available data to that form.

If x is known or can be assumed (to be ~ for example), C can be computed

from the inverse of equation (C.S.l)

:jp (flP)X c = w /

(C.5.2)

For example, this can be used with the component leakage area formulation

which has been used to characterize openings for infiltrations calculations

(ASHRAE, 1985, p 22.16). The effective leakage area is given by

L = 10000 Q / JUP7P (C.5.3)

where

L = the effective leakage area (cm2),

Q = the volumetric flow (m3/s), \

AP = standard pressure difference (4 Fa), and

p = density of air (1.2 kg/m3).

This corresponds to a value of C = 0.0001414 L for use with equation (C.S.1).

When two values, wl(~Pl) and W2(~P2)' are known, both C and x can be

computed:

x= In(wl) -In(wz) ] / In(!!Pl -In(APz) (C.5.4)

with C then computed from equation (C.5.2). This can be used to convert

detailed duct calculations to powerlaw formulation. Tests indicate that the

characteristics of ducts with friction modeled by the Colebrook equation can

be replaced by the powerlaw representation with an accuracy of better than 2%

for flows that vary by up to a factor of four (a typical variation in a VAV

system). Equation (C.5.4) would also be useful for reducing a collection of

flow elements which are in series to a single flow element.

74
Constant Flow Element

One particularly simple but useful airflow element sets a constant flow

between two nodes. Since the flow is constant, the partial derivatives of

flow with respect to the node pressures must be zero. The constant flow

element element does not contribute to the Jacobian, [A], but it does add to

the right side vector, {B}.

Constant flow elements do not mathematically link the pressures of the

adjacent nodes. It is necessary that all node in the network be linked to

constant pressure nodes in order to have a unique solution. Violation of this

restriction will produce a division by zero somewhere in the solution of the

equations. Consider the following simple network:

.=c 1 =.=F 1 =.=C2 =.=F 2 =.=C3 =.

Pl P2 P3 P4 Ps P6

where P1 and P6 are known pressures, and Fl and F2 are known flows. Since the

flow through C2 is determined by P3-P4, and no other flow is related to those

pressures, there are not enough equations to determine P3 SD4 P4 uniquely. In

addition, it is possible for Fl and F2 to be assigned different flows, which

produces a physically impossible condition.

Constant Power Fan Element

One simple fan model is applicable for an assumption of near constant

operating conditions. This is the constant power fan model, which assumes

that a constant power is delivered by the fan to the airstream. The basic

equation is:

H = Q .A.P = w .A.p / p (-G.7.1

where

H == power (W),

Q == volume flow rate (m3/s) ,

bP == pressure rise (Pa),

w == mass flow rate (kg/s) , and

p == air density (kg/m3).

Equation (C.7.! rearranges to

w = H .p / (Pm -Pn) (C.7.2)

for which the partial derivatives are:

(C.l.3) .£.!i 8P = -Hp/bp2

n

~i 8P = Hp/f.p2 (C.7.4)

m

75
RelationshiI1s for Cracks

Baker, Sharples, and Ward (1987) indicate that infiltration openings can

be more accurately modeled by a quadratic relationship of the form

t,.P = A Q + B Q2 (C.B.l)

They give theoretical relationships between A and B and the physical

characteristics of the openings. These are

A = 12jjz/Ld3 (C.8.2a)

B = pC/2d2L2 (C.8.2b)

where

I" = viscosity,

P = density,

z = distance along the direction of flow,

d = crack width,

L = crack length, and

C = 1.5 + number of bends in the flow path.

Equations (C.8.2a,b) are easily converted to mass flow form for the quadratic

flow element:

a = 12jJz/pLd3 (C.8.3a)

and

b = C/2pd2L2 (C.8.3b)

Clarke (1985) also indicates a relationship for flow through cracks that

can be converted directly into a powerlaw airflow element. This relationship

is:

Q = k a (!J.P)X (C.8.4a)

x = 0.5 + 0.5 exp(-Wj2) (C.8.4b)

k = O.OO97(O.OO92)X (C.8.4c)

where

W = crack width (rom), and

a = crack length (m).

The powerlaw equation is

Wi = C JPn (6P)X (C.B.S)

Therefore, x is given by (C.8.4b) and

c = JPn a O.OO97(O.OO92)X (C.B.G)

76
C.9 References

ASHRAE (1985) ASHRAE Handbook -1985 Fundamentals, ASHRAE, Atlanta GA

Baker, P.R., Sharples, S., & Ward, I.C. (1987) "Air Flow Through Cracks"

Building and Environment, Vol 22, No 4, Pergamon, pp 293-304.

Clarke, J.A., 1985.

Bristol and Boston.

Energy Simulation in Building Design. Adam Hilger Ltd.

Conte, S.D. & de Boor, C., 1972. Elementary Numerical Analysis: an

Algorithmic Auuroach, McGraw-Hill, pp 84-87.

Dhatt, G., & Touzot, G. (translated by G. Cantin), 1984.

Method Displayed, John Wiley & Sons.

The Finite Element

Kernighan, B.W. & Ritchie, D.M., 1978.Hall, 

Inc.

The C PrograrnrninQ: Language, PrenticePissanetzky, S., 1984. Sparse Matrix Technologx. Academic Press.

Stoecker, W.F., 1971. Design of Thermal Systems, McGraw-Hill, pp 92,93.

77