On Linear and Nonlinear Acoustics in Stratified, Variable-Area Ducts and Atmospheres, and Lighthill’s Proposition

The propagation of finite-amplitude acoustic waves, and especially the creation of shocks within them, are of interest for physical problems across a wide range of scales – from the collapse of sonoluminescent bubbles (Lin & Szeri, 2001) to the eruption of stellar envelopes in pre-supernova outbursts (Smith, 2013). We shall consider the case of quasi-one-dimensional oscillations of an otherwise hydrostatic, stratified background state, allowing for variations in the cross-sectional area of the channel or atmosphere.

An important point of reference is the simplest limit in which the cross-sectional area is constant and the background state is uniform. Riemann (1860) solved this case exactly, up to the moment of shock formation within the flow, in terms of invariants conserved along sound fronts (characteristic trajectories). For the subsequent evolution, Whitham (1975, building on the work of Landau 1945) developed an elegant analytical formalism to very accurately predict the strength and trajectory of a weak shock. Specifically, Whitham’s ‘equal area’ construction applies to shocks made by a disturbance traveling in one direction only, i.e., a simple wave.

Can these solutions be adapted to the general problem? In chapter 2 of Waves in Fluids (Lighthill, 1978), M. J. Lighthill proposed that they can. Considering the equivalent of simple waves, Lighthill proposed an invariant invariant should be approximately conserved along characteristic trajectories – approximately, because the argument relies on low-amplitude perturbations and on the high-frequency limit in which wave reflection is negligible. Coupled with the pressure-velocity relation appropriate to simple waves, this conservation law allows one to transform a traveling-wave problem in the general problem to its analog in the simple limit, for which Riemann’s and Whitham’s approaches provide a solution.

This method is appealing as an analytical technique even though it is not exact. Low-amplitude, high-frequency waves are precisely the type that sample variations in their environment before forming a shock. Furthermore, such waves are challenging to treat with direct numerical simulations because accuracy demands hundreds or thousands of numerical zones per wavelength. A shock forms only after many wavelengths of propagation if the wave amplitude is low.

While Lighthill’s proposal has not, to our knowledge, been re-examined, similar approaches have certainly been considered. High-frequency scalings, familiar from WKB theory, motivate Subrahmanyam et al. (2001)’s proposed exact solutions for specific problems in linear acoustics. In our own previous work (Ro & Matzner, 2017), we reinvented aspects of Lighthill’s proposal, including his criterion for shock formation (his eq. 254), by positing that the high-frequency scalings apply to individual characteristics. Other authors, such as Lin & Szeri (2001), Tyagi & Sujith (2003), and Tyagi & Sujith (2005), also derive shock formation criteria that coincide with Lighthill’s, as Tyagi & Sujith note.

To develop the theory further, we shall express linear acoustical motions in terms of wave amplitudes that are conserved along characteristics except for the effects of reflection. These amplitudes are modified versions of Riemann’s invariants; they are equivalent to Lighthill’s proposed invariant in a linear travelling wave; and they are related in a simple way to the wave luminosity. When the linearized amplitude equation is applied to a nonlinear wave, errors are suppressed by factors of the gradients in the background. This reflects the fact that wave amplitudes reduce to Riemann’s invariants in the uniform case. Together, these points prove the validity of Lighthill’s proposal and the exactness of his shock formation criterion in certain contexts, while also providing a convenient way to evaluate wave reflection.

We start with adiabatic, quasi-one-dimensional fluid equations in Lagrangian form, and obtain Riemann’s invariants $J_{\pm}$ for the simple case of planar motion of a uniform fluid. We then derive the variation of $J_{\pm}$ in the general problem. Linearizing this result and introducing wave amplitudes ${\cal R}_{\pm}$, we arrive at the linear amplitude evolution equation. We relate these amplitudes to the wave luminosity, and consider the dynamics of wave reflection. We then examine the nonlinear errors incurred when the linearized equation is applied to waves of finite amplitude, showing that these are small in the high-frequency limit. Moreover, we show that Lighthill’s shock formation condition is exact in certain cases, and essentially exact for shocks forming within high-frequency waves. Then, we demonstrate that nonlinear solutions are available for a certain class of problems. Finally, we conduct numerical experiments to demonstrate the application to traveling waves in spherical and strongly stratified environments.

