DEVELOPMENT OF A VISUAL Tool
FOR HVACSIM+
By
JI ADITYA VARANASI
Bachelor of Engineering
Osmania University
Hyderabad, India
1998
Submitted to the Faculty of the
Graduate College of the
Oklahoma State University
in partial fulfillment of
the requirements for
the Degree of
MASTER of SCIENCE
December, 2002
DEVELOPMENT OF A VISUAL TOOL
FOR HVACSlM+
Thesis Approved:
::t:>> ~!d(JVi=::;1__
11
ACKNOWLEDGE:MENTS
This thesis grew out of a series of dialogues with my advisor Dr. Dan Fisher and
my committee member Dr. Ron Delahoussaye.
I am indebted to Dr. Fisher for all his support during slack times. He always
'showed me a way out. His comments on chapter drafts are themselves a course in critical
thought upon which I will always draw.
Through his Socratic questioning, Dr. Delahoussaye brought me closer to the
reality I had initially perceived, eventually enabling me to grasp its rich complexity.
I am grateful to my committee member Dr. Jeff Spitler who was at the center of
our friendships, turning water into wine with his music and incredible imagination.
I am also thankful to the Federal Highway Administration (FHWA) for providing
continuous support for the work.
Many thanks to Xioabing Uu, Ray Young and Mahadevan Ramamoorthy for
providing decisive and energetic support all the way.
I am forever grateful to my parents whose foresight and values paved the way for
a privileged education, and to my brother who gently offered unconditional support at
each tum of the road.
III
TABLE OF CONTENTS
Chapter P.age
1. INTRODUCTION AND LITERATURE REVIEW 1
1.1. Introduction 1
1.2. Thesis objective and scope 2
1.3. Background 3
1.4. ,HVACSIlYI+ 5
1.5. Simulation Methodology 8
2. DETAILED REVIEW OF HVACSIM+ 14
2.1. Introduction 14
2.2. Architecture of HyACSIM+ 14
2.3. Modular Simulation Program (MODS1M) 16
2.3.1. Hierarchical and Modular Approach 17
2.3.2. Hybrid Simulation Time Steps 20
2.3.3. Time Dependent Boundary Conditions 20
2.3.4. Reports Generation 21
2.4. Numerical methods in MODSIlYI 22
2.4.1. Solving Systems of Simultaneous Nonlinear Algebraic Equations , 22
2.4.2. Mathematical description of the solver 25
2.4.2.1. Powell's Hybrid Method 28
iv
Chapter Page
2.4.3. Convergence Criteria in the Equation Solver ~ 33
2.4.4. Integration of Stiff Ordinary Differential Equations 36
2.4.5. Interpolation of data 38
3. HVACSIM+ ENIIANCEMENTS 40
3.1. Development of the userinterface 40
3.1.1. Visual Modeling Tool 42
3.1.2. Residual Reporting ~ 43
3.2. Development of new convergence metrics 47
3.2.1. Design of new convergence metrics 48
3.2.2. Performance evaluation of the new convergence metrics 50
3.2.3. Simulation results for the Smart Bridge system 53
3.2.4. Simulation of Single zone air handler system with a cooling coiL 60
4. TECHNICAL OVERVIEW OF THE USER INTERFACE 64
4.1. Concept 64
4.2. Architecture 64
4.3. Data Format 65
4.4. Configuration 68
v
Chapter Page
4.5. Workspace "' 68
4.6. Boundary File Editor : 72
4.7. Simulation 73
4.8. Plotting 14
4.9. Installation and maintenance 75
4.10. Code Description 76
5. EXAMPLES 82
5.1. Simulation of Examples. from HVACSIM+ documentation 82
5.1.1. Example 1: Inlet pipe and linear valve system 82
5.1.1.1. System Description 83
5.1.1.2. Creating the system modeL 85
5.1.1.3. Preparing the Boundary Value File 94
5.1.1.4. Running the simulation 97
5.1.1.5. Analyzing the simulation results 99
5.1.2. Example 2: Heating Coil Control Loop, Variable Inlet Water Temperature
Control 102
5.1.2.1. Model Description 103
5.1.2.2. Creating the simulation 104
5.1.2.3. Preparing the Boundary Value File..: 106
vi
Chapter Page
5.1.2.4. Running the simulation 107
5.1.2.5. Analyzing·the simulation results 108
5.1.3. Example 3: Cooling Coil Control Loop with Variable Water Flow Rate..... 110
5.1.3.1. Model Description 110
5.1.3.2. Creating the simulation 1I 1
5.1.3.3. Preparing the Boundary Value File 112
5.1.3.4. Running the simulation 113
5.1.3.5. Analyzing the simulation results 113
5.1.4. Example 4: Air handling umt without heat exchangers 115
5.1.4.1. System Description 116
5.1.4.2. Creating the simulation 117
5.1.4.3. Preparing the Boundary Value File 118
5.1.4.4. Running the simulation 119
5.1.4.5. Analyzing the simulation results 119
5.1.5. Example 5: Single Zone Air Handler with Cooling CoiL 120
5.1.5.1. Model Description 121
5.1.5.2. Creating the simulation 122
5.1.5.3. Preparing the Boundary Value File 123
5.1.5.4. Running the simulation 124
5.1.5.5. Analyzing the sim.ulation results . 125
vii
Chapter Page
5.2. Simulation of Smart Bridge Systems 126
5.2.1. Design Procedure Step 1 126
5.2.1.1. Model Description 127
5.2.1.2. Creating the system model.. 128
5.2.1.3. Preparing the Boundary Value File 132
5.2.1.4. RuniIing the simulation 133
5.2.1.5. Analyzing the simulation results 135
5.2.2. Design Procedure Step 2 136
5.2.2.1. Model Description 136
5.2.2.2. Creating the system modeI. 137
5.2.2.3. Preparing the Boundary Value File 140
5.2.2.4. Analyzing the simulation results 141
5.2.3. Design Procedure Step 3 144
5.2.3.1. Model Description 144
5.2.3.2. Creating the system model.. 145
5.2.3.3. Preparing the Boundary Value File 149
5.2.3.4. Analyzing the simulation results 149
viii
Chapter r Page
6. CONCLUSIONS AND RECOMMENDATIONS 151
6.1. Conclusions ~ .' 151
6.2. Recommendations ~ :.: 152
REFERENCES 154
APPENDIX A  DESCRIPTION OF HVACSIM+ MODELS 160
APPENDIX B  SIMULATION TASK TREE 194
APPENDIX C  HVACSIM+ EXAMPLE SIMULATIONS 196
APPENDIX D  SMART BRIDGE SIMULATIONS 219
APPENDIX E  MODIFIED MODSIM CODE 236
IX
..
LIST OF TABLES
Table 1 'u Page
Table 1.31. List of few Building Energy Simulation tools 4
Table 1.41. Table showing the differences between the work done by Yosuke et aI. and
this ,thesis 7
• I
Table 2.3.11. Categories of variables 17
Table 2.3.12. Component models arranged in functional categories 18
Table 2.4.11. List of parameters to SNSQ 23
Table 3.2.31. Results of simulation performed while varying only the error tolerance
applied to the solution vector 56
Table 3.2.32. Results of simulation performed while varying the error tolerance for the
residual vector 58
Table 3.2.33. Results of simulation performed while varying the error tolerances ........ 60
Table 3.2.41. Performance of the simulation at selected time steps 61
Table 3.2.42. Results of simulation performed while varying only the error tolerance
applied to the solution vector 62
Table 3.2.43. Results of simulation performed while varying only the error tolerance for
the residual vector '" 63
Table 5.1.1.31. Example 1: Boundary values for the simulation 95
Table 5.1.2.31. Example 2: Boundary values for the simulation 106
Table 5.1.3.31. Example 3: Boundary values for the simulation 113
x
Table Page
Table 5.1.4.31. Example 4: Boundary values US
,
Table 5.1.5.31. Example 5: Boundary values 123
I
Table 5.1.5.51. Example 5: Temperature response 125
Table 5.1.5.52. Example 5: Control signal response 126
r
"
...
Xl
LIST OF FIGURES
Figure ,n Page
Figure 2.3.11. Flow diagram of programs and data files within HVACSIM+ 15
Figure 2.4.2.11. Powell's hybrid step 30
Figure 2.4.2.12. DOGLEG Algoritlun 32
Figure 3.1.21. Simulation Requirements c¥alog used to initiate residual reporting ........ 45
Figure 3.1.22. Sample residual reporting screen (LOG scale) 45
F~gu~e 3.1.23. Sample Diagnostic information file : 46
Figure 3.1.24. Residual plot generated by the visual modeling tool 47
Figure 3.2.21. Geothermal Smart Bridge (Design Procedure Step 2) 51
Figure 3.2.31. Plot of residuals when the snow event begins 54
Figure 3.2.32. Plot showing the time step (216000 seconds) when the solution did not
converge 55
Figure 3.2.33. Plot of residual for different error tolerance for the solution vector ........ 57
Figure 3.2.34. Plot of the average bridgesurface temperature against time when the error
tolerance for the solution vector is varying 57
Figure 3.2.35. Plot of residual for different error tolerance for the residual vector ........ 59
Figure 3.2.36. Plot of average bridge,surface temperature against time when error
tolerance for the residual vector is varying 59
Figure 4.31. Input file structure 67
Figu~e 4.32. Format of additional infonnation in the raw output file 67
Xli
Figure Page
Figure 4.51. Interface menu and toolbars 1 •• 69
Figure 4.52. Figure shQwing the description of component icons ~ 70
Figure 4.53. Dialog used to set the parameters of units " 71
Figure 4.54. Dialog used to set the properties of boundary conditions 71
Figure 4.55. Error message generated while making a bad connection 72
Figure 4.61. Section of a boundary file editor 72
Figure 4.71. Figure showing the concurrent working of DOS shell and the interface 73
Figure 4.72. Error message while running a deficient simulation 74
Figure 4.81. Forms used for plotting simulation results ,and residuals 75
Figure 5.1.1.11. Example 1: System schematic diagram 82
Figure 5.1.1.21. Example 1: Dialog showing the simulation requirements 87
Figure 5.1.1.22. Example 1: Dialog showing the simulation requirements 88
Figure 5.1.1.23. Example 1: Dialog showing the names for the boundary conditions 90
Figure 5.1.1.24. Example 1: Setting the properties of boundary conditions 90
Figure 5.1.1.25. Example 1: Arrangement of entities in the simulation 91
Figure 5.1.1.26. Example 1: Initial values for unit 1 92
Figure 5.1.1.27. Example 1: Initial values for unit 2 93
Figure 5.1.1.28. Example 1: Setting variable names 93
Figure 5.1.1.29. Example 1: Screenshot showing the reported variables selected for the
simulation , 94
xili
~~re h~
Figure 5.1.1.31. Example 1: Boundary me editor 96
Figure 5.1.1.41. Example 1: Initial screen before the simulation begins 97
Figure 5.1.1.42. Example 1: Screen showing the end of simulation , 98
Figure 5.1.1.43. Ex.ample 1: Partial listing of the Definition file generated during
simulation ~ _ 98
Figure 5.1.1.44. Example 1: Plot of reported variables 99
Figure 5.1.1.51. Example 1: Temperature response 101
Figure 5.1.1.52. Example 1: Control valve response 102
Figure 5.1.21. Example 2: System setup 103
Figure 5.1.2.21. Example 2: Dialog showing the simulation requirement 105
Figure 5.1.2.22. Example 2: Dialog showing the simulation requirements 105
Figure 5.1.2.41. Example 2: Screen showing the simulation requirements 107
Figure 5.1.2.42. Example 2: Plot of reported variables 108
Figure 5.1.2.51. Example 2: Plot between coil outlet air temperature and time 109
Figure 5.1.31. Example 3: System schematic diagram 110
Figure 5.1.3.21. Example 3: Dialog showing the simulation requirements 111
Figure 5.1.3.51. Example 3: Controlled variable (air temperature) response 115
Figure 5.1.41. Example 4: System Schematic 116
Figure 5.1.4.21. Example 4: Simulation requirements dialog 117
Figure 5.1.4.41. Example 4: Plot of reported variables 119
xiv

fl~~ h~
Figure 5.1.51. Example 5: System schematic diagram ~ ( 121
Figure 5.1.5.21. Example 5: Simulation requirements dialog "" 122
Figure 5.2.11. Smart Bridge  Design procedure step 1: System schematic 127
Figure 5.2.1.21. Smart Bridge  Design procedure step 1: Simulation requirements 128
Figure 5.2.1.22. Smart Bridge  Design procedure step 1: Error tolerances ., 129
Figure 5.2.1.23. Smart Bridge  Design procedure step 1: Dialog showing the names for
the boundary conditions 130
Figure 5.2.1.24. Smart Bridge  Design procedure step 1: Screenshot showing the
various units, boundary conditions and reported variables in the simulation 131
Figure 5.2.1.25. Smart Bridge  Design procedure step 1: Initial values for unit 1 131
Figure 5.2.1.26. Smart Bridge  Design procedure step 1: Initial values for unit 2 132
Figure 5.2.1.27. Smart Bridge  Design procedure step 1: Initial values for unit 3 132
Figure 5.2.1.31. Smart Bridge  Design procedure step 1: Boundary file editor 133
Figure 5.2.1.41. Smart Bridge  Design procedure step 1: Residual plot generated by
MODSThi during simulation 134
Figure 5.2.1.42. Smart Bridge  Design procedure step 1: Plot of simulation results 134
Figure 5.2.1.51. Smart Bridge  Design procedure step 1: Plot of simulation results when
heater is operating in mode 1 with the maximum outlet temperature set at 47 degrees
centigrade 135
xv
Figure Page
Figure 5.2.1.52. Smart Bridge  Design procedure step 1: Plot of simulation results when
heater is operating in mode 2 with the maximum outlet temperature set at 47 degrees
centigrade 136
Figure 5.2.2.21. Smart Bridge  Design procedure step 2: Screenshot showing the
various units, boundary conditions and reported variables in the simulation......... 138
Figure 5.2.2.22. Smart Bridge  Design procedure step 2: Initial values for unit 1 138
Figure 5.2.2.23. Smart Bridge  Design procedure step 2: Initial values for unit 2 139
Figure 5.2.2.24. Smart Bridge  Design procedure step 2: Initial values for unit 3 139
Figure 5.2.2.25. Smart Bridge  Design procedure step 2: Initial values for unit 4 140
Figure 5.2.2.31. Smart Bridge  Design procedure step 2: Boundary file editor. 141
Figure 5.2.2.41. Smart Bridge  Design procedure step 2: Plot of simulation results with
maximum number of boreholes set at eight.. 142
Figure 5.2.2.42. Smart Bridge  Design procedure step 2: Plot of simulation results with
maximum number of boreholes set at four 143
Figure 5.2.31. Smart Bridge  Design procedure step 3: System schematic 144
Figure 5.2.3.21. Smart Bridge  Design procedure step 3: Screenshot showing the
various units, boundary conditions and reported variables in the simulation......... 146
Figure 5.2.3.22. Smart Bridge  Design procedure step 3: Initial values for unit I 146
Figure 5.2.3.23. Smart Bridge  Design procedure step 3: Initial values for unit 2 147
Figure 5.2.3.24. Smart Bridge  Design procedure step 3: Initial values for unit 3 147
xvi
NO.MENCLATURE J •
n = Number of scalars in a vector
l(.) = Jacobian matrix
, r .
x • Solution vector
•• I I'
x k Intennediate value of x
1·1 = Absolute value
II, II Euclidean norm
e, = Error tolerance in the equation solver
er = Relative error tolerance
eo = Absolute error tolerance
t g = Time interval
esol vec = Error tolerance applied to a solution vector
eruidlUJl vec = Error tolerance applied to the residual vector
p = Pressure variable
C. = Control signal
T = Temperature variables
W Mass flow rate
xviii
Figure Page
Figure 5.2.3.25. Smart Bridge  Design procedure step 3: Initial values for unit 4 148
Figure 5.2.3.26. Smart Bridge  Design procedure step 3: Initial values for unit 5 148
Figure 5.2.3.27. Smart Bridge  Design procedure step 3: Initial values for unit 6 149
Figure 5.2.3.41. Smart Bridge  Design procedure step 3: Plot of simulation results 150
xvii

I CHAPTERl
INTRODUCTION AND LITERATURE REVIEW
1.1.1ntroductitJn
'. .f
The quest for higher building energy efficiency in the world has encouraged
1 •
engineers to focus on the relationship between design variables and energy perfonnance.
~
Assessment of building energy perfonnance is fundamental in making decisions
regarding energyefficient design of buildings and in quantifying the impact of energy
conservation measures. Evaluation of energy characteristics of buildings serves as a
critical base for developing building energy standards and assessing their effectiveness.
Simulation of building heating and cooling systems has become a very useful tool
in designing energyefficient and costeffective systems. However, as described by
Spitler (2001), system simulation is an essential design tool for groundsource heat pump
(GSHP) systems, particularly hybrid GSHP systems that incoxporate additional heat
rejection components, such as cooling towers, ponds, and heated pavement systems. It is
fortuitous that the same simulation tools useful for designing building heating and
cooling systems are also useful for designing similar thennal systems, such as "Smart
Bridge" (Chiasson and Spitler 2001), which are not connected to buildings.
HVACSIM+, which stands for 'HVAC SIMulation PLUS other systems', is a
nonproprietary simulation package developed at the National Institute of Standards and
Technology (NIST), Gaithersburg, Maryland, U.S.A. It is capable of modeling HVAC
(heating, ventilation and airconditioning) systems, HVAC controls, the building model,
1
and energy management systems. In HVACSIM+ both the specifications by the user and
the internal representation of HVAC elements are represented !in terms of individual
components like fan, duct, heating coil, boiler, pump, pipe, etc., which are connected to
fonn complete systems. This allows the users to develop new models and introduce them
in the package to simulate them in various configurations. This capability has been
extended by the development of component models for ground loop heat exchangers
(Yavuzturk and Spitler 2000), heated pavements (Ohiasson, et aI. 2000a), watertowater
heat pumps (Jin and Spitler 2002, Ramamoorthy 2(01), and ponds (Chiasson, et al.
2000b; Ramamoorthy, et a1. 2001). The modular component based modeling
environment, coupled with the available component models makes this an ideal tool for
design of heated bridge deck systems. Unfortunately, the existing HVACSIM+ user
interface makes the program unnecessarily difficult to use, even for experienced users.
1.2. Thesis objective and scope
This study aims at the development of a Visual Modeling Tool for HVACSJM+.
The study is coupled with the implementation of new convergence metrics in the tool.
The main objectives of this study can be summarized as follows:
• To design, implement and validate a new tool that provides a better userinterface
for HVACSIM+.
• To design, implement and validate new convergence criteria in order to
achieve improved convergence speed and minimized residuals in the outputs.
Chapter 2 of this thesis conducts a detailed review of HVACSIM+ and the various
programs in the package. The main program, MODSIM, is explained in detail. The
numerical methods used in the package are also discussed.
2
Chapter 3 introduces the problems in HVACSIM+ and discusses the solutions
implemented to tackle the shortcomings of the userinterface and the convergence criteria
in the equation solver. This chapter explains the design and implementation of the new
convergence metrics. The performance evaluation results of the new criteria are also
presented in this chapter.
Chapter 4 provides the technical overview of the userinterface in the visual tool.
It discusses the architecture, file formats, configuration and features of the tool. It
presents a detailed description of the source code structure in the visual tool.
Chapter 5 is intended to validate the visual tool by subjecting it to specific test
cases. The examples given in the HVACSIM+ Users Guide are selected for the study.
These example simulations are reproduced using the tool and the obtained results are
compared with the actual output.
1.3. Background
Simulation programs, intended for research, have been around for many years but
their development into usable tools has not kept pace with other software commonly used
by engineers  spreadsheets, CAD, etc. The large gap between what is available in the
market and what is actually used and comprehended by the design professionals suggests
that implementation, rather than just technical capability, is key to increasing the use of
simulation tools and, hence, energy efficiency.
Many buildingenergy analysis programs were developed in USA and Europe.
There are more than two hundred programs available in USA and over a hundred
programs in Europe and elsewhere. Table 1.31 presents a list of few programs currently
being used.
3
ProeraDl Reference Source .~ Country of use
APECESPll Wickham, 1985 USA
ASEAM Ohadi, Meyer, and Pollington, 1989 USA
BESA BESA,1993 Canada
BLAST BLAST,1991 USA
BUNYIP MoHer and Wooldridge, 1985 l\ustralia
DOE2 Birdsall, et a!., 1990; LBL, 1981 USA
EnergyPlus Crawley, et ai, 2(x)() USA
ESPr ESRU,2000 UK
HVACSIM+' Clark,1985 . (I . I USA
SPARK BuhI, et al., 1993 USA
TRACE 600 Trane Company, 1992a & b USA
TRNSYS TRNSYS,1988 USA
Table 1.31. List offew Building Energy Simulation tools
These packages hav~ different levels of ~tail and input requirements. They come
from research institutions, equipment manufacturers and private consulting finns. Some
of the programs are public domain programs and others are proprietary programs. Apan
from energy analysis, many of them also allow for standard HVAC design load
calculations.
Buildingenergy analysis programs have undergone a slow evolution since arrival
nearly three decades ago. The simulation techniques are rapidly changing with decreasing
cost and increasing flexibility of compute, systems. In the 1970's, the cost of conducting
an energy analysis study was high. Energy programs in those «;fays were developed on
mini or mainframe computers and were inaccessible to the vast majority of potential
users. In the 1980's and 1990's, introduction of microcomputers and emergence of
microcomputer versions of the simulation programs made it more affordable and
accessible to carry out detailed energy studies. Many of the current simulation programs
are succeeding generations of some previous ones.
4
Considerable amount of literature is available regarding studies perfonned on
existing energy analysis tools. Some of the important works are presented below.
1.4. HVACSIM+
Park, Clark and Kelly (1985) give an overview of the HVACSIM+ package. The
I
architecture and the main programs in the package are discussed. The paper discusses the
features of MODSIM and introduces the numerical methods implemented in it. The
techniques to solve stiff differential equations and data interpolation are discussed. The
paper also presents a detailed discussion of the building shell and zone models used in
HVACSIM+. The discussion on HVACSIM+ in this paper includes the solving
techniques used, convergence criteria and numerical instability of the methods used in the
equation solver.
Hensen (1995) gives an overview and examples of various approaches to system
simulation in buildings. Advantages and disadvantages of different methods are
described. The author classifies the existing software tools in four levels, as LEVEL AD.
In the case of LEVEL A, design specification and representation of plant systems is
purely conceptual and only the room processes are considered. In the case of LEVEL B,
the specification by the user is in terms of real systems like variableairvolume, constantvolume
zone reheat system, residential wet central heating, etc. Examples of simulation
systems operating on this level are DOE2 and BLAST. In the case of LEVEL C both the
specification by the user and the internal representation is in terms of individual plant
components like fan, duct, heating coil, boiler, pump, pipe, etc., which are connected to
form complete systems. Two main approaches can be distinguished in terms of individual
component models. They are the inputoutput based approach and the conservation
5
equation based approach. The inputoutput based approach repre ents each part of the
system as an equivalent inputoutput relationship. These are connected to comprise the
whole system in such a way that the output from one component is fed into the next as an
input. TRNSYS is a well known example implementing this approach. In the
conservation equation based approach HVAC system modeling is achieved by a modular,
componentwise approach, involving representation of parts of a system by discrete nodal
schemes and by the derivation of energy and mass flow equation sets that are solved
simultaneously for each time step. An ,example of a conservation equation based system
is HVACSIM+. In the case of LEVEL D, the specification is in the same form as in
LEVEL C but the internal representation is based on the category of the variables. ESPr
is an example of a LEVEL D software tool.
Sahlin (1996) reviews many en~rgy analysis tools like TRNSYS, HVACSIM+.
and SPARK and describes the design of a generalpurpose simulation environment, IDA,
and the Neutral Model Format (NMF), a program independent language for modeling of
dynamical systems using differentialalgebraic equations with discrete events. IDA and
NMF are used to effectivdy develop special purpose GUIbased tools that (1) are easy to
use for end users without simulation expertise, and that (2) have good prospects for long
term maintenance and reuse. IDA and NMF aliso serve as a general modeling and
development environment for the sophisticated user.
Yosuke, Xiangyang, and Nobuo (1999) describe the development of a
configuration tool for HVACSIM+. The tool is based on Object Modeling Technique
(OMT), which is used to analyze the information of a HVACSIM+ simulation. The
purpose of the tool is to replace HVACGEN, the preprocessing program of HVACSIM+
6
and provide an integrated environment to run simulation . The software architecture
consists of the system model framework and the user interface framework. These
frameworks are object class group.s for describing the design information, inoluding
connections between components and the hierarohy of components, and the screens
which allow users to edit this infOlmation, respectively. Some of the aspects in which the
work differs from the implementation strategies followed in this thesis are presented in .
the Table 1.41. , I • \
Work done by Yosuke et ala Current thesis study
Programming Uses Object Oriented Uses structured programming
technique Programming technique which is technique which is easy to
difficult to understand and code. understand and code.
Programming
Small;talk Visual Basic
lan2Ua2e .. .. l.
Description Incorporates description of a Single diagram describes the
.of Blocks component connection diagram for system and displays all
each block. This procedure is components and connections
difficult to maintainand between them. This procedure is
understand and causes the users to easy to maintain and understand.
work with a pool of diagrams.
Component New component models can be Component models cannot be
models registered and existing ones can be added or edited.
edited.
Table 1.41. Table showing the differences between the work done by Yosuke et al. and
this thesis.
Sowell and Haves (2000) present a discussion on the efficiencies of methods
employed in solution of building simulation models and provide means for benchmarking
them. Several simulation packages such as HVACSIM+, TRNSYS and SPARK are
compared with each other. Specifically, HVACSIM+ and SPARK. are studied and
comparison analysis is presented in tenns of simulation execution times.
Ramamoorthy (2001) presents a discussion on two energy simulation packages,
TRNSYS and HVACSIM+. The author discusses the structure of the two packages, their
7
relative merits and demerits, and steps to convert component models from TRNSYS to
HVACSIM+. According to the author, the simultaneou nonlinear equation solver is
supposed to obtain a selfconsistent solution for each time step in HVACSIM+, as
compared to a relatively simple successive substitution algorithm in TRNSYS. The userinterface
of TRNSYS is better than HVACSlM+ making it easier to work with. TRNSYS
is more suited for using hourly time steps, whereas HVACSIM+ is suited for both shorter
and longer time steps if the variable step algorithm in HVACSIM+ can. be effectively
exploited.
1.5. Simulation Methodology "
A variety of mathematical algorithms are implemented in existing building energy
analysis tools to solve systems of nonlinear algebraic equations. It is important to study
them in order to develop efficient diagnostic tools, specifically, when applying these
algorithms to nontraditional applications. A few selected works are presented below.
Powell (1970) describes his hybrid method, which combines the features.of
Newton and Steepest descent methods. The method is tested and techniques for the
approximation of the Jacobian are suggested. Readers are referred to Chapter 2 for a
r ,
detailed description of Powell's hybrid method.
Hiebert (1982) performed an evaluation of the existing mathematical software that
solves systems of nonlinear equations. Eight FORTRAN applications were considered for .
the study: COSNAF, BRENT, HYBRD, NSOIA, QN, SOSNLE, ZONE, and ZSYSTM.
These programs implemented three different methods: Brown or Brent, quasiNewton,
and Powell's hybrid. The aim of the work was to evaluate the code written and to find the
relati've merits and demerits of the algorithms implemented. A defined set of test cases
8
was used to evaluate and make conclusions regarding the suitabj.lity of :the software for
specific problems. The following observations were made by the author: .. I
• The initial guess has little effect on the performance of the different code .
• Poor function scaling has essentially no effect on the perfonnance of the
Brown and Brent codes. This is due to the fact that the functions value are
supplied one'at a time.
• The hybrid method has not necessarily been an improvement over the quasiNewton
method, especially when poor variable or function scaling is present.
• Work needs to be done on how to handle poor variable scaling.
Dennis and Schnabel (1983) explain various aspects of solving nonlinear system
of equations in their book,. titled 'Numerical Methods for Unconstrained Optimization
and Nonlinear Equations'. The book is a compilation of works done by several pioneers
in the field. It explains the fundamentals of the calculus behind solving nonlinear
equations and describes the Newton's method and several of its variants in detail. In
addition, Jacobian approximation methods such as the Secant method and the Broyden's
method are discussed in detail. The work, designed to be a text book, mentions examples
and gives brief performance analysis for each method discussed.
Shacham (1990) describes a variable order method for solving nonlinear algebraic
equations. The variable order method is a modified fonn of the Improved Memory
Method (IMM), which was earlier proposed by the same author. The IMM method uses
continued fractions to pass an inverse interpolating polynontial through all the previously
calculated points, to fmd a new estimate of the solution. The modified form allows the
number of previously calculated points, which governs the order of the method, to be
9
changed. The algorithms discussed in the paper focus on the implementation of the
method, the approach to attain global convergence and the stopping criteria. The method
was subjected to ninety test cases and solutions for most of them were obtained in less
than twenty iterations, proving that it is robUSt and reliable. It was also observed that by
selecting more previously calculated points, the Fate of convergence of the method oan be
increased. However, as the order increases, the improvement in the solution is moderate.
BuzziFerraris and Tmnconi (1993) describe an improved convergence criterion
in the solution of nonlinear algebraic equations. The inadequacies in the existing
convergence criteria, which requixethat either the sum of squares or the weighted sum of
squares of the residuals decreases after each iteration, motivate the development of the
new criterion. The new criterion introduces scaling of residuals to overcome the problem
in existing criteria that is related to the fact that, if the residuals have widely differing
orders of magnitudes, they tend to appear in the convergence criterion with incorrect and
unbalanced weights. The improved criterion determines convergence using a modified
merit function which equals exactly to the squared distance between the current values of
the unknown and the solution vector in the case of linear system of equations. The
authors provide a geometrical interpretation of the criterion which shows that the method
introduces natural scaling or weighting factors for the residuals in the objective function,
although it cannot be proven to guarantee global convergence. In addition, the evaluation
of the criterion is shown to be computationally inexpensive. The authors end their
discussion recommending the implementation of the criterion in generalpurpose
equation solvers.
10
Paloschi (1994) describes a hybrid continuation algorithm to olve algebraic
nonlinear equations that combines Newton or QuasiNewton and Continuatlion methods
in such a way that the continuation trategy will be used only when the other convergence
criteria fail. The continuation strategy is used only as a last resort because it is
computationally expensive compared to QuasiNewton method. The author refers to
several papers and indicates that the existing methods lack in robustness and proposes the
use of a continuation or homotopy method that is specially designed to improve
robustness. The strategy is governed by the following homotopy function
h(x,H) =H f(x) + (1  H) g(x)
where g(x) is a function such that g(x) = 0 has a known solution. It is selected as g(x) =
f(x)  f(Xt). H is the homotopy parameter that allows defining a family of functions where
h(x,H) is solved easily and hex,] )=f(x). The equation h(x,H) = 0 defines a path x(H)
which can lead to the solution. Three hundred and twenty seven examples that cover
many different problems and reflect different scaling conditions are used to the test the
method. It is found that the algorithm solves almost every problem and is proved to be
robust. Finally, the author compares the perfonnance of the method with the HYBRD
(MINPACK) code, which is the previous version of the routine SNSQ used in MODS1M,
and finds out that the HYBRD code is capable of solving only two hundred test cases.
Schonwalder and Sans (1995) present a method for the steady state analysis and
optimization of nonlinear autonomous electrical circuits having unknown final oscillation
periods in their steady state response. The approach involves discretizing the response
and solving the associated system of nonlinear algebraic equations. Most iterative
algorithms prove useless since autonomous circuits have infrnite solutions leading to
11
difficulties in convergence. The authors discus the use of two robust solution
approaches, namely, the Newton's method and the Fast Simulated Diffusion (FSD)
method. Since the Newton's method is known to be locally quadratically convergent, the
authors also suggest the use of line searoh algorithms or trust region methods that help in
maintaining the Newton's directions or the radii of regions around the solutions, which
decide the maximum step lengths. Better. performance can be attained from. trust region
methods since they provide multidimensional solution search. The motivation behind the
'use of FSD method can be attributed to the fact that many techniques reject any tentative
steps that result in an increase in the associated objective function, thus, making it
difficult for the techniques to escape from local minima. The FSD method generates two
kinds of steps called the descendent steps, which help find local minima and random
steps, which ensure that the whole design space may eventually be reached.
Appropriately reducing the length of.random steps to zero, using simulated annealing
technique, ensures that the algorithm will converge. The authors fmally present
application examples to validate the proposed approach.
Spedicato (1997) investigates the performance of several Newton like methods,
namely Newton's method, the ABS Huang method, the ABS row update method and six
QuasiNewton methods. Thirty one families of problems with dimensions n=10, 50, 100
and two starting points are selected to test the algorithms. Newton's method appears to be
the best in terms of number of problems solved, followed closely by the ABS Huang
method. In addition, Newton's method has the smallest nonconvergence region. The
performance of the QuasiNewton can be summarized as follows:
• Broyden's method and Greenstadt's second method show very poor performance.
12
• Greenstadt's first method, Martinez' column method, nonlinear Huang's method
perfonn similarly.
• Thomas' method appears to be marginally more robust and fast and provides a
better approximation to the Jacobian. In addition, the Thomas' method shows the
smallest nonconvergence region.
MacMillan (1999), in his doctoral thesis, suggests a relaxed convergence criterion
to improve the convergence rate. The existing convergence criteria are compared to find
out the reasons for their failure. It is shown that relaxing the conditions required to prove
global convergence can improve the performance of an algorithm. It is also shown that
minimizing an estimate of the distance to the minimum relaxes the convergence
conditions in such a way as to improve an algorithm's convergence rate. A new
linesearch algorithm based on these ideas is presented that does not force a reduction in
the objective function at each iteration, yet it allows the objective function to increase
during an iteration only i~ this will result in faster convergence. Perfonnance of the
algorithm on some standard test functions is presented to illustrate that it is well defined
and globally convergent.
Gould, Orban, Sartenear and Toint (200 1) present a parameterized variant of the
Newton's method. The Newton's method is perturbed by adding a tenn with a scalar
parameter that is driven to zero as the iteration proceeds. The exact local solutions to the
perturbed systems then fonn a path leading to a solution of the original system, the scalar
parameter determining the progress along the path. It is shown that asymptotic
convergence rate may be obtained in finding the solution and the residual vectors.
Numerical experiments are used to validate the method.
13
CHAPTER 2
DETAILED REVIEW OF HVACSIM+
2.1. Introduction
The HVACSIM+ package has been developed primarily as a research tool for
complete building system modeling. It is written in ANSI standard Fortran 77 and has
enough flexibility to perform simulations of HVAC components, control systems,
building models, or any combination thereof.
Documentation for HVACSIM+ consists primarily of three publications: a
Reference Manual (Clark, 1985), a Users Guide (Clark and May, 1985) and a Building
loads calculation Manual (park, Clark and Kelly, 1986). The Reference Manual explains
the component models distributed with HVACSIM+ and accompanying support routines.
The Users Guide describes the procedure to build and run simulations. It also presents six
examples of system simulations using HVACSIM+. The Building loads calculation
Manual describes component models and support routines and presents sample
simulations for building loads calculation. A description of the numerical methods used
in HVACSIM+ is also discussed in this document.
2.2. Architecture ofHVACSIM+
Operations performed by INACSIM+ can be categorized as preprocessing,
simulation and postprocessing. Separate programs in HVACSIM+ perform these
operations. Figure 2.21 shows how HVACSIM+ programs and data files interact.
14
. Preproces~g 
Thermal
Properties of
Building
Weather Data Materials
File
RDTAPE I Weather Tape l Weather data·
CRWDATA
.. _,
~'l........., Boundary Value File & i
Simulation Control Data ,
,
Simulation
.:_ Summary. Raw output and Initialization files ,
: • • I
Postprocessing                  .    .  .             :
,
Raw output
filE'
,,
:__.__._~
I .
• to
Figure 2.3.11. Flow diagram ofprograms and dlJtafiles within HVACSIM+
The programs that perfonn pre"processing operations are HVACGEN,
SLIM:CON, RDTAPE, CRWDATA and CTFGEN. The main program MODSIM takes
input from the model definition file that is created by SLIM:CON from a work file
generated by HVACGEN. HVACGEN is a frontend program supplied with HVACSIM+
that is used to edit a work file interactively. In generating the simulation work me,
HVACGEN employs a data file containing component model infonnation. Programs
RDTAPE and CRWDATA are intended to read raw weather files and to create weather
data flies in a fonn recognized by MODSIM. The RDTAPE program can recognize
unprocessed weather data files only if they are in one of the following fonnats:
• NOAA Test Reference Year (TRY)
• NOAA Typical Meteorological Year (TMY)
• NOAA SOLMET
15
I ~ I • Weather Year for Energy Calculations (WYEC)
The first three fonnats were proposed by National Oceanic and Atmospheric
Administration (NOAA). If raw weather data is not available, the CRWDATA program
can produce a designdayweather data file. T,be CTFGEN program generates the
conduction transfer functions of multilayered building elements; such as roofs and walls.
The MODSIM program is the main program in HVACSIM:+. It consists of a main
driver program and many subprograms that perform the following functions:
I'J • • Input/output operations
• Block and state variable status control
• Integration of stiff ordinary differential equations
• Solving of a system of simultaneous nonlinear algebraic equations
MODSIM reads the model definition and boundary data files and executes the
simulation with user specified parameters. After a successful simulation, three data files
are generated. These are the summary, raw output and initialization data meso
The postprocessing program, SORTSB, is used for sorting the raw output data.
The output of this program may then be used for plotting with a usersupplied graphic
routine.
2.3. Modular Simulation Program (MODSIM)
MODSIM stands for MODular SIMulation. Since its inception in 1983 at the
University of Wisconsin Solar Energy Laboratory, MODSIM has been enhanced
significantly to include new algorithms for performance and efficiency. Important
features of the current version of MODSIM are described below.
16
2.3.1. Hierarchical and Modular Approach
MODSIM consists of many subroutines called TYPE ubroutines. These
subroutines represent system components that the program can recognize. When a
simulation is constructed, the components modeled by the TYPE subroutines are linked
together. A UNIT in HVACSlM+ represents a component model of HVAC systems or
controls, or a building element. To the user, a unit is presented as a black box that takes a
set of inputs and produces a set of outputs. Using a modular approach, each UNIT is
modeled in the subroutine TYPEn, where n is the index number assigned to the specific
component. Components are also given UNIT numbers so that similar types can be
distinguished. Each unit has distinct inputs, outputs and parameters. Each input and
output variable must belong to a known category. Table 2.3.11 shows a list of available
categories, their index numbers and default units.
Index number Cate20ry name Catel!ory units
1 Pressure kPa
2 Flow Rate kws
3 Temperature °C
4 Control Signal Nondimensional
5 Rotation rate Revolutions per second
6 Ener~y .kJ
7 Power kW
8 Humidity kg of waterlkg of dry air
Table 2.3.11. Categories of variables
A library of component models, suitable for dynamic simulation, was developed
for use with HVACSIM+. This library is called Types component library. It consists of
thirty TYPE subroutines for models of HVAC components and controls. Table 2.3.12
shows the list of models categorized based on their functionality.
17
Functional cateeory TYPE Model Description
Fan or Pumps TYPE 1 Fan or Pump
Conduits TYPE 2 Conduit (Duct or Pipe)
TYPE 3 Inlet Conduit (Duct or Pipe)
Flow Splits and Mer~es TYPE 4 Flow Mer~e
TYPE 6 Flow Split
TYPE 13 Threeway Valve with Actuator
TYPE 17 Mixing Dampers and Merjte
TYPE 18 Plenum
TYPEll 'Grounded' Flow Split
Valves, Dampers and Flow Restrictors TYPES Damper or Valve
TYPE 9 Linear Valve with Actuator
TYPE 13 Threeway Valve with Actuator
I', TYPE 17 Mixing Dampers and Merge
TYPE 23 Steam Nozzle
TYPE 24 Ideal Gas Nozzle
Sensors and Controller TYPE 7 Temperature Sensor
TYPE 8 ProportionalIntegral Controller
TYPE 16 'Sticky' Proportional Controller
TYPE 20 Hijth or Low Limit Controller
TYPE 26 Control Signal Inverter
Heat Exchan,gers TYPE 10 Hot Water to Air Heatin~ Coil (Simple)
TYPE 11 Hot Water to Air Heating Coil (Detailed)
TYPE 12 Cooling or Dehumidifying Coil
TYPE 25 Steam to Air Heat Exchanger
Humidifiers and Dehumidifiers TYPE 12 Cooling or Dehumidifying Coil
TYPE 14 Evaporative Humidifier
TYPE 22 Steam Spray Humidifier
Conditioned Zone TYPE IS Room
Superblock Couplin,g Components TYPE 16 'Sticky' Proportional Controller
TYPE 19 Flow Balance Control
Building Shell TYPE 50 Zone Envelope
TYPESI Building Surface
TYPES3 Weather Input
Building Zone TYPE 52 Zone Model
Smart Bridge deck heatinj]; system TYPE 700 Slab Model: Explicil Slab Model
TYPE 711 Heat Pump: Gang Of Water To Water Heat Pump
TYPE 712 Heal Pump: Single Water To Waler Heat Pump
TYPE 721 GLHE Model: Hierarchical Load Agltlegation Model
TYPE 730 Flow Control: Linear Proportional Control (Heating)
TYPE 740 Heater: 2Mode Ideal Steady Stale Heater
TYPE 750 Pump Model: Simple Pump Model
Table 2.3.12. Component models arranged in functional categories
Simulating very large systems involving many algebraic and differential
equations takes considerable time and greater effort by the routines that implement the
numerical methods to solve them. This problem is tackled by grouping sets of equations
together into smaller subsets. MODSlM[ uses a hierarchical simulation structure
consisting of BLOCKS and SUPERBWCKS.
18
A number of functionally related units (or a single unit) form a block. The
connections between the units define a set of simultaneous equations that are solv~d by
the equation solver in MODSlM. Input variables that do not participate in the set of
equations are boundary conditions. "
A number of blocks (or a single block) make up a superblock. Superblocks (or a
single superblock) comprise a simulation. Connections between the blocks in a
superblock define a set of simultaneous equations that are solved by the equation solver.
MODSIM assumes that the couplings between, superblocks in a simulation are weak
enough so that each superblock can be treated as an indepeij.dent subsy,stem. No
simultaneous equations ~ defined between superblocks. Since superblocks evolve
independently in time, the inputs to a superblock may change at a time when that
superblock is not scheduled for processing. To ease this problem, a superblock input
scanning option is provided. When this option is selected, all superblock inputs are
scanned after each time step. If the inputs to a superblock that was not processed during
the time step have changed, the superblock is called and its state is calculated based on
the new values of its inputs. When the input scanning option is not selected, each
superblock is called at time determined by the timestep control algorithm, regardless of
any changes in its inputs at intennediate times.
In MODSIM, when state variables reach steady state, they are frozen, and
removed from the system of equations until deviations in the steadystate values are
encountered. Similarly, a block is deactivated if all the input variables to the block are
frozen and activated back when any of them change. The user is allowed to enable or
disable variable freezing in the simulation.
19
2.3.2. Hybrid Simulation Time Steps
Many ~componentmodels for HVAC and control systems involve differential
equations. When thesystem is unsteady, a large time step invites numerical instability,
especially when the components in the system have a wide range of time constants. To
prevent this instability, smaller intervals are necessary at simulation startup or dUring a
period when sudden change occurs. After the system has stabilized, the use of long time
steps is more appropriate.
The MODSIM program incorporates two different types of time steps, namely,
the fixed time step and the variable time step. While defining the simulation, the user
chooses the minimum and maximum time step values. If these values are different, it
indicates that the simulation will use variable time steps.
Variable time steps are best suited for solving stiff differential equations that have
widely varying time constants. Building models may use fixed time steps since they are
based on unifonnly distributed time sampling.
MODSIM employs variableorder Gear algorithms (Gear, 1971) that are proven
to work efficiently in integrating stiff differential equations. The order of integration and
the size of the time steps are dynamically detennined during the calculation to minimize
the computation required for a desired degree of accuracy. The time step and integration
order are determined independently for each superblock since they are considered to be
independent subsystems that proceed independently in time.
2.3.3. Time Dependent Boundary Conditions
A state variable that is external to the system being simulated can be designated as
a boundary variable when the simulation work file is generated. The boundary variables
20
may be constant or time dependent. Data for timedependent boundary variables are
stored in the boundary value file and read as the simulation progresses.
The first column in the boundary value file must always contain time values.
Time intervals in the boundary value file need not be equal. since a third order
Lagrangian interpolation method. explained in the Section 2.4.6, is used to find the
boundary variable values at simulation time steps. The remaining columns contain values
of boundary variables at the specified times. MODSIM begins all simulations at time zero
.regardless of the first time entry in the boundary value file.
Including in the boundary data file two different data values of a boundary
variable at a given time signals a discontinuity in a boundary condition. When the change
in a boundary variable is discontinuous, the routine that integrates the differential
equations is reset at the time of discontinuity to bring the simulation time step to a
minimum value.
2.3.4. Reports Generation
MODSIM produces three files on output, namely the Final State File (.FIN), the
Output File (.OUT) and the Summary File (.SUM). The Final State File contains numbers
defming the state of the simulation at the simulation stopping time. The Output File
contains the raw data at each time step in the simulation. The Summary File contains a
summary of the simulation including information relating to configuration, diagnostics
performed and listing of reported variables at equal time intervals.
21
2.4. Numerical methods in ltfODSIM
The numerical methods employed in the MODSIM program involve techniques
for solving systems of simultaneous nonlinear algebraic equations, integrating stiff
ordinary differential equations, and interpolating data sampled in either a fixed period or
variable time intervals. A large number of subprograms in MODSIM are related to these
numerical algorithms.
2.4.1. Solving Systems ofSimultaneous Nonlinear Algebraic Equations
The nonlinear equation solver is implemented in the subroutine SNSQ. The
purpose of this subroutine is to find a minimum of a system of N nonlinear functions in N
variables using the Powell's hybrid method.
The parameters to SNSQ are classified as input and output. Input parameters must
be specified on entry to SNSQ and are not changed on exit. Output parameters need not
be specified on entry and are set to appropriate values on exit from SNSQ. Table 2.4.11
presents a list of input and output parameters and their descriptions. Apart from the listed
parameters, the SNSQ routine also takes the name of the usersupplied subroutine that
calculates the functions as a parameter.
The user has the option of either providing a subroutine that calculates the
Jacobian or of allowing the code to calculate it by a forwarddifference approximation
method. The input parameter roPT is used to identify which method the user prefers to
use. A numerical value of one indicates that the user wishes to use the former method. In
this case, the name of the subroutine that calculates the Jacobian should also be passed to
SNSQ as a parameter.
22
Parameter Type Input! Output Description
IOPT Integer Input IOPT pecifie how the Jacobian will be calculated.
IfIOPT = 1, then the u er must supply the Jacobian
through the subroutine JAC. If IOPT = 2, then the
code will approximate the Jacobian using the forward
i difference method.
N Integer Input It is a positive integer variable set to the number of
functions and variables.
, X Real Input and It is an array of length N. On input; X must contain an
Output initial estimate of the solution vector. On output, X
contains the final estimate of the solution vector.
FVEC Real Output FVEC is an array of length N that contains the
functions evaluated at the output X.
FJAC Real Output FJAC is an NbyN array that contains the orthogonal
matrix Q produced by the QR factorization of the final
approximate Jacobian.
LDFJAC Integer Input LDFJAC is a positive integer variable not less than N
that specifies the leading dimension of the array FJAC.
XTOL Real Input It is a nonnegative number. Tennination occurs when
the relative error between two consecutive iterations is
at most XTOL. Therefore. XTOL measures the
relative error desired in the approximate solution.
MAXFEV Integer Input It is a positive integer variable. Termination occurs
when the number of calls 'to FCN is at least MAXFEV
by the end of iteration.
ML Integer Input ML is a nonnegative integer variable that specifies the
number of subdiagonals within the band of the
Jacobian matrix. If the Jacobian is not banded or
IOPT=I. set ML to at least N  1.
MU Integer Input MU is a nonnegative integer variable that specifies the
number of superdiagonals within the band of the
Jacobian matrix. If the Jacobian is not banded or
IOPT=1, set MU to at least N  I.
EPSFCN Real Input It is used in determining a suitable step for the forwarddifference
approximation. This approximation
assumes that the relative errors in the fu.nctions are of
the order of EPSFCN. IfEPSFCN is less than the
I machine precision, it is assumed that the relative errors
in the functions are of the order of the machine
precision. If IOPT. = 1, then EPSFCN can be ignored.
DrAG Real Input It is an array of length N. IfMODE = I (see below),
DIAG is internally set. If MODE = 2, DIAG must
contain positive entries that serve as implicit scale
factors for the variables.
Table 2.4.11. List ofparameters to SNSQ
23
Parameter Type Input! OutPut ~riptiOD
MODE Integer Input It is an integer variable. IfMODE =I, the variables
J [;"
will be scaled internally. IfMODE =2, the scaling is
I' I specified by the input DIAG. Other values ofMODE
are equivalent to MODE =1.
FACTOR' Real Input It is a positive variable used in determining the initial
step bound. This bound is set to the product of
FACTOR and the Euclidean nonn of DlAG*X if
nonzero, or else to FACTOR itself. In most cases,
FACTOR should lie in the interval (0.1,100.0). The
I generally recommended value is 100.0.
NPRINT Integer Input It is an integer variable that enables controlled Printing
of iterates if it is positive. In this case, FCN is called
with IFLAG =0 at the beginning of the first iteration
and every NPRlNT iteration thereafter and
 immediately prior to return, with X and FVEC
available for printing. Appropriate print statements
, must be added to FCN. If NPRINT is not positive, no
special calls of FCN with IFLAG =0 are made.
NFEV Integer Outeut It is an integer variable set to the number of calls to
FCN.
NJEV Integer Output It is an integer variable set to the number of calls to
JAC. (If IOPT=2, then NJEV is set to zero.)
LR Integer 'Input It is a positive integer variable not less than
(N*(N+ I »12.
QTF Real Output It is an output array of length N which contains the
vector (QT)*FVEC.
INFO Integer Output It is an integer variable. If the user has terminated
execution, INFO is set to the (negative) value of
IFLAG. Otherwise, INFO is set as follows.
• INFO = 0, improper input parameters.
• INFO =I, relative error between two consecutive
iterates is at most XTOL.
• INFO = 2, number of calls to FCN has reached or
exceeded MAXFEV.
• INFO = 3, XTOL is too small. No further
improvement in the approximate solution X is
possible.
• INFO =4, iteration is not making good progress, as
measured by the improvement from the last five
Jacobian evaluations.
• INFO =5, iteration is not making good progress, as
measured by the improvement from the last ten
iterations.
R Real Output R is an output array of length LR that contains the
I upper triangular matrix produced by the QR
factorization of the final approximate Jacobian, stored
rowwise.
WAI, WAl, Real Output • They are work arrays oflength N.
WA3, WA4
Table 2.4.11. (Contd.) List ofparameters to SNSQ
24
2.4.2. Mathematical description ofthe solver
A detailed description of the solver was given by Powell (1970); This section
explains the mathematical description of the methods implemented in the'SNSQ
subroutine. . !
The system of nonlinear equations can be written in vector form as:
(2.4.21 )
Expanding the system in a Taylor series, and neglecting the high order tenns, the
linearized, approximate system bepo~es:
(2.4.22)
where, x' is the solution vector and J (xk
) is a Jacobian evaluated at xlc
. • .
The Newton step of the nonlinear system, &, can be expressed as:
(2.4.23)
The major advantages of Newton's method are its rapid (quadratic) convergence
and its suitability with highly interacting systems because interactions between variables
are fully taken into account by the Jacobi~ matrix. However, the algorithm requires a
good initial guess for successful convergence. In addition, more computational effort is
required to calculate the inverse of the Jacobian matrix. In order to reduce the
computational effort, a QuasiNewton method (Dennis and Schnabel, 1983) is used. 'The
main idea of this method is to approximate the Jacobian at each iteration using only the
current and previous values of x and f(x).
The QuasiNewton method uses Broyden's rankone update (Broyden, 1965)
instead of calculating the full Jacobian at each iteration. For the first iteration, the
Jacobian is calculated either by a usersupplied subroutine or by a forwarddifference
25
..
approximation. The Jacobian is not recalculated until the rankone method fails to give
satisfactory progress. The method is explained below.
Let the k+1 step give the solution for the problem, in the Equation 2.4.21.
Equation 2.4.22 can then be modified as:
f (X k
+
1) =f (xk) +J (x k )(Xk+1 x k) =0
Rearranging the tenns,
J(xk)(Xk+
1
 x k ) = f(xk
)  f(X k+1
)
Let Ak represent the approximation of J(xk
). This equation can be rewritten as:
Ak(xk+[ _xk ) = f(x k) _ f(X k+1)
Equation 2.4.26 is referred to as the Secant equation.
Revisiting Equation 2.4.22,
f(x·)=f(xk)+J(Xk)(x· xk)=O
Rewriting the Equation 2.4.27 in a model form gives,
f(x) = f(x k
) +At (x x k )
Let mk (x) represent the right side of the Equation 2.4.22. We have,
mk (x) =f (xk
) + Ak (x  x k
)
(2.4.24)
(2.4.25)
(2.4.26)
(2.4.27)
(2.4.28)
(2.4.29)
For the next iteration, Xk+l, in the above equation can be modified as:
mk+l (x) =f(X k+1 ) + Ak+L(x _ Xk+I) (2.4.210)
Taking the difference between the equations 2.4.210 and 2.4.29 and rearranging gives,
mk+1(x) _mk(x) = f(X k+1) f(x k ) _ Ak+l(Xt+1_ x k)+ (Ak+J  Ak)(x x k) (2.4.211)
Substituting the Equation 2.4.26 in 2.4.211 and rearranging gives,
26
mJc
+1(x)  m" (x) = (A"+1  A" )(x  x") , l (2.4.212)
Considering the Equation 2.4.212 as a minimizing function and defining step, x, to be
taken in order to reach the global minimum, as
xx" =a(xH1_xk.)+t (2.4.213)
where, a is a small value. Variable t is defined by the equation t T
Sc = O. Substituting
•
the Equation 2.4.212 in Equation 2.4.212 gives,
(2.4.214)
The behavior of the fIrst term on the right side in the Equation 2.4.214 cannot be
predicted since it involves f(x). However, the second term can be made zero by
choosing Ak+1 such that (AIe +1  Ak )t =0 for all t orthogonal to (x Jc
+
1
 x le
). This is
possible only when (AHI  AIe) is in the form u(xle
+
1
 xk l , U E ~n • This can be
represented as:
(2.4.215)
Applying the new form of (A k+1  Ale) to the secant equation in the Equation 2.4.26, we
have,
(2.4.216)
The above equation gives an expression for u. Substituting the Equation 2.4.216 in the
Equation 2.4.215, we obtain,
rP( k) f( Ie+l)l Jc+l k)T
Ak+I=A" +If x  x JX x
(xle+1 _ Xk)T (Xk+l _ xk )
(2.4.217)
The Equation 2.4.217 gives the Broyden's update. The word 'update' indicates that
l(Xk+1) is not updated from scratch, rather the approximation Ale to l(x") into an
27
approximation Ak+l. to J (X
k
+
l
). For the first iteration, A0 is chosen to be eitherthe
negative identity matrix or the finite difference method is applied to the Jacobian. matrix
and the matrix is inverted. I I
The QuasiNewton method onlt requires a single pass through the system at each
iteration. It accurately reflects the interaction of variables and does not need to perform
inversions of matrices at each step. The calculation of second derivatives is also
eliminated.
The convergence of the algorithm is only superlinear (not quadratic), because of
which more iterations are needed for the algorithm to converge. However, the cost per
step of Broyden's method is substantially less.
2.4.2.1. Powell's Hybrid Method
The major shortcoming of the QuasiNewton method is the necessity of a good
initial guess. To improve this property Powell's Hybrid method (Powell, 1970) is
employed in the SNSQ routine. The hybrid step, also called the "dogleg" step, is a
combination of quasiNewton and steepest gradient step.
In the Powell's Hybrid method, Equation 2.4.23 is modified as:
(2.4.2.11 )
where, I is a unit matrix and Ii' is a nonnegative p~ameter. As Il is increased the
• I
algorithm approaches the steepest descent method with small steps:
(2.4.2.12)
while as J1.k is decreased to zero the algorithm becomes GaussNewton.
28
t The value of pi. is based on an objective function, also called the minimization
function. The objective function used in the solver is the sum of squares of residual : A
residual may be defined as the difference between the function evaluated at the calculated
solution and the function evaluated at the actual solution. From the Equation 2.4.21, we
can observe that the function evaluated at the actual solution is zero. The objective
function can be expressed as follows:
n
F(x) = L [f (x)Y
i=,l
(2.4.2.13)
The value of j./' is calculated in such a way that the following inequality is satisfied.
F(Xk.+1
) < F(xk
) (2.4.2.14)
The algorithm begins with f.l set to some small value. If a step does not yield a
smaller value for F(x), then the step is'repeated with j.l' multiplied by some factor v > 1.
Eventually F(x) should decrease, since the algorithm takes small steps in the direction of
steepest descent. If a step produces a smaller value of F(x) , then Jlk is divided by v for
the next step. This enables the algorithm to approach GaussNewton, thus providing
faster convergence. The method thus provides a good compromise between the speed of
Newton's method and the guaranteed convergence of steepest descent.
Figure 2.4.2.11 shows the Powell's Hybrid step. The figure presents an example
performance surface, between two parameters, on which the method is applied. The black
dot represents the global minimum of the surface, which is (10, 3). The gray dot
represents the initial guess, which is (3, 3). The dashed arrow represents the direction
taken for small,l/ , which corresponds to the GaussNewton direction. The solid arrow
29
represents the direction taken for large p!' , which corresponds to the steepest descent
direction. The gray curve represents the hybrid step for all intermediate values of I1 Ie
•
15
10
x
2 5
o
o 5 10 15
• r
Figure 2.4.2.11. Powell's hybrid step
The HVACSIM+ equation solver implements Powell's hybrid method by using a
specialized algorithm called the DOGLEG algorithm that utilizes a combination of
steepestdescent and Newton steps in the process of minimizing a function. As long as
the gradient steps are relatively large, they are used. However, since gradient steps tend
to perform poorly in valleys on the performance surface, Newton steps are also used.
Newton steps, however, perform poorly with data points far away from the minimum.
Hence, the algorithm uses a bound on the minimum step size and provides a compromise
step, called a "dogleg" step, that combines the gradient and Newton steps. Depending on
the performance of the DOGLEG algorithm the bound is either increased or decreased for
future iterations.
30
The DOGLEG algorithm requires the following inputs:
• Function f (x)
• Starting values Xo
d I
• An initial radius (DELTA) to provide an upper bound on step size. The radius is
an estimate of a region about the current solution in which the algorithm can
predict the behavior of the function reasonably well. A default value of 1.0 is used
as an initial estimate.
• The maximum number of iterations
• Convergence tolerance for the gradient
• Relative coefficient change
The algorithm returns the step that it decided to take (Newton, steepestgradient or
dogleg). The acceptability of the Newton step is based on Euclidean nonns. The
Euclidean norm is defined below.
The Euclidean nonn, also called leastsquares nonn, of a vector
(2.4.2.15)
The dogleg algorithm is shown in the Figure 2.4.2.12. At the beginning of each
iteration, the GaussNewton direction is calculated. The stepbound limitation is then
checked. If the Newton's step is not acceptable, a gradient step is employed. If the
gradient step also fails to improve the minimization function, shown in the Equation
2.4.2.14, a dogleg step is attempted. If all the three steps fail, the algorithm repeats with
31
the same solution. Thebound radius (DELTA) is then increased or decreased and a new
iteration begins.
Calculate the Newton' s Step. t5N
No
Calculate the gradient Step. t5SG
If
l\t5sal L DELTA
Calculate the dogleg Step. t3D
Figure 2.4.2.12. DOGLEG Algorithm
The initial choice of DELTA can profoundly affect the perfonnance of the
algorithm  different values sometimes lead to finding different local minimizers. Too
small or too large a value can cause the algorithm to spend several function evaluations in
the first iteration to decide whether to decrease or increase DELTA.
32
2.4.3. Convergence Criteria in the Equation Solver
The convergence test implemented in the solver selects xK as a solution vector if
the following condition is satisfied: . .
(2.4.31)
where, e, is the error tolerance usually specified by the user, f (x) is a vector of residual
functions, and the double bars denote the Euclidean norms. The value of the error
tolerance depends upon the simulation setup and its initial values. As a rule of thumb. its
value may be greater than or equal to the sum of the absolute and relative error
tolerances.
Convergence is particularly a problem at the beginning of a simulation. when the
equation solver must work from the set of usersupplied initial conditions. If the initial
conditions are a large distance from the final solution, the equation solver may fail to
converge. The choice of the error tolerance for the equation solver is an important factor.
Even if the solutions are convergent, an improperly selected error tolerance can result in
unacceptably large errors.
The problem of slow and uncertain convergence to a selfconsistent set of initial
conditions is unavoidable, but the consequences can be minimized by the use of the
initialization option. This permits a series of simulations to be run from the same initial
state, without repeating an inaccurate and computationally expensive startup transient at
the beginning of each run.
The convergence test in the solver can be rewritten as Ilxk+1  xk II ~ e,lIxlc II '
33
where, Xi is the solution vector and e, is the error tolerance. This test is based on
Euclidean nonns. Another kind of nonn, called the Infinity norm, is used in the field of
optimization. The Infmity nonn, represented as II. t 'may be defined on a vector
Ilvll =maxlv I DO ISiS" I
If, for example, X k+1  Xi = [r, ...,rY ,where, r is any scalar quantity, then
,
Ilxk+1  Xl II = ~Irl and Ilxk+1  Xi t = Irl '
,
(2.4.2.12)
where, n is the number of variables. If n is large, the Euclidean norm of the gradient can
be large even if r is small. This can distort the convergence test and so it is preferable to
use the Infinity norm when large problems are solved. The use of norms in the
convergence test can have other side effects. Suppose that:
X
l

1 =[1.55564 0.00082 0.0000078Y
Xl = [1.55552 0.00011 O.OOOOOIIY
If we had chosen e, =103 then
!IX
k

1
 xk II.., =11[0.00012 0.00071 0.0000067Yll =0.0007201 S; 1.5555.103
,
and xk would pass this test. If, for an application, it were important that the significant
digits of all variables be accurate, then x k would not be satisfactory since its second and
third components are still changing in their first significant digit. The use of norms
emphasizes the larger components in a vector and so the smaller components may have
poor relative accuracy. This effect can be rectified by nondimensionalizing the solution
vector. This is achieved by scaling the solution vector, which can ensure that all the
34
variables in the transformed problem would have approximately the arne relative
accuracy.
A scale vector. represented here as d. where d = (dll d2 .d3
••• .dn
) , can
profoundly affect the performance of the optimization algorithms. Vector, d , should be
such that the elements Idjxi I, 1 ~ i ~ n • are iii comparable units. A reasonable choice of d
is obtained by guessing the upper bounds, represented as qj' on Ixj·', where x· is the
,solution vector, and setting dj equal to 11 qj' Several methods have been suggested to
estimate the scale vector. Gay (1990) summarized the available methods and concluded
that selection of a scaling method is problem specific.
Scaling in the equation solver is computed by estimating the upper bound of the
solution vector. q , and computing the scale as 1/q .There is no user interaction
required to perform scaling.
In order to solve the system of equations we need a starting value and, since f(x)
is not scale invariant, we also need a way to measure the scale of the residuals
(r(x) = y f(x)) at the beginning of the computations. We also may need to remeasure
the scale as the computation continues.
Since the initial conditions supplied by the user are usually not good enough, the
convergence is particularly a problem at the beginning of the simulation. If the initial
conditions are too far from the actual solution, the equation solver may not converge. In
some cases, errors made by the solver can cause erratic solutions for several time steps,
even if the solutions are convergent. The choice of the error tolerance for the equation
solver is an important factor that influences this behavior.
35
Whenever the solution does not converge, MODSIM.shows. a warning message
"Iteration not making good progress". This message is shown when the steps taken in the
gradient direction, while searching for a global minimum in the domain, exceed 10. The
number of iterations, which limit the gradient step bounds, are used by the algorithm to
check the progress of simulation. For certain simulations, a closer investigation of the
residuals, printed in the diagnostic report,. revealed that the solution failed to converge
even though it had satisfied all the stringent convergence criteria. In addition, the
program did not show any warning messages. I.'
The convergence criteria used in MODSlM do not directly consider the
magnitude of the residuals. In a system that simulates mainly the temperatures (which are
mostly in the range of 10 to 60 degrees centigrade), a residual greater than 1 degree
affects the results significantly. Moreoyer, in computing the heat transfer rates, where the
temperature differential is multiplied by high flow rates, the effect of high residuals is
very significant. The only limiting criterion of the solver, to check whether the optimum
solution is reached or not, is the step size used by the solution algorithm.
The convergence p"roperties of the equation solver are also dependent on the
Block/Superblock structure. The user must carefully defme blocks and superblocks to
obtain good convergence since they determine the equations sets to be solved during a
simulation. This fact limits the flexibility of simulation construction.
2.4.4. Integration ofStiffOrdinary Differential Equations
The algorithm, implemented in MODS1M, to integrate differential equations is
based on a method developed by Brayton, Gustavson and Hachtel (1972). This method
uses variable time steps and variableorderintegration mchniques to solve sets of
36
differential equations and hence can significantly reduce the amount of computer tim.e
required for dynamic simulations. This section highlights the important aspects of this
method.
A general system. of differential equations may contain higher order derivatives,
but such a system may be reduced to firstorder fonn by introducing new variables. A
system of firstorder differential algebraic equations can be expressed as:
f(x,x,t) =0 (2.4.41 )
where, x is a state vector that is a function of time, t, and x is the derivatives of x.
If the solution vector x(t) of the above equation had been obtained at previous discrete
times, t = tn' t = (nl , ... , and t = tn+J_ p ' then the solution xn+1 at the current time, t =t"+I'
satisfies:
(2.4.42)
For stiff differential equations, having widely varying time constants, the pth
order backward differentiation formula is given as:
(2.4.43)
where, aj are constants and h is the step size (In+!  tn ) •
Substituting equation 2.4.53 in 2.4.52 yields a set of nonlinear algebraic
equations in X"+l at time ("+1' The obtained system of equations is solved by the equation
solver. If the solution is converged, the iteration is tenninated and the order, p, and step
size, h, are updated as follows.
The truncation error of pth order backward differential formula is given as,
37
r (2.4.44)
where, x:+1 is the guess value of x,,+} . The guess value is formulated using the same
regressor expression in Equation 2.4A3. The guess value is expres ed as,
p+l
x g ,,+1 = "L.'Jr;xn+1;
;=0
where, ri are constants.
If the truncation error satisfies the condition,
(2.4.45)
(2.4.46)
where, tif is the time interval, er is the relative error tolerance and eo is the absolute
error tolerances, the step size and the order are acceptable. The time interval and the error
tolerances are specified by the user.
To update the step size and the order, the above condition is applied for orders
p +1, P and p 1 . The order corresponding to the largest time step size is seleCted. The
computation of (p 1) th order differential is skipped if p is I, in which case, only the
step sizes at p and p + 1 are compared. The maximum order of p is limited to six.
2.4.5. Interpolation ofdata
In the Section 2.3.3, it was mentioned that HVACSIM+ interpolates the data in
the boundary value file to find the boundary values at simulation times. This section
describes the application of Lagrangian interpolation techniques to perform interpolation
of data.
38
During a simulation, the MODSIM program might take time step values for
which boundary values are not specified in the boundary value fIle. In order to fmd the
data values for the boundary variables the third order Lagrangian interpolation method is
used. This method uses the simulation time to interpolate the data. It is expressed as
follows.
x(t) =L4[IT4(~ttJlXi
;=0 j=\ t; t j
j~j
where, x(t) is the interpolated boundary variable at time t, and Xi =x(tj ) •
39
(2.4.51)
CHAPTER 3
HVACSIM:+ ENHANCEMENTS
Although HVACSIM+ was developed over two decades ~go, a number of
problems exist that make the program difficult to use. These problems, as documented by
Clark (1985) include:
• complexity of the user interface
• difficulty in debugging the simulation
• improper convergence of the equation solver
3.1. Development ofthe userinterface
The user interface of HVACSIM+ may be considered primitive. The program
interacts with the user through a combination of command and menudriven systems. The
menu driven program control has a hierarchy of menus that the user can select by typing
one or more characters to direct the program. However, if the user makes an error while
configuring a system, there is no way that it can be corrected immediately. The user has
to either abort the process or supply arbitrary values to get to the next command prompt.
Aborting the process discards all changes that have been made. Once the user returns to
the main menu, the EDIT command can be chosen to edit the changes in the system
configuration. This is usually time consuming and requires patience.
Each input and output variable in a simulation is given an index number. The user
uses this number to identify the variable that needs to be processed for setting
connections or setting initial values. Keeping track of the indices is required in order to
40
represent the system correctly in a physical sense. Assigning numbers to variable in
large systems is tedious and time consuming. Repetition of indices, or 90nnectivity of one
category with a different category (e.g. control signal with pressure) distorts the whole
system configuration. The simulation work file, created using HVACGEN, and the model
definition file, generated by SLIMCON, are also difficult to understand.
Using an eventdriven approach instead of the questionandanswerbased logic
can rectify the problems with the userinterface. In the questionandanswer method, a
central controlling program detennines the sequence of operations for an application.
Eventdriven programming determines the sequence of operations for an application by
the user's interaction with the application interface (forms, menus, buttons, etc.). For
example, rather than having a main procedure that executes a series of processes, an
eventdriven application remains in the backgrou.nd until the user picks the processes to
be performed.
An event is an action recognized by an element of a graphical user interface.
Some examples include a 'mouse click' on a button or a 'key press' on the keyboard.
Events can also relate to actions that trigger processes. For example, a 'simulation ended'
event at the end of a simulation could, among other things, close any open files.
Eventdriven progranuning and graphica user interfaces (GUIs) are related since
fonns and the graphical interface objects on the fonns serve as the skeleton for the entire
application. In an eventdriven application, small programs are attached to events
associated with objects. In this way, the behavior of the application is determined by the
interaction of a number of small manageable programs rather than one large program.
The primary advantages of eventdriven programs are the following:
41
• Flexibility  Events rather than a equential program control the flow of the
application. Thus, the user does not have to conform to the programmer's
understanding ofhow tasks should be executed.
• Robustness  Eventdriven applications tend to be more robust since they are
less sensitive to the order in which users perform activities. In conventional
programming, the programmer has to anticipate virtually every sequence of
activities the user might perform and defme responses to these sequences.
As discussed in Chapter 2, the main programs in HVACSIM+ that provide an
interface for building and executing a simulation are the HVACGEN and MODSIM
programs. Both these programs use the conventional questionandanswer based
approach discussed above to accept input from the user. HVACGEN is used to create an
initial work file that is processed by the SLIMCON program to create a definition file.
The MODSIM program reads the definition and the boundary value files and simulates
the system. The user should use all three programs to run a simulation.
The user interface for HVACSlM+ has been improved by creating a Visual
Modeling Tool to replace HVACGEN and adding residual reporting feature in MODSIM.
3.1.1. Visual Modeling Tool
The visual modeling tool has been developed primarily for research purpose.s and
is intended to carry out dynamic simulation studies,. A detailed description of the tool is
given in Chapter 4. The features of the tool can be summarized as follows:
• The tool uses icons and pictures to represent units, boundary variables and
connections between them.
42
• The tool imitates the functionality of HVACGEN u ing an eventdriven
approach but also goes a step forward and makes calls to SLIMCON and
MODSIM programs in the background.
• The tool eliminates the need for the user to use three different program to
build and run a simulation.
• The user is not required to keep track of variable indices.
• The tool provides facilities to plot the simulation output.
• An online help system supports the tool.
3.1.2. Residual Reporting
During a simulation, HVACSIM+ does not use any validation tests to determine if
the results being generated are reasonable. For simulations that run several hours such
tests would allow the users to decide whether to continue with the simulation or not.
Since no such test is available, the user is forced to wait for the simulation to end, analyze
the results and rerun the simulation with any corrections.
The solution algorithm and its associated subroutines are more difficult to debug,
since the code lacks structure and logical flow. The GOTO statement has been used
extensively in the code, which makes the code difficult to comprehend. Variables are
named inconsistently. The source code also lacks necessary documentation.
The only reprieve during the debugging phase is the diagnostic report that can be
printed out for a specified duration in the simulation. The diagnostic report can be made
to print information such as the Jacobian, residual information, "hybrid step" size of the
optimization algorithm (indicating the progress of the iteration), limiting convergence
43
criteria and intermediate solution vector. This information, although useful, is in ·ufficient
to find the source of most problems.
In order to give a useful picture of the simulation progress, the interface is capable
of plotting residuals while the simulation is running. This has been achieved by adding a
child window in which the plotting is performed. The user can use the menu to switch
from the simulation status window to the plot window.
The residual plot feature can only be activated if MODSIM is called from the
visual tool. To activate the feature the user is required to interact with the form shown in
Figure 3.1.21. Complete details about this dialog are presented in the Chapters 4 and 5.
The user can decide not to show the plot only at this stage. Once the plot window is
activated, there is no way to close it.
The application plots the Euclidean nonn of the residual vector against time on a
bar graph. The scale on the yaxis, which shows the norms, can be linear or log. The scale
on the xaxis, which shows the time step values and iterations within them, i always
linear. Each time step is allotted a random color that remains the same for all iterations
belonging to that time step. A filled circle represents convergence of the algorithm over a
bar. Figure 3.1.22 shows a sample plot of residuals generated during a simulation. The
plot shows that during the initial stages of the simulation the residuals are high. This is
because of poor initialization of variables. After a good guess is made, the algorithm
performs better. After the first two time steps, the residuals in the plot are zero. This is
probably because of no loads in the simulation or because the algorithm is performing
very efficiently. The peaks start appearing when there is considerable effort required by
the equation solver to find a solution. As the iterations for a time step proceed. the
44
residuals should reduce. Some time steps indicate convergence although the residual
values are not zero.
Okabl. F'"""I!VeriebIe ...... 7
 .......rn Tim. Slop (second.
s__
OOl"'U__ID .... elQet><afC_1D be .........d ....... ol ...._ .... ?
Figure 3.1.21. Simulation Requirements dialog used to initiate residual reporting
Figure 3.1.22. Sample residual reporting screen (LOG scale)
45
The program also generates an output file that contains the plotting data,
including the xaxis and yaxis values from the residuals plot discussed above. This file
has the file extension' . DIG'. It will be referred to as the Diagnostic Information file.
Figure 3.1.23 shows a sample of the Diagnostic Information file.
***********************************************************•••****.**•••••*****
DIAGNOSTIC INFORMATION (RESIDU PLOTTING) FILE for
WeatherFord Bridge Deck With Temperature Difference Control for Recharge
.*.***.*.****••*.*.*****••*********••••*••••••••••••***•••••••••••••••••**.*••
TIME IlRESIDUAL VECTORII CONvERGED IN (ITERATIONS)
••**•••••**.*.*•••••*••••••**••*••••••••••*••••**•••*.*.**.**••••*••••**••***.*
3600.000
3600.000
3600.000
7200.000
7200.000
7200.000
7200.000
7200.000
10800.00
10800.00
10800.00
10800.00
10800.00
24.42056
24.42056
24.42056
17.2526~
17.25268
17.25268
17.25268
17.25268
356.6238
106.7234
63.81583
59.72607
59.28603
oo3o
ooo5
ooo
o
0.
Figure 3.1.23. Sample Diagnostic infonnation file
The visual tool c0nsists of a routine that can read the Diagnostic Infonnation file
and regenerate the residuals plot. Figure 3.1.24 shows a sample residual plot generated
by visual tool. The residual plot is a barplot of Euclidean nonn of the residual vector
versus time. The plot features log scale on the yaxis and a linear scale on the xaxis.
Each bar represents an iteration for a time step, between two tick marks, shown on the xaxis.
46
I Semple Plot of "ldUII." I
"t':!::"",.: :~.' ~
~_.
...... " .....•.. _._..
1.CIlOOOO
0.10aa00
:0
I
iJ·Ol0c00
J!
i
'80.001000 Ic
r·oooloo
~w
D.000010 !'aln:.
Figure 3.1.24. Residual plot generated by the visual modeling tool
3.2. Development ofnew convergence metrics
Since the initial conditions supplied by the user are usually poor, the convergence
is often a problem particularly at the beginning of the simulation. If the initial conditions
are too far from the actual solution, the equation solver may not converge. In some case,
errors made by the solver can cause erratic solutions for several time steps, even if the
solutions are convergent. The choice of the error tolerance for the equation solver is an
important factor that influences this behavior.
Whenever the solution does not converge, MODSIM shows the warning message
"Iteration not making good progress". This message is shown when the steps taken in the
gradient direction, while searching for a global minimum in the domain, exceed 10. The
number of iterations limits the gradient step bounds. For certain simulations, a closer
investigation of the residuals, printed in the diagnostic report, revealed that the solution
failed to converge even though it had satisfied all the stringent convergence criteria and
did not show any warning messages.
47
In order to ensure that HVACSIM+ would not consider a simulation time tep
converged when unacceptably large residuals were present, it was necessary to establish
new convergence metrics. This section discusses the development, implementation and
testing of new metrics in the equation solver. The new metrics aim to provide good
convergence criteria and increase the simulation speed.
3.2.1. Design ofnew convergence metrics
The convergence criteria used in MODSIM do not directly consider the absolute
value of the residuals. In a system that simulates mainly temperatures (which are mostly
in the range of 10 to 60 degrees centigrade), a residual greater than 1 degree affects the
results significantly. Moreover, in computing the heat flux, where the temperature
differential is multiplied by high flow rates, the effect of high residuals is very
significant. The only limiting criterion of the solver is the step size.
The major shortcoming of the existing HVACSIM+ convergence criteria is that it
does not consider the absolute value of the residuals. The effect of not taking m~gnitudes
was explained in the Section 3.1.3. This fact motivated the development of new
convergence criteria.
The convergence criteria, described in the Section 2.4.4 of Chapter 2, are
represented as:
(3.2.11)
where, x is the solution vector, f(x) is the residual vector, and er is the error tolerance.
For complex systems, the residual functions will rarely be equal to a numerical
value of zero. It is unlikely that the calculated value of the functions would ever be
48
exactly zero because of rounding errors in computer calculation . Even if there were no
rounding errors, no algorithm is guaranteed to find such a solution in a finite amount of
time.
As an alternative, a tolerance value can be used for the residual functions. This
modifies f(x) in the Equation 3.3.11 as:
(3.2.12)
where, ereridual is a small value. However, identifying an approximately small value for
the tolerance may not be easy. Equation 3.3.12 faces another problem when the
magnitude of the solution vector can vary over a wide range. Suppose that the
measurement units were changed. For example, suppose that instead of measuring the
objective in terms of kilometers it was measured in tenns of millimeters. This would
cause the function to be multiplied by 106 and hence would cause the norm to be
multiplied by 106
. Simply changing the units would make the convergence test much
more difficult to satisfy unless the error tolerance was also changed to reflect the change
in units.
During the simulation, the solution vector holds values for different kinds of
variables irrespective of their units. For example, a solution vector might contain
variables that are temperature and pressure variables. Having a single tolerance value that
is applied to all variables may lead to inaccurate results. Using separate tolerance values
for different categories of variables could solve this problem by allowing the user to
assign different convergence criteria to different types of variables. For example, mass
flow rates in a system may be less important than temperatures. A large tolerance may be
applied to variables that represent massflowrates and smaller tolerance to temperature
49
variables. This way the algorithm would spend mQre time on important variable . Let
e,ol_vu be an array that contains the tolerance values to be applied on the olution vector.
Each element in erol _vee is based on the category of the variable in the solution vector. Let
eresid/iQl_wc be an array that contains the tolerance values to be applied to the residual
vector. Each element in eresidU/ll_vu is based on the category of the variable in the residual
vector. The new convergence criteria may be represented as:
(3.2.13)
From the Equation 3.3.13, it can be observed that the residual vector is no longer
normalized. The absolute value of each scalar in the residual vector is compared with the
corresponding tolerance value in the eresiduaJ _we vector. In addition, if the vector, erol_vec'
contains the same values arid the vector, eresidual _"lC ' contains zeros then the Equation
3.3.13 represents the Euclidean nonn, given in the Equation 3.3.11.
3.2.2. Performance evaluation ofthe new convergence metrics
In order to test the new convergence metcics two systems were selected: the
Geothermal Smart Bridge system and the Single zone air handler system with a cooling
coil. These systems were simulated using the new convergence criteria and the output
generated was analyzed.
The Geothermal Smart Bridge system was developed as a result of ongoing
research at Oklahoma State University. The simulation models a bridge deck heating
system that will eliminate preferential icing on bridges and reduce bridge maintenance
costs. This system makes use of ground source heat pump technology and hydronic
circuits embedded in the concrete bridge deck. The system is comprised of two parts,
50
namely, the loadside and the sourceside. The load ide ubsystem consists of the
bridge deck. The sourceside consists of a ground loop heat exchange.r (GLHE) usually in
the fonn of a vertical borehole field. Heat transfer from the earth to the fluid takes place
in the ground loop heat exchanger. 1\ gang (bank) of watertowater heat pump connects
the loadside and sourceside subsystems. A controller is used to determine the number
of heat pumps to be operated in th.e system. The controller also sends a designed value of
mass flow rates of the fluid to the loadside and sourceside subsystems. In addition, the
GLHE model supplies the average temperature of the fluid in the exchanger to the
controller.
During the development of Smart Bridge technology, several configurations were
simulated. One such configuration, called "Design Procedure Step 2", will be used as a
test case. In this procedure, only the loadside subsystem is simulated with constant
sourceside conditions. A constant value of  3°C is supplied to the heat pump as the
sourceside entering fluid temperature. The average fluid temperature in the GLHE
model, supplied to the controller, is set at  l()(}OC. The Figure 3.2.21 shows the system
schematic for this configuration.
GANGOF
WAffiR.TQ..WATER
CONTROLLER HEAT RJMPS
WEATHER
[)\TA
r,
I I
II o BRIDGE DECK
Figure 3.2.21. Geothennal Smart Bridge (Design Procedure Step 2)
51
The components and their functions are listed below:
• TYPE 700:' Slab model  This model represents the bridge deck, the
temperature of which is maintained at a specific value.
• TYPE 711: Gang ofWater to Water Heat Pumps model  This model
simulates a series of heat pumps. Th,e number of heat pumps in the series is
detennined by the TYPE 730, the linear proportional controller.
• typE 750: Simple Pump model This type is similar to TYPE 1 Pump that is
distributed with HVACSIM+.
• TYPE 730: Linear Proportional Controller for heating  This model
determines the number of heat pump pairs that should be operated based on
the conditions of the bridge deck.
The design consists of two superblocks. The flow controller that decides the
number of heat pumps is the only unit in Superblock #2. All the other components are
grouped together in Superblock #1. A design data set, consisting of one snow event, was
created for the simulation to be used as boundary values. The equation sets in both
Superblocks consist of only temperature variables. The average bridgesurface
temperature was used for analyzing the results.
The single zone air handler with a cooling coil is discussed in Section 5.1.5 of
Chapter 5. The system consists of multiple blocks. Block 1 represents a cooling coil
control loop, block 2 represents an airhandling unit, and block 3 represents a room
configuration with a supply zone, a conditioned zone and a return duct.
52
The new convergence criteria were applied to the selected te t cases. The overall
execution time and residuals were inspected during the simulations. The test results and
inferences made from them are presented below.
Since the new convergence oriteria are conditioned by the error tolerances fa the
solution and the residual vectors, only they were changed in the tests. The primary aim
was to show that the results generated are acceptable when different limits are applied to
the residuals. The improvement in the simulation speed is automatically achieved when
the convergence criteria are not rigid. Comparative analysis was performed using the
results obtained with the original criteria.
All the simulations were run on a computer using an AMDAthlon Processor
running at 750 MHz and having 128 MB RAM. The simulation execution times were
estimated by making two calls to the F:ortran function, SECNDS. The fIrst call sets the
timer to zero seconds and the second call estimates the time taken since the first call. The
code to initialize, simulate and report results is placed between the two call statements.
The visual modeling tool was used to plot the residuals generated during the
simulations. However, the plots have been modifIed so that each time step is drawn with
a different color shade making the plots more readable.
3.2.3. Simulation results for the Smart Bridge system
In order to conduct a comparison analysis, the Smart Bridge model was fIrst
simulated using the original convergence criteria. The error tolerance applied to the
solution vector was selected as 0.0002. The simulation time step was fixed at 3600
seconds (l hour). A simulation stopping time of 338400 seconds (approximately 4 days)
was selected.
53
The simulation execution time was 45.2 econds. The re iduals in the simulation
are zero until the snow event begins at 169200 econds, as hown in Figure 3.2.31. Thi
is because the heat pumps are not running when there is no snow or ice formation on the
bridge and hence, there are no equations to solve in the system. The residuals are high
during the first few iterations of every time step taken during the 'snow event. As the
iterations for a time step proceed, the equation solver is getting closer to finding the
solution for that time step. This can be observed by the reduction of the residuals. The
solution does not converge when the time is 216000 seconds, as shown in Figure 3.2.32.
The solutions converge for all other time steps.
Figures 3.2.31 and 3.2.32 show barplots of Euclidean nonn of the residual
vector versus time when the snow event begins and during the time step when the
algorithm does not converge, respectively. Th.e plots feature log scale on the yaxis and a
linear scale on the xaxis. Each bar represents an iteration for a time step, between two
tick marks, shown on the xaxis. Time steps can be distinguished by the color of the bars.
Smart Bridge Syllem Cllign Procedure Slep 2
, 00000o
ii
11
~ ounm
!
I ODlOOllO
~
Ej 0001000
Jo~oo
0JlOOO10
00000D1
~::.~:.;.:~::::::.:.:: ~::::::.::::
. _._ ..

,,..00
Figure 3.2.31. Plot of residuals when the snow event begins
54
I Smart Bridge Sy~em Oestgn Procedure Step 2 I
Figure 3.2.32. Plot showing the time step (216000 seconds) when the solution did not
converge
The Smart Bridge model was then simulated using various test conditions. Three
test environments were selected, which include:
• varying only the error tolerance applied to the solution vector
• varying only the error tolerance applied to the residual vector
• varying both the error tolerance for the equation solver and the error tolerance
for the residual vector
Since the equation sets defined for the model include only temperature variables,
the error tolerances were specified only for temperatures. The error tolerances for all
other categories of variables were set at 0.0002 and at 0.0 for the equation solver and for
the residual vector, respectively.
The Table 3.2.31 presents the results of simulations obtained by varying only the
error tolerance applied to the solution vector. The error tolerance for the residual vector
was maintained at zero for all categories of variables.
55
Test #
Error tolerance applied Simulation time
to the solution vector (seconds)
סס. 1 0 OO1 58
2 , 0.0001 48
3 1, 0.002 39
4 T .~i 0.02 36
5 " .. ., 0.2 34 . ,
Table 3.2.31. Results of simulation performed while varying only the error tolerance
applied to the solution vector
Convergence was not achieved when the error tolerance was set to 0. סס OO1. As
the error tolerance was reduced to 0.002, the simulation output matched the output
generated by using the original criteria with the error tolerance for the solution vector was
set at 0.0002. When the error tolerance was further reduced, the simulation run time was
also reduced since convergence was achieved within the fIrst three iterations for all time
steps. Figure 3.2.33 shows three residual plots for different error tolerances for the time
steps, 212400 and 216000 seconds. The error tolerance values for plots (A), (B), and (C)
are 0.00001, 0.002, and 0.2, respectively. Figure 3.2.34 shows the plot of the average
bridgesurface temperature against time for the original criteria (error tolerance of
0.0002) and the new criteria with the error tolerances value set at 0.2. As shown in the
figure, a less restrictive error tolerance results in a larger error band around the deck
temperature. The fluctuations begin at the time 280800 seconds and end at about 334800
seconds. At all other times, the simulation converges to the same solution for both
criteria.
56
,axao
oamiO
Ol/JOJ. fUJL......,
(A) Error tolerance = 0. סס OO 1 (B) Error tolerance =0.002
(C) Error tolerance =0.2
Figure 3.2.33. Plot ofresidual for different error tolerance for the solution vector
4 ._.
3lleOO 75eOO 111800 147llOO
2 t
~  
 
21ueoo 26l58OO 2lI1tlOO 327110O
.10 l. '
Error toIerllnC8 =0.0002 (OrlglnaJ Crtlerla) Error tole<ance = 0.2 I
Figure 3.2.34. Plot ofthe average bridgesurface temperature against time when the
error tolerance for the solution vector is varying
57
The Table 3.2.32 presents the results of the imulation obtained by varying the
error tolerance for the residual vector. The error tolerance for the solution vector wa
maintained at 0.0002 for all categories of variables.
Test #
Error tolerance applied Simulation time
to the residual vector (seconds)
סס. 1 0 OO1 45
2 ' <, 0.0001 41
3 0.001 40
4 0.01 36
Table 3.2.32. Results ofsimulation performed while varying the error tolerance for the
residual vector
From the Table 3.2.32, it can be observed that the simulation with error tolerance
set at 0.00001 required the same time as the simulation with a tolerance of 0.0 (the
original criteria).. This suggests that even a tolerance value of 0. סס OO1 is rigid to achieve
convergence. Convergence was not achieved when the error tolerance for the residual
vector was set at 0. סס OO1. The simulation output showed larger surface temperature
fluctuations when the error tolerance was set at 0.01, suggesting that the value selected
was high. Figure 3.2.35 shows two residual plots for different error tolerances at the time
steps taken at 212400 and 216000 seconds. The error tolerance values for the plot (A),
and (B) are 0. סס OO1, and 0.01, respectively. Figure 3.2.36 shows the plot of the average
bridgesurface temperature against time for the error tolerances, 0. סס OO1 and 0.01.
58
Smart Bridge Syatem Design Procedtro SI p 2 I
(A) Error tolerance = 0. סס OO1
i ..... :
jI

'Ii"",..
jj
] ...... w
(B) Error tolerance =0.01
Figure 3.2.35. Plot ofresidual for different error tolerance for the residual vector
4r,
10 ~ ____J
IEmrtdennle =0.0lXXl1  Error toIenn:e ,,0.011
Figure 3.2.36. Plot ofaverage bridgesurface temperature against time when error
tolerance for the residual vector is varying
It can be observed from the Figures 3.2.34 and 3.2.36 that setting a tolerance
value of 0.2 for the solution vector or a tolerance value of 0.01 for the residual vector
produces similar results. Using unreasonably high values for the error tolerances makes
59
the convergence criteria loose, which cau e the solutions to con erge during the first few
iterations even before the algorithm has approached closer to the correct solution.
The Table 3.2.33 presents the results of the simulations obtained by varying both
the error tolerances applied to the residual vector and the solution vector.
Test #
Error tolerance applied Error tolerance applied Simulation time
to the solution vector to the residual vector (seconds)
1 0.002 0.0001 42
2 0.002 0.001 40
Table 3.2.33. Results ofsimulation performed while varying the error tolerances
The output generated through the tests shown in the above table matched the
output generated by using the original criteria having very rigid error tolerances. The
simulation running time was reduced. Similar to the previous test conditions, the
algorithm did not converge for the time step at 216000 seconds.
The saving in the simulation run time is significant for annual simulations. Using
the original criteria, the simulation would take 4210 seconds (1.17 hours, approximately)
to end. Using the conditions for Test #2 shown in the Table 3.2.33, the simulation would
need only 3728 seconds (1.04 hours, approximately).
3.2.4. Simulation ofSingle zone air handler system with a cooling coil
In order to conduct a comparison analysis, the system was first simulated using
the original convergence criteria. The error tolerance for the solution vector was selected
as 0.0005. Variable time step is used in the simulation with a minimum time step of 0.1
seconds and a maximum time step of 300 seconds. A simulation stopping time of 86400
seconds was selected.
60
The simulation execution time was recorded as 37.2 seconds for 7554 overall
iterations. During many time steps, the equation solver demonstrated its characteri tic of
reporting convergence although large residuals were present in the solution . The Table
3.2.41 shows a list of selected time steps, the Euclidean nonn of the residual vector
when the solver reported convergence and the number of iterations taken for that time
step. The data in the table was extracted from the Diagnostic Information file generated
during the simulation. The data in the table illustrates the convergence problem
mentioned above.
Time (seconds) Norm of Residual vector Number of iterations
7858.3 0.11028 19
54940.8 0.00759 270
56140.8 0.00376 24
71740.8 , 0.007326 29
80440.8 0.031 37
Table 3.2.41. Performance ofthe simulation at selected time steps
The solution vector consists of only pressure, mass flow rate and temperature
variables. The system was simulated using various combinations of error tolerances
applied to these variable categories. The first set of tests involves varying the error
tolerance applied to the solution vector for pressures, mass flow rates and temperatures.
The second set of tests involves varying the error tolerance applied to the residual vector
for these categories.
In the first set of tests, the error tolerance applied to the residual vector was
maintained at zero and the error tolerance applied to the solution vector was changed. A
tolerance of 0.005 was applied to each category in the solution vector and the
performance of the solver was analyzed while maintaining the tolerance for other
61
categories at 0.0005. A summary of the test results is shown in Table 3.2.42. The table
displays the error tolerance value, category on which the test is applied and imulation
run time for the test. The table also displays the total number of iterations performed in
the simulation.
Test Error tolerance applied Categories Simulation Number of
# to the solution vector selected time (seconds) iterations
1 0.005 Pressure 27.8 5453
2 0.005 Flow Rate 30 6052
3 0.005 Temperature 25.1 4934
4 0.05
Pressure and
27.4 5460
Flow rate
5 0.05
Flow rate and
24 4579
Temperature
Table 3.2.42. Results ofsimulation performed while varying only the error tolerance
applied to the solution vector
From the above table, it can be observed that by reducing the error tolerance on
pressures and temperatures we achieve faster convergence compared to reducing
tolerance on flow rates. The solution vector contains more pressures than any other
categories of variables. Thus, by reducing the tolerance on the pressure variables, the
solver required less effort to perform computations. In addition, the presence of only one
temperature in the solution vector influenced faster convergence in the simulation when
the tolerance value for temperatures was increased.
The simulation output matched the original output for the fIrst three test cases.
Test #4 and Test #5 generated results that showed fluctuations in output at irregular time
intervals.
In the second set of tests, the error tolerance applied to the solution vector was
maintained at 0.0005 for all categories of variables and the error tolerance applied to the
62
residual vector was changed. A summary of the te tresults is shown in Table 3.2.43.
The table displays the error tolerance value, category on which the test was applied and
simulation run time for that test. The table also displays the total number of iteration
performed in the simulation. The error tolerance values for categories not part of the tests
were maintained at 0.01.
Test Error tolerance applied Categories Simulation Number of
# to the residual vector selected time (seconds) iterations
1 0.01 All categories. 24.34 4721
2 0.1 Pressure 17.26 3111
3 0.1 Row rate 24.73 . 4721
4 0.1 Temperature 24.33 4677
5 0.1
Pressure and
20.33 3729
Flow rate
6 0.1
Flow rate and
23.9 4677
Temperature
Table 3.2.43. Results ofsimulation performed while varying only the error tolerance for
the residual vector
Convergence was achieved in all the test cases shown in Table 3.2.43. However,
the solver reported convergence when the residuals were large. When the error tolerance
was set at 0.01, as in Test #1, the simulation generated output that matched the output
generated when the tolerance was set at 0.0. The simulation run time was the shortest
when the tolerance on the pressure variables was set at 0.1. When the tolerance of both
the pressure and flow rate variables was controlled, the simulation showed fastest
convergence with only 3729 overall iterations.
The two test sets show that the new convergence criteria cannot rectify the
convergence problem in the equation solver. However, they provide means to control the
simulation to yield similar results with shorter run times.
63
CHAPTER 4
TECHNICAL OVERVIEW OF THE USER INTERFACE
4.1. Concept
The Visual Modeling Tool, discussed in the Chapter 3, provides a powerful
graphical user interface for HVACSIM+. The interface provides:
• a workspace for representing infonnation about a HVAC system.
• a means to define the environment of such a system in terms of boundary
values.
• facilities to simulate the system.
• a means to plot the results.
4.2. Architecture
The core components of the application include graphical administration tools, the
boundary file editor, the system simulator, and the plotting tool.
The graphical administration tools create components in the workspace, provide
sufficient graphical infonnation on the component icons, perform various editing
operations like cut, copy and paste, and remove icons from the workspace. These tools
also handle operations like dragging components around the workspace, to provide a
visually wellorganized system configuration.
The boundary file editor tries to mimic a traditional spreadsheet organizer having
a matrix of rows and columns, where columns represent boundary conditions and rows
64
represent data values for all the boundary conditions at specific times. It provides
facilities for creating boundary value files, which supply input data to the simulation.
Additional support allows for associating columns of data to boundary conditions in th
simulation, and handling editing operations on the data between applications (like
Microsoft Excel) using the windows clipboard as an intermediate buffer.
The simulation manager analyzes the workspace for basic errors and ets up an
environment for calling backend applications like SLIMCON and MODSIM. The
simulation manager identifies errors from MODSIM and interprets these errors for the
user. It also modifies the output files to add additional information, like variable names
assigned by the user, to clearly differentiate between them in these files.
The major purpose of the plotting tool is to plot the data in the output files. The
user is allowed to select reported variables and the tool reads the variable data and the
associated time values from the output file.
. l
4.3. Data Format
All the information entered by the user is written into an input file with the default
file extension of ' . DAT'. The user can save his work any time while using the visual tool
using the 'Save' or 'Save As' options under the 'File' menu. This input file has a special
fLle fonnat, shown in the Figure 4.31, that include the following:
• Header: A header identifying that the file belongs to the application.
• Block and Superblock details: Information about the number of blocks and
superblocks in the workspace.
• Unit Information: Information about the various units involved including the
location of the icons in the workspace, superblock and block numbers, type
65
•
number, initial and parameter values, reporting status of all the input and
output ports, and variables names.
Boundary condition information: Infonnation about the timedependent and
constant boundary conditions involved including their location in the
workspace, name, and constant value.
• Connection information: Information about the connections involved
including the originating port identities and ending port identities.
• Simulation information: Simulation details including the title, error tolerances,
simulation times, and answers to questions related to diagnostics.
• Boundary value information: Information relating to the boundary file such as
the header names of columns and the data values in the grid.
After execution of the simulation, MODISM generates three output files,
discussed in the sections 2.3.4. The application appends the raw output file with
additional information required for plotting the results. The format of the additional
information is shown in the Figure 4.32. The other two files generated are not modified
by the application. The plotting routine can also read output files without having the
additional information. Each reported variable is identified with its Unit and TYPE
numbers. For example, an inlet temperature for a fan or pump may be represented as Unit
1, TYPE 1: Inlet Temperature.
MODSIM also generates the Diagnostic information file, mentioned in Section
3.2.2. This file is not modified by the visual tool and is used only to plot the residuals.
66
Header }
Block
&
Superblock
Details
Unit
Information
Boundary
Condition
Information
Connection
Information
Simulation
Information
Boundary
Value
Information
••*••••**••••••••••••••••••••••••••••••••••*
• VISUAL MODEUNG TOOL FOR HVACSIM+ •
• In~tffie • ••**•••*•••*••••••••••••••••••••*•••••*••••*
<Blank line>
<Number of Superblocks>
<Number of Blocks in each Superblock on separate lines>
<Active Superblock number>
<Active Block number>
<Number of units involved in the simulation>
[List of all units involved in the following fonnat for each unit]
<XpositiOn>, <Yposition>, <Superblock>.<Block>,<fype number>. <Initial value >, <Parameter
values>, <Reporting status; 0 or 1>
<Variables names, each on a separate line>
<Number of boundary conditions in the simulation>
[List of all boundary conditions involved in the following format for each condition)
<XPosition>, <YPosition>. <Is it a constant condition; 0 or 1>, <Category index>, <Constant
value>, <Name>
[For constant conditions. Name is •. ]
<Number of connections in the simulation>
[List of all connections involved in the following format for each connection]
<From unit>, <From port>, <From category>. <To unit>, <To port>, <To category>
.[Connections involving boundary conditions have 1 for the From port and From category values.)
<Simulation title>
<Relative error tolerance>
<Absolute error tolerance>
<Tolerance for simultaneous equations>
<Error tolerance for time interval>
<Reporting interval>
<Variable freezing option>
<Variable scanning option>
<Minimum time step value>
<Maximum time step value>
<Simulation stopping time>
<Answers. one in each line, to questions related to diagnostics>
<Total number of headers defined>
[List of header names in the foHowing format:1
<Header name>
<Total number of rows>
<Total number of columns>
[Boundary value data in the following format:)
<Value in each column; comma delimited>
Figure 4.31. Input file structure
[The raw output listing)
••••••••••••••••••••••••••••••••••••••••••••••••••
<Simulation title>
<Date and time of creation>
<Location of the corresponding input file>
<Number of Superblocks>
[For each Superblock.. the reported variable information is represented in the following fonnat:]
<Number of Reported variables in the Superblock>
<Variable name, one on each line>
Figure 4.32. Format ofadditional information in the raw output file
67
4.4. Configuration
Configuration of the product is carried out by identifying the corrict Types
component library ftle (TYPAR.DAT) and copying the file to the application directory
before starting the application. The Types component library contains information about
all the types which SLIMCON and MODSIM can identify. The application reads this file
while starting and prepares a list containing the type numbers and descriptions of each
TYPE in the file.
. Since the application is based on an eventdriven approach, discussed in the
Section 3.2, there is no sequence of steps to follow to configure a simulation. However, a
task tree is presented in the Appendix B to help build a simulation. The task tree may be
used as a checklist to ensure that the simulation is completely configured. The tree lists
the processes to be performed and shows options, in square brackets, on how to perform
them.
4.5. Workspa£e
The main window consists of a comprehensive menu that supports the important
operations. The Figure 4.51 shows the available menu and toolbar options. Operations
such as opening, saving and printing a workspace can be performed using the File menu.
The Edit menu options provide facilities to perform edition operations such as cut, copy,
paste, delete and find. The Edit menu options work with units, boundary conditions and
connections selected by the user. For large and complicated simulations, involving many
components, the Hide All Connections and Show All Connections menu options can be
used to hide or show the connections in the workspace. The Settings menu can be used to
set. the simulation requirements, the locations of MODSIM and SLIMCON and the serial
68
number for the workspace. The Boundary File Editor can be acces ed through the
BoundaryConditions menu. A simulation can be initiated and executed using the
Simulation menu. The simulation output can be plotted using the options in the Reports
menu.
BoundaryConditions
Se1tin . Menu
S~ Requirements.•.
status Irtorrnatton
File Locations .
5ave... ], SerieI Nl.IOOer .
Save As...
Create5IM Fie, ..
'~~:""':~=I
ptilt Workspace (In!aQe)
Prri Workspace. (Text)
File Menu
Edit Menu
elt
~opy
p~
Delete
SimulaliOIl
RI.I;J....
Reports Menu
PkJI; Silnl,etian Reds...
Plot Res1dJals""
Flri!:l CQrnppnent••• ClrI+F
Hde AI CorY,'l«tlons
SIlO\lI AJ CorneetiOnSFigure
4.51. Inteiface menu and toolbars
The Types component library and the Boundary conditions library shown in the
Figure 4.51 can be used to insert HVACSIM+ models and boundary conditions in the
workspace, respectively. These components are visually represented as icons in the
workspace. The Figure 4.52 shows the various components and their icon descriptions
using an example.
69
: Connection line
Green  Pressure
Blue  Flow
Magenta  Temperature &
Power
Btack Control signal,
Rotation Rate &
Humidity
Cyan  Energy
Variable Category
p~ Pressure····
WFlow
T  Temperature
C  Control Signal
R  Rotation Rate
E  Energy
Q  Power
H  Humidity /
/
/
Output Port t .... ... .
._"
Reported Variable;
, Contains the TYPE
number and the
description of the model
Figure 4.52. Figure showing the description ofcomponent icons
The user can view or change properties of units in the workspace. Properties that
can be set include parameter values, and initial values for output and input variables.
Initial values for the input variables can be set using constant boundary conditions. The
Figure 4.53 shows the dialog where the parameter values and initial values of an
example unit have been set.
70
.. PlOperhe. of INLET CONDUIT [DUCT Oil PIP[J £3
...... v............
INlET FLUJD OUTl.ET FLUID
'RaSURE 'RESSURE
INLEl'FLUlD
.ERAlVlIE
.............,..,
I .. 25.2
_I8lTIlIR
~E
ORa)
0.015
Figure 4.53. Dialog used to set the parameters ofunits
Two kinds of boundary conditions can be inserted in the workspace, namely, the
constant boundary conditions and timedependent boundary conditions. The interface
provides an efficient way of attaching timedependent boundary variables to columns in
the boundary file. The Figure 4.54 shows the dialog where the properties of a boundary
condition can be set or modified.
• ~.!'den!~~&1d~
Constltlll Bounde/y Conaltion
Figure 4.54. Dialog used to set the properties of boundary conditions
71
The workspace features an efficient and effective error handling. Suitable error
messages are displayed when the user makes an error. The Figure 4.55. shows an error
message generated while making a bad connection.
TtlfllPerdture bour1dery, condtioils shoUd be clllYlected only to Temperature input ports.
Figure 4.55. Error message generated while making a bad connection
4.6. Boundary File Editor
This section describes the features of the boundary file editor and its subcomponents.
The boundary file editor represents boundary files as spreadsheets. It provides
functions to create and extract data from boundary value files into the workspace. Edition
operations like cut, copy and paste can be perfonned on the data in the editor. The editor
supports copy and paste operations between external applications like Microsoft Excel. In
addition, several mathematical functions, such as sine, cosine, average, maximum, and
minimum, are available to be applied on the data. The Figure 4.61 shows a section of the
boundary file editor containing sample data.
Figure 4.61. Section ofa boundary file editor
72
4.7. Simulation
This section describes the features of the simulation manager and its subcomponents.
The program efficiently manages the execution of SLIMCON and MODS1M in
the DOS shell and tries to interpret the errors generated. The Figure 4.71 shows a
simulation while it is running. It shows the DOS shell overlaid on the simulation
manager. Efficient error handling is performed to interpret any errors generated. The
Figure 4.72 shows an error message generated when boundary values are not identified
while running a simulation. There 'is also no need for the user to open the output fIles and
study them manually since the interface can plot the graphs immediately aft.er the
simulation ends.
Runnmg SImulation  Smart Brtdge 5y"tem DesIgn p,'oledure step I : ".; ;rr
,,:). .
"Tbe toDOwint S~P$;are being executed. as the
simulatio.ns l1ro en s. Click tit STOP button to
stop the simuJ,atiall at au time.
Checlci,ng envitonmept selbDgs.
Figure 4.71. Figure showing the concurrent working ofDOS shell and the interface
73
Running SImulation  Smart Bndge System Design Procedure Step I . ,~;;;rf'
The (onome steps are beine: pucat d. th
simulations pro se. Click the STOP button to
stop the simulation at any timP.
Checking enviromnent settings.
Generaticg the simu1atI.on (SIM) &le.
Setting environment for SLIMCON.
Executing SLIMCON.
Checlcing for errors in the definition file.
Setting en'ilironmentf'i r MODSllv1
ExecutingMODSIM (to stop close MODSIM).
Ending the simulation (mocDfying the output file).
/No b01Dldin)' value data in the wOl·klipace.
not create tile b01Uuhll)' fill'.
Figure 4.72. Error message while running a deficient simulation
4.8. Plotting
This section lists the features of the plotting routines and its subcomponents.
• The tool allows the user to plot multiple reported variables in one graph. This
feature can be used to study the relationships between variables in the system.
• The tool can read output files that were created from a previous simulation. It
can also plot data from these files.
• The interface facilitates plotting residuals both through a simulation as well as
from a previously created Diagnostic information file (. DIG).
• Routines are provided in the tool to print the plots. These routines also allow
the user to scale the plots on the paper before printing them.
74
The Figures 4.81 shows the forms used for plotting simulation results. from the
output files and residuals, from the Diagnostic Information files.
~ Plol SImulation Res,dlJdl. ..
, to •
~
•• _ i
Semple Tme 1 I Sllmple Tme 2 I Sa.TIlTlll 3 i
Time (Seconds)
Visual Modeling Tool For HVACSIM+ •
Sample Residual Plot
ri!
iIIII
I
II
III!
 Sample 1
 Sample 2
15
10
.Il' S
" .:(
> 0
5
10
X·Alds
Visual Modeling Tool For HVACSIM+·
Sample Graph
I~ Plol, Analyze Slmulatlon Rewlh ;:;.
Eie
Open a Diagnostic Information (.DIG) file from 'he menu to I
plot lI1e residuals. NOTE: Each bar represents an tieratio.n..._._.
J1
for a particular time step. Also. use the scroll bar to scroll
through lime values
..._._...__. .._...__ .._..__.
Figure 4.81. Forms usedfor plotting simulation results and residuals
4.9. Installation and maintenance
An installation script was developed using the Packaging wizard supplied with
Microsoft