# Theoretical and Numerical Investigation
of Nonlinear Thermoacoustic, Acoustic,
and Detonation Waves

Prateek Gupta

###### Abstract

Finite amplitude perturbations in compressible media are ubiquitous in
scientific and engineering applications such as gas-turbine engines, rocket
propulsion systems, combustion instabilities, inhomogeneous solids, and
traffic flow prediction models, to name a few. Small amplitude waves in
compressible fluids propagate as sound and are very well described by linear
theory. On the other hand, the theory of nonlinear acoustics, concerning high-
amplitude wave propagation (Mach$<$2) is relatively underdeveloped. Most of
the theoretical development in nonlinear acoustics has focused on wave
steepening and has been centered around the Burgers’ equation, which can be
extended to nonlinear acoustics only for purely one-way traveling waves. In
this dissertation, theoretical and computational developments are discussed
with the objective of advancing the multi-fidelity modeling of nonlinear
acoustics, ranging from quasi one-dimensional high-amplitude waves to
combustion-induced detonation waves.

We begin with the theoretical study of spectral energy cascade due to the
propagation of high amplitude sound in the absence of thermal sources. To this
end, a first-principles-based system of governing equations, correct up to
second order in perturbation variables is derived. The exact energy corollary
of such second-order system of equations is then formulated and used to
elucidate the spectral energy dynamics of nonlinear acoustic waves. We then
extend this analysis to thermoacoustically unstable waves – i.e. amplified as
a result of thermoacoustic instability. We drive such instability up until the
generation of shock waves. We further study the nonlinear wave propagation in
geometrically complex case of waves induced by the spark plasma between the
electrodes. This case adds the geometrical complexity of a curved, three-
dimensional shock, yielding vorticity production due to baroclinic torque.
Finally, detonation waves are simulated by using a low-order approach, in a
periodic setup subjected to high pressure inlet and exhaust of combustible
gaseous mixture. An order adaptive fully compressible and unstructured Navier
Stokes solver is currently under development to enable higher fidelity studies
of both the spark plasma and detonation wave problem in the future.

Prateek Gupta Doctor of PhilosophyPhDAugust2019 Carlo Scalo West Lafayette

Dr. Carlo Scalo, ChairSchool of Mechanical Engineering Dr. Ivan ChristovSchool
of Mechanical Engineering Dr. Stuart BoltonSchool of Mechanical Engineering
Dr. Sally BaneSchool of Aeronautics and Astronautics Dr. Jay P. GoreAssociate
Head for Graduate Studies

###### Acknowledgements.

I would like to express my sincere gratitude to my advisor Prof. Carlo Scalo
for his help and support. His passion to research upon a fundamental idea all
the way up to application and strong advocacy of scientific rigor have been
constant motivations. I would also like to thank Prof. Ivan Christov, Prof.
Sally Bane, and Prof. Stuart Bolton for serving on my defense committee and
for their feedback on this dissertation. I wish to thank my colleagues at
Compressible Flow and Acoustics Lab for creating a vibrant atmosphere in the
workplace. A special thanks to my colleagues and friends, Zongxin Yu, Yongkai
Chen, Emmanuel Gill Torres, Karl Jantze, Julien Brillon, and Joel Redmond for
their company. A very special thanks to my friends at Purdue, including Kumar
Akash, Kunal Pardikar, Akash Patil, Yeshaswi Menghmalani, Rohil Jain, Viplove
Arora, Vishwanath Ganesan, and Vishal Anand. I would also like to thank Prof.
Guido Lodato at INSA de Rouen for his guidance regarding the SD3DvisP solver
in my early years of PhD. Besides people directly involved in my PhD, I am
extremely grateful to my bachelor thesis advisor, Prof. Supreet Singh Bahga at
IIT Delhi, Mechanical Engineering Department, for his support during tough
times. His guidance played a pivotal role in conditioning me into a patient
and thorough researcher. I am grateful to my parents, brother, sister-in-law,
and the extended family for their continued support and unconditional love.
This dissertation is a culmination of their support, time, and energy put into
me. To them, I dedicate this dissertation.

###### Contents

  1. 1 Introduction
  2. 2 Planar second order nonlinear acoustics
    1. I Introduction
    2. II Governing equations and scaling analysis
      1. II.1 Fully compressible 1D Navier-Stokes and entropy scaling in weak shocks
      2. II.2 Perturbation variables and non-dimensionalization
      3. II.3 Second order nonlinear acoustics equations
  3. 3 Energy cascade and decay of nonlinear acoustic waves
    1. I Introduction
    2. II Second order perturbation energy
    3. III High Fidelity Simulations with Adaptive Mesh Refinement
      1. III.1 Numerical approach
      2. III.2 Initial conditions
    4. IV Scales of acoustic energy cascade and dissipation
      1. IV.1 Spectral energy flux and dissipation rate for periodic perturbations
      2. IV.2 Taylor microscale
      3. IV.3 Integral length scale
      4. IV.4 Kolmogorov length scale
    5. V Scaling of spectral quantities
    6. VI Concluding remarks
  4. 4 Thermoacoustic nonlinear energy cascade
    1. I Introduction
    2. II Problem Formulation
      1. II.1 Design of the minimal-unit model
      2. II.2 Navier-Stokes calculations
    3. III Regimes of Thermoacoustic Amplification
      1. III.1 Spectral energy density
      2. III.2 Regime classification
    4. IV Harmonic growth analysis
      1. IV.1 Linear thermoviscous quasi-planar wave equations
      2. IV.2 Acoustic energy budgets for quasi-planar wave perturbations
      3. IV.3 Effects of temperature gradient and geometry on growth rates
    5. V Formulation of a Nonlinear Thermoacoustic Model
      1. V.1 Governing equations for nonlinear thermoviscous perturbations
      2. V.2 Free-shear nonlinear wave propagation
      3. V.3 Wall-shear and wall-heat flux
      4. V.4 Nonlinear thermoviscous wave equations
    6. VI Nonlinear Spectral Energy Dynamics
    7. VII Scales of thermoacoustically sustained spectral energy cascade
    8. VIII Summary
  5. 5 Spark-plasma generated shock waves
    1. I Introduction
    2. II Experimental study and motivation (Credit: Ms. Bhavini Singh, Mr. Lalit Rajendran, and Dr. Sally Bane)
    3. III Computational setup
    4. IV Results
  6. 6 Numerical modeling of detonation waves
    1. I Introduction
    2. II Experimental study and motivation (Credit: Mr. Kyle Schwinn and Dr. Carson Slabaugh)
    3. III Numerical validation
      1. III.1 Validation with Shock-to-detonation transition
      2. III.2 Pre-mixed laminar flame
    4. IV One-dimensional problem formulation
    5. V Results
    6. VI Conclusions
  7. 7 Development of an order adaptive Navier Stokes solver
    1. I Introduction
    2. II Spectral difference methods : Overview
      1. II.1 Governing equations
      2. II.2 Mapping and spatial discretization
    3. III Adaptive order refinement
    4. IV Test cases
      1. IV.1 1D Shock tube
      2. IV.2 Convecting isentropic vortex
      3. IV.3 Shock Vortex interaction
    5. V Summary and future work
  8. 8 Contributions, conclusions, and future research
    1. I Major contributions
      1. I.1 Spectral energy dynamics in nonlinear acoustics
      2. I.2 Spectral energy cascade in thermoacoustic shock waves
      3. I.3 Development of an order adaptive Navier-Stokes solver
    2. II Conclusions and recommendations for future research
      1. II.1 Spectral energy cascade and decay in planar nonlinear acoustic waves
      2. II.2 Spectral energy cascade in thermoacoustic shock waves
      3. II.3 Mesh adaptive Navier-Stokes solver
  9. 9 Status of publications
    1. I Refereed journal articles
    2. II Presentations and non-refereed conference publications
  10. 10 Derivation of second order equations
  11. 11 Derivation of Spectral Energy Transfer
  12. 12 Experimental Validation of linear model
  13. 13 Scaling of governing equations
  14. 14 Time domain model and convergence
  15. 15 Windowed shock capturing

###### List of Tables

  1. 3.1 Simulation parameter space for TW, SW, and AWT cases listing base state dimensionless viscosity $\nu_{0}$ (cf. Eq. (2.18)), initial characteristic perturbation amplitude $A_{\mathrm{rms},0}$ (cf. Eq. (3.16)), and dimensional characteristic perturbation in velocity $u^{*}_{\mathrm{rms}}$ and pressure $p^{*}_{\mathrm{rms}}$ fields (Eq. (3.17)).
  2. 3.2 Initial spectral compositions for traveling wave (TW), standing wave (SW), and acoustic wave turbulence (AWT). $\delta\left(\cdot\right)$ is the Dirac delta function.
  3. 3.3 Summary of the three length scales $\ell$, $\lambda$, and $\eta$, respective definitions, and the range of spectrum characterized by them. The integral length scale characterizes the energy containing range $(k_{0},k_{E})$. The Taylor microscale is the characteristic of the energy transfer and dissipation range $(k_{E},k_{\delta})$. The Kolmogorov length scale corresponds to the highest wavenumber generated as a result of nonlinear acoustic energy cascade.
  4. 4.1 Thermodynamic parameters for base state. Base pressure and cold-side values of density and temperature are set to reference values: $P_{0}=P_{\mathrm{ref}}$, $\rho_{c}=\rho_{\mathrm{ref}}$, $T_{C}=T_{\mathrm{ref}}$.
  5. 6.1 Injection stagnation pressures $p_{0,\mathrm{inj}}$ and the weighting factors $w_{i}$ considered.

###### List of Figures

  1. 1 Introduction
  2. 1.2 Physical setup of a looped thermoacoustically unstable resonator.
  3. 1.3 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.
  4. 1.4 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.
  5. 1.5 Advection of an isentropic vortex in a subsonic base flow (top) and adaptive order refinement of the underlying mesh (bottom).
  6. 2.1 Variation of entropy and pressure within a weak shock wave for varying viscosity values.
  7. 2.2 Scaling of entropy jump across the shock wave and maximum entropy inside the shock wave with pressure jump.
  8. 3.1 Comparison of perturbation energy function for nonlinear acoustic waves $E^{(2)}$ with the linear acoustic energy $E^{(1)}$ in the case of $p^{\prime}=u^{\prime}$ (assumed for illustrative purpose). The correction $f(p^{\prime})$ is independent of $u^{\prime}$.
  9. 3.2 Spatial profile of finite amplitude waves for a traveling wave, a standing wave, and an acoustic wave turbulence case and temporal evolution of the average perturbation energy ($\left\langle E^{(2)}\right\rangle$ (–); $\left\langle E^{(1)}\right\rangle$ ($--$)) evaluated from the DNS data
  10. 3.3 Illustration of the binary tree implementation of Adaptive Mesh Refinement (AMR) technique.
  11. 3.4 Schematic illustrating the comparison of various length scales associated with spectral energy cascade in nonlinear acoustics in both spatial $(a)$ and spectral $(b)$ space. $(a)$ shows the perturbation velocity $u^{\prime}$ ($--$) and pressure $p^{\prime}$ (–) fields in AWT obtained from the DNS data for $\nu_{0}=$1.836$\times$10-7 and $A_{\mathrm{rms}}=10^{-1}$. $(b)$ shows the corresponding spectral energy $\widehat{E}_{k}$ in log-log space. The integral length scale $\ell$ corresponds to the characteristic distance between the shock waves traveling in the same direction. The Kolmogorov length scale $\eta$ corresponds to the shock wave thickness. The Taylor microscale $\lambda$ is the diffusive length scale and satisfies $\ell\gg\lambda\gg\eta$. $L$ corresponds to the length of the domain.
  12. 3.5 Schematic illustrating the comparison of various length scales associated with spectral energy cascade in nonlinear acoustics in the spectral space, and corresponding variation of spectral energy flux $\hat{\Pi}_{k}$ ($a$) and spectral dissipation $\mathcal{D}_{k}$ $(b)$. The integral length scale $\ell$ corresponds to the characteristic distance between the shock waves traveling in the same direction. The Kolmogorov length scale $\eta$ corresponds to the shock wave thickness. The Taylor microscale $\lambda$ is the diffusive length scale and satisfies $\ell\gg\lambda\gg\eta$. $L$ corresponds to the length of the domain.
  13. 3.6 Spectro-temporal evolution of $\widehat{E}_{k}$ (top) and spectral flux $\widehat{\Pi}_{k}$ (bottom) for a TW ($a$), SW ($b$) and AWT ($c$). Spectral flux $\widehat{\Pi}_{k}$ for a traveling wave simply increases towards high wavenumbers. For a standing wave, $\widehat{\Pi}_{k}$ oscillates at low wavenumbers cyclically due to collisions of oppositely traveling shock waves while high wavenumber behaviour resembles that of a traveling shock. For AWT, the spectral broadening occurs for $k>k_{E}$ with small fluctuations in time for $k<k_{E}$.
  14. 3.7 Temporal evolution of scaled total energy <$E^{(2)}$> ${A^{-2}_{\mathrm{rms},0}}$ (–) (top), dissipation rate $\epsilon{A^{-3}_{\mathrm{rms},0}}$ ($--$) (mid) and normalized Taylor microscale $\lambda/\sqrt{\delta\tau}$ (bottom) for TW ($a$), SW ($b$) and AWT ($c$) against the scaled time $t/\tau$ for varying perturbation Reynolds number $\mathrm{Re}_{L}$. The time $t_{0}$ signifies fully broadened spectrum of the perturbation field.
  15. 3.8 Evolution of the integral length scale $\ell$ $(a)$ and the Reynolds number $\mathrm{Re}_{\ell}$ ($b$) defined in Eqs. (3.36) and (3.45), respectively, for all the cases of AWT considered. For small thermoviscous diffusivity, $\ell$ increases approximately as $t^{1/3}$ before saturating to the dimensionless domain length $L=1$ and $\mathrm{Re}_{\ell}$ remains approximately constant.
  16. 3.9 Fully developed spectra of compensated energy $(a)$, spectral energy flux $(b)$, and cumulative dissipation $(c)$ for TW at time instant $t_{0}\approx 0.03$. Harmonics with wavenumbers such that $k\eta<1$ contain all the energy. The spectral energy flux vanishes at $k\eta\approx 1$ thus indicating numerical resolution of all the energy containing harmonics. The marked regime $0.1<k\eta<1$ signifies the dissipation range. The constant $C\approx 0.075$. (–) $A_{\mathrm{rms},0}=10^{-1}$; ($--$) $A_{\mathrm{rms},0}=10^{-2}$; ($\cdots$) $A_{\mathrm{rms},0}=10^{-3}$
  17. 3.10 Fully developed spectra of compensated energy $(a)$, spectral energy flux $(b)$, and cumulative dissipation $(c)$ for SW averaged over one time cycle after $t_{0}\approx 0.04$. Harmonics with wavenumbers $k\eta<1$ contain all the energy. The spectral energy flux vanishes at $k\eta\approx 1$ thus indicating numerical resolution of all the energy containing harmonics. The constant $C\approx 0.075$. (–) $A_{\mathrm{rms},0}=10^{-1}$; ($--$) $A_{\mathrm{rms},0}=10^{-2}$; ($\cdots$) $A_{\mathrm{rms},0}=10^{-3}$
  18. 3.11 Fully developed spectra of compensated energy, $(a)$, spectral energy flux $(b)$, and cumulative dissipation $(c)$ against scaled wavenumber $k\eta$ for the randomly initialized broadband noise (AWT) cases with $A_{\mathrm{rms},0}$ and $\nu_{0}$ listed in Table 3.1 at dimensionless time $t/\tau=\tau_{0}\approx 6\times 10^{-4}$. The marked regime $0.1<k\eta<1$ signifies the dissipation range. (–) $A_{\mathrm{rms},0}=10^{-1}$; ($--$) $A_{\mathrm{rms},0}=10^{-2}$; ($\cdots$) $A_{\mathrm{rms},0}=10^{-3}$
  19. 4.1 Two-dimensional axially periodic computational setup for minimal-unit simulations (top, not to scale), geometrical parameters (bottom left), and qualitative illustration of the equivalent experimental build-up of a variable-area looped thermoacoustic resonator (bottom right). Top figure: (– –), adiabatic/sip conditions; (–), isothermal/no-slip conditions. The minimal unit is traced in the bottom right with dashed lines (– –). The listed geometrical parameters provide sub-optimal growth rates across all the values of $T_{H}$, the boxed quantities have been determined via optimization (cf. figure  and table 4.1).
  20. 4.3 Illustration of selection portions of the unstructured mesh near duct interfaces for Grid A with total number of elements $N_{\mathrm{el}}=2004$.
  21. 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.
  22. 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}$.
  23. 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. ()). ($a$): (–), $\Phi^{T}_{\mathcal{P}}$; (- -), $-\Phi^{S}_{\mathcal{P}}$.
  24. 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}$.
  25. 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.
  26. 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*}$.
  27. 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}$.
  28. 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.
  29. 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.
  30. 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).
  31. 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).
  32. 5.2 Shock wave induced due to the spark plasma heat deposition.
  33. 5.3 Computational setup of the electrodes. The setup models axisymmetric flow conditions.
  34. 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.
  35. 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.
  36. 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.
  37. 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.
  38. 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.
  39. 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 )
  40. 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 )
  41. 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 )
  42. 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}$.
  43. 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)$
  44. 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.
  45. 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$.
  46. 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.
  47. 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$).
  48. 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}$.
  49. 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 ($--$).
  50. 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}$).
  51. 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.
  52. 7.3 Mortar element connecting elements with different local orders with non-conforming flux points.
  53. 7.4 Spatio-temporal evolution of the one dimensional Sod shock tube problem and the adaptive order refinement.
  54. 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$.
  55. 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.
  56. 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).
  57. 12.1 Neutral stability curves obtained from linear analysis (§ IV), Yazaki_PhysRevLet_1998 ’s data ($a$), and BiwaTY_JAcousSocAm_2011 ’s reported instability limit ($b$). $r$ is the hydraulic radius of regenerator pore (square cross-section). ($a$): (–), Y-1998 in 3D ; (- -), 3D minimal unit corresponding to Y-1998; ($\cdot\cdot\cdot$), guedra_use_2012 ’s results; ($\circ$), Yazaki_PhysRevLet_1998 ’s reported data. ($b$): (–), B-2011 in 3D (with $\ell_{buffer}=0.4~{}$m); (- -), 3D minimal unit corresponding to B-2011 without buffer length ; ($\circ$), reported instability limit.

