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

1 Introduction

In computational aeroacoustics (CAA), hybrid simulation approaches are state of the art [1]. The advantage of the classical hybrid simulation approach is that the flow field and the acoustic propagation are simulated independently. First, the flow field is simulated, and acoustic sources are extracted. Subsequently, the acoustic sources are used in an acoustic propagation simulation. In doing so, the challenging multi-scale character of aeroacoustics is tackled by separating the scales of hydrodynamics and acoustics. Schoder and Kaltenbacher [2] present a overview of the recent developments in hybrid modeling of aeroacoustic sound.

Depending on the degree of modeling of the acoustic sources, steady or unsteady Reynolds-averaged Navier–Stokes (RANS) equations, detached eddy simulation (DES), or large-eddy simulations (LES) are used to simulate the flow field. Using the steady RANS equations, acoustic sources can be modeled, e.g., with the stochastic noise generation and radiation (SNGR) model of Bechara et al. [3]. Unsteady simulations like DES or LES promise to give more accurate acoustic sources, but the required computational effort increases considerably. Here, the acoustic sources are extracted from the unsteady flow field.

In complex flow regimes like separation, high-fidelity methods, such as wall resolved LES, are beneficial. With increased available computing power, these methods become more and more relevant for industrial CAA application. However, wall resolved LES of large configurations are often still computational not feasible. Zonal RANS–LES approaches are used more often e.g., by Zhangg et al. [4], Erbig and Maihöfer [5], and Kuhn et al. [6] and are also used in the CAA community. Terracol [7], Bernicke et al. [8], and Satcunanathan et al. [9] use a zonal LES approach to predict the trailing edge noise of airfoils. The idea is to extract the acoustic sources inside the zonal domain and apply hybrid CAA methods to predict the far-field noise.

Traditionally, in hybrid CAA of low Mach number cases, the acoustic sources are extracted from an incompressible flow simulation. For example, Nusser and Becker [10] use this approach to predict vehicle interior noise. In such a case, specific aeroacoustic effects like acoustic feedback cannot be predicted due to a neglected interaction of hydrodynamics and acoustics. However, acoustic feedback can be relevant in developing technical products. Solving the compressible Navier–Stokes equations (NSE) can capture these effects in a direct noise computation (DNC), with the downside of high computational effort. An overview of DNC can be found in Bailly et al. [11]. Using DNC, Frank and Munz [12] predicted aeroacoustic feedback on a side-view mirror, and Kuhn et al. [6] and Erbig and Maihöfer [5] predicted Rossiter feedback on cavity configurations relevant for the automotive industry.

When applying the zonal approach in combination with DNC, the acoustic far-field propagation outside this domain cannot be calculated due to the spatial limitation of the domain. In this case, the hybrid CAA approach is suitable to predict the far-field noise. Moreover, one obtains a solution that describes the acoustics. This is advantageous since the solution of the compressible NSE includes both the hydrodynamic and acoustic solution, and the separation of both is costly and still subject to active research, see Schoder et al. [13].

We extend the discontinuous Galerkin framework FLEXI¹ to a zonal direct-hybrid CAA framework within this work. Here, the flow field and the acoustic propagation are solved with the same underlying discontinuous Galerkin scheme. Within the framework, the compressible NSE are solved to predict the flow field in a zonal manner. Solving the compressible NSE allows depicting aeroacoustic feedback in the flow simulation using DNC. Here, high-order schemes, such as discontinuous Galerkin, are well suited due to the beneficial dissipation and dispersion properties, see Flad et al. [14]. The chosen LES approach to determine the acoustic sources is the following: In the first, step a RANS simulation of the whole domain is simulated. Subsequently, a restricted area of interest is defined where a zonal LES using the data provided by the RANS simulation is performed. Unlike DES, the zonal LES are carried out wall resolved as an implicit LES without an explicit LES model. At the RANS–LES interface, we use the recycling–rescaling anisotropic linear forcing (RRALF) introduced by Kuhn et al. [15] to generate the inflow turbulence. This turbulent inflow method requires the time-averaged mean velocities and the full Reynolds stress tensor in the inflow region as input data. Usually, these data are not known beforehand. We rely on data obtained from a RANS simulation as input data. However, in commonly used turbulence models, the distribution of the Reynolds stresses is not modeled. Therefore, we developed a modeling approach to model the full Reynolds stress tensor based on a RANS simulation.

The framework offers two possibilities to perform CAA simulations. First, the flow simulation is performed independently, and the extracted acoustic sources are stored. Here, the necessary I/O operations quickly become the bottleneck, especially in large-scale simulations done in high-performance computing. Second, the flow and acoustic propagation simulations are performed simultaneously in large-scale simulations, and the acoustic sources are communicated via the messaging passing interface (MPI). This reduces the amount of data stored and thus makes the simulations easier to handle. The bottleneck of the I/O operations was also identified in the work of Schlottke-Lakemper [16]. He used a similar direct-hybrid approach. The imbalance between the acoustic and flow solver is taken into account by a static load balancing.

The new proposed model of the Reynolds stresses is validated by comparing results of a zonal LES based on RANS simulation data and a zonal LES based on wall resolved LES data. Here, we chose the example of turbulent flow over a flat plate. The direct-hybrid acoustic framework is validated against the DNC of a NACA 0012 airfoil. The NACA 0012 airfoil at Re = 10⁵ is chosen since it is known to exhibit tonal self-noise. We show that the direct-hybrid method can reproduce tonal self-noise in the acoustic propagation simulation. The applicability of the new model to a zonal direct-hybrid simulation is demonstrated using the zonal simulation of a NACA 64418 trailing edge at Re = 10⁶. This example is used to validate the whole framework. First, the zonal hydrodynamic flow field is validated against a wall resolved LES of the entire airfoil. Second, the acoustic propagation simulation results are compared to the full airfoil’s DNC and the zonal LES’s DNC.

2 Numerical method

