## Abstract

The propagation of acoustic waves in a fluid that saturates a Darcy-type porous medium is considered under finite-amplitude theory. The equation of motion is derived, an acceleration wave analysis is carried out, and a travelling wave solution (TWS) is obtained. In addition, analytical findings are supported with numerical work generated by a simple, but effective, finite-difference scheme and results obtained are compared with those of the nonporous and linear cases. Most notably, this analysis reveals the following: (i) that the equation of motion is a new, hyperbolic form of Kuznetsov's equation; (ii) that finite-time blow-up of the wave amplitude is possible even if dissipation is present; (iii) the presence of a porous medium increases, with respect to the nonporous case, the rate at which amplitude growth/decay occurs; (iv) in the case of porous media propagation, not all compressive acceleration waves suffer blow up; and (v) that there exists a connection between acceleration waves and TWSs.

## 1. Introduction

For over a century the study of propagating singular surfaces, and in particular, acceleration waves, has been of great interest to many areas of the physical sciences, most notably continuum mechanics (e.g. Truesdell & Toupin 1960; Chen 1973; McCarthy 1975 and references therein). As opposed to a shock wave, which in the present context is defined as a propagating jump discontinuity (or jump) in at least one of the acoustic field variables, an acoustic acceleration wave is defined as a propagating jump in at least one of the first derivatives of the acoustic field variables, with the field variables themselves all being continuous. One of the main reasons why these waves are so interesting is that, in certain cases, the wave amplitude can suffer finite time blow-up, also known as ‘gradient catastrophe’ (Logan 1994) in the mathematical literature.

Thomas (1957) appears to have been the first to determine exact expressions describing the growth and/or decay of acceleration waves (or ‘sonic discontinuities’ as they are referred to by Thomas) in flows of inviscid ideal gases. Coleman & Gurtin (1967) studied acceleration waves and mild discontinuities in fluids that exhibit mechanical dissipation via the relaxation of internal state variables. These authors went on to conjecture that acceleration wave blow-up implied that a shock wave had, in fact, formed. They noted, however, that a rigorous mathematical proof of their supposition was lacking. Lindsay & Straughan (1978) examined the evolutionary behaviour of acceleration waves of arbitrary shape in a perfect fluid. More recently, Catheline *et al*. (2003) reported the results of an experimental study in which a transverse shock wave (technically referred to as a vortex sheet; McCarthy 1975) was observed forming in a sample of the soft tissue model *Agar-gelatin phantom* as the result of applying a sinusoidal boundary signal. Still more recently, Jordan & Christov (2005), who considered lossless finite-amplitude acoustic waves in organ pipes, have provided numerical support for the ‘shock-conjecture’ of Coleman & Gurtin (1967). It should also be mentioned that acceleration wave analysis has been, and remains, an important tool in theoretical investigations of the various formulations of thermoelasticity that have been proposed since the 1960s (e.g. Lindsay & Straughan 1976; Quintanilla & Straughan 2004).

In this work, we examine the propagation of finite-amplitude acoustic waves in a fluid that saturates a fixed, rigid, homogeneous and isotropic porous media. To incorporate the effects of the pores on the acoustic field within the fluid, we make use of Darcy's law, which is given by (e.g. Nield & Bejan 1999)(1.1)where * v* is the velocity vector,

*P*is an intrinsic pressure, and the positive constants

*μ*,

*K*and

*Χ*(<1), respectively, denote the dynamic viscosity, permeability and porosity. Here, we note that the filtration (or Darcy) velocity vector is related to the intrinsic (i.e. volume-averaged over a volume element consisting of fluid only) velocity vector

*by the Dupuit–Forchheimer relationship =*

**v***Χ*

*and it should also be noted that Darcy's law does not, in fact, represent a balance of forces (Nield & Bejan 1999). (Mathematically, Darcy's law is equivalent to a body force term.)*

