2.1 Design of the minimal-unit (or single-pore) thermoacoustic resonator model

The proposed computational set-up (figure 1, top) is a straight, two-dimensional, axially periodic minimal-unit (or single-pore) thermoacoustic device composed of four constant-area sections (a–d). Such configuration represents an idealization of a looped thermoacoustic resonator (figure 1, bottom) similar – but not identical – to the one adopted by Yazaki et al. (Reference Yazaki, Iwata, Maekawa and Tominaga1998). Adiabatic slip conditions are applied everywhere, with the exception of the thermoacoustic regenerator, or core (section $b$ ), where isothermal no-slip walls are used to impose a linear wall temperature distribution $T_{w}(x)$ from the cold, $T_{C}$ , to the hot-side temperature, $T_{H}$ , resulting in the base temperature distribution, $T_{0}(x)$ . A body force is added to suppress Gedeon streaming (§ 2.2.1), which would otherwise cause convective heat transport away from the hot end of section $b$ and require the introduction of a thermal buffer tube and a secondary ambient heat exchanger (Penelet et al. Reference Penelet, Gusev, Lotton and Bruneau2005), thus introducing further complications in the proposed canonical set-up (discussed in § 2.2.1). The resulting relaxation of base state from $T_{H}$ to $T_{C}$ outside the regenerator is due to molecular diffusion. Consequently, a diffusive thermal layer develops which is very thin compared to the acoustic wave length and is neglected in the subsequent analysis. The uniform base pressure and cold-side values of density and temperature are set to be equal to the reference thermodynamic quantities $P_{0}=P_{ref}$ , $\unicode[STIX]{x1D70C}_{C}=\unicode[STIX]{x1D70C}_{ref}$ , $T_{C}=T_{ref}$ (table 1), chosen for air (§ 2.2.1). Stacking of any number of thus-conceived single-pore models in the $y$ direction, i.e. preserving the area ratio and increasing the number of pores in the regenerator, would yield the same numerical results.

The minimal-unit choice is dictated by the need to maximize the numerical resolution in the propagation direction of the captured shocks. Although full resolution of the propagating shocks is still not computationally feasible due to the very large length of the set-up (of the order of 1 m) compared to the typical shock-thickness scale (of the order of $100~\unicode[STIX]{x03BC}\text{m}$ ). Performing a fully resolved three-dimensional simulation of an equivalent experimental set-up would be (even more) unfeasible: for instance, the set-up studied by Yazaki et al. (Reference Yazaki, Iwata, Maekawa and Tominaga1998) consists of two heat exchangers and a ceramic catalyst with approximately 1000 pores. By design, the minimal-unit configuration neglects the thermoviscous losses outside the regenerator that attenuate the thermoacoustic amplification. The two-dimensional effects of the curvature of the resonator walls are also neglected. The former become less important as the ratio of the pore diameter to resonator diameter and number of pores increases (42 and 1000, respectively, in Yazaki’s set-up), and the latter are only relevant for very small ratios of the curvature radius to acoustic wavelength ( ${\sim}0.3$ in Yazaki’s set-up, assuming curvature radius of $L/6$ where $L$ is the total length). In conclusion, performing a full three-dimensional simulation taking into account approximately 1000 thermoacoustic pores, viscous losses in the resonator, and the curvature effects, would not lead to any significant additional insights into the physics of thermoacoustically generated shock waves, especially within the theoretical scope of the present study.

Table 1.

Thermodynamic parameters for base state. Base pressure and cold-side values of density and temperature are set to reference values: $P_{0}=P_{ref}$ , $\unicode[STIX]{x1D70C}_{C}=\unicode[STIX]{x1D70C}_{ref}$ , $T_{C}=T_{ref}$ .

The total length of the device is fixed to $\ell _{a}+\ell _{b}+\ell _{c}+\ell _{d}=2.58~\text{m}$ , taken from the experimental set-up of Yazaki et al. (Reference Yazaki, Iwata, Maekawa and Tominaga1998). The height of the regenerator has been chosen such that $h_{b}\sim 2\unicode[STIX]{x1D6FF}_{k}$ , where

(2.1) $$\begin{eqnarray}\unicode[STIX]{x1D6FF}_{k}=\sqrt{\frac{2\unicode[STIX]{x1D708}}{\unicode[STIX]{x1D714}Pr}},\end{eqnarray}$$

$\unicode[STIX]{x1D708}$ is the kinematic viscosity, $Pr$ is the Prandtl number and $\unicode[STIX]{x1D714}$ is the angular frequency of the unstable mode. This results in a porosity $h_{b}/h_{c}=0.91$ . The area ratio $h_{d}/h_{c}$ and the length $\ell _{c}$ have been chosen to yield sub-optimal values of growth rates across all hot-side temperature settings (figure 2), assuring high enough thermoacoustic instability to achieve rapid shock-wave formation. This optimization has been carried out with the system-wide numerical approach developed by Lin et al. (Reference Lin, Scalo and Hesselink2016), using the governing equations outlined in § 4.1. Further justification of the set-up choice and a comparison of the linear numerical analysis with Yazaki et al.’s (Reference Yazaki, Iwata, Maekawa and Tominaga1998) stability data are given in appendix A.

Figure 2.

Iso-contours of thermoacoustic growth rate, $\unicode[STIX]{x1D6FC}$ , of the second harmonic versus resonator area ratio, $h_{d}/h_{c}$ , and length $\ell _{c}$ (see figure 1) for hot-side temperatures $T_{H}=400~\text{K}$ (a), $T_{H}=450~\text{K}$ (b), $T_{H}=500~\text{K}$ (c) and $T_{H}=550~\text{K}$ (d). (——), $\unicode[STIX]{x1D6FC}>0$ ; (- - - -), $\unicode[STIX]{x1D6FC}<0$ ; (●), sub-optimal values $\ell _{c}=0.59~\text{m}$ and $h_{d}/h_{c}=1.28$ chosen for present investigation (figure 1). Growth rate values for each $T_{H}$ at sub-optimal geometry are listed in the table.

2.2 Navier–Stokes calculations

2.2.1 Governing equations

Fully compressible Navier–Stokes simulations are carried out by solving the conservation laws for mass, momentum and total energy in two dimensions, given by

(2.2a ) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}t}(\unicode[STIX]{x1D70C})+\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}x_{j}}(\unicode[STIX]{x1D70C}u_{j})=0, & \displaystyle\end{eqnarray}$$

(2.2b ) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}t}(\unicode[STIX]{x1D70C}u_{i})+\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}x_{j}}(\unicode[STIX]{x1D70C}u_{i}u_{j})=-\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}x_{i}}p+\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}x_{j}}\unicode[STIX]{x1D70F}_{ij}+\unicode[STIX]{x1D6FF}_{1i}f_{B}, & \displaystyle\end{eqnarray}$$

(2.2c ) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}t}(\unicode[STIX]{x1D70C}E)+\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}x_{j}}[u_{j}(\unicode[STIX]{x1D70C}E+p)]=\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}x_{j}}(u_{i}\unicode[STIX]{x1D70F}_{ij}-q_{j}), & \displaystyle\end{eqnarray}$$

respectively, where $x_{i}$ ( $x_{1}$ , $x_{2}$ or equivalently, $x$ , $y$ ) are the axial and cross-sectional coordinates, $u_{i}$ are the velocity components in each of those directions and $p$ , $\unicode[STIX]{x1D70C}$ , $T$ and $E$ are the instantaneous pressure, density, temperature and total energy per unit mass, respectively.

Due to the propagation of finite amplitude nonlinear acoustic waves and the periodic nature of the set-up, Gedeon streaming (Gedeon Reference Gedeon and Ross1995) is expected. The mean flow caused by the Gedeon streaming results in the transport of heat away from the hot side of the regenerator into the whole device affecting the mean temperature distribution outside the regenerator. Such thermal leakage is usually mitigated by a secondary cold heat exchanger and a thermal buffer tube to achieve steady-state conditions. Due to the inhomogeneous mean temperature outside the regenerator, such a set-up would exhibit a very large design parameter space. Moreover, the mean flow caused by the streaming would further disperse the acoustic waves, affecting the wave steepening and thus delaying the steady state further. Such effects fall beyond the scope of the present investigation. Thus, acoustic streaming in the current work is purposefully suppressed. To this end, a uniform mean pressure gradient $f_{B}$ in the axial direction is dynamically adjusted to relax the net axial mass flow rate to zero. The relevant expression for $f_{B}$ is:

(2.3) $$\begin{eqnarray}f_{B}=\frac{\unicode[STIX]{x1D6FC}}{\unicode[STIX]{x0394}t}({\dot{m}}-{\dot{m}}_{0}),\end{eqnarray}$$

where $\unicode[STIX]{x1D6FC}=0.3$ is a relaxation coefficient, $\unicode[STIX]{x0394}t$ is the time step, ${\dot{m}}$ is the instantaneous volume-averaged mass flow rate and ${\dot{m}}_{0}=0$ is the target value. The viscous stress tensor, $\unicode[STIX]{x1D70F}_{ij}$ , and heat flux, $q_{i}$ , are formulated based on the Stokes and the Fourier laws as

(2.4a,b ) $$\begin{eqnarray}\unicode[STIX]{x1D70F}_{ij}=2\unicode[STIX]{x1D707}\left[\unicode[STIX]{x1D61A}_{ij}-\frac{1}{3}\frac{\unicode[STIX]{x2202}u_{k}}{\unicode[STIX]{x2202}x_{k}}\unicode[STIX]{x1D6FF}_{ij}\right],\quad q_{j}=-\frac{\unicode[STIX]{x1D707}C_{p}}{Pr}\frac{\unicode[STIX]{x2202}T}{\unicode[STIX]{x2202}x_{j}},\end{eqnarray}$$

respectively, where $\unicode[STIX]{x1D61A}_{ij}$ is the strain-rate tensor, given by $\unicode[STIX]{x1D61A}_{ij}=(\unicode[STIX]{x2202}u_{j}/\unicode[STIX]{x2202}x_{i}+\unicode[STIX]{x2202}u_{i}/\unicode[STIX]{x2202}x_{j})/2$ , $Pr$ is the Prandtl number and $C_{p}$ is the specific heat capacity at constant pressure. The working fluid is air, assumed to be calorically and thermally perfect. The dynamic viscosity, $\unicode[STIX]{x1D707}(T)$ , is varied with the temperature according to the Sutherland’s law, $\unicode[STIX]{x1D707}(T)=\unicode[STIX]{x1D707}_{ref}(T/T_{S,ref})^{1.5}(T_{S,ref}+S)/(T+S)$ where $S=120~\text{K}$ is the Sutherland constant, $T_{S,ref}=291.15~\text{K}$ and $\unicode[STIX]{x1D707}_{ref}=1.827\times 10^{-5}~\text{kg}~\text{m}^{-1}~\text{s}^{-1}$ . The values of the other unspecified fluid parameters, valid for air, are $Pr=0.72$ , the ratio of isobaric to isochoric specific heat capacities $\unicode[STIX]{x1D6FE}=1.4$ , the reference density $\unicode[STIX]{x1D70C}_{ref}=1.176~\text{kg}~\text{m}^{-3}$ , the pressure $p_{ref}=101\,325~\text{Pa}$ , the temperature $T_{ref}=300~\text{K}$ and the gas constant $R=p_{ref}/(\unicode[STIX]{x1D70C}_{ref}\,T_{ref})$ .