High-order methods are beneficial when dealing with turbulence and acoustics due to their low dissipation and dispersion errors [17]. In this work, we use the discontinuous Galerkin method, implemented in our framework FLEXI [18]. This framework has been successfully applied to DNC simulations [6, 12] in the past. The extension towards a zonal direct-hybrid framework is described in the following.

2.1 Discontinuous Galerkin spectral element method

Within this work, we consider the compressible NSE, which intrinsically include hydrodynamics and acoustics, and the APE-4 equation system proposed by Ewert and Schröder [19] for acoustic propagation. Both systems of equations are solved with the discontinuous Galerkin spectral element method (DGSEM).

To numerically solve the system of equations, the physical domain is discretized with three-dimensional, non-overlapping hexahedral elements. For a better representation of the geometry, we use curved elements, and for flexible meshing, we allow an unstructured mesh topology. Each element is transformed from the physical space to the unit reference element E = [−1; 1]³, where represents the reference coordinates. The transformation follows Kopriva [20]. The variational form is obtained by multiplying a basis function by a test function ϕ, which is chosen in the Galerkin approach identical to the basis function. Integrating over the reference element E and integration by parts yields the weak formulation of the DGSEM:

(1)

where and denote the numerical flux function normal to the surface for the inviscid and the viscous term, respectively. In the case of the NSE the volume flux is and in the inviscid APE-4 case the viscous flux vanishes, and the volume flux is .

As basis function, we choose a tensor product ψ^N of 1-D Lagrange polynomials of degree N,

(2)

where is the solution vector and the time-dependent nodal degree of freedom. Using the (N + 1)³ Legendre–Gauss–Lobatto quadrature points for integration and interpolation yields an efficient collocation approach following Kopriva [21].

The Roe and the Lax–Friedrichs Riemann flux [22] is used to determine the inviscid surface flux at the cell interface in the case of the NS and the APE-4 equation system. The viscous flux of the NSE system is approximated by the lifting procedure by Bassi and Rebay [23]. To stabilize the DG scheme in the NS case, the kinetic energy preserving two-point split flux form proposed by Gassner et al. [24] and implemented in the DG framework by Flad and Gassner [25] is used. It is based on the kinetic energy preserving flux introduced by Pirozzoli [26]. As time integration scheme, the low storage fourth-order explicit Runge–Kutta method of Carpenter and Kennedy [27] is applied.

The treatment of the boundaries is crucial when dealing with acoustics. To reduce artificial reflections of acoustic waves and the turbulent flow structures we use a sponge zone proposed by Pruett et al. [28] and boundary conditions of Dirichlet type in weak form. Further details about the acoustic properties of the framework FLEXI can be found in Flad et. al [14].

2.2 Turbulent inflow

In order to perform a zonal LES, a turbulent inflow method is necessary to generate the required turbulence in the inflow region. In this work, we apply the recycling–rescaling anisotropic linear forcing introduced by Kuhn et al. [15]. This method is a combination of an anisotropic linear forcing (ALF) introduced by de Laage de Meux et al. [29] and a recycling rescaling method introduced by Lund et al. [30]. Recycling the fluctuating part of the solution at the interface in combination with the volumetric, anisotropic linear forcing produces good quality inflow turbulence. Furthermore, the restriction of self-similarity by the recycling rescaling methods is weakened by this combination. In the case of non-self-similar problems the ALF corrects the introduced error at the inflow.

The combination of these two methods gives us a flexible turbulent inflow method. Well established turbulent inflow methods like digital filtering by Klein et al. [31] require target statistics, namely the mean flow, the Reynolds stresses, and length scales. The quality of the inflow turbulence and the relaxation length towards natural turbulence after the inflow directly depend on the quality of the turbulent statistics. Here, the RRALF method offers better flexibility. It is not required to provide the Reynolds stresses used by the ALF in the entire inflow region since the ALF is not required to be used over the entire inflow region. It can be partially applied in the inflow region by using it primarily as a correcting mechanism. Here one would choose a region of good quality target data e.g., measurements or another simulation. At the same time, the recycling of the turbulent fluctuations ensures high-quality turbulence over the whole inlet, even in the near wall region.

Similar to the original recycling rescaling, the scaling function is split into an inner and outer scaling. Here, the scaling is done by a static scaling of the recycling plane. It is based on the inflow plane node distribution y_in and the provided target data. The recycling plane is build up in the following way:

(3)

where is the time-averaged target velocity gradient at the wall of the inlet and at the recycling plane. The dimensionless wall distance reads as and is a smooth blending function defined as,

(4)

To prevent a drift of the solution only the fluctuating part of the solution is recycled, according to Sagaut et al. [32]. At the inflow plane the instantaneous solution is reconstructed with the provided fixed mean target profile. The reconstruction follows:

(5)

where Δz accounts for an optional spanwise shift in order to reduce the periodicity of the recycling.

The volume forcing follows de Laage de Meux et al. [29]. It is a volumetric forcing of the velocity field towards a prescribed mean velocity profile as well as a prescribed Reynolds stress distribution . The volume force f_i = A_iju_j + b_i, added to the momentum equation, results from solving the following linear system of equations for each time step. Equations (6) and (7) are solved sequentially for A_ij and b_i. Equation (6) is solved using Crout’s method:

(6)

(7)

The estimates of the time-averaged mean velocities and the Reynolds stresses are computed using an exponential time filter which is updated in every time step. Due to the recycling, the turbulent fluctuations can overshoot the prescribed target data at the inflow. The original ALF was modified by introducing an explicit damping that corrects the overshoot rapidly. This damping acts on the momentum equation in the spatial direction where the target data is overshot.

