A Numerical Method for Predicting Acoustical Wave Propagation in Open Spaces

Abstract

We present a numerical methodology for evaluating wave propagation phenomena in two dimensions in the time domain with focus on the linear acoustic second-order wave equation. An outline of the higher-order compact discretization schemes followed by the time discretization technique is first presented. The method is completed with the addition of spatial filtering based on the same compact schemes' principles. The important role of boundary conditions is subsequently addressed. Two popular ways to truncate the computational domain in the near field are presented and compared here: first the formulation of “absorbing conditions” in the form of partial differential equations especially for the origin and second the construction of an absorbing layer surrounding the domain, in which waves (after they have exited the domain) are attenuated and decayed exponentially. Subsequently, the method is assessed by recalling three benchmark problems. In the first where a Gaussian pulse is generated and propagated in a 2D rectangular domain, the accuracy and absorbability of the boundary conditions are compared. In the second, a similar situation is investigated but under curvilinear coordinates and under the presence of a solid body which scatters the pulse. Finally the sound field inducted by the flow of corotating vortex pair is calculated and compared with the corresponding analytical solution.

1. Introduction

In many computational wave propagation problems, including, for example, acoustics, electromagnetic, and fluid dynamic, low-order schemes (second or lower) are not accurate enough. The advantage of high-order finite-difference schemes is twofold: they allow one either to increase accuracy while keeping the number of mesh points fixed or to reduce the computational cost by decreasing the grid dimension while preserving accuracy. And although they require more work per node, the fact that fewer points need to be stored and computed makes them more efficient than low-order methods.

High-order finite difference schemes can be classified into two main categories, with respect to the way derivative calculation is made: explicit schemes and Padé type or compact schemes. Explicit schemes compute the derivative directly at each mesh point, while compact schemes compute the derivatives implicitly usually by solving a large system of equations. Compact schemes are more accurate than the corresponding explicit schemes of the same order as many investigators have shown [1]. Compact schemes can be optimized over a specific wavenumber range by minimizing their dispersion relation or resolving efficiency, thus providing more accuracy [1, 2]. The nondissipative nature of these schemes may give rise to spurious oscillations with higher wavenumber outside the working region of the scheme. These oscillations may grow during iterations, affect solution’s accuracy, and at worst drive the scheme to instability. Many authors propose the use of spatial filters to overcome this [3], while others avoid applying the filter directly to the solution and add damping terms to the discretized equation/s, which act implicitly and have zero contribution on small wave numbers [4, 5]. A similar sixth-order compact scheme with filtering is recalled in this work with a special boundary treatment for the imposition of von Neumann pure reflection conditions.

Apart from discretization, the most sensitive and important issue when dealing with numerical simulation of hyperbolic phenomena is the proper choice of boundary conditions. Radiation in the far field, reflection from obstacles, changes in the propagation medium must be accurately predicted, while keeping the solution bounded and the whole scheme stable. Radiation in the far field however seems to be the most difficult and challenging issue in numerical simulation of wave propagation. Since 30 years it has been and continues to be a topic of active research common in several application fields. Attempts to cope with the reflections generated from the truncation of the computational domain in the near field first took place by Engquist and Majda [6, 7] who developed theoretical perfectly absorbing boundary conditions (ABCs) for general classes of wave equations. Subsequently, by making use of Padé approximations of the solutions to the wave equation representing plane waves travelling to the left, they derived a hierarchy of stable highly absorbing local BCs, which approximate the theoretical nonlocal (exact) initially obtained. As the order of the BC increases the absorbing properties become successively better. Later Giles [8] derived local BCs for the Euler equations by using low-order Taylor series approximation, and Bayliss and Turkel [9] localized the exact BC by taking an asymptotic expansion of the solution for large distances from the source region. This formulation was also proposed by Tam and Webb [2] and has been used successively later by other authors like Bogey and Bailly [10] in three-dimensional cases. Another approach to the problem was used by Thompson [11, 12], who developed exact ABCs for the Euler equations based on characteristic analysis.

The absorbing properties of the above schemes depend on the incident angle and the boundary distance from the source. The further is the boundary location from the source, the better the absorbing properties become. Perfect absorption incurs only for plane waves. This means all that these formulations have a small degree of reflection which can affect the numerical solution. Another approach to the problem of finding good absorbing BCs are the so-called absorbing layers. These are used in conjunction with an ABC to enhance its efficiency or even without an ABC. First introduced by Israeli and Orszag [13] absorbing layer treatments typically provide damping of disturbances prior to their interaction with an ABC. This is done by introducing artificial dissipation to the equations by adding special terms or most simply by adding linear friction terms as done later by R. Kosloff and D. Kosloff [14] for the Schrödinger and the acoustic wave equations. Whatever disturbances are reflected by the boundary are similarly damped again as they propagate back through the layer before reaching the domain as reflection error. The basic problem with this approach is that the internal boundary of the layer is itself reflective, so the only way to obtain satisfactory results is to gradually increase damping from zero over a relatively long distance. This results in thick layers inefficient layers. The perfectly matched layer (PML) first developed by Berenger [15] for solving electromagnetism problems with the FDTD method (finite difference time domain) is an alternative way to construct layers with zero theoretical reflection at any frequency and any incident angle (perfect matching). The key idea of Berenger's formulation was an unphysical split of the radiated variable into its geometrical components (e.g., 𝑝^𝑥,𝑝^𝑦,𝑝^𝑧 for acoustic pressure) then enforcement of damping only in the direction perpendicular to the layer and then choice of the damping functions (𝜎_𝑥,𝜎_𝑦,𝜎_𝑧) at each direction by requiring a small reflection factor. Well posedness of this formulation was weak due to the unphysical field splitting, and implementation in the time domain was under debate due to the requirement of writing second-order wave equations as a first-order system (Euler equations for acoustics) and thus introducing more auxiliary variables (e.g., velocities for acoustics). A later more elegant, general, and strongly posed approach known as the stretched-coordinate PML was used first by Chew and Weedon [16] again in electromagnetism problems and later by Hu [17, 18] in acoustics. According to this formulation equations inside the PML are first transformed into the complex plane where damping is imposed only into the imaginary part of each coordinate. This causes waves to attenuate inside the PML in each direction no matter their incident angle. Physical coordinates equations are received by grouping together all imaginary components and then identifying them as auxiliary variables. To deal with second-order wave equations later in this work, the Laplace transform is used in the same manner. PMLs and absorbing layers have to be properly tuned in order to operate right. This is done only by experimentation and the tuneable parameters are the width of the layer and the damping profiles inside.

