© The Author(s), Published by EDP Sciences, 2024

1 Introduction

Many numerical simulations of flows around land or air vehicles or isolated components are carried out for acoustic purposes. These numerical simulations aim to meet the ever-growing need for accurate predictions of the noise generated by aerodynamic flows. Articles on the state of the art of calculation methods are regularly published (see for example [1, 2]). Most of the time, these aeroacoustic simulations are performed in two stages for efficiency: an aerodynamic simulation of the flow concentrating the computational effort (grid density, order of the resolution scheme, etc.) in a volume as small as possible around the object in order to limit the numerical cost, then an acoustic analysis of the field of disturbances. The most widely used tools for this acoustic analysis are still currently based on the integral surface formulations of Kirchhoff [3] or of Ffowcs Williams and Hawkings (FW-H) [4]. Starting from the flow disturbances on a surface including the noise sources, assuming uniform flow and acoustic propagation in free field outside this surface, these integral methods make it very easy to deduce the noise radiated outside the integration surface, up to the far field. Taking into account the reflection or shadow effects of solid bodies outside the integration surface – such as a wing or a fuselage in the noise radiated by a propeller or a jet, for example – can be a significant improvement. A solution may be to include these reflecting bodies in the aerodynamic calculation and to extend the integration surface beyond these bodies, but the cost of the aeroacoustic simulation is considerably increased. Indeed, it is then necessary to adapt the grid density and/or the order of the calculation for a correct propagation of acoustic signals of low amplitude (compared to the aerodynamic perturbations) in a much larger domain. Coupling an aerodynamic calculation with a propagation calculation by solving Euler’s equations (CFD/CAA coupling, [5–8]) is more economical for the aerodynamic part but remains quite expensive overall. A coupling with boundary element methods (BEM) in the frequency domain is also possible [9–11], but this coupling can also turn out to be quite expensive if the acoustic signals are broadband and extend to high frequencies, such as with blade vortex interaction (BVI) noise of a helicopter rotor, high-speed impulsive (HSI) noise of a transonic rotor or propeller, or broadband noise of a jet. Ray-tracing methods can also be used to predict noise scattering from solid bodies [12–15], but they require the noise sources behind the incident sound field to be modeled. In an application to a numerical aerodynamic simulation, the validity of the result is therefore closely linked to the source model.

The idea is to be able to evaluate the main reflection and shadow effects of a body, by a fast “first-order” method, without weighing the aerodynamic computation down. This approximate calculation can also show, depending on the results, whether a more precise but more expensive calculation is necessary before performing it. For example, a more precise calculation may be unnecessary if the reflected levels are deemed negligible. Conversely, it may be necessary if the levels found show a strong gradient in a critical region of interest. The simplified method can also speed up the optimization of the position of a reflective surface with the aim of reducing noise levels in a given direction, and save an elaborate calculation for the most promising configuration, to ensure an accurate level prediction. However, the method should not be limited to academic geometries or noise sources. It must work for bodies of any geometry (fuselage, wing, etc.) and noise sources whose spatial and spectral extent can be relatively large (jet noise, for example).

The solution proposed in this study is based on the Kirchhoff method in the time domain. It is therefore directly applicable to the flow disturbances provided by unsteady numerical simulations, without prior spectral analysis. It consists of a first calculation of free-field radiation to determine the incident acoustic pressure on a fictitious surface representing the skin of the reflecting body, then in a second radiation calculation to deduce the effects of this reflecting body on the complete acoustic field. This approach relies on two simplifying assumptions: a uniform flow around the body and short wavelengths compared to the local radius of curvature of the body surface. It can be noted that the Kirchhoff method has already been used for simplified reflection calculations [16–20], generally in the frequency domain, but the approach presented in this article is original in the simplicity and efficiency of the calculation algorithm that does not require any information on the origin of the incident acoustic field.

The article is organized as follows. In Section 2, the expression of the Kirchhoff surface integral is first reviewed and the method for calculating reflection and shadow effects is described. In Section 3, the consequences of the assumption of short wavelength compared to the local radius of curvature of the reflecting body are quantified by comparison with analytical solutions. An application example is then given for a wing in a Mach 0.5 flow, in Section 4. The calculation cost and a first assessment of the method are provided in the conclusion, together with leads for possible improvements.

2 Formulation

Let us define a reference frame Y in which the unperturbed fluid is in translation at a uniform velocity $U_0{\vec{Y}}_1$ and S, the solid surface of a body fixed in this reference frame. Let x_i and y_i respectively be the space coordinates of an observation point $\vec{x}$ and a point $\vec{y}$ of surface S (Fig. 1). Let p_inc be the pressure radiated on S by noise sources (that do not need to be determined for the reflection calculation).

Figure 1 Schematic of acoustic waves towards the solid body and reflection towards the observation point.

The radiation of the pressure $p\left(\vec{x},t\right)$ reflected by S at the observation point $\vec{x}$ can then be calculated by a Kirchhoff integral starting from the reflected pressure p_ref on S and its gradient $\vec{\nabla }{p}_{\mathrm{ref}}$, that can expressed as (see Appendix A):

$$\begin{array}{c}p\left(\vec{x},t\right)=\underset{\tau }{\overset{}{\int }}\underset{S}{\overset{}{\int }}{F}_{\mathrm{ref}}\left(\vec{x},\vec{y},\tau \right)\mathrm{GdSd\tau }\\ {\mathrm{with} F}_{\mathrm{ref}}\left(\vec{x},\vec{y},\tau \right)=p_{\mathrm{ref}}\frac{{\beta }^2}{d^2}\left(\vec{x}-\vec{y}\right)·\vec{n}-\frac{\partial p_{\mathrm{ref}}}{\partial n}+M_0^2\frac{\partial p_{\mathrm{ref}}}{\partial y_1}{n}_1+\frac{1}{a_0}\left(\frac{\left(\vec{x}-\vec{y}\right)·\vec{n}}{d}+M_0{n}_1\right)\frac{\partial p_{\mathrm{ref}}}{\partial \tau }\end{array}$$(1)

