FINITE-AMPLITUDE ACOUSTICS UNDER THE CLASSICAL THEORY OF PARTICLE-LADEN FLOWS

. We consider acoustic propagation in a particle-laden fluid, specif- ically, a perfect gas, under a model system based on the theories of Mar-ble (1970) and Thompson (1972). Our primary aim is to understand, via analytical methods, the impact of the particle phase on the acoustic velocity field. Working under the finite-amplitude approximation, we investigate singu- lar surface and traveling wave phenomena, as admitted by both phases of the flow. We show, among other things, that the particle velocity field admits a singular surface one order higher than that of the gas phase, that the particle-to-gas density ratio plays a number of critical roles, and that traveling wave solutions are only possible for sufficiently small values of the Mach number.

1. Introduction. The basis of the present study is Marble's [27] theory of particleladen flows (PLF)s, which he put forth in 1970. Here, as in Ref. [27], the fluid phase is taken to be a perfect gas; i.e., one that obeys what is usually referred to as the "ideal gas law" subject to the proviso that c p > c v > 0 are both constant [34]. In Eq. (1), ℘(> 0) is the thermodynamic pressure, ϱ(> 0) is the mass density of the gas, ϑ(> 0) is the absolute temperature of the gas, and c p and c v denote the specific heats at constant pressure and volume, respectively, of the gas. Our focus here shall be on finite-amplitude acoustic phenomena in such "dusty gases", specifically, traveling waves and singular surfaces generated by compressive piston motion. Employing the tools and methods of classical analysis, our primary aim is to understand how the particle phase impacts the acoustic field, as modeled under the finite-amplitude approximation, vis-à-vis these two classes of waveforms. To simplify the analysis, however, we do not consider Marble's [27] theory in its most general form; instead, the governing system adopted in the present study stems from the (1D) special case of the former whose distinguishing assumption is i.e., the particle specific heat (at constant pressure) is negligibly small; see Thompson [34, pp. 553-556]. The immediate (practical) consequence of this assumption is τ T = 0 (⇒ the absence of thermal inertia), where τ T ∝ c pp is the thermal relaxation time.
It is appropriate to mention that the model formulated below is also a special case of the (1D) weakly-nonlinear dusty gas model developed by Crighton [15], who took the fluid phase to be a classical thermoviscous perfect gas, which we note is a class of (perfect) gases more general than that considered in the present study. In Ref. [15], however, the focus is limited to the gas phase and solutions are presented for only a single special case of the relaxation parameters, specifically (see Ref. [15, p. 75 where τ M is the momentum relaxation time (see Eq. (7) below). Hence, the present investigation will be seen to complement Crighton's analysis by providing solutions/results corresponding to the case wherein τ M > 0, τ T = 0 hold simultaneously, which of course is not a special case of Eq. (3), and for both phases of the flow. Here, we observe that our assumption c pp ≃ 0 corresponds to c ≃ 0 in Ref. [15].
The remainder of this exposition is organized as follows. In Sect. 2, what we term the Marble-Thompson model-1 1 is stated and the corresponding finite-amplitude system of equations (see Sys. (11), below) is presented. In Sect. 3, a singular surface analysis of the both the bi-directional-linear and right-running weakly-nonlinear versions this system is carried out. Then, in Sect. 4, we seek traveling wave solutions (TWS)s of the aforementioned system and examine their analytical properties. And lastly, in Sect. 5, some implications of our findings are reviewed and connections to other fields are noted.

Remark 1.
When the fluid phase involves a thermoviscous perfect gas, the thermal relaxation time can be expressed as where Pr is the Prandtl number and Nu is the Nusselt number. It is noteworthy that in his dusty gas model, Crighton assumes Nu = 2 from the outset; see Ref. [15, p. 71]. This "unforced" limitation is, however, easily circumvented through the use of Eq. (4). Hence, for the general case Nu > 0, the assumption stated in Eq. (3) implies that c pp assuming a thermoviscous perfect gas; see also the footnote on Ref. [34, p. 553].

The Marble-Thompson model-1.
Assuming that the fluid phase is a perfect gas of negligible thermal conductivity and that the particles are tiny, uniform, rigid spheres that are composed of a material for which Eq. (2) is a valid assumption, 1D acoustic propagation in such a dual-phase medium is approximately described by the following system: where the absence of both external body forces and internal heat sources has been assumed. Here, u = (u(x, t), 0, 0) and v = (υ(x, t), 0, 0) are the gas and particle velocity vectors, respectively; η = η(x, t) is the specific entropy of the gas; n = n(x, t) is the number of particles per unit volume; m p (> 0), the mass of each particle, is assumed constant; γ = c p /c v , where γ ∈ (1, 5/3] in the case of perfect gases; A = ϱ 0 c 2 0 denotes the adiabatic bulk modulus [31] of the clean gas, where c 0 = √ γ℘ 0 /ϱ 0 is the sound speed in the undisturbed (clean) gas; and τ M , which we recall is the momentum relaxation time, is given by where the density of the material that constitutes the dust particles, ϱ p , the shear viscosity of the (clean) gas, µ, and the particle radius, r p , are all positive and constant-valued. Also, in this communication the ambient state of both phases is assumed to be homogeneous-quiescent [31, p. 14]; i.e., while ℘ 0 , ϱ 0 , ϑ 0 , η 0 , and n 0 are constantvalued and positive, υ 0 = u 0 = 0, where a "0" subscript attached to a dependent variable denotes the ambient state value of that variable.
At this point the following observations regarding Sys. (6) are in order: (a) From Eq. (6b), the momentum equation for the gas phase, we see that the viscous contribution to the total stress tensor has been neglected; however, viscous effects still impact the flow, in the form of Stokes drag, via the parameter τ M . (b) Since it implies a flow free of the effects of thermal inertia (i.e., τ T ≃ 0), Eq. (2) also implies that where ϑ p (x, t) denotes the particle absolute temperature; see Ref. [34, p. 556]. (c) The right-hand side (RHS) of Eq. (6c) has been simplified using Eq. (8); see Remark 2 below. (d) Eq. (6d), the non-isentropic equation of state (EoS) for a perfect gas, stems, in part, from Eq. (1); see Ref. [20, p. 129]. And (e), our assumption of a homogeneous-quiescent ambient state is consistent with the requirement that the ambient state values of the field variables must, themselves, satisfy Sys. (6).
Observing that the 1D nature of the flow renders it irrotational, it follows that u(x, t) = ϕ x (x, t), where ϕ is the scalar velocity potential 2 corresponding to the gas phase. Also, anticipating our use of the finite-amplitude approximation, we now neglect both the RHS and the term uη x in Eq. (6c); consequently, the flow is now homentropic 3 and our EoS is reduced to a barotropic expression.

PEDRO M. JORDAN
Thus, on making the substitution u = ϕ x , observing that homentropic flow implies η(x, t) = η 0 , and introducing the following dimensionless variables: Sys. (6) is reduced to (10f) In Sys. (10), s = s(x, t) is the condensation; we term n = n(x, t) the particle condensation; U piston (> 0) and L(> 0), both constant, respectively denote the piston speed and a characteristic (macro-scale) length associated with the particular problem/flow under consideration; ϵ = c −1 0 U piston denotes the acoustic Mach number; τ = c 0 L −1 τ M is the dimensionless momentum relaxation time; κ, the mass fraction of particles, is given by κ = m p n 0 /ϱ 0 [34, Eq. (11.78)]; and all diamond (⋄) superscripts are here and henceforth suppressed, but should remain understood.
Under the finite-amplitude approximation, the basis of which are the assumptions ϵ ≪ 1 and |s| ∼ O(ϵ) (see, e.g., Jordan [20, §2.3]), Sys. (10) can be reduced to a weakly-nonlinear system involving only two PDEs. Omitting the details 4 , it is a straightforward matter to show that this approximating system readŝ whereĉ e =ĉ 0 / √ 1 + κ is the dimensionless version of the "equilibrium sound speed" [34, p. 557],ĉ 0 (= 1) is the dimensionless version of c 0 , and we observe thatĉ e < c 0 = 1. In deriving this system the following were also assumed: and, in accordance with the finite-amplitude scheme, only nonlinear terms of O(ϵ 2 ) have been neglected.

Remark 2.
It is the following, i.e., not Thompson's [34, Eq. (11.80)], that is the "exact" form of the entropy production equation vis-à-vis the assumptions of Ref. [34, §11.5]: which we have expressed in terms of the dimensional variables and parameters of the present communication. Here, we note that τ T (i.e., τ t in terms of Thompson's notation) is given, most generally, by where h(> 0) denotes the convective heat transfer coefficient. (In the thermoviscous case, this expression can, of course, be re-expressed as Eq. (4).) 3. Singular surface results: The piston problem in the half-space x > 0.

Mathematical preliminaries.
In this section our primary object is to investigate how the particle phase impacts shock propagation in the gas phase, in particular, the evolution of [[u]]. Employing the usual notation of singular surface theory, we define the amplitude of the jump discontinuity in a function F = F(x, t) across a propagating surface x = Σ(t) as where are assumed to exist, and where a "+" superscript corresponds to the region into which Σ is advancing while a "−" superscript corresponds to the region behind Σ. If F is a velocity component and represents the shock amplitude and x = Σ(t) the wave-front, or, in this case, the shock-front, across which F is discontinuous. The half-space version of the piston problem considered in this section has been investigated for over a century now; see, e.g., Lamb [25, §63] and Moran and Shen [28] for variants of this problem involving dust-free gases. Under the formulation considered, the piston, located at x = 0, impulsively begins "moving" to the right at time t = 0, thus compressing the dusty gas that fills the half-space x > 0; here, however, as in most of the previous studies of this version of the problem, the displacement of the piston in not explicitly taken into account 5 . Known in the acoustics literature as a signaling problem [16, p. 189], this version of the piston problem is, since it involves initial conditions (IC)s, also a start-up problem.
As Eq. (11a) is a quasilinear PDE, it does not directly lend itself to the study of wave phenomena in which [[u]] ̸ = 0, at least not by analytical methods. As such, we are left with two options. The first is to consider the linearized (i.e., ϵ := 0) version of Sys. (11), which we do in the next subsection. The second is to turn to what is generally referred to as the unidirectional approximation. As its name implies, this scheme is based on the assumption of unidirectional propagation; it is the methodology used to derive the KdV equation from the Boussinesq equations (see, e.g., Whitham [38, §8]) and the well known Burgers equation from the PDE which some authors refer to as the "Blackstock-Lesser-Seebass-Crighton" (BLSC) acoustic model (see, e.g., Crighton [14, p. 16]). While its disadvantage is clear (i.e., reflected waves are not allowed), the advantage of the unidirectional approximation is that it yields a simplified EoM whose highest-order time derivative is one order less than that of the original EoM.
As our attention in this section is limited to right-running waves, propagating along the +x-axis, with no possibility of reflection, reducing Eq. (11a) to a second order PDE via the unidirectional approximation is readily accomplished 6 . Thus, in terms of this new (simplified) EoM for the gas phase Sys. (11) becomes where we have made use of both the relation u = ϕ x and the fact thatĉ 0 = 1. Here, we have introduced β(> 1), the coefficient of nonlinearity [4], where β = (γ + 1)/2 in the case of perfect gases. It is noteworthy that Eq. (16a) is a special case of the PDE known today as the hyperbolic Burgers equation (HBE), which first arose as a kinematic wave-based traffic flow model; see, e.g., Refs. [12,19] and those cited therein.
Anticipating the singular surface analysis that shall be performed in Sect. 3.3, we close the present subsection by recasting Eq. (16a) as the kinematic-type system where, as in Ref. [19, §4.3], j has been introduced here to serve in the role of a flux. As its characteristics, which are defined by dx/dt = ±1, are real-valued and unequal, Sys. (17), and thus Eq. (16a), are actually strictly hyperbolic [26], meaning that their solutions satisfy the requirements of causality. Equally important, vis-àvis analyzing shock phenomena using singular surface theory, is the following: Since its eigenvalues are constants (i.e., ±1) and all nonlinearity is confined to its source vector, Sys. (17), and thus Eq. (16a), are also semilinear.

Remark 3.
Since we have assumed only right-running waveforms, it is clear that Σ(t) = t + x 0 under Sys. (17), where the constant x 0 denotes the initial location of Σ(t) on the x-axis. That is, in the case of Sys. (17) Σ(t) is a planar wave-front, across which u = u(x, t) may suffer a jump discontinuity, that propagates (to the right) along, and perpendicular to, the x-axis with speed dΣ(t)/dt = 1.
3.2. Velocity field shocks: Linearized system. In this subsection we consider the following initial-boundary value problem (IBVP) involving the linearized version of Sys. (11): Here, Θ(ζ) denotes the Heaviside unit step function; we have made use of both the relation u = ϕ x and the fact thatĉ 0 = 1; and we note that the first and fourth ICs imply that u + and υ + , the values of u and υ immediately ahead of Σ, are both zero.
On applying the Laplace transform [11] to Eqs. (18a), (18b), and the boundary conditions, IBVP (18) is reduced to, after making use of the ICs, the following system of subsidiary equations: which is to be solved subject tō Here, s is the Laplace transform parameter and a bar over a quantity denotes the image of that quantity in the Laplace transform domain. Solving Sys. (19) is not difficult; after simplifying, the exact (transform-domain) solutions can be expressed asū υ(x, s) =τ −1 where we have set σ := (τ c 2 e ) −1 . We now expand Eqs. (21) and (22) for large-s: Since the coefficient of sx in the factor e −sx , which each of these expansions contains, is −1, it follows 7 that dΣ(t)/dt = 1; therefore, in the case of IBVP (18) where the fact that the signal is inserted at the boundary x = 0 (⇒ Σ(0) = 0) means that the resulting constant of integration is zero. Now applying the theorem given in Ref. [7, §4] to Eq. (23), it is readily established that, across x = t, where we note that σ − 1/τ = κ/τ , and it follows that Eq. (25) is a shock-front in the case of the gas phase. In contrast, applying the aforementioned theorem to Eq. (24) reveals that, across see, e.g., Ref. [6, p. 182]. From Eqs. (27) we find that while υ(x, t) ∈ C 0 (R + , R + ), its first derivatives suffer jump discontinuities across x = t; i.e., with respect to the particle phase Eq. (25) describes an acceleration wave 8 in the case of IBVP (18). In closing this subsection we compute and plot the solutions of IBVP (18) by numerically inverting the Laplace domain solutions (i.e., Eqs. (21) and (22)) using the following formula: where M ≫ 1 is an integer and Eq. (29), we note, is a modified version of Tzou's [35] Riemann-sum inversion approximation, which Keiffer et al. [24] recently put forth, wherein the sinc(mπ/M ) factors have been introduced in order to reduce the Gibbs phenomenon that occurs when approximating our (discontinuous) solution u(x, t) using Tzou's formula. The sequences shown in Fig. 1 clearly illustrate the analytical results presented in this subsection, e.g., the fact that dΣ(t)/dt = 1. Comparing the (discontinuous) u vs. x profiles with those of the Type II case in Ref. [3, Fig. 3] reveals that the most obvious effect of the particle phase (i.e., takingτ > 0) on the acoustic velocity field is to "round-off" the corner that would otherwise exist on the leading side of the propagating Heaviside function; in contrast, comparison with those of the Type III case in Ref. [3, Fig. 3] makes clear that without the term u ttt , solutions of Eq. (18a) do not satisfy causality. Lastly, Fig. 1 also illustrates the fact that the acceleration wave amplitude (of the particle phase) is simply the slope of the tangent to the (continuous) υ vs. x profile at the acceleration wave(-front), which is located at the point (t, 0).

Remark 4.
Eq. (18a), which appears to have been first derived by Stokes [32] in his 1851 paper on acoustic propagation in radiating gases, is a special case of the 1D version of the PDE that has been termed the "Moore-Gibson-Thompson" (MGT) equation [23]; see also Ref. [22,Eq. (64)]. In this regard it is noteworthy that Eq. (18a) satisfies the exponential stability condition discussed by Kaltenbacher [22, p. 464]; specifically, her critical parameter (i.e., γ) is equal to κ, and κ > 0.
A phase plane analysis of Eq. (32) reveals the stability characteristics of its (two) equilibria, which we denote using a superposed tilde (i.e.,S). As α • > 0, since u + = 0, it follows thatS = {0, α • } are stable and unstable, respectively. Here, where we note for later reference the approximation where we observe that S(0) = 1. Eq. (35) indicates that, from the mathematical standpoint, the evolution of S(t) can occur in any one of the following three ways: Here, t ∞ , the time at which S(t) exhibits blow-up, is given by where we observe that t ∞ > 0 holds only under Case (III). From these results it is clear that Sys. (16) Turning our attention to the velocity field of the particle phase, we observe that taking jumps of Eq. (16b) yieldŝ On setting [[υ]] = 0, based on the findings of the previous subsection, it is easy to show that where we have once again made use of Eq. (28).

Remark 5.
Eq. (26) is recovered from Case (I) by first replacing α • with the first term of the approximation given in Eq. (34) and then letting ϵ → 0.

Traveling wave analysis:
The piston problem on the real line. We begin our study of this well known problem 10 with the following observation: since Eq. (11a) is invariant under the transformation x → −x, we shall, without loss of generality, hereafter restrict our attention to only right-running waves. That is, waveforms that travel with constant speed, a > 0, and whose single argument is the similarity variable ξ, which takes the specific form ξ := x − at.
The formulation of this version of the piston problem has the piston located at x = −∞, and moving to the right with constant speed U piston ; it also involves the assumption that at x = +∞ the medium is in its equilibrium state, i.e., at rest with all thermodynamic variables assuming their ambient state values.
And lastly, in this study we employ the following definition for the waveform metric known as the shock thickness: 4.1. Velocity field of gas phase. Introducing the ansatz ϕ(x, t) = (ξ) and substituting it into Eq. (11a), the latter is reduced to, after integrating once with respect to ξ, the second order ODÊ where a prime denotes d/dξ and K 1 is the constant of integration. On setting f(ξ) = ′ (ξ) and then imposing and enforcing the asymptotic conditions of the piston problem, i.e., f → 1, 0 as ξ → ∓∞, the following results are readily established: K 1 = 0, which follows directly from the right-asymptotic condition, and which we note is the positive root ofĉ −2 e a 2 − 1 = ϵβa; and thus, our associated ODE is found to beτ On integrating this (Bernoulli type) ODE, which we do subject to the usual condition f(0) = 1/2, we find that our TWS assumes the form of a Taylor shock 11 , specifically, where ϵ * := κ/β is a critical value of the Mach number (recall Eq. (37)), with shock thickness Here, we observe that this restriction on the Mach number must be satisfied if f(ξ) is to represent a physically admissible solution; in particular, f(ξ) does not satisfy the imposed asymptotic conditions, because it becomes an increasing function, when ϵ > ϵ * is taken.

Particle velocity field.
Under the traveling wave assumption, i.e., on setting v(ξ) = υ(x, t), Eq. (11b) becomes where in simplifying we have set Z = 4ξ/ℓ and defined As we expect v(Z) to be bounded for all Z ∈ R, we see from Eq. (46) that this implies lim respectively. 11 Note that a Taylor shock is a special case of a kink-type traveling wave profile.
With exp(−λZ) as the integrating factor, and employing the identity Eq. (46) yields the quadrature Here, K 3 , the constant of integration, has been set to zero to suppress this ODE's complementary solution, which is given by v c (Z) = K 3 exp(λZ). As shown in the Appendix, this integral can be evaluated for arbitrary values of λ(> 0). For the purposes of the present section, however, the following series representation of v(Z), which we obtained by repeatedly applying integration by parts to Eq. (50), shall prove adequate: From Eq. (51) it is clear that v is a perturbation of f. As such, satisfying the assumption given in Eq. (12) 1 requires that we now impose the restriction λ ≫ 1.
Note, however, that satisfying both λ ≫ 1 and 0 < ϵ < ϵ * (see Eq. (44)) requires us to impose the following (comprehensive) restriction: where That is, if the inequality in Eq. (52) is satisfied, then the velocity TWSs of both phases are physically admissible. From Eq. (51) we see that, as is true in the case of the gas phase, the particle velocity profile assumes the form of a kink, to which a shock thickness value can be assigned. Denoting by ℓ p the shock thickness corresponding to the particle velocity field, we now revert back to the (independent) variable ξ and introduce the function Noting that v ′′ 1 (ξ m ) = 0, where it is not difficult to establish, based on Eqs. (40) and (48), that From this it is clear that ℓ p < ℓ. Here, we observe that ξ m < 0; more specifically, And lastly, from Eq. (51) it is easily shown that v(0) = 1 2 from which we see that 0 < v(0) < f(0) = 1/2. For completeness we note that Eqs. (A.66) and (A.71) (both) give which is an exact result. Here, ψ(ζ) is the digamma function [1, §6] and we observe that for v(0) given by Eq. (59), where 1 − ln(2) ≈ 0.3069.

Observations and final remarks.
(i) For both singular surface and traveling wave phenomena, the model considered demands that we restrict the value of the Mach number; in both cases, the ratio ϵ * := κ/β appears as a factor in the corresponding upper bound. (ii) Observing that ϵ * can be expressed as where it should be noted that the numerator factors relate only to the particle phase while the denominator factors relate only to the gas phase, the condition ϵ < ϵ * , e.g., is seen to imply that i.e., we have placed an upper bound on a parameter over which an experimentalist has full and direct control. Such relations can, in principal, allow us to experimentally determine properties of the particle phase when the type of (perfect) gas that constitutes the gas phase is known. (iii) Under the present model, the particle velocity field admits a singular surface one order higher than that of the gas phase. (iv) On letting m p → 0, Eq. (11a) reduces to the 1D (dimensionless) case of the PDE some authors refer to as the lossless BLSC equation (see, e.g., Ref. [20, Eq. (51)]) while Eq. (11b) becomes simply υ = u; i.e., in this limit the particles move with the same velocity as the gas, but their vanishing mass does not allow them to influence the motion of the latter. (v) With regard to the gas phase, 1/τ plays the role of a Reynolds number, specifically, one based on Lighthill's diffusivity of sound, δ(> 0) [34, §4.13]. What is more, within this analog to thermoviscous gas dynamics we find that m p corresponds to δ, as highlighted by the following observation: Since m p → 0 (⇒τ → 0) causes ℓ → 0, the f(ξ) vs. ξ profile "shocks-up" in the same way that the (traveling wave) profile described by Ref. [16,Eq. (22)] does when we let δ → 0. (vi) The fact that ℓ p < ℓ indicates that, in the context of traveling waves, the particle velocity field experiences less dissipation than that of the gas phase; it should be noted, however, that ℓ p → ℓ (⇒ particle phase dissipation increases) as λ → ∞. Since v(Z) must be physically well-defined for all Z ∈ R, we now employ the analytic continuation 13 of Φ given by Guillera  where F (ς 1 , ς 2 ; ς 3 ; ζ) denotes the Gauss hypergeometric series [1, §15].