© The Author(s), published by EDP Sciences, 2023

1 Introduction

Since a direct computation of aeroacoustic noise by a turbulent scale resolving simulation of the flow field including the generation and propagation of the sound waves is computationally expensive, two-step or hybrid, i.e., computational fluid dynamics (CFD) – computational aeroacoustics (CAA) methods are usually preferred. In these methods the flow field is computed first and subsequently, the acoustic field is determined using the results of the flow field solution. A well established hybrid CFD/CAA approach consists of a large-eddy simulation (LES) to compute the unsteady turbulent flow field and a subsequent CAA step, in which, e.g., the acoustic perturbation equations (APE) are solved to predict the propagation of the acoustic waves through the inhomogeneous near field up to the far field. In such a hybrid approach, the acoustic source terms for the CAA step are obtained from the instantaneous flow field. A major drawback of using two disjoint solvers is the necessity to transfer the volume source terms from the CFD to the CAA solver. If this is done via disk I/O, huge disk storage is required for large-scale problems. On high-performance computing hardware, the overall parallel scalability is often limited by the available I/O bandwidth, which can render a coupling over disk I/O inefficient. Therefore, a new fully-coupled direct-hybrid method has been developed, in which the LES and CAA simulations are conducted in a shared framework, i.e., both solvers are implemented in a single software framework [1]. This allows for the concurrent execution of both solvers with an in-memory data exchange of the acoustic source terms thus fully avoiding disk I/O. Furthermore, by employing a joint hierarchical Cartesian grid, which is partitioned on a coarse level via a space-filling curve, both solvers are efficiently coupled and parallelized, whereas the use of solution adaptive meshes with dynamic load balancing is possible [2]. Thus, the direct-hybrid method combines the advantages of the direct and hybrid solution strategies, i.e., minimizes I/O and exploits the different length scales in the acoustic and flow fields, such that it is a good candidate to perform large-scale simulations of three-dimensional turbulent aeroacoustic problems.

The lattice Boltzmann (LB) method is known to be an efficient CFD solver for low Mach number flows. Therefore, a direct-hybrid CFD/CAA method is implemented, in which an LB method is used for the flow field prediction being coupled to a discontinuous Galerkin scheme for the solution of the APE. The coupling strategy is similar as for a finite-volume solver for the Navier–Stokes equations, which has been developed previously in [1, 3]. For the first time, the accuracy is discussed for the two directly coupled CFD/CAA methods, i.e., the LB-DG and the FV-DG method.

Since all solvers are formulated on Cartesian meshes, boundary conditions for the immersed surfaces of solid bodies have to be formulated. While well know methods exist for the CFD methods, such boundary conditions are more difficult to formulate for the higher-order DG method. Different approaches exist for the immersed boundary (IB) treatment of solid bodies. In [4], a Cartesian ghost-cell IB method for high-order finite-difference schemes was presented. Chung and Morris [5] modeled the solid body as a region of higher impedance, which naturally scatters impinging acoustic waves. In the compressible Brinkman penalization method by Liu [6], equations in porous media modeling have been utilized to formulate a solid boundary treatment with a region of high impedance. In [7], the method was used on a boundary-fitted Cartesian-like grid and a slip wall boundary is enforced by an anisotropic penalization. Various modifications have been proposed in the literature, e.g., in [8], a Galilean invariant formulation was proposed, and in [9] the penalization method was extended to account for more general Neumann and Robin boundary conditions. In the present work, an IB method motivated by porous media modeling is applied, which is particularly designed for a conservative numerical scheme.

This paper has the following structure. First, the computational methodologies, i.e., the DG method for solving the APE extended by the IB formulation, the CFD methods, and the direct-hybrid coupling procedure, are introduced in Section 2. Subsequently, the numerical methods are validated in various setups explained in Section 3. Finally, conclusions are drawn in Section 4.

2 Computational methodologies

As previously stated, the flow and acoustic fields are predicted by direct-hybrid coupled CFD/CAA methods. The flow field is simulated by a CFD method, which is based either on a LB method or a FV method solving for the Navier–Stokes equations. Both methods can be directly coupled to a discontinuous Galerkin method, which is used to solve the APE [10] to predict the acoustic field. The noise generating source terms, i.e., the fluctuating Lamb vector for the spinning vortex case and the subsonic cylinder flow, are obtained from the instantaneous flow field and are communicated to the CAA solver. All methods are formulated on hierarchical Cartesian grids, where the coupled solvers share a joint mesh with an independent grid refinement level for each solver as described in [2]. Here, the grid spacing Δ _n for the grid refinement level n is given by Δ _n  = Δ₀/2 ⁿ , where Δ₀ is a predefined length of the level n = 0. In the following, all numerical methods are briefly discussed.

2.1 Computational aeroacoustics method

2.1.1 Discontinuous Galerkin spectral element method

The generation and propagation of acoustic waves is simulated by solving the APE-4 variant of the acoustic perturbation equations$$\frac{\mathrm{\partial }{\underline{u}}^{\mathrm{\prime }}}{\mathrm{\partial }t}+\nabla \left(\underline{\bar{u}}\cdot {\underline{u}}^{\mathrm{\prime }}+\frac{p^{\mathrm{\prime }}}{\bar{\rho }}\right)={\underline{q}}_{m,}$$(1) $$\frac{\mathrm{\partial }{p}^{\mathrm{\prime }}}{\mathrm{\partial }t}+{\bar{c}}_s^2\nabla \cdot \left(\bar{\rho }{\underline{u}}^{\mathrm{\prime }}+\underline{\bar{u}} \frac{p^{\mathrm{\prime }}}{{\bar{c}}_s^2}\right)=0,$$(2)where the time average is denoted by $\overline{(\cdot )}$ and $(\cdot {)}^{\mathrm{\prime }}=(\cdot )-\overline{(\cdot )}$ is a perturbed quantity. For aerodynamically induced noise at low Mach number, only the perturbed Lamb vector ${\underline{L}}^{\mathrm{\prime }}=(\underline{\omega }\times \underline{u}{)}^{\mathrm{\prime }}$ is retained as an acoustic source ${\underline{q}}_m\approx -{\underline{L}}^{\mathrm{\prime }}$. Equations (1) and (2) are discretized in space by a high-order nodal discontinuous Galerkin spectral element method (DGSEM). Only a brief outline of the DGSEM is given. For more details the reader is referred to [11–13] or, for the particular implementation in the simulation framework m-AIA used for this paper, to [1].

