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

1 Introduction

Acoustic measurements for aircraft certification must be acquired by microphones placed at 1.2 m above the ground, as specified by the International Civil Aircraft Organization in the Annex 16 Environmental Protection – Aircraft Noise. Because of the source motion and interference between the direct wave and reflected waves at the ground, this geometry leads to a time-dependent comb filtering on the noise measurements [1]. The interference pattern is significantly sensitive to ground effects, meteorological effects, and source trajectory. As a result, large uncertainties are associated with sound pressure levels (SPL) measurements, particularly due to single tones such as tonal components emitted by the aircraft engines [2]. Ongoing efforts are thus being made by the Committee on Aviation Environmental Protection to establish enhanced noise measurement methods for certification of civil aircraft (see, e.g. [2–5]).

The use of sound propagation models is promising to help design new standards for aircraft certification. Analytical solutions for moving sources above an absorbing ground surface have been obtained in very simple cases, i.e. for a monopole source following a rectilinear motion at constant speed and at constant altitude above a flat ground in a homogeneous atmosphere [6–9]. In particular, an asymptotic formulation in far-field and at grazing angle, called the “Dopplerized” Weyl Van der Pol formulation, has been proposed [10]. For the case of an inhomogeneous atmosphere in motion, it is necessary to rely on numerical simulations. Regarding outdoor acoustic propagation, numerous methods are available in the literature. They can be sorted into geometrical approaches, such as ray-tracing [11, 12] or Gaussian beam methods [13], and wave-based approaches, that solve either paraxial wave equations [14, 15] or linearized fluid mechanics equations [16, 17]. In particular, time-domain approaches [18, 19] are well suited for studying radiation of moving sources, as they can consider arbitrary trajectories.

As the aircraft noise impact studies need to be carried out over large scale, the current models employed in the industry are often based on faster and simplified approaches. They use either analytical models or geometrical approaches. The former usually assume a homogeneous atmosphere and a quasi-stationary source, neglecting the difference between the source position at the emission and reception times. These assumptions are however invalid for sources traveling at a non-negligible Mach number in an outdoor environment, which is a characteristic scenario of aircraft noise. Enhancing current prediction methods is crucial to simulate reliable aircraft noise predictions and to propose new experimental procedures. Recent works focused on the use of ray-tracing (RT) algorithms (e.g. [20–23]) as it is an efficient sound propagation method that can take into account meteorological effects. The RT methods however present some limitations, such as the presence of caustics, as well as shadow zones where rays can not propagate. In addition, the inclusion of source motion effects (convective amplification, Doppler effect, consideration of emission time) on sound propagation is particularly tricky in RT methods. Moreover, validation of the existing approaches often rely on test cases for a nonmoving source.

In order to improve the prediction of the sound pressure levels due to a moving source for application to aircraft noise, this work proposes an innovative methodology based on the coupling of a heuristic formulation for a monopole in arbitrary motion [10, 24] with a RT model [11]. The heuristic formulation can be seen as an extension of the analytical “Dopplerized” Weyl Van der Pol formulation. It accounts for source motion effects, as well as ground effects with the ground surface admittance determined at the Doppler frequency [25, 26]. The RT model accounts for refraction and is used to determine meteorological-dependent parameters, such as curved ray length, incidence angle for reflection, or propagation time, that are used as inputs for the heuristic model. Validation of the coupled model is performed for several test cases with different source trajectories and atmospheric conditions using a reference solution obtained from a three-dimensional (3D) Finite Difference Time Domain (FDTD) solver of the linearized Euler equations [19].

We first review analytical models for moving sources above an absorbing ground and present the coupled RT heuristic model in Section 2. The validation of the model is presented in Section 3. Then, the SPL variability due to the ground properties and meteorological conditions is illustrated in the context of aircraft noise in Section 4. Finally, concluding remarks are given in Section 5.

2 Sound field formulation for a monopole in motion

Let us consider sound radiation by a harmonic monopole moving above a flat ground in an inhomogeneous and moving atmosphere. In the Cartesian coordinate system x = (x, y, z), the source trajectory is denoted x_s(t) = (x_s(t), y_s(t), z_s(t)) with t the time. The equations governing sound propagation are the reduced linearized Euler equations [27]:

$$\frac{\partial p}{\partial t}+{\mathbf{V}}_{\mathbf{0}}·\nabla p+{\rho }_0{c}_0^2\nabla ·\mathbf{v}=c_0^2{S}_0 {\mathrm{e}}^{-\mathrm{i}{\omega }_{0}t}\delta [\mathbf{x}-{\mathbf{x}}_s(t\left)\right],$$(1)

$$\frac{\partial \mathbf{v}}{\partial t}+{\mathbf{V}}_{\mathbf{0}}·\nabla \mathbf{v}+\frac{1}{{\rho }_0}\nabla p+\mathbf{v}·\nabla {\mathbf{V}}_{\mathbf{0}}=0.$$(2)

where p and v are respectively the acoustic pressure and velocity, ρ₀ is the air density, c₀ is the sound speed, V₀ is the base flow velocity, and ω₀ = 2πf₀ with f₀ the source frequency. The source amplitude S₀ is set to 1 kg s⁻¹ and is omitted in the following for simplicity. The boundary condition at z = 0 is given by:

$${\rho }_0{c}_0{v}_z+\left[\tilde{\beta }\mathrm{*}p\right]=0,$$(3)

where $\tilde{\beta }$ is the inverse Fourier transform of the ground normalized surface admittance β(ω), and the symbol * denotes the convolution operator.

From the source trajectory, we define the source speed vector V_s = dx_s/dt, the Mach number vector M_s = V_s/c₀, the Mach number M_s = |M_s|, and the time derivative of the Mach number vector ${\dot{\mathbf{M}}}_{\mathbf{s}}=\mathrm{d}{\mathbf{M}}_{\mathbf{s}}/\mathrm{d}t$. We introduce also the Mach number vector associated to the mean flow as M₀ = V₀/c₀. Finally, we define the source-receiver distance at the time t as R(t) = |R(t)| with R(t) = x − x_s(t).

Unless indicated, the mean flow, the air density and sound speed can vary spatially. For the test cases in Section 3 and the application in Section 4, the air density and sound speed are constant and are equal to ρ₀ = 1.2 kg m⁻³ and c₀ = 340 m s⁻¹. In addition, atmospheric absorption is not considered in Sections 2 and 3 for simplicity, but is accounted for in the application in Section 4.

2.1 Analytical solution for a monopole moving at constant speed and height in homogeneous atmosphere

An analytical solution for the pressure field is known for a source moving at constant speed and height above an absorbing ground in a homogeneous atmosphere at rest [9]. The sound speed c₀ is constant in this section and V₀ = 0. Without loss of generality, the source position is chosen as x_s(t) = (M_sc₀t, 0, z_s), where M_s is constant. Figure 1 illustrates the geometry considered. Note that the subscripts 1 and 2 refer respectively to the direct and reflected waves. The angles φ and θ are respectively the elevation angle and the azimuthal angle of the sound wave propagation path. The acoustic path length at the corresponding emission time t_e (Eqs. (11) and (13)) is denoted by R_e (Eq. (12)).

Figure 1 Geometry considered for a moving monopole in rectilinear motion above the ground.

We use in the following an approximate solution in far-field [24, 28, 29], which is often called the “Dopplerized” Weyl Van der Pol formula. The pressure field is written as the sum of a direct wave and a reflected wave:

$$p(\mathbf{x},t)=-\frac{\mathrm{i}{\omega }_0{\mathrm{e}}^{-\mathrm{i}{\omega }_{0}t}}{4\pi }\left[\frac{{\mathrm{e}}^{\mathrm{i}{k}_0{R}_{\mathrm{e}1}}}{R_{e1}(1-M_{{\mathrm{r}}_1}{)}^2}+\mathcal{Q}\frac{{\mathrm{e}}^{\mathrm{i}{k}_0{R}_{\mathrm{e}2}}}{R_{e2}(1-M_{r_2}{)}^2}\right],$$(4)

where k₀ = ω₀/c₀ is the wave number and

$$M_{r_1}={\mathbf{M}}_s·{\mathbf{R}}_{e1}/R_{e1}=M_s\cos {\phi }_1\cos {\theta }_1,$$(5)

$$M_{r_2}={\mathbf{M}}_s·{\mathbf{R}}_{e2}/R_{e2}=M_s\cos {\phi }_2\cos {\theta }_2,$$(6)

are the projections of the Mach number vector in the source-receiver directions. Even if the source signal is monochromatic, the pressure at the receiver is not monochromatic anymore as R_e1 and R_e2 are varying with time due to source motion. However, an instantaneous frequency at the receiver can be obtained by performing a Taylor expansion of the phase around a reference time within a characteristic period of oscillation (see, e.g. Ref. [30] for details). For the direct arrival, this yields an instantaneous frequency at observer ${\omega }_{e_1}={\omega }_0/(1-M_{r_1})$. Similarly, the frequency associated to the reflected wave is ${\omega }_{e_2}={\omega }_0/(1-M_{r_2})$. It differs from that associated to the direct wave ${\omega }_{e_1}$, which implies that the direct and reflected waves are not perceived at the same frequency at the receiver.

The reflected wave is weighted by the spherical wave reflection coefficient:

$$\mathcal{Q}=\mathcal{R}+(1-\mathcal{R})F\left(d\right),$$(7)

where

$$F\left(u\right)=1+\mathrm{i}u\sqrt{\pi }w\left(u\right)$$(8)

is the boundary loss factor, w is the Faddeeva function,

$$\mathcal{R}=\frac{\sin {\phi }_2-\beta \left({\omega }_{e_2}\right)}{\sin {\phi }_2+\beta \left({\omega }_{e_2}\right)},$$(9)

is the plane wave reflection coefficient, and

$$d=\frac{1}{2}(1+\mathrm{i})\sqrt{\frac{k_0{R}_{e2}}{1-M_{r_2}}}[\sin {\phi }_2+\beta ({\omega }_{e_2}\left)\right],$$(10)

is the numerical distance. Note that the ground surface admittance is evaluated at the Doppler frequency associated to the reflected wave ${\omega }_{e_2}$ and not at the source frequency. In Reference [18], it was shown that evaluating the ground admittance at the source frequency ω₀ instead of the Doppler frequency ${\omega }_{e_2}$ can lead to a significant error, especially for grounds modelled as a hard-backed layer of a porous material, for which the ground admittance can show large variations with the frequency.

In addition, the amplification factors (1 − M_r)⁻² for the direct and reflected waves are induced by the source motion, which is referred as to the convective amplification effect in the literature (see, e.g. [31]). The exponent in the convective amplification factor depends on the type of source considered: it is two for monopoles and dipoles but it is three for quadripoles.

The geometric quantities in equation (4) must be expressed as a function of the emission time t_e as the sound received by an observer at time t has been emitted at time t_e:

$$t_e=t-R_e/c_0,$$(11)

where R_e = |R_e| = R(t_e) is the source-receiver distance at the emission time and R_e = x − x_s(t_e). One has:

$$R_e=\sqrt{(x-M_s{c}_0{t}_e{)}^2+y^2+(z\pm z_s{)}^2}.$$(12)

where the sign − has to be chosen for the direct wave and the sign + for the reflected wave. Introducing equation (12) into equation (11), the emission time is then determined by solving a polynomial equation of degree two in t_e. The physical solution leads to:

$$t_e=\frac{c_{0}t-M_{s}x-\sqrt{x^2+(1-M_s^2)[y^2+(z\pm z_s{)}^2]}}{c_0(1-M_s^2)}.$$(13)

Note that, for a given value of t, the emission time is different for the direct and the reflected wave. Finally, equation (13) is used to evaluate the position of the direct and image sources at the emission time and determine the geometric quantities appearing in equation (4).

2.2 Heuristic formulation for a monopole in arbitrary motion in homogeneous atmosphere

In the context of aircraft noise, we need to consider a source in arbitrary motion, as sketched in Figure 2. An analytical solution is available for a monopole in arbitrary motion in free-field (see, e.g. [32]), but not above an absorbing ground.

Figure 2 Geometry considered for a moving monopole in arbitrary motion above the ground.

Nevertheless, an heuristic formulation that extends equation (4) has been proposed by Attenborough and Van Renterghem [10] for this case:

$$p=-\frac{\mathrm{i}{\omega }_0{\mathrm{e}}^{-\mathrm{i}{\omega }_{0}t}}{4\pi }\left[\left(1-M_{r_1}+\frac{\mathrm{i}{\dot{M}}_{r_1}}{{\omega }_0}\right)\frac{{\mathrm{e}}^{\mathrm{i}{k}_0{R}_{e1}}}{R_{e1}(1-M_{r_1}{)}^3}+\mathcal{Q}\left(1-M_{r_2}+\frac{\mathrm{i}{\dot{M}}_{r_2}}{{\omega }_0}\right)\frac{{\mathrm{e}}^{\mathrm{i}{k}_0{R}_{e2}}}{R_{e2}(1-M_{r_2}{)}^3}\right].$$(14)

In this formulation, M_r = M_s · R_e/R_e is still the projection of the Mach number vector in the source-receiver direction at the emission time. Similarly, ${\dot{M}}_r={\dot{\mathbf{M}}}_{\mathbf{s}}·{\mathbf{R}}_{\mathrm{e}}/R_e$ is the projection of the vector ${\dot{\mathbf{M}}}_{\mathbf{s}}$ in the source-receiver direction and is related to the acceleration of the source. Attention has to be paid to the z-component of the Mach number vector associated to the image source, as the image source travels along the z-direction in the opposite direction to the direct source. In addition, the reflection coefficient has the same expression than in Section 2.1.

For an arbitrary motion, there is no explicit formula giving the emission time t_e as a function of t, contrary to the case of rectilinear uniform motion. Thus, the emission time must be determined numerically. To do so, for each value of t, one solves the equation t_e − t + R(t_e)/c₀ = 0. From the emission time, the geometric quantities in equation (14) can be determined.

It is noteworthy that ${\dot{M}}_{\mathrm{r}}$ becomes negligible when it is small compared to ω₀(1 − M_r). Denoting by τ_a a characteristic time of the source acceleration, this occurs if M ≪ ω₀τ_a(1 − M). Therefore, except for high subsonic source motion, for source trajectories with rapid variations, or for very low frequencies, the term ${\dot{M}}_r$ can be safely removed (see details in [30], Sect. 5.1.2). This is especially expected for scenarios of aircraft noise certification. We observed that the ${\dot{M}}_r$ term has a negligible influence on the various test cases outlined in the paper. When ${\dot{M}}_r$ is removed from the solution, equation (14) has the same form than equation (4).

2.3 Heuristic formulation for a monopole in arbitrary motion in inhomogeneous atmosphere

