Validation of a Discontinuous Galerkin Implementation of the Time-Domain Linearized Navier–Stokes Equations for Aeroacoustics

Abstract

The propagation of small perturbations in complex geometries can involve hydrodynamic-acoustic interactions, coupling acoustic waves and vortical modes. A propagation model, based on the linearized Navier–Stokes equations, is proposed. It includes the mechanism responsible for the generation of vorticity associated with the hydrodynamic modes. The linearized Navier–Stokes equations are discretized in space using a discontinuous Galerkin formulation for unstructured grids. Explicit time integration and non-reflecting boundary conditions are described. The linearized Navier–Stokes (LNS) model is applied to two test cases. The first one is the time-harmonic source line in an incompressible inviscid two-dimensional mean shear flow in an infinite domain. It is shown that the proposed model is able to capture the trailing vorticity field developing behind the mass source and to represent the redistribution of the vorticity. The second test case deals with the analysis of the acoustic propagation of an incoming perturbation inside a circular duct with a sudden area expansion in the presence of a mean flow and the evaluation of its scattering matrix. The computed coefficients of the scattering matrix are compared to experimental data for three different Mach numbers of the mean flow, M${}_0$ = 0.08, 0.19 and 0.29. The good agreement with the experimental data shows that the proposed method is suitable for characterizing the acoustic behavior of this kind of network.

1. Introduction

When developing a model for aeroacoustic problems, it is necessary to establish the level of physical approximation ranging from the wave equation, without or with convective effects, to the linearized Euler or Navier–Stokes equations [1]. For problems in complex geometries, significant hydrodynamic-acoustic interactions, coupling acoustic waves and vortical modes, may occur. For example, in the case of acoustic propagation in ducts with sudden changes of area, where flow separation occurs in correspondence of sharp edges, there is a consequent generation of vorticity due to viscous effects. To correctly capture this coupling, the mechanisms responsible for the generation of vorticity associated with the hydrodynamic modes must be included in the linearized model. Presently, in computational aeroacoustics (CAA), the main focus is on the solution of the linearized Euler equations (LEE), which model both the acoustic propagation, as well as the vorticity transport. The role of viscosity, as a source of vorticity, is introduced by means of some variants of the Kutta condition. However, some ambiguities remain on the vorticity generation process, which is ultimately due to the viscous effects. To have a more exact, and less empirical, model, it is necessary to resort to the linearized Navier–Stokes equations (LNS).

In addition, for the complete definition of the model, it is necessary to define the numerical scheme. For an accurate simulation of the weak acoustic perturbations over long distances, the numerical scheme must have a low degree of numerical dissipation and dispersion. Such properties are typical of high-order numerical schemes. A wide range of high-order numerical schemes have been proposed, among them, the DRP [2] and compact finite difference schemes [3] are the most popular in CAA for structured grids. The occurrence of geometrical complexities makes it necessary to employ unstructured grids and to resort to high-order finite element/volume formulations, such as the spectral element method [4] or the discontinuous Galerkin method (DGM). The first DGM was proposed in 1973 by Reed and Hill [5], but only after 30 years was rediscovered and applied to compressible flows by Bassi and Rebay [6] and Cockburn and co-workers [7]. Many variations based on the discontinuous discretization introduced by the DGM have been proposed since then. Despite all its advantages, the DGM did not have a significant impact, up to now, in practical applications of computational fluid dynamics, because of its high computational cost. The presence of extra nodes with respect to a standard continuous Galerkin method introduces severe stability limitations. However, in CAA, the need for a numerical scheme with low dispersive and dissipative properties in unstructured grids makes DGM more appealing despite its computational cost, and today, it is one of the most popular schemes for the numerical simulation of wave propagation.

The objective of the present work is to include the viscous effects into the acoustic wave propagation model obtained linearizing the Navier–Stokes equations, for a compressible flow, with respect to a representative mean flow. The inclusion of the viscous terms not only solves the problem of the correct generation of vorticity, but may also contribute to solving the well-known problem of representing the Kelvin–Helmholtz instabilities, connected with the propagation of the vortical modes, which make the current LEE models highly unstable. The explicit inclusion of the viscous terms removes the necessity of adding artificial viscosity or to resort to other fictitious mechanisms to stabilize the numerical solution of the LEE. The LNS equations can be a valid alternative to the LEE model for acoustic propagation problems. In this work, an efficient numerical algorithm for the solution of the LNS equations is proposed. To our knowledge, there are only a few works dealing with the solution of the LNS for aeroacoustics. In [8,9], the LNS equations are solved in the frequency domain for acoustic propagation in internal flows at a low Mach number. A method for the solution of the LNS in the time domain is described in [10,11] for the study of a slit resonator under grazing flow. The present method solves the LNS equations in the time domain on unstructured grids. The present work is an extension of previous work on the DGM applied to the solution of the LEE in the frequency and time domains [12].

The paper is organized as follows. In Section 2, the LNS equations are presented for general coordinates and for axisymmetric problems, with the assumption of no tangential velocity and the remaining quantities independent of the tangential coordinate. The DGM formulation for the LNS equations is described in Section 3, as well as the time integration algorithm, based on a low dissipation formulation of a fourth-order accurate Runge–Kutta scheme [13]. Explicit time integration, more appropriate for acoustic wave propagation, avoids the inversion of a large algebraic system, and it is well suited for parallel computation. Along numerical boundaries, to avoid incoming spurious reflections, a sponge-layer boundary condition is used [14]. As the first test case, in Section 4, the harmonic perturbation from a line source in an incompressible inviscid two-dimensional mean flow with linear shear is studied. The LNS equations, as well as the LEE, obtained from the LNS equations by simply dropping the viscous terms, are able to correctly represent the hydrodynamic wake caused by the shed vorticity of the pulsating source and compared to the analytical results of Rienstra et al. [15]. As a second test, in Section 5, the evaluation of the scattering matrix for a sudden area discontinuity in a cylindrical duct in the presence of mean flow is reported. The Mach number ranges from $0.08$ to $0.29$. The numerical results are validated using the experimental data of Ronneberger [16]. Finally, some concluding remarks are reported in Section 6.

2. Linearized Navier–Stokes Equations

2.1. General Coordinate Formulation

The propagation of sound waves in a medium with non-uniform mean velocity is governed by the LNS equations [1]. Denoting with ${(·)}_0$ the mean flow variables and with ${(·)}^{\prime }$ the acoustic perturbations, the governing equations for a general coordinate system read:

$${\displaystyle \frac{\partial \mathbf{U}}{\partial t}}+\nabla ·{\mathbf{F}}^c=\nabla ·{\mathbf{F}}^d+\mathbf{G}+\mathbf{H}$$

(1)

where $\mathbf{U}={\left({\rho }^{\prime },{\rho }_0{\mathbf{u}}^{\prime },p^{\prime }\right)}^T$ is the vector of the acoustic variables. ${\mathbf{F}}^c$ and ${\mathbf{F}}^d$ are the convective and viscous fluxes, defined as follows, $\mathbf{I}$ being the identity matrix,

$${\mathbf{F}}^c=\left(\begin{array}{c}{\rho }_0{\mathbf{u}}^{\prime }+{\rho }^{\prime }{\mathbf{u}}_0\\ {\rho }_0{\mathbf{u}}_0{\mathbf{u}}^{\prime }+p^{\prime }\mathbf{I}\\ \gamma p_0{\mathbf{u}}^{\prime }+p^{\prime }{\mathbf{u}}_0\end{array}\right)$$

(2)

