2.1 Engine model design

The TAP engine model (figure 1) is 510 mm in length and is divided into two cylindrical, constant-area sections: one of 19.5 mm in diameter, enclosing an axisymmetric thermoacoustic stack (table 1) and the other of 71 mm in diameter, capped by a piezoelectric diaphragm tuned to the thermoacoustically amplified mode (388 Hz) for maximization of acoustic energy extraction. The TAP engine design is inspired by the design and experimental work of Smoker et al. (Reference Smoker, Nouh, Aldraihem and Baz2012).

An axisymmetric model cannot account for three-dimensional flow effects. The scope of this study is instead focused on the accurate modelling of thermoacoustic acoustic energy production (§ § 3 and 4), nonlinear thermoacoustic transport (§ 5) and energy extraction (§ 7), for which three-dimensional flow effects are secondary. Moreover, at the highest acoustic amplitude achieved in the present computations ( $\simeq 6000~\text{Pa}$ ), the Stokes Reynolds numbers based on the maximum centreline velocity amplitude in the device (at approximately $x=245~\text{mm}$ ) is $Re_{\unicode[STIX]{x1D6FF}_{\unicode[STIX]{x1D708}}}<100$ , where $\unicode[STIX]{x1D6FF}_{\unicode[STIX]{x1D708}}$ is the Stokes boundary layer thickness (3.4), falling well within the fully laminar regime of oscillatory boundary layers (Jensen, Sumer & Fredsøe Reference Jensen, Sumer and Fredsøe1989). Even at significantly higher Reynolds numbers and acoustic amplitudes, such as the ones achieved in the three-dimensional calculations of a large travelling-wave engine by Scalo et al. (Reference Scalo, Lele and Hesselink2015b ), hydrodynamic nonlinearities – such as Reynolds stresses associated with transition turbulence – were found to be negligible with respect to acoustic nonlinearities such as streaming and thermoacoustic transport.

In the experiments by Smoker et al. (Reference Smoker, Nouh, Aldraihem and Baz2012), a square-weave mesh-screen regenerator is used with porosity and hydraulic radius of $\unicode[STIX]{x1D719}=0.25$ and $r_{h}=0.34~\text{ mm}$ , respectively. The regenerator is heated on one side (in the hot cavity) by a resistive filament sustaining a hot temperature of $T_{h}=790~\text{K}$ , without a cold heat exchanger on the opposing side (Nouh, personal communication). As a result, the mean axial temperature gradient weakens throughout the course of the experiment due to conduction in the metal and thermoacoustic transport in the pore volume.

The thermoacoustic stack in our theoretical axisymmetric TAP engine model is composed of coaxial cylindrical annuli (table 1) with a linear axial wall temperature profile (from $T_{h}$ to $T_{c}$ ) imposed via isothermal boundary conditions. This choice allows for the direct application of Rott’s theory for verification of growth rates and frequencies observed during the start-up phase of the Navier–Stokes simulations, along with a clear definition of the geometrical parameter space for the exploration of the optimal stack design.

Figure 1. Illustration of the axisymmetric TAP engine model (not to scale) inspired by the experimental work of Smoker et al. (Reference Smoker, Nouh, Aldraihem and Baz2012). All lengths are given in millimetres. The dashed lines in the hot cavity indicate the original experimental design. The axial distribution of mean temperature, $T_{0}(x)$ , is qualitatively sketched. Different stack designs (table 1) and temperature settings (table 3) have been considered in the simulations.

Three stack configurations have been investigated (table 1) and were obtained by varying two parameters: (i) the number of coaxial solid annuli, $n_{s}$ , surrounding a central solid rod of $h_{s}/2$ in radius; and (ii) the ratio of the solid annulus thickness to the annular gap width, $h_{s}/h_{g}$ . Following the geometrical constraint

where $R_{stk}=19.5~\text{ mm}$ (radius of the small-diameter section enclosing the thermoacoustic stack), unique values for $h_{s}$ and $h_{g}$ were determined for given values of $n_{s}$ and $h_{s}/h_{g}$ . Volume porosity, $\unicode[STIX]{x1D719}$ , and hydraulic radius, $r_{h}$ , are calculated as

where $V_{g}$ is the total gas-filled volume in the stack, $S_{heat}$ is the gas–solid contact surface through which wall-heat transfer occurs, $A_{g}$ is the cross-sectional area available to the gas and $A_{g}+A_{s}=A_{stk}=\unicode[STIX]{x03C0}R_{stk}^{2}$ where $A_{s}$ is the cross-sectional area occupied by the solid.

Stack I has been designed by selecting a combination of $n_{s}$ and $h_{s}/h_{g}$ resulting in a porosity and hydraulic radius close to the values of the mesh-wire regenerator of  Smoker et al. (Reference Smoker, Nouh, Aldraihem and Baz2012). Stack II is characterized by a higher porosity with respect to Stack I, without significant differences in the hydraulic radius. Stack III has been designed by imposing $h_{s}=h_{g}$ and $n_{s}=3$ , resulting in a more porous and regularly spaced stack, and allowing for the formation of an inviscid acoustic core in the annular gap (missing in Stack I and II, see table 1), at the expense of thermal contact (see discussion in § 4.3).

Five different temperature settings have been considered in the Navier–Stokes simulations (table 3), bracketing values observed in the experiments (Smoker et al. Reference Smoker, Nouh, Aldraihem and Baz2012; Nouh, Aldraihem & Baz Reference Nouh, Aldraihem and Baz2014), ranging from a close-to-critical (case 1) to a very strong thermoacoustic response (case 5), the latter corresponding to a temperature gradient that might be challenging to sustain experimentally. A linear temperature profile, ranging from hot, $T_{h}$ , to cold, $T_{c}$ , is imposed on the thermoacoustic stack walls; no-slip isothermal boundary conditions corresponding to ambient conditions $T_{a}=300~\text{K}$ are imposed everywhere else, with the exception of the left and right end (including the piezoelectric diaphragm, if applicable), which are kept adiabatic.

2.2 Governing equations

The conservation equations for mass, momentum and total energy, solved in the fully compressible Navier–Stokes simulations of the TAP engine model are, respectively,

