ABSTRACT
 Title of thesis: NUMERICAL AND EXPERIMENTAL STUDIES
 ON FREE PISTON STIRLING ENGINES
 Dibesh Shrestha, 2012
 Thesis directed by: Professor Balakumar Balachandran
 Department of Mechanical Engineering
 Free piston Stirling engine (FPSE) is a closed cycle engine that converts ther-
 mal energy into mechanical energy. The focus of this thesis is on understanding
 limit-cycle motions in FPSEs. First, making use of reduced-order models, paramet-
 ric studies are carried out to understand what FPSE parameters a ect the creation
 of oscillatory motions. It is shown that quasi-static variations of the sti ness and
 the damping terms on the power piston can lead to conditions for Hopf instabilities
 in the system. The e ect of the inclusion of a nonlinear spring term to the system is
 also investigated through numerical studies. The nonlinear springs include harden-
 ing springs and magnetic springs. The results, which include  rst results for FPSEs
 with hardening springs, show that nonlinear springs can help these systems realize
 limit-cycle motions. Finally, preliminary experimental studies conducted to realize
 oscillatory motions in a FPSE with a magnetic spring are also reported.
NUMERICAL AND EXPERIMENTAL STUDIES
 ON FREE PISTON STIRLING ENGINES
 by
 Dibesh Shrestha
 Thesis submitted to the Faculty of the Graduate School of the
 University of Maryland, College Park, in partial ful llment
 of the requirements for the degree of
 Master of Science
 2012
 Advisory Committee:
 Professor Balakumar Balachandran, Chair and Advisor
 Professor Amr Baz, Department of Mechanical Engineering
 Assistant Professor Nikhil Chopra, Department of Mechanical Engineering
c Copyright by
 Dibesh Shrestha
 2012
Acknowledgments
 I would like to thank Professor Balachandran, my advisor and committee
 chair, for his support and guidance during my undergraduate and graduate studies.
 I would also like to thank my committee members, Professor Baz and Professor
 Chopra, for their assistance in completing this thesis.
 I want to thank Vince Nguyen for aiding me in my research. I also want to
 thank all of my other colleagues in the Dynamics and Vibrations Group for their
 support and help. Finally, I would like to thank my family and friends for their
 support and love.
 ii
Table of Contents
 List of Figures iv
 1 Introduction 1
 1.1 Con gurations of Stirling Engines . . . . . . . . . . . . . . . . . . . . 2
 1.2 Existing Literature on FPSE Dynamics . . . . . . . . . . . . . . . . . 6
 1.3 Contributions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
 1.4 Thesis Organization . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
 2 Governing Equations 11
 2.1 Dynamic Equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
 2.2 Schmidt Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
 2.3 System with Gas Spring . . . . . . . . . . . . . . . . . . . . . . . . . 17
 2.4 Pressure Loss Term . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
 3 Parametric Studies 22
 3.1 Parametric Investigations using Root Locus Plots . . . . . . . . . . . 25
 3.2 E ect of Nonlinear Spring . . . . . . . . . . . . . . . . . . . . . . . . 29
 4 Governing Equations for System with Magnetic Springs 40
 4.1 Nondimensional Equations . . . . . . . . . . . . . . . . . . . . . . . . 43
 4.2 Repulsive Magnetic Forces . . . . . . . . . . . . . . . . . . . . . . . . 45
 4.3 Results for System Response . . . . . . . . . . . . . . . . . . . . . . . 49
 5 Experimental Studies 54
 5.1 Design and Fabrication . . . . . . . . . . . . . . . . . . . . . . . . . . 55
 5.2 Preliminary Experimental Results . . . . . . . . . . . . . . . . . . . . 58
 6 Summary and Concluding Remarks 61
 6.1 Parametric Studies . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62
 6.2 System with Magnetic Springs . . . . . . . . . . . . . . . . . . . . . . 63
 6.3 Experimental Studies . . . . . . . . . . . . . . . . . . . . . . . . . . . 63
 6.4 Future Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64
 A Equations of Motion and Parameters from Ulusoy (1994) 65
 B Experimental Parameters 69
 C Matlab Codes 72
 iii
