  
 
 
ABSTRACT 
 
 
 
 
Title of Dissertation: AN OBJECT ORIENTED SIMULATION 
FRAMEWORK FOR STEADY-STATE ANALYSIS 
OF VAPOR COMPRESSION REFRIGERATION 
SYSTEMS AND COMPONENTS 
 
  
 David Richardson, Doctor of Philosophy, 2006 
  
Dissertation Directed By: Professor Reinhard Radermacher 
Department of Mechanical Engineering 
 
 
Thermo fluid energy system simulation has shown to be a useful tool for 
engineers; encompassing component design, system design and with increasing 
interest, system optimization.  Thermo fluid energy systems, be they for comfort 
cooling, comfort heating, power generation, or any other purpose typically possess a 
unique composition and function.  This has resulted in simulations for individual 
  
rather narrowly defined energy systems, each customized for the particular system of 
interest. However, it is impossible to ignore that the majority of thermo fluid energy 
systems share, among others, the common characteristics of fluid flow, mechanical 
work input/output and energy input/output via heat transfer.  This dissertation exploits 
this similarity, and develops an object oriented methodology for modeling 
components and solving systems created from such components, operating in steady-
state.  The technique is novel in that it discriminates between systems, and their sub-
systems, referred to as components.  This methodology serves as a functional starting 
point which will appeal to the objectives of individual research groups, such as 
industrial sponsors, academic professionals, and students.  The dissertation then 
presents several examples highlighting the major points in the analysis, and a 
complex example that demonstrates where such a tool may be usefulness in a product 
design environment.  Lastly, the dissertation presents a component based, user-
friendly interface specifically for vapor compression refrigeration systems.  Several 
examples are used to validate the component models, reproducing experimental data 
reasonably well within a range of 5% for most performance variables.   
 
 
 
 
 
 
 
  
 
 
 
 
 
 
 
 
 
 
 
 
 
AN OBJECT ORIENTED SIMULATION FRAMEWORK FOR STEADY-
STATE ANALYSIS OF VAPOR COMPRESSION REFRIGERATION SYSTEMS 
AND COMPONENTS 
 
By 
David H. Richardson 
 
 
 
Dissertation submitted to the Faculty of the Graduate School of the  
University of Maryland, College Park, in partial fulfillment 
Of the requirements for the degree of 
Doctor of Philosophy 
2006 
 
 
 
 
 
 
 
Advisory Committee: 
Professor Reinhard Radermacher, Dissertation Chair/Advisor 
Professor Shapour Azarm  
Doctor Yunho Hwang 
Professor Gregory Jackson 
Professor Herbert Rabin 
 
 
  
 
 
 
 
 
 
 
 
 
 
? Copyright by 
David H. Richardson 
2006 
 
 
ii 
Dedication 
This work is dedicated to my very special friends Glenn Campopiano and Adam 
Carter.  Rest peacefully; knowing that regardless of my state of mind, you are always 
there to assure me that I can do this.  Understand that even though you could not stick 
around to see it done, you have been here every step of the way. 
 
iii 
Acknowledgements 
I would like to acknowledge Dr. Reinhard Radermacher, and the entire Mechanical 
Engineering Department at the University of Maryland College Park.   
I first thank Dr. Radermacher first and foremost, for affording me the opportunity to 
pursue my studies at the University of Maryland, and for his infinite patience as I found 
my way through the process. 
I also thank the other members of my dissertation committee, Dr. Azarm, Dr. Hwang, 
Dr. Jackson and Dr. Rabin.   
I would like to thank Dr. Ugo Piomelli, for listening to me above and beyond the call 
of Fluid Mechanics, and for lending advice of a personal nature. 
Lastly, I would like to thank Ronnie Brown, who helped me find my way through a 
very dark place. 
 
iv 
Nomenclature 
Symbols 
A 
Cross-Sectional / 
Surface Area 
m
2
 
C Heat Capacity Rate W/K 
c* 
Heat Capacity Rate 
Ratio 
- 
CF Correction Factor - 
COP 
Coefficient of 
Performance 
- 
C
p
 Specific Heat J/kg-K 
DP
Factor
 
Pressure Drop 
Correction Factor 
- 
d Diameter m 
E
 
 
Energy Flow W 
E
v
 Viscous Dissipation W 
e Electrical Charge C 
F
 G
 
Force N 
FR Fin Ratio - 
f Friction Factor - 
g 
Gravitational 
Acceleration 
m/s
2
 
G Gibbs Free Energy W 
h Specific Enthalpy J/kg 
h 
Heat Transfer 
Coefficient 
W/m
2
-K 
k Thermal Conductivity W/m-K 
L Length m
m Mass Kg 
m   Mass Flow Rate kg/s 
? Viscosity Pa-s 
? Chemical Potential J 
NTU 
Number of Transfer 
Units 
- 
Nu Nusselt Number - 
P Property -
p,P Pressure Pa, Psia 
P0 Environmental Pressure Pa, Psia 
PLF Pseudo Length Factor - 
Q  
Characteristic Energy 
Flow Rate 
W 
,QQ
 
 
Heat Transfer/ Heat 
Transfer Rate 
J/W 
r residual W,kg/s 
R Thermal Resistance 
m
2
-
K/W 
Re Reynold?s Number - 
s Specific Entropy J/kg-K 
S Surface Area m
2
 
SEER 
Seasonal Energy 
Efficiency Ratio 
Btu/W-
hr 
SH Superheat K 
T Temperature K,?C,?F
T0 
Environmental 
Temperature 
K 
u 
Specific Internal 
Energy 
J/kg 
U 
Overall Heat Tranfer 
Coefficient 
W/m
2
-K 
V Electrical Potential V 
V Volume m
3
 
v
 G
 Velocity m/s
,WW
 
 Work/Power J/W 
x 
Thermodynamic 
Quality 
- 
z Elevation m 
 
Greek  Letters 
?  Deviation W,kg/s 
? Change - 
?  Relative Roughness m 
?  Convergence Criterion W,kg/s 
?  Effectiveness -
?
s 
Isentropic Efficiency - 
?
v
 Volumetric Efficiency - 
? Viscosity Pa-s 
? Chemical Potential J 
? 
Gravitational Potential 
Energy 
J 
? Humidity Ratio g/kg 
 
Subscripts/Superscripts 
1,2 Weak/Strong Fluid 
Air 
Air 
Cond Condenser
Des Desuperheater 
i/o In/Out
J Junction 
LHT Liquid Heat Transfer 
PD Pressure Drop
TPHT Two-Phase Heat Transfer 
p Port
ref 
Refrigerant 
Sub Subcooler 
t Transported 
* Characteristic Value
 
 
v 
Table of Contents 
Dedication..................................................................................................................... ii 
Acknowledgements......................................................................................................iii 
Nomenclature............................................................................................................... iv 
Table of Contents.......................................................................................................... v 
List of Tables .............................................................................................................. xii 
List of Figures............................................................................................................ xiv 
1. Introduction....................................................................................................... 1 
1.1 Vapor Compression Refrigeration and Simulation and Research Motivation 1 
1.2 Literature Review............................................................................................ 2 
1.2.1 Generalized Simulation........................................................................... 2 
1.2.2 Air-Conditioning/Refrigeration Simulation............................................ 4 
1.2.3 Existing Software Packages.................................................................. 10 
1.2.4 Summary............................................................................................... 11 
1.3 Research Objective ....................................................................................... 11 
1.4 References..................................................................................................... 13 
2. General Thermo-Fluid Simulation Structure .................................................. 18 
2.1 Introduction................................................................................................... 18 
2.2 A Note on Thermo Fluid Properties.............................................................. 19 
2.3 Conservation Equations ................................................................................ 19 
2.4 Simulation Data Structures ........................................................................... 20 
2.4.1 Selection of Independent Variables ...................................................... 20 
2.4.2 Components .......................................................................................... 22 
 
vi 
2.4.3 Ports ...................................................................................................... 23 
2.4.4 Junctions ............................................................................................... 23 
2.4.5 Junction Evaluation............................................................................... 24 
2.5 Constructing a Thermo Fluid System ........................................................... 25 
2.6 Controlling the System ................................................................................. 29 
2.7 Simulation Hierarchy.................................................................................... 31 
2.8 Simulation Framework Design Rules ........................................................... 31 
2.9 Summary....................................................................................................... 32 
2.10 References................................................................................................. 33 
3. The Junction Solver ........................................................................................ 34 
3.1 The Junction Component Port Matrix (JCP) ................................................ 35 
3.2 Component Specification.............................................................................. 37 
3.3 Independent Variables and Residual Equations............................................ 37 
3.4 Solution of the Model ................................................................................... 39 
3.4.1 Simultaneous Approach........................................................................ 39 
3.4.2 Fractional Step ...................................................................................... 40 
3.4.3 Solving Hydraulic Relations................................................................. 44 
3.4.4 Solving Thermal Relations ................................................................... 44 
3.5 Conditioning of the Model............................................................................ 44 
3.6 The Solution Seed......................................................................................... 46 
3.7 Summary....................................................................................................... 46 
3.8 References..................................................................................................... 47 
4. Simulation of a Heat Exchanger ..................................................................... 49 
 
vii 
4.1 Model Components....................................................................................... 49 
4.1.1 Tubes..................................................................................................... 49 
4.1.1.1 Hydraulic Relations ......................................................................... 51 
4.1.1.2 Thermal Relations............................................................................ 52 
4.1.1.3 Refrigerant Side Heat Transfer Coefficients.................................... 53 
4.1.1.4 Air Side Heat Transfer Coefficients ................................................ 54 
4.1.1.5 Air Surface Area Calculations ......................................................... 55 
4.1.1.6 Distributed Models........................................................................... 55 
4.1.2 A Refrigerant Circuit ............................................................................ 56 
4.1.3 An Air Circuit ....................................................................................... 56 
4.2 The Actual Heat Exchanger.......................................................................... 58 
4.3 The Heat Exchanger Model .......................................................................... 60 
4.3.1 Wrapping the Junction Solver............................................................... 63 
4.3.2 Seeding the Junction Solver for the Heat Exchanger............................ 64 
4.3.3 Scaling Parameters for the Model......................................................... 64 
4.4 Validation...................................................................................................... 67 
4.4.1 Thermophysical Property Distributions................................................ 67 
4.4.2 Mass Balances and Energy Balances.................................................... 69 
4.5 Verification ................................................................................................... 71 
4.6 Summary....................................................................................................... 75 
4.7 References..................................................................................................... 76 
5. VapCyc ? Vapor Compression Refrigeration System .................................... 77 
5.1 VapCyc Component Specifications .............................................................. 78 
 
viii 
5.2 The Test System............................................................................................ 80 
5.3 Simulation of Vapor Compression Refrigeration ......................................... 81 
5.4 Premise.......................................................................................................... 82 
5.5 VapCyc?s Component Library...................................................................... 83 
5.5.1 A Generic Compressor Model .............................................................. 83 
5.5.2 The Alpha Compressor ......................................................................... 84 
5.5.3 The Beta and Gamma Compressor ....................................................... 86 
5.5.4 The Delta Compressor .......................................................................... 86 
5.5.5 ARI 10 COEFFICIENT MODEL......................................................... 87 
5.5.6 The Generic Orifice .............................................................................. 87 
5.5.7 Generic Heat Exchanger ....................................................................... 88 
5.5.8 The Alpha Heat Exchangers ................................................................. 90 
5.5.9 The Eta and Epsilon Condensers .......................................................... 91 
5.5.10 Heat Exchanger Summary ................................................................ 93 
5.6 VapCyc Utility.............................................................................................. 93 
5.6.1 Charge Management ............................................................................. 94 
5.6.2 Component Inter-Changeability............................................................ 94 
5.6.3 Multiple Refrigerants............................................................................ 94 
5.6.4 Multiple Unit Systems .......................................................................... 95 
5.6.5 Physical Variables................................................................................. 95 
5.6.6 Economic Variables.............................................................................. 95 
5.6.7 Seasonal Performance........................................................................... 96 
5.6.8 Interactive Mode ................................................................................... 96 
 
ix 
5.6.9 Automated Mode: Parametric Analysis................................................ 98 
5.7 VapCyc System Model ................................................................................. 98 
5.7.1 Scaling the Independent Variables and Balance Equations.................. 98 
5.7.1.1 The System Context....................................................................... 100 
5.7.1.2 The VapCyc Air-Conditioning Context......................................... 105 
5.7.1.3 The Multiple VapCyc Contexts ..................................................... 106 
5.7.1.4 Scaling With the Context............................................................... 108 
5.7.2 Control of VapCyc.............................................................................. 108 
5.8 VapCyc Examples of Capabilities and Verification ................................... 109 
5.8.1 Verification of the Model.................................................................... 110 
5.8.2 Verification of the Scaling.................................................................. 112 
5.8.3 Examining the VapCyc Context ......................................................... 114 
5.8.4 Comparison with Test Data ................................................................ 117 
5.8.5 Optimization of Charge....................................................................... 122 
5.8.6 SEER Example.................................................................................... 124 
5.8.7 Application Examples: Expandability of the Simple System............. 125 
5.8.7.1 Evaporators in Parallel................................................................... 125 
5.8.7.2 Compressors in Series.................................................................... 127 
5.8.8 Parametric Study of COP With Respect to Charge and Air Velocity. 129 
5.9 Summary..................................................................................................... 130 
5.10 References............................................................................................... 133 
6. Conclusions................................................................................................... 135 
6.1 Thermal Energy System Simulation Framework........................................ 135 
 
x 
6.2 Basic Components ...................................................................................... 136 
6.3 Component Models..................................................................................... 136 
6.4 Compound Energy Systems........................................................................ 137 
6.5 VapCyc ....................................................................................................... 138 
6.6 Future Work................................................................................................ 141 
6.7 Major Accomplishments............................................................................. 141 
6.8 Major Engineering Accomplishments ........................................................ 143 
A1. TEST System ................................................................................................ 147 
A1.1 The TEST Compressor Model................................................................ 147 
A1.2 The TEST Orifice Model........................................................................ 148 
A1.3 The TEST Evaporator Model.................................................................. 149 
A1.4 The TEST Condenser Model .................................................................. 151 
A1.5 Raw Data................................................................................................. 152 
A1.6 References............................................................................................... 156 
A2. VapCyc ? The Program ................................................................................ 157 
A2.1 Basic Description.................................................................................... 157 
A2.2 The Interface ........................................................................................... 158 
A2.2.1 Choosing Components........................................................................ 160 
A2.2.2 Editing Component Independent Variables ........................................ 161 
A2.2.3 The Parametric Interface..................................................................... 162 
A2.3 The Solution Routine .............................................................................. 163 
A2.3.4 The Pre Processor Phase ..................................................................... 164 
A2.3.5 The Solution Phase ............................................................................. 164 
 
xi 
A2.3.6 The Post Processor Phase.................................................................... 164 
  
 
xii 
List of Tables 
Table 2-1: Equations used to Model Generalized Thermo Fluid Energy Systems........... 28 
Table 2-2: Equations used to Model Generalized Thermo Fluid Energy Systems with 
Negligible Kinetic and Potential Energies................................................................ 29 
Table 3-1: JCP for a Vapor Compression Refrigeration System...................................... 36 
Table 3-2: A List of Necessary Properties for all Components Utilized in the Junction 
Solver ........................................................................................................................ 37 
Table 4-1: Test Data for R-404a Condenser..................................................................... 59 
Table 4-2: Test Data for R-290 Condenser....................................................................... 60 
Table 4-3: Correction Factors for Modeled Heat Exchanger (R-404a) ............................ 62 
Table 4-4: Correction Factors for Modeled Heat Exchanger (R-290).............................. 62 
Table 5-1: VapCyc Component Properties....................................................................... 79 
Table 5-2:Alpha Compressors Rated For Air-Conditioning............................................. 85 
Table 5-3: Alpha Compressors Rated For Low-Temperature Refrigeration.................... 85 
Table 5-4: Eta Condenser Nominal Operating Point for Air-Conditioning...................... 91 
Table 5-5: A Summary of VapCyc Heat Exchangers....................................................... 93 
Table 5-6: VapCyc Evaporating and Condensing Temperatures ................................... 107 
Table 5-7: Scaled Independent Variables and Balance Equations.................................. 109 
Table 5-8: Mass and Energy Balances for Test Configurations ..................................... 110 
Table 5-9: Air Temperatures of A/C and Refrigeration Context.................................... 115 
Table 5-10: Measured vs. Simulated Data for R-290 Test System................................. 121 
Table 5-11: Measured vs. Simulated Data for R-404a Test System............................... 122 
Table 5-12: Alpha 3600 System  Performance at ASHRAE "A" and "B" Conditions... 124 
 
xiii 
Table 5-13: Results for Nominal 20 kW Alpha Split Air Conditioning System ............ 127 
 
 
xiv 
List of Figures 
Figure 2-1: Energy System Components .......................................................................... 22 
Figure 2-2: A Junction ...................................................................................................... 25 
Figure 2-3: Open (a) and Closed (b) Energy Systems ...................................................... 26 
Figure 3-1: A Modified Vapor Compression Refrigeration System................................. 36 
Figure 3-2: Flow Chart for Simultaneous Equations ........................................................ 41 
Figure 3-3: Fractional Step Method of Jiang (2003) ........................................................ 42 
Figure 3-4: A Modified Fractional Step for the Junction Solver...................................... 43 
Figure 3-5: Successive Substitution Energy Routine........................................................ 44 
Figure 4-1: A Refrigerant Circuit ..................................................................................... 56 
Figure 4-2: An Air Stream: Air Flows from Left to Right, Through Three Air Segments, 
Simulating Air Flow Over Three Tube Segments .................................................... 57 
Figure 4-3: Cross Section of Heat Exchanger Configuration.  Air Flows Across Banks of 
Tubes from Left to Right .......................................................................................... 58 
Figure 4-4: Heat Exchanger Physical Parameters............................................................. 59 
Figure 4-5: Geometry Modification For Heat Exchanger Model.  The Real Condenser on 
the Left is Approximated with Circuitry on Right.................................................... 61 
Figure 4-6: Heat Exchanger Solution Flow Chart ............................................................ 63 
Figure 4-7: Algorithm for a Linear Pressure Distribution ................................................ 64 
Figure 4-8: Condition Number of Jacobian Matrix versus Hydraulic (Inner) Iteration 
Number.  The Iterative Hydraulic Routine is Wrapped in a Successively Substituted 
Outer Thermal Routine.  Scaling of the Independent Variables and Residual 
 
xv 
Equations Reduces the Condition Number of the Jacobian Matrix of the Hydraulic 
Routine...................................................................................................................... 66 
Figure 4-9: Refrigerant (a) Pressure and (b) Quality Distribution.  Air Flow, Not Shown 
In the Figure, Flows from Left To Right Over the Tube Banks, Cooling Them. ..... 69 
Figure 4-10: Heat Exchanger Temperature Distribution ? ............................................... 70 
Figure 4-11: Heat Exchanger (a) Mass Balance and (b) Energy Balance ........................ 70 
Figure 4-12: Mass Flow Rate Comparison for R-404a Condenser (See Table 4-1)......... 72 
Figure 4-13: Heat of Rejection Comparison for R-404a Condenser (See Table 4-1) ...... 73 
Figure 4-14: Subcooling Comparison for R-404a Condenser (See Table 4-1) ................ 73 
Figure 4-15: Mass Flow Rate Comparison for R-290 Condenser (See Table 4-2) .......... 74 
Figure 4-16: Heat of Rejection Comparison for R-290 Condenser (See Table 4-2) ........ 74 
Figure 4-17: Subcooling Comparison for R-290 Condenser (See Table 4-2) .................. 75 
Figure 5-1: ?
s
 versus Compression Ratio for Delta Compressor...................................... 86 
Figure 5-2: The Alpha 3600 System in Air-Conditioning Mode...................................... 97 
Figure 5-3: Changing of the Compressor.......................................................................... 97 
Figure 5-4: Simple Vapor-Compression Refrigeration System........................................ 99 
Figure 5-5: A Compressor Model Number (Tecumseh, 2004)....................................... 100 
Figure 5-6: ARI 540-99 Rating Points for Compressors ................................................ 101 
Figure 5-7: ARI 540-99 Heating and Cooling Compressor Operating Points................ 102 
Figure 5-8: Performance Summary of a Copeland CS10K6E-PFV Compressor (Copeland 
Web Site) ................................................................................................................ 103 
Figure 5-9: Performance Summary of a Copeland ARJ11K1E-PFV Compressor......... 104 
Figure 5-10: A Nominal Air-Conditioning Operating Point for VapCyc....................... 106 
 
xvi 
Figure 5-11: A Nominal Air-Conditioning Operating Point for VapCyc....................... 107 
Figure 5-12: P-h Diagrams for Various Alpha 3600 Configurations ............................. 111 
Figure 5-13: P-h Diagrams of Numerical Experiment Systems ..................................... 113 
Figure 5-14: Jacobian Condition Number vs. Iteration Number .................................... 113 
Figure 5-15: Jacobian Condition Number vs. Iteration Number .................................... 114 
Figure 5-16: Performance of Alpha3600 Components:.................................................. 115 
Figure 5-17: Low Temperature Solution to an Alpha Hybrid System............................ 116 
Figure 5-18: Simulated vs. Measured Evaporator Heat Load......................................... 118 
Figure 5-19: Simulated vs. Measured COP .................................................................... 118 
Figure 5-20: Simulated vs. Measured Mass Flow Rate .................................................. 119 
Figure 5-21: Simulated vs. Measured Refrigerant Charge ............................................. 120 
Figure 5-22: Suction and Discharge Pressure for R290 Test System............................. 120 
Figure 5-23: Measured and Simulated COP vs. Condenser Sub-cooling....................... 123 
Figure 5-24: ASHRAE "A" and "B" Rating Conditions for Air-Conditioning Systems 124 
Figure 5-25: Alpha 3600 System SEER vs. System Charge .......................................... 125 
Figure 5-26: A Nominal 20 kW Alpha Split Air Conditioning System ......................... 126 
Figure 5-27: System COP versus Ambient Temperature ............................................... 129 
Figure 5-28: COP and Capacity versus Charge and Air Velocity .................................. 130 
 
 
1 
1. Introduction 
1.1 Vapor Compression Refrigeration and Simulation and Research 
Motivation 
In the global economy of the twenty first century, energy production, distribution and 
conversion are of paramount importance.  Environmental acts, such as the Montreal 
Protocol and the Kyoto Protocol have demonstrated a global initiative toward protecting 
the environment from man made products which may cause damage.  A large part of the 
world?s energy usage relates to heating, ventilation, air-conditioning and refrigeration 
(HVAC/R) systems.  Consequently, research is motivated to meeting simulation needs 
for component manufacturers, system manufacturers and researchers. 
A simulation framework, catered toward vapor compression refrigeration and its 
components is a useful tool for research.  Such a Framework that is general in nature 
would be very helpful.  Ideally, the framework facilitates the development and use of 
libraries of individual component models that integrate into systems seamlessly, and 
without a priori knowledge of the system to which they are being inserted. 
The generality of the vapor compression refrigeration framework allows future 
generations to utilize it without specific knowledge of the inner workings of the 
framework.  This prevents generation after generation of re-invention of the simulation 
structure, which hampers the development of the research of individual components 
and/or system configurations. 
 
2 
1.2 Literature Review 
1.2.1 Generalized Simulation 
A variety of engineering problems encountered in practice exhibit behavior that is 
predicated upon the interaction of complex physical phenomena, i.e. their response is 
determined not by a single phenomena, but by the multidisciplinary interaction (Gupta 
and Meek, 2000). Considering individual components as a ?discipline?, combines 
different components and estimating their performance as a system is multidisciplinary in 
nature. 
Object oriented simulation allows for increased programming productivity in code 
maintenance and re-usability.  A single generalized simulation technique can be used to 
simulate many different individual energy systems, which can be assembled into 
coordinating multidisciplinary energy systems.  Applying object oriented simulation to 
the standard vapor compression refrigeration system, each individual system component; 
compressor, condenser, expansion device and evaporator can be modeled independently 
from the coordinating system simulation tool which will integrate them to a complete 
system. 
This approach is beneficial, as modular component models which can utilize a 
majority of the same code.  A secondary benefit is realized when a generalized solution 
routine is used within a coordinating system model. 
Along with these object oriented techniques comes the potential benefit of 
components conforming to a single standard for their use in the generalized simulation; 
namely an industry wide component library.  This industry wide component library 
 
3 
coupled with a simulation tool to utilize can then estimate the performance of systems 
comprised of these individual components. 
Multidisciplinary optimization, or MDO, (Alexandrov and Hussaini, 1997) may be 
defined as a methodology together with a set of tools for assisting in the design of 
complex coupled systems, that is systems whose behavior is governed by many distinct 
but interacting physical phenomena.  An object oriented, simulation tool for thermo fluid 
energy systems is a logical foundation for an industry wide analysis and optimization tool 
which may be used for MDO. 
Although moderate work has been performed in generalized energy system 
simulation, little has focused in on vapor compression refrigeration.  In 1994, Paulus et 
al. (1994) describes a technique using vectors and matrices to specify a system.  The 
technique is lacking because it fails to discriminate between the component equations, 
and system property balances.  This in effect creates a single system of equations which 
must be solved simultaneously.  Performance relations are hard coded into the program, 
and there is no flexibility to change components. 
A generalized approach (Majumdar et.al, 1995), initially developed to simulate the 
fluid flow systems found in rocket propulsion became a very comprehensive program for 
solving generalized flow networks.  The interface for the program was somewhat 
cumbersome, in that the user interacts with the program through many input and output 
files, as well as through an interactive mode.  As with Paulus, the program is configured 
for a single system.  A change in make-up of the system is achieved through the input of 
a different set of input files, or an entirely new run in the interactive mode.  The use of 
 
4 
input/output files, coupled with component models which are integrated into the 
simulation itself are limiting factors in achieving a flexible modeling environment. 
In 2000, Lindsay (2000) offers ?Network Builder?, which is a foundation for the 
generalized energy system simulation technique given here.  The formulation is based on 
a series of components connected by junctions, and serves as a starting point for this 
work. 
Jiang (2003) creates an object oriented heat exchanger simulation, based upon 
Lindsay?s techniques.  The simulation is comprehensive and flexible, allowing many 
parameters to be specified via a user interface.  In this simulation, tubes act as individual 
components, which can be combined into different configurations by the user.  Although 
general in nature, the only components which could be changed were tubes. 
1.2.2 Air-Conditioning/Refrigeration Simulation 
The simulation of the vapor compression refrigeration system is clearly a benefit to 
the engineer, and it is therefore not surprising that great deal of work has been done 
toward creating good simulations.  Meaningful simulations of refrigeration systems can 
be traced back to as early as 1976, when Hiller and Glicksman (1976) as well as Davis 
and Scott (1976) offered the first modern looks at simulation, which were limited in 
scope to air-to-air heat pump simulations.  Ellison and Creswick in 1978 then used the 
Hiller and Glicksman model to examine the change in system performance due to 
changes in components.  Also appearing around the same time was a simulation for 
optimizing the heating function of heat pumps by Carrington (1978).  In 1979, Ellison 
and Rice gave a report on the Oak Ridge heat pump model, a model that continued to 
evolve over the years, and is indeed still existing and evolving today. 
 
