We present a physical model and a numerical method based on a space- and time-dependent Galilean-type coordinate transformation to simulate acoustic waves in the presence of an accelerating background flow field with sonic transition. Kinematically, the coordinate transformation is designed so as to maintain the well-posedness of the transformed wave equation, which is solved in a fixed computational domain using standard finite differences. Considering an acoustic black hole analogy, we analyze the nonlinear dynamics of acoustic waves in a stationary but non-uniformly accelerating flow field under the assumption of spherical symmetry. The choice of the acoustic black hole analogy is motivated by the fact that the steady-state spherical sonic horizon allows us to parameterize the wave-flow configuration in terms of a Helmholtz number He=c2/(λagh), which is expressed as a function of the speed of sound c, the emitted wavelength λa, and the flow acceleration at the sonic horizon, that is, the acoustic surface gravity gh. The results of the numerical simulations show that He describes geometrically similar sets of wave characteristics for different combinations of gh and λa. However, we also observe nonlinear variations of the wave amplitude along the wave characteristics, which are attributed to nonlinear Doppler modulations. It appears that these amplitude modulations depend on the acceleration of the flow field and can, therefore, differ for geometrically similar characteristics.

The nonlinear behavior of acoustic waves emitted or scattered by accelerating boundaries has recently attracted attention as nonlinear amplitude modulations and frequency side-bands can reveal distinct information about the motion of a wave emitter or scatterer. A comprehensive review of different (potential) applications is given in the work of Gasperini et al.,1 who emphasize the difficulty of predicting nonlinear Doppler shifts and propose a computational method based on a suitable coordinate transformation that can account for the accelerating motion of wave scattering boundaries. Using a similar approach, we recently confirmed the analytical predictions by Christov and Christov2 regarding the nonlinear Doppler waveform of an oscillating emitter by the means of numerical simulations.3 However, the background medium is assumed to be at rest in the aforementioned studies, whereas it is known that Doppler shifts are also caused by the motion of the background medium itself. This is harnessed, for instance, in medical Doppler ultrasonography to estimate the velocity of blood flows or moving tissue and to identify pathological patterns in the Doppler spectrum.4,5 While considerable research efforts have been made to improve Doppler signal processing techniques for medical applications,6 a more fundamental understanding of nonlinear Doppler modulations in moving flow fields is still missing.

In the present work, we derive a physical model and present a computational method to simulate acoustic waves in moving background flow fields. We then employ the acoustic black hole analogy in order to investigate the dynamics of an acoustic wave propagating in an accelerating background flow field.

An acoustic black hole is a convergent fluid flow with sonic transition, where the point of sonic transition forms an acoustic event horizon from which small signal sound waves cannot escape.7 Unruh8 discovered that the propagation of sound waves in an irrotational and inviscid fluid flow, and in particular, the associated acoustic metric, corresponds exactly to the motion of a scalar field in a curved (3 + 1)-dimensional Lorentzian space–time geometry.9 Since then, the acoustic black hole analogy and other back hole analogue systems have gained considerable attention because the kinematic connection to curved space–time geometry allows for the investigation of certain kinematic aspects of gravitational black holes, such as Hawking radiation8,10–13 or rotational super-radiance,14,15 in laboratory experiments.16,17 Apart from the investigation of astrophysical phenomena, the study of acoustic black holes may also contribute to an improved understanding of the acoustic wave behavior in accelerating background flow fields, which may, in turn, help to improve relevant technical applications such as medical Doppler ultrasonography.

In principle, the presence of a sonic transition is not essential to the wave dynamics discussed in this work, which merely depend on the flow acceleration. However, as it is further discussed in Sec. II A, a steady-state sonic horizon has the beneficial feature that it provides a fixed reference point of the propagating waveform that is subject to a constant background flow acceleration, even if the background flow field is non-uniform in space. This allows to parameterize the wave-flow configuration by combining a suitable Helmholtz number with the acceleration at the sonic horizon, which again depends on the size of the steady-state acoustic black hole. We will then make use of this parametrization to study the amplitude modulations along geometrically similar wave characteristics.

Starting from the work of Unruh,8 Visser,9 and Balbinot et al.,18 we first derive a quasi-one-dimensional wave equation for a moving and compressible background medium based on first-order perturbation theory (see Sec. II B). As we aim to study the dynamics of acoustic waves in a smooth background flow field, it is briefly motivated in Sec. II C how, under certain conditions, a smooth sonic transition is supported in de Laval nozzle flows19 and Bondi–Parker-type flows,20,21 where the latter describe the mass accretion and plasma wind flows around compact stellar objects.

In Sec. III, we derive a suitable coordinate transformation in physical space, conveyed by a generic change of variables, which is instrumental in maintaining the well-posedness of the wave equation and used as basis for a numerical solution. While this approach is well-established in the context of wave emitters moving relative to a quiescent background medium,1–3 the well-posedness of the governing equation is still a particular problem when the wave emitter approaches or even exceeds the sonic speed.2 Despite the fact that we consider a wave emitter of zero velocity relative to the moving background medium, so that no sonic shock can form, the governing wave equation becomes ill-posed for supersonic flows if the wave equation is formulated in a fixed Eulerian reference frame, which is passed by the supersonic flow. Therefore, a space- and time-dependent coordinate transformation recovers the proper two-sided wave propagation behavior at any point in physical space. This idea evolves along similar lines as the utilization of Galilean and Lorentz invariants of the Navier–Stokes equations and the wave equation as presented by Gregory et al.22 and Wang.23 Gregory et al.22 pointed out that the description of Doppler phenomena in moving fluids can greatly simplify under the appropriate coordinate transformation.

Finally, we consider the canonical case of a spherically symmetric steady-state acoustic black hole, for which quasi-planar and linear wave propagation is obtained in the limit of very large radii of the sonic horizon. The focus is on the geometrical similarity of the wave characteristics as described by the Helmholtz number, and the effect of the sonic horizon's surface gravity, that is, the background flow acceleration, on the acoustic wave dynamics.

In the following, we discuss important parameters that govern the acoustic wave dynamics in a stationary spherical flow with sonic transition. Figure 1 (left) shows the schematic of a spherical sonic horizon at radial position r=rh, where the magnitude of the inward directed radial flow velocity u coincides with the sonic speed c. In the context of the present work, R and Ṙ are the instantaneous position and velocity, respectively, of a virtual spherical wave emitting boundary moving along with the flow as indicated in Fig. 1 (left). Thus, there is no relative motion between the background medium and the wave emitter, so that any nonlinear Doppler shifts are exclusively imparted by the non-uniform motion of the background medium. For the following derivations, it is also advantageous to write the radial velocity distribution as

u(r,t)=Ṙ(t)(R(t)r)n(r,t).
(1)
FIG. 1.

Schematic of a steady-state spherical acoustic black hole. At the sonic horizon as indicated by the radial position r=rh on the left, the velocity magnitude of the inward directed spherically symmetric background flow u0 is equal to the speed of sound c so that sound waves cannot escape to infinity from inside this region. The total flow velocity is u=u0+εu1, where εu1 is a small acoustic perturbation. On the right, the wave characteristics are depicted in a space–time diagram. The characteristics are curved due to the non-uniform background flow field, where the supersonic flow region is indicated by the gray color. The background flow can be thought of as a potential sink with singularity at the flow center. Far away from the center, the flow velocity gradient becomes small so that flat (non-divergent) characteristics are recovered.

FIG. 1.

Schematic of a steady-state spherical acoustic black hole. At the sonic horizon as indicated by the radial position r=rh on the left, the velocity magnitude of the inward directed spherically symmetric background flow u0 is equal to the speed of sound c so that sound waves cannot escape to infinity from inside this region. The total flow velocity is u=u0+εu1, where εu1 is a small acoustic perturbation. On the right, the wave characteristics are depicted in a space–time diagram. The characteristics are curved due to the non-uniform background flow field, where the supersonic flow region is indicated by the gray color. The background flow can be thought of as a potential sink with singularity at the flow center. Far away from the center, the flow velocity gradient becomes small so that flat (non-divergent) characteristics are recovered.

