Spectral energy cascade in thermoacoustic shock waves

TL;DR

This paper investigates the nonlinear development of thermoacoustic waves in a resonator, revealing a transition from harmonic growth to shock-dominated regimes through spectral energy cascade, using theoretical and numerical methods.

Contribution

It introduces a comprehensive model of thermoacoustic wave regimes based on first principles, highlighting the nonlinear spectral cascade leading to shock formation.

Findings

Identification of three distinct wave regimes

Observation of nonlinear spectral broadening

Demonstration of shock-wave limit cycle saturation

Abstract

We have performed a theoretical and numerical investigation of thermoacoustically amplified quasi-planar nonlinear waves driven to the limit of shock-wave formation in a variable-area 2.58 m long looped resonator designed via eigenvalue analysis to maximize the growth rate of the second harmonic (265 Hz). High-order unstructured fully compressible Navier-Stokes simulations reveal the presence of three regimes: (i) Monochromatic harmonic growth, governed by linear thermoacoustics; (ii) Hierarchical spectral broadening, characterized by nonlinear energy cascade; (iii) Shock-wave dominated limit cycle, where energy production is balanced by dissipation occurring at the captured shock-thickness scale. We have modeled these regimes of thermoacoustic wave amplification based on first-principle governing equations and elucidated the physics of energy density production and saturation to a limit cycle.

Tables

Table 1. 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 subscript 𝑃 0 subscript 𝑃 ref P_{0}=P_{\mathrm{ref}} , ρ c = ρ ref subscript 𝜌 𝑐 subscript 𝜌 ref \rho_{c}=\rho_{\mathrm{ref}} , T C = T ref subscript 𝑇 𝐶 subscript 𝑇 ref T_{C}=T_{\mathrm{ref}} .

$P_0$	${\rho }_C$	$T_C$	$T_H$
101325 Pa	1.176 kg/m3	300 K	400 K, 450 K,
			500 K, 550 K

Table 2. Table 2 : Discretization order used per element, p 𝑝 p , total number of elements, N el subscript 𝑁 el N_{\mathrm{el}} , degrees of freedom, N dof = p 2 ​ N el subscript 𝑁 dof superscript 𝑝 2 subscript 𝑁 el N_{\mathrm{dof}}=p^{2}N_{\mathrm{el}} , and minimum and maximum linear element size Δ h ( min ) subscript subscript Δ ℎ min {\Delta_{h}}_{\mathrm{(min)}} , Δ h ( max ) subscript subscript Δ ℎ max {\Delta_{h}}_{\mathrm{(max)}} (cf. ( 9 )). Reported values are the ones used at the limit cycle. In all cases, simulations are carried out with p = 3 𝑝 3 p=3 and Grid A throughout the transient.

$p=3$
$p=5$
$p=7$

Table 3. Table 3 : Frequency and growth rates of unstable mode and its first seven overtones as obtained from Navier-Stokes simulations (N-S) and time integration of the nonlinear model (Model).

N-S

Model

Equations

$${\delta }_k=\sqrt{\frac{2\nu }{\omega Pr}},$$

$$\frac{\partial }{\partial t}\left(\rho \right)$$

$$\frac{\partial }{\partial t}\left(\rho u_i\right)$$

$$\frac{\partial }{\partial t}\left(\rho E\right)$$

$$f_B=\frac{\alpha }{\Delta t}\left(\dot{m}-{\dot{m}}_0\right),$$

$${\tau }_{ij}=2\mu \left[S_{ij}-\frac{1}{3}\frac{\partial u_k}{\partial x_k}{\delta }_{ij}\right],q_j$$

$$\xi =c_{s}t-x,$$

$$\frac{\Delta s}{R}=\sum _i\frac{\Delta s_i}{R}\approx \frac{\gamma +1}{12{\gamma }^2}\sum _i{\left(\frac{\Delta p_i}{p_i}\right)}^3,$$

$$\frac{\Delta s}{R}\approx \frac{2\gamma }{(\gamma +1{)}^2}\frac{(M^2-1{)}^3}{3}.$$

$$M(r_{\Delta })=m_0+m_1{r}_{\Delta }+m_2{r}_{\Delta }^2+m_3{r}_{\Delta }^3,$$

$$r_{\Delta }={\overline{\Delta }}_h/{\overline{\Delta }}_h{|}_{p=7}^B,{\overline{\Delta }}_h=\sqrt{\Omega /N_{\mathrm{dof}}},$$

$$\frac{\partial \mathbf{X}}{\partial t}=\mathbf{f}(\mathbf{X}),\mathrm{with}\mathbf{X}={\left(\frac{u^{\prime }}{a_0},\frac{p^{\prime }}{{\rho }_0{a}_0^2}\right)}^{\mathrm{T}},$$

$$E=\frac{1}{2}{\rho }_0{a}_0^2{\mathbf{X}}^{\mathrm{T}}\mathbf{X}=\frac{1}{2}{\rho }_0{u}^{\prime 2}+\frac{p^{\prime 2}}{2{\rho }_0{a}_0^2}.$$

$$u^{\prime }(x,t)=\sum _{\begin{array}{c}k=-\infty \\ k\ne 0\end{array}}^{+\infty }{\hat{u}}_k(x,ϵt)e^{i\frac{k}{2}{\omega }_{0}t},$$

$${\hat{u}}_{-k}={\hat{u}}_k^{\ast },$$

$$\overline{E}=2\sum _{k=1}^{\infty }{E}_k,E_k=\frac{1}{2}{\rho }_0|{\hat{u}}_k{|}^2+\frac{|{\hat{p}}_k{|}^2}{2{\rho }_0{a}_0^2},$$

$$\overline{(\cdot )}=\frac{1}{T_0}{\int }_{ϵt}^{ϵt+T_0}(\cdot )dt,T_0=\frac{2\pi }{{\omega }_0},$$

$${\alpha }_k=\alpha k/2,k\in \{4,6,8,\dots \}.$$

$$\frac{\partial {\rho }^{\prime }}{\partial t}+{\rho }_0\frac{\partial u^{\prime }}{\partial x}+u^{\prime }\frac{d{\rho }_0}{dx}+{\rho }_0\frac{\partial v^{\prime }}{\partial y}=0,$$

$$\frac{\partial u^{\prime }}{\partial t}+\frac{1}{{\rho }_0}\frac{\partial p^{\prime }}{\partial x}-{\nu }_0\frac{{\partial }^2{u}^{\prime }}{\partial y^2}=0,$$

