A numerical validation of a high-order finite-difference compact scheme for computational aeroacoustics

Computing sound directly from unsteady flows featuring large-scale instabilities requires numerical methods that can resolve the small amplitude perturbations of sound embedded in the large-amplitude hydrodynamic flow unsteadiness. The perturbations due to sound need to be propagated with an acceptable dispersion and dissipation over significant distances to adequately resolve in space the noise near field. In this paper, a high-order compact scheme is detailed that promise to deliver the small dispersion and dissipation characteristics to directly model the noise and the unsteadiness in low Mach number flows. The performance of the scheme is assessed by cross-comparison against selected benchmark problems from the first workshop on benchmark problems in Computational Aeroacoustics (CAA) and additional test cases, including the reflection of an acoustic pulse into a corner, which is a challenging application due to the confluence of two wall conditions at the corner. Analytical solutions for the test cases are also produced to validate the predictions. The numerical results show that the scheme achieves comparable performances in computing these relatively simple two-dimensional test problems, with levels of dispersion and dissipation comparable to that in the published literature from similar computational aeroacoustic schemes. The results give confidence in developing these scheme to tackle more complex, three-dimensional noise producing flows, of interest for practical engineering applications.


I. Introduction
This work presents the validation of a high-order finite-difference scheme to model the wave propagation through homogeneous media, single and multiple wall reflection.The numerical method has been developed at the Institut de Mécanique des Fluides de Toulouse (IMFT), France, and it consists of an in-house sixth order compact finite-difference scheme.Fourth-order accuracy in time is obtained from a usual Runge-Kutta scheme.The predictions show that the numerical results are in very good agreement with the reference analytical solution.In performing such tests the guidelines of the first workshop on benchmark problems in CAA 1 have been followed as rigorously as possible, to provide a solid validation statement.The absence of spurious numerical features in areas of multiple wall reflections is a novel test to verify the absence of numerical artefacts where a more non-linear interaction of larger amplitudes occurs at the corner of the computational domain, which is usually a difficult area to treat with adequate numerical conditions.

II. Numerical Method
The numerical method solves the compressible, non-linear Euler equations with conservative variables in Cartesian coordinates. 2The method has been extended to solve the Navier-Stokes equations but, for the present work, the Euler code is used.

II.A. Spatial and temporal discretization
The spatial discretization is a compact sixth-order finite-difference scheme optimized in a five-point stencil.A progressive-regressive formulation by Kloker 2,3 generates numerical fluxes of alternate sign that are used in a classical fourth order Runge-Kutta scheme to time-march the flow and the acoustic field.At each temporal iteration, the forward-backward sequence is inverted to implicitly eliminate high-frequency numerical oscillations.Neither filtering nor the introduction of an artificial viscosity term are necessary to stabilize the scheme when coupled with the computational boundary conditions.This scheme is stable up to a Courant number of 1.11. 3 The forward and backward formulations to compute the derivative of the flux F in the x direction are given in table I, (Ia,Ib): ∂Fi−1,j ∂x Close to the computational domain boundaries, a one-sided five-point stencil sixth-order compact scheme is used at the second and the last but one points.The stencil of the one-sided scheme is from Table I The conservative variables fluxes F at the interior domain and at the boundaries are estimated by solving the ordinary differential equation : where A is a tridiagonal matrix containing the left-hand-side coefficients of equations 1-3, B is a pentadiagonal matrix containing the right-hand-side coefficients of equations 1-3 and F are the fluxes arranged as a matrix.The derivatives of the fluxes, DF , are determined by LU decomposition.
In order to obtain the derivatives of the fluxes in the whole computational domain, the vector F requires information from the the computational boundaries, i.e. the first and last points.As a consequence, the matrices A and B are completed with a boundary stencil, which is an explicit fifth-order scheme in a sixpoint stencil [cf.Table I, (IV) 3 ].The numerical values of the fluxes at the computational boundaries are not a direct output of equation 4, but are computed at later stage using specific boundary conditions.These schemes are an extension of the classical sixth-order centered schemes proposed by Lele.In the present computation, the schemes are applied in the i and j directions in operator split form. 4

