NLAS provides a more sophisticated subgrid treatment that allows the extraction of acoustic sources from the temporal variation within the modeled subgrid structures. The quasi-steady near-wall RANS solution is obtained a priori so that the grid requirements can be relaxed and reduced in the near-wall region during the NLAS transient calculation, compared to the LES solvers. At the same time, the dissipation effects of a subgrid eddy viscosity model are avoided; Thus, the NLAS solver proves less dissipation than the classic LES or hybrid RANS/LES simulation on course meshes (Figure 2). One of the most important advantages of the NLAS is able to account for both the turbulence-related broadband noise and the discrete tones produced from coherent structures or resonance [24].

The NLAS control equation is derived from the N-S equations by dividing each term in the equation into the mean and the fluctuation terms, ϕ =   ϕ ̄ +   ϕ ′ substituting into the N-S equations, and rearranging the mean and the fluctuation terms to obtain the set of non-linear perturbation Eq. 1 [23, 24]: ∂ Q ′ ∂ t + ∂ F i ′ ∂ x i − ∂ F i v ′ ∂ x i = − ∂ Q ̄ ∂ t − ∂ F i ̄ ∂ x i + ∂ F i v ̄ ∂ x i (1) where.

Q′ is the transient fluctuation term;

Q ̄ is the transient mean term;

F i ′ is the linear inviscid fluctuation term;

F i ̄ is the inviscid averaging term;

( F i v ) ′ is the viscous fluctuation term;

F i v ̄ is the viscous averaging term.

The specific expressions are: Q ′ =   ρ ′ ρ ̄ u j ′ + ρ ′ u j ̄ + ρ ′ u j ′ e ′ , Q ̄ =   ρ ̄ ρ ̄ u j ̄ e ̄ F i ′ =   ρ ̄ u i ′ + ρ ′ u i ̄ ρ ′ u i ̄ u j ̄ + ρ u i ̄ u j ′ + ρ ̄ u i ′ u j ̄ + p ′ δ i j u i ′ e ̄ + p ̄ + u i ̄ e ′ + p ′ +   ρ ′ u i ′ ρ ̄ u i ′ u j ′ + ρ ′ u i ′ u j ̄ + ρ ′ u i ̄ u j ′ + ρ ′ u i ′ u j ′ u i ′ e ′ + p ′ F i v ′ =   0 τ ′ i j − θ ′ i + u ′ k τ ̄ k i + u ̄ k τ ′ k i F i ̄ =   ρ ̄ u i ̄ ρ ̄ u i ̄ u j ̄ + p ̄ δ i j u i ̄ e ̄ + p ̄ , F i v ̄ =   0 τ i j ̄ − θ i ̄ + u ̄ k τ k i ̄

Neglecting the density fluctuation and taking a time average of the Eq. 2: L H S ̄ = R H S ̄ = ∂ R i ∂ x i (2) R _i is a correlation of the standard Reynolds stress tensor and turbulent heat flux: R i =   0 ρ ̄ u ′ i u ′ j ̄ c p ρ T ′ u ′ i ̄ + ρ ̄ u ′ i u ′ k ̄ u k ̄ + 1 2 ρ ̄ u ′ k u ′ k u ′ i ̄ + u ′ k τ k i ̄

The key step in NLAS is to obtain these unknown terms in advance from the classical RANS method. The small-scale quantities which cannot be solved are obtained by turbulent reconstruction generated by the subgrid source terms. The turbulent reconstruction method is proposed by Batten [25]. After both the mean levels and subgrid sources are established, the time-dependent calculations can then be carried out to determine the transmitted perturbations around the mean flows by using the above non-linear disturbance equations. The sound pressure in far-field observers can be calculated by the FW-H equations given by Farassat et al. [31, 32]. The flow field information is obtained by NLAS.