Another important aspect of wave propagation numerics is the choice of time discretization scheme. The well-known Runge-Kutta (RK) scheme and its variations have been well studied, improved, and widely used in acoustical wave propagation applications. Even recently Kondaxakis and Tsangaris [19] employ the RK-Nyström method to solve the resulting from a spectral discretization ordinary differential equation system to accurately model the acoustical wave propagation inducted from a pulsatile internal flow. Seo and Mittal [20] use the classical four-stage RK method to march in time the discretized perturbed Euler equations to successfully predict sound propagation in a number of complex geometry configurations with the immersed boundary method. Among other analysts, Bogey et al. [21] demonstrate the use of a class of optimized RK schemes to perform noise computations by solving directly the compressible Navier-Stokes equations without using decoupled analogy approaches. A similar time integration scheme firstly introduced by Hu et al. [22] is described later in this work.

The paper is organized as follows: Section 2 presents the considered initial boundary value wave problem along with its boundary conditions. Section 3 presents the numerical implementation regarding spatial and temporal discretization. Section 4 is devoted to numerical assessment of the method by recalling 3 well-known benchmark problems in wave propagation.

2. Governing Equations

We consider a time-dependent acoustical pressure field 𝑝^ propagating through unbounded two-dimensional space and assume that all sources and initial disturbances are confined to the curvilinear domain Ω. We further assume the speed of propagation 𝑐₀>0 to be constant. Inside Ω the field 𝑝(𝑥,𝑦,𝑡) satisfies𝑝𝑝^𝑡𝑡−𝑐²⁰∇²𝑝^=𝑆,𝑡>0,(1)_=𝑝𝑝₀(𝑥,𝑦),𝑡=0,(2)_𝑡=0,𝑡=0,(3) where 𝑆=𝑆(𝑥,𝑦) represents bounded sources inside the domain.

2.1. Lighthill's Acoustic Analogy

For the purposes of numerical assessment of the method using the spinning vortex pair benchmark problem later in this paper, Lighthill's acoustic analogy is recalled [23], which states that under isentropic propagation, pressure fluctuations can be determined form sources inducted by fluid flow by solving a second-order inhomogeneous wave equation of the following form:𝑝^𝑡𝑡−𝑐²⁰∇²𝑝^=𝑐²⁰𝜕²𝑇_𝑖𝑗𝜕𝑥_𝑖𝜕𝑥_𝑗,(4) where 𝑇_𝑖𝑗  represents the Lighthill's source tensor which has the form:𝑇_𝑖𝑗=𝜌𝑢_𝑖𝑢_𝑗−𝜏_𝑖𝑗.(5)

In (5), term 𝜌𝑢_𝑖𝑢_𝑗 represents nonlinear acoustic contributions due to convective forces, whereas term 𝜏_𝑖𝑗 represents acoustic contributions due to fluid's viscous stresses.

However, if the flow is incompressible it can be shown by taking the divergence of the momentum equation and then using the continuity equation to cancel terms, that the incompressible pressure satisfies the following equality:∇²𝑝^inc=𝜌𝜕²𝑢_𝑖𝑢_𝑗𝜕𝑥_𝑖𝜕𝑥_𝑗.(6) Under the assumption that in an incompressible flow pressure fluctuations can be split into the acoustic and incompressible part 𝑝′=𝑝^𝑎+𝑝^inc, Lighthill's analogy becomes upon substitution to (4) and using(6)𝑝^𝑎𝑡𝑡−𝑐²⁰∇²𝑝^𝑎=−𝑝^inc_𝑡𝑡.(7) Equation (7) provides an alternative way of computing acoustic pressure using only the second derivative of incompressible pressure as input and not the divergence of Lighthill's tensor which is computationally more demanding to calculate. Incompressible pressure data can be available from a fluid flow solver.

2.2. Boundary Conditions

In the following two paragraphs, a brief discussion about the two considered boundary condition formulations is attempted.

2.2.1. Local ABC's

Following the exact radiation boundary condition for the 1D homogenous wave equation and using pseudodifferential operators, we can easily construct ABCs for the 2D case [13]. The 2D homogenous wave equation is rewritten by taking 𝑐₀=1 as𝜕^2𝑥𝑥−𝐿²𝑝^=0,(8) where 𝐿² is the operator𝐿²=𝜕^2𝑡𝑡−𝜕^2𝑦𝑦.(9) The 2D analog of the exact 1D radiation condition would be written as𝜕_𝑥+𝐿𝑝^=0,(10) where 𝐿 is one of the “square roots” of the operator 𝐿². For plane waves traveling in the 𝑦 direction 𝑝^=𝑝^(𝑥)⋅𝑒^{𝑖𝜔𝑡+𝑖𝑘𝑦}, we can obtain an expression for the operator 𝐿² utilizing (9).𝐿𝐿²==±−√𝜔²+𝑘²,−𝜔²⋅1−𝑘²𝜔²𝐿=±𝑖𝜔,or1−𝑘²𝜔².(11) Now expanding the square root quantity in powers of 𝑘/𝜔;𝐿≈±𝑖𝜔−𝑘²2𝜔.(12) Changing back to time domain using 𝑖𝜔↔𝜕_𝑡 and 𝑖𝑘²↔𝜕_𝑦𝑦, substituting to (10), and taking one more time derivative, we obtain the boundary conditions for the 𝑥 direction:𝑝^𝑥𝑡±𝑝^𝑡𝑡∓12𝑝^𝑦𝑦=0.(13) Equivalently for 𝑦 direction:𝑝^𝑦𝑡±𝑝^𝑡𝑡∓12𝑝^𝑥𝑥=0.(14) The first of (13) and (14) stand for the right and top boundary and the second for the left and bottom boundary, respectively. (12) is the second approximation of the square root quantity in (11). Taking the first approximation (i.e., 𝐿≈±𝑖𝜔) would had led to the 1D exact radiation condition for each direction:𝑝𝑝^𝑥_𝑦=±𝑝_𝑡,=±𝑝_𝑡.(15) A comparison of different approximations of (11) is available in [24] where it is shown that the Chebychev-Padé-type approximation is leaving less residual energy in the 2D channel test. In [25] Trefethen and Halpern formulate theorems for determining the well posedness for these type of approximations. (15) are easier to implement but lacking in accuracy for 2D problems. Equations (13) and (14) are more difficult to implement, especially in variational formulation methods, because of the presence of mixed derivatives thus requiring usage of auxiliary variables. A more general form of (15) is𝑝^𝐧=𝑝^𝑡,(16) where 𝐧 denotes the normal direction.