$${\mathbf{F}}^d=\left(\begin{array}{c}0\\ {\tau }^{\prime }\\ (\gamma -1)\kappa \nabla ·T^{\prime }\end{array}\right)$$

(3)

γ is the specific heat ratio for a perfect gas and κ the thermal conductivity coefficient. The viscous stresses ${\tau }^{\prime }$, μ being the dynamic viscosity coefficient and ${\mu }_B$ the bulk viscosity coefficient, are linearly related to the velocity fluctuation gradients:

$${\tau }^{\prime }=\mu \left[(\nabla {\mathbf{u}}^{\prime }+\nabla {\mathbf{u}}^{\prime T})-{\displaystyle \frac{2}{3}}(\nabla ·{\mathbf{u}}^{\prime })\mathbf{I}\right]+{\mu }_B(\nabla ·{\mathbf{u}}^{\prime })\mathbf{I}$$

In the present applications, the bulk viscosity term is neglected, because it affects the acoustic propagation only over extremely long distances, when the accumulation of its influence over each cycle becomes important and eventually dissipates the acoustic perturbation.

The vector $\mathbf{G}$, containing the zeroth order terms and the mean flow derivatives, reads:

$$\mathbf{G}=\left(\begin{array}{c}0\\ -{\rho }^{\prime }\left({\displaystyle \frac{\partial {\mathbf{u}}_0}{\partial t}}+{\mathbf{u}}_0·\nabla {\mathbf{u}}_0\right)-{\rho }_0{\mathbf{u}}^{\prime }·\nabla {\mathbf{u}}_0\\ (\gamma -1)\left({\mathbf{u}}^{\prime }·\nabla p_0-p^{\prime }\nabla ·{\mathbf{u}}_0\right)+(\gamma -1)\left[\left({\tau }_0·\nabla \right)·{\mathbf{u}}^{\prime }+\left({\tau }^{\prime }·\nabla \right)·{\mathbf{u}}_0\right]\end{array}\right)$$

The source vector $\mathbf{H}$ contains explicit source terms, like a mass source or a force source. The linearized state equation for a perfect gas reads:

$${\displaystyle \frac{p^{\prime }}{p_0}}={\displaystyle \frac{{\rho }^{\prime }}{{\rho }_0}}+{\displaystyle \frac{T^{\prime }}{T_0}}$$

The LNS equations can be simplified assuming an isoentropic relation between pressure and density fluctuations. This is not a general assumption, for example not being valid for problems involving combustion, but it is applicable to the problems studied in the present work [8]. In this case, pressure and density are related by the relation $p^{\prime }=c_0^2{\rho }^{\prime }$, where $c_0^2=\gamma p_0/{\rho }_0$ is the velocity of sound. With this assumption, one dependent variable is removed from the system, and therefore, the energy equation can be omitted from the system Equation (1).

The LEE equations can be obtained from Equation (1) by simply dropping out the viscous fluxes ${\mathbf{F}}^d$ and setting to zero ${\tau }^{\prime }$ and ${\tau }_0$ in the expression of $\mathbf{G}$. LEE support hydrodynamic modes; therefore, it is possible to further simplify the governing equations, removing from the LEE the terms responsible for the transport of the non-acoustic modes of vorticity and entropy. The resulting equations are called acoustic perturbation equations (APE) [17]. It must be pointed out that the wave operator in the left-hand side of the APE is exact only in the case of irrotational mean flow fields; consequently, the presence of mean vorticity may cause errors in the computed sound propagation. These errors are assumed to be small for low mean vorticity levels. There are several formulations for the APE; in Section 4, some results are obtained with the APE-4 formulation [12,17].

2.2. Axisymmetric Formulation

For axisymmetric problems, it is convenient to rewrite the governing equations in a cylindrical coordinate system, namely $\left(r,\theta ,z\right)$. If both the geometry and the mean flow can be assumed axisymmetric, i.e.,

$$v_{0\theta }=0,{\displaystyle \frac{\partial }{\partial \theta }}\left[{\left(·\right)}_0\right]=0,{\displaystyle \frac{\partial }{\partial \theta }}\left[{\left(·\right)}^{\prime }\right]=0$$

where $v_{0\theta }$ is the θ-component of the mean flow velocity, the LNS Equation (1) can be written in cylindrical coordinates as:

$${\displaystyle \frac{\partial \mathbf{U}}{\partial t}}+{\displaystyle \frac{1}{r}}{\displaystyle \frac{\partial {\mathbf{F}}_r^c}{\partial r}}+{\displaystyle \frac{\partial {\mathbf{F}}_z^c}{\partial z}}={\displaystyle \frac{1}{r}}{\displaystyle \frac{\partial {\mathbf{F}}_r^d}{\partial r}}+{\displaystyle \frac{\partial {\mathbf{F}}_z^d}{\partial z}}+{\mathbf{F}}_p+{\mathbf{G}}_{axy}+\mathbf{H}$$

(4)

where $\mathbf{U}={\left({\rho }^{\prime },{\rho }_0{v}_r^{\prime },{\rho }_0{v}_z^{\prime },p^{\prime }\right)}^T$ is the acoustic perturbation vector. Here, $(v_r^{\prime },v_z^{\prime })$ are the velocity components in $(r,z)$ directions, respectively. ${\mathbf{F}}_r^c$ and ${\mathbf{F}}_z^c$ contain part of the inviscid fluxes, the pressure gradient term being included in the term ${\mathbf{F}}_p$; ${\mathbf{F}}_r^d$ and ${\mathbf{F}}_z^d$ are the viscous fluxes along r and z directions, respectively:

$${\mathbf{F}}_r^c=\left(\begin{array}{c}{\rho }^{\prime }{v}_{0r}+{\rho }_0{v}_r^{\prime }\\ {\rho }_0{v}_{0r}{v}_r^{\prime }\\ {\rho }_0{v}_{0r}{v}_z^{\prime }\\ r\left(\gamma p_0{v}_r^{\prime }+p^{\prime }{v}_{0r}\right)\end{array}\right),{\mathbf{F}}_r^d=\left(\begin{array}{c}0\\ r{\tau }_{rr}^{\prime }\\ r{\tau }_{rz}^{\prime }\\ (\gamma -1)\kappa r{\displaystyle \frac{\partial T^{\prime }}{\partial r}}\end{array}\right)$$

$${\mathbf{F}}_z^c=\left(\begin{array}{c}{\rho }^{\prime }{v}_{0z}+{\rho }_0{v}_z^{\prime }\\ {\rho }_0{v}_{0z}{v}_r^{\prime }\\ {\rho }_0{v}_{0z}{v}_z^{\prime }\\ \gamma p_0{v}_z^{\prime }+p^{\prime }{v}_{0z}\end{array}\right),{\mathbf{F}}_z^d=\left(\begin{array}{c}0\\ {\tau }_{rz}^{\prime }\\ {\tau }_{zz}^{\prime }\\ (\gamma -1)\kappa {\displaystyle \frac{\partial T^{\prime }}{\partial z}}\end{array}\right)$$

and:

$${\mathbf{F}}_p=-{\left(0,{\displaystyle \frac{\partial p^{\prime }}{\partial r}},{\displaystyle \frac{\partial p^{\prime }}{\partial z}},0\right)}^T$$

The viscous stress components read:

$$\begin{array}{ccc}\hfill {\tau }_{rr}^{\prime }& =& 2\mu {\displaystyle \frac{\partial v_r^{\prime }}{\partial r}}+\left({\mu }_B-{\displaystyle \frac{2}{3}}\mu \right)\left[{\displaystyle \frac{1}{r}}{\displaystyle \frac{\partial }{\partial r}}\left(r{v}_r^{\prime }\right)+{\displaystyle \frac{\partial v^{\prime }}{\partial z}}\right]\hfill \\ \hfill {\tau }_{rz}^{\prime }& =& \mu \left({\displaystyle \frac{\partial v_r^{\prime }}{\partial z}}+{\displaystyle \frac{\partial v_z^{\prime }}{\partial r}}\right)\hfill \\ \hfill {\tau }_{zz}^{\prime }& =& 2\mu {\displaystyle \frac{\partial v_z^{\prime }}{\partial z}}+\left({\mu }_B-{\displaystyle \frac{2}{3}}\mu \right)\left[{\displaystyle \frac{1}{r}}{\displaystyle \frac{\partial }{\partial r}}\left(r{v}_r^{\prime }\right)+{\displaystyle \frac{\partial v^{\prime }}{\partial z}}\right]\hfill \end{array}$$

Vector $\mathbf{H}$ represents the acoustic sources, and ${\mathbf{G}}_{axy}$ contains terms of the mean flow due to axisymmetry and the mean flow derivatives:

$${\mathbf{g}}_{axy}=\left(\begin{array}{c}0\\ -\left[{\rho }^{\prime }{\displaystyle \frac{\partial v_{0r}}{\partial t}}+\left({\rho }^{\prime }{v}_{0r}+{\rho }_0{v}_r^{\prime }\right){\displaystyle \frac{\partial v_{0r}}{\partial r}}+\left({\rho }^{\prime }{v}_{0z}+{\rho }_0{v}_z^{\prime }\right){\displaystyle \frac{\partial v_{0r}}{\partial z}}\right]\\ -\left[{\rho }^{\prime }{\displaystyle \frac{\partial v_{0z}}{\partial t}}+\left({\rho }^{\prime }{v}_{0r}+{\rho }_0{v}_r^{\prime }\right){\displaystyle \frac{\partial v_{0z}}{\partial r}}+\left({\rho }^{\prime }{v}_{0z}+{\rho }_0{v}_z^{\prime }\right){\displaystyle \frac{\partial v_{0z}}{\partial z}}\right]\\ (\gamma -1)g_{axy}^z\end{array}\right)$$

where:

$$\begin{array}{ccc}\hfill g_{axy}^z& =& \left(v_{0r}{\displaystyle \frac{\partial p^{\prime }}{\partial r}}+v_{0z}{\displaystyle \frac{\partial p^{\prime }}{\partial z}}\right)+\left(v_r^{\prime }{\displaystyle \frac{\partial p_0}{\partial r}}+v_z^{\prime }{\displaystyle \frac{\partial p_0}{\partial z}}\right)+\hfill \\ & & \left({\displaystyle \frac{{\tau }_{0rr}}{r}}{\displaystyle \frac{\partial }{\partial r}}\left(r{v}_r^{\prime }\right)+{\tau }_{0rz}{\displaystyle \frac{\partial v_r^{\prime }}{\partial z}}\right)+\left({\displaystyle \frac{{\tau }_{0rz}}{r}}{\displaystyle \frac{\partial }{\partial r}}\left(r{v}_z^{\prime }\right)+{\tau }_{0zz}{\displaystyle \frac{\partial v_z^{\prime }}{\partial z}}\right)+\hfill \\ & & \left({\displaystyle \frac{{\tau }_{rr}^{\prime }}{r}}{\displaystyle \frac{\partial }{\partial r}}\left(r{v}_{0r}\right)+{\tau }_{rz}^{\prime }{\displaystyle \frac{\partial v_{0r}}{\partial z}}\right)+\left({\displaystyle \frac{{\tau }_{rz}^{\prime }}{r}}{\displaystyle \frac{\partial }{\partial r}}\left(r{v}_{0z}\right)+{\tau }_{zz}^{\prime }{\displaystyle \frac{\partial v_{0z}}{\partial z}}\right)\hfill \end{array}$$

3. Discontinuous Galerkin Method

The LNS Equation (1) can be written in a more compact form as:

$${\displaystyle \frac{\partial \mathbf{U}}{\partial t}}+\nabla ·\mathbf{F}=\mathbf{S}$$

(5)

where the vector $\mathbf{F}$ is defined as the sum of the convective and viscous contributions:

$$\mathbf{F}={\mathbf{F}}^c+{\mathbf{F}}^d$$

3.1. Spatial Discretization

The standard DGM is formulated to approximate first order derivatives of the primitive variables, while the viscous components of the fluxes of the LNS Equation (5) show a dependence up to the second-order derivatives of the conservative variables. Among the several extensions of the DGM to convection-diffusion problems with higher-order derivatives, following Yan and Shu [18], it is possible to introduce an auxiliary variable $\mathbf{Q}=\nabla \mathbf{U}$. The global system of equations for the variables $\mathbf{U}$ and $\mathbf{Q}$ is:

$$\begin{array}{ccc}& & \mathbf{Q}-\nabla \mathbf{U}=0\hfill \end{array}$$

(6)

$$\begin{array}{ccc}& & \frac{\partial \mathbf{U}}{\partial t}+\nabla ·\mathbf{F}(\mathbf{U},\mathbf{Q})=\mathbf{S}\hfill \end{array}$$

(7)

The domain Ω is approximated by a partition of non-overlapping elements ${\Omega }^k$:

$$\Omega \approx {\Omega }_h=\bigcup _{k=1}^M{\Omega }^k$$

and within each element ${\Omega }^k$, it is assumed that the solution is approximated by the nodal expansion [19]:

$${\mathbf{U}}_h^k(\mathbf{x},t)=\sum _{l=0}^N{\mathbf{u}}^k({\mathbf{x}}_l,t)L_l^k\left(\mathbf{x}\right)=\sum _{l=0}^N{\mathbf{u}}_l^k\left(t\right)L_l^k\left(\mathbf{x}\right)$$

(8)

where N is the degree of freedom inside the element, ${\Omega }^k$ and $L_l^k\left(\mathbf{x}\right)$ are the multivariate Lagrange interpolation polynomials, defined by the points ${\mathbf{x}}_l$ in the element ${\Omega }^k$, and ${\mathbf{u}}_l^k\left(t\right)$, the expansion coefficients, are the values of the solution at the nodal points. The nodal expansion is preferred with respect to the modal expansion. Both expansions are characterized by an exponential rate of convergence, but in the nodal expansion, the solution is evaluated at the nodal points, while in the modal expansion, the solution is known at the Gauss quadrature points, and interpolation is necessary for recovering the solution at the quadrature points on the elemental edges. The nodal points are obtained with the method proposed in [20], resulting in fully-unstructured nodal sets with a large degree of symmetry, including the element vertices and with complete one-dimensional polynomials supported by the nodes on each edge.

Applying a Galerkin projection to Equation (7) onto each member of the basis set $L_i^k$, the weak form of the problem can be written for each element ${\Omega }^k$ as:

$${\int }_{{\Omega }^k}{L}_i^k\left({\displaystyle \frac{\partial {\mathbf{U}}_h^k}{\partial t}}+\nabla ·{\mathbf{F}}_h^k\right)d{\Omega }^k={\int }_{{\Omega }^k}{L}_i^k{\mathbf{S}}_h^{k}d{\Omega }^k$$

(9)

where the flux and the source terms are approximated using the basis sets, likewise for the conservative variables Equation (8),

$${\mathbf{S}}_h^k=\sum _{l=1}^N{\mathbf{s}}_l^k{L}_l^k,{\mathbf{F}}_h^k=\sum _{l=1}^N{\mathbf{f}}_l^k{L}_l^k$$