List of Figures
 1.1 An example of Stirling engine con guration. (Source: Urieli, 2010) . . 2
 1.2 Di erent types of Stirling engines: (a) Kinematically linked Stirling
 engine (Source:http://en.wikipedia.org/wiki/Stirling_engine).
 (b) Free Piston Stirling engine (Source: http://www.bekkoame.ne.
 jp/\protect\unhbox\voidb@x\penalty\@M\{}khirata/english/fpse.
 htm). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2
 1.3 Cycles for kinematic Stirling engines (Source: http://en.wikipedia.
 org/wiki/Stirling_engine). . . . . . . . . . . . . . . . . . . . . . . 3
 1.4 The cyclic process of free piston Stirling engines (Source: Ulusoy, 1994). 4
 1.5 Di erent con gurations of Stirling engines: (a) Alpha type con gu-
 ration. (b) Beta type con guration. (c) Gamma type con guration
 (Source: Urieli, 2010). . . . . . . . . . . . . . . . . . . . . . . . . . . 6
 2.1 Schematic of a beta type FPSE. . . . . . . . . . . . . . . . . . . . . . 12
 2.2 Assumed temperature pro le of regenerator (Source: Urieli, 2010). . . 14
 2.3 Schematic of a gas spring beta type FPSE. . . . . . . . . . . . . . . . 18
 3.1 Schematic of a gas spring beta type FPSE. . . . . . . . . . . . . . . . 22
 3.2 Root locus plot with spring parameter, ~kp, set to zero while linear
 load parameter, ~fl, is varied from -1.00 to 0.00. . . . . . . . . . . . . 27
 3.3 Root locus plot with linear load parameter, ~fl, set to zero while spring
 parameter, ~kp, is varied from 1.00 to 0.00. . . . . . . . . . . . . . . . 28
 3.4 Root Locus Plot with varying ~fl and ~kp. (a) ~fl: -1.00 to 0.00 and
 ~kp=0.20. (b) ~fl: -1.00 to 0.00 and ~kp=0.30. (c) ~fl: -1.00 to 0.00 and
 ~kp=0.60. (d) ~fl: -1.00 to 0.00 and ~kp=0.80. . . . . . . . . . . . . . . . 29
 3.5 Displacement response of the displacer and the power piston ex-
 cluding nonlinear terms. The chosen parameters are ~fl = 0:00 and
 ~kp = 0:10. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
 3.6 Displacement response of the displacer and the power piston ex-
 cluding nonlinear terms. The chosen parameters are ~fl = 0:10 and
 ~kp = 0:10. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32
 3.7 Displacement responses for the displacer and the power piston due to
 additional nonlinear spring added to the power piston. The parame-
 ters are ~kp = 0:10, ~kc = 0:10, and ~fl = 0:00. (a) ICs: (1.00 0.00 1.00
 0.00) and (b) ICs: (7.00 0.00 7.00 0.00). . . . . . . . . . . . . . . . . . 33
 3.8 Phase portraits with chosen initial conditions of (1.00 0.00 1.00 0.00)
 and (7.00 0.00 7.00 0.00). (a) Response of displacer and (b) Response
 of power piston. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34
 3.9 Displacement response comparison of cases with low sti ness nonlin-
 ear spring and without nonlinear spring. (a) Response of displacer
 and (b) Response of power piston. . . . . . . . . . . . . . . . . . . . . 36
 iv
3.10 Displacement response comparison of cases with high sti ness non-
 linear spring and without nonlinear spring. (a) Response of displacer
 and (b) Response of power piston. . . . . . . . . . . . . . . . . . . . . 37
 3.11 Frequency spectrum comparison for cases with and without nonlinear
 spring. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38
 3.12 Amplitude response variation with respect to linear sti ness param-
 eter ~kp, and nonlinear sti ness parameter ~kc, with linear load param-
 eter, fl, set to -0.02. (a) Response of displacer and (b) Response of
 power piston. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39
 4.1 Schematic of magnetic spring FPSE. . . . . . . . . . . . . . . . . . . 40
 4.2 Magnetic spring force curve compared to various types of spring forces. 42
 4.3 Magnets modeled as charges. . . . . . . . . . . . . . . . . . . . . . . . 46
 4.4 Repulsive force between two similar magnets. . . . . . . . . . . . . . 48
 4.5 Displacement response of magnetic spring FPSE without any damp-
 ing or head loss. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50
 4.6 Phase portraits for the case with no damping and with ICs: (0.01 0.0
 0.0 0.0). (a) Displacer and (b) Power Piston. . . . . . . . . . . . . . 50
 4.7 Displacement response of magnetic spring FPSE with low damping
 and head loss. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51
 4.8 Phase portraits for the case with low damping and head loss and with
 ICs: (0.01 0.0 0.0 0.0). (a) Displacer and (b) Power Piston. . . . . . 51
 4.9 Displacement response of magnetic spring FPSE with high damping
 and head loss. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52
 4.10 Displacement response of magnetic spring FPSE with high damping
 and head loss with temperature di erence of 500:00 Celcius. . . . . . 52
 5.1 Kinematic Stirling engine to be modi ed into FPSE. . . . . . . . . . 54
 5.2 Schematic of the experimental setup. . . . . . . . . . . . . . . . . . . 55
 5.3 Modi ed FPSE. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56
 5.4 Increased height for the power piston. . . . . . . . . . . . . . . . . . . 57
 5.5 Increased volume for the expansion space. . . . . . . . . . . . . . . . 58
 5.6 Power piston displacements for a kinematically linked Stirling engine. 59
 5.7 Power piston displacements for  rst FPSE arrangement. . . . . . . . 59
 5.8 Power piston displacements for second FPSE arrangement. . . . . . . 60
 A.1 Schematic of FPSE from Ulusoy (1994). . . . . . . . . . . . . . . . . 65
 v
Chapter 1
 Introduction
 The Stirling engine, invented by Robert Stirling in 1816, is a device that can
 be used to convert thermal energy into mechanical energy by taking advantage of
 the pressure variations within the system developed through a temperature gradi-
 ent. This is a closed cycle regenerative engine, which compresses a  xed amount
 of gas through cooling and then expands the gas again through heating; this pro-
 cess is repeated cyclically. A schematic of an alpha type Stirling engine is shown
 in Figure 1.1, which shows the components of a Stirling engine. In order to create
 the temperature gradient, the working space of the engine is separated into two
 di erent compartments: the compression space where the working gas is cooled,
 and the expansion space where the working  uid is heated. Other necessary compo-
 nents for a Stirling engine are a cooling source to cool the  uid in the compression
 space, a heating source to heat the expansion space, a displacer piston, and a power
 piston. The main purpose of the pistons are to move the working  uid back and
 forth between the two compartments, which is further explained later in this chap-
 ter. Adding a regenerator, which acts as a heat exchanger, is said to increase the
 e ciency of a Stirling engine. A regenerator retains heat from the working  uid and
 then redelivers the heat again to the  uid as the  uid travels back and forth between
 the two compartments (Organ, 1992). The placement of the various components of
 1
Figure 1.1: An example of Stirling engine con guration. (Source: Urieli, 2010)
 the Stirling engine de ne the type of Stirling engine, which will be also discussed
 further in this chapter.
 1.1 Con gurations of Stirling Engines
 Stirling engines can be classi ed on the basis of their driving mechanisms,
 namely, constrained drive systems such as kinematic Stirling engines and uncon-
 strained drive systems like Free Piston Stirling engines (FPSE). Schematics for the
 kinematic Stirling engine and FPSE are shown in Figure 1.2. As shown in Figure
 1.2(a), the motions of the kinematic Stirling engines are de ned by the crankshafts
 and the rods that connect the power piston and the displacer. The power piston
 (a) (b)
 Figure 1.2: Di erent types of Stirling engines: (a) Kinematically linked Stirling en-
 gine (Source:http://en.wikipedia.org/wiki/Stirling_engine). (b) Free Piston
 Stirling engine (Source: http://www.bekkoame.ne.jp/~khirata/english/fpse.
 htm).
 2
and the displacer are linked so that there is a phase di erence with respect to one
 another, and the resulting system is a single degree-of-freedom (SDOF) system. On
 the other hand, FPSEs are unconstrained such that the piston and the displacers
 move solely due to the pressure variation of the working gas caused by the temper-
 ature di erence. This arrangement results in independent movement of the piston
 and displacer, allowing the system to be a multi-degree-of-freedom (MDOF) sys-
 tem, as depicted in Figure 1.2(b). A FPSE has the following advantages over the
 kinematic Stirling engine: (i) can be produced at a low cost, (ii) does not require
 external high pressure seal, (iii) long operating life without the need for lubrication,
 and (iv) quiet operating conditions.
 The cyclic process for the movement of a kinematic Stirling engine is depicted
 in Figure 1.3. For a kinematic Stirling engine, initial conditions are given to the
 power piston to compress the working  uid in the compression space (Step 1). The
 working  uid is heated, increasing the pressure within the system. The heated
  uid increases the pressure and pushes the power piston to its upstroke (Step 2).
 (a) Step 1 (b) Step 2 (c) Step 3 (d) Step 4
 Figure 1.3: Cycles for kinematic Stirling engines (Source: http://en.wikipedia.
 org/wiki/Stirling_engine).
 3
The displacer, which is kinematically linked with the power piston with a phase
 di erence, moves the heated working  uid to the compression space (Step 3). Finally,
 the cooled gas is compressed by the power piston due to its momentum along with
 the decreased pressure (Step 4).
 On the other hand, a FPSE relies on the pressure variations for its operation,
 as depicted in Figure 1.4. When heat is applied, the power piston of FPSE is pushed
 outwards by the expanding gas. The increase in volume decreases pressure in the
 compression space, which creates a pressure di erence across the displacer. The
 Figure 1.4: The cyclic process of free piston Stirling engines (Source: Ulusoy, 1994).
 4
working  uid moves from the expansion space to the compression space, forcing
 the displacer towards the expansion space (Step 1). The power piston is almost
 stationary; hence, the volume within the system can be considered constant. The
 working gas is cooled at a constant volume, reducing the pressure of the gas. The
 reduced pressure causes the power piston to move towards the hot end by its own
 inertia, compressing the cold working  uid in the compression space (Steps 1 and
 2). This increases the pressure in the compression space causes the working  uid
 to  ow to the expansion space. This results in displacer moving towards the power
 piston and thus collapsing the compression space (Steps 2 and 3). As the working
  uid heats up, the pressure of the system increases and causes the power piston to
 move outward (Steps 3 and 4), leading to a cyclic process (Walker, 1980).
 Stirling engines can be further classi ed on the basis of their con guration as
 alpha, beta, and gamma type engines (Urieli and Berchowitz, 1984; Urieli, 2010).
 The three aforementioned con gurations are illustrated in Figure 1.5. The alpha
 engine, shown in Figure 1.5(a), contains two pistons in two di erent cylinders at
 the opposite ends that are connected by a heater, regenerator, and cooler in series.
 This is the simplest type of a Stirling engine con guration but the pistons have to
 be sealed for proper functionality. Unlike the alpha con guration, the beta con gu-
 ration, shown in Figure 1.5(b), only has one cylinder with a hot end and a cool end.
 The displacer sits inside the cylinder while the power piston is above the compres-
 sion space. The gamma con guration, shown in Figure 1.5(c), is similar to the beta
 type except that the power piston and displacer are in two di erent cylinders. The
 cycles for all three con gurations will depend upon whether it is a kinematic or free
 5
(a)
 (b)
 (c)
 Figure 1.5: Di erent con gurations of Stirling engines: (a) Alpha type con guration.
 (b) Beta type con guration. (c) Gamma type con guration (Source: Urieli, 2010).
 piston Stirling engine, as described previously.
 1.2 Existing Literature on FPSE Dynamics
 Various studies on the FPSE are based on the Schmidt?s analysis, as this ap-
 proach can lead to a closed form solution for the engine performance. The Schmidt
 model helps relate the thermodynamics to the dynamics of the FPSE. The assump-
 tions made in the model are as follows: (i) isothermal compression and expansion,
 (ii) ideal gas behavior, (iii) spatially constant pressure, and (iv) closed system en-
 6
gine. By using the aforementioned assumptions, the working gas pressure can be
 determined in terms of volume variations.
 In 1985, Berkowitz and Redlich did a linear analysis of FPSE using the Schmidt
 isothermal model to determine the requirements for oscillation and general behavior
 under load. The linear analysis of the system showed that the equilibrium point
 near the origin had one pair of complex eigenvalues with a zero or small positive
 real part and another pair with a negative real part. The growing oscillations need
 to be controlled by introducing nonlinearities to the system, such as gas springs.
 It is suggested to change the damping coe cients by using a control mechanism to
 push the positive real part of the conjugate pair towards the imaginary axis. It is
 also found that there is a minimum hot end temperature for which oscillations can
 be expected. By making the frequency of the displacer oscillations to match the
 frequency of the power piston oscillations, the engine can start at a lower hot end
 temperature and power can be maximized (Berchowitz and Redlich, 1985).
 In 1990, Benvenuto, de Monte, and Farina studied a methodology for design
 optimization and performance evaluation of FPSE dynamics for space applications.
 The e ects of temperature variations of the working spaces on the spring terms
 along with gas hysteresis losses for the bu er space and for the gas spring space are
 included as linear damping terms. The bu er space is the nominal space under the
 power piston, and the spring space is the air space that acts as a spring attached
 to the displacer. The developed method allows analytical relations for determining
 dynamic behavior for their modeled systems. The model also allows for minimizing
 hysteresis losses in the gas springs and viscous losses in the heat exchangers. Finally,
 7
an analytical relationship was constructed for the phase angle between the displacer
 and the piston motions in order to optimize engine e ciency (Benvenuto and de
 Monte, 1995).
 In 1994, Ulusoy investigated the nonlinear e ects on FPSE through isothermal
 and nonisothermal modeling. The thermodynamics and dynamics of the system are
 coupled via isothermal and nonisothermal methods. The e ects of nonlinear damper
 load, nonlinear pressure loss, and a gas spring acting on the displacer were studied.
 By using nonlinear analysis, it was shown that periodic motions of the piston and
 displacer can be attained and controlled by using the coe cient of the nonlinear load
 term. The e ect of the temperature variation due to nonisothermal behavior of the
 working gas on FPSE was also investigated by using numerical means. The results
 showed that variation in the working gas temperature does not have a signi cant
 e ect on the dynamics of the engine, leading to the conclusion that an isothermal
 assumption is appropriate for a qualitative dynamic analysis (Ulusoy, 1994).
 In 2009, at the University of Maryland, in the same research group as this
 thesis author, Choudhary investigated how to engineer a Hopf bifurcation of an
 equilibrium solution in a FPSE via nonlinear analysis and simulations. Reduced-
 order models were developed on the basis of Schmidt and nodal analysis. Through
 eigenvalue analysis, it was shown that introducing an appropriate nonlinearity into
 the system can lead to a Hopf bifurcation, the result of which can be an attracting
 limit cycle. The method of multiple scales was used to study the weakly nonlinear
 system analytically in order to develop limit-cycle motions in a beta FPSE with
 cubic damping on the power piston side. The analytical prediction was compared to
 8
numerical results, which show that the analytical and numerical solutions are close
 to one another near the bifurcation point. Finally, a simpli ed nodal method was
 developed to reduce the order of FPSE model that captures more of the thermody-
 namic behavior than the Schmidt model, while making the analysis of the dynamics
 possible. However, the simpli ed nodal methods presented in the thesis study did
 not yield physically valid models (Choudhary, 2009). This thesis builds upon the
 studies conducted by Choudhary. Here, further parametric tests are conducted to
 evaluate various parameters that can introduce Hopf instabilities. Additional non-
 linearities that can produce oscillatory motions are introduced into the system and
 studied. Finally, experimental studies are done.
 1.3 Contributions
 A main purpose of this study is to better understand oscillatory motions in
 a FPSE by introducing nonlinearities into the system. A FPSE system with gas
 spring for the displacer and a nonlinear spring for the power piston are introduced.
 The introduction of the spring parameters gives one the ability to change the design
 frequency of the FPSE without much trouble when the masses of the pistons are
 constrained. This system is linearized to observe which parameters, in particular, the
 spring parameter and the damping parameter, a ect the Hopf instabilities. Through
 root locus analysis, the in uences of the sti ness of the system and the damping
 are assessed. Next, the e ect of cubic nonlinear spring addition to the power piston
 is studied through numerical studies to see if limit-cycle motions can be attained.
 9
This  nding can be applied to FPSEs without cubic damping in order to introduce
 limit-cycle motions. Once this was completed, magnetic springs are introduced in
 the displacer and power pistons to study if the FPSE system can exhibit limit-cycle
 motions. A study of the e ects of the magnetic springs on the FPSE is carried out
 to see if this is a viable option for application. After the magnetic spring model
 showed an oscillatory motion, an experimental setup is fabricated and the motions
 of the displacer and piston are studied. A recent study conducted on FPSEs is the
 study of Formosa (2009).
 1.4 Thesis Organization
 In the next chapter, the governing equations are presented. First, the Schmidt
 analysis is provided and the general form of the FPSE governing equations is givem.
 Second, a spring parameter is introduced into the gas spring FPSE studied by Ulusoy
 (Ulusoy, 1994). The pressure losses within the FPSE are also derived. In Chapter
 3, a parametric study is presented for the various cases of the modi ed gas spring
 FPSE. Afterwards, the governing equations for the magnetic spring system are de-
 rived in Chapter 4 and numerical studies carried out for the magnetic spring FPSE
 are presented. The fabrication details and experimental studies for the FPSE are
 presented in Chapter 5. Summary and future work are provided in Chapter 6 along
 with concluding remarks. To close this thesis, in the appendices, some details of
 Ulusoy?s work , the experimental parameters, and the codes used in the work are
 provided.
 10
Chapter 2
 Governing Equations
 The governing equations of motions for the displacer and power piston are
 obtained by using force balance for a beta type FPSE. However, this methodology
 can be used for any other type of FPSE as well. By using Schmidt analysis, the
 thermodynamics of the system is described as algebraic functions of the displace-
 ments and velocities of the displacer and power piston. Subsequently, a nonlinear
 spring term is added to the power piston of the gas spring FPSE studied by Ulusoy
 (Ulusoy, 1994). The equations of motion for the modi ed case are provided along
 with the parameters used in the studies carried out in Chapter 3. The pressure loss
 for a general type of FPSE is also derived in this chapter.
 2.1 Dynamic Equations
 The dynamic equations for the schematic of FPSE con guration shown in
 Figure 2.1 are derived. The FPSE can be considered as a two DOF system governed
 by force balance of pressure, spring, and external load forces. The net force acting
 on the displacer piston and the power piston is a resultant of the pressure di erence
 between the two sides of the pistons. Hence, applying Newton?s second law, the
 11
equations of motion for the given FPSE can be written as:
 md xd = Ad(Pc  Pe) Fspring1 (2.1)
 mp xp = Ap(Pb  Pc) + Fload _xp  Fspring2 (2.2)
 where Ad is the area of the displacer, Ap is the area of the piston, md is the mass
 of the displacer, mp is the mass of the piston, Vb is the nominal volume under the
 power piston, Pb is the nominal pressure under the power piston, and Fload is an
 external load that is attached to the power piston. The numerical values given for
 Figure 2.1: Schematic of a beta type FPSE.
 12
Fload should be negative because since this load parameter represents energy being
 taken out of the system by equipment such as an alternator.
 To solve equations (2.1) and (2.2), one requires information about the pres-
 sures in the expansion, compression, and bu er space. The working gas pressure is
 dependent on the displacements of the piston and the displacer. In the following
 section, pressure is related to the volume via Schmidt analysis.
 2.2 Schmidt Analysis
 In 1871, an isothermal analysis of the Stirling engine was completed by Gus-
 tav Schmidt. Major assumptions of this analysis are isothermal compression and
 expansion, perfect regeneration, and spatially constant instantaneous working gas
 pressure. The Schmidt analysis provides a simple way to relate the dynamics and
 thermodynamics of the system.
 The pressure terms are modeled as a function of the volume and the temper-
 ature. For this, the conservation of mass and the ideal gas law are used. By using
 the conservation of mass, the total mass of working gas inside the engine can be
 expressed as (Urieli and Berchowitz, 1984)
 mt = me +mh +mR +mk +mc (2.3)
 where mt is the total mass of the working gas, me is the mass of the gas inside the
 expansion space, mh is the mass of the gas inside the heater, mR is the mass of the
 gas inside the regenerator, mk is the mass of the gas inside the cooler, and mc is the
 13
mass of the gas in the compression space.
 An assumption of instantaneous working gas pressure is made, which leads to
 P  Pe  Ph  PR  Pk  Pc (2.4)
 It is also assumed that the temperature of expansion space is equal to the heater
 temperature and the temperature of the compression space is equal to the cooler
 temperature:
 Te = Th and Tc = Tk (2.5)
 Hence,using the ideal gas relation: mt =
 X PV
 RT
 , the total mass of the working gas
 can be written as
 mt =
 P
 R
  
Ve
 Th
 +
 Vh
 Th
 +
 VR
 TR
 +
 Vk
 Tk
 +
 Vc
 Tk
  
(2.6)
 where Ve, Vh, VR, Vk, and Vc are the respective volumes. The temperature along the
 regenerator is assumed to be linear as shown in Figure 2.2. Hence, the temperature
 Figure 2.2: Assumed temperature pro le of regenerator (Source: Urieli, 2010).
 14
in the regenerator can be extrapolated as
 T (x) =
 (Th  Tk)x
 Lr
 + Tk (2.7)
 where Lr is the regenerator length. Now, the total mass of the gas in the regenerator
 can be found by using ideal gas relationship:
 mr =
 VrP
 R
 Z Lr
 0
 dx
 (Th  Tk)x+ TkLr
 (2.8)
 mr =
 VrP
 R
 ln(Th=Tk)
 (Th  Tk)
 (2.9)
 Hence, the regenerator temperature is found to be
 TR =
 Th  Tk
 ln (Th=Tk)
 (2.10)
 Next, the compression volume and expansion volume is de ned by the geom-
 etry of the system; that is,
 Vc = Vcm  Apxp + Adxd (2.11)
 Ve = Vem  Adxd (2.12)
 where Vcm is the mean compression volume and Vem is the mean expansion volume.
 15
Next, new parameters are introduced for the simpli cation:
 S =
 Vem
 Th
 +
 Vh
 Th
 +
 VR ln
  
Th
 Tk
  
Th  Tk
 +
 Vk
 Tk
 +
 Vcm
 Tk
 (2.13)
 a =
 Ap
 Tk
 1
 S
 (2.14)
 b =
  
Ad
 Tk
  
Ad
 Th
  
1
 S
 (2.15)
 On substituting equations (2.10), (2.11), and (2.12) into equation (2.6) and using
 equations (2.13), (2.14), and (2.15), the pressure term can be determined as
 P =
 mtR
 S (1 axp + bxd)
 (2.16)
 It is assumed that the pressure of the working gas is equal to the mean pressure,
 P = Pm, when the displacer and the piston are at mid-stroke; that is xp = 0 and
 xd = 0. Hence, the mean pressure can be expressed as
 Pm =
 mtR
 S
 (2.17)
 By using equation (2.17), the working space pressure can be expressed as
 P =
 Pm
 1 axp + bxd
 (2.18)
 It is assumed that the compression pressure is equal to the working gas pressure,
 and the expansion pressure is the compression pressure minus pressure drop ( P )
 16
across the heater, regenerator, and cooler. It is also assumed that the bu er pressure
 is equal to the mean pressure (Ulusoy, 1994).
 Pc = P (2.19)
 Pe = Pc   P (2.20)
 Pb = Pm (2.21)
 Hence, the equations of motions become
 md xd = Ad P  Fspring1 (2.22)
 mp xp = ApPm
  
1 
1
 F
  
+ Fload _xp  Fspring2 (2.23)
 F = 1 axp + bxd (2.24)
 2.3 System with Gas Spring
 In 1994, Ulusoy studied the e ects of nonlinearity in a FPSE. One such study
 was the modeling of gas spring to introduce nonlinearity into the system, as shown
 in Figure 2.3. In this section, the equations of motions are presented after adding a
 spring parameter to the power piston of the gas spring FPSE studied by Ulusoy.
 The gas spring is modeled by assuming ideal gas relation for an adiabatic
 17
process and de ning the gas spring volume (Ulusoy, 1994)
 Ps = Pm
  
Vsm
 Vs
   
(2.25)
 Vs = Vsm + Arxd (2.26)
 where Vs is the gas spring volume, Vsm is its average value, and Ar is the area of the
 gas spring rod. Introducing a new quantity
 c =
 Ar
 Vsm
 (2.27)
 Figure 2.3: Schematic of a gas spring beta type FPSE.
 18
the gas spring pressure term can be written as:
 Ps = Pm
  
1
 1 + cxd
   
(2.28)
 Also, the b parameter in equation (2.24) should have the area of gas spring rod (Ar)
 taken into account as well. This means
 b =
  
Ad  Ar
 Tk
  
Ad
 Th
  
1
 S
 (2.29)
 After the spring parameter is added to the power piston as shown in Figure
 2.3, the equations of motion for the gas spring FPSE become
 md xd = Ad (Pc  Pe) Ar (Pc  Ps) (2.30)
 mp xp = Ap (Pb  Pc) + Fload _xp  Fspring (2.31)
 2.4 Pressure Loss Term
 The pressure drop term,  P , can be expressed as
  P =  pcooler +  pregenerator +  pheater (2.32)
 In this thesis study, including the experimental setup, the FPSE does not include
 a regenerator; hence, the pressure loss in the regenerator is ignored. Without the
 regenerator, the pressure drop can be modeled as a continuous pressure drop through
 19
one channel assuming the geometry of the channel does not vary. Thus, the pressure
 drop is modeled only as a drop through the cooler area. It is assumed that the
 working  uid acts turbulently inside the FPSE. Hence, the pressure drop derived by
 Urieli is used (Urieli and Berchowitz, 1984)
  Pcooler =
 1
 2
  (ft + kh)ucoolerjucoolerj (2.33)
 where  is the density of the working  uid, kh is the head loss, u is the velocity of
 the working  uid, and ft is the turbulent friction. This factor is given by
 ft =
 4ffL
 dh
 (2.34)
 where ff is the Fanning friction coe cient factor, L is the length of the  ow passage,
 and dh is the hydraulic diameter. The Fanning friction coe cient factor is de ned
 as
 ff = CfRe
 n (2.35)
 where Re is the Reynolds number de ned as
 Re =
  Udh
  
(2.36)
 and  is the viscosity of the working  uid. Some commonly used working gases?
 viscosities can be calculated by using relations given by Martini (1983). The rec-
 ommended values for Cf and n depend on the details of the  ow.
 20
Next, the gas  ow velocities need to be related to those of the reciprocating
 elements. This is done by de ning the volumetric  ow rate of the compression and
 expansion space through
 _Vc =
 dVc
 dt
 (2.37)
 _Ve =
 dVe
 dt
 (2.38)
 Taking the respective term derivatives of equations (2.10) and (2.11), the following
 relationships are found:
 _Vc = Ad _xd  Ap _xp (2.39)
 _Ve =  Ad _xd (2.40)
 Next, the velocity of gas with respect to the velocities of the pistons are ap-
 proximated as
 ucooler =
 _V
 Acooler
 (2.41)
 where _V = _Ve _Vc because positive values of _Ve and _Vc indicate increasing volumes,
 which results in
 _V = Ap _xp  2Ad _xd (2.42)
 Finally, inserting equation (2.42) into equation (2.41), and then into equation
 (2.32), the pressure drop can be written as
  P =
 1
 2
  (ft + kh)
  
Ap _xp  2Ad _xd
 Ak
 j
 Ap _xp  2Ad _xd
 Ak
 j
  
(2.43)
 21
Chapter 3
 Parametric Studies
 In this chapter, parametric studis are pursued by using the model presented in
 the previous chapter. This model is shown in Figure 3.1. The numerical studies help
 understand which parameters push the Stirling engine to exhibit Hopf instabilities.
 The equations of motions reported in Ulusoy?s study are given in Appendix A. In
 that study, the equations of motion are rewritten by  rst nondimensionalizing the
 Figure 3.1: Schematic of a gas spring beta type FPSE.
 22
equations and then expanding the nonlinear terms by using Taylor series. Also, the
 pressure loss terms are separated as a combination of linear and nonlinear terms
 (Ulusoy, 1994). For the purposes of this study, a spring parameter is introduced
 into the power piston to study Hopf instabilities. The signi cance of the new spring
 parameter becomes apparent for FPSEs without cubic damper, as this nonlinearity
 can help realize oscillatory motions. Additionally, one can change the oscillation
 frequency of the system through the spring parameter. The frequency response of
 the system is determined by the displacer and piston masses. However, in practice,
 these parameters cannot be easily changed. On the other hand, one can easily change
 the frequency of the sytem through springs to get a desired frequency associated with
 a periodic motion..
 The equations of motions of the modi ed syste have been reported in Chapter
 2.3 and are presented in equations (2.30) and (2.31). Next, the equations are put in
 state-space form for the system with the gas spring:
 _x =
 2
 6
 6
 6
 6
 6
 6
 6
 6
 6
 6
 4
 0 1 0 0
 k1 c1 k2 c2
 0 0 0 1
 k3 0 k4  
~kp
 ~mp
 ~fl
 ~mp
 3
 7
 7
 7
 7
 7
 7
 7
 7
 7
 7
 5
 2
 6
 6
 6
 6
 6
 6
 6
 6
 6
 6
 4
 x1
 x2
 x3
 x4
 3
 7
 7
 7
 7
 7
 7
 7
 7
 7
 7
 5
 +
 2
 6
 6
 6
 6
 6
 6
 6
 6
 6
 6
 4
 0
  1HOT1 +  2HOT2
 0
  
~fc
 ~mp
 x23x4 +  3HOT3
 3
 7
 7
 7
 7
 7
 7
 7
 7
 7
 7
 5
 (3.1)
 23
Here,
 ~xd = x1 (3.2)
 _~xd = x2 (3.3)
 ~xp = x3 (3.4)
 _~xp = x4 (3.5)
 ~kp =
 lpkp
 2ApPm
 (3.6)
 ~fl =
 Fload
 ApPm
 (3.7)
 For the parametric studies, the parameters given by Ulusoy are used. The
 pressure loss model is given by (Ulusoy, 1994)
  P =  0:028x2  0:511x
 3
 2 + 0:021x4 + 1:197x
 2
 2x4  0:979x2x
 2
 4 + 0:269x
 3
 4
 Finally, the derived state-space model is given by
 _x =
 2
 6
 6
 6
 6
 6
 6
 6
 6
 6
 6
 4
 0:00 1:00 0:00 0:00
  0:38  1:60  1:30 1:20
 0:00 0:00 0:00 1:00
 0:62 0:00  1:06 
~kp
 0:26
 ~fl
 0:26
 3
 7
 7
 7
 7
 7
 7
 7
 7
 7
 7
 5
 2
 6
 6
 6
 6
 6
 6
 6
 6
 6
 6
 4
 x1
 x2
 x3
 x4
 3
 7
 7
 7
 7
 7
 7
 7
 7
 7
 7
 5
 +
 2
 6
 6
 6
 6
 6
 6
 6
 6
 6
 6
 4
 0
  1HOT1 +  2HOT2
 0
  
~fc
 0:26
 x23x4 +  3HOT3
 3
 7
 7
 7
 7
 7
 7
 7
 7
 7
 7
 5
 (3.8)
 where
 HOT1 =  26:21x
 3
 2 + 66:38x
 2
 2x4  53:21x2x
 2
 4 + 15:31x
 3
 4
 24
HOT2 = 0:002x
 2
 1 0:001x
 3
 1 0:007x1x3 0:002x
 2
 1x3 0:006x
 2
 3 0:003x1x
 2
 3 0:002x
 3
 3
 HOT3 = 0:025x
 2
 1  0:004x
 3
 1  0:085x1x3  0:02x
 2
 1x3  0:074x
 2
 3  0:035x1x
 2
 3  0:02x
 3
 3
 3.1 Parametric Investigations using Root Locus Plots
 The system of equations is solved for the equilibrium points, which is found to
 be at the origin (0.00,0.00,0.00,0.00). In order to study the behavior of the system,
 the Jacobian of the system at the equilibrium point is investigated:
 [Dif ] =
 2
 6
 6
 6
 6
 6
 6
 6
 6
 6
 6
 4
 0:00 1:00 0:00 0:00
  0:38  1:60  1:30 1:20
 0:00 0:00 0:00 1:00
 0:62 0:00  1:06 
~kp
 0:26
 ~fl
 0:26
 3
 7
 7
 7
 7
 7
 7
 7
 7
 7
 7
 5
 (3.9)
 The characteristic equation is derived from the Jacobian and used to study the
 instabilities of the system are studied. A parametric study is done by varying ~kp,
 the sti ness of spring, and ~fl, the linear damping parameter associated with power
 piston. The goal is to realize a system that exhibits limit-cycles through a Hopf
 bifurcation. The criteria for a Hopf bifurcation of an equilibrium point are given by
 (Nayfeh and Balachandran, 1995)
 i) The system has an equilibrium point x = x0 at the critical point  =  c.
 ii) The Jacobian has a pair of purely imaginary eigenvalues while all of the other
 eigenvalues have non-zero real parts at (x0;  c).
 25