**v**Notable early works in this area, which is often referred to as ‘poroacoustics’, include the book by Morse & Ingard (1968) in which the equation of motion for linear acoustic propagation in porous media is derived,1 and the paper by Emanuel & Jones (1968) in which the flow of a compressible fluid through a porous plate is considered. Later, Pascal (1986) studied linear acoustic propagation in porous media, along with deriving the equation of motion under Forchheimer's (nonlinear) flow law, and Nield (1994) proposed a modification to the Forchheimer–Brinkman equation, which describes the high-speed flow of compressible fluids in saturated porous media. More recently, Fellah & Depollier's (2000) investigated acoustic propagation in rigid porous media based on a fractional calculus model. (See also Fellah *et al*. 2003 and references therein for industrial applications involving poroacoustics.) For a listing of recent books and review articles on conventional (i.e. nonporous) finite-amplitude acoustics, see Jordan (2004).

The primary aim of the present work is to carry out, for the first time that the author is aware of, an in-depth study of acoustic acceleration waves in Darcy-type porous media based on finite-amplitude theory. In addition, we compare/contrast our findings with those of the linear theory and the nonporous case. To this end, the present article is arranged as follows: in §2, the equation of motion under the finite-amplitude approximation is determined. In §3, an acceleration analysis is carried out. In §4, numerical results are presented to illustrate how the findings of §3 relate to a model system. In §5, travelling wave solutions (TWSs) are investigated. And finally, in §6 results are summarized and a discussion follows.

## 2. Equation of motion under finite-amplitude theory

The equations of continuity, momentum, energy and state for a viscous, homogeneous, compressible fluid that is flowing homentropically in a fixed and rigid, non-thermally conducting, homogeneous and isotropic porous medium of permeability *K* and porosity *Χ* are(2.1)(2.2)(2.3)(2.4)where is the mass density, *℘* is the thermodynamic pressure, *η* is the entropy per unit mass, a superposed dot denotes the material derivative and all body forces have been neglected. Moreover, we note the homentropic nature of the flow, i.e. the fact that *η*=*η*_{0}(=const.) for all *t*≥0, implies that it is isentropic as well (Ockendon & Ockendon 2004). By equilibrium state, we mean the unperturbed state in which * v*=(0,0,0), ,

*℘*=

*℘*

_{0}and

*η*=

*η*

_{0}, where and

*℘*

_{0}are positive constants.

Restricting our attention to irrotational (i.e. ∇×* v*=0) flows, from which

*=∇*

**v***ϕ*follows, we find, on dividing through by and shifting the pressure gradient term to the left-hand side, that equation (2.2) reduces to(2.5)where

*ϕ*is the scalar velocity (or acoustic) potential.

In the case of finite-amplitude acoustic disturbances, it is assumed that(2.6)Here, *ς*∈(0,1) is a (sufficiently) small constant; is the condensation, where is the acoustic density; *p*=*℘*−*℘*_{0} is the over pressure; and the constant *c*_{0}(>0), known as the adiabatic sound speed, denotes the sound speed in the undisturbed fluid. (In perfect gases, , where *γ*=*c*_{p}/*c*_{v} is the adiabatic index and the constants *c*_{p}>*c*_{v} denote the specific heats at constant pressure and volume, respectively.) Now, the Taylor series expansion of the (isentropic) equation state about the equilibrium value is (e.g. Beyer 1997)(2.7)where the ratio *B*/*A* is a positive constant known as the nonlinearity parameter (Beyer 1997) and is equal to the adiabatic (bulk) elastic modulus. Substituting this expansion into equation (2.5), and using the fact that , yields, after simplifying,(2.8)In equation (2.8), is the fluid's kinematic viscosity and the coefficient of nonlinearity *β* is defined as (Beyer 1997)(2.9)Next, expanding both (1+*s*)^{−1} factors and employing the following (non-dimensional) transformations:(2.10)where *V*(>0) and *L*(>0), respectively, denote a characteristic flow speed and length, equation (2.8) becomes, after simplifying,(2.11)where all primes have been suppressed but are understood. Here, is the flow's Mach number, where we recall that , and we note that . What is more, the Darcy term's coefficient has become *δ*=*Χ*/Re, where Re=*c*_{0} *l*/*ν* is a Reynolds number and we have defined *l*≡*K*/*L* for convenience. We also note that the continuum assumption demands Kn<0.01, where is the flow's Knudsen number (Zucrow & Hoffman 1976). Assuming, in addition, that , we find that under the finite-amplitude approximation, whereby terms of order are neglected, equation (2.11) reduces to(2.12)the direct integration of which gives(2.13)where is an arbitrary function. Regarding equation (2.13) as a quadratic in *s*, solving for the root corresponding to the (physically relevant) ‘+’ case, and then expanding the result using the binomial theorem and the fact that is small, gives(2.14)where we have taken .

Returning now to equation (2.1), recasting it in terms of *s* and *ϕ* and employing the dimensionless transformations given in equation (2.10), yields, after rearranging terms and expanding the resulting (1+*s*)^{−1} factor,(2.15)Finally, replacing *s* in equation (2.15) with the right-hand side of equation (2.14) and then simplifying, we obtain, after neglecting terms of order , the equation of motion(2.16)which we observe is the porous media generalization of the lossless version of Kuznetsov's equation (e.g. Kuznetsov 1971; Chester 1991; Hamilton & Morfey 1997; Naugolnykh & Ostrovsky 1998). Now, when the velocity vector is of the form * v*=(

*u*(

*x*,

*t*),0,0), equation (2.16) reduces to(2.17)Since