The above formulations consider propagation paths as straight lines between source and receiver. In order to take into account the refraction of acoustic waves due to meteorological effects, we use a ray-tracing approach [11, 33] to calculate the parameters that depend on meteorological conditions (e.g. ray path, propagation time, angle of incidence of the ray on the ground, etc.). These parameters are then used as input directly in the heuristic formulation described in equation (15). Figure 3 presents the geometry considered.

Figure 3 Geometry considered for a moving monopole in arbitrary motion in a inhomogeneous atmosphere above a ground.

For simplicity, we consider thereafter that there are only two eigenrays reaching the receiver during the motion of the source. We also assume that the sound rays do not pass through any caustic. We propose the following heuristic formula for the pressure field, which could be extended to multipath propagation:

$$p(\mathbf{x},t)=-\frac{\mathrm{i}{\omega }_0 {\mathrm{e}}^{-\mathrm{i}{\omega }_{0}t}}{4\pi }\left[A_1{\mathrm{e}}^{\mathrm{i}{\omega }_0{\tau }_1}+\mathcal{Q}{A}_2{\mathrm{e}}^{\mathrm{i}{\omega }_0{\tau }_2}\right].$$(15)

with τ the travel time from the source at the emission time to the receiver, i.e. τ = t − t_e, and A the amplitude, that accounts for propagation effects and convection effects, due to the wind and the source motion. In this formulation, the plane-wave reflection coefficient and the numerical distance d involved in the spherical wave reflection coefficient $\mathcal{Q}$ are expressed as:

$$\mathcal{R}=\frac{\sin {\phi }_g-\beta \left({\omega }_{e_2}\right)}{\sin {\phi }_g+\beta \left({\omega }_{e_2}\right)},$$(16)

$$d=\frac{1}{2}(1+\mathrm{i})\sqrt{{\omega }_{e_2}{\tau }_2}[\sin {\phi }_g+\beta ({\omega }_{e_2}\left)\right],$$(17)

with ${\omega }_{e_2}$ the Doppler frequency associated to the reflected wave (see Eq. (20)), and φ_g the angle of incidence of the ray on the ground. It should be remarked that φ_g differs from φ₂ as long as the reflected ray path is curved. Without source motion, the expressions for the plane-wave reflection coefficient and the numerical distance correspond to those reported in Reference [10], Eq. (11.7) for an inhomogeneous and moving atmosphere. In addition, it can be shown that equation (15) reduces to equation (14) in the case of a homogeneous atmosphere at rest (neglecting ${\dot{M}}_r$).

Denoting by x_r the ray position and by n the local unit vector normal to the wavefront, the six ray equations are written as (see, e.g. [11, 12, 33]):

$$\frac{\mathrm{d}{\mathbf{x}}_r}{\mathrm{d}{t}_r}=c_0\mathbf{n}+{\mathbf{V}}_{\mathbf{0}}$$(18)

$$\frac{\mathrm{d}\mathbf{n}}{\mathrm{d}{t}_r}=(\mathbf{q}·\mathbf{n})\mathbf{n}-\mathbf{q},$$(19)

with t_r the travel time along the ray and q = ∇c₀ + (∇V₀) · n, where (∇V₀) · n is the directional derivative of the flow velocity vector. To avoid any ambiguity, the ray equations are given with Einstein notation in Appendix A. The first three equations indicate that the ray travels at the group velocity, which accounts for the propagation at the local sound speed and the convection at the local wind speed. The other three equations are associated to the curvature of the ray by the gradients of the sound speed and wind speed. At t_r = 0, the eigenray is located at the source, i.e. x_r(t_r = 0) = x_s(t_e) and at t_r = τ, the eigenray reaches the receiver, i.e. x_r(t_r = τ) = x. The ray equations are solved using the standard fourth order Runge-Kutta algorithm.

Note that the Doppler frequency ω_e for an inhomogeneous medium in motion is defined as [30], Eq. (5.68):

$${\omega }_e={\omega }_0\frac{1+{\mathbf{M}}_{\mathbf{0}}\left(0\right)·\mathbf{n}\left(0\right)}{1+\left[{\mathbf{M}}_{\mathbf{0}}\right(0)-{\mathbf{M}}_{\mathbf{s}}(0\left)\right]·\mathbf{n}\left(0\right)},$$(20)

where n(0) has to be understood as the unit vector normal to the wavefront of the emitted wave along the ray path connecting the source and the receiver, i.e. the eigenray, M₀(0) is the Mach number vector associated to the flow at the source position at the emission time and M_s(0) is the source Mach number vector at the emission time.

We have also sinφ_g = −n_g · e_n, with n_g is the unit vector normal to the incident wavefront on the ground and e_n the unit vector normal to the ground surface directed towards the air (for a flat ground, e_n = [0, 0, 1]).

As we consider a monopole, the rays are parameterized using two angles φ and θ, corresponding respectively to the elevation and azimuthal angle. The amplitude along a ray is determined by calculating the ray tube cross-section S(t_r). To do so, we consider the variations of the ray position with respect to the angles θ and φ; denoting x_φ = dx_r/dφ and x_θ = dx_r/dθ, one has S(t_r) = |x_φ × x_θ|. The evolution of x_φ and x_θ along the ray requires to solve 12 additional ordinary differential equations. For completeness, these equations are presented in Appendix A. Finally, the amplitude at the receiver is calculated with:

$$A\left(\tau \right)=\sqrt{\frac{{\rho }_0\left(\tau \right)c_0\left(\tau \right)E(\phi ,\theta )}{S\left(\tau \right)}}\frac{1}{1+{\mathbf{M}}_{\mathbf{0}}\left(\tau \right)·\mathbf{n}\left(\tau \right)}$$(21)

with E(φ, θ) the energy flux, which has a constant value along each ray but depends on the ray launch angles. Its expression has been obtained by matching the ray solution with an analytical solution for a monopole moving at constant speed in a homogeneous moving atmosphere. This gives:

$$E\left(\phi ,\theta \right)=\frac{\left|\cos \phi \right|}{{\rho }_0\left(0\right)c_0\left(0\right)}\frac{[1+{\mathbf{M}}_{\mathbf{0}}(0)·\mathbf{n}(0){]}^2}{[1+({\mathbf{M}}_{\mathbf{0}}\left(0\right)-{\mathbf{M}}_{\mathbf{s}}\left(0\right))·\mathbf{n}(0){]}^4}$$(22)

Details on the analytical derivation are given in Appendix B. Note that the effect of source acceleration on the source directivity is not accounted for in (22).

The determination of the rays and of the amplitude at the receiver is done in several steps. First, a set of emission time values t_e is prescribed. The position of the source at these emission times is obtained from the source trajectory. Second, the eigenrays connecting the source and the receiver are sought. To do so, a shooting method is employed. To avoid launching a large number of rays for each position of the source and thus reduce the computational cost, the rays are launched from the receiver using the reciprocity principle; from the flow reversal theorem, this requires to reverse the flow in the ray equations, i.e. to replace V₀ by −V₀ in equations (A.18) and (A.19). A coarse discretization of the ray parameters θ and φ is used. Once the ray trajectory has been computed, we determine the minimal distance d_min between each source position and the ray. This gives a map d_min(θ, φ) for each source position. We then determine the local minima of d_min(θ, φ). If the minimum is smaller than a threshold set to 1 m, we consider that the minimum is associated to an eigenray. The identification of eigenrays is repeated for each source position. In a third step, we optimize the launch angles (θ, φ) for each eigenray in order to reduce the minimal distance between the source and the receiver. We employ for that the MATLAB function “fminsearch”, using a reduced threshold of 0.01 m. In the last step, the pressure amplitude along the eigenrays is calculated. For that, we launch again the eigenrays from the source to the receiver for each source position solving the full set of the 18 ray equations. Finally, the amplitude is calculated with equation (21) and is used with the travel time to determine the pressure field at the receiver with equation (15).