##  Chapter 1 Introduction

Nonlinear wave processes are observed in a variety of engineering and physics
applications such as acoustics Hamilton_NLA_1998 ; Lighthill_JSV_1978 ;
Whitham2011 , thermoacoustics ScaloLH_JFM_2015 ; GuptaLS_JFM_2017 ;
Culik_AGARD_2006 ; Mcmanus1993review ; Schwinn_CandF_2018 ;
PoinsotV_numComb_2011 , surface waves AblowitzClarkson1991 , plasma-physics
Gurbatov2012 , and inhomogeneous solids Nazarov2015 . In acoustics, nonlinear
wave propagation (finite amplitude acoustics) has been an active area of
research since more than five decades now and can be classified as both
hyperbolic and dispersive, based on the medium of propagation Whitham2011 .
Planar nonlinear waves in compressible fluids primarily exhibit non-dispersive
hyperbolic nature. Propagation of such waves causes two main nonlinearities
known, acoustic streaming Lighthill_JSV_1978 ; Gedeon_1997_Cryocoolers and
wave steepening Hamilton_NLA_1998 ; Naugolnykh1998 . While acoustic streaming
is a kinematic nonlinear effect, wave steepening is caused by the variation of
speed of propagation (local speed of sound) due to variation in pressure. For
very high pressure variations, nonlinear acoustic waves form shock waves,
which can also cause combustion of the medium due to intense temperature
variation, depending on the ignition properties of the medium. Moreover, non-
planar nonlinear waves can cause intense hydrodynamic events such as vorticity
generation Truesdell_1952_JAS , which can be eventually utilized for changing
scalar and momentum mixing properties of flows.

The scope of the present work includes theoretical and numerical study of such
nonlinear wave propagation effects, in particular, mathematically quantifying
the planar nonlinear acoustic waves identifying the extent of second order
nonlinear theory, effect of boundary heating and combustion heating on such
nonlinear waves, followed by the hydrodynamic effects of such waves on flows.
We begin with the basic definition and mathematical identification of such
nonlinearities. As the high amplitude planar waves propagate, multiple length
scales are generated which modify the energy dissipation dynamics
significantly. Moreover, such nonlinear waves can be generated from a variety
of heat sources such as thermoacoustic instability, heat release in detonation
waves, spark plasma discharges etc. We outline the analysis of such nonlinear
wave-heat source interactions for compressible fluids. Throughout, theoretical
developments are supported by high-fidelity canonical numerical simulations.

This dissertation is organized in five chapters, each summarized briefly
below.

$(a)$$(b)$$(c)$$(d)$$(e)$$(f)$

### Summary of Chapters 2 and 3: Decaying nonlinear acoustics

Nonlinear acoustic wave steepening occurs due to gradients in the wave speed
associated with thermodynamic nonlinearities, and it entails generation of
smaller length scales (harmonics) via a nonlinear energy cascade, which can be
realized mathematically via the multiplication of two truncated Fourier
series,

$\displaystyle\left(\sum^{n}_{k=-n}a_{k}e^{ikx}\right)\left(\sum^{m}_{l=-m}b_{l}e^{ilx}\right)=={\sum_{k}a_{k}b_{-k}}+{\underset{k+l\neq
0}{\sum_{k}\sum_{l}}a_{k}b_{l}e^{i(k+l)x}}.$ (1.1)

The left hand side of the Eq. (1.1) represents any quadratic nonlinear term
appearing in a governing equation. Continued nonlinear evolution results in
further generation of smaller length scales, as depicted by the second term on
the right hand side of Eq. (1.1), ultimately leading to spectral broadening,
as shown in figure LABEL:fig:_Introfig1. In the specific case of steepening of
nonlinear acoustic waves, the shock thickness denotes the smallest length
scale generated in the flow, governed by the viscous dissipation. The latter
causes saturation of energy spectra, hence establishing an energy flow among
various scales.

In this work, we focus on the characterization of the spectral energy cascade
in nonlinear acoustics. In particular, we study the spatio-temporal and
spectro-temporal evolution of decaying finite amplitude planar nonlinear
acoustic waves in three canonical configurations, traveling waves (TW),
standing waves (SW), and randomly initialized Acoustic Wave Turbulence (AWT).
Utilizing the second order nonlinear acoustics approximation, we derive
analytical expressions for the spectral energy, energy transfer function, and
dissipation. Analogous to the study of small scale generation in hydrodynamic
turbulence, well quantified by the K41 theory kolmogorov1941a ;
kolmogorov1941b ; kolmogorov1941c , we define relevant length scales
associated with fully developed nonlinear acoustic waves elucidating the
scaling features of the energy spectra.

### Summary of Chapter 4: Thermoacoustic shock waves

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

Figure 1.2: Physical setup of a looped thermoacoustically unstable resonator
$(a)$, pressure probed inside the resonator $(b)$, and the isentropic acoustic
energy spectrum of the probed signal $(c)$ depicting spectral energy cascade.

The acoustic energy cascade in the spectral space due to nonlinear acoustic
wave propagation is analogous to the turbulent spectral energy cascade in
which vortex stretching (hydrodynamic nonlinearity) causes breakdown of large
length scale eddies into smaller eddies up to a point where the viscous
dissipation dominates and prevents generation of further smaller scales
pope2000turbulent . In equilibrium energy cascade, mean shear causes injection
of turbulent kinetic energy into large eddies which breakdown into smaller
eddies thus establishing an equilibrium energy spectrum tennekes1972first .
Usually, hydrodynamic instabilities are associated as the mechanisms of energy
injection into large length scales by mean shear, such as Tollmien-Schlichting
instability in low speed hydrodynamic boundary layers schlichting1960boundary
. Analogously, acoustic energy cascade due to nonlinear acoustic wave
propagation can be sustained via acoustic instabilities. We show such
existence of equilibrium spectral energy cascade in which energy in large
length scales (harmonics) is injected by thermoacoustic instabilities
resulting in formation of thermoacoustically sustained shock waves (see figure
1.2).

Thermoacoustic amplification of waves in a compressible flow is the result of
a fluid dynamic instability emerging from the favourable coupling between
pressure and heat-release fluctuations. The wavelength of thermoacoustically
unstable waves is set by the size of the enclosing resonant chamber, while the
heat release, providing the energy source for the amplification, is confined
to a compact region. The heat release rate is a function of local velocity and
pressure fluctuations, affecting for example the instantaneous flame surface
area in a combustion chamber or the rate of convective heat extraction from a
hot wire-mesh screen in a Rijke tube. The resulting fluctuations in the heat
release rate drive a cycle of compressions and dilatations, which act as a
source for pressure fluctuations that, if within a quarter phase from the heat
release itself, become thermoacoustically amplified, as identified by the
Rayleigh’s criterion Rayleigh_Nature_1878 ),

$\int_{\Omega}p^{\prime}q^{\prime}d\Omega>0.$ (1.2)

Variables $p^{\prime}$ and $q^{\prime}$ denote the pressure and heat release
rate perturbations in a confined domain $\Omega$.

In this work, we show the application of the previously discussed theory on
nonlinear acoustic energy cascade. We investigate the high amplitude (or
macrosonic) limit of thermoacoustically driven nonlinear waves characterized
by the formation of self-sustaining resonating shock waves and inter-scale
energy transfer. A comprehensive nonlinear theoretical and high-fidelity
modeling approach is adopted to accurately describe macrosconic thermoacoustic
waves. To this end, a canonical travelling-wave looped resonator, inspired by
Yazaki _et al._ ’s Yazaki_PhysRevLet_1998 experimental setup but geometrically
optimized via linear theory Rott_ZAMP_1969 ; Rott_ZAMP_1973 ; LinSH_JFM_2016 ,
has been designed to maximize the growth rate of the quasi-travelling-wave
second harmonic and thus achieve rapid shock wave formation. Yazaki _et al._
’s Yazaki_PhysRevLet_1998 looped configuration allows quasi-travelling-wave
acoustic phasing which facilitates faster nonlinear energy cascade compared to
standing wave resonators BiwaEtAl_JASA_2014 . It is shown that the energy
content in spectral domain resembles the equilibrium energy cascade observed
in turbulence, similar to the spectral energy distribution of an ensemble of
acoustic waves interacting nonlinearly among each other Nazarenko_2011_WT ;
Zakharov_2012_WT .

### Summary of Chapter 5: Spark Plasma induced shock waves

Figure 1.3: 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.

For planar weak shocks, the acoustic nonlinearities constitute the primary
nonlinearities in the flow. However, for weak (or strong) shocks in 2D or 3D,
more interesting features appear. One such interesting feature is the
generation of vorticity behind a curved shock. As first derived by Truesdell
Truesdell_1952_JAS , a steady 2D curved shock induces a vorticity jump across
the shock directly proportional to the curvature of the shock. Later, this
calculation was generalized and analyzed by Kevlahan Kevlahan_1997_JFM for an
unsteady 2D shock, who showed the baroclinic moment applied by the curved
shock also causes vorticity generation.

In general, the baroclinic moment acting on the flow is given by,

$\Gamma=\frac{\nabla p\times\nabla\rho}{\rho^{2}}.$ (1.3)

For an isentropic perturbation, the pressure gradient $\nabla p$ and density
gradient $\nabla\rho$ are parallel. Consequently, isentropic waves can not
apply baroclinic moments on the fluid. However, waves generated due to intense
heat release may exhibit entropy gradients (even though perturbations are
weak) not aligned with pressure perturbations. Such entropy variations result
in baroclinic moments which cause vorticity changes in a flow.

In this work, we analyze the flow field generated by the heat deposition
between two electrodes due to dielectric breakdown. Such setups are utilized
in controlling flow fields and combustion mixing due to coherent structures
generated near the region of dielectric breakdown. We show that the baroclinic
moment plays an important role in generating coherent structures near conical
shaped electrodes. We model the spark plasma heat deposition be a spato-
temporally varying heating kernel.

### Summary of Chapter 6: High pressure detonation waves

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

$(a)$$(b)$

Figure 1.4: 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.

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. We numerically investigate the effect of fuel-oxidizer
injection rates on the dynamics of detonation waves and sustenance in a
periodic one-dimensional domain. To this end, 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.

### Summary of Chapter 7: Mesh adaptive Navier Stokes solver

In the previous two topics investigated, an extremely wide spectrum of length
scales is observed. In vortex breakdown behind a curved shock wave induced by
spark-plasma heat deposition, vortex rings have been observed in experiments
BsinghEtAl_2019_AIAA , which are generated by the curved shock wave. These
vortex rings collide and break down generating very small scale structures in
the flow. For detonation wave resonance in a periodic domain, the detonation
waves propagate in a domain $~{}10^{6}$ times longer than the thickness
Schwer_AIAA_2010 ; Schwer_CandF_2011 . As a result, very high numerical
resolution is required in only a very small region inside the domain which
mandates the mesh adaptive approach in multiple dimensions for efficient
computations.

To this end, a high fidelity mesh adaptive Navier-Stokes solver is developed
and presented in this work, as a tool for future development and study of the
topics discussed above. The solver adapts the local order of interpolation and
derivative calculation upon detecting shock waves and vorticity, thus enabling
higher degrees of freedom locally for higher resolution. Test results on cases
involving high gradients in flow fields due to shock waves and vorticity (see
figure 1.5) are presented.