iii) For   =  c, let the analytic continuation of the pair of imaginary eigenvalues
 be   i!. Then
 d 
d 
6= 0 at   =  c.
 In this thesis work, the parameters ~fl and ~kp are synonymous with the two param-
 eter. Condition (i) is met by solving for the equilibrium point. After using the
 approximate values of -0.10 and 0.10 for ~fl and ~kp, respectively. It is found that the
 system is "close" to satysfying the Condition (ii); the associated eigenvalues of the
 Jacobian matrix are
  1;2 =  0:99 0:76i
  3;4 =  0:94i
 Condition (iii) requires the eigenvalues to cross the imaginary axis with nonzero
 speed, called the transversality condition. When the eigenvalues cross the imaginary
 axis, then, there must be a point when the Jacobian matrix has a pair of purely
 imaginary eigenvalues. In order to check the transversality condition, the root locus
 plots of the system are observed as the parameters ~kp and ~fl are varied. First, the
 value of the spring parameter, ~kp, is set to zero and the value of the linear load
 parameter, ~fl, is varied. Next, ~fl is set to zero and ~kp is varied. Finally, both ~kp and
 ~fl are varied to observe how the root locus plots change. Through the numerical
 studies, one can examine which parameters are important and can be tuned to
 introduce Hopf instabilities.
 Since it can be assumed that the load is dissipating energy from the system,
 the linear load parameter ( ~fl) is varied from -1.00 to 0.00 and spring term is set to
 26
