Comparison of optimized high-order finite-difference schemes for Computational Aeroacoustics

The aerospace sector is actively engaged in developing quieter aircraft to meet progressively more stringent aircraft noise legislation. The European Union has set the target of reducing aircraft noise by 10 dB over 20 years, which is driving the research agenda towards break-through aeroacoustic technology. The low-noise design methods developed by the aerospace sector area of interest to the wider engineering community, such as the automotive, marine and the railways, where low noise adds to the commercial value of their products. A key element to achieve this challenging aircraft noise target is upstreaming quiet technology in the aircraft design process. While computational aeroacoustics is a relatively new development of computational fluid dynamics, its maturity is actively encouraged by the aerospace sector, driven by its need to predict the acoustic performance of virtual prototypes from computer based design. In computational aeroacoustic methods that use the direct approach, the aerodynamic near-field fluctuations that generate noise and the near-field acoustic radiation are solved as a single problem. The hydrodynamic pressure fluctuations of the acoustically active flow can be of the order of the free stream dynamic head, while near-field noise levels down to at least 80 dB are typically required for design purposes. The challenge of computational aeroacoustics is to be able to resolve accurately on the same computational domain these perturbations that differ by four or more orders of magnitude, albeit in different regions. This requires numerical methods with low dispersion and low dissipation characteristics. The seminal work of Lele on high-order finite-difference schemes for computational aeroacoustics encouraged active work in this field, documented by four computational aeroacoustic workshops on benchmark problems and related literature. The early work focussed on achieving a high-order formal accuracy, as determined by the truncation error of the approximate Taylor series expansion of the differential operators. This was done by evaluating the coefficients of the finite-difference stencil to cancel the residual from the Taylor series expansion up to the desired order of accuracy. In later schemes, the emphasis shifted towards achieving a low dispersion and low dissipation error over the wavenumber range of interest for a specific application. In explicit time-marching methods, separate optimizations for the spatial differentiation and the time integration are performed. Pirozzoli proposed an analysis of the global discretization error that takes into account two approximately independent sources of error, associated with the discretization of space and time. Rigorous error and cost metrics have been defined for the one-dimensional linear advection


I. Introduction
The aerospace sector is actively engaged in developing quieter aircraft to meet progressively more stringent aircraft noise legislation. 1The European Union has set the target of reducing aircraft noise by 10 dB over 20 years, 2 which is driving the research agenda towards break-through aeroacoustic technology.The low-noise design methods developed by the aerospace sector area of interest to the wider engineering community, such as the automotive, marine and the railways, where low noise adds to the commercial value of their products.
A key element to achieve this challenging aircraft noise target is upstreaming quiet technology in the aircraft design process.While computational aeroacoustics is a relatively new development of computational fluid dynamics, its maturity is actively encouraged by the aerospace sector, driven by its need to predict the acoustic performance of virtual prototypes from computer based design.
In computational aeroacoustic methods that use the direct approach, the aerodynamic near-field fluctuations that generate noise and the near-field acoustic radiation are solved as a single problem.The hydrodynamic pressure fluctuations of the acoustically active flow can be of the order of the free stream dynamic head, while near-field noise levels down to at least 80 dB are typically required for design purposes.The challenge of computational aeroacoustics is to be able to resolve accurately on the same computational domain these perturbations that differ by four or more orders of magnitude, albeit in different regions. 3This requires numerical methods with low dispersion and low dissipation characteristics.
The seminal work of Lele 4 on high-order finite-difference schemes for computational aeroacoustics encouraged active work in this field, documented by four computational aeroacoustic workshops on benchmark problems [5][6][7][8] and related literature. 9The early work focussed on achieving a high-order formal accuracy, as determined by the truncation error of the approximate Taylor series expansion of the differential operators.This was done by evaluating the coefficients of the finite-difference stencil to cancel the residual from the Taylor series expansion up to the desired order of accuracy.In later schemes, 10 the emphasis shifted towards achieving a low dispersion and low dissipation error over the wavenumber range of interest for a specific application.irozzoli 12,13 proposed an analysis of the global discretization error that takes into account two approximately independent sources of error, associated with the discretization of space and time.Rigorous error and cost metrics have been defined for the one-dimensional linear advection equation, which yield appropriate choices for the grid spacing and the time step by means of a cost minimization procedure for a given level of error.This analysis led to stating relatively simple criteria for building space and time optimized computational aeroacoustic simulations.Pirozzoli 12 applied this technique to the classical penta-diagonal sixth order space scheme of Lele 4 coupled with the fourth order Runge-Kutta time integration, obtaining cost-based optimal coefficients for this scheme.
Starting from the implicit finite-difference cost-optimized stencils, the objective of this paper is to obtain analytically equivalent stencils for the prefactored finite-difference compact scheme of Hixon, 14 coupled with the classical fourth order Runge-Kutta time marching scheme.The compact notation works with shorter stencils.This reduces the extent of the near-boundary region where special boundary stencils are used. 14In multi-block, multi-processor applications, the compact formulation should reduce the depth of the computational block overlaps across inter-block boundaries, reducing the data exchange among processors.By combining the simpler near boundary data treatment of a compact scheme with a cost-optimized interior formulation, this paper contributes to developing a computational aeroacoustic scheme that is suitable to the multi-block parallel computations that are required to model the aeroacoustic of complex geometries, representative of real aerospace components.
Section II reviews the analytical findings of Ref. 12 and extends the optimization theory to the prefactored schemes.Section III presents the numerical methods used at the University of Roma 'La Sapienza' and at the University of Leicester for the present comparison.Section IV shows the validation of the optimized compact formulation against the optimized implicit multi-diagonal formulation and the reference analytical benchmark solution for the one-dimensional advection of a monochromatic sinusoidal wave.The variation of the L 2 error norm with resolution in space and time is explored analytically and numerically.