$$\frac{\partial s^{\prime }}{\partial t}+u^{\prime }\frac{d{s}_0}{dx}-\frac{k_0}{{\rho }_0{T}_0}\frac{{\partial }^2{T}^{\prime }}{\partial y^2}=0,$$

$$\frac{\partial p^{\prime }}{\partial t}+\frac{{\rho }_0{a}_0^2}{h}\frac{\partial U^{\prime }}{\partial x}=\frac{{\rho }_0{a}_0^2}{h}\frac{q^{\prime }}{{\rho }_0{C}_p{T}_0},$$

$$\frac{1}{h}\frac{\partial U^{\prime }}{\partial t}+\frac{1}{{\rho }_0}\frac{\partial p^{\prime }}{\partial x}=\frac{1}{h}\frac{{\tau }_w^{\prime }}{{\rho }_0},$$

$$U^{\prime }={\int }_{-h/2}^{+h/2}{u}^{\prime }(x,y,t)dy,$$

$$q^{\prime }=2k_0{\frac{\partial T^{\prime }}{\partial y}|}_{y=+h/2},{\tau }_w^{\prime }=2{\mu }_0{\frac{\partial u^{\prime }}{\partial y}|}_{y=+h/2}.$$

$$p^{\prime }(x,t)=\hat{p}(x)e^{\sigma t},U^{\prime }(x,t)=\hat{U}(x)e^{\sigma t},$$

$$\sigma \hat{p}=\frac{{\rho }_0{a}_0^2}{h_b}\left[\frac{1}{1+\left(\gamma -1\right)f_k}\left(\frac{\Theta (f_k-f_{\nu })}{(1-f_{\nu })(1-Pr)}-\frac{d}{dx}\right)\right]\hat{U},\Theta =\frac{1}{T_0}\frac{d{T}_0}{dx},$$

$$\sigma \hat{U}=-\frac{h_b}{{\rho }_0}\left(1-f_{\nu }\right)\frac{d\hat{p}}{dx},$$

$$f_{\nu }=\frac{\tanh (\eta h_b/2)}{\eta h_b/2},f_k=\frac{\tanh (\eta h_b\sqrt{Pr}/2)}{\eta h_b\sqrt{Pr}/2},\mathrm{with}\eta =\sqrt{i\omega /{\nu }_0}.$$

$$\frac{\partial E}{\partial t}+\frac{\partial \mathcal{I}}{\partial x}=\mathcal{P}-\mathcal{D},$$

$$E=\frac{1}{2}\frac{p^{\prime 2}}{{\rho }_0{a}_0^2}+\frac{1}{2}{\rho }_0{\left(\frac{U^{\prime }}{h_b}\right)}^2,$$

$$\mathcal{I}=\frac{p^{\prime }{U}^{\prime }}{h_b},\mathcal{P}-\mathcal{D}=\frac{p^{\prime }{q}^{\prime }}{h_b{\rho }_0{C}_p{T}_0}+\frac{{\tau }_w^{\prime }{U}^{\prime }}{h_b^2}.$$

$$\frac{d}{d(ϵt)}{\int }_L\overline{E}(x,ϵt)dx={\int }_L\left(\overline{\mathcal{P}}-\overline{\mathcal{D}}\right)dx=\mathcal{R},$$

$${\hat{\tau }}_w=h_b\frac{\partial \hat{p}}{\partial x}{f}_{\nu },$$

Full text

1 Introduction

Thermoacoustic amplification of waves in a compressible flow is the result of a fluid dynamic instability emerging from the favourable coupling between pressure and heat-release fluctuations. The wavelength of thermoacoustically unstable waves is set by the size of the enclosing resonant chamber, while the heat release, providing the energy source for the amplification, is confined to a compact region. The heat release rate is a function of local velocity and pressure fluctuations, affecting for example the instantaneous flame surface area in a combustion chamber or the rate of convective heat extraction from a hot wire-mesh screen in a Rijke tube. The resulting fluctuations in the heat release rate drive a cycle of compressions and dilatations, which act as a source for pressure fluctuations that, if within a quarter phase from the heat release itself, become thermoacoustically amplified (Rayleigh 1878’s criterion). In combustion applications, efforts are made to suppress thermoacoustic amplification of waves (Poinsot & Veynante, 2011), which is, on the other hand, desired in thermoacoustic waste heat recovery devices (Swift, 1988).

The goal of the present manuscript is to investigate the high amplitude (or macrosonic) limit of thermoacoustically driven nonlinear waves characterized by the formation of self-sustaining resonating shock waves and inter-harmonic, or inter-scale, energy transfer. In thermoacoustically unstable combustors, modelling based on a modally truncated analysis of nonlinearly coupled oscillators has been shown to successfully describe the energy transfer among harmonics (Culick, 2006; Ananthkrishnan et al., 2005). Culick (1971) also formalized the procedure for estimating the nonlinear growth and decay of the thermoacoustically unstable mode and a few companion overtones (modally truncated analysis). However, in the limit of self-sustaining resonating shock waves, harmonics with very short wavelengths are generated as a result of energy cascade, rendering modally truncated analysis insufficient. While the problem of shock wave resonance sustained by near-resonant frequency excitation has been analysed in detail (Saenger & Hudson, 1960; Chester, 1964), self-sustaining shocks generated via modal instability remain to be analysed in detail to the best of authors’ knowledge.

Multiple authors (Swift, 1992; Karpov & Prosperetti, 2000; Hamilton et al., 2002; Penelet et al., 2005) have discussed the higher harmonic generation in non-combustion-driven thermoacoustically unstable waves, restricting their analysis to only a first few overtones of the unstable mode. Biwa et al. (2011) reported the first experimental observation of energy cascade in thermoacoustically sustained shock waves in a looped thermoacoustic resonator. The stability limits of a similar experimental setup were studied previously by Yazaki et al. (1998), who generated second-mode quasi-travelling wave instability, although no shock formation was reported. Olivier et al. (2015) performed numerical modelling of a looped thermoacoustic resonator similar (but not identical to) Biwa et al. (2011)’s apparatus, combining approximate nonlinear propagation equations (Menguy & Gilbert, 2000) in the resonator with linear thermoviscous equations (Rott, 1969, 1973) in the thermoacoustic regenerator, i.e. the heat exchanger used to sustain the wall temperature gradient required to generate the instability. The approximate equations by Menguy & Gilbert (2000), however, assume low frequencies, hence are invalid over the entire spectrum generated by a nonlinear energy cascade driven up to shock wave formation, and neglect thermoacoustic nonlinearities in the regenerator, which, as shown later, need to be accounted for.

