Rodger W. Dyson

Glenn Research Center, Cleveland, Ohio

Scott D. Wilson

Sest, Inc., Middleburg Heights, Ohio

Roy C. Tew

Glenn Research Center, Cleveland, Ohio

Rikako Demko

Sest, Inc., Middleburg Heights, Ohio

Fast Whole-Engine Stirling Analysis

NASA/TM—2005-213960

October 2005

AIAA–2005–5558
The NASA STI Program Office . . . in Profile

Since its founding, NASA has been dedicated to

the advancement of aeronautics and space

science. The NASA Scientific and Technical

Information (STI) Program Office plays a key part

in helping NASA maintain this important role.

The NASA STI Program Office is operated by

Langley Research Center, the Lead Center for

NASA’s scientific and technical information. The

NASA STI Program Office provides access to the

NASA STI Database, the largest collection of

aeronautical and space science STI in the world.

The Program Office is also NASA’s institutional

mechanism for disseminating the results of its

research and development activities. These results

are published by NASA in the NASA STI Report

Series, which includes the following report types:

• TECHNICAL PUBLICATION. Reports of

completed research or a major significant

phase of research that present the results of

NASA programs and include extensive data

or theoretical analysis. Includes compilations

of significant scientific and technical data and

information deemed to be of continuing

reference value. NASA’s counterpart of peerreviewed formal professional papers but

has less stringent limitations on manuscript

length and extent of graphic presentations.

• TECHNICAL MEMORANDUM. Scientific

and technical findings that are preliminary or

of specialized interest, e.g., quick release

reports, working papers, and bibliographies

that contain minimal annotation. Does not

contain extensive analysis.

• CONTRACTOR REPORT. Scientific and

technical findings by NASA-sponsored

contractors and grantees.

• CONFERENCE PUBLICATION. Collected

papers from scientific and technical

conferences, symposia, seminars, or other

meetings sponsored or cosponsored by

NASA.

• SPECIAL PUBLICATION. Scientific,

technical, or historical information from

NASA programs, projects, and missions,

often concerned with subjects having

substantial public interest.

• TECHNICAL TRANSLATION. Englishlanguage translations of foreign scientific

and technical material pertinent to NASA’s

mission.

Specialized services that complement the STI

Program Office’s diverse offerings include

creating custom thesauri, building customized

databases, organizing and publishing research

results . . . even providing videos.

For more information about the NASA STI

Program Office, see the following:

• Access the NASA STI Program Home Page

at http://www.sti.nasa.gov

• E-mail your question via the Internet to

help@sti.nasa.gov

• Fax your question to the NASA Access

Help Desk at 301–621–0134

• Telephone the NASA Access Help Desk at

301–621–0390

• Write to:

 NASA Access Help Desk

 NASA Center for AeroSpace Information

 7121 Standard Drive

 Hanover, MD 21076
Rodger W. Dyson

Glenn Research Center, Cleveland, Ohio

Scott D. Wilson

Sest, Inc., Middleburg Heights, Ohio

Roy C. Tew

Glenn Research Center, Cleveland, Ohio

Rikako Demko

Sest, Inc., Middleburg Heights, Ohio

Fast Whole-Engine Stirling Analysis

NASA/TM—2005-213960

October 2005

National Aeronautics and

Space Administration

Glenn Research Center

Prepared for the

Third International Energy Conversion Engineering Conference

sponsored by the American Institute of Aeronautics and Astronautics

San Francisco, California, August 15–18, 2005

AIAA–2005–5558
Acknowledgments

We'd like to thank Kumud Ajmani, Jim Heidman, Vijay Garg (deceased), and Isaac Lopez for their suggestions and

guidance with high performance computing and simulation. The work described in this paper was performed for

the Science Mission Directorate (SMD) and the Radioisotope Power System (RPS) Program, which provided

funding for this project.

Available from

NASA Center for Aerospace Information

7121 Standard Drive

Hanover, MD 21076

National Technical Information Service

5285 Port Royal Road

Springfield, VA 22100

This report is a formal draft or working

paper, intended to solicit comments and

ideas from a technical peer group.

Trade names or manufacturers’ names are used in this report for

identification only. This usage does not constitute an official

endorsement, either expressed or implied, by the National

Aeronautics and Space Administration.

This report contains preliminary

findings, subject to revision as

analysis proceeds.

Available electronically at http://gltrs.grc.nasa.gov
An experimentally validated approach is described for fast axisymmetric Stirling engine

simulations. These simulations include the entire displacer interior and demonstrate it is

possible to model a complete engine cycle in less than an hour. The focus of this effort

was to demonstrate it is possible to produce useful Stirling engine performance results in

a time-frame short enough to impact design decisions. The combination of utilizing the

latest 64-bit Opteron computer processors, fiber-optical Myrinet communications, dynamic

meshing, and across zone partitioning has enabled solution times at least 240 times faster

than previous attempts at simulating the axisymmetric Stirling engine. A comparison of

the multidimensional results, calibrated one-dimensional results, and known experimental

results is shown.

This preliminary comparison demonstrates that axisymmetric simulations can be very

accurate, but more work remains to improve the simulations through such means as modifying the thermal equilibrium regenerator models, adding fluid-structure interactions, including radiation effects, and incorporating mechanodynamics.

Nomenclature

γ Porosity of the Medium

τ eff Stoke’s Tensor

ρf Fluid Density

ρs Solid Density

~v Fluid Velocity Vector

Ef Fluid Energy

Es Solid Energy

h Sensible Enthalpy

J Diffusion Flux of Helium

kef f Effective Thermoconductivity

p Pressure

Sh

f Fluid Enthalpy Source Term

T Temperature

I. Introduction

Fast Whole-Engine Stirling Analysis 

Rodger W. Dyson 

National Aeronautics and Space Administration 

Glenn Research Center 

Cleveland, Ohio 44135 

Scott D. Wilson 

Sest, Inc. 

Middleburg Heights, Ohio 44130 

Roy C. Tew 

National Aeronautics and Space Administration 

Glenn Research Center 

Cleveland, Ohio 44135 

Rikako Demko 

Sest, Inc. 

Middleburg Heights, Ohio 44130 

NASA/TM—2005-213960 1

Great strides have been made in the last decade in expanding the designer’s analysis options 

for building improved free-piston Stirling engines1

 ranging from thermal and structural dynamic 

analysis,2,3 Computer-assisted-design (CAD) software for complete three-dimensional part 
• Structural

– linear or nonlinear static or transient deformation

– modal, spectrum, harmonic, and random vibration dynamics

– linear or nonlinear buckling

• Thermal

– steady-state, transient, conduction, convection, and radiation

• CFD

– Steady-state/transient, compressible/incompressible, laminar/turbulent, forced/natural convection, conjugate and radiation heat transfer, laminar to turbulent transition modeling, scalable

parallel performance

• Electromagnetics

– electrostatics, magnetostatics, low-frequency electromagnetics, current conduction, circuit analysis, harmonic, transient, modal, scattering, perfect electric and magnetic conductors, impedence

boundaries, near and far electromagnetic field extension, frequency selective surface, antenna

radiation patterns, and specific absorption rate.

• Coupled Physics

– acoustics-structural, electric-magnetic, fluid-structural, fluid-thermal, electromagnetic-thermal,

piezoelectric, piezoresitive, thermal-electric, thermal-structural, electromagnetic-thermal-structural,

and eletrostatic-structural

Recently, attempts at utilizing these commercially available options have focused on only the thermal

dynamic cycle simulation because it is difficult to achieve a useful solution in a reasonable time, let alone

attempt a more comprehensive simulation. For example, previous attempts at whole engine simulation

with these multiphysics commercial codes ended in unconverged solutions.31, 41 Although a recent simplified

axisymmetric solution75 shows reasonable energy balance convergence, the results are not validated yet by

experiment.

The first reported successful three-dimensional whole Stirling engine simulation was reported by Mahkamov42, 43, 44

in the United Kingdom, but due to export control regulations, limited information and verfication is available. Apparently a gamma-type Stirling engine was simulated using commercially available software with

good experimental agreement. It is not clear how much detail was simulated nor how practical the results

are to engine design. The second reported successful three-dimensional Stirling engine simulation was by

design,4

 one-dimensional thermodynamic cycle analysis5,6,7,8,9,10,11,12,13,14,15,16,17,18, and optimization, 

and one-dimensional system mechanodynamic simulation.19,20,21,22,23,24,25 A review of the progress 