Similar expressions were derived in [4, 6] but following the asymptotic solution of the Euler equations and/or using the wave equation in polar coordinates, respectively. The boundary expression for pressure can be written as𝜕_𝐧𝑝^+𝜕_𝑡𝑝^+12𝑅𝑝^=0.(17) Another expression for circle boundaries provided in [26] is𝜕_𝐧𝑝^+𝜕_𝑡𝑝^+12𝑅𝑝^=1+𝑅𝜕_𝑡⁻¹18𝑅𝑝^+12𝑅𝜕^2𝜃𝑝^.(18) The improvement over (17) however is not significant.

2.2.2. Perfectly Matched Layer (PML)

A brief presentation of the coordinate stretched PML equations derivation following [27, 28] is given below, considering wave propagation in a rectangular 2D domain Ω, surrounded by a layer of thickness 𝐿_𝑖,𝑖=1,2.

The Laplace transform of the pressure fluctuations is𝐿(𝑝′)=𝑝^̃(𝐱,𝑠)=^∞0𝑒^𝑠𝑡𝑝^(𝐱,𝑡)𝑑𝑡.(19) The Laplace transform of (1) gives the modified Helmholtz equation for 𝑝^̃, which after introducing the coordinate transformation (20) takes the form of (21).𝑥_𝑖⟶𝑥̃_𝑖∶=𝑥_𝑖+1𝑠^𝑥𝑖₀𝑠𝑍_𝑖(𝑥)𝑑𝑥,𝑖=1,2,(20)₂𝑝^̃=𝜕𝑥̃₁𝑐²⁰𝜕𝑥̃₁𝑝^̃+𝜕𝑥̃₂𝑐²⁰𝜕𝑥̃₂𝑝^̃.(21) The damping profiles 𝑍_𝑖 are positive inside the layer but zero in Ω. If 𝑝^̃ satisfies the modified Helmholz (21) under the new stretched coordinates, then the actual pressure fluctuations will remain unaltered inside Ω, but vanish exponentially inside the layer. Transforming (21) back to time domain and introducing auxiliary variables, one gets the partial differential equations which must be solved inside the layer. From (20) the partial differentiation relationship between coordinates is as follows:𝜕𝑥̃_𝑖=1+𝑍^𝑖𝑠⋅𝜕^𝑥𝑖=𝑘_𝑖⋅𝜕_𝑥𝑖.(22) After replacing derivatives in (21), we get𝑠²𝑘₁𝑘₂𝑝^̃=𝜕_𝑥1𝑐²⁰𝑘₂𝑘₁𝜕_𝑥1𝑝^̃+𝜕^𝑥2𝑐²⁰𝑘₂𝑘₁𝜕_𝑥2𝑝^̃.(23) Or equivalently𝑠²+𝑠𝑍₁+𝑍₂+𝑍₁𝑍₂𝑝^̃=+𝜕𝜕_𝑥1_𝑥1𝑐²⁰𝜕_𝑥1𝑝^̃+𝜕_𝑥2𝑐²⁰𝜕_𝑥2𝑝^̃𝑐²⁰𝑍₂−𝑍₁𝑠+𝑍₁𝜕_𝑥1𝑝^̃+𝜕^𝑥2𝑐²⁰𝑍₁−𝑍₂𝑠+𝑍₂𝜕_𝑥2𝑝^̃.(24) Introducing the auxiliary vector variable 𝐟̃=(𝑓₁,𝑓₂)=(𝑐²⁰(𝑍₂−𝑍₁)/(𝑠+𝑍₁)𝜕_𝑥1𝑝^̃,𝑐²⁰(𝑍₁−𝑍₂)/(𝑠+𝑍₂)𝜕_𝑥2𝑝^̃), and reverting to time domain:𝑝^𝑡𝑡+𝑍₁+𝑍𝐟₂𝑝^𝑡+𝑍₁𝑍₂𝑝^=∇⋅𝑐²⁰∇𝑝^+𝐟_𝑡=𝐀₁𝐟+𝑐²⁰𝐀₂∇𝑝^,(25) where matrices 𝐀₁ and 𝐀₂ are defined as follows:𝐀₁=diag−𝑍₁,−𝑍₂,𝐀₂=diag𝑍₂−𝑍₁,𝑍₁−𝑍₂.(26) The key difference from other stretched coordinate formulations was the use of the Laplace transform, which suits better here, because the considered pde is the second-order wave equation and we did not wish to split it to a first-order system (Euler equations) introducing more unknowns (velocities). Note that if 𝑍_𝑖=0 and 𝐟=𝟎 (25) takes the form of (1). Equation (25) must be solved inside the layer using the same time marching scheme as with (1). Auxiliary variable 𝐟 increases the calculation load as it is an additional unknown in the system. The layer can be tuned by properly selecting the widths 𝐿_𝑖 and the damping functions 𝑍_𝑖. This must be done only by experimentation. A general rule of thumb is that damping functions must be gradually increased from zero avoiding steep changes. It should be noted here that (25) was derived for a rectangular domain. However, it must be possible to use (25) in general curvilinear coordinates, if under an appropriate transformation one substitutes the gradient and divergence operators with the proper expressions. Nevertheless, this has not be done in this work. Benchmark problems where the PML BC was employed were solved in rectangular domains only.