(10)

Integrating by parts Equation (9), one obtains:

$${\int }_{{\Omega }^k}{L}_i^k{\displaystyle \frac{\partial {\mathbf{U}}_h^k}{\partial t}}d{\Omega }^k-{\int }_{{\Omega }^k}\nabla L_i^k·{\mathbf{F}}_h^{k}d{\Omega }^k+{\int }_{{\Sigma }^k}{L}_i^k{\mathbf{F}}_h^k·\mathbf{n}d{\Sigma }^k={\int }_{{\Omega }^k}{L}_i^k{\mathbf{S}}_h^{k}d{\Omega }^k$$

(11)

where ${\Sigma }^k$ indicates the boundary limiting ${\Omega }^k$ and $\mathbf{n}$ is the outward pointing unit normal referred to each edge. This term can be written as:

$${\int }_{{\Sigma }^k}{L}_i^k{\mathbf{F}}_h^k·\mathbf{n}d{\Sigma }^k={\int }_{{\Sigma }^k}\left(L_i^{k+}{\mathbf{F}}_h^{k+}·{\mathbf{n}}^{+}+L_i^{k-}{\mathbf{F}}_h^{k-}·{\mathbf{n}}^{-}\right)d{\Sigma }^k$$

(12)

where “+” and “−” distinguish the values of the discontinuous quantities on either side of an interface. The + sign corresponds to the element ${\Omega }^k$, with the normal direction consistent with the edge orientation. Defining the jump operator $\left[\left[\mathbf{A}\right]\right]={\mathbf{A}}^{+}·{\mathbf{n}}^{+}+{\mathbf{A}}^{-}·{\mathbf{n}}^{-}=\left({\mathbf{A}}^{+}-{\mathbf{A}}^{-}\right)·{\mathbf{n}}^{+}$ and the average operator $\langle \mathbf{A}\rangle =\left({\mathbf{A}}^{+}+{\mathbf{A}}^{-}\right)/2$, Relation Equation (12) can be rewritten as:

$${\int }_{{\Sigma }^k}{L}_i^k{\mathbf{F}}_h^k·\mathbf{n}d{\Sigma }^k={\int }_{{\Sigma }^k}\left(\left[\left[L_i^k\right]\right]\langle {\mathbf{F}}_h^k\rangle +\left[\left[{\mathbf{F}}_h^k\right]\right]\langle L_i^k\rangle \right)d{\Sigma }^k$$

(13)

In the DGM, the argument of Integral Equation (13) is replaced by an interface flux function Γ depending on the test function, the local normal and the values of the vector $\mathbf{U}$ and its derivatives of both elemental edges delimiting the interface: $\Gamma \left({\mathbf{U}}^{+},{\mathbf{U}}^{-},L_i^{k+},L_i^{k-},\mathbf{n}\right)$. To ensure consistency, Γ must satisfy the condition, with h a measure of the grid element size,

$$\underset{h\to 0}{lim}\Gamma \left({\mathbf{U}}^{+},{\mathbf{U}}^{-},L_i^{k+},L_i^{k-},\mathbf{n}\right)=\left[\left[{\mathbf{L}}_{\mathbf{i}}^{\mathbf{k}}\right]\right]\langle {\mathbf{F}}_h^k\rangle$$

To make the DGM element-wise conservative, it must be:

$$\Gamma \left({\mathbf{U}}^{+},{\mathbf{U}}^{-},L_i^{k+},0,\mathbf{n}\right)=-\Gamma \left({\mathbf{U}}^{+},{\mathbf{U}}^{-},0,L_i^{k-},\mathbf{n}\right)$$

The convective and diffusive fluxes are discretized separately. The convective discretization is stabilized by using an approximate Riemann solver $\mathcal{H}$ for Γ:

$$\Gamma \left({\mathbf{U}}^{+},{\mathbf{U}}^{-},L_i^{k+},0,\mathbf{n}\right)=L_i^{k+}\mathcal{H}\left({\mathbf{U}}^{+},{\mathbf{U}}^{-},\mathbf{n}\right)$$

The approximate Riemann solver $\mathcal{H}$ is based on a local Lax–Friedrichs splitting of the form:

$$\mathcal{H}\left({\mathbf{U}}^{+},{\mathbf{U}}^{-},\mathbf{n}\right)={\displaystyle \frac{1}{2}}\left[{\mathbf{F}}^c\left({\mathbf{U}}^{+}\right)+{\mathbf{F}}^c\left({\mathbf{U}}^{-}\right)-\vartheta {\parallel a\parallel }_{max}\left({\mathbf{U}}^{-}-{\mathbf{U}}^{+}\right)\right],\vartheta \ge 0$$

where ${∥a∥}_{max}$ is the maximum (absolute value) of the eigenvalues of the Jacobian matrix associated with the normal flux ${\mathbf{F}}^c\left(\mathbf{U}\right)·\mathbf{n}$ and ϑ is an upwind parameter. The numerical solution of the weak formulation may correspond to either a physical mode or a spurious mode depending on the value of ϑ, as demonstrated by Ainsworth [21]. The local Lax–Friedrichs formula, with $\vartheta =1$, one of the simplest Riemann flux formulations, is commonly used in many DGM implementations because of its low operational cost. In this case, the spurious mode is damped so that it seldom has an influence.

For the diffusive flux, the splitting is based on the interior penalty method [22], evaluating the residual using only the neighboring points of the element.

With Expansions Equations (8) and (10), setting:

$${\mathbf{F}}_h^k·\mathbf{n}=\sum _{l=1}^{N_{edge}}{\mathbf{f}}_l^k{\overline{L}}_l^k·\mathbf{n}$$

where $N_{edge}$ is the number of the nodes on the element edges and ${\overline{L}}_l^k$ the edge polynomial basis, the weak formulation Equation (11) reads:

$$\begin{array}{ccc}\hfill {\int }_{{\Omega }^k}{L}_i^k\left(\sum _{l=1}^N{\displaystyle \frac{\mathrm{d}{\mathbf{u}}_l^k}{\mathrm{d}t}}{L}_l^k\right)d{\Omega }^k& -& {\int }_{{\Omega }^k}\nabla L_i^k\left(\sum _{l=1}^N{\mathbf{f}}_l^k{L}_l^k\right)d{\Omega }^k+{\int }_{{\Sigma }^k}{L}_i^k\left(\sum _{l=1}^{N_{edge}}{\mathbf{f}}_l^k{\overline{L}}_l^k·\mathbf{n}\right)d{\Sigma }^k\hfill \\ & =& {\int }_{{\Omega }^k}{L}_i^k\left(\sum _{l=1}^N{\mathbf{s}}_l^k{L}_l^k\right)d{\Omega }^k\hfill \end{array}$$

(14)

which holds for every element of the domain.

Similarly, Equation (6) can be written in a weak formulation using the same spatial discretization and basis functions:

$${\int }_{{\Omega }^k}{L}_i^k\left(\sum _{l=1}^N{\mathbf{q}}_l^k{L}_l^k\right)d{\Omega }^k+{\int }_{{\Omega }^k}\nabla L_i^k\left(\sum _{l=1}^N{\mathbf{u}}_l^k{L}_l^k\right)d{\Omega }^k-{\int }_{{\Sigma }^k}\left(\sum _{l=1}^{N_{edge}}{\mathbf{u}}_l^k{\overline{L}}_l^k\right)L_i^k·\mathbf{n}d{\Sigma }^k=0$$

(15)