II.B. Non-reflecting boundary conditions
Two different non-reflecting boundary conditions are implemented in the code and have been used in this study.The first one is a formulation based on the method of characteristics proposed by Giles. 5This formulation decouples the conservative Euler equations into characteristic Euler that represent vorticity, entropy and acoustic waves travelling at their respective speeds.To make the boundaries non-reflecting, the incoming waves are suppressed.Characteristic boundary conditions are essentially one-dimensional and are good for absorbing waves travelling normal to the boundary.However, the formulation presents problems when the travelling waves are oblique with respect to the boundary and corners are especially difficult to treat.The second method used to implement non-reflecting boundary conditions was developed by Tam. 6,7  uses the acoustic far-field asymptotic expression of the Euler equations to simulate a free-field acoustic radiation.This approach is formulated in cylindrical coordinates, where the origin is usually the centre of the noise source.The main advantage of this approach is that it is multi-dimensional and the problems at the corners are reduced due to the use of cylindrical coordinates.

II.C. Wall boundary conditions
Solid boundaries in this investigation are modelled as isothermal and no-slip.The formulation of an appropriate wall boundary condition in high-order schemes is still an open problem, since common formulations for low-order schemes might lead to numerical instabilities when applied to high-order schemes.In this study, two different boundary conditions are implemented and tested.The first one is a robust formulation proposed by Gloerfelt. 8The isothermal condition is directly imposed by calculating the density at the given temperature at the wall.The no-slip condition is applied to estimate the convective fluxes in conservative form normal to the wall: where subscripts (ρ, ρu, ρv, ρe) indicate the estimation of the conservative fluxes for the conservation of mass, x-momentum, y-momentum, and energy, respectively.In equation( 5) ∂v/∂y is obtained from the second-order approximation where v i,j = 0 is the flow-normal velocity at a horizontal wall.
The second method is a ghost cell method with no point on the wall.A ghost cell is placed outside the computational domain, at the same distance from the wall as the first interior point.The isothermal and no-slip conditions are used to estimate the flow state at the ghost cells, denoted by g.
where subscript w denotes a condition at the wall and r is the temperature recovery factor.The condition p g = p 1 is equivalent to imposing ∂p/∂y = 0 at the wall by a second order centred scheme.

III. Validation Path
A hierarchy of benchmark problems are selected to test the implementation of the high-order quadrature scheme.These test cases are of increasing complexity.The benchmark problems are taken from the first workshop on benchmark problems in CAA 1 and an additional test case, including the reflection of an acoustic pulse into a corner has been added.For all the simulations, the computational domain is an equidistant cartesian grid of 201 × 201 points, and the space increment in x and y is 1 m.