2.2.3. Solid Body Condition

When a acoustic wave incidents with a solid surface, then the following homogenous von Neumann boundary condition is valid: 𝜕_𝐧𝑝^=0.(27) Equation (27) is adequate to provide total reflection for incident waves, and no other expressions need to be formulated like in Euler equations where the normal acoustic velocity must also be set to zero (𝐮⋅𝐧=0). This is an additional benefit of modeling directly the second-order wave equation and not the solving for all the primitive acoustic variables.

3. Numerical Implementation

3.1. Coordinate Transformation

In order to deal with complex geometries, the governing equation first have to be transformed into curvilinear coordinates using a Jacobian transformation. Metrics must be calculated by the same spatial discretization scheme to retain the same order of accuracy.𝐽=𝜕(𝜉𝜉,=𝑛)𝜉(𝑥,𝑦),𝜂=𝜂(𝑥,𝑦),𝜕(𝑥,𝑦)=det𝑥^𝜉𝑥_𝑛𝑦_𝜉𝑦_𝑛=𝑥^𝜉𝑦_𝑛−𝑥_𝑛𝑦_𝜉.(28) Using the differentiation chain rule (1) transforms to𝑝^𝑡𝑡=𝑐²⁰+𝑝𝜉^2𝑥_𝑛𝑛+𝜉^2𝑦𝑝^𝜉𝜉+2𝜉_𝑥𝑛_𝑥+𝜉_𝑦𝑛_𝑦𝑝^𝜉𝜂+𝑛^2𝑥+𝑛^2𝑦+𝜉_𝑥𝑥+𝜉_𝑦𝑦𝑝^𝜉+𝑛_𝑥𝑥+𝑛_𝑦𝑦𝑝^𝑛+𝑆.(29) The normal derivative which appears in boundary conditions expressions transforms to𝜕_𝐧𝑝^=𝑝^𝜉𝜉_𝑥𝑛_𝑥+𝜉_𝑦𝑛_𝑦±𝑝^𝑛𝑛^2𝑥+𝑛^2𝑦𝑛^2𝑥+𝑛^2𝑦^−1/2.(30) The above stands for the top/bottom boundary, and if gridlines are vertical at that region (most grid generating techniques meet well this requirement) then 𝑥_𝜉𝑥_𝑛+𝑦_𝜉𝑦_𝑛=−𝑛_𝑦𝜉_𝑦−𝑛_𝑥𝜉_𝑥=0 and therefore the normal derivative can be calculated in the transformed coordinate system simply by𝜕_𝐧𝑝^=±𝑝^𝑛𝑛^2𝑥+𝑛^2𝑦.(31) The expression for left/right boundaries is similar:𝜕_𝐧𝑝^=±𝑝^𝜉𝜉^2𝑥+𝜉^2𝑦.(32) So the wall conditions take the form:𝑝^𝜉=0,𝑝^𝑛=0.(33)

3.2. Spatial Discretization—Compact Schemes

Compact schemes (CS) have become popular in finite difference numerical applications due to their simplicity and easiness in implementation. In this work, a sixth-order CS was employed and a special formulation for the solid boundary was derived. The general form of tridiagonal compact approximations for the first derivative [1] is (subscript here denotes discrete value at that point and not partial differentiation)𝑎𝑝^𝑖−1^+𝑝^𝑖^+𝑎𝑝^𝑖+1^=𝑏𝑝^𝑖+1−𝑝^𝑖−12ℎ+𝑐𝑝^𝑖+2−𝑝^𝑖−24ℎ.(34) And for the second derivative,𝑎𝑝^𝑖−1^+𝑝^𝑖^+𝑎𝑝^𝑖+1^=𝑏𝑝^𝑖+1−2𝑝^𝑖+𝑝^𝑖−1ℎ²+𝑐𝑝^𝑖+2−2𝑝^𝑖+𝑝^𝑖−24ℎ².(35) Unknown coefficients are calculated by taking Taylor expansions of the left- and right-hand side terms around 𝑝^𝑖, then gathering together like terms and finally equating quantities of the same power of ℎ. By requiring sixth-order accuracy, the resulting linear system has 3 unknowns and 2 equations. Thus a family of one parameter schemes can be generated. As this parameter goes to zero, the family merges to the fourth-order central difference scheme.

The arithmetic values of the coefficients for the first derivative are [1]𝑎=13,𝑏=194,𝑐=91.(36) And for the second derivative,𝑎=211,𝑏=1112,𝑐=131.(37) Implementation to the curvilinear system is straightforward by taking ℎ=1.

3.3. Boundary Treatment

To calculate derivatives when nonperiodic conditions are concerned, then one-sided formulas need to be constructed, to avoid using points outside the domain. They must be constructed such that the banded (tridiagonal) form of the interior scheme is maintained. Also their stencil size can be adjusted to match the formal order of the accuracy of the corresponding interior scheme. Using the same method as above, the following schemes for the 1st derivative were constructed

point 1: 𝑝^1^+5𝑝^2^=19760𝑝^1−512𝑝^2+5𝑝^3−53𝑝^4+512𝑝^5−120𝑝^6ℎ1,(38) point 2:18𝑝^1^+𝑝^2^+34𝑝^3^=−4396𝑝^1−56𝑝^2+98𝑝^3−16𝑝^4+196𝑝^5ℎ1.(39) And for the 2nd derivative,

point 1:𝑝^1^+11216𝑝^2^=2077157𝑝^1−2943110𝑝^2+57344𝑝^3−16799𝑝^4+1811𝑝^5−57110𝑝^6+1311980𝑝^7ℎ1₂,(40) point 2:211𝑝^1^+𝑝^2^−13122𝑝^3^=17788𝑝^1−50744𝑝^2+78344𝑝^3−20122𝑝^4+8188𝑝^5−344𝑝^6ℎ1₂.(41) The 1st derivative coefficients for the right boundary formulas are the negative of the coefficients of (38) and (39), whereas the 2nd derivative coefficients are the same as in (40) and (41).