5 
The early 1980?s gave birth to many refrigeration simulations.  Rice et al. (1981) used 
computer simulation to explore optimization.  This was achieved by combining the Oak 
Ridge model with an optimization program to optimize system level performance as a 
function of specific component level variables.  Also in the early 1980?s were McMullan 
et al. (1981) used simulation to compare field performance with (laboratory) test 
performance.  Domanski and Didion in 1983 created a model (HPSIM) for the United 
States National Bureau of Standards (NBS-Now known as the National Institute for 
Testing and Standards, NIST).  The HPSIM model is still in existence today, and has 
undergone many enhancements over the years which allow it to be very flexible.   
Refrigerant mixtures were introduced to vapor-compression models in the mid 
1980?s.  Connon (1984) presents a simple simulation model in 1984, utilizing refrigerant 
mixtures.   
Parise (1986) offers a simulation model that is based upon the governing equations 
for its component models, rather than performance maps and empirical equations.  Parise 
claims at this point that empirical performance maps are responsible for deviations of 
simulations from actual equipment behavior. 
Welsby et al. (1988) creates a simulation model for water-to-water heat pump 
systems, and gives a literature review as the starting point, detailing the growth of the 
field.  In 1987, simulation was put to use comparing the performance of pure refrigerants 
versus mixtures, as given in McLinden and Radermacher (1987). 
Simulation literature becomes more prevalent in the 1990?s, as the personal computer 
became affordable enough to be used by more researchers 
 
6 
In 1990, Jolly et al. used simulation to present work on heat pump assisted drying.  
Also in 1990 the Oak Ridge model had evolved into CYCLEZ, or the Mark III system.  
Rice and Sand (1990) also present a modified version of CYCLEZ for a Lorentz-
Meutzner cycle refrigerator-freezer model.  In 1992 Stafanuk et al. present a model 
specifically addressing the need for a so called charge inventory, or and accounting of the 
refrigerant mass contained within the system  That work addresses the affects of the 
charge inventory, specifically the effect it has on prediction of off-design performance.  
Prior to Stafanuk only the Oak Ridge model and the HPSIM (NIST) simulation model 
accounted for system charge.  Bare (1993) presents a model designed to address 
alternative refrigerants and their performance potential.  It is noted that the paper presents 
a very simple model, as the focus is on refrigerant effect on the system, not accurate 
simulation of the system itself.  The topic of alternate refrigerants is also the focus of a 
study by Ouazia and Snelson (1994).  This is also a simple model, but is typical of system 
simulations involving alternative refrigerants at this time, which is co-incident with the 
phase out of CFCs.  Around this time the Oak Ridge model had grown to the Mark V, 
also named PUREZ (Rice et al. 1994).  
A ?toolkit? for chiller simulation is presented by Bourdouxh et al. (1994), which 
represents an ASHRAE initiative for simple simulation of HVAC components for 
building simulation.  Levins et al. (1996) used simulation to explore the effects of over-
sizing of residential air conditioning units.  The Oak Ridge Model (Mark II) was 
modified to reflect performance data from an experimental apparatus, and the results used 
to estimate energy savings by reducing compressor size in existing units.  The DOE/Oak 
 
7 
Ridge model history and overview, as well as some applications to R-22 alternatives were 
presented in 1997 at a symposium on heat pumps in cold climates (Rice et. al. 1993) 
Water chillers are a small part of the vapor compression simulation family, and oft 
ignored, but in 1998 Richardson as well as Browne and Bansal (1998) address them.  The 
study by Brown and Bansal addresses many sources of simulation literature to date, and 
concludes that although there is (to that date) a large amount of simulation research in 
existence, chillers are different in operation and therefore the field of chiller simulation 
was immature.  Browne and Bansal conclude that the two main differences between 
water chillers and typical vapor compression heat pumps are the heat exchanger 
configurations and the compressors.  It is concluded that simulations that accurately 
accounted for these differences would be very similar in nature to existing simulations.   
Anand  and Tandon  (1999) discuss two separate works, both of which focus on 
simulation of refrigeration systems.  The works involve using modeling techniques 
separating the component models from the system models. 
LeRoy (2000) uses PUREZ to predict the performance of unitary air conditioners.  
The concept of ?tuning? the PUREZ model is introduced to help the simulation match 
actual data.   
Browne and Bansal (2001) continue their modeling of chillers with an improved heat 
exchanger model.  Updates to the Oak Ridge model (now available for on-line simulation 
over the internet) can be found from the Oak Ridge web site (Rice 2001). 
An indicator of the direction is the move to component based modeling, Grossman et 
al. (2001) present ABSIM, a simulation of absorption systems, presents a modular format 
 
8 
to simulations.  Here Grossman discusses the advantages of using modular simulation 
components for the sharing of work as well as for development of the tool. 
Several key issues emerge from the study of the literature, each of which needs to be 
included into a comprehensive simulation, namely: 
? Charge Inventory 
? Architecture and solution technique of the system model 
? (Lack of) Optimization 
 
Few simulations to date account for ?charge inventory?.  Stefanuk et al. clearly 
presents (1994) that for off-design simulation the charge inventory is mandatory.  
Consequently, the majority of the existing models are only representative of a single 
operating point, not a uniquely charged system at different operating points. 
The majority of the system models in existence today integrate the components and 
the system into a single system of equations to be solved.  Nearly every cited reference 
implements this technique.  
Many solution techniques exist for the system of equations developed in simulations 
of vapor compression refrigeration systems.  The large majority of the existing 
simulations use a mixture of iteration and successive substitution.  The overwhelmingly 
popular solution technique involves solving systems of equations in successive 
substitution technique created by combining the component models with the system 
conservation laws.  Typically satisfying some of the equations solves a ?high side? 
pressure, and the results of the ?high side? outputs are then fed to the ?low side?, which 
must in turn satisfy the remaining equations.  Results from the ?low side? are fed to the 
?high side?.  The process continues until neither the high nor low sides change.  A second 
successive substitution technique (Annand, Tandon 1999) involves separate models for 
 
9 
the components.  Utilizing the charge inventory, a single component can be solved for its 
operating point and inlet and outlet flows (given an initial guess for its independent 
variables and charge).  Any discrepancy between inlet and outlet flow results in a net 
change in charge for the component.  This charge difference is passes along to the next 
component, along with the thermodynamic state of the fluid entering the component.  
This process is repeated until there is no change in the system, i.e. there are no 
discrepancies in charge or between mass flows into or out of any components.  This 
process is difficult to implement in complex systems (components in parallel for 
example) and requires fundamental modeling differences between components that store 
significant amounts of refrigerant and those that do not.  Another popular technique is to 
formulate the equations and solve them using an iterative solver such as Newton-Raphson 
(Jolly et al. 1990, Paulus et al. 1994, Bourdoux et al. 1994, Richardson 1998, Brown and 
Bansal 1998, LeRoy et al. 2000). 
Some attempts to optimize refrigeration systems as a whole, have been made (Rice et 
al. 1981, Fisher and Rice 1986,Levins et al. 1996).  However, in all of these cases the 
system simulation model is typically customized to begin with and then is either (a) 
integrated with the optimization routine, making the program even more customized, or 
(b) the optimization is performed manually with simulation program results.   
Object oriented simulation allows more flexibility to customize individual 
components, which are separated from the solution routines, which themselves are 
separated from optimization routines. 
 
10 
1.2.3 Existing Software Packages 
Many software packages currently exist, which may be used for energy system 
simulation.  Refrigeration systems definitely fit within the scope of these simulation 
packages.  The current software comes in two groups and hybridizations of the groups.  
One group is general equation solvers, which allow the user to specify the system model 
in terms of governing equations, and the software solves the set of equations.  A second 
group is advanced energy system software. 
Clearly general equation solvers are valid for any conceivable energy system and 
therefore refrigeration systems.  However, the equation solvers are extremely dependent 
upon the user to formulate the model and properly specify it, as the software has no 
intelligence with regard to the specific problem.  Examples of general equation solvers 
(or programming environments conducive to general equation solving) are: The 
Engineering Equation Solver (F-Chart), MATLAB (Mathworks), and Mathematica 
(Wolfram). 
Advanced energy system software on the other hand is often very specific to the task 
at hand (determining the operation of a refrigeration system, steam power plant, gas 
turbine power plant, etc).  Often the user is ?locked? into the specific system, and 
modification or generalization of the system is difficult or impossible.  Examples of 
advanced energy system software are: Aspen Plus (Aspen Technologies) and Gate Cycle 
(Stork). 
Hybridizations of the two groups offer fixed systems, but user defined component 
models.  An example of hybridization is Sinda/Fluint (C&R). 
 
11 
Each of the three classes of software has advantages and drawbacks, often leading a 
user to create their own custom software to specifically suit their individual needs. 
1.2.4 Summary 
Many simulations for refrigeration systems exist, and there are many software 
packages that exist to facilitate custom system models.  It is clear from the literature, that 
most often, simulations are created to fill a specific need, rather than to serve as a general 
tool.  The main output of each simulation is the performance of the specific system being 
simulated.  A drawback of previous simulations is the rigid definition of the system.  
Often existing simulations offer generality through a wide range of component 
independent variables.  For the most part, component models themselves are integrated 
into the simulation itself and therefore not subject to any change not built into the 
simulation at its creation. 
Two main challenges are inherent in current modeling of vapor compression in regard 
to creating generalized modeling tools, namely (1) a fixed systems configuration and (2) 
specific, usually fixed, component configurations. 
1.3 Research Objective 
The findings from the literature demonstrate that although the fields of energy system 
simulation in general, and vapor compression refrigeration specifically, are mature, many 
opportunities for research and progress exist today. 
A significant engineering challenge is the creation of a general vapor compression 
refrigeration system simulation tool, which can serve many individual constituent 
objectives.  This tool best serves the individual user by allowing all the detail required by 
 
12 
an individual constituent, but does not force a complete understanding of the inner 
workings of disciplines beyond the scope of this individual constituent group.  
A true general-purpose vapor compression simulation program includes all decided 
features in a set defined by components, refrigerants, refrigerant charge, system boundary 
conditions and controls.  
The primary objective of the dissertation is the development of a general-purpose 
framework for simulating fluid flow systems, such as those found in refrigeration and air-
conditioning applications.  The cornerstone of the general-purpose nature is the 
prediction of the performance of the complete set of components, working fluids, 
working fluid charge, system boundary conditions and controls, while each component of 
the system has its own independent, stand-alone simulation model. 
This thesis accomplishes the following objectives: 
1. Develop and create a simulation framework, and subsequent simulation tool in 
an object oriented manner, which provides the ability to simulate experiments 
in a laboratory.  The framework is developed with the following fundamental 
objectives: 
a. Serve as a functional starting point which will appeal to the individual 
objectives of each individual user (industrial sponsors (development of 
models and/or real equipment), academic professionals (utility of the 
simulation tool without an intimate understanding of its inner 
workings), and students (re-using of existing software to facilitate the 
growth of the project(s))) 
 
13 
b. Serve as a simulation tool such that each component within the system 
is a stand-alone simulation model, enabling the analysis and if desired 
optimization of individual simulation models and/or the system as a 
whole 
2. Verify the simulation tool with experimental data, specifically a condenser, 
and air conditioning system utilizing this condenser. 
1.4 References 
Alexandrov, N., Hussaini, M. Y., 1997, Multidisciplinary Design Optimization State of 
the Art, Society for Industrial and Applied Mathematics, 1997 
 
Anand, G., 1999, ?Transient and Steady State Model of a Household Refrigerator?, 
Masters Thesis, Department of Mechanical Engineering, University of Maryland, 1999 
 
Bare, J.C., 1993.  ?Simulation Of Performance Of Chlorine-Free Fluorinated Ethers And 
Flourinated Hydrocarbons To Replace CFC-11 And CFC-114 In Chillers?, ASHRAE 
Transactions 99 (1): pp 397-407 
 
Bourdouxh, J.H., Grodent, M, Lebrun, J.J, Saavedra, C, Silva, K, 1994, ?A Toolkit For 
Primary HVAC System Energy Calculation-Part 2: Reciprocating Chiller Models?.  
ASHRAE Transactions 100 (2): pp. 774-786 
 
Browne, M., Bansal, P., 1998, ?Steady-State Model Of Centrifugal Liquid Chillers?, 
International Journal of Refrigeration, 1998, Vol. 21, pp. 343-358 
 
Browne, M, Bansal, P., 1998, ?Challenges in Modeling Vapor-Compression Liquid 
Chillers?, ASHRAE Transactions, 1998, V. 104, Pt. 1 
 
Browne, M., Bansal, P., 2001, ?An elemental NTU-? model for vapour-compression 
liquid chillers?, International Journal of Refrigeration, 2001, Vol. 24, pp. 612-627 
 
Carrington, C.G., 1978, ?Optimizing a heat pump for heating purposes?, International 
Journal of Energy Research, Vol. 2, pp. 153-170 
 
CEEE web site, http://www.eng.umd.edu/Department_gateways/CEEE.html 
 
Clark, D. R. and May, W. R, Jr., 1985, HVACSIM+ Building Systems and Equipment 
Simulation Program Users Guide. NBSIR 85-3243. National Bureau of Standards, 
September, 1985. 
 
14 
 
Davis, G.L., Scott, T.C., 1976, ?Component Modeling Requirements For Refrigeration 
Systems Simulation?, Proceedings of the 1976 Purdue Compressor Technology 
Conference, West Lafayette, USAA, July 1976, pp. 401-408 
 
Domanski, P., Didion, D, 1983.  ?Computer Modeling of the Vapor Compression Cycle 
with Constant Flow Area Expansion Device?, NBS Building Science Series 155, May 
1983 
 
Ellison, R.D. and Creswick, F.A. 1978.  ?A Computer Simulation Of Steady-State 
Performance Of Air-To-Air Heat Pumps?, Report ORNL/CON-16, Department of 
Energy, Oak Ridge National Laboratory, Tennessee, USA 
 
Ellison, R.D. and Rice, C.K. 1979.  ?ORNL Heat Pump Model Update ? May 1979?, 
Department of Energy, Oak Ridge National Laboratory, Tennessee, USA 
 
Glicksman, L.R., 1976, ?Improving Heat?Pump Performance Via Compressor Capacity 
Control ? Analysis And Test?.  Volumes I and II, Energy Laboratory Reports MIT-EL 
76-0001 and MIT-EL 76-0002, Department of Mechanical Engineering, Massachusetts 
Institute of Technology, Massachusetts, USA (January 1976) 
 
Gupta, K.K., Meek, J.L., 2000, Finite Element Multidisciplinary Analysis, AIAA 
Education Series, American Institute of Aeronautics and Astronautics, Inc, 2000 
 
Grossman, G., Zaltash, A., 2001, ?ABSIM ? Modular Simulation Of Advanced 
Absorption Systems?, International Journal of Refrigeration, 2001, Vol. 24, pp 531-543 
 
Jiang, H., 2003, Development Of A Simulation And Optimization Tool For Heat 
Exchanger Design, Ph.D. Thesis, Department of Mechanical Engineering, University of 
Maryland, 2003 
 
Jolly, P., Jia, X. and Clements, S., 1990 ?Heat Pump Assisted Continuous Drying Part I: 
Simulation Model?, International Journal of Energy Research, 1990, Vol 14, pp. 757-
770 
 
LeRoy, J, Groll, E., Braun, J, 2000, ?Evaluating The Accuracy Of PUREZ In Predicting 
Unitary Equipment Performance?, ASHRAE Transactions, 106 (1)  
 
Levins, W., Rice, C.K., Baxter, V.D. 1996, ?Modeled And Measured Effects Of 
Compressor Downsizing In An Existing Air Conditioner/Heat Pump In The Cooling 
Mode?, ASHRAE Transactions, V 102, Pt 2 
 
Lindsay, D., 2000, ?Network Builder Notes, University of Maryland ? CEEE, Internal 
Document, July, 2000 
 
 
15 
Paulus, D. Jr., Gaggioli, R., Park, H. and Dunbar, W, 1994, ?Development of Personal 
Computer Software for Energy System Simulation?, Thermodynamics and the Design, 
Analysis and Improvement of Energy Systems, Advanced Energy Systems Division 
(Publication) AES-Vol. 33, Proceedings of the 1994 International Mechanical 
Engineering Congress and Exposition, p. 177-184 
 
Majumdar, A. K., Bailey, J. W., Schallhorn, P. A. and Steadman, T., 1999, ?Generalized 
Fluid System Simulation Program (GFSSP) Version 3.0? Sverdrup Technology MSFC 
Group Report No. MG-99-290, Nov, 1999, MSFC Science and Engineering Contract 
NAS 8-40836. 
 
McLinden, M.O. and Radermacher, R., 1987.  ?Methods for Comparing the Performance 
of Pure and Mixed Refrigerants in the Vapor Compression Cycle?, International Journal 
of Refrigeration, 10, pp. 318-325 
 
McMullan, J.T., Morgan, R., Hughes, D.W., 1981.  ?The Discrepancy Between Heat-
Pump Field And Test Performance: A Simulation Study?, International Journal of 
Energy Research, Vol. 5, pp. 83-94 
 
Metcalf, R.R., 1989, Development of a Generalized Fan System Simulation for the 
Building Loads Analysis and System Thermodynamics Program (BLAST), Masters 
Thesis, University of Illinois, 1989 
 
Ouazia, B., Snelson, P,  1994. ?Predicting System Performance Of Alternative 
Refrigerants Using A Water-Water Heat Pump?.  ASHRAE Transactions, 100(2): pp. 
140-147 
 
Parise, J, 1986.  ?Simulation of Vapour-Compression Heat Pumps? Simulation 46(2): pp. 
71-76 
 
Rice, C.K., 2001, ORNL-HPDM-MarkV-VI-Update.pps, Power Point Presentation on 
Oak Ridge Heat Pump Design Model for the web.  
http://www.ornl.gov/~wlj/hpdm/doehpdm.html, January, 2001 
 
Rice, C.K. and Jackson, W.L., 1994. ?PUREZ?The Mark V ORNL Heat Pump Design 
Model for Chlorine-Free Pure and Near-Azeotropic Refrigerant Alternatives?, Document 
Package, Oak Ridge National Laboratory 
 
Rice, C.K., Sand, J.R., 1990. ?Initial Parametric Results Using CYCLEZ?An LMTD 
Specified Lorentz-Meutzner Cycle Refrigerator-Freezer Model?, West Lafayette, PA: 
Proc. International Refrigeration Conference at Purdue University, July 17-20.  pp. 448-
458. 
 
Richardson, D. H., 1998.  ?Performance Models For Use in Computer Simulation of 
Energy Systems?, Masters Thesis, Department of Mechanical Engineering, Marquette 
University, 1998 
 
16 
 
Resource Dynamics Corporation, 2002, ?Integrated Energy Systems (IES) for Buildings:  
A Market Assessment?, Prepared for: Energy Efficiency and Renewable Energy U.S. 
Department of Energy ,Washington, DC, and Oak Ridge National Laboratory Oak Ridge, 
TN, August 2002 
 
Rice, C.K., Fischer, S. K., Ellison, R.D., Jackson, W.L. 1981.  ?Design Optimization of 
Conventional Heat Pumps: Application to Steady-State Heating Efficiency,? ASHRAE 
Transactions. Vol. 87, pp. 1037-1054 
 
Rice, C.K., Wright, L. and Bansal, P.K., 1993, ?Thermodynamic Cycle Evaluation Model 
For R-22 Alternatives In Heat Pumps-Initial Results And Comparisons?.  Moneton, 
Canada: Proc. 2
nd
 Int. Conference on Heat Pumps in Cold Climates, Aug 16-17 
 
Stefanuk, N.B.M., Aplevich, J.D., Renksizbulut, 1992. ?Modeling And Simulation Of A 
Superheat-Controlled Water-To-Water Heat Pump?, ASHRAE Transactions, 98(2) 
 
Taylor, R. D., Pedersen, C. O., Fisher, D., Liesen, R. and Lawrie, L., 1991. ?Impact Of 
Simultaneous Simulation Of Buildings And Mechanical Systems In Heat Balance Based 
Energy Analysis Programs On System Response And Control?, Proc. of the Building 
Simulation ?91 Conference, August 20-22, 1991, Sophia-Antipolis, Nice, France, pp. 
227-234, International Building Performance Simulation Association. 
 
Tandon, A., 1999, Object Oriented Modeling of Vapor Compression Systems and 
Components, Masters Thesis, Department of Mechanical Engineering, University of 
Maryland, 1999 
 
University of Wisconsin, Solar Lab, TRNSYS. 
 
Welsby, P., Devotta, S., Diggory, P.J., 1988.  ?Steady and Dynamic-State Simulations of 
Heat-Pumps. Part I: Literature Review?, Applied Energy, Vol. 31, pp. 189-203 
 
Welsby, P., Devotta, S., Diggory, P.J., 1988.  ?Steady and Dynamic-State Simulations of 
Heat-Pumps. Part II: Modeling of a Motor Driven Water-to-Water Heat-Pump?, Applied 
Energy, Vol. 31, pp. 239-262 
 
Commercial Software 
 
Aspen Plus, Aspen Technologies 
 
Engineering Equation Solver, F-Chart Software 
 
Gate Cycle, Stork Engineers & Contractors, Radarweg 60, 1043 NT Amsterdam, The 
Netherlands 
 
Mathematica, Wolfram Software 
 
17 
 
MATLAB, Mathworks 
 
SINDA/FLUINT, C&R Technologies, 9 Red Fox Lane, Littleton CO 80127, USA 
 
18 
2. General Thermo-Fluid Simulation Structure 
2.1 Introduction 
The primary objective of the dissertation is the development and creation of a thermo 
fluid simulation framework primarily for use in refrigeration applications.  A secondary 
goal is to build a simulation tool utilizing this framework which appeals to a diverse user 
group.  The framework is to be general in nature, such that the tool built upon it will 
provide the ability to simulate experiments as they may be performed in a laboratory.  
This simulation tool is designed, such that each component within the system functions 
independently, enabling the analysis and if desired optimization of component models 
and/or the system as a whole.  Creation of an object oriented tool is crucial for allowing 
users to utilize a single tool, despite having their own simulation goals, objectives and 
resources. 
By virtue of the component independence, the simulation tool is to be 
multidisciplinary in nature.  Although the components of a thermo fluid system, such as a 
vapor-compression refrigeration system are joined together, each represents a distinct 
discipline.  While a physical system is combined using pipes and a working fluid, the 
virtual system, or simulation, is joined by governing equations.  This chapter develops 
data structures and governing equations for a hierarchic thermo fluid system simulation.  
This undertaking results in a novel presentation for a generalized, multidisciplinary, 
thermo fluid system simulation. 
 
19 
2.2 A Note on Thermo Fluid Properties 
The scope of this thesis is limited to the modeling of steady state thermo fluid energy 
systems.  A further simplification is that of single substance mixtures.  That is to say 
although multi-component mixtures, such as R404a, can be used as a working fluid, the 
composition remains constant.   
  In this context, a thermo fluid energy system is a system which may allow fluid 
flow, heat transfer and/or work. 
2.3  Conservation Equations  
For a non-isothermal system, there are three conservation equations that describe the 
relations between the inlet and outlet conditions of the stream(s) (Bird et al. 1960)
1
: 
a. The mass balance (given in equation (2.1)) obtained by integrating the 
equation of continuity over the volume of the flow system. 
b. The momentum balance (given in equation(2.2)) obtained by integrating the 
equation of motion over the volume of the flow system. 
c. The energy balance (given in equation (2.3)) obtained by integrating the 
energy equation of continuity over the volume of the flow system. 
The steady state continuity equation: 
0
in out
mm=?
??
      (2.1) 
describes the (lack of) mass accumulation within the system. 
The steady state momentum equation: 
                                                 
1
 All equations are given here in their steady state form 
 
20 
( )0
tot
mp m=?? + ?v g  SF+    (2.2) 
governs the change between fluid kinetic energy, mv   , force on the fluid flow boundaries 
(product of pressure and flow surface areas), force of the control volume upon the fluid 
and gravity. 
The steady state energy equation: 
2
?
0
2
pv
umQW
?
??
??
=?? + + +? + ?
????
??
??
     (2.3) 
governs the change between fluid internal energy, m  u, flow work, m 
p
?
, kinetic 
energy, m 
2
2
v
 , gravitational potential energy, 
?
m?  , heat transfer out of the control 
volume, Q, and work done by the fluid, W. 
2.4 Simulation Data Structures 
The generalized simulation framework created in this thesis is premised upon 
individual component models communicating via a single coordinating simulation 
mechanism.  The simulation is constructed via several object oriented data structures 
which facilitate communication of information by the coordinating, or parent system.  
The data structures presented here are components, ports and junctions. 
2.4.1 Selection of Independent Variables 
The framework created in this thesis is based upon the premise that components 
satisfy their own mass, momentum and energy balances, the system, comprised of these 
components is resolved by satisfying the balances at the points in between them, later 
termed ?junctions? in this thesis. 
 
21 
Typical simulation approaches often simulate components differently from one 
another.  For example, heat exchangers and pipes are often provided independent 
variables of mass flow rate and the thermodynamic state of the inlet fluid.  In these cases, 
the component dependent variables are the thermodynamic state of the outlet fluid, such 
that the mass, momentum and energy balances are satisfied.  Contrastingly, compressors 
and valves are often provided independent variables of thermodynamic state of the inlet 
fluid and outlet pressure, while the component dependent variables are the mass flow rate 
through the component and the thermodynamic state of the exit flows, which satisfy the 
component mass, momentum and energy balances. 
Allowing components to possess different independent and dependent variables 
places a burden on an algorithm attempting to resolve a system comprised of these 
components.  Often the algorithm requires additional information about which 
components require which specific variables.  Additionally, often the evaluation order of 
the components has an affect on the solution algorithm. 
In this thesis, the generality of the framework requires that all components have the 
same set of independent variables, and calculate the same set of dependent variables. 
Choosing mass flow rate for component independent variables is problematic for a 
general routine, where a priori knowledge of system size and flow directions (which may 
be operating point dependent, rather than fixed) does not exist. 
Therefore, this thesis utilizes the working fluid pressure at all flow points, and the 
fluid enthalpy at the fluid inlet points as component independent variables.  For the single 
component mixture, the fluid pressure and enthalpy are adequate to define the fluid 
thermodynamic state.  For multi-component mixtures, a composition is also required. 
 
22 
It is noted here, that a challenge exists here, due to the mandate.  As mentioned 
previously, relations for many components typically exist in the form ? P = ? (mass flow).  
This makes the evaluation of the component mass flow rate implicit. 
2.4.2 Components 
Components are represented by appropriate engineering models which satisfy the 
macroscopic balances given in (2.1) - (2.3).  A fluid flow enters or leaves a component 
through a port, which in turn communicates with a junction. 
Figure 2-1 illustrates a typical ?black box? component as an example, and the detailed 
component, which it represents.  In this case, the component to be simulated is a 
combination of two compressors and an inter-cooler.  This component communicates 
with its environment pressure and temperature (P
0
,T
0
) through work (W
 
) and heat 
transfer ( Q
 
).  The component communicates with other components in the system 
through their respective junctions. 
Junction 1
Pressure 1
Port 1 Port 2
Junction 2
Pressure 2
P, T
00
(b)(a)
WQ
  
,
Q
 
W
 
Junction 2
Pressure 2
Junction 1
Pressure 1
P, T
00
 
Figure 2-1: Energy System Components 
Physical components are subject to certain conditions in which they operate.  
Conditions imposed upon the component at its system boundaries are referred to in this 
 
23 
thesis as boundary conditions.  Components combined into a system, are subject to both 
internal and external boundary conditions.  
 In this context, external boundary conditions are boundary conditions imposed on the 
