Timothy F. Regan and Edward J. Lewandowski

Sest, Inc., Middleburg Heights, Ohio

Development of a Linear Stirling System

Model With Varying Heat Inputs

NASA/CR—2007-215012

December 2007
NASA STI Program . . . 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 plays a key part in helping NASA maintain

this important role.

The NASA STI Program operates under the auspices

of the Agency Chief Information Officer. It collects,

organizes, provides for archiving, and disseminates

NASA’s STI. The NASA STI program provides access

to the NASA Aeronautics and Space Database and its

public interface, the NASA Technical Reports Server,

thus providing one of the largest collections of

aeronautical and space science STI in the world.

Results are published in both non-NASA channels and

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 counterpart of peer-reviewed 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 also include creating custom

thesauri, building customized databases, organizing

and publishing research results.

For more information about the NASA STI

program, 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 STI Help Desk

at 301–621–0134

• Telephone the NASA STI Help Desk at

301–621–0390

• Write to:

 NASA Center for AeroSpace Information (CASI)

 7115 Standard Drive

 Hanover, MD 21076–1320
National Aeronautics and

Space Administration

Glenn Research Center

Cleveland, Ohio 44135

Prepared under Contract NNC07TA38T

Timothy F. Regan and Edward J. Lewandowski

Sest, Inc., Middleburg Heights, Ohio

Development of a Linear Stirling System

Model With Varying Heat Inputs

NASA/CR—2007-215012

December 2007

Prepared for the

Fifth International Energy Conversion Engineering Conference and Exhibit (IECEC)

sponsored by the American Institute of Aeronautics and Astronautics

St. Louis, Missouri, June 25–27, 2007
Acknowledgments

This work was performed for NASA Headquarters, Science Mission Directorate under the Radioisotope Power Systems

program and was supported by the NASA Glenn Research Center.

Available from

NASA Center for Aerospace Information

7115 Standard Drive

Hanover, MD 21076–1320

National Technical Information Service

5285 Port Royal Road

Springfield, VA 22161

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

Trade names and trademarks are used in this report for identification

only. Their usage does not constitute an official endorsement,

either expressed or implied, by the National Aeronautics and

Space Administration.

Level of Review: This material has been technically reviewed by expert reviewer(s).
NASA/CR—2007-215012 1

Development of a Linear Stirling System Model 

With Varying Heat Inputs 

Timothy F. Regan and Edward J. Lewandowski 

Sest, Inc. 

Middleburg Heights, Ohio 44130 

Abstract 

The linear model of the Stirling system developed by NASA Glenn Research Center (GRC) has been extended to 

include a user-specified heat input. Previously developed linear models were limited to the Stirling convertor and 

electrical load. They represented the thermodynamic cycle with pressure factors that remained constant. The 

numerical values of the pressure factors were generated by linearizing GRC’s nonlinear System Dynamic Model 

(SDM) of the convertor at a chosen operating point. The pressure factors were fixed for that operating point, thus, 

the model lost accuracy if a transition to a different operating point were simulated. Although the previous linear 

model was used in developing controllers that manipulated current, voltage, and piston position, it could not be used 

in the development of control algorithms that regulated hot-end temperature. This basic model was extended to 

include the thermal dynamics associated with a hot-end temperature that varies over time in response to external 

changes as well as to changes in the Stirling cycle. The linear model described herein includes not only dynamics of 

the piston, displacer, gas, and electrical circuit, but also the transient effects of the heater head thermal inertia. The 

linear version algebraically couples two separate linear dynamic models, one model of the Stirling convertor and one 

model of the thermal system, through the pressure factors. The thermal system model includes heat flow of heat 

transfer fluid, insulation loss, and temperature drops from the heat source to the Stirling convertor expansion space. 

The linear model was compared to a nonlinear model, and performance was very similar. The resulting linear model 

can be implemented in a variety of computing environments, and is suitable for analysis with classical and state 

space controls analysis techniques. 

Nomenclature 

Ad Displacer area (m2) 

Ap Piston area (m2) 

Arod Displacer rod area (m2

) 

Ch Heater head thermal capacitance (J/K) 

Ct Tuning capacitance (μF) 

Dd Displacer damping (N-s/m) 

Dp Piston damping (N-s/m) 

f Operating frequency (Hz) 

Ialt Alternator current (A) 

Kbounce Bounce space spring rate (N/m) 

Kd Displacer spring rate (N/m) 

Kmagnet Alternator magnet space spring rate (N/m) 

Kp Piston spring rate (N/m) 

Lalt Alternator inductance (H) 

Md Effective displacer mass (kg) 

Mp Effective piston mass (kg) 

MW Gas molecular weight 

N Number of turns on the alternator winding 

P Stirling cycle dynamic pressure (Pa) 

Pc Compression space pressure (Pa) 

Pd Displacer pressure factor (Pa/m) 

Pe Expansion space pressure (Pa) 

Pexp Expansion space PV power (W) 