Figure 1.5: Advection of an isentropic vortex in a subsonic base flow (top)
and adaptive order refinement of the underlying mesh (bottom).

##  Chapter 2 Planar second order nonlinear acoustics

The contents of this chapter have been published in _Physical Review E_
GuptaScalo_PRE_2018 (see Section II) and have been reported here in abridged
form with minor modifications.

###  I Introduction

High amplitude planar nonlinear acoustic waves exhibit two main nonlinearities
:- acoustic streaming Lighthill_JSV_1978 ; Gedeon_1997_Cryocoolers and wave
steepening Hamilton_NLA_1998 ; Naugolnykh1998 . Usually, problems in
nondispersive nonlinear wave propagation are studied utilizing the model
Burgers equation Whitham2011 ; Gurbatov_1981_JETP ; Naugolnykh1998 ;
Gurbatov2012 . However, in this work, we focus on quantifying the spatio-
temporal and spectro-temporal evolution of finite amplitude planar nonlinear
acoustic waves exactly via second order nonlinear governing equations, as
derived from the continuum gas dynamics (compressible Navier Stokes) equations
in this chapter. In later chapters, we study the evolution of finite amplitude
planar nonlinear acoustic waves in three canonical configurations, traveling
waves (TW), standing waves (SW), and randomly initialized Acoustic Wave
Turbulence (AWT) utilizing these equations.

###  II Governing equations and scaling analysis

In this section, we derive the governing equations for nonlinear acoustics
truncated up to second order (in the acoustic perturbation variables) for a
single-component ideal gas. We begin with fully compressible one-dimensional
Navier-Stokes equations for continuum gas dynamics and analysis of entropy
scaling with pressure jumps in weak shocks formed due to the steepening of
nonlinear acoustic waves (Section II.1). We then briefly discuss the variable
decomposition and non-dimensionalization in Section II.2, followed by the
derivation of second order governing equations for nonlinear acoustics in
Section II.3.

####  II.1 Fully compressible 1D Navier-Stokes and entropy scaling in weak
shocks

One-dimensional governing equations of continuum gas dynamics (compressible
Navier-Stokes) for an ideal gas are given by,

$\displaystyle\frac{\partial\rho^{*}}{\partial
t^{*}}+\frac{\partial(\rho^{*}u^{*})}{\partial x^{*}}=0,$ (2.1)
$\displaystyle\frac{\partial}{\partial
t^{*}}\left(\rho^{*}u^{*}\right)+\frac{\partial}{\partial
x^{*}}\left(\rho^{*}u^{*2}\right)=-\frac{\partial p^{*}}{\partial
x^{*}}+\frac{\partial}{\partial
x}\left(\left(\frac{4}{3}\mu^{*}+\mu^{*}_{B}\right)\frac{\partial
u^{*}}{\partial x^{*}}\right),$ (2.2)
$\displaystyle\rho^{*}T^{*}\left(\frac{\partial s^{*}}{\partial
t^{*}}+u^{*}\frac{\partial s^{*}}{\partial
x^{*}}\right)=\frac{\partial}{\partial
x^{*}}\left(\frac{\mu^{*}C^{*}_{p}}{Pr}\frac{\partial T^{*}}{\partial
x^{*}}\right)+\left(\frac{4}{3}\mu^{*}+\mu^{*}_{B}\right)\left(\frac{\partial
u^{*}}{\partial x^{*}}\right)^{2},$ (2.3)

which are closed by the ideal gas equation of state,

$\displaystyle p^{*}=\rho^{*}R^{*}T^{*},$ (2.4)

where $p^{*},u^{*},\rho^{*},T^{*},s^{*}$, respectively, denote total pressure,
velocity, density, temperature, and entropy of the fluid, $x^{*}$ and $t^{*}$
denote space and time, and $\mu^{*}$ denotes dynamic viscosity. Results in
this chapter and subsequent chapters have been obtained by DNS of Eqs.
(2.1)-(2.3), to resolve all the length scales of planar nonlinear acoustic
waves. For our simulations (see Section III), we choose the gas specific
constants for air at standard temperature and pressure (STP),

$R^{*}=287.105~{}\frac{\mathrm{m}^{2}}{\mathrm{s}^{2}\cdot
K},\quad\quad\mu^{*}_{B}=0,\quad\quad Pr=0.72.$ (2.5)

$(a)$$(b)$

Figure 2.1: Weak shock wave structure $(a)$ pressure $p^{*}$ and $(b)$ entropy
$s^{*}$ propagating with a speed $a^{*}_{s}>a^{*}_{0}$ obtained from DNS (see
Section III). $\Delta s^{*}/R^{*}$ and $s^{*}_{\mathrm{max}}/R^{*}$ are the
entropy jump and maximum entropy respectively. With increasing viscosity, the
peak in entropy remains constant. The DNS data has been obtained for base
state viscosity values given in Table 3.1. Figure 2.2: Entropy jump $\Delta
s^{*}=s^{*}_{2}-s^{*}_{1}$ and maximum entropy generated
$s^{*}_{\mathrm{max}}$ versus pressure jump $\Delta p^{*}$ across a planar
shock wave. In the labeled region ($\Delta p^{*}/\gamma p^{*}_{1}<1$, referred
as ‘weak shocks’ hereafter), the entropy jump $\Delta s^{*}$ scales as
$\mathcal{O}\left(\Delta p^{*3}\right)$, whereas the maximum entropy generated
$s^{*}_{\mathrm{max}}$ scales as $\mathcal{O}\left(\Delta p^{*2}\right)$,
approximately. Markers denote DNS data (see Section III), ( , )
$\mu^{*}=7.5\times 10^{-3}~{}$kg$\cdot$m${}^{-1}\cdot$s-1; ( , ),
$\mu^{*}=7.5\times 10^{-4}~{}$kg$\cdot$m${}^{-1}\cdot$s-1; ( , ),
$\mu^{*}=7.5\times 10^{-5}~{}$kg$\cdot$m${}^{-1}\cdot$s-1 for varying values
of $\Delta p^{*}$. Solid lines correspond to Eqs. (2.6) and (2.8).

Planar nonlinear acoustic waves steepen and form weak shocks. For weak shocks,
the smallest length scale (shock-thickness) is also significantly larger than
the molecular length scales. Hence, in this work, we neglect the molecular
vibrational effects in the single component ideal gas ($\mu^{*}_{B}=0$) ,
typically modeled via bulk viscosity effects Cramer_PoF_2012 . Across a freely
propagating planar weak shock (figure 2.1), the entropy jump ($\Delta
s^{*}=s^{*}_{2}-s^{*}_{1}$) is given by the classical gas-dynamic relation
LiepmannRoshko ,

$\displaystyle\frac{\Delta s^{*}}{R^{*}}$
$\displaystyle=\frac{1}{\gamma-1}\ln\left(1+\frac{2\gamma}{\gamma+1}\left(M^{2}-1\right)\right)-\frac{\gamma}{\gamma-1}\ln\left(\frac{\gamma+1}{\gamma-1+2/M^{2}}\right),$
(2.6)

where $M$ is the Mach number, given by,

$\frac{\Delta p^{*}}{\gamma p^{*}_{1}}=\frac{p^{*}_{2}-p^{*}_{1}}{\gamma
p^{*}_{1}}=\frac{2}{\gamma+1}\left(M^{2}-1\right),$ (2.7)

and $\Delta p^{*}=p^{*}_{2}-p^{*}_{1}$ is the pressure jump with $p^{*}_{1}$
and $p^{*}_{2}$ being the pre-shock and post-shock pressures, respectively.
Near the inflection point of the fluid velocity profile, the entropy reaches a
local maximum $(s^{*}=s^{*}_{\mathrm{max}})$. According to Morduchow and Libby
Morris_JASc_1949 , maximum entropy $s^{*}_{\mathrm{max}}$ assuming
$\mu^{*}_{B}=0$ and $Pr=3/4$, can be obtained as,

$\frac{s^{*}_{\mathrm{max}}}{R^{*}}=\frac{1}{\gamma-1}\ln\left(1+\frac{\gamma-1}{2}M^{2}\left(1-\xi\right)\xi^{\frac{\gamma-1}{2}}\right),$
(2.8)

where,

$\xi=\frac{\gamma-1}{\gamma+1}+\frac{2}{\gamma+1}\frac{1}{M^{2}}.$ (2.9)

For weak shock waves, ($\Delta p^{*}/\gamma p^{*}_{1}<1$), the entropy jump
$\Delta s^{*}$ and maximum entropy $s^{*}_{\mathrm{max}}$ scale with pressure
jumps as (cf. figure 2.2),

${\Delta s^{*}}=\mathcal{O}\left(\Delta p^{*3}\right),\quad
s^{*}_{\mathrm{max}}=\mathcal{O}\left({\Delta p}^{*2}\right),$ (2.10)

independent of $\mu^{*}$ (cf. Eqs. (2.6) and (2.8)). The overall entropy jump
$\Delta s^{*}$ is due to irreversible thermoviscous losses occurring within
the shocks. However, the overshoot in entropy ($s^{*}_{\mathrm{max}}>\Delta
s^{*}$) is due to both reversible and irreversible processes, and is not in
violation of the second law of thermodynamics Morris_JASc_1949 . Moreover, in
the range of pressure jumps considered in the DNS (see Section III), the
maximum Mach number of the shock is around $M\approx 1.4$, which is well
within the limits of validity of the continuum approach
Bird_MolecularGasDynamicsBook_1994 . Hence, it is physically justified to draw
conclusions regarding the smallest length scales through the governing
equations based on continuum approach and assuming thermodynamic equilibrium.

####  II.2 Perturbation variables and non-dimensionalization

In this section, we utilize the previous consideration on the second order
scaling of the maximum entropy $s^{*}_{\mathrm{max}}$ inside a weak shock wave
to derive second order nonlinear acoustics equations. To this end, we
decompose the variables in base state and perturbation fields and derive
equations containing only linear and quadratic terms in perturbation fields.
Denoting the base state with the superscript $()_{0}$ and the perturbation
fields with the superscript $({})^{\prime}$, we obtain,

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

where 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}.$ (2.12)

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

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

where $L^{*}$ is the length of the one-dimensional periodic domain. As also
typically done in classical studies of homogeneous isotropic turbulence
Ishihara_ARFM_2009 ; Batchelor1953theory ; pope2000turbulent ;
monin1971statistical ; Hosokawa_PRE_2008 ; Burattini_PRE_2006 , periodic
boundary conditions represent a common (yet not ideal) way to approximate
infinite domains; as such, a spurious interaction between the flow physics
that one wishes to isolate and the periodic box size may occur. For the TW and
SW test cases analyzed herein, $L^{*}$ corresponds to the initial (and hence
largest) reference length scale of the acoustic perturbation; in the AWT case,
the value of $L^{*}$ should be chosen as much larger than the integral length
scale $\ell$ or Taylor microscale $\lambda$ (see Section IV), which truly
define the state of turbulence.

Due to thermodynamic nonlinearities, wave propagation velocity increases
across a high-amplitude compression front, resulting in wave-steepening
Whitham2011 and hence generation of small length scales associated with
increasing temperature and velocity gradients responsible for thermoviscous
dissipation. The increase in thermoviscous dissipation results in positive
entropy perturbations peaking within the shock structure. For pressure jumps
$\Delta p^{*}/\gamma p^{*}_{1}<1$, the maximum entropy scales approximately as
$\mathcal{O}\left(\Delta p^{*2}\right)$ (cf. figure 2.2). Moreover, as we
discuss in a later section (see Section II), the second order nonlinear
acoustic equations impose a strict limit of $|p^{\prime}|<1/\gamma$ ($\simeq
0.714$ for $\gamma=0.72$) for base state normalized (Eq. (2.13)) (not pre-
shock state normalized (Eq. (2.7))) perturbations. Hence, in our simulations
(see Section III), we consider a suitable range of
$10^{-3}<p^{\prime}<10^{-1}$, which satisfies the aforementioned constraints.
Thus, the second order scaling of entropy holds in our simulations.

Below, we utilize this entropy scaling to derive the correct second order
nonlinear acoustics equations governing the spatio-temporal evolution of
dimensionless perturbation variables $p^{\prime}$ and $u^{\prime}$, as defined
in Eq. (2.13).

####  II.3 Second order nonlinear acoustics equations

For a thermally perfect gas, the differential in dimensionless density $\rho$
can be related to differentials in dimensionless pressure $p$ and
dimensionless entropy $s$ as,

$\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-\frac{\rho(\gamma-1)}{\gamma}ds.$
(2.14)

Nondimensionalizing the continuity Eq. (2.1) and substituting Eq. (2.14), we
obtain,

$\displaystyle\frac{\partial\rho}{\partial t}+\frac{\partial\rho u}{\partial
x}=0,$ $\displaystyle\frac{\partial p}{\partial t}+u\frac{\partial p}{\partial
x}+\gamma p\frac{\partial u}{\partial x}=(\gamma-1)p\left(\frac{\partial
s}{\partial t}+u\frac{\partial s}{\partial x}\right).$ (2.15)

Substituting the dimensionless forms of Eqs. (2.3) and (2.4) and utilizing the
decomposition given in Eqs. (2.13), we obtain the following truncated equation
for pressure perturbation $p^{\prime}$,

$\displaystyle\frac{\partial p^{\prime}}{\partial t}+\frac{\partial
u^{\prime}}{\partial x}+\gamma p^{\prime}\frac{\partial u^{\prime}}{\partial
x}$ $\displaystyle+u^{\prime}\frac{\partial p^{\prime}}{\partial
x}=\nu_{0}\left(\frac{\gamma-1}{Pr}\right)\frac{\partial^{2}p^{\prime}}{\partial
x^{2}}$ $\displaystyle+\mathcal{O}\left(p^{\prime}s^{\prime},s^{\prime
2},p^{\prime 3},\left(\frac{\partial u^{\prime}}{\partial
x}\right)^{2}\right).$ (2.16)

Similarly, the truncated equation for velocity perturbation $u^{\prime}$ is
obtained as,

$\displaystyle\frac{\partial u^{\prime}}{\partial t}+\frac{\partial
p^{\prime}}{\partial x}+\frac{\partial}{\partial x}\left(\frac{u^{\prime
2}}{2}-\frac{p^{\prime 2}}{2}\right)$
$\displaystyle=\frac{4}{3}\nu_{0}\frac{\partial^{2}u^{\prime}}{\partial
x^{2}}$ $\displaystyle+\mathcal{O}\left(\rho^{\prime 2}p^{\prime},\rho^{\prime
3}p^{\prime}\right).$ (2.17)

In Eqs. (2.16) and (2.17), $\nu_{0}$ is the dimensionless kinematic viscosity
given by,

