parallelized Fortran 90 code employing the spectral difference local spatial
reconstruction for hexahedral elements on unstructured grids
kopriva1996conservative ; kopriva1996conservative2 ; sun2007high ;
Jameson_JSC_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 roe1981approximate ’s flux with the entropy correction by
harten1983self . 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 lodato:14b ; ChapelierLJ_CF_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:06 ; lodato:16 . 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. The same solver
has been used and validated in a wide variety of flow configurations,
including turbulent channel flow lodato2013discrete ; lodato2014structural and
unsteady shock-wavy wall interaction problems lodato:16 ; lodato:17 .

Results shown at shock-dominated limit cycle hereafter correspond to the
finest computational grid simulated. The reader is referred to
GuptaLS_JFM_2017 for the detailed grid-sensitivity analysis.

###  III Regimes of Thermoacoustic Amplification

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

Figure 4.4: Time series of pressure fluctuations at $x=1.54$ m (figure 4.1)
for $T_{H}=$450 K with insets showing regimes ($i$) modal growth, ($ii$)
hierarchical spectral broadening, and ($iii$) the shock-dominated limit cycle.

####  III.1 Spectral energy density

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

$\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^{2}_{0}}\right)^{\mathrm{T}},$
(4.10)

where $\mathbf{X}$ is the state vector containing the dimensionless
perturbation variables. Throughout, we utilize the isentropic acoustic energy
norm (see Section II), which is the squared $L_{2}$ norm, as the perturbation
energy density Naugol1975spectrum :

$E=\frac{1}{2}\rho_{0}a^{2}_{0}\mathbf{X}^{\mathrm{T}}\mathbf{X}=\frac{1}{2}\rho_{0}u^{\prime
2}+\frac{p^{\prime 2}}{2\rho_{0}a^{2}_{0}}.$ (4.11)

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,

$\displaystyle u^{\prime}(x,t)=\sum_{\begin{subarray}{c}k=-\infty\\\ k\neq
0\end{subarray}}^{+\infty}\hat{u}_{k}(x,\epsilon
t)\,e^{i\frac{k}{2}\omega_{0}t},$ $\displaystyle\quad
p^{\prime}(x,t)=\sum_{\begin{subarray}{c}k=-\infty\\\ k\neq
0\end{subarray}}^{+\infty}\hat{p}_{k}(x,\epsilon
t)e^{i\frac{k}{2}\omega_{0}t},$ (4.12)
$\displaystyle\hat{u}_{-k}=\hat{u}^{*}_{k},$
$\displaystyle\quad\hat{p}_{-k}=\hat{p}^{*}_{k},$

into Eq. (4.11) and cycle averaging yields:

$\overline{E}=2\sum_{k=1}^{\infty}E_{k},\quad
E_{k}=\frac{1}{2}\rho_{0}|\hat{u}_{k}|^{2}+\frac{|\hat{p}_{k}|^{2}}{2\rho_{0}a^{2}_{0}},$
(4.13)

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

$\overline{(\cdot)}=\frac{1}{T_{0}}\int_{\epsilon t}^{\epsilon
t+T_{0}}(\cdot)dt,\quad T_{0}=\frac{2\pi}{\omega_{0}},$ (4.14)

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 4.4. 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.

####  III.2 Regime classification

Figure 4.5: Evolution of dimensionless spectral energy density of the unstable
mode $k=2$ and its first seven overtones, scaled by the second mode angular
frequency $\omega_{0}$, base speed of sound $a_{0}$, and rate of spectral
energy transfer $\mathcal{E}$ Eq. (4.66). (–), $T_{H}=450~{}\mathrm{K}$; (-
-), $T_{H}=500~{}\mathrm{K}$; ($\cdot\cdot\cdot$), $T_{H}=550~{}\mathrm{K}$.

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

  1. 1.

Modal growth: Only the thermoacoustically unstable mode amplifies
exponentially and all the other modes excited by the initial perturbation
field decay. Higher harmonics begin to grow after $\alpha t\approx 3.2$
(figure 4.5), 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 Section V and
Section VI), $\alpha$, that is

$\alpha_{k}=\alpha k/2,\quad\,k\in\\{4,6,8,\dots\\}.$ (4.15)

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}$
(Section VII).

The harmonic growth analysis corresponds to the linear stability analysis of
the resonator by formulation of a system-wide differential eigenvalue problem.
The reader is referred to GuptaLS_JFM_2017 for further details on eigenvalue
analysis and thermoacoustic energy budget analysis.

###  IV 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 LinSH_JFM_2016 in (Section IV.1). Utilizing the eigenvalue
analysis, the computational setup has been optimized (Section II.1) and the
acoustic energy budgets are derived (Section IV.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 Section IV.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.

####  IV.1 Linear thermoviscous quasi-planar wave equations

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

$\displaystyle\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,$ (4.16a) $\displaystyle\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,$ (4.16b)
$\displaystyle\frac{\partial s^{\prime}}{\partial
t}+u^{\prime}\frac{ds_{0}}{dx}-\frac{k_{0}}{\rho_{0}T_{0}}\frac{\partial^{2}T^{\prime}}{\partial
y^{2}}=0,$ (4.16c)

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

$\displaystyle\frac{\partial p^{\prime}}{\partial
t}+\frac{\rho_{0}a^{2}_{0}}{h}\frac{\partial U^{\prime}}{\partial
x}=\frac{\rho_{0}a^{2}_{0}}{h}\frac{q^{\prime}}{\rho_{0}C_{p}T_{0}},$ (4.17a)
$\displaystyle\frac{1}{h}\frac{\partial U^{\prime}}{\partial
t}+\frac{1}{\rho_{0}}\frac{\partial p^{\prime}}{\partial
x}=\frac{1}{h}\frac{\tau^{\prime}_{w}}{\rho_{0}},$ (4.17b)

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

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

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

$q^{\prime}=2k_{0}\left.\frac{\partial T^{\prime}}{\partial
y}\right|_{y=+h/2},\quad\tau^{\prime}_{w}=2\mu_{0}\left.\frac{\partial
u^{\prime}}{\partial y}\right|_{y=+h/2}.$ (4.19)

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 Eq. (4.17),
namely,

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

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
LinSH_JFM_2016 :

$\displaystyle\sigma\hat{p}=\frac{\rho_{0}a^{2}_{0}}{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{dT_{0}}{dx},$
(4.21a)
$\displaystyle\sigma\hat{U}=-\frac{h_{b}}{\rho_{0}}\left(1-f_{\nu}\right)\frac{d\hat{p}}{dx},$
(4.21b)

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

$f_{\nu}=\frac{\tanh(\eta h_{b}/2)}{\eta h_{b}/2},\quad f_{k}=\frac{\tanh(\eta
h_{b}\sqrt{Pr}/2)}{\eta
h_{b}\sqrt{Pr}/2},~{}~{}\mathrm{with}~{}~{}\eta=\sqrt{i\omega/\nu_{0}}.$
(4.22)

####  IV.2 Acoustic energy budgets for quasi-planar wave perturbations

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

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

for the one-dimensional acoustic energy density

$E=\frac{1}{2}\frac{p^{\prime
2}}{\rho_{0}a^{2}_{0}}+\frac{1}{2}\rho_{0}\left(\frac{U^{\prime}}{h_{b}}\right)^{2},$
(4.24)

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

$\mathcal{I}=\frac{p^{\prime}U^{\prime}}{h_{b}},\quad\mathcal{P}-\mathcal{D}=\frac{p^{\prime}q^{\prime}}{h_{b}\rho_{0}C_{p}T_{0}}+\frac{\tau^{\prime}_{w}U^{\prime}}{h_{b}^{2}}.$
(4.25)

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

$\frac{d}{d(\epsilon t)}\int_{L}\overline{E}(x,\epsilon
t)dx=\int_{L}\left(\overline{\mathcal{P}}-\overline{\mathcal{D}}\right)dx=\mathcal{R},$
(4.26)

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

$\displaystyle\hat{\tau}_{w}=h_{b}\frac{\partial\hat{p}}{\partial x}f_{\nu},$
(4.27a)
$\displaystyle\hat{q}=h_{b}(i\omega)C_{p}T_{0}\left(\frac{\Theta}{\left(1-Pr\right)\omega^{2}}\frac{\partial\hat{p}}{\partial
x}\left(f_{k}-f_{\nu}\right)-\frac{\gamma-1}{a^{2}_{0}}\hat{p}f_{k}\right).$
(4.27b)

$(a)$$(b)$$(c)$

Figure 4.6: Variation of thermoviscous functionals affecting the cycle
averaged thermoacoustic production $\overline{\mathcal{P}}$ ($a$), and
dissipation $\overline{\mathcal{D}}$ ($b$ and $c$) of acoustic energy density
versus the ratio of the regenerator half-width $h_{b}/2$ to the Stokes
boundary layer thickness $\delta_{k}$ (cf. Eq. (4.28)). ($a$): (–),
$\Phi^{T}_{\mathcal{P}}$; (- -), $-\Phi^{S}_{\mathcal{P}}$.

Combining Eq. (4.25), Eq. (4.27), and Eq. (4.21b), an analytical expression is
obtained for the cycle-averaged production $\overline{\mathcal{P}}$ and
dissipation $\overline{\mathcal{D}}$ of the acoustic energy density Eq.
(4.24):

$\displaystyle\overline{\mathcal{P}}=\frac{\Theta}{2\left(1-Pr\right)h_{b}}\left[\Phi^{T}_{\mathcal{P}}\;\mathrm{Re}\left(\hat{p}^{*}\hat{U}\right)-\Phi^{S}_{\mathcal{P}}\;\mathrm{Im}\left(\hat{p}^{*}\;\hat{U}\right)\right],$
(4.28a)
$\displaystyle\overline{\mathcal{D}}=\mathrm{Re}\left(if_{k}\right)\frac{\omega\left(\gamma-1\right)}{2\rho_{0}a^{2}_{0}}|\hat{p}|^{2}-\mathrm{Im}\left(\frac{f_{\nu}}{1-f_{\nu}}\right)\,\frac{\rho_{0}\omega}{2h^{2}_{b}}|\hat{U}|^{2}.$
(4.28b)

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:

$\Phi^{T}_{\mathcal{P}}=\mathrm{Re}\left(\frac{f_{k}-f_{\nu}}{1-f_{\nu}}\right),\quad\Phi^{S}_{\mathcal{P}}=\mathrm{Im}\left(\frac{f_{k}-f_{\nu}}{1-f_{\nu}}\right).$
(4.29)

Figure 4.7: Normalized thermoacoustic production
$\Phi^{T}_{\mathcal{P}}\cos(\psi_{Up})-\Phi^{S}_{\mathcal{P}}\sin(\psi_{Up})$
(cf. Eq. (4.28a)) versus the phase angle difference between $\hat{p}$ and
$\hat{U}$, $\psi_{Up}=\angle{\hat{U}}-\angle{\hat{p}}$, for increasing values
of $h_{b}/2\delta_{k}=~{}$0.5, 1.0, 1.5, 2.0. (–); (- -), Optimum phasing
maximizing the thermoacoustic production for continually varying
$h_{b}/2\delta_{k}$.

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
4.6$a$ and 4.7). 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 4.7), for fixed temperature settings,
dissipation also diverges (figure 4.6$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 4.7). For the
temperature settings considered here (table 4.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.

####  IV.3 Effects of temperature gradient and geometry on growth rates

$(a)$$(b)$

Figure 4.8: Dimensionless thermoacoustic growth rate, $\alpha\tau$ (cf. Eq.
(4.32)), versus the natural logarithm of the temperature ratio $T_{H}/T_{C}$
for $T_{C}=300~{}$ K ($a$) and resonator area ratio $h_{d}/h_{c}$ ($b$).
($a$): (–), Linear stability analysis; (- -), logarithmic estimate of
$\alpha\tau$ Eq. (4.34) fitted using values at $T_{H}=450\,\mathrm{K}$ and
$T_{H}=550\,\mathrm{K}$; ($\circ$), Navier-Stokes simulations.

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

$\overline{\mathcal{P}}_{*}=\frac{\overline{\mathcal{P}}\tau}{\rho_{0}a^{2}_{0}},\quad\overline{\mathcal{D}}_{*}=\frac{\overline{\mathcal{D}}\tau}{\rho_{0}a^{2}_{0}},\quad\mathcal{R}_{*}=\frac{\mathcal{R}\tau}{h_{b}\rho_{0}a^{2}_{0}}.$
(4.30)

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

$\displaystyle\mathcal{R}_{T*}=\frac{\tau}{2(1-Pr)\rho_{0}a^{2}_{0}h^{2}_{b}}\int_{b}\Theta\,\mathrm{Re}\left(\hat{p}^{*}\hat{U}\right)\Phi^{T}_{\mathcal{P}}dx,$
(4.31a)
$\displaystyle\mathcal{R}_{S*}=-\frac{\tau}{2(1-Pr)\rho_{0}a^{2}_{0}h^{2}_{b}}\int_{b}\Theta\,\mathrm{Im}\left(\hat{p}^{*}\hat{U}\right)\Phi^{S}_{\mathcal{P}}dx,$
(4.31b)

respectively, are also analysed, where

$\tau=h^{2}_{b}/\nu_{0}$ (4.32)

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

$\overline{\mathcal{P}}_{*}\sim\Theta=\frac{d}{dx}\ln[T_{0}(x)],$ (4.33)

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

$\mathcal{R}_{*}\sim\alpha\tau\approx A\ln(T_{H}/T_{C})-B,$ (4.34)

where $A$ and $B$ are geometry dependent fitting coefficients.

$(a)$$(b)$

Figure 4.9: Dimensionless Rayleigh index $\mathcal{R}_{*}$ versus the natural
logarithm of the temperature ratio $T_{H}/T_{C}$ for $T_{C}=300\,$K ($a$) and
the resonator area ratio $h_{d}/h_{c}$ ($b$). ($b$): (–), $\mathcal{R}_{*}$ ;
(- -), $\mathcal{R}_{T*}$ ; ($\cdot\cdot\cdot$), $\mathcal{R}_{S*}$.

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 4.8$b$). Also,
the net production of acoustic energy density due to the travelling wave
component $\mathcal{R}_{T*}$ (figure 4.9$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. Eq. (4.29)).

###  V 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 BiwaEtAl_JASA_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 Eq. (4.49). In Section V.1 mass, momentum, and energy
(combined with the second law of thermodynamics) equations correct up to
second order are introduced. Furthermore, in Section V.2 and Section V.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 Section V.3.

####  V.1 Governing equations for nonlinear thermoviscous perturbations

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

