Applications of the discontinuous Galerkin method to propagating acoustic wave problems

Abstract

The purpose of this article is to demonstrate that the discontinuous Galerkin method is efficient and suitable to solve linearized Euler equations, modelling sound propagation phenomena. Several benchmark problems were chosen for this purpose. We studied the effect of the underlying computational mesh on the convergence rate and showed the importance of high-quality meshes in order to achieve the theoretical convergence rates. Various acoustic boundary conditions were examined. Perfectly matched layer was used as a non-reflecting boundary condition.

Keywords

Computational aeroacoustics linearized Euler equations discontinuous Galerkin method effect of mesh perfectly matched layer

Introduction

The linearized Euler equations (LEE) are often used as the mathematical model to simulate acoustic propagation. An alternative to the LEE is the acoustic perturbation equation (APE).¹ High-order methods are necessary to achieve the desirable accuracy of acoustic simulations. The discontinuous Galerkin method (DGM) therefore seems ideal for the spatial discretization of the LEE. It allows one to reach high-order accuracy while maintaining the ability to deal with complex geometries. This is achieved by combining features of the finite element method (high-order solution) and finite volume method (local element-based solution). The DGM was first introduced in Cockburn and Shu² and has been since developed for hyperbolic systems. The first successful application of the DGM was presented in Reed and Hill³ to solve the neutron transport equation. Examples of the successful application of DGM to acoustic problems are Reymen et al.⁴ or Bauer et al.⁵ In Reymen et al.,⁴ the DGM was applied to three-dimensional (3D) acoustic pulse propagation problem to analyse the convergence rate of the method. In Bauer et al.,⁵ the DGM was applied to solve two-dimensional (2D) APE-4 for monopole in boundary layer problem. Our aim is to test the DGM for several 2D benchmark problems with various boundary conditions (BCs). Additionally, different meshing algorithms were tested to analyse the effect of the computational mesh on the convergence rate of the method.

Outline

The article is organized as follows. In section ‘LEE’, the 2D LEE are derived from the more general Euler equations (EE). Section ‘Discretization’ describes the nodal DGM, which is used for the spatial discretization of the LEE, as well as the specification of numerical flux. Numerical flux is one of the most important aspects of this method and a common feature of the DGM and the finite volume method. It is the numerical flux that allows us to describe the dynamics of the system under consideration. We then turn our attention to time discretization, where we utilize the explicit Runge–Kutta method to perform time integration. The effect of the underlying mesh used for computation is also discussed. Experimenting with different meshes revealed that the theoretical convergence rate $N + 1$ , where $N$ is the order of the approximation, is not always reached and that this is closely related to the quality of the computational mesh. We therefore experimented with various meshing software and algorithms, and we present our findings before moving on to acoustic applications.

The second part of this article is devoted to applying the DGM to acoustic simulations. We chose a simple acoustic pulse propagation problem as a benchmark problem to demonstrate the high-order nature of this algorithm. We also test the effectiveness of the BC used, namely, reflecting BC and non-reflecting BC (NRBC). The latter is used to simulate the propagation of acoustic waves in a truncated domain to ensure that the outgoing waves do not reflect back to the computational domain. One of the most promising approaches to implement NRBC is the perfectly matched layer (PML), which was first formulated in Hu⁶ for LEE. It was found that this formulation had stability issues, and stable formulation for uniform mean flow was presented in Hu⁷ and extended to non-uniform mean flow in Hu.⁸ This approach was later extended to full EE and Navier–Stokes equations (cf. Hu et al.⁹).

LEE

This section is devoted to the derivation of the 2D LEE. They are obtained by linearizing the EE written in primitive variables (density $ρ$ , velocities $u$ , $v$ and pressure $p$ )

\begin{array}{l} \frac{\partial w}{\partial t} + A_{1} (w) \frac{\partial w}{\partial x} + A_{2} (w) \frac{\partial w}{\partial y} = 0, \\ x = (x, y) \in ℝ^{2}, t > 0 \end{array}

(1)

where

\begin{array}{l} w (x, t) = (\begin{matrix} ρ (x, t) \\ u (x, t) \\ v (x, t) \\ p (x, t) \end{matrix}), A_{1} (w) = (\begin{matrix} u & ρ & 0 & 0 \\ 0 & u & 0 & \frac{1}{ρ} \\ 0 & 0 & u & 0 \\ 0 & κ p & 0 & u \end{matrix}), \\ A_{2} (w) = (\begin{matrix} v & 0 & ρ & 0 \\ 0 & v & 0 & 0 \\ 0 & 0 & v & \frac{1}{ρ} \\ 0 & 0 & κ p & v \end{matrix}) \end{array}

(2)

and $κ$ is the ratio of specific heats (adiabatic index), $κ = 1.4$ in the case of dry air.

We assume that the variables can be decomposed into two parts – the reference state (sometimes called mean value) and fluctuations (or perturbations) – in the following way

w (x, t) = w_{0} (x) + w' (x, t)

(3)

where variables with the subscript ‘ $0$ ’ denote the mean values and those with the superscript ′ denote the fluctuations (cf. Ali et al.¹⁰) Furthermore, we suppose that the reference states are known functions of $x$ or constants, which are independent of time $t$ .