3 Validation of the heuristic formulation with a FDTD numerical solver

This section presents four test cases for the validation of the formulations in equations (4), (14), and (15) with a reference finite-difference time-domain solution of the linearized Euler equations (Eqs. (1) and (2)). The validation is performed at several space scale and for several Mach number. However, the maximum propagation distance is of the order of 100 m given that the 3D FDTD simulations are costly in terms of memory and calculation duration. In the following, the ground admittance is specified using the slit-pore model [34] for a porosity of 0.38 and an effective flow resistivity of 514 kPa s m⁻² which correspond to a bare type of ground.

3.1 Principles of the FDTD numerical solver

The linearized Euler equations in equations (1) and (2) are solved using FDTD methods. More specifically, optimized fourth order finite-difference schemes [35] are used. The numerical methods employed in the code are detailed in Reference [19] and recent updates are described in Reference [36]. A broadband impedance boundary condition is implemented at the ground [37]. The sides and the top of the domain feature a PML layer [38] to prevent unwanted reflections from the limits of the domain.

Instead of the Dirac delta function in the right-hand side of equation (1), the source in the FDTD model has a Gaussian spatial distribution:

$$S\left(\mathbf{x}\right)=\frac{1}{{\pi }^{3/2}{B}^3}\exp \left(-\frac{|\mathbf{x}{|}^2}{B^2}\right),$$(23)

with B the Gaussian width. Note that S(x) tends to the Dirac delta function as B tends to zero. As discussed in Reference [18], the directivity of the Gaussian source is modified by its motion. In free-field and in a homogeneous and quiescent atmosphere, the acoustic pressure for a moving source in uniform rectilinear motion with a Gaussian spatial support, denoted p_B, is related to that with a Dirac delta function, denoted p₀, by:

$$p_B(\mathbf{x},t)=\exp \left(-\frac{{\omega }_{\mathrm{e}}^2{B}^2}{4c_0^2}\right)p_0(\mathbf{x},t).$$(24)

As a Gaussian source in motion does not behave like a monopole, the directivity has to be taken into account for comparison with the heuristic solution. To do so, comparison with the FDTD solution is done for a modified heuristic solution, in which the direct and reflected waves are multiplied by the directivity of the Gaussian source. For instance, for a homogeneous atmosphere at rest and for a source in arbitrary motion, equation (14) becomes:

$$p(\mathbf{x},t)=-\frac{\mathrm{i}{\omega }_0{\mathrm{e}}^{-\mathrm{i}{\omega }_{0}t}}{4\pi }[\exp \left(-\frac{k_0^2{B}^2}{4(1-M_{r_1}{)}^2}\right)\left(1-M_{r_1}+\frac{\mathrm{i}{\dot{M}}_{r_1}}{{\omega }_0}\right)\frac{{\mathrm{e}}^{\mathrm{i}{k}_0{R}_{e1}}}{R_{e1}(1-M_{r_1}{)}^3}+\mathcal{Q} \exp \left(-\frac{k_0^2{B}^2}{4(1-M_{r_2}{)}^2}\right)\left(1-M_{r_2}+\frac{\mathrm{i}{\dot{M}}_{r_2}}{{\omega }_0}\right)\frac{{\mathrm{e}}^{\mathrm{i}{k}_0{R}_{e2}}}{R_{e2}(1-M_{r_2}{)}^3}].$$(25)

Note that the directivity correction of the Gaussian source has been derived for a source moving in uniform rectilinear motion in an inhomogeneous and quiescent atmosphere and is also assumed to be valid for a source in arbitrary motion in a homogeneous and moving atmosphere.

For all simulations, the mesh is uniform with a spatial step defined as Δx = Δy = Δz = 0.1 m. The time step Δt is determined from the Courant-Friedrichs-Lewy (CFL) number, so that CFL = c₀ Δt/Δx = 0.5. The width of the Gaussian source is set to $B=2\Delta x/\sqrt{\log 2}$.

The size of the computational domain varies depending on the test case. For the most computationally demanding configuration, detailed in Section 3.2.4, the computational domain is [−10 m; 210 m] × [−10 m; 18 m] × [0; 15 m]. The number of spatial points in the 3D domain is 2201 × 281 × 151, corresponding to about 93 million points. The simulation time is 4 s; with Δt = 1.47 × 10⁻⁴ s, this leads to 27000 time iterations. The code is parallelized using OpenMP. The simulations have been carried out using 16 core nodes of Intel 6142 Skylake. The CPU time is 1350 h, which corresponds to 84 h of elapsed time.

3.2 Test cases

3.2.1 Influence of source motion effects

The first case study considers a monopole in rectilinear motion at a constant Mach number M_s = 0.3 along the x axis and at a constant height z_s = 2 m above the ground. The ground admittance is specified using the slit-pore model [34] for a porosity of 0.38 and an effective flow resistivity of 514 kPa s m⁻².

Figure 4 compares the instantaneous sound pressure level (SPL), defined by:

$$\mathrm{SPL}(\mathbf{x},t)=20{\log }_{10}\left(\frac{\left|p\right(\mathbf{x},t\left)\right|}{p_{\mathrm{ref}}}\right)$$(26)

with p_ref = 2 × 10⁻⁵ Pa, for several formulations to highlight the importance of the source motion effects. The first one is the FDTD reference solution. The second one is the quasi-static formulation, denoted QS. It neglects the convective amplification and it assumes that the source position at the emission time is that at the reception time. As a consequence, the Doppler effect is not accurately accounted for. The third formulation, noted $\tilde{H}$, also neglects convective amplification, but accounts for the source motion between the emission and reception time. Finally, the fourth formulation is the heuristic formulation that considers both convective amplification and the correct source position at the emission time (see Eq. (25)). It can be observed that the heuristic solution fits perfectly with the reference FDTD solution. Moreover, the QS and $\tilde{H}$ formulations both show large error (up to several dB). This level difference is due to convective amplification that is not taken into account. Indeed, convective amplification tends to increase the SPL as the source approaches the receiver, and reduce the SPL as the source moves away. Furthermore, shifts in interference patterns are observed with the QS formulation. They are caused by phase errors because of the inaccurate position of the source at the emission time in the QS formulation.

Figure 4 Sound pressure level as a function of the time obtained for several formulations for a source in rectilinear uniform motion and for a source frequency of f0 = 500 Hz, considering two receivers located at (40, 0, 0) m (left) and (40, 0, 5) m (right).

We have carried out additional simulations for source frequencies of 100 Hz and 200 Hz and we have also compared the numerical and analytical solutions for a receiver height of 2 m. We found excellent agreement between the two solutions for all frequencies and receiver positions. The comparisons are not presented here for conciseness.

3.2.2 Validation for a 2D trajectory

The second case study is carried out for a monopole moving with a constant acceleration at a constant height z_s = 2 m above the ground. The monopole moves in the (x, y) plane following the equation of motion:

$${\mathbf{x}}_{\mathbf{s}}\left(t\right)={\mathbf{x}}_{\mathbf{s}}\left(t=0\right)+{\mathbf{M}}_s\left(t=0\right)c_{0}t+\frac{1}{2}{\dot{\mathbf{M}}}_s{c}_0{t}^2,$$(27)

