The effects of coupling on finite-amplitude acoustic traveling waves in thermoviscous gases: Blackstock's models

We consider the propagation of acoustic and thermal waves in classical perfect gases under a coupled, weakly-nonlinear system first derived by Blackstock. Our primary aim is to ascertain the usefulness of Blackstock's system as an approximate model of nonlinear acoustic phenomena. Working in the context of the piston problem, and using a solvable special case of the Navier--Stokes--Fourier system as our benchmark, we compare Blackstock's system against a simpler weakly-nonlinear model whose constitute equations are not coupled. In particular, traveling wave solutions (TWS)s are determined, the structure of the solution profiles is analyzed, numerical comparisons are presented, and follow-on studies are suggested.

1. Introduction. Blackstock [6], in a 1963 report written for the General Dynamics Corp., sought simplifications of the Navier-Stokes-Fourier (NSF) system by assuming that disturbances applied to the fluid cause the field variables to undergo small, but finite-amplitude, fluctuations about their equilibrium state values.
One of the many results presented by Blackstock [6] is a coupled, weakly-nonlinear, system that he derived as a model of (finite-amplitude) thermal and mechanical disturbances in thermoviscous perfect gases. Here, "thermoviscous" refers to fluids that exhibit not only (Newtonian) viscosity but also the ability to conduct heat, while a "perfect gas" [30] is one that obeys the equation of state (EoS) under the assumption that c p and c v are constants. In Eq. (1), ℘(> 0) is the thermodynamic pressure, (> 0) is the mass density, ϑ(> 0) is the absolute temperature, and c p > c v > 0 denote the specific heats at constant pressure and volume, respectively. Hereafter restricting our attention to acoustic propagation (i.e., irrotational flow), and employing a slightly altered version of his notation, we express Blackstock's [6, Eqs. (5c), (6)] as respectively, where the irrotational assumption means that ∇ × u = 0. Here, u = (u, υ, w) is the velocity vector; φ denotes the scalar velocity potential, where u = ∇φ since the flow is irrotational; θ = (ϑ − ϑ 0 )/ϑ 0 ; the adiabatic exponent is defined as γ = c p /c v , where γ ∈ (1, 5/3] under the perfect gas assumption [30]; b = 4/3 + µ B /µ is the viscosity number, where µ(> 0) and µ B (≥ 0) are the coefficients of shear and bulk viscosity, respectively; ν = µ/ 0 is the kinematic viscosity coefficient; κ = K/( 0 c p ) is the thermal diffusivity, where K(> 0) is the coefficient of thermal conductivity; c 0 , the speed of sound in the undisturbed gas, is given by c 0 = γ℘ 0 / 0 ; and a "0" subscript attached to a quantity denotes the (constant) equilibrium state value of said quantity. While Blackstock's [6] report was published over 50 years ago, the present author is not aware of any works wherein the utility of Sys. (2) as an approximation to the compressible NSF system has been investigated. It is therefore the aim of the present article to probe this topic in-depth; specifically, in the context of the classical piston problem, to which analytical solutions are sought in the form of traveling waves. Our primary focus shall be on determining the structure and the accuracy of the resulting solution profiles; the latter we accomplish by comparing the profiles of Sys. (2) with those of a simpler (weakly-nonlinear) model system (see Sect. 2 below), whose constitute equations are not coupled, using a solvable special case of the NSF system as our benchmark. Accordingly, we hereafter limit our focus to 1D planar propagation along, and perpendicularly to, the x-axis, i.e., to a flow geometry in which Sys. (2) reduces to where, in Sys. (3), u = (u(x, t), 0, 0), meaning of course that u(x, t) = φ x (x, t), and θ = θ(x, t).