Close modal

The parameter n(r,t) in Eq. (1) can formally be both time- and space-dependent, and it can be specified in such a way that the continuity equation for a compressible spherically symmetric flow is satisfied. In the present work, we perform numerical solutions for the particular case of a stationary spherically symmetric flow, for which the continuity equation 4πr2ρ(r)u(r)=4πR2ρ(R)Ṙ yields

u(r)=ρ(R)ρ(r)Ṙ(t)(R(t)r)2,
(2)

which is equivalent to Eq. (1) for

n(r)=2+log(R/r)(ρ(R)ρ(r)).
(3)

Hence, we have n =2 for an incompressible flow and n2 for a compressible flow with monotonically increasing density toward the center.

For a given stationary velocity field u(r), R may be obtained as solution of Eq. (1) subject to initial conditions for R and Ṙ. In the context of the present work, this solution is required as a boundary condition for the governing wave equation and must, therefore, be obtained analytically. In order to facilitate the analytical solution of Eq. (1), we take n as spatial-temporally constant. This assumption is justified as we are only interested in a confined region around the sonic horizon and the principle effect of the magnitude of n. Furthermore, n may also be chosen in such a way that Eq. (1) represents a non-spherical flow geometry, as for instance the quasi-one-dimensional flow in a de Laval nozzle. This allows us to apply the following considerations to any quasi-one-dimensional flow with sonic transition.

With n=const, it can be inferred from the work of Visser9 that the acoustic surface gravity gh at the sonic horizon of a steady-state acoustic black hole associated with a quasi-one-dimensional flow follows as

gh=nc2rh+ccr|rh.
(4)

The term “surface gravity” reflects the analogy to the acceleration experienced by a freely falling observer at the event horizon of a gravitational Schwarzschild black hole.24 For a nearly constant speed of sound in the proximity of the steady-state sonic horizon, gh represents the acceleration of the stationary background flow field. The change of the fluid acceleration with respect to r is proportional to the acceleration itself and determines how quickly neighboring wave characteristics diverge away from each other as indicated in Fig. 1 (right). This suggests that the acoustic surface gravity is an important parameter that governs the rate at which an outward propagating wave is red-shifted.

The variation of the background flow field over a waveform further depends on the wavelength itself. Therefore, we adopt the Helmholtz number25 to characterize the size of the acoustic black hole relative to the emitted wavelength. With fa being the emitted wave frequency and λa=c/fa the emitted wavelength, the Helmholtz number is given by He=2πrh/λa. Substituting rh from Eq. (4) and dropping the dimensionless parameter n as well as the constant 2π, the Helmholtz number can be rewritten as

He=c2λa(1ghccr|rh)1,
(5)

meaning that the ratio c2/λa is uniquely defined by the values of He,gh, and cc/r|rh, where the latter contribution is negligible if c varies slowly over the acoustic wavelength. This relation suggests a geometrical similarity of the wave characteristics for a particular value of He and a given flow field as parameterized by n. A particular value of He can be obtained from different combinations of gh and λa. As the Helmholtz number depends, in this particular form, on the flow acceleration at a distinct reference point rather than the radius of a spherical wave emitter, its applicability is not restricted to spherical symmetry, but extends to any quasi-one-dimensional geometry, where the corresponding flow field can be parameterized by the scaling factor n in Eq. (1).

Because it is of particular relevance to the present work, the derivation of a wave equation describing the propagation of acoustic waves in a moving fluid is briefly repeated here, drawing on the works of Unruh,8 Visser,9 and Balbinot et al.18 We start from the three-dimensional representation, which is then simplified to the case of spherical symmetry at the end of this section. Neglecting viscous stresses, the dynamics of the fluid motion and the acoustic waves are governed by the Euler equations, that is, the conservation equations of mass and momentum, given by

ρt+·(ρu)=0,
(6)
ut+(u·)u=1ρpψ.
(7)

In Eqs. (6) and (7), u, ρ, and p are the fluid velocity, density, and pressure, respectively, and ψ is the gravitational potential. Further assuming the flow to be irrotational (×u=0), a velocity potential ϕ exists so that

u=ϕ.
(8)

Integration of Eq. (7) and using Eq. (8) gives the transient Bernoulli equation

ϕt+12(ϕ)2=p0pdpρ(p)ψ=h(p)ψ,
(9)

where h(p) is the difference of the fluid's enthalpy between the total pressure p and some undisturbed reference pressure p0. By introducing the first-order perturbations ϕ=ϕ0+εϕ1,p=p0+εp1,ρ=ρ0+ερ1, and u=u0+εu1 about the background quantities,9 where the subscript 0 indicates the background quantities and the terms scaled by ε are small amplitude acoustic perturbations, substituting into Eqs. (6) and (7), and further linearizing the enthalpy difference according to h(p0+εp)h0+εp/ρ0, we obtain

ρ0t+·(ρ0u0)=0,
(10)
ρ1t+·(ρ0u1+ρ1u0)=0,
(11)
ϕ0t+12(ϕ0)2+h0=ψ,
(12)
ϕ1t+u0·ϕ1+p1ρ0=0.
(13)

Solving Eq. (13) for p1 gives

p1=ρ0Dϕ1Dt,
(14)

where

DDt=t+u0·
(15)

is the material or Lagrangian derivative operator. Inserting Eq. (14) into the linearized barotropic relation dp/dρ=c2p1/ρ1 of the acoustic perturbation and then substituting for ρ1 in Eq. (11), one arrives at the wave equation9 

t(ρ0c2Dϕ1Dt)·(ρ0ϕ1u0ρ0c2Dϕ1Dt)=0.
(16)

Similar derivations based on first-order perturbations are found in the works of Ewert and Schröder26 and Godin,27 where Ewert and Schröder26 derived an additional term that represents the effect of rotational velocity components in addition to the irrotational flow field. In spherical coordinates and under the assumption of spherical symmetry, Eq. (11) is rewritten as ρ1/t+(1/r2)[r2(ρ0u1+ρ1u0)]/r=0, where r and u0 are the radial coordinate and the radial background flow velocity, respectively. Repeating the above steps analogously, the wave equation becomes

t(ρ0c2Dϕ1Dt)1r2r[r2(ρ0ϕ1ru0ρ0c2Dϕ1Dt)]=0,
(17)

where D/Dt=/t+u0/r is the Lagrangian/material derivative operator analogously to Eq. (15). Expanding the derivative operators in Eq. (17) gives

2ϕ1t2+2u02ϕ1rt+Du0Dtϕ1r2c2rϕ1r(c2u02)2ϕ1r2=c2ρ0ρ0rϕ1r+2cDcDtDϕ1Dt.
(18)

As the focus of the present work is on the effect of the background flow acceleration on the acoustic wave dynamics in a continuous flow field, we make the assumption of a stationary and smooth (shock-free) sonic transition. In this section, we briefly indicate under which conditions a smooth sonic transition is physically possible. In principle, a smooth sonic transition can be established in de Laval nozzles19 or in spherically symmetric Bondi–Parker-type flows.20,21 As pointed out by Williams and Dyson19 and briefly outlined in the following, the prototypical Bondi–Parker equations are very similar to the equations describing the quasi-one-dimensional steady-state flow in a de Laval nozzle.

We assume that the barotropic relation dp/dρ=c2 holds for the background fluid state, which is indeed justified under the assumption of a shock-free flow. With the surface area A(r) passed by the flow in the normal direction, a stationary, inviscid, and quasi-one-dimensional flow with mass flow rate ρ0u0A and gravitational potential ψ is governed by the continuity and momentum equations28,29

1AAr+1ρ0ρ0r+1u0u0r=0,
(19)
u0u0r=c2ρ0ρ0rψr.
(20)

Equations (19) and (20) can be combined to give the radial velocity gradient

u0r=c2u0u02c2(1AAr1c2ψr).
(21)

A singular velocity gradient (shock) can only be avoided if the parenthesized term on the right-hand side of Eq. (21) vanishes. The gradient of the gravitational potential in a Bondi–Parker-type flow is given by29 