$\displaystyle\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}=\left[-\rho^{\prime}\frac{\partial
u^{\prime}}{\partial x}-u^{\prime}\frac{\partial\rho^{\prime}}{\partial
x}\right],$ (4.35a) $\displaystyle\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}}-\frac{1}{\rho_{0}}\frac{\partial}{\partial
x}\left[\mu_{0}\left(\xi_{B}+\frac{4}{3}\right)\frac{\partial
u^{\prime}}{\partial
x}\right]=\left[-\frac{\rho^{\prime}}{\rho_{0}}\frac{\partial
u^{\prime}}{\partial t}-\frac{1}{2}\frac{\partial u^{\prime 2}}{\partial
x}\right],$ (4.35b) $\displaystyle\frac{\partial s^{\prime}}{\partial
t}+u^{\prime}\frac{ds_{0}}{dx}-\frac{Rk_{0}}{p_{0}}\frac{\partial^{2}T^{\prime}}{\partial
y^{2}}-\frac{R}{p_{0}}\frac{\partial}{\partial x}\left(k_{0}\frac{\partial
T^{\prime}}{\partial
x}\right)=\left[-\frac{p^{\prime}}{p_{0}}\left(\frac{\partial
s^{\prime}}{\partial
t}+u^{\prime}\frac{ds_{0}}{dx}\right)-u^{\prime}\frac{\partial
s^{\prime}}{\partial x}\right],$ (4.35c)

where, primed variables $\left(\cdot\right)^{\prime}$ represent the
fluctuations in corresponding quantities, the subscript $0$ 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$. Higher harmonics with correspondingly shorter
wavelengths are generated due to the nonlinear spectral energy cascade.
Consequently, the axial diffusion terms in Eq. (4.35b) and Eq. (4.35c), which
have been neglected in the linear regime, become significant and act as the
primary sink of energy at large harmonic scales in the spectral space.