The solution procedure in d-dimensions consists of the following major steps. The equations in each grid cell are mapped from physical space $\underline{x}$ to reference space $\underline{\xi }\in [-\mathrm{1,1}{]}^d$ by $\underline{x}\left(\underline{\xi }\right)$. The equations are then multiplied by test functions and integrated over the reference element E before integration by parts is applied to get the weak formulation. The solution space in each element is approximated by the tensor product ansatz of piecewise polynomials of degree at most N, i.e., the function space ${\mathbb{P}}^N\left(E\right)$. A full set of basis functions is obtained by the Lagrange interpolating polynomials, $\left\{l_i\right(\xi ){\}}_{i=0}^N$, defined by a set of interpolation nodes $\{{\xi }_i{\}}_{i=0}^N$, i.e., $l_i\left({\xi }_j\right)={\delta }_{\mathrm{ij}}$. Hence, the solution vector is expanded as$$\underline{U}\left(\underline{\xi },t\right)\approx {\underline{U}}_h\left(\underline{\xi },t\right)=\sum _{i,j=0}^N\mathrm{}{\underline{\widehat{U}}}_{\mathrm{ij}}\left(t\right){\psi }_{\mathrm{ij}}\left(\underline{\xi }\right) \mathrm{with}\hspace{0.5em}{\psi }_{\mathrm{ij}}\left(\underline{\xi }\right)=l_i\left(\xi \right)l_j\left(\eta \right),$$(3)for d = 2 with the nodal degrees of freedom ${\underline{\widehat{U}}}_{\mathrm{ij}}$. The fluxes and the sources are approximated using the same nodal expansion as in equation (3). Due to the discontinuity of the solution space across element boundaries ∂E, the surface flux is substituted by a numerical flux, here the local Lax-Friedrichs flux. Adhering to the Galerkin ansatz, test functions $\underline{\phi }$ are also taken from ${\mathbb{P}}^N\left(E\right)$, in particular $\underline{\phi }=\underline{\psi }$. The space integrals are replaced by a numerical quadrature rule. A collocation-type discretization with integration nodes $\{{\xi }_i{\}}_{i=0}^N$ and weights $\{w{\}}_{i=0}^N$ is applied, resulting in an orthogonal approximation space and a diagonal discrete mass matrix due to $l_i\left({\xi }_j\right)={\delta }_{\mathrm{ij}}$. Here, either Legendre–Gauss nodes with an integration precision of 2N + 1 or Legendre–Gauss–Lobatto nodes with an integration precision of 2N − 1 are used. The obtained semi-discrete ordinary differential equation is integrated in time by an appropriate time integration scheme.

In [14] it was shown that for most combinations of polynomial order and accuracy bounds Gauss DGSEM is superior in terms of the required points per wavelength (PPW) over Gauss–Lobatto DGSEM. Therefore, the Gauss DGSEM is the preferred choice in this paper. To render the overall space discretization fourth-order in space in accordance with the time integration third-order basis functions are used in the following unless otherwise stated.

2.1.2 Immersed boundary formulation of acoustic perturbation equations

Since the CAA solver is based on a Cartesian mesh, all curved boundaries of the domain are immersed in the computational mesh. The accurate, efficient and robust numerical formulations for the immersed boundaries in the high-order DGSEM [1, 11, 12] will be explained in the following. The ansatz to incorporate the acoustic scattering effect of the solid surfaces is to consider solid bodies as permeable media in the limiting case of vanishing porosity ϕ → 0 and permeability κ → 0. From these assumptions physically motivated model terms can be derived, which are included in the APE. For the limiting case of a solid, non-permeable body, the form of the porous model terms can be adapted to yield a robust discrete system.

The approach presented in the following closely follows the ideas of the Brinkman penalization method (BPM) and the Impedance Mismatch Method (IMM) [5, 15]. The former, originally applied to incompressible flows with rigorous error estimates [16], has been extended to the compressible framework in [6] and also applied to acoustic scattering problems, e.g., in [7]. The principle in IMM is to simulate the solid body as a region of higher characteristic impedance Z = ρc _s , which according to acoustics theory will naturally reflect waves impinging on this region. For numerical reasons Chung [5] solved an auxiliary problem with Z ∝ 1/ρ, where the impedance is controlled by the mean density ratio, keeping the speed of sound constant to maintain the same CFL condition. In the compressible BPM formulation by Liu [6] in addition to penalizing the momentum and energy equation, the porosity is introduced through the continuity equation in porous media, effectively modeling a high impedance region, which overlaps with the ideas in IMM.

