4.1. A Lagrangian model of the thermodynamic cycle

Inspired by Swift (Reference Swift1988)’s theoretical TAE (figure 6, top), a one-dimensional analytical model is derived in the following to rigorously explain the interaction between a mean temperature gradient and a Lagrangian fluid parcel in a generic planar acoustic wave field. The latter can be expressed as the linear superposition of a left- and a right-travelling wave, which in complex form reads

(4.1) $$\begin{eqnarray}\displaystyle & \hat{p}(x)=p^{-}\text{e}^{\text{i}[kx+{\it\phi}^{-}]}+p^{+}\text{e}^{\text{i}[-kx-{\it\phi}^{+}]}, & \displaystyle\end{eqnarray}$$

(4.2) $$\begin{eqnarray}\displaystyle & {\it\rho}_{0}a_{0}\hat{u} (x)=-p^{-}\text{e}^{\text{i}[kx+{\it\phi}^{-}]}+p^{+}\text{e}^{\text{i}[-kx-{\it\phi}^{+}]}, & \displaystyle\end{eqnarray}$$

$p^{+/-}$

${\it\phi}^{+/-}$

$k$

$p^{\prime }=\hat{p}(x)\text{e}^{\text{i}{\it\omega}t}$

$u^{\prime }=\hat{u} (x)\text{e}^{\text{i}{\it\omega}t}$

$\text{i}$

${\it\rho}_{0}a_{0}$

$\{{\it\rho}_{0},T_{0},P_{0}\}$

Let $x_{p}=\overline{x}+x_{p}^{\prime }$ be the instantaneous position of a fluid parcel oscillating with small displacements $x_{p}^{\prime }$ about the position $\overline{x}$ , where a linear temperature gradient is locally imposed (figure 6). The fluid parcel velocity can be approximated, based on the assumption $k\max (x_{p}^{\prime })\ll 1$ , as

(4.3) $$\begin{eqnarray}\frac{\text{d}}{\text{d}t}x_{p}^{\prime }=u^{\prime }(\overline{x}+x_{p}^{\prime },t)=\text{Re}\{\hat{u} (\overline{x}+x_{p}^{\prime })\text{e}^{\text{i}{\it\omega}t}\}\simeq \text{Re}\{\hat{u} (\overline{x})\text{e}^{\text{i}{\it\omega}t}\},\end{eqnarray}$$

yielding the relation

(4.4) $$\begin{eqnarray}\hat{x}_{p}=\hat{u} (\overline{x})/\text{i}{\it\omega}.\end{eqnarray}$$

Introducing a Lagrangian base state ${\it\rho}_{L,0},P_{L,0},T_{L,0}$ (specified later) for the fluid parcel density, entropy and temperature,

(4.5) $$\begin{eqnarray}\displaystyle & {\it\rho}_{L}={\it\rho}_{L,0}+\text{Re}\{\hat{{\it\rho}}_{L}\text{e}^{\text{i}{\it\omega}t}\}, & \displaystyle\end{eqnarray}$$

(4.6) $$\begin{eqnarray}\displaystyle & s_{L}=s_{L,0}+\text{Re}\{{\hat{s}}_{L}\text{e}^{\text{i}{\it\omega}t}\}, & \displaystyle\end{eqnarray}$$

(4.7) $$\begin{eqnarray}\displaystyle & T_{L}=T_{L,0}+\text{Re}\{\hat{T}_{L}\text{e}^{\text{i}{\it\omega}t}\}, & \displaystyle\end{eqnarray}$$

$x_{p}$

(4.8) $$\begin{eqnarray}{\it\rho}_{L}T_{L}\frac{\text{d}s_{L}}{\text{d}t}=-{\it\alpha}_{T}[T_{L}-T_{0}(x_{p})],\end{eqnarray}$$

where the (total) time derivative on the left-hand side, applied to the Lagrangian fluid parcel’s specific entropy, has replaced the material derivative of the Eulerian entropy field.

Substituting Gibbs’ relationship linearized about the Lagrangian base state $T_{L,0},p_{L,0}$ ( $=p_{0}$ ),

(4.9) $$\begin{eqnarray}\text{d}s_{L}=\frac{c_{p}}{T_{L,0}}\text{d}T_{L}-\frac{R}{p_{0}}\text{d}p_{L},\end{eqnarray}$$

into (4.8), which is also linearized assuming ${\it\rho}_{L}T_{L}\simeq {\it\rho}_{L,0}T_{L,0}$ , yields

(4.10) $$\begin{eqnarray}{\it\rho}_{L,0}c_{p}\frac{\text{d}T_{L}}{\text{d}t}-\frac{\text{d}p_{L}}{\text{d}t}=-{\it\alpha}_{T}[T_{L}-T_{0}(x_{p})].\end{eqnarray}$$

The imposed mean temperature at the particle location, $T_{0}(x_{p})$ , can be expressed via the Taylor expansion

(4.11) $$\begin{eqnarray}T_{0}(x_{p})=T_{0}(\bar{x})+\left.x_{p}^{\prime }\frac{\text{d}T_{0}}{\text{d}x}\right|_{x=\bar{x}}\end{eqnarray}$$

which is exact in the case of a linear base temperature distribution. Substituting (4.11) into (4.10) yields

(4.12) $$\begin{eqnarray}{\it\rho}_{L,0}c_{p}\frac{\text{d}T_{L}}{\text{d}t}-\frac{\text{d}p_{L}}{\text{d}t}=-{\it\alpha}_{T}\left[T_{L}-T_{0}(\bar{x})-\left.x_{p}^{\prime }\frac{\text{d}T_{0}}{\text{d}x}\right|_{x=\bar{x}}\right].\end{eqnarray}$$

Defining the Lagrangian-based state for temperature such that $T_{L,0}=T_{0}(\overline{x})$ , assuming $p_{L}=p$ and switching to the complex form,

(4.13) $$\begin{eqnarray}{\it\rho}_{L,0}c_{p}\text{i}{\it\omega}\hat{T}_{L}-\text{i}{\it\omega}\hat{p}=-{\it\alpha}_{T}\left[\hat{T}_{L}-\left.\hat{x}_{p}\frac{\text{d}T_{0}}{\text{d}x}\right|_{x=\bar{x}}\right],\end{eqnarray}$$

where the Lagrangian base density is ${\it\rho}_{L,0}=p_{0}/(RT_{L,0})$ .

If the acoustic field is assigned, so are the complex pressure amplitude $\hat{p}$ and particle displacement $\hat{x}_{p}=(\text{i}{\it\omega})^{-1}\hat{u} (\overline{x})$ , which then allows one to solve for $\hat{T}_{L}$ from (4.13). The density of the Lagrangian parcel can then be calculated from the linearized equation of state

(4.14) $$\begin{eqnarray}\hat{{\it\rho}}_{L}=\frac{\hat{p}-{\it\rho}_{L,0}R\hat{T}_{L}}{RT_{L,0}}.\end{eqnarray}$$

The work done by the fluid parcel on the surrounding ambient per unit time (generated acoustic power) is

(4.15) $$\begin{eqnarray}{\dot{w}}=-\overline{\frac{p}{{\it\rho}}\frac{\text{D}{\it\rho}}{\text{D}t}}\simeq -\frac{1}{{\it\rho}_{L,0}}\frac{1}{2}\text{Re}\{\hat{p}(\text{i}{\it\omega}\hat{{\it\rho}}_{L})^{\ast }\}.\end{eqnarray}$$