?3.5 ?3 ?2.5 ?2 ?1.5 ?1 ?0.5 0 0.5?1
 ?0.8
 ?0.6
 ?0.4
 ?0.2
 0
 0.2
 0.4
 0.6
 0.8
 1
 Root Locus, kp=0
 Real
 Imaginar
 y
 Figure 3.2: Root locus plot with spring parameter, ~kp, set to zero while linear load
 parameter, ~fl, is varied from -1.00 to 0.00.
 0.00 for generating the plot shown in Figure 3.2. The values of ~fl are only varied upto
 0.00 since positive values would mean that energy is being added to the system. It
 can be seen that one pair of the eigenvalues crosses the imaginary axis, whereas the
 other pair has a negative real part. The system has only one pair of eigenvalues that
 becomes purely imaginary whereas the second pair has real and imaginary parts,
 which meets the Hopf instability condition (ii). It can also be seen from the graph
 that the transversality condition for a Hopf instability is met. The value for ~fl at
 which there is a pair of purely imaginary eigenvalues is approximately -0.22, which
 is the same as what Ulusoy (1994) had determined. This is expected since when ~kp
 is 0.00, both systems are identical.
 Next, the spring parameter is varied from 1.00 to 0.00, while the linear load
 parameter is set to 0.00, as shown in Figure 3.3. It can be seen that a pair of
 eigenvalues cross the imaginary axis while the other pair has a negative real part.
 27
?1.2 ?1 ?0.8 ?0.6 ?0.4 ?0.2 0 0.2 0.4?2.5
 ?2
 ?1.5
 ?1
 ?0.5
 0
 0.5
 1
 1.5
 2
 2.5
 Root Locus, fl=0
 Real
 Imaginar
 y
 Figure 3.3: Root locus plot with linear load parameter, ~fl, set to zero while spring
 parameter, ~kp, is varied from 1.00 to 0.00.
 Like the  rst case, this case meets both the Hopf instability conditions (ii) and
 (iii). The value of ~kp at which there is a pair of purely imaginary eigenvalues is
 approximately 0.40.
 Finally, both ~fl is varied for various values of ~kp to determine their combined
 e ect on the system. From Figure 3.4, it can be seen that as ~kp increases, the
 conditions for a Hopf instability are not met. As the sti ness increases, the linear
 load parameter needs to decrease in order to meet the Hopf instability criteria.
 The value found by Ulusoy (1994) for ~fl without any spring parameter was -0.22,
 whereas when the spring parameter with value of 0.10 is added, the value for ~fl has
 to be approximately -0.10 in order to meet the criterion for Hopf instability. Beyond
 a certain linear spring parameter value, one pair of eigenvalues does not cross the
 imaginary axis. This observation points to the important role played by the sti ness
 parameter in determining Hopf instabilities.
 28
?3.5 ?3 ?2.5 ?2 ?1.5 ?1 ?0.5 0 0.5?1.5
 ?1
 ?0.5
 0
 0.5
 1
 1.5
 Root Locus, kp=0.2
 Real
 Imaginar
 y
 (a)
 ?3 ?2.5 ?2 ?1.5 ?1 ?0.5 0 0.5?1.5
 ?1
 ?0.5
 0
 0.5
 1
 1.5
 Root Locus, kp=0.3
 Real
 Imaginar
 y
 (b)
 ?2.5 ?2 ?1.5 ?1 ?0.5 0?2
 ?1.5
 ?1
 ?0.5
 0
 0.5
 1
 1.5
 2
 Root Locus, kp=0.6
 Real
 Imaginar
 y
 (c)
 ?2.5 ?2 ?1.5 ?1 ?0.5 0?2
 ?1.5
 ?1
 ?0.5
 0
 0.5
 1
 1.5
 2
 Root Locus, kp=0.8
 Real
 Imaginar
 y
 (d)
 Figure 3.4: Root Locus Plot with varying ~fl and ~kp. (a) ~fl: -1.00 to 0.00 and
 ~kp=0.20. (b) ~fl: -1.00 to 0.00 and ~kp=0.30. (c) ~fl: -1.00 to 0.00 and ~kp=0.60. (d)
 ~fl: -1.00 to 0.00 and ~kp=0.80.
 3.2 E ect of Nonlinear Spring
 The stability of periodic solutions in the presence of nonlinear pressure loss and
 nonlinear damping load terms have been studied by Ulusoy (1994). He concluded
 that the nonlinear pressure loss and the nonlinear damping load terms both help
 produce limit-cycle motions. Here, the e orts are focused on studying the e ect
 of the nonlinear spring term added to the power piston. For a nonlinear spring
 function, a series expansion of the function can be written as a combination of
 29
linear and nonlinear spring components. In this study, a nonlinear spring with a
 linear spring element and a cubic nonlinear spring element is chosen. The force-
 displacement relationship is described as
 F (x) = kx+  kx3 (3.10)
 where  is the nonlinear spring coe cient. A positive value of  results in a harden-
 ing spring while a negative value of  results in a softening spring (Balachandran and
 Magrab, 2003). The use of softening spring can cause the system equilibrium posi-
 tion to be unstable; hence, a hardening spring is used in this study. The state-space
 model for the system is
 _x =
 2
 6
 6
 6
 6
 6
 6
 6
 6
 6
 6
 4
 0:00 1:00 0:00 0:00
  0:38  1:60  1:30 1:20
 0:00 0:00 0:00 1:00
 0:62 0:00  1:06 
~kp
 0:26
 ~fl
 0:26
 3
 7
 7
 7
 7
 7
 7
 7
 7
 7
 7
 5
 2
 6
 6
 6
 6
 6
 6
 6
 6
 6
 6
 4
 x1
 x2
 x3
 x4
 3
 7
 7
 7
 7
 7
 7
 7
 7
 7
 7
 5
 +
 2
 6
 6
 6
 6
 6
 6
 6
 6
 6
 6
 4
 0
  1HOT1 +  2HOT2
 0
  
~fc
 0:26
 x23x4 +  3HOT3  ~kcx
 3
 3
 3
 7
 7
 7
 7
 7
 7
 7
 7
 7
 7
 5
 (3.11)
 With the linear spring parameter value, ~kp, of 0.10 and the linear load term set
 to 0.00, the system is expected to be unstable when all the nonlinearities are ex-
 cluded. The expected outcome is veri ed in Figure 3.5, with initial displacements
 for displacer and power piston set at a value of 1.00 unit, while the velocities are
 set to 0.00 units; the initial conditions (ICs) are presented as (xd _xd xp _xp). The
 30
0 10 20 30 40 50 60 70?6000
 ?4000
 ?2000
 0
 2000
 4000
 6000
 Time (units)
 Amplitude (units
 )
 Displacer and Power Piston Movements
  
 
Displacer
 Piston
 Figure 3.5: Displacement response of the displacer and the power piston excluding
 nonlinear terms. The chosen parameters are ~fl = 0:00 and ~kp = 0:10.
 eigenvalues for the linearized system are given by
  1;2 = 0:14 1:01i
  3;4 =  0:94 0:66i
 where the eigenvalues  1;2 dominate; the associated oscillation frequency is approxi-
 mately 1 rad/s. Next, both the linear spring parameter, ~kp, and the linear damping
 parameter, ~fl are set to 0.10. The system response is still unstable with the linear
 damping term, as shown in Figure 3.6.
 Next, all of the nonlinear terms shown in equation (3.11) except the cubic
 spring term are held  xed at 0.00. The values of ~kp = 0:10, ~kc = 0:10, and ~fl = 0:00
 are chosen to determine the e ect of nonlinear spring added to the power piston. The
 31
0 10 20 30 40 50 60 70 80 90 100?3
 ?2
 ?1
 0
 1
 2
 3
 4 x 10
 7
 Time (units)
 Amplitude (units
 )
 Displacer and Power Piston Movements
  
 
Displacer
 Piston
 Figure 3.6: Displacement response of the displacer and the power piston excluding
 nonlinear terms. The chosen parameters are ~fl = 0:10 and ~kp = 0:10.
 initial condition to generate Figure 3.7(a) is (1.00 0.00 1.00 0.00), while the initial
 condition for Figure 3.7(b) is (7.00 0.00 7.00 0.00). The results are presented in
 terms of nondimensional values, represented by units. Comparing Figure 3.7(a) and
 Figure 3.7(b), the amplitude for the oscillating motions for di erent intial conditions
 are about the same: about 3.89 units for the power piston, and 2.82 units for the
 displacer. This is further analyzed through phase portraits as illustrated in Figure
 3.8, for the two di erent initial conditions. Figure 3.8(a) is the phase portrait for
 the displacer piston and Figure 3.8(b) is the phase portrait for the power piston.
 It can be seen that by adding a cubic spring term to power piston, the result is a
 limit-cycle motion when all of the other nonlinearities are not included.
 Next, the cubic spring system is compared to the system studied by Ulusoy
 (1994). The nonlinear pressure loss terms and nonlinear gas spring terms are in-
 cluded in the cubic spring system ( 1 = 1,  2 = 1, and  3 = 1 in equation (3.11)).
 32