where the initial position of the source is x_s(t = 0) = (0, 0, 2 m), the initial Mach number vector is M_s(t = 0) = (0.1, 0, 0), and ${\dot{\mathbf{M}}}_s=(0, 0.2 {\mathrm{s}}^{-1},0)$. Figure 5 presents the source Mach number as a function of the time and the corresponding trajectory.

Figure 5 Mach number (left), and trajectory (right) of the monopole. The receivers are placed at the positions represented by the colored dots.

Figure 6 presents the time evolution of the sound pressure levels at the receivers, for source frequencies of 100 Hz and 500 Hz. A perfect agreement is also found between the FDTD and heuristic solutions at the four receivers and for the two source frequencies. Note that the dips observed at 500 Hz are due to interference between direct and reflected fields.

Figure 6 Comparison of sound pressure level obtained with the FDTD solution and the heuristic solution at m1 = (40, 0, 0) m, m2 = (40, 25, 0) m, m3 = (40, 0, 2) m and m4 = (40, 25, 2) m, for both f0 = 100 Hz (left) and f0 = 500 Hz (right).

3.2.3 Validation for a 3D trajectory

The third test case corresponds to a scenario in which the monopole moves along a 3D helical trajectory. This test case is chosen in order to test the validity of the heuristic formulation for a context beyond that of aircraft noise. The position of the monopole is governed by the following system:

$$x_s=r_s\left(t\right)\cos (2\pi t/T),$$(28)

$$y_s=r_s\left(t\right)\sin (2\pi t/T),$$(29)

$$z_s=v_{z0}t,$$(30)

where v_z0 = 1 m/s, T = 1 s, r_s(t) = a_rt² is the instantaneous radius, and a_r = 0.4 m/s². With such parameters, the Mach number is varying with time (see Fig. 7).

Figure 7 Mach number (left), and trajectory (right) of the monopole for the 3D helical test case. The two receivers m1 = (0, 10, 0) m and m2 = (0, 10, 5) m are represented by the colored dots.

Figure 8 shows the time evolution of the SPL at two receivers located at (0, 10, 0) m and (0, 10, 5) m, for source frequencies of 100 Hz and 500 Hz. The FDTD and heuristic solutions match closely, whether for the receiver at the ground or for the receiver placed at z = 5 m. In the latter case, slight offsets are visible in the interference patterns dips for t > 4 s. Note that the shifts are on the order of hundredths of a second which validates the robustness of the heuristic approach, especially given the spatial scale considered and the fast variation of the Mach number in time of this case study.

Figure 8 Comparison of sound pressure level obtained with the FDTD solution (plain lines) and the heuristic solution (dashed lines) at m1 = (0, 10, 0) m and m2 = (0, 10, 5) m, for both f0 = 100 Hz (left) and f0 = 500 Hz (right).

3.2.4 Validation in presence of flow

This last test case aims at validating the heuristic formulation (15) with refraction effects. A comparison of the FDTD solution and the heuristic solution is presented in the following, considering a monopole moving at constant speed and height in the presence of flow. The wind is directed along the x-direction, i.e. V₀ = (V₀, 0, 0). The flow profile is defined by:

$$V_0\left(z\right)=A_u\log \left(1+\frac{z}{z_0}\right),$$(31)

with A_u = 2 m s⁻¹ and z₀ = 0.5 m. The monopole is moving in a rectilinear motion along the x axis at a constant height z_s = 2 m and Mach number M_s = 0.1. It is initially located at x_s(t = 0) = (0, 0, 2) m.

As an illustration, Figure 9 shows two snapshots of the instantaneous sound pressure levels obtained from the FDTD solution for f₀ = 500 Hz. The effect of acoustic refraction can be observed on the two snapshots. In Figure 9a, focusing of acoustic waves near the ground is noticed for x > 140 m, because of downwind conditions. In Figure 9b, a shadow zone appears near the ground for x < 40 m, because of upwind conditions.

Figure 9 Snapshots of the instantaneous sound pressure levels in the planes y = 0 and z = 2 m obtained from the FDTD solution at time (a) t = 0.706 s and (b) t = 3.971 s. The monopole with a source frequency of f0 = 500 Hz is moving at constant height zs = 2 m and Mach number Ms = 0.1.

Figure 10 presents the eigenrays for some positions of the source obtained with the RT-heuristic formulation (see Eq. (15)) along with the vertical wind profile V₀(z) considered. The source position at different instants in time is represented by the blue points and the receiver position is shown with a red point. The direction of the source motion x_s is indicated by the black arrow. Note that the eigenrays are located in the plane y = 0. The effect of acoustic refraction on the eigenrays can be observed. When the source approaches the receiver (x_s < 80 m), downwind propagation occurs and the rays are bent downwards. Conversely, when the source recedes from the receiver (x_s > 80 m), upwind propagation occurs and the rays are bent upwards. For such propagation conditions, the receiver is in the shadow zone when the source is located at x_s > 140 m.

Figure 10 (left) Eigenrays determined for some positions of the source and for a moving atmosphere. The source positions are indicated by blue dots and the receiver position by a red dot. The corresponding wind vertical profile V0(z) is shown in the right figure.

Figure 11 presents a comparison of the SPL obtained with the FDTD reference solution, with the RT-heuristic approach, and with the heuristic solution in homogeneous atmosphere. Results are shown for f₀ = 200 Hz and f₀ = 500 Hz for a receiver placed at (40, 0, 5) m. The results of the FDTD solution and the RT-heuristic approach match closely, both in terms of amplitude and position of interference patterns. The heuristic solution in homogeneous atmosphere is only valid when the source is close to the receiver (in this case for 2 s < t < 2.6 s), where refraction can be neglected due to the short propagation distance. It is important to note that the RT-heuristic formulation does not allow for estimating SPL in the shadow zone since no rays would reach the receiver under these conditions.

Figure 11 Comparison of the FDTD solution with the RT-heuristic solution and the heuristic solution in a homogeneous atmosphere for the sound pressure level at (80, 0, 5) m above a bare ground and in presence of flow. The pressure field is simulated for f0 = 200 Hz (left) and f0 = 500 Hz (right).

4 Aircraft noise application

This section presents simulations carried out in the context of aircraft noise certification in order to highlight the influence of outdoor propagation effects on the results, using the formulations detailed in Section 2.3. The source follows a climbing trajectory at constant speed, which is described by equation (27) with x_s(t = 0) = (0, 0, 250) m, M_s(t = 0) = (0.295, 0, 0.05) and ${\dot{\mathbf{M}}}_{\mathrm{s}}=(0, 0, 0)$. A receiver that represents the certification measurement point is located at (400, 0, 1.2) m, for which the acoustic pressure is calculated with equation (15). Figure 12 shows the corresponding geometry.

Figure 12 Trajectory of the monopole for the aircraft noise certification scenario. The receiver (red point) is located at (400, 0, 1.2) m and represents the aircraft certification measurement point.

The emission spectrum is composed of broadband pink noise (−3 dB by octave band) along with 3 single tones emerging by 12 dB respectively at f₀ = 63, 250, and 1500 Hz [2]. A first reference simulation is performed assuming a homogeneous atmosphere and a rigid ground. A second simulation is computed for a homogeneous atmosphere and an absorbing grassy ground in order to quantify sound level variability due to a change in the ground admittance. The ground admittance is specified using the slit-pore model [34] for a porosity of 0.76 and an effective flow resistivity of 71.7 kPa s m⁻². A third simulation is performed for a moving atmosphere with a log wind profile with A_u = 2 m s⁻¹ (see Eq. (31)), over the same grassy ground.