In principle, all of the integrals in Equations (14) and (15) may be evaluated over each single element ${\Omega }^k$ of the domain. However, it is convenient to map every element to a master element ${\Omega }_M$ defined on the Cartesian reference system $(\xi ,\eta )$. For triangular elements, the master element is the unit right triangle with vertices $({\xi }_1,{\eta }_1)=(0,0)$, $({\xi }_2,{\eta }_2)=(1,0)$ and $({\xi }_3,{\eta }_3)=(0,1)$. For quadrangular elements, the master element is a square centered in the origin of the $(\xi ,\eta )$ system, with sides parallel to the reference axes and length $l_{side}=2$. The mapping between ${\Omega }^k$ and ${\Omega }^M$ can be expressed in the form:

$$\begin{array}{ccc}\hfill x(\xi ,\eta )& =& \sum _{m=0}^N\sum _{n=0}^{N-m}{\vartheta }_{mn}{\xi }^m{\eta }^n\hfill \end{array}$$

(16)

$$\begin{array}{ccc}\hfill y(\xi ,\eta )& =& \sum _{m=0}^N\sum _{n=0}^{N-m}{\psi }_{mn}{\xi }^m{\eta }^n\hfill \end{array}$$

(17)

the coefficients ${\vartheta }_{mn}$ and ${\psi }_{mn}$ are solution of a linear system, which is always well conditioned, because the terms ${\sum }_{m=0}^N{\sum }_{n=0}^{N-m}{\xi }^m{\eta }^n$ form a linearly-independent set of bases. The mapping expressions Equations (16) and (17) enable a fast evaluation of the transformation Jacobian J and its determinant $\left|\right|J\left|\right|$. The mapping does not impose any limit on the geometric shape of the physical cell, and therefore, it is suited for the description of curvilinear elements.

Equations (14) and (15) are transformed on the master element ${\Omega }^M$ as follows:

$$\begin{array}{ccc}\hfill {\int }_{{\Omega }^M}{L}_i^k\left(\sum _{l=1}^N{\mathbf{q}}_l^k{L}_l^k\right)\left|\right|J^k\left|\right|d{\Omega }^M& +& {\int }_{{\Omega }^M}\nabla L_i^k\left(\sum _{l=1}^N{\left(J^k\right)}^{-1}{\mathbf{u}}_l^k{L}_l^k\right)\left|\right|J^k\left|\right|d{\Omega }^M\hfill \\ \mathrm{(18)}& & -& {\int }_{{\Sigma }^M}{L}_i^k\left(\sum _{l=1}^{N_{edge}}{\left(J^k\right)}^{-1}{\mathbf{u}}_l^k{\overline{L}}_l^k\right)·\mathbf{n}\left|\right|\partial J^k\left|\right|d{\Sigma }^M=0\hfill \end{array}$$

$$\begin{array}{ccc}\hfill {\int }_{{\Omega }^M}{L}_i^k\left(\sum _{l=1}^N{\displaystyle \frac{\mathrm{d}{\mathbf{u}}_l^k}{\mathrm{d}t}}{L}_l^k\right)\left|\right|J^k\left|\right|d{\Omega }^M& +& {\int }_{{\Sigma }^M}{L}_i^k\left(\sum _{l=1}^{N_{edge}}{\left(J^k\right)}^{-1}{\mathbf{f}}_l^k·\mathbf{n}\right){\overline{L}}_l^k\left|\right|J^k\left|\right|d{\Sigma }^M\hfill \\ & -& {\int }_{{\Omega }^M}\nabla L_i^k\left(\sum _{l=1}^N{\left(J^k\right)}^{-1}{\mathbf{f}}_l^k{L}_l^k\right)\left|\right|J^k\left|\right|d{\Omega }^M\hfill \\ \mathrm{(19)}& & =& {\int }_{{\Omega }^M}{L}_i^k\left(\sum _{l=1}^N{\mathbf{s}}_l^k{L}_l^k\right)\left|\right|J^k\left|\right|d{\Omega }^M\hfill \end{array}$$

For each element, the following matrices are evaluated:

$$\begin{array}{ccc}\hfill {\int }_{{\Omega }^k}{L}_i^k{L}_l^{k}d{\Omega }^k& =& {\int }_{{\Omega }^M}{L}_i^k{L}_l^k\left|\right|J^k\left|\right|d{\Omega }^M={\mathbf{M}}_{\mathbf{e}}\hfill \\ \hfill {\int }_{{\Omega }^k}\nabla L_i^k{L}_l^{k}d{\Omega }^k& =& {\int }_{{\Omega }^M}\nabla L_i^k{L}_l^k{\left(J^k\right)}^{-1}\left|\right|J^k\left|\right|d{\Omega }^M={\mathbf{D}}_{\mathbf{e}}\hfill \\ \hfill {\int }_{{\Sigma }^k}{L}_i^k{\overline{L}}_l^{k}d{\Sigma }^k& =& {\int }_{{\Sigma }^M}{L}_i^k{\overline{L}}_l^k\left|\right|\partial J^k\left|\right|d{\Sigma }^M={\mathbf{B}}_e\hfill \end{array}$$

where ${\mathbf{M}}_{\mathbf{e}}$ is the mass matrix, of size $N\times N$, ${\mathbf{D}}_{\mathbf{e}}$ the stiffness matrix of size $N\times N\times N_{var}$ and $N_{var}$ the number of acoustic variables. The matrix ${\mathbf{B}}_e$ is of size $N\times N_{edge}$. For the particular choice of the basis and test functions, the mass matrix is reduced to a diagonal matrix; the terms of the mass matrix consist simply of the geometric Jacobian of the element multiplied by the weight function of the integration method.

The base functions are evaluated over the master element. They are the Lagrangian polynomials defined on the node set $T_p=\{x_j,j=1,N\}$, where N is the number of nodes. For rectangular elements, the bases are obtained as tensor product of one-dimensional Lagrangian polynomials defined on the Gauss-l-Lobatto nodes. ${\varphi }_l\left(\xi \right)$, with $l=1,N_{\xi }$, and ${\varphi }_r\left(\eta \right)$, $l=1,N_{\eta }$, being one-dimensional Lagrangian polynomials, the two-dimensional polynomials are defines as:

$${\mathsf{\Phi }}_i(\xi ,\eta )={\varphi }_l\left(\xi \right){\varphi }_r\left(\eta \right),i=1,N_{\xi }{N}_{\eta }$$

For triangular elements, the Lagrangian polynomials are constructed on a set of nodes, whose locations are the minimum energy solution of a steady electrostatic problem [20]. The nodes along the edges correspond to the one-dimensional Gauss–Lobatto quadrature points.

3.2. Time Integration

Time integration is performed using a fourth-order, six-stage Runge–Kutta scheme, which has low dispersion and dissipation errors [13].

Classical third- and fourth-order Runge–Kutta schemes provide relatively large stability limits, but for acoustic calculations, the stability consideration alone is not sufficient, since the Runge–Kutta schemes retain both dissipation and dispersion errors. Hu et al. [23] have shown that to get time accurate solutions in wave propagation problems, time steps much smaller than those allowed by the stability limit of the classical Runge–Kutta schemes must be used. This constraint certainly undermines the efficiency of the classical integration schemes. Instead of choosing the coefficients of the Runge–Kutta scheme to optimize the maximum order of accuracy, it is possible to select coefficients, so as to minimize the dissipation and the dispersion errors. Moreover, this optimization does not introduce additional stability constraints, and sufficiently large time steps can be used, which therefore increase the efficiency of the computation.

The Runge–Kutta scheme is implemented using the low-storage Willamson’s formulation, which only requires two storage locations per variable [24].

3.3. Boundary Conditions