0 10 20 30 40 50 60 70 80 90 100
 ?6
 ?4
 ?2
 0
 2
 4
 6
 Time (units)
 Amplitude (units
 )
 Displacer and Power Piston Movements
  
 
Displacer
 Piston
 (a)
 0 10 20 30 40 50 60 70 80 90 100
 ?6
 ?4
 ?2
 0
 2
 4
 6
 Time (units)
 Amplitude (units
 )
 Displacer and Power Piston Movements
  
 
Displacer
 Piston
 (b)
 Figure 3.7: Displacement responses for the displacer and the power piston due to
 additional nonlinear spring added to the power piston. The parameters are ~kp =
 0:10, ~kc = 0:10, and ~fl = 0:00. (a) ICs: (1.00 0.00 1.00 0.00) and (b) ICs: (7.00 0.00
 7.00 0.00).
 Reasonable values for ~kp and ~kc also must be chosen. By using the parameters given
 by Ulusoy (1994) and substituting them into the nondimensional form of spring
 33
?4 ?2 0 2 4 6 8?10
 ?8
 ?6
 ?4
 ?2
 0
 2
 4
 6
 8
 xd (units)
 x d
 ?vel (units
 )
 Phase Portrait
  
 
IC: [1 0 1 0]
 IC: [7 0 7 0]
 (a)
 ?8 ?6 ?4 ?2 0 2 4 6 8?15
 ?10
 ?5
 0
 5
 10
 15
 xp (units)
 x p
 ?vel (units
 )
 Phase Portrait
  
 
IC: [1 0 1 0]
 IC: [7 0 7 0]
 (b)
 Figure 3.8: Phase portraits with chosen initial conditions of (1.00 0.00 1.00 0.00)
 and (7.00 0.00 7.00 0.00). (a) Response of displacer and (b) Response of power
 piston.
 parameters, ~kp is approximately equal to 0:000001k and ~kc is approximately equal
 to 0:000001 k, where k is the spring constant and  is the nonlinear coe cient.
 Next, choosing a spring constant value of 100.00 N/m, and  of 100.00, the values
 34
of ~kp and ~kc are 0.0001 and 0.001, respectively. Finally, the values of ~fl =  0:02,
 and ~fc = 0:05 are used for comparing the Ulusoy?s model results to the results of
 the model studied in this chapter. Examining Figure 3.9, there does not seem to
 be much di erence due to the addition of a low sti fness spring because the mo-
 tions and amplitudes are close to each other in the two cases. The amplitude of
 the displacer with cubic spring and Ulusoy?s model parameters is about 1.14 units.
 Likewise, the amplitude of the power piston with cubic spring and Ulusoy?s model
 parameters is about 1.54 units. It can also be seen that the oscillation frequencies
 for the two models are close to each other since the motions for the two models
 almost overlap with each other.
 However, the oscillation frequency and the amplitude can be changed by in-
 creasing or decreasing the sti ness of the system. Increasing the sti ness should
 decrease the amplitudes of motions of the pistons while increasing the frequency.
 Hence, the values for ~kp and ~kc are chosen to be 0.01. The amplitudes of the dis-
 placer and power piston, along with the oscillation frequencies are compared to the
 case without any spring attached to the power piston. The obtained motions for
 the di erent cases are presented in Figure 3.10. It can be seen that the amplitudes
 without any spring added to the power piston have higher values. The amplitudes
 for the displacer and power piston without any spring parameter are 1.12 units and
 1.53 units, respectively. The amplitudes for the displacer and power piston with the
 added spring parameter are 1.07 units and 1.46 units, respectively, which are about
 5.00% less. While increasing the sti ness decreased the amplitude, the oscillation
 frequency of the system increased. This is shown in Figure 3.11, where the frequency
 35
0 10 20 30 40 50 60 70 80 90 100?1.5
 ?1
 ?0.5
 0
 0.5
 1
 1.5
 Time (units)
 Amplitude (units
 )
 Displacer Movement
  
 
With Cubic Spring
 Without Cubic Spring
 (a)
 0 10 20 30 40 50 60 70 80 90 100?2
 ?1.5
 ?1
 ?0.5
 0
 0.5
 1
 1.5
 2
 Time (units)
 Amplitude (units
 )
 Power Piston Movement
  
 
With Cubic Spring
 Without Cubic Spring
 (b)
 Figure 3.9: Displacement response comparison of cases with low sti ness nonlinear
 spring and without nonlinear spring. (a) Response of displacer and (b) Response of
 power piston.
 spectrum for the two di erent cases are plotted. The oscillation frequency of the
 system without any spring parameter added to the power piston is about 24.20 Hz
 while the oscillation frequency for the added spring parameter case is 25.50 Hz; there
 36
10 20 30 40 50 60 70 80 90
 ?1.5
 ?1
 ?0.5
 0
 0.5
 1
 Time (units)
 Amplitude (units
 )
 Displacer Movement
  
 
Without Cubic Spring
 With Cubic Spring
 (a)
 10 20 30 40 50 60 70 80 90 100
 ?2
 ?1.5
 ?1
 ?0.5
 0
 0.5
 1
 1.5
 Time (units)
 Amplitude (units
 )
 Power Piston Movement
  
 
Without Cubic Spring
 With Cubic Spring
 (b)
 Figure 3.10: Displacement response comparison of cases with high sti ness nonlinear
 spring and without nonlinear spring. (a) Response of displacer and (b) Response of
 power piston.
 is an increase of about 5.00%.
 Finally in Figure 3.12, the amplitude responses of the displacer and power
 pistons as the linear, ~kp, and nonlinear, ~kc, spring parameters are varied are shown.
 37
23.5 24 24.5 25 25.5 26 26.5 27 27.5
 0.5
 1
 1.5
 2
 2.5
 x 104
 Frequency (Hz)
 m/H
 z
  
 
Without Cubic Spring
 With Cubic Spring
 Figure 3.11: Frequency spectrum comparison for cases with and without nonlinear
 spring.
 The linear load parameter, fl, is held constant at -0.02. It can be seen that when the
 linear sti ness of the system is 0.00, the system is unstable. However, when the linear
 sti ness is introduced to the system, the response amplitude decreases. Increasing
 the linear sti ness of the system quickly decreases the response amplitude of the
 system to zero. However, with only a large cubic sti ness, the amplitude response
 variation of the system does not go to zero quickly.
 38
(a)
 (b)
 Figure 3.12: Amplitude response variation with respect to linear sti ness parameter
 ~kp, and nonlinear sti ness parameter ~kc, with linear load parameter, fl, set to -0.02.
 (a) Response of displacer and (b) Response of power piston.
 39
Chapter 4
 Governing Equations for System with Magnetic Springs
 Introducing nonlinearities causes the FPSE to exhibit limit-cycle motions. Pre-
 viously, it was shown that nonlinear springs, especially hardening springs, also cause
 the system to exhibit periodic motions. Hence, magnetic springs are introduced into
 the FPSE, which display hardening spring behavior in one direction and softening
 spring behavior in another direction. The schematic for the magnetic spring FPSE
 is shown in Figure 4.1. Magnetic springs can be more robust than conventional
 springs or elastic bands because the magnetic springs would not be in contact with
 one another. Hence, using magnetic springs can increase the longevity of a FPSE
 and may make the FPSE more reliable. Also, the property of elastic bands depend
 Figure 4.1: Schematic of magnetic spring FPSE.
 40
on temperature, displacement, and time, which may give very inconsistent results.
 Unlike gas springs, the tolerances to use magnetic springs can be larger since they
 do not require the FPSE to be perfectly sealed to prevent gas leakage. Finally,
 the usage of magnetic springs can also reduce the operational noise in the system.
 Gravity must be introduced to the system in order to have restoring force for proper
 analysis. The magnetic spring force can be described as
 Fmagneticspring = Mg  Fmagnet (4.1)
 where M is the mass of the top magnet, g is gravity, and Fmagnet is the repulsive
 magnetic force. Adding gravity to the system is critical, since as the magnets are
 repulsed away from each other, the weight of the top magnet acts as the restoring
 force. A generic force curve of magnetic spring, derived in the next section, is shown
 in Figure 4.2. The black square in the  gure represents a magnet; it can be seen
 that the magnetic spring behaves as a hardening spring when the upper magnet
 comes in close proximity with the lower magnet and as a softening spring when it
 goes further away from the other magnet. The equilibrium point of the magnetic
 spring can be varied by changing the mass of the top magnet.
 The equations of motions with the magnetic spring can be written as
  md zd = Ad P + Fmagnet1   mdg (4.2)
  mp zp = ApPm
  
1 
1
 F
  
+ Fload _zp + Fmagnet2   mpg (4.3)
 41
?0.05 ?0.04 ?0.03 ?0.02 ?0.01 0 0.01 0.02 0.03 0.04 0.05
 ?500
 ?400
 ?300
 ?200
 ?100
 0
 100
 200
 300
 Distance (m)
 Force (N
 )
 Spring Force?Displacement Curve
  
 
Magnetic
 Linear
 Hardening
 Softening
 Figure 4.2: Magnetic spring force curve compared to various types of spring forces.
 where
 zd = xd + e1
 zp = xd + e2
  md = md +M1
  mp = mp +M2
 Fload = ( fl  fcz
 2
 p)
 Furthermore, a and b are the equilibrium distances between the two repulsive mag-
 nets, M1 and M2 are the masses of the magnets, fl is the linear damping parameter,
 and fc is the nonlinear damping term. The equilibrium distances can be found by
 42
setting all the velocity and acceleration terms to zero; that is
  mdg = Fmagnet1 (4.4)
  mpg = ApPm
  
1 
1
 F
  
+ Fmagnet2 (4.5)
 4.1 Nondimensional Equations
 In this section, the variables for the governing equations are nondimensional-
 ized. The di erent quantities are as follows:
 Masses:
 ~md =
  mdld
 AdPm
 (4.6)
 ~mp =
  mpld
 AdPm
 (4.7)
 Displacements:
 ~zd =
 zd
 ld
 (4.8)
 ~zp =
 zp
 ld
 (4.9)
 where ld is the length of the displacer rod.
 Area:
 ~Ad =
 Ad
 Ad
 (4.10)
 ~Ap =
 Ap
 Ad
 (4.11)
 43
Pressure:
  ~P =
  P
 Pm
 (4.12)
 ~Pm =
 Pm
 Pm
 (4.13)
 Forces:
 ~Fmagnet1 =
 Fmagnet1
 AdPm
 (4.14)
 ~Fmagnet2 =
 Fmagnet2
 AdPm
 (4.15)
 ~Fload =
 Fload
 AdPm
 (4.16)
 ~md~g =
  mdg
 AdPm
 (4.17)
 ~mp~g =
  mpg
 AdPm
 (4.18)
 Other quantities:
 ~a =
 a
 ld
 (4.19)
 ~b =
 b
 ld
 (4.20)
 ~F = 1 ~a~xp + ~b~xd (4.21)
 The nondimensional equations of motions are of the form
 ~md~ zd = ~Ad ~P + ~Fmagnet1  ~md~g (4.22)
 ~mp~ zp = ~Ap ~Pm
  
1 
1
 ~F
  
+ ~Fload~_zp + ~Fmagnet2  ~mp~g (4.23)
 44
4.2 Repulsive Magnetic Forces
 In order to understand how the repulsive magnetic forces behave, the principle
 of superposition is applied, which states that the interaction between any two charges
 is completely una ected by the presence of others (Gri ths, 1999). Hence, in order
 to  nd a force on a particular charge, Q, among several charges, qi, one must compute
 force F1 due to q1 while ignoring all other charges, then, compute force F2 due to
 q2 alone, and so forth. Finally, the total force felt by charge Q is determined as
 the summation of all forces. The force between any two point charges are given by
 Coulomb?s law, which reads as
 F =
 1
 4  0
 qQ
 r2
 r^ (4.24)
 where  0 is the permittivity of free space, q is the source charges, Q is the test charge,
 r is the distance between the two charges, and r^ is the associated unit vector. The
 force is repulsive if q and Q have the same sign, and attractive if they have opposite
 signs. Finally, for several point charges q1, q2,...,qn at distances r1, r2,...,rn from Q,
 the total force on Q is given by (Gri ths, 1999)
 F = F1 + F2 + ::: =
 1
 4  0
  