To set up the turbulent inflow region, the recycling plane is typically placed 5–10 boundary layer thicknesses δ_99,in inside the fluid domain. The region where the ALF is applied ends about 2–3 δ_99,in in the streamwise direction in front of the recycling plane to allow for natural relaxation of the turbulence. The overall performance impact of this method depends on the simulated case. Typically, there is a performance overhead of 5–20%. Computing the volumetric forcing term inside the ALF region and evaluating the turbulent statistics at the recycling plane are efficient point local operations and only partially contribute to the overhead. Despite hiding the communication of the turbulent statistics from recycling plane to the inflow plane, using non-blocking communication, the recycling shows the greatest performance impact.

2.3 Target data generation

The ALF and the determination of the scaling of the recycling plane rely on provided target data, namely the time-averaged mean velocities and the Reynolds stress distribution in the turbulent inflow region. In this work, we use data provided by a simulation solving the RANS equations using Menter’s Shear Stress Transport turbulence model. The RANS simulations are computed with the open-source CFD toolbox OpenFOAM. This toolbox uses a finite volume method. We use a generalized coupling approach to couple the finite volume solution to the discontinuous Galerkin framework. Based on scattered data, we construct a Delaunay triangulation and reconstruct the solution on the nodal points of the high-order mesh. Using such a triangulation, we are independent of the data structure of the RANS solver, and it is easily interchangeable.

The time-averaged mean velocities are obtained directly from the RANS simulation results. The distribution of each single Reynolds stresses is not modeled in the turbulence model, only the turbulent kinetic energy k and the turbulent dissipation rate ω. The obvious option is to use a Reynolds stress transport model, but these are not very commonly used. Therefore, we propose a simple model approach to approximate the full Reynolds stress tensor based on the distribution of the turbulent kinetic energy k. The approach is based on the findings of Bradshaw et al. [33]. In this one equation turbulence model, the shear stress is modeled by a constant scalar scaling of the turbulent kinetic energy k. This model is very accurate between the viscous sublayer and the outer layer. We extended this approach to the full Reynolds stress tensor. Therefore, instead of a single scalar, we define a symmetric scaling tensor S_ij with an individual scaling scalar dedicated for each Reynolds stress component. Here, the trace of tensor S_ij must be exactly two, to conserve the turbulent kinetic energy k. To generalize the model, the sign of the Reynolds stresses approximated by Boussinesq’s assumption is used to specify the appropriate sign of the Reynolds stresses. The model is described in Equation (8):

(8)

In Figure 1 the Reynolds stresses normalized by the turbulent kinetic energy k of a wall resolved LES is plotted in the case of a turbulent flow over a flat plate. Figure 1 indicates that the assumption of a constant scaling factor is valid for the outer region of the boundary layer. These factors are also only slightly varying for different attached flow configurations.

Figure 1 Reynolds stresses normalized by the local turbulent kinetic energy k of a turbulent boundary layer.

Figures 4 and 10 display the modeled Reynolds stresses compared to a wall resolved LES at the same position of the cases considered in this work, namely the turbulent flow over a flat plate and a NACA 64418 airfoil. As expected, the modeled Reynolds stresses match the reference well in the boundary layer’s outer region. To handle the inaccuracy at the wall-near region, we follow the strategy that we skip the wall-near region when applying the ALF. Typically we apply the ALF to the outer 2/3 of the boundary layer. Therefore, we skip the region with the largest modeling errors and rely on the natural development of turbulence in the near-wall region. This approach is applicable due to the recycling of the RRALF method, which allows for rapid natural development of near-wall turbulence.

2.4 Direct-hybrid simulation framework

Within the scope of this work, the discontinuous Galerkin framework FLEXI was extended towards a direct-hybrid acoustic simulation framework. Besides solving the compressible NSE, the APE-4 equation system proposed by Ewert and Schröder [19] is now supported within this framework. In the vortex noise dominated problems at low Mach number considered in this work, we use the instantaneous perturbed Lamb vector L′ = [ω × u]′ as the acoustic source, where ω is the vorticity and u the velocity vector. In the region where acoustic production is assumed, e.g., at the trailing edge of an airfoil, the acoustic source is extracted from a zonal LES solving the compressible NSE. This is done by masking beforehand the region where acoustic source terms are extracted. This offers the flexibility to reduce the data required to communicate as well as to skip specific acoustic source regions.

There are two possibilities to perform the hybrid acoustic simulation implemented within this framework. First, the simulation to determine the acoustic sources is carried out independently and the acoustic sources are extracted and stored. Here, the acoustic sources do not have to be stored at every explicit time step to achieve good results. The acoustic sources are interpolated to the required time step in the acoustic propagation simulation. Such a hybrid simulation approach requires a lot of storage space and the required I/O operations are getting the bottleneck of such simulations.

The second approach is to compute the flow simulation and the acoustic propagation simulation simultaneously in a direct-hybrid manner. This offers the advantage of directly communicating the extracted acoustic sources to the acoustic solver. This approach omits to save the acoustic sources. The resulting reduction of required I/O operations makes this approach well suited for large-scale simulations done in high-performance computing. Further, an exchange of the acoustic sources at each time step or even at each Runge–Kutta stage is possible without significant impact on the performance compared to the approach of storing the acoustic sources. This also prevents the uncertainty of a negative impact due to the temporal resolution of the source term data. Both solvers within the framework are explicit solver. Therefore, within each timestep, both solver synchronize the minimal timestep. Thereby, and due to the larger computational effort to solve the NSE, a static load balancing is introduced between both solvers. Solving the NSE is approximately twice as expensive as solving the APE-4 equation system within this DG framework. All communications between the two solvers are implemented in a non-blocking manner allowing for a staggered execution and results in an improved performance. Despite the synchronization in each timestep, it is possible to exchange the acoustic sources each nth timestep. Here, an improvement in performance is noticeable up to synchronizing each 10th timestep. Beyond that, a performance improvement is not noticeable. In such a scenario, the acoustic solver buffers the source data of four time instances and performs a third-order interpolation in time. All communications between the two solvers are implemented in a non-blocking manner which allows for a staggered execution and results in improved performance. At the moment, there is still a restriction to use the same mesh for both solvers. To capture the convection and refraction in the acoustic simulation and for computing the perturbed Lamb vector, we need to provide mean flow data. This can be done by time-averaging a precursor simulation or, in case of hybrid simulations without precursor simulation, e.g., in the direct-hybrid case, by computing a moving average in the flow solver and communicating it to the acoustic solver.

