Nonlinear Acoustics and Shock Formation in Lossless Barotropic Green--Naghdi Fluids

The equations of motion of lossless compressible nonclassical fluids under the so-called Green--Naghdi theory are considered for two classes of barotropic fluids: (\textit{i}) perfect gases and (\textit{ii}) liquids obeying a quadratic equation of state. An exact reduction in terms of a scalar acoustic potential and the (scalar) thermal displacement is achieved. Properties and simplifications of these model nonlinear acoustic equations for unidirectional flows are noted. Specifically, the requirement that the governing system of equations for such flows remain hyperbolic is shown to lead to restrictions on the physical parameters and/or applicability of the model. A weakly nonlinear model is proposed on the basis of neglecting only terms proportional to the square of the Mach number in the governing equations, without any further approximation or modification of the nonlinear terms. Shock formation via acceleration wave blowup is studied numerically in a one-dimensional context using a high-resolution Godunov-type finite-volume scheme, thereby verifying prior analytical results on the blowup time and contrasting these results with the corresponding ones for classical (Euler) fluids.


Introduction.
Recently, there has been significant interest in the mathematics of nonlinear acoustics [42] and, specifically, in proving abstract mathematical results for model nonlinear acoustic equations, including well-posedness and control [43,44,45,46,9,10,11]. To this end, it is important to examine the model equations and understand their origin and applicability, in order for the mathematical results to have clear implications for a diverse set of physical contexts [12]. Broadly speaking, model equations of nonlinear acoustics [22,32,35] are derived with the goal of including only those physical effects important to a given application, while discarding physical effects deemed secondary or unimportant.
In the present work, we focus on a class of models for nonlinear acoustics in lossless fluids that conduct heat, building on a general theory [26,27] based on abstract considerations [28,29,30] proposed by Green and Naghdi (GN) two decades ago. 1 Since Green and Naghdi's initial work, there has been interest in understanding both acoustic ("first sound") [40,62,3] and thermal ("second sound") [3,4,70,63,71,2] nonlinear wave propagation under such nonclassical continuum theories, which motivates the present study.
Specifically, in Section 2 of the present work, we consider the governing equations of motion of a lossless compressible Green-Naghdi (GN) fluid that conducts heat. The equations of motion are made dimensionless in Section 3. Then, in Section 4, the discussion is specialized to two specific classes of barotropic fluids (i.e., ones for which the thermodynamic pressure is a function of the density alone): (i) perfect gases and (ii) liquids obeying a quadratic equation of state. In Section 5, an exact reduction to a coupled system of model nonlinear acoustic equations is achieved, using a scalar acoustic potential and the (scalar) thermal displacement, and these equations' properties and simplifications are also noted. Additionally, in Section 5.3, weakly nonlinear acoustic equations for lossless GN fluids are proposed on the basis of neglecting only terms proportional to the square of the Mach number in the governing equations, without any further approximation or modification of the nonlinear terms. Then, in Section 6, shock formation via acceleration wave blow up is studied numerically in a one-dimensional context using a high-resolution Godunov-type finite-volume scheme. In particular, prior exact results on the blow up time are verified and contrasted with the corresponding ones for classical (Euler) fluids. Finally, Section 7 summarizes our results and notes potential future work.
2. Governing equations. A lossless compressible GN fluid [26] in homentropic flow [Eq. (1d)] is described by the Euler equations [Eqs. (1a) and (1b)] of inviscid compressible hydrodynamics [73,67], which are augmented by a term involving the gradient of the thermal displacement from GN theory [40] and coupled to an energy equation involving the latter [Eq. (1c)]: ∇ · (̺δ) = 0, (1c) Here, v(x, t) is the velocity vector, ̺(x, t)(> 0) is the mass density, ℘(x, t)(> 0) is the thermodynamic pressure, η(x, t) is the specific entropy, δ(x, t) is the gradient of the thermal displacement from GN theory, 2 m is the (constant) "GN parameter" [40], which carries units of (L/T) 4 /K 2 with L and T denoting length and time, respectively, all body forces have been omitted, and t subscripts denote partial differentiation with respect to time t. To reiterate: the contributions of GN theory manifest themselves in the form of a (new) term on the right-hand side of the momentum equation (1b) and the additional solenoidality condition (1c). We assume that, initially, the fluid is in its equilibrium state, i.e., ̺ = ̺ 0 , v = 0, ℘ = ℘ 0 , δ = δ 0 and η = η 0 at time t = 0, where all the quantities demarcated with the subscript "0" are constant. Finally, we henceforth restrict to longitudinal (i.e., acoustic) waves, which implies irrotational flows, namely ∇ × v = 0 at the initial instant of time t = 0 and, hence, for all t > 0 [73]. The system (1) is not closed until ℘ is specified in terms of the other thermodynamic variables through an equation of state. In Section 4, we introduce the two classes of barotropic (meaning ℘ = ℘(̺) only [73,59]) equations of state that we restrict our discussion to.