component from outside the thermo fluid energy system which surrounds the component.  
External boundary conditions could be, but are not limited to: heat exchanger secondary 
fluid flow rates and temperature, an environmental temperature or a compressor?s speed.   
In this context, internal boundary conditions refer to the boundary conditions imposed 
on a component by the coordinating thermo fluid energy system of which the component 
is a part.  Internal boundary conditions are limited to the fluid properties from fluid flow 
interactions with other components (via the component ports and system junctions).  
Internal boundary conditions include, but are not limited to: pressure, temperature, 
enthalpy and/or a fluid mass flow. 
2.4.3 Ports 
Ports are a component?s interface to junctions in the system where fluid flow is the 
mechanism of energy transport.  The fluid flow through a port has distinct values for flow 
rate, velocity, position (in a gravitational field) and thermodynamic state. 
2.4.4 Junctions 
Junctions are points, which allow for interaction between ports via fluid flow.  Within 
the scope of this thesis, a junction is considered to be adiabatic and does not allow for any 
work interaction with its environment.  Unlike components, junctions do not posses a 
unique geometry, and therefore are characterized by a single pressure. 
 
24 
Each fluid flow into or out of a junction has a unique mass flow rate, thermodynamic 
state, velocity and height.  
The steady state solution of the simulation thus requires the mass and energy balances 
of equations (2.1) and (2.3) to be satisfied at each junction 
For a special case, where kinetic and potential energies of the individual flows are 
ignored, junctions are assumed to have perfect mixing and thus the thermodynamic state 
of all fluid flows leaving a junction is identical.  In this case, the energy balance can be 
used to resolve the appropriate value for the fluid enthalpy for each exiting fluid flow.  In 
this instance, the single value is characteristic of the junction itself and can be referred to 
as the junction?s mixing enthalpy.  
2.4.5 Junction Evaluation 
As components implicitly satisfy the conservation equations (2.1) - (2.3), the 
junctions can be used to simulate the thermo fluid system.  For example, Figure 2-2, 
illustrates a general thermo fluid energy system junction.  This particular junction has 
fluid flowing into and out of it from four different components.  The individual fluid flow 
rates (
4321
,,, mmmm      ), fluid flow velocity (v
1
, v
2
, v
3
, v
4
) and fluid flow height (z
1
, z
2
, z
3
, 
z
4
) are determined from ports connected to the junction, and are determined by the 
solution of the component?s models connected to the junction.  Both 
21
 and mm    flow into 
the junction, and thus their enthalpy values, h
1
 and h
2
 are determined by the exit ports of 
the components to which they are connected.  Both 
43
 and mm    are flowing to 
components.  The corresponding enthalpy values, h
3
 and h
4
 are independent and must be 
determined to satisfy the energy equation applied at the junction.   
 
25 
m
h
1
1
v
z
1
1
m, h, 
33
v
33
, z
m, h, 
44
v
44
, z
m
h
2
2
v
z
2
2
 
Figure 2-2: A Junction 
The independent variables of a junction are therefore its pressure and enthalpy values 
for each exiting flow.  Other junction variables, namely: fluid flow rates, velocities, 
elevation values and values of fluid enthalpy for entering streams are provided by the 
components connected to the junction, are then used to evaluate the junction mass and 
energy balances. 
As determined previously, in the absence of kinetic and potential energies, the values 
for junction exit enthalpies are identical, thus this mixed enthalpy is a property of the 
junction.  The independent variables of a junction are then its pressure and single mixed 
enthalpy value. 
2.5 Constructing a Thermo Fluid System 
Figure 2-3 shows two general energy systems; one is an open system, while the other 
is closed.  Each system contains six components, connected via junctions.  All component 
ports not communicating across a system boundary must be connected to an internal 
junction. 
 
26 
Comp onent
Component
Component
Component
Component
Compon ent
Compon ent
P
1
P
2
P
3
Component
Component
Component
Component
Component
J
1
J
2
J
3
J
4
J
5
J
6
J
1
J
2
J
3
J
4
J
5
(a)
(b)
System Boundary
 
Figure 2-3: Open (a) and Closed (b) Energy Systems 
Assuming thermo fluid energy flow, i.e. internal, kinetic and potential energy only, an 
?energy of the flow? at each port 
p
E
 
is expressed as 
2
()
2
p
pp p
v
mh zg++  , where 
p
m   is the 
mass flow rate at a given port, and h
p
 is the enthalpy of the fluid flowing through the port 
at a velocity v, and z is the port elevation in a gravitational field g. (Bird, Steward and 
Lightfoot, 1960) 
Solution of this energy system simulation requires knowing the energy flow rates at 
each component port.  The pressure difference of the fluid between ports is the driving 
potential for the flow, and will thus determine the flow rate of fluid through the 
component. 
The pressure of each port connected to a junction is identical; therefore there are N
J
 
(number of junctions) unknown pressures.  Each port may have a unique value for 
enthalpy, mass flow, position (elevation relative to a reference) and velocity, adding an 
additional 4 N
p
 (number of ports) unknown variables.  This results in a total of 4 N
p
 + N
J
 
unknown variables which must be solved for.   
 
27 
The governing equations of conservation of mass and conservation of energy at the 
junction level provide 2N
J
 equations.  The remaining equations are provided by the 
components, which provide their port position (as a component property) and determine 
their port mass flow and velocity as well as the energy state for outlet ports based upon 
boundary conditions.  These boundary conditions are the thermodynamic state at the 
ports (internal boundary conditions) and component specific independent variables. 
As mentioned previously, a special circumstance exists when internal energy plus 
flow energy represent the only energy accounted for by a flowing fluid, i.e. kinetic 
energy and potential energies are considered negligible.  In this case, the unknown 
variables for the simulation are then N
J
 unknown pressures, as well as a unique value for 
enthalpy and mass flow at each port.  This results in a total of 2 N
p
 + N
J
 unknown 
variables which must be solved for (Lindsay 2000).   
Table 2-1 summarizes the system level equations for generalized thermo fluid energy 
system simulation for a single component fluid. 
The important differences between open systems, where mass is allowed to cross the 
system boundary, and closed systems, where mass is conserved within the system 
boundary, are the boundary conditions imposed on the system, and the effect they have 
on the system charge. 
In the open system, the system?s boundaries are specified using thermodynamic 
properties, and allow mass to flow freely through them.  A change in system boundary 
conditions changes the flow through the system and consequently, the mass within.  In a 
closed system, mass is conserved not only at junctions, but for the system as a whole 
under all conditions. 
 
28 
Consequently, for the closed system, mass must be accounted for in each of the 
components, at all operating points.  A change in system boundary conditions (via the 
components) changes the operating point of the system, but not its system charge.  For 
the closed system the junction mass flow rates are no longer independent.  A mass 
balance over each junction results in a redundant equation, so one junction mass balance 
must be removed.  The number of system variables in the simulation remains the same 
however.  A mass balance over the entire system provides the equation to close the 
system of equations. 
Property Balance 
Number of 
Equations 
Equation 
Junction Mass Balance
2
 N
J
  or  N
J
-1 
0
1
=
?
=
Nports
p
p
m   
Junction Energy Balance N
J
 
2
1
()0
2
ports
N
p
pp p
p
mh z
=
+ +=
?
 
v
g  
Mass flow through port p of  
component i 
N
p
 m
p
 = ?
1
(P
in
, h
in
, P
out
) 
Thermodynamic state of each 
exit port in component 
N
p
 Out of 
Components 
h
3
 = ?
2
(P
in
, h
in
, P
out
)  
Velocity of fluid through port 
p of  component i 
N
p
 v
p
 = ?
3
(P
in
, h
in
, P
out
) 
Position (height) of port p of  
component i 
N
p
 z
p
 = Constant 
System Mass Balance
3
 1 or 0 
1
Components
N
System i
i
Mass Mass
=
=
?
 
Table 2-1: Equations used to Model Generalized Thermo Fluid Energy Systems 
 
Table 2-2 summarizes the changes in Table 2-1 for a simplified thermo fluid energy 
system simulation, where kinetic and potential energies are neglected. 
                                                 
2 The number of these equation depends upon the system boundary conditions:
 
An open system requires N
J
 mass balances, while a 
closed system requires only N
J
-1 mass balances 
3 A system mass balance is only required for closed systems. 
 
29 
 
Property Balance 
Number of 
Equations 
Equation 
Junction Mass Balance
4
 
N
J
 (Open systems) or 
N
J
-1 (Closed Systems) 
0
1
=
?
=
Nports
p
p
m   
Junction Energy Balance N
J
 
2
1
()0
2
ports
N
p
pp p
p
mh z
=
+ +=
?
 
v
g  
System Mass Balance
5
 1 or 0 
1
Components
N
System i
i
Mass Mass
=
=
?
 
Table 2-2: Equations used to Model Generalized Thermo Fluid Energy Systems with Negligible 
Kinetic and Potential Energies 
2.6 Controlling the System 
The general energy system simulation technique detailed here will result in a solution 
of the components that make up the energy system and for the internal boundary 
conditions supplied to the components.  However, energy systems can have more than a 
single operating point.  When a component has adjustable parameters, a change in one 
may affect the operation of the system.  The changing of system parameters, including 
component parameters, to achieve results is controlling of the system.  The system can be 
controlled by altering system level parameters, or component level parameters. 
For an open system, one method of control, which does not involve manipulation of 
any component parameters, is achieved by modifying a system boundary condition.  For 
example: changing the throttle pressure of a steam turbine, changing the back pressure on 
a compressor and changing the inlet temperature on a heat exchanger.  For a closed 
system, the only non-component, system boundary condition that can be utilized to 
                                                 
4 The number of these equation depends upon the system boundary conditions:
 
An open system requires N
J
 mass balances, while a 
closed system requires only N
J
-1 mass balances 
5 A system mass balance is only required for closed systems. 
 
30 
control the system is the working fluid charge contained within the system.  Increasing or 
decreasing the system charge alters the system operating point. 
System control through component manipulation can be, but is not limited to, an 
opening of a valve orifice to change its flow characteristic, the increase of a condenser 
fan speed, to provide additional cooling of the condenser, or the bleeding of hot gas from 
the compressor discharge to the evaporator inlet to accommodate low load conditions. 
It is desirable for systems to have certain system performance variables set to 
particular values.  This is done via active control of the components independent 
variables.  This is a challenge, as components do not have access to system level 
information that is not communicated through a junction.  This can be addressed by 
giving the system a control.  A control is aware of the following: a current value for the 
objective dependent variable, the desired value of that objective variable, the component 
being controlled, the independent variable to control in the component and the value of 
the independent variable. 
An example of an effective control scheme for a vapor compression refrigeration 
system using the above described technique is that of the thermostatic expansion valve 
(TXV).  A system simulation that has a TXV and uses the value of evaporator superheat 
to solve the system cannot be modeled with the generalized technique, as there is no 
physical flow between the evaporator outlet and the TXV.  Rather, a control is added to 
the system.  For this example system, the system can be aware of several system-level 
variables which may have set-point values.  For instance evaporator superheat, condenser 
sub-cooling, compressor head pressure, or compressor suction pressure.  The control can 
evaluate these system variables after the completion of a simulation.  The control can 
 
31 
then perturb the 'control variable' (TXV orifice opening) and have the system solve again, 
and update the new value of the control variable until the objective is obtained. 
2.7 Simulation Hierarchy 
The primary objective of the dissertation is the development of a general-purpose 
simulation framework, for analysis and design of the vapor-compression system and its 
components.  The thermo fluid energy system itself is the parent, or coordinating system.  
The coordinating system manipulates values of independent variables at junctions, which 
share information with ports, and thus components, which are sub-systems. 
2.8 Simulation Framework Design Rules 
A generalized framework is developed here which allows for the construction of 
thermo fluid energy systems with character similar in nature to refrigeration systems.  
Within the context of this thesis, the framework is simplified down to a steady state 
model, with negligible kinetic and potential energy.  A set of design rules defines the 
framework, and is presented as: 
1. The system independent variables consist of the value for the fluid pressure at 
each junction, as well as the value for fluid enthalpy for all flows leaving 
junctions, flowing to components. 
2. All components are individual units, and are subject to identical inputs and 
outputs.  In this case, all components must evaluate their dependent variables 
when provided values for working fluid, fluid pressure at each port, and fluid 
enthalpy at each port where fluid flows into the device. 
 
32 
3. All components must provide identical outputs.  In this case, all components 
provide the value for fluid mass flow rate at each outlet port, and the value for 
fluid enthalpy at each outlet port. 
4. The system solution must satisfy an energy balance at each junction and mass 
balance at N-1 junctions for systems with no open boundaries and N junctions 
for systems with open boundaries, where N is the number of junctions in the 
system. 
5. Components have complete control of their independent variables, other than 
port pressures and enthalpies.  Note: A component can, upon request provide a 
list of independent variables, which may be of use to a program implementing 
the framework, and serve as a system control.  This is demonstrated in chapter 
6, where VapCyc is explained. 
2.9 Summary 
This chapter details the creation of a thermo fluid simulation framework.  This 
framework is such that each component within the system is an independent, stand alone 
program.   
The framework is created upon a system of data structures and property balances.  
The structures consist of components, ports and junctions which are then combined to 
form specific energy systems. 
The components are similar to real world devices in that when boundary conditions 
are applied, a steady state solution results.  The component has several ports, to which it 
communicates the flows through the component and the thermodynamic state of fluids 
flowing through them. 
 
33 
A coordinating system, which contains sub-systems, is resolved through conservation 
equations of mass and energy.  The coordinating system has independent variables, which 
are intensive thermo-fluid properties: junction pressures and enthalpy.  The sub-systems 
then provide values for their dependent variables: mass flow and exit specific enthalpy 
for a given set of junction properties.  The coordinating system is resolved when the 
junction pressures and enthalpies provided to the sub-system result in mass flow and exit 
specific enthalpy values which satisfy the coordinating system conservation laws. 
2.10 References 
Rice, C.K., Wright, L. and Bansal, P.K., 1993, ?Thermodynamic Cycle Evaluation Model 
For R-22 Alternatives In Heat Pumps-Initial Results And Comparisons?.  Moneton, 
Canada: Proc. 2
nd
 Int. Conference on Heat Pumps in Cold Climates, Aug 16-17 
 
Bird, Steward and Lightfoot, 1960, ?Transport Phenomena?, John Wiley and Sons, 1960 
 
Hatsopolous, G., Keenan, J, 1965, ?Principles of General Thermodynamics?, John Wiley 
and Sons, 1965 
 
Obert, E., Gaggioli, R., 1963, ?Thermodynamics?, McGraw Hill, 1963 
 
Lindsay, D., 2000, ?Network Builder Notes, University of Maryland ? CEEE, Internal 
Document, July, 2000 
 
34 
3. The Junction Solver 
In this chapter, a general purpose thermo-fluid simulation solver is created to serve as 
the coordinating system in an object oriented multidisciplinary tool.  This general 
purpose solver, referred to in this thesis as the Junction Solver, provides specific 
functionality to VapCyc, a vapor compression refrigeration simulation tool.  The Junction 
Solver feature list is as follows: 
1. Single solver many thermo fluid energy systems 
2. Specifically applicable to:  
? Vapor compression refrigeration (closed system) 
? Pipe network (open or closed system) 
? Heat exchanger (open system) 
3. Allow user defined component connectivity. 
4. The functionality to change, or swap, component models used in a fixed-
component simulation at run-time thus replacing components without the 
need to re-compile the simulation code itself. 
5. The functionality to add or subtract components at run time, without the need 
to re-compile the simulation code itself. 
6. The ability to constrain and specify the total mass of working fluid for the 
system simulation. 
7. Allow for the selection from several different working fluids/refrigerants. 
8. Flexibility to alter component independent variables, independently from the 
main simulation program. 
 
35 
9. Arbitrary order of component evaluation, eliminating the need for a run-time 
specified order for component evaluation 
This chapter details necessary structure the Junction Solver possesses to fulfill the 
objectives given in the above feature list. 
3.1 The Junction Component Port Matrix (JCP) 
The foundation of the Junction Solver is the connectivity, defined here as the manner 
in which individual component ports are connected to the coordinating system junctions.  
Connectivity is described to the solution algorithm at run time, via a data structure called 
a Junction Component Port matrix, or, JCP.  
As an example, Figure 3-1 illustrates a basic vapor compression refrigeration system, 
with a hot gas bypass valve installed between the compressor discharge and evaporator 
inlet.  This valve takes the compressor discharge flow and meters it into the evaporator 
inlet to keep the evaporator pressure up during low load situations.  In the figure, 
components are numbered 1-5, junctions are numbered 1-4 and component ports are 
numbered from 1-n, with n being the maximum number of ports in an individual 
component. 
 Table 3-1 represents the JCP for the thermo fluid system depicted in Figure 3-1, a 
two dimensional matrix.  In the JCP, rows are numbered 1-N
J
, the columns are numbered 
1-N
C
, where Nj is the number of junctions in the coordinating system and N
C
 is the 
number of components (sub-systems) in the system.  The matrix elements, JCP[j][c] are 
the port number corresponding to component C which connects to junction J.   
 
36 
Q
low
Q
High
W
Net
Component
1
P
1
P
2
P
1
P
2
P
2
P
2
P
1
P
1
J
2
J
1
J
3
J
4
Component
2
Component
4
Component
3
Component
5
P
1
P
2
 
Figure 3-1: A Modified Vapor Compression Refrigeration System 
The JCP is a flexible data structure which serves as the starting point for the Junction 
Solver solution.  The JCP enables a simulation program utilizing the Junction Solver as 
its main solution engine, providing the utility given in 1-4 of the above feature list given 
in the introduction to this chapter. 
Component 
Junction 
1 2 3 4 5 
1 1 0 0 2 0 
2 2 1 0 0 1 
3 0 2 1 0 0 
4 0 0 2 1 2 
 Table 3-1: JCP for a Vapor Compression Refrigeration System 
 
37 
3.2 Component Specification 
Features 5-7 from the above feature given at the beginning of this chapter are 
achieved through the object oriented multidisciplinary nature of the simulation.  
Components, regardless of function, are to a degree ?all the same? to the coordinating 
system, in this case the Junction Solver.  Consequently, there exists a list of properties 
and methods the components must possess.  Table 3-2 displays a list of properties the 
each component utilized by the Junction Solver must posses to resolve the simulation. 
Property Description 
Charge 
The mass contained within the component 
at the current operating point. 
Port Enthalpy(p) 
The enthalpy of the working fluid, at port 
p, at the current operating point 
Port Massflow(p) 
The mass flow of the working fluid, at 
port p, at the current operating point 
Port Pressure(p) 
The pressure of the working fluid, at port 
p, at the current operating point 
Refrigerant The working fluid 
SystemUnits 
The unit system which the component 
must comply with when receiving/giving 
numeric values 
Table 3-2: A List of Necessary Properties for all Components Utilized in the Junction Solver 
 
After the Junction Solver sets the values for the component thermodynamic boundary 
conditions, the component is then run via the Run method, and values for the component 
dependent variables can be assembled and used to evaluate the Junction Solver residual 
equations. 
3.3 Independent Variables and Residual Equations 
The Junction Solver, configured by an interface, formulates a system of equations to 
resolve, subject to a set of independent variables. 
 
38 
Chapter 2 of this thesis developed the necessary property balances, namely mass and 
energy, which when properly resolved at each junction; result in a solution sufficient to 
the energy system as a whole.  This is predicated by the component models themselves, 
solving the respective mass, momentum, and energy balances.  This thesis neglects the 
effects of kinetic and potential energies. 
Recalling from chapter 2 that a component, with a number of ports numbered 1 to N
p
, 
and a number of inlet ports numbered 1 to N
p,in
, has independent variables, P
i
 and h
j
, 
where i numbers from 1 to N
p
, and j numbers from 1 to N
p,in
.  The component then 
calculates its dependent variables ,,
ii
mz 
i
v  and h
k
, where k refers to ports where mass 
flow exits, and numbers from 1 to N
p
- N
p,in
. 
The system simulation independent variables are then: P
i
, and h
j
, where, 1 < i < 
N
Junctions
, and 1 < j < N
Junction Exit Streams
. 
The Junctions Solver then uses these dependent variables to calculate the following 
residual equations: 
 
1
Nports
ip
p
rm
=
=
?
      (3.1) 
 
1
ports
N
jpp
p
rmh
=
=
?
      (3.2) 
 
1
NComponents
lSystem C
C
r Mass Charge
=
=?
?
   (3.3) 
The residual equations given here are the junction mass balance, junction energy 
balance and the overall system mass balance. 
When the junction solver operates in ?open system mode?, a mass balance is 
necessary for each junction.  When the junction solver operates in ?closed system mode?, 
 
39 
a mass balance is only necessary N
Junctions
-1 junctions, and equation (3.3), the system 
mass balance, is necessary. 
In all cases, an energy balance is required at each junction. 
3.4 Solution of the Model 
Organizing the residual equations, given in equations (3.1)-(3.3) in a vector?
 G
, the 
vector of residual values r
 G
 can be evaluated through the use of the state (solution) vector 
S
 G
.   
The resulting system of equations,  
() () 0
() ( )
aS
bSr
? =
? =
 G G
 G G
 G
     (3.4) 
is a non-linear set.  Equation (3.4)(a) is the actual set of equations, and (3.4)(b) is the 
approximated set, where r
 G
 represents the residual, or deviation from zero.  In this case, 
?
i
 represent the mass and energy balances to be solved.  S
i
 are the individual state 
variables which make up the independent variables, namely P
j 
and
 
h
j
.   
Solution of the system of equations given by equations (3.1) - (3.3) can be 
accomplished in many ways.  Two methods are currently implemented in the junction 
solver, namely (1) a simultaneous approach, where the full set of equation is specified to 
a non-linear equation solver, and (2) a modified fractional step method, based upon the 
work of Jiang (2003).  
3.4.1 Simultaneous Approach 
A general approach which does not require any a priori knowledge of system 
configuration defines the state vector S as S
T
 = [P
1
,P
..
,P
nj
, h
1
,h
?
, h
Nh
], where N
J
 is the 
 
40 
number of system junctions, and N
h
 is the number of independent enthalpy values. Here S 
is a vector of thermodynamic properties only.  The state vector also has no implications 
about system configuration.  This approach gives a ?Junction Solver?, which takes a 
candidate state vector, and populates the system junctions with pressure and enthalpy 
values.  After propagating these values of pressure and enthalpy to the proper 
components, the components are then executed and the resulting balance equations 
evaluated.  Figure 3-2 illustrates the calculation procedure for this method. 
This simultaneous approach allows employing any non-linear equation solver, but has 
the drawbacks of reliance on numerical derivatives and a need for a ?sufficient? initial 
guess. 
The current version of the Junction Solver achieves the simultaneous solution of the 
non-linear equation set via the Newton-Raphson or Broyden?s method. (Nash and Sofer, 
1996) 
3.4.2 Fractional Step 
Jiang offers a so called fractional step technique for solving heat exchangers modeled 
in a component/junction manner, as illustrated in Figure 3-3  This two step technique is 
characterized by (1) a solution for the system junction pressures through resolution of the 
system mass balances at fixed system (thermal) energy, and (2) a solution for the 
unknown junction enthalpy values through resolution of the system energy balances at 
fixed (port) flow rate. 
Although the system pressure field must always be solved simultaneously, Jiang 
exploits the fact that the energy equations can often be solved successively, reducing the 
 
41 
need to compute numeric partial derivatives.  Figure 3-3 gives a flow chart for Jiang?s 
method of solving heat exchangers. 
Seed Candidate State Vector
S[1] = P
S[2... ] = 
S[N ]  = 
S
*
1
Junct io ns
P
P
[N ] = h
S[ ] = h
S[2N ] = h
i..
NJunctions
Junct io ns1
i...
Junctions NJunctions
*
*
*
*
*
N
N..N
Junctions-1
+1
Junct io ns+ 2 Junc tions-1
Conver ged?
De-code S, and populate junctions with P  and h.
j j
Feed P and h to Components
Run  Comp onents
Evaluate Residual Equations
Upd a te S tate Ve ctor
Done
No
Yes
 
Figure 3-2: Flow Chart for Simultaneous Equations 
Jiang utilizes his fractional step technique to solve heat exchanger systems, 
interacting with a secondary fluid, namely air.  Furthermore, it is seen in Figure 3-3 that 
the solution of the component models is not completely separated from the solution of the 
refrigerant side; which is evident from the updating of the air side enthalpy values. 
 
42 
Recalling (Balling  and Sobieszczanksi-Sobieski, 1994), the system is considered 
?non-hierarchic? when coordinating and sub-system interactions are not forbidden.  This 
represents a more complex case than the ?hierarchic? one, and is beyond the scope of the 
Junction Solver. 
Additional methods are therefore necessary for the Junction Solver to utilize a 
fractional step technique, namely a (1) the ability to determine the component mass flow 
rate at the given internal (thermodynamic) boundary conditions without changing the 
component?s energy state, and (2) the ability to determine the component energy state 
given specified mass flow rate values at the ports. 
 
Figure 3-3: Fractional Step Method of Jiang (2003) 
It is noted here, updating the ?enthalpy field? for a constant component mass flow 
rate values, may result in a density change within one or more components.  This density 
change may result in a deviation between the component?s current specified mass flow 
 
43 
rate and one calculated with the current enthalpy field.  It is therefore necessary for 
components to have additional properties, namely values for their (deviations in) mass, 
and energy balances, at a specified set of port pressure and enthalpy values. 
Figure 3-4 shows a modified fractional step approach which can be used to in the 
hierarchic structure of the Junction Solver. 
Guess Initial State Vector
Component Mass Balance Converged?
Put Guess Vector Into Indi vidual Components
Solve for Junction Pr essure Values
Subject To
Junction Mass Flo w Residuals
At Constant E nergy
Done
Solve for Junction Enthalpy Values
With Implicit Component Energy Solu tion
Subject To
Junction Energy Flow Residuals
Junction Energy Balances Converged?
Component Energy Balance Converged?
Yes
Ye s
No
No
Elim inates
premature convergence
 
Figure 3-4: A Modified Fractional Step for the Junction Solver 
The modified fractional step method will iterate until the junction mass, and energy 
balances are satisfied, and the component mass, momentum, energy balances are satisfied 
too. 
 
44 
3.4.3 Solving Hydraulic Relations 
When the Junction Solver resolves the junction mass balances only, as given in 
equation (3.1) a non-linear equation solution technique such as Newton-Raphson or 
Broyden?s method is employed. (Nash and Sofer, 1996).  
3.4.4 Solving Thermal Relations 
For the special case of the Junction Solver resolving the junction energy balances 
only, as given in equation (3.2), when kinetic and potential energies are ignored, a 
successive substitution method is employed.  This technique consists of the repetition of 
two steps, namely (1) mixing of junctions, where the junction?s mixture enthalpy is 
determined and fed to all tubes that receive flow from the junction and (2) running of the 
tubes energy routines at constant mass flow.  The routine is presented graphically in 
Figure 3-5. 
Has Junction Mixture
Enthalpy C hanged?
Mix Junctions and Distribute to Componen ts
Done
Yes
No
Run Components @ Constant Mass flow
 
Figure 3-5: Successive Substitution Energy Routine 
3.5 Conditioning of the Model 
The iterative solution of nonlinear systems of equation often requires the solution of 
auxiliary linear systems; hence the latter has basic influence on the efficiency of the 
 
