Implementation of a high-order finite difference scheme to model wave propagation

The large disparity between the time and length scales of an acoustically active flow field, and the ones of the resulting generated acoustic field, is an issue in computational aeroacoustics (CAA). Numerical schemes used to calculate the time and space derivatives in CAA should exhibit a low dispersion and dissipation error. This paper presents the implementation of a high-order finite difference scheme to model CAA problems. The numerical scheme consists of a sixth-order prefactored compact scheme coupled with the 4-6 Low Dispersion and Dissipation Runge-Kutta (LDDRK) time marching scheme. Onesided explicit boundary stencils are implemented and a buffer zone is used to eliminate spurious numerical waves. Results are in good agreement with a 1-D advection equation benchmark problem. This constitutes a preliminary validation for the scheme to explore and investigate the acoustic propagation of sound generated aerodynamically.


I. Introduction
The field of computational aeroacoustics (CAA) has grown rapidly during the last decade due to a rising interest in aeroacoustic phenomena driven by more stringent noise legislation and increasing environmental awareness.CAA focuses on the generation and propagation of aerodynamic noise.The inherent unsteady nature of aeroacoustic phenomena, the disparity in magnitude between mean and acoustic flow variables, and the necessity to model the high-frequency noise components place stringent demands on the numerical methods. 12Moreover, different aeroacoustic problems often exhibit totally different characteristics and behaviours.As a result, there is no single algorithm that can be used to simulate all problems with adequate resolution and accuracy. 3High-order finite difference schemes have low dispersion and dissipation properties that are required to reduce the number of grid points per wavelength, while still ensuring tolerable levels of numerical error.There are two main classes of high-accuracy finite difference schemes: explicit schemes and compact schemes. 4xplicit schemes employ large computational stencils for a given level of accuracy, while compact schemes use smaller stencils by solving for the spatial derivatives as independent variables at each grid point.While compact schemes are more accurate than the explicit stencil-equivalent ones, they have two disadvantages: first a matrix must be inverted to obtain the spatial derivative at each point; second, the boundary stencil has a larger effect on the stability and accuracy of the scheme.Hixon derived a class of prefactored compact schemes 5 that retain a high-order accuracy while using a even shorter stencil.These schemes reduce the computational effort while offering an easier boundary stencil specification than the non-prefactored schemes. 5he specific objective of this work is to implement and analyze a high-order accurate scheme for computational aeroacoustic applications.In this context, a scheme is regarded as high-order if it has a formal spatial and temporal accuracy higher than three.This scheme aims to resolve the acoustic near-field with adequate dispersion and dissipation properties for selected subsonic near-incompressible aerospace applications. 6The high-order finite-difference numerical scheme consists of a sixth-order compact scheme coupled with a 4-6 Low Dispersion Low Dissipation Runge-Kutta time marching scheme.The paper is organized as follows: Section II describes the numerical method and section III validates the scheme against a 1D advection equation problem. 7

II.A. Spatial discretization
A general compact derivative of a function f may be written as where D is the spatial derivative of the function f , and [B] is a square non-singular diagonally dominant matrix.Expanding this equation for a general pentadiagonal system at point i, For sixth-order accuracy: There are several interesting features of these compact schemes.The first is that the linear system contains tridiagonal matrices, the second is that the scheme requires two boundary stencils due to the tridiagonal nature of the scheme. 8

II.B. Sixth-order prefactored compact scheme
In a prefactored compact scheme, the non-prefactored derivative operator D i is split into a forward D F i and a backward D B i component, so that 8 where Let Substituting equations ( 7),( 8) in ( 5),( 6) and premultiplying and substituting in (4), it is obtained that Equation ( 9) requires that 8 which is true for constant-coefficient tridiagonal matrices.By comparing equations ( 1) and ( 9), the general form of bi-diagonal prefactored compact schemes can be derived: For sixth-order accuracy, The stencil on the right hand-side of the equation ( 3) is reduced from five points to three points, with the tridiagonal matrix on the left hand-side replaced by two independent bi-diagonal matrices.Figure 1 shows the dispersive characteristics and the relative amplitude error of the sixth-order prefactored compact scheme against the classical sixth-order 7-point stencil scheme (ST7) and the explicit fourth order 7-point Tam DRP scheme. 9The graph shows that the six-order prefactored compact scheme exhibits a better wavenumber resolution performance and has a lower error for a given number of grid points per wavelength.