where G = δ(g)/4πd the free-space Green’s function in R₃ for the convective wave equation, δ(f) the Dirac distribution, g = τ – t + σ/a₀, σ = (d − M₀ (x₁ − y₁))/β², $d=\sqrt{C_i{\left(x_i-y_i\right)}^2}$, C_i = (1 − β²) δ_1i + β², δ_ij the Kronecker delta symbol (1 if i = j, zero otherwise), ${\beta }^2=1-M_0^2$, M₀ = U₀/a₀, a₀ the speed of sound in the unperturbed medium, $\vec{n}$ the normal unit vector on the surface S, pointing outward the solid body, $\frac{\partial }{\partial n}=\frac{\partial }{\partial y_i}{n}_i$, τ an emission time and t a reception time.

The surface of the reflecting body is assumed to be rigid and parallel to the mean flow. The latter assumption should have limited consequences for aeronautical applications, where the mean flow approaches most aircraft surfaces with a low angle of incidence.

If the body’s local radius of curvature is large enough compared to the wavelength, S can be locally considered as a plane surface from the acoustic point of view. The reflected pressure p_ref on S and its gradient $\vec{\nabla }{p}_{\mathrm{ref}}$ can then be simply deduced from the incident pressure p_inc and its gradient $\vec{\nabla }{p}_{\mathrm{inc}}$. The locally plane reflection hypothesis leads to the following relations:

$$\begin{array}{c}{p}_{\mathrm{ref}}=p_{\mathrm{inc}}\\ \vec{\nabla }{p}_{\mathrm{ref}}=\vec{\nabla }{p}_{\mathrm{inc}}-2\frac{\partial p_{\mathrm{inc}}}{\partial n}\vec{n}\end{array}$$

which lead to:

$$\begin{array}{c}\frac{\partial p_{\mathrm{ref}}}{\partial \tau }=\frac{\partial p_{\mathrm{inc}}}{\partial \tau }\\ \frac{\partial p_{\mathrm{ref}}}{\partial n}=\frac{\partial p_{\mathrm{inc}}}{\partial n}-2\frac{\partial p_{\mathrm{inc}}}{\partial n}\\ \frac{\partial p_{\mathrm{ref}}}{\partial y_1}{n}_1=\frac{\partial p_{\mathrm{inc}}}{\partial y_1}{n}_1-2\frac{\partial p_{\mathrm{inc}}}{\partial n}{n}_1^2.\end{array}$$

Expressing p_ref, $\frac{\partial p_{\mathrm{ref}}}{\partial n}$ and $\frac{\partial p_{\mathrm{ref}}}{\partial y_1}{n}_1$ in formula (1) by their expressions above and rearranging leads to:

$$F_{\mathrm{ref}}\left(\vec{x},\vec{y},\tau \right)=p_{\mathrm{inc}}\frac{{\beta }^2}{d^2}\left(\vec{x}-\vec{y}\right)·\vec{n}-\frac{\partial p_{\mathrm{inc}}}{\partial n}+M_0^2\frac{\partial p_{\mathrm{inc}}}{\partial y_1}{n}_1+2\left(1-M_0^2{n}_1^2\right)\frac{\partial p_{\mathrm{inc}}}{\partial n}+\frac{1}{a_0}\left(\frac{\left(\vec{x}-\vec{y}\right)·\vec{n}}{d}+M_0{n}_1\right)\frac{\partial p_{\mathrm{inc}}}{\partial \tau }$$(2)

The source term $F_{\mathrm{ref}}\left(\vec{x},\vec{y},\tau \right)$ of the reflected pressure can thus be expressed entirely as a function of the incident pressure and its gradient, which are input data assumed to be known.

In practice, in order to best model the reflection and shadow effects of the solid body, four cases must be distinguished for each surface element ΔS of its surface S. These four cases are illustrated in Table 1:

Case (a): the incident wave is incoming with respect to the surface element, and the observation point is on the same side as the local normal;

Case (b): same configuration for the incident wave, but the observation point is on the opposite side to the local normal;

Case (c) and case (d): the incident wave is outgoing with respect to the surface element.

Formula (2) only applies to an observer on the reflection side (case (a)). For an observer on the other side (shadow side, case (b)), the shadow effect of the body can be modeled in the calculation by a term creating a radiation opposite to the incident wave, so that the sum with the direct (free-field) pressure is zero. By taking into account the orientation of the local normal $\vec{n}$, the corresponding source term is simply formula (1) applied to p^inc:

$$F_{\mathrm{shadow}}\left(\vec{x},\vec{y},\tau \right)=p_{\mathrm{inc}}\frac{{\beta }^2}{d^2}\left(\vec{x}-\vec{y}\right)·\vec{n}-\frac{\partial p_{\mathrm{inc}}}{\partial n}+M_0^2\frac{\partial p_{\mathrm{inc}}}{\partial y_1}{n}_1+\frac{1}{a_0}\left(\frac{\left(\vec{x}-\vec{y}\right)·\vec{n}}{d}+M_0{n}_1\right)\frac{\partial p_{\mathrm{inc}}}{\partial \tau }.$$(3)

Finally, there is no reflection or shadow effect to calculate in the case of the incident wave coming out of the body surface (outgoing wave, cases (c) and (d)). The source term corresponding to these cases is therefore zero.

For each point on the surface S and for each observation point, the direction of propagation of the incident wave and the relative position of the observation point with respect to the reflecting element can be determined by the sign of (∂p/∂n)(∂p/∂τ) and of $\left(\vec{x}-\vec{y}\right)·\vec{n}$ as shown in Table 1.

All cases can therefore be easily managed by two coefficients κ₁ and κ₂ in the following expression of the source term:

$$F\left(\vec{x},\vec{y},\tau \right)={\kappa }_1\left(\begin{array}{c}{p}_{\mathrm{inc}}\frac{{\beta }^2}{d^2}\left(\vec{x}-\vec{y}\right)·\vec{n}\\ -\frac{\partial p_{\mathrm{inc}}}{\partial n}+M_0^2\frac{\partial p_{\mathrm{inc}}}{\partial y_1}{n}_1\\ +\frac{1}{a_0}\left(\frac{\left(\vec{x}-\vec{y}\right)·\vec{n}}{d}+M_0{n}_1\right)\frac{\partial p_{\mathrm{inc}}}{\partial \tau }\\ +2{\kappa }_2\left(1-M_0^2{n}_1^2\right)\frac{\partial p_{\mathrm{inc}}}{\partial n}\end{array}\right).$$(4)