in Stirling analysis shows a steady relaxation of assumptions/limitations imposed on the analysis 

over the past century.26

Initial two-dimensional computational fluid dynamic (CFD) simulation work began over 

twenty years ago and focused on specific regions of the engine but were not time-efficient and the 

projects were put on hold.6,27,28,29,30 And other researchers attempted to build specific in-house 

research codes with mixed success at modeling specific components.31,32,33,34,35,36

More recently, multidimensional multiphysics simulations that allow for realistic Stirling 

engine system solutions now exist and include the following capabilities that are potentially 

useful for Stirling convertor design:37,38,39,40 

NASA/TM—2005-213960 2
Zhang45 in China in which it took approximately 3 months to reach a steady-harmonic converged solution

on a simplified, academic geometry that did not include the appendix gap, solid walls, or internal displacer.

A multi-dimensional analysis has an important place in Stirling engine design as is discussed at length

in Ref.46 and outlined below:

• verify the one-dimensional results,

• properly simulate inherently three-dimensional turbulence–including transition,

• provide empirical heat transfer and friction factor coefficients for complex geometries before the hardware is built

• integrate all the parts into a CAD system and test for structural and relative motion issues in the

overall design

• assist experiments with determining hard to reach data

• provide a fluid-structure interaction capability

• generate linear reduced order models for controller development

• solve high-power Stirling applications in which the one-dimensional flow assumption breaks down

• identify areas of excessive flow losses due to unintended dead zones, recirculation zones, dissipative

turbulence and other losses such as:47, 48, 49

1. Inefficient heat exchange and pressure loss in the regenerator, heater and cooler,

2. Gas spring and working space loss due to hysterisis and turbulence,

3. Appendix gap losses due to pumping and shuttle effects,

4. Mixing gas losses from nonuniform temperature and flow distributions perpendicular to primary

engine flow axis,

5. Conduction losses from the hot to cold regions

6. Losses due to combined radiation, conduction and convection in hollow spaces

7. And in general, inaccurate loss representations due to use of 1-D flow design codes to account

for flow and heat transfer through area changes (between components) where phenomena such as

flow separation and jetting from tube into regenerator may occur.

The current multiphysics simulation tools are not optimized for Stirling analysis as discussed in Ref.50

For example, recently developed oscillating time advance strategies51, 52, 53, 54, 55, 56, 57 and high order spatial

differencing could make modeling oscillating turbulent transition a reality. But given the cost of developing

new codes, it is desirable to fully utilize current technology first.

This paper will present the first successful U.S. axisymmetric Stirling engine simulation incorporating an

actual engine’s geometry including appendix gap, displacer seals (inner and outer), solid walls, and displacer

interior. Only the flexures, piston seals, bounce space, radiation shields, and heat exchanger fins were not

included for this initial work, but reasonable agreement with experiment has been achieved despite these

omissions.

II. Description of the Problem

The dual opposed configuration shown in figure 1

58, 59, 60, 61 is being developed for multimission use

(i.e., for use in atmospheres and space), including providing electric power for potential missions such as

unmanned Mars rovers and deep space missions.62, 63

NASA/TM—2005-213960 3
Figure 1. Dual Opposed Stirling Convertors Reduce

Vibration (Schreiber)

Only the Stirling engine part of the convertor

(Fig. 2) is simulated multi-dimensionally although

one could anticipate the entire convertor may one

day be simulated using modified multiphysics software.

An example of the desired full Stirling engine

simulation is shown in Fig. 3. It is anticipated that

full 3D simulations will provide a level of geometric

and flow detail necessary for further design improvements. For example, hardware experiments have shown that large performance gains can be made by varying

manifolds and heat exchanger designs to improve flow distributions in the heat exchangers.59, 64, 65, 66

Figure 2. Stirling Engine Part of Convertor Simulated Only

Figure 3. Example Transient Simulation of Stirling

Convertor-Colored by Temperature

The kinds of flow and geometry that occur in a

Stirling engine for which a modeling technique must

address are as follows:49, 48

1. Oscillating flow which changes the effective

conduction, flow friction coefficients and heat

diffusion.67

2. Low mach number flow (no shocks),

3. Compressible flow due to enclosed varying volumes and heat transfer effects,

4. Laminar, Transitional, and Turbulent flow

with Reynolds numbers from 100 to 10000

(based on various length scales pertinent to

the flow region),

5. Conjugate heat transfer, thermal dispersion

and local thermal nonequilibrium in the

porous media regenerator,

6. Micron to Millimeter scale geometry6

NASA/TM—2005-213960 4
7. Sliding Interfaces in the appendix gap

8. Deforming Flow Regions Due to Compression

and Expansion

III. Current Multi-Dimensional Stirling Cycle Analysis Tools

At the level of multidimensional analysis, relatively little has been completed because the third order

analysis, including Sage,5 and to some extent, DeltaE,13 and MS*2 Stirling Cycle Code11 are faster and for

the most part have been adequate engineering design tools. However, to improve efficiency, reliability, and

supplement experimental testing (to understand and reduce losses, identify component temperatures, assist

in sensor calibration, correlation, and placement) it will likely require a better understanding of the actual

flow features and heat transfer throughout the engine. Some of the primary multidimensional analysis tools

are shown below.

1. Modified Computer Aided Simulation of Turbulence (CAST)

The modified CAST code69 is based upon the Semi-Implicit Method for Pressure-Linked Equations (SIMPLE)70 method but is restricted to two dimensions.71, 72, 73 It was modified to include oscillatory boundary

conditions and conjugate heat transfer but does not allow for complicated geometry nor sliding interfaces.

It has been used to model Stirling components but has not been extended to a whole engine simulation tool.

A pressure-splitting technique was added74 to reduce the computational requirements. It was based

on separating the thermodynamic and hydrodynamic pressures so that these widely varying scales could

be solved with less round-off error and better efficiency. And the non-acoustic form of the equations was

later incorporated into the code since the fast acoustic waves and small dimensions of the engine essentially

resulted in ”incompressible” like flow behaviour and removing the acoustical pressure disturbance could

enhance computational speed and accuracy.32

2. CFD-ACE

This commercial code has been used to model a two-dimensional representative Stirling engine31, 75, 33, 45 and

it has been used to model an actual engine.41 Full validation of these simulations has not been completed.

It is also based upon the SIMPLE technique. The regenerator is not currently modeled correctly since

thermal equilibrium is assumed between the gas and solid, but as shown later, this error may not be large

compared to the overall engine performance. The code does support sliding interfaces on parallel computers

in three-dimensions and the company is introducing sliding, overlapping Chimera grids which should enable

better appendix gap modeling. Moreover, the highly compressed volumes occuring at both ends of the

displacer may be more easily modeled with Chimera grids. This finite volume code can utilize both structured

and unstructured grids.

3. Fluent

This commercial code is based upon the SIMPLE and PISO methods. It currently has similar regenerator

modeling limitations as CFD-ACE. It does however have a sliding interface for two-dimensional axisymmetric

flows that can be used for appendix gap modeling on parallel computers and will be discussed in more detail

later. It is also finite volume based and can utilize both structured and unstructured grids.

It has been used with mixed success by industry for Stirling engine simulations76 by several commercial

manufacturers. It was also used by Mahkamov43 in what appears to be the first full three-dimensional

simulation, but details are embargoed. Also, in the related business of cryocooler modeling it has been used

with reasonable success on simplified problems77, 78, 79, 81, 82, 83, 84, 85, 86, 87, 88 And this report has utilized it to

provide a detailed axisymmetric simulation.

NASA/TM—2005-213960 5
4. STAR-CD

The Simulation of Turbulent Flow in Arbitrary Regions (STAR)38 code also uses SIMPLE and PISO methods.

Its companion product (STAR-HPC) is the parallel computer version. It also has sliding interfaces and

deforming mesh capability. It has been used in the related field of internal combustion engine piston modeling,

and some Stirling engines and regenerators have been modeled with it.89, 90

5. CFX

This code also uses SIMPLE and PISO methods on unstructured grids. It also has sliding interfaces implemented, but no stirling engine modeling with this software has been publicly published.37 It currently seems

to have the greatest potential for full multiphysics Stirling convertor simulations because of it’s fluid-structure

interaction, magnetic flux analysis, and structural/thermal stress analysis capabilities. It’s dynamic meshing