We seek to collapse Eq. (4.35a)–(4.35c) to obtain a set of equations similar
to Eq. (4.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:

$\rho^{\prime}=\alpha_{s}p^{\prime}+\alpha_{p}s^{\prime}+\frac{1}{2}\left(\beta_{s}p^{\prime
2}+\beta_{p}s^{\prime 2}+2\beta_{sp}s^{\prime}p^{\prime}\right),$ (4.36)

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

$\displaystyle\alpha_{s}=\left(\frac{\partial\rho}{\partial
p}\right)_{s}=\frac{1}{a^{2}_{0}},\quad\alpha_{p}=\left(\frac{\partial\rho}{\partial
s}\right)_{p}=-\frac{\rho_{0}}{C_{p}},$ (4.37a)
$\displaystyle\beta_{s}=\left(\frac{\partial^{2}\rho}{\partial
p^{2}}\right)_{s}=-\frac{\gamma-1}{\rho_{0}a^{4}_{0}},\quad\beta_{p}=\left(\frac{\partial^{2}\rho}{\partial
s^{2}}\right)_{p}=\frac{\rho_{0}}{C^{2}_{p}},$ (4.37b)
$\displaystyle\beta_{sp}=\left[\frac{\partial}{\partial
s}\left(\frac{\partial\rho}{\partial
p}\right)_{s}\right]_{p}=\left[\frac{\partial}{\partial
s}\left(\frac{\rho}{\gamma p}\right)\right]_{p}=-\frac{1}{C_{p}a^{2}_{0}}.$
(4.37c)

Figure 4.10: Thermodynamic cycle in $p^{\prime}-\rho^{\prime}$ plane for non-
isentropic thermoacoustic wave amplification. (–), Nonlinear; (- -), Linear.
Sample perturbation fields are: $p^{\prime}=0.5p_{0}\cos\omega t$ and
$s^{\prime}=-Rp^{\prime}/p_{0}+250\sin\left(\omega(t+\tau)\right)$ for
$\tau=h^{2}_{b}/\nu=0.0107~{}\mathrm{s}$.

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 Eq. (4.36) account for the
corresponding nonlinear effects. HedbergR_JAP_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 4.10 shows the hysteresis cycle for the second order
constitutive relation in Eq. (4.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_NLA_1998 . Consequently, the
constitutive relation Eq. (4.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 Eq. (4.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}$ can
be derived and the final results are presented in the following sections.

####  V.2 Free-shear nonlinear wave propagation

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

$\displaystyle\frac{\partial p^{\prime}}{\partial t}=-\frac{\gamma
p_{0}}{h}\left(1+\frac{1+\gamma}{\gamma}\frac{p^{\prime}}{p_{0}}\right)\frac{\partial
U^{\prime}}{\partial
x}+\frac{k}{\rho_{0}}\left(\frac{1}{C_{v}}-\frac{1}{C_{p}}\right)\frac{\partial^{2}p^{\prime}}{\partial
x^{2}},$ (4.38) $\displaystyle\frac{\partial U^{\prime}}{\partial
t}=-\frac{h}{\rho_{0}}\frac{\partial p^{\prime}}{\partial
x}+\nu_{0}\left(\xi_{B}+\frac{4}{3}\right)\frac{\partial^{2}U^{\prime}}{\partial
x^{2}},$ (4.39)

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

####  V.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 ChuK_JFM_1958 ; Pierce_1989 , the following decomposition for
the entropy and velocity fluctuations inside the regenerator is considered:

$u^{\prime}(x,y,t)=\tilde{u}(x,t)+u^{\prime}_{\nu}(x,y,t),\quad
s^{\prime}=\tilde{s}(x,t)+s^{\prime}_{q}(x,y,t),$ (4.40)

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

$\frac{\partial u^{\prime}_{\nu}}{\partial
t}=\nu_{0}\frac{\partial^{2}u^{\prime}_{\nu}}{\partial y^{2}},\quad
u^{\prime}_{\nu}(x,y=\pm h_{b}/2,t)=-\tilde{u}(x,t).$ (4.41)

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

$\frac{\partial s^{\prime}_{q}}{\partial
t}+u^{\prime}\frac{ds_{0}}{dx}=\frac{\nu_{0}}{Pr}\frac{\partial^{2}s^{\prime}_{q}}{\partial
y^{2}},\quad s^{\prime}_{q}\left(x,y=\pm
h_{b}/2,t\right)=-\tilde{s}+s^{\prime}_{\mathrm{w}},$ (4.42)

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

$s^{\prime}_{\mathrm{w}}=-\frac{1}{\rho_{0}T_{0}}p^{\prime}=-\frac{R}{p_{0}}p^{\prime}.$
(4.43)

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

$\displaystyle
u^{\prime}_{\nu}=-\tilde{u}+\sum_{j=0}^{\infty}\check{u}_{j}(x,t)\cos(\zeta_{j}y),\quad
s^{\prime}_{q}=-\tilde{s}+s_{\mathrm{w}}+\sum_{j=0}^{\infty}\check{s}_{j}(x,t)\cos\left(\zeta_{j}y\right),$
(4.44a)
$\displaystyle\quad\text{with}\quad\zeta_{j}=\left(2j+1\right)\frac{\pi}{h_{b}}.$
(4.44b)

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

$\displaystyle\frac{\partial\check{u}_{j}}{\partial
t}+\nu_{0}\zeta^{2}_{j}\check{u}_{j}=(-1)^{j+1}\frac{2}{\zeta_{j}h_{b}}\left(\frac{1}{\rho_{0}}\frac{\partial
p^{\prime}}{\partial x}\right),$ (4.45)
$\displaystyle\frac{\partial\check{s}_{j}}{\partial
t}+\check{u}_{j}\frac{ds_{0}}{dx}+\frac{\nu_{0}}{Pr}\zeta^{2}_{j}\check{s}_{j}=(-1)^{j}\frac{2R}{\zeta_{j}h_{b}p_{0}}\frac{\partial
p^{\prime}}{\partial t}.$ (4.46)

Equations (4.45) and (4.46) determine the evolution of the transverse modes of
the longitudinal velocity $\check{u}_{j}$ and entropy $\check{s}_{j}$
fluctuations.

####  V.4 Nonlinear thermoviscous wave equations

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

$\displaystyle\frac{\partial U^{\prime}}{\partial
t}+\frac{h_{b}}{\rho_{0}}\frac{\partial p^{\prime}}{\partial
x}=\tau^{\prime}_{w}+\frac{1}{\rho_{0}}\frac{\partial}{\partial
x}\left[\mu_{0}\left(\xi_{B}+\frac{4}{3}\right)\frac{\partial
U^{\prime}}{\partial x}\right],$ (4.47)

where

$\displaystyle\tau^{\prime}_{w}=2\nu_{0}\sum^{\infty}_{j=0}(-1)^{j+1}\check{u}_{j}(x,t)\zeta_{j}.$
(4.48)

Equation (4.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, Eqs. (4.35a), (4.35c), and
(4.36) are combined to obtain

$\displaystyle\underbrace{\frac{\partial p^{\prime}}{\partial
t}+\frac{\rho_{0}a^{2}_{0}}{h_{b}}\frac{\partial U^{\prime}}{\partial
x}}_{\text{wave
propagation}}=\frac{\rho_{0}a^{2}_{0}}{h_{b}}\Bigg{(}\frac{q^{\prime}}{C_{p}\rho_{0}T_{0}}+\underbrace{q_{2}+\mathbb{T}-\mathbb{Q}}_{\begin{subarray}{c}\text{thermodynamic}\\\
\text{nonlinearities}\end{subarray}}+\mathbb{D}_{s}\Bigg{)}-\mathbb{C},$
(4.49)

where

$q^{\prime}=\frac{2\nu_{0}\rho_{0}T_{0}}{Pr}\sum^{\infty}_{j=0}(-1)^{j+1}\check{s}_{j}(x,t)\zeta_{j},$
(4.50)

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

$\mathbb{Q}=\frac{\gamma p^{\prime}q^{\prime}}{C_{p}p_{0}\rho_{0}T_{0}},$
(4.51)

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 Eq. (4.35b)
and Eq. (4.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 (4.49) and expressions for the terms
$q_{2}$, $\mathbb{D}_{s}$, $\mathbb{T}$, and $\mathbb{C}$ are given in
GuptaLS_JFM_2017 . Equations (4.45), (4.46), (4.47), and (4.49) constitute the
governing equations for the spatio-temporal evolution of large amplitude
acoustic perturbations inside the regenerator. The macrosonic thermoacoustic
interaction Eq. (4.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:

$\frac{\partial p^{\prime}}{\partial
t}+\frac{\rho_{0}a^{2}_{0}}{h_{b}}\frac{\partial U^{\prime}}{\partial
x}\approx\frac{\rho_{0}a^{2}_{0}}{h_{b}}\left\\{\frac{1}{C_{p}}\left[\left(1-\frac{\gamma
p^{\prime}}{p_{0}}\right)\frac{q^{\prime}}{\rho_{0}T_{0}}\right]\right\\}.$
(4.52)

In the above equation, terms $\mathbb{T}$, $\mathbb{D}_{s}$, and $q_{2}$ are
neglected for simplicity. Equations (4.45) and (4.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
Eqs. (4.45) and (4.46) yields:

$\displaystyle\check{u}_{j}=(-1)^{j+1}\frac{2}{\zeta_{j}h_{b}\rho_{0}}\int^{t}_{-\infty}e^{-\frac{t-\eta}{\tau_{j}}}\frac{\partial
p^{\prime}}{\partial x}(x,\eta)d\eta,$ (4.53)
$\displaystyle\check{s}_{j}=-\frac{ds_{0}}{dx}\int^{t}_{-\infty}e^{-\frac{t-\eta}{Pr\tau_{j}}}\check{u}_{j}(x,\eta)d\eta+(-1)^{j}\frac{2R}{\zeta_{j}hp_{0}}\int^{t}_{-\infty}e^{-\frac{t-\eta}{Pr\tau_{j}}}\frac{\partial
p^{\prime}}{\partial\eta}(x,\eta)d\eta.$ (4.54)

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,

$\mathcal{G}_{j}(\phi,\tau_{j})=\int^{t}_{-\infty}e^{-\frac{t-\eta}{\tau_{j}}}\phi(x,\eta)d\eta,$
(4.55)

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

$\displaystyle\tau^{\prime}_{w}=\frac{4\nu_{0}}{\rho_{0}h_{b}}\sum^{\infty}_{j=0}\mathcal{G}_{j}\left(\frac{\partial
p^{\prime}}{\partial x},\tau_{j}\right),$ (4.56) $\displaystyle
q^{\prime}=\frac{2\rho_{0}\nu_{0}T_{0}}{Pr}\sum^{\infty}_{j=0}\left[(-1)^{j}\frac{ds_{0}}{dx}\mathcal{G}_{j}\left(\zeta_{j}\check{u}_{j},\tau_{j}Pr\right)-\frac{2R}{h_{b}p_{0}}\mathcal{G}_{j}\left(\frac{\partial
p^{\prime}}{\partial t},\tau_{j}Pr\right)\right].$ (4.57)

Equations (4.56) and (4.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 (4.47) and (4.52), complete the nonlinear wave
propagation model equations. However, in the present work,Eqs. (4.45) and
(4.46) have been considered for time integration for simplicity.

The solution technique for model Eqs. (4.38), (4.39), (4.45), (4.46), (4.47),
and (4.49) is discussed in detail in GuptaLS_JFM_2017 .

###  VI Nonlinear Spectral Energy Dynamics

$(a)$$(b)$$(c)$

Figure 4.11: Comparison of filtered thermoacoustic signal as obtained from
time integration of nonlinear model Eqs. (4.38), (4.39), (4.45), (4.46),
(4.47), and (4.49) probed at $x=1.54~{}\mathrm{m}$ in harmonic growth regime
($a$), hierarchical spectral broadening regime ($b$), and limit cycle ($c$).
(–), Model with macrosonic thermoacoustic interaction $\mathbb{Q}$; (- -),
Model without $\mathbb{Q}$, ($\circ$), Navier-Stokes simulations.

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 (4.52) are in fairly good
agreement with the fully compressible Navier-Stokes simulations (see figure
4.11). Since several nonlinearities are neglected in calculations via (4.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 (4.51), the accuracy of the
predicted limit cycle amplitude gets reduced to $60\%$.

Figure 4.12 shows the time evolution of the spectral energy density $E_{k}$ of
the unstable mode (cf. Eq. (4.13)) and its first seven overtones, as obtained
from the nonlinear model discussed in Section V, and compares it to the
results obtained from the fully compressible Navier-Stokes simulations for the
signal shown in figure 4.11. 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$, 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 4.12
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 modeling. Assuming propagation of purely travelling waves in the
system and eliminating $U^{\prime}$ from Eqs. (4.38) and (4.39), the following
Burgers equation for the pressure fluctuations $p^{\prime}$ is obtained:

$\frac{\partial p^{\prime}}{\partial
t}-\frac{(\gamma+1)}{4\rho_{0}a_{0}}\frac{\partial p^{\prime
2}}{\partial\xi}=\frac{\delta}{2}\frac{\partial^{2}p^{\prime}}{\partial\xi^{2}},$
(4.58)

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

$\delta=\nu_{0}\left(\frac{4}{3}+\xi_{B}\right)+\frac{k}{\rho_{0}}\left(\frac{1}{C_{v}}-\frac{1}{C_{p}}\right).$
(4.59)

In general, the wall-shear and the wall-heat flux expressions 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 Eq. (4.58) yields:

$\frac{\partial p^{\prime}}{\partial
t}-\frac{(\gamma+1)}{4\rho_{0}a_{0}}\frac{\partial p^{\prime
2}}{\partial\xi}=\frac{\delta}{2}\frac{\partial^{2}p^{\prime}}{\partial\xi^{2}}+\alpha_{\mathrm{th}}p^{\prime},$
(4.60)

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

$p^{\prime}=\sum_{k}p_{k}(t)\sin\left(\frac{k\omega_{0}}{2a_{0}}\xi\right),~{}~{}\mathrm{where}~{}~{}k=2,4,6,8,\cdots$
(4.61)

the following modal evolution equation is obtained:

$\frac{dp_{k}}{dt}=\alpha_{\mathrm{th}}p_{k}+Q(p_{k})-\frac{\delta}{16}\left(\frac{\omega_{0}}{a_{0}}\right)^{2}k^{2}p_{k},$
(4.62)

where

$Q(p_{k})=\frac{(\gamma+1)\omega_{0}}{8\rho_{0}a^{2}_{0}}\Bigg{(}\sum^{n\leq
k-2}(k-n)p_{n}p_{k-n}-k\sum_{n\geq k+2}p_{n}p_{n-k}\Bigg{)},$ (4.63)

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

Figure 4.12: Comparison of evolution of the spectral energy density $E_{k}$
for the unstable mode and its first seven overtones as obtained from the
signal shown in figure 4.11. (–), Nonlinear model; ($\circ$), Navier-Stokes
simulations; (- -), Linear interpolation of spectral energy density evolution
in hierarchical spectral broadening regime.

$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 Eq. (4.62) 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 Eq. (4.62) can be approximated
such that

$\frac{dp_{2}}{dt}\approx\alpha
p_{2},\quad\frac{dp_{4}}{dt}\approx\frac{(\gamma+1)}{4\rho_{0}a^{2}_{0}}\omega_{0}p^{2}_{2},\quad\frac{dp_{6}}{dt}\approx\frac{3(\gamma+1)}{4\rho_{0}a^{2}_{0}}\omega_{0}p_{2}p_{4},$
(4.64)

and so on for higher harmonics. Equation (4.64) 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.

###  VII Scales of thermoacoustically sustained 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 modeling,
analogous to the turbulent energy cascade pope2000turbulent ;
Nazarenko_2011_WT . 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

$E^{(1\mathrm{D})}=\int^{\frac{2\pi
a_{0}}{\omega_{0}}}_{0}E\,d\xi=\sum_{k}E^{(1\mathrm{D})}_{k},\quad
E^{(1\mathrm{D})}_{k}=\frac{\pi a_{0}}{\omega_{0}}E_{k},$ (4.65)

where $E$ is the squared $L_{2}$ norm defined in Eq. (4.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 Eq. (4.59). Thus, $\mathcal{E}$
is purely governed by thermoacoustic wave amplification at smaller harmonics
and viscous dissipation at higher harmonics and scales as

$\frac{\mathcal{E}}{\alpha^{2}_{\mathrm{eff}}\delta}=\mathrm{const.},$ (4.66)

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

Figure 4.13: 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 ($\circ$),
$T_{H}=450~{}$K; ($+$), $T_{H}=500$ K; ($\square$), $T_{H}=550$ K; (- -)
compares the variation of the energy with harmonic frequency with the power
law derived in Eq. (4.69).

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

$\alpha_{\mathrm{eff}}\approx(1-\gamma|p^{\prime}|/p_{0})\alpha,$ (4.67)

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:

$\frac{E^{(1\mathrm{D})}_{k}\omega^{5/2}_{k}}{\rho_{0}\mathcal{E}^{3/2}}=\mathrm{const.}.$
(4.68)

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

$\frac{E_{k}\omega^{5/2}_{k}}{\mathcal{E}^{3/2}}=\mathrm{const.}.$ (4.69)

The scaling derived in Eq. (4.69) 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 4.13. 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 4.5) 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
GuptaLS_SCITECH_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
Eq. (4.69) can be further improved accounting for propagation of very large
harmonics inside the regenerator.

###  VIII Summary

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._ Yazaki_PhysRevLet_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 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.

##  Chapter 5 Spark-plasma generated shock waves

The contents of this chapter were presented in the AIAA-Scitech presentation
in January 2019 BsinghEtAl_2019_AIAA and have been reported here in abridged
form. The work presented here was motivated by the experimental study of wave
propagation due to spark-plasma heat deposition by Singh  _et al._
Bsingh2018_EinF

###  I Introduction

Plasma discharges are a growing field of interest with numerous applications
across a wide range of areas in science and engineering including flow and
combustion control, biomedicine, materials processing, nanotechnology, and
environmental engineering Corke_ARFM_2010 . In the area of aerodynamic flow
control, considerable work over the past couple of decades has focused on
using plasma actuators as active flow control devices on aerodynamic bodies.
Nanosecond pulsed plasmas have recently been attracting great interest due to
their extremely efficient generation of excited and ionized species and
relatively low power consumption. Recent studies have investigated the effect
of various parameters of the nanosecond pulse driven plasma discharges on the
nature of the plasma Pai_PSST_2010 . The aim of this study is to analyze the
flow field generated as a result of high-frequency repeated plasma between two
electrodes with the help of high-fidelity simulations, low order modeling, and
experiments.

###  II Experimental study and motivation (Credit: Ms. Bhavini Singh, Mr.
Lalit Rajendran, and Dr. Sally Bane)

Figure 5.1 shows the snapshots of repeated spark discharge between two
electrodes as seen in PIV measurements. Multiple bursts of pulses were sent to
the pulse-generator which creates the repeated spark discharge. The
instantaneous heat deposition due to the spark discharge generates a heat-
induced cylindrical shock wave which propagates outward and is deflected by
the electrode geometry.

Figure 5.1: PIV measurements of repeated spark discharges at 1 kHz of
frequency between the electrodes 5 mm apart (Courtesy of Ms. Bhavini Singh,
Mr. Lalit Rajendran, and Dr. Sally Bane). Figure 5.2: Shock wave induced due
to the spark plasma heat deposition.

Behind the shock wave, vortical structures are generated due to grazing of the
shock wave over hard no-slip walls of electrodes. Such vortical structures are
utilized in the boundary layer re-attachment over aerodynamic objects to
improve lift characteristics. Figure 5.2 shows the preliminary simulation
results obtained via modeling the heat deposited by a single spark event by a
simplistic isothermal boundary condition.

###  III Computational setup

A high-order unstructured, fully compressible Navier-Stokes solver (sd3DvisP)
is used to simulate the shocks induced by the heat deposition by the spark.
The sd3DvisP solver is an MPI parallelized fortran 90 code for compressible
flows based on the high-order Spectral Difference (SD) scheme for unstructured
hexahedral elements. kopriva:96 ; sun:07 The solver is capable of running with
arbitrary preselected orders of accuracy and provides minimal numerical
dissipation. lodato:14b For the computation of high-speed compressible flows
with shocks and discontinuities, a self-calibrating shock capturing
methodology is available in the solver. Originally developed for the
Discontinuous Galerkin scheme, persson:06 this method, enables sub-cell
resolution of discontinuities which are detected through a modal sensor. The
sd3DvisP solver has been successfully applied to the computation of flows with
shocks, including shock/vortex interaction and shock/wavy-wall interactions.
lodato:15 ; lodato:16

Figure 5.3: Computational setup of the electrodes. The setup models
axisymmetric flow conditions.

In the scope of the current work, we investigate the flow field generated in
an axisymmetric computational setup (figure 5.3). The governing equations
solved are conservation laws for mass, momentum, and total energy in three
dimensions, given by,

$\displaystyle\frac{\partial}{\partial t}\left(\rho\right)$
$\displaystyle+\frac{\partial}{\partial x_{j}}\left(\rho
u_{j}\right)=\sigma\left(\rho_{\infty}-\rho\right),$ (5.1a)
$\displaystyle\frac{\partial}{\partial t}\left(\rho u_{i}\right)$
$\displaystyle+\frac{\partial}{\partial x_{j}}\left(\rho
u_{i}u_{j}\right)=-\frac{\partial}{\partial x_{i}}p+\frac{\partial}{\partial
x_{j}}\tau_{ij}+\sigma\left(\left(\rho u\right)_{\infty}-\left(\rho
u\right)\right),$ (5.1b) $\displaystyle\frac{\partial}{\partial
t}\left(\rho\,E\right)$ $\displaystyle+\frac{\partial}{\partial
x_{j}}\left[u_{j}\left(\rho\,E+p\right)\right]=\frac{\partial}{\partial
x_{j}}\left(u_{i}\tau_{ij}-q_{j}\right)+\dot{\mathcal{Q}}_{\mathrm{spark}}+\sigma\left(\left(\rho
E\right)_{\infty}-\rho E\right),$ (5.1c)

respectively, are solved where $x_{i}$($x_{1}$, $x_{2}$, $x_{3}$ or
equivalently, $x$, $y$, $z$) are the cartesian 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. We model the energy deposition due to the spark
plasma utilizing the $\dot{\mathcal{Q}}_{\mathrm{spark}}$ term in the energy
equation given by,

$\dot{\mathcal{Q}}_{\mathrm{spark}}=\mathcal{Q}_{0}\exp\left(-\alpha\left({t-\tau}\right)^{2}\right),$
(5.2)

where,

$\mathcal{Q}_{0}=\rho_{0}C_{p}{\alpha\tau T_{s}}g(x,y,z),$ (5.3)

$T_{s}$ is the characteristic spark temperature, and $\alpha$ and $\tau$ are
the spark energy deposition rate parameters. Function $g(x,y,z)$ determines
the shape of the spark. For simplicity, we model the spark as a straight
channel with heat deposition such that,

$g(x,y,z)=\exp\left(-8\times 10^{8}\left(R-R_{cut}\right)^{2}\right),$ (5.4)

where $R_{cut}=0.6~{}\mathrm{mm}$. The exponential relaxation terms
(characterized by $\sigma$ in Eqs. 4.5a-4.5c) are used for non-reflective
sponge layers (borders distinguished by red line in figure 5.3) which model
the outflow boundaries. Within the domain $|x|<x_{s}$ and $R<R_{s}$,
$\sigma=0$ is used and outside the domain,

$\displaystyle\sigma=\frac{f(x,y,z)}{\Delta t}$
$\displaystyle=\frac{1}{4\Delta
t}\left(1+\tanh\left(\beta|x-x_{s}|\right)\right)\left(1+\tanh\left(\beta\left(R(y,z)-R_{s}\right)\right)\right),$
(5.5) $\displaystyle\mathrm{for}~{}~{}|x|>x_{s},~{}~{}R>R_{s},$

is used. Parameter $\beta$ determines the spatial variation of the sponge
layers. In the current simulations $x_{s}=12.7~{}\mathrm{mm}$ and
$R_{s}=10~{}\mathrm{mm}$ are chosen.

Below we discuss two simulation results, one corresponding to no-slip
adiabatic walls of the electrodes (viscous), and the other corresponding to
hard-slip walls (inviscid). These simulations enable comparison of viscous and
inviscid vorticity generation mechanisms. The spark parameters used in the
current simulations are,

$T_{s}=2000~{}\mathrm{K},~{}\tau=20~{}\mathrm{ns}~{}\mathrm{and}~{}\alpha=1\times
10^{17}~{}\mathrm{s}^{-2}.$ (5.6)

###  IV Results

Figure 5.4 illustrates the basic flow features observed behind the curved
shock induced by the heat deposition. As the temperature of the discharge
region drops after reaching a maximum, an expansion wave is generated behind
the shock wave which creates a backward flow (flow towards the electrode
center) in the domain. Moreover, due to the shock curvature, the vorticity
generated also enhances the flow towards the domain. Below, we quantitatively
discuss the effect of baroclinic torque on the vorticity.

Figure 5.4: Schematic illustrating the flow features generated by the heat
deposition induced shock wave and the curvature in the shock wave caused by
the electrode geometry.

Figure 5.5 shows the initial stages of vorticity $\omega_{n}$ normal to $x-y$
plane. As the heat deposition peaks, a cylindrical shock wave emanates from
the spark. Due to the slant boundaries, the shock wave diffracts and develops
curvature. At the initial stages, the shock wave diffraction from the walls
causes high localized baroclinic torque which results in production of
vorticity. Moreover, behind the curved shocks, boundary layer develops in case
of no-slip electrode walls (figure 5.5). The vorticity field near the
electrode walls is modified considerably due to no-slip walls (figure 5.5)
compared to slip walls (figure 5.6).

Behind the shock wave, due to sudden decrease in heat rate (after reaching
maximum), an expansion wave is generated. As a result, the flow decelerates
and is directed inward (towards the center of the spark) by the electrode
walls. Since, the current simulations model spark boundary as a symmetric
boundary, the flow is re-directed upwards, hence driving the large scale
vortices generated. While the magnitude of vorticity is higher, the sense of
vorticity is similar to the one observed in experiments. Due to the
simulations being axisymmetric, the vortices do not breakdown. Thus, no
smaller scale vortices are observed and the large scale vortices simply decay
due to fluid viscosity.

Figure 5.5: Evolution of the vorticity $\omega_{n}$ normal to the $x-y$ plane
in the initial stages of the shock wave propagation for no-slip electrode
walls.

Figure 5.6: Evolution of the vorticity $\omega_{n}$ normal to the $x-y$ plane
in the initial stages of the shock wave propagation for slip electrode walls.

Figure 5.6 shows the the vorticity evolution for a simulation with identical
$\alpha$, $\tau$, and $T_{s}$ but only hard-wall symmetric condition imposed
on the electrode walls. Due to fluid slip, the walls do not act as a source of
vorticity. However, the large scale vortices can still be observed with
similar magnitude of $\omega_{n}$ thus suggesting that inviscid baroclinic
terms result in initial generation of vorticity.

Figure 5.7: Evolution of the normalized baroclinic torque $\Gamma_{n}$ normal
to the $x-y$ plane in the initial stages of the shock wave propagation for no-
slip electrode walls.

Figures 5.7 and 5.8 show the evolution of baroclinic torque $\Gamma_{n}$
normal to the x-y plane, given by,

$\Gamma_{n}=\frac{\left(\nabla\rho\times\nabla
p\right)_{n}}{\rho^{2}}\left(\frac{\rho_{0}^{2}l^{2}}{\Delta
p\Delta\rho}\right),$ (5.7)

where $l$ is the characteristic length scale of pressure and density gradients
($l\sim 1~{}\mathrm{mm}$), $\Delta p$ and $\Delta\rho$ are the pressure and
density jumps across the shock wave. In the current simulations (both no-slip
wall and slip wall), the pressure jump varies from $\Delta p\sim
110~{}\mathrm{kPa}$ in the center to $\Delta p\sim 50~{}\mathrm{kPa}$ near the
electrode walls. Similarly, the density jump varies from $\Delta\rho\sim
0.8~{}\mathrm{kg/m}^{3}$ in the center to $\Delta\rho\sim
0.5\sim\mathrm{kg/m}^{3}$ near the electrode walls. Values near the electrode
walls have been used for normalization in Eq. 5.7.

Figure 5.8: Evolution of the normalized baroclinic torque $\Gamma_{n}$ normal
to the $x-y$ plane in the initial stages of the shock wave propagation for
slip electrode walls.

In figures 5.7 and 5.8, we note that the baroclinic torque $\Gamma_{n}$ in
both the slip wall and no-slip wall cases evolves similarly. Moreover, the
sense of the torque matches the vorticity alignment shown in figures 5.5 and
5.6 thus indicating that baroclinic torque indeed plays an important role in
initial generation of the large scale vortices observed in both the
simulations and the experiments.

##  Chapter 6 Numerical modeling of detonation waves

The contents of this chapter were presented in the AIAA-Aviation meeting in
June 2018 GuptaSLSS_Aviation_2018 and have been reported here in abridged
form. The work presented here was motivated by the experimental study by
Schwinn  _et al._ Schwinn_CandF_2018 at Maurice J. Zucrow Laboratories at
Purdue.

###  I Introduction

Combustion waves are primarily classified as deflagration or detonation waves,
the latter being high amplitude pressure waves propagating at supersonic
speeds. Heating caused by a propagating adiabatic shock-compression results in
ignition of fuel-oxidizer mixture. The fluid containing reaction products
expands behind the shock wave and depending on downstream boundary conditions,
further accelerates the wave front sustaining its propagation JHlee_2008 .
Classical Zel’dovich, von Neumann and Döring (ZND) model postulates
equilibrium one-dimensional detonation waves. In experiments and detailed
theoretical studies, unstable detonation wave propagation with complex
reactive chemistry and compressible flow physics interactions are often
observed JHlee_2008 ; SHEPHERD2009 . Consequently, unsteady dynamics of
unstable one-dimensional pulsating detonation waves have received wide-spread
attention Erpenbeck_PoF_1964 ; Lee_JFM_1990 ; Sharpe_PRSA_1999 .

Detonation waves result in very high pressure gains, which exhibit higher
thermodynamic efficiencies when utilized for mechanical work Wolanski_RDE_2015
; WintenbergerShepherd_JPP_2006 . Rotating Detonation Engines (RDEs) are
propulsion devices, which utilize such pressure gains through continuously
spinning detonation waves for generating thrust Bykovskii_JPP_2006 . Fuel-
oxidizer mixture is injected axially into an annular shaped device which is
undergoing combustion due to the rotating detonation waves. High pressure
combustion expels products axially from the opposite end thus generating
thrust. However, the complex combustion wave propagation dynamics in such
devices require further attention and careful analysis for effective design
Schwinn_CandF_2018 . In this work, we perform numerical investigations of such
sustained detonation dynamics in a periodic domain, inspired by the recent
experimental study of Schwinn _et al._ Schwinn_CandF_2018 in which a straight-
line detonation chamber was analyzed experimentally exhibiting sustained
resonance of detonation waves.

Numerical studies of RDEs usually focus on the general flow features exhibited
by sustained detonation waves Schwer_AIAA_2010 ; Schwer_CandF_2011 ;
Lietz_AIAA_2018 . In this work, we numerically investigate the effect of fuel-
oxidizer injection rates on the dynamics of detonation waves and sustenance in
a periodic one-dimensional domain. We model fuel-oxidizer injection rates
utilizing a one-dimensional model reduction which represents solving the
governing equations close to the injector plate in an experimental setup.
Moreover, we adopt a device-scale dynamical system perspective in the current
study to elucidate the wave propagation dynamics.

###  II Experimental study and motivation (Credit: Mr. Kyle Schwinn and Dr.
Carson Slabaugh)

Schwinn _et al._ Schwinn_CandF_2018 conducted an experimental study on
sustained detonation wave resonance in a linear semi-bounded channel with
CH4+O2 mixture. Figure 6.1 shows a schematic of the experimental setup of
Schwinn _et al._ Schwinn_CandF_2018 with the instrumentation points and the
two possible window locations for obtaining optical measurements. The
experimental setup provided modular capability between making single-point
measurements with pressure transducers and ion probes, and obtaining high
frame rate images of the chamber dynamics through large, fused quartz windows.
Pressure fluctuations in the chamber were measured using piezoelectric
pressure transducers (PCB 113B26).

Figure 6.1: Schematic diagram of instrumentation port nomenclature for the
combustion chamber (CC), the fuel manifold (FM), and the oxidizer manifold
(OM). White panels in the combustion chamber indicate regions of optical
access.(Courtesy of Schwinn _et al._ Schwinn_CandF_2018 )

Using different diagnostics on opposites sides of the narrow chamber,
simultaneous pressure and chemiluminescence measurements of the chamber
combustion dynamics were performed. Figure 6.2$(a)$ shows a detonation wave
(6.2A) on the left side of the image, propagating in the positive x-direction
and approaching the instrumentation port location CC-09 (see figure 6.1). The
average intensity at probe location CC-09, which is proportional to the local
combustion heat release, is shown in figure 6.2$(b)$ along with the in-phase
local pressure fluctuation.

$(a)$$(b)$

Figure 6.2: a) Detonation imaged using OH*-chemiluminescence (Window B), with
instrumentation port CC-09 marked with red circle. b) Overlay of intensity
data with pressure data from instrument at location CC-09. (Courtesy of
Schwinn _et al._ Schwinn_CandF_2018 )

Schwinn _et al._ Schwinn_CandF_2018 observed wave steepening and amplification
caused by the fuel-oxidizer injection through a time series of
OH*-chemiluminescence, as shown in figure 6.3. In the first frame (with
arbitrary time assignment $+0\,\mu s$), a combustion front (figure 6.3A)
propagates from an auto-ignition source at the interface between the reactants
and product gases. In this same frame, a high-intensity compact combustion
wave (6.3B) from the previous cycle is traveling in the positive x-direction.
The recently-ignited combustion front propagates upstream into the reactant
jets and reaches the injector face, and then travels along the injector face
in the positive x-direction (6.3C). Between two successive waves, deflagrative
burning in the reactant fill is observed as a region of low intensity (6.3D).
As the combustion front continues to propagate in the positive x-direction
into the reactants, a distinct contact surface forms between the reactants and
product gases (6.3E). The combustion wave accelerates due to the heat release
behind the shock, intensifying and becoming more spatially compact (6.3F). As
the steep-fronted combustion wave continues to propagate through the channel,
the expanding combustion front from an auto-ignition event appears (6.3G),
beginning the process for the next cycle.

Figure 6.3: Time series of OH*-chemiluminescence images (imaged in Window A)
detailing the combustion zones in the development region of the combustion
chamber (Scale is $10\,mm$; see first frame).(Courtesy of Schwinn _et al._
Schwinn_CandF_2018 )

By inspection of images in figure 6.3, the combustion waves are determined to
be supersonic and the nature of the limit-cycle behavior of the chamber
dynamics was linked to periodic injection, ignition, amplification, and
acceleration processes. The high-speed chemiluminescence images obtained in
Window B, as seen in figure 6.2 enable the identification of structures in the
reaction zone that confirm the presence of self-excited, periodic detonation
waves.

Motivated by the experimental study conducted by Schwinn _et al._
Schwinn_CandF_2018 , we study the possibility of a reduced order one-
dimensional detonation wave resonance model in this work to elucidate the
nature of sustained detonation wave resonance from the perspective of
dynamical systems. We note that the fuel-oxidizer injection and coupling of
injector nozzles is crucial to allow the spatial variation of fuel injection
coupled with the local pressure variations caused by the detonation waves. In
our study, we assume premixed fuel-oxidizer injec human being tion of H2+O2+Ar
(due to ease of achieving detonations) mixture through sonic (converging)
nozzles in the chamber. Detailed model for hydrogen-oxygen-argon reaction is a
9-species and 19-step reaction scheme based on Jachimowski’s mechanisms for
hydrogen-air combustion Jachimowski1988_Hydrogen . In the following sub-
section, we discuss the validation of detailed chemistry model in a high-order
spectral difference solver, sd3DvisP lodato2013discrete ; lodato2014structural
; lodato:14b ; GuptaLS_JFM_2017 used in the previous chapter.

###  III Numerical validation

Detonation systems are usually modeled with single step reactive chemistry
models based on induction times in which progress variable of reaction
replaces the species concentrations in dynamical equations. The dynamics of
resonating detonations are elucidated in this work based on detailed chemistry
models utilizing the state-of-the-art reaction mechanisms and spectral
difference solver utilized in previous chapters.

####  III.1 Validation with Shock-to-detonation transition

$(a)$$(b)$$(c)$

Figure 6.4: Validation test setup $(a)$ identical to the one studied by Oran
_et al._ Oran_CandF_1982 and Yungster and Radhakrishnan
YungsterRadha_AIAA_1996 . The combustion front due to shock-heating at the
adiabatic (symmetric) wall accelerates $(b)$ eventually merging with the
reflected shock wave. The resulting steady state detonation wave propagates
steadily at $v_{D}\approx 1600~{}\mathrm{m/s}$ $(c)$. The pressure peak in the
detonation wave corresponds to the Von-Neumann state
$p_{vn}=560156.5~{}\mathrm{Pa}$.

The numerical model is validated against one-dimensional steady detonation
wave propagation, as shown in figure 6.4. The validation setup is the one
defined by Yungster and Radhakrishnan YungsterRadha_AIAA_1996 . In the
validation study, the simulation is initialized with a shock wave of mach
number $M_{s}=3$, propagating towards an adiabatic (symmetric) wall and fluid
mixture at $p_{1}=6666.12~{}\mathrm{Pa},T_{1}=300K$. On the right hand side,
an inlet is defined driving the shock wave in a mixture of H2+O2+$Ar$ with
molar ratio $2:1:7$. The reflection of the shock wave results in instantaneous
shock-heating at the adiabatic wall, which ignites the mixture. The combustion
wave generated due to the instantaneous heating initially lags behind the
reflected shock wave (figure. 6.4$(a)$). The combustion front accelerates
(figure. 6.4$(b)$) eventually merging with the reflected shock wave and
steadily propagates with the velocity $v_{D}$ thus resulting in a steady
propagating detonation wave. The steady state velocity $v_{D}$ saturates to
$1600~{}\mathrm{m}/\mathrm{s}$ as approximated by the steady state ZND theory.
The merging of combustion wave front and the shock wave increases the peak
pressure (figure 6.4$(c)$) as well. The peak pressure corresponds to the Von-
Neumann state in the steady state detonation propagation theory JHlee_2008 and
matches the value $p_{vn}=560156.5~{}\mathrm{Pa}$. The temperature of the Von-
Neumann state is $T_{vn}=2459.33~{}\mathrm{K}$ and the maximum heat release
rate behind the Von-Neumann state $\omega_{T}=2.58\times
10^{11}~{}\mathrm{W}/\mathrm{m}^{3}$.

The shock-to-detonation transition validation study (figure 6.4) validates the
numerical model for detailed chemistry implemented. We note that the shock-
capturing viscosity is approximately $100$ times higher than the maximum
physical diffusivity in the system and hence only supports stable numerical
computations. The inviscid and chemical reaction based characteristics of the
detonation waves are captured accurately. In the section below, we outline the
one-dimensional problem formulation based on previous numerical studies
Schwer_CandF_2011 ; Lietz_AIAA_2018 and experimental studies Fotia_AIAA_2014 ;
Schwinn_CandF_2018 for sustained detonation waves.

####  III.2 Pre-mixed laminar flame

The pre-mixed laminar CH4-air flame was simulated at various equivalence
ratios $\phi$ for validating the combined implementation of transport and
combustion terms in the governing equations. Since the formulation in the
current work is fully compressible, the results differ from the usual low-Mach
number implementation results of the energy conservation equation. The
computation is initialized with a fully converged steady state solution
obtained from CANTERA goodwin2009cantera . Due to fully compressible
formulation and fixed Lewis numbers for all the species, the spatial
distribution of velocity and temperature adjusts to the steady state solution
shown in figure 6.5 with a corresponding change in pressure. Relaxation
sponges are used at the extremes of the domain to eliminate thermoacoustic
interactions and convergence to steady state solution. The
$\mathrm{CH}_{4}-\mathrm{air}$ mixture (left) flows towards a stationary
flame. Flow speed before the flame corresponds to the flame velocity $u_{f}$
in an unsteady case.

$(a)$$(b)$

Figure 6.5: Spatial distribution of fluid velocity and temperature (top) and
mass fractions of $\mathrm{CH}_{4}$ and $\mathrm{CH}^{.}_{3}$ (bottom) for a
laminar flame testcase for equivalence ratio $\phi=1$ $(a)$. Comparison of
CH4-air laminar flame speeds for various equivalence ratios with experimental
data vagelopoulos1994 , low-mach number unity lewis number calculation ($-$)
PoinsotV_numComb_2011 , and current calculations ($\circ$) $(b)$

Smaller intermediate species like H$\cdot$ and H2 have higher diffusivity
coefficients which would tend to increase the flame speed
($u_{f}\propto\sqrt{K\mathcal{D}}$ where $K$ is some measure of rate of
reaction and $\mathcal{D}$ is the species diffusivity). Hence, the simplified
assumption of $Le_{k}=1$ for all species results in under-prediction of flame
speeds compared to the experiments.

###  IV One-dimensional problem formulation

Detonation propagation can be sustained via fuel injection and exhaust purge
as shown in various experiments Fotia_AIAA_2014 ; Schwinn_CandF_2018 . Various
numerical studies of such sustained detonations in two dimensions have also
been conducted Schwer_CandF_2011 ; Lietz_AIAA_2018 . However, two-dimensional
studies focus more on the flow characteristics induced by sustained
detonations and are computationally expensive. In this study, particular focus
is given to establishing which parameters affect the sustenance of detonation
waves. To the best of our knowledge, no previous attempt has been made to
derive first-principle based 1D model capable of simulating sustained
detonation waves.

$(a)$$(b)$

Figure 6.6: Schematic illustrating the one-dimensional model reduction based
on the OH∗ chemiluminescence visualization of experimental study by Schwinn
_et al._ Schwinn_CandF_2018 $(a)$ and sketch of the spatial profile of
pressure corresponding to the visualized detonation wave $(b)$. Dashed line (-
-) shows the stagnation injection pressure.

We reduce the problem to a one-dimensional problem as shown in figure 6.6$(a)$
by considering a small region of height $h$ in which gradients in $y$
direction can be neglected. Two dimensional governing equations can be written
as,

$\displaystyle\frac{\partial\rho}{\partial t}+\frac{\partial\left(\rho
u\right)}{\partial x}+\frac{\partial\left(\rho v\right)}{\partial y}=0,$ (6.1)
$\displaystyle\frac{\partial\rho u}{\partial t}+\frac{\partial\rho
u^{2}}{\partial x}+\frac{\partial\rho vu}{\partial y}=-\frac{\partial
p}{\partial x}+\frac{\partial}{\partial x}\left(\frac{4\mu}{3}\frac{\partial
u}{\partial x}\right)+\frac{\partial}{\partial y}\left(\mu\frac{\partial
u}{\partial y}+\mu\frac{\partial v}{\partial x}\right),$ (6.2)
$\displaystyle\frac{\partial\rho E}{\partial t}+\frac{\partial}{\partial
x}\left(u(\rho E+p)\right)+\frac{\partial}{\partial y}\left(v(\rho
E+p)\right)=\frac{\partial}{\partial x}\left(k\frac{\partial T}{\partial
x}\right)+\frac{\partial}{\partial y}\left(k\frac{\partial T}{\partial
y}\right)$ $\displaystyle+\frac{\partial}{\partial
x}\left(u\frac{4\mu}{3}\frac{\partial u}{\partial x}+v\mu\left(\frac{\partial
u}{\partial y}+\frac{\partial v}{\partial
x}\right)\right)+\frac{\partial}{\partial y}\left(u\mu\left(\frac{\partial
u}{\partial y}+\frac{\partial v}{\partial
x}\right)+v\frac{4\mu}{3}\frac{\partial v}{\partial y}\right)+$
$\displaystyle\omega_{T}-\frac{\partial}{\partial
x}\left(\rho\sum^{N}_{k=1}h_{s,k}Y_{k}V_{k,x}\right)-\frac{\partial}{\partial
y}\left(\rho\sum^{N}_{k=1}h_{s,k}Y_{k}V_{k,y}\right),$ (6.3)
$\displaystyle\frac{\partial\rho Y_{k}}{\partial t}+\frac{\partial}{\partial
x}\left(\rho(u+V_{k,x})Y_{k}\right)+\frac{\partial}{\partial
y}\left(\rho(v+V_{k,y}+V^{C}_{k,y})Y_{k}\right)=\omega_{k}.$ (6.4)

Integrating the Eqs. (6.1)-6.4 in $y$ within the region of height $h$ shown in
figure 6.6$(a)$ and neglecting stresses, diffusion velocities, and temperature
gradients at inlet and outlet, we obtain the following system of equations in
1D,

$\displaystyle\frac{\partial\rho}{\partial t}+\frac{\partial\left(\rho
u\right)}{\partial x}=\frac{1}{h}\left(\left(\rho
v\right)_{\mathrm{inlet}}-\left(\rho v\right)_{\mathrm{outlet}}\right),$ (6.5)
$\displaystyle\frac{\partial(\rho u)}{\partial t}+\frac{\partial(\rho
u^{2})}{\partial x}=-\frac{\partial p}{\partial x}+\frac{\partial}{\partial
x}\left(\frac{4\mu}{3}\frac{\partial u}{\partial
x}\right)+\frac{1}{h}\left(\left(\rho uv\right)_{\mathrm{inlet}}-\left(\rho
uv\right)_{\mathrm{outlet}}\right),$ (6.6) $\displaystyle\frac{\partial(\rho
E)}{\partial t}+\frac{\partial}{\partial x}\left(u(\rho
E+p)\right)=\frac{\partial}{\partial x}\left(k\frac{\partial T}{\partial
x}\right)+\frac{\partial}{\partial x}\left(u\frac{4\mu}{3}\frac{\partial
u}{\partial x}+v\mu\left(\frac{\partial u}{\partial y}+\frac{\partial
v}{\partial x}\right)\right)$
$\displaystyle+\omega_{T}-\frac{\partial}{\partial
x}\left(\rho\sum^{N}_{k=1}h_{s,k}Y_{k}V_{k,x}\right)+\frac{1}{h}\left(\left(v(\rho
E+p)\right)_{\mathrm{inlet}}-\left(v(\rho
E+p)\right)_{\mathrm{outlet}}\right),$ (6.7) $\displaystyle\frac{\partial(\rho
Y_{k})}{\partial t}+\frac{\partial}{\partial
x}\left(\rho(u+V_{k,x})Y_{k}\right)=\omega_{k}+\frac{1}{h}\left(\left(\rho
vY_{k}\right)_{\mathrm{inlet}}-\left(\rho
vY_{k}\right)_{\mathrm{outlet}}\right).$ (6.8)

Equations (6.5)-6.8 represent the one-dimensional model equations accounting
for fuel/oxidizer injection and outlet close to the injector plate in the full
experimental setup. We assume $h=1~{}\mathrm{mm}$ in the current study based
on experimental visualizations. We note that values of $h$ should be much
smaller than the overall transverse dimension of the rig to account for
negligible $y$ gradients in the domain. Furthermore, we assume that
$v_{\mathrm{outlet}}=v_{\mathrm{inlet}}=v_{\mathrm{inj}}$ where
$v_{\mathrm{inj}}$ is evaluated assuming isentropic expansion through injector
nozzles, as introduced by Schwer and Kailasnath Schwer_AIAA_2010 . The
assumption of $v_{\mathrm{outlet}}=v_{\mathrm{inlet}}$ will be corrected in
future studies to account for the flux of reaction products out of the domain
correctly. The outlet quantities $\rho$, $p$, $\rho u$, $\rho E$, and $\rho
Y_{k}$ correspond to the quantities in the domain. Consequently, the model
equations become,

$\displaystyle\frac{\partial\rho}{\partial t}+\frac{\partial\left(\rho
u\right)}{\partial
x}=\frac{v_{\mathrm{inj}}}{h}\left(\rho_{\mathrm{inlet}}-\rho\right),$ (6.9)
$\displaystyle\frac{\partial(\rho u)}{\partial t}+\frac{\partial(\rho
u^{2})}{\partial x}=-\frac{\partial p}{\partial x}+\frac{\partial}{\partial
x}\left(\frac{4\mu}{3}\frac{\partial u}{\partial
x}\right)+\frac{v_{\mathrm{inj}}}{h}\left(-\left(\rho u\right)\right),$ (6.10)
$\displaystyle\frac{\partial(\rho E)}{\partial t}+\frac{\partial}{\partial
x}\left(u(\rho E+p)\right)=\frac{\partial}{\partial x}\left(k\frac{\partial
T}{\partial x}\right)+\frac{\partial}{\partial
x}\left(u\frac{4\mu}{3}\frac{\partial u}{\partial x}+v\mu\left(\frac{\partial
u}{\partial y}+\frac{\partial v}{\partial x}\right)\right)$
$\displaystyle+\omega_{T}-\frac{\partial}{\partial
x}\left(\rho\sum^{N}_{k=1}h_{s,k}Y_{k}V_{k,x}\right)+\frac{v_{\mathrm{inj}}}{h}\left(\left(\rho
E+p\right)_{\mathrm{inlet}}-(\rho E+p)\right),$ (6.11)
$\displaystyle\frac{\partial(\rho Y_{k})}{\partial t}+\frac{\partial}{\partial
x}\left(\rho(u+V_{k,x})Y_{k}\right)=\omega_{k}+\frac{v_{\mathrm{inj}}}{h}\left(\left(\rho
Y_{k}\right)_{\mathrm{inlet}}-\rho Y_{k}\right).$ (6.12)

To calculate $v_{\mathrm{inj}}$ and $\rho_{\mathrm{inlet}}$ in the model Eqs.
(6.9)-6.12, we specify the stagnation conditions at the injectors
$p_{0,\mathrm{inj}},T_{0,\mathrm{inj}}$. Through isentropic nozzle expansion
relations, the mass flow rate of the injected fluid is given by,

$\rho_{\mathrm{inlet}}v_{\mathrm{inj}}=p_{0,\mathrm{inj}}\sqrt{\frac{\gamma}{RT_{0,\mathrm{inj}}}}M\left(1+\frac{\gamma-1}{2}M^{2}\right)^{-\frac{1}{2}\left(\frac{\gamma+1}{\gamma-1}\right)},$
(6.13)

where $M$ is the exit Mach number of the premixed fuel/oxidiser mixture
governed by the pressure in the domain as,

$M=\sqrt{\frac{2}{\gamma-1}\left(\left(\frac{p_{0,\mathrm{inj}}}{p}\right)^{\frac{\gamma-1}{\gamma}}-1\right)},$
(6.14)

and $\gamma$ is the adiabatic exponent of the mixture. Where
$p>p_{0,\mathrm{inj}}$, the mixture injection does not take place. We note
that, isentropic expansion relations also yield the density injected
$\rho_{\mathrm{inj}}$ inside the domain as,

$\rho_{\mathrm{inj}}=\rho_{0,\mathrm{inj}}\left(1+\frac{\gamma-1}{2}M^{2}\right)^{\frac{1}{1-\gamma}}.$
(6.15)

For $p>p_{0,\mathrm{inj}}$, there is no fuel injection. For
$p_{s,\mathrm{inj}}<p<p_{0,\mathrm{inj}}$ the fuel injection rates are
determined by Eqs. (6.13)-6.15, where $p_{s,\mathrm{inj}}$ is the sonic
pressure corresponding to $M=1$ in equation (6.14). For $p<p_{s,\mathrm{inj}}$
in the domain, all the quantities are calculated at the sonic state, $M=1$.

Figure 6.7: Initial propagation of detonation waves in a periodic domain at
$t_{1}=0$, $t_{2}=10~{}\mu\mathrm{s}$, and $t_{3}=20~{}\mu\mathrm{s}$ for
$p_{0,\mathrm{inj}}=0.3~{}\mathrm{MPa}$ and $w_{i}=0.2$.

###  V Results

Table 6.1: Injection stagnation pressures $p_{0,\mathrm{inj}}$ and the weighting factors $w_{i}$ considered. $p_{0,\mathrm{inj}}$ | 0.2 MPa | 0.3 MPa | 0.4 MPa
---|---|---|---
| 0 | 0 | 0
| 0.1 | 0.1 | 0.1
| 0.2 | 0.2 | 0.2
$w_{i}$ | - | 0.4 | -
| - | 0.6 | -
| - | 0.8 | -
| - | 1.0 | -

In the current study, we restrict the analysis to stoichiometric H2+O2+Ar
mixture with molar fractions $2:1:3$ and consider three injection pressures
$p_{0,\mathrm{inj}}$ to elucidate the dynamics of detonation wave propagation
with premixed fuel injection. We note that increasing stagnation pressure
$p_{0,\mathrm{inj}}$ causes higher density of the mixture injected in the
domain which results in higher heat release rates due to combustion. Since
high heat release rates cause very large gradients in species concentrations,
the numerical simulations become unstable. In this study, to achieve stable
computations, we assume the inlet density $\rho_{\mathrm{inlet}}$ as an
average between the injected density $\rho_{\mathrm{inj}}$ and the initialized
density $\rho_{\infty}$ as,

$\rho_{\mathrm{inlet}}=w_{i}\rho_{\mathrm{inj}}+(1-w_{i})\rho_{\infty}.$
(6.24)

$(a)$$(b)$$(c)$$(d)$

Figure 6.8: Spatio-temporal evolution of pressure $p$ (–) and heat release
$\omega_{T}$ (–) (left) ; mass fraction of Hydrogen $Y_{H_{2}}$ (–) and
temperature $T$ (–) (right) for times $t=0$ $(a)$, $12~{}\mu\mathrm{s}$ $(b)$,
$30~{}\mu\mathrm{s}$ $(c)$, and $58~{}\mu\mathrm{s}$ $(d)$ for
$p_{0,\mathrm{inj}}=0.3~{}\mathrm{MPa}$ and $w_{i}=0.2$. The initialized pulse
propagates to the right and steepens thus forming a shock wave of $M\approx
6$. The rapid temperature rise behind the shock wave causes the ignition of
the fuel thus resulting in detonation wave formation.

We note that such choice of averaged inlet density is akin to a numerical
under-relaxation and is un-physical. However, it yields preliminary physical
insights into the detonation waves reproducible by one-dimensional model
reduction. We consider a periodic domain of length $L=5~{}\mathrm{cm}$ (see
figure 6.6) and initialize the system with a stoichiometric mixture of
H2+O2+Ar with molar fractions $2:1:3$. The pressure $p$, density $\rho$, and
velocity $u$ in the domain are initialized with an isentropic Gaussian pulse
(see figure 6.7) propagating in $+x$ direction. The initial field is given as,

$\displaystyle p(x,t=0)$ $\displaystyle=p_{\infty}+\Delta pf_{x},$ (6.25)
$\displaystyle\rho(x,t=0)$
$\displaystyle=\rho_{\infty}\left(\frac{p}{p_{\infty}}\right)^{1/\gamma_{\infty}},$
(6.26) $\displaystyle u(x,t=0)$
$\displaystyle=\frac{2a_{\infty}}{\gamma_{\infty}-1}\left(\left(\frac{\rho}{\rho_{\infty}}\right)^{\frac{\gamma_{\infty}-1}{2}}-1\right),$
(6.27)

where $p_{\infty}=0.02~{}\mathrm{MPa}$, $\Delta p$ is the initial pressure
jump and $f_{x}$ is the Gaussian given by,

$\Delta
p=0.82~{}\mathrm{MPa},~{}f_{x}=\exp\left(-10^{5}\left(x-0.02\right)^{2}\right),$
(6.28)

respectively. The pulse propagates and steepens thus forming a shock wave of
$M\approx 6$. The temperature jump thus caused ignites the mixture of H2+O2+Ar
such that a combustion front is created behind the shock wave. Unlike the
validation case shown in figure 6.4, the detonation wave is not driven by an
external flow as a result of which, the pressure falls behind the wave instead
of saturating to a constant value. As the pressure falls below
$p_{0,\mathrm{inj}}$ the fuel injection is activated thus reducing the
temperature and increasing the fuel/oxidizer fractions in the mixture. Table
6.1 lists the different values of $p_{0,\mathrm{inj}}$ and $w_{i}$ considered
in the current study. Throughout the computations, we assume constant
injection stagnation temperature $T_{0,\mathrm{inj}}=700~{}\mathrm{K}$.

Figure 6.8 below shows the spatio-temporal evolution of pressure $p$, heat
release rate $\omega_{T}$, mass fraction of Hydrogen $Y_{H_{2}}$, and the
temperature $T$ in the domain for $p_{0,\mathrm{inj}}=0.3~{}\mathrm{MPa}$ and
$w_{i}=0.2$. The shock wave formation due to steepening of the initialized
pulse causes rapid heating of the mixture behind the shock wave up to
temperature $T\approx 1000~{}\mathrm{K}$. As a result, the mixture in the
domain ignites resulting in detonation waves. Due to no driving support (for
instance, a moving piston), the pressure downstream of the detonation wave
decreases generating an expansion wave. As the pressure decreases below the
injection stagnation pressure $p_{0,\mathrm{inj}}$, the fuel injection is
activated locally. As the detonation wave continues to propagate, the fuel
injection regime follows the wave. However, the expansion wave behind the
detonation front flattens in space thus causing the fuel injection zone to lag
behind the detonation front.

$(a)$$(b)$$(c)$$(d)$$(e)$$(f)$$w_{i}=0$$w_{i}=0.1$$w_{i}=0.2$

Figure 6.9: Evolution of the maximum heat release $\omega_{T}$ behind the
detonation wave in the system for $w_{i}=0$ $(a)$, $w_{i}=0.1$ $(b)$, and
$w_{i}=0.2$ $(c)$ for injection stagnation pressures
$p_{0,\mathrm{inj}}=0.2~{}\mathrm{MPa}$ (–), 0.3 MPa ($--$), and 0.4 MPa
($\cdots$).

The flattening of the expansion wave downstream of the detonation wave results
in widening of the gap between the fuel injection zone and the detonation
wave. Consequently, detonation wave eventually merges with the fuel injection
zone resulting in quenching of the reaction and failure of the detonation wave
sustenance. Figure 6.9 shows the evolution of maximum heat release
$\omega_{T}$ behind the detonation wave as the wave propagates in the periodic
domain for $w_{i}=0,~{}0.1,~{}0.2$ and increasing values of injection
stagnation pressures
$p_{0,\mathrm{inj}}=0.2~{}\mathrm{MPa},~{}0.3~{}\mathrm{MPa},~{}0.4~{}\mathrm{MPa}$.
With increasing stagnation pressures, the density of the injected fuel
increases. Since the chemical composition of the mixture is fixed, the total
amount of fuel and oxidizer increases in the system thus causing higher heat
release behind the detonation wave.

As shown in figure 6.9, detonations exhibiting higher heat release sustain for
more time compared to those which exhibit lower heat release. However, for
$w_{i}=0$, the heat release is independent of the injection stagnation
pressure since the density of the fluid remains constant in the domain. For
such conditions, lower injection stagnation pressures increase the separation
between the detonation front and the fuel injection zone. Hence, the waves are
quenched by the fuel injection at later times for lower injection stagnation
pressures sustaining waves for longer. Figure 6.10 shows the trajectory of the
combustion front in time compared with the corresponding C-J velocities,
obtained by CEA CEAcite . The combustion front velocity is very close to the
C-J velocity initially. However, the detonation wave decelerates to velocity
slightly below C-J velocity slowly and eventually gets quenched.

$(a)$$(b)$$(c)$$w_{i}=0$$w_{i}=0.1$$w_{i}=0.2$

Figure 6.10: Evolution of the detonation wave trajectory $x_{D}$ in time for
$w_{i}=0$ $(a)$, $w_{i}=0.1$ $(b)$, and $w_{i}=0.2$ $(c)$ for injection
stagnation pressures $p_{0,\mathrm{inj}}=0.2~{}\mathrm{MPa}$ (–), 0.3 MPa
($--$), and 0.4 MPa ($\cdots$). (–) are the C-J velocities for
$p_{0,\mathrm{inj}}=0.4$ for various $w_{i}$. Figure 6.11: Evolution of the
maximum heat release $\omega_{T}$ behind the detonation wave in the system for
$w_{i}=0.4,~{}0.6,~{}0.8$ and $w_{i}=1.0$ for injection stagnation pressure
$p_{0,\mathrm{inj}}=$ 0.3 MPa ($--$).

Similar to increasing the injection stagnation pressure $p_{0,\mathrm{inj}}$,
with increasing the weighting factor $w_{i}$ also, the density of the injected
mixture increases thus causing high heat release rates, as shown in figure
6.11. With increasing release rates, the detonation velocity becomes higher.
This results in the merging of detonation front and the fuel injection zone.
Consequently, the detonation wave fails and the pressure becomes greater than
the injection stagnation pressure everywhere in the domain. Such over-pressure
in the domain highlights the unphysical nature of the assumption
$v_{\mathrm{outlet}}=v_{\mathrm{inlet}}$. In reality, high pressure generated
due to combustion of reactants expels the products at very high speeds
compared to the inlet speeds. Such increase in outlet momentum, in-fact,
results in generation of useful thrust.

###  VI Conclusions

A numerical investigation of planar detonation wave sustenance due to fuel
injection in a periodic domain, inspired by the experimental investigation of
Schwinn _et al._ Schwinn_CandF_2018 , was performed utilizing spectral
difference scheme. The detonation waves were modeled using detailed chemical
reaction scheme and the results were analyzed for H2+O2+Ar mixture in a
periodic domain of $5~{}\mathrm{cm}$ length. The fuel injection was modeled
with one-dimensional model reduction representative of a region close to the
injector plates in an original experimental setup. The assumptions in the
fuel-injection model result in failure of detonation wave sustenance.

The initiation of detonation waves was done by a high pressure pulse which
steepens and forms shock waves. Rapid heating due to the shock wave causes the
fuel/oxidizer mixture in the domain to ignite. The combustion front thus
formed follows the shock wave resulting in a detonation wave. Downstream of
the detonation wave, an expansion wave is generated. The pressure profile
becomes flattened as the expansion and the detonation waves propagate. The
fuel injection is activated as the pressure falls below the injection
stagnation pressure. The fuel injection zone travels with the expansion wave,
lagging behind the detonation wave. As a result, the detonation wave merges
with the fuel-injection zone, thus causing quenching of the wave. We reported
the wave propagation and failure for various densities and stagnation pressure
of injected fuel.

We note that, the model assumptions are highly restrictive in representing the
fuel-injection and exhaust phenomenon observed in experiments
Schwinn_CandF_2018 , or multi-dimensional simulations Lietz_AIAA_2018 ;
Schwer_AIAA_2010 thus causing failure of detonation wave sustenance. Moreover,
the length of the domain, restricted by the computational cost in the current
study, also plays an important role in wave sustenance due to the time taken
by the detonation wave to merge with the fuel-injection zone.

##  Chapter 7 Development of an order adaptive Navier Stokes solver

The contents of this chapter summarize the development of a mesh adaptive
compressible flow solver, undertaken as a part of summer internship in 2018 at
HySonic Technologies, LLC. The numerical method discussed here-in is central
to the numerical results presented previously in this dissertation and all the
computations showed were performed utilizing the presented numerical method.
However, the particular mesh adaptive solver discussed here is developed as a
future tool to conduct more in-depth and efficient analysis of the problems
presented in previous chapters, in particular, resonating detonation waves.

###  I Introduction

High order discontinuous finite element methods such as spectral element
methods karniadakis2005spectral , spectral volume methods Wang_2002_JCP , and
multidomain spectral difference methods kopriva1996conservative provide
flexibility of obtain highly accurate solutions to partial differential
equations in geometrically complex domains GuptaLS_JFM_2017 ; LodatoJ_2019_LES
; Zhou_2010_JCP . In particular, spectral difference methods are characterized
by very less numerical diffusion for increasing orders of accuracy lodato:14b
; ChapelierLJ_CF_2016 thus providing highly accurate solutions to nonlinear
wave propagation problems GuptaScalo_PRE_2018 ; GuptaLS_JFM_2017 ; lodato:16
and simulation of compressible turbulent flows lodato2013discrete ;
LodatoJ_2019_LES . In such methods, each subdomain is solved utilizing high
order polynomial interpolations local to the subdomain. Consequently, the
methodology exhibits compactness of the discretization stencil and can be
optimally parallelized.

In this work, we develop a methodology to adapt the local orders of the
subdomains to obtain a high degree of polynomial reconstruction near the flow
features exhibiting high gradients utilizing the spectral difference methods.

###  II Spectral difference methods : Overview

In this section, an overview of spectral difference methods for fully
compressible Navier Stokes equations for an ideal gas is given. We begin with
presenting the governing equations, followed by discussing the discretized
solution evaluation methodology specific to unstructured hexahedral meshes
with straight-edged elements. However, we note that, such methods can be
easily extended to generic shaped meshes containing tetrahedra and prisms
Taylor_2000_SIAM ; Luo_2006_Luo ; karniadakis2005spectral . In further
sections, we discuss the adaptive order refinement strategy for multidomain
spectral difference methods.

####  II.1 Governing equations

Fully compressible three dimensional Navier-Stokes equations can be written in
the following conservative form,

$\displaystyle\frac{\partial\mathbf{Q}}{\partial
t}+\frac{\partial\mathbf{F}_{i}}{\partial
x_{i}}=\frac{\partial\mathbf{F^{v}}_{i}}{\partial x_{i}},$ (7.1)

where $t$ is time, $x_{i}$ represent spatial directions for $i=1,2,3$,
$\mathbf{Q}$ is the vector of conservative variables, namely, mass, momentum
in $1,2,3$ directions, and total energy, with $\mathbf{F}_{i}$ and
$\mathbf{F^{v}}_{i}$ being the inviscid and viscous flux vectors in $i^{th}$
direction respectively, given by,

$\mathbf{Q}=\left(\begin{matrix}\rho\\\ \rho u_{1}\\\ \rho u_{2}\\\ \rho
u_{3}\\\ \rho e\\\
\end{matrix}\right),\quad\mathbf{F}_{i}=\left(\begin{matrix}\rho u_{i}\\\ \rho
u_{1}u_{i}+p\delta_{i1}\\\ \rho u_{2}u_{i}+p\delta_{i2}\\\ \rho
u_{3}u_{i}+p\delta_{i3}\\\ \left(\rho e+p\right)u_{i}\\\
\end{matrix}\right),\quad\mathrm{and}\quad\mathbf{F^{v}}_{i}=\left(\begin{matrix}0\\\
\tau_{1i}\\\ \tau_{2i}\\\ \tau_{3i}\\\ \tau_{ij}u_{j}-q_{i}\\\
\end{matrix}\right).$ (7.2)

For an ideal gas, the total specific energy $e$ is given by,

$e=\frac{p}{\rho\left(\gamma-1\right)}+\frac{u_{i}u_{i}}{2}.$ (7.3)

The viscous stress tensors $\tau_{ij}$, assuming a Newtonian fluid, are given
as,

$\tau_{ij}=2\mu\left(\frac{\partial u_{i}}{\partial x_{j}}+\frac{\partial
u_{j}}{\partial x_{i}}+\left(\beta-\frac{2\mu}{3}\right)\frac{\partial
u_{k}}{\partial x_{k}}\delta_{ij}\right),$ (7.4)

where $\mu$ and $\beta$ are the coefficients of dynamic (shear) viscosity and
bulk viscosity, respectively. The heat flux $q_{i}$ is given by,

$q_{i}=-k\frac{\partial T}{\partial x_{i}},$ (7.5)

where $k$ is the coefficient of conduction, and related to the dynamic shear
viscosity via the Prandtl number,

$Pr=\frac{\mu C_{p}}{k},$ (7.6)

where $C_{p}$ is the heat capacity at constant pressure.

####  II.2 Mapping and spatial discretization

In the spectral difference method, the domain is assumed to be divided into
non-overlapping subdomains or elements. The elements can be hexahedral,
prismatic, or even tetrahedral. However, in this work, we focus on development
for hexahedral elements only. Each hexahedral element is mapped on a standard
computational domain, with the mapping function given by,

$x_{i}=\sum^{N}_{n=1}\chi^{n}\left(\xi_{j}\right)\xi^{n}_{i},$ (7.7)

where $\chi^{n}$ is the mapping function corresponding to the $n^{th}$ vertex
used for defining the physical element, $x_{i}$ are the physical coordinates,
and $\xi_{i}$ are the computational coordinates. The standard computational
domain or element is assumed to extend from $\left[-1,1\right]$ such that,

$\left(\xi_{1},\xi_{2},\xi_{3}\right)\in[-1,1]\times[-1,1]\times[-1,1].$ (7.8)

We note that the total number of vertices $N$ depends on the physical shape of
the element. For quadratic or higher order elements, isoparametric mapping
functions of same order can be utilized for mapping. However, in this work, we
focus only on trilinear elements, i.e., linear in all directions, which
correspond to $N=8$ for hexahedra.

Figure 7.1: Mapping $\chi$ between a standard hexahedral domain (left) defined
in computational space ($\xi_{1},~{}\xi_{2},~{}\xi_{3}$) and a skewed
trilinear hexahedral defined in physical space ($x_{1},~{}x_{2},~{}x_{3}$).

For the mapping function $\chi$ in Eq. 7.7, the jacobian of transformation is
given by the determinant,

$J=\left|\frac{\partial x_{i}}{\partial\xi_{j}}\right|,$ (7.9)

which plays an important role in transformation of the conservation Eqs. 7.2,
which take the following form in the mapped computational coordinate system
Landmann_2010_AIAA ,

$\frac{\partial\tilde{\mathbf{Q}}}{\partial
t}+\frac{\partial\tilde{\mathbf{F}}_{i}}{\partial\xi_{i}}=\frac{\partial\tilde{\mathbf{F^{v}}}_{i}}{\partial\xi_{i}},$
(7.10)

where $\tilde{\cdot}$ denotes quantities in the computational coordinate
system and are given by,

$\tilde{\mathbf{Q}}=J\mathbf{Q},~{}\mathrm{and}~{}\tilde{\mathbf{F}}_{i}=J\frac{\partial\xi_{i}}{\partial
x_{j}}\mathbf{F}_{j},~{}\mathrm{and}~{}\tilde{\mathbf{F^{v}}}_{i}=J\frac{\partial\xi_{i}}{\partial
x_{j}}\mathbf{F^{v}}_{j}.$ (7.11)

$(a)$$(b)$

Figure 7.2: Solution and flux point distribution in 3D ($a$) and 2D projection
on the $\xi_{1}-\xi_{3}$ plane ($b$) for $P=5$. Solution points are shown by
black circles and red crosses denote flux points.

Inside the standard hexahedron, two sets of points are defined namely,
solution points (figure 7.2) and flux points (figure 7.2), which are staggered
w.r.t. each other. The conservative variables $\mathbf{Q}$ are stored on the
solution points inside the hexahedron and the fluxes are computed on the flux
points. The solution points and flux points are the tensor product of one-
dimensional Gauss-Legendre and Gauss-Lobatto-Legendre points, respectively
kopriva1996conservative . Here onwards, we denote the set of solution points
as,

${X_{i}},~{}\mathrm{for}~{}i=1,2,..P,~{}\mathrm{and}~{}{X_{i+1/2}},~{}\mathrm{for}~{}i=0,1,2,..P,$
(7.12)

respectively. With $P$ solution points and $P+1$ in each direction inside an
element, following Lagrange basis functions are used for interpolations and
differentiations,

$h_{i}(X)=\prod^{P}_{j=1,j\neq
i}\left(\frac{X-X_{j}}{X_{i}-X_{j}}\right),~{}\mathrm{and}~{}l_{i+1/2}(X)=\prod^{P}_{j=0,j\neq
i}\left(\frac{X-X_{j+1/2}}{X_{i+1/2}-X_{j+1/2}}\right).$ (7.13)

Polynomials in Eq. 7.13 yield the following reconstructions of the vectors
$\tilde{\mathbf{Q}}$ and $\tilde{\mathbf{F}}_{i},\tilde{\mathbf{F^{v}}}_{i}$
in the computational space,

$\displaystyle\tilde{\mathbf{Q}}$
$\displaystyle=\sum^{P}_{k=1}\sum^{P}_{j=1}\sum^{P}_{i=1}\tilde{\mathbf{Q}}_{ijk}h_{i}(\xi_{1})h_{j}(\xi_{2})h_{k}(\xi_{3}),$
(7.14) $\displaystyle\tilde{\mathbf{F}}_{1},\tilde{\mathbf{F^{v}}}_{1}$
$\displaystyle=\sum^{P}_{k=1}\sum^{P}_{j=1}\sum^{P}_{i=0}\tilde{\mathbf{F}}_{1,ijk},\tilde{\mathbf{F^{v}}}_{1,ijk}l_{i+1/2}(\xi_{1})h_{j}(\xi_{2})h_{k}(\xi_{3}),$
(7.15) $\displaystyle\tilde{\mathbf{F}}_{2},\tilde{\mathbf{F^{v}}}_{2}$
$\displaystyle=\sum^{P}_{k=1}\sum^{P}_{j=0}\sum^{P}_{i=1}\tilde{\mathbf{F}}_{2,ijk},\tilde{\mathbf{F^{v}}}_{2,ijk}h_{i}(\xi_{1})l_{j+1/2}(\xi_{2})h_{k}(\xi_{3}),$
(7.16) $\displaystyle\tilde{\mathbf{F}}_{3},\tilde{\mathbf{F^{v}}}_{3}$
$\displaystyle=\sum^{P}_{k=0}\sum^{P}_{j=1}\sum^{P}_{i=1}\tilde{\mathbf{F}}_{3,ijk},\tilde{\mathbf{F^{v}}}_{3,ijk}h_{i}(\xi_{1})h_{j}(\xi_{2})l_{k+1/2}(\xi_{3}),$
(7.17)

which are used to advance the conservative variable vector
$\tilde{\mathbf{Q}}$ in time at solution points. Normal flux at the element
faces is calculated utilizing approximate Riemann solver (HLLC)
Toro_2013_Riemann and averaged diffusion flux. For time integration, we
utilize the strong stability preserving, five step RK3 ruuth2006global scheme,
however, we note that any order time integration can be allowed by the spatial
discretization methodology described above.

###  III Adaptive order refinement

Local order of each computational subdomain can be adapted according to the
order of accuracy required. For steady state problems, Kopriva
Kopriva_1996_Mortar deviced a way utilizing mortar elements to calculate the
interface fluxes for non-conforming set of points for steady state problems.
In this work, we extend the mortar element method to the unsteady problems
which require high spatial resolution locally.

Due to local order refinement, interface of various computational subdomains
(called mortars hereafter) might be connected across non-conformal set of flux
points, as shown in figure 7.3. As a result, local interpolation along the
face of computational elements is required.

Figure 7.3: Mortar element connecting elements with different local orders
with non-conforming flux points.

To this end, we consider the reconstruction of left and right states along the
interface, $\tilde{\mathbf{Q}}^{L}$, $\tilde{\mathbf{Q}}^{R}$ respectively,

$\displaystyle\tilde{\mathbf{Q}}^{L}(X)=\sum^{P_{L}}_{i=1}\sum^{P_{L}}_{j=1}\tilde{\mathbf{Q}}^{L}_{ij}h^{L}_{i}(\xi_{1})h^{L}_{j}(\xi_{2}),\quad\tilde{\mathbf{Q}}^{R}(X)=\sum^{P_{R}}_{i=1}\sum^{P_{R}}_{j=1}\tilde{\mathbf{Q}}^{R}_{ij}h^{R}_{i}(\xi_{1})h^{R}_{j}(\xi_{2}),$
(7.18)

where $P_{L}$ and $P_{R}$ denote the orders of computational elements on the
left and right of the mortar and $\xi_{1},\xi_{2}$ are the local coordinates
on the mapped face. On the mortar, we consider total number of points to be
$P_{M}=\mathrm{max}\left(P_{L},P_{R}\right)$ in each direction. To
interpolate, we utilize the least squares approximation, as done by Kopriva
Kopriva_1996_Mortar . To this end, we assume left and right interpolated
states on the mortar to have the following reconstruction,

$\tilde{\mathbf{Q}}^{L_{M},R_{M}}=\sum^{P_{M}}_{i=1}\sum^{P_{M}}_{j=1}\tilde{\mathbf{Q}}^{L_{M},R_{M}}_{ij}h^{M}_{i}(\xi_{1})h^{M}_{j}(\xi_{2}).$
(7.19)

The least squares approximation yields the following relation for the
interpolation from element to mortar on the left side Kopriva_1996_Mortar ,

$\displaystyle\sum^{P_{M}}_{i=1}\sum^{P_{M}}_{j=1}\tilde{\mathbf{Q}}^{L_{M}}_{ij}\int_{\mathrm{mortar}}h^{M}_{i}(\xi_{1})h^{M}_{j}(\xi_{2})h^{M}_{p}(\xi_{1})h^{M}_{q}(\xi_{2})d\xi_{1}d\xi_{2}=$
$\displaystyle\sum^{P_{L}}_{i=1}\sum^{P_{L}}_{j=1}\tilde{\mathbf{Q}}^{L}_{ij}\int_{\mathrm{mortar}}h^{L}_{i}(\xi_{1})h^{L}_{j}(\xi_{2})h^{M}_{p}(\xi_{1})h^{M}_{q}(\xi_{2})d\xi_{1}d\xi_{2},$
(7.20)

and similar for the right side. Solving Eq. 7.20 for every mortar yields the
interpolated values $\tilde{\mathbf{Q}}^{L_{M}}_{ij}$ and
$\tilde{\mathbf{Q}}^{R_{M}}_{ij}$ from the left and right sides of the mortar
on the mortar points. Utilizing these, fluxes
$\tilde{\mathbf{F}}_{i,pq}\left(\tilde{\mathbf{Q}}^{L_{M}}_{ij},\tilde{\mathbf{Q}}^{R_{M}}_{ij}\right)$
on the mortar points are evaluated and interpolated back on the left and right
elements utilizing the Eq. 7.20.

To adapt the local element orders based on the flow features, such as vortices
and shock waves, we detect the local flow velocity gradients. For shock waves,
we also add a localized shock capturing bulk viscosity KawaiLele_2008_JCP ;
Premasuthan_2014_I ; Premasuthan_2014_II based on the local magnitude of
divergence of the velocity field, which indicates the degree of
compression/rarefaction in the flow.

To detect the shocks, we calculate the local normalized divergence of the
velocity field

$\epsilon_{s}=\frac{1}{2}\left(\tanh\left(c_{1}-c_{2}\frac{\Delta
x}{a}\left|\sum^{3}_{j=1}\frac{\partial
u_{i}}{\partial\xi_{j}}\frac{\partial\xi_{j}}{\partial
x_{i}}\right|\right)-\tanh\left(c_{1}\right)\right),$ (7.21)

where $c_{1}>0$ and $c_{2}>0$ are the tuning parameters, $a$ is the local
speed of sound, and $\Delta x$ is the characteristic element thickness. The
error is bounded by $0$ and $1$ $0<\epsilon_{s}<1$, and based on predefined
tolerances, the local order of the elements is increased. Furthermore, a shock
capturing bulk viscosity and corresponding increase in thermal conductivity
are also added in the stress tensor calculation (see Eq. 7.4 ) and heat flux
calculation (see Eq. 7.5),

$\displaystyle\beta^{*}=c_{\beta}\rho\epsilon_{s}\left(\Delta
x\right)^{2}\left|\sum^{3}_{j=1}\frac{\partial
u_{i}}{\partial\xi_{j}}\frac{\partial\xi_{j}}{\partial x_{i}}\right|,$ (7.22)
$\displaystyle k^{*}=C_{p}\beta^{*}/Pr.$ (7.23)

Figure 7.4: Spatio-temporal evolution of the one dimensional Sod shock tube
problem and the adaptive order refinement.

To detect the local vorticity, we utilize the ratio of vorticity to velocity
magnitude as the sensing parameter, and define the vorticity sensor as,

$\epsilon_{v}=\frac{1}{2}\left(\tanh\left(c_{1}-c_{2}\Delta
x\frac{\omega_{i}\omega_{i}}{u_{i}u_{i}}\right)-\tanh\left(c_{1}\right)\right).$
(7.24)

Local resolution error is defined as the
$\epsilon=\mathrm{max}\left(\epsilon_{s},\epsilon_{v}\right)$ and refinement
or coarsening is done comparing $\epsilon$ with predefined tolerances.

###  IV Test cases

In this section, we discuss the three test cases which test the shock and
vorticity detecting sensors for adaptive order refinement.

####  IV.1 1D Shock tube

The first test case that we consider is the 1D Sod shock tube problem
SOD_1978_JCP with the dimensionless left and right states defined as,
$\left(\rho_{L},u_{1,L},p_{L}\right)$ = $\left(1,0,1\right)$ and
$\left(\rho_{R},u_{1,R},p_{R}\right)=\left(0.125,0,0.1\right)$ in a 1 unit
length domain with 10$\times$2 elements. Order refinement is also donw
corresponding to the initial conditions so as to capture the initial
conditions as accurately as possible. Figure 7.4 shows the spatio-temporal
evolution of the dimensionless pressure and adaptive order refinement of the
mesh. The test case was initialized with $P=2$ throughout and a maximum order
of 7.

Figure 7.5: Comparison of the solution obtained from the adaptive
$p$-refinement strategy with the analytical solution of the Sod shock tube
case at dimensionless time $t=0.25$.

We note that the sensor tracks the right propagating shock and the left
propagating expansion wave very effectively. Figure 7.5 shows the comparison
of the solution at $t=0.25$ units with the analytical solution of the inviscid
problem.

####  IV.2 Convecting isentropic vortex

The second test case considered is a convecting isentropic vortex in a base
flow of 0.5 Mach number in domain with $19\times 10$ elements. The vortex is
defined through the following dimensionless velocity and temperature values,

$\displaystyle u_{1}$
$\displaystyle=U_{\infty}+\varepsilon\frac{r}{r_{0}}e^{\alpha\left(1-\left(r/r_{0}\right)^{2}\right)}\sin\theta,$
(7.25) $\displaystyle u_{2}$
$\displaystyle=-\varepsilon\frac{r}{r_{0}}e^{\alpha\left(1-\left(r/r_{0}\right)^{2}\right)}\cos\theta,$
(7.26) $\displaystyle T$
$\displaystyle=T_{\infty}-\frac{\left(\gamma-1\right)\varepsilon^{2}e^{2\alpha\left(1-\left(r/r_{0}\right)^{2}\right)}}{4\alpha\gamma},$
(7.27)

where $U_{\infty}=0.5$ and $T_{\infty}=1$. Initial pressure and density are
evaluated through the isentropic relation $p/\rho^{\gamma}=\mathrm{const.}$
and the ideal gas equation of state. The vortex related parameters considered
were $\varepsilon=0.3,r_{0}=0.05,\alpha=0.204$. Jiang_1996_WENO

The vortex core propagates at the base flow velocity. Moreover, acoustic waves
are emanated from the vortex initialization which propagate due to temperature
initialization and the corresponding pressure and density changes.

Figure 7.6: Spatio temporal evolution of the isocontours of $u_{1}$ of the
vortex defined by velocity and temperature field in Eq. 7.27 (top), adaptive
order refinement as the vortex core moves (mid), and comparison of the
vertical velocity $u_{2}$ at $x_{2}=0.5$ with the analytical solution.

The adaptive order refinement follows the vortex core as well as the acoustic
waves propagating away from the vortex, as shown in figure 7.6. The domain is
initialized with $P=1$ in all the cells and then pre-refined before the first
time integration step based on the initial conditions. The maximum allowed
order is $P=7$. In figure 7.6, the comparison of $u_{2}$ with the analytical
result of pure advection with a velocity of $U_{\infty}$ is also shown.

Figure 7.7: Spatio temporal evolution of pressure isocontours for the shock
vortex interaction test case and the corresponding dynamic mesh refinement
based on both the shock sensor (see Eq. 7.21) and the vorticity sensor (see
Eq. 7.24).

####  IV.3 Shock Vortex interaction

The third and final test case considered in this work is the shock-vortex
interaction case with a shock of Mach number 1.1. Figure 7.7 shows the
adaptive order refinement and iso-contours of pressure field as the vortex
passes through the shock in a 2D domain with $32\times 16$ elements. The
domain is initialized with $P=1$ in all the cells and then pre-refined before
the first time integration step based on the initial conditions. The maximum
allowed order is $P=7$. We note that when the vortex passes through the shock,
the shock bifurcates thus creating triple points in the domain. Moreover, due
to the shock-vortex interaction, intense pressure oscillations at the element
boundaries are generated due to high order polynomial interpolation. However,
the hydrodynamic features are captured very well by the high-order refinement.

###  V Summary and future work

In this work, we developed and tested an unstructured fully compressible
Navier Stokes solver which adaptively refines the mesh by increasing the local
order of refinement of the interpolation and derivative calculation. The
solver is based on the spectral difference technique which allows definition
of order of interpolation local to the mesh elements. We tested the solver on
three test cases. Results for one-dimensionl shock wave and a vortex advection
show very good agreement with analytical results. However, for the complex
test case of shock vortex interaction, we note that element boundaries exhibit
large oscillations of pressure. Future research from the developed solver
includes further development of adaptive element refinement ($h-$refinement)
which breaks a given element into subelements thus allowing higher resolution
with low order interpolation as well.

##  Chapter 8 Contributions, conclusions, and future research

In this chapter, major contributions and conclusions of this dissertation are
summarized. We also present recommendations for future research.

###  I Major contributions

####  I.1 Spectral energy dynamics in nonlinear acoustics

  * •

Second order scaling of maximum entropy generated within weak shock waves with
pressure jump for an ideal gas was shown.

  * •

Utilizing the second order scaling of maximum entropy, the correct second
order equations for planar nonlinear acoustics were derived which hold for any
general nonlinear acoustic system (within second order approximation), unlike
Burgers equation which holds only for a planar propagating wave.

  * •

The second order equations yielded the correct perturbation energy corollary
for nonlinear acoustics and the analytical expression of the correct
perturbation energy norm for nonlinear acoustics.

  * •

From the correct perturbation energy corollary, the spectral energy, spectral
energy flux, and spectral energy dissipation expressions were obtained and
analyzed utilizing DNS of fully compressible one-dimensional Navier-Stokes
equations.

####  I.2 Spectral energy cascade in thermoacoustic shock waves

  * •

First high fidelity numerical and analytical modeling of thermoacoustic shock
waves was done utilizing the fully compressible Navier-Stokes simulations,
pushing the thermoacoustic amplification of acoustic waves to the formation of
shock waves.

  * •

In the linear regime of thermoacoustic amplification of waves, thermoacoustic
energy budgets were analyzed yielding the optimum phasing required for
thermoacoustic growth.

  * •

In the nonlinear regime, novel framework of nonlinear acoustic waves
propagating near solid no-slip isothermal walls was developed, yielding the
macrosonic interaction of wall-heat flux and wall-shear on the nonlinear
acoustic waves.

  * •

Scaling of spectral acoustic energy was shown
($\widehat{E}_{k}\propto\omega^{-5/2}_{k}$), utilizing the computational
models and theoretical dimensionless scaling arguments.

####  I.3 Development of an order adaptive Navier-Stokes solver

  * •

Motivated by the study of spark plasma shock wave induced flow field and
sustained planar detonation wave dynamics, a high fidelity order adaptive
unstructured Navier-Stokes solver was developed, which can be further utilized
for simulating even more complex flow phenomena.

###  II Conclusions and recommendations for future research

####  II.1 Spectral energy cascade and decay in planar nonlinear acoustic
waves

  * •

The maximum entropy perturbations $s^{\prime}$ in an ideal gas scale as
$p^{\prime 2}$ for normalized pressure perturbation
$p^{\prime}\sim\mathcal{O}\left(10^{-3}-10^{-1}\right)$.

  * •

In traveling waves (TW) and standing waves (SW), spatial average of
perturbation energy decays as $\left\langle E^{(2)}\right\rangle\sim t^{-2}$
and dissipation rate as $\epsilon\sim t^{-3}$ in time. Moreover, the spectral
energy for both traveling and standing waves assumes the self-similar form:
$\widehat{E}_{k}k^{2}\epsilon^{-2/3}\ell^{1/3}\sim 0.075f(k\eta)$, where
$\ell$ is the integral length scale and $\eta$ is the Kolmogorov length scale.

  * •

In acoustic wave turbulence (AWT), due to gradual increase of the integral
length scale $\ell$ caused by the shock coalescence, the approximate decay
laws are $\left\langle E^{(2)}\right\rangle\sim t^{-2/3}$ and $\epsilon\sim
t^{-5/3}$.

  * •

The correct energy corollary can be utilized to evaluate nonlinear
interactions of high amplitude sound and heat release due to combustion or
thermal gradients.

  * •

Spectral energy cascade due to the nonlinear acoustic effects plays an
important role in the study of high Reynolds number compressible turbulence
Lele_AnnRev_1994 . Hence, importance of nonlinear compressible effects can be
quantified via spectral energy transfer terms derived in this work.

####  II.2 Spectral energy cascade in thermoacoustic shock waves

  * •

Three regimes of thermoacoustic wave amplification were 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.

  * •

The existence of an equilibrium _thermoacoustic energy cascade_ was shown. The
spectral energy density at the limit cycle, in particular, was found to decay
as $\omega^{-5/2}$ in spectral space, the relevant intensity scaling with
growth rate as $\alpha^{3}$.

  * •

Nonlinear wave propagation equations in the presence of wall-shear and wall-
heat transfer can be extended to higher dimensions establishing a generalized
framework for understanding nonlinear wave propagation and effects of wall-
shear and wall-heat transfer on the acoustic variables.

####  II.3 Mesh adaptive Navier-Stokes solver

  * •

Axisymmetric simulations of flow field behind a spark induced shock wave were
carried out, revealing the importance of baroclinic torque in the generation
of vorticity field behind the shock wave.

  * •

Utilizing the mesh adaptive solver, full three-dimensional simulations of the
flow field behind the spark-plasma induced shock wave can be carried out, thus
enabling the study of low Reynolds number turbulence observed in the
experiments BsinghEtAl_2019_AIAA .

  * •

Further development of the mesh adaptive solver (implementation of detailed
chemistry models) would enable the study of detonation wave propagation in
three dimensions and resonant dynamics sustained by fuel injection and
exhaust.

##  Chapter 9 Status of publications

###  I Refereed journal articles

  * •

Prateek Gupta, Guido Lodato, Carlo Scalo, _Spectral energy cascade in
thermoacoustic shock waves_ , _Journal of Fluid Mechanics_ , 2017, 831;
358-393

  * •

Danish Patel, Prateek Gupta, Carlo Scalo, _Acoustic Impedance Calculation via
Numerical Solution of the Inverse Helmholtz Problem_ , In review _Journal of
Sound and Vibration_ , 2018

  * •

Prateek Gupta, Carlo Scalo, _Spectral energy cascade and decay in nonlinear
acoustic waves_ , Submitted to _Physical Review E_ , 2018

###  II Presentations and non-refereed conference publications

  * •

Prateek Gupta, Guido Lodato, Carlo Scalo, _Nonlinear modeling of
thermoacoustically driven energy cascade_ , 69th APS DFD (Portland, Oregon,
USA), Nov. 2016

  * •

Prateek Gupta, Guido Lodato, Carlo Scalo, _Numerical investigation of
thermoacoustic shock waves_ , 55th AIAA Scitech (Grapevine, Texas, USA), Jan.
2017

  * •

Danish Patel, Prateek Gupta, Thomas Rothermel, Markus Kuhn, Carlo Scalo,
_mpedance Characterization of Ultrasonically Absorptive Coatings via Numerical
Solution of the Inverse Helmholtz Problem_ , 55th AIAA Scitech (Grapevine,
Texas, USA), Jan. 2017

  * •

Mario Tindaro Migliorino, Prateek Gupta, Carlo Scalo, _Real fluid effects on
thermoacoustic standing-wave resonance in supercritical $CO_{2}$_, 55th AIAA
Aviation (Denver, Colorado, USA), June 2017

  * •

Prateek Gupta, Carlo Scalo, _Glimpses of Kolmogorov energy cascade in
nonlinear acoustic waves_ , 70th APS DFD (Denver, Colorado, USA), Nov. 2017

  * •

Prateek Gupta, Kyle Schwinn, Guido Lodato, Carson Slabaugh, Carlo Scalo,
_Numerical investigation of sustained planar detonation waves in a periodic
domain_ , 56th AIAA Aviation (Atlanta, Georgia, USA), June 2018

##  Chapter 10 Derivation of second order equations

Decomposing the variables in base state and perturbation fields, we obtain,

$\displaystyle\rho^{*}=\rho^{*}_{0}+{\rho^{*}}^{\prime},\quad
p^{*}=p^{*}_{0}+{p^{*}}^{\prime},$ (10.1a) $\displaystyle
u^{*}={u^{*}}^{\prime},\quad s^{*}={s^{*}}^{\prime},\quad
T^{*}=T^{*}_{0}+{T^{*}}^{\prime},$ (10.1b)

where the superscript $()_{0}$ denotes the base state, the superscript
$({})^{\prime}$ denotes the perturbation values, no mean flow $u^{*}_{0}=0$ is
considered, and $s^{*}_{0}$ is arbitrarily set to zero. We neglect the
fluctuations in the dynamic viscosity as well, i.e.,

$\mu^{*}=\mu^{*}_{0}.$ (10.2)

While in classic gas dynamics, pre-shock values are used to normalize
fluctuations or jumps across the shock (e.g. see Eq. (2.7)), here we choose
base state values to non-dimensionalize the variables,

$\displaystyle\rho=\frac{{\rho^{*}}}{\rho^{*}_{0}}=1+\rho^{\prime},~{}~{}p=\frac{p^{*}}{\gamma
p^{*}_{0}}=\frac{1}{\gamma}+p^{\prime},$ (10.3a) $\displaystyle
u=\frac{u^{*}}{a^{*}_{0}}=u^{\prime},~{}~{}s=\frac{s^{*}}{R^{*}}=s^{\prime},~{}~{}T=\frac{T^{*}}{T^{*}_{0}}=1+T^{\prime},$
(10.3b) $\displaystyle
x=\frac{x^{*}}{L^{*}},~{}~{}t=\frac{a^{*}_{0}t^{*}}{L^{*}}.$ (10.3c)

where $L^{*}$ is the length of the one-dimensional periodic domain. For a
chemically inert generic gas, infinitesimal changes in dimensionless density
$\rho(p,s)$ in terms of pressure $p$ and entropy $s$ are given by,

$\displaystyle d\rho$ $\displaystyle=\left(\frac{\partial\rho}{\partial
p}\right)_{s}dp+\left(\frac{\partial\rho}{\partial s}\right)_{p}ds,$
$\displaystyle=\frac{\rho}{\gamma
p}dp-\left(\frac{\rho^{*}_{0}T^{*}_{0}R^{*}}{\gamma
p^{*}_{0}}\right)\frac{\rho^{2}T}{p}\left(\frac{\gamma-1}{\gamma}\right)ds.$
(10.4)

Substituting the above relation in the dimensionless continuity Eq. (2.15), we
obtain,

$\displaystyle\frac{\partial p}{\partial t}$ $\displaystyle+u\frac{\partial
p}{\partial x}+\gamma p\frac{\partial u}{\partial x}$
$\displaystyle=\left(\frac{\rho^{*}_{0}T^{*}_{0}R^{*}}{p^{*}_{0}}\right)\left(\frac{\gamma-1}{\gamma}\right){\rho
T}\left(\frac{\partial s}{\partial t}+u\frac{\partial s}{\partial x}\right).$
(10.5)

Non-dimensionalizing the entropy Eq. (2.3) utilizing the Eq. (10.3), we
obtain,

$\rho T\left(\frac{\partial s}{\partial t}+u\frac{\partial s}{\partial
x}\right)=\frac{\nu_{0}}{Pr}\frac{C^{*}_{p}}{R^{*}}\frac{\partial^{2}T}{\partial
x^{2}}+\frac{4\nu_{0}}{3}\frac{a^{*2}_{0}}{R^{*}T^{*}_{0}}\left(\frac{\partial
u}{\partial x}\right)^{2}.$ (10.6)

Substituting the above equation in Eq. (10.5), we obtain,

$\displaystyle\frac{\partial p}{\partial t}+u\frac{\partial p}{\partial
x}+\gamma p\frac{\partial u}{\partial x}$
$\displaystyle=\left(\frac{\rho^{*}_{0}T^{*}_{0}R^{*}}{p^{*}_{0}}\right)\frac{\gamma-1}{\gamma}\left(\frac{\nu_{0}}{Pr}\frac{C^{*}_{p}}{R^{*}}\frac{\partial^{2}T}{\partial
x^{2}}+\frac{4\nu_{0}}{3}\frac{a^{*2}_{0}}{R^{*}T^{*}_{0}}\left(\frac{\partial
u}{\partial x}\right)^{2}\right).$ (10.7)

Substituting the decomposition of variables (cf. Eq. (10.3)) in the above Eq.
(10.7), we obtain the pressure perturbation equation for a generic gas,

$\displaystyle\frac{\partial p^{\prime}}{\partial t}+\frac{\partial
p^{\prime}}{\partial x}+u^{\prime}\frac{\partial p^{\prime}}{\partial
x}+\gamma p^{\prime}\frac{\partial u^{\prime}}{\partial x}$
$\displaystyle=\left(\frac{\rho^{*}_{0}T^{*}_{0}R^{*}}{p^{*}_{0}}\right)\frac{\gamma-1}{\gamma}\left(\frac{\nu_{0}}{Pr}\frac{C^{*}_{p}}{R^{*}}\frac{\partial^{2}T^{\prime}}{\partial
x^{2}}+\frac{4\nu_{0}}{3}\frac{a^{*2}_{0}}{R^{*}T^{*}_{0}}\left(\frac{\partial
u^{\prime}}{\partial x}\right)^{2}\right).$ (10.8)

As shown in Section II.2, the entropy perturbations are atmost
$2^{\mathrm{nd}}$ order in pressure, independent of viscosity. Consequently,
the first and second term on right hand side of Eq. (10.6) are second and
third order in pressure perturbations, respectively. Truncating the Eq. (10.8)
up to second order, we obtain the second order equation for pressure
perturbations for a generic fluid as,

$\displaystyle\frac{\partial p^{\prime}}{\partial t}+u^{\prime}\frac{\partial
p^{\prime}}{\partial x}+\frac{\partial u^{\prime}}{\partial x}+\gamma
p^{\prime}\frac{\partial u^{\prime}}{\partial x}=$
$\displaystyle\frac{\nu_{0}}{Pr}\left(\frac{\rho^{*}_{0}T^{*}_{0}R^{*}}{p^{*}_{0}}\right)\left(\frac{\gamma-1}{\gamma}\right)\left(\frac{\partial
T}{\partial
p}\right)_{s,0}\frac{C^{*}_{p}}{R^{*}}\frac{\partial^{2}p^{\prime}}{\partial
x^{2}}$ $\displaystyle+\mathcal{O}\left(p^{\prime}s^{\prime},s^{\prime