Our quasi-one-dimensional equations incorporate a couple assumptions: first, that motions transverse to the gradient of $r$ can be ignored; and second, that displacements are small enough ($\delta r\ll r_{0}$) that smooth functions of $r$ can be replaced with their values at $r_{0}$: for instance, $g=g_{0}$ and $A=A_{0}$. In addition, we neglect dissipative effects such as thermal diffusion; see Ro & Matzner (2017) for a justification in the context of stellar outbursts.

We start with the conservation of mass,

$$d_{t}\ln\rho+A^{-1}\partial_{r}Av=0,$$

momentum,

$$d_{t}v+\frac{1}{\rho}\partial_{r}P=-g,$$

as well as background hydrostatic equilibrium,

$$\frac{1}{\rho}_{0}\partial_{r_{0}}P_{0}=-g_{0}.$$

We consider the fluid to be adiabatic and for simplicity we adopt an ideal gas equation of state,

$$P/P_{0}=(\rho/\rho_{0})^{\gamma}=x^{\gamma},$$

where $\gamma$ may be a function of $r_{0}$.

We convert the spatial gradient ($\partial_{r}$) into gradients with respect to the initial radius ($\partial_{r_{0}}$) using $A\rho\,dr=A_{0}\rho_{0}\,dr_{0}$, which implies

$$\rho^{-1}\partial_{r}=(A/A_{0})\rho_{0}^{-1}\partial_{r_{0}}\simeq\rho_{0}^{-1}\partial_{r_{0}}$$

using $A\simeq A_{0}$. At the same time, we switch spatial variables from $r$ to $r_{0}$, allowing $d_{t}\rightarrow\partial_{t}$. Our equation for mass conservation becomes (using $r^{-1}\simeq r_{0}^{-1}$),

$$\partial_{t}\ln(\rho/\rho_{0})+(\rho/\rho_{0})\partial_{r_{0}}v+v\,\partial_{r_{0}}\ln A=0.$$

The momentum equation becomes

$$\partial_{t}v+\frac{1}{\rho_{0}}\partial_{r_{0}}\delta P=0$$

or, using $P=x^{\gamma}P_{0}$ and writing out the derivative,

$$\partial_{t}v+\frac{\gamma P_{0}}{\rho_{0}}x^{\gamma-1}\partial_{r_{0}}x+\frac{\partial_{r_{0}}P_{0}}{\rho_{0}}(x^{\gamma}-1)+c_{0}^{2}(\ln x)\partial_{r_{0}}\ln\gamma=0,$$

which we re-write, using $\gamma P_{0}=\rho_{0}c_{0}^{2}$, $\partial_{r_{0}}P_{0}=-\rho_{0}g_{0}$, and $c_{0}^{2}x^{\gamma-1}=c^{2}$, as

$$\partial_{t}v+c^{2}\partial_{r_{0}}x=(x^{\gamma}-1)g_{0}-c_{0}^{2}(\ln x)\partial_{r_{0}}\ln\gamma.$$

The mass equation takes in the form

$$\partial_{t}x+x^{2}\partial_{r_{0}}v=-vx\,\partial_{r_{0}}\ln A,$$

where we have multiplied through by $x$ to put the time derivatives in the same form as the momentum equation. Combining the momentum and mass equations,

$$\partial_{t}\mathbf{u}+\mathsfbi{B}\,\partial_{r_{0}}\mathbf{u}=\mathbf{s},$$

(1)

where

$$\mathbf{u}=\left(\begin{array}[]{l}v\\ x\end{array}\right),~{}~{}\mathsfbi{B}=\left[\begin{array}[]{ll}0&c^{2}\\ x^{2}&0\end{array}\right],~{}~{}{\rm and}~{}~{}\mathbf{s}=\left(\begin{array}[]{c}(x^{\gamma}-1)g_{0}-c_{0}^{2}(\ln x)\partial_{r_{0}}\ln\gamma\\ -vx\,\partial_{r_{0}}\ln A\end{array}\right).$$

Note that our approximation of small radial perturbation only affects the form of $\mathbf{s}$.

In the simplest limit of a homogeneous background fluid and constant area, $g_{0}=0$ and $\partial_{r}\ln A=0$ so that $\mathbf{s}=0$. To derive the Riemann invariants, we seek a function $J(x,v)$ and a Lagrangian speed $\lambda$ that satisfy $\partial_{t}J+\lambda\partial_{r_{0}}J=0$.

