NBSIR 87-3661
Indoor Air Quality Modeling
Phase II Report
James Axley
U.S. DEPARTMENT OF COMMERCE
National Bureau of Standards
National Engineering Laboratory
Center for Building Technology
Building Environment Division
Gaithersburg, MD 20899
October 1987
Prepared for:
U.S. Environmental Protection Agency
U.S. Department of Energy
-------
NBSIR 87-3661
INDOOR AIR QUALITY MODELING
PHASE II REPORT
James Axley
U.S. DEPARTMENT OF COMMERCE
National Bureau of Standards
National Engineering Laboratory
Center for Building Technology
Building Environment Division
Gaithersburg, MD 20899
October 1987
Prepared for:
U.S. Environmental Protection Agency
U.S. Department of Energy
U.S. DEPARTMENT OF COMMERCE, C. William Verity, Acting Secretary
NATIONAL BUREAU OF STANDARDS, Ernest Ambler, Director
-------
NBS: Indoor Air Quality Model ABSTRACT
Phase II Report
ABSTRACT
This interim report presents the results of Phase II of the NBS General Indoor Air
Pollution Concentration Model Project. It describes the theoretical basis of a
general-purpose nonreactive contaminant dispersal analysis model for
buildings, the computational implementation of a portion of this model in the
program CONTAM86, and examples of the application of this model to practical
problems of contaminant dispersal analysis. Presently the model is being
extended to handle problems of reactive contaminant dispersal analysis and full
computational implementation of all portions of the model is being completed.
The contaminant dispersal analysis model is based upon the idealization of
building air flow systems as an assemblages of flow elements connected to
discrete system nodes corresponding to well-mixed air zones within the
building and its HVAC system. Equations governing the air flow processes in
the building (e.g., infiltration, exfiltration, HVAC system flow, & zone-to-zone
flow) and equations governing the contaminant dispersal due to this flow,
accounting for contaminant generation or removal, are formulated by
assembling element equations so that the fundamental requirement of
conservation of mass is satisfied in each zone. The character and solution of
the resulting equations is discussed and steady and dynamic solution methods
outlined.
KEY WORDS: contaminant dispersal analysis , flow simulation, building
simulation, building dynamics, computer simulation techniques, discrete
analysis techniques,
-------
NBS: Indoor Air Quality Model ACKNOWLEDGEMENTS
Phase II Report
ACKNOWLEDGEMENTS
Although the author of this report assumes full responsibility for the contents of
the report it is important to acknowledge the contribution made by Richard Grot
and George Walton who together with the author acted, in effect, as a project
team.
Dr. Richard Grot of the Indoor Air Quality and Ventilation Group, Building
Environment Division, National Bureau of Standards closely supervised all
research reported in this document, providing essential critical evaluation and
guiding the direction of the work by applying his considerable experience in the
field and keen intellect to the task at hand. This he accomplished with his
always engaging sense of humor and tireless enthusiasm.
The indoor air quality model presented in this report is based largely upon the
work of George Walton of the Mechanical Systems and Controls Group,
Building Environment Division, National Bureau of Standards. In fact, the
present model should, properly, be presented as an extension of his earlier
work. George was involved in Phase I and the early part of Phase II of this
project and continued thereafter to provide his invaluable insight in the model
development effort.
This work was supported by the Department of Energy, Contract IAG: DE-AI01-
86CE21013, the Environmental Protection Agency, Contract IAG: DW-
13931103-01-2, and the Department of Architecture, Cornell University.
iii
-------
NBS: Indoor Air Quality Model CONTENTS
Phase II Report
CONTENTS
ABSTRACT ii
ACKNOWLEDGEMENTS Hi
CONTENTS iv
PREFACE . vii
PART I - Theoretical Basis
1. General Considerations 1-1
1.1 Definition of Problem 1-2
1.2 Modeling Approaches 1-4
1.3 The Well-Mixed Macroscopic Model 1-6
2. Contaminant Dispersal Analysis . . . 2-1
2.1 Element Equations 2-1
2.2 System Equations 2-3
2.3 Boundary Conditions . 2-5
2.4 Elimination of Massless DOFs 2-7
2.5 Qualitative Analysis of System Equations 2-8
2.6 Solution of System Equations 2-17
2.6.1 Steady State Behavior 2-18
2.6.2 Free Response Behavior 2-18
2.6.3 Dynamic Behavior 2-19
3. Air Flow Analysis 3-1
3.1 Pressure Variation within Zones 3-1
3.2 Element Equations 3-3
3.2.1 Flow Resistance Element Equations 3-3
3.2.2 Fan/Pump Element Equations 3-14
3.3 System Equations 3-18
3.4 Simple Examples 3-20
3.5 Solution of Flow Equations 3-26
3.5.1 Successive Substitution 3-26
3.5.2 Newton-Raphson Iteration 3-30
3.5.3 Incremental Formulation 3-33
4. Summary and Directions of Future Work 4-1
iv
-------
NBS: Indoor Air Quality Model CONTENTS
Phase II Report
PARTI References Refl-1
PART II - CONTAM86 Users Manual
5. General Instructions 5-1
6. Command Conventions 6-1
7. Introductory Example 7-1
8. Command Reference 8-1
8.1 Intrinsic Commands 8-1
8.1.1 HELP 8-1
8.1.2 ECHO 8-1
8.1.3 LIST 8-1
8.1.4PRINTA= 8-1
8.1.5 DIAGRAM A= 8-1
8.1.6 SUBMIT F= 8-1
8.1.7 RETURN 8-2
8.1.8 QUIT 8-2
8.2 CONTAM86 Commands 8-3
8.2.1 FLOWSYS 8-3
8.2.2 FLOWELEM 8-4
8.2.3 STEADY 8-4
8.2.4 TIMECONS 8-5
8.2.5 Dynamic Analysis 8-7
8.2.5.1 FLOWDAT 8-7
8.2.5.2 EXCITDAT 8-9
8.2.5.3 DYNAMIC 8-10
8.2.5.4 Dynamic Analysis Example 8-11
8.2.6 RESET 8-13
9. Example Problems 9-1
9.1 Single Zone Examples 9-1
9.1.1 Case 1: Contaminant Decay under
Steady Flow Conditions 9-2
9.1.2 Case 2: Contaminant Decay under
Unsteady Flow Conditions 9-7
9.1.3 Case 3: Contaminant Dispersal Analysis of an
-------
NBS: Indoor Air Quality Model CONTENTS
Phase II Report
Experimental Test 9-10
9.2 Two Zone Example 9-14
9.3 Full-Scale Multi-zone Residential Example .... 9-19
PARTH References Refll-1
Appendix- FORTRAN 77 Source Code
vi
-------
NBS: Indoor Air Quality Model PREFACE
Phase II Report
PREFACE
The work reported here is a product of the Gieneral Indoor Air Pollution
Concentration Model Project initiated in 1985 at the National Bureau of
Standards with the support of the U. S. Environmental Protection Agency and
the U.S. Department of Energy. The fundamental objective of this project is to
develop a comprehensive validated computer model to simulate dynamic
pollutant movement and concentration variation in buildings. The scope of the
project is ambitious; a full-scale, multi-zone building contaminant dispersal
model that simulates flow processes (e.g., infiltration, dilution, & exfiltration) and
contaminant generation, reaction, and removal processes is being developed.
During the planning stage of this project it was decided to organize efforts into
three distinct phases:
Phase I: formulation of a general framework for the development of general
indoor air quality analysis models (see [1] for report of Phase I work),
Phase II: development of a residential-scale model, based on the simplifying
assumption that air is well-mixed within each building zone, providing
simple simulation of HVAC system interaction, and
Phase III: extension of modeling capabilities to allow more complete simulation
of HVAC system interaction and consideration of rooms that are not
well-mixed.
This report presents a model that satisfies the scope and objectives set for
Phase II of the "General Indoor Air Pollution Concentration Model" Project and,
as such, completes Phase II efforts. The report is organized in two parts. In the
first part of the report the theoretical basis of the model is presented;
Section 1: outlines the general aspects of indoor air quality simulation making
the distinction between contaminant dispersal analysis and air flow
analysis,
Section 2: presents the theoretical basis of contaminant dispersal analysis,
Section 3: presents the theoretical basis of air flow analysis.
VII
-------
NBS: Indoor Air Quality Model PREFACE
Phase II Report
The second part of the report presents the practical implementation of the
contaminant dispersal analysis model in the program CONTAM86;
Sections 5 -8: provide a users manual for the program CONTAM86, and
Section 9: gives examples of application of CONTAM86, and its underlying
theory, to problems of contaminant dispersal analysis.
The complete source code for CONTAM86 is listed in the appendix.
VIII
-------
NBS: Indoor Air Quality Model
Phase II Report
PART I Theoretical Basis
1. General Considerations
1. General Considerations
Airborne contaminants introduced into a building disperse throughout the
building in a complex manner that depends on the nature of air movement in-to
(infiltration), out-of (exfiltration), and within the building system, the influence of
the heating ventilating and air conditioning (HVAG) systems on air movement,
the possibility of removal, by filtration, or contribution, by generation, of
contaminants by the HVAC system, and the possibility of chemical reaction or
physical-chemical reaction (e.g., adsorption or absorption) of contaminants with
each other or the materials of the buildings construction and furnishings.
Fig. 1.1 Contaminant Dispersal in a Residence
Our immediate objective, here, is to develop a model of this dispersal process
for residential-scale building systems that comprehensively accounts for all of
these processes that affect the actual contaminant dispersal phenomena. We
shall, however, attempt, to develop this residential-scale modeling capability
within a more general context so that techniques developed here may be
extended to more complex problems of indoor air quality analysis. To this end,
in this section, the problem is given a general definition and the basic modeling
strategy used to address this problem is outlined.
1 -1
-------
NBS: Indoor Air Quality Model
Phase II Report
PART I Theoretical Basis
1. General Considerations
1.1 Definition of Problem
The building air flow system may be considered to be a three dimensional field
within which we seek to completely describe the state of infinitesimal air
parcels. The state of an air parcel will be defined by its temperature, pressure,
velocity, and contaminant concentration (for each species of interest) - the state
variables of the indoor air quality modeling problem.
Air Parcel State Variables
Temperature: T(x,y,Z,t)
Pressure; P(x,y,Z,t)
Flow Velocity: v(x,y,z,t)
Concentration: C(x,y,Z,t), C(x,y,Z,t), ...
for species
-------
NBS: Indoor Air Quality Model
Phase II Report
PART I Theoretical Basis
1. General Considerations
[=] mass of species/mass of air
a, p = species type indices
x, y, z = spacial coordinates
t = time
and shall refer to the process of determining the spacial and temporal variation
of these species concentrations as contaminant dispersal analysis .
Contaminant dispersal analysis, for a single nonreactive species "a", depends
on the air velocity field and its variation with time;
"CCx.y.Z.t) = "(X v(x,y,Z,t) ) & B.C. : Contam. Dispersal Anal.
(1.1)
But the air velocity field depends on the pressure field which is affected by the
temperature field through buoyancy and, completing the circle, the temperature
field is dependent on the velocity field;
v(x,y,z,t) = v( P(x,y,z,t)) & B.C.
|
P(x,y,z,t) = P( T(x,y,z,t)) & B.C.
T(x,y,z,t) = T( v(x,y,z,t)) & B.C.
Flow Analysis
Bucyancy Effects
Thermal Analysis
(1.2)
(1.3)
(1.4)
where;
B.C = boundary conditions
v = air flow velocity
P = air pressure
T = air temperature
Thus, in general, contaminant dispersal analysis, for a single nonreactive
species, is complicated by a coupled nonlinear flow-thermal analysis problem.
Therefore, a comprehensive indoor air quality model will eventually have to
address the related flow and thermal problems.
For cases of reactive contaminants, contaminant dispersal analysis, itself, will
1-3
-------
NBS: Indoor Air Quality Model PART I Theoretical Basis
Phase II Report 1. General Considerations
become a coupled (and, generally, nonlinear) analysis problem as individual
species' concentrations will depend on other species' concentrations in
addition to the air velocity field;
"(Xx.y.Z.t) = "(Xv, *C, *C, ...) : Species o< Dispersal Analysis (l.5a)
C(x,y,Z,t) = C(V, "C, C, ...) •• Species fl Dispersal Analysis (1.5b)
In this report we shall focus on single, nonreactive species dispersal analysis
and the associated problem of flow analysis, for a completely defined thermal
field and its variation. The approach taken, however, has been formulated to be
compatible with thermal analysis modeling techniques developed earlier [2].
Presently, we are addressing the reactive, multiple species dispersal analysis
problem and see no difficulty with extending the approach to this more complex
situation.
1.2 Modeling Approaches
We shall attempt to solve the general field problems posed above by attempting
to determine the state of air at discrete points in the building air flow system. It
will be shown that this spacial discretization allows the formulation of systems
of ordinary differential equations that describe the temporal variation of the state
fields. Two basic approaches may be considered, one based upon the
microscopic equations of motion (i.e., continuity, motion, and energy equations
for fluids) and the other based upon a "well-mixed" zone simplification of
macroscopic mass, momentum, and energy balances for flow systems (for a
concise and complete review of these basic approaches see [3]).
1 -4
-------
NBS: Indoor Air Quality Model
Phase II Report
PART I Theoretical Basis
1. General Considerations
n
For 3D: 0(1000 nodes/room )
Representative Flow Elements
N
Representative
Nodes
3D: 0(1-10 nodes/room)
Microscopic Model
Macroscopic Model
Fig. 1.3 Basic Spacial Discretization Approaches
In the microscopic modeling approach one of several techniques of the
generalized finite element method, which includes the finite difference method
[4], could be used to transform the systems of governing partial differential
equations into systems of ordinary differential equations that then can be solved
using a variety of numerical methods. The macroscopic modeling approach
leads directly to similar systems of ordinary differential equations.
In both approaches the building air flow system is modeled as an assemblage
of discrete flow elements connected at discrete system nodes . Systems of
ordinary differential equations governing the behavior of elements are then
formed and assembled to generate systems of ordinary differential equations
that describe the behavior of the system as a whole (i.e., in terms of the spacial
and temporal variation of the discrete state variables). These systems of
equations may then be solved — given system excitation, initial conditions, and
boundary conditions — to complete the analysis.
Virtually all computational procedures, except those used to form the element
equations, would be practically identical for both approaches. From a practical
point of view, however, microscopic modeling will involve on the order of 1000
nodes per room while the macroscopic model will involve on the order of only
10 nodes/room to realize acceptably accurate results. With six state variables
1 -5
-------
NBS: Indoor Air Quality Model
Phase II Report
PART I Theoretical Basis
1. General Considerations
for a single species - temperature, pressure, three velocity components and
species concentration - the microscopic modeling approach can lead to
extremely large systems of equations that therefore limit its use, at this time, to
research inquiry. The macroscopic approach, resulting in systems of equations
that are on the order of two magnitudes smaller than the microscopic approach,
is a reasonable candidate for practical analysis, although it can not provide the
detail of the microscopic approach.
Within this report we shall limit consideration to the macroscopic approach,
although the specific techniques employed to implement this approach have
been formulated to be compatible with the microscopic approach and it is
expected that one may, in the future, be able to use both approaches in analysis
to gain the benefits of detail in specific areas of the building system and yet
account for full-system interaction.
Fia. 1.4 Possible Hybrid Micro-Macro Discretization
1.3 The Well-Mixed Macroscopic Model
Here, the building air flow system shall be modeled as an assemblage of flow
elements connected to discrete system nodes corresponding to well-mixed
air zones.
1-6
-------
NBS: Indoor Air Quality Model
Phase II Report
PART I Theoretical Basis
1. General Considerations
• — Zone Node
2 — Node Number
Flow Element-
Element Flow
Actual Building
Air-Flow System Idealization
Fig. 1.5 Well-Mixed Macroscopic Model
Limiting our attention to the contaminant dispersal and flow analysis problems
we associate with each system node the discrete variables or degrees of
freedom (DOFs) of pressure, air mass generation (typically zero), species
concentration, species mass generation, and temperature;
(P) = (Pi, P2, P3, - )
{W} = {W,, W2, W3, ... }
{"O = {"c,, txc2, %,
( ^3>
(T) = {T,, T2, T3> ... }
: Pressure DOFs
-. Air Mass Generation DOFs
-. Species
-------
NBS: Indoor Air Quality Model PART I Theoretical Basis
Phase II Report 1. General Considerations
element air mass flow rate, we. The element mass flow rates will be related to
the nodal state variables through specific properties associated with each
particular element to form element equations.
In the formulation of both the contaminant dispersal model, presented in
Section 2, and the flow model, presented in Section 3, we will assemble the
governing element equations to form equations governing the behavior of the
building system - the system equations - by demanding conservation of mass
flow at each system node.
1 -8
-------
NBS: Indoor Air Quality Model
Phase II Report
PART I-Theoretical Basis
2. Contaminant Dispersal Analysis
2. Contaminant Dispersal Analysis
In this section contaminant dispersal element equations are formulated.
Demanding continuity of mass flow at each system node these element
equations are then assembled to form contaminant dispersal equations
governing the behavior of the full building system. Finally, methods for solution
of the system equations are presented.
2.1 Element Equations
Two nodes2'1 and a total mass flow rate, we , will be associated with each
flow element, where flow from node i to j is defined to be positive. An element
species concentration, "C®, and an element species mass flow rate, aw£ , will
be associated with each element node, k=i, j. The element species mass flow
rate is defined so that flow from each node into the element is positive.
nodei
W
jlement "e"
node j
Fig. 2 .1 Contaminant Dispersal Element DOFs
It follows from fundamental considerations that these element variables are
related directly to the element total mass flow rate as;
2'1 The distinction between element nodes and systems nodes must be made because the
element species concentration vector, {aCe}, is taken as a subset of the system species
concentration vector, {aC}.
2-1
-------
NBS: Indoor Air Quality Model
Phase II Report
PART I - Theoretical Basis
2. Contaminant Dispersal Analysis
or
_
;forwe>0
(2.1 a)
-1~l ,v Q
{aCG} ;forwe<0
{"w8} =
where;
\ = (V , X >T
{ace> = {acr, ace}T
(2.
; element species mass flow rate vector
; element species concentration vector
[f *] = element total mass flow rate matrix
;forwe>0
;forwe<0
For the purposes here, element nodes will be selected to correspond to specific
system nodes, consequently, the element nodal species concentrations will
have a one-to-one correspondence with the corresponding system node
species concentrations.
If the element acts as a filter and removes a fraction, T\ , of the contaminant
passing through the filter then the element flow rate matrix becomes;
[f *] = element total mass flow rate matrix
= w
1 0
L(T1-D Oj
; for we > 0
= w
"0 (Ti-1)l .
.0 1 J'
for we < 0
2-2
-------
NBS: Indoor Air Quality Model PART I - Theoretical Basis
Phase II Report 2. Contaminant Dispersal Analysis
The fraction, TJ , is commonly known as the "filter efficiency" and may have
values in the range of 0.0 to 1.0.
2.2 System Equations
System equations that relate the system concentration DOFs, {aC}, to the
system generation DOFs,{aG}, may be assembled from the element equations
by first transforming the element equations to the system DOFs and then
demanding conservation of species mass flow at each system node.
There exists a one-to-one correspondence between each element's
concentration DOFs, {aCe}, and the system concentration DOFs, {aC}, that may
be defined by a simple Boolean transformation;
{aCe} = [aB8]{aC} (2.2)
where;
[ B ] is an m x n Boolean transformation matrix consisting of zeros
and ones; m = the number of element nodes (here, m=2); n =
the number of system nodes
For example, an element with nodes i & j (or 1 & 2 ) connected to system nodes
5 & 9, respectively, of a 12-node system would have ones in the 1st row, 5th
column and the 2nd row, 9th column and all other elements of the 2x12
Boolean transformation matrix would be set equal to zero.
In a similar manner, we may define a "system-sized vector" to represent the net
species mass flow rate from the system node into an element "e", {aWe} , and
relate it to the corresponding element species mass flow rate using the same
transformation matrix, as;
(2 3)
For an arbitrary system node n, with connected elements "a", "b",... as indicated
below in Fig. 2.2, we then demand conservation of species mass as;
2-3
-------
NBS: Indoor Air Quality Model
Phase II Report
PART I - Theoretical Basis
2. Contaminant Dispersal Analysis
(elem. species mass flow) +
connected
elements
'rateofchange"!
•of
species mass J
generation
of
species mass
(2.4)
system node n
or,
... + vr
d°C
n
(2.5)
or, for the system as a whole;
X^ a
(2-6)
where;
M
Vi
= diag(V1, V2, ...) ; the system volumetric mass matrix
= the volumetric mass of node i
system node "n"
direct species mass generation
.element "c" species mass flow
Fia. 2.2 Conservation of Species a Mass Flow at System Node n
2-4
-------
NBS: Indoor Air Quality Model PART I - Theoretical Basis
Phase II Report 2. Contaminant Dispersal Analysis
Substituting relations (2.2) and (2.3) we obtain the final result;
(2.7a)
where;
[F] = X [aBe]T[fe][aBe] (2.7b)
e = a,b, ...
= the system mass flow matrix
= A[fe] ; the direct assembly sum of element flow matrices
Equation (2.7a) defines the contaminant dispersal behavior of the system as a
whole and is said to be assembled from the element equations through the
relation given by equation (2.7b). The assembly process, as formally
represented in equation (2.7b), has found widespread application in the
simulation of systems governed by conservation principles and is, therefore,
often represented by the so-called assembly operator A as indicated above. It
should be noted that while the formal representation of the assembly process is
important from a theoretical point of view it is generally far more efficient,
computationally, to assemble the element equations directly, without explicitly
transforming them ( see, for example, the "LM Algorithm" in [24]).
2.3 Boundary Conditions
The variation of concentration or generation rate, but not both, may be specified
at system nodes. Concentration or generation conditions in the discrete model
are equivalent to boundary conditions in the corresponding continuum model
and will, therefore, be referred to as such.
Formally then, we may distinguish between those DOFs for which concentration
will be specified, {aCc}, and those for which generation rate will be specified,
{aCg}, and partition the system of equations accordingly;
2-5
-------
NBS: Indoor Air Quality Model
Phase II Report
PART I - Theoretical Basis
2. Contaminant Dispersal Analysis
FCC Fqj
LFgc Fggj
ac
cc
d°Cc
dt
daC
g
dt
aGgJ
(2.8)
Using the second equation and simplifying we obtain;
[Faa]{aCa> + [Vaa]
gg
dt
= {aGa} - [Fac]{aCc}
or
[F]{aC}
= {aE}
(2.9a)
(2.9b)
where;
A
[F] = [Fgg] ; the generation driven mass flow matrix
{aC} = {aCg} ; the generation driven nodal concentration vector
{"£} s {aGg} - [Fgc]{aCc} ; the system excitation (2.9c)
It should be noted that the response of the system is driven by the system
excitation involving both specified contaminant mass generation rates and
contaminant concentrations which may, in general, vary with time.
Equation (2.9b), written in the standard form of a set of first order differential
equations similar to the form of equation (2.7a), most directly defines the
contaminant dispersal behavior of the system. The formation and solution of
equation (2.9b^ will be considered the central task of contaminant dispersal
analysis.
The response of the system is defined by the solution of equation (2.9b) for the
generation rate specified DOFs, {aCg}. The generation rates, {aGc}, required to
maintain the specified concentrations, {aCc}, may be determined from the
response of the system to the specified excitation using the first equation of
2-6
-------
NBS: Indoor Air Quality Model
Phase II Report
PART I - Theoretical Basis
2. Contaminant Dispersal Analysis
equation (2.8) as;
{aGc> =
+ [Fcg]{aCg}
d
(2.10)
Alternatively, one may numerically imposed specified concentration conditions
by directly modifying equation (2.7a). The effect of an infinite source or sink, of
the desired concentration, may be effected by scaling the appropriate diagonal
terms of the system matrices by a large number and setting the corresponding
generation rates equal to the product of the specified concentration and the
scaled diagonal term. (The current version of CONTAM uses this strategy.)
2.4 Elimination of Massless DOFs
Often the analyst will define flow nodes within a complex building airflow system
to model zones having negligibly small volumetric masses (e.g., junctions in
HVAC system ductworks) and the analyst may prefer to model theses zones as
if their nodal volumetric masses were zero. Additionally, the response at such
nodes may be of little interest and the analyst may prefer to eliminate these
nodal DOFs from consideration.
If the system of equations (2.9b) is partitioned into those DOFs having zero
nodal volumetric masses, {aCz}, and those having non-zero volumetric masses,
{aCn}, as;
^zz rz
A A
'
nnj
0 0 1
o vnnj
d°Cz
~dt~
daCn
. dt .
(2.11)
we may eliminate the massless DOFs from consideration by first solving for
these DOFs using the upper equation;
("EZ) - [Fzn]{aCn}
(2.12)
2-7
-------
NBS: Indoor Air Quality Model PART I - Theoretical Basis
Phase II Report 2. Contaminant Dispersal Analysis
and substituting this result in the lower equation to obtain;
[FKaC>
(2.13a)
where;
[F] = [Fnz][Fzz]~1[Fzn] ; the reduced system flow matrix (2.13b)
{"£} ={aEn}-[Fzzr1{aEz} ; the effective system excitation (2.13c)
("C} ^ {aCn}
[V] ^ [Vnn]
Equation (2.13a) is simply a reduced form of equation (2.9b); being a system of
smaller size it may be solved more efficiently. In addition, the elimination of
massless DOFs should help to avoid some numerical problems associated with
round-off error. Eventhough the massless DOFs have been eliminated from
consideration in equation (2.13a) their values may be recovered, at any time,
using equation (2.12). (The current version of CONTAM does not eliminate
massless DOFs.)
2.5 Qualitative Analysis of System Equations
It is important to keep in mind that we have developed equations that described
the contaminant dispersal behavior of building idealizations, based upon
assemblages of ideal flow elements, and have not, strictly speaking, developed
equations that govern the behavior of the actual buildings being considered.
Although it is hoped that these building idealizations will accurately describe
the behavior of the actual buildings being modeled it is possible that they will
not. In fact, it is quite possible to create idealizations that result in equations that
have no solution, at all.
In this section, therefore, we shall consider the conditions that must be met to
yield contaminant dispersal equations that have solutions and in so doing we
shall also learn something about the general qualitative character of the
solutions that are possible.
2-8
-------
NBS: Indoor Air Quality Model PART I - Theoretical Basis
Phase II Report 2. Contaminant Dispersal Analysis
It should come as no surprise that building idealizations that satisfy
conservation of total mass flow (i.e., as distinguished from species mass flow)
will lead to system of equations that do, in fact, have solutions, but to get to this
seemingly obvious conclusion we shall have to consider the details of the
system flow and mass matrices and their impact upon the dynamic character of
the system as a whole.
System Flow Matrix
The system flow matrix [F], being a direct assembly sum of nonsymmetric
element matrices, will also, in general, be nonsymmetric. The details of the
assembly process reveal that the diagonal elements of the flow matrix are
always positive and the off-diagonal elements negative. Furthermore, if the total
mass flow into a system node is equal to the total mass flow out of a system
node, then the diagonal elements of the flow matrix will be less than or equal to
the "row sum" or the "column sum" of the corresponding off-diagonal elements.
More specifically, for a given system node i the diagonal element, FJJ , is simply
equal to the total mass flow out of a node, theow sum of row i equals the sum
of total mass flow inla the node weighted by the filter efficiency factors (TI -1);
row sum of row i = > | Fy | = weighted total mass flow into node i (2.14)
j* i
and the column sum equals the sum of total mass flow out of the node weighted
by the filter efficiency factors (TI -1);
n
column sum of col. i = /, | FJJ = weighted total mass flow out of node i (2.15)
J=1
Therefore, if total mass flow is conserved at each node, we may assert;
2-9
-------
NBS: Indoor Air Quality Model PART I • Theoretical Basis
Phase II Report 2. Contaminant Dispersal Analysis
y
JJ > F = row sum of row i (2.1 6)
and
Fii ^ Fi = column sum of col. i (2.17)
where the equality is strict when filter efficiencies of the elements connected to
node i are zero (i.e., all TI = 0) and the inequality holds if any of the connected
outflow elements (for the row sum) or inflow elements (for the column sum) have
nonzero filter efficiencies.
If all elements of a flow system idealization have nonzero filter efficiencies then
the system flow matrix will be strictly diagonally dominant (i.e., for all i the
inequalities above will hold); a condition that insures, by itself, the possibility of
solution; that is to say, a sufficient condition to prove that the flow matrix would
be nonsingular. For the (unlikely) limiting case where all elements have filter
efficiencies equal to 1.0 the flow matrix becomes diagonal and, therefore, all
zones act as independent (i.e. .uncoupled) single zone systems.
At the other (more likely) extreme where all elements have filter efficiencies
equal to 0.0 the equalities of equations (2.16) and (2.17) hold for all nodes and
the flow matrix is no longer strictly diagonally dominant and, therefore, may not
be assumed to be nonsingular. We may show, however, that the important
submatrix of the flow matrix identified earlier as the generation driven mass flow
matrix is, in fact, nonsingular by demanding conservation of total mass flow of
all subassemblages of system nodes and their inter-connecting elements and
using some relatively esoteric theorems relating to the general class of matrices
known as M-matrices.
An M-matrix may be defined in a number of alternative, but equivalent ways.
Using the alternative employed by Funderlic and Plemmons [5] an M-matrix is a
square nonzero real matrix with all off-diagonal elements nonpositive that has
2-10
-------
NBS: Indoor Air Quality Model PART I - Theoretical Basis
Phase II Report 2. Contaminant Dispersal Analysis
eigenvalues with nonnegative real parts. It may be shown [6] that a real square
matrix [A], with positive diagonal elements and nonpositive off-diagonal
elements;
a) is an M-matrix (possibly singular) if and only if it can be shown that
[ [A] + £[l] ] is a nonsingular M-matrix for all scalars £ > 0 and
b) is a nonsingular M-matrix if [A] is strictly diagonally dominant
In the case at hand, clearly [ [F] + £[l] ] is strictly diagonally dominant, and
therefore a nonsingular M-matrix, for all scalars £ > 0, (if, of course, total mass
flow is conserved at all nodes). Thus we can conclude that [F] is an M-matrix,
although it will be singular for the limiting case when all filter efficiencies are
zero.
It has also been shown that each principal submatrix of an irreducible M-matrix
(other than the M-matrix itself) is a nonsingular M-matrix [7]. The flow matrix
would be said to be reducible if it is possible, using an appropriate numbering
of the system nodes, to assemble the flow matrix in the form;
[F] =
FH F,
0 F
22.
(2.18)
where F-n and F22 are square matrices, otherwise [F] would be said to
irreducible. Recalling thatsuperdiagonal term, FJJ ; j > i, corresponds to flow
from node j to node i and a subdiagonal term, FJJ ; j > i, corresponds to flow from
node i to node j, a flow matrix of the form of equation (2.18) would correspond to
a flow system idealization having a total mass flow from subassembly 2 to
subassembly 1, without a return flow from 1 to 2, and, therefore, conservation of
total mass flow would be violated.
We may conclude, then, that;
a) the flow matrix, [F], will be an irreducible M-matrix and, therefore,
b) the generation driven mass flow matrix, [F] , a principal submatrix of the
flow matrix will be a nonsingular M-matrix.
2-11
-------
NBS: Indoor Air Quality Model PART I - Theoretical Basis
Phase II Report 2. Contaminant Dispersal Analysis
if they are formed based upon a flow idealization that satisfies conservation of
total mass flow
Inasmuch as the solution of the generation driven contaminant dispersal
equations (equation (2.9b)) is the central task of contaminant dispersal analysis
and the nonsingularity of the generation driven flow matrix is a necessary
perequisite to assure the possibility of solution of these equations, the
conclusion that the generation driven flow matrix will be nonsingular when the
flow system idealization satisfies the condition of total mass conservation is of
paramount importance. An additional property of nonsingular M-matrices
provides the additional benefit of allowing efficient numerical solution strategies
to be employed in the solution of these equations.
^
Nonsingular M-matrices, and therefore, properly formed [F] matrices, have the
important additional property that they may be factored into the product of lower,
A
[L], and upper, [U], triangular matrices, [F] = [L][U] , by Gauss elimination
without the need of pivoting in an efficient and numerically stable manner (i.e.,
resulting in no more accumulation of error that that which would result if pivoting
were employed) [8]. Therefore, not only may we be certain that a properly
formed flow matrix will lead to the possibility of solution but it will also allow the
advantage of the use of very efficient methods of solution associated with LU
decomposition.
System Volumetric Mass Matrix
By definition the system volumetric mass matrix, [V], is diagonal and
nonnegative. In those instances when some nodal volumetric masses are so
small that the analyst prefers to modeled them with zero values the system of
contaminant dispersal equations may be reduced, by eliminating the massless
equations (see section 2.4), to a form having an all positive, and therefore,
nonsingular, volumetric mass matrix. The inversion of the positive volumetric
mass matrix is trivial;
[V]-1 = diagd/VLlA/2,... 1/Vn) ;V^O (2.19)
2-12
-------
NBS: Indoor Air Quality Model PART I - Theoretical Basis
Phase II Report 2. Contaminant Dispersal Analysis
System Equations - Steady Flow
The generation driven contaminant dispersal equations, equation (2.9b), may
now be rewritten in the form;
«
fl LJ I A H f* A
i = [V]-1[aE] (2.20)
A
where, in general, the [F] will vary with time.
A 4 A
The product matrix [V]~ [F] contains the essential dynamic character of the
system being studied. For properly formed idealizations (being the product of a.
positive diagonal matrix and a nonsingular M-matrix [9]) it will be a nonsingular
M-matrix and, therefore,
a) solutions to equation (2.20) will exist, and
b) the product matrix may also be factored into the product of lower, [L], and
upper, [U], triangular matrices, [V]'1[F] = [L][U] , by Gauss elimination
without the need of pivoting in an efficient and numerically stable manner.
We may gain some insight into the general character of solutions to equation
(2.20) by considering the case of steady flow ( [F] constant) without excitation
(i.e., the homogeneous case);
[Vr1[F]{aC} + = [0] (2.21)
Anticipating the result we try solutions of the form;
{«C} = {"Ole^ (2.22)
where;
T = decay time constant
2-13
-------
NBS: Indoor Air Quality Model
Phase II Report
PART I-Theoretical Basis
2. Contaminant Dispersal Analysis
{"O } = vector of unknown magnitudes
which, when substituted into equation (2.21) lead to the standard eigenvalue
problem;
[[V]-1[F] -
= {0}
(2.23)
The solution of this standard eigenvalue problem and its relation to the first
order system of differential being considered is discussed elsewhere [10], [11]
and is well beyond the scope of this report. Suffice it to say, for a properly
formed flow system idealization of n nodes there will be n solutions to this
eigenvalue problem consisting of n pairs of time constants, T, (or equivalently
their inverses, 1/i - the system eigenvalues) and their associated eigenvectors,
In some cases it may be possible to transform the product matrix [V]~ [F], by
similarity transformations, to diagonal form leaving the eigenvalues on the
diagonal as;
[Sr1[[Vf1[F]][S] =
where;
'(1/T,) 0 ... 0
0 (1/x2) ... 0
0 0 ... (1/tn)
(2.24)
[S] = the similarity transformation
For these cases it will be possible to express the general solution to the
homogeneous problem, equation (2.21), as a linear combination of simple
exponential decay terms;
(2.25)
where the scalar coefficients, a-|, a2, ... an, are determined from the initial
conditions using the similarity transformation employed as;
2-14
-------
NBS: Indoor Air Quality Model PART I - Theoretical Basis
Phase II Report 2. Contaminant Dispersal Analysis
a2
• = [sr1-
aC2(t=0)
(2.26)
The n pairs of time constants and associated eigenvectors are often referred to
as the system modes and the response of the system is often described in
terms of the degree to which each mode participates. From the form of the free
response, equation (2.25), it is clear that as time passes the contribution of
those modes with larger time constants will dominate the character of the
response until, eventually, the response, in all zones, will be dominated by the
mode with the largest time constant and therefore will appear to be a simple
exponential decay.
The similarity transformation [S] may be chosen as a matrix whose columns
equal the eigenvectors, in this case, and, therefore, by equation (2.26) we can
see that we may trigger a decay response in any single mode if we simply set
the initial conditions equal to the corresponding eigenvector (or a scalar
multiple of it), although, for some modes the eigenvectors will have negative
components that, for contaminant dispersal problems, would not be physically
admissible.
In general, the solution of the eigenvalue problem will be computationally
demanding. However, for the limiting case discussed earlier, when all flow
elements have filter efficiencies equal to 1.0, eigenanalysis is trivial. For this
case the product matrix [V]~1[F] will be diagonal, therefore;
a) the time constants, TJ, will be simply equal to (V/Fjj),
b) the similarity transformation will be equal to the identity matrix,
c) the eigenvectors will be equal to the unit vector corresponding to each
DOF (i.e., the columns of the identity matrix), and
d) the scalar coefficients will equal the initial conditions corresponding to
each DOF { a1f a2,... an} = {aC1(t=0)> aC2(t=0),... aCn(t=0)}.
2-15
-------
NBS: Indoor Air Quality Model PART I - Theoretical Basis
Phase II Report 2. Contaminant Dispersal Analysis
For this limiting case all zones act independently as single zone "systems" and,
therefore, these results follow directly from the more familiar single-zone theory.
For general contaminant dispersal systems we may apply the Gerschgorin
Theorem [10], given the volumetric mass matrix is diagonal, to obtain a poorly
bounded, but computationally inexpensive, estimate of the (real part of) system
time constants as;
:foralli (2-27)
'
This expression simplifies, exactly, to the values obtained for the limiting case
discussed above, when all filter efficiencies equal 1.0, while at the other
extreme, when all filter efficiencies are 0.0, it assures only that the system time
constant will fall within the range;
al1 f'lter efficiencies = 0.0 (2.28)
as, in these cases the off-diagonal row sum will be equal to the diagonal value
of the flow matrix.
In some cases it will not be possible to diagonalize the product matrix [V]~1[F],
but in these cases it will always be possible to transform the product matrix to a
form known as the Jordan canonical form, an upper block-triangular matrix with
the eigenvalues (inverse time constants) on the diagonal. For these cases, it
will still be possible to express the general solution to the homogeneous
problem, equation (2.21), as a combination of exponential decay terms, but now
some of these decay terms will have factors equal to powers of time (i.e., in
addition to terms like eW we will have to include terms like teW , t2e-(t/T) ,
t3e-(t/T) i6tc.).
In all cases the system time constants will have positive real parts, as the
product matrix is a nonsingular M-matrix, and therefore all components making
up the general solution will approach zero with time. That is to say, the
2-16
-------
NBS: Indoor Air Quality Model PART I - Theoretical Basis
Phase II Report 2. Contaminant Dispersal Analysis
homogeneous contaminant dispersal equations are stable; the concentration
at all nodes will (eventually) approach zero. Furthermore, following the
argument similar to that presented earlier in the discussion of the flow matrices,
we may show that the sum of the product matrix and its transpose;
[[VT1[F] + [[Vr1[F]]T]
is also a nonsingular M-matrix with positive (real parts of) eigenvalues and,
therefore, the sum of the squares of the system concentrations (i.e., the
Euclidean norm of the concentration vector) will decay at every instant of time
[12];
;t,0 (2.29)
.
where;
H{aC(t)}||2 = (I^COI2 + fC2(t)|2 + ... |aCn(t)|2)
These results are consistent with experience (and intuitive expectation) that
while some nodal concentrations may at first increase with time (e.g., due to
zone-to-zone mixing) in the long run all concentrations will diminish toward the
zero level and at all times (some reasonable measure of) the mean
concentration will also be diminishing.
The response of steady flow systems to nonzero excitation (i.e., the
inhomogeneous case) may also be expressed in terms of linear combination of
— 1 A
the eigenvectors of the product matrix [V] [F]. For practical contaminant
dispersal analysis, however, it is more convenient to solve the system equations
directly using numerical integration techniques that are not limited to steady
flow cases.
2.6 Solution of System Equations
The governing system of equations, equation (2.9b), have the form of a system
of first order linear differential equation with constant coefficients. In many
practical situations, however, the mass flow rates will not be constant in time,
and thus, in general, we may consider equation (2.9b) to be a system of first
order differential equations with nonconstant coefficients. Here we shall
consider the solution of these equations for;
2-17
-------
NBS: Indoor Air Quality Model PART I - Theoretical Basis
Phase II Report 2. Contaminant Dispersal Analysis
1) Steady State: steady contaminant generation rates under conditions of
steady element mass flow,
2) Free Response: transient decay of contaminant concentration under
conditions of steady element mass flow,
3) Dynamic Response: to steady flow with unsteady generation rates, to
unsteady flow with steady generation rates, or to unsteady flow with unsteady
generation rates.
In the discussion below, equation (2.9b) will be written dropping the hat, A, to
simplify notation.
2.6.1 Steady State Behavior
For systems with steady element mass flows driven by steady contaminant
generation rates and/or specified concentrations the response of the system
will, eventually, come to a steady state (i.e., {daC/dt} = 0 ) given by the solution
of;
[F]{a.C> = {"£} (2.30)
As discussed in section 2.5 above this equation may be solved by LU
decomposition without pivoting in an efficient and numerically stable manner.
2.6.2 Free Response Behavior
The free response behavior of steady flow systems has been discussed above
and shown to be closely related to the solution of the eigenproblem given by
equation (2.23) that yields system time constants and associated eigenvectors.
For steady flow systems knowledge of the system time constants provides
invaluable insight into the dynamic character of the system yet eigenanalysis is
computationally time consuming. It is, therefore, tempting to estimate the
system time constants, after single-zone theory, by the ratio of the volumetric
2-18
-------
NBS: Indoor Air Quality Model PART I - Theoretical Basis
Phase II Report 2. Contaminant Dispersal Analysis
mass of each zone to the total air flow out of the zone. This estimate of system
time constants will be designated as the nominal system time constants and,
from the discussion in section 2.5, may be represented as;
TJ ~ •=- ;the nominal system time constants (2.31)
hi
For typical situations, however, the error bound on this estimate is very large
(see section 2.5) and this estimate of the actual system time constants is likely to
be a very poor estimate.
A variety of techniques exist that will provide better solutions to the governing
eigenvalue problem and thereby provide better estimates of the actual system
time constants [13]. The program CONTAM uses a relatively simple, published
procedure, based on Jacobi iteration, that transforms the product matrix,
[V]~1[F], to upper triangular form leaving the eigenvalues on the diagonal [14].
(The command TIMECONS in the program CONTAM reports both nominal and
actual time constants for comparative purposes.)
2.6.3 Dynamic Behavior
The governing systems of equations, equation (2.9b), may be solved for cases
of steady flow with general unsteady contaminant generation using any number
of different finite difference solution schemes. Here we shall employ a general
form predictor-corrector method.
For cases of unsteady flow it is likely that this same predictor-corrector solution
scheme will prove useful, providing, of course, the system flow matrix, [F], is
updated appropriately, although for cases of rapidly changing flow rates small
time steps may be required to control error. If difficulties arise, an iterative
scheme may have to be nested within the predictor-corrector time integration
scheme.
A finite difference scheme for the approximate integration of the semidiscrete
equation (2.9b) may be developed by dividing time domain into discrete steps;
tn+1 = *n + 5t ; n = 0,1,2,3... (2.15)
2-19
-------
NBS: Indoor Air Quality Model PART I - Theoretical Basis
Phase II Report 2. Contaminant Dispersal Analysis
t0 = initial time
where;
8t = integration time step (often constant but may be variable)
demanding the satisfaction of equation (2.9b) at each of these steps;
[F]{aC}n+1 + Mp = (aE}n+i (2.33)
I at J n+1
where;
{aC}n+1 =
Jdacl f daC(tn+1)|
1 dt Jn+1 1 dt J
(aE}n+1 = {aE(tn+1)}
Substituting into this equation the consistent difference approximation
represented by;
{aC}n + (1-e)5t + e (2.34)
I at J n I at J n+1
where;
0<6<1
6 = 0 corresponds to the Forward Difference scheme
9=1/2 corresponds to the Crank-Nicholson scheme
0 = 2/3 corresponds to the Galerkin scheme
6 = 1 corresponds to the Backward Difference scheme
a general implicit finite difference scheme is formulated;
65t[F] + [V] ] p - (aE}n+1 - [F] {aC}n + (1+9)5t ' (2.35a)
2-20
-------
NBS: Indoor Air Quality Model PART I - Theoretical Basis
Phase II Report 2. Contaminant Dispersal Analysis
or, equivalently;
ek)VIf°C)" + (1-6)5t{ lif} „ <2'35b)
Computationally it is useful to implement this general finite difference scheme,
equation (2.35), as a three step predictor-corrector algorithm;
(aC}n+1 E={aC}n + (l-9)5t ^- -predictor (2.36a)
L «i J p|
[ (65t)[F] + [V] ] | Ijp | - {aE}n+1 - [F]{aC}n ; (i.e. eqn (2.35a) ) (2.36b)
I at J n+i
{aC}n+1 - {aC}n+1 + (66t){^p| ; corrector (2.36c)
I at J n+1
It should be noted that;
a) this algorithm is self-starting; given initial conditions, {aC(t0)} , equation
(2.33) may be solved to obtain an estimate of the initial rate of change of
nodal temperatures, {daC(t0)/dt} , and the first predictor step, equation
(2.36a) may then be computed, and
b) equation (2.36b) may also be solved by LU decomposition, without the
need of pivoting; importantly then, the matrix [ (95t)[F] + [V] ] may first be
factored into the L and U product matrices and need not be refactored again
until there is a change in the system flow matrix (i.e., due to unsteady
element flows) and equation (2.36b) may then be solved, at minimum
computational cost by back and forward substitution using the LU factors, for
the first and each subsequent time step.
This predictor-corrector scheme has been analyzed by Taylor [15] and Huebner
[16] and a more general predictor-multicorrector scheme that includes this
implicit scheme has been analyzed by Hughes [1 7] for systems with constant
coefficient matrices (i.e.. [F] and [V] constant) . For 6 > 1/2 this scheme leads to
an unconditionally stable solution; 9 > 3/4 (approximately) leads to an
unconditionally stable non-oscillatory solution; beyond this, Taylor makes some
2-21
-------
NBS: Indoor Air Quality Model PART I - Theoretical Basis
Phase II Report 2. Contaminant Dispersal Analysis
recommendations regarding selection of 0 and step size, 8t, to limit error while
minimizing computational effort. (In the program CONTAM the default value of 6
is set to 0.75, and may be reset by the user, and an estimate of the time step
needed to limit error is reported (for the given initial conditions) using a method
developed by Taylor [15].)
2-22
-------
NBS: Indoor Air Quality Model
Phase II Report
PART I - Theoretical Basis
3. Air Row Analysis
3. Air Flow Analysis
In this section air flow element equations are formulated that relate mass flow
rate through flow elements to pressure differences across the elements, the
assembly of these element equations to form equations governing the flow
behavior of the building air flow system is discussed, and methods of solving
these equations are presented. The formulation of the air flow equations
presented herein is based, in large part, on the work of Walton [18], an example
presented by Carnahan et. al. [19], and Chapter 33 of the ASHRAE Handbook
1985 Fundamentals [20].
3.1 Pressure Variation within Zones
A general model of building airflow systems, the "well-mixed macroscopic
model", and system DOFs relating to this model were defined in Section 1.3 of
this report. For this model, fluid density within any zone i, pj( will be assumed
constant and thus the variation of static pressure within a zone, pj(z), will be
given by;
P.(Z) = P.
where;
Z.
Z
g
g
. -Z)
..
= the elevation of node i relative to an arbitrary datum
= elevation relative to an arbitrary datum
= the acceleration due to gravity
= dimensional constant (1.0 (kg m)/(N s2) )
(3.1)
3-1
-------
NBS: Indoor Air Quality Model
Phase II Report
PART I-Theoretical Basis
3. Air Flow Analysis
• — Zone Node
2 — Node Number
Flow Element-
Element Flow -
Fig. 3.1 Elevations Defined Relative to a Datum
Static pressures (i.e., under still conditions) acting on exterior surfaces may be
approximated as;
•P z
r
; on exterior surfaces, calm conditions
(3.2)
where Pa and pg are the atmospheric pressure and air density at the level of the
outdoor datum.
To account for pressures due to wind effects the pressure on any exterior
surface may be approximated using published wind pressure coefficients [21]
as;
P I
p(2) = P +C -2-
a p 2
; on exterior surfaces, windy conditions
(3.3)
where Cp is a dimensionless pressure coefficient associated with the position
on the exterior surface and the characteristics of the wind and UH is the wind
speed at the roof level of the building. Usually, local wind data will not be
3-2
-------
NBS: Indoor Air Quality Model PART I - Theoretical Basis
Phase II Report 3. Air Flow Analysis
available; reference [21] suggests one modification of equation (3.3) to allow
use of airport wind speed data.
(Strictly speaking equation (3.2) is exact for only a homogeneous atmosphere,
i.e., of constant density. Typically, however, the lower atmosphere, at the scale
of even the tallest buildings, has characteristics that fall between that of an
isothermal atmosphere and a homogeneous atmosphere and equation (3.2)
provides a very good estimate of air pressure for this range of conditions.
Equation (3.3), on the other hand, provides only very approximate estimates of
surface pressures. This is due to the great uncertainty of both pressure
coefficients and the local wind speeds.)
3.2 Element Equations
Two classes of elements will be developed here; the first class, flow resistance
elements, is a very general class that may be used to model a large variety of
flow paths that provide passive resistance to flow (e.g., conduits, ducts,
ductwork assemblies, small orifices such as cracks, etc.); the second class is
developed to model fan-driven air flow. These two classes of elements should
allow modeling of a large variety of complex and complete building airflow
systems. It is anticipated, however, that special elements may need to be
developed, in the future, to provide better models of some flow paths (e.g., flow
through large openings such as doors and windows). Special elements may be
developed using the resistance and fan/pump element formulations as
examples of the general approach of element formulation.
3.2.1 Flow Resistance Element Equations
Resistance to flow will be modeled by flow elements having a single entry and
exit (e.g., simple ducts, openings between zones, orifices, etc.). Flow
components with multiple entries, exits, or both may be modeled as
assemblages of these simpler elements.
Flow resistance elements shall be two-node elements. With each node we
associate element pressure, Pf , temperature, Tf, and flow rate, w®, DOFs (i.e.,
for flow from the node into the element). Element nodes are selected to have
3-3
-------
NBS: Indoor Air Quality Model
Phase II Report
PART I-Theoretical Basis
3. Air Flow Analysis
the same elevation as the zone nodes they connect3'1.
,
m m
m
Zone I
ZoneJ
W
n n
DATUM
Fig. 3.2 Flow Resistance Element DOFs
Fluid flow within each flow resistance element is assumed to be incompressible,
isothermal, and governed by the Bernoulli equation as applied to duct design
[20];
pv2 pv2
(P,
2g,
-1) + -S-pCzf! -
2gc gcp '
Where, for the purposes of developing the general element equations, the more
conventional flow variables, indicated below, have been used;
P-i , P2 = entry and exit pressures, respectively
^1 , V2 = entry and exit mean velocities, respectively
gc = dimensional constant, 1.0 (kg-m)/(N-sec2)
g = the acceleration of gravity (e.g., 9.80665 m/sec2)
p = density of fluid flowing through the element
z1 , z2 = elevations of entry and exits, respectively
we = mass flow rate through the element
3'1 The distinction between element nodes and system nodes must be made because the
element pressure vector, {P8}, is taken as a subset of the system pressure vector, {P}.
3-4
-------
NBS: Indoor Air Quality Model PART I - Theoretical Basis
Phase II Report . 3. Air Flow Analysis
% Ap0 = the sum of all frictional and dynamic losses in the elements
reference
cross section o ...
Fig. 3.3 Conventional. Flow Variables
o
The losses, £ Ap0, are commonly related to the velocity pressure, pV0 / 2gc , of
the fluid flow at reference cross sections "o";
Pv2
where; C0 = loss coefficient
- For conduits of constant cross-section:
= f LVDeq
with;
f = dimensionless friction factor (see Chapter 33
equation (22) and/or Chapter 2 equations (16), (17), &
(1 8) of ASHRAE 1985 Handbook of Fundamentals)
~ constant for turbulent flow (i.e. Re > 2 x 1 03)
=64/Re for laminar flow (i.e. Re < 2 x 1 03)
Re = V0Deq/n
ji = the fluid viscosity
L = length of conduit
Deq = equivalent diameter of conduit
= 4A/PW = 4(flow area)/(wetted perimeter)
- For "fittings" of air flow systems see Appendix B, Chapter 33,
ASHRAE Handbook 1985 Fundamentals.
- For flow through an orifice (Chapter 2, ASHRAE 1985
Fundamentals):
3-5
-------
NBS: Indoor Air Quality Model PART I - Theoretical Basis
Phase II Report 3. Air Flow Analysis
C(j = orifice coefficient
* constant for turbulent flow (0.6 typically)
= (constant)/Re for laminar flow
D = diameter of approach to orifice
d = diameter of orifice opening
Thus the loss sum takes the form;
ZAP = (—) (c pv2 + c pv2 + c pv2 +...) (3.6)
o 2a ° ° ' P P q q
y
Recognizing that the mass flow rate, we, at each of these sections must be
equal;
= pVA = ... = pV A = pV A = pV A ='... = pV0A0 (3.7)
r \ \ roorpprqq r 2 2
equation (3.6) may be rewritten in terms of mass flow rate as;
= (1/2g p)(C /A2 + C /A2 + C /A 2 + ... )(we )2 (3.8)
o ^cr o o P p q q
and equation (3.4) then simplifies to;
2 - Ze) - C6(we) (3.9)
where;
C6 = (1/2g p)(-1/A2 + ... C /A2 + C /A2 + C /A 2 ... +1/A2) (3.10)
ycr 1 o o p p q q 2
Equation (3.9) may now be rewritten in terms of the element pressure DOFs,
using equation (3.1), as;
3-6
-------
NBS: Indoor Air Quality Model
Phase II Report
PART I - Theoretical Basis
3. Air Flow Analysis
It may be seen from equation (3.11) that mass flow through element e is driven
by the absolute pressure differences between zones (Pm - Pp) modified by
buoyancy effects created by density differences that are, in turn, due to zone
temperature differences.
Introducing a new variable, Be, for the buoyancy induced pressure component;
en
K - V
equation (3.11) may be rewritten as;
w
or
we =
- Pne) * aV
where; a = (C
Be )'
(3.133)
(3.13c)
where the second form, equations (3.13b) and (3.13c), will provide the correct
sign for we.
Variation of Flow With Zone Pressure
It is useful, at this point, to develop analytical expressions for the variation of
mass flow with zone pressure. This expressions will be seen to be useful for
solving the nonlinear flow system equations using schemes based upon the
classical Newton-Raphson iteration method. Therefore, from equations (3.13b)
and (3.13c) we obtain;
3-7
-------
NBS: Indoor Air Quality Model PART I - Theoretical Basis
Phase II Report 3. Air Flow Analysis
gpe aps
m m
( | R -Pe+Be | ) - l(cV( | Pe-Pe+Be | r (3.,4b)
9pe 2 9pe ' m n ' 2 ' m n '
n n
and from equation (3.10) we obtain;
9 3C _o 3C 9 3C
= (1/2g p)(A"2^ - A~2—£ + A"2—9- + ... ) (3.15a)
8Pe C ° 3Pe P 3Pe q 3Pe
m m m m
? 9C 9 3C _9 3C
r2-^- - A"2^ + A"2—S. - ... ) (3.15W
3Pe C ° 3Pe P 3Pe q 3Pe
n n n n
that is, the variation of Ce with pressure is simply a weighted sum of the
variation of individual pressure loss coefficients contributing to the total
pressure loss along the element. Analytical expressions for these partial
derivatives of the pressure loss coefficients are not easily formulated, but by
considering limiting cases of flow we can gain some insight.
In general, the loss coefficients depend, in a rather complex and poorly
understood way, upon the nature of flow, as indicated by the Reynolds number,
Re, and detailed characteristics of the flow geometry (e.g., roughness,
constrictions, etc.). For many situations, however, the loss coefficients are
practically constant for the limiting case of fully turbulent flow (i.e., Re > 106), at
one extreme, and proportional to 1/Re for laminar flow (i.e., Re < 2 x 103) at the
other;
C0 = constant (3.16)
for fully developed turbulent flow
C0 - C*0/Re = C*0 u/pD0V0 (3.17)
3-8
-------
NBS: Indoor Air Quality Model PART I - Theoretical Basis
Phase II Report 3. Air Row Analysis
fully developed laminar flow
where; C*0 = constant
- For conduits of constant cross-section;
= 64 L/Deq
- For "fittings" values of C*0 are not available; it may be
reasonable to estimate values based upon equivalent
lengths of conduits used in turbulent flow calculations
(e.g. see ASHRAE 1985 Handbook of Fundamentals
Chptr 34).
- For flow through an orifice;
= ??
|i = fluid viscosity
D0 = a characteristic dimension of the flow geometry
In fully developed turbulent flow, with each of the pressure loss coefficients
constant, the partial derivatives of equations (3.15) become zero and
consequently the first term of equations (3.14) becomes zero and, using
equations (3.13), may be simplified to;
3we 1 e
^-LL- = — a ; for fully turbulent flow (3
8P8 2
m
3we 1 e
X^— =- — a ; for fully turbulent flow (3.18b)
8P6 2
n
Limiting consideration to flow resistance elements of constant cross-section, we
may formulate a modified expression for laminar flow in an element, in a
manner similar to that used to formulate equations (3.13). We obtain;
we ~ a,e (Pe - Pe) + a,e Be (3.i9a)
L m n L
K K X
c c c
where; a. = (2g p/ji)( —2_ + —2_ + —2_ + ... ) (3.19b)
L 3cr ^ D A DA DA
oo p p q q
3-9
-------
NBS: Indoor Air Quality Model
Phase II Report
PART I - Theoretical Basis
3. Air Row Analysis
for which the evaluation of the variation of flow with pressure is straightforward;
8w _
8P
8w
2
m
e
= a.
= - a.
laminar flow, constant cross section
laminar flow, constant cross section
(3.20a)
(3.20b)
It is instructive to compare the fully turbulent flow equation, equation (3.13) with
Ce constant, with this particular case (i.e., constant cross section) fully laminar
flow equation;
Fully Laminar Flow
i • • •. J^L
L m n
m n
m n
n
Fig. 3.4 Limiting Case Flow Relations- Elements of Constant Cross-Section
It is seen that ae, the tangent slope of the fully turbulent curve, becomes
unbounded as flow approaches zero-flow conditions while a® does not.
If the variations of the pressure loss coefficients, C0 , Cp , Cq with flow are
well defined (i.e., for conduits: if the friction factor relations are reliable) then the
flow defined by equations (3.13) should asymptotically approach these two
curves at the upper and lower limits of flow. (Note: this is not to say that these
3-10
-------
NBS: Indoor Air Quality Model PART I - Theoretical Basis
Phase II Report 3. Air Flow Analysis
two curves provide an upper or lower bound to flow magnitude, in fact, they do
not (e.g., orifice flow: see reference [22] Fig. 18)).
Our purpose, here, is not to use these limiting-case flow relations in place of the
more general relation of equations (3.13), but rather to use these limiting cases
to provide an estimate of the variation of element flow with zone pressure to be
used in nonlinear solution algorithms. Specifically, we shall only employ
equations (3.19) and (3.20) for very low flow conditions, when the more general
expression for flow, equation (3.13b), and the approximation for the variation of
flow with pressure, equations (3.18), will tend to become unbounded.
Matrix Formulation of the Element Flo'w Equations
The element equations may be recast into matrix form, using the element DOFs
defined above, by first noting;
W = W = -W • (3.21)'
m n
thus;
(3.22a)
where;
} = [vv W } (3.22b)
net m ' n
= the element net mass flow rate vector
>e} = {P2 P6} (3.22c)
m ' n
= the element pressure vector
3-11
-------
NBS: Indoor Air Quality Model PART I - Theoretical Basis
Phase II Report 3. Air Row Analysis
r i n
P PI i ^ i i
[a ] = a i i I ; for all but very low flow conditions (3.22d)
= a _i i I ; for very low flow conditions (3.22e)
= matrix of pressure-flow coefficients
{w_} = a B {1-1} ; for all but very low flow conditions (3.220
D
= a B { 1' -1 } ; for very low flow conditions (3.22g)
= bouyancy-induced mass flow rate vector
and;
}
net - u ' ' ' ' : for all but very low flow conditions (3.23a)
e r 1 r
- !_ ] ~]
"2 L -i i.
-ae[ ]-jl
•\ L -i u
; for very low flow conditions (3.23b)
The element pressure-flow coefficients ae and a® are defined in such a way that
they are always positive, therefore, the matrix of pressure-flow coefficients will
be positive semi-definite.
Some complicating details deserve special note;
a) the direction of flow will be determined by the sign of (pm - pn + BS); if
positive, the flow will be from m to n,
b) the density p, of the fluid flowing through the element, will depend on the
direction of flow;
3-12
-------
NBS: Indoor Air Quality Model PART I - Theoretical Basis
Phase II Report 3. Air Flow Analysis
= P ; for flow from m to n
rm
= p ; for flow from n to m
rn
c) the flow coefficient, C®, will also depend on the direction of flow due to the
dependency of p on direction and the dependency of the pressure loss
coefficients C0 that also, in general, depend on the direction of flow,
d) the pressure-flow coefficient matrix [aej will also be flow-direction
dependent due to the flow-direction dependency of Ce and Be,
e) equation (3.22a) is highly nonlinear due to the flow-direction
dependencies, noted above, the dependency of the pressure-flow
coefficient matrix [ae] and the buoyancy-induced mass flow rate vector {WQ}
on the pressure, and the dependency of density on fluid temperatures which
are, in turn, dependent on the rate of flow.
3-13
-------
NBS: Indoor Air Quality Model
Phase II Report
PART I - Theoretical Basis
3. Air Row Analysis
3.2.2 Fan/Pump Element Equations
General operating characteristics of fans are discussed in the ASHRAE
Handbook and Product Directory: 1979 Equipment [23]. Flow behavior of fans
is generally described in terms of performance curves that have the following
typical form;
shutoff; ;
/•........!....;.....: {....}
Volumetric Flow Rate = AV
Actual Performance Curve
Mass Flow Rate = pAV
Performance Curve Represented
as a Familu of Linear Functions
W
o
Fia. 3.5 Fan Performance Curves
Performance curves may be easily converted to pressure-mass flow curves, by
scaling by the fluid density, and represented by the family of equations of the
general form;
e e e
w = w - a AP
o
where; w
(3.24)
= the free delivery mass flow rate of the fan
w
e e/ eN
a = a -(w J ; the fan pressure-flow coefficient
AP = the effective pressure drop across the fan .
This family of equations has the advantage of being able to represent saddle
shaped performance curves but, unfortunately, the members of the family used
3-14
-------
NBS: Indoor Air Quality Model PART I - Theoretical Basis
Phase II Report 3. Air Row Analysis
to "capture" the saddle shape portion of the performance curve provide very
poor representations of the change of mass flow with changes of pressure and,
therefore, should not be expected to perform well when used with Newton-
Raphson type nonlinear solution strategies.
For nonlinear solution techniques that require the determination of the change
of mass flow with changes of pressure we shall have to resort to a more
restricted form of representation having;
ae = ae(AP) (3.25)
Unfortunately, a true saddle shape may not be represented with this form.
An attractive candidate for this more restricted form is offered by the following
polynomial form;
ae = a^ + a^AP + a* AP2 + ... (3.26)
or
we = w6 -(aVp + a^AP2 + a* AP3 + ...) . (3.27)
o V 2 3
where the coefficients, a®, a| would be determined by a best fit to published
or measured performance curve data.
Defining fan element degrees of freedom consistent with flow resistant element
degrees of freedom, as shown below, Fig. 3.6, and accounting for buoyancy
effects, as in the development of the flow resistant element equations, equation
(3.27) may be rewritten as;
we = we -ae(Pe-Pe + BB) (3.28)
o m n
or in terms of element flow rate DOFs as;
= [ae]{Pe} + {w°} + {w }
u O
(3.29a)
3-15
-------
NBS: Indoor Air Quality Model
Phase II Report
PART I-Theoretical Basis
3. Air Flow Analysis
AP
w
m
m
m
m
Zone
w
,
n n
ZoneJ
DATUM
Fig. 3.6 Fan Element DOFs
where, now;
[ae] =
'[ •'. -i]
= aeBe( -i
i }
{w ) = we{ I -I)
o o
(3.29D)
(3.29c)
(3.29d)
Typically, the fan pressure-flow coefficient will be positive and therefore the
matrix of fan pressure-flow coefficients, [ae], will be negative semi-definite. To
account for the possibility of a system driving a fan beyond the shut off pressure
- free delivery range (i.e., to account for the possibility of back flow or pressure
assisted forward flow) the fan performance curve must be defined outside the
conventional range of flows.
Using the polynomial form of fan performance curve, equation (3.27), we may
develop analytical expressions for the variation of flow with zone pressures;
3-16
-------
NBS: Indoor Air Quality Model PART I - Theoretical Basis
Phase II Report 3. Air Flow Analysis
e eeee eeee2
= - a, - 2a0(P - P + B ) - 3a,(P - P + B ) - ... (3.30a)
_,_R 1 2 m n 3 m n
3P
m
= * a" + 2aV- Pe + Be) + 3a'(Pe- Pe+ Be)2 + ... (3.30b)
,-^e 1 2 m n 3 m n
3P
n
or, in terms of the element mass flow rate DOFs;
p. ? r -1 11
(3.31)
r -i 11
3-17
-------
NBS: Indoor Air Quality Model PART I - Theoretical Basis
Phase II Report 3. Air Flow Analysis
3.3 System Equations
Requiring conservation of mass at each flow-related node we demand;
( mass generation^ _ V1 ( net mass flow "\
I I — / j I I
V rate / A~. Arate into element/
connected J system node
elements
the element equations may be assembled to form a system of equations, that
govern the flow behavior of the system, of the form;
(W) = [AKPJ + {WJ + (W }
D O
(3.33a)
where;
(w) = {w,, w2,... wn}T (3i33b)
(P) = {P, , P2 , ... Pn)T ' (3.33c)
NR NF
[A] = A[ae] + A[ae] (3.33d)
9=1 9=1
NU NE
{WJ = A (WQ) + A{wn} (3.33e)
B e=i B e=1 B
NF
(W } = A{W6} (3.330
0 e=1 °
NR , NF = the number of flow resistance and fan elements respectively
A = the element assembly operator; a combination Boolean
transformation and matrix summation (see section 2.2, [2] or [24]
for details)
The system flow matrix, [A], is" the sum of positive semi-definite flow resistance
element matrices and negative semi-definite fan/pump element equations
3-18
-------
NBS: Indoor Air Quality Model PART I - Theoretical Basis
Phase II Report 3. Air Flow Analysis
and, therefore, may become negative definite! Given the "1 , -1 ,1 , -1" form of
the flow resistance element equations and the "-1 ,1 , -1 ,1" form of the
fan/pump element equations we need only check the diagonal elements of the
[A] matrix - if any are negative then [A] will be negative semi-definite otherwise
it will be positive semi-definite. As will be seen in the following examples,
transformation from a semi-definite to a definite matrix results upon the
specification of a single nodal pressure.
3-19
-------
NBS: Indoor Air Quality Model
Phase II Report
PART I - Theoretical Basis
3. Air Flow Analysis
3.4 Simple Examples
Two two-zone air flow examples are considered below. For these examples
the density of air will be estimated using the ideal gas law as;
p = (M/RXP/T) = ( 28.9645 kq/kqmolej(p/T) = 0.00348365 (P/T)
K 8314.41 N-m/kqmole-°K
where;
-m/kgmole
p = density [=] kg/m3
M = the mean molecular weight per mole of dry air
R = the universal gas constant
P = the absolute pressure [=] Pa (i.e., N/m2)
T = the absolute temperature [=] °K
Example 1
In the first example, illustrated below, two zones are linked by two flow
resistance elements, conduits in this case.
1 m
1-4—M D = 0.25 m
4
T,,P,
1 1
t
2m
LI*
k
3.5m
X
r
elema
O* O
^
elemb
. A
4
2m
0.5 m
Fig. 3.7 Examle 1 Flow Idealization
The temperature in zone 1 is maintained at 10 °C and that of zone 2 at 20 °C
or;
^ =(10 + 273.15) = 283.15 °K
273.15) = 293.15°K
3-20
-------
NBS: Indoor Air Quality Model PART I - Theoretical Basis
Phase II Report 3. Air Flow Analysis
and we seek to determine the mass flow rates through these elements and the
zone pressures that will be induced by buoyancy-driven flow induced by these
zone temperature differences.
Zone Densities:
- assume sea level pressure 101.325 kPa
- p10oc = 0.00348365 * (101.3257(10° + 273.15°) ) = 0.0012466 kg/m3
- p20°c = 0.00348365 * (101.3257(20° + 273.15°) ) = 0.0012041 kg/m3
Element Equations:
- Relative roughness = e/D = 0.00015/0.25 = 0.0006
- Friction factor: from ASHRAE Fundamentals, Chapter 2, Fig. 13; f = 0.032
- Cross sectional area: A = TTD2/4 = rr 0.252/4 = 0.049 m2
- Pressure loss coefficient: Co = f L/D = 0.032 * 1.0 - 0.25 = 0.128
- Element a : connectivity 1-2
- Assume flow is from zone 2 to zone 1 thus p = p20°C
Ca = (1/2gcp)( C0/A20) = (1/(2 * 1 * 0.0012041)(0.128/0.0492) = 22137
= (9.81/1.0X0.0012466(2-3.5) +0.0012041(3.5-3.5) +
0.0012041(3.5-2))
= -0.00062576
- Initial element matrices (from equations (16) ): (assume Pam = pan)
(1/Ca | pam- pan+ ea J )'/2 = (i/(22137 * 1 0 + (-0.00062576) | )'/2 =
0.268679
f} = { -0.00016813 0.00016813}
D
- T
~ L
0.268679 -0.268679]
-0.268679 0.268679 J
Element b: connectivity 1-2
- Assume flow is from zone 1 to zone 2 thus p = pio°c
3-21
-------
NBS: Indoor Air Quality Model PART I - Theoretical Basis
Phase II Report 3. Air Flow Analysis
C& = (1/2gcp)( C0/A20) = (1/(2 * 1 x 0.0012466))(0.128/0.0492) = 21382
Bb
= (9.81/1.0)(0.0012466(2-0.5) +0.0012466(0.5-0.5)
+0.0012041(0.5-2))
= 0.00062576
- Initial element matrices (from equations (16) ): (assume Pbm = pbn)
(1/Cb | pbm- pbn+ Bb | )'/2 = (17(257923 x 1 0 + (0.00062576) | )'/2 = 0.27338
j = { 0.00017107 -0.00017107 }T
D
- [°-:
• L-o.
.27338262 -0.27338262
27338262 0.27338262 J
System Equations:
The system equations may be assembled from the element equations; in this
case we obtain, assuming no mass generation in the zones;
/ 0 | _ f 0.54206194 -0.54206194] f P\\ f 0. 00000294
\. Oj ~ L -0.54206194 0.54206194 J \.P2J + I -0.00000294
J
As they stand this set of equations is singular - they describe only the pressure
difference between zones. If we specify the pressure in one zone, say P! =
101.325, then a first estimate of P2 may be determined; P2 = 101.32500543.
The element arrays may then be recomputed with these new estimates of P-\ &
P2 and the system equations formed and solved. By repeating this process
until the results converge to acceptable accuracy a solution is obtained. For
this problem we obtain, upon convergence;
PI = 101.3250000 Pa (i.e., as specified)
P2 = 101.3250814 Pa
3-22
-------
NBS: Indoor Air Quality Model
Phase II Report
PART I-Theoretical Basis
3. Air Row Analysis
Wa = -0.00016922 kg/sec
Wb = 0.00016995 kg/sec
For comparison, the system equations at convergence are;
_ f 0.54216990 -0.54216990~| / P, \ f 0.00000515
~ L-0.54216990 0.54216990 \ \ P2J 1-0.00000515
Example 2
In this example, illustrated below, two zones are linked by a flow resistance
element, identical to element "a" used in the example above, and a fan
element.
1 m
UN D = 0.25m
Trp,
i
3.5m
2m
i
elema
t
elemb
~"tv/vt~ A
T2,P2
2m
0 0.0015
Fan F'erf ormance Curve
^b = Wo ~ cbAP
= 0.0015 - 0.00001 AP
w [=] kg/sec
0.5m
Fig. 3.8 Example 2 Flow Idealization
Again the temperature in zone 1 is maintained at 10 °C and that of zone 2 at
20 °C and we seek to determine the mass flow rates through the elements and
the zone pressures that will be induced by the combined effects of buoyancy-
driven and fan-driven flow.
Element Equations:
- Element a: connectivity 1-2 (as above)
- Initial element matrices (from example 1): (assume pam = pan)
{W0} = { -0.00016813 0.00016813}
D
3-23
-------
NBS: Indoor Air Quality Model PART I - Theoretical Basis
Phase II Report 3. Air Row Analysis
- r °-:
" L -o.
.268679 -0.268679"]
268679 0.268679 J
- Element b: connectivity 1-2
- From the fan performance curve, above, we obtain Ct>= 0.00001 and wt»o =
0.0015
- B& is equal to that calculated for the resistance element b above; B& =
0.00062576
- Initial Element Matrices (from equations (18) > (assume pam = pan)
{W^} = { -0.00000001 0.00000001 }
B
{Wb} = { 0.0015 -0.0015 }
o
r &1 _ [ -0.00001 O.OOOOll
La J " L 0.00001 -0.00001 J
System Equations:
The system equations may be assembled from the element equations; in this
case we obtain, assuming no mass generation in the zones;
0\ _ f 0.26866932 -0.268669321 /P, \ ^ / -0.00016814\ ^ / 0.0015\
OJ ~ L -0.26866932 0.26866932 J \P2J I 0.00016814J \-0.0015j
Again we obtain a singular set of equations that describe the pressure
difference between zones. By specifying one zone pressure, say Pi =
101.325, a first estimate of P2 may be determined, in this iteration P2 =
101.32995727. Again, we iteratively update element matrices with these
estimates of zone pressures until results converge to acceptable accuracy.
Reasonably convergent results are;
P, = 101.32500000 Pa (i.e., as specified)
P2 = 101.37379028 Pa
3-24
-------
NBS: Indoor Air Quality Model PART I - Theoretical Basis
Phase II Report a Air Row Analysis
wa = -0.00149407 kg/sec
wt> = 0.00150048 kg/sec
For comparison, the system equations at "convergence" are;
( 0\ - F 0.03022450 -0.03022450 1 . / -0.00001893 \ . / 0.0015
I Oj ~ L-0.03022450 0.03022450 J I 0.00001893J X -O.OOIS
3-25
-------
NBS: Indoor Air Quality Model PART I - Theoretical Basis
Phase II Report 3. Air Row Analysis
3.5 Solution of Flow Equations
Two classic nonlinear solution strategies and their variations;
a) Method of Successive Substitutions or Fixed-Point Iteration
Direct
Jacobi Iteration
Zeid's Modified Jacobi Iteration
Gauss-Seidel Iteration
Successive Overrelaxation Method
b) Newton-Raphson Method
Classic Newton-Raphson Method
Modified Newton-Raphson Method
and incremental formulations of these methods will be considered as
candidates for solving the system of nonlinear flow equations, equations (3.33).
To set the stage for a discussion of these solution methods we rewrite the
system equations, equations (3.33), in two alternate forms:
(F(P)} = [AKP} + {WJ + {W } - {W} = {0} (3.35)
D O
and
[AKP} = {g} = {W} - {W } - {WJ (3.36)
O D
where, it is important to be mindful that [A] and {Wg} are both dependent on the
state of the system pressure variables, {P}, and may also vary with time if the
flow prblem is embedded in a dynamic thermal response problem.
3.5.1 Successive Substitution
A class of nonlinear solution techniques have been developed and studied for
equations of the form of equation (3.36) with {a} not a function of the dependent
variable {P} that are based upon the use of an approximate inverse [C]. By
3-26
-------
or
NBS: Indoor Air Quality Model PART I - Theoretical Basis
Phase II Report 3. Air Flow Analysis
adding the vector [C]{P} to both sides of equation (5.2);
[AKP} + [CKP} = {g} + [C](P} (3.37a)
(P) = (P) + [C]H{ {g} - [A]{P} } (3.37b)
the governing equation is recast in a form that suggest a general iterative
scheme;
(Pk+1) = (Pk) - [cVi {gk} - [Ak]{Pk}} (3.38)
where k is an iteration index. The term {{gk} - [Ak]{pk}} may be thought of as a
residual or error that could be monitored to evluate the convergence of the
method.
The choice of the [C] matrix is key to the success of this approach. Clearly [C]
must be nonsingular. Zeid shows, furthermore, that to ensure convergence [C]
must satisfy the following condition [25],[26];
[I]-[C]"'[A] < 1 (3.39)
where the double bars || indicate any appropriate norm (e.g., maximum norm or
Euclidean norm).
We shall consider the following alternatives, based on those developed for
systems with {g} not a function of the dependent variable {P};
Direct Iteration
The most straigtforward approach simply sets [C] = [A];
3-27
-------
NBS: Indoor Air Quality Model PART I - Theoretical Basis
Phase II Report 3. Air Flow Analysis
[Akr1{ {gk> - [
- or
{Pk+1} = [Akr1{gk>
Computationally, it is efficient to avoid inversion and instead successively solve
the system of equations;
+ (3.40C)
For systems with {g} # (g(P)} this method often does not converge [27] and,
therefore, will not be considered further.
Jacobi Iteration
Splitting the [A] matrix into upper and lower components as;
[A] = [D][ [L] + [I] + [U] ] ; [D] = diag (A^) (3.41)
we set [C] = [D] to obtain;
{pk+i} = {pk} + [Dk-,-i{ {gk}_[Ak](pk} } (342a)
or
{Pk+1} = [Dk]~1{gk} - [ [Lk] + [Uk] ]{Pk} (3.42b)
For systems with {g} * (g(P)} this method converges if [Ak] is strictly diagonally
dominant [25],[26]. In general, [A] will not be strictly diagonally dominant, thus.
this method is not useful here.
Zeid's Modified Jacobi Iteration
Zeid has developed a modified form of Jacobi iteration that does not require
strict diagonal dominance [25],[26]. In this method we set;
3-28
-------
NBS: Indoor Air Quality Model PART I - Theoretical Basis
Phase II Report 3. Air Flow Analysis
[C ] = diag(oO ; - [Uk]{Pk} + [Dk]~1{gk} (3.44b)
For systems with {g} * (g(P)} the rate of convergence of this method is linear. In
indicial notation this method is;
k\f \f
A—. N N
.v. - z^ A.P. + g.
K j = 1 >J J H 'J J '
r = -" -1 (3.44c)
1 Ak.
n
P.+ = P. + r. ; i = 1, 2, ... n (3.44d)
where r is the residual that may conveniently be monitored to evaluate
convergence.
Successive Overrelaxation Method
A variant of of Gauss-Seidel iteration, commonly know as the successive
overrelaxation or SOR method, attempts to to accellerate convergence by
scaling the residual by a relaxation factor, o>, as;
3-29
-------
NBS: Indoor Air Quality Model PART I - Theoretical Basis
Phase II Report 3. Air Row Analysis
{pk+l} = {pk} + t| + Lkj-1^-1^ {gk} _ [Ak]{pk} }
or
{Pk+1} =-[l_k]{Pk+1} + (1-co){Pk} + [[Lk] + co[Lk]]{Pk}
-co[Uk]{Pk} + co[Dkr1{gk} (3.45b)
where for co=1.0 this reduces to Gauss-Seidel iteration. In indicial notation this
method is;
k _
i
\V |< |< +
-|,AUPJ
' " 2Ajjpj *
Ak.
k
,
(3.45c)
_k+1 nk k
P. = P. + ur. I = 1, 2, ... n (3i45d)
This method can only converge for 0 < co< 2 [28].
For the governing flow equations, equations (3.36), the forcing vector {g} will, in
general, depend upon the dependent variable {P} and thus the convergence
rates and conditions on convergence noted above can, at best, provide only
guidelines; we are not in a position at this time to say much about the
convergence of these adaptations of classical fixed-point methods.
Upon closer examination, however, we note that {g} = (WB(P)} + {W0}, the sum
of a bouyancy-related flow vector, that is pressure dependentand a fan-related
flow vector that is not. If the flow is largely forced (i.e., by fans or wind-induced
pressure), so that the bouyancy-related flow is relatively small, then we should
expect these adapted methods to behave as theory predicts.
3.5.2 Newton-Raphson Iteration
The following development of the Newton-Raphson Method and its variants is
based largely on the formulation presented by Bjork and Anderson [28].
3-30
-------
NBS: Indoor Air Quality Model PART I - Theoretical Basis
Phase II Report 3. Air Row Analysis
Using Taylor's formula, generalized for a system of n equations, we may
approximate the function {F(P)}( from equation (3.35), from its value at a nearby
vector {Pk} as;
} + [F'(PN)]{ {P} - {PN} } + 0( {PHP*} \ ') (3.46)
where F' is the Jacobian defined as;
3{P) '{P}={Pk}
or
k
F|H(P)="~ k (3'47b)
IJ 9Pj {P)=(Pk}
Equation (3.46) leads naturally to the general form of the popular Newton-
Raphson iterative method;
[F'(Pk)]{APk+1} = - (F(Pk)} (3.48a)
} (3.48b)
where, again, k is the iteration index. Given an initial guess {Po} sufficiently
close to the solution the method will converge at a quadratic rate.
The high rate of convergence has made this approach popular, but the method
involves the formation of the n x n entries of the Jacobian and the solution of an
n x n system of equations at each iteration - tasks that become computationally
prohibitive as n increases.
Evaluation of the Jacobian
For the problem at hand, equation (3.35), the Jacobian involves the evaluation
of;
3-31
-------
NBS: Indoor Air Quality Model PART I - Theoretical Basis
Phase II Report 3. Air Flow Analysis
{W_}}
— I ,. (3.49a)
3{P)
or
e
{p}={P} e=1 8{P} {p(P)} 2=} 3{pe} {P(p)}
e
G49W
that is, the Jacobian is simply evaluated as an element assembly sum of the
element Jacobians evaluated at the element pressures {Re} corresponding to
the current iterative estimate of the system pressure {Pk}.
Modified Newton-Raphson Iteration
To avoid some of the computational expense of forming and solving the
Newton-Raphson equations, (3.48), at each iteration, one may reform [A] , [F*]f
and {WB} only occaissonally, say every m steps, as;
[F'(PP)KAPk+1} = -[AP]{Pk} - {WP} - (W } - {W} ; k=p, ... p*i (3.50)
D O
This modified Newton-Raphson method saves computation but at a cost of
convergence. Attempts have been made to compensate for a lower
convergence rate by using an overtaxation factor , co, applied to the residual,
AP, calculated at each step, as;
(3.51)
the choice of co is likely to be "problem dependent and the experience of the
analyst will be crucial" [29]; values of co = 2.0 are often used.
3-32
-------
NBS: Indoor Air Quality Model PART I - Theoretical Basis
Phase II Report 3. Air Row Analysis
3.5.3 Incremental Formulation
When flow is driven primarily by fans (i.e., when the buoyancy-related flow is
relatively small) it may prove useful to approach a solution incrementally by
considering incremental increases in fan free delivery flow {W0}. After
Zienkkiewicz [27] we rewrite equation (3.35) in the form;
[A]{P} + {WJ + X{ {W } - {W} } = {0} (3.52)
D o
and solve a series of nonlinear problems, incrementally increasing \ to 1.0; the
solution at each increment may then be used as the initial guess of the solution
at the next increment. For increments of X suitably small we may be assured
that the initial guesses of the incremental solutions will be sufficiently close to
the solution to guarantee convergence, if the solution to;
[A]{P} + {W J = {0} (3.53)
D
is available (e.g., if {WB} is a zero vector) or can be computed.
For m increments of A,;
X = 1/m , 2/m ,... m/m (3.54)
the Gauss-Seidel method, with overrelaxation becomes, in indicial notation;
i-1 . , i n . ' .
v1 m . k m_ k+i v"1 m . k m_k m, ,k m. /, , , . x
- X A p ~ ,2^ A P - W + X(W - W )
•y' = -±LJLJ—^-^J — L_fl (3.55a)
V.
M
m k"*~ 1 m k m k
P. = P. + 00 r. ; i=l, 2,... n, (3.55b)
and the modified Newton-Raphson method, also with overrelaxation, becomes;
3-33
-------
NBS: Indoor Air Quality Model PART I - Theoretical Basis
Phase II Report 3. Air Flow Analysis
} - %{ (W } - {W} } (3.56a)
D O
{V*'> = {V> * wteV*'} U = p, p.l,... p.l (3.560)
with updating of system arrays every 1+1 steps. In both cases, at each
increment, m, one iterates on k.
3-34
-------
NBS: Indoor Air Quality Model PART I - Theoretical Basis
Phase II Report 4. Summary and Directions of Future Work
4. Summary and Directions of Future Work
Summary
The theoretical basis of a building indoor air quality model has been presented
that provides for;
a) contaminant dispersal analysis of nonreactive contaminants, and
b) mechanical, wind, and thermally-driven air flow analysis
in multi-zone buildings of arbitrary complexity. It has been shown that both
contaminant dispersal analysis and air flow analysis equations may be
assembled from element equations that govern the behavior of discrete flow
elements in the building airflow system. The general, qualitative character of
these equations has been discussed and efficient numerical methods have
been presented for their solution.
This theoretical work extends the work of others (e.g., [18], [30],[31]) in that;
a) for both contaminant dispersal and flow analysis;
- the governing equations are assembled from element equations so that
systems of arbitrary complexity may be considered, existing
computational strategies based upon element assembly methods may be
employed, and formal analysis of the system equations is possible from
the new perspective of the element assembly operation,
- efficient numerical methods have been identified for the practical solution
of the governing equations, and
b) for contaminant dispersal analysis;
- filtering of contaminants has been accounted for,
- practical methods of accounting for unsteady flow conditions have been
identified,
4-1
-------
NBS: Indoor Air Quality Model PART I - Theoretical Basis
Phase II Report 4. Summary and Directions of Future Work
- the qualitative analysis of the multi-zone contaminant dispersal equations
has been extended demonstrating, importantly, that the conservation of
total air flow, alone, in a building idealization (without the need to place
special qualifications on zones isolated from exterior air infiltration, e.g.,
[31] p. 225) leads to nonsingular M-matrices that may be efficiently
factored to LU form, and
c) for flow analysis;
- element equations governing passive resistance air flow paths has been
extended to allow consideration of a variety of simple and complex air
flow paths,
- element equations governing fan-driven air flow have been developed
that may readily be assembled, with the general resistance element, to
allow analysis of building air flow systems of arbitrary complexity, and
- low-flow conditions have been modeled consistently with existing flow
theory in such a way that should help to avoid convergence problems
experienced by others (some preliminary computational studies indicate
success here).
In PART II of this report a program, CONTAM86, is presented that implements
the contaminant dispersal portion of the theory and examples of its application,
that provide preliminary validation, are discussed.
Directions of Future Work
In the near future, work will be directed toward the two general areas
considered thus far - contaminant dispersal analysis and air flow analysis. In
addition, the inverse contaminant dispersal problem will be considered (i.e., the
determination of airflows, in a multi-zone building system, from knowledge of
zonal concentrations due to known excitations). In the distant future, hopefully,
the coupled multi-zone building flow and thermal analysis problem and its
integration with the contaminant dispersal analysis problem will be considered
by integrating the building thermal analysis methods developed earlier [2] with
4-2
-------
NBS: Indoor Air Quality Model PART I • Theoretical Basis
Phase II Report • 4. Summary and Directions of Future Work
the methods introduced here.
In the area of contaminant dispersal analysis the present "theory will be
extended;
a) through the development of reaction elements , to allow modeling of the
dispersal of single and multiple reactive contaminants, and
b) through the development of one-dimensional convection-diffusion flow
elements , to allow modeling of the details of contaminant dispersal for flow
in duct-type flow passages.
In addition, an attempt will be made to develop elements to model the dynamics
of contaminant adsorption and absorption into the building fabric and
furnishings.
The flow analysis theory will be implemented to provide computational tools that
may be used in an integrated manner with the contaminant dispersal analysis
tools presently available in CONTAM86. An attempt will be made to evaluate
the several nonlinear solution strategies, discussed in section 3.5, so that
guidelines for their use may be formulated.
The inverse problem of determining multi-zone air flow rates from measured
contaminant concentration and generation rate data (e.g., as used in tracer gas
flow measuring techniques) will, also, be addressed. That the inverse problem
is inherently an ill-conditioned problem (i.e., small errors in concentration and
generation rate data typically result in large errors in estimated airflow
quantities) is not well appreciated, therefore, this effort will place an emphasis
on determination of the conditioning of the inverse problem, for specific
applications, and identification of strategies of formulating the inverse problem
to minimize ill-conditioning. Coupling the formulation and solution of specific
inverse analysis problems with the determination of their conditioning provides,
as an additional benefit, a means to place error bounds on the estimates of
airflows. Again, the inverse problem will be formulated using an element
assembly approach, to allow consideration of systems of arbitrary complexity,
and implemented so as to augment the computational tools available and
presently under development for dispersal and flow analysis.
4-3
-------
NBS: Indoor Air Quality Model PART I - Theoretical Basis
PhaseIIReport PARTI References
PARTI References
[1] McNall, P., Walton, G., Silberstein, S., Axley, J., Ishiguro, K., Grot, R., &
Kusuda, T., Indoor Air Quality Modeling Phase I Report: Framework for
Development of General Models. NBSIR 85-3265, U. S. Dept. of Commerce,
National Bureau of Standards, October 1986
[2] Axley, J.W., DTAM1: A Discrete Thermal Analysis Method for Building Energy
Simulation: Part I Linear Thermal Systems with DTAM1 Users' Manual.
(presently under review for publication by the National Bureau of Standards,
Building Environment Division, Center for Building Technology)
[3] Bird, R.B., Stewart, W.E., & Lightfoot, E.N., Transport Phenomena. John Wiley
& Sons, lnc.,N.Y.,1960
[4] Zienkiewicz, O.C. & Morgan, K., Finite Elements and Approximation. John
Wiley & Sons, N.Y., 1983, pp. 154-157
[5] Funderlic, R.E. & Plemmons, R.J., "LU Decomposition of M-Matrices by
Elimination Without Pivoting", Linear Algebra and Its Applications, Vol 41,
pp. 99-110, Elsevier, North Holland.1981
[6] Plemmons, R.J., Nonegative Matrices in the Mathematical Sciences. Chapter
6: M-Matrices, Academic Press, 1979, p. 137 and 147
[7] Ibid, p 156
[8] Funderlic, R.E., Neumann, M., & Plemmons, R.J., "LU Decomposition of
Generalized Diagonally Dominant Matrices", Numer. Math., Springer-Verlag,
Vol 40, p 57-69, (1982) Theorem 4
[9] Graybill, F.A., Matrices with Applications in Statistics. Wadsworth, Belmont
CA., 1983 Section 11.4
[10] Noble, B., Applied Linear Algebra. Prentice-Hall, N.J., 1969
[11] Strang, G., Linear Algebra and Its Application. Academic Press, N.Y., 1980
Ret I -1
-------
NBS: Indoor Air Quality Model PART I - Theoretical Basis
PhaseIIReport PARTI References
[12] Ibid, p 213
[13] Wilkinson, J.H. & Reinsch, C., Handbook for Automatic Computation:
Volume II: Linear Algebra. Springer-Verlag, 1971 - Part II: The Algebraic
Eigenvalue Problem
[14] ibid., Contribution 11/12: "Solution to the Eigenproblem by a Norm Reducing
Jacobi Type Method", by P.J. Eberlein & J. Boothroyd
[15] Taylor, Robert L, *HEAT*. A Finite Element Computer Program for Heat-
Conduction Analysis. Report 75-1, Prepared for: Civil Engineering
Laboratory, Naval Construction Battalion Center, Port Hueneme, California,
May 1975
[16] Huebner, K.H., & Thornton, E., The Finite Element Method for Engineers:
Second Edition, John Wiley & Sons, New York, 1982
[17] Hughes, T.J., "Analysis of Some Fully-Discrete Algorithms for the One-
Dimensional Heat Equation", International Journal for Numerical Methods in
Engineering, Vol. 21, John Wiley & Sons, 1985
[18] Walton, G., Thermal Analysis Research Program Reference Manual. NBSIR
83-2655, U.S. Dept. of Com., National Bureau of Standards, Gaithersburg,
MD, March 83
[19] Carnahan, B., Luther, H.A., Wilkes, J.O., Applied Numerical Methods. John
Wiley & Sons, 1969
[20] ASHRAE, ASHRAE Handbook: 1985 Fundamentals. Chapter 33 Duct
Design,, ASHRAE, Atlanta, GA, 1985
[21] ASHRAE, ASHRAE Handbook: 1985 Fundamentals. Chapter 14 Airflow
Around Buildings, ASHRAE, Atlanta, GA, 1985
[22] ASHRAE, ASHRAE Handbook: 1985 Fundamentals. Chapter 2 Fluid Flow,
ASHRAE, Atlanta, GA, 1985
[23] ASHRAE, ASHRAE Handbook & Product Directory: 1979 Equipment,
Chapter 3 Fans, ASHRAE, New York, 1979
Refl-2
-------
NBS: Indoor Air Quality Model PART I • Theoretical Basis
PhaseIIReport PARTI References
[24] Bathe, K.J., Finite Element Procedures in Engineering Analysis. Second
Edition. John Wiley & Sons, New York, 1982
[25] Zeid, I., "Fixed-Point Iteration to Nonlinear Finite Element Analysis. Part I:
Mathematical Theory and Background," International Journal for Numerical
Methods in Engineering, Vol.21, No.11, John Wiley & Sons, New York, 1985
[26] Zeid, I., "Fixed-Point Iteration to Nonlinear Finite Element Analysis/Part II:
Formulation and Implementation," International Journal for Numerical
Methods in Engineering, Vol.21, No.11, John Wiley & Sons, New York, 1985
[27] Zienkiewicz, O.C., The Finite Element Method. 3rd Edition. McGraw-Hill,
London, 1977, p. 452
[28] Bjork, A., Anderson, N., Numerical Methods. Prentice-Hall, Inc., New Jersey,
1974
[29] Chow, Y.K. & Kay, S., "On the Aitken Acceleration Method for Nonlinear
Problems", Computers & Structures, Vol. 19, No. 5/6,1984, pp. 757-761
[30] Sinden, F.W., "Multi-Chamber Theory of Air Infiltration", Building and
Environment, Vol. 13, Pergamon Press, Great Britain, 1978, pp. 21-28
[31] Sandberg, M, "The Multichamber Theory Reconsidered from the Viewpoint
of Air Quality Studies", Building and Environment, Vol. 19, No. 4, Pergamon
Press, Great Britain, 1984, pp. 221-233
Refl-3
-------
NBS: Indoor Air Quality Model
Phase II Report
PARTII - CONTAM86 Users Manual
5. General Instructions
5. General Instructions
The program CONTAM86 is a command processor5'1; it responds to
commands in the order that they are presented and processes data associated
with each command. Commands may be presented to the program
interactively, using keyboard and monitor, or through the use of command/data
input files; that is to say, it offers two modes of operation - interactive and batch
modes.
For most practical problems of contaminant dispersal analysis the batch mode
of operation will be preferred. For these problems, analysis involves three basic
steps;
Step 1: Idealization of the Building System and Excitation
rap
• — Zone Node
2-Node Number
Flow Element-
Element Flow-
Actual Building Air-Flow System Idealization
Fig. 5.1 Idealization of the Building System and Excitation
Idealization of the building flow system involves
a) discretization of the system as an assemblage of appropriate flow
elements connected at system nodes,
b) identification of boundary conditions, and
c) numbering of system nodes optimally (i.e., to minimize the bandwidth -
5'1 CONTAM86 is written in FORTRAN 77. The complete source code for the program may be
found in the attached appendix.
5-1
-------
NBS: Indoor Air Quality Model
Phase II Report
PARTII - CONTAM86 Users Manual
5. General Instructions
node number difference - of system equations).
The excitation (i.e., specified contaminant concentrations and generation rates)
may be modeled to be steady or defined in terms of arbitrary time histories. For
the latter case initial conditions of nodal contaminant concentration will have to
also be specified.
Step 2: Preparation of Command/Data Input File
Edit
Fig. 5.2 Preparation of Input Command/Data File
In the batch mode, the program reads ASCII text files of commands and
associated data, collected together in distinct data groups, that define the
building flow idealization and excitation. The command/data input file may be
prepared with any available ASCII text editing program and given a file name,
, specified by the user. The must, however, consist of 8
or less alphanumeric characters and can not include an extension (i.e.,
characters separated from the filename by a period,".").
5-2
-------
NBS: Indoor Air Quality Model
Phase II Report
PART1I - CONTAM86 Users Manual
5. General Instructions
Step 3: Execution of CONTAM86
Plot Output File
.PLT
Binaru Data Files
Command/data
Input File
CONTAM
00101
01100
00011
0010'
01100
00011
.WDT
.FEL
0010
01100
00011
. EOT
Results Output File s-
.OUT
Fig. 5.3 Execution of CONTAM86
CONTAM86 is then executed. Initially CONTAM86 will be in the interactive
mode. To enter the batch mode the command "SUBMIT F=" may be
used to "submit" the command/data input file to the program. The program will
then proceed to form element and system arrays and compute the solution to
the posed problem. CONTAM86 reads the ASCII command/data input file and
creates an ASCII (i.e., printable) output file .OUT. The results of an
analysis, .OUT, may be conveniently reviewed using an ASCII
editor and, from the editor, portions or all of the results may be printed out. Key
response results are also written to the ASCII file .PLT in a format
that may easily be transferred to some spreadsheet and plotting programs (i.e.,
data values within each line are separated by the tab character) for plotting or
subsequent processing.
File Summary
Depending upon the commands processed, CONTAM86 will also create a
variety of binary files for out-of-core storage needed for subsequent processing.
A summary of files read and created includes;
Files Read
an ASCII input file specified by the user that contains
5-3
-------
NBS: Indoor Air Quality Model
Phase II Report
PARTll - CONTAM86 Users Manual
5. General Instructions
commands and associated data
Files Created
.OUT a printable ASCII output file that contains analysis
results
.PLT
.FEL
an ASCII output file that contains key analysis results in
a form that may be transferred to spreadsheet and/or
plotting programs
a binary file used for out-of-core storage of flow element
data
.WDT a binary file used for out-of-core storage of element flow
time history data
.EDT a binary file used for out-of-core storage of excitation
time history data
In the interactive mode is set to the default value of "CONTAM86"
and commands are read from the keyboard. A help command, "HELP" or "H",
will produce a screen listing of all available commands.
5-4
-------
NBS: Indoor Air Quality Model PARTII • CONTAM86 Users Manual
Phase II Report 6. Comand Conventions
6. Command Conventions
Commands and their associated data (if any) may be single-line or multiple-line
command/data groups.
Single-Line Commands
Single line command/data groups begin with the command keyword and may
have any number of associated data items identified by data identifies of the
typical form;
COMMAND A=n1,n2,n3 B=n4 C=n5,n6 D=c1c2c3
where n1,n2,n3,... is numeric data and c1c2c3 is character data. In this
example the keyword COMMAND is the command keyword and the data
identifiers are A=, B=, C=, and D=.
Multiple-Line Commands
Multiple-line command/data groups are delimited by the command keyword and
the keyword END and may have any number of data subgroups terminated by
the symbol "<" within. They have the typical form of;
COMMAND A=n1,n2
n1 I=n2,n3,n4 B=n5 C=c1c2c3c4
n1 I=n2,n3,n4 B=n5 C=c1c2c3c4
n1 I=n2,n3,n4 B=n5 C=c1c2c3c4
<
n1,n2,n3 D=n4,n5,n6 E=n7 F=c1c2c3
n1,n2,n3 D=n4,n5,n6 E=n7 F=c1c2c3
n1,n2,n3 D=n4,n5,n6 E=n7 F=c1c2c3
<
C1C2c3c4c5c6
END
6-1
-------
NBS: Indoor Air Quality Model PARTII - CONTAM86 Users Manual
Phase II Report 6. Comand Conventions
Classes of Commands
Two general groups of commands are available, the "Intrinsic Commands" and
the "CONTAM86 Commands". The "Intrinsic Commands" are useful, primarily, in
the interactive mode allowing the user to examine system arrays generated by
the "CONTAM86 commands" and save them for further processing by the CAL-
80 command processor or other command processors based on the CALSAP
in-core management routines [1]. The "CONTAM86 Commands" provide
contaminant dispersal analysis operations.
Command/data Lines
Normally the line length (i.e., the number of character and spaces on a line) is
limited to 80. A backslash T at the end of information on any line will, however,
allow the next line to be interpreted as a continuation of the first line providing
an effective line length of 160.
Use of the symbol "<" within in any line indicates the end of information on that
line. Information entered to the right of this symbol is ignored by the program
and may, therefore, be used to annotate a command/data input file.
An asterisk "*" at the beginning of any line will cause the line to be echoed as a
comment on the console and to the output file. Lines marked in this way may,
then, be used to annotate the output file and help indicate the progress of
computation when using the batch mode of operation.
Data Identifiers
Data identifiers and their associated data may be placed in any order within
each line of the command/data group with the exception that the first line of a
command/data group must begin with the command keyword. In some
instances data may not be associated with a data identifier, such data must be
placed first in a line.
Data
6-2
-------
NBS: Indoor Air Quality Model PART1I - CONTAM86 Users Manual
Phase II Report 6. Comand Conventions
Decimal points are not required for real numeric data. Scientific notation of the
form nnE+nn or nn.nnE+nn (e.g., 5.79E-13) may be used. Simple arithmetic
expressions employing the conventional operators +, -, *, and / may be used.
The order of evaluation is sequential from left to right - unlike FORTRAN or other
programing languages where other "precedence" rules are used.
If fewer data values are supplied than required the missing data will assumed to
be zero, blank, or set to default values as appropriate.
6-3
-------
NBS: Indoor Air Quality Model
Phase II Report
PARTII • CONTAM86 Users Manual
7. Introductory Example
7. Introductory Example
For purposes of contaminant dispersal analysis the specific command/data
groups that need to be included in a command/data input file will depend upon
the details of the flow system idealization, the nature of the excitation, and the
type of analysis to be computed. A specific introductory example, should
however, provide some useful insight into the more general aspects of
contaminant dispersal analysis using CONTAM86
Consider the two-story residence with basement shown, in section, below. In
this residence interior air is circulated by a forced-air furnace and exterior air
infiltrates the house through leaks around the two first floor windows. The flow
system may be idealized using flow elements to model the ductwork, room-to-
room, and infiltration flow paths as shown below.
CO2
Generation
0.55
kg/hr
70 m3/hr
• — Zone Node
2 — Node Number
70 m3/hr
Flow Element —
Element Flow —
Actual Building
Air-Flow System Idealization
Fia 7.1 Hypothetical Residential Example
For this building idealization we shall consider the hypothetical problem of
determining the steady state distribution of CO;? generated by a kerosene
heater placed in room "2", distributed by the furnace flow system operated at
constant conditions, and diluted by infiltration at a constant rate. The CO2
7-1
-------
NBS: Indoor Air Quality Model PART1I • CONTAM86 Users Manual
Phase II Report 7. Introductory Example
generation rate is assumed to be 0.55 kg/hr, exterior C02 concentration is
assumed to be 760 u.g C02/ g air, and the assumed air volumetric flow rates are
indicated on the drawings above.
The CONTAM86 command/data file to complete this steady state analysis is
listed below. Command/data groups needed to complete a time constant
analysis and dynamic analysis for this building idealization are presented as
examples in the reference section of this manual.
7-2
-------
NBS: Indoor Air Quality Model
Phase II Report
PART1I - CONTAM86 Users Manual
7. Introductory Example
Command/data File for Residential Example
Note: CONTAM86 keywords and identifies are displayed in boldface below.
Description
Column
Comments:
Comments
Comments
Comments
Comments
Comments
System Definition:
Boundary Conditions
Flow Element Data:
Element Number & Connectivity,
Command/data File
1
Steady State Solution:
Flow Element Mass Flow Rates
*
* Six-Zone (7-Node)
Units: kg, m
Example
, hr
* Concentration [=] kg-CO2/kg-air
* Generation
*
FLOWSYS N=7
7 BC=C
END
FLOWELEM
1 1=1,2
2 1=1,3
3 1=7,2
4 1=2,7
5 1=7,3
6 1=3,7
7 1=2,4
8 1=3,5
9 1=4,6
10 1=5,6
11 1=6,1
END
STEADY
1,2 W=70*1.2
3,6 W=20*1.2
7,10 W=70*1.2
11 W-140*1.2
rate [=] kg-CO2/hr
< System has 7 Nodes
< Ext. "Zone" Cone. Spec.
< Flow Element 1
< Flow Element 2
< Flow Element 3
< Flow Element 4
< Flow Element 5
< Flow Element 6
< Flow Element 7
< Flow Element 8
< Flow Element 9
< Flow Element 10
< Flow Element 1 1
< (Air Density 1. 2 kg/m3)
< Supply Ducts
< Infiltration
< Return Loop
< Main Return Duct
Contaminant Excitation
Return to Interactive Mode
2 CG=0.55
7 CG=0.000760
END
RETURN
< Node 2: Generation Rate
< Node 7: Ext. CO2 Cone.
Details are given on the following pages for each of CONTAM86's
command/data groups.
7-3
-------
NBS: Indoor Air Quality Model PARTII - CONTAM86 Users Manual
Phase II Report 8. Command Reference
8. Command Reference
8.1 Intrinsic Commands
8.1.1 HELP
The command HELP, or simply H, will produce a list of all available commands,
in abbreviated form.
8.1.2 ECHO
The.command ECHO-ON acts to cause computed results normally directed to
the results output file to be echoed to the screen. The command ECHO-OFF
turns this feature off. At start-up CONTAM86 is set to ECHO-ON. Selective use
of ECHO-ON and ECHO-OFF can speed computation as writing results to the
screen consumes a significant amount of time.
8.1.3 LIST
The command LIST, or simply L, will produce a list of all arrays currently in the
in-core array database.
8.1.4 PRINT A=
The command PRINT A= or simply IP A= will "print"
array named , a one-to four character name, to the screen.
8.1.5 DIAGRAM A=
The command DIAGRAM A= will "print" a diagram of array named
, a one-to four character name, to the screen indicating position of
zero and nonzero terms. (Character arrays can not be diagramed.)
8.1.6 SUBMIT F=
The command SUBMIT F= will cause the program to switch to
batch mode and read all subsequent commands from the file .
8-1
-------
NBS: Indoor Air Quality Model PARTII - CONTAM86 Users Manual
Phase II Report 8. Command Reference
8.1.7 RETURN
The command RETURN returns the operation of the program from batch mode
to interactive mode. RETURN or QUIT will normally be the last line of batch
command/data input files.
8.1.8 QUIT
The command QUIT or simply Q terminates execution of the program and
returns the user to the control of the operating system.
8-2
-------
NBS: Indoor Air Quality Model PARTII • CONTAM86 Users Manual
Phase II Report 8. Command Reference
8.2 CONTAM86 Commands
The following conventions will be used for the command definitions presented
in this section;
- an ellipses,'...', indicates unlimited repetition of similar data items or data
lines within a data subgroup
- square brackets, [...], indicate optional data,
- numeric data is indicated by lower case n, as n1 ,n2,..., and
- character data by lower case c, as c1.
8.2.1 FLOWSYS
The size of the flow system and boundary conditions of system nodes are
defined with the following command/data group;
FLOWSYS N=n1
n2,n3,n4 BC=c1
END
where; ri1 = the number of flow nodes
n2,n3,n4 = first node, last node, node increment of a series of nodes
with identical boundary conditions
c1 = boundary condition code; C for concentration prescribed
nodes; G for generation prescribed nodes; (default = C)
The direct species mass generation rate or the species concentration - but not
both - may be specified at each node to establish boundary conditions of
prescribed contaminant generation or concentration.
If this boundary condition data is omitted all nodes will be assumed to be
species mass generation rate DOFs. Typically, nodes associated with outdoor
environmental conditions will be assigned specific contaminant concentrations
8-3
-------
NBS: Indoor Air Quality Model PART1I - CONTAM86 Users Manual
Phase II Report 8. Command Reference
and nodes associated with indoor air zones will be assigned specific species
generation rates although zero generation rates will often be appropriate for
these nodes.
See the introductory example presented earlier for an example of the use of this
command.
8.2.2 FLOWELEM
Two-node flow elements may be added to the flow system assemblage with the
following command/data group;
FLOWELEM
n1 I=n2,n3 GEN=n4 E=n5
• • •
END
where; n1 = the element number
n2, n3 = the element node numbers
n4 = generation increment (default = 1)
n5 = the element filter efficiency (default = 1.0)
Element data must be supplied in numerical order. Omitted data is automatically
generated by incrementing the preceding node numbers by the current
generation increment. Generated elements will have the properties of the
current element.
See the introductory example presented earlier for an example of the use of this
command.
8.2.3 STEADY
The response of the system to steady contaminant generation with steady
element mass flow may be computed with the following command/data group;
STEADY
n1,n2,n3 W=n4
8-4
-------
NBS: Indoor Air Quality Model PARTII - CONTAM86 Users Manual
Phase II Report . 8. Command Reference
n5,n6,n7 CG=n8
END
where; n1 ,n2,n3 = first element, last element, element number increment of a
series of elements with identical mass flow rates
n4 = element total mass flow rate; (default = 0.0)
n5,n6,n7 = first node, last node, node increment of a series of nodes
with identical excitation
n8 = contaminant concentration or contaminant generation
rate, as appropriate to the boundary condition of the node;
(default = 0.0)
Net total mass flow rate at each system node will be reported, but computation
will not be aborted if net mass flow is nonzero. The analyst must assume the
responsibility to check continuity of mass flow from these reported values.
See the introductory example presented earlier for an example of the use of this
command.
8.2.4 TIMECONS
System time constants, nominal and actual, may be computed with the following
command/data group;
TIMECONS [E=n1]
n2,n3,n4 W=n5
n6,n7,n8 V=n9
END
where; n1 = optional convergence parameter, epsilon ; (default =
machine precision)
n2,n3,n4 = first element, last element, element number increment of a
8-5
-------
NBS: Indoor Air Quality Model
Phase II Report
PARTII - CONTAM86 Users Manual
8. Command Reference
series of elements with identical mass flow rates
n5 = element total mass flow rate; (default = 0.0)
n6,n7,n8 = first node, last node, node increment of a series of nodes
with identical volumetric masses
n9 = nodal volumetric mass; (default = 0.0)
The nominal time constants are computed for each node as the quotient of the
nodal volumetric mass divided by the total air flow out of a zone. The actual
time constants are computed using an eigenanalysis routine that is a variant of
Jacobi iteration adapted for nonsymmetric matrices [2]. It should be noted that
the actual time constants are likely to be very different from the nominal time
constants for systems having well-coupled zones. Be advised: eigenanalysis of
the flow system matrices is a time consuming task.
Example
To determine the time constants associated with the building idealization
presented earlier, in the introductory example, the following command/data
group would have to be added to the command/data file.
TIMECONS
1,2 W=70*1.2
3,6 W=20*1.2
7,10 W=70*1.2
11 W=140*1.2
< (Air Density 1.2 kg/m3)
< Supply Ducts
< Infiltration
< Return Loop
< Main Return Duct
1 V=1.2*1.0
2.3 V«1.2*40.0
4.5 V-1.2*30.0
6 V-1.2*0.1
7 V=1.2*1.0E+06
END
< Node 1 Vol. Mass
< Nodes 2 & 3 Vol. Mass
< Nodes 4 & 5 Vol. Mass
< Node 6 Vol. Mass
< Node 7 Ext. Vol. Mass
8-6
-------
NBS: Indoor Air Quality Model PARTII - CONTAM86 Users Manual
Phase II Report 8. Command Reference
8.2.5 Dynamic Analysis
The response of the system, including transients, to general dynamic excitation,
may be computed using the command DYNAMIC. The dynamic solution
procedure used is driven by discrete time histories of excitation and element
mass flow data that must fjrsl be generated with the commands FLOWDAT and
EXCITDAT. (In future releases of CONTAM element mass flow data may also
be generated by a detailed flow analysis of the flow system.)
8.2.5.1 FLOWDAT
Discrete time histories of element mass flow rate may be defined, in step-wise
manner, from given element mass flow data, as illustrated below;
Data
Value
1 Given Data Value
-Time
TIME(I) TIME(2) TIME(3) TIME(4)
Fig. 8.1 Arbitrarily Defined Time History Data
px alternatively, discrete time histories of element mass flow data, defined in a
step-wise manner at equal time-step intervals along piece-wise linear
segments, may be generated from given element mass flow data over a time
range defined by an initial time, Tj, a final time, Tf, and a generation time
increment, AT, as illustrated below;
8-7
-------
NBS: Indoor Air Quality Model
Phase II Report
PARTII - CONTAM86 Users Manual
: 8. Command Reference
Data
Value
—Given Data Value
Generated Value
I ' I I ' '
AT
' ' ' I ' ' I
f
•Time
TIME(1) TIME(2) TIME(3) ' TIME(4)
Fia. 8.2 Eaual-Time-Step-GeneratQd Time History Data
using the following command/data group;
FLOWDAT [T=n1,n2,n3]
TIME=n4
n5,n6,n7 W=n8
TIME=n4
complete
n5,n6,n7 W=n8
[additional TIME data, as necessary, to define the
excitation time history]
END
where;
n1 ,n2,n3 = initial time, final time, time step increment used for
the piece-wise linear generation option
n4 = time value for subsequent data subgroups
n5,n6,n7 = first element, last element, element number increment of a
series of elements with identical mass flow data
n8 = prescribed element mass flow: (default = 0.0)
If data values n1.n2.n3 are specified, step-wise time histories will be generated
from the given data, along piece-wise linear segments as illustrated in Fig. 8.2
above, otherwise the given data will be used directly, as illustrated in Fig. 8.1
above.
At least two "TIME" data subgroups must be provided. FLOWDAT writes the
generated time history to the file .WDT so that this data may
subsequently be accessed by the command DYNAMIC.
8-8
-------
NBS: Indoor Air Quality Model
Phase II Report
PARTII - CONTAM86 Users Manual
8. Command Reference
8.2.5.2 EXCITDAT
Discrete time histories of excitation data may be defined in the two ways
discussed above for the FLOWDAT command using the following
command/data group;
EXCITDAT [T=n1(n2,n3]
T!ME=n4
n5,n6,n7 CG=n8
TIME=n4
n5,n6,n7 CG=n8
[additional TIME data, as necessary, to define the
complete excitation time history]
END
where; n1,n2,n3
n4
n5,n6,n7
n8,
= initial time, final time, time step increment used for
the piece-wise linear generation option
= time value for subsequent data subgroups
= first node, last node, node number increment of a series
of nodes with identical excitation data
= prescribed contaminant concentration or'prescribed
contaminant generation rate (as appropriate to node
boundary condition): (default = 0.0)
If data values n1,n2,n3 are specified, step-wise time histories will be generated,
from the given data, along piece-wise linear segments as illustrated in Fig. 8.2
above, otherwise the given data will be used directly, as illustrated in Fig. 8.1
above.
At least two "TIME" data subgroups must be provided. EXCITDAT writes the
generated time history to the file .EDT so that it may subsequently be
accessed by the command DYNAMIC.
8-9
-------
NBS: Indoor Air Quality Model PART1I - CONTAM86 Users Manual
Phase II Report 8. Command Reference
8.2.5.3 DYNAMIC
The response of the system to excitation defined by the EXCITDAT command,
using the prescribed element flow data defined by the FLOWDAT command,
may be computed using the following command/data group;
DYNAMIC
T=n1,n2,n3 [THETA=n4] [Pl=n5] [PS=n6]
n7,n8,n9 V=n10
n7,n8,n9 IC=n11
END
where; n1 ,n2,n3 = initial time, final time, time step increment
n4 = integration parameter, 9, where 0 < 0 <, 1 ; (default =
0.75) instability may result for 0 < 0.5,
n5 = response results print interval; (default = 1 )
n6 = plot file results scale factor; if not equal to 0.0, an ASCII
file, .PLT, of concentration response results will
be created with values scaled by the factor n6
n7,n8,n9 = first node, last node, node increment of a series of
nodes with identical data
n1 0 = nodal volumetric mass; (default = 0.0)
n1 1 = initial nodal concentration; (default = 0.0)
The response is computed using the predictor-corrector method discussed in
PART I of this report. With this method, the system flow matrix is updated at the
discrete times used to define element flow rate time histories and the system
excitation is updated at the discrete times used to define excitation time
histories, as illustrated below;
8-10
-------
NBS: Indoor Air Quality Model
Phase II Report
PART1I - CONTAM86 Users Manual
8. Command Reference
Flow
Data
Value
t
t
,
i
1'
,
Generated by
FLOWDAT
)
Flow Time Step
Excitation
Data
Value
<
*
1 ]
1
1
j
1
Generated by
EXCITDAT
Excitation Time Step
Respons
e
Value
Computed
by
DYNAMIC
Time
Integration Time Step
Fig. 8.3 Flow and Excitation Driven Dynamic Solution Procedure
The accuracy of the computed response is, therefore, dependent upon the
choice of the flow data time step, the excitation data time step, and the
integration time step chosen by the analyst. Furthermore, the flow data and
excitation data time steps may be nonconstant. The analyst should, therefore,
consider investigating the effects of the choice of these time step variables to
gain a sense of the error they induce.
8.2.5.4 Dynamic Analysis Example
To provide an example of a command/data sequence needed for dynamic
analysis we may consider an extension to the introductory example presented
earlier; the analysis of the dynamic response of the given building system,
under conditions of constant air flows, to a step change in CO2 generation.
Specifically, to consider the case where the kerosene heater is turned on and
then turned off 133 minutes later the following command/data group would have
8-11
-------
NBS: Indoor Air Quality Model PART1I • CONTAM86 Users Manual
Phase II Report 8. Command Reference
to be added to the command/data file used in the introductory example.
FLOWDAT
*
* Element flow rates modeled as constant.
*
TIME=0
1,2 W=70*1.2 < Supply Ducts
3,6 W=20*1.2 < Infiltration
7,10 W=70*1.2 < Return Loop
11 W=140*1.2 , < Main Return Duct
<
TIME=5
1,2 W=70*1.2 < Supply Ducts
3,6 W=20*1.2 < Infiltration
7,10 W=70*1.2 < Return Loop
11 W=140*1.2 < Main Return Duct
END
EXCITDAT < Nodal Excitation
TIME=0
*
* Kerosene heater turned on at time = 0 mins.
*
2 CG=0.55 < Node 2: Generation Rate
7 CG=0.000760 <.Node 7: Ext. C02Conc.
<
TIME=133/60
*
* Kerosene heater turned off at time = 133 mins.
*
2 OG-0.0 < Node 2: Generation Rate
7 CG=0.000760 < Node 7: Ext. CO2 Cone.
<
TIME=5
2 CG=0.0 < Node 2: Generation Rate
7 CG=0.000760 < Node 3: Ext. CO2 Cone.
END
DYNAMIC
1=0,4,0.1 PS=1 .OE+6 < Time-step; Plot Scale
1 V-1.2*1.0
-------
NBS: Indoor Air Quality Model PARTII - CONTAM86 Users Manual
Phase II Report 8. Command Reference
8.2.6 RESET
The command RESET resets the system in preparation for a new analysis
problem (i.e., key internal variables are re-initialized, contaminant dispersal
analysis system arrays are deleted from memory, and existing binary files are
deleted from disk storage). The system is automatically reset, if necessary,
upon execution of the FLOWSYS command.
RESET may be used to delete binary files that would otherwise be left on disk at
the termination of the program.
8-13
-------
NBS: Indoor Air Quality Model
Phase II Report
PART1I - CONTAM86 Users Manual
9. Example Problems
9. Example Problems
9.1 Single Zone Examples
It is useful to first consider a single zone building air flow system that exchanges
indoor air with the exterior environment. Such a single zone system may be
modeled as an assemblage of two flow elements, corresponding to inlet and
exhaust flow paths, connected to two system nodes, corresponding to the inside
air zone and the exterior environment "zone" as illustrated below;
C2 , G2, V2
Interior Zone
C,,G,,V,
. flow element
Element 1 number
Flow Rate w,
Element 2
Flow Rate w2
Fia. 6.1 A Single Zone Building and Corresponding Flow Model
The equations governing this simplest flow system have the following general
form;
-w2
w2 _
v, o"
0 N/2_
dC,
dt
dC2
dt
(9.1)
where;
Wi ,W2
Ci , C2
Vi,V2
G1,G2
= intake and exhaust element flow rates, respectively
= interior and exterior contaminant concentrations, respectively
= interior and exterior volumetric masses, respectively
= interior and exterior contaminant generation rates,
respectively.
9-1
-------
NBS: Indoor Air Quality Model PARTII - CONTAM86 Users Manual
Phase II Report 9. Example Problems
From a consideration of mass continuity we require w^ = w2 = w and therefore
equations (9.1) may be rewritten in expanded form as;
w C, - w C2 + V = G, (9 -2a)
-w C, * w C2 + V2 = GZ (9 .2b)
at
With these equations in hand we shall proceed to consider three cases;
Case 1 : Contaminant Decay under Steady Flow Conditions
Case 2: Contaminant Decay under Unsteady Flow Conditions
Case 3: Contaminant Dispersal Analysis of an Experimental Test
In all three cases, system characteristics will be based on those of an
experimental test reported by Traynor, et. al [3] involving measurements of
pollutant emissions from portable kerosene heaters.
9.1 .1 Case 1 : Contaminant Decay under Steady Flow Conditions
Consider the particularly simple, and familiar, case of contaminant decay from
some initial value, ^(1=0), under steady flow conditions, w = constant, with
concentration in the exterior environment maintained at the zero level, C2 = 0.
Under these conditions equation (9.2a) simplifies to;
W C, - V = 0 (9.3)
dt
whose exact solution is;
C, . C,(t=0) / (9.4)
9-2
-------
NBS: Indoor Air Quality Model
Phase II Report
PART1I - CONTAM86 Users Manual
9. Example Problems
(the quotient (V^w) is commonly know as the time constant of the system).
This exact solution is compared, below, to approximate solutions generated
with the program CONTAM using integration time steps of At = 2.0, 1.0, and 0.5
hrs with 0^=0) = 1.0 x 10'6 kg / kg air, V, = 31.87 kg, and w = 12.75 kg/hr (i.e.,
0.4 air changes per hour).
i.o
0.9
0.8
0.7
0.6
Contaminant
Concentration 0.5
0.4
0.3
0.2
0.1
0.0
— Exact
0.5hr
— At = 1.0 hr
— At = 2.0hr
Fig. 9.2 Single Zone Model: Contaminant Decay under Steady Flow Conditions
The accuracy of the general predictor-corrector method used to approximate the
response of this system is related to the time constant of the system being
studied. In this case the time constant is (31.87 kg/12.75 kg/hr) = 2.5 hr. From
the results of this single study, then, it appears that using an integration time
increment equal to a fraction of the system time constant will assure practically
accurate results.
Case 1: Command/data Input File for At = 0.5
The CONTAM command/data file and resulting results output file are listed
below. It should be noted that a large number was used for the volumetric mass
9-3
-------
NBS: Indoor Air Quality Model
Phase II Report
PARTII - CONTAM86 Users Manual
9. Example Problems
of the exterior "zone" to affect a model of a practically infinite contaminant sink.
FLOWSYS N=2
2 BC=C
END
FLOWELEM
1 1=1,2
2 1=2,1
END
FLOWDAT
TIME=0
1 W=12.75
2 W=12.75
<
TIME=15
1 W=12.75
2 W=12.75
END
EXCITDAT
TIME=0
1 CG=0.0
2 CG=0.0
<
TIME=15
1 CG=0.0
2 CG=0.0
END
DYNAMIC
T=0,10,0.5
1 V=31.87
2 V=1.0E+9
<
1' IC=1.0E-06
2 IC=0.0
END
RETURN
< Single-Zone (2-Node) Example
< Flow Element 1
< Flow Element 2
< Element Mass Flow Rates [=] kg/hr
< Nodal Excitation
< Node 1: Zero Generation Rate [=] kg/hr
< Node 2: Zero Concentration [=] kg C02/kg
< Node 1: Zero Generation Rate [=] kg/hr
< Node 2: Zero Concentration [=] kg C02/kg
< Initial Time, Final Time, Time Step Increment
< Node 1: Volumetric Mass [=] kg
< Node 2: Volumetric Mass [=] kg
< Node 1: Initial Concentration [=] kg C02/kg
< Node 2: Initial Concentration [=] kg C02/kg
Case 1: Results Output File
I CONTAM: Contaminant Dispersal Analysis for Building Systems
Jim Axley
===== FLOWSYS: FLOW SYSTEM CONTROL VARIABLES
Number of flow system nodes 2
== Node Boundary Conditions
Negative Eqtn-# = concentration-prescribed boundary.
Positive Eqtn-# = generation-prescribed boundary.
I
—Ver-10-86 —
Cornell & NBS
MTOT: 50000
9-4
-------
NBS: Indoor Air Quality Model
Phase II Report
PARTII • CONTAM86 Users Manual
9. Example Problems
Node Eqtn-# Node Eqtn-# Node Eqtn-f Node Eqtn-# Node Eqtn-#
11 2-2
===== FLOWELEM: FLOW ELEMENTS
Elem
• 1
2
I-Node
1
2
J-Node
2
1
Filter Efficency
.000
.000
== FLOWDAT: ELEMENT FLOW TIME HISTORY DATA
=== Generation Control Variables
Initial time
Final time
Time step increment
.000
15.0
15.0
Element Mass Flow Time History Data
Time':
.000
Elem Value
1 12.7
== Time: 15.0
Elem Value Elem Value
2 12.7
Elem Value Elem Value
Elem Value Elem Value Elem Value Elem Value
1 12.7 2 12.7
EXCITDAT: EXCITATION TIME HISTORY DATA
Generation Control Variables
Elem Value
Initial time
Final time
Time step increment
.000
15.0
15.0
== Nodal Excitation Time History Data
== Time: .000
"*" = independent DOFs
Node Value Node Value
1 .000 2* .000
Node
== Time: 15.0
"*" = independent DOFs
Node Value Node Value
1 .000 2* .000
== DYNAMIC: DYNAMIC SOLUTION
Node
"U" = undefined DOFs.
Value Node Value Node Value
"U" = undefined DOFs.
Value Node Value Node Value
9-5
-------
NBS: Indoor Air Quality Model PARTII - CONTAM86 Users Manual
Phase II Report 9. Example Problems
== Solution Control Variables
Initial time 000
Final time • 10.0
Time step increment 500
Integration parameter: alpha ... .750
Results print interval 1
== Nodal Volumetric Mass
"*" = independent DOFs "U" = undefined DOFs.
Node Value Node Value Node Value Node Value Node Value
1 31.9 2* 0.100E+10
== Initial Conditions: Nodal Concentrations
"*" = independent DOFs "U" = undefined DOFs.
Node Value Node Value Node Value Node Value Node Value
1 0.100E-05 2* .000
== Element Flow Rate Update ================•================== Time: .000
Elem Value Elem Value Elem Value Elem Value Elem Value
1 12.7 2 12.7
== Net Total Mass Flow
"*" = independent DOFs "U" = undefined DOFs.
Node Value Node Value Node Value Node Value Node Value
1 .000 2* .000
"*" = independent DOFs "U" = undefined DOFs.
• Node Value Node Value Node Value Node Value Node Value
1 .000 2* .000
== Time Step Estimate for Initial Conditions
— NOTE: Estimated time step to limit error to approx. 5.00% is: .925
Specified time step is: .500
"*" = independent DOFs "U" = undefined DOFs.
Node Value Node Value Node Value Node Value Node Value
1 0.826E-06 2* -0.593E-30
"*" = independent DOFs "U" = undefined DOFs.
9-6
-------
NBS: Indoor Air Quality Model PART1I - CONTAM86 Users Manual
Phase II Report 9. Example Problems
Node Value Node Value Node Value Node Value Node Value
1 0.682E-06 2* -0.187E-29
"*" = independent DOFs "U" = undefined DOFs.
Node Value Node Value Node Value Node Value Node Value
1 0.564E-06 2* -0.372E-29
== Response ==============================<==================== Time: 2.00
"*" = independent DOFs "U" = undefined DOFs.
Node Value Node Value Node Value Node Value Node Value
1 0.466E-06 2* -0.603E-29
"*" = independent DOFs "U" = undefined DOFs.
Node Value Node Value Node Value Node Value Node Value
1 0.385E-06 2* -0.874E-29
== Response ================================================== Time: 3.00
"*" = independent DOFs "U" = undefined DOFs.
Node Value Node Value Node Value Node Value Node Value
1 0.318E-06 2* -0.118E-28
( et cetera )
"*" = independent DOFs "U" = undefined DOFs.
Node Value Node Value Node Value Node Value Node Value
1 0.219E-07 2* -0.686E-28
9.1.2 Case 2: Contaminant Decay under Unsteady Flow Conditions
To investigate the consequence of unsteady flow on the nature of the behavior
of the "real" system and the numerical characteristics of its simulation we shall
extend Case 1 by considering the decay of a contaminant under conditions of
linearly increasing flow rates, that is to say with;
w = w° t ; t > 0.0 (9.5)
The decay problem is now governed by the equation;
9-7
-------
w° t C, + V^ = 0
dt
or
w°tdt = Vr^
C,(t=0) = 1.0
C,(t=0) = 1.0
NBS: Indoor Air Quality Model PARTII - CONTAM86 Users Manual
Phase II Report 9, Example Problems
(9.6a)
(9.6b)
The second form, with variables t and C1 separated, may be integrated directly
to obtain the exact solution;
. t2
C, = 1.0 e (2V,/w°) (6.7)
Again this exact solution is compared to approximate solutions generated with
the program CONTAM86, below. For this case, however, the numerical
consequences of both integration time step, At, and step-wise approximation of
the unsteady flow, Atw, (i.e., the flow approximation time step) can be
considered. (The solution was generated for V-i = 31.87 kg, and w° = 3.187
kg/hr2.)
In this case, using an integration time step equal to the flow approximation time
step, At = Atw, (i.e., updating the system flow matrix at each time step) provides
practically accurate results for even the relatively large time step of 2.0 hr (see
Figure 9.3). Updating the system flow matrix every other time step introduces
an offset error equal to the flow approximation time step (when compared to
results obtained with updating at each time step) for the first time step that is
gradually diminished with each successive time step (see Figure 9.4). This
initial offset error results because of the initial zero flow condition; in other cases
the initial error would not be expected to be as great.
9-8
-------
NBS: Indoor Air Quality Model
Phase II Report
PART1I - CONTAM86 Users Manual
9. Example Problems
OQ .
.7
OQ .
.0 '
O~l .
.(
Of .
Contaminant
OA .
.*T.
0-r .
,9
OO .
.£
01 .
.1 '
n n .
N^
•Jk
"^
N
»»
01 23456789 10
Time (hr)
Fig. 9.3 Single Zone Contaminant Decay under Unsteady Flow Conditions
with Flow Updating at Each Integration Time Step
1 .U •
OQ .
.7 '
0 8 •
0~j .
. ( '
Of .
.O
Contaminant
Concentration 0.5 •
(ug/g)
f\ A .
U.*t •
OTt .
.O
Of\ .
.4. •
O\ .
.1
On .
Xs
V
\
\
^N
^
v\
V
\
v\
"^ ^
— Exact
— At, Atw= 1.0, 1.0 hr
— At, Atw = 1 .0, 2.0 hr
V.
VSN
^>
^>
,w
^^^^fa-^ ^
==•;--
a^a^
01234567891
Time (hr)
Fig. 9.4 Sinale Zone: Contaminant Decay under Unsteadv Flow Conditions
With Flow Updating at Even/ Other Integration Time Step
9-9
-------
NBS: Indoor Air Quality Model
Phase II Report
PARTII - CONTAM86 Users Manual
9. Example Problems
Case 2: Command/data Input File for At = 1.0 and Atw = 2.0
The CONTAM command/data file used for one of these studies is listed below.
It should be noted that a large number was used for the volumetric mass of the
exterior "zone" to affect a model of a practically infinite contaminant sink.
FLOWSYS N=2
2 BC=C
END
FLOWELEM
1 1=1,2
2 1=2,1
END
FLOWDAT
TIME=0
1 W=0.0
2 W=0.0
<
TIME=12
1 W=38.244
2 W=38.244
END
EXCITDAT
TIME=0
1 CG=0.0
2 CG=0.0
<
TIME=15
1 CG=0.0
2 CG=0.0
END
DYNAMIC
T=0,10,1.0
1 V=31.87
2 V=1.0E+9
<
1 IC=1.0E-06
2 IC=0.0
END
RETURN
T=0,12,2
< Single-Zone (2-Node) Example
< Flow Element 1
< Flow Element 2
< Element Mass Flow Rates [=] kg/hr
< t=0 : w = 3.187 X 0.0 = 0.0
< t=12 : w = 3.187 X 12 = 38.244
< Nodal Excitation
< Node 1 : Zero Generation Rate
< Node 2: Zero Concentration
[=0 kg/hr
[«=] kg C02/kg
< Node 1: Zero Generation Rate [=] kg/hr
< Node 2: Zero Concentration [=] kg C02/kg
< Initial Time, Final Time, Time Increment
< Node 1: Volumetric Mass [=] kg
< Node 2: Volumetric Mass [=] kg
< Node 1: Initial Concentration [=] kg CO2/kg
< Node 2: Initial Concentration [=] kg C02/kg
9.1.3 Case 3: Contaminant Dispersal Analysis of an Experimental Test
As noted above Traynor, et.al. reported the time variation of contaminant
concentrations in a single zone system generated by portable kerosene
heaters. In this example the variation of NO concentration, C1 , in a single zone
system is computed, using measured properties of the system and NO
generation rate, and compared to experimental results. The properties of the
9-10
-------
NBS: Indoor Air Quality Model
Phase II Report
PARTII - CONTAM86 Users Manual
9. Example Problems
system and excitation used in the model are as follows;
V^ single zone volumetric mass = 31.87 kg (based on the reported volume of
27 m3 and an assumed air density of 1.18 03 kg/m3 corresponding to 26 °C
and 1 atm) .
GT : NO generation rate = 0.000186 kg/hr constant for one hour, zero
thereafter (based on the product of the reported emission rate of 23.7 |ig/kJ
times the fuel consumption of 7830 kJ/hr)
V2 : exterior "zone" volumetric mass = 1.0 - 109 kg (infinite sink modeled as a
large number)
C2 : exterior "zone" ambient concentration = 0.0 kg NO/kg air (based on
reported initial conditions)
w: air mass flow rate = 12.43 kg/hr (based on reported air change rate of 0.39
ACH)
Experimental results are compared below, Figure 9.66, to analytical results
using two integration time steps. The reported generation rate time history is
shown in Figure 9.5.
01 ^ .
.10
NO
Generation rt ,_
Rate °-10'
(g/hr)
On^ .
.uo •
Onn .
0.0 0.5 1
0 1 .5 2.
Time (hr)
Fig. 9.5 NO Generation Rate Time History Models
9-11
-------
NBS: Indoor Air Quality Model
Phase II Report
PARTII - CONTAM86 Users Manual
9. Example Problems
NO
Concentration 2.5
•*• Analytical: At = 0.1 hr
•°- Analytical: At = 0.01 hr
0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0
Time (hr)
Fig. 9.6 Single Zone: NO Contaminant Dispersal Analysis of an Experimental
Test
Traynor, et. al. also studied the time variation of CO2 concentration generated
by portable kerosene heaters in the same single zone system. Experimental
results for one of these studies are compared to analytical results below, Figure
9.7. Again, the predicted results agree well with measured data.
9-12
-------
NBS: Indoor Air Quality Model
Phase II Report
PART1I - CONTAM86 Users Manual
9. Example Problems
1 OUUU •
i Anon .
1 'rUUU
i ?nnn •
1 rtArtrt .
1 UUUU
C02
uoncenTrdTion ouuu
((Jlg/g)
-------
NBS: Indoor Air Quality Model
Phase II Report
PARTII • CONTAM86 Users Manual
9. Example Problems
EXCITDAT
TIME=0.0
1 CG=0.000186
2 CG=0.0
<
TIME=1.0
1 CG=0.0
2 CG=0.0
<
TIME=3.5
1 CG=0.0
2 CG=0.0
<
END
DYNAMIC
T=0,2,0.1
1 V=31.87
2 V=1.0E+9
<
1 IC=0.0
2 IC=0.0
END
RETURN
< Nodal Excitation
< Node 1: Generation Rate [=] kg/hr
< Node 2: Concentration [=] kg NO/kg
< Node 1: Generation Rate [=] kg/hr
< Node 2: Concentration [=] kg NO/kg
< Node 1: Generation Rate [=] kg/hr
< Node 2: Concentration [=] kg NO/kg
< Initial Time, Final Time, Time Increment
< Node 1: Volumetric Mass [=] kg
< Node 2: Volumetric Mass [=] kg
< Node 1: Initial Concentration [=] kg NO/kg
< Node 2: Initial Concentration [=] kg NO/kg
9.2 Two Zone Example
In another study Traynor et. al. [4] studied the variation of contaminant
concentration generated by portable kerosene heaters in a multi-room
residence that was modeled as a two-zone flow system. In this study a
kerosene heater was placed in a master bedroom that was allowed to
exchange air with the rest of the house and the exterior environment under a
variety of test conditions. Here we shall attempt to model one of these tests that
allowed relatively large flow rates between the master bedroom and the rest of
the house.
For this test Traynor et. al. report the time history of the flow rate between the
master bedroom and the rest of the house, the whole-house infiltration rate, and
the volumes of the master bedroom and the rest of the house. The contaminant
generation rate produced by the kerosene heater was reported in the earlier
study discussed above. The heater was operated for a period of 133 minutes.
Based on these reports a two-zone building and its corresponding flow model
may be formulated as illustrated below (Figure 9.8).
9-14
-------
NBS: Indoor Air Quality Model
Phase II Report
PART1I - CONTAM86 Users Manual
9. Example Problems
Exterior "Zone"
\
Vol, = 205m3 Vol2 = 31m3
Fia. 9.8 A Two Zone Building and Corresponding Flow Model
It will be assumed that infiltration will be equal to exfiltration for each zone (i.e.,
w3 = w4 and ws = w6) given by the product of the reported whole-house
infiltration rate (0.35 ACH) and the respective volumetric masses. The average
indoor air temperature of 16 °C will be used to compute volumetric mass
quantities and mass flow rates from the reported values (i.e., a constant density
of 1.22 kg/m3 is assumed for air).
The "inter-room" mass flow rate time histories (i.e., w^t) or equivalently w2(t)),
based on the reported volumetric flow rate histories, are plotted below along
with the computed variation of CO2 concentration in each zone, figures 9.9 and
9.10.
5000
4000
Mass Flow Rate
(kg/hr)
3000
0.0
1.0 1.5
Time (hr)
Fia. 9.9 Two Zone Example: Inter-Room Mass Flow Rate
9-15
-------
NBS: Indoor Air Quality Model
Phase II Report
PART1I - CONTAM86 Users Manual
9. Example Problems
4000
3500
Peak Value: 3769 Hg/g (2480 ppm),
C02
Concentration
fog/g)
0.0
1.0 1.5
Time (hr)
2.0
2.5
Fia. 9.10 Two Zone Example: Response Based on Measured Flow and COo
Generation Data
The peak CO2 concentration measured during the test was 3709 u.g/g (2440
ppm) that compares very well with the predicted concentration of 3769 u,g/g
(2480 ppm). It should be noted, however, that the reported flow rates were
determined to an accuracy of only ± 33 % so the close agreement of
experimental and analytical peak values must be considered to be largely
fortuitous.
Traynor et. al. also reported inter-room temperature differences for the test
considered above which suggested thermal equilibrium had been achieve by
the time the heater was shut off (i.e., the temperature difference between the
master bedroom and the rest of the house remained relatively steady. Based
on this observation the inter-room mass flow rate was assumed to have also
reached steady state (i.e., the rightmost extrapolated portion of Figure 9.9
above) for the purposes of analysis.
It is interesting, then, to consider a hypothetical extension of this test - How
would CO2 concentration vary under these (apparently) steady conditions? To
answer this question an additional analysis was computed using the flow time
history reported above (Figure 9.9), with flow assumed constant after 1.7 hours,
and a constant generation rate (i.e., without shutting off the heater). The results
9-16
-------
NBS: Indoor Air Quality Model
Phase II Report
PARTII - CONTAM86 Users Manual
9. Example Problems
of this study are plotted below. The program CONTAM, in this instance, was
used to estimate both the steady state and the dynamic response of the system.
7000
6000
sooo
C02 4000
Concentration
3000
2000
1000
Stea
/
/,
'-
dy Sta1
f
t~
e Asyr
..-''
nptotes
^^
4^
tff"^
—
r*^
•^•^
Node 1
Node 2
Node 1 - Steady State
Node 2 - Steady State
0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0
Time (hr)
Fig..9.11 Two Zone Example: Hypothetical Constant CO? Generation Rate
Response
Command/data Input File •
The CONTAM command/data file used for the first study is listed below. It
should be noted that a large number was used for the volumetric mass of the
exterior "zone" to affect a model of a practically infinite contaminant sink.
FLOWSYS N=3
3 BC=C
END
FLOWELEM
1 1=2,1
2 1=1,2
3 1=1,3
4 1=3,1
5 1=2,3
6 1=3,2
END
FLOWDAT
TIME=0
1 W=0
2 W=0
< Two-Zone (3-Node) Example
< Exterior "Zone" (Node 3) Will Have Cone. Specified
< Flow Element 1
< Flow Element 2
< Flow Element 3
< Flow Element 4
< Flow Element 5
< Flow Element 6
T=0,180/60, 0.1 < Element Ma.ss Flow Rates [=] kg/hr
< Inter-Room Flow
< Inter-Room Flow
9-17
-------
NBS: Indoor Air Quality Model
Phase II Report
PARTII - CONTAM86 Users Manual
9. Example Problems
3 W=0.35*205*1.22
4 W=0.35*205*1.22
5 W=0.35*31*1.22
6 W=0.35*31*1.22
TIME=28/60
1 W=250*1.22
2 W=250*1.22
3 W=0.35*205*1.22
4 W=0.35*205*1.22
5 W=0.35*31*1.22
6 W=0.35*31*1.22
TIME=52/60
1 W=500*1.22
2 W=500*1.22
3 W=0.35*205*1.22
4 W=0.35*205*1.22
5 W=0.35*31*1.22
6 W=0.35*31*1.22
TIME=76/60
1 W=1205*1.22
2 W=1205*1.22
3 W=0.35*205*1.22
4 W=0.35*205*1.22
5 W=0.35*31*1.22
6 W=0.35*31*1.22
TIME=101/60
1 W=3375*1.22
2 W=3375*1.22
3 W=0.35*205*1.22
4 W=0.35*205*1.22
5 W=0.35*31*1.22
6 W=0.35*31*1.22
TIME=210/60
W=3375*1.22
W=3375*1.22
W=0.35*205*1.22
35*205*1.22
35*31*1.22
35*31*1.22
W=0
W=0
W=0
END
EXCITDAT
TIME=0.0
2 CG=0.549
3 CG=0.000760
TIME=133/60
2 CG=0.0
3 CG=0.000760
TIME=210/60
2 CG=0.0
< 0.35 ACH
< 0.35 ACH
< 0.35 ACH
< 0.35 ACH
< Inter-Room Flow
< Inter-Room Flow
< 0.35 ACH
< 0.35 ACH
< 0.35 ACH
< 0.35 ACH
< Inter-Room Flow
< Inter-Room Flow
< 0.35 ACH
< 0.35 ACH
< 0.35 ACH
< 0.35 ACH
< Inter-Room Flow
< Inter-Room Flow
< 0.35 ACH
< 0.35 ACH
< 0.35 ACH
< 0.35 ACH
< Inter-Room Flow
< Inter-Room Flow
< 0.35 ACH
< 0.35 ACH
< 0.35 ACH
< 0.35 ACH
< Inter-Room Flow
< Inter-Room Flow
< 0.35 ACH
< 0.35 ACH
< 0.35 ACH
< 0.35 ACH
< Nodal Excitation
< Node 2: Generation Rate [=] kg/hr
< Node 3: Exterior C02 Concentration [=] kg C02/kg
< Kerosene heater turned off at 133 minutes.
< Node 2: Generation Rate [=] kg/hr
< Node 3: Exterior C02 Concentration [=] kg C02/kg
< Node 2: Generation Rate [=] kg/hr
9-18
-------
NBS: Indoor Air Quality Model
Phase II Report
PARTII - CONTAM86 Users Manual
9. Example Problems
3 CG=0.000760 < Node 3: Exterior C02 Concentration [=] kg C02/kg
END
DYNAMIC
1=0,150/60,0.1
1 V=205*1.22
2 V=31*1.22
3 V=1.0E+09
<
1 100.000760
2 IC=0.000760
3 IC=0.000760
END
RETURN
< Initial Time, Final Time, Time Increment
< Node 1: Volumetric Mass [=] kg
< Node 2: Volumetric Mass [=] kg
< Node 3: Exterior Volumetric Mass [=] kg
< Node 1: Initial Concentration [=] kg C02/kg
< Node 2: Initial Concentration [=] kg C02/kg
< Node 3: Initial Concentration [=] kg C02/kg
9.7 Full-Scale Multi-zone Residential Example
To provide an example of a more complex multi-zone problem consider the
hypothetical full-scale residential flow system illustrated below. In this example,
CO2 generated in one room of a two story four room residence is dispersed
throughout the building by the hot-air system and diluted by outside air
infiltration at the rate of 0.5 ACH in the two lower rooms. The CO2 is generated
by a portable kerosene heater, whose generation characteristics are assumed
to be the same as that used above in the single zone examples, is operated for
133 minutes and then turned off. The results of the analysis are plotted below
illustrating the detailed dynamic variation of pollutant concentration in the
building air flow system.
9-19
-------
NBS: Indoor Air Quality Model
Phase II Report
PART1I - CONTAM86 Users Manual
9. Example Problems
m-
0.55
kg/hr
760 pg C02/g air
20 mVhr
70 mVhr
• — Zone Node
2 — Node Number
70 mVhr
Flow Element—
Element Flow -
Actual Building
Air-Flow Sustem Idealization
Fio. 9.12 Full-Scale Residence and Corresponding Flow Model
9000
8000
7000
6000
C02 5000
Concentration
g/g 4000
3000
2000
1000
0
V^ G-'
-------
NBS: Indoor Air Quality Model
Phase II Report
PARTII - CONTAM86 Users Manual
9. Example Problems
Command/data Input File
The CONTAM command/data input file used for this study is listed below.
FLOWSYS N=7
7 BC=C
END
FLOWELEM
1 1=1,2
2 1=1,3
3 1=7,2
4 1=2,7
5 1=7,3
6 1=3,7
7 1=2,4
8 1=3,5
9 1=4,6
10 1=5,6
11 1=6,1
END
TIMECONS
1,2 W=70*1.2
3,6 W=20*1.2
7,10 W=70*1.2
11 W=140*1.2
1 V=l. 2*1.0
2,3 V=l. 2*40.0
4,5 V=l. 2*30.0
6 V=l. 2*0.1
7 V=1.2*1.0E+06
END
FLOWDAT
TIME=0
1,2 W=70*1.2
3,6 W=20*1.2
7,10 W=70*1.2
11 W=140*1.2
< Six-Zone (7-node) Example
< Exterior "Zone" (Node 7) Will Have Cone. Specified
< Flow Element 1
< Flow Element 2
< Flow Element 3
< Flow Element 4
< Flow Element 5
< Flow Element 6
< Flow Element 7
< Flow Element 8
< Flow Element 9
< Flow Element 10
< Flow Element 11
< 0.50 Building ACH each
< 0.25 Room ACH each
< 0.50 Building ACH each
< 1.00 Building ACH
< Node 1: Volumetric Mass [=] kg
< Nodes 2 & 3: Volumetric Mass [=] kg
< Nodes 4 & 5: Volumetric Mass [=] kg
< Node 6: Volumetric Mass [=] kg
< Node 7: Exterior Volumetric Mass [=] kg
< Element Mass Flow Rates [=] kgm/hr
< 0.50 Building ACH each
< 0.25 Room ACH each
< 0.50 Building ACH each
< 1.00 Building ACH
TIME=5
1,2 W=70*1.2
3,6 W=20*1.2
7,10 W=70*1.2
11 W=140*1.2
END
EXCITDAT
TIME=0
2 CG=0.549
7 CG=0.000760
TIME=133/60
2 CG=0.0
7 CG=0.000760
< 0.50 Building ACH each
< 0.25 Room ACH each
< 0.50 Building ACH each
< 1.00 Building ACH
< Nodal Excitation
< Node 2: Generation Rate [=] kg/hr
< Node 7: Exterior C02 Concentration [=] kg C02/kg
< Kerosene Heater Turned Off at 133 minutes
< Node 2: Generation Rate [=] kg/hr
< Node 7: Exterior C02 Concentration [=] kg C02/kg
9-21
-------
NBS: Indoor Air Quality Model
Phase II Report
PARTII - CONTAM86 Users Manual
9. Example Problems
TIME=5
2 CG=0.0
3 CG=0.000760
<
END
DYNAMIC
T=0,4,0.5
1 V=l.2*1.0
2,3 V=l.2*40.0
4,5 V-l.2*30.0
6 V-l.2*0.1
< Node 2: Generation Rate [=] kg/hr
< Node 3: Exterior CO2 Concentration [=] kg C02/kg
< Initial Time, Final Time, Time Increment
< Node 1: Volumetric Mass [=] kg
< Nodes 2 & 3: Volumetric Mass [=] kg
< Nodes 4 & 5: Volumetric Mass [=] kg
< Node 6: Volumetric Mass [=] kg
7 V=1.2*1.0E+09 < Node 7: Exterior Volumetric Mass [=] kg
1,7 IC=0.000760 < Initial Concentration [=] kg C02/kg
END
RETURN
It will be noticed that, in this case, system time constants were to be computed.
The results of the time constants analysis are listed below;
TIMECONS: TIME CONSTANTS - CONTAMINANT DISPERSAL SYSTEM
Convergence parameter, epsilon,
Element Mass Flow Rates
0.100E-15
Elem
1
6
11
Value
84.0
24.0
168.
Elem
2
7
Value
84.0
84.0
Elem
3
8
Value
24.0
84.0
Elem
4
9
Value
24.0
84.0
Elem
5
10
Value
24.0
84.0
Net Total Mass Flow
= independent DOFs
Node
1
6
Value
.000
.000
Node
2
7*
Value
.000
.000
Node
3
== Nodal Volumetric Mass
"*" = independent DOFs
Node Value Node Value Node
1 1.20 2 48.0 3
6. .120 7* 0.120E+07
Nominal Time Constants
Node Value Node
1 0.714E-02 2
6 0.714E-03 7
Value Node
.444 3
0.250E+05
-U" = undefined DOFs.
Value Node Value Node Value
.000 4 .000 5 .000
"U" = undefined DOFs.
Value Node Value Node Value
48.0 4 36.0 5 36.0
Value Node Value Node Value
.444 4 .429 5 .429
Actual Time Constants
9-22
-------
NBS: Indoor Air Quality Model
Phase II Report
PARTII. CONTAM86 Users Manual
9. Example Problems
Num. Value Num.
1 0.714E-03 2
6 3.73 . 7
Value Num. Value
0.714E-02 3 .230
-0.852E+16
Num. Value Num. Value
4 .429 5 .444
Number of iterations used
11
9-23
-------
NBS: Indoor Air Quality Model PARTII - CONTAM86 Users Manual
PhaseIIReport PARTII References
PARTII References
[1] Wilson, E. L & Hoit, M. I., "A Computer Adaptive Language for the
Development of Structural Analysis Programs," Computers & Structures, Vol.
19, No. 3, pp 321-338, 1984
[2] Eberlein, PJ. & Boothroyd, J., "Contribution 11/12: Solution to the
Eigenproblem by a Norm Reducing Jacobi Type Method," Handbook for
Automatic Computation: Volume II: Linear Algebra, Wilkinson, J.H. &
Reinsch,& Reinsch, C. - editors, Springer-Verlag, 1971
[3] Traynor, G.W., Allen, J.R., Apte, M.G., Girman, J.R., & Hollowell, C.D.,
"Pollution Emissions from Portable Kerosene-Fired Space Heaters",
Environmental Science & Technology, Vol. 17, June 1983, pp.369-371
[4] Traynor, G.W., Apte, M.G., Carruthers, A.R., Dillworth, J.F., Grimsrud, D.T., &
Thompson, W.T., "Indoor Air Pollution and Inter-Room Pollutant Transport
Due to Unvented Kerosene-Fired Space Heaters,", Lawrence Berkeley
Laboratory - University of California, Applied Science Division, LBL-17600,
Feb., 1984
Ref II -1
-------
NBS: Indoor Air Quality Model
Phase II Report
Appendix - FORTRAN 77 Source Code
Appendix - FORTRAN 77 Source Code
The program CONTAM86 is listed below. In
this listing you will note that compiler directives
to "includVcode stored in separate "include
files" are used. These "include files" contain
common block data specifications that are
shared by many subroutines. The contents of
these include files are listed on the last page of
this appendix.
C— 3.0 COMMAND PROCESSOR LOOP
PROGRAM CONTAM
TORSION FY»6
Developed by JAKES AXLE*
Dapt. of Architecture. Cornell University
Building Environment Division, NBS
Fall, 1986
Uiing;
A) CAL-SAP Library at subroutines developed by ED WILSON.
U.C. BERKELEY
B) Microsoft FORTRAN V2.2 Compiler for Apple Macintosh
For Nee
1. Set logical unit numbers, in SUBROUTINE INITIO, es>
NTR - 9 ; NTW - 9 ; NCMD - 9
2. INCLUDE statements use <(ilenam».INC (i.e., without •)
I. In SUBROUTINE PROMPT nsei WRITE(NTN, • (A,\) 'I STRING
C) IBM PC Professional FORTRAN (Ryen-McFarland)
1. Set logical unit numbers, in SUBROUTINE INITIO, as;
NTR - 5 I NTW - ( ; NCMD - 9
2. INCLUDE statements use '.INC' (i.e., uith ')
3. In SUBROUTINE PROMPT usei WRITE (NTW. • (A) •) STRING
Memory for dynamically allocated/defined array! i« located in
vector lA(MTOT) in blank eenRDn. To increase or decrease this
area alter the dimension of IA, in the section 0.0 balov, set
MTOT, in section 1.0 balo«. equal to this new dinsnaion, and
reconcile the code. As integers are 4 bytes wide, memory
dedicated to IA(MTOT) is equal to MTOT*4 bytes.
C - 3.1 CHECK BLANK COMMON STORAGE
C
30 NSTOR - (IDIF.-NEXT-20)>IP(1)/IP(2)
IF(NSTOR.LE.IOO) WRITE (MTV, 2100) NSTOR
2300 FORMAT!
*• •••• WARNING! Array storage available -•,!>,• real numbers.')
C
C - 3.2 GET COMMAND LINE
C
IF (MODE. EC..' INTER ') CALL PROMPT C CMND>')
CALL FREE
IFIMOOE.EO.'llATCH'l CALL FREEWR(NTN)
C
C - 3.3 INTERPRET COMMAND LINE
C
C - GET COMMAND i ARRAY NAMES, IF ANY
C
CALL FREECC •,NCMND,e,l)
CALL FREECC*', N1<1>, 4, 7)
IMPLICIT REAL>< (A-H.O-Z)
C—0.0 DATA SPECIFICATIONS I COMMON STORAGE
C
COMMON MTOT.NP,IA(20000)
INCLUDE ARYCOM.INC
INCLUDE lOCOM.INC
INCLUDE CMDCOM.INC
INCLUDE CNTCON86.INC
LOGICAL ERR •
C—1.0 INITIALIZE INTERNAL VARIABLES
C
MTOT - 20000
CALL INITAR(HTOT)
CALL INITIO
CALL INITCN
ERR - .FALSE.
C—2.0 WRITE BANNER
C
CALL BANNER (NTW)
CALL BANNER (NOT)
WRITE (NOT, 2200) (FNAMEdiLTNAME)//' .OUT1)
2200 FORMAT!/' — RESULTS OUTPUT FILEi ',(A»
INTRINSIC COMMANDS
IF((NNCMKD.EQ.'a').OR.(NNCMKD.EQ.'BELP')) THEN
IFIMOOE.EO.'BATCa'l TEEN
WRITE(NTH,2310)
WRITE(NOT.2310)
CALL RETRN
ELSE
CALL EEIJ>
ENDIF
ELSEIF(NNCMKD.EO.'ECBO-ON') TEEN
ECHO - .THUS.
ELSEIF(NNCMHD.EQ.'ECBO-OPP') TBEN
ECHO - .FJOSE.
ELSEIF((NNCMND.EO.'Lt>.OR.(NNCMND.EO..'LIST')) TBEN
IF(MOOE.EO.'BATCH') TBEN
WRITE (NTH, 2310)
WRITE(NOT,2310)
CALL RETRN
ELSE
CALL LIST
ENDIF
ELSEIFUNNCMND.EQ.'P'j.OR.INNCMND.EQ.1 PRINT')) TBEN
CALL PRINT
ELSEIFKNNCMND.EO.'D') .OR. INNCMHD.EQ.'DIAGRAM')) THEN
CALL DIAG8M
ELSEIFINNCmD.EO. 'SUBMIT') THEN
IFIMODE.EO.'BATCH') THEN
WRITE(MTW,2310)
WRITE(MOT, 2310)
CALL RETRN
ELSE
CALL SUBMIT
ENDIF
ELSEIF(NNCI
-------
NBS: Indoor Air Quality Model
Phase II Report
Appendix • FORTRAN 77 Source Code
CALL FLOSYS
ELSEIF(NNCMND.EO. 'PLOVELEH') THEN
CALL FLOELM
ELSEIF(NNCMND.EO..'STEADY') MEN
CALL STEADY
ELSEIFINNCKND.Ea.'TIMECONS') THEN
CALL TIHCON
ELSEIFINNCMKD.EO.'FLOMDAT') TEEN
CALL FLODAT
ELSEIF(NNCHND.EO.'EK:iTDAT') TBEN
CALL EXCDAT
BLSEIFINNCMND.EO..'DYNAMIC1) TEEN
CALL DYNAH
SUBROUTINE INITCN
C—SUBiINITCN - INITIALIZES CONTAM LABELED COMMON /CNTCOM/
INCLUDE CNTCOMX.INC
NFNOD - 0
NFEQN - 0
NFBAM - 0
NFELM - 0
MPV - 0
MFF - 0
KPC - 0
MFC - 0
MPXXQ - 0
EP - l.OD-16
RETURN
END
ELSEIF(NNCMMD.EO. 'BESET') TBEN
CALL RESET
SUBROUTINE BANNEB(LUN)
C—SUBt BANNER - WRITES PROGRAM BANNER TO LOGICAL UNIT LUN
ELSE
NRITE(NTH. 2310)
IFIMODE.EO.'BATCB'l TBEN
CALL RETRN
ENDIF
HTOT.NP,IA(1)
•RITE(LUN.2000) MTOT
2000 FORMAT(//.1X,78(1B-)./,
.'I CONTAM9 «',T7»,'|1,/,
.' I Contaminant Diaparaal Analyai* for Building Syctarra •
ENDIF
GO TO 30
C
2310 FORMAT) ' •••• ERROR: Comnand net dafinad in BATCB moda.')
2320 FORMAT C •••• ERROR I Comand not dafinad in INTERACTIVE nx
2330 FORMAT!' •••• ERROKi Conmand not dafinad. •)
SUBROUTHIE INITARIHTOT)
C SUBtrNITAR *• INZTXALIZXS DYNAMIC ARRAY MANAGER VARIABLES
C IN BLANK COMMON AND LABELED COMMON /ARYCOH/
INCLUDE ARYCOM.INC
NUMA - 0
IDIR - NTOT
IP(1) - •
ZP(2) • 0
IP (3) - 1
RETURN
END
SUBROUTINE INITIO
C— SUB I INITIO - INITIALIZES LABELED COMMON /IOCOM/
C OPENS DEFAULT RESULTS OUTPUT FILE
INCLUDE IOCOM. INC
LOGICAL FOUND
NTR - »
NTH - »
NCMD - »
NIN • 10
NOT - 11
ND1 - 12
N01 - 13
ND3 • 14
ND1 • 15
FNAME - 'CONTAM'
LFNAME - C
EXT - '
CALL NOPEN(NOT, (FNAME(ltLFNAME)//' .OUT'), 'FORMATTED')
MODE - 'IHTEB-
ECHO - .TRUE.
RETum
END
.,D»,'I',/,
. ' 1 V«r«ion tltt - Jin Axl«y - Cornell < NBS ' ,
.T7», ' 1 ',/,lX,71(lB-),/,65X, 'MTOTi MS)
to.') RETURN
END
C
C
SUBROUTINE BBLP
C— SUBl EILP - PROVIDES ON-SCREEN HELP
C
C
C .' EELP (8) Li to continu*.'
•RITE (NTH, 2010)
PAUSE ' — Enur to eontinu*. '
NRITE (NTM, 2020)
PAUSE ' — Ent.r to eontinu.. •
•RITE (NTR, 2030)
PAUSE ' — Eater to continu*.'
•RITE (NTN, 2040)
RETURN
C
C QELF LI3T3 • • • •
C
2000 FORMAT(/, ' — - INTRINSIC COMMANDS',//,
.' HELP IB) Li Print array n«n»d .',/.
.' DIAGRAM (D) A- Diagram array nanad Ml.',/,
.' SUBMIT F- Raad conmndi from batch .
.' mnuRN La
-------
NBS: Indoor Air Quality Model
Phase II Report
Appendix • FORTRAN 77 Source Code
.' FLOWSYS N-nl
.' n2,n3,n4 BC-cl
.' END
.• FLOVELEM
.' nl I-n2,n3 E-n4
.' END
. ' STEADY
.' nlf n2,n3 W-n4
.' n5,n(,n7 CG-n8
.• END')
2020 FORMAT!/,
.• TIMECONS I-nl
. • n2.n].M W-nS
.• nC,n7,ne V-n9
.' END')
2030 FORMAT!/,
.• FLOWDAT [T-nl.n2,n3]
.' TIME-nl
.' nl,n2,n3 W-n4
.' i',/,
.' END',//,
.' EXCITDAT (T-nl,n2,n3)
.' TIME-nl
.' nl,n2,nj CG-n4
.' I',/,
.' END')
2040 FORMAT!/.
.' DYNAMIC
.' T-nl,nJ,nJ [A-n4] [PI'
.' n1.ne,n9 V-nlO
.' t
.• n7,n>,n9 IC-nll
.' END ',//.
. • RESET
Flott*y*t*m control variabl**.',/,
nl - number of {Ion nod**'./.
n2,n3,n4 • nod*i fir*t, la*t. iner.',/,
el - boundary condition; G or C',//,
Flov *l*m*nt coraoand/data group.',/,
nl • *l*ia*nt number',/,
n2,n3 - *l*m*nt *nd nod**',/,
n4 - filt*r *ffiei*ney',//.
Steady *tata colution.'./,
nl,n2,n3 » *l*mi fir*t, la*t, incr.',/,
n4 - *l*n*nt flow rat*',/,
n5,n(,n7 - nod*t fir*t. l**t, incr.',/,
n8 - pr**crib*d cone, or g*n. rat*',/,
Tim* constant *olution, Al • *p*ilon',/,
n2,n3,n4 - *l*mi Cirat, la*t, iner.',/,
n5 • *l*n*nt flow rat*',/,
nfi,n7,na • nod*t fir*t, last, iner.',/.
n9 - nodal volumetric nan',/,
G*n*rat* «l*ra*nt flow tin* histori**.',/,
nl - tin*',/,
nl,n2,n3 • *l*mt fir*t, la*t, incr.',/,
n4 • *l*a*nt oa** flow rat*.',/.
G«n*rat* vxeitation tin* hiitori**.',/,
nl - tin*',/,
nl,n2,n3 • nod*t fir*t, la*t, iner.',/,
n4 • cxeitationt eonc. or g*n. rat*.',/.
Dynamic *olution.',/,
>nS] (PS-nC)',/,
nl.n2,n3 - init,final,iner; n4 -alpha',/,
n5 - print interval; n( - plot leal*',/,
n7,nB,n9 - nod*. fir*t, la>t, iner.',/,
nlO - nodal voluaetrie na*«',/,
nil - initial nodal concentration',/,
R***t CONTAM for nev problem.')
1C - IDIR
DO 100-I-l.NUMA
II - IL + 1
ILOC - 1
1ST - 0
IA6 - IAIIC+6)
IA7 - IAIIC+7)
IA9 - IAIIC+9)
•CHECK FOR LOCATION AND STORAGE TYPE
IFIIA9.CT.O) ILOC-2
iriIA7.LT.01 ILOC-2
iriIA7.EO.-l) IST-1
ir(IA7.EO.-J) IST-2
IFIIA9.GT.O) IST-3
IPN - 1C - 1
DO 10 J-1,4
IPN - I7N f 1
NAM(J) > CH)iRIIA(IPN))
DATA TO TERMINAL
IFIIST.EQ.OI WRITE(NTH,1100) (NAM(J),J-l,4),
• IAIIOM),IA{IC+5), !TYPB!K, IA6) ,K-1, 9).
* (LOCIL,ILO<:),L-1,4)
IFIIST.EO.l) WRITE(NTH.1100) (NAM|J),J-1,4),
• IA(IO4),Ul(IC+5), (TYPEIK,IAC),X-1,9),
• ILOCIL, ILOi:). L-l, 4), ISTORIM, 1), M-l, 13)
IFIIST.EQ.2I WRITE INTW, 1300) (NAM(J), J-l, 4),
• IAIIC-M),(LOCIL.ILOO.L-1.4).(STORIM,2) ,M-1,13)
IFIIST.EQ.3) WRITE(NTW, 1200) INAMIJ), J-1,4),
• IAIIC«4).UtIC+S),IAIIC»«). |LOC(L,ILOC).L-1,4),
• (STOR(M,2),M-1,13)
1C - 1C * 10
-CHECK FOR NUMBER OF LINES PRINTED
IP(IL.LT.20) GO TO 100
If fl.EQ.NimA) CO TO 100
CALL PROMPT(• •• Do you «ant mar* 7 IY/N) •)
READ(NTR,2200)
IF!ICHK.EO.'n').OR. ICBJt.EQ.'N'» GO TO 900
II - 0
WRITEINTW.2000)
I CONTIHUE
SUBROUTINE LIST
C—SUBiLIST - LIST DIRECTORY OF ALL ARRAYS IN BLANK COMMON
C
C—HELP LIST
C
C
C-
. • LIST (L)
Li*t th. diraetory of all array*. ',/,
COMMON MTOT.NP.IAU)
INCLUDE ARYCOM.INC
INCLUDE IOCOM.INC
CHARACTER-1 NAM(4) ,LOC(4,2) .TYPEI9, 3),STORI13,2)
CHARACTER*! CHK
DATA TYPE/'I','N','T','E','G','E','R',' ',' ',
1 'R','E','A'.'L',' ',' ',' ',' ',' ',
2 'C'.'H'.'A'.'R'.'A'.'C'.'T'.'E'.'R'/
DATA LOC/'C', '0', 'R', 'E', 'D', 'I', 'S', 'KV
DATA STOR/'S','E','O'.'U'.'E','N'.'T','I','A'.'L',' ',' ',' ',
-LIST DIRECTORY OF ALL ARRAYS IN DATA BASE
IF(NUMA.EO.O) CO TO 900
-WRITE HEADER FOR SCREEN LISTING OF FILE DATA
WRITE(NTW, 1000)
START COUNT OF LINES TO SCREEN
1000 FORMAT!' — — LISTi ARRAY LIST',//,
- • Nam*',2)t,'Number',2X,'Number',5X,'Data',SX,
• 'Location'.SX, 'Storage',/,8X,'Ro»*',2X,
• 'Column«',5X,'Type',19X,'Typ«',/>
1100 FORMAT(IX, 'IA1, JX,14,4X,14,SX,9A1,4X,4A1,4X,13A1)
1200 FORMAT!1X, Print array naimd .',/.
C
COMMON MTOT.NP.IAIl)
INCLUDE ARYCOM.INC
INCLUDE IOCOM.INC
INCLUDE CHDCOM.INC
-PRINT OF PEAL OR INTEGER ARRAY
CALL PROMO 11)
-LOCATE MATRIX TO BE PRINTED
IF (ECHO) WRITE INTW, 2000) Ml
WRITE (NOT, 2000) Ml
Append-3
-------
NBS: Indoor Air Quality Model
Phase II Report
Appendix - FORTRAN 77 Source Code
CALL LOCATE (Ml,NA.NR.NC)
IF(NA.EQ.O) THIN
URITI (Hilt, 2010) Ml
WRITE(NOT,J010) m
CALL ABORT
RETURN
ELSEIFfNA.LT.O) TEEN
WRIT!(DTK, 2020) Ml
WRITS(NOT, 2020) Ml
CALL ABORT
RETURN
ELSE
IKNP.EO.l) CALL IPRT(IA(NA) ,NR,NC)
IFIKP.E0..2) CALL RPRT(IA(NA),NR.NC)
ENDI?
RETURN
C
2000 FORMAT!/'
2010 FORMAT! •
2020 FORMAT I•
END
— PRINT OF ARRAY ••,«Al,t")
••• ERROR. Array ••.4A1,"' do»> net miat.')
••• ERRORt Array ••,4X1," it out of cor..')
MUTE (NOT, 2002) J, (A(J,K).K-I,IN)
IF(ECBO) WRITS(NTH,2002) J. (A(J,X). K-I. IN)
ENDIF
100 CONTINUE
C RETURN
2000 rORMATI/1 COLI -',6112)
2001 FORMATC ROW, 14, IE12.S)
2002 FORMAT)' RON',14,«F12.5)
END
SUBROUTINE DIAGRH
C— SUBiDIACRM - COMMAND TO -DIAGRAM" ARRAY TO RESULTS OUTPUT FILE
C
C—BELP LIST
C
C .' DIAGRAM (D) A- Diagram array nan»d Ml.',/,
C
SUBROUTINE IIRT(N,NR,NC)
C—SUB: IPRT - PRINTS INTEGER ARRAY TO RESULTS OUTPUT FILE
DIMENSION N(NR,NC)
INCLUDE IOCOM.INC
NUMC - 14
DO 100 I-l.NC.NUMC
IN - I + NUMC - 1
IPIIN.CT.NC) IN - HC
KRITE(NOT,2000) (K, K-I. IN)
IF (ECHO) WRITE (NTH, 2000) (K. K-I, IN)
DO 100 J-l.NR
WRITE(NOT,2001) J,(N(J.K),K-I,IN)
IFIECBO) WRITE(NTW,2001) J, (N(J.K) .K-I.IN)
100 CONTINUE
2000 FORKkTI/' COL* -',1415)
2001 FORMATC ROW'. 14,HIS)
END
COMMON MTOT.NP.IA(l)
INCLUDE IOCOM.INC
INCLUDE CMDCOM.INC
-PRINT OF REAL OR INTEGER ARRAY
CALL PROMS(II
-LOCATE MATRIX TO BE PRINTED
IP (ECHO) WRITE (Km, 2000) HI
KRITE (NOT, 2000) Ml
CALL LOCATE (HI, NA.NR.NC)
IF(NA.EO.O) TBSN
WRITS(NTW,2010) m
WRITE (NOT, 2010) HI
CALL ABORT
RETURN
SLSEIF(NA.LT.O) .TEEN
•RITE (MTV. 2020) n
WRITE(NOT,20JO) HI
CALL ABORT
RETURN
ELSE
IFIN7.E0.1) CALL IDIASR(IA(NA) ,NR,NC>
IFIHP.IQ.2) CALL RDIAGR(IA(NA),NR,NC)
ENDIF
RETURN
SUBROUTINE RPRT(A,NR,NC)
C—SUB. RPRT - PRINTS REAL ARRAY TO RESULTS OUTPUT FILE
IMPLICIT REAL*< (A-B.O-I)
DIMENSION AINR.NC)
INCLUDE IOCOM.INC
2000 FORMAT(/'
2010 FORMATC
2020 FORMAT C
END
—DIAGRAM OF ARRAY •'.4A1,'-')
»• ERRORi Array •',4A1,'" doaa not .Ki.t.'l.
•" ERRORt Array "'.4A1," ia out of oora.')
SUBROUTINE IDIACR(N,NR,NC)
C—SUB I IDIAGR - •DIAGRAMS* INTEGER ARRAY TO RESULTS OUTPUT FILE
XMAX - 0.00
DO 90 I-l.NR
DO 90 J-l.NC
XX - DABS(A(I,J))
IF(XX.CT.XMAX) XMAX - XX
50 CONTINUE
H - 1
IP(XMAX.LT.99939.) M - 2
IP(XMAX.LT.0.1000) N - 1
IF(XMAX.EO.O.O) M - 2
NUMC - 6
DO 100 I-1,NC,NUMC
IN - I » HUME - 1
IP(IN.GT.NC) IN - NC
WRITE (NOT. 2000) (K, K-I, IN)
IF(ECBO) WRITE (NTW, 2000) (K, K-I, IN)
DO 100 J-l.NR
IF(N.Ea.l) TEEN
WRITE (NOT, 2001) J. (A( J. X). K-I. IN)
IP(ECHO) WRITE(NTW,20011 J,(A(J,K) ,K-I, IN)
ELSEIF(M.E0.2> TBEN
INTEGER N(NR.NC)
CBARACTER-1 ICONO6)
INCLUDE IOCOM.INC
: DIAGRAM INTEGER ARRAY
NUMC - 16
DO 200 I-l.NC.NUMC
IN - I + NUMC - 1
IF(IN.CT.NC) IN - NC
•RITE(NOT,2000) (INT(K/10).K-I.IN)
•RITE (NOT, 2010) «X-INT(K/10)»10).K-I.IN>
IF(ECBO) WRITE(NTW,2000) |INT(K/10),X-I, IN)
IF(ECBO) WRITE(NTW,2010) ((K-INTIK/10)-10),K-I.IN)
DO 200 J-l.NR
DO 100 K-I, IN
ICON(K) - '•'
IF(N(J,K).EQ.O> ICON(K) - ' '
100 CONTINUE
WRITE(NOT,2020) J,(ICON(K),K-I,IN)
IF(ECHO) WRITE(NTW, 2020) J,(ICON(K),K-I,IN)
200 CONTINUE
Append-4
-------
NBS: Indoor Air Quality Model
Phase II Report
Appendix - FORTRAN 77 Source Code
RETURN
2000 FORMAT(/• COLI -',36(1X,II))
2010 FORHAT(7X,3S(1X.I1))
2020 FORMATC ROW1,13,36 (IX, AD)
END
SUBROUTINE RDIACR (A, NR, NC)
C—SUB, TOUCH - -DIAGRAMS" REAL ARRAY TO RESULTS OUTPUT FILE
SUBROUTINE RETRN
C—SUBtXKTRN - RETURNS TO INTERACTIVE MODE
C
C—HELP LIST
C
C .' RETURN Laat command in batch .', /,
C
INCLUDE IOCOH.INC
REAL-8 A(KR.KC)
CHARACTER-1 ICON(36)
INCLUDE IOCOM.INC
C DIAGRAM INTEGER ARRAY
NUMC - 36
DO 200 I-1,NC,NUMC
IN - I * NUMC - 1
IF(IN.CT.KC) IN - NC
•RITE (NOT. 2000) (INTIK/10) .K-I.IN)
nun (NOT, 2010) ((K-INTIK/IOI-IOI.K-I.IN)
IF (ECHO) WRITS (HTII, 2000) (INT(K/10) ,K-I. IN)
I? (ECHO) mUTE I (K-INTIK/10)-10). K-I, IN)
DO 200 J-l.NR
BO 100 K-I.IN
ICON(X) - ••'
tFIA(J.It).EO.O.OOO) ICON IK) - • •
100 CONTINUE
KSITB(NOT. 2020) J. (ICON(lt), K-I, IN)
IF(ECBO) WRITE(NTM,2020) J, (ICON(K).K-I.IN)
200 CONTINUE
C
RETURN
2000 FORMAT(/' COL* -'.36(IX, ID I
2010 FORMATC7X,36(1X.I1»
2020 FORMAT C ROW', U, 3«<1X,A1))
END
SUBROUTINE SUBMIT
C—SUB] SUBMIT - SWITCHES TO BATCH MODE AND OPENS BATCH COMMAND FILE
C
C—HELP LIST
C
C .' SUBMIT P- Read cotmuidi from batch .',/,
C
CLOSE(NCMD)
CLOSE(NOT)
FNAME - 'CONTAM'
LFNAME - 6
OPEN (NOT, FIIJXFNAMEdrLFNAME)//'. OUT') .STATUS-'OLD'.
•f FORM-' FORMATTED')
REWIND NOT
10 READ(NOT,«,EHD-20I
CO TO 10
20 BACKSPACE(NOT)
NCMD - NTR
NOOK - 'INTER*
WRITE(NTW,2010)
WRITE(NOT,2010)
2010 FORMAT!' •••• CONTAM raturnad to INTERACTIVE rod..')
RETURN
END
CONTAM COMMANDS
-FLOSYS
SUBROUTINE fLOSYS
C—SUB t FLOSYS - COMMAND TO READ t PROCESS FLOW SYSTEM CONTROL VARIABLES
C EilTABLISEES FLOW SYSTEM EQUATION NUMBERS < B.C.
C
C—HELP LIST •-
C
INCLUDE IOCOM.INC
LOGICAL FOUND
FLOWSYS N-nl
n2,nl,n4 BC-cl
n2,n3,n4 V-nS
.' END1,//,
Flovayatam control variablat.',/,
nl - numbar of flow nodaa',/,
n2,nl,n4 - nodal flrat, la
OPEN (NCMD, FILE-FNAME (ItLTNAME). STATUS-'OLD')
REWIND NCMD
CLOSE(NOT)
CALL NOPEN(NOT, (FNAME(liLFNAME)//- .OUT1), 'FORMATTED')
CALL BANNER (NOT)
WRITE(NOT,2020) (FNAME(liLFNAME)//'.OUT')
2020 FORMAT!/' — RESULTS OUTPUT FILEi '.(AD
ELSE
WRITE (NTW, 2030)
2030 FORMAT C — NOTEl Submit fila not found.')
CALL ABORT
ENDIF
COMMON HTOT.HP, IAI1)
INCLUDE IOCON.INC
INCLUDE CNTCOM86.INC
LOGICAL ERR
INTEGER UK 13)
EXTERNAL BCDATO
ERR - .FALSE.
C
C—1.0 READ NUMBER OF FLOW SYSTEM NODES
C
CALL FREEII'N'.NFNOD.l)
IF(NFNOD.LH.O) THEN
WRITE (NTII, 2100)
WRITE(NOT, 2100)
ERR - .THUS.
GO TO 40(1
ENDIF
2100 FORMATC •"•* ERRORl Nunbar of flow tyttam nodaa mutt ba graatar',
+ • than O.'t
RETURN
END
IF(MOOE.EO. 'INTER') THEN
WKITE(NT1«,2110)
WRITE (NTH, 2120) NFNOD
Append-5
-------
NBS: Indoor Air Quality Model
Phase II Report
Appendix - FORTRAN 77 Source Code
WRITS(NTH,1130}
ENDIP
WRITE {NOT, J110)
WRITE(NOT,2120) NFNOD
•RITE(NOT,2130)
2110 FORMAT!/, • — PLOTS YSI FLOW SYSTEM CONTROL VARIABLES')
2120 FORMAT]/
+' Nunbar of flow lytten node* ', 15)
2130 FORMAT!/, • — Nod* Boundary Condition!')
NTMN - NFNOD
C—HELP LIST-
C
7LOMELEH Flow element eoranend/data group.',/,
nl I-n2,n3 E-n4 nl • element number',/,
... n2,n3 - element and nodea',/,
DO n4 - filter efficiency',//,
COMMON MTOT.NP.IAIl)
INCLUDE tOCOH.INC
INCLUDE CNTCOMM.INC
C—2.0 DEFINE KZO ARRAY AND NUMBER EQUATIONS IN NODE ORDER
C
CALL DELETE('KEO ')
CALL DEPINICKEO ' .MPKEO.NFNOD, 1)
NN - 0
DO 30 N-MPKEO,MPKEO+NFNOD-1
NN • NNtl
20 IAIN) - NN
C
C—3.0 PROCESS BOUNDARY CONDITION DATA
C
CALL DATEEN(BCDATO,0,ERR)
C
C—4.0 RETORT BC IF NO ERROR ENCOUNTERED, ELSE ABORT
C
400 IF(ERR) THEN
CALL DELETECKEQ ')
MFKEO - 0
ERR - .FALSE.
CALL ABORT
ELSE
IF (ECHO) WRITE (NTH, 2400)
WRITE (NOT, 2400)
ITIECBO) WRITE(NTW,2410) (|N),IA(N+MPXEO-l),N-l.NPNOD)
NRITE(NOT,2410) ( (N), IA(N+MPKEQ-1> ,N-1,NFNOO)
ENDIF
RETURN
2400 FORMAT!/.
. 6X,'Negative Eqtn-t • concentration-preacritoed boundary.',/,
. fX, 'Poaitiva EQtfi~v m oeneration-preacribed boundary.',//,
.4X.V Nona Eqtn',2X»
2410 FORMAT((4X,3(II.1X,I«,2X)))
END
SUBROUTINE BCOATO (N, ERR)
C—SUBlBCDATO - READS FLOH B.C. DATA
C
COMMON MTOT.NP.IAU)
INCLUDE ICCOM.INC
INCLUDE CNTCOMK.INC
LOGICAL ERR
CHARACTER BC*1
CALL FREECCC'.BC.l.l)
IF((BC.NE.'C').AND.(Be.NE.'e'» THEN
WRITE (NTH, 2000) BC
WRITE (NOT. 2000) BC
ERR - .TRUE.
RETURN
ELSEIF(BC.EQ.'C') THEN
IA(N»MPKEO-1> - -IA(N*MPRM-D
ENDIF
RETURN
2000 FORMATC •••• ERRORt Boundary condition ',A1,• not available.')
END
REAL'S EFT
LOEICAL ERR
EXTERNAL FLOWED
ERR - .FALSE.
WRITE(NOT, 2000)
WRITS(NTH, 2000)
2000 FORMAT!/,' — FLONELEMi FLOW ELEMENTS')
C
C—1.0 CHECK TO SEE IF SYSTEM NODES t EQUATION NUMBERS ARE DEFINED
C
IF(NFHOD.EQ.O) TEEN
WRITE(NTW,2100)
WRITE(HOT,2100)
2100 FORMAT I
+ • •••• ERRORt Nunbar of flo. iy>tan noda> - O.',/,
+ ' FLOWSYS eonmuid nuat ba axaeutad.')
CALL ABORT
RETURN
ENDIF
C
C—2.0 OPEN .FEL
C
IF(NFBLM.EO.O)
* CALL NOPEN(NDl,(FNAHE(liLFNAME)//'.FEL'),'UNFORMATTED')
IF(NFELM.CT.O) THEN
WRITE(NTW, 2200)
WRITE(NOT,2200)
CALL ABORT
RETURN
ENDIF
2200 FORMAT!' •••• ERRORl Flow alananta nava already been defined.')
C
C—3.0 GET ELEMENT DATA
C
NELDOF - 2
CALL tLCEN(FLOWEO.IA(MPXEO>.NELDOF,NFNOD.MFBAN,ERR)
C IF ERR ABORT COMMAND
30 IF (ERR) TEEN
NFELM - 0
CALL ABORT
CLOSE (HDD
RETURN
ENDIF
C
C—4.0 REPORT ELEMENT DATA
C
REWIND (HDD
WRITE (NOT,2400>
IF(ECBO) KRITEINT*.5400)
2400 FORMAT!/, ' Elan I-Node J-Nod. Filter Effic.ncy'1
DO 40 N-l,NFELM
READ(NDl) LM1,LM2,EFF
IF (ECHO) WRITE (NTW, 2410) N, LM1, LM2, EFF
40 WRITE (NOT. 2410) N. LM1, LM2, EFF
2410 FORMAT(3(3X.IS).SX, CIO.})
C
C—5.0 CLOSE ELEMENT DATA FILE
C
CLOSE (HDD
RETURN
END
SUBROUTINE FLOELM
C—SUBlFLOELM - COMMAND TO READ I PROCESS FLOW ELEMENT DATA
C
-FLOELM c-
SUBROUTINE FLOWED (NIL. LM. ERR)
C—SUB I FLOWED - READS ADDITIONAL ELEMENT DATA
C WRITES FLOW ELEMENT DATA TO LOGICAL UNIT ND1
Append-6
-------
NBS: Indoor Air Quality Model
Phase II Report
Appendix - FORTRAN 77 Source Code
INCLUDE IOCOM.INC
INCLUDE CNTCOmC.INC
REAL'S EFF
INTEGER LM(2),NEL
LOGICAL ERR
C
C—1.0 GET FILTER EFFICIENCY
C
EF? - 0.0
CALL FREERf'E-.EFF.l)
IF(EFF.LT.O.ODO) TEEN
WRITE (NTH, 2100)
WRITE(NOT,2100)
2100 FORMAT)
• t .... EUROS, Filter •fflci«nci«> niit b» grratar Chan 0.')
ERR - .TRUE.
RETURN
ENOIF
C
C—2.0 WRITE ELEMENT INFORMATION TO ND1 -
C
WRITE(ND1) LM(l). LM(2). EFF
NFELM - NEL
C—2.0 DEFINE AND INITIALIZE ARRAYS
C
CALL DELETE)'*E
CALL DELETE CC
CALL DELETE (T
CALL DEFINRCF
CALL DEFINRCC
CALL OEFimrWB
',«PF,NFEQN,2*MPBAN-1)
'.MPG.NFEON. 1)
•.HPWE.NFELM.l)
CALL tERO!l(LMMK),NFEQN, 1)
CALL ZEROR(IMMPWE).NFELM. 1)
C
C—1.0 GET ELEMENT FLOW RATES (WE)
C
CALL READWE(KRR)
IF(ERR) THEN
CALL ABORT
RETURN
ENOIF
•4.0 FORM IF)
ALLOH *END* BEFORE EXCITATION DATA TO JUST FORM COMPACT (C)
OPEN(ND1.FIUS-(FRAME(11LFNANg) //'.FSL'), STATUS-• OLD'.
•••FORM- • UNFORMATTED •)
REWIND ND1
RETURN
END
SUBROUTINE STEADY
C—SUB I STEADY - COMMAND TO FORM STEADY PROBLEM [F] 1C) - (C) t SOLVE
C SOLUTION (C) IS WRITTEN OVER (G)
C
C—SKIP LIST
C
C .' STEADY Study «t»t« lolution.',/,
C .' nl,n2,n3 W»n4 nl.n2,n3 - •lam fixet, la*t. inex.1,/.
C .' ... n4 • •l«mnt flex xat«',/.
CC*n8
nS.ni.n7 • i
oona. ox a«n.
.• END1.//.
IMPLICIT RIAL'8 (A-fl. 0-Z)
COMMON MTOT.NP.IA(l)
INCLUDE IOCOM.INC
INCLUDE CMDCOM.INC
INCLUDE CNTCOMM.INC
LOGICAL ERR
CHARACTER ENDFLAC*]
ERR - .FALSE.
WRITE(NOT,2000)
WRITS(NTW.2000)
2000 FORMAT)/, ' — STEADY! STEADY STATE SOLUTION1)
C--1.0 CEECX IF FLOW SYSTEM AND ELEMENT DATA ARE DEFINED
C
IF(NFEON.EQ.O) TEEN
WRITE(NTW,2100)
WRITE(NOT.2100)
2100 FORMAT(
+ • .... ERROR! Nuatex of flow •yBtra DOF* • O.',/,
+ ' FLOWSYS eoRffland m«t b* ameutad.')
RETURN
ILSEIF (NFELM. EO.O) THEN
WRITE(NTW.2110)
WRITE(NOT, 2110)
2110 FORMAT!
» ' •••• ERROR! Nontax of flo« Cloo «l«n»nt« . O.',/,
+ ' FLOWELEM eonmand muct b> •x*eut*d.')
RETURN
ENDIF
C
CALL FORMFILMMPKEOl.IAIMPFj.IAIMPWE),1 BAND1,ERJt)
IF (ERR) THEN
CALL ABORT
RETURN
ENDIF
CLOSE (NOD
CALL FREECC ',ENDFLAC, 3,1)
IFIENDFLAG.EQ.fEND') RETURN
C—S.O FORM (C)
C
CALL FORME(ERR)
inex.1,/, IF(ERR) THEN
r.t.',/, CALL ABORT
RETURN
ENDIF
C
C—(.0 MODIFY (G) AND [F) FOR PRESCRIBED CONCENTRATIONS
C
CALL MODIF(HKMPKEQ), IA(MPF), IA(MPC) .NFNOD.NFEON, MFBAN)
C
C—7.0 SOLVE
C
CALL FACTCA(IA(MPF).NFEQN,MFBAN,ERR)
IF (ERR) THEM
CALL ABORT
RETURN
ENDIF
CALL SOLVCAIIA(MPF).IA(HPG).NFEON,HFBAN,ERR)
IF (ERR) THElf
CALL ABORT
RETURN
ENDIF
C
C—8.0 REPORT SOLUTION
C
IFjECBO) WRITE)NTW,2800)
WRITE(NOT,2800)
2800 FORMAT!/,' — Ruponui Nod* Concentration.•)
CALL REPRTN(IA(MPG).IA(MPXEO),NFEaN.NFNOD)
C
C—9.0 DELETE ARRAYS
C
CALL DELETE ('ME •)
CALL DELETECC •)
CALL DELETE I'F •)
RETURN
SUBROUTINE REAOWE(ERR)
Append-7
-------
NBS: Indoor Air Quality Model
Phase II Report
Appendix - FORTRAN 77 Source Code
C—SUB IREADME - READS I REPORTS ELEMENT TOTAL MASS FLOW RATE DATA
C
COMMON MTOT.NP.IA(l)
INCLUDE IOCOM.INC
INCLUDE CNTCOM06.INC
LOGICAL ERR
EXTERNAL WEDATO
WRITE(NTW,2000)
WRITE(NOT,2000)
2000 FORMAT!/, • — El*n»nt H»l Flow Rat*!-)
CALL DATGENIWEDATO.NrELM.ERR)
IF(ERR) RETURN
CALL REPRTE(IA(MPWE),HFELM)
RETURN
END
SUBROUTINE HEOATO (N. ERR)
C—SUB t NEDATO - CALLS KEOAT1 PASSING ARRAYS
C
COMMON MTOT.NP.LMl)
INCLUDE CNTCOMaC.INC
LOGICAL ERR
"" IIEDAn(IA,IA
-------
NBS: Indoor Air Quality Model
Phase II Report
Appendix - FORTRAN 77 Source Code
LOGICAL ERR
CHARACTER ENOFLAC'3
WRITE(NOT. 2500)
2500 FORMAT I/.' •— Nominal Tin. Conitanti')
ERR - .FALSE.
CALL REPRTT(Lk(HPP),IA(MPV),NFEON,l)
C—0.0 WRITE BEAOER AND READ PRECISION
C
WRITE (NOT. 2000)
WRITE (NTH, 2000)
2000 FORMAT!/,
+ • — TIMECONSl TIME CONSTANTS - CONTAMINANT DISPERSAL SYSTEM <)
EP1 - E?
CALL TOER('B',BP1,1)
WRITE(NOT, 2010) EF1
WRITE (NTW, 2010) EP1
2010 FORMAT)/' Cojivarganca pazan»taz, apiilon, ...', 010.3)
C
C—1.0 CHECK IF FLOW SYSTEM AND ELEMENT DATA ARE DEFINED
C
IF(NFEQN.EO.O) TEEN
WRITE(NTW.2100)
WRITE(NOT, 21001
2100 FORMAT (
+ • •••• ERROR: Nunbar of Clef iy>tam DOFi - O.',/,
+ ' FLOWSYS eoramnd taumt ba axacutad.')
RETURN
ELSEIF(NFELM.EO.O) TEEN
WRITE (NTW, 2110)
WRITE(NOT,2110)
2110 FORMAT!
* • ••" ERROR: Nunbar of fKM floo alamntl . O.',/,
* ' FLOVELEM conftand nuct b* axaoutad.')
RETURN
ENDIF
C
C—2.0 DEFINE AND INITIALIZE ARRAYS
C
CALL DEUTECWE •)
CALL DELETE I'V •)
CALL DELETECF •)
CALL DEFINRCF •.MPF.NFEON.NFEaN)
CALL DEFINRCV '.MPV.NFEON,!)
CALL DEFnairWE ',MPWE,NFELN,1I
CALL ZEROR(IA(MPV),NFEON,1)
CALL ZEROR(IA(MPWE),NFELM,1)
C
C—3.0 GET ELEMENT FLOW RATES (WE)
C
CALL READWE(ERR)
IF (ERR) TEEN
CALL ABORT
RETURN
ENDIF
C
C—4.0 FORM [F] (ALLOW "END- BEFORE VOL. MASS DATA TO JUST FORM [F])
C
OPEN(NDl,FILE-OLD<,
•FORM-' UNFORMATTED')
REWIND ND1
CALL FORMF(IA(MPKEQ),IA(MPF),IA)MPNE), 'FULL',ERR)
IF(ERR) TEEN
CALL ABORT
RETURN
ENDIF
CLOSE(ND1)
CALL FREXCC '.ENDFLAG, 3.1)
IF)ENDFLAC.EO..'END'> RETURN
C
C—5.0 GET NODAL VOLUMETRIC MASS AND REPORT NOMINAL TIME CONSTANTS
C
CALL READV(ERR)
IF(ERR) THEN
CALL ABORT
RETURN
ENDIF
IF(ECHO) WRITE(NTH,2500)
C—(.0 PREKULTIPLY [F] BY [V] INVERSE
C
CALL VINVF(IA(MPF),IA(MPV>,NPEQN,EP,ERR)
IF(ERR) TSIN
CALL ABORT
RETURN
ENDIF
C •
C—1.0 SOLVE EIGENVALUE PROBLEM
C
IF(ECHO) WRITE(NTW, 2700)
WRITE (NOT,2700)
2700 FORMAT)/, ' — Actual Tin. Coflltantl •)
WRITE (NTW, 2710)
2710 FORMAT!/, • — NOTEt Confutation of actual tin con.tanti •,
+ 'aay taka eofiaidarabla time.')
NIT - 0
CALL EIGEN2(U(MPF),IA(M?F>,NFEON,NIT,EP1)
C
C—a.O REPORT TIME CONSTANTS ( ITERATION INFORMATION
C
CALL REPRTT(Uk(MPF),IA(MPV),NFEQN,2)
WRITE(NTW,26()0) ABS(NIT)
WRITE (NOT. MOO) ABS(NIT)
2*00 FORMAT!/• Nunfaar of it.r.tion. u>ad ...-,15)
IF((NIT.LT.O).OR.(NAIT.EQ.50» THEN
WRITE (NTW, :i>10)
WRITE (NOT, :|610)
2110 FORMAT I' •••• WARNING. Procaaduia did net conv.rg..')
ENDIF
C
C—9.0 DELETE ARRAYS
C
"M- DELETECWI •)
CALL DELETECV ')
CALL DSLETECF •)
RETURN
END
SUBROUTINE VTNVF(F,V,NFEgN,EP,ERR)
C—SUBt VTNVFi BVJiUATBS [V)-1[F] I CALLED BY TIMCON
INCLUDE IOCOM.INC
REAL«I F(NFiaN.l), V(NFEON), EP, EPZERO
LOGICAL ERR
C
C—1.0 FIND MAX VOLUMETRIC MASS TO ESTABLISH RELATIVE MACHINE ZERO
C
VMAX - O.ODO
DO 10 I-l.Nl'EON
IF(Vd).GT.VMAX) VMAX-V(I)
10 CONTINUE
EPtERO - EP'VMAX
C
C—2.0 EVALUATE PKOOUCT [V]-1(F]I ERR IF DIV BY MACHINE ZERO
C
DO 20 I-l.NFEQN
VII - V(I)
IF(VII.LE.EPtERO) THEN
WRITE (NTH, 2000) I
WRITE (NOT, 2000) I
ERR - .TRUE.
RETURN
ENDIF
2000 FORMAT) .
+ • •••• BRRORi Volumtrie naa» !••• than rvlativ* naehina taco.',/,
+' Equation nunfeari ',15)
DO 20 J-l.NFEQN
F(l.J) - F(I,J)/VII
20 CONTINUE
Append-9
-------
NBS: Indoor Air Quality Model
Phase II Report
Appendix - FORTRAN 77 Source Code
RETURN
END
SUBROUTINE REPRTT(F,V,NFEON,OPT)
C— SUBjRSPRTT - RETORTS TIME CONSTANTS! CALLED BY TIHCON
C
C— 1
C
INCLUDE IOCOM.INC
REAL'S F(NFEQN.l).V(NFEaN>
INTEGER OPT
IF(OPT.EQ.l) THEN
.0 RETORT NOMINAL TIME CONSTANTS V(I,I)/F(I,I>
nun (NOT, 2010)
Tviwnni HDTTVfimf 9nin\
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
r-
TIMEI3)
KPKE WEINFELM)
DAT(l) 1 •
1 1
1 •-
1 1
1 *-
1 1
DAT(J) 1- - •-
TM(2) TM(l)
KHDAT TDATI2)
KPNDAT WDAT(NFELM,2)
: START TIKE, ENDTIME, TIMESTEP
: CURRENT ELEMENT MASS FLOW VALUES
Tina hiatoriai of axeitation data ara
dafinad •• atap-viaa function* of tin»
aaino; arbitrary valuat or, optionally,
ganaratad intarnadiata valuai of
aqual atap aiaa.
i CURRENT ARBITRAItT TIKE VALUES
i CORRESPONDING ELEM. FLO* DATA
HRITE (NOT. 2020) (N. V(N)/F(N.N), N-l.NFEON)
IF(ECBO) WRITE (Km, 20201 (N, V(N)/F(N,N), N-l.NFEON)
COMMON man i MPTDAT.KPNDAT
REAL»I TIKE(31
LOGICAL ERR
CHARACTER ENDFLAG*3
C—2.0 REPORT ACTUAL TIME CONSTANTS
C
WRITE(NOT.2040)
IF(ECEO) NRITE(NT»,2040)
WRITE(NOT,2020) (N, 1.000/FIN.N), N-l.NPEON)
IF(ECEO) mm (NT*. 2020) (N, 1.0DO/F(N,N), N-l.NFEON)
END IF
2010 FORMAT!/,6X,4(2X,'Noda Valua',]X»
2020 P010IAT«tX..IIDT
OFTIONALLT EQUAL STEP TIME HISTORIES MAY BE GENERATED
.' FLOHDAT [T»nl.n2,n3J Ganarata alanant flow tina historias.', /,
. • TIKE-nl nl - tina',/,
.• nl.nl,nJ N-n4 nl.n2.n3 - noda: flrit, la«t, incr.',/,
.' ... B4 • alanant naaa flow rata.',/,
.' i1,/,
.' END1,//,
IMPLICIT REAL-e(A-H.O-I)
• CAL-SAPt DATA i COMMON STORAGE
COMMON MTOT.NP,IA(1)
INCLUDE IOCOM.INC
INCLUDE CNTCOmC.INC
• FLOOAT) DATA I COMMON STORAGE
— DICTIONA RY OP VARIABLES
POINTER VARIABLE DESCRIPTION
ERR - .FALSE.
WRITE(NOT,2000)
NRITE (DTK, 2000)
2000 FORMAT!/, ' — FLONDATt ELEMENT FLOW TIME HISTORY DATA')
C
C—1.0 CHECK TO SEE IF ELEMENTS HAVE BEEN DEFINED
C
IF(NFELM.EO.O) THEN
•RITE(NTH, 2100)
WRITE (NOT, 2100)
2100 FORMAT)
» • •••• ERRORt Nuntoar of flo* alanantl - O.'./,
» ' FLOWELEN eooHnd wit ba axacutad.'l
CALL ABORT
RETURN
EKDIF
C
C—2.0 GET DATA GENERATION CONTROL DATA
C
TIMS(l) - O.ODO
TIME(2) - O.ODO
TIME(3) - O.ODO
CALL FREERCr,TIME(l),3)
IF(TIMEO).LT.O.ODO) TBEN
WRITE (NTN, 2200)
•RITE(NOT,2200)
2200 FORMAT|' •••• ERROR) Tina (tap nay not ba nagativa.')
"". ABORT
RETURN
ELSEIF(TIMEO).GT.O.ODO) THEN
IF(TIMEI2).LT.TIME(1)) THEN
WRITE(NTH,2210)
WRITE (NOT, 2210)
2210 FORMAT)
+' •••• ERROR: Final tina mitt ba graatar than initial tima.')
CALL ABORT
RETURN
ENDIF
IF (ECHO) WRITE (NTN, 2220)
WRITE (NOT, 2220)
2220 FORMAT!/, ' — Oanazation Control Vaziablai')
IFIECHO) WRITB.mrr
C
CALL NOPEN(ND1, (FNAMB(liLFNAME)//' .WOT'), 'UNFORMATTED')
C
C—4.0 READ ( GENERATE FLOW DATA
C
Append-10
-------
NBS: Indoor Air Quality Model
Phase II Report
Appendix - FORTRAN 77 Source Code
IRITE (NOT, 2400}
mtlTE (HTM, 2400)
2400 FORMATI/, ' — Slmuit HIII Flo. Tin Biitory Data')
c
C 4.1 DEFINE 1 INITIALISE ARRAYS
C
ERR - .THUS.
RETURN
ENDIF
>"BT-T- GETDDTIHDAT.ERR)
IF (ERR) RETURN
CALL DELETE ('TDAT')
CALL DELETE I'NDAT')
NXDAT - 1
iriTIME(l).CT.O.ODO) TEEN
CALL DELETE CUE •)
CALL DEFINRCNE •.MPME.NFELM, 1)
CALL tEROR(IA(MP«E),NFELM.l)
KNDAT - 2
ENDIF
CALL DETINRt'WDAT'.MPTOAT.NFELM.NHDAT)
CALL DEFINRrTDAT',MPTDAT,l,2>
CALL ZEROR(IA(HP*DAT),NFELM,NK)AT)
CALL ZEROR(LMMPTDAT>,1.2>
.2 CENERATE VALUES i RRITE TO .*DT
IP(TIME(3).CT.O.ODO) THEN
>•"' CENKDl(LMNPIIX),IA(MPTDAT),IA(MniDAT),TIME,ERR)
17(ERR) TEKN
CALL ABORT
RETURN
END IT
ELSE
CALL CENND2(IA(MPTDAT),IA(MPNDAT).ERR)
I?(ERR) TEEN
CALL ABORT
RETURN
ENDIF
ENDIF
C—5.0 DELETE ARRAYS, CLOSE ELEMENT FLO* DATA FILE, SKIP TO -END'
C
CALL DELETE('TDAT')
CALL DELETE CWDAT-)
CALL DELETE CUE •)
CLOSE(ND1)
IFIMODE.EO.•BATCB-) TEEN
900 IFIEOC) RETURN
CALL FREE
SO TO 300
ENOIF
RETURN
END
CALL GETTDT(TDAT)
IFIEOC) TEEN
•KITE(NTH. 2100)
•RITE(NOT.2100)
ERR - .TRUE.
RETURN
ELSEIF(TDATC1).LT.TDAT(2» TEEN
WRITE (NTH, 2110)
WRITE (NOT, 2110)
2110 FORMAT)' "••• ERROR. Tin data out of ••q
ERR - .TRUE.
RETURN
ENDIF
CALL EETNDT(WDAT,ERR)
IF(ERR) RETURN
C
C—2.0 GENERATION TIME LOOP
C
DO 100 T-TIMI(1).TINE(2),TIME(I>
C
C 2.1 UPDATE EXCITATION FUNCTION DATA IF NEEDED
C
20 IF(T.CT.TDAT(1» THEN
CALL CETTDT(TDAT)
IFIEOC) TBEN
«RITE(Nn,2100)
•RITE(NOT,2100)
ERR - .TRUE.
RETURN
EMEIF|TDAT(1).LT.TDAT(2I) TEEN
•RITE(HIV,2110)
mun (NOT, 2110)
ERR - .TRUE.
RETURN
ENDIF
CALL CETHIT (NDAT, ERR)
IF (ERR) RinURN
CO TO 20
ENDIF
C
C 2.2 COMPUTE INTERPOLATION FRACTION
C
XT- (T-TDAT(2))/(TDAT(1)-TDAT(2))
C
C 2.] COMPUTE I«E(T))
C
CALL IEROR(»E,NFELN,1>
SUBROUTINE GEN*D1(KE,TDAT,NDAT,TIME,ERR)
C—SUB: CENND1 - GENERATES ELEMENT MASS FLOW DATA, AT EQUAL TIME STEP
C INTERVALS, FROM GIVEN ARBITRARY DISCRETE TIME DATA
c
IMPLICIT REAL'S(A-fl.O-I)
00 2> N-l.NFELM
VEIN) - KDATIN. 2) * XT*(NDATIN.1)-KDATIN,J))
2} CONTINUE
C
C 2.4 NUTE TIKE, COMMON STORAGE
C
COMMON /rum/ MPTDAT, MPNDAT
LOGICAL ERR
c
C CENHDlt DATA > COMMON STORAGE
C
.REAL'S •E(NFELM),TDAT(2),IIDAT(NFELM,2),TIMEO)
C
C—1.0 GET FIRST TKO TIME HISTORY RECORDS I TDAT(l), NDATINFELM, 1) )
C
CALL GETTDT(TDAT)
IFIEOC) TEEN
•RITE(NTN,2100)
•RITE(NOT,2100)
2100 FORMAT!' •••• ERRORt In«»fIici«nt data.')
200 CONTINUE
C
C—3.0 NRITE ONE ADDITIONAL TIME VALUE TO DISK
C
•RITE (ND1) T
RETURN
END
SUBROUTINE GETTDT(TDAT)
C—SUB: CETTDO - UPDATES TIME DATA VALUES
C
INCLUDE IOCOM.INC
Append-11
-------
NBS: Indoor Air Quality Model
Phase II Report
Appendix - FORTRAN 77 Source Code
REAL'S TDATI2)
C
C—1.0 UPDATE OLD VALUES
C
THAT 12) - TDAT(l)
C
C—2.0 READ HEN VALUE
C
C CHECK FOR END-OP-COMMAND -DID-
ir(EOC) THEN
EOD • .TRUE.
RETURN
END IF
IF (MOOS. EO.'INTER') CALL PROMPT <' TIME>')
CALL FREE
IFIMOOB.EQ. 'BATCH') CALL FREENR(NTN)
C CHECK FOR END-OP-COMMAND -END"
IF(EOC) THEN
EDO - .TRUE.
RETURN
ENDIF
CALL FREERCS'.TDAT(l).l)
C REPORT
IP(ECBO) HRITE(NTH. 2020) TO»T(1)
WRITE (NOT, 2020) TDAT(l)
2020 FORMAT!/, • — Tin»t ',010.3)
RETURN
END
SUBROUTINE CETNDT(NDAT.ERR)
-SUB. CETNDT - UPDATES ELEMENT FLON DATA VALUES
INCLUDE CNTCOmi.INC
RETURN
END
SUBROUTINE CENHD2 (TOAT, NDAT, ERR)
C—SUB I CENMD2 - GENERATES ELEMENT MASS FLON DATA, AT SIVEN TIKE STEP
C INTERVALS, FROM GIVEN DISCRETE TIME DATA
C
IMPLICIT REAL»I(A-B,0-I)
INCLUDE IOCOM.INC
INCLUDE CNTCOMM.INC
• FLOVDATl DATA I COMMON STORAGE
COMMON /FLODT/ MPTDAT, MPVDAT
LOGICAL ERR
EXTERNAL NDATO
• GEWD2t DATA t COMMON STORAGE
REAL-8 TDAT(2),M)AT(NrELM.l)
C—1.0 GET FIRST TIME BISTORT RECORD ( TDAT(l), KDATINFELM, 1) )
C
CALL CETTDT(TDAT)
IF(IOC) RETURN
TOAT(J) - TDAT(l)
CALL DATSEN|VDATO,NFELM,ERR)
IF (ERR) RETURN
CALL RKPKTE(IDAT(1,1),NFELM>
HUTE (NOD TDAT(l)
WRITE(ND1) (KDAT(I,1).I-1,NFELM)
LOGICAL ERR
REAL'S «DAT|NFELM,2)
EXTERNAL ODATO
C
C—1.0 UPDATE 'OLD' DATA VALUES; INITIALItS 'HEN' DATA VALUES
C
DO 10 N-1,NFXU<
KDAT(N,2) - NDAT(N.l)
10 WHAT(N.I) - O.ODO
C
C—2.0 READ NEN VALUES
C
CALL DATCEN(«DATO,NFELM.ERR>
IF(ERR) RETURN
CALL RXPRTE(WDAT(1,1>,N7ELM)
RETURN
END
SUBRODTIHI «DATO(N, IRK)
-SUBtNDATO - rAT,T,» KDAT11 PASSING ARRAYS
COMMON MTOT,NP,IA(1)
INCLUDE CNTCOMaC.INC
COMMON /FLODT/ MPTDAT, KPNDAT
LOGICAL ERR
CALL NDATHIA(MPNDAT).NTELM.N>
RETURN
END
REAL*« «DAT(NTELM,1>
CALL 7REER('V,NDAT(N,1),1)
C—2.0 GET ADDITIONAL TIME HISTORY RECORDS
C
20 CALL EETTDT(TDAT)
IFIEOC) CO TO 300
IP(TDAT|1>.LT.TDAT(2» THEN
•RITE(NTH.J100)
WRITE(NOT.2100)
aiOO roRMATC **** ERRORl Tine data out of ••qu«ne*.')
xm - .TRUE.
RETURN
ENDIF
TDAT(2) - TDAT(l)
CALL DATGEN(KDATO,NFELM,ERR)
IF (ERR) RETURN
CALL REPRTE(«DAT(1,1),NFELM)
VRITE(NDl) TDAT(l)
WRIT! (HDD (NDAT(I,1),I-1,NFELM)
GO TO 20
C
C—3.0 KRTTE ONE ADDITIONAL TIME VALUE TO DISK
-KDATO C
300 KRITI(NDl) TDAT(l)
RETURN
END
SUBROUTINE HDAT1 (HDAT, NFELM, N)
-SUBtNDATl - READS ELEMENT MASS FLO* RATE TIME HISTORY DATA
SUBROUTINE EXCDAT
-SUB:EXCDAT - COMMAND TO READ EXCITATION DATA t GENERATE STEPHISE
TIKE HISTORIES OF EXCITATION VALUES, E(NfEQN),AND
WRITES TIME HISTORIES IN FORMAL-
TIME
(E(I),I-1,NFEQN)
TIME
(E(I).I-l.NFEON)
TO FILE .EDT
.' EXCITDAT IT-nl,n2,n)) C«r..r«t. ..citation tin* hi.tori...',/,
.' TIME-nl nl - ti«»',/,
.' nl,n2,nl CC-n< nl,n2,n] - nodai tint, lait, Incr.1,/,
Append-12
-------
NBS: Indoor Air Quality Model
Phase II Report
Appendix - FORTRAN 77 Source Code
C .'....
C .' :',/.
C .' END',//,
C
n4 • •xeitationt oone. or gmn. rat*.',/, ENDIF
IF(ECBO) WRITE (NTW, 22201
WRITE (NOT, 1220)
IMPLICIT REAL*>(A-B,0-Z)
COMMON moT.NP.IA(l)
INCLUDE IOCOM.INC
INCLUDE CNTCOM88.INC
I bAlUATI UAJA t Wnnwt slVMAbA — —
C
C
C POINTER VARIABLE DESCRIPTION
C
C TIME)}) : START TIKE, ENDTIME, TIKESTEP
C KPE E(N?ELM) t CURRENT EXCITATION VALUES
C
C TIME
C
C DATI1)
C
C
C
C
C
C DATI2)
C
C
BISTORT DATA
• Tin* histaria* of axeitation data ar*
1 oafinvd aa atap-vi** function* of tin*
•- ncing
arbitrary valua* or, optionally.
1 gvnaratad int*rm*diat* valu** of
*- *qnal
1 ,
- - •-
TMI2) TH(1>
•tap tita.
C MPTOAT TDATI2) t CURRENT ARBITRARY TIME VALUES
C~ MPEDAT EDAT(NFNOD,2> : CORRESPONDING EXCITATION DATA
COMMON /EXCDT/ MPTDAT,MPEOAT
REAL'S TIME (3)
CHARACTER ENDFIAC']
LOGICAL ERR
ERR - .FALSE.
«RITE (NOT, 2000)
WRITE (NT», 20001
2000 FORMAT!/, • — EXCITDAT: EXCITATION TIME BISTORY DATA')
C—1.0 CBZCX TO SEE IF FLOW SYSTEM HAS BEEN DEFINED
C
IF(NFEOX.EQ.O) THEN
WRITE(NTH,2100)
WRITE (NOT, 2100)
2100 FORMAT(
* • •••• ERROR. Hunter of fie* TEEN
WRITE(NTW,2210)
WRITE(NOT, 2210)
2210 FORMAT(
»' ••" ERROR! Final tina> unit b* grcatar than initial tin*.')
CALL ABORT
RETURN
IF(ECBO) WRITE (NTH, 2230) (TIKE(I), 1-1, 3)
WRITE(NOT,22JO) (TIME(I), 1-1, 3)
2230 FORMAT!/,
Initial tina ',C10.3,/,
Final tima '.G10.3,/,
. • Tin* itap inoranint ',C10.3)
ENDIF
C
C—3.0 OPEN
2400 FORMAT!/,' — Nodal E«oltation Tin* Biitozy Data')
C
1.1 DEFINE t INITIALIZE ARRAYS
CALL DELETE I'THAT')
CALL DELETE CEDAT')
CALL DELETE!'E ')
CALL DEFINRCE '.MPE.KFEOjN.l)
CALL .EDT
IP(TIME(3).«T.O.ODO) TBIN
CALL CENtDl (lA(MPKEO), IA(K?E), lA(MPTDAT). IAIMPEDAT) .TIME, ERR)
IF (ERR) TI1KN
CALL ABORT
RETURN
ENDIF
ELSE
CALL GENEI)2(IA(MPKEO),IA(MPE),IA(MPTDAT),IA(MPEDAT),ERR)
IF (ERR) TI3EN
CALL AB3RT
RETURN
ENDIF
ENDIF
C—5.0 DELETE ARRAYS, CLOSE ELEMENT FLOW DATA FILE, SKIP TO "END"
C
"M- DELETE CTDAT1)
CALL DELETE CEDAT'l
CALL DELETE I'E ')
IF(MODE.EO.'BATCH') TBEN
500 IF(EOC) fITURN
CALL FREE
CO TO 900
ENDIF
RETURN
END
SUBROUTINE GENEDKKEO.E.TDAT.EDAT, TIME, ERR)
C—SUB I CENXD1 - GENERATES EXCITATION DATA, AT EQUAL TIME STEP
C INTERVALS, FROM GIVEN ARBITRARY TIME DATA
C
IMPLICIT REAL't(A-B, C-Z)
INCLUDE IOCOM.INC
INCLUDE CNTCOH8C.INC
Append-13
-------
NBS: Indoor Air Quality Model
Phase II Report
Appendix - FORTRAN 77 Source Code
LOGICAL ERR
c
C GENEDli DATA > COMMON STORAGE
C
RIAL'S E
INTEGER XEO(NFNOD)
C—1.0 GET FIRST TMO TIKE BISTORT RECORDS ( TDAT(2), EDATIN7NOO, 2) )
C
CALL GETTDT(TOAT)
IF(gOC) THEN
•RITE(NT», 2100)
BUTE (NOT, 2100)
2100 FORMATC •••• ERRORi In.uftici.nt dau.M
ERR - .TRUE.
RETURN
ENDIF
CALL CETEDTIEDAT.ERX)
IF(ERR) RETURN
CALL CETTDHTOAT)
IFIEOC) THIN
nUTEINTR.2100)
•RITE (NOT. 2100)
ERR - .TRUE.
RETURN
ELSEIF(TDAT(1).LT.TDAT(2)) THEN
•RITE(NTN,2110)
•RITE(NOT,2110)
2110 FORMAT!' •••• ERROR> Tim data oat of ••queiic*. •)
ERR - .TRUE.
RETURN
ENDIF
CALL GXTEDT (EDAT, ERR)
IF(ERR) RETURN
200 CONTINUE
C
C—3.0 mm ONE ADDITIONAL TIME VALUE TO DISK
C
•RmiNDl) T
RETURN
END
SUBROUTINE GETEDT(EDAT.ERK)
C—SUB I GETEDT - UPDATES EXCITATION DATA VALUES
C
COttON MTOT,NP,IA(1)
INCLUDE IOCOH.INC
INCLUDE CNTCOMM.INC
C
C— CETEDTi DATA t COMMON STORAGE
C
LOGICAL ERR
REALM SDAT(NmOD.2>
EXTERNAL EDATO
C
C—1.0 UPDATE 'OLD' DATA VALUES; INITIALIZE 'NEW' DATA VALUES
C
DO 10 N-l.NFHOO
EDAT(N,2) - EDAT(N.l)
10 EDAT(N.l) - O.ODO
C
C—2.0 READ NEK VALUES
C
CALL DATGEN (EDATO, NFNOD, ERR)
IF(ERR) RETURN
C—2.0 GENERATION TIME LOOP
C
DO 200 T-TIME(1),TIME<2>,TIME(J>
C
C 2.1 UPDATE EXCITATION FUNCTION DATA IF NEEDED
C
20 IFIT.GT.TDAT(l)) TEEN
CALL CETTDT(TOAT)
IFIEOC) TEEN
•RITE(NTH,1100)
•RITE(NOT,2100)
ERR - .TRUE.
RETURN
EUEIF(TDAT(1).LT.TDAT(2» THIN
«RITE(NTN,2110)
•RITE (NOT, 2110)
ERR - .TRUE.
RETURN
ENDIF
CALL GETEDT (EDAT, ERR)
IF(ERR) RETURN
GO TO 20
ENDIF
-2.2 COMPUTE INTERPOLATION FRACTION
XT-
DO 21 N-l,NFNOD
NEO - ABS(KEO(N»
IF(NEO.NE.O) E(NEO) - EDATIN.2) + XT* (EOATIN, 1)-EDAT(N,2) )
23 CONTINUE
C
C 2.4 WRITS TIKE, (E(T)l TO ND1
C
«*".'. REPRTN(EDAT(1,1),IA(MPKEQ),NFNOD,NFNOD)
RETURN
SUBROUTINE EDATO (N, ERR)
C—SUBiEDATO - CALLS EDAT1 PASSING ARRAYS
C
COMMON HTOT,NP,IA(1)
INCLUDE CNTCOMM.INC
COMMON /EXCDT/ MPTDAT.KPEDAT
LOGICAL EM
<*".'- EDAT1(IA(KPEDAT),NFHOD,N)
RETURN
END
SUBROUTINE EDAT1 (EDAT. NFNOO, N)
C—SUBlEOATO - READS EXCITATION TIME HISTORY DATA
C
REAL'S EDAT(NFNOD,1)
CALL FREERCG',EDAT(N,1),1)
RETURN
END
SUBROUTINE GENED2(KEQ,E,TDAT. EDAT. ERR)
C—SUB I GENED2 - GENERATES EXCITATION DATA FROM CIVEN TIME DATA
C
IMPLICIT REAL*KA-H,0-t>
COMMON MTOT.NP.IA(l)
INCLUDE IOCOM.INC
INCLUDE CNTCOMaC.INC
VRITE(NDl) T
•RITE(NDl) (E(I),I-1,NFEON)
LOGICAL EM
EXTERNAL EDATO
Append-14
-------
NBS: Indoor Air Quality Model
Phase II Report
Appendix - FORTRAN 77 Source Code
CENED2I DATA 4 COMMON STORAGE
REAL'S TDAT(2),EDAT(NFKOD,1), E(NFEQN)
INTECER 1
C—1.0 GET FIRST TIME HISTORY RECORD ( TDAT(l). EDAT(NFNOD, 1) )
C
CALL CETTDTITDAT)
IF(EOC) RETURN
TDATI2) - TOAT(l)
CALL DATGENIEDATO.NFNOD.ERR)
IF(ERR) RETURN
DO 10 N-l.NFNOD
NEO - ABS(KEO(N»
10 IF(NEQ.NE.O) E(NEQ) - EDAT(N.l)
CALL REPRTN(E,IA(MPKEO).NFEQN,NTNOD)
WRITE (NOD TDAT(l)
WRITE(NDl) (E(I).1-1.NFEQN)
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c —
VARIABLE
TIME!))
TWDAT
TEDAT
PINT
P SCALE
POINT
START TIME, END TIME, TIMS INCREMENT
TIME OF NEXT ELEMENT FLOW RATE RECORD
TIME OF NEXT EXCITATION RECORD
RESPONSE RESULTS PRINT INTERVAL
RESULTS PLOT FILE SCALE FACTOR
ERS TO BLANK COMMON LOCATIONS
MPFS FS (NFEQN. 2«MFBAN-1) I [F>] DYNAM ALG. MATRIX (ASYM-COHPACT)
NPC
MPCD
MPCDD
MPG
C(NFEQN) i CURRENT (C)
CD(NFEQN) i CURRENT d|C)/dt
CDD(NFEQN) ' CURRENT d/dt ldlC)/dt|
G (NFEQN) I CURRENT |G)
ERR - .FALSE.
C—2.0 GET ADDITIONAL TIKE HISTORY RECORDS
C
20 CALL CERDT(TDAT)
IF(EOC) CO TO 300
IF(TDATU).LT.TDATm) TEEN
WRITE(NTH,2100)
WRITE(NOT,2100)
2100 FORMAT!1 •••• ERROR) Tin data out of iwiuanca.')
ERR - .TRUE.
RETURN
ENDIF
TDAT(J) - TDAT(l)
CALL DATGENIEDATO.NFNOO.nUl)
DO 22 N-l.NTNOD
NEO - ABSIKEOIN))
22 IF(NEQ.NE.O) E(NEQ) • gOAT(N.l)
CALL REPRTN(E.IA(KPKIQ>.NFION.NTNOD>
•RITE (HDD TDAT(l)
WRITS(ND1) (E(I).I-1,OTEQN)
GO TO 20
C
C—j.O VRITE ONE ADDITIONAL TIKE VALUE TO DISK
C
300 VRITE (HDD TOAT(l)
MRITE(NOT,2000)
MRITE(NTH,2000)
2000 FORMAT!/, • -— DYNAMIC I DYNAMIC SOLUTION')
C—1.0 CHECK IF SYSTEM, ELEMENT, AND EXCITATION DATA ARE DEFINED I AVAIL
C
IFtNFEON.EQ.'JI THEN
WRITE(NTW,1100)
WRITE (NOT,1100)
2100 FORMAT!
••• ERRORi Hunter of Iloo lyitra DOFi - O.'./.
+ • FLOWSYS coraBud nu.t ba axaeutad.')
CALL ABORT
RETURN
ELSEIFINFELH.EO.O) THEN
WRITE (NTH, 2110)
WRITE (NOT. 2110)
2110 FORMAT!
+ ' •••• ERUORt NuntMr of flow alanntl - O.',/,
+ ' FLOWELEM eonnuid n)it b* •x«cuc*d.t)
CALL ABOR1'
RETURN
ENDIF
RETURN
END
SUBROUTINE DYNAM
C—SUBtDYNAM - COKNAND TO FORM 1 SOLVE DYNAMIC PROBLEM
C
C [F(t)](C| t [V)d(C)/dt - (C|t)l
c
C • EXCITATION, (C) AND PRESCRIBED (C). UPDATED AT DISCRETE
C TIMES USED TO DEFINE EXCITATION (READ FROM HDD
C • FLOW MATRIX, [Fl. UPDATED AT DISCRETE TIKES USED TO
C DEFINED ELEMENT FLOW RATES (READ FROM ND2)
C
c
c
c
c
c
c —
. ' DYNAMIC
. ' T-nl.n3.n3 A-n«
.' n9,n(,n1 IC-ne
.' ...
.' i1,/.
.' END1,//.
Dynanie •olution. ' ,/,
nl,n2,n3 • init. final, iner; n4 •alpha',/.
n5,n6,nT - nod«! fint, la«t. Inez.'./,
n8 » nodal initial concentration*'./.
IKPLICIT REAL'S(A-fl.O-Z)
COMMON MTOT.NP.IA(l)
INCLUDE lOCON.mC
INCLUDE CNTCOM86.INC
COMMON /OYNM/ TWDAT.TEOAT
LOGICAL ERR, FOUND
REAL'S TIKE(3), PSCALE
INTECER PINT
INQUIRE (FILIXFNAKE (11 LFNAKt)//'. F8L'), EXIST-FOUND)
IFI.NOT.FOUHD) THEN
WRITE(NTW, 2120) (FNAME(ltLFNAME) //' .FEL')
WRITS(NOT, 2120) (FNAME(llLFNAME)//'.FEL')
2120 FORMAT!' "••• ERRORI Elanant data flla '.A,' not found.',/,
+ ' FLOWELEM eoanand rauit b* axacutad.')
CALL ABORT
RETURN
tNDIF
INQUIRE (FIU1- (FNAME (1 tLFNAME) //' .WDT'). EXIST-FOUND)
IF(.NOT.rOUllD) TEEN
WRITE(NTW,2130) (FNAMEIIiLFNAME)//'.WDT')
WRTTE(NOT.2130) IFNAMEIItLFNAME)//' .WDT')
2130 FORMATI• ••" ERROBi Flov data fil» '.A,' not found.',/,
+ ' FLOWDAT eonnBnd wit b* axaeutvd.')
CALL ABORT
RETURN
ENDIF
INQUIRE (Fill-(FNAKE (1 tLFNAME)//'.EDT') .EXIST-FOUND)
IF (.NOT.FOUND) THEN
WRITE (NTH', 2140) (FNAME (1:LTNAHEI//'.EOT'I
KRITEINOT,2HOI (FNAHE(liLFNAKE)//' .EOT')
2140 FORMAT!' •••• ERRORI Excitation data film ',A,' not found.',/,
* ' BXCITDAT connand njct b« axacutad.')
CALL ABOUT
RETURN
ENDIF
C
C—2.0 GET DYNAMIC SOLUTION CONTROL VARIABLES
C
WRITE (NTW, IIJOO)
WRITE (NOT. :1JOO)
Append-15
-------
NBS: Indoor Air Quality Model
Phase II Report
Appendix - FORTRAN 77 Source Code
2900 FORMATS1 — Solution Control Variablai')
IP (MODE.EO.'INTER1) CALL PROMPT!' DATAV)
CALL FREE
IFIMODE.EQ. 'BATCH') CALL FREEWR(NTII)
TIHE(l) - 0.000
TIKE(2) - O.ODO
TIKE(3) - 0.000
CALL FREERCT',TIKE(11,31
IF(TIME(3).LE.0.000) TEEM
WRIT! (NTH, 2210)
HUTE (NOT, 2210)
2210 FORMAT!' •••• ERRORJ Tin itap mo.t b. graatar than 0.0.')
CALL ABORT
RETURN
ELSEIF(TIMB(2).LT.TIMEI1)) THIN
WRITE !Nn, 2220)
•RITE (NOT, 2220)
2220 FORMAT!
+ ' •••• ERROR] Final tin» mat b* graatar than initial tin.')
CALL ABORT
RETURN
ENOIF
IF(ERR) TBEN
CALL ABORT
RETURN
EKDIF
C
C—5.0 SET NODAL INITIAL CONCENTRATIONS
C
CALL READIC(ERR)
IF (ERR) TEEN
CALL ABORT
RETURN
EHOIF
C
C—(.0 OPEN ELEMENT, FLON AND EXCITATION DATA FILES, ( PLO
C
OPEN(ND1,FILE-(FNANE(liLFNAME)//'.FEL'>.STATUS-'OLD'
+FORM-'UNFORMATTED')
REWIND ND1
OPEN(ND2, FILE-(FNAKE111LPNAME) //• .NOT'), STATUS-• OLD'
+FORM-'UNFORMATTED')
REMIND ND2
READ(NDJ) TWOAT
ALPEA - 0.75DO
CALL FREER!'A',ALPEA,1)
IPHALPBA.LT.0.000).OR. IALPBA.CT.1.0DO)) TEEN
WRITE (HTM,2230)
VRITE(NOT,2230)
2230 FORMATC —• ERRORi Alpha nut to in rang* 0.0 to 1.0.')
CALL ABORT
RETURN
ENDIF
PINT - I
CALL FREEICI'.PINT.l)
IFIPINT.LT.O) TEEN
WRITE(HTM.2240)
•RITE(NOT,2240)
2240 FORMATC •••• ERROR I Raaulta print intaryal mat to > 0.')
CALL ABORT
RETURN
ENDIF
OPEN(ND3,FILE-(FNAME(ltLFNAME)//'
+1 ORM>' UNFORMATTED')
REWIND ND3
READIND3) TSDAT
IFIPSCALE.HE.O.ODO) TEEN
CALL MOPEN(RD4, (niAMEUlLFNAMS)//'
ENDIF
1.0 DEFINE ADDITIONAL SOLUTION ARRAYS
CALL DELSTECFS •)
CALL DELETE('CO •)
CALL DELETE('COD ')
CALL DEFINRCCDD •.MPCDD.NFEON.l)
CALL DEFINRCCD
CALL DEFINRCFS
EOT1),STATUS--OLD1,
>PLT'),'FORMATTED')
.KPCD.NFEQN.l)
, MPFS.NFEON. 2*MFBAN-1)
CALL IEROR(IA(MPCD),NrEON,l)
CALL IEROR(IA|MPCDD).NFEON,1)
PSCALE - O.ODO
CALL FREER!'S',PSCALE,1)
IF (ECHO) WRITE (NTN. 2230) (TIKE (II, 1-1, 3) .ALPHA, PINT
WRITE (NOT, 2250) (TIME(I),1-1.3). ALPBA.PINT
22SO FORMAT!/,
Initial tin '.C10.3./,
.' Final tin* ',510.3,/,
.' Tin* atap ineramnt '.G10.3,/,
.' Integration paramatar: alpha .... '.C10.3,/,
Raiolta print intarval Mi)
IFtPSCALE.NE.O.ODO) TEEN
IF (ECHO) WRITE(in*,2260) PSCALE
WRITE(NOT,22(0) PSCALE
ENDIF
2260 FORMAT!' Raaulta plot-tila icala factor .. ',C10.3)
C—».0 CALL PREOIC TO DO TBS WORK
C
CALL PREDICIIA(KPXEO), IA(KPF), lAIMPFS), IAIMPV),IAIMPC), IA(KPC),
•flA(HPCD), IAIKPCDD).TIME. ALPHA.NFNOO,NFEQN,MFBAN. PINT, PSCALE, ERR)
IF (ERR) CALL ABORT
C
C—10.0 DELETE UNHEEDED ARRAYS t CLOSE FILES
C
CALL DELETECFS ')
CALL DELETE ('CD •)
CALL OELETECCDD ')
CLOSE (HDD
CLOSE (NO J)
CLOSE(ND3)
CLOSE (ND4)
C—3.0 DEFINE AND INITIALIZE SYSTEM ARRAYS
C
CALL DELETE ('WE •)
CALL DELETECC •)
CALL DELETE ('C ')
CALL DELETE('F ')
CALL DELETE I'V ')
CALL DEFINRCV
CALL DEFINRCF
CALL DEFINRVS
CALL DEFINRCC
CALL DEFINRCHE
CALL ZEROR(IA(MPV),NFEON,1)
i-E" tERORIIAIMPO.NFEOM.l)
.HPV.NFEON.D
, MPF, NFEQN, 2*MFBAN-1>
,MPC,NFEON.l)
,MPC,NFEQN,1)
.MPWE.NFELM, 1)
C—4.0 CET NODAL VOLUMETRIC MASS
C
CALL READV(ERR)
C—11.0 SKIP TO END-OF-COMKAND DELIMITER 'END'
C
IF (MODE. EO.' INTER') RETURN
IF(MODE.EO.'BATCH') TEEN
1100 IFIEOC) RETURN
CALL FREE
EO TO 1100
ENDIF
END
SUBROUTINE READ1C(ERR)
C—SUBlREADIC - READS 4 REPORTS INITIAL CONCENTRATION CONDITIONS DATA
C
I MTOT.NP.LM1)
INCLUDE IOCOM.INC
Append-16
-------
NBS: Indoor Air Quality Model
Phase II Report
Appendix - FORTRAN 77 Source Code
INCLUDE CNTCOMS..INC
LOGICAL ERR
mure (NTM, 20001
MUTE (NOT, 2000)
2000 FORMAT!/, ' — Initial Condition.: Nodal Conc.ntr.tiom')
CALL DATCEN(ICDATO,NFNOD.ERR>
IF(HUI) RETURN
CALL REPRTN(IA(MPC),IA(MPKEO.),NFEON,NFNOD)
RETURN
END
SUBROUTINE ICDATO(N.ERR)
C— SUBlICDATO - CALLS ICDAT1 PASSING ARRAYS
C
COMMON MTOT.NP, IA(1)
INCLUDE CNTCOmC.INC
LOGICAL IRR
C
C
£
tnATmnTB HBBnnTnfltVfHi
C-DUMHY
C
C
C
C
C
C
C
C
C
C
C
— ICDATO C
C
C
C
C
C
C
C
C
C—
ID(NNOD)
K(NEQN,2>MBAN-1)
KS (NEON. 2*NBAN-] )
CINEON)
E(NEON)
T(NEON)
TO (NEON)
TDD (NEON)
TIME(J)
ALPBA
NNOO
NEON
MBAN
PINT
PSCALE
ERR
t EQUATION NUMBER/CODE (ORDERED BY EOTN 1)
EOTN « OF NODE N - ABS(ID(N»
ID(N| - 0 t NODE IS NOT DEFINED DOF
ID(N) - PCS t NODE IS E-PRESCRIBED DOF
ID(N) - NEE I NODE IS T-PRESCRIBED DOF
t [K) MATRIX) ASTM-BANDEO COMPACT-STORED
> (K*) - 1C] * aDT(X) MATRIX (SCALED FOR NEC ID)
t CURRENT [C] (ORDERED BY EOTN «)
1 CURRENT (El (ORDERED BY EQTN <)
> CURRENT IT) (ORDERED BY EOTN «)
> CURRENT (dT/dt) (ORDERED BY EOTN t)
1 INITIAL !d/dt ( dT/dt ) ) TO EST TIME STEP
> START TIKE, END TIME, TIME INCREMENT
t INTEGRATION PARAMETER
I NUMBER OF SYSTEM NODES
t NUMBER OF EQUATIONS
t BALF BANDWIDTH OF SYSTEM
1 OUTPUT RESULTS PRINT INTERVAL
t RESULTS PLOT-FILE SCALE FACTOR
t ERROR FLAG
CALL ICOATHIA!MPKEO),IAtHPC),NrNOD,NFEON,N,ERR>
RETURN
SUBROUTINE ICDAT1(10W,C,NFNOO,NFZQN,N.ERR>
C— SUBiICOATl - READS INITIAL CONCENTRATION CONDITIONS DATA
C
INCLUDE IOCOM.INC
INTEGER KKQ(NFNOD)
REAL'I C(NFEQN),CDAT
LOGICAL ERR
COAT - O.ODO
CALL FREER! 'C', COAT. II
IFICDAT.LT. O.ODO) TBEN
NRITE (NTH, 2000)
mm (NOT, 2000)
ERR - .TRUE.
RETURN
ENDIF
2000 FORMAT) • •••• ERROR. Nodal
ICDAT1
concentration! nay not ba negative.')
NEO - ABSIKEQ(N))
C(NEQ) - COAT
RETURN
END
PREDIC
SUBROUTINE PREDIC (ID, K.KS,C,E,T, TO, TDD, TIME, ALPBA, NNOO, NEON, MBAN,
*PINT, PSCALE, ERR)
•SUB: PREDIC - PREDICTOR-CORRECTOR 1ST O.D.E. EQUATION SOLVER
TIME STEP ESTIMATE BASED ON METHOD IN -BEAT*
BY R.L.TAYLOR - U.C. BERKELEY
SOLVES EQUATION!
[K(t)J(T| * (C)(dT/dt| - (E(t)l
NBERE; [K(t>] - STORED IN COMPACT ASYMMETRIC BANDED FORM
[C] - DIAGONAL; STORED AS VECTOR
IE It)} - EXCITATION; DEFINED PIECE-WISE LINEAR
BASED ON DIFFERENCE APPROXIMATION;
(T)n«l - |T)n * (l-a)DT(dT/dt)n * (a)DT(dT/dtln»l
NBBR2; a • "alpha*, an integration paxanatar
• 0 corresponds to Forward Difference method
• 1 corresponds to Backward Difference method
- 1/2 corresponds to Crank-Nicholson method (unstable)
DT • tiro* step increment
IMPLICIT REAL'6 (A-fl,O-Z)
INCLUDE IOCOM.INC
- PREDIC t DATA I COMMON STORAGE
RBAL-6 MNEQN,2*KBAN-l),K8 (NEON, 2>MBAN-1),C(NEQN),E (NEON), T(NEON),
•fTD (NSQH), TDD (NEON) , TIME (», ALPBA, PSCALE
INTEGER PINT, ID(NNOD)
LOGICAL ERR, TDOF, KUPDAT, EUPDAT
C—1.0 FORM INITIAL [K)
C
CALL UPDATKIK.TIMEd),KUPDAT,IRR)
IF (ERR) RETVRII
C
C—2.0 COMPOTE INITIAL TEMPERATURE BATES: (dT(0)/dt) FROM
C
C (C)|dT(0)/dt) - (E(0)| - (K)[T(OI)
C
C 2.1 GET INITIAL EXCITATION
C
CALL UPDATE(E,TIME(1).EUPDAT,ERR)
IF(ERR) RETURN
C
C 2.2 FORM RBSt (dT/dtl-0 FOR 'T'-DOF, IE)-(K] (Tl FOR 'E'-DOF
C
DO 22 I-l.NECM
C -T'-DOFi SET (oT/dt)-0
IF(TDOF(I,ID, NNOO) I THEN
C -T'-OOFi CBECX FOR dT/dt INFINITE
IF(T(I).NE.E(D) THEN
MRITEINTll, 2S20) I
NRITE(NOT, 2220) I
2220 FORMATC •"• ERROR: Can not compute for step change',
+ • in dependent variable number I',15)
ERR - .TIIUE.
RETURN
ELSE
TD(I) - O.ODO
ENDIF
C •E'-OOFi FORM (E)-(K)(TI HBERE (K) IS IN COMPACT STORAGE
ELSE
TEMP - 8(1)
Fa - MAX(l.MBAN-I-fl)
K2 - MN(3'KBAN-1,KBAN*NEQN-I)
DO 20 KX-IU.1U
J - I » KK - MBAN
TEMP - TEK> - K(I,IOX)>T(J)
20 CONTINUE
C SOLVE
- TENP/CID
Append-17
-------
NBS: Indoor Air Quality Model
Phase II Report
Appendix - FORTRAN 77 Source Code
c—j.o COMPUTE TAYLOR'S TIMESTBP CHICK
c
IP (ECHO! WRITE(NTK.J300)
WRITE(NOT.2300)
2100 PORMAT(/.' — Tim Step E>tin»te for Initi.l Condition.')
C
C I.I COMPUTE INITIAL RATE Or TEMP RATES
C FORM AMD SOLVE. (C)d(dT/dt)/dt - -(XHdT/dtl
C
00 32 I-l.NEQN
W(TDW(i,iD.NiiOD)) TEEN
TOD(I) - O.ODO
ELSE
TEMP - O.ODO
1U . MAX(l.MBAN-I-M)
K2 - MrN(2«MBAN-l,MBAN»NEOJtl-I)
00 30 KX-K1.K2
J . I * XX - MBAN
TEMP - TEMP - K(I.KX,)*TD(J)
30 CONTINUE
TDD(I) - TEKP/CII)
32 ENDIP
C
C 3.2 COMPOTE NORMS! M(T(0»|I, I HdT(O)/dt) II, I Id/dt (dT(O)/dt | M
C
TN • O.ODO
TON - O.ODO
TDDN - O.ODO
DO 34 N-l.NEON
TN - TN + TIN)"2 '
TON - TON * TD(N)"2
94 TDDN - TDDN » TDD IN) "2
TN - SORT(TN)
TON • SORT (TON)
TDDN - SORT (TDDN)
C
C 3.3 EVALUATE TAILORS EXPRESSION TOR TIME STEP ESTIMATE
C
B - O.OSDO
II-ITDON.NE.O.ODO) TEXN
DTEST - IB-TUN » SOKT(B«B«TON-TTJN * 2.000-B-TN-TDDN)) /TOON
IP (ECHO) WRITE(NTN,2320) B'lOO.ODO, OTCST.TTKtO)
NRITI(NOT.2320) B'100.000, DTEST.TIKE(3)
2320 FORMAT)/' — NOTSt Iitineted tin* itep to Unit error to',
.' tppro«.',F9.2,'% i»i',C10.3,/
. • Specified tin* lt«p i«i",S10.3>
ELSE
IP(ECBO) NRITE(NTH, 2340)
WRITS (NOT,2340)
2340 FORMAT (/' — NOTE I Unable to eltinate tin itep to limit ',
.' error for the given «y>teia.')
ENDIF
C
C—4.0 FORM AND FACTOR (X*)
C
CALL PORHKSIID.K.KS.C,ALPHA.TI«E(3),NNOD.NEON,HBAN)
CALL FACTCA (US, NEON, MBAN, ERR)
IF(ERR) RETURN
C
C—5.0 TIME STEP THRU SOLUTION
C
ADT - ALPHA-TIME (3)
DTA - (1.000 - ALPHA) *TIME(3)
ISTEP - 0
DO 500 TM-TTME(1)+TIMI(3>,TINE(I>.TIME(3)
ISTEP - ISTEP « 1
C
C 5.1 UPDATE [K], FORM AND FACTOR (K-)
C
CALL U?OATK(K,TM,XUPDAT,ERR)
IF(ERR) RETURN
IF(KUPDAT) THEN
'•"•' FORMKS(ID,X,IIS,C,ALPBA,TIMI(3),NNOD,NEON,MBAN)
CALL FACTCAdtS.NEON.MBAN.ERR)
IF(ERR) RETURN
ENDIF
C - S.2 FORM IE)
C
CALL UPDATE (I,TM,EUPOAT, ERR)
IF (ERR) RETURN
C
C - 5.3 PREDICT Tl (T) - (T) t (l-»)DT(dT/dt)
C
DO 51 N-l.NSON
IF(TDOF(N,ID.NNOD» TEEN
T(N) - BIN)
ELSE
TIN) -TIN) t DTA'TD(N)
ENDIF
51 CONTINUE
C
C - 5.4 FORM RHS:(E)-[X)(T) FOR FLUX-DOF, IdT/dt I
C
CALL RBS(ID,T,TD,E,X,KS,NNOD,NEON,MB»I)
_ C
C - 5.5 SOLVE FOR IdT/dt)
C
CALL SOLVCA(XS.TD,NEaN.NBAN,ERR)
IF (ERR) RETURN
C
C - 5.C CORRECT Ti (T) - (T) » «DT(dT/dt)
C
DO 55 N-l.NEON
IF(TDOF(N,ID,NNOD)> THEN
T(N) - I(N)
ELSE
TIN) - T(N) » ADT*TD(N)
ENDIF
55 CONTINUE
C
C - 5.1 REPORT RESULTS
C
IFIMODIISTEP.PINTI.EO.O) THEN
IF (ECHO) NRITEINTN.2510) TN
•RITE (NOT, 2510) TM
2510 FORMAT!/, • — Re.pon.. '.
CALL REFRTN(T.ID,NEaN.NNOD)
FOR TEMP-OOF
Tim '.S10.3)
TO FILE .PLT for plotting
IF(PSCALE.NE. O.ODO) TEEN
WRITE (ND4, 2530) TM, (CHARUI , Til) 'PSCAL!, 1-1, NIONI
2530 FORHATIF10.3, (10(A1,E10.4) I )
ENDIF
ENDIF
500 CONTINUE
RETURN
SUBROUTINE UPDATKIX.TM.XUPDAT.ERR)
C—SUBiUPDATX - UPDATES [X]-[F] IF ELEMENT MASS FLOW RATES CHANCE
COMMON NTOT.NP.IAIl)
INCLUDE IOCOM.INC
INCLUDE CNTCOMai.INC
COMMON /DYNM/ TUDAT.TEDAT
RIAL** K(NFEQN.**HPBAN-1), TM, TNDAT, TEDAT
LOGICAL ERR, KD7DAT
C
C—1.0 UPDATE ELEMENT FLOW RATES IF(TM.eE.TWDAT)
C
CALL UPDAT(ND2,TM,TWDAT,IA(MPNE),NFELM,KUPDAT,ERR)
IF(KUPDAT) THEN
IF (ECHO) WRITE (NTW, 2000) TM
WRITS (NOT. 2000) TM
2000 FORHAT(/, • — El*mMt Flo» lUt. UpeUt. '.30(IB-),
•f • Tim ',£10.31
CALL REPRTEIIA(HPWE).NFELM)
CALL FORMF(IA(MPKEQ),X, lA(KPWE), 'BAND', ERR)
Append-18
-------
NBS: Indoor Air Quality Model
Phase II Report
Appendix - FORTRAN 77 Source Code
ENDIF
RETURN
END
SUBROUTINE UPDATE (E.TM. EUPDAT, ERR)
.C—SUBlUPDATE - UPDATES 1EI-IO IF EXCITATION CHANCES
COMMON MTOT.NF.IA(l)
SUBROUTINE FORMXS(ID.K,KS,C.ALFBA.DT,NNOO,NEON,MBAN)
C—SUBlFORNKS - FOIIMS;
-UPDATE C (K«) - [C] * aDTIK]
C
C SCALES (*•] - IK«)*1.0D15 FOR 'T'-DOF
C
IMPtlCIT REAL«B(A-H,0-Z)
INCLUDE IOCOM.IHC
INCLUDE CNTCOMat.IKC
COMMON /DYNH/ TNDAT.TEDAT
RIAL'S E(NFEON), TK, TNPAT, TEDAT
LOGICAL ERR. EUPDAT
CALL UPDAT(ND3.TM,TEDAT,E.NFEQN,EUPDAT,ERR>
IP(EUPDAT) THEN
IF(ECEO) WRITE (DTK, 2000) TM
WRITE (NOT. 2000) TM
2000 FORMAT!/, ' — EielUtion Updat. ',31 (IB-),' Tint '.CIO.3)
CALL REPRTN(E,IA(MPKEa>.NPXQN,NrNOD>
ENDIF
RETURN
SUBROUTINE UPDAT(LUN,T,TD,D,ND,UPDATE,ERR)
C—SUB I UPDAT
C SEARCHES A SEQUENTIAL DATA RECORD, ON UNIT LUN, OF TEE FORM;
TD
. (D(I).I-l.ND)
TD
(0(11,1-1.(TO)
REAL'S K(NE01».2«MBAN-1>,KS(NEON.2-MBAN-1),C(NEON>
INTEGER ID(NNOD)
LOGICAL TDOF
ADT - ALPBA-OT
DO 10 N-l.NEQN
DO 10 H-1,2*1QAN-1
10 KS(N.H) -ADT'K(N,M>
DO 20 N-1.NE9N
20 XS(N.MBAN) - KS(N.MBAN)
C(N)
DO 30 N-l.NEON
30 ir(TDOF(N, ID.NNOOI) KS(N.MBAN) - ItS(N,MBAN)-1.0D15
TO UPDATE DATA VALUES TO CURRENT TIME. 'TV IF DATA VALUES ARE
UPDATED LOGICAL 'UPDATE- IS SET TO TRUE.
t DISCRETE TIME VALUE
I UPDATED TO NEXT VALUE
t CORRESPONDING DISCRETE DATA VALUES
TD
DID
UPDAT MUST BE •PRIMED- BY READING FIRST TD VALUE TO MEMORY
INCLUDE ICCOM.INC
REAL'S D(ND).T,TD
LOGICAL ERR. UPDATE
UPDATE - .FALSE.
10 IF(T.CE.TD) TEEN
C UPDATE DISCRETE DATA VALUES
READ(LUN, ERR-800, END-tOO) (0(I). I-l.ND)
IF(ERR) RETURN
UPDATE - .TRUE.
C GET NEXT DISCRETE TIME
READ(LUN, ERR-800, END-tOO) TD
IF(ERR) RETURN
GO TO 10
ELSE
RETURN
ENDIF
600 ERR - .TRUE.
WRITE (NTH, 1000)
WRITE(NOT,«000)
8000 FORMAT!' "" ERROR: Tin history data til* r««d error.')
RETURN
900 ERR - .TRUE.
WRITE(NTH,9000)
MUTE (NOT, 9000)
9000 FORMAT!
+ ' *••* ERJtORi EOF •ncoqntarvd on tin* history data fila.',/.
+ ' Insufficient tina history data.')
RETURN
RETURN
END
SUBROUTINE RHS(ID,T,TD,B,X.KS,HNOD,NEQN.MBAN>
C-SUBiRBS - FORMS RES OF [X>)(dT/dt| - IE')
C
C
C
C
C
C
C
IE*(t» -
lEMt)) -
IE(t»-[K](T(tl) ; FOR 'E'-OOF
(dT(t)/dt)*DIAE OF IK'l ; FOR • T'-DOF
IE-) IS WRITTEN OVER (TO)
IK] 4 IK*]
ADI AYSM-BANOED COMPACT STORED
IMPLICIT REAL'S (A-H.O-S)
REAL'S T.TD(NEQN),E(NEON>.K(NEaN.2*MBAN-l>,
•ltS(NtON.2'MBAN-l)
INTEGER ID(tMOO)
LOGICAL TD01
DO 20 I-l.NNN
SCALE BY DLIiSONAL FOR TEMP PRESCRIBED NODES
IP(TDOF(I,in.NNOO)> TEEN
TD(I) - TtKI)'KS(I.KBAN)
FORM (B]-[K](T| WHERE IK] IS IN COMPACT STORAGE
ELSE
TEMP - Ed)
Kl - KAX(l,MBAN-I+l>
X2 - HINP'MBAN-l.KBAFMTCON-I)
DO 10 XK-ia,K2
J - I » IDC - MBAN
TEMP - TEHP - K(X,KX)»T(J)
10 CONTINUE
TD(I) - TIIKP
ENDIF
20 CONTINUE
RETURN
END
FUNCTION TDOF(NEQ,ID.NNOD)
C-FUNlTDOF - DETERMINES IF EQUATION NUMBER NEQ IS A TEMPERATURE DOF
LOGICAL TOO?
INTEGER ID(SNOO)
TDOF - .FALSE.
DO 10 N-l.KXOD
IF I (ID(N).LT.O).AND.(ABSIID(N)).EQ.NEQ)I TBEN
TDOF - .TRUE.
RETURN
ENDIF
10 CONTINUE
RETURN
Append-19
-------
NBS: Indoor Air Quality Model
Phase II Report
Appendix - FORTRAN 77 Source Code
LOGICAL ERR
-RESET CALL VDATULMMPKEQ) , LMMPV) ,NFNOD, WECW.N.ERRl
SUBROUTINE RESET
C— SUB:RESET - COMMAND TO RESET CONTAM BY RE-IMTTIALI2INC POINTERS AND
C COUNTERS AND DELETES ARRAYS LEFT BY CONTAM IN BLANK COMMON
C — HELP LIST
C - C
RETURN
END
CALL INITCN
SUBROUTINE VDAT1 (KEO.,V,NFNOD,NFEQN,N,IRR)
C — SUBlVDATl - READS NODE VOLUME DATA
RETURN
END
C
c
c
c
c
INCLUDE IOCOM.INC
REAL«« VINFEQNI.VD
IF(VDAT.LT.O.ODO)
SUBROUTTNX NDCHK
C — SUBlNDCHK - CHECKS FOR OUT-OF-RANCE ELIMENT NODE NUMBERS
INCLUDE IOCOM.INC
LOGICAL ERR
DIMENSION ND(NDIM)
WRITE (NOT, 2000)
ERR - .TRUE.
RETURN
ENDIF
2000 FORMAT C •«•• ERROR. Hotel volun»tric «•• nay not b. n.gotiv..')
NEO - ABS(KEO REPORTS NODE VOLUME DATA
C
COMMON MTOT.NP, IA<1)
INCLUDE IOCOM.INC
INCLUDE CNTCOMS6.INC
LOGICAL ERR
EXTERNAL VDATO
WRITE (NTW, 2000)
WRITE (NOT, 2000)
2000 FORMAT!/, ' — Nodal Voluntric Ma»')
CALL DATGEN (VDATO, NFELM. ERR)
IF (ERR) RETURN
CALL REPRTN(IA(MPV).IA(NPXEQ).N7EaN.NrNOD)
RETURN
END
SUBROUTINE VDATO (N, ERR)
C-SHB.VDATO - CALLS VDAT1 PASSIKC ARRAYS
COMMON MTOT.NP,IA(1)
INCLUDE CNTCOmC.INC
ERR - .FALSE.
FIRSTL - .TRUE.
ADV C
C—1.0 GIT LINE OF DATA
C
100 IF(MODE.EO.'INTER') CALL PROMPT!' DATA>')
rar.T. FREE
IFIMOOE.EO.'BATCH') '•"' FREENR(NTW)
C
C—2.0 CHECK FOR •END"
C
IF (IOC) TEEN
IF(FIRSTL) THEN
WRITE(NTW,2200)
WRITE(NOT.2200)
2200 FORMAT)' •••• ERRORI Data anpactad; 'END* found.')
ERR - .TRUE.
RETURN
RETURN
ENDIF
ENDIF
C—1.0 GET INCREMENTING RULE; RETURN IF IJK(l) .EO.O
C
ATO IJK<1> * °
CALL FREEH1' '.IJKID.D
IFIIJXID.EQ.O) RETURN
IF(IJX(2).EQ.O) IJK(2)-IJK(1>
IF(IJKIl).EQ.O) IJKI31-1
DO 100 N-IJK(l), IJK(2), IJK(3)
IF(MAXNO.CT.O) CALL NDCHKfN,MAXNO, 1,ERR)
IF(ERR) RETURN
CALL DATAO (N, ERR)
Append-20
-------
NBS: Indoor Air Quality Model
Phase II Report
Appendix - FORTRAN 77 Source Code
IF(ERR) RETURN
100 CONTINUE
FIRST! - .FALSE.
CO TO 100
SUBROUTIN1
C— SUBiELEMIN -
C
C
C
C
C
C
C
C
C VARIABLE
C INPUT
C ELEHO
C NELDOF
C KEQ
C NSYNOD
C OUTPUT
C MSTBAN
C ERR
C LOCAL
C LMNEN.LHO!
C NOLD.NNEN
C INCR
C
I ELGENIELEMO.KEQ.NELDOF, NSYNOD, MSYBAN.ERR)
READS ELEMENT NUMBER, CONNECTIVITY, i GENERATION DATA
GENERATES HISSING ELEMENTS, UPDATES SYSTEM BANDNIDTH
PMLS -ELEMO- TO READ ELEMENT PROPERTY DATA
RETURNS KHEN DATA LINE IS BLANK, IS -|-. OR IS "END-
CHECKS ALL GENERATED NODE NUMBERS .LE. NSYNOD
•• CURRENTLY LIMITED TO FOUR-NODE ELEMENTS OR LESS ••
PROCEDURE NAME TO READ ELEMENT PROPERTY DATA
NUMBER OF ELEMENT DECREES OF FREEDOM
SYSTEH EQUATION NUMBERS (BY NODE NUMBER)
NUMBER OF SYSTEM NODES
SYSTEH BAND WIDTB
ERROR FLAG
U) ELEMENT LOCATION/CONNECTIVITY DATA
ELEMENT NUMBERS
GENERATION INCREMENT
RETURN
ENDIF
C GENERATE MISSING ELEMENTS
IFINNEN.CT.NOLD+1) TEEN
DO 24 N-NOtO-H.NNEM-1,1
DO 22 I-l.NELDOF
22 LMOLD(I) - LMOLD(I) f INCH
CALL NOCBX(LHOLD,NSTNOD,NELDOF,ERR)
IF(ERR) RETURN
CALL ELEHOLO,NSYNOD.NELOOF,ERR)
IF(ERR) RETUPN
CALL ELBAN(KIO,LHOLD.MSYBAN,NELDOP,NSYNOD)
CALL KLEMOtNOLD.LMOLO.ERR)
IP (ERR) RETURN
SO TO 20
2200 FORMAT!• •••• ERRORi El.rn.nt nunter ',15,' it out of ord«r.'1
END
SUBROUTINE ELBAN( KEO,LM,MSYBAN,NELOOF,NSTNOD>
C—SUBlELBAN - CONFUTES ELEMENT BANNIDTB t UPDATES SYSTEM BANDWIDTH
DIMENSION LM(NELDOF),XEQ(NSTNOD)
INCLUDE IOCOM.INC
LOGICAL ERR
INTEGER NEUX)F,LMNEN(4).UCLO(4),NOLD.NNEH,KEO(NSYNOO>
EXTERNAL ELEMO
C
C—1.0 art FIRST LINE OF ELEMENT DATA
C
INCH - 0
IF(MODE.EQ.'INTER') CALL PROMPTC DATA>'|
f»t,T- FREE
IFIMODE.EQ.'BATCH') CALL FREENR(NTN)
C CHECK FOR -END-
IFIEOC) RETURN
HOLD • 0
CALL FREEIC >,NOLD,1)
IF(NOLD.EO.O) RETURN
CALL FREEI('I',LMOLD(1),NELDOF)
CALL NDCHKILNOLD.NSYNOD,NELDOF,ERR)
IF(ERR) RETURN
<•"'- ELBANtKEQ.LMOLD.MSYBAN.NELDOF,NSYNOD)
CALL ELEHO (HOLD, LHOLD, ERR)
IF(ERR) RETURN
C
C—2.0 GET NEXT LINE OF ELEMENT DATA
C
20 IF(MODE.EQ.'INTER') CALL PROMPT!' DATA>')
CALL FREE
IF(MODE.EO.'BATCH') «*"-'- FREENR(NTN)
C CHECK FOR "END"
IFIEOC) RETURN
C GET NEH ELEMENT INFORMATION
NNEM - 0
CALL FREEIC ',NNEN,1|
IF(NNEN.EQ.O) RETURN
CALL FREEI('I',LMNEN(1),NELDOF)
CALL FREEH'N',INCR,1>
IF(INCR.EQ.O) INCR'l
C CHECK NUMERICAL ORDER
IF(NNEN.LE.NOLD) THEN
WRITEINTN.2200) NNEN
WRITE(NOT,2200) NNEN
ERR-.TRUE.
C
C INPUT
C LM
C' NELDOF
C HSYBAN
C KEO
C NSYNOD
C OUTPUT
C MSYBAN
e
ELEMENT LOCATION/CONNECTTVTTY
NUMBER OF ELEMENT DECREES OF
CURRENT SYSTEM BANDWIDTH
ARRAY
FREEDOM
SYSTEH EQUATION NUMBERS (BY NODE NUMBER)
NUMBER OF 8XSTKN NODES
UPDATED SYSTEM BAND WIDTH
MAX - ABS(XEfl(LM(l))>
KIN - ABS(KE()(LN(2)))
DO 10 1-1,NELDOF
NN - ABS(!a:0(LM(I»)
IF(NN.CT.MM) MAX-NN
IF(NN.LT.nH) MIN-NN
10 CONTINUE
MELBAN - MAX-MIN*!
IFIKELBAN.CT.HSYBAN) MSYBAN-MELBAN
RETURN
END
SUBROUTINE R!!PKTN(X.KEQ,NX,NK£Q>
C—SUBtREPRTN - RETORTS VECTOR
IMPLICIT REAL'S (A-H,O-Z>
INCLUDE IOCOH.INC
REAL'S X(NX),XX(5)
INTEGER KEO(lflCEQ)
CHARACTER*! irLG(S)
WRITE(NOT.2000)
IF(ECHO) WRITE INT",2000)
Append-21
-------
NBS: Indoor Air Quality Model
Phase II Report
Appendix - FORTRAN 77 Source Code
2000 FORMAT(/.
.13X,•••• - ind>p.nd»nt OOF.
.CX,4(2X, 'Mod* V.lvi«',JX»
DO 100 N-1.NKEQ.4
NN > Mni(N+J,NKEO.>
DO 10 I-N.NN,1
NEQ - KEQ(I)
NNEO - ABS(NEQ)
ITINEQ.LT.O) THIN
XXd-N+1) - X(KNEQ)
rLSd-N+1) - •••
ELSIir(NSQ.la.0) THIN
XX(I-N+1) - 0.000
ncd-N+i) - -o-
ELS!
xxd-N-fi) - x(NHEO>
FLSd-H*!) - ' '
ENDIF
10 CONTINUE
ir(ECBO) nun (NTH. 2010)
•U- . undefined DOFi.',//,
d,PLB(I-N+l(,XX(I-N+1),I-N.NN)
100 mUTE(NOT.2010) (I,FLC(I-Ntl),XXd-N+l).I-N.NN)
2010 FORMAT((«X,4(I6,lAl,eil.J)))
RETOIW
SUBROUTINE REPRTE(X,NX)
C—SUBlREPRTE - REPORTS VECTOR (X).IN ELEMENT ORDER SEQUENCE
C X(HX) - VECTOR OF VALUES ORDERED BY ELEMENT NUMBER
IMPLICIT RIAL'S(A-B.O-Z)
INCLUDE IOCOM.INC
DIMENSIOH X(NX)
WRITE (NOT. 2000)
IF(ECHO! WRITE (NTM, JOOO)
•RITE (NOT. 2010) IN, X(N). N-l.NX)
iriECBO) WRITE(NTH,2010) IN, X(N). N-l.NX)
2000 rORHAT(/.CX.4(2X.'Elra Vain'.IX))
2010 FORMAT((6X,4(I(,1X.C11.3)»
RETURN
END
SUBROUTINE FORMF(KXQ,r,VE,rORM.ERR)
C—SUBiFORHT - CALLS FORHFO TO FORM SYSTEM FLOW MATRIX
C ARRAY GOUT USED TO CHICK NODAL MASS FLOW CONTINUITY
COMMON MTOT.NP.IA(l)
INCLUDE CNTCOMSC.INC
IMPLICIT REAL'S(A-B.O-Z)
INCLUDE IOCOM.INC
INCLUDE CNTCOHSt.INC
REAL-9 F(NPEON.l), NE(NFELM), KLF<2,2), CONT(NFEQN),EFF
INTEGER KEQ(NFNOD), LM(2)
LOGICAL ERR
CBARACTER FORM-4
C
C—1.0 FOR EACH ELEMENT FORK ELEMENT PLOH MATRIX AND ADD TO [F]
C ACCUMULATE TOTAL MASS FLOW (CONTINUITY) AT EACH NODE
REMIND ND1
DO 10 N-l.NFELM
READINDl, ERR->00, END-900) LMdl, LMI2), EFF
N - «E(N)
Nl - ABS(KEOILMIl)))
N2 - ABS|KEO(LH(2)))
IF(W.ST.O.ODO) TBIN
ELFd.l) - •
ELF(1,21 - O.ODO
ELT(2,1) - -*•(!.ODO-BFFI
ELF (1,2) - 0.000
CONT(Nl) - CONT(Nl) * «
CONT(N2) - CONT(N2) - K
ELSIirdl.LT.O.ODO) TBEN
ELril.l) - O.ODO
BLF(1,2) - **(1.0DO-Err)
ELF(2,1) - O.ODO
ELT(1,2) - -N
CONT(N1) - CONT(Nl) * II
CONT(N2) - CONTIN2) - N
ELSE
CO TO 10
ENDir
IT (FORM.KQ.'BAND'I CALL ADDCA(KEQ,NreOD.ELr.r.2,KrEQN.MFBAN,LM)
IT (FORM.10.'FULL') CALL ADDA(lfflQ,NFNOO,SU,F,2,NTEW.LM)
10 CONTINUE
C
C—2.0 REPORT NIT TOTAL MASS FLOIt
C
VRXTS (NOT, 2200)
IF(ICBO) URITI(NTH,2200)
2200 rORMATI/,' — N*t Tot.l K>» Flo.')
CALL REPRTNICONT.KIQ.NrEON.NrNOO)
RETURN
>00 WRITE (NTM, 2100)
NRITE (NOT, 2*00)
2tOO FORMAT)
+ • •••• ERROR I Rwd or EOr urar on fln< «l«n«nt data fil«'l
ERR - .TRUI.
RETURN
RIAL'S riNraQN.l), NE(NrSLM)
INTEGER KEO(NmOD), MPCONT
LOGICAL ERR
CBARACTER FORH*4
CALL DELETE('CONT')
CALL OErTNRCCONT'.MPCONT.NrEQN.l)
CALL IEROR(IA(MPCOHT),NPEON.l)
IF (FORM. EQ. 'BAND') CALL ZERORIF, NPF.OX, 2'HFBAN-l)
ir(FORM.EQ.'FULL') CALL lERORIF.NTEON.NrEON)
CALL rORMrO(KEQ.r.NE,IA(MPCONT),rORM,ERR)
CALL DELETE C CONT')
RETURN
END
SUBROUTINE FORMFOIKEQ.P.NE. CONT. FORM. ERR)
C—SUBtFORMTO - FORMS SYSTEM PLOW MATRIX
C ARRAY CONT USED TO CHECK NODAL MASS FLON CONTINUITY
END
SUBROUTINE ADDCA(XEQ,NSYNOD,ELA,SYA,NELDOF,NSYDOF,MSYBAN, LM)
-SUBiAODCA - ADDS ILEMENT ARRAY TO COMPACT ASYMMBTRIC SYSTEM ARRAY
' REAL'S ELA(NELDOF.NELDOr), SYA(NSYDOF,1)
INTECER KEO(NSYNOa),LM(HELDOr)
C
C
C
C
C
C
C
C
C
C
C
C
VARIABLE
XEO(NSYNOD) t
NSYNOD I
ELA(NELDOF.NELDOF) i
SYA(NSYDOF, 2*HSYBAN-1) I
NILDOF i
NSYDOF i
MSYBAN i
LM(NELDOF) t
SYSTEM NODAL EQUATION NUMBERS
NUMBER OF SYSTEM NODES
ELEMENT ARRAY
COMPACTED ASYM. SYSTEM ARRAY
NUMBER Or ELEMENT DECREES OF FREEDOM
NUMBER OF SYSTEM DECREES OF FREEDOM
EALF BANDNIDTB OF SYSTEM ARRAY
ELEMENT LOCATION/CONNECTIVITY
00 20 I-l.NELDOF
Append-22
-------
NBS: Indoor Air Quality Model Appendix - FORTRAN 77 Source Code
Phase II Report
II ' ABS(KEO.(LM(I))) INCLUDE IOCOH.INC
DO 10 J-l.NELDOF CHARACTER STRING* (M
JJ - MSYBAN - II * ABS(KEO(LM(J))> NRITEINTII. ' (Jk,\) •) STRING
SYA(II.JJ) - SYA(II,JJ) * ELA(I.J) RETURN
10 CONTINUE END
20 CONTINUE
RETURN C
END SUBROUTINE PKOMB(N)
C—SUBtPXOMS - -HOLLERITH PROMPT*
SUBROUTINE ADDA(KEO,NSYNOD.ELA,SYA,NELDOF,NSYDOP,LM) COMMON MTOT.UP,IA(1)
INCLUDE IOCOH.INC
CHARACTIR*! ItCMND.M
COMMON /CMND/ NCMNDI8) ,M(4,7)
C PROMPT FOR A1UIAY NAMES
IF(MODE.EO.'BATCH') CO TO 900
DO 200 I-l.N
100 IF(M(1,N).NE.' ') CO TO 200
WRITE(NTW, 2000) N
CALL FREE
CALL FREECC ',Mil.N),8,1)
CO TO 100
200 CONTINUE
C
900 RETURN
2000 FORMATC •• Entar array nana *',1I1,'*
DO 20 I-1,NELDOF END
II - ABS(XEO(LM(I»)
DO 10 J-1,NZLDOF C
C— SUBiADDCA - ADDS ELEMENT ARRAY TO FULL ASYMMETRIC SYSTEM ARRAY
C
REALM ELAINELDOP.NELDOF)
, SYA(NSYDOF.l)
INTEGER KEQ(NSYNOO),LM(NELDOF)
C
C
C VARIABLE
C XEOINSYNOD) :
C NSYNOD >
C ELA(NELDOF,NXLDO?) I
C SYA(NSYDOF,2*MSYBAN-1) :
C NELBOF :
C NSYDOF >
C MSYBAN 1
C LM(NELDOF) t
C
DESCRIPTION -----..— _——----. i ._
SYSTEM NODAL EQUATION NUMBERS
NUMBER OP SYSTEM NODES
ELEMENT ARRAY
COMPACTED ASYM. SYSTEM ARRAY
NUMBER' OF ELEMENT DECREES OF FREEDOM
NUMBER OF SYSTEM DECREES OF FREEDOM
HALF BANDWIDTH OF SYSTEM ARRAY
ELEMENT LOCATION/CONNECTIVITY
JJ - ABS(KEO(LM(J») SUBROUTINE PROM (NR.NC)
SYA(II.JJ) - SYAIII.JJ) * ELA(I,J) C—SUB: PROM - 'INTEGER PROMPT*
10 CONTINUE
20 CONTINUE INCLUDE IOCON.IHC
RETURN
END C ASK FOR NUMBER OF ROWS AND COLUMNS
IF(MODE.SO.'BATCH') CO TO 900
C C 100 IF(NR.CT.O) SO TO 200
C , C «•"-'- PROMPT!' •• Entar nunfaar of rotui •)
C COMMAND PROCESSOR UTILITIES C «•"' FREE
C C CALL FREEH' '.NR.ll
C C CO TO 100
C
C NOPEN 200 IF(NC.CT.O) CO TO 900
SUBROUTINE NOPEN (LUN.niAME.FRM) CALL PROMPT C •• Ent.r nunbu of column
C—SUB. NOPEN - OPENS A FILE AS A HEN FILE WHETHER IT EXISTS OR NOT CALL FREE
C LUN - LOGICAL UNIT NUMBER CALL FRESIC ',NC,1)
C FNAME - FILENAME CO TO 200
C FRM - FORM; 'UNFORMATTED' OR 'FORMATTED' C
900 RETURN
INTEGER LUN END
CHARACTER rNAME'O, FRM* (•)
LOGICAL FOUND C
SUBROUTINE HBORT
INQUIRE (FILE-FNAME.EXIST-FOUND) C—SUB I ABORT - ABORTS COMMAND AND RETURNS TO INTERACTIVE MODE
IF (FOUND) THEN C
OPEN(LUN,FILE-FNAME.STATUS-'OLD',FORM-FRM) INCLUDE ICCON.INC
IFIFRM.EO.'FORMATTED') THEN
WRITE(LUN,2000) LUN WRITE(NTH,2COO)
i FORMAT(K) WRITE (NOT, 2000)
ELSEiriFRH.EQ. 'UNFORMATTED') TEEN 2000 FORMAT!' •••• COMMAND ABORTED')
WRITE(LUN) LUN IF(MODE.BO.'BATCH') CALL RETRN
EHDIF
CLOSE(LUN,STATUS-'DELETE') RETURN
OPEN(LUN,FILE-FNAME,STATUS-'NEN',FORM-FRM) END
ELSE
OPEN (LUN, FILE-FNAME, STATUS- • NEW', FORM-FRH)
ENDIF
RETURN
END
SUBROUTINE PROMPT(STRINC)
C—SUB t PROMPT - INLINE PROMPT
C
C CALSAPX. LIBRARY
C
C AN EXTENSION OF "CAL-SAP" LIBRARY OF SUBROUTINES
C DEVELOPED BY ED WILSON, U.C. BERKELEY
C
C 1.0 FREE-FIELD TOPUT SUBROUTINES
C
r
C
C
C
C
C
FREE
Append-23
-------
NBS: Indoor Air Quality Model
Phase II Report
Appendix - FORTRAN 77 Source Code
SUBROUTINE FREE
C—SUBiFREE - READ LIMI 07 FREE FIELD DATA
C COWffiVTS LINES ECHOED TO SCREEN
WRITE(LUN,2000) (LINE (I). 1-1, JJ)
2000 FORMAT (IX, SOU)
INCLUDE IOCOH.INC
INCLUDE FRECOM.INC
RETURN
END
C-0.0-INITIALIZE VARIABLES
C
EOD - .FALSE.
EOC - .FALSE.
DO 5 1-1,160
5 LINE(I)-' '
C
C-1.0 GET LINE OF DATA
C
10 I - 1
II- eo
READ(NCND, 1000,ERR-100) (LINE(K) ,K-I.III
SUBROUTINE FRXEFNISEP.NC, FOUND)
C—SUBlFREEFN - FINDS NEXT NC-CBARACTER SEPARATOR IN INPUT FILE
C SEP(NC)«1 - CBARACTER STRING
INCLUDE IOCOH.INC
INCLUDE FRECOH.INC
CHARACTER'1 SEP(NC)
LOGICAL FOUND
FOUND - .FALSE.
C CHECK FOR ADDITIONAL LINE
JJ - LEHTRH(LLINS)
DO 12 R-I.JJ
IKLINKKl.EO.'V) TEEN
I - K
II- K*79
READ(NCMD, 1000,ERR-100) (LXNE(XX) ,KX-I,II)
1000 FORMAT I «OA1)
60 TO 14
EKDIP
12 CONTINUE
SO CALL FREE
IF(NC.LE.II) THEN
DO iO N-l.NC
(0 IF(SEF(N).NE.LINE(N)) GO TO SO
FOUND • .TRUE.
RETURN
ELSE
eo TO so
ERDIF
RETURN
END
C CHECK FOR COMMENT
14 IFILHiEID.EQ.'•') THEN
IF (HODS.EQ.1 BATCH1) CALL FREENR(NTH)
CALL FRZEWR(NOT)
GO TO 10
END IF
C
C-2.0 DETERMINE LENGTH-OF-INFORMATION
C
JJ - LENTRH(LLHR)
C
C-3.0 DETERMINE LENCTH-OF-DATA AND CONVERT DATA TO UPPER CASE
C
ISP - ICHARC •)
LA - ICHAR('l')
DO 10 I-l.JJ
IF(LINE(I).EO.•<') GO TO 32
NN - ICBARILINEID)
IF(NN.GE.IA) LINE(I) - CBAR(NN-ISP)
10 CONTINUE
32 II - I - 1
C
C-«.0 CHECK FOR END-OF-DATACROUP t END-OF-COMMAND
C
IFILINE(l).E0.'<'> EOD - .TRUE.
IF(LINE(1)//LINE(2)//LINE(3).EO.'END<) EOC- .TRUE.
RETURN
C ERROR IN READ
100 WRITE(NOT,2000)
WRITE(NTN,2000)
2000 FORMATC •••• ERRORi Error in r«.ding input line.')
CALL ABORT
SUBROUTINE FREENR(LUN)
C—SUB 1 FREENR - WRITE COMMAND/DATA LINE TO FILE LUN
C LUN - LOGICAL UNIT NUMBER TO WRITE TO
INCLUDE lOCOM.INC
INCLUDE FRECOM.INC
SUBROUTINE FRIERIIC,DATA,HUM)
C—SUB i FREER - FIND AND INTERPRET REAL DATA
C IC'l - DATA IDENTIFIER CBARACTER
C DATA - REAL DATA RETURNED
C HUM - NUMBER OF DATA VALUES TO EXTRACT
IMPLICIT RIAL'S (A-B.O-l)
DIMENSION DATA(IO)
CHARACTER IC'l
INCLUDE FRECOM.INC
C FIND REAL STRING
10 1-0
IFIIC.EO.' ') GO TO 230
DO 100 1-1,11
IF((LINE .AND. (LINEd-H) .EO.'-'» CO TO 2SO
100 CONTINUE
RETURN
C EXTRACT REAL DATA
290 DO 2CO J-l.NUH
2tO DATA(J)-0.0
DO 300 J-l.NUM
JJ-0
270 IF(I.CT.II) GO TO 300
CALL FREERIII.XX.NN)
IF(JJ.NE.O) CO TO 27S
DATA(J) - XX
CO TO 290
C ARITHMETRIC STATEMENT
279 IF(JJ.EO.l) DATA(J)-DATA(J)*XX
IFIJJ.EQ.2) DATA (J)-DATA (JI/XX
IFIJJ.E0.3) DATA(J)-DATA(J)*XX
IFIJJ.EQ.4) DATA (J)-DATA (J)-»
IFIJJ.NI.9) CO TO 290
C EXPONENTIAL DATA
JJ - DABS(XX)
IF(JJ.EQ.O) CO TO 290
DO 2<0 K-l.JJ
IF(XX.LT.O.O) DATAIJ) - OATA(J)/10.
IF(XX.CT.O.O) DATA(J) - DATA(J)'10.
260 CONTINUE
C SET TYPE OF STATEMENT
290 JJ-0
IF(LINE(I).EQ.'"I JJ-1
Append-24
-------
NBS: Indoor Air Quality Model
Phase II Report
Appendix - FORTRAN 77 Source Code
IF (LINE (I) .IQ. V) JJ-2
IF(LINEd) .EO.'*') JJ-]
IF(LINEd) .EQ.'-'> JJ-4
IP(LINE(I).EO. •£•> JJ-5
IFtLINE JJ-0
IF(JJ.NE.O) GO TO 210
IFINN.CT.9) RETURN
300 CONTINUE
RETURN
END
SUBROUTINE FREER1(I.XX.NN)
C-SUBtFREERl - INTERPRETS X SINGLE REAL VALUE
IMPLICIT REAL'S (A-H,C~!>
INCLUDE FRECOM.INC
C CONTORT STRING TO REAL FLOATING POINT NUMBER —
IF (LINE (1*1) .EO. •-'( 1-1*1
T-0
IS-1
XX-0.0
IFillNEtl+D.EO.'-') TEEN
IS—1
I-I+l
ELSEIF(LINE(I+l).EO.t + <) THEN
IS-1
I-I+l
ELSE
CONTINUE
ENDI7
261 IF(LINE.NE.' 'I CO TO 210
I-I+l
17(1.CT.II) CO TO 300
GO TO 2(1
210 I-I+l
IFd.ST.II) CO TO 300
IFKLIKSdl.EO.' •).AND.(LINE(I+1).EQ.> ')) CO TO 210
NN - ICSAKI LINE (I) ) - ICHARCO1)
XN-ISIGN(NN,IS)
IFILINSd) .ME.'.') CO TO 215
t-1.0
GO TO 210
215 IF(LINE(II.ED.' •) CO TO 300
IF(LINE(I).EQ.', ') CO TO 300
IFKNN.LT.Oj.OR. (NN.CT.9» CO TO 300
IFIY.EO.O) GO TO 2SO
Y-T/10.
XN-Xtfl
XX-XX+XM
CO TO 210
280 XX-10.>XX+XN
CO TO 210
300 RETURN
END
SUBROUTINE FREEIIIC.IDATA.NUM)
C—SUBiFREEI - FIND AND INTERPRET INTEGER DATA
C IC>1 - DATA IDENTIFIER CHARACTER
C IDATA - INATECER DATA RETURNED
C HUM - NUMBER OF DATA VALUES TO EXTRACT
CHARACTER*! IC.LNE
DIMENSION IDATAI12)
INCLUDE FRECOM.INC
DO 250 J-l.NUM
ISICN - 1
C SKIP BLANKS BETWEEN INTEGERS
215 IFILINEd+D.NE.' ' I CO TO 220
I-I+l
IF(I.CT.II) GO TO tOO
GO TO 215
220 I-I+l
IF(I.CT.II) GO TO 230
C CHECK FOR SIGN
LNE - LINE (I)
IFtlNE.NE.1-1) CO TO 225
ISICN - -1
CO TO 220
C EXTRACT INTEGER
225 IPONE.EQ.' ') CO TO 230
IFUNE.EO.',') CO TO 230
IF(LNE.EO.'I') GO TO 230
NN - ICBAR(LKE) - ICBAR('O')
IF((NN.LT.O).OR.(NN.CT.9)> CO TO 900
IOATA (J)-10* I DATA (J)+NN
CO TO 220
C SET SICN
230 IDATA(J) - IDATA(J)•ISICN
ISO CONTINUE
900 RETURN
END
C FIND INTEGER STRING
90 1-0
IFIIC.EQ.' ') CO TO 200
DO 100 1-1,11
IF(ILINE(I).EO.IC).AND.(LINE(I+1).EO.'-')) GOTO 200
100 CONTINUE
RETURN
C IERO INTEGER STRING
200 DO 210 J-l.NUM
210 IDATA(J)-0
IFILINE(Itl).EO.'-'I I-I+l
SUBROUTINE HIEEC (1C, IDATA. NC, HUM)
C—SUBtFREEC - FIND AND INTERPRET CHARACTER DATA
C IC<1 - DATA IDENTIFIER CHARACTER
C IDATA - CHARACTER DATA RETURNED
C DC - NUMBER OF CHARACTERS PER DATA VALUE
C HUM - NUMBER OF DATA VALUES TO EXTRACT
CHARACTER-1 1C,IDATA
DIMENSION IDJlTAINC.NUH)
INCLUDE rRECOM.INC
C FIND DATA IDIWTIFIER
90 1-0
IFdC.EQ. • •! CO TO 200
DO 100 1-2, II
IF((LINE(I-1| .EO.IO.AND.ILINEID.EO.'-')) GOTO 200
100 CONTINUE
RETURN
C EXTRACT CHARACTER DATA
200 DO 210 J-l.NIIM
DO 210 N-1,W:
210 IDATA(N.J)-' >
C
DO 300 J-l.NUM
2(0 I - I + 1
IF(I.CT.II) IX) TO 400
IFUINEdl.EO.',') CO TO 2(0
IF(LINEd) .E(J. ' •) CO TO 2(0
DO 290 N-l.NC
IFILINEd) .E(J.'i') GO TO 300
IF(LINE(I).EU.' •) CO TO 300
IF (LINE (I). Ell.',' I CO TO 300
IDATA(N,J) - LINE (I)
IFIN.EO.KC) GO TO 290
I-I + l
290 CONTINUE
300 CONTINUE
400 RETURN
END
FUNCTION LENTRH(STRING)
C—FUNiLENTRN - DETERMINES LENGTH OF TRIMMED STRING - A STRING KITH
C mULINC BLANKS REMOVED
C
C LENTOT I THE TOTAL LENGTH OF THE STRING
C LENTRM I THE LENGTH OF THE TRIMMED STRING
C
CHARACTER STRING* <*>
INTEGER LENTOT, LENTRM
Append-25
-------
NBS: Indoor Air Quality Model
Phase II Report
Appendix • FORTRAN 77 Source Code
LENTOT - LEN (STRING)
DO 10 I-LENTOT, l.-l
IFISTRING(IlI).NE.'
10 CONTINUE
20 LENTRM - I
RETURN
END
C 2.0 DYNAMIC ARRAY MANAGEMENT
C
SUBROUTINE DEFINR(NAMB,NA.NR.NC)
C—SUBiOEFINR - DEFINE DIRECTORY AND RESERVE STORAGE
C TOR REAL ARRAY IN DATABASE
C NAME - NAME OF ARMY
C NA - BLANK COMMON POINTER TO ARRAY (RETURNED)
C UK - NUMBER OF ROTS
C NC - NUMBER OF COLUMNS
C
COMMON MTOT.NF, IA(1)
CHARACTER*! NAME<4)
NP - 2
CALL DZFIN(NAME.NA.NR.NC)
RETURN
SUBROUTINE DEFINI(NAME,NA,NR,NC)
C—SUBiDEFINI - DEFINE DIRECTORY AND RESERVE STORAGE
C FOR INTEGER ARRAY IN DATABASE
C NAME - NAME OF ARRAY
C NA - BLANK COMMON POINTER TO ARRAY (RETURNED)
C MR - NUMBER OF RONS
C NC - NUMBER OF COLUMNS
C
COMMON MTOT.NF, XAU)
CHARACTER*! HAM*(4)
NP - 1
CALL DEFIN(NAME,NA.NR,NC)
RETURN
END
SUBROUTINE DEFINC(NAKE,NA,NR,NC)
C—SUBiDEFINC - DEFINE DIRECTORY AND RESERVE STORAGE
C FOR CHARACTER*! (HOLLERITH) ARRAY IN DATABASE
C NAME - NAME OF ARRAY
C NA - BLANK COMMON POINTER TO ARRAY (RETURNED)
C NR - NUMBER OF RONS
C NC - NUMBER OF COLUMNS
C
CHARACTER'1 NAME(4)
COMMON MTOT.NP.IAd)
NP - 1
CALL DEFININAME.NA.NR.NC)
RETURN
END
SUBROUTINE DEFIN(NAKE,NA,NR,NC)
-DEFINE AND RESERVE STORAGE FOR ARRAY
C
C
C
C
C
C
C
C
c-
c
C
C
C
C
C
C
C
C
c-
c
C
C
C
C
C
C
C - EVALUATE STORAGE REQUIREMENTS -
NSIIE -
N8IIX - NSIIE*2 * 2
NA - NEXT
NEXT - NEXT -I- NSIZE
C - SET UP NEN DIRECTORY -
HUMA - NUMA * 1
IDIR - IDIR - 10
I - IDIR
C - CHECK STORAGE LIMITS -
IF(I.GE.NXXT) GO TO 100
I - NEXT - I » MTOT - 1
KRITE (NTH, 2000) I. MTOT
•RITE (NOT. 2000) I. MTOT
PAUSE
STOP
100 CALL ICOH(NAMX.IA(I))
»(!«•«> - NR
IAII+S) - NC
IAII+4) - NP
LMI+7) - NA
IA(I+» - NSIZE
IA(It») - 0
RETURN
NEXT - NEXT AVAILABLE STORAGE LOCATION
IDIR - START OF DIRECTORY IN BLANK COMMON
IP - NUMBER OF LOCICALS CONTAINED IN DATA TYPE
LENR ' NUMBER OF LOCICALS IN PHYSICAL RECORD
NP ' TYPE OF DATA
• 1 INTEGER DATA
- 2 REAL DATA
- 3 LOGICAL DATA
-DIRECTORY DEFINITION FOR CORE OR SEQUENTIAL FILES
IDIR(l.N) - NAME OF ARRAY - INANE (4 CBAR.)
IDIR(S.N) - NUMBER OF RONS - NR
IDIR«,N) - NUMBER OF COLUMNS - NC
IDIRP.N) - TYPE OF DATA - NP
IDIRO.N) - INCORE ADDRESS - NA
- -1 IF SEQUENTIAL FILE ON DISK
- -2 IF DIRECT ACCESS ON DISK
IDIRO.N) - SHE OF ARRAY
IDIR(IO.N) - 0 IF IN CORE STORAGE
-DIRECTORY DEFINITION FOR DIRECT ACCESS FILES
IDIR(S.N) - NUMBER OF INTEGERS
IDIRK.N) - NUMBER OF REAL WORDS
IDIRO.N) - NUNBER OF LOCICALS
IDIR(a.N) - NUMBER OF LOGICAL RECORDS
IDIR(I.N) • LOGICAL RECORD NUMBER
IDIR(IO.N) - LUN IF ON LOGICAL UNIT LUN
fOO
2000 FORMAT)
ERROR i Insufficient blank COMMON itor«5..'./,
Storage nquind NTOT -',n./.
Stongi «»«il«tl« MTOT -'.IT)
SUBROUTINE DEFDIR(NAMX.NR,NC,ISTR)
C—SUBiDEFDIR - DEFINE DIRECTORY FOR OUT-OF-CORE FILE
C NANS - NAME OF ARRAY
C NR NUMBER OF RONS
C NC - NUMBER OF COLUMNS
C ISTR - OUT OF CORE FLAG (—1)
C :
COMMON MTOT,NF,IA(1>
INCLUDE ARYCOM.IMC
INCLUDE IOCOM.INC
COMMON MTOT.NP.IAd)
INCLUDE ARYCOM.INC
INCLUDE IOCOM.INC
CHARACTER*! NAME (4)
•OEFIN VARIABLES
NAME - NAME OF ARRAY - 4 LOCICALS MAXIMUM
NA - LOCATION OF ARRAY IF IN BLANK COMMON
NR - NUMBER OF RONS
NC - NUMBER OF COLUMIS
MTOT - END OF DIRECTORY
NUMA - NUMBER OF ARRAYS IN DATA BASE
CHARACTER*! NAME(4)
C EVALUATE STORAGE REQUIREMENTS
IF(NP.EO.O) NP - 2
C SET UP NEM DIRECTORY
NUMA - NUMA » 1
IDIR - IDIR - 10
I - IDIR
C CHECK STORAGE LIMITS
IFII.GE.NEXT) GO TO 100
I - NEXT - I + MTOT - 1
NRITE(NTN.2000) I,MTOT
Append-26
-------
NBS: Indoor Air Quality Model
Phase II Report
Appendix - FORTRAN 77 Source Code
MUR (NOT, 2000) I.KIOT
PAUSE
STOP
100 <*"-'- ICON (NAME. IA (I))
IA(I+4) - NK
IA(I+5> - KC
IA(I»«> - N?
1AII+7) - ISTR
IA(I*8> - 0
IAII+9) - 0
900 RETURN
2000 FORMAT!
.. .... ERRORi inefficient blank COMMON •toragi.1,/,
•' Storag* raquind MTOT -',I7,/,
•' Storag* available MTOT -',17)
END
SUBROUTINE LOCATE(NAME, NA,NX, NC>
C—SUBlLOCATE - LOCATE ARRAY -NAME- AND RETURN
C NA - POINTER TO LOCATION IN BLANK COMMON
C NX - NUMBED OF ROWS
C NC - NUMBER OF COLUMNS
C
COMMON MTOT.N7.IAm
CHARACTER-1 NAME
DIMENSION NAME(4),INANE(1)
C LOCATE AND RETURN PROPERTIES ON ARRAY
NA - 0
CALL ICON (NAME, INAME)
I - IFIHD (INAME, 0)
IF(I.EQ.O) GO TO >00
C RETURN ARRAY PROPERTIES
NA » lA(I-fl)
NR - IAII-H)
NC - IAII+S)
NP - IA(I+«)
900 RETURN
DO «SO J-1,10
IAIII) - IA(It-10)
«50 II - II - 1
IF(IA(I»7).LE.O> CO TO 840
IF(IA(I+9>.EQ.O) IA(I+7> - IA(I*7) - NSIZE
• (0 I - I - 10
C
900 RETURN
1000 FORMAT!' — Nan> '.4A1, • it biing ui.d for an'.
. t OUT-OF-COR! 111*.',/)
END
SUBROUTINE ICON (NAME. INAME)
CHARACTER-1 NUMB HI
DIMENSION INAKEI4)
C^^—CONVERT LOGICALS TO INTECER DATA
-LOCATE DO 100 I - 1, «
100 INAME(I) - ICBARI NAHS(I) )
C
RETURN
END
FUNCTION IFINO(INAME, LUN)
C—FUNtiriND - FIND
COMMON MTOT.NP, IA(1)
INCLUDE ARYCOM.INC
DIMENSION INAMK4)
-FIND ARRAY LOCATION
I - IDIR
DO 100 N-l.NUKA
SVBROUTINS DELETE(NAME)
C—SUBtDELZTE - DELETE ARRAY -NAME- FROM DATABASE
IT(LUN.NE.IA(I+9H CO TO 100
IF (INAME(l).NI.IA(I » GO TO 100
IF (INAME(2).NI.IA(I*1» CO TO 100
IF (INANE|l>.NE.IA(I+2)> CO TO 100
IF (INAME(4> .IQ.IA(IO)) GO TO 200
100 I - I » 10
I - 0
200 IFIND - I
RETURN
END
COMMON MTOT.NP.IA(l)
INCLUDE ARYCOM.INC
INCLUDE IOCOM.INC
C 3.0 MATRIX OPERATION UTILITIES
C
CHARACTER-1 NAME
DIMENSION NAME(4),INAME(4)
C DELETE ARRAY FROM STORAGE
100 "'-'- ICON (NAME, INAME)
I - IFIND(INAME,0)
IF(I.EQ.O) CO TO 900
C CHECK ON STORAGE LOCATION
200 NSIIE - IA(I+8)
C SET SItB OF ARRAY
NEXT - NEXT - NSIZE
NUMA - NUMA - 1
NA - IA(I+7>
C CHECK IF OUT OF CORE OR DIRECT ACCESS
IF(NA.CT.O) CO TO SOO
NRITEflRlf.lOOO) NAME
XRITE (NOT,1000) NAME
CO TO SOO
SOO IT (NA.EG.NEXT) CO TO 100
C COMPACT STORAGE
II - NA » NSIZE
NNXT - NEXT - 1
DO 700 J'NA.NNXT
IAIJ) - lA(II)
700 II - II « 1
C—COMPACT AND UPDATE DIRECTORY
800 NA > I - IDIR
IDIR - IDIR * 10
IF(NA.EQ.O) GO TO 900
NA • NA/10
DO 8«0 K-l.NA
II - I » 9
SUBROUTINE ZEROI(IA,NR,NC)
C—SUBitlRORI - SET ARRAY IA(NR,NC) TO 0
DIMENSION IA(NR,NC)
DO 10 I-l.NR
DO 10 J-l.NC
IAII.J) - 0
10 CONTINUE
RETURN
SUBROUTINE ZERORIA.NR.NC)
C—SUBlZEROR - SET ARRAY A(NR,NC> TO O.O
REAL** A(NR.NC)
DO 10 I-l.NR
DO 10 J-l.NC
A(I.J) - O.ODO
10 CONTINUE
RETURN
END
SUBROUTINE F»CTCA(A,NEO,MBAND,ERR>
C—SUBlFACTCA - FACTORS COMPACT ASYMMETRIC MATRIX
C FACTORS (A] - IL](U)
C [L] [U] IS WRITTEN OVER [A]
C (A] MAYBE SYM OR ASYM, POSITIVE DEFINITE
C (Al BAS SEMI-BANDMIDTB MBAND 1 IS STORED COMPACTLY
C FROMl BUEBHER I THORNTON "TBE FINITE ELEMENT METHOD FOR ENCRS.
Append-27
-------
NBS: Indoor Air Quality Model
Phase II Report
Appendix - FORTRAN 77 Source Code
IMPLICIT REAL'S (A-B.O-Z)
INCLUDE IOCOH.INC
DIMENSION A|NBQ,2-MBAND-1)
LOGICAL ERR
NCOLS > 2-MBAND-l
mm - KBAND * i
DO 50 N-l.NEO
IF(X(N,HBAND).EQ.O.ODO) SO TO CO
IFIAIN,KBAND).BO.l.ODO) CO TO 20
C - l.ODO/AIN.MBAND)
DO 10 K-KMIN, NCOLS
IFIA(N,K>.EQ.0.000) CO TO 10
A(N.K) - C*A(N,K>
10 CONTINUI
20 CONTINUE
DO 40 L-2.MBAND
LL - NBAND + 1
DO 50 M-l.NEO
N - NEQ + 1 - M
DO 40 L-LL, NCOLS
ir(A(N,L).EO.O.ODO) CO TO 40
K - N * L - NBAND
B(N) - B(N) - A IN, LI «B IK)
40 CONTINUE
50 CONTINUE
RZTURN
(0 ERR - .TRUE.
BUTE (NTH, 2000) N
RETURN
2000 FORMAT(• "" ERROR) SUBiSOLVCA - Equation! may ba (Ingular.',/,
+ ' Diagonal of aquation numbar ',15,' ia xaro.'}
END
SUBROUTINE EICEN2(A,T,N.TKX.IP)
C—SUB I EICEN2 - Unaymnatrie Eigan Analyaia Boutin.
JJ - HBAND - L + 1
I - N « L - 1
IF(I.CT.NEO) CO TO 40
IF(A|I,JJ).EO.O.ODO) CO TO 40
KI ' KBAND » 2 - L
XT - NCOLS + 1-1
J ' MBAKO
DO 30 K-KI.KT
J - J + 1
ir71
Solv.i aiganproblan for raal matrix A(N,N). lyn. or uniyn., by
a laqnanca of Jaeobi-liKa tranafomationi [T]-1(A][T] «h.r. (T)-
[T1][T2][T3] Each [Ti] ii of tha fern IRI] [Si] wh.r«;
Rt Rk.k - Rm,n - CM(«) ; Rn,k - -R)t,n - •in(x)
Ri,i - 1 ; Ri,j - 0 l (i,J - x,nl
Si Sk.k • Sm.iB - eoah(y) i Sm,k - Sk,n - ->inh|y)
Si. i - 1 1 Si, j - 0 1 (i.J - k.m)
in which x, y ara datarminad by tha alanant* of [Ai] .
In tha limiting matrix raal aiganvaluaa occupy tha diagonal whila
raal and imaginary part! of eomplax aigan valua* occupy tha
diagonal and off-diagonal eornart of 2x2 bloeka cantarad on diog.
Array T(N,N) mult ba providad to raeaiva aiganvaetor! .
TMX*0 i aiganvaetora not ganaratad and A(N,N) may ba
paaiad aa T(N,N)
TMX<0 : ganarata laft, (TJ-1, tranaformationi
TMX>0 t ganarata right, [T), tranafarmationt
Eiganvactora of raal aiganvaluaa occurr aa row! Icoli) of [T]— 1
I(T]). Eiganvaetori for a eonplax aiganvalua pair aj,i 1 iaj,j+l
nay ba formad by tj 1 ltj+1 Hhara tj. tj+1 ara th. corraiponding
mm (cola) of [Tl-1 I(T)>
Itarationa ara limitad to 50 maxinum. On axit from tha proeadura
TKX raeorda tha numbar of itarationa parformad. Failur. to
eonvarga ia indicatad by TMX-50 or. if all tranaformationi in
Tha maehina dapandant variabla EP ia aat to 1E-06 and ahould ba
raaat for maehina praoiaion availabla.
—INPUT
A(N.N) Array to ba analyiad.
N . Syatam «iaa
TMS Control paramatar
C^— OUTPUT
C
C
C—
C
T|N,N) Array to raeaiva aiganvactora .
TMX Itaratlon count/it. ration flag
—LOCAL
EP Praeiaion
DO 20 L-2,KBAND
JJ - HBAND - L + 1
I; - N » L -1
IP(I.CT.NZQ) CO TO 20
IF(Ad.JJ) .EO.O.ODO) CO TO 20
BID - BID - A|I.JJ)*B|N)
20 CONTINUE
30 CONTINUE
C—2.0 BACIUUBSTITVnON
IMPLICIT REAL'S (A-B.O-t)
REAL** A(N.N),TIN,N),EP
INTEGER N.TKX
LOCICAL MARK, LEFT, RIGHT
C
C—0.0 INITIALIZE CONTROL VARIABLES
C
IF(EP.LE.O.ODO) EP - l.OD-e
EPS - SORTIE?)
LEFT - .FALSE.
RIGHT - .FALSE.
Append-28
-------
NBS: Indoor Air Quality Model
Phase II Report
Appendix - FORTRAN 77 Source Code
IFITMX.LT.O) THEN
LEFT - .TRUE.
ELSEIF(TMX.CT.O) TBEN
RICH! - .TRUE.
CNOIF
MARK - .FALSE.
C
C— 1.0 INITIALIZE [T] AS IDENTITY MATRIX
C
IF(THX.NE.O) THEN
DO 10 I-l.N
T(I,I) - l.ODO
DO 10 J-m,N
T(I.J) - 0.000
T(J,I) - O.ODO
10 CONTINUE
ENDIF
C
C— 2.0 MAIN LOOP
C
DO 26 IT-1,30
C
C— 2.1 IF MARX IS SET
C TRANSFORMATIONS OF PREVIOUS ITERATION HERE OMITTED
C PROCEIOURE HILL NOT CONVERGE
C
IF (MARK) TBEN
TMX - 1-IT
RETURN
ENDIF
C
C— 2.2 COMPUTE CONVERGENCE CRITERIA
C
DO 20 I-l.N-1
All - Ad, II
DO 20 J*I+1»N
AIJ - A(I,J)
AJI - A(J.I)
IF ((ABS (AIJ+AJI) .CT.EPS) .OR.
* ((ABS(AIJ-AJI). CT.EPS). AND. (ABSIAII-AIJ, J) ) . CT.EPS))) THIN
GOTO 21
ENDIF
20 CONTINUE
TMX - IT -1
RETURN
C
C~~2 . 3 BEGIN NEXT TRANSFORMATION
C
21 MARX - .TRUE.
DO 29 K-l.N-1
DO 29 M-K+l.N
B - O.ODO
C - O.ODO
BJ - O.ODO
YB - O.ODO
DO 22 I-l.N
AIK - AIX.K)
AIM - A(I,M)
TE - AIX-AIK
TEE m AIM*AIM
YB - YB + TE - TEE
IF((I.NE.X).AND.(I.NE.M)I TEEN
AKI > A(K, II
AMI - A(M, II
B * B * AKI-AHI - AIK'AIM
TEP - TB + AMI-AMI
TEM - TEE * AKI-AKI
G - G » TEP t TEM
BJ - BJ - TEP + TEM
ENDIF
22 CONTINUE
B - 8 * B
D - AIK, HI - AIM. Ml
AKM - A(X.M)
AKK - A(M. HI
C - AXM * AMX
E - AKM - AMX
C
r rfwDtirv vmtnrrv nv re* i
C
IFIABS(C).IE.EP) TBEN
CX - 1.000
SX - O.ODO
ELSE
COT2X - D/C
SIC - SICNI1.0.COT2X)
COTX - COT2X * (SIC-SQRTI1.0DO » COT2X»COT2X) I
SX - SIC/SORT(1.0DO * COTX'COTX)
CX - SX-COTX
ENDIF
IF(YB.LT.O.ODO) TBEN
TEM - CX
CX - SX
SX - -TEM
ENDIF
COS2X - CX«CX - SX'SX
SIN2X - 2.000'SX«CX
0 - D-COS2X » C-SIN2X
B - B*CO32X - BJ*SIN2X
DEN - G » 2.0DOMI-! * D-D)
TANBY - (I'D - 8/2. 0001 /DEN
C
C COMPUTE ELEMENTS OF [Sil
C
IFIABS(TANBY).LE.EP) TBEN
CBY - l.ODO
SBY • O.ODO
ELSE
CBY - 1.0DO/SQXTI1.0DO - TANBY*TANBY)
SBY - Cfll-TANBY
ENDIF
C
C — ••• COHPUT& KUMENTS OF [Ti] • [Ri} [Si]
C
Cl - CBY«Ot - SBY'SX
C2 - CBY'CX » SBY>SX
SI - CflY'SX •• SBY*CX
82 - -CBY'SX •>• SBY-CTt
C
C APPLY TRANSFORMATION IF MARRANTED
C
IFUABS(Sl) .CT.EP) .OR. (ABS(S2) .CT.EP) ) TBEN
HARK - .IALSE.
DO 23 I-l.N
AKI - AIK. I)
AMI - AIM, I)
A|K, I) - Cl-AKI * SI-AMI
AIM, I) - S2*AKI + C2«AMI
IF (LEFT) TBEN
TKI •• TIX.I)
TKI - T(M,I)
T(K,I) - C1»TKI + S1«TMI
TIM, I) - S2-TKI * CJ*TMI
ENDIF
2] CONTINUE
DO 24 I-l.N
AIK - JMI.K)
AIM - AII.M)
AII.K) - C2-AIK - S3 -AIM
AII.M) - -Sl'AIK + Cl'AIM
IF(RICBT) TBEN
TIK • T(I,K)
TIM • T(I,M)
Til, I) - C2-TIK - SJ'TIM
TII.B) - -Sl-TIK + Cl-TIM
ENDIF
24 CONTINUE
ENDIF
25 CONTINUE
2t CONTINUE
TMX - 90
RETURN
END
Append-29
-------
NBS: Indoor Air Quality Model
Phase II Report
Appendix • FORTRAN 77 Source Code
Include Rles:
REALM EP
COMMON /CNTCOM/NFNOD.NPEON.MFBAN.NFELM.EP.
*M7V, KPF, MFC. MFC. MPKEQ. MPWS
C CALSAPX ARRAY MANAGEMENT
lARYCOM.INC
COMMON /ARYCON/ NUMA.NEXT. IDIR,IP»>
C VARIABLE-
-DESCRIPTION-
C HIOT SIZE OP BLANK COMMON VECTOR IA
C HI CURRENT DATA TYPE: 1-INTECER: 2-REALi 3-CHAR.
C IA(MTOT) BLANK COMMON VECTOR
C NUMA NUMBER OF ARRAYS IN BLANK COMMON DATA BASE
C NEXT NEXT AVAILABLE STORAGE LOCATION IN BLANK COMMON
C IDIR START OF DIRECTORY IN BLANK COMMON
C IPO) NUMBER OF BYTES IN INTEGER, REAL, CHARACTER DATA
C
C CALSAPX I/O FILE MANAGEMENT
NFEQN
NFBAN
NFELM
EP
POINTERS
—POINTER—ARRAY-
INTEGER LFNAME
LOGICAL ECEO.EOD.EOC
CHARACTER*! FNAME*12, EXT*], MODE'S
COMMON /IOCC*a/NTR,NTN,NCMD.NIN.NOT,HDl,ND2,ND3,ND4,
+LFNAME,ECHO,ZOO,EOC
COMMON /ICCON2/ NODE,EXT,FNAME
C VARIABLE DESCRIPTION
C /IOCOM/
C NTR LOGICAL UNIT NUMBER FOR TERMINAL-READ (KEYBOARD)
C NT* LOGICAL UNIT NUMBER FOR TERMINA-WRITE (SCREEN)
C NCMD LOGICAL UNIT NUMBER FOR COMMAND/DATA INPUT
C NIN LOGICAL UNIT NUMBER FOR INPUT DATA ASCII FILE
C NOT LOGICAL UNIT NUMBER FOR OUTPUT DATA ASCII FILE
C ND1 thru ND« LOGICAL UNIT NUMBERS FOR GENERAL USE
C FNAME*12 RESULTS OUTPUT FILE NAME
C LFNAMS LENGTH OF FILENAME KITS TRAILING BLANKS REMOVED
C EXT*J RESULTS OUTPUT FILE EXTENSION
C MODE COMMAND HODSt -INTER••INTERACTIVE, 'BATCH'-BATCH
C ECBO WHEN .TRUE. ECBO RESULTS OUTPUT TO NTH (SCREEN)
C BOD EHD-OF-DATA LOGICAL
C EOC EMD-OF-COMMAND LOGICAL
C
C VARIABLE-
C ERR
C
C
C
C
C
C
C
C
c-
c
C
C
C
C
C
C
C
C
C
C
c-
-DESCRIPTION-
DO-NIILE TERMINATOR FLAG
NUMBER OF FLOW SYSTEM NODES
NUMBER OF FLON SYSTEM EQUATIONS
- NFNOD (CURRENT VERSION)
(HALF) BANDWIDTH OF FLON SYSTEM EQUATIONS
NUMBER OF FLON ELEMENTS
MACHINE PRECISION
TO BLANK COMMO!
LOCATIONS
HPV V(NSNOD) I VOLUMETRIC MASSES
MOT F(NPEQN.2*MSBAN-1)I FLOW MATRIX (UNSYMMETRIC)
MPC C (NFEQN) > CONTAMINANT CONCENTRATION
MPG G(NFEON) i CONTAMINANT GENERATION
MPXZQ KEO(NFNOO) I SYSTEM EQUATION NUMBERS
t 0 • UNDEFINED
t NEC - CONCENTRATION PRESCRIBED DOF
i PCS - GENERATION PRESCRIBED DOF
i ELEMENT MASS FLO* RATES
MPVE ME (NFELM)
C CALSAPX FREE-FIELD
c
IFRECON.INC
CHARACTER LtNS*l. LLINE*1«0
/CLIMS1/ LTNI(ICO)
/CLINX2/ II,JJ
EQUIVALENCE (UINE.LINEd))
SAVE /CLINEl/,/CLINE2/
C LINE (1601
C II
C JJ
c
C CALSAPX CONN
COMMAND LINE BUFFER
END-OF-OATA IN LINE BUFFER
END-OF-INFORMATION IN LINE BUFFER
INFORMATION - DATA 1 COMMENTS
AND MANAGEMENT (CMDCOM.INC
CBARACTER'l NCMND,M1,M2,M3,M4,MS,MC,K>,NNCMHD'<
COMMON /CMND/ NCMND(8).MIf4>.K2(4).MJ(4>,m«>,MS(4),MC(4),M1(4)
EQUIVALENCE (NNCMND.NCMND(l))
-VARIABLE-
-OESCRIPTIOI
NCMND(«)>1 CURRENT COMMAND
NNCKND-B CURRENT COMMAND
to M7«) CURRENT ARRAY NAMES
C CONTAM COMMON STORAGE
SCNTCOM86.INC
Append-30
-------
NB5-114A IREV. 2-eo
u.s. oe.fr. OF COHM.
BIBLIOGRAPHIC DATA
SHEET fSee instructions)
1. PUBLICATION OR
REPORT NO.
NBSIR 87 3661
2. Performing Organ. Report No.
3. Publication Date
4. TITLE AND SUBTITLE
Indcfcr Air^Quality Modeling - Phase II Report
Residential Indoor Air Quality Simulation
5. AUTHOR(S)
James W. Axley
6. PERFORMING ORGANIZATION (If joint or other than NBS. see instructions;
NATIONAL BUREAU OF STANDARDS
DEPARTMENT OF COMMERCE
WASHINGTON, D.C. 20234
7. Contract/Grant No.
1. Type of Report & Period Covered
9. SPONSORING ORGANIZATION NAME AND COMPLETE ADDRESS (Street. City. State. ZIP)
U.S. Environmental Protection Agency U.S. Department of Energy
1000 Independence Ave., SW
Washington, DC 20585
Washington, DC 20460
10. SUPPLEMENTARY NOTES
| | Document describes a computer program; SF-I8S, PIPS Software Summary, is attached.
11. ABSTRACT (A 200-word or less factual summary of most significant information. If document includes a significant
bibliography or literature survey, mention it here)
This interim report presents the results of Phase II of the NBS General Indoor Air
Pollution Concentration Model Project. It describes; the theoretical basis of a general-
purpose nonreactive contaminant dispersal analysis model for buildings, the computation-
al implementation of a portion of this model in the program CONTAM86, and examples of
the application of this model to practical problems of contaminant dispersal analysis.
Presently the model is being extended to handle problems of reactive contaminant dis-
persal analysis and full computational implementation of all portions of the model is
being completed.
The contaminant dispersal analysis model is based upon the idealization of building
air flow systems as an assemblages of flow elements connected to discrete system nodes
corresponding to well-mixed air zones within the building and its HVAC system. Equation
governing the air flow processes in the building (e.g., infiltration, exfiltration, HVAC
system flow, and zone-to-zone flow) and equations governing the contaminant dispersal
due to this flow, accounting for contaminant generation or removal, are formulated by
assembling element equations so that the fundamental requirement of conservation of mass
is satisfied in each zone. The character and solution of the resulting equations is
discussed and steady and dynamic solution methods outlines.
12. KEY WORDS (Six to twelve entries; alphabetical order; capitalize only proper names; and separate key words by semicolons)
contaminant .dispersal analysis; flow simulation; building simulation; building dynamics;
computer simulation techniques; discrete analysis techniques.
13. AVAILABILITY
XS2 Unlimited
| | For Official Distribution. Do Not Release to NTIS
| | Order From Superintendent of Documents, U.S. Government Printing Office, Washington, D.C.
20402.
Order From National Technical Information Service (NTIS), Springfield, VA. 22161
14. NO. OF
PRINTED PAGES
156
15. Price
$18.95
USCOMM-OC 6043-PBO
------- |