On the structure of isothermal acoustic shocks under classical and artificial viscosity laws: selected case studies*

Assuming Newton’s law of cooling, the propagation and structure of isothermal acoustic shocks are studied under four different viscosity laws. Employing both analytical and numerical methods, 1D traveling wave solutions for the velocity and density fields are derived and analyzed. For each viscosity law considered, expressions for both the shock thickness and the asymmetry metric are determined. And, to ensure that isothermal flow is achievable, upper bounds on the associated Mach number values are derived/computed using the isothermal version of the energy equation.


Introduction
In a celebrated paper published in 1849, Stokes [1] appears to have been the first to consider the impact of shear viscosity on the 1D propagation of linear acoustic waves in gases.Two years later, Stokes [2] went on to examine the propagation of linear, time-harmonic, acoustic plane waves in a gas whose only loss mechanism is its ability to radiate heat to its surroundings, the process of which he modeled using Newton's law of cooling; see Appendix A. In 1910, Rayleigh [3, pp. 270-271] presented a partial analysis of the isothermal propagation of an infinite (1D) "wave of condensation" in a constant viscosity, but thermally non-conducting, gas, the process of radiation being invoked to maintain the flow's isothermal nature.Later, Lamb (see, e.g., Ref. [4, § 360]) and, in 1953, Truesdell [5, p. 687] generalized Stokes' [2] radiant propagation model to include the effects of both viscosity and thermal conduction.In 2004, LeVeque [6] employed the assumption of isothermal flow to investigate 1D "delta shocks" in lossless perfect gases; because he did not consider the energy equation, however, LeVeque's analysis, like Rayleigh's [3, pp. 270-271], must be considered as only partially complete.
The primary aim of the present study is to not only complete, but also extend Rayleigh's analysis, wherein only the equation of motion (EoM) for the density field was derived and integrated, by considering the effects of various viscosity laws on the propagation and structure of acoustic shocks in the setting of the isothermal piston problem.Specifically, we determine and analyze 1D traveling wave solutions (TWS)s, under the assumption that the temperature of the gas in question is held constant, based on the following four viscosity laws: (i) constant shear viscosity (i.e., the case considered in Ref. [3, pp. 270-271]), (ii) shear viscosity proportional to mass density, (iii) von Neumann-Richtmyer artificial viscosity, and (iv) Evans-Harlow-Longley artificial viscosity.These particular laws were selected primarily because all lead to model systems that are amiable to study by analytical means.
To provide a mechanism for achieving isothermal flow, we, like Rayleigh, invoke the process of radiation, which like Stokes [2] we model via Newton's law of cooling.To ensure that isothermal flow is not only mathematically but also physically possible, under each of the four cases considered, we also determine, based on the full form of the isothermal energy equation, the upper bound on the range of (piston) Mach number values corresponding to each case.
In the next section, we begin our investigation with a review of the thermodynamics of perfect gases and the formulation of our governing system of equations.
2 Formulation of mathematical model

Thermodynamical aspects of perfect gases
As defined by Thompson [8], a perfect gas is one in which p(> 0), the thermodynamic pressure, ρ(> 0), the mass density, and ϑ(> 0), the absolute temperature, obey the following special case of the ideal gas law [8, § 2.5]: Here, c p > c v > 0 are the specific heats at constant pressure and constant volume, respectively, and γ = c p /c v , where γ ∈ (1, 5/3] in the case of perfect gases.Furthermore, a zero ("0") subscript attached to a field variable identifies the uniform equilibrium state value of that variable; i.e., in the present investigation, the gas is assumed to be homogeneous when in its equilibrium state.
With regard to the perfect gas assumption we note that the equilibrium state values of the adiabatic and isothermal sound speeds in the gas are given by and respectively, where we observe that b 0 ∈ (0, c 0 ).That is, c 0 and b 0 are, respectively, the speeds of infinitesimal-amplitude acoustic signals under the adiabatic and isothermal assumptions; see, e.g., Ref. [4, § 278], wherein b 0 is referred to as the "Newtonian velocity of sound".
In this study we will also make use of the thermodynamic axiom known as the Gibbs relation [8, p. 58], which in the case of a perfect gas can be written as where η is the specific entropy and E is the specific internal energy.Here, we note for later reference the relations where H is the specific enthalpy.