where $x_{1}$ , $x_{2}$ and $x_{3}$ (equivalently, $x$ , $y$ and $z$ ) are axial and cross-sectional coordinates, $u_{i}$ are the velocity components in each of those directions and $p$ , $\unicode[STIX]{x1D70C}$ and $E$ are respectively pressure, density and total energy per unit mass. The gas is assumed to be ideal, with equation of state $p=\unicode[STIX]{x1D70C}\,R_{gas}\,T$ and a constant ratio of specific heats, $\unicode[STIX]{x1D6FE}$ . The gas constant is fixed and calculated as $R_{gas}=p_{ref}(T_{ref}\,\unicode[STIX]{x1D70C}_{ref})^{-1}$ , based on the reference thermodynamic density, pressure and temperature, $\unicode[STIX]{x1D70C}_{ref}$ , $p_{ref}$ and $T_{ref}$ , respectively. The viscous and conductive heat fluxes are: where $\unicode[STIX]{x1D61A}_{ij}$ is the strain-rate tensor, given by $\unicode[STIX]{x1D61A}_{ij}=(1/2)(\unicode[STIX]{x2202}u_{j}/\unicode[STIX]{x2202}x_{i}+\unicode[STIX]{x2202}u_{i}/\unicode[STIX]{x2202}x_{j})$ ; $Pr$ is the Prandtl number; and $\unicode[STIX]{x1D707}$ is the dynamic viscosity, given by $\unicode[STIX]{x1D707}=\unicode[STIX]{x1D707}_{ref}(T/T_{ref})^{n}$ , where $n$ is the viscosity power-law exponent and $\unicode[STIX]{x1D707}_{ref}$ is the reference viscosity. Simulations have been carried out with the following gas properties: $\unicode[STIX]{x1D6FE}=1.4$ , $\unicode[STIX]{x1D70C}_{ref}=1.2~\text{kg}~\text{m}^{-3}$ , $p_{ref}=101\,325~\text{Pa}$ , $T_{ref}=300~\text{K}$ , $\unicode[STIX]{x1D707}_{ref}=1.98\times 10^{-5}~\text{kg}~\text{m}^{-1}~\text{s}^{-1}$ , $Pr=0.72$ , and $n=0.76$ , valid for air (De-Yi & Bu-Xuan Reference De-Yi and Bu-Xuan1990).

No-slip and isothermal boundary conditions are used on all axial boundaries in the model. Direct acoustic energy extraction is only allowed from the piezoelectric diaphragm (figure 1), modelled with impedance boundary conditions

formulated in the time domain following the numerical implementation by Scalo et al. (Reference Scalo, Bodart and Lele2015a ), summarized in appendix B. The broadband (dimensional) impedance $Z(\unicode[STIX]{x1D714})$ is derived by collapsing the experimentally determined two-port electromechanical admittance matrix for the piezoelectric element and fitting the resulting impedance with a multi-oscillator approach (Fung & Ju Reference Fung and Ju2001) as discussed in detail in § 6. In the present work, the characteristic specific acoustic impedance

is absorbed within the value of the impedance in (2.5); hence, (2.5) is treated as dimensional in implementing both single- and multi-oscillator impedance boundary conditions (§ 6). As described in appendix B, the impedance boundary conditions (2.5) are implemented via imposition of the complex wall softness coefficient, $\widehat{\widetilde{W}}_{\unicode[STIX]{x1D714}}$ , defined as

which is related to the complex reflection coefficient, $\widehat{W}_{\unicode[STIX]{x1D714}}$ , via

Hard-wall (purely reflective) conditions correspond to the limit of infinite impedance magnitude $|Z|\rightarrow \infty$ and therefore can be imposed by setting $\widehat{\widetilde{W}}_{\unicode[STIX]{x1D714}}=0$ . The subscript $\unicode[STIX]{x1D714}$ is introduced to avoid ambiguity when (2.7) and (2.8) are extended to the Laplace space, via the transformation $s=i\,\unicode[STIX]{x1D714}$ , yielding

It is important to stress that $\widehat{\widetilde{W}}_{\unicode[STIX]{x1D714}}(\cdot )$ and $\widehat{\widetilde{W}}_{s}(\cdot )$ are different functional forms, and the latter is convenient for the implementation of the TDIBC, as in § 6.

2.3 Computational set-up

The three different stack types (I, II and III) required different computational grids. For stack type I (figure 2), three different levels of grid resolution were considered (A, B and C, from coarse to fine). Stack types II and III were only meshed at the highest grid resolution level (table 2). Simulations with temperature settings 1–4 (table 3) have been performed on the finest grid resolution level C and only for stack type I. The viscous and thermal Stokes thicknesses at 300 K and 388 Hz (frequency of the thermoacoustically amplified mode) are $\unicode[STIX]{x1D6FF}_{\unicode[STIX]{x1D708}}\sim 0.12~\text{mm}$ and $\unicode[STIX]{x1D6FF}_{\unicode[STIX]{x1D705}}\sim 0.14~\text{mm}$ , respectively, and are resolved on all grids considered. The coarsest near-wall grid resolution considered is $\unicode[STIX]{x0394}r_{w}=0.06~\text{mm}$ (table 2). While the full three-dimensional Navier–Stokes equations are solved, azimuthal gradients are not captured on the adopted computational grid (figure 2), which is extruded azimuthally with $1^{\circ }$ increments for a total of 5 cells, with rotational periodicity imposed on the lateral faces. The results from the numerical computations are, in practice, axisymmetric.

Figure 2. Computational grid for resolution/stack-type A/I (see tables 2 and 3).

The governing equations are solved using $\text{CharLES}^{X}$ , a control-volume-based, finite-volume solver for the fully compressible Navier–Stokes equations on unstructured grids, developed as a joint effort among researchers at Stanford University. $\text{CharLES}^{X}$ employs a three-stage, third-order Runge–Kutta time discretization and a grid-adaptive reconstruction strategy, blending a high-order polynomial interpolation with low-order upwind fluxes (Ham et al. Reference Ham, Mattsson, Iaccarino and Moin2007). The code is parallelized using the message passing interface protocol and highly scalable on a large number of processors (Bermejo-Moreno et al. Reference Bermejo-Moreno, Bodart, Larsson, Barney, Nichols and Jones2013).

3.1 Hot cavity, pulse tube and resonator

In the hot cavity, pulse tube and resonator, a constant axial mean temperature distribution is assumed (figure 1), yielding the linearized equations

enforcing the (combined) conservation of mass and energy (3.1a ) and momentum (3.1b ), respectively. In these sections, the total cross-sectional area corresponds to the area available to the gas, $A=A_{g}$ . The complex thermoviscous functions $f_{\unicode[STIX]{x1D708}}$ and $f_{\unicode[STIX]{x1D705}}$ in (3.1) are

where $\text{J}_{n}(\cdot )$ are Bessel functions of the first kind and $\unicode[STIX]{x1D702}$ is the dimensionless complex radial coordinate

where $\unicode[STIX]{x1D708}_{0}=\unicode[STIX]{x1D707}(T_{0})/\unicode[STIX]{x1D70C}_{0}$ is the kinematic viscosity based on mean values of density and temperature, and $\unicode[STIX]{x1D702}_{w}$ in (3.2) is the dimensionless coordinate (3.3) calculated at the radial location of the isothermal, no-slip wall. The viscous, $\unicode[STIX]{x1D6FF}_{\unicode[STIX]{x1D708}}$ , and thermal, $\unicode[STIX]{x1D6FF}_{\unicode[STIX]{x1D705}}$ , Stokes thicknesses are

and are related via the Prandtl number, $\unicode[STIX]{x1D6FF}_{\unicode[STIX]{x1D708}}=\sqrt{Pr}\;\unicode[STIX]{x1D6FF}_{\unicode[STIX]{x1D705}}$ . The effective laminar boundary layer thickness is approximately three times the Stokes thickness. For a Prandtl number below unity, $Pr<1$ , the thermal boundary layer is thicker than the viscous layer.

3.2 Thermoacoustic stack

The analytical expression for the radial profile of the complex axial velocity amplitude within the $m\text{th}$ annular gap of the stack (table 1) has been derived for a generic axial location $x$ by neglecting radial variations of pressure, i.e. $\hat{p}^{(m)}(x,r)=\hat{p}^{(m)}(x)$ (appendix A), yielding