Writing out these partial derivatives in subscript notation, $J_{t}=J_{x}x_{t}+J_{v}v_{t}=(J_{v},J_{x})\mathbf{u}_{t}$ and likewise $J_{r_{0}}=(J_{v},J_{x})\mathbf{u}_{r_{0}}$. The equation we wish to solve, $J_{t}+\lambda J_{r_{0}}=0$, becomes $(J_{v},J_{x})\mathbf{u}_{t}+\lambda(J_{v},J_{x})\mathbf{u}_{r_{0}}=0$; using (1) to eliminate the time derivative,

$$(J_{v},J_{x})(\mathsfbi{B}-\lambda\mathsfbi{I})\mathbf{u}_{r_{0}}=0.$$

For this to be true for any $\mathbf{u}_{r_{0}}$ requires $(J_{v},J_{x})(\mathsfbi{B}-\lambda\mathsfbi{I})=0$. Taking the transpose,

$$(\mathsfbi{B}^{T}-\lambda\mathsfbi{I})(J_{v},J_{x})^{T}=0,$$

which is solved if $\lambda$ and $(J_{v},J_{x})^{T}$ are an eigenvalue and corresponding eigenvector of $\mathsfbi{B}^{T}$, respectively.

The eigenvalues are $\lambda_{\pm}=\pm cx$. These represent the Lagrangian speeds of sound fronts moving to the right or left though space at speed $v\pm c$, because $A_{0}\rho_{0}(dr_{0}/dt)=A\rho(dr/dt-v)$ along any trajectory, and $A\rightarrow A_{0}$ for small perturbations. (We refer to $\lambda_{\pm}$ variously as the wave speed, propagation speed, or characteristic speed, and to trajectories $C^{\pm}$ obeying $\dot{r}_{0}=\lambda_{\pm}$ as characteristics or sound fronts.)

The corresponding eigenvectors are $(J_{v},J_{x})_{\pm}=(1,\pm c/x)$. The functions that have these partial derivatives are the usual Riemann quantities,

$$J_{\pm}=v\pm\int c\frac{dx}{x}=v\pm\frac{c}{q},$$

which are invariants when $\mathbf{s}=0$. Having reconstructed the conservation of the Riemann quantities for this restricted problem, we now turn to how $J_{\pm}$ vary in the more general context.

Lifting the restrictions of planar flow and a uniform fluid, Riemann’s quantity $J$ is no longer invariant. Its time variation along a characteristic is

$$\dot{J}={\sum_{f\in\{v,x,c_{0},q\}}}J_{f}(f_{t}+\lambda f_{r_{0}})=(J_{v},J_{x})\cdot\mathbf{s}+\lambda J_{c_{0}}\partial_{r_{0}}c_{0}+\lambda J_{q}\partial_{r_{0}}q.$$

The first term on the right hand side arises from non-zero gravity and changes of area (components of $\mathbf{s}$), while the second and third terms account for gradients of the fluid properties; these contain only radial derivatives because $c_{0}$ and $\gamma$ are constant within each fluid element.

We prefer to work with radial derivatives for later convenience, which we compute as $d_{r_{0}}J=\dot{J}/\lambda$. Written out for outward and inward waves,

$$d_{r_{0}}^{\pm}{J_{\pm}}=\pm\frac{x^{\gamma}-1}{x^{(\gamma+1)/2}}\frac{g_{0}}{c_{0}}\mp\frac{c_{0}^{2}\ln x}{cx}\partial_{r_{0}}\ln\gamma-\frac{v}{x}\partial_{r_{0}}\ln A\pm\frac{c}{q}\partial_{r_{0}}\ln c_{0}\mp\frac{c}{q}(1-q\ln x)\partial_{r_{0}}\ln q.$$

Part of the change in $J_{\pm}$ is acoustical, but part is due to the radial gradient of the unperturbed state. To isolate the acoustical portion, we subtract $\partial_{r_{0}}J_{0\pm}$ to obtain