Figure 3.

Illustration of selected portions of the unstructured mesh near the duct interfaces for Grid A with total number of elements $N_{el}=2004$ .

Table 2.

Discretization order used per element, $p$ , total number of elements, $N_{el}$ , degrees of freedom, $N_{dof}=p^{2}N_{el}$ , and minimum and maximum linear element size $\unicode[STIX]{x1D6E5}_{h_{(min)}}$ , $\unicode[STIX]{x1D6E5}_{h_{(max)}}$ (cf. (2.9)). Reported values are the ones used at the limit cycle. In all cases, simulations are carried out with $p=3$ and Grid A (figure 3) throughout the transient.

The fully compressible high-fidelity Navier–Stokes calculations have been carried out with the discontinuous finite element sd3DvisP solver, an MPI parallelized Fortran 90 code employing the spectral difference local spatial reconstruction for hexahedral elements on unstructured grids (Kopriva Reference Kopriva1996; Kopriva & Kolias Reference Kopriva and Kolias1996; Sun et al. Reference Sun, Wang and Liu2007; Jameson Reference Jameson2010). The solver reconstructs the local solution inside each element as the tensor product of polynomials up to the user-specified order $p=N_{p}-1$ , where $N_{p}$ is the number of solution points per dimension inside the element. Inter-element discontinuities in the solution are handled utilizing Roe’s (Reference Roe1981) flux with the entropy correction by Harten & Hyman (Reference Harten and Hyman1983). The numerical dissipation at the element interfaces scales as $\unicode[STIX]{x1D6E5}_{h}^{N_{p}+1}$ where $\unicode[STIX]{x1D6E5}_{h}$ is the characteristic length scale of neighbouring elements (Jameson & Lodato Reference Jameson and Lodato2014; Chapelier, Lodato & Jameson Reference Chapelier, Lodato and Jameson2016). Sub-cell shock capturing is enabled through a Laplacian artificial diffusion term applied in regions of steep gradients, which are detected by means of a modal sensor based on a Legendre polynomial expansion (Persson & Peraire Reference Persson and Peraire2006; Lodato, Vervisch & Clavin Reference Lodato, Vervisch and Clavin2016). The time integration is carried out explicitly with a third-order Runge–Kutta scheme and discretization order $p=3$ during the transient, whereas a fifth-order Runge–Kutta scheme is adopted at the shock-dominated limit cycle where a systematic grid sensitivity study has been carried out (cf. table 2 and figure 4), as discussed in the following section. The same solver has been used and validated in a wide variety of flow configurations, including turbulent channel flow (Lodato, Castonguay & Jameson Reference Lodato, Castonguay and Jameson2013, Reference Lodato, Castonguay and Jameson2014) and unsteady shock–wavy wall interaction problems (Lodato et al. Reference Lodato, Vervisch and Clavin2016; Lodato, Vervisch & Clavin Reference Lodato, Vervisch and Clavin2017).

2.2.2 Grid sensitivity analysis

A grid sensitivity analysis is carried out at the shock-dominated limit cycle spanning the simulation parameter space outlined in table 2 only for the $T_{H}=450~\text{K}$ temperature setting. The pressure waveform, expressed as a function of the travelling-wave coordinate,

(2.5) $$\begin{eqnarray}\unicode[STIX]{x1D709}=c_{s}t-x,\end{eqnarray}$$

where $c_{s}\simeq a_{0}$ is the (weak) shock speed, and $a_{0}$ is the base speed of sound, monotonically converges to a step function (figure 4 a) as the numerical degrees of freedom, $N_{dof}$ , of the high-order Navier–Stokes calculations are increased. The Mach number of the captured shock waves (figure 4 c) is estimated, for each grid resolution level in table 2, by calculating the total entropy jump, $\unicode[STIX]{x0394}s$ , as the sum of the point-to-point entropy increments:

(2.6) $$\begin{eqnarray}\frac{\unicode[STIX]{x0394}s}{R}=\mathop{\sum }_{i}\frac{\unicode[STIX]{x0394}s_{i}}{R}\approx \frac{\unicode[STIX]{x1D6FE}+1}{12\unicode[STIX]{x1D6FE}^{2}}\mathop{\sum }_{i}\left(\frac{\unicode[STIX]{x0394}p_{i}}{p_{i}}\right)^{3},\end{eqnarray}$$

where $\unicode[STIX]{x0394}p_{i}=p_{i+1}-p_{i}$ , $\unicode[STIX]{x0394}s_{i}=s_{i+1}-s_{i}$ are the pressure and the entropy jumps from point $i$ to $i+1$ in the $\unicode[STIX]{x1D709}$ space (figure 4 b). Thus the Mach number of the captured wave is obtained upon solving

(2.7) $$\begin{eqnarray}\frac{\unicode[STIX]{x0394}s}{R}\approx \frac{2\unicode[STIX]{x1D6FE}}{(\unicode[STIX]{x1D6FE}+1)^{2}}\frac{(M^{2}-1)^{3}}{3}.\end{eqnarray}$$

Evaluating the overall entropy jump, and therefore, the shock Mach number, based on the total pressure jump, would include spurious contributions from isentropic compression waves which do not participate in the coalescence of waves forming the shock. Finally, the shock Mach number $M_{r_{\unicode[STIX]{x1D6E5}}\rightarrow 0}$ at the continuum limit is estimated via a polynomial extrapolation,

(2.8) $$\begin{eqnarray}M(r_{\unicode[STIX]{x1D6E5}})=m_{0}+m_{1}r_{\unicode[STIX]{x1D6E5}}+m_{2}r_{\unicode[STIX]{x1D6E5}}^{2}+m_{3}r_{\unicode[STIX]{x1D6E5}}^{3},\end{eqnarray}$$

where $m_{0}$ , $m_{1}$ , $m_{2}$ and $m_{3}$ are fitting parameters and $r_{\unicode[STIX]{x1D6E5}}$ is the grid refinement ratio,

(2.9a,b ) $$\begin{eqnarray}r_{\unicode[STIX]{x1D6E5}}=\overline{\unicode[STIX]{x1D6E5}}_{h}/\overline{\unicode[STIX]{x1D6E5}}_{h}|_{p=7}^{B},\quad \overline{\unicode[STIX]{x1D6E5}}_{h}=\sqrt{\unicode[STIX]{x1D6FA}/N_{dof}},\end{eqnarray}$$

in which $\unicode[STIX]{x1D6FA}$ denotes total area of the computational domain.

Figure 4.

Grid sensitivity analysis at the shock-dominated limit cycle for $T_{H}=450~\text{K}$ . Pressure versus travelling-wave coordinate $\unicode[STIX]{x1D709}$ for increasing numerical degrees of freedom (table 2) from time series extracted at $x=-2~\text{mm}$ (figure 1) (a), point-to-point dimensionless entropy increments (b), and shock Mach number $M$ versus grid refinement ratio $r_{\unicode[STIX]{x1D6E5}}$ (c). (c): (- - - -), polynomial extrapolation (cf. (2.8)); (○), shock Mach number; (▫), Mach number in the continuum limit $r_{\unicode[STIX]{x1D6E5}}\rightarrow 0$ .

3.1 Spectral energy density

Any acoustic dynamical model (linear or nonlinear) can be written as

(3.1) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}\boldsymbol{X}}{\unicode[STIX]{x2202}t}=\boldsymbol{f}(\boldsymbol{X}),\quad \text{with }\boldsymbol{X}=\left(\frac{u^{\prime }}{a_{0}},\frac{p^{\prime }}{\unicode[STIX]{x1D70C}_{0}a_{0}^{2}}\right)^{\text{T}},\end{eqnarray}$$

where $\boldsymbol{X}$ is the state vector containing the dimensionless perturbation variables. Throughout, we define the following squared $L_{2}$ norm as the perturbation energy density (Naugol’Nykh & Rybak Reference Naugol’Nykh and Rybak1975):

(3.2) $$\begin{eqnarray}E=\frac{1}{2}\unicode[STIX]{x1D70C}_{0}a_{0}^{2}\boldsymbol{X}^{\text{T}}\boldsymbol{X}=\frac{1}{2}\unicode[STIX]{x1D70C}_{0}{u^{\prime }}^{2}+\frac{{p^{\prime }}^{2}}{2\unicode[STIX]{x1D70C}_{0}a_{0}^{2}}.\end{eqnarray}$$

In the nonlinear growth and limit cycle regimes, velocity and pressure fluctuations, $u^{\prime }$ and $p^{\prime }$ , respectively, are composed of higher harmonics of the linearly unstable mode. Substituting velocity and pressure fluctuations, expressed as complex Fourier expansions,

(3.3a,b ) $$\begin{eqnarray}u^{\prime }(x,t)=\mathop{\sum }_{\substack{ k=-\infty \\ k\neq 0}}^{+\infty }\hat{u} _{k}(x,\unicode[STIX]{x1D716}t)\text{e}^{\text{i}(k/2)\unicode[STIX]{x1D714}_{0}t},\quad p^{\prime }(x,t)=\mathop{\sum }_{\substack{ k=-\infty \\ k\neq 0}}^{+\infty }\hat{p}_{k}(x,\unicode[STIX]{x1D716}t)\text{e}^{\text{i}(k/2)\unicode[STIX]{x1D714}_{0}t},\end{eqnarray}$$

(3.4a,b ) $$\begin{eqnarray}\hat{u} _{-k}=\hat{u} _{k}^{\ast },\quad \hat{p}_{-k}=\hat{p}_{k}^{\ast },\end{eqnarray}$$

into (3.2) and cycle averaging yields:

(3.5a,b ) $$\begin{eqnarray}\overline{E}=2\mathop{\sum }_{k=1}^{\infty }E_{k},\quad E_{k}=\frac{1}{2}\unicode[STIX]{x1D70C}_{0}|\hat{u} _{k}|^{2}+\frac{|\hat{p}_{k}|^{2}}{2\unicode[STIX]{x1D70C}_{0}a_{0}^{2}},\end{eqnarray}$$

where $\unicode[STIX]{x1D716}\sim \unicode[STIX]{x1D6FC}/\unicode[STIX]{x1D714}\ll 1$ is the smallness parameter such that $t$ and $\unicode[STIX]{x1D716}t$ correspond to fast and slow time scales, respectively, $\overline{(\cdot )}$ denotes the cycle averaging operator defined as