where

is the inviscid acoustic velocity, which varies with the axial direction $x$ ; $H_{n}^{(1)}(\cdot )$ are Hankel functions of the first kind; and

Rott’s wave equations can be written for the $m\text{th}$ annular flow passage of cross-sectional area $A_{g}^{(m)}$ (where $m\in \{1,\ldots ,n_{s}+1\}$ ), in the diagonalized form:

where the complex thermoviscous functions, $f_{\unicode[STIX]{x1D708}}^{(m)}$  and $f_{\unicode[STIX]{x1D705}}^{(m)}$ , are, in this case (appendix A)

Assuming that all annular flow passages share the same instantaneous pressure field, $\hat{p}^{(m)}(x)=\hat{p}(x)$ , and mean density and temperature axial distribution, and considering that the thermoviscous functions $f_{\unicode[STIX]{x1D708}}^{(m)}$ and $f_{\unicode[STIX]{x1D705}}^{(m)}$ differ at most by $2\,\%$ over all values of $m$ , it is possible to collapse the $n_{s}+1$ equations in (3.8a ) via area-weighted averaging and to take the arithmetic sum of (3.8b ) over $m$ , yielding a new set of approximate wave equations for the thermoacoustic stack,

where the total cross-sectional area available to the gas, $A_{g}$ , and flow rate, $\hat{U}$ , are

and an area-weighted equipartitioning of the flow rates, $\hat{U} ^{(m)}=A_{g}^{(m)}/A_{g}\;\hat{U}$ , has been assumed in (3.10a ).

3.3 Discretization, boundary and inter-segment conditions

Isolated-component eigenvalue problems for the cavity ( $c$ ), thermoacoustic stack ( $s$ ), pulse tube ( $t$ ) and resonator ( $r$ ) control volumes (figure 3) are first assembled in the form

with $\boldsymbol{v}=\{\boldsymbol{u}_{c},\boldsymbol{u}_{s},\boldsymbol{u}_{t},\boldsymbol{u}_{r}\}$ where $\boldsymbol{u}_{l}=\{\hat{\boldsymbol{p}}_{l},\hat{\boldsymbol{U}}_{l}\}$ is the collection of the discrete complex amplitudes of pressure and flow rate for the $l\text{th}$ segment where $l\in \{c,s,t,r\}$ , $\unicode[STIX]{x1D644}$ is the identity matrix and $\boldsymbol{B}$ is an operator discretizing the equations (3.1) for the hot cavity, pulse tube and resonator, and (3.10) for the thermoacoustic stack.

Equations for each segment are discretized on a staggered, uniform grid (figure 3). The mass and energy equation is written at each pressure node while the one for momentum is written at each flow rate location, both with a second-order central discretization scheme. Inter-segment conditions are

where subscripts $-2$ and $-1$ indicate the second last and last point of a segment, subscripts $+1$ and $+2$ indicate the first and second point of the following one and $\unicode[STIX]{x0394}\hat{p}_{ml}$ is the pressure drop due to minor losses (7.8), if applied. The extrapolation (3.14a ) does not constrain the axial derivative of pressure to be continuous. The continuity and jumps in the acoustic power are correctly captured. Zero flow-rate conditions (hard walls), $\hat{U} =0$ , are imposed on both ends of the device. The corresponding zero-Neumann condition for pressure $\text{d}\hat{p}/\text{d}x=0$ does not need to be explicitly enforced numerically, as it is a natural outcome of the solution of the eigenvalue problem. Inter-segment and boundary conditions are inserted into (3.13), yielding the complete eigenvalue problem. Several analytical results for variable-area duct acoustic systems have been reproduced to machine-precision accuracy (Dowling & Williams Reference Dowling and Williams1983) and excellent agreement with the Navier–Stokes calculations of the TAP engine model is found in both the linear and low-acoustic-amplitude nonlinear regime (as discussed later).

4.1 Engine start-up

Navier–Stokes simulations are carried out first without piezoelectric energy absorption (i.e. with hard-wall boundary conditions on the right end of the resonator) for all cases in table 3. Initial conditions are that of zero velocity, ambient pressure and temperature matching the expected mean axial distribution at equilibrium (figure 1). No initial velocity or pressure perturbations are prescribed. As also observed in Scalo et al. (Reference Scalo, Lele and Hesselink2015b ), for a sufficiently large background temperature gradient, the simple activation of the heat source triggers a disturbance that is thermoacoustically amplified, initiating a transient exponential growth, followed by a saturation of the pressure amplitude (figure 4). During the late stages of energy growth the pressure amplitude overshoots its limit cycle value, especially for close-to-critical values of the temperature gradient. This behaviour was not observed in the travelling-wave engine investigated by Scalo et al. (Reference Scalo, Lele and Hesselink2015b ).

Figure 5. Frequency (a) and growth rate (b) versus temperature difference, $\unicode[STIX]{x0394}T=T_{h}-T_{c}$ , for grid resolution/stack-type C/I. Predictions from linear model (– – –) and Navier–Stokes simulations at start-up (○); limit cycle frequencies (

); square-root fit (4.1) (——) of limit cycle pressure amplitudes (▾) in the hot cavity versus temperature difference; and linear fit of growth rates from Navier–Stokes calculations (——).

Growth rates and frequencies predicted by the linear model developed in § 3 are in good agreement with the nonlinear simulations (figure 5). A linear fit of the growth rates extracted from the Navier–Stokes simulations against the temperature difference, $\unicode[STIX]{x0394}T$ , yields a critical temperature difference of $\unicode[STIX]{x0394}T_{cr}=305~\text{K}$ . Linear theory predicts $\unicode[STIX]{x0394}T_{cr}=315.7~\text{K}$ , while fitting the limit cycle pressure amplitude, $p_{lc}$ , with the equilibrium solution from a supercritical Hopf bifurcation model with dissipation term scaling as $p_{lc}^{2}$ ,

yields $\unicode[STIX]{x0394}T_{cr}=332.7~\text{K}$ . While small discrepancies between linear theory and Navier–Stokes simulations are expected, a difference of 30 K between the $\unicode[STIX]{x0394}T_{cr}$ calculated from the growth rates and the limit cycle pressures suggests that hysteresis effects associated with subcritical bifurcation may be present. This phenomenon was not observed in the numerical simulations of a travelling-wave engine by Scalo et al. (Reference Scalo, Lele and Hesselink2015b ).

Very good agreement is also found between pressure and flow rate eigenfunctions, and pressure and flow-rate amplitudes extracted from the Navier–Stokes calculations via least squares fitting during the start-up phase (figure 6). Results were confirmed with peak finding and windowed short-time Fourier transform (STFT). Minor discrepancies are present at locations of abrupt area change, where assumptions of quasi-one-dimensionality break down. Amplitude and phase detection were applied to time series of cross-sectionally averaged pressure and cross-sectionally integrated axial flow velocity components.

Figure 6. Axial distribution of pressure (a) and flow-rate (b) amplitudes of the thermoacoustically unstable mode for temperature setting 5 and grid resolution/stack-type C/I, as predicted by linear theory (——), rescaled to match amplitudes extracted from Navier–Stokes simulations (○) during the start-up phase. The frequency predicted by linear theory is $f=378.9~\text{Hz}$ , while the frequency extracted from the Navier–Stokes calculations yields $f=381.8~\text{Hz}$ (see figure 5). Vertical dashed lines indicate locations of abrupt area change (figure 1).