In the present work, a comprehensive nonlinear theoretical and high-fidelity modelling approach is adopted to accurately describe macrosconic thermoacoustic waves. To this end, a canonical travelling-wave looped resonator, inspired by Yazaki et al. (1998)’s experimental setup but geometrically optimized via linear theory (Rott, 1969, 1973; Lin et al., 2016), has been designed to maximize the growth rate of the quasi-travelling-wave second harmonic and thus achieve rapid shock wave formation. Yazaki et al. (1998)’s looped configuration allows quasi-travelling-wave acoustic phasing which facilitates faster nonlinear energy cascade compared to standing wave resonators (Biwa et al., 2014). It is shown that the energy content in spectral domain resembles the equilibrium energy cascade observed in turbulence, similar to the spectral energy distribution of an ensemble of acoustic waves interacting nonlinearly among each other (Nazarenko, 2011; Zakharov et al., 2012). As demonstrated by the numerical simulation data and companion low-order nonlinear modelling, thermoacoustically sustained shock waves exhibit inter-scale energy transfer dynamics analogous to Kolmogorov’s equilibrium hydrodynamic turbulent energy cascade (Kolmogorov, 1941). Throughout, the results obtained via the proposed nonlinear model are verified and compared with fully compressible high-fidelity Navier-Stokes simulations.

The development of an accurate nonlinear thermoacoustic wave propagation theory and modelling framework warrants the support of high-fidelity numerical simulations. A high-order spectral difference numerical framework (Kopriva & Kolias, 1996; Kopriva, 1996; Sun et al., 2007; Jameson, 2010) combined with an artificial Laplacian viscosity (Persson & Peraire, 2006) shock-capturing scheme has been adopted for the present study. Moreover, the computational setup has been reduced to a minimal-unit (or single-pore) configuration, as done by El-Rahman et al. (2017), to reduce the computational cost and ensure the maximum possible numerical resolution in the direction of shock propagation for a given number of discretization points, or degrees of (numerical) freedom. In spite of this choice, full resolution of the propagating shocks was still not attainable with the available resources.

The paper is organized as follows. Details of the computational setup are presented in § 2, followed by a grid convergence study. In § 3, the various regimes of thermoacoustic amplification and energy cascade are identified. In § 4, the dependence of instability growth rates on acoustic phasing is presented utilizing the acoustic energy budgets. A time-domain nonlinear acoustic model is developed in § 5, highlighting the importance of accounting for thermodynamic nonlinearities to match the results from the Navier-Stokes simulations presented in § 6. Finally, in § 7, the distribution of the spectral energy density at steady state is discussed and modelled heuristically based on dimensional analysis statements, establishing the existence of an equilibrium thermoacoustic energy cascade analogous to hydrodynamic turbulent cascade.

2 Problem Formulation

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

The proposed computational setup (figure 1, top) is a straight, two-dimensional, axially periodic minimal-unit (or single-pore) thermoacoustic device composed of four constant-area sections ($a$, $b$, $c$, and $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. (1998). 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., 2005), thus introducing further complications in the proposed canonical setup (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_{\mathrm{ref}}$, $\rho_{C}=\rho_{\mathrm{ref}}$, $T_{C}=T_{\mathrm{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 setup (of the order of $1$ m) compared to the typical shock thickness scale (of the order of $100\,\mu$m). Performing a fully resolved three-dimensional simulation of an equivalent experimental setup would be (even more so) unfeasible: for instance, the setup studied by Yazaki et al. (1998) 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 setup), and the latter are only relevant for very small ratios of the curvature radius to acoustic wavelength ($\sim 0.3$ in Yazaki’s setup, 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.

The total length of the device is fixed to $\ell_{a}+\ell_{b}+\ell_{c}+\ell_{d}=2.58~{}\mathrm{m}$, taken from the experimental setup of Yazaki et al. (1998). The height of the regenerator has been chosen such that $h_{b}\sim 2\delta_{k}$, where

[TABLE]

$\nu$ is the kinematic viscosity, $Pr$ is the Prandtl number, and $\omega$ 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. (2016), using the governing equations outlined in § 4.1. Further justification of the setup choice and a comparison of the linear numerical analysis with Yazaki et al. (1998)’s stability data are given in appendix A.

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

[TABLE]

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$, $\rho$, $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 setup, Gedeon streaming (Gedeon, 1995) 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 setup 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:

[TABLE]

where $\alpha=0.3$ is a relaxation coefficient, $\Delta 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, $\tau_{ij}$, and heat flux, $q_{i}$, are formulated based on the Stokes and the Fourier laws as

[TABLE]

respectively, where $S_{ij}$ is the strain-rate tensor, given by $S_{ij}=\frac{1}{2}\left(\partial u_{j}/\partial x_{i}+\partial u_{i}/\partial x_{j}\right)$, $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, $\mu(T)$, is varied with the temperature according to the Sutherland’s law, $\mu(T)=\mu_{\textrm{ref}}(T/T_{S,\textrm{ref}})^{1.5}(T_{S,\textrm{ref}}+S)/(T+S)$ where $S=120$ K is the Sutherland constant, $T_{S,\textrm{ref}}=291.15$ K and $\mu_{\textrm{ref}}=1.827\times 10^{-5}$  kg$\cdot$m${}^{-1}\cdot$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 $\gamma=1.4$, the reference density $\rho_{\textrm{ref}}=1.176$  kg$\cdot$m*-3*, the pressure $p_{\textrm{ref}}=101\,325$  Pa, the temperature $T_{\textrm{ref}}=300$  K, and the gas constant $R=p_{\textrm{ref}}/(\rho_{\textrm{ref}}\,T_{\textrm{ref}})$.

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 & Kolias, 1996; Kopriva, 1996; Sun et al., 2007; Jameson, 2010). 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 the Roe (1981)’s flux with the entropy correction by Harten & Hyman (1983). The numerical dissipation at the element interfaces scales as $\Delta_{h}^{N_{p}+1}$ where $\Delta_{h}$ is the characteristic length scale of neighbouring elements (Jameson & Lodato, 2014; Chapelier et al., 2016). 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, 2006; Lodato et al., 2016). The time integration is carried out explicitly with a 3rd order Runge-Kutta scheme and discretization order $p=3$ during the transient, whereas a 5th 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 et al., 2013, 2014) and unsteady shock-wavy wall interaction problems (Lodato et al., 2016, 2017).

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$ K temperature setting. The pressure waveform, expressed as a function of the travelling wave coordinate