A system based on the Blackstock-Lesser-Seebass-Crighton equation.
After obtaining Sys. (2), Blackstock goes on to derive, among other things, the following third order equation of motion (EoM): where the coefficient δ, termed the diffusivity of sound by Lighthill [18, §3.1], is given by and where P r = ν/κ denotes the Prandtl number. Here, we note that this PDE, which will play an important role in the analysis carried out below, appears as Blackstock's [6,Eq. (13)], wherein a 1 = γ − 1 in the case of a perfect gas, and is referred to by some authors as the Blackstock-Lesser-Seebass-Crighton (BLSC) equation; see, e.g., Ref. [13].
It should be noted that under the finite-amplitude scheme, Eq. (4) can also be derived directly from the following version of the (compressible) NSF system: which is the continuity equation, the (irrotational) momentum equation, the (recast) EoS for a perfect gas, which is readily derived from Eq. (1) using the Gibbs equation [30]. Here, s = ( − 0 )/ 0 is known as the condensation; we have set e := (η − η 0 )/c v , where η denotes the specific entropy; a superposed dot denotes the material derivative; Φ(≥ 0) is the dissipation function [30]; and the absence of both external body forces and internal heat sources has been assumed. Now introducing the following (additional) dimensionless variables: where L 0 (> 0) represents a characteristic (macro-scale) length, we recast Sys. (3) in dimensionless form, viz.: Likewise, the dimensionless, uncoupled, system formed by the 1D versions of Eqs. (4) and (7) reads (11b) In this communication, Pr = bP r is the longitudinal Prandtl number 1 , we have set κ := (Re Pr) −1 , and all diamond ( ) superscripts have been suppressed but should remain understood. 4. Traveling wave analysis. We begin this section with the following observation: since the equations of Sys. (10), as well as Eq. (11a), are all invariant under the transformation x → −x, we shall, without loss of generality, hereafter confine our attention to only right-running waves. That is, waveforms that travel with constant speed a > 0 and whose (single) argument, ξ, which some authors refer to as the wave variable, is of the specific form ξ := x − at.
To distinguish between expressions/quantities relating to Sys. (10) and those relating to Sys. (11), we employ, beginning in Sect. 4.1, the index n, where a subscript n = 1, 2 corresponds to the former and latter systems, respectively.
And, for later reference, we note that all shock thickness expressions given below are based on the following definition 2 : Eq. (12), we observe, applies only to those traveling wave profiles F(ξ) vs. ξ that take the form of a kink [2].

4.1.
Ansatzes and associated ODEs. We now introduce the ansatzes φ(x, t) = Π 1 (ξ) and θ(x, t) = θ p,1 T 1 (ξ). Here, θ p,1 = (ϑ p,1 − ϑ 0 )/ϑ 0 , where ϑ p,1 denotes the absolute temperature of the gas in contact with the face of the piston under Sys. (10). On making these substitutions and then integrating each of the resulting ODEs once with respect to ξ, Sys. (10) is reduced to where a prime denotes d/dξ, we have set U 1 = Π 1 , and A and B are the constants of integration. On imposing and enforcing the asymptotic conditions of the piston problem, i.e., T 1 , U 1 → 1 as ξ → −∞ and T 1 , U 1 → 0 as ξ → ∞, the following results are readily established: A = B = 0, which follows directly from the right-asymptotic condition; which we note is the positive root of a 2 − 1 = βa; and from which it is easy to see that θ p,1 > 0 (i.e., ϑ p,1 > ϑ 0 ). Here, for convenience, we have introduced β, the coefficient of nonlinearity [5,20], which in the case of a perfect gas is given by β = (γ + 1)/2. And in the case of Eq. (11a), applying the traveling wave assumption yields the simple first order ODE which we observe is of the Bernoulli type.

4.2.
Velocity field results. On eliminating T between the ODEs in Sys. (13), we are led to consider the single, second order, ODE [6], however, we neglect all terms (i.e., both nonlinear and linear) of O( 2 ); this, then, reduces Eq. (17) to the first order ODE which, like Eq. (16), is a simple Bernoulli equation.
And since, as in the case of Sys. (13), a is given by Eq. (14), Eq. (16) can be simplified to where we recall that a is the positive root of a 2 − 1 = βa.
On integrating Eqs. (18) and (19) via separation of variables, it is not difficult to establish that the integral curves of these ODEs assume the form of a Taylor shock, specifically, (20), we have imposed and enforced the wavefront conditions U 1,2 (0) = 1/2 and 1,2 (> 0), the shock thicknesses corresponding to the stated velocity field expressions, are given by Remark 1. In recent years some authors have begun to re-examine Blackstock's [6] models under a modified version of his approximation scheme in which only nonlinear terms of O( 2 ) are neglected; see Kaltenbacher [14] and the references cited therein. Under this less restrictive scheme, Eq. (17) remains second order but reduces to an ODE which is identical in form to the associated ODE of the Fisher-KPP equation. Here, we have introduced the re-scaling ξ = ζ γκ/( βaRe) and K 1 (> 0) is given by As shown in a number of works 3 , if the condition K 1 = 5/ √ 6 can be satisfied, then Eq. (22) admits the exact, kink-type, solution where r 0 = 1 − √ 2 for the wavefront condition U 1 (0) = 1/2. Unfortunately, K 1 = 5/ where we have set m := a/(γκ), we recall that A = 0, and we let C denote the second constant of integration corresponding to the temperature field. On evaluating the two remaining integrals using the software package Mathematica (ver. 9), and then setting C = 0, so that the resulting integral curve satisfies the left-asymptotic condition, Eq. (25) assumes the more tractable form Here, F (h, i; j; ς) denotes the Gauss hypergeometric series [1, §15]; we have set λ := 1 + m 1 /4; and we observe that where ψ(ς) is the digamma function [1, §6].