Nondimensionalization.
Let us now rewrite the dimensional variables using the following dimensionless ones denoted by a ⋆ superscript: where the positive constants v m := max x,t |v|, L and δ 0 := |δ 0 | ≡ √ δ 0 · δ 0 denote a characteristic speed, length and thermal displacement gradient, respectively, and c 0 is the sound speed in the undisturbed fluid. 3 In terms of the dimensionless variables in Eq. (5), and recalling that we have assumed an irrotational flow, the governing system of equations (1) becomes where we have left the ⋆ superscripts understood without fear of confusion. Here, we have introduced two dimensionless groups: ǫ := v m /c 0 is the Mach number, and ς 2 := (mδ 2 0 )/(c 0 v m ) could be termed the "GN number" (squared), a quantity measuring the relative effects of the thermal displacement field to the acoustic field.
where, as before, we can formally introduce the effective (dimensionless) pressurẽ Here, it is evident that the (gradient of the) thermal displacement results in a O(ǫ) or "smaller" correction to the pressure ℘, even if ς = O(1).
The linearization of the system of equations (7) about some state where ℘ ′ ≡ d℘/d̺. The eigenvalues λ of the coefficient matrix in Eq. (9) satisfy For the system (7) to remain strictly hyperbolic (see, e.g., [51, Chap. 2]), we must require that the eigenvalues of the coefficient matrix of the linearized system (9) remain real and nonzero or, equivalently, that the discriminant of the polynomial (10) be positive: which simplifies to In the limit of ς → 0 + , the latter reduces to the usual assumption, i.e., ℘ ′ (̺ i ) > 0, under which the Euler system remains strictly hyperbolic. This assumption is, of course, equivalent to requiring that the speed of sound in the fluid remains a real number [see also Eq. (18) below and the discussion thereof]. More importantly, however, the hyperbolicity condition given in Eq. (12) can be rewritten in terms of the effective pressure℘ introduced in Eq. (8) to yield Therefore, under a redefinition of the functional form of the barotropic pressure, as in Eq. (8), the 1D lossless compressible GN fluid retains the same (hyperbolic) mathematical structure [49,76] as the 1D inviscid compressible Euler fluid! 4 4. Barotropic equation(s) of state. At this point in the discussion, we have not yet specified the equation of state, i.e., the relationship between the pressure ℘ and the remaining field variables (̺, v, δ, etc.), for the lossless compressible GN fluid under consideration. For the purposes of the present work, we only consider barotropic fluids, i.e., ones for which ℘ = ℘(̺) [73,59]. It is worth reminding the reader that the flow of a barotropic fluid must be barotropic, however, the converse does not necessarily have to be true.
In particular, following [17,18], we consider two example of barotropic equations of state. The first is the case of a perfect gas, 5 for which, under the above assumption of homentropic flow, it can be shown that the equation of state assumes a polytropic form [67], i.e., ℘ is a power law of ̺ only. The second is the case of a barotropic fluid with a quadratic equation of state. For fluids in general, especially liquids, it is not always possible to determine the mathematical relationship between ℘ and the other thermodynamic variables [54], hence expanding ℘ to second order in ̺ (at constant entropy η) about its equilibrium value is an approximation that is often made [54,48].

4.1.
Homentropic flow of a perfect gas. As mentioned above, for the homentropic flow of a perfect gas, the equation of state assumes a (dimensionless) polytropic from: where we remind the reader that the constant γ(> 1) is the adiabatic index, and we have used the dimensionless variables from Eq. (5).

