Table of Contents
ISRN Mechanical Engineering
Volume 2011, Article ID 174031, 15 pages
http://dx.doi.org/10.5402/2011/174031
Research Article

A Numerical Method for Predicting Acoustical Wave Propagation in Open Spaces

Fluids Section, School of Mechanical Engineering, National Technical University of Athens, 9 Heroon Polytechniou Avenue, Zografou, 15780 Athens, Greece

Received 31 January 2011; Accepted 17 March 2011

Academic Editors: S.-H. Chuang and D. Zhou

Copyright © 2011 Johnny Papageorgakopoulos and Sokrates Tsangaris. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

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>0 to be constant. Inside Ω the field 𝑝(𝑥,𝑦,𝑡) satisfies𝑝𝑡𝑡𝑐202𝑝𝑝=𝑆,𝑡>0,(1)=𝑝0𝑝(𝑥,𝑦),𝑡=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:𝑝𝑡𝑡𝑐202𝑝=𝑐20𝜕2𝑇𝑖𝑗𝜕𝑥𝑖𝜕𝑥𝑗,(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:2𝑝inc𝜕=𝜌2𝑢𝑖𝑢𝑗𝜕𝑥𝑖𝜕𝑥𝑗.(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)𝑝𝑎𝑡𝑡𝑐202𝑝𝑎=𝑝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 𝑐0=1 as𝜕2𝑥𝑥𝐿2𝑝=0,(8) where 𝐿2 is the operator𝐿2=𝜕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 𝐿2. For plane waves traveling in the 𝑦 direction 𝑝=𝑝(𝑥)𝑒𝑖𝜔𝑡+𝑖𝑘𝑦, we can obtain an expression for the operator 𝐿2 utilizing (9).𝐿2=𝜔2+𝑘2,𝐿=±𝜔2𝑘12𝜔2,or𝐿=±𝑖𝜔𝑘12𝜔2.(11) Now expanding the square root quantity in powers of 𝑘/𝜔;𝑘𝐿±𝑖𝜔22𝜔.(12) Changing back to time domain using 𝑖𝜔𝜕𝑡 and 𝑖𝑘2𝜕𝑦𝑦, 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𝜕𝐧𝑝+𝜕𝑡𝑝+1𝑝2𝑅=0.(17) Another expression for circle boundaries provided in [26] is𝜕𝐧𝑝+𝜕𝑡𝑝+1𝑝2𝑅=1+𝑅𝜕𝑡11𝑝8𝑅+1𝜕2𝑅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𝑠𝑥𝑖0𝑍𝑖𝑠(𝑥)𝑑𝑥,𝑖=1,2,(20)2̃𝑝=𝜕̃𝑥1𝑐20𝜕̃𝑥1̃𝑝+𝜕̃𝑥2𝑐20𝜕̃𝑥2̃𝑝.(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𝑠2𝑘1𝑘2̃𝑝=𝜕𝑥1𝑐20𝑘2𝑘1𝜕𝑥1̃𝑝+𝜕𝑥2𝑐20𝑘2𝑘1𝜕𝑥2̃𝑝.(23) Or equivalently𝑠2𝑍+𝑠1+𝑍2+𝑍1𝑍2̃𝑝=𝜕𝑥1𝑐20𝜕𝑥1̃𝑝+𝜕𝑥2𝑐20𝜕𝑥2̃𝑝+𝜕𝑥1𝑐20𝑍2𝑍1𝑠+𝑍1𝜕𝑥1̃𝑝+𝜕𝑥2𝑐20𝑍1𝑍2𝑠+𝑍2𝜕𝑥2̃𝑝.(24) Introducing the auxiliary vector variable ̃𝑓𝐟=(1,𝑓2)=(𝑐20(𝑍2𝑍1)/(𝑠+𝑍1)𝜕𝑥1̃𝑝,𝑐20(𝑍1𝑍2)/(𝑠+𝑍2)𝜕𝑥2̃𝑝), and reverting to time domain:𝑝𝑡𝑡+𝑍1+𝑍2𝑝𝑡+𝑍1𝑍2𝑝𝑐=20𝑝𝐟+𝐟𝑡=𝐀1𝐟+𝑐20𝐀2𝑝,(25) where matrices 𝐀1 and 𝐀2 are defined as follows:𝐀1=diag𝑍1,𝑍2,𝐀2𝑍=diag2𝑍1,𝑍1𝑍2.(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𝑝𝑡𝑡=𝑐20𝜉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𝑝𝑖1𝑝2+𝑐𝑖+2𝑝𝑖24.(34) And for the second derivative,𝑎𝑝𝑖1+𝑝𝑖𝑝+𝑎𝑖+1𝑝=𝑏𝑖+12𝑝𝑖+𝑝𝑖12𝑝+𝑐𝑖+22𝑝𝑖+𝑝𝑖242.(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]1𝑎=3,𝑏=1491,𝑐=9.(36) And for the second derivative,2𝑎=11,𝑏=12311,𝑐=11.(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𝑝+52=197𝑝6015𝑝122+5𝑝353𝑝4+5𝑝1251𝑝2061,(38) point 2:18𝑝1+𝑝2+34𝑝3=43𝑝96156𝑝2+98𝑝316𝑝4+1𝑝9651.(39) And for the 2nd derivative,

point 1:𝑝1+126𝑝112=2077𝑝15712943𝑝1102+573𝑝443167𝑝994+18𝑝11557𝑝1106+131𝑝1980712,(40) point 2:2𝑝111+𝑝2131𝑝223=177𝑝881507𝑝442+783𝑝443201𝑝224+81𝑝8853𝑝44612.(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,𝑏2=36,𝑏3=16,𝑏4=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=12𝑖=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 are1𝑎=37,𝑏=1,𝑐=6,𝑑1=1,𝑑2=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,1𝑎=43,𝑏=2,(46)1𝑎=610,𝑏=5.(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.

fig1
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.

fig2
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=𝐔𝑛+𝐊4,(50) where the optimized weight values are 𝑤1=0,𝑤2=1/4,𝑤3=1/3,𝑤4=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𝑝𝑡𝑞=𝑞,𝑡=𝑐202𝑝.(51) The same practice was followed for the ABCs. So equations (13)-(14) become, after reintroducing 𝑐0,𝑞𝑡=±𝑐012𝑝𝑦𝑦,𝑝𝑡=𝑞±𝑝𝑥𝑐0.(52) And the PML equation (25) becomes𝑝𝑡𝑞=𝑞,𝑡+𝑍1+𝑍2𝑝𝑡+𝑍1𝑍2𝑝𝑐=20𝑝,𝐟+𝐟𝑡=𝐀1𝐟+𝑐20𝐀2𝑝.(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.

tab1
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.

tab2
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.

174031.fig.003
Figure 3: Pressure history—PML BC.
174031.fig.004
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.1𝐸=2Ω𝑐20𝜕𝑥𝑝2+𝑐20𝜕𝑦𝑝2+𝜕𝑡𝑝2𝑑Ω.(54)

174031.fig.005
Figure 5: Absorbing ability comparison of ABC and PML.
174031.fig.006
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𝑟0 with circulation intensity Γ. A schematic is presented in Figure 7. The vortices rotate around each other with a period 𝑇=8𝜋2𝑟20/Γ. Each vortex induces an angular velocity 𝑣𝜃=Γ/(4𝜋𝑟0). The configuration results in a rotating speed 𝜔=Γ/(4𝜋𝑟20) and rotating Mach number 𝑀𝑟=𝑣𝜃/𝑐0=Γ/(4𝜋𝑟0𝑐0)=2𝜋𝑟0/𝑇𝑐0.

174031.fig.007
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𝑧2𝑏12𝑧2ΓΓ𝜋𝑖ln𝑧𝑏2𝜋𝑖2𝑧2,(55) where 𝑧=𝑟𝑒𝑖𝜃 and 𝑏=𝑟0𝑒𝑖𝜔𝑡. 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.𝑢𝑥𝑖𝑢𝑦=𝜕Φ=Γ𝜕𝑧𝑧𝑖𝜋𝑧2𝑏2,𝑝(56)inc=𝑝0𝜌0𝜕(Re(Φ))1𝜕𝑡2𝜌0𝑢2𝑥+𝑢2𝑦.(57) Using (57) and differentiating twice, it is easy to calculate the acoustic sources in (7). The result is𝑝inc𝑡𝑡=𝜌0ΦRe+̇𝑢2𝑥+𝑢𝑥̈𝑢𝑥+̇𝑢2𝑦+𝑢𝑦̈𝑢𝑦,(58) where 𝜌0 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 𝑟/𝑟01.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).

tab3
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]𝑝𝜌(𝑟,𝜃,𝑡)=0Γ464𝜋3𝑟20𝑐20𝐽2(𝑘𝑟)cos(2𝜃2𝜔𝑡)𝑌2(,𝑘𝑟)sin(2𝜃2𝜔𝑡)(59) where 𝑘=2𝜔/𝑐0 and 𝐽2,𝑌2 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 𝑟/𝑟080. 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.

fig8
Figure 8: Acoustic sources (a) contour and (b) surface plot at 𝑡=120.
fig9
Figure 9: Predicted pressure fields at (a) 𝑡=120 and (b) 𝑡=280.
174031.fig.0010
Figure 10: Pressure decay along the 𝑥 axis at 𝑡=280.
174031.fig.0011
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,𝜃=900), 𝐵(𝑟=5,𝜃=1350), and 𝐶(𝑟=5,𝜃=1800) and compare them with distributions from an analytical expression.

174031.fig.0012
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.

tab4
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.

174031.fig.0013
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.𝜙(𝑟,𝜃,𝑡)=0𝑒𝐴(𝑟,𝜃,𝜔)sin(𝜔𝑡)𝑑𝜔,𝐴(𝑟,𝜃,𝑡)=𝑏𝜔/4𝐽2𝑏0𝜔𝑟2+𝑥2𝑠2𝑟𝑥𝑠+cos𝜃𝑘=0𝜀Re𝑘𝐻𝑘(𝑟𝜔)cos(𝑘𝜃)2𝑘𝜋𝐻𝑘(𝜔/2)𝜋𝜔𝐻𝑘+1(𝜔/2)𝜋0𝜔0.5𝑥𝑠cos(𝜂)0.25+𝑥2𝑠𝑥𝑠cos(𝑛)×𝐽1𝜔0.25+𝑥2𝑠𝑥𝑠,cos(𝑛)cos(𝑘𝑛)𝑑𝑛(60) where 𝐽0,𝐽1 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.

174031.fig.0014
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).

References

  1. S. K. Lele, “Compact finite difference schemes with spectral-like resolution,” Journal of Computational Physics, vol. 103, no. 1, pp. 16–42, 1992. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet · View at Scopus
  2. C. K. W. Tam and J. C. Webb, “Dispersion-relation-preserving finite difference schemes for computational acoustics,” Journal of Computational Physics, vol. 107, no. 2, pp. 262–281, 1993. View at Publisher · View at Google Scholar · View at MathSciNet · View at Scopus
  3. M. R. Visbal and D. V. Gaitonde, “On the use of higher-order finite-difference schemes on curvilinear and deforming meshes,” Journal of Computational Physics, vol. 181, no. 1, pp. 155–185, 2002. View at Publisher · View at Google Scholar · View at MathSciNet · View at Scopus
  4. C. K. W. Tam, “Computational aeroacoustics: an overview,” in Proceedings of the RTO AVT Symposium on Ageing Mechanisms and Control—Part A: Developments in Computational Aero- and Hydro-Acoustics, Manchester, UK, October 2001.
  5. C. Bogey and C. Bailly, “A family of low dispersive and low dissipative explicit schemes for flow and noise computations,” Journal of Computational Physics, vol. 194, no. 1, pp. 194–214, 2004. View at Publisher · View at Google Scholar · View at Scopus
  6. B. Engquist and A. Majda, “Absorbing boundary conditions for numerical simulation of waves,” Applied Mathematical Sciences, vol. 74, pp. 1765–1766, 1977. View at Google Scholar
  7. B. Engquist and A. Majda, “Radiation boundary conditions for acoustic and elastic wave calculations,” Communications on Pure and Applied Mathematics, vol. 32, pp. 313–357, 1979. View at Google Scholar
  8. M. B. Giles, “Nonreflecting boundary conditions for Euler equation calculations,” AIAA Journal, vol. 28, no. 12, pp. 2050–2058, 1990. View at Google Scholar · View at Scopus
  9. A. Bayliss and E. Turkel, “Radiation boundary conditions for wave like equations,” Communications on Pure and Applied Mathematics, vol. 33, pp. 707–725, 1980. View at Google Scholar
  10. C. Bogey and C. Bailly, “Three-dimensional non-reflective boundary conditions for acoustic simulations: far field formulation and validation test cases,” Acta Acustica United with Acustica, vol. 88, no. 4, pp. 463–471, 2002. View at Google Scholar · View at Scopus
  11. K. W. Thompson, “Time dependent boundary conditions for hyperbolic systems,” Journal of Computational Physics, vol. 68, no. 1, pp. 1–24, 1987. View at Google Scholar · View at Scopus
  12. K. W. Thompson, “Time-dependent boundary conditions for hyperbolic systems, II,” Journal of Computational Physics, vol. 89, no. 2, pp. 439–461, 1990. View at Publisher · View at Google Scholar · View at Scopus
  13. M. Israeli and S. A. Orszag, “Approximation of radiation boundary conditions,” Journal of Computational Physics, vol. 41, no. 1, pp. 115–135, 1981. View at Google Scholar · View at Scopus
  14. R. Kosloff and D. Kosloff, “Absorbing boundaries for wave propagation problems,” Journal of Computational Physics, vol. 63, no. 2, pp. 363–376, 1986. View at Google Scholar · View at Scopus
  15. J. P. Berenger, “A perfectly matched layer for the absorption of electromagnetic waves,” Journal of Computational Physics, vol. 114, no. 2, pp. 185–200, 1994. View at Publisher · View at Google Scholar · View at MathSciNet · View at Scopus
  16. W. C. Chew and W. H. Weedon, “3D perfectly matched medium from modified Maxwell's equations with stretched coordinates,” Microwave and Optical Technology Letters, vol. 7, no. 13, pp. 599–604, 1994. View at Google Scholar · View at Scopus
  17. F. Q. Hu, “On absorbing boundary conditions for linearized euler equations by a perfectly matched layer,” Journal of Computational Physics, vol. 129, no. 1, pp. 201–219, 1996. View at Publisher · View at Google Scholar · View at MathSciNet · View at Scopus
  18. F. Q. Hu, “A stable, perfectly matched layer for linearized Euler equations in unsplit physical variables,” Journal of Computational Physics, vol. 173, no. 2, pp. 455–480, 2001. View at Publisher · View at Google Scholar · View at MathSciNet · View at Scopus
  19. D. Kondaxakis and S. Tsangaris, “Aerodynamic and aeroacoustic analysis of a pulsating internal flow by a multidomain weak collocation spectral method,” Fluid Dynamics Research, vol. 40, no. 5, pp. 311–342, 2008. View at Publisher · View at Google Scholar
  20. J. H. Seo and R. Mittal, “A high-order immersed boundary method for acoustic wave scattering and low-Mach number flow-induced sound in complex geometries,” Journal of Computational Physics, vol. 230, pp. 1000–1019, 2011. View at Google Scholar
  21. C. Bogey, C. Bailly, and O. Marsden, “Progress in direct noise computation,” International Journal of Aeroacoustics, vol. 9, pp. 123–143, 2010. View at Google Scholar
  22. F. Q. Hu, M. Y. Hussaini, and J. L. Manthey, “Low-dissipation and low-dispersion Runge-Kutta schemes for computational acoustics,” Journal of Computational Physics, vol. 124, no. 1, pp. 177–191, 1996. View at Publisher · View at Google Scholar · View at MathSciNet
  23. M. J. Lighthill, “On sound generated aerodynamically: I. General theory,” Proceedings of the Royal Society of London, vol. 211, pp. 564–587, 1952. View at Google Scholar
  24. R. A. Renaut and J. Petersen, “Stability of wide-angle absorbing boundary conditions for the wave equation,” Geophysics, vol. 54, no. 9, pp. 1153–1163, 1989. View at Google Scholar
  25. L. N. Trefethen and L. Halpern, “Well-posedness of one-way wave equations and absorbing boundary conditions,” Mathematics of Computation, vol. 47, pp. 421–435, 1986. View at Google Scholar
  26. A. Bayliss, M. Gunzburger, and E. Turkel, “Boundary conditions for the numerical solution of elliptic equations in exterior regions,” SIAM Journal on Applied Mathematics, vol. 42, pp. 430–451, 1982. View at Google Scholar
  27. S. Abarbanel, D. Gottlieb, and J. S. Hesthaven, “Long time behavior of the perfectly matched layer equations in computational electromagnetics,” Journal of Scientific Computing, vol. 17, no. 1–4, pp. 405–422, 2002. View at Google Scholar
  28. D. Appelö and G. Kreiss, “Application of a perfectly matched layer to the nonlinear wave equation,” Wave Motion, vol. 44, no. 7-8, pp. 531–548, 2007. View at Publisher · View at Google Scholar · View at MathSciNet
  29. D. V. Gaitonde, J. S. Shang, and J. L. Young, “Practical aspects of higher-order numerical schemes for wave propagation phenomena,” International Journal for Numerical Methods in Engineering, vol. 45, no. 12, pp. 1849–1869, 1999. View at Publisher · View at Google Scholar
  30. M. Popescu, W. Shyy, and M. Garbey, “Finite volume treatment of dispersion-relation-preserving and optimized prefactored compact schemes for wave propagation,” Journal of Computational Physics, vol. 210, no. 2, pp. 705–729, 2005. View at Publisher · View at Google Scholar
  31. R. Ewert and W. Schröder, “Acoustic perturbation equations based on flow decomposition via source filtering,” Journal of Computational Physics, vol. 188, no. 2, pp. 365–398, 2003. View at Publisher · View at Google Scholar · View at MathSciNet
  32. D. J. Lee and S. K. Koo, “Numerical study of sound generation due to a spinning vortex pair,” AIAA Journal, vol. 33, no. 1, pp. 20–26, 1995. View at Google Scholar
  33. J. A. Ekaterinaris, “New formulation of Hardin-Pope equations for aeroacoustics,” AIAA Journal, vol. 37, no. 9, pp. 1033–1039, 1999. View at Google Scholar
  34. B. E. Mitchell, S. K. Lele, and P. Moin, “Direct computation of the sound from a compressible co-rotating vortex pair,” Journal of Fluid Mechanics, vol. 285, pp. 181–202, 1995. View at Google Scholar
  35. D. J. Bodony, “Analysis of sponge zones for computational aeroacoustics,” in Proceedings of the 2nd Computational Aeroacoustics (CAA) Workshop on Benchmark Problems, NASA CP-3352, 1997.