Linear and Nonlinear Modeling of a Traveling-Wave Thermoacoustic Heat Engine
UUnder consideration for publication in J. Fluid Mech. Linear and Nonlinear Modeling of aTraveling-Wave Thermoacoustic Heat Engine
C A R L O S C A L O † , S A N J I V A K. L E L E L A M B E R T U S H E S S E L I N K Center for Turbulence Research, Stanford, CA, 94305, USA Dept. of Aeronautics and Astronautics and Mechanical Eng., Stanford, CA, 94305, USA Dept. of Aeronautics and Astronautics and Electrical Eng., Stanford, CA, 94305, USA(Received ?; revised ?; accepted ?. - To be entered by editorial office)
We have carried out three-dimensional Navier-Stokes simulations, from quiescent con-ditions to the limit cycle, of a theoretical traveling-wave thermoacoustic heat engine(TAE) composed of a long variable-area resonator shrouding a smaller annular tube,which encloses the hot (HHX) and ambient (AHX) heat-exchangers, and the regenerator(REG). Simulations are wall-resolved, with no-slip and adiabatic conditions enforced atall boundaries, while the heat transfer and drag due to the REG and HXs are modeled.HHX temperatures have been investigated in the range 440K – 500K with AHX temper-ature fixed at 300K. The initial exponential growth of acoustic energy is due to a networkof traveling waves thermoacoustically amplified by looping around the REG/HX unit inthe direction of the imposed temperature gradient. A simple analytical model demon-strates that such instability is a localized Lagrangian thermodynamic process resemblinga Stirling cycle. An inviscid system-wide linear stability model based on Rott’s theoryis able to accurately predict the operating frequency and the growth rate, exhibitingproperties consistent with a supercritical Hopf bifurcation. The limit cycle is governedby acoustic streaming – a rectified steady flow resulting from high-amplitude nonlinearacoustics. Its key features are explained with an axially symmetric incompressible modeldriven by the wave-induced stresses extracted from the compressible calculations. Thesefeatures include Gedeon streaming, Rayleigh streaming in the resonator, and mean re-circulations due to flow separation. The first drives the mean advection of hot fluid fromthe HHX to a secondary heat exchanger (AHX2), in the thermal buffer tube (TBT),necessary to achieve saturation of the acoustic energy growth. The direct evaluation ofthe nonlinear energy fluxes reveals that the efficiency of the device deteriorates with thedrive ratio and that the acoustic power in the TBT is balanced primarily by the meanadvection and thermoacoustic heat transport.
Key words:
Authors should not enter keywords on the manuscript
1. Introduction
Thermoacoustic heat engines (TAE) are devices that can convert available thermal en-ergy into acoustic power with very high efficiencies. This potential is due to the absenceof moving parts and relative simplicity of the components. This results in low manufac-turing and maintenance costs making these systems an attractive alternative for clean † Email address for correspondence: [email protected] a r X i v : . [ phy s i c s . f l u - dyn ] A ug Carlo Scalo, Sanjiva K. Lele, Lambertus Hesselink and effective energy generation or waste-energy reutilization. The core energy conversionprocess occurs in the regenerator – a porous metallic block, placed between a hot anda cold heat-exchanger, sustaining a mean temperature gradient in the axial direction.Acoustic waves propagating through it (under the right conditions) can be amplified viaa thermodynamic process resembling a Stirling cycle. Most designs explored up to themid 1980’s were based on standing waves and had efficiencies typically less than 5%. Asignificant breakthrough was made by Ceperley (1979) who showed that traveling-wavescan extract acoustic energy more efficiently, leading to the design concept for traveling-wave TAEs (Backhaus & Swift 2000; de Waele 2009). In this configuration the generatedacoustic power is in part resupplied to the regenerator via some form of feedback andin part directed towards a resonator for energy extraction. A secondary ambient heatexchanger is typically needed to contain the heat leaking (due to nonlinear effects) fromthe hot-heat exchanger. This design is the focus of the present study.Improving the technology behind TAEs has been of particular interest in the lastdecade with research efforts being made worldwide (see Garret (2004) for a review). Arecent breakthrough, for example, has been made by Tijani & Spoelstra (2011) who de-signed a traveling-wave TAE achieving a remarkable overall efficiency of 49% of theCarnot limit. The state-of-the-art prediction capabilities and technological design ofTAEs can, however, significantly benefit from a high-fidelity description of the under-lying fluid mechanics. The potential of such approach to fill important modeling gapshas inspired the present study which relies on full-scale three-dimensional flow simula-tions to gain insight into the linear and nonlinear processes occuring in a theoreticaltraveling-wave thermoacoustic engine.A comprehensive theoretical analysis of thermoacoustic effects in ducts is provided inthe seminal publications by Rott and co-workers (Rott 1969, 1973, 1974, 1975 a , b , 1976;Zouzoulas & Rott 1976; Rott 1980; M¨uller & Rott 1983; Rott 1984) where a predictiveanalytical framework (restricted to simple configurations) is derived improving upon pre-existing theories by Kirchhoff (1868) and Kramers (1949). Issues addressed include theonset of instability, thermoacoustic heating, transport due to acoustic nonlinearities andeffects of variable cross-sectional area. Later, Swift and co-workers used Rott’s work as thebasis for the development of semi-empirical low-order models for the acoustics in variouscomponents found in real thermoacoustic engines (Swift 1988, 1992; Swift & Ward 1996).This resulted in the development of the prediction software package DeltaE (Ward &Swift 1994) (replaced now by
DeltaEC ), which, together with similar modeling toolssuch as
SAGE and
REGEN , is still actively used in the academic literature as well as inindustry. Other examples of advanced low-order modeling relying on Rott’s theory andsystematic asymptotic approximations (Bauwens 1996; In ’T Panhuis et al. et al. odeling of a Thermoacoustic Heat Engine et al. et al. a , b , 2012)identified in the inadequate modeling of such nonlinear effects the primary reason forthe failure of low-order models to correctly capture wave-amplitude saturation, even insimple geometries. It is therefore necessary to adopt a direct modeling approach as doneby Boluriaan & Morris (2003) who performed simulations solving the fully compressibleviscous flow equations in an idealized two-dimensional configuration modeling traveling-wave streaming suppressed by a jet pump. To properly account for the viscous interactionswith the solid wall, the Stokes thickness needs to be resolved. Three-dimensional flowsimulations of similar configurations, fully resolving thermo-viscous effects and transportdue to streaming, are yet to be attempted. In the present work we take on the challengeof studying streaming occurring in a three-dimensional flow with geometric complex-ities by analyzing its effects on the engine performance but also directly modeling itwith a vorticity-streamfunction formulation. A fluid dynamic analogy can be exploitedto model the streaming (Rudenko & Soluyan 1977) as an incompressible flow driven bywave-induced Reynolds stresses. By following this approach, we developed a simplifiednumerical model to gain insight into the spatial structure of the acoustic stresses and theirrelationship with complex geometrical features. The model reproduces the key nonlineareffects such as Rayleigh and Gedeon streaming (Gedeon 1997).A very high computational cost, however, is associated with the direct resolution ofthe governing equations. In a full TAE the range of temporal and spatial scales can span4 orders of magnitude (Hamilton et al. δ ν typically of the order of 10 − mm, to the resonator length, typically of the orderof the acoustic wavelength λ (cid:39) et al. (2005), extending it to a three-dimensional setup and adding a secondary ambient heat-exchanger necessary to achieve a limit cycle. The heat-exchangers and regenerators aremodeled using semi-empirical source terms available in the literature. This allows usto focus on the full-scale resolution of the high-amplitude acoustics interacting withcomplex geometrical features. Details omitted in `a Nijeholt et al. (2005) regarding thegeometry and the modeling of the heat-exchangers and regenerators are reconstructedto the best of the authors’ ability with the help of Ray Hixon (pers. comm., 2012).The present work can be regarded as the first step towards a simple benchmark casefor validation of computational modeling of thermoacoustic devices. The spatially andtemporally resolved data will be used to gain insight into the aforementioned linearand nonlinear governing physical processes and contribute to their modeling. The lack of Carlo Scalo, Sanjiva K. Lele, Lambertus Hesselink
Figure 1.
Illustration of the model traveling-wave thermoacoustic engine inspired by `a Nijeholt et al. (2005). Full system (top), regenerator/heat-exchanger area (bottom), drawn to scale. Adash-dotted line indicates the line of symmetry. The different components of the engine arethe resonator (R), the annular feedback inertance loop (i), the compliance (c), the hot heatexchanger (HHX), the regenerator (REG), the ambient heat exchanger (AHX), secondary heatexchanger (AHX2), and the thermal buffer tube (TBT). The annular tube enclosing the AHX2,HHX, REG, and AHX is also refered to as pulse tube, of which the TBT is only a section. experimental data available for the proposed theoretical device has required the adoptionof several companion lower-order models to verify the results obtained from the three-dimensional simulations. The levels of modeling adopted range from zero-dimensionaland purely analytical to axially symmetric and nonlinear.In the following, the computational setup is first introduced, discussing the adoptedmeshing strategy and the semi-empirical heat-transfer and drag models for the heat-exchangers and the regenerator. Results follow investigating first the start-up phase,where linear models are adopted to describe the nature of the wave propagation through-out the system and amplification via thermoacoustic instability. Instantaneous data isthen collected at the limit-cycle where acoustic streaming is investigated. Results from anincompressible numerical model used to directly solver for streaming flow are discussed.Finally, thermal energy budgets in the TBT are analyzed.
2. Problem formulation
The full-scale approach advocated in the present study requires the resolution of thecomplete set of compressible viscous flow equations. The conservation equations for mass,momentum and total energy are reported here in index notation ∂ρ∂t + ∂ρu i ∂x i = 0 (2.1 a ) ∂ρu i ∂t + ∂ρu i u j ∂x j + ∂τ ij ∂x j = − ∂p∂x i + D i (2.1 b ) ∂ρE t ∂t + ∂∂x j [( ρE t + p ) u j ] + ∂β j ∂x j = S E (2.1 c ) odeling of a Thermoacoustic Heat Engine x , x and x (or x , y and z ) are the axial and cross-sectional coordinates, u i thevelocity components in those directions, p , ρ and E t are the pressure, density and totalenergy. The viscous stress tensor and energy flux, respectively, τ ij and β i , are given by β i = − u j τ ij + q i (2.2 a ) τ ij = − µ (cid:20) S ij − ∂u k ∂x k δ ij (cid:21) (2.2 b )where q i is the molecular heat-flux, µ the dynamic viscosity, S ij the deviatoric part ofthe strain-rate tensor. The fluid is assumed to be an ideal gas with reference state givenby ρ ref = 1 . , p ref = 101325 Pa and T ref = 300 K. The dynamic viscosity varieswith temperature based on the power law µ ( T ) = µ ref ( T /T ref ) . . The Prandtl numberis P r = 0 . et al. (2005) we choose to model the heat transfer and drag in suchcomponents via the source terms D i and S E on the right-hand side of (2.1b) and (2.1c).They are expressed as (cid:40) D i = − (cid:104) R C + R F ( u j u j ) / (cid:105) u i (2.3 a ) S E = − u i D i + S h (2.3 b )where drag term D i is modeled following the parametrizations R C = C sf µ (1 − φ )4 d w φ (2.4 a ) R F = ρC fd d w (2.4 b )where φ and d w are the characteristic porosity and mesh wire size of the component, and C sf and C fd are dimensionless fitting constants taken from the ILK Dresden and K¨uhlmetal felts correlations (Thomas & Pittman 2000). These are specific to TAE regeneratorsand are derived under different oscillating flow conditions. Unfortunately, `a Nijeholt et al. (2005) does not provide values for the wire size d w , which have been estimated by usingthe correlation suggested by Organ (1992) r h = d w φ − φ ) , (2.5)where r h is the pore hydraulic radius.The source term S h in (2.3b) accounts for heat-transfer in the REG/HX unit and ismodeled as S h = − α T (cid:2) T − T ( x ) (cid:3) , (2.6)where T is the instantaneous fluid temperature and T ( x ) is the target mean temperatureprofile which is equal to T h and T a , respectively, in the hot and ambient heat exchangers,and varies linearly in the regenerator between the two values. No information is provided Carlo Scalo, Sanjiva K. Lele, Lambertus Hesselink
Table 1.
Parameters for regenerator/heat-exchanger model (2.3). Hydraulic radius r h , charac-teristic wire diameter d w (mm), porosity φ , and drag coefficients C sf , C fd (Thomas & Pittman2000). Values for r h and φ are taken from `a Nijeholt et al. (2005).Parameter Heat Exchangers Regenerator r h (mm) 0.1 0.041 d w (mm) 0.2667 0.0670 φ C sf C fd in `a Nijeholt et al. (2005) for the proportionality constant α . We propose to model it as(Ray Hixon, pers. comm, 2012) α T = α h ρRγ − R is the gas constant and γ is the ratio of specific heats. The ratio ρR/ ( γ − ∂ρE t /∂T ) derived using the equation of state. The constant α h is defined as α h = 1 /τ h (2.8)where τ h is the characteristic time scale for heat-transfer in the void spaces of the heatexchangers and regenerator. An estimate for τ h can be derived by modeling such com-ponents as stacks of parallel plates with spacing b matching the given hydraulic radiusof the pores, b = 2 r h (table 1). This results in (Bejan 2004) τ h = ( b/ k/ρC p = b ρP r µ (2.9)where k is the thermal conductivity, C p the specific heat capacity of the gas.This simplified model is expected to predict the intensity of the heat transfer rate tothe pore flow within an order of magnitude. It is based on the assumption of perfectthermal contact and is in quantitative agreement with a similar model used in DeltaEC (Ward & Swift 1994). While the thermal regime resulting from (2.6) is not affected by theflow velocity, its linear dependency from the temperature facilitates the lower-order mod-eling efforts made in the present manuscript. Overall, despite their simplicity and coarseapproximation, the application of the source terms (2.3) reproduces the essential ther-modynamic and hydrodynamic processes that occur in regenerators and heat-exchangersin TAEs, as discussed in the following.
3. Numerical Model
The governing equations are discretized on an unstructured hexahedral mesh (figure 2)and solved in a 90 ◦ sector with rotational periodicity applied in the azimuthal direction.The infinitely thin annular tube wall is introduced by breaking the mesh connectivity,creating two overlapping boundary surfaces with opposite orientation. Three concentric O-grids in the cross-section, two in the annular tube and one in the pulse tube (notshown), are required in this region to map the polar mesh at the resonator walls to aquasi-uniform Cartesian block at the center. High resolution is retained near the sharp odeling of a Thermoacoustic Heat Engine (a)(b) Figure 2.
Cross-section of three-dimensional computational grid A for the full length resonator(a), zoom on the right end (b) corresponding to view in figure 1. The computational grid hasbeen mirrored about the centerline for illustrative purposes. Computations are performed in a90 ◦ sector. edges of the annular tube walls to properly capture the shear layer caused by periodicflow separation. Visual inspection of the flow in previous calculations does not reveal asignificant vorticity magnitude away from the wall for x < − . x =-0.146 m (figure 2b), resulting in a coarser radial grid distribution for x < -0.146m. Points have also been concentrated in the AHX2 ( x = 0.031 m), due to intenseinstantaneous temperature and velocity gradients at the limit cycle created by hot fluidstreaming away from the HHX (discussed later). Overall, a significant effort has beenmade to retain a high-quality structured grid when possible.A preliminary grid refinement study has been carried out to ensure adequate resolutionof the axially symmetric components of the acoustic field and accurate prediction of thegrowth rate. The latter is sensitive primarily to the resolution in the axial direction(main direction of propagation of the acoustic waves), both in the resonator and in thethermal buffer tube. The viscous boundary layers are resolved with resolution of 0.1mm at the wall for all cases. These considerations have lead to the design of a baselinegrid distribution, grid A, (figure 2) used to rapidly advance in time through the initialtransient. Simulations have been carried out for hot-heat exchanger temperatures of T h = 440K, 460K, 480K and 500K and ambient heat-exchanger temperature of T a = 300Kin all cases. The acoustic perturbation is initialized with a standing wave of 0.5 kPa ofamplitude and the source terms (2.3) are applied from the beginning. The former is onlyused to reduce the duration of the transient and is not required to achieve acoustic energygrowth (see section 4). Once a limit cycle is reached, two successive grid refinement stepsare carried out, resulting in grid B and C. At each step the resolution was increasedin the axial direction, especially around the sharp edges, and systematically doubled inthe azimuthal direction. The mesh size in the radial direction is then adjusted and/orincreased to optimize the cells aspect ratio. The sensitivity of the wave-induced Reynoldsstresses to these changes (shown later) is used as a metric for grid convergence and onlycarried for calculations at T h =500 K (table 2).The governing equations for mass, momentum and total energy are solved in thefinite-volume unstructured code CharLES X developed as a joint-effort project amongresearchers at Stanford University. The flux reconstruction method is grid-adaptive atthe preprocessing stage and solution-adaptive at run-time. It blends a high-order poly-nomial interpolation scheme (up to fourth-order on uniform meshes) with a lower-orderscheme to ensure numerical stability in areas of low grid quality (Ham et al. Carlo Scalo, Sanjiva K. Lele, Lambertus Hesselink
Table 2.
Parameter space for numerical simulations carried out to limit cycle. Total number ofcontrol volumes N cv, control volumes in the azimuthal direction N θ , temperature in the HHX T h . For all cases T a = 300K. Three meshes with increasing level of resolution and quality (fromGrid A to C) are considered. Computations performed (x) and not performed ( · ).Grid Type N cv N θ T h = 460K T h = 480K T h = 500KA 0.46m 20 (x) (x) (x)B 1.25m 40 ( · ) ( · ) (x)C 5.08m 80 (x) ( · ) (x) .
00 0 .
05 0 .
10 0 .
15 0 .
20 0 .
25 0 .
30 0 .
35 0 .
40 0 .
45 0 .
50 0 .
55 0 . t (s) − − U ( m / s ) .
00 0 .
05 0 .
10 0 .
15 0 .
20 0 .
25 0 .
30 0 .
35 0 .
40 0 .
45 0 .
50 0 .
55 0 . − − P − P a t m ( k P a ) (a)(b) Figure 3.
Time series of pressure (a) along the centerline, r = 0, at x = − . x = 0 . x = − .
92 m (—) for case T h = 500 K on grid A. unwanted oscillations in the solution caused by the application of the source terms (2.3).The discretized system of equations is integrated in time with a fully-explicit, third-orderRunge-Kutta scheme. The code is parallelized using the Message Passing Interface (MPI)protocol and highly scalable on a large number of processors. The adoption of computa-tionally intensive discretizations such as ENO in a limited portion of the domain has leadto a load-balancing problem that required a volume-based dual-constrained partitioning(Karypis & Kumar 1998; Schloegel et al.
4. Engine Start-Up
In all of the numerical trials performed, the abrupt activation of the source terms (2.3)alone in a quiescent flow provides a sufficiently intense initial disturbance ( ∼ ∼ odeling of a Thermoacoustic Heat Engine
150 155 160 165 170 175 180 185 x (mm) T ( K ) AHXREGHHX . . . . . . . t (s) x L ( mm ) A H X RE G HH X − . − .
04 0 .
00 0 .
04 0 . v L (m /kg) − − p L ( k P a ) ◦ ◦ ◦ ◦ . . . . . . . ρ ( k g / m ) − . − . . . . u L ( m / s ) ◦ ◦ ◦ ◦ (a) (b)(c) Figure 4.
Base temperature (—) and density profiles (- -) taken along the centerline in theREG/HX unit at t = 0 .
55s for T h = 500 K, grid A (see table 2) (a); Lagrangian fluid parcelaxial position x L starting at x = 0 .
171 m at t = 0 (b); fluctuating pressure, p (cid:48) L , (—) andLagrangian velocity, u (cid:48) L , (- -) plotted against specific volume fluctuation, v (cid:48) L , over one acousticcycle t = 0 .
55s (c) with heat release fluctuation, S h (cid:48) (2.6) shown with arrows (oriented upwardsfor heat addition S (cid:48) h > S (cid:48) h < x < Thermodynamic Cycle
The driver of the thermoacoustic instability, converting heat into acoustic power, is themean temperature gradient imposed in the REG/HX unit (figure 4a). Insight into thefundamental energy conversion mechanisms can be gained by looking at the evolutionof a Lagrangian fluid parcel in the regenerator interacting with the acoustic field. Theslight drift in the direction of the mean temperature gradient (figure 4b) is a nonlineareffect known as acoustic streaming (discussed later in section 5), which can be ignoredat this stage.The fluid parcel in the regenerator experiences a thermodynamic cycle converting heatinto acoustic power, which is neither the ideal Stirling or Carnot cycle (figure 4c). Forexample, purely isochoric transformations, present in the ideal Stirling cycle, are not0
Carlo Scalo, Sanjiva K. Lele, Lambertus Hesselink possible due to the sinusoidal waveform of the acoustic velocity and pressure. Moreover,isentropic transformations, present in the ideal Carnot cycle, are not possible due tothe heat-exchange and viscous losses in the REG/HX unit. However, analogously tothe ideal Stirling cycle, most of the heating and cooling occurs, respectively, during theexpansion and compression stages. The heat transfer model in (2.6) assumes perfectthermal contact and is therefore likely that its adoption leads to an overestimation of thereal thermoacoustic response under comparable conditions.Due to the orientation of the background temperature gradient ( d T ( x ) /dx <
0, fig-ure 4a), a fluid parcel in the regenerator that is displaced towards the hot heat exchanger( u (cid:48) <
0) will experience heating. This occurs with a given phase-lag with respect to thevelocity depending, in particular, on the nature of the heat-transfer. In our case the heataddition, S (cid:48) h >
0, peaks ∼ ◦ after the maximum negative velocity (when the parcelcomes to rest) and so does the positive pressure fluctuation, p (cid:48) > S (cid:48) h p (cid:48) >
0, the ob-served phase differences suggest the presence of positive acoustic energy production. Thephasing between velocity and pressure fluctuations is consistent with a standing wave,short of approximately 15 ◦ . This slight phase difference contributes to a negative corre-lation between u (cid:48) and p (cid:48) , i.e. the left-traveling wave propagating through the REG/HX ismore intense than the right-traveling wave. Acoustic power is therefore being produced.Inspired by Swift (1988)’s theoretical thermoacoustic engine (figure 5, top), a one-dimensional analytical model is derived in the following to rigorously explain the inter-action between a mean temperature gradient and a Lagrangian fluid parcel in a genericplanar acoustic wave field. The latter can expressed as the linear superposition of a left-and a right-traveling wave, which in complex form readsˆ p ( x ) = p − e i [ k x + φ − ] + p + e i [ − k x − φ + ] (4.1) ρ a ˆ u ( x ) = − p − e i [ k x + φ − ] + p + e i [ − k x − φ + ] , (4.2)where p + / − is the amplitude of the right/left- traveling wave and k is the wave number.The acoustic pressure and velocity in time are given based on the convention p (cid:48) a =ˆ p ( x ) e i ω t and u (cid:48) a = ˆ u ( x ) e i ω t , where i is the imaginary unit and the base impedance ρ a is given by the state { ρ , T , P } .Let x p = x + x (cid:48) p be the instantaneous position of a fluid parcel oscillating with smalldisplacements x (cid:48) p about the position x , where a linear temperature gradient is locally im-posed (figure 5). The fluid parcel velocity can be approximated, based on the assumption k max( x (cid:48) p ) <<
1, as ddt x (cid:48) p = u (cid:48) ( x + x (cid:48) p , t ) = R{ ˆ u ( x + x (cid:48) p ) e i ω t } (cid:39) R{ ˆ u ( x ) e i ω t } (4.3)yielding the relation ˆ x p = ˆ u ( x ) /iω. (4.4)Introducing a Lagrangian base state ρ L, , P L, , T L, (specified later) for the fluid parceldensity, entropy and temperature, ρ L = ρ L, + R (cid:8) ˆ ρ L e i ω t (cid:9) (4.5) s L = s L, + R (cid:8) ˆ s L e i ω t (cid:9) (4.6) T L = T L, + R (cid:110) ˆ T L e i ω t (cid:111) , (4.7) odeling of a Thermoacoustic Heat Engine x p , yields ρ L T L ds L dt = − α T [ T L − T ( x p )] (4.8)where the (total) time derivative on the l.h.s., applied to the Lagrangian fluid parcel’sspecific entropy, has replaced the material derivative of the Eulerian entropy field.Substituting Gibbs’ relationship linearized about the Lagrangian base state T L, , p L, (= p ), ds L = c p T L, dT L − Rp dp L (4.9)into (4.8), which is also linearized assuming ρ L T L (cid:39) ρ L, T L, , yields ρ L, c p dT L dt − d p L dt = − α T [ T L − T ( x p )] . (4.10)The imposed mean temperature at the particle location, T ( x p ), can be expressed viathe Taylor expansion T ( x p ) = T ( x ) + x (cid:48) p ∇ T , (4.11)which is exact in the case of linear temperature profile. Substituting (4.11) into (4.10)yields ρ L, c p dT L dt − d p L dt = − α T (cid:2) T L − T ( x ) − x (cid:48) p ∇ T (cid:3) (4.12)Defining the Lagrangian based state for temperature such that T L, = T ( x ), assuming p L = p and switching to complex form ρ L, c p iω ˆ T L − iω ˆ p = − α T (cid:104) ˆ T L − ˆ x p ∇ T (cid:105) (4.13)where the Lagrangian base density is ρ L, = p / ( R T L, ).If the acoustic field is assigned, so are the complex pressure amplitude ˆ p and particledisplacement ˆ x p = ( iω ) − ˆ u ( x ), which then allows to solve for ˆ T L from (4.13). The densityof the Lagrangian parcel can then be calculated from the linearized equation of stateˆ ρ L = ˆ p − ρ L, R ˆ T L R T L, . (4.14)The work done by the fluid parcel on the surrounding ambient per unit time (generatedacoustic power) is ˙ w = − pρ DρDt (cid:39) − ρ L, R{ ˆ p ( i ω ˆ ρ L ) ∗ } . (4.15)Depending on the nature of acoustic wave being imposed (figure 5a), ranging from purelyleft-traveling to purely right-traveling, a different phasing between ˆ p and ˆ u is achieved,determining ˙ w (figure 5b). In the case of a standing wave, acoustic energy will be absorbed( ˙ w <
0) or generated ( ˙ w >
0) depending on the location of the temperature gradient, x . For a sufficiently high amplitude of the right-, p + >> p − , (or left-, p − >> p + )traveling wave, acoustic power is ultimately only absorbed (or generated) for any x (figure 5b). This shows that an acoustic wave traveling in the same or opposite directionof an imposed mean temperature gradient (applied over a region small compared to thewavelength) will be amplified or absorbed. Moreover, the acoustic power associated withthe energy conversion occurring in an (almost) purely traveling wave is remarkably higher(see area enclosed by the p − v diagrams in figure 5c) than the one of a standing wave of2 Carlo Scalo, Sanjiva K. Lele, Lambertus Hesselink π/ π/ π/ πkx | ˆ p | / P × π/ π/ π/ πkx − − ˙ w / P ω × − . − .
15 0 .
00 0 .
15 0 . v L /v L, − − − − p L / P × (a)(b) (c) Figure 5.
Results from the linear Lagrangian model (section 4.1) of a one-dimensional the-oretical thermoacoustic engine composed of a generic plane wave interacting with a negativetemperature gradient located at x (top). Imposed pressure distribution as a function of k x (a),overall cycle-averaged work (i.e. generated acoustic power) as a function of kx (temperature gra-dient location) (b), and Lagrangian thermodynamic cycle extracted in the temperature gradientregion for k x = π/ p − = 3 × − P , p + = 1 × − P (––); p + = p − = 2 . × − P (– · –); p − = 1 × − P , p + = 3 × − P (– –) for P = ρ a .For all cases φ − = φ + = − π/ comparable amplitude. This confirms Ceperley (1979)’s seminal intuition that led to therevolutionary concept of traveling-wave energy conversion, trumping thereafter designsbased on standing waves.4.2. Acoustic Network of Traveling Waves
In spite of the finite amplitude of the initial perturbation (exceeding 1% of the base pres-sure) and the immediate establishment of nonlinear effects, the exponentially growingacoustic amplitude (with uniform growth rate in the entire system) suggests that thesystem-wide behavior in the start-up phase can be analyzed by invoking linear acous-tics. The low frequencies observed in the numerical simulations and the high-aspect ratioof the resonator (with lowest cut-on frequency ∼ odeling of a Thermoacoustic Heat Engine − . − . − . − . − . − . − . − .
25 0 .
00 0 . . . . . . . . p − , p + ( k P a ) − . − . − . − . − . − . − . − .
25 0 .
00 0 . x (m) − π - π π π φ − , φ + ( r a d ) (a)(b) Figure 6.
Left ( − , (cid:47) , (cid:74) ) and right (+, (cid:46) , (cid:73) ) traveling-wave amplitudes, p ± , (a) and phases, φ ± , (b) calculated based on the linear approximation (4.16) for time t = 0 .
55 s for T h = 500 Kcase on grid A. The locations where data is extracted are shown in the top figure. White andblack symbols correspond to values extracted on the centerline and in the feedback inertance,respectively. (figure 5c) and propagating freely in the latter. The acoustic power propagating throughthe inertance is looped back into the REG/HX unit via the compliance, hence creatinga network of self-amplified traveling-waves.This picture can be confirmed with the aid of the instantaneous numerical data. Anexact local decomposition in terms of left ( − ) and right (+) traveling waves (cid:26) p (cid:48) ( t ) = p − f ( ω t + φ − ) + p + f ( − ω t + φ + ) , (4.16 a ) u (cid:48) ( t ) = − u − f ( ω t + φ − ) + u + f ( − ω t + φ + ) , (4.16 b )can allow for the direct evaluation of the amplitudes p ± and u ± and phases φ ± of purelytraveling waves of a given generic waveform f () (periodic function of period 2 π ). In ourcase, the angular frequency ω is much larger than the growth rate (discussed in section4.3), allowing to ignore the variations of pressure and velocity amplitudes over one acous-tic period as well as the change in base impedance, ρ a , due to a DC mode in pressure(discussed above). The following analysis is restricted to the start-up phase and for lo-cations in the engine outside of the REG/HX unit where the strong mean temperaturecauses non-negligible spatial gradients of the base impedance and the isentropic wavepropagation assumption to be violated.The acoustic perturbation (4.16) can be rewritten in the form of a complex Fourier4 Carlo Scalo, Sanjiva K. Lele, Lambertus Hesselink t (s) (b)(c) (a) P a m p ( k P a ) . . . . P a m p ( d B ) Figure 7.
Time series of pressure amplitude (kPa) on the left end of the resonator and corre-sponding sound pressure level (+dB) for T h = 460K (a), T h = 480K (b) and T h = 500K (c) ongrid A (see table 2). Decaying case for T h = 440K not shown. series + ∞ (cid:88) k = −∞ ˆ p k e i k ω t = p − + ∞ (cid:88) k = −∞ ˆ f k e i k [ ω t + φ − ] + p + + ∞ (cid:88) k = −∞ ˆ f ∗ k e i k [ ω t − φ + ] , (4.17 a ) ρ o a o + ∞ (cid:88) k = −∞ ˆ u k e i k ω t = − p − + ∞ (cid:88) k = −∞ ˆ f k e i k [ ω t + φ − ] + p + + ∞ (cid:88) k = −∞ ˆ f ∗ k e i k [ ω t − φ + ] . (4.17 b )where p ± = ± ρ a u ± and the superscript ∗ indicates the complex conjugate. Lettingthe m -th mode be any non-zero Fourier component, the unknowns p ± and φ ± are easilydetermined by isolating such mode from (4.17), yielding ˆ p m = ˆ f m (cid:104) p − e i m φ − (cid:105) + ˆ f ∗ m (cid:104) p + e − i m φ + (cid:105) , (4.18 a ) ρ o a o ˆ u m = − ˆ f m (cid:104) p − e i m φ − (cid:105) + ˆ f ∗ m (cid:104) p + e − i m φ + (cid:105) . (4.18 b )The amplitudes and phases of the modes in (4.16) can be obtained from the magnitudeand phases of the complex unknowns grouped into squared brackets above, yielding (cid:40) p − = | ˆ p m − ρ a ˆ u m | / | f m | (4.19 a ) p + = | ρ a ˆ u m + ˆ p m | / | f ∗ m | . (4.19 b )The phases φ − and φ + are then readily extracted from (4.18) given the pressure ampli-tudes p − and p + .This procedure is applied to the discrete set of points in figure 6(top) located alongthe resonator axis and around the pulse tube. Results show that left-traveling wavesleaving the REG/HX unit propagate into the resonator and are reflected back withslightly lower amplitude due to losses in the resonator (figure 6a). Consistently with theone-dimensional approximation in (4.16) the acoustic power can be expressed as˙ E + a = (cid:90) A p (cid:48) u (cid:48) dA = Aρ a (cid:104) p +2 − p − (cid:105) + ∞ (cid:88) k = −∞ ˆ f ∗ k ˆ f k . (4.20)providing an energetic interpretation to the imbalance p − ≷ p + . The data (figure 6a,b)confirms that, as anticipated earlier in this section, the acoustic power generated in the odeling of a Thermoacoustic Heat Engine p +2 >>p − (directly responsible for nonlinearities such as Gedeon streaming, analyzed later insection 5). The spatial distribution of φ − and φ + (figure 6b) is consistent with a standingwave at the resonator scale with a slight deviation in the REG/HX and in the inertance.Such difference is sufficiently small to suggest that it is in fact the same planar wavepropagating through the REG/HX and the inertance in each direction at once but withpart of its acoustic power amplified or damped in the former depending on the direction ofpropagation. Overall, the present results confirm the qualitative picture outlined earlier inthis section, explaining the instability mechanisms leading to the acoustic energy growth(figure 7). 4.3. System-wide Linear Modeling
Results shown so far suggest that nonlinearities do not play an important role in explain-ing the acoustic energy propagation and amplification mechanisms during the start-upphase. However, given the high amplitude of the initial perturbation ( ∼ ddt p (cid:48) c = w c (cid:104) U (cid:48) i − U (cid:48) t (cid:105) (4.21)where p (cid:48) c , U (cid:48) i and U (cid:48) t are, respectively, the instantaneous fluctuating pressure in thecompliance and flow rates exchanged with the inertance through surface i , and with thepulse tube through t (figure 8), and w c = γP /V c where P is the base pressure and V c the volume of the compliance.The very small variation among the growth rates and frequencies extracted from thenumerical simulations at the locations in figure 6 (not shown) suggests that normalmodes can be assumed for all fluctuating quantities. By adopting the e + iσ t conventionwith σ = − iα + ω , where α and ω are, respectively, the growth rate and the angularfrequency, (4.21) becomes iσ ˆ p c = w c (cid:104) ˆ U i − ˆ U t (cid:105) . (4.22)The same modeling approach is adopted for convenience at the junction where the con-servation of mass reads iσ ˆ p J = w J (cid:104) ˆ U R − ˆ U i + ˆ U t (cid:105) (4.23)with w J = γP /V J , where V J is the volume of the control volume modeling the junction(figure 8).Phase variations along the axial coordinate, x , are significant for the other componentsof the engine and the direct application of the complete set of linearized Euler equationsis necessary. In all cases, the fluctuating field and the base state, defined by ρ , T , and P , are assumed to be exclusively a function of x .6 Carlo Scalo, Sanjiva K. Lele, Lambertus Hesselink
For the resonator ( R ) and feedback inertance ( i ) isentropic wave propagation is as-sumed, yielding a simplified set of linearized equations for mass and momentum, iσ ˆ p = − ρ a A ( x ) d ˆ Udx (4.24 a ) iσ ˆ U = − A ( x ) ρ d ˆ pdx (4.24 b )valid for a variable cross-sectional area distribution A ( x ), where a = √ γ R T is thespeed of sound based on the base temperature. By introducing a spatial discretization,the set of equations (4.24) can be recast in algebraic form, as( iσ I − B R ) · u R = 0 (4.25)for the resonator, and ( iσ I − B i ) · u i = 0 (4.26)for the inertance, where I is the identity matrix, B is an operator discretizing the r.h.s.of (4.24), and u is a discrete collection of complex amplitudes for pressure and flowrate, specifically, u R = { ˆ p R , ˆ U R } for the resonator, and u i = { ˆ p i , ˆ U i } for the feedbackinertance.The systems of equations (4.25) and (4.26) are isolated component eigenvalue problems(with boundary conditions to be specified) and their resolution, in the context of asystem-wide linear stability analysis, is only meaningful if coupled with all of the othercomponents in the engine, as discussed in the following.The heat-transfer and drag in the pulse tube, and the presence of gradients of basedensity and temperature in the REG/HX unit, require variations of entropy to be ex-plicitly accounted for. Replacing the conservation equation for the total energy with thetransport equation for entropy, expressed in terms of temperature and pressure usingGibbs’ relation, yields iσA ˆ ρ + d ˆ Udx ρ + ˆ U dρ dx = 0 (4.27 a ) iσ ˆ U + Aρ d ˆ pdx = − R C ρ ˆ U (4.27 b ) ρ C p (cid:34) iσ ˆ T + ˆ UA dT dx (cid:35) − iσ ˆ p = − α T ˆ T (4.27 c )where the source term on the r.h.s of (4.27 b ) is obtained by linearizing the drag-model(2.3 a ), due to Organ (1992) and the heat-transfer model on the r.h.s of (4.27c) is the sameone used in (2.6) (Bejan 2004). Such terms are only activated in the heat-exchangers andthe regenerator. The spatial distribution of the base state in the pulse tube is taken fromthe numerical data at t = 0 .
55 s (figure 4a), at the early stages of the start-up phase(figure 7).Recasting the system of equations in (4.27) in diagonalized form yields iσ ˆ p = (cid:20) ρ R B ˆ T ˆ U − P A ddx − R T A dρ dx (cid:21) ˆ U + (cid:2) ρ R B ˆ T ˆ T (cid:3) ˆ T (4.28 a ) iσ ˆ U = (cid:20) − Aρ ddx (cid:21) ˆ p + (cid:20) − R c ρ (cid:21) ˆ U (4.28 b ) iσ ˆ T = B ˆ T ˆ U ˆ U + B ˆ T ˆ T ˆ T (4.28 c ) odeling of a Thermoacoustic Heat Engine B ˆ T ˆ U = − T A (cid:20) T d T dx + ( γ − ddx (cid:21) (4.29) B ˆ T ˆ T = − γ α T ρ C p (4.30)which, discretized in space, yields ( iσ I − B t ) · u t = 0 (4.31)where u t = { ˆ p t , ˆ U t , ˆ T t } .The complete eigenvalue problem can finally be solved by coupling the isolated com-ponent eigenvalue problems (4.22), (4.23), (4.25), (4.26), and (4.31) via the followingconditions ˆ U R = 0 (4.32 a ) ddx ˆ p R (cid:12)(cid:12)(cid:12)(cid:12) = 0 (4.32 b )ˆ p J = ˆ p R = ˆ p i = ˆ p t (4.32 c )ˆ p c = ˆ p i = ˆ p t (4.32 d )ˆ T t = 1 ρ R γ − γ ˆ p J (4.32 e )ˆ T t = 1 ρ R γ − γ ˆ p c (4.32 f )representing, respectively, the hard-wall condition on the left end of the resonator, conti-nuity of pressure at the junction and in the compliance, and an isentropic closure for thetemperature fluctuations at the two ends of the pulse tube. The final eigenvalue prob-lem can now be built by first combining (4.25), (4.26), and (4.31), into one system ofequations, iσ I − B t B R
00 0 B i · v = 0 , (4.33)where v = { u t ; u R ; u i } , and then incorporating the conditions (4.32) and the equations(4.22) and (4.23) to close the problem. Each of the conditions in (4.32), (4.22), (4.23)replaces one corresponding equation in (4.33), therefore, not affecting the rank of thesystem. The eigenvalue structure is finally recovered by absorbing the equations that donot contain σ (i.e. the ones deriving from (4.32)) via Gaussian elimination.The linear modeling framework composed of equations (4.22) and (4.23), (4.24), and(4.27) is first tested against the three-dimensional numerical simulation data for verifica-tion purposes. First, the acoustic impedance at different axial positions in the resonator,obtained by numerically integrating (4.24), has been quantitatively verified against thesimulation data, as well as the constants w c in (4.22) and w J in (4.23), which are equalto 135 × N/m and 312 × N/m , respectively. Good agreement is also obtainedby directly integrating (4.27) from section t to t and i to i (figure 8) using datafrom the numerical simulations as initial conditions (table 3). The integration has beencarried out with the exact value of the frequency and growth rate extracted from thesimulations (figure 9). Numerical trials have shown that, for a given i.c. and base state,the direct integration of (4.27) is much more sensitive to the angular frequency, ω than8 Carlo Scalo, Sanjiva K. Lele, Lambertus Hesselink
Figure 8.
Sketch illustrating the subdivision of the engine into control volumes representingdifferent components (cfr figure 1): compliance ( c ), pulse tube ( t ), feedback inertance ( i ), junction( J ), and resonator ( R ). Control surfaces are shown with dashed lines marking the beginning (0)and the end (1) of one-dimensional flow segments (oriented according to arrows) modeling thepulse tube, inertance and resonator. Complex pressure, volumetric flow rate and temperatureare indicated with ˆ p , ˆ U , and ˆ T respectively. References to the thermal buffer tube are droppedin the context of the analysis of the start-up phase since nonlinear effects can be neglected. Theleft end of the control volume representing the junction is located at x = -23 mm. Table 3.
Comparison between linear theory and simulation data for T h = 460K and T h = 500Kextracted at locations shown in figure 8 at t = 0 .
55 s. The linearized equations (4.27) areintegrated from t to t and i to i (figure 8) with initial conditions, at t and i , respectively,and α , ω , T ( x ) and ρ ( x ) taken from the numerical simulations (figures 4a,9). T h = 460 K T h = 500 K simulation data linear theory simulation data linear theory pulse tube: ˆ p t (1035.43 Pa, 0 ◦ ) used as i.c. (1566.54 Pa, 0 ◦ ) used as i.c.ˆ U t (0.00159 m /s, -39.0 ◦ ) used as i.c. (0.00229 m /s, -36.9 ◦ ) used as i.c.ˆ U t (0.00586 m /s, -68.6 ◦ ) (0.00651 m /s, -70.2 ◦ ) (0.00895 m /s, -66.9 ◦ ) (0.01046 m /s, -69.0 ◦ )ˆ p t (952.30 Pa, 3.28 ◦ ) (955.75 Pa, 3.11 ◦ ) (1437.71 Pa, 3.37 ◦ ) (1442.83 Pa, 3.03 ◦ ) inertance: ˆ U i (0.00349 m /s, 68.8 ◦ ) used as i.c. (0.00537 m /s, 69.0 ◦ ) used as i.c.ˆ p i (946.04 Pa, 3.01 ◦ ) used as i.c. (1427.41 Pa, 3.09 ◦ ) used as i.c.ˆ U i (0.00216 m /s, 55.0 ◦ ) (0.00218 m /s, 53.5 ◦ ) (0.00337 m /s, 56.2 ◦ ) (0.00337 m /s, 54.2 ◦ )ˆ p i (1031.70 Pa, 0.15 ◦ ) (1024.34 Pa, 0.07 ◦ ) (1561.24 Pa, 0.17 ◦ ) (1563.52 Pa, -0.23 ◦ ) in the growth rate α , which suggests that the prediction of the latter, in thermoacousticsystems, is potentially problematic, especially within the framework of linear modelingin the spectral domain.The eigenvalues σ are finally calculated by directly solving the eigenvalue problem ( ?? )for operating conditions ranging between τ = 1 .
25 and τ = 1 .
85, where τ = T h /T c isthe hot-to-cold temperature ratio. The segments representing the pulse tube, inertanceand resonator were discretized with 256, 64 and 32 points respectively, with a forth-orderspatial polynomial reconstruction. The frequency of instability and its variation with τ arepredicted within a ∼ ∼ odeling of a Thermoacoustic Heat Engine .
25 1 .
40 1 .
55 1 .
70 1 . τ f r e q u e n c y , ω / π ( / s ) .
25 1 .
40 1 .
55 1 .
70 1 . τ − . − . . . . . . g r o w t h r a t e , α ( / s ) − − li m i t - c y c l e p r e ss . a m p l . p l c ( k P a ) (a) (b) Figure 9.
Frequency, ω/ π (a), growth rate, α , and limit-cycle pressure amplitude p (cid:48) lc (b) versushot-to-cold temperature ratio, τ = T h /T c . Linear stability model (- - -), numerical simulations(symbols) with corresponding fitting (—) yielding critical temperature ratio τ cr = 1 .
505 and p (cid:48) lc | δτ =1 = 41 ,
000 Pa. in the complete variable area resonator alone (without the pulse tube), in accordancewith phase distribution shown in figure 6b, which is consistent with simple standingwave resonance. The growth rate is slightly over-predicted (figure 9b) having neglectedviscous and nonlinear losses. Overall, the quantitative agreement is very encouraging,serving both as a verification step for the full Navier-Stokes calculations and to gaininsight into the nature of the instability, also briefly discussed in the following section.The complete eigenvalue problem could also be solved by iteratively integrating (4.27)and (4.28) in space starting from a given set of initial conditions or guesses. This ap-proach is adopted in
DeltaEC (Ward & Swift 1994) to predict the limit-cycle pressureand velocity distributions in TAEs, which is a valid approximation in the case of rel-atively low pressure amplitudes, limited waveform distortion and simple geometries. Inprevious numerical trials this approach has proven to be unsuccessful in the context of thepresent device, especially in predicting the correct growth rate. A fully implicit spatialformulation, similar to a Helmholtz solvers used in reactive flows (Poinsot & Veynante2011), which has been adopted in the present context, is the only one that has proven tobe robust, cost effective and reliable.4.4.
Supercritical Hopf Bifurcation
A linear fit of the growth rates, α , extracted from the numerical simulation data versusthe temperature ratio τ = T h /T c (figure 9b) suggests a critical value of τ cr = 1 . p (cid:48) lc = p (cid:48) lc | δτ =1 √ τ − τ cr , (4.34)to the limit-cycle pressure amplitudes p (cid:48) lc versus τ . A similar result is obtained via nonlin-ear modeling by Mariappan & Sujith (2011). The fitting parameters in (4.34) are p (cid:48) lc | δτ =1 (dimensional) and τ cr (dimensionless). Moreover, two assumptions required by the Hopfbifurcation theorem are also satisfied: the non-hyperbolicity condition, α = 0 and ω (cid:54) = 0at τ cr (figure 9a), and the transversality condition, d α /d τ (cid:54) = 0 at τ cr (figure 9b). Theseresults have important implication on the parametrization of nonlinear fluxes, discussedlater in section 5.2.0 Carlo Scalo, Sanjiva K. Lele, Lambertus Hesselink
5. Nonlinear Regime
The analysis carried out so far has been exclusively based on the assumption of linearacoustic perturbations and therefore limited to the start-up phase. Nonlinear effects,however, are already detectable after only a few cycles of operation. These include thedeparture from exponential growth of the pressure amplitude (figure 7), the presence ofbroadband fine-scale flow structures associated with transitional turbulence, and a driftin the fluid parcels’ velocity, already noticable in the REG/HX during the startup phase(figure 4b). The latter phenomenon is known as acoustic streaming, which is the focuson this section.The most dramatic manifestation of acoustic streaming in traveling-wave TAEs isthe advective heat leakage from the hot heat exchanger (figure 10) which requires theintroduction of a secondary ambient heat exchanger (AHX2) (figure 1) to remove theexcess heat and achieve a limit cycle. Acoustic streaming occurs everywhere in the engineand its prediction and its suppression is one of the main technological challenges for thedesign of efficient TAEs.5.1.
Direct Modeling of Acoustic Streaming
A triple decomposition can be invoked to separate the streaming flow (rigorously definedbelow) from the acoustic field and the small-scale high-frequency fluctuations, startingwith the Reynolds decomposition ρ = ρ + ρ (cid:48) (5.1 a ) p = p + p (cid:48) (5.1 b ) u i = u ,i + u (cid:48) (5.1 c )where the subscript ‘0’ indicates a sharp-spectral-filtered quantity (also used before toindicate mean quantities), such as u ,i = u i = (cid:90) ∞−∞ u i ( x , t + τ ) sin( π f c τ ) π τ dτ (5.2)where f c is the cut-off frequency. The filtering operation (5.2) is, in practice, carriedout over 6 acoustic periods by adopting Simpson’s quadrature rule on the discrete datasampled at 2.2 kHz and f c = 0 . f where f = ω/ π is the acoustic frequency (figure 9a).The remainder of the filtering operation in (5.1) can be further decomposed into a purelyacoustic (subscript ‘ a ’) and a small-scale component (subscript ‘ t ’), ρ (cid:48) = ρ (cid:48) a + ρ (cid:48) t (5.3 a ) p (cid:48) = p (cid:48) a + p (cid:48) t (5.3 b ) u (cid:48) i = u (cid:48) a,i + u (cid:48) t,i (5.3 c )where the acoustic component can be isolated by applying another filtering operationthat removes frequencies higher than f while preserving the full acoustic amplitude. Dueto the truncation in time of the filter kernel, this was achieved, in practice, with a cut-offfrequency of f c > f .Vorticity contours from instantaneous visualizations (figure 10) suggest that the mostintense small-scale fluctuations occur in the feedback inertance. The unsteadiness ofthe (larger-scale) acoustic fluctuations does not allow turbulence to reach a fully (oreven partially) developed state. The Reynolds number based on the Stokes thickness δ ν = (cid:112) ν/ω and the maximum velocity amplitude at the center of the feedback inertanceis approximately Re δ ν = 480 (disturbed laminar regime, Jensen et al. (1989)), suggestingthat the turbulent kinetic energy generated from the break-up of the vortices rolling-up odeling of a Thermoacoustic Heat Engine Figure 10.
Instantaneous visualizations of temperature contours (colorbar in bottom figure)showing streaming of hot fluid in the thermal buffer tube and vorticity magnitude (white)showing intense vortex shedding and transitional turbulence. Data is shown over one completeacoustic cycle with 90 ◦ phase increments from (a) to (d). (Color online) Carlo Scalo, Sanjiva K. Lele, Lambertus Hesselink − −
250 0 250 500 750 1000 1250 1500 1750 2000 2250 2500 2750 3000 h u t,x u t,x i , h u a,x u a,x i (m /s ) r ( mm ) Figure 11.
Profiles of cycle- and time-averaged variance of acoustic (cid:104) u (cid:48) a,x u (cid:48) a,x (cid:105) (—)and small-scale (cid:104) u (cid:48) t,x u (cid:48) t,x (cid:105) axial velocity fluctuations (– –) extracted in the inertance, for x/L i = 0 . , . , . , . , . , .
99, from left to right, where L i = 0 .
204 m (shifted by500 m /s for clarity) for T h = 500 K and grid C (table 2). (a)(b) Figure 12.
Contours of axial component of time-averaged streaming velocity (cid:104) ˜ u x (cid:105) (5.7) for T h = 500K and grid C (table 2). Full-scale visualization (a) and zoom on the right end (b).Results have been mirrored about the centerline and streamlines have been only numericallyapproximated for illustrative purposes. (Color online) from the edges of the annular tube is not sustained. Moreover, turbulent stresses extractedfor the T h = 500 K case (highest drive ratio) and for the finest grid available (table 2)are approximately two orders of magnitude smaller than the acoustic stresses (figure 11)and will be neglected in the following analysis. This choice also accommodates the needto devise a simple predictive modeling framework for the streaming velocity (5.7), whichis discussed in the following.Substituting the decomposition (5.1) into the time-filtered conservation of mass, ig-noring temporal and spatial variations of the filtered density field (a strong assumption,particularly for the regions around the sharp edges and in the thermal buffer tube), as-suming that second-order quantities in the small-scale fluctuations are negligible withrespect to their acoustic counterpart (figure 11), u (cid:48) t,i u (cid:48) t,i << u (cid:48) a,i u (cid:48) a,i (5.4) ρ (cid:48) t u (cid:48) t,i << ρ (cid:48) a u (cid:48) a,i , (5.5) odeling of a Thermoacoustic Heat Engine ∂ ˜ u i ∂x i = 0 , (5.6)where ˜ u i is the density-weighted velocity field, ˜ u i = ρu i /ρ , defined based on the filteringoperation (5.2), which, under the assumptions made, can be expressed as˜ u i (cid:39) u ,i + ρ (cid:48) a u (cid:48) a i ρ . (5.7)In the present manuscript the density-weighted velocity field, ˜ u i , is adopted as the def-inition of the streaming velocity based on (5.7), which is second-order accurate in waveamplitude.The streaming velocity field in our case is axially-symmetric and spans the full extentof the engine (figure 12). Very large and elongated recirculations, of the order of a quarterof the acoustic wavelength, are visible in the resonator and are driven by the wall-normalgradient of the wave-induced shear stresses (not shown). Large recirculations of the orderof the resonator radius near the sharp edges of the pulse tube and a mean flow circulatingaround the pulse tube (following the direction of the amplified waves) are also observed.The latter is called Gedeon (or DC) streaming (Gedeon 1997) and is responsible forthe advective heat-leakage of the type shown in figure 10, which limits the efficiency ofmost traveling-wave thermoacoustic engines (as also discussed later in the context of thepresent engine).Assuming that time scales of variation of the filtered quantities ρ and u ,i are muchlonger than the acoustic period, and under the same assumptions underlying the deriva-tion of (5.7), it can be shown that ˜ u i satisfies the incompressible Navier-Stokes equations(Rudenko & Soluyan 1977), ∂ ˜ u i ∂t + ∂∂x j ˜ u i ˜ u j + ∂p ∂x i − ν ∇ ˜ u i = F a,i (5.8)where p = P /ρ and the forcing term F a,i is the divergence of the wave-inducedReynolds stresses, which can be expressed to second-order accuracy in wave-amplitudeas F a,i = − ∂∂x j u (cid:48) a,j u (cid:48) a,i + ∂∂x j (cid:26) − νρ (cid:20) ∂∂x j ρ (cid:48) a u (cid:48) a,i + ∂∂x i ρ (cid:48) a u (cid:48) a,j (cid:21) + 23 νρ ∂∂x k ρ (cid:48) a u (cid:48) a,k δ ij (cid:27) . (5.9)For large Reynolds numbers based on the streaming velocity, the terms containing themolecular diffusion in (5.9) can be neglected finally yielding, F a,i (cid:39) − ∂∂x j u (cid:48) a,j u (cid:48) a,i . (5.10)The maximum value of the stresses u (cid:48) a,j u (cid:48) a,i is expected to be found in the feedbackinertance where the acoustic power is maximum in the system. A grid-sensitivity studyon all of the components of the stresses in the inertance (figure 13) shows monotonicgrid-convergence of the stresses from grid A to C (table 2).The direct evaluation of (5.10) from the numerical data reveals very high values of F a,i near the sharp edges of the annular tube (figure 14) which locally drive the largeaforementioned recirculations. On the other hand, Gedeon Streaming is driven by theviscous decay of the wave amplitude in the annular inertance. This results in a negativeaxial gradient of normal stress u (cid:48) a,x u (cid:48) a,x visible in both figure 14 and 13b.4 Carlo Scalo, Sanjiva K. Lele, Lambertus Hesselink − −
10 0 10 20 30 40 50 60 70 80 90 100 110 120 h ˜ u x i (m/s) r ( mm ) (a) − −
10 0 10 20 30 40 50 60 70 80 90 100 110 120 q h u a,x u a,x i (m/s) r ( mm ) (b) − −
10 0 10 20 30 40 50 60 70 80 90 100 110 120 q h u a,r u a,r i (m/s) r ( mm ) (c) − −
75 0 75 150 225 300 375 450 525 600 675 750 825 900 h u a,x u a,r i (m /s ) r ( mm ) (d) Figure 13.
Profiles of time-averaged streaming velocity (a) and wave-induced Reynolds stresses(b),(c),(d) in the inertance, for x/L i = 0 . , . , . , . , . , .
99, respectively from left toright, where L i = 0 .
204 m (shifted by 20 m/s or 150 m /s for clarity). Data for T h = 500 K andgrid A (– · –), grid B (– –) and grid C (—) (table 2). The time-averaged volumetric flux throughthe inertance (intensity of Gedeon streaming) for this drive ratio is approximately 0.0033 m /s,corresponding to a mean streaming velocity of 1.2 m/s. odeling of a Thermoacoustic Heat Engine (a)(b) Figure 14.
Contour plots of the divergence of wave-induced Reynolds stresses. Axial, F a,x (a)and radial, F a,r (b) components extracted for T h = 500K on grid C (table 2). An axially symmetric numerical model,
Stream X (for more details see Appendix A),has been developed to directly simulate the streaming velocity field as the solution of theincompressible equations (5.8) driven by the divergence of the wave-induced stresses (fig-ure 14) extracted from the three-dimensional fully compressible calculations. Secondaryfeatures such as steady large-scale recirculations near the sharp edges of the annular tubeare only qualitatively reproduced (figure 15) with the appearance of a second recircu-lation in the compliance, which is not observed in the calculations. The actual targetof the present low-order modeling effort is the prediction of the intensity of the Gedeonstreaming. In spite of numerical challenges involved in solving incompressible flow in thepresence of a sharp edge, the latter is predicted fairly accurately (figure 17a), especiallyfor low drive ratios. For high drive ratios, resulting in very high limit-cycle acoustic am-plitudes, the errors associated with the assumptions made in deriving (5.7), (5.8) and(5.10) become too severe. ‘Slow’ streaming (Rudenko & Soluyan 1977) never actuallyoccurs in our engine, where the maximum intensity of ˜ u is comparable to the acousticvelocity amplitude in all cases.5.2. Efficiency and Energy Fluxes in the Thermal Buffer Tube
As the acoustic energy grows during the initial transient, so does the intensity of theGedeon streaming, increasing the rate of advective transport of hot fluid away from theHHX towards the AHX2. This results in unwanted heat leakage, which lowers the overallenergy conversion efficiency. The gradual expansion of the gas in the TBT determinesa slow increase of the background pressure in the system, which stops only when thehot temperature front reaches the AHX2. At this point, a rapid increase of the growthrate is observed, as shown by the kink in the time series in figure 7. A limit cycleis only reached later, with a constant background pressure, when the acoustic energyproduction is balanced by the losses in the system, which include streaming in resonatorand dissipation associated with the turbulent vortex shedding from the pulse tube walls.The exact conservation equation for the density-averaged internal energy, ˜ e = ρ e/ρ ,6 Carlo Scalo, Sanjiva K. Lele, Lambertus Hesselink
Figure 15.
Positive (—) and negative (– –) iso-levels of Stokes-streamfunction from the incom-pressible solver
Stream X driven by the wave-induced Reynolds stresses (figure 14) extractedfrom the fully compressible three-dimensional calculations for T h = 500K on grid C (table 2). reads (Lele 1994) ∂∂t ( ρ ˜ e ) = − ∂∂x j (cid:16) ρ (cid:104) ˜ u j ˜ e + (cid:93) h (cid:48)(cid:48) u (cid:48)(cid:48) j (cid:105) − q j (cid:17) + ∂∂x j (cid:16) p (cid:48) u (cid:48)(cid:48) j (cid:17) + u (cid:48)(cid:48) i ∂p∂x i − p (cid:48) ∂u (cid:48)(cid:48) i ∂x i (5.11)where u (cid:48)(cid:48) j = u i − ˜ u i and h (cid:48)(cid:48) = h i − ˜ h are the fluctuations of the density-weighted averagesof velocity and enthalpy, and q j is the time filtered molecular heat flux. Applying (5.11)to the flow in the TBT approximated as quasi one-dimensional, neglecting small termsand assuming equilibrium conditions, yields (cid:104) ρ ˜ u ˜ e (cid:105) (cid:124) (cid:123)(cid:122) (cid:125) Advective H.T. + (cid:104) ρ (cid:93) h (cid:48)(cid:48) u (cid:48)(cid:48) (cid:105) (cid:124) (cid:123)(cid:122) (cid:125) Thermoacoustic H.T. − (cid:104) p (cid:48) u (cid:48)(cid:48) (cid:105) (cid:124) (cid:123)(cid:122) (cid:125) Acoustic Energy Flux (cid:39) const (5.12)which is verified with fairly good approximation in the simulations (figure 16b). Theintensity of the advective heat transport in the TBT is proportional to the mean tem-perature profile, since the streaming velocity is uniform in this region (figure 12). Theadjustment length back to ambient temperature of the mean temperature distributionincreases with the drive ratio (figure 16a). This is due to thermoacoustic heat transportmechanisms. As the hot fluid front is transported with stronger intensity towards theAHX2 by the high-amplitude velocity fluctuations, steeper instantaneous temperaturegradients form at the interface between the AHX2 and the TBT (figure 10). The resultis a net cycle-average conductive heat flux in the positive axial direction creating a tem-perature buffer region. The intensity of the conductive heat flux in (5.11) is, however,neglibile compared to the quantitites in (5.12), which dominate the energy transportbudget in the TBT.An accurate evaluation of the other terms in (5.11) is, unfortunately, made impracticalby the (necessary) application of a second-order ENO reconstruction in the TBT andthe smoothing associated with azimuthally averaging the three-dimensional unstructureddata. The energy balance expressed by (5.12) is, however, satisfied to a sufficient degreeof accuracy to gain insight into the role of the Gedeon streaming in determining theoverall efficiency of the device and the scaling of the energy fluxes in (5.12).While no direct energy extraction component (e.g. an acoustic load such as a piezoelectric element or a linear alternator) has been included in the setup investigated, a odeling of a Thermoacoustic Heat Engine A H X RE G HH X A H X ∼ TBT ∼ .
00 0 .
05 0 .
10 0 .
15 0 . h T i ( K ) .
00 0 .
05 0 .
10 0 .
15 0 . x (m) − − − p o w e r ( k W ) (a)(b) Figure 16.
Axial distribution of mean temperature for T h = 460K, T h = 480K and T h = 500K(a) and surface integrated energy fluxes in (5.12) only for grid C and T h = 500 K at limit cycle.Acoustic power (- -), thermoacoustic heat transport, (– –), advective heat transport (– · –), andoverall sum ( ◦ ) (b). .
50 1 .
55 1 .
60 1 .
65 1 . τ . . . . . . G e d eo n s t r e a m i n g ( m / s ) .
50 1 .
55 1 .
60 1 .
65 1 . τ − . − . − . − . . . p o w e r ( k W ) . . . . . . η H (a) (b) Figure 17.
Intensity of Gedeon streaming versus hot-to-cold temperature ratio for fully com-pressible Navier-Stokes simulations ( (cid:3) ) with corresponding fit ( ), incompressible model( (cid:4) ), order-of-magnitude analysis (5.16) ( (cid:46) ), efficiency η H (5.13) ( (cid:63) ). (a); average nonlinear en-ergy fluxes at limit cycle, full Navier-Stokes simulations (symbols) and corresponding fit basedon (5.18), (5.20), (5.22) (lines), respectively for advective heat transport ( (cid:68) ,– · –), acoustic power( (cid:52) ,---), and thermoacoustic heat transport ( (cid:3) ,– –). Carlo Scalo, Sanjiva K. Lele, Lambertus Hesselink metric for efficiency can still be defined as η H = (cid:104) p (cid:48) u (cid:48)(cid:48) (cid:105)(cid:104) ρ (cid:93) h (cid:48)(cid:48) u (cid:48)(cid:48) (cid:105) + (cid:104) ρ ˜ u ˜ e (cid:105) , (5.13)which is the ratio between acoustic energy produced and total heat dissipated by thesecond-ambient heat exchanger. For example, our theoretical device produces ∼ T h = 500 K while losing ∼ ∼ > >
30% ifbuilt in a cascaded configuration (Gardner & Swift 2003). The efficiency directly evalu-ated from the simulation data (figure 17a) decreases rapidly with the temperature ratio,which could be expected in thermoacoustic engines with excessive Gedeon streamingcontrolling the energy balance in the TBT (G. Swift, pers. comm., 2014). However, theuncertainties in the integral quantities in figure 17b, due to averaging over a limitednumber (approximately 25) of acoustic cycles, makes the direct metric for the efficiency(5.13) not very reliable.A more robust estimate for η H can be derived by investigating the scaling of thevolume-averaged energy fluxes in (5.12) (figure 17b) with a simple order-of-magnitudeanalysis and curve fitting. The advective flux in the TBT, driven by the Gedeon stream-ing, is expected to scale as (cid:104) ρ ˜ u ˜ e (cid:105) ∼ ρ ,tbt u dc C v T h , (5.14)where ρ ,tbt ∼ p / ( R T c τ ) is the average density in the TBT and the hot temperature issimply T h = T c τ . The analytical expression of the Stokes drift due to a freely propagatingtraveling wave suggests that the intensity of the streaming velocity scales as u dc ∼ | u (cid:48) lc | a , (5.15)where u (cid:48) lc ∼ p (cid:48) lc / ( ρ a ) is the limit-cycle velocity amplitude. While the quadratic scalingof the streaming velocity is an expected result, a more quantitative estimate for u dc can beobtained by further simplifying the analysis in section 5.1, where the streaming velocityis directly modeled based on an hydrodynamic analogy. In fact, by roughly measuring theaverage spatial decay rate of the axial acoustic stresses in the inertance, (cid:104) F a,x (cid:105) i (figure18,top), and equating that to the linearized viscous losses in the pulse tube (see (2.3a)),the intensity of the Gedeon streaming can be estimated as u dc ∼ ρ A i L i (cid:104) F a,x (cid:105) i V AHX2 R c, AHX2 + V HHX R c, HHX + V REG R c, REG + V AHX R c, AHX , (5.16)where V and R c are, respectively, the volume and drag coefficients (obtained by lineariz-ing (2.3)) of the heat-exchanger/regenerators in the pulse tube. The estimate (5.16) isin good quantitative agreement with the results from section 5.1 (figure 17a). The casefor T h = 480K is an outlier simply due to the lack of grid convergence of the acousticstresses for this particular case (table 2). The quadratic scaling adopted in (5.15) is,therefore, further justified by (5.16) where (cid:104) F a,x (cid:105) i ∼ u (cid:48) lc . By invoking the previouslyderived scaling for the limit-cycle pressure (4.34), (5.15) becomes u dc = A dc
12 1 ρ a p (cid:48) lc (cid:12)(cid:12)(cid:12) δτ =1 ( τ − τ cr ) (5.17)where the fitting coefficient is A dc = 0 .
52 (figure 17b). Substituting (5.17) into (5.14) odeling of a Thermoacoustic Heat Engine (cid:104) ρ ˜ u ˜ e (cid:105) = A adv (cid:20) C v P R A dc
12 1 ρ a p (cid:48) lc (cid:12)(cid:12)(cid:12) δτ =1 (cid:21) ( τ − τ cr ) , (5.18)where the fitting coefficient is A adv = 0 .
51. The same procedure can be applied to theacoustic energy flux, which scales as (cid:104) p (cid:48) u (cid:48)(cid:48) (cid:105) ∼ p (cid:48) lc u (cid:48) lc
12 cos(∆ φ ) (5.19)where ∆ φ is the phase difference between pressure and velocity observed in the REG/HX(figure 4c), leading to (cid:104) p (cid:48) u (cid:48)(cid:48) (cid:105) = A ap
12 cos(∆ φ ) 1 ρ a p (cid:48) lc (cid:12)(cid:12)(cid:12) δτ =1 ( τ − τ cr ) (5.20)where the fitting constant is A ap = 0 .
3. Finally, the thermoacoustic heat flux is expectedto scale as (cid:104) ρ (cid:93) h (cid:48)(cid:48) u (cid:48)(cid:48) (cid:105) ∼ ρ ,tbt C p ∂T∂P p (cid:48) lc p (cid:48) lc ρ a
12 cos(∆ φ ) (5.21)where ∆ φ is needed to account for the correlation between pressure and velocity (like in(5.19)), resulting in (cid:104) ρ (cid:93) h (cid:48)(cid:48) u (cid:48)(cid:48) (cid:105) = A ta (cid:20) ρ a p (cid:48) lc (cid:12)(cid:12)(cid:12) δτ =1
12 cos(∆ φ ) (cid:21) ( τ − τ cr ) (5.22)with the fitting coefficient A ta = 0 .
93. With all fluxes scaling as ∼ ( τ − τ cr ), this leadsto a robust estimate for the efficiency, η H, avrg = 0 .
13 (5.23)effectively averaged over the range of temperature ratios investigated.Linear acoustic solvers applied at the limit cycle can directly estimate second-orderquantities in the acoustic amplitudes such as (5.19) and (5.22), or even the mean wave-amplitude decay (cid:104) F a,x (cid:105) i due to viscous losses in a duct, without prior knowledge of thecritical temperature ratio, τ cr . However, process-based parametrizations for the Gedeonstreaming, and therefore the advective transport (5.18), are still currently missing despitehaving a first-order impact on the acoustic energy budgets and the efficiency. We haveshown that, for slow streaming,commonly found in realistic traveling-wave thermoacous-tic engines, a hydrodynamic analogy (Lighthill 1978) can be invoked and, further sim-plified leading to a very low-order modeling approach (5.16), which is very amenable inthe context of simple linear solvers used for engineering prediction of TAEs. Moreover,estimates such as (5.16) do not require the knowledge of the critical temperature ratio τ cr , making the hydrodynamic analogy investigated in section 5.1 an attractive modelingparadigm for streaming.
6. Conclusions
We have carried out three-dimensional numerical simulations of a theoretical traveling-wave thermoacoustic heat-engine (TAE). This is the first step in a broader researcheffort aimed at building multi-fidelity, full-scale prediction tools for TAEs. The goal isto assist technological design by directly simulating, under realistic operating conditions,the physical processes controlling the overall efficiency of such devices. These include thethermoacoustic instability, wave propagation and amplification in the startup phase, the0
Carlo Scalo, Sanjiva K. Lele, Lambertus Hesselink x (mm) h u a , x u a , x i ( m / s ) (c)(b)(a) Figure 18.
Illustration of balance (5.16) between traveling-wave decay in feedback inertanceand drag due to heat-exchangers and regenerator (top). Streamwise profile of axial wave-inducedstresses extracted at r = 56 mm for T h = 460K (a), T h = 480K (b), and T h = 500K (c) withlinear fitting (- -) approximating the mean decay rate. nonlinear effects at the limit cycle (mainly acoustic streaming and turbulence) and theeffects of geometrical complexities. The last two are not directly captured in state-of-the-art predictive tools for TAE.Inspired by the work of `a Nijeholt et al. (2005), we have devised a simple traveling-waveTAE model that could serve as a benchmark case for high-fidelity numerical simulationsof similar devices. We have extended such setup to three-dimensions and introduced asecond ambient heat-exchanger to achieve a limit cycle, modeling typical fluid dynamicconditions found in thermal buffer tubes. Details omitted in `a Nijeholt et al. (2005)regarding the geometry and the modeling of the heat-exchangers and regenerators havebeen reconstructed to the best of authors’ ability and reported in detail. In spite of itstheoretical nature, the model retains all the essential critical components, features andcomplexities of real traveling-wave TAEs.The time integration is carried out from initial quiescent conditions to the limit cy-cle. It is shown that the mechanisms responsible for the acoustic energy generationand propagation in the system during the start-up phase can be explained with lin-ear acoustics, despite the high amplitude ( ∼
1% of the mean) of the initial perturbationresulting from the activation of the source terms modeling the heat transfer. An ana-lytical linear Lagrangian model shows that the thermoacoustic instability occurring inthe regenerator/heat-exchanger (REG/HX) unit intensifies plane waves traveling in thedirection of the imposed temperature gradient via a process resembling a thermodynamicStirling cycle. The result is the establishment of a network of self-amplifying travelingwaves looping around the REG/HX unit. A system-wide linear stability model based onRott’s theory accurately predicts the frequency of the (only) unstable mode as well asthe the critical temperature ratio, despite not accounting for viscous and other nonlinear odeling of a Thermoacoustic Heat Engine
Carlo Scalo, Sanjiva K. Lele, Lambertus Hesselink
Appendix A.
Stream X : an axially symmetric incompressible flowsolver model Stream X solves the incompressible Navier-Stokes equations in cylindrical coordinates, ∂u x ∂t + 1 r ∂∂r ( ru r u x ) + ∂∂x u x = − ∂p∂x + ν (cid:20) r ∂∂r (cid:18) r ∂u x ∂r (cid:19) + ∂ u x ∂x (cid:21) + F x (A 1) ∂u r ∂t + 1 r ∂∂r (cid:0) r u r (cid:1) + ∂∂x u x u r = − ∂p∂r + ν (cid:20) r ∂∂r (cid:18) r ∂u r ∂r (cid:19) + ∂ u r ∂r − u r r (cid:21) + F r (A 2)where p = P/ρ and ( F x , F r ) is a body force, relying on a Stokes-streamfunction -vorticity, Ψ - ζ , formulation. A stagghered collocation for u x and u r , and nodal for Ψand ζ is adopted (figure 19a). Spatial derivatives are approximated with a second-ordercentral-difference scheme.The non-simply connected computational domain (figure 19a) requires a special time-advancement strategy to directly solve for (cid:126) Ψ n +1 = (0 , , Ψ n +1 ), given its unknown valueat the boundary ∂ Ω (while the value on ∂ Ω is fixed and arbitrary). The velocity fieldis first predicted at the next time step ( u ∗ , v ∗ ) with an explicit Runge Kutta integrationof the r.h.s. of (A 1), carried out without the pressure terms. This guarantees the exactprediction of vorticity and circulation at time t n +1 , respectively (cid:126)ζ n +1 = (cid:126) ∇ × ( u ∗ x , u ∗ r ) (A 3)Γ n +1 = (cid:73) ( u ∗ x , u ∗ r ) ˙ d(cid:126)l. (A 4)Knowing the exact boundary conditions for (cid:126) Ψ n +1 on ∂ Ω and ∂ Ω would allow to directlysolve (cid:126) ∇ × (cid:32) (cid:126) ∇ × (cid:126) Ψ n +1 r (cid:33) = (cid:126)ζ n +1 (A 5)and complete the time advancement yielding ( u n +1 x , u n +1 r ). A straightforwad workaroundis to express the Stokes-streamfunction at time t n +1 as the linear combinationΨ n +1 = Ψ A + α Ψ B , (A 6)where α is an unknown coefficient, and Ψ A and Ψ B are the solutions to (cid:126) ∇ × (cid:16) (cid:126) ∇ × (cid:126) Ψ B /r (cid:17) = 0 , for Ψ B | ∂ Ω = 0 , Ψ B | ∂ Ω = 1 (A 7) (cid:126) ∇ × (cid:16) (cid:126) ∇ × (cid:126) Ψ A /r (cid:17) = (cid:126)ζ n +1 , for Ψ A | ∂ Ω = 0 , Ψ A | ∂ Ω = 0 (A 8)this allows to calculate the circulations Γ A and Γ B and α = (Γ − Γ A ) / Γ B , (A 9)where Γ is known from (A 3). While Ψ A can be calculated in preprocessing, Ψ B needs tobe re-evaluated at every time step. Finally, Ψ n +1 is calculated from (A 6). Acknowledgments
The authors acknowledge the support of the Precourt Institute for Energy Seed Grantat Stanford and the computational time provided by the NSF-MRI grant on the StanfordCertainty cluster. The authors thank Dr. Gregory Swift for his very useful comments onthe draft and acknowledge the help of Prof. Ray Hixon for assisting in the development of odeling of a Thermoacoustic Heat Engine (a) (b) Figure 19.
Illustration of domain with same degree of connectivity as computational domain infigure 1 with ∂ Ω and ∂ Ω representing the resonator walls and the annular tube, respectively,and Γ the clock-wise circulation calculated around ∂ Ω (a). Variable collocation in Stream X foraxial velocity u ij , radial velocity v ij , vorticity ζ ij and Stokes-streamfunction Ψ ij (b). the source terms to model the heat-transfer and drag in heat-exchanger and regenerators.Carlo Scalo would like to thank, in particular, Dr. Julien Bodart and Dr. Ivan Bermejo-Moreno for their precious technical help and Jeffrey Lin providing the authors with acomplete DeltaEC model of the engine, which has lead to the choice of introducing asecondary ambient heat exchanger.
REFERENCESBackhaus, S & Swift, G. W.
J.Acoust. Soc. Am. , 3148–3166.
Bauwens, L
J. Fluid Mech. (135–161).
Bejan, A.
Convective Heat Transfer . John Wiley & Sons Inc.
Boluriaan, S. & Morris, P. J. . Boluriaan, S. & Morris, P. J.
Intl. Journalof Aeroacoustics (3-4), 255–292. Ceperley, P.H.
J. Acoust.Soc. Am. , 1239–1244. Gardner, D.L. & Swift, G. W.
J. Acoust. Soc. Am. (4).
Garret, S.L.
Am. J.Phys. , 11–17. Gedeon, D.
Cryocoolers , 385–392. Ham, F., Mattsson, K., Iaccarino, G. & Moin, P.
Towards Time-Stable and AccurateLES on Unstructured Grids , Lecture Notes in Computational Science and Engineering ,vol. 56, pp. 235 – 249. Springer Berlin Heidelberg.
Hamilton, M. F., Ilinksii, Yu. A. & Zabolotskaya, E. A.
J. Acoust. Soc. Am. (5).
Hireche, O., Weisman, C., Baltean-Carles, D., Quere, P. L., Francois, M. & Bauwens,L.
Comptes Rendus Mecanique (1),18–23.
In ’T Panhuis, P. H. M. W., Reinstra, S. W., Molenaar, J. & Slot, J. J. M.
J.Fluid Mech. , 41–70.
Jensen, B. L., Sumer, B. M. & Fredsøe, J.
Journal of Fluid Mechanics , 265–297. Carlo Scalo, Sanjiva K. Lele, Lambertus Hesselink
Karypis, G. & Kumar, V.
Proc. Supercomputing ’98 . Kirchhoff, G.
Pogg. Ann. (177).
Kramers, H. A.
Physica (971). Lele, Sanjiva K.
Annu. Rev. Fluid Mech. . Lighthill, J.
J. Sound Vib (3), 391–418. Mariappan, S. & Sujith, R. I.
J. Fluid. Mech. , 511 – 533.
M¨uller, U. A. & Rott, N.
Z. Angew. Math. Phys. (609). `a Nijeholt, J.A. Lycklama, Tijani, M.E.H. & Spoelstra, S. J. Acoust. Soc. Am. (4), 2265–2270.
Olson, J. R. & Swift, G. W.
Cryogenics , 769 – 776. Organ, A. J.
Thermodynamics and Gas Dynamics of the Stirling Cycle Machine . Penelet, G., Guedra, M., Gusev, V. & Devaux, T.
Int. J. Heat Mass Tran. , 6042–6053. Penelet, G., Gusev, V., Lotton, P. & Bruneau, M. a Experimental and theoreticalstudy of processes leading to steady-state sound in annular thermoacoustic engines.
Phys.Rev. E. (016625). Penelet, G., Gusev, V., Lotton, P. & Bruneau, M.
Phys. Lett. A ,268–273.
Penelet, G., Job, S., Gusev, V., Lotton, P. & Bruneau, M. b Dependance of soundsamplification on temperature distribution in annular thermacoustic engines.
Acta. Acust. , 567–577. Poinsot, T. & Veynante, D.
Theoretical and Numerical Combustion , third edition edn.R.T. Edwards, Inc.
Rott, N.
Z. Angew. Math. Phys. (230). Rott, N.
Z.Angew. Math. Phys. (54). Rott, N.
Z. Angew. Math. Phys. (417). Rott, N. a Thermally driven acoustic oscillations, part III: Second-order heat flux.
Z.Angew. Math. Phys. (43). Rott, N. b Thermally driven acoustic oscillations, part IV: Tubes with variable cross-section.
Z. Angew. Math. Phys. (43). Rott, N.
Neue Zurecher Ztg. (210).
Rott, N.
Adv. Appl. Mech. , 135–175.
Rott, N.
J. FluidMech. (1).
Rudenko, O. V. & Soluyan, S. I.
Schloegel, K., Karypis, G. & Kumar, V.
Swift, G. W.
J. Acoust. Soc. Am. (4), 1145–1181. Swift, G. W.
J. Acoust. Soc.Am. (3). Swift, G. W. & Ward, W.C.
J. Thermophys.Heat Tr. (4). Thomas, B. & Pittman, D.
Energy Conversion odeling of a Thermoacoustic Heat Engine Engineering Conference and Exhibit, 2000. (IECEC) 35th Intersociety , , vol. 1, pp. 76–84vol.1.
Thompson, M. W., Atchley, A. A. & Maccarone, M. J.
J. Acoust. Soc. Am. pp. 1939–1849.
Tijani, M.E.H. & Spoelstra, S.
J. Appl.Phys. , 093519. de Waele, A.T.A.M.
J. Sound Vib. , 974–988.
Ward, W.C. & Swift, G. W.
J. Acoust.Soc. Am. , 3671. Zouzoulas, G. & Rott, N.
Z. Angew. Math. Phys.27