*Δ*is always positive under the finite-amplitude approximation, where is the discriminant of equation (2.17), it follows that our equation of motion is strictly hyperbolic.

## 3. Acceleration waves

### (a) Hadamard's lemma and derivation of jump amplitude equation

In this section, we employ the theory of singular surfaces to obtain expressions for the speed and amplitude of propagating jump discontinuities in the first derivatives of the acoustic field variable *u*. The power of singular surface theory, in cases where it can be applied, is that this information can be determined without need of solving the partial differential equation (PDE) in question.

As the first step in this process, we employ the well-known relations(3.1)where the latter follows from equation (2.14) on neglecting terms of order , to recast equation (2.17) as the first order system(3.2)where *x* and/or *t* subscripts are now used to denote partial differentiation and(3.3)Clearly, the eigenvalues of *A* are(3.4)where *λ*_{1,2} are the roots of the characteristic equation det(*A*−*λI*_{2})=0 and *I*_{2} denotes the 2×2 identity matrix. Since and *λ*_{1}≠*λ*_{2}, it is evident that equations (3.2) comprise a strictly hyperbolic system, with characteristics defined by d*x*/d*t*=*λ*_{1,2} (Logan 1994). (Note that this is consistent with the fact that *Δ*>0.)

Consider now a smooth planar surface *x*=*Σ*(*t*) that is propagating along the *x*-axis of a Cartesian coordinate system, say in the positive *x*-direction. Let *Σ* be advancing into a region occupied by a fluid-saturated porous medium of the type described in §3. (Following the usual convention, we use a ‘+’ superscript to denote to the region into which *Σ* is advancing and a ‘−’ superscript to denote the region behind *Σ*.) Suppose, further, that *u* and *s* are both continuous functions of *x*, *t* within in the regions and , but that one of their first derivatives, say *u*_{t}, suffers a jump across *Σ*, i.e. [*u*]=[*s*]=0 but [*u*_{t}]≠0, where the amplitude of the jump in a function *F*=*F*(*x*,*t*) across *Σ* is defined as(3.5) are assumed to exist, and we observe that [*F*] and *F*^{+} and *F*^{−} are at most functions of *t* only. If *F*^{−}≠*F*^{+}, then *Σ* is termed a singular surface with respect to *F*. Since, in our particular case, *F*=*u*_{t}, it follows that *Σ* is classified as an acceleration wave (Chen 1973). Hence, given the value of [*u*_{t}] at *t*=0, and assuming that the fluid in is in its equilibrium state, we set ourselves the task of determining the behaviour of [*u*_{t}] for all *t*>0.

To this end, we begin by stating Hadamard's lemma (e.g. Chen 1973; Bland 1988), a result of fundamental importance in singular surface theory. For the problem under consideration, Hadamard's lemma takes the form(3.6)where |*U*|(≠0) is the speed of *Σ* with respect to the fluid ahead and denotes the (one-dimensional) displacement derivative with respect to *Σ* (Chen 1973). (Here, we should note that Hadamard's lemma is also known as Maxwell's theorem when [*F*]=0). Now, since [*u*]=[*s*]=0, then by equation (3.6) we have(3.7)(3.8)Moreover, since equation (3.2) holds on both sides of *Σ*, we also have(3.9)(3.10)Here, we note that in deriving the last two jump equations, the formula for the jump of a product (e.g. Lindsay & Straughan 1976)(3.11)was used along with the assumption that the fluid occupying the region is in its equilibrium state, i.e. *u*^{+}=*s*^{+}=0, to show that [*uu*_{t}]=[*ss*_{t}]=0.