45 
overall method (Farago, I. and Karatson, J, 2002).  The residual equations presented in 
(3.1) - (3.3) are nonlinear in nature due to the thermo physical properties of the fluids. 
The Junction Solver developed in this thesis incorporates iterative solution techniques 
which require the solution of auxiliary linear systems.  Rather than investigate elaborate 
preconditioning techniques, the junction solver incorporates the scaling of the 
independent variables and residual equations.  The scaling is basic, scaling all 
independent variables between a maximum and a minimum value, and all residual 
equations by a characteristic value for mass flow rate or energy. 
This leaves the determination of the scaling values to the specific energy system 
simulation which utilizes the junction solver, to determine the value for scaling 
parameters. 
The scaling process is described here, while the specific methods for determining the 
value of the scaling parameters is left for future chapters, where the Junction Solver is 
utilized to simulate specific systems. 
As the independent variables of the Junction Solver are pressure and enthalpy values, 
the dimensionless independent variables are defined as: 
 
* ilow
i
high low
PP
P
P P
?
=
?
     (3.5) 
 
* ilow
i
high low
hh
h
hh
?
=
?
     (3.6) 
The Junction Solver must be provided values for P
high
, P
low, 
h
high
, h
low
, which ideally 
represent the highest and lowest values for pressure and enthalpy encountered in the 
system. 
 
46 
The residual equations given in (3.1) - (3.3) represent the mass and energy balances 
of the system.  The Junction Solver resolves scaled version of these equations, give as: 
1
*
Nports
p
p
i
m
r
m
=
=
?
 
 
     (3.7) 
 
1
*
ports
N
p p
p
j
mh
r
E
=
=
?
 
 
     (3.8) 
1
*
NComponents
System C
C
l
Mass Charge
r
Mass
=
?
=
?
    (3.9) 
The values of 
**
,mE
 
   and Mass
*
, represent characteristic values for the system mass 
flow, energy flow and total working fluid mass. 
3.6 The Solution Seed 
The Junction Solver attempts to resolve a state vector, S
 G
, subject to equations (3.7) ? 
(3.9).  The junction solver itself is completely independent of system configuration, thus 
the determination of the seed is given to the specific energy system simulation which 
utilizes the junction solver. 
3.7 Summary 
This chapter describes the general function of the Junction Solver, a general purpose 
thermo fluid simulation solver.  In addition to the JCP matrix, which describes the system 
connectivity to the Junction Solver, a limited set of component properties and methods is 
introduced.  The Junction Solver solution methodologies are developed and presented 
here as well.  The solution methodologies are dependent upon the component level 
properties and methods presented in this chapter. 
 
47 
Lastly, a general method for scaling the problem for numeric stability is presented.  
The method is broad, introducing scaling factors which must be provided values by the 
thermo fluid system simulation utilizing the Junction Solver. 
The simulation framework presented is able to serve many of the diverse needs of the 
constituents, as presented previously, namely: 
1. Single solver for many thermo fluid energy systems. 
2. Allow user defined component connectivity. 
3. The functionality to change, or swap, component models used in a fixed-
component simulation at run-time thus replacing components without the 
need to re-compile the simulation code itself. 
4. The functionality to add or subtract components at run time, without the need 
to re-compile the simulation code itself. 
5. The ability to constrain and specify the mass of working fluid for the 
simulation of system simulation.  For a vapor-compression refrigeration 
application, this is referred to as ?Charge Management?. 
6. Allow for the selection from several different working fluids/refrigerants. 
7. Flexibility to alter component independent variables, independently from the 
main simulation program. 
8. Arbitrary order of component evaluation, eliminating the need for a run-time 
specified order for component evaluation. 
3.8 References 
Balling, R.J. and Sobieszczanksi-Sobieski, J., ?Optimization Of Coupled Systems: A 
Critical Overview Of Approaches?, Proceedings of the 5
th
 AIAA/NASA/USAF/ISSMO 
 
48 
Symposium on Multidisciplinary Analysis and Optimization, Panama City, Florida, 
September 7-9, 1994.  Paper No. AIAA-94-4330-CP 
 
Farago, I. and Karatson, J., 2002, Numerical Solution of Nonlinear Elliptic Problems Via 
Preconditioning Operators, Nova Science Publishers, Inc, New York, 2002 
 
Jiang, H., 2003, Development Of A Simulation And Optimization Tool For Heat 
Exchanger Design, Ph.D. Thesis, Department of Mechanical Engineering, University of 
Maryland, 2003 
 
49 
4. Simulation of a Heat Exchanger 
A secondary objective of the dissertation is the development of a heat exchanger 
model built upon the Junction Solver.  A heat exchanger of the type simulated in this 
chapter, represents a compound thermo fluid energy system, where two systems interact, 
specifically a refrigerant system and an air system.  Many heat exchanger models exist, 
and the objective of this simulation is the illustration of a specific incarnation of the 
Junction Solver.  The example demonstrates the Junction Solver, and accomplishes three 
goals: 
1. Validation of the Junction Solver general purpose thermo fluid simulation 
tool 
2. Utilization of the Junction Solver to solve a compound system where two 
dependent thermo fluid systems interact 
3. Verification of a simulation model which will be used in a system simulation 
model, VapCyc, and facilitate its validation 
The actual heat exchanger models are designed to represent physical pieces of 
hardware which have been instrumented and tested in the laboratory (Hwang et al. 2004).  
4.1 Model Components 
4.1.1 Tubes 
Tubes are finned conduits for fluid flow.  As are all other thermo fluid energy 
systems, they are characterized by mass, momentum, and energy conservation.  For the 
scope of this heat exchanger study, kinetic and potential energies are considered 
 
50 
negligible, and the steady state balance equations necessary to simulate the tube are: 
(Bird et al., 1960) 
Continuity: 
0
in out
mm=?
??
      (4.1) 
For the pipe geometry, the continuity equation simply sates the mass flow into the 
pipe is equivalent to the mass flow out of the pipe. 
Momentum: 
( )0
tot
mp m=?? + ?vg  SF+    (4.2) 
The momentum equation is the relationship between the flow, and the pressure 
difference required to drive the fluid through the resistance the pipe walls exert upon the 
fluid.  The unknown quantity in the momentum equation is this force F.  Typically, 
correlations created from empirical data are used to estimate this force. 
Energy: 
0
p
umQW
?
????
=?? + + ?
????
????
     (4.3) 
The energy equation governs the change between fluid internal energy, m  u, flow 
work, m 
p
?
, heat transfer out of the control volume, Q, and work done by the fluid, W.  
For a pipe, there is no (shaft) work done on the working fluid. 
The tube has three distinct modes of operation.   
1. When supplied appropriate internal boundary conditions of pressures and 
enthalpy values for all inlet streams, and external boundary conditions (air 
temperature, flow rate), the tube will resolve its mass flow rate and exit 
 
51 
enthalpy state which meets both the thermal and hydraulic relations for the 
model. 
2. When supplied an inlet pressure and outlet pressure, a mass flow rate, inlet 
enthalpy value, and external boundary conditions, the tube will resolve an 
exit enthalpy which meets the thermal relations.  Note: The resulting solution 
may violate hydraulic relations for the model. 
3. When supplied with boundary pressures and enthalpy values, the tube will 
resolve a mass flow rate which meets the hydraulic relations.  Note: The 
resulting solution may violate the thermal relations for the model. 
4.1.1.1 Hydraulic Relations 
White (1986), under the assumptions of incompressible flow, gives a method for 
obtaining the pressure drop through a pipe as 
 
351
44 4
7/4
0.158P LdV??
?
??     (4.) 
For this calculation, density, ? , viscosity, ?, tube diameter d and fluid flow velocity, 
V are representative of average values.  This formula can be rearranged to an explicit 
form for mass flow 
5
4
31
44
2.8702
dP
mA
L
?
??
??
?
?
??
??
      (4.5) 
It is noted here, that this correlation is valid in the range 4000 < Re
d
 < 10
5
.  While the 
simulation may have values for Re
d
 outside of this range, the relationship captures the 
characteristic behavior of the flow for most applications, and serves the purposes of this 
work. 
 
52 
Additionally, a correction term, DP
Factor
 is included in the model.  This term is a 
constant that can be adjusted to fine tune the tube performance. 
5
4
31
44
Factor
DP  
2.8702
P
d
mA
L
?
??
???
??
?
??
??
      (4.6) 
It is stressed here that this model is dependent upon average values for the fluid 
thermophysical properties. 
For a pipe, with a specified mass flow rate and thermodynamic state, the pipe model 
can then calculate a mass flow rate based upon this thermodynamic state.  If the 
calculated mass flow rate does not agree with a specified mass flow rate, a mass flow 
residual, representing the tube?s violation of its mass balance is given by 
mass specified calculated
rm m=?      (4.7) 
4.1.1.2 Thermal Relations 
Relations between the heat transfer and the heat exchanger geometry can be 
calculated using an effectiveness-NTU method (White, 1984). 
 
()
1, 1,
2, 1,
1
1
oi
ii Max
p
TT
Q
TT Q
UA UA
NTU
c mc
?
?
==
?
==
 
 
 
     (4.8) 
The effectiveness (? ) is a ratio of heat transfer to maximum possible heat transfer, 
while the number of transfer units (NTU) is a ?thermal size? relating heat transfer area 
(A) and heat transfer coefficient (U) to heat capacity rate (
p
mc  ).  The heat transfer kinetic 
relations for the two fluids are then used to complete the effectiveness-NTU relationship. 
 
53 
 
1
11
Contact Fouling
oi
U
RR
hhFR
=
+++
    (4.9) 