While good agreement is also retained in the nonlinear regime, and used to extract the axial distribution of acoustic power at the limit cycle with piezoelectric energy extraction (discussed later, figure 19 in § 7), a frequency shift in the range 5–20 Hz, from high to low temperature settings, is observed in the transient leading to the limit cycle (figures 5, 11). This frequency change, as discussed in § 7, is enough to alter significantly the rate of energy extraction from a piezoelectric diaphragm – indicating that, especially for complex geometries, its fine tuning should be performed based on actual measurements or nonlinear calculations, and not by solely relying on linear theory. The observed frequency shift is also discussed in § 5 in the context of acoustic streaming.

4.2 Grid sensitivity study

Figure 7. Time series of pressure amplitudes within the hot cavity ( $\circ$ ) with semi-logarithmic fit (– – –) over initial start-up phase (a) and growth rates (b) for grid resolution levels A, B and C, temperature setting 5 and stack type I (tables 1, 2 and 3). An estimate of the growth rate (▫) at zero-grid spacing is derived via Richardson extrapolation.

Nonlinear calculations for stack type I and temperature setting 5 have been carried out at all three available grid resolution levels – A, B and C (table 3) – with a successive linear grid refinement factor of approximately $r=\sqrt{N_{cv}^{(C)}/N_{cv}^{(B)}}\simeq \sqrt{N_{cv}^{(B)}/N_{cv}^{(A)}}\simeq \sqrt{2}$ . The order of grid convergence estimated from the growth rates extracted from the Navier–Stokes calculations (figure 7 b) is

where $\unicode[STIX]{x1D6FC}_{A}$ , $\unicode[STIX]{x1D6FC}_{B}$ and $\unicode[STIX]{x1D6FC}_{C}$ are the growth rates associated with the grid resolution levels $A$ , $B$ and $C$ , respectively (table 3). Using Richardson extrapolation, the predicted growth rate in the limit of zero-grid spacing is $\unicode[STIX]{x1D6FC}_{h=0}=75.64~\text{ s}^{-1}$ with an error band of $5.6\,\%$ .

Grid-convergent values of the growth rate show that the amount of acoustic energy lost to the numerical scheme is related to its truncation error. The estimated order of grid convergence (4.2) is, however, only slightly above first order, lower than the nominal order of spatial accuracy of the solver. This result demonstrates the inherent difficulties associated with the extraction of the growth rates from Navier–Stokes simulations of full-scale thermoacoustic devices on unstructured grids. Issues include: the arbitrary choice of the time window used for the semi-logarithmic fit of the pressure amplitude time series (figure 7 a); defining a systematic grid refinement criteria for complex unstructured grids (figure 2) that compensates for changes in the effective order of the numerical discretization scheme due to regions of intense skewness and stretching; and the nature of the growth rate itself, which is an accumulation over several cycles of relatively small amounts of energy per cycle – modelling and/or numerical errors, which would otherwise be deemed negligible, accumulate in the same way.

Despite the aforementioned technical and conceptual issues related to the exact definition of the growth rate, results from Navier–Stokes simulations on the highest grid resolution level (grid C) are in very good agreement with linear theory (figures 5, 6, 8 and 9) and thus will be used for the remainder of the manuscript.

Figure 8. Frequency (a) and growth rate (b) versus stack porosity during start-up phase. Results from the nonlinear Navier–Stokes simulations ( $\circ$ ) for temperature setting 5 and grid resolution/stack-type C/I, C/II and C/III; predictions from linear theory for $n_{s}=3,5,7,9$ and 11 (——). Vertical arrows denote the difference between the linear theory prediction and Navier–Stokes simulations.

4.3 Optimal stack design

The frequency of the thermoacoustically amplified mode decreases with increasing porosity (figure 8 a), i.e. as a larger fraction of the cross-sectional area in the stack is made available to gas flow ( $A_{g}=\unicode[STIX]{x1D719}\,A_{stk}$ ). For example, only by adopting stack I, which matches the porosity of the regenerator adopted in Smoker et al.’s (Reference Smoker, Nouh, Aldraihem and Baz2012) experiments, and by carrying out the simulations to a limit cycle (figure 11), is it possible for the TAP engine model to operate close to the experimentally observed frequency of 388 Hz, at which the piezoelectric diaphragm was tuned.

The growth rate is dramatically affected by the stack porosity (figure 8 b). For any given number of solid annuli $n_{s}$ , decreasing the porosity reduces the volume of gas available to thermoacoustic energy production ( $V_{g}\rightarrow 0$ ), while increasing viscous blockage ( $f_{\unicode[STIX]{x1D708}}\rightarrow 1$ ) leads to negative growth rates in the limit of $\unicode[STIX]{x1D719}\rightarrow 0$ . For example, for $n_{s}=11$ , reducing the porosity from $\unicode[STIX]{x1D719}=0.42$ to $0.2$ reduces the growth rate from $\unicode[STIX]{x1D6FC}_{max}=174.4~\text{s}^{-1}$ to zero. On the other hand, increasing the porosity increases the cross-sectional area available to the gas ( $A_{g}\rightarrow A$ ) at the expense of thermal contact, ultimately leading to an attenuation of the growth rate. For example, for $n_{s}=5$ , a maximum growth rate of $\unicode[STIX]{x1D6FC}_{max}=117.73~\text{s}^{-1}$ is obtained at $\unicode[STIX]{x1D719}=0.308$ ; this declines to $\unicode[STIX]{x1D6FC}=90.24~\text{s}^{-1}$ for $\unicode[STIX]{x1D719}=0.443$ and to $\unicode[STIX]{x1D6FC}=50~\text{ s}^{-1}$ for $\unicode[STIX]{x1D719}=0.6$ . Negative growth rates for $\unicode[STIX]{x1D719}>0.9$ are reached for $n_{s}=3$ . The degree of thermal contact of the $m\text{th}$ annular flow passage is accounted for by the thermoacoustic gain term $(f_{k}^{(m)}-f_{\unicode[STIX]{x1D708}}^{(m)})/(1-f_{\unicode[STIX]{x1D708}}^{(m)})$ in (3.10a ) multiplying the mean temperature gradient, which is the driver of thermoacoustic instability. This term decays to a very small (but non-zero) value for $\unicode[STIX]{x1D719}\rightarrow 1$ (Swift Reference Swift2002, p. 95). The optimal value of porosity, $0<\unicode[STIX]{x1D719}_{opt}<1$ , is therefore the result of a trade-off between thermal contact and available pore volume for thermoacoustic energy production.

Increasing the growth rate for fixed values of porosity is possible by increasing $n_{s}$ . This results in an increased solid-to-gas contact surface $S_{heat}$ and greater thermal contact without increasing flow obstruction. A higher stack density, however, requires a higher porosity to maintain the optimal growth rate, $\unicode[STIX]{x1D6FC}_{max}$ , to compensate for the increased viscous blockage. Moreover, the achievable $\unicode[STIX]{x1D6FC}_{max}$ increases with $n_{s}$ , demonstrating the importance of available surface area $S_{heat}$ : for $n_{s}=3$ , the maximum growth rate achievable is $\unicode[STIX]{x1D6FC}_{max}=89.1~\text{s}^{-1}$ , while for $n_{s}=11$ it increases to $\unicode[STIX]{x1D6FC}_{max}=174.4~\text{s}^{-1}$ .