A. Efficiency analysis of finite-difference schemes
This section reviews the theory and the main results of Ref. 12.The model problem is the linear propagation (at constant phase speed a) of sinusoidal disturbances with wavelength λ (and wavenumber w = 2π/λ) in a one-dimensional unbounded domain ∂u ∂t + a ∂u ∂x = 0, u(x, 0) = û0 e i w x , uniformly discretized both in space (with grid spacing h) and time (with time step k).Ref. 12 assumed s-stage Runge-Kutta time integration and (either explicit or implicit) finite difference discretization of the spatial derivatives of the type (denoted in the following as CP p) where v j is the approximate solution of Eq. 1 at node j.
Let E be the relative L 2 error norm of the computed solution at time T = nk, and let C be the computational cost, defined as where C is assumed to be proportional to the total number of points in the computational domain (L/h), the number of operations/node required by Eq. (2) (say N op ), the number of Runge-Kutta stages (s) and the number of time steps (T /k).It is possible to derive, with a few approximations, 12 the following expressions for normalized cost and error metrics where ϕ = wh is the reduced wavenumber, σ = ak/h is the Courant number, (awT ) is a measure of the number of wavelengths traveled by the wave in a time interval T , (wL) is a measure of the number of wavelengths contained in the computational domain, and is the amplification factor of the difference scheme. 15The spatial discretization enters the analysis through the modified wavenumber Φ(ϕ), 4 defined as and time integration through the coefficients γ m , which are related to the coefficients a ij and b i of the underlying Runge-Kutta algorithm (Eq.( 7)).Linearly r-th order RK accurate schemes (i.e.r-th order accurate for linear ODEs) are obtained by setting for which ( 7) is a r-th order Taylor series expansion of the "exact" amplification factor.When dealing with multidimensional nonlinear propagation problems of broadband signals, the primary effect is to have a whole range of spatial scales, say |w| ≤ w, and propagation velocities, say |a| ≤ a (which imply ϕ = wh and σ = ak/h), and a different weight of ϕ in the cost metric (the number of points are scaled with h d , where d is the number of space dimensions).The formulas for the normalized cost and error can then be arranged as The changes with respect to the single-scale, one-dimensional case are: i) replacement of the "local" normalized error function e in Eq. ( 5) with the "global" error e, representing the maximum of e in the entire range of relevant wavenumbers and Courant numbers; and ii) replacement of ϕ 2 in Eq. ( 6) with ϕ d+1 in the cost metric.Optimizing the performance of a given scheme (i.e. for given values of s, N op ) for a given problem amounts to requiring that the computational cost be minimum for a given error level; this can be achieved by specifying a target level (say ε) for the relative error, and finding a pair of values ϕ * (ε), σ * (ε) that minimize the cost metric.To clearly illustrate the optimization procedure, Fig. 1 reports the iso-lines of the normalized global error e and the two-dimensional cost c 2 in the (ϕ, σ) plane for the C23 spatial discretization (corresponding to α 1 = 1/2, α 2 = 1/20, a 1 = 17/24, a 2 = 101/600, a 3 = 1/600), coupled with the classical fourth-order Runge-Kutta time discretization (RK4, corresponding to γ 1 = 1, γ 2 = 1/2, γ 3 = 1/6, γ 4 = 1/24).Optimality requires tangency of the two families of curves.One (sample) optimal working point is indicated with a bullet in Fig. 1.
Ref. 12 also showed that the global error can be expressed with good approximation as e(ϕ, σ) ≈ max (e s (ϕ), e t (σ ϕ)) , where e s is the spatial error (i.e. the error in the case of exact time integration), and e t is the temporal error (i.e. the error in the case of exact space integration), with z = σ ϕ, z = σ ϕ.The condition of tangency of the iso-error and iso-cost curves for a normalized error level ε is therefore realized with good approximation when The validity of the latter condition is checked in Fig. 2, where the error map reported in Fig. 1 is shown together with the corresponding approximation (12).The figure confirms that the "true" optimal points lie close to the points where Eq. ( 15) is satisfied.The condition (12) can be used to determine the approximate optimal working conditions of a given scheme for a given problem, by separately considering the contributions of space-and time-discretization, as follows: i) determine the optimal reduced wavenumber according to ii) determine the optimal Courant number from The normalized cost corresponding to the optimal working conditions is (see Eq. ( 10)) B. Optimization of finite-difference schemes Equation (18) suggests that the approximate analysis can also be exploited to develop new schemes.Indeed, optimized finite difference schemes for a specified target error level ε can be designed trying to maximize ϕ * (ε) and z * (ε) in Eq. ( 18), which amounts to separately optimize the spatial and time discretization schemes for the same error level.