For this example, atmospheric absorption is taken into account. For that, the amplitude of the direct and reflected arrivals are multiplied by exp[−α(ω_e)s], with α the absorption coefficient evaluated at the Doppler frequency and s the distance travelled by the acoustic wave from the source at the emission time to the receiver along the direct or reflected ray. The absorption coefficient α is calculated according to ISO-9613 standard [39] for a temperature of 15 °C and a relative humidity of 50%. Finally, the distance s corresponds to the arc length of the eigenray given by:

$$s={\int }_0^{\tau }\mathrm{}\left|\frac{\mathrm{d}{\mathbf{x}}_r}{\mathrm{d}{t}_r}\right|\mathrm{d}{t}_r$$(32)

where the integrand has to be evaluated from the right hand side of equation (18) (see Sect. 2.3).

The spectrograms for the three cases are presented in Figure 13. They are displayed with respect to the frequency at the observer. As the RT-heuristic formulation in equation (15) gives the acoustic pressure as a function of the source frequency, we use equation (20) to express the acoustic pressure as a function of the observer frequency. Note that the relation between the frequency at the observer and the source frequency depends on the propagation path, and thus differs for the direct and the reflected waves. To facilitate the comparison between spectrograms, fixed white dashed lines have been added, and the sound levels are normalized by the maximum of pressure at the receiver of the first simulation. The results highlight a comb-filtering effect that is induced by the interference between the direct wave and the reflected wave. This comb-filtering effect is time dependent: indeed, as the source is in motion, the source-receiver geometry changes over time. The interference patterns are more pronounced in the case of rigid ground, as it is the case for a nonmoving source. It is observed that the levels are higher when the source approaches the receiver (t < 4 s) than when the source moves away from it (t > 4 s) which is due to convective amplification. The positions of the interference patterns are shifted between the reference simulation (a) and the simulations (b) and (c). Regarding simulation (b), the slight shifts are due to a phase difference induced by the change in ground admittance properties. Regarding simulation (c), the interference patterns shifts are also due to ground effects, as well as to refraction by the wind vertical gradient that induces a change in the travel time and Doppler effect. Note that the interference shifts due to the wind are mainly observed when the source is far from the receiver (t < 2 s and t > 6 s). Indeed, for a stratified atmosphere, refraction effect is negligible for rays propagating vertically and tends to be stronger as rays propagate horizontally. Finally, as expected, the effect of atmospheric absorption is mainly visible at high frequencies and large distances, for which a significant reduction of SPL is observed.

Figure 13 Spectrograms of the acoustic pressure at the receiver position (400, 0, 1.2) m for a homogeneous atmosphere and a rigid ground (a), a homogeneous atmosphere and a grassy ground (b), and a moving atmosphere with a log wind profile and a grassy ground (c). Fixed white dashed line are added in order to facilitate the comparison. The levels are normalized by the maximum of pressure at the receiver for the first simulation (a).

This SPL variability can exert a significant impact on the results of aircraft noise acoustic certification because of the change in the overall level, and because the emitted tonal components may alternate between a maximum or a minimum of pressure depending on the outdoor propagation conditions.

5 Conclusion and future work

A coupled approach based on a ray-tracing model and a heuristic formulation has been introduced for predicting the sound radiation by a monopole in arbitrary motion in an inhomogeneous and moving atmosphere. It provides an efficient formulation that includes source motion effects (Doppler effect, convective amplification), as well as outdoor propagation effects (ground reflection, refraction due to temperature and wind speed gradients). This heuristic approach has been validated against a reference FDTD solution for several case studies, with different meteorological conditions and source trajectories. Results showed a perfect agreement with the FDTD solution. Finally, an application focusing on aircraft noise certification conditions was proposed. An acoustic source emitting a broadband noise and three single tones was considered. The study highlighted the comb-filtering effect resulting from interference between direct and reflected rays, which is time-dependent due to the motion of the source. Shifts in interference patterns leading to variability of SPL were observed and were attributed to outdoor propagation effects such as wind gradient and ground absorption. Convective amplification was also noticed, leading to a SPL higher when the source approaches the receiver than when it moves away from it.

The variability of SPL can significantly influence the measurements for aircraft noise certification, particularly if tonal components alternate between interference maxima and minima. Thus, the heuristic formulation presented in this paper can be employed to guide the design of new aircraft noise control strategies. Indeed, it enables modeling accurately outdoor sound propagation for a moving source, while maintaining a low computational cost compared to reference wave-based approaches. However, the formulation still has some improvement points that could be addressed in subsequent works. Regarding noise propagation, considering atmospheric turbulence is needed to improve the accuracy of the results in the interference patterns region as well as in shadow zones. Phase decorrelation methods which leads to a coherence loss in ground effects can be employed to simulate atmospheric turbulence effects in the interference patterns region. As no ray propagates in shadow zones, engineering solutions, such as the one proposed in the Harmonoise project [40], can provide a first-order approximation for modeling the decrease in SPL. Regarding noise emission, taking into account several sources placed on the aircraft, each of them with a specific directivity and spectrum, would allow for a better consideration of aircraft noise emission characteristics.

Acknowledgments