Pp Piston pressure factor (Pa/m) 

Qe Heat into the Stirling cycle (W) 

Qin Heat delivered by the heat source (W) 

Qout Heat rejected (W) 

Ralt Alternator resistance (Ω) 

Rc Compression space resistance (W/m-K) 

Re Expansion space resistance (W/m-K) 

Rgas Gas constant (J/kg-K) 

Rh Heater head resistance (W/m-K) 

Rhh Heater head conduction loss resistance (W/mK) 

Rins Insulation resistance (W/m-K) 

t time (sec) 

Ta Ambient temperature (K) 

Talt Alternator temperature (K) 

Tc Compression space temperature (K) 

Te Expansion space temperature (K) 

Th Hot heat exchanger temperature (K) 

Tk Cold heat exchanger temperature (K) 

Tr Regenerator space temperature (K) 

Tsource Heat source temperature (K) 

Vc Compression space volume (m3

) 

Vco Mean compression space volume (m3

) 

Ve Expansion space volume (m3

) 

Veo Mean expansion space volume (m3

) 

Vh Hot heat exchanger volume (m3

) 

Vk Cold heat exchanger volume (m3

) 
NASA/CR—2007-215012 2

Vr Regenerator volume (m3

) 

xd Displacer position (m) 

Xd Displacer position amplitude (m) 

xp Piston position (m) 

Xp Piston position amplitude (m) 

ΔP Pressure drop across regenerator, heater and 

cooler (Pa) 

ηmag Alternator magnetic efficiency 

Φ Flux (Wb) 

Φd Displacer phase angle (deg) 

Φp Pressure phase angle (deg) 

ω Operating frequency (rad/sec) 

I. Introduction 

NASA’s Glenn Research Center is using its System Dynamic Model (SDM) to develop linear models of freepiston Stirling convertors for use in real-time simulations of Stirling convertors. The outputs of such a simulation 

could be used in place of actual convertor hardware in tests of controllers and other convertor support hardware. The 

SDM is a non-linear model that simulates Stirling convertor system dynamics for systems of arbitrary complexity. 

SDM has been used to simulate a dual-convertor system with controllers based on pulse-width modulated (pwm) 

power electronics, however the maximum time step allowed for a 50 kHz pwm controller simulation is on the order 

of nano-seconds. Consequently, no more than a few cycles of such a system can be simulated on SDM running on a 

typical desk-top computer. On the other end of the spectrum, when the slow variation of system temperatures 

requires study, the simulation times can become so long as to make the simulation approach impractical. To address 

these extremes, SDM was used to linearize the model at an operating point and the linearized coefficients were used 

to build a linear model. Linear models run faster than SDM and they can run in a variety of computing hardware 

environments. Linear models of Stirling convertor systems, including a dual-opposed system and a single convertor 

system with a balancer were analyzed in detail as part of NASA GRC’s Stirling convertor modeling effort1

. The 

models reported in the literature have included only the mechanical and mounting dynamics. The pressure wave was 

modeled by constants—called pressure factors—that were calculated by SDM to represent the pressure wave at the 

point at which the model was linearized. The models described in Ref. 1 were not equipped to track the relationship 

between output power, hot-end temperature and piston amplitude. Convertor conditions which cause an increase or 

decrease of hot-end temperature such as a change in piston amplitude, a change in the input heat or the rejection 

temperature were not tracked successfully by those linear models. The models lost some accuracy whenever the 

operating point moved away from the operating point at which the model was linearized. That is, the more the piston 

amplitude, hot-end temperature or output power changed during a simulation using the models described in Ref. 1, 

the less accurate the results became. Such a model was not useful in controlling against disturbances in thermal 

environment or user load since disturbances in these areas invalidate the model. 

The linear models developed so far have followed the state-variable formulation in which the system matrix—

the A-matrix—was constrained to have constant coefficients. If the system changes with operating conditions, then 

the coefficients in the model cannot remain constant. In the present work, a method is described for re-computing 

the hot-end temperature based on the calculated expansion space PV power and re-computing certain entries in the 

A-matrix as the system temperatures change. 

The validation data for the models described herein were published by NASA in 1999 from operation of the 

Component Test Power Convertor (CTPC). The CTPC effort published design details as well as test data. The final 

report2,3 has been used for validation of Stirling simulation codes and design codes. A model of the CTPC was 

constructed in SDM and used to produce linearized A-matrix coefficients as in the models described in Ref. 1. 

Section II begins with a description of the CTPC and its SDM model. Section II-A describes the linear electromechanical model. Section II-B describes the thermal system model. Section II-C describes how the models are 

coupled through re-computation of the pressure factors. Section III summarizes the results. 

II. CTPC Dynamic Model 

A high-level schematic of the SDM model of the CTPC is shown in figure 1. Starting at the bottom of the figure, 

the heat input to the Stirling cycle was represented by the heat input Qin, shown with two flow paths emanating from 