After testing various variants regarding accuracy, robustness, and ease of implementation, the following set of equations are proposed, which are valid in the whole CAA domain$$\frac{\mathrm{\partial }{\underline{u}}^{\mathrm{\prime }}}{\mathrm{\partial }t}+\nabla \left(\underline{\bar{u}}\cdot {\underline{u}}^{\mathrm{\prime }}+{\bar{c}}_s^2\frac{{\rho }_s^{\mathrm{\prime }}}{\bar{\rho }}\right)+\underline{\underline{R}} {\underline{u}}^{\mathrm{\prime }}={\underline{q}}_m,$$(4) $$\frac{\mathrm{\partial }{\rho }_s^{\mathrm{\prime }}}{\mathrm{\partial }t}+\frac{1}{\varphi }\nabla \cdot \left(\varphi \bar{\rho }{\underline{u}}^{\mathrm{\prime }}+\varphi \underline{\bar{u}}{\rho }_s^{\mathrm{\prime }}\right)=0.$$(5)These equations can be understood to be formally derived by volume averaging equations (1) and (2) over a porous medium with porosity ϕ. For the use in a conservative numerical scheme, the equivalent density fluctuation ${\rho }_s^{\mathrm{\prime }}=p^{\mathrm{\prime }}/{\bar{c}}_s^2$ in an isentropic acoustic medium corresponding to the pressure fluctuation is introduced. Similarly, since ϕ is independent of time the porosity is lumped together to form a superficial isentropic density perturbation, i.e., equations (4) and (5) are solved for $({\underline{u}}^{\mathrm{\prime }},\varphi {\rho }_s^{\mathrm{\prime }}{)}^T$. This prevents the need for a special numerical treatment of the interface jump of the porosity and allows a straight forward space discretization. The quantity $\underline{\underline{R}}=\frac{1}{\kappa }\underline{n}\otimes \underline{n}$ is the “anisotropic” permeability tensor with a unit normal vector of the solid body $\underline{n}$, introduced in [7] to enforce the slip boundary condition, compliant with the first-order APE system.

Note that it does not contradict the acoustics theory, that similar to [5] in the present method the effective density is smaller in the solid body region. This is because in the APE a part of the linearized terms are shifted to the right-hand side source vector, which inverts the effect of density variations on the wave operator. The present IB treatment degenerates to replacing the density by an effective density ${\bar{\rho }}_{\mathrm{eff}}=\varphi \bar{\rho }$ during the initialization of the mean field and applying a penalty term. Indeed, two contributions can be identified. The assumption of the mean density to be uniform across the fluid-solid interface yields the theoretical reflection coefficient for plane waves$$R\left(\omega ,\varphi ,\kappa \right)=\frac{{\tilde{Z}}_{\mathrm{solid}}-\varphi }{{\tilde{Z}}_{\mathrm{solid}}+\varphi }\in \mathbb{C}\hspace{0.5em}\mathrm{with}\hspace{1em}{\tilde{Z}}_{\mathrm{solid}}=\sqrt{1-\frac{i}{\omega \kappa }},$$(6)where $i=\sqrt{-1}$, ω is the angular frequency, and ${\tilde{Z}}_{\mathrm{solid}}$ is the relative impedance of the solid domain. Note that also a frequency independent formulation with vanishing phase error $\mathfrak{I}\left(R\right)$ can be easily derived by adding a pressure relaxation term and modifying the flux Jacobians. However, tests showed this to be numerically less robust and the poor frequency dependence for linear acoustics is expected to have a negligible effect on the acoustic field. Numerical experiments showed $\kappa \propto \mathcal{O}\left(10^{-2}\right)$ to be a good compromise between the ambiguities about the exact boundary interface for κ → ∞ and wave reflection.

2.1.3 Discussion of stability

In the fluid domain, i.e., where ϕ = 1 and κ → 0, the semi-discretized equations are advanced in time using an explicit five-stage fourth-order low-storage Runge–Kutta (ERK) scheme [17]. However, the stability for ϕ → 0 and vanishing $\kappa$ is restricted by two numerical phenomena. The first restriction emerges from the IB treatment modifying the flux vector to be a rational function and therefore, the quadrature is not exact. In the literature, different stabilization techniques for this kind of aliasing error exist, e.g. solving a split form of the equations, applying a modal filtering, which increases the numerical dissipation of the scheme or applying an overintegration technique to those terms, which for rational flux functions will not fully avoid the error. Using Gauss DGSEM neither filtering, nor overintegration have shown unconditional stability, i.e., resulted in positive eigenvalues of the space operator in the present case. The cure is to switch to Gauss–Lobatto DGSEM in a small layer around solid bodies, i.e., for ϕ < 1, which has been observed to result in an unconditionally stable semi-discrete system and apply Gauss DGSEM elsewhere. Gauss–Lobatto DGSEM includes the boundary points and satisfies the summation-by-parts (SBP) property, which is the discrete counterpart of integration-by-parts. Together with a diagonal mass matrix, which results for the collocation-type approach with $l_i\left({\xi }_j\right)={\delta }_{\mathrm{ij}}$, this enables energy and entropy proofs to be applied to the discrete scheme. Smaller porosity values ϕ result in a space operator with a larger spectral radii, resulting in a more severe time step restriction for explicit integration methods. With the requirement that the eigenvalue spectrum lies within the linear stability envelope of the ERK method, in the subsequent computations the porosity is set to ϕ = 1/15. The second restriction arises from the stiffness of the penalty terms. Here the system is split into the linear stiff part and the remainder, such that the ERK method is combined with a sixth-stage formally 2nd-order explicit, singly-diagonally implicit Runge–Kutta (ESDIRK) method also in 2N-storage format to an implicit–explicit Runge–Kutta (IMEX-RK) scheme. The low-storage algorithm and the respective coefficients of the ESDIRK scheme are given in Table 1. For the steady permeability tensor $\underline{\underline{R}}$, which contains only local geometric information, this results in a single inversion of a local coefficient matrix $\underline{\underline{I}}-b\Delta t\underline{\underline{R}}$, precomputed during initialization and reused in the RK sub-stages 2–6.

