`ISRN Mechanical EngineeringVolume 2011, Article ID 174031, 15 pageshttp://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 to be constant. Inside the field satisfies 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: where   represents the Lighthill's source tensor which has the form:

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: Under the assumption that in an incompressible flow pressure fluctuations can be split into the acoustic and incompressible part , Lighthill's analogy becomes upon substitution to (4) and using(6) 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 as where is the operator The 2D analog of the exact 1D radiation condition would be written as 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). Now expanding the square root quantity in powers of ; Changing back to time domain using and , substituting to (10), and taking one more time derivative, we obtain the boundary conditions for the direction: Equivalently for direction: 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: 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 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 Another expression for circle boundaries provided in [26] is 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 .

The Laplace transform of the pressure fluctuations is The Laplace transform of (1) gives the modified Helmholtz equation for , which after introducing the coordinate transformation (20) takes the form of (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: After replacing derivatives in (21), we get Or equivalently Introducing the auxiliary vector variable , and reverting to time domain: where matrices and are defined as follows: 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 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: 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 (). 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. Using the differentiation chain rule (1) transforms to The normal derivative which appears in boundary conditions expressions transforms to 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 and therefore the normal derivative can be calculated in the transformed coordinate system simply by The expression for left/right boundaries is similar: So the wall conditions take the form:

##### 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) And for the second derivative, 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] And for the second derivative, Implementation to the curvilinear system is straightforward by taking .

##### 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: point 2: And for the 2nd derivative,

point 1: point 2: 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: Then the one-sided first derivative formulas (38)-(39) can be used to form the tridiagonal system of equations, where the value of 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 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: In (44), one may notice the presence of the derivative at the first point   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 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, 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.

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]: If a component of the solution vector is denoted by , filtered values are obtained by solving again a tridiagonal system. (48) corresponds to a -order filter and a 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  , 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  . The absolute value of the complex amplification factor represents the dissipation errors, whereas dispersion errors are represented by the phase.

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: update: where the optimized weight values are [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 The same practice was followed for the ABCs. So equations (13)-(14) become, after reintroducing , And the PML equation (25) becomes 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 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 .

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.

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 with circulation intensity . A schematic is presented in Figure 7. The vortices rotate around each other with a period . Each vortex induces an angular velocity . The configuration results in a rotating speed and rotating Mach number .

Figure 7: Spinning vortex pair schematic.

The incompressible inviscid flow field can de determined numerically [31] by the evaluation of the complex potential function : where and . The approximation in (55) is valid only in the far field that is for . 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. Using (57) and differentiating twice, it is easy to calculate the acoustic sources in (7). The result is where is the fluid’s density and 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 , 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 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] where 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 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 and 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 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 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 . 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.

Figure 8: Acoustic sources (a) contour and (b) surface plot at .
Figure 9: Predicted pressure fields at (a) and (b) .
Figure 10: Pressure decay along the axis at .
Figure 11: Pressure time series at .
##### 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 , and the noise source is located at . The sound source can be either periodic, or nonperiodic where in the later case a pulse is released at and the initial pressure distribution is given as . One must find the pressure time series at points , , and 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 . 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 . 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 .

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. where are the zeroth and first-order Bessel functions of the first kind and is the th order Hankel function of the first kind. for and for , and is the source location in the axis. The pressure can then be found as the first time derivative: Care must be taken when numerically evaluating (60)-(61), since Bessel and Hankel functions go to infinity at zero. A value of 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).

#### 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.
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.
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.
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.
6. B. Engquist and A. Majda, “Absorbing boundary conditions for numerical simulation of waves,” Applied Mathematical Sciences, vol. 74, pp. 1765–1766, 1977.
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.
8. M. B. Giles, “Nonreflecting boundary conditions for Euler equation calculations,” AIAA Journal, vol. 28, no. 12, pp. 2050–2058, 1990.
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.
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.
11. K. W. Thompson, “Time dependent boundary conditions for hyperbolic systems,” Journal of Computational Physics, vol. 68, no. 1, pp. 1–24, 1987.
12. K. W. Thompson, “Time-dependent boundary conditions for hyperbolic systems, II,” Journal of Computational Physics, vol. 89, no. 2, pp. 439–461, 1990.
13. M. Israeli and S. A. Orszag, “Approximation of radiation boundary conditions,” Journal of Computational Physics, vol. 41, no. 1, pp. 115–135, 1981.
14. R. Kosloff and D. Kosloff, “Absorbing boundaries for wave propagation problems,” Journal of Computational Physics, vol. 63, no. 2, pp. 363–376, 1986.
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.
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.
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.
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.
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.
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.
21. C. Bogey, C. Bailly, and O. Marsden, “Progress in direct noise computation,” International Journal of Aeroacoustics, vol. 9, pp. 123–143, 2010.
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.
23. M. J. Lighthill, “On sound generated aerodynamically: I. General theory,” Proceedings of the Royal Society of London, vol. 211, pp. 564–587, 1952.
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.
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.
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.
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.
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.
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.
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.
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.
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.
33. J. A. Ekaterinaris, “New formulation of Hardin-Pope equations for aeroacoustics,” AIAA Journal, vol. 37, no. 9, pp. 1033–1039, 1999.
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.
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.