Clearly, non-zero values of *U* that satisfy the above system of four jump equations exist only if the determinant of the corresponding coefficient matrix is zero. As a result, we are led to consider the Fresnel–Hadamard propagation condition(3.12)which is equivalent to the quadratic equation *U*^{2}−1=0, and it is therefore easy to see that *U*=±1. However, for the purposes of the present study, we need only consider the right-travelling case (i.e. the positive solution). From this, it is clear that *Σ*(*t*) corresponds to the moving plane *x*=*t*+*x*_{0}, where *x*_{0} is a constant.

Having determined the propagation speed, our next step is to derive the equation governing [*u*_{t}]. To this end, we first record the following equalities(3.13)all of which follow from equations (3.7)–(3.10). Next, we differentiate equation (3.2) with respect to *t* and then take jumps. After some rearranging, this results in(3.14)(3.15)In deriving equation (3.15), we used the fact that and , which respectively denote the values of *s*_{x} and *s*_{t} ahead of *Σ*, are both zero in the equilibrium state. Omitting the details, it is a relatively straightforward process to show, using Hadamard's lemma and equations (3.13)–(3.15), that the jump amplitude is governed by following Bernoulli-type ordinary differential equation (ODE):(3.16)

### (b) Evolutionary behaviour of acceleration waves

The exact solution to this version of Bernoulli's equation is well known and has been exhaustively studied by Chen (1973). Consequently, it is not difficult to show that(3.17)where the positive constant is known as the critical amplitude of the acceleration wave and [*u*_{t}]_{0}(≠0) denotes the value of [*u*_{t}] at *t*=0.

Based on equation (3.17), we find that [*u*_{t}] can evolve over time in any one of the following four possible ways:

If [

*u*_{t}]_{0}<0, then [*u*_{t}]∈[−|[*u*_{t}]_{0}|,0), for all*t*≥0, and [*u*_{t}]→0 monotonically from below as*t*→∞.If 0<[

*u*_{t}]_{0}<*α**, then [*u*_{t}]∈(0,[*u*_{t}]_{0}], for all*t*≥0, and [*u*_{t}]→0 monotonically from above as*t*→∞.If [

*u*_{t}]_{0}=*α**, then [*u*_{t}]=*α** for all*t*≥0.If [

*u*_{t}]_{0}>*α**, then , where*t*_{∞}(>0) is given by

(3.18)By contrast to the nonporous (i.e. *δ*→0) case (see §3*c*(ii)), where [*u*_{t}] is an algebraic function of *t* that can only grow or decay, it is clear that the amplitude behaviour is exponential and that four evolutionary paths are possible, of which two correspond to [*u*_{t}]→0, in the poroacoustics context. Hence, we see that one effect of the presence of a porous medium is to increase the rate at which the growth/decay of the amplitude occurs, as compared with that of the nonporous case. It should also be mentioned that case (iii), where [*u*_{t}] maintains a constant amplitude, is not possible in the nonporous context.

Finally, an acceleration wave is expansive (respectively, compressive) if [*u*_{t}]<0 (respectively, [*u*_{t}]>0) (see Chen 1973). From this, it is not difficult to establish the following:

If [

*u*_{t}]_{0}<0, then*Σ*is expansive for all*t*≥0.If [

*u*_{t}]_{0}∈(0,*α**], then*Σ*is compressive for all*t*≥0.If [

*u*_{t}]_{0}>*α**, then*Σ*is compressive for all*t*∈[0,*t*_{∞}).

From case (B) we see that cases (ii) and (iii) describe compressive acceleration waves with bounded amplitudes. This is in contrast to the nonporous case where *all* compressive acceleration waves suffer finite-time blow-up (again, see §3*c*(ii)).

### (c) Special/limiting cases

#### (i) Linear case