2.2 Computational fluid dynamics methods

2.2.1 Cumulant lattice Boltzmann method

The Boltzmann equation describes the evolution in time t of the momentum distribution function $f(\underline{x},\underline{v},t)$. This function represents the density of particles at the position $\underline{x}$ and time t with a velocity $\underline{v}$. The Boltzmann equation reads$$\frac{\mathrm{\partial }f}{\mathrm{\partial }t}+\underline{v}\cdot \nabla f=\mathrm{\Omega }\left(f\right),$$(7)with the collision operator $\mathrm{\Omega }$ as its source term accounting for the effect of the momentum exchange of particle collisions on the distribution function f. Enskog demonstrated [18] through a series expansion of the momentum distribution function with the perturbation parameter $ϵ$, i.e., $f=f^{\left(0\right)}+ϵf^{\left(1\right)}+\mathcal{O}\left({ϵ}^2\right)$, that a first-order approximation of the Boltzmann equation recovers the Navier–Stokes equations. Discretizing equation (7) yields the lattice Boltzmann equation$$f_i(\underline{x}+c_i\Delta t,t+\Delta t)=f_i^{\mathrm{*}}(\underline{x},t)=f_i(\underline{x},t)+{\mathrm{\Omega }}_i\left(f\right),$$with $c_i=v_i/\sqrt{3}$ denoting the discrete particle velocity in the discrete direction i, and the asterisk ($\mathrm{*}$) indicating the post-collision state. In the present study, the particle velocity space is discretized in a three-dimensional Cartesian lattice featuring 27 discrete velocity directions (D3Q27). Macroscopic flow quantities such as the density ρ and the flow velocity $\underline{u}$ are obtained from the moments of the momentum distribution function. In discrete form, the integrals read$$\rho \left(\underline{x},t\right)=\sum _i\mathrm{}{f}_i\left(\underline{x},t\right)\hspace{0.5em}\mathrm{and}\hspace{0.5em}\underline{u}(\underline{x},t)=\frac{1}{\rho }\sum _i\mathrm{}{\underline{c}}_i{f}_i(\underline{x},t) .$$The collision step can be performed, e.g., by using the Bhatnagar–Gross–Krook (BGK) operator [19] or by a cumulant collision operator [20]. While the former is widely used for low Reynolds number flow and is computationally less expensive, the later has been recently developed for a stable and accurate prediction of high Reynolds number flow. The BGK operator reads$${\mathrm{\Omega }}_i={\omega }_{\mathrm{BGK}}[f_i^{\mathrm{eq}}-f_i],$$with $f_i^{\mathrm{eq}}$ being the Maxwell equilibrium distribution function and ${\omega }_{\mathrm{BGK}}$ being the non-dimensional relaxation frequency given by$${\omega }_{\mathrm{BGK}}=\frac{\Delta t{c}_s^2}{{\nu }_{\mathrm{eff}}+\frac{1}{2}\Delta t{c}_s^2},$$(8)where $\Delta t$ denotes the time step and ${\nu }_{\mathrm{eff}}$ the effective viscosity. It represents the sum of the fluid viscosity and an artificially introduced viscosity, which is only used in a sponge region close to the domain boundaries to damp reflected waves. The equilibrium distribution function formulated using an isothermal assumption reasonable for low Mach number flows is$$f_i^{\mathrm{eq}}=w_i\rho \left[1+\frac{{\underline{c}}_i\cdot \underline{u}}{c_s^2}+\frac{\underline{u} \underline{u}:\left({\underline{c}}_i{\underline{c}}_i-c_s^2\underline{\underline{I}}\right)}{2c_s^4}\right],$$(9)where the weighting factors $w_i$ are equal to 8/27, 2/27, 1/54, and 1/216 for the rest, the six Cartesian, the twelve cubic edge-diagonal, and the eight cubic space-diagonal directions, respectively.

For the cumulant collision operator the momentum distribution functions is transformed into so called countable cumulants c _α , which are observable quantities that are Galilean invariant and statistically independent [20]. The collision in cumulant space reads$$c_{\alpha }^{\mathrm{*}}=c_{\alpha }+{\omega }_{\alpha }[c_{\alpha }^{\mathrm{eq}}-c_{\alpha }],$$where $c_{\alpha }^{\mathrm{eq}}$ is the Maxwell equilibrium in cumulant space and ${\omega }_{\alpha }$ is the frequency to relax the cumulants towards equilibrium. In the present study, following the nomenclature of [20], all relaxation rates but ${\omega }_1$ are set to unity$${\omega }_1={\omega }_{\mathrm{BGK}},\hspace{0.5em}{\omega }_{\alpha }=1:\hspace{1em}\alpha \ne 1 .$$The propagation step, i.e., scattering the post-collision distribution into the neighboring lattices is always performed in momentum space. Therefore, the post-collision cumulants are transformed back into momentum space. Local grid refinement is implemented by using different relaxation frequencies for each level of refinement following the method of Dupuis and Chopard [21]. To account for compressibility effects an acoustic scaling is applied such that the computational time step scales with the grid spacing $\Delta t\sim \Delta$. From equation (8) the expression for the viscosity follows:$${\nu }_{\mathrm{eff}}=\Delta t{c}_s^2\left(\frac{1}{{\omega }_{\mathrm{BGK}}}-\frac{1}{2}\right) .$$In order to keep the viscosity constant across the grid refinement, the relaxation frequency is scaled correspondingly to $\Delta t$. For the no-slip boundary condition an interpolated bounce-back scheme [22] is implemented allowing for a second order accurate representation of arbitrary shaped geometries. Inlet (outlet) boundary condition are realized by setting an equilibrium state (eq. (9)) based on a prescribed velocity (density) and on the second order accurate extrapolated density (velocity) from the interior of the domain.