$\nu_{0}=\frac{\mu^{*}_{0}}{\rho^{*}_{0}a^{*}_{0}L^{*}},$ (2.18)

and quantifies viscous dissipation of waves relative to propagation. Equations
(2.16) and (2.17) constitute the nonlinear acoustics equations truncated up to
second order, governing spatio-temporal evolution of finite amplitude acoustic
perturbations $p^{\prime}$ and $u^{\prime}$. The entropy scaling
($s^{*}_{\mathrm{max}}=\mathcal{O}\left(\Delta p^{*2}\right)$) discussed
previously results in the dissipation term on the right hand side of Eq.
(2.16). The left hand side of Eqs. (2.16) and (2.17) contains terms denoting
linear and nonlinear isentropic acoustic wave propagation. The Detailed
derivation of Eqs. (2.16) and (2.17) is given in Appendix 10, where we also
show that the nonlinear terms on the left hand side (LHS) of Eqs. (2.16) and
(2.17) are independent of the thermal equation of state. The functional form
of the second order perturbation energy norm ($E^{(2)}$, Eq (3.10)) – being
exclusively dictated by such terms (see Section II) – is independent of the
thermal equation of state of the gas. The results shown in this work focus on
ideal-gas simulations merely for the sake of simplicity, with no loss of
generality pertaining to inviscid nonlinear (up to second order) spectral
energy transfer dynamics.

We note that, Eq. (2.16) consists of the velocity derivative term ($\gamma
p^{\prime}{\partial u^{\prime}}/{\partial x}$), and is different from those
obtained by Naugol’nykh and Rybak Naugol1975spectrum , which in dimensionless
form read,

$\displaystyle\frac{\partial p^{\prime}}{\partial
t}-(\gamma-1)p^{\prime}\frac{\partial p^{\prime}}{\partial t}+\frac{\partial
u^{\prime}}{\partial x}+p^{\prime}\frac{\partial u^{\prime}}{\partial x}=0,$
(2.19) $\displaystyle\frac{\partial u^{\prime}}{\partial t}+\frac{\partial
p^{\prime}}{\partial x}+\frac{\partial}{\partial x}\left(\frac{u^{\prime
2}}{2}-\frac{p^{\prime 2}}{2}\right)=0.$ (2.20)

We adopt Eqs. (2.16) and (2.17) throughout the study since they represent the
truncated governing equations exactly. Unlike Naugol’nykh and Rybak
Naugol1975spectrum , we do not approximate the density $\rho$ using the Taylor
series and only use the total differential form given in Eq. (2.14).

Additionally, we note that Eqs. (2.16) and (2.17) can be combined into
Westervelt’s equation Hamilton_NLA_1998 only if the Lagrangian defined as,

$\mathcal{L}=\frac{u^{\prime 2}}{2}-\frac{p^{\prime 2}}{2},$ (2.21)

is zero, which holds only for linear pure traveling waves. The derivation of
Burgers’ equation in nonlinear acoustics follows from the Westervelt’s
equation Hamilton_NLA_1998 . Hence, it is inadequate in modeling general
nonlinear acoustics phenomena involving mixed phasing of nonlinear waves which
occurs in the Standing Wave (SW) and Acoustic Wave Turbulence (AWT) cases
analysed here.

##  Chapter 3 Energy cascade and decay of nonlinear acoustic waves

The contents of this chapter have been published as a regular article in
_Physical Review E_ GuptaScalo_PRE_2018 (see Sections III - VII) and have been
reported here in abridged form with minor modifications.

###  I Introduction

In this chapter, we analyze the energy dynamics of decaying nonlinear acoustic
waves utilizing three canonical configurations, traveling waves (TW), standing
waves (SW), and randomly initialized Acoustic Wave Turbulence (AWT) utilizing
the second order equations derived in the previous chapter and DNS of the
fully compressible Navier-Stokes equations. The spectral energy and decay
dynamics of 1D Burgers’ turbulence have been studied extensively by Kida
Kida_1979_JFM , Gurbatov et al. Gurbatov_JFM_1997 ; Gurbatov1991 , Woyczynski
Woyczynski_Gottingen_2006 , Fournier and Frisch Fournier1983Burgers , and
Burgers Burgers1974nonlinear . However, the equations of second order
nonlinear acoustics can be reduced to Burgers’ equations only with the
restrictive assumption of pure traveling waves (TW). Consequently, we utilize
the continuum gas dynamics governing equations to elucidate the spectral
energy cascade and decay dynamics of nonlinear acoustics. The decay of
nonlinear perturbations is governed by viscosity and thermal conductivity
(thermoviscous diffusion). The decay exhibits a power law in time due to the
gradual increase of the length scale over which the thermoviscous diffusion
acts. Such decay dynamics occur due to the separation of energy containing and
diffusive length scales and resemble those of decaying homogeneous isotropic
turbulence (HIT), which also have been studied extensively in the literature
both theoretically and experimentally batchelor1948decay ; batchelor1949nature
; hinze1975turbulence ; Ishihara_ARFM_2009 .

We develop theoretical arguments based on the second order nonlinear acoustics
equations derived in the previous chapter, and compare our results with DNS
results of fully compressible Navier Stokes Eqs. (2.1)-(2.3). We analyze
various length scales associated with the decaying nonlinear acoustic waves,
as well as the spectra of waveforms.

###  II Second order perturbation energy

In this section we derive a new perturbation energy function for nonlinear
acoustic waves utilizing Eqs. (2.16) and (2.17). To this end, we derive the
perturbation energy conservation relation (energy corollary) for high
amplitude acoustic perturbations. We show that the spatial average of the
perturbation energy function satisfies the definition of the Lyapunov function
for high amplitude acoustic perturbations and evolves monotonically in time
(cf. figure 3.2). Utilizing the energy corollary, we derive spectral energy
transport relations in further sections.

Multiplying Eqs. (2.16) and (2.17) with $p^{\prime}$ and $u^{\prime}$,
respectively, and adding, we obtain,

$\displaystyle\frac{\partial}{\partial t}\left(\frac{p^{\prime
2}}{2}+\frac{u^{\prime 2}}{2}\right)+\frac{\partial}{\partial
x}\left(u^{\prime}p^{\prime}+\frac{u^{\prime 3}}{3}\right)+\gamma p^{\prime
2}\frac{\partial u^{\prime}}{\partial
x}=\nu_{0}\left(\frac{\gamma-1}{Pr}\right)p^{\prime}\frac{\partial^{2}p^{\prime}}{\partial
x^{2}}+\frac{4}{3}\nu_{0}u^{\prime}\frac{\partial^{2}u^{\prime}}{\partial
x^{2}}.\quad\quad$ (3.1)

Spatial averaging of Eq. (3.1) over a periodic domain $[0,L]$ yields,

$\displaystyle\frac{d\left\langle
E^{(1)}\right\rangle}{dt}=-\left\langle\gamma p^{\prime 2}\frac{\partial
u^{\prime}}{\partial x}\right\rangle$
$\displaystyle-\nu_{0}\left(\frac{\gamma-1}{Pr}\right)\left\langle\left(\frac{\partial
p^{\prime}}{\partial
x}\right)^{2}\right\rangle-\frac{4}{3}\nu_{0}\left\langle\left(\frac{\partial
u^{\prime}}{\partial x}\right)^{2}\right\rangle,$ (3.2)

where $\left\langle.\right\rangle$ is the spatial averaging operator,

$\left\langle.\right\rangle=\frac{1}{L}\int^{L}_{0}\left(.\right)dx,$ (3.3)

and

$E^{(1)}=\frac{u^{\prime 2}}{2}+\frac{p^{\prime 2}}{2},$ (3.4)

is the first order isentropic acoustic energy. Equation (3.2) suggests that,
in a lossless medium ($\nu_{0}\rightarrow 0$), $\left\langle
E^{(1)}\right\rangle$ would exhibit spurious non-monotonic behavior in time
due to the first term on right hand side. Such non-monotonic behaviour is
confirmed by the DNS results shown in figure 3.2. Consequently, the linear
acoustic energy norm $E^{(1)}$ does not quantify the perturbation energy
correctly for high amplitude perturbations since the spatial average
$\left\langle E^{(1)}\right\rangle$ supports spurious growth and decay in the
absence of physical sources of energy.

Figure 3.1: Comparison of perturbation energy function for nonlinear acoustic
waves $E^{(2)}$ with the linear acoustic energy $E^{(1)}$ in the case of
$p^{\prime}=u^{\prime}$ (assumed for illustrative purpose). The correction
$f(p^{\prime})$ is independent of $u^{\prime}$.

The corrected perturbation energy function can be obtained upon recursively
evaluating the velocity derivative term ($\gamma p^{\prime 2}\partial
u^{\prime}/\partial x$ in Eq. (3.1)) utilizing Eq. (2.16) as,

$\displaystyle\gamma p^{\prime 2}\frac{\partial u^{\prime}}{\partial
x}=-\frac{\partial}{\partial t}\left(\frac{\gamma p^{\prime
3}}{3}\right)-\frac{\partial}{\partial x}\left(\frac{\gamma
u^{\prime}p^{\prime 3}}{3}\right)-$
$\displaystyle\gamma\left(\gamma-\frac{1}{3}\right)p^{\prime 3}\frac{\partial
u^{\prime}}{\partial x}-\frac{\nu_{0}(\gamma-1)}{Pr}\gamma p^{\prime
2}\frac{\partial p^{\prime}}{\partial x}.$ (3.5)

Furthermore, the third term in above Eq. (3.5) on the right can be further
evaluated as,

$\displaystyle\gamma\left(\gamma-\frac{1}{3}\right)p^{\prime 3}\frac{\partial
u^{\prime}}{\partial x}=-\frac{\partial}{\partial
t}\left(\frac{\gamma}{4}\left(\gamma-\frac{1}{3}\right)p^{\prime 4}\right)$
$\displaystyle-\frac{\partial}{\partial
x}\left(\frac{\gamma}{4}\left(\gamma-\frac{1}{3}\right)u^{\prime}p^{\prime
4}\right)-\gamma\left(\gamma-\frac{1}{3}\right)\left(\gamma-\frac{1}{4}\right)p^{\prime
4}\frac{\partial u^{\prime}}{\partial
x}-\frac{\nu_{0}(\gamma-1)}{Pr}\gamma\left(\gamma-\frac{1}{3}\right)p^{\prime
3}\frac{\partial^{2}p^{\prime}}{\partial x^{2}},$ (3.6)

and so on. Continued substitution according to Eqs. (3.5) and (3.6) yields the
closure of the system and the following energy corollary,

$\frac{\partial E^{(2)}}{\partial t}+\frac{\partial I}{\partial
x}=\nu_{0}\left(\frac{\gamma-1}{Pr}\right)h(p^{\prime})\frac{\partial^{2}p^{\prime}}{\partial
x^{2}}+\frac{4}{3}\nu_{0}u^{\prime}\frac{\partial^{2}u^{\prime}}{\partial
x^{2}},$ (3.7)

where,

$I(p^{\prime},u^{\prime})=p^{\prime}u^{\prime}+\frac{u^{\prime
3}}{3}+u^{\prime}f(p^{\prime}),$ (3.8)

is the intensity (energy flux) of the field, $h(p^{\prime})$ is given by,

$h(p^{\prime})=p^{\prime}+\frac{\partial f(p^{\prime})}{\partial
p^{\prime}}=\frac{\partial E^{(2)}}{\partial p^{\prime}}.$ (3.9)

and $E^{(2)}$ is given by,

$E^{(2)}(p^{\prime},u^{\prime})=\frac{u^{\prime 2}}{2}+\frac{p^{\prime
2}}{2}+f(p^{\prime})=E^{(1)}+f(p^{\prime}),$ (3.10)

and defines the second order perturbation energy for high amplitude acoustic
perturbations. The energy corollary Eq. (3.7) is mathematically exact for the
governing Eqs. (2.16) and (2.17).

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

Figure 3.2: Spatial profile of finite amplitude waves (top) for TW($a$),
SW($b$), and AWT($c$). Evolution of the average perturbation energy
($\left\langle E^{(2)}\right\rangle$ (–); $\left\langle E^{(1)}\right\rangle$
($--$)) evaluated from the DNS data (bottom) scaled by the initial value
against scaled time $t/\tau$ (cf. Eq. (3.31)) for increasing values of
perturbation amplitude $A_{\mathrm{rms}}$ defined in Eq. (3.16) at $\nu_{0}=$
1.836$\times$10-7 (see Table 3.1). The curves are shifted vertically by 0.25
for illustrative purpose only. With increasing perturbation amplitude
$A_{\mathrm{rms}}$, the variation of linear acoustic energy norm < $E^{(1)}$ >
becomes increasingly non-monotonic. The vertical dashed line (bottom)
highlights the end of approximately inviscid spectral energy cascade regime.
In this regime, the energy is primarily redistributed in the spectral space
due to the nonlinear propagation ($\epsilon\simeq 0$).

The correction term $f(p^{\prime})$ in $E^{(2)}$ appears due the thermodynamic
nonlinearities and can be derived in the closed form as,

$f(p^{\prime})=\sum^{\infty}_{n=2}T_{n}=\sum^{\infty}_{n=2}(-1)^{n+1}\frac{\gamma
p^{\prime n+1}}{n+1}\prod^{n}_{i=3}\left(\gamma-\frac{1}{i}\right),$ (3.11)

where $T_{2}$ and $T_{3}$ can be identified in Eqs. (3.5) and (3.6),
respectively. Upon isolating the $n^{th}$ term of the above infinite series
as,

$T_{n}=(-1)^{n+1}\frac{\gamma p^{\prime
n+1}}{n+1}\underbrace{\left(\gamma-\frac{1}{3}\right)\left(\gamma-\frac{1}{4}\right)\cdots\left(\gamma-\frac{1}{n}\right)}_{n-2~{}\mathrm{terms}}.$
(3.12)

Multiplied fractions in the Eq. (3.12) above yield the $n^{th}$ term as,

$T_{n}=-\frac{2\gamma}{\left(\gamma-1\right)\left(2\gamma-1\right)}\left(\gamma
p^{\prime}\right)^{n+1}\binom{1/\gamma}{n+1}.$ (3.13)

Finally, the energy correction $f(p^{\prime})$ can be recast as,

$\displaystyle f(p^{\prime})=\sum^{\infty}_{n=2}T_{n}=$
$\displaystyle-\frac{2\gamma}{\left(\gamma-1\right)\left(2\gamma-1\right)}\Big{(}\left(1+\gamma
p^{\prime}\right)^{1/\gamma}-1-p^{\prime}+\frac{\left(\gamma-1\right)p^{\prime
2}}{2}\Big{)}.$ (3.14)