(3.6) $$\begin{eqnarray}\overline{(\cdot )}=\frac{1}{T_{0}}\int _{\unicode[STIX]{x1D716}t}^{\unicode[STIX]{x1D716}t+T_{0}}(\cdot )\,\text{d}t,\quad T_{0}=\frac{2\unicode[STIX]{x03C0}}{\unicode[STIX]{x1D714}_{0}},\end{eqnarray}$$

and $(\cdot )^{\ast }$ denotes the complex conjugate. Here $\unicode[STIX]{x1D714}_{0}$ is the angular frequency of the unstable second harmonic, and $E_{k}$ is the spectral energy density of the $k$ th mode. The pressure and velocity amplitudes of the $k$ th harmonic ( $|\hat{p}_{k}|$ and $|\hat{u} _{k}|$ , respectively) are functions of the $x$ coordinate and the slow time $\unicode[STIX]{x1D716}t$ , and are extracted via a short time-windowed Fourier transform (over 8 cycles of time period $T_{0}$ ) of the time series shown in figure 5. In the nonlinear growth regime, the energy cascades from the unstable second mode ( $k=2$ ) into its overtones only ( $k=4,6,8,\ldots$ ) with no energy content in the odd-numbered harmonics.

3.2 Regime classification

Based on a scale-by-scale analysis of the growth of spectral energy density (figure 6), the aforementioned three regimes of thermoacoustic wave amplification are identified as follows.

1. (i) Modal growth. Only the thermoacoustically unstable mode amplifies exponentially and all the other modes excited by the initial perturbation field decay (see § 4). Higher harmonics begin to grow after $\unicode[STIX]{x1D6FC}t\approx 3.2$ (figure 6), setting the end of a purely harmonic growth. In the modal growth regime, the system is well approximated by the linear system of equations.

2. (ii) Hierarchical spectral broadening. Energy cascades down to higher harmonics hierarchically: the $k$ th harmonic ( $k>2$ ) grows at a rate equal to $k/2$ -times the modal growth rate of the second harmonic (see §§ 5 and 6), $\unicode[STIX]{x1D6FC}$ , that is

   (3.7) $$\begin{eqnarray}\unicode[STIX]{x1D6FC}_{k}=\unicode[STIX]{x1D6FC}k/2,\quad k\in \{4,6,8,\ldots \}.\end{eqnarray}$$
   The saturation of the spectral energy density $E_{k}$ occurs at $\unicode[STIX]{x1D6FC}t\approx 10$ followed by the formation of resonating shock waves at limit cycle.

3. (iii) Shock-dominated limit cycle. In this regime ( $\unicode[STIX]{x1D6FC}t>10$ ), the continued injection of energy in the second-mode harmonic is balanced by the cascade of the spectral energy density into the overtones of the second mode and terminates by viscous dissipation at very high overtones ( $k\approx 300$ ). The maximum number of overtones generated is a function of the acoustic phasing of the unstable mode. Moreover, at the limit cycle, the spectral energy density scales with instability growth rate $\unicode[STIX]{x1D6FC}$ approximately as $\unicode[STIX]{x1D6FC}^{3}$ (§ 7.2).

4.1 Linear thermoviscous quasi-planar wave equations

The time-domain linear thermoviscous governing equations for a two-dimensional perturbation are:

(4.1a ) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}^{\prime }}{\unicode[STIX]{x2202}t}+\unicode[STIX]{x1D70C}_{0}\frac{\unicode[STIX]{x2202}u^{\prime }}{\unicode[STIX]{x2202}x}+u^{\prime }\frac{\text{d}\unicode[STIX]{x1D70C}_{0}}{\text{d}x}+\unicode[STIX]{x1D70C}_{0}\frac{\unicode[STIX]{x2202}v^{\prime }}{\unicode[STIX]{x2202}y}=0, & \displaystyle\end{eqnarray}$$

(4.1b ) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}u^{\prime }}{\unicode[STIX]{x2202}t}+\frac{1}{\unicode[STIX]{x1D70C}_{0}}\frac{\unicode[STIX]{x2202}p^{\prime }}{\unicode[STIX]{x2202}x}-\unicode[STIX]{x1D708}_{0}\frac{\unicode[STIX]{x2202}^{2}u^{\prime }}{\unicode[STIX]{x2202}y^{2}}=0, & \displaystyle\end{eqnarray}$$

(4.1c ) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}s^{\prime }}{\unicode[STIX]{x2202}t}+u^{\prime }\frac{\text{d}s_{0}}{\text{d}x}-\frac{k_{0}}{\unicode[STIX]{x1D70C}_{0}T_{0}}\frac{\unicode[STIX]{x2202}^{2}T^{\prime }}{\unicode[STIX]{x2202}y^{2}}=0, & \displaystyle\end{eqnarray}$$

$(\cdot )^{\prime }$

$0$

$y$

(4.2a ) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}p^{\prime }}{\unicode[STIX]{x2202}t}+\frac{\unicode[STIX]{x1D70C}_{0}a_{0}^{2}}{h}\frac{\unicode[STIX]{x2202}U^{\prime }}{\unicode[STIX]{x2202}x}=\frac{\unicode[STIX]{x1D70C}_{0}a_{0}^{2}}{h}\frac{q^{\prime }}{\unicode[STIX]{x1D70C}_{0}C_{p}T_{0}}, & \displaystyle\end{eqnarray}$$

(4.2b ) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{1}{h}\frac{\unicode[STIX]{x2202}U^{\prime }}{\unicode[STIX]{x2202}t}+\frac{1}{\unicode[STIX]{x1D70C}_{0}}\frac{\unicode[STIX]{x2202}p^{\prime }}{\unicode[STIX]{x2202}x}=\frac{1}{h}\frac{\unicode[STIX]{x1D70F}_{w}^{\prime }}{\unicode[STIX]{x1D70C}_{0}}, & \displaystyle\end{eqnarray}$$

$h$

$U^{\prime }$

(4.3) $$\begin{eqnarray}U^{\prime }=\int _{-h/2}^{+h/2}u^{\prime }(x,y,t)\,\text{d}y,\end{eqnarray}$$

whereas $q^{\prime }$ and $\unicode[STIX]{x1D70F}_{w}^{\prime }$ are the wall-heat flux and the wall shear, respectively:

(4.4a,b ) $$\begin{eqnarray}q^{\prime }=2k_{0}\left.\frac{\unicode[STIX]{x2202}T^{\prime }}{\unicode[STIX]{x2202}y}\right|_{y=+h/2},\quad \unicode[STIX]{x1D70F}_{w}^{\prime }=2\unicode[STIX]{x1D707}_{0}\left.\frac{\unicode[STIX]{x2202}u^{\prime }}{\unicode[STIX]{x2202}y}\right|_{y=+h/2}.\end{eqnarray}$$

In sections $a,c$ and $d$ , linear wave propagation is assumed to be inviscid and adiabatic, hence isentropic, resulting in $q^{\prime }=0$ and $\unicode[STIX]{x1D70F}_{w}^{\prime }=0$ . Applying the normal mode assumption to (4.2), namely,

(4.5a,b ) $$\begin{eqnarray}p^{\prime }(x,t)=\hat{p}(x)\text{e}^{\unicode[STIX]{x1D70E}t},\quad U^{\prime }(x,t)=\hat{U} (x)\text{e}^{\unicode[STIX]{x1D70E}t},\end{eqnarray}$$

where $\unicode[STIX]{x1D70E}=\unicode[STIX]{x1D6FC}+\text{i}\unicode[STIX]{x1D714}$ is the complex eigenvalue of the system with growth rate $\unicode[STIX]{x1D6FC}$ and angular frequency $\unicode[STIX]{x1D714}$ , leads to the thermoviscous set of quasi-planar wave equations in the frequency domain (Lin et al. Reference Lin, Scalo and Hesselink2016):

(4.6a ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D70E}\hat{p}=\frac{\unicode[STIX]{x1D70C}_{0}a_{0}^{2}}{h_{b}}\left[\frac{1}{1+(\unicode[STIX]{x1D6FE}-1)f_{k}}\left(\frac{\unicode[STIX]{x1D6E9}(f_{k}-f_{\unicode[STIX]{x1D708}})}{(1-f_{\unicode[STIX]{x1D708}})(1-Pr)}-\frac{\text{d}}{\text{d}x}\right)\right]\hat{U} ,\quad \unicode[STIX]{x1D6E9}=\frac{1}{T_{0}}\frac{\text{d}T_{0}}{\text{d}x}, & \displaystyle\end{eqnarray}$$

(4.6b ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D70E}\hat{U} =-\frac{h_{b}}{\unicode[STIX]{x1D70C}_{0}}(1-f_{\unicode[STIX]{x1D708}})\frac{\text{d}\hat{p}}{\text{d}x}, & \displaystyle\end{eqnarray}$$

$f_{\unicode[STIX]{x1D708}}$

$f_{k}$

(4.7a,b ) $$\begin{eqnarray}f_{\unicode[STIX]{x1D708}}=\frac{\tanh (\unicode[STIX]{x1D702}h_{b}/2)}{\unicode[STIX]{x1D702}h_{b}/2},\quad f_{k}=\frac{\tanh (\unicode[STIX]{x1D702}h_{b}\sqrt{Pr}/2)}{\unicode[STIX]{x1D702}h_{b}\sqrt{Pr}/2},\quad \text{with }\unicode[STIX]{x1D702}=\sqrt{\text{i}\unicode[STIX]{x1D714}/\unicode[STIX]{x1D708}_{0}}.\end{eqnarray}$$

4.2 Acoustic energy budgets for quasi-planar wave perturbations

Multiplying (4.2a ) by $p^{\prime }/(\unicode[STIX]{x1D70C}_{0}a_{0}^{2})$ and (4.2b ) by $\unicode[STIX]{x1D70C}_{0}U^{\prime }/h_{b}$ and adding them, yields the conservation equation

(4.8) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}E}{\unicode[STIX]{x2202}t}+\frac{\unicode[STIX]{x2202}{\mathcal{I}}}{\unicode[STIX]{x2202}x}={\mathcal{P}}-{\mathcal{D}},\end{eqnarray}$$

for the one-dimensional acoustic energy density

(4.9) $$\begin{eqnarray}E=\frac{1}{2}\frac{{p^{\prime }}^{2}}{\unicode[STIX]{x1D70C}_{0}a_{0}^{2}}+\frac{1}{2}\unicode[STIX]{x1D70C}_{0}\left(\frac{U^{\prime }}{h_{b}}\right)^{2},\end{eqnarray}$$

consistent with the definition (3.2). The instantaneous acoustic flux ${\mathcal{I}}$ and the net energy production ${\mathcal{P}}-{\mathcal{D}}$ therein are given by

(4.10a,b ) $$\begin{eqnarray}{\mathcal{I}}=\frac{p^{\prime }U^{\prime }}{h_{b}},\quad {\mathcal{P}}-{\mathcal{D}}=\frac{p^{\prime }q^{\prime }}{h_{b}\unicode[STIX]{x1D70C}_{0}C_{p}T_{0}}+\frac{\unicode[STIX]{x1D70F}_{w}^{\prime }U^{\prime }}{h_{b}^{2}}.\end{eqnarray}$$