2.2.2 Finite-volume method for the solution of the Navier–Stokes equations

A finite-volume method is used to discretize the non-dimensional conservation equations of mass, momentum, and energy formulated in integral form for an arbitrary control volume V $$\frac{\mathrm{d}}{\mathrm{d}t}{\int }_{V\in \mathrm{\Omega }}^{}\underline{Q} \mathrm{d}V+{\oint }_{\mathrm{\partial }V\in \mathrm{\Omega }}^{}\left({\underline{\underline{H}}}^{\mathrm{inv}}-{\underline{\underline{H}}}^{\mathrm{vis}}\right)\cdot \underline{n} \mathrm{d}S=0 .$$(10)The vector of conservative variables is defined as $\underline{Q}={\left[\rho ,\rho \underline{u},\rho E\right]}^T$. The inviscid and viscous fluxes are given by$${\underline{\underline{H}}}^{\mathrm{inv}}-{\underline{\underline{H}}}^{\mathrm{vis}}=\left[\begin{array}{l}\rho \underline{u}\\ \rho \underline{u}\otimes \underline{u}+p\underline{\underline{I}}\\ \rho E\underline{u}+p\underline{u}\end{array}\right]-\frac{1}{\mathrm{Re}}\left[\begin{array}{l}0\\ \underline{\underline{\tau }}\\ \underline{\underline{\tau }}\cdot \underline{u}+\frac{1}{\mathrm{Pr}}\underline{q}\end{array}\right] .$$(11)The quantities E, $\underline{q}$, $\underline{\underline{I}}$, Re, and Pr denote the total energy, heat-flux vector, unit tensor, Reynolds number, and Prandtl number. The discrete approximation is formulated for a Cartesian mesh with immersed boundaries. A fully conservative cut-cell method is used for the formulation of boundary conditions at the domain boundaries [23]. An explicit low storage 5-stage Runge–Kutta scheme is applied to advance the solution in time. The overall method is second-order accurate in time and space. For more details of this method, the reader is referred to [24].

2.3 Coupling algorithm in the direct-hybrid CFD/CAA method

In the direct-hybrid method, a joint hierarchical Cartesian grid is used for the CFD and the CAA solver [1]. Mesh refinement is performed during the grid generation process individually according to constraints specified for each solver. All cells are tagged according to their assignment to the solvers. Cells are organized in an octree data structure with parent, child, and neighbor relationships [25]. Thus, each coarse partition cell becomes the root cell of a subtree of the whole grid. The partitioning of the computational grid is conducted on this coarse level, where a Hilbert space-filling curve (SFC) is used to obtain a one-dimensional ordering of all partition cells. A computational workload is assigned to each cell of the joint grid depending on its use by one or both solvers. By traversing the subtrees of the grid, the accumulated workload for each coarse partition cell is determined. A domain decomposition is performed by splitting the Hilbert curve such that equal-sized partition with respect to the workload are obtained. Each partition, i.e., a spatially compact amount of partition cells with its corresponding subtrees of cells, is then allotted to a dedicated parallel process.

The partitioning strategy, which operates on the joint Cartesian mesh, facilitates an efficient coupling between the solvers, since all CFD and CAA cells contained in the same physical volume of the three-dimensional domain are assigned to the same process allowing in-memory exchange of the data containing the acoustic source terms. Spatial interpolation from the CFD to the CAA solution representation is performed by local Galerkin projection [26], which can handle arbitrary spatial mappings while requiring only data that is locally available. Intra-process spatial compactness of cells is implicitly ensured by the locality property of the Hilbert SFC, which generates compact surfaces of the subdomains and thus reduces the communication overhead. In this study, the time step used in the FV CFD simulation is chosen equal to the time step used in the DG CAA simulation. Thus, a temporal interpolation of source term data, e.g., by quadratic interpolation [27], is not required. The parallel coupled algorithm alternately advances the Runge–Kutta time integration stages of both solvers. A more detailed description of the coupling approach can be found in [1]. For the coupling of the LB-DG method the procedure slightly differs. First, the solvers alternately advance the solution by a full time step, since the LB method is not based on a multi-stage Runge–Kutta time integration scheme, see Section 2.2.1, making it more difficult to hide the communication of the CAA solver behind the CFD computation. Second, the time step of the CAA solver is chosen to match the time step of a specific LB refinement level. Therefore, the time steps used for the LB method on coarser levels than this specific grid level are multiples of the CAA time step. In the present study, the specific refinement level is set to the finest level in the LB method and acoustic source terms calculated by the LB method on the coarser levels are kept constant for the DG method until a matching time level for the corresponding cells is reached.

The interleaved execution pattern prevents additional overhead due to differing load compositions generated by a varying number of CFD and CAA cells on each partition and allows to overlap communication with computations. However, this requires the accumulated workload of the CFD and the CAA solver to be identical on each parallel subdomain. In general, the parallel efficiency of the overall computation is determined by the initial domain decomposition, which is based on estimates regarding the different computational costs for CFD and CAA cells.

3 Results

The solid body treatment of the CAA in the DG solver and the two-step direct-hybrid FV-DG and LB-DG methods are evaluated by a number of generic well-documented benchmark cases and where possible compared against analytical solutions. The complexity of the cases increases from a CAA simulation alone with a prescribed source towards fully-coupled simulations, where sources are extracted from the simultaneously running CFD simulation.

3.1 Diffraction of acoustic waves around a thin flat plate