Optimized space discretization
The optimization of a scheme for spatial discretization is achieved by maximizing the 'spatial resolving efficiency' ϕ * for a given value of the normalized error.For the sake of the analysis, let us consider a family of fourth-order compact difference approximations that satisfy Eqn. 2, with P = 1, p = 2, and where α 1 is a free parameter.The specific choice of α 1 = 1/3, yields the only sixth-order scheme in the family (C12).The authors have attempted to find members of this family of schemes that maximize the resolving efficiency as a function of the normalized error level ǫ.The optimal value of the coefficients α 1 is reported in tabular form for representative values of ǫ in Table 1, that reports also the data on the resolving efficiency of the (non-optimized) C12 scheme.Optimized space discretization schemes have greater values of the spatial coefficients α 1 and (α 2 ) than the non-optimized one.It is evident that optimized spatial discretizations tailored for a specific error level in principle outperform C12, yielding 25 ÷ 30% increase in the spatial resolving efficiency.

Optimized time discretization
The optimization of schemes for time integration is achieved by maximizing z * , which can be interpreted as a 'temporal resolving efficiency', for a given value of the normalized error, under the stability constraint in this analysis it is considered as representative example a four-stage, second-order Runge-Kutta scheme, i.e. set γ 1 = 1, γ 2 = 1/2, and left two free parameters γ 3 , γ 4 .The results of the analysis are reported in Table 2, which also reports the data on the resolving efficiency of the non-optimized RK4 scheme.Optimal time integration schemes have smaller values of γ 3 , γ 4 than the non-optimized Runge-Kutta scheme, which has γ 3 = 1/6, γ 4 = 1/24, and outperform RK4 yielding 20÷30% increase in the temporal resolving efficiency.