ψr=GMr2,
(22)

where G is the gravitational constant and M the mass of a compact stellar object. The critical/sonic point is at rh=GM/(2c2),29–31 while (A/r)/A=2/r for a spherical flow, so that a smooth sonic transition at r=rh is indeed supported by the above equations.19 

The stationary quasi-one-dimensional flow in a de Laval nozzle is governed by the same set of equations, however, with ψ/r=0. Here, the parenthesized term vanishes at the point of sonic transition under the appropriate boundary conditions since A/r=0 at the nozzle throat, effectively resulting in non-convergent/non-divergent streamlines.19 Kovalenko and Eremin32 demonstrate the stability of the spherically symmetric transonic Bondi-type flow against spherically symmetric acoustic perturbations. They further suggest that the breakdown of spherical symmetry is responsible for the onset of instability and shock formation. By linear stability analysis, Ray and Bhattacharjee31 show that while the Bondi–Parker equations admit a stationary transonic solution associated with the state of lowest energy of the system,20 this solution is physically not realizable. Ray and Bhattacharjee31 further show that a transonic solution is only supported by linear stability analysis if the flow is transient. However, while the radial velocity profile becomes time-dependent under the assumption of a transient flow, the instantaneous velocity profiles are still qualitatively similar to the solution of the stationary Bondi–Parker equations.31 Therefore, one may assume that the timescale of the variation of the background flow velocity u0 is negligibly small as compared to the timescale of the acoustic wave problem. Maicke et al.33 discuss various distinct three-dimensional effects leading to the departure from the shock-free condition in a supersonic de Laval nozzle flow. Moreover, the realization of nearly quasi-one-dimensional conditions, precluding the occurrence of perturbations of the background flow field, requires superfluidic conditions as can be achieved with Bose–Einstein condensates34,35 or superfluidic helium II.36 

As is common in the context of acoustic black holes,7,8,18,37 we take the speed of sound as a constant in our analysis. This simplifying assumption is associated with a first-order approximation of the background flow compressibility, as it is shown by expanding p0(ρ0) around some reference density ρ0,ref under the assumption of a barotropic flow. With dp0/dρ0=c2, the Taylor series expansion yields

p0(ρ0,ref+Δρ0)=p0(ρ0,ref)+c2|ρ0,refΔρ0+cdcdρ0|ρ0,ref(Δρ0)2+O((Δρ0)3),
(23)

showing that the variation of c with respect to ρ0 is a second-order effect. With c=const, the second term on the right-hand side of Eq. (18) vanishes, c/r=0 in Eq. (4), and Eq. (5) reduces to

He=c2ghλa.
(24)

Again employing the critical point rh=GM/(2c2) of a spherical Bondi–Parker-type flow,29–31 as well as Eq. (4) for c=const, Eq. (22) is conveniently expressed as

ψr=2nc4r2gh.
(25)

Equation (20) readily gives the density gradient

ρ0r=ρ0c2(u0u0r+ψr).
(26)

Substituting Eq. (26) into Eq. (18) yields

2ϕ1t2+2u02ϕ1rt+Du0Dtϕ1r2c2rϕ1r(c2u02)2ϕ1r2=(u0u0r+ψr)(c2/ρ0)ρ0/rϕ1r.
(27)

We will use Eqs. (24), (25), and (27) in the further course of this study. A general advantage of the above formulation is that the effect of a variable and non-uniform background density field is absorbed in the right-hand side term that exclusively depends on the background flow field. The appropriate velocity distribution can then be treated as an input to the coefficients of the wave equation, and it can be obtained either analytically, numerically, or experimentally. However, the coefficient c2u02 of the Laplacian 2ϕ1/r2 in Eq. (27) changes its sign at the point of sonic transition so that the problem is ill-posed for regions where u02>c2. In Sec. III, we will introduce a time- and space-dependent coordinate transformation in physical space in order to extend the validity range of Eq. (27) to supersonic background flow fields.

Regarding the numerical procedure and the uniqueness of the solution, it is a particular problem that the coefficient of the Laplacian on the left-hand side of Eq. (27) switches its sign at the point of sonic transition, rendering the problem ill-posed. The same problem is pointed out by Christov and Christov,2 however, for the wave equation describing a moving wave emitter approaching the sonic shock condition in a quiescent background medium. In the present work, the acoustic wave does not form a shock as the wave emitter motion is aligned with the background flow, and the problem merely originates in the circumstance that the wave equation is formulated in a fixed Eulerian reference frame, which the fluid passes at the speed of sound at the sonic horizon. Hence, at the point of sonic transition, an outward propagating wave “freezes” and appears to have a zero frequency (infinite wavelength). In a reference frame moving with the background flow, by contrast, proper two-sided wave propagation at finite wavelength is observed. This behavior can be recovered by a suitable Galilean-type coordinate transformation, which emphasizes the fact that the well-posedness is not primarily a property of the governing equation itself, but rather of the reference frame in which it is formulated.

In order to switch to a more suitable reference frame to describe the wave dynamics for transonic background flows, we introduce the moving physical domain Ω(t)=[R(t),r] with coordinates (r, t), moving wave emitting boundary R(t), and fixed boundary r, and the fixed computational domain Θ=[XR,X] with coordinates (ξ(r,t),t) and fixed boundaries XR and X. Next, we invoke the time-dependent coordinate transformation r:[XR,X][R(t),r],(ξ,t)r(ξ,t) to implement the generic change of variables,

ϕ1(r,t)=Φ1(ξ(r,t),t),
(28)
u0(r,t)=U0(ξ(r,t),t),
(29)
ψ(r,t)=Ψ(ξ(r,t),t).
(30)

A similar approach for a wave emitter moving relative to a quiescent background medium is presented by Gasperini et al.1 The derivative operators conveying the change of variables follow as

ddr=ξrξ,
(31)
ddt=ξtξ+t.
(32)

Applying the derivative operators to Eqs. (28) and (29) and substituting into Eq. (27), the transformed wave equation,

ξr[2U02ξξt+U0t+(ξt+U0ξr)U0ξ]Φ1ξ+ξr[ξr(Ψξ+U0U0ξ)2rc2]Φ1ξ[(c2U02)ξr2ξrξ2ξt2ξt2ξtξ]Φ1ξ[(c2U02)(ξr)22U0ξtξr(ξt)2]2Φ1ξ2+2(ξt+U0ξr)2Φ1ξt=2Φ1t2,
(33)

is obtained, where the spherical Laplacian gives rise to the additional contribution (2/r)Φ1/ξ. Only the derivatives of Φ1 need to be discretized. The derivatives of ξ, ψ, and U0 depend on the specification of the coordinate transformation and the background flow field and can, therefore, be seen as space- and time-dependent coefficients.

In particular, the coordinate transformation can be chosen in such a way that a non-negative coefficient of the Laplacian 2Φ1/ξ2 is obtained, thereby guaranteeing the well-posedness of the problem. It follows from Eq. (33) that the Laplacian coefficient is non-negative if the condition

c|U0+ξt(ξr)1|
(34)

is satisfied. The term ξ/t is the velocity of a fixed point in Θ as viewed from a moving point in Ω(t).

We opt for a linear mapping between the moving physical domain Ω(t) and the fixed computational domain Θ as described by the linear function

ξ(r,t)=X+(rr)XXRrR(t),
(35)

and choose XR=0 and X=1. One is, in principle, free in the choice of XR and X. However, we will choose the sonic horizon radius on the order of unity and scale the emitted wavelength accordingly, so that the Jacobian of the coordinate transformation ξ/r is on the order of unity with the above choices of XR and X. The derivatives of ξ follow as

ξr=XXRrR(t),
(36)
ξt=XξrR(t)Ṙ(t),
(37)
2ξtξ=Ṙ(t)rR(t),
(38)
2ξrξ=0,
(39)
2ξt2=2ξt2ξtξR¨ξrXξXXR.
(40)

From the above relations and Eqs. (1) and (25), the derivatives of ψ and U follow as