3.4. Boundary Conditions

Dirichlet boundary conditions are easy and obvious to implement with CS's. Neumann boundary conditions require more care. The homogenous wall boundary condition (33) was applied using the following one-sided explicit formula:𝑎𝑝^1=_𝑖=2𝑏_𝑖𝑝^𝑖.(42) Then the one-sided first derivative formulas (38)-(39) can be used to form the tridiagonal system of equations, where the value of 𝑝^1 must be substituted from (42) in order to calculate the right-hand side expressions. The Taylor matching procedure is used again to determine the unknown coefficients in (42). Each term is expanded about node 1 and like terms are equated to arrive at certain orders of accuracy. However, in this particular case the sum of coefficients of ℎ can be set to any nonzero constant (here we used 1 to keep calculations simple), due to the presence of the homogenous condition for the first derivative. For fourth-order accuracy, the coefficients of (42) are𝑎=48,𝑏₂=−36,𝑏₃=16,𝑏₄=−3.(43) In the case of the second derivative, however, a different approach has been followed in order to incorporate the Neumann boundary condition. The following relation was formed:𝑎⋅𝑝^1^+𝑏ℎ𝑝^1^+𝑐⋅𝑝^2^=1ℎ²^𝑖=1𝑑_𝑖⋅𝑝^𝑖.(44) In (44), one may notice the presence of the derivative at the first point (𝑝^1)^  in the left-hand side. (44) is more general and can be used even for nonhomogenous Neumann conditions. Requiring fourth-order accuracy again (thus 2 terms on the right-hand side of (44)), the arithmetic values of coefficients are𝑎=13,𝑏=1,𝑐=76,𝑑¹=−1,𝑑₂=1.(45) For the second node and the node before the last, the fourth-order central scheme was used with coefficients for the first and second derivative, respectively,𝑎=14,𝑏=32,(46)𝑎=110,𝑏=65.(47) No special treatment was needed for the other types of boundary conditions (13)–(17), where a transient term is included. The derivatives were calculated using the one-sided formulas (38)–(41).

3.5. Error Analysis

Taking the Fourier approximation of the depended variable, substituting into (34), and introducing the modified wavenumber 𝜔=ℎ𝜔, we can build the spectral transfer function for the considered CS. It is well known [1] that CS's are nondissipative; therefore, the spectral function for the interior points has no imaginary part. In Figure 1(a), a plot of the 1st derivative transfer function for 5 different in accuracy schemes is depicted. Diagonal corresponds to exact differentiation. The superiority over the classical second-order differentiation is obvious.

(a)

(b)

Figure 1 

1st derivative spectral functions (a) interior and (b) boundary points (imaginary part).

However, as can be seen in Figure 1(b) some dissipation is added in the boundary one-sided formulas at the high bandwidth. This dissipative error is increased as the order of the scheme and the deviation from noncentered differentiation increases. This is the reason why we dropped the order of our scheme as we approached the boundaries to 4, although high wavenumbers never came into consideration throughout this study.

3.6. Linear Equations System Solver

The classical tridiagonal matrix algorithm (TDMA) was employed to solve the resulting from discretization sparse linear equations system. The main benefit was the low storage requirements of this method, as it only requires storage of the three main diagonals of the sparse matrix.

3.7. Spatial Filtering

Stability analysis of a CS is not always a straightforward procedure. Some analysis can be done by considering the spectral function, but for numerical purposes stability is achieved by adding damping to the scheme. This can be done either by adding terms into the partial differential equation, or by directly filtering the solution after each time advancement. For the CS presented here, this was done by using the following filter formula [3]: 𝑎_𝑓𝑝̂^𝑖−1+𝑝^𝑖̂+𝑎_𝑓𝑝̂^𝑖+1=^𝑁_𝑛=0𝑎_𝑛2𝑝^𝑖+𝑛+𝑝^{𝑖−𝑛}.(48) If a component of the solution vector is denoted by 𝑝^𝑖, filtered values 𝑝^𝑖̂ are obtained by solving again a tridiagonal system. (48) corresponds to a 2𝑁-order filter and a 2𝑁+1 point stencil. Coefficients are calculated again using Taylor series analysis and can be found along with error analysis in [29]. The adjustable parameter 𝑎_𝑓 satisfies the inequality  −0.5≤𝑎_𝑓≤0.5, with higher values corresponding to a less dissipative filter, whereas analysts suggest [3] that values between 0.4 and 0.5 are appropriate for a wide range of applications. As far as the order of the filter is concerned, a general rule is to choose the filter order to be at least 2 orders higher than the CS it is coupled with, meaning that an eighth-order filter should coupled with our sixth-order scheme. Special formulas are required for the boundary points due to the relatively large stencil of the filter. In this work, a more suitable approach was followed: the filter order was gradually decreased upon approaching the boundary and the outer points were left unfiltered. In multidimensional problems, the filter must be applied one time for each spatial direction. However, there is no rule regarding the filter application frequency. This can found be experimentation and taking into consideration that applying the filter after every time advancement increases heavily the computational load.

3.8. Time Integration

The 4-stage Low-Dispersion Low-Dissipation Runge Kutta (LDDRK4) scheme has been employed to march the wave equation, the ABCs, and the PML equations in time [22]. The main property of this modified RK scheme is that it preserves the dispersion relation between the ordinary differential equation and its numerical approximation by minimizing a quadratic norm of the amplification factor difference over a wide range of wavenumbers. For reasons of completeness, we present the dispersion and dissipation errors in conjunction with the classic RK4 method in Figures 2(a) and 2(b) along with the stability (𝑅) and accuracy (𝐿) limits for a given precision of  𝜀=0.001. The absolute value of the complex amplification factor 𝑛 represents the dissipation errors, whereas dispersion errors are represented by the phase.

(a)

(b)

Figure 2 

Magnitude (a) and phase (b) of the relative amplification factor for the RK schemes.

For the a linear first-order system of equations 𝐔′=𝐅(𝐔), the algorithmic part of the method can be written as

