Citation
Variable-density oscillatory flow and heat transfer in 2D models of a Stirling microrefrigerator

Material Information

Title:
Variable-density oscillatory flow and heat transfer in 2D models of a Stirling microrefrigerator
Creator:
Guo, Guobao, 1963-
Publication Date:
Language:
English
Physical Description:
x, 104 leaves : ill. ; 29 cm.

Subjects

Subjects / Keywords:
Boundary conditions ( jstor )
Cooling ( jstor )
Energy conversion ( jstor )
Engineering ( jstor )
Helium ( jstor )
Pistons ( jstor )
Pressure ( jstor )
Refrigerators ( jstor )
Thermodynamics ( jstor )
Velocity ( jstor )
Aerospace Engineering, Mechanics and Engineering Science thesis, Ph. D
Dissertations, Academic -- Aerospace Engineering, Mechanics and Engineering Science -- UF
Fluid mechanics -- Mathematical models ( lcsh )
Heat -- Transmission -- Mathematical models ( lcsh )
Genre:
bibliography ( marcgt )
non-fiction ( marcgt )

Notes

Thesis:
Thesis (Ph. D.)--University of Florida, 1996.
Bibliography:
Includes bibliographical references (leaves 98-103).
General Note:
Typescript.
General Note:
Vita.
Statement of Responsibility:
by Guobao Guo.

Record Information

Source Institution:
University of Florida
Holding Location:
University of Florida
Rights Management:
Copyright [name of dissertation author]. Permission granted to the University of Florida to digitize, archive and distribute this item for non-profit research and educational purposes. Any reuse of this item in excess of fair use or other copyright exemptions requires permission of the copyright holder.
Resource Identifier:
023165152 ( ALEPH )
35010235 ( OCLC )

Downloads

This item has the following downloads:


Full Text









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.




Full Text

PAGE 1

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 1996

PAGE 2

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 i11 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 CFO 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 tharll
PAGE 3

I am truly grateful to my wife Ying for her love patience and encourageme11t 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. 111

PAGE 4

TABLE OF CONTENTS page ACKNOWLEDGMENTS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 LIST OF TABLES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . VI LIST OF FIGURES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Vll ABSTRACT . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . lX 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 VALIDA TI ON . . . . . . . 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 IV

PAGE 5

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 Phas e Angle Effect . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86 4 5 Impact of Varying R ef rigeration Temperatures on Performance s . 87 4 6 Changing Channel Length . . . . . . . . . . . . . . . . . . . . . . . . . 89 4 7 Effects of Changing the Ratio of Chamber Width to Channel Width 90 4 .8 Dis cussio n on Regeneration .. ... .... ..... ....... ... .. ... ............... 91 5 CONCLUDING REMARKS...... ........... ... .... ....... . . .... ........ ..... 95 REF ERE NC E S . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98 BIOGRAPHICAL SKETCH . .... .... ...... ............................. .. ... .. ....... . 104 V

PAGE 6

Table Table 2.1 Table 4.1 Table 4.2 Table 4.3 Table 4.4 Table 4.5 LIST OF TABLES Terms in the governing equation (2. 34) ............................. . List of the cases investigated ........................................... Comparison of Models A 1 and A2 ........................... ..... ... . Comparison of Models A2 and B2 .................................... Comparison of helium and air .......................................... Comparison of cases with different channel length s ...... .......... VI page 35 66 83 84 85 89

PAGE 7

Figur e Fi g. 1.1 F ig 1. 2 F i g 1 .3 F i g 1 4 F i g. 2.1 F i g 3 1 F i g. 3. 2 F i g 3 .3 F i g. 3 4 Fi g. 3.5 F i g. 3 .6 F i g. 3. 7 F i g 3 8 F i g. 3.9 F i g 3.10 F i g. 3 .11 LIST OF FIGURES page St i rling cy c le and a three component Stirling refrigerator c on f i g uration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 Schematic of a five component Stirling refrigerator ............... Singl echannel models. ( a ) Model Al ; ( b ) Model A 2 ... .. ... .. Mu I ti c hannel mod e ls. ( a ) Model B 1 ; (b ) Model B 2 . ... ... . . . Illu s trativ e 1D dis c r e tization . .. .... .. . ......... . ... . ..... ...... Two typical grids. ( a ) wt = 1r / 4 ; ( b ) wt = 1r . ... ... . .. .... .. . .. . 1D g r i ds . .. ............ .. .... .... . . ..... ..... . . . . ... . .. .. . . ... C ontrol volum e s for scalar quantities ..... .. . .... ... ............ .. S ta g ge red u and v c ontrol volum es .. .. .. .. ... ... ... .. .. .. .. .. .. ... S c 11 e matic o f a general cell for at point P ........ .. .. . . . .... . A s outh boundary ce ll with Diri c hl e t boundary condition .. . . . A s ou t h boundar y ce ll with N e umann boundary condition .. ..... S ol uti o n pro ce dur es f or variabl ed e nsity time periodi c flow and h e at transfer probl e m s . .... .... ..... ...... .. .. .. . . . . ... ..... ..... A piston -c ylinder s e t up . ............. .......... . ....... ... ... ... .. Co mput e d hydrodynamic pre s sure distributions along th e x ax is fo r th e pi st on cylind e r set up . .. ... ............ .... ... . ...... ..... .. T h e rmodynamic pre s sur e in th e third c y c l e for the pi s ton Vll 5 1 8 19 2 6 39 3 9 4 2 4 2 43 4 8 48 55 58 58

PAGE 8

Fig. 3.12 Fig. 3.13 Fig. 4.1 Fig. 4.2 Fig. 4.3 Fig. 4.4 Fig. 4.5 Fig. 4.6 Fig. 4. 7 Fig. 4. 8 Fig 4.9 Fig. 4.10 Fig. 4.11 cylinder set-up . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61 Temperature distributions along the x axis for the pistoncylinder set-up . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61 PV diagram of the second test case ................................... Transient processes of pressures. (a) Thermodynamic pressure; (b) Pressure drop between pistons. Convergence history to time-periodic state ........................... Velocity fields of Case 1 63 69 71 (a) wt=O; (b) wt=1r/2; . . . . . . . . . . . . . . . . . . . . . . .. . .. 72 (c) wt=1r; (d) wt=31r/2. . . . . . . . . . . . . . . . . . . . . . . .. . 73 Hydrodynamic pressure contours of Case 1 (a) wt=O; (b) wt=1r/2; ................................................... 74 (c) wt=1r; (d) wt=31r/2. . . . .. .. . .. . . ... .. . . . . . . . .. . .. .. 75 Temperature contours of Case 1 (a) wt=O; (b) wt=1r/2; . . . . . . . . . . . . . . . . . . . . .. . . . . 77 (c) wt=1r; (d) wt=31r/2. . . . . . . . . . . . . . . . . . . . . . .. . . 78 Density contours of Case 1 (a) wt=O; (b) wt=1r/2; . . . . . . . . . . ... . . . . .. . . . . .. . . 79 (c) wt=1r; (d) wt=31r/2. . . . . . . . . . . .. . . . . .. . .. . .. . . . 80 P V diagram of Case 1 .................................................. p hase angle effects ........................................................ Refrigeration temperature effects ...................................... Effects of changing the ratio of chamber width to channel width Temperature distributions along x axis for two different {3 's (a) {3 = 0. 02 9 8; (b) {3 = 0. 2 9 8 ....................................... Vlll 81 87 88 91 94

PAGE 9

Abstract of Dissertation Presented to the Graduate School of the University of Florida in Partial Fulftlment 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 pha s e and conne c ted by narrow channels with a linear wall temperature distribution simulatin g the regenerator. Navier Stokes equations governing the variable density os c illatory flow and heat transf e r within two dimensional domains with moving boundaries ar e numerically solved by using the SIMPLE algorithm incorporated i nto 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 tran s fer problems in Stirling type devices. The pressure splitting technique separates th e hydrodynami c and thermodynamic pressures which are different by ord e rs o f magnitude for gaseous flows with heat transfer like those in reciprocating machines lX

PAGE 10

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 i11corporated 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 90 For the multichannel 1nodel, the ratio of the chamber width to the channel width also affects the cooling performance, and the optimum ratio depends on specific operation conditions. X

PAGE 11

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), alter11ative 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 1

PAGE 12

Cold Chamber 2 p 3 R egenerator : :: ;: : : :: : :: : ::: :: : ::::::: :: : ::: ::: ::: ::::::: : : : : :; : ::~: :;:: :::::::: :;: :: :::: :::::::: :: : : : '""'""' . .. .. ....... .:::-:,;,:,:::: :::::::::::::::: : : ::::::::::::::::: !=:~ :: ::: : ::~:: :-: ::: ::: ,::: :::: :::::: ::::: :: :: . ....... ............ .................... ,. . ...... .... ........... :: :: ::: ::: ::: : :::: ..... .. :-::::: ::: ::: ::: ::: : :-: ~=~ : ~:i:;: i:i :~: !:: .:t:::: : ;:;:i:~:~: i:~:;:;: I I I I I I I I I 2 / Q re j ected ......... T h Q stored 1 Q res t ored ""Q absorbed T c 4 I V H ot Chamber 1 Compression of gas at T h ; Heat rejected out 2 Gas moves th r ough R at V 2 and is precooled. 3 E x pans i on of gas a t T c ; Gas absorbs heat. 4 Gas moves through R a t V 1 and i s preheated. 1 I Fig 1.1 Stirling cycle and a three component Sti r I i ng refrige r ator configuration

PAGE 13

3 porous materials so it has an extremely high thern1al 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 c and the left piston expands. Owi11g 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

PAGE 14

4 cl1amber during the compression and all is in tl1e 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 so1newhat 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 tirne 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

PAGE 15

Left Piston Expansion Chamber 5 Cooler Regenerator :,:,:,:,:.:-:,:,:,:,::::::::::::::::: :::::::::::;::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: ..... .. ......... ................. ... i ... ................................ ............................... ........ ... ....................... .. ................. ... ........ :::::::::::::::::::::::::::::::::::::::::;:;::::::::::::::::::::::::::::::::::~:::::::::;.;::::::.:::::::::::::::::.:: ............. ................................... ........ .. ... ............................... ::::::::::::::::::::::::::::;:;::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: :::::::::::::::::::::::~::::::::::::::::::::::::::::::::::::!::::::::::::::::::::::::::::::::::::~:::::::::::::::::::: .. ............................... Heater Compression Chamber Fig. 1.2 Schematic of a five-component Stir l ing r ef r igerator Right Piston (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 Stir l ing 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 inve11tion almost two centuries ago, many different Stirling devices have been developed and found their various potential and actual applications: automotive engines s1nall electric generators, underwater power units, solar thermal conversion space power artificial heart power heat pumping, cryogenic cooling engines (cryocoolers) near ambient-temperatu r e refrigerating machines (refrigerators) The early and primary developments focused only on machines providing mechanical driving power (i.e., prime

PAGE 16

6 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.

PAGE 17

7 W11en 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 witl1 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 p r ovider for this purpose. A small Stirling engine in such devices uses radioisotopes as the heat source for a long and continuous I ife. 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

PAGE 18

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 Swee 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

PAGE 19

9 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

PAGE 20

10 temperature superconductors and their increasing applications in electronics the current need for Stirling cryocoolers is more demandin g 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 ele c tronic components that need cryogenic cooling for higher resolution and sensitivity Thes e 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 th e 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 spac e missions (especially the long-term ones) it is not feasible to carry a large amount o f 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 si z e and li g htwei g ht mechanical quietness and expected very long life without maintenance As r e ported 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 Stirlin g cryocoolers as mentioned before and to its higher efficiency than others in the temperature range of interest

PAGE 21

11 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). A c cording to Walker et al ., at near ambient temperatures the coefficient of performance o f 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 1nanned 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 lightweigh t ( 101 lbs including the heat pipes) quiet and efficient ( heat lift of 90.6 watts at 4 C with a power input of 60 watts and heat lift of 88.5 watts at 22 C with power inpu t of 70 watts). It was successfully flown on board a space shuttle in F e bruary 1994 for three days as a freezer and four days as a refrigerator. Sunpower Inc.

PAGE 22

12 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 -26 C 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 50 C 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

PAGE 23

13 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 (1D) 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 ;

PAGE 24

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 l1igher order terms Another way of treating the nonlinearity is to express the nonlinear ter1ns 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

PAGE 25

15 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

PAGE 26

16 i11 the present work is below lxl0 4 .) 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 on1y within one end chamber of a Stirling microrefrigerator with various orifice configurations and they recommended a better arrangement of the orifices. They also used tl1e 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.