it. One flow path represents the heat lost through the insulation to ambient and has the thermal resistance Rins. The 

other flow path represents the heat conducted to the heater. There is a temperature drop from the heat source to the 

heater. The temperature drop was modeled by the resistance Rh. The thermal time constant of the heat input system 

was modeled by the thermal capacitance Ch. This parameter was determined based on the mass of the CTPC heater 

and its heat capacity. 
NASA/CR—2007-215012 3

Not all of the heat entering the heater in the Stirling convertor 

goes into the Stirling cycle. Some is conducted to the cold end 

through paths including the heater head wall, inner cylinder, 

helium gas, regenerator matrix and the displacer. This loss was 

modeled by the heat flow through the thermal resistance Rhh. 

The Stirling cycle thermodynamics were modeled based on 

the Schmidt model4, which assumes an isothermal Stirling cycle, 

but in which pumping losses through heat exchangers and the 

regenerator were considered. The thermodynamic portion of the 

model determines the internal expansion space and compression 

space gas pressures Pe and Pc. These pressures generate the 

driving forces on the displacer mass Md and the piston mass Mp. 

In the model, gas spring and damping forces on the displacer were 

represented by Kd and Dd. Displacer damping is a non-linear 

function of fluid flow through the porous regenerator material and 

through the heat exchanger tubes. The piston is subjected to 

conservative forces from the bounce space gas spring and the 

magnet spring. In the model, these were represented by Kbounce

and Kmagnet. The electrical load circuit also contributes spring 

force if the reactance of the alternator winding is not completely 

balanced by the capacitive reactance of the tuning capacitor. . 

The alternator produces a damping force on the piston based 

on the alternator current. Alternator electrical dynamics were 

determined by the alternator resistance, inductance, motor 

constant, and the tuning capacitor. 

The Stirling convertor was assumed to be rigidly connected 

to the ground. Casing motion, dual-opposed dynamics, and the 

effect of a dynamic balancer could be added to this model if 

desired. 

A. Linear Dynamic Stirling Model 

The ultimate design of a model is determined by a number of factors, including fidelity required, simulation 

time, and model platform. Linearized parameters were used to create a model in Simplorer (Ansoft Corporation) and 

other dynamic modeling platforms. The equations for the electro-mechanical portion of the model were developed in 

Reference 1 and are reprinted as eqs. (1) through (4) below:. 

d

d x

dt

dx =  (1) 

p

p d

p d

p d

d r

d d

d d

d d

d r d x x M

P x A x M

P x A x M

P x A x M

P K A

dt

dx    

 1 1 1 1

⎟

⎟

⎠

⎞

⎜

⎜

⎝

⎛

∂

∂Δ

+

∂

∂ − ⎟

⎟

⎠

⎞

⎜

⎜

⎝

⎛

∂

∂Δ

+ ⎟

⎟

⎠

⎞

⎜

⎜

⎝

⎛

∂

∂ = − + (2) 

p

p x

dt

dx =  (3) 

alt

p mag p

p

p

p

p

p p

d p p

d p

p

p I

dx M

d x N

M

D x x M

P x K A x MP A

dt

dx 1 1 1 1

η

Φ − + ⎟

⎟

⎠

⎞

⎜

⎜

⎝

⎛

∂

∂ − +

∂

∂ = −   (4) 

Equations (1) through (4) represent the dynamic equations for the displacer and the piston. In eq. (2), the force 

due to the pressure wave’s action on the displacer area Ad has been accounted for in two parts. This is because the 

pressure is expressed as a linear combination of the piston and displacer positions. To reconstruct the pressure wave, 

Figure 1.—Schematic of the CTPC 

dynamic model. 
NASA/CR—2007-215012 4

the contribution of both components must be considered. The same discussion applies to the pressure force on the 

piston in eq. (4). In state-variable format, 

X = A⋅ X + B ⋅ u  (5a) 

with the state vector defined as 

[ ]

T X xd xd xp xp =   , (5) 

the A-matrix entries are seen to be the coefficients of the state variables in eqs. (1) through (4), above. In eq. (4), the 

right-most term expresses the piston damping due to alternator current. In a state-variable formulation, Ialt is the 

input, u, and the coefficient of Ialt is the B-matrix entry that contributes to eq. (4). Thus, 

T

dx p mag M p

d B N

⎥

⎥

⎦

⎤

⎢

⎢

⎣

⎡

η

Φ = 1 1 0 0 0 . (6)

B. Thermal System Model 

Figure 2 shows the thermal system modeled as a series of 

thermal resistances with a heat capacitance, heat input, and 

temperature sinks. It is analogous to figure 1 with the 

electromechanical components removed. The schematic 

shows the system as it is configured in the Simplorer 

modeling platform 

Thermal energy from the heater Qin is shown to be 

divided between two flow paths: 1) heat lost through the 

insulation to ambient was determined from Rins and 2) heat 

transferred to the heater head was determined from Rh. The 