C. Extension to prefactored schemes
To derive the optimized prefactored compact schemes, the authors follow the previous work of Hixon 14 and Ashcroft. 10Let us consider the forward and backward operators v where v ′ i is the discrete approximation of the spatial derivative of the function v, at node i.The generic stencils for the forward and backward derivative operators are respectively defined as: where the coefficients must be chosen so that by adding together the two biased stencils, the original central compact scheme (see Eq. ( 2)) is recovered.To achieve this, it is important to keep in mind a peculiar property of the compact MacCormack schemes: 10,16 the real (dispersive) components of the modified wavenumbers of the forward and backward stencils are equal and identical to the modified wavenumber of the original central compact scheme, whilst the imaginary (dissipative) components of the modified wavenumbers are equal and opposite.Such property is used to develop the new class of optimized prefactored compact schemes.Combining Eq. ( 8) and Eq. ( 19), the modified wavenumber of the compact C12 scheme can be written as: 3(1+2α1) sin (ϕ) + (4α1−1) 6(1+2α1) sin (2ϕ) The modified wavenumber of the forward and backward operators may be determined in a similar manner, using the Fourier analysis.The real and imaginary component of the reduced wavenumber, of the generic forward stencil defined by Eq. ( 22), are given by: and for the backward stencil from Eq. ( 23): .
In the previous equations, by imposing the following restrictions on the coefficients of the backward stencil, it is ensured that the imaginary components of forward and backward operators are equal and opposite and also, to ensure that in the regions of zero gradient the computed derivatives vanish, the following relation is introduced Finally, by matching the various terms of Eq. ( 25) with those of Eq. ( 24), the following system of equations is obtained, considering a simplified form of the generic forward biased stencil (with a F = b F = 0), to replicate the characteristics of the scheme: Due to the quadratic term in the third of these relations, the coefficients are not uniquely determined by the above system.To choose among the range of possible coefficients, according to Ashocroft 10 and Hixon, 14 the coefficients are selected to minimize the ratio α F /β F , so that the influence of errors at the boundaries on the interior scheme is minimized.The new prefactored optimized coefficients are given in Table 3.

III. Numerical method
A. Numerical method of the University of Rome 'La Sapienza' The baseline numerical method used at the University of Rome 'La Sapienza' is the tenth-order compact pentadiagonal scheme coupled with the classical fourth order Runge-Kutta time-integration. 4,12 his scheme is referred to as the c23rk4 scheme.The spatial and temporal coefficients are given in Tables 1 and 2. The error and cost maps corresponding to c23rk4 are reported in Fig. 3a.Optimized versions of the baseline c23rk4 scheme have then been considered.The spatial discretization has been selected using the criterion of the maximization of the 'spatial resolving efficiency' φ * for a given value of normalized error.The optimized coefficients for the spatial discretization are given in Table 1 for 10 −3 ≤ ē ≤ 10 −5 .The optimization of the time integration has been achieved starting from a four-stage (s = 4) second-order Runge-Kutta scheme, with γ 1 = 1 and γ 2 = 1/2, to preserve formal second order accuracy.The two free parameters γ 3 , γ 4 have been determined so as to maximize the 'temporal resolving efficiency' z * for a given value of normalized error.The results of this optimization of the γ coefficients are reported in Table 2. Collecting together the coefficients for optimized space discretization for a given level of error with the optimized time integration ones for the same level of error, a space-and time-optimized scheme is obtained with design error of 10 −n .This class of schemes is denoted as epsmn-c23rk4.The error maps obtained with the epsm4 -c23rk4 are shown in Fig. 3b.

B. Numerical method of the University of Leicester
The numerical method used at the University of Leicester is the sixth-order prefactored tridiagonal compact scheme of Hixon, 14,17 coupled with the classical fourth order time marching scheme, labelled as prefc12rk4.The error and cost maps corresponding to prefc12rk4 are reported in Fig. 4a.Adopting the previous technique (i.e.section II C), the authors have optimized the classical prefactored prefc12rk4 for the same levels of error of section A. The new optimized prefactored scheme is labelled as epsmn-prefc12rk4.The spatial and temporal coefficients are respectively given in Table 3 and 2, and the corresponding iso-error and iso-map are shown in Fig. 4b.