Ψξ=2nc4r2gh(ξr)1,
(41)
U0ξ=nU0r(ξ)(ξr)1,
(42)
U0t=(Rr)n[R¨+nR(Ṙ)2+nṘrξt(ξr)1].
(43)

Finally, R, Ṙ, and R¨ need to be specified. For a stationary spherical flow as described by Eq. (1), these quantities follow by solving Eq. (1) for R under the boundary condition u(rh)=const=c (steady-state flow) and the initial condition R(t=0)=R0, which gives

R(t)=[R0n+1(n+1)crhnt]1n+1,
(44)
Ṙ(t)=crhn[R0n+1(n+1)crhnt]nn+1,
(45)
R¨(t)=nc2rh2n[R0n+1(n+1)crhnt]2n+1n+1,
(46)

where it is assumed that n=const.

An explicit finite difference time domain (FDTD) method based on central differences in space is used to solve the transformed wave equation [Eq. (33)] in the fixed computational domain Θ. Without special treatment, the explicit FDTD solution may be contaminated by dispersive numerical noise,38,39 which can be amplified by the frequency side-bands that are generated by the accelerating background flow field. In order to counteract the development and amplification of dispersive numerical noise, we opt for sixth-order accurate central differences in space, where this relatively high order helps to delay the onset of dispersive numerical noise.40 The ξ-discrete solution is explicitly advanced in time using the predictor–corrector method by Dey and Dey.41 Originally developed for ordinary differential equations, the predictor–corrector method was later applied to the linear wave equation in the context of seismic wave propagation by Nascimento and Pestana,42 who showed that the method effectively mimics an attenuation term weighted by a factor γ. This term predominantly acts on the higher wave frequencies so that the diffusive effect on the numerical solution vanishes with increasing grid resolution while preserving the stabilizing effect. In the present work, we use the recommended value γ=0.5, for which the predictor–corrector method evolves into a second-order Runge–Kutta scheme when applied to an ordinary differential equation.43 However, on application to Westervelt's wave equation,44 we previously found a convergence rate between first and second orders, which is mainly attributed to the treatment of the mixed derivatives as source terms as described in our previous work,3 in which a more detailed description of the numerical method and rigorous convergence studies are provided. For a detailed discussion of the stability properties of the predictor–corrector method, we refer to the work of Dey.43 

We consider a spherically symmetric steady-state acoustic black hole for a Bondi-type (inward directed) flow as described by Eq. (27), with the gradient of the gravitational potential given by Eq. (25). The input flow field is given by Eq. (1), where the boundary motion is specified by Eqs. (44) and (45). We choose n =2 (incompressible flow field) for our base configuration, but we also present results for n <2.

The acoustic waves are emitted from a virtual spherical boundary with instantaneous radius R(t) collapsing along with the spherically symmetric fluid flow (see Fig. 1). We arbitrary choose a sound speed of c=1500m/s. Given the dispersion relation c=λafa, the values of He [see Eq. (24)] and gh [see Eq. (4) for c/r=0] uniquely define the emitted wave frequency fa and the background flow velocity distribution. The simulations comprise a dimensionless time of fat=27.5 emitted wave periods. The sonic horizon position rh is implicitly specified by Eq. (1), where u0(rh)=c. The initial wave emitter position R0>rh (the wave emitter is “falling” into the acoustic black hole) is chosen in such a way that the sonic horizon approximately intersects the 14th wave period for a slowly varying background velocity u0. The position r of the far field boundary is chosen in such a way that the waves do not reach the boundary. Apart from the grid convergence study in Sec. IV B, the grid resolution is Nppw=400, where Nppw is the number of grid points per emitted wavelength λa. The Courant–Friedrich–Lewy number is always CFL=cΔt/Δr=0.1, where Δt is the time step size and Δr the spatial step size as specified by Nppw and λa.

Figure 2 depicts an instantaneous acoustic waveform intersected by the sonic horizon at r/rh=1 for five different grid resolutions. The acoustic perturbation potential ϕ1 is normalized by the excitation amplitude Δϕ1,a. The acoustic surface gravity is gh=7.5×106m/s, and the Helmholtz number as given by Eq. (24) is He=20. This configuration involves the largest surface gravity and the smallest Helmholtz number considered in this study. Hence, it involves the largest gradients of the background flow velocity relative to the acoustic wavelength and constitutes, therefore, the most challenging case for the numerical method. The waveform in Fig. 2 is recorded at fat=27.5. The figure on the right shows a detail of the wave profile close to the outer domain boundary r, which corresponds to the longest travel distance of the wave. It can be seen that for Nppw=400, a sufficiently well-converged solution is obtained. This required resolution agrees with the findings of our previous work,3 in which a more rigorous analysis of the accuracy of the numerical method is presented.

FIG. 2.

Instantaneous acoustic wave profile intersected by the sonic horizon at r/rh=1 (left) with a detailed view (right) recorded at fat=27.5 for gh=7.5×106m/s [see Eq. (4) for c/r=0] and He=20 [see Eq. (24)]. The acoustic perturbation potential ϕ1 is normalized by the excitation amplitude Δϕ1,a. The wave profile is depicted for five different grid resolutions, where Nppw is the number of grid points per emitted wavelength λa.

FIG. 2.

Instantaneous acoustic wave profile intersected by the sonic horizon at r/rh=1 (left) with a detailed view (right) recorded at fat=27.5 for gh=7.5×106m/s [see Eq. (4) for c/r=0] and He=20 [see Eq. (24)]. The acoustic perturbation potential ϕ1 is normalized by the excitation amplitude Δϕ1,a. The wave profile is depicted for five different grid resolutions, where Nppw is the number of grid points per emitted wavelength λa.

Close modal

Figure 3 shows the distribution of the acoustic perturbation potential ϕ1 in the space–time plane. The spatial axis is normalized by the sonic horizon radius rh and the time axis by the emitted wave frequency fa. The distribution is depicted for different values of the surface gravity gh [see Eq. (4) for c/r=0] and the Helmholtz number He [see Eq. (24)]. Starting from subfigure (a), which represents the smallest acoustic black hole relative to the fixed emitted wavelength, the distributions in (b)–(d) are obtained by increasing the sonic horizon radius rh while leaving the other parameters unchanged, so that gh decreases and He increases accordingly. Initially, the wave emitting boundary with instantaneous radius R(t) is located outside the sonic horizon as indicated by the black solid line. Subsequently, the wave emitting boundary collapses and eventually passes the sonic horizon. It can be seen that the wave characteristics diverge away from the sonic horizon over time, which corresponds to a gradual red-shift of the outgoing wave train.

FIG. 3.

Distribution of the acoustic perturbation potential ϕ1, normalized by the excitation amplitude Δϕ1,a, in the space–time plane for different values of the surface gravity gh [see Eq. (4) for c/r=0] and the Helmholtz number He [see Eq. (24)]. The different combinations in the subfigures (a)–(d) are obtained by varying the sonic horizon radius rh for a fixed emitted wavelength λa. A decreasing value of the acoustic black hole radius rh is associated with an increasing magnitude of the surface gravity, associated with faster divergence of neighboring wave characteristics as well as a more pronounced curvature of both the wave emitter trajectory and the wave characteristics. (a) He = 20; gh = 7.500× 106 m/s2. (b) He = 40; gh = 3.750× 106 m/s2. (c) He = 60; gh = 2.500× 106 m/s2. (d) He = 80; gh = 1.875× 106 m/s2.

FIG. 3.

Distribution of the acoustic perturbation potential ϕ1, normalized by the excitation amplitude Δϕ1,a, in the space–time plane for different values of the surface gravity gh [see Eq. (4) for c/r=0] and the Helmholtz number He [see Eq. (24)]. The different combinations in the subfigures (a)–(d) are obtained by varying the sonic horizon radius rh for a fixed emitted wavelength λa. A decreasing value of the acoustic black hole radius rh is associated with an increasing magnitude of the surface gravity, associated with faster divergence of neighboring wave characteristics as well as a more pronounced curvature of both the wave emitter trajectory and the wave characteristics. (a) He = 20; gh = 7.500× 106 m/s2. (b) He = 40; gh = 3.750× 106 m/s2. (c) He = 60; gh = 2.500× 106 m/s2. (d) He = 80; gh = 1.875× 106 m/s2.