resulting heat flow to the heater head causes an increased 

convertor hot-end temperature Th. 

Thermal energy from Th is shown to be divided among 

three flow paths: 1) heat lost via conduction through the 

convertor was determined from Rhh, 2) heat transferred to the 

gas in the expansion space is determined from Re, and 3) heat 

transferred to the heater head material was determined from 

the heater head thermal capacitance Ch. Thermal capacitance 

only becomes apparent during transient analysis; at steady-state operation the heater head material temperature 

equals Th, and heat flow is divided between the conduction loss and heating of the gas expansion space, with the 

bulk of the heat being transferred to the latter. 

Heat is rejected from the convertor to the cold-sink temperature Tk. Most of the rejected heat flows from the 

working fluid. The remainder of the heat rejected is from the conduction losses. The circuit model shown in figure 2 

displays the operation of the thermal system. Circuit solvers such as the one included with Simplorer offer a 

convenient way to calculate the temperatures and heat flows throughout the thermal system. Simplorer will not 

simulate the solution in real time however. To implement a Stirling convertor simulator in real time, the circuit 

solution must be implemented in more direct mathematical language. The method presented here follows a state 

variable formulation. There is one energy storage element in the thermal system: the heater head thermal mass, Ch, 

therefore, there is only one state variable in the formulation. 

1. State Variable Formulation of the Thermal System 

In the circuit diagram of figure 3, as in the circuit of figure 2, there is one energy storage element, the heater 

head thermal mass, Ch. The single state variable is the temperature across it, Th. There are four inputs—the heat in, 

the expansion space PV power, the ambient temperature, and the rejection temperature. The inputs are Qin, Qpv, Ta, 

and Tk, respectively. Heat source Qin models a fixed heat source such as a radioisotope heat source. Qpv is a 

RhhCondLoss

Ch

Rh

H Pexp

Tambient

Rins

Re

Rc Tc

Th

Te

Qsource, 

Tsource

Qout, Tk

Figure 2.—The thermal system was modeled as a 

circuit in SDM. 

Qin, 

Tsource

 Rhh
NASA/CR—2007-215012 5

continuous calculation of the expansion space work over the cycle. Observe that when Qin is greater than Qpv, the 

excess heat flows into the capacitor, Ch. The temperature across Ch will rise as a result. The heat through a thermal 

mass (capacitance) is related to the temperature across it as shown in eq. (7). 

dt

dT Q Q Q C h c = in − pv = h . (7) 

If the heat flow into the capacitor increases, then Qin will be greater than Qpv, and the temperature across the 

capacitor will increase. If the heat flow into the capacitor decreases, then Qin is less than Qpv, and the temperature 

across the capacitor will decrease. 

The circuit of figure 3 is more simple than the circuit of figure 2 in that the thermal resistor, Rh has been 

eliminated. This was done to simplify the math used to calculate the expansion space temperature. The process of 

reducing the thermal circuit of figure 3 to state-variable form was motivated by the observation that the heat 

supplied to the convertor by Qin is the sum of the insulation loss, the heater head conduction, the expansion space 

PV power, and the heat into the thermal capacitance 

dt

dT Q Q Q Q Q Q Q Q C h

in = ins + hh + pv + c = ins + hh + pv + h (8) 

Additionally, the insulation loss and the heater head conduction can be expressed in terms of the inputs and the 

state variable as shown in eqs. (9) and (10). 

ins

h A

ins R

T T Q − = (9) 

hh

h k

hh R

T T Q − = . (10) 

The expansion and compression space temperatures, Te and Tc can be expressed in terms of input, Qpv and the 

state variable, Th, thus: 

Te Th Re Qpv = − ⋅ and Tc Th Rc Qpv = + ⋅ (11) 

By substituting eqs. (9) and (10) into eq. (8) and simplifying, the following expression was obtained: 

h ins

A

h hh

k

h

PV

h

in

h

h h hh ins

C R

T

C R

T

C

Q

C

Q

C

T R R

dt

dT

⋅ − ⋅ + − − − + =

− − ( ) 1 1

 (12) 

Equation (12) is in the basic form of the state equation, 

Qins

Qhh

Qc

Figure 3.—A simplified circuit was used to derive the state-variable equation 

used in the linear model. 

Th Te

Rc

Ch Qpv

Ta

Tk

Rhh Qin

Rins Re

Tc
NASA/CR—2007-215012 6

X = AX + Bu  (13) 

where A is the single coefficient 

h

hh ins

C

R R A a ( ) 1 1

11

− − − + = = (14) 

and B is the vector [b1 b2 b3 b4]

T

 where 

Ch

b 1

1 = (15a) 

Ch

b 1

2

− = (15b) 

Ch Rins

b ⋅ = 1

3 (15c) 

Ch Rhh

b ⋅ = 1

4 (15d) 

Thus, there is a B-vector coefficient for each of the four inputs, 

Qin, QPV, TA and Tk. When the product Bu is evaluated, one number is 