4.2.
Homentropic flow with a quadratic equation of state. As mentioned above, for fluids in general and liquids in particular, an equation of state can be obtained by expanding ℘ to second order in ̺ about its equilibrium value, namely, where B/A is known as the nonlinearity parameter [5,48] and is related to the coefficient of nonlinearity β via β = 1 + B/(2A) [48]. In Eq. (15), ℘ 0 is the (dimensionless) pressure in the equilibrium state, and the expansion (15) has been performed at constant entropy η, which is allowed because of our assumption of homentropic flow [recall Eq. (1d)]. In this work, we take Eq. (15) as exact, i.e., we consider a (hypothetical) fluid that satisfies it exactly for some value of B/A (or, equivalently, β). Note that ℘ 0 = 1 because we made the pressure dimensionless using ̺ 0 c 2 0 in Eq.  14) about ̺ = 1, it follows that β ∈ (1, 4/3] based on the restrictions on the value of the adiabatic index γ for perfect gases [35]. Therefore, even though the equation of state (15) maintains the hyperbolicity of the Euler system for perfect gases, for the (hypothetical) GN liquids, for which we have assumed this equation of state to be exact, we must impose the further restriction that β ∈ (1, 3/2]. 6 Values of β as high as ≈ 7 [54, p. 582] and ≈ 10 [48, Table 8.4] have been measured, hence this restriction is nontrivial.   It is clear from Fig. 1 that the quadratic approximation [Eq. (15)] of the polytropic law [Eq. (14)] is quite good for air at room temperature, across a range of densities. Of course, for ̺ ≫ 1, we expect the quadratic approximation to become worse; indeed, the discrepancy between the solid and dashed curves is becoming noticeable near ̺ = 3 in Fig. 1. The most prominent feature of Fig. 1, however, is the singular behavior of the effective GN fluid pressure near the vacuum state (̺ = 0). Whereas for an Euler fluid, both Eqs. (14) and (15) predict a finite value for ℘(0), namely 0 and ℘ 0 + β − 2, respectively, the effective pressure℘ given by Eq. (8) blows up, i.e.,℘(̺) → +∞ as ̺ → 0 + .
At first glance, the latter observation might suggest that the most prominent difference in the acoustic behavior of classical (Euler) and nonclassical (GN) lossless compressible fluids might be at low densities, and this could be a regime wherein they can be experimentally differentiated. [From Eq. (8), one can show that at higher densities the difference diminishes and℘ − ℘ = O(̺ −1 ) as ̺ → ∞.] However, the low density regime is problematic as, just from a visual inspection of the gray curves in Fig. 1, it is evident that℘ ′ (̺) changes sign for some ̺ ∈ (0, 1) (for the chosen values of γ, β, ǫ and ς). When we verify the condition (12), which is necessary to maintain the hyperbolicity of the 1D governing system of equations, we find that℘ ′ (̺) changes sign at ̺ • = ǫς 2 1/(γ+1) for a perfect gas and at the real root 7  In other words, this GN theory appears unsuitable at low gas densities. These restrictions are, of course, related to the more specific ones derived by Jordan and Straughan [40,Section 3] between the dimensional GN parameter m and the state of the GN gas ahead of an acceleration wave.

5.
Model equations of nonlinear acoustics in GN fluids. In this section, following [40, Section 2], we discuss the reduction of the governing systems of equations from Section 1 using a potential function for the acoustic field and the scalar thermal displacement, namely φ and α such that v = ∇φ and δ = ∇α, respectively. For the acoustic field, we are allowed to introduce a scalar potential since we have assumed that ∇ × v = 0, i.e., that the flow is irrotational. Meanwhile, δ is, by definition, the gradient of the (scalar) thermal displacement α [26,40].
Due to the assumption of a barotropic equation of state, it is possible to introduce the (dimensionless) specific enthalpy function (see also [54,32,58]) through the definition where̺ is a "dummy" integration variable and h 0 ≡ h(1) is the (dimensionless) specific enthalpy of the fluid in its equilibrium state. Using the chain rule of differentiation, one can then show that At this point in the derivation, we introduce the (dimensionless) local speed of sound (see, e.g., [73,48]): which is a thermodynamic variable. Note that the derivative in Eq. (18) is, in fact, a partial derivate taken at fixed η, however, due our assumption of homentropic flow [recall Eq. (1d), which implies η = const. everywhere], it becomes a total derivative.