II.C. Boundary stencils and buffer zone
1][12] The reason for this is that the error from the boundary stencil derivative can propagate well into the computational domain. 5o compute the state variables at the computational domain boundaries and wall points, explicit sixth-order one-sided derivative stencils are defined for the sixth-order compact scheme.To accomplish this, the Taylor series for the forward and backward interior derivatives was matched to the sixth order.The resulting boundary stencils for the backward sweep are and for the forward sweep: where the coefficients are given in the table 1.The numerical scheme possesses low dissipation and dispersion characteristics.Hence any inconsistencies due to numerical treatments at the boundaries will introduce errors or spurious wave reflections in the computation, which will eventually degrade the solution.To overcome this problem, a buffer zone is introduced in the numerical model.A buffer zone is a computational domain enclosing the physical domain in which damping is directly applied to the numerical solution vector Q at each time step: where Q is the damped solution vector and Q target is a given reference.The damping coefficient σ is defined as where l is the distance from the outer boundary of the buffer zone and L is the buffer zone width.Parameters α 2 (=1) and β 1 (=1.5) are used to control the damping coefficient σ which varies smoothly from zero at the interface between the buffer zone and the propagation region to a value, in this case 1.0, at the outer boundary.This specific damping coefficient provides a sufficient amount of damping to filter spurious reflections. 13

II.D. Time marching scheme
The time marching scheme employed is the 4-6 Low Dispersion and Dissipation Runge-Kutta (LDDRK) optimized scheme of Hu et al. 14 based on an optimization that minimizes the dissipation and dispersion errors for wave propagation.The 4-6 notation signifies a two-step marching cycle, where a four stage is used for the first time step and a six stage for the second time step in the cycle.Given the equation: where U t is the derivative in time and F x is the operator containing the discretization of the spatial derivatives, the LDDRK is given by: Step 1 Step 2 U n+1 = U n + h (6) Using this time marching method, the sixth-order prefactored compact scheme is stable up to a Courant number of 1.266.This scheme is fourth-order accurate in time for linear problems and second-order accurate for non-linear problems.This method has been previously tested and found to be accurate and stable for unsteady problems. 15The advantage of the LDDRK 4-6 method is that the optimization used reduces the error from one to three orders of magnitude for a given time step compared to the classical fourth-order accurate Runge-Kutta scheme. 14

III.A. 1-D Linear wave propagation
The implementation of the sixth-order prefactored compact scheme was tested against the benchmark problem of Category 1.1 taken from the First CAA Workshop. 7This problem asks for the solution at a nondimensional time=400 of the 1-D linear convection equation: where The boundary conditions of this problem are: at the inflow boundary, x = −20, U is equal to zero, at the outflow boundary, x = −450, U is calculated by the code using the one-sided explicit boundary stencils.When the buffer zone is applied (L = 10 on both sides), the boundary conditions are shifted to the outer region of the buffer zone. Figure 2(a) shows the predictions at the non-dimensional time ta ∞ /∆x = 200, obtained by using the sixthorder prefactored compact scheme without a buffer zone at the boundaries.The analytical solution at both boundaries is U = 0.The numerical solution displays a gradual departure from the analytical solution as U rises above zero approaching the left hand boundary.Figure 2(b) shows the same configuration but with a buffer zone at the boundaries.The introduction of the buffer zone has essentially eliminated the spurious rise in U close to the computational domain boundary.Figure 3 shows the predictions at the non-dimensional time ta ∞ /∆x = 400, calculated at two different Courant numbers.The numerical solution is in good agreement with the analytical solution sampled at every ∆x.The stability limit of the prefactored compact scheme allows to march the solution in time with a Courant number of 1.25, as shown in figure 3(b), without any significant degradation of the predictions as compared to a smaller time step, the results from which are shown in figure 3(a).

IV. Conclusion
A high-order finite difference scheme with low dispersion and dissipation characteristics has been coded up and tested.The scheme consists of a sixth-order prefactored compact scheme coupled with the 4-6 Low Dispersion and Dissipation Runge-Kutta (LDDRK) time marching scheme.One-sided explicit boundary stencils are implemented and the use of the buffer zone is described to eliminate spurious numerical waves.Results are in good agreement with a 1-D advection equation benchmark problem.This is a preliminary validation for the scheme.Future work will validate the scheme against test cases of increasing complexity to explore and investigate the acoustic propagation of sound generated aerodynamically.

Table 1 .
Coefficients sj and ej for one-sided boundary stencils.