To avoid spurious reflections from far-field boundaries, sponge layers are added to the computational domain. In the current formulation, sponge layers act by both damping the solution and applying a virtual stretch to the computational mesh [14].

The damping is applied by multiplying the solution by a smoothing function, which gradually decreases from one to zero:

$${\mathbf{u}}_{\text{sponge}}=\zeta \mathbf{u}$$

where ${\mathbf{u}}_{\text{sponge}}$ is the solution in the sponge layer and ζ is the damping function. The function ζ is defined as [25]:

$$\zeta =(1-C_1{x}_l^2)\left(1-{\displaystyle \frac{1-e^{C_2{x}_l^2}}{1-e^{C_2}}}\right)$$

where $C_1=0$ and $C_2=13$. The quantity $x_l$ is the normalized distance from the inner border of the sponge layer. Thus, $x_l$ ranges from zero to one, marking, respectively, the beginning and the end of the sponge layer. The virtual stretching has the effect to gradually slow down waves in the layer; this can be achieved with a coordinate transformation. For a vertical layer, the transformation $x=x\left(\xi \right)$ is defined as the inverse solution from the virtually-stretched coordinates $\xi =\left[{\xi }_1,{\xi }_2\right]$, with ${\xi }_1$ and ${\xi }_2$ being the starting and the ending coordinates of the virtually-stretched layer, of the ordinary differential equation:

$${\displaystyle \frac{\mathrm{d}x}{\mathrm{d}\xi }}=\eta \left(x\left(\xi \right),\xi \right),x\left({\xi }_2\right)=x_2$$

where $x_2$ is the ending coordinate of the physical layer. The stretching function $\eta \left(x\right)$ can be chosen as:

$$\eta =1-\left(1-{\epsilon }_l\right){\left[1-{\left(1-x_l\right)}^p\right]}^q$$

where $p=3.25$, $q=1.75$ and ${\epsilon }_l={10}^{-4}$.

The combination of the two strategies leads to an efficient approach that requires short buffer layers and is computationally efficient.

Walls are assumed impermeable and acoustically rigid; this means that no flow passes through the boundary and that acoustic waves are totally reflected. The flux normal to a wall is evaluated using the relation:

$${\left({\mathbf{F}}^c+{\mathbf{F}}^d\right)}_{w,n}=\left({\mathbf{F}}^c+{\mathbf{F}}^d\right)·\mathbf{n}$$

where $\mathbf{n}$ is the unit vector normal to the wall and ${\mathbf{F}}^c$ is evaluated with Equation (2) and ${\mathbf{F}}^d$ with Equation (3) imposing wall conditions for pressure and velocity fluctuations. Two different boundary conditions can be imposed at the walls: slip and no-slip boundary conditions. The slip flow boundary condition forces the velocity to be tangent to the wall. The velocity on the wall at time step $j+1$ can be therefore evaluated from the values at the previous time step as:

$${\left({\mathbf{u}}_w^{\prime }\right)}^{j+1}={\left({\mathbf{u}}_w^{\prime }\right)}^j-\left[{\left({\mathbf{u}}_w^{\prime }\right)}^j·\mathbf{n}\right]\mathbf{n}$$

whereas the pressure fluctuations at the wall are evaluated linearizing the exact solution of the Riemann problem for a reflective wall [26]:

$${\left(p_w^{\prime }\right)}^{j+1}={\left(p_w^{\prime }\right)}^j+p_0\frac{\gamma }{c_0}{\left(u_w^{\prime }\right)}^j$$

In a similar way, the no slip boundary condition can be imposed forcing all of the components of the velocity to be zero:

$${\left({\mathbf{u}}_w^{\prime }\right)}^{j+1}=0$$

For the LNS equations, an additional boundary condition should be imposed on the temperature fluctuations $T^{\prime }$. Assuming there is no heat flux normal to the wall (adiabatic walls), the spatial gradient of the temperature on a wall can be evaluated as:

$${\left(\nabla T_w^{\prime }\right)}^{j+1}={\left(\nabla T_w^{\prime }\right)}^j-\left[{\left(\nabla T_w^{\prime }\right)}^j·\mathbf{n}\right]\mathbf{n}$$

4. Line Source in an Incompressible Linear Shear Flow

As the first numerical test, the case of a line source in a low speed mean flow with shear is studied. As pointed out by Brambley et al. [27] and Rienstra et al. [15], if a source is located inside a mean flow with shear, a non-modal contribution appears. This field is not present when the mass source is within a uniform flow; therefore, the existence of this non-modal mode is independent of the source. The non-modal mode can be identified as a wake, of hydrodynamic nature, due to the presence of the mass source inducing harmonic isentropic perturbations. This contribution is analogous to the vortex stretching. In the case of a two-dimensional parallel sheared flow, with constant density ${\rho }_0$ and sound speed $c_0$, no vortex stretching may occur. Therefore, in the absence of an external force, the vorticity is conserved, and the particle vorticity changes, due to the mass source, can only be a redistribution of the vorticity, because there is no vorticity production.

Considering the effect of a time-harmonic monopole line (along the third dimension) source on an incompressible inviscid two-dimensional mean shear flow in an infinite domain, the correct representation of its associated trailing vorticity is a challenging test for the LEE and LNS equations.

It can be shown [15] that in an unbounded parallel linearly-sheared mean flow with a point monopole located at the origin and time dependence e${}^{i\omega t}$, the pressure field, the solution of the incompressible form of the Pridmore–Brown equation [1], has a singularity at wave number $k_0=\omega /U_0$, where $U_0>0$ is the velocity of the mean flow at the origin. This singularity, being a pole, does not represent a mode. It is a semi-infinite vorticity sheet waving with wavenumber $k_0$. Since $k_0$ is real and positive, the contribution is of constant amplitude in xand is found to exist for $x>0$. It corresponds to the trailing vorticity behind the mass source.

For the numerical test, the harmonic monopole is a Gaussian function centered at the origin:

$$S=ϵ\exp \left[-\ln \left(2\right)·\left(\frac{x^2+y^2}{b^2}\right)\right]\text{sin}\left(\omega t\right)$$

where $ϵ=1000$ and $b=0.125$. The domain of integration is a rectangle $\left[-4,6\right]\times \left[-3,3\right]\mathrm{m}$, $U_0=3\mathrm{m}/\mathrm{s}$, and $\sigma =\mathrm{d}U/\mathrm{d}y=61/\mathrm{s}$. The angular frequency is $\omega =8\mathrm{rad}/\mathrm{s}$. The corresponding hydrodynamic wavenumber is $k_0=2.667$.

In Figure 1, the snapshot of the pressure field, computed with the LEE, is shown. The field is the sum of the acoustic modes and of the hydrodynamic component. The trailing vorticity solution is obtained subtracting from the complete field the acoustic component, which is obtained as a solution of the APE equations. In Figure 2, the velocity components and the vorticity are shown. From a visual inspection, it can be seen that the computed wavenumber of the vortical wake is approximately $k=2.35$. The small difference with respect to the analytical value $k_0$ may be due to the difference in the spatial extension of the monopole source. The analytical point source has been approximated with a Gaussian distribution.

The same results have been obtained solving the LNS equations. This is not surprising, because in this case, there is no vorticity production for viscous effects. However the test case enables us to compare the computational effort in the two cases. The computational domain, discretized with 960 elements, is surrounded by sponge layers. The complete computational domain, the $\left[-4,6\right]\times \left[-3,3\right]$ rectangle surrounded by sponge layers of 3 m in depth, is discretized with 3300 elements. Both computations ended at time $t=5$ s, with elements of order five and a CFLnumber of 0.25. The computational time for the LEE was about 290 min, while the LNS computation took about 480 min on eight processors: Intel(R) Xeon(R) CPU E5440 (2.83 GHz).