5.
Exact TWS of the NSF system: Becker's assumption. It is not, generally speaking, possible to exactly solve the un-approximated, compressible version of the NSF system, even, as Lord Rayleigh [26] recognized in 1910, if one's analysis is limited to the 1D piston problem involving a perfect gas. In 1922, however, Becker [4] obtained an exact solution to this problem for the special case P r = 3/4 and µ B = 0, where we note that monatomic gases are the only gases for which µ B = 0 [25, p. 550]. In 1944 Roy [27] revisited Rayleigh's [26] study. Implicit 4 in Roy's analysis was the fact that Becker's assumption could be generalized to (the equivalent of) Pr = 1, with µ B = 0 reverting back to the usual restriction µ B ≥ 0. And although he does not refer to it as such, Roy [27] appears to be the first to point out that the total enthalpy exhibited by the (piston problem) solution corresponding to Pr = 1 is constant; see also Morduchow and Libby [22] 5 , who like Roy recognize this fact but do not use the term "enthalpy", and Hayes [10], who refers to Becker's solution as the "constant total enthalpy" special case.
In this section we follow the more recent study of Margolin et al. [21], which revisits and extends that of Morduchow and Libby [22], but wherein the problem formulation is one of traveling waves with the same asymptotic conditions as in the present communication. And while they are valid for such values of the Mach number, none of the equations, expressions, etc., presented in this section are based on the assumption 1.