[TABLE]

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_{\mathrm{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, $\Delta s$, as the sum of the point-to-point entropy increments:

[TABLE]

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

[TABLE]

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_{\Delta}\rightarrow 0}$ at continuum limit is estimated via a polynomial extrapolation,

[TABLE]

where $m_{0}$, $m_{1}$, $m_{2}$, and $m_{3}$ are fitting parameters and $r_{\Delta}$ is the grid refinement ratio,

[TABLE]

in which $\Omega$ denotes total area of the computational domain.

3 Regimes of Thermoacoustic Amplification

Three different regimes of thermoacoustic wave amplification can be identified by visual inspection of the pressure time series in figure 5. We attempt a rigorous classification here based on the dimensionless collapse of the nonlinear growth regime of the spectral energy density (figure 6), derived in more detail in § 7.2.

3.1 Spectral energy density

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

[TABLE]

where $\mathbf{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, 1975):

[TABLE]

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,

[TABLE]

into (11) and cycle averaging yields:

[TABLE]

where $\epsilon\sim\alpha/\omega\ll 1$ is the smallness parameter such that $t$ and $\epsilon t$ correspond to fast and slow time scales, respectively, $\overline{(\cdot)}$ denotes the cycle averaging operator defined as

[TABLE]

and $(\cdot)^{*}$ denotes the complex conjugate. Here $\omega_{0}$ is the angular frequency of the unstable second harmonic, and $E_{k}$ is the spectral energy density of the $k^{\mathrm{th}}$ mode. The pressure and velocity amplitudes of the $k^{\mathrm{th}}$ harmonic ($|\hat{p}_{k}|$ and $|\hat{u}_{k}|$, respectively) are functions of the $x$ coordinate and the slow time $\epsilon 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,\dots$) 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:

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 $\alpha 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. 2.