The correction function $f(p^{\prime})$ defined in Eq. (3.14) accounts for
second order isentropic nonlinearities and is not a function of entropy
perturbation. Hence, $E^{(2)}$ accounts for the effect of high amplitude
perturbations on perturbation energy isentropically. We note that this
separates $E^{(2)}$ fundamentally from generalized linear perturbation energy
norms, such as the ones derived by Chu Chu1965EnergyNorm for small amplitude
non-isentropic perturbations, and by Meyers myers1986exact for acoustic wave
propagation in a steady flow. Moreover, as discussed in the previous section
(and shown in Appendix 10), since the isentropic nonlinearities on the LHS of
Eqs. (2.16) and (2.17) are independent of the thermal equation of state, the
functional form of $E^{(2)}$ and $I$ are also independent of the equation of
state. However, the dissipation term on the right hand side of the energy
corollary Eq. (3.7) may change with the thermal equation of state.

The energy correction $f(p^{\prime})$ is infinite order in pressure
perturbation $p^{\prime}$ and converges only for perturbation magnitude
$|p^{\prime}|<1/\gamma$ thus naturally yielding the strict limit of validity
of second order acoustic equations in modeling wave propagation and wave
steepening. Figure 3.1 shows the newly derived second order perturbation
energy $E^{(2)}$ compared against the isentropic acoustic energy $E^{(1)}$.
Both $E^{(2)}$ and $E^{(1)}$ are non-negative in the range
$|p^{\prime}|<1/\gamma$ ($p^{\prime}=u^{\prime}$ is assumed for illustrative
purpose). Furthermore, $E^{(2)}$ is asymmetric in nature, with larger energy
in dilatations compared to compressions of same magnitude, as shown in figure
3.1. Such asymmetry signifies that the medium (compressible ideal gas in the
present study) relaxes towards the base state faster for finite dilatations
compared to compressions.

For compact supported or spatially periodic perturbations, the energy
conservation Eq. (3.7) shows that the spatially averaged energy $\left\langle
E^{(2)}\right\rangle$ decays monotonically in time (in the absence of energy
sources) accounting for the nonlinear interactions i.e.,

$\displaystyle\dot{\mathscr{V}}=\frac{d\left\langle E^{(2)}\right\rangle}{dt}$
$\displaystyle=-\nu_{0}\left(\frac{\gamma-1}{{Pr}}\right)\left\langle\frac{\partial^{2}E^{(2)}}{\partial
p^{\prime 2}}\left(\frac{\partial p^{\prime}}{\partial
x}\right)^{2}\right\rangle-\frac{4}{3}\nu_{0}\left\langle\left(\frac{\partial
u^{\prime}}{\partial x}\right)^{2}\right\rangle$
$\displaystyle=\left\langle\mathcal{D}\right\rangle=-\epsilon\leq 0,$ (3.15)

where $\mathcal{D}$ is the perturbation energy dissipation and $\epsilon$ is
the negative of its spatial average. The spatial average $\left\langle
E^{(2)}\right\rangle$ is non-negative ($E^{(2)}\geq 0$), and Eq. (3.15) and
figure 3.2 confirm that $\left\langle E^{(2)}\right\rangle$ evolves
monotonically in time in the absence of physical energy sources. Hence, the
spatial average of the perturbation energy function $\left\langle
E^{(2)}\right\rangle$ defines the Lyapunov function $\mathscr{V}$ of the
nonlinear acoustic system governed by the set of second order governing Eqs.
(2.16) and (2.17) exactly. The spatial average $\left\langle
E^{(2)}\right\rangle$ should be used for studying the stability of nonlinear
acoustic systems George2012 ; Strogatz2018NonlinearDynamics , which, however,
falls beyond the scope of this work.

Wave-front steepening entails cascade of perturbation energy into higher
wavenumbers thus broadening the energy spectrum. A fully broadened spectrum of
acoustic perturbations exhibits energy at very small length scales which
causes high thermoviscous energy dissipation. We analyse the separation of
length scales and energy decay caused by nonlinear wave steepening and
thermoviscous energy dissipation in the following sections. To this end, we
utilize the direct numerical integration of Navier-Stokes Eqs. (2.1)-(2.4)
resolving all the length scales (DNS) and the exact energy corollary Eq. (3.7)
for second order truncated Eqs. (2.16) and (2.17).

###  III High Fidelity Simulations with Adaptive Mesh Refinement

We perform shock-resolved numerical simulations of 1D Navier-Stokes (DNS) Eqs.
(2.1)-(2.4) with Adaptive Mesh Refinement (AMR). We use the perturbation
energy $E^{(2)}$ defined in Eq. (3.10) to define the characteristic
dimensionless perturbation amplitude $A_{\mathrm{rms}}$ as,

$A_{\mathrm{rms}}=\sqrt{\left\langle E^{(2)}\right\rangle},$ (3.16)

which is varied in the range $10^{-3}-10^{-1}$. The dimensionless kinematic
viscosity at base state $\nu_{0}$ is also varied from 1.836$\times$10-5 to
1.836$\times$10-7. The base state conditions in the numerical simulations
correspond to STP, i.e. $p^{*}_{0}=101325~{}\mathrm{Pa}$ and
$T^{*}_{0}=300~{}\mathrm{K}$.

Table 3.1: Simulation parameter space for TW, SW, and AWT cases listing base state dimensionless viscosity $\nu_{0}$ (cf. Eq. (2.18)), initial characteristic perturbation amplitude $A_{\mathrm{rms},0}$ (cf. Eq. (3.16)), and dimensional characteristic perturbation in velocity $u^{*}_{\mathrm{rms}}$ and pressure $p^{*}_{\mathrm{rms}}$ fields (Eq. (3.17)). $\nu_{0}$ |  1.836$\times$10-5 |  1.836$\times$10-6 |  1.836$\times$10-7
---|---|---|---
$A_{\mathrm{rms},0}$ | $10^{-3}$ | $10^{-2}$ | $10^{-1}$
$u^{*}_{\mathrm{rms}}~{}\left(\mathrm{m}/\mathrm{s}\right)$ | 0.347 | 3.472 | 34.725
$p^{*}_{\mathrm{rms}}~{}\left(\mathrm{kPa}\right)$ | 0.142 | 1.419 | 14.185

The goal of spanning $A_{\mathrm{rms}}$ and $\nu_{0}$ over three orders of
magnitude is to achieve widest possible range of energy cascade rate and
dissipation within computationally feasible times. Equation (3.16) yields the
definitions of the perturbation Reynolds number $\mathrm{Re}_{L}$ ,
characteristic perturbation velocity field $u^{*}_{\mathrm{rms}}$, and
pressure field $p^{*}_{\mathrm{rms}}$ as,

$\mathrm{Re}_{L}=\frac{A_{\mathrm{rms}}a^{*}_{0}L^{*}}{\nu_{0}^{*}},~{}\quad
u^{*}_{\mathrm{rms}}=a^{*}_{0}A_{\mathrm{rms}},~{}\quad
p^{*}_{\mathrm{rms}}=\rho^{*}_{0}{a_{0}^{*}}^{2}A_{\mathrm{rms}},$ (3.17)

where $\mathrm{Re}_{L}$ denotes ratio of diffusive to wave steepening time
scale over the length $L$. In the simulations, we keep $\mathrm{Re}_{L}\gg 1$,
which corresponds to very fast wave steepening rates compared to diffusion. In
further sections (see Section IV), we define the wave turbulence Reynolds
number $\mathrm{Re}_{\ell}$ based on the integral length scale $\ell$. Below,
we briefly discuss the numerical scheme utilized for shock-resolved
simulations and outline the initialization of the three configurations (TW,
SW, and AWT) for numerical simulations.

####  III.1 Numerical approach

Figure 3.3: Illustration of the binary tree implementation of Adaptive Mesh
Refinement (AMR) technique (top left). The mesh is refined based on the
resolution error in pressure field in each cell acting as a node of a binary
tree. The pressure field shown (middle) corresponds to the randomly
initialized AWT case with $A_{\mathrm{rms}}=10^{-1}$,
$\nu_{0}=1.836\times$10-6 (Table 3.1) at $t/\tau=0.04$. The inset shows the
resolved shock wave with $(+)$ denoting the cell interfaces. The mesh
refinement levels (bottom) show the depth $d$ of the binary tree.

We integrate the fully compressible 1D Navier-Stokes Eqs. (2.1)-(2.3) in time
utilizing the staggered spectral difference (SD) spatial discretization
approach kopriva1996conservative . In the SD approach, the domain is
discretized into cells. Within each cell, the orthogonal polynomial
reconstruction of variables allows numerical differentiation with spectral
accuracy. We refer the reader to the work by Kopriva and Kolias
kopriva1996conservative for further details.

To accurately resolve spectral energy dynamics at all length scales, i.e. for
resolved weak shock waves, we combine the SD approach with the adaptive mesh
refinement (AMR) approach as first introduced by Mavriplis
mavriplis1994adaptive for spectral methods. The SD-AMR approach eliminates the
computational need of very fine grid everywhere for resolving the propagating
shock waves. To this end, we expand the values of a generic variable $\phi$
local to the cell in the Legendre polynomial space as,

$\phi=\sum^{N}_{i=1}\hat{\phi}_{i}\psi_{i}(x)$ (3.18)

where $\psi_{i}(x)$ is the Legendre polynomial of $(i-1)^{\mathrm{th}}$
degree. The polynomial coefficients $\hat{\phi}_{i}$ are utilized for
estimating the local resolution error $\varepsilon$ defined as
mavriplis1994adaptive ,

$\varepsilon=\left(\frac{2\hat{\phi}^{2}_{N}}{2N+1}+\int^{\infty}_{N+1}\frac{2f_{\varepsilon}^{2}(n)}{2n+1}dn\right)^{1/2},\quad
f_{\varepsilon}(n)=ce^{-\sigma n},$ (3.19)

where $f_{\varepsilon}$ is the exponential fit through the coefficients of the
last four modes in the Legendre polynomial space. When the estimated
resolution error $\varepsilon$ exceeds a pre-defined tolerance, the cell
divides into two subcells, which are connected utilizing a binary tree (shown
in figure 3.3). The subcells merge together if the resolution error decreases
below a pre-defined limit.

####  III.2 Initial conditions

Table 3.2: Initial spectral compositions for traveling wave (TW), standing wave (SW), and acoustic wave turbulence (AWT). $\delta\left(\cdot\right)$ is the Dirac delta function. | TW | SW | AWT
---|---|---|---
$k_{0}$ | 1 | 1 | 1
$k_{E}$ | 1 | 1 | 100
$b_{0}(k)$ | 0 | 0 | $e^{-(|k|-k_{E})^{2}}$
$\widehat{E}_{k}$ | ${A}_{\mathrm{rms}}^{2}\delta(k_{0})$ | ${A}_{\mathrm{rms}}^{2}\delta(k_{0})$ | ${A}_{\mathrm{rms}}^{2}$

We utilize the Riemann invariants for compressible Euler equations to
initialize the propagating traveling and standing wave cases in the numerical
simulations. The Riemann invariants in terms of perturbation variables
assuming nonlinear isentropic changes are given by,

$\displaystyle
R_{-}=\frac{2}{\gamma-1}\left(\left(1+{\rho^{\prime}}\right)^{\frac{\gamma-1}{2}}-1\right)-u^{\prime},$
(3.20) $\displaystyle
R_{+}=\frac{2}{\gamma-1}\left(\left(1+{\rho^{\prime}}\right)^{\frac{\gamma-1}{2}}-1\right)+u^{\prime},$
(3.21)

where $R_{-}$ and $R_{+}$ are the left and right propagating invariants,
respectively, and $u^{\prime}$ and ${\rho^{\prime}}$ are normalized velocity
and density perturbations, as defined in Eq. (2.13). Initial conditions for TW
and SW cases correspond to $R_{-}=0$ and $R_{-}=R_{+}$ respectively.

To initialize the broadband noise case, we first choose $p^{\prime}$ and
$u^{\prime}$ pseudo-randomly from a uniform distribution for the whole set of
discretization points in $x$. Low-pass filtering of $p^{\prime}$ and
$u^{\prime}$ yields,

$\displaystyle\widetilde{\widehat{p}_{k}}(t$
$\displaystyle=0)=\widehat{p}_{k}b_{0}(k),\quad\widetilde{\widehat{u}_{k}}(t=0)=\widehat{u}_{k}b_{0}(k)$
$\displaystyle b_{0}(k)$ $\displaystyle=\begin{cases}1&k_{0}\leq|k|\leq
k_{E}\\\ e^{-(|k|-k_{E})^{2}}&|k|>k_{E}\end{cases}.$ (3.22)

where ${\widehat{p}_{k}}$ and ${\widehat{p}_{k}}$ are the Fourier coefficients
of pseudo-random fields $p^{\prime}$ and $u^{\prime}$, respectively.
$\widetilde{\widehat{p}_{k}}$ and $\widetilde{\widehat{u}_{k}}$ are the low-
pass filtered coefficients. The inverse Fourier transform of Eq. (3.22) yields
smooth initial conditions with the initial spectral energy
${\widehat{E}_{k}}$, as defined in Section IV (cf. Eq. (3.27)). For TW and SW,
only the single harmonic ($k=1$ in the current work) contains all of the
initial energy. However, for AWT, $\widehat{E}_{k}$ is governed by the
correlation function of velocity and pressure fields. In Table 3.2, we
summarize the initial spectral energy for all three cases based on Eq. (3.22).

###  IV Scales of acoustic energy cascade and dissipation

In this section, we derive the analytical expressions of spectral energy,
energy cascade flux, and spectral energy dissipation utilizing the exact
energy corollary Eq. (3.7) (see Section II). We then identify the integral
length scale $\ell$, the Taylor microscale $\lambda$, and the Kolmogorov
length scale $\eta$ for TW, SW, and AWT cases in a periodic domain utilizing
the DNS data (see figures 3.4 and 3.5 and Table 3.3). Temporal evolution laws
of these length scales yield energy decay laws, which are used for
dimensionless spectral scaling relations (see Section V).

####  IV.1 Spectral energy flux and dissipation rate for periodic
perturbations

The exact perturbation energy conservation equation is given by (cf. Eq.
(3.7)),

$\frac{\partial E^{(2)}}{\partial t}+\frac{\partial I}{\partial
x}=\mathcal{D}$ (3.23)

Integrating over the periodic domain, the above energy corollary can be
converted into the following statement of conservation of perturbation energy
in the spectral space,

$\frac{d}{dt}\sum_{|k^{\prime}|\leq
k}\widehat{E}_{k^{\prime}}+\widehat{\Pi}_{k}=\sum_{|k^{\prime}|\leq
k}\widehat{\mathcal{D}}_{k^{\prime}},$ (3.24)

$(a)$$(b)$

