Table of Contents Author Guidelines Submit a Manuscript
Journal of Applied Mathematics
VolumeΒ 2012, Article IDΒ 783101, 30 pages
http://dx.doi.org/10.1155/2012/783101
Research Article

The Technique of MIEELDLD in Computational Aeroacoustics

Department of Mathematics and Applied Mathematics, University of Pretoria, Pretoria 0002, South Africa

Received 3 January 2012; Accepted 14 February 2012

Academic Editor: M. F.Β El-Amin

Copyright Β© 2012 A. R. Appadu. 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

The numerical simulation of aeroacoustic phenomena requires high-order accurate numerical schemes with low dispersion and low dissipation errors. A technique has recently been devised in a Computational Fluid Dynamics framework which enables optimal parameters to be chosen so as to better control the grade and balance of dispersion and dissipation in numerical schemes (Appadu and Dauhoo, 2011; Appadu, 2012a; Appadu, 2012b; Appadu, 2012c). This technique has been baptised as the Minimized Integrated Exponential Error for Low Dispersion and Low Dissipation (MIEELDLD) and has successfully been applied to numerical schemes discretising the 1-D, 2-D, and 3-D advection equations. In this paper, we extend the technique of MIEELDLD to the field of computational aeroacoustics and have been able to construct high-order methods with Low Dispersion and Low Dissipation properties which approximate the 1-D linear advection equation. Modifications to the spatial discretization schemes designed by Tam and Webb (1993), Lockard et al. (1995), Zingg et al. (1996), Zhuang and Chen (2002), and Bogey and Bailly (2004) have been obtained, and also a modification to the temporal scheme developed by Tam et al. (1993) has been obtained. These novel methods obtained using MIEELDLD have in general better dispersive properties as compared to the existing optimised methods.

1. Introduction

Computational aeroacoustics (CAA) has been given increased interest because of the need to better control noise levels from aircrafts, trains, and cars due to increased transport and stricter regulations from authorities [1]. Other applications of CAA are in the simulation of sound propagation in the atmosphere to the improved design of musical instruments.

In computational aeroacoustics, the accurate prediction of the generation of sound is demanding due to the requirement for preservation of the shape and frequency of wave propagation and generation. It is well known [2, 3] that, in order to conduct satisfactory computational aeroacoustics, numerical methods must generate the least possible dispersion and dissipation errors. In general, higher order schemes would be more suitable for CAA than the lower-order schemes since, overall, the former are less dissipative [4]. This is the reason why higher-order spatial discretisation schemes have gained considerable interest in computational aeroacoustics.

The field of CAA has grown rapidly during the last decade and there has been a resurgence of interest in aeroacoustic phenomena characterised by harsher legislation and increasing environmental awareness. CAA is concerned with the accurate numerical prediction of aerodynamically generated noise as well as its propagation and far-field characteristics. CAA involves mainly the development of numerical methods which approximate derivatives in a way that better preserves the physics of wave propagation unlike typical aerodynamic computations [1].

Since multidimensional finite-volume algorithms are generally more expensive in terms of numerical cost than finite-difference algorithms, the majority of CAA codes are based on the finite-difference methodology [5]. The trend within the field of CAA has been to employ higher-order accurate numerical schemes that have in some manner been optimized for wave propagation in order to reduce the number of grid points per wavelength while ensuring tolerable levels of numerical error. Apart from acoustics and aeroacoustics, low amplitude wave propagation takes place over distances characterized by large multiples of wavelength in other areas such as [6](1)electromagnetics, for microcircuit design,(2)elastodynamics, for nondestructive testing,(3)seismology, for oil exploration,(4)medical Imaging, for accurate diagnosis,(5)hyperthermia, for noninvasive surgery.

Aerodynamics and other areas of fluid mechanics have benefitted immensely from the development of CFD [7]. The numerical analysis of flows around full aircraft configurations has become feasible with advances in both numerical techniques and computing machines. The temptation to apply effective CFD methods to aeroacoustic problems has been unavoidable and has been met with some success, but, in some cases, it has been observed that there is a necessity for some numerical protocols specific to problems involving disturbance propagation over long distances. The difference between aerodynamic and aeroacoustic problems lies mainly in the fact that for aeroacoustic computations, the solution is desired at some large distance from the aerodynamic source, but, in the case of aerodynamic problems, flow properties are required accurately only on the body itself [7]. Most aerodynamics problems are time independent, whereas aeroacoustics problems are, by definition, time dependent [8]. There are computational issues that are unique to aeroacoustics. Thus, computational aeroacoustics requires somewhat independent thinking and development.

At specific Courant numbers and angles of propagation, the perfect-shift property can be obtained, leading to exact propagation for all wavenumbers [9]. The perfect shift property refers to the situation when the error from the spatial discretisation precisely cancels that from the temporal discretisation. Several numerical schemes which combine the spatial and temporal discretisation produce the perfect shift property at specific Courant numbers [9]. Often this perfect cancellation of temporal and spatial errors occurs at cfl = 1.0. For such methods, the sum of the spatial and temporal error increases when the cfl number is decreased as the temporal error no longer cancels the spatial error. As the cfl number tends to zero, so does the temporal error and thus only spatial error remains. For most schemes, a low cfl represents the worst case associated with large dispersion or large dissipation errors as there is no cancellation of temporal and spatial errors [9]. Thus it is important to assess numerical methods over a range of Courant numbers [9]. However, this is not an issue for schemes built up from a high-accuracy spatial discretisation with a high-accuracy time-marching method. These schemes generally do not rely on cancellation to achieve high accuracy and thus the error does not increase as the Courant number is reduced.

The imaginary part of the numerical wavenumber represents numerical dissipation only when it is negative [10]. Due to the difference between the physical and numerical wavenumber, some wavenumbers propagate faster or slower than the wave speed of the original partial differential equation [11]. This is how dispersion errors are induced. The real part of the modified wavenumber determines the dispersive error while the imaginary part determines the dissipative error [9]. The group velocity of a wavepacket governs the propagation of energy of the wavepacket. The group velocity is characterised by Re((𝑑/(𝑑(πœƒβ„Ž)))(πœƒβˆ—β„Ž))βˆ’1.0 which must be almost one in order to reproduce exact result [12]. Otherwise, dispersive patterns appear. When Re((𝑑/(𝑑(πœƒβ„Ž)))(πœƒβˆ—β„Ž))=1.0, the scheme has the same group velocity or speed of sound as the original governing equations and the numerical waves are propagated with correct wave speeds.

2. Organisation of Paper

This paper is organised as follows. In Section 3, we briefly describe the technique of Minimised Integrated Exponential Error for Low Dispersion and Low Dissipation (MIEELDLD) when used to optimise parameters in numerical methods. We also describe how this technique can be extended to construct high order methods with low dispersive and low dissipative properties in computational aeroacoustics. In Sections 4–8, we use MIEELDLD to obtain some optimized spatial methods based on a modification of the methods constructed by Tam and Webb [3], Lockard et al. [13], Zingg et al. [14], Zhuang and Chen [15], Bogey and Bailly [16]. Section 9 introduces an optimised temporal scheme which is obtained using MIEELDLD and based on a modification of the temporal discretisation method constructed by Tam et al. [17]. In Section 10, we construct numerical methods based on blending each of the five new spatial schemes with the new time discretisation scheme when used to discretise the 1-D linear advection equation and obtain rough estimates of the range of stability of these methods. Section 12 highlights the salient features of the paper.

3. The Concept of Minimised Integrated Exponential Error for Low Dispersion and Low Dissipation

In this section, we describe briefly the technique of Minimized Integrated Exponential Error for Low Dispersion and Low Dissipation (MIEELDLD). This technique have been introduced in Appadu and Dauhoo [18] and Appadu and Dauhoo [19]. We now give a resume of how we have derived this technique of optimisation.

Suppose the amplification factor of the numerical scheme when applied to the 1-D linear advection equation:πœ•π‘’πœ•π‘‘+π›½πœ•π‘’πœ•π‘₯=0(3.1) isπœ‰=𝐴+𝐼𝐡.(3.2) Then the modulus of the Amplification Factor (AFM) and the relative phase error (RPE) are calculated as||πœ‰||,1AFM=RPE=βˆ’π‘Ÿπ‘€tanβˆ’1𝐡𝐴,(3.3) where π‘Ÿ and 𝑀 are the cfl number and phase angle, respectively.

For a scheme to have Low Dispersion and Low Dissipation, we require||||1βˆ’RPE+(1βˆ’AFM)⟢0.(3.4)