Averaging (4.8) over one acoustic cycle and integrating axially over the periodic domain, $L$ , yields:

(4.11) $$\begin{eqnarray}\frac{\text{d}}{\text{d}(\unicode[STIX]{x1D716}t)}\int _{L}\overline{E}(x,\unicode[STIX]{x1D716}t)\,\text{d}x=\int _{L}(\overline{{\mathcal{P}}}-\overline{{\mathcal{D}}})\,\text{d}x={\mathcal{R}},\end{eqnarray}$$

where ${\mathcal{R}}$ is the Rayleigh index and $\unicode[STIX]{x1D716}t$ is the slow time scale (cf. (3.5)). Relation (4.11) allows us to unambiguously identify the onset of an instability via the criterion  ${\mathcal{R}}>0$ . This expression also accounts for wall shear and wall-heat flux losses outside the regenerator (if present), which attenuate the thermoacoustic instability. Such thermoviscous losses are captured in the heat flux $q^{\prime }$ and shear stress $\unicode[STIX]{x1D70F}_{w}^{\prime }$ terms (cf. (4.10)) in the respective duct sections. Utilizing the frequency-domain linear equations (4.6), the wall-heat flux $\hat{q}$ and wall shear $\hat{\unicode[STIX]{x1D70F}}_{w}$ in the frequency domain are given by

(4.12a ) $$\begin{eqnarray}\displaystyle & \displaystyle \hat{\unicode[STIX]{x1D70F}}_{w}=h_{b}\frac{\unicode[STIX]{x2202}\hat{p}}{\unicode[STIX]{x2202}x}f_{\unicode[STIX]{x1D708}}, & \displaystyle\end{eqnarray}$$

(4.12b ) $$\begin{eqnarray}\displaystyle & \displaystyle \hat{q}=h_{b}(\text{i}\unicode[STIX]{x1D714})C_{p}T_{0}\left(\frac{\unicode[STIX]{x1D6E9}}{(1-Pr)\unicode[STIX]{x1D714}^{2}}\frac{\unicode[STIX]{x2202}\hat{p}}{\unicode[STIX]{x2202}x}(f_{k}-f_{\unicode[STIX]{x1D708}})-\frac{\unicode[STIX]{x1D6FE}-1}{a_{0}^{2}}\hat{p}f_{k}\right). & \displaystyle\end{eqnarray}$$

Figure 7.

Variation of thermoviscous functionals affecting the cycle-averaged thermoacoustic production $\overline{{\mathcal{P}}}$ (a), and dissipation $\overline{{\mathcal{D}}}$ (b,c) of acoustic energy density versus the ratio of the regenerator half-width $h_{b}/2$ to the Stokes boundary layer thickness $\unicode[STIX]{x1D6FF}_{k}$ (cf. (4.13)). (a) (——), $\unicode[STIX]{x1D6F7}_{{\mathcal{P}}}^{T}$ ; (- - - -),  $-\unicode[STIX]{x1D6F7}_{{\mathcal{P}}}^{S}$ .

Combining (4.10), (4.12), and (4.6b ), an analytical expression is obtained for the cycle-averaged production $\overline{{\mathcal{P}}}$ and dissipation $\overline{{\mathcal{D}}}$ of the acoustic energy density (4.9):

(4.13a ) $$\begin{eqnarray}\displaystyle & \displaystyle \overline{{\mathcal{P}}}=\frac{\unicode[STIX]{x1D6E9}}{2(1-Pr)h_{b}}[\unicode[STIX]{x1D6F7}_{{\mathcal{P}}}^{T}\,\text{Re}(\hat{p}^{\ast }\hat{U} )-\unicode[STIX]{x1D6F7}_{{\mathcal{P}}}^{S}\,\text{Im}(\hat{p}^{\ast }\hat{U} )], & \displaystyle\end{eqnarray}$$

(4.13b ) $$\begin{eqnarray}\displaystyle & \displaystyle \overline{{\mathcal{D}}}=\text{Re}(\text{i}f_{k})\frac{\unicode[STIX]{x1D714}(\unicode[STIX]{x1D6FE}-1)}{2\unicode[STIX]{x1D70C}_{0}a_{0}^{2}}|\hat{p}|^{2}-\text{Im}\left(\frac{f_{\unicode[STIX]{x1D708}}}{1-f_{\unicode[STIX]{x1D708}}}\right)\frac{\unicode[STIX]{x1D70C}_{0}\unicode[STIX]{x1D714}}{2h_{b}^{2}}|\hat{U} |^{2}. & \displaystyle\end{eqnarray}$$

$\unicode[STIX]{x1D6F7}_{{\mathcal{P}}}^{T}$

$-\unicode[STIX]{x1D6F7}_{{\mathcal{P}}}^{S}$

$\text{Re}(\hat{p}^{\ast }\hat{U} )$

$\text{Im}(\hat{p}^{\ast }\,\hat{U} )$

(4.14a,b ) $$\begin{eqnarray}\unicode[STIX]{x1D6F7}_{{\mathcal{P}}}^{T}=\text{Re}\left(\frac{f_{k}-f_{\unicode[STIX]{x1D708}}}{1-f_{\unicode[STIX]{x1D708}}}\right),\quad \unicode[STIX]{x1D6F7}_{{\mathcal{P}}}^{S}=\text{Im}\left(\frac{f_{k}-f_{\unicode[STIX]{x1D708}}}{1-f_{\unicode[STIX]{x1D708}}}\right).\end{eqnarray}$$

Figure 8.

Normalized thermoacoustic production $\unicode[STIX]{x1D6F7}_{{\mathcal{P}}}^{T}\cos (\unicode[STIX]{x1D713}_{Up})-\unicode[STIX]{x1D6F7}_{{\mathcal{P}}}^{S}\sin (\unicode[STIX]{x1D713}_{Up})$  (cf. (4.13a )) versus the phase angle difference between $\hat{p}$ and $\hat{U}$ , $\unicode[STIX]{x1D713}_{Up}=\angle \hat{U} -\angle \hat{p}$ , for increasing values of $h_{b}/2\unicode[STIX]{x1D6FF}_{k}=0.5$ , 1.0, 1.5, 2.0 (——); (- - - -), optimum phasing maximizing the thermoacoustic production for continually varying $h_{b}/2\unicode[STIX]{x1D6FF}_{k}$ .

For $h_{b}/2\unicode[STIX]{x1D6FF}_{k}\leqslant 1.13$ , $|\unicode[STIX]{x1D6F7}_{{\mathcal{P}}}^{T}|>|\unicode[STIX]{x1D6F7}_{{\mathcal{P}}}^{S}|$ , which implies that the regenerator half-width $h_{b}$ must remain comparable to or smaller than the Stokes boundary layer thickness $\unicode[STIX]{x1D6FF}_{k}$ , to achieve higher thermoacoustic amplification of travelling waves ( $\unicode[STIX]{x1D713}_{Up}\simeq 0^{\circ }$ ) (figures 7 a and 8). However, to maximize thermoacoustic energy production for standing waves ( $\unicode[STIX]{x1D713}_{Up}\simeq \pm 90^{\circ }$ ), a larger regenerator half-width ( $h_{b}/2\unicode[STIX]{x1D6FF}_{k}>1.13$ ) is required. While production alone for a purely travelling wave ( $\unicode[STIX]{x1D713}_{Up}=0$ ) is maximized in the limit $h_{b}/2\unicode[STIX]{x1D6FF}_{k}\rightarrow 0$ (figure 8), for fixed temperature settings, dissipation also diverges (figure 7 c). Therefore, pure travelling-wave phasing, if at all achieved, always results in smaller net production of acoustic energy density compared to an optimal combination of standing and travelling waves. Moreover, varying the temperature inside the regenerator results in a local variation of the ratio $h_{b}/2\unicode[STIX]{x1D6FF}_{k}$ which, in turn, causes the optimum phase to vary along the regenerator (figure 8). For the temperature settings considered here (table 1), the optimum phasing angle, averaged over the regenerator length changes from $43.17^{\circ }$ to $36.41^{\circ }$ as the temperature $T_{H}$ is increased.

4.3 Effects of temperature gradient and geometry on growth rates

The variation of the instability growth rates with the temperature and the geometry is further analysed utilizing the acoustic energy budget formulation developed in the previous section. To this end, the production, dissipation and Rayleigh index (normalized by the pressure amplitude) are plotted in the convenient dimensionless forms

(4.15a-c ) $$\begin{eqnarray}\overline{{\mathcal{P}}}_{\ast }=\frac{\overline{{\mathcal{P}}}\unicode[STIX]{x1D70F}}{\unicode[STIX]{x1D70C}_{0}a_{0}^{2}},\quad \overline{{\mathcal{D}}}_{\ast }=\frac{\overline{{\mathcal{D}}}\unicode[STIX]{x1D70F}}{\unicode[STIX]{x1D70C}_{0}a_{0}^{2}},\quad {\mathcal{R}}_{\ast }=\frac{{\mathcal{R}}\unicode[STIX]{x1D70F}}{h_{b}\unicode[STIX]{x1D70C}_{0}a_{0}^{2}}.\end{eqnarray}$$

Cycle-averaged production of acoustic energy density due to travelling-wave and standing-wave components given by

(4.16a ) $$\begin{eqnarray}\displaystyle & \displaystyle {\mathcal{R}}_{T\ast }=\frac{\unicode[STIX]{x1D70F}}{2(1-Pr)\unicode[STIX]{x1D70C}_{0}a_{0}^{2}h_{b}^{2}}\int _{b}\unicode[STIX]{x1D6E9}\,\text{Re}(\hat{p}^{\ast }\hat{U} )\unicode[STIX]{x1D6F7}_{{\mathcal{P}}}^{T}\,\text{d}x, & \displaystyle\end{eqnarray}$$

(4.16b ) $$\begin{eqnarray}\displaystyle & \displaystyle {\mathcal{R}}_{S\ast }=-\frac{\unicode[STIX]{x1D70F}}{2(1-Pr)\unicode[STIX]{x1D70C}_{0}a_{0}^{2}h_{b}^{2}}\int _{b}\unicode[STIX]{x1D6E9}\,\text{Im}(\hat{p}^{\ast }\hat{U} )\unicode[STIX]{x1D6F7}_{{\mathcal{P}}}^{S}\,\text{d}x, & \displaystyle\end{eqnarray}$$

(4.17) $$\begin{eqnarray}\unicode[STIX]{x1D70F}=h_{b}^{2}/\unicode[STIX]{x1D708}_{0}\end{eqnarray}$$

is a reference viscous time scale in the regenerator with $\unicode[STIX]{x1D708}_{0}$ evaluated at $T_{C}=300~\text{K}$ .

Figure 9.