Figure 3.4: Schematic illustrating the comparison of various length scales
associated with spectral energy cascade in nonlinear acoustics in both spatial
$(a)$ and spectral $(b)$ space. $(a)$ shows the perturbation velocity
$u^{\prime}$ ($--$) and pressure $p^{\prime}$ (–) fields in AWT obtained from
the DNS data for $\nu_{0}=$1.836$\times$10-7 and $A_{\mathrm{rms}}=10^{-1}$.
$(b)$ shows the corresponding spectral energy $\widehat{E}_{k}$ in log-log
space. The integral length scale $\ell$ corresponds to the characteristic
distance between the shock waves traveling in the same direction. The
Kolmogorov length scale $\eta$ corresponds to the shock wave thickness. The
Taylor microscale $\lambda$ is the diffusive length scale and satisfies
$\ell\gg\lambda\gg\eta$. $L$ corresponds to the length of the domain.

where the first term corresponds to the temporal rate of change of cumulative
spectral energy density,

$\frac{d\widehat{E}_{k}}{dt}\approx\frac{d}{dt}\left(\frac{|\widehat{u}_{k}|^{2}}{2}+\frac{|\widehat{p}_{k}|^{2}}{2}\right)+\Re\left(\widehat{p}_{-k}\frac{d\widehat{g}_{k}}{dt}\right),$
(3.25)

and $\widehat{g}$ is the Fourier transform of $g(p^{\prime})$ given by,

$g(p^{\prime})=\frac{\gamma}{\gamma-1}\left(\left(1+\gamma
p^{\prime}\right)^{1/\gamma}-1-p^{\prime}\right).$ (3.26)

The spectral energy $\widehat{E}_{k}$ is given by,

$\widehat{E}_{k}=\frac{|\widehat{u}_{k}|^{2}}{2}+\frac{|\widehat{p}_{k}|^{2}}{2}+\Re\left(\widehat{p}_{-k}\left(\widehat{\frac{f(p^{\prime})}{p^{\prime}}}\right)_{k}\right).$
(3.27)

$(a)$$(b)$

Figure 3.5:  Schematic illustrating the comparison of various length scales associated with spectral energy cascade in nonlinear acoustics in the spectral space, and corresponding variation of spectral energy flux $\hat{\Pi}_{k}$ ($a$) and spectral dissipation $\mathcal{D}_{k}$ $(b)$. The integral length scale $\ell$ corresponds to the characteristic distance between the shock waves traveling in the same direction. The Kolmogorov length scale $\eta$ corresponds to the shock wave thickness. The Taylor microscale $\lambda$ is the diffusive length scale and satisfies $\ell\gg\lambda\gg\eta$. $L$ corresponds to the length of the domain. Table 3.3:  Summary of the three length scales $\ell$, $\lambda$, and $\eta$, respective definitions, and the range of spectrum characterized by them. The integral length scale characterizes the energy containing range $(k_{0},k_{E})$. The Taylor microscale is the characteristic of the energy transfer and dissipation range $(k_{E},k_{\delta})$. The Kolmogorov length scale corresponds to the highest wavenumber generated as a result of nonlinear acoustic energy cascade. Length scale | Integral | Taylor | Kolmogorov
---|---|---|---
length scale | Microscale | length scale
| $\ell$ | $\lambda$ | $\eta$
Definition | $\sqrt{\frac{\sum_{k}{\widehat{E}_{k}}/{k^{2}}}{\sum_{k}\widehat{E}_{k}}}$ | $\sqrt{\frac{2\delta\left\langle E^{(2)}\right\rangle}{\epsilon}}$ | $\frac{\delta}{\sqrt{\left\langle E^{(2)}\right\rangle}}$
Characteristic | $(k_{0},k_{E})$ | $(k_{E},k_{\delta})$ | $(k_{\delta},\infty)$
spectral range

It is noteworthy that the correction in spectral energy does not follow
directly from the nonlinear correction function $f(p^{\prime})$ derived in the
physical space. In Eq. (3.25), we have made the following approximation,

$\frac{d}{dt}\left(\Re\left(\widehat{p}_{-k}\left(\widehat{\frac{f(p^{\prime})}{p^{\prime}}}\right)_{k}\right)\right)\approx\Re\left(\widehat{p}_{-k}\frac{d\widehat{g}_{k}}{dt}\right)$
(3.28)

The second term $\widehat{\Pi}_{k}$ in Eq. (3.24) is the flux of spectral
energy density from wavenumbers $|k^{\prime}|\leq k$ to $|k^{\prime}|>k$ and
is given by, Finally, the spectral dissipation $\widehat{\mathcal{D}}_{k}$ is
given by,

$\displaystyle\widehat{\mathcal{D}}_{k}=\nu_{0}\frac{\gamma-1}{Pr}$

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

Figure 3.6: Spectro-temporal evolution of $\widehat{E}_{k}$ (top) and spectral
flux $\widehat{\Pi}_{k}$ (bottom) for a TW ($a$), SW ($b$) and AWT ($c$).
Spectral flux $\widehat{\Pi}_{k}$ for a traveling wave simply increases
towards high wavenumbers. For a standing wave, $\widehat{\Pi}_{k}$ oscillates
at low wavenumbers cyclically due to collisions of oppositely traveling shock
waves while high wavenumber behaviour resembles that of a traveling shock. For
AWT, the spectral broadening occurs for $k>k_{E}$ with small fluctuations in
time for $k<k_{E}$.

Detailed derivation of Eqs. (3.24)-(LABEL:eq:_SpectralDissipation) is given in
appendix 11. Figures 3.4 and 3.5 summarize the typical shape of the spectral
energy $\widehat{E}_{k}$, spectral energy flux $\widehat{\Pi}_{k}$, and the
spectral dissipation $\widehat{\mathcal{D}}_{k}$ along with the relative
positions of the three relevant length scales, the integral length scale,
$\ell$, the Taylor microscale, $\lambda$, and the Kolomogorov length scale
$\eta$ in the spectral space. The spectro-temporal evolution of any
configuration of nonlinear acoustic waves can be quantified utilizing these
length scales and the respective evolution in time which is discussed in
detail in the subsections below. Table 3.3 summarizes these length scales and
the characteristic spectral range. In further sections, we discuss all the
spectral quantities as functions of absolute value of wavenumbers and drop the
$|.|$ notation for convenience.

The spectral energy flux $\widehat{\Pi}_{k}$, defined in Eq.
(LABEL:eq:_SpectralFlux), is in terms of interactions of the Fourier
coefficients of the pressure $\widehat{p}_{k}$ and velocity $\widehat{u}_{k}$
perturbations. For compact support or periodic perturbations,
$\widehat{\Pi}_{k}$ approaches zero in the limit of very large wavenumbers
$k\rightarrow\infty$,

$\lim_{k\rightarrow\infty}\widehat{\Pi}_{k}=\left\langle\frac{\partial
I}{\partial x}\right\rangle=0.$ (3.29)

The last two terms in Eq. (LABEL:eq:_SpectralFlux) result in
$\widehat{\Pi}_{k}\rightarrow 0$ for large $k$ for general acoustic phasing.
Hence, they are most relevant in SW and AWT cases. In a pure traveling wave
(TW), $u^{\prime}=p^{\prime}$ at first order due to which the last two terms
in Eq. (LABEL:eq:_SpectralFlux) become negligible. Furthermore, the sequence
of $\widehat{\Pi}_{k}$ also converges monotonically i.e.,

$\lim_{k\rightarrow\infty}\left(\widehat{\Pi}_{k-1}-\widehat{\Pi}_{k}\right)\rightarrow
0^{+},$ (3.30)

as shown in Figs. 3.5$(a)$ and 3.6. The flattening of the spectral energy flux
$\widehat{\Pi}_{k}$ (Eq. (3.30)) begins at a specific wavenumber $k_{\delta}$
associated to the Kolmogorov length scale, $\eta$, as shown in figures 3.4 and
3.5. The spectral energy $\widehat{E}_{k}$ deviates off the $k^{-2}$ decay
near the wavenumber $k_{\delta}$. Figure 3.6 shows the spectro-temporal
evolution of the spectral energy $\widehat{E}_{k}$ and the flux
$\widehat{\Pi}_{k}$ for TW, SW, and AWT prior to formation of shock waves.

For TW, $\widehat{\Pi}_{k}$ increases in time due to spectral broadening. In
SW, $\widehat{\Pi}_{k}$, while increasing, also oscillates at low wavenumbers
due to the periodic collisions of oppositely propagating shocks. A combination
of these processes takes place in a randomly initialized smooth finite
amplitude perturbation, which at later times develops into AWT. At later
times, nonlinear waves in all three configurations fully develop in to shock
waves. Up to the shock formation, the spectral dynamics of all configurations
simply involve increase of the spectral flux $\widehat{\Pi}_{k}$. The
dimensionless shock formation time $\tau$ can be estimated as,

$\tau=\frac{2}{(\gamma-1)A_{\mathrm{rms},0}}.$ (3.31)

Upon shock formation, the dynamic evolution of TW and SW remains
phenomenologically identical. The isolated shocks propagate and the total
perturbation energy of the system decays due to thermoviscous dissipation
localized around the shock wave. However, for AWT, along with collisions of
oppositely propagating shocks, those propagating in the same direction
coalesce due to differential propagating speeds. As we discuss below, this
modifies the energy decay and spectral energy dynamics in AWT significantly
compared to TW and SW.

In the sub-sections below, we elucidate the energy dynamics before and after
shock formation for TW, SW, and AWT. To this end, we define and discuss the
relevant length scales as mentioned above, namely: the Taylor microscale
$\lambda$, the integral length scale $\ell$, and the Kolmogorov length scale
$\eta$. Particular focus is given to the AWT case due to modified dynamics
caused by shock coalescence.

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

Figure 3.7: Temporal evolution of scaled total energy <$E^{(2)}$>
${A^{-2}_{\mathrm{rms},0}}$ (–) (top), dissipation rate
$\epsilon{A^{-3}_{\mathrm{rms},0}}$ ($--$) (mid) and normalized Taylor
microscale $\lambda/\sqrt{\delta\tau}$ (bottom) for TW ($a$), SW ($b$) and AWT
($c$) against the scaled time $t/\tau$ for varying perturbation Reynolds
number $\mathrm{Re}_{L}$. The time $t_{0}$ signifies fully broadened spectrum
of the perturbation field.

####  IV.2 Taylor microscale

In hydrodynamic turbulence, the Taylor microscale $\lambda$ separates the
inviscid length scales from the viscous length scales tennekes1972first ;
pope2000turbulent . Due to the spectral energy cascade in planar nonlinear
acoustics, we note that the spectral energy varies as $\widehat{E}_{k}\sim
k^{-2}$ due to the formation of shocks and the spectral dissipation due to
thermoviscous diffusion varies as $\widehat{\mathcal{D}}_{k}\sim
k^{2}\widehat{E}_{k}$. Consequently, the dissipation acts over most of the
length scales with $k>k_{E}$ (figure 3.5$(b)$), unlike hydrodynamic turbulence
where the viscous dissipation dominates only the smaller length scales
tennekes1972first ; pope2000turbulent . As shown in figure 3.5$(a)$, length
scales in the range $(k_{E},k_{\delta})$ exhibit both dissipation
$\widehat{\mathcal{D}}_{k}$ and energy transfer $\widehat{\Pi}_{k}$. For
$k>k_{\delta}$, $\widehat{\Pi}_{k}$ begins to converge monotonically to 0 and
the interval $(k_{\delta},1/\eta)$ primarily exhibits dissipation
$\widehat{\mathcal{D}}_{k}$ only. The Taylor microscale $\lambda$ quantifies
the length scale associated to the whole dissipation range.

Utilizing the definition of the total perturbation energy $\left\langle
E^{(2)}\right\rangle$ and the dissipation rate $\epsilon$ (cf. Eq. (3.15)),
the microscale $\lambda$ can be defined as,

$\lambda(t)=\sqrt{\frac{2\delta\left\langle E^{(2)}\right\rangle}{\epsilon}},$
(3.32)

where $\delta$ is the thermoviscous diffusivity, given by,

$\delta=\nu_{0}\left(\frac{4}{3}+\frac{\gamma-1}{\mathrm{Pr}}\right).$ (3.33)

Equation (3.32) indicates that the Taylor microscale can be identified as the
geometrical centroid of full energy spectrum, i.e.

$\lambda\sim\sqrt{\frac{\sum_{k}\widehat{E}_{k}}{\sum_{k}k^{2}\widehat{E}_{k}}}.$
(3.34)

As the smaller length scales (higher harmonics) are generated, the dissipation
rate $\epsilon$ tends to increase reaching a maximum in time. The increase of
dissipation rate $\epsilon$ implies decrease of the length scale $\lambda$ in
time. Minima of $\lambda$ indicates the fully-broadened spectrum of energy
limited by the thermoviscous diffusivity at very large wavenumbers. Further
spatio-temporal evolution of the system is dominated by dissipation thus
indicating the purely diffusive nature of the Taylor microscale, i.e.,

$\lambda\rightarrow C\sqrt{\delta t}.$ (3.35)

The temporal evolution of $\lambda$ is qualitatively similar for TW, SW, and
AWT, the constant $C$ in Eq. (3.35) differs for TW and SW compared with AWT
due to the different spatial structure of perturbations. The time $t_{0}$ at
which $\lambda$ reaches minimum signifies fully developed nonlinear acoustic
waves. In case of AWT, it signifies fully developed acoustic wave turbulence.

Figure 3.7 shows the decay of scaled total perturbation energy $\left\langle
E^{(2)}\right\rangle A_{\mathrm{rms},0}^{-2}$ and total dissipation rate
$\epsilon A_{\mathrm{rms},0}^{-3}$ for the TW, SW, and AWT. We note that the
total energy decays as a power law $t^{-2}$ for both TW and SW, whereas, for
AWT, the initial decay law is $t^{-2/3}$. Asymptotic evolution (at large $t$)
of the Taylor microscale follows from the decay laws as $\lambda=\sqrt{\delta
t}$ and $\lambda=\sqrt{3\delta t}$, respectively. Since energy decay law of a
single harmonic traveling and standing waves is rather trivial, we focus
primarily on the AWT case for further discussion.

####  IV.3 Integral length scale

$(a)$$(b)$

Figure 3.8: Evolution of the integral length scale $\ell$ $(a)$ and the
Reynolds number $\mathrm{Re}_{\ell}$ ($b$) defined in Eqs. (3.36) and (3.45),
respectively, for all the cases of AWT considered. For small thermoviscous
diffusivity, $\ell$ increases approximately as $t^{1/3}$ before saturating to
the dimensionless domain length $L=1$ and $\mathrm{Re}_{\ell}$ remains
approximately constant.

We identify the integral length scale $\ell$ as the characteristic length
scale of the energy containing scales. In general, random smooth broadband
noise (AWT) develops into an ensemble of shocks, propagating left and right in
a one-dimensional system. For an ensemble of shock waves distributed spatially
along a line, $\ell$ corresponds to the characteristic distance between
consecutive shock waves traveling in the same direction, as shown
schematically in figures 3.4 and 3.5. Formally, we define $\ell$ as,

