# A very fast high-order flux reconstruction for Finite Volume schemes for Computational Aeroacoustics

Luis Ramirez, Javier Fernández-Fidalgo, José Paris, Michaël Deligant, Sofiane Khelladi, Xesús Nogueira

<figure><img src="image_1.png" /></figure>

## ► To cite this version:

Luis Ramirez, Javier Fernández-Fidalgo, José Paris, Michaël Deligant, Sofiane Khelladi, et al. A very fast high-order flux reconstruction for Finite Volume schemes for Computational Aeroacoustics. Engineering with Computers, 2024, 41, pp.667-680. 10.1007/s00366-024-02039-2. hal-05444387

HAL Id: hal-05444387

https://hal.science/hal-05444387v1

Submitted on 6 Jan 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# A very fast high-order flux reconstruction for Finite Volume schemes for Computational Aeroacoustics

Luis Ramírez<fnref n="1" /> · Javier Fernández-Fidalgo<fnref n="2" /> · José París<fnref n="1" /> · Michael Deligant<fnref n="3" /> · Sofiane Khelladi<fnref n="3" /> · Xesús Nogueira<fnref n="1" />

Received: 6 February 2024 / Accepted: 24 July 2024

© The Author(s) 2024

## Abstract

Given the small wavelengths and wide range of frequencies of the acoustic waves involved in Aeroacoustics problems, the use of very accurate, low-dissipative numerical schemes is the only valid option to accurately capture these phenomena. However, as the order of the scheme increases, the computational time also increases. In this work, we propose a new high-order flux reconstruction in the framework of finite volume (FV) schemes for linear problems. In particular, it is applied to solve the Linearized Euler Equations, which are widely used in the field of Computational Aeroacoustics. This new reconstruction is very efficient and well suited in the context of very high-order FV schemes, where the computation of high-order flux integrals are needed at cell edges/faces. Different benchmark test cases are carried out to analyze the accuracy and the efficiency of the proposed flux reconstruction. The proposed methodology preserves the accuracy while the computational time relatively reduces drastically as the order increases.

**Keywords** High-order methods · Finite volume · Mean preserving moving least squares · Computational aeroacoustics

## 1 Introduction

Aeroacoustics is concerned with predicting the noise generated by unsteady flows. Computational Aeroacoustics (CAA) deals with the numerical simulation of these phenomena, which is a very complex problem [1–4]. Numerical simulation of aerodynamic flows is a mature technology which is widely used nowadays in multiple scientific and industrial environments. However, the application of numerical methods developed for aerodynamics to CAA is

not straightforward and presents a number of problems. Most of these problems are related with the low magnitude of acoustic waves and the importance of low amplitude pressure waves to the accurate representation of the aeroacoustic field, which is wider than in aerodynamics. Thus, the use of high-resolution, low-dissipation numerical schemes is mandatory in the numerical simulation of this kind of computations.

Historically, high-resolution finite differences schemes [5–8] have been the predominant methods in CAA. However, high-order Finite Volume (FV) solvers [9–11] and even meshless methods based on Riemann solvers [12] have also

Javier Fernández-Fidalgo, José París, Michael Deligant, Sofiane Khelladi and Xesús Nogueira have contributed equally to this work.

<figure><img src="image_3.png" /></figure>

Luis Ramírez
luis.ramirez@udc.es

Javier Fernández-Fidalgo
j.fernandez.fidalgo@upm.es

José París
jose.paris@udc.es

Michael Deligant
michael.deligant@ensam.eu

Sofiane Khelladi
sofiane.khelladi@ensam.eu

Xesús Nogueira
xesus.nogueira@udc.es

<footnote n="1">¹ Group of Numerical Methods in Engineering-GMNI, Center for Technological Innovation in Construction and Civil Engineering-CITEEC, Civil Engineering School, Universidade da Coruña, Campus de Elviña, 15071 A Coruña, Spain</footnote>

<footnote n="2">² Departamento de Ingeniería Geológica y Minera, ETSI de Minas y Energía, Universidad Politécnica de Madrid, Calle Ríos Rosas 21, 28003 Madrid, Spain</footnote>

<footnote n="3">³ Arts et Métiers Institute of Technology CNAM, LIFSE, 75013, Paris, France</footnote>
been proposed for CAA computations. The interest of Finite Volume or meshless methods relies on its geometrical flexibility, which is a very interesting feature in order to address real engineering problems, since the use of finite difference schemes in such problems is cumbersome when complex geometries are involved.

Finite volume schemes are the current standard in industry for numerical simulation of compressible flows. Very high-order Finite Volume methods fulfill the requirements for CAA, since they are very low-dissipative schemes and with a relatively low dispersion error. Thus, they are suitable candidates for CAA computations of engineering problems. In order to achieve a high order finite volume scheme it is mandatory to use a high order reconstruction and also an appropriate quadrature for the integrals. So the order of the numerical scheme will be lowest between the order of the reconstruction and the numerical integration. However, using a pointwise representation of the variables the formal order is limited to order 3 for unsteady problems unless use a pseudo-mass matrix formulation, which relates the meanwise and pointwise framework [10] at the price of an increase of computational cost. A different approach is to use a meanwise formulation where arbitrary order can be also achieved [23] avoiding the increase of cost and complexity when dealing with the pseudo-mass matrix. Unfortunately, increasing the convergence order of a finite volume scheme implies increasing the cost of the numerical scheme both in terms of memory requirements and computational time. A great amount of this increase is related to the number of integration points at the interfaces between cells. This increase in computational cost is especially important in the case of three-dimensional problems, and it is a limitation of the convergence order of the numerical scheme in practical applications.