Close modal

Figure 4 shows the instantaneous wave profiles corresponding to the configurations in Fig. 3, recorded at fat=27.5. For reference, the black dashed line indicates a 1/r-decay of the waveform, starting from the sonic horizon. The reference value for the wave amplitude at the sonic horizon is obtained by linear interpolation between the two neighboring wave peaks. The general observation is that the entire wave is increasingly amplified with increasing magnitude of gh. In all cases, the rate of amplification is most pronounced in the innermost region close to the wave emitting boundary so that the maximum wave amplitude is observed at some distance away from the location of wave emission, in contrast to the 1/r-decay dictated by the spherically symmetric geometry. It appears that the dimensionless distance r/rh of the maximum wave amplitude increases with a decreasing magnitude of gh.

FIG. 4.

Instantaneous wave profiles corresponding to the space–time diagrams in Fig. 3 recorded at fat=27.5. In the subfigures (a)–(d), the sonic horizon radius rh is varied systematically while leaving the emitted wavelength λa constant, resulting in different values of the Helmholtz number He [see Eq. (24)] and the surface gravity gh [see Eq. (4) for c/r=0]. For reference, the black dashed line represents the 1/r-decay of the waveform, starting from the sonic horizon. (a) He = 20; gh = 7.500× 106 m/s2. (b) He = 40; gh = 3.750× 106 m/s2. (c) He = 60; gh = 2.500× 106 m/s2. (d) He = 80; gh = 1.875× 106 m/s2.

FIG. 4.

Instantaneous wave profiles corresponding to the space–time diagrams in Fig. 3 recorded at fat=27.5. In the subfigures (a)–(d), the sonic horizon radius rh is varied systematically while leaving the emitted wavelength λa constant, resulting in different values of the Helmholtz number He [see Eq. (24)] and the surface gravity gh [see Eq. (4) for c/r=0]. For reference, the black dashed line represents the 1/r-decay of the waveform, starting from the sonic horizon. (a) He = 20; gh = 7.500× 106 m/s2. (b) He = 40; gh = 3.750× 106 m/s2. (c) He = 60; gh = 2.500× 106 m/s2. (d) He = 80; gh = 1.875× 106 m/s2.

Close modal

Figure 5 depicts the wave amplitude modulation at the sonic horizon over time, where the amplitude modulation factor is defined as AM(rh,t)=ϕ1(rh,tth)/ϕ1(rh,th) for tth. The time th marks the time instant at which the wave emitting boundary with instantaneous position R(t) passes the sonic horizon at rh. It can be seen that the growth rate of AM increases with gh. Again, this illustrates the overall amplification of the waveform imparted by the acceleration of the background flow. Even though we are lacking an analytical expression for AM at this time, it is interesting to note that the curvatures of the AM lines appear to be qualitatively related to the curvatures of the corresponding wave emitter trajectories in Fig. 3.

FIG. 5.

Amplitude modulation factor AM(rh,t)=ϕ1(rh,tth)/ϕ1(rh,th) for tth evaluated at the sonic horizon, where th marks the time instant at which the wave emitting boundary (instantaneous position R(t)) passes the sonic horizon at rh. The different combinations correspond to Figs. 3 and 4, where the sonic horizon radius rh is varied systematically while leaving the emitted wavelength λa constant, resulting in the different values of the Helmholtz number He [see Eq. (24)] and the surface gravity gh [see Eq. (4) for c/r=0].

FIG. 5.

Amplitude modulation factor AM(rh,t)=ϕ1(rh,tth)/ϕ1(rh,th) for tth evaluated at the sonic horizon, where th marks the time instant at which the wave emitting boundary (instantaneous position R(t)) passes the sonic horizon at rh. The different combinations correspond to Figs. 3 and 4, where the sonic horizon radius rh is varied systematically while leaving the emitted wavelength λa constant, resulting in the different values of the Helmholtz number He [see Eq. (24)] and the surface gravity gh [see Eq. (4) for c/r=0].

Close modal

For larger values of gh in Fig. 3, the wave characteristics possess a distinct curvature, whereas they flatten out as gh decreases. At the same time, the characteristics tend to evolve more and more parallel to the sonic horizon line, which means that for very small gh and very large He, the entire wave train degenerates into a standing wave as viewed from the reference frame of the sonic horizon. This is demonstrated by the space–time diagram of ϕ1 in Fig. 6(a), in which the sonic horizon rh is further increased while again leaving the other parameters unchanged. Figure 6(b) depicts the corresponding instantaneous wave profile at fat=27.5, indicating that the waveform becomes planar and exhibits increasingly linear behavior as rh increases. The increasingly linear behavior manifests in a reduction of the frequency side-bands and amplitude modulations.

FIG. 6.

Distribution of the acoustic perturbation potential ϕ1, normalized by the excitation amplitude Δϕ1,a, in the space–time plane (a) and the corresponding instantaneous wave profile recorded at fat=27.5 (b) for a large sonic horizon radius rh as compared to the emitted wavelength λa, where He=104/3 and gh=4.500×104m/s2. (a) Space–time diagram. (b) Instantaneous wave profile.

FIG. 6.

Distribution of the acoustic perturbation potential ϕ1, normalized by the excitation amplitude Δϕ1,a, in the space–time plane (a) and the corresponding instantaneous wave profile recorded at fat=27.5 (b) for a large sonic horizon radius rh as compared to the emitted wavelength λa, where He=104/3 and gh=4.500×104m/s2. (a) Space–time diagram. (b) Instantaneous wave profile.

Close modal

As discussed in Sec. II A, the same Helmholtz number He [see Eq. (24)] can be obtained for different combinations of the emitted wavelength λa and the surface gravity gh. Figure 7 shows such a variation for the smallest (a) and the largest (d) acoustic black hole in Fig. 3, which are represented by subfigures (a) and (b) in Fig. 7. Figures 7(c) and 7(d) are based on the same Helmholtz numbers as (a) and (b), respectively. However, the surface gravity in (c) is equal to the surface gravity in (b), and the surface gravity in (d) is equal to the surface gravity in (a), which means that in both (c) and (d), the emitted wavelength is adjusted to obtain the corresponding values for He. It follows from Eq. (24) that the required scaling factor of λa is the inverse of the surface gravity ratios, which is 4 in this particular case. The emitted wave frequency changes according to fa=c/λa. The comparison between Figs. 7(a) and 7(c) and Figs. 7(b) and 7(d), respectively, shows that for identical values of He, the wave characteristics evolve along virtually identical trajectories, however in a space–time plane normalized by rh and fa. This means that the characteristics corresponding to the same Helmholtz number He for different combinations of gh and λa are not identical in physical space–time, albeit geometrically similar.

FIG. 7.

Space–time diagram of the acoustic perturbation potential ϕ1, normalized by the excitation amplitude Δϕ1,a, for a variation of gh and λa to obtain geometrically similar wave characteristics for different Helmholtz numbers He as given by Eq. (24). Subfigures (c) and (d) are based on the same Helmholtz numbers as (a) and (b), respectively, but the surface gravity in (c) is equal to the surface gravity in (b), and the surface gravity in (d) is equal to the surface gravity in (a). In (c) and (d), λa is adjusted by multiplying/diving by a factor of 4 as compared to (a) and (b) to obtain the corresponding values for He. (a) He = 20; gh = 7.500× 106 m/s2. (b) He = 80; gh = 1.875× 106 m/s2. (c) He = 20; gh = 1.875× 106 m/s2. (d) He = 80; gh = 7.500× 106 m/s2.

FIG. 7.