Velocity field. Assuming traveling waves and setting
where α < U(ξ) < 1, it is readily established that the velocity field under the aforementioned generalization of Becker's assumption, which for convenience we shall hereafter refer to as simply "Becker's assumption", is described by the following special case of Abel's equation: where = 2γ a(γ + 1)Re and and where α ∈ (0, 1) can also be expressed as α = 1 − /a. Here, we observe that Eq.
6.1. Shock thickness comparisons: Sys. (10) vs. Sys. (11). In Fig. 1 we have plotted the velocity field shock thickness expressions given in Eq. (21), for both n = 1, 2, and Eq. (38). From these curves it is clear that, for all perfect gases, the n = 2 case of the former, which corresponds to Sys. (11), best approximates that of the NSF system under Becker's assumption. As such, it follows that the n = 2 case of Eq. (20) better approximates the (exact) kink solution profile described by Eq. (37). Similarly, in Fig. 2, we see that the temperature field shock thickness expression given in Eq. (31), which corresponds to Sys. (11), is a much closer approximation to Eq. (40), the temperature field shock thickness expression corresponding to the NSF system under Becker's assumption, than that given in Eq. (30) for the full range of perfect gases.
6.2. Temperature field profile. We begin this subsection by introducing the waveform metric (see Ref. [13]): i.e., D(T 1,2 ) denotes the distance, with respect to the L 2 norm, between the (exact) temperature profile given in Eq. (39) and those given in Eqs. (26) and (28), respectively. From Fig. 3, wherein we have plotted D(T 1,2 ) vs. γ, it is obvious that for all γ ∈ (1, 5/3]. Thus, from this inequality, the plots presented in Fig. 2, and the fact that θ p,2 = θ p (see Sect. 5.2), it is clear that the temperature profile of Sys. (11) (i.e., T 2 (ξ)) provides, vis-à-vis that of Sys. (10), the closest approximation to the temperature profile of the NSF system under the assumption Pr = 1. We note in passing that Fig. 3 also indicates that the approximation offered by T 2 (ξ) improves, albeit very slightly, as γ is decreased from γ = 5/3, its maximum value for perfect gases. In particular, T 2 (ξ) is a better approximation for polyatomic gases at high temperatures, for which γ ≈ 1.1 [19, p. 8], than it is for diatomic gases (e.g., N 2 , O 2 ), for which γ = 7/5 [30]; in contrast, T 2 (ξ) is least accurate in the case of monatomic gases (see Remark 2), for which of course γ = 5/3. 6.3. Entropy approximation. We begin this subsection by introducing the version of Eq. (8c) used in the derivation of Eq. (11a), specifically, the (1D) linearized approximation to the former expressed in dimensionless form: see, e.g., Chester [7, p. 46]. Now setting η(x, t) = c v S 2 (ξ), and recalling that T 2 (ξ) = θ(x, t)/θ p,2 represents the temperature field TWS of Sys. (11), we integrate Eq. (47) once with respect to ξ. Thereafter solving for the ensuing constant of integration by assigning to η and θ their equilibrium state values, we obtain While computing T 2 (ξ) is a relatively simple matter, the derivation of Eq. (11a) involves use of an additional approximation. Referring the reader to Ref. [7, pp. 46-47] for the details, it follows that Eq. (48) can be re-expressed as an expression which, it should be noted, is consistent with Blackstock's approximation scheme. Here, we have set ∆S 2 (ξ) := S 2 (ξ) − S 0 , where we recall that S 0 := η 0 /c v . Thus, on calculating U 2 (ξ) from the n = 2 case of Eq. (20) and simplifying, we find that a result which has been known since (at least) the late 1950s; see Landau and Lifshitz [15], who derive it in the case of Burgers' equation. Thus, from Eq. (50) it is immediately clear that Thus, which under Becker's assumption reduces to Although ∆S(ξ) tends to constant, but unequal, asymptotic limits it is not a kink; see the Gold curves plotted in Fig. 4. This means that ∆S 2 (ξ) is incapable of matching the left-asymptotic limit of the ∆S(ξ) vs. ξ profile, as the plots in Fig. 4 clearly illustrate. It should be noted, however, that while ∆S p > 0, and ξ * is strictly negative 7 , further reducing ( 1) improves the overall approximation offered by ∆S 2 (ξ); indeed, one can get a sense of this by observing that And finally, by comparing Figs. 4(a,b) it is clear that, unlike the velocity and temperature profiles, ∆S(ξ) is quite sensitive to changes in γ.
Remark 3. Like Eqs. (50) and (54), the fact that ∆S(ξ) also attains a global maximum within the shock transition region for values of Pr other than one has been known for many years; see, e.g., Serrin and Whang [28].
which of course is the KdV equation.

7.
Closure. While counter intuitive, and established only for the special case Pr = 1, the superiority of Sys. (11) over Sys. (10) as a model of nonlinear acoustic phenomena is clearly welcome news to theoreticians and numerical modelers alike. Of course, the TWSs admitted by these systems should now be compared against each other using those of other (analytically) solvable special cases of the NSF system (e.g., those given in Refs. [8, §5.2] and [11, §3]), as well as those of arbitrary cases of the NSF system which can only be determined numerically (see Refs. [8,Eq. (22)] and [21,Sys. (9), (10), (13)] 8 ), as benchmarks.
Finally, the fact that its solutions do not satisfy the requirements of causality has prompted some researchers to seek generalizations of the NSF system which are free of this defect. And while they are naturally more complicated than the NSF system, these generalized systems are necessarily hyperbolic. Consequently, under the finiteamplitude assumption some such "causality-capturing" formulations would appear to be reducible to quasilinear hyperbolic systems. Cases where this simplification is possible naturally carry considerable advantage since a variety of analytical and numerical tools can then be brought to bear on them; see, e.g., Refs. [17,29]. For more on this branch of thermoviscous compressible flow, whose roots can be traced back to Maxwell, see, e.g., Refs. [9,16,23,24] and [29, § §3.2,3.3], as well as those cited therein.    numerical results and figures presented herein were generated using the software package Mathematica (ver. 9).