In this work, we propose a new high-order flux reconstruction method in the context of finite volume schemes applied to linear problems. LR: However, it is important to note that the proposed methodology is also applicable to other schemes such as Discontinuous Galerkin methods. This methodology allow us to obtain large reductions in computation time preserving the accuracy of the original schemes. The savings in computation time obtained are greater as the order of the scheme increases. In particular, we apply the new schemes to the Linearized Euler Equations, which are widely used in the field of Computational Aeroacoustics. Different benchmark test cases are carried out to analyze the accuracy and the efficiency of the new flux reconstruction. It is shown that the proposed reconstruction method is very efficient and well suited in the context of very high-order schemes, where high-order numerical integration of the flux is needed at cell edges/faces.

The paper is structured as follows. In Sect. 2 the governing equations are presented. In Sect. 3 the formulation of

classical finite volume methods is presented, whereas the new formulation is presented in Sect. 4. Some benchmark CAA test cases are analyzed in Sect. 5, to show the accuracy and efficiency gain of the proposed methodology. Finally, some conclusions are drawn.

## 2 Governing equations

The three-dimensional Linearized Euler Equations written in conservative form read

$$ \partial_t U + \nabla \cdot \mathcal{F}(U) + \mathcal{H}(U) = S(U) \quad (1) $$

where $U$ and $\mathcal{F} = (\mathcal{F}_x, \mathcal{F}_y, \mathcal{F}_z)$ are the vector of the conservative variables and the inviscid flux vector, namely

$$ U = \begin{pmatrix} \rho' \\ \rho_0 u' \\ \rho_0 v' \\ \rho_0 w' \\ p' \end{pmatrix} \quad (2) $$

$$ \mathcal{F}_x = \begin{pmatrix} \rho' u_0 + \rho_0 u' \\ \rho_0 u' u_0 + p' \\ \rho_0 v' u_0 \\ \rho_0 w' u_0 \\ p' u_0 + \gamma p_0 u' \end{pmatrix}, \quad \mathcal{F}_y = \begin{pmatrix} \rho' v_0 + \rho_0 v' \\ \rho_0 u' v_0 \\ \rho_0 v' v_0 + p' \\ \rho_0 w' v_0 \\ p' v_0 + \gamma p_0 v' \end{pmatrix}, \quad (3) $$

$$ \mathcal{F}_z = \begin{pmatrix} \rho' w_0 + \rho_0 w' \\ \rho_0 u' w_0 \\ \rho_0 v' w_0 \\ \rho_0 w' w_0 + p' \\ p' w_0 + \gamma p_0 w' \end{pmatrix} $$