Space–time diagram of the acoustic perturbation potential ϕ1, normalized by the excitation amplitude Δϕ1,a, for a variation of gh and λa to obtain geometrically similar wave characteristics for different Helmholtz numbers He as given by Eq. (24). Subfigures (c) and (d) are based on the same Helmholtz numbers as (a) and (b), respectively, but the surface gravity in (c) is equal to the surface gravity in (b), and the surface gravity in (d) is equal to the surface gravity in (a). In (c) and (d), λa is adjusted by multiplying/diving by a factor of 4 as compared to (a) and (b) to obtain the corresponding values for He. (a) He = 20; gh = 7.500× 106 m/s2. (b) He = 80; gh = 1.875× 106 m/s2. (c) He = 20; gh = 1.875× 106 m/s2. (d) He = 80; gh = 7.500× 106 m/s2.

Close modal

However, as shown by the instantaneous wave profiles in Fig. 8, the radial distribution of the acoustic perturbation ϕ1 is different for identical values of He but different values of gh. While the difference is negligible for He=80, it becomes more pronounced for He=20. Furthermore, Fig. 8 echoes our previous observations; a larger surface gravity induces a stronger amplification of the waveform and is accompanied by a shift toward closer radii r/rh of the wave amplitude maximum.

FIG. 8.

Instantaneous wave profiles corresponding to the space–time diagrams in Fig. 7 recorded at fat=27.5. Subfigures (c) and (d) are based on the same Helmholtz numbers as (a) and (b), respectively, but the surface gravity in (c) is equal to the surface gravity in (b), and the surface gravity in (d) is equal to the surface gravity in (a). (a) He = 20; gh = 7.500× 106 m/s2. (b) He = 80; gh = 1.875× 106 m/s2. (c) He = 20; gh = 1.875× 106 m/s2. (d) He = 80; gh = 7.500× 106 m/s2.

FIG. 8.

Instantaneous wave profiles corresponding to the space–time diagrams in Fig. 7 recorded at fat=27.5. Subfigures (c) and (d) are based on the same Helmholtz numbers as (a) and (b), respectively, but the surface gravity in (c) is equal to the surface gravity in (b), and the surface gravity in (d) is equal to the surface gravity in (a). (a) He = 20; gh = 7.500× 106 m/s2. (b) He = 80; gh = 1.875× 106 m/s2. (c) He = 20; gh = 1.875× 106 m/s2. (d) He = 80; gh = 7.500× 106 m/s2.

Close modal

In order to investigate the effect of the flow field scaling parameter n in Eq. (1), we start from the case in Fig. 3(a), where n =2, and reduce the scaling parameter to n =0.5 while leaving the remaining parameters unchanged. The magnitude of the radial flow velocity then increases toward the center with r1/2 [see Eq. (1)], which corresponds to the free-fall condition in a Bondi-type flow.45Figure 9 shows both the space–time distributions of ϕ1 and the instantaneous wave profiles at fat=27.5 for both configurations. A comparison of (a) and (c) shows that a reduction of n reduces the curvature of the wave characteristics, which is due to the reduction of the background flow velocity gradients. Equation (4) implies that the surface gravity gh reduces linearly with decreasing n. A comparison of (b) and (d) shows that the smaller magnitude of gh for n =0.5 manifests in a less pronounced amplification of the waveform, which is in agreement with the previous observations.

FIG. 9.

Space–time diagrams of ϕ1/Δϕ1,a and corresponding instantaneous wave profiles recorded at fat=27.5 for a variation of the background velocity scaling factor n in Eq. (1) to mimic the effect of a compressible background fluid. Starting from the base configuration in (a) and (b), n is reduced in (c) and (d) while leaving the remaining parameters unchanged, where {He=20;gh=7.500×106m/s2} for n =2 and {He=80;gh=1.875×106m/s2} for n =0.5. (a) Space–time diagram of ϕ1∕Δϕ1,a for n = 2. (b) Wave profile for n = 2. (c) Space–time diagram of ϕ1∕Δϕ1,a for n = 0.5. (d) Wave profile for n = 0.5.

FIG. 9.

Space–time diagrams of ϕ1/Δϕ1,a and corresponding instantaneous wave profiles recorded at fat=27.5 for a variation of the background velocity scaling factor n in Eq. (1) to mimic the effect of a compressible background fluid. Starting from the base configuration in (a) and (b), n is reduced in (c) and (d) while leaving the remaining parameters unchanged, where {He=20;gh=7.500×106m/s2} for n =2 and {He=80;gh=1.875×106m/s2} for n =0.5. (a) Space–time diagram of ϕ1∕Δϕ1,a for n = 2. (b) Wave profile for n = 2. (c) Space–time diagram of ϕ1∕Δϕ1,a for n = 0.5. (d) Wave profile for n = 0.5.

Close modal

Finally, it is noted that we do not observe geometrical similarity for identical Helmholtz numbers and different values of the flow field parameter n. This can be seen by comparing Fig. 9(c) with Fig. 3(d). Both figures are based on the same Helmholtz number and also the same surface gravity, but the latter exhibits a significantly smaller expansion of the wave characteristics around the sonic horizon region. Hence, it appears that geometrical similarity of the wave characteristics strongly relies on a given background flow field.

The acoustic black hole analogy employed to analyze acoustic waves propagating in non-uniform background flow fields reveals several interesting features. Most importantly, the acoustic black hole analogy allows for a parametrization of the wave-flow configuration in terms of a Helmholtz number He that is expressed as a function of the acoustic black hole's surface gravity [see Eq. (24)]. The numerical simulations demonstrate that He describes geometrically similar wave characteristics, which can, however, be subject to different amplitude modulations depending on the magnitude of the surface gravity. In particular, an increase in the surface gravity and, hence, the overall magnitude of the background flow acceleration cause an increasing amplification of the entire waveform, which again becomes increasingly noticeable over time as the wave propagates through the accelerating flow field. Consequently, the maximum wave amplitude occurs at some distance away from the wave emitting boundary. A qualitatively similar effect is reported by Christov and Christov2 for a wave emitter accelerating through a quiescent background medium in a one-dimensional Cartesian case, where the nonlinear Doppler amplification becomes more noticeable with increasing distance from the source. In principle, the above observations are expected to apply to subsonic flows in a very similar fashion as the amplitude modulations merely depend on the acceleration of the flow field and not on the magnitude of the flow velocity. In fact, Fig. 6 shows that an acoustic wave exhibits quasi-linear behavior across the sonic horizon if the surface gravity is small enough and He is large. Loosely speaking, this scenario is associated with very large acoustic black holes, having in mind that the acoustic black hole's size should be assessed in conjunction with the acoustic wavelength.

In addition to the merits regarding the study of nonlinear wave dynamics in accelerating flow fields, the acoustic black hole analogy may also prove useful for the verification and benchmarking of numerical methods for acoustic waves in moving fluids. On the one hand, the system is straightforward to describe geometrically, with the sonic horizon as a distinct reference point, and it is well defined in terms of a small number of parameters. On the other hand, it is rich in phenomena and complexity, which can be increased arbitrary by considering a non-stationary sonic horizon, circumferential flow components, and other three-dimensional effects, nonlinear wave steepening, or discontinuities in the background flow field.

A general advantage of the presented model formulation is that the effects of a compressible background medium are absorbed in terms that exclusively depend on the stationary background flow velocity [see Eq. (27)]. Hence, the problem reduces to the knowledge of the appropriate background flow field that satisfies the continuity equation for a compressible flow of a barotropic fluid. In particular, we can study the principle effect of a space-dependent distribution of ρ0 by varying the background velocity field using the parameter n in Eq. (1). While n =2 corresponds to the incompressible flow assumption, n <2 is expected in the region of sonic transition as it can be inferred from Eq. (3) and the works of Treves et al.45 and Ray and Bhattacharjee.31 As shown by Fig. 9(a), a decrease in n reduces the curvature of the wave characteristics as a result of the smaller velocity gradients. Equation (1) implies that a reduction of n moves the sonic horizon toward larger radii for given values of R and Ṙ. The reduced curvature of the wave characteristics manifests in a less pronounced amplification of the waveform [see Fig. 9(b)]. We may further anticipate the qualitative effect of a variable speed of sound, which is a second-order effect as discussed in Sec. II C. Due to the compression of the fluid toward the center point, the speed of sound should increase toward the center accordingly. Taking the speed of sound at the sonic horizon radius rh as a reference, c increases toward the inner region and decreases toward the outer region. For an outward propagating wave, this has the qualitative effect that the characteristics in the inner region evolve slower toward the center point, whereas the characteristics in the outer region escape at a lower speed. Therefore, the variable speed of sound is expected to straighten the wave characteristics, which is qualitatively similar to the effect of the variable density field. It is further noted that a variable speed of sound also alters the surface gravity [see Eq. (4)].