5. Scattering Matrix for a Sudden Area Discontinuity in a Cylindrical Duct in the Presence of a Mean Flow

The second test case deals with the analysis of the acoustic propagation of an incoming perturbation inside a circular duct with a sudden area expansion in the presence of a mean flow and the evaluation of its scattering matrix. In this case, a strong interaction between acoustic and vortical modes occurs. At the corners, the acoustic oscillations are strongly affected by viscous effects. Vorticity is generated at the wall corners and convected into the duct by the mean flow field.

No analytical method is able to predict the coefficients of the scattering matrix of such an acoustic network with a satisfying level of accuracy, especially for mean flows characterized by Mach numbers greater than 0.1 [28,29]. Therefore, the development of an accurate computational model for this kind of problem is of great practical interest. Previous numerical simulations of this test case [30], made with the LEE, failed to correctly predict the vorticity production. In the absence of the viscous terms, the Kutta condition at the corner was generated by numerical viscosity enhanced by the geometric singularity. Moreover, for high Mach numbers of the mean flow, Kelvin–Helmholtz instabilities prevented the convergence of the calculation. To have a realistic representation of the involved phenomena, it is necessary to include the viscous terms, as made in [31] in the frequency domain, solving the LNS equations. Experimental data have been provided by Ronneberger [16] for different mean flow Mach numbers, namely $M_0=0.08$, $0.19$ and $0.29$. The measurements were made on a cylindrical area discontinuity with an upstream diameter $h_1=50\mathrm{mm}$ and a downstream diameter $h_2=85\mathrm{mm}$, corresponding to an area ratio equal to $\eta =0.346$. In the computational domain, the length of the upstream and downstream ducts are $l_1=1.5\mathrm{m}$ and $l_2=2.75\mathrm{m}$, respectively.

In the upstream and downstream regions of the area expansion, a quadrilateral structured grid is used, while around the area discontinuity, the mesh is finer and unstructured with the exception of the region near the corner, where no-slip boundary conditions are imposed. In this part of the domain, the mesh is kept structured. Since the acoustic field is axisymmetric, only half of the domain is considered. A close-up of the mesh around the area enlargement is shown in Figure 3.

The LNS equations are solved using elements of degree $p=4$. High-order elements enable one to perform the computations on coarser grids with respect to low-order methods [31].

To absorb the outgoing wave and avoid spurious reflections, two sponge layers with a thickness of $0.75\mathrm{m}$ are placed at the inflow and outflow boundaries.

To calculate the perturbed field, it is necessary to perform a preliminary evaluation of the mean flow field in the duct. The fluid-dynamic analysis is performed with the open source software OpenFoam, solving the steady-state incompressible Reynolds-averaged Navier–Stokes (RANS) equations with a $k-\omega$ model. The simulations of the mean flow are made imposing velocity profiles at the inflow corresponding to the Mach numbers $0.08$, $0.19$ and $0.29$, according to the above-mentioned experimental data. In all simulations, the inlet turbulent length scale is set equal to 1 mm and the turbulent intensity to $10\%$, while a fixed ambient pressure boundary condition is imposed at the outlet boundary.

5.1. Time Domain Wave Packet

In frequency-domain calculations, the acoustic sources are modeled as sinusoidal functions with a specified frequency. In the time-domain approach, the drawback of the single frequency analysis is that, before starting to collect numerical data, the simulation should reach a time periodic state. This could take a long time, leading to computationally-expensive calculations. However, for time-domain integrations, it is possible to reduce the computational cost of the simulations; the monochromatic waves can be substituted with a time domain wave packet $\mathsf{\Phi }\left(t\right)$ [32]. This approach, suitable only for linear problems, replaces the monochromatic sinusoidal sources with a single temporally-compact broadband pulse. Such a source will generate a wave packet that contains a broad range of frequencies in a short time duration (Figure 4). With this approach, it is possible to solve with one computation all frequencies within numerical resolution. The wave packet can be generated using the expression, with $\overline{t}=t-t_0$,