We assume that the magnitudes of the fluctuations are very small compared to those of the mean values, that is

| ρ' | << | ρ_{0} |, | u' | << \sqrt{u_{0}^{2} + v_{0}^{2}}, | v' | << \sqrt{u_{0}^{2} + v_{0}^{2}}, | p' | << | p_{0} |

Therefore, we can approximate the matrices $A_{j} (w_{0} + w')$ as

A_{j} (w_{0} + w') \approx A_{j} (w_{0})

(4)

Inserting approximation (4) into EE (1), we obtain the LEE

\begin{array}{l} \frac{\partial w'}{\partial t} + A_{1} (w_{0}) \frac{\partial w'}{\partial x} + A_{2} (u_{0}) \frac{\partial w'}{\partial y} + h = 0, \\ x \in ℝ^{2}, t > 0 \end{array}

(5)

where $h$ is a vector containing derivatives of the mean flow

h = A_{1} (w_{0}) \frac{\partial w_{0}}{\partial x} + A_{2} (w_{0}) \frac{\partial w_{0}}{\partial y}

(6)

If we consider uniform flow, that is, the mean values are constant in space, the term $h$ in equation (5) vanishes since the spatial partial derivatives in equation (6) become zero vectors.

Discretization

This section deals with discretization techniques used to solve the LEE. The DGM is a suitable method for solving the LEE because the LEE represents conservation laws of mass, momenta and energy, and the DGM is capable of respecting these laws on individual mesh elements with high order. In contrast to the standard dispersion-relation-preserving finite difference schemes for computational acoustics (cf. Tam and Webb¹¹), the DGM is able to handle complex geometries.

Spatial discretization

The algorithm we use is based on the nodal DGM described in Hesthaven and Warburton¹² extended to the case of the LEE. We use the quadrature-free implementation as proposed in Atkins and Shu¹³ to avoid numerical integration. The algorithm is presented for the LEE with $h = 0$ . We have left out the superscript ′ in the notation of the unknown perturbations in this section to simplify the notation.

We suppose that the computational domain $Ω$ can be approximated by the union of non-overlapping elements $D^{k}$

Ω ≃ Ω_{h} = ⋃_{k = 1}^{K} D^{k}

(7)

We consider straight-sided triangles although other options are possible. For example, in the case of a domain with curved boundary, elements of higher order can improve accuracy by describing the boundary more accurately (cf. Krivodonova and Berger¹⁴ and Toulorge and Desmet¹⁵). Within each triangle $D^{k}$ , we assume that the local solution is well approximated by a 2D polynomial of order $N$

w_{h}^{k} (x, t) = \sum_{i = 1}^{N_{p}} w_{h}^{k} (x_{i}^{k}, t) ℓ_{i}^{k} (x)

(8)

where $ℓ_{i}^{k} (x)$ are 2D fundamental Lagrange interpolation polynomials of order $N$ on the set of nodes ${x_{i}^{k}}_{i = 1}^{N_{p}}$ and $N_{p} = 1 / 2 (N + 1) (N + 2)$ . The global solution is approximated by direct sum of local solutions

w (x, t) ≃ w_{h} (x, t) = \oplus_{k = 1}^{K} w_{h}^{k} (x, t)

(9)

Local solutions ${w_{h}^{k}}_{k = 1}^{K}$ are defined by the relation

\int_{D^{k}} (\frac{\partial w_{h}^{k}}{\partial t} + A_{1} \frac{\partial w_{h}^{k}}{\partial x} + A_{2} \frac{\partial w_{h}^{k}}{\partial y}) ℓ_{i}^{k} d x = 0, i = 1, \dots, N_{p}

Applying Gauss’s divergence theorem, we obtain

\begin{array}{l} 1 \int_{D^{k}} \frac{\partial w_{h}^{k}}{\partial t} ℓ_{i}^{k} d x = \int_{D^{k}} (A_{1} \frac{\partial ℓ_{i}^{k}}{\partial x} + A_{2} \frac{\partial ℓ_{i}^{k}}{\partial y}) w_{h}^{k} d x \\ - \int_{\partial D^{k}} (A_{1} n_{1}^{k} + A_{2} n_{2}^{k}) w_{h} ℓ_{i}^{k} d S, i = 1, \dots, N_{p} \end{array}

(10)

where $\partial D^{k}$ denotes the boundary of $D^{k}$ and $n^{k} = (n_{1}^{k}, n_{2}^{k})$ is the unit outward normal to $D^{k}$ . Let us denote by $S (k)$ the set of indices of all triangles adjacent to $D^{k}$ . Furthermore, let $\partial D^{kl} = {\bar{D}}^{k} \cap {\bar{D}}^{l}$ , $l \in S (k)$ , denote the common edge of neighbouring triangles $D^{k}$ and $D^{l}$ , and $n^{kl} = (n_{1}^{kl}, n_{2}^{kl})$ is the unit outward normal on $\partial D^{kl}$ . Then, we approximate

\begin{matrix} \int_{\partial D^{kl}} (A_{1} n_{1}^{k} + A_{2} n_{2}^{k}) w_{h} ℓ_{i}^{k} d S \\ = \sum_{l \in S (k)} \int_{\partial D^{kl}} (A_{1} n_{1}^{kl} + A_{2} n_{2}^{kl}) w_{h} ℓ_{i}^{k} d S ≃ \\ \sum_{l \in S (k)} \int_{\partial D^{kl}} H (w_{h}^{k} |_{\partial D^{kl}}, w_{h}^{l} |_{\partial D^{kl}}, n^{kl}) ℓ_{i}^{k} d S \end{matrix}

(11)

where $H (w_{h}^{k} |_{\partial D^{kl}}, w_{h}^{l} |_{\partial D^{kl}}, n^{kl})$ is the numerical flux from $D^{k}$ to $D^{l}$ through $\partial D^{kl}$ . We chose Lax–Friedrichs numerical flux

H (u, v, n) = \frac{1}{2} (A_{1} n_{1} + A_{2} n_{2}) (u + v) + \frac{1}{2} C (n) (u - v)

(12)

where

C (n) = max_{λ} | λ (A_{1} n_{1} + A_{2} n_{2}) |

and $λ (A_{1} n_{1} + A_{2} n_{2})$ are eigenvalues of $A_{1} n_{1} + A_{2} n_{2}$ . From equations (10), (11) and (12), we obtain the system of ordinary differential equations

\begin{array}{l} M^{k} \frac{d {(W_{h}^{k})}^{T}}{d t} = S_{x}^{k} {(A_{1} W_{h}^{k})}^{T} + S_{y}^{k} {(A_{2} W_{h}^{k})}^{T} \\ - \sum_{l \in S (k)} M^{k l} (H (W_{h}^{k}, W_{h}^{l}, n^{k l}))^{T} \end{array}

(13)

where $k = 1, \dots, K$

\begin{matrix} M^{k} = \int_{D^{k}} (ℓ^{k})^{T} ℓ^{k} d x, S_{x}^{k} = \int_{D^{k}} \frac{\partial (ℓ^{k})^{T}}{\partial x} ℓ^{k} d x, \\ S_{y}^{k} = \int_{D^{k}} \frac{\partial (ℓ^{k})^{T}}{\partial y} ℓ^{k} d x \\ M^{kl} = \int_{\partial D^{kl}} (ℓ^{k})^{T} ℓ^{k} d S \end{matrix}

and

ℓ^{k} = (ℓ_{1}^{k} (x), \dots, ℓ_{N_{p}}^{k} (x)), W_{h}^{k} = (w_{h}^{k} (x_{1}, t), \dots, w_{h}^{k} (x_{N_{p}}, t))

The advantage of straight-sided triangles is that we can map each triangle to a reference triangle by means of linear transformation with a constant Jacobian matrix. This way, the above integrals need to be evaluated only once. We refer the reader to Hesthaven and Warburton¹² for more information on how to compute the matrices defined above efficiently.

Time discretization

As for time discretization, we used the explicit Runge–Kutta method of order 4. We shall focus on low-storage schemes due to memory requirements. Generally speaking, we want to solve an initial value problem for the system of ordinary differential equations of order 1, written in matrix form as

y' = f (t, y), y (t_{0}) = y_{0}

The low-storage explicit Runge–Kutta method reads as

p^{(0)} = y_{n}

for i = 1, \dots, s {\begin{matrix} k^{(i)} = a_{i} k^{(i - 1)} + Δ t f (t_{n} + c_{i} Δ t, p^{(i - 1)}) \\ p^{(i)} = p^{(i - 1)} + b_{i} k^{(i)} \end{matrix}

y_{n + 1} = p^{(s)}

where $s$ is the number of stages of a particular scheme and $Δ t$ denotes timestep. We chose the coefficients $a_{i}$ , $b_{i}$ and $c_{i}$ from Carpenter and Kennedy.¹⁶ In computational aeroacoustics (CAA) applications, it is useful to employ special schemes, called low-dissipation and low-dispersion schemes, which minimize dissipation and dispersion errors. Examples are given in Hu et al.¹⁷ or Berland et al.¹⁸

Effect of the mesh

This section is devoted to studying the effect of the underlying mesh on the convergence rate of the numerical method presented in the previous section. The ideal theoretical convergence rate is $N + 1$ , where $N$ is the order of the polynomial approximating the local solution.

Our first step for applying the DGM to acoustic simulations was to test this method for solving the advection equation

\frac{\partial u}{\partial t} + \frac{\partial u}{\partial x} + \frac{\partial u}{\partial y} = 0

for $x \in [- 5, 15]^{2}$ and $t \in [0, 10]$ . We first used Gmsh,¹⁹ which is a free-meshing and post-processing tool, to create the computational meshes. The results were satisfactory for moderate values of $N$ , but we experienced a drop in the convergence rate for high order of the local solution, namely, $N = 6$ and $N = 8$ . We used two of the meshing algorithms available in Gmsh – Delaunay and Frontal. We also compared these two algorithms with Gambit meshing software included in Fluent simulation software.

When talking about mesh quality, we have to define the shape quality parameter for an element. We used the one from Geuzaine and Remacle,¹⁹ defined as

γ_{k} = \frac{4 \sin α \sin β \sin γ}{\sin α + \sin β + \sin γ}

where $α$ , $β$ and $γ$ are the inner angles of the kth-triangle. This parameter ranges from 0 to 1, where the equilateral triangle is assigned the value 1 and the flat triangle the value 0. We have then computed this quality parameter for all triangles in the resulting mesh. A testing domain $Ω = [- 1, 1]^{2}$ was chosen to study the quality of the mesh produced by the three meshing algorithms. The results are presented in Figures 1 –3. The figures show the generated mesh (on the left side), as well as a histogram of element quality parameter (on the right side). For the histogram plots, x-axis denotes the element quality parameter and y-axis denotes the percentage of elements falling between given quality bounds.

Figure 1.

Gmsh mesh generator (Delaunay algorithm) – mesh (left) and quality parameter histogram (right).

Figure 2.

Gmsh mesh generator (Frontal algorithm) – mesh (left) and quality parameter histogram (right).

Figure 3.

Gambit mesh generator – mesh (left) and quality parameter histogram (right).

Based on our results, we can state that Gambit produces meshes of the highest quality. For a very fine mesh, we can see that Gambit produces mesh with the highest percentage of high-quality elements. We also observed that the Frontal algorithm produces a fairly structured mesh with most of the elements oriented in the same way. In contrast to that, the Delaunay algorithm implemented in Gmsh produces totally unstructured grids while element size varies more significantly.

When studying the convergence rate of the numerical method applied to the advection equation mentioned above, the following convergence rates were experimentally observed for the case $N = 8$ . We summarize our findings in Table 1. We see that only Gambit meshing software can produce meshes that allow us to reach theoretical convergence rates for high orders of local approximation.

Table 1.

Convergence rates for different algorithms.

Algorithm	Delaunay	Frontal	Gambit
Convergence rate	$8.545$	$8.784$	$9.197$

BCs

Next, we describe acoustic BCs. Reflecting BCs are used to reflect acoustic waves from the rigid wall in examples presented in subsections ‘Verification of reflecting BC’ and ‘Duct acoustics’. NRBCs are used to truncate the computational domain for the simulations corresponding to the sound propagation in an unbounded domain.

The reflecting BC is modelled by the slip condition, which is prescribed for the unknown perturbations (cf. Toulorge²⁰). Let us denote $u' = (u', v')^{T}$ as the velocity perturbation. The slip BC is prescribed as

u' \cdot n = 0

(14)

where $n$ denotes the outward normal to the boundary. This is implemented by direct evaluation of the flux

(n_{1} A_{1} + n_{2} A_{2}) w' = (n_{1} A_{1} + n_{2} A_{2}) {w'}_{t}

where

{w'}_{t} = (ρ', u'_{t}, v'_{t}, p')

and

(u'_{t}, v'_{t}) = u' - (u' \cdot n) n

(15)

As for the NRBC, there are many possibilities to implement it. It can be modelled using characteristics (cf. Thomson²¹) or various absorbing zones (cf. Mani²²). An overview can be found in Colonius,²³ Hu²⁴ and Richards et al.²⁵ In general, absorbing zones are designed in such a way that the waves are attenuated, yielding as little reflections as possible. We chose the PML approach and utilized the formulation from Hu.⁷ Denoting

{\bar{w}}^{'} = (\begin{matrix} w' \\ q \end{matrix})

The modified LEE for uniform mean flow aligned with x-axis and with $M$ for the Mach number can be written in matrix form as

\frac{\partial {\bar{w}}^{'}}{\partial t} + {\bar{A}}_{1} \frac{\partial {\bar{w}}^{'}}{\partial x} + {\bar{A}}_{2} \frac{\partial {\bar{w}}^{'}}{\partial y} = \bar{s}

(16)

where

{\bar{A}}_{1} = (\begin{matrix} A_{1} & σ_{y} A_{1} \\ O_{4} & O_{4} \end{matrix}), {\bar{A}}_{2} = (\begin{matrix} A_{2} & σ_{x} A_{2} \\ O_{4} & O_{4} \end{matrix})

and the right-hand side term is given by

\bar{s} = - (\begin{matrix} (σ_{x} + σ_{y}) I_{4} + \frac{σ_{x} M}{1 - M^{2}} A_{1} & σ_{x} σ_{y} (I_{4} + \frac{M}{1 - M^{2}} A_{1}) \\ - I_{4} & O_{4} \end{matrix}) {\bar{w}}^{'}

Here, $O_{4}$ is the zero matrix of order $4$ , $I_{4}$ is the unit matrix of order $4$ , and $q$ is the auxiliary vector defined by

\frac{\partial q}{\partial t} = w'

Coefficients $σ_{x}$ and $σ_{y}$ are often considered to be power functions

σ_{x} (x) = σ_{m} {| \frac{x - x_{0}}{D} |}^{β}, σ_{y} (y) = σ_{m} {| \frac{y - y_{0}}{D} |}^{β}

(17)

where $D$ denotes the thickness of the PML, $β$ is an integer, and $x_{0}$ and $y_{0}$ denote the beginning of PML in the direction of x-axis and y-axis, respectively.

To terminate the PML, we impose the following BCs, which were proposed in Hu.⁷ These are

\begin{array}{l} At outflow x = X_{m a x}, y = Y_{m i n} and y = Y_{m a x} : p = 0 \\ At inflow x = X_{m i n} : ρ = v = p = 0 \end{array}

(18)

where $[X_{\min}, X_{\max}] \times [Y_{\min}, Y_{\max}]$ denotes the whole domain as shown in Figure 4.

Figure 4.

PML domain – absorption coefficients.

Discretization of equation (16) is carried out similarly to the LEE because these two sets of equations are formally similar. Furthermore, it was discovered that eigenvalues of the matrix $n_{1} {\bar{A}}_{1} + n_{2} {\bar{A}}_{1}$ are the same as those of the matrix $n_{1} A_{1} + n_{2} A_{2}$ with the only difference that the additional four eigenvalues are equal to zero.

Acoustic simulations

In this section, we test the DGM on several benchmark problems of acoustics. The benchmark problems presented in this article are taken from the Proceedings of Computational Aeroacoustics Workshop sponsored by the National Aeronautics and Space Administration (NASA) and the Institute for Computer Applications in Science and Engineering (ICASE), published in 1995 as ICASE/LaRC Workshop on Benchmark Problems in Computational Aeroacoustics (cf. Tam²⁶). Our implementation in MATLAB environment was used to solve these test problems (based on DGM implementation presented in Hesthaven and Warburton¹²).

For these test problems, we define the reference speed of sound, $a_{0}$ , and the reference Mach number of the mean flow, $M_{0}$ , as

a_{0} = \sqrt{\frac{γ p_{0}}{ρ_{0}}}, M_{0} = \frac{\sqrt{u_{0}^{2} + v_{0}^{2}}}{a_{0}}

(19)

It is convenient to use length, time, and density scales to further define velocity and pressure scales. We may then use a dimensionless version of LEE using these scales (cf. Tam²⁶). The velocity of the mean flow is then specified by $M_{x}$ and $M_{y}$ , which denote Mach numbers of the mean flow in the direction of x-axis and y-axis, respectively.

Propagation of an acoustic pulse

We begin with an initial value problem designed to simulate the propagation of acoustic waves. We consider a medium at rest, that is, $M_{x} = M_{y} = 0$ . As for the initial conditions (IC), we consider an acoustic Gaussian pulse, which is done by prescribing the following IC at $t = 0$

ρ' = ε \exp (- κ r^{2}), u' = v' = 0, p' = ε \exp (- κ r^{2})

(20)

where $r^{2} = (x - x_{a})^{2} + (y - y_{a})^{2}$ with $x_{a} = y_{a} = 0$ , that is, the pulse is located at the origin. Furthermore, we prescribe the amplitude $ε = 1$ and $κ = \ln 2 / b^{2}$ , where $b = 2$ is the half-width of this pulse.

The analytical solution was derived in Tam and Webb.¹¹ We used it to evaluate the error of the numerical method. We considered a finite domain $Ω = [- 20, 20]^{2}$ and chose computation time $T = 10$ to avoid the acoustic wave reaching the boundary of the domain $Ω$ since our immediate aim is to evaluate the propagation of acoustic waves and not to verify any BC. We only need to prescribe BC to ensure the stability of the scheme, so we chose the reflecting BC for this purpose. We solved this problem using the element size parameter $h = 1$ (resulting in $K = 4272$ elements when using the Delaunay algorithm in Gmsh) and $N = 6$ , the degree of the polynomial representing the local solution within an element. Figure 5 shows the solution for pressure perturbations: on the left side, we present the computed numerical solution, and on the right side, comparison with the analytical solution along y-axis is plotted. Next, we present the numerical solution and comparisons with the analytical solution along y-axis and x-axis, respectively, for horizontal velocity fluctuations $u'$ and vertical velocity perturbations $v'$ in Figures 6 and 7, respectively. We observe a high level of agreement between the numerical solution and the analytical solution in all three cases.

Figure 5.

Solution for pressure perturbations – top view (left) and cut at $y = 0$ (right).

Figure 6.

Solution for velocity perturbations $u^{'}$ – top view (left) and cut at $y = 0$ (right).

Figure 7.

Solution for velocity perturbations $v^{'}$ – top view (left) and cut at $x = 0$ (right).

Seeing that the DGM is able to reproduce propagating acoustic waves very accurately, we would like to verify the convergence rate of this method. To this end, we solved this problem for different combinations of mesh size parameter $h$ (and a corresponding number of elements $K$ ) and order $N$ of approximating polynomial. We used $L^{1}$ -norm to measure the error of our numerical solution (cf. Rinaldi et al.²⁷). The computed $L^{1}$ -norm of errors for pressure perturbations are shown in Table 2. Specifications of the used meshes are given in Table 3.

Table 2.

$L^{1}$ -errors for the acoustic pulse propagation problem.

N	h
N	3.50	3.00	2.50	2.00	1.50	1.00	0.75
2	2.49E–03	2.02E–03	9.84E–04	5.15E–04	2.02E–04	6.02E–05	2.41E–05
4	7.25E–05	4.06E–05	1.73E–05	5.83E–06	1.33E–06	1.98E–07	4.44E–08
6	1.99E–06	7.59E–07	2.74E–07	6.13E–08	7.77E–09	5.18E–10	6.51E–11
8	4.76E–08	1.23E–08	3.54E–09	5.09E–10	3.51E–11	1.32E–11	7.00E–12

Table 3.

Specification of meshes.

h	3.50	3.00	2.50	2.00	1.50	1.00	0.75
K	344	460	618	942	1734	3726	6750

The observed errors are plotted in Figure 8. The scale of both axes is set to be logarithmic and a line is fitted to these values in terms of the least squares method. The obtained convergence rates are presented in Table 4.

Figure 8.

Convergence analysis for acoustic pulse propagation problem.

Table 4.

Convergence rates for acoustic pulse propagation.

Solution order $N$	2	4	6	8
Convergence rate	3.17	5.02	6.96	6.50

We can see that almost optimal convergence rates are obtained in all cases except for $N = 8$ , where accuracy deteriorates when using very fine meshes. If we omit the two finest meshes, we obtain the convergence rate $8.90$ . These results were obtained using meshes generated by the Frontal algorithm in Gmsh. Our experiments revealed that meshes generated by Gambit yielded lower convergence rates compared to the Frontal algorithm in Gmsh for this test case, in contrast to the advection problem discussed in section ‘Effect of the mesh’. The computed convergence rates are shown in Figure 8.

Verification of reflecting BC

This is an initial value boundary problem to test the implemented reflecting BC. We again consider a medium at rest, that is, $M_{x} = M_{y} = 0$ , and assume simple IC in the form of an acoustic Gaussian pulse

ρ' = ε \exp (- κ r^{2}), u' = v' = 0, p' = ε \exp (- κ r^{2})

(21)

where this $r^{2} = x^{2} + (y - 25)^{2}$ , that is, the position of the pulse is $(x_{a}, y_{a}) = (0, 25)$ . Furthermore, we prescribe the amplitude $ε = 1$ and $κ = \ln 2 / b^{2}$ , where $b = 5$ is the half-width of this pulse.

The analytical solution, presented in Tam,²⁶ was used to evaluate the numerical results. This time, our computational domain was $Ω = [- 100, 100] \times [0, 150]$ bounded by a wall at $y = 0$ . Our aim here is to verify the implemented reflecting BC, so we chose $T = 75$ as the final time of the computation so that the acoustic wave does not reach the artificial boundaries at $y = 150$ and $x = \pm 100$ . The computational domain, together with the IC, is depicted in Figure 9.

Figure 9.

Domain and initial conditions.

We solved this problem using an element size parameter $h = 3$ (resulting in $K = 8900$ elements) with $N = 6$ being the degree of the polynomial representing the local solution within an element. We first show the computed solution for density (on the left side) and velocity $u'$ (on the right side) perturbations in Figure 10 and for velocity $v'$ (on the left side) and pressure (on the right side) perturbations in Figure 11. We can see that the solutions for density and pressure perturbations are the same, which is in agreement with the analytical solution.

Figure 10.

Solution for density (left) and velocity $u^{'}$ (right) fluctuations.

Figure 11.

Solution for velocity $v^{'}$ (left) and pressure (right) fluctuations.

To show more clearly that the implemented reflecting BC yields satisfactory results, we present the solution for pressure perturbations along y-axis (on the left side) and pressure perturbations in time at point $(x, y) = (0, 0)$ on the wall (on the right side) in Figure 12. There is a high level of agreement with the analytical solution in both cases.

Figure 12.

Pressure fluctuations – cut at $y = 0$ (left) and solution at point $(x, y) = (0, 0)$ on wall in time (right).

Verification of NRBCs

This is an initial value boundary problem designed to test the implemented NRBCs. The mean flow is assumed to be in the direction of the x-axis with its velocity equal to half of the reference speed of sound, that is, $M_{x} = 0.5$ and $M_{y} = 0$ . The IC is slightly more general in contrast to the acoustic pulse propagation problem because we assume all three types of pulses –acoustic, entropy and vorticity (LEE support three types of waves, cf. Roeck²⁸). Therefore, the IC at $t = 0$ are

\begin{matrix} ρ' = ε_{1} \exp (- κ_{1} r_{a}^{2}) + ε_{2} \exp (- κ_{2} r_{e}^{2}) \\ u' = ε_{3} (y - y_{v}) \exp (- κ_{3} r_{v}^{2}) \\ v' = - ε_{3} (x - x_{v}) \exp (- κ_{3} r_{v}^{2}) \\ p' = ε_{1} \exp (- κ_{1} r_{a}^{2}) \end{matrix}

(22)

where $r_{a}^{2} = (x - x_{a})^{2} + (y - y_{a})^{2}$ , $r_{e}^{2} = (x - x_{e})^{2} + (y - y_{e})^{2}$ , $r_{v}^{2} = (x - x_{v})^{2} + (y - y_{v})^{2}$ and $κ_{j} = \ln 2 / b_{j}^{2}$ , $j = 1, 2, 3$ . We chose $(x_{a}, y_{a}) = (0, 0)$ to specify the position of the acoustic pulse and similarly $(x_{e}, y_{e}) = (67, 0)$ and $(x_{v}, y_{v}) = (67, 0)$ are positions of the entropy and the vorticity pulses, respectively. The amplitudes and half-widths of these pulses are summarized in Table 5.

Table 5.

Parameters $ε_{j}$ and $b_{j}$ in equation (22).

$j$	$ε_{j}$	$b_{j}$
1	1.00	3
2	0.10	5
3	0.04	5

We compared our results with the analytical solution provided in Tam.²⁶ We assumed a computational domain $Ω = [- 100, 100]^{2}$ , where the whole boundary was treated with NRBC. The positions of the initial pulses, together with the computational domain, are shown in Figure 13.

Figure 13.

Domain and initial conditions.

We chose $N = 6$ and $h = 3$ , resulting in $K = 11, 750$ elements in the domain of interest and 2292 elements in the absorbing zones. In this case, NRBC based on PML was implemented. We set $D = 10$ as the width of the PML zones, $σ_{m} = 2$ as the magnitude of the absorption coefficients and $β = 2$ , that is, quadratic profiles of the absorption coefficients. We let the computation run until $T = 300$ to test whether all the waves are correctly attenuated and no waves reflect back to the computational domain of interest. To verify this, we stored solutions at three points $(x_{R}, y_{R}) = (100, 0)$ , $(x_{T}, y_{T}) = (0, 100)$ and $(x_{L}, y_{L}) = (- 100, 0)$ during the computation and then compared them to the provided analytical solution in the very same points in time interval $[0, 300]$ . A high level of agreement was observed in all three points. We present results for the point located on the right side of our computational domain in Figures 14 and 15.

Figure 14.

Density (left) and velocity $u^{'}$ (right) fluctuations in time at point $(x, y) = (100, 0)$ on boundary.

Figure 15.

Velocity $v^{'}$ (left) and pressure (right) fluctuations in time at point $(x, y) = (100, 0)$ on boundary.

Duct acoustics

Duct acoustics is an important part of acoustic widely used in industry. Flow is often in one direction and the computational domain is bounded by walls in the other directions. The problem we now present does not have a known analytical solution. To verify the implemented NRBCs, we solved the same problem using a significantly larger computational domain terminated by reflecting BC on both ends. The secondary domain was chosen large enough so that the reflected waves are not able to reach our domain of interest. The solution obtained by simulation using the larger domain, truncated to the original computational domain, can then be used as a reference solution to the one obtained by simulation using the original domain with NRBCs.

Here, we attempt to solve a simple testing problem presented in Hu.²⁴ We assumed a uniform mean flow with $M_{x} = 0.5$ and $M_{y} = 0$ and the IC in the form of an acoustic Gaussian pulse

ρ' = ε \exp (- κ r^{2}), u' = v' = 0, p' = ε \exp (- κ r^{2})

(23)

where $r^{2} = (x + 50)^{2} + y^{2}$ , the amplitude of the pulse is $ε = 1$ and $κ = \ln 2 / b^{2}$ , where $b = 2$ is the half-width of this pulse.

Our computational domain was a rectangular domain $Ω = [- 100, 100] \times [- 50, 50]$ bounded by walls at $y = \pm 50$ . The boundary at $x = \pm 100$ was modelled with NRBC. The computational domain is shown in Figure 16.

Figure 16.

Domain and initial condition.

We solved this problem using the element size parameter $h = 3$ (resulting in $K = 5984$ elements in the actual domain of interest and 920 elements in the absorbing zones) with $N = 4$ being the degree of the polynomial representing the local solution within an element. We run the simulation until $T = 300$ to test the efficiency of the implemented NRBC. PML was used again as the NRBC with parameters set to $D = 15$ , $β = 2$ and $σ_{m} = 0.5$ . The black vertical lines represent where PML zones start. We present a pressure field at times $t = 30$ , $t = 100$ , $t = 180$ and $t = 270$ in Figures 17 –20, respectively.

Figure 17.

Pressure fluctuations at $t = 30$ .

Figure 18.

Pressure fluctuations at $t = 100$ .

Figure 19.

Pressure fluctuations at $t = 180$ .

Figure 20.

Pressure fluctuations at $t = 270$ .

Figure 21 shows the comparison with the solution obtained using domain $Ω = [- 325, 375] \times [- 50, 50]$ with an order of approximation $N = 6$ , which serves as a reference solution to validate the implemented NRBC. We compare the solution at point $(x, y) = (100, 0)$ , which lies on the right boundary of the truncated domain, in time interval $[0, 300]$ . The results agree closely with the reference solution.

Figure 21.

Comparison for pressure perturbations with reference solution in time at point $(x, y) = (100, 0)$ on boundary.

We see that all outgoing waves are correctly attenuated in the PML zones (compared to the reference solution obtained by the previous simulation on larger computational domain). These results are in agreement with those presented in Hu²⁴ and Nogueira et al.²⁹

Conclusion

We presented the nodal DGM as a suitable method for acoustic simulations. We performed simulations of acoustic pulse propagation to demonstrate that the DGM can reach a high convergence rate. We also tested typical acoustic BC, namely, reflecting BC and NRBC, and we verified that the DGM enables their effective implementation. In particular, we have also shown how to use the PML approach in the DGM. In addition, we demonstrated the importance of high-quality computational mesh. Our experiments revealed that the actual convergence rate of the method is strongly affected by quality of the underlying mesh.

All the experiments were done using MATLAB code. The algorithm was then ported to C++ programming language, which yielded approximately 3× shorter execution times. Additionally, work has begun on implementation of distributed simulation for computational cluster. We would like to extend the algorithm to 3D.

Footnotes

Academic Editor: Yi Wang

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship and/or publication of this article.

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship and/or publication of this article: This study was financially supported by the project no. GA 13-27505S, funded by the Czech Grant Agency.

References

Ewert

Schröder

. Acoustic perturbation equations based on flow decomposition via source filtering. J Comput Phys 2003; 188: 365–398.

Cockburn

Shu

C-W

. TVB Runge-Kutta local projection discontinuous Galerkin finite element method for conservation laws II: general framework. Math Comput 1989; 52: 411–435.

Reed

Hill

. Triangular mesh methods for the neutron transport equation. Technical Report LA-UR-73-479, Los Alamos Scientific Laboratory, Los Alamos, NM, 1973.

Reymen

Roeck

Rubio

. A 3D discontinuous Galerkin method for aeroacoustic propagation. In: Proceedings of the 12th international congress on sound and vibration, 11–14 July 2005, pp.761–770. Lisbon, Portugal: International Institute of Acoustics and Vibration.

Bauer

Dierke

Ewert

. Application of a discontinuous Galerkin method to discretize acoustic perturbation equations. AIAA J 2011; 49: 898–908.

. On absorbing boundary conditions for linearized Euler equations by a perfectly matched layer. J Comput Phys 1996; 129: 201–219.

. A stable perfectly matched layer for linearized Euler equations in unsplit physical variables. J Comput Phys 2001; 175: 455–480.

. A perfectly matched layer absorbing boundary condition for linearized Euler equations with a non-uniform mean flow. J Comput Phys 2005; 208: 469–492.

Lin

. Absorbing boundary conditions for nonlinear Euler and Navier-stokes equations based on the perfectly matched layer technique. J Comput Phys 2008; 227: 4398–4424.

10.

Ali

Becker

Utzmann

. Aeroacoustic study of a forward facing step using linearized Euler equations. Physica D 2008; 237: 2184–2189.

11.

Tam

Webb

. Dispersion-relation-preserving finite difference schemes for computational acoustics. J Comput Phys 1993; 107: 262–281.

12.

Hesthaven

Warburton

. Nodal discontinuous Galerkin methods: algorithms, analysis, and applications. New York: Springer, 2008.

13.

Atkins

Shu

C-W

. Quadrature-free implementation of discontinuous Galerkin method for hyperbolic equations. AIAA J 1998; 36: 775–782.

14.

Krivodonova

Berger

. High-order accurate implementation of solid wall boundary conditions in curved geometries. J Comput Phys 2006; 211: 492–512.

15.

Toulorge

Desmet

. Curved boundary treatments for the discontinuous Galerkin method applied to aeroacoustic propagation. AIAA J 2010; 48: 479–489.

16.

Carpenter

Kennedy

. Fourth-order 2N-storage Runge–Kutta schemes. Technical Report, NASA Langley Research Center, Hampton, June 1994.

17.

Hussaini

Manthey

. Low-dissipation and low-dispersion Runge-Kutta schemes for computational aeroacoustics. J Comput Phys 1996; 124: 177–191.

18.

Berland

Bogey

Bailly

. Low-dissipation and low-dispersion fourth-order Runge–Kutta algorithm. Comput Fluids 2005; 35: 1459–1463.

19.

Geuzaine

Remacle

. Gmsh: a three-dimensional finite element mesh generator with built-in pre- and post-processing facilities. Int J Numer Meth Eng 2009; 79: 1309–1331.

20.

Toulorge

. Efficient Runge-Kutta discontinuous Galerkin methods applied to aeroacoustics. PhD Thesis, Department of Mechanical Engineering, Faculty of Engineering, Katholieke Universiteit Leuven, 2012.

21.

Thomson

. Time dependent boundary conditions for hyperbolic systems. J Comput Phys 1987; 68: 1–24.

22.

Mani

. Analysis and optimization of numerical sponge layers as a nonreflective boundary treatment. J Comput Phys 2012; 231: 704–716.

23.

Colonius

. Modeling artificial boundary conditions for compressible flow. Annu Rev Fluid Mech 2004; 36: 315–345.

24.

. Absorbing boundary conditions. Int J Comput Fluid Dynamics 2004; 18: 513–522.

25.

Richards

Zhang

Chen

. The evaluation of non-reflecting boundary conditions for duct acoustic computation. J Sound Vib 2004; 270: 539–557.

26.

Tam

. Benchmark problems and solutions. In: Proceedings of the ICASE/LaRC workshop on benchmark problems in computational aeroacoustics (CAA): NASA conference publication, Hampton, 24–26 October 1994, pp.1–13. Hampton: National Aeronautics and Space Administration, Langley Research Center.

27.

Rinaldi

RDR

Iob

Arina

. An efficient discontinuous Galerkin method for aeroacoustic propagation. Int J Numer Meth Fl 2012; 69: 1473–1495.

28.

Roeck

. Hybrid methodologies for the computational aeroacoustic analysis of confined, subsonic flows. PhD Thesis, Department of Mechanical Engineering, Faculty of Engineering, Katholieke Universiteit Leuven, 2007.

29.

Nogueira

Colominas

Cueto-Felgueroso

. Resolution of computational aeroacoustics problems on unstructured grids with a higher-order finite volume scheme. J Comput Appl Math 2010; 234: 2089–2097.