Abstract
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.
1. 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.
A recent review of high-order method can be found in [1, 2]. Four different families of high-order methods can be encountered in the literature: spectral and pseudospectral [3], Discontinuous Galerkin (DG) [4], Weighted Essentially Nonoscillatory (WENO) [5], and optimized finite difference (FD) methods [6–9].
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–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 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 . 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 two-dimensional problems in Section 8. The conclusions are made in the last section.
2. Spatial Discretization
Consider a uniform mesh size and the values of a function at nodes . The first-order spatial derivative of the function at can be approximated by a central, -point stencil, finite difference scheme aswhere the coefficients are such as and , providing a scheme without dissipation. For standard schemes, the coefficients are fully determined from matching the Taylor series to 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 OrderFourth OrderSixth OrderEighth OrderTenth OrderTwelfth Order
For standard centered finite difference schemes, the coefficients are fully derived from the above (2)–(7). When , 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: Std7p6o: (2)–(4), . Std9p8o: (2)–(5), . Std11p10o: (2)–(6), . Std13p12o: (2)–(8), .
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 wide-band 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.
3. 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 arewhere , is the wave number, and represents the Fourier transform of function .
By comparing the two sides of (9), it is clear that the modified wave number can be written aswhere 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
4. 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 and from the numerical dispersion relation , we obtain
If , 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 (group velocity) and (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 1where is the group velocity error, is group velocity, and the phase velocity in (13) in the one-dimension wave equation is assumed to be 1 in this equation.
5. Minimization Objective Function
As shown above, a good finite difference scheme should have minimum phase velocity (dispersive) and group velocity error at the same time. 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.
Tam and Webb’s Integrated Absolute Dispersion Error
Bogey and Bailly’s Integrated Relative Dispersion Error
Combined Integrated Absolute Dispersion Error and Group Velocity Errorwhere , 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 . 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 , the combined phase and group velocity integrated error in (17) will degenerate to Tam’s absolute error.
6. 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, spectral-like optimized schemes may not be well-suited, since DRP schemes sacrifice formal order of truncation in favor of wide-band 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 . 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 , and for maximum resolution characteristics is as follows.
7-Point Schemes
Second Order. Solve (2) for where and then substitute in (19) where .
Fourth Order. Solve (2)-(3), for , where , and then substitute in (19), where .
Sixth Order. Solve (2)–(4). This gets the standard scheme.
9-Point Schemes
Second Order. Solve (2) for , where , and then substitute in (19), where .
Fourth Order. Solve (2)-(3) for , where , and then substitute in (19), where .
Sixth Order. Solve (2)–(4) for , where , and then substitute in (19), where .
Eighth Order. Solve (2)–(5). This gets the standard scheme.
11-Point Schemes
Second Order. Solve (2) for , where , and then substitute in (19), where .
Fourth Order. Solve (2)-(3) for , where , and then substitute in (19), where .
Sixth Order. Solve (2)–(4) for , where , and then substitute in (19), where .
Eighth Order. Solve (2)–(5) for , where , and then substitute in (19), where .
Tenth Order. Solve (2)–(6). This gets the standard scheme.
13-Point Schemes
Second Order. Solve (2) for , where , and then substitute in (19), where .
Fourth Order. Solve (2)-(3) for , where , and then substitute in (19), where .
Sixth Order. Solve (2)–(4) for , where , and then substitute in (19), where .
Eighth Order. Solve (2)–(5) for , where , and then substitute in (19), where .
Tenth Order. Solve (2)–(6) for , where , and then substitute in (19), where .
Twelfth Order. Solve (2)–(7). This gets the standard scheme.
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 6-order truncation error, which has two free parameters that can be optimized.
7. Optimization Ranges
in (15)–(17) is a multivariable function of , , and , where . The value of is the shading area in the Figure 1; minimization of is to minimize this shading area in the Figure 1.
At most time, we choose the integral lower limit value ; 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 , when , since ; 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 , , optimize the problem using the corresponding (2)–(6), (19) conditions of the scheme from Sections 2 and 5; get an initial range and the corresponding coefficient vector .(2)Loop in range (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 , 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 for each row; the phase velocity error is specified to in this work, which is the phase velocity error tolerance. Add to the list of obtained in previous steps; at last, we get a list of , .(4)Find the best and the corresponding in the matrix , the obtained is the best integral upper limit, and the best optimization coefficients of the finite difference scheme are .
Maximum Integral Upper Limit for Opt11p6o. Follow the algorithm above, take the 11-point 6-truncation error order scheme as a demo, we can easily find , and the corresponding and the coefficients are
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 (). From the optimization result, we get that the resolution range is , and the coefficients are
Case 2 (). From the optimization result, we get that the resolution range is , and the coefficients areModified 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 , 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.
8. 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.
8.1. 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 aswhere 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 . The time integration is numerically handled by 6-stage Runge-Kutta with a time step of and the . 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 norm error of solution is defined aswhere is the numerical solution and is the exact solution of (25).
The solution solved for the disturbance has been advected over (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 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.
8.2. 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 and . An initial pressure pulse is released at and an acoustic pulse is generated. The wave front of the acoustic pulse expanded radially. It involves the two-dimensional linearized Euler equation in the form below: wherewhere , , , 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 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 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 are shown in Figures 7 and 8.
The pressure waveform along the axis of at 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.
8.3. 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 and two cylinders with diameter 0.75 are placed at and . The spatially distributed axisymmetric acoustic source is placed at the origin, , and source parameters are set to , . The source is defined by
The computational domain is , and a uniform Cartesian grid with 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 profile along the line is compared with the analytical solution [25] in Figures 11(a) and 11(b). We note that the present results coincide well with the analytical solution thereby demonstrating the ability of the proposed scheme to provide an accurate solution with relatively complex geometries.
(a)
(b)
9. Conclusion
High-order finite difference scheme for computational aeroacoustics was presented in the last few decades [10, 11]. There is a large body contributing to this method due to its simplicity and efficiency. Yet quite a few problems still exist relevant to this method. First of all, the determination of the proper threshold of error tolerance which was used to determine the maximum resolution in the optimization procedure will be dependent on specific application problems; and the integral upper limit value affects the optimization result significantly. Besides, the choice of the minimization integral error objective functions is also varied from different authors. More recently, as Zingg [13] and Cunha and Redonnet [14] pointed, for long wave propagation, large stencil with higher truncation error order should be better than the traditional truncation error order of fourth order.
The explicit centered finite difference schemes with various truncation error orders are optimized based upon the DRP concept in this paper. An algorithm for determining the best integral upper limit in the optimization procedure is given, followed by an intuitive explanation. Three types of minimum error objective functions are discussed for the optimization procedure. It can be concluded that Tam’s integral 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.
Appendix
The coefficients of the central finite difference schemes studied and related to this work are shown in Tables 1–3. The term STD in the tables means the last columns in Tables 1–3 which is the standard maximum-order finite difference schemes fully derived from the Taylor series expansion. All the data in the appendix has been uploaded to the web [27].
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgments
This work received the support of the National Natural Science Foundation of China (No. 51175195). The authors wish to thank Dr. Jung Hee Seo for his guidance on the analytical solutions of the sound wave scattering by three cylinders. The authors would also like to thank Dr. Xiaolu Zhang (University of Southampton) for his helpful comments on the revision of this paper.