Here the overall heat transfer coefficient is described in terms of the outer (air) and 
inner (refrigerant) heat transfer coefficients (h) as well as thermal resistances for 
conduction, surface fouling and an outside to inside surface area ratio (FR). 
The tubes are modeled as air to refrigerant heat exchanger in cross-flow, with the air 
side mixed.  The effectiveness-NTU relationship for this is (White, 1984) 
unmixedmixed,)]);exp(1[exp(1(
1
unmixedmixed,];
)exp(1
exp[1
minmax
*
*
maxmin
*
*
CCNTUC
C
CC
C
CNTU
????=
??
??=
?
?
 (4.10) 
The heat capacity rate of a fluid 
p
Cm   is denoted by C, and C* is the ratio of two heat 
capacity rates, such that C* ? 1.   
Although this relationship is fine for super-heated and sub-cooled liquids, this 
relationship does not apply to fluids in a two-phase state.  A condensing or evaporating 
fluid can be thought to have an infinite heat capacity since its temperature does not 
change with an increase in energy.  Therefore when a single phase fluid is in thermal 
contact with an evaporating or condensing fluid, the effectiveness is expressed as: 
 )exp(1 NTU??=?     (4.11) 
4.1.1.3 Refrigerant Side Heat Transfer Coefficients 
For simplicity, all flows are modeled as turbulent.  Single phase heat transfer 
coefficients are calculated using Gnielinski?s correlation (Kakac, 2002, p.96), given as 
 
54 
 
()
2
1/2
1.58ln( ) 3.28
(/2)( )
112.7( /2) ( 1)
/
d 2/3
d
fRe
f Re -1000 Pr
Nu
fPr
hNukd
?
=?
=
+?
=
   (4.12) 
In this case, the correlation is valid for fully developed turbulent flow with a constant 
wall heat flux, Pr > 0.7.   
Condensation heat transfer coefficients are calculated using Shah?s (1989) 
relationship 
 ()
()
0.04
0.76
0.8
0.38
3.6 1
1
lo
Cr
xx
hh x
P
P
??
??
?
=?+
??
??
????
   (4.13) 
where h
lo
 is the heat transfer coefficient were the fluid only liquid, and in this case is 
calculated by Dittus and Boelter (1930) as 
 ( )
0.8 0.3
0.023 /hkdRePr=    (4.14) 
The above correlation is valid in the range Re
lo
 > 350, Re
vo
 > 35000. 
4.1.1.4 Air Side Heat Transfer Coefficients 
Likewise Kim et. al (1999) provide heat transfer calculations in terms the fin 
geometry (height, depth, spacing), maximum air velocity and air thermophysical 
properties. 
-0.369 0.106 0.0138 2/3
max
0.163 ( / ) ( / ) Pr
air d Fin Fin Fin air air
h Height Depth Spacing d u Cp?
?
=
Re
 (4.15) 
For the geometry and air flow given in Figure 4-4, this results in h
air
 = 68.99 W/m
2
-K. 
 
55 
4.1.1.5 Air Surface Area Calculations 
The air side calculations are independent of the refrigerant side of the tube.  A 
specified tube geometry and air flow rate facilitates the calculation of the air heat transfer 
coefficient, fin ratio and fin efficiency.   
Utilizing relations for fin ratio and fin efficiency given in Chang et al. (1997)  
()
Tube Fins Fins
A d L N Thickness?=?    (4.16) 
The fin surface area is the cross-sectional area of the fin parallel to the air flow. 
,
2( )
Fins Fins Fin Fin cross tube
A N Depth Height A=?   (4.17) 
The fin ratio is defined as the ratio of fin surface area to tube surface area were no 
fins present. 
Tube Fins
A A
FR
dL?
+
=      (4.18) 
4.1.1.6 Distributed Models 
The effectiveness-NTU relationship provided by equation (4.8) is only adequate for 
constant values of the overall heat transfer coefficient U.  In situations where a tube has 
both heat transfer between air and a single phase fluid, as well as heat transfer between 
air and a condensing fluid, the value of U will change radically along the length of the 
tube.  In this thesis, a distributed heat transfer model, where each heat transfer mode 
(desuperheating, condensing and subcooling) is accounted for separately is implemented 
to help account for the radical change in thermophysical properties. 
 
56 
4.1.2 A Refrigerant Circuit 
For the purposes of heat exchanger modeling, a refrigerant circuit is defined here as 
the flow path of refrigerant from an entrance to an exit.  A refrigerant circuit may branch 
off into two or more paths, re-combine and/or exit at one or more locations.  The 
refrigerant circuit is simply all tubes in hydraulic contact with one another between inlet 
and outlet.  Figure 4-1 illustrates a possible complex refrigerant circuit. 
 
Figure 4-1: A Refrigerant Circuit 
4.1.3 An Air Circuit 
The tubes modeled earlier in this are prescribed a certain air velocity, or flow rate.  In 
the tube models, air properties are attributes of the tube itself.  When given the 
appropriate air boundary conditions, inlet temperature and flow rate, the outlet air 
temperature is calculated.  Realistically, the air flow is an energy system distinct from the 
refrigerant.  Although in thermal contact with the refrigerant circuit(s), the air flows 
along its own path, and is therefore hydraulically independent from a tube. 
For simplicity in this thesis, the equations governing the hydraulic flow of the air are 
ignored, and a flow rate is imposed upon the air circuit.  Consequently, continuity is 
implicit.  Figure 4-2 illustrates a single air circuit flowing over three refrigerant segments.   
In this particular case, the air is flowing in the +y direction, over three separate tube 
segments, where the refrigerant flows along the x-axis.  Note there is no assumption that 
the tube segments be part of the same refrigerant flow stream. 
 
57 
 
Figure 4-2: An Air Stream: Air Flows from Left to Right, Through Three Air Segments, Simulating 
Air Flow Over Three Tube Segments 
Each air circuit has one or more tube segments over which it flows.  The three 
individual ?cells? are air segments, combined to make the single air stream.   
Using the energy equation given in Table 2.1, an energy balance, the governing 
equation of energy for the air stream is: 
,,
i,Ref
1i1
()0
air out air in
NSegments NSegments
i
i
mh h Q
QQ Q
==
?+=
==?
??
 
 
    
    (4.19) 
In the case where moist air has water vapor condense on the tube segment, 
conservation of species (water vapor/liquid water) is no longer implicit, and species 
balances must also be included.  The mass and energy equations for moist air condensing 
on heat exchanger segments are (McQuiston and Parker, 1988): 
Comment: In figure caption, explain in 
more detail the goemetry.
 
58 
2
22
,
,,,,
i,Ref
1i1
()
0
air out in H O out
air air out H O out air in H O in fg
NSegments NSegments
i
i
mm
mh h h h QQ
QQ Q
? ?
??
==
?=
+??++=
==?
??
  
  
 
    
  (4.20) 
Recalling that junctions have no interaction with the fluid, the junction energy 
balance is: 
22
2
,,,,
,,
()0
Or, by defining: 
0
air air out H O out air in H O in
air out H O out
out in
mh h h h
hh h
hh
? ?
?
+??=
=+
?=
 
   (4.21) 
4.2 The Actual Heat Exchanger 
The actual heat exchanger modeled is a cross flow tube-fin heat condenser.  Three 
distinct circuits enter the heat exchanger, and combine into a single stream, which serves 
as a sub-cooler.   
Figure 4-3 illustrates the configuration of the heat exchanger modeled in this chapter, 
while the physical parameters for the heat exchanger are given in Figure 4-4. 
 
Figure 4-3: Cross Section of Heat Exchanger Configuration.  Air Flows Across Banks of Tubes from 
Left to Right 
 
59 
 
Figure 4-4: Heat Exchanger Physical Parameters 
The condenser given in (Hwang, et. al 2004) was part of a system, and thus subject to 
various system test points.  Table 4-1 details the operating point of the condenser under 
the various system test points for the R-404a case.  Similarly, Table 4-2 details the 
operating point of the condenser for the various system test points for the R-290 case. 
Test Parameter Test1 Test2 Test3 Test4 
Pin (kPa) 2422.5 2515.1 2709.6 1743.3 
Tin (K) 353.0 355.2 360.8 337.0 
Pout (kPa) 2285.5 2384.9 2589.3 1606.9 
Subcooling (K) 8.2 13.5 17.8 12.0 
Air Inlet 
Temperature (K) 
308.7 308.3 308.3 291.5 
Mass flow rate (g/s) 100.0 99.9 98.5 100.0 
Q (kW) 16.04 16.72 17.02 17.65 
Table 4-1: Test Data for R-404a Condenser 
 
 
 
 
 
60 
Test Parameter Test1 Test2 Test3 Test4 Test5 
Pin (kPa) 1764.9 1786.1 1820.8 1877.3 1238.9 
Tin (K) 353.3 354.2 355.1 357.7 331.2 
Pout (kPa) 1667.0 1695.4 1735.5 1795.2 1130.8 
Subcooling (K) 8.1 10.2 12.5 15.4 6.1 
Air Inlet 
Temperature (K) 
308.0 307.9 307.9 307.8 292.0 
Mass flow rate 
(g/s) 
41.7 41.6 41.7 41.2 42.7 
Q (kW) 15.89 16.1 16.35 16.45 16.77 
Table 4-2: Test Data for R-290 Condenser 
4.3 The Heat Exchanger Model 
A heat exchanger is constructed by combining refrigerant circuits, as defined in 
section 4.1.2, with air streams, as defined in section 4.1.3.  The two energy systems are 
then coupled by associating the air properties of individual air segments with the air 
properties of the appropriate refrigerant segment.  
The heat exchanger represented by Figure 4-3 and Figure 4-4 is modeled here in with 
a slight geometric modification.  For a  given refrigerant circuit path, when the circuitry 
snakes back upon itself, while staying in the same plane, perpendicular to the air flow, the 
tubes in that plane are combined into a single tube, thus combining multiple ?passes? into 
a single ?pass?.  This process is demonstrated in Figure 4-5, where the total number of 
tubes is reduced from 66 to 15.  As a consequence, a ?pseudo length factor? (PLF) is 
introduced.  This factor multiplies a given tube?s surface area and length for the purposes 
of calculating heat transfer and pressure drop for individual tubes. 
 
61 
 
Figure 4-5: Geometry Modification For Heat Exchanger Model.  The Real Condenser on the Left is 
Approximated with Circuitry on Right 
Some factors which may contribute to the inaccuracy of the model to represent the 
actual heat exchanger include: 
? The reduction of many tube passes down to a single tube pass 
? The use of a single segment to represent a relatively long length of tube (0.99 
m) with lumped pressure drop and heat transfer correlations 
? The possibility of the flow characteristics being beyond the acceptable range 
for the given correlations 
As the scope of this work is more the use of the heat exchanger to validate a 
simulation framework, rather than an accurate heat exchanger model, three additional 
correction factors (CF) are introduced to the model, for the purposes agreement between 
experimental data and the model.  Equations (4.22) - (4.25) show how the pressure drop 
correction factor (CF
PD
), liquid heat transfer correction factor (CF
LHT
) and two phase 
heat transfer correction factor (CF
TPHT
) are used in the modeling equations. 
351
44 4
7/4
PD
0.158 CF  P LPLF d V??
?
??    (4.22) 
 
62 
LHT
CF /
Liquid d
hNuk=     (4.23) 
T
CF /
Two Phase PHT d
hNuk=    (4.24) 
()1
1
p
UA UA PLF
NTU
c mc
==
 
    (4.25) 
For the purpose of modeling the heat exchanger, the values for the correction factors 
are determined using a design test point (test 1) from Table 4-1 and Table 4-2.  The 
values are determined by solving the following equations: 
 0
calc test
mm? =      (4.26) 
 
,,
0
out calc out test
hh? =     (4.27) 
 
,subcooler inlet calc Sat
hh=     (4.28) 
Equations (4.26) - (4.27) represent an agreement between experimental and simulated 
mass flow rate and outlet enthalpy (thus heat transfer), while equation (4.28) attempts to 
meet a design criteria that the final three tubes in the refrigerant circuit be used as a sub-
cooler circuit. 
Results for the correction factors with the condenser using R-404a are given in Table 
4-3, while results for the condenser using R-290 are given in Table 4-4. 
CF
PD
 CF
LHT
 CF
TPHT
 
4.25 0.7 0.9 
Table 4-3: Correction Factors for Modeled Heat Exchanger (R-404a) 
 
CF
PD
 CF
LHT
 CF
TPHT
 
6.0 0.2 2.5 
Table 4-4: Correction Factors for Modeled Heat Exchanger (R-290) 
 
63 
The relatively high value of the correction factors indicate that the modeling 
assumptions and correlations used are not the best for simulating this particular heat 
exchanger.  They do however allow this model to accurately simulate this heat exchanger 
at the design point, and allow for the model to represent the actual heat exchanger in the 
vicinity of the design operating point.  
4.3.1 Wrapping the Junction Solver 
The Junction Solver is a non-hierarch system, meaning the sub-systems; in this case, 
tubes are not allowed to communicate.  Consequently, the coordinating system, in this 
case, the heat exchanger itself, must facilitate communication between sub-systems, 
(components) such that convergence criteria are met.   
Configure heat exchanger
Set air side temperature and enthalpy values
Are air junction energy
balances satisfied?
Propagate air properties from
outle t of t ubes t o inlet of
corre sponding  next tube.
This implicitly satisfies air energy balance
Done
Solve Junction Solver, such that:
Refrigerant Junction Energy Balances
Component Energy B alance Converged
Yes
No
 
Figure 4-6: Heat Exchanger Solution Flow Chart 
Figure 4-6 illustrates the process of wrapping the Junction Solver within a routine 
which updates the heat exchanger air side.  The Junction Solver is responsible for 
 
64 
resolving the refrigerant side, while the coordinating heat exchanger is responsible for 
providing air property values (the tube external boundary conditions) to the Junction 
Solver sub-systems. 
4.3.2 Seeding the Junction Solver for the Heat Exchanger 
The Junction Solver can be filled with any potential solution.  When no a priori guess 
value is known, the junction solver initializes to a linear pressure distribution and a 
constant thermal distribution.   
The algorithm for determining the linear junction pressure distribution is given in 
Figure 4-7. 
Has Average Junction
Pressure C hanged?
Average Port Pressures at Each Junction
Set Each Port Pressure to This Pressure
Done
Yes
No
Average Tube Pressure and
Apply to Each Port
 
Figure 4-7: Algorithm for a Linear Pressure Distribution 
The value for the junction enthalpies is set to the heat exchanger inlet enthalpy.  The 
air temperature is set to the heat exchanger inlet air temperature at all air junctions. 
4.3.3 Scaling Parameters for the Model 
Recall, that for conditioning purposes, the Junction Solver has the parameters built in 
for scaling: P
high
, P
low, 
h
high
, h
low
, 
**
,mE
 
  .  Recall the junction solver scales independent 
variables in the form of 
* ilow
i
high low
PP
P
P P
?
=
?
     (4.29) 
 
65 
* ilow
i
high low
hh
h
hh
?
=
?
     (4.30) 
Similarly, the Junction Solver balance equations are scaled as 
1
*
Nports
p
p
i
m
r
m
=
=
?
 
 
     (4.31) 
 
2
1
*
()
2
ports
N
p
p pp
p
j
mh z
r
E
=
++
=
?
 
 
v
g
    (4.32) 
For heat exchanger simulation, the values of P
high
 and P
low
 come from the heat 
exchanger boundary pressures, specifically the inlet and outlet pressure respectively.  For 
a condenser, the refrigerant inlet enthalpy is the highest value, and represents h
high
, while 
h
low
 is calculated as the minimum enthalpy the refrigerant can achieve, that which 
corresponds to the air inlet temperature at the heat exchanger outlet pressure.  As the heat 
exchanger models are designed after a specific heat exchanger, values for mass flow and 
heat of rejection at the design point are used for
**
,mE
 
  .  Taking the Test 1 condition from 
Table 4-1 the independent variables are scales as 
*
2285.5
2422.6 2285.5
ilow i
i
high low
PP P
P
PP
??
==    (4.33) 
,
*
( 2422.6, 352.95) 419.830 /
( 2285.5, 308.65) 250.741 /
250.741
419.830 250.741
high in in
low out air in
ilow i
i
high low
hhP T kJkg
hhP T kJkg
hh h
h
hh
== = =
== = =
??
==
  (4.34) 
While the residual equations are scaled as 
1
0.100 /
Nports
p
p
i
m
r
kg s
=
=
?
 
    (4.35) 
 
66 
 
1
16.04
ports
N
p p
p
j
mh
r
kW
=
=
?
 
    (4.36) 
Recall from Figure 4-6 that a fractional step solution is employed.  In this particular 
case, a simultaneous solution of the hydraulic equations is utilized, while a successive 
substitution routine is employed for the thermal step.  Conditioning of the equations for 
use in iterative linear algebra routines is therefore only necessary for the hydraulic step.   
Figure 4-8 illustrates the effect the scaling of the Junction Solver independent 
variables and residual equation has on the solution process.  Each curve in Figure 4-8 
represents an individual hydraulic step of the algorithm shown in Figure 4-6, which itself 
has one or more iterations to solve. 
Iteration Number
Ja
co
b
i
a
n
C
o
n
d
i
ti
o
n
N
u
m
b
e
r
(
P
r
e
ssu
r
e
R
e
si
d
u
a
l
k
g
/
s
-
P
a
)
1 2 3 4 5 6 7 8
0
20
40
60
80
100
Iteration #12
Iteration #20
Iteration #19
Iteration #1
(a)
Iteration Number
J
a
c
o
b
i
a
n
C
o
n
d
i
t
i
o
n
N
umbe
r
(
D
i
m
e
ns
i
o
nl
e
s
s
P
r
e
s
s
u
r
e
R
e
s
i
dua
l
)
1 2 3 4 5 6 7 8
0
20
40
60
80
100
Iteration #20
Iteration #19
Iteration #1
(b)
 
Figure 4-8: Condition Number of Jacobian Matrix versus Hydraulic (Inner) Iteration Number.  The 
Iterative Hydraulic Routine is Wrapped in a Successively Substituted Outer Thermal Routine.  
Scaling of the Independent Variables and Residual Equations Reduces the Condition Number of the 
Jacobian Matrix of the Hydraulic Routine.   
Figure 4-8 shows curves representing the condition number of the Jacobian matrix 
used for solving the pressure calculations in the fractional step technique.  Each curve 
represents the value of the Jacobian condition number at individual iterations during the 
solution of the pressure routine. 
Comment: This needs more detail.  
What does ?Inner? mean?  Or refr to text 
in fgure caption. 
 
67 
Figure 4-8(a) has curves showing the condition number of the Jacobian when 
unscaled independent variables and residual equations are used, while the curves shown 
in Figure 4-8(b) has curves representing the Jacobian when independent variables and 
residual equations are scaled in the manner given in equations (4.33) and (4.35). 
It is seen that due to the nature of the equations, most of the pressure calculations are 
converged to tolerance within two iterations. 
It is evident from the plots that performance improvements are seen.  For the scaled 
variables and equations, the initial iteration has a reduced condition number, and 
completed in one less iteration.  Improvements in later iterations are seen, as after the 
first iteration, no pressure routine for the scaled variables and equations requires more 
than four iterations, where some of the iterations for the unscaled variables and equations 
do. 
4.4 Validation 
For verification, the heat exchanger model is run for the R-404a ?test 1? case, as 
given in Table 4-1.  Graphical output is used to verify that the model is indeed behaving 
in a manner similar to that of a real heat exchanger. 
For test purposes, convergence tolerances are set to 0.01 W and 10
-8
 kg/s for the heat 
exchanger simulation.  Plots of the heat exchanger thermophysical properties are made to 
verify that the model is solved properly, when convergence occurs. 
4.4.1 Thermophysical Property Distributions 
Figure 4-9 displays the refrigerant pressure distribution and quality distribution.  The 
refrigerant pressure and quality behave in a predictable manner.  Figure 4-9(a) illustrates 
 
68 
a gradual decrease in refrigerant pressure along the flow path(s), while Figure 4-9(b) 
shows a steady drop in the refrigerant quality, along the condenser flow path(s). 
Figure 4-10 displays the temperature distribution for the heat exchanger.  The 
temperature distribution of the heat exchanger is as expected.  Figure 4-10(a) shows the 
refrigerant cooling off along the length of the heat exchanger path(s).  Due to the 
condensation, the temperature remains near constant, while it condenses.  The figure is 
consistent with Figure 4-9(b), in that the nearly constant temperature section corresponds 
to the sections of the heat exchanger where the refrigerant quality is greater than 0, but 
less than 1.  Figure 4-10(b) (note the heat exchanger is viewed from the opposite side) 
illustrates the air temperature distribution superimposed upon the refrigerant tubes.  In 
this figure, air is flowing from left to right, and the temperature rises along the air flow 
path. 
 
69 
 
Figure 4-9: Refrigerant (a) Pressure and (b) Quality Distribution.  Air Flow, Not Shown In the 
Figure, Flows from Left To Right Over the Tube Banks, Cooling Them. 
4.4.2 Mass Balances and Energy Balances 
Figure 4-11 illustrates the refrigerant junction and tube mass and energy balances.  
Figure 4-11(a) illustrates that at convergence, the junction and tube mass balances are 
less than the convergence criterion 10
-8
 kg/s (in absolute value).  Similarly, Figure 
4-11(a) illustrates that the refrigerant tubes and junction energy balances are less than the 
specified 0.01 W. 
 
Comment: Better show airflow too
 
70 
 
Figure 4-10: Heat Exchanger Temperature Distribution ? 
 (a) Refrigerant Only and (b) Refrigerant and Air 
 
Figure 4-11: Heat Exchanger (a) Mass Balance and (b) Energy Balance 
 
71 
4.5 Verification 
The performance results predicted by the simulation are compared against the 
experimental data provided and displayed in Table 4-1, Table 4-2, and Figure 4-12 - 
Figure 4-17.  Table 4-1 gives a numerical comparison of the R-404a condenser, while 
Table 4-2 gives a numerical comparison of the R-290 condenser.  See Hwang et al 
(2004), for a detail on the actual experiment.  
The experimental data represents the single heat exchanger geometry under a variety 
of test conditions, with two working fluids.  Given the specified boundary conditions on 
the condenser, namely inlet and outlet pressures, refrigerant inlet state, air inlet 
temperature and flow rate, the relevant performance properties, which correspond to 
model outputs are: refrigerant mass flow rate, heat rejection and refrigerant exit 
thermodynamic state represented by the degree of subcooling. 
Figure 4-12 graphically displays the comparison of simulated to measured mass flow 
rate for the R-404a condenser for the four test points.  The first three test points, which 
correspond to the same inlet air temperature, but different refrigerant pressures and inlet 
energy state have good agreement (within 5%) to the test data. 
The final test point, which corresponds to a test point with a lower inlet air 
temperature, shows a deviation of approximately 6%. 
Figure 4-13 shows that under all test points, the heat of rejection reported by the 
simulation is within 5%. 
Figure 4-14 illustrates agreement within 5% for all test points except low inlet air 
temperature test point, which deviates by approximately 11%.  
 
72 
Figure 4-15 graphically displays the comparison of simulated to measured mass flow 
rate for the R-404a condenser for the four test points.  The first three test points, which 
correspond to the same inlet air temperature, but different refrigerant pressures and inlet 
energy state have good agreement (within 5%) to the test data. The final test point, which 
corresponds to a test point with a lower inlet air temperature, shows a deviation of 
approximately 10%. 
Figure 4-16 shows that under all test points, the heat of rejection reported by the 
simulation is within 5%. 
Figure 4-17 illustrates subcooling agreement within 5% for all test points except low 
inlet air temperature test point, which deviates by approximately 11%.  
 
Figure 4-12: Mass Flow Rate Comparison for R-404a Condenser (See Table 4-1) 
 
73 
 
Figure 4-13: Heat of Rejection Comparison for R-404a Condenser (See Table 4-1) 
 
 
Figure 4-14: Subcooling Comparison for R-404a Condenser (See Table 4-1) 
 
74 
 
Figure 4-15: Mass Flow Rate Comparison for R-290 Condenser (See Table 4-2) 
 
 
Figure 4-16: Heat of Rejection Comparison for R-290 Condenser (See Table 4-2) 
 
75 
 
Figure 4-17: Subcooling Comparison for R-290 Condenser (See Table 4-2) 
4.6 Summary 
In this chapter, the secondary objective of the dissertation, the development of a heat 
exchanger model built upon the Junction Solver is accomplished.  The model facilitates 
(a) validation of the Junction Solver as a general purpose thermo fluid simulation tool, 
and (b) a heat exchanger model for use in VapCyc to facilitate its validation. 
In this chapter the necessary simulation models of tube and air segment are combined 
via the Junction Solver creating a condenser, which is suitable for use in VapCyc.  The 
actual heat exchanger models are designed to represent physical pieces of hardware 
which have been instrumented and tested in the laboratory (Hwang, et al. 2004).  
The tube models allow for three correction factors (pressure drop correction factor 
(CF
PD
), liquid heat transfer correction factor (CF
LHT
) and two phase heat transfer 
correction factor (CF
TPHT
)), which are determined numerically to match a specific test 
case for the heat exchanger, for each of the two tested fluids. 
 
76 
Due to the factors, the simulation data produced from the simulation results matches 
the experimental data to within a few percent for flow rate and heat of rejection. 
This chapter achieves the goal of validating the Junction Solver as a general purpose 
thermo fluid simulation tool by its success in the application of modeling a specific 
thermo fluid system, namely an air cooled condenser.  In addition, the R-134a and R-290 
condensers are general purpose components, suitable for use in VapCyc, another specific 
application of the Junction Solver. 
4.7 References 
Chang, Y., Wang, C., 1997, ?A Generalized Heat Transfer Correlation for Louver Fin 
Geometry?, International Journal of Heat and Mass Transfer, Vol.40, pp. 533- 544. 
 
Dittus, S.J., Boelter, L.M.K., 1930, ?Heat Transfer In Automobile Radiators Of The 
Tubular Type?, University of California Publications on Engineering, Vol.2, pp.443-461. 
 
Hwang, Y., Dea-Huan, J., Radermacher, R., ?Comparison Of Hydrocarbon R-290 And 
Two Hfc Blends R-404a And R-410a For Medium Temperature Refrigeration 
Applications?, ARI: Global Refrigerant Environmental Evaluation Network (GREEN) 
Program, Final Report, March 2004 
 
Kakac, S., Liu, H., 2002, Heat Exchangers: Selection, Rating, And Thermal Design, 
2th edition, CRC Press. 
 
Kim, N.H., Young, B., Webb, R.L., 1999, ?Air-Side Heat Transfer and Friction  
Correlations for Plain Fin-and-Tube Heat Exchangers with Staggered Tube 
Arrangements?, Journal of Heat Transfer, Transactions of ASME, Vol.121, pp.662-667. 
  
Nash, G. S., and Sofer, A., 1996.  Linear and Nonlinear Programming, McGraw ?Hill 
publishing 
 
Shah, M.M., 1989, ?A General Correlation for Heat Transfer during Film Condensation 
inside Pipes?, International Journal of Heat and Mass Transfer, Vol.22, pp.547-556. 
 
White, F. M., Fluid Mechanics, Second Edition, McGraw-Hill, 1986 
 
 
77 
5. VapCyc ? Vapor Compression Refrigeration System 
To this point, this dissertation has focused on the creation of a generalized thermo 
fluid energy system simulation tool, the Junction Solver.  The power of this simulation 
tool is realized when sub-systems of different character are combined into a larger 
system.  
One objective of this section is the development of a general-purpose comprehensive 
vapor-compression refrigeration simulation tool, VapCyc, for analysis and design of the 
vapor-compression system and its components.  Building VapCyc upon the Junction 
Solver allows for the following objectives to be met: 
? The ability to constrain and specify the mass of working fluid for the system 
simulation.  For a vapor-compression refrigeration application, this is referred 
to as ?Charge Management?. 
? The functionality to change, or swap, component models used in a fixed-
component simulation at run-time thus replacing components without the need 
to re-compile the simulation code itself. 
? The functionality to add or subtract components at run time, without the need 
to re-compile the simulation code itself. 
? A (limited) user defined component connectivity. 
? The selection from several different working fluids/refrigerants. 
? Flexibility to alter component independent variables, independently from the 
main simulation program. 
? Incorporation of  both ?interactive? and ?automated? modes of operation 
 
78 
? Off design simulation 
This chapter presents the VapCyc tool, as a specific implementation of the Junction 
Solver.  Examples are used to validate and verify against data from a research project in 
the CEEE heat pump laboratory (Hwang et al. 2004), as well as demonstrate the 
overcoming of the numerical issues of conditioning and providing a necessary initial 
guess.  This method of providing an initial guess to the system solution is provided via 
the introduction of a novel concept called the ?context?.  
5.1 VapCyc Component Specifications 
Section 3.2 detailed component specification for use within the Junction Solver 
framework.  VapCyc is a simulation which incorporates the Junction Solver.  While the 
components used within VapCyc must have all the specifications of a component used 
within the Junction Solver, VapCyc requires specifications beyond those of the basic 
Junction Solver specification.  For example, since VapCyc is a vapor compression 
refrigeration simulation, it must report the system?s coefficient of performance (COP).  
This is calculated by dividing the system?s cooling capacity by the power consumption.  
Thus all VapCyc components must provide both their cooling capacity (typically only the 
evaporator would have a non-zero value) and their power consumption.   
The object oriented manner in which components are designed allows application 
programs, such as VapCyc to incorporate component properties which are non-traditional 
in standard simulations.  For example, a component property, ?Lead Time?, gives an idea 
of how quickly a system could be built.  Utilizing the properties of ?gravity? 
(gravitational constant) and component ?mass? allows the study of systems in different 
applications, such as earth and in space. 
 
79 
Table 5-1 illustrates the component properties that all VapCyc components are 
required to provide, for this introductory version of VapCyc.  The manner in which 
VapCyc is written allows additional properties to be added later. 
Property Description 
Charge 
The mass contained within the component at 
the current operating point. 
Cost 
The cost to purchase/manufacture the 
component 
Gravity 
The value of the gravitational potential 
experienced by the component 
HeatOut 
The energy rejected to the environment via 
heat transfer at the current operating point 
Lead Time 
Time (in weeks) required to 
obtain/manufacture the component 
Manufacturability 
A numeric value representing the ability/cost 
to manufacture the component 
Mass 
The mass of the component (devoid of 
working fluid) 
Port Enthalpy(p) 
The enthalpy of the working fluid, at port p, at 
the current operating point 
Port Massflow(p) 
The mass flow of the working fluid, at port p, 
at the current operating point 
Port Pressure(p) 
The pressure of the working fluid, at port p, at 
the current operating point 
Port Velocity(p) 
The velocity of the working fluid, at port p, at 
the current operating point 
Port Position(p) 
The position in a gravitational field of the 
working fluid, at port p, at the current 
operating point 
Power Consumption 
Power consumed by the component at the 
current operating point 
Refrigerant The working fluid 
Reliability 
A number representing the reliability of the 
component at the current operating point 
SystemUnits 
The unit system which the component must 
comply with when receiving/giving numeric 
values 
Volume Volume of the component 
WorkOut 
The energy rejected to the environment via 
work (rate) at the current operating point 
Table 5-1: VapCyc Component Properties 
 
80 
5.2 The Test System 
Due to growing environmental awareness and resulting concerns, refrigerants, the 
working fluids for refrigeration systems, heat pumps and air conditioners, have attracted 
considerable attention. Following policies to reduce global warming, industry is 
developing technologies that can reduce emissions and improve energy efficiency 
(Hwang, et al. 2004).  The scope of the project is to evaluate the performance of a 
medium temperature refrigeration unit with respect to the working fluid, R-404a, R-290 
and R-410a. 
This test project is a useful vehicle for the validation and verification of VapCyc.  
Rather than completely reproduce the results of the test project in its entirety, only the R-
404a and R-290 test results are reproduced. 
The tests performed for the test project include a ?charge optimization? which 
includes the design point of the systems, and a ?part load? test.  The charge optimization 
consists of a fixed system configuration with environmental boundary conditions (air 
conditions at both the condenser and evaporator) are held constant while the refrigerant 
charge is varied to maximize the system COP.  The part load test consists of the fixed 
system configuration at the ?optimum? charge, evaluated at a second set of air conditions. 
Appendix I contains the data points taken from the test project, as well as the process 
for the creation of the TEST components.  The TEST components are used in VapCyc to 
verify the simulation against the test data. 
 
81 
5.3 Simulation of Vapor Compression Refrigeration  
Vapor compression (v-c) refrigeration represents a large part of the United States 
energy consumption, and therefore a practical target for energy savings.   
Experimental analysis of v-c refrigeration is difficult, as even a simple system can 
have many independent variables; making complete control of the system during 
experimentation difficult if not impossible.  Logically, simulation is a good vehicle to 
analyze system performance over a rigidly controlled set of independent system 
variables.  Simulation also offers an opportunity to save time and effort when 
optimization of a system is the goal. 
The study of vapor compression refrigeration has been facilitated by many computer 
simulations created over the years, dating back to approximately 1976.  These 
simulations are created for the purposes of design and analysis, but rarely optimization of 
the system.  The large majority of the simulations created are very specific to a single 
system configuration.  Typically, to make a very specific simulation more flexible a large 
number of independent variables are allowed in the simulation.   
The literature demonstrates moderate efforts to utilize optimization in vapor 
compression refrigeration at the system level.  Previous optimization work focuses on (a) 
simulation programs designed for and integrated with optimization routines and (b) 
optimization results obtained from manually retrieving data from simulation programs.  
In previous work, programs designed for optimization (Rice, et al. 1981, Levins et al. 
1996) are limited in nature due to the complex nature of the optimization/simulation 
integration.  Manual optimization techniques (Fischer and Rice, 1986) are (i) tedious to 
perform and (ii) not conducive to expansion.  Thus, an opportunity exists to create a 
 
82 
simulation structured such that a great deal of the functionality of the optimization 
routine is determined at run-time, thus separating the simulation and optimization 
routines, adding flexibility and utility to the optimization process.   
5.4 Premise  
VapCyc is a tool for a steady-state vapor compression refrigeration system 
simulation.  The simulation is designed to be comprehensive and robust, as the goal is its 
utilization in optimization routines.  
VapCyc is built around four independent component models, namely compressor, 
condenser, expansion device and evaporator.  VapCyc simulates the performance of 
systems comprised of at least this basic set.  The simulation consists of a set of 
independent variables, such as: system charge, component model numbers and 
component independent variables, and results in a set of dependent system variables, such 
as: COP, capacity, weight and volume.  Using the general system simulation modeling 
technique outlined in chapter 2, VapCyc also facilitates the insertion of additional 
components, beyond the four basic components into the system.  The additional 
components can be positioned in parallel; between any two junctions, or in series; 
between any two components. 
Although the basic configuration of the system is fixed, the models and independent 
variables, the boundary conditions, for each component can be variable.  Once the 
individual component models, the component boundary conditions and a value for the 
refrigerant charge are chosen, the complete set of independent variables for the system is 
specified.  The simulation is achieved through the solution of the system level 
conservation laws of mass, and energy, as described in chapter two.  The evaluation of 
 
83 
the components subject to their boundary conditions, allows for the solution of the system 
level conservation laws. 
Component inter-changeability, coupled with a solver for the run-time determined set 
of independent variables is achieved through a two level object oriented scheme.  The 
system level object is responsible for resolving the system conservation laws.  The 
component level objects are responsible for the reporting of their mass flow, refrigerant 
charge and energy interactions with their environment to the system level when given 
boundary conditions. 
5.5 VapCyc?s Component Library 
5.5.1 A Generic Compressor Model 
A simple model for a compressor accounts for all compressor workings with two 
compressor properties; namely an isentropic efficiency and a volumetric efficiency.  
When subject to thermodynamic boundary conditions, the volumetric displacement per 
revolution and compressor speed, the operating point of the compressor is determined.  In 
this context, the operating point refers to the flow rate through the compressor (both mass 
and volumetric flow rates), as well as the thermodynamic state of the refrigerant at the 
compressor discharge port. 
The model inputs are,  
? Volumetric Efficiency,  
RPMV
m
m
m
actualactual
v
?
?
 
 
 
==
max
 
? Isentropic Efficiency,  
Actual
Isentropic
s
W
W
 
 
=?  
? RPM (Compressor speed in revolutions per minute) 
 
84 
? Displacement Volume (per revolution, V) 
The results of the model, which solves for mass flow rate and outlet enthalpy are: 
 ()
in v
s in
out in
s
mVRPM
hh
hh
??
?
=
?
=+
 
    (5.1) 
The model is now ready for a general energy system simulation as described in 
chapter 2.  With the energy system providing thermodynamic boundary conditions to the 
model, the compressor is able to provide a mass flow at each of its ports, and a value for 
the enthalpy at each of its outlet ports. 
It is stressed here that this generic compressor model is a specific compressor model 
designed to execute quickly, to facilitate the validation of VapCyc.  Although this model 
is fairly simply in nature, models of any level of sophistication can be used within 
VapCyc. 
5.5.2 The Alpha Compressor 
VapCyc is equipped with a library of components, which have a range of nominal 
capacities at the VapCyc air-conditioning context.  The components are created in 
families, where a family of a certain component has the same operating characteristics as 
all others in the family.  For example, using the techniques of section 5.5.1 a family of 
compressors called the Alpha Compressor was created.  The Alpha Compressor is based 
upon the Generic Compressor model of section 5.5.1.  The model has four independent 
variables, namely displacement per revolution, speed (rev/min), volumetric efficiency 
(?
v
) and isentropic efficiency (?
s
).  The Alpha Compressor Family fixes RPM at 1000 
rev/min, ?
v
 at 0.95 and ?
s
 at 0.65.  Based upon R-134a performance, a family of 
 
85 
compressors is sized for displacement/revolution to give nominal capacities of 1.24, 2.44, 
3.6, 5.0, 10.0, 11.0 and 20 kW. 
Table 5-2 shows the Alpha Compressor Family, rated with different refrigerants, and 
Table 5-3 shows the Alpha Compressor Family, rated for low temperature refrigeration 
applications. 
 
Compressor Model Compressor Capacity
6
 (kW) 
 R-134a R-22 R-12 R-404a Ammonia 
AlphaCompressor1240 1.24 1.93 1.20 1.95 2.17 
AlphaCompressor2440 2.44 3.81 2.37 3.84 4.28 
AlphaCompressor3600 3.60 5.63 3.50 5.67 6.32 
AlphaCompressor5000 5.00 7.81 4.86 7.88 8.78 
AlphaCompressor10000 10.0 15.6 9.71 15.7 17.5 
AlphaCompressor11000 11.0 17.7 10.7 17.3 19.3 Co
mp
res
s
o
r
 Nam
e
 
AlphaCompressor20000 20.0 31.2 19.4 31.5 35.1 
Table 5-2:Alpha Compressors Rated For Air-Conditioning 
 
Compressor Model Compressor Capacity
7
 (kW) 
 R-134a R-22 R-12 R-404a Ammonia 
AlphaCompressor1240 0.177 0.359 0.210 0.334 0.336 
AlphaCompressor2440 0.349 0.707 0.414 0.658 0.662 
AlphaCompressor3600 0.515 1.04 0.611 0.970 0.977 
AlphaCompressor5000 0.716 1.50 0.849 1.35 1.36 
AlphaCompressor10000 1.43 2.90 1.70 2.70 2.71 
AlphaCompressor11000 1.57 3.18 1.86 2.96 2.98 
Co
m
p
re
ss
or 
Na
m
e
 
AlphaCompressor20000 2.86 5.80 3.40 5.40 5.43 
Table 5-3: Alpha Compressors Rated For Low-Temperature Refrigeration 
 
                                                 
6 Capacity based upon 40?F (4.4?C) / 115?F (46.1?C) evaporating / condensing temperature, 7?F superheat 
/ 10?F (3.9-5.6?C) subcooling 
7 Capacity based upon -40?F (-40.0?C) / 105?F (40.6?C) evaporating / condensing temperature, 7?F 
superheat / 10?F (3.9-5.6?C) subcooling 
 
86 
5.5.3 The Beta and Gamma Compressor 
The Beta and Gamma compressors are essentially the same as the Alpha 
Compressors.  The difference between models is the isentropic efficiency value.  Beta 
compressor?s have an isentropic efficiency of 0.85, while the Gamma Compressor has an 
isentropic efficiency of 0.55.  The three similar compressors are provided to offer high, 
medium and low efficiency compressors of similar capacity. 
5.5.4 The Delta Compressor 
The Delta compressor is built upon the Alpha Compressors.  To make the model 
more realistic, the isentropic efficiency of the compressor varies as a function of the 
compression ratio the compressor pumps against.  The isentropic efficiency versus 
compression ratio is dependent upon the compressor?s working fluid, and therefore 
currently the Delta compressor only supports simulation for R-134a applications. 
 
Figure 5-1: ?
s
 versus Compression Ratio for Delta Compressor 
 
87 
Figure 5-1 shows the value of the Delta compressor?s isentropic efficiency as a 
function of pressure ratio.  The Delta compressor is designed to have an isentropic 
efficiency identical to that of the Alpha compressor at the air conditioning design point, 
?
s
 = 0.65 @ r
p
 = 3.4. 
5.5.5 ARI 10 COEFFICIENT MODEL 
The compressor model is based upon the 10-coefficient modeling method given in the 
Air Conditioning and Refrigeration Institute standard ARI 540-99 (ARI, 1999).  For this 
model, a compressor?s performance is described by polynomials of the form 
223223
01 2 3 4 5 6 7 8 9
XCCSCDCSCSDCDCSCDSCSDCD=+ + + + + + + + + . (5.2) 
In this equation, X, represents a performance value, C
i
 are coefficients, specific to the 
performance value of interest, S is the Suction Temperature (in ?F) and D is the 
Discharge Temperature (in ?F).  For this particular model, two sets of coefficients are 
used to represent the compressor mass flow rate and power requirement at a given 
operating point.  Note, suction and discharge temperature are the dew temperature 
corresponding to the compressor?s suction pressure and bubble temperature 
corresponding to the compressor?s discharge pressure respectively. 
Compressor manufacturer?s supply the coefficients for their specific compressors by 
testing the compressor under the conditions given in the standard and fitting the 
coefficients to the test data. 
5.5.6 The Generic Orifice 
Stoecker (1958) gives a relation for the flow of refrigerant through an orifice as 
2
mKD p?=?       (5.3) 
 
88 
where K is a dimensionless, p?  has units of Pa, ?  has units of kg/m
3
 and m  has units of 
kg/s.  The generic orifice uses this formula to compute the mass flow rate through the 
orifice in terms of the fluid thermodynamic state and the specified pressure boundary 
conditions on the orifice. 
5.5.7 Generic Heat Exchanger 
A generic heat exchanger model is one which can offer a great deal of flexibility.  
One type of model is a single tube, with refrigerant flowing on the inside, and air flowing 
in a cross-flow configuration on the outside.  If the model were to experience either 
condensation or evaporation, it divides into sub-segments, keeping heat transfer and 
momentum relations separate for the different heat transfer and momentum transport 
modes. 
The single tube model is not realistic in practice.  To keep a generic heat exchanger 
model as simple as possible the length of the heat exchanger is an equivalent length for a 
realistic heat exchanger with the given duty of that modeled heat exchanger. 
The following model assumptions which apply to the condensing and evaporation 
sections, as well as single phase flow can be made for a computationally efficient heat 
exchanger model. 
? Constant air flow rate and thermo-physical properties 
? Heat transfer coefficients are functions of fin ratio and fin efficiency 
? Refrigerant mass flow is a function of component boundary pressures 
? All of the (refrigerant) pressure drop for the component occurs in the 
thermodynamic regime (vapor, two phase or liquid) of the entering fluid, 
and over the entire length of the heat exchanger 
 
89 
 
Using a classic relation of pressure drop to flow is (Bird, et al., 1960),  
?
?
?
?
?
?
?
?
?
=
2
2
1
4
1
V
P
L
D
f
?
.    (5.4) 
This relation is valid for laminar flow only, which leads to inaccuracies under 
conditions which result in turbulent flow.  It does however qualitatively capture the 
behavior of the pressure ? resistance characteristic.  
An explicit form of mass flow can be determined, which is an advantage 
computationally.  
fL
P
Dm
24
2
5
?? ?
=  .    (5.5) 
Taking the heat exchanger inlet properties, and using the entire heat exchanger length 
as the parameter L, allows for a simple approximation of the heat exchanger mass flow 
rate as a function of pressure difference, inlet density (enthalpy) and heat exchanger 
geometry. 
This radical simplification may differ greatly from the actual behavior one may 
expect in the laboratory, but it captures qualitatively the behavior looked for from the 
model, namely a decrease in flow with decreased pressure drops and/or increased inlet 
superheats.  The main benefits of the simplification are the reduced calculation time and 
elimination of the need for iteration in the heat exchanger model. 
It is stressed here that this generic heat exchanger model is a specific compressor 
model designed to execute quickly, to facilitate the validation of VapCyc.  Although this 
model is fairly simple in nature, models of any level of sophistication can be used within 
VapCyc. 
 
90 
Heat transfer between air and refrigerant during the phase change is modeled with a 
simple relation, 
( )
()
?
?
?
?
?
?
?
? +
?=
?=
2
,,
,ref
,,
outrefinref
Inair
InrefOutrefrefref
TT
TUAQ
TTCpmQ  
    (5.6) 
Combining the two expressions in (5.6), and rearranging terms, allows the refrigerant 
outlet temperature to be expressed as, 
UACm
TTUATCm
T
p
inrefinAirinrefp
outref
+
?+
=
 
 
2
)2(2
,,,
,
  (5.7) 
Heat transfer between air and refrigerant (either evaporating or condensing) is 
modeled through a simple relation, 
( )
)(
,,Ref
/,
inrefoutrefref
CondEvapInairAir
hhmQ
TTUAQ
?=
?=
 
    (5.8) 
A caveat of this approach is a potential violation of the 2
nd
 Law of Thermodynamics.  
This is guarded against by checking the outlet temperature of the refrigerant.  If it is less 
than the saturation temperature, than the model is sub divided into sub-sections. 
5.5.8 The Alpha Heat Exchangers 
The Alpha heat exchangers are based upon the same model as the generic heat 
exchangers described in 5.5.7.  In this case values for length and diameter of the heat 
exchanger tube are held constant, thus simulating a particular heat exchanger. 
Using the VapCyc air-conditioning context, the Alpha Condenser and evaporator are 
sized by determining the appropriate length and diameter, which allows them to match a 
desired capacity at the nominal rating point. 
 
91 
5.5.9 The Eta and Epsilon Condensers 
Although the Alpha Condenser family calculates quickly and qualitatively behaves 
well, it is not sensitive to air conditions. In reality, air conditions may affect the heat 
exchanger performance and therefore the system as a whole.  The Eta and Epsilon 
condenser are constructed from distributed segments, and sensitive to air conditions. 
ASHRAE (1988) recommends air face velocities around 6.5 ft/s (2 m/s) and air flow 
rates in the range of 600 to 1200 CFM /ton (80 to 160 l/s-kW), resulting in air 
temperature changes of 15 to 40?F (7 to 22?C). 
The Eta Condenser family was designed to meet not only the refrigerant design 
specifications of the VapCyc air-conditioning context, but also the ASHRAE 
recommendations for the air side.  As was revealed in chapter five, heat exchangers built 
from many segmented tubes can be very costly computationally.  Therefore, a single 
segment which meets the above requirements was modeled.   
Although the single segment may not represent a realistic model, it represents a 
scaled version of the desired condenser.  It is therefore easy to imagine several of these 
segments in parallel.  To get a desired capacity, the segment values need only be 
multiplied by a performance factor.  Table 7-4 shows the dependent variables provided 
by the Eta Condenser at the VapCyc air-conditioning context. 
Condenser Property Nominal Value
Mass flow (g/s) 0.152 
Heat Transfer (W) 29.04 
Charge (g) 0.108 
Inlet Air Temperature (?C) 35.0 
Outlet Air Temperature (?C) 42.0 
Table 5-4: Eta Condenser Nominal Operating Point for Air-Conditioning 
 
 
92 
Similar in nature to the Eta condenser family is the Epsilon condenser family.  Rather 
than a single segment representing the heat exchanger, a single tube represents the heat 
exchanger.  This tube can be split into an arbitrary number of segments, allowing the heat 
transfer and pressure drop to be evaluated with greater accuracy, with respect to the 
physical phenomena.   
The Eta and Epsilon condensers also calculate air pressure drop.  According to Shah 
in (Kreith,1999), the pressure drop through the core of an array of tubes covered with fins 
can be approximated as: 
2
41
2
ch m
fLG
p
gD ?
??
?=
??
??
   
 (5.9) 
using the heat exchanger friction factor, f , the fin depth, L, the mass flux, defined as 
m
G
A
=
 
, the hydraulic diameter, D
h
 and the mean air density.  The pressure drop is 
calculated from equation (5.9), assuming a constant density and a friction factor given by 
(Webb, 1994) as: 
 ( )
1.318
0.521
0
0.508 /
dt
fReXd
?
=     (5.10) 
In this equation, X
t
, is the fin depth and d
0
 is the tube diameter. 
The fan power requirement is given as the thermodynamic power (sometimes referred 
to as Brake Horse Power) required to increase the given flow through the calculated 
pressure drop.  The actual power requirement is the thermodynamic power divided by a 
fan efficiency and motor efficiency (White, 1986). 
 
fan motor
gQH
P
?
??
=      (5.11) 
 
93 
In this case the air volume flow rate, Q is multiplied with air density, ? , the 
gravitational acceleration constant and the head, defined as 
p
H
?
?
= .  The Eta and 
Epsilon condensers use ?
fan
 = ?
motor
 = 1.0; 
5.5.10 Heat Exchanger Summary 
A summary the heat exchangers presented in this chapter is detailed in, Table 5-5 
Name Features 
Generic Heat Exchanger 
? Constant air flow rate and thermo-physical 
properties 
? Heat transfer coefficients are functions of 
fin ratio and fin efficiency 
? Refrigerant mass flow is a function of 
component boundary pressures 
? All of the (refrigerant) pressure drop for the 
component occurs in the thermodynamic 
regime (vapor, two phase or liquid) of the 
entering fluid, and over the entire length of 
the heat exchanger 
Alpha 
Condenser/Evaporator 
? Based upon Alpha Heat Exchanger 
? Fixed Geometry 
Eta Condenser 
? Single distributed segment (see section 4.2) 
? Constant Geometry 
? Air pressure drop/power requirement 
Epsilon Condenser 
? A series of distributed segments  
? Constant Geometry 
? Air pressure drop/power requirement 
Table 5-5: A Summary of VapCyc Heat Exchangers 
5.6 VapCyc Utility 
As an engineering tool, VapCyc provides a certain amount of utility, as outlined next.  
The tool is designed to offer a comprehensive set of features, which are currently lacking 
in existing software.  The ability to add functionality is naturally facilitated via the 
generalized energy system modeling approach detailed in this thesis. 
 
94 
5.6.1 Charge Management 
Recalling the modeling equations of chapter 2, as presented in Table 2-2; systems are 
categorized into two categories: open and closed.  With regards to mass, an open system 
is one where mass is allowed to flow across the system boundaries, while a closed system 
contains a constant amount of mass at all times.  Charge management is an accounting of 
the mass within a system, and thus utilizing the closed system mass balance, the charge 
can be specified, and accounted for at all system operating points. 
5.6.2 Component Inter-Changeability 
The generalized energy system modeling technique calls for a set of components, 
which operate in concert with one another via the solution of system level mass and 
energy balances.  Thus an advanced tool, such as VapCyc, facilitates component inter-
changeability by employing this technique.  The components are distinct engineering 
models.  The resolution of the modeling equations presented in Table 2-2 is accomplished 
independently from components. 
Furthermore, by allowing the user to interface with the component models 
themselves, individual component boundary conditions can be modified independently of 
the solver resolving the system level equations of Table 2-2.   
5.6.3 Multiple Refrigerants 
Requiring the component models themselves to resolve their operating point given 
specified boundary conditions allows for a system level solver which is independent of 
choice of refrigerant.  In VapCyc each component utilizes RefProp (NIST, 2001).  
VapCyc can poll components at run time to see which refrigerants from RefProp the 
 
95 
components are capable of utilizing.  By taking all components into account, a list of 
refrigerants available to the user to utilize in the simulation is presented. 
5.6.4 Multiple Unit Systems 
Requiring component models to operate with a specified set of unit systems, allows 
the VapCyc interface to specify the units of the boundary conditions specified to 
components.  This allows the user of the simulation to operate in a unit system of choice. 
5.6.5 Physical Variables 
A vapor compression refrigeration system simulation is used to estimate many 
physical variables over a range of operating conditions.  The current list of physical 
variables calculated in the simulation includes: 
? COP 
? Capacity 
? Weight 
? Volume 
? SEER ? Seasonal Energy Efficiency Ratio 
5.6.6 Economic Variables 
A meaningful simulation also supplies economic variables of interest.  A current list 
of variables includes: 
? System cost 
? System reliability 
? System manufacturability 
 
The system level variables are calculated using values obtained from the individual 
components, which are dependent upon individual models.  The usefulness of this is 
realized when component values, such as component reliability are functions of 
 
96 
component performance.  An example is a compressor which calculates a reliability value 
which is dependent upon the compressor?s pressure ratio.   
5.6.7 Seasonal Performance 
ARI standard 210/240-94 ?Unitary Air Conditioning and Air-Source Heat Pump 
Equipment? (ARI, 1996) provides a testing standard, which results in a seasonal figure of 
merit.  This figure of merit is the Seasonal Equipment Energy Efficiency Ratio (SEER), a 
ratio of the equipment cooling or heating performance to the system?s electrical 
requirements.  The calculation procedure given in the standard is straightforward, and 
implemented in VapCyc. 
5.6.8 Interactive Mode 
The user of VapCyc is able to obtain meaningful results without a thorough 
understanding of the programming and/or mathematical modeling involved in the 
simulation.  To this end, a ?user-friendly? interface allows the user to modify the 
simulation interactively and view the results immediately. 
Figure 5-2 is a screen shot of the VapCyc Cycle Calculation page, with the Alpha 
3600 system loaded.  This page represents the main interface for interactive use of 
VapCyc.  Visible in the page are: the four icons, representing the four main VapCyc 
components, model numbers for the currently selected components, the selected 
refrigerant, the system charge, system performance variables (not all are available before 
the simulation is run), and the nominal capacity and charge of the four main components 
in the current context.  For demonstration purposes, a figure of merit called ?Final Energy 
 
97 
Balance? is given on this screen.  This represents the total energy balance of the system 
after solution (in W). 
 
Figure 5-2: The Alpha 3600 System in Air-Conditioning Mode 
The pictures representing models can be clicked with the mouse, allowing the user to 
(a) choose a component, (b) edit the current component or (c) view the component?s 
output variables (after run). 
Figure 5-3 illustrates the changing of the compressor, in this particular case; the 
Alpha 3600 compressor is replaced with the Alpha 20000 compressor. 
 
(a)       (b) 
Figure 5-3: Changing of the Compressor 
 
98 
Note the change in the system variables that accompany the compressor change.  
Upon the change of any variable, the nominal system charge and/or nominal capacity 
may change.  The component nominal capacities and charges can be seen in Figure 
5-3(b).  Here, it is seen that changing the compressor has not changed the nominal system 
charge, but the compressor?s nominal capacity is now quite different than the other three 
components. 
5.6.9 Automated Mode: Parametric Analysis 
To be a useful tool for more than a single operating point, parametric analysis is 
offered through an automated mode.  This automated mode, operating with Microsoft 
Excel, allows system configurations to be specified by text only.  Multiple runs are then 
performed without monitoring the simulation. 
5.7 VapCyc System Model 
Figure 5-4 shows a simple four component vapor compression refrigeration system.  
The components are connected to one another via junctions.  The components and their 
connectivity are specified to the Junction Solver to form the vapor compression 
refrigeration system. 
5.7.1 Scaling the Independent Variables and Balance Equations 
Recall, that for conditioning purposes, the Junction Solver has the parameters built in 
for scaling: P
high
, P
low, 
h
high
, h
low
, 
**
,mE
 
  .  Recall the junction solver scales independent 
variables in the form of 
* ilow
i
high low
PP
P
P P
?
=
?
    (5.12) 
 
99 
 
* ilow
i
high low
hh
h
hh
?
=
?
    (5.13) 
Q
low
Q
High
W
Net
Condenser
Evaporator
Expansion Device Compressor
J
P
h
1
1
1
J
P
h
2
2
2
J
P
h
3
3
3
J
P
h
4
4
4
Port 1Port 8
Port 2
Port 4
Port 5
Port 6
Port 7
Port 3
 
Figure 5-4: Simple Vapor-Compression Refrigeration System 
Similarly, the Junction Solver balance equations are scaled as 
1
*
Nports
p
p
i
m
r
m
=
=
?
 
 
    (5.14) 
 
2
1
*
()
2
ports
N
p
p pp
p
j
mh z
r
E
=
++
=
?
 
 
v
g
   (5.15) 
The choice of scaling parameters should represent some physical meaning to the 
system.  This often offers meaningful values to convergence criteria and independent 
variable bounds.  The system depicted in Figure 5-4 however, does not give any 
 
100 
indication of characteristic pressures (P
high
, P
low
), enthalpies (h
high
, h
low
), mass flow rate, 
or energy flow rate. 
5.7.1.1 The System Context 
Scaling the variables and equations of a simulation when the result of the simulation 
is unknown is achieved using the ?real life? concept of equipment rating.  Figure 5-5 
shows an example of a compressor manufacturer?s model number explanation 
(Tecumseh, 2004).  Included in the compressor model number is a code for a nominal 
rating point, and nominal capacity at that rating point. 
 
Figure 5-5: A Compressor Model Number (Tecumseh, 2004).   
The compressor model number shown in Figure 5-5  represents a compressor for 
?high temperature application? (45?F / 7.2?C) and provides a cooling capacity of 4000 
Btu/h (1.17 kW) at its rating point.  The rating point, although not published in Figure 5-5 
is defined in ARI standard 540-99 (ARI, 1999), and consists of suction and discharge 
pressures, as well as in inlet gas temperature.  The compressor?s flow rate can be 
 
101 
correlated to an evaporator capacity by multiplying it by a nominal heat transfer per unit 
mass.  This nominal heat transfer is given in terms of an evaporating pressure, 
condensing temperature condenser exit liquid temperature (subcooling) and evaporator 
superheating, or compressor return gas temperature.  The rating point and the 
compressor?s cooling capacity at that point provide enough information to derive the 
compressor?s suction and discharge pressures, as well as the suction temperature and 
mass flow rate. 
Figure 5-6 shows rating points which vary with compressor application, as published 
in the 540-99 standard. 
 
Figure 5-6: ARI 540-99 Rating Points for Compressors 
A similar set of rating points is given in the standard for compressors which are to be 
used as both cooling and heating compressors in a heat pump application.   These rating 
points are shown in Figure 5-7. 
Even though compressors have different operating characteristics, having similar 
ratings at a common operating point allows for an educated guess that the compressors 
should behave similarly, at least, in the region of this defined operating point. 
 
102 
Consider the two compressors represented by published manufacturer?s data shown in 
Figure 5-8 and Figure 5-9, the Copeland CS10K6E-PFV and ARJ11K1E-CAV 
respectively (Copeland, 2004).   
 
Figure 5-7: ARI 540-99 Heating and Cooling Compressor Operating Points 
At the ASHRAE high temperature rating point (marked in red) in Figure 5-8, these 
compressors are very similar, giving a cooling capacity of 11,600 and 11,300 Btu/hr (3.4 
and 3.3 kW) respectively.   The capacity difference at that particular rating point is only 
2.6%.  However, comparing the two displacement rates, it is seen that these are not as 
mechanically similar as the capacity at the ASHRAE high temperature rating point would 
suggest; with displacement rates differing by 7.5%.  The capacity of the machines, when 
rated at the lower temperature rating point, corresponding to the second ASHRAE rating 
point in Figure 5-9 (marked in blue) differ by approximately 10%. 
In this thesis, a distinct thermodynamic operating point, which is applied to an 
individual piece of equipment, is deemed to be a context.  Using this nomenclature, each 
of the five entries in Figure 5-6 is a unique context at which a given compressor will have 
a distinct operating point and performance. 
 
103 
A compressor can be described with reasonable accuracy simply by denoting a 
nominal operating condition.  Choosing either the high temperature (45?F / 7.2?C) or 
medium temperature (20?F / -6.7?C) as a context, a nominal operating point is 
established.  Provided the solution of the system model does not force the compressor 
model to vary from this range significantly, its mass flow rate is fairly close to the 
nominal value.  Thus, the nominal flow rate (or nominal capacity), as defined by the 
nominal rating point is a good choice for scaling the mass flow rate and mass balance 
equations. 
 
Figure 5-8: Performance Summary of a Copeland CS10K6E-PFV Compressor (Copeland Web Site) 
Similar to the compressor, the remaining components have standards for their rating.  
Evaporators are rated in such standards as ARI 420-2000, Unit Coolers for Refrigeration 
and ARI 440-1998, Room Fan-Coils.  Condensers are specified in standards, such as ARI 
450-99, Water-Cooled Refrigerant Condensers, Remote Type, and ARI 460-2000, 
Remote Mechanical-Draft Air Cooled Refrigerant Condensers.  Expansion devices are 
rated according to ARI 750-2001: Thermostatic Refrigerant Expansion Valves.  Each 
 
104 
standard defines a context in which the given piece of equipment can be, at least partially, 
characterized. 
 
Figure 5-9: Performance Summary of a Copeland ARJ11K1E-PFV Compressor  
(Copeland Web Site) 
Although these standards were chosen because they are all refrigeration/air-
conditioning components, the standards themselves are for individual components.  There 
is no insinuation that these devices will work in concert with others.  In practice, when 
assembled into a refrigeration system, the eventual operating point will likely not result in 
the specified cooling capacity of any of the individual components.  However, if the 
individual pieces of equipment have similar nominal capacities, at the same nominal 
operating point, the system will operate near this nominal operating point. 
Individual component models do not necessarily relate to one another in any intuitive 
fashion.  The challenge to the simulation program is to establish if one piece of 
equipment (model) is compatible with another, and if so, how might it go about solving 
the system model. 
 
105 
VapCyc utilizes a system context to help overcome the challenges.  From a simulation 
stand point, it is not unreasonable for VapCyc to supply components a context, or 
operating point, and request its flow rate and/or heat transfer in that context.  This offers a 
first approximation of component compatibility for the system as a whole.  
5.7.1.2 The VapCyc Air-Conditioning Context 
Fixing at least four components in the cycle, the compressor, condenser, expansion 
device and evaporator, VapCyc deviates from a general energy system to a specific one.  
VapCyc allows for the addition of components and the ability to change model and 
independent variables of components, which makes it general, but ultimately, VapCyc is 
some sort of vapor compression refrigeration cycle. 
Taking advantage of this allows a context to be specified.  This context is used to 
scale the simulation independent variables and balance equations.  One particular context 
employed by VapCyc is a standard air-conditioning system.  The compressor pressures 
are set by determining the nominal evaporating and condensing temperatures.  A typical 
air-conditioning application would be approximately 40?F (4.4?C) evaporating and 115?F 
(46.1?C) condensing.  Nominal pressure drop through the condenser and evaporator are 
assumed to be 3%.  Superheating (evaporator) and sub cooling (condenser) are typically 
in the 7-10?F (3.9-5.6?C) ranges. 
Figure 5-10 and Figure 5-11 graphically depict the system context of this simple 
vapor compression cycle.  Figure 5-10 shows the set of nominal values that make up the 
standard air-conditioning context for VapCyc, while Figure 5-11 is a Pressure-enthalpy 
diagram of the nominal cycle.  Components rated at this context provide mass flow, 
outlet energy state and refrigerant charge. 
 
106 
Q
low
Q
High
W
Net
Condenser
Evaporator
Expansion Device Compressor
P = Psa t(40?F) = Psat(4.4?C )
P = 342.954 kPa
SH = 7?F = 3.89?C
P = Psat(115?F) = Psat(46.1?C)
P = 1193.70 kPa
SH = 35?F = 20?C
Subcooling = 10?F = 5.6?C
P = 0.97 Pin
P = 1157.919 kPa
Port 1Port 8
Port 2
Por t 4
Port 5
Port  6
Port  7
Por t 3
P =1.03 Pevap
P = 353.561 kPa
 
Figure 5-10: A Nominal Air-Conditioning Operating Point for VapCyc 
Consider a context given to four components; were the condenser to have a rated 
mass flow value 20 times that of the other components, this is a sign that the system may 
be poorly matched.  Poorly matched systems in practice may function to some degree, but 
in simulation often lead to a failure to find a solution (if it exists).  The performance value 
of individual components at the given context allows the system to be analyzed for 
possible failure of the solution routine. 
5.7.1.3 The Multiple VapCyc Contexts 
VapCyc has five built in contexts: Residential Air-Conditiong (RAC), Residential 
Heat Pump (RHP), High Temperature Refrigeration (HT), Medium Temperature (MT) 
Refrigeration and Low Temperature Refrigeration (LT).  Currently, all contexts have 
fixed values for condenser subcooling, 10?F/5.56?C; evaporator superheat, 7?F/3.88?C; 
compressor discharge superheat,  37.5?F/20.8?C; and heat exchanger pressure drop of 3% 
 
107 
of the inlet value.  The evaporating and condensing pressures are given in Table 5-6 in 
the form of evaporating and condensing temperatures. 
 
Figure 5-11: A Nominal Air-Conditioning Operating Point for VapCyc 
Context 
Evaporating 
Temperature (?F/?C) 
Condensing 
Temperature (?F/?C) 
Residential Air-Conditioning 45.0 / 7.2 130.0 / 54.4 
Residential Heat Pump 30.0 / -1.1 110.0 / 43.3 
High Temperature 
Refrigeration 
20.0 / -6.7 120.0 / 48.9 
Medium Temperature 
Refrigeration 
-25.0 / -31.7 120.0 / 48.9 
Low Temperature 
Refrigeration 
-40.0 / -40.0 105.0 / 40.6 
Table 5-6: VapCyc Evaporating and Condensing Temperatures 
 
108 
5.7.1.4 Scaling With the Context 
The simulation equations and variables are scaled to the context.  The nominal 
evaporating and condensing temperatures ( (), ()
Evap sat Evap Cond sat Cond
PPTPPT==) provide 
good choices to replace P
High
 and P
Low
 in equation (5.12), as they represent the highest 
and lowest pressures in the air-conditioning context.  Likewise the compressor discharge 
and evaporator inlet enthalpy values are the highest and lowest of the air-conditioning 
context.  The air-conditioning context is for steady state, and thus each component has 
the same flow rate.  Therefore the compressor nominal flow rate is chosen as the 
characteristic flow rate to scale mass flow residuals.  The evaporator?s nominal capacity 
is a good choice for a characteristic heat transfer value to scale energy balances with.  
The condenser?s nominal heat rejection could be used just as easily.  Also, at the given 
context each component will have a nominal charge.  The nominal system charge is used 
to scale the system mass balance.  Table 5-7 gives a full set of the scaled independent 
variables and equations used in VapCyc. 
5.7.2 Control of VapCyc 
The power of a simulation tool over many operating conditions for many components 
is useful, but the need to specify the system completely is problematic when design is the 
objective.  VapCyc implements a control, which allows a user to specify an objective 
variable, and a variable to control it.  For example, if one wants to size an expansion 
orifice to give a desired level of superheat at the evaporator outlet, VapCyc allows the 
user to specify the orifice diameter as an independent variable and the evaporator 
superheat as a dependent variable. 
 
109 
Variable / Equation Scaled Variable / Equation 
Junction Pressures, P* 
EvapCond
Evapj
j
PP
PP
P
?
?
=
*
 
Junction Enthalpy 
inEvapedisch
inEvapj
j
hh
hh
h
,arg
,*
?
?
=  
Junction Mass Balance 
0
1
=
?
=
Compl
Nports
p
p
m
m
 
 
 
Junction Energy Balance 
0
,
1
,
1
=
?
?
?
=
=
Evap
poOutletNports
po
po
InletNports
pi
ppi
Q
h
m
hm
 
 
 
System Mass Balance 
0
1
=
?
?
=
Charge
ChargeSystemCharge
sNComponent
i
i
 
Table 5-7: Scaled Independent Variables and Balance Equations 
5.8 VapCyc Examples of Capabilities and Verification 
With the capability to change component models, as well as expand upon the basic 
system configuration, VapCyc represents a generalized energy system modeling approach 
to the standard vapor compression refrigeration cycle.  VapCyc offers more utility than 
existing v-c simulations due to its set of features, namely, charge management, 
component inter-changeability, multiple refrigerants and seasonal performance. 
Several examples are presented to validate the techniques described in this chapter 
which are used within VapCyc.  Examples include: the affect scaling the independent 
variables and balance equations, the usefulness of the context as seen in the interactive 
mode, simulation and verification against actual experimental data, and seasonal 
performance of a system.  The flexibility of the tool is demonstrated by modifications of 
the system configuration, by utilizing several evaporators in parallel and two compressors 
 
110 
in series.  Lastly, VapCyc is used to perform a simple system optimization versus system 
charge and condenser air flow velocity. 
5.8.1 Verification of the Model 
The model is ?solved? when the Junction Solver resolves all junction level mass and 
energy balances, in addition to the system mass balance.  To confirm that the solution is 
indeed a physically meaningful one, P-h diagrams are generated and mass and energy 
balances are calculated.  Three cases are examined for verification 
1. The standard Alpha 3600 Air Conditioning System 
2. The standard Alpha 3600 Air Conditioning System with the condenser inlet 
air set to 40?C 
3. The standard Alpha 3600 Air Conditioning System, with a hot gas bypass 
valve installed between the compressor discharge and the evaporator inlet 
In all cases, the convergence tolerance is set to 10
-4
 x compressor nominal capacity 
for energy (3600 W x 10
-4
 = 0.36 W) and 10
-4
 x compressor nominal flow rate (0.024 
kg/s x 10
-4
 = 2.4 x 10
-6
 kg/s). 
 
Junction Balance (Mass kg/s x 10
-6
, Energy W) 
 
1 2 3 4 
Test # 
Mass Energy Mass Energy Mass Energy Mass Energy 
I 
-0.565 -0.15 0.00934 0.01 0.00644 0.11 0.407 0.05 
II -1.67 0.11 -0.059 -0.03 0.180 -0.69 1.55 0.14 
III -1.50 0.29 0.150 -0.0021 0.059 0.0013 -1.74 -0.18 
Table 5-8: Mass and Energy Balances for Test Configurations 
 
111 
Table 5-8 displays the final junction mass and energy balances for the three tests, 
numerically verifying the solution. 
Figure 5-12 illustrates the results of the three cases used to validate VapCyc in 
graphical form using a P-h diagram.  The shape of the P-h diagrams validates the model, 
being the expected shape for the four junction systems.  The behavior of the system is as 
expected, and is seen increase in system condensing pressure for the increase in 
condenser air temperature, and the increase in suction pressure for the hot gas bypass 
configuration.  The shape of the P-h diagrams, coupled with the mass and energy 
balances shown in Table 5-8 validate the simulation approach. 
 
Figure 5-12: P-h Diagrams for Various Alpha 3600 Configurations 
 
112 
5.8.2 Verification of the Scaling 
To verify the scaling is appropriate and validate the use of the context, a simple 
numerical experiment is performed.  Two systems are simulated using the nominal 
charge, as provided by the component for the Air Conditioning context.  The first system 
is the standard Alpha 3600 system (each of the four components is of the nominal 3.6 kW 
variety), while the second system is the Alpha 3600 system, but with an Alpha 10000 
condenser.  This second system is referred to here as the ?hybrid system?.  The 
simulation is performed using a NEWTON type non-linear equation.  The significance of 
the NEWTON type solver is that a Jacobian matrix is calculated and a correction vector 
is calculated using this matrix.  Large condition numbers for this matrix can result in 
erroneous solutions to this correction vector.  The simulations are carried out with (a) 
with no scaling at all and (b) utilizing the Air-Conditioning context for scaling the 
independent variables and residual equations.  In both cases the initial guess is provided 
to the system from the Air-Conditioning context. 
Figure 5-13 gives a P-h diagram with three cycles overlaid.  The first represents the 
VapCyc Air-Conditioning context, which serves as the initial guess for the two 
simulations, the second, is the solution of the Alpha 3600 system, with its nominal 
charge, and the third is the Alpha 3600 system with an Alpha 10000 condenser in place 
of the original.  It is noted here, that as expected the Alpha 3600 system solution lies very 
near the VapCyc Air-Conditioning context, and except for the condensing temperature 
(pressure), the hybrid system does as well. 
 
113 
 
Figure 5-13: P-h Diagrams of Numerical Experiment Systems 
Figure 5-14 shows a result of the Alpha 3600 system simulation.  Two main points 
are evident, namely (1) The condition number of the Jacobian Matrix used for the 
NEWTON type solver at each iteration is of order 10
6
 (O(10
6
)) for the non-scaled case 
and of O(10) for the scaled case, and (2) the non-scaled case requires 6 iterations, while 
the scaled case requires only 1. 
 
Figure 5-14: Jacobian Condition Number vs. Iteration Number 
 
114 
Figure 5-15 illustrates a result of the hybrid system simulation, with the larger 
condenser. 
 
Figure 5-15: Jacobian Condition Number vs. Iteration Number 
Points evident from this figure are: (1) The condition number of the Jacobian Matrix 
used for the NEWTON type solver at each iteration is again of order 10
6
 (O(10
6
)) for the 
non-scaled case and of O(10) for the scaled case, and (2) the non-scaled case requires 12 
iterations, while the scaled case requires only 7.  Clearly, the numerical solver must work 
harder to solve this system of equations provided by this different set of components, and 
the use of the Context makes the computations easier. 
The power of the Context as both a vehicle for an initial guess and a scaling tool is 
demonstrated by the ability to reduce the computational power to solve the problem and 
its ability to find a solution where a non-scaled problem may not be able to. 
5.8.3 Examining the VapCyc Context 
Figure 5-2 illustrates the main VapCyc interface.  The screen shot shows the VapCyc 
calculation page, in the VapCyc RAC context.  Recalling the ?Alpha? series components 
 
115 
were designed for air-conditioning application, using the Alpha3600 series for each 
component results in a nominal 3.6 kW capacity at the design point, which corresponds to 
the VapCyc RAC context. 
Note, the option to use the default nominal charge is selected, so the sum of the 
nominal charge for each component (shown in the table in the bottom left corner) is equal 
to the specified system charge (in the upper right hand corner). 
Table 5-9 shows the changes applied to the system when low temperature context 
(LT) is utilized.  In this particular case, where both the evaporator and condenser have air 
as the secondary fluid, the only changes are the two secondary air temperatures. 
 
Air-Conditioning Context 
Low Temperature 
Refrigeration Context 
Condenser  Air 
Temperature (?C) 
35.00 35.00 
Evaporator Air  
Temperature (?C) 
26.67 -23.00 
Table 5-9: Air Temperatures of A/C and Refrigeration Context 
Figure 5-16, taken from the main VapCyc screen, details the Alpha3600 components, 
subject to the LT context.  Each of the components behaves differently than when subject 
to the residential air conditioning (RAC) context. 
 
Figure 5-16: Performance of Alpha3600 Components: 
Low Temperature Refrigeration Context 
The compressor, operating at a much lower suction pressure, has a drastically reduced 
capacity.  The condenser, at nearly the same nominal operating point changes little in 
nominal capacity, and nominal charge contained within the condenser is virtually 
 
116 
indistinguishable.  The changes in evaporating and condensing temperatures results in a 
slightly increased expansion device capacity.  The evaporator experiences large 
differences in both nominal capacity and nominal charge, due to the large density 
differences between the two operating points.  Comparing the nominal capacities of the 
components, the disparity between them is evident, and a solution to the model is 
unlikely; as expected, VapCyc cannot find a valid solution.  Figure 5-17 shows a hybrid 
cycle of Alpha components, in the LT context.  The Alpha 3600 orifice is replaced with 
the Alpha 500 orifice and the Alpha 3600 condenser is replaced with the Alpha 1240 
condenser (smallest).   
 
Figure 5-17: Low Temperature Solution to an Alpha Hybrid System 
The VapCyc context has allowed the system to re-specify its operating condition from 
air-conditioning to low temperature refrigeration, as well as provide an initial guess value 
to the solver.  Utilizing the nominal capacity allows for the identification of components 
which may offer difficulty in the cycle operating in an acceptable range, due to mis-
 
117 
matched sizes.  In this case the components are replaced manually to match the nominal 
capacities at the LT context better.  Using the nominal capacity (or deviation in nominal 
capacity from other components) as a figure of merit, the process may be automated, 
allowing VapCyc to exchange equipment automatically.  The ability of the context to 
offer insight to component compatibility as well as seeding of an initial guess to a solver 
gives VapCyc great power. 
5.8.4 Comparison with Test Data 
The comparison of the test data (Hwang, et al. 2004), and VapCyc simulations using 
TEST component models (see Appendix 1) is presented here for verification purposes.  A 
main focus of the test project is the performance of the entire system, for a given 
refrigerant, as a function of the system charge, at constant evaporator superheat (TXV 
controlled). 
The only degree of freedom in the test system experiment is the system charge.  Two 
issues make comparison of refrigerant charge inventory between simulation and 
experiments difficult, namely (1) simulation models (Rice, 1987) can vary greatly in 
estimation of charge depending upon the correlations used to calculate them, and (2) as 
no information regarding piping dimensions is included in the report, the physical system 
piping connecting the various components are not accounted for in the simulation.  The 
liquid line between condenser and expansion device in particular can hold a considerable 
amount of refrigerant. 
Consequently, the VapCyc simulations are carried out at constant superheat, and 
refrigerant charge is varied to match the value of sub cooling at the condenser exit from 
 
118 
the given experimental test point.  Results are given graphically in Figure 5-18 - Figure 
5-21. 
 
Figure 5-18: Simulated vs. Measured Evaporator Heat Load 
Figure 5-18 demonstrates that the VapCyc simulation predicts the evaporator heat 
load within 5%, even for the part load test.  The point corresponding to the ?test 1? 
(Appendix I) data point agrees within 5%.  This is not surprising, as the components were 
designed based upon this point.  All other points slightly over predict the system?s 
cooling performance. 
 
Figure 5-19: Simulated vs. Measured COP 
Figure 5-19 illustrates results similar to that of the system?s cooling performance, 
namely the system COP is predicted within or near 5%.  This is expected, as the cooling 
 
119 
load results are predicted well, and the compressor?s off-design performance (see 
appendix 1) is very similar to the measured compressor off-design performance. 
The good agreement of measured versus simulated compressor performance is also 
evident in Figure 5-20, where the system mass flow rate agrees with measured data to 
within 5%. 
 
Figure 5-20: Simulated vs. Measured Mass Flow Rate 
Figure 5-21 illustrates how difficult it is to simulate refrigerant charge inventory.  To 
compare refrigerant charge in the system, a crude ?charge correction? was made to the 
measured data.  The charge correction is taken as the difference between the measured 
system charge at ?test 1? from the test data (appendix 1) and the virtual charge from each 
component at ?test 1?.  This charge correction is assumed to be constant for each test 
point. 
 
120 
 
Figure 5-21: Simulated vs. Measured Refrigerant Charge 
As expected, the ?test 1? data point has perfect agreement between simulation and 
experiment; however, all remaining points disagree significantly.  Clearly, if 
experimental and simulated charge management is to be in agreement, the components 
themselves must be modeled more accurately with respect to charge inventory, and 
piping must be accounted for. 
To gain some insight to the difference between the simulation results and 
experimental results, Figure 5-22 presents the measured and simulated system suction 
and discharge pressures for the R-290 test data. 
 
Figure 5-22: Suction and Discharge Pressure for R290 Test System 
 
121 
The figure reveals that the experimental pressures are not varying from one test to the 
next in the same manner as the simulation pressures.  This phenomenon is consequence 
of two main issues: (1) refrigerant charge inventory models different from reality, and (2) 
off design component performance differing from experimental performance. 
Table 5-10 and Table 5-11 provide numerical comparisons of both the R-290 system 
and the R-404a system  Comparisons of the simulated evaporator heat load, Q
Evap
, COP, 
mass flow rate, corrected system charge, suction and discharge pressures to the measured 
values from the Test Project are provided. 
  
Q
Evap
 
(kW) 
COP 
Massflow 
(g/s) 
Charge
8
 
(g) 
P
Suction
 
(kPa) 
P
Discharge
 
(kPa) 
Measured 11.44 1.935 41.70 365 407.2 1764.9 
Simulated 11.48 1.927 42.67 365 405.7 1781.1 
Test 1
 
% Deviation -0.35 0.41 -2.32 0.0 0.37 -0.92 
Measured 11.59 1.943 41.62 456 407.6 1786.1 
Simulated 11.50 1.890 42.65 385 406.1 1829.6 
Test 2
 
% Deviation 0.78 2.73 -2.48 15.6 0.37 -2.44 
Measured 11.74 1.940 41.75 546 409.3 1820.8 
Simulated 11.47 1.832 42.22 400 406.3 1884.4 
Tes
t
 3
 
% Deviation 2.30 5.57 -1.13 26.74 0.74 -3.49 
Measured 11.59 1.870 41.62 637 407.6 1786.1 
Simulated 11.43 1.739 42.35 423 406.6 1985.0 
Test 4
 
% Deviation 0.52 7.00 -1.75 33.60 0.25 -11.1 
Measured 13.47 3.911 42.73 456 396.4 1238.9 
Simulated 12.97 3.800 42.70 365 396.1 1250.5 
Tes
t
 5
 
% Deviation 3.71 2.84 0.07 19.96 6.89 -0.94 
Table 5-10: Measured vs. Simulated Data for R-290 Test System 
 
 
 
 
 
 
                                                 
8
 Corrected Charge 
 
122 
  
Q
Evap
 
(kW) 
COP 
Massflow 
(g/s) 
Charge
9
 
(g) 
P
Suction
 
(kPa) 
P
Discharge
 
(kPa) 
Measured 10.94 1.669 100.1 866 503.1 2422.5 
Simulated 10.94 1.668 100.0 866 405.7 1781.1 
Test 1
 
% Deviation 0.0 0.06 0.10 0.0 19.36 26.48 
Measured 11.46 1.709 99.9 1320 504.2 2515.1 
Simulated 11.10 1.670 101.2 920 406.1 1829.6 
Test 2
 
% Deviation 3.14 2.28 -1.30 30.30 19.46 27.26 
Measured 11.41 1.621 98.5 1773 504.3 2709.6 
Simulated 11.22 1.690 99.2 1020 406.3 1884.4 
Te
st 3
 
% Deviation 1.67 -4.26 -0.71 42.47 19.43 30.45 
Measured 13.75 2.558 100.0 1320 482.9 1743.3 
Simulated 13.25 2.480 100.9 920 406.6 1985.0 
Test 4
 
% Deviation 3.64 3.05 -0.90 30.30 15.80 -13.86 
Table 5-11: Measured vs. Simulated Data for R-404a Test System 
5.8.5 Optimization of Charge 
A goal of the test project is the determination of the optimal charge for a given set of 
components at constant superheat at the evaporator exit, with COP being the objective 
function.  In the test project, several data points are taken for different system charges.  
As demonstrated previously, meaningful comparisons between measured and simulated 
charge are difficult to achieve with the current data.  Figure 5-23 offers comparison of 
simulated versus measured data by illustrating the optimum charge as a function of 
condenser exit sub-cooling. 
                                                 
9
 Corrected Charge 
 
123 
 
Figure 5-23: Measured and Simulated COP vs. Condenser Sub-cooling 
The values of COP are similar in both cases, but the optimum COP is shifted to a 
lower value of sub-cooling for the simulation.  This again is a consequence of both the 
component models off design performance and charge inventory behavior.  The 
refrigerant charge inventory of a system is strongly related to the system?s internal 
volume, which was not a focus of the of the current component models.  To achieve 
better agreement, more data with respect to component charge management is needed, 
and the components need to be modeled to perform more accurately off the design point.  
The evaporator capacity, mass flow, suction and discharge pressure agree within 5% for a 
given evaporator superheat and condenser subcooling.  Also, the character of the 
refrigerant charge curve with respect to condenser subcooling is consistent with test data.  
This demonstrates an opportunity to use VapCyc to tune the component models to help 
accurately simulate component charge as a function of operating point. 
 
124 
5.8.6 SEER Example 
Determining the SEER of a given system using the standard requires the 
determination of the system?s operation at two steady state operating points, namely the 
ASHRAE ?A? and ASHRAE ?B? rating points as defined in the standard.  For reference 
Figure 5-24 illustrates the rating points. 
 
Figure 5-24: ASHRAE "A" and "B" Rating Conditions for Air-Conditioning Systems 
The performance of the Alpha 3600 air-conditioning system, with the nominal system 
charge (0.467 kg) at the ASHRAE ?A? and ?B? conditions is given in Table 5-12.  Using 
these values for system performance, in accordance with ASHRAE standard 210/240-94 
yields a system SEER of 13.76 
ASHRAE Test Q (Btu/hr / kW) Power (Btu/hr / kW) COP 
A 12,490 / 3.66 3,316 / 0.92 3.77 
B 12,523 / 3.67 2,771 / 0.82 4.52 
Table 5-12: Alpha 3600 System  Performance at ASHRAE "A" and "B" Conditions 
 
Figure 5-25 illustrates the variation in the Alpha 3600 SEER as a function of system 
charge.  It is noted here that as system charge increases, system superheat at the 
compressor suction port decreases, as sub-cooling at the condenser exit increases.  The 
increase in SEER is due to this increase in sub-cooling.  The increase in COP is seen to 
 
125 
level off at higher charges, as the condensing temperature of the system decreases and the 
condenser (refrigerant) exit temperature approaches the ambient air temperature cooling 
it.   
 
Figure 5-25: Alpha 3600 System SEER vs. System Charge 
For a given operating point, the charge has a maximum value, as some degree of 
suction superheat must be maintained.  The Alpha 3600 system must not exceed 0.57 kg 
of refrigerant charge for the ASHRAE ?B? in order to maintain some degree of suction 
superheat. 
5.8.7 Application Examples: Expandability of the Simple System 
5.8.7.1 Evaporators in Parallel 
The object oriented structure of the VapCyc simulation allows for the simple four 
component vapor compression system to be expanded at run time.  The object oriented 
nature of the Junction Solver, which is the underlying solution engine allows a single 
solver to be utilized to solve many different systems.   
 
126 
A common practice in refrigeration and air conditioning is to server several 
evaporators off of a single compressor and condenser.  Figure 5-26 illustrates a ?split? 
system, where five evaporators of varying size are connected in parallel to a single 
compressor, condenser and expansion orifice in series. 
 
Figure 5-26: A Nominal 20 kW Alpha Split Air Conditioning System 
Such a system may serve as an air conditioning system, where each evaporator would 
serve as the cooling coil for individual rooms requiring different cooling loads.  The 
 
127 
system illustrated in Figure 5-26 is a nominal 20 kW system, consisting of the Alpha 
Compressor 20000, Alpha Condenser 20000, Alpha Orifice 20000 and a set of 
evaporators, Alpha 1240, 2440, 3600, 5000, and 10000.  The system is configured in 
VapCyc, and the system charge is chosen to achieve acceptable levels of compressor 
suction superheat and condenser subcooling. 
Table 5-13 displays the mass flow rate, power consumption and heat rejection by the 
individual components and the total system. 
Component 
Mass Flow 
Rate (g/s) 
Power Out (kW) Heat Out (kW) 
Evaporator 1240 7.35 0.0 -1.16 
Evaporator 2440 14.46 0.0 -2.29 
Evaporator 3600 21.40 0.0 -3.39 
Evaporator 5000 29.67 0.0 -4.70 
Evaporator 10000 59.35 0.0 -9.41 
Compressor 20000 132.22 -5.45 0.0 
Condenser 20000 132.22 0.0 26.41 
Orifice 20000 132.22 0.0 0.0 
System Total 132.22 -5.45 5.46 
Table 5-13: Results for Nominal 20 kW Alpha Split Air Conditioning System 
As expected, the total flow rate through the evaporators is equivalent to that of the 
other three components.  The evaporators each deliver different amounts of cooling to 
their specific zones, and the overall system energy balance is obeyed.  This example 
demonstrates the expandability of VapCyc, and the underlying Junction Solver, and thus 
the flexibility of the simulation tool and simulation technique. 
5.8.7.2 Compressors in Series 
The Delta compressor presents a more realistic behavior of a compressor?s isentropic 
efficiency versus its compression ratio.  Figure 5-27 illustrates the performance of two 
systems as the ambient air temperature for the condenser varies.  Each system contains an 
 
128 
Alpha3600 condenser, Alpha3600TXV (set for 10K superheat) and an Alpha3600 
evaporator.  One system contains an Alpha3600 compressor and the other a Delta3600 
compressor.  Recall, the Alpha compressor has a constant isentropic efficiency of 0.65, 
and the Delta compressor has a parametric isentropic efficiency, which varies with the 
compressor compression ratio, and has a value of 0.65 at the compression ratio of 3.42.  It 
is shown in the figure that the system COP varies with the ambient air temperature, and 
this effect is more pronounced when the (more realistic) Delta compressor is used. 
A combination of two Delta compressors in series can be used, rather than a single 
Delta compressor, when the rotational speed of the second compressor is reduced.  
Simulation demonstrates an increase in system COP by utilizing two delta compressors in 
series, rather than one.   
Figure 5-27 illustrates the COP versus ambient air temperature of a system utilizing 
the single Delta compressor and two Delta compressors in series.  The figure 
demonstrates that two Delta compressors in series can provide improvements in system 
COP at higher than design ambient temperatures. 
It is illustrated in the figure that the system COP curve shifts up or down, depending 
upon the second compressor?s rotational speed.  The design speed of the second 
compressor would therefore be decided based upon the off design operation of the 
system. 
 
129 
 
Figure 5-27: System COP versus Ambient Temperature 
5.8.8 Parametric Study of COP With Respect to Charge and Air Velocity 
Consider a situation where a set of vapor compression refrigeration components are 
supplied for a given application, such as the ASHRAE ?A? condition, per the ARI 
standard 210/240.  Specifically, consider the case of an Alpha 3600 compressor, Alpha 
3600 Evaporator, Alpha 3600 Orifice and an Eta 3600 condenser. 
VapCyc can be used to optimize the system.  Figure 5-28 is a representation of the 
performance of such a system. 
The performance is viewed as a function of system charge and condenser air velocity.  
The figure has lines of constant COP flooded with color and lines of constant cooling 
 
130 
capacity.  Evident from the figure is the two dimensional maximum in COP.  
3
7
7
5
3
7
5
0
3
6
0
0
365
0
3
7
0
0
3
7
1
8
.
2
5
System Charge (g)
Ai
r
t
V
e
l
o
c
i
t
y
(
m
/
s
)
50 60 70 80
1
2
3
COP
3.5
3.4
3.3
3.2
3.1
3
2.9
2.8
2.7
2.6
2.5
2.4
2.3
2.2
2.1
Capacity (W)
 
Figure 5-28: COP and Capacity versus Charge and Air Velocity 
This system has a maximum COP of 3.53 when the condenser air velocity is 2.19 m/s 
and the system charge is 64.7g. 
5.9 Summary 
This chapter summarizes VapCyc, a comprehensive and robust tool for steady-state 
vapor compression refrigeration system simulation.  The objective of the tool is for 
analysis, design and optimization of vapor compression refrigeration systems and their 
components.  To mimic the behavior of real vapor compression refrigeration systems, a 
list of features including: charge management, component inter-changeability, a library of 
available refrigerants, a variety of unit systems, a series of physical variables, and the 
ability to utilize economic variables are integral to VapCyc.  The tool can be used in an 
 
131 
interactive mode for simulating a single operating point as well as an automated mode for 
parametric analysis. 
The VapCyc calculation engine is an object oriented entity, based upon the 
generalized energy system simulation principals given in chapter 2 of this thesis.  The 
component library is built upon the component models given in chapters 4-6 of this 
thesis. 
A major simulation challenge is the solution of the mathematical model.  This thesis 
develops a simulation framework, and tool, the Junction Solver, to create systems in a 
general manner.  The challenges of providing an initial guess and proper scaling of the 
conservation equations and independent variables fall upon the simulation program 
employing the Junction Solver.  In this chapter, a novel concept called the system context, 
which is a logical extension of rating equipment, is developed.  The system context is a 
link between the specific system being simulated and the general nature of the Junction 
Solver.  Several contexts are defined for several different vapor compression refrigeration 
applications.  The context allows for communication between the system level solver, and 
the component level models.  The system level solver can utilize information provided by 
components for the specified context to intelligently develop an initial guess, and 
diagnose potential disastrous system configurations, composed of mis-matched 
component models. 
VapCyc allows the research concerns of individual users to be addressed.  These 
concerns include: 
? The analysis and performance of particular pieces of hardware, over a range 
of operating conditions 
 
132 
? A forum for sharing of performance information, while maintaining propriety 
of such information 
? A tools to access/publish/share this performance information 
? An assembled libraries of components which can be used to assemble full 
systems 
? The ability to explore the performance of specific systems over a range of 
operating conditions to which the system may be subject 
? Tools to access/publish/share this system performance information 
? An environment for re-usable experiments/simulations 
? Component level experiments/simulation 
? System level experiments/simulation 
? Flexible and adjustable system design for experiments/simulations 
? Interchangeable components 
? Complete separation of component models from system models, enabling 
o Research to be conducted on systems, with individual component 
models being independent variables to the experiment/simulation 
o Research to be conducted on individual components 
o Research to be conducted on simulation solvers independently from 
components 
? User friendly software 
? Common simulation framework over multiple system simulations 
? A component/system simulation environment that closely mimics that of the 
physical experimental environment 
 
133 
5.10 References 
Air-Conditioning and Refrigeration Institute, 1994, Unitary Air Conditioning and Heat 
Pump, ANSI/ARI Standard 210/240. 
 
Air-Conditioning and Refrigeration Institute, 1999, Positive Displacement Refrigerant 
Compressors and Compressor Units, ANSI/ARI Standard 540. 
 
ASHRAE Handbook: Equipment, American Society of Heating, Refrigerating and Air-
Conditioning Engineers, 1988 
 
Copeland Corporation, Division of Emerson Climate Technologies, On-Line Catalog, 
2004, http://www.copeland-corp.com/americas/entercataloginfo.htm 
 
Davis, G.L., Scott, T.C., 1976.  ?Component Modeling Requirements For Refrigeration 
Systems Simulation?, Proceedings of the 1976 Purdue Compressor Technology 
Conference, West Lafayette, USAA, July 1976, pp. 401-408 
 
Fischer, S.K., and Rice, C.K., 1986  ?System Design Optimization and Validation for 
Single-Speed Heat Pumps?, ASHRAE Transactions, Vol. 91, Part 2 
 
Hiller, C.C., Glicksman, L.R., 1976, ?Improving Heat?Pump Performance Via 
Compressor Capacity Control ? Analysis And Test?.  Volumes I and II, Energy 
Laboratory Reports MIT-EL 76-0001 and MIT-EL 76-0002, Department of Mechanical 
Engineering, Massachusetts Institute of Technology, Massachusetts, USA (January 1976) 
 
Hwang, Y., Dea-Huan, J., Radermacher, R., ?Comparison Of Hydrocarbon R-290 And 
Two Hfc Blends R-404a And R-410a For Medium Temperature Refrigeration 
Applications?, ARI: Global Refrigerant Environmental Evaluation Network (GREEN) 
Program, Final Report, March 2004 
 
Kreith, F., 1999, Mechanical  Engineering Handbook, Boca Raton: CRC Press LLC, 
1999 
 
Levins, W., Rice, C.K., Baxter, V.D. 1996. ?Modeled And Measured Effects Of 
Compressor Downsizing In An Existing Air Conditioner/Heat Pump In The Cooling 
Mode?, ASHRAE Transactions, V 102, Pt 2 
 
NIST, 2001, NIST Thermodynamic and Transport Properties of Refrigerants and 
refrigerant Mixtures-REFPROP Version 7.0, National Institute of Standards and 
Technology, Boulder, Colorado, USA. 
 
Rice, C.K., Fischer, S. K., Ellison, R.D., Jackson, W.L. 1981.  ?Design Optimization of 
Conventional Heat Pumps: Application to Steady-State Heating Efficiency,?  ASHRAE 
Transactions. Vol. 87, pp1037-1054 
 
 
134 
Rice, C.K., 1987 ?The Effect of Void Fraction Correlation and Heat Flux Assumption on 
Refrigerant Charge Inventory Predictions?, ASHRAE Transactions. Vol. 93, part 1, pp 
341-367 
 
Richardson, 1998.  Performance Models For Use in Computer Simulation of Energy 
Systems, Masters Thesis, Department of Mechanical Engineering, Marquette University, 
1998 
 
Stoeker, W., Refrigeration and Air Conditioning, McGraw-Hill, 1958 
 
Tecumseh Products, Model Chart, 2004, http://www.tpc-nacg.com/mdlnmbr1.htm 
 
Tandon, A., Object-Oriented Modeling of Vapor Compression Systems and Components, 
M.S. Thesis,Department of Mechanical Engineering, University of Maryland, College 
Park, 1999. 
 
Webb, R.L. 1994, Principals of Enhanced Heat Transfer, John Wiley & Sons, New York. 
 
 
135 
6. Conclusions 
The objectives of this thesis are: 
1. Develop and create a simulation framework, and subsequent simulation tool in 
an object oriented manner, which provides the ability to simulate experiments 
in a laboratory.  The framework is developed with the following fundamental 
objectives: 
a. Serve as a functional starting point which will appeal to the individual 
objectives of each individual user with a wide range of backgrounds: 
industrial uses (development of models and/or real equipment), 
academic professionals (utility of the simulation tool without an 
intimate understanding of its inner workings), and students (re-using of 
existing software to facilitate the growth of the project(s))) 
b. Serve as a simulation tool such that each component within the system 
is a stand-alone simulation model, enabling the analysis and if desired 
optimization of individual simulation models and/or the system as a 
whole 
2. Verify the simulation tool with experimental data 
6.1 Thermal Energy System Simulation Framework 
The generalized methodology is built upon the structures of: components, ports and 
junctions.   
Components are individual energy systems themselves, which are connected by the 
user at run time to create a larger system.  Each component has ports with which 
 
136 
communication exist via junctions.  Similar in nature to real world devices, boundary 
conditions are applied to a component, and a steady state solution results. 
Through the iteration on junction independent variables, such as pressure and 
enthalpy, system level conservation equations are solved, consequently simulating the 
system. 
6.2 Basic Components 
In this thesis, the energy interactions for the simulation of the vapor compression 
refrigeration system are limited to thermo fluid interactions, such as thermal energy, heat 
transfer and work.  Components are developed as sub-systems to be combined into larger 
systems. 
Each type of component has unique characteristics which govern its specific 
relationships between boundary conditions; pressure and inlet energy state, to mass flow, 
heat transfer and work output.  Individual components determine outlet conditions as well 
as heat/work interactions via a common set of boundary conditions.  A generalized 
energy system simulation using such components can be constructed utilizing system 
level intensive variables as the component boundary conditions.  In this particular case, 
the system level intensive variables; junction pressure and junction exit specific enthalpy 
are employed and communicated to components, via ports. 
6.3 Component Models 
Building energy systems from components and utilizing these energy systems as 
components in a larger system presents mathematical challenges for simulation.  Two 
main challenges are solving the sub-system models to achieve a desired solution for a 
 
137 
system, and the mathematical character of the partial differential equations which 
describe the physics being modeled. 
The main rule to focus upon when modeling components for use in larger energy 
systems are determined to be: proper scaling of independent variables and balance 
equations.  The scaling is shown to affect the solution process, particularly when gradient 
based solution methods are employed. 
The independent variable scales must be indicative of the physical phenomena being 
modeled.  This thesis exploits component level boundary conditions for the scaling of 
independent variables. 
6.4 Compound Energy Systems 
Chapter four of this thesis broaches the concept of a compound energy system, which 
is the combination of two or more hydraulically independent systems.  In chapter four, a 
refrigerant circuit and an air circuit are combined to make a heat exchanger.  Later, in 
chapter five, a compressor, condenser, expansion device and evaporator are combined 
into a vapor compression refrigeration system.  Also in chapter five, additional 
components are added, demonstrating the ability to create energy systems in an arbitrary 
fashion by combining individual components. 
Compound energy systems have the same scaling issues as basic components.  The 
scaling parameters for both independent variables and balance equations should be 
indicative of both independent systems.  In the case of the air to refrigerant cross flow 
heat exchanger example, the air entrance enthalpy, coupled with the boundary conditions 
on the refrigerant side allow scaling values to be calculated. 
 
138 
Uncertainty propagation between components and the system generally needs to 
begin at the lowest levels, such that when propagated to the highest level, the 
convergence criterion of the highest system is met. 
The character of the mass flow relations for the compound energy system is elliptic, 
while the energy relations are essentially hyperbolic.  The flow arrangement of the two 
fluids however need not be hyperbolic in the same direction.  Consequently, some form 
of de-coupled mass and energy relations are necessary.  A fractional step method 
provided by Jiang (2003) is examined.  Additionally, more complex methods which 
couple the mass and energy relations of the two systems to varying degrees are 
developed.  Ultimately a robust method for solving the example compound system is 
developed. 
The heat exchanger model developed in this thesis has been validated and verified 
against actual test data. 
6.5 VapCyc 
This thesis reaches a culmination with VapCyc, a comprehensive and robust tool for 
steady-state vapor compression refrigeration system simulation.  The objective of the tool 
is for analysis, design and optimization of vapor compression refrigeration systems and 
their components.  VapCyc serves as a simulation tool such that each component within 
the system is an individual program, enabling the analysis and if desired optimization of 
individual simulation models and/or the system as a whole 
This simulates the behavior of real vapor compression refrigeration systems, where a 
list of relevant component and/or system features includes: charge management, 
component inter-changeability, a library of available refrigerants, a variety of unit 
 
139 
systems, a series of physical variables, and the ability to utilize economic variables are 
integral to VapCyc.  The tool can be used in an interactive mode for simulating a single 
operating point as well as an automated mode for parametric analysis. 
The VapCyc calculation engines, as well as its component library are object oriented 
entities, based upon the generalized energy system simulation principals.   
The components themselves are completely unrelated unless connected to one another 
via junctions.  Therefore, neither the system level balance equations, nor the independent 
variables can be scaled relative to the component models.   
Overcoming the complexity of formulating an initial guess is accomplished through a 
novel concept called the system context, which is a logical extension of rating equipment.  
Several contexts are defined for several different vapor compression refrigeration 
applications.  The context allows for communication between the system level solver, and 
the component level models.  The system level solver can utilize information provided by 
components for the specified context to intelligently develop an initial guess, and 
diagnose potential disastrous system configurations, composed of mismatched 
component models. 
Several examples demonstrate how the solution technique is employed in a flexible 
setting.  Single point simulations, parametric analyses and optimizations are possible.  
Independent variables for the simulations are not limited to component independent 
variables, but the refrigerant, system charge and even component model. 
VapCyc allows the research concerns of individual users to be addressed.  These 
concerns include: 
 
140 
? The analysis and performance of particular pieces of hardware, over a range 
of operating conditions 
? A forum for sharing of performance information, while maintaining propriety 
of such information 
? A tools to access/publish/share this performance information 
? An assembled libraries of components which can be used to assemble full 
systems 
? The ability to explore the performance of specific systems over a range of 
operating conditions to which the system may be subject 
? Tools to access/publish/share this system performance information 
? An environment for re-usable experiments/simulations 
? Component level experiments/simulation 
? System level experiments/simulation 
? Flexible and adjustable system design for experiments/simulations 
? Interchangeable components 
? Complete separation of component models from system models, enabling 
o Research to be conducted on systems, with individual component 
models being independent variables to the experiment/simulation 
o Research to be conducted on individual components 
o Research to be conducted on simulation solvers independently from 
components 
? User friendly software 
? Common simulation framework over multiple system simulations 
 
141 
? A component/system simulation environment that closely mimics that of the 
physical experimental environment 
In this thesis, VapCyc has been validated and verified against actual test data. 
6.6 Future Work 
Future work opened up by this research takes several different forms.  
Mathematically, there is much to be done.  Although the complex nature of the 
momentum and energy relations cannot be determined before run time, there is much 
room in a generalized ?junction solver? to interpret the solution process. 
The generalized energy system process can be either applied to different systems 
and/or expanded.  With energy transport being similar in nature from energy form to 
energy form, it is conceivable that many types of energy be incorporated.  The 
generalized solver which is used to solve the simplified system given in VapCyc, is easily 
altered to include kinetic energy, potential energy electro-magnetic energy.  
VapCyc has many potential uses.  The component swapping capability, combined 
with the evolutionary optimization techniques make it an ideal candidate to create in 
industry wide simulation tool.  Creating a component library from every major 
manufacturer, would allow a systems engineer combining equipment into systems to 
explore thousands of potential systems very quickly. 
6.7 Major Accomplishments 
This thesis accomplished the following objectives: 
1. Develop and create a simulation framework in an object oriented manner, 
which provides the ability to simulate experiments in a laboratory 
 
142 
2. Develop a solution algorithm for the simulation tool where the system 
solution is separate from individual components 
3. Create a heat exchanger simulation, operating with this tool 
4. Verify the simulation tool with experimental data 
5. Create VapCyc, a comprehensive (Graphical User Interface + Junction Solver 
+ Components) steady state vapor-compression refrigeration simulation tool 
The thesis developed and created the Junction Solver.  The Junction Solver is an 
object oriented structure of premised upon the data structures of components, ports and 
junctions which are used to virtually construct thermo thermal energy systems.  As 
detailed in chapter 3, the methodology allows for components to be treated as black box 
entities, which conform to specific input/output standards.  The input/output parameters 
to the component models are intensive thermodynamic variables.  This allows the 
individual component function to be separate from the overall system governing 
equations.  The simulation is then deemed to be multidisciplinary in nature, as each 
component can be considered as an individual discipline, of which the Junction Solver 
(coordinating system) is independent.  The Junction Solver provides the abilities to 
simulate experiments as in a laboratory, such as: 
1. The ability to constrain and specify the mass of working fluid for the 
simulation of system simulation.  For a vapor-compression refrigeration 
application, this is referred to as ?Charge Management?. 
2. The functionality to change, or swap, component models used in a fixed-
component simulation at run-time thus replacing components without the need 
to re-compile the simulation code itself. 
 
143 
3. The functionality to add or subtract components at run time, without the need 
to re-compile the simulation code itself. 
4. Allow user defined component connectivity. 
5. Allow for the selection from several different working fluids/refrigerants. 
6. Flexibility to alter component independent variables, independently from the 
main simulation program. 
The heat exchanger simulation detailed in chapter 4 was used to verify and validate 
the Junction Solver.  The heat exchanger model created to simulate a specific heat 
exchanger, and the model was verified against the experimental data. 
VapCyc, a vapor compression refrigeration system simulation tool was created in 
chapter 5.  As a specific application built upon the generality of the Junction Solver, the 
system context concept was created to overcome the numerical challenges of scaling and 
initial guess values.  Examples are used to demonstrate the validity of the scaling 
provided by the system context and demonstrate the flexibility and expandability the 
Junction Solver provides to a specific application, such as VapCyc. 
6.8 Major Engineering Accomplishments 
This thesis advanced the state of thermo-fluid energy system simulation, specifically 
thermo fluid systems such as heat exchangers, pipe networks and closed loop systems 
such as vapor compression refrigeration.  A generalized framework and simulation tools 
were developed.  Engineering principles were applied to the framework and tools to 
demonstrate its use and develop an advanced vapor compression refrigeration simulation 
tool.  Key engineering accomplishments differentiate this work from previous work, 
namely: 
 
144 
1. The framework has scaling of independent variables and conservation 
equations at the solver level.  This frees the solution algorithm from 
determining application specific variables, allowing development of the 
mathematics separately from engineering. 
2. The selection of dependent and independent variables.  The Junction Solver 
has a set of independent variables, consisting of thermodynamic properties, 
which are intensive variables.  These independent variables are the input to 
sub-system components, which then provide a set of dependent variables, of 
mass and energy flow rate through its ports, as well as the mass of working 
fluid contained within the component at the given operating point.  The 
Junction Solver is resolved using the mass and energy flow rates flowing 
between components, and for closed systems, a total system mass balance.  
The benefits of these choices are: 
a. While typical simulation programs utilize working fluid flow rates as 
independent variables.  This is problematic: 
i. When flows split, and/or flow direction is not known a priori.  
ii. The order of magnitude of the flow throughout the system is 
not known.  For example, a 1kW air conditioner and a 1000kW 
air conditioner have radically different flow rates, yet are 
essentially the same system. 
b. Pressure and Enthalpy are logical choices for system independent 
variables, as they are independent of system size, mass flow rate and 
energy flow rate. 
 
145 
3. A comprehensive, yet flexible simulation framework, which also allows for 
charge management. 
4. With consistent independent variables, the framework does not require 
differentiation of components 
a. Traditional simulations had different dependent/independent variables for 
different types of components.  For example, compressors/valves may 
have boundary pressures and internal state as input and gave mass flow 
rate and exit state as output, while heat exchangers may have mass flow 
rate and inlet state as input and provided outlet state as output.  This 
required a solution algorithm with system configuration knowledge. 
b. The Junction Solver allows components to be connected in any order, and 
allows for a solution algorithm independent of component evaluation 
order. 
c. The system configuration can change arbitrarily at run time, with no need 
to re-compile. 
5. Implementation of the System Context: a vehicle for communication between 
components and the system which they create prior to attempting to resolve 
the simulation. 
a. Simulation program can determine potential component mismatches. 
b. Simulation program can provide appropriate scaling values to the junction 
solver based upon potential system operation. 
c. Simulation program can provide an initial guess to the system solution. 
d. Simulation program can provide an estimate of the system charge. 
 
146 
6. An intuitive interface that is useful to engineers with a diverse set of 
objectives 
a. Component simulation 
b. System simulation 
c. Variation of working fluid 
d. Solver development 
 
147 
A1. TEST System 
A1.1 The TEST Compressor Model 
The TEST Compressor models are modified versions of the ARI 540-99 model, as 
presented in chapter 5.  The modification is slight, correcting for an ?off design? 
operation.  The test points given in the standard all correspond to a certain return gas 
temperature to the compressor.  In reality, a change in temperature at constant pressure 
will change the gas density and thus affect the mass flow rate through the compressor.  
Tandon (1999) introduces a scale factor of the form 
 
tan
(,)
(, )
actual actual
actual s dard
PT
s
PT
?
?
= ,    (A1.1)  
where the actual density (? ) is scales by the design (standard) density for the given 
pressure.  This scale factor, s, is used to correct the model mass flow rate for the 
difference between the actual temperature and the design temperature given in the 
standard. 
Figure A1- 1 and Figure A1- 2 illustrate the comparison of the Green Compressor models 
(both R-404a and R-290) and the measured performance data from the report (Hwang et. 
al 2004).  It is clear from both figures that, at least at the measured data points, the model 
predicts the compressor mass flow rate and power requirement within +/- 5%. 
 
148 
 
Figure A1- 1: TEST Compressor R-404a 
 
 
Figure A1- 2: TEST Compressor R-290 
A1.2 The TEST Orifice Model 
The TEST orifice is based upon the generic orifice, as presented in chapter 5.  The orifice 
model has two degrees of freedom, namely the constant K and diameter D from equation 
(5.3).  Using the first test point from test project data (Hwang, et. al 2004), and assuming 
a value of K=1, the TEST orifice is sized to give the corresponding mass flow rate. 
Figure A1- 3 graphically displays the agreement between experimentally measured and 
calculated mass flow rate for the TEST orifice. 
 
149 
 
Figure A1- 3: TEST Orifice Performance 
The figure shows some scatter in the data, with most points falling within 5% of the 
measured data.  In both the R-404a and R-290 orifice, the deviation for the off-design test 
point is greater than for the others. 
A1.3 The TEST Evaporator Model 
The TEST evaporator is based upon the generic evaporator, as presented in chapter 5.  
The length and diameter of the generic evaporator are sized such that the mass flow rate 
and refrigerant energy increase (refrigeration capacity) agree with the test data at test 
point 1 (see appendix I). 
 
Figure A1- 4: Simulated Versus Measured Mass Flow Rate for TEST Evaporator (a) R-290 and (b) 
R-404a 
 
150 
Figure A1- 4 - Figure A1- 6 show the agreement of the simulated results versus the 
measured test points. 
 
Figure A1- 5: Simulated Versus Measured Superheat for TEST Evaporator 
R-290 and (b) R-404a 
 
 
Figure A1- 6: Simulated Versus Measured Heat Load for TEST Evaporator (a) R-290 and (b) R-404a 
The TEST evaporator simulation agrees with measured data to within 5% for all 
cases, except the R-290 part load case.  The simulated superheat, a reflection of the outlet 
condition for the evaporator does not agree very well.  The large deviation in superheat, 
however does necessarily represent a large deviation in overall performance, as 
demonstrated by the fairly good representation of evaporator heat load.  In all cases 
 
151 
except the R-404a part load test, the simulated heat load agrees with the measured heat 
load within 5%. 
A1.4 The TEST Condenser Model 
The TEST condenser is the condenser modeled in chapter 4.  Results for the condenser?s 
performance are given in that chapter. 
  
 
152 
A1.5 Raw Data 
R290 3.8 lb 
Tair 
Evap 1.59
Tair 
Cond 34.79
Mass 
flow Rate 41.72 Charge* 1.00 
 Pin Pout Tin Tout hin hout Q W 
Power 
Consumption
Compressor 400.87 1764.94 3.19 82.16 583.54 695.74 0.00 4680.47 5.91 
Condenser 1764.94 1667.04 80.07 40.86 690.98 310.16
-
15886.29 0.00 0.00 
Expansion 1595.33 424.16 38.67 -4.32 303.83 303.83 0.00 0.00 0.00 
Evaporator 424.16 407.21 -4.32 0.13 303.83 577.96 11435.75 0.00 0.91 
      
Energy 
Balance 229.92 COP 1.94 
    
R290 4.0 lb 
Tair 
Evap 1.76
Tair 
Cond 34.69
Mass 
flow Rate 41.62 Charge* 1.05 
 Pin Pout Tin Tout hin hout Q W 
Power 
Consumption
Compressor 401.42 1786.14 3.42 83.18 583.91 697.47 0.00 4726.15 5.96 
Condenser 1786.14 1695.40 81.04 39.34 692.58 305.73
-
16099.41 0.00 0.00 
Expansion 1622.57 424.07 37.34 -4.37 299.99 299.99 0.00 0.00 0.00 
Evaporator 424.07 407.60 -4.37 0.39 299.99 578.39 11586.20 0.00 0.90 
          
      
Energy 
Balance 212.94 COP 1.94 
    
      
R290 4.2 lb 
Tair 
Evap 1.67
Tair 
Cond 34.66
Mass 
flow Rate 41.75 Charge* 1.11 
 Pin Pout Tin Tout hin hout Q W 
Power 
Consumption
 
153 
Compressor 403.15 1820.76 3.12 84.08 583.32 698.55 0.00 4810.23 6.05 
Condenser 1820.76 1735.47 81.89 38.01 693.50 301.87
-
16349.45 0.00 0.00 
Expansion 1662.73 425.76 36.14 -4.29 296.54 296.54 0.00 0.00 0.00 
Evaporator 425.76 409.32 -4.29 0.05 296.54 577.73 11738.99 0.00 0.91 
          
      
Energy 
Balance 199.77 COP 1.94 
    
      
R290 4.4 lb 
Tair 
Evap 1.76
Tair 
Cond 34.69
Mass 
flow Rate 41.62 Charge* 1.16 
 Pin Pout Tin Tout hin hout Q W 
Power 
Consumption
Compressor 401.42 1786.14 3.42 83.18 583.91 697.47 0.00 4726.15 6.20 
Condenser 1786.14 1695.40 81.04 39.34 692.58 305.73
-
16099.41 0.00 0.00 
Expansion 1622.57 424.07 37.34 -4.37 299.99 299.99 0.00 0.00 0.00 
Evaporator 424.07 407.60 -4.37 0.39 299.99 578.39 11586.20 0.00 0.90 
          
      
Energy 
Balance 212.94 COP 1.87 
    
      
R290 4.0 lb - PL 
Tair 
Evap 1.59
Tair 
Cond 18.76
Mass 
flow Rate 42.73 Charge* 1.05 
 Pin Pout Tin Tout hin hout Q W 
Power 
Consumption
Compressor 389.41 1238.89 1.53 59.84 581.16 662.54 0.00 3478.05 4.71 
Condenser 1238.89 1130.82 57.95 25.26 658.56 266.25
-
16765.38 0.00 0.00 
Expansion 1080.01 413.11 23.89 -5.26 262.51 262.51 0.00 0.00 0.00 
Evaporator 413.11 396.35 -5.26 -0.36 262.51 577.59 13464.99 0.00 0.91 
 
154 
          
      
Energy 
Balance 177.66 COP 2.86 
 
R404a 10 lb Tair Evap 1.61 Tair Cond 35.50 Mass flow Rate 100.14 Charge* 1.00 
 Pin Pout Tin Tout hin hout Q W 
Power 
Consumption
Compressor 503.06 2422.46 2.16 81.28 368.86 421.72 0.00 5293.63 6.56 
Condenser 2422.46 16.00 79.83 40.94 419.87 259.74 -16035.44 0.00  
Expansion 2175.13 544.01 39.41 -3.67 257.28 257.28 0.00 0.00 
Evaporator 544.01 514.16 -3.67 0.05 257.28 366.56 10943.32 0.00  
         
Energy Balance 201.51 COP 1.67 
          
R404a 11 lb Tair Evap 1.76 Tair Cond 35.15 Mass flow Rate 99.87 Charge* 1.10 
 Pin Pout Tin Tout hin hout Q W 
Power 
Consumption
Compressor 504.19 2515.08 2.03 83.59 368.70 423.26 0.00 5448.74 6.70 
Condenser 2515.08 2384.86 82.00 37.42 421.23 253.78 -16723.82 0.00  
Expansion 2276.97 544.52 36.09 -3.70 251.71 251.71 0.00 0.00 
Evaporator 544.52 515.12 -3.70 -0.04 251.71 366.46 11460.31 0.00  
         
Energy Balance 185.23 COP 1.71 
          
R404a 12 lb Tair Evap 1.84 Tair Cond 35.08 Mass flow Rate 98.54 Charge* 1.20 
 Pin Pout Tin Tout hin hout Q W 
Power 
Consumption
Compressor 504.33 2709.60 2.43 89.40 369.08 427.88 0.00 5794.44 7.04 
Condenser 2709.60 2589.27 87.62 36.96 425.57 252.87 -17017.81 0.00  
Expansion 2481.26 543.50 35.63 -3.75 250.81 250.81 0.00 0.00 
Evaporator 543.50 514.76 -3.75 0.14 250.81 366.64 11413.61 0.00  
         
 
155 
      Energy Balance 190.25 COP 1.62 
   
       
R404a 11 lb - PL Tair Evap 1.66 Tair Cond 18.34 Mass flow Rate 99.98 Charge* 1.10 
 Pin Pout Tin Tout hin hout Q W 
Power 
Consumption
Compressor 482.88 1743.29 0.27 61.27 367.54 407.98 0.00 4042.97 5.38 
Condenser 1743.29 1606.85 59.78 21.70 406.21 229.66 -17650.94 0.00  
Expansion 1527.59 523.96 20.99 -5.20 228.61 228.61 0.00 0.00 
Evaporator 523.96 495.25 -5.20 -0.84 228.61 366.17 13752.73 0.00  
         
Energy Balance 144.75 COP 2.56 
 
156 
A1.6  References 
Hwang, Y., Dea-Huan, J., Radermacher, R., ?Comparison Of Hydrocarbon R-290 And 
Two Hfc Blends R-404a And R-410a For Medium Temperature Refrigeration 
Applications?, ARI: Global Refrigerant Environmental Evaluation Network (GREEN) 
Program, Final Report, March 2004
 
157 
A2. VapCyc ? The Program 
A2.1 Basic Description 
VapCyc is a general-purpose, steady state, vapor compression refrigeration system 
simulation tool.  ISOC aspires to make VapCyc the most advanced vapor compression 
system modeling tool available.  In this instance, most advanced, entitles a combination 
of simulation ability, coupled with functionality, currently not offered (at least in the 
public domain) in the field of vapor compression refrigeration simulation. 
The desired applicability of VapCyc is system simulation, for the purposes of: 
? Analysis of the system, as a function of component performance 
? Analysis of a component, as a part of a system 
? Optimization of system level boundary conditions 
? Optimization of component level parameters 
 
To achieve these goals VapCyc is a simulation of a fixed configuration vapor 
compression system, consisting of: 
? 1 Compressor 
? 1 Condenser 
? 1 Expansion Device 
? 1 Evaporator 
? Expandable Components, both in Parallel and in Series with Four Main 
Components 
? Piping Between Components 
? Suction Line Heat Exchanger 
? Accumulator 
 
This system is general purpose in that each of the components can be selected and/or 
changed, and (if applicable) modified. 
 
158 
VapCyc has a library of components, consisting of ?fixed? models, as well as 
?generic? components (components which allow some or all of their properties to be 
modified).  Fixed components mimic ?real life? components, in that their performance is 
dictated by their geometry, which cannot be changed. 
An example of a fixed versus generic component is a compressor.  A purchased 
compressor has a performance that is dictated by its boundary conditions and 
unchangeable geometry.  A generic compressor allows the user to specify certain 
parameters.  The Generic Compressor supplied with VapCyc has inputs for: 
Displacement Volume, RPM, Volumetric Efficiency and Isentropic efficiency. 
It is desirable to offer a limited number of system level variables (those which would 
controllable with a ?real? machine) that can be modified by the user.  The system level 
variables available in VapCyc are: 
? The refrigerant type 
? System charge 
? System boundary conditions 
o Condenser coolant temperature (If applicable) 
o Condenser coolant flow rate (If applicable) 
o Evaporator heating fluid temperature (If applicable) 
o Evaporator heating fluid flow rate (If applicable) 
A2.2  The Interface 
The main VapCyc screen gives the user a snapshot of the basic system, with pictures 
representing the four main components, as shown in Figure A2- 1.  Evident from the 
main interface page is the ability to change the refrigerant of the system, via a drop-
down-box in the upper right hand corner.  Also in this section, the user can specify the 
refrigerant charge, contained within the system.  The check box to the right of the 
 
159 
refrigerant charge, allows VapCyc to calculate a nominal, or default charge, based upon 
the specified operating context and default charges for the currently selected components. 
 
Figure A2- 1: VapCyc Main Interface Screen 
In the middle of the right hand side of the main screen are some performance 
variables calculated by VapCyc when the simulation is run.  Currently the system?s 
cooling capacity, coefficient of performance (COP), superheat value at the evaporator 
outlet, condenser subcooling and refrigerant massflow are reported.  
The bottom right hand corner of the screen shows performance variables which may 
be of interest.  Physical variables of system cost, weight, volume are calculated as the 
sum of individual component values.  System reliability and manufacturability have place 
holders in the current version of the program, but are not yet fully functional. 
 
160 
A2.2.1 Choosing Components 
Component models are pre-loaded, by default, the ?Alpha3600? series is loaded.  
These components are designed to give a 3.6kW cooling capacity under air-conditioning 
conditions.  Clicking the right mouse button while the cursor hovers over a component 
opens a dialog box allowing the user to select additional components.  Figure A2- 2 
demonstrated the selection of a different compressor from a list of available models. 
 
Figure A2- 2: Changing of a Component in VapCyc 
When the cycle is defined, in terms of component models, refrigerant charge value 
and refrigerant, the system can be simulated by selecting ?Run Current Cycle? from the 
Calculate menu.  Alternatively, the user can press Cntrl-R as a short cut.  At this point the 
solution values are displayed.  Illustrates the result of running the default Alpha3600 
system. 
 
161 
 
Figure A2- 3: The Solution of the Alpha3600 System at the Default Operating Point 
A2.2.2 Editing Component Independent Variables 
Component models may have independent variables which can alter their 
performance.  An example of this is the air temperature of a condenser.  The components 
themselves are responsible for providing the means to alter these values, which can be 
accessed by right clicking the right button, while hovering the mouse over the component 
of interest and selecting the ?edit? function.  Figure A2- 4 illustrates the edit window for 
the Alpha3600 condenser. 
 
162 
 
Figure A2- 4: Edit Menu for the Alpha 3600 Condenser 
A2.2.3 The Parametric Interface 
VapCyc has a parametric interface, which utilizes Microsoft Excel.  Through this  
interface, the system is completely specified.  Component models as well as their 
independent variables are accessed through this interface.  Figure A2- 5 is a shot of the 
Excel interface for VapCyc.  From the figure, it is seen that the component models are 
specified using their name and the system charge is varied for the individual run. 
 
163 
 
Figure A2- 5: Excel Interface for VapCyc's Parametric Input Mode 
A2.3  The Solution Routine 
VapCyc is an application program to simulate the performance of a specific vapor 
compression refrigeration system.  As such, it has a Graphical User Interface (GUI) 
which implements a calculation engine.  The calculation engine is a set of routines which 
receives information from the GUI and give the necessary information to the Junction 
Solver.  After the Junction Solver has resolved the system, the calculation engine then 
retrieves information from the Junction Solver and feeds it to the GUI.  This is achieved 
through three distinct phases of operation, namely a Pre Processor Phase, Solution Phase 
and Post Processor Phase. 
 
164 
A2.3.4 The Pre Processor Phase 
The Pre Processor Phase is an interactive phase that allows the user to specify the 
system configuration and interpret information.  VapCyc employs the specified system 
context to evaluate the current set of components and reports the results to the user.  The 
pre processor is responsible for determining the nominal capacity and charge of a 
component for the specified system context.  The capacity is a rough indicator of the flow 
rate the device will have at the operating component?s operating point for the specified 
system context.  If component?s nominal capacities differ by a large degree, then the 
specified system context will be far from the solution point, if a solution for the specified 
system exists.  The sum of the nominal charge contained in all components gives an 
indication as to the appropriateness of the specified system charge. 
A2.3.5 The Solution Phase 
The solution phase of the program begins when the user tells VapCyc to begin 
solving the currently configured cycle.  At this point, the VapCyc GUI forms the JCP and 
provides the connectivity information to the calculation engine.  The system context is 
then used to provide the calculation engine with the high and low pressure and enthalpy 
as well as the nominal system charge values to scale the Junction Solver independent 
variables and residual equations. 
A2.3.6 The Post Processor Phase 
The Post Processor Phase is the final communication between the calculation engine 
and VapCyc GUI.  The calculation engine returns a value to the GUI communicating 
information about the success of the solver.  The results are on of the following: 
 
165 
? The numeric solver in the calculation engine solved resolved to a specified 
tolerance. 
? The numeric solver in the calculation engine solved a spurious convergence, 
or local minimum, but not to a specified tolerance. 
? The numeric solver had solution problems and the returned solution is invalid. 
At this point, the GIU queries the calculation engine for solution variables, including: 
evaporator cooling capacity, condenser heat of rejection, power consumed by each 
component, system COP, evaporator superheat, condenser subcooling, and compressor 
mass flow rate. The values of several components variables are available by hovering the 
mouse over a given component.  Figure A2- 6 shows the screen after solution, and 
illustrates the presentation of compressor variables. 
 
Figure A2- 6: VapCyc Main Screen After Solution, Showing Compressor Variables