The advantage of the direct-hybrid approach is relevant for high-performance computing. This approach is beneficial in the case of large-scale simulations with only few executions. Here, we do not greatly profit from the possible reusability of the acoustic source date by storing them. Also, large-scale simulations likely produce a huge amount of data to be handled and stored. In the case of high reusability of the acoustic sources, e.g., when there is no significant impact of changes in the geometry or flow condition on the acoustic sources, the traditional approach is the more economical approach.

Comparing the performance of both implementations at the example of writing each 10th timestep shows an overall improvement of 40% for the direct-hybrid approach, which highly depends on the available hardware and I/O performance. More frequent exchange of the acoustic source date favours the second approach where as a less frequent favours the first approach. In this work, the second approach is chosen.

3 Simulation results

The above described direct-hybrid simulation framework consists of two major building blocks. First, the turbulent inflow method is used for the zonal LES approach based on the RRALF method. For coupling with a RANS simulation, the model approach described in Section 2.3 was introduced. Second, the hybrid acoustic simulation framework itself. It consists of extracting acoustic source terms, the communication of acoustic sources to the acoustic propagation solver, and the acoustic propagation solver itself.

To validate the turbulent inflow based on the new model approach, in Section 3.1 a turbulent flow over a flat plate is considered. Here, a zonal LES using the modeled target data is validated against a zonal LES based on data from an extended tripped LES.

The second building block is validated using two test cases. First, in Section 3.2 the hybrid simulation framework is validated by simulation of a NACA 0012 airfoil. The acoustic propagation simulation results are compared to the results of DNC. Afterward, in Section 3.3 results of a direct-hybrid zonal trailing edge simulation of a NACA 64418 airfoil are presented. Here, the complete zonal direct-hybrid framework is validated against target data obtained from a compressible wall resolved LES of the whole airfoil. Also the hybrid acoustic results are compared to the DNC results within the zonal region.

3.1 Turbulent flow over a flat plate

Three simulations of a weakly compressible turbulent flat plate were carried out with the same numerical scheme and mesh resolution. First, a tripped boundary layer simulation was conducted, in the following referred to as full LES. Second, these simulation results were used as target data for a zonal LES of the rear half of the flat plate (ZLES I). Third, a second zonal LES of the rear half of the flat plate using RANS simulation data as target data (ZLES II). In this case, a two-dimensional RANS with Menter’s SST turbulence model was carried out. Based on the RANS simulation, the Reynolds stresses were modeled using the model proposed in Section 2.3.

The tripped boundary layer used the incompressible, turbulent velocity profile provided by Eitel-Amor et al. [34] at Re_θ = 750 and a Mach number of Ma = 0.3 as inflow data. The quadratic trip is placed one boundary layer thickness behind the inflow and has the size of y⁺ = 50 in wall units referring to the inflow profile. The domain size is 205δ_99,in × 12δ_99,in × 16δ_99,in, δ_99,in being the boundary layer thickness at the inflow of the tripped case. The domain resolves the boundary layer up to about Re_θ = 2800.

The mesh sketched in Figure 2 comprises 459 150 elements, which results in about 235 million DOF for a polynomial degree of N = 7. The mesh resolution is improving towards higher Re_θ, from x⁺ = 25–20.5, y⁺ = 1–0.8 and z⁺ = 12–9.8. The simulation was carried out for 531 convective time units T* = δ_99,z/U_∞ and was averaged over the last 115T* to obtain the first- and second-order turbulence statistics. Here, U_∞ denotes the free-stream velocity and δ_99,z the boundary layer thickness at the inflow of the zonal domain at Re_θ = 1800 and is from now on the reference length scale. The simulation was carried out on 256 Nodes with 128 cores each, and computational cost of the simulation was about 445 000 CPUh, including I/O.

Figure 2 Sketch of the mesh used for the tripped and zonal simulations of the turbulent flow over a flat plate. Zonal region reaching from Reθ = 1800–2800.

The two zonal LES had the same mesh resolution as the full LES and resolved the boundary layer from Re_θ = 1800–2800. In Figure 2 the zonal region is highlighted. Hereby, the mesh was reduced to about 112 million DOF. The ALF region starts at the inflow and extends over 6δ_99,z in streamwise direction. In the case of ZLES I, the first wall-normal high-order cell was skipped. In ZLES II the 14 wall-normal cells were skipped. The recycling plane was placed at 8 δ_99,z in streamwise direction.

At the inflow of ZLES II the mean velocity profile of the RANS simulation is enforced. In Figure 3, the mean velocity profile at the inlet of ZLES II and the results of the full LES are plotted. The agreement with the full LES is good in the viscous sublayer, the outer part of the log-law region, and the outer layer. There is a noticeable deviation from the full LES at the buffer layer. In Figure 4 the Reynolds stresses of the full LES, ZLES II and the modelled Reynolds stresses are plotted at the inflow. The good agreement between the full LES and ZLES II demonstrates the advantage of the recycling approach. As stated by Bradshaw et al. [33] the turbulent shear stress is captured well over the whole boundary layer. The normal stresses show good agreement with the full LES in the outer region of the boundary layer. Here, the model overestimates the Reynolds stresses slightly. In the wall-near region, the model is not capable of reproducing the behavior of the full LES. This is due to the constant scalar scaling within this model approach. As described in Section 2.3 the wall-near region is skipped in the ALF forcing. The region where the ALF forcing is active is indicated in blue.

Figure 3 Comparison of the mean velocity profile between the full LES (tripped), the RANS simulation, and ZLES II (zonal LES based on the RANS simulation data) at the inflow of the zonal LES. The mean velocity profile of the RANS simulation and the ZLES are identical.