q1Q
 r21
 r^1 +
 q2Q
 r22
 r^2 + :::
  
(4.25)
 F =
 Q
 4  0
  
q1
 r21
 r^1 +
 q2
 r22
 r^2 + :::
  
(4.26)
 45
F = QE (4.27)
 where E is the electric  eld of the charges and is de ned as:
 E = C
  
nX
 i=1
 qi
 r2i
 r^i
 !
 (4.28)
 and
 C =
 1
 4  0
 (4.29)
 Next, each magnet used here is modeled as magnetic dipole, or a set of two charges
 as shown in Figure 4.3. As such, the charge is de ned as
 Z
 S
  0MdA =  qi (4.30)
 where S is a closed surface area, qi is the charge enclose within that area, and M
 is the magnetization. The charge is de ned as negative because the magnetic  eld
 moves from positive to negative charge (Gri ths, 1999). Integrating over a circular
 surface area leads to
 qi =   0M R
 2 (4.31)
 Figure 4.3: Magnets modeled as charges.
 46
which is negative for a negative charge source, but positive for a positive charge
 source.
 Assuming constant radius, the forces between each charge from one magnet to
 another are found and then summed as follows
 Fq1;q3 = C
 q1q3
 (x+ L)2
 (4.32)
 Fq1;q4 = C
 q1q4
 (x+ 2L)2
 (4.33)
 Fq2;q3 = C
 q2q3
 (x)2
 (4.34)
 Fq2;q4 = C
 q2q4
 (x+ L)2
 (4.35)
 After using equation (4.31), an approximation for the repulsive magnetic force can
 be found as
 X
 F =
   0
 4
 M2R4
  
1
 x2
 +
 1
 (x+ 2L)2
  
2
 (x+ L)2
  
(4.36)
 The force-displacement curve of the repulsive magnetic spring is shown in Figure
 4.4. The dotted points are the representative force given by the manufacturer KJ-
 Magnetics, while the other line is the model prediction curve. Both curves follow a
 similar pro le; also, the range of motion for the numerical and experimental studies
 is greater than 1.00 cm. Large values for the constant will make the magnetic spring
 sti er, while small values will make the magnetic spring softer.
 Equation (4.36) is compared to the magnetic force derived by Vokoun, Beleg-
 gia, Heller, and Sittner (2009). They derive the force between sets of magnets from
 47
the total magnetostatic interaction energy E of the system, which is given by:
 Fmag =  grad(E) (4.37)
 The assumptions made are that the magnets are cylindrical and that they are made
 of the same material characterized by saturation magnetization M . The cylindrical
 magnets are of equal radius R and are magnetized uniformly along the cylinder axis
 of symmetry. The force between two cylindrical permanent magnets with a common
 axis is derived to be (Vokoun et al., 2009)
 Fmagnet =
   0
 4
 M2R4
  
1
 x2
 +
 1
 (x+ 2t)2
  
2
 (x+ t)2
  
(4.38)
 where M is the magnetization, R is the radius of the magnet, t is the height of the
 magnet, and  0 is the permeability of vacuum. Comparing equations (4.36) and
 0 0.005 0.01 0.015 0.02 0.025 0.030
 10
 20
 30
 40
 50
 60
 70
 80
 90
 100
 Distance (m)
 Force (N
 )
  
 
Representative fit from the manufacturer
 Model prediction
 Figure 4.4: Repulsive force between two similar magnets.
 48
(4.38), it can be seen that the equations are similar.
 By using Equations (2.43), (4.22), and (4.23), the nondimensionalized equa-
 tions of motion for the system with the magnetic spring can be written as
 ~md~ xd =
 ~Ad
 Pm
  
1
 2
  (ft + kh)
 2Ad _zd  Ap _zp
 Ak
 j
 2Ad _zd  Ap _zp
 Ak
 j+
   0M2dR
 4
 d
 4AdPm
  
1
 z2d
 +
 1
 (zd + 2t)2
  
2
 (zd + t)2
  
 ~md~g (4.39)
 ~mp~ xp = ~Ap
  
1 
1
 F
  
+ ~Fload~_zp+
   0M2pR
 4
 p
 4AdPm
  
1
 z2p
 +
 1
 (zp + 2t)2
  
2
 (zp + t)2
  
 ~mp~g
 (4.40)
 Here Ak is the area of the annulus between the displacer and the wall in Figure 4.1
 through which the working  uid is displaced.
 4.3 Results for System Response
 In this section, various  ndings for the modi ed system with modeled pressure
 loss is presented. The parameters used for the modeling are presented in Appendix
 B. First, the magnetic spring case without any damping is shown. Then, the head
 loss within the annulus and damping is introduced.
 The movements of the displacer and power piston without any damping can
 be seen in Figure 4.5. The initial condition given to the system is (0.01 0.0 0.0 0.0).
 Although damping is not present, the system shows oscillatory motions. This is
 further veri ed through the phase portrait diagram shown in Figure 4.6. The power
 49
1 2 3 4 5 6 7 8 9 10
 ?0.06
 ?0.05
 ?0.04
 ?0.03
 ?0.02
 ?0.01
 0
 0.01
 0.02
 0.03
 Time (s)
 Displacement (m
 )
  
 
Displacer
 Power Piston
 Figure 4.5: Displacement response of magnetic spring FPSE without any damping
 or head loss.
 piston does not follow the movements of the displacer, which shows that there is a
 weak coupling between the two pistons.
 Next, a low level of damping and head loss is introduced into the power piston.
 The linear damping term fl is chosen to be 0.01, the cubic nonlinear term fc is set
 at 0.02, and the head loss kh is set at 0.05. The initial condition is held constant at
 ?0.01 ?0.005 0 0.005 0.01 0.015?0.2
 ?0.15
 ?0.1
 ?0.05
 0
 0.05
 0.1
 0.15
 0.2
 xd (m)
 x
 d
 ?vel (m/s
 )
 (a)
 ?0.08 ?0.06 ?0.04 ?0.02 0 0.02 0.04?0.8
 ?0.6
 ?0.4
 ?0.2
 0
 0.2
 0.4
 0.6
 0.8
 xp (m)
 x
 p
 ?vel (m/s
 )
 (b)
 Figure 4.6: Phase portraits for the case with no damping and with ICs: (0.01 0.0
 0.0 0.0). (a) Displacer and (b) Power Piston.
 50
0 10 20 30 40 50 60 70 80 90 100?0.06
 ?0.04
 ?0.02
 0
 0.02
 0.04
 Time (s)
 Displacement (m
 )
  
 
Displacer
 Power Piston
 Figure 4.7: Displacement response of magnetic spring FPSE with low damping and
 head loss.
 their previous values. The amplitudes for the displacer and power piston are shown
 in Figure 4.7. As expected, the amplitudes decrease with damping. Once again, the
 phase portrait is shown to verify oscillatory motions in Figure 4.8.
 Finally, the damping parameter is increased to verify that the amplitudes die
 out when the system is highly damped. The linear damping term fl is chosen to
 ?0.01 ?0.005 0 0.005 0.01 0.015?0.2
 ?0.15
 ?0.1
 ?0.05
 0
 0.05
 0.1
 0.15
 0.2
 xd (m)
 x
 d
 ?vel (m/s
 )
 (a)
 ?0.06 ?0.05 ?0.04 ?0.03 ?0.02 ?0.01 0 0.01 0.02 0.03 0.04?0.8
 ?0.6
 ?0.4
 ?0.2
 0
 0.2
 0.4
 0.6
 0.8
 xp (m)
 x
 p
 ?vel (m/s
 )
 (b)
 Figure 4.8: Phase portraits for the case with low damping and head loss and with
 ICs: (0.01 0.0 0.0 0.0). (a) Displacer and (b) Power Piston.
 51
0 20 40 60 80 100 120 140 160 180 200?0.03
 ?0.02
 ?0.01
 0
 0.01
 0.02
 Time (s)
 Displacement (m
 )
  
 
Displacer
 Power Piston
 Figure 4.9: Displacement response of magnetic spring FPSE with high damping and
 head loss.
 be 0.50, the cubic nonlinear term fc is set at 0.50, and the head loss kh is set at
 5.00. From Figure 4.9, it can be seen that the amplitudes for both the displacer
 and the power piston die out with large damping. However, when the temperature
 82 84 86 88 90 92 94 96 98
 ?0.1
 ?0.05
 0
 0.05
 0.1
 Time (s)
 Displacement (m
 )
  
 
Displacer
 Power Piston
 Figure 4.10: Displacement response of magnetic spring FPSE with high damping
 and head loss with temperature di erence of 500:00 Celcius.
 52
di erence between the expansion space and compression space is increased to 500:00 
Celsius, the system response shows oscillatory motion as shown in Figure 4.10. The
 numerical studies of this chapter indicate that with a magnetic spring, a FPSE may
 not show oscillatory motions due to high levels of damping. However, this may be
 overcome with an appropriate temperature gradient. In addition, further work is
 needed to realize limit-cycle motions in a FPSE with magnetic springs.
 53
Chapter 5
 Experimental Studies
 In this chapter, experimental studies conducted in an attempt to realize a
 FPSE system are presented. Various con gurations of the FPSE are designed for
 the possibilities of having oscillatory motions. The kinematic Stirling engine shown
 in Figure 5.1 is modi ed and turned into a FPSE. The  ywheel connecting the power
 piston and the displacer are disconnected and replaced with magnets to levitate the
 power piston and the displacer. The power piston is fabricated from graphite, and
 the compression volume around the power piston is made of borosilicate glass in
 order to reduce friction.
 The schematic of the experimental setup is shown in Figure 5.2. The power
 source is connected to the temperature controller, which controls the thermoelectric
 Figure 5.1: Kinematic Stirling engine to be modi ed into FPSE.
 54
Figure 5.2: Schematic of the experimental setup.
 cooler (TEC) and the thermistor. Both, the TEC and the thermistor are attached to
 the bottom of the FPSE. The TEC can be controlled to vary the temperature of the
 expansion space of FPSE while the thermistor is used to measure the temperature.
 At the top, a  ber optic sensor is used to measure the movement of the power piston.
 5.1 Design and Fabrication
 The  rst iteration of the experimental setup is shown in Figure 5.3. The rod
 connecting the power piston to the magnet, along with the magnet holders were
 fabricated by this thesis author. The top magnets are attached to the magnet
 holders via a press  t. The holders are designed to be adjustable so that it can slide
 up and down. This way, the equilibrium point can be adjusted by adding weights
 55
such that the displacer and power piston sit at the middle of the expansion and
 compression space, respectively. In this setup, the magnetic interference between
 the two magnetic springs are found to be high because of their close proximity; thus,
 introducing high damping into the system due to eddy currents. The repulsion and
 attraction between the top magnet for the displacer and the magnets for the power
 piston cause the displacer rod to bend slightly either towards or away from the the
 power piston. This increases the friction between the displacer rod and its casing.
 Similarly, interactions between the displacer and the power piston magnets increase
 the friction between the power piston and its casing.
 In order to reduce the magnetic interference and decrease piston rod related
 damping, the magnetic spring for the power piston is elevated. The modi ed
 schematic is shown in Figure 5.4. A tube that encased the power piston was fab-
 ricated, and the rod connecting the power piston was made longer. However, the
 dimensions and the placement of the power piston were kept the same as that of the
 Figure 5.3: Modi ed FPSE.
 56
previous con guration. Increasing the height of the power piston rod had adverse
 a ects. The weight of the magnet and the holder for the power piston caused the
 power piston rod to tilt. This increased the normal force between the power piston
 and the casing, which in turn increased the friction. Another problem that was also
 noticed was that the displacer had very little room for movement. Hence, when the
 displacer moved, it tended to stick to the top or the bottom of the casing. When
 the displacer sticks to the top, it closes the path for the air to get into compres-
 sion space; hence, the air is prevented from being compressed. Likewise, when the
 displacer sticks to the bottom, the air is prevented from being expanded.
 Finally, the expansion space was doubled so that the displacer could have
 more room for oscillations. This con guration is shown in Figure 5.5. Again, this
 con guration did not produce oscillatory motions. One of the reasons might be that
 the compression space is small compared to that of the expansion space. Thus, when
 Figure 5.4: Increased height for the power piston.
 57