calculate:𝐊_𝑖=𝑑𝑡⋅𝐅𝐔^𝑛+𝑤_𝑖𝐊_𝑖−1,𝑖=1,…,4,(49) update:𝐔^𝑛+1=𝐔^𝑛+𝐊₄,(50) where the optimized weight values are 𝑤₁=0,𝑤₂=1/4,𝑤₃=1/3,𝑤₄=1/2 [30].

An additional benefit of this formulation over the traditional RK4 method is the lower storage requirement. The loop implementing (49) needs to store only two values of 𝐤_𝑖 (the current and previous ones) and not all four as in the traditional RK4. This becomes more beneficial if the number of stages of the method is increased.

To use this method with (1), we have to create first a first-order system by incorporating an auxiliary variable 𝑞. Note that this not the same as solving the Euler equations where an additional auxiliary variable would be needed. Thus the differential system becomes𝑞𝑝^𝑡=𝑞,_𝑡=𝑐²⁰∇²𝑝^.(51) The same practice was followed for the ABCs. So equations (13)-(14) become, after reintroducing 𝑐₀,𝑞_𝑡=±𝑐₀12𝑝^𝑦𝑦𝑝_𝑡,=∓𝑞±𝑝^𝑥𝑐₀.(52) And the PML equation (25) becomes𝑞𝑝^𝑡=𝑞,_𝑡+𝑍₁+𝑍₂𝐟𝑝^𝑡+𝑍₁𝑍₂𝑝^=∇⋅𝑐²⁰∇𝑝^+𝐟_𝑡,=𝐀₁𝐟+𝑐²⁰𝐀₂∇𝑝^.(53) We note here that no special memory management was followed for the PML equations. This means that the vector variable 𝐟 was defined over the whole computational domain, despite that it was solved only for the layer nodes. Nodes outside the layer were assigned with zero values and were never updated. This implementation is simply and straightforward but sacrifices computational resources.

4. Test Cases

4.1. Simple Gaussian Pulse Propagation

With this numerical experiment, we compare the accuracy, absorbing ability, and long time stability of the ABCs (52) and the PML (53). Equations (51) were solved in a rectangular domain with a Gaussian pulse released at its center, and a reference solution (RS) used for comparisons was computed in a 5 times larger domain keeping the same discretization density. The sixth-order CS ((34)-(35) and (38)-(39)) and the LDDRK4 method (49)-(50) were used for spatial and temporal discretization, respectively. Stable time step has been calculated using the LDDRK4 stability limits from Figure 2. Because of the grid uniformity, no filtering scheme was necessary. More details can be found in Table 1 where SI units are assumed, but have been omitted for simplicity.

Table 1 

Numerical details for Case  1.

PML details are summarized in Table 2, where 𝑏_𝑖 and 𝐿_𝑖 are the PML's inner interface coordinate components and width, respectively. Among other choices for damping profiles, we found that the profile in Table 2 provides very smooth transition from the inner PML interface. The width and damping constant choices are adequate to provide enough attenuation of the wave front before it even reaches the outer boundary, thus making the choice of the outer boundary condition trivial. For testing purposes however, a reflective BC (Neumann) was imposed on the outer boundary.

Table 2 

PML details for Case  1.

In Figure 3, pressure distributions were drawn at times where the wave front reaches and exits the inner PML interface and the outer computational boundary, where the Neumann BC is valid. No spurious reflections were observed. However, this was not the case with the ABC. In Figure 4, a comparison between the RS and the ABC solution is plotted. Pressure profiles are shown at two discrete time steps. At time 𝑡=30 when the pulse has already left the domain, one can clearly observe reflection waves traveling in the opposite direction. However, the magnitude of these waves is just 0.89% of the incident wave and can range up to 1.0% depending on the incident angle. This shows a good but not reflectionless operation of the ABCs.

Figure 3 

Pressure history—PML BC.

Figure 4 

Pressure profiles along 𝑥 axis at 𝑦=0.

A more extensive analysis of the absorbability of the two boundary condition types was further attempted. First the normalized wave energy defined as in (54) was calculated over the whole domain for each type of boundary condition and then the L2 norm of the error from the RS was calculated for a long period of time. Both types showed good stability characteristics over time, but the PML clearly outperforms the ABCs as far as the accuracy and the residual energy are concerned. In Figures 5 and 6, one can observe that the normalized energy for the PML case is at least 4 orders of magnitude smaller, whereas the deviation from the reference solution is at least 2 orders of magnitude smaller.𝐸=12_Ω𝑐²⁰𝜕_𝑥𝑝^²+𝑐²⁰𝜕_𝑦𝑝^²+𝜕_𝑡𝑝′²𝑑Ω.(54)

Figure 5 

Absorbing ability comparison of ABC and PML.

Figure 6 

Solution accuracy comparison of ABC and PML.

4.2. Spinning Vortex Pair

An important motivation for investigating this type of vortex sound problem is the existence of an analytical solution for the acoustic far field, which allows direct validation of the acoustic results. The corotating vortex pair consists of two point vortices separated by a fixed distance of 2𝑟₀ with circulation intensity Γ. A schematic is presented in Figure 7. The vortices rotate around each other with a period 𝑇=8𝜋²𝑟²⁰/Γ. Each vortex induces an angular velocity 𝑣_𝜃=Γ/(4𝜋𝑟₀). The configuration results in a rotating speed 𝜔=Γ/(4𝜋𝑟²⁰) and rotating Mach number 𝑀_𝑟=𝑣_𝜃/𝑐₀= Γ/(4𝜋𝑟₀𝑐₀)=2𝜋𝑟₀/𝑇𝑐₀.

Figure 7 

Spinning vortex pair schematic.