The mean stationary solution is $U_0 = (\rho_0, \rho_0 u_0, \rho_0 v_0, \rho_0 w_0, p_0)^T$. In the acoustic perturbations variables $U$, $\rho'$ is the fluctuating density, $u' = (u', v', w')^T$ the fluctuating velocity and $p'$ the fluctuating pressure. Moreover, $\mathcal{H}(U)$ contains terms related to the gradients of the mean flow, which are zero when the mean flow is uniform

$$ \mathcal{H}(U) = \begin{pmatrix} 0 \\ (\rho_0 \mathbf{u}' + \rho' \mathbf{u}_0) \cdot \nabla u_0 \\ (\rho_0 \mathbf{u}' + \rho' \mathbf{u}_0) \cdot \nabla v_0 \\ (\rho_0 \mathbf{u}' + \rho' \mathbf{u}_0) \cdot \nabla w_0 \\ (\gamma - 1) p' \nabla \cdot \mathbf{u}_0 - (\gamma - 1) u' \nabla p_0 \end{pmatrix}. \quad (4) $$

## 3 Finite volume discretization

LR: Given a system of conservation laws
$$ \partial_t U + \nabla \cdot \mathcal{F} = S, \qquad (5) $$

**LR:** the FV discretization integrates it inside each control volume $\Omega_I$

$$ \int_{\Omega_I} \partial_t U \, d\Omega + \int_{\Omega_I} \nabla \cdot \mathcal{F} \, d\Omega = \int_{\Omega_I} S \, d\Omega. \qquad (6) $$

**LR:** After applying Gauss theorem and the mean value of the conservative variable, $\hat{U}_I$, we obtain

$$ \Omega_I \partial_t \hat{U}_I + \int_{\Gamma_I} \mathcal{F} \cdot n \, d\Gamma = \int_{\Omega_I} S \, d\Omega \qquad (7) $$

**LR:** where $\hat{U}_I$ is the mean value of the conservative variables $U(x, t)$ at control volume $I$, $\Gamma_I$ is the boundary of control volume $I$ and $n$ is the unitary outward normal of the boundary.

**LR:** For an arbitrary control volume $I$ we can split the face/edge integral of equation (7) in $n_f$ integrals, so $\sum_{j=1}^{n_f} \Gamma_{I_j} = \Gamma_I$, where $n_f$ is the total number of edges/faces in 2D/3D of the control volume $I$. Therefore, Eq. (7) reads

$$ \Omega_I \partial_t \hat{U}_I + \sum_{j=1}^{n_f} \int_{\Gamma_{I_j}} \mathcal{F} \cdot n_j \, d\Gamma = \int_{\Omega_I} S \, d\Omega. \qquad (8) $$

In a Finite Volume framework, the hyperbolic/inviscid flux vector of Eq. (8) is computed using a suitable numerical flux, denoted as $\Theta(U^+, U^-, n_j)$, where $U^+$ and $U^-$ are the right and left Riemann states of the variables on face $j$. Introducing the numerical flux in Eq. (8) we obtain

$$ \Omega_I \partial_t \hat{U}_I + \sum_{j=1}^{n_f} \int_{\Gamma_{I_j}} \Theta(U^+, U^-, n_j) \, d\Gamma = \int_{\Omega_I} S \, d\Omega \qquad (9) $$

Standard high-order Finite Volume discretization approximates the first integral of Eq. (9) by means of an appropriate numerical quadrature. Here, we use a Gauss-Legendre quadrature. Then, the discrete approximation of the flux integral of (9) reads as

$$ \sum_{j=1}^{n_f} \int_{\Gamma_{I_j}} \Theta(U^+, U^-, n_j) \, d\Gamma \approx \sum_{j=1}^{n_f} \sum_{ig=1}^{n_g} (\Theta(U_K^+, U_I^-, n_j))_{ig} L_j W_{ig} \qquad (10) $$

where $L_j$ denotes the length in 2D or area in 3D of edge/face $j$, $n_g$ is the number of quadrature points for each edge/face and $W_{ig}$ is the corresponding quadrature weight for the quadrature point (ig) where the function $\Theta$ is evaluated. In Eq. (10) subindexes $K$ and $I$ refer to the cells that share the edge/face (see Fig. 1). To obtain a high-order accuracy method, it is necessary to estimate the flux at each integration point using a high-order reconstruction of the variables, such as

<figure><img src="image_4.png" /><figcaption>**Fig. 1** Schematic representation of the integration points over edge $\Gamma_I$ of control volumes $I$ and $K$</figcaption></figure>

the MLS reconstruction [9, 10, 13] k-exact reconstruction [14, 15], ADER schemes [16, 17], WENO reconstruction [18–21] or recently the CWENO reconstruction [22], among others.

One of the key points in the development of a high-order Finite Volume method is that the order of the quadrature of the numerical integration must be, at least, of the same order of the FV method. For example, a $n^{th}$-order FV method using a Gauss-Legendre quadrature requires, at least, $(n+1)/2$ integration points for each edge in 2D and $((n+1)/2)^2$ integration points for each face in 3D.

In Fig. 1, we plot, as an example, three quadrature points along the edge $\Gamma_I$ for control volumes $I$ and $K$. For each quadrature point we need to compute the reconstructed variables $U_I^-(x_{ig})$ and $U_K^+(x_{ig})$ to calculate the numerical flux of Eq. (9). As explained earlier, the order of the reconstruction must be, at least, of the same order of the desired FV method.

In this work, we use the Mean Preserving Moving Least Squares, proposed by the authors in [23] which employs the information of the surrounding control volumes to high-order reconstruct the variables inside each control volume. For all control volumes, each conservative variable is approximated by

$$ U_I(x) = \hat{p}^T(x)\alpha = \sum_{k=1}^{n_x} N_k^I(x) \hat{U}_k \qquad (11) $$

where $N_k^I(x)$ are the MP-MLS shape functions, $n_x$ is the number of cells of the stencil and $\hat{p}^T(x)$ is the polynomial basis and is defined in 1D as

$$ \hat{p}_j(x) = \frac{1}{\Omega_I} \int_{\Omega_I} \hat{p}_j(x) d\Omega = \frac{1}{\Omega_I} \int_{\Omega_I} \left( \frac{x-x_I}{h} \right)^{j-1} d\Omega \qquad (12) $$

where $\alpha$ is the vector of reconstruction coefficients. These coefficients are calculated by minimizing an error functional of weighted residual taking into account the geometric information and mean values of the surrounding control volumes. The functional reads
$$J(\alpha(\mathbf{x})) = \int_{\mathbf{x}_k \in \Omega_l} W(\mathbf{x} - \mathbf{x}_k, h) \left[ \hat{U}(\mathbf{x}_k) - \frac{1}{\Omega_k} \int_{\Omega_k} \mathbf{p}^T(\mathbf{x}_k) \alpha(\mathbf{x}) d\Omega_k \right]^2 d\Omega_k. \quad (13)$$

where the $W(\mathbf{x} - \mathbf{x}_k, h)$ refers to a kernel or weight function which is centered on a point $\mathbf{x}$, and weights the importance in the approximation of the points inside the stencil. The stencil is determined by $h$ (the smoothing length), leading to a compact point cloud and must be at least the number of elements of the polynomial basis. The resulting approximation, $U_l(\mathbf{x})$ can be written in terms of the MP-MLS shape functions, $N(\mathbf{x})$ as seen in Eq. (11). It is important to remark that the number of points in the stencil increases as the order of the reconstruction is increased. Full details of the Mean Preserving Moving Least Squares can be found in [23].

Using this approach, for each integration point, $ig$ in Eq. (10), we compute the high-order reconstruction of the Riemann states as

$$\hat{U}_I^-(\mathbf{x}_{ig}) = \sum_{k=1}^{n_x} N_k^I(\mathbf{x}_{ig}) \hat{U}_k \quad (14)$$

$$\hat{U}_K^+(\mathbf{x}_{ig}) = \sum_{k=1}^{n_x} N_k^K(\mathbf{x}_{ig}) \hat{U}_k \quad (15)$$

Then, the discrete approximation of Eq. (9) reads

$$\Omega_l \partial_t \hat{U}_l + \sum_{j=1}^{n_f} L_j \sum_{ig=1}^{n_g} [\Theta(\hat{U}_K^+, \hat{U}_I^-, n_j)]_{ig} W_{ig} = \int_{\Omega_l} S \, d\Omega \quad (16)$$

Notice that Eq. (16) is the state of the art discretization to obtain a high-order finite volume scheme, where $\hat{U}_K^+$ and $\hat{U}_I^-$ are the high-order reconstruction of the variables.

## 4 Flux reconstruction through face average integration method

From Eq. (16), it becomes evident that as the order of the FV scheme is increased, the computational cost largely increases, specially in 3D. This is mainly due to the increase in the number of integration points needed to preserve high-order accuracy. As previously commented, for a Gauss-Legendre numerical face-integral in 3D, we need at least $((n+1)/2)^2$ integration points to preserved the order of a method of $n_{th}$ order. To reduce this large increase in cost, we propose to compute the face average values of the variables using a face averaged shape function as

$$\hat{U}_I = \frac{1}{L_j} \int_{\Gamma_j} U(\mathbf{x}, t) \, d\Gamma = \frac{1}{L_j} \int_{\Gamma_j} N^I(\mathbf{x}) \hat{U} \, d\Gamma = (\hat{N}^I)^T \hat{U} \quad (17)$$

where $\hat{N}^I$ is the mean value over the edge $\Gamma_j$ of the MP-MLS shape functions that reconstruct the variables inside the control volume $I$. Each component of the mean shape function is computed as

$$\hat{N}^I_j = \frac{1}{L_1} \int_{\Gamma_1} N_j^I(\mathbf{x}) \, d\Gamma \approx \sum_{ig=1}^{n_g} [N_j^I]_{ig} W_{ig} \quad (18)$$

Notice that the edge/face integral is computed for the MP-MLS shape functions instead of the numerical flux in Eq. (16). Since MP-MLS shape functions are independent of the value of the variables and, therefore, for fixed grids, they can be computed only once at the preprocessing stage. Therefore, we can pre-compute the value of $\hat{N}$ for each edge/face of the mesh. With this procedure, we are able to compute face-averaged values of the variables by using the following equation

$$\hat{U}_I = (\hat{N}^I)^T \hat{U}; \quad \hat{U}_I^+ = (\hat{N}^I)^T \hat{U} \quad (19)$$

This integration procedure is specially relevant for linear equations, such as the Linearized Euler Equations, where it is verified that

$$\Theta(\hat{U}_K^+, \hat{U}_I^-, n_j) = \frac{1}{L_j} \int_{\Gamma_j} \Theta(\hat{U}_K^+, \hat{U}_I^-, n_j) \, d\Gamma \quad (20)$$

And therefore, the discrete equation can be written as

$$\Omega_l \partial_t \hat{U}_l + \sum_{j=1}^{n_f} L_j \Theta(\hat{U}_K^+, \hat{U}_I^-, n_j) = \int_{\Omega_l} S \, d\Omega \quad (21)$$

From the point of view of coding, solving Eq. (21) is more efficient than solving Eq. (16), as schematically represented in the comparison between algorithms 1 and 2. Note that with this procedure we avoid the loop at integration points for each time step, since we have previously computed that loop at the preprocessing stage. LR: The total number of flux evaluations is related to the total number of integration points in the domain, that is for example, on a structured domain with a discretization of $N_x \times N_y$ elements, $(N_x + 1) \times (N_y + 1) \times \text{ngauss for the traditional integration scheme and } (N_x + 1) \times (N_y + 1) \text{ evaluations for the Face Average Integration Method. Thus, since the number of integration points is related to the order of our numerical method, the expected time savings will become larger as we increase the order of our numerical method.}$
**Algorithm 1** Computation of the flux using the traditional integration procedure reconstruction through Face Average Integration Method.

```
for icell = 1 to ncell
    flux(icell) = 0
    for if = 1 to nfaces(icell)
        for ig = 1 to ngauss
            flux(icell) = flux(icell) + Θ(Û<sub>K</sub><sup>+</sup>, Û<sub>I</sub><sup>-</sup>, n(if)) · W(ig) · L(if)
        end for
    end for
end for
```

**Algorithm 2** Computation of the flux using the Face Average Integration Method.

```
for icell = 1 to ncell
    flux(icell) = 0
    for if = 1 to nfaces(icell)
        flux(icell) = flux(icell) + Θ(Û<sub>K</sub><sup>+</sup>, Û<sub>I</sub><sup>-</sup>, n(if)) · L(if)
    end for
end for
```

We remark that the proposed Face Average Integration Method can be also applied to other mean-preserving reconstruction procedures. The key point is to be able to reformulate the reconstruction procedure, as expressed in Eq. (11), using shape functions.

# 5 Numerical results

This section presents the numerical results for several benchmark problems aimed at assessing the efficiency and performance of the proposed Face Average Integration Method, while keeping the same level of accuracy than the traditional scheme. In all the examples, the chosen time integration scheme is a TVD third-order Runge-Kutta [24]. **LR:** In this work, the numerical flux, $\Theta(Û_K^+, Û_I^-, n(if))$, is discretized using Roe's flux difference splitting [25] or the Rusanov's HLL flux [26], but we want to remark that the proposed approach is also valid for other numerical fluxes.

## 5.1 One-dimensional linear wave equation

In the first test case we analyze the benchmark problem of the First ICASE/LaRC Workshop on Benchmark Problems

in Computational Aeroacoustics [27]. The aim of this example is to test the accuracy and the order of convergence of the method. The following one-dimensional linear advection equation is solved

$$ \frac{\partial u}{\partial t} + a \frac{\partial u}{\partial x} = 0 \qquad (22) $$

where $a$ is the phase velocity. For simplicity, in this example, it has been taken $a = 1$. The initial condition is

$$ u(x, t = 0) = 0.5 \exp \left[ -\ln(2) \left( \frac{x}{3} \right)^2 \right] \qquad (23) $$

Equation (22) with initial condition (23) has the following analytic solution

$$ u(x, t = 0) = 0.5 \exp \left[ -\ln(2) \left( \frac{x - at}{3} \right)^2 \right] \qquad (24) $$

The domain is $-20 \le x \le 450$, and periodic boundary conditions are imposed. The domain is discretized with a set of uniform meshes of: 256, 1024, 2048, 4096, 8192, 16384, 32768 and 65536 control volumes. Once the solution is obtained the $L_2$-norm of the error is computed. **LR:** Since the time discretization is third order and the spatial discretization is arbitrary, following [28] we scale the $\Delta t$ to ensure that errors due the temporal integration converge
<figure><img src="image_5.png" /><figcaption>**Fig. 2** One dimensional linear wave problem. $L_2$ error norm of the variable for different convergence orders and spatial discretizations</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Fig. 3** One dimensional linear wave problem. Comparison of the numerical solutions $t = 400$ (left) and $t = 5100$ (right) on a 1024 control volume mesh obtained using different convergence orders</figcaption></figure>

at least with the same order as the spatial discretization. The numerical flux employed is an arbitrary-order upwind scheme, which in this test case is equivalent to the Rusanov flux. The obtained numerical results for different orders of accuracy are shown at $t = 400$ in Fig. 2.

The formal order of accuracy is achieved for the different orders tested. In Fig. 3 a) the numerical solution on a 1024 control volume mesh at $t = 400$ and b) at ($t = 5100$, which corresponds to 10 periods in the domain, are plotted. As expected,

the numerical method converges to the analytical solution as the order is increased, even for long time simulations.

## 5.2 2D-Gaussian acoustic pulse

The second test case is devoted to analyze the accuracy and efficiency of the proposed methodology in two dimensions. The two-dimensional Gaussian Acoustic pulse test case [5] is

<figure><img src="image_7.png" /></figure>
<figure><img src="image_8.png" /><figcaption>Fig. 4 2D-Gaussian acoustic pulse. Pressure field at $t = 0, t = 5$ and $t = 10$ with a first order finite volume method (left) and a tenth order FV- MP-MLS method (right)</figcaption></figure>

commonly used in the literature as a benchmark case [10–12, 29, 30] for CFD and CAA codes.

The computational domain is a square $\Omega = [-20, 20] \times [-20, 20]$ with periodic boundary conditions. The initial condition for the fluctuating variables is defined as

$$ \begin{aligned} \rho'(x, y, t = 0) &= e^{-\kappa r^2}, & u' &= v' = 0, \\ p'(x, y, t = 0) &= e^{-\kappa r^2} \end{aligned} \qquad (25) $$
<figure><img src="image_9.png" /><figcaption>**Fig. 5** 2D-Gaussian acoustic pulse. Pressure and velocity perturbation profiles at $x = 0$ obtained at $t = 10$ with a first order finite volume method and a tenth order FV-MP-MLS method compared with the analytical solution</figcaption></figure>

<figure><img src="image_10.png" /></figure>

<figure><img src="image_11.png" /><figcaption>**Fig. 6** 2D-Gaussian acoustic pulse. $L_2$ norm of the density for different orders of accuracy and spatial discretizations</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Fig. 7** 2D-Gaussian acoustic pulse. Computational times for both the face averaged flux integration and the classical Gauss Legendre integration</figcaption></figure>

where $\kappa = \ln(2)/b^2$, $b = 2$ and $r$ denotes the distance to the initial location of the pulse, which is located at $(x_0, y_0) = (0, 0)^T$. In this test case, the mean velocity is set to zero. The numerical solution is computed using the Rusanov numerical flux at $t = 10$.

First, in Fig. 4 the pressure field solution at times $t = 0$, 5 and 10 is plotted for a first-order finite volume method (left column) and a tenth-order FV-MP-MLS method (right column) on a mesh with 4096 control volumes. Notice that

the solution obtained with the first-order dissipates as the time increases.

The pressure and velocity field profiles at $x = 0$ obtained at time $t = 10$ are plotted in Fig. 5 and compared with the analytical solution. As expected, it becomes evident the increase of accuracy as the order of the numerical method is increased.

Second, we use this test case to check and compare the formal order of accuracy using both flux reconstruction

<figure><img src="image_13.png" /></figure>
approaches. To do so, the computational domain is discretized with a set of structured meshes, (256, 1024, 2048, 4096, 8192, 16384, 32768 and 65536 cells) and the solution is computed for each mesh with different orders of accuracy. **LR:** In order to avoid the influence of the time integration scheme, we have chosen a fixed time step of $\Delta t = 0.001$ to compare the computational times. Once the solution is obtained at $t = 10$ the $L_2$-norm of the density error is computed. The numerical results are plotted in Fig. 6, where no difference in terms of accuracy have been found between both flux reconstruction approaches. For convenience, we only plot the results obtained with the new flux reconstruction approach.

Once the accuracy and formal order of the schemes are verified, we compare the efficiency of the proposed Face Averaged Flux Integration method with the classical numerical integration method for different orders of accuracy. For a given mesh resolution of 4096 control volumes, the obtained results are summarized in Fig. 7, where the times are normalized by the time required for the second order scheme with the proposed methodology.

As the convergence order of the scheme is increased, the classical numerical integration of the flux requires more integration points for every edge of each control volume and therefore the computational time increases. However, the Face Averaged Flux Integration requires, for the same accuracy and order, lower computing times than the classical integration approach.

<figure><img src="image_14.png" /><figcaption>**Fig. 8 LR:** 2D-Gaussian Acoustic pulse. $L_2$ norm of the density for different orders of accuracy and the computational times for both the Face Averaged Flux Integration and the classical Gauss Legendre integration using the a mesh resolution of 4096</figcaption></figure>

Moreover, we note that, using the Face Averaged Flux Integration method, the computational time is only slightly increased as the convergence order is increased. This increasing is due to the larger stencil required for higher-order approaches. Thus, with the proposed methodology, the computational required for the seventh order scheme is less than two times larger than that required for the second order scheme, whereas using the classical approach, the time required is around eighth times larger. The tenth order scheme requires around four times more computing time than the second order scheme. With the classical scheme, the time required is 20 times larger. **LR:** For a fixed mesh of 4096 control volumes, we can plot the $L_2$-norm of the density error and the computational times required to obtain them with different order of accuracy. It can be observed in Fig. 8 the large computational time savings achieve with the proposed methodology.

We have verified that the new methodology keeps the accuracy of the original scheme with large savings of computational time, which becomes greater as the order increases.

## 5.3 3D-Gaussian acoustic pulse

In this example, we extend the previous numerical test case to 3D. The computational domain is $\Omega = [-20, 20] \times [-20, 20] \times [-20, 20]$. The initial position of the pulse is $(x_p, y_p, z_p) = (0, 0, 0)^T$ and it is defined as

<figure><img src="image_15.png" /><figcaption>**Fig. 9 3D-Gaussian acoustic pulse.** Computational times for both the face averaged flux integration and the classical Gauss Legendre integration</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>Fig. 10 Monopole source radiation in uniform (M = 0.5). Computational Mesh (left) and the absorbing layer surrounding the computational domain (right)</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>Fig. 11 Monopole source radiation in uniform (M = 0.5) mean flow at time t = 5400. Comparison of the pressure field on the computational domain using the proposed scheme for second (left), fourth (middle) and sixth (right) convergence orders</figcaption></figure>

$$ \rho'(x, y, z, t = 0) = e^{-\kappa r^2}, \quad u' = v' = w' = 0, \\ p'(x, y, z, t = 0) = e^{-\kappa r^2} \qquad (26) $$

where $\kappa = \ln(2)/b^2$, $b = 2$ and r denotes the distance to the initial location of the pulse. The mean velocity is set to zero. In order to compare the efficiency of the Face Averaged Flux Integration, the mesh resolution is fixed for all orders in 32³. The computational times for the different orders of accuracy are represented in Fig. 9, where the times are normalized by the time required for the second order scheme with the proposed methodology.

As expected, the computational time savings increases as the order of convergence increases. Notice the in this test case, the eight order method with the new approach is more than 20 times faster than the classical integration procedure.

<figure><img src="image_18.png" /><figcaption>Fig. 12 Monopole source radiation in uniform (M = 0.5) mean flow at time t = 5400. Pressure contours obtained using the proposed scheme for second, fourth and sixth convergence orders</figcaption></figure>
## 5.4 Monopole source radiation in uniform flow

This test case is commonly used to check the ability of the method to simulate long distance wave propagation. Here, we also use this test case to check if the Face Averaged Flux Integration method is able to keep the accuracy and efficiency shown in Cartesian meshes, when unstructured meshes are used to discretize the computational domain. We use the setup proposed in [31, 32]. A monopolar radiation source is placed in an uniform subsonic ($M = 0.5$) flow. It is located at $(x_p, y_p)^T = (0, 0)^T$ and is defined as

$$S_p = \frac{1}{2} \exp\left(-\ln(2)\frac{r^2}{2}\right) \sin(\omega t) \times [1, 0, 0, 1]^T \quad (27)$$

where $\omega = 2\pi/30$ is the angular frequency, $r$ denotes the distance to the location of the source and $t$ is the time coordinate. The wave length is $\lambda = 30$, and the computational domain is a circle with radius $r = 100$. The source term is made dimensionless with $[\rho_\infty c_\infty / \Delta x, 0, 0, \rho_\infty c_\infty^3 / \Delta x]^T$. An absorbing layer has been placed surrounding the computational domain to avoid spurious reflections at boundaries [31]. LR: The Roe numerical flux is used with a fixed time step of 0.5 for each order. Figure 10, shows the mesh resolution used in this problem and the absorbing layer around the computational domain.

In order to evaluate the long term accuracy of the method, we let the computations running until $t = 5400$, that corresponds to 180 periods of the source. Two acoustic waves propagate upstream and downstream of the source. Due to the effect of the mean flow, the apparent wavelength is modified. Thus, it is different upstream ($\lambda_1 = (1 - M)\lambda$) and downstream ($\lambda_2 = (1 + M)\lambda$) of the source (Fig. 11).

Figure 12 shows the results for the pressure profile along the axis $y = 0$. It is observed that the solution obtained for the proposed method is in agreement with the analytical solution, reproducing the different apparent wavelength upstream and downstream of the source.

Concerning the efficiency of the Face Averaged Integration Method, in Table 1 the CPU time is compared for different orders of reconstruction with the classical face

<figure><img src="image_19.png" /><figcaption>Fig. 13 Monopole source radiation in shear mean flow at time $t = 190$. Pressure contours obtained using the proposed scheme for sixth-order</figcaption></figure>

integration. The same accuracy has been obtained using both integration approaches.

For the given mesh, the fourth-order method using the Face-Averaged Integration Method proposed in this work obtains a solution that matches the analytical one 12, requiring a computation time comparable to that used by the second-order scheme using the standard integration method 1.

## 5.5 LR: Monopole source radiation in shear mean flow

LR: In this test case we analyze the proposed methodology when it is applied in presence of a non-uniform mean flow [32–35]. The initial condition is a monopolar radiation source is located at $(x_p, y_p)^T = (0, 0)^T$ and is defined as

<table><caption>**Table 1** Monopole source radiation in uniform flow: CPU time comparison between the Standard Face Integration Method using a Gauss Legendre quadrature and the Face Averaged Integration Method</caption><thead><tr><th></th><th>Standard integration method</th><th>Face averaged integration method</th></tr></thead><tbody><tr><th>Order</th><td>CPU time (s)</td><td>CPU time (s)</td></tr><tr><th>2th order</th><td>5.83</td><td>5.80</td></tr><tr><th>4th order</th><td>12.50</td><td>5.97</td></tr><tr><th>6th order</th><td>36.95</td><td>9.64</td></tr></tbody></table>

<table><caption>**Table 2** Monopole source radiation in shear mean flow: CPU time comparison between the standard face integration method using a Gauss Legendre quadrature and the face averaged integration method</caption><thead><tr><th></th><th>Standard integration method</th><th>Face averaged integration method</th></tr></thead><tbody><tr><th>Order</th><td>CPU time (s)</td><td>CPU time (s)</td></tr><tr><th>2th order</th><td>16.18</td><td>15.45</td></tr><tr><th>4th order</th><td>35.23</td><td>16.36</td></tr><tr><th>6th order</th><td>76.84</td><td>28.00</td></tr></tbody></table>
<figure><img src="image_20.png" /><figcaption>Fig. 14 Vortex convection in a uniform mean flow. Schematic representation of the test case setup</figcaption></figure>

$$S_p = \frac{1}{2}\exp\left(-\ln(2)\frac{r^2}{2}\right)\sin(\omega t) \times [1, 0, 0, 1]^T \quad (28)$$

**LR:** where $\omega = 2\pi/30$, $r$ denotes the distance to the location of the source and $t$ is the time coordinate. The wave length is $\lambda = 30$, and the computational domain and mesh are the same as the previous test case, a circle with radius $r = 100$. The mean flow is defined as

$$U_0 = (1, M_x \tanh(y/2), 0, 1/\gamma) \quad (29)$$

**LR:** where $M_x = 0.125$. The numerical solution is computed with both approaches at time $t=190$ using the 6th-order

scheme. The obtained pressure field is shown in Fig. 13 obtained with the Face Average Integration Method. The obtained results are in agreement with those presented in [34, 35].

**LR:** No discrepancies has been observed with the solution obtained with the standard integration scheme. The computational times for both methods are shown in Table 2.

## 5.6 LR: Vortex convection in a uniform mean flow

**LR:** In this last test case, we analyze the convection of a vortex in a uniform mean flow [31]. The aim of this test case is to analyze the proposed integration scheme in acoustic problems involving fluctuations in the velocity. The initial condition is defined as:

$$ \begin{aligned} \rho' &= 0, & u' &= y\epsilon\exp(-\ln(2)\frac{r^2}{b^2}), \\ v' &= -x\epsilon\exp(-\ln(2)\frac{r^2}{b^2}), & p' &= 0 \end{aligned} \quad (30) $$

<table><caption>Table 3 Vortex convection in a uniform mean flow: CPU time comparison between the standard face integration method using a Gauss Legendre quadrature and the face averaged integration method</caption><thead><tr><th></th><th>Standard integration method</th><th>Face averaged integration method</th></tr></thead><tbody><tr><td>Order</td><td>CPU time (s)</td><td>CPU time (s)</td></tr><tr><td>2th order</td><td>5.30</td><td>5.04</td></tr><tr><td>4th order</td><td>15.44</td><td>6.94</td></tr><tr><td>6th order</td><td>42.39</td><td>12.22</td></tr></tbody></table>

<figure><img src="image_21.png" /><figcaption>Fig. 15 Vortex convection in a uniform mean flow. Comparison of the residual fluctuating pressure for the standard and the proposed scheme using the 6th-order method</figcaption></figure>
**LR:** where $b = 5$ and $\epsilon = 0.03$. The vortex defined in (30) is convected with a Mach $M_x = 0.5$ on a domain of size $200 \times 200$ with periodic boundary conditions. An schematic representation of the test case is shown in Fig. 14. **LR:** Following [31], we evaluate the magnitude of the spurious noise generated by the vortex to verify that the proposed integration method does not introduce additional spurious noise and that it obtains the same solution as the traditional integration method. To measure it, we will use the residual fluctuating pressure, defined as

$$L_p = \sqrt{\frac{1}{N} \sum_{i=1}^{N} (p - p_0)^2} \qquad (31)$$

where $N$ is equal to the number of cells. In this example we have used a structured mesh of 4096 cells to discretize the computational domain. The vortex is convected for a complete period until $t = 400$. The obtained residual fluctuating pressure is plotted in Fig. 15, for a sixth-order method with both schemes. No discrepancies have been observed between both solutions. Concerning the computational time, the obtained CPU times are shown in Table 3 for both approaches and different orders. The time step is fixed for all orders in $\Delta t = 0.05$ and the Rusanov numerical flux is used.

## 6 Conclusions

In this work we have proposed a new and very efficient discrete flux reconstruction method for linear equations, which is well suited in the context of very high-order schemes, where high-order flux integrals are needed at cell edges/ faces. We have applied the proposed methodology to the Linearized Euler Equations using the MP-MLS reconstruction procedure on a FV context. **LR:** We remark that the proposed Face Average Integration Method is not restricted to the MP-MLS reconstruction in FV schemes. In fact, it can be applied to other reconstruction procedures and other numerical schemes to reduce the required edges/face integration points, such as Discontinuous Galerkin Methods.

Different benchmark test cases in the field of Computational Aeroacoustics are carried out to analyze the accuracy and the efficiency of the proposed flux reconstruction. It is seen that the proposed methodology greatly reduces the computing times as the convergence order of the numerical method is increased, specially in three-dimensional computations and it is able to deal with unstructured meshes. In particular, using the proposed methodology, we are able to reduce up to 20 times the computing time of the eight-order scheme in three dimensional computations. This has important consequences for the application of the method to industrial problems, which will be the next step of our research.

**Acknowledgements** L. Ramírez and X. Nogueira acknowledge the support provided by the [Grant PID2021-125447OB-I00] funded by MCIN/AEI/ 10.13039/501100011033 and by “ERDF A way of making Europe”, and the funds by [Grant TED2021-129805B-I00] funded by MCIN/AEI/ 10.13039/501100011033 and by the “European Union NextGenerationEU/PRTR”. They also acknowledge the funding provided by the Xunta de Galicia (Grant #ED431C 2022/06).

**Funding** Open Access funding provided thanks to the CRUE-CSIC agreement with Springer Nature.

**Data availability** Data can be obtained by reasonable request to the corresponding author.

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

## References

1. Tam CKW (1995) Computational acoustics: issues and methods. AIAA J 33(10):1788-1796

2. Tam CKW (2004) Computational aeroacoustics: an overview of computational challenges and applications. Int J Comput Fluid Dyn 18(6):547-567

3. Colonius T, Lele SK (2004) Computational aeroacoustics: progress on nonlinear problems of sound generation. Prog Aerosp Sci 40:345-416

4. Lele SK, Nichols JW (2014) A second golden age of aeroacoustics? Philos Trans R Soc A 372:20130321

5. Tam CKW, Webb JC (1993) Dispersion-relation-preserving finite difference schemes for computational acoustics. J Comput Phys 107:262-281

6. Lele SK (1992) Compact finite difference schemes with spectral-like resolution. J Comput Phys 103:16-42

7. Bogey C, Bailly C (2004) A family of low dispersive and low dissipative explicit schemes for flow and noise computations. J Comput Phys 194:194-214

8. Daude F, Berland J, Emmert T, Lafon P, Crouzet F, Bailly C (2012) A high-order finite-difference algorithm for direct computation of aerodynamic sound. Comput Fluids 61:46-63

9. Nogueira X, Khelladi S, Colominas I, Cueto-Felgueroso L, París J, Gómez H (2011) High-resolution finite volume methods on unstructured grids for turbulence and aeroacoustics. Arch Comput Methods Eng 18:315-340

10. Khelladi S, Nogueira X, Bakir F, Colominas I (2011) Toward a Higher-Order Unsteady Finite Volume Solver Based on Reproducing Kernel Particle Method. Comput Methods Appl Mech Eng 200:2348-2362

11. Foulquié C, Khelladi S, Deligant M, Ramírez L, Nogueira X, Mardjono J (2020) Numerical assessment of fan blades screen effect on fan/OGV interaction tonal noise. J Sound Vib 481:115428
12. Ramírez L, Nogueira X, Khelladi S, Krimi A, Colominas I (2018) A very accurate arbitrary Lagrangian–Eulerian mesh-less method for computational aeroacoustics. Comput Methods Appl Mech Eng 342:116–141

13. Cueto-Felgueroso L, Colominas I, Nogueira X, Navarrina F, Casteleiro M (2007) Finite volume solvers and moving least-squares approximations for the compressible Navier–Stokes equations on unstructured grids. Comput Methods Appl Mech Eng 196:4712–4736

14. Ollivier-Gooch CF, Van Altena M (2002) A high-order accurate unstructured mesh finite-volume scheme for the advection-diffusion equation. J Comput Phys 181(2):729–752

15. Ollivier-Gooch CF, Nejat A, Michalak K (2007) On obtaining high-order finite-volume solutions to the Euler equations on unstructured meshes, AIAA 2007-4464. In: 18th AIAA computational fluid dynamics conference

16. Toro E, Hidalgo A (2009) ADER finite volume schemes for nonlinear reaction diffusion equations. Appl Numer Math 59(1):73–100

17. Titarev VA, Toro EF (2005) ADER schemes for three-dimensional non-linear hyperbolic systems. J Comput Phys 204:715–736

18. Hu CQ, Shu CW (1999) Weighted essentially non-oscillatory schemes on triangular meshes. J Comput Phys 150(1):97–127

19. Friedrich O (1998) Weighted essentially non-oscillatory schemes for the interpolation of mean values on unstructured grids. J Comput Phys 144(1):194–212

20. Dumbser M, Kaser M (2007) Arbitrary high order non-oscillatory finite volume schemes on unstructured meshes for linear hyperbolic systems. J Comput Phys 221(2):693–723

21. Tsoutsanis P, Antoniadis A, Drikakis D (2014) WENO schemes on arbitrary unstructured meshes for laminar, transitional and turbulent flows. J Comput Phys 256:254–276

22. Tsoutsanis P, Dumbser M (2021) Arbitrary high order central non-oscillatory schemes on mixed-element unstructured meshes. Comput Fluids 225:104961

23. Ramírez L, Edreira L, Couceiro I, Ouro P, Nogueira X, Colominas I (2023) A new mean preserving moving least squares method for arbitrary order finite volume schemes. Appl Math Comput 443:127768

24. Shu CW, Osher S (1998) Efficient implementation of essentially nonoscillatory shock-capturing schemes. J Comput Phys 77:439–471

25. Roe PL (1981) Approximate Riemann solvers, parameter vectors and difference schemes. J Comput Phys 43:357–372

26. Rusanov VV (1962) The calculation of the interaction of non-stationary shock waves and obstacles. USSR Comput Math Math Phys 1(2):304–320

27. Hardin JC, Ristorcelli JR, Tam CKW (1995) ICASE/LaRC workshop on benchmark problems in computational aeroacoustics, NASA Conference Publication 3300

28. Ghosh D, Baeder JD (2012) Compact reconstruction schemes with weighted ENO limiting for hyperbolic conservation laws. SIAM J Sci Comput 34(3):A1678–A1706

29. Moguen Y, Delmas S, Perrier V, Bruel P, Dick E (2015) Godunov-type schemes with an inertia term for unsteady full Mach number range flow calculations. J Comput Phys 281:556–590

30. Moguen Y, Bruel P, Dick E (2019) A combined momentum-interpolation and advection upstream splitting pressure-correction algorithm for simulation of convective and acoustic transport at all levels of Mach number. J Comput Phys 384:16–41

31. Nogueira X, Cueto-Felgueroso L, Colominas I, Khelladi S, Navarrina F, Casteleiro M (2010) Resolution of computational aeroacoustics problem on unstructured grids with a high-order finite volume scheme. J Comput Appl Math 234(7):2089–2097

32. Bailly C, Juvé D (2000) Numerical solution of acoustic propagation problems using linearized Euler equations. AIAA J 38(1):22–29

33. Bogey C, Bailly C, Juvé D (2002) Computation of flow noise using source terms in linearized Euler's equations. AIAA J 40(2):235–243

34. Foulquié C, Mardjono J, Khelladi S, Deligant M, Henner M (2016) High-order aeroacoustics propagation solver with sliding-mesh capabilities for subsonic turbomachinery, AIAA 2016-2723

35. Cand M (2005) A 3D high-order aeroacoustics model for turbomachinery fan noise propagation, PhD thesis, Imperial College London, Department of Mechanical Engineering