The transformation of the governing wave equation for a moving background medium conveyed by a generic change of variables is shown to be instrumental to maintain the well-posedness of the governing wave equation in the presence of a transonic flow field. In the present work, we opt for a time-dependent linear transformation in physical space. The linearity of the transformation facilitates the specification of the involved derivatives as a function of the boundary motion and the motion of the background flow field [see Eqs. (36)–(43)]. However, the linearity of the transformation also imposes a constraint on the possible combinations of the acoustic black hole size and the size of the moving/deforming physical domain Ω(t). This is because the speed of any point in Ω(t) relative to the moving background flow must remain below the sonic speed in order to preserve the well-posedness of the problem as discussed in Sec. III A. Drawing on the method of characteristics, this constraint can be relaxed by specifying the transformation in such a way that any point of the moving physical domain Ω(t) moves exactly with the background fluid. For a uniform background flow field, the convective wave equation [Eq. (27) with zero right-hand side] would then reduce to the standard wave equation under this particular Galilean transformation as it readily follows from the transformed wave equation [Eq. (33)]. Any potential non-uniformity of the flow field generates additional terms, which are then exclusively related to nonlinear Doppler modulations, that is, the generation of frequency side-bands and amplitude modulations. Thus, the appropriate coordinate transformation is expected to greatly facilitate the study of Doppler shifts and nonlinear amplitude modulations in complicated flow situations, as has also been concluded by Gregory et al.22 For instance, the method of characteristics may permit full analytical solutions for certain flow situations, in which separations of length and/or time scales can be identified. An example is found in the work by Christov and Christov,2 where an asymptotic analytical solution for the wave amplitude and frequency modulations caused by an oscillating wave emitting boundary in a quiescent background medium is given, under the assumption that the timescale of the boundary motion is large as compared to the emitted wave period. While Fig. 6 demonstrates linear wave behavior in the limit of a very small surface gravity, a more detailed analysis is needed to understand the asymptotic behavior for very large magnitudes of the surface gravity.

Furthermore, the change of variables may be extended to the time variable, which would admit Lorentz transformations as described by Wang.23 While the transformation in physical space is sufficient as long as the wave emitter moves along with the background fluid flow as it is the case in the present work, Lorentz transformations may indeed help to preserve the well-posedness of the governing wave equation in the case of relative motion between the wave emitter and the background fluid. Finally, the coordinate transformation may also be harnessed to redistribute grid points in physical space such that sharp features of the waveforms are accurately resolved. In such a moving grid method, a fixed, uniform grid is implemented in transformed space and the coordinate transformation governing the grid point motion is linked to a measure of targeted accuracy.46,47

Within the scope of an acoustic black hole analogy, the acoustic wave-flow configuration can be parameterized in terms of the Helmholtz number He=c2/(ghλa), where c, gh, and λa are the speed of sound, the acoustic black hole's surface gravity,9 and the emitted wavelength, respectively. The Helmholtz number in this particular form is shown to describe geometrically similar wave characteristics for different combinations of gh and λa. Its applicability is not restricted to spherically symmetric flows, but extends to any quasi one-dimensional geometry. However, the wave amplitude modulations along geometrically similar wave characteristics depend on the magnitude of gh. While the acoustic black hole analogy is a useful means to characterize the acoustic wavelength relative to the acceleration of a stationary but non-uniform background flow field, the observations regarding the nonlinear wave amplitude modulations are expected to apply to subsonic flows in a very similar fashion as they merely depend on the flow acceleration. For small gh, the linear wave behavior is recovered in the proximity of the sonic horizon. Therefore, the present work promotes the acoustic black hole analogy for the study of nonlinear acoustic wave dynamics and Doppler modulations in accelerating flow fields. Regarding the computational method, it is shown that a suitable Galilean-type coordinate transformation from a moving/deforming physical domain to a fixed computational domain can maintain the well-posedness of the governing wave equation in the presence of a transonic background flow. Following the suggestions by Gregory et al.22 and Wang,23 the flexibility afforded by the coordinate transformation can be further exploited in order to isolate terms that are exclusively related to nonlinear Doppler effects from those terms that constitute the invariants of the wave equation. This is expected to greatly simplify the analysis of nonlinear Doppler modulations in moving flow fields.22 

This research was funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation), Grant Nos. 441063377 and 443546539.

The authors have no conflicts to disclose.

Soeren Schenke: Conceptualization (lead); Data curation (lead); Formal analysis (lead); Investigation (lead); Methodology (lead); Software (lead); Validation (lead); Visualization (lead); Writing – original draft (lead). Fabian Sewerin: Formal analysis (supporting); Funding acquisition (supporting); Methodology (supporting); Writing – original draft (supporting); Writing – review & editing (equal). Berend van Wachem: Funding acquisition (equal); Project administration (equal); Resources (equal); Supervision (supporting); Writing – review & editing (equal). Fabian Denner: Funding acquisition (equal); Methodology (supporting); Project administration (equal); Resources (equal); Software (supporting); Supervision (lead); Writing – review & editing (equal).

The data that support the findings of this study are available from the corresponding author upon reasonable request.