$$d_{r_{0}}^{\pm}\delta J_{\pm}=\pm\frac{x^{\gamma}-1}{x^{(\gamma+1)/2}}\frac{g_{0}}{c_{0}}\mp\frac{c_{0}^{2}\ln x}{cx}\partial_{r_{0}}\ln\gamma-\frac{v}{x}\partial_{r_{0}}\ln A\pm\frac{\delta c}{q}\partial_{r_{0}}\ln c_{0}\mp\left(\frac{\delta c}{q}-c\ln x\right)\partial_{r_{0}}\ln q.$$

(2)

Equation (2) is exact within the quasi-one-dimensional approximation we have adopted. It therefore provides a point of comparison for the linearized wave amplitude equation derived below.

At this point we focus on linear perturbations and expand to leading order in $\delta x=x-1$:

$$d_{r_{0}}^{\pm}\delta J_{\pm}=\pm\gamma\,\delta x\frac{g_{0}}{c_{0}}\mp c_{0}\delta x\,\partial_{r_{0}}\ln\gamma-v\partial_{r_{0}}\ln A\pm\frac{\delta c}{q}\partial_{r_{0}}\ln c_{0}\mp\left(\frac{\delta c}{q}-c_{0}\,\delta x\right)\partial_{r_{0}}\ln q.$$

The first term, due to gravity, can be rewritten using $\gamma g_{0}/c_{0}=-c_{0}\partial_{r_{0}}\ln P_{0}$.

We now express the perturbations in terms of $\delta J_{\pm}$ using $v=(\delta J_{+}+\delta J_{-})/2$ and $(\delta c)/q=(\delta J_{+}-\delta J_{-})/2=c_{0}\delta x+{\cal O}(\delta x^{2})$. The term involving $\partial_{r_{0}}\ln q$ is zero to leading order, and the rest can be written compactly, using $\gamma P_{0}=\rho_{0}c_{0}^{2}$ and $\delta J_{+}-\delta J_{-}=\pm(\delta J_{\pm}-\delta J_{\mp})$, as

$$d_{r_{0}}^{\pm}\delta J_{\pm}=\delta J_{\pm}\,\partial_{r_{0}}\ln\frac{1}{\sqrt{\rho_{0}c_{0}A}}+\delta J_{\mp}\,\partial_{r_{0}}\ln\sqrt{Z}$$

(3)

where $Z=\rho_{0}c_{0}/A$ is the characteristic acoustic impedance. The first term represents the inward or outward propagation of an acoustic disturbance, while the second term encapsulates a conversion between it and the counter-propagating disturbance in the process of reflection, which is catalyzed by gradients of $Z$.

Equation (3) is simplified in terms of a new variable, the wave amplitude:

$${\cal R}_{\pm}\equiv\sqrt{A\rho_{0}c_{0}}~{}\delta J_{\pm}.$$

(4)

In linear theory this quantity evolves along $C^{\pm}$ only due to reflection:

$$d_{r_{0}}^{\pm}{\cal R}_{\pm}={\cal R}_{\mp}\,\partial_{r_{0}}\ln\sqrt{Z},$$

(5)

so it is effectively conserved so long as reflection is negligible. This wave amplitude equation provides a means to solve linear acoustics in the general problem by the method of characteristics. As we shall see, it also provides a connection to Lighthill’s approximate invariant and to the wave luminosity. Moreover, it provide a useful means to approximate nonlinear acoustics as well.

Our applications of interest involve the nonlinear effects of wave self-distortion, shock formation, and shock evolution. So, we now consider wave motion in the nonlinear regime. We begin by demonstrating the existence of nonlinear solutions to certain problems in the high-frequency limit. We then evaluate the nonlinear error associated with Lighthill’s proposal that the linear invariant may be used to approximate nonlinear problems. Finally, we consider the validity of Lighthill’s shock formation criterion.

We present three numerical tests covering the phases before, at, and after shock formation. We use the hydrodynamics code FLASH4.6 (Fryxell et al., 2000) set to use the HLLC Riemann solver with fifth-order reconstruction. The spherical simulation domain ($A=4\pi r^{2}$) spans $1\leq r/r_{0,i}\leq 3$ with uniform radial resolution. The initial fluid pressure is constant, but the density increases as $\rho_{0}\propto r^{4}$ such that $Z$ is constant. We adopt an ideal gas equation of state with $\gamma=5/3$. Variations in pressure, density, and velocity related by $v=\delta c/q$ are applied at the inner boundary to generate an outward traveling wave.