Figure 9. Axial distribution of pressure (a) and flow rate (b) amplitudes for the second (natural) resonant mode predicted by linear theory (——) rescaled to match pressure and flow-rate amplitudes extracted from Navier–Stokes simulations in figure 10 (○) at 30 cycles. Vertical dashed lines indicate locations of abrupt area change (figure 1).

Higher growth rates lead to higher limit cycle acoustic amplitudes (figure 4); for example, limit cycle pressure amplitudes of approximately $6000~\text{ Pa}$ are obtained for stack type I and temperature setting 5, while stack type II reaches $P_{amp}\simeq 11\,500~\text{ Pa}$ (not shown) for the same imposed temperature gradient. This result reflects the increased thermal contact in stack II, which has almost twice the available solid-to-gas contact surface area of stack I. Stack III exhibits the lowest growth rate due to poor thermal contact, as suggested by the presence of an inviscid core (table 1). Higher growth rates, however, may not straightforwardly be associated with higher thermal-to-acoustic efficiencies; in the case of increased $S_{heat}$ , a higher external thermal energy input will be required.

Figure 10. Time series of pressure in the hot cavity for grid-resolution/stack-type C/I at $\unicode[STIX]{x0394}T=0$ with initial quarter-wavelength pressure distribution of 6000 Pa in amplitude. Time is expressed in cycles of mode 2 with frequency $f=633.9~\text{ Hz}$ , in agreement with linear theory (table 4).

Figure 11. Temporal evolution of frequency of thermoacoustically amplified mode for temperature settings 1–5 and grid resolution/stack-type C/I (table 3). Frequency is obtained via peak finding and windowed over two acoustic periods. Higher temperature differences correspond to lower limit cycle frequencies, as shown by the arrow.

4.4 Unexcited acoustic modes

Deactivating the temperature gradient in the stack allows for the analysis of the natural, unexcited acoustic modes of the TAP engine. Initiating the Navier–Stokes calculations with a large-amplitude quarter-wavelength pressure distribution allows the observation of the simultaneous decay of the first two resonant modes (figure 10). The second mode (633.0 Hz) decays slower than the first mode (335.4 Hz), which is thermoacoustically amplified for $\unicode[STIX]{x0394}T>\unicode[STIX]{x0394}T_{cr}$ . This is due to the structure of the second mode (figure 9), which exhibits relatively low flow-rate amplitudes in the stack, where the most intense viscous losses are concentrated.

The second mode is weakly thermoacoustically sustained by the temperature gradient, as shown by the increase in the growth rate with respect to the unexcited case (table 4). The persistence of a negative growth rate indicates that the associated thermoacoustic energy production (made inefficient by a pressure amplitude minima in the stack), is insufficient to overcome viscous dissipation.

In preliminary numerical trials at piezoelectric energy extraction, mode switching from the first mode to the second mode was mistakenly triggered. This was due to the erroneous application of a physically inadmissible impedance with negative resistance ( $\text{Re}(Z)<0$ ) at frequencies close to 633 Hz. While the second mode is not prone to being thermoacoustically amplified, the erroneously assigned impedance forced the device to operate at a frequency different from the fundamental one, effectively controlling the thermoacoustic response. Admissibility issues arise, in particular, due to the fact that an impedance with negative resistance represents an active boundary element, i.e. it injects acoustic power into the system (Rienstra Reference Rienstra2006).

6.1 One-port electromechanical impedance model

The experimental characterization of the electromechanical frequency response of a PZT-5A (lead zirconate titanate) piezoelectric diaphragm has been carried out by Smoker et al. (Reference Smoker, Nouh, Aldraihem and Baz2012), resulting in the system of equations

where $\hat{x}_{c}$  [m], $\hat{q}$  [C], $\hat{p}$ [Pa] and $\hat{V}$ [V] are, respectively, the complex amplitudes of the fluctuating centreline displacement (positive along the $x$ direction), electric charge, pressure and voltage. The electromechanical admittances, $T_{mn}$ , in (6.1) have been measured for a broadband range of frequencies and fitted with the rational function

where the fitting coefficients $a_{mn}$ and $b_{mn}$ are reported in appendix C. Expressing the centreline displacement, $\hat{x}_{c}$ , and charge, $\hat{q}$ , in terms of velocity, $\hat{u} _{c}$ , and current, $\hat{I}$ , respectively, yields

In the experiments, the piezoelectric diaphragm drives a load of resistance $R_{L}=3170~\unicode[STIX]{x03A9}$ , relating voltage and current via $\hat{V}=R_{L}\hat{I}$ . This allows (6.3) to be collapsed into a one-port model

which corresponds to the (purely mechanical) impedance

and wall softness

where $Z_{0}=\unicode[STIX]{x1D70C}_{0}a_{0}$ is the characteristic specific acoustic impedance of the gas. While the impedance (6.5) is based on experimental measurements of the broadband frequency response of the diaphragm (measured at the centreline only), it is not necessarily computationally stable and/or physically admissible (as discussed below) and therefore cannot be applied directly in time-domain nonlinear Navier–Stokes simulations.

The Navier–Stokes calculations were carried out with a computationally and physically admissible impedance approximating (6.5) (approximations derived below), uniformly applied over a circular area scaled in size to preserve the surface-averaged displacement amplitude of the experimentally measured deflection profile by Smoker et al. (Reference Smoker, Nouh, Aldraihem and Baz2012). This technique allows the matching of overall acoustic power output for the same pressure amplitude levels (figure 13). Impedance boundary conditions impose a specific relationship between the Fourier transforms of velocity and pressure at a stationary boundary and should not be confused with the imposition of a moving boundary. The nature of the resulting power extraction is an acoustic-to-mechanical energy conversion, since it is associated with a mechanical deflection of a membrane driven by acoustic excitation. Mechanical-to-electric energy conversion is not directly accounted for.

Figure 13. Illustration of the PZT-5A piezoelectric diaphragm installation by Smoker et al. (Reference Smoker, Nouh, Aldraihem and Baz2012) capping the resonator of the TAP engine (see figure 1). All lengths are given in millimetres. An aluminium substrate backs the piezoelectric diaphragm with an added weight at the centreline used for tuning. Rational polynomial fit of impulse response measurements of the electromechanical admittances (6.1) has been performed solely based on centreline measurements. To model the piezoelectric diaphragm in the Navier–Stokes simulations, a patch of uniformly distributed impedance is used, with size scaled to match the acoustic power output of the actual PZT-5A diaphragm for the same pressure amplitude levels.

6.2 Single-oscillator approximation

A simple approach towards constructing a computationally admissible impedance approximating the experimental value (6.5) is to use a damped Helmholtz oscillator model (Tam & Auriault Reference Tam and Auriault1996), expressed as the three-parameter impedance

where $R$ , $X_{+1}$ and $X_{-1}$ are the resistance, acoustic mass and stiffness, respectively. Only one undamped resonant frequency,