A. Sinusoidal monochromatic acoustic wave
Problem A considers the monochromatic acoustic wave propagation in the positive-x direction subject to the following initial conditions: and periodic boundary conditions, on a computational domain [−10; 10] × [−10; 10], on cartesian uniformly spaced grid with h x = h y = h.The simulations have been time-advanced up to the time T = 1.Extensive simulations have been carried out by varying both n, which control the wavenumber of the disturbances, and h (the mesh spacing), noticing that for the present problem the relevant wavenumber is w = 2πn, and the speed of propagation of the disturbances is a = 1.The analytical solution and the L 2 error norm of Ref. 12 are used to measure the relative error at time T between analytical and numerical solution.Figure 5(a) and figure 5(b) show the computed normalized error, as a function of ϕ = wh, σ = ak/h, for the classical c23rk4 and the optimized epsm4 -c23rk4.The computed normalized error for the classical prefactored prefc12rk4 and for the optimized version epsm4 -prefc12rk4 is plotted respectively in Fig. 6(a) and 6(b).Excellent agreement has been observed between analytical and theoretical iso-error maps, for both schemes, in the classical and optimized version.This confirms that the results derived for the scalar, one-dimensional case, including optimal choice of mesh spacing and time step, is entirely applied to the present, single-scale vector problem.

B. Comparison of the schemes
In Fig. 7 the classical scheme c23rk4 is compared against the optimized version epsm4 -c23rk4 for a level of design error of 10 −4 .Fig. 7 shows the nominal performance of the optimized epsm4 -c23rk4 scheme in terms of optimal error, reduced wavenumber and Courant number for a given cost against that of the c23rk4 scheme.The space and time optimized scheme epsm4 -c23rk4 is able to reduce the cost by 40 − 50% in a decade around the design error level of 10 −4 , compared to the corresponding non-optimized c23rk4 scheme.This can be readily observed from Fig. 7(a) by tracing a horizontal line at the design error level ǫ = 10 −4 .The intercept of this line with the epsm4 -c23rk4 curve gives a substantial decrease of the optimal cost c compared to the intercept with the classical c23rk4 curve.The direct consequence is the increase in the range of well resolved wavenumbers φ * in Fig. 7(b), and the decrease of the optimal Courant number σ * in Fig. 7(c).Fig. 7(a) shows that the epsm4 -c23rk4 curve intercepts the c23rk4 curve at cost c ≃ 4 × 10 2 .Therefore using the epsm4 -c23rk4 scheme on very fine meshes becomes disadvantageous as the relative error ǫ is higher than for the c23rk4 scheme.By using the epsm5 -c23rk4 for such computations, the intercept point with the c23rk4 curve in Fig. 7(a) is shifted to the right so that using epsm5 -c23rk4 remains beneficial for fine mesh computations.Figure 8 shows the comparison of the classical prefactored scheme prefc12rk4 against the the optimized version epsm4-prefc12rk4 for a level of design error of 10 −4 .In this case, the space and time optimized scheme epsm4-prefc12rk4 gives a cost reduction of 60 − 70% in a decade around the design error level of 10 −4 , compared to the non-optimized prefactored version prefc12rk4.

V. Conclusions
In this paper we have extended a recently developed procedure 12 to determine the optimal performance of finite-difference numerical schemes for computational aeroacoustic applications based on quantitative normalized error/cost analysis.Theoretical results based on the one-dimensional linear propagation of monochromatic waves are described, and the procedure to develop cost optimized schemes for a given application is illustrated.Starting from the original compact finite-difference formulation, we have presented the formulation of optimized schemes in the framework of the class of prefactored finite-difference compact schemes /bin/bash: :q: command not found The comparison of the performance of the two different formulations of optimized high-order finite-difference schemes (with Runge-Kutta time integration) has been presented for a one-dimensional test case involving the linear propagation of monochromatic waves, which has confirmed that theoretical predictions immediately apply to real-world aeroacoustic problems.In particular, numerical experiments confirm that reductions of the computational errors (for given computational cost) up to 60−70% are possible.The extension of the prefactored version of the optimized schemes to multi-dimensional problems is under way.

Table 1 .
Coefficients for the spatial discretization of the classical and optimized scheme

Table 2 .
Coefficients for the time integration of the classical RK4 and optimized scheme

Table 3 .
Coefficients for the spatial discretization of the prefactored classical C12 and optimized scheme