In the general case, the direction of the incident wave is not constant, and the coefficients κ₁ and κ₂ must obviously be determined for each time step of the reflection calculation. The discontinuities created by this coefficient switching during a time-domain calculation have no major impact on the results, mainly because the influence of these discontinuities is linked to the size of the elementary surface ΔS of the reflecting surface S, and the contribution decreases on finer meshes. The space-time shifting of these discontinuities on the surface S, during the simulation, also tends to spread their effects over time at the observation point.

In summary, a reflecting body can be taken into account in an approximate way in a radiation calculation in free field starting from the knowledge of the incident pressure, its time derivative and its gradient on the body surface. The incident pressure on the body surface is provided by the radiation calculation in free field and its time derivative can be easily deduced by finite differences. So, only the pressure gradient remains to be determined. In the general case, this can also be calculated by finite differences of the pressure on a double-layered mesh on either side of the body skin. This approach is used in the calculations that follow. It has the disadvantage of requiring to define the distance between the two layers so that the calculation of the pressure gradient is precise enough for the frequency band of interest (thickness small enough compared to the wavelengths while remaining compatible with the precision of the incident radiation calculation). A more accurate and safer way to proceed is possible. It consists in calculating the pressure gradient directly by the same integral formulation as that used for the pressure. Expressions of this gradient are provided in Appendix B, for fixed integration surfaces, for both Kirchhoff and FW-H integral formulations.

Remarks:

• For a correct summation of direct and reflected pressures (Fig. 2), it is of course necessary for both pressures to be calculated with the same time reference point, so that the causal link between both pressures at the observer point is maintained. The calculation of the reflection of each elementary surface ΔS of S must therefore be carried out from the reception times of the direct radiation on this elementary surface, so that the final reception times at the observer point include all propagation times.

  Figure 2 Illustration of the effect of propagation times on the duration of the valid total signal (case of non-periodic signals).

• In the case of non-periodic acoustic fields (for which the signals cannot be folded back on a single period), the different propagation times must be taken into account to determine the part of the total signal that is valid, i.e. the part that has received both the direct signal and the entire contribution of the reflecting surface, as shown in Figure 2.

• The reflection calculation can include several reflective surfaces, but, for each surface, it only takes into account the direct radiation, not the radiation of other reflective surfaces. The method is therefore not intended for acoustic propagation in ducts or for the prediction of stationary waves that can form a standing wave pattern between two reflective surfaces.

3 Comparison with analytical solution

In this section, reflection calculation using the Kirchhoff surface integral is applied to the prediction of the scattering of a point source by a rigid sphere, and the results are compared with the analytical solution. The configuration is as follows: a sphere of radius R = 0.17 m in a fluid at rest (M₀ = 0), a monopole point source at frequency f placed at x_S = 0.34 m from the center of the sphere in its horizontal median plane z = 0 (Fig. 3), 180 observation points evenly distributed on a circle of radius r = 0.6 m around the sphere, in the horizontal plane. The calculations are performed for a speed of sound a₀ = 340 m/s and for four frequencies: f = 1 kHz (λ/R = 2), f = 2.5 kHz (λ/R = 0.8), f = 5 kHz (λ/R = 0.4) and f = 7.5 kHz (λ/R = 0.27).

Figure 3 Calculation of the scattering of a point source by a rigid sphere.

The analytical solution of the scattered pressure in complex notations, at the observation points (r, θ) is [21]:

$$\begin{array}{c}{P}_T\left(r,\theta \right)=\frac{\mathrm{ik}}{4\pi }\sum _{m=0}^{\infty }\left(2m+1\right)P_m\left(\cos \theta \right)\mathrm{\times }\left[j_m\left(\mathrm{kr}\right)h_m\left(k{x}_S\right)-\frac{j_m^{\mathrm{\prime }}\left(\mathrm{kR}\right)}{h_m^{\mathrm{\prime }}\left(\mathrm{kR}\right)}{h}_m\left(k{x}_S\right)h_m\left(\mathrm{kr}\right)\right]\hspace{1em}\mathrm{r}\le {\mathrm{x}}_{\mathrm{S}}\\ P_T\left(r,\theta \right)=\frac{\mathrm{ik}}{4\pi }\sum _{m=0}^{\infty }\left(2m+1\right)P_m\left(\cos \theta \right)\mathrm{\times }\left[j_m\left(k{x}_S\right)h_m\left(\mathrm{kr}\right)-\frac{j_m^{\mathrm{\prime }}\left(\mathrm{kR}\right)}{h_m^{\mathrm{\prime }}\left(\mathrm{kR}\right)}{h}_m\left(k{x}_S\right)h_m\left(\mathrm{kr}\right)\right]\hspace{1em}\mathrm{r}>{\mathrm{x}}_{\mathrm{S}}\end{array}$$

with i² = −1, k = 2π/λ and where P_m is the Legendre polynomial of degree m, j_m the first kind spherical Bessel function of order m, $j_m^{\mathrm{\prime }}$ its derivative, h_m = j_m + iy_m the Hankel function of the first kind and order m, $h_m^{\mathrm{\prime }}$ its derivative and y_m the second kind spherical Bessel function of order m. Below, these series are calculated for a number of terms that ensure a error of less than 10⁻¹⁰.

For the Kirchhof method, the surface of the sphere is represented by a uniform unstructured mesh composed of triangles. This mesh is duplicated on either side of the sphere surface (radial distance ± 0.00034 m) for the calculation of the pressure gradient. The complex pressure radiated by the monopole on a node $\vec{x}$ of the double-layer mesh and at an instant t, is [22, 23]:

${\widehat{p}}^{\mathrm{\prime }}\left(\vec{x},t\right)={\rho }_0\left(\frac{\partial \varphi }{\partial t}+U_0\frac{\partial \varphi }{\partial x_1}\right)$ where $\varphi \left(\vec{x},\vec{y},t\right)=\frac{A}{4\pi d}{\exp }^{i\left(2\pi f\right)(t-\frac{\sigma }{a_0})}$ is the complex potential of the monopole, with A its strength and $\vec{y}=(0.34, 0, 0)$. The input to Kirchhoff’s calculation is the real part of this pressure:

$$p^{\mathrm{\prime }}\left(\vec{x},t\right)=\frac{A{\rho }_0}{4\pi d}\left[2\pi f\left(1-\frac{U_0}{a_0}\frac{\partial \sigma }{\partial x_1}\right)\sin \left(2\pi f\left(t-\frac{\sigma }{a_0}\right)\right)+\frac{U_0}{d}\frac{\partial d}{\partial x_1}\cos (2\pi f\left(t-\sigma /a_0)\right)\right].$$

The computation time step is equal to 1/100 of the period T of the monopole.

As a first step, calculations are carried out for the highest frequency (f = 7.5 kHz), for three meshes of different densities, to check the stability of the results. The number of cells, the average distance L between two nodes of a cell, and the λ/L ratio are given in Table 2, for the three meshes and for the four frequencies. These values are given for guidance only. Since the angle of incidence of the acoustic waves on the sphere varies from 0° to 90°, choosing a uniform mesh based on the λ/L criterion alone would be unnecessarily penalizing. The criterion considered is simply the minimum mesh density at which the results no longer change significantly.

The results for the three meshes are compared in Figure 4. The quantity plotted is the attenuation factor Gamma Total (GT) defined as γ_Τ(f) = <p_T(f)>/<p_I(f)>, with < . > indicating the ensemble average operator, and p_I and p_T, the incident and the total pressure fluctuations, respectively.

Figure 4 Influence of the mesh density on the Kirchhoff results (f = 7.5 kHz).

There are no significant differences between the directivities obtained with the medium and fine meshes. On the other hand, the coarse mesh leads to significantly different levels for observation angles between 45° and 315°. However, the differences remain moderate, indicating a certain robustness to the method. Given these results for the highest frequency, comparisons with analytical solutions are made with the medium mesh.

For all four frequencies, the Kirchhoff method provides realistic acoustic directivities with the correct number of lobes (Fig. 5). For the lowest frequency (f = 1 kHz), the wavelength is not small compared to the surface’s radius of curvature (λ/R = 2) and the locally plane reflection assumption distorts the levels. In particular, the shadow effect in the half-space x < 0 is overestimated due to the incomplete modeling of refraction. The levels are correct only in a small angular sector in front of the source (−10° < θ < 10°). The accuracy of calculations by the Kirchhoff method logically improves with increasing frequency, i.e. with decreasing the ratio λ/R (Figs. 5b–5d). For the highest frequency (f = 7.5 kHz, λ/R = 0.27), the levels match the analytical solution throughout all the half-space x > 0 on the source side.

Figure 5 Directivity of the scattering of a point source by a sphere. (a) f = 1 kHz, λ/R = 2; (b) f = 2.5 kHz, λ/R = 0.8; (c) f = 5 kHz, λ/R = 0.4 and (d) f = 7.5 kHz, λ/R = 0.27.

Each of these calculations carried out for 11,968 reflective cells, 100 time steps, and 180 observation points, takes 12.3 s on a processor running at 3 GHz.

4 Comparisons with Euler calculations

In this section, reflection calculations by the Kirchhoff surface integral are carried out for an airfoil in a flow, as an application example. The results are compared with the numerical solution of the linearized Euler equations, in the presence of a mean flow.

The reference solution is a 2D Euler calculation in the frequency domain, based on a high-order discontinuous Galerkin (DG) method [24, 25]. The reflecting surface is a NACA 0012 profile of chord C = 1 m, without incidence in a flow at Mach number 0.5. The computational domain is a disk of radius 8 m around the profile. A noise source near the leading edge of the profile is modeled by a sinusoidal pressure disturbance at frequency f = 2 kHz, with a Gaussian distribution of width 0.02 m at the midpoint of the Gaussian. The center of the source is 0.1 m upstream of the leading edge of the profile, 0.2 m below. Two unstructured meshes made of triangles are designed, one without the profile, the other with it, to highlight the effect of the profile on the pressure field radiated by the acoustic source (Fig. 6). Particular care is taken to ensure that these two meshes are as similar as possible in the common part of the fluid domain. The number of cells is respectively 10,324 for the mesh without profile and 8580 for the mesh with profile. The mean flow is uniform in the case without profile. It is determined by solving Euler’s equations in the case with the profile. DG acoustic calculations are performed at order 8 in both cases. The harmonic solution of the density, velocity and pressure fields is projected onto 90 points distributed in a circle of radius 0.1 m around the source for FW-H free-field radiation calculations (the source term F of this formulation is recalled in Appendix B).

Figure 6 Zoom in the unstructured meshes used for the 2D DG calculations in free field (left) and with the NACA profile (right). Circle of input data for FW-H calculations (in red).

This 2D configuration is best modeled in the 3D FW-H and Kirchhoff calculations presented below, with a wing and a cylindrical “source” surface obtained by duplicating the profile and the “source” circle in the third direction. This duplication is carried out over 50 m on either side of the DG calculation plane, with a step of 0.1 m, i.e. a total of 1001 vertical planes (Fig. 7). For FW-H calculations of the source free-field radiation, the “source” surface is therefore described by a structured grid of 1001 × 90 nodes. On each of the 90 lines of this grid, the harmonic solution (amplitude and phase) provided by the DG calculation on the “source” circle, is converted into the time domain with a time step equal to 1/50 of the period T of the noise source. For reflection calculations using the Kirchhoff method, the wing is described by a structured double-layer grid, with 1001 points in span and 99 points in chord, many of which concentrated at the leading edge for a precise description of its geometry. The thickness e of the double layer on either side of the wing skin is 0.001 m (λ/e = 170) for an accurate calculation of the pressure gradients at its leading edge. The calculations are carried out for a speed of sound a₀ = 340 m/s.

Figure 7 Configuration of the reflection calculation by Kirchhoff method: control surface for free-field acoustic radiation using FW-H formulation (in red), reflecting surface (in green), and observation plane (in blue).

The calculation procedure is as follows. First, FW-H calculations of free-field radiation are performed for 109 × 90 observation points describing a 4 m long and 3 m high rectangle around the wing, in its median vertical plane (Fig. 7) and for the nodes of the wing’s double-layer mesh. Then, using the results on the wing, the Kirchhoff calculations of the pressure reflected by the wing is carried out for the observation points. Finally, the results of both calculations in the observation plane are summed. It should be noted that this procedure, combining FW-H calculations (with an integration surface surrounding the noise sources) and Kirchhoff calculations (to take into account the reflections of solid bodies outside the FW-H control surface), can be applied to most numerical aerodynamic simulations.

