Poroacoustic Traveling Waves under the Rubin–Rosenau–Gottlieb Theory of Generalized Continua

: We investigate linear and nonlinear poroacoustic waveforms under the Rubin–Rosenau– Gottlieb (RRG) theory of generalized continua. Working in the context of the Cauchy problem, on both the real line and the case with periodic boundary conditions, exact and asymptotic expressions are obtained. Numerical simulations are also presented, von Neumann–Richtmyer “artiﬁcial” viscosity is used to derive an exact kink-type solution to the poroacoustic piston problem, and possible experimental tests of our ﬁndings are noted. The presentation concludes with a discussion of possible follow-on investigations.


Introduction
What is known today as the "RRG theory" was put forth by Rubin et al. [1] in 1995. This phenomenological-based theory of generalized continua is thought capable of modeling dispersive effects caused by the introduction of a medium's characteristic length, which Rubin et al. denote as α. Under RRG theory, α is regarded as an inherent material property. From the modeling standpoint, this theory exhibits a number of appealing features, the two most important of which are the following: (i) it is only the pressure stress (i.e., isotropic) part of the Cauchy (i.e., total) stress tensor and the specific Helmholtz free energy that are modified, but these modifications are achieved by adding perburtative terms, which must satisfy certain constraint equations, to the constitutive relations of the former and latter; and (ii), no additional boundary nor initial conditions, beyond those required to solve classically formulated problems, are needed ( [1], p. 4063).
To date, RRG theory has only been applied to single-phase media; see, e.g., Ref. [2] and those cited therein. Hence, there is an obvious need to investigate the nature of the solutions, e.g., those of the traveling wave type, predicted by this theory in the case of multi-phase media.
Accordingly, the aim of this communication is to carry out a preliminary investigation of RRG theory in the context of acoustic problems involving propagation in dual-phase (specifically, fluid + solid) media-dual-phase media being, of course, the simplest case of multi-phase media. Employing both analytical and numerical methodologies, we consider linear and finite-amplitude poroacoustic propagation under the RRG-based generalization of what some refer to as the Brinkman poroacoustic model (BPM) (Although he does not refer to it as such, the general, multi-D, version of the BPM follows on setting C = 0 in Burmeister [3].). Here, it should be noted that the original version of the drag law on which the BPM is based reads (see, e.g., Refs. [4,5]) ∇P =μχ∇ 2 u − (µχ/K)u. (1) In System (3), ℘(> 0) is the thermodynamic pressure; (> 0) is the mass density of the gas; n is the specific entropy of the gas; the parameter γ denotes the ratio of specific heats, where γ ∈ (1, 5/3] in the case of perfect gases; we take α(> 0), which carries the unit of length, to be a constant (That is, we have assumed the simplest version of RRG theory; see (Ref. [1], Equation (20)).); the problem geometry dictates that, here and henceforth, u = (u(x, t), 0, 0), ℘ = ℘(x, t), and = (x, t); and a zero ("0") subscript attached to a dependent variable denotes the (constant) equilibrium state value of that variable, where we note that u 0 = (0, 0, 0).
Here, we observe that since the flow has been assumed homentropic, our RRG-based poroacoustic model is obtained by perturbing only the pressure tensor term in the BPM. Also, we record for later reference that c 0 = γ℘ 0 / 0 is the (constant) equilibrium state value of the sound speed, i.e., the speed of sound in the undisturbed gas; see, e.g., (Ref. [7], §4.3).