Hierarchical spectral broadening: Energy cascades down to higher harmonics hierarchically: the $k^{\mathrm{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), $\alpha$, that is

[TABLE]

The saturation of the spectral energy density $E_{k}$ occurs at $\alpha t\approx 10$ followed by the formation of resonating shock waves at limit cycle. 3. 3.

Shock dominated limit cycle: In this regime ($\alpha 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 $\alpha$ approximately as $\alpha^{3}$ (§ 7.2).

4 Harmonic growth analysis

In this section, the amplification of acoustic waves in the linear regime in the proposed minimal-unit model is discussed. A system-wide differential eigenvalue problem is formulated and solved numerically utilizing the strategy adopted by Lin et al. (2016) in (§ 4.1). Utilizing the eigenvalue analysis, the computational setup has been optimized (§ 2.1) and the acoustic energy budgets are derived (§ 4.2). This provides analytical expressions for the cycle-averaged thermoacoustic production and dissipation, and elucidates the role of the acoustic phasing on thermoacoustic instability. Finally, the effects of varying the hot-side temperature and the geometry on the thermoacoustic growth rates are discussed in § 4.3. While the eigenvalue analysis is restricted to the minimal-unit model, the physical conclusions and analysis presented in this section hold for any thermoacoustically unstable device operating in the low-acoustic amplitude regime.

4.1 Linear thermoviscous quasi-planar wave equations

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

[TABLE]

where primed variables $\left(\cdot\right)^{\prime}$ represent the fluctuations in the corresponding quantities and the subscript [math] denotes the base state. Axial diffusion terms and the $y$-momentum equation in (16) have been neglected based on the scaling analysis reported in Appendix B. Combining the cross-sectionally averaged (16a) and (16c) and cross-sectionally averaging (16b), accounting for isothermal and no-slip boundary conditions, yields

[TABLE]

respectively, where $h$ denotes the cross-sectional width of the duct, $U^{\prime}$ denotes the fluctuations in the flow rate,

[TABLE]

whereas $q^{\prime}$ and $\tau^{\prime}_{w}$ are the wall-heat flux and the wall-shear, respectively:

[TABLE]

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

[TABLE]

where $\sigma=\alpha+i\omega$ is the complex eigenvalue of the system with growth rate $\alpha$ and angular frequency $\omega$, leads to the thermoviscous set of quasi-planar wave equations in the frequency domain (Lin et al., 2016):

[TABLE]

where the thermoviscous functions $f_{\nu}$ and $f_{k}$ are given by

[TABLE]

4.2 Acoustic energy budgets for quasi-planar wave perturbations

Multiplying (17a) by $p^{\prime}/(\rho_{0}a^{2}_{0})$ and (17b) by $\rho_{0}U^{\prime}/h_{b}$ and adding them, yields the conservation equation

[TABLE]

for the one-dimensional acoustic energy density

[TABLE]

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

[TABLE]

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

[TABLE]

where $\mathcal{R}$ is the Rayleigh index and $\epsilon t$ is the slow time scale (cf. (13)). Relation (26) allows 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 $\tau^{\prime}_{w}$ terms (cf. (25)) in the respective duct sections. Utilizing the frequency domain linear equations (21), the wall-heat flux $\hat{q}$ and wall-shear $\hat{\tau}_{w}$ in the frequency domain are given by

[TABLE]

Combining (25), (27), and (21b), an analytical expression is obtained for the cycle-averaged production $\overline{\mathcal{P}}$ and dissipation $\overline{\mathcal{D}}$ of the acoustic energy density (24):

[TABLE]

In the above relations, $\Phi^{T}_{\mathcal{P}}$ and $-\Phi^{S}_{\mathcal{P}}$ weigh the contributions to the thermoacoustic energy production by the travelling-wave, $\mathrm{Re}(\hat{p}^{*}\hat{U})$, and the standing-wave, $\mathrm{Im}(\hat{p}^{*}\,\hat{U})$, components respectively. Their expressions read:

[TABLE]

For $h_{b}/2\delta_{k}\leq 1.13$, $|\Phi^{T}_{\mathcal{P}}|>|\Phi^{S}_{\mathcal{P}}|$, which implies that the regenerator half-width $h_{b}$ must remain comparable to or smaller than the Stokes boundary layer thickness $\delta_{k}$ to achieve higher thermoacoustic amplification of travelling waves ($\psi_{Up}\simeq 0^{\circ}$) (figures 7$a$ and 8). However, to maximize thermoacoustic energy production for standing waves ($\psi_{Up}\simeq\pm 90^{\circ}$), a larger regenerator half-width ($h_{b}/2\delta_{k}>1.13$) is required. While production alone for a purely travelling wave ($\psi_{Up}=0$) is maximized in the limit $h_{b}/2\delta_{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\delta_{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

[TABLE]

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

[TABLE]

respectively, are also analysed, where

[TABLE]

is a reference viscous time scale in the regenerator with $\nu_{0}$ evaluated at $T_{C}=300\,$K.

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

[TABLE]

and more rapidly than the dissipation $\overline{\mathcal{D}}_{*}$, yielding positive values of the Rayleigh index $\mathcal{R}_{*}$ (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 (figure 9$a$, 10$a$) as

[TABLE]

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*}$ (figure 10$b$) peaks at $h_{d}/h_{c}\approx 1.28$. Moreover, the high mechanical impedance of the section $d$ ($\rho_{0}a_{0}h_{d}$) for high values of $h_{d}$ increases the standing wave component of the acoustic power $\mathcal{R}_{S*}$ and decreases the travelling wave component $\mathcal{R}_{T*}$. The variation of the transverse geometrical parameters does not significantly alter the frequency. Therefore, changes in the ratio $(f_{k}-f_{\nu})/(1-f_{\nu})$ are also negligible and do not influence the growth rates significantly (cf. (29)).

5 Formulation of a Nonlinear Thermoacoustic Model

As a result of the modal thermoacoustic instability, large pressure amplitudes ($\sim 160~{}\mathrm{dB}$) are generated, which result in the nonlinear steepening of the waveform. In the spectral space, the nonlinear steepening can be viewed as the cascade of energy from the unstable mode into higher harmonics with correspondingly shorter wavelengths. Moreover, inside the regenerator, large amplitude perturbations in thermodynamic quantities are responsible for thermoacoustic nonlinearities. As a result of nonlinear wave propagation, thermoacoustically sustained shock waves propagate in the system. While the analysis above highlights that quasi-travelling wave phasing is essential for high thermoacoustic growth rates, nonlinear steepening is also favoured by such phasing (Biwa et al., 2014).

In this section, a first-principles-based theoretical framework accounting for acoustic and thermoacoustic nonlinearities up to second order is developed, and a quasi one-dimensional evolution equation is obtained for nonlinear thermoacoustic waves (49). In § 5.1 mass, momentum, and energy (combined with the second law of thermodynamics) equations correct up to second order are introduced. Furthermore, in § 5.2 and § 5.4, cross-sectionally averaged nonlinear spatio-temporal evolution model equations are derived with the time-domain approximations of wall-shear and wall-heat flux outlined in § 5.3.

5.1 Governing equations for nonlinear thermoviscous perturbations

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

[TABLE]

where, again, primed variables $\left(\cdot\right)^{\prime}$ represent the fluctuations in corresponding quantities, the subscript [math] denotes the base state, whereas $\xi_{B}=2/3$ is the ratio of the bulk viscosity coefficient $\mu_{B}$ to the shear viscosity coefficient $\mu$. The terms on the left-hand side are linear in the perturbation variables while those on the right-hand side are nonlinear. The entropy generation due to viscous dissipation is neglected, as well are the pressure gradients and velocity in $y$ direction (boundary layer assumption), and the fluctuations in the diffusivity coefficients $\mu$ and $k$ (see appendix B). Higher harmonics with correspondingly shorter wavelengths are generated due to the nonlinear spectral energy cascade. Consequently, the axial diffusion terms in (35b) and (35c), which have been neglected in § 4, become significant and act as the primary sink of energy at large harmonic scales in the spectral space.

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

[TABLE]

where the thermodynamic coefficients $\alpha$ and $\beta$ are given by

[TABLE]

The coefficients $\alpha_{p}$ and $\alpha_{s}$ contribute to first order wave propagation and thermoacoustic effects while the second-order coefficients $\beta_{s},~{}\beta_{p},$ and $\beta_{sp}$ in (36) account for the corresponding nonlinear effects. Hedberg & Rudenko (2011) have demonstrated the hysteretic effects of nonlinear wave propagation retaining only the $\beta_{s}$ nonlinear term and $\alpha_{p}$ term accounting for irreversible entropy changes. Figure 11 shows the hysteresis cycle for the second order constitutive relation in (36) 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 et al., 1998). Consequently, the constitutive relation (36) can be approximated with the terms corresponding to $\alpha_{s},\alpha_{p}$, and $\beta_{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 (36). 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 setup (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 (35b) and $T^{\prime}$ in (35c) can be neglected and entropy fluctuations remain second order in the nonlinear regime as well (Hamilton et al., 1998). Up to second order, the wave propagation in the duct sections $a$, $c$, and $d$ is governed by

[TABLE]

where $h=h_{a},~{}h_{c},$ or $h_{d}$. Second-order nonlinearities in (38) 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\delta_{\nu}\sim 1$, $h_{b}/2\delta_{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, 1958; Pierce, 1989), the following decomposition for the entropy and velocity fluctuations inside the regenerator is considered:

[TABLE]

where $u^{\prime}_{\nu}$ 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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

Equations (45) and (46) 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 (44a) for $j\rightarrow\infty$ are discussed.

5.4 Nonlinear thermoviscous wave equations

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

[TABLE]

where

[TABLE]

Equation (38) 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, (35a), (35c), and (36) are combined to obtain

[TABLE]

where

[TABLE]

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

[TABLE]

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 (35b) and (35c) 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 (49) and expressions for the terms $q_{2}$, $\mathbb{D}_{s}$, $\mathbb{T}$, and $\mathbb{C}$ are given in Appendix C. Equations (45), (46), (47), and (49) constitute the governing equations for the spatio-temporal evolution of large amplitude acoustic perturbations inside the regenerator. The macrosonic thermoacoustic interaction (51) 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 ($\rho^{\prime}>0$), compared to dilatations ($\rho^{\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:

[TABLE]

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

[TABLE]

where $\tau_{j}=1/\nu_{0}\zeta^{2}_{j}$ defines the viscous relaxation time for the $j^{\mathrm{th}}$ viscous mode. Hence, writing the relaxation functional for some function $\phi(x,t)$, namely,

[TABLE]

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

[TABLE]

Equations (56) and (57) 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 (47) and (52), complete the nonlinear wave propagation model equations. However, in the present work, (45) and (46) have been considered for time integration for simplicity.

Sugimoto (2010) 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 equation (55) 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 (2016) 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 (56) and (56) 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 (38), (39), (45), (46), (47), and (49) are integrated in time utilizing a $4^{\mathrm{th}}$ 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 ($\Delta p_{\mathrm{ml}}$) due to nonlinear vortex generation which results in second order losses. Hence, minor losses can be modeled as quadratic functions of the first order interface flow rate, i.e., $\Delta p_{\mathrm{ml}}\sim U^{2}_{i}$. In analogy with the approach followed by Lin et al. (2016), 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:

[TABLE]

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

[TABLE]

Hence, the second order flow rate $U^{(2)}_{i}$ is calculated imposing the condition in (58). 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.

6 Nonlinear Spectral Energy Dynamics

In order to elucidate the physics of the hierarchical spectral broadening regime, the results from the nonlinear model derived in the previous section are discussed here and compared to the high-order Navier-Stokes calculations. Pressure time series from the time integration of (52) are in fairly good agreement with the fully compressible Navier-Stokes simulations (see figure 13). Since several nonlinearities are neglected in calculations via (52), time integration results in spurious temporal variations of the time averaged pressure fluctuations which are removed in further discussions. The time integration of the nonlinear model up to the limit cycle offers a significant reduction in computational cost (about $500$ times faster than the fully compressible Navier-Stokes simulations) and predicts the captured limit cycle amplitudes within $80\%$ accuracy. Additionally, upon excluding the macrosonic thermoacoustic interaction term from the model (51), the accuracy of the predicted limit cycle amplitude gets reduced to $60\%$.

Figure 14 shows the time evolution of the spectral energy density $E_{k}$ of the unstable mode (cf. (13)) and its first seven overtones, as obtained from the nonlinear model discussed in § 5, and compares it to the results obtained from the fully compressible Navier-Stokes simulations for the signal shown in figure 13. In the spectral broadening regime of thermoacoustic wave amplification, the growth of spectral energy density of the $k^{\mathrm{th}}$ harmonic obtained from linear interpolation is approximately $k\alpha_{2}/2$ (table 3), where $\alpha_{2}=2\alpha$ (since $E\propto p^{\prime 2}$) is the growth rate of the spectral energy density of the unstable mode, i.e., $E_{k}\sim e^{k\alpha_{2}t/2}$. For instance, figure 14 corresponds to the case $T_{H}=450$ K, for which the growth rate is $\alpha=8.64\,\mathrm{s}^{-1}$; therefore, the growth rate of the spectral energy density is $\alpha_{2}=2\alpha=17.28\,\mathrm{s}^{-1}$ for the unstable mode.

Nonlinear energy cascade in the spectral space can be further explained using reduced order modelling. Assuming propagation of purely travelling waves in the system and eliminating $U^{\prime}$ from equations (38) and (39), the following Burgers equation for the pressure fluctuations $p^{\prime}$ is obtained:

[TABLE]

where $\xi=a_{0}t-x$ is the travelling wave coordinate and $\delta$ is the axial dissipation coefficient given by

[TABLE]

In general, the wall-shear and the wall-heat flux expressions derived in (56) and (57) can be used as forcing functions in the above Burgers equation. However, the abrupt area changes present in the setup under study make the generalization and further time domain simplification seemingly challenging. Hence, we seek to only qualitatively explain the temporal evolution of the nonlinear cascade utilizing the Burgers equation. Introducing thermoacoustic amplification by adding a simple linear forcing term in (60) yields:

[TABLE]

where $\alpha_{\mathrm{th}}$ accounts for the thermoacoustic growth rate. Substituting a Fourier expansion for acoustic pressure $p^{\prime}$, namely,

[TABLE]

the following modal evolution equation is obtained:

[TABLE]

where

[TABLE]

$\alpha_{\mathrm{th}}$ determines the rate of thermoacoustic amplification of the $k^{\mathrm{th}}$ mode ($\alpha_{\mathrm{th}}=8.64\,\mathrm{s}^{-1}$ for $k=2$), and $\omega_{0}$ is the angular frequency of the unstable mode.

$Q(p_{k})$ is the nonlinear cascade function quantifying the scale-by-scale flux of energy in the spectral space from the unstable mode to the harmonics which are dissipated by the molecular diffusion effects (momentum and thermal diffusivity). The third term on the right-hand side of (64) signifies dissipation of the $k^{\mathrm{th}}$ mode through molecular diffusion. Assuming that the subsequent overtones of the unstable mode are characterized by pressure amplitudes which are an order of magnitude smaller ($p_{k+2}/p_{k}\ll 1$), the modal growth rate (64) can be approximated such that

[TABLE]

and so on for higher harmonics. Equation (66) shows that the growth rates of the overtones of the unstable mode due to the energy cascade are proportional to the ratio of the oscillation frequencies, that is to say, $p_{k}\sim e^{k\alpha t/2}$ in the hierarchical spectral broadening regime. It is important to note that the growth of higher overtones is not strictly exponential since they are not subject to thermoacoustic instability. Energy is cascaded into the overtones of the thermoacoustically unstable mode via nonlinear energy cascade due to the high acoustic wave amplitude. Nonetheless, it is possible to assume modal growth for a small time interval in the spectral broadening regime to quantify the exponential growth rate for each harmonic (table 3).

7 Thermoacoustically sustained shock waves

The continued nonlinear energy cascade results in thermoacoustically sustained propagating shock waves at the limit cycle. While shock wave propagation outside the regenerator is quasi-planar, inside the regenerator, shock wave propagation is significantly affected by the wall-shear stresses and the wall-heat flux. In this section, time domain nonlinear model results are used to quantify the production of the energy density at the limit cycle (§ 7.1). Furthermore, utilizing the Navier-Stokes simulations and the nonlinear model results, a spectral scaling for the spectral energy density is derived; such a scaling, in turn, suggests the existence of thermoacoustic energy cascade analogous to turbulent energy cascade (§ 7.2).

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. Figures 15$a-d$ show the phasing of $p^{\prime}$ and $u^{\prime}$ extracted at $x=20$ mm, $y=0$ mm and compare the evolution of the relative phase angle $\psi_{up}$ obtained from the time domain nonlinear model and the fully compressible Navier-Stokes simulations using the correlation function expression

[TABLE]

In general, the above relation can be used to study the evolution of the phase difference between close to harmonically oscillating quantities (Balasubramanian & Sujith, 2008). The relative phase $\psi_{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 (51) 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 (47) and (52) such that

[TABLE]

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

Furthermore, figure 15$f$ shows that the phase angle $\psi_{pq}$ saturates at a value $\psi_{pq}\approx 33.57^{\circ}$ compared to $\psi_{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 $\psi_{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, 2000; Nazarenko, 2011). Assuming travelling wave propagation, the total energy density of the $k^{\mathrm{th}}$ harmonic associated to planar wave propagation, $E^{(1\mathrm{D})}_{k}$, can be defined as

[TABLE]

where $E$ is the squared $L_{2}$ norm defined in (11) where $p^{\prime}$ and $u^{\prime}$ are defined as a function of the travelling wave coordinate $\xi$. 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 $\delta$ in (61). Thus, $\mathcal{E}$ is purely governed by thermoacoustic wave amplification at smaller harmonics and viscous dissipation at higher harmonics and scales as

[TABLE]

where $\alpha_{\mathrm{eff}}$ is the effective energy amplification rate at the limit cycle.

Based on the macrosonic thermoacoustic interaction, $\alpha_{\mathrm{eff}}$ can be estimated as

[TABLE]

where $\alpha$ 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^{(1\mathrm{D})}_{k}/\rho_{0}$ only depends on the rate of energy cascade $\mathcal{E}$ and on the angular frequency of the harmonic $\omega_{k}$, the following scaling is obtained:

[TABLE]

Substituting into (69) and eliminating $\rho_{0}$, the scaling for the spectral energy density $E_{k}$ remains the same such that

[TABLE]

The scaling derived in (73) shows that the distribution of the energy density at limit cycle in the spectral space decays as $\omega^{-5/2}_{k}$—where $\omega_{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 $\alpha^{3}_{\mathrm{eff}}$ (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 $\alpha^{3}_{\mathrm{eff}}$) 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 et al., 2017), 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 (73) can be further improved accounting for propagation of very large harmonics inside the regenerator.

8 Conclusions

The linear and nonlinear regimes of thermoacoustic wave amplification have been modelled up to the formation of shock waves in a minimal unit looped resonator with the support of high-fidelity fully compressible Navier-Stokes simulations. The computational setup is inspired by the experimental investigations conducted by Yazaki et al. (1998) and geometrically optimized to maximize growth rates for the quasi-travelling wave mode. Three regimes of thermoacoustic wave amplification have been identified: ($i$) a monochromatic or modal growth regime, ($ii$) a hierarchical spectral broadening or nonlinear growth regime and ($iii$) a shock-dominated limit cycle. The modal growth regime is characterized by exponential amplification of thermoacoustically unstable modes. An acoustic energy budget formulation yielding a closed form expression of the Rayleigh index has been developed and the effect of variations in geometry and hot-to-cold temperature ratios on the thermoacoustic growth rates have been elucidated. The limit cycle regime exhibits many features of Kolmogorov’s equilibrium turbulence, where energy, steadily injected at the integral length scale (wavelength of the second-mode harmonic), cascades towards higher wave numbers via inviscid mechanisms (wave steepening) and is finally dissipated at the Kolmogorov’s length scale (of the order of the shock thickness). A grid sensitivity analysis has been carried out at the limit cycle to ensure that the entropy jump across the captured shock waves is grid convergent, hence assuring the same fidelity typically attributed to direct numerical simulations of turbulent flows, with the caveat that shocks in the present study are not fully resolved by the computational mesh.

The existence of an equilibrium thermoacoustic energy cascade has thus been shown. The spectral energy density at the limit cycle, in particular, has been found to decay as $\omega^{-5/2}$ in spectral space, the relevant intensity scaling with growth rate $\alpha$ as $\alpha^{3}$. Such findings are confirmed by a novel theoretical framework to model thermoacoustic nonlinearities, which has lead to the formulation of a one-dimensional time-domain nonlinear acoustic model. The model is correct up to second order in the perturbation variables and addresses the fundamental problem of high amplitude wave propagation in the presence of wall-shear and wall-heat flux, accounting for thermodynamic nonlinearities such as the second-order interactions between the pressure fluctuations and the wall-heat flux, namely the macrosonic thermoacoustic interaction term. The model also confirms the dynamics of energy transfer across scales in the nonlinear growth regime: the growth of higher harmonics is hierarchical in nature and higher harmonics are amplified at faster rates. in particular the $k^{\mathrm{th}}$ harmonic grows at the rate of $\alpha k/2$, where $\alpha$ is the rate of growth of energy in the unstable mode.

Acknowledgments

Prateek Gupta acknowledges the support of the Lynn Fellowship at Purdue University, Ford Motor Company and Rolls-Royce Corporation (Indianapolis). The use of the spectral difference solver originally developed by Antony Jameson’s group at Stanford University, and financial support of Guido Lodato by CNRS under the INSA Turbulence and Simulation Chair are also gratefully acknowledged. Computations have been run on the high-performance computing resources provided by the Rosen Center for Advanced Computing (RCAC) at Purdue University and by the Institut du Développement et des Ressources en Informatique Scientifique (IDRIS-CNRS) under the allocation i2015-2a7361 and of the Haute Normandie Computing center CRIANN.

Appendix A Experimental Validation of linear model

Figure 17$a$ shows the comparison between the neutral stability curve evaluated from the linear eigenvalue analysis discussed in § 4, experimental data obtained by Yazaki et al. (1998) (setup referenced as Y-1998), and the numerical predictions obtained by Guedra & Penelet (2012). The lengths of smaller cross-section ducts ($\ell_{a}$ and $\ell_{c}$ in figure 1) have been estimated from the laser Doppler velocimetry (LDV) data reported by Yazaki et al. (1998) ($\ell_{a}=0.153$ m and $\ell_{c}=0.5$ m). The rest of the missing geometrical details have been taken from Guedra & Penelet (2012). The looped geometry of  Biwa et al. (2011) (setup referenced as B-2011), is shown to be thermoacoustically unstable at the reported critical $T_{H}$ value, ($T_{cr}/T_{C}=1.52$) only if a mean negative temperature gradient from $\ell_{b}<x<\ell_{\mathrm{buffer}}+\ell_{b}$, is assumed in the linear analysis. Moreover, the unstable mode is a quasi-travelling wave ($-90^{\circ}<\psi_{Up}<90^{\circ}$). However, by considering an abrupt temperature change from $T_{H}$ to ambient temperature $T_{C}$, the ratio $T_{cr}/T_{C}$ is very large ($\sim$2.3) and the unstable mode is a standing wave ($\psi_{Up}\simeq\pm 90^{\circ}$). Since the nonlinear cascade leading to shock formation is inhibited by the standing wave phasing compared to the quasi-travelling wave phasing (Biwa et al., 2014), a geometry similar to Y-1998 has been chosen in the current study which exhibits quasi-travelling wave unstable mode due to a wider duct section.

Appendix B Scaling of governing equations

The governing equations for compressible flows in two dimensions read:

[TABLE]

where $\Phi$ accounts for entropy generation due to viscous gradients. These equations are written in terms of perturbations, denoted by primed symbols, in the relevant variables $u$, $v$, $\rho$, $p$, $T$, and $~{}s$. It shall be noted that, for a nonlinear acoustic field of velocity amplitude scale $\mathcal{U}$, the acoustic Mach number $\mathcal{M}=\mathcal{U}/a_{0}$ is $\mathcal{O}(10^{-1})$ and the aspect ratio of the regenerator $h_{b}/l_{b}$ is $\mathcal{O}(10^{-2})$. As a result, the perturbations in the $y$ component of velocity (namely, $v^{\prime}/u^{\prime}\sim h_{b}/l_{b}$) from the momentum and entropy equations can be neglected and the following equations, including second order terms in $\mathcal{M}$, are obtained:

[TABLE]

In the above equations, the subscript [math] denotes the base state. Also, the terms on the left-hand side of the equations are linear in the perturbation variables, whereas the terms on the right-hand side are non-linear. Assuming $\omega^{-1}$ as the characteristic time scale of the acoustic field, the viscous dissipation terms on the left-hand side of (75b) and (75c), and the entropy generation due to viscous gradients $\Phi$ scale relatively, as

[TABLE]

where the following scaling relations have been used:

[TABLE]

Due to the nonlinear energy cascade, higher harmonics are generated, as a result of which, the characteristic time scale of the acoustic field $\omega^{-1}$ decreases. Consequently, the characteristic Stokes layer thickness, $\delta_{\nu}$, decreases and the ratio $h_{b}\omega/a_{0}$ increases. Accordingly, the terms which scale as $\mathcal{O}((\delta_{\nu}/h_{b})^{2}\mathcal{M})$ can be dropped and terms corresponding to the local perturbations in viscosity and thermal conductivity, as well as the entropy generation term due to viscous gradients, are neglected. Finally, the equations governing the spatio-temporal evolution of the acoustic field, correct up to the second order, are given as

[TABLE]

Appendix C Time domain model and convergence

Integrating (79a) and (79c) in $y$, the following equations are obtained:

[TABLE]

[TABLE]

where,

[TABLE]

Also, the density constitutive equation, up to second order, is given by

[TABLE]

and the entropy perturbations $s^{\prime(1)}$, up to first order, are

[TABLE]

Notice that the superscripts $(1)$ and $(2)$ are dropped hereafter for convenience of notation. Integrating equation (84) along $y$, differentiating in time, and combining with (80) and (81) to eliminate $\tilde{\rho}$, the following equation is obtained:

[TABLE]

where

[TABLE]

and,

[TABLE]

The above equations account for thermodynamic, as well as, convective nonlinearities, with $\mathbb{D}_{s}$ denoting the axial conduction term. The treatment of quadratic nonlinearities gets significantly simplified due to the Fourier expansions of the viscous and entropic modes. For simplicity, only the nonlinear macrosonic thermoacoustic interaction $\mathbb{Q}$ is retained in the present one-dimensional computations. However, the complete model equations (86)–(88) should be considered in the case of relatively large regenerators.

In order to show the convergence of the viscous and entropic modes, $\check{u}_{j}$ and $\check{s}_{j}$, respectively, (45) and (46) are pre-multiplied by $(-1)^{j+1}\zeta_{j}$ and summed over $j$ to obtain:

[TABLE]

For the wall-shear $\tau^{\prime}_{w}$ and the wall-heat flux $q^{\prime}$ to converge, (89) and (90) must yield the same values of $q^{\prime}$ and $\tau^{\prime}_{w}$ for $N$ and $N+1$ in the limit of $N\rightarrow\infty$. Hence, the following conditions are obtained:

[TABLE]

The above relations show that, in order for convergence to be ensured, the magnitudes of the viscous modes $\check{u}_{j}$ and $\check{s}_{j}$ must decay as $\zeta^{-3}_{j}$ for large values of the index $j$. The relaxation functional used in (53) and (54) for very large $j$ yields:

[TABLE]

for some function $\phi(x,t)$. After substituting the above approximation in (53) and (54), it is shown that the magnitudes of the viscous modes $\check{u}_{j}$ and $\check{s}_{j}$ decay as $\zeta^{-3}_{j}$ in the limit of $j$ such that the approximation in (92) holds valid, which, in turn, guarantees the convergence of the infinite series approximation for the viscous and entropic modes.

Appendix D Windowed shock capturing

At the limit cycle, the perturbation fields exhibit shock wave propagation causing the formation of large gradients at the limit of grid resolution. Spatially windowed Legendre polynomial expansions are therefore used in order to add an artificial viscosity term to ensure a proper resolution of longitudinal gradients in the model equations. The implemented strategy, which is briefly summarized below, is similar to the one proposed by Persson & Peraire (2006) for discontinuous Galerkin methods. For a spatial window of size $N_{p}$, the local pressure field is reconstructed using Legendre polynomial basis expansions up to orders $N_{p}$ and $N_{p}-1$, such that

[TABLE]

where $\psi_{n}$ denotes the Legendre polynomial of order $n$ and $p_{n}$ is the relevant $n^{\mathrm{th}}$ mode. A modal smoothness indicator $S_{e}$ is hence constructed from the ratio of the inner products of $p^{\prime}-\tilde{p}^{\prime}$ and $p^{\prime}$, such as to detect the onset of excessively high high-frequency modes, which is typical of insufficiently unresolved signals:

[TABLE]

The sensor is then used to trigger the shock capturing artificial diffusivity $\epsilon_{e}$, which is evaluated as

[TABLE]

where $s_{e}=\log S_{e}$ and the corresponding parameters are $s_{0}=3.2$, $\kappa=0.9(\mathrm{max}(s_{e})-s_{0})$, and $\epsilon_{0}=100\Delta x/N_{p}$, $\Delta x$ being the grid spacing of the discretized one-dimensional domain.

Finally, the artificial viscosity augmented equations to capture the nonlinear acoustic fields become:

[TABLE]