is associated with (6.7), with corresponding wall softness, expressed in the Laplace domain,

where $s=\text{i}\unicode[STIX]{x1D714}$ , with $\unicode[STIX]{x1D714}$ here being extended to the complex domain (via an abuse of notation). Computational admissibility requires the time-domain equivalent of (6.9) to be causal, that is, the poles of the wall softness must lie in the left half of the $s$ -plane (negative real part) or, equivalently, in the upper half of the complex $\unicode[STIX]{x1D714}$ -plane (positive imaginary part). Poles of $\widehat{\widetilde{W}}_{s}(s)$ in the $s$ -domain are in biunivocal correspondence with the poles of $\widehat{\widetilde{W}}_{\unicode[STIX]{x1D714}}(\unicode[STIX]{x1D714})$ in the complex $\unicode[STIX]{x1D714}$ -domain.

The wall softness of a generic oscillator with a single resonant frequency can be expressed via a decomposition in partial fractions in the Laplace domain,

with one set of complex conjugate residues ( $\unicode[STIX]{x1D707},\unicode[STIX]{x1D707}^{\ast }$ ) and poles ( $p,p^{\ast }$ ), where $\unicode[STIX]{x1D707}=a+\text{i}b$ and $p=c+\text{i}d$ with $a,b,c,d\in \mathbb{R}$ .

In order for ( $\unicode[STIX]{x1D707},\unicode[STIX]{x1D707}^{\ast }$ ) and ( $p,p^{\ast }$ ) to represent a single damped Helmholtz oscillator in the form of the three-parameter model (6.7), the following conditions, derived by comparing (6.10b ) with (6.9), must be satisfied:

where $C=bd+ac$ is the phase parameter. Physical admissibility (boundary is a passive acoustic absorber) and causality require $R=-(1+2\,c/a)>0$ and $c=\text{Re}(p)<0$ , respectively. It is important to stress that a generic oscillator of the form (6.10b ) cannot be equivalent to the single damped Helmholtz oscillator (6.9) unless its phase parameter is zero (6.11a ). Scalo et al. (Reference Scalo, Bodart and Lele2015a ) have demonstrated that it is possible to perform turbulent flow simulations with imposed wall impedance of type (6.9) without encountering numerical stability issues, confirming that (6.9) is in fact physically and computationally admissible.

Fung & Ju (Reference Fung and Ju2001) have suggested that it is not necessary for a single-oscillator model such as (6.10) to have a zero phase parameter for its use in time-domain computations. However, in preliminary numerical trials, it was found that leaving the phase parameter unconstrained ( $C\neq 0$ ) leads to unstable numerical simulations and causing, in our case, spurious mode switching and near-DC (near zero-frequency) acoustic power extraction.

As seen from (6.10b ), the phase parameter is dominant in the low-frequency limit ( $s\rightarrow 0$ ), thus influencing the phase of $Z(\unicode[STIX]{x1D714})$ over a broad range of near-DC frequencies. A non-zero phase parameter yields a purely real, non-zero and finite $Z(\unicode[STIX]{x1D714})$ at zero frequency. Because the experimentally measured wall softness (6.6) has a zero magnitude (infinite impedance magnitude) in the DC limit, a zero phase parameter $C=0$ is necessary in both the single- and multi-oscillator impedance approximations to (6.6) (the latter discussed below) to retain physical admissibility.

Following the aforementioned considerations, the impedance (6.5) was first approximated by the three-parameter impedance model (6.7) (guaranteeing a zero phase parameter) with $R$ , $X_{+1}$ and $X_{-1}$ determined directly via least squares fitting of $\text{Re}(Z_{exp})$ and $\text{Im}(Z_{exp})$ , where $Z_{exp}$ is the impedance corresponding to the collapsed two-port model (6.5). The fitting window used is $f=388~\text{ Hz}\pm 10~\text{ Hz}$ with resulting parameters reported in table 5. As expected, good agreement is found only for frequencies close to $f=388~\text{ Hz}$ (figure 14). The largest discrepancies are in the values of resistance $R$ (not constant in the experiments), which is responsible for differences in the location of the minima of $|Z|$ . The latter is an attractor for the thermoacoustically unstable mode at the limit cycle. Negative values of resistance in the experimentally measured impedance are observed for frequencies above 450 Hz, which is unphysical for a passive acoustic element and therefore are a challenge in the context of deriving a multi-oscillator impedance approximation.

Figure 14. Magnitude (a), phase (b), real part (c) and imaginary part (d) of the experimentally measured impedance (6.5) (○), single-oscillator impedance model (6.7), (6.9) (– – –) and multi-oscillator impedance fit (6.12b ) with $\bar{\unicode[STIX]{x1D6FC}}=0.06$ and $n_{o}=18$ oscillators (——). The single-oscillator model is fitted in the range $388\pm 10~\text{ Hz}$ while the multi-oscillator model is fitted over the entire frequency range and with wall softness $\widehat{\widetilde{W}}$ constrained to have a positive real part, resulting in $\text{Re}(Z)\geqslant -Z_{0}$ (see text). The shaded area highlights the frequency interval $f>450~\text{Hz}$ of negative resistance $\text{Re}\{Z_{exp}(\unicode[STIX]{x1D714})\}<0$ for the experimentally determined impedance (6.5), deemed unphysical.

6.3 Multi-oscillator approximation

In order to fit (6.6) over a broader frequency range, a linear superposition of the wall softness coefficients of $n_{o}$ oscillators, each decomposed in partial fractions with one conjugate pair of residues ( $\unicode[STIX]{x1D707}_{k},\unicode[STIX]{x1D707}_{k}^{\ast }$ ) and poles ( $p_{k},p_{k}^{\ast }$ ), is used, yielding

where (6.12b ) is an alternative form to (6.10b ) adopted by Fung & Ju (Reference Fung and Ju2004), where $\unicode[STIX]{x1D714}_{0,k}$ is the resonant (or basis) frequency (6.8) of the $k\text{th}$ oscillator, $\bar{\unicode[STIX]{x1D6FC}}$ is a damping parameter (common to all oscillators) and $A_{k}$ and $B_{k}$ are fitting coefficients corresponding to $2a$ and $-2\,C$ in the single-oscillator model in (6.10b ). The experimentally measured wall softness is expected to approach zero for $\unicode[STIX]{x1D714}\rightarrow \infty$ and $\unicode[STIX]{x1D714}\rightarrow 0$ (with implications on $B_{k}=0$ , discussed below), making its functional form better suited for fitting than the impedance itself, which, in the case of a damped Helmholtz resonator, diverges for the same extremes. Moreover, fitting the wall softness as a linear superposition of oscillators is consistent with the numerical implementation in the time domain, whereas linearly superimposing impedances is not. Note that the linear superimposition of wall softnesses (as in (6.12a )), which is the approach used in the present work, is not equal to the wall softness resulting from the linear superposition of the corresponding single-oscillator impedances, that is

where

The damping parameter $\bar{\unicode[STIX]{x1D6FC}}$ in (6.12b ) – common to all $n_{o}$ oscillators – controls the bandwidth of the frequency response of each oscillator centred about its basis frequency; for low (high) values of $\bar{\unicode[STIX]{x1D6FC}}$ , each oscillator will exhibit a narrowband (broadband) response. Therefore, for a given fitting frequency window, a low (high) value of $\bar{\unicode[STIX]{x1D6FC}}$ will require a larger (smaller) number of oscillators to approximate a given wall softness. A large number of narrowband oscillators results in a more accurate fit – requiring, however, a closer spacing of basis frequencies.