Navier-Stokes-Fourier system in 1D
For the 1D flow of a perfect gas along the x-axis of a Cartesian coordinate system, the Navier-Stokes-Fourier (NSF) system [9, p. 513] can, assuming the absence of all body forces and that the ratio µ b /µ is constant, be written in the following form: where p = p(x, t), ρ = ρ(x, t), ϑ = ϑ(x, t), and η = η(x, t) under this flow geometry.In Sys. ( 6), u = (u(x, t), 0, 0) and q = (q(x, t), 0, 0) are, respectively, the velocity and heat flux vectors; µ(> 0) is the shear (or dynamic) viscosity; 3 + µ b /µ is the viscosity number, where µ b (≥ 0) is the bulk viscosity; K(> 0) is the thermal conductivity; and r, which carries units of W/kg (i.e., m 2 /s 3 ), represents the external rate of supply of heat per unit mass.Moreover, Φ, the dissipation function [8,9], here takes the form Φ = µυ (∂u/∂x) 2 ; Eq. (6e) follows from integrating Eq. ( 4) after substituting in E from Eq. ( 5); and we note that the pressure gradient term that would normally have appeared in Eq. (6b) was eliminated using the relevant 1D special case of the thermodynamic relation [8, p. 71] In what follows we consider the compressive version of the classic (1D) piston problem.That is, the piston, whose face we take to be thermally insulated, is located at x = −∞ and moving to the right with constant speed u p (> 0) while the gas at x = +∞ is in its equilibrium state.In this regard we also observe that, since the piston's motion is strictly compressive, it follows that u ∈ [0, u p ] and ∂u/∂x ≤ 0 will always hold under the assumed flow geometry.
Additionally, we will invoke the assumption which of course is Stokes' hypothesis [8,9].Here, we stress that Eq. ( 9), which in the case of air, e.g., is only an approximation [8, Table 1.1], will not be a limitation on our analysis.This is because our numerical results will focus exclusively on the (monatomic) gas Ar -one for which Eq. ( 9) has been shown, by both theory and experiment, to hold exactly; see, e.g., Refs.[11,12].
In Eq. (13c), we have taken r to be given by the isothermal special case of Newton's law of cooling, viz.: where κ 0 (> 0), the "velocity of cooling" [2, p. 307], carries units of 1/s; see also Ref. [5, pp. 651, 687].Here, ϑ e = ϑ e (x, t), the temperature of the surrounding environment, is an additional dependent variable that is controllable by the experimenter based on "output" from Eq. (13c).This, of course, is necessary because the gas must be able to radiate away the resulting heat of its compression, at a rate that cannot be assumed constant, if our assumption of isothermal flow is to be satisfied + maintained.(Note that the ability of the gas to conduct heat plays no role under the isothermal assumption.)In this regard, we have adopted the assumptions stated by Rayleigh [10, p. 28] regarding the physical and thermal characteristics of the enclosing cylinder.

Viscosity laws: Classical and artificial
In this subsection we give the precise statement of the four viscosity laws mentioned in Section 1.Note that the first two stem from classical continuum theory while the latter two are of the artificial type.
(i) Constant shear viscosity: recall, this is the case considered in Ref. [3, pp. 270-271]; it is also the exact form of µ for the current (i.e., isothermal) problem under the kinetic theory of gases.
Here, µ 0 and ν 0 = µ 0 /ρ 0 represent, respectively, the (constant) equilibrium state values of the shear and kinematic viscosity coefficients; in the case of hard-sphere molecules we have, according to the kinetic theory of gases [8, § 2.7], an expression which also follows on setting "δ" in Ref. Lastly, since our investigation is to be carried out primarily by analytical methodologies, and we seek an approach that would allow one to compare/contrast these four cases in a consistent manner, we have taken ∆x = λ 0 , where ∆x is the spatial mesh increment used in the usual statements of both the vNR and EHL laws.And with regard to λ 0 we record here, for later reference, the expression which is easily obtained after eliminating c0 between Eqs. ( 19) and ( 20), and we note that 16/ √ 50π ≈ 1.277.
3 Traveling wave reduction

Ansatzs and wave variable
Invoking the traveling wave assumption, we set where ζ := x − vt is the wave (i.e., similarity) variable, and where the parameter v(> 0) will be seen to represent the speed of the resulting shocks.On substituting the above ansatzs into Sys.( 13), the latter is reduced to the following system of ODEs: where a prime denotes d/dζ.This system is to be integrated subject to the asymptotic conditions which of course correspond to a shock moving to the right; recall that v > 0.

Shock speed and associated ODE
Using the fact that Eq. (23a) integrates to allows us to, in turn, integrate Eq. (23b): where the resulting constants of integration were found to be K 1 = −ρ 0 v and K 2 = ρ 0 b 2 0 , respectively.Now eliminating g between the former and latter equations yields i.e., a Riccati type equation.Employing the asymptotic conditions once again leads us to consider a quadratic whose only positive root is where Ma = u p /b 0 is the piston Mach number.With these results in-hand, we can reduce Sys.(23) to the two-equation system Eq. (30a) is the associated ODE of our traveling wave analysis; in the next four sections, we shall integrate it, under each of the aforementioned cases of µ, subject to the wave-front condition f (0) = 1 2 u p .