In order to reduce the amount of mesh during the numerical simulation, the numerical models are usually modified or simplified from the experiment models. For example, the span of the experiment airfoil L is much bigger than the chord c. In order to use the lower mesh number, the span of the numerical simulation can be reduced from L to L _s , which is smaller than the chord c. Then, the aeroacoustic calculation can be speeded up extensively; However, the aeroacoustic sound pressure level (SPL) obtained from the numerical results, and the experimental results cannot be compared directly. In such cases, some corrections have to be introduced to the numerical sound pressure level (SPL). In the paper, we use the correction method first proposed by Kato et al [33, 34].

When L C ≤ L S S P L = S P L S + 10 ⁡ log L / L S (3) when L S < L C ≤ L S P L = S P L S + 20 ⁡ log L C / L S + 10 ⁡ log L / L C (4) when L C > L S P L = S P L S + 20 ⁡ log L / L S (5)

The validation benchmark is given by the experiment data of Jacob [35]. In the experiment (Figure 3), the NACA0012 airfoil (with chord length c = 0.1 m) is located at one-chord length after the cylinder (with diameter d = 0.01 m), both extending by L = 0.3 m in the spanwise direction. The acoustic far-field receiver is at 1.85 m from the airfoil center. The incoming flow velocity is 72 m/s and the Reynolds number of the cylinder is 48,000. The experiment was conducted in an anechoic room (10 m × 8 m × 8 m). The air was supplied by a high-speed subsonic anechoic wind tunnel at Mach numbers ranging up to 0.34.

In order to reduce the computational cost, the span of the numerical model is set as 0.05 m which is smaller than the experiment model (L = 0.3 m). The simulation domain was X (−0.2 m, 0.3 m), Y (−0.2 m, 0.2 m), and Z (−0.05 m, 0 m). After the grid convergence check, the multiblock structure meshes with 3 million (mesh cell number) computational grids were used for the RANS/NLAS simulation. The first interior point was located at y+ < 1 from the airfoil surface, yielding a sufficient resolution of the viscous sublayer. In the RANS simulation, the cubic k-e turbulence model was used. NLAS provided a more sophisticated subgrid treatment that allowed the extraction of acoustic sources from the temporal variation within the modeled subgrid structures. In this paper, the subgrid for NLAS was x (−0.15 m, 0.15 m), y (−0.1 m, 0.1 m), and z (−0.05 m, 0 m) with the resolution ratio 0.002.

Four offset distances (D = 0d, 0.5d, 1d, and 3d) are used to reveal the effect of offset distances. The simulation domain and offset distances can be found in Figure 4.

RANS calculations were carried out to obtain a statistical average of the local Reynolds stress tensor and heat flux. Then they were interpolated into the NLAS noise calculation grid. The turbulence flow is reconstructed based on this statistical average. The velocity inlet was applied in the left boundary. The outflow boundary was used in the right boundary. The far-field and periodic boundary are used in the up/down and the spanwise direction. The viscous non-slip adiabatic wall condition was applied for cylinders and airfoils.

The flow velocity (on the line x/c = 0.25 and 0 < y/c < 0.6, as show in Figure 5A.), obtained from the RANS, was compared with the experimental results in Jacob’s literature ([35]), as shown in Figure 5B. It can be seen that the velocity variation trend obtained from the RANS calculation in this paper is the same as the URANS method used in the literature ([35]). Both methods differ significantly from the experimental variation trend within 0.06 < y/c < 0.26; however, for the region y/c > 0.26, the results of the RANS calculations in this paper are closer to the experimental data than the URANS method.

The instantaneous Cl for the cylinder is given in Figure 6. The Cl presents a sinusoidal oscillation mode which means the convergence.

Figure 7 shows the vorticity contours obtained from the URANS calculations, the large-scale vortices at the trailing edge of the cylinder are continuously shedding and spreading downstream to the leading edge of the wing. The wake of the cylinder interacts with the leading edge of airfoils and generates unsteady pressure fluctuation across the airfoil surface.

The RANS calculation results were interpolated to the NLAS grid with turbulent reconstruction for acoustic calculations. The numerical time step is Δt = 2 × 10⁻⁵ s, and the simulation is performed 4,000 steps in total, the simulation is performed using AMD EPYC 7452 with 64 cores for about 20 h. In this paper, the numerical model span is smaller than the experimental model span, and the correction of Eq. 5 is used to correct the acoustic results accordingly.