Dimensionless thermoacoustic growth rate, $\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D70F}$ (cf. (4.17)), versus the natural logarithm of the temperature ratio $T_{H}/T_{C}$ for $T_{C}=300~\text{K}$ (a) and resonator area ratio $h_{d}/h_{c}$ (b). (a) (——), Linear stability analysis; (- - - -), logarithmic estimate of $\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D70F}$  (4.19) fitted using values at $T_{H}=450~\text{K}$ and $T_{H}=550~\text{K}$ ; (○), Navier–Stokes simulations.

Figure 10.

Dimensionless Rayleigh index ${\mathcal{R}}_{\ast }$ versus the natural logarithm of the temperature ratio $T_{H}/T_{C}$ for $T_{C}=300~\text{K}$ (a) and the resonator area ratio $h_{d}/h_{c}$ (b). (b) (——), ${\mathcal{R}}_{\ast }$ ; (- - - -), ${\mathcal{R}}_{T\ast }$ ; ( $\cdots \cdots$ ),  ${\mathcal{R}}_{S\ast }$ .

Increasing the hot-side temperature $T_{H}$ , the thermoacoustic production $\overline{{\mathcal{P}}}_{\ast }$ increases monotonically, approximately as

(4.18) $$\begin{eqnarray}\overline{{\mathcal{P}}}_{\ast }\sim \unicode[STIX]{x1D6E9}=\frac{\text{d}}{\text{d}x}\ln [T_{0}(x)],\end{eqnarray}$$

and more rapidly than the dissipation $\overline{{\mathcal{D}}}_{\ast }$ , yielding positive values of the Rayleigh index ${\mathcal{R}}_{\ast }$ (figure 10 a) for $T_{H}/T_{C}>1.29$ (figure 9 a). The Rayleigh index can thus be used to quantify the thermoacoustic growth rate (figures 9 a, 10 a) as

(4.19) $$\begin{eqnarray}{\mathcal{R}}_{\ast }\sim \unicode[STIX]{x1D6FC}\unicode[STIX]{x1D70F}\approx A\ln (T_{H}/T_{C})-B,\end{eqnarray}$$

where $A$ and $B$ are geometry dependent fitting coefficients. With increasing resonator area ratio $h_{d}/h_{c}$ , keeping $h_{c}$ fixed, the frequency decreases monotonically by approximately 4 % in the range of $h_{d}/h_{c}$ considered. However, the growth rates vary non-monotonically, reaching a local maximum at $h_{d}/h_{c}\approx 1.28$ (figure 9 b). Also, the net production of acoustic energy density due to the travelling-wave component ${\mathcal{R}}_{T\ast }$ (figure 10 b) peaks at $h_{d}/h_{c}\approx 1.28$ . Moreover, the high mechanical impedance of the section $d$ ( $\unicode[STIX]{x1D70C}_{0}a_{0}h_{d}$ ) for high values of $h_{d}$ increases the standing-wave component of the acoustic power ${\mathcal{R}}_{S\ast }$ and decreases the travelling-wave component ${\mathcal{R}}_{T\ast }$ . The variation of the transverse geometrical parameters does not significantly alter the frequency. Therefore, changes in the ratio $(f_{k}-f_{\unicode[STIX]{x1D708}})/(1-f_{\unicode[STIX]{x1D708}})$ are also negligible and do not influence the growth rates significantly (cf. (4.14)).

5.1 Governing equations for nonlinear thermoviscous perturbations

The nonlinear governing equations, correct up to second order, for a two-dimensional perturbation read:

(5.1a ) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}^{\prime }}{\unicode[STIX]{x2202}t}+\unicode[STIX]{x1D70C}_{0}\frac{\unicode[STIX]{x2202}u^{\prime }}{\unicode[STIX]{x2202}x}+u^{\prime }\frac{\text{d}\unicode[STIX]{x1D70C}_{0}}{\text{d}x}+\unicode[STIX]{x1D70C}_{0}\frac{\unicode[STIX]{x2202}v^{\prime }}{\unicode[STIX]{x2202}y}=\left[-\unicode[STIX]{x1D70C}^{\prime }\frac{\unicode[STIX]{x2202}u^{\prime }}{\unicode[STIX]{x2202}x}-u^{\prime }\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}^{\prime }}{\unicode[STIX]{x2202}x}\right], & \displaystyle\end{eqnarray}$$

(5.1b ) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}u^{\prime }}{\unicode[STIX]{x2202}t}+\frac{1}{\unicode[STIX]{x1D70C}_{0}}\frac{\unicode[STIX]{x2202}p^{\prime }}{\unicode[STIX]{x2202}x}-\unicode[STIX]{x1D708}_{0}\frac{\unicode[STIX]{x2202}^{2}u^{\prime }}{\unicode[STIX]{x2202}y^{2}}-\frac{1}{\unicode[STIX]{x1D70C}_{0}}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}x}\left[\unicode[STIX]{x1D707}_{0}\left(\unicode[STIX]{x1D709}_{B}+\frac{4}{3}\right)\frac{\unicode[STIX]{x2202}u^{\prime }}{\unicode[STIX]{x2202}x}\right]=\left[-\frac{\unicode[STIX]{x1D70C}^{\prime }}{\unicode[STIX]{x1D70C}_{0}}\frac{\unicode[STIX]{x2202}u^{\prime }}{\unicode[STIX]{x2202}t}-\frac{1}{2}\frac{\unicode[STIX]{x2202}{u^{\prime }}^{2}}{\unicode[STIX]{x2202}x}\right],\qquad & \displaystyle\end{eqnarray}$$

(5.1c ) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}s^{\prime }}{\unicode[STIX]{x2202}t}+u^{\prime }\frac{\text{d}s_{0}}{\text{d}x}-\frac{Rk_{0}}{p_{0}}\frac{\unicode[STIX]{x2202}^{2}T^{\prime }}{\unicode[STIX]{x2202}y^{2}}-\frac{R}{p_{0}}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}x}\left(k_{0}\frac{\unicode[STIX]{x2202}T^{\prime }}{\unicode[STIX]{x2202}x}\right)=\left[-\frac{p^{\prime }}{p_{0}}\left(\frac{\unicode[STIX]{x2202}s^{\prime }}{\unicode[STIX]{x2202}t}+u^{\prime }\frac{\text{d}s_{0}}{\text{d}x}\right)-u^{\prime }\frac{\unicode[STIX]{x2202}s^{\prime }}{\unicode[STIX]{x2202}x}\right],\qquad & \displaystyle\end{eqnarray}$$

$\left(\cdot \right)^{\prime }$

$0$

$\unicode[STIX]{x1D709}_{B}=2/3$

$\unicode[STIX]{x1D707}_{B}$

$\unicode[STIX]{x1D707}$

$y$

$\unicode[STIX]{x1D707}$

$k$

We seek to collapse (5.1a–c ) to obtain a set of equations similar to (3.1). To this end, the following quadratic thermodynamic constitutive equation relating the density fluctuations $\unicode[STIX]{x1D70C}^{\prime }$ with the pressure and entropy fluctuations ( $p^{\prime }$ and $s^{\prime }$ , respectively) is considered:

(5.2) $$\begin{eqnarray}\unicode[STIX]{x1D70C}^{\prime }=\unicode[STIX]{x1D6FC}_{s}p^{\prime }+\unicode[STIX]{x1D6FC}_{p}s^{\prime }+{\textstyle \frac{1}{2}}(\unicode[STIX]{x1D6FD}_{s}{p^{\prime }}^{2}+\unicode[STIX]{x1D6FD}_{p}{s^{\prime }}^{2}+2\unicode[STIX]{x1D6FD}_{sp}s^{\prime }p^{\prime }),\end{eqnarray}$$

where the thermodynamic coefficients $\unicode[STIX]{x1D6FC}$ and $\unicode[STIX]{x1D6FD}$ are given by

(5.3a,b ) $$\begin{eqnarray}\unicode[STIX]{x1D6FC}_{s}=\left(\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}}{\unicode[STIX]{x2202}p}\right)_{s}=\frac{1}{a_{0}^{2}},\quad \unicode[STIX]{x1D6FC}_{p}=\left(\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}}{\unicode[STIX]{x2202}s}\right)_{p}=-\frac{\unicode[STIX]{x1D70C}_{0}}{C_{p}},\end{eqnarray}$$

(5.3c,d ) $$\begin{eqnarray}\unicode[STIX]{x1D6FD}_{s}=\left(\frac{\unicode[STIX]{x2202}^{2}\unicode[STIX]{x1D70C}}{\unicode[STIX]{x2202}p^{2}}\right)_{s}=-\frac{\unicode[STIX]{x1D6FE}-1}{\unicode[STIX]{x1D70C}_{0}a_{0}^{4}},\quad \unicode[STIX]{x1D6FD}_{p}=\left(\frac{\unicode[STIX]{x2202}^{2}\unicode[STIX]{x1D70C}}{\unicode[STIX]{x2202}s^{2}}\right)_{p}=\frac{\unicode[STIX]{x1D70C}_{0}}{C_{p}^{2}},\end{eqnarray}$$

(5.3e ) $$\begin{eqnarray}\unicode[STIX]{x1D6FD}_{sp}=\left[\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}s}\left(\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}}{\unicode[STIX]{x2202}p}\right)_{s}\right]_{p}=\left[\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}s}\left(\frac{\unicode[STIX]{x1D70C}}{\unicode[STIX]{x1D6FE}p}\right)\right]_{p}=-\frac{1}{C_{p}a_{0}^{2}}.\end{eqnarray}$$

Figure 11.

Thermodynamic cycle in $p^{\prime }{-}\unicode[STIX]{x1D70C}^{\prime }$ plane for non-isentropic thermoacoustic wave amplification. (——), Nonlinear; (- - - -), linear. Sample perturbation fields are: $p^{\prime }=0.5p_{0}\cos \unicode[STIX]{x1D714}t$ and $s^{\prime }=-Rp^{\prime }/p_{0}+250\sin (\unicode[STIX]{x1D714}(t+\unicode[STIX]{x1D70F}))$ for $\unicode[STIX]{x1D70F}=h_{b}^{2}/\unicode[STIX]{x1D708}=0.0107~\text{s}$ .