produced. It has dimensions of Kelvins per second. An 

implementation of the state equation in block diagram form is shown 

in figure 4. 

The thermal model determines the hot-end temperature (Th) for 

the Stirling cycle, but does not calculate cycle performance or 

convert heat energy into mechanical energy. This is accomplished by 

coupling the linear Stirling model to the thermal system model. The 

next section will discuss how the model coupling was performed. 

C. Coupling Stirling Linear Model and Thermal System 

Model 

To complete the linear model of the CTPC, the Stirling dynamic 

model has to be coupled to the thermal system model. One approach would be to expand the Stirling cycle model of 

eqs. (1) to (4) to include the thermal system. The concern with this approach is that the response time of the Stirling 

cycle is on the order of milliseconds, while the response time of the thermal system are on the order of seconds or 

minutes. By keeping the electro-mechanical system separate from the thermal system, the likelihood of formulating 

an ill-conditioned matrix is reduced. 

An alternative approach is to algebraically couple two separate linear dynamic models, one model of the 

Stirling cycle and one model of the thermal system, through the pressure factors. The pressure factors are the partial 

derivatives of the pressure wave P with respect to displacer position xd and piston position xp. The equations to 

couple the systems can be derived from the Stirling cycle equations as derived by Urieli and Berchowitz4

. 

−1

⎟

⎟

⎠

⎞

⎜

⎜

⎝

⎛ = ⋅ ⋅ + + + + c

c

e

e

k

k

r

r

h

h w gas T

V

T

V

T

V

T

V

T

V P M R (16) 

Figure 4.—Thermal model block diagram 

used to solve the state equation.

I

Integrator

CONST

Input

KI := A

b1*Qin+b2*Qpv+b3*Ta+b4*Tk

ThDot

Th
NASA/CR—2007-215012 7

In eq. (16) the volumes Ve and Vc are functions of xd and xp. The other variables are assumed to be fixed over a given 

cycle. 

Ve =Veo − Ad xd (17) 

Vc Vco Ap x p Ad Arod xd = − + ( − ) (18) 

The temperatures Tc and Te can be expressed in terms of the heat into the heater head Qh based on the thermal circuit 

shown in figure 2, where the thermal resistance Re is used to calculate the temperature drop from Th to Te, and Rc is 

used to calculate the temperature drop from Tc to Tk. 

Te = Th − QhRe (19) 

Tc = QhRc + Tk (20) 

The regenerator temperature can be expressed as the log mean temperature between Th and Tk: 

ln( / ) h k

h k r T T

T T T − = (21) 

Substituting eqs. (17) to (20) into eq. (16) and differentiating yields the following expressions for the pressure 

factors 

2 ( ) ln( / )

( )

−

⎥

⎦

⎤ ⎢

⎣

⎡

−

− + + − + +

+

− + − ⋅ + = ∂

∂

h e e

eo d d

h

h

h k

r h k

k

k

e c k

co p p d rod d

e c k

p w gas

p T Q R

V A x

T

V

T T

V T T

T

V

Q R T

V A x A A x

Q R T

A M R

x

P (22) 

2 ( ) ln( / )

( )

( ) −

⎥

⎥

⎦

⎤

−

− + + − ⎢

⎣

⎡

+ +

+

− + − ⋅ ⎥

⎦

⎤ ⎢

⎣

⎡

− − +

− = − ∂∂

h e e

eo d d

h

h

h k

r h k

k

k

e c k

co p p d rod d

h e e

d

e c k

d rod w gas

d T Q R

V A x

T

V

T T

V T T

T

V

Q R T

V A x A A x

T Q R

A

Q R T

A A M R xP (23) 

The pressure factors are now functions of Th, Tk, Qe, which can be determined based on the dynamics of the thermal 

system. The pressure factors are also functions of xp and xd. In theory, the pressure factors should be the same 

regardless of the point in the cycle at which xp and xd were measured, but analysis of simulation data has shown that 

the calculated values of the pressure factors do vary over the cycle. It has been observed that the calculated values of 

the pressure factors are closest to the actual value as determined from analysis of the pressure wave when xp = xd and 

are both negative. Further analysis is needed to understand the reason for this and to understand the limitations. 

Given displacer and piston amplitudes Xd and Xp and a displacer phase angle φd, xp and xd are equal and negative 

when 

⎟

⎟

⎠

⎞

⎜

⎜

⎝

⎛

⎥

⎥

⎦

⎤

⎢

⎢

⎣

⎡

− φ

φ = − −

p d d

d d

p p X X

X x X

cos

sin sin tan 1 (24) 

⎟

⎟

⎠

⎞

⎜

⎜

⎝

⎛

+ φ ⎥

⎥

⎦

⎤

⎢

⎢

⎣

⎡

− φ

φ = − − d

p d d

d d

d d X X

X x X

cos

sin sin tan 1 (25) 

These values for xp and xd should be used in eqs. (22) and (23). 