capabilities are not ideal for highly compressed volumes as occur in the Stirling engine since only the spring

analogy mesh adaption is currently available.

6. Others

While there are other research codes, they are usually limited to modeling only specific regions of the Stirling

engine such as the regenerator, or the displacer.80, 81, 15, 91, 92

IV. Regenerator Modeling

A very important specific area of modeling difficulty is the regenerator. Since the regenerator (depending

upon one’s definition of effective) has roughly 3 to 40 times more effective heat transfer than the heater,93

any inefficiency of the regenerator represents a significant loss for the entire Stirling engine. Hence, any

numerical losses or inaccuracies in this region will disproportionately influence the entire Stirling simulation.

Also, the geometrical shape of the matrix elements in the regenerator is important and as shown in

Figs. 4(b-f), many shapes are possible. For practical calculations a Representative Elementary Volume

(REV) is often used as shown in Fig. 4(a). An REV requires some closure assumptions that can be partially

derived from experiment and potentially from microscale (pore level) numerical simulations.

Since the size and shape of the regenerator fibers do affect Stirling engine performance, numerically

simulating a regenerator has been pursued at both microscopic (direct) and macroscopic (REV) levels of

fidelity.6, 94, 95, 96, 97, 98, 99, 100 An example simulation utilizing Fluent by Harvey from Georgia Tech is shown

in Fig. 5.

However, our initial approach at whole engine simulation simply adds source terms to the governing

equations and does not take into account the different temperatures of the solid matrix and the fluid. The

source terms serve to enforce a pressure drop according to the Darcy-Forcheimer67 equation. And only a

single energy equation is used, although it does include an effective solid/fluid average conductivity68 and

an effective solid/fluid average energy.

The continuity equation is not modified, but the following source term is added to the momentum

equations (for i = x, y, orz):

Si = −



μ

K

vi + C2

1

2

ρvmagvi



(1)

where permeability K = 4.08e

−10m2

, C2 =

2

√

Cf

K

and inertial coefficient, Cf = 0.178. These coefficients

were determined experimentally101 and used without modification here. And the fluid/solid averaged energy

equation can be written as:

NASA/TM—2005-213960 6
(a) Representative Elementary

Volume

(b) Random (c) Wire Mesh

(d) Perforated Disk (e) Foam Metal (f) Sintered Glass (Failed)

Figure 4. Regenerator Geometry

∂

∂t (γρfEf + (1 − γ)ρsEs) + ∇ · (~v(ρfEf + p)) = ∇ · "

kef f∇T −

 X

i

hiJi

!

+ ( ̄hτ · ~v)

#

+ S

h

f

(2)

V. Dangers of Component Modeling

Modeling the entire engine, as demonstrated herein, is preferable to modeling one component at a time due

to the oscillatory nature of the flow. A component requires known unsteady boundary conditions for an entire

cycle. Since a given boundary may have flow going in both directions, the direction of the characteristics

at these boundaries are difficult to define. Incorrect specifications of boundary conditions result in solution

instability, nonconvergence of solutions, and/or convergence to inaccurate results. It is well known that even

for steady harmonic flows, non-physical reflections can occur at inflow and outflow locations.102 In essence,

by artificially truncating the domain to model a single engine component an artificial boundary condition

is required. However, the exact specification of artificial/non-reflecting boundary conditions in nonlinear

turbulent flow is an important and as yet unresolved area of research.103

For example, a subsonic inflow has four characteristics entering and hence four physical boundary conNASA/TM—2005-213960 7
(a) Schematic (b) Pressure (c) Pathlines

Figure 5. Flow Through Wire Mesh Screen Regenerator

ditions need to be specified, typically total pressure, total temperature, and two flow angles, allowing only

one of the velocity components to be free (which is extrapolated from the interior flow solution). A subsonic outflow has only one characteristic entering and hence only one physical boundary condition should be

specified, for example, static pressure.

Simultaneous inflow/outflow along each engine component’s interface is clearly demonstrated by examining flow velocities near zero (from -1 m/s to 1 ms). Regions of blue and red in Fig. 6 show borders of

mixed inflow/outflow in the axial direction. For example, the piston compression space in Fig. 6(a), the

seal exit in Fig. 6(b), the cooler exit in Fig. 6(c), and the appendix gap entrance in Fig. reffig:oscillating(d)

all show mixed inflow/outflow borders at this moment in time. Moreover, these mixed boundaries moving

during the cycle, further complicating the numerical boundary treatment. If a component boundary crosses

any of these indicated areas, the boundary conditions will not be properly defined. Of course in the case of

commercial software users the code will attempt to stabilize the solution with essentially guessed boundary

conditions, but the solutions will be inherently non-physical.

Repeating the same velocity test, but this time examining velocity close to zero in the y-direction, we

see a new set of boundary lines that will pose difficulties resulting in non-physical component solutions in

Fig. 7.

A. Adjust Components within a Full Engine Simulation

Instead of pursuing the mathematically ill-posed problem of component modeling in complicated oscillating

flow, it is easier, more efficient, and more reliable to determine a whole engine steady-harmonic solution

with the procedures described herein. The engine components can be easily examined and it only takes a

few hours to test the effect of component changes since the overall engine has the correct steady-harmonic

temperature distribution.

At the same time, it is possible to numerically determine the one-dimensional equivalent friction factors

and heat transfer terms in the new components required by one-dimensional analysis tools (e.g. Sage). In

this way, a one-dimensional analysis be made more accurate while not requiring the ill-posed boundary conditions that multidimensional component modeling entails.46 Finally, it is worth noting that one-dimensional

analysis does not have difficulty with component modeling because a boundary is always either an inflow or

an outflow, never both.

NASA/TM—2005-213960 8
(a) Entire Engine (b) Seal Region

(c) Compression Space (d) Heater Manifold and Appendix Gap

Figure 6. Oscillating Axial-Direction

NASA/TM—2005-213960 9
(a) Entire Engine (b) Heater Region

(c) Appendix Gap Region (d) Seal Region

Figure 7. Oscillating Radial-Direction

NASA/TM—2005-213960 10
VI. Approach for Fast Simulations

Achieving practical computer simulations of entire Stirling engines depends upon judicial use of computer

hardware, software engineering, grid generation, and algorithm design. A combination of all four areas is

presented next that enables axisymmetric whole engine simulations to be accurately completed in less than

a week.

A. Computer Design

An important area new to Stirling analysis is the use of parallel computers for faster simulations. NASA

Glenn was one of the first organizations to pioneer computer clusters back in 1992. Many universities have

assembled their own low cost computers,104 and many states have supercomputer consortiums105 to assist

in designing and operating small computer clusters that utilize commercial off-the-shelf (COTS) hardware

and software for low cost systems. These systems are referred to as Linux Beowulf clusters106, 107 and now

represent the majority of the top 500 supercomputer systems in the world.108 Recently, many commercial organizations have begun assembling ”Beowulf”-like clusters including Dell,109 Angstrom,110 LinuxNetworx,111

Appro,112 and Microway113 as shown in Fig. 8. They cost a bit more, but provide comprehensive health

monitoring and controls that significantly reduce the adminstration and integration costs. The latest systems

can house 260 processors in a single 7 foot tall tower.

(a) Angstrom Cluster (b) Appro Cluster (c) Linux Networx Cluster

Figure 8. Appro, Angstrom, and Linux Networx Clusters

Our sixteen dual-Opteron (32) processor Microway cluster is shown in Fig. 9. This was selected because

last year it represented the best value for performance and it integrated well with the commercial software,

Fluent. The Opteron chip, at the time, was the only chip that could run 32-bit and 64-bit applications with

an exceptional price/performance ratio. All the whole Stirling engine results shown here were calculated

using it.

A very important consideration when designing your computer cluster is the communication bandwidth

and latency since commercial fluid simulation software utilizes older numerical approaches that require a large

NASA/TM—2005-213960 11
(a) Front View (b) Dual CPU Board

Figure 9. Microway Cluster

number of frequent communications between the compute nodes. In Fig. 10, two of the best communications

network interface card (NIC) and network switches are shown. The Infiniband system can perform 10 Gb/s

bi-directional communication rates or about 20 times faster than the Myrinet in our cluster. Our simulation

utilized Myrinet communications because last year Infiniband was not a proven technology.

B. Software Decision

We chose Fluent as our computational platform because, at the time, it was the only commercial vendor

supplying an axisymmetric sliding interface for parallel computers. Presently, as mentioned in Sec. III,