Finite-Amplitude Equation of
Motion: The Case µ := const.
Hence, on invoking the finite-amplitude approximation, and introducing the following dimensionless variables: where the positive constants L and U p respectively denote a macro-length scale characteristic of the propagation domain and the magnitude of the peak particle velocity in the gas, it is not difficult to establish (See, e.g., the derivation performed in (Ref. [8], §2), and note that (Ref. [8], Equation (10)) is the σ, δ := 0 special case of Equation (5) where here and henceforth all diamond ( ) superscripts have been suppressed for convenience. In Equation (5), which we note reduces to the corresponding EoM for the (1D) BPM on setting a 0 := 0, = U p /c 0 is the Mach number, where 1 is assumed; δ = νχL/(c 0 K) is the dimensionless Darcy coefficient, where ν = µ/ 0 is the kinematic viscosity of the gas; a 0 , the dimensionless version of α, is given by a 0 = α √ 2/L; we have set σ := χ/Re B , where Re B = c 0 L/ν is a Reynolds number, and whereν =μ/ 0 ; and β(> 1) denotes the coefficient of nonlinearity [9], which in the case of a perfect gas is given by In deriving Equation (5) we have assumed that δ, σ, a 0 , |s| ∼ O( ) and, in accordance with the finite-amplitude approximation, only nonlinear terms O( 2 ) have been neglected. Although derived under the finite-amplitude approximation, Equation (5) is still too complicated for treatment by analytical means. Fortunately, however, the nature of the problems to be considered below is such that we may employ the uni-directional approximation to reduce the order of Equation (5) by one and confine its nonlinearity to a single (quadratic) term. Omitting the details, we find that under, say, the right-running case of this approximation (See, e.g., Crighton's ( [10], p. 16) derivation of the acoustic version of Burgers' equation.), which in the present setting reads φ x −φ t , our EoM becomes, after making use of the relation u(x, t) = φ x (x, t) and simplifying, which on switching to the variables X = x − t and T = t is further reduced to If we once again make use of the right-running approximation, which now takes the form u T −u X , to re-express only the third order term in Equation (8), which is justified since a 0 ∼ O( ) (i.e., (a 2 0 /2)u TXX is a "small" term), then Equation (8) assumes its final form, specifically, a PDE which we term the damped Burgers-KdV (dBKdV) equation.
In closing this sub-section we stress that Equations (7)-(9) apply only to right-running waveforms; i.e., to problems wherein reflection (to the left) is not possible.

Comparison of Linearized EoMs: The Cauchy Problem
In this section we compare the BPM with its RRG-based counterpart under the linear approximation, which at the EoM level corresponds to setting := 0. We do so in the context of what is perhaps the best known problem from classical PDE theory.
To this this end, we consider the linearized version of Equation (9) in the setting of the following initial value problem (IVP), i.e., in the setting of the classical Cauchy problem: Here, we take f (X), our initial condition (IC), to be defined on the real line and such that its Fourier transform exists.
On applying the Fourier transform to both Equation (10a) and the IC, and then solving the resulting (first order) ODE, it is readily shown that where k is the Fourier transform parameter and a hat over a quantity denotes the Fourier transform of that quantity. In turn, applying F −1 (·), the inverse Fourier transform, to Equation (11) gives 3.1. The RRG Case: a 0 > 0 Using the convolution theorem, and letting Ai(·) denote the Airy function of the first kind, the RRG (i.e., a 0 > 0) case of Equation (12) can be recast in the more explicit form For obvious reasons, the following two special cases of f (X) are of particular interest: Here, d(·) denotes the Dirac delta function and b(> 0) is a (dimensionless) constant.
3.2. The BPM Case: a 0 := 0 If we assume instead the BPM, then the solution of IVP (10) is readily obtained on setting a 0 := 0 in Equation (12); for the two aforementioned cases of f (X), we find that

Remarks: RRG vs. BPM
With regard to the Gaussian IC, the primary difference between the linearized RRG and BPM cases is that the pulse profile corresponding to the former instantly becomes oscillatory about the X-axis, due to the Airy function in its integrand, while that of the latter maintains, for all T > 0, the shape and strict positivity of the initiating Gaussian. The clearly contrasting behaviors exhibited by these two models should, therefore, allow researchers to experimentally determine which of the two best describes propagation in a given poroacoustic system. Before examining it in its most general form, and for the benefit of those readers who are not well acquainted with the intricacies of nonlinear evolution equations, it is instructive to first review selected special cases of Equation (9). The right-running models discussed in the next three sub-sections, all of which, it should be noted, have applications beyond poroacoustics, will each have a role to play in the analysis performed in Section 4.4.
Since we have assumed δ 1, applying the Kryloff-Bogoliubov asymptotic expansion method to the dKdV equation yields, as Ott and Sudan [11] have shown, the large-T expression where we have taken N(0) = 1 (see Ref. [11]), and we let see also Ref. [12], as well as Ref. [13] (Ref. [13] contains a number of recently identified typographical errors; see Appendix A below.) and those cited therein. Equation (17) represents a damped, and decelerating, solitary waveform, and as such it cannot be a soliton in the classical sense [14]. Note, however, that the acoustic version of the classic soliton solution of the KdV equation (see Ref. [8]) is recovered as the limiting case

Case (II): Damped Burgers' Equation
This case, which corresponds to setting a 0 := 0, reads Equation (20) is the right-running EoM stemming from the BPM, and in this context it has recently been investigated by Rossmanith and Puri [15].
As shown by Nimmo and Crighton [16], this generalization of Burgers' equation does not admit a linearizing (i.e., Cole-Hopf type) transform. As shown by Malfliet [17], however, its TWS, which assumes the form of a damped kink, is readily approximated. To the order expressed explicitly in Ref. [17], the TWS of Equation (20) is given by Here, where λ B = 4σ/( β) is the shock thickness exhibited by the TWS given below in Equation (25); In Ref. [17], the parameter c, which herein has the value c = 2/λ B , is defined so that Equation (21) yields the limiting case Equation (25) and λ B are the TWS, which we note takes the form of a Taylor shock, and corresponding shock thickness, which was determined using Equation (2), respectively, admitted by the classic Burgers equation.

Numerical Results
Inspired by, and closely following, Zabusky and Kruskal's [14] analysis of the classic KdV equation, in this subsection we perform numerical experiments on Equation (9), and its special cases listed above as Cases (I) and (II), in the setting of the following initial-boundary value problem (IBVP) with periodic boundary conditions: u(X, 0) = cos(2πX), |X| < 1.
In (Ref. [14], Figure 1), snapshots of the evolution of the KdV's solution profile were displayed in units of (dimensionless) time T B , where Zabusky and Kruskal used T B to denote the "breakdown time" (i.e., the time at which finite-time gradient catastrophe occurs) of the solution to the Cauchy problem involving the classic (i.e., undamped) Riemann equation. In our analysis of IBVP (30), T * B shall play the role of T B .
The graphs presented in Figures 1-3 were computed and plotted using MATHEMATICA (ver. 11.2). Except for the value of β(= 1.2), which corresponds to diatomic gases (e.g., air) [9], all other parameter values were selected based on our desire to produce clear, illustrative, graphs and the need to satisfy the assumptions under which Equation (9) was derived.
From Figure 1 it is easy to see that, except for attenuation of the profile (caused by the Darcy term) and a slight phase shift, the dKdV profiles are qualitatively similar to those of the classic KdV equation in the setting of IBVP (30). And, as is also true in the case of the latter, reducing (resp. increasing) a 0 increases (resp. decreases) the number of pulses seen in Figure 1b.
In contrast, the plots shown in Figures 2 and 3 highlight the fact that, like that of the damped Burgers equation, the dBKdV profile suffers attenuation, again due to the Darcy term, and it also develops a "dull sawtooth" appearance as it shocks-up (to the right), but never breaks since σ > 0. More interesting, however, is the fact that for large-T, both the damped Burgers equation and dBKdV profiles are seen to re-assume the periodic form of the IC. As Figures 2c and 3c illustrate, both profiles evolve to become damped, and in the dBKdV case slightly phase-shifted (to the left), versions of the IC given in Equation (30c). This suggests that for sufficiently large values of T, one may employ the approximations u(X, T) ≈ u 1,2 (X, T), where and where we require 1,2 (T) > 0. Comparing the blue-broken curve in Figure 2c with its counterpart in Figure 3c we see that 2 here, for simplicity, we have assumed 1,2 (T) and ψ 1,2 (T) to be linear functions of T. In the setting of IBVP (30), then, the presence of the third order (i.e., RRG) term in the dBKdV equation gives rise to both a phase shift and slightly less attenuation vis-à-vis the damped Burgers equation.
While their usefulness may be limited to certain "windows" of T-values, the functions 1,2 (T) and ψ 1,2 (T) should be constructible based on Equation (31) and numerically generated, large-T, data sets using one of the many data-fitting methodologies found in the literature.

The RRG Case with "Artificial" µ
In 1950, von Neumann and Richtmyer (vNR) [23] introduced their artificial viscosity coefficient. In this section, we make use of this celebrated artifice not to regularize numerical schemes used to calculate shock profiles, as was vNR's aim, but rather to obtain an analytical solution to the poroacoustic version of the piston problem (Unlike Ref. [23], wherein Lagrangian coordinates were used, in this communication we employ the Eulerian description; see, e.g., (Ref. [24], §V-D-1) wherein vNR's system is recast under the latter.).
To this end, we return to the RRG case of System (3) and assume that but continue to regardμ as a positive constant. Here, we have expressed the length-scale factor in vNR's artificial viscosity coefficient as α, instead of some grid spacing ∆x.
For simplicity, we now assume that the porous solid in question is comprised of packed beds of rigid solid spheres, all of radius r(> 0), which are fixed in place. For such a configuration, the permeability is given by the well known Kozeny-Carmen relation [4]: As these spheres are scatters of acoustic waves, we take α to be proportional to the characteristic length now associated with our dual-phase medium; i.e., we take α = b 1 r, where b 1 (> 0) is an "adjustable" (dimensionless) constant ( [24], p. 233).
If, moreover, we limit our focus to kink-type waveforms, as physical intuition suggests, and have the piston located at x = −∞ and moving to the right along the x-axis, then u x < 0 and Equation (32) becomes where we have used the relation u = φ x . Returning to our dimensionless variables, and once again applying the finite-amplitude approximation, it is readily established that, under the aforementioned assumptions, the following (simpler) weakly-nonlinear PDE replaces Equation (5) as our bi-directional EoM: In Equation (35), which we observe applies only to the RRG case, we have set a 1 := b 1 r √ 2/L, where we require that a 1 ∼ O( ), and where the requirement δ 1 ∈ (0, 1) implies that b 1 must satisfy the inequality Assuming the gas at x = +∞ is in its equilibrium state, and thus motionless, and observing that in the present context U p is the dimensional speed of the piston, we let φ(x, t) = G(η), where η = x − v 1 t and the (dimensionless) shock speed v 1 is taken to be a positive constant, and then substitute into Equation (35). On integrating once with respect to η and then imposing/enforcing the asymptotic conditions g → 1, 0 as η → ∓∞, respectively, Equation (35) is reduced to the ODE where we note that the resulting constant of integration is zero. In Equation (38), g(η) = G (η), where a prime denotes d/dη; we have defined recalling that β is the coefficient of nonlinearity (see Equation (6)); and which we observe is the positive root of To apply the solution methodology employed in (Ref. [2], §2) to Equation (38), the following condition must be satisfied: In (Ref. [2], §2), satisfying Equation (42) required that the value of the Mach number be fixed, a constraint which clearly limits the usefulness of the TWSs presented in that article. Here, however, we shall use this restriction to our advantage; specifically, in the following sense: Since the value ofμ for a given poroacoustic flow is, in general, not known, and we are seeking a kink-type TWS, then the only possible solution of Equation (42) in the present context is where we observe that v 1 = v 1 c 0 is the dimensional shock speed and, moreover, that v 1 > 1 implies v 1 > c 0 . On imposing the usual wave front condition g(0) = 1/2, but otherwise referring the reader to (Ref. [2], §2) for details regarding its derivation, the TWS we seek is where K = tanh −1 −1 + √ 2 . Letting λ 1 = 1 /L denote the dimensionless shock thickness (Recall Equation (2)) admitted by Equation (44), it is easily established that where 1 is the corresponding dimensional shock thickness. Also, with regard to computing λ 1 , it should be noted that g (η * ) = 0, where η * (< 0) is given by and where it should also be noted that g(η * ) = 5/9. The usefulness of Equation (45) might be ascertained as follows. Assume that v 1 and 1 can both be determined, either directly or indirectly, from experimental measurements and, moreover, that both are (at most) slowly varying functions of time. With v 1 known, b 1 can, of course, be computed using Equations (36), (39) and (40). If this (inferred) value of b 1 satisfies the inequality in Equation (37), a 1 ∼ O( ) is also satisfied, and the measured value of 1 is in agreement with that computed from Equation (45) over, say, some span of time t ∈ T , then we can expect Equation (45) to prove useful as an approximation within the transition region of our kink-type traveling wave profile for t ∈ T .

Discussion: Possible Follow-On Studies
In addition to gaining a better understanding of how the solution of IBVP (30) behaves for large-T, in particular, determining to what extent (if any) the recurrence behavior seen in Figures 2c  and 3c is related to the functional form of a given IC, future work on poroacoustic RRG theory could included the use of homogenization methods in problems wherein K and/or χ vary with position. Other possible extensions include the poroacoustic generalization of the study carried out in Ref. [25], wherein α was taken to be a function of (u x ) 2 , and also the case in whichμ is a power-law function of the shear rate tensor. Follow-on work might also include the study of poroacoustic signaling problems involving sinusoidal and/or shock input signals, as well as problems in which changes in entropy and temperature are taken into account.