IV.A. 2-D acoustic, entropy and vorticity waves
Problem IV.A tests the effectiveness of the radiation boundary conditions, the inflow and outflow boundary conditions and the isotropy property of the computational algorithm in a complex problem.This is problem 1, category 3 of the first workshop on benchmark problems in Computational Aeroacoustics. 1 It consists on an acoustic wave, a vorticity wave and an entropy wave propagating in a M = 0.5 uniform flow with no solid boundary, as shown in Fig. 1(a).At t = 0 the initial flow conditions are Dimensional variables are used throughout.The reference values are p 0 = 10 5 Pa, T 0 = 298 K and u 0 = M x c 0 , c 0 is the ambient speed of sound, M x = 0.5 and ε = 0.01.The non-reflecting boundary conditions detailed in section II.B are applied to the computational domain boundaries.Using the boundary conditions by Tam 6,7 gave similar results to those from using the formulation by Giles. 5For clarity, only the results obtained with the Giles characteristics boundary conditions are included here.
Non-dimensional values for density ρ, pressure p, and time t are defined.The normalization is done using the following scales: length scale ∆x = ∆y, time scale ∆x/c 0 and density scale ερ 0 /γ.From these, the velocity and pressure are scaled by c 0 and ερ 0 c 2 0 .From these, the non-dimensional values of the perturbations ρ and p are: ρ = ρ − γ/ε and p = pγ − 1/ε.
Even though the computations are done in dimensional form, the final results are normalized and compared against the non-dimensional analytical solution given in the first workshop on benchmark problems in Computational Aeroacoustics: where 1/2 and J 0 () is a Bessel function of order 0. The evaluation of the integral has been done numerically at the University of Leicester with MATLAB R 7.3.0using the adaptive Lobatto quadrature with an absolute error tolerance of 1.0e − 9.The computed results for non-dimensional density ρ are shown in figure 2. All results are given at the non-dimensional times t = 30, 60 and 120. Figure 2(a) shows the non-dimensional density contours for the combined convected acoustic, vorticity and entropy waves of the benchmark problem IV.A at a non-dimensional computational time t = 30.Four non-dimensional contour levels are shown, in accordance to the benchmark problem solution guideline.These are ρ = (−0.02,0.01, 0.02, 0.04).The waves are shown approaching the right computational boundary.At this time, these waves have yet to interact with the boundary and are fully contained within the computational domain.The numerical prediction, shown by the black dashed line, overlaps the analytical solution, which is shown by the continuous grey line.There is no appreciable azimuthal distortion of any of the two waves, that is, the waves are essentially circular, showing that the numerical solution does not suffer from any appreciable degradation where the waves propagate at an angle with respect to the Cartesian mesh.A cross-section of the non-dimensional density distribution along the x axis is given to the right of the contour plot.The location of this cross-section is as specified in the first workshop on benchmark problems in CAA. 1 The numerical solution, shown by the the squares (2), is well-matched to the analytical solution, shown by the solid line.There is no appreciable departure between the two solutions at the right boundary of the plot, where the predictions are closest to the right computational boundary and are therefore more likely to start to suffer from numerical boundary effects.Figure 2(b) shows the predicted non-dimensional density distribution at the non-dimensional time t = 60.At this time, the acoustic, entropy, and vorticity waves impinge simultaneously against the right computational boundary.The same normalized contour levels are used as in Fig. 2(a).The numerical prediction (−−) overlaps the analytical solution (−).The cross-section of the density distribution along the x axis to the right of the contour map shows more details of the interaction between the three waves.As the waves coalesce onto one another at the right computational boundary, the density perturbation constructively interfere, creating a normalized density peak of about ρ ∼ 0.125.The numerical solution reports a maximum that is just a fraction lower than the analytical result, due to the maximum lying between two nodes of the discretized computational domain.
Figure 2(c) shows the non-dimensional density contours at the non-dimensional computational time t = 120.The same contour levels of Figs.2(a) and 2(b) are used.At this non-dimensional time, the waves have left the computational domain through the right, top and bottom computational boundaries.No reflection is shown inside the domain at the lowest non-dimensional density contour level of 0.01.The non-dimensional density distribution along the x axis, plotted to the right of the contour map, shows a good agreement between the numerical prediction and the analytical solution.
The right hand side plots of figures 2 (a,b,c) show a small density rise at the centre of the pulse.At t = 30, this rise is located at x = 15∆x.This is due to the numerical representation of the Gaussian pulse.Specifically, the pulse maximum lies between two nodes of the discretized computational domain.This feature advects rightwards with the Mach 0.5 flow, remaining reasonably centred with respect to the expanding Gaussian pulse.