To algebraically couple the thermal system with the Stirling model, it was also necessary to calculate the 

expansion space PV power. The expansion space PV power is approximately equal to the heat flow into the Stirling 

cycle 
NASA/CR—2007-215012 8

∫ Qh ≈ f PdV (26) 

For the expansion space, 

Ve =Veo − Ad xd (27) 

Since 

( ) sin( ) d d d x t = X ωt + φ (28) 

x (t) X sin( t) p = p ω (29) 

then 

cos( ) d

d d d e A X t

dt

V = − ω ω + φ (30) 

Also, the pressure in the expansion space is a function of the pressure factors 

P(t) P x (t) P x (t) = d d + p p (31) 

Substituting eqs. (30) and (31) into eq. (26) and integrating gives the PV power: 

pv d Xd X pPp d Q (32) = π A f sin φ

Determining the PV power is straightforward by use of eq. (32) if the phase angle and amplitudes are known. 

SDM makes use of Simplorer’s state diagrams to calculate them from the piston and displacer position signals. The 

same approach can be taken with Matlab/Simulink (The MathWorks, Inc.) using the Stateflow (The MathWorks, 

Inc.) add-on package. Implementations that make use of the Very-high-speed integrated circuit Hardware Definition 

Language (VHDL) are able to use the state machine definition features of VHDL to determine amplitudes and 

angles. The amplitude of a sinusoidal signal can be determined by observing when the derivative of the signal 

crosses zero. Since the derivative of a sinusoid is a sinusoid lagging by 90°, the zero-crossing of the derivative 

coincides with the maximum or minimum of the original signal. Detection of zero-crossings can be detected by state 

machines programmed to change states when the control signal changes signs. 

An even more direct solution was used in the Simplorer implementation of the linear model that did not require a 

state machine. The PV power was calculated by the product of frequency and the output of a power measurement 

block. The two inputs to the power measurement block were expansion space volume and the pressure wave instead 

of voltage and current. 

If temperature changes in the hot end, then the working fluid redistributes itself within the working space and the 

non-working spaces of the convertor. The non-working spaces comprise the bounce space and the interior of the 

displacer. The mass of fluid in the working space is computed from the total charge, Mtot, as follows: 

⎟

⎟

⎠

⎞

⎜

⎜

⎝

⎛

+ + + + +

⎟

⎟

⎠

⎞

⎜

⎜

⎝

⎛

+ + + +

= ⋅

bounce

bounce

c

co

e

eo

k

k

r

r

h

h

c

co

e

eo

k

k

r

r

h

h

w tot

T

V

T

V

T

V

T

V

T

V

T

V

T

V

T

V

T

V

T

V

T

V

M M (33) 

The pressure factors as calculated by eqs. (22) and (23) were used to re-compute the entries for a21, a23, and a41

in eqs. (34) through (36). 
NASA/CR—2007-215012 9

d d

d rod

x M

P

a K A 1

21 ⋅ ⎟

⎟

⎠

⎞

⎜

⎜

⎝

⎛

∂

∂ = − + ⋅ (34) 

p d

rod x M

P a A 1

23 ⋅

⎟

⎟

⎠

⎞

⎜

⎜

⎝

⎛

∂

∂ = − ⋅ (35) 

d p

p x M

P

a A 1

41 ⋅ ⎟

⎟

⎠

⎞

⎜

⎜

⎝

⎛

∂

∂ = − ⋅ (36) 

p p

p p x M

P

a K A 1

43 ⋅ ⎟

⎟

⎠

⎞

⎜

⎜

⎝

⎛

∂

∂ = − + ⋅ (37) 

The only other coefficients to be calculated are 

those responsible for the displacer damping. Their 

values depend on the pressure drop between the 

expansion and compression spaces. The pressure drop 

was modeled as a sinusoidal pressure expressed as a 

linear combination of the piston velocity and the 

displacer velocity. Values for these pressure drop 

factors were determined by performing an SDM model 

simulation through a temperature range and fitting the 

model’s recorded output to a linear curve with hot-end 

temperature as a parameter. 

1 +β1 = α ⋅ ∂

∂Δ

h

d

T

x

P

 (38) 

2 + β2 = α ⋅ ∂

∂Δ

h

p

T

x

P

 (39) 

Finally, a22 and a24 were calculated: 

d

d

d M

x

P

a A 1

22 ⋅

∂

∂Δ = ⋅ ⋅ (40) 

d

p

d M

x

P

a A 1

24 ⋅

∂

∂Δ = ⋅ ⋅ (41) 

D. Implementation 

The linear models and their coupling equations were implemented in Simplorer and in Matlab/Simulink. These 

two high-level modeling platforms both offer convenient means to compute the integral of the state equations and to 

evaluate the output power and the expansion space PV power. In Simplorer, the integration of a state equation is 

accomplished with the block diagram of figure 5. Simulink offers a ready-made block for state equations which 