The linearized version of equation (2.17) is the PDE(3.19)which is a special case of the telegraph equation known as the damped wave equation (DWE). Like equation (2.17), the DWE is hyperbolic. Unlike equation (2.17), however, the DWE arises in the mathematical analysis of diverse physical phenomena ranging from the motion of DNA strands to ‘second sound’ propagation in thermally conducting media (e.g. Mickens & Jordan 2004 and the references cited therein). The following expression for the acceleration wave amplitude corresponding to the DWE is obtained by letting in equation (3.17)(3.20)a result that has been noted in the literature of many fields (e.g. Jordan *et al*. 2000 and references therein). Here, we observe that blow-up in this case is impossible and that *Σ* is expansive (respectively, compressive) if [*u*_{t}]_{0}<0 (respectively, [*u*_{t}]_{0}>0).

#### (ii) Nonporous case

In this case, the acceleration wave amplitude, which is obtained by letting *δ*→0, is given by the algebraic expression (e.g. Jordan & Christov 2005)(3.21)Thus, we see that if [*u*_{t}]_{0}<0, then the following are all true: (a) [*u*_{t}]<0 for all *t*≥0, (b) *Σ* is expansive for all *t*≥0 and (c) [*u*_{t}]→0 monotonically from below as *t*→∞. If, on the other hand, [*u*_{t}]_{0}>0, then , where(3.22)and *Σ* is compressive for all *t*∈[0,†_{∞}) since [*u*_{t}]>0 for all *t*∈[0,†_{∞}).

#### (iii) A new form of the Lighthill–Westervelt equation

Since *Σ* is a planar wavefront propagating with unit speed, it follows that equation (2.17) is, in the present acceleration wave context, equivalent to the PDE(3.23)which is obtained from the former by employing the approximation *ϕ*_{x}≈−*ϕ*_{t} only on the right-hand side (i.e. affecting only the ‘small’ terms). Equation (3.23) is a dissipative version of the hyperbolic Westervelt, or Lighthill–Westervelt, equation (e.g. Jordan & Christov 2005 and references therein). More importantly, however, with the elimination of the mixed derivative term from equation (2.17), we now have a much simpler equation of motion. Consequently, in §4, where acceleration waves are numerically simulated, we can justifiably use equation (3.23) in place of equation (2.17), greatly simplifying our tasks. (See also Hamilton & Morfey 1997, pp. 52–56, and references therein.)

## 4. Model system: analytical and numerical results

We now apply the theoretical results developed thus far to the study of a physical system. In addition, computational method/tools are employed to numerically solve the resulting initial boundary-value problem (IBVP) and to illustrate the analytical findings given in §3. In particular, we construct a simple finite-difference scheme that is surprisingly effective in capturing the shocking-up behaviour described in case (iv).

### (a) Problem formulation

To this end, we recast equation (3.23) in terms of the non-dimensional over-pressure and consider the following IBVP(4.1a)(4.1b)(4.1c)where *H*(.) is the Heaviside unit step function, *ω* is a dimensionless frequency, *n*=0,1 as demanded, *x*_{0} is now zero, and all primes are again suppressed. Furthermore, so as to avoid changes in the waveform due to reflection, we have limited our study to *Σ*'s initial transit of *x*∈(0,1), hence the restriction *t*<1. Here, we note that the (amplitude-dependent) wave speed predicted in equation (4.1*a*) is(4.2)from which the well-known result is easily obtained (e.g. Naugolnykh & Ostrovsky 1998).

Perhaps the simplest physical system modelled in this IBVP is that of acoustic propagation in a muffler pipe caused by the input of a sinusoidal signal at one end. (See Crighton & Scott (1979) for a discussion of this boundary condition in the context of acoustics.) Specifically, suppose the muffler's length is *ℓ* and that it is filled with a porous material of the kind described in §2. Initially, the air inside is in its equilibrium state. Beginning at *t*=0, however, a sinusoidal signal of constant amplitude and frequency is inserted at the end *x*=0, while the end at *x*=*ℓ* is open to the atmosphere. Assuming that the open end is maintained in the equilibrium state, it follows that , where is the amplitude of the input signal; *L*=*ℓ*; and , where we recall that ∈(0,1) is small.

### (b) Acceleration wave analysis of IBVP