1.
D.
Gasperini
,
H.-P. P.
Beise
,
U.
Schroeder
,
X.
Antoine
, and
C.
Geuzaine
, “
A frequency domain method for scattering problems with moving boundaries
,”
Wave Motion
102
,
102717
(
2021
).
2.
I. C.
Christov
and
C. I.
Christov
, “
On mechanical waves and Doppler shifts from moving boundaries
,”
Math. Methods Appl. Sci.
40
,
4481
4492
(
2017
).
3.
S.
Schenke
,
F.
Sewerin
,
B.
van Wachem
, and
F.
Denner
, “
Explicit predictor-corrector method for nonlinear acoustic waves excited by a moving wave emitting boundary
,”
J. Sound Vib.
527
,
116814
(
2022
).
4.
A.
Heimdal
and
H.
Torp
, “
Ultrasound Doppler measurements of low velocity blood flow: Limitations due to clutter signals from vibrating muscles
,”
IEEE Trans. Ultrason. Ferroelectr. Freq. Control
44
,
873
881
(
1997
).
5.
A. A.
Oglat
,
M. Z.
Matjafri
,
N.
Suardi
,
M. A.
Oqlat
,
M. A.
Abdelrahman
, and
A. A.
Oqlat
, “
A review of medical Doppler ultrasonography of blood flow in general and especially in common carotid artery
,”
J. Med. Ultrasound
26
,
3–13
(
2018
).
6.
A. H.
Khandoker
,
F.
Marzbanrad
, and
Y.
Kimura
, “
Editorial: Recent advances in Doppler signal processing and modeling techniques for fetal monitoring
,”
Front. Physiology
9
,
691
(
2018
).
7.
Q.-B.
Wang
and
X.-H.
Ge
, “
Geometry outside of acoustic black holes in (2 + 1)-dimensional spacetime
,”
Phys. Rev. D
102
,
104009
(
2020
).
8.
W. G.
Unruh
, “
Experimental black-hole evaporation?
,”
Phys. Rev. Lett.
46
,
1351
1353
(
1981
).
9.
M.
Visser
, “
Acoustic black holes: Horizons, ergospheres and Hawking radiation
,”
Classical Quantum Gravity
15
,
1767
1791
(
1998
).
10.
S. W.
Hawking
, “
Particle creation by black holes
,”
Commun. Math. Phys.
43
,
199
220
(
1975
).
11.
H.
Furuhashi
,
Y.
Nambu
, and
H.
Saida
, “
Simulation of an acoustic black hole in a Laval nozzle
,”
Classical Quantum Gravity
23
,
5417
5438
(
2006
).
12.
P.
Jain
,
A. S.
Bradley
, and
C. W.
Gardiner
, “
Quantum de Laval nozzle: Stability and quantum dynamics of sonic horizons in a toroidally trapped Bose gas containing a superflow
,”
Phys. Rev. A
76
,
023617
(
2007
).
13.
R.
da Rocha
, “
Black hole acoustics in the minimal geometric deformation of a de Laval nozzle
,”
Eur. Phys. J. C
77
,
355
(
2017
).
14.
V.
Cardoso
,
A.
Coutant
,
M.
Richartz
, and
S.
Weinfurtner
, “
Detecting rotational superradiance in fluid laboratories
,”
Phys. Rev. Lett.
117
,
271101
(
2016
).
15.
T.
Torres
,
S.
Patrick
,
A.
Coutant
,
M.
Richartz
,
E. W.
Tedford
, and
S.
Weinfurtner
, “
Rotational superradiant scattering in a vortex flow
,”
Nat. Phys.
13
,
833
(
2017
).
16.
B.
Demirkaya
,
T.
Dereli
, and
K.
Güven
, “
Acoustic superradiance from a Bose–Einstein condensate vortex with a self-consistent background density profile
,”
Phys. Scr.
95
,
055001
(
2020
).
17.
S.
Patrick
,
H.
Goodhew
,
C.
Gooding
, and
S.
Weinfurtner
, “
Backreaction in an analogue black hole experiment
,”
Phys. Rev. Lett.
126
,
041105
(
2021
).
18.
R.
Balbinot
,
S.
Fagnocchi
, and
A.
Fabbri
, “
Quantum effects in acoustic black holes: The backreaction
,”
Phys. Rev. D
71
,
064019
(
2005
).
19.
R. J. R.
Williams
and
J. E.
Dyson
, “
On sonic transitions in astrophysical flows
,”
Mon. Not. R. Astron. Soc.
270
,
L52
L56
(
1994
).
20.
H.
Bondi
, “
On spherically symmetrical accretion
,”
Mon. Not. R. Astron. Soc.
112
,
195
(
1952
).
21.
E. N.
Parker
, “
Dynamics of the interplanetary gas and magnetic fields
,”
Astrophys. J.
128
,
664
(
1958
).
22.
A.
Gregory
,
S.
Sinayoko
,
A.
Agarwal
, and
J.
Lasenby
, “
An acoustic space-time and the Lorentz transformation in aeroacoustics
,”
Int. J. Aeroacoust.
14
,
977
1003
(
2015
).
23.
S.
Wang
, “
Extensions to the Navier–Stokes equations
,”
Phys. Fluids
34
,
053106
(
2022
).
24.
T. L.
Chow
, “
The physics of black holes
,” in
Gravity, Black Holes, and the Very Early Universe: An Introduction to General Relativity and Cosmology
(
Springer
,
New York
,
2008
), pp.
81
109
.
25.
D.
Long
and
R.
Arndt
, “
The role of Helmholtz number in jet noise
,” in
AIAA 22nd Aerospace Sciences Meeting
(
AIAA
,
1984
), p.
403
.
26.
R.
Ewert
and
W.
Schröder
, “
Acoustic perturbation equations based on flow decomposition via source filtering
,”
J. Comput. Phys.
188
,
365
398
(
2003
).
27.
O.
Godin
, “
An exact wave equation for sound in inhomogeneous, moving, and non-stationary fluids
,”
OCEANS'11—MTS/IEEE Kona, Program Book
(
IEEE
,
2011
), pp.
1
5
.
28.
E.
Keto
, “
Stability and solution of the time-dependent Bondi–Parker flow
,”
Mon. Not. R. Astron. Soc.
493
,
2834
2840
(
2020
).
29.
A.
Aguayo-Ortiz
,
E.
Tejeda
,
O.
Sarbach
, and
D.
López-Cámara
, “
Spherical accretion: Bondi, Michel, and rotating black holes
,”
Mon. Not. R. Astron. Soc.
504
,
5039
5053
(
2021
).
30.
S. K.
Chakrabarti
, “
Accretion processes on a black hole
,”
Phys. Rep.
266
,
229
390
(
1996
).
31.
A.
Ray
and
J.
Bhattacharjee
, “
Static and dynamic aspects of transonicity in Bondi accretion
,”
Indian J. Phys.
80
,
1123
(
2007
).
32.
I. G.
Kovalenko
and
M. A.
Eremin
, “
Instability of spherical accretion—I. Shock-free Bondi accretion
,”
Mon. Not. R. Astron. Soc.
298
,
861
870
(
1998
).
33.
B. A.
Maicke
,
J.
Majdalani
, and
R. L.
Geisler
, “
Characterization of the startup and pressure blowdown processes in rocket nozzles
,”
Aerosp. Sci. Technol.
25
,
273
282
(
2013
).
34.
P.-H.
Chavanis
and
T.
Harko
, “
Bose-Einstein condensate general relativistic stars
,”
Phys. Rev. D
86
,
064011
(
2012
).
35.
C. G.
Böhmer
and
T.
Harko
, “
Can dark matter be a Bose–Einstein condensate?
,”
J. Cosmol. Astropart. Phys.
2007
,
025
.
36.
M.
Anacleto
,
F.
Brito
, and
E.
Passos
, “
Superresonance effect from a rotating acoustic black hole and Lorentz symmetry breaking
,”
Phys. Lett. B
703
,
609
613
(
2011
).
37.
B.
Zhang
and
E. N.
Saridakis
, “
Thermodynamics of acoustic black holes in two dimensions
,”
Adv. High Energy Phys.
2016
,
5710625
.
38.
A.
Kyriakou
,
E.
Neufeld
,
B.
Werner
,
G.
Székely
, and
N.
Kuster
, “
Full-wave acoustic and thermal modeling of transcranial ultrasound propagation and investigation of skull-induced aberration correction techniques: A feasibility study
,”
J. Ther. Ultrasound
3
,
11
(
2015
).
39.
R. K.
Dubey
, “
Data dependent stability of forward in time and centred in space (FTCS) scheme for scalar hyperbolic equations
,”
Int. J. Numer. Anal. Model.
13
,
689
704
(
2016
).
40.
R.
Hixon
, “
On increasing the accuracy of MacCormack schemes for aeroacoustic applications
,” in
3rd AIAA/CEAS Aeroacoustics Conference
(
AIAA
,
1997
), p.
1586
.
41.
S.
Dey
and
C.
Dey
, “
An explicit predictor-corrector solver with applications to Burgers' equation
,”
NASA Technical Memorandum, Technical Report No. 84402
(
National Aeronautics and Space Administration
,
1983
).
42.
W. C. R.
Nascimento
and
R. C.
Pestana
, “
An anti-dispersion wave equation based on the predictor-corrector method for seismic modeling and reverse time migration
,” in
SEG Technical Program Expanded Abstracts
(
Society of Exploration Geophysicists
,
2010
), pp.
3226
3230
.
43.
S.
Dey
, “
A novel explicit finite difference scheme for partial differential equations
,”
Math. Modell. Anal.
4
,
70
78
(
1999
).
44.
P. J.
Westervelt
, “
Parametric acoustic array
,”
J. Acoust. Soc. Am.
35
,
535
537
(
1963
).
45.
A.
Treves
,
L.
Maraschi
, and
M.
Abramowicz
, “
Basic elements of the theory of accretion
,”
Publ. Astron. Soc. Pac.
100
,
427
(
1988
).
46.
W.
Huang
and
R.
Russell
,
Adaptive Moving Mesh Methods
(
Springer
,
2011
), Vol.
174
.
47.
F.
Sewerin
and
S.
Rigopoulos
, “
An explicit adaptive grid approach for the numerical solution of the population balance equation
,”
Chem. Eng. Sci.
168
,
250
270
(
2017
).