IV.B. 2D acoustic wave against a wall
Problem IV.B tests the implementation of no-slip boundary conditions over a single solid wall to reflect a convecting pressure wave.This the problem 1, category 4 of the first workshop on benchmark problems in Computational Aeroacoustics. 1 It consists of an acoustic wave convected by a Mach 0.5 uniform flow that impacts against a wall located to the south of the perturbation source as shown in Fig. 1(b).
In this problem, the computational domain is −100 ≤ x ≤ 100, 0 ≤ y ≤ 200.The equations solved are the Euler equations and the initial flow conditions are where the reference values are p 0 = 10 5 Pa, T 0 = 298K and u 0 = M x c 0 where c 0 is the ambient speed of sound, the Mach number M x = 0.5 and ε = 0.01.Two boundary conditions have been used for the wall: Gloerfelt's boundary conditions and the ghost cells method.For the other boundaries, the non-reflecting characteristics boundary conditions of Giles have been used.Both solid boundary conditions gave similar results.For clarity, only the results obtained using ghost cells are shown.
The results have been normalized as in section IV.A and compared against the non-dimensional analytical solution given in the first workshop on benchmark problems in Computational Aeroacoustics 1/2 and J 0 () is a Bessel function of order 0. Figure 3 shows the results for density pertubations at the non-dimensional times t = 30, 60 and 120.At the earliest non-dimensional time t of 30, fig.3(a), the acoustic pulse has just impacted against the south wall.The early stages of reflected wave are described by the acoustic density contours map.Contour levels 0.01 and 0.05 have been selected in accordance to the guidelines of the first workshop on benchmark problems in CAA. 1 The analytical reference solution (−) has been overlaid on the numerical predictions (−−) and there is a consistent good overlap between the two.From the non-dimensional density field, a cross-section of the density perturbation distribution has been extracted along x = y line.The selection of the cross-section matches with the directives of this test case. 1 The density perturbation is non monotonic along this line, yet the numerical scheme is able to predict this (2) and follow the analytical benchmark solution (−) pretty accurately.
Figure 3(b) shows the evolution of the incident reflected acoustic waves before any wave goes through any of the non-reflecting computational boundaries.The agreement between the predicted density perturbations (−−) and the reference analytical solution (−) is good.The density perturbations cross-section to the right of the contour map shows two similar peaks, the peak at s = 50 is the pulse reflection from the wall and its shape is similar to the incident wave, which is centred at s ∼ 100.The numerical predictions (2) follow the analytical benchmark solution (−) accurately.
Figure 3(c) shows the acoustic pulse after it has crossed the right computational boundary.The overall agreement between the contours of predicted density perturbations (−−) and the reference analytical solution (−) is reasonably good.There is a confined area where the reflected pulse crosses the right computational boundary in which the predictions seem to display some departures from the reference analytical solution.This feature is further highlighted by the density perturbation cross-section along x = y to the right of fig.3(c).At s > 100, the predicted density perturbation minimum is slightly under-estimated by the numerical method and small-amplitude waves with the characteristic wave-number higher than that of the acoustic pulse are shown around the density perturbation minimum.The density perturbation maximum is instead well-captured by the numerical prediction (2).