The results of the FW-H calculations are compared with the reference solution in Figures 8 and 9. The discrepancies between the two calculations are not discernible for instantaneous fields (Fig. 8). The agreement of the integrated levels (Fig. 9) confirms both the validity of the FW-H calculations starting from the DG fields on the “source” surface, and the sufficient span of the wing to find the 2D solution in the median plane. The slight loss of accuracy of the DG calculations in the bottom-left corner is due to the slightly less dense mesh in this area (Fig. 6). The small level oscillations in the FW-H results are due both to the rather coarse discretization of the “source” cylinder in span (Δz = 0.1 m) and to its finite length. Overall, these FW-H calculations appear to provide sufficiently accurate input data to the Kirchhoff reflection calculation to enable valid comparisons with the reference solution.

Figure 8 Instantaneous pressure radiated by the noise source in free-field (without the profile, direct pressure). DG calculation (left) and FW-H calculation (right).

Figure 9 Same comparison as above for the integrated pressure levels.

The results of the Kirchhoff reflection calculation are shown in Figure 10. It can be seen that the profile radiates reflected waves on the source side and waves in phase opposition with the incident pressure on the opposite side. Naturally, the integrated levels do not allow us to differentiate the constructive contributions of this radiation from the destructive ones.

Figure 10 Reflected pressure. Kirchhoff calculation. Instantaneous pressure (left) and integrated levels (right).

The total pressure, resulting from the sum of the direct and reflected pressures, is plotted in Figures 11 and 12. For the frequency studied, they show that the profile strongly modifies the acoustic field of the noise source, with fairly marked areas of amplification or silence. The integrated levels predicted by the simplified method are close to those obtained by Euler’s calculations. The two main deviations from the reference solution are: less pronounced directivity lobes upstream of the profile, and a less destructive interference at 45° below the profile. The shadow zone above the profile is close to the reference solution, despite the assumption of a large local radius of curvature compared to the wavelength. This is a very satisfactory result, given the small radius of curvature of the profile trailing edge. The deviations in the bottom-left corner, due to slight numerical dissipation in the DG calculation, show, conversely, the interest of Kirchhoff calculations for far-field acoustic predictions, especially when the mesh of aerodynamic simulation cannot be designed for acoustic propagation over long distances for reasons of computational cost.

Figure 11 Instantaneous pressure with the profile (total pressure). DG calculation (left) and FW-H calculation (right).

Figure 12 Same comparison as above for the integrated pressure levels.

The reflection calculation carried out with the Kirchhoff method for 99,000 reflective cells, 50 time steps, and 9810 observation points, take about 40 min on a processor running at 3 GHz.

5 Conclusion

Assuming locally plane reflections and a uniform mean flow parallel to the reflecting surface, a calculation of the reflection and shadow effects of a solid body (wing, fuselage, building, ground relief) in an acoustic field was obtained rather simply with a Kirchhoff surface integral. The calculation time of the reflected pressure field is about 5. 10⁻⁸ s/reflective cell/time step/obs on a processor running at 3 GHz. For example, for a reflective surface described by 20,000 cells and signals described by 5000 instants, the calculation time is about 5 s per observation point. In practice, a few hundred observation points are enough to characterize acoustic radiation in the far field. Calculating the reflected pressure field is therefore inexpensive in itself. In fact, what can be most costly is to determine the incident pressures on the body surface, but this input data is necessary whatever method is used.

Test cases show that the proposed model can provide realistic and useful results for a first analysis of acoustic scattering, for example as part of an optimization process. The results are all the more accurate as the assumption of short wavelength compared to the local radius of curvature of the reflective surface is satisfied. In the present study, the incident pressure gradient on the reflective surface was calculated by finite differences on two layers on either side of the surface, spaced by a prescribed value. When the incident pressure is calculated using the Kirchhoff or Ffowcs Williams and Hawkings formulation, it is preferable to determine this pressure gradient directly using the same integral formulation, with the corresponding source terms. This prevents the inaccuracies of the calculation by finite differences, especially for broadband noise sources, where the distance between the two layers cannot logically be optimal for all frequencies.

Apart from its computational cost, another strength of the method is that it can be applied to unsteady aerodynamic numerical simulations without any sound source model. It works directly starting the incident pressure and its gradient on the reflective surfaces and requires no information on the origin of the incident acoustic field. Furthermore, time-domain calculation offers a significant advantage for periodic signals rich in harmonics or broadband signals. This method is therefore particularly attractive for BVI or HSI noise sources of rotors or propellers, and for noise sources of wakes or jets (recent example of application to a jet in [26]). For all these applications, BEM methods in the frequency domain remain considerably more expensive (when applicable).

The drawbacks of the present approach are mainly its lack of accuracy for high wavelengths (with respect to the body’s radius of curvature), and its inability to handle sharp angles correctly. A possible improvement could be to use this high-frequency approximation in conjunction with a BEM method for the low frequencies. Once the spectral analysis of the results of the Kirchhoff calculation has been carried out, this would consist of replacing the imprecise levels at low frequencies with those obtained with the BEM calculations. It should also be possible to enrich the method to deal with sharp angles, although, from our point of view, at the cost of a level of complexity that would probably preclude interest of this approach for industrial applications.

Conflicts of interest

The authors declare no conflicts of interest.

Data availability statement

The data are available from the corresponding author on request.

Appendix

Expression of the Kirchhoff surface integral for a fixed integration surface and fixed observation points

The starting point of the formulation is the generalized Green’s formula derived by Goldstein ([3], pp. 59–61):

$$p\left(\vec{x},t\right)=\underset{\tau }{\overset{}{\int }}\underset{S}{\overset{}{\int }}\left[\left(p\frac{\partial G}{\partial n}-G\frac{\partial p}{\partial n}\right)+\frac{M_0{n}_1}{a_0}\left(G\frac{\mathrm{dp}}{\mathrm{d\tau }}-p\frac{\mathrm{dG}}{\mathrm{d\tau }}\right)\right]\mathrm{dSd\tau }$$

with $\frac{d}{\mathrm{d\tau }}=\frac{\partial }{\partial \tau }+U_0\frac{\partial }{\partial y_1}$ and the notations given in Section 2.