$\ell=\sqrt{\frac{\sum_{k}\frac{\widehat{E}_{k}}{k^{2}}}{\sum_{k}\widehat{E}_{k}}},$
(3.36)

which is identical to the integral length scale defined in Burgers’ turbulence
Gurbatov_JFM_1997 . Definition in Eq. (3.36) yields the centroid wavenumber of
the initial energy spectrum (unlike the Taylor microscale, which corresponds
to the full energy spectrum) and hence is characteristic of the large length
scales of fully developed AWT. To elucidate the evolution of the total
perturbation energy $\left\langle E^{(2)}\right\rangle$ utilizing the integral
length scale, we assume the following model spectral energy density
$\widehat{E}_{k}$,

$\widehat{E}_{k}=\begin{cases}C_{1}k^{n}&k_{0}\leq k\leq k_{E}\\\
C_{2}k^{-2}&k_{\delta}>k>k_{E}\end{cases},$ (3.37)

where $k^{n}$ corresponds to the shape of initialized energy spectrum in the
range $(k_{0},k_{E})$ (figures 3.4 and 3.5). In this work, we only utilize the
white noise initialized AWT cases which correspond to $n=0$ (see Table 3.2).
Moreover,

$C_{2}=C_{1}k_{E}^{n+2}.$ (3.38)

The wave numbers $k_{E}$ and $k_{\delta}$ vary in time due to decaying energy.
By definition, the mean of perturbations is zero. Hence, the smallest
wavenumber containing energy $k_{0}$ (cf. figures 3.4 and 3.5) is the
reciprocal of the domain length $L$, i.e.,

$k_{0}=1/L.$ (3.39)

We note that the above model spectral energy $\widehat{E}_{k}$ holds for two
primary reasons. Firstly, the energy cascade results in the $k^{-2}$ decay of
the spectral energy $\widehat{E}_{k}$ due to formation of shock waves
Gurbatov_JFM_1997 . In the limit of vanishing viscosity $\delta\rightarrow 0$,
such decay extends up to $k\rightarrow\infty$ in which case the developed
shock waves render the system $C_{0}$ discontinuous. Secondly, the shape of
the spectral energy $\widehat{E}_{k}$ for $k\rightarrow k_{0}$ corresponds to
$k^{n}$, which is also the shape of initial energy spectral at time $t=0$.
Such argument corresponds to the concept of _permanence of large eddies_ in
hydrodynamic turbulence Batchelor1953theory , which in spectral space can be
written as,

$\widehat{E}_{k}(t)\approx\widehat{E}_{k}(t=0),~{}~{}\mathrm{as}~{}~{}k\rightarrow
k_{0}.$ (3.40)

Gurbatov _et al._ Gurbatov_JFM_1997 utilized a similar argument in the context
of Burgers’ turbulence. Combining the Eqs. (3.36)-(3.40), the integral length
scale $\ell$ is given by,

$\displaystyle\ell$
$\displaystyle\approx\begin{cases}\sqrt{\frac{n+1}{n-1}\left(\frac{k^{n-2}_{E}+k^{n-3}_{E}k_{0}+\cdots
k^{n-2}_{0}}{k^{n}_{E}+k^{n-1}_{E}k_{0}+\cdots k^{n}_{0}}\right)}&n\neq 1\\\
\sqrt{\frac{2\ln(k_{E}/k_{0})}{k^{2}_{E}-k^{2}_{0}}}&n=1.\end{cases},$ (3.41)

where we have used the simplifying approximation of $k_{\delta}\gg k_{E}$. We
note that the Eq. (3.41) indicates the dependence of $\ell$ and consequently
the energy decay law on $n$. In the present work, we perform numerical
simulations for an uncorrelated white noise (filtered) which corresponds to
$n=0$, and

$\ell\approx\frac{1}{\sqrt{k_{0}k_{E}}}.$ (3.42)

As a result of _permanence of large eddies_ , the decay of energy in the
initial regime of AWT is associated only with the decreasing $k_{E}$ or
increasing integral length scale $\ell$. Integrating Eq. (3.37) in the
spectral space and differentiating in time yields (for $n=0$ in the current
simulations),

$\displaystyle\frac{d\left\langle E^{(2)}\right\rangle}{dt}$
$\displaystyle=C_{1}\left(2\frac{dk_{E}}{dt}\left(1-\frac{k_{E}}{k_{\delta}}\right)+\left(\frac{k_{E}}{k_{\delta}}\right)^{2}\frac{dk_{\delta}}{dt}\right)$
(3.43) $\displaystyle\approx-\frac{2C_{1}}{k_{0}\ell^{3}}\frac{d\ell}{dt}.$
(3.44)

Above relation shows that derivation of the energy decay power law amounts to
finding the kinetic equations of the integral length scale $\ell$ and the
limiting wavenumber $k_{\delta}$. For TW and SW, $\ell$ remains constant by
definition. Consequently, the energy decay rate only depends on decrease of
wavenumber $k_{\delta}$ and the coefficient $C_{2}$ due to the thermoviscous
diffusion (cf. Eq. (3.48) ). However, for an ensemble of shock waves in AWT,
$\ell$ increases monotonically in time, as shown in figure 3.8$(a)$ due to the
coalescence of shock waves propagating in the same direction. At large times,
the domain consists of only two shock waves propagating in opposite
directions.

In the context of Burgers’ turbulence in an infinite one-dimensional domain,
Burgers Burgers1974nonlinear and Kida Kida_1979_JFM have derived the
appropriate asymptotic evolution laws for the integral length scale $\ell$
based on the dimensional arguments. However, in the present work, the
finiteness of the domain renders the asymptotic analysis infeasible. Our
numerical results indicate that $\ell\sim t^{1/3}$ ($k_{E}\sim t^{-2/3}$) for
randomly distributed shock waves at various $\mathrm{Re}_{L}$ values
considered, as shown in figure 3.8$(a)$. Equations (3.44) and (3.42) show that
such scaling is consistent with the observed energy decay law $\left\langle
E^{(2)}\right\rangle\sim t^{-2/3}$ thus validating the result in Eq. (3.44).
It is noteworthy that decay $k_{E}\sim t^{-2/3}$ is a result analogous to the
one discussed in Burgers’ turbulence Burgers1974nonlinear ; Kida_1979_JFM ;
Gurbatov_JFM_1997 considered in an infinite one-dimensional domain. Due to
infinitely long domain, the average distance between the shocks approaches
$1/k_{E}$ (not $\ell$) simply due to larger number of shocks in the domain
separated by the distance $1/k_{E}$ since $k_{E}$ corresponds to the largest
wavenumber carrying initial energy, thus implying that mean distance between
the shocks increases as $t^{2/3}$ as noted by Burgers Burgers1974nonlinear .

Based on the integral length scale, the Reynolds number $\mathrm{Re}_{\ell}$
can be defined as,

$\mathrm{Re}_{\ell}=\mathrm{Re}_{L}\ell,$ (3.45)

which captures the ratio of the diffusive time scale to the wave turbulence
time. Upon formation of shock waves, the perturbation energy decays due to
coalescence. Shock waves coalesce locally thus increasing the characteristic
separation between the shock waves thus causing $\ell$ to increase. In this
regime, the Reynolds number $\mathrm{Re}_{\ell}$ remains constant (figure
3.8$(b)$) which denotes that the ratio of shock coalescence time scale $(\ell
L^{*})/(a^{*}_{0}A_{\mathrm{rms}})$ and the diffusive time scale $(\ell
L^{*})^{2}/\nu^{*}_{0}$ remains constant. As the wave turbulence decays
further, $\ell\rightarrow L$ with continued decay of energy. Consequently,
$\mathrm{Re}_{\ell}$ also begins to decay.

####  IV.4 Kolmogorov length scale

For spectral energy $\widehat{E}_{k}\sim k^{-2}$ over the intermediate range
of wavenumbers, $k\in(k_{E},k_{\delta})$ (cf. figures 3.4 and 3.5), the Taylor
microscale can be estimated as,

$\lambda\sim\frac{1}{\sqrt{k_{E}k_{\delta}}},$ (3.46)

utilizing the Eq. (3.34). Equation (3.46) shows that $\lambda$, despite being
a dissipative scale, is not the smallest scale generated due to the energy
cascade. Analogous to the hydrodynamic turbulence, we define the Kolmogorov
length scale $\eta$ tennekes1972first as the smallest length scale generated
as a result of the acoustic energy cascade. The length scale $\eta$ can be
approximated by the balance of nonlinear steepening and energy dissipation,
i.e.,

$\frac{{A}^{2}_{\mathrm{rms}}}{\eta}\sim\delta\frac{A_{\mathrm{rms}}}{\eta^{2}},\quad\eta\sim\frac{\delta}{{A}_{\mathrm{rms}}},$
(3.47)

where $A_{\mathrm{rms}}$ is defined in Eq. (3.16). Figures 3.4 and 3.5
illustrates the integral length scale $\ell$ and the Kolmogorov length scale
$\eta$ in a typical AWT field. Visual inspection indicates $\ell\gg\eta$ which
is as expected. We note that $\eta$ and $1/k_{\delta}$ evolve in time
similarly, differing only by a constant value. For AWT, this is immediately
understandable since, Eq. (3.46) shows that $k_{\delta}\sim t^{-1/3}$ and Eq.
(3.47) shows that $\eta\sim t^{1/3}$ which implies $k_{\delta}\eta$ remains
constant when the energy decays.

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

Figure 3.9: Fully developed spectra of compensated energy $(a)$, spectral
energy flux $(b)$, and cumulative dissipation $(c)$ for TW at time instant
$t_{0}\approx 0.03$. Harmonics with wavenumbers such that $k\eta<1$ contain
all the energy. The spectral energy flux vanishes at $k\eta\approx 1$ thus
indicating numerical resolution of all the energy containing harmonics. The
marked regime $0.1<k\eta<1$ signifies the dissipation range. The constant
$C\approx 0.075$. (–) $A_{\mathrm{rms},0}=10^{-1}$; ($--$)
$A_{\mathrm{rms},0}=10^{-2}$; ($\cdots$) $A_{\mathrm{rms},0}=10^{-3}$

For TW and SW, the spectral energy given by Eq. (3.37) corresponds to the
degenerate case of $k_{0}=k_{E}=1$. For such a form of spectral energy, the
energy evolution (cf. Eq. (3.44)) changes to,

$\frac{d\left\langle
E^{(2)}\right\rangle}{dt}=\frac{1}{k_{0}}\frac{dC_{2}}{dt}\left(1-\frac{k_{0}}{k_{\delta}}\right)+\frac{C_{2}}{k^{2}_{\delta}}\frac{dk_{\delta}}{dt}.$
(3.48)

As shown in figure 3.7, the Taylor microscale $\lambda\rightarrow\sqrt{\delta
t}$. Consequently, for $k_{E}=k_{0}$ constant, Eq. (3.46) shows that
$k_{\delta}\sim t^{-1}$. Equation (3.48) shows that the decay of perturbation
energy is due to decay in $C_{2}$ and $k_{\delta}$. Our numerical results (cf.
figure 3.7) show that for TW and SW, $\left\langle E^{(2)}\right\rangle\sim
t^{-2}$ which suggests that $C_{2}\sim t^{-2}$ for $k_{\delta}\gg 1$ from Eq.
(3.48). Hence, the compensated energy spectrum $k^{2}\widehat{E}_{k}\sim
t^{-2}$ for both TW and SW indicating that dissipation
$\widehat{\mathcal{D}}_{k}$ remains active over all the length scales
$k>k_{0}$ while the energy decays.

Equation (3.47) shows that the Reynolds number based on the Kolmogorov length
scale or the shock thickness $Re_{\eta}=\eta\mathrm{Re}_{L}$ remains constant
in time,

$Re_{\eta}=\frac{\rho^{*}_{0}a^{*}_{0}L^{*}\eta
A_{\mathrm{rms}}}{\mu^{*}_{0}}=\frac{4}{3}+\frac{\gamma-1}{Pr}.$ (3.49)

Above relation shows that $Re_{\eta}=\mathcal{O}\left(1\right)$ indicating
that $\eta$ is the length scale at which diffusion dominates the nonlinear
wave steepening.

###  V Scaling of spectral quantities

In this section, we discuss the variation and scaling of the energy
$\widehat{E}_{k}$, the spectral energy flux $\widehat{\Pi}_{k}$, and the
cumulative dissipation $\sum_{k^{\prime}<k}\widehat{\mathcal{D}}_{k^{\prime}}$
for high amplitude TW, SW, and AWT cases utilizing the length scale analysis
presented in the previous sections. We show that the spectral energy
$\widehat{E}_{k}$ and the cumulative dissipation
$\sum_{k^{\prime}<k}\widehat{\mathcal{D}}_{k^{\prime}}$ for all the cases can
be collapsed on to a common structure versus the reduced wavenumber $k\eta$
however, the flux $\widehat{\Pi}_{k}$ lacks such a universality.

As discussed in previous section (cf. Eq. (3.48)), the decay of total energy
$\left\langle E^{(2)}\right\rangle$ and dissipation rate $\epsilon$ for TW is
given by,

$\left\langle E^{(2)}\right\rangle\sim t^{-2},~{}\mathrm{and}~{}\epsilon\sim
t^{-3},$ (3.50)

which are well known results for the Burgers’ equation as well Bec_PhyRep_2007
.

While the results in Eq. (3.50) are well known, we note that such power law
decay results in a universally constant structure of shock waves in the
spectral space, as shown in Figs. 3.9 and 3.10. Utilizing the estimate of
Kolmogorov length scale $\eta$ given in Eq. (3.47), the energy dissipation
rate $\epsilon$ and the Kolmogorov length scale $\eta$ can be related as,

$\epsilon\sim\frac{{A_{\mathrm{rms}}^{3}}}{\ell},~{}\mathrm{and}~{}\eta\sim\frac{\delta}{\left(\epsilon\ell\right)^{1/3}}.$
(3.51)

Hence, the energy spectrum $\widehat{E}_{k}$ can be written in the following
collapsed form (figure 3.9$a$).

$\widehat{E}_{k}k^{2}\epsilon^{-2/3}\ell^{1/3}\sim CF(k\eta).$ (3.52)

In Eq. (3.52), the integral length scale $\ell$ is used for making the left
hand expression dimensionless. For TW and SW, the integral length scale $\ell$
remains constant by definition ($\ell=L$). Hence, $C$ in Eq. (3.52) is
constant and can be attributed to the Kolmogorov’s universal equilibrium
theory for hydrodynamic turbulence. $F(.)$ is a function which decays as the
reduced wavenumber $k\eta$ increases to 1. From the numerical simulations for
cases listed in Table 3.1 we obtain,

$C\approx 0.075.$ (3.53)

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