$$\mathsf{\Phi }\left(t\right)=\left\{\begin{array}{cc}{\displaystyle \frac{sin\left(2\pi f_{max}\overline{t}\right)}{2\pi f_{max}\overline{t}}}{e}^{\left(ln\delta \right){\left(\overline{t}/t_0\right)}^2}\hfill & \text{if }|\overline{t}|\le t_0\hfill \\ 0\hfill & \text{if }|\overline{t}|>t_0\hfill \end{array}\right.$$

(20)

where $f_{max}$ is the maximum frequency of the packet, δ is a constant equal to 0.01 and $t_0$ is half of the temporal duration of the wave packet. To generate the wave packet inside a duct, a source term is specified in the continuity with the form, for a horizontal duct,

$$h_{\text{WP}}(x,t)=\mathsf{\Phi }\left(t\right)e^{-\left(ln2\right)\left(x^2-x_0^2\right)/b^2}$$

(21)

where $x_0$ is the location of the source mid-plane inside the duct and b is the half-width of the source region. The parameters for all of the simulations are: $f_{max}=4000\mathrm{Hz}$, $b=0.025\mathrm{m}$ and $t_0=5\times {10}^{-4}\mathrm{s}$. The plane $x_0$, for both left- and right-running wave packets, is placed inside the duct near the interior boundary of the corresponding sponge layer.

5.2. Acoustic and Vortical Modes

The mean flow can be divided into three parts: the unperturbed channel flows upstream and downstream of the discontinuity and the region close to the discontinuity characterized by the presence of a jet-like flow.

In presence of sharp edges, part of the energy associated with the acoustic modes is redistributed to the vortical ones, because of the strong viscous effects due to the large velocity gradients present near the corner. The acoustic Reynolds number being high, the vortical motion developing past the corner is unstable, subject to the Kelvin–Helmholtz instability. In Figure 5, the case of a wave packet, entering from the inflow boundary, is shown at the time when it is passing on the corner. It is possible to see, behind the section enlargement, the rise of a pattern of vortical structures convected by the mean flow.

The effect of the mean flow can be seen in Figure 6, where the instantaneous vorticity perturbation is plotted at the same time, $t={5.2510}^{-3}$ s, for the three Mach numbers, $0.08$, $0.19$ and $0.29$. With increasing mean flow velocity, the convection speed of the eddies increases, but the vorticity intensity levels are almost the same, as expected, because the hydrodynamic mode is convected with the mean flow velocity, but the mechanism of vorticity generation depends only on the imposed perturbation, which is the same in the three cases.

5.3. Scattering Matrix

The relationship between the transmission and the reflection of incoming sound waves through an acoustic element can be described using the scattering-matrix formalism. A general network element can be thought of as a black box to which acoustic waves enter and exit through the so-called ports. For a two-port element, the scattering-matrix formalism could be written as [33,34]:

$$\left(\begin{array}{c}{p}_{a-}\\ p_{b+}\end{array}\right)=\underset{\mathit{Passive}\mathit{part}}{\underbrace{\stackrel{\mathbf{S}}{\overbrace{\left(\begin{array}{cc}{R}^{+}& T^{-}\\ T^{+}& R^{-}\end{array}\right)}}\left(\begin{array}{c}{p}_{a+}\\ p_{b-}\end{array}\right)}}+\underset{\mathit{Active}\mathit{part}}{\underbrace{\left(\begin{array}{c}{p}_{a+}^s\\ p_{b-}^s\end{array}\right)}}$$

(22)

where a and b represent the inflow and outflow ports, respectively, and $p_{+}$ and $p_{-}$ are the Fourier coefficients of the transmitted and reflected acoustics waves, defined as in Figure 7. The matrix $\mathbf{S}$ is the scattering matrix and takes into account the passive transmission and reflection of sound waves. The active processes that generates sound within the element are considered in the second term of Equation (22). In the case of the duct with area expansion, all elements are assumed passive, i.e., they do not generate sound on their own; therefore, the scattering-matrix formalism of a two-port becomes:

$$\left(\begin{array}{c}{p}_{a-}\\ p_{b+}\end{array}\right)=\left[\mathbf{S}\right]\left(\begin{array}{c}{p}_{a+}\\ p_{b-}\end{array}\right)$$

To evaluate the scattering matrix $\mathbf{S}$, the two-source location method could be applied: first, the problem is solved with a time-harmonic wave entering from the inflow, Configuration I (Figure 7(a)), and then, the problem is solved with the same wave entering from the outflow port, Configuration $II$ (Figure 7(b)). With this approach, it is possible to obtain two linearly-independent complex equations for the four unknown quantities of the scattering matrix:

$$\mathbf{S}=\left[\begin{array}{cc}{R}^{+}& T^{-}\\ T^{+}& R^{-}\end{array}\right]=\left[\begin{array}{cc}{p}_{a-}^I& p_{a-}^{II}\\ p_{b+}^I& p_{b+}^{II}\end{array}\right]{\left[\begin{array}{cc}{p}_{a+}^I& p_{a+}^{II}\\ p_{b-}^I& p_{b-}^{II}\end{array}\right]}^{-1}$$

To calculate the elements of the scattering matrix, it is necessary to know the Fourier coefficients of the transmitted and reflected acoustic waves at each port. These quantities are not directly available from the numerical solution, because the solution represents the time-domain acoustic field inside the tested element. Therefore, it is necessary to post-process the solution to obtain the quantities needed for the evaluation of the scattering matrix: this could be done through a method called the plane wave decomposition method. Assuming that only plane waves propagate through the element, the Fourier coefficients of the acoustic pressure can be written as ([31,35] and [36]),

$$p\left(x,f\right)=p_{+}(x,f)+p_{-}(x,f)=|p_{+}|e^{i{\varphi }^{+}}{e}^{-i{\kappa }^{+}x}+\left|p_{-}\right|e^{i{\varphi }^{-}}{e}^{-i{\kappa }^{-}x}$$

(23)

where $|p_{\pm }|$ and ${\varphi }^{\pm }$ are real quantities, which represent the amplitudes and the phases of up- and down-stream propagating waves. Neglecting the viscosity effects on the acoustic propagation, the wave numbers ${\kappa }_{\pm }$ are real quantities. For plane waves, the wave numbers ${\kappa }^{\pm }$ can be therefore written as ${\kappa }^{\pm }=\mp \omega /\left(c\pm u_0\right)$, where $u_0$ is the mean flow velocity along the duct axis.

Once the time-domain acoustic field is known, the time history of the pressure is extracted on N points along the center-line of each port and converted in the frequency domain using the Fourier transform. For every input frequency, a plane wave decomposition Equation (23) is written in all of the point locations, obtaining an overdetermined system:

$$\left\{\begin{array}{c}p\left(x_1\right)=|p_{+}|e^{i{\varphi }^{+}}{e}^{-i{\kappa }^{+}{x}_1}+\left|p_{-}\right|e^{i{\varphi }^{-}}{e}^{-i{\kappa }^{-}{x}_1}\hfill \\ ⋮\hfill \\ p\left(x_j\right)=|p_{+}|e^{i{\varphi }^{+}}{e}^{-i{\kappa }^{+}{x}_j}+\left|p_{-}\right|e^{i{\varphi }^{-}}{e}^{-i{\kappa }^{-}{x}_j}\hfill \\ ⋮\hfill \\ p\left(x_N\right)=|p_{+}|e^{i{\varphi }^{+}}{e}^{-i{\kappa }^{+}{x}_N}+\left|p_{-}\right|e^{i{\varphi }^{-}}{e}^{-i{\kappa }^{-}{x}_N}\hfill \end{array}\right.$$

(24)

which can be solved using an iterative non-linear least-square method to obtain the amplitudes and the phases of the plane waves. The iterative method needs an initial guess of the solution. This initial solution can be evaluated solving the system obtained considering the first and the last equations of System Equation (24):

$$\begin{array}{cc}\hfill {\left(|p_{+}|e^{i{\varphi }^{+}}\right)}_{\text{guess}}& ={\displaystyle \frac{p\left(x_1\right)·e^{-i{\kappa }^{-}{x}_N}-p\left(x_N\right)·e^{-i{\kappa }^{-}{x}_1}}{e^{-i{k}^{+}{x}_1-i{k}^{-}{x}_N}-e^{-i{k}^{+}{x}_N-i{k}^{-}{x}_1}}}\hfill \\ \hfill {\left(|p_{-}|e^{i{\varphi }^{-}}\right)}_{\text{guess}}& ={\displaystyle \frac{p\left(x_1\right)·e^{-i{\kappa }^{+}{x}_N}-p\left(x_N\right)·e^{-i{\kappa }^{+}{x}_1}}{e^{-i{k}^{-}{x}_1-i{k}^{+}{x}_N}-e^{-i{k}^{-}{x}_N-i{k}^{+}{x}_1}}}\hfill \end{array}$$

The scattering coefficients are reported for the three Mach numbers, $0.08$, $0.19$ and $0.29$, in Figure 8, Figure 9, Figure 10, Figure 11, Figure 12 andFigure 13. The numerical results are compared to the experimental ones [16]. A good agreement between the numerical and experimental data could be observed. This is true in particular for the lowest Mach number $M_0=0.08$, where there is an almost perfect superimposition of the results. For higher Mach numbers, small discrepancies could be noticed, especially in the module of the transmitted component $|T^{-}|$. The errors between the numerical and the experimental data observed at higher Mach numbers could be partially explained with the assumption of incompressible flow adopted during the simulation of the mean flow.

6. Conclusions

The LNS equations for the propagation of small perturbations in complex geometries have been derived for general coordinates and for axisymmetric geometries. The model is able to represent the hydrodynamic-acoustic interactions, coupling acoustic waves and vortical modes, and the mechanism responsible for the generation of vorticity associated with the hydrodynamic modes.

The LNS equations are discretized in space using a discontinuous Galerkin formulation, which is consistent and stable. It can deal with unstructured grids, and it is compact, making the method ideally suited for parallel computing.

The LNS model has been tested for two representative cases. The first one is the time-harmonic source line in an incompressible inviscid two-dimensional mean shear flow in an infinite domain. The comparison with the analytical solution shows that the proposed model is able to correctly capture the trailing vorticity field developing behind the mass source and to represent the redistribution of the vorticity.

In the second case, the acoustic propagation of an incoming perturbation inside a circular duct with a sudden area expansion in the presence of a mean flow has been studied. A method for computing the coefficient of the scattering matrix has been discussed. The numerical computations clarify the mechanism of the vorticity generation at corners, where the acoustic oscillations are strongly affected by viscous effects, and vortical perturbations are generated at the wall and convected into the duct by the mean flow field. The computed coefficients of the scattering matrix were compared to experimental data, for three different Mach numbers of the mean flow, M${}_0$ = 0.08, 0.19 and 0.29. The good agreement with the experimental data shows that the proposed method is suitable for characterizing the acoustic behavior of this kind of network.