Figure 8 shows the comparison between the NLAS used in this paper and several calculation methods (LES [35], URANS [13], DES [9]) in the literature. It is found that the sound pressure levels obtained by all four calculation methods are higher than the experimental values in the low-frequency range of 300 Hz–900 Hz. The peak sound pressure obtained by the NLAS method is almost the same as the experimental data and is more accurate than the other methods. Above 1700 Hz, the results obtained by EASM DES and NLAS method are closer to the experimental values. It shows that the NLAS method is able to capture the sound pressure accurately in most frequencies and with an acceptable computational cost.

Four different sets of computational mesh are used to test the grid independence. As shown in Figure 9, the peak sound pressure level is 91 dB at 1,365 Hz in the experiment [35]. The difference for different computational meshes is quite small. The peak sound pressure level and its frequency are compared in Table 1. The maximum difference of peak sound pressure level is less than 4 dB, and the peak frequency is less than 6 dB, which indicates the calculated results are less correlated with the computational grid. In this paper, the number of meshes is chosen to be 2.41 million to give an error of less than 2 dB for both the peak sound pressure level and its frequency.

In order to analyze the influence of wakes at the trailing edge of the cylinder, different offset distances (0d, 0.5d, 1d, and 3d) are used in this section. The sound pressure is collected in 1.85 m from the center of the wing at zero offset distance. This location is not changed in later research.

Figure 10 shows the variation of vortex intensity at different offset distances calculated by URANS, where the vertical coordinate indicates the vortex intensity ω and the horizontal coordinate indicates distances behind the rod center. The figure shows that the vortex intensities in the wake of the cylinder show a non-linear decreasing trend, and the offset distance has less effect on them. The vortex in the wake of the cylinder dissipates due to the viscous effect of the fluid. The dissipation rate of the vortex is found to be: y = 2615 + 432 e − 4.8 x (6) where y is the biggest vortex intensity in the wake of the rod and x is the distance behind the rod center.

Figure 11 shows the pressure distribution of the airfoil along the chord in the middle span. It is found that only when the offset distance D = 0d the pressure on the upper and lower surfaces of the airfoil shows symmetric changes. With the increase of offset distance in the region near the leading edge (0.014 < x/c < 0.3), the pressure gradually decreases and with a asymmetric distribution. This may be caused by the fact that with the increase of the offset distance, the higher the flow velocity interacts with the leading edge of the airfoil. When the offset distance increases, the airfoil gradually get out of the vortex street influence area, and the flow field around the body is dominated by the velocity field without cylindrical disturbance.

The pressure distribution at the leading edge of airfoils along the spanwise is presented in Figure 12. The pressure in the center of the leading edge is affected by the vortex shedding dramatically. The pressure difference in the center of the leading edge gets its maximum for about 1,300 Pa by comparing D = 0d and D = 3d. This difference decreases with the increase of offset distance D.

The size and mesh resolution for the integral surface of sound sources are the same as the benchmark case. After 4,000 steps of NLAS calculation, Figures 13, 14 show the time domain of sound pressure in the time domain and frequency domain, respectively. The variation of peak sound pressure and its frequency are plotted in Figure 15. It can be seen that the sound pressure level decreases gradually with the increase of offset distance, from 89 dB to 80.6 dB. The peak sound frequency decreased significantly when the offset distance changed from 0d to 0.5d. Its frequency did not change significantly when the offset distance changed from 0.5d to 1d and increased slightly from 0.5d to 3d. With the increase of offset distance, the sound pressure level decreases almost linearly from 0d to 1d; From 1d to 3d, the decrease rate slows down significantly.

The Q-criteria contour surfaces (Figure 16) obtained from NLAS show that there are both large-scale vortices and small-scale vortices generated in the wake of the cylinder. With the increase of the offset distance, the interaction between the wake of the cylinder and the leading edge of the airfoil becomes smaller. So the sound pressure decreases, as predicted in Figure 15.