Figure 3.10: Fully developed spectra of compensated energy $(a)$, spectral
energy flux $(b)$, and cumulative dissipation $(c)$ for SW averaged over one
time cycle after $t_{0}\approx 0.04$. Harmonics with wavenumbers $k\eta<1$
contain all the energy. The spectral energy flux vanishes at $k\eta\approx 1$
thus indicating numerical resolution of all the energy containing harmonics.
The constant $C\approx 0.075$. (–) $A_{\mathrm{rms},0}=10^{-1}$; ($--$)
$A_{\mathrm{rms},0}=10^{-2}$; ($\cdots$) $A_{\mathrm{rms},0}=10^{-3}$

Scaling of $\widehat{\Pi}_{k}$ with the energy dissipation rate $\epsilon$
shows the relative magnitude of spectral energy flux compared to the energy
dissipation. For increasing Reynolds numbers $\mathrm{Re}_{L}$, we note that
$\widehat{\Pi}_{k}/\epsilon$ increases but still remains less than 1 in the
energy transfer and dissipation range, as shown in figure 3.9$(b)$. This
highlights the primary difference between energy spectra of nonlinear acoustic
waves and hydrodynamic turbulence, in which, the energy transfer range does
not exhibit viscous dissipation pope2000turbulent . However, in nonlinear
acoustics, the dissipation occurs over all the smaller length scales which do
not contain energy initially (figure 3.5$(b)$). Moreover, for $k\eta\approx
0.1$, the flux $\widehat{\Pi}_{k}$ rapidly approaches to zero. In the regime
$k\eta>0.1$, scaled cumulative dissipation
$\sum_{k^{\prime}<k}\widehat{\mathcal{D}}_{k^{\prime}}/\epsilon\rightarrow 1$
as $k\eta\rightarrow 1$.

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

Figure 3.11: Fully developed spectra of compensated energy, $(a)$, spectral
energy flux $(b)$, and cumulative dissipation $(c)$ against scaled wavenumber
$k\eta$ for the randomly initialized broadband noise (AWT) cases with
$A_{\mathrm{rms},0}$ and $\nu_{0}$ listed in Table 3.1 at dimensionless time
$t/\tau=\tau_{0}\approx 6\times 10^{-4}$. The marked regime $0.1<k\eta<1$
signifies the dissipation range. (–) $A_{\mathrm{rms},0}=10^{-1}$; ($--$)
$A_{\mathrm{rms},0}=10^{-2}$; ($\cdots$) $A_{\mathrm{rms},0}=10^{-3}$

Such functional forms of spectral energy, spectral energy flux, and cumulative
dissipation can also be realized for the SW case. At later times, the
nonlinear evolution results in two opposite traveling shock waves which
collide with each other twice in one time period. Such collisions cause
instantaneous peaks in the dissipation rate $\epsilon$ and corresponding
oscillations in the Taylor microscale $\lambda$, as shown in figure 3.7.
However, the total energy $\left\langle E^{(2)}\right\rangle$ decays
monotonically by definition. In the spectral space, such collisions generate
periodic oscillations in the spectral energy flux $\widehat{\Pi}_{k}$, as
shown in figure 3.6. Averaging over one such time cycle yields the energy
spectra forms similar to that for TW, as shown in figure 3.10. Such cycle
averaging is allowed since the total energy $\left\langle
E^{(2)}\right\rangle$ and the dissipation rate $\epsilon$ decay such that
averaged behavior is identical to the one of traveling waves. Furthermore, the
value of the constant $C$ is identical for SW. We further note that for the
case with lowest Reynolds number $\mathrm{Re}_{L}$ ($\nu_{0}=1.836\times
10^{-5}$ and $A_{0,\mathrm{rms}}=10^{-3}$), the spectra exhibit energy for
$k\eta>1$ (figure 3.10$(a)$ since the Eq. (3.47) underpredicts $\eta$. This
suggests that the nonlinear spectral energy transfer is small compared to the
spectral dissipation, as shown by figure 3.10$(b)$.

As discussed in previous sections, the decay phenomenology of AWT is different
from that of TW and SW. Typical acoustic field
$u^{\prime}(x,t),p^{\prime}(x,t)$ for a randomly initialized perturbation at a
time after shock formation is shown in figure 3.4$(a)$. The velocity field
corresponds to randomly positioned shocks connected with almost straight slant
lines (expansion waves) and the pressure field with identical distribution of
shocks but connected with horizontal lines. Shocks traveling in the same
direction collide inelastically and coalesce, while those traveling in
opposite directions pass through. As discussed previously, the integral length
scale $\ell$ defines the average distance between the adjacent shock traveling
in the same direction. Due to gradual coalescence of the shocks, $\ell$
increases in time. Moreover, as $t\rightarrow\infty$, it is obvious that two
opposite traveling shocks remain in the domain and $\ell\rightarrow L$. We
note that such behaviour is similar to the Burgers’ turbulence Kida_1979_JFM .
Figure 3.11 shows the fully developed compensated spectra at scaled
dimensionless time $t/\tau=t_{0}\approx 6\times 10^{-4}$. For AWT, the
compensated energy spectrum $\widehat{E}_{k}k^{2}\epsilon^{-2/3}\ell^{1/3}$
defined in Eq. (3.52) does not remain constant in the energy transfer range of
wavenumbers due to decay laws of energy and dissipation derived in the
previous section. Moreover, for lowest $\mathrm{Re}_{L}$ case, the spectra
exhibit energy for $k\eta>1$ (figure 3.11$(a)$) due to underprediction of
$\eta$ obtained via balancing of nonlinear wave propagation and thermoviscous
dissipation effects. The dissipation acts at large length scales also in the
lowest $\mathrm{Re}_{L}$ case. Consequently, the spectral energy flux
$\widehat{\Pi}_{k}$ is very small compared to dissipation $\epsilon$ and the
length scale $\eta$ is primarily governed by diffusion only.

###  VI Concluding remarks

We have studied the spectral energy transport and decay of finite amplitude
planar nonlinear acoustic perturbations governed by fully compressible 1D
Navier-Stokes equations through shock-resolved direct numerical simulations
(DNS) focusing on propagating single harmonic traveling wave (TW), standing
wave (SW), and randomly initialized Acoustic Wave Turbulence (AWT). The
maximum entropy perturbations scale as $p^{\prime 2}$ for normalized pressure
perturbation $p^{\prime}\sim\mathcal{O}\left(10^{-3}-10^{-1}\right)$.
Consequently, the second order nonlinear acoustic equations are adequate to
derive physical conclusions about spectral energy transfer in the system.
Utilizing the second order equations, we derived the analytical expression for
corrected energy corollary for finite amplitude acoustic perturbations
yielding infinite order correction term in the perturbation energy density. We
have shown that the spatial average of the corrected perturbation energy
density can be classified as a Lyapunov function for the second order
nonlinear acoustic system with strictly monotonic behaviour in time.

Utilizing the corrected energy corollary, we derived the expressions for
spectral energy, spectral energy flux, and spectral dissipation, analogous to
the spectral energy equation studied in hydrodynamic turbulence. Utilizing the
spectral expressions, we performed theoretical study of three possible length
scales characterizing a general nonlinear acoustic system, namely, the
integral length scale $\ell$, the Taylor microscale $\lambda$, and the
Kolmogorov length scale $\eta$.

In traveling waves (TW) and standing waves (SW), $\ell$ remains constant in
the decaying regime. 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. The Kolmogorov scale increases linearly in time
($\eta\sim t$) in the decaying regime. 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)$.

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}$, similar to the Burgers’ turbulence Burgers1974nonlinear . While,
various cases for AWT qualitatively collapse with the scaling
$\widehat{E}_{k}k^{2}\epsilon^{-2/3}\ell^{1/3}$, quantitative scaling can only
be obtained utilizing a statistically stationary ensemble of shock waves
combined with random forcing, which falls beyond the current scope.

##  Chapter 4 Thermoacoustic nonlinear energy cascade

The contents of this chapter were previously published by Gupta, Lodato, Scalo
in _Journal of Fluid Mechanics_ GuptaLS_JFM_2017 and have been reported here
in abridged form with minor modifications.

###  I Introduction

In previous chapter, we discussed the spectral energy cascade and decay in
initialized nonlinear acoustic waves. In this chapter, we discuss
thermoacoustically sustained spectral energy cascade in a canonical
thermoacoustically unstable resonator exhibiting high amplitude (macrosonic)
thermoacoustic waves at limit cycle. A comprehensive nonlinear theoretical and
high-fidelity modeling approach is adopted to accurately describe macrosconic
thermoacoustic waves. The travelling-wave looped resonator, inspired by Yazaki
_et al._ Yazaki_PhysRevLet_1998 ’s experimental setup but geometrically
optimized via linear theory developed by Rott Rott_ZAMP_1969 ; Rott_ZAMP_1973
; LinSH_JFM_2016 , has been designed to maximize the growth rate of the quasi-
travelling-wave second harmonic and thus achieve rapid shock wave formation.
Yazaki _et al._ ’s Yazaki_PhysRevLet_1998 looped configuration allows quasi-
travelling-wave acoustic phasing which facilitates faster nonlinear energy
cascade compared to standing wave resonators BiwaEtAl_JASA_2014 . It is shown
that the energy content in spectral domain resembles the equilibrium energy
cascade observed in turbulence, similar to the spectral energy distribution of
an ensemble of acoustic waves interacting nonlinearly among each other
Nazarenko_2011_WT ; Zakharov_2012_WT . As demonstrated by the numerical
simulation data and companion low-order nonlinear modeling, thermoacoustically
sustained shock waves exhibit inter-scale energy transfer dynamics analogous
to Kolmogorov’s equilibrium hydrodynamic turbulent energy cascade
kolmogorov1941local . Throughout, the results obtained via the proposed
nonlinear model are verified and compared with fully compressible high-
fidelity Navier-Stokes simulations.

The development of an accurate nonlinear thermoacoustic wave propagation
theory and modeling framework warrants the support of high-fidelity numerical
simulations. A high-order spectral difference numerical framework
kopriva1996conservative ; kopriva1996conservative2 ; sun2007high ;
Jameson_JSC_2010 combined with an artificial Laplacian viscosity persson:06
shock-capturing scheme has been adopted for the present study. Moreover, the
computational setup has been reduced to a minimal-unit (or single-pore)
configuration, as done by Rahman _et al._ RahmanEtAl_IntJHMT_2017 , to reduce
the computational cost and ensure the maximum possible numerical resolution in
the direction of shock propagation for a given number of discretization
points, or degrees of (numerical) freedom. In spite of this choice, full
resolution of the propagating shocks was still not attainable with the
available resources.

###  II Problem Formulation

####  II.1 Design of the minimal-unit model

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

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

Figure 4.1: Two-dimensional axially periodic computational setup for minimal-unit simulations (top, not to scale), geometrical parameters (bottom left), and qualitative illustration of the equivalent experimental build-up of a variable-area looped thermoacoustic resonator (bottom right). Top figure: (– –), adiabatic/sip conditions; (–), isothermal/no-slip conditions. The minimal unit is traced in the bottom right with dashed lines (– –). The listed geometrical parameters provide sub-optimal growth rates across all the values of $T_{H}$, the boxed quantities have been determined via optimization (cf. figure II.2 and table 4.1). Table 4.1: Thermodynamic parameters for base state. Base pressure and cold-side values of density and temperature are set to reference values: $P_{0}=P_{\mathrm{ref}}$, $\rho_{c}=\rho_{\mathrm{ref}}$, $T_{C}=T_{\mathrm{ref}}$. $P_{0}$ |  $\rho_{C}$ |  $T_{C}$ |  $T_{H}$
---|---|---|---
101325 Pa |  1.176 kg/m3 |  300 K |  400 K, 450 K,
|  |  |  500 K, 550 K

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

$\delta_{k}=\sqrt{\frac{2\nu}{\omega Pr}},$ (4.4)

$\nu$ is the kinematic viscosity, $Pr$ is the Prandtl number, and $\omega$ is
the angular frequency of the unstable mode. This results in a porosity
$h_{b}/h_{c}=0.91$. The area ratio $h_{d}/h_{c}$ and the length $\ell_{c}$
have been chosen to yield sub-optimal values of growth rates across all hot-
side temperature settings (figure II.2), assuring high enough thermoacoustic
instability to achieve rapid shock wave formation. This optimization has been
carried out with the system-wide numerical approach developed by Lin _et al._
LinSH_JFM_2016 . Further details of the linear stability analysis are given in
GuptaLS_JFM_2017 .

####  II.2 Navier-Stokes calculations

##### Governing Equations

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

$\displaystyle\frac{\partial}{\partial t}\left(\rho\right)$
$\displaystyle+\frac{\partial}{\partial x_{j}}\left(\rho u_{j}\right)=0,$
(4.5a) $\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}+\delta_{1i}f_{B},$ (4.5b)
$\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),$ (4.5c)

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

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

$T_{H}$ |  400 K |  450 K |  500 K |  550 K
---|---|---|---|---
$\alpha$ |  1.79 $\mathrm{s}^{-1}$ |  8.64 $\mathrm{s}^{-1}$ |  14.51 $\mathrm{s}^{-1}$ |  19.43 $\mathrm{s}^{-1}$
Figure 4.2: Iso-contours of thermoacoustic growth rate, $\alpha$, of the
second harmonic versus resonator area ratio, $h_{d}/h_{c}$, and length
$\ell_{c}$ (see figure 4.1) for hot-side temperatures $T_{H}=400~{}\mathrm{K}$
($a$), $T_{H}=450~{}\mathrm{K}$ ($b$), $T_{H}=500~{}\mathrm{K}$ ($c$), and
$T_{H}=550~{}\mathrm{K}$ ($d$). (–), $\alpha>0$; (- -), $\alpha<0$; ( ), sub-
optimal values $\ell_{c}=0.59$ m and $h_{d}/h_{c}=1.28$ chosen for present
investigation (figure 4.1). Growth rate values for each $T_{H}$ at sub-optimal
geometry are listed in the table.

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

$f_{B}=\frac{\alpha}{\Delta t}\left(\dot{m}-\dot{m}_{0}\right),$ (4.8)

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

$\displaystyle\tau_{ij}=2\mu\left[S_{ij}-\frac{1}{3}\frac{\partial
u_{k}}{\partial x_{k}}\delta_{ij}\right],\quad q_{j}$ $\displaystyle=$
$\displaystyle-\frac{\mu\,C_{p}}{Pr}\frac{\partial T}{\partial x_{j}},$ (4.9a)

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

Figure 4.3: Illustration of selection portions of the unstructured mesh near
duct interfaces for Grid A with total number of elements
$N_{\mathrm{el}}=2004$.

The fully compressible high-fidelity Navier-Stokes calculations have been
carried out with the discontinuous finite element sd3DvisP solver, an MPI