Figure 15. Magnitude of experimentally measured impedance (6.5) (a) and (b) wall softness magnitude (○) compared with multi-oscillator model for $\bar{\unicode[STIX]{x1D6FC}}=0.2237$ ( $n_{o}=1$ ) (—  —), $\bar{\unicode[STIX]{x1D6FC}}=0.12$ ( $n_{o}=4$ ) (– –), $\bar{\unicode[STIX]{x1D6FC}}=0.10$ ( $n_{o}=6$ ) (- - -), $\bar{\unicode[STIX]{x1D6FC}}=0.08$ ( $n_{o}=13$ ) (

), to $0.06$ ( $n_{o}=18$ ) (——). The shaded area highlights the frequency interval $f>450~\text{Hz}$ of negative resistance $\text{Re}\{Z_{exp}(\unicode[STIX]{x1D714})\}<0$ for the experimentally determined impedance (6.5), deemed unphysical.

The impedance $Z_{exp}$ has been fitted with the following numbers of oscillators: $n_{o}=1$ , $n_{o}=4$ , $n_{o}=6$ , $n_{o}=13$ , $n_{o}=18$ (see (6.12b ) and figure 15). For the single-oscillator case, values of $f_{0,1}(=\unicode[STIX]{x1D714}_{0,1}/2\unicode[STIX]{x03C0})$ , $\bar{\unicode[STIX]{x1D6FC}}$ , $A_{1}$ , corresponding to the single-oscillator model in (6.9), are reported in the third row of table 5. For the multi-oscillator case, $n_{o}>1$ , values of $f_{0,k}(=\unicode[STIX]{x1D714}_{0,k}/2\unicode[STIX]{x03C0})$ , $\bar{\unicode[STIX]{x1D6FC}}$ , $A_{k}$ are reported in table 6. For a given $\bar{\unicode[STIX]{x1D6FC}}$ , basis frequencies were selected through a gradient descent-based iterative method such that the approximate impedance $Z_{fit}$ is the least squares minimizer of $\log \left|Z_{exp}\right|-\log |Z_{fit}|$ and $\arg Z_{exp}-\arg Z_{fit}$ .

As $n_{o}$ increases, the basis frequencies are more closely spaced, corresponding to a decrease in $\bar{\unicode[STIX]{x1D6FC}}$ and an increase in the accuracy of the fit in the frequency domain (figure 15). For each $\bar{\unicode[STIX]{x1D6FC}}$ and thus for a particular $n_{o}$ , the impedance, as a function of frequency and basis frequencies, is fitted with least squares over frequencies $f\in [1,440]~\text{ Hz}$ , with $5.5$ -fold weighting on $f\in [360,440]~\text{ Hz}$ and $22$ -fold weighting on $f\in [378,398]~\text{ Hz}$ .

As seen in figure 14, at higher frequencies, the real component of the experimentally measured impedance becomes negative, which is not consistent with a passive acoustic element and may be the spurious result of the sampling rate used for the eigensystem realization algorithm as reported by Smoker et al. (Reference Smoker, Nouh, Aldraihem and Baz2012) or simply the extrapolation of the rational polynomial fit beyond the tuned frequency of the piezoelectric diaphragm. To avoid unphysical values of the reconstructed impedance at high frequencies, $A_{k}$ in (6.12b ) is constrained to be positive, since negative $A_{k}$ can lead to unbounded negative resistance. By combining the constraints

where $B_{k}$ is the phase parameter (see discussion in § 6.2), with (6.6) and (6.12b ), the real part of the resulting impedance $Z$ has a lower bound, i.e. $\text{Re}(Z)>-Z_{0}$ . Without such a constraint, a multi-oscillator impedance model could fit, with arbitrary accuracy, the given experimentally measured impedance (6.5), but would cause the impedance to inject energy into the system for $f>450~\text{Hz}$ , hence exciting its second mode (at $\simeq 640~\text{Hz}$ , where $\text{Re}(Z_{exp})$ is large and negative). This causes unphysical mode switching (see § 4.4), with the piezoelectric diaphragm no longer acting as a passive element but a (spurious) driver of oscillations.

In the following, results from Navier–Stokes calculations with piezoelectric energy absorption with the multi-oscillator model with $n_{o}=18$ and $\bar{\unicode[STIX]{x1D6FC}}=0.06$ are shown, since this model provides the highest level of fidelity over the frequency range of interest.

7.1 Thermal-to-mechanical efficiency

Acoustic energy extraction, modelled via the TDIBCs designed in § 6, is applied once a limit cycle without energy absorption is achieved, and only for the TAP engine model with stack type I. The latter most closely matches the porosity and hydraulic radius of the regenerator used in the experiments (table 1) and, as a result, the operating frequency at which the piezoelectric diaphragm is tuned ( ${\sim}388~\text{Hz}$ ).

Figure 16. Time series of pressure in the hot cavity for temperature setting 5 and grid resolution/stack-type C/I from start-up to limit cycle, without ( $t<0.544~\text{ s}$ ) and with ( $t>0.544~\text{ s}$ ) piezoelectric energy absorption, as modelled by the multi-oscillator impedance model (6.12b ) with $n_{o}=18$ (table 6).

The imposition of the impedance boundary conditions designed in § 6 results in a decrease in the pressure amplitude (figure 16), corresponding to an extraction of acoustic energy (figure 17), following an initial assessment phase with spurious high-frequency oscillations due to the abrupt initialization of the convolution integral (B 5). A new limit cycle is rapidly obtained with a slight frequency shift due to the resonance tuning of the piezoelectric diaphragm.

For temperature setting 5, acoustic energy absorption results in a pressure amplitude decrease of 10 %. The same acoustic energy absorption with temperature setting 1 (the close-to-critical temperature gradient) suppresses the thermoacoustic instability. The net power output per cycle (figure 17),

is extracted via sharp spectral filtering of the instantaneous acoustic power output,

where $p^{\prime }$ and $U^{\prime }$ are the pressure and surface-averaged volumetric flow-rate amplitudes at the diaphragm location. The convolution integral in (7.1) is, in practice, limited to $\pm 17$ acoustic cycles with a cutoff frequency of $f_{c}=22.5~\text{ Hz}$ , lower than half that of the thermoacoustically amplified mode. The power extracted at the boundary is, at most, 111.25 mW, corresponding to temperature setting 5. Thermal-to-mechanical efficiency $\unicode[STIX]{x1D702}$ is calculated for each case as the ratio of $\overline{P}_{out}$ and the cycle-averaged heat transfer rate through the stack walls, and is at most 1.3 % (table 7).

Figure 17. Time series of instantaneous acoustic power (7.2) (——) extracted at the limit cycle (left axis), cycle-averaged power (7.1) shown with the shaded area (right axis) for temperature setting 5, grid resolution/stack-type C/I. The beginning of the time series in this figure corresponds to the vertical dashed line in figure 16.