PAGE 27

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 (Bl 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 c and that of the right is at Th. The isothermal surfaces are extended into part of the connecting channels to a le11gth "d" in Models A2 and B2. As indicated in the figures, the channel half width is "a", the total channel lengtl1 is "L" the half chamber width is "b". In either the left

PAGE 28

18 (a) Single-channel model: A 1 T c Th y I I I I Tlin ear C I I b C ~1 I I I t I I X I I I I I I I I I< =I IE =I E I I a I I I I I I I I I I I I I I x, ett xrigh t x,eft = -L / 2 c + Sin(wt) xright = L/2 + c + Sin(wt+
PAGE 29

(a) Multichannel model: 81 x e tt x e 11 C 0 U) a. O> C ro ::, rn C Left chamber = -L / 2 c + e sin ( (l.)t) T o b T o T lln ea r \ ........... ... -.. .-. lll~llllli~l!lll~llil~lll ......... ., ......... ~~~t~~J~;I;ti~i~~~tt :::: :::::: ::: ::: :: ; ::: :; : ;: ::;;;: ::: 1 -::::: ::: :: : :: ~;:~:: ::::: f.:: :: ::: .; ;,,,;: ~;:;:-;,;,; : 1 :-: ::-:.-: ;~~li~~~~J~lilit y i~immrt~i~~~~iif J !~l~il~l~l~~~!~l~ll~l~~~lI ll~!Il!I~illll! ................. ll!!l~!!ll~!!!l!!!f il~ll!f ill~!ll!lll!~ll!~!lll l!~!ll!I~ i .... . . . . . . .. ....... ... .,.. .. .. .. . . . ... . . . . ......... .... ....... = : : :-:. ::-: : ::: :; : :-:-:-:-::-:l~!l!llilllli!l~lj~j!j~ll~~ : ... ... . , . .. . ,., ... 1 : :;:;:: :::: :;: ~:::: :: ::: ::: :: ::: : I .............. . .............. : :::::::::::~:::::::::::::::::::::? :;:: : :: : :-:: ........ . . . . .. . ... . . . . -:. .. .... . . L ; : : :: :.: : ::: ::.:: :::: : : ::::: : : : : :=:: ::: :: : : ::::::: :: : : :::::: :::: ... ...... .... .... . ::;::: ::: ::::: ::: : : : : : :: :: :;: : :::: :: : : ::: :;: ;: : : : : ::::: :: :: :: ::::: :: :: : :;:;:;::::::::::::::;:;:;:;::::: : :: ., .. . .... : :: : :: ;:: : :: ;:::: =:~ -: : :-: .: :; :: ;-;~: : 1 T h a T h C 0 U) a. O> C ro ::, rn C Right chamber X rig ht X rfg ht X U2 + c + e sin ( wt+

C ro ::, rn C Left chamber -U2 c + e sin(wt) T c b T c T llnear :,. -:-:-:-:-:Itrrm rrrrr : ; : : :: : ::: :::::::: :-::::-:-:-:-:, ::::~: '. :'.: '. :i:'.: i: :i: i :/ I/ii t::i:: :-:-:-:-:-:-:-:-:ti\!Ii\ :::;::::::::::::::: :-:::-:-:-:::-: :i::I::::t:l: t??{} :::::::::: : ::::::: : !/:III .-.-....: ,-,-.-. :::;::: ;!!: ::::::: :-::-:-:-:.;.:.:lilil!!ll!l!llli :-:-:-:-:-:-:-:-:. .. .-. . :: :: ::::::: :: ::: :: ::::::::::=:::::: : ~Ittf : :-:-:-:-:-:!: : ~: :: ::: ; :::;: :: l!!l!!!liii!!!! ::::::::::::::::: JI~~ttf ===-==-:-:-:: : :-:-:-:-:-:-:, L : ~lllt tl~I~~! 11:1. d (a) Model Bl; (b) Model B2 y -------:-:-:-:-:-:-:-:, !:!i!ltl!l!ll:!i j ::::::::::::::::~ :-:-:-:-:::-:-: ::::::::::::::::: : :: : ::::: :::: :: :: . . .. -. : : :: :-:: : :::-:-:-:-:: ::::::::::::::::: :-:: -:-:::: :: : : : : ::::;:: ::::: ft/:;}: 1111 '. '. :-: '.'. '.'.: T h C: 0 ... x right rn a. O> C ro ::, en C ll~!~~!!l!l1 ~t~=t~:~i 1.: X !lilIIil!liI~ : ~~@J~il l ~lll~~ll!lll : ;-: : : : : : ,;-; iII:II :;::::::::::::::::: ffttI -:-:-:-::-:-:... trtrt ::::::::::::::::::: i:J:III :!1111111 1 1 1 1:l!l!I :::::::::'.:'.:::::~: T h JI~r~~t ;: : :; : := ::;:;:::: tlll :~: ;: ~::: :: ::::: -::: . i)?\?i ::::::::::::: : :::: )}{ff : ((}tf :-:: Right chamber x ri9ht U2 + c + e sin(wt+cp o ) \0

PAGE 30

20 or the right chamber, a piston oscillates sinusoidally with an amplitude of "e" and an angular frequency of "w". The piston motion in the right (hot) chamber leads that in the left by a phase angle "cp 0 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 ther e 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, i11 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 pr ess ure splitting tech11ique. Since the domain of interest is bounded by moving pistons a movin g grid te c hnique as used by Demirdzic and Perie (1990) is also incorporat e d 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 devi ce s.

PAGE 31

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 ex = a where w is the oscillation angular frequency and v the kinetic viscosity.) This dime11sionless 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 90 The effects of varying the temperature ratios between expansio11 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 th e 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 arrangeme11t of tl1is 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.

PAGE 32

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 Stirli11g 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 22

PAGE 33

23 d f f pdV + p(V VJftdS dt V(t) S(t) = 0 (2.1) Conservation of momentum, d f pVdV + f pV(V VJftdS = f (-pi+i)ftdS, (2.2) dt V(t) S(t) S(t) Conservation of energy d f f f f~ pEdV + pE(V VJftdS =qftdS+ V ( pl+-'t).ftdS. dt \f(t ) S(t ) S(t) S(t) (2.3) In the above equations v(t) represents the control volume S(t) the moving control surface y the velocity of the control surface and ft the outward unit normal to the b 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 i the deviatoric stress tensor and q the conduction heat flux vector. Here the assumption of N ewto nian 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 tl1e velocity of the control surface is zero, then the Eulerian approach is taken. Otherwise the viewpoint is neither Lagrangian nor E ulerian as happens in this study. Since the body force is negligible potential energy is negligible too. The total e n ergy is then given as E = e 1 + VV 2 (2. 4)

PAGE 34

24 where _!y.y 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 = C v T. ( 2.5) Here c v is the specific heat at constant volume and T the absolute temperature. According to Stokes hypothesis the deviatoric stress tensor is 2 ~ -(VV)I 3 ( 2.6 ) where the supers c ript "T" represents the transposition of a tensor and 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 q = k VT, ( 2.8 ) wher e k is the coefficient of thermal conductivity. Viscosity is determined by using Sutherland s formula = (2 .9 ) in which T r is some reference temperature r the viscosity of the fluid at the reference te1nperature and T a constant that depends on the fluid. Thermal conductivity k can be determined from viscosity and specific heat at

PAGE 35

constant pressure 25 k = cP Pr ( 2.10 ) where Pr is the Prandtl number which is usually a co11stant for a specific fluid. For a perfect gas and 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 (2.11 ) (2.12 ) The pressure p(x y z t) consists of hydrodynamic pressure denoted by ph and thermodynamic pressure denoted by p t, i.e. p(x y z t) = Ph + p t (2.13 ) The hydrodynamic pressure is related to the inertial and viscous forces and is usually a function of both time and space, so (2.14) a11d the thermodynamic pressure is related to the volumetric changes and temperature variations. At v e ry low Mach numbers the thermodynamic pressure is a function of t ime only 1. e. (2.15) Consequently

PAGE 36

In numerically solving the momentum equation (2. 2) pressure gradient terms like ap/ax need to be approximated. As a demonstration consider a 1 D problem where the pressure gradient is dp / dx. As shown in Fig. 2.1, the pressure gradient at xi+ 1 12 estimated by where 26 (2.16) p j P i+1 I I I X x I X ; +1 12 x i+1 Fig 2.1 Illustrative lD discretization (xi + x i+i )/2, denoted as (dp/dx)i + 1 12 is (2.17) (2.18) (2 .19) Fro1n the above discussion p t is independent of x and is several orders lar ger 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 nu1nerical solutions are erroneous (Golub and Ortega 1992). Although using double-precision arithmetic in

PAGE 37

27 computations can increase the accuracy, a price of more computer memory storage and C PU 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 accura cy and efficiency could be increased The rationale follows. Substituting Eqs. (2. 18 ) and (2. 19 ) into Eq. (2 .1 7) yields ---dph ( dx )i + l/2 (2 20) Then the problem of catastrophic cancellation does not exist anymore and therefore the c omputing ac c uracy is enhanced without using double precision arithmetic. Th e c omputer CPU time memory as well as storage are saved. Now the question is how to determine p h and p t Such a m et hod is developed in the prese11t s tudy First consider how to determin e the hydrodynami c pressur e ph. I t can be seen from Eq. (2 20) that th e pressure p appearing in the momentum equation (2.2) can now be replaced by the hydrodynamic pr ess ure p h and therefore p h is solved by using the pressure based SIMPL E algorithm which will be dis c ussed 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 pl. It can be determined by using the "' P e rsonal communication 1994

PAGE 38

28 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 +pt p RT RT (2. 21) Tl1e total mass of the entire system is known at any time instant, and Total Mass = J p(x ,y,z, t)d'v' = J p1xx,y,z,t) +pt(t) d'v'. ~(t) ~(t) R T{x,y,z, t) (2.22) Therefore -1 -[Total Mass J Ph d'v'] f 1 dV ~RT ~(t) RT (2.23) 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

PAGE 39

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 = Xi e ft 0 ::s y ::s b left piston wall no slip: U = U1 eft = -e wSin(wt) V = 0 x 1 eft < x < -L / 2 y = b, imaginary boundary symmetry: au / ay = 0 v = O. x = L / 2 a < y ::s b vertical wall of the left chamber no slip: U = 0 V = 0. L / 2 ::s x ::s 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 < x right, y = b imaginary boundary symmetry: au l ay = 0 v = O. x = x rig h t, 0 ::s y ::s b, right piston, no slip: u = u rig h t = ewsin(wt+<,0 0 ), v = 0. X 1eft < x < x right, y = 0 symmetry: au / ay = 0 v = O The temperature boundary conditions are as follows: X = X1 eft 0 ::S Y ::S b: X 1eft < X < L / 2 Y =b: aT ax = 0. aT 1 ay = o. (2.24) (2 25)

PAGE 40

x = L / 2 a < y :s b: L / 2 :s x < L / 2+d, y=a: -L/2+ d :s x :s L / 2-d, y=a: L / 2-d < x :s L / 2, y=a: x = L / 2 a < y :s b: L / 2 < X < X rig ht y = b: x = Xri g h t 0 :s y :s b: X1 eft < X < Xright, Y = 0: 30 T = T c + (Th T c ) (x + L / 2 d) /( L 2d). aT1ay = o. aT!ax = 0. aT1ay = o. Here the boundary conditions for velocity and temperature only include two types, Dirichlet and Neumann and all the Neuman11 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

PAGE 41

31 nondimensionalized as follows. First reference quantities are properly chosen: Reference length: L r = a (channel half width) ; Reference time : tr = 1 / w (w = angular frequency of the piston oscillation) ; Reference velocity: Ur = L r/ t r = aw. The reference velocity is defined in terms of the reference length and reference time and cannot be arbitrarily chosen ; Reference temperature: Tr. It should be chosen such that it is as close as possible to the temporal and spatial average temperature ; Reference density: P r Like Tr the reference density is also chosen so that it is as close as possible to the temporal and spatial average density ; Reference viscosity: r = r) T11en th e nondim e nsional quantities which are denoted by bars can be defined as x 1 a r t t = e e 2' Ur E U, 1 aw E 1~-=. = e+ VV u2 2 r p = p' Pr 't 't Ur r L r T T = T' r Sp e cial care should be take11 in nondimensionalizing the hydrodynamic and thermodynamic pressures because of their inherent different scales. The corresponding nondimensional quantities are defined as

PAGE 42

32 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 df f ...... p d'v' + p (V Vb) ft dS dt 'v(t) S(t) = 0 (2.26) Conservation of momentum f phftdS + 1 2 f i ftdS S(t) a S(t) (2.27) Conservation of energy d f pTdV + f pT(V VJftdS = y 2 fV(T)ftdS dt 'v(t) S(t) Pr S(t) (y-l)fpTViidS + y(y-!)M 2 fV'riidS (2.28) sro a sro + y(y l)M 2 [ d J .!pVVdV + f .!pVV(V Vb)ftdS], dt 'v(t) 2 S(t) 2 Equation of state

PAGE 43

33 (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: Ratio of specific heats: Mach number: = (reference speed) I (speed of sound at T r), Womersley number: (2 .30 ) (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 latt e r is defined as Va = a 2 w/11 = cx 2 = L U I v r r r r (2 .32 ) It is called the Valensi nu1nber 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 c onditions are b /a, c/a, d / a e/a, L/a, T c / Tr Th / T r and 'P o Since T is used as the dependent variable to be solved for in the energy equation (2.2 8) quantities p e and consequently E are all expressed in terms of T there. The sec ond term including the minus sign "-", on the right hand side (RHS) of the equation

PAGE 44

34 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 th is study the dissipation term and the terms related to kinetic energy are neglected because they all include the factor M 2 and M here is always smaller than 10 4 Therefore, the energy equation is simplified without loss of computational accuracy to d f pTd\i + f pT(V-VJftdS = y 2 f V{T)fldS-(y 1) f pTV ftdS. dt...,. ( t) set ) Pra set> S(t) ( 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 f pdV + f p(V Vb)11dS = f (rV)ftdS + f Cd\f dt 'v'(t ) S(t) S(t ) v'(t) (2.34) where r is the diffusion coefficient r the source per unit volume. Each specific integral equation corresponds to a specific set of rand r 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 ; the second term on the LHS is the rate of efflux across the control surface

PAGE 45

35 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 a11d energy equations. Table 2.1 Terms in the governing equation (2.34) r r Continuity 1 0 0 aph 1 a au 2 a av a av Ol. 2 ax +3 ax(ax)-3 ax(ay)+ ay(ax) u-Momentum u aph 1 a av 2 a au a au Ol. 2 ay +3 ay( ay) 3 ay( ax) + ax( ay) v-Momentum V '}' -(y-l)[ a(puT) + a(pvT)] Pr 0t. 2 ax ay Energy T 2.6 Dimensional Analysis The physical quantities of rnuch 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

PAGE 46

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 Lr tr Pr and T r 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 Lr 5 tr 3 PrTr 0 i e ., p~ 5 w 3 From the IT theorem (Fox and McDonald 1978) we have P I (p~ 5 w 3 ) = .9{Pr 'Y M, a,


PAGE 47

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 Demirdzic and Perie (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 witl1 grids into a system of subdomains or cells. At every time step new 37

PAGE 48

38 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 = 1r / 4 and 1r, 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 1 ine segment parallel to the x-axis starting at x = x 1 and ending at x = xN with the first and last grid spacings being specified as llx 0 and dx r respectively (Fig. 3.2). Grid points are distributed according to the hyperbolic tangent distribution given by Thompson et al. (1985) as follows. First define (3 .1 ) and H = ( 3.2 ) Then solve one of the following nonlinear equations for the parameter o by using the Newton Raphson method (Hornbeck 1982):

PAGE 49

(a) (b) 39 11111111 11111111 111111 11 ill I ll llllllllllllllllt 1111 1 111: 11 1111111111111111111 11111t1 I ll ll ll llllllllllllll lDIDI :::::;:~:::::~::::::::::::::::::~~::~;*::::;:~::::;::::~:::::;:~::::~:::::~:::;:::::;:::::;::~~::::::::::;:::::~~: :0-:''' .... 11 I 11 111 111111111 1 11 II 1111 111 1111111 11 I Ill I Jilli II '111111 111111111 11 1111111 111 1!111111 I I II 11111111 1 111 mr : I lllll lllllllll !l lllllll lll : III I UI IIIIIIIIIIII III III III : I 11 11 11 l j l lllllllllllll llll i 1 111111 111111111 11 111 1111 i 11 11 11 111111111111111111 111 1 111 1 1111 1111 1111 111 11 1111 : 1111111 11111111111 11 11111 111 I llllllll lllllllll ll llll llll : I 1111111111 1111111111111 1111 l 'l lllUI III I IIIIIII 1 1 1 11 1111 ,1 I HMI I U11J H1t111 11 f ll!t1 ,,,111m,~1, 'ti U"I' : ~:;-~:;;;m!::~:i!:::: 11111 1111111111111111 1 11 m I : 1 ll j l1 I 1 l ll ll l l 11l I l l 1 11 1 11 1 1 i 1l l ll ll l llllll, I I I l l ll l l 1 i i ...... , .. ..... , ......... ..... ' ........ ........ , ,: ,: '\: : .... .... , .............. .: ... .. ' .. .. ::: .. ... ... ~:i:~:~:~:;:~:~:~:i~~:~~:;:;:~:i:~~:~:~~~=~:i:i::$:::::::::::::::::~::::::::::::::::::::::::::::;::::::f::;::::i::;~~:; ::::::::::::::::::::::::;:::::::::::::::::::::::::::~:::::::::::::::::;:::::::::::::::::::::::::::::::::::;::, 111 1 11 1111111111 1 11 11 11 11 i: 11111 1111 111 1 11 11111 : llllll mmm I 111 11111, ,: IDP IIIIUUilllllll lll l ll : ... .. ... .............. . . . . . . .. ............... ......... IIUIJ IIUU IIII IIIII UIIU IU l! ff fltltllllll!I IIUIOI : 1uuu m111 1 11 1 u1 1r1nl1+ u,no~ .. ,f"l!"ffi''"'" :!: :,,~ 1 i ~n. ~~ti: : -:,i:,:;:I !,!,,,: ::::,,: JJ -'_......_ _ .................................. ~~ gsi!j ~ : : .... ::1gii;:,:L,::1:i.:; : i:i; ;: i; :ii ,;0,, 1 ; 1 ; ... ;; '. :I== = a:==__:_ :::1:;:lUil;U:1 !!::::.: t ;===-= ,.._.. ---------------""". ........ .. ,t .. ,t:n,t.n.nr1 1 .. 1 .. j" =-= -..... :::::::::!:n:11:11:1:r::1, r --=5=:i~ ;;;m11ajw11; ,::;;, ,, .. ====-:::::::::: ::, !,.1,,.:: I: ~=====~~ ---=i .... ,.,,1 .......... ,.,,,1 =-a = ........................ -=-!!;:,. :.,~,::;::,:;;:;:;;1.;i .. ; == ;I,:;: ,::::::::::::111 I === ==. ,..,,.,,,,,,;;,,,,,1 -------~-....... .... ' .. .......... J ;,; ........... ,. :: ----= ==~ WHiHHHHH:i:::: i 11e= ;:;;~ = --=~ aaa, 1 ::::;frnmff1:1 i-i == 9=:sihii : : -:,::::::::::1:: .: ---: ...... .. r .. .. ;:a::: : m:;:;r,:f:p.i:m::. : 1::: : a: :: :m.-:: : :::::;:::. : m u 1 m lnj11 u 111111 un,1 1 11 11111 "'11 1 '"1'' 111111111 1111 I I 'I 11111 11111111 1 111111111 11" 1 1111 1111111 1 I II 11111 t J1I II 111111 Ill' 11 1111, 1 i i l l l lll ll1l11 l 1 ll 1 1 '1I I I U l lllllllllllllllll II l[I 1 1 11 11111111111111 11111 1' 111 1111111111111 1 u 111111 1 111 1 111111111 1 11 111 1 11 1 1111111111111111111 1111111 II II II UIIIIIIIIII I I II III I 1111 II Ill 111111111 111111111 : :::::!::1~1 ::.u: !I ::~ I ,,.,..,,tH ttn ; j"I "' 1,nu11tu lfU U I l I ,uj : l:I I Hl ltl U I UII II f ti II) lilllllllllll lll l l ll I I Iii : t:111111 111 111111 I 1 1 1 Ill ... .. ..... ..... l:IIIII IUIIIIIIUI 11 1111 1111n 11 1 m 11111111 111 11 111: Ill Ill Ill llll 11 1 11111111 111 i "i; .................................. ... :!:::::::::::!::!::::::::::~::~::::::::::::::::::::::::::::::::~::::;:!:::::::::;:::::::::::;::::::::::::::::::::::::~:~::: ........................ ::::::::::::::::::::::::::::::::::;::::::::::::::::::~:::::;::::,:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: ... .. .~ ............................. ... ... . . . . . . . . .. . .. ................... ....... :::::::::::::::::::::~:::;:;:~::::::::::::::::::::::::::::::;::::::::::::::~:::~::::::::::::~:::::::~::::::::::::::::::~: ........................................ .. .... ........ 111111 1111 1 1111111111111 111 I 1:111111 11 11 11 11111 11 1111 :;::::::::::~:::::::: : Jllllllll lllllllll 1 1111 1111 i 1 11n11n 111111111 1 1111 1111 : 1 11111111 111111 111 11 111111, ... 1 1111111111 11111111 1 1 111111 1 111111 111 11 111 1 1111111 1111 .... ... ::::::::::;::::::::::::~:::!::::::::::::::::::::::;:::::::::::::::::::::::;:::::;:::;::::::::::::::::::::::::::::::::::::::::: 11 11 11 111 111111 Ill II i ll ::::::~:::::::::;:::::::::::::::::::::::::::::::::;::::::::::::::::::::::::::::::::;::~:::~:::::::;:::::::::::::::::::::::::: ~iJ~l~~~~~it;;!;~;J~if~J~~i~~!J;Ji~)l;~~~~f 1il]i~;~f ,~ ........... .. ........................................ .. ........ .. .. ...... .. ....... ... ... ... :;:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::;:::::::;:::::::~=~:;:::::::::::::::::::;::::::::::::;::::::::: ... ... ~imtt~~~~J~~tJ~~~i~~~~jJttJtttt~r~tt~I];f f {~t~;l~~t~;~ !::::::~::::::::::::::;;:~:::::::::::::::::::::::::::::::::::::::::::::::::::::~:.::::::~:::~::::::::::.::::::::,::::,:f~:~:i .................. . . . .. . . .. .. . . .. . . . . . :::;:::::!::i:::::::~::::::::::~::::::::::::::::::::::::~:::::::::::~::::::::;~:~~:~~::;:::::::::::::::::;:::;:;::::::::: ::::::::::::~::::::!:::::::~:::::::::::::::::::::::::::~:::::::;::::::;::::::::::::::::::::::::::::::::;:::::~::::::::::: ~J~l~~~;J~~ltflltttlt~~ttf(;~;t(;!tiJftt:ttt~:i:~:~:\ ;::::::::::::::::::::::::::::::;:::::~::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::;::::::::::::::::::::: ........... ............. ............................ : . ::. ;. ,:. : :: :.: :. :: :: ; ,:.,.: @f i;~~Jt r tlJttJttrtitIJtm~~~Jm -+-f--l--+-+-+-f--f--lf--lf--lf--lf--lf--lf--lf--1--1--+--+-+-!--lf--+-+-+--+-+--+-+--t: : .................. : .... : :.: : :. ::: : :. =.: : . Fig. 3.1 d)( o 0 0 X 1 X 2 .. .. Two typical grids. (a) wt 1r / 4 (b) wt 0 0 0 X 3 x I Fig 3.2 1D grids 7r d)(f 0 0 X

PAGE 50

sinho 0 sino 0 1 G 1 G 40 for G < 1 for G > 1 The i th grid is distributed at where g + --1 --{X X ) xi = x l H + ( 1 H) gi N 1 1 gi = { 1 + 2 tanh[ o c i i l ) l N 1 2 ------} tanh! 2 (3 3) (3.4) ( 3 5 ) (3.6) For each chamber block, either xN or x 1 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 tange11t 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

PAGE 51

41 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 si1nilarly v control volumes also shift in they 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 Sand 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 g ) e-at the east interface and (u g ) w-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 over cell P is

PAGE 52

0 0 0 Fig. 3.3 0 0 i:i:~:f~:1t:i:~~ ::: ::: ::~::: *: :::::: ,: ,: :-: :::::;::::::: ... '. h ;: : ;::::;;:;~::::: :: : }{~11]} 42 0 0 V u 0 0 0 Contro l vo l umes fo r scalar q u a n tities 0 p p 0 :~!~~:li~~:j:~il~lil~II!t . . . ....... .. .. . . ' ....... '' . . .. . . ....... ... . H ... I . ....... . i~!:l~l:~~~~il~l~l~~l~!:!:l~ll ;i:~;~:~i~;;~~I~;~~~;1~~~;;~1 [ ~::::1~;l~~!~~~lj~~~i~~ v =~t1itl~tl r -~ : .. ~.: i~mwmr~ . . . . . . . .. . . . . .... ., ... ........ :-:-: ::: : :::: ::: :~: :: :=:: :::: :: :: :: lltli1Jt}tlt / IlIJIIIlll T ~::;::~:;,:,:::,:,:-:~-:=:::::::::,:,:-:::::-:::: 0 V c ontr ol lume V O u co ntr o l vo lum e F i g. 3.4 Staggered u and v co n tro l vo lu mes

PAGE 53

43 d + f Sp( t ) ftdS + J 'v'p(t) ( 3.7 ) where the s u bscript p d e no tes t h e qu a n tity of ce ll P 1---~xw I ~Xp ~ E I N YN 0 0 0 n v n I~iil~l~~~l~ t ~}~! i !~! i~i 8 y n ...... ............ .. ....... :::::::::::::::::::::::::::~:::::::::::::::: :::::-:;::;::: ::::::::::::::=:::=:::::::::::::::::::::::::::: .... .. ....... .. .... ...... ........ .... ........ . . ......... .. ............... I I ...... .. ....... '..... : . ... ..... ............. : :-:-:-:,:-:-:-:-: :-. ,:,:,:::::::::::::. :: ~ : :: e :W 1 Yp w e 0 E 0 . .......... .. ....... -:,.,,:-:,:,:-:-:-:-:,:-:-:, YI ........... ..... ,,..,., ... 1 .v,v.,., .,_.,. ,., .......... .. .. ,. ( u g ) e :,.,,:.;.: . : ,: :;::: :::::::::::::::::;::::~::::::::::::::::::::::::::::::::::::::::::::::: .... ,, u ... :W: ,:,:-:-:,:. -::. ' ... ,,,,,,,:,:,:-:-1,:-:,:-:,,:-:-:-:-:-:,:,:,:-:,:,:,:,:-:-::::: ...... .. .. .. .... .... ......... :;:::::::::::::::::::::::::;: : ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::,:::::,::; .. i :-:.:,:-:-:::-:,: :-:,:, :,:,:,:::::: ,. .. ... '.' ....... ....... : 8 y s .......................... .......... ............ . . . .. .............. ''' , . . ........... .. ..... :-::-: s Y s 0 0 0 1 s I I F i g. 3 5 Schematic of a genera l ce ll for at p oi n t P Consi d er the ra t e of change te rm f ir st. I n r eference t o F i g. 3 5 th i s term can be discretized as d J dt 'v'p(t) pd V ~t ( 3.8) where ~t is the time step size V p the ce ll vo lum e an d t h e s up erscri p t II O II re p resen t s a q u antity at the l ast time l eve l Here o nl y a 2D case i s cons id e r e d so

PAGE 54

44 (3.9) and consequently (3 10) The convection term is discretized as f p(V Vb) fidS = Pee(u ug) e ayp Pww(u ug)wayp Sp(Q (3.11) + Pnn V n axp P s s v s axp Define mass fluxes through the cell faces F O = (pv) 0 dXp F s = (pv) s Lixp then (3 12) (3.13) f p(j)(V-Vb)fidS = Fee Fww + Fnn Fss (3.14) Sp(t ) The diffusion term is then discretized as (3.15) The first derivatives in the above expressions are approximated by first order finite differencing, e.g .,

PAGE 55

45 (3.16) Thus one has ayP ayP RHS of Eq.(3.15) = re -( B P) rwo (p-w) oxe dxP dxP (3. 17 ) + rn O ( N p) rs O ( p s) Yn Y s = De ( E P) D w ( P w) + D n ( N p) D s ( p S ) where the diffusion fluxes are defined as The discretization of the source term is straightforward, J ,dv = vp,p = Ap'p 'v'p(t) (3 18) (3.19) Substituting all the discretized terms into the integral equation (3. 7) gives ,h o o ,ho ppAp'f'p ppAP'f'P + F ,h F ,h + F ,h F ,h = at e'f'e w'f'w n'f'n s 'f' s (3.20) Noti ce that the fully implicit first order backward Euler scheme is used. In the above equation the val u es at cell faces i. e e, n w and s have to be expressed in te rm s of the values at cell centers i.e. 8 N w and s since the latter ar e regarded as unknowns to be solved. Convection schemes are needed to do this (f or

PAGE 56

46 details see Patankar 1980 Shyy 1994). Applying a convection scheme to the above equation and then collecting terms gives the algebraic equation for in the final form: (3 21 ) where for compactness the subscript "nb" denotes evaluations at the neighboring points of point P and the coefficients are defined by a E = D e C(I Pec e l) + max{-F e' O} a w= D w C(IPec w l ) + max{F w, O}, a N = DnC(IPecnl ) + max{-Fn O} a s = D s C ( I Pee s I) + max{F s, O} 0 ap = (pA)~ at (3.22) ( 3 23) (3.24) (3 25) ( 3 26 ) (3 27 ) ( 3 28) The symbol "Pee ' in these expressions is the cell Peclet number which is the ratio of the convection flu x to the dif f usion flux. The Peclet number at the east face Pec e, is defined as ( 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 ( I Pee I ) = max{O l ( 1 / 2) I Pee I}. ( 3.30 )

PAGE 57

47 This scl1eme is essentially central difference for I Pec e I 2 and switches to first order upwind for I Pec e I > 2. For the problem under consideration the source term is only linear in at mo s t and so for simplicity it is just lumped into the term B (E q. (3.26)). Special care should be taken in the discretization of boundary cells. As discussed befor e there ar e 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 s is given In this case, let point S coincide with point s then (3.31) H e r e also note the special definition of oy s 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 a<1> ay is specified is taken as an example to demonstrate the boundary treatment. As shown in Fig. 3. 7 an imaginary cell S i s added to th e south of cell P such that and so Ay 8 = Ay p, oy 5 =Ay p Applying the boundary condition with central difference approximation yields (3. 32 ) (3. 33 )

PAGE 58

48 I I I I I!! ru


PAGE 59

49 (3.34) then (3.35) 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 s = s for such special cases. Throughout the discretization process, linear interpolation is frequently used to evaluate values whenever they are not readily available, except 's themselves at cell faces which are determined by convection schemes instead. For instance, in evaluating the mass flux Fn = PnV 0 ~Xp (refer to Fig. 3.7) if p 0 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 involveme11t 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) a11d (3.13)). Since the grid velocity in this work is neither zero 11or tl1e 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

PAGE 60

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 e quation 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 , (pt) (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 f inal form of the discretized u momentum equation for an u -c ontrol volume Pis given here. (3.36) 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

PAGE 61

51 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, (3.37) where Wu is an underrelaxation parameter (0 < wu < 1). Owing to the nonlinear coupling 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 (3.38) 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 co ntinuity equation because the solution is based on a guessed hydrodynamic pressure field (ph) There fore the d iscretizing the continuity equation over a pressure control vo lume at point P (refer to Fig. 3.5) yields {p*A)p-(pA)~ +F F p F* = ( 0) At e w + n s Sm t,. (3 39) The pressure needs to be corrected so that the resultant velocity field satisfies the co ntinuity equation. If the pressure correction is denoted by p ', then from the momentum equations (3.36) and (3.38), the corresponding velocity corrections at tl1e i nterface s of the pressure-control volume are given by

PAGE 62

and I u = e 52 (3.40) (3.41) Since there is a half cell shift between the velocity and pressure -c ontrol volumes, an index shift is also applied when deriving the above equations from the momentum e quations. The coeff icient~ in Eqs. (3.40) differ from those in (3.41) because the u and v-control volumes are not collocated. The corresponding flux corrections ar e (3 .42 ) It is required that the corrected pressure and velocity satisfy (p A)p (pA)p 0 , 1 -----+ F* F +F* F* +F F F F = 0 dt e w n s e w + n s (3. 43 ) Substituting Eqs. ( 3.42) into Eq. (3.43) and collecti11g terms give the pressure correction e quation

PAGE 63

53 I = E bniJ)nb Sm (3.44) nb where b 's are the coefficients, and s"'m 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 consideratio11. 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 w' in Eqs. (3.42) is zero. The F w 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 (3 .45) where underrelaxation is used again to maintain convergence (0 < wP < 1). The correspondi11g 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 ' v=v +v'. (3.46) (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

PAGE 64

54 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 pt 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=O, the left piston velocity is zero (Eq. (2.24)) but the right piston velocity is not if


PAGE 65

55 L1 0 In :>ut Pr Re =a M,~ ', Geometr1 I.C. & B.C. at t= Q) E C ..---------,;~ t = t + flt Generate Grid. U :>date B.C. u-Momentum u~=-----, v-Momentum v Co ntinuit', Ee. :> U :>date P n u, v ' Ener~ { Ee. T U odate P t Ee. of State D f Conver: ed at t? ___ N_o ____ Yes No .___ ____ Ti me-P erio die? f Yes ..... I I Q) E ..... ..... Q) ..... II,... Q) ..... F i g. 3. 8 S ol u t ion pr oce dur es fo r v ari a bl ed e nsi ty t im ep e ri o d ic flow an d h eat tr an sfe r probl e m s

PAGE 66

56 A quantity cp is defined by the following expression, 'P = 'P o [ 1 -cos(2wt)] 2 0 ~wt
PAGE 67

57 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 a11d the code developed were validated for variable-density time-periodic flows and heat transfer in regions with moving boundaries. Tl1e 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 p 0 = 1 std atm (=1.013x1Q5 N / m 2 ) temperature T 0 = 288.15 K, and density p 0 =1.225 kg / m 3 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 X pt s o n e 2 [ cos(wt) -1] ecos(wt) 0 s wt< 1r 1r s wt. (3 51) The oscillation frequency is f= 1 Hz (w=21rf =21r rad/s). The geometric parameters are a = lx10 3 m e = la and L = 3a. The nondimensional parameters are -y=l.4 Pr=0.72 ~= 0.6559 and M = l.846xl0 5 It is known from thermodynamics that for such small dimensions and slow oscillations the thermodynamic quantities should be uniform within tl1e cylinder and the quantities can be determined analytically by the adiabatic relationships:

PAGE 68

15 I C\I I ,........_ 8 en ........... a. I I 10 -I I ,........_ -+,J .. 0 .. en LO T""" ........... 5 .c Cl.. I ,........_ -+,J .. 0 .. >< ........... .c Cl.. I I 0 -1 58 y ~ Insulating wall I I I I ... I I :: I : ..., ___ ,.... -L-----+.: :; l 0 I I I I I I I I I I I I I a :,:,:,:-::: X plston = COS(wt) Fig. 3.9 A piston -c ylinder set up No.1 Single precision, Without pressure-splitting No 2 Double precision, Without pressure-splitting No 3 Single precision, With pressure-splitting f ... ... IHI lHI lHI .,_ .,_ f I I I I I w t:: n/ 4, No. 1 I I I I I I I \ I I ,, I I II I I 1, I I I II ,, I I 11 I Symbol w t Computation I ., II I 11 I I I I I -xrt/ 4 No 1 I I I I I I I I I I I I I I I I I I I I + 2 7r/ 3 No. 1 I I I I I I I I I I I I I n/ 4 No. 2 or3 I I I I I I I I I 0 I I 1 1 I I I I I I I /),. 2 rt/ 3 No. 2 or3 I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I JI I 11 11 1' 11 I / I I w t::2 7r/ 3 I ' A No. 1 I ........................... w t::2 7r/ 3 ooo ooo ooo No. 2 or3 o ooo 0 000000 w t:: n/ 4, No. 2 or 3 0 1 x / a 00000000000,0 2 Fig. 3.10 Computed hydrodynamic pressure distributions along the x axis for the piston -cy linder set-up 3

PAGE 69

59 t \/ T y p t = ( o) y = ( p ) y = (-) y -1 Po \;/ Po To (3.52) 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)=(l.5a,O), therefore the pressure is the hydrodynamic pressure ph. As indicated by the figure at time wt=1r / 4, the hydrodynamic pressure from Computation No. 1 oscillates wildly along the x axis but at wt=21r / 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 observatio11 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

PAGE 70

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( t) in the third oscillation cycle is shown in Fig. 3.11. To help compare the magnitudes of ph and pt, 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 tl1er1nodynamic 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 i s presented by the horizontal line in the figure. The temperature from Cornputation 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 -s plitting technique not only gives tl1e 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. Geo1netric parameters are (refer to Fig. 1.3a): a= lxlQ3 m b=4a c=3a, e=2a and L=6a The wall temperature (including the piston s) is kept at

PAGE 71

61 3.Sx 10 9 I I 3.0x10 9 C\I 3 (U a. 2.5x10 9 I I --+-' 0.. 2 0x10 9 1.5x10 9 0 100 200 300 wt (degrees) Fig. 3.11 Thermo d ynamic pressure in the third cyc l e for the piston -cy linder set-up 0.98 0.96 0 I......... +-' 0.94 0 X .......... I0.92 0.90 -1 Symbol w t Computation X n/ 4 No. 1 + 2 '1t/ 3 No.1 0 n/ 4 No. 2 or 3 2 n/ 3 No. 2 or3 Analytical No.1 S i ngle precision W i thout pressure-splitting No.2 Double precision Without pressure-splitting No.3 Single precis i on, W i thout pressure-splitting Analytical ; No .2, or 3 + N 1 + + + o 0 1 x/a 2 3 Fig. 3.12 Temperature distributions along the x a.xis for the piston-cylinder set-up

PAGE 72

62 constant T w =288.15 K (the r efore T c =Th=T w ) Pistons oscillate sinusoidally with the right piston leading the left by 90 in phase (<,0 0 =90) with an oscillation frequency of f = 2.32 Hz. Air is the working medium. At t = O p 0 =l atm, T 0 = 288.15 K p 0 = 1.225 kg / m 3 The corresponding nondimensionaJ parameters are -y=l.4 Pr=0.72 cx=l and M = 4.29xl0 5 In this calculation, the pressure splitting algorithm is used. For the computational domain which is only the upper half geometry a 9lx31 grid system is used. The grid distribution is: 31x31 in the left chamber, 31xll in the channel and 31x31 in the right chamber. The time step is chosen s uch that each cyc le is discretized into 360 uniform intervals i.e. At = l / (360f) = l.2xl0 3 sec. It takes about four cycles for the solution to reach the tirne periodic state The P V diagram for this ca l culation is shown in Fig. 3.13. This figure shows that the minimum pressure is P mwl p 0 =0. 7099 and the maximum is Pma x / p 0 = 1. 6420 which can be checked by examining two extreme conditions namely isothermal and adiabatic. For an isothermal condition p / p o = v o/ v so For an adiabatic condition p / p o = ( V o/ V ) "'' so P mwl P o = ( V 0 / V m aJ"' = 0.6389 and pfll1U( / p 0 = ( V 0 / V mw> "' = 1.9402 The numerical maximum should be lar ge r than the adiabatic counterpart and less than the isothermal one whereas the numerical minimum shou ld be larger than the isothermal one

PAGE 73

1.4 0 c.. -1.2 c.. 1.0 0 .8 0.6 0 8 63 a= 1 mm, b=4a, C=3a,e-=2a L = 6a, f=2.32Hz, cp 0 = 90 Air Pr = 1atm, T r =288.15K T c = T h = T r I O;= 1 1 2 1. 4 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 tl1an to the adiabatic ones because of the boundary condition of isothermal wall s. That is indeed the case: m1n1m um: maximum: 0.6389 (adiabatic) < 0.7099 ( numerical) < 0 7262 (isothermal) 1. 6035 (isothermal) < 1. 6420 (numerical) < 1. 9402 (adiabatic) The c orrectness of the numerical solution is also supported by checking the overall energy balance. First define the power inputthe mean rate at which work is done 011 the gas by some external f orce -as

PAGE 74

64 't' P = _!_ f { f pd'ef} dt t O v'(t) (3.53) where 7 is the period of piston oscillations, i.e. T = 21rlw. The mean rate at which heat is transferred out of some domain is 't' Q f I f kVTMS}dt, o swc:tCt ) (3.54) where S we t(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 1 er 1 =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 Q c hann e i = 0.5516D watts, and that through the right chamber walls plus the piston surface is Q r i g ht = 2.7297D watts. The power input is P=0.8415D watts. According to the first law of thermodynamics, the following expression holds, Q e ft Q c bann e l Qrigh t = 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.

PAGE 75

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 (-y = l.4 for air and -y=l.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 Pr, which can be obtained via the equation of state as P r= P r / (RT r ) with R being 287 J / (KgK) for air and 2077.2 J / (KgK) for helium For the meanings of the geometric parameters listed in Table 4.1 please refer to Figs. 1.3 and 1.4. 65

PAGE 76

Tabel 4.1 List of the cases investigated Case Model a b / a c / a d / a e / a L / a 'P o Gas P r T r T c/ T r T b/ T r f Ci M No. [m] [ o ] [attn] [K] [Hz] xl0 3 xl0 6 1 Bl 1 8 8 N I A 7 12 90 Heli11m 1 288 2 0.9 1.1 1 0.229 6.291 2 Al 1 8 8 N I 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 H e li11m 1 288.2 0.9 1.1 1 0.229 6.291 4 B 2 1 8 8 6 7 24 90 H e lium 1 288.2 0.9 1 .1 1 0.229 6.291 5 B 2 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 I A 2 10 45 Heli11m ...... 1 288.2 0.9 1.1 1 0.229 6 291 7 Al 1 6 3 N I A 2 10 75 H e lium 1 288 2 0.9 1 1 1 0.229 6.291 8 Al 1 6 3 N I A 2 10 90 H e li11m 1 288.2 0 9 1.1 1 0.229 6.291 9 Al 1 6 3 N I A 2 10 105 H e li11m 1 288 2 0.9 1 .1 1 0.229 6.291 10 Al 1 6 3 N I A 2 10 135 Helium 1 288.2 0 .9 1.1 1 0.229 6.291 11 Al 1 6 3 N I A 2 10 90 H e li11m ....... 1 288.2 0.95 1 .05 1 0.229 6.291 1 2 Al 1 6 3 N I A 2 10 90 Helium 1 288.2 0.85 1.15 1 0 .229 6.291 13 Al 1 6 3 N I A 2 10 90 H e li11m 1 288 2 0.8 1.2 1 0.229 6.291 14 Al 1 6 3 N I A 2 10 90 H e lium l 288.2 0.7716 1.2284 1 0.229 6.291 15 Al 1 6 3 N I A 2 10 90 h e li11m 1 288.2 0.75 1.25 1 0.229 6.291 16 Al 1 8 8 N I A 7 12 90 Heli11m 1 288.2 0.9 1.1 1 0.229 6.291 1 7 Bl 1 4 8 N I A 7 12 90 Heli11m 1 288.2 0.9 1 1 1 0.229 6.291 N I A d enote s not applicable for the case

PAGE 77

Table 4.1--Continued Case Model a b l a c l a d / a e l a L i a 'P o Gas P r T r T c/ T r T h/ T r f a M No. [m] ( 0 ] [attn] [K] [Hz] xl0 3 xl0 6 18 Bl 1 12 8 N I A 7 12 90 Heli11m ,,_.. 1 288.2 0 .9 1 1 1 0.229 6.291 19 Bl 1 15 8 N I A 7 12 90 Heli11m 1 288.2 0.9 1.1 1 0.229 6.291 20 Al 1 6 3 N I A 2 10 90 H e li1rrn 1 288.2 0.9 1.1 10 0.72 4 62 .9 1

PAGE 78

68 4.1 Velocity and Thermodynamic Fields of a Typical 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 frrst 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 (3lx31 in the left chamber 21xl 1 in the channel and 3lx31 in the right channel) and a time step size of ~t = 7 / 360 ( = 2. 78x10 3 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 = 1r l 4 and wt=1r 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 54 72 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 a11d temperature conditions in the initial quarter cycle discussed at the end of the last chapter Figure 4.1 (b) shows that the (hydrodynamic) 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 2xl0 6 to 2xl0 6 atms with the latter being only about one millionth of the

PAGE 79

69 (a) Thermodynamic pressure a:1 mm, b:8a, c:8a, e: 7a, L=12a, f=1 Hz,

< .._ a. I ...... 0 ..c 0) t: >< a. a:1mm, b:Sa, c:Sa, e:7a, L=12a, f=1Hz, q, 0 : 90 Helium, Pr=1 atm, T,:288.15K, T c =0.9Tr, T h=1.1Tr, a:0. 229 4 X 1 06 .......,,---,,.....,..--,--,--,---,---,..-.,..--,......,,---,,---,.,--,--,--,---,---,..--,--, 2x1 o 6 0 -2x1 o 6 0 1 2 3 4 t I 't cycle Fig. 4.1 Transient processes of pressures. (a) Thermodynamic pressure; (b) Pressure drop between pistons

PAGE 80

70 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. I,: I I <1>:J I + I <1>:J' I l lJ (4.1) where is a general notation representing any specific dependent variable of u, v ph Tor 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 au a v aph and aT 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 loga "' 's have exactly the same slope. Figures 4.3 (a) (b) (c) and (d) give the velocity vectors at time instants wt=O, 1r / 2 1r and 31r / 2 respectively. Note that to l1ave 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

PAGE 81

71 a=1 mm, b=8a, c=8a, e =7a, L=12a, f=1 Hz,


PAGE 82

2 0 10 0 -10 20 -10 ' I I (a) wt=O ,. ---I ~:;.;::::;:~::::::::::::;:::::;:::::::;:::::;:;:;:::;:;:;:;:;:::;:~:= ... : l , I r ,. ::::::::::::::::::~:::::::::::::::::::::: ; ::::::::~::::::::::::~::: , ,. --! ==========-==..,,.. .... ' ' ' I ,. I I ' ' ~i~l~l!l!!!~ll~ l!!!l!l!!~!ll!!!lll~ll!l!!!l!l!~~l!!!!f: f Ifiififilflflttttlf , ,. ---' -:-:-:-:-:-:-:,:-:-:-:-:-:-:-:-:-:-:-:-:-:-:-:-:-:-:-:-:-:-:-:-:-:-:-: 0 x / a , , ' , , ... , ' 10 Fig. 4.3 Velocity fields of Case 1 1 OO(aw)=0 628m /s > .. 20 .................. -----... .. ., --------, _.,__...,. _..,,._., .. ,.,,. 1 0 .......................... ..... .............. ... --0 .... _,.,.,, ---,---..... -1 0 __ .._ .................... ._ .. -, .. ,, __ ______ -20 -....... _._. .... 20 -20 -10 (a) wt O (b) wt (c) wt (b) wt=n/2 i~l~!!llll~ll!!!!!!l!!ll!!lill!!!ll!l!!~!f !ll!!!~l~l!l~l~il 1~ 1 1 1 1 !IIIJI[li{ J .... -:::::::;::::::::~:::::::::~::::~::::;::~:::::::::::::::::::::::::: .......... ............... ., .......... -::: 1 !l!l~!!!llllllllll!ll!lllll!ll!l~ill!l!l~l~l!ll~lilil;; :!lll!!lll!!ll!l~llll!~llllll!lllll!l!lf ~!lllll!ll~!!llll!l~l~l! 0 x/a 7r (d) wt=31r / 2. 10

PAGE 83

20 10 0 -10 -20 -20 .. .. .. ......... . .. .. .. .. .. .. . ................ \, ... .. .......... "''''' .. .. .. .. ................ ,,,, .. ...................... ,..._, \ ---... .. ...... __ .. .. .. __ ......... ,,.,.,.,,,,,,,. ......... ,,,,,,,,,,, .......... ,,,,,,, ...... .... ,,,, ,.,1,,,,,,, .. ..... ,, . .. ............... . . .. .. ... ... .. .. .. .. .. .. .. .. ........... ..... .. . ................. ,,,, .. . ................. .... ..... \ .. ..... ..... -... .... ......... -.......... ----...... --,.,,.,,,,,,. .... ,,,,,,,,. ....... ,,,,,,, .... ,,,,, .. . . .... .. .. .. ................ ........ ,,, .... .. ................. ,,,, .. .. ............. ,, ,,...... \ .. ---. .... ___ ,,,.,,,,,.,,,,,,,. ..... ,,,,,,,,,, .......... ,,,,,,, ..... ,,,,,, ... .. -10 (c) wt=n .. .. .. .. 0 x/a 1 O O ( a cu) =0. 628 m /s . , , , , f ; , .,,;.,. ,,,,,., :t~,,~~ / k :: ,,. ,,. _. _. 0 ~~ ::::=~ \, '' ........... .. ,,,, .... ,,,, ...... ,,,, .... ."'-'' ,. ....... ' ... ...... .. .... . , .,,. , ., , ., , , ,,,,, ,,.,,,,., ; .,, ..,,. .., ~" ,, ,,,, .., .,,. .,,,,._ --~-. ' \ \ -----, ... ... ... ... ,, ... ,, ... ' .. .. ... ..... ... ... ... . , , , , , , , I I ; ,' , ; ,,,.,,..,,. t, ~-:. ~~;-: :: ~-. ' ' ' -, ... .... '' ... ' ' ' ' ' . 10 ,, ., ., _. _. .. ... ... ... ... > 20 10 -10 -20 20 -20 Fig. 4.3 --c ontinued _____ ..,..,_ .......... ---........................ -... ..................... --.......... ,,,,, ". ...................... ,,,.,, ..._,,,,,, \ -~,,,,,, \' \ __ _ _....,.._,. I, x,,, _,,,.,,,, _,,,,,.,,. _.,.,,.,,,,,,,, _____ ,,_,,,,I# I .... ................... ................ ..... ................... .... ...... ,,,,, .... ....... ,,,,,' ....._.._.. .... ,,, , \. .,. .. ......... ............. ,, \ _.._, 'l'lf ,...,.,.,.,,' ~,,,,.,,,,,, _,,,,,,,, _ __ ,,,.,,,_,,II _, ,,,,I# _, ,,,,,, .. --.... ............... ... ... ......... ,,,. _, ............... ,,,. ...... ... ,,,. .... ........ ,,, .. ........,,,,, \' __ ...._...._,,,,\\I .... ,-. .... .......... ~,, \ :,,; '/I, .,_,,,,,, I .,,,,,,, _.,,.,,,,,.,,,,., ,,,, __ ,,,,. ,,,,.,,,,.,,,#, .. . -10 (d) wt=31t / 2 lllllllllllljJllllll~ll~ll~ l ll l !l l llllllllll!l l ll l lll l ~ ~llf 0 x/a , # , I I I I I I I I # ,,,,,,,,,,,,,,,. .,,,,,,,,,,,,,,, .. . ,,,,,.,,,,,,,,,,,, .. ~ .,,,, .,,,,,. ,,_~,,--.I --- "' ---. --_ .. ...... ----, ................ ':'to,,,,,....,..._ ................. .,,,,,,,,, .................. ,, ... . ,,,, ......... ,,,,,, ..... ' ' ' ' ' ' ' . . ,.,., ........ ,, #, ,,,,,,,, .... .,,,,,,,,,, ..... .,,,,,,,,,,,., .. . ,,,,,,,,,,,,,.,,,,.,,,, ... . -;,,;.,,..,,.,_,,, ........ ----.. -------... _:::::: : ,, ,,, ,....,, ................................. .. ,,,,,,,,, .............. ,, ...... .. ,, ,,,,,, ... ,,,, .... .. .,,,,,,,,,, ...... .. ,,,, .. .......... . .. . . f f f , , , . .. ..... .,..,...,. __ ,... . --...... -------._. ............ .. ~ ,, ,,,_ ......... ......... , ,,,,,, .............................. ,,,,,, ... ................... .. ., .. ,,,,,, ...... ..... .. . . . ' .. .. .. . .. .. . .. 10 20

PAGE 84

(a) wt=O (b) wt=n/2 20 20 10 10 0 ctS ........... 0 !l!l!l!~l!l!l~~llll!~l~~!ll!~llllll~lll!llll~l[~!I -10 -10 -20 -20 -10 0 10 20 -20 -10 0 10 x/a x/a Fig. 4.4 Hydrodynamic pressure contours of Case 1. (a) wt =0 (b) wt 1r / 2; (c) wt 7r (d) wt 31r / 2.

PAGE 85

(c) wt=n 20 10 0 -10 -20 -20 -10 0 10 x/a 20 10 -10 -20 20 -20 Fig. 4.4--continued -10 (d) wt=31t/2 0 x/a 10 20

PAGE 86

76 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=O (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=1r (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 cha1nber 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

PAGE 87

(a) wt=O 1 6 20 20 I I 10 10 5 I A 0 0 -10 6 -10 -20 -20 -10 0 10 20 -20 -10 x/a Fig. 4 5 Te mp era tur e c o nto ur s of Case 1 (a) wt O ) (b) wt (b) wt=n / 2 !~l~lm~l~ll~~!ll~ll~lll~llll~~l~llllll~llll!ll!ll ll~illlllllllillll~lr~,11111r~llll~lllllll~llllll ) 0 x/a (c) wt 7r 10 ( d ) wt T /T r A6 1 625 AS 1 6 A4 1 575 A3 1 55 A2 1 525 A1 1 5 AO 1 475 z 1 45 y 1 425 X 1 4 w 1 375 V 1 35 u 1 325 T 1.3 s 1 275 R 1 25 Q 1 225 p 1 2 0 1 175 N 1 15 M 1 125 L 1 1 K 1 075 J 1 05 I 1 025 ....J ....J H 1 G 0 975 F 0.95 E 0 925 0 0.9 C 0 8 7 5 B 0 85 A 0 825 9 0 8 8 0 7 7 5 7 0 7 5 6 0 7 25 5 0.7 4 0 & 7 5 3 0.65 2 0 6 2 5 1 0.6 3 1r /2

PAGE 88

(c) wt=rc 20 I \ 20 ) I ) 10 10 K 0 1) 1 .. ) llll~llllllilililil~llill~ilililllillilillilililillillllllli -10 -10 it~~~1f t@iJ;~J~;~~IIJii~~J~II~ ;:;:;:;:;:;::::: : :::::;:;:;:;::::::~:;::::::::::::::::::::::::::::: \ \ -20 r -20 -20 -10 0 10 20 -20 -10 x/a Fig. 4.5--continued ( d) wt=3rc/2 ) ) I I I l!l!l!!lll!l!l!l~!l!l!l!l!l!lll!lll~~!ijlJ~!lll!~l~~llllll lfj 1 flllll 1 I 0 x/a 10 p 20 T/T r A8 1 825 AS 1 8 A4 1 575 A3 1 .55 A2 1 525 A1 1 5 AO 1 475 z 1 45 y 1 425 X 1.4 w 1 375 V 1 35 u 1 325 T 1 3 s 1 275 R 1 25 Q 1 225 p 1 2 0 1 175 N 1 15 M 1 125 L 1.1 K 1 .0 75 J 1 05 I 1 025 H 1 -.l G 0 9 7 5 00 F 0 95 E 0 925 D 0.9 C 0 .8 75 B 0 85 A 0 825 9 0 .8 8 0. 775 7 0 75 8 0. 725 5 0.7 4 0 875 3 0.85 2 0 825 1 0.8

PAGE 89

(a) wt=O 20 20 10 10 0 rn -0 p -10 H N -10 -20 -20 10 0 10 20 20 -10 x/a Fig. 4.6 Density co ntour s of Case 1. (a) wt O ) (b) wt (b) wt=n/2 !~lll!~!l!l~!~llllll!ll!~llll!ll~!ll~lf lllf lllll!ll!ll!~l~! ~ll~llill!!~ll!!~!!l!~ll!llll!llll~~!l~l!~l!~l!l!!!l!I 0 x/a (c) wt 7r' ) (d) 10 wt p/pr A2 2.4 A1 2.35 AO 2.3 z 2 .25 y 2 2 X 2 15 w 2 1 V 2 05 u 2 T 1 95 s 1 9 R 1 85 Q 1 8 p 1 75 0 1 7 N 1 65 M 1 6 L 1 55 K 1 5 J 1 4 5 I 1 4 H 1 35 G 1 3 F 1 25 E 1 2 -.J l,O 0 1 15 C 1 1 B 1 05 A 1 9 0 95 8 0.9 7 0.85 6 0 8 5 0 75 4 0 7 3 0,65 2 0.6 1 0 55 31r / 2.

PAGE 90

(c) wt=1t 20 2 0 10 10 0 I I I -10 10 20 -20 -20 -10 0 10 20 20 -10 x/a Fig 4. 6--continued (d) wt=31t / 2 111111}11 I 111,11111 I 0 x/a 10 A2 A 1 AO z y X w V u T s R Q p 0 N M L K J H G F E D C 8 A 9 8 7 8 5 4 3 2 1 20 pip 2 4 2 35 2 3 2 25 2 2 2 15 2 1 2 05 2 1 95 1 9 1 85 1 8 1 7 6 1 7 1 85 1 8 1 55 1.5 1 45 1 4 1 35 1 3 1 25 1 2 1 15 00 1 1 0 1 05 1 0 95 0 9 0 85 0.8 0.75 0 7 0 85 0 6 0 55

PAGE 91

... 2 .0 0 0... -1.5 ,__ 0... 1.0 81 E xpa nsion process a=1 mm, b=8a, c=8a, e=7a L = 12a, f = 1 Hz


PAGE 92

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 c ) is defined as the cooling capacity and is denoted by Q c, while the rate at which heat is rejected from the gas through the solid surfaces having hot temperature (Th ) is defined as the heating rate and denoted by (1. 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 eak. The quantities here are: Q c = ll 574D watts (1 = 21.643D watts and Q 1 eak = 3.990D watts. So the total rate of energy into the system is P +QC= 14 075D + ll 574D = 25.649D and the total rate of energy out of the system is (1 + Q 1 eak = 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 c (1 Q 1 eak 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 f or these cases are the same except that in Case 2 the temperature along the entire chann e l 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

PAGE 93

either end segment. Table 4.2 c ompares the performances of the two cases. The negative value of Q ea1c in Case 3 means that heat is added into the gas through the channe] walls having a linear temperature distribution. In the table COP stands for the coefficient of performance which is defined as (Walker 1973) 83 Q C Q e ak p COP 1/ r Table 4.2 Comparison of Models Al and A2 Case 2 Case 3 Model Al Model A2 12.7403 11.8207 17.4084 19.3578 2.4482 0.3929 7.1340 7 .1626 1.7859 1.6503 39 69% 36.67 % COP = Cooling capacity / Power input = Q c / 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) ( 4.3) In Table 4.2 the symbol r,r represents the relative efficiency defined by (Walker 1973 ) 'Y/ r = Actual COP / Carnot COP = COP / (COP)carn ot ( 4.4) The relative effi c iency 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

PAGE 94

84 performance under the given conditions. Next compare Models A2 and B2. Case 3 uses the single-channel model A2 and Case 4 uses the multicl1annel 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 requires a power input of 11. 3809 while Case 3 requires only 7.1626. As a result Case 4 has a COP as low as 0.9947. The use of the multichannel model in Case 4 makes the compression and expansion processes closer to adiabatic than in Case 3 where single channel model is used so Case 4 should have a larger range of pressure variations despite the same Table 4.3 Comparison of Models A2 and B2 Case 3 Case 4 Model A2 Model B2 Q C 11.8207 11.3202 19.3578 22.9956 Q,ea1c -0.3929 -0.3421 p 7.1626 11.3809 COP 1.6503 0.9947 1/ r 36.67% 22.10% 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 excl1anging area of Case 4 is much smaller than that of Case 3. Even though

PAGE 95

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 a11d hydrogen The selection of a gas for a specific Stirling cooling device not only depends on the thermodynamic performance of the gas, but also depends on other factors such as cost, availability safety and so on (Reader and Hooper 1983). With the computer code it is very convenient to compare the performances of working gases. Consider Table 4.4 Comparison of helium and air Case 4 Case 5 Helium Air Q C 11.3202 5.2169 22.9956 14.1773 Q,eak 0.3421 1.2477 p 11.3809 10.2219 COP 0.9947 0.5104 1] r 22.10% 11.34 % 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 wl1ile 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

PAGE 96

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 'P o = 45 75 90 105 and 135 respectively with other parameters being maintained constant Figure 4.8 plots the cooling capacity (Q c ) rate of heat rejected (~ ) power input (P ) and COP versus the phase angle 'P 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 o 110 the cooling capacity a maximum at 'P o 96 and COP a maximum at cp 0 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 'P o 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).

PAGE 97

87 a=1 mm b=6a C=3a e=2a L=1 Oa f=1 Hz Helium, P r =1 atm T r =288.15K, T c =0.9T r T h= 1.1 T r a=0.229 2.0 _______ ,......,.. _________ ....-.._--,--.,.......,. __ --,--.,.......,.---,---,--.,.......,.----.---,--.,.......,--, 2.0 COP 1. 5 Rate of heat rejected Q h 1 5 ..... ca 0... 0... 1 .0 C oo ling capacity, Q c 1 .0 0 0 ..c a 0 a 0.5 Power input, P 0 .5 0.0 ............ ___.__........_..........___._....__....__._ ......... ....__ .................... ...._ ................... ...._ .......... ____ ......__..__.__...__.. ___ 0 0 45 60 75 90 105 120 1 35 Phase angle


PAGE 98

,;-:cu ......., fl .l:! 0 a.. 0 0 made by expansion. 4 3 2 1 0 0.6 88 a=1 mm b=6a C=3a e=2a, L=1 Oa f=1 Hz


PAGE 99

89 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. Tl1at is indicted by the negative cooling capacity in the figure. If the hot -e nd temperature Th is 300 K th e 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 length in the former is 24a and that in the latter is 12a and other parameters in both cases are the same. The loss with the longer channel is 0.6131 which is less than 1 .2 443 --t he counterpart with the shorter channel. However the cooling capacity with the longer channel is Table 4.5 Comparison of cases with different channel lengths Case 2 Case 16 L=24a L=12a Q oss 0.6131 1.2443 Q C 12.7403 13.7450 p 7.1340 8.6155 COP 1. 7859 1.6000 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

PAGE 100

90 dead volume--is also increased therefore the cooling produced by expansion decreases. The net effect (i.e., the available cooling QJ 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 11eeded. Either situation may degrade the cooling performance. For the multicl1annel 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. Q c /( 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 c/ (2bD) versus b / a. The figure shows that both COP and Q c /( 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 c/( 2bD) reaches a maximum value of about 750 Watts / m 2

PAGE 101

0 84 0 82 0... 8 0.80 0.78 0 76 4 91 a=1 mm C=8a e=7a L=12a


PAGE 102

a2 tdif = 1C 92 = -V (4.5) where K is the thermal diffusivity. The residence time during which the gas passes through the channel and exchanges heat with the channel is t = L res U m (4.6) in which Um is the mean magnitude of the gas velocity within the channel. This quantity is estimated by using the ratio of the chamber width (2b) to the channel width (2a) and the mean magnitude of the piston velocity Upm Since the piston velocity is in the form of ewcos(wt+<,0) the quantity U p m is Then one has b b ew U = -U = ma pm a[i It is necessary for es to be greater than tdif, i.e. t >tdif. res Substituting Eqs. (4 .5 )-(4.7) into Eq. (4.8) and then rearranging terms yield p = 1 a 2 Pr < 1 fi a L (4.7) (4.8) (4.9) The dimensionless parameter {3 defined in Eq. (4.9) is a combination of two other parameters related to the geometry (i.e. b / a and e / L), flow oscillation (i.e. a) and heat transfer (i e. Pr ) For effective regeneration {3 should be kept far below unity. Compare Cases 8 and 20. The parameter {3 in the former is 0.0298 and that in the latter is 0.298. Such an increase in {3 is made possible by increasing the oscillation

PAGE 103

93 frequency by nine times. Figures 4.11 (a) and (b) plot the temperature distributions along the x axis at eight time instants for the two cases respectively. For the case with the smaller {3 at all the time instants and within almost the entire channel, the gas temperature is always very close to the local channel wall temperature. This indicates that the channel functions well as a regenerator. However for the case with the larger {3 the te1nperature of the gas within the channel deviates significantly from that of the local wall. As a result the case with {3=0.0298 has a COP of 1.8606 but the case with {3 = 0.298 has merely a COP of 0.9265. The comparison made here clearly shows the role played by the parameter {3 and the need to keep this well below unity in a typical channel regenerator.

PAGE 104

94 (a ) a=1 mm, b=6a C=3a, =2a, L=1 Oa, 'P a =90 f=1 Hz, U=0.229 Helium P r =1 atm T r =288.15K T c =0.9T r T h =1.1 T r ... ._ -0 >< ...__.. ._ 1 .2 ------------------...-..--.----1.0 T c Cold 0.8 Chamber -10 -5 ,~ = 0.02981 Channel 0 x / a Hot Chamber 5 10 T h a=1 mm b=6a c=3a =2a L=1 Oa < ...__.. ._ 1.4 ----------------I~= 0.2981 1.2 Cold Chamber 1 0 0.8 -10 Hot ~-Channel ~ ~ ;--~ ~ -5 0 x / a Chambe 5 10 F i g 4. 11 Te mp e ratur e d is tributions alon g x axi s for t wo di ffe r e n t {3 s ( a ) {3 = 0.0 2 98 ; ( b ) {3 = 0 298.

PAGE 105

CHAPTER 5 CONCLUDING REMARKS A numerical technique was devised to separate the hydrodynamic and thermodynamic pressures for low Mach number variable density flows with heat transfer which occur in reciprocating thermodynamic machines. Separation of these two pressures differ e nt in orders of magnitudes avoids the catastrophic cancellations in numerically evaluating the pressure gradients, and so the technique can increase computational accuracy and efficiency substantially (potential saving on CPU time is about 30% and that on storage and memory is 50 %) It can be easily incorporated to any existing pressur ebased numerical algorithms Two dimensional models of three -c omponent Stirling microrefrigerators were developed in this study. The models consist of all three heat exchangers: heater cooler and regenerator Tl1e heater and cooler are ended with pistons oscillating out of phase and connected by narrow channels with a linear wall temperature distribution simulating the reg e nerator Numerical computations were made for the time periodic variable density flow and heat transfer within domains with moving boundaries by using the SIMPL E algorithm incorporated with a moving grid technique and with the newly devised pressur esplitting method. Such a solution methodology was proven to be feasible--accurate and efficient--for the complex flow and heat transfer problems in the 95

PAGE 106

96 Stirling-type devices. This is the first reported attempt to systematically discuss the numerical aspects pertaining to the simulations of Stirling devices. Results based on the microrefrigerator models agreed very well with the experimental data obtained from actual devices. The model s regeneration and cooling performance are sensitive to the dimensionless parameter that is the product of Womersley number squared, Prandtl number the relative piston strokes and the ratio of chamber width to channel width. This parameter should be kept far below unity for the connecting channels to function well as a regenerator. As the cold chamber wall temperature decreases the cooling capacity decreases linearly but the power input increases linearly and consequently the performance degrades For a set of given operation conditions there is a lowest temperature below which no cooling is possible. This is caused by the increasing cooling loss through the channel as the cold chamber 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. The phase angle that results in the maximum cooling capacity may not necessarily result in the maximum efficiency. The optimum phase angle for either cooling capacity or efficiency is near 90 The ratio of the chamber width to the channel width also affects the cooling

PAGE 107

97 performance. Under given operation conditions there is an optimum ratio for the cooling capacity and efficiency respectively. Future work on this topic should be in the following directions. First turbulence models for oscillatory internal flows should be incorporated into the computations so that simulations can also be made for larger devices and for higher frequencies. Second model the regenerator as porous medium to simulate the conjugate heat transfer mor e accurately Third use a five-component model instead of the three component model Finally the dimensions of the model should be increased from two to three. These extensions will need a very large increase in computing resources of both computer memory and CPU time and it will be necessary to use parallel computations on Cray s upercomputers.

PAGE 108

REFERENCES Ahn K. H. and Ibrahim M.B. "A 2-D Oscillating Flow Analysis in Stirling Engine Heat Exchanger NASA TM 103781 ICOMP-91 04 1991. Andraka C.E. Goods S.H. Bradshaw, R.W. Moreno J.B. Moss T.A. and Jone s S.A. "NAK Pool-Boiler Bench Scale Receiver Durability Test: Test Results and Materials Ananysis Proceedings of the 29th Intersociety Energy Conversion Engineering Conference AIAA 1994 Vol. 4 pp. 1942 1949. Beam R.M. and Warming R.F. "An Implicit Factored Scheme for the Compressible Navier Stokes Equations AIAA Journal Vo l 16 No. 4 1978, pp. 393-402. Berchowitz D.M. "Free Piston Stirling Coolers for Intermediate Lift Temperatures Proceedings of the 27th Intersociety Energy Conversion Engineering Conference SAE 1992 Vol 5 pp. 115-121. Cairelli J.E. Geng S.M. and Skupinski R.C., "Results From Basline Tests of the SPRE I and Comparison with Code Model Predications Proceedings of the 24th Intersociety Energy Conversion Engineering Conference IEEE 1989 Vol 5 pp. 2 249 2256. Cairelli J.E. Geng S.M. Wong W A. Swee, D.M. Doeberling T.J. Madi F.J and Vora D. C., "1993 Update on Results of SPRE Testing at NASA Lewis Proceedings of the 28th Intersociety Energy Conversion Engineering Conference ACS 1993 Vol.2 pp. 633 638. Cairelli J.E. Swee D.M and Skupinski R.C. and Rauch, S. "Update on Results o f SPRE Testing at NASA Proceedings of the 25th Intersociety Energy Conversion Engineering Conference AIChE 1990, Vol 5 pp. 237 244. Cairelli, J.E. Swee D.M. Wo11g W.A ., Doeberling T.J. and Madi F.J. "Update on Results of SPRE Testing at NASA Lewis Proceedings of the 26th Intersociety Energy Conversion Engineering Conference ANS 1991 Vol.5 pp. 217-222. 98

PAGE 109

99 Chan C.K., Tward E. and Burt, W.W. "Overview of Cryocooler Technologies for Space-Based Electronics and Sensors," in Advances in Cryogenic Engineering, Vol. 35 (Editor: R.W. Fast), Plenum Press New York, 1990, pp. 1239-1250. Chase V. "Here s a Chilling Thought: How to Keep Cool Without CFCs R&D Magazine Vol. 37 No. 10, 1995 pp. 39 40. Demirdzic I. and Perie M. "Finite Volume Method for Prediction of Fluid Flow in Arbitrarily Shaped Domains with Moving Boundaries International Journal for Numerical Methods in Fluids, Vol. 10 1990, pp. 771 790. Diver R.B. Andraka C.E., Moreno J.B. Adkins D.R and Moss T.A. "Trends in Dish Stirling Receiver Designs," Proceedings of the 25th Intersociety Energy Conversion Engineering Conference, AIChE, 1990 Vol.5, pp. 303 310. Dochat G. "Free-Piston Space Stirling Technology Program: An Update Proceedings of the 25th Intersociety Energy Conversion Engineering Conference, AIChE 1990 Vol.5 pp. 219-233. Dochat G. and Dhar M. "Free-Piston Stirling Component Test Power Converter Proceedings of the 26th Intersociety Energy Conversion Engineering Conference ANS 1991 Vol.5 pp 239-244. Fox, R. W. and McDonald A. T ., Introduction to Fluid Mechanics Second Edition John Wiley & Sons New York 1978 pp. 307 308. Gedeon D ., "Manifest: A Computer Program for 2 D Flow Modeling in Stirling Machines NASA Contractor Report 182290 1989. Geng S.M. and Tew, R C. "Comparisons of GLIMPS and HFAST Stirling Engine Code Predictions with Experimental Data Proceedings of the 27th Intersociety Energy Conversion Engineering Conference, SAE 1992, Vol.5 pp. 53 58. Golub G.H. and Ortega, J.M. Scientific Computing and Differential Equations: An In troduction to Numerical Methods Academic Press San Diego CA 1992 pp. 6-8. Hansen A.G. Fluid Mechanics Wiley New York, 1967. Hirsch C. Numerical Computation of Internal and External Flows, John Wiley & Sons New York 1988 Vol. 1 pp 474 476. Hornbeck R.W ., Numerical Methods, Prentice-Hall, Englewood Cliffs NJ 1982 pp 66 89.

PAGE 110

100 Huang S.C. "HFAST--A Harmonic Analysis Program for Stirling Cycles," Proceedings of the 27th Intersociety Energy Conversion Engineering Conference, SAE 1992 Vol.5 pp. 47 52. Ibral1im M. Hashim,W. Tew R.C. and Dudenl1oefer, J.E., "Heat Transfer in Oscillating Flows with Sudden Change in Cross Section," Proceedings of the 27th Intersociety Energy Conversion Engineering Conference, SAE, 1992 Vol. 5 pp. 503-508. Ibrahim, M ., Kannapareddy M., Tew R.C. and Dudenhoefer, J.E., "Instantaneous Heat Transfer Coefficient Based upon Two Dimensional Analysis of Stirling Space Engine Components Proceedings of the 26th Intersociety Energy Conversion Engineering Conference ANS 1991, Vol. 6 pp. 149 159. Isshiki N. Kohira Y. Sato N. and Kinoshita, H., "Experimental Studies on Alpha Type Stirling Cooler Proceedings of the 29th Intersociety Energy Conversion Engineering Conference AIAA 1994 Vol. 4 pp. 1813-1816. Isshiki N. Watanabe H. Shishido K. and Watanabe, K. R&D on Small Solar Stirling Engines TNT 2 & 3 NAS 1," Proceedings of the 26th Intersociety Energy Conversion Engineering Conference ANS 1991 Vol.5 pp. 382 387. Jones, D. "Free Piston Stirling Engine Scaling Study Proceedings of the 25th lntersociety Energy Conversion Engineering Conference AIChE, 1990 Vol.5 pp. 245 249. Karki K. C. "A Calculation Procedure for Viscous F low s at All Speeds in Complex Geometries '' Ph D. Thesis Univers ity of Minnesota June 1986 pp. 39-45. Keck T. Schiel W. and Benz, R ., "An Iru1ovative Dish / Stirling System," Proceedings of the 25th Intersociety Energy Conversio11 Engineering Conference, AIChE 1990 Vol.6 pp. 317-322. Kim S.-T. Cho, K. S. and Chung W. S., "Study of Stirling Cycle Refrigerator," Proceedings of the 28th Intersociety Energy Conversion Engineering Conference ACS 1993 Vol.2 pp. 615 619. Kohler J.W.L. "The Gas Refrigerating Machine and Its Position in Cryogenic Technique in Progress in Cryogenics Vol.2 (Editor: K. Mendelssohn) Academic Press New York 1960, pp 41 67. Kormhauser A A. and Smith J .L. Jr. "Heat Transfer with Oscillating Pressure and Oscillating Flow Proceedings of the 24th Intersociety Energy Conversion Engineering Conference IEEE 1989 Vol. 5 pp. 2347-2353.

PAGE 111

101 Kurzweg U.H. and Zhao A.X. "The Stirling Micro-refrigerator: An Innovative Cooling Device for Space Applications," ISRP 92 93 Final Report-Subcontract No. NGT40015 1993. Loktionov Y.V. "Conceptual Design of Free-Piston Stirling Conversion System for Solar Power Units Proceedings of the 26th lntersociety Energy Conversion Engineering Conference ANS 1991 Vol.5 pp. 370-375. McDonald, K. Berchowitz D., Rosenfeld J. and Lindemuth, J. "Stirling Refrigerator for Space Shuttle Experiments," Proceedings of the 29th Intersociety Energy Conversion Engineering Conference, AIAA, 1994, Vol. 4, pp. 1807-1812. Nambudripad N. "Small refrigerators in Cryogenic Engineering (Editor: B.A. Hands) Academic Press London 1986, pp. 431 455. Organ A.J. Thermodynamics and Gas Dynamics of the Stirling Cycle Machine Cambridge University Press, 1992. Otaka T., Saito K. Saito T. Araoka K. and Kagawa N. "Study on a lOOW-Class Stirling Cycle Cooler -Experimental results of Operating Characteristics," Proceedings of the 28th Intersociety Energy Conversion Engineering Conference, ACS, 1993 Vol.2 pp. 621 626. Patankar S. V. Numerical Heat Transfer and Fluid Flow Hemisphere Washington D.C. 1980. Reader G.T. and Hooper C. Stirling Engines E. & F. N. Spon New York 1983 Shaltens, R.K. and Schreiber, J.G. "Comparison of Conceptual Designs for 24 kWe Advanced Stirling Conversion Systems for Dish Electric App l ications Proceedings of the 24th Intersociety Energy Conversion Engineering Conference, IEEE, 1989 Vol. 5 pp. 2305 2315. Sherman, A. "History, Status and Future Applications of Spaceborne Cryogenic Systems in Advances in Cryogenic Engi11eering, Vol. 27 (Editor: R. W. Fast) Plenum Press New York 1982 pp. 1007-1029. Shimko M.A. Wallis P.N. and Martin C.M. "Diaphragm Defined Stirling Refrigerator Proceedings of the 29th Intersociety Energy Conversion Engineering Conference AIAA 1994 Vol 4 pp. 1796 1801. Shyy W., Computational Modeling for Fluid Flow and lnterfacial Transport Elsevier Amsterdam the Netherlands, 1994 pp. 156-163.

PAGE 112

102 Smith J.L. Jr. Lienhard, J.H. V Tziranis A.K. and Ho, Yung "M.I.T. Stirling Cycle Heat Transfer Apparatus," Proceedings of the 27th Intersociety Energy Conversion Engineering Conference, SAE 1992, Vol. 5 pp. 509-516. Tang X. and Cheng P. "Correlations of the Cycle Averaged Nusselt Number in a Periodically Reversing Pipe Flow Int. Comm. Heat Mass Transfer Vol. 20 No. 2 1993 pp. 161 172. Tew R.C. "Overview of Heat Transfer and Fluid Flow Problem Areas Encountered in Stirling Engine Modeling in Fluid Flow and Heat Transfer in Reciprocating Machinery (Editors: T. Morel J.E. Dudenhoefer T. Uzkan and P.J. Singh), ASME, New York 1987, pp. 77-88. Tew R. C. "Status of Several Stirling Loss Characterization Efforts ans Their Significance for Stirling Space Power Development," Proceedings of the 23rd Intersociety Energy Conversion Engineering Conference ASME, 1988 Vol.I pp 113 119. Tew R.C. and Geng, S.M. "Overview of NASA Supported Stirling Thermodynamic Loss Research Proceedings of the 27th Intersociety Energy Conversion Engineering Conference, SAE 1992, Vol.5 pp. 489-494. Thie1ne L.G. and Swee D.M. "Overview of the NASA Lewis Component Technology Program for Stirling Program for Stirling Power Converters," Proceedings of the 27th Intersociety Energy Conversion Engineering Conference SAE, 1992 Vol.5 pp 283-288. Thompson J.F., Warsi Z.U.A. and Mastin C.W. Numerical Grid Generation: Foundations and Applications, North-Holland, 1985 pp. 307-308. Timmerhaus K.D. and Flynn, T.M. Cryogenic Process Engineering Plenum Press New York 1989. Walker G. Stirling-Cycle Machines Oxford University Press London 1973. Walker G. Fauvel O.R. Reader G., Ellison W. Zylstra S. and Scott M.J. "Commercial Applications for Stirling Refrigerators Proceedings of the 23rd Intersociety Energy Conversion Engineering Conference ASME, 1988 Vol. l pp. 15 19. Walker G. Reader G. Fauvel R. and Bingham E.R. "Stirling Near Ambient Tempetrature Refrigerators: Innovative Compact Designs Proceedings of the 27th Intersociety Energy Conversion Engineering Conference SAE, 1992 Vol.5 pp. 93 96

PAGE 113

103 West, C.D. Principles and Applications of Stirling Engines Van Nostrand Reinhold New York 1986 White M.A. "Implantable Energy Source for Artificial Hearts," Proceedings of the 1st International Symposium on Current Problems for Future Development of Artificial Heart and Assis Device Springer-Verlag 1985 pp. 36 48. Wong W.A., Cairelli J.E. Swee D.M. Doeberling T.J. Lakatos T.F. and Madi F J. "NASA Lewis Stirling Testing and Analysis with Reduced Number of Cooler Tubers Proceedings of the 27th Intersociety Energy Conversion Engineering Co11ference SAE 1992 Vol.5 pp. 443 448

PAGE 114

BIOGRAPHICAL SKETCH Guobao Guo was born on January 17 1963 in Nanyang County, Henao Province the People's Republic of China. In 1980, he enrolled in the Aerospace Engineering Department of the Northwestern Polytechnical University Xi'an, Shaanxi Province, P.R China. He received his bachelor s degree in 1984 and master s in 1987 and then he continued his study as a Ph.D. candidate at that university till 1989. In that year he came to the United States and enrolled as a graduate student in the Department of Aerospace Engineering, Mechanics and Engineering Science of the University of Florida. 104

PAGE 115

I certify that I have read this study and that in my opinion it conforms to acceptable standards of scholarly presentation and is fully adequate, in scope and quality as a dissertation for the degree of Doctor of Philosophy. Ulrich H. Kurzweg, Chairperson Professor of Aerospace Engineering Mechanics and Engineering Science I certify that I have read this study and that in my opinion it conforms to acceptable standards of scholarly presentation and is fully adequate in scope and quality as a dissertation for the degree of Doctor of Philosophy. Chen-Chi Hsu Professor of Aerospace Engineering Mechanics and Engineering Science I certify that I have read this study and that in my opinion it conforms to acceptable standards of scholarly presentation and is fully adequate in scope and quality as a dissertation for the degree of Doctor of Philosophy. Wei Shyy Professor of Aerospace Engineering Mechanics and Engineering Science

PAGE 116

I certify that I have read this study and that in my opinion it conforms to acceptable standards of scholarly presentation and is fully adequate in scope and quality as a dissertation for the degree of Doctor of Philosophy. David W. Mikolaitis Associate Professor of Aerospace Engineering Mechanics and Engineering Science I certify that I have read this study and that in my opinion it conforms to acceptable standards of scholarly presentation and is fully adequate in scop e and quality as a dissertation for the degree of Doctor of Philosophy. e mit N. Sigmon Associate Professor of Mathematics This dissertation was submitted to the Graduate Faculty of th e College o f E ngineering and to the Graduate School and was accepted as partial fulfillment of th e requirements for the degree of Doctor of Philosophy. May 1996 C---->Winfred M. Phillips Dean College of Engineering Karen A. Holbrook Dean Graduate School

PAGE 117

j I 1 LO 1780 199 0 ....._..,, q _, ,_


xml version 1.0 encoding UTF-8
REPORT xmlns http:www.fcla.edudlsmddaitss xmlns:xsi http:www.w3.org2001XMLSchema-instance xsi:schemaLocation http:www.fcla.edudlsmddaitssdaitssReport.xsd
INGEST IEID EEPC6TROU_QSP89A INGEST_TIME 2013-01-23T14:33:44Z PACKAGE AA00012919_00001
AGREEMENT_INFO ACCOUNT UF PROJECT UFDC
FILES