Figure 5.5: Increased volume for the expansion space.
 the compression space collapses, the decrease in volume does not lead to a su cient
 pressure increase in the compression area. Hence, the pressure di erence across the
 displacer does not increase very much.
 5.2 Preliminary Experimental Results
 The di erent con gurations were heated to about 315:00 Celsius while the
 cold end of the engine was at room temperature of about 25:00 Celsius. The pa-
 rameters for the experimental setup are given in Appendix B. The response obtained
 for the kinematically linked Stirling engine is shown in Figure 5.6. For this system,
 oscillatory motions are observed. The frequency of the system is about 6.20 Hz.
 For the  rst modi ed system depicted in Figure 5.3, the displacer and power
 piston motions decay right away. The motion of the power piston was captured
 using a Philtec optical sensor, and the observed response is shown in Figure 5.7.
 58
0 1 2 3 4 5 6 7 8 9 10?4
 ?3
 ?2
 ?1
 0
 1
 2
 3
 4 x 10
 ?3
 Time (s)
 Distance (m
 )
 Figure 5.6: Power piston displacements for a kinematically linked Stirling engine.
 Due to high damping, the decaying motion is as expected. By using the data along
 with the logarithmic decrement method, the damping ratio is approximated to be
 0.31, and the undamped frequency is 1.59 Hz.
 For the second modi ed system depicted in Figure 5.8, the power piston os-
 0.2 0.4 0.6 0.8 1 1.2 1.4
 ?10
 ?5
 0
 5
 x 10?3
 Time (s)
 x p
  (m
 )
 Figure 5.7: Power piston displacements for  rst FPSE arrangement.
 59
0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8
 ?0.01
 ?0.005
 0
 0.005
 0.01
 Time (s)
 x p
  (m
 )
 Figure 5.8: Power piston displacements for second FPSE arrangement.
 cillates very little and decays again. The oscillation of the power piston is shown in
 Figure 5.8. Again, the decaying oscillations are due to the high damping inherent
 in the system due to the design. The damping ratio approximated for the con g-
 uration is 0.21, and the undamped frequncy is 5.60 Hz. The last con guration did
 not provide any positive results either as the decaying motions were very similar to
 the previous con guration. The numerical studies of the previous chapter support
 the observations made in this chapter; that is, a high damping level can inhibit os-
 cillatory motions. Future experiments can be conducted with a higher temperature
 gradient to see if that would help in realizing oscillatory motions.
 60
Chapter 6
 Summary and Concluding Remarks
 A FPSE has various advantages over kinematically linked Stirling engines; as
 such, FPSEs would be better for practical purposes. A kinematically linked Stirling
 engine relies on the linkage between the displacer and power piston for its motions.
 A FPSE rely solely on the pressure di erence created by the temperature gradient
 and volume variations. Various parameters within the FPSE are critical for its func-
 tionality. Hence, parametric studies have been conducted to study which parameters
 can help the system exhibit oscillatory motions. Additionally, a study of a nonlinear
 system with a spring attached to the power piston is conducted. The considered
 beta engine exhibits Hopf instabilities with a hardening spring on the power piston
 side and the introduction of nonlinear spring leads to limit-cycle motions. Next,
 magnetic springs are introduced into the FPSE system, and the motions of the dis-
 placer and the power piston are studied. The system with magnetic springs shows
 oscillatory motions. Finally, experimental studies are conducted in an attempt to
 show oscillatory responses for a FPSE. The signi cance of this thesis work is the
 understanding developed on Hopf instabilities in FPSEs. Also, the experimental
 studies help understand di culties in realizing a FPSE from a kinematic Stirling
 engine. This thesis study complements the analytical-numerical work conducted in
 the group by Choudhary (2009).
 61
6.1 Parametric Studies
 The model studied by Ulusoy (1994) is modi ed to include spring term (~kp)
 for the power piston. With the introduction of the spring parameter, the Jacobian
 matrix of the system is investigated. By using the characteristic equation found from
 the Jacobian, a parametric study is conducted by varying the spring parameter, ~kp,
 and the damping parameter, ~fl. It is shown that a Hopf instability can be attained
 by varying the damping parameter, while the spring parameter is zero. Additionally,
 one can have a Hopf instability in the absence of the system damping due to the
 spring parameter as well. The in uences of the damping parameter and the spring
 parameter are established. With increasing sti ness in the system, it is found that
 the system damping needs to be reduced for oscillatory motions to occur.
 Next, a nonlinear hardening spring is introduced into the power piston. The
 system response is unstable for a case with linear damping. However, the addition
 of the nonlinear spring leads to oscillatory motions. Introducing the spring to the
 system a ected the amplitude and the oscillatory frequency of the system. Compar-
 ing the motions of FPSE without nonlinear spring added to the power piston to the
 system with a low sti ness nonlinear spring added to the power piston showed that
 the amplitude between the two does di er by much. By increasing the sti ness of
 the spring, the response amplitude decreases but the oscillation frequency increases.
 The advantage realized by adding a nonlinear spring to the power piston is to realize
 limit-cycle motion in a FPSE without a cubic damper, while also being able to vary
 the oscillatory frequency of the system.
 62
6.2 System with Magnetic Springs
 Once it was found that nonlinear spring helps attain oscillatory motions, mag-
 netic springs were introduced into the system. The advantages of the inclusion of
 magnetic springs in a FPSE are less friction leading to higher durability and less
 operational noise. The repulsive forces between two magnets are found by modeling
 the charges as dipole moments along with Coulomb?s law. The numerical results
 obtained for the magnetic spring model showed that the displacer and the power pis-
 ton are weakly linked. The results showed that oscillatory motions are possible with
 magnetic springs. However, with very high damping, the motions die out. In order
 to overcome large damping, a high temperature di erence between the expansion
 and the compression space is necessary.
 6.3 Experimental Studies
 Various con gurations for FPSE were designed and fabricated in an attempt
 to realize oscillatory motions. The  rst modi ed con guration shown in Figure 5.7
 had very high damping due to the proximity of the magnets used in the displacer and
 power pistons. Hence, the motions decay quickly. This con guration was modi ed
 as shown in Figure 5.8. The new design created high friction between the power
 piston and the casing along with the displacer sticking to the hot end and cold end of
 the FPSE. Therefore, this con guration did not exhibit oscillatory motions either.
 Finally, the expansion volume was expanded to overcome the sticking problem.
 This modi cation did not result in oscillatory motions either. It is suspected that
 63
the damping levels in the experimental system may not be appropriate to realize
 oscillatory motions; a large temperature di erence between the expansion and the
 compression space would be needed to realize oscillatory motions. The experimental
 studies also helped realize the di culties in fabricating a FPSE from a kinematically
 linked Stirling engine.
 6.4 Future Work
 The present thesis work supports the promise of realizing low power FPSEs.
 However, further analytical and numerical studies are needed to understand the
 nature of Hopf instabilities. NASA has previously considered Stirling engines for
 space applications. For these applications, the current system with magnetic springs,
 which relies on gravity will not work. In order to overcome this, the design of the
 magnetic springs can be changed to have a magnet repel in between two di erent
 magnets, which also makes the magnetic spring sti from both directions. Addi-
 tionally, the experiments were not successful in realizing oscillatory motions. The
 experimental setup should be redesigned to reduce the a ect of damping in order
 to realize oscillatory motions. Also, various other designs for FPSEs need to be
 considered and appropriate parametric studies need to be conducted.
 64
Appendix A
 Equations of Motion and Parameters from Ulusoy (1994)
 Equations of Motion:
 The schematic of the FPSE that Ulusoy (1994) studied is shown in Figure A.1.
 The equations of motions for this system are given by
 md xd = Ad(Pc  Pe) Ar(Pc  Ps) (A.1)
 mp xp = Ap(Pb  Pc) + Fload _xp (A.2)
 Figure A.1: Schematic of FPSE from Ulusoy (1994).
 65
The equations of motions resulting after Schimdt analysis are as follows
  xd =
 Ad
 md
  P  
ArPm
 md
  
F  
 
1
 1 + asxd
    
(A.3)
  xp =
 ApPm
 mp
 (1 F ) +
 Fload
 mp
 (A.4)
 where
 F =
 1
 1 apxp + adxd
 (A.5)
 Fload =
  
fL  fcx
 2
 p) _xp
  
(A.6)
 Next the parameters are nondimensionalized. These parameters are given by
 ~xp =
 xd
 lp=2
 (A.7)
 ~xp =
 xp
 lp=2
 (A.8)
 ~P =
 P
 Pm
 =
 1
 1 ~ap~xp + ~ad~xd
 (A.9)
 ~ap = ap
 lp
 2
 (A.10)
 ~ad = ad
 lp
 2
 (A.11)
 ~Pc = Pc=Pm = ~P (A.12)
 4 ~P =
 4P
 Pm
 (A.13)
 ~Pe = Pc=Pm = ~P  4 ~P (A.14)
 ~Ps = Ps=Pm =
  
1
 1 + ~as~xd
   
(A.15)
 66
~as = as
 lp
 2
 (A.16)
 ~md =
 !2 lp2md
 PmAp
 (A.17)
 ~mp =
 !2 lp2mp
 PmAp
 (A.18)
 ~Ad = Ad=Ap (A.19)
 ~Ar = Ar=Ap (A.20)
 ~Ap = Ap=Ap = 1 (A.21)
 ~Fload =
  
~fL  ~fcx
 2
 3
  
x4 (A.22)
 ~fL =
 fL!
 lp
 2
 PmAp
 (A.23)
 ~fc =
 fc!
  lp
 2
  3
 PmAp
 (A.24)
 ~kp =
 kp
 lp
 2
 PmAp
 (A.25)
 where lp is the maximum stroke of the power piston. Hence, the nondimensional
 equations become:
  ~xd =
 ~Ad
 ~md
  ~P  
~Ar
 ~md
  
~F  
 
1
 1 + ~as~xd
    
(A.26)
  ~xp =
 1
 ~mp
  
1 ~F
  
+
 ~Fload
 ~mp
 (A.27)
 Ulusoy?s Parameters:
 The parameters used in his work are provided next.
 Temperatures (K):
 Th = 900
 67
Tk = 300
 Masses (kg):
 Power Piston Mass (mp) = 6.2
 Displacer Mass (md) = 0.426
 Volumes (m3):
 Displacer Gas Spring Mean Volume (Vso) = 31:8  (10 6)
 Compression Space Mean Volume (Vco)= 103:6  (10 6)
 Expansion Space Mean Volume (Veo) = 63:6  (10 6)
 Heater Volume (Vh)= 27:309  (10 6)
 Cooler Volume (Vk)= 20:422  (10 6)
 Regenerator Volume (Vr)= 186:534  (10 6)
 Design Frequency (Hz)= 30
 Piston length stroke (cm) =4.20
 Mean Pressure (Pa) = 7,000,000
 68
Appendix B
 Experimental Parameters
 The parameters corresponding to the FPSEs studied in the Vibrations Labo-
 ratory are provided below.
 Masses (kg)
 md = 0.001
 mp = 0.001
 mda = 0.1
 mpa = 0.1
 mdt = md+mda
 mpt = mp+mpa
 Length (m)
 ld= 21/1000
 Area (m2)
 Ad= pi  (140=2000)2
 Ap= pi  (35=2000)2
 Ak= pi  (150=2000)2
 Ah = pi  (150=2000)2
 69
Temperatures (K)
 Th = 315 + 273
 Tk = 25 + 273
 Volumes (m3)
 Vcm = pi  (35=2000)2  (20=1000) + pi  (152=2000)2  (10=1000)
 Vem = pi  (152=2000)2  (13=1000)
 Pressures (Pa)
 Pm = 101500
 Magnets
 mu =4  pi  (10 7)
 Bo = 1.3200 (T)
 R = 0.0254 (m)
 t = 0.003175 (m)
 M = (2  Bo=mu)
 Other Parameters
 g=9.81 (m=s2)
 rho=1
 Do= .152 (m)
 70