From the results of §3*b*, it is not difficult to show that the jump amplitude and breakdown time are, in the case of IBVP (4.1*a**–c*), given by(4.3)where we have used the fact that [*p*_{t}]=−[*p*_{x}],(4.4)and the amplitude of the jump in the time-derivative of the boundary condition at *x*=0 across the plane *t*=0 is(4.5)In addition, we observe that(4.6)This clearly illustrates the fact that not only must *Σ* be compressive, but *ω*>*α*^{*} must hold as well if finite-time blow-up is to occur. (Actually, there is no other possible way that blow-up can occur here other than in finite time.)

### (c) Construction of finite-difference scheme

Since an analytical solution of IBVP (4.1*a**–c*) appears to be out of the question, we therefore must turn to numerical methods if further progress is to be made. Hence, we consider, guided by the work of Jordan & Christov (2005) on the *δ*≡0 case, the following simple discretization of equation (4.1*a*):(4.7)where(4.8)Since equation (4.7) is linear in the most advanced time-step approximation , we can solve for this term to obtain, after some manipulation, the (explicit) scheme:(4.9)where *R*=(Δ*t*)/(Δ*x*) and the truncation error is of order .

### (d) Linearized IBVP

In the linearized case, equation (4.1*a*) reduces to the DWE *p*_{xx}−*p*_{tt}−*δp*_{t}=0 and the corresponding exact solution can be determined using the Laplace transform. For the purposes of the present study, however, the following simple expression, which very closely approximates the former during *Σ*'s initial transit of the muffler pipe, will be adequate (e.g. Jordan *et al*. 2000):(4.10)

### (e) Numerically generated pressure profile plots

To illustrate the evolution and behaviour of *Σ*, including its apparent transition from acceleration to shock wave, we have, in figures 1–3, plotted the numerically computed solution of IBVP (4.1*a**–c*). Here, we have taken the saturating fluid to be air at 20 °C. As such, it follows that *γ*=1.402 (⇒*β*=1.201), *c*_{0}=343 m s^{−1} and *ν*=1.496×10^{−5} m^{2} s^{−1} (see appendix 10, table (c), of Kinsler *et al*. 1982).

Let us now consider the porous material that fills the muffler pipe. For simplicity, we take it to consist of packed beds of rigid solid spheres, all of diameter *d*, which are fixed in place. In this configuration, the permeability is given by the well-known Kozeny–Carmen relation (Nield & Bejan 1999)(4.11)If we further suppose that the spheres are packed in a cubical fashion, then *Χ*=*Χ*_{c}, where *Χ*_{c}≡0.4764 (Nield & Bejan 1999), and the damping coefficient becomes(4.12)If we also take *ℓ*=0.5 m, which is about the length of a typical muffler found on small off-road recreational vehicles, and *d*=4.5 mm, then *δ*≈0.2341 and it follows that the continuum assumption is valid since Kn≈4.21×10^{−4}. Additionally, let us choose *ω*=3.7, which corresponds to a dimensional frequency of approximately 2.54 kHz, so that the sinusoidal nature of the input signal will be apparent.

The curves shown in bold, which correspond to , were produced from datasets computed by a simple algorithm, which implemented our scheme (equation (4.9)) on a desktop PC running Mathematica (version 5.0). Interpolations between the points were then accomplished using the built-in cubic splin routine that is part of this software package. To help illustrate the decay and growth of the acceleration wave amplitudes, we have also included the tangents to the solution profiles at *x*=*t*, which are generated from equation (4.3) and plotted as broken line segments. And to allow for comparisons of the nonlinear versus linear versions of IBVP (4.1*a**–c*), we have also graphed equation (4.10) as the thin-solid curves.

The sequences shown in figures 1–3 illustrate the evolution of the solution profiles during the acceleration wave's initial transit of the muffler and correspond to cases (i), (iii) and (iv), respectively, of §3*b*. (A figure corresponding to case (ii) has not been given since, qualitatively, cases (i) and (ii) describe the same behaviour.) Figures 1 and 3 and figures 2 and 3 illustrate the effects of changing *n* and , respectively. Here, it should be noted that in producing figure 3, which depicts most of the finite-time blow-up of [*p*_{x}], we required *t*_{∞}=1 so that *Σ* would live out its entire life during the course of its initial passage of *x*∈(0,1). Consequently, a Mach number of was needed. While this value of may be a bit large, considering that we are working under the finite-amplitude approximation, it does place the flow within the range where compressibility effects are becoming important (Liggett 1994) and, as we will soon see, provides for clearly distinguishable plots of the nonlinear versus linear solution profiles appearing in figures 1 and 3.