A challenging benchmark case for the IB condition described in Section 2.1.2, is the scattering and diffraction of acoustic waves by a thin plate as presented in [28] to test their wall boundary condition in the framework of high-order finite difference-schemes. The geometric arrangement is depicted in Figure 1. The origin of the coordinate system is in the center on the upper surface of the flat plate of width L. The waves are generated by a harmonic source of the form$$S_p=A_p\exp \left\{-\left(\ln 2\right)\frac{x^2+(y+L{)}^2}{(L/8{)}^2}\right\}\cos \left(\omega t\right)$$(12)applied to the perturbation pressure equation with A _p  = 0.01 and wavelength λ = L/4, i.e., $\omega =0.32\pi$. The computational domain is defined by $(x,y)\in [-4L,4L]\times [-3.5L,4.5L]$. Data is post-processed along the line labeled “output”, which is located 4L above the plate. The domain is shifted towards the upper side to avoid direct influences on the output line from the boundary conditions. The overall setup is symmetric about the y-axis.

Figure 1 Sketch of the computational domain. Data is sampled along the line labeled “output”.

Without more sophisticated modifications, the IB method is restricted to finite sized bodies, hence the plate thickness is set to 0.0625L, which corresponds to half the coarsest grid spacing. Simulations have been conducted on uniform grids of refinement levels n = {6, 7, 8, 9, 10}, corresponding to grid spacings Δ _n  = 200/2 ⁿ . Taking the definition $\mathrm{PPW}=\lambda (N+1)/\Delta$ this corresponds to $\mathrm{PPW}\ge 8$. For the flat plate with sides aligned to the coordinate axes a body-fitted grid can be generated for the grids with n ≥ 7, where solid wall boundary conditions can be incorporated by particular boundary fluxes. In the following, these boundary treatments will be termed “exact” or labeled by the letter “e” in the legend. The reference solution for the grid convergence study is obtained by the exact boundary treatment on refinement level n = 10.

The mean values in equations (4) and (5) are initialized for a quiescent medium with homogeneous fluid properties. At the far field, radiation boundary conditions (RBC) are applied to minimize reflections from outgoing acoustic waves. Additionally, to exclude any errors from the time integration the time step is set to 20 of the maximum possible time step, dictated by stability. To evaluate the long-time stability and accuracy, solutions are compared for t = 250. In Figure 2, the instantaneous perturbation pressure fields for the exact wall treatment on the left half and the IB method on the right half of the figure for n = 7 are shown. Qualitatively the agreement between both fields is very good. Only in the immediate vicinity of the plate minor differences are visible. These are due to the sharp boundary treatment compared to the IB method, where the boundary region for this test case extends over $\approx 0.25\Delta$. A more detailed quantitative comparison among the different grids is shown by the pressure signal distributions along the output line in Figure 3. The solution obtained with the exact boundary treatment on the finest grid (n = 10e) is taken as reference. A very good match of the amplitude and phase is achieved. However, for n = 7 differences are visible in the close-up view. The convergence is shown in Figure 4 by means of the L ² error along the output line. For reference, the errors obtained with the exact boundary treatment and n = {7, 8, 9} are also plotted. Imperfections from the far-field boundaries are expected to deteriorate the theoretical convergence order of four. Clearly, the absolute errors for IB method lie above those for the exact boundary treatment, but the convergence among both treatments is consistent. Note that the solution obtained with n = 10e is also only an approximation to the exact analytical solution.

Figure 2 Instantaneous perturbation pressure fields at t = 250 for n = 7. Left half corresponds to exact boundary treatment and right half to the IB method.

Figure 3 Acoustic pressure along the output line, shown in Figure 1 for varying mesh refinement levels n = 7 to n = 10, where level n = 10e denotes the solution at level n = 10 with an exact boundary formulation.

Figure 4 L 2-norm of the error of the acoustic pressure on the output line over cell length Δ for the exact boundary formulation (exact) and the IB method (IBM).

3.2 Sound generation of a co-rotating vortex pair

The two direct-hybrid CFD/CAA methods are applied to a case for which an analytic solution is available. The acoustic field generated by a co-rotating vortex pair is computed and the results are compared to each other with respect to the accuracy and the solution effort. In this two-dimensional test case, shown in Figure 5, two vortices with identical circulation $\mathrm{\Gamma }$ are placed in a distance of 2r ₀ and move along a circular path with radius r ₀. The induced rotational Mach number is $M_r=\mathrm{\Gamma }/\left(4\pi r_0{c}_s\right)$. The period of rotation is $T_r=8{\pi }^2{r}_0^2/\mathrm{\Gamma }$ corresponding to an angular frequency of ${\omega }_r=\mathrm{\Gamma }/\left(4\pi r_0^2\right)$. For an inviscid and incompressible flow field Müller and Obermeier [29] determined an analytical solution for the excited acoustic field, which is given by a complex potential function. The expression for the normalized perturbed pressure field then reads$$\frac{p^{\mathrm{\prime }}\left(r,\theta ,t\right)}{{\rho }_0{c}_s^2}=-\frac{1}{64{\pi }^3}{\left(\frac{\mathrm{\Gamma }}{r_0{c}_s}\right)}^4\left[J_2\left(\mathrm{kr}\right)\sin \left(\mathrm{\Psi }\right)+Y_2\left(\mathrm{kr}\right)\cos \left(\mathrm{\Psi }\right)\right],$$(13)with the wave number $k=2{\omega }_r/c_s$ and the Bessel functions of the first and second kind J ₂ and Y ₂. The time dependent angle $\mathrm{\Psi }=2(\theta -{\omega }_{r}t)$ shows that the acoustic frequency is twice the angular frequency of the vortices. To accurately represent the vortices without singularity in the vortex center, the vortex model proposed by Scully [30] is used$$u_{\theta }\left(r^{\mathrm{\prime }}\right)=\frac{\mathrm{\Gamma }}{2\pi }\frac{r^{\mathrm{\prime }}}{\left(r_c^2+r^{\mathrm{\prime }2}\right)},$$(14)in which u _θ is the circumferential velocity, $r\mathrm{\prime }$ is the distance from the vortex origin, and r _c is the radius of a specified vortex core region (Fig. 6). Thus, the initial flow field is given by$$u\left(r,\theta \right)=-\frac{\mathrm{\Gamma }}{2\pi }\left[\frac{r_{+}}{\left(r_c^2+r_{+}^2\right)}\sin \left({\theta }_{+}\right)+\frac{r_{-}}{\left(r_c^2+r_{-}^2\right)}\sin \left({\theta }_{-}\right)\right],$$ $$v(r,\theta )=\frac{\mathrm{\Gamma }}{2\pi }\left[\frac{r_{+}}{(r_c^2+r_{+}^2)}\cos \left({\theta }_{+}\right)+\frac{r_{-}}{(r_c^2+r_{-}^2)}\cos \left({\theta }_{-}\right)\right],$$where the subscripts + and − indicate the two vortices. The rotational Mach number is M _r  = 0.1 and the vortex-core radius is r _c  = 0.2r ₀. The computational domain for the CFD is $(x,y)\in [-64r_0,64r_0{]}^2$ and for the CAA $(x,y)\in [-128r_0,128r_0{]}^2$ The grid spacing is determined by the grid refinement level n via Δ _n  = Δ₀/2 ⁿ with the definition of Δ₀ = 256r ₀. The grid of the CAA domain is uniformly refined at level n = 9, whereas the CFD grid features four circular refinement patches around the origin with the radii r/r ₀ = 15, 10, 5, and 4 consecutively increasing the grid refinement level from n = 10 to n = 14. The CFD and CAA domains are periodic in z-direction with an extension of two cell length of the coarsest refinement level applied, i.e., z ∈ [−Δ₉, Δ₉]. The acoustic source terms are determined in a circular region within the radius r/r ₀ = 3.8 around the origin. A polynomial degree of two is used in the DG solver such that the acoustic wave is resolved by 18 PPW.

Figure 5 Sketch of the co-rotating vortex pair. The vortices are moving along a circular path with radius r 0.

Figure 6 Comparison of the tangential velocity u θ around a point vortex between the analytical and the vortex-core model proposed by Scully [30].

Contours of the perturbed pressure field obtained by the LB-DG method are depicted in Figure 7 showing the typical double spiral pattern. In Figure 8 the distribution of the analytical perturbed pressure (eq. (13)) along y = x is compared to the results obtained from the FV-DG and the LB-DG simulations. The amplitude is scaled by the pressure value of the first peak $p_{\max }^{\text{'}}$. The scaling accounts for the reduced energy introduction into the acoustic system when comparing a vortex core model to an analytically defined vortex. The numerical results of the FV-DG method match very well with the analytical solution. The amplitude obtained by the LB-DG method slightly over-predicts the maximum of the acoustic pressure.

Figure 7 Contours of the perturbed pressure scaled by the magnitude of the first peak $p_{\max }^{\text{'}}$ of the co-rotating vortex setup.

Figure 8 Perturbation pressure distribution over the radial distance r along the line y = x.

3.3 Flow around a circular cylinder

The flow and acoustic fields around a circular cylinder are predicted for a free-stream Mach number of M = 0.3 and a Reynolds number of Re _D  = 200 based on the cylinder diameter D. The CFD and CAA domains are defined by $(x,y)\in [-10D,30D]\times [-10D,10D]$ and $(x,y)\in [-80D,80D]\times [-80D,80D]$. The base grid refinement level in the CFD domain is n = 10 with Δ₀ = 160D. Four refinement patches decreasing the grid spacing consecutively up to a level of n = 14 as shown in Figure 9. The domain of the CAA solver is discretized on a mesh which is one level coarser than the CFD mesh and where the three finest patches of refinement are omitted. That is, the grid starts from a level n = 9 and features a region of level n = 10 around the cylinder and its near wake region. The vortex shedding with a period of $T{c}_{\mathrm{\infty }}/D=17$ [10] excites acoustic waves with a wave length of λ/D = 17. The acoustic source terms are calculated in the region $(x,y)\in [-2.5D,25D]\times [-4.0D,4.0D]$. To reduce the generation of spurious noise at the boundary of the source term region, the source terms are smoothly decreased by applying a cosine shaped filter function. The filtering is applied in a layer of 1D upstream, 15D downstream, and 1.5D at the top and bottom of the source region, as depicted in Figure 9. The current grid spacing and polynomial degree of three results in a resolution of the acoustic source term region by 19 points per cylinder diameter in the CAA. In Figure 10 a snapshot of the perturbed Lamb vector and the instantaneous perturbed pressure field are shown at the same time level. It can be seen that the source terms are smoothly reduced towards the downstream end of the source term region by the cosine filter. The minimum acoustic wave length in front of the cylinder is ${\lambda }_{\min }/D=11.9$, which is resolved by 152 PPW. Scattering and reflection effects at the solid surface are clearly observed upstream and downstream of the cylinder.

Figure 9 Computational grid (bottom) and Q-criterion (top) to visualize vortical structures are shown for the flow over the circular cylinder. The region to determine the acoustic sources is defined by the light red line (). The source terms are smoothed in the region between the light and dark red box () by a cosine shaped filter function.

Figure 10 Contours of the x- and y-component of the perturbed Lamb vector $\mathbf{L}\mathrm{\prime }$ (left) and the perturbed pressure (right). The source term region is defined by the light red lines (). The distribution along the white line (right) is shown in Figure 11.

The distribution of the perturbed pressure along the line x/D = 0 outside of the source term region is shown in Figure 11 for the FV-DG and the LB-DG methods and a direct numerical simulation (DNS) of the acoustic field from [10]. In the DNS the perturbed pressure signal is directly determined from the flow field by $p\mathrm{\prime }=p-\bar{p}$. Taking the DNS as the reference solution, the wavelength is slightly underpredicted by the FV-DG solution and the LB-DG distribution matches better. The amplitude is slightly underpredicted by the FV-DG method, while the LB-DG amplitude fits for the negative peaks and is overpredicted at the positive peaks. All deviations decrease significantly at increasing distance from the cylinder. A similar deviation of the perturbed pressure field as observed between the FV-DG and the DNS reference results is reported for a hybrid CFD/CAA simulation [10], in which the APE-4 system for the acoustic and the Navier–Stokes equations for the flow field are solved on a body-fitted O-type grid with a similar grid spacing in circumferential direction on the cylinder’s surface as in the present setup. For the identical CAA setups, the discrepancy in the acoustic pressure fields results from different predictions of the vortex noise sources by the LB and the FV solver. The better agreement of the LB-DG solution may originate from the different lower dissipation error of the LB method even compared to Navier–Stokes high order schemes [31].

Figure 11 Distribution of the perturbed pressure of the flow over a circular cylinder at $\mathrm{R}{\mathrm{e}}_D=200$ over the vertical line x/D = 0.

4 Conclusion

In a direct-hybrid CFD/CAA simulation framework two different CFD methods were coupled to a high-order discontinuous Galerkin spectral element CAA method to determine aeroacoustic noise. In the CAA, the APE have been modified to include an IB formulation to mimic a solid wall, based on porous media modeling in the limit of vanishing porosity and permeability. The results for the validation case of the diffraction of acoustic waves interacting with a plate showed that this boundary formulation reproduces the analytical solution with good agreement. For decreasing grid spacing the L ₂-norm of the error of the acoustic pressure shows a similar slope for the error reduction as an exact boundary formulation, i.e., for the considered case approximately third-order accuracy.

The two CFD methods, i.e., a FV method for the compressible Navier–Stokes equations and a LB method, were coupled to the same DG method and applied to two benchmarks cases. The results of both direct-hybrid methods for the acoustic pressure field of a co-rotating vortex pair showed likewise good agreement with the analytical solution. For the low Reynolds number flow around a circular cylinder the perturbed pressure field is well estimated by the combination of the proposed immersed boundary formulation of the acoustic perturbation equation and the two hybrid CFD/CAA methods in comparison with a DNS reference solution. A comparable solution quality is obtained by both coupled methods confirming the findings from the spinning vortex pair.

Since the LB method is a one-stage method in contrast to the multi-stage Runge–Kutta method used for the FV and DG solver, the communication needed for the parallelized DG method cannot be easily hidden by an interleaved execution pattern for the LB method as it can be done for the FV method. Therefore, a general conclusion on the performance of the two coupled solvers is difficult to draw, since the parallel efficiency obtained for the direct-hybrid methods depends to a large extent on the implementation of the solvers, the problem size, its domain partitioning, and the number of computing cores involved.

Conflict of interest

The authors declare no conflict of interest.

Acknowledgments

The funding of the European Union’s Horizon 2020 research and innovation programme within the project INVENTOR (innovative design of installed airframe components for aircraft noise reduction) listed under the grant agreement No. 860538 and of the Deutsche Forschungsgemeinschaft (DFG) grant number SCHR 309/54 is gratefully acknowledged. The authors also gratefully acknowledge the Gauss Centre for Supercomputing e.V. (www.gauss-centre.eu) for funding this project by providing computing time on the GCS Supercomputer HAWK at Höchstleistungsrechenzentrum Stuttgart (www.hlrs.de).

References

All Tables

All Figures

	Figure 1 Sketch of the computational domain. Data is sampled along the line labeled “output”.
In the text

	Figure 2 Instantaneous perturbation pressure fields at t = 250 for n = 7. Left half corresponds to exact boundary treatment and right half to the IB method.
In the text

	Figure 3 Acoustic pressure along the output line, shown in Figure 1 for varying mesh refinement levels n = 7 to n = 10, where level n = 10e denotes the solution at level n = 10 with an exact boundary formulation.
In the text

	Figure 4 L 2-norm of the error of the acoustic pressure on the output line over cell length Δ for the exact boundary formulation (exact) and the IB method (IBM).
In the text

	Figure 5 Sketch of the co-rotating vortex pair. The vortices are moving along a circular path with radius r 0.
In the text

	Figure 6 Comparison of the tangential velocity u θ around a point vortex between the analytical and the vortex-core model proposed by Scully [30].
In the text

	Figure 7 Contours of the perturbed pressure scaled by the magnitude of the first peak $p_{\max }^{\text{'}}$ of the co-rotating vortex setup.
In the text

	Figure 8 Perturbation pressure distribution over the radial distance r along the line y = x.
In the text

	Figure 9 Computational grid (bottom) and Q-criterion (top) to visualize vortical structures are shown for the flow over the circular cylinder. The region to determine the acoustic sources is defined by the light red line (). The source terms are smoothed in the region between the light and dark red box () by a cosine shaped filter function.
In the text

	Figure 10 Contours of the x- and y-component of the perturbed Lamb vector $\mathbf{L}\mathrm{\prime }$ (left) and the perturbed pressure (right). The source term region is defined by the light red lines (). The distribution along the white line (right) is shown in Figure 11.
In the text

	Figure 11 Distribution of the perturbed pressure of the flow over a circular cylinder at $\mathrm{R}{\mathrm{e}}_D=200$ over the vertical line x/D = 0.
In the text