VARIABLE-DENSITY OSCILLATORY FLOW AND HEAT TRANSFER
IN 2D MODELS OF A STIRLING MICROREFRIGERATOR
By
GUOBAO GUO
A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL
OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILMENT
OF THE REQUIREMENTS FOR THE DEGREE OF
DOCTOR OF PHILOSOPHY
UNIVERSITY OF FLORIDA
ACKNOWLEDGMENTS
I would like to thank my supervisor, Professor Ulrich H. Kurzweg, for his
guidance, constant encouragement and patience. I have greatly benefited not only from
his teaching in engineering but also from his integrity, which has had a positive influence
on my life. I sincerely thank Dr. Chen-Chi Hsu for giving me the opportunity to join
the graduate program here and for teaching me many of the fundamental, yet important,
aspects of CFD during my first phase of graduate course work. Likewise, I am very
grateful to Dr. Wei Shyy for our discussions and his fruitful suggestions in my work.
I also express my appreciation to my other supervisory committee members, Dr. David
M. Mikolaitis and Dr. Kermit N. Sigmon, for their interest and suggestions.
My special thanks go to my fellow graduate students of the department here for
their help and companionship. I especially thank Dr. Ming-Hsiung Chen for his warm
friendship, help and encouragement over the years.
I am thankful to the Department of Aerospace Engineering, Mechanics and
Engineering Science for providing financial assistance. I extend my special thanks to Dr.
Michael A. DeLorenzo of the Dairy and Poultry Science Department and Dr. Ranganatha
Narayanan of the Chemical Engineering Department for offering me part-time jobs
during my study, which exposed me to some other interesting engineering and
management problems. They have certainly enriched my education.
ii
I am truly grateful to my wife Ying for her love, patience and encouragement
during this work, which were more than I could ever have requested.
The Pittsburgh Supercomputing Center supported this work through the allocation
of 183 Service Units on the Cray C90 supercomputer for three and half years, and is
acknowledged here.
TABLE OF CONTENTS
cage
ACKNOWLEDGMENTS ........................................ ............... ii
LIST OF TABLES ...................................... ..... ............... vi
LIST OF FIGURES ........................... .... ...................... vii
ABSTRACT ........................................ ........................... ix
CHAPTERS
1 INTRODUCTION ........................................ ............. 1
1.1 Background of Stirling Refrigerators ................................... 1
1.2 Literature Review ............................ ............................ 13
1.3 2D Computational Models and Scope of the Present Work .......... 17
2 GOVERNING EQUATIONS AND
PRESSURE-SPLITTING TECHNIQUE ................................... 22
2.1 Governing Equations ....................................................... 22
2.2 Separation of Hydrodynamic and Thermodynamic Pressures ........ 25
2.3 Boundary Conditions ....................................................... 28
2.4 Nondimensionalization of Governing Equations and Boundary
Conditions ....................................... ........ ...... 30
2.5 General Form of the Governing Equations ............................. 34
2.6 Dimensional Analysis ................... .... ..... ............. 35
3 NUMERICAL TECHNIQUES AND CODE VALIDATION ............. 37
3.1 Numerical Grid Generation ................... .......................... 37
3.2 Discretization of the Governing Equations ............................. 40
3.3 Solution Procedures ............................. .................... 49
3.4 Initial and Slow-Start Boundary Condations ........................... 54
3.5 Code Validation ....................................... .............. 57
4 RESULTS AND DISCUSSION ................... ............................ 65
4.1 Velocity and Thermodynamic Fields for a Typical Case ............ 68
4.2 Model Comparison ........................... ................ 82
4.3 Comparison of Using Helium and Air ................................. 85
4.4 Phase Angle Effect ......................................................... 86
4.5 Impact of Varying Refrigeration Temperatures on Performances ... 87
4.6 Changing Channel Length ................................................ 89
4.7 Effects of Changing the Ratio of Chamber Width to Channel Width 90
4.8 Discussion on Regeneration ............................................ 91
5 CONCLUDING REMARKS.................................................. 95
REFERENCES............................. .................................... 98
BIOGRAPHICALSKETCH...................................... 104
LIST OF TABLES
Table page
Table 2.1 Terms in the governing equation (2.34) .............................. 35
Table 4.1 List of the cases investigated ....................... .............. 66
Table 4.2 Comparison of Models Al and A2 ..................................... 83
Table 4.3 Comparison of Models A2 and B2 ...................................... 84
Table 4.4 Comparison of helium and air .......................................... 85
Table 4.5 Comparison of cases with different channel lengths ................ 89
LIST OF FIGURES
Figure page
Fig. 1.1 Stirling cycle and a three-component Stirling refrigerator
configuration ................................. .............. ... 2
Fig. 1.2 Schematic of a five-component Stirling refrigerator ............... 5
Fig. 1.3 Single-channel models. (a) Model Al; (b) Model A2 ............ 18
Fig. 1.4 Multichannel models. (a) Model Bl; (b) Model B2 .............. 19
Fig. 2.1 Illustrative ID discretization .................... .................... 26
Fig. 3.1 Two typical grids. (a) tt=T/4; (b) wt=w ........................... 39
Fig. 3.2 1D grids ....................................................... 39
Fig. 3.3 Control volumes for scalar quantities ................................. 42
Fig. 3.4 Staggered u and v control volumes .................................... 42
Fig. 3.5 Schematic of a general cell for 4 at point P ......................... 43
Fig. 3.6 A south boundary cell with Dirichlet boundary condition .......... 48
Fig. 3.7 A south boundary cell with Neumann boundary condition ......... 48
Fig. 3.8 Solution procedures for variable-density time-periodic flow and
heat transfer problems .................... .. .................... 55
Fig. 3.9 A piston-cylinder set-up .................. ................... 58
Fig. 3.10 Computed hydrodynamic pressure distributions along the x axis
for the piston-cylinder set-up .......................................... 58
Fig. 3.11 Thermodynamic pressure in the third cycle for the piston-
cylinder set-up ........................................ ............... 61
Fig. 3.12 Temperature distributions along the x axis for the piston-
cylinder set-up ........................................ ............. 61
Fig. 3.13 P-V diagram of the second test case .................................. 63
Fig. 4.1 Transient processes of pressures.
(a) Thermodynamic pressure; (b) Pressure drop between pistons. 69
Fig. 4.2 Convergence history to time-periodic state ........................... 71
Fig. 4.3 Velocity fields of Case 1
(a) wt=0; (b) wt=r/2; ........... ...... .................. 72
(c) wt=r; (d) wt=3r/2. ................ ................... 73
Fig. 4.4 Hydrodynamic pressure contours of Case 1
(a) ct=0; (b) wt=7/2; ..................... .... ................ 74
(c) ot=ir; (d) wt=3ir/2. ................. ................... 75
Fig. 4.5 Temperature contours of Case 1
(a) wt=0; (b) wt= r/2; .............................. ... ....... 77
(c) Wt= r; (d) ct=3r/2. ............................................. 78
Fig. 4.6 Density contours of Case 1
(a) wt=0; (b) wt=r/2; ................. ................... 79
(c) ct= r; (d) wt=3r/2. ............................................. 80
Fig. 4.7 P-V diagram of Case 1 ........... ............. .................. 81
Fig. 4.8 Phase angle effects ............................... .............. 87
Fig. 4.9 Refrigeration temperature effects ....................................... 88
Fig. 4.10 Effects of changing the ratio of chamber width to channel width 91
Fig. 4.11 Temperature distributions along x axis for two different s's
(a) # = 0.0298; (b) 0 = 0.298 ....................................... 94
Abstract of Dissertation Presented to the Graduate School
of the University of Florida in Partial Fulfilment of the
Requirements of the Degree of Doctor of Philosophy
VARIABLE-DENSITY OSCILLATORY FLOW AND HEAT TRANSFER
IN 2D MODELS OF A STIRLING MICROREFRIGERATOR
By
Guobao Guo
May 1996
Chairperson: Dr. Ulrich H. Kurzweg
Major Department: Aerospace Engineering, Mechanics and Engineering Science
Two-dimensional models of three-component Stirling microrefrigerators are
examined. The models consist of heater and cooler ended with pistons oscillating out of
phase and connected by narrow channels with a linear wall-temperature distribution,
simulating the regenerator. Navier-Stokes equations governing the variable-density,
oscillatory flow and heat transfer within two-dimensional domains with moving
boundaries are numerically solved by using the SIMPLE algorithm incorporated into a
moving grid technique and a pressure-splitting technique that is newly developed in this
study. Such a solution methodology is accurate and efficient for the complex flow and
heat transfer problems in Stirling-type devices. The pressure-splitting technique separates
the hydrodynamic and thermodynamic pressures, which are different by orders of
magnitude for gaseous flows with heat transfer like those in reciprocating machines.
Separation of the pressures increases the computing accuracy and efficiency substantially
(with saving on memory by 50% and on CPU time by 30%), and the technique can be
easily incorporated to any pressure-based numerical algorithms. The model's cooling
performance is sensitive to a dimensionless parameter that is the product of Prandtl
number, Womersley number squared, and two dimensionless geometric parameters. This
parameter should be kept small enough for the regenerator to function well. As the
refrigeration temperature decreases, the cooling capacity decreases linearly, but the
power input increases linearly, and consequently the performance degrades. It is found
that for a set of given operation conditions, there exists a lowest temperature below
which no cooling becomes possible. This is caused by the increasing cooling loss
through the channel as the refrigeration temperature decreases. Lengthening the channels
decreases such loss, but meanwhile increases the dead space, which in turn reduces the
cooling made by the expansion. The net effect depends on these two competing factors
under specific conditions. The phase angle affects the cooling capacity and efficiency,
and the optimum phase angle for either cooling capacity or efficiency is near 900. For
the multichannel model, the ratio of the chamber width to the channel width also affects
the cooling performance, and the optimum ratio depends on specific operation conditions.
CHAPTER 1
INTRODUCTION
1.1 Background of Stirling Refrigerators
1.1.1 The Stirling Principle
With the ban on the production of chlorofluorocarbon (CFC) refrigerants in the
U.S. after January 1, 1996 (and eventually the ban on the use of them worldwide),
alternative cooling devices such as Stirling refrigerators, pulse tube refrigerators and
thermoacoustic refrigerators have received more and more attention (Chase 1995). The
Stirling refrigerator is named after Scottish clergyman Robert Stirling who invented the
Stirling engine in 1816. A Stirling-type device uses a common gas (not ozone-depleting
CFCs), usually air or helium, as the working medium and works in a closed cycle. By
compressing and expanding the gas in two separated chambers at different temperatures,
either refrigeration is obtained by doing work on the gas (as in a refrigerator), or
mechanical power is produced by consuming heat (as in an engine). Figure 1.1
illustrates the working principle of a Stirling cooling device. The device consists of three
basic heat-exchanging components: the cold chamber, the hot chamber and the
regenerator, so it is called a three-component configuration. As shown in the figure, the
chambers are connected by a regenerator. In each chamber, there is a piston oscillating
out of phase with the piston in the other chamber. The regenerator is made of metallic
2
P
Rejected
Th
3 stored
rtored
Qabsobed T 4
V2 V, V
Cold Chamber Regenerator Hot Chamber
R 1
SCompression of gas at T,;
Heat rejected out.
R 2
SGas moves through R at V2
and is precooled.
SExpansion of gas at T.;
-I I Gas absorbs heat.
RR 4
SGas moves through R at V,
Sand is preheated.
R TUL__E
Fig. 1.1 Stirling cycle and a three-component Stirling refrigerator configuration
porous materials, so it has an extremely high thermal capacity compared with the
working gas and meanwhile has excellent thermal contact with the gas in it. It functions
as a heat sponge, periodically absorbing and releasing heat from and to the gas flowing
through it. The Stirling cycle consists of four idealized processes indicated in Fig. 1.1.
Process 1-2--Isothermal Compression:
All the gas is located in the hot chamber at temperature Th and the right piston
compresses the gas. As a result, heat is rejected from the gas through the chamber wall.
Process 2-3--Constant Volume Precooling:
The two pistons move in unison to keep the total volume constant and to move
the hot gas from the hot chamber to the cold chamber through the regenerator. During
this process, heat is stored in the regenerator and the gas is precooled.
Process 3-4--Isothermal Expansion:
All the gas is in the cold chamber at temperature T, and the left piston expands.
Owing to the expansion, heat is absorbed by the gas through the chamber wall. Thus,
refrigeration is achieved.
Process 4-1--Constant Volume Preheating:
The pistons move in unison so the total volume is constant and the cold gas is
moved from the cold chamber to the hot chamber through the regenerator. The
previously stored heat in the regenerator is now released back to the cold gas. Then the
closed system is ready for the next cycle.
The above discussion is based on highly idealized conditions under which the
piston motions are discontinuous as indicated in Fig. 1.1, and all the gas is in the hot
4
chamber during the compression and all is in the cold chamber during expansion. The
seidealized conditions also imply perfect heat transfer and infinite heat capacity of the
regenerator so that gas temperature in either end chamber is maintained constant and the
compression or expansion is isothermal, and the amount of heat absorbed by the
regenerator is exactly the same as that restored. It is not so in reality. First, since the
discontinuous piston motions are mechanically difficult to implement, the pistons in a real
Stirling cooling device usually oscillate sinusoidally such that the hot-chamber volume
variation leads the cold-chamber volume variation by a phase angle of about 90.
Second, it is not possible for all the gas to be in either end chamber for compression or
expansion because there is always some gas in the regenerator or in other pockets (the
volume that can not be swept by pistons is called dead space). And finally, both the
thermal capacity of the regenerator and the rate of heat transfer between the gas and the
three exchangers are finite. Therefore, the gas in the chambers does not exactly have
the wall temperature and the regenerator may not function perfectly either. As a result
of these complicating factors, the thermodynamic cycle is not precisely composed of two
isothermal and two equal-volume processes, and the actual P-V diagram will be
somewhat different from the one shown in Fig. 1.1 (Chapter 4 will show a realistic P-V
diagram).
In the three-component configuration, the expansion or compression takes place
in the big end chambers, so there may not be enough time and area for heat exchange
between the gas and the chamber walls. For better performance, usually the compression
and expansion chambers are separated from the heating area (heater) and the cooling area
Compression Right
Fig. 1.2 Schematic of a five-component Stirling refrigerator
(cooler), respectively. The heater or cooler is made up of many metal tubes of small
diameters such that the heat transfer is adequate. Figure 1.2 demonstrates a five-
component configuration.
1.1.2 Brief History of Stirling Machine Development and Applications
Details on the development history of Stirling-type devices can be found in books
by Walker (1973), Reader and Hooper (1983), West (1986) and Organ (1992). Since
Stirling's invention almost two centuries ago, many different Stirling devices have been
developed and found their various potential and actual applications: automotive engines,
small electric generators, underwater power units, solar thermal conversion, space
power, artificial heart power, heat pumping, cryogenic cooling engines (cryocoolers),
near-ambient-temperature refrigerating machines (refrigerators). The early and primary
developments focused only on machines providing mechanical driving power (i.e., prime
movers). After the middle of the nineteenth century, because of the inventions of the
internal combustion engine and electric motor, the Stirling prime movers gradually lost
their commercial market to these new inventions. However, the appealing characteristics
of the Stirling engine still attracted engineers to continue research and development.
When working as a prime mover, a Stirling engine can use essentially all kinds of heat
sources: solid, liquid or gaseous fossil fuel combustion, solar heat, nuclear heat, etc.
The combustion with a Stirling engine is external and the process can be easily
controlled, producing a lower air pollution level than an internal combustion engine. The
Stirling engine does not require any valves to control the flow or pressure, and the
combustion (if needed) does not result in periodic explosions. So the engine is very quiet
and the machine vibration is very small in comparison with other engines. When
working as a refrigerator a Stirling refrigerator uses no CFCs, so it is environmentally
safe.
One of the leading institutions for the Stirling devices is the Philips Research
Laboratories in Eindhoven, Holland. Back to the late 1930s, Philips developed portable
Stirling electric generators for radio sets. Later they studied and applied new
technologies to the then existing Stirling engines. By the mid-1940s, they had increased
the Stirling engine's power per unit weight by a factor of 50, increased its size per unit
power by a factor of 125 and improved its efficiency by a factor of 15 (Reader and
Hooper 1983). It is Philips who brought the Stirling engine into its modern era.
From the late 1950s, Philips along with several European and American
companies (including Ford and GM) began to develop Stirling automotive engines.
7
When the oil crisis hit the world in the 1970s, the interest in such engines increased
suddenly because of the wide variety of energy (heat) sources that could be used by
Stirling engines. In 1978, sponsored by the U.S. Department of Energy, several U.S.
companies including the AM General Corporation, Mechanical Technology Inc. (MTI)
and the United Stirling carried out intensive researches in hopes of replacing the internal
combustion diesel and gasoline engines in automobiles. The main difficulty with Stirling
automotive engines becoming commercialized is that they cannot compete with existing
diesel and gasoline engines in manufacture and maintenance costs because the latter have
already been well developed and are under continuous improvement. As a result, all
Stirling automotive programs stopped about a decade ago.
One of the active areas in the Stirling research and development is the artificial
heart. Since the 1960s the U.S. National Institutes of Health and some companies (e.g.,
Stirling Technology Company) have been developing totally implantable artificial heart
devices (West 1986, White 1985). Its merits of quietness, very small vibration and high
reliability make the Stirling engine the best power provider for this purpose. A small
Stirling engine in such devices uses radioisotopes as the heat source for a long and
continuous life.
Since the 1960s the Stirling engine has been identified as the key to the future
space energy conversion systems by NASA Lewis Research Center (LeRC) (Shaltens and
Schreiber 1989). The heat sources for the systems can be nuclear reactor heat,
radioisotopic heat or solar heat. Stirling systems have the potential of meeting space
power's requirements: high power capacity, high thermal efficiency, high reliability, low
8
vibration and very long life. NASA initiated the free-piston Stirling space power
program in the mid-1980s, and a research engine was built by MTI under contract with
NASA (Tew 1988). Some main design goals are: output electric power 12.5 kW,
thermal efficiency 25%, life 60,000 hours (seven years), hot-end temperature 1050 K,
cold-end temperature 525 K, frequency 70 Hz and mean pressure 150 atms. Initial test
results with the originally built engine indicated that the engine performance did not
reach the design goal (Tew 1988). NASA LeRC, MTI and several other U.S. companies
have done intense experimental and analytical work in the past several years (Cairelli et
al. 1989-1991 & 1993, Dochat 1990, Jones 1990, Dochat and Dhar 1991, Thieme and
Swec 1992, Wong et al. 1992), and some significant improvements on the engine have
been made and novel component technologies have been obtained.
The commercialization of Stirling solar energy conversion systems has been
pursued for almost two decades in several countries, e.g., Germany (Keck et al. 1990),
Japan (Isshiki et al. 1991), Russian (Loktinonov 1991), Sweden (Diver et al. 1990) and
the U.S. (Andraka et al. 1994). In such systems the solar heat is collected by a parabolic
dish and used as the heat source of a Stirling engine, then the engine can drive a water
pump for irrigation or drive an electrical generator for electricity supply. The Stirling
solar system is more efficient than conversion systems using photovoltaic cells. In the
early solar Stirling engines, collected solar rays were aimed directly at the engine heater
tubes. This produced very hot spots in the tubes and so the temperature distribution
there was not uniform, resulting in a too restrictive requirement on materials and thus
limiting the engine life and performance. To solve this problem, a liquid-metal reflux
receiver is used to transport the collected solar energy to the heater tubes. In the reflux
receiver, a liquid metal (usually sodium or sodium-potassium alloy NaK-78) evaporates
at the solar absorber due to the very high temperature there, and then condenses at the
engine heater tubes. The latent heat released by the condensing vapour is utilized to heat
the tubes. Most of the current technology advances are made in the U.S., primarily due
to the Department of Energy's Solar Thermal Program and the studies made by some
private companies.
The most successful commercial application of the Stirling engine is cooling.
Seventeen years after Stirling's invention, John Herschel in 1834 realized that the Stirling
machine could also work as a refrigerator, and the first Stirling refrigerator was built by
Gorrie in 1845 (Kohler 1960). A review on the applications of Stirling refrigerators was
given by Walker et al. (1988). In 1954 Philips first commercialized a Stirling cryocooler
with a cooling capacity of 1 kW at 80 K for air liquefaction. Since then Stirling
refrigerators have been commercially developed for various applications with a cooling
capacity from a few milliwatts to several kilowatts and with refrigeration temperatures
from as low as a few degrees Kelvin to near ambient. At very low temperatures, usually
below 160 K, the efficiency of a Stirling refrigerator is much higher than all other
cooling devices such as the ones using the Joule-Thomson effect, or the ones using the
evaporation of a liquid (i.e., the commonly used refrigerators). So Stirling refrigerators
have dominated the world market of gas liquefiers for half a century. Liquified gases
are mainly used for cryogenic experiments, in which the temperature is usually in the
cryogenic range, namely below 120 K. With the findings of more and more high-
10
temperature superconductors and their increasing applications in electronics, the current
need for Stirling cryocoolers is more demanding than ever.
Beside the large cryocoolers for gas liquefaction, small and miniature Stirling
cryocoolers have also been widely used to cool some special small sensors and electronic
components that need cryogenic cooling for higher resolution and sensitivity. These
include infrared detectors, X-ray and y-ray solid state-detectors, etc. Increasing numbers
of such detectors are carried by airplanes, space shuttles, space stations, satellites and
deep space probes for civilian and military purposes (e.g., Earth observation,
atmospheric science, astronomy and space surveillance). Since the 1960s, many efforts
have been made in the U.S. to develop spaceborne cryogenic systems (Sherman 1982).
The required cooling power for the systems is usually between a few milliwatts and tens
of watts at temperatures below 80 K. Due to the restriction on the size and weight of
the loads in space missions (especially the long-term ones), it is not feasible to carry a
large amount of solid or liquid cryogens with spacecraft. Therefore miniature
mechanical cryocoolers (of the Stirling, Joule-Thompson absorption, pulse tube, acoustic,
turbo Brayton and magnetic types) have been pursued because of their small size and
lightweight, mechanical quietness, and expected very long life without maintenance. As
reported by Chan et al. (1990), of all the eight known space flights (between 1970 and
1989) using mechanical cryocoolers, Stirling cryocoolers were flown six times. The lead
of the Stirling cryocooler over others is attributed to the mature technologies for the
large-size Stirling cryocoolers as mentioned before and to its higher efficiency than others
in the temperature range of interest.
In search of alternatives to the existing CFC-containing vapour-compression
refrigerators exclusively used at near ambient temperatures, the Stirling refrigerator has
recently been considered as a promising contender (Walker et al. 1992, Berchowitz 1992,
Chase 1995). According to Walker et al., at near ambient temperatures, the coefficient
of performance of a vapour-compression refrigerator is about 60-70% of the Carnot
value, whereas that of the Stirling refrigerator is only about 45% of the Carnot value,
and that is why the former has long dominated the refrigeration market at near ambient
temperatures. Researchers in many countries are now developing Stirling refrigerators
for domestic refrigeration, air conditioning, electronic cooling and so on. In the late
1980s, the NASA Manned Spacecraft Center at Houston, Texas, sponsored a company
to build a Stirling refrigerator for manned spacecraft. Air was used as the working fluid
in the device, which makes it safe for astronauts in case of refrigerant leakage and simple
for gas replenishment since the working gas is just the air in the spacecraft. To replace
the existing refrigerator or freezer used in the space shuttle mid-deck for the preservation
of experiment samples, primarily blood and urine, NASA's Life Sciences Projects
Division at the Johnson Space Center in 1992 also initiated a project to design and
develop a Stirling refrigerator (McDonald et al. 1994). The Stirling refrigerator was
designed by Sunpower Inc. of Athens, Ohio. This unit used helium as the working fluid,
and was lightweight (101 Ibs, including the heat pipes), quiet and efficient (heat lift of
90.6 watts at 4C with a power input of 60 watts, and heat lift of 88.5 watts at -22C
with power input of 70 watts). It was successfully flown on board a space shuttle in
February 1994 for three days as a freezer and four days as a refrigerator. Sunpower Inc.
has done pioneer work in developing free-piston Stirling refrigerators for intermediate
and near ambient temperatures and in commercializing them (Berchowitz 1992). The
company has developed two prototypes of Stirling refrigerators: "Super Fridge" and
"Cool Junior." The "Super Fringe" is aimed at replacing the existing domestic vapour-
compression refrigerators, and is capable of lifting 250 watts at -26C with an input of
215 watts. The "Cool Junior" is very small, and so is ideal for cooling electronics,
vaccines, etc. It has a cooling capacity of 35 watts at -50C and 54 watts at -26"C.
Work on the development of Stirling refrigerators for domestic applications has also
recently been done by Gold-Start Co. of Korea (Kim et al. 1993) and by several Japanese
companies and institutions like Toshiba (Otaka et al. 1993, Isshiki et al. 1994). Creare
Inc. of Hanover, New Hampshire, has developed a proof-of-principle diaphragm defined
Stirling refrigerator for near ambient applications (Shimko et al. 1994). The design
innovation in this machine is the application of diaphragms to sweep and seal the
compression and expansion volumes. That completely avoids the serious sealing and
lubricating problems associated with the pistons in the conventional design.
1.1.3 The Challenges in Analyzing and Simulating Stirling-Type Devices
In reciprocating machines, such as internal combustion engines, Stirling type-
devices (prime mover, refrigerators and heat pumps) and pulse tube refrigerators, the
flow and heat transfer are time periodic, multidimensional with moving boundaries and
with abrupt cross-section changes. Besides, owing to large temporal volumetric
oscillations and spatial temperature gradients, density varies significantly in time and in
space, and pressure level also oscillates to a large degree. Moreover, periodic transition
to turbulence may also occur for large enough machine dimensions or at high enough
operation frequencies, and the time-periodic flow and heat transfer in the porous medium
of the regenerator complicates the problem even further. Therefore analysing or
numerically simulating such flow and heat transfer remains a challenging task. In the
design and analysis of such machines, usually one-dimensional (ID) models are used
along with empirical formulas for gas-wall heat transfer coefficients, skin friction
coefficients, pressure drops and so on (Reader and Hooper 1983, West 1986). These
formulas, however, are correlated only from the experimental data of incompressible
(constant-density) steady unidirectional flows; and it has been shown that the formulas
and consequently the 1D models are not able to accurately predict the performance of
Stirling devices (Geng and Tew 1992, Smith et al. 1992, Cairelli et al. 1993), although
they meet such design needs as simplicity and short design cycle.
1.2 Literature Review
A classical thermodynamic analysis method for the Stirling engine was first
presented by G. Schmidt in 1861 (Walker 1973). The differences between the Schmidt
theory and the highly idealized Stirling cycle are that the Schmidt theory allows for the
more realistic sinusoidal piston motions, and that it also takes the dead space into
consideration. This theory, however, still retains some of the highly idealized Stirling-
cycle assumptions. For example, the gas temperature distribution in each heat exchanger
is uniform, therefore the thermodynamic process isothermal; the regeneration is perfect;
14
pressure drops across engine components are negligible and thus the instantaneous
pressure is the same through out the system. Although the method can give closed form
solutions based on these assumptions, it is still of very limited use. This is because some
assumptions it uses are still far from the real situation, and so the predicated engine
characteristic values must be multiplied by a factor of 0.3 to 0.4 in order for them to
match those of a real engine.
The Schmidt method was later replaced by newer analytical methods like linear
harmonic analysis, which was introduced to the Stirling community after 1960s. West
(1984) thoroughly reviewed the linear harmonic analysis method. This method regards
all the major component-averaged quantities (volume, pressure, temperature, etc.) as
unknown sinusoidal functions that are linear combinations of sin(wt) and cos(wt). The
governing equations used in the method are algebraic and ordinary differential equations
(ODEs) written for each individual heat exchanger with the component-averaged
quantities as the dependent variables. Examples of the algebraic equations are the
equation of state, and the equation of continuity, namely the summation of the mass in
all the three (sometimes five) components equals constant. The sinusoidal functions are
substituted into the governing equations; the resulting nonlinear terms of sin(wt) and
cos(wt) are linearized by using binominal expansion first and then simply neglecting
higher order terms. Another way of treating the nonlinearity is to express the nonlinear
terms of the major quantities themselves in the equations as Fourier series truncated after
the sine and cosine terms. In either way, the linearized equations are simply a few (at
most less than dozens) simultaneous algebraic equations. The latter approach is able to
account for the effects of the non-isothermal or transient heat transfer process and
pressure drops, which the Schmidt analysis omits. Huang (1992) reported a one-
dimensional computer code based on such harmonic analysis method.
In the past decade, fundamental problems of oscillating internal flow and heat
transfer with Stirling applications have been experimentally studied to improve analysis
and design accuracy. Tang and Cheng (1993) carried out experiments on incompressible
periodically reversing pipe flows, and correlated a cycle-averaged Nusselt number in
terms of the Reynolds number, dynamic Reynolds number and the nondimensional fluid
displacement. Smith and co-workers (Kormhauser and Smith 1989, Smith et al. 1992)
did experiments on an apparatus including a piston-cylinder space connected to a
dead-end heat exchanger space, and proposed a "complex Nusselt number" model in
which heat flux consists of one part proportional to the gas-wall temperature difference
plus a second part proportional to the rate of change of the difference. Tew and Geng
(1992) gave a review on related experimental projects supported by NASA.
Several researchers have recently begun using two-dimensional (2D) models to
simulate Stirling devices numerically. Gedeon (1989) developed a 2D computer code
modelling the unsteady laminar flow and heat transfer in Stirling components. He used
Darcy's law to simulate the flow and heat transfer in the regenerator, and applied the
density-based Beam-Warming scheme (Beam and Warming 1978) as the solution
algorithm. This scheme was originally designed for solving compressible (high Mach
number) flows, but the Mach number of the flows in a Stirling device is quite low,
typically below 0.1. (As will be seen later, the Mach number for the micro-configuration
16
in the present work is below lx10".) Density-based numerical schemes, e.g., the Beam-
Warming schemes, are inefficient and inaccurate in solving low Mach number or
incompressible flows (Karki 1986), and Gedeon reported similar problems encountered
in his applications. Ibrahim and co-workers (1991, 1992) also used 2D laminar models
to study Stirling components numerically, and suggested a way to correlate the gas-wall
heat transfer coefficient under oscillating flow conditions, and they investigated the
effects of sudden cross-section changes on heat transfer and friction coefficients. The
solution method employed in Ibrahim and co-workers' work was the pressure-based
algorithm, SIMPLE, which was initially developed by Patankar (1980) for incompressible
flows. By utilizing a 2D model, Kurzweg and Zhao (1993) numerically examined the
velocity and pressure fields only within one end chamber of a Stirling microrefrigerator
with various orifice configurations, and they recommended a better arrangement of the
orifices. They also used the SIMPLE algorithm. Density in the chamber was regarded
by them as a function of time only.
In all the above-mentioned 2D numerical studies, the heat-exchanger components
were separately considered as open-ended channels (at least at one end), and the velocity
boundary conditions (B.C.) at the component entrance or exit were also simply
approximated by slug flow and the temperature and density there were taken as uniform.
Therefore, the effects of large volume and pressure oscillations were not properly taken
into account.
17
1.3 2D Computational Models and Scope of the Present Work
The main purposes of this work are 1) to propose 2D models for Stirling
microrefrigerators, 2) to give a systematic solution methodology suitable for solving the
variable-density time-periodic laminar flow and heat transfer within confined regions with
moving boundaries, and 3) to numerically investigate the flow and heat transfer
phenomena in such configurations.
For computational simplicity, only the three-component configuration is
considered in this study. The model adopted here integrates all the three heat exchangers
(cooler, heater and regenerator) together, allowing the effect of large pressure and
volume variations to be considered properly. To reduce the computational complexity,
the regenerator is replaced by narrow channels with linear wall temperature distributions.
Since the channels are narrow, there is enough time for the gas flowing in them to adjust
its temperature to the desirable wall temperature. Two types of models are studied:
single-channel models Al (Fig. 1.3 (a)) and A2 (Fig. 1.3 (b)); and multichannel models
Bl (Fig. 1.4 (a)) and B2 (Fig. 1.4 (b)). In the single-channel models (Al and A2), the
two chambers are connected only by one channel, and all the chamber walls are
isothermal. In the multichannel models (B1 and B2), the two chambers are connected
by more than one channel. In Models A and B, the piston surface is regarded as
insulating; and the wall temperature of the left chamber is set at T, and that of the right
is at Th. The isothermal surfaces are extended into part of the connecting channels to
a length "d" in Models A2 and B2. As indicated in the figures, the channel half width
is "a", the total channel length is "L", the half chamber width is "b". In either the left
(a) Single-channel model: Al
XI.ft L- L 1 xI right
xen = -L/2 c + Esin((ot) x,g,, = L/2 + c + Esin(ot+(po)
(b) Single-channel model: A2
xl. -- L ---- Xgt
x,ia = -L/2 c + esin(ot) xh = L/2 + c + esin(cot+()
Fig. 1.3 Single-channel models. (a) Model Al; (b) Model A2
x4
uos!d 6u!lnsul
H -el
uols!d Bu!lelnsul
0-
E
a-
C-
C:
2C
I I
I 1
F-
E
0
'U
-
m I
Eg
.o, E
o
20
or the right chamber, a piston oscillates sinusoidally with an amplitude of "e" and an
angular frequency of "o". The piston motion in the right (hot) chamber leads that in the
left by a phase angle "po" for cooling. The models assume a small clearance of "c-E"
between the piston and the vertical wall. The reasons for this are 1) to avoid difficulties
in numerical computations if there is no clearance and 2) to model the actual situation
where there is always a clearance and some dead space. Since the models are intended
to simulate microrefrigerators, the models are in the centimeter range.
Even though the flow within the 2D models is laminar because of the very small
sizes used, numerically solving the pertaining heat and transfer problem still constitutes
a difficult task--the challenging elements mentioned before remain the same. Besides,
in the variable-density flows with heat transfer, there are significant volumetric changes.
The associated thermodynamic pressure, therefore, is many orders larger than the
hydrodynamic pressure, and that could result in low computational accuracy and
efficiency, and even physically unrealistic solutions. To overcome this problem, a
pressure-splitting technique is developed to separate the thermodynamic and
hydrodynamic pressures in this study. The Navier-Stokes and energy equations
governing the variable-density, oscillatory flow and heat transfer are solved numerically
by using the pressure-based algorithm, SIMPLE (Patankar 1980), incorporated with the
pressure-splitting technique. Since the domain of interest is bounded by moving pistons,
a moving grid technique as used by Demirdzic and PeriC (1990) is also incorporated into
the algorithm. This study proves that such a methodology is accurate and efficient in
solving flow and heat transfer problems like those occurring in Stirling type devices.
21
By using a code based on the above method, the 2D models of Stirling
microrefrigerators are studied. It is found that the model's regeneration and cooling
performance are very sensitive to a dimensionless parameter that is the product of the
Womersley number squared, Prandtl number, and two geometric parameters.
(Womersley number is defined as c = a ,/wl/v, where w is the oscillation angular
frequency and v the kinetic viscosity.) This dimensionless parameter should be kept
small enough to let the gas passing through the channel adequately adjust the temperature
to that of the wall's. A parametric study of the phase angle change is also made and the
optimum phase angle is shown to be close to 900. The effects of varying the temperature
ratios between expansion and compression chambers are also examined in this study.
The cooling capacity is found to vary linearly with the temperature ratio. Other issues
investigated in this study are the effects of changing the channel length and changing the
ratio of chamber width to channel width. Some of the findings from the current
numerical study agree well with experimental data with actual Stirling refrigerators, and
that suggests the proposed 2D models are valuable in helping gain insights into the
Stirling microrefrigerators.
The arrangement of this dissertation is as follows. The governing equations and
the pressure-splitting technique will be discussed in Chapter 2. Chapter 3 will describe
in detail the numerical techniques solving these equations, and the validation of the
solution methodology will also be made in Chapter 3. Chapter 4 will focus on the
numerical results and discussions, and finally Chapter 5 contains the concluding remarks.
CHAPTER 2
GOVERNING EQUATIONS AND PRESSURE-SPLITTING TECHNIQUE
In this chapter, the governing equations of fluid mechanics and heat transfer
appropriate to the problem under consideration are first given. A method is devised to
split the pressure into hydrodynamic and thermodynamic parts so that the computational
efficiency and accuracy for incompressible (very low Mach number) but variable-density
flows with heat transfer are enhanced. The equations and their boundary conditions are
nondimensionalized. Then the governing equations are cast into a form that is suitable
for the numerical solution algorithm to be used. At the end of this chapter, dimensional
analysis is carried out for the cooling power, heating power and the rate of work input
in the Stirling microrefrigerator models.
2.1 Governing Equations
The governing equations of viscous, variable-density, unsteady flows and heat
transfer with negligible body forces in domains with deformable moving boundaries can
be written in an integral form as (Hansen 1967):
Conservation of mass,
d fpdv + fp(V9- Vb)dS = 0 (2.1)
SV(t) s(t)
Conservation of momentum,
d- fpdV + fpV(9V-V fidS = f(-pi+t)-dS, (2.2)
v(t) s(t) s(t)
Conservation of energy,
d fpEdV+ pE(9- -idS= fIdS+ 9-(-pi+T)-fdS. (2.3)
V(t) S(t) S(t) (t)
In the above equations, V(t) represents the control volume, S(t) the moving control
surface, Vb the velocity of the control surface, and fl the outward unit normal to the
control surface; p is the density of the fluid, V the velocity of the flow, E the total
energy of the fluid per unit mass, p the pressure, i the identity stress tensor, j the
deviatoric stress tensor, and i the conduction heat flux vector. Here the assumption of
Newtonian fluid is used.
Notice that if the velocity of the control surface is the same as the fluid velocity
there, then the Lagrangian approach is taken; if the velocity of the control surface is
zero, then the Eulerian approach is taken. Otherwise the viewpoint is neither Lagrangian
nor Eulerian, as happens in this study.
Since the body force is negligible, potential energy is negligible, too. The total
energy is then given as
1e +
E=e+ ,
2
where 1-.V is the kinetic energy per unit mass, and e the internal energy per unit
2
mass. The internal energy for a perfect gas is given by,
e = cT. (2.5)
Here c, is the specific heat at constant volume, and T the absolute temperature.
According to Stokes' hypothesis, the deviatoric stress tensor is
S= I[VV +(VV)'] t2-(VV)I (2.6)
where the superscript "T" represents the transposition of a tensor, and p is the molecular
viscosity of the fluid.
The equation of state gives the relationship among p, p and T, and reads
p = pRT. (2.7)
where R is the gas constant.
In the energy equation, the conduction heat flux is given by Fourier's law
[ = -k VT (2.8)
where k is the coefficient of thermal conductivity.
Viscosity is determined by using Sutherland's formula
Tr +TI T
p = T ( )3 (2.9)
T+T Tr
in which T, is some reference temperature, p, the viscosity of the fluid at the reference
temperature, and T, a constant that depends on the fluid.
Thermal conductivity k can be determined from viscosity 1 and specific heat at
constant pressure
k (2.10)
Pr
where Pr is the Prandtl number, which is usually a constant for a specific fluid. For a
perfect gas,
p c, = R, (2.11)
and CP = Y7c, (2.12)
where y is called the ratio of specific heats, and is a constant for a given gas.
2.2 Separation of Hydrodynamic and Thermodynamic Pressures
The pressure p(x,y,z,t) consists of hydrodynamic pressure, denoted by ph, and
thermodynamic pressure, denoted by p', i.e.,
p(x,y,z,t) = ph + p'. (2.13)
The hydrodynamic pressure is related to the inertial and viscous forces and is usually a
function of both time and space, so
ph = p(x,y,z,t); (2.14)
and the thermodynamic pressure is related to the volumetric changes and temperature
variations. At very low Mach numbers, the thermodynamic pressure is a function of
time only, i.e.,
p' = P(t), (2.15)
Consequently,
p O(M2) c .
p'
In numerically solving the
momentum equation (2.2),
pressure gradient terms like ap/ax
need to be approximated. As a
demonstration, consider a ID
problem where the pressure Fig.
Fig.
gradient is dp/dx. As shown in
Fig. 2.1, the pressure gradient at xi+l,, = (x, + x,+
estimated by
(2.16)
Pi P,.1
I I I x
--I--i-- --^x
X, Xi+1/2 Xi+1
2.1 Illustrative ID discretization
,)/2, denoted as (dp/dx),+,,,, is
(dp i Pi -P
dx )lr2 x.l-xi,
(2.17)
where
Pi = phi + pt, (2.
Pi+l = Phi+ + pt, (2.19)
From the above discussion, pt is independent of x and is several orders larger than ph.
So what is done in Eq. (2.17) is subtracting two very close numbers from one another
with the difference between them being many orders smaller than the two numbers
themselves. If the computer used does not have enough digits, large round-off errors and
even catastrophic cancellation may occur in the calculation and the numerical solutions
are erroneous (Golub and Ortega 1992). Although using double-precision arithmetic in
27
computations can increase the accuracy, a price of more computer memory, storage and
CPU time must be paid.
Wei Shyy suggested that effort should be made in order to find a way of
separating hydrodynamic and thermodynamic pressures so that the computing accuracy
and efficiency could be increased. The rationale follows. Substituting Eqs. (2.18) and
(2.19) into Eq. (2.17) yields
p (Pt:h h+p) h h
).2 (P Pt(Pi +P ) Pil-Pi" d 1 (2.20)
dx) ixi1 _xi xi-x dx i
Then the problem of catastrophic cancellation does not exist anymore, and therefore the
computing accuracy is enhanced without using double-precision arithmetic. The
computer CPU time, memory as well as storage are saved. Now the question is how to
determine ph and p'.
Such a method is developed in the present study. First consider how to determine
the hydrodynamic pressure ph. It can be seen from Eq. (2.20) that the pressure p
appearing in the momentum equation (2.2) can now be replaced by the hydrodynamic
pressure ph and therefore p' is solved by using the pressure-based SIMPLE algorithm,
which will be discussed in Chapter 3. It should be noted that the total pressure appearing
in the energy equation (2.3) has nothing to do with the pressures splitting procedures
because the pressure there is replaced with pRT via the equation of state (2.7).
Next consider the thermodynamic pressure p'. It can be determined by using the
'Personal communication, 1994.
equation of state and the total mass of the whole system, and the procedure follows.
From Eqs. (2.7) and (2.13), the density is given by
p ph+p t (2.21)
RT RT
The total mass of the entire system is known at any time instant, and
Total Mass = p(x,y,z,t)dV = f Phxyzt)+pt)dV (2.22)
vt) vt) RT(x,y,z, t)
Therefore,
1 h
pt(t) 1 [Total Mass dV]
1-dV t)RT (2.23)
v(t RT
In this way, the hydrodynamic and thermodynamic pressures, which are of
different physical meaning and of different order of magnitudes, are determined
separately, and consequently the computing accuracy and/or efficiency is increased.
2.3 Boundary Conditions
The symmetry of the geometry and the boundary conditions under consideration
makes the solution also symmetric about the x-axis (see Figs. 1.3 and 1.4). Numerical
calculation takes the advantage of the symmetry by considering only half the domain of
interest as the computational domain with symmetry boundary conditions being applied
at the symmetry line, as a result the computing effort needed is greatly reduced. For
example, the upper half domain of interest in the present work is taken as the
29
computational domain. In what follows, the boundary conditions for Model B2 (Fig. 1.4
(b)) are listed.
Velocity boundary conditions are as follows:
x = x,,a, 0 y y b, left piston wall, no-slip:
u = u, f = -awsin(0t), v = 0. (2.4)
xif < x < -L/2, y = b, imaginary boundary, symmetry:
au/ly =0, v=0.
x = -L/2, a < y b, vertical wall of the left chamber, no-slip:
u = 0, v = 0.
-L/2 -< x L/2, y = a, upper channel wall, no-slip:
u=0, v = 0.
x = L/2, a < y : b, vertical wall of the right chamber, no-slip:
u = 0, v = 0.
L/2 < x < xg, y = b, imaginary boundary, symmetry:
aulay =0, v=O.
x = xrh, 0 y b, right piston, no-slip:
u = uht = -eosin(wt+ .o), v = 0. (2.25)
xija < x < x,,, y = 0, symmetry:
au/ay =0, v=0.
The temperature boundary conditions are as follows:
x = xIa, 0 y s b: aT/ax = 0.
xfn < x < -L/2, y=b: aT/Oy = 0.
30
x = -L/2, a < y b: T = T,.
-L/2 x < -L/2+d, y=a: T = T,.
-L/2+d x L/2-d, y=a: T = T, + (Th T.) (x + L/2 d)/(L 2d).
L/2-d < x L L/2, y=a: T = Th.
x = L/2, a < y b: T =Th.
L/2 < x < x,,,, y=b: MT/ay = 0.
x = xgh, 0 y b: aT/ax = 0.
x,,a < x < xr, y =0: aT/ay = 0.
Here the boundary conditions for velocity and temperature only include two types,
Dirichlet and Neumann, and all the Neumann boundary conditions here have zero
gradients. Around the entire boundary, the velocity component normal to the boundary
is always known.
Pressure does not requires explicit boundary conditions here. As will be
discussed in the next chapter, the numerical algorithm adopted in this study does not need
explicit pressure boundary conditions when the velocity normal to the boundary is given.
2.4 Nondimensionalization of Governing Equations and Boundary Conditions
In the analysis or numerical computation of a physical problem, the governing
equations and boundary conditions are usually nondimensionalized. The reasons for this
are twofold: first, by doing so the nondimensional parameters governing the physical
problem are obtained; second, computation convenience is also achieved by using
nondimensional quantities. The governing equations and boundary conditions are
nondimensionalized as follows.
First reference quantities are properly chosen:
Reference length: L, = a (channel half width);
Reference time: t, = 1/w (w = angular frequency of the piston oscillation);
Reference velocity: U, = L/t, = aw. The reference velocity is defined in terms of the
reference length and reference time and cannot be arbitrarily chosen;
Reference temperature: T,. It should be chosen such that it is as close as possible to
the temporal and spatial average temperature;
Reference density: p,. Like T,, the reference density is also chosen so that it is as
close as possible to the temporal and spatial average density;
Reference viscosity: r, = i(T).
Then the nondimensional quantities, which are denoted by bars, can be defined as
Xi xi t ui T
x t = ot, u, = T
Lr a tr Ur ao Pr T,
e E 1 T
-i e=--, E= e2+- I V
1r Ur2 U 2 U,
L,
Special care should be taken in nondimensionalizing the hydrodynamic and
thermodynamic pressures because of their inherent different scales. The corresponding
nondimensional quantities are defined as
h h
P ,
prU2'
-t_ Pt
PrRTr
For convenience from now on the bars over the nondimensional quantities will be
dropped with the understanding that nondimensional quantities are used, unless otherwise
mentioned.
Finally, the above nondimensional quantities are substituted into the governing
equations and boundary conditions. The nondimensional equations are:
Conservation of mass,
d fpdv + fp(V-V0-dS = 0, (2.26)
V(t) S(t)
Conservation of momentum,
d fpVdV+fp V(V--V'-fdS = -fphfdS+ -fidS, (2.27)
v(t) s(t) s(t) s(t)
Conservation of energy,
d fpTdV + pT( -V.)ArdS = s V(PlT) fdS
to st) Pr t)
-( -1) pTV-fidS + (Y (--)M2 f -f dS (2.28)
s(t) a S()
+ y(y-1)M2[ df lpV -VdV + flpV*V(V-9-b)dS]
dt)2 S(t)2
Equation of state,
pt + yM2ph = pT. (2.29)
Here one should notice the following two points: 1) the pressure in the momentum
equation (2.27) is just the hydrodynamic pressure ph, and 2) the form of the
nondimensional equation of state (2.29).
The nondimensional parameters in these equations are:
Prandtl number: Pr = (,p /k), = /K,,
Ratio of specific heats: y = Cp /c,,
Mach number: M = aw/(yRT,)/2 = U,/(7RT,)1/2
= (reference speed) / (speed of sound at T,), (2.30)
Womersley number: a = a (Pl//i,)1/2 = a (I/r)l12. (2.31)
This parameter characterizes oscillatory flows, and is also known to be Stokes number.
Instead of the above parameter, another parameter is also used in the literature. The
latter is defined as
Va = a '/v, = a2 = L, U,/,. (2.32)
It is called the Valensi number or nondimensional frequency (Ahn and Ibrahim 1991),
or oscillating Reynolds number (Tew 1987).
For simplicity, the nondimensional velocity and temperature boundary conditions
will not be listed here. The nondimensional parameters resulting from the boundary
conditions are b/a, c/a, d/a, E/a, L/a, T, IT,, Th /T, and oP,.
Since T is used as the dependent variable to be solved for in the energy equation
(2.28), quantities p, e and consequently E are all expressed in terms of T there. The
second term, including the minus sign "-", on the right-hand side (RHS) of the equation
is the rate at which work is done to compress the gas and therefore the rate at which the
internal energy increases owing to compression; the third term on the RHS is the rate of
dissipation of mechanical energy due to viscosity; and the last terms in the square
brackets on the RHS are related to kinetic energy. In this study the dissipation term and
the terms related to kinetic energy are neglected because they all include the factor M2
and M here is always smaller than 10'. Therefore, the energy equation is simplified,
without loss of computational accuracy, to
d f pTdV+ pT(V-V0-"idS= f V(pT)fidS-(y-1) fpTV-fdS.
v(t) S(t) P ,(t) s(
(2.33)
2.5 General Form of the Governing Equations
Any of the integral equations presented in the preceding section can be cast into
the following general form
SfpdV + ) p(V-~ b)dS = (f(rV4)'fdS + f(dv, (2.34)
dvt) S(t) S(W) V()
where F is the diffusion coefficient, the source per unit volume. Each specific integral
equation corresponds to a specific set of 4, P and [, which are tabulated in Table 2.1 for
two dimensions.
The first term on the left-hand side (LHS) of Eq. (2.34) is the rate of change of
the extensive property inside the control volume, that corresponds to the intensive
property 4d; the second term on the LHS is the rate of efflux across the control surface,
and called the convection term. The first term on the RHS of the equation is called the
rate of diffusion, and the second term on the RHS the source term. Anything that cannot
be put into the three terms mentioned above is put into the source term, as is the case
in the momentum and energy equations.
Table 2.1 Terms in the governing equation (2.34)
r r
Continuity 1 0 0
ph1 a au 2a av a av
Ca ax 3ax ax 3a x y ay ax
u-Momentum u
2 -ph 1 a av 2a au, a au
ay 3ay ay a3 y ax ax ay
v-Momentum v
Pr -( )[ a(puT) a(pvT)
Energy T x ay
2.6 Dimensional Analysis
The physical quantities of much interest with the computational models are:
cooling power, heating power and rate of work input, the definitions of which will be
given in the next chapter. Since they all have the dimension of power, they are generally
36
denoted as P for the dimensional analysis purpose. For the flow and heat transfer
problem under consideration, there are four basic dimensions, [length], [time], [mass]
and [temperature]. Here "[ ]" means "the dimension of". The four basic quantities,
which include all these basic dimensions and at the same time never make a
nondimensional quantity themselves, are chosen as L,, t,, p, and T,. They were used
as the references in the previous nondimensionalization of the governing equations and
boundary conditions. The quantity P is nondimensionalized by dividing it by Lrt;,3pT,
i.e., paV3. From the I theorem (Fox and McDonald 1978), we have
P/I (pa5) = {Pr, y, M, a, (o, T,/T,, Th/T,, b/a, c/a, d/a, E/a, L/a}. (2.35)
CHAPTER 3
NUMERICAL TECHNIQUES AND CODE VALIDATION
The integral equations discussed in the last chapter need to be solved numerically
for the specified boundary conditions. The gas flow under consideration is still
incompressible even though the density of the gas changes considerably with space and
time. This is because the Mach number here is very small, and therefore the density
variation is not caused by pressure gradients, but by temperature gradients and
volumetric changes. So the flow should be called incompressible variable-density flow,
and the SIMPLE algorithm (Patankar 1980), which was originally developed for
incompressible flows and heat transfer, can be adopted as the proper solution algorithm.
To solve the problem with moving boundaries, the moving grid technique by Demird&iC
and PeriC (1990) is incorporated into the SIMPLE algorithm. The pressure-splitting
method proposed in the last chapter and the code developed accordingly are validated
with some cases whose solutions are known beforehand so that the code can be used for
production runs.
3.1 Numerical Grid Generation
In order for the integral equations to be solved numerically, the entire domain is
discretized with grids into a system of subdomains, or cells. At every time step, new
grids are generated according to the oscillating-piston locations. Owing to the
rectangular nature of the boundaries of the models under consideration, Cartesian
rectangular grids are used. Figures 3.1 (a) and 3.1 (b) show two typical grid
distributions at time instants wt = ir/4 and i, respectively. The domain is decomposed
into three blocks with each in the left chamber, the channel and the right chamber. The
total number of grid points in each direction within a particular block is kept constant and
the grids move and deform only in the x-direction but not in the y-direction.
Consider the grid distribution in the x-direction. Suppose there are N grid points
to be distributed along a line segment parallel to the x-axis, starting at x=x, and ending
at x=XN with the first and last grid spacings being specified as Axo and Ax,,
respectively (Fig. 3.2). Grid points are distributed according to the hyperbolic tangent
distribution given by Thompson et al. (1985) as follows.
First define
G -x (3.1)
XN XI
N-l
and
H Ax= (3.2)
S Axo
Then solve one of the following nonlinear equations for the parameter 6 by using the
Newton-Raphson method (Hornbeck 1982):
Ij '_ i i l-. .
"'" I i T i :: '"i .i i..
777-1--l I
= : : : .:: :; : :
I .
: I : Iilt -ii: : : : :,,::
Fig. 3.1 Two typical grids. (a) wt=r/4; (b) wt=wr.
Ax, AX
X1 X2 X ,,. Xi
S XN- XN- XN X
Fig. 3.2 ID grids
sinha 1 for G < 1, (3.3)
a G
sin = 1 for G > 1 (3.4)
a G
The i-th grid is distributed at
xi x + i (3.5)
H+(1-H)gi
where
i-i 1
tanh[ ( 1)]
gi 1 + N-1 (3.6)
2 6
tanh-8
2
For each chamber block, either XN or x, varies with time due to piston
oscillations. The two end spacings are given in terms of the instantaneous average grid
spacing at each time instant so that a good moving grid distribution is achieved. The
hyperbolic tangent distribution is also used for generating the grid in the channel block
with each end spacing being the same as the respective adjacent end spacing of the
chamber block, resulting in a smooth grid distribution near the chamber and channel
interfaces. Similar procedures are applied to generate grids along the y-direction.
3.2 Discretization of the Governing Equations
As mentioned earlier, the governing equations need to be discretized on a finite
number of cells (i.e., finite volumes), so that they can be solved numerically. The
staggered grid layout is adopted here to avoid the pressure-velocity decoupling (see
Patankar 1980). Figure 3.3 shows that the control volumes for scalar quantities p, p and
T are the same and the quantities are located at the centers of the control volumes. A
typical scalar control volume is highlighted in the figure. Velocity components (u and
v) are located at the interfaces of the scalar control volumes. As a result, the u control
volumes shift in the x-direction by half cell distance from the scalar control volumes, and
similarly v control volumes also shift in the y-direction by half cell distance. Figure 3.4
shows a typical u control volume and a typical v control volume. Each of the quantities
indicated there is regarded as the average value over the relevant control volume. It can
be seen from the above discussion that for the 2D staggered grid layout, there are three
sets of control volumes: one set for all scalar quantities, and the other two for the two
velocity components.
The general equation (2.34) is taken here as an example to demonstrate the
procedures for discretizing the governing equations. Consider a cell centered at point
P with neighboring cells centered at points N, W, S and E (as shown in Fig. 3.5). The
midpoints of the interfaces between cell P and each of its neighboring cells are denoted
by lower case letters n, w, s and e. As mentioned in the last section, the grid points
move in the x-direction. So Fig. 3.5 also shows the grid velocities: (u,),-at the east
interface, and (ug),--at the west interface. The grid velocities are evaluated by using first
order finite difference based on the grid coordinates at current and last time levels.
From Eq. (2.34), the integral equation for the general quantity 0 over cell P is
Fi-. 3.3 Control volumes for scalar quantities
Fig. 3.3 Control volumes for scalar quantities
P t. T
u control volume
Fig. 3.4 Staggered u and v control volumes
v control
volume
d f p4dV + f p(9V-Vb,)adS = VfrV A idS + f cdV,
vd(t) sTt) s(t) Vp(t)
where the subscript "P" denotes the quantity of cell P.
.-Axw ---x AXp A---
N
AyN *
n vn
A, V\ w P e 'E
s 1
Ay
II s___
(3.7)
Fig. 3.5 Schematic of a general cell for 0 at point P
Consider the rate-of-change term first. In reference to Fig. 3.5, this term can be
discretized as
df P,V = pppp RN
dt At
where At is the time step size, Vp the cell volume, and the superscript "o" represents a
quantity at the last time level. Here only a 2D case is considered, so
44
V, = Ap = AXpayp, (3.9)
and consequently,
d f pOdV = P _- p;A (3.10)
dt v At
The convection term is discretized as
fp( -Vb)-fdS = pOe(u -u,),YP -p (U -uAyp-),wAy
st) (3.11)
+ pnV Axp- PsVs AXp
Define mass fluxes through the cell faces
F, = [p(u-u,)]. Aye, F. = Lo(u-u,)]w Ayp, 0.12)
F. = (pv) AXp, F, = (pv), Ax, (3.13)
then
fp0)(V- ) -fidS = Fe FP + F n F.< (3.14)
Sp(t)
The diffusion term is then discretized as
frVo-adS = (T )y-(rp )Ay,
t) ax(3.15)
+(r ).Axp-(r),Ax .
ay 0ay
The first derivatives in the above expressions are approximated by first order finite
differencing, e.g.,
( -) = OE (3.16)
ax 8x
Thus one has,
Ay, Ayp
RHS of Eq.(3.15) = rA (E ) l w- )
Ax Axp (3.17)
bay, y,
=Do( E PO -Dw P *W)+Dn(oN-P) -D8( p-$) ,
where the diffusion fluxes are defined as
Ay Ayp Axp AX
D r D = D =r- D. = F- (3.18)
w8x, ayn' s 7
The discretization of the source term is straightforward,
f CdV = VCp = Ap( (3.19)
Substituting all the discretized terms into the integral equation (3.7) gives
pp#p -poAPO~
At +F, +-F +FA- F (3.20)
D.(-E -P)-D,((P--w) +D-(On -4P)-D.(
Notice that the fully implicit first order backward Euler scheme is used. In the above
equation, the 0 values at cell faces, i.e. 0,, 0,, 0, and 4,, have to be expressed in
terms of the values at cell centers, i.e. E,, ~ w w and Os, since the latter are
regarded as unknowns to be solved. Convection schemes are needed to do this (for
46
details, see Patankar 1980, Shyy 1994). Applying a convection scheme to the above
equation and then collecting terms gives the algebraic equation for 4 in the final form:
app = aEbE+awW +a + +as +B = ab+B (3.21)
nb
where, for compactness, the subscript "nb" denotes evaluations at the neighboring points
of point P, and the coefficients are defined by
aE = D.C(IPec,) + max{-F, 0}, (3.22)
aw = D,C(IPec,)) + max{F,, 0}, (3.23)
a = D.C([Pec,,) + max{-F,, 0}, (3.24)
a = DC(IPec,|) + max{F,,, 0}(3.25)
B = CAp + apoo (3.26)
S(pA)t (3.27)
At
ap = aE+aW+aN+as+ap = E aa (3.28)
nb
The symbol "Pec" in these expressions is the cell Peclet number, which is the ratio of
the convection flux to the diffusion flux. The Peclet number at the east face, Pec,, is
defined as
Pece = F, / De. (3.29)
The function C depends on the convection scheme used, and in this work the hybrid
scheme is used. Then C takes the form
C(|PecJ) = max{0,1-('/)|Pec|}. (3.30)
47
This scheme is essentially central difference for I Pec, I< 2, and switches to first-order
upwind for |Pece > 2. For the problem under consideration, the source term is only
linear in 4 at most, and so for simplicity it is just lumped into the term B (Eq. (3.26)).
Special care should be taken in the discretization of boundary cells. As discussed
before, there are only two kinds of boundary conditions involved: Dirichlet and
Neumann. Figure 3.6 shows cell P with the south face being the boundary where the
Dirichlet boundary condition is specified there, namely, 4, is given. In this case, let
point S coincide with point s, then
4s=4,. (3.31)
Here also note the special definition of 6y, at the boundary cell, as could be seen clearly
by comparing it with that in Fig. 3.5. With these two special definitions, the procedures
and formulations to discretize the boundary cell are exactly the same as those for an
interior cell as discussed before. By the same token, all other boundary cells, whether
they are at the north or the east boundary, are discretized.
The Neumann boundary condition is treated in another way. Again, a cell P with
the south face being the boundary where the a4/8y is specified is taken as an example
to demonstrate the boundary treatment. As shown in Fig. 3.7, an imaginary cell S is
added to the south of cell P such that
Ays=Ayp, (3.32)
and so y,= Ayp. (3.33)
Applying the boundary condition with central difference approximation yields
48
|< -- 8x --- 41. 8x" -
Fig. 3.6 A south boundary cell with Dirichlet boundary condition
I--- xw I- 5",----
0
S
-------
Fig. 3.7 A south boundary cell with Neumann boundary condition
(A) _P -s (3.34)
ay 8y,
then,
4s =
Now the boundary cell P can be treated as if it were an interior cell. As mentioned in
the last chapter, all the Neumann boundary conditions are with zero gradients in this
study, consequently the last equation gives 4s = 4, for such special cases.
Throughout the discretization process, linear interpolation is frequently used to
evaluate values whenever they are not readily available, except O's themselves at cell
faces, which are determined by convection schemes instead. For instance, in evaluating
the mass flux F. = pvAxp (refer to Fig. 3.7), if p, is not readily available, as is the
case when the control volume is for any quantity but v, then the density there should be
approximated by linearly interpolating the values at the neighboring points.
It has been seen that the involvement of moving grids doesn't complicate the
algorithm much. It is the velocity relative to the control surface that is used in defining
the convection mass flux (refer to Eqs. (3.12) and (3.13)). Since the grid velocity in this
work is neither zero nor the fluid velocity, the grid is neither Eulerian nor Lagrangian.
3.3 Solution Procedures
With the SIMPLE algorithm, the discretized governing equations are solved one
at a time in a sequential manner. Besides, pressure is used as the basis, i.e., the
50
hydrodynamic pressure is determined by letting the velocity field satisfy the continuity
equation, as is in contrast to the compressible algorithms (also called density-based
algorithms) that use density as the basis by solving the density from the continuity
equation and then determine pressure from the equation of state.
At the very beginning of the calculation, initial conditions must be properly
specified. Then time is increased by At; at the new time level grids are generated based
on the new boundaries and boundary conditions are specified accordingly. The flow and
thermal fields need to be calculated, and the procedures are described as follows.
First, all the unknown quantities are guessed, and denoted as u*, v', (p'), (ph)',
T', p', etc. Good initial guesses should correspond to the values at the last time level.
With these guessed quantities, the u-momentum equation is discretized first. Since the
way to discretize a general equation was already given in detail in the last section, only
the final form of the discretized u-momentum equation for an u-control volume P is given
here.
apup = rab yp(Pl -PPh) +s (3.36)
nb
where the source term is intentionally written into two parts such that the part solely
related to the driving force--pressure--is separated from all other contributions, as this
will help the discussion of pressure-velocity coupling.
Equation (3.36) is linearized already because the coefficients involved are
evaluated in terms of the guessed u* values but not the current u values. Every cell
results in an equation like Eq. (3.36), and thus a system of linear algebraic equations for
u is formed, which is solved by an iterative method--the line successive relaxation
method (Hirsch 1988). The u* values can now be updated by the newly obtained u
values,
u' = WOu + (1-oW)u, (3.37)
where w, is an underrelaxation parameter (0
between the u-momentum equation and other equations, underrelaxation would be needed
in updating u' values so that the solution procedure converges.
The discretized v-momentum equation is given as
ap = b -Axp(pNh -p ) +, (3.38)
nb
In a similar way, v is solved and v' is updated with underrelaxation.
The velocity field determined this way, however, may not necessarily satisfy the
continuity equation because the solution is based on a guessed hydrodynamic pressure
field (ph)*. Therefore, the discretizing the continuity equation over a pressure-control
volume at point P (refer to Fig. 3.5) yields
(p'A)-(pA)? (3.39
(pA)1-(pA), +F: -F +F -F. = sF (= 0) .(3.39)
At
The pressure needs to be corrected so that the resultant velocity field satisfies the
continuity equation. If the pressure correction is denoted by p', then, from the
momentum equations (3.36) and (3.38), the corresponding velocity corrections at the
interfaces of the pressure-control volume are given by
s -Y (PE-PP) u, = -A (p -Pw) (3.40)
ap ap
and
n = --A N-P ), v. (Ax( -P's) (3.41)
ap ap
Since there is a half-cell shift between the velocity- and pressure-control volumes, an
index shift is also applied when deriving the above equations from the momentum
equations. The coefficient a, in Eqs. (3.40) differ from those in (3.41) because the u-
and v-control volumes are not collocated. The corresponding flux corrections are
Ayp ,
Fe = Pu'Ayp = -PeA PEY(P-)Ayp
Ay
F, = pu'Ayp = -p, (pp-pw)Ay ,
ap
F = p.vAxp = -p.- (p-ps)A xp
ap
It is required that the corrected pressure and velocity satisfy
(p'A)p-(pA) +F*-F-+F.-F,+F-FF-F = 0. (3.43)
At
Substituting Eqs. (3.42) into Eq. (3.43) and collecting terms give the pressure correction
equation
b,pp, = b E'S, (3.44)
ab
where b's are the coefficients, and s', is defined in Eq. (3.39). The latter quantity
serves as a measure of how closely the continuity equation is satisfied.
It was mentioned in the last chapter that there are no explicit pressure boundary
conditions for the problems under consideration. Pressure correction equations like Eq.
(3.44), however, do need boundary conditions. Boundary conditions for p' need to be
derived from the given velocity boundary conditions. For example, consider a cell with
the west face as the boundary where the normal velocity component, u, is predefined,
then there should be no velocity (and therefore mass flux) correction there, i.e., F,' in
Eqs. (3.42) is zero. The F,' term will not appear in Eq. (3.43), so the coefficient bw
in Eq. (3.44) takes the value of zero (bw = 0).
After pressure corrections are solved from the pressure correction equation, the
hydrodynamic pressure is updated by
ph (h)* + pp', (3.45)
where underrelaxation is used again to maintain convergence (0< ,< 1).
The corresponding velocity corrections u' and v' are calculated from Eqs. (3.40)
and (3.41), respectively, and the velocity components are then updated as
u=u'+u', (3.46)
v=v'+v'. (3.47)
One should note that the energy equation is coupled with the continuity and
momentum equations, and the reasons are 1) viscosity is a function of temperature, and
2) for the variable-density problem, density is also a function of temperature. So the
discretized energy equation is solved with the guessed density p' and the above updated
velocity components. The thermodynamic pressure p' can then be computed by using the
method described in the last chapter via Eq. (2.23), and density is also updated from Eq.
(2.21). Underrelaxation is also used in updating density.
With all the updated quantities, one goes back to the momentum equation again
and repeats the solution procedures until the solutions are fully converged at the given
time level. The time level is next increased by one when solutions have converged.
Since the time-periodic problem is being simulated, it would take several cycles (or
periods) for the numerical solution to reach a time-periodic state. Figure 3.8 gives a
summary of the overall solution procedures in a flowchart.
3.4 Initial and Slow-Start Boundary Conditions
Because of the iterative feature of the solution procedure, good initial guesses for
the unknown quantities can make the solution converge fast, as is especially the case at
the very first time step. At t=0, the left piston velocity is zero (Eq. (2.24)), but the
right piston velocity is not if op, d 0 (Eq. (2.25)). Meanwhile, if TTh, it would be
difficult to specify the precise initial velocity and temperature conditions for the solution
process to start with because the initial velocity and temperature fields are not uniform
in such conditions. To circumvent such a difficulty, a slow-start procedure was utilized
in specifying temperature boundary conditions, and the phase angle lag p, and therefore
the velocity boundary conditions as well.
Ilnput Pr, Re(=az), M,y, Geometry, I.C. & B.C. at t=0
Generate Grid. Update B.C.
IGuess u, v, p', p", T and p at t
u-Momentum Eq. u -
0)
E
.c Iv-Momentum Eq. = v
-C
(0 II
2 IContinuity Eq. p Update p", u, v
lEnergy Eq. =-; T
lUpdate p'. Eq. of State -pi
Converged at t? No
SYes
No
o Time-Periodic?
Yes
Stop
Fig. 3.8 Solution procedures for variable-density time-periodic
flow and heat transfer problems
A quantity ip is defined by the following expression,
[ 1 -cos(2wt) 0
S' (3.48)
where o, is the desirable phase angle. The newly defined f will take the place of (p in
the expression for Xrih, in Figs. 1.3 and 1.4, and also in Eq. (2.25) for u,. Then at
t=0, both pistons have zero velocity, and so the velocity in the whole computational
domain is also zero. It takes a quarter of the first cycle for p to adjust smoothly to its
final value p,.
For temperature boundary conditions, suppose that the desirable chamber wall
temperatures are: T, for the left chamber and Th for the right, then define
1 (Th-T)
S {(Th+T)- ( [1 -cos(2ct)]} 0
'le
T. ct ,
and
1 (Th -T) r
{ ~((Th+Tc)+ -cos(2wt)]}, 0 wt < (3
h 22 (3.50)
Trght 2
Th ,
where T,, and Tr, are the temperatures to be specified as the boundary conditions for
the corresponding chamber walls. At t=0, Tft=Tr,h, hence a uniform initial
temperature field is expected. Over a quarter of the first cycle, the boundary
temperatures gradually develop to the desirable values.
3.5 Code Validation
The proposed pressure-splitting method and the code developed were validated for
variable-density, time-periodic flows and heat transfer in regions with moving
boundaries. The first test case is a 2D piston-cylinder set-up with all the boundaries
being insulated (Fig. 3.9). Air is present in the cylinder with initial pressure po= 1 std
atm (= 1.013x105 N/m2), temperature T,=288.15 K, and density p,= 1.225 kg/m3. At
the very beginning, the piston stands still at x=0, so there is no fluid motion. Then the
piston oscillates, and its location is given by
I cost() 1, O wt < r, (3.51)
ecos(Wt), 7< 5Ct.
The oscillation frequency is f=l Hz (w=2rf=2r rad/s). The geometric parameters are
a=lxl03 m, e=la, and L=3a. The nondimensional parameters are ,y=1.4, Pr=0.72,
a=0.6559 and M=1.846x105. It is known from thermodynamics that for such small
dimensions and slow oscillations the thermodynamic quantities should be uniform within
the cylinder, and the quantities can be determined analytically by the adiabatic
relationships:
Y -Insulating wall
SL--
a
* a
x
xpisn = cos(ot)
Fig. 3.9 A piston-cylinder set-up
No.1 Single preclson, Without pressure-spliting
No.2 Double precision, Without presure-spllttlng
No.3 Single precision, With pressure-spltting
15 co t=7i/4, No.1 I
It
Symbol at Computation I
( -ix- 7/4 No. 1 I I
0 II I I
S10 + 21J3 No. 1 I I
o ,/4 No. 2 or3 I Ii
A 2n/3 No.2or3
SI i Ill
iI j I I
1 i I
S ot=27r3 k
No. 1 .1
0 ..*...1 --- *-* ****
wt=2n/3
No.2or3 o*'O DO wt=4, No.2or3
-1 0 1 2
x/a
Fig. 3.10 Computed hydrodynamic pressure distributions along
the x axis for the piston-cylinder set-up
3
t V ( T
SP (o)Y = (_)y =( )i (3.52)
o Po T
where v is the total volume and the subscript "o" denotes initial values.
This case was computed numerically on a machine with significant digits of seven
for single precision and of 16 for double precision. It was first computed with single
precision and without split pressures (denoted as Computation No. 1), next it was
computed with double precision and still without split pressures (Computation No. 2),
and finally the same case was computed with single precision and with split pressures
(Computation No. 3). The numerical pressure distributions along the x axis are plotted
against each other in Fig. 3.10. The pressure shown in the figure is actually the
pressure difference with respect to the point (x,y)=(1.5a,0), therefore the pressure is the
hydrodynamic pressure p'. As indicated by the figure, at time wt=ir/4, the
hydrodynamic pressure from Computation No. 1 oscillates wildly along the x axis, but
at ot=2r/3, it is just uniform. None of these two distributions is physically realistic.
When the computation accuracy is increased by using double precision (Computation No.
2), the hydrodynamic pressure distribution is neither oscillating nor uniform. Such an
observation reveals that for problems like this, the pressure field cannot be correctly
resolved by using the computational algorithm of Computation No. 1, and using double
precision in the calculations (Computation No. 2) is indeed one remedy. Figure 3.10
clearly shows that Computation No. 3 yields the same pressure as Computation No. 2
does. Since Computation No. 3 uses single precision, it saves computer CPU time by
about 30%, memory and storage by almost 50% over Computation No. 2. This proves
60
the significant benefits given by using the pressure-splitting algorithm, and verifies the
analysis presented in the last chapter.
The instantaneous thermodynamic pressure, p'(t), in the third oscillation cycle is
shown in Fig. 3.11. To help compare the magnitudes of ph and p', pt in this figure is
normalized with the same quantity as that for ph. By comparing Figs. 3.10 and 3.11, it
can be seen that the magnitude of the thermodynamic pressure is greater than that of the
hydrodynamic pressure by nine orders. It is for this reason that Computation No.1 with
only seven significant digits can not accurately resolve the pressure field when the
thermodynamic and hydrodynamic pressures are not treated separately in the solution.
In Fig. 3.12, the numerical temperature distribution along the x axis is plotted
against the analytical solution. At a given time, the analytical temperature distribution
is uniform, as is presented by the horizontal line in the figure. The temperature from
Computation No. 1 is not uniform, but varies along the x axis. The incorrect
temperature is produced by the erroneous velocity field resulting from the pressure field
that is not accurately resolved numerically. It can also be seen that the numerical
temperature distribution either from Computation No. 2 or from Computation No. 3 is
in perfect agreement with the analytical solution, which again shows that the proposed
pressure-splitting technique not only gives the right solutions, but also saves computer
CPU time, memory and storage substantially.
In the second test case, the set-up is similar to Model Al except that the piston
surfaces are isothermal. Geometric parameters are (refer to Fig. 1.3a): a=lxl0'3 m,
b=4a, c=3a, e=2a and L=6a. The wall temperature (including the piston's) is kept at
3.5x10'
3.0x10"
2.5x10*
S2.0x10"
1.5x10,
0 100 200 300
cot (degrees)
Fig. 3.11 Thermodynamic pressure in the third cycle for the piston-cylinder set-up
0.98 .. ...
Mt=l/4
Symbol ot Computation
x s/4 No.1
0.96 2ni3 No.1
o 'l4 No.2or3
A 2J/3 No.2or3
-Analytical
0.94
No.1 Single preclaon, Without pressure-aptting
No.2 Double precision, Without pressure-splitting
No.3 Single precision, Without pressure-splitting
0.92
0.90 "yloal; No.2, or 3 t
It .. .** *No.1
-1 0 1 2
Fig. 3.12 Temperature distributions along the x axis for the piston-cylinder set-up
62
constant T,=288.15 K (therefore T,=Th=T, ). Pistons oscillate sinusoidally with the
right piston leading the left by 900 in phase (tpo=90) with an oscillation frequency of
f=2.32 Hz. Air is the working medium. At t=, po= 1 atm, T,= 288.15 K, p= 1.225
kg/m3. The corresponding nondimensional parameters are y=1.4, Pr=0.72, a=1 and
M=4.29x10-5.
In this calculation, the pressure-splitting algorithm is used. For the computational
domain, which is only the upper half geometry, a 91x31 grid system is used. The grid
distribution is: 31x31 in the left chamber, 31x11 in the channel, and 31x31 in the right
chamber. The time step is chosen such that each cycle is discretized into 360 uniform
intervals, i.e., At= 1/(360f)= 1.2x 10 sec. It takes about four cycles for the solution to
reach the time-periodic state.
The P-V diagram for this calculation is shown in Fig. 3.13. This figure shows
that the minimum pressure is p,/po=0.7099 and the maximum is p,, /Po= 1.6420,
which can be checked by examining two extreme conditions, namely isothermal and
adiabatic. For an isothermal condition,
p/p. = v,/v,
so, Pl,/Po = V/V, = 0.7262, and p,/p, = VJVo = 1.6055.
For an adiabatic condition,
p/p. = (v,/v),
so, m,J/P = (VJv.,) = 0.6389, and p,,/po = (Vo/Vd = 1.9402.
The numerical maximum should be larger than the adiabatic counterpart and less than the
isothermal one, whereas the numerical minimum should be larger than the isothermal one
0.6 0.8 1.0 1.2 1.4
V / V,
Fig. 3.13 P-V diagram of the second test case
and less than the adiabatic one; moreover, numerical values should be closer to the
isothermal values than to the adiabatic ones because of the boundary condition of
isothermal walls. That is indeed the case:
minimum: 0.6389 (adiabatic) < 0.7099 (numerical) < 0.7262 isothermall),
maximum: 1.6035 isothermall) < 1.6420 (numerical) < 1.9402 (adiabatic).
The correctness of the numerical solution is also supported by checking the
overall energy balance. First define the power input--the mean rate at which work is
done on the gas by some external force-as
p = f -pdV)dt, (3.53)
0 V(t)
where r is the period of piston oscillations, i.e. r = 2r/lu.
The mean rate at which heat is transferred out of some domain is
Q f 1 f -kVT-fdS}dt (3.54)
Q 0 osct)
where S,,(t) represents the solid surface that is in contact with the gas, and the surface
may change with time.
The rate at which the gas "absorbs" heat through the left chamber walls plus the
piston surface is Q,,a=2.4589D watts (here D is the depth in the third dimension in
meters); the rate at which heat is rejected out of the gas through the channel walls is
Qd =0.5516D watts, and that through the right chamber walls plus the piston surface
is Q,,h=2.7297D watts. The power input is P=0.8415D watts. According to the first
law of thermodynamics, the following expression holds,
Qeat Qcha-l Qnth = -P.
Then, LHS = 2.4589 0.5516 2.7297 = -0.8224,
RHS = -0.8415.
Therefore, there is only a very small relative error of 2.3% in the total energy balance.
The preceding two test cases have verified that the pressure-splitting algorithm is
computationally accurate and efficient, and in the mean time they have also validated the
computer code.
CHAPTER 4
RESULTS AND DISCUSSION
By using this code, investigations were made on the models' cooling capacity and
efficiency in terms of working gases, oscillation frequencies, channel length, ratios of
the chamber width to the channel width, etc. The velocity and temperature fields were
also examined. Discussion is made on the criterion for the model to result in good
regeneration. Table 4.1 lists the cases that were studied. Note that to make the number
of parameters appearing in the table as small as possible, yet adequate to present each
case, those parameters that can be obtained directly, or via auxiliary relations if needed,
are omitted there. For instance, parameters that do not need such relations are: the ratio
of specific heats (7= 1.4 for air and 7= 1.7 for helium), Prandtl number (Pr=0.72 for
air and Pr=0.67 for helium), and the constant T, used in the Southerland's formula (Eq.
(3.9)) (T=110.4 K for air and T,=61.67 K for helium). One of the parameters that
need such relations is the reference density p,, which can be obtained via the equation
of state as p,=p, /(RT,) with R being 287 J/(Kg-K) for air and 2077.2 J/(Kg-K) for
helium. For the meanings of the geometric parameters listed in Table 4.1, please refer
to Figs. 1.3 and 1.4.
Tabel 4.1 List of the cases investigated
Case Model a b/a c/a d/a c/a L/a ,p. Gas p, T, T,/T, T,/T, f a M
No. [m] [] [atm] [K] [Hz]
x103 xl0'
1 B1 1 8 8 N/A' 7 12 90 Helium 1 288.2 0.9 1.1 1 0.229 6.291
2 Al 1 8 8 N/A 7 24 90 Helium 1 288.2 0.9 1.1 1 0.229 6.291
3 A2 1 8 8 6 7 24 90 Helium 1 288.2 0.9 1.1 1 0.229 6.291
4 B2 1 8 8 6 7 24 90 Helium 1 288.2 0.9 1.1 1 0.229 6.291
5 B2 1 8 8 6 7 24 90 Air 1 288.2 0.9 1.1 1 0.656 18.46
6 Al 1 6 3 N/A 2 10 45 Helium 1 288.2 0.9 1.1 1 0.229 6.291
7 Al 1 6 3 N/A 2 10 75 Helium 1 288.2 0.9 1.1 1 0.229 6.291
8 Al 1 6 3 N/A 2 10 90 Helium 1 288.2 0.9 1.1 1 0.229 6.291
9 Al 1 6 3 N/A 2 10 105 Helium 1 288.2 0.9 1.1 1 0.229 6.291
10 Al 1 6 3 N/A 2 10 135 Helium 1 288.2 0.9 1.1 1 0.229 6.291
11 Al 1 6 3 N/A 2 10 90 Helium 1 288.2 0.95 1.05 1 0.229 6.291
12 Al 1 6 3 N/A 2 10 90 Helium 1 288.2 0.85 1.15 1 0.229 6.291
13 Al 1 6 3 N/A 2 10 90 Helium 1 288.2 0.8 1.2 1 0.229 6.291
14 Al 1 6 3 N/A 2 10 90 Helium 1 288.2 0.7716 1.2284 1 0.229 6.291
15 Al 1 6 3 N/A 2 10 90 helium 1 288.2 0.75 1.25 1 0.229 6.291
16 Al 1 8 8 N/A 7 12 90 Helium 1 288.2 0.9 1.1 1 0.229 6.291
17 Bl 1 4 8 N/A 7 12 90 Helium 1 288.2 0.9 1.1 1 0.229 6.291
* N/A denotes not applicable for the case.
Table 4.1--Continued
Case Model a b/a c/a d/a e/a L/a 'Po Gas p, T, T /T, Th ,T, f a M
No. [m] [ ] [atm] [K] [Hzl
xl03 x10l
18 B1 1 12 8 N/A 7 12 90 Helium 1 288.2 0.9 1.1 1 0.229 6.291
19 B1 1 15 8 N/A 7 12 90 Helium 1 288.2 0.9 1.1 1 0.229 6.291
20 Al 1 6 3 N/A 2 10 90 Helium 1 288.2 0.9 1.1 10 0.724 62.91
4.1 Velocity and Thermodynamic Fields of a Tyvical Case
In this section, the computations and results of Case 1 (refer to Table 4.1) are
examined. Computations were made on the Cray C90 supercomputer of the Pittsburgh
Supercomputer Center. Different grids and time step sizes were first tested to make sure
that they were fine enough for the numerical solutions to be accurate enough. It was
identified that for this particular case a reasonable choice was: a grid system of 2153
points (31x31 in the left chamber, 21x11 in the channel, and 31x31 in the right channel)
and a time step size of At = r/360 (= 2.78xl03 sec.), for further refining either the grid
or the time step size did not result in any noticeable changes in the numerical solutions.
The grids at wt= T/4 and ut=r for this case can also be found in Figs. 3.1 (a) and (b)
respectively, though these figures originally served another purpose there. Solution
iteration numbers used were: 8 iterations for solving the u-momentum equation, 8 for the
v-momentum equation, 40 for the p' equation and 8 for the energy equation; and 300
repetitions of the above procedures at a give time level. Under these given conditions,
it took about 5472 second CPU time to run one cycle on the Cray C90 machine.
Figures 4.1 (a) and (b) plot the transient processes of the thermodynamic pressure
and the pressure drop between the left and right pistons. These figures reflect the slow-
start velocity and temperature conditions in the initial quarter cycle discussed at the end
of the last chapter. Figure 4.1 (b) shows that the hydrodynamicc) pressure drop adjusts
itself to time-periodic state only in few cycles. The thermodynamic pressure oscillates
in the range of 0.5 to 2.5 atms, whereas the hydrodynamic pressure drop is only in the
range of -2x10-6 to 2x0I atms, with the latter being only about one millionth of the
(a) Thermodynamic pressure
a=lmm, b=8a, c=8a, e=7a, L=12a, f=1Hz, ,.=90
Helium, p,=latm, T,=288.15K, T=0.9T,, Th11T, a=0.229
t / cycle
(b) Pressure drop
4x10"
2x10"6
0
-2x 106
a=lmm, b=8a, c=8a, =7a, L=12a, f=1Hz, p.=90*
Helium, p,=latm, T,=288.15K, T.=0.9T,, T,=1.1T,, a=0.229
t / T cycle
Fig. 4.1 Transient processes of pressures.
(a) Thermodynamic pressure; (b) Pressure drop between pistons.
former.
A way to check if numerical solutions have reached the time-periodic state is to
compare the current solutions with those at the time one cycle earlier. The quantity a
defined below is such a measure.
o = i1 (4.1)
where 4 is a general notation representing any specific dependent variable of u, v, ph,
T or p. Such a measure is objective, because it takes the overall domain into account
by the summation and meanwhile is a relative measure by adding the denominator in the
definition. Figure 4.2 plots the quantities a,, ao, ph and or in the logarithmic scale
vs. time. At any given time, they are of the same order and they all decrease linearly
as time increases. The lines of logo's have exactly the same slope.
Figures 4.3 (a), (b), (c) and (d) give the velocity vectors at time instants wt=0,
7/2, 7r and 3r/2, respectively. Note that to have a better view of the field quantity in
the configuration, the domain shown in the figures is actually as large as six times the
computational domain. No significant circulation region is observed for this case because
of the small Womersley number used.
The corresponding hydrodynamic pressure contours are plotted in Figs. 4.4 (a),
(b), (c) and (d). These figures show that when neither piston is extremely close to the
channel exit, pressure drops mainly occur in the channels, otherwise the pressure drop
in the transverse direction within the chamber is comparable to that in the channel. The
a=1 mm, b=8a, c=8a, e=7a, L=12a, f=1Hz, p,=90
Helium, p,=1 atm, T,=288.15K, T,=0.9T,, T,=1.1T,, a=0.229
2.0 2.5 3.0
t / T cycle
3.5 4.0
Fig. 4.2 Convergence history to time-periodic state
100(a))=0.628m/s
(a) ot=O (b) ot=1/2
20 20
10 10 -
-10 -10
-20 -20
-10 0 10 20 -20 -10 0 10
x/a x/a
Fig. 4.3 Velocity fields of Case 1. (a) ut=O; (b) wt=r/2; (c) wt=T; (d) wt=3r/2.
100(at)=0.628m/s
(c) At=7I
:. .
------ ,,:-- :
,zc2 I ..
-10 0
x/a
10 20 -20
-10 0
x/a
Fig. 4.3-continued
10 20
' i ,
_li:
_
~-~.. .._........._._.1 ~4~. ~
..~.~..~.~.~.~..--' L4~C~L .
(d) cot=37/2
(a) ot=O
20
10
C 0
-10
-20
Aph=500p,(ao)2=3.34x1 0"N/m2=3.30x1 Oatm
Fig. 4.4 Hydrodynamic pressure contours of Case 1. (a) wt=0; (b) jt=r/2; (c) ot=r; (d) wt=3-/2.
I TfI r Ti7IT,'
T i1 iiI
i
- -IU 1
-20 '
20 -20 -10 0 10
(b) ot=d/2
i I i
,"lrr mnrr~ij~
(c) (Ot=l7
10 20 -20
Aph=500p,(ao)2=3.34x10- N/m2=3.30x1 O-atm
Fig. 4.4--continued
-10 0
x/a
-10 0
x/a
10 20
(d) oAt=3n/2
S-------------
: .^ ^ .li _)
1 i-
total pressure drop in space is only a few millionths of an atmosphere at most. It will
be seen later that this is in contrast with the thermodynamic pressure, which changes with
time by an amount in the order of one atmosphere.
Figures 4.5 (a), (b), (c) and (d) demonstrate the temperature contours. At all the
times, the temperature of the gas within the channel is very close to the local channel
wall temperature, indicating that the channel functions as a regenerator properly. At
wt=0 (Fig. 4.5 (a)), which is at the compression stage, the temperature in the right
chamber is much higher than that of the right chamber wall (1.1 nondimensional), and
a value as high as 1.625 is observed. Inversely, at wt=a- (Fig. 4.5 (c)), which is at the
expansion stage, the temperature in the left chamber is much lower than that of the left
chamber wall (0.9), and a value of 0.625 is observed there. On average over time and
space the mean gas temperature in the right chamber is higher than that of the chamber
wall, so heat is rejected out of the gas through the right chamber wall; on the other
hand, the mean gas temperature in the left chamber is lower than that of the chamber
wall, and heat is absorbed by the gas there.
Density contours, shown in Fig. 4.6 (a), (b), (c) and (d), indicate that density
changes significantly with time and space. For the time instants chosen, the
nondimensional density takes values from 0.55 to 2.4. At any time, the patterns of
density contours look like those of temperature contours due to the reciprocal relationship
between density and temperature at a given time (Eq. (3.9)).
Figure 4.7 plots the P-V diagram for this case. The direction of the loop is
counterclockwise, indicating that some external force does work on the gas. The power
(a) ot=O (b) ot= n/2 T/T
| j|A| 1.625
20 20
/'W 1 475
10 10
uiT "'
-10 -- ,:jR I T TlT I r '-
to a P 1.125
0F 0.65
-10 0 10 20 -20 -10 0 10
x/a-20 -20 o
x/a x/a **
Fig. 4.5 Temperature contours of Case 1. (a) wt=0; (b) wt=r/2; (c) wt=w; (d) wt=3-/2.
(c) (t=in
-20 -10 0
x/a
10 20 -20 -10 0
x/a
20
10
S0
-10
-20
Fig. 4.5--continued
10 20
I l l l lr/t l l l L L
i/I/Ill l/ l)I
(d) wt=37/2
(a) At=O
-10 0 1(
x/a
20
10
-a
-10
-10
20 -20 -10 0 1(
x/a
Fig. 4.6 Density contours of Case 1. (a) wt=0; (b) tt=r/2; (c) wt=ir; (d) tt=3r/2.
(b) ot=n/2
(d) ct=37/2
20 -20 -10
10 20
Fig. 4.6--continued
I
ii
: I
20
10
-10
-20
-20 -10
' '
I I
(c) wt=n
) )
k
2.5
2.0
0
- 1.5
1.0
0.5
Fig. 4.7 P-V diagram of Case 1
0.5 1.0 1.5
V I V
82
input is P = 14.075D watts (as mentioned before, the symbol "D" represents the third
dimension in meters). The rate at which the gas absorbs heat through the solid surfaces
having cold temperature (T,) is defined as the cooling capacity and is denoted by Q,
while the rate at which heat is rejected from the gas through the solid surfaces having hot
temperature (T,) is defined as the heating rate and denoted by Qh. Heat can also be
conducted into or out of the channel walls along which the temperature distribution is
linear. This is an undesirable leakage, the rate of which is denoted by Q,. The
quantities here are: Q, = 11.574D watts, Qh = 21.643D watts, and Q,1 = 3.990D
watts. So the total rate of energy into the system is
P + Qc = 14.075D + 11.574D = 25.649D,
and the total rate of energy out of the system is
Q, + Q,, = 21.643D + 3.990D = 25.633D.
Hence the overall energy balance is excellent.
In the above discussions, the factor "D" and unit "watt" always go with the
quantities: Q,, Q ,Q, and P. The factor and units will be omitted in the following for
simplicity.
4.2 Model Comparison
To compare Models Al and A2, consider Cases 2 and 3. As listed in Table 4.1,
all the conditions for these cases are the same except that in Case 2 the temperature along
the entire channel wall is linear, while in Case 3 the channel wall temperature is linear
only along the middle segment of the channel and the wall temperature is constant along
either end segment. Table 4.2 Table 4.2 Comparison of Models
Al and A2
compares the performances of the two
Case Case3
cases. The negative value of Q1, in Model Al ase 3
Model Al Model A2
Case 3 means that heat is added into Q. 12.7403 11.8207
the gas through the channel walls Qh 17.4084 19.3578
Q1, 2.4482 -0.3929
having a linear temperature P 7.1340 7.1626
distribution. In the table, COP stands COP 1.7859 1.6503
for the coefficient of performance, 7r 39.69% 36.67%
which is defined as (Walker 1973)
COP = Cooling capacity / Power input = Q/ P. (4.2)
The COP of Case 2 is 1.7859, and that of Case 3 is 1.6503; the latter is only slightly
lower than the former. According to the Carnot theorem, of all the thermodynamic
cycles with common values for the maximum and minimum temperatures, pressures and
volumes, the Carnot cycle has the maximum COP. For a cooling device, the COP of
a Carnot cycle is (Walker 1973)
(COP). = T, / (T T.). (4.3)
In Table 4.2, the symbol i1, represents the relative efficiency defined by (Walker 1973)
i, = Actual COP / Carnot COP = COP / (COP) (4.4)
The relative efficiency of an actual device usually varies between 40% and a few percent
depending on the size of the device and the operation conditions. Here the relative
efficiency of either Case 3 or 4 is only slightly below 40%. In terms of the cooling
capacity, power input and hence COP, Models Al and A2 have almost the same
performance under the given conditions.
Next compare Models A2 and B2. Case 3 uses the single-channel model A2 and
Case 4 uses the multichannel model B2, with all other conditions being the same. It can
be seen from Table 4.3 that Case 4 gives a cooling capacity of 11.3202, which is quite
close to Case 3's 11.8207. Case 4
Table 4.3 Comparison of Models
requires a power input of 11.3809, while A2 and B2
Case 3 requires only 7.1626. As a result,Case 3 Case 4
Case 3 Case 4
Case 4 has a COP as low as 0.9947. The Model A2 Model B2
Q, 11.8207 11.3202
use of the multichannel model in Case 4 19.3 2.
Qh 19.3578 22.9956
makes the compression and expansion Q1,, -0.3929 -0.3421
P 7.1626 11.3809
processes closer to adiabatic than in Case P .1626 11.3809
COP 1.6503 0.9947
3 where single-channel model is used, so
n 36.67% 22.10%
Case 4 should have a larger range of
pressure variations despite the same
volumetric changes for the two cases. Computed data (not tabulated here) indicate
pressure limits of [0.64, 2.12] (atm) for Case 3 and limits of [0.57, 2.22] for Case 4.
That explains why the latter case needs a larger power input than the former. The
compression and expansion processes in Case 4 produce lower temperature in the left
chamber and higher temperature in the right chamber, respectively, than in Case 3; the
cooling or heating of Case 4, however, is still close to the counterpart of Case 3, because
the heat exchanging area of Case 4 is much smaller than that of Case 3. Even though
85
the multichannel model has a lower COP than the single-channel model does, the former
is more practical from an engineering viewpoint.
4.3 Comparison of Using Helium and Air
The most widely used working gases for Stirling-type machines are air, helium
and hydrogen. The selection of a gas for
Table 4.4 Comparison of helium and air
a specific Stirling cooling device not only
depends on the thermodynamic Case4 Case 5
Helium Air
performance of the gas, but also depends Q. 11.3202 5.2169
on other factors such as cost, availability, Qh 22.9956 14.1773
QI. -0.3421 1.2477
safety and so on (Reader and Hooper P 11.3809 10.2219
1983). With the computer code, it is COP 0.9947 0.5104
very convenient to compare the ?Ir 22.10% 11.34%
performances of working gases. Consider
Cases 4 and 5, in which geometric parameters, temperature boundary conditions, piston
oscillation frequencies and amplitudes are the same, but helium is the working gas in the
former and air in the latter. From Table 4.4 it can be seen that helium produces a
cooling capacity twice as large as air does, while the power input needed for one case
is close to the other. As a result, helium is about twice as efficient as air under the
given conditions.
There are mainly two reasons contributing to this. First, the monatomic gas
helium has a ratio of specific heats of 1.67, which is higher than air's 1.40. So, for the
86
same amount of expansion, more cooling is produced by helium than by air. Second,
helium has a smaller Prandtl number than air (0.67 versus 0.72) and a kinematic
viscosity eight times as large as air. As result, the channel's heat regeneration for
helium is better than for air. And a detailed discussion on the regeneration will be found
in Section 4.8.
4.4 Phase Angle Effect
In order to examine the effect of phase angle on the performances of the models,
five runs (Cases 6 through 10 in Table 4.1) were made with phase angles of o, = 45,
75, 90, 1050 and 1350, respectively, with other parameters being maintained constant.
Figure 4.8 plots the cooling capacity (Q,), rate of heat rejected (Q,), power input (P)
and COP versus the phase angle o. The figure shows that each of the four quantities
takes its maximum in the range of phase angles covered, and that the phase angle
corresponding to the maximum value of one quantity may differ from that of another.
For instance, the power input has a maximum at (p, = 1100, the cooling capacity a
maximum at op = 960, and COP a maximum at op, = 75. The optimum phase angle
for the cooling power is not optimum for COP. Generally speaking, the phase angles
corresponding to the maxima are not too far off 90, and each of the four quantities
shown in the figure is not very sensitive to changes in po near the optimum phase angle.
These observations based on the numerical simulations are in agreement with
experimental data (Reader and Hooper 1983, pp. 75-77).
a=1 mm, b=6a, c=3a, e=2a, L=1 Oa, f=l Hz
Helium, p,=1atm, T,=288.15K, T,=0.9T,, Th=1.1T,, a=0.229
COP
.1.0 Cooling capacity, Q 1.0 0
.Power input, P 0.5
0.0 0.0
45 60 75 90 105 120 135
Phase angle 9, ()
Fig. 4.8 Phase angle effects
4.5 Impact of Varying Refrigeration Temperatures on Performances
To understand the effects of changing the refrigeration temperature on the
performances, several runs (Cases 2 and 11 through 15) were made, in which the
temperature difference between the cold and hot chambers was changed to a large extent
while all other parameters were kept unchanged. In Fig. 4.9, the cooling capacity,
power input and COP are plotted versus the ratio of the cold to hot-end temperature,
T,/T Also plotted in the figure is the mean rate at which heat is conducted and
convected into the cold chamber through the gas at the cross section of the channel exit.
This rate is regarded as a loss (denoted by Qo,,,), because it will consume some cooling
a=1 mm, b=6a, c=3a, e=2a, L=1 Oa, f=1 Hz, (p,=900
Helium, p,=latm, a=0.229
To I Th
Fig. 4.9 Refrigeration temperature effects
made by expansion.
The cooling capacity, power input, the loss and the COP all vary monotonously
with the temperature ratio. When the refrigeration temperature rises, the cooling
capacity increases, but the power input and the loss decrease, therefore the COP
increases, too. The closer the cold chamber temperature is to the hot-end temperature
(usually the ambient), the better is the performance. For instance, when T, /T, equals
0.76, the COP is about 1. However, as T, /Th is increased to 0.9, the COP is close to
4. Such phenomena were also reported by McDonald and coworkers (McDonald et al.
1994) based on their design analysis and experimental data with an actual Stirling
refrigerator.
When the cold chamber wall temperature is set too low, cooling is no longer
possible. Instead, heat is rejected out of the gas in the left chamber through the walls.
That is indicted by the negative cooling capacity in the figure. If the hot-end temperature
Th is 300 K, the lowest obtainable refrigeration temperature will be about 189 K under
the operation conditions considered. The large heat loss at low refrigeration temperatures
can consume all the cooling produced by expansions.
4.6 Changing Channel Length
Increasing the channel length (or the regenerator length) will decrease the
temperature gradient along the channel, and therefore can reduce the conduction loss.
Comparing Case 2 and Case 16 (Table
4.5) will demonstrate that. The channel Table 4.5 Comparison of cases with
different channel lengths
length in the former is 24a and that in
the latter is 12a, and other parameters in Case 2 Case 16
L=24a L=12a
both cases are the same. The loss with 0.6131 1.2443
the longer channel is 0.6131, which is Qc 12.7403 13.7450
less than 1.2443--the counterpart with the P 7.1340 8.6155
COP 1.7859 1.6000
shorter channel. However, the cooling
capacity with the longer channel is
12.7403, which is surprisingly even less than the value (13.745) of the cooling capacity
with the shorter channel. The reason for this is that although the loss is reduced by
lengthening the channel, meanwhile the volume that cannot be swept by the pistons--the
90
dead volume--is also increased, therefore the cooling produced by expansion decreases.
The net effect (i.e., the available cooling Q) depends on these two competing factors
under specific operation conditions. Despite such complexity, it is always true that
increasing the channel length also increases the COP because the loss is reduced. Table
4.5 does show that the COP with the longer channel is 1.7859 but that with the shorter
one is only 1.6.
4.7 Effects of Changing the Ratio of Chamber Width to Channel Width
Under given conditions, if the ratio of chamber width (2b) to channel width (2a),
i.e. b/a, is too small, the gas in the cold chamber is not effectively separated from that
in the hot; besides, the effective cooling surface is also too small. On the other hand,
if the ratio is too large, the gas in one chamber cannot be promptly transferred to the
other when needed. Either situation may degrade the cooling performance. For the
multichannel models, in addition to the COP, another appropriate measure of the
performance is the cooling capacity per unit area of the chamber cross section, i.e. Qc
/(2bD), which indicates how economically the device uses space. Four cases were run
with the following values of b/a: 4 (Case 18), 8 (Case 1), 12 (Case 18) and 15 (Case
19). Figure 4.10 plots COP and Q /(2bD) versus b/a. The figure shows that both COP
and Q. /(2bD) vary with b/a, and there exists an optimum value of b/a for either quantity.
When b/a is approximately 6.4, the efficiency of performance takes its maximum value
of about 0.84. When b/a is about 11.5, the quantity Q,/(2bD) reaches a maximum value
of about 750 Watts/m2.