In figure 1, where |*p*| vs. *x* is plotted since *n*=0, we observe, as predicted in case (i), that the magnitude of the slope of the profile at the wavefront decreases as *t* increases and that the propagation speed of *Σ* is unity, in agreement with the fact that *U*=1. We also observe that the nonlinear profile lags behind that of the linear case, a result of the fact that *c*(*p*)<1 when *n*=0 (see equation (4.2)), and that this lag becomes more pronounced as *t* increases.

On the other hand, from figure 2, in which *p* versus *x* is shown for *n*=1 and *ω*=*α**, the fact that is a constant (i.e. −*α**) is very evident. We also see that now the linear and nonlinear profiles are much closer together, essentially overlapping for small *t*. This is, of course, due to the fact that figure 2 was generated using the smaller Mach number value of . However, for larger values of *t*, we see that the nonlinear profile has pulled ahead, albeit slightly, of its linear counterpart, a reflection of the fact that *c*(*p*)>1 when *n*=1.

Figure 3, which like figure 2 plots *p* versus *x* for *n*=1, but for *ω*>*α** captures about 95% of the life span of *Σ* under case (iv). The behaviour shown is due to the amplitude dependence of the wave speed, where as in figure 2 *c*(*p*)>1, and involves a rapid increase in |[*p*_{x}]| as *t*→*t*_{∞}. Consequently, it appears from the last frame of figure 3 that the solution profile is preparing to ‘break’ as *Σ* nears *x*=1, thus suggesting that a jump in *p* itself (i.e. a shock wave) is forming.

In closing this section, we observe that figures 1–3 exemplify how the effects of nonlinearity are cumulative in nature, even for small (see figure 2). This phenomena, clearly a result of the amplitude dependence of *c*, is especially apparent in figures 1 and 3, where the relatively large value of used gives rise to a rapid departure from the linear profile. Finally, we note that our scheme has accurately captured both the attenuation of the profile caused by the Darcy term *δp*_{t} and the corner at *x*=*t* corresponding to *Σ*.

## 5. Travelling wave solution

We begin our study of TWSs with the following observation: since equation (2.17) is invariant under the transformation *x*→−*x*, it is clear that we need only consider, without loss of generality, the case of right-travelling waves. Thus, we seek solutions of the specific form(5.1)where *ξ*≡*x*−*v*_{0}*t* and the positive constant *v*_{0} denotes the (non-dimensional) propagation speed. To this end, we substitute equation (5.1) into (2.17), integrate once with respect to *ξ*, rearrange terms, and simplify to arrive at the nonlinear ODE(5.2)where is a positive constant, , and *C*_{1} is the constant of integration. While equation (5.2) can be integrated exactly as it stands, for the purpose of the present investigation we will limit our attention to the special case *v*_{0}=*U*(=1), and thus consider the simpler Chrystal equation . This ODE is easily solved to yield the exact, but unbounded, non-trivial TWS(5.3)where *C*_{2} is the second constant of integration. Equation (5.3) describes an undamped, undistorted, inverted parabolic waveform which propagates to the right at a speed of *v*_{0}=1, i.e. at the (non-dimensional) sound speed in the undisturbed fluid. Now, while it appears from this solution that |*f*| and |d*f*/d*ξ*| both tend to infinity, as , when *v*_{0}=1, we note with interest that(5.4)That is, we have found a connection between acceleration wave theory and TWSs; specifically, for the special case *v*_{0}=*U*, the acceleration expression derived from equation (5.3) is equivalent to the jump amplitude solution cited in case (iii) of §3*b*, but only for *x*<*Σ*. The resolution of this apparent anomaly, along with a detailed analysis of equation (5.2), will be published elsewhere.

## 6. Summary and final remarks