Developing the derivatives d/dτ gives,

$$p\left(\vec{x},t\right)=\underset{\tau }{\overset{}{\int }}\underset{S}{\overset{}{\int }}\left[p\left(n_i\frac{\partial G}{\partial y_i}-\frac{M_0{n}_1}{a_0}\frac{\partial G}{\partial \tau }-M_0^2{n}_1\frac{\partial G}{\partial y_1}\right)-G\left(\frac{\partial p}{\partial n}-\frac{M_0{n}_1}{a_0}\frac{\partial p}{\partial \tau }-M_0^2{n}_1\frac{\partial p}{\partial y_1}\right)\right]\mathrm{dSd\tau }$$

Knowing that, for a fixed surface S and a fixed observation point $\vec{x}$, $\partial g/\partial \tau =1$ and $\partial d/\partial \tau =0$,

$$\begin{array}{c}\underset{\tau }{\overset{}{\int }}p\frac{\partial G}{\partial \tau }\mathrm{d\tau }=\underset{g}{\overset{}{\int }}p\frac{\partial G}{\partial g}\mathrm{dg}=\underset{g}{\overset{}{\int }}\frac{p}{4\pi d}\frac{\partial \left(\delta \right(g\left)\right)}{\partial g}\mathrm{dg}\\ =-\underset{g}{\overset{}{\int }}\frac{\partial p}{\partial g}\frac{\delta \left(g\right)}{4\pi d}\mathrm{dg}=\underset{\tau }{\overset{}{\int }}\frac{\partial p}{\partial \tau }\mathrm{Gd\tau }\end{array}$$

Developing $\frac{\partial G}{\partial y_i}=\frac{1}{4\pi d}\left(\frac{\partial \left(\delta \right(g\left)\right)}{\partial g}\frac{\partial g}{\partial y_i}-\frac{\delta \left(g\right)}{d}\frac{\partial d}{\partial y_i}\right)$ gives:

$$\begin{array}{c}\underset{\tau }{\overset{}{\int }}p{n}_i\frac{\partial G}{\partial y_i}\mathrm{d\tau }=\underset{g}{\overset{}{\int }}\frac{p}{4\pi d}\frac{\partial g}{\partial y_i}{n}_i\frac{\partial \left(\delta \right(g\left)\right)}{\partial g}\mathrm{dg}-\underset{\tau }{\overset{}{\int }}\frac{p}{d}\frac{\partial d}{\partial y_i}{n}_i\frac{\delta \left(g\right)}{4\pi d}\mathrm{d\tau }\\ =-\underset{\tau }{\overset{}{\int }}\frac{\partial }{\partial \tau }\left(p\frac{\partial d}{\partial y_i}{n}_i\right)\frac{\delta \left(g\right)}{4\pi d}\mathrm{d\tau }-\underset{\tau }{\overset{}{\int }}\frac{p}{d}\frac{\partial d}{\partial y_i}{n}_i\frac{\delta \left(g\right)}{4\pi d}\mathrm{d\tau }\\ =-\underset{\tau }{\overset{}{\int }}\frac{\partial p}{\partial \tau }\frac{\partial g}{\partial y_i}{n}_i\mathrm{Gd\tau }-\underset{\tau }{\overset{}{\int }}\frac{p}{d}\frac{\partial d}{\partial y_i}{n}_i\mathrm{Gd\tau }\end{array}$$

$p(\vec{x},t)$ can therefore write: $p\left(\vec{x},t\right)=\underset{\tau }{\overset{}{\int }}\underset{S}{\overset{}{\int }}F\left(\vec{x},\vec{y},\tau \right)\mathrm{GdSd\tau }$

$$\begin{array}{c}\mathrm{with} F\left(\vec{x},\vec{y},\tau \right)=-\frac{\partial g}{\partial y_i}{n}_i\frac{\partial p}{\partial \tau }-\frac{\partial d}{\partial y_i}{n}_i\frac{p}{d}+\frac{M_0{n}_1}{a_0}\frac{\partial p}{\partial \tau }\\ +M_0^2\frac{\partial g}{\partial y_1}{n}_1\frac{\partial p}{\partial \tau }+M_0^2\frac{\partial d}{\partial y_1}{n}_1\frac{p}{d}\\ -\frac{\partial p}{\partial n}+\frac{M_0{n}_1}{a_0}\frac{\partial p}{\partial \tau }+M_0^2{n}_1\frac{\partial p}{\partial y_1}\end{array}$$

Developing $\frac{\partial d}{\partial y_i}=\frac{-C_i\left(x_i-y_i\right)}{d}$ and $\frac{\partial g}{\partial y_i}=\frac{1}{a_0{\beta }^2}\left(\frac{-C_i\left(x_i-y_i\right)}{d}+M_0{\delta }_{1i}\right)$ and rearranging terms eventually leads to:

$$F\left(\vec{x},\vec{y},\tau \right)=p\frac{{\beta }^2}{d^2}\left(x_i-y_i\right)n_i-\frac{\partial p}{\partial y_i}{n}_i+M_0^2\frac{\partial p}{\partial y_1}{n}_1+\frac{1}{a_0}\left(\frac{\left(x_i-y_i\right)n_i}{d}+M_0{n}_1\right)\frac{\partial p}{\partial \tau }$$

Appendix

Expression of pressure gradient

The integration surface and the observation points are assumed to be fixed in the calculation reference frame. The flow is uniform at velocity U₀ = a₀M₀ in the direction of the x axis.

General expression (see notations in Sect. 2)

Starting from the expression of the pressure for the Kirchhoff or FW-H formulation in the form [27]:

$$p\left(\vec{x},t\right)\underset{\tau }{\overset{}{\int }}\underset{S}{\overset{}{\int }}F\left(\vec{x},\vec{y},\tau \right)\mathrm{GdSd\tau } \hspace{0.5em}\mathrm{with} G=\frac{\delta \left(g\right)}{4\pi d}$$

the gradient can write:

$$\frac{\partial p\left(\vec{x},t\right)}{\partial x_j}=\underset{\tau }{\overset{}{\int }}\underset{S}{\overset{}{\int }}\frac{\partial F}{\partial x_j}\mathrm{GdSd\tau }+\underset{\tau }{\overset{}{\int }}\underset{S}{\overset{}{\int }}F\frac{\partial G}{\partial x_j}\mathrm{dSd\tau }$$