Shock thicknesses, Q-metric, and jump amplitudes
Employing the notation we define the shock thickness of the velocity profile by a definition which Morduchow and Libby [18, p. 680] attribute to Prandtl, and that of the density profile by Here, ζ = ζ • j and ζ = ζ * j denote the relevant stationary points of f and g , respectively, i.e., f (ζ • j ) = 0 and g (ζ * j ) = 0; the subscript j represents the number [i.e., (i)-(iv)] of the case under consideration; and with regard to Eq. ( 33) we note the following: and We also define the Q-(or asymmetry) metric where respectively.As Schmidt [19, p. 369]2 points out, not only is Q a "sensitive measure of asymmetry", but it also complements the shock thickness as a characterizing metric since the latter "fails to give sufficiently detailed information about [shock] structure; . .." [19, p. 361].In Section 8.3 (below), we use Q to quantify the degree of asymmetry exhibited by each of the four g vs. ζ profiles studied below.Lastly, we define, following Morro [20] and Straughan [21], the amplitude of the jump in a function where it is assumed that both limits exist and that they are different.Here, we call attention to the fact that [[F]] is positive (resp.negative) when the jump in F is from higher (resp.lower) to lower (resp.higher) values.

Constant shear viscosity case
Recall that under Case (i), µ = µ 0 ; consequently, Eq. (30a) becomes which on setting This ODE, like the former, is a particularly simple special case of Abel's equation [24]; as such, it is easily integrated and yields the exact, but generally implicit3 , TWS: Here, the constant of integration is zero, by way of the fact that f (0) = u p /2 implies F(0) = 0; we have set and we note for later reference (in Appendix B) that Eq. ( 41) can also be expressed as Because k < −1, the integral curves described by Eq. (41) can only take the form of fully dispersed shocks, also referred to by some as kinks; see, e.g., Ref. [25, § 5.2.2].For this velocity traveling wave profile, therefore, we can show that the shock thickness is given by to which corresponds the stationary point here, we observe that which when related back to f yields In order to determine l i , we must first determine the (only) positive root of Π(Y ) = 0, where This quadratic arises when we attempt to solve g (ζ) = 0 after expressing it in terms of F and simplifying, i.e., when seeking the solution of which we do subject to the constraint |F| ∈ (0, 1).Denoting the aforementioned root by Y = F * i , it is not difficult to establish that which when related back to f yields where the stationary point of the corresponding g vs. ζ profile is given by With Eq. (49) in hand and observing that, under Case (i), Eq. (34) becomes the density profile is seen to admit the shock thickness Remark 2: Small-|ζ| and large-|ζ| approximations to F can easily be determined from the corresponding expressions for "U" given in Ref. [25,Remark 9], wherein "ξ" plays the role of ζ.Mention should also be made of the explicit, but approximate, result [25, Remark 10] for which the shock thickness is ˆ i = 16ν 0 /(3b 0 Ma), and we note that Ma 1 implies |k| 1.

The case µ ∝ ρ
Under Case (ii), µ = ν 0 × Eq. ( 25); as such, Eq. (30a) is reduced to the following Bernoulli-type equation: which is easily integrated and yields the exact TWS where the shock thickness for this case of f is given by With the aid of Eq. ( 56), the density profile for this case is found to be and this profile admits the shock thickness with corresponding stationary point From the latter it is readily established that ζ * ii < 0, which follows from the fact that Ma > 0, and, moreover, that which follows from the Taylor expansion of the ln-term about Ma = 0.