Depending on the nature of acoustic wave being imposed (figure 6 a), ranging from purely left-travelling to purely right-travelling, a different phasing between $\hat{p}$ and $\hat{u}$ is achieved, determining ${\dot{w}}$ (figure 6 b). In the case of a standing wave, acoustic energy will be absorbed ( ${\dot{w}}<0$ ) or generated ( ${\dot{w}}>0$ ) depending on the location of the temperature gradient, $\overline{x}$ . For a sufficiently high amplitude of the right-, $p^{+}\gg p^{-}$ , (or left-, $p^{-}\gg p^{+}$ ) travelling wave, acoustic power is ultimately only absorbed (or generated) for any $\overline{x}$ in the case of a negative mean temperature gradient, $\text{d}T_{0}/\text{d}x<0$ (figure 6 b). This shows that an acoustic wave travelling in the same direction as or opposite to the imposed mean temperature gradient (applied over a region small compared to the wavelength) will be amplified or absorbed, respectively. Moreover, the acoustic power associated with the energy conversion occurring in an (almost) purely travelling wave is remarkably higher (see area enclosed by the $p{-}v$ diagrams in figure 6 c) than that of a standing wave of comparable amplitude. This confirms Ceperley (Reference Ceperley1979)’s seminal intuition that led to the revolutionary concept of travelling-wave energy conversion, trumping thereafter designs based on standing waves.

4.2. Acoustic network of travelling waves

In spite of the finite amplitude of the initial perturbation (exceeding 1 % of the base pressure) and the immediate establishment of nonlinear effects, the exponentially growing acoustic amplitude (with uniform growth rate in the entire system) suggests that the system-wide behaviour in the startup phase can be analysed by invoking linear acoustics. The low frequencies observed in the numerical simulations and the high aspect ratio of the resonator (with lowest cut-on frequency ${\sim}$ 1.7 kHz) allows us to neglect radial or azimuthal acoustic modes at the resonator scale. This suggests that, as a first approximation, our analysis can be restricted to planar waves. The nature of the thermoacoustic instability in the REG/HX (analysed in § 4.1) suggests that the pulse tube acts as an amplifier of left-travelling waves. These are expected to propagate into the resonator, be reflected back and, upon returning to the REG/HX unit as right-travelling waves, propagate both through the pulse tube and in the feedback inertance (figure 1), being absorbed in the former case (figure 6 c) and propagating freely in the latter. The acoustic power propagating through the inertance is looped back into the REG/HX unit via the compliance, hence creating a network of self-amplified travelling waves. This scenario, typical of travelling-wave TAEs built in a looped configuration (Backhaus & Swift Reference Backhaus and Swift2000), is confirmed in the following analysis.

An exact local decomposition in terms of left $(-)$ and right $(+)$ travelling waves

(4.16a ) $$\begin{eqnarray}\displaystyle & p^{\prime }(t)=p^{-}f({\it\omega}t+{\it\phi}^{-})+p^{+}f(-{\it\omega}t+{\it\phi}^{+}), & \displaystyle\end{eqnarray}$$

(4.16b ) $$\begin{eqnarray}\displaystyle & u^{\prime }(t)=-u^{-}f({\it\omega}t+{\it\phi}^{-})+u^{+}f(-{\it\omega}t+{\it\phi}^{+}), & \displaystyle\end{eqnarray}$$

$p^{\pm }$

$u^{\pm }$

${\it\phi}^{\pm }$

$f(\,)$

$2{\rm\pi}$

${\it\omega}$

${\it\rho}_{0}a_{0}$

Figure 7. Left (–, ◃, ◂) and right ( $+$ , ▹, ▸) travelling-wave amplitudes, $p^{\pm }$ , (a) and phases, ${\it\phi}^{\pm }$ , (b) extracted from the Navier–Stokes calculations based on the linear approximation (4.16) for time $t=0.55~\text{s}$ for the $T_{h}=500~\text{K}$ case on grid A. The locations where data is extracted are shown in the top subfigure. White and black symbols correspond to values extracted on the centreline and in the feedback inertance, respectively.

The acoustic perturbation (4.16) can be rewritten in the form of a complex Fourier series

(4.17a ) $$\begin{eqnarray}\displaystyle & \displaystyle \mathop{\sum }_{k=-\infty }^{+\infty }\hat{p}_{k}\text{e}^{\text{i}k{\it\omega}t}=p^{-}\mathop{\sum }_{k=-\infty }^{+\infty }\hat{f}_{k}\text{e}^{\text{i}k[{\it\omega}t+{\it\phi}^{-}]}+p^{+}\mathop{\sum }_{k=-\infty }^{+\infty }\hat{f}_{k}^{\ast }\text{e}^{\text{i}k[{\it\omega}t-{\it\phi}^{+}]}, & \displaystyle\end{eqnarray}$$

(4.17b ) $$\begin{eqnarray}\displaystyle & \displaystyle {\it\rho}_{o}a_{o}\mathop{\sum }_{k=-\infty }^{+\infty }\hat{u} _{k}\text{e}^{\text{i}k{\it\omega}t}=-p^{-}\mathop{\sum }_{k=-\infty }^{+\infty }\hat{f}_{k}\text{e}^{\text{i}k[{\it\omega}t+{\it\phi}^{-}]}+p^{+}\mathop{\sum }_{k=-\infty }^{+\infty }\hat{f}_{k}^{\ast }\text{e}^{\text{i}k[{\it\omega}t-{\it\phi}^{+}]}, & \displaystyle\end{eqnarray}$$

$p^{\pm }=\pm {\it\rho}_{0}a_{0}u^{\pm }$

$^{\ast }$

$m$

$p^{\pm }$

${\it\phi}^{\pm }$

(4.18a ) $$\begin{eqnarray}\displaystyle & \displaystyle \hat{p}_{m}=\hat{f}_{m}[p^{-}\text{e}^{\text{i}m{\it\phi}^{-}}]+\hat{f}_{m}^{\ast }[p^{+}\text{e}^{-\text{i}m{\it\phi}^{+}}], & \displaystyle\end{eqnarray}$$

(4.18b ) $$\begin{eqnarray}\displaystyle & \displaystyle {\it\rho}_{o}a_{o}\hat{u} _{m}=-\hat{f}_{m}[p^{-}\text{e}^{\text{i}m{\it\phi}^{-}}]+\hat{f}_{m}^{\ast }[p^{+}\text{e}^{-\text{i}m{\it\phi}^{+}}]. & \displaystyle\end{eqnarray}$$

(4.19a ) $$\begin{eqnarray}\displaystyle & \displaystyle p^{-}=|\hat{p}_{m}-{\it\rho}_{0}a_{0}\hat{u} _{m}|/|2\hat{f}_{m}|, & \displaystyle\end{eqnarray}$$

(4.19b ) $$\begin{eqnarray}\displaystyle & \displaystyle p^{+}=|{\it\rho}_{0}a_{0}\hat{u} _{m}+\hat{p}_{m}|/|2\hat{f}_{m}^{\ast }|. & \displaystyle\end{eqnarray}$$

${\it\phi}^{-}$

${\it\phi}^{+}$

$p^{-}$

$p^{+}$

This procedure is applied to the discrete set of points in figure 7(top) located along the resonator axis and around the pulse tube. Results show that left-travelling waves leaving the REG/HX unit propagate into the resonator and are reflected back with a slightly lower amplitude due to losses in the resonator (figure 7 a). Consistently with the one-dimensional approximation in (4.16), the acoustic power can be expressed as

(4.20) $$\begin{eqnarray}{\dot{E}}_{a}^{+}=\int _{A}\overline{p^{\prime }u^{\prime }}\text{d}A=\frac{A}{{\it\rho}_{0}a_{0}}[{p^{+}}^{2}-{p^{-}}^{2}]\mathop{\sum }_{k=-\infty }^{+\infty }{\hat{f}_{k}}^{\ast }\hat{f}_{k},\end{eqnarray}$$

providing an energetic interpretation to the imbalance $p^{-}\gtrless p^{+}$ . The data (figure 7 a,b) confirm that, as anticipated earlier in this section, the acoustic power generated in the REG/HX is fed back to it by the compliance after being channelled through the inertance, where ${p^{+}}^{2}\gg {p^{-}}^{2}$ .