The quantity, |1βˆ’RPE| measures dispersion error while (1βˆ’AFM) measures dissipation error. Also when dissipation neutralises dispersion optimally, we have||||||||1βˆ’RPEβˆ’(1βˆ’AFM)⟢0.(3.5) Thus on combining these two conditions, we get the following condition necessary for dissipation to neutralise dispersion and for low dispersion and low dissipation character to be satisfied:||||||βˆ’||+ξ€·||||+ξ€Έeldld=1βˆ’RPE(1βˆ’AFM)1βˆ’RPE(1βˆ’AFM)⟢0.(3.6)

Similarly, we expectξ€·||||||βˆ’||ξ€Έξ€·||||+ξ€Έeeldld=exp1βˆ’RPE(1βˆ’AFM)+exp1βˆ’RPE(1βˆ’AFM)βˆ’2⟢0,(3.7) in order for Low Dispersion and Low Dissipation properties to be achieved.

The measure, eeldld, denotes the exponential error for low dispersion and low dissipation. The reasons why we prefer eeldld over eldld is because the former is more sensitive to perturbations.

We next explain how the integration process is performed in order to obtain the optimal parameter(s).

Only One Parameter Involved
If the cfl number is the only parameter, we compute ξ€œπ‘€10eeldld𝑑𝑀,(3.8) for a range of π‘€βˆˆ[0,𝑀1], and this integral will be a function of π‘Ÿ. The optimal cfl is the one at which the integral quantity is closest to zero.

Two Parameters Involved
We next consider a case where two parameters are involved and whereby we would like to optimise these two parameters.
Suppose we want to obtain an improved version of the Fromm’s scheme which is made up of a linear combination of Lax-Wendroff (LW) and Beam-Warming (BW) schemes. Suppose we apply BW and LW in the ratio πœ†βˆΆ1βˆ’πœ†. This can be done in two ways.
In the first case, if we wish to obtain the optimal value of πœ† at any cfl, then we compute ξ€œπ‘Ÿ10ξ€œπ‘€10eeldldπ‘‘π‘€π‘‘π‘Ÿ,(3.9) which will be in terms of πœ†.
The value of π‘Ÿ1 is chosen to suit the region of stability of the numerical scheme under consideration while 𝑀1 is chosen such that the approximated RPEβ‰₯0 for π‘Ÿβˆˆ[0,π‘Ÿ1]. In cases where phase wrapping phenomenon does not occur, we use the exact RPE instead of the approximated RPE and in that case, π‘€βˆˆ[0,πœ‹].
The second way to optimise a scheme made up of a linear combination of Beam-Warming and Lax-Wendroff is to compute the IEELDLD as βˆ«π‘€10eeldld𝑑𝑀 and the integral obtained in that case will be a function of π‘Ÿ and πœ†. Then a 3-D plot of this integral with respect to π‘Ÿβˆˆ[0,π‘Ÿ1] and πœ†βˆˆ[0,1] enables the respective optimal values of π‘Ÿ and πœ† to be located. The optimised scheme obtained will be defined in terms of both a cfl number and the optimal value of πœ† to be used.

Considerable and extensive work on the technique of Minimised Integrated Exponential Error for Low Dispersion and Low Dissipation has been carried out in Appadu and Dauhoo [18], Appadu and Dauhoo [19], Appadu [20–22].

In Appadu and Dauhoo [18], we have obtained the optimal cfl for some explicit methods like Lax-Wendroff, Beam-Warming, Crowley, Upwind Leap-Frog, and Fromm’s schemes. In Appadu and Dauhoo [19], we have combined some spatial discretisation schemes with the optimised time discretisation method proposed by Tam and Webb [3] in order to approximate the linear 1-D advection equation. These spatial derivatives are a standard 7-point and 6th-order central difference scheme (ST7), a standard 9-point and 8th-order central difference scheme (ST9) and optimised spatial schemes designed by Tam and Webb [3], Lockard et al. [13], Zingg et al. [14], Zhuang and Chen [15] and Bogey and Bailly [16]. The results from some numerical experiments were quantified into dispersion and dissipation errors and we have found that the quality of the results is dependent on the choice of the cfl number even for optimised methods, though to a much lesser degree as compared to standard methods.

Moreover, in Appadu [20], we obtain the optimal cfl of some multilevel schemes in 1-D. These schemes are of high order in space and time and have been designed by Wang and Liu [23]. We have also optimised the parameters in the family of third-order schemes proposed by Takacs [24]. The optimal cfl of the 2-D CFLF4 scheme which is a composite method made up of Corrected Lax-Friedrichs and the two-step Lax-Friedrichs developed by Liska and Wendroff [25] has been computed and some numerical experiments have been performed such as 2-D solid body rotation test [26], 2-D acoustics [27], and 2-D circular Riemann problem [26]. We have shown that better results are obtained when the optimal parameters obtained using MIEELDLD are used.