von Neumann-Richtmyer artificial viscosity
Under this formulation [i.e., Case (iii)], µ is given by Eq. ( 17); Eq. (30a), therefore, becomes A phase plane analysis of Eq. ( 63) reveals that its equilibrium solutions {0, u p } are rendered stable and unstable, respectively, and thus consistent with the compressive version of the piston problem, only when the "−" sign case is selected.
On rejecting the "+" sign case, it is a straightforward matter to show (see, e.g., Refs.[13,14]) that the velocity TWS under the vNR case is given by the following piecewise defined integral curve of Eq. (63): the shock thickness of which is [recall Eq. ( 32) We also find, on substituting Eq. (64) into Eq.( 26) and making use of Eq. ( 29), that the density TWS for this case is given by to which corresponds the (density) shock thickness [recall Eq. ( 33)] The stationary point corresponding to l iii is given by which we note can also be expressed as where we observe that − π 4 iii < ζ * iii < 0.
Remark 3: The f vs. ζ profile admits a pair of weak discontinuities4 , both of second order; i.e., the f vs. ζ profile exhibits two jumps, the amplitudes and locations of which are: Likewise, the present g vs. ζ profile also exhibits (two) weak discontinuities of order two; in the case of these jumps we have respectively.

Evans-Harlow-Longley artificial viscosity
Under this formulation [i.e., Case (iv)], for which µ is given by Eq. ( 18), Eq. (30a) reduces to which like Eq. ( 39) is a special case of Abel's equation.The integration of this ODE is not difficult; omitting the details, we find that Here, the corresponding shock thickness is given by which we computed by evaluating the limit Similarly, we have for the density profile which admits the shock thickness with Remark 4: As one differentiation of Eq. ( 73) reveals, the f vs. ζ profile under this case exhibits an acoustic acceleration wave6 ; in other words, the f vs. ζ profile suffers a jump discontinuity, the amplitude and location of which are: From Eq. ( 76) it is clear that the g vs. ζ profile also exhibits an acceleration wave, whose amplitude and location are 8 Numerical results

Parameter values: Ar
In the case of Ar at ϑ 0 = 300 K and p 0 = 50 mTorr, the gas/conditions on/under which Alsmeyer's [26] shock experiments were performed, we have the following: from which we find that b 0 ≈ 249.88 m/s, ν 0 ≈ 0.212 m 2 /s, λ 0 ≈ 1.084 mm.(82) Here, we have also made use of the result γ = 5/3, which according to the kinetic theory of gases holds for all monatomic gases; see Refs.[8,11,12].The values of ρ 0 , µ 0 , and c 0 were obtained from the NIST Chemistry WebBook, SRD 69 (see: https://webbook.nist.gov/chemistry/form-ser/);those of b 0 , ν 0 , and λ 0 , in contrast, were computed using Eq. ( 3), the defining relation ν 0 = µ 0 /ρ 0 , and Eq. ( 21), respectively.And so as to achieve iii = iv = 3λ 0 , which follows from Ref. [15, p. 233], based on our assumption ∆x = λ 0 and the shock thickness definition we have adopted, b iii = 1 4 √ 27 and b iv = 9/2 shall henceforth be taken., which displays data for non-isothermal shock propagation in Ar.In doing so, we have introduced the dimensionless reciprocal shock thickness parameter, viz.:  2], i.e., the fact that the former predicts a smaller shock thickness than the latter, is not unexpected.This is because under the isothermal assumption, K is absent from the NSF system [recall Sys. ( 13)]; i.e., the dissipation associated with the gas's ability to conduct heat, which tends to increase the shock's thickness, does not occur in isothermal propagation.

Determination of Q-metric values
In Table 1, values of Q based on the parameter values given above for Ar are presented.While those for Cases (ii)-(iv), as well as the M s = √ 2 case of Case (i) (see Appendix B), were computed using expressions obtained from evaluating the integrals in Eq. (36) analytically, the entries corresponding to the remaining M s values under Case (i) were computed from the numerical solution of Eq. (39), the process of which involved use of both the NDSolve[ ] and NIntegrate[ ] commands provided in Mathematica.
From Table 1 it is clear that, in all cases shown, Q is a strictly increasing function of the shock Mach number; i.e., profile asymmetry increases with M s .This behavior, we observe, is in agreement with the data presented in Refs.[19,26] for non-isothermal shock propagation in Ar.It is also clear from Table 1 that Cases (i)-(iii) all yield values of Q that are quite close to unity (i.e., quite close to being perfectly symmetric), and to each other, for Ma = 0.01.In this regard, we point out that Ma 1 defines the "realm" of finite-amplitude (or weakly-nonlinear) acoustics theory, within which velocity traveling wave profiles   36) have been changed to ∓ 1 4 π iii , respectively.‡ For this case, +∞ in Eq. ( 36) has been changed to iv ln (2).can be expected to have the form of Eqs. ( 54) and (56); see, e.g., Refs.[25,27], as well as those cited therein.
Finally, the last three entries in the first row of Table 1 highlight the following noteworthy feature associated with Case (i): for