Figure 4 Comparison of the Reynolds stresses between the full LES (tripped), the RANS simulation (modeled), and zonal LES II (based on the RANS data) at the inflow of the zonal LES. The used values for the scaling tensor Sij are S11 = 0.9, S22 = 0.45, S33 = 0.65 and S12 = 0.3. The area where the ALF is applied is indicated in blue.

ZLES I was carried out for 277T* and averaged over 69T* and ZLES II was carried out for 238T* and averaged over 66T* to obtain the turbulence statistics. The zonal simulation was carried out on 128 nodes with 128 cores each, and the computational cost was about 227 000 CPUh, including I/O. Figures 5 and 6 depict the first- and second-order turbulence statistics at Re_θ = 2240 and Re_θ = 2536 within the zonal region, for which reference data from Eitel-Amor et al. [34] and Schlatter and Örlü [35] are available. In Figure 5 the time-averaged streamwise velocity of ZLES I and II, the full LES, and reference data from literature is shown in viscous wall units. The agreement of the velocity profiles is good at both Reynolds numbers. The natural development of the boundary layer can equalize the introduced error of the mean velocity field as depicted in Figure 3. The second-order turbulent statistics in form of the normal stresses and the shear stress are illustrated in Figure 6. There is an overall good agreement between ZLES I and II and the references for both Reynolds numbers under investigation. ZLES I is able to reproduce the turbulent statistics of the full LES over the whole boundary layer height very well. However, ZLES I and the full LES tend to underestimate the Reynolds stresses compared to the literature slightly. In the wall-near region, ZLES II reproduces the turbulence well with the exception of the normal stress in the spanwise direction. Here, it is underestimated, especially at the higher Reynolds number. In the outer region of the boundary layer, the Reynolds stresses of ZLES II slightly overestimate the reference. This is a result of the slightly overestimated Reynolds stresses of the model and the pure natural development in the near wall region.

Figure 5 Comparison of the mean velocity profiles of the turbulent boundary layer between ZLES I (target data based on LES data), ZLES II (target data based on a RANS simulation), the full LES (tripped), and reference data from the literature.

Figure 6 Comparison of the Reynolds stresses of the turbulent boundary layer between ZLES I (target data based on LES data), ZLES II (target data based on a RANS simulation), the full LES (tripped), and reference data from the literature.

The zonal simulation of the turbulent flow over a flat plate showed that the turbulent inflow method RRALF is capable to reproduce the solution of a full LES when using the same numerical scheme and the same mesh resolution of the simulation providing the target data. Therefore the deviations of ZLES II to the reference ZLES I can be dedicated to the chosen modeling approach. These deviations are small and the results obtained from the zonal LES based on a two-dimensional RANS simulation as target data are still very good. This capability motivates us to use this approach in zonal hybrid acoustic simulations since it is able to predict turbulence very well, which is essential for hybrid acoustic simulations.

3.2 NACA 0012 Airfoil

A two-dimensional NACA 0012 airfoil is simulated in a direct-hybrid manner to validate the direct-hybrid simulation framework. For the flow and the acoustic propagation simulations, the same mesh of the complete airfoil and the acoustic sources of the whole fluid domain are used. This is done to compare the acoustic solution of both simulations. The airfoil has an angle of attack of 0.5°, the free-stream Mach number is Ma = 0.3 and the Reynolds number is Re = 10⁵. The setup is chosen according to Jones and Sandberg [36]. This configuration is well known for tonal airfoil self-noise. The used mesh consists of 40 934 high-order Elements, the polynomial degree is chosen to N = 5, and a circular sponge region is applied around the airfoil to avoid artificial reflections at the far-field boundaries. This test case was chosen to emphasize the advantages of a compressible simulation as an underlying flow simulation. Here, it is possible to capture tonal self-noise in the acoustic propagation simulation. Due to the interaction of the acoustics and the hydrodynamics in the flow simulation, the hydrodynamic instabilities, excited by the acoustics, form tonal self-noise generating vortices which are extracted as acoustic sources. In Figure 7 the pressure fluctuation is plotted one chord length above the trailing edge. The agreement of the pressure fluctuation between the DNC and the acoustic propagation simulation is good. In the acoustic solver, the amplitude of the pressure fluctuation is slightly underestimated but there is no noticeable phase error. This good agreement is also visible in the acoustic spectra displayed in Figure 8. Here, the acoustic solver is able to reproduce the results of the DNC very well. It is worth mentioning that a tonal frequency at f_tonal = 4.2 is captured in both simulations. Besides the agreement in the tonal frequency, the agreement in the overall spectrum is good again.

Figure 7 Time-dependent dimensionless pressure recorded at an arbitrary time one chord length above the trailing edge. Comparing the DNC results (blue) with the direct-hybrid simulation results (red).

Figure 8 Comparison of the dimensionless PSD of the pressure between the DNC results (blue) and the direct-hybrid simulation results (red).

3.3 NACA 64418 – zonal direct-hybrid simulation

A zonal LES of a NACA 64418 airfoil was carried out. The airfoil has an angle of attack of 6°, which leads to a small separation close to the trailing edge, a free-stream Mach number of Ma = 0.2, and a Reynolds number of Re = 10⁶. As a reference, a tripped LES of the full airfoil was computed according to Krais et al. [18]. The mesh of the simulation consists of 229 620 elements. The simulation was advanced for 15T*, T* being the convective time unit T* = c/u_∞, and averaged over the last 5T* to get the statistics. The computational cost to advance the simulation by one convective time unit was about 43 300 CPUh using 256 nodes.