The incompressible inviscid flow field can de determined numerically [31] by the evaluation of the complex potential function Φ(𝑧,𝑡):Φ(𝑧,𝑡)=Γ2𝜋𝑖ln𝑧²1−𝑧𝑏²²≈𝜋Γ𝑖ln𝑧−2Γ𝜋𝑖𝑧𝑏²²,(55) where 𝑧=𝑟𝑒^𝑖𝜃 and 𝑏=𝑟₀𝑒^𝑖𝜔𝑡. The approximation in (55) is valid only in the far field that is for |𝑧/𝑏|≫1. The first term on the right-hand side of (55) represents a steady vortical flow and the second term represents the contribution due to the vortex motion [32]. From (55), it is easy to derive the hydrodynamic quantities such as velocities and pressure field, by differentiating with respect to 𝑧 and using the unsteady Bernoulli equation, respectively.𝑢_𝑥−𝑖𝑢_𝑦=𝜕Φ𝜕𝑧=𝑖Γ𝜋𝑧²−𝑧𝑏²𝑝,(56)_inc=𝑝₀−𝜌₀𝜕(Re(Φ))𝜕𝑡−12𝜌⁰𝑢^2𝑥+𝑢^2𝑦.(57) Using (57) and differentiating twice, it is easy to calculate the acoustic sources in (7). The result is𝑝^inc_𝑡𝑡=−𝜌₀Φ⃛_Re+𝑢̇^2𝑥+𝑢_𝑥𝑢̈_𝑥+𝑢̇^2𝑦+𝑢_𝑦𝑢̈_𝑦,(58) where 𝜌₀ is the fluid’s density and Φ⃛_Re represents the third time derivative of the real part of the complex potential function and 𝑢_𝑥,𝑢_𝑦 are the velocity components from (56). Care must be taken when evaluating (58) since the hydrodynamic quantities have very steep and large gradients near the point vortices. Recall that the complex potential from (55) is valid only in the far field. To avoid this singularity, a cut-off practice [33] was applied. Thus for nodes located at distances 𝑟/𝑟₀≤1.5, no acoustic source was calculated. As in the previous test case, the sixth-order CS ((34)-(35) and (38)-(39)) and the LDDRK4 method (49)-(50) were used for spatial and temporal discretization, respectively. No spatial filtering was necessary. Numerical details for the simulation can be found in Table 3 (SI units are omitted).

Table 3 

Numerical details for Case  2.

The PML technique was used again to truncate the computational domain with the same profile and damping constant as in case  1, except for the width, which is 𝐿_𝑖=80,𝑖=1,2 and the additional number of nodes per side which is 50. For comparison with the numerical results, the analytical solution of the acoustic pressure fluctuations from the vortex pair obtained with the matched asymptotic expansion method was used [34]𝑝^(𝑟,𝜃,𝑡)=𝜌⁰Γ⁴64𝜋³𝑟²⁰𝑐²⁰−𝑌𝐽₂₂(𝑘𝑟)⋅cos(2𝜃−2𝜔𝑡)(𝑘𝑟)⋅sin(2𝜃−2𝜔𝑡),(59) where 𝑘=2𝜔/𝑐₀ and 𝐽₂,𝑌₂ are the second-order Bessel functions of the first and second kind, respectively. Figures 8(a) and 8(b) show the quadrupole source term obtained form the computation at 𝑡=120 as a contour and a surface plot, respectively. The singularity issue can clearly be seen. As can be observed, the source region is large enough to avoid significant truncation of the sources. The amplitudes at the boundary of the source region are about 1% of their maximum amplitude reached within the source domain. In Figure 9, predicted pressure distributions at 𝑡=120 and 𝑡=280 are depicted (PML region not shown). Again the PML worked pretty well absorbing the outgoing waves, without affecting the solution's accuracy in the inner domain, as no spurious reflections were observed in the distributions. Further the predicted solution accuracy was tested against the analytical solution. To this end, the pressure values along the half 𝑥 axis at 𝑡=280 were plotted and compared with those from (59) in Figure 10, whereas in Figure 11 a time series comparison between the predicted pressure values at (𝑥,𝑦)=(0,80) and the ones from (59) is attempted. As can be seen, the results are in very good agreement with the analytical solution. Note the analytical solution is valid for the far field nodes, which here is minimally defined as two acoustic wavelengths away from the vortices, which is 𝑟/𝑟₀≥80. The discrete evaluation of the acoustic sources can be one reason for the observed lower amplitudes in the results, corresponding in general to around 95% of the analytical values. Another reason would be the applying of cutoff values near the vortices region when computing the sources.

(a)

(b)

Figure 8 

Acoustic sources (a) contour and (b) surface plot at 𝑡=120.

(a)

(b)

Figure 9 

Predicted pressure fields at (a) 𝑡=120 and (b) 𝑡=280.

Figure 10 

Pressure decay along the 𝑥 axis at 𝑡=280.

Figure 11 

Pressure time series at (𝑥,𝑦)=(0,80).

4.3. Acoustic Scattering from a Cylindrical Solid Body

The objective of this problem first presented in the 2nd CAA Workshop on Benchmark problems [35] is to find the sound field generated by propeller scattered off by the fuselage of on aircraft (see Figure 12). The pressure loading on the fuselage is an input to the interior noise problem. The fuselage is represented as a circular cylinder of diameter 𝐷 at (𝑥,𝑦)=(0,0), and the noise source is located at (𝑥,𝑦)=(4𝐷,0). The sound source can be either periodic, or nonperiodic where in the later case a pulse is released at 𝑡=0 and the initial pressure distribution is given as 𝑝(𝑥,𝑦,0)=𝑒^{(−ln2/0.04)((𝑥−4𝐷)2+𝑦2)}. One must find the pressure time series at points 𝐴(𝑟=5,𝜃=90⁰), 𝐵(𝑟=5,𝜃=135⁰), and 𝐶(𝑟=5,𝜃=180⁰) and compare them with distributions from an analytical expression.

Figure 12 

Acoustic scattering problem schematic.