IVAN C. CHRISTOV
Using Eqs. (17) and (18) together with v = ∇φ, the conservation of mass equation (6a) can be recast as Meanwhile, rewriting the conservation of momentum equation (6b) in terms of φ and α and employing Eq. (17) 2 , we obtain Integrating Eq. (20) over space and enforcing the equilibrium conditions at t = 0, a relation between h, φ and α emerges: which is an exact expression of the conservation of momentum for a lossless compressible GN fluid in terms of the scalars φ and α. Equation (22) is not closed because, at this point, we have not specified c in terms of φ and α. In addition, Eq. (22) is still coupled to Eq. (6c) through the thermal displacement gradient ∇α ≡ δ. In terms of α and h, Eq. (6c) can be rewritten as c 2 ∇ 2 α + ∇h · ∇α = 0 (23) or, upon introducing Eq. (21) into the latter, Next, we consider two cases in which c can be written explicitly in terms of φ: homentropic flow of (i) a perfect gas and (ii) a fluid with a quadratic equation of state.
5.1. Homentropic flow of a perfect gas. Referring to [17] for the details, it can be shown that, in this case, and

Homentropic flow with a quadratic equation of state.
Referring to [17] for the details, it can be shown that, in this case, where W 0 (·) is the principal branch of the Lambert W -function [21], a special function function with a surprising number of applications in the physical sciences (see, e.g., [38]). Finally, using Eq. (21) to eliminate h, Eq. (29) furnishes the closure relation Note that by obtaining the closed-form expression (30) which, using Eq. (21) to eliminate h, becomes Thus, by recalling that β = 1 + (γ − 1)/2 for a perfect gas, we see that Eqs. (32) and (27) are equivalent, and that the sound speed in a perfect gas and in a liquid with a qudratic barotropic equation of state coincide in the weakly nonlinear regime, i.e., for ǫ ≪ 1 (see also the discussion in [17]).