requires only the specification of a 4 x 4 system matrix (A-matrix) and a 4 x 2 input matrix (B-matrix) to implement 

the same system shown in figure 5. Implementation of the linear model also required the equations for re-computing 

the six A-matrix entries, and evaluating the expansion space PV power. 

Figure 5.—Simplorer implementation of 4-state, 2-input 

state equation uses equation block and four integrators.

I

I

I

I

x1

x2

x3

x4

State Equations x1Dot:=a11*x1+a12*x2+a13*x3+a14*x4+b11*u1+b12*u2

x2Dot:=a21*x1+a22*x2+a23*x3+a24*x4+b21*u1+b22*u2

x3Dot:=a31*x1+a32*x2+a33*x3+a34*x4+b31*u1+b32*u2

x4Dot:=a41*x1+a42*x2+a43*x3+a44*x4+b41*u1+b42*u2
NASA/CR—2007-215012 10

There are a number of software tools available for translating Matlab/Simulink models into code that can be 

executed on a microprocessor or a field-programmable gate array6

 (FPGA). The advantage of running the model on 

such hardware rather than on a desktop computer is that the model can simulate convertor performance in real time. 

With real-time simulation, the actual thermal time constants can be used in the model and more realistic response 

obtained. Simplorer offers a feature whereby the model shown in figure 4 can be translated into VHDL code. The 

feature can be used to generate a first cut of the VHDL code used to program the FPGA to execute the linear model. 

The code generated is not directly synthesizable—the user must first substitute the non-synthesizable variable 

declarations with appropriate fixed-point declarations. 

By implementing the linear model on a readily-available demo board for a microprocessor or FPGA, and connecting 

the EMF output to a power amplifier, a convertor simulator such as that shown in figure 6 can be constructed. 

The State variable representing the piston velocity is multiplied by the alternator constant and the result is used 

as the main output of the model. That output is equal to the alternator EMF. It is amplified in the power amplifier 

and the output is wired in series with an inductance and a resistance. The values used are equal to the inductance and 

resistance of the alternator winding. The power measured by the watt meter is the actual power output of the 

modeled convertor. In the case of the CTPC, it is 12.5 kW, so the amplifier in the notional schematic of figure 6 is 

not one that is readily available. The voltage on the load end of Ralt is the same as the voltage at the terminals of the 

convertor. Controller hardware and software can be connected and exercised. Various parameters can be changed in 

such a system, perhaps even while the system is running. The changes can be introduced during a controller test to 

verify that the controller performs as expected when the convertors endure radiator temperature changes as well as 

the expected radioisotope decay. The real-time simulator based on the linear model could also find use as a system 

for demonstrating design margins in convertor matching for dual-opposed pairs. For example, two simulators and 

two controllers could be used to compare system performance with spring rates matched to within 0.1, 1, and 10%. 

III. Summary 

The linear model described extends the four-state electro-mechanical linear model previously developed from 

outputs of SDM. The extended model includes a user-specified heat input and a hot-end temperature that varies with 

operating conditions. 

Seven mathematical constructs were required in order to implement the model 

1. A four-state, single-input, single-output state space model of the piston and displacer motion and 

damping as shown in eqs. (1) through (4). 

2. A single-state, four-input, single-output state space model of the thermal processes in the hot end as 

shown in eqs. (12) and (13). 

3. Recalculation of the mass of the working fluid in the working spaces as shown in eq. (33) whenever a 

change in temperature is detected. 

4. Recalculation of the pressure factor of eq. (20) and associated A-matrix entries, a23 (eq. (33)) and a43

(eq. (35)) whenever a change in either temperature, amplitude or phase angle is detected. 

5. Recalculation of the pressure factor of eq. (21) and associated A-matrix entries, a21 (eq. (32)) and a41

(eq. (34)) whenever a change in either temperature, amplitude or phase angle is detected. 

6. Recalculation of the pressure drop factor for displacer damping with respect to displacer velocity from 

eq. (36) and the associated A-matrix entry a22 from eq. (38). 

7. Recalculation of the pressure drop factor for displacer damping with respect to piston velocity from 

eq. (37) and the associated A-matrix entry a24 from eq. (39). 

Lalt Ralt

FPGA

ADC

ADC

DAC

V

+

-

EMF

A CONTROLLER

HeatInSet

Figure 6.—Notional convertor simulator is based on the linear model running in the FPGA and 

producing an EMF signal in real time.

Qin Set 
NASA/CR—2007-215012 11

A possible use of the model is in Stirling convertor simulators. The simulator will generate realistic convertor 

electrical terminal voltages and currents so that convertor systems including controllers can be exercised under a 

variety of conditions. 

References 

1. T.F. Regan and E.J. Lewandowski, “Stirling System Modeling for Linear Dynamics Analysis,” Proceedings of the Third 

International Energy Conversion Engineering Conference (IECEC 2005), San Francisco, CA, 2005. 