The coefficients $\unicode[STIX]{x1D6FC}_{p}$ and $\unicode[STIX]{x1D6FC}_{s}$ contribute to first-order wave propagation and thermoacoustic effects while the second-order coefficients $\unicode[STIX]{x1D6FD}_{s},\unicode[STIX]{x1D6FD}_{p}$ and $\unicode[STIX]{x1D6FD}_{sp}$ in (5.2) account for the corresponding nonlinear effects. Hedberg & Rudenko (Reference Hedberg and Rudenko2011) have demonstrated the hysteretic effects of nonlinear wave propagation retaining only the $\unicode[STIX]{x1D6FD}_{s}$ nonlinear term and $\unicode[STIX]{x1D6FC}_{p}$ term accounting for irreversible entropy changes. Figure 11 shows the hysteresis cycle for the second-order constitutive relation in (5.2) in comparison with the first-order approximation. Nonlinear wave propagation outside the regenerator (sections $a$ , $c$ , and $d$ ) is not affected by the no-slip and isothermal boundary conditions. Only higher order irreversible entropy fluctuations are generated due to the axial conduction terms (Hamilton, Blackstock & others Reference Hamilton and Blackstock1998). Consequently, the constitutive relation (5.2) can be approximated with the terms corresponding to $\unicode[STIX]{x1D6FC}_{s},\unicode[STIX]{x1D6FC}_{p}$ and $\unicode[STIX]{x1D6FD}_{s}$ retained. However, inside the regenerator, wall shear and wall-heat flux from the no-slip isothermal boundaries generate first-order reversible entropy fluctuations; hence, in order to capture the nonlinear thermoacoustic wave amplification inside the regenerator, nonlinear terms in entropy perturbations need to be included, as shown in (5.2). Starting from the general second-order governing equations discussed above, the spatio-temporal evolution equations for the fluctuations in pressure $p^{\prime }$ and flow rate $U^{\prime }$ are derived in the following sections.

5.2 Free-shear nonlinear wave propagation

Waves outside the regenerator (sections $a$ , $c$ and $d$ ) in the minimal-unit set-up (figure 1) propagate in the absence of wall shear and wall-heat flux. As a result, terms involving the transverse gradient of $u^{\prime }$ in (5.1b ) and $T^{\prime }$ in (5.1c ) can be neglected and entropy fluctuations remain second order in the nonlinear regime as well (Hamilton et al. Reference Hamilton and Blackstock1998). Up to second order, the wave propagation in the duct sections $a$ , $c$ and $d$ is governed by

(5.4) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}p^{\prime }}{\unicode[STIX]{x2202}t}=-\frac{\unicode[STIX]{x1D6FE}p_{0}}{h}\left(1+\frac{1+\unicode[STIX]{x1D6FE}}{\unicode[STIX]{x1D6FE}}\frac{p^{\prime }}{p_{0}}\right)\frac{\unicode[STIX]{x2202}U^{\prime }}{\unicode[STIX]{x2202}x}+\frac{k}{\unicode[STIX]{x1D70C}_{0}}\left(\frac{1}{C_{v}}-\frac{1}{C_{p}}\right)\frac{\unicode[STIX]{x2202}^{2}p^{\prime }}{\unicode[STIX]{x2202}x^{2}}, & \displaystyle\end{eqnarray}$$

(5.5) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}U^{\prime }}{\unicode[STIX]{x2202}t}=-\frac{h}{\unicode[STIX]{x1D70C}_{0}}\frac{\unicode[STIX]{x2202}p^{\prime }}{\unicode[STIX]{x2202}x}+\unicode[STIX]{x1D708}_{0}\left(\unicode[STIX]{x1D709}_{B}+\frac{4}{3}\right)\frac{\unicode[STIX]{x2202}^{2}U^{\prime }}{\unicode[STIX]{x2202}x^{2}}, & \displaystyle\end{eqnarray}$$

where $h=h_{a},h_{c}$ or $h_{d}$ . Second-order nonlinearities in (5.4) cause the waveform distortion and steepening.

5.3 Wall shear and wall-heat flux

By design, the regenerator width is comparable to the local viscous and thermal Stokes layer thickness ( $h_{b}/2\unicode[STIX]{x1D6FF}_{\unicode[STIX]{x1D708}}\sim 1$ , $h_{b}/2\unicode[STIX]{x1D6FF}_{k}\sim 1$ ). Due to the wall shear, the velocity fluctuations inside the regenerator vary in the $y$ direction as well. Moreover, the wall-heat flux contributes to the first-order entropy fluctuations. Noting that nonlinear acoustic waves can be decomposed into acoustic, viscous and entropic modes (Chu & Kovásznay Reference Chu and Kovásznay1958; Pierce Reference Pierce1989), the following decomposition for the entropy and velocity fluctuations inside the regenerator is considered:

(5.6a,b ) $$\begin{eqnarray}u^{\prime }(x,y,t)=\tilde{u} (x,t)+u_{\unicode[STIX]{x1D708}}^{\prime }(x,y,t),\quad s^{\prime }=\tilde{s}(x,t)+s_{q}^{\prime }(x,y,t),\end{eqnarray}$$

where $u_{\unicode[STIX]{x1D708}}^{\prime }$ is the viscous velocity fluctuation and $\tilde{u}$ is the nonlinear acoustic wave field. The former is diffused by viscosity and is governed by the unsteady diffusion equation

(5.7a,b ) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}u_{\unicode[STIX]{x1D708}}^{\prime }}{\unicode[STIX]{x2202}t}=\unicode[STIX]{x1D708}_{0}\frac{\unicode[STIX]{x2202}^{2}u_{\unicode[STIX]{x1D708}}^{\prime }}{\unicode[STIX]{x2202}y^{2}},\quad u_{\unicode[STIX]{x1D708}}^{\prime }(x,y=\pm h_{b}/2,t)=-\tilde{u} (x,t).\end{eqnarray}$$

Similarly, $\tilde{s}$ accounts for the entropy changes due to the nonlinear acoustic wave propagation and $s_{q}^{\prime }$ corresponds to the first-order entropy changes due to the wall-heat flux inside the regenerator and is governed by the equation

(5.8a,b ) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}s_{q}^{\prime }}{\unicode[STIX]{x2202}t}+u^{\prime }\frac{\text{d}s_{0}}{\text{d}x}=\frac{\unicode[STIX]{x1D708}_{0}}{Pr}\frac{\unicode[STIX]{x2202}^{2}s_{q}^{\prime }}{\unicode[STIX]{x2202}y^{2}},\quad s_{q}^{\prime }(x,y=\pm h_{b}/2,t)=-\tilde{s}+s_{w}^{\prime },\end{eqnarray}$$

where $s_{w}^{\prime }$ corresponds to the entropy fluctuations at the isothermal walls driven by pressure fluctuations

(5.9) $$\begin{eqnarray}s_{w}^{\prime }=-\frac{1}{\unicode[STIX]{x1D70C}_{0}T_{0}}p^{\prime }=-\frac{R}{p_{0}}p^{\prime }.\end{eqnarray}$$

Equations (5.7) and (5.8) suggest the following infinite series solution forms for the viscous and the entropic fields:

(5.10a ) $$\begin{eqnarray}\displaystyle & & \displaystyle \displaystyle u_{\unicode[STIX]{x1D708}}^{\prime }=-\tilde{u} +\mathop{\sum }_{j=0}^{\infty }{\check{u}}_{j}(x,t)\cos (\unicode[STIX]{x1D701}_{j}y),\quad s_{q}^{\prime }=-\tilde{s}+s_{w}+\mathop{\sum }_{j=0}^{\infty }{\check{s}}_{j}(x,t)\cos (\unicode[STIX]{x1D701}_{j}y),\qquad\end{eqnarray}$$

(5.10b ) $$\begin{eqnarray}\displaystyle & & \displaystyle \qquad \text{with }\unicode[STIX]{x1D701}_{j}=(2j+1)\frac{\unicode[STIX]{x03C0}}{h_{b}}.\qquad\end{eqnarray}$$

Performing eigenfunction expansions along the $y$ direction (5.10a,b ) yields the following evolution equations for the Fourier coefficients corresponding to the viscous and entropic modes:

(5.11) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}{\check{u}}_{j}}{\unicode[STIX]{x2202}t}+\unicode[STIX]{x1D708}_{0}\unicode[STIX]{x1D701}_{j}^{2}{\check{u}}_{j}=(-1)^{j+1}\frac{2}{\unicode[STIX]{x1D701}_{j}h_{b}}\left(\frac{1}{\unicode[STIX]{x1D70C}_{0}}\frac{\unicode[STIX]{x2202}p^{\prime }}{\unicode[STIX]{x2202}x}\right), & \displaystyle\end{eqnarray}$$

(5.12) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}{\check{s}}_{j}}{\unicode[STIX]{x2202}t}+{\check{u}}_{j}\frac{\text{d}s_{0}}{\text{d}x}+\frac{\unicode[STIX]{x1D708}_{0}}{Pr}\unicode[STIX]{x1D701}_{j}^{2}{\check{s}}_{j}=(-1)^{j}\frac{2R}{\unicode[STIX]{x1D701}_{j}h_{b}p_{0}}\frac{\unicode[STIX]{x2202}p^{\prime }}{\unicode[STIX]{x2202}t}. & \displaystyle\end{eqnarray}$$

Equations (5.11) and (5.12) determine the evolution of the transverse modes of the longitudinal velocity ${\check{u}}_{j}$ and entropy ${\check{s}}_{j}$ fluctuations. In appendix C, the convergence properties of the infinite series in (5.10a,b ) for $j\rightarrow \infty$ are discussed.

5.4 Nonlinear thermoviscous wave equations

The axial velocity fluctuations $u^{\prime }$ are governed by (5.1b ) up to second-order accuracy. However, nonlinearities in (5.1b ) result in acoustic streaming, which is suppressed in the current analysis, and is therefore neglected (Hamilton et al. Reference Hamilton and Blackstock1998). Integrating the resulting momentum equation in $y$ and substituting (5.10a,b ) yields

(5.13) $$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x2202}U^{\prime }}{\unicode[STIX]{x2202}t}+\frac{h_{b}}{\unicode[STIX]{x1D70C}_{0}}\frac{\unicode[STIX]{x2202}p^{\prime }}{\unicode[STIX]{x2202}x}=\unicode[STIX]{x1D70F}_{w}^{\prime }+\frac{1}{\unicode[STIX]{x1D70C}_{0}}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}x}\left[\unicode[STIX]{x1D707}_{0}\left(\unicode[STIX]{x1D709}_{B}+\frac{4}{3}\right)\frac{\unicode[STIX]{x2202}U^{\prime }}{\unicode[STIX]{x2202}x}\right], & & \displaystyle\end{eqnarray}$$

where

(5.14) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D70F}_{w}^{\prime }=2\unicode[STIX]{x1D708}_{0}\mathop{\sum }_{j=0}^{\infty }(-1)^{j+1}{\check{u}}_{j}(x,t)\unicode[STIX]{x1D701}_{j}. & & \displaystyle\end{eqnarray}$$

Equation (5.4) governs the evolution of the pressure fluctuations up to second order in the free-shear/adiabatic ducts. In order to derive an analogous governing equation for the regenerator, (5.1a ), (5.1c ) and (5.2) are combined to obtain

(5.15) $$\begin{eqnarray}\displaystyle \underbrace{\frac{\unicode[STIX]{x2202}p^{\prime }}{\unicode[STIX]{x2202}t}+\frac{\unicode[STIX]{x1D70C}_{0}a_{0}^{2}}{h_{b}}\frac{\unicode[STIX]{x2202}U^{\prime }}{\unicode[STIX]{x2202}x}}_{wave~propagation}=\frac{\unicode[STIX]{x1D70C}_{0}a_{0}^{2}}{h_{b}}\left(\frac{q^{\prime }}{C_{p}\unicode[STIX]{x1D70C}_{0}T_{0}}+\underbrace{q_{2}\,+\,\mathbb{T}\,-\,\mathbb{Q}}_{\substack{ thermodynamic \\ nonlinearities}}\,+\,\mathbb{D}_{s}\right)-\mathbb{C}, & & \displaystyle\end{eqnarray}$$