Some more interesting features of MIEELDLD are detailed in Appadu [21]. In that paper, we extend the measures used by Tam and Webb [3], Bogey and Bailly [16], Berland et al. [28] in a computational aeroacoustics framework to suit them in a computational fluid dynamics framework such that the optimal cfl of some known numerical methods can be obtained. Thus, we define the following integrals: integrated error from Tam and Webb [3], (IETAM), integrated error from Bogey and Bailly [16] ((IEBOGEY), and integrated error from Berland et al. [28] (IEBERLAND) as follows:ξ€œIETAM=𝑀10||||1βˆ’RPE2ξ€œπ‘‘π‘€,IEBOGEY=𝑀10||||ξ€œ1βˆ’RPE𝑑𝑀,IEBERLAND=𝑀10||||(1βˆ’π›Ό)1βˆ’RPE+𝛼(1βˆ’AFM)𝑑𝑀.(3.10)

The optimal cfl is obtained by plotting the respective integral with respect to the cfl number and locating the cfl at which the integral is least. The techniques used to obtain the quantities IETAM, IEBOGEY, and IEBERLAND are named Minimised Integrated Error from Tam and Webb [3] (MIETAM), Minimised Integrated Error from Bogey and Bailly [16] (MIEBOGEY), and Minimised Integrated Error from Berland et al. [28] (MIEBERLAND) respectively. It is shown that MIEELDLD has an upper hand over the other techniques of optimisation: MIETAM, MIEBOGEY, and MIEBERLAND.

The work in Appadu [22] helps us to understand why not all composite schemes can be effective to capture shocks with minimum dispersion and dissipation. The findings concluded are that some efficient composite methods to approximate the 1-D linear advection equation are as follows: composite methods using Lax-Wendroff and Beam-Warming as either predictor or corrector steps, a linear combination of either Lax-Wendroff and Beam-Warming schemes or MacCormack and two-step Lax-Friedrichs and the composite MacCormack/Lax-Friedrichs schemes [29].

4. Modification to Space Discretisation Scheme Proposed by Tam and Webb [3]

Tam and Webb [3] constructed a 7-pt and 4th-order central difference method based on a minimization of the dispersion error.

They approximated πœ•π‘’/πœ•π‘₯ at π‘₯0 asπœ•π‘’=1πœ•π‘₯β„Ž3𝑖=βˆ’3π‘Žπ‘–π‘’ξ€·π‘₯0ξ€Έ+π‘–β„Ž,(4.1) where β„Ž is the spacing of a uniform mesh and the coefficients π‘Žπ‘– are such that π‘Žπ‘–=βˆ’π‘Žβˆ’π‘–, providing a scheme without dissipation. On applying spatial Fourier Transform to (4.1), the effective wavenumber πœƒβˆ—β„Ž of the scheme is obtained and it is given byπœƒβˆ—β„Ž=23𝑖=1π‘Žπ‘–sin(π‘–πœƒβ„Ž).(4.2)

Taylor expansion of πœƒβˆ—β„Ž about πœƒβ„Ž from (4.2) gives the following:2π‘Ž1ξ‚€1πœƒβ„Žβˆ’6(πœƒβ„Ž)3+1120(πœƒβ„Ž)5+2π‘Ž2ξ‚€12πœƒβ„Žβˆ’6(2πœƒβ„Ž)3+1120(2πœƒβ„Ž)5+2π‘Ž3ξ‚€13πœƒβ„Žβˆ’6(3πœƒβ„Ž)3+1120(3πœƒβ„Ž)5+β‹―.(4.3)

To obtain a 4th-order accurate method, we must have2π‘Ž1+4π‘Ž2+6π‘Ž3π‘Ž=1,1+8π‘Ž2+27π‘Ž3=0.(4.4)

Since, we have 2 equations and 3 unknowns, we can choose, for instance, say π‘Ž1 as a free parameter. Thus,π‘Ž2=9βˆ’4205π‘Ž1,π‘Ž3=15ξ‚€π‘Ž1βˆ’23.(4.5)

Hence, the numerical wavenumber can be expressed asπœƒβˆ—β„Žβ‰ˆ2π‘Ž1ξ‚€9sin(πœƒβ„Ž)+2βˆ’4205π‘Ž11sin(2πœƒβ„Ž)+25π‘Ž1βˆ’215sin(3πœƒβ„Ž).(4.6)

The optimisation procedure used by Tam and Webb [3] is to find π‘Ž1 which minimizes the integrated error, 𝐸 defined asξ€œπΈ=01.1||πœƒβˆ—||β„Žβˆ’πœƒβ„Ž2𝑑(πœƒβ„Ž).(4.7) The value obtained for π‘Ž1 is 0.7708823806. The corresponding values for π‘Ž2 and π‘Ž3 are βˆ’0.1667059045 and 0.0208431428.

Hence, the scheme developed by Tam and Webb [3] has numerical wavenumber, πœƒβˆ—β„Ž, and group velocity given byπœƒβˆ—β„Ž=1.5417647612sin(πœƒβ„Ž)βˆ’0.3334118090sin(2πœƒβ„Ž)+0.0416862856sin(3πœƒβ„Ž),(4.8)groupvelocity=1.5417647612cos(πœƒβ„Ž)βˆ’0.6668236180cos(2πœƒβ„Ž)+0.1250588568cos(3πœƒβ„Ž)(4.9) and is termed as β€œTAM” method.

We next consider the numerical wavenumber in (4.2) and use the technique of MIEELDLD to find optimal values of π‘Ž1, π‘Ž2, and π‘Ž3. The integrated error using MIEELDLD is given byξ€œπΈ=01.1ξ€·ξ€·||β„œξ€·πœƒexpβˆ—β„Žξ€Έ||+||β„‘ξ€·πœƒβˆ’πœƒβ„Žβˆ—β„Žξ€Έ||ξ€Έ||||β„œξ€·πœƒ+expβˆ—β„Žξ€Έ||βˆ’||β„‘ξ€·πœƒβˆ’πœƒβ„Žβˆ—β„Žξ€Έ||||ξ€Έβˆ’2𝑑(πœƒβ„Ž).(4.10) Since we are considering a 7-point and 4th-order central difference method, the numerical wavenumber, πœƒβˆ—β„Ž, does not have an imaginary part, that is, β„‘(πœƒβˆ—β„Ž)=0. Hence, (4.10) simplifies toξ€œπΈ=01.1ξ€·||β„œξ€·πœƒ2expβˆ—β„Žξ€Έ||ξ€Έβˆ’πœƒβ„Žβˆ’2𝑑(πœƒβ„Ž),(4.11) and on minimising this integral using the function NLPSolve in maple, we obtain π‘Ž1 as 0.7677206709. Corresponding values for π‘Ž2 and π‘Ž3 are 0.1641765367 and 0.0202108009, respectively.

Hence we have obtained a modified method which is 7-point and of 4th-order which we term as β€œTAM-NEW” method. Expressions for the numerical wavenumber and the group velocity of the β€œTAM-NEW” method are given byπœƒβˆ—β„Ž=1.5354413418sin(πœƒβ„Ž)βˆ’0.3283530734sin(2πœƒβ„Ž)+0.0404216018sin(3πœƒβ„Ž),(4.12)groupvelocity=1.5354413418cos(πœƒβ„Ž)βˆ’0.6567061468cos(2πœƒβ„Ž)+0.1212648054cos(3πœƒβ„Ž).(4.13) We next perform a spectral analysis of the two methods. We compare the variation of numerical wavenumber versus the exact wavenumber in Figure 1. A plot of the dispersion error versus the exact wavenumber is depicted in Figure 2. The dispersion error for TAM-NEW is slightly less than that for TAM for 0<πœƒβ„Žβ‰€1, but for 1β‰€πœƒβ„Žβ‰€πœ‹/2, the dispersion error from TAM is slightly less than that for TAM-NEW.

783101.fig.001
Figure 1: Plot of the variation of numerical wavenumber versus exact wavenumber for the methods: TAM, TAM-NEW, ZINGG, and ZINGG-NEW.
783101.fig.002
Figure 2: Plot of the dispersion error on a logarithmic scale versus exact wavenumber for the methods: TAM, TAM-NEW, ZINGG, and ZINGG-NEW.

We now compare quantitatively these two methods: TAM and TAM-NEW. We use four accuracy limits [5, 16] defined as follows:||πœƒβˆ—||β„Žβˆ’πœƒβ„Žπœ‹β‰€5Γ—10βˆ’3,||πœƒβˆ—||β„Žβˆ’πœƒβ„Žπœ‹β‰€5Γ—10βˆ’4,||||π‘‘ξ€·πœƒπ‘‘(πœƒβ„Ž)βˆ—β„Žξ€Έ||||βˆ’1.0≀5Γ—10βˆ’3,||||π‘‘ξ€·πœƒπ‘‘(πœƒβ„Ž)βˆ—β„Žξ€Έ||||βˆ’1.0≀5Γ—10βˆ’4(4.14) and compute the minimum number of points per wavelength needed to resolve a wave for each of the four accuracy limits. The results are summarised in Table 1.

tab1
Table 1: Comparing the dispersive and group velocity properties for two spatial discretisation methods β€œTAM” and β€œTAM-NEW” in terms of number of points per wavelength (to 4 d.p).

It is seen that the scheme β€œTAM-NEW” is not superior to the TAM method as for a given accuracy it requires more points per wavelength in regard to the dispersive and group velocity properties. This is because the technique of MIEELDLD aims to minimize both dispersion and dissipation in numerical methods but here our aim is to construct a 7-point and 4th-order central difference method with no dissipation.

5. Modification to Space Discretisation Scheme Developed by Lockard et al. [13]

Lockard et al. [13] constructed a 7-point and 4th-order difference method by approximating πœ•π‘’/πœ•π‘₯ at π‘₯0 asπœ•π‘’=1πœ•π‘₯β„Ž3𝑖=βˆ’4π‘Žπ‘–π‘’ξ€·π‘₯0ξ€Έ+π‘–β„Ž,(5.1) and therefore the real and imaginary parts of the numerical wavenumber are obtained as follows:β„œξ€·πœƒβˆ—β„Žξ€Έ=βˆ’π‘Žβˆ’4sin(4πœƒβ„Ž)βˆ’π‘Žβˆ’3sin(3πœƒβ„Ž)βˆ’π‘Žβˆ’2sin(2πœƒβ„Ž)βˆ’π‘Žβˆ’1sin(πœƒβ„Ž)+π‘Ž1sin(πœƒβ„Ž)+π‘Ž2sin(2πœƒβ„Ž)+π‘Ž3β„‘ξ€·πœƒsin(3πœƒβ„Ž),(5.2)βˆ—β„Žξ€Έξ€·π‘Ž=βˆ’βˆ’4cos(4πœƒβ„Ž)+π‘Žβˆ’3cos(3πœƒβ„Ž)+π‘Žβˆ’2cos(2πœƒβ„Ž)+π‘Žβˆ’1cos(πœƒβ„Ž)+π‘Ž0+π‘Ž1cos(πœƒβ„Ž)+π‘Ž2cos(2πœƒβ„Ž)+π‘Ž3ξ€Έ.cos(3πœƒβ„Ž)(5.3)

To obtain a 4th-order method, we require 4 conditions based on the real and imaginary parts of πœƒβˆ—β„Ž, namely, π‘Ž1+2π‘Ž2+3π‘Ž3βˆ’4π‘Žβˆ’4βˆ’3π‘Žβˆ’3βˆ’2π‘Žβˆ’2βˆ’π‘Žβˆ’1=1,βˆ’π‘Ž1βˆ’8π‘Ž2βˆ’27π‘Ž3+64π‘Žβˆ’4+27π‘Žβˆ’3+8π‘Žβˆ’2+π‘Žβˆ’1π‘Ž=0,0+π‘Ž1+π‘Ž2+π‘Ž3+π‘Žβˆ’4+π‘Žβˆ’3+π‘Žβˆ’2+π‘Žβˆ’1=0,βˆ’π‘Ž1βˆ’4π‘Ž2βˆ’9π‘Ž3βˆ’16π‘Žβˆ’4βˆ’9π‘Žβˆ’3βˆ’4π‘Žβˆ’2βˆ’π‘Žβˆ’1=0.(5.4)

The coefficients obtained by Lockard et al. [13] areπ‘Žβˆ’4=0.0103930209,π‘Žβˆ’3=βˆ’0.0846974943,π‘Žβˆ’2π‘Ž=0.3420311831,βˆ’1=βˆ’1.0526812838,π‘Ž0=0.2872741244,π‘Ž1π‘Ž=0.5861624738,2=βˆ’0.0981442817,π‘Ž3=0.0096622576.(5.5)

Hence, the real and imaginary parts of πœƒβˆ—β„Ž for LOCKARD scheme are given as follows:β„œξ€·πœƒβˆ—β„Žξ€Έ=1.63884375sin(πœƒβ„Ž)βˆ’0.44017538sin(2πœƒβ„Ž)+0.09433201sin(3πœƒβ„Ž)βˆ’0.01039020sin(4πœƒβ„Ž),(5.6)β„‘(πœƒβˆ—β„Ž)=βˆ’0.28727412+0.46651881cos(πœƒβ„Ž)βˆ’0.24388682cos(2πœƒβ„Ž)+0.07500749cos(3πœƒβ„Ž)βˆ’0.01039020cos(4πœƒβ„Ž), (23)respectively.

We now obtain a modification to the scheme developed by Lockard et al. [13]. We consider the numerical wavenumber in (5.2) and (5.3) and replace π‘Žβˆ’1, π‘Ž0, π‘Ž1, π‘Ž2, and π‘Ž3 in terms of π‘Žβˆ’2, π‘Žβˆ’3, π‘Žβˆ’4, and πœƒβ„Ž. Our aim is to minimise the following integral:ξ€œπΈ=01.1ξ€·ξ€·||β„œξ€·πœƒexpβˆ—β„Žξ€Έ||+||β„‘ξ€·πœƒβˆ’πœƒβ„Žβˆ—β„Žξ€Έ||ξ€Έ||||β„œξ€·πœƒ+expβˆ—β„Žξ€Έ||βˆ’||β„‘ξ€·πœƒβˆ’πœƒβ„Žβˆ—β„Žξ€Έ||||ξ€Έβˆ’2𝑑(πœƒβ„Ž).(5.8)

The integral is a function of π‘Žβˆ’2, π‘Žβˆ’3, and π‘Žβˆ’4. We use the function NLPSolve and obtain optimal values for π‘Žβˆ’4, π‘Žβˆ’3, and π‘Žβˆ’2 as 0.0113460667, βˆ’0.0891980000, and 0.3499980000. Then the values of the other unknowns can be obtained and we are out withπ‘Žβˆ’1=βˆ’1.0582666667,π‘Ž0=0.2866010000,π‘Ž1π‘Ž=0.5895196001,2=βˆ’0.1,π‘Ž3=0.01.(5.9)

The modified method is termed as β€œLOCKARD-NEW” and has real and imaginary parts of its numerical wavenumber described byβ„œξ€·πœƒβˆ—β„Žξ€Έβ„‘ξ€·πœƒ=1.6477862670sin(πœƒβ„Ž)βˆ’0.4499980000sin(2πœƒβ„Ž)+0.0991980000sin(3πœƒβ„Ž)βˆ’0.0113460667sin(4πœƒβ„Ž),(5.10)βˆ—β„Žξ€Έ=βˆ’0.2866010000+0.4687470669cos(πœƒβ„Ž)βˆ’0.2499980000cos(2πœƒβ„Ž)+0.0791980000cos(3πœƒβ„Ž)βˆ’0.0113460667cos(4πœƒβ„Ž),(5.11) respectively.

We next perform a spectral analysis of the two methods: LOCKARD and LOCKARD-NEW. We compare the variation of numerical wavenumber versus the exact wavenumber in Figure 3 and in Figure 4, we have the plot of the dispersion error versus the exact wavenumber.

783101.fig.003
Figure 3: Plot of the variation of numerical wavenumber versus exact wavenumber for the methods LOCKARD and LOCKARD-NEW.
783101.fig.004
Figure 4: Plot of the variation of dispersion error in logarithmic scale versus exact wavenumber for LOCKARD and LOCKARD-NEW schemes.

We now compare quantitatively the two methods by computing the minimum number of points per wavelength needed to resolve a wave for each of the four accuracy limits and the results are summarized in Table 2.

tab2
Table 2: Comparing the dispersive and group velocity properties for two spatial methods LOCKARD and LOCKARD-NEW in terms of number of points per wavelength (to 4 d.p).

Clearly, LOCKARD-NEW has appreciably better phase and group velocity properties as compared to LOCKARD scheme.

6. Modification to Spatial Discretisation Scheme Developed by Zingg et al. [14]

Zingg et al. [14] constructed a 4-point and 4th-order difference method. They approximated πœ•π‘’/πœ•π‘₯ at π‘₯0 byπœ•π‘’=1πœ•π‘₯β„Ž3𝑖=1π‘Žπ‘–ξ€·π‘’ξ€·π‘₯0ξ€Έξ€·π‘₯+π‘–β„Žβˆ’π‘’0+1βˆ’π‘–β„Žξ€Έξ€Έβ„Žξƒ©π‘‘0𝑒π‘₯0ξ€Έ+3𝑖=1𝑑𝑖𝑒π‘₯0ξ€Έξ€·π‘₯+π‘–β„Ž+𝑒0ξƒͺ.βˆ’π‘–β„Žξ€Έξ€Έ(6.1)

The real and imaginary parts of the numerical wavenumber are obtained asβ„œξ€·πœƒβˆ—β„Žξ€Έξ€·π‘Ž=21sin(πœƒβ„Ž)+π‘Ž2sin(2πœƒβ„Ž)+π‘Ž3ξ€Έβ„‘ξ€·πœƒsin(3πœƒβ„Ž),(6.2)βˆ—β„Žξ€Έξ€·π‘‘=βˆ’0+2𝑑1cos(πœƒβ„Ž)+2𝑑2cos(2πœƒβ„Ž)+2𝑑3ξ€Έcos(3πœƒβ„Ž).(6.3)

The conditions to have a 4th-order difference method are as follows.

(i)If we consider β„œ(πœƒβˆ—β„Ž), then 2π‘Ž1+4π‘Ž2+6π‘Ž3βˆ’1=1,6π‘Ž1βˆ’86π‘Ž2βˆ’276π‘Ž3=0,(6.4) and these two conditions give π‘Ž2=9βˆ’4205π‘Ž1,π‘Ž(6.5)3=15ξ‚€π‘Ž1βˆ’23.(6.6)(ii)If we consider β„‘(πœƒβˆ—β„Ž), then 𝑑0+2𝑑1+2𝑑2+2𝑑3=0,βˆ’π‘‘1βˆ’4𝑑2βˆ’9𝑑3=0,(6.7) and this gives 𝑑0=6𝑑2+16𝑑3,𝑑(6.8)1=βˆ’4𝑑2βˆ’9𝑑3.(6.9)

Based on the optimisation performed by Zingg et al. [14], the following values are obtained:π‘Ž1=0.75996126,π‘Ž2=βˆ’0.15812197,π‘Ž3=0.01876090,𝑑0𝑑=0.1,1=βˆ’0.07638462,𝑑2=0.03228962,𝑑3=βˆ’0.00590500.(6.10)

We now obtain a modification to the scheme proposed by Zingg et al. [14] using MIEELDLD. We considerβ„‘ξ€·πœƒβˆ—β„Žξ€Έξ€·π‘‘=βˆ’0+2𝑑1cos(πœƒβ„Ž)+2𝑑2cos(2πœƒβ„Ž)+2𝑑3ξ€Έcos(3πœƒβ„Ž).(6.11)

Since Im(πœƒβˆ—β„Ž) must be negative and the method must have sufficient dissipation, we can choose 𝑑0=0.1 and hence obtain𝑑15=βˆ’8𝑑2π‘‘βˆ’0.05625,33=0.00625βˆ’8𝑑2.(6.12)

We next plot Im(πœƒβˆ—β„Ž) versus 𝑑2 versus πœƒβ„Žβˆˆ[0,2πœ‹] and obtain the range of 𝑑2 such that Im(πœƒβˆ—β„Ž)<0. The maximum value of 𝑑2 is 0.0323. Having fixed the values of 𝑑0 as 0.1 and 𝑑2 as 0.0323, now we can compute the values of 𝑑1 and 𝑑3. We are out with 𝑑1=βˆ’0.0764375 and 𝑑3=βˆ’5.8625Γ—10βˆ’3. Hence, we minimize the following integral:ξ€œ01.1ξ€·||β„œξ€·πœƒexpβˆ—β„Žξ€Έ||+||β„‘ξ€·πœƒβˆ’πœƒβ„Žβˆ—β„Žξ€Έ||ξ€Έ||||β„œξ€·πœƒ+expβˆ—β„Žξ€Έ||βˆ’||β„‘ξ€·πœƒβˆ’πœƒβ„Žβˆ—β„Žξ€Έ||||βˆ’2.0𝑑(πœƒβ„Ž)(6.13) which is a function of π‘Ž1, using NLPSolve.

We obtain π‘Ž1=0.7643155206, and therefore, using (6.5) and (6.6), we obtain π‘Ž2=βˆ’0.1614524165 and π‘Ž3=0.0195297708.

Hence, the real and imaginary parts of the real and imaginary parts of the numerical wavenumber of the scheme ZINGG-NEW are as follows:β„œξ€·πœƒβˆ—β„Žξ€Έβ„‘ξ€·πœƒ=1.5286310410sin(πœƒβ„Ž)βˆ’0.3229048330sin(2πœƒβ„Ž)+0.0390595416sin(3πœƒβ„Ž),(6.14)βˆ—β„Žξ€Έ=βˆ’0.1+0.1528750000cos(πœƒβ„Ž)βˆ’0.0646000000cos(2πœƒβ„Ž)+0.0117250000cos(3πœƒβ„Ž).(6.15)

Plots of β„œ(πœƒβˆ—β„Ž) versus πœƒβ„Ž and also for β„‘(πœƒβˆ—β„Ž) versus πœƒβ„Ž for ZINGG and ZINGG-NEW schemes are depicted in Figures 1 and 6, respectively. It is observed based on Figure 6 that the two methods have almost the same dissipation error for πœƒβ„Žβˆˆ[0,πœ‹]. Based on (Figure 1), we observe that for πœƒβ„Ž<0.2 and 0.8<πœƒβ„Ž<πœ‹/2, the dispersion error from ZINGG-NEW is less than that for ZINGG. For 0.2<πœƒβ„Ž<0.8, the dispersion error from ZINGG is less than ZINGG-NEW.

Based on Table 3, for the four accuracy limits tested, we can conclude that the new scheme developed is superior to the ZINGG method in terms of both dispersive and group velocity properties as it requires less points per wavelength in all the four cases.

tab3
Table 3: Comparing the dispersive and group velocity properties for two spatial methods ZINGG and ZINGG-NEW in terms of number of points per wavelength (to 4 d.p).

7. Modification to Spatial Scheme Developed by Zhuang and Chen [15]

Zhuang and Chen [15] constructed a 7-point and 4th-order difference method by approximating πœ•π‘’/πœ•π‘₯ at π‘₯0 asπœ•π‘’=1πœ•π‘₯β„Ž2𝑖=βˆ’4π‘Žπ‘–π‘’ξ€·π‘₯0ξ€Έ+π‘–β„Ž,(7.1) and therefore the real and imaginary parts of the numerical wavenumber are obtained asβ„œξ€·πœƒβˆ—β„Žξ€Έ=π‘Ž1sin(πœƒβ„Ž)+π‘Ž2sin(2πœƒβ„Ž)βˆ’π‘Žβˆ’4sin(4πœƒβ„Ž)βˆ’π‘Žβˆ’3sin(3πœƒβ„Ž)βˆ’π‘Žβˆ’2sin(2πœƒβ„Ž)βˆ’π‘Žβˆ’1β„‘ξ€·πœƒsin(πœƒβ„Ž),(7.2)βˆ—β„Žξ€Έξ€·π‘Ž=βˆ’βˆ’4cos(4πœƒβ„Ž)+π‘Žβˆ’3cos(3πœƒβ„Ž)+π‘Žβˆ’2cos(2πœƒβ„Ž)+π‘Žβˆ’1cos(πœƒβ„Ž)+π‘Ž0+π‘Ž1cos(πœƒβ„Ž)+π‘Ž2ξ€Έ.cos(2πœƒβ„Ž)(7.3)

To obtain a 4th-order method, we require 4 conditions based on the real and imaginary parts of πœƒβˆ—β„Ž:π‘Ž1+2π‘Ž2βˆ’4π‘Žβˆ’4βˆ’3π‘Žβˆ’3βˆ’2π‘Žβˆ’2βˆ’π‘Žβˆ’1=1,βˆ’π‘Ž1βˆ’8π‘Ž2+64π‘Žβˆ’4+27π‘Žβˆ’3+8π‘Žβˆ’2+π‘Žβˆ’1π‘Ž=0,0+π‘Ž1+π‘Ž2+π‘Žβˆ’4+π‘Žβˆ’3+π‘Žβˆ’2+π‘Žβˆ’1=0,βˆ’π‘Ž1βˆ’16π‘Žβˆ’4βˆ’4π‘Žβˆ’2βˆ’π‘Ž1βˆ’4π‘Ž2βˆ’9π‘Žβˆ’3=0.(7.4)

These simplify to the following if we let π‘Žβˆ’4,π‘Žβˆ’3,π‘Žβˆ’2 as free parameters:π‘Ž2=10π‘Žβˆ’4+4π‘Žβˆ’3+π‘Žβˆ’2βˆ’16,π‘Žβˆ’1=βˆ’20π‘Žβˆ’4βˆ’10π‘Žβˆ’3βˆ’4π‘Žβˆ’2βˆ’13,π‘Ž0=45π‘Žβˆ’4+20π‘Žβˆ’3+6π‘Žβˆ’2βˆ’12,π‘Ž1=βˆ’36π‘Žβˆ’4βˆ’15π‘Žβˆ’3βˆ’4π‘Žβˆ’2+1.(7.5)

On plugging π‘Ž2, π‘Žβˆ’1, π‘Ž0, and π‘Ž1 in terms of functions of π‘Žβˆ’4, π‘Žβˆ’3, π‘Žβˆ’2 in (7.2) and (7.3), we getβ„œξ€·πœƒβˆ—β„Žξ€Έξ‚€=βˆ’5165π‘Žβˆ’4+π‘Žβˆ’3βˆ’4115sin(πœƒβ„Ž)+6ξ€·60π‘Žβˆ’4+24π‘Žβˆ’3ξ€Έβˆ’1sin(2πœƒβ„Ž)βˆ’π‘Žβˆ’3sin(3πœƒβ„Ž)βˆ’π‘Žβˆ’4β„‘ξ€·πœƒsin(4πœƒβ„Ž),(7.6)βˆ—β„Žξ€Έ=12βˆ’45π‘Žβˆ’4βˆ’20π‘Žβˆ’3βˆ’6π‘Žβˆ’2+16ξ€·cos(πœƒβ„Ž)336π‘Žβˆ’4+150π‘Žβˆ’3+48π‘Žβˆ’2ξ€Έβˆ’4βˆ’π‘Žβˆ’3cos(3πœƒβ„Ž)βˆ’π‘Žβˆ’41cos(4πœƒβ„Ž)+6ξ€·βˆ’60π‘Žβˆ’4βˆ’24π‘Žβˆ’3βˆ’12π‘Žβˆ’2ξ€Έ+1cos(2πœƒβ„Ž).(7.7)

The coefficients obtained by Zhuang and Chen [15] areπ‘Žβˆ’4=0.0161404967,π‘Žβˆ’3=βˆ’0.1228212790,π‘Žβˆ’2π‘Ž=0.4553322778,βˆ’1=βˆ’1.2492595883,π‘Ž0=0.5018904380,π‘Ž1π‘Ž=0.4399321927,2=βˆ’0.0412145379,(7.8) and, therefore, the real and imaginary parts of πœƒβˆ—β„Ž are given as follows:β„œξ€·πœƒβˆ—β„Žξ€Έβ„‘ξ€·πœƒ=1.689191781sin(πœƒβ„Ž)βˆ’0.4965468157sin(2πœƒβ„Ž)+0.1228212790sin(3πœƒβ„Ž)βˆ’0.0161404967sin(4πœƒβ„Ž),(7.9)βˆ—β„Žξ€Έ=βˆ’0.5018904390+0.8093273950cos(πœƒβ„Ž)βˆ’0.4141177399cos(2πœƒβ„Ž)+0.1228212790cos(3πœƒβ„Ž)βˆ’0.0161404967cos(4πœƒβ„Ž),(7.10) respectively.

We now obtain a modification to the scheme developed by Zhuang and Chen [15]. We consider the numerical wavenumber in (7.6) and (7.7) and minimise the following integralξ€œπΈ=01.1ξ€·ξ€·||β„œξ€·πœƒexpβˆ—β„Žξ€Έ||+||β„‘ξ€·πœƒβˆ’πœƒβ„Žβˆ—β„Žξ€Έ||ξ€Έ||||β„œξ€·πœƒ+expβˆ—β„Žξ€Έ||βˆ’||β„‘ξ€·πœƒβˆ’πœƒβ„Žβˆ—β„Žξ€Έ||||ξ€Έβˆ’2𝑑(πœƒβ„Ž).(7.11)

The integral is a function of π‘Žβˆ’4, π‘Žβˆ’3, and π‘Žβˆ’2. We use the function NLPSolve and obtain optimal values for π‘Žβˆ’4, π‘Žβˆ’3, and π‘Žβˆ’2 as 0.01575, βˆ’0.122, and 0.4553 respectively. Corresponding values for π‘Ž2, π‘Žβˆ’1, π‘Ž0, and π‘Ž1 are then obtained as βˆ’0.0418666600, βˆ’1.2495333300, 0.5005500000, and 0.4418000000, respectively.

The modified method is termed as ZHUANG-NEW and has real and imaginary parts of its numerical wavenumber described byβ„œξ€·πœƒβˆ—β„Žξ€Έβ„‘ξ€·πœƒ=1.6913333333sin(πœƒβ„Ž)βˆ’0.4971666667sin(2πœƒβ„Ž)+0.1220000000sin(3πœƒβ„Ž)βˆ’0.0157500000sin(4πœƒβ„Ž),(7.12)βˆ—β„Žξ€Έ=βˆ’0.5005500000+0.8077333330cos(πœƒβ„Ž)βˆ’0.4134333333cos(2πœƒβ„Ž)+0.1220000000cos(3πœƒβ„Ž)βˆ’0.0157500000cos(4πœƒβ„Ž),(7.13) respectively.

We next perform a spectral analysis of the two methods: ZHUANG and ZHUANG-NEW. We compare the variation of real part and imaginary parts of the numerical wavenumber versus the exact wavenumber in Figures 7 and 5, respectively. We have the plot of the dispersion error versus the exact wavenumber in Figure 8 and we observe that, for 0<πœƒβ„Ž<1, ZHUANG-NEW is slightly better than ZHUANG in terms of dispersive properties.

783101.fig.005
Figure 5: Plot of the imaginary part of numerical wavenumber versus exact wavenumber for the methods: LOCKARD, LOCKARD-NEW, ZHUANG, and ZHUANG-NEW.
783101.fig.006
Figure 6: Plot of the imaginary part of numerical wavenumber versus exact wavenumber for ZINGG and ZINGG-NEW.
783101.fig.007
Figure 7: Plot of the variation of numerical wavenumber versus exact wavenumber for the methods: ZHUANG, ZHUANG-NEW, BOGEY, and BOGEY-NEW.
783101.fig.008
Figure 8: Plot of the variation of dispersion error in logarithmic scale versus exact wavenumber for the methods ZHUANG and ZHUANG-NEW.

We now compare quantitatively these two methods. We compute the minimum number of points per wavelength needed to resolve a wave for each of the four accuracy limits. The results are summarized in Table 4.

tab4
Table 4: Comparing the dispersive and group velocity properties for two spatial methods ZHUANG and ZHUANG-NEW in terms of number of points per wavelength (to 4 d.p).

ZHUANG-NEW requires fewer points per wavelength than ZHUANG scheme for |(πœƒβˆ—β„Žβˆ’πœƒβ„Ž)/πœ‹|≀0.005.

8. Modification to Spatial Discretisation Scheme Developed by Bogey and Bailly [16]

Bogey and Bailly [16] modified the measure used by Tam and Webb [3] by minimizing the relative difference between πœƒβ„Ž and πœƒβˆ—β„Ž. They define the integrated error, 𝐸, asξ€œπΈ=πœ‹/2πœ‹/16||πœƒβˆ—||β„Žβˆ’πœƒβ„Žπœƒβ„Žπ‘‘(πœƒβ„Ž).(8.1)

Bogey and Bailly [3] use a 9-point stencil with coefficients π‘Žβˆ’4, π‘Žβˆ’3, π‘Žβˆ’2, π‘Žβˆ’1, π‘Ž0, π‘Ž1, π‘Ž2, π‘Ž3, π‘Ž4 and choose π‘Ž0=0, π‘Žβˆ’1=βˆ’π‘Ž1, π‘Žβˆ’2=βˆ’π‘Ž2, π‘Žβˆ’3=βˆ’π‘Ž3, and π‘Žβˆ’4=βˆ’π‘Ž4 and therefore the numerical wavenumber can be written asπœƒβˆ—ξ€·π‘Žβ„Ž=21sin(πœƒβ„Ž)+π‘Ž2sin(2πœƒβ„Ž)+π‘Ž3sin(3πœƒβ„Ž)+π‘Ž4ξ€Έsin(4πœƒβ„Ž).(8.2)

To obtain a 4th-order method, π‘Ž1 and π‘Ž2 are chosen such asπ‘Ž1=23+5π‘Ž3+16π‘Ž4,π‘Ž21=βˆ’6ξ‚€12+24π‘Ž3+60π‘Ž4,(8.3) respectively.

The coefficients π‘Ž3 and π‘Ž4 are chosen to minimize the integrated error defined in (8.1), and the values which Bogey and Bailly [16] have obtained are as follows:π‘Ž1=0.841570125,π‘Ž2=βˆ’0.2446786318,π‘Ž3=0.0594635848,π‘Ž4=βˆ’0.0076509040.(8.4)

We now construct a method based on a 9-point stencil using MIEELDLD. The wavenumber is set as follows:πœƒβˆ—ξ‚€2β„Ž=23+5π‘Ž3+16π‘Ž4ξ‚ξ‚€βˆ’1sin(πœƒβ„Ž)+212βˆ’4π‘Ž3βˆ’10π‘Ž4sin(2πœƒβ„Ž)+2π‘Ž3sin(3πœƒβ„Ž).(8.5)

The integrated error using MIEELDLD is defined asξ€œπœ‹/2πœ‹/16ξ€·||β„œξ€·πœƒ2expβˆ—β„Žξ€Έ||ξ€Έβˆ’πœƒβ„Žβˆ’2𝑑(πœƒβ„Ž),(8.6) which is a function of π‘Ž3 and π‘Ž4. Using NLPSolve, we obtain the optimal values of π‘Ž3 and π‘Ž4 as 0.0613000000 and βˆ’0.0080500000, respectively. Hence, we obtain π‘Ž1 and π‘Ž2 as 0.8443666667 and βˆ’0.2480333333, respectively.

Using MIEELDLD, a new scheme is obtained and is termed as BOGEY-NEW with its numerical wavenumber given byπœƒβˆ—β„Ž=1.6887333332sin(πœƒβ„Ž)βˆ’0.4960666667sin(2πœƒβ„Ž)+0.1226000000sin(3πœƒβ„Ž)0.0161000000sin(4πœƒβ„Ž).(8.7)

We next perform a spectral analysis of the two methods: BOGEY and BOGEY-NEW. We compare the variation of numerical wavenumber versus the exact wavenumber in Figures 7 and 9; we have the plot of the dispersion error versus the exact wavenumber.

783101.fig.009
Figure 9: Plot of the variation of dispersion error in logarithmic scale versus exact wavenumber.

We now compare quantitatively these two methods. We compute the minimum number of points per wavelength needed to resolve a wave for each of the four accuracy limits.

Table 5 indicates that BOGEY-NEW has appreciably better phase and group velocity properties as compared to BOGEY scheme.

tab5
Table 5: Comparing the dispersive and group velocity properties for two spatial methods BOGEY and BOGEY-NEW in terms of number of points per wavelength (to 4 d.p).

9. Optimized Time Discretisation Schemes

9.1. Time Discretisation Scheme by Tam et al. [17]

Tam et al. [17] have developed a time-marching scheme which is four-level and accurate up to π‘˜3. They expressedπ‘ˆ(𝑛+1)βˆ’π‘ˆ(𝑛)β‰ˆπ‘˜3𝑗=0π‘π‘—ξ‚€π‘‘π‘ˆξ‚π‘‘π‘‘(π‘›βˆ’π‘—).(9.1) We next summarize how the coefficients have been obtained.

The effective angular frequency of the time discretisation method is obtained asπœ›=𝐼(exp(βˆ’πΌπœ”π‘˜)βˆ’1)π‘˜βˆ‘3𝑗=0𝑏𝑗exp(πΌπ‘—πœ”π‘˜).(9.2)

For πœ›π‘˜ to approximate πœ”π‘˜ to order (πœ”π‘˜)4, we must have𝑏0+𝑏1+𝑏2+𝑏3𝑏=1,1+2𝑏2+3𝑏31=βˆ’2,𝑏1+4𝑏2+9𝑏3=13.(9.3)

Since we have 4 equations and 3 unknowns, we can choose 𝑏0 as a free parameter, and hence we have𝑏1=5312βˆ’3𝑏0,𝑏2=3𝑏0βˆ’163,𝑏3=2312βˆ’π‘0.(9.4)

Hence, we can express πœ› as follows:πœ›π‘˜=𝐴𝐢+𝐡𝐷1ξ€·+πΌπ΅πΆβˆ’π΄π·1𝐢2+(𝐷1)2,(9.5) where𝐴=sin(πœ”π‘˜),𝐡=cos(πœ”π‘˜)βˆ’1,𝐢=𝑏0+ξ‚€5312βˆ’3𝑏0cos(πœ”π‘˜)+3𝑏0βˆ’163cos(2πœ”π‘˜)+2312βˆ’π‘0𝐷cos(3πœ”π‘˜),1=ξ‚€5312βˆ’3𝑏0sin(πœ”π‘˜)+3𝑏0βˆ’163sin(2πœ”π‘˜)+2312βˆ’π‘0sin(3πœ”π‘˜).(9.6)

The weighted integral error incurred by using πœ› to approximate πœ”, used by Tam et al. [17], is computed as𝐸𝑇=ξ€œ0.5βˆ’0.5ξ€ΊπœŽ(β„œ(πœ›π‘˜)βˆ’πœ”π‘˜)2+(1βˆ’πœŽ)(β„‘(πœ›π‘˜))2𝑑(πœ”π‘˜),(9.7) and 𝜎 is chosen as 0.36.

On minimizing 𝐸𝑇, the value of 𝑏0 is obtained as 2.30255809 and therefore the corresponding values for 𝑏1, 𝑏2, and 𝑏3 are βˆ’2.49100760, 1.57434094, and βˆ’0.38589142, respectively.

9.2. Modified Temporal Discretisation Scheme Using MIEELDLD

We consider the equation in (9.5) which expresses πœ›π‘˜ in terms of πœ”π‘˜ and define the quantity, eeldld asξ€·||β„œ||+||β„‘||ξ€Έξ€·||β„œ||βˆ’||β„‘||ξ€Έexp(πœ›π‘˜)βˆ’πœ”π‘˜(πœ›π‘˜)+exp(πœ›π‘˜)βˆ’πœ”π‘˜(πœ›π‘˜)βˆ’2.(9.8)

We minimizeξ€œ0.5βˆ’0.5eeldld𝑑(πœ”π‘˜)(9.9) and this integral is a function of 𝑏0. Using NLPSolve, we obtain the value of 𝑏0 as 2.2796378228. A plot of 𝐸𝑇 versus 𝑏0 is shown in Figure 10.

783101.fig.0010
Figure 10: Plot of IEELDLD versus 𝑏0.

Corresponding values of 𝑏1,𝑏2,𝑏3 are obtained as βˆ’2.4222468020, 1.5055801360, and βˆ’0.3629711560. This modified temporal discretisation scheme obtained by modifying the temporal scheme of Tam et al. [17] is termed as β€œTAM-MODIFIED” scheme. Plots of β„œ(πœ›π‘˜) versus πœ”π‘˜ and β„‘(πœ›π‘˜) versus πœ”π‘˜ for the TAM-MODIFIED scheme are shown in Figures 11 and 12, respectively.

783101.fig.0011
Figure 11: Plot of β„œ(πœ›π‘˜) versus πœ”π‘˜.
783101.fig.0012
Figure 12: Plot of β„‘(πœ›π‘˜) versus πœ”π‘˜.

For |β„‘(πœ›π‘˜)|≀3Γ—10βˆ’3, we require πœ”π‘˜β‰€0.42.

9.3. Comparison between Temporal Discretisation Schemes: TAM and TAM-MODIFIED

Plots of β„œ(πœ›π‘˜) versus πœ”π‘˜ for the two methods are shown in Figure 13. We also compare their dispersive properties at two different levels of accuracy in terms of number of points per wavelength and the results are tabulated in Table 6. Clearly, TAM-MODIFIED is more superior as it requires less points per wavelength for the same accuracy.

tab6
Table 6: Comparing the dispersive properties for two temporal discretisation methods TAM and TAM-MODIFIED in terms of number of points per wavelength (to 4 d.p).
783101.fig.0013
Figure 13: Plot of β„œ(πœ›π‘˜) versus πœ”π‘˜ for TAM and TAM-MODIFIED schemes.

10. Stability of Some Multilevel Optimized Combined Spatial-Temporal Finite Difference Schemes

The stability of the combined spatial and temporal finite difference scheme developed by Tam and Webb [3] and Tam et al. [17], which is 7-point in space and 4-point in time and which is referred to as the Dispersion-Relation-Preserving (DRP) scheme, satisfies the stability condition, π‘Ÿβ‰€0.229 [3]. The condition on the spatial discretisation is that |(πœƒβˆ—β„Žβˆ’πœƒβ„Ž)/πœ‹|≀0.05 and this gives πœƒβ„Žβ‰€1.76. The interval 0<πœ›π‘˜β‰€0.4 has been chosen in order to maintain satisfactory temporal resolution and this interval is obtained by requiring the condition: β„‘(πœ›π‘˜)≀0.003.

Sinceπœ›π‘˜=(π›½πœƒ)π‘˜,(10.1) we also haveπœ›π‘˜=π‘Ÿπœƒβ„Ž.(10.2)

Since, we require πœ›π‘˜β‰€0.4, this implies that π‘Ÿ(πœƒβ„Ž)≀0.4. Also, we have πœƒβ„Žβ‰€1.76 and thus π‘Ÿβ‰€0.4/1.76.

The stability of the DRP scheme therefore satisfies the condition: π‘Ÿβ‰€0.23.

Using the approach just described in the preceding paragraph, the ranges of stability of some methods are obtained, namely, TAM-NEW, ZINGG-NEW, ZHUANG-NEW, LOCKARD-NEW, and BOGEY-NEW when combined with TAMMODIFIED. We also obtain the range of stability for the methods: ZINGG, ZINGG, ZHUANG, LOCKARD, and BOGEY when they are combined with the temporal discretisation scheme of Tam et al. [17]. The results are tabulated in Table 7. It is seen that the new combined spatial-temporal methods constructed using MIEELDLD have a slightly greater region of stability than the existing combined spatial-temporal methods.

tab7
Table 7: Region of stability for some combined spatial-temporal discretisation schemes.

11. Comparison of Some Metric Measures

Spatial Scheme of Tam and Webb [3]
The integrated error is defined as ξ€œ01.1||πœƒβˆ—||β„Žβˆ’πœƒβ„Ž2𝑑(πœƒβ„Ž).(11.1)
The quantity, |πœƒβˆ—β„Žβˆ’πœƒβ„Ž|2 is equivalent to |1βˆ’RPE|2 in a computational fluid dynamics framework. A plot of |1βˆ’RPE|2 versus RPE∈[0,2] is shown in Figure 14(a).

fig14
Figure 14: Plot of different metrics from Tam and Webb [3], Bogey and Bailly [16] and Appadu and Dauhoo [18].

Spatial Scheme of Bogey and Bailly [16]
In this case, the integrated error is defined as ξ€œπœ‹/2πœ‹/16||πœƒβˆ—||β„Žβˆ’πœƒβ„Žπœƒβ„Žπ‘‘(πœƒβ„Ž),(11.2) or ξ€œlnπœ‹/2lnπœ‹/16||πœƒβˆ—||π‘‘β„Žβˆ’πœƒβ„Ž(ln(πœƒβ„Ž)).(11.3) The quantity (|πœƒβˆ—β„Žβˆ’πœƒβ„Ž|)/πœƒβ„Ž is equivalent to (|1βˆ’RPE|)/RPE while |πœƒβˆ—β„Žβˆ’πœƒβ„Ž| is equivalent to |1βˆ’RPE|. Plots of (|1βˆ’RPE|)/RPE and |1βˆ’RPE|, both versus RPE∈[0,2], are shown in Figures 14(b) and 14(c).

Spatial Scheme Using MIEELDLD
A plot of eeldld=exp((||1βˆ’RPE|βˆ’(1βˆ’AFM)|)+exp(|1βˆ’RPE|+(1βˆ’AFM))βˆ’2) versus RPE∈[0,2] versus AFM∈[0,1] is shown in Figure 14(d).
We observe from Figures 14(a), 14(b), and 14(c) that the measure is zero when RPE=1 whereas, in Figure 14(d), the measure is zero provided RPE=1 and AFM=1.

12. Conclusions

In this work, we have used the technique of Minimised Integrated Exponential Error for Low Dispersion and Low Dissipation (MIEELDLD) in a computational aeroacoustics framework to obtain modifications to optimized spatial schemes constructed by Tam and Webb [3], Zingg et al. [14], Lockard et al. [13], Zhuang and Chen [15], and Bogey and Bailly [16], and also a modification to the optimized temporal scheme devised by Tam et al. [17] is obtained. It is seen that, in general, improvements can be made to the existing spatial discretisation schemes, using MIEELDLD. The new temporal scheme obtained using MIEELDLD is superior in terms of dispersive properties as compared to the one constructed by Tam et al. [17]. The region of stability has also been obtained. In a nutshell, we conclude that MIEELDLD is an efficient technique to construct high order methods with low dispersion and dissipative properties. An extension of this work will be to use the new spatial discretisation schemes and the novel temporal discretisation method constructed to solve 1-D wave propagation experiments and quantify the errors into dispersion and dissipation. Moreover, MIEELDLD can be used to construct low dispersive and low dissipative methods which approximate 2-D and 3-D scalar advection equation suited for computational aeroacoustics applications.

Nomenclature

𝐼=√(βˆ’1)
π‘˜:Time step
β„Ž:Spatial step
𝑛:Time level
𝛽:Advection velocity
πœƒβˆ—β„Ž:Numerical wavenumber
πœƒβ„Ž:Exact wavenumber
π‘Ÿ:cfl/courant number
π‘Ÿ=π›½π‘˜/β„Ž
𝑀:Phase angle in 1-D
𝑀=πœƒβ„Ž
πœ”:Exact angular frequency
πœ›:Effective angular frequency of time discretization scheme
RPE:Relative phase error per unit time step
AF:Amplification factor
AFM=|𝐴𝐹|
LDLD:Low Dispersion and Low Dissipation
IEELDLD:Integrated Exponential Error for Low Dispersion and Low Dissipation
MIEELDLD:Minimised Integrated Exponential Error for Low Dispersion and Low Dissipation.

Acknowledgment

This work was funded through the Research Development Programme of the University of Pretoria.

References

  1. S. Johansson, High Order Finite Difference Operators with the Summation by Parts Property Based on DRP Schemes, Division of Scientific Computing, Department Information Technology, Uppsala University, Sweden, 2007.
  2. J. Hardin and M. Y. Hussaini, Computational Aeroacoustics, Springer, New-York, NY, USA, 1992.
  3. C. K. W. Tam and J. C. Webb, β€œDispersion-relation-preserving finite difference schemes for computational acoustics,” Journal of Computational Physics, vol. 107, no. 2, pp. 262–281, 1993. View at Publisher Β· View at Google Scholar
  4. R. Hixon, β€œEvaluation of high-accuracy MacCormack-type scheme using benchmark problems,” NASA Contractor Report 202324, ICOMP-97-03-1997.
  5. G. Ashcroft and X. Zhang, β€œOptimized prefactored compact schemes,” Journal of Computational Physics, vol. 190, no. 2, pp. 459–477, 2003. View at Publisher Β· View at Google Scholar Β· View at Zentralblatt MATH Β· View at Scopus
  6. P. Roe, β€œLinear bicharacteristic schemes without dissipation,” SIAM Journal on Scientific Computing, vol. 19, no. 5, pp. 1405–1429, 1998. View at Publisher Β· View at Google Scholar Β· View at Zentralblatt MATH
  7. V. L. Wells and R. A. Renault, β€œComputing aerodynamically generated noise,” Annual Review Fluid Mechanics, vol. 29, pp. 161–199, 1997. View at Google Scholar
  8. C. K. W. Tam, β€œComputational aeroacoustics: issues and methods,” AIAA Journal, vol. 33, no. 10, pp. 1788–1796, 1995. View at Publisher Β· View at Google Scholar Β· View at Zentralblatt MATH Β· View at Scopus
  9. D. W. Zingg, β€œComparison of high-accuracy finite-difference methods for linear wave propagation,” SIAM Journal on Scientific Computing, vol. 22, no. 2, pp. 476–502, 2001. View at Publisher Β· View at Google Scholar Β· View at Zentralblatt MATH
  10. T. K. Sengupta, G. Ganeriwal, and S. De, β€œAnalysis of central and upwind compact schemes,” Journal of Computational Physics, vol. 192, no. 2, pp. 677–694, 2003. View at Publisher Β· View at Google Scholar Β· View at Zentralblatt MATH Β· View at Scopus
  11. W. De Roeck, W. Desmet, M. Baelmans, and P. Sas, β€œAn overview of high-order finite difference schemes for computational aeroacoustics,” in Proceedings of the International Conference on Noise and Vibration Engineering, pp. 353–368, Katholieke Universiteit Leuven, Belgium, September 2004.
  12. M. Popescu, W. Shyy, and M. Garbey, β€œFinite-volume treatment of dispersion-relation- preserving and optimized prefactored compact schemes for wave propagation,” Tech. Rep. UH-CS-05-03.
  13. D. P. Lockard, K. S. Brentner, and H. L. Atkins, β€œHigh-accuracy algorithms for computational aeroacoustics,” AIAA Journal, vol. 33, no. 2, pp. 246–251, 1995. View at Publisher Β· View at Google Scholar Β· View at Zentralblatt MATH Β· View at Scopus
  14. D. W. Zingg, H. Lomax, and H. Jurgens, β€œHigh-accuracy finite-difference schemes for linear wave propagation,” SIAM Journal on Scientific Computing, vol. 17, no. 2, pp. 328–346, 1996. View at Publisher Β· View at Google Scholar Β· View at Zentralblatt MATH
  15. M. Zhuang and R. F. Chen, β€œApplications of high-order optimized upwind schemes for computational aeroacoustics,” Aiaa Journal, vol. 40, no. 3, pp. 443–449, 2002. View at Publisher Β· View at Google Scholar Β· View at Scopus
  16. C. Bogey and C. Bailly, β€œA family of low dispersive and low dissipative explicit schemes for flow and noise computations,” Journal of Computational Physics, vol. 194, no. 1, pp. 194–214, 2004. View at Publisher Β· View at Google Scholar Β· View at Zentralblatt MATH Β· View at Scopus
  17. C. K. W. Tam, J. C. Webb, and Z. Dong, β€œA study of the short wave components in computational acoustics,” Journal of Computational Acoustics, vol. 1, no. 1, pp. 1–30, 1993. View at Publisher Β· View at Google Scholar
  18. A. R. Appadu and M. Z. Dauhoo, β€œThe concept of minimized integrated exponential error for low dispersion and low dissipation schemes,” International Journal for Numerical Methods in Fluids, vol. 65, no. 5, pp. 578–601, 2011. View at Publisher Β· View at Google Scholar Β· View at Scopus
  19. A. R. Appadu and M. Z. Dauhoo, β€œAn overview of some high order and multi-level finite difference schemes in computational aeroacoustics,” Proceedings of World Academy of Science, Engineering and Technology, vol. 38, pp. 365–380, 2009. View at Google Scholar Β· View at Scopus
  20. A. R. Appadu, β€œSome applications of the concept of minimized integrated exponential error for low dispersion and low dissipation,” International Journal for Numerical Methods in Fluids, vol. 68, no. 2, pp. 244–268, 2012. View at Publisher Β· View at Google Scholar Β· View at Scopus
  21. A. R. Appadu, β€œComparison of some optimisation techniques for numerical schemes discretising equations with advection terms,” International Journal of Innovative Computing and Applications, vol. 4, no. 1, pp. 12–27, 2012. View at Google Scholar
  22. A. R. Appadu, β€œInvestigating the shock-capturing properties of some composite numerical schemes for the 1-D linear advection equation,” International Journal of Computer Applications in Technology, vol. 43, no. 2, pp. 79–92, 2012. View at Google Scholar
  23. J. Wang and R. Liu, β€œA new approach to design high-order schemes,” Journal of Computational and Applied Mathematics, vol. 134, no. 1-2, pp. 59–67, 2001. View at Publisher Β· View at Google Scholar Β· View at Zentralblatt MATH
  24. L. L. Takacs, β€œA two-step scheme for the advection equation with minimized dissipation and dispersion errors,” Monthly Weather Review, vol. 113, no. 6, pp. 1050–1065, 1985. View at Publisher Β· View at Google Scholar Β· View at Scopus
  25. R. Liska and B. Wendroff, β€œComposite schemes for conservation laws,” SIAM Journal on Numerical Analysis, vol. 35, no. 6, pp. 2250–2271, 1998. View at Publisher Β· View at Google Scholar Β· View at Zentralblatt MATH
  26. M. Lukacova, β€œFinite volume schemes for multi-dimensional hyperbolic systems based on the use of bicharacteristics,” Applications of Mathematics, vol. 51, no. 3, pp. 205–228, 2006. View at Publisher Β· View at Google Scholar
  27. C. Kim, Multi-dimensional Upwind Leapfrog Schemes and their Applications, Ph.D. thesis, Aerospace Engineering, University of Michigan, 1997.
  28. J. Berland, C. Bogey, O. Marsden, and C. Bailly, β€œHigh-order, low dispersive and low dissipative explicit schemes for multiple-scale and boundary problems,” Journal of Computational Physics, vol. 224, no. 2, pp. 637–662, 2007. View at Publisher Β· View at Google Scholar Β· View at Zentralblatt MATH
  29. A. R. Appadu, M. Z. Dauhoo, and S. D. D. V. Rughooputh, β€œControl of numerical effects of dispersion and dissipation in numerical schemes for efficient shock-capturing through an optimal Courant number,” Computers & Fluids, vol. 37, no. 6, pp. 767–783, 2008. View at Publisher Β· View at Google Scholar