Di= .139 (m)
 L=10=1000 (m)
 71
Appendix C
 Matlab Codes
 Representative programs used in this thesis study are included in this Ap-
 pendix.
 72
% Parametric study 
 
close all 
clear all 
clc 
  
format short 
  
b=zeros(100,5); 
t=1; 
n=0; 
for kp=0 
    n=n+1; 
    figure(n); 
    tic 
    for i=-1:.01:0 
        A=[0 1 0 0; -.389 -1.598 -1.302 1.182; 0 0 0 1; 0.614 0 -1.056-
 (kp/0.257) i/0.257]; 
%       A=[0 1 0 0; -.389 -1.598 -1.302 1.182; 0 0 0 1; 0.614 0 -1.0171-
 (kp/0.257) i/0.257]; 
        c=eig(A); 
%         figure(n) 
        scatter(real(c),imag(c),'.') 
        hold on 
        b(t,1)=i; 
        b(t,2)=c(1,1); 
        b(t,3)=c(2,1); 
        b(t,4)=c(3,1); 
        b(t,5)=c(4,1); 
        t=1+t; 
    end 
    ti=title(['Root Locus, kp=' num2str(kp)]); 
    x=xlabel('Real'); y=ylabel('Imaginary'); 
    set(ti,'Fontsize',12) 
    set(x,'Fontsize',12) 
    set(y,'Fontsize',12) 
    axis on 
    grid on 
    toc 
end 
 
 
%%Solving for magnetic equilibriums  
 
u=4*Pi*10^ -7; 
R=.0254; 
t=0.003174; 
B=1.32;  
M=(2*B/u);  
Ap=Pi*(45/1000)^2;  
Ad=Pi*(140/1000)^2;  
Vcm=Pi*(45/1000)^2*(8/1000);  
Vem=Pi*(145/1000)^2*(10/1000);  
Pm= 500; 
Tk=273+0;  
Th=500+100;  
S=Vem/Th+Vcm/Tk;  
a=Ap/Tk(1/S);  
b=(Ad/Tk -Ad/Th)(1/S);  
m1=.1001; 
 
Solve[m1*9.81 -(Pi*M^2*u*R^4)/4*(1/(z1)^2+1/(z1+2t)^2 -2/(z1+t)^2) 0,z1] 
 
 
u=4*Pi*10^ -7; 
R=.0254; 
t=0.003174; 
B=1.32;  
M=(2*B/u);  
Ap=Pi*(45/1000)^2;  
Ad=Pi*(140/1000)^2;  
Vcm=Pi*(45/1000)^2*(8/1000);  
Vem=Pi*(145/1000)^2*(10/1000);  
Pm= 500; 
Tk=273+0;  
Th=500+100;  
S=Vem/Th+Vcm/Tk;  
a=(Ap/Tk)(1/S);  
b=(Ad/Tk -Ad/Th)(1/S);  
m1=.1001; 
F=1-a*(z2)+b*(0.09964981303487692`);  
Solve[m1*9.81 -(Pi*M^2*u*R^4)/4*(1/(z2)^2+1/(z2+2t)^2 -2/(z2+t)^2) -Ap*Pm (1 -
 1/F) 0,z2] 
 
function xp=Exp1(t,x)  
  
%%% Masses %%%  
  
md = 0.001; % mass of displacer 
mp = 0.001; % mass of piston 
  
mda = .1; % mass added to displacer 
mpa = .1; % mass added to piston 
  
mdt = md+mda; % total displacer mass 
mpt = mp+mpa; % total piston mass 
  
%%% Length %%% 
  
ld= 21/1000; % length of displacer rod 
  
%%% Area %%% 
  
Ad= pi*(140/1000)^2; % Area of displacer 
Ap= pi*( 15/1000)^2; % Area of piston 
Ak= pi*(150/1000)^2; %Area of cooler 
Ah = pi*(150 /1000)^2; %Area of heater 
  
%%% Temperatures %%%  
  
Th = 500+273; % Heater temperature  
Tk = 0+273; % Cooler temperature  
  
%%% Volumes %%%  
  
Vcm = pi*(45/1000)^2*(10/1000); %Mean compression volume  
Vem = pi*(145/1000)^2*(10/1000); %Mean expansion volume  
  
S = Vem/Th+ Vcm/Tk; % Introducing new term  
a = Ap/Tk * (1/S); % Placeholder 
b = (Ad/Tk - Ad/Th) * (1/S); % Placeholder 
  
%%% Pressures %%% 
  
Pm = 101500; %Mean pressure;  
  
  
%%% Magnets %%%  
  
mu = 4*pi*(10^ -7);  
Bo = 1.3200; % Magnetic flux density value  
R = 0.0254; % Radius of magnet 
h = 0.003175; % Thickness of magnet  
M = (2*Bo/mu); % Magnetization of magnet  
mfield =(pi*mu*M^2*R^4)/4;  
  
  
%%% Other %%%  
  
g=9.81; %gravity  
rho=1; 
fl=0.05; %linear damping 
fc=0.05; %cubic damping 
kh=.05; %head loss 
  
%%% EOMS %%%  
  
xp=zeros(2,1); 
  
Cf= 0.0457; % Fanning friction coeff; 
n1= -0.2; %recommended by Martini;  
Do= 152.4; %Outer diameter;  
Di= 139.7; %Inner diameter; 
dh = 4*0.25*pi*(Do^2 -Di^2)/(pi*(Do+Di))/1000; %Hydraulic diameter for 
annulus; 
vis1=(181.94+0.536*(Th -293)+1.22*Pm)*10^ -6; %Viscosity of air  
vis2=(181.94+0.536*(Tk -293)+1.22*Pm)*10^ -6; %Viscosity of air  
L=10/1000; 
  
dd=(pi*((Do/2)^2 -(Di/2)^2));  %Area of annulus 
  
e1=.0996755; % equilibrim pt for mass 1.  
e2=.501; % equilibrim pt for mass 2.  
  
utemp2 = (-2*Ad*x(2)+Ap*x(4))/dd;  
utemp1 = (-2*Ad*x(2)+Ap*x(4))/dd;  
  
%%%Model Fit%%%  
  
xp(1)=x(2); 
xp(2)=(Ad/mdt*((.5*rho*(ffcoef2*utemp2 ^n1 +kh)*utemp2*abs(utemp2)) ... 
    + (.5*rho*(ffcoef1*utemp1 ^n1 +kh)*utemp1*abs(utemp1)))) ... 
    +(mfield/mdt)*(1/(x(1)+e1)^2+1/(x(1)+e1+2*h)^2 -2/(x(1)+e1+h)^2) -g; 
xp(3)=x(4); 
xp(4)=(Ap*Pm/mpt*(1 -1/(1-a*(x(3)+e2)+b*(x(1)+e1) )))-
 (fl+fc*(x(3)+e2)^2)*x(4)/mpt+ ... 
    (mfield/mpt)*(1/(x(3)+e2)^2+1/(x(3)+e2+2*h)^2 -2/(x(3)+e2+h)^2) -g; 
 
  
Bibliography
 Banduric, R. D., and Chen N. J. C. (1984). Nonlinear Analysis of Stirling
 Engine Thermodynamics. Oak Ridge National Laboratory Report, 154.
 Balachandran, B. and Magrab, E. (2003). Vibrations. Cengage Learning.
 Benvenuto, G., and de Monte, F. (1995). Analysis of Free-Piston Stirling
 Engines/Linear Alternator Systems Part 1: Theory. Journal of Propulsion and
 Power, 11: 1036-1046.
 Benvenuto, G., and de Monte, F. (1995). Analysis of Free-Piston Stirling
 Engines/Linear Alternator Systems Part II: Results. Journal of Propulsion and
 Power, 1:1047-1055.
 Benvenuto, G., and de Monte, F. (1996). The E ect of Nonlinear Thermo-
 Fluid-Dynamic Terms on Free-Piston Stirling Machine Stability. Proceedings of
 the 31st Energy Conversion Engineering Conference, Washington, DC.
 Benvenuto, G., de Monte, F., and Farina, F. (1990). Dynamic Behaviour Pre-
 diction of Free-Piston Stirling Engines. Proceedings of the 25th Energy Conversion
 Engineering Conference.
 Berchowitz, D. M., and Redlich, R. W. (1985). Linear Dynamics of Free-
 Piston Stirling Engines. Proceedings of the Institution of Mechanical Engineers,
 199: 203-213.
 Chen, N. J. C., and Gri n, F. P. (1983). A Review of Stirling Engine Mathe-
 matical Models. Oak Ridge National Laboratory Report, 135.
 Chen, N. J. C., and Gri n, F. P. (1983). Linear Harmonic Analysis of Free
 Piston Stirling Engines. Oak Ridge National Laboratory Report, 172.
 Chen, N. J. C., Gri n, F. P., and West, C. D. (1984). Linear Harmonic Anal-
 ysis of Stirling Engine Thermodynamics. Oak Ridge National Laboratory Report,
 155.
 Choudhary, F. (2009). Dynamics of Free Piston Stirling Engines. Master?s Thesis,
 Department of Mechanical Engineering, University of Maryland, College Park,
 MD.
 77
Formosa, F. (2009). Nonlinear Dynamics Analysis of a Membrane Stirling Engine:
 Starting and Stable Operation. Journal of Sound and Vibration, 326: 794-808.
 Gri ths, D. J. (1999). Introduction to Electrodynamics. Prentice Hall, Upper
 Saddle River, NJ.
 Ibrahim, M. B., and Tew, R. C. (2001). Study of Two-Dimensional Compress-
 ible Non Acoustic Modeling of Stirling Machine Type Components. NASA Tech-
 nical Memorandum 211066.
 Kankam, M. D., and Rauch, J. S. (1991). Comparative Survey of Dynamic
 Analyses of Free-Piston Stirling Engines. NASA Technical Memorandum 104491.
 Kankam, M. D., and Rauch, J. S. (1993). Controllability of Free-Piston Stirling
 Engine/Linear Alternator Driving a Dynamic Load. NASA Technical Memoran-
 dum 106497.
 Kankam, M. D., Madi, F. J., Rauch, J. S., and Santiago, W. (1992). A
 Free Piston Stirling Engine/Linear Alternator Controls and Load Interaction Test
 Facility. NASA Technical Memorandum 105825.
 Martini, W. R. (1983). Stirling Engine Design Manual. NASA Technical Memo-
 randum 168088.
 Martini, W. R. (2001). Stirling Technology Development at GRC. NASA Techni-
 cal Memorandum 211315.
 Nayfeh, A. H. (1981). Introduction to Perurbation Techniques. John Wiley and
 Sons, New York.
 Nayfeh, A. H. and Balachandran, B. (1995). Applied Nonlinear Dynamics:
 Analytical, Computational and Experimental Methods. John Wiley and Sons,
 New York.
 Organ, A. J. (1992). Thermodynamics and Gas Dynamics of the Stirling Cycle
 Machine. Cambridge University Press, UK.
 Organ, A. J. (1997). The Regenerator and the Stirling Engine. Mechanical Engi-
 neering Publications Limited, London and Bury St Edmunds, UK.
 78
Senft, J. R., and Walker, G. (1985). Free Piston Stirling Engines. Springer-
 Verlag, NY.
 Ulusoy, N. (1994). Dynamic Analysis of Free Piston Stirling Engines. Ph.D. the-
 sis, Department of Mechanical and Aerospace Engineering, Case Wester Reserve
 University, OH.
 Urieli, I., and Berchowitz, D. M. (1984). Stirling Cycle Engine Analysis. Adam
 Hilger LTD., Bristol, NY.
 Urieli, I. (2010). Stirling Cycle Machine Analysis.
 http://www.ohio.edu/mechanical/stirling/me422.html.
 Vokoun, D., Beleggia, M., Heller, L., and Sittner, P. (2009). Magnetostatic
 Interactions and Forces between Cylindrical Permanent Magnets. Journal of Mag-
 netism and Magnetic Materials. 321: 3758-3763.
 Walker, G. (1980). Stirling Engines. Clarendon Press, Oxford University Press,
 NY.
 79