The adopted coaxial configuration is not efficient as it allows part of the acoustic power that is returned by the resonator to be absorbed by the pulse tube. This results in a spatial distribution of the phases ${\it\phi}^{-}$ and ${\it\phi}^{+}$ almost resembling a standing wave (figure 7 b). However, the clear deviation in the REG/HX and in the inertance from such a pattern, both in the phases, ${\it\phi}^{\pm }$ , and in the wave amplitudes, $p^{\pm }$ , is due to the acoustic energy production mechanisms illustrated in figure 5 c and modelled in § 4.1, which are based on the travelling-wave concept (Ceperley Reference Ceperley1979).

The presence of acoustic power being looped around the REG/HX unit (the feedback inertance), which is a component not found in standing-wave TAEs, effectively qualifies the proposed model as a travelling-wave TAE and, as discussed later in § 5, is directly responsible for nonlinear processes such as Gedeon streaming. If the same model had been rearranged in a looped configuration, the departure from a pure standing wave in the REG/HX unit and inertance would only be (quantitatively) more significant, while retaining the same (qualitative) structure as found in the present configuration.

4.3. System-wide linear modelling

Results shown so far suggest that nonlinearities do not play an important role in explaining the acoustic energy propagation and amplification mechanisms during the startup phase. However, given the high amplitude of the initial perturbation ( ${\sim}1$  kPa) and the presence of complex geometrical features such as the sharp edges of the pulse tube (inducing vortex shedding from the start), it is important to assess to what extent a system-wide linear model is able to quantitatively explain the observed instability.

Building upon well-established linear modelling approaches (Rott Reference Rott1969; Ward & Swift Reference Ward and Swift1994; de Waele Reference de Waele2009), the engine is divided into a collection of control volumes (figure 8), representing different components, each modelled as a one- or zero-dimensional element, exchanging acoustic power and mass with adjacent components. Data from the Navier–Stokes simulations suggests that the pressure field is uniform within the compliance, consistent with de Waele (Reference de Waele2009)’s modelling choices. Such a component is, in fact, compact with respect to the fundamental wavelength. By imposing the conservation of mass and assuming an isentropic relation between volume-averaged density and pressure variations one obtains, in the time domain,

(4.21) $$\begin{eqnarray}\frac{\text{d}}{\text{d}t}p_{c}^{\prime }=w_{c_{0}}[U_{i_{1}}^{\prime }-U_{t_{0}}^{\prime }],\end{eqnarray}$$

where $p_{c}^{\prime }$ , $U_{i_{1}}^{\prime }$ and $U_{t_{0}}^{\prime }$ are, respectively, the instantaneous fluctuating pressure in the compliance and flow rates exchanged with the inertance through surface $i_{1}$ , and with the pulse tube through $t_{0}$ (figure 8), and $w_{c_{0}}={\it\gamma}P_{0}/V_{c_{0}}$ , where $P_{0}$ is the base pressure and $V_{c_{0}}$ the volume of the compliance.