the other vendors are now also providing similar capability. It is important to always have at least two

commercial solvers available because of the invariably large number of bugs in their software. Generally, a

bug is more easily found if the same simulation is repeated with two different codes. For example, Fluent

version 6.1.22 properly handles heated moving solids, but the new release (version 6.2.16) did not as identified

by Ibrahim.114 The latest version 6.2.19 has been corrected.

Another reason for utilizing at least two commercial codes is the lack of control the user has in the

solvers utilized. Generally, the commercial codes all utilize slight variants of older, well established lower

order numerical techniques. The slight variants are the only independent way to examine the effect of a

given solver and to look for possible defects. Of course developing one’s own software offers more freedom,

but this can be risky in the short-term.

Our second platform utilizes CFD-ACE and reasonable solutions have been achieved with it. It’s primary

limitation has been the need to work in three-dimensions when solving cases including axisymmetry sliding

interfaces. Another strength of this software is it includes a wide range of multiphysics capabilities.

A third very interesting platform is CFX because of it’s integrated multiphysics capability enhanced by

NASA/TM—2005-213960 12
(a) Myrinet Card (b) Myrinet Switch (c) Infiniband Card (d) Infiniband 144 port

Switch

Figure 10. Communications

the merger between ANSYS and CFX. The only apparent limitation of the software is grid layering is not

yet available. Grid layering is defined by the automatic addition or removal of a row of cells as the mesh

domain shrinks or expands and will be discussed later.

Other packages such as STAR-HPC are very capable as well. However, the commercial cost of these

packages are very expensive and it is not practical to test all of them.

Finally, full three-dimensional simulations are still quite challenging and the parallal scalability of the

commercial software is somewhat limited due to the implicit nature of their steady flow solvers. Better

approaches are currently available50 but have not yet been applied to the Stirling cycle analysis.

C. Grid Design

When using commercial codes, it is primarily the grid quality that determines the speed and accuracy

of the solution. Ideally, one could simply read in a three-dimensional CAD geometry and a grid would

automatically be created. That is in fact the implied goal of computational fluid dynamics code developers.

But for complicated geometry as occurs inside a Stirling engine with all it’s moving parts, the technology

currently available still requires a great deal of human intervention. Taking several weeks to develop a grid

is not unusual and is often required to get a reasonable solution. For example, the rate at which conjugate

heat transfer can be solved is directly related to the smoothness and quality of the grid at the interface. This

is particularly true in Stirling engines when the heat capacity of the solid is vastly diffent than the Helium.

If short-cuts are taken in developing the grid, the price is often vastly longer solution times, if a solution is

produced at all.

All the grids produced in this work were done with Gambit, which is the Fluent,Inc. proprietary grid

generator. There are other grid generators available with their own distinct advantages, but despite the bugs

encountered in Gambit, we found it to be a useful tool.

First, the domain is scaled up by 1000 to compensate for the 32-bit precision available in the grid generator.

The Stirling engine, when entirely modeled, has a fairly wide range of length scales ranging from a

thousandth of an inch in the narrow seal and appendix gap region to an inch or more in the heat

exchanger regions. This scaling has the effect of making very small numbers larger so that smoother

grids can be created. Once the mesh is created satisfactorly, the entire grid domain is scaled back down

by a 1000 before utilizing it.

Second, try to use structured or paver unstructured grids wherever possible to reduce the number of bad

cells and to reduce the number of cells required. Avoid unstructured tetrahedrals because of the

considerably reduced rate of solution convergence.

NASA/TM—2005-213960 13
Third, be sure a sliding interface has essentially matching cell sizes on both sides of the interface to prevent

interpolation errors and avoid instabilities.

Fourth, be sure to place grid points in the boundary layer to properly model friction losses.

Fifth, use sizing functions to smoothly expand the cell size from the boundaries at a rate of 10 percent.

Sixth, use layering wherever possible, but avoid grid adaption and remeshing because this slows down the

solution time considerably. Choose a layering size large enough to permit a reasonable 160 time steps

per cycle, but small enough to allow for smooth transitioning to boundary layers.

Seventh, utilize moving boundary layers on all walls that are moving to enhance conjugate heat transfer.

And place a boundary layer on both sides of the interface (on the solid side and the fluid side).

Eighth, include as much geometry as is practical to avoid having to redo the grid later when simplifying

assumptions are found to matter.

Ninth, clean up the geometry first in a CAD package to avoid creating virtual geometry that is not as

robust.

Tenth, test your grid by forcing the piston and displacer to be stationary and simulate the heat transfer in

the engine. Once a converged solution is observed, continue running the solver and look for grid induced

heat sources. It should be noted that Cleveland State University found producing a steady-state result

first provided a better starting condition for the transient solution. But this last suggestion goes a bit

further in providing a means for improving the grid by ”over-converging” the steady solution.

For the successful case presented here, we used about 700,000 grid points and all grid boundary layers

start at a distance of .001 in from the fluid/solid interface.

In Fig. 11, we show the entire engine, piston region and spider region grid utilized in this simulation.

Notice the piston has a moving boundary layer and grid layering is used to avoid small cells when the piston

fully compresses. Also, notice in the spider region that all helium/wall interfaces have boundary layers on

both sides to enhance conjugate heat transfer. The overall engine can be seen to not include any grid inside

the heater or cooler since the heat exchanger wall temperture is specified as constant. Also notice there is

some extra wasted grid in the spider region near the compression space. This is necessary to maintain grid

quality in the other nearby boundary layer regions. Grid quality is more imporant than grid count for faster

conjugate heat transfer calculations.

In Fig. 12, the hot half of the engine is shown, including the gridding in the displacer and appendix gap.

Notice that despite the complexity, tetrahedral grids were not required. The layering inside the displacer

was done to avoid moving all the grid points inside the displacer for faster run-times. For example, without

grid layering (the addition or subtraction of rows of cells), all the grid cells inside the displacer would need to

move. Instead, only the single layer of cells located at the top of the displacer (right side in picture) moves

at a time in Fig. 12(b). Also, the expansion space uses layering to avoid skewed cells when the displacer is

at top dead center (TDC).

In Fig. 13, the gridding utilized in the seals and compression space are shown. The challenge in the seal

region is to avoid overly skewed cells, overly mismatched cells across the interface, and reasonable cell size

variation near the boundary layer. It is quite difficult to balance all these constraints. Layering was used

in the compression region to avoid crushed cells and the axial direction spacing was chosen to maximize the

time-step.

In Fig. 14 the seal inside the displacer is shown. Notice the use of unstructured quadrilaterals throughout

to avoid the inefficient tetrahedral gridding. Also notice just how small the seal gap is even compared to

the thickness of the displacer wall. This region was very difficult to grid and required utilizing suggestion 10

in Section C. Velocities are high in this region and this grid may not be adequate here, but since this seal

region is not being studied this was left as seen.

NASA/TM—2005-213960 14
(a) Entire Engine (b) Piston Region

(c) Spider Space

Figure 11. Grid Strategy

NASA/TM—2005-213960 15
(a) Right Half (b) Layering Inside Dome

(c) Layering in Expansion Space (d) Appendix Gap

(e) Appendix Gap Close

Figure 12. Hot End Grid Strategy

NASA/TM—2005-213960 16
(a) Top and Bottom Seal, Rod, Compression Space (b) Top Seal and Compression

(c) Top Seal Close-up (d) Bottom Seal and Compression

Figure 13. Seal Grid Strategy

NASA/TM—2005-213960 17
(a) Left Half (b) Inside Seal

(c) Inside Seal Close-up (d) Inside Seal Close-up2

Figure 14. Inner Seal Grid Strategy

NASA/TM—2005-213960 18
D. Solver Optimization

While one’s choices are limited with commercial codes, they do offer a variety of low-order solvers to choose

from, including low Mach number variable density flows.115 A comparison of two commercial codes with the

latest high-order numerical techniques demonstrates some of the differences a solver decision can make.50

One of the limitations currently in Fluent is only first order accuracy is available in time when the mesh

moves. CFD-ACE permits a higher accuracy in time with moving meshes, but does not have a coupled

solver. The results shown in the next section utilize the segregrated implicit solver in Fluent, although it

was later discovered the coupled implicit solver is potentially significantly better. Of course, the grid used for

the segregrated solver may need further improvement if a coupled solver is used. It was found that the best