and developing $\frac{\partial G}{\partial x_j}$ gives: $\frac{\partial G}{\partial x_j}=\frac{1}{4\pi d}\left(\frac{\partial \left(\delta \right(g\left)\right)}{\partial g}\frac{\partial g}{\partial x_j}-\frac{\delta \left(g\right)}{d}\frac{\partial d}{\partial x_j}\right)$

For fixed integration surfaces and observation points ∂g/∂τ = 1 and thus:

$$\underset{\tau }{\overset{}{\int }}\frac{F}{4\pi d}\frac{\partial g}{\partial x_j}\frac{\partial \left(\delta \right(g\left)\right)}{\partial g}\mathrm{d\tau }=\underset{g}{\overset{}{\int }}\frac{F}{4\pi d}\frac{\partial g}{\partial x_j}\frac{\partial \left(\delta \right(g\left)\right)}{\partial g}\mathrm{dg}=\underset{g}{\overset{}{\int }}\frac{\partial }{\partial g}\left(\frac{\partial g}{\partial x_j}\frac{F}{4\pi d}\right)\delta \left(g\right)\mathrm{dg}=-\underset{\tau }{\overset{}{\int }}\frac{\partial g}{\partial x_j}\frac{\partial F}{\partial \tau }\frac{\delta \left(g\right)}{4\pi d}\mathrm{d\tau }$$

Therefore, the pressure gradient writes:

$$\frac{\partial p(\vec{x},t)}{\partial x_j}=\underset{\tau }{\overset{}{\int }}\underset{S}{\overset{}{\int }}\left(\frac{\partial F}{\partial x_j}-\frac{\partial g}{\partial x_j}\frac{\partial F}{\partial \tau }-\frac{\partial d}{\partial x_j}\frac{F}{d}\right)\mathrm{GdSd\tau }$$

Writing $\widehat{r}=\frac{x_i-y_i}{d}$:

$$\frac{\partial d}{\partial x_j}=C_j{\widehat{r}}_j \hspace{0.5em}\mathrm{and}\hspace{0.5em} \frac{\partial g}{\partial x_j}=\left(\frac{C_j{\widehat{r}}_j-M_0{\delta }_{1j}}{a_0{\beta }^2}\right).$$

Application to Kirchhoff formulation

The source term F is: $F= M_0^2{n}_1\frac{\partial p}{\partial y_1}-\frac{\partial p}{\partial n}+\frac{1}{a_0}\left(M_0{n}_1+\frac{\left(x_i-y_i\right)n_i}{d}\right)\frac{\partial p}{\partial \tau }+{\beta }^2\frac{\left(x_i-y_i\right)n_i}{d^2}p$

As the pressure p does not depend on x_j:

$$\frac{\partial F}{\partial x_j}=\frac{1}{a_0}\left(\frac{n_j}{d}-\frac{\left(x_i-y_i\right)n_i}{d^2}\frac{\partial d}{\partial x_j}\right)\frac{\partial p}{\partial \tau }+{\beta }^2\left(\frac{n_j}{d^2}-\frac{2\left(x_i-y_i\right)n_i}{d^3}\frac{\partial d}{\partial x_j}\right)p$$

Writing ${\widehat{r}}_n={\widehat{r}}_i{n}_i$:

$$\begin{array}{c}\frac{\partial F}{\partial x_j}=\frac{1}{a_{0}d}\left(n_j-{\widehat{r}}_n{C}_j{\widehat{r}}_j\right)\frac{\partial p}{\partial \tau }+\frac{{\beta }^2}{d^2}\left(n_j-{2\widehat{r}}_n{C}_j{\widehat{r}}_j\right)p\\ \frac{\partial g}{\partial x_j}\frac{\partial F}{\partial \tau }=\frac{\left(C_j{\widehat{r}}_j-M_0{\delta }_{1,j}\right)}{a_0{\beta }^2}\times \left[M_0^2{n}_1\frac{{\partial }^{2}p}{\partial {\tau \partial y}_1}-\frac{{\partial }^{2}p}{\partial \tau \partial n}+\frac{1}{a_0}\left(M_0{n}_1+{\widehat{r}}_n\right)\frac{{\partial }^{2}p}{\partial {\tau }^2}+{\beta }^2\frac{{\widehat{r}}_n}{d}\frac{\partial p}{\partial \tau }\right]\\ \frac{\partial d}{\partial x_j}\frac{F}{d}=\frac{C_j{\widehat{r}}_j}{d}\left[M_0^2{n}_1\frac{\partial p}{{\partial y}_1}-\frac{\partial p}{\partial n}+\frac{1}{a_0}\left(M_0{n}_1+{\widehat{r}}_n\right)\frac{\partial p}{\partial \tau }+{\beta }^2\frac{{\widehat{r}}_n}{d}p\right]\end{array}$$

Grouping terms finally leads to the following expression of the pressure gradient:

$$\begin{array}{c}{\alpha }_{1j}=\frac{{\beta }^2}{d^2}\left(n_j-3{\widehat{r}}_n{C}_j{\widehat{r}}_j\right)\\ {\alpha }_{2j}={-M}_0^2{n}_1\frac{C_j{\widehat{r}}_j}{d}\\ \begin{array}{c}{\alpha }_{3j}=\frac{C_j{\widehat{r}}_j}{d}\\ {\alpha }_{4j}=\frac{1}{a_{0}d}\left(n_j-\left(M_0{n}_1+3{\widehat{r}}_n\right)C_j{\widehat{r}}_j+M_0{\widehat{r}}_n{\delta }_{1j}\right)\\ \begin{array}{c}{\alpha }_{5j}=-\frac{M_0^2{n}_1}{a_0{\beta }^2}\left(C_j{\widehat{r}}_j-M_0{\delta }_{1j}\right)\\ {\alpha }_{6j}=\frac{1}{a_0{\beta }^2}\left(C_j{\widehat{r}}_j-M_0{\delta }_{1j}\right)\\ {\alpha }_{7j}=-\frac{1}{a_0^2{\beta }^2}\left(C_j{\widehat{r}}_j-M_0{\delta }_{1j}\right)\left(M_0{n}_1+{\widehat{r}}_n\right)\end{array}\end{array}\end{array}$$