Figure 8. Sketch illustrating the subdivision of the engine into control volumes representing different components (cf. 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 segments (oriented according to arrows) modelling the pulse tube, inertance and resonator. Complex pressure, volumetric flow rate and temperature are indicated with $\hat{p}$ , $\hat{U}$ and $\hat{T}$ respectively. References to the TBT are dropped in the context of the analysis of the startup phase during which nonlinear effects can be neglected. The left end of the control volume representing the junction is located at $x=-23~\text{mm}$ .

The very small variation among the growth rates and frequencies extracted from the numerical simulations at the locations in figure 7 (not shown) suggests that normal modes can be assumed for all fluctuating quantities. By adopting the $\text{e}^{+\text{i}{\it\sigma}t}$ convention with ${\it\sigma}=-\text{i}{\it\alpha}+{\it\omega}$ , where ${\it\alpha}$ and ${\it\omega}$ are, respectively, the growth rate and the angular frequency, (4.21) becomes

(4.22) $$\begin{eqnarray}\text{i}{\it\sigma}\hat{p}_{c}=w_{c_{0}}[\hat{U} _{i_{1}}-\hat{U} _{t_{0}}].\end{eqnarray}$$

The same modelling approach is adopted for convenience at the junction, where the conservation of mass reads

(4.23) $$\begin{eqnarray}\text{i}{\it\sigma}\hat{p}_{J}=w_{J_{0}}[\hat{U} _{R_{1}}-\hat{U} _{i_{0}}+\hat{U} _{t_{1}}]\end{eqnarray}$$

with $w_{J_{0}}={\it\gamma}P_{0}/V_{J_{0}}$ , where $V_{J_{0}}$ is the volume of the control volume modelling the junction (figure 8).

Phase variations along the axial coordinate, $x$ , are significant for the other components of the engine and the direct application of the complete set of linearized Euler equations is necessary. In all cases, the fluctuating field and the base state, defined by ${\it\rho}_{0}$ , $T_{0}$ and $P_{0}$ , are assumed to be exclusively a function of $x$ .

For the resonator (R) and feedback inertance ( $\text{i}$ ) isentropic wave propagation is assumed, yielding a simplified set of linearized equations for mass and momentum,

(4.24a ) $$\begin{eqnarray}\displaystyle & \displaystyle \text{i}{\it\sigma}\hat{p}=-\frac{{\it\rho}_{0}a_{0}^{2}}{A(x)}\frac{\text{d}\hat{U} }{\text{d}x}, & \displaystyle\end{eqnarray}$$

(4.24b ) $$\begin{eqnarray}\displaystyle & \displaystyle \text{i}{\it\sigma}\hat{U} =-\frac{A(x)}{{\it\rho}_{0}}\frac{\text{d}\hat{p}}{\text{d}x}, & \displaystyle\end{eqnarray}$$

$A(x)$

$a_{0}=\sqrt{{\it\gamma}RT_{0}}$

(4.25) $$\begin{eqnarray}(\text{i}{\it\sigma}\unicode[STIX]{x1D644}-\boldsymbol{B}_{R})\boldsymbol{\cdot }\boldsymbol{u}_{R}=0\end{eqnarray}$$

for the resonator, and

(4.26) $$\begin{eqnarray}(\text{i}{\it\sigma}\unicode[STIX]{x1D644}-\boldsymbol{B}_{i})\boldsymbol{\cdot }\boldsymbol{u}_{i}=0\end{eqnarray}$$

for the inertance, where $\unicode[STIX]{x1D644}$ is the identity matrix, $\boldsymbol{B}$ is an operator discretizing the right-hand side of (4.24), and $\boldsymbol{u}$ is a discrete collection of complex amplitudes for pressure and flow rate, specifically, $\boldsymbol{u}_{R}=\{\hat{\boldsymbol{p}}_{R},\hat{\boldsymbol{U}}_{R}\}$ for the resonator, and $\boldsymbol{u}_{i}=\{\hat{\boldsymbol{p}}_{i},\hat{\boldsymbol{U}}_{i}\}$ for the feedback inertance.

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 a system-wide linear stability analysis, is only meaningful if coupled with all of the other components in the engine, as discussed in the following.

The heat transfer and drag in the pulse tube, and the presence of gradients of base density and temperature in the REG/HX unit, require variations of entropy to be explicitly accounted for. Replacing the conservation equation for the total energy with the transport equation for entropy, expressed in terms of temperature and pressure using Gibbs’ relation, yields

(4.27a ) $$\begin{eqnarray}\displaystyle & \displaystyle \text{i}{\it\sigma}A\hat{{\it\rho}}+\frac{\text{d}\hat{U} }{\text{d}x}{\it\rho}_{0}+\hat{U} \frac{\text{d}{\it\rho}_{0}}{\text{d}x}=0, & \displaystyle\end{eqnarray}$$

(4.27b ) $$\begin{eqnarray}\displaystyle & \displaystyle \text{i}{\it\sigma}\hat{U} +\frac{A}{{\it\rho}_{0}}\frac{\text{d}\hat{p}}{\text{d}x}=-\frac{R_{C}}{{\it\rho}_{0}}\hat{U} , & \displaystyle\end{eqnarray}$$

(4.27c ) $$\begin{eqnarray}\displaystyle & \displaystyle {\it\rho}_{0}C_{p}\left[\text{i}{\it\sigma}\hat{T}+\frac{\hat{U} }{A}\frac{\text{d}T_{0}}{\text{d}x}\right]-\text{i}{\it\sigma}\hat{p}=-{\it\alpha}_{T}\hat{T}, & \displaystyle\end{eqnarray}$$

$t=0.55~\text{s}$

Recasting the system of equations in (4.27) in diagonalized form yields

(4.28a ) $$\begin{eqnarray}\displaystyle & \displaystyle \text{i}{\it\sigma}\hat{p}=\left[{\it\rho}_{0}R\mathscr{B}_{\hat{T}\hat{U} }-\frac{P_{0}}{A}\frac{\text{d}}{\text{d}x}-\frac{RT_{0}}{A}\frac{\text{d}{\it\rho}_{0}}{\text{d}x}\right]\hat{U} +[{\it\rho}_{0}R\mathscr{B}_{\hat{T}\hat{T}}]\hat{T}, & \displaystyle\end{eqnarray}$$

(4.28b ) $$\begin{eqnarray}\displaystyle & \displaystyle \text{i}{\it\sigma}\hat{U} =\left[-\frac{A}{{\it\rho}_{0}}\frac{\text{d}}{\text{d}x}\right]\hat{p}+\left[-\frac{R_{c}}{{\it\rho}_{0}}\right]\hat{U} , & \displaystyle\end{eqnarray}$$

(4.28c ) $$\begin{eqnarray}\displaystyle & \displaystyle \text{i}{\it\sigma}\hat{T}=\mathscr{B}_{\hat{T}\hat{U} }\hat{U} +\mathscr{B}_{\hat{T}\hat{T}}\hat{T}, & \displaystyle\end{eqnarray}$$

(4.29) $$\begin{eqnarray}\displaystyle & \displaystyle \mathscr{B}_{\hat{T}\hat{U} }=-\frac{T_{0}}{A}\left[\frac{1}{T_{0}}\frac{\text{d}T_{0}}{\text{d}x}+({\it\gamma}-1)\frac{\text{d}}{\text{d}x}\right], & \displaystyle\end{eqnarray}$$

(4.30) $$\begin{eqnarray}\displaystyle & \displaystyle \mathscr{B}_{\hat{T}\hat{T}}=-\frac{{\it\gamma}{\it\alpha}_{T}}{{\it\rho}_{0}C_{p}}, & \displaystyle\end{eqnarray}$$

(4.31) $$\begin{eqnarray}(\text{i}{\it\sigma}\unicode[STIX]{x1D644}-\boldsymbol{B}_{t})\boldsymbol{\cdot }\boldsymbol{u}_{t}=0,\end{eqnarray}$$

where $\boldsymbol{u}_{t}=\{\hat{\boldsymbol{p}}_{t},\hat{\boldsymbol{U}}_{t},\hat{\boldsymbol{T}}_{t}\}$ .

The complete eigenvalue problem can finally be solved by coupling the isolated component eigenvalue problems (4.22), (4.23), (4.25), (4.26) and (4.31) via the following conditions

(4.32a ) $$\begin{eqnarray}\displaystyle & \displaystyle \hat{U} _{R_{0}}=0, & \displaystyle\end{eqnarray}$$

(4.32b ) $$\begin{eqnarray}\displaystyle & \displaystyle \left.\frac{\text{d}}{\text{d}x}\hat{p}_{R}\right|_{0}=0, & \displaystyle\end{eqnarray}$$

(4.32c ) $$\begin{eqnarray}\displaystyle & \displaystyle \hat{p}_{J}=\hat{p}_{R_{1}}=\hat{p}_{i_{0}}=\hat{p}_{t_{1}}, & \displaystyle\end{eqnarray}$$

(4.32d ) $$\begin{eqnarray}\displaystyle & \displaystyle \hat{p}_{c}=\hat{p}_{i_{1}}=\hat{p}_{t_{0}}, & \displaystyle\end{eqnarray}$$

(4.32e ) $$\begin{eqnarray}\displaystyle & \displaystyle \hat{T}_{t_{1}}=\frac{1}{{\it\rho}_{0}R}\frac{{\it\gamma}-1}{{\it\gamma}}\hat{p}_{J}, & \displaystyle\end{eqnarray}$$

(4.32f ) $$\begin{eqnarray}\displaystyle & \displaystyle \hat{T}_{t_{0}}=\frac{1}{{\it\rho}_{0}R}\frac{{\it\gamma}-1}{{\it\gamma}}\hat{p}_{c}, & \displaystyle\end{eqnarray}$$

where $\boldsymbol{v}=\{\boldsymbol{u}_{t};\boldsymbol{u}_{R};\boldsymbol{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) and (4.23) replaces one corresponding equation in (4.33), therefore, not affecting the rank of the system. The eigenvalue structure is finally recovered by absorbing the equations that do not contain ${\it\sigma}$ (i.e. the ones deriving from (4.32)) via Gaussian elimination.

The linear modelling framework composed of equations (4.22)–(4.24) and (4.27) is first tested against the three-dimensional numerical simulation data for verification purposes. First, the acoustic impedance at different axial positions in the resonator, obtained by numerically integrating (4.24), has been quantitatively verified against the simulation data, as well as the constants $w_{c_{0}}$ in (4.22) and $w_{J_{0}}$ in (4.23), which are equal to $135\times 10^{6}~\text{N}~\text{m}^{-5}$ and $312\times 10^{6}~\text{N}~\text{m}^{-5}$ , respectively. Good agreement is also obtained by directly integrating (4.27) from section $t_{0}$ to $t_{1}$ and $i_{0}$ to $i_{1}$ (figure 8) using data from the numerical simulations as initial conditions (table 3). The integration has been carried out with the exact value of the frequency and growth rate extracted from the simulations (figure 9). Numerical trials have shown that, for a given initial condition and base state, the direct integration of (4.27) is much more sensitive to the angular frequency ${\it\omega}$ than the growth rate ${\it\alpha}$ , which suggests that the prediction of the latter, in thermoacoustic systems, is potentially problematic, especially within the framework of linear modelling in the spectral domain.

Figure 9. (a) Frequency, ${\it\omega}/2{\rm\pi}$ , and (b) growth rate, ${\it\alpha}$ , and limit-cycle pressure amplitude, $p_{lc}^{\prime }$ , versus hot-to-cold temperature ratio, ${\it\tau}=T_{h}/T_{c}$ . Linear stability model (- - -), numerical simulations (symbols) with corresponding fitting (—–) yielding the critical temperature ratio ${\it\tau}_{cr}=1.505$ and $\left.p_{lc}^{\prime }\right|_{{\it\delta}{\it\tau}=1}=41\,000~\text{Pa}$ .

Table 3. Comparison between linear theory and simulation data for $T_{h}=460~\text{K}$ and $T_{h}=500~\text{K}$ extracted at locations shown in figure 8 at $t=0.55~\text{s}$ . The linearized equations (4.27) are integrated from $t_{0}$ to $t_{1}$ and $i_{0}$ to $i_{1}$ (figure 8) with initial conditions (i.c.) at $t_{0}$ and $i_{0}$ , respectively, and ${\it\alpha}$ , ${\it\omega}$ , $T_{0}(x)$ and ${\it\rho}_{0}(x)$ taken from the numerical simulations (figures 5 a, 9).

The eigenvalues ${\it\sigma}$ are finally calculated by directly solving the complete eigenvalue problem for operating conditions ranging between ${\it\tau}=1.25$ and ${\it\tau}=1.85$ , where ${\it\tau}=T_{h}/T_{c}$ is the hot-to-cold temperature ratio. The segments representing the pulse tube, inertance and resonator were discretized with 256, 64 and 32 points, respectively, with a fourth-order spatial polynomial reconstruction. The frequency of instability and its variation with ${\it\tau}$ are predicted within a ${\sim}0.2~\text{Hz}$ error (figure 9 a). The operating frequency of the system could also be predicted (with a ${\sim}1~\text{Hz}$ error) by simply solving the eigenvalue problem (4.24) in the complete variable-area resonator alone (without the pulse tube), in accordance with the phase distribution shown in figure 7(b), which is consistent with simple standing-wave resonance. The growth rate is slightly overpredicted (figure 9 b), having neglected viscous and nonlinear losses. Overall, the quantitative agreement is very encouraging, serving both as a verification step for the full Navier–Stokes calculations and to gain insight into the nature of the instability, also briefly discussed in the following section.

One could also think of solving the complete eigenvalue problem by iteratively integrating (4.27) and (4.28) in space, starting from a given set of initial conditions or guesses. This approach is adopted in DeltaEC (Ward & Swift Reference Ward and Swift1994) to predict the limit-cycle pressure and velocity distributions in TAEs, which is a valid approximation in the case of relatively low-pressure amplitudes, limited waveform distortion and simple geometries. Several numerical trials, however, have shown that applying the same approach to the eigenvalue problem above (valid for the startup phase only) leads to numerically unstable results, especially when trying to predict the growth rate. Accurately predicting the correct instability frequency, on the other hand, appears to be possible regardless of the quality of the specific strategy adopted. A fully implicit spatial formulation, similar to Helmholtz solvers used in investigations of thermoacoustically unstable reactive flows (Poinsot & Veynante Reference Poinsot and Veynante2011), which has been adopted in the present context, is the only one that has proved to be reliable.

4.4. Supercritical Hopf bifurcation

A linear fit of the growth rates, ${\it\alpha}$ , extracted from the numerical simulation data versus the temperature ratio ${\it\tau}=T_{h}/T_{c}$ (figure 9 b) suggests a critical value of ${\it\tau}_{cr}=1.505$ . This is in perfect agreement with the same value obtained by fitting the functional form suggested by the supercritical Hopf bifurcation model

(4.34) $$\begin{eqnarray}p_{lc}^{\prime }=\left.p_{lc}^{\prime }\right|_{{\it\delta}{\it\tau}=1}\sqrt{{\it\tau}-{\it\tau}_{cr}},\end{eqnarray}$$

to the limit-cycle pressure amplitudes $p_{lc}^{\prime }$ versus ${\it\tau}$ . A similar result is obtained via nonlinear modelling by Mariappan & Sujith (Reference Mariappan and Sujith2011). The fitting parameters in (4.34) are $\left.p_{lc}^{\prime }\right|_{{\it\delta}{\it\tau}=1}$ (dimensional) and ${\it\tau}_{cr}$ (dimensionless). Moreover, two assumptions required by the Hopf bifurcation theorem are also satisfied: the non-hyperbolicity condition, ${\it\alpha}=0$ and ${\it\omega}\neq 0$ at ${\it\tau}_{cr}$ (figure 9 a), and the transversality condition, $\text{d}{\it\alpha}/\text{d}{\it\tau}\neq 0$ at ${\it\tau}_{cr}$ (figure 9 b). These results have important implications on the parametrization of nonlinear fluxes, discussed later in § 5.2.

5.1. Direct modelling of acoustic streaming

A triple decomposition can be invoked to separate the streaming flow (rigorously defined below) from the acoustic field and the small-scale high-frequency fluctuations, starting with the Reynolds decomposition

(5.1a ) $$\begin{eqnarray}\displaystyle & {\it\rho}={\it\rho}_{0}+{\it\rho}^{\prime }, & \displaystyle\end{eqnarray}$$

(5.1b ) $$\begin{eqnarray}\displaystyle & p=p_{0}+p^{\prime }, & \displaystyle\end{eqnarray}$$

(5.1c ) $$\begin{eqnarray}\displaystyle & u_{i}=u_{0,i}+u^{\prime }, & \displaystyle\end{eqnarray}$$

(5.2) $$\begin{eqnarray}u_{0,i}=\overline{u}_{i}=\int _{-\infty }^{\infty }u_{i}(\boldsymbol{x},t+{\it\tau})\frac{\text{sin}({\rm\pi}f_{c}{\it\tau})}{{\rm\pi}{\it\tau}}\text{d}{\it\tau},\end{eqnarray}$$

where $f_{c}$ is the cutoff frequency. The filtering operation (5.2) is, in practice, carried out over six acoustic periods by adopting Simpson’s quadrature rule on the discrete data sampled at 2.2 kHz and $f_{c}=0.9f$ , where $f={\it\omega}/2{\rm\pi}$ is the acoustic frequency (figure 9 a). The remainder of the filtering operation in (5.1) can be further decomposed into a purely acoustic (subscript ‘ $a$ ’) and a small-scale component (subscript ‘ $t$ ’),

(5.3a ) $$\begin{eqnarray}\displaystyle & {\it\rho}^{\prime }={\it\rho}_{a}^{\prime }+{\it\rho}_{t}^{\prime }, & \displaystyle\end{eqnarray}$$

(5.3b ) $$\begin{eqnarray}\displaystyle & p^{\prime }=p_{a}^{\prime }+p_{t}^{\prime }, & \displaystyle\end{eqnarray}$$

(5.3c ) $$\begin{eqnarray}\displaystyle & u_{i}^{\prime }=u_{a,i}^{\prime }+u_{t,i}^{\prime }, & \displaystyle\end{eqnarray}$$

$f$

$f_{c}>10f$

No evidence of significant turbulent activity has been found in the resonator. The Stokes Reynolds number in the neck is approximately 750, barely entering the intermittently turbulent regime (Jensen, Sumer & Fredsøe Reference Jensen, Sumer and Fredsøe1989). The low Stokes Reynolds number, together with other complex geometrical features, does not allow turbulence to be fully developed during the acoustic cycle.

Vorticity contours from instantaneous visualizations (figure 10) suggest that the most intense small-scale fluctuations occur in the feedback inertance. The unsteadiness of the (larger-scale) acoustic fluctuations does not allow turbulence to reach a fully (or even partially) developed state. The Reynolds number based on the Stokes thickness ${\it\delta}_{{\it\nu}}=\sqrt{2{\it\nu}/{\it\omega}}$ and the maximum velocity amplitude at the centre of the feedback inertance is approximately $Re_{{\it\delta}_{{\it\nu}}}=360$ (far below the upper limit of the disturbed laminar regime described in Jensen et al. (Reference Jensen, Sumer and Fredsøe1989)), suggesting that the turbulent kinetic energy generated from the breakup of the vortices rolling up from the edges of the annular tube is not sustained by the wall-induced shear. Moreover, turbulent stresses extracted for the $T_{h}=500~\text{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 need to devise a simple predictive modelling framework for the streaming velocity (5.7), which is discussed in the following.

Figure 11. Profiles of cycle- and time-averaged variance of acoustic $\langle \overline{u_{a,x}^{\prime }u_{a,x}^{\prime }}\rangle$ (—–) and small-scale $\langle \overline{u_{t,x}^{\prime }u_{t,x}^{\prime }}\rangle$ axial velocity fluctuations (– –) extracted in the inertance, for $x/L_{i}=0.01,0.2,0.40,0.60,0.80,0.99$ , from left to right, where $L_{i}=0.204~\text{m}$ (shifted by $500~\text{m}^{2}~\text{s}^{-2}$ for clarity) for $T_{h}=500~\text{K}$ and grid C (table 2).

Substituting the decomposition (5.1) into the time-filtered conservation of mass, ignoring temporal and spatial variations of the filtered density field (a strong assumption, particularly for the regions around the sharp edges and in the TBT), assuming that second-order quantities in the small-scale fluctuations are negligible with respect to their acoustic counterpart (figure 11),

(5.4) $$\begin{eqnarray}\displaystyle & \overline{u_{t,i}^{\prime }u_{t,i}^{\prime }}\ll \overline{u_{a,i}^{\prime }u_{a,i}^{\prime }}, & \displaystyle\end{eqnarray}$$

(5.5) $$\begin{eqnarray}\displaystyle & \overline{{\it\rho}_{t}^{\prime }u_{t,i}^{\prime }}\ll \overline{{\it\rho}_{a}^{\prime }u_{a,i}^{\prime }}, & \displaystyle\end{eqnarray}$$

(5.6) $$\begin{eqnarray}\frac{\partial \tilde{u} _{i}}{\partial x_{i}}=0,\end{eqnarray}$$

where $\tilde{u} _{i}$ is the density-weighted velocity field, $\tilde{u} _{i}=\overline{{\it\rho}u_{i}}/{\it\rho}_{0}$ , defined based on the filtering operation (5.2), which, under the assumptions made, can be expressed as

(5.7) $$\begin{eqnarray}\tilde{u} _{i}\simeq u_{0,i}+\frac{\overline{{\it\rho}_{a}^{\prime }u_{a_{i}}^{\prime }}}{{\it\rho}_{0}}.\end{eqnarray}$$

In the present manuscript the density-weighted velocity field, $\tilde{u} _{i}$ , is adopted as the definition of the streaming velocity based on (5.7), which is second-order accurate in wave amplitude.

The streaming velocity field in our case is axially symmetric and spans the full extent of the engine (figure 12). Very large and elongated recirculations, of the order of a quarter of the acoustic wavelength, are visible in the resonator and are driven by the wall-normal gradient of the wave-induced shear stresses (not shown). Large recirculations of the order of the resonator radius near the sharp edges of the pulse tube and a mean flow circulating around the pulse tube (following the direction of the amplified waves) are also observed. The latter is called Gedeon (or DC) streaming (Gedeon Reference Gedeon1997) and is responsible for the mean transport of hot fluid away from the HHX, in the direction of the AHX2. As a result of such advective process, the TBT experiences a significant thermal load, with hot fluid occupying, on average, a significant portion of the TBT’s length at the limit cycle (figure 10). Such a process, leading to heat leakage, limits the efficiency of most travelling-wave TAEs, as also discussed later in the context of the present engine.

Figure 12. Contours of axial component of time-averaged streaming velocity $\langle \tilde{u} _{x}\rangle$ (5.7) for $T_{h}=500~\text{K}$ and grid C (table 2). Full-scale visualization (a) and zoom on the right end (b). Results have been mirrored about the centreline and streamlines have been only numerically approximated for illustrative purposes.

Assuming that time scales of variation of the filtered quantities ${\it\rho}_{0}$ and $u_{0,i}$ are much longer than the acoustic period, and under the same assumptions underlying the derivation of (5.7), it can be shown that $\tilde{u} _{i}$ satisfies the incompressible Navier–Stokes equations (Rudenko & Soluyan Reference Rudenko and Soluyan1977),

(5.8) $$\begin{eqnarray}\frac{\partial \tilde{u} _{i}}{\partial t}+\frac{\partial }{\partial x_{j}}\tilde{u} _{i}\tilde{u} _{j}+\frac{\partial p_{0}}{\partial x_{i}}-{\it\nu}{\rm\nabla}^{2}\tilde{u} _{i}=F_{a,i},\end{eqnarray}$$

where $p_{0}=P_{0}/{\it\rho}_{0}$ and the forcing term $F_{a,i}$ is the divergence of the wave-induced Reynolds stresses, which can be expressed to second-order accuracy in wave amplitude as

(5.9) $$\begin{eqnarray}F_{a,i}=-\frac{\partial }{\partial x_{j}}\overline{u_{a,j}^{\prime }u_{a,i}^{\prime }}+\frac{\partial }{\partial x_{j}}\left\{-\frac{{\it\nu}}{{\it\rho}_{0}}\left[\frac{\partial }{\partial x_{j}}\overline{{\it\rho}_{a}^{\prime }u_{a,i}^{\prime }}+\frac{\partial }{\partial x_{i}}\overline{{\it\rho}_{a}^{\prime }u_{a,j}^{\prime }}\right]+\frac{2}{3}\frac{{\it\nu}}{{\it\rho}_{0}}\frac{\partial }{\partial x_{k}}\overline{{\it\rho}_{a}^{\prime }u_{a,k}^{\prime }}{\it\delta}_{ij}\right\}.\end{eqnarray}$$

For large Reynolds numbers based on the streaming velocity, the terms containing the molecular diffusion in (5.9) can be neglected, finally yielding

(5.10) $$\begin{eqnarray}F_{a,i}\simeq -\frac{\partial }{\partial x_{j}}\overline{u_{a,j}^{\prime }u_{a,i}^{\prime }}.\end{eqnarray}$$

The maximum value of the stresses $\overline{u_{a,j}^{\prime }u_{a,i}^{\prime }}$ is expected to be found in the feedback inertance where the acoustic power is maximum in the system. A grid-sensitivity study on all of the components of the stresses in the inertance (figure 13) shows monotonic grid convergence of the stresses from grid A to C (table 2).

Figure 13. Profiles of time-averaged streaming velocity (a) and wave-induced Reynolds stresses (b–d) in the inertance, for $x/L_{i}=0.01,0.2,0.40,0.60,0.80,0.99$ , respectively from left to right, where $L_{i}=0.204~\text{m}$ (shifted by $20~\text{m}~\text{s}^{-1}$ or $150~\text{m}^{2}~\text{s}^{-2}$ for clarity). Data for $T_{h}=500~\text{K}$ and grid A (— ⋅ —), grid B (– –) and grid C (—–) (table 2). The time-averaged volumetric flux through the inertance (intensity of Gedeon streaming) for this drive ratio is approximately $0.0033~\text{m}^{3}~\text{s}^{-1}$ , corresponding to a mean streaming velocity of $1.2~\text{m}~\text{s}^{-1}$ .

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 large aforementioned recirculations. On the other hand, Gedeon streaming is driven by the viscous decay of the wave amplitude in the annular inertance. This results in a negative axial gradient of normal stress $\overline{u_{a,x}^{\prime }u_{a,x}^{\prime }}$ , visible in both figures 14 and 13(b).

Figure 14. Contour plots of the divergence of wave-induced Reynolds stresses. (a) Axial, $F_{a,x}$ , and (b) radial, $F_{a,r}$ , components extracted for $T_{h}=500~\text{K}$ on grid C (table 2).

An axially symmetric numerical model, $\text{Stream}^{\text{X}}$ (for more details see appendix A), has been developed to directly simulate the streaming velocity field as the solution of the incompressible equation (5.8) driven by the divergence of the wave-induced stresses (figure 14) extracted from the three-dimensional fully compressible calculations. Secondary features such as steady large-scale recirculations near the sharp edges of the annular tube are only qualitatively reproduced (figure 15), with the appearance of a second recirculation in the compliance, which is not observed in the calculations. The actual target of the present low-order modelling effort is the prediction of the intensity of the Gedeon streaming. In spite of the numerical challenges involved in solving incompressible flow in the presence of a sharp edge, the latter is predicted fairly accurately (figure 17 a), especially for low drive ratios. For high drive ratios, resulting in very high limit-cycle acoustic amplitudes, the errors associated with the assumptions made in deriving (5.7), (5.8) and (5.10) become too severe. ‘Slow’ streaming (Rudenko & Soluyan Reference Rudenko and Soluyan1977) never actually occurs in our engine, where the maximum intensity of $\tilde{u}$ is comparable to the acoustic velocity amplitude in all cases.

Figure 15. Positive (—–) and negative (– –) iso-levels of the Stokes streamfunction from the incompressible solver $\text{Stream}^{\text{X}}$ driven by the wave-induced Reynolds stresses (figure 14) extracted from the fully compressible three-dimensional calculations for $T_{h}=500~\text{K}$ on grid C (table 2).

5.2. Efficiency and energy fluxes in the TBT

As the acoustic energy grows during the initial transient, so does the intensity of the Gedeon streaming, increasing the rate of advective transport of hot fluid away from the HHX towards the AHX2. This results in unwanted heat leakage, which lowers the overall energy conversion efficiency. The gradual expansion of the gas in the TBT determines a slow increase of the background pressure in the system (not shown), which plateaus only when the hot temperature front enters into contact with the AHX2. At this point, a rapid increase of the growth rate is also observed, as shown by the kink in the time series in figure 4. This is due to a change in the structure of the velocity field in the TBT after the first contact between the hot gas front and the AHX2 occurs. As shown in figure 10, and in time series of instantaneous axial velocity in the TBT (not shown), when the hot gas is displaced towards the AHX2, the hot flow front is abruptly stopped by the sudden density increase caused by the cooling in the AHX2. This results in a reduction of the cycle-averaged velocity and, therefore, a temporary containment of the intensity of the Gedeon streaming and heat leakage from the HHX.

A limit cycle is only reached later, with a constant background pressure, when the losses in the system balance the enhanced acoustic energy production. These losses include streaming in the resonator and dissipation associated with the turbulent vortex shedding from the pulse tube walls.

The exact conservation equation for the density-averaged internal energy, $\tilde{e}=\overline{{\it\rho}e}/\overline{{\it\rho}}$ , reads (Lele Reference Lele1994)

(5.11) $$\begin{eqnarray}\frac{\partial }{\partial t}(\overline{{\it\rho}}\tilde{e})=-\frac{\partial }{\partial x_{j}}\left(\overline{{\it\rho}}\left[\tilde{u} _{j}\tilde{e}+\widetilde{h^{\prime \prime }u_{j}^{\prime \prime }}\right]-\overline{q}_{j}\right)+\frac{\partial }{\partial x_{j}}\left(\overline{p^{\prime }u_{j}^{\prime \prime }}\right)+\overline{u_{i}^{\prime \prime }}\frac{\partial \overline{p}}{\partial x_{i}}-\overline{p^{\prime }\frac{\partial u_{i}^{\prime \prime }}{\partial x_{i}}},\end{eqnarray}$$

where $u_{j}^{\prime \prime }=u_{i}-\tilde{u} _{i}$ and $h^{\prime \prime }=h_{i}-\tilde{h}$ are the fluctuations of the density-weighted averages of velocity and enthalpy, and $\overline{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 terms and assuming equilibrium conditions, yields

(5.12) $$\begin{eqnarray}\underbrace{\langle \,\overline{{\it\rho}}\tilde{u} \tilde{e}\rangle }_{\text{Advective H.T.}}+\underbrace{\left\langle \,\overline{{\it\rho}}\widetilde{h^{\prime \prime }u^{\prime \prime }}\right\rangle }_{\text{Thermoacoustic H.T.}}-\underbrace{\left\langle \,\overline{p^{\prime }u^{\prime \prime }}\right\rangle }_{\text{Acoustic energy flux}}\simeq \text{const.},\end{eqnarray}$$

which is verified with fairly good approximation in the simulations (figure 16 b). The intensity of the advective heat transport in the TBT is proportional to the mean temperature profile, since the streaming velocity is uniform in this region (figure 12). The adjustment length back to ambient temperature of the mean temperature distribution increases with the drive ratio (figure 16 a). This is due to thermoacoustic heat transport mechanisms. As the hot fluid front is transported with stronger intensity towards the AHX2 by the high-amplitude velocity fluctuations, steeper instantaneous temperature gradients form at the interface between the AHX2 and the TBT (figure 10). The result is a net cycle-average conductive heat flux in the positive axial direction, creating a temperature buffer region. The intensity of the conductive heat flux in (5.11) is, however, negligible compared to the quantities in (5.12), which dominate the energy transport budget in the TBT.

Figure 16. (a) Axial distribution of mean temperature for $T_{h}=460~\text{K}$ , $T_{h}=480~\text{K}$ and $T_{h}=500~\text{K}$ and (b) surface integrated energy fluxes in (5.12) for grid C and $T_{h}=500~\text{K}$ at limit cycle. Acoustic power (- -), thermoacoustic heat transport, (– –), advective heat transport (— ⋅ —) and overall sum (○).

An accurate evaluation of the other terms in (5.11) is, unfortunately, made impractical by the (necessary) application of a second-order ENO reconstruction in the TBT and the smoothing associated with azimuthally averaging the three-dimensional unstructured data. The energy balance expressed by (5.12) is, however, satisfied to a sufficient degree of accuracy to gain insight into the role of the Gedeon streaming in determining the overall 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 piezo electric element or a linear alternator) has been included in the set-up investigated, a metric for efficiency can still be defined as

(5.13) $$\begin{eqnarray}{\it\eta}_{H}=\frac{\left\langle \,\overline{p^{\prime }u^{\prime \prime }}\right\rangle }{\left\langle \,\overline{{\it\rho}}\widetilde{h^{\prime \prime }u^{\prime \prime }}\right\rangle +\langle \,\overline{{\it\rho}}\tilde{u} \tilde{e}\rangle },\end{eqnarray}$$

which is the ratio between acoustic energy produced and total heat dissipated by the AHX2. For example, our theoretical device produces ${\sim}0.2~\text{kW}$ of acoustic power at $T_{h}=500~\text{K}$ while losing ${\sim}2.0~\text{kW}$ , mostly due to mean advection caused by Gedeon streaming, achieving a modest efficiency ( ${\sim}10\,\%$ ). Realistic travelling-wave TAEs can reach overall efficiencies of ${>}20\,\%$ , or even ${>}30\,\%$ if built in a cascaded configuration (Gardner & Swift Reference Gardner and Swift2003). The efficiency directly evaluated from the simulation data (figure 17 a) decreases rapidly with the temperature ratio, which could be expected in TAEs with excessive Gedeon streaming controlling the energy balance in the TBT (G.W. Swift, Personal communication, 2014). However, the uncertainties in the integral quantities in figure 17 b, due to averaging over a limited number (approximately 25) of acoustic cycles, makes the direct metric for the efficiency (5.13) not very reliable.

Figure 17. (a) Intensity of Gedeon streaming versus hot-to-cold temperature ratio for fully compressible Navier–Stokes simulations (▫) with corresponding fit (——), incompressible model (▪), order-of-magnitude analysis (5.16) (▹), efficiency ${\it\eta}_{H}$ (5.13) (

). (b) Average nonlinear energy fluxes at limit cycle, full Navier–Stokes simulations (symbols) and corresponding fit based on (5.18), (5.20), (5.22) (lines), respectively, for advective heat transport (

, — ⋅ —), acoustic power (▵, - - - ), and thermoacoustic heat transport (♢, – –).

Figure 18. (a) Illustration of balance (5.16) between travelling-wave decay in feedback inertance and drag due to HXs and REG. (b) Streamwise profile of axial wave-induced stresses extracted at $r=56~\text{mm}$ for $T_{h}=460~\text{K}$ (a), $T_{h}=480~\text{K}$ (b) and $T_{h}=500~\text{K}$ (c) with linear fitting (- - -) approximating the mean decay rate.

A more robust estimate for ${\it\eta}_{H}$ can be derived by investigating the scaling of the volume-averaged energy fluxes in (5.12) (figure 17 b) with a simple order-of-magnitude analysis and curve fitting. The advective flux in the TBT, driven by the Gedeon streaming, is expected to scale as

(5.14) $$\begin{eqnarray}\langle \overline{{\it\rho}}\tilde{u} \tilde{e}\rangle \sim {\it\rho}_{0,tbt}u_{dc}C_{v}T_{h},\end{eqnarray}$$

where ${\it\rho}_{0,tbt}\sim p_{0}/(RT_{c}{\it\tau})$ is the average density in the TBT and the hot temperature is simply $T_{h}=T_{c}{\it\tau}$ . The analytical expression of the streaming velocity due to a freely propagating travelling wave,

(5.15) $$\begin{eqnarray}u_{dc}=\frac{1}{2}\frac{|u_{lc}^{\prime }|^{2}}{a_{0}},\end{eqnarray}$$

suggests that the intensity of the streaming velocity scales quadratically with the limit-cycle velocity amplitude $u_{lc}^{\prime }\sim p_{lc}^{\prime }/({\it\rho}_{0}a_{0})$ . While the quadratic scaling of the streaming velocity is an expected result, a more quantitative estimate for $u_{dc}$ can be obtained by further simplifying the analysis in § 5.1, where the streaming velocity is directly modelled based on an hydrodynamic analogy. In fact, by roughly measuring the average spatial decay rate of the axial acoustic stresses in the inertance, $\langle F_{a,x}\rangle _{i}$ (figure 18 a), 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

(5.16) $$\begin{eqnarray}u_{dc}\sim \frac{{\it\rho}_{0}A_{i}L_{i}\langle F_{a,x}\rangle _{i}}{V_{\mathit{AHX2}}R_{c,\mathit{AHX2}}+V_{\mathit{HHX}}R_{c,\mathit{HHX}}+V_{\mathit{REG}}R_{c,\mathit{REG}}+V_{\mathit{AHX}}R_{c,\mathit{AHX}}},\end{eqnarray}$$

where $V$ and $R_{c}$ are, respectively, the volume and drag coefficients (obtained by linearizing (2.3)) of the REG/HX unit in the pulse tube. The estimate (5.16) is in good quantitative agreement with the results from § 5.1 (figure 17 a), showing that, in this case, the viscous wave-amplitude decay in the inertance, leading to the generation of strong divergence of wave-induced Reynolds stresses, dominates over mean pressure gradient ( $\partial p_{0}/\partial x_{i}$ in (5.8)) effects. The case for $T_{h}=480~\text{K}$ is an outlier simply due to the lack of grid convergence of the acoustic stresses for this particular case (table 2). The quadratic scaling adopted in (5.15) is, therefore, further justified by (5.16), where $\langle F_{a,x}\rangle _{i}\sim {u_{lc}^{\prime }}^{2}$ . By invoking the previously derived scaling for the limit-cycle pressure (4.34), (5.15) becomes

(5.17) $$\begin{eqnarray}u_{dc}=A_{dc}\frac{1}{2}\frac{1}{{\it\rho}_{0}^{2}a_{0}^{3}}\left.{p_{lc}^{\prime }}^{2}\right|_{{\it\delta}{\it\tau}=1}({\it\tau}-{\it\tau}_{cr}),\end{eqnarray}$$

where the fitting coefficient is $A_{dc}=0.52$ (figure 17 b). Substituting (5.17) into (5.14) yields

(5.18) $$\begin{eqnarray}\langle \,\overline{{\it\rho}}\tilde{u} \tilde{e}\rangle =A_{adv}\left[\frac{C_{v}P_{0}}{R}A_{dc}\frac{1}{2}\frac{1}{{\it\rho}_{0}^{2}a_{0}^{3}}\left.{p_{lc}^{\prime }}^{2}\right|_{{\it\delta}{\it\tau}=1}\right]({\it\tau}-{\it\tau}_{cr}),\end{eqnarray}$$

where the fitting coefficient is $A_{adv}=0.51$ . The same procedure can be applied to the acoustic energy flux, which scales as

(5.19) $$\begin{eqnarray}\langle \,\overline{p^{\prime }u^{\prime \prime }}\rangle \sim p_{lc}^{\prime }u_{lc}^{\prime }{\textstyle \frac{1}{2}}\,\text{cos}({\rm\Delta}{\it\phi}),\end{eqnarray}$$

where ${\rm\Delta}{\it\phi}$ is the phase difference between pressure and velocity observed in the REG/HX (figure 5 c), leading to

(5.20) $$\begin{eqnarray}\langle \,\overline{p^{\prime }u^{\prime \prime }}\rangle =A_{ap}\frac{1}{2}\,\text{cos}({\rm\Delta}{\it\phi})\frac{1}{{\it\rho}_{0}a_{0}}\left.{p_{lc}^{\prime }}^{2}\right|_{{\it\delta}{\it\tau}=1}({\it\tau}-{\it\tau}_{cr}),\end{eqnarray}$$

where the fitting constant is $A_{ap}=0.3$ . Finally, the thermoacoustic heat flux is expected to scale as

(5.21) $$\begin{eqnarray}\langle \,\overline{{\it\rho}}\widetilde{h^{\prime \prime }u^{\prime \prime }}\rangle \sim {\it\rho}_{0,tbt}C_{p}\frac{\partial T}{\partial P}p_{lc}^{\prime }\frac{p_{lc}^{\prime }}{{\it\rho}_{0}a_{0}}\frac{1}{2}\,\text{cos}({\rm\Delta}{\it\phi}),\end{eqnarray}$$

where ${\rm\Delta}{\it\phi}$ is needed to account for the correlation between pressure and velocity (like in (5.19)), resulting in

(5.22) $$\begin{eqnarray}\langle \,\overline{{\it\rho}}\widetilde{h^{\prime \prime }u^{\prime \prime }}\rangle =A_{ta}\left[\frac{1}{2{\it\rho}_{0}a_{0}}\left.{p_{lc}^{\prime }}^{2}\right|_{{\it\delta}{\it\tau}=1}\frac{1}{2}\,\text{cos}({\rm\Delta}{\it\phi})\right]({\it\tau}-{\it\tau}_{cr})\end{eqnarray}$$

with the fitting coefficient $A_{ta}=0.93$ . With all fluxes scaling as ${\sim}({\it\tau}-{\it\tau}_{cr})$ , this leads to a robust estimate for the efficiency,

(5.23) $$\begin{eqnarray}{\it\eta}_{H,\mathit{avrg}}=0.13,\end{eqnarray}$$

effectively averaged over the range of temperature ratios investigated.

Linear acoustic solvers applied at the limit cycle can directly estimate second-order quantities in the acoustic amplitudes such as (5.19) and (5.22), or even the mean wave-amplitude decay $\langle F_{a,x}\rangle _{i}$ due to viscous losses in a duct, without prior knowledge of the critical temperature ratio, ${\it\tau}_{cr}$ . However, process-based parametrizations for the Gedeon streaming, and therefore the advective transport (5.18), are still currently missing despite having a first-order impact on the acoustic energy budgets and the efficiency. We have shown that, for slow streaming, commonly found in realistic travelling-wave TAEs, a hydrodynamic analogy (Lighthill Reference Lighthill1978) can be invoked, and further simplified, leading to a very low-order modelling approach (5.16), which is very amenable in the 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 ${\it\tau}_{cr}$ , making the hydrodynamic analogy investigated in § 5.1 an attractive modelling paradigm for streaming.