2. M. Dhar, “Stirling Space Engine Program, Volume 1—Final Report,” NASA/CR—1999-209164/VOL1. 

3. M. Dhar, “Stirling Space Engine Program, Volume 2—Appendixes A, B, C, and D,” NASA/CR—1999-209164/VOL2. 

4. I. Urieli, and D.M. Berchowitz, Stirling Cycle Engine Analysis, p. 21, Adam Hilger Ltd., Bristol, 1984. 

5. S. Qiu and A. Peterson , “Linear Dynamic Modeling and Numerical Simulation of an STC Stirling Convertor,” Proceedings 

of the 1st International Energy Conversion Engineering Conference (IECEC), Portsmouth, VA, 2003. 

6. D. Pellerin, S. Thibault, Practical FPGA Programming in C, Prentice Hall, Upper Saddle River, NJ, 2005. 
REPORT DOCUMENTATION PAGE Form Approved 

OMB No. 0704-0188 

The 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 Department of Defense, Washington Headquarters Services, Directorate for Information Operations and Reports (0704-0188), 1215 Jefferson Davis Highway, Suite 1204, Arlington, VA 22202-4302. 

Respondents should be aware that notwithstanding any other provision of law, no person shall be subject to any penalty for failing to comply with a collection of information if it does not display a currently valid OMB 

control number. 

PLEASE DO NOT RETURN YOUR FORM TO THE ABOVE ADDRESS. 

1. REPORT DATE (DD-MM-YYYY)

01-12-2007 

2. REPORT TYPE

Final Contractor Report 

3. DATES COVERED (From - To)

4. TITLE AND SUBTITLE

Development of a Linear Stirling System Model With Varying Heat Inputs 

5a. CONTRACT NUMBER

NNC07TA38T 

5b. GRANT NUMBER

5c. PROGRAM ELEMENT NUMBER

6. AUTHOR(S)

Regan, Timothy, F.; Lewandowski, Edward, J. 

5d. PROJECT NUMBER

5e. TASK NUMBER

5f. WORK UNIT NUMBER

WBS 138494.04.01.01 

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

Sest, Inc. 

18000 Jefferson Parkway 

Middleburg Heights, Ohio 44130 

8. PERFORMING ORGANIZATION 

 REPORT NUMBER

E-16188 

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

National Aeronautics and Space Administration 

Washington, DC 20546-0001 

10. SPONSORING/MONITORS 

 ACRONYM(S)

NASA 

11. SPONSORING/MONITORING 

 REPORT NUMBER

NASA/CR-2007-215012 

12. DISTRIBUTION/AVAILABILITY STATEMENT

Unclassified-Unlimited 

Subject Categories: 20 and 44 

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

This publication is available from the NASA Center for AeroSpace Information, 301-621-0390

13. SUPPLEMENTARY NOTES

14. ABSTRACT

The linear model of the Stirling system developed by NASA Glenn Research Center (GRC) has been extended to include a user-specified heat input. Previously developed 

linear models were limited to the Stirling convertor and electrical load. They represented the thermodynamic cycle with pressure factors that remained constant. The numerical 

values of the pressure factors were generated by linearizing GRC’s nonlinear System Dynamic Model (SDM) of the convertor at a chosen operating point. The pressure factors 

were fixed for that operating point, thus, the model lost accuracy if a transition to a different operating point were simulated. Although the previous linear model was used in 

developing controllers that manipulated current, voltage, and piston position, it could not be used in the development of control algorithms that regulated hot-end temperature. 

This basic model was extended to include the thermal dynamics associated with a hot-end temperature that varies over time in response to external changes as well as to 

changes in the Stirling cycle. The linear model described herein includes not only dynamics of the piston, displacer, gas, and electrical circuit, but also the transient effects of 

the heater head thermal inertia. The linear version algebraically couples two separate linear dynamic models, one model of the Stirling convertor and one model of the thermal 

system, through the pressure factors. The thermal system model includes heat flow of heat transfer fluid, insulation loss, and temperature drops from the heat source to the 

Stirling convertor expansion space. The linear model was compared to a nonlinear model, and performance was very similar. The resulting linear model can be implemented in 

a variety of computing environments, and is suitable for analysis with classical and state space controls analysis techniques.

15. SUBJECT TERMS

Stirling engines; Design analysis; Block diagrams; Electromechanics; Dynamics models; Simulation; AC generators; Stirling cycle; 

FPGA; Simulink; VHDL 

16. SECURITY CLASSIFICATION OF: 17. LIMITATION OF 

 ABSTRACT

UU 

18. NUMBER 

 OF 

 PAGES

11 

19a. NAME OF RESPONSIBLE PERSON

STI Help Desk (email:help@sti.nasa.gov) 

a. REPORT

U 

b. ABSTRACT

U 

c. THIS 

PAGE

U 

19b. TELEPHONE NUMBER (include area code)

301-621-0390 

Standard Form 298 (Rev. 8-98)

Prescribed by ANSI Std. Z39-18