Figure 18. Limit cycle frequency, $\unicode[STIX]{x1D714}/2\unicode[STIX]{x03C0}$ , versus temperature difference (a) with active piezoelectric energy extraction; cycle-averaged power output (7.1) (○) and the square of limit cycle pressure amplitudes in the hot cavity $p_{lc,cav}^{2}$ (▾) versus temperature difference (b). Results for grid resolution/stack-type C/I. Note that limit cycle frequency values reported here differ from the ones in figure 5(a), as the latter are obtained without imposition of piezoelectric energy extraction.

In the experiments, an acoustic power output of $1.32~\text{mW}$ is reported for conditions nominally meant to match the temperature setting 5 used in the present TAP engine model (Smoker et al. Reference Smoker, Nouh, Aldraihem and Baz2012). However, results in Nouh et al. (Reference Nouh, Aldraihem and Baz2014) from the same engine show that thermoacoustic heat leakage and natural relaxation of the thermal gradient in the stack leads to unsteady temperature distributions in the stack, approaching temperature setting 1. Due to differences in the regenerator/stack and uncertainties in the actual temperature gradient used in the experiments, numerical simulations with (steady) isothermal conditions cannot reproduce the experimentally observed limit cycle acoustic pressure amplitude. A normalized power output can be defined by compensating for the differences in pressure amplitude,

where $A_{piezo}$ is the area of the equivalent uniform impedance patch used in the present simulations and the superscript $(\text{exp})$ indicates experimental values. The good matching observed between the two non-dimensional powers (table 7) confirms that the impedance boundary conditions are imposing the correct phasing between pressure and velocity.

After the application of the TDIBC, the limit cycle operating frequency shifts (not shown) towards the frequency corresponding to the minimum impedance magnitude (maximum acoustic energy absorption). This is due to the increased compliance of the piezoelectric diaphragm at higher frequencies, corresponding to a reduction in the value of the resistance at higher frequencies (as seen in figure 14 and discussed in § 6.3). In the case of the single-oscillator impedance model (6.7) with a constant value of resistance, the limit cycle frequency is controlled exclusively by the reactance. In all cases, an excessively large shift in frequency would disrupt the thermoacoustic phasing in the stack, leading to a suppression of the instability.

Figure 19. Pressure (a) and flow-rate (b) amplitudes of the thermoacoustically amplified mode predicted by linear theory (——) rescaled to match amplitudes extracted from companion Navier–Stokes simulations (○) for temperature setting 5, grid resolution/stack-type C/I, with active energy extraction at the limit cycle. Minor losses have been incorporated; however, the exclusion of minor losses (not shown) does not significantly alter the amplitudes predicted by linear theory. (——). (c) Axial distribution of acoustic power (7.4) from eigenfunctions predicted by linear theory (§ 3), without minor losses (– – –), with linearized minor losses ((7.6) and (7.7)) (

), and with minor losses calibrated from the Navier–Stokes simulations (——). Also shown are the acoustic power values extracted from the simulations, using only centreline (○) and using full cross-sectionally averaged (●) values of axial velocity and pressure. Inter-segment locations are numbered above, referenced in table 8(b). Results correspond to temperature setting 5 and grid resolution/stack-type C/I, with active energy extraction. Vertical dashed lines indicate locations of abrupt area change (figure 1).

The linear thermoacoustic model developed in § 3 has been augmented with minor losses (§ 7.2) and is used here to reconstruct the axial distribution of acoustic power by applying a (constant) impedance, $Z(\unicode[STIX]{x1D714}_{0})$ where $\unicode[STIX]{x1D714}_{0}/2\unicode[STIX]{x03C0}=388.0~\text{ Hz}$ at $x=0.51~\text{m}$ . The axial distribution of acoustic power ${\dot{W}}$ can then be calculated as

where $\hat{p}$ and $\hat{U}$ are the eigenvectors predicted by the linear model and rescaled such that pressure and volume flow-rate amplitudes match the Navier–Stokes calculations with TDIBCs, at limit cycle. The resulting eigenfunctions and axial power distribution are shown in figure 19. Results from the Navier–Stokes calculations were collapsed into axial time series of volume flow rate and pressure via surface integration, and the resulting quantities of $\hat{U} (x)$ and $\hat{p}(x)$ were fitted to complex phasors. Acoustic power along the axis is then calculated using (7.4). The linear model predicts an acoustic power extraction of $111.09~\text{mW}$ , in good agreement with the result of $111.25~\text{ mW}$ from the Navier–Stokes calculations, as reported in table 7.

7.2 Acoustic energy budgets

As expected, a positive slope in the acoustic power is present in the stack, while a negative slope in the pulse tube and the resonator volume indicates acoustic power dissipation due to viscous dissipation and thermal relaxation. The acoustic power distribution is consistent with that predicted by engineering design software such as DeltaEC in other standing-wave engines in literature (Swift Reference Swift1992; Ward et al. Reference Ward, Clark and Swift2012).

The balance of acoustic energy at the limit cycle can be heuristically expressed as

where the cycle-averaged acoustic energy per unit length, $E_{a}$ , is assumed to be steady. The divergence of the acoustic energy flux is balanced by thermoacoustic source terms, $S_{ta}$ , and viscous dissipation, $D_{\unicode[STIX]{x1D707}}$ , which are accurately predicted by Rott’s theory, as made evident by the predicted slope matching that of the data extracted from nonlinear calculations. Hydrodynamic minor losses due to abrupt area changes, $D_{m}$ , manifest themselves as jumps in the value of ${\dot{W}}$ and will therefore be formally incorporated in the budgets above via Dirac functions (table 8 a,b).

Incorporating minor losses is necessary to reproduce first-order features in the limit cycle acoustic power distribution for the thermoacoustic piezoelectric engine. In the case of steady flow, minor losses due to abrupt area changes can be parameterized via the Borda–Carnot formula,

in the case of expansions or via the formula presented by Idelchik (Reference Idelchik2003),

in the case of sudden contractions, where $A_{0}$ and $A_{1}$ are the smaller and larger areas. In the present study, the two losses are combined, $\unicode[STIX]{x1D701}=\unicode[STIX]{x1D701}_{e}+\unicode[STIX]{x1D701}_{c}$ , and the pressure drop condition

linearized about the limit cycle axial velocity amplitude distribution $u_{lc}$ , is incorporated in the inter-segment condition (3.14a ) in the linear thermoacoustic model in § 3. While a similar approach to modelling minor losses in oscillatory flows has been adopted with great success by Ward et al. (Reference Ward, Clark and Swift2012), a more accurate parameterization of $\unicode[STIX]{x1D701}$ should be derived with ad hoc numerical or experimental investigations. Despite the strong assumptions made in deriving and introducing minor losses in the linear model, the agreement with the Navier–Stokes data is remarkable (figure 19 c). Axial amplitude profiles of pressure and flow rate are, however, not visibly altered by the addition of minor losses; the condition (7.8) primarily affects the pressure–velocity phasing in locations of sudden area change.

It is possible, however, to determine the value of $u_{lc}$ in (7.8) a priori by imposing a zero growth-rate condition in the eigenvalue solver in § 3. This is indirectly done in modelling software such as DeltaEC, where the limit cycle pressure and velocity amplitudes are found iteratively by assuming that acoustic energy budgets are balanced at steady-state conditions.