Determination of Ma sup values
Implicit in Eq. (30b) is the thermodynamic requirement Θ e > 0. A simplified approach-which usually yields only approximate7 results, however-of ensuring that this inequality is everywhere satisfied begins with recasting Eq. (30b) as the inequality which is easily broken down into the individual cases where, in deriving these inequalities, we have made use of Eq. (32).It is noteworthy that, because f is a semi-compact8 function under our two artificial viscosity laws, both inequalities in Relation (87) become simply 0 < ϑ 0 on those intervals of the ζ-axis over which f (and therefore f ) is identically zero.
In the remainder of this subsection, we seek the value of Ma sup for each case appearing in Relations ( 86) and (87), assuming one exists, where Ma ≥ Ma sup means that isothermal propagation is not possible because the speed of the piston (i.e., u p ) is too great relative to b 0 .The numerical values given below were computed based on γ = 5/3 and the simplifying assumption M r = 1, where we define the radiative Mach number as M r := b −1 0 √ κ 0 ν 0 .

Case (i)
In dimensionless form, this case of Relation (86) can be recast as 0 < Π i (Ma), where and where we have set Ω i (M r ) = (8/3)M 2 r /(γ − 1).On solving Π i (Ma 0 sup ) = 0, which is readily accomplished due to its biquadratic form, we find that Ma sup ≈ Ma 0 sup , where here, , where we observe that Ω i (M r ) = 4 yields Y i = 2, and we have introduced the "dummy" variable Ma 0 sup .Unfortunately, Ma 0 sup is not, generally speaking, a very accurate approximation to Ma sup under this case.This is due to the fact that, in general, i.e., ζ e i = ζ • i , where Θ e (ζ e i ) = 0.One can, in principle, determine Ma sup for this case exactly by setting Θ e (ζ) = 0, after first using Eq.(39) to eliminate f from Eq. (30b) and then applying d/dζ to the result.The need to employ which is the only positive root of Π iii (Ma sup ) = 0.
Thus, under Case (iii): Remark 5: The Θ e (ζ) vs. ζ profile exhibits two jumps under Case (iii), the amplitudes and locations of which are: (99)

Case (iv)
Because the Case (iv) version of Eq. (30b) exhibits explicit dependence on not only f but also f , the corresponding case of Eq. (87) will only yield an approximation to Ma sup -one that can be shown to be quite poor.Hence, let us specialize Eq. (30b) to Case (iv) (see Section 7), with b iv = 9/2, and then replace Θ e (ζ) with its dimensionless counterpart T e (ζ); doing so yields, after simplifying, the exact expression where T e (ζ) ∝ Θ e (ζ) and Ω iv (M r ) = Ω iii (M r ).Although doing so is not simply a matter of setting T e (ζ) = 0, it can be shown that ζ e iv , the critical point of T e (ζ) under this case, is exactly given by where we have set Here, we observe that T e (ζ e iv ) = 0 only for Ma > 1, where the corresponding stationary point of the T e (ζ) vs. ζ profile is also its absolute minimum; for Ma . Consequently, we are led to consider the inequality 0 < Π iv (Ma), where Eq. ( 103) follows from Eq. (100) on setting Π iv (Ma) = T e (ζ e iv ), and as such contains no approximations.
It is now possible to determine Ma sup exactly by solving Π iv (Ma sup ) = 0, which can be simplified to read where we note the critical value and we observe that Ω iv (M r ) = 1/3.While it is a trivial matter to establish that suffice it to say that the process of determining Ma sup for the cubic case of Eq. ( 104) is algebraically intensive.As such, we set our analytical efforts aside and return to our numerical tools; hence, following an approach similar to that used in Section 8.4.1, we find that under Case (iv): • While f, g ∈ C ∞ (R) under both Cases (i) and (ii), the same is not true under Cases (iii) and (iv); specifically, f, g ∈ C 1 (R) and f, g ∈ C 0 (R), respectively, under the former and latter.
• By averaging the (1D) Euler system under finite-scale theory (FST), one obtains a system containing an "effective viscosity" that is given by Eq. ( 17) with b iii = 1/ √ 12 and λ 0 replaced by L, where L(> 0) is the length scale that appears in the averaging transform of FST; see Refs.[13,30] and those cited therein.
• While they are quite distinct from those presented for Ar in Refs.[19,26], the TWS profiles generated under Case (iv) are very similar to those given in Ref. [33,Fig. 2], which depict non-isothermal shocks in CO 2 ; see also Ref. [34].

Fig. ( 1 )
Fig. (1) has been generated in accordance with Ref. [26, Fig. 2], which displays data for non-isothermal shock propagation in Ar.In doing so, we have introduced the dimensionless reciprocal shock thickness parameter, viz.:

, 2 .
Interestingly, as shown in Appendix B, these values are three of the five shock Mach number values for which Case (i) yields explicit TWSs.

Table 1 :
Values of Q in Ar corresponding to Cases (i)-(iv) for selected Mach number values Case M s ≈ 1.005 For this case, ∓∞ in Eq. (