The setup of the zonal LES is displayed in Figure 9. The simulation setup contains the trailing edge of the airfoil. On the suction side the domain starts at x/c = 0.5 and on the pressure side at x/c = 0.6. The spanwise extend of the domain is 0.05c and is reduced by half compared to the full airfoil. In the background, the mesh, consisting of 58 800 high-order cells, is displayed. The mesh resolution is chosen similar to the reference LES and the polynomial degree is set to N = 7. The resulting mesh resolution at the wall in wall units, normalized by the factor N + 1, is the following: x⁺ declines towards the separation on the suction side from 45 to 5, y⁺ from 3 to 1, and z⁺ from 20 to 2. On the pressure side, x⁺ declines from 45 to 30, y⁺ and z⁺ is in the range of 3–15. Furthermore, in Figure 9 the scaled inflow planes on the suction and pressure side and the corresponding recycling planes are displayed. On the suction and pressure side, the recycling planes are shifted 0.2x/c and 0.11x/c, respectively into the domain which correspondents to 16δ_99,in and 11δ_99,in, δ_99,in being the boundary layer thickness at the corresponding inflow plane. Between the recycling and the inflow planes, the region where the ALF forcing is active is displayed in blue. It extends from the inflow plane to 0.15x/c and 0.08x/c into the domain on the suction and pressure side, respectively which correspondents to 13δ_99,in and 9δ_99,in. Behind the recycling planes, the turbulent vortex structures are visualized by isosurfaces of the Q-criterion, the second invariant of the velocity gradient tensor, of value Q = 200 and colored by the streamwise velocity component. Values of Q > 0 indicate areas dominated by rotation rate compared to the strain rate. Figure 10 displays the modeled target Reynolds stresses within the ALF region on the suction side compared to the full LES. The area where the ALF is applied is indicated in blue.

Figure 9 Trailing edge of the NACA 64418 airfoil. Vortex structure visualized by the Q-criterion (Q = 200) colored by the streamwise velocity in the range [−0.3, 1.3]. The displayed planes indicate the inflow and recycling planes. The forcing region is highlighted in blue.

Figure 10 Comparison of the modelled Reynolds stresses and the results of a wall resolved LES. Plotted in wall-normal direction on the suction side of a NACA 64418 airfoil. The used values for the scaling tensor Sij are S11 = 0.8, S22 = 0.55, S33 = 0.65 and S12 = 0.3. The area where the ALF is applied is indicated in blue.

In a first step, the zonal LES of the trailing edge was simulated and validated against the wall resolved LES of the whole airfoil. The simulation was advanced for 7T* and the flow field was averaged over 4T*. The computational cost to advance the simulation by one convective time unit was about 9300 CPUh. The cost reduces to 20% in the zonal approach. The pressure distribution obtained from the time-averaged flow field is displayed in Figure 11. Here, the pressure distribution of the zonal LES, the full airfoil and the RANS simulation, used for the target data generation, are displayed. Looking at the pressure distribution of the RANS simulation, we see quite good agreement with the result of the full airfoil up to the trailing edge region. As expected from a RANS simulation, the prediction of the separation forming at the trailing edge is not as good as in a wall resolved LES case. Comparing the result of the zonal LES, based on the input of the RANS data, with the results of the full airfoil, we can see an overall good agreement, even at the trailing edge. Therefore, by using the zonal LES approach, based on RANS simulation data in the attached region, we are able to reproduce the separation behavior.

Figure 11 Pressure distribution on the surface of the tripped NACA 64418 airfoil. Comparison of the LES and the RANS simulation of the complete airfoil and the zonal LES based on the RANS simulation data.

In Figure 12 the distribution of the time-averaged mean velocity profiles at four positions on the suction side, pressure side, and the wake are displayed. The positions are displayed in Figure 12 (top right). At the investigated positions, there is good agreement between the zonal LES and the full LES except in the inflow plane. Here, there is a slight deviation to the full LES in the near wall region. In Figure 13 the distribution of the Reynolds stresses is displayed at the same positions. The velocity fluctuations show larger deviations from the reference compared to the mean velocity profiles. Looking at the distribution at the inflow on the suction and pressure side, we can see the largest deviations. Due to the recycling of the fluctuations from the recycling plane to the inflow and the pressure gradient on the suction and pressure side, the shape of the distribution of the velocity fluctuations differs from the reference. In streamwise direction, the distribution trends towards the reference. This is due to the correction by the forcing as well as the natural development of the turbulence. The good agreement of the distributions at the trailing edge leads to a good agreement in the wake. This agreement is quite important because the region at the trailing edge greatly impacts the generation of acoustic emissions.

Figure 12 Mean velocity profiles in comparison to the full LES at the displayed positions. Profiles staggered in streamwise direction. Top left: Suction side. Top right: Mean streamwise velocity, including the positions of the velocity fluctuation profiles. Bottom left: Pressure side. Bottom right: Wake region.

Figure 13 Reynolds stress distributions in comparison to the full LES at the displayed positions. Profiles staggered in streamwise direction. Top left: Suction side. Top right: Mean streamwise velocity, including the positions of the velocity fluctuation profiles. Bottom left: Pressure side. Bottom right: Wake region.

To compute the acoustic spectrum, the zonal LES and the acoustic propagation simulation are advanced for 0.5T* to propagate the acoustic waves in the acoustic solver. Afterward, the pressure signal is recorded for 1.4T* in both simulations. To balance the load, the zonal LES was computed using 42 nodes and the acoustic propagation solver using 22 nodes, 128 CPUs each node. To advance both solver, the simulation of one convective time unit required about 18 100 CPUh. This is approximately twice as expensive as the standalone simulation of the zonal LES. Here, a coupling of the acoustic sources in each Runge–Kutta stage does not allow for an optimal communication hiding, and the higher computational efficiency of the acoustic propagation solver is not fully utilized. To convert the recorded pressure signal into a PSD, a FFT and averaging over three blocks with a Hanning window function applied and 50% overlap is used. In the DNC of the zonal LES as well as the full LES corresponding time intervals are chosen. In Figures 14 and 15 the acoustic spectrum of the direct-hybrid, the DNC of the zonal LES, and the DNC of the full airfoil is plotted.

Figure 14 Dimensionless PSD of the pressure above the trailing edge. Comparison of the direct-hybrid simulation, the zonal DNC, and the DNC of the full LES.