where

(5.16) $$\begin{eqnarray}q^{\prime }=\frac{2\unicode[STIX]{x1D708}_{0}\unicode[STIX]{x1D70C}_{0}T_{0}}{Pr}\mathop{\sum }_{j=0}^{\infty }(-1)^{j+1}{\check{s}}_{j}(x,t)\unicode[STIX]{x1D701}_{j},\end{eqnarray}$$

defines the fluctuating wall-heat flux and couples the pressure evolution (5.15) with the entropic mode evolution (5.12), whereas $\mathbb{Q}$ denotes the nonlinear interaction of pressure and wall-heat flux fluctuations,

(5.17) $$\begin{eqnarray}\mathbb{Q}=\frac{\unicode[STIX]{x1D6FE}p^{\prime }q^{\prime }}{C_{p}p_{0}\unicode[STIX]{x1D70C}_{0}T_{0}},\end{eqnarray}$$

hereafter referred as macrosonic thermoacoustic interaction. The term denoted by $q_{2}$ corresponds to the second-order heat flux which is a quadratic function of the entropy gradient in $y$ . The nonlinear terms denoted by the double faced $\mathbb{T}$ correspond to the constitutive (thermodynamic) nonlinearities which account for the second-order density fluctuation due to first-order entropic modes. The terms denoted by the double faced $\mathbb{C}$ correspond to the convective nonlinearities in (5.1b ) and (5.1c ) and those denoted by $\mathbb{D}_{s}$ account for the axial diffusion of gradients in highly nonlinear regimes of thermoacoustic wave amplification. A detailed derivation of (5.15) and expressions for the terms $q_{2}$ , $\mathbb{D}_{s}$ , $\mathbb{T}$ , and $\mathbb{C}$ are given in appendix C. Equations (5.11)–(5.13) and (5.15) constitute the governing equations for the spatio-temporal evolution of large amplitude acoustic perturbations inside the regenerator. The macrosonic thermoacoustic interaction (5.17) breaks the thermodynamic symmetry between the interactions of compressions and dilatations with the wall-heat flux inside the regenerator, thus highlighting that the entropy of a Lagrangian parcel of fluid changes by a small amount under high amplitude compressions ( $\unicode[STIX]{x1D70C}^{\prime }>0$ ), compared to dilatations ( $\unicode[STIX]{x1D70C}^{\prime }<0$ ), for the same amount of heat input or output.

In general, for thermoacoustic devices in looped configuration, the length of the regenerator is very short compared to the total length of the device. As a result, higher-order terms affecting only the propagation of the acoustic perturbations, such as convective nonlinearities, can be neglected inside the regenerator. Under such hypotheses, the following approximate nonlinear governing equation for the pressure fluctuations $p^{\prime }$ inside a short regenerator is obtained:

(5.18) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}p^{\prime }}{\unicode[STIX]{x2202}t}+\frac{\unicode[STIX]{x1D70C}_{0}a_{0}^{2}}{h_{b}}\frac{\unicode[STIX]{x2202}U^{\prime }}{\unicode[STIX]{x2202}x}\approx \frac{\unicode[STIX]{x1D70C}_{0}a_{0}^{2}}{h_{b}}\left\{\frac{1}{C_{p}}\left[\left(1-\frac{\unicode[STIX]{x1D6FE}p^{\prime }}{p_{0}}\right)\frac{q^{\prime }}{\unicode[STIX]{x1D70C}_{0}T_{0}}\right]\right\}.\end{eqnarray}$$

In the above equation, terms $\mathbb{T}$ , $\mathbb{D}_{s}$ and $q_{2}$ are neglected for simplicity. Equations (5.11) and (5.12) can also be integrated in time analytically to express the wall shear $\unicode[STIX]{x1D70F}_{w}^{\prime }$ and wall-heat flux $q^{\prime }$ in terms of acoustic variables. The time integration of (5.11) and (5.12) yields:

(5.19) $$\begin{eqnarray}\displaystyle & \displaystyle {\check{u}}_{j}=(-1)^{j+1}\frac{2}{\unicode[STIX]{x1D701}_{j}h_{b}\unicode[STIX]{x1D70C}_{0}}\int _{-\infty }^{t}\text{e}^{-(t-\unicode[STIX]{x1D702})/\unicode[STIX]{x1D70F}_{j}}\frac{\unicode[STIX]{x2202}p^{\prime }}{\unicode[STIX]{x2202}x}(x,\unicode[STIX]{x1D702})\,\text{d}\unicode[STIX]{x1D702}, & \displaystyle\end{eqnarray}$$

(5.20) $$\begin{eqnarray}\displaystyle & \displaystyle {\check{s}}_{j}=-\frac{\text{d}s_{0}}{\text{d}x}\int _{-\infty }^{t}\text{e}^{-(t-\unicode[STIX]{x1D702})/Pr\unicode[STIX]{x1D70F}_{j}}{\check{u}}_{j}(x,\unicode[STIX]{x1D702})\,\text{d}\unicode[STIX]{x1D702}+(-1)^{j}\frac{2R}{\unicode[STIX]{x1D701}_{j}hp_{0}}\int _{-\infty }^{t}\text{e}^{-(t-\unicode[STIX]{x1D702})/Pr\unicode[STIX]{x1D70F}_{j}}\frac{\unicode[STIX]{x2202}p^{\prime }}{\unicode[STIX]{x2202}\unicode[STIX]{x1D702}}(x,\unicode[STIX]{x1D702})\,\text{d}\unicode[STIX]{x1D702},\qquad & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D70F}_{j}=1/\unicode[STIX]{x1D708}_{0}\unicode[STIX]{x1D701}_{j}^{2}$ defines the viscous relaxation time for the $j$ th viscous mode. Hence, writing the relaxation functional for some function $\unicode[STIX]{x1D719}(x,t)$ , namely,

(5.21) $$\begin{eqnarray}{\mathcal{G}}_{j}(\unicode[STIX]{x1D719},\unicode[STIX]{x1D70F}_{j})=\int _{-\infty }^{t}\text{e}^{-(t-\unicode[STIX]{x1D702})/\unicode[STIX]{x1D70F}_{j}}\unicode[STIX]{x1D719}(x,\unicode[STIX]{x1D702})\,\text{d}\unicode[STIX]{x1D702},\end{eqnarray}$$

and summing (5.19) and (5.20) over $j$ , the following expressions for the wall shear and wall-heat flux are obtained:

(5.22) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D70F}_{w}^{\prime }=\frac{4\unicode[STIX]{x1D708}_{0}}{\unicode[STIX]{x1D70C}_{0}h_{b}}\mathop{\sum }_{j=0}^{\infty }{\mathcal{G}}_{j}\left(\frac{\unicode[STIX]{x2202}p^{\prime }}{\unicode[STIX]{x2202}x},\unicode[STIX]{x1D70F}_{j}\right), & \displaystyle\end{eqnarray}$$

(5.23) $$\begin{eqnarray}\displaystyle & \displaystyle q^{\prime }=\frac{2\unicode[STIX]{x1D70C}_{0}\unicode[STIX]{x1D708}_{0}T_{0}}{Pr}\mathop{\sum }_{j=0}^{\infty }\left[(-1)^{j}\frac{\text{d}s_{0}}{\text{d}x}{\mathcal{G}}_{j}(\unicode[STIX]{x1D701}_{j}{\check{u}}_{j},\unicode[STIX]{x1D70F}_{j}Pr)-\frac{2R}{h_{b}p_{0}}{\mathcal{G}}_{j}\left(\frac{\unicode[STIX]{x2202}p^{\prime }}{\unicode[STIX]{x2202}t},\unicode[STIX]{x1D70F}_{j}Pr\right)\right]. & \displaystyle\end{eqnarray}$$

Equations (5.22) and (5.23) provide first-order expressions for the wall shear and the wall-heat flux as a function of a generic acoustic field near the walls and, together with (5.13) and (5.18), complete the nonlinear wave propagation model equations. However, in the present work, (5.11) and (5.12) have been considered for time integration for simplicity.

Sugimoto (Reference Sugimoto2010) systematically derived the functional form of ${\mathcal{G}}_{j}$ for extremely thick and thin diffusion layers using the linear acoustic field approximation. The functional ${\mathcal{G}}_{j}$ in (5.21) approximates the wall shear and the wall-heat flux in terms of any acoustic field, linear or nonlinear, though within the restriction of linear decomposition of the field into viscous and entropic modes. Recently, Sugimoto (Reference Sugimoto2016) also developed a theoretical framework to elucidate high amplitude nonlinear wave propagation in a shear-dominated duct with the restriction of very thick diffusion layers and focusing on streaming. However, the wall shear and wall-heat flux expressed in terms of the acoustic field variables by the relaxation functional ${\mathcal{G}}_{j}$ in (5.22) and (5.23) hold true irrespective of the relative thickness between the diffusion layers and the channel width.

5.5 Solution technique, minor losses and shock capturing

The model equations (5.4), (5.5), (5.11), (5.12), (5.13) and (5.15) are integrated in time utilizing a fourth-order explicit Runge–Kutta advancement with a second-order staggered spatial discretization. Appropriate interface conditions for ducts with abrupt area jumps are provided as well. In general, abrupt area jumps cause minor losses in the pressure ( $\unicode[STIX]{x0394}p_{ml}$ ) due to nonlinear vortex generation which results in second-order losses. Hence, minor losses can be modelled as quadratic functions of the first-order interface flow rate, i.e.  $\unicode[STIX]{x0394}p_{ml}\sim U_{i}^{2}$ . In analogy with the approach followed by Lin et al. (Reference Lin, Scalo and Hesselink2016), in order to match the difference between the left-hand pressure ( $p_{i}^{-}$ ) and the right-hand pressure ( $p_{i}^{+}$ ) limits at the interface to these minor losses, the linearly extrapolated pressures from the closest discretization points are used, as schematically shown in figure 12. Accordingly, the interface condition imposed reads:

(5.24) $$\begin{eqnarray}U_{i}^{(2)}:p_{i}^{-}-p_{i}^{+}=({\textstyle \frac{3}{2}}p_{1}-{\textstyle \frac{1}{2}}p_{0})-({\textstyle \frac{3}{2}}p_{2}-{\textstyle \frac{1}{2}}p_{3})=\unicode[STIX]{x0394}p_{ml}\approx {\textstyle \frac{1}{2}}\unicode[STIX]{x1D70C}_{0}KU_{i}^{(1)^{2}},\end{eqnarray}$$

where $K$ represents an empirically determined minor loss coefficient for an incompressible steady flow and $U_{i}^{(1)}$ corresponds to the first-order interface flow rate calculated such that