Since problem is axisymmetric, computations are conducted only for the upper half of the field. Computational domain is chosen as a semicircular region, which is bounded by an outer semicircle, the cylinder wall and the 𝑥 axis. An O-type mesh with 𝜉 lines as concentric circles, with the first and last circles being, respectively, the cylinder and the outer boundary, and the 𝑛 lines directed radically from the cylinder to the outer boundary with diameter equal to 10.5𝐷. The absorbing boundary condition (18) was posed on the outer boundary, whereas the solid wall condition (27) was posed on the cylinder surface and on the symmetry axis. The initial grid density corresponding to (201 𝜉 lines) × (91 𝑛 lines) was not adequate to resolve the problem well. Additionally spurious oscillations were observed as the pulse propagated towards the cylinder. The polar grid nonuniformity was the major reason for the appearance of spurious oscillations. Doubling the grid density to (401 𝜉 lines) × (181 𝑛 lines) and introducing a eighth-order spatial filter every 10 iterations proved a necessity towards better quality solutions and elimination of the spurious oscillations. As in the previous test case, the sixth-order CS ((34)-(35) and (38)-(39)) and the LDDRK4 method ((49)-(50)) were used for spatial and temporal discretization, respectively. More numerical details can be found in Table 4.

Table 4 

Numerical details for Case  3.

Figure 13 shows a picture of the acoustic wave pattern at 𝑡=8.20. There are three wave fronts. The one farthest from the cylinder is the wave front created by the initial condition. The next front is the wave reflected off the right surface of the cylinder directly facing the initial pulse. The wave front closest to the cylinder is generated when the two parts of the initial wave front (split by the cylinder) collided and merged together to the left side of the cylinder. the later wave front is weaker relative to the other two. ABCs worked well despite their small reflections (Case  1), as no serious disruption can be seen in the figure. In addition, the reflection condition formula (44)–(47) was also worked well, keeping the algorithm stable.

Figure 13 

Acoustic wave pattern at 𝑡=8.20.

Using the method of superposition for the incident and the reflected wave, the analytical solution for the velocity potential, after correcting some errors found in [35], is given by the following equation.𝐴=(−𝑟,𝑒𝜃,𝑡)_{−𝑏𝜔/4}𝜙(𝑟,𝜃,𝑡)=^∞0𝐴(𝑟,𝜃,𝜔)sin(𝜔𝑡)𝑑𝜔,2𝑏𝐽⁰𝜔+𝑟²+𝑥^2𝑠−2𝑟𝑥_𝑠cos𝜃_∞_𝑘=0Re𝜀^𝑘𝐻_𝑘(𝑟𝜔)cos(𝑘𝜃)⋅2𝑘𝜋⋅𝐻_𝑘(𝜔/2)−𝜋𝜔⋅𝐻_𝑘+1(𝜔/2)_𝜋0𝜔0.5−𝑥^𝑠cos(𝜂)×𝐽₁0.25+𝑥^2𝑠−𝑥_𝑠cos(𝑛)𝜔0.25+𝑥^2𝑠−𝑥_𝑠cos(𝑛)cos(𝑘𝑛)𝑑𝑛,(60) where 𝐽₀,𝐽₁ are the zeroth and first-order Bessel functions of the first kind and 𝐻_𝑘 is the 𝑘th order Hankel function of the first kind. 𝜀_𝑘=2 for 𝑘>0 and 𝜀_𝑘=1 for 𝑘=0, 𝑏=ln2/0.04 and 𝑥_𝑠 is the source location in the 𝑥 axis. The pressure can then be found as the first time derivative:𝑝(𝑟,𝜃,𝑡)=−𝜕𝜙𝜕𝑡=−^∞0𝐴(𝑟,𝜃,𝑡)𝜔cos(𝜔𝑡)𝑑𝜔.(61) Care must be taken when numerically evaluating (60)-(61), since Bessel and Hankel functions go to infinity at zero. A value of 𝑘=55 is adequate to terminate the summation in (60).

Figure 14 shows the time history of pressure variations at the prescribed measurement points along with exact solutions. There is excellent agreement between the computed and the exact time histories.

Figure 14 

Pressure time histories at monitor points.

5. Conclusions

A numerical method for solving the second-order wave equation by combining a sixth-order compact spatial discretization scheme with the Low dispersion—dissipation Runge—Kutta time advancement scheme has been proposed. A slightly different formulation of coordinate stretching based on the Laplace transform for the Perfectly Matched Layer has been derived. This formulation was applied to rectangular domains only. The incorporation of Neumann boundary conditions with the sixth-order compact scheme has been discussed and a formula for the second derivative has been derived. Local ABCs formulations have been recalled and presented for comparison reasons. Three benchmark problems have been presented where the accuracy of the method and the boundary conditions have been tested. It was showed that local ABCs have small reflections and leave four orders of magnitude more residual energy inside the computational domain, whereas a PML with proper tuning can outperform them. However, the PML's auxiliary variable increases the computational load and memory requirements. Filtering was proved to be a necessity especially in problems where grids with a high level of nonuniformity are used. Finally the method gave excellent result agreement in two benchmark problems with analytical solutions with and without the presence of a source term in the right-hand side.

Nomenclature

Symbols

𝑝:	Pressure fluctuations
𝑝𝑎:	Acoustic pressure
𝑝inc:	Pressure due to incompressible fluid flow
𝑐0:	Speed of sound at normal conditions
𝜌:	Medium density
𝑢𝑖:	Velocity components
𝑇𝑖𝑗:	Lighthill tensor components
𝜏𝑖𝑗:	Viscous stresses tensor components
𝑘:	Wavenumber
𝜔:	Angular frequency
𝑠:	Laplace space variable
𝑍𝑖, 𝑘𝑖:	Damping profile in the 𝑖  direction, damping constant in the 𝑖  direction
𝑎,𝑏,𝑐,𝑑,𝑒,…:	Compact scheme coefficients
𝑎𝑓:	Filtering constant
𝑆:	Any kind of source distribution
Φ:	Complex potential function.

Subscripts

𝑡𝑡,𝑥𝑥,𝑡,𝑥:	First- or second-order partial differentiation with respect to time or space variables
𝐧:	Normal derivative
𝜉,𝑛:	Partial differentiation with respect to curvilinear space coordinates
𝑖:	Variable value at the discrete point 𝑖
𝑜:	Initial variable distribution.

Acknowledgments

This paper is part of the 03ED-866 research project, implemented within the framework of the “Reinforcement of Human Research Manpower” (PENED) and cofinanced by National and Community Funds (20% from Greek Ministry of Development—General Secretariat of Research and Technology and 80% from E.U.—European Social Fund).