In the far field, the integration can be limited to the last three terms.

Application to Ffowcs Williams and Hawkings formulation

The source term F is: $F= \frac{{\beta }^2}{d^2}{\mathrm{\Sigma }}_1+\frac{\partial }{\partial \tau }\left(\frac{{\mathrm{\Sigma }}_1}{a_{0}d}+\frac{{\mathrm{\Sigma }}_2}{{\beta }^2}\right)$

$$\begin{array}{c}{\mathrm{\Sigma }}_1=\left(A_i\left(x_i-y_i\right)+\frac{\left(M_0^2{A}_1-U_{0}B\right)}{{\beta }^2}\left(x_1-y_1\right)\right) \mathrm{and} {\mathrm{\Sigma }}_2=B-\frac{M_0{A}_1}{a_0}\\ A_i=p{n}_i+p{u}_i\left(u_n+U_0{n}_1\right) \mathrm{and} B=-{\rho }_0{U_{0}n}_1+\rho \left(u_n+{U_{0}n}_1\right)\end{array}$$

where $\vec{u}$ is the velocity of the fluid in a reference frame in which the unperturbed fluid is at rest, and $u_n=\vec{u}·\vec{n}$

$$\begin{array}{c}\frac{\partial F}{\partial x_j}=\frac{{\beta }^2}{d^2}\left(\frac{\partial {\mathrm{\Sigma }}_1}{\partial x_j}-\frac{2{\mathrm{\Sigma }}_1}{d}{C}_j{\widehat{r}}_j\right)+\frac{\partial }{\partial \tau }\left(\frac{1}{a_{0}d}\frac{\partial {\mathrm{\Sigma }}_1}{\partial x_j}-\frac{{\mathrm{\Sigma }}_1}{a_0{d}^2}{C}_j{\widehat{r}}_j+\frac{1}{{\beta }^2}\frac{\partial {\mathrm{\Sigma }}_2}{\partial x_j}\right)\\ \frac{\partial g}{\partial x_j}\frac{\partial F}{\partial \tau }=\left(\frac{C_j{\widehat{r}}_j-M_0{\delta }_{1j}}{a_0{\beta }^2}\right)\left[\frac{{\beta }^2}{d^2}\frac{\partial {\mathrm{\Sigma }}_1}{\partial \tau }+\frac{{\partial }^2}{\partial {\tau }^2}\left(\frac{{\mathrm{\Sigma }}_1}{a_{0}d}+\frac{{\mathrm{\Sigma }}_2}{{\beta }^2}\right)\right]\\ \frac{\partial d}{\partial x_j}\frac{F}{d}=C_j{\widehat{r}}_j\left[\frac{{\beta }^2}{d^3}{\mathrm{\Sigma }}_1+\frac{\partial }{\partial \tau }\left(\frac{{\mathrm{\Sigma }}_1}{a_0{d}^2}+\frac{{\mathrm{\Sigma }}_2}{{\beta }^{2}d}\right)\right]\end{array}$$

As A_i et B do not depend on x_j: $\frac{\partial {\mathrm{\Sigma }}_2}{\partial x_j}=0$.

Grouping terms and writing ${\mathrm{\Sigma }}_j^{\text{'}}=\frac{\partial {\mathrm{\Sigma }}_1}{\partial x_j}-3C_j{\widehat{r}}_j\frac{{\mathrm{\Sigma }}_1}{d}=A_j+\frac{\left(M_0^2{A}_1-U_{0}B\right)}{{\beta }^2}{\delta }_{1j}-3C_j{\widehat{r}}_j\frac{{\mathrm{\Sigma }}_1}{d}$ finally leads to the following expression of the pressure gradient:

with:

$$\begin{array}{c}{F}_{1j}={\beta }^2\frac{{\mathrm{\Sigma }}_j^{\text{'}}}{d^2}\\ F_{2j}=\frac{{\mathrm{\Sigma }}_j^{\text{'}}}{a_{0}d}+M_0{\delta }_{1j}\frac{{\mathrm{\Sigma }}_1}{a_0{d}^2}-\frac{C_j{\widehat{r}}_j}{{\beta }^2}\frac{{\mathrm{\Sigma }}_2}{d}\\ F_{3j}=-\frac{\left(C_j{\widehat{r}}_j-M_0{\delta }_{1j}\right)}{{\beta }^2}\left(\frac{{\mathrm{\Sigma }}_1}{a_0^{2}d}+\frac{{\mathrm{\Sigma }}_2}{a_0{\beta }^2}\right)\end{array}$$

In the far field, the integration can be limited to the third term.

References

All Tables

All Figures

	Figure 1 Schematic of acoustic waves towards the solid body and reflection towards the observation point.
In the text

	Figure 2 Illustration of the effect of propagation times on the duration of the valid total signal (case of non-periodic signals).
In the text

	Figure 3 Calculation of the scattering of a point source by a rigid sphere.
In the text

	Figure 4 Influence of the mesh density on the Kirchhoff results (f = 7.5 kHz).
In the text

	Figure 5 Directivity of the scattering of a point source by a sphere. (a) f = 1 kHz, λ/R = 2; (b) f = 2.5 kHz, λ/R = 0.8; (c) f = 5 kHz, λ/R = 0.4 and (d) f = 7.5 kHz, λ/R = 0.27.
In the text

	Figure 6 Zoom in the unstructured meshes used for the 2D DG calculations in free field (left) and with the NACA profile (right). Circle of input data for FW-H calculations (in red).
In the text

	Figure 7 Configuration of the reflection calculation by Kirchhoff method: control surface for free-field acoustic radiation using FW-H formulation (in red), reflecting surface (in green), and observation plane (in blue).
In the text

	Figure 8 Instantaneous pressure radiated by the noise source in free-field (without the profile, direct pressure). DG calculation (left) and FW-H calculation (right).
In the text

	Figure 9 Same comparison as above for the integrated pressure levels.
In the text

	Figure 10 Reflected pressure. Kirchhoff calculation. Instantaneous pressure (left) and integrated levels (right).
In the text

	Figure 11 Instantaneous pressure with the profile (total pressure). DG calculation (left) and FW-H calculation (right).
In the text

	Figure 12 Same comparison as above for the integrated pressure levels.
In the text