IV.C. 2D acoustic wave against a corner
Problem IV.C tests the ability of the solid boundary conditions of section II.C to handle multiple pressure wave reflections into a corner.This test case is inspired by the previous ones, with the innovation of the presence of two walls, on the south and east boundaries, forming the corner.This problem has no incoming flow.The computational domain and initial conditions are shown in fig.1(c).The computational domain is −100 ≤ x ≤ 100, −100 ≤ y ≤ 100, the governing equations are the discrete Euler equations and the initial flow conditions are where the reference values are p 0 = 10 5 Pa, T 0 = 298K, u 0 = M x c 0 , c 0 is the ambient speed of sound, the Mach numbers are M x = M y = 0, and ε = 0.01.At the walls, the ghost cell method of section II.C is used.At the remaining boundaries, the non-reflecting characteristic boundary conditions of Giles have been used.The alternative formulation by Gloerfelt for a no-slip boundary condition was tested and similar results were obtained.
The results have been normalized as in section IV.A.An analytical solution has been derived by mirroring the solution in section IV.B about the right wall and by adding the original and the mirrored pressure and density fields.This gives 1/2 and J 0 () is a Bessel function of order 0.
Figure 4 shows the propagation of the acoustic wave at three different non-dimensional times t = 60, 120 and 180. Figure 4(a) shows the contours of non-dimensional density perturbation at the nondimensional time t = 60.Contour levels have been selected consistently with test case IV.A.These are ρ = (−0.02,0.01, 0.02, 0.04).The numerical predictions (−−) overlap the analytical reference solution (−).At this time, the acoustic wave has not yet reached any boundary and the solution simply reflects the good dispersive and dissipative properties of the numerical scheme in the interior of the computational domain.
To the right of the contours, the density perturbation distribution along the x axis is shown.The pattern is symmetric with respect to the origin of the Cartesian reference system, x = 0.The predicted distribution (2) matches the analytical reference result (−).
Figure 4(b) shows the non-dimensional density perturbation at the later time t = 120.The acoustic wave has impacted against both the south and the east walls, causing local reflections.These reflections are apart from one-another and the prediction behaves similarly to test case IV.B.The overlap between the predicted density perturbation contours (−−) and the analytical reference solution (−) is good, just like for test case IV.B.The density perturbation cross-section along the y = 0 plane to the right of the contour plot shows the profile of the reflected pulse.It is noticeable that whereas in Fig. 4(a) the right hand side normalized density perturbation features a minimum followed by a maximum, close to the east boundary, in Fig. 4(b) a density perturbation maximum is followed by a minimum, confirming that this is a reflection from incident wave.The overlap between the predicted (2) and the reference (−) distributions is very good.
Figure 4(c) shows the normalized density perturbation distribution at t = 180.At this time, the reflected waves from the south and east walls have intersected, resulting in the wave interference pattern shown in Fig. 4(c).The wave branch close to the corner (x, y) = (100, 0) results from the combined wall conditions at the corner.There is no evidence that such condition generated any spurious numerical feature as the numerical predictions (−−) well match the reference solution (−) along this wave branch.The reflected waves also constructively interfere close to the origin of the Cartesian reference system, at the centre of the computational domain.This interference seems to have been handled well by the numerical scheme, which displays no high-wavenumber rippling effect at the cross-over point.Where the reflected waves cross the farfield boundaries, the normalized density fluctuation contours exhibit a small departure from the analytical

V. Conclusion
The study presented the validation of a high-order finite-difference scheme to model the wave propagation through homogeneous media and single and multiple wall reflection.The predictions show that the numerical results are in very good agreement with the reference analytical solution.In performing such tests the guidelines of the the first workshop on bechmark problems in CAA 1 have been followed as rigorously as possible, to provide a solid validation statement.The absence of spurious numerical features in areas of multiple wall reflections is a novel test to verify the absence of numerical artefacts where a more non-linear interaction of larger amplitudes occurs at the corner of the computatinal domain, which is usually a difficult area to treat with adequate numerical conditions.Again the performance of the code is good.
The Navier-Stokes extension of the Euler code will be used for sound control purpose on a two-dimensional compressible cavity flow following the work carried out by Spagnoli 9, 10 and Wei and Freund 11 on an unstable shear layer.The optimal open loop control will consist on defining an area where sound is targeted to decrease, and an area where the actuators will act (on one of the cavity walls for wall control, inside the cavity for control modelling).The gradients of the functional to minimize are given by continuous adjoint Navier-Stokes equations solved on approximately the same computational domain.Then a conjugate gradient algorithm will provide the optimal control.The feedback control will be carried out for the same test case using reduced order modelling from DNS simulations and Galerkin projections.Different alternatives are to be explored (Proper Orthogonal Decomposition, instability eigen modes).Feedback control has been done experimentally, and it is aimed to do it computationally. 12