The present work is part of the program MAMBO “Méthodes avancées pour la modélisation du bruit moteur et avion” (Advanced methods for engine and aicraft noise modelling” coordinated by Airbus SAS. It was supported by the Direction Générale de l’Aviation Civile (DGAC) under the Grant n° 2021-50. This publication was performed within the framework of the LABEX CeLyA (ANR-10-LABX-0060) of Universite de Lyon, within the program “Investissements d’Avenir” (ANR-16-IDEX-0005) operated by the French National Research Agency (ANR). This work was granted access to the HPC resources of PMCS2I (Pole de Modelisation et de Calcul en Sciences de l’Ingenieur et de l’Information) of Ecole Centrale de Lyon and PSMN (Pole Scientifique de Modelisation Numerique) of ENS de Lyon, members of FLMSN (Federation Lyonnaise de Modelisation et Sciences Numeriques), partner of EQUIPEX EQUIP@MESO.

Conflicts of interest

The authors declare that they have no conflicts of interest in relation to this article.

Data availability statement

The data are available from the corresponding author on request.

Appendix A

Ray equations

In this appendix, we detail the implementation of the ray approach. In addition to the six ray equations indicated in equations (18) and (19), 12 equations are solved for x_φ, x_θ, n_φ and n_θ in order to determine the amplitude along the rays. The equations for x_φ and n_φ write:

$$\frac{\mathrm{d}{\mathbf{x}}_{\phi }}{\mathrm{d}{t}_r}=({\mathbf{x}}_{\phi }·\nabla c_0)\mathbf{n}+c_0{\mathbf{n}}_{\phi }+({\mathbf{x}}_{\phi }·\nabla ){\mathbf{V}}_0,$$(A.1)

$$\frac{\mathrm{d}{\mathbf{n}}_{\phi }}{\mathrm{d}{t}_r}=({\mathbf{q}}_{\phi }·\mathbf{n})\mathbf{n}+(\mathbf{q}·{\mathbf{n}}_{\phi })\mathbf{n}+(\mathbf{q}·\mathbf{n}){\mathbf{n}}_{\phi }-{\mathbf{q}}_{\phi },$$(A.2)

with q_φ = (x_φ · ∇)∇c₀ + [(x_φ · ∇)∇ V₀] · n + (∇V₀) · n_φ. The equations for x_θ and n_θ are similar to equations (A.1) and (A.2) with the subscript φ replaced by θ.

Using Einstein notation, the full system of ray equations writes:

$$\frac{\mathrm{d}{x}_{\mathrm{ri}}}{\mathrm{d}{t}_r}=c_0{n}_i+V_{0i}$$(A.3)

$$\frac{\mathrm{d}{n}_i}{\mathrm{d}{t}_{\mathrm{r}}}=q_j{n}_j{n}_i-q_i,$$(A.4)

$$\frac{\mathrm{d}{x}_{\phi i}}{\mathrm{d}{t}_{\mathrm{r}}}=x_{\phi j}\frac{\partial c_0}{\partial x_j}{n}_i+c_0{n}_{\phi i}+x_{\phi j}\frac{\partial V_{0i}}{\partial x_j},$$(A.5)

$$\frac{\mathrm{d}{n}_{\phi i}}{\mathrm{d}{t}_{\mathrm{r}}}=q_{\phi j}{n}_j{n}_i+q_j{n}_{\phi j}{n}_i+q_j{n}_j{n}_{\phi i}-q_{\phi i},$$(A.6)

with the following expressions for the components of vectors q and q_φ:

$$q_i=\frac{\partial c_0}{\partial x_i}+n_j\frac{\partial V_{0j}}{\partial x_i},$$(A.7)

$$q_{\phi i}=x_{\phi j}\frac{{\partial }^2{c}_0}{\partial x_i\partial x_j}+n_j{x}_{\phi k}\frac{{\partial }^2{V}_{0j}}{\partial x_i\partial x_k}+n_{\phi j}\frac{\partial V_{0j}}{\partial x_i}.$$(A.8)

The ray variables are initialized with the relations:

$${\mathbf{x}}_r(t_r=0)={\mathbf{x}}_s\left(t_e\right),$$(A.9)

$$\mathbf{n}(t_r=0)=[\cos \phi \cos \theta ,\cos \phi \sin \theta ,\sin \phi {]}^T,$$(A.10)

$${\mathbf{n}}_{\phi }(t_r=0)=[-\sin \phi \cos \theta ,-\sin \phi \sin \theta ,\cos \phi {]}^T,$$(A.11)

$${\mathbf{n}}_{\theta }(t_r=0)=[-\cos \phi \sin \theta ,\cos \phi \cos \theta ,0{]}^T,$$(A.12)

$${\mathbf{x}}_{\phi }(t_r=0)=\mathbf{0},$$(A.13)

$${\mathbf{x}}_{\theta }(t_r=0)=\mathbf{0}.$$(A.14)

Finally, an extra attention has to be taken for the reflection of the rays from the ground. In particular, the reflection conditions are of crucial importance for the determination of the amplitude along the ray. We introduce the unit vector normal to the ground oriented towards the air e_n. For a flat horizontal ground, one has e_n = [0, 0, 1]. Following Gainville [41], once the ray impinges on the ground, the ray variables have to be reinitialized. By noting the updated ray variables after reflection with the superscript ′, one has:

$${\mathbf{x}}_r^{\mathrm{\prime }}={\mathbf{x}}_r,$$(A.15)

$${\mathbf{n}}^{\mathbf{\text{'}}}=\mathbf{n}-2(\mathbf{n}·{\mathbf{e}}_{\mathbf{n}}){\mathbf{e}}_{\mathbf{n}},$$(A.16)

$${\mathbf{x}}_{\phi }^{\mathbf{\text{'}}}={\mathbf{x}}_{\phi }-2({\mathbf{x}}_{\phi }·{\mathbf{e}}_{\mathbf{n}}){\mathbf{e}}_{\mathbf{n}},$$(A.17)

$${\mathbf{x}}_{\theta }^{\mathrm{\prime }}={\mathbf{x}}_{\theta }-2({\mathbf{x}}_{\theta }·{\mathbf{e}}_{\mathbf{n}}){\mathbf{e}}_{\mathbf{n}},$$(A.18)

$${\mathbf{n}}_{\phi }^{\text{'}}={\mathbf{n}}_{\phi }-2({\mathbf{n}}_{\phi }·{\mathbf{e}}_{\mathbf{n}}){\mathbf{e}}_{\mathbf{n}}-\frac{{\mathbf{x}}_{\phi }·{\mathbf{e}}_{\mathbf{n}}}{c_0\mathbf{n}·{\mathbf{e}}_{\mathbf{n}}}\left[\frac{\mathrm{d}\mathbf{n}}{\mathrm{d}t}-\frac{\mathrm{d}\mathbf{n}\mathrm{\prime }}{\mathrm{d}t}-2\left(\frac{\mathrm{d}\mathbf{n}}{\mathrm{d}t}·{\mathbf{e}}_{\mathbf{n}}\right){\mathbf{e}}_{\mathbf{n}}\right],$$(A.19)

$${\mathbf{n}}_{\theta }^{\mathrm{\prime }}={\mathbf{n}}_{\theta }-2\left({\mathbf{n}}_{\theta }·{\mathbf{e}}_{\mathbf{n}}\right){\mathbf{e}}_{\mathbf{n}}-\frac{{\mathbf{x}}_{\theta }·{\mathbf{e}}_{\mathbf{n}}}{c_0\mathbf{n}·{\mathbf{e}}_{\mathbf{n}}}\left[\frac{\mathrm{d}\mathbf{n}}{\mathrm{d}t}-\frac{\mathrm{d}{\mathbf{n}}^{\mathrm{\prime }}}{\mathrm{d}t}-2\left(\frac{\mathrm{d}\mathbf{n}}{\mathrm{d}t}·{\mathbf{e}}_{\mathbf{n}}\right){\mathbf{e}}_{\mathbf{n}}\right].$$(A.20)

Note that dn/dt and dn′/dt have to be evaluated from the right-hand side of equation (19), using either the ray variables before or after reflection.

Appendix B

Expression of the energy density flux

In this appendix, we detail the derivation of the energy flux in equation (22) for a monopole moving at uniform speed in an uniformly moving medium.

B.1 Analytical solution for a monopole moving at uniform speed in an uniformly moving medium

We briefly present the derivation of the analytical solution for a monopole moving at a constant Mach number vector M_s in a uniformly moving atmosphere with Mach number vector M₀. First, we use the change of coordinates x′ = x − M₀c₀t and t′ = t. Equations (1) and (2) write:

$$\frac{\partial p}{\partial t\mathrm{\prime }}+{\rho }_0{c}_0^2\nabla \mathrm{\prime }·\mathbf{v}=c_0^2{S}_0 {\mathrm{e}}^{-\mathrm{i}{\omega }_0{t}^{\mathrm{\prime }}}\delta [{\mathbf{x}}^{\mathbf{\text{'}}}-({\mathbf{M}}_{\mathrm{s}}-{\mathbf{M}}_0\left)c_{0}t\mathrm{\prime }\right],$$(B.1)

$$\frac{\partial \mathbf{v}}{\partial \mathrm{t\text{'}}}+\frac{1}{{\rho }_0}\nabla \mathrm{\text{'}p}=0.$$(B.2)

where ∇′ = ∂/∂x′. In the system of coordinates (x′, t′), the problem is equivalent to that of a monopole moving at Mach number M_s − M₀ in an homogeneous atmosphere at rest, for which an analytical solution is available. Back in the original system of coordinates (x, t), the acoustic pressure can be written as [31]:

$$p(\mathbf{x},t)=-\frac{\mathrm{i}{\omega }_0{\mathrm{e}}^{-\mathrm{i}{\omega }_{0}t}}{4\pi }\frac{{\mathrm{e}}^{\mathrm{i}{k}_0{R}_e}}{R_e(1-M_{s0}\cos {\psi }_e{)}^2}$$(B.3)

where R_e and cosψ_e are related to the coordinates at the reception time by the relations:

$$R_e=\frac{R}{{\beta }^2}\left(M_{s0}\cos \psi +\sqrt{M_{s0}^2{\cos }^2\psi +{\beta }^2}\right)$$(B.4)

$$\cos {\psi }_e=M_{s0}+\frac{R}{R_e}\cos \psi$$(B.5)

with R = x − M_sc₀t, R = |R|, cosψ = [R · (M_s − M₀)]/(RM_s0), M_s0 = |M_s − M₀|, and $\beta =\sqrt{1-M_{s0}^2}$.

B.2 Pressure amplitude in the ray-tracing approach

The conservation of energy flux along a raytube writes [30, 42]:

$$\frac{A(t_r{)}^2}{{\rho }_0\left(t_r\right)c_0\left(t_r\right)}[1+{\mathbf{M}}_{\mathbf{0}}(t_r)·\mathbf{n}(t_r\left){]}^{2}S\right(t_r)=\mathrm{cste}=E(\phi ,\theta ).$$(B.6)

We can then deduce the pressure amplitude with:

$$A\left(t_r\right)=\sqrt{\frac{E\left(\phi ,\theta \right){\rho }_0\left(t_r\right)c_0\left(t_r\right)}{S\left(t_r\right)}}\frac{1}{\left[1+{\mathbf{M}}_{\mathbf{0}}\left(t_r\right)·\mathbf{n}\left(t_r\right)\right]}.$$(B.7)

B.3 Ray solution for a monopole moving at uniform speed in an uniformly moving medium

For a homogeneous medium with a uniform velocity, we can integrate the 18 ray equations from t_r = 0 to τ to obtain the ray solution:

$${\mathbf{x}}_r\left(\tau \right)={\mathbf{x}}_S\left(t_e\right)+\left({\mathbf{M}}_0+\mathbf{n}\right)c_0\tau ={\mathbf{x}}_S\left(\tau \right)+\left({\mathbf{M}}_0-{\mathbf{M}}_{\mathrm{s}}+\mathbf{n}\right)c_0\tau ,$$(B.8)

$$\mathbf{n}\left(\tau \right)=\mathbf{n}\left(0\right),$$(B.9)

$${\mathbf{n}}_{\phi }\left(\tau \right)={\mathbf{n}}_{\phi }\left(0\right),$$(B.10)

$${\mathbf{n}}_{\theta }\left(\tau \right)={\mathbf{n}}_{\theta }\left(0\right),$$(B.11)

$${\mathbf{x}}_{\phi }\left(\tau \right)={\mathbf{n}}_{\phi }\left(0\right)c_0\tau ,$$(B.12)

$${\mathbf{x}}_{\theta }\left(\tau \right)={\mathbf{n}}_{\theta }\left(0\right)c_0\tau .$$(B.13)

One can then deduce the raytube cross-section: S(τ) = |x_φ × x_θ| = (c₀τ)²|cosφ|.

The ray solution writes:

$$p(\mathbf{x},t)=-\frac{\mathrm{i}{\omega }_0{\mathrm{e}}^{-\mathrm{i}{\omega }_{0}t}}{4\pi }A\left(\tau \right){\mathrm{e}}^{\mathrm{i}{\omega }_0\tau }.$$(B.14)

By identification with equation (B.3), one deduces that R_e = c₀τ. The ray trajectory in equation (B.8) can be rewritten as R = (M₀ − M_s + n)R_e. By projecting R on M₀ − M_s, one finds with equation (B.5) the equivalence M_s0cosψ_e = n · (M_s − M₀).

Finally, in order that equation (B.14) matches with equation (B.3), one has the relation:

$$A\left(\tau \right)=\frac{1}{c_0\tau [1+\mathbf{n}·({\mathbf{M}}_{\mathbf{0}}-{\mathbf{M}}_{\mathbf{s}}){]}^2}.$$(B.15)

One deduces from equation (B.6), the expression for the energy density flux:

$$E(\phi ,\theta )=\frac{\left|\cos \phi \right|}{{\rho }_0{c}_0}\frac{[1+{\mathbf{M}}_{\mathbf{0}}·\mathbf{n}{]}^2}{[1+\mathbf{n}·({\mathbf{M}}_{\mathbf{0}}-{\mathbf{M}}_{\mathbf{s}}){]}^4}.$$(B.16)

In the case of a moving source in arbitrary motion or/and in an inhomogeneous moving medium, the parameters in equation (B.16) are not constant and are evaluated at the source position and at the emission time, which gives equation (22).

References

All Figures

	Figure 1 Geometry considered for a moving monopole in rectilinear motion above the ground.
In the text

	Figure 2 Geometry considered for a moving monopole in arbitrary motion above the ground.
In the text

	Figure 3 Geometry considered for a moving monopole in arbitrary motion in a inhomogeneous atmosphere above a ground.
In the text

	Figure 4 Sound pressure level as a function of the time obtained for several formulations for a source in rectilinear uniform motion and for a source frequency of f0 = 500 Hz, considering two receivers located at (40, 0, 0) m (left) and (40, 0, 5) m (right).
In the text

	Figure 5 Mach number (left), and trajectory (right) of the monopole. The receivers are placed at the positions represented by the colored dots.
In the text

	Figure 6 Comparison of sound pressure level obtained with the FDTD solution and the heuristic solution at m1 = (40, 0, 0) m, m2 = (40, 25, 0) m, m3 = (40, 0, 2) m and m4 = (40, 25, 2) m, for both f0 = 100 Hz (left) and f0 = 500 Hz (right).
In the text

	Figure 7 Mach number (left), and trajectory (right) of the monopole for the 3D helical test case. The two receivers m1 = (0, 10, 0) m and m2 = (0, 10, 5) m are represented by the colored dots.
In the text

	Figure 8 Comparison of sound pressure level obtained with the FDTD solution (plain lines) and the heuristic solution (dashed lines) at m1 = (0, 10, 0) m and m2 = (0, 10, 5) m, for both f0 = 100 Hz (left) and f0 = 500 Hz (right).
In the text

	Figure 9 Snapshots of the instantaneous sound pressure levels in the planes y = 0 and z = 2 m obtained from the FDTD solution at time (a) t = 0.706 s and (b) t = 3.971 s. The monopole with a source frequency of f0 = 500 Hz is moving at constant height zs = 2 m and Mach number Ms = 0.1.
In the text

	Figure 10 (left) Eigenrays determined for some positions of the source and for a moving atmosphere. The source positions are indicated by blue dots and the receiver position by a red dot. The corresponding wind vertical profile V0(z) is shown in the right figure.
In the text

	Figure 11 Comparison of the FDTD solution with the RT-heuristic solution and the heuristic solution in a homogeneous atmosphere for the sound pressure level at (80, 0, 5) m above a bare ground and in presence of flow. The pressure field is simulated for f0 = 200 Hz (left) and f0 = 500 Hz (right).
In the text

	Figure 12 Trajectory of the monopole for the aircraft noise certification scenario. The receiver (red point) is located at (400, 0, 1.2) m and represents the aircraft certification measurement point.
In the text

Figure 13 Spectrograms of the acoustic pressure at the receiver position (400, 0, 1.2) m for a homogeneous atmosphere and a rigid ground (a), a homogeneous atmosphere and a grassy ground (b), and a moving atmosphere with a log wind profile and a grassy ground (c). Fixed white dashed line are added in order to facilitate the comparison. The levels are normalized by the maximum of pressure at the receiver for the first simulation (a).

In the text