In this paper, a study of nonlinear acoustic acceleration waves in Darcy-type porous media has been presented. Our investigation, which extends the work of Jordan & Christov (2005) to the area of poroacoustics, has revealed how acceleration waves propagate and evolve over time in such media and, via case (iii), has uncovered a link between singular surfaces and TWSs. We have also compared/contrasted our findings with those of the nonporous and linear cases and, employing both analytical and numerical methods, we have applied our results to a model system that approximately describes acoustic propagation in a muffler pipe. The numerical work presented here not only provides additional support for the shock conjecture of Coleman & Gurtin (1967), but also demonstrates how effective a simple finite-difference scheme and relatively inexpensive computational hardware can be in simulating transient phenomena associated with nonlinear, hyperbolic PDEs.

Of the findings reported, case (iii) is perhaps the most mathematically intriguing. In particular, we observe that under the conditions of case (iii), the effects of dissipation and nonlinearity cancel each other out in the context of acceleration waves. This is based on the observation that if equation (4.1*a*) was replaced with the linear wave equation *p*_{xx}−*p*_{tt}=0, then the same jump amplitude expression would result when *ω*=*α**. It must be stressed, however, that case (iii) is unstable in the sense that any discrepancy, however small, in the initial value of the jump amplitude will lead to either case (ii) or (iv).

Let us briefly return to the TWS analysis presented in §5. If we differentiate equation (5.2) with respect to *ξ* and recast it in terms of , where , then requiring that , where are constants, leads to the condition . From this, we conclude that equation (2.17) is an example of a dissipative, integrable PDE that does not admit a TWS for in the form of a Taylor shock profile.2 This is an interesting result in light of the fact that Kuznetsov's equation(6.1)where Re_{d} is a Reynolds number involving Lighthill's diffusivity of sound, has recently been shown to possess, in terms of , exact TWSs of the tanh [.] type (Jordan 2004). Now, equations (2.17) and (6.1) are clearly equivalent when their respective dissipation terms are neglected. Note, however, that Kuznetsov's equation, which describes acoustic disturbances in thermo-viscous fluids, is not a wave equation in the strict sense; it has a diffusive character due to the third-order damping term. In contrast, dissipation in equation (2.17) arises from the first-order term *δϕ*_{t}, and thus the hyperbolicity of equation (2.17) is left intact (see §2). Clearly, the form of the dissipation term can have a major impact on the nature of the solution(s).

Lastly, we must emphasize that hyperbolicity does not preclude a Taylor shock-type profile TWS. For example, consider the following undamped, nonlinear PDEs:(6.2)where *m* and *g*(>0) are constants and *φ* is a dummy variable. The first is the famous sine-Gordon equation of soliton theory and the other is the *φ*^{4} (or phi-four) equation of particle physics. While both are obviously hyperbolic, they nevertheless admit TWSs of the form tan^{−1} [.] and tanh [.], respectively (e.g. Dodd *et al*. 1982). However, we should point out that, unlike equation (2.17), for which *λ*_{1,2} are amplitude-dependent (see equation (3.4)), the characteristics of the two PDEs in equation (6.2) are simply the straight lines defined by d*x*/d*t*=±1. Consequently, they, like the linear wave equation and the DWE, are incapable of exhibiting the type of acceleration wave blow-up that has been discussed here.

## Acknowledgments

The author would like to express his sincere thanks to Professors S. A. Chin-Bing, M. J. Ablowitz, C. I. Christov, B. Straughan, and H. and J.R. Ockendon for the helpful discussions and for sharing their insight. Some of the findings reported in this article were part of a lecture that the author delivered both at the Oxford Centre for Industrial and Applied Mathematics (OCIAM), Oxford University, Oxford, UK and at Grey College, University of Durham, Durham, UK on 13 and 20 February 2004, respectively. P.M.J. was supported by ONR/NRL funding (PE 061153N).

## Footnotes

↵Since it is expressed as a function of the porosity, the wave velocity coefficient in their equation (6.2.24) is, unfortunately, incorrect.

↵Broadly speaking, any strictly monotonic, bounded function of

*ξ*that approaches constant, but unequal, limits as ; also known as a ‘diffusive soliton’ (see Jordan 2004 and references therein).- Received November 17, 2004.
- Accepted March 21, 2005.

- © 2005 The Royal Society