ABSTRACT Title of Document: DYNAMICS OF FREE PISTON STIRLING ENGINES Farhan Choudhary, M.S., 2009 Directed By: Professor Balakumar Balachandran, Department of Mechanical Engineering Free piston Stirling engines (FPSEs) are examples of closed cycle regenerative engines, which can be used to convert thermal energy into mechanical energy. FPSEs are multi-degree-of- freedom dynamical systems that are designed to operate in a periodic manner. Traditionally, the designed periodic orbits are meta-stable, making the system operation sensitive to disturbances. A preferred operating state would be an attracting limit cycle, since the steady-state dynamics would be unique. In this thesis, it is investigated as to how to engineer a Hopf bifurcation of an equilibrium solution in a FPSE. Through a combination of weakly nonlinear analysis and simulations, it is shown that it is possible to engineer a Hopf bifurcation in a FPSE system. Through the analyses, reduced-order-models are developed on the basis of Schmidt formulations and nodal analysis. This thesis effort could serve as a platform for designing FPSEs which take advantage of nonlinear phenomena in either the beta or double acting alpha configuration. DYNAMICS OF FREE PISTON STIRLING ENGINES By Farhan Choudhary Thesis submitted to the Faculty of the Graduate School of the University of Maryland, College Park, in partial fulfillment of the requirements for the degree of Master of Science 2009 Advisory Committee: Professor Balakumar Balachandran, Chair Associate Professor Gregory Jackson Assistant Professor Nikhil Chopra ? Copyright by Farhan Choudhary 2009 ii Acknowledgements I would like to offer thanks to Dr. Balachandran for his support and guidance throughout my graduate school career. I also thank the committee members, Dr. Jackson and Dr. Chopra, for their insights and assistance in completing this thesis. I also thank all of the members of the Dynamics and Vibrations Group for providing assistance and guidance. Finally, my gratitude goes to my family for their unconditional love and encouragement. iii Table of Contents Acknowledgements......................................................................................................................... ii Table of Contents ..........................................................................................................................iii List of Tables ...................................................................................................................................v List of Figures ............................................................................................................................... vi 1 Introduction............................................................................................................................ 1 1.1 Stirling Engines..........................................................................................................................1 1.2 Configurations of Stirling Engines...........................................................................................3 1.3 Free Piston Stirling Engine.......................................................................................................5 1.4 Existing Literature on FPSE Dynamics...................................................................................6 1.5 Advantages of Nonlinear Design ..............................................................................................8 1.6 Contributions .............................................................................................................................8 1.7 Organization of Thesis ..............................................................................................................9 2 Governing Equations ........................................................................................................... 10 2.1 Nomenclature...........................................................................................................................11 2.2 Comparison with Existing Literature....................................................................................12 2.3 Dynamic Equations .................................................................................................................14 2.4 Schmidt Analysis .....................................................................................................................14 2.4.1 Beta Engine..........................................................................................................................................15 2.4.2 Modified Beta Engine..........................................................................................................................19 2.4.3 Single Cylinder Alpha Engine .............................................................................................................21 2.4.4 Double Acting Alpha Engine: Two Cylinder Case..............................................................................23 2.4.5 Double Acting Alpha Engine: Three Cylinder Case............................................................................25 2.5 Simplified Nodal Analysis.......................................................................................................27 2.5.1 Single Cylinder Alpha Engine .............................................................................................................27 2.5.2 Single Cylinder Case ...........................................................................................................................31 2.5.3 Beta Engine..........................................................................................................................................32 3 Analytical and Numerical Results....................................................................................... 35 3.1 Criterion for Self Starting Engine..........................................................................................35 3.1.1 Linearized Schmidt Model: Root Loci.................................................................................................35 3.2 Modal Coordinates, Method of Multiple Scales, and Numerical Results...........................46 3.2.1 Beta Engine..........................................................................................................................................47 3.2.2 Double Acting Alpha Engine...............................................................................................................58 3.3 Numerical Results from Simplified Nodal Analysis .............................................................64 3.3.1 Single Cylinder Case ...........................................................................................................................65 3.3.2 Beta Engine..........................................................................................................................................68 4 Concluding Remarks and Suggestions for Future Work................................................... 72 iv 4.1 Limit Cycle Motions................................................................................................................72 4.2 Reduced Order Models ...........................................................................................................73 4.3 Future Work ............................................................................................................................73 Appendix A: Computer Programs............................................................................................... 75 Root Locii ..............................................................................................................................................75 Root Locus for Beta Engine...............................................................................................................................75 Root Locus for Double Acting Alpha Engine: Three Cylinder Case ................................................................76 Numerical Integration of Schmidt Model...........................................................................................77 Beta Configuration Free Piston Stirling Engine.................................................................................................77 Parameter Selection............................................................................................................................................79 Transformation of Cubic Damping Terms .........................................................................................................80 Three Cylinder Double Acting Engine...............................................................................................................81 Nodal Methods ......................................................................................................................................83 Isothermal Single Cylinder ................................................................................................................................83 Non-Isothermal Single Cylinder ........................................................................................................................85 Isothermal Simulation of Sunpower RE-1000 Engine .......................................................................................87 Non-Isothermal Simulation of Sunpower RE-1000 Engine...............................................................................89 Bibliography................................................................................................................................. 91 v List of Tables Table 2.1 Comparative table of dynamic analyses ...................................................................... 13 Table 2.2 Comparative table of nodal analyses ........................................................................... 14 Table 3.1 Comparison of analytical and numerical results; ?=0.281........................................... 54 Table 3.2 Comparison of analytical and numerical results; ?=0.15............................................. 54 vi List of Figures Figure 1.1 A simplified schematic displaying the essential components of a Stirling engine. Adapted from [5]............................................................................................................................. 1 Figure 1.2 The ideal Stirling cycle is depicted along with the P-V and T-S diagrams. In moving from 1 to 2, the fluid is compressed isothermally; going from 2 to 3, the fluid is heated at a constant volume. Between states 3-4 the fluid is expanded at a constant temperature. Finally, moving from state 4 to 1 involves constant volume cooling. The cycle is then repeated. Adapted from [5]. .......................................................................................................................................... 2 Figure 1.3 A simplified schematic of an alpha type Stirling engine. The compression and expansion spaces both contain a piston. Adapted from [5]. .......................................................... 3 Figure 1.4 A simplified schematic of a beta type Stirling engine. The displacer serves to force the gas through the heat exchangers; the power piston is contained in the compression space. Both the displacer and power piston are contained within a single cylinder. Adapted from [5]... 4 Figure 1.5 A simplified schematic of a gamma type Stirling engine. The displacer again serves the purpose of forcing the working fluid through the heat exchangers; however, the power piston is contained within a separate cylinder. Adapted from [5]. .......................................................... 5 Figure 1.6 A simplified schematic of a double acting alpha type Stirling engine; note that the compression space and expansion space of each individual unit are located in adjacent cylinders. An equivalent single acting engine would require 8 pistons rather than 4, giving this configuration a higher specific power............................................................................................. 5 Figure 1.7 This is a simple schematic of a free piston Stirling engine. Note that the reciprocating components are mounted on springs rather than being kinematically linked. This engine is one example of a beta configuration...................................................................................................... 6 Figure 2.1 This schematic depicts the well known beta type FPSE for which the governing equations are derived; this configuration is one of a working engine. ......................................... 10 Figure 2.2 A schematic showing the characteristics of a beta configuration FPSE. This configuration is known to work. ................................................................................................... 15 Figure 2.3 This is a variation of the beta engine shown in Figure 2.2. The displacer rod is removed and the displacer is attached to the piston, through a spring, rather than to the casing. 20 Figure 2.4 This configuration is a free piston version of the single acting alpha engine. ........... 22 Figure 2.5 This configuration is of a two cylinder double acting FPSE. There are two thermodynamic cycles that are dynamically coupled. The expansion and compression spaces of each cycle are located in adjacent cylinders. Notice that the bounce volume is now eliminated and the need for only one reciprocating element per thermodynamic cycle. ............................... 23 Figure 2.6 This configuration is a three cylinder double acting FPSE. Although it is not shown in the figure, the expansion space located in the third cylinder connects to the compression space in the first cylinder........................................................................................................................ 25 Figure 2.7 The division of the engine into smaller control volumes is shown. The bounce spaces are each assumed to have a constant mass and constant pressure with respect to time; this assumption is valid for engines in which the change in bounce space volume is small compared to the mean bounce space volume. The heat exchanger volume is assumed to have a spatial gradient in temperature and pressure, although this distribution is assumed to remain relatively constant in time. The expansion and compression spaces are assigned spatially averaged temperature and pressure values that vary in time........................................................................ 28 vii Figure 2.8 A schematic of the single cylinder device analyzed.................................................... 31 Figure 2.9 The nodal analysis is applied to the beta configuration. The expansion and compression spaces are assumed to have time varying thermodynamic properties, while all other parts of the engine are assumed to have constant properties in time............................................ 32 Figure 3.1 A plot of the eigenvalues of the undamped beta engine as the hot side temperature, Th, is increased. There are four eigenvalues that migrate away from the imaginary axis with increasing temperature, two of which have positive real parts..................................................... 36 Figure 3.2 A plot of the eigenvalues of the beta engine with damping applied only to the piston. Rather than cross the imaginary axis, a pair of eigenvalues jumps to the negative real axis, indicating that purely imaginary eigenvalues are not possible in this configuration................... 37 Figure 3.3 A plot of the eigenvalues of the beta engine with damping applied to both the piston and displacer. The pair of right half-plane eigenvalues now cross the imaginary axis, which indicates a possible bifurcation point............................................................................................ 38 Figure 3.4 A plot of the undamped modified beta eigenvalues as the temperature is increased; all four system eigenvalues are purely imaginary and show no significant variation with the heater temperature; the small real part can be attributed to numerical error. .......................................... 39 Figure 3.5 A plot of the modified beta eigenvalues as the damping is increased. Damping is applied to both the piston and displacer; the eigenvalues remain in the left-half plane for any loading........................................................................................................................................... 40 Figure 3.6 A plot of the undamped single acting alpha eigenvalues as the temperature is increased; all four system eigenvalues are purely imaginary and show no significant variation with the heater temperature. The small real part can be attributed to numerical error. ............... 41 Figure 3.7 A plot of the modified single cylinder alpha as the damping is increased. Damping is applied to both pistons; the eigenvalues remain in the left-half plane for any loading. ............... 42 Figure 3.8 A plot of the undamped two cylinder double acting alpha system eigenvalues as the temperature is increased; all four system eigenvalues are purely imaginary and show no significant variation with respect to the heater temperature. The small real part can be attributed to numerical error.......................................................................................................................... 43 Figure 3.9 A plot of the two cylinder double acting system eigenvalues as the damping is increased. Damping is applied to all pistons; the eigenvalues remain in the left-half plane for any loading........................................................................................................................................... 44 Figure 3.10 A plot of the eigenvalues of the undamped, three cylinder double acting alpha engine system as the hot side temperature is increased. The case considered is for three identical cylinders. There are two eigenvalues that migrate away from the imaginary axis with positive real part. ........................................................................................................................................ 45 Figure 3.11 A plot of the eigenvalues of the three cylinder double acting alpha engine system as the damping is increased. The same damping is applied to all three elements. The right half- plane eigenvalues cross the imaginary axis, indicating a possible bifurcation point.................... 46 Figure 3.12 Numerical solution to the equations of motion when the damping is too high. As predicted, the system exhibits decaying motions.......................................................................... 53 Figure 3.13 Numerical solution of the steady-state motion near the bifurcation point; the piston motion is described by the sinusoid with large amplitude whereas the displacer has the smaller amplitude....................................................................................................................................... 53 Figure 3.14 A parametric plot of the steady-state piston and displacer displacements. As predicted, in the limit, the system migrates to a limit cycle. The analytical approximation (inner ellipse) is shown for comparison with the numerical solution (outer ellipse). ............................. 53 viii Figure 3.15 The system still possesses a limit cycle far from the bifurcation point.................... 54 Figure 3.16 A parametric plot of the piston and displacer displacements far from the bifurcation point. The analytical approximation (inner ellipse) is not valid here. ........................................ 54 Figure 3.17 Parametric plot of solution curves for different initial conditions when the system has purely imaginary eigenvalues and no nonlinear damping. As discussed in [5] and [19], the origin is a center, about which there is a continuum of periodic orbits. ....................................... 55 Figure 3.18 Numerical solution of the steady-state motion near the bifurcation point; the piston motion is described by the sinusoid with large amplitude whereas the displacer has the smaller amplitude....................................................................................................................................... 56 Figure 3.19 A parametric plot of the steady-state piston and displacer displacements. As predicted, in the limit, the system migrates to a limit cycle. The analytical approximation (inner elipse) is shown for comparison with the numerical solution (outer ellipse). ............................. 56 Figure 3.20 A parametric plot of the piston and displacer displacements for the parameter values specified in this section. The outer ellipse is the analytical approximation, whereas the inner ellipse is the numerical solution.................................................................................................... 58 Figure 3.21 (a) A plot of the steady-state displacements of the three pistons with respect to time; all three are identical except for a phase shift. (b) A three-dimensional parametric plot of the steady-state displacements............................................................................................................ 62 Figure 3.22 A plot of mode one with respect to time. The maximum value is negligible; thus, the approximation that mode one decays with increasing time is valid. ...................................... 63 Figure 3.23 Plots of the real part of (a) mode two and (b) mode three with respect to time; they are the same................................................................................................................................... 63 Figure 3.24 Plots of the imaginary part of (a) mode two and (b) mode three with respect to time; they are of opposite sign. .............................................................................................................. 64 Figure 3.25 Plots of various quantities versus time obtained from an isothermal nodal analysis of the single cylinder device: (a) piston displacement, (b) working space pressure, and (c) mass of working fluid................................................................................................................................. 66 Figure 3.26 Plots of various quantities versus time obtained from a non-isothermal nodal analysis of the single cylinder device: (a) piston displacement, (b) working space pressure, (c) fluid temperature, and (d) mass of working fluid................................................................................. 67 Figure 3.27 Plots of various quantities versus time obtained from an isothermal nodal analysis of the Beta configuration Sunpower RE-1000 engine: (a) piston displacement, (b) displacer displacement, (c) expansion space pressure, and (d) pressure differential. ................................. 69 Figure 3.28 Plots of various quantities versus time obtained from a non-isothermal nodal analysis of the Beta configuration Sunpower RE-1000 engine: (a) piston displacement, (b) displacer displacement, (c) compression space pressure, and (d) mass of working fluid in working space volumes......................................................................................................................................... 71 1 1 Introduction In this chapter a general discussion of Stirling engines is provided. The author begins by explaining the underlying physical principles at work within the engine and then proceeds to elaborate on different types of engine configurations. Then, the free piston Stirling engine (FPSE) is introduced, followed by a description on the analyses on FPSEs in the literature. The chapter concludes with the motivation for applying nonlinear analysis to FPSEs and a brief description of this thesis. 1.1 Stirling Engines A Stirling engine is a mechanical device used to generate power or to transfer heat from a cool region to a hot region [1]. It is an external combustion engine, and as a whole, it is a closed system; however, the working fluid can flow between various compartments within the engine. In fact, the temperature gradient and volume variations within the engine cause the working fluid to flow throughout the engine in a cyclic manner. This cyclic flow can be used to convert heat into work or to transfer heat against a temperature gradient. A simplified schematic of a Stirling engine is given in Figure 1.1. As illustrated in Figure 1.2, the ideal Stirling cycle consists of the following four reversible processes: (i) isothermal compression, (ii) constant volume regenerative heating, (iii) isothermal expansion, and (iv) constant volume regenerative cooling. In the ideal cycle, the largest possible amount of heat is converted to work; in other words, the ideal efficiency of the Stirling engine is equivalent to the Carnot efficiency. In reality, a real Stirling engine operates with many non- idealities. These include losses due to viscous flow, the continuous motion of the pistons (as opposed to the discontinuous motion required to achieve the ideal cycle), non-isothermal expansion/compression, fluid leakage, and heat transfer losses [2]. All of these serve to reduce the efficiency of a real machine below the maximum theoretical value. It becomes the task of the designer to minimize these losses. Figure 1.1 A simplified schematic displaying the essential components of a Stirling engine. Adapted from [5]. 2 Figure 1.2 The ideal Stirling cycle is depicted along with the P-V and T-S diagrams. In moving from 1 to 2, the fluid is compressed isothermally; going from 2 to 3, the fluid is heated at a constant volume. Between states 3-4 the fluid is expanded at a constant temperature. Finally, moving from state 4 to 1 involves constant volume cooling. The cycle is then repeated. Adapted from [5]. Historically, the Stirling engine was invented in 1816 by Robert Stirling. The regenerative heat exchanger was an innovative feature of this engine. Along with other hot air engines, the Stirling engine found use mostly as a low power engine until the gasoline powered spark-ignition engine and electric motor came into widespread use [1]. By World War I, Stirling engines were scarce; however, in the mid-1930s, Philips Research Laboratory in Eindhoven, the Netherlands, sought to develop a Stirling engine that could provide power in remote locations for radios. After World War II, improvements in batteries and radios had eliminated the original need for the Stirling engine; however, Philips was able to apply its research to develop Stirling engines capable of generating hundreds of horsepower. Since then, research carried out on the Stirling engine has demonstrated its potential for naval, automotive, and space power applications. One of the notable studies is NASA?s continuing efforts to improve free piston Stirling engine design for use in Mars rovers and deep space applications [3], [4]. The advantages that the Stirling engine offers are low noise, reduced emissions (through cleaner combustion), multi-fuel capacity, and long life [1]. The principal limitations of the Stirling 3 engine lies in the availability of appropriate materials to withstand the constant high temperature imposed by the engine. 1.2 Configurations of Stirling Engines Stirling engines are typically configured in one of three ways. The respective configurations are the alpha, beta, and gamma type configurations [5]. In Figure 1.3, the alpha type configuration is illustrated. In this configuration, one uses separate piston cylinder combinations as the expansion and compression spaces; the two cylinders are connected by heat exchangers. The beta configuration makes use of two reciprocating elements within the same cylinder to achieve the thermodynamic cycle; the piston is used to derive power while the displacer motion forces the fluid through heat exchangers. As illustrated in Figure 1.4, the beta configuration is more compact than an equivalent alpha configuration. Figure 1.3 A simplified schematic of an alpha type Stirling engine. The compression and expansion spaces both contain a piston. Adapted from [5]. 4 Figure 1.4 A simplified schematic of a beta type Stirling engine. The displacer serves to force the gas through the heat exchangers; the power piston is contained in the compression space. Both the displacer and power piston are contained within a single cylinder. Adapted from [5]. As illustrated in Figure 1.5, the gamma configuration makes use of two separate cylinders; in one cylinder, the reciprocating motion of the displacer forces the working fluid through the heat exchangers. The piston moves in the adjacent cylinder, in which the working fluid cyclically undergoes compression and expansion. 5 Figure 1.5 A simplified schematic of a gamma type Stirling engine. The displacer again serves the purpose of forcing the working fluid through the heat exchangers; however, the power piston is contained within a separate cylinder. Adapted from [5]. Figure 1.6 A simplified schematic of a double acting alpha type Stirling engine; note that the compression space and expansion space of each individual unit are located in adjacent cylinders. An equivalent single acting engine would require 8 pistons rather than 4, giving this configuration a higher specific power. It is also important to distinguish between single acting engines and double acting engines. A single acting engine represents a complete engine within itself and it is independent of any other engine that may also be contributing power to drive the same load. All of the previously described configurations in this thesis are single acting engines. By contrast, in a double acting engine, multiple piston-cylinder configurations are connected so that the expansion space and compression space of each subsystem are located in adjacent cylinders. A double acting engine can be constructed by connecting several alpha type units together; this configuration reduces the number of reciprocating components by a factor of two and, as a result, this engine has a higher specific power than the single acting counterpart [2]. In Figure 1.6, a schematic of a double acting engine is provided. 1.3 Free Piston Stirling Engine The free piston Stirling engine (FPSE) is a variation of the Stirling engine invented by William Beale in the 1950?s. An illustration of a FPSE system is provided in Figure 1.7. In a kinematic engine, the various moving parts are linked to force the engine components to move in a particular way relative to one another; these constraints serve to make the kinematic engine a single degree-of-freedom dynamical system. In a FPSE, this is not the case; the pistons are not kinematically linked, but rather their motion is induced by fluid pressure variations within the engine. Thus, the pistons can move independently of one another, making a FPSE a multi-degree-of-freedom dynamical system. This arrangement results in a strong coupling between the engine thermodynamics and the engine dynamics. There are several arrangements of a FPSE; for an extensive list, the reader is referred to [2]. 6 Figure 1.7 This is a simple schematic of a free piston Stirling engine. Note that the reciprocating components are mounted on springs rather than being kinematically linked. This engine is one example of a beta configuration. The principal advantages of an FPSE over a kinematic engine are mechanical simplicity, low cost, no leakage, low wear, and self-starting capability [5], [2]. However, the dynamic behavior of a FPSE is much more complicated than in a kinematic engine; this can be a challenge to the designer. 1.4 Existing Literature on FPSE Dynamics In most basic analysis of FPSE dynamics, the Schmidt isothermal model is utilized to relate the thermodynamics and dynamics of the engine [5]. The Schmidt model leads to closed-form expressions for the working gas pressure. Its assumptions are as follows: 1) spatially constant pressure 2) isothermal compression and expansion 3) ideal gas behavior 4) the engine is a closed system These assumptions can be used to determine the working gas pressure in terms of the volume variations. In a single degree-of-freedom kinematic engine, it is then easy to analyze the behavior of the engine as a function of the crank angle. However, in a multi-degree-of-freedom free piston engine, the volume variations are unknown functions of time which must be solved for. Thus, in the case of free piston Stirling engines, the pressure equation is linearized, in addition to any other non-linear terms; in the linearization procedure, one either uses a truncated Taylor series expansion and/or replaces non-linear springs/dampers/loads with linear elements 7 that do the same work as the non-linear elements over some specified portion of the engine cycle. The result is a set of linear differential equations with constant coefficients, which can be analytically solved. Most analyses of Stirling engines begin with the Schmidt model; however, as the Schmidt model is overly simplified, the model is modified to include the dissipative effects of the pressure gradient. This modification utilizes the assumption of quasi-steady flow to relate the pressure drop to the fluid velocity; the fluid velocity is then approximated as a function of the volume variations within the engine [5]. By using the linearized set of governing equations, one can assess the dynamic stability of the FPSE. For a linear system of differential equations, stability relies on assuring that none of the system poles have a positive real part. However, in order to generate power, the engine must have one pair of purely imaginary poles; thus, the system will oscillate without growing amplitude. Several methods have been employed to assure that this is the case; these include the Nyquist criterion, Routh Criterion, and the direct use of the characteristic polynomial [3]. It is often reported that if the non-linearity in an FPSE is considered, it tends to stabilize the motion of the engine. In one study, the effects of the non-linearity in the working gas pressure, pressure gradient, and loading terms on the engine dynamics have been investigated [6]. The analytical relations are again modified versions of the Schmidt model, but they are not linearized. By using center manifold theory and normal form analysis, the parameter values resulting in a stable limit cycle are determined. While the above techniques are analytical in nature, there is also another class of techniques used to study the dynamics of the FPSE known as nodal methods [7]. These methods break the internal working fluid into a number of control volumes; each control volume has a separate pressure and temperature associated with it. The mass flow rate between the control volumes is determined as a function of the pressure difference between those volumes. Applying conservation of mass, momentum, and energy to each of the control volumes, a system of differential equations describing the FPSE thermodynamics is obtained. The equations governing the motion of the engine components are obtained by using force balance methods, as before. The governing equations are then numerically integrated in time in order to determine the dynamics associated with a particular set of parameter values and initial conditions. Since the system of governing equations does not have a closed-form solution, nodal methods do not offer the same insights as the modified versions of the Schmidt model; on the other hand, this method does not overly simplify the engine thermodynamics. In recent work computational fluid dynamics has been used to study multi-dimensional effects present within the working fluid, with the goal of creating computationally efficient two-dimensional models that exceed existing one- dimensional models in accuracy [8]. Alternatively, the same nodal equations can be solved by using Harmonic Analysis, which assumes that all variables vary harmonically about some mean value. After substituting solutions of this form into the governing equations, the products of harmonic functions are replaced with the corresponding truncated Fourier series representations. The final set of equations contain the unknown amplitudes of oscillation, which are solved for [9],[10]. 8 Experimental studies of FPSE behavior are scarce in the open literature. However, NASA has made use of an experimental setup to study the dynamic behavior of the FPSE; their study evaluates the effectiveness of various control mechanisms, the variation of output power with system parameters, and transient response of the engine [11]. NASA has also conducted another experimental study by using the Sunpower RE-1000 engine to collect several experimental measurements of engine behavior [12], and another study in which a hydraulic load is used to extract power from the engine [13]. In addition, there are several cases in which analytical models are applied to working FPSEs in order to compare the analytical results with the known engine behavior [5], [14]. 1.5 Advantages of Nonlinear Design As mentioned, almost all of the previous analytical studies on the free piston Stirling engine have relied on linearizing the governing equations. The solutions to a system of linear ordinary differential equations will either grow in an unstable manner, converge to a point, or include an infinite number of periodic solutions. However, it is important to point that the desired behavior of the FPSE is for it to oscillate indefinitely while producing power. Thus, at least one pair of the system poles for the linear part of the system must be purely imaginary. However, the oscillatory solution of the linearized system of differential equations is neutrally stable rather than asymptotically stable. In other words, the amplitude is not unique and becomes sensitive to perturbations; if the amplitude is sensitive to perturbations, then linear analysis does not rule out the possibility of internal collisions within the engine. The ideal design of a Stirling engine will be to have an asymptotically stable periodic orbit; thus, even when perturbed, the engine would have a tendency to return to the design orbit. Such a motion cannot occur in a linear system; therefore linear analysis cannot provide any guidance on how to achieve a design with a stable periodic orbit. If the design utilizes nonlinearity within the engine, it may then be possible to determine a set of parameters that will produce a stable periodic orbit for a free piston Stirling engine. 1.6 Contributions In this study, it is sought to investigate the stabilizing characteristics of nonlinearities on FPSE dynamics. In the analysis, analytical expressions based on the Schmidt model and the techniques of nodal analysis are used. An attempt is made to simplify the nodal methods in order to construct a method of analysis that is simpler than a full nodal discretization, but also captures more thermodynamic complexity than the Schmidt model. Furthermore, alternative configurations of FPSEs are considered in order to assess what modifications can be made to FPSEs. A self-starting criterion for the FPSE is also developed. In addition, nodal methods are applied to a single piston device in order to examine some of the physical phenomena present in a FPSE system in a simplified manner. 9 1.7 Organization of Thesis The rest of this thesis is organized as follows. In Chapter 2, the development of models used in this study is presented. Both, the Schmidt model and a simplified nodal scheme are discussed. In Chapter 3, the numerical and analytical results obtained by using these models are presented. In Chapter 4, the significance of the thesis contributions are discussed and areas of future work are identified. Finally, the Appendices contain the Mathematica and Matlab programs used to conduct the study. 10 2 Governing Equations The author opens this chapter by providing a brief comparison of this study with similar studies in the literature. This comparison is intended to provide the reader a sense of where this work fits in relation to the broader body of literature pertaining to FPSEs. Then, the author proceeds to explain the derivation of equations used to model the FPSE. The dynamic equations governing the motion of the pistons are obtained by using force balance. It is important to keep in mind that for any oscillatory element, positive displacements are taken to stretch the attached spring. Two models were used to describe the thermodynamics of the engine; the first is the based on the known Schmidt model. The second model is based on nodal methods. In Figure 2.1, a beta type FPSE to be studied is illustrated. Figure 2.1 This schematic depicts the well known beta type FPSE for which the governing equations are derived; this configuration is one of a working engine. 11 2.1 Nomenclature x displacement c damping coefficient P pressure T temperature ? fluid density m mass V volume k spring stiffness R gas constant A area Q heat u fluid velocity KL minor loss coefficient lo initial height of deforming volumes Cp specific heat at constant pressure Cv specific heat at constant volume Subscripts avg averaged quantity e pertaining to expansion space k pertaining to compression space CS pertaining to cross section i arbitrary control volume b pertaining to bounce space k pertaining to cooler h pertaining to heater st storage term flux flux quantity p pertaining to the piston of a beta type configuration d pertaining to the displacer of a beta type configuration 1,2,3 used to indicate cylinders in alpha type engines 2.2 Comparison with Existing Literature Before proceeding to describe the models used in this study, a description is provided in Table 2.1 and Table 2.2 comparing this study with others. There have been several different models applied to FPSEs, some emphasizing either the dynamics or thermodynamics, and still others giving equal emphasis to both; in order to make clear where this study deviates from existing models, these tables have been constructed. In Table 2.1, different levels of analyses emphasizing the dynamics of FPSEs are compared; the analyses include linear analysis [5], nonlinear analysis [6], and the analysis used in this study. In Table 2.2, different levels of analyses using nodal methods are described; the analyses are a full nodal discretization [15], a three node discretization [6], and the two node discretization used in this study. 13 Table 2.1 Comparative table of dynamic analyses Urieli [5] Ulusoy [6] Current Study Dynamic Equations -force balance -force balance -force balance Thermodynamic Equations - Schmidt model for pressure - net mass flow rate approximated by net volumetric flow rate -quasi-steady flow assumed -pressure gradient approximated with equivalent linear dampers - Schmidt model for pressure -net mass flow rate approximated by net volumetric flow rate -quasi-steady flow assumed -pressure gradient approximated by a cubic function -Schmidt model for pressure -pressure gradient not explicitly accounted for Springs -gas springs -gas springs -linear mechanical springs Load -linear damper -damping coefficient varies with displacement amplitude -linear damper Nonlinearity - all nonlinearities linearized - pressure drop, working space pressure, load - cubic damper Analysis -constraints for periodic motions determined -center manifold and normal form transformations used to determine parameters for which Hopf bifurcations arise -root loci constructed to ensure a pair of eigenvalues transversely cross the imaginary axis -method of multiple scales applied to determine the criteria for an attracting limit cycle -analytical approximation of the limit cycle near the bifurcation point developed 14 Table 2.2 Comparative table of nodal analyses Schock [15] Ulusoy [6] Current Study Dynamic Equations -force balance -force balance -force balance Thermodynamic Equations -full nodal discretization to emphasize fluid behavior -3 nodes used to analyze an adiabatic expansion/compression model -2 node simplification used in order to increase thermodynamic accuracy while emphasizing dynamic behavior -fluid flow correlations used to describe behavior in heat exchangers Load -general -linear alternator -linear dampers Analysis -numerical integration -numerical integration -numerical integration 2.3 Dynamic Equations Assuming the engine casing is fixed, the Free Piston Stirling Engine (FPSE) becomes a two degree-of-freedom system. The motion of each of the components is governed by a balance of pressure, spring, and any external load forces. If present, gravitational forces can be taken into account by shifting the static equilibrium position of the system. The force balance results in the following equations: Error! Objects cannot be created from editing field codes. ( 2.1) d d d d d dAm x PdA k x= ??&& ( 2.2) In general, the pressure is a function of spatial location within the FPSE; however, since it is very difficult to know this function, it will be approximated as discussed in the next section. The FLoad term represents the force transmitted to the engine by any device powered by the FPSE; for the time being, it is treated as a generic force. 2.4 Schmidt Analysis The Schmidt analysis offers the simplest representation of the working fluid behavior; it can be used to obtain the pressure variation as a function of the working space volume variations. Furthermore, it is often reported that the spatial pressure variation within the engine is small and that the idealized Schmidt pressure is a reasonable approximation to the working fluid pressure [5]. Thus, the Schmidt model will be applied to the well known Beta configuration; the governing equations will be expressed in nondimensional form and subsequently linearized. The same approach is applied to all of the potential engine configurations; thus, it is shown in detail only for the Beta engine, which is illustrated in Figure 2.1 and Figure 2.2. 15 2.4.1 Beta Engine In a well designed FPSE, the spatial pressure variation is negligible. In such a machine, the spatially averaged pressure of the working fluid can be approximated by using the ideal gas law, isothermal approximation, and the fact that the FPSE is a closed system. It is important to note that the bounce volume is, ideally, sealed off from the compression space; therefore, it is treated as a separate closed system. In addition, the bounce volume is generally very large, meaning that its pressure can be approximated as constant and equivalent to the charge pressure within the engine. Dividing the working spaces and heat exchangers of the engine into zones of different temperature, computing the mass of working fluid within those zones, and adding, one obtains the total mass of working fluid within the engine as: ?= i i ifluid fluid T V R Pm ( 2.3) where the summation is over the different temperature zones. It is important to note that the expansion and compression spaces will both have time varying volumes. By evaluating the summation at equilibrium, one can define the following parameter: Figure 2.2 A schematic showing the characteristics of a beta configuration FPSE. This configuration is known to work. h dd h h kh khr k k k pp T lA T V TT TTV T V T lAS ++ ?++= )( )/ln( ( 2.4) 16 In this analysis, the regenerator does not have a constant spatial temperature; therefore, an effective temperature is computed assuming a linear temperature distribution [5]. With this new parameter S, the previous equations can be solved for the spatially averaged pressure in a beta free piston Stirling engine, leading to 1)( 1 ? ??? ? ??? ? ??? ? ??? ? ??+?= d h d k rd p k pfluid fluid xST A ST AAx ST A S RmP ( 2.5) Linearizing the fluid pressure for small displacements about the static-equilibrium position by using a truncated Taylor series expansion, one obtains ?? ? ? ??? ? ??? ? ??? ? ???+? d h d k rd p k p mfluid xST A ST AAx ST APP )(1 ( 2.6) where the engine charge pressure has been substituted into the equation. By using the linearized form of the pressure, a set of linear ordinary differential equations with constant coefficients can be obtained for studying the behavior of a FPSE system. After substituting the approximate pressure distribution into the dynamic equations the following equation is obtained, for the beta case: ?? ? ? ??? ?= ??? ? ??? ? ? ? ? ? ? ? ? ? ? ? ? ? +?? +??+ +?? ? ? ??? ? ?? ? ?? ?+ ??? ? ??? ? ?? ? ?? ? 0 0 )( )( 0 0 0 0 2 1 2 2 1 2 1 x x ST AAP ST AAAPk ST AAP ST AAP ST AAAP ST APk x x c c x x m m h drm k rdrm d h prm h pdm k rdpm h pm p d p d p & & && && ( 2.7) In this equation, the author has included a damping load on each element in order to study the effects that dissipation has on the system dynamics. The damping matrix in the above equation is, in general, a simplification of the true dissipation occurring within the engine. Next, the dynamic equations are transformed to the state space form; this is done by introducing a vector of length 2n, where n is the number of dynamic equations; the first n elements of the vector represent the displacements and the last n elements represent the corresponding velocities. The transformation is ?? ? ? ? ? ? ?? ? ? ? ? ? = ?? ? ? ? ? ? ?? ? ? ? ? ? d p d p x x x x z z z z & & 4 3 2 1 ( 2.8) 17 which yields the following state space representation: ?? ? ? ? ? ? ?? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ???+?? ?????= ? ?? ? ? ? ? ? ?? ? ? ? ? ? 4 3 2 1 2 4 3 2 1 0)( 0)( 1000 0100 w w w w m c STm AAP STm AAAP m k STm AAP m c STm AAP STm AAAP STm AP m k w w w w d d hd drm kd rdrm d d hd prm p p hp pdm kp rdpm hp pm p p ( 2.9) where the prime is used to indicate differentiation with respect to time. Next, the governing equations are nondimensionalized; the advantage of doing so will become apparent by comparing the dimensional and nondimensional forms of the equations. The governing equations can be nondimensionalized by choosing a characteristic time scale and a characteristic length. The compression space clearance and the piston natural frequency are used for this purpose, and the nondimensionalization is as follows: pii lx ?= ( 2.10) p p k mt ?= ( 2.11) After substituting these new variables and performing the appropriate algebra, the equations of motion become 0)1)1(( 1 =+?+?+?+?? ? pdppppp ????????? ( 2.12) 1((1 ) 1) 0d d d d p d d? ? ? ? ?? ?? ????? ?+ + ? + ? + = ( 2.13) pp p p km c=? ( 2.14) ppd dp d kmm cm=? ( 2.15) pp mp p lk PA=? ( 2.16) ppd mrp d lkm PAm=? ( 2.17) dp pd mk mk=? ( 2.18) ST lA k pp=? ( 2.19) 18 ST lAST lAA h pd k prd ??= )(? ( 2.20) The nondimensional parameters help to compare different effects within the engine; to this end, a physical interpretation of them is provided. The parameter ?p is simply twice the viscous damping factor associated with the piston spring-mass combination. For a given system mass, the viscous damping factor of a damped oscillator can be viewed as a ratio of the dissipative and spring forces present in the system. Thus, the parameter ?d can be viewed as a scaled ratio of the dissipative forces acting on the displacer to the spring forces acting on the piston; the scaling factor is the ratio of the piston mass to the displacer mass. In short, the nondimensional damping parameters compare the dissipative forces with the spring forces in the system. The parameter ?p is the ratio of the static pressure forces to the static spring forces that can act on the piston. On the other hand, ?d is a scaled ratio of the static pressure forces acting on the displacer to the static spring forces acting on the piston, where the scaling factor is the ratio of the piston mass to the displacer masses. Thus, the nondimensional pressure parameters help to measure the relative importance of pressure and spring forces acting on the elements of the system. The parameter ? is quite simple to understand; it is the square of the ratio of the displacer natural frequency to the piston natural frequency. The parameter ? is the weighted volume fraction of the compression space to the total engine volume, where the weight assigned to each volume is the inverse of the temperature associated with that volume. Since the Schmidt model assumes a spatially constant pressure and temporally constant temperatures, this ratio is equivalent to the mass fraction of the working fluid inside the compression space. A cold volume will have to have more working fluid inside of it than a hot volume; furthermore, other things being equal, the pressure of a fluid changes more in response to a volume change if the mass of that fluid is increased. Differentiating the ideal gas law with respect to volume, one can see that the magnitude of this partial derivative increases linearly with respect to mass; that is 2P mRTV V? = ?? ( 2.21) Tying these facts together, the parameter ? measures how much the mass distribution, and thus the working fluid pressure, will change in response to the piston motion. The parameter ? is also a weighted volume ratio; however, the volumes of interest are those effected by the displacer motion, namely the compression space (excluding the displacer rod volume) and the expansion space. Thus, it measures how much the pressure will change in response to the displacer motion. On examining the parameters in more detail, it is noted that the amplitudes have been scaled so that they are less than unity (provided that the expansion space clearance is not much larger than the compression space clearance); in addition, the parameters ? and ? also have magnitudes which are less than unity. Linearizing the nondimensional governing equations about the static- equilibrium position, one obtains 19 (1 ) 0p p p p p p d? ? ? ? ? ? ? ???? ?+ + + ? = ( 2.22) ( ) 0d d d d p d d? ? ? ? ?? ? ? ? ??? ?+ + + ? = ( 2.23) Applying the substitution 1 2 3 4 p d p d xz xz xz xz ? ?? ? ? ?? ? ? ?? ? = ? ?? ? ? ? ?? ? ?? ? ? ? ( 2.24) the state space form of the linearized system is obtained. 1 1 2 2 3 3 4 4 0 0 1 0 0 0 0 1 (1 ) 0 0 p p p d d d z z z z z z z z ? ? ? ? ? ? ? ? ? ? ? ? ?? ? ? ? ? ? ? ?? ? ? ? ? ?? ?= ? ? ? ?? ?? + ? ? ? ? ?? ?? ? + ?? ? ? ?? ? & & & & ( 2.25) where the overdots indicate differentiation with respect to time. In the subsequent sections, the same analysis is applied to different engine configurations; since the analysis is presented in detail for the Beta engine configuration, less detail is shown for other engine configurations. However, it is important to note that the key steps are the same and that the nondimensional parameters have similar physical meanings. 2.4.2 Modified Beta Engine The beta engine shown in Figure 2.2 has a rod that passes through the piston; in order to prevent significant leakages, this interface must have as close a tolerance as possible while preventing physical contact between the piston and displacer rod. A natural improvement maybe to remove the rod and spring the displacer to the piston, as opposed to the casing. Such a configuration is shown in Figure 2.3; the Schmidt model is applied to this configuration to first determine the fluid pressure as 1 1 1p pd d d dfluid m p d m p d k k h k k h A AA A A AP P x x P x x ST ST ST ST ST ST ?? ? ? ?? ? ? ? = ? + ? ? + ? ?? ? ? ?? ? ? ?? ? ? ? ? ? ? ?? ? ? ? ( 2.26) h dd h h kh khr k k k pp T lA T V TT TTV T V T lAS ++ ?++= )( )/ln( ( 2.27) 20 Figure 2.3 This is a variation of the beta engine shown in Figure 2.2. The displacer rod is removed and the displacer is attached to the piston, through a spring, rather than to the casing. The governing equations can be nondimensionalized again by choosing a characteristic time scale and a characteristic length. To this end, the compression space clearance and the piston natural frequency are used for this purpose, leading to pii lx ?= ( 2.28) p p k mt ?= ( 2.29) On substituting these new variables and performing the appropriate algebra, the equations of motion become 1((1 ) 1) 0p p p p d d? ? ? ? ?? ?? ???? ?+ + ? + ? + = ( 2.30) ( ) 0d d d d p? ? ? ? ? ??? ?+ + ? = ( 2.31) pp p p km c=? ( 2.32) ppd dp d kmm cm=? ( 2.33) p m p p A P k l? = ( 2.34) 21 p p k A l T S? = ( 2.35) d p d p k h A l A l T S T S? = ? ( 2.36) dp pd mk mk=? ( 2.37) Upon linearizing the nondimensional equations about the static-equilibrium position and applying the state-space representation, the result is 1 1 2 2 3 3 4 4 0 0 1 0 0 0 0 1 (1 ) 0 0 p p p d z z z z z z z z ? ? ? ? ? ? ? ? ? ?? ? ? ? ? ? ? ?? ? ? ? ? ?? ?= ? ? ? ?? ?? + ? ? ? ? ?? ?? ?? ? ? ?? ? & & & & ( 2.38) 2.4.3 Single Cylinder Alpha Engine It is often the case that FPSEs are of a beta configuration. In Figure 2.4, an alpha configuration is shown. This configuration is examined to determine what advantages it may offer; applying the Schmidt model, the pressure of the working fluid is obtained as follows: 1 1 1 2 2 1 1 2 21 1 fluid m m h k h k A x A x A x A xP P P ST ST ST ST ?? ?? ? ? ? = ? + ? + +? ?? ? ? ?? ? ? ? ? ?? ? ( 2.39) hh h kh khr k k k T lA T V TT TTV T V T lAS 1122 )( )/ln( ++ ?++= ( 2.40) As in the beta case, it is assumed that bounce volume is large enough that it can be treated as having a constant pressure equivalent to the charge pressure of the engine. For the single cylinder alpha engine, a simple case in which both sets of pistons and springs have the same mass, area, and stiffness is examined. Thus, one can adopt the following substitution: i ix l?= ( 2.41) mt k?= ( 2.42) 22 Figure 2.4 This configuration is a free piston version of the single acting alpha engine. On substituting these new variables and performing the appropriate algebra, the equations of motion become 11 1 2 1 1((1 ) 1) 0? ?? ? ?? ?? ???? ?+ + ? ? ? + = ( 2.43) 12 2 2 1 2((1 ) 1) 0? ?? ? ?? ?? ???? ?+ + ? ? ? + = ( 2.44) cmk? = ( 2.45) mAPkl? = ( 2.46) k Al T S? = ( 2.47) h Al T S? = ( 2.48) Linearizing about the static-equilibrium position and transforming into the state-space representation, one obtains: 1 1 2 2 3 3 4 4 0 0 1 0 0 0 0 1 (1 ) 0 (1 ) 0 z z z z z z z z ?? ?? ? ?? ?? ? ? ? ? ?? ? ? ? ? ?? ? ? ? ? ?? ?= ? ? ? ?? ?? + ? ? ? ? ? ?? ?? ? + ?? ?? ? ? ? & & & & ( 2.49) 23 2.4.4 Double Acting Alpha Engine: Two Cylinder Case Figure 2.5 This configuration is of a two cylinder double acting FPSE. There are two thermodynamic cycles that are dynamically coupled. The expansion and compression spaces of each cycle are located in adjacent cylinders. Notice that the bounce volume is now eliminated and the need for only one reciprocating element per thermodynamic cycle. For the double acting alpha engine shown in Figure 2.5, there are two separate Stirling engines with coupled dynamics. Thus, each Stirling cycle is denoted by the letter a or b. The dynamic components are denoted by 1 or 2. A double acting configuration offers higher specific power over the single acting equivalent; thus, it is natural to ask if a double acting FPSE can be constructed. The Schmidt model yields the following for the pressure variation in the different thermodynamic cycles. ?? ? ? ??? ? +?? ??? ? ??? ? ?+= ? aaaa kaha aech kaha aechafluid TS xA TS xAP TS xA TS xAPP 2211 arg 1 2211 arg 11 ( 2.50) a e a a aa aaa a a a c h a h h kh khr k k k a a T lA T V TT TTV T V T lAS 21 )( )/ln( ++ ?++= ( 2.51) ?? ? ? ??? ? ?+? ??? ? ??? ? +?= ? bbbb kbhb bech kbhb bechbfluid TS xA TS xAP TS xA TS xAPP 2211 arg 1 2211 arg 11 ( 2.52) 2 1ln( / )( )c b b b b b e b b b b b b b k r h k h b b k k h k h h A l V V T T V AlS T T T T T T= + + + +? ( 2.53) 24 Realistically, the mean pressure within both separate Stirling cycles will be the same because of leakages. However, if one assumes that both cylinders have identical parameter values, then the stiffness matrix becomes symmetric with restoring forces; this system will have purely imaginary eigenvalues. Here, the author nondimensionalized the system without this assumption in order to investigate the more general case. The following substitution is implemented: 1 ci i x l?= ( 2.54) 1 1 mt k?= ( 2.55) The governing equations read as: 1 11 1 1 1 2 1 1 2 1((1 ) (1 ) ) 0b b a a? ? ? ? ? ? ? ? ? ? ? ? ?? ??? ?+ + + ? ? + ? + = ( 2.56) 1 12 2 2 2 1 2 2 1 2((1 ) (1 ) ) 0a a b b? ? ? ? ? ? ? ? ? ? ? ? ??? ??? ?+ + + ? ? + ? + = ( 2.57) 11 1 1 c m k? = ( 2.58) 1 22 2 1 1 m c m m k? = ( 2.59) 11 1 1c mA P k l? = ( 2.60) 1 22 2 1 1c mm A P m k l? = ( 2.61) 1 1c a a k a Al T S? = ( 2.62) 2 1c b b k b A l T S? = ( 2.63) 2 1c a a h a A l T S? = ( 2.64) 1 1c b b h b A l T S? = ( 2.65) 2 1 1 2 k m k m? = ( 2.66) Linearizing about the static equilibrium position and applying the state space transformation, one obtains 25 1 1 2 2 1 1 13 3 2 2 24 4 0 0 1 0 0 0 0 1 1 ( ) ( ) 0 ( ) ( ) 0 b a a b b a a b z z z z z z z z ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?? ? ? ? ? ? ? ?? ? ? ? ? ?? ?= ? ? ? ?? ?? ? + + ? ? ? ? ?? ?+ ? ? + ?? ? ? ?? ? & & & & ( 2.67) 2.4.5 Double Acting Alpha Engine: Three Cylinder Case A three cylinder double acting configuration Figure 2.6 is also examined in order to determine whether the system dynamics is altered by the addition of more cylinders. The pressure variation in each thermodynamic system is derived by using the Schmidt model, and this reads as ?? ? ? ??? ? +?? ??? ? ??? ? ?+= ? aaaa kaha aech kaha aechafluid TS xA TS xAP TS xA TS xAPP 2211 arg 1 2211 arg 11 ( 2.68) Figure 2.6 This configuration is a three cylinder double acting FPSE. Although it is not shown in the figure, the expansion space located in the third cylinder connects to the compression space in the first cylinder. 26 a e a a aa aaa a a a c h a h h kh khr k k k a a T lA T V TT TTV T V T lAS 21 )( )/ln( ++ ?++= ( 2.69) 1 3 3 3 32 2 2 21 1 b b b b b b fluid b m m b k b h b k b h A x A xA x A xP P P S T S T S T S T ?? ? ? ? = + ? ? ? +? ? ? ?? ? ? ? ? ? ? ? ( 2.70) b e b b bb bbb a a b c h b h h kh khr k k k b b T lA T V TT TTV T V T lAS 32 )( )/ln( ++ ?++= ( 2.71) 1 3 3 3 31 1 1 11 1 c c c c c c fluid c m m c k c h b k c h A x A xA x A xP P P S T S T S T S T ?? ? ? ? = + ? ? ? +? ? ? ?? ? ? ? ? ? ? ? ( 2.72) c e b c cc ccb c a c c h c h h kh khr k k k c c T lA T V TT TTV T V T lAS 13 \ )( )/ln( ++ ?++= ( 2.73) Again, the author considers that the leakage will force the long term operation of the device to have the same mean pressure within the different cycles. For simplicity, it is assumed that all three cylinders are identical, in that the stiffnesses, masses, clearances, temperature differentials, and so on are the same for each cylinder. Next, a similar scheme is applied to nondimensionalize the three cylinder double acting alpha engine. Doing so allows the author to reduce the number of parameters. Thus, setting i ix l?= ( 2.74) mt k?= ( 2.75) the following equations of motion, are obtained: 1 11 1 3 1 1 2 1((1 ) (1 ) ) 0? ?? ? ?? ?? ?? ?? ?? ??? ?+ + + ? ? + ? + = ( 2.76) 1 12 2 1 2 2 3 2((1 ) (1 ) ) 0? ?? ? ?? ?? ?? ?? ?? ??? ?+ + + ? ? + ? + = ( 2.77) 1 13 3 2 3 3 1 3((1 ) (1 ) ) 0? ?? ? ?? ?? ?? ?? ?? ??? ?+ + + ? ? + ? + = ( 2.78) cmk? = ( 2.79) mAPkl? = ( 2.80) k Al T S? = ( 2.81) h Al T S? = ( 2.82) Linearizing about the static-equilibrium position and using a state-space representation, one obtains 27 1 1 2 2 3 3 4 4 5 5 6 6 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 1 ( ) 0 0 1 ( ) 0 0 1 ( ) 0 0 z z z z z z z z z z z z ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?? ? ? ? ? ?? ? ? ? ? ?? ? ? ? ? ?? ?= ? ? ? ?? ?? ? + ? ? ? ? ?? ? ? ? ? ?? ?? ? + ? ? ? ? ?? ?? ? ? ?? ? + ? ? ?? ? ? ? & & & & & & ( 2.83) 2.5 Simplified Nodal Analysis Nodal techniques are based on applying conservation laws to finite control volumes within the engine. Typically, the expansion, compression, and bounce spaces will constitute one node each while the heat exchangers contain several nodes. While such a discretization offers precision, it demands a great deal of computational effort. Here, an effort is made to eliminate all nodes aside from the expansion and compression spaces by using algebraic expressions to represent the heat exchangers. The potential advantages of such a formulation are computational simplicity in comparison to the full nodal techniques, while capturing more of the thermodynamic complexity than the Schmidt model. Various formulations for some select engines are shown below. The single cylinder alpha engine is shown in detail, whereas less detail is presented for the subsequent configurations. 2.5.1 Single Cylinder Alpha Engine The working fluid inside of an FPSE is typically air, hydrogen, or helium, all of which approximate ideal gas behavior very well. Thus, assuming the engine operation to be quasi- equilibrium, one can apply the ideal gas equation of state at any point within the working fluid. P RT?= ( 2.84) If the ideal gas law is applied to a finite control volume, then it can be used to assign averaged thermodynamic properties to that control volume. PV mRT= ( 2.85) Thus, the bounce volumes, expansion space, heat exchanger volume, and compression space are all divided into control volumes with an average pressure and temperature. 28 Figure 2.7 The division of the engine into smaller control volumes is shown. The bounce spaces are each assumed to have a constant mass and constant pressure with respect to time; this assumption is valid for engines in which the change in bounce space volume is small compared to the mean bounce space volume. The heat exchanger volume is assumed to have a spatial gradient in temperature and pressure, although this distribution is assumed to remain relatively constant in time. The expansion and compression spaces are assigned spatially averaged temperature and pressure values that vary in time. Next, the conservation of mass principle is applied to the expansion and compression spaces. For an arbitrary, deforming control volume with a spatially constant density, the change in stored mass is given by 2st P PT PVm V V VRT RT RT? ? ? ?= + = ? +? ? ? ? & & &&&& ( 2.86) The expression on the right hand-side is obtained by substituting the ideal gas law for density. It is assumed that no control volumes other than the expansion and compression spaces have any net change in stored mass; thus, any mass that leaves the expansion space must end up in the compression space. Assuming the pressure and temperature profiles within the dead volume (thus the fluid density) do not change significantly throughout the cycle, this approximation is valid. Mathematically, this implies the following: 1 2st st m m= ?& & ( 2.87) In order to determine the rate at which mass leaves/enters a control volume, a quasi-steady flow with negligible entry length is assumed. If this assumption applies and the minor losses are the dominant loss mechanism, the pressure drop between two points is related to the average fluid velocity by the following equation 29 2 ( ) 2 b a b L a uP P Sign u K ?? = ? ( 2.88) where the fluid velocity is positive if going from a to b. The mass flux is related to the average fluid velocity by the following equation: csm A u?=& ( 2.89) On substituting and solving for the mass flux, one obtains the following relation 2( ) avga b cs a b a bb L a m A Sign P P P P K ? ? = ? ??& ( 2.90) Applying these equations to the expansion and compression space mass balances while solving for the time derivative of pressure, one obtains the following differential equations: 1 2 1 1 11 1 1 1 1 m RT T VP P V T V ? ? ??= + ?? ? ? ? & &&& ( 2.91) 1 2 2 2 22 2 2 2 2 m RT T VP P V T V ? ? ?= + ?? ? ? ? & &&& ( 2.92) 1 2 1 2 1 22 1 2( ) avg cs L m A Sign P P P P K ? ? = ? ??& ( 2.93) To determine the average temperature of a control volume, an energy balance is applied. The work and heat input/output of a control volume along with net enthalpy convection are balanced by the energy stored within the control volume. Assuming constant specific heats with both enthalpy and internal energy defined to be zero at zero Kelvin, the following equations are obtained: ( )v p flux inmC T mC T PV Q? = ? + &&& ( 2.94) The flux temperature depends on whether the mass flux enters or leaves the control volume. A heating term is included but there may be no external heating applied to a control volume. In order to simplify the analysis, it is assumed that the heat exchange process is perfect. In other words, the exit temperature from the heater will always be some value Th while the exit temperature from the cooler will always be some Tk. Applying the above equation, the previous assumption, and solving for the time derivative of temperature, one obtains the following equations: 30 v invfluxp Cm QVPTCTCmT 1 11121 1 11 )( &&&& +???= ? ( 2.95) v invfluxp Cm QVPTCTCmT 2 22221 2 22 )( &&&& +??= ? ( 2.96) where the flux temperatures are given by 1 1 1 2 1 2 for 0 or for 0flux hT T m T m? ?= > <& & ( 2.97) 2 1 2 2 1 2 for 0 or for 0flux kT T m T m? ?= > <& & ( 2.98) In the flux temperature equation, Th and Tk represent the heater and cooler temperatures, respectively. The above equations are assembled to obtain a final set of equations used to analyze the FPSE; they are repeated below: 1 1 1 1 1 1( )p b Loadm x P P A k x F= ? ? +&& ( 2.99) 2 2 2 2 2 2( )p bm x P P A k x= ? ?&& ( 2.100) 1 2 1 1 11 1 1 1 1 m RT T VP P V T V ? ? ??= + ?? ? ? ? & &&& ( 2.101) 1 2 2 2 22 2 2 2 2 m RT T VP P V T V ? ? ?= + ?? ? ? ? & &&& ( 2.102) 1 2 1 2 1 22 1 2( ) avg cs L m A Sign P P P P K ? ? = ? ??& ( 2.103) v invfluxp Cm QVPTCTCmT 1 11121 1 11 )( &&&& +???= ? ( 2.104) v invfluxp Cm QVPTCTCmT 2 22221 2 22 )( &&&& +??= ? ( 2.105) 1 1 1 2 1 2 for 0 or for 0flux hT T m T m? ?= > <& & ( 2.106) 2 1 2 2 1 2 for 0 or for 0flux kT T m T m? ?= > <& & ( 2.107) i ii i PVm RT= ( 2.108) ( ) ii i o i V A l x= ? ( 2.109) 31 2.5.2 Single Cylinder Case In order to gain some insight as to the interaction between the dynamics and thermodynamics within an FPSE, a single cylinder system, depicted in Figure 2.8, is also analyzed. The equations governing this system can be derived in the same way as the full FPSE equations. The single cylinder equations arise as a special case of the equations governing the alpha configuration; the constraint is that the second ?cylinder? (which is in reality the ambient air) has constant properties. The governing set of equations is listed below: ( )p b Loadm x P P A k x F= ? ? +&& ( 2.110) m RT T VP PV T V? ?= + ?? ? ? ? & &&& ( 2.111) 1 12 1 2( ) flux cs o o L m A Sign P P P P K ?= ? ? ? & ( 2.112) ( ) ( )p flux v in w o v m C T C T PV Q hA T TT m C ? ? + ? ?= &&&& ( 2.113) for 0 or for 0flux oT T m T m= < >& & ( 2.114) for 0 or for 0flux oP m mRT? ?= < >& & ( 2.115) PVm RT= ( 2.116) ( )oV A l x= ? ( 2.117) Figure 2.8 A schematic of the single cylinder device analyzed. 32 2.5.3 Beta Engine Figure 2.9 The nodal analysis is applied to the beta configuration. The expansion and compression spaces are assumed to have time varying thermodynamic properties, while all other parts of the engine are assumed to have constant properties in time. The simplified nodal analysis outlined for the single acting alpha configuration can be easily adapted to the beta configuration. The real difference lies in the equations governing the volume variations, as well as in the pressure force acting on each reciprocating element. ( )p p b c p p p Loadm x P P A k x F= ? ? +&& ( 2.118) ( ) ( )d d c e d b c r d dm x P P A P P A k x= ? + ? ?&& ( 2.119) ( ) ( )c p p p d r dV A l x A A x= ? + ? ( 2.120) ( )e d d dV A l x= ? ( 2.121) While this is the only modification necessary, the author makes a few further alterations; an isothermal assumption is used and the behavior is compared to that of a non-isothermal model. Also, a simpler expression for the mass flow rate between different control volumes is used. These modifications are based on information from references [10] and [16]. 2.5.3.1 Isothermal Case In the isothermal model, the energy equation is not needed. Each control volume is assigned a temperature that is assumed to remain constant in time; thus, the expansion space is at the heater 33 temperature and the compression space is at the cooler temperature. Differentiating the ideal gas law in time leads to iiiiii VPVPmRT &&& += ( 2.122) Then, assuming a quasi-steady flow, from [10], one has the following equation for mass flow- rate ( ) a b a b a b P Pm k ? ? ?= ?& ( 2.123) which is in contrast to the equation used for the single acting alpha engine; this equation is expected to have greater numerical stability and for an appropriately selected constant, may not sacrifice too much accuracy. Thus, the isothermal model only requires the following equations: ( )p p b c p p p Loadm x P P A k x F= ? ? +&& ( 2.124) ( ) ( )d d c e d b c r d dm x P P A P P A k x= ? + ? ?&& ( 2.125) ( ) ( )c p p p d r dV A l x A A x= ? + ? ( 2.126) ( )e d d dV A l x= ? ( 2.127) iiiiii VPVPmRT &&& += ( 2.128) ( ) a b a b a b P Pm k ? ? ?= ?& ( 2.129) 2.5.3.2 Non-Isothermal Case Retaining the simplified mass flow rate equation from above, the author removes the assumption of isothermal behavior. In the case of the single acting alpha engine, the ideal gas law is substituted into the conservation of mass equation; here, based on [16], the author instead substituted it into the energy equation. Thus, the result is ( )p p b c p p p Loadm x P P A k x F= ? ? +&& ( 2.130) ( ) ( )d d c e d b c r d dm x P P A P P A k x= ? + ? ?&& ( 2.131) ( ) ( )c p p p d r dV A l x A A x= ? + ? ( 2.132) ( )e d d dV A l x= ? ( 2.133) ( ) cc c flux c c c P m RT PVV?= ?& && ( 2.134) ( ) ee e flux e e e P m RT PVV?= ?& && ( 2.135) 34 p v C C? = ( 2.136) ( ) a b a b a b P Pm k ? ? ?= ?& ( 2.137) ( ) c e c e c P Pm k ? ?= ?& ( 2.138) ( ) e c e c e P Pm k ? ?= ?& ( 2.139) for 0 or for 0 eflux h e e e T T m T m= > <& & ( 2.140) for 0 or for 0 cflux k c c c T T m T m= > <& & ( 2.141) i ii i PVT m R= ( 2.142) In the flux temperature equation, Te and Tc represent the expansion space and compression space temperatures, respectively; they are computed using the ideal gas law. In this chapter, equations governing a FPSE system have been presented. The development of these equations on the basis of the Schmidt model and nodal analysis has been discussed. The different differential system models presented in this chapter are investigated in the next chapter. 35 3 Analytical and Numerical Results In this chapter, the different numerical and analytical results obtained in this study are presented. First, a conceptual approach based on the Schmidt model and the theory of dynamical systems is presented to determine which of the alternative engine configurations developed might exhibit a Hopf bifurcation [17], leading to a limit cycle. Next, cubic damping is introduced into those configurations that may have a Hopf bifurcation and the method of multiple scales is used to determine an analytical criterion for a limit cycle motion. This analytical relation is then compared with numerical results. At the end of the chapter, the results obtained from the simplified nodal methods are presented. The programs used to obtain the numerical results are provided in the Appendices. 3.1 Criterion for Self Starting Engine 3.1.1 Linearized Schmidt Model: Root Loci Qualitatively, if all dissipative effects within an engine are neglected, including the external load, then the engine should become dynamically unstable. In other words, there is a net buildup of energy within the engine in the absence of dissipative effects, resulting in motions with increasing amplitudes. Mathematically, this means that the linearized, undamped equations of motion have one or more eigenvalues with positive real part. Therefore, it is instructive to examine the equations of motion when dissipation is neglected; if the parameters can be varied to obtain eigenvalues with positive real parts, then the configuration may be capable of overcoming dissipative effects. Otherwise, the configuration does not hold promise. Thus, the essential features of the undamped equations of motion for each engine configuration are captured and parameters are varied to see if the above criterion can be met. Then, linear dissipation is introduced into the system in order to see if a pair of complex conjugate eigenvalues transversely cross the imaginary axis, which implies the possibility of a Hopf bifurcation [17]. 3.1.1.1 Beta Engine The Beta configuration is known to work; thus the aforementioned qualitative analysis is applied to this configuration. Two parameters depend on the heater temperature; as the heater temperature increases, the limiting behavior of those parameters is obtained as lim lim ln( / ) ( ) h h p p p p T T p p p p k k r h k h d d k k k h k h h A l A l A l A l VV V T T V A lT T T T T T T ? ?? ?? = = +? ? + + + +? ?? ? ? ( 3.1) ( ) ( )lim lim h h d r p d p d r p T T k h p p k A A l A l A A l T S T S A l V? ??? ?? ? ?? ?= ? = ? ? ? +? ? ( 3.2) 36 Figure 3.1 A plot of the eigenvalues of the undamped beta engine as the hot side temperature, Th, is increased. There are four eigenvalues that migrate away from the imaginary axis with increasing temperature, two of which have positive real parts. Noting, also, that the magnitudes of both of these parameters must always be less than unity, a root locus is constructed as shown in Figure 3.1. Arrows have been used in this figure and the following figures to show the directions in which the eigenvalues move with respect to the considered parameter variation. For simplicity variation in ? is neglected and ? is made to approach ? from below, which is the limiting behavior computed above. As expected, the undamped engine is an unstable system; furthermore, it is possible to obtain roots with a positive real part and nonzero imaginary part, resulting in growing oscillations. In order to see whether the engine can be stabilized with dissipation, linear damping is applied just to the piston; this would correspond to a device being powered by the engine. In other words, this damping is used to represent the dissipative effect of a load. 37 Figure 3.2 A plot of the eigenvalues of the beta engine with damping applied only to the piston. Rather than cross the imaginary axis, a pair of eigenvalues jumps to the negative real axis, indicating that purely imaginary eigenvalues are not possible in this configuration. As shown in Figure 3.2, a pair of eigenvalues approaches the imaginary axis from the right; however, rather than cross the imaginary axis, the eigenvalues jump to the negative real axis. Since, this is not the desired behavior, a damper is now applied to both the piston and displacer. The corresponding results are shown in Figure 3.3. 38 Figure 3.3 A plot of the eigenvalues of the beta engine with damping applied to both the piston and displacer. The pair of right half-plane eigenvalues now cross the imaginary axis, which indicates a possible bifurcation point. Now, the roots can all be forced into the left half plane, and a pair can be made to cross the imaginary axis. These root loci indicate that both the displacer and piston must have a dissipative force acting upon them in order for a Hopf bifurcation to arise. While in reality, the viscous effects within the engine will provide this dissipation to the displacer, the exact form of that dissipation is unknown; at best, it can be approximated, perhaps as in [5]. However, the device cannot be stabilized without it. 3.1.1.1.1 Potential Function When the Beta engine is represented by the Schmidt model and dissipation is neglected, the system can be described by a potential function. If the equation governing piston motions is multiplied by the integrating factor d p ??? ??= ? ( 3.3) the following potential function can be derived: 2 2 ( , ) ln 1 2 2pd dp d p d p p d dV ?? ?? ? ?? ?? ?? ? ? ? ??= ? + ? ? + + ( 3.4) 39 Figure 3.4 A plot of the undamped modified beta eigenvalues as the temperature is increased; all four system eigenvalues are purely imaginary and show no significant variation with the heater temperature; the small real part can be attributed to numerical error. Although it is not used further in this study, this potential function could be used in constructing a Lyapunov function to further analyze the dynamics of the Beta engine. 3.1.1.2 Rodless Beta Engine Here, the effects of removing the displacer rod in the Beta configuration are examined. Again, by increasing the temperature, one obtains lim lim ln( / ) ( ) h h p p p p T T p p p p k k r h k h d d k k k h k h h A l A l A l A l VV V T T V A lT T T T T T T ? ?? ?? = = +? ? + + + +? ?? ? ? ( 3.5) lim lim h h d p d p d p T T k h p p k A l A l A l T S T S A l V? ??? ?? ? ?= ? = ? ? ? +? ? ( 3.6) Again, ? is assumed to be constant and ? approaches ? from below. 40 As shown in Figure 3.4, the undamped configuration has four purely imaginary roots as the temperature is varied; thus, the undamped linearized system cannot be made unstable. While in theory, the nonlinearities can cause growing or decaying oscillations, the undamped case is not Figure 3.5 A plot of the modified beta eigenvalues as the damping is increased. Damping is applied to both the piston and displacer; the eigenvalues remain in the left-half plane for any loading. realistic because some dissipation exists within the device. Therefore, this configuration is not promising. Applying a load on both, the piston and displacer, the results shown in Figure 3.5 are obtained. As expected, the damped system has eigenvalues only in the left half-plane, meaning a loaded version of this device would have decaying motions. 3.1.1.3 Single Cylinder Alpha Engine Next, the single cylinder alpha configuration is considered. In the case that the expansion space and compression space have identical dimensions, the limiting behavior of the parameters becomes lim lim ln( / ) ( ) h hT T kk r h k h k k k h k h h Al Al Al VV V T T VAl AlT T T T T T T ? ?? ?? = = +? ? + + + +? ?? ? ? ( 3.7) 41 lim lim 0 h hT T h Al T S??? ??= = ( 3.8) As shown in Figure 3.6, the undamped equations of this configuration also have four purely imaginary roots; therefore, applying the same reasoning as in the modified beta configuration, it is not a promising configuration. Figure 3.6 A plot of the undamped single acting alpha eigenvalues as the temperature is increased; all four system eigenvalues are purely imaginary and show no significant variation with the heater temperature. The small real part can be attributed to numerical error. 42 Figure 3.7 A plot of the modified single cylinder alpha as the damping is increased. Damping is applied to both pistons; the eigenvalues remain in the left-half plane for any loading. As shown in Figure 3.7, once linear damping is added to the system, it forces the eigenvalues into the left half plane. Thus, there is no balance point where the system has purely imaginary eigenvalues in the presence of dissipation. 3.1.1.4 Double Acting Alpha Engine: Two Cylinder Configuration Next, the two cylinder alpha configuration is considered. For generality, the two cylinders are not assumed to be the same; in this case, the limiting behavior of the parameters becomes 1 1 11 2 lim lim ln( / ) ( ) c c h ha a c ac a a a a a e a a a a a a a a a aT T a ka k r h k h a k k k h k h h Al A l Al VAl V V T T V A lT T T T T T T ? ?? ?? = = +? ? + + + +? ?? ?? ? ? ( 3.9) 2 2 22 1 lim lim ln( / ) ( ) c c h hb b c ac b b b b b e b b b b b b b a a bT T b kb k r h k h b k k k h k h h A l A l A l VA l V V T T V AlT T T T T T T ? ?? ?? = = +? ? + + + +? ?? ?? ? ? ( 3.10) 43 2 1lim lim 0c h ha a a aT T h a A l T S??? ??= = ( 3.11) 1 1lim lim 0c h hb b b bT T h b A l T S??? ??= = ( 3.12) Allowing the hot space temperatures of the two cylinders to increase, the results shown in Figure 3.8 are obtained. All four eigenvalues are purely imaginary; thus, the undamped engine cannot be made unstable. While the nonlinearity in the engine might drive the system away from equilibrium, any linear dissipation present in the system will make the motions decay in time. Figure 3.8 A plot of the undamped two cylinder double acting alpha system eigenvalues as the temperature is increased; all four system eigenvalues are purely imaginary and show no significant variation with respect to the heater temperature. The small real part can be attributed to numerical error. 44 Figure 3.9 A plot of the two cylinder double acting system eigenvalues as the damping is increased. Damping is applied to all pistons; the eigenvalues remain in the left-half plane for any loading. As expected, one can see from Figure 3.9 that linear dissipation drives the eigenvalues into the left half-plane; thus, a two cylinder alpha configuration is not promising. 3.1.1.5 Double Acting Alpha Engine: Three Cylinder Configuration Finally, the three cylinder alpha configuration is considered. In the case that all three cylinders have identical dimensions, the limiting behavior of the parameters becomes lim lim lim ln( / ) ( ) h h hT T Tk kk r h k h k k k h k h h Al Al Al T S Al VV V T T VAl AlT T T T T T T ? ?? ?? ?? = = = +? ? + + + +? ?? ? ? ( 3.13) lim lim 0 h hT T h Al T S??? ??= = ( 3.14) Increasing the temperature, the results shown in Figure 3.10 are obtained. 45 Figure 3.10 A plot of the eigenvalues of the undamped, three cylinder double acting alpha engine system as the hot side temperature is increased. The case considered is for three identical cylinders. There are two eigenvalues that migrate away from the imaginary axis with positive real part. There are, in fact, eigenvalues with positive real parts. Now applying damping to all three cylinders, the results shown in Figure 3.11 are obtained. 46 Figure 3.11 A plot of the eigenvalues of the three cylinder double acting alpha engine system as the damping is increased. The same damping is applied to all three elements. The right half-plane eigenvalues cross the imaginary axis, indicating a possible bifurcation point. A pair of complex conjugate eigenvalues can be made to cross the imaginary axis, while the remaining eigenvalues are in the left-half plane, which is the desired behavior. 3.2 Modal Coordinates, Method of Multiple Scales, and Numerical Results The displacement magnitudes of the nondimensionalized governing equations are always less than unity; as a consequence the system may be considered to be weakly nonlinear. This fact allows the application of perturbation analysis to the system; more specifically, the method of multiple scales will be used in order to determine the steady-state behavior of the system. However, in order to expedite the analysis, a change of coordinates is used in order to decouple the linear part of the system. Nonlinear dissipation is introduced into the system because it was found to be necessary to produce limit cycles in the FPSE. Cubic dissipation is used for simplicity; cubic damping can physically appear in the system in the form of a nonlinear feedback controller. Alternatively, viscous interactions in turbulent flow can also be modeled with a cubic function [6]; such dissipation can appear through the external load (a pump) or through the pressure drop across the heat exchangers. In general, so long as the nonlinear dissipation produces a secular term [18], it will serve the same purpose as the cubic damper, which is to produce an attracting limit cycle; for example, higher order odd nonlinearities will produce the same effect. The external load 47 serves as a means to physically introduce this nonlinearity into the system; for example, a linear (motion) alternator is a nonlinear load that is often used with the FPSE [6]. 3.2.1 Beta Engine 3.2.1.1 Modal Coordinates In order to express the governing equations in terms of the modal coordinates, the simplifying assumption of proportional damping is imposed. More specifically, the nondimensional damping on both the piston and displacer are assumed to be equal; this means that the damping matrix is a scalar multiple of the identity matrix. Then the eigenvalues and eigenvectors of the linear stiffness matrix are determined. Under these assumptions the equations become 1 2 1 2 1 ( , )1 0 1 0 0 ( , )0 1 0 1 0 p p p p p p d d d d d d HOT HOT ? ? ? ? ? ? ? ? ?? ? ? ? ? ? ? ? ? ? ? ? ? + ?? ? ? ? ? ?? ? ? ?? ? ? ? ? ? + + + =? ? ? ? ? ? ? ? ? ?? ?? ? ? ? ? ? ? ? ? ? ?? ? ? ? ? ?? ? ? ? ( 3.15) where HOT represents higher order terms. Solving the eigenvalue problem ?? ? ? ??? ?= ??? ? ??? ? ?? ? ?? ? ?+ ??+ 0 01 d p dd pp ? ? ?????? ????? ( 3.16) the following eigenvalues are obtained ( )?????????? pdiib +??== 4020202,1 mm ( 3.17) provided that ( ) 040 <+?? ??????? pd ( 3.18) where 210 ?????? dp ?++= ( 3.19) The eigenvectors are given by ?? ? ? ??? ? ?=?? ? ? ? ?? ? ? ? ?+ = ?? ?? ??? iv p p 11 1 2,12,1 ( 3.20) In order to transform the equations into the modal coordinates, the following change of variables is used: 48 [ ] ?? ? ? ??? ?= ??? ? ??? ? ?? ? ?? ? ?+=??? ? ??? ? 2 1 2 111 ? ? ? ? ????? ? V iid p ( 3.21) [ ] ? ? ? ?? ? ?? ??=? 1 1 2 1 ?? ?? ? i iiV ( 3.22) Left multiplying by the inverse of V, in the new coordinates, the equations become ??? ? ??? ?= ??? ? ??? ? ?? ? ?? ? ?? ??+ ??? ? ??? ?? ? ? ?? ? + ?+? ??? ? ??? ? ?? ? ?? ?+? ??? ? ??? ? ?? ? ?? ? 0 0 ),( ),( 1 1 2 0 0 10 01 10 01 21 21 2 1 2 0 2 0 2 1 2 1 ?? ?? ?? ?? ? ? ? ? ? ? ?? ? ? d p HOT HOT i ii ib ib ( 3.23) where the higher order terms have also been transformed. The advantage of this transformation becomes clear once the method of multiple scales is applied. 3.2.1.2 Method of Multiple Scales Analysis In the decoupled equations, the modal stiffness quantities are complex, with the imaginary part of the stiffness quantities representing energy transfer terms in the system. In order to investigate how these terms affect the amplitude growth, the method of multiple scales is applied [18]. To begin with, a bookkeeping parameter ? is introduced, and the damping, complex part of the stiffness, and the higher order terms are rescaled. When the damping parameter and imaginary part of the eigenvalues are close in magnitude, the system is near a bifurcation point; in other words, the eigenvalues of the state-space representation cross the imaginary axis in this region. Rescaling the difference between the damping and imaginary part of the eigenvalue corresponds to rescaling the real part of the eigenvalues in the state-space representation. This difference must be small in order to obtain a valid approximation. ??? ? ??? ?= ??? ? ??? ? ?? ? ?? ? ?? ??+ ??? ? ??? ?? ? ? ? ? ? + ?+? ??? ? ??? ? ?? ? ?? ?+? ??? ? ??? ? ?? ? ?? ? 0 0 ),( ),( 1 1 2 0 0 10 01 10 01 21 21 2 1 2 0 2 0 2 1 2 1 ?? ?? ?? ?? ?? ? ? ?? ?? ? ??? ? ? d p HOT HOT i ii bi bi ))) ( 3.24) The following time scales are defined, with ?n? being a nonnegative integer ?? nnT = ( 3.25) which turns the time differentiation into 49 0 1 d d T T?? ? ?? + + ? ? K ( 3.26) The solution is approximated as follows K+++? ),(),(),(),( 12211101 ?????????? tttt ( 3.27) K+++? ),(),(),(),( 22221202 ?????????? tttt ( 3.28) where ?1 and ?2 are referred to as the first and second modes, respectively. After substituting the above equations into the original problem, the following equations are obtained, where a subscript on the inverse of V indicates a row number ( ) ( ) ( ) [ ] 2 10 11 10 11 0 1 0 1 12 0 10 11 1 21( ) ( , ) 0p T T T T i b V HOT ? ? ?? ?? ? ? ?? ? ? ? ?? ? ? ?? ? ? ? ?? ? ? ?? ? ? ? ? ?+ + + + + + +? ?? ? ? ?? ?? ?? ? ? ? ? ? ? ?? ?? ? + ? + + + = )K K ) K ( 3.29) ( ) ( ) ( ) [ ] 2 20 21 20 21 0 1 0 1 12 0 12 12 1 22( ) ( , ) 0d T T T T i b V HOT ? ? ?? ?? ? ? ?? ? ? ? ?? ? ? ?? ? ? ? ?? ? ? ?? ? ? ? ? ?+ + + + + + +? ?? ? ? ?? ?? ?? ? ? ? ? ? ? ?? ?? ? + + + + + = )K K ) K ( 3.30) Grouping terms with like powers of ?, one obtains the following hierarchy of equations ( ) 0?O ( 3.31) 2 2 10 0 102 0 0T ? ? ?? + =? 2 2 20 0 202 0 0T ? ? ?? + =? The general solution to these equations in the complex domain is 0 0 0 0 0 0 0 01 1 2 1( ) ( )10 0 1 1 1 2 1 1 1 2 11 1( , ) ( ) ( ) ( ) ( )2 2i T i T i T i Ti T i TT T A T e A T e a T e e a T e e? ? ? ?? ?? ? ?= + = + ( 3.32) 0 0 0 0 3 1 0 0 0 04 1( ) ( )20 0 1 3 1 4 1 3 1 4 11 1( , ) ( ) ( ) ( ) ( )2 2i T i T i T i T i Ti TT T A T e A T e a T e e a T e e? ? ? ? ??? ? ?= + = + ( 3.33) Not considering the higher order terms for now, one has 50 ( )1?O ( 3.34) 22 2 10 10 11 0 11 102 0 1 0 0 2 ibT T T T? ?? ? ? ? ?? ?? + = ? ? +? ? ? ? )) 22 2 20 20 21 0 21 202 0 1 0 0 2 ibT T T T? ?? ? ? ? ?? ?? + = ? ? ?? ? ? ? )) Grouping the sources of secular terms with positive frequency, and considering solutions of positive and negative frequencies separately, the following complex modulation equations are obtained for positive frequencies 0 1 0 1 12 0i A i A ibA? ? ??? ? + = ( 3.35) 0 3 0 3 32 0i A i A ibA? ? ??? ? + = ( 3.36) and for negative frequencies, one has 0 2 0 2 22 0i A i A ibA? ? ?? + + = ( 3.37) 0 4 0 4 42 0i A i A ibA? ? ?? + + = ( 3.38) Applying the polar form of the slowly varying amplitudes and separating the real and imaginary parts, the amplitude and phase of the positive frequency solutions becomes 1 11 02 2 ba aa ? ?? = ? ( 3.39) 1 0??= ( 3.40) 3 33 02 2 ba aa ? ?? = ? ? ( 3.41) 3 0?? = ( 3.42) and those for the negative frequencies are obtained as 2 22 02 2 ba aa ? ?? = ? ? ( 3.43) 2 0?? = ( 3.44) 4 44 02 2 ba aa ? ?? = ? ( 3.45) 4 0?? = ( 3.46) The amplitude of the first mode dies out for negative frequencies while the amplitude of the second mode dies out for positive frequencies. Based on this fact, a simplifying assumption is made; the steady-state solutions will contain positive frequency solutions for the first mode and negative frequency solutions for the second mode. Furthermore, since the coefficients of the 51 equations governing modes 1 and 2 are complex conjugates of one another; the solutions are assumed to be complex conjugates of one another. Thus, the steady-state solution is assumed to have the form 0 0 0 01( )1 0 1 1 11( , ) ( ) ( )2j T j Ti TT T A T e a T e e? ??? = = ( 3.47) 0 0 0 01( )2 0 1 1 11( , ) ( ) ( )2j T j Ti TT T A T e a T e e? ??? ? ??= = ( 3.48) Next, the effect of nonlinear damping is considered; nonlinear dissipation is capable of canceling the growth of the solutions for large amplitudes, resulting in a limit cycle. For example, applying a cubic damper to the piston, the following contributions to the governing equations in modal form are obtained: [ ] [ ] ( ) ( ) 3 31 1 1 2 0 1 3 1 2 0 1 0 0 2 pV V T T ii i T T ?? ? ? ?? ? ??? ? ? ? ? ?? ? ? ? ?? ?? ?? ? ? ?? ? ? + +? ?? ?= = ? ?? ?? ? ? ?? ?? ? ? ? ? ?? ? ? ? ? ?? ? ?? ? ? ?+ + ? ?? ?? ?? ? ? ?? ?? ? ? ? ? ? ) ( 3.49) Considering either mode 1 or mode 2, this will contribute a secular term to the amplitude and phase equations to yield: 2 20 0 ? ? ?3 2 2 16 ba a a???? ? ?? ??= ? ?? ?? ? ? ?? ?? ?? ? ( 3.50) 3 30 0?3 016a a??? ? ???+ = ( 3.51) A nonzero equilibrium solution is given by ( )03 03 8 ba ?? ?? ?= ( 3.52) 2 2 2 20 1 0 0 0?3 316 16a T a?? ??? ? ? ? ? ?? ?= + = + ( 3.53) For the original displacements, one has: 52 [ ] [ ] 0 01 11 2 1 1 0 02 1 cos(( ) ) cos(( ) ) p d av v v v ra ? ? ? ? ?? ? ? ? ? ? ? ?? ? + +? ? ? ?? ? ? ?= = = ? ? ? ?? ? ? ? + + +? ? ? ?? ? ? ? ( 3.54) where 2 2r ? ?= + ( 3.55) 1tan ?? ?? ? ?= ? ? ? ? ( 3.56) and the value of ?0 will depend on initial conditions and possibly transients. 3.2.1.3 Numerical Results In this section, the analytical results from the previous section are compared to numerically computed solutions of the nonlinear equations governing the Beta engine; these equations are derived based on the Schmidt model. All of the analysis is carried out in terms of the nondimensional equations; several numerical examples are presented in this section. 3.2.1.3.1 Case 1: No Frequency Correction As mentioned, the approximation is valid near the bifurcation point. When the product of the frequency and damping is larger than the imaginary part of the complex eigenvalues, the solutions decay. However, when the situation is reversed, a limit cycle is observed. The following parameter values are used in the numerical example: 20p? = 13d? = .4? = .3? = 9.1? = 0.1? = The associated eigenvalues and eigenvectors are 1,2 9 0.894i? = m 1,2 10.149v i? ?= ? ?? ? ? First, consider 01 / 0.298b? ?= > = . The corresponding results are shown in Figure 3.12. 53 Figure 3.12 Numerical solution to the equations of motion when the damping is too high. As predicted, the system exhibits decaying motions. As predicted, the solutions decay. Next, consider 0.281 /b? ?= < . The corresponding results are shown in Figure 3.13 and Figure 3.14. Figure 3.13 Numerical solution of the steady-state motion near the bifurcation point; the piston motion is described by the sinusoid with large amplitude whereas the displacer has the smaller amplitude. Figure 3.14 A parametric plot of the steady-state piston and displacer displacements. As predicted, in the limit, the system migrates to a limit cycle. The analytical approximation (inner ellipse) is shown for comparison with the numerical solution (outer ellipse). 54 Table 3.1 Comparison of analytical and numerical results; ?=0.281 Analytical Approximation Numerical Integration Angular Frequency 3 3 Piston Amplitude 0.22 0.22 Amplitude Ratio 0.149 0.16 Displacer Phase Angle 90? 95? However, the approximation breaks down far from the bifurcation point. Also, one can observe that as the amplitudes become larger, a bias can be seen, likely caused by neglected quadratic terms. For instance, for 0.15 /b? ?= < , the corresponding results are as shown in Figure 3.15 and Figure 3.16. Figure 3.15 The system still possesses a limit cycle far from the bifurcation point. Figure 3.16 A parametric plot of the piston and displacer displacements far from the bifurcation point. The analytical approximation (inner ellipse) is not valid here. Table 3.2 Comparison of analytical and numerical results; ?=0.15 Analytical Approximation Numerical Integration Angular Frequency 3 2.73 Piston Amplitude .66 .8 Amplitude Ratio 0.149 .375 Displacer Phase Angle 90? 95? 55 Figure 3.17 Parametric plot of solution curves for different initial conditions when the system has purely imaginary eigenvalues and no nonlinear damping. As discussed in [5] and [19], the origin is a center, about which there is a continuum of periodic orbits. It is also worth noting that at the bifurcation point, if the cubic damper is removed, the result is that the system has a continuum of periodic orbits around the origin. This situation is analogous to the linear analysis presented in references [5] and [19]. Setting 0/b? ?= and varying the initial conditions shows that the origin becomes a center. The corresponding results are given in Figure 3.17. 3.2.1.3.2 Case 2: Frequency Correction Now, a second numerical example in which the frequency correction is nonzero is constructed. The analysis is again carried out in terms of the nondimensional equations. The following parameter values are used in the numerical example: 20p? = ( 3.57) 10d? = ( 3.58) .4? = ( 3.59) .3? = ( 3.60) 8? = ( 3.61) 0.1? = ( 3.62) This results in the following eigenvalues and eigenvectors: 1,2 7 4.472i? = m ( 3.63) 1,2 10.333 0.745v i? ?= ? ?? ? ? 56 Figure 3.18 Numerical solution of the steady-state motion near the bifurcation point; the piston motion is described by the sinusoid with large amplitude whereas the displacer has the smaller amplitude. Figure 3.19 A parametric plot of the steady-state piston and displacer displacements. As predicted, in the limit, the system migrates to a limit cycle. The analytical approximation (inner elipse) is shown for comparison with the numerical solution (outer ellipse). Considering 01.6 / 1.69b? ?= < = , the results shown in Figure 3.18 and Figure 3.19 are obtained. 3.2.1.3.3 Dimensional Parameters for a Prototype Engine Although an analytical criterion has been developed to engineer Hopf bifurcations in Beta configuration engines, the results are in the form of nondimensional parameters. In order to test the feasibility of obtaining such parameters, dimensional quantities of an engine are specified here. More specifically, the parameters of the Sunpower RE-1000 are modified in order to satisfy the following set of nondimensional relationships. 57 The parameters of the Sunpower RE-1000 as reported in [5] are given by 2 2 2 7.1 814.3 322.8 4.073 * / 25.69 2.176 18.61 18.30 6.20 .426 295.7 / 14.76 / 461.5 / 215.7 / m h k p d e c p d p d l d P Mpa T K T K S cm mm K A cm A cm l mm l mm m kg m kg k kN m k kN m c Ns m c Ns m = = = = = = = = = = = = = = In order to remain in the neighborhood of feasible values, the heat exchanger parameters are not varied; in addition, the displacer damping, which is largely due to viscous losses, is not altered. The piston mass, displacer mass, piston spring stiffness, displacer spring stiffness, load, and cubic damping are varied in order to satisfy the necessary criterion for a limit cycle, which is 0 0b ??? > . By using Mathematica, the control parameters are varied to satisfy the necessary condition with the following parameter modifications: 3 3 5.33 3.2 20 / 244 / 359 / 2.6 / p d p d p cubic m kg m kg k kN m k kN m c Ns m c Ns m = = = = = = The nondimensional parameters corresponding to this choice of dimensional parameters are 49.7p? = 7.03d? = 0.358? = 0.186? = 58 Figure 3.20 A parametric plot of the piston and displacer displacements for the parameter values specified in this section. The outer ellipse is the analytical approximation, whereas the inner ellipse is the numerical solution. 20.26? = 1.1? = 0.01? = where the nondimensional and dimensional cubic dampers are related as follows 2 p p cubic p p k l c m m? = ( 3.64) With this choice of parameters, the steady-state trajectory shown in Figure 3.20 is obtained. The parameter modifications required to obtain the limit cycle shown in Figure 3.20 increase the system mass by 1.9 kg, reduce the load coefficient by 22%, and introduce a very small cubic damper into the system. Lastly, the parameter modifications very nearly switch the displacer spring stiffness and the piston spring stiffness. All of these changes are, in fact, physically realizable; indicating that a machine designed for a limit cycle motion is also physically realizable. 3.2.2 Double Acting Alpha Engine 3.2.2.1 Modal Coordinates In order to express the governing equations in terms of the modal coordinates, the simplifying assumption of proportional damping is applied. More specifically, the nondimensional damping associated with all three pistons is equal; this means that the damping matrix is a scalar multiple of the identity matrix. Proceeding to determine the eigenvalues and eigenvectors of the linear stiffness matrix, the following relations are obtained: 59 1 1 1 2 2 2 3 3 3 1 2 3 1 0 0 1 0 0 1 ( ) 0 1 0 0 1 0 1 ( ) 0 0 1 0 0 1 1 ( ) 0 0 0 HOT HOT HOT ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? + + ? ?? ?? ? ? ?? ? ? ?? ? ? ? ? ? ? ?? ? ? ? ? ?+ + ? + + ? ? ? ? ? ? ?? ? ? ? ? ?? ? ? ? ? ? ? ? ? ? ? ?? ? + +? ?? ? ? ?? ? ? ?? ? ? ? ? ? ? ? ? ?+ = ? ? ? ?? ? ? ? ? ? ? ? ( 3.65) Solving the eigenvalue problem ? ? ? ? ? ? ? ? ? ? = ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?++?? ??++? ???++ 0 0 0 )(1 )(1 )(1 3 2 1 ? ? ? ?????? ?????? ?????? ( 3.66) the following eigenvalues are obtained: 1 1 ( )( 1)? ? ? ?= + + ? ( 3.67) 2,3 2 ( )(2 1) 3( )2 i? ? ? ? ?? + + + ?= m ( 3.68) where 22 1? ?= ( 3.69) 21 2,3Re( )? ?= ( 3.70) 2,3Im( )b ?= ( 3.71) The eigenvectors are given as follows [ ] [ ]1 2 3 1 3 1 31 2 2 1 3 1 31 2 2 1 1 1 i i i iV v v v ? ?? ? ? + ? ? ? ? ? ?? + ? ?= = ? ? ? ? ? ? ? ?? ? ? ? ( 3.72) The inverse of [V] is 60 [ ] 1 1 1 1 1 1 3 1 3 1 3 2 2 1 3 1 3 1 2 2 i iV i i ? ? ? ? ? ? ?? + ? ? ? ?= ? ? ? ?? ? ? + ? ? ? ?? ? ( 3.73) Based on these results, if the three cylinders are identical, the eigenvectors do not depend on the parameters. Furthermore, the complex entries in V have a phase of ?120?, indicating symmetry in the motion. More specifically, if all three cylinders are identical, the three oscillating elements are always one third of a cycle apart from one another. In order to transform the equations to the modal coordinates, the following change of variables is used [ ] ? ? ? ? ? ? ? ? ? ? = ? ? ? ? ? ? ? ? ? ? 3 2 1 321 3 2 1 ? ? ? ? ? ? vvv ( 3.74) and the equations are left multiplied by the inverse of [V]. In the new coordinates, the equations become 2 1 1 2 1 2 2 2 1 2 2 3 3 1 3 1 2 3 1 0 0 1 0 0 0 0 0 1 0 0 1 0 0 0 0 0 1 0 0 1 0 0 1 1 1 0 1 1 3 1 3 1 0 3 2 2 0 1 3 1 3 1 2 2 ib ib HOTi i HOT HOTi i ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?? ?? ? ? ?? ? ? ? ? ?? ? ? ? ? ?? ? ? ?+ + ? ? ?? ? ? ? ? ?? ? ? ?? ? ? ? ? ? ? ?? ? ? ? +? ?? ? ? ?? ? ? ?? ? ? ? ? ? ? ? ?? ?? + ? ? ? ?? ?+ = ? ?? ?? ? ? ?? ? ?? ? ? + ? ? ? ?? ? ? ? ? ? ?? ? ? ( 3.75) where the higher order terms have also been transformed. 3.2.2.2 Simplifying Approximations No analytical approximation was developed to describe limit cycles in the three cylinder double acting FPSE. However, based on the form of the equations, some simplifying assumptions are made. First, since there is no imaginary part in the first modal stiffness, the first mode is assumed to decay. This assumption is based on the behavior of the beta configuration engine, where the imaginary part of the modal stiffness causes amplitude growth. 61 Second, since the coefficients of the equations governing modes 2 and 3 are complex conjugates of one another, the solutions to these equations are also assumed to be complex conjugates of one another. Finally, the higher order terms in the pressure are neglected. These assumptions are stated as 1lim 0 t ? ?? = ( 3.76) 2 3? ?= ( 3.77) Based on this, the following substitution can be made to determine the real and imaginary parts of modes 2 and 3: 2 32,3Re( ) 2? ?? ? += = ( 3.78) 2 32,3Im( ) 2i? ?? ? ?= = ( 3.79) If cubic damping is added and the above assumptions are applied, the modal equations reduce to 2 2 2 20 2 3 3 30 2 2 3 2 2 3 1 0 1 0 0 0 1 0 1 0 03 0 ib ib ? ? ??? ? ? ?? ? ?? ? ? ? ? ? ??? ? ? ? ? ?? ? ? ? + + ? ?? ? ? ? ? ?? ? ? ? + ? ? ? ?? ? ? ? ? ?? ? ? ?? ? ? ?+ = ? ? ? ?? ? ? ?? ? ( 3.80) First adding and then subtracting the equations for mode 2 and 3, then applying the substitution, the following equations are obtained governing the real and imaginary parts of mode 2 and 3 2 2 21 3 ( ) 0b? ?? ? ? ? ? ? ? ??? ? ? ? ?+ + + + + = ( 3.81) 2 2 21 3 ( ) 0b? ?? ? ? ? ? ? ? ??? ? ? ? ?+ + ? + + = ( 3.82) An analytical approximation can be used to solve these equations and the displacements in the physical coordinates can be determined by reversing the coordinate transformations. 3.2.2.3 Numerical Results If the following parameter values are chosen, then an attracting limit cycle is observed in the three cylinder FPSE, as shown in Figure 3.21. 0.3? = 5? = 0.4? = 0.2? = 62 (a) (b) Figure 3.21 (a) A plot of the steady-state displacements of the three pistons with respect to time; all three are identical except for a phase shift. (b) A three-dimensional parametric plot of the steady-state displacements. 0.1? = As predicted, the motions of the three pistons are one third of a cycle apart. Next, the assumptions applied to the solutions are investigated. 63 Figure 3.22 A plot of mode one with respect to time. The maximum value is negligible; thus, the approximation that mode one decays with increasing time is valid. (a) (b) Figure 3.23 Plots of the real part of (a) mode two and (b) mode three with respect to time; they are the same. 64 (a) (b) Figure 3.24 Plots of the imaginary part of (a) mode two and (b) mode three with respect to time; they are of opposite sign. As shown in Figure 3.22, mode 1 is observed to have a negligible contribution to the motion. Based on the results shown in Figure 3.23 and Figure 3.24 modes 2 and 3 are observed to be complex conjugates of one another, as predicted. Finally, the largest quadratic term in the expansion of the pressure is 2 2 .0256i? ? = , which is only 6.4% of the amplitude. Thus, the simplifying assumptions are valid for the numerical case considered. 3.3 Numerical Results from Simplified Nodal Analysis The models obtained on the basis of simplified nodal methods proved to be unstable. In other words, the numerical solutions would often exhibit growing characteristics. Furthermore, the solution to these equations required stiff solvers and long computation times for short time intervals. These problems are attributed to invalid simplifying assumptions. For instance, the quasi-steady flow assumption is likely not satisfied. Indeed, previous studies have shown that this assumption is not always valid. Another assumption that may cause these problems is assuming that only the expansion and compression spaces have time varying masses. Finally, the method used to define flux quantities is also of importance. More specifically, the properties 65 assigned to the fluid flowing in or out of a control volume can be specified in more than one way. Some of the results are included as well as the programs used to obtain them; however, the effort to create a reduced-order model of FPSEs that captures more thermodynamic complexity than the Schmidt model is an area that requires more work. All dimensional quantities are specified in fundamental SI units. 3.3.1 Single Cylinder Case Applying an isothermal assumption to the single cylinder equations, the plots shown in Figure 3.25 are obtained. (a) (b) 66 (c) Figure 3.25 Plots of various quantities versus time obtained from an isothermal nodal analysis of the single cylinder device: (a) piston displacement, (b) working space pressure, and (c) mass of working fluid. The numerical solution tends to grow in time; the root cause of this is attributed to inconsistent assumptions. Of particular interest is the net buildup of fluid mass within the system, which is a clear indicator that the solution is not physically realizable. In the next simulation, the isothermal assumption was removed. Furthermore, periodic heat was applied to the cylinder such that heat addition would occur during an expansion process; a conduction path was also allowed from within the cylinder to the outside air. If this conduction path was not included, then the solution would again grow in time. This additional complexity within the model required the use of stiff solvers, and depending on the parameter values chosen, the computation time was drastically increased. The plots obtained in one case are shown in Figure 3.26. (a) 67 (b) (c) (d) Figure 3.26 Plots of various quantities versus time obtained from a non-isothermal nodal analysis of the single cylinder device: (a) piston displacement, (b) working space pressure, (c) fluid temperature, and (d) mass of working fluid. Although the system properties oscillate in time, the thermodynamic properties do not vary smoothly in time; much of this can be attributed to the discontinuous heating applied to the cylinder. 68 3.3.2 Beta Engine The simplified nodal analysis was also applied to the Beta configuration FPSE. Similar problems plagued these simulations, despite specifying the parameters based on a working engine. The parameter values for the Sunpower RE-1000 were obtained from [12] and [5]. 3.3.2.1 Isothermal Case The isothermal equations again had growing solutions. Furthermore, increasing the load coefficients did not remedy the situation. A stiff solver was used to obtain the solution; but the solutions were found to grow to nonphysical values in a fraction of a second; the results are shown in Figure 3.27. (a) (b) 69 (c) (d) (e) Figure 3.27 Plots of various quantities versus time obtained from an isothermal nodal analysis of the Beta configuration Sunpower RE-1000 engine: (a) piston displacement, (b) displacer displacement, (c) expansion space pressure, and (d) pressure differential. 70 3.3.2.2 Non-Isothermal Case In contrast to the previous isothermal simulation, the non-isothermal simulation results tended to exhibit characteristics of decaying motions. The parameters were again chosen from [12] and [5]; the results are shown in Figure 3.28. (a) (b) (c) 71 (d) Figure 3.28 Plots of various quantities versus time obtained from a non-isothermal nodal analysis of the Beta configuration Sunpower RE-1000 engine: (a) piston displacement, (b) displacer displacement, (c) compression space pressure, and (d) mass of working fluid in working space volumes. In this chapter, analytical and numerical results have been presented to illustrate the possibility of engineering a Hopf bifurcation of an equilibrium solution and ensuring limit cycle motions in a FPSE. The obtained numerical results suggest that the Hopf bifurcations are supercritical in nature. 72 4 Concluding Remarks and Suggestions for Future Work In this chapter, the work done in this thesis is summarized and the significance of the work in relation to FPSE research as a whole is explained. Subsequently, areas of future work are identified. 4.1 Limit Cycle Motions This study has applied a strategy towards engineering a Hopf bifurcation in a weakly nonlinear dynamical system. Specifically, a plot of the system eigenvalues is obtained as a control parameter is varied; if the transversality condition [17] is met, then introduction of the appropriate nonlinearity into the system can lead to a Hopf bifurcation, leading to an attracting limit cycle. An appropriate simplification method can be used to study the weakly nonlinear system analytically, such as the Method of Multiple Scales or Center Manifold Reduction. This method was applied to develop an analytical basis for limit cycle motion in beta FPSEs with cubic damping. The Schmidt model was used to obtain the governing equations of motion for the FPSE. After nondimensionalizing the equations of motion, they were observed to be weakly nonlinear. Then, cubic damping was added to the system and the Method of Multiple scales was used to develop an analytical basis to obtain attracting limit cycles in beta FPSEs. The system was assumed to be proportionally damped. By using this assumption, the system dynamics could be described by the eigenvalues of the stiffness matrix, a damping parameter, and a cubic damping parameter. Since the stiffness matrix is a non-symmetric matrix, it is possible for it to have complex eigenvalues, which occur in complex conjugate pairs. It is the imaginary part of these eigenvalues which were found to cause growing motions in the system; the linear dissipation competed with this effect. If the system is configured to have growing motions from the static equilibrium, then the addition of a cubic damper limits the growth, leading to an attracting limit cycle. The obtained analytical relation was compared to the numerical results in order to evaluate its validity. The analytical and numerical solutions are found to be close to one another near the bifurcation point; although the approximation breaks down far from the bifurcation point, at the bifurcation point itself, and for large enough amplitudes. In addition, a set of dimensional parameters for a feasible engine design incorporating these results is provided; the parameter selection was guided by the parameters reported in the literature for the Sunpower RE-1000 engine. It is also worth noting that although this study considered cubic damping, there are other types of damping that can lead to attracting limit cycles; the method presented here can be applied to other forms of nonlinear damping as well. This result demonstrates the advantage of using nonlinearity in the design of FPSEs because such an approach can yield a unique attracting periodic trajectory. A linear analysis yields periodic motions whose amplitudes are sensitive to disturbances. Thus, the nonlinearity can be used to stabilize the motion of the device. 73 In addition, attracting limit cycles were also observed in double acting three cylinder FPSEs by using numerical integration. However, no analytical result was obtained describing the limit cycle. Even so, a foundation has been laid to find a similar analytical basis for limit cycles in three cylinder double acting engines. Based on the form of the equations, simplifying assumptions were made, and then verified using numerical integration. The simplifying assumptions lead to a system of two real valued second order differential equtations; these equations can be used to develop an analytical relation to describe limit cycles in double acting engines. Such an analysis could be generalized to double acting FPSEs with even more cylinders. 4.2 Reduced Order Models The simplified nodal methods were an attempt at developing a reduced order model of the FPSE that captured more of the thermodynamic behavior than the Schmidt model, but at the same time, still made dynamic analysis possible. Specifically, the nodal discretization schemes that have been applied to Stirling engine analysis were simplified to replace all of the heat exchanger nodes with algebraic expressions. Thus, the simplified nodal methods made use of two deforming control volumes and fluid flow correlations. The simplified nodal methods proved to be difficult to integrate numerically. Even stiff solvers would take a long time to solve relatively small time intervals; furthermore, over long time periods, the solvers would break down. Furthermore, the solutions were, in many cases, not physical. Thus, the simplified nodal methods as presented in this study did not yield physically valid models. There are several factors which can contribute to the failure of these models. First, there are a number of simplifying assumptions applied to the fluid flow that are not necessarily true. More specifically, the flow is assumed to be quasi-steady. Thus, the fluid is modeled with steady state pipe flow correlation; this approximation is only valid under certain conditions [15]. Then, there is the possibility of numerical instability [20]. Depending on the types of finite differences used, the stability of the solution is affected. Furthermore, the stability of the solution depends on how the properties of the fluxes are defined. Any one of the above factors, or several of them, can be at work in causing the simplified nodal methods to fail. 4.3 Future Work This study has shown that limit cycles are possible in FPSEs when nonlinear damping is present; however, there is still the issue of how this nonlinear dissipation comes into the system. On the one hand, the viscous losses associated with the working fluid introduce nonlinear dissipation if the flow is turbulent; however, it needs to be described in a manner such that the pressure gradient can be expressed in terms of the motion of the engine components. Alternatively, it is possible to introduce a nonlinear load. A comprehensive study on the various nonlinearities that can be introduced into an FPSE to cause limit cycle motion would be of value. Another property of the limit cycle that was not analytically examined is its domain of attraction. If one considers the modulation equations associated with the slowly varying amplitudes as a dynamical system, then an appropriate Lyapunov function associated with this system can be used to estimate the domain of attraction for the limit cycle [17]. 74 In addition, there are applications in which the FPSE would be subject to disturbances in the environment. For example, if it were used on a rover for space applications, as the rover moved over a rough planetary surface, the attached FPSE would experience a time varying external force. A study of how the limit cycle is affected by external disturbances would prove valuable for such applications. Furthermore, this study neglected any direct effects of viscous dissipation within the system; rather, linear and cubic dampers were used to capture the dissipation in the system. An experimentally validated model of how the viscous losses relate to the system dynamics would facilitate dynamic analysis of FPSEs. Another area of further research is to develop an analytical basis for limit cycles in double acting engines, for three cylinders or possibly more. Such a relation could prove valuable in the design of double acting FPSEs. Finally, although the simplified nodal methods did not provide valid solutions, the underlying goal of those methods is still of interest. That is, to create a reduced-order model of FPSEs that, while emphasizing the dynamics of the engine, provides greater accuracy of the fluid thermodynamics than the Schmidt model. Further work using the nodal methods should incorporate the Riemann problem in fluid mechanics [21]. This problem serves as a prototype for flow problems in which the boundary data contains a discontinuity. In the case of the Stirling engine, the discontinuity comes from the flux quantities because the flow oscillates. The Riemann problem provides information about how flux quantities should be specified and the maximum allowable grid size required for convergence of the numerical scheme. Appendix A: Computer Programs In this appendix, the different Matlab and Mathematica programs used to generate the root locus plots and carry out the numerical integrations are presented. Root Locii Root Locus for Beta Engine %The eigenvalues of the beta engine are plotted as the hot side temperature %or the damping are varied. The equations are in nondimensional form. clc clear for i=1:201 alphap = 2; alphad=1; beta = .4; kapa=1.7; %The dimensionless hot side temperature is varied. gamma(i)=.4-.4/(((i-1)*.1+1)); c(i)=0; A=[0, 0, 1, 0 ; 0, 0, 0, 1;-(1+alphap*beta),alphap*gamma(i), -c(i), 0; -alphad*beta, -(kapa-alphad*gamma(i)), 0,0]; p=poly(A); r=roots(p); r plot(r,'bx') title('Undamped Beta') xlabel('Re[s]') ylabel('Im[s]') hold on clear end figure for i=1:201 alphap = 2; alphad=1; beta = .4; kapa=1.7; %Damping is applied to only the piston. gamma(i)=.4-.4/(((201-1)*.1+1)); c(i)=(i-1)*.02; A=[0, 0, 1, 0 ; 0, 0, 0, 1;-(1+alphap*beta),alphap*gamma(i), -c(i), 0; -alphad*beta, -(kapa-alphad*gamma(i)), 0,0]; p=poly(A); r=roots(p); r plot(r,'bx') title('Beta: Damped Piston') xlabel('Re[s]') ylabel('Im[s]') hold on clear end figure for i=1:201 alphap = 2; alphad=1; beta = .4; kapa=1.7; %Damping is applied to both the piston and displacer. gamma(i)=.4-.4/(((201-1)*.1+1)); c(i)=(i-1)*.02; A=[0, 0, 1, 0 ; 0, 0, 0, 1;-(1+alphap*beta),alphap*gamma(i), -c(i), 0; -alphad*beta, -(kapa-alphad*gamma(i)), 0,-c(i)]; p=poly(A); r=roots(p); r plot(r,'bx') title('Beta: Damped Piston and Displacer') xlabel('Re[s]') ylabel('Im[s]') hold on clear end 76 Root Locus for Double Acting Alpha Engine: Three Cylinder Case %Plot the eigenvalues of the %equations of motion as the temperature or loading are varied. The %equations are in nondimensional form. clc clear alpha = 3, beta = .4 for i=1:201 alpha=3; beta=.4; %The dimensionless hot side temperature is varied. gamma(i)=1/((i-1)*.2+1); c(i)=0; A=[0, 0, 0, 1, 0, 0 ; 0, 0, 0, 0, 1, 0; 0, 0, 0, 0, 0, 1;-1-alpha*(beta+gamma(i)),gamma(i), beta, -c(i), 0, 0; beta, -1- alpha*(beta+gamma(i)), gamma(i), 0, -c(i), 0; gamma(i), beta,-1-alpha*(beta+gamma(i)), 0, 0, -c(i)] p=poly(A); r=roots(p); r plot(r,'bx') title('Undamped Three Cylinder Alpha') xlabel('Re[s]') ylabel('Im[s]') hold on clear end figure for i=1:201 alpha=3; beta=.4; %The damping is varied. gamma(i)=1/((201-1)*.2+1); c(i)=(i-1)*.01; A=[0, 0, 0, 1, 0, 0 ; 0, 0, 0, 0, 1, 0; 0, 0, 0, 0, 0, 1;-1-alpha*(beta+gamma(i)),gamma(i), beta, -c(i), 0, 0; beta, -1- alpha*(beta+gamma(i)), gamma(i), 0, -c(i), 0; gamma(i), beta,-1-alpha*(beta+gamma(i)), 0, 0, -c(i)] p=poly(A); r=roots(p); r plot(r,'bx') title('Damped Three Cylinder Alpha') xlabel('Re[s]') ylabel('Im[s]') hold on clear end 77 Numerical Integration of Schmidt Model Beta Configuration Free Piston Stirling Engine This program is used to numerically solve the equations of motion for the Beta configuration Stirling engine. The equations of motion are based on the Schmidt model and are in nondimensional form. The numerical solution is also compared to the analytical approximation generated by using the method of multiple scales. The first cell is used to define the dimensionless parameters in the equations of motion. alphap=20 alphad=1/3 beta=.4 gamma=.3 kappa=9.1 c=.281 mu=.1 The dimensionless stiffness matrix is defined here; information required for the analytical approximation is extracted. The analytical approximation is defined. K={{1+alphap*beta, -alphap*gamma},{alphad*beta, kappa-alphad*gamma}} lambda=Eigenvalues[K] V0=Eigenvectors[K] V=V0/Part[V0,1,1] w=Sqrt[Re[Part[lambda,1]]] b=Abs[Im[Part[lambda,1]]] kisei=Re[Part[V,1,2]] psi=Abs[Im[Part[V,1,2]]] a=Sqrt[(b-c*w)/(3*mu*w^3/8)] beta1=3*mu*kisei*w^2*a^2/(16*psi) r=Sqrt[kisei^2+psi^2] theta=If[kisei==0,ArcTan[Infinity], ArcTan[psi/kisei]] Axp[t_]:=a*Cos[(w+beta1)*t] Axd[t_]:=r*a*Cos[(w+beta1)*t+theta] The governing equations are numerically solved. In order to verify that the limit cycle is attracting, a second solution with varying intial conditions is also obtained. q=NDSolve[{xp''[t]+c*xp'[t]+mu*xp'[t]^3+xp[t]+alphap((1-beta*xp[t]+gamma*xd[t])^-1- 1)==0,xd''[t]+c*xd'[t]+kappa*xd[t]+alphad((1-beta*xp[t]+gamma*xd[t])^-1-1)==0, xp[0]==1, xd[0]==0, xp'[0]==0, xd'[0]==0},{xp,xd},{t,0,300}, MaxSteps->10^6] qq=Table[NDSolve[{xp''[t]+c*xp'[t]+mu*xp'[t]^3+xp[t]+alphap((1-beta*xp[t]+gamma*xd[t])^- 1-1)==0,xd''[t]+c*xd'[t]+kappa*xd[t]+alphad((1-beta*xp[t]+gamma*xd[t])^-1-1)==0, xp[0]==i, xd[0]==0, xp'[0]==0, xd'[0]==0},{xp,xd},{t,0,300}, MaxSteps->10^6],{i,0,1,.15}] 78 Various plots of transient and steady state behavior are generated. Also, the function values over a time interval are placed in an array in order to determine the period of motion and phase angle between the displacer and piston responses. Plot[Evaluate[{xp[t],xd[t]}/.q],{t,0,20},PlotRange->All, AxesLabel->{\[Tau],Subscript[x,i]}] Plot[Evaluate[{xp[t],xd[t]}/.q],{t,100,120},PlotRange->All, AxesLabel- >{\[Tau],Subscript[x,i]}] Plot[{Axp[t],Axd[t]},{t,100,120},PlotRange->All, AxesLabel->{\[Tau],Subscript[x,i]}] ParametricPlot[{Evaluate[{{xp[t],xd[t]}/.q}]},{t,0,30},PlotRange->All,AxesLabel - >{Subscript[x,p],Subscript[x,d]}] ParametricPlot[{{Evaluate[{{xp[t],xd[t]}/.q}]},{Axp[t],Axd[t]}},{t,260,280},PlotRange->All, AxesLabel ->{Subscript[x,p],Subscript[x,d]}] ParametricPlot[{Evaluate[{{xp[t],xd[t]}/.qq}]},{t,280,300},PlotRange->All, AxesLabel - >{Subscript[x,p],Subscript[x,d]}] Table[Evaluate[xp[i]/.q],{i,105,109,.05}] Table[Evaluate[xd[i]/.q],{i,105,109,.05}] 79 Parameter Selection The Sunpower RE-1000 engine is modified so that the system produces an attracting limit cycle. The displacer damping is a reflection of the viscous dissipation that occurs in the heat exchangers; thus, it is not changed. Also, the heat exchanger parameters are not changed. These parameters are defined first. cd=215.7 beta=.358 gamma=.1855 The remaining parameters are varied in order to determine where the quantity b-c*w>0, which is necessary for an attracting limit cycle to appear in the system. Recall that b is the imaginary part of the stiffness matrix eigenvalue, c is the nondimensional daming, and w is the square root of the real part of the stiffness matrix eigenvalue. Manipulate[{mp=cp*md/cd,k=cp*kd/(cd*kp), c=cp/(Sqrt[mp*kp]),ad=cp*.01189*Pm/(cd*kp),ap=.140*Pm/kp, K={{1+ap*beta, -ap*gamma},{ad*beta,k-ad*gamma}}, s=Eigenvalues[K], b=Abs[Im[Part[s,1]]],w=Sqrt[Abs[Re[Part[s,1]]]]},{kd,200000,500000},{cp,100,1000}, {Pm,1000000,10000000},{md,.3,20}] Manipulate[Plot[b-c*w,{kp,10000,500000}, AxesLabel-> {kp,y}],{kd,200000,500000},{cp,100,1000}, {Pm,1000000,10000000}, {md,.3,20}] 80 Transformation of Cubic Damping Terms The following program was used to determine the transformed cubic terms when the modal transformation is applied. The terms are left in symbolic form and substituted in at the end. In addition, the assumption that mode 1 is zero is applied at the end. h=-r-I*i k=-r+I*i P={{1,h,k },{1,k,h},{1,1,1}} Cubic={{With[{x1=n1+k*n2+h*n3,x2=n1+h*n2+k*n3, x3=n1+n2+n3},x1^3]},{With[{x1=n1+k*n2+h*n3,x2=n1+h*n2+k*n3,x3=n1+n2+n3},x2^3]}, {With[{x1=n1+k*n2+h*n3,x2=n1+h*n2+k*n3,x3=n1+n2+n3},x3^3]}} Inverse[P].Cubic Collect[%,{n2,n3}] FullSimplify[%] Evaluate[%,n1=0,r=1/2, i=Sqrt[3]/2] 81 Three Cylinder Double Acting Engine The following program is used to obtain a numerical solution to the equations governing the motion of a three cylinder double acting free piston Stirling engine. The equations are derived based on the Schmidt model. The program begins with the specification of the nondimensional parameters. a=5 b=.4 g=.2 c=.3 mu=.1 The equations of motion are numerically solved. A second solution with varying initial conditions is obtained in order to verify that the system converges to a unique limit cycle. q=NDSolve[{x1''[t]+c*x1'[t]+mu*x1'[t]^3+x1[t]+a((1+b*x3[t]-g*x1[t])^-1-(1+b*x1[t]- g*x2[t])^-1)==0,x2''[t]+c*x2'[t]+mu*x2'[t]^3+x2[t]+a((1+b*x1[t]-g*x2[t])^-1-(1+b*x2[t]- g*x3[t])^-1)==0, x3''[t]+c*x3'[t]+mu*x3'[t]^3+x3[t]+a((1+b*x2[t]-g*x3[t])^-1-(1+b*x3[t]- g*x1[t])^-1)==0, x1[0]==1, x2[0]==0, x3[0]==0, x1'[0]==0, x2'[0]==0, x3'[0]==0},{x1,x2,x3},{t,0,300}, MaxSteps->10^6] qq=Table[NDSolve[{x1''[t]+c*x1'[t]+mu*x1'[t]^3+x1[t]+a((1+b*x3[t]-g*x1[t])^-1-(1+b*x1[t]- g*x2[t])^-1)==0,x2''[t]+c*x2'[t]+mu*x2'[t]^3+x2[t]+a((1+b*x1[t]-g*x2[t])^-1-(1+b*x2[t]- g*x3[t])^-1)==0, x3''[t]+c*x3'[t]+mu*x3'[t]^3+x3[t]+a((1+b*x2[t]-g*x3[t])^-1-(1+b*x3[t]- g*x1[t])^-1)==0, x1[0]==i, x2[0]==0, x3[0]==0, x1'[0]==0, x2'[0]==0, x3'[0]==0},{x1,x2,x3},{t,0,300}, MaxSteps->10^6],{i,0,1,.15}] The transient and stead state behavior are plotted. An array of discrete time values is also formed in order to determine the frequency and relative phase between the oscillating components. Plot[Evaluate[{x1[t],x2[t],x3[t]}/.q],{t,0,20},PlotRange->All, AxesLabel- >{\[Tau],Subscript[x,i]}] Plot[Evaluate[{x1[t],x2[t],x3[t]}/.q],{t,280,300},PlotRange->All, AxesLabel- >{\[Tau],Subscript[x,i]}] ParametricPlot3D[{Evaluate[{{x1[t], x2[t],x3[t]}/.q}]},{t,280,300},PlotRange->All,AxesLabel- > {Subscript[x,1],Subscript[x,2],Subscript[x,3]}] ParametricPlot3D[{Evaluate[{{x1[t], x2[t],x3[t]}/.qq}]},{t,280,300},PlotRange- >All,AxesLabel->{Subscript[x,1],Subscript[x,2],Subscript[x,3]}] Table[Evaluate[x1[i]/.q],{i,285,293,.05}] Table[Evaluate[x2[i]/.q],{i,285,293,.05}] Table[Evaluate[x3[i]/.q],{i,285,293,.05}] Max[Table[(Evaluate[x1[i]/.q]+Evaluate[x2[i]/.q]+Evaluate[x3[i]/.q])/3,{i,280,295,.05}]] The modal solutions are determined by using the coordinate transformation. 82 Plot[(Evaluate[x1[t]/.q]+Evaluate[x2[t]/.q]+Evaluate[x3[t]/.q])/3,{t,280,300}, AxesLabel->{t, Subscript[n,1]}] Realn2=Plot[Re[((.5(-1+I*Sqrt[3]))Evaluate[x1[t]/.q]+(.5(-1- I*Sqrt[3]))Evaluate[x2[t]/.q]+Evaluate[x3[t]/.q])/3],{t,280,300}, AxesLabel- >{t,Subscript[Re_n,2]}] Realn3=Plot[Re[((.5(-1-I*Sqrt[3]))Evaluate[x1[t]/.q]+(.5(- 1+I*Sqrt[3]))Evaluate[x2[t]/.q]+Evaluate[x3[t]/.q])/3],{t,280,300}, AxesLabel- >{t,Subscript[Re_n,3]}] Imaginaryn2=Plot[Im[((.5(-1+I*Sqrt[3]))Evaluate[x1[t]/.q]+(.5(-1- I*Sqrt[3]))Evaluate[x2[t]/.q]+Evaluate[x3[t]/.q])/3],{t,280, 300}, AxesLabel- >{t,Subscript[Im_n,2]}] Imaginaryn3=Plot[Im[((.5(-1-I*Sqrt[3]))Evaluate[x1[t]/.q]+(.5(- 1+I*Sqrt[3]))Evaluate[x2[t]/.q]+Evaluate[x3[t]/.q])/3],{t,280,300}, AxesLabel- >{t,Subscript[Im_n,3]}] 83 Nodal Methods Isothermal Single Cylinder The arrangement of a spring mounted piston inside of a cylinder is modeled by using a finite control volume scheme. The fluid within the cylinder is modeled with the ideal gas law; the fluid flow is assumed to be quasi-steady; upwinding is used to specify flux quantities. The cylinder is assumed to remain at some constant elevated temperature, whereas the ambient conditions are assumed to be constant. Bounded solutions are not obtained for this formulation, which is likely caused by physically incorrect simplifying assumptions. The parameters are defined first. All parameters are specified in SI units; the working fluid is air. The properties of the spring mass system are defined. mp=5 k=1000 c=.5 The ambient fluid conditions are defined. Pb=100000 The geometry of the cylinder and hole are defined as well as the loss coefficient. A=.3 lo=.2 Ah=.01 KL=10000 Fluid properties are defined. R=285 T=373 rho[t]=P[t]/(R*T) rhoa=1.205 An influx and outflux function are defined; then the density of the flux is defined. i[t]=-(Sign[P[t]-Pb]-1)/2 o[t]=(Sign[P[t]-Pb]+1)/2 rhof[t]=rho[t]*o[t]+rhoa*i[t] Various algebraic expressions are defined in order to simplify the input to the solver. V[t]=A*(lo-y[t]) Vdot[t]=-A*y'[t] m[t]=P[t]*V[t]/(R*T) mdot[t]=Ah*Sign[P[t]-Pb]*Sqrt[2*rhof[t]*Abs[P[t]-Pb]/KL] 84 The governing equations are numerically integrated. s=NDSolve[{mp*y''[t]==-k*y[t]-c*y'[t]+Ah*(Pb-P[t]),P'[t]==mdot[t]*R*T/V[t]- P[t]*(Vdot[t]/V[t]),y[0]==-.5, y'[0]==0, P[0] ==100000},{y ,P},{t,0,10}, MaxSteps->10^6] Various plots are obtained. Plot[Evaluate[{y[t]}/.s],{t,0,10}, AxesLabel->{t,y}] Plot[Evaluate[P[t]/.s],{t,0,10},PlotRange->All, AxesLabel->{t,P}] Plot[Evaluate[{m[t]}/.s],{t,0,10}, AxesLabel->{t,m}] 85 Non-Isothermal Single Cylinder The arrangement of a spring mounted piston inside of a cylinder is modeled by using a finite control volume scheme. The fluid within the cylinder is modeled with the ideal gas law; the fluid flow is assumed to be quasi-steady; upwinding is used to specify flux quantities. In addition, periodic heating is applied to the cylinder and a conduction path from the interior of the cylinder to the ambient air is allowed. The parameters are defined first. All parameters are specified in fundamental SI units. The working fluid is air. The properties of the spring mass system are defined. mp=6 k=100 c=.5 The geometry of the cylinder and hole, heat transfer parameters, and the loss coefficient are defined. lo = .4 A=.015 Ah=.05A Qdot = 100 h=20 KL= 1000 Fluid properties are defined as well as the ambient conditions. Cp = 1005 Cv = 718 R= 287 To = 300 Po = 100000 rhoo=1.161 Influx and outflux functions are defined and used to specify flux quantities. A function is also defined such that it is unity during cylinder expansion and zero during compression. i[t]=(Sign[Po-P[t]]+1)/2 o[t]=(1-Sign[Po-P[t]])/2 rho[t]=P[t]*V[t]/(R*T[t]) rhoflux[t]= i[t]*rhoo+o[t]P[t]*rho[t] Tflux[t] = i[t]To+o[t]T[t] Ex[t]=(Sign[Vdot[t]]+1)/2 Various algebraic expressions are defined in order to simplify the input to the solver. 86 mf[t] = Sign[Po-P[t]]*Ah*Sqrt[2*rhoflux[t]*Abs[P[t]-Po]/KL] V[t]=A*(lo-x[t]) Vdot[t]=-A*v[t] m[t] = P[t]*V[t]/(R*T[t]) The governing equations are numerically integrated. s=NDSolve[{x'[t]==v[t], v'[t] ==(A*(Po-P[t])-k*x[t]-c*v[t])/mp,T'[t]==(mf[t]*(Cp*Tflux[t]- Cv*T[t])+Qdot*Ex[t]-h(T[t]-To)-P[t]Vdot[t])/(Cv*m[t]), P'[t]==(mf[t]+m[t]*((mf[t]*(Cp*Tflux[t]-Cv*T[t])+Qdot*Ex[t]-h(T[t]-To)- P[t]Vdot[t])/(Cv*m[t]))/T[t]-P[t]Vdot[t]/(R*T[t]))*R*T[t]/V[t], x[0]==0, v[0]==-.0009, T[0]==300, P[0]==Po-.25},{x, v, T, P},{t,0,10},Method->StiffnessSwitching, MaxSteps-> 10^7] Various plots are obtained. Plot[Evaluate[x[t]/.s], {t, 0,10}, PlotRange ->All, AxesLabel->{t,x}] Plot[Evaluate[P[t]/.s], {t, 0,10}, PlotRange ->All, AxesLabel->{t,P}] Plot[Evaluate[T[t]/.s], {t, 0,10}, PlotRange ->All, AxesLabel->{t,T}] Plot[Evaluate[m[t]/.s], {t, 0,10}, PlotRange ->All, AxesLabel->{t,m}] 87 Isothermal Simulation of Sunpower RE-1000 Engine The Sunpower RE-1000 is a beta configuration free piston Stirling engine; the parameters are available in the literature, as well as experimentally measured performance data. The simplified nodal analysis is applied to this engine. However, the springs in this model are treated as linear mechanical springs rather than gas springs. All dimensional quantities are specified in fundamental SI units. The working fluid is helium. The spring stiffnesses, masses, loads, areas, and clearances are defined. Kp=190107 Kd=14760 Mp=6.2 Md=.426 cl=461.5 cd=35.34 Ap=25.69*10^-4 Ar=2.176*10^-4 Ad=Ap lp=18.3*10^-3 ld=18.61*10^-3 The engine charge pressure, absolute heater and cooler temperatures, and fluid properties are defined. Pm=71*10^5 Th=814.3 Tk=322.8 gamma=1.4 R=2077 The proportionality constant between the pressure drop and mass flow rate is defined. Initial mass contained within the compression and expansion spaces are defined. k=810000 mco=Pm*Ap*lp/(R*Tk) meo=Pm*Ad*ld/(R*Th) Several algebraic expressions are defined in order to simplify the input to the solver. Vc[t]=Ap*(lp-xp[t])+(Ad-Ar)*xd[t] Ve[t]=Ad*(ld-xd[t]) Vcdot[t]=-Ap*vp[t]+(Ad-Ar)*vd[t] Vedot[t]=-Ad*vd[t] Te[t]=Pe[t]*Ve[t]/(me[t]*R) 88 Tc[t]=Pc[t]*Vc[t]/(mc[t]*R) The equations are numerically integrated. s=NDSolve[{xp'[t]==vp[t], xd'[t]==vd[t],Mp*vp'[t]==-Kp*xp[t]-cl*vp[t]+Ap*(Pm- Pc[t]),Md*vd'[t]==-Kd*xd[t]-cd*vp[t]+Ad*(Pc[t]-Pe[t])+Ar*(Pm-Pc[t]),Pe'[t]==(- Pe[t]*Vedot[t]-k*R*Th((Pe[t]-Pc[t])))/Ve[t], Pc'[t]==(-Pc[t]*Vcdot[t]-k*R*Tk((Pc[t]- Pe[t])))/Vc[t], me'[t]==-(Pe[t]-Pc[t])/k, mc'[t]==-(Pc[t]-Pe[t])/k, xp[0]==0,xd[0]==0, vp[0]==.5,vd[0]==.001, Pe[0]==Pm, Pc[0]==Pm, me[0]==meo, mc[0]==mco},{xp,xd, vp, vd, Pe, Pc, me, mc},{t,1}, Method ->StiffnessSwitching, MaxSteps->10^7] Various plots of the output data are obtained. Plot[Evaluate[{xp[t]}/.s],{t,0,.387},AxesLabel->{t,xp}] Plot[Evaluate[{xd[t]}/.s],{t,0,.387},AxesLabel->{t,xd}] Plot[Evaluate[{Pe[t]}/.s],{t,0,.387},AxesLabel->{t,Pe}] Plot[Evaluate[{Pe[t]-Pc[t]}/.s],{t,0,.387},AxesLabel->{t,Pe-Pc}] Plot[Evaluate[{me[t]}/.s],{t,0,.387},AxesLabel->{t,me}] 89 Non-Isothermal Simulation of Sunpower RE-1000 Engine The Sunpower RE-1000 is a beta configuration free piston Stirling engine; the literature contains the engine parameters, as well as experimentally measured performance data. The simplified nodal analysis is applied to this engine, removing the isothermal assumption. However, the springs in this model are treated as linear mechanical springs rather than gas springs. All dimensional quantities are specified in fundamental SI units. The working fluid is helium. The spring stiffnesses, masses, loads, areas, and clearances are defined. Kp=190107 Kd=14760 Mp=6.2 Md=.426 cl=461 cd=35 Ap=25.69*10^-4 Ar=2.176*10^-4 Ad=Ap lp=18.3*10^-3 ld=18.61*10^-3 The engine charge pressure, absolute heater and cooler temperatures, and fluid properties are defined. Pm=71*10^5 Th=814.3 Tk=322.8 gamma=1.4 R=2077 The proportionality constant between the pressure drop and mass flow rate is defined. Initial mass contained within the compression and expansion spaces are defined. k=81000 mc0=Pm*Ap*lp/(R*Tk) me0=Pm*Ad*ld/(R*Th) Various algebraic expressions are defined to simplify the input to the solver. Vc[t]=Ap*(lp-xp[t])+(Ad-Ar)*xd[t] Ve[t]=Ad*(ld-xd[t]) Vcdot[t]=-Ap*vp[t]+(Ad-Ar)*vd[t] Vedot[t]=-Ad*vd[t] Te[t]=Pe[t]*Ve[t]/(me[t]*R) 90 Tc[t]=Pc[t]*Vc[t]/(mc[t]*R) Tfluxe[t]=.5*(Th*(Sign[Pc[t]-Pe[t]]+1)-Te[t](Sign[Pc[t]-Pe[t]]-1)) Tfluxc[t]=.5*(Tk*(Sign[Pe[t]-Pc[t]]+1)-Tc[t](Sign[Pe[t]-Pc[t]]-1)) The equations are numerically integrated. s=NDSolve[{xp'[t]==vp[t], xd'[t]==vd[t],Mp*vp'[t]==-Kp*xp[t]-cl*vp[t]+Ap*(Pm- Pc[t]),Md*vd'[t]==-Kd*xd[t]-cd*vd[t]+Ad*(Pc[t]-Pe[t])+Ar*(Pm-Pc[t]),Pe'[t]==gamma*(- Pe[t]*Vedot[t]-k*R*Tfluxe[t]((Pe[t]-Pc[t])))/Ve[t], Pc'[t]==gamma*(-Pc[t]*Vcdot[t]- k*R*Tfluxc[t]((Pc[t]-Pe[t])))/Vc[t],me'[t]==-(Pe[t]-Pc[t])/k, mc'[t]==-(Pc[t]- Pe[t])/k,xp[0]==0,xd[0]==0, vp[0]==.005,vd[0]==.001, Pe[0]==Pm, Pc[0]==Pm, me[0]==me0, mc[0]==mc0},{xp,xd, vp, vd, Pe, Pc, me, mc},{t,1}, Method->{StiffnessSwitching, Method- >{ExplicitRungeKutta,Automatic}}, MaxSteps->10^7] Various plots of the solution are obtained. Plot[Evaluate[{xp[t]}/.s],{t,0,1}, AxesLabel->{t,xp}, PlotRange->All] Plot[Evaluate[{xd[t]}/.s],{t,0,1}, AxesLabel->{t,xd},PlotRange->All] Plot[Evaluate[{Pc[t]}/.s],{t,0,1}, AxesLabel->{t,Pc},PlotRange->All] Plot[Evaluate[{mc[t]+me[t]}/.s],{t,0,1}, AxesLabel->{t,mc+me},PlotRange->All] 91 Bibliography [1] Walker, G., 1980, Stirling Engines, Clarendon Press, Oxford University Press, NY. [2] Senft, J. R., and Walker, G., 1985, Free Piston Stirling Engines, Springer-Verlag, NY. [3] Kankam, M. D., and Rauch, J. S., 1991, ?Comparative Survey of Dynamic Analyses of Free-Piston Stirling Engines,? NASA Technical Memorandum 104491. [4] Schreiber, J. G., 2001, ?Assessment of the Free Piston Stirling Converter as a Long Life Power Converter for Space,? NASA Technical Memorandum 210604. [5] Urieli, I., and Berchowitz, D. M., 1984, Stirling Cycle Engine Analysis, Adam Hilger LTD., Bristol, NY. [6] Ulusoy, N., 1994, ?Dynamic Analysis of Free Piston Stirling Engines,? Ph.D. thesis, http://www.ohiolink.edu/etd/view.cgi?case1061217408. [7] Chen, N. J. C., and Griffin, F. P., 1983, ?A Review of Stirling Engine Mathematical Models,? Oak Ridge National Laboratory Report, 135. [8] Ibrahim, M. B., and Tew, R. C., 2001, ?Study of Two-Dimensional Compressible Non- Acoustic Modeling of Stirling Machine Type Components,? NASA Technical Memorandum 211066. [9] Chen, N. J. C., Griffin, F. P., and West, C. D., 1984, ?Linear Harmonic Analysis of Stirling Engine Thermodynamics,? Oak Ridge National Laboratory Report, 155. [10] Chen, N. J. C., and Griffin, F. P., 1986, ?Linear Harmonic Analysis of Free Piston Stirling Engines,? Oak Ridge National Laboratory, 172. [11] 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. [12] Schreiber, J. G., Geng, S. M., and Lorenz, G. V., 1986, ?RE-1000 Free Piston Stirling Engine Sensitivity Test Results,? NASA Technical Memorandum 88846. [13] Schreiber, J. G., and Geng, S. M., 1987, ?RE-1000 Free Piston Stirling Engine Hydraulic Output System Description,? NASA Technical Memorandum 100185. [14] 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, pp 1036-1046. 92 [15] Schock, A., 1978, ?Stirling Engine Nodal Analysis Program,? Journal of Energy, 2, pp 354-362. [16] Organ, A. J., 1992, Thermodynamics & Gas Dynamics of the Stirling Cycle Machine, Cambridge University Press, NY, NY. [17] Nayfeh, A. H., and Balachandran, B., 1995, ?Applied Nonlinear Dynamics,? John Wiley & Sons, NY, NY. [18] Nayfeh, A. H., 1981, ?Introduction to Perturbation Techniques,? John Wiley & Sons, NY, NY. [19] Berchowitz, D. M., and Redlich, R. W., 1985, ?Linear Dynamics of Free-Piston Stirling Engines,? Proceedings of the Institution of Mechanical Engineers, 199, pp. 203-213. [20] Larsson, S., and Thomee, V., 2005, ?Partial Differential Equations with Numerical Methods,? Springer, Berlin, Germany. [21] Ben-Artzi, M., and Falcovitz, J., 2003, ?Generalized Riemann Problems in Computational Fluid Dynamics,? Cambridge University Press, NY, NY [22] Benvenuto, G., de Monte, F., and Farina, F., 1990, ?Dynamic Behaviour Prediction of Free-Piston Stirling Engines,? Proc. 25th Energy Conversion Engineering Conference, 5, pp. 346-351. [23] Benvenuto, G., and de Monte, F., 1996, ?The Effect of Nonlinear Thermo-Fluid- Dynamic Terms on Free-Piston Stirling Machine Stability,? Proc. 31st Energy Conversion Engineering Conference, Washington, DC, 2, pp. 1224-1231. [24] 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, pp 1047-1055. [25] Martini, W. R., 1983, ?Stirling Engine Design Manual,? NASA Technical Memorandum 168088. [26] Kankam, M. D., and Rauch, J. S., 1993, ?Controllability of Free-Piston Stirling Engine/Linear Alternator Driving a Dynamic Load,? NASA Technical Memorandum 106497. [27] Mason, S. L., Schreiber, J. G. and Thieme, L. G, 2001, ?Stirling Technology Development at GRC,? NASA Technical Memorandum 211315. [28] Banduric, R. D., and Chen, N. J. C., 1984, ?Nonlinear Analysis of Stirling Engine Thermodynamics,? Oak Ridge National Laboratory Report, 154. 93 [29] Rogdakis, E.D., Bormpilas, N.A., and Koniakos, I.K., 2003, ?A Thermodynamic Study for the Optimization of Stable Operation of Free Piston Stirling Engines,? Energy Conversion & Management, 45, pp. 575-593. [30] Schock, A., 1978, ?Nodal Analysis of Stirling Cycle Devices,? Proceedings of Intersociety Energy Conversion Engineering Conference, 3, pp 1771-1779. [31] Cheng, H., 2007, ?Advanced Analytic Methods in Applied Mathematics, Science, and Engineering,? LuBan Press, Boston, MA.