(5.25) $$\begin{eqnarray}U_{i}^{(1)}:p_{i}^{-}=p_{i}^{+}.\end{eqnarray}$$

Hence, the second-order flow rate $U_{i}^{(2)}$ is calculated imposing the condition in (5.24). In the limit of shock-wave formation, the spectral energy cascade results in the formation of large gradients, which reach the limit of grid resolution. Local Legendre polynomial expansions of the pressure field are thus used to introduce an artificial bulk viscosity which increases the effective diffusion coefficient in the direction of propagation. Details of such a windowed artificial viscosity implementation are discussed in appendix D.

Figure 12.

Solution technique for integrating quasi-one-dimensional governing equations in time across abrupt area changes. Interface bulk velocity at $n$ time step, $U_{i}^{n}$ is calculated such that difference of left and right extrapolations of pressure ( $p_{i}^{-}$ and $p_{i}^{+}$ respectively) equals minor loss $\unicode[STIX]{x0394}p_{ml}$ at $n+1$ time step which is a function of $U_{i}^{n}$ .

7.1 Phase evolution

As discussed in § 4.1, the cycle-averaged production of the acoustic energy density $\overline{{\mathcal{P}}}$ is governed by the relative phasing between pressure and velocity perturbations $p^{\prime }$ and $u^{\prime }$ , respectively. Figure 15(a–d) shows the phasing of $p^{\prime }$ and $u^{\prime }$ extracted at $x=20~\text{mm}$ , $y=0~\text{mm}$ and compares the evolution of the relative phase angle $\unicode[STIX]{x1D713}_{up}$ obtained from the time-domain nonlinear model and the fully compressible Navier–Stokes simulations using the correlation function expression

(7.1) $$\begin{eqnarray}\cos (\unicode[STIX]{x1D713}_{up})=\frac{\displaystyle \int _{0}^{t}p^{\prime }u^{\prime }\,\text{d}\unicode[STIX]{x1D70F}}{\displaystyle \sqrt{\int _{0}^{t}{p^{\prime }}^{2}\,\text{d}\unicode[STIX]{x1D70F}}\sqrt{\int _{0}^{t}{u^{\prime }}^{2}\,\text{d}\unicode[STIX]{x1D70F}}}.\end{eqnarray}$$

In general, the above relation can be used to study the evolution of the phase difference between close to harmonically oscillating quantities (Balasubramanian & Sujith Reference Balasubramanian and Sujith2008). The relative phase $\unicode[STIX]{x1D713}_{up}$ saturates at a value of approximately $40^{\circ }$ , which is also the value in the harmonic growth regime indicating production of the energy density at the limit cycle inside the regenerator. However, due to large amplitude perturbations at the limit cycle, the production is attenuated due to the macrosonic thermoacoustic interaction term (5.17) discussed in the previous section. The net production of the energy density inside the regenerator in the nonlinear regime can be approximated by combining equations (5.13) and (5.18) such that

(7.2) $$\begin{eqnarray}{\mathcal{P}}-{\mathcal{D}}\approx \frac{p^{\prime }}{h_{b}}\left(\frac{q^{\prime }}{C_{p}\unicode[STIX]{x1D70C}_{0}T_{0}}-\mathbb{Q}\right)+\frac{\unicode[STIX]{x1D70F}_{w}^{\prime }U^{\prime }}{h_{b}^{2}},\end{eqnarray}$$

which shows a decrease, compared to the linear regime, due to the macrosonic thermoacoustic interaction term $\mathbb{Q}$ .

Figure 15.

Evolution of pressure perturbation $p^{\prime }$ (a), axial velocity perturbation $u^{\prime }$ at $y=0~\text{mm}$ (b), phase portrait in the $p^{\prime }{-}u^{\prime }$ plane at the limit cycle (c), acoustic phasing $\unicode[STIX]{x1D713}_{pu}$ at $x=20~\text{mm}$ inside the regenerator, phase portrait in the $p^{\prime }{-}q^{\prime }$ plane at the limit cycle (e) and phase $\unicode[STIX]{x1D713}_{pq}$ evolution inside the regenerator at $x=20~\text{mm}$ (f). (d) (——), Nonlinear thermoacoustic model; (○), Navier–Stokes simulations.

Figure 16.

Scaled spectral energy density $E_{k}{\mathcal{E}}^{-3/2}$ at the limit cycle against frequency of the harmonics. (——), Time-domain nonlinear model at $T_{H}=450$ ; Navier–Stokes simulations at (○), $T_{H}=450~\text{K}$ ; ( $+$ ),  $T_{H}=500~\text{K}$ ; (▫), $T_{H}=550~\text{K}$ ; (- - - -) compares the variation of the energy with harmonic frequency with the power law derived in (7.7).

Furthermore, figure 15(f) shows that the phase angle $\unicode[STIX]{x1D713}_{pq}$ saturates at a value $\unicode[STIX]{x1D713}_{pq}\approx 33.57^{\circ }$ compared to $\unicode[STIX]{x1D713}_{pq}\approx 7.41^{\circ }$ in the linear regime, thus showing the decrease in the production of the energy density at the limit cycle. However, since the phase angle $\unicode[STIX]{x1D713}_{pq}\neq 90^{\circ }$ at the limit cycle, the production of acoustic energy density continues and shock waves are therefore sustained.

7.2 Scales of acoustic spectral energy cascade

At the limit cycle, the energy of the unstable mode continues to increase and further cascades into higher harmonics on account of the nonlinear wave propagation. Higher harmonics have correspondingly shorter wavelengths due to which gradients in the longitudinal direction $x$ become large. Through the bulk viscosity and the thermal conductivity, the energy density is dissipated at higher harmonics, thus establishing a steady flow of energy from the unstable mode to the higher harmonics. The distribution of the spectral energy density in the harmonics can be derived utilizing an energy cascade modelling, analogous to the turbulent energy cascade (Pope Reference Pope2000; Nazarenko Reference Nazarenko2011). Assuming travelling-wave propagation, the total energy density of the $k$ th harmonic associated with planar wave propagation, $E_{k}^{(1\text{D})}$ , can be defined as

(7.3a,b ) $$\begin{eqnarray}E^{(1\text{D})}=\int _{0}^{2\unicode[STIX]{x03C0}a_{0}/\unicode[STIX]{x1D714}_{0}}E\,\text{d}\unicode[STIX]{x1D709}=\mathop{\sum }_{k}E_{k}^{(1\text{D})},\quad E_{k}^{(1\text{D})}=\frac{\unicode[STIX]{x03C0}a_{0}}{\unicode[STIX]{x1D714}_{0}}E_{k},\end{eqnarray}$$

where $E$ is the squared $L_{2}$ norm defined in (3.2) where $p^{\prime }$ and $u^{\prime }$ are defined as a function of the travelling-wave coordinate $\unicode[STIX]{x1D709}$ . At the limit cycle, the rate of cascade of energy ${\mathcal{E}}$ in the spectral space balances the thermoacoustic wave amplification. As a result, higher harmonics (namely, the overtones of the unstable mode) are generated and later dissipated by molecular dissipation, as indicated by the presence of the molecular dissipation factor $\unicode[STIX]{x1D6FF}$ in (6.2). Thus, ${\mathcal{E}}$ is purely governed by thermoacoustic wave amplification at smaller harmonics and viscous dissipation at higher harmonics and scales as

(7.4) $$\begin{eqnarray}\frac{{\mathcal{E}}}{\unicode[STIX]{x1D6FC}_{eff}^{2}\unicode[STIX]{x1D6FF}}=\text{const.},\end{eqnarray}$$

where $\unicode[STIX]{x1D6FC}_{eff}$ is the effective energy amplification rate at the limit cycle. Based on the macrosonic thermoacoustic interaction, $\unicode[STIX]{x1D6FC}_{eff}$ can be estimated as

(7.5) $$\begin{eqnarray}\unicode[STIX]{x1D6FC}_{eff}\approx (1-\unicode[STIX]{x1D6FE}|p^{\prime }|/p_{0})\unicode[STIX]{x1D6FC},\end{eqnarray}$$

where $\unicode[STIX]{x1D6FC}$ is the growth rate of the unstable mode in the linear regime, and $|p^{\prime }|$ is the amplitude of the pressure fluctuation at the limit cycle. Assuming that the total energy per unit mass $E_{k}^{(1\text{D})}/\unicode[STIX]{x1D70C}_{0}$ only depends on the rate of energy cascade ${\mathcal{E}}$ and on the angular frequency of the harmonic $\unicode[STIX]{x1D714}_{k}$ , the following scaling is obtained:

(7.6) $$\begin{eqnarray}\frac{E_{k}^{(1\text{D})}\unicode[STIX]{x1D714}_{k}^{5/2}}{\unicode[STIX]{x1D70C}_{0}{\mathcal{E}}^{3/2}}=\text{const}.\end{eqnarray}$$

Substituting into (7.3) and eliminating $\unicode[STIX]{x1D70C}_{0}$ , the scaling for the spectral energy density $E_{k}$ remains the same such that

(7.7) $$\begin{eqnarray}\frac{E_{k}\unicode[STIX]{x1D714}_{k}^{5/2}}{{\mathcal{E}}^{3/2}}=\text{const}.\end{eqnarray}$$

The scaling derived in (7.7) shows that the distribution of the energy density at limit cycle in the spectral space decays as $\unicode[STIX]{x1D714}_{k}^{-5/2}$ – where $\unicode[STIX]{x1D714}_{k}$ is the frequency of the harmonic – as shown in figure 16. Moreover, in the spectral broadening regime, the energy density in the spectral space varies with the effective growth rate as $\unicode[STIX]{x1D6FC}_{eff}^{3}$ (figure 6) which, in turn, allows the comparison among cases with varying hot-side temperature $T_{H}$ . Such a scaling arises purely from the mechanism of the nonlinear saturation resulting from the spectral energy cascade. It is however challenging to accurately estimate the acoustic energy production (quantified by $\unicode[STIX]{x1D6FC}_{eff}^{3}$ ) at the limit cycle utilizing the time-domain nonlinear model, or even the fully compressible Navier–Stokes simulations, due to the shock-capturing artificial viscosity that has to be added to ensure numerical stability in the presence of shocks. Overall, the spectral cascade of the energy density is balanced by the thermoacoustic wave amplification and the viscous dissipation, i.e. the shock waves are thermoacoustically sustained. However, as shown previously (Gupta, Scalo & Lodato Reference Gupta, Scalo and Lodato2017), at the location of the steepest gradient in pressure, the wall-heat flux from the acoustic field inside the isothermal walls and wall shear stress are maximum inside the regenerator. Thus, the propagation of sharp gradients in the acoustic field $p^{\prime }$ and $U^{\prime }$ inside the regenerator are counteracted by the wall shear stresses and the wall-heat flux resulting in smoothing of shock waves. Consequently, the scaling argument (7.7) can be further improved accounting for propagation of very large harmonics inside the regenerator.