method in Fluent and the best method in CFD-ACE performed comparably on simpler heat transfer test

problems.50 And with sufficient grid quality and density, the answers are comparable to the latest numerical

methods.

One of the challenges to utilizing actual engine geometry is the difficulty of creating a high quality mesh

as discussed previously. However, certain switches in the solvers may be adjusted to attempt to compensate

for poor grid quality in regions of conjugate heat transfer.116 In what follows are specific recommendations

for improving conjugate heat transfer and they contain Fluent specific switches. The other commercial

packages should be modified with their own specific format.

First, skewed meshes require secondary temperature gradients to be turned off to improve stability and

to reduce undershoot/overshoot of local cell temperatures. In Fluent, this is accomplished by the

command, ”(rpsetvar ‘temperature/secondary-gradient? #f)”

Second, in conjugate heat transfer problems with vastly different thermal conductivities between the fluid

and solid zones, cells near the interface can experience round-off errors. To mitigate this effect, we

want to ensure that the coarse levels of the AMG are visited. It is recommended that the default AMG

cycle for energy be changed from Flexible to either F, V, or W-cycle. The F and V-cycles are more

efficient in parallel simulations. The termination criteria should also be reduced to either 0.01 or 0.001

in order to tighten up convergence.

Third, although a reasonable temperature distribution can be achieved using single precision, in order to

ensure that the total heat transfer is balanced the use of double precision is recommended.

Fourth, if convergence still remains problematic with secondary gradients turned off, it is recommended

that explicit under-relaxation of temperature be used to stabilize the solution. When this is activated,

the energy equation is solved with a URF of 1, but the solution is explicitly relaxed with respect

to the previous iteration. This allows for the equations to be efficiently solvered (URF equal 1) and

instabilities due to poor quality mesh and non-linear material properties to be mitigated (relaxation

factor). The GUI entry of URF for energy will specify the explicit relaxation factor to be used. To

active this, the scheme command is: (rpsetvar ‘temperature/explicit-relax? #t)

NASA/TM—2005-213960 19
VII. Comparison with Analytical and Experimental Data

The approach just described was applied to simulating a Stirling engine built by Infinia Corp.117 (formerly

Stirling Technology Company) being considered for use in space missions. The engine is currently being tested

for long-duration missions and hence a good deal of experimental performance information is available for

comparison. The geometry is export controlled for these engines and so limited information may be shared

with the public. The entire simulation was performed first, with no knowledge of the ”correct” answer either

from experiment or calibrated one-dimensional analysis. The operating conditions used in the simulation

are shown below (Table 1):.

Hot-End Temperature, K 923

Cold-End Temperature,K 353

Ambient Temperature, K 293

Frequency, Hz 80

Mean Pressure, Pa 2.429E+6

Power Piston Amplitude, mm 6.0e-3

Table 1. Operating Conditions

The pressures and heat transfer occuring in various parts of the engine over two cycles are shown in

Fig. 15(a). The pressure in the left (cooler end), middle, and right (heater end) of the regenerator is shown

in Fig. 15(b). The average pressure in the expansion and compression spaces are shown to be about equal

in Fig. 15(c). And, the temperature after 412.376 cycles of simulation are shown in Fig. 16 for an instant in

time (135.36 degree Crank angle with reference to 0 degree displacer located at mean position) in the cold

end of the machine and for the hot end are shown in Fig. 17. Notice the effects of radial heat transfer on

the regenerator.

An initial comparison with Sage suggested the axisymmetric simulation was over-predicting the indicated

power as shown in Table 2. And the overall heat transfer in and out fell within a band predicted by Sage in

which heat loss to the environment was included and then assumed excluded (adiabatic).

Axisymmetric Modified Skupinski Sage 1D Modified Skupinski Sage 1D

Simulation With Ambient Heat Loss No Ambient Heat Loss

PV Power, W 79.65 68.77 70.14

Heat In, W 247.315 260.94 193.8

Heat Out, W 168.313 191.9 123.7

PV Efficiency .322 .264 0.362

Pressure Ratio 1.209 1.187 1.187

Regenerator Delta Press., Pa 10466.3 15,888 15790

Heater Delta Press., Pa 682 194 192

Cooler Delta Press., Pa 896 316 315

Pressure Amplitude., Pa 225269 203100 203200

Mean Pressure, Pa 2429130 2378000 2378000

Table 2. Comparison With Sage 1D Results

NASA/TM—2005-213960 20
(a) Pressure and Heat Transfer over Two Cycles (b) Pressure Along Regenerator over Two Cycles

(c) Comparison of Expansion and Compression Pressure

Figure 15. Pressure and Heat Transfer over Two Cycles

NASA/TM—2005-213960 21
(a) Left Half (b) Compression and Cooler

(c) Cooler and Seal (d) Seal Region

(e) Seal Pretty Close (f) Seal Very Close

Figure 16. Cold End Temperature

NASA/TM—2005-213960 22
(a) Right Half (b) Below Heater

(c) Around Heater (d) Regenerator

Figure 17. Hot End and Regenerator Temperature

NASA/TM—2005-213960 23
But then we gathered the actual performance of two identical engines identified as Technology Demonstrator Convertor (TDC) #13 and #14 from over many thousands of hours of operation since August 2003.

As shown in Table 3, the axisymmetric solution which assumed no ambient heat loss to the environment,

predicted not only the higher indicated power more accurately, but the efficiency was essentially within the

variability of each built machine.

A few words of caution about this encouraging result are in order. First, due to experimental limitations

the piston amplitude is not known for certain. The original designs call for a piston amplitude of 6+/- .08 mm.

And the mean operating pressure is 2.517 Mpa compared to the simulated mean pressure of 2.42913 Mpa.

If the piston amplitude in the experiment is exactly 6.00mm, then this mean pressure difference implies the

simulated power would need to be corrected up to 82.5874W. But since the engines are tuned for maximum

performance and the energy measurements are considered more reliable than the piston amplitudes, it seems

appropriate to conclude this simulation did in fact predict the performance and efficiency within the bounds

of variability of the two engines tested. Since a slightly lower piston amplitude would counter-balance the

corrected higher power just shown, it is important to improve the experimentally derived piston amplitude

for future comparisons.

Th Tc freq net Qin PV Pout net Qout

TDC #13 646.0 92.4 81.4 242.1 78.2 163.9

TDC #14 646.5 94.4 81.4 250.4 79.6 170.8

Simulation 650.0 80.0 80.0 247.3 79.7 168.3

%Err-AVE 100.6% 100.3% 98.2% 100.4% 100.9% 100.6%

%Err-TDC#13 100.5% 100.2% 98.2% 102.1% 101.9% 102.7%

%Err-TDC#14 100.6% 100.3% 98.2% 98.8% 100.0% 98.5%

Table 3. Comparison With Experimental Results

VIII. Conclusion

For the first time in the U.S., a fully converged axisymmetric simulation of an actual Stirling engine was

successfully performed. And moreover, by utilizing the latest advances in computer processors, communication hardware, parallel decomposition strategy, dynamic meshing, solution algorithms, and careful grid

design, it is possible to compute one complete Stirling cycle in less than an hour.

While not conclusive until additional testing is completed, the full axisymmetric simulation predicted

the engine performance within the bounds of both engine test’s data. And predicted the overall efficiency

and indicated power considerably better than our initial efforts with using a leading one-dimensional Stirling

analysis prediction on the first attempt. The one-dimensional analysis can certainly be improved by modifying various friction and heat transfer factors, either from experimental data or from this axisymmetric

solution. But the real value of the axisymmetric solution is no adjustments were required because the engine

is simulated from first principles. Of course the one-dimensional analysis is valuable for other reasons such

as fast optimizations and performance mapping.

These results and observations lead to the conclusion that one, two, and three-dimensional modeling

should all be employed and all three paradigms provide important capabilities that when combined provide

a

potent combination of initial design, empirical coefficient adjustment, optimization, and final prototype

demonstration before the first metal is cut.

NASA/TM—2005-213960 24
References

1Beale, W.T.:”Free Piston Stirling Engines - Some Model Tests and Simulations”, [Preliminary investigation of free-piston

Stirling engines.], SAE International Automtive Engineering Conference, Jan. 1969

2Geng, S.M. : ”A 3-D magnetic analysis of a linear alternator for a Stirling power system.”, NASA/TM-2000-210249,

2000

3Geng, S.M. : ”A 3-D magnetic analysis of a Stirling convertor linear alternator under load”, NASA/TM-2001-211084,

2001

4Chua, C.K. : ”Rapid prototyping: principles and applications”, World Scientific, 2003, 420p.

5Gedeon, D. : ”Sage - Object-oriented software for Stirling machine design”, AIAA Paper 94-4106. 1994

6Gedeon, D. : ”Manifest: A Computer Program for 2-D Flow Modeling in Stirling Machines”, NASA CR-182290, 1989

7Gedeon, D. : ”Sage:Stirling-Cycle Model-Class Reference Guide, Third Ed.”, Gedeon Assocatiates, 1999

8Gedeon, D. : ”Sage:User’s Guide, Third Edition” Gedeon Assocatiates, 1999

9Gedeon, D. : ”GLIMPS Version 4: User’s Manual”, Gedeon Associates, July 1992.

10Bauwens, L., ”Consistency, Stability, Convergence of Stirling Engine Models”, Mitch/Stirling Machines/Systems, Inc.

Berkeley, CA 94710, Dept. Mech. Engineering, Univ. of Calif. at Berkeley, Berkely, CA 94720

11Mitchell, M. P. : ”Modelling Pulse Tube Coolers with the MS*2 Stirling Cycle Code, 10th Int. Cryocooler Conference,

Monterey, 1998.

12http://www.bell.mma.edu/ jmcent/laser.html, ”LASER - Linear Analysis of Stirling Cycle Engines and Refrigerators”.

13Ward, B., Swift, G. ”DeltaE- Design Environment for Low-Amplitude ThermoAcoustic Engines” LA-CC-93-8, Feb.

1996, June, 2001

14”ARCOPTR : Ames Research Center Orifice Pulse Tube Refrigerator”, http://arcoptr.arc.nasa.gov/ARCOPTRoptions.html, 2003

15Kuriyama, T. (Japan) Lewis, M. A. Radebaugh, R. Xiao, J. H. (PHYSICAL AND CHEMICAL PROPERTIES DIVISION - 838) Effect of Regenerator Geometry on Pulse Tube Refrigerator Performance Advances in Cryogenic Engineering ,

July, 1997 , Portland, OR - June 01, 1998

16Schock, A. : ”Stirling Engine Nodal Analysis Program”, J. Energy, Vol. 2, No. 6, 1978

17Swift, G.W., Ward, W.C. : ”Simple Harmonic Analysis of Regenerators”, J. of Thermophysics and Heat Transfer, Vol.

10, No. 4, Oct.-Dec. 1996

18Lu, G.Q., Cheng, P. : ” Numerical and experimental study of a Gifford-McMahon-type pulse tube refrigerator” JOURNAL OF THERMOPHYSICS AND HEAT TRANSFER 17 (4): 457-463 OCT-DEC 2003

19Regan, T.F., Lewandowski, E.J.: ”Application of the GRC stirling convertor system dynamic model : final report”

NASA-CR-213364, 2004

20Regan, T.F.: ”Free-Piston Stirling Convertor Simulation Model Includes Case Motion”, AIAA-2003-5931, Aug. 2003

21: Regan, T.F., Gerber, S.S., Roth, M.E. : ”Development of a Dynamic, End-to-End Free Piston Stirling Convertor

Model”, NASA/TM-2004-212941

22Lewandowski, E.J., Regan, T. : ”Overview of the GRC Stirling Convertor System Dynamic Model”, AIAA-2004-5671,

Aug. 2004

23Ulusoy, N. : ”Dynamic Analysis of Free-Piston Stirling Engines”, Ph.D. Thesis, Case Western Reserve University, 1994.

24Benvenuto, G., Monte, F. de : ”Analysis of Free-Piston Stirling Engine/Linear Alternator Systems Part 2: Results”, J.

of Prop. and Power, Vol. 11, No. 5, Sept.-Oct. 1995

25: Qui, S., Peterson, A.A. : ”Linear Dynamic Modeling and Numerical Simulation of an STC Stirling Convertor”, 1st

Int. Energy Conv. Eng. Conf., Aug. 2003

26Dyson, R.W., Wilson, S.D., Tew, R.C.: ”Review of Computational Stirling Analysis Methods”, Second IECEC, Providence, RI, Aug. 2005

27Gedeon, D. : ”Computational Techniques for the two-dimensional Gasdynamic Equations in Stirling Engine Regenerators

and Associated Manifolds”, 1985.

28Gedeon, D. : ”Developments in two-dimensional regenerator modeling”, 23rd IECEC, 1988, Vol. 1, pp39-44.

29Ibrahim, M.B., ”Development of Multi-Dimensional Stirling Engine Code”, August Progress Report on NASA Grant

NAG 3-955., 1993.

30Hall, C.A., Porshing, T.A., Medley, J, Tew, R.C. : ”Multidimensional computer simulation of Stirling cycle engines”,

25th IECEC, 1990, Vol. 5, pp. 407-411

31Ibrahim, M.B., Mittal, M., Simon, T., and Gedeon, D., ”A 2-D CFD model of a Free Piston Stirling Engine for Space

Applications with Annular Heat Exchangers”, AIAA 2nd Int. Energy Conversion Engineering Conference, Providence, RI, Aug.

16-19, 2004. AIAA-2004-5583.

32Tew, R.C., Ibrahim, M.B. : ”Two-dimensional Compressible Non-Acoustic Modeling of Stirling Machine-Type Components”, J. of Prop. and Power, Vol. 19, No. 5, Sept.-Oct. 2003

33Zhang, Z., ”2D Numerical Model of Flow and Heat Transfer for Stirling Engine Components”, Doctoral Dissertation,

Cleveland State University, May 2003.

NASA/TM—2005-213960 25
34Guo, G. : ”Variable-density Oscillatory Flow and Heat Transfer in 2D models of a Stirling Microrefrigerator”, Ph.D.

Thesis, University of Florida, 1996

35Goldberg, L.F. : ”One- and Two-Dimensional Stirling Machine Simulation Using Experimentally Generated Reversing

Flow Turbulence Models”, NASA CR-185285, Aug. 1990

36Tang, H.Z., Cheng, P., Xu, K. ”Numerical simulations of resonant oscillations in a tube ” NUMERICAL HEAT TRANSFER PART A-APPLICATIONS 40 (1): 37-54 JUL 2001

37Computational Fluids Dynamics Company website at http://www-waterloo.ansys.com/cfx.

38Star CD, ”Methology”,Computational Dynamics Limited, 2002 http://www.fnb.maschinenbau.tudarmstadt.de/de/software/starcd/manuals/method/star method.pdf

39Fluent, Inc., website at http://www.fluent.com/

40ESI, Inc., website at http://www.esi-cfd.com/

41Wilson, S.D., Dyson R.W., Tew, R.C., and Ibrahim, M.B. : ”Multi-D CFD Modeling of a Free-Piston Stirling Convertor

at NASA GRC”, AIAA 2nd Int. Energy Conversion Engineering Conference, Providence, RI, Aug. 16-19, 2004, AIAA-2004-5573

42Makhkamov, K., Ingham, D.B. : ”Theoretical Investigations on the Stirling Engine Working Process”, A00-37715,

AIAA-2000-2815, 2000

43Makhkamov, K., Djumanov, D. : ”Three-dimensional CFD modeling of a Stirling engine”, ST TA1 15, 11th Int. Stirling

Engine Conf., Rome, Italy, Nov. 2003

44Mahkamov, K., Djumanov, D., Orunov, B., Korobkov, A. : ”A γ - Type Stirling Engine: Some Results of Second-Order

Numerical Simulations and Experimentals Tests”, 1st Int. Energy Conversion Eng. Conf., AIAA 2003-5933, Aug. 2003

45Zhang, Z., and Ibrahim, M. ”Development of CFD Model for Stirling Engine and its Components”, AIAA 2nd IECEC,

Providence, RI, August 16-19, 2004. AIAA-2004-5674

46Dyson, R.W., Wilson, S.D., Tew, R.C., Demko, R. : ”On the Need for Multidimensional Analysis”, Third IECEC, San

Francisco, Aug. 2005

47Tew, R.C., Cairelli, J.E., Ibrahim, M., Simon, T.W., Gedeon, D.: ”Overview of NASA Multi-Dimensional Stirling

Convertor Code Development and Validation Effort”, NASA/TM-2002-211997, Dec. 2002.

48Hirata, K., Kagawa, N., Takeuchi, M., Yamashita, I., Isshiki, N., Hamaguchi, K. : ”Test Results of Applicative 100 W

Stirling Engine”, 96204.

49Ibrahim, M.B., Tew, R.C., Zhang, Z., Gedeon, D. Simon, T. : ”CFD Modeling of Free-Piston Stirling Engines”,

NASA/TM-2001-211132, IECEC2001-CT-38, Sept. 2001.

50Dyson, R.W., Wilson, S.D., Tew, R.C., Demko, R. : ”Stirling Analysis Comparison of Commercial vs. High-Order

Methods”, Third IECEC, San Francisco, Aug. 2005

51Zuijlen, A., Bijl, H. : ”Implicit and explicit higher order time integration schemes for structural dynamics and fluidstructure interaction computations”, Comp. and Structures, 83, pp. 93-105, 2005

52Tsitouras, Ch., Simos, T.E. : ”Optimized Runge-Kutta pairs for problems with oscillating solutions”, J. Comp. and

Applied Math., 147, pp. 397-409, 2002.

53Franco, J.M. : ”Runge-Kutta methods adapted to the numerical integration of oscillatory problems”, Applied Num.

Math., 50, pp. 427-443, 2004

54Iserles, A. : ”Think globally, act locally: Solving highly-oscillatory ordinary differential equations”, Applied Num. Math.,

43, pp. 145-160, 2002.

55Tsitouras, C. : ”A High-Order Explicit Runge-Kutta Pair for Initial Value Problems with Oscillating Solutions”, Computers and Mathematics with Applications, 37, pp. 31-36, 1999.

56Franco, J.M. : ”Runge-Kutta-Nystrom methods adapted to the numerical integration of perturbed oscillators”, Comp.

Phys. Comm., 147, pp. 770-787, 2002.

57Calvo, M., Franco, J.M., Montijano, J.I., Randez, L. ”Explicit Runge-Kutta methods for initial value problems with

oscillating solutions”, J. Comp. and Applied Mathematics, 76, pp. 195-212, 1996.

58Thieme, L.G.; Schreiber, J.G. : ”NASA GRC Stirling Technology Development Overview”, NASA/TM-2003-212454,

Aug. 2003

59Thieme, L.G., Schreiber, J.G., Mason, L.S. : ”Stirling Technology Development at NASA GRC”, NASA/TM-2001-

211315/REV1, Jan. 2002

60Schreiber, J. : ”Tutorial on Free-Piston Stirling Power Conversion Technology (How Does it Work?), Space Technology

and Applications International Forum, Feb. 2001

61Schreiber, J.G., Thieme, L.G., Mason, L.S.: ”Status of Free-piston Stirling Power Conversion in the United States Space

Program”, 10th International Stirling Engine Conference, Osnabruk, Germany, Sept. 2001.

62Furlong, R., Shaltens, R. : ”Technology Assessment of DOE’s 55-We Stirling Technology Demonstrator Convertor

(TDC)”, NASA/TM-2000-210509, Oct. 2000

63Abelson, R.D., Balint, T.S., Coste, K., Elliott, J.O.,Randolph, J.E., Schmidt, G.R., Schriener, T., Shirley, J.H., Spilker,

T.R.: ”Expanding Frontiers with Standard Radioisotope Power Systems”, JPL D-28902, PP-266 0332, Jan. 12, 2005

64Tew, R., Ibrahim, M., Simon, T., Mantell, S., Gedeon, D., Qui, S., Wood, G.: ”Overview 2003 of NASA Multi-D

Stirling Convertor Code Development and DOE & NASA Stirling Regenerator R & D Efforts”, NASA TM-2004-212908, Space

Technology Applications Int. Forum 2004.

NASA/TM—2005-213960 26
65Tew, R., Dyson, R.W., Wilson, S.D., Demko, R. ”Overview 2004 of NASA-Stirling Convertor CFD Model Development

and Regenerator R & D Efforts”, STAIF 2005, Albuquerque, NM, Feb. 16-17

66Tew, R.C., Thieme, L.G., Dudenhoefer, J.E. : ”Recent Stirling Engine Loss-Understanding Results”, 25th IECEC,

NASA TM 103122, 1990

67Kaviany, M. : ”Performance of a Heat Exchanger Based on Enhanced Heat Diffusion in Fluids by Oscillation: Analysis”,

J. of Heat Transfer, Vol. 112 (49), Feb. 1990

68Kurzweg, U.H. : ”Enhanced heat conduction in oscillating viscous flow within parallel-plate channels”, J. Fluid Mech.,

Vol. 156, 1985, pp. 291-300

69Peric, M., Scheuerer, G.: ”CAST - A Finite Volume Method for Predicting Two-dimensional Flow and Heat Transfer

Phenomena”, GRS-Technische Notiz, SSR-89-01, Sept. 1989.

70Patankar, S.V., Spalding, D.B. : ”A Calculation Procedure for Heat, Mass, and Momentum Transfer in ThreeDimensional Parabolic Flows”, Int. J. Heat Mass Transfer, Vol. 15, p. 1787, 1972

71Ibrahim, M., Tew, R.C. Dudenhoefer, J.G., ”Two dimensionsal numerical simulation of a Stirling engine heat exchanger”,

24th IECEC, 1989, Vol.6, pp2795-2802.

72Ibrahim, M., Tew, R.C., Dudenhoefer, J.G., ”Further two-dimensional code development for Stirling space engine

components”, 25th IECEC, 1990, Vol. 6, pp. 329-335.

73Ibrahim, M., Kannapareddy, M. : ”Computational Heat Transfer Analysis for Oscillatory Channel Flows”, NASA-CR199238, 1992

74Guo, G., Kurzweg, U.H. : ”2D Numerical Simulation of Heat Transfer in a Stirling Microrefrigerator Model”, AIAA

Paper 96-1807, 31st AIAA Thermophysics Conference, New Orleans, LA, June, 1996

75Ibrahim, M.B., Mittal, M., Jiang, N., Simon, T. : ”Validation of Multi-Dimensional Stirling Engine Codes: Modeling

of the Heater Head”, IECEC, San Francisco, Aug. 2005

76Qui, S. ”Preliminary Computational Fluid Dynamics Modeling of STC Stirling Engines”, IECEC, Providence, Rhode

Island, Aug. 2004

77Flake, B., Razani A., ”Modeling Pulse Tube Cryocoolers With CFD”, 2004

78Flake, B. Razani, A., ”Phase Shift and Compressible Fluid Dynamics in Inertance Tubes”, 2004

79Cha, J. ”CFD Simulation of Multi-Dimensional Effects in Inertance Tube Pulse Tube Cryocoolers”, Master’s thesis,

Georgia Inst. of Tech., May, 2004

80Yuan, Z.S. : ”Oscillatory Flow and Heat Transfer in a Stirling Engine Regenerator”, Case Western Reserve University,

Ph.D. Thesis, 1992.

81Harvey, J.P. ”Oscillatory Compressible Flow and Heat Transfer in Porous Media – Application to Cryocooler Regenerators”, Doctoral Dissertation, Georgia Inst. of Tech., Dec. 2003

82Wang, C., Wu, P.Y., Chen, Z.Q. : ”NUMERICAL MODELING OF AN ORIFICE PULSE TUBE REFRIGERATOR”

CRYOGENICS 32 (9): 785-790 1992

83Wang C ”Numerical analysis of 4K pulse tube coolers .1. Numerical simulation ” CRYOGENICS 37 (4): 207-213 APR

1997

84Ju Y.L., Wang C., Zhou Y. ”Numerical simulation and experimental verification of the oscillating flow in pulse tube

refrigerator”, CRYOGENICS 38 (2): 169-176 FEB 1998

85Alexeev A, Gutfinger C ”Resonance gas oscillations in closed tubes: Numerical study and experiments ”, PHYSICS OF

FLUIDS 15 (11): 3397-3408 NOV 2003

86Zhang X.J., Liu H.X., Wang R.Z., Shi F. ”Numerical simulation of heat transfer in regenerator of solid adsorption

refrigeration system” RENEWABLE ENERGY 26 (4): 599-610 AUG 2002

87Zhu S.W., Matsubara Y. ”A numerical method of regenerator ”, CRYOGENICS 44 (2): 131-140 FEB 2004

88Ataer O.E. ”Numerical analysis of regenerators of free-piston type Stirling engines using Lagrangian formulation ”,

INTERNATIONAL JOURNAL OF REFRIGERATION-REVUE INTERNATIONALE DU FROID 25 (5): 640-652 AUG 2002

89Zhang, Z. Posted on the web at http://www.cfdreview.com/jobs.pl?op=view res&rid=264, No specific details since

proprietary, 2003

90Noh D.S., Hong S.K., Ryou H.S., Lee S.H. ”An experimental and numerical study on thermal performance of a regenerator system with ceramic honeycomb ”, KSME INTERNATIONAL JOURNAL 15 (3): 357-365 MAR 2001

91Bauwens, L., ”Entropy Balance and Performance Characterization of the Narrow Basic Pulse-Tube Refrigerator”, AIAA97-0814, 1997

92Bauwens, L., ”Oscillating Flow of a Heat-Conducting Fluid in a Narrow Tube”, J. Fluid Mechanics, Vol. 324, 1996, pp.

135-161

93Gedeon, D., Tew, R. : Personal Communication, 2004

94Ibrahim, M., Rong, W., Simon, T., Gedeon, D. : ”Simulations of Flow and Heat Transfer Inside Regenerators Made of

Stacked Welded Screens Using Periodic Cell Structures”, 2nd IECEC, Aug. 2004

95Ibrahim, M.B., Rong, W., Simon, T.W., Tew, R., Gedeon, D. ”Microscopic Modeling of Unsteady Convective Heat

Transfer in a Stirling Regenerator Matrix”, 1st IECEC, 2003

NASA/TM—2005-213960 27
96Howse J.W., Hansen G.A., Cagliostro D.J., Muske K.R. ”Solving a thermal regenerator model using implicit NewtonKrylov methods”, NUMERICAL HEAT TRANSFER PART A-APPLICATIONS 38 (1): 23-44 JUL 2000

97Nakayama, A., Kuwahara, F. : ”Numerical Modeling of Convective Heat Transfer in Porous Media Using Microscopic

Structures”, Handbook of Porous Media

98Bastian, P. : ”Higher Order Discontinuous Galerkin Methods for Flow and Transport in Porous Media”, 2003

99Nam, K., Jeong, S. : ”Novel flow analysis of regenerator under oscillating flow with pulsating pressure ”, CRYOGENICS

45 (5): 368-379 MAY 2005

100Tanaka, M., Yamashita, I., Chisaka, F. ”FLOW AND HEAT-TRANSFER CHARACTERISTICS OF THE STIRLING

ENGINE REGENERATOR IN AN OSCILLATING FLOW ” JSME INTERNATIONAL JOURNAL SERIES II-FLUIDS ENGINEERING HEAT TRANSFER POWER COMBUSTION THERMOPHYSICAL PROPERTIES 33 (2): 283-289 MAY 1990

101Wood et. al, Sunpower Oscillating Flow Test Rig

102Saxer A.P., Giles M.B.: ”Quasi-3-Dimensional Nonreflecting Boundary-Conditions For Euler Equations Calculations ”,

JOURNAL OF PROPULSION AND POWER 9 (2): 263-271 MAR-APR 1993

103Hagstrom T, Warburton T: ”A new auxiliary variable formulation of high-order local radiation boundary conditions:

corner compatibility conditions and extensions to first-order systems ”, WAVE MOTION 39 (4): 327-338 APR 2004

104Modi, A., ”Unsteady Separated Flow Simulations Using a Cluster of Workstations”, Master’s Thesis, Pennsylvania State

University, May 1999

105http://www.osc.edu/

106Sterling, Thomas Lawrence. ”Beowulf cluster computing with Linux / [edited by] Thomas Sterling ” QA76.58 .B46 2002

2002 Owned by: GODDARD IVV

107David HM Spector 1st ed. ”Building Linux clusters : [scaling Linux for scientific and enterprise applications]” QA76.76

.O63 S6647 2000 2000

108http://www.top500.org/

109http://www.dell.com/

110http://www.angstrom.com/

111http://www.linuxnetworx.com/

112http://www.appro.com/

113http://www.microway.com/

114Ibrahim, M. ”Personal Communication”, April 2005.

115Schneider T., Botta N., Geratz K.J., Klein R. ”Extension of finite volume compressible flow solvers to multi-dimensional,

variable density zero Mach number flows”, JOURNAL OF COMPUTATIONAL PHYSICS 155 (2): 248-286 NOV 1 1999

116Fluent User’s Group Meeting 2005.

117http://www.stirlingtech.com/

NASA/TM—2005-213960 28
This publication is available from the NASA Center for AeroSpace Information, 301–621–0390.

REPORT DOCUMENTATION PAGE

2. REPORT DATE

19. SECURITY CLASSIFICATION

 OF ABSTRACT

18. SECURITY CLASSIFICATION

 OF THIS PAGE

Public reporting burden for this collection of information is estimated to average 1 hour per response, including the time for reviewing instructions, searching existing data sources,

gathering and maintaining the data needed, and completing and reviewing the collection of information. Send comments regarding this burden estimate or any other aspect of this

collection of information, including suggestions for reducing this burden, to Washington Headquarters Services, Directorate for Information Operations and Reports, 1215 Jefferson

Davis Highway, Suite 1204, Arlington, VA 22202-4302, and to the Office of Management and Budget, Paperwork Reduction Project (0704-0188), Washington, DC 20503.

NSN 7540-01-280-5500 Standard Form 298 (Rev. 2-89)

Prescribed by ANSI Std. Z39-18

298-102

Form Approved

OMB No. 0704-0188

12b. DISTRIBUTION CODE

8. PERFORMING ORGANIZATION

 REPORT NUMBER

5. FUNDING NUMBERS

3. REPORT TYPE AND DATES COVERED

4. TITLE AND SUBTITLE

6. AUTHOR(S)

7. PERFORMING ORGANIZATION NAME(S) AND ADDRESS(ES)

11. SUPPLEMENTARY NOTES

12a. DISTRIBUTION/AVAILABILITY STATEMENT

13. ABSTRACT (Maximum 200 words)

14. SUBJECT TERMS

17. SECURITY CLASSIFICATION

 OF REPORT

16. PRICE CODE

15. NUMBER OF PAGES

20. LIMITATION OF ABSTRACT

Unclassified Unclassified

Technical Memorandum

Unclassified

National Aeronautics and Space Administration

John H. Glenn Research Center at Lewis Field

Cleveland, Ohio 44135–3191

1. AGENCY USE ONLY (Leave blank)

10. SPONSORING/MONITORING

 AGENCY REPORT NUMBER

9. SPONSORING/MONITORING AGENCY NAME(S) AND ADDRESS(ES)

National Aeronautics and Space Administration

Washington, DC 20546–0001

Available electronically at http://gltrs.grc.nasa.gov

October 2005

NASA TM—2005–213960

AIAA–2005–5558

E–15276

WBS–22–972–30–01

34

Fast Whole-Engine Stirling Analysis

Rodger W. Dyson, Scott D. Wilson, Roy C. Tew, and Rikako Demko

Stirling cycle analysis

Unclassified -Unlimited

Subject Category: 20

Prepared for the Third International Energy Conversion Engineering Conference sponsored by the American Institute of

Aeronautics and Astronautics, San Francisco, California, August 15–18, 2005. Rodger W. Dyson and Roy C. Tew, NASA

Glenn Research Center; Scott D. Wilson and Rikako Demko, Sest, Inc., 18000 Jefferson Park Road, Middleburg Heights,

Ohio 44130. Responsible person, Rodger W. Dyson, organization code RPT, 216–433–9083.

An experimentally validated approach is described for fast axisymmetric Stirling engine simulations. These simulations

include the entire displacer interior and demonstrate it is possible to model a complete engine cycle in less than an hour.

The focus of this effort was to demonstrate it is possible to produce useful Stirling engine performance results in a

time-frame short enough to impact design decisions. The combination of utilizing the latest 64-bit Opteron computer

processors, fiber-optical Myrinet communications, dynamic meshing, and across zone partitioning has enabled solution

times at least 240 times faster than previous attempts at simulating the axisymmetric Stirling engine. A comparison of the

multidimensional results, calibrated one-dimensional results, and known experimental results is shown. This preliminary

comparison demonstrates that axisymmetric simulations can be very accurate, but more work remains to improve the

simulations through such means as modifying the thermal equilibrium regenerator models, adding fluid-structure

interactions, including radiation effects, and incorporating mechanodynamics.