First, we test the prediction that there is no wave reflection in a region of uniform impedance, so wave luminosity should be conserved along characteristics. Here, we introduce one period of a sinusoidal wave (which starts in compression), and set its velocity amplitude to $10^{-3}c_{0}$ at the inner boundary. We choose $\omega$ so the acoustic wavelength is twice the density scale height (twice $r/4$) at the inner boundary; however, the drop in $c_{0}$ causes the wavelength to be significantly smaller than $r/4$ at the outer boundary. Figure 2 shows the normalized instantaneous wave luminosity $L_{w+}$ from four simulations with resolution that increases relative to a base resolution of $N=8192$ cells. Wave maxima and minima both display peaks of $L_{w+}$, which are slightly higher at the maxima. We see that the deviations from constant luminosity in the wave maxima (minima) decreases with increasing spatial resolution to a maximum of $0.02\%$ ($0.4\%$) for the highest resolution simulation. The numerical results support the conclusion that luminosity is conserved along characteristics, across nearly two orders of magnitude variations in background density, for a wave traveling in a fluid with constant impedance.

Next, we check Lighthill’s proposition that mildly nonlinear wave distortion can be accurately predicted using the linear adiabatic invariant. Using the same background profile and initialization strategy, and adopting the highest resolution, we increase the amplitude of the sinusoidal wave to $10^{-2}c_{0}$ while doubling the wave frequency, so that shock formation occurs within the domain. In Figure 4, we compare numerical and analytical predictions for the form of the wave pulse at the time of shock formation. The two agree very well, at least at the level estimated in § 6.2, despite the fact that the wave has traversed over two orders of magnitude in background density from the point of initialization.

For a final demonstration, we show in Figure 3 that the adiabatic invariance of ${\cal R}_{\pm}$ can be used to adapt Whitham’s equal-area method to weak shock evolution within a varying background. Here, a wave train is introduced for which the initial waveform is a ‘reverse saw-tooth’ with a linear rise and sudden decline in ${\cal R}_{+}(t_{0})$. The waveform converts to a ‘forward saw-tooth’ at the point of shock formation. After shock formation, the wave dissipates in a manner we discuss in detail in a companion paper (Matzner & Ro, 2020). Notably, the leading or ‘head’ shock decays distinctly more slowly than the ‘internal’ shocks that form within the wave train.

Our primary findings are as follows:

• -

  The linearized equations of hydrodynamics are compactly expressed (eq. 5) in terms of wave amplitudes ${\cal R}_{\pm}$, which interact in the presence of variations of the impedance.

• -

  In the context of ‘quasi-simple’ waves (§ 5.1) – those in which one wave family is absent and not regenerated by reflection – the conservation of ${\cal R}_{\pm}$ is consistent with the invariant proposed by Lighthill (1978), and also with the conservation of wave luminosity (§ 5.2).

• -

  Familiar results about the dynamics of reflection follow directly from the amplitude equation (§5.3). In particular, wave reflection is strongly suppressed in the high-frequency limit.

• -

  Lighthill’s proposition amounts to the statement that linear conservation laws may be extrapolated into the mildly nonlinear regime, providing a simple means to calculate wave self-distortion and shock formation, as well as weak shock dissipation by Whitham’s equal-area method. We validate this method and delimit its validity. In the high-frequency limit, the nonlinear error is minimal at the point of shock formation (§ 6.2). Because there is no nonlinear error at the wave node, Lighthill’s criterion for shock formation is exact so long as the shock forms at a node and reflection is negligible, as it is for ‘head’ shocks, and for any shocks generated within high-frequency waves (§ 6.3).

• -

  Going beyond linearized theory, we find that nonlinear solutions are available in a restricted class of problems (§ 6.1). In the example considered, the novel element is a nonlinear correction to the wave amplitude, which leads to a modified phase speed.

In a companion work (Matzner & Ro, 2020), we use these findings to predict the evolution of shock strength and shock dissipation in the context of super-Eddington outbursts of evolved massive stars. There, we provide useful formulae for the shock-deposited heat, and we evaluate the potential of head shocks to eject matter from the stellar surface.

This work was funded by an NSERC Discovery Grant (CDM) and by the Gordon and Betty Moore Foundation through Grant GBMF5076 (SR).