Explicit High Accuracy Maximum Resolution Dispersion Relation Preserving Schemes for Computational Aeroacoustics

A set of explicit finite difference schemes with large stencil was optimized to obtain maximum resolution characteristics for various spatial truncation orders. The effect of integral interval range of the objective function on the optimized schemes’ performance is discussed. An algorithm is developed for the automatic determination of this integral interval. Three types of objective functions in the optimization procedure are compared in detail, which show that Tam’s objective function gets the best resolution in explicit centered finite difference scheme. Actual performances of the proposed optimized schemes are demonstrated by numerical simulation of three CAA benchmark problems.The effective accuracy, strengths, and weakness of these proposed schemes are then discussed. At the end, general conclusion on how to choose optimization objective function and optimization ranges is drawn.The results provide clear understanding of the relative effective accuracy of the various truncation orders, especially the trade-off when using large stencil with a relatively high truncation order.


Introduction
In a wide range of industrial fields such as aircraft, automobile, and turbomachinery flow induced noise and noise often interact with the turbulent flow strongly.These aeroacoustic problems usually involve sound waves of many octave bands, due to the large disparity in the characteristic scales of the acoustic and the flow fluctuations; high-order accurate methods with minimal dissipation and dispersion errors are required.
Finite difference schemes are considered in the present work.Owing to the simplicity, a significant body of research has been devoted to optimize the finite difference method.
Several finite difference schemes [6,[10][11][12] have been designed or optimized to resolve the waves with high accuracy and good spectral resolution.Lele [10] constructed such schemes by constraining some of the coefficients in the compact finite difference to achieve a certain truncation error and then determining the remainder of the coefficients by requiring the modified wave number to be equal to exact differentiation at certain wave numbers.Following a similar procedure, Tam and Webb [11] increased the resolution of a central finite difference approximation by minimizing the integrated dispersive (phase) errors in the wave number domain and proposed the dispersion-relation-preserving (DRP) scheme.The coefficients of the central discretization are determined by the truncation order and the minimization of errors.The fourth-order spatial central scheme of an optimized seven-point stencil shows better resolution characteristics than the standard maximum-order central finite difference schemes.However, as Zingg [13] pointed out, if the schemes are optimized too aggressively, they perform poorly 2 Mathematical Problems in Engineering for longer distances of travel.Furthermore, if a waveform has significant low wave number content, as in the case of a Gaussian pulse, the benefits of an optimized scheme can be minimal, and maximum-order schemes can even be superior.Cunha and Redonnet [14] also expressed a similar opinion on comparison of the optimization scheme and standard maximum-order scheme.
For the optimization of the finite difference scheme, there are two key factors that can affect the optimization result mostly.The first one is the optimization ranges for the integral modified wave number error.In Tam and Webb's [11] first paper, this integral range is arbitrarily determined to [0, /2].As Lockard et al. [15] pointed out, this was an aggressive choice.Later, in another paper by Tam et al. [8], they gave a more reasonable optimization range.However, they did not explain these values in detail.In Kim and Lee's [16] paper on optimized compact scheme, this optimization range is determined by a trial and error method.In Cheong and Lee's [17] paper on grid-optimized DRP scheme on General Geometries for Computational Aeroacoustics they just simply follow Tam's integral ranges.More recently, Liu et al. [18] develop an algorithm to determine this range, but it is a little complex and did not give intuitive explanation; in this work, we will give a simpler algorithm for the optimization ranges; intuitive explanation is followed in discussion part.This simple algorithm can be readily adopted for optimized implicit compact scheme, nonuniform grid DRP scheme, upwind optimized DRP scheme, and curvilinear version of the above scheme too.
The second one is the truncation error order.Tam and Webb [11] choose the 4-order formal accuracy order for 7-point central scheme and larger stencil finite difference scheme too.However, as Zingg [13] and Cunha and Redonnet [14] pointed out, this scheme was not a good scheme for long distance travel wave; it is not better than the maximum-order schemes if there are many low frequency contents.So using a larger stencil scheme combined with higher truncation error order will be a natural compromised way.Kim and Lee [16] consider various truncation error orders for compact finite schemes.The similar various truncation error orders for explicit central finite difference schemes will be considered in this work.
The organization of the paper is as follows.In Section 2, brief introduction to spatial discretization of the central finite difference schemes is recalled; phase velocity and group velocity errors of the finite difference schemes in the wave number domain are given via Fourier analysis in Sections 3 and 4. Section 5 defines three types of minimum integral error objective functions for optimization.Section 6 describes the optimization strategies and shows a detailed procedure of optimization of the various truncated error order schemes.An intuitive simple algorithm was proposed in Section 7 to achieve a reasonable integral range.Additional discussion on optimization objective function is given in Section 7 as well.Then, the effectiveness of the various error orders optimized schemes is applied to one-and twodimensional problems in Section 8.The conclusions are made in the last section.

Spatial Discretization
Consider a uniform mesh size Δ and the values of a function   = (  ) at nodes  (=0,1,...,) .The first-order spatial derivative of the function / at   can be approximated by a central, (2 + 1)-point stencil, finite difference scheme as where the coefficients   are such as  0 = 0 and  − = −  , providing a scheme without dissipation.For standard schemes, the coefficients   are fully determined from matching the Taylor series to 2 order accuracy, so that the maximum order is reached.
The relations between the coefficients of   are derived by matching the Taylor series coefficients of various truncation orders.These relations are as follows.
Second Order ( Fourth Order Sixth Order Eighth Order Tenth Order Twelfth Order For standard centered finite difference schemes, the coefficients are fully derived from the above (2)- (7).When  = 3, 4, 5, 6, we get standard schemes using 7, 9, 11, and 13 points, which are referred to as Std7p6o, Std9p8o, Std11p10o, and Std13p12o in this work and are of orders 6, 8, 10, and 12th, respectively.The Taylor series conditions for these schemes are listed as following: For convenience, their coefficients are listed in Tables 1-3 (the last column) of Appendix.These standard schemes have unique coefficients, and these are the highest order ones obtainable with these schemes.
For optimized dispersion-relation-preserving schemes, which sacrifice formal order of truncation in favor of wideband performance, they can provide significantly better wave propagation characteristics in the high wave number range, that is, high resolution.Following Tam and Webb's [11] DRP way, the other lower-order schemes must have free coefficients that are not determined until more constraints are imposed and these can be used to improve the resolution characteristics.The additional constraint conditions for the coefficients will be discussed later in Section 6.The nature of these constraints is the minimization of dispersive (phase) errors and group velocity in the wave number domain.

Dispersion Error
The finite difference equation ( 1) is of a central difference; it has no dissipative errors.In this section, the differencing errors of (1) are analyzed in terms of the dispersive errors.A Fourier analysis provides an effective tool to quantify the dispersive errors and resolution characteristics of a differencing approximation and so this quantification will be used further to guide the optimization of the differencing scheme.
Setting  = Δ, take the Fourier transform of the left and right sides of (1) that are where  = √ −1,  is the wave number, and f represents the Fourier transform of function .
By comparing the two sides of ( 9), it is clear that the modified wave number can be written as where  is the modified (effective) wave number and Δ represents the exact wave number.The numerical dispersive errors come from the deviation between the effective and exact wave number, which can be defined as

Group Velocity Error
For the DRP (dispersion-relation-preserving) scheme, the dispersion relation of the linear one-dimension wave equation can be written as [19]  () =  ( ()) .
Assuming that waves propagate in the -direction only, the group velocity of the wave is With an acceptable approximation / ≃ 1 and from the numerical dispersion relation / ≃ , we obtain If / = 1, the scheme is then to reproduce the same group velocity of the original partial differential equation.Finally, the requirements on the numerical method for correct predictions should be / = 1 (group velocity) and / = 1 (phase velocity), so combining the integrated absolute dispersion error and group velocity error as the objective function for optimized finite difference is a natural way, which will be discussion in Section 5.
The group velocity error can be given as the difference between the group velocity and the value 1 where   is the group velocity error, V  is group velocity, and the phase velocity  in (13) in the one-dimension wave equation is assumed to be 1 in this equation.

Minimization Objective Function
As shown above, a good finite difference scheme should have minimum phase velocity (dispersive) and group velocity where   ,   are defined in ( 13), (14).(Δ)/(Δ) are the group velocity, which should be close to 1; the  term is a weighting function that allows emphasis to be put on either the phase velocity or group velocity component. can be function of Δ or constant; here, for simplicity,  is set to be constant in range 0 <  < 1.This form was essentially presented by Tam and Webb [11] in the context of phase velocity error but is equally valid group velocity error.Then, (Δ)/(Δ) can be derived from (11) easily The second type (Bogey and Bailly's) error objective function will not be considered in this work since the relative error gives more weight on the high wave number, which is not a good choice for many cases.When  = 1, the combined phase and group velocity integrated error in (17) will degenerate to Tam's absolute error.

Optimization Procedure
Just as Cunha and Redonnet [14] and Zingg [13] argued, for a very low error tolerance or a long transfer distance, especially when there are many low frequency components, spectrallike optimized schemes may not be well-suited, since DRP schemes sacrifice formal order of truncation in favor of wideband performance.They are intended to more adequately represent high-frequency waves at the expense of allowing more error in low-frequency content.
To solve this problem, we can use standard schemes as Cunha and Redonnet [14] suggested; another compromise method is to use a larger stencil with the truncation error order more than 4.
In this paper, analytic and systematic constraints for determination of the free coefficients for various truncation error orders are considered.The nature of these constraints is the minimization of dispersive (phase) errors in the wave number domain, that is, the stencil wave number optimization introduced by Tam and Webb [11].
As initially proposed in [11], spatial derivative operators using explicit centered finite difference schemes can be optimized so that they have minimum integrated dispersion errors.The optimization method is to combine the traditional truncated Taylor series finite difference approximation and the wave number space approximation.
The integrated error in ( 15) is a function of the coefficients   ,  = 1, 2, . . ., .It is necessary to find the optimum values of the coefficients that minimize the integrated error.The conditions that make  a local minimum value are proposed by Tam and Webb [11] as follows: where the error  can be one of ( 15)- (17).For simplicity, we choose (15) as a demo to get the following conditions.Equations ( 2)-( 6) and (19) provide a system of linear algebraic equations by which the optimum coefficients can be determined.The optimization procedure to determine the coefficients   ,  = 1, 2, . . ., , and (Δ) ℎ for maximum resolution characteristics is as follows.
For convenience, the schemes are referred to as Opt11p6o and so on in this work, where Opt means that this is a DRP optimization scheme, the first number is the point number of the scheme, and the second number is the truncation order of the scheme.Therefore, the name Opt11p6o refers to a DRP optimized scheme based on an 11-point spatial scheme of 6order truncation error, which has two free parameters that can be optimized.
At most time, we choose the integral lower limit value (Δ)  = 0; however, the upper limit value of (Δ) ℎ is very important; this section will discuss this key term further, using 11-point 6-truncation order scheme (Opt11p6o) as a demo.The value of this upper limit should be in the range [0, ], when Δ = , Δ = 0 since sin() = 0,  = 1, 2, . . ., ; this means that the modified wave number defined in ( 9) always falls to a value of zero when the exact wave number is .Thus, most dispersive errors exist in the high wave number ranges, and it is preferable that this wave number range should be omitted in the minimum objective function in ( 15)- (17).
In this section, we give a simple algorithm to determine this upper limit value in ( 15)- (17).We take a small step through the entire ranges from 0 to  since the optimization is not time consuming; at first, we will specify an error to get the well-resolved domain; then we can find a maximum resolution ranges and the corresponding integral upper limit value (Δ) ℎ ; in this meaning, we get a maximum resolution for the schemes.
Based on the procedure mentioned above, a simple algorithm is summarized to implement the proposed method.In the algorithm,  = Δ,  ℎ * = Δ.
(1) Setting  =  0 ,  1 =  0 + , optimize the problem using the corresponding (2)-( 6), (19) conditions of the scheme from Sections 2 and 5; get an initial range  1 and the corresponding coefficient vector (2) Loop in range 0 <  < 2 (2 is enough for the explicit central finite difference scheme upper to 13 points), get a matrix of ranges and coefficients where  is the number of   in ranges [0, 2],  is the number of the coefficients, and each row represents an integral upper limit value and the corresponding coefficients at that integral upper limit.(3) From the matrix   , find  ℎ * , with |Δ − Δ| <  ph for each row; the phase velocity error  ph is specified to 1 × 10 The  ℎ * , changes as   =   Δ increases; this is shown in Figure 2. At first, the best resolution range  ℎ * , increases linearly as the integral upper limit  ℎ , until  ℎ reaches the value 1.37 and  ℎ * , reaches the maximum value 1.35.The difference between the  ℎ * , and  ℎ is 0.02, which is the length we adapt to step through the entire ranges.Then there is a big jump in the curve, which shows that the phase velocity error will increase quickly in the high frequency domain  when using aggressive integral ranges.These explained why the proper integral upper limit is so important.

Comparison of Minimization Objective
Function.Now we turn back to the discussion on error objective function.As stated in Section 5, we will only compare the Tam's error objection function and the combined error objective function.We again take Opt11p6o scheme as a demo.Two values of  are adapted in the comparison.
Using the optimization procedure in Section 6, Tam's minimum objective function, the coefficients for the two cases are list as below.
Case 1 ( = 0.5).From the optimization result, we get that the resolution range is  max = 1.27, and the coefficients are Case 2 ( = 0.1).From the optimization result, we get that the resolution range is  max = 1.32, and the coefficients are Modified wave numbers of the Opt11p6o-0.5, Opt11p6o-0.1, and Opt11p6o are shown in Figure 3, and the error is shown in Figure 4. From the figure, it is shown that when  → 0, the phase velocity behavior will tend to the tam's error function case, and using the combined error function will get a conservative optimization.The resolution will be less than Tam's objective function cases.When using combined error as minimization objective function and emphasis on group velocity, the optimization result will be conservative since the resolution decreases as the weight function  increases.However, we believe this may be investigate further for other schemes; in fact, in paper of Venutelli [21], he chooses variations of the phase speed with the wave number for compact scheme optimization, which is directly related to the group velocity error.

Illustration and Discussion
In this section, three problems were considered to investigate the performances of the proposed various truncation error orders finite difference scheme.The first problem involved the one-dimensional scalar wave convection, the second one was concerned about two-dimensional acoustic pulse convection, and the last benchmark problem is two-dimensional sound wave scattering by multiple cylinders.During the computation, various truncation error order schemes were implemented for the spatial discretization, and the six-stage Runge-Kutta optimized scheme of Hu et al. [22] is applied to temporal integration.

Convection of a Harmonic Signal.
To illustrate the strengths of the large stencil with a higher truncation order, the one-dimensional convection equation was considered; the solution consists of an initial disturbance which travels along the -axis at a unitary dimensionless speed: The initial condition is given as where  is the amplitude of the disturbance and Δ is the wavelength.The latter one, which equals the number of points-per-wavelength (PPW), is here considered Δ ∈ [4,16].The time integration is numerically handled by 6stage Runge-Kutta with a time step of Δ = 0.1 and the CFL = Δ/Δ = 0.1.The temporal approximation induces only negligible errors with a low CFL value.The total error is approximately equal to the spatial discrete differential operators.
For quantitative comparison, the  2 norm error of solution is defined as where  num is the numerical solution and  exact is the exact solution of (25).
The solution solved for the disturbance has been advected over 50Δ (50 PPW).With the increase of the wavelength (or, equivalently, PPW), the error level of the DRP schemes with truncation error less than or equal to 4 (Opt11p4o) did not decrease (Figure 5).On the contrary, the error level of their standard counterparts (Opt7p6o, Opt9p8o, and Opt11p10o) decayed rapidly due to the high truncation error order.And the DRP schemes with truncation order greater than 4 (Opt11p6o and Opt11p8o) also have similar error level to the standard counterparts.With the comparison of Opt11p8o and Opt9p8o, when the PPW is larger than about 6, these two schemes have almost the same decaying ratio.Furthermore, the error of Opt9p8o is always less than the error of Opt9p8o, especially in high frequency range.This indicates that Opt11p8o has a better accuracy for short wave components.
There are two key values of PPW in Figure 5, which are PPW = 7.6 and PPW = 9.2.Before the value of 7.6 in Figure 5, Opt11p4o has the lowest error, which is about 1e-2 for a convected distance of 50 PPW, while Std11p10o has the quickest convergence, yet the greatest error in these ranges.After this value, the performance of these three schemes is inverted.Std11p10o has the minimum error level; however, the behavior of Opt11p8o scheme is also good in this 50 PPW distance.It is found that one can conclude that, provided that 1% is set as the acceptable error level, the Opt11p8o has similar accuracy with the maximum-order Std11p10o scheme in low and mid frequency ranges but has a marginal advantage in the higher frequency range; this could increase the robust of the scheme when applied to real problems.This indicates that the Opt11p8o scheme is better than the standard counterpart Std11p10o scheme.
Before the second value PPW = 9.2, the standard Std9p8o has the largest error, while Opt11p6o and Opt11p4o have the best accuracy in this high frequency range.However, after this value, the error level of Opt11p4o does not decay substantially.Thus, it can be concluded that Opt11p6o scheme has better accuracy in a wider range when comparing to the Opt11p4o and Std9p8o schemes.
Figure 6 shows the  2 norm error of solution of (24) decreased with the increase of the truncation error order.It has been shown that the Std11p10o scheme has the least error, while the Opt11p4o has the greatest error.The Opt11p8o and Std11p10o have very little error, showing the effectiveness of this scheme.

Two Dimension Pulse.
A two-dimension acoustic pulse convection is simulated to show the effectiveness of the proposed various truncation error order schemes presented in this work.The computation domains are −50 ≤  ≤ 50 and 0 ≤  ≤ 110.An initial pressure pulse is released at  = 0 and an acoustic pulse is generated.The wave front of the acoustic pulse expanded radially.It involves the twodimensional linearized Euler equation in the form below: where where , , V, and  are the density, velocities, and pressure, respectively.  and   are set as 0 in the computation, which are the Mach number of the mean flow in the and -directions.The initial conditions are given as The simulation is done over a domain with uniform grids (Δ = Δ = 1) and the value of  is set as 0.01.CFL number of 0.1 is employed.The spatial derivatives are discretized using Opt11p6o scheme.During the simulation, Perfectly Matched Layers (PML) of 10 grid points are applied to absorb outgoing waves in the left, the right, and the upper of computation domain [23].Wall boundary V = 0 is applied at the lower domain [24].In order to keep the solution stable, filters are used to damp the spurious waves of the four variables of  after the final stage of the 6-stage Runge-Kutta time integration.Each time the filter is firstly applied in direction and then in -direction.
The pressure contours at  = 50 are shown in Figures 7  and 8.
The pressure waveform along the axis of  = 0 at  = 50 is shown in Figure 8.The density waveform of the scheme Opt11p6o is approximately equal to the exact curve.
For a quantitative comparison, the relative rate of different schemes' numerical result along the -axis is defined as The relative rate of Opt11p6o is 0.0342 and that of Std7p6o is 0.1231.This demonstrates the better performance of Opt11p6o scheme.

Acoustic Wave Scattering by Multiple Cylinders.
The last benchmark problem is two-dimensional sound wave scattering in complex geometry.This problem simulates a sound field scattered by three rigid circular cylinders of different sizes from a monopole acoustic source located between the cylinders in a zero mean flow [25,26].The geometry configuration is shown in Figure 9.The first cylinder with diameter 1 is placed at (−3, 0) and two cylinders with diameter 0.75 are placed at (3, −4) and (3,4).The spatially distributed axisymmetric acoustic source is placed (31) The computational domain is −9 ≤ (,) ≤ 9, and a uniform Cartesian grid with Δ = 0.02 is used.In the computation, the Opt11p6o scheme is again used for spatial derivatives; temporal derivatives are discretized by 6-stage LDDRK [22].
Figure 10 shows the instantaneous fluctuating pressure field, and   rms profile along the  = 0 line is compared error is the proper objective function for explicit scheme in most cases.Yet further study is suggested for other schemes.The performances of the proposed various truncation orders finite difference schemes are demonstrated when applying to both one and two-dimensional test problems.

Figure 1 :
Figure 1: Optimization region for central finite difference schemes.

Figure 3 :
Figure 3: Modified wave number of Opt11p6o and Opt11p6o optimized with combined integral error objective function.

Figure 4 :
Figure 4: Error of modified wave number of Opt11p6o and Opt11p6o optimized with combined integral error objective function.

Figure 5 :Figure 6 :
Figure 5: Maximum  2 -norm error made on disturbance by various truncation error order finite difference schemes, for a propagation distance of 50 PPW.

Table 1 :
Optimized coefficients for 7 point schemes.

Table 2 :
Optimized coefficients for 9 point schemes.

Table 3 :
[20]mized coefficients for 11 point schemes.There are two common error definitions in literatures as the optimization objective function to get a good finite difference scheme.These are Tam and Webb's[11]integrated absolute dispersion error definition and Bogey and Bailly's[20]integrated relative dispersion error.In this work, a third type error, which combines the group velocity error and the phase error, is also proposed for comparison.