Weakly nonlinear model equations.
In this section, we derive consistent weakly nonlinear (also known as "finite amplitude") approximations by neglecting terms of O(ǫ 2 ), in the spirit of Blackstock [7], Lesser and Seebass [50] and Crighton [22], who pioneered similar approaches for the thermoviscous case. In [17], it was shown that the consistent weakly nonlinear approximation, which does not involve "unnecessary" further modifications of the nonlinear terms, results in a solution closest to the reference Euler solution for a model shock tube problem in the ǫ ≪ 1 regime. In the context of GN theory, other weakly nonlinear model acoustic equations have also been proposed [40] based on further (approximate) manipulations of the nonlinear terms. As noted in Section 5.2, the exact equations governing the homentropic flow of a perfect GN gas are identical to the approximate equations governing the homentropic flow of a barotropic GN liquid with a quadratic equation of state, in the weakly nonlinear (ǫ ≪ 1) regime. Thus, to derive the consistent (also known as "straightforward" [17]) weakly nonlinear approximation, we drop all terms explicitly of O(ǫ 2 ) in Eqs. (22), (24) and (32), arriving at Clearly, Eq. (33c) can be introduced into Eqs. (33a) and (33b) to yield a system of two coupled partial differential equations (PDEs) for the scalars φ and α. At this point in the discussion, we have not specified the order of magnitude of the GN number ς. 8 Jordan and Straughan [40] argued that the nonclassical effects must necessarily be small, therefore it may be appropriate to take ς = O(ǫ). In this regime, Eqs. (33) simplify further, since we are neglecting all terms of O(ǫ 2 ): Equations (34)  Since only the square of ς appears in Eqs. (33), one could consider the regime ς = O(ǫ 1/2 ), which still corresponds to a "small" GN number ς, however, now the GN number is (asymptotically) larger than the Mach number ǫ. In such a weakly nonlinear regime, Eqs. (33) simplify to This appears to be, in a sense, the "simplest" straightforward weakly nonlinear model for a lossless compressible GN fluid in which there is a direct coupling between the acoustic and thermal displacement fields.
6. Shock formation. In this section, we present a numerical study of 1D shock formation in the homentropic flow of a lossless compressible GN fluid. In particular, we are interested in shock formation via acceleration wave blow up (see, e.g., [39,20]). This type of common scenario, first considered in [74] and later extended in [77,24,53], is one of a general class of problems of discontinuity formation from smooth initial data [1,76,25,52] with applications throughout continuum mechanics [13], and even in social and biological systems [6,72]. In the acoustics context, shock formation via acceleration wave blow up is, in particular, relevant to the theory of shock tubes and resonators [66,16], of which organ pipes are one example [60]. For the numerical study presented in this section, we use the so-called MUSCL-Hancock scheme [75] to solve the hyperbolic system of conservation laws (7); for implementation details, the reader is referred to [20] and [17, Appendix A], while noting the typographical correction regarding [17, Appendix A] given in [18]. The spatial grid size used is ∆x = 7.5 × 10 −4 , while the temporal step size is chosen adaptively [20,17] to satisfy the Courant-Friedrichs-Lewy stability condition (see [75,51]) in every computational cell at every time step. The spatial interval on which the system is solved is chosen to be large enough so that no reflections occur from the downstream boundary. Now, let us define the dimensionless acoustic density (also known as the condensation) ρ ≡ ̺ − 1. Then, we consider the 1D system of conservation laws (7) subject to the initial conditions and the boundary condition where H(·) is the Heaviside unit step function. The two conditions in Eq. (36) reflect the fact that the medium is initially in its equilibrium state, while Eq. (37) describes a pulse of finite duration t w (equivalently, finite width) being introduced at the x = 0 boundary of the domain at time t = 0 + . Specifically, we consider a sinusoidal pulse f (t) = ǫ sin(πt), which is particularly relevant to acoustics (see, e.g., [8, Section III.B]). For the initial and boundary conditions given in Eqs. (36) and (37), although ̺ and v are initially continuous, ̺ t (and, through the coupled nature of ̺ and v, v t ) suffers a jump discontinuity across the surface x = 0. Let x = Σ(t) be the location of this jump discontinuity at later times with Σ ′ (t) being its velocity relative to the fluid [13]. Then, we define the jump of any variable, say F, across Σ as where we have made the expression from [40] dimensionless using the variables introduced in Eq. (5) and, for a medium ahead at rest, used the identity α + x = δ + = 1/̺ + = 1, which is valid for unidirectional flow [recall Eq. (3)]. Furthermore, it can be shown that, if one of the first derivatives of a field variable suffers a jump discontinuity, then so do the rest, and all the jumps can be related to A(t) [40,Eqs. (3.9) and (3.19)]. Specifically, noting that V in [40] denotes Σ ′ , which we make dimensionless via V = c 0 V ⋆ , and recalling that v + = 0 and ̺ + = 1 for a wavefront advancing into a medium at rest, [40,Eqs. (3.9)] become Note the appearance of Σ ′ in the latter relations, which is not the case for an acceleration wave advancing into an Euler fluid at rest (see, e.g., [20,Eq. (18) Finally, from Eqs. (39)-(41), we find an expression for the slope of ̺(x, t) at the wavefront, i.e., at x = Σ(t), Since, for the problem posed above, the acceleration wave is compressive [i.e., A(0) > 0] [13], we expect acceleration wave blow up. In this case, it is easy to see that Eq. (39) implies a blow up time, i.e., t ∞ such that lim t→t − ∞ A(t) = +∞, of The expressions in Eqs. (39), (41), (42) and (43) are also valid for the Euler fluid by simply taking the limit ς → 0 + . Specifically, lim ς→0 + [[̺ x ]](t) = −ǫπ/(1 − ǫa 0 πt), lim ς→0 + Σ(t) = t and lim ς→0 + a 0 = β with β = (γ + 1)/2 for a perfect gas, which are the corresponding expressions for an acceleration wave advancing into a perfect (Euler) gas at rest [39,20,17]. Moreover, for unidirectional flow, owing to Eq. (3) and [40,Eqs. (3.9) and (3.19)], we can immediately conclude that an acceleration wave necessarily implies a temperature-rate wave [31,14,57,71], as was also shown in [40]. For the purposes of illustrating shock formation via acceleration wave blow up, let us take the medium to have similar properties to air at 20 • C, i.e., γ = 1.4 (⇒ β = 1.2), while ǫ = 0.26503 and t w = 1. Then, for these parameter values, the predicted blow up time of the acceleration wave's amplitude is t ∞ ≈ 1.00 for ς = 0 (i.e., the Euler fluid). The shock formation process is illustrated in Fig. 2 for (a) ς = 0 and (b) ς = 0.5, which corresponds to the scaling ς = O(ǫ 1/2 ) for ǫ = 0.26503. We observe from the dashed lines in Fig. 2 that the theoretical predictions (from singular surface theory) for the locations of and slopes at the wavefront agree very well with the numerical simulations of both the classical (Euler) and nonclassical (GN) fluids for all t ∈ [0, t ∞ ]. In particular, we see that shock formation occurs earlier and the wavefront moves slower in the GN fluid than in the Euler fluid [as predicted by Eqs. (43) and (41), respectively]. The numerical simulation additionally provides the full wave profile, even behind the wavefront Σ, where singular surface theory is not applicable. We observe that the classical and nonclassical wave profiles are quite similar, consistent with the assumption that the nonclassical effects are expected to be "weak," i.e., GN theory is a small correction to the classical theory.
7. Conclusion. In this paper, motivated by the final remarks in [17, p. 491], we re-examined the model equations of nonlinear acoustics for a class of lossless but thermally conducting nonclassical compressible fluids described by Green-Naghdi (GN) theory [26,27], in which a thermal displacement variable α is introduced in addition to the usual field variables of the fluid (i.e., density ̺, velocity v and specific entropy η). Specifically, the homentropic flow of GN fluids obeying a barotropic equation of state was studied; two representative barotropic equations of state were considered: (i) the polytropic equation of state of a perfect gas and (ii) a quadratic approximation to a general barotropic equation state (actually applicable to both gases and liquids), which was taken as exact.
A reduction in terms of a scalar acoustic potential and the scalar thermal displacement was achieved, yielding two systems of coupled nonlinear PDEs: (I) Eqs. (22), (24) and (27), which is exact for perfect GN gases; and (II) Eqs. (22), (24) and (30), which is exact for GN fluid with a quadratic barotropic equation of state. Additionally, under the assumption of a small Mach number (ǫ ≪ 1), some consistent weakly nonlinear models were noted. If the GN number is asymptotically of the same order as the Mach number, i.e., ς = O(ǫ), then the thermal displacement's evolution does not affect the acoustic field, and a one-way coupled systems of PDEs (34) is obtained. Alternatively, if the GN number is taken to be small but asymptotically larger than the Mach number, e.g., ς = O(ǫ 1/2 ), then the acoustic and thermal fields become coupled, as evidenced by Eqs. (35).
In addition, the exact nonlinear acoustic equations were studied for unidirectional flows, and a significant simplification was achieved through the introduction of an effective pressure. In this way, the system of conservation laws for mass and momentum in 1D takes an identical form for Euler (classical) and GN (nonclassical) fluids, as shown by Eqs. (7) and (8). Using this formalism, bounds were derived for the validity of GN theory. Specifically, unlike the Euler case, it was shown in Section 4.3 that the system of conservation laws (7) loses hyperbolicity for ̺ ∈ (0, ̺ • ) for some ̺ • . Thus, GN theory appears to not be applicable at low fluid densities.
In the same 1D context of unidirectional flow, shock formation via acceleration wave blow up was examined numerically. A high-resolution Godunov-type (i.e., shock-capturing) finite-volume scheme was used to solve Eqs. (7) numerically for the canonical problem of the injection of a sinusoidal density signal into a quiescent fluid. Jordan and Straughan's exact results [40] (based on singular surface theory) for acceleration wave blow up in this context were thus confirmed numerically, specifically showing that blow up occurs earlier in the nonclassical fluid and the wavefront's velocity is slowed down in comparison to the classical fluid. Additionally, the numerical simulations allowed for a comparison of the full wave profiles (at and behind the wavefront), showing that the wave profiles of the Euler and GN fluids are qualitatively similar for this model problem.