Figure 15 Dimensionless PSD of the pressure below the trailing edge. Comparison of the direct-hybrid simulation, the zonal DNC, and the DNC of the full LES.

The spectrum in Figure 14 is evaluated 0.2 chord length above the trailing edge. Here, we want to focus on the lower frequencies to discuss the influence of the turbulent inflow. At low frequencies, the zonal DNC overpredicts the reference. This is the influence of the turbulent inflow method which generates artificial noise, also partly due to the recycling. As seen in Figure 13, the recycled turbulence deviates from the target data at the inflow. The relaxation causes artificial noise. At those low frequencies, the direct-hybrid simulation is closer to the reference. This is due to the fact that the extraction of the acoustic sources begins after the turbulent inflow region and the artificial noise of the inflow is reduced. This is also one motivation to use this kind of hybrid approach. At frequencies around f = 100 the reference is overpredicted, both by the zonal DNC as well as by the hybrid simulation. In the high frequencies, there is good agreement with the reference. The decline and the amplitude of the acoustic level are comparable, both for the zonal DNC as well as for the direct-hybrid simulation.

In Figure 15 the acoustic spectrum is evaluated 0.1 chord length below the trailing edge focusing on the higher frequencies. We investigate a similar behavior as discussed before. At lower frequencies, we have a dominant influence of the artificial noise, which can be omitted in the hybrid simulation. In the midrange of frequencies, the hybrid simulations over predicts the result compared to the reference. However, in this case, the zonal DNC matches quite well with the reference DNC. Again we have matching results for both three simulations at high frequencies.

The results of the direct-hybrid simulation in combination with a zonal LES prove that this approach is justified. We were able to show that we can reproduce the turbulent statistics on the flow simulation side with a reduced computational expense. Using the hybrid approach, we can filter, at least a dominant part, of the artificial inflow noise generated by the turbulent inflow method, as well as the acoustic propagation to at least half the cost. Further, it allows to simulate the acoustic far-field propagation outside of the zonal domain. The downside of this approach is that it still relies on the quality of the underlying RANS simulation, and there is a reduction in the quality of the acoustic results.

4 Conclusion

A method for direct-hybrid aeroacoustic analyses was implemented in the high-order discontinuous Galerkin framework FLEXI. In large-scale simulations on high-performance computers, the commonly used approach in hybrid aeroacoustic simulations by consecutively simulating the flow field and the acoustic propagation simulation reaches its limits. The storage of the extracted acoustic sources leads to I/O operations that are getting the simulation’s bottleneck. Therefore we follow a direct-hybrid simulation approach, where the flow solver and the acoustic solver are executed parallel and exchange the extracted acoustic sources via the messaging passing interface.

We solve the compressible Navier–Stokes equations (NSE) with an explicit high-order discontinuous-Galerkin scheme to simulate the flow field. The solution of the compressible NSE equations already inherently includes the acoustics and resolves acoustic flow interactions. Due to the high computational cost of such a direct noise computation, our approach is to restrict the highly resolved computation of the flow field to the relevant acoustic source region using a zonal large eddy simulation (LES). In the domain outside of the zonal region, the low fidelity Reynolds-averaged Navier–Stokes (RANS) equations are solved. The turbulent inflow method we use is the recycling rescaling anisotropic linear forcing (RRALF), a combination of a traditional recycling rescaling approach and an anisotropic linear forcing. This method produces a high-quality turbulent flow as needed to predict the acoustic sources. However, this method requires a time-averaged mean velocity field and the full Reynolds stress tensor as input values. In RANS simulations, the distribution of the Reynolds stresses is not modeled by the commonly used turbulence models. Therefore, we developed a simple, robust approach to model the distribution of the Reynolds stresses based on scaling of the turbulent kinetic energy obtained from a RANS simulation.

In the acoustic propagation solver, we also solve the acoustic perturbation equations (APE) in version 4 with the high-order discontinuous Galerkin scheme. Due to the beneficial dissipation and dispersion properties, high-order schemes, such as discontinuous Galerkin, are well suited. In vortex noise-dominated problems, we use the perturbed Lamb vector as an acoustic source. In the region where acoustic production is assumed, e.g., at the trailing edge of an airfoil, the Lamb vector is extracted from the zonal LES, solving the compressible NSE. The imbalance between the acoustic and flow solver is taken into account by a static load balancing.

The new proposed model of the Reynolds stresses was validated using the example of a turbulent flow over a flat plate in the first step. As a reference, a wall resolved tripped boundary layer over a flat plate was computed. In the following, the rear part of the boundary layer was simulated with a zonal LES using the RRALF method with the target data from the tripped LES. We were able to show that this method is capable of almost exactly reproducing the tripped LES when using the same numerical scheme and the same mesh resolution. A second zonal LES based on target data derived from a RANS simulation applying the new proposed model is validated against this reference zonal LES. Here, the differences can be accounted to the errors introduced by the model. In the wall near region, the trend is to underestimate the Reynolds stresses, whereas in the outer region of the boundary layer, the Reynolds stresses are slightly overestimated. Overall, the results of the zonal LES using the new model approach are in good agreement with the references. This example demonstrates that the new proposed model can generate a realistic turbulent field based on two-dimensional RANS simulation data.

Next, the direct-hybrid framework was validated at the example of a NACA 0012 airfoil. The investigated configuration is well known from the literature to exhibit tonal self-noise. We showed that the compressible flow solver can resolve the interaction of the acoustics and the hydrodynamics. The hydrodynamic instabilities excited by the acoustics form the vortices generating self-noise at the trailing edge. These vortices are extracted as an acoustic source, and we demonstrated that the direct-hybrid acoustic solver is capable reproducing the tonal component in the acoustic spectrum. This is the advantage of computing a compressible simulation compared to an incompressible one. In the case of an incompressible flow simulation, hydrodynamics and acoustics are strictly separated, and tonal self-noise can not be captured.

The turbulent inflow based on RANS simulation data and the acoustic solver, both building blocks of the zonal direct-hybrid simulation, were applied to the example of the simulation of a NACA 64418 airfoil’s trailing edge. The presented framework is capable of reproducing the mean velocity field nearly identical compared to the reference. In the inflow region, the distribution of the Reynolds stresses differs from the reference due to the recycling approach and the pressure gradient. However, after the turbulent inflow region and in the region at the trailing edge and in the wake, the distribution and amplitude match the reference well. This proves the proposed model approach for the Reynolds stress distribution combined with the turbulent inflow method to be well suited, even at more complex cases with pressure gradient. The acoustic emissions were compared between a direct acoustic simulation of the LES of the full airfoil, the DNC results of the zonal LES, and the results of the acoustic solver. At high frequencies, the results of the three simulations show similar behavior. In the medium range of the considered frequencies, the acoustic solver overpredicts the acoustic level. At low frequencies, the acoustic results of the zonal LES show the influence of the turbulent inflow. Here, artificial noise is produced, and a higher noise level is predicted compared to the full LES. Based on the extracted acoustic source terms from the zonal LES, the direct-hybrid simulation omits the acoustic sources in the inflow region. This can strongly reduce the amount of artificial noise in the acoustic simulation, and the acoustic behavior of the full LES can be approximated.

The acoustic results of the zonal direct-hybrid simulation approach show great potential. The zonal computation of the acoustic sources opens the possibility to apply high-fidelity methods for the calculation of acoustic sources efficiently. Combined with an acoustic propagation solver, the propagation of the acoustics beyond the restricted LES domain to predict far-field noise is feasible. Furthermore, solving the compressible NSE in the flow solver allows depicting acoustic feedback in the acoustic solution. Applying the direct hybrid approach further enables the efficient execution of direct-hybrid aeroacoustic simulations on high-performance computers.

Conflicts of interest

Authors declared no conflict of interests.

Data availability statement

Data are available on request from the authors.

Acknowledgments

The authors gratefully acknowledge the support and computing resources granted by the High Performance Computing Center Stuttgart (HLRS) on the national supercomputer HPE Apollo (Hawk) under the grant hpcdg. Also we want to acknowledge the support of the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany’s Excellence Strategy – EXC 2075 – 390740016. We further thank Nico Krais for providing the baseline setup of the NACA 64418 profile.

References

All Figures

	Figure 1 Reynolds stresses normalized by the local turbulent kinetic energy k of a turbulent boundary layer.
In the text

	Figure 2 Sketch of the mesh used for the tripped and zonal simulations of the turbulent flow over a flat plate. Zonal region reaching from Reθ = 1800–2800.
In the text

	Figure 3 Comparison of the mean velocity profile between the full LES (tripped), the RANS simulation, and ZLES II (zonal LES based on the RANS simulation data) at the inflow of the zonal LES. The mean velocity profile of the RANS simulation and the ZLES are identical.
In the text

	Figure 4 Comparison of the Reynolds stresses between the full LES (tripped), the RANS simulation (modeled), and zonal LES II (based on the RANS data) at the inflow of the zonal LES. The used values for the scaling tensor Sij are S11 = 0.9, S22 = 0.45, S33 = 0.65 and S12 = 0.3. The area where the ALF is applied is indicated in blue.
In the text

	Figure 5 Comparison of the mean velocity profiles of the turbulent boundary layer between ZLES I (target data based on LES data), ZLES II (target data based on a RANS simulation), the full LES (tripped), and reference data from the literature.
In the text

	Figure 6 Comparison of the Reynolds stresses of the turbulent boundary layer between ZLES I (target data based on LES data), ZLES II (target data based on a RANS simulation), the full LES (tripped), and reference data from the literature.
In the text

	Figure 7 Time-dependent dimensionless pressure recorded at an arbitrary time one chord length above the trailing edge. Comparing the DNC results (blue) with the direct-hybrid simulation results (red).
In the text

	Figure 8 Comparison of the dimensionless PSD of the pressure between the DNC results (blue) and the direct-hybrid simulation results (red).
In the text

	Figure 9 Trailing edge of the NACA 64418 airfoil. Vortex structure visualized by the Q-criterion (Q = 200) colored by the streamwise velocity in the range [−0.3, 1.3]. The displayed planes indicate the inflow and recycling planes. The forcing region is highlighted in blue.
In the text

	Figure 10 Comparison of the modelled Reynolds stresses and the results of a wall resolved LES. Plotted in wall-normal direction on the suction side of a NACA 64418 airfoil. The used values for the scaling tensor Sij are S11 = 0.8, S22 = 0.55, S33 = 0.65 and S12 = 0.3. The area where the ALF is applied is indicated in blue.
In the text

	Figure 11 Pressure distribution on the surface of the tripped NACA 64418 airfoil. Comparison of the LES and the RANS simulation of the complete airfoil and the zonal LES based on the RANS simulation data.
In the text

	Figure 12 Mean velocity profiles in comparison to the full LES at the displayed positions. Profiles staggered in streamwise direction. Top left: Suction side. Top right: Mean streamwise velocity, including the positions of the velocity fluctuation profiles. Bottom left: Pressure side. Bottom right: Wake region.
In the text

	Figure 13 Reynolds stress distributions in comparison to the full LES at the displayed positions. Profiles staggered in streamwise direction. Top left: Suction side. Top right: Mean streamwise velocity, including the positions of the velocity fluctuation profiles. Bottom left: Pressure side. Bottom right: Wake region.
In the text

	Figure 14 Dimensionless PSD of the pressure above the trailing edge. Comparison of the direct-hybrid simulation, the zonal DNC, and the DNC of the full LES.
In the text

	Figure 15 Dimensionless PSD of the pressure below the trailing edge. Comparison of the direct-hybrid simulation, the zonal DNC, and the DNC of the full LES.
In the text