SCIENTIFIC ARTICLE

OPEN ACCESS

# Heuristic solution for the acoustic radiation of a moving monopole in an inhomogeneous and moving atmosphere. Application to aircraft noise

Bill Kayser*, Didier Dragna, and Philippe Blanc-Benon

Ecole Centrale de Lyon, CNRS, Université Claude Bernard Lyon 1, INSA Lyon, LMFA, UMR5509, 36 avenue Guy de Collongue, 69130 Ecully, France

Received 19 April 2024, Accepted 5 August 2024

**Abstract** – A coupled approach is proposed for predicting sound radiation from a monopole in arbitrary motion in a moving and inhomogeneous atmosphere. It is based on a heuristic model proposed in the literature for sound radiation by a moving source in a homogeneous atmosphere at rest above an absorbing ground and a ray-tracing code, which takes into account meteorological effects. Validation of the model is performed with a reference three-dimensional finite-difference time-domain solution of the linearized Euler equations for several test-cases with different source trajectories and atmospheric conditions. We show that neglecting convective amplification or the source motion between the emission and reception times can lead to significant errors in the predictions. Finally, an application case for aircraft noise is presented. The importance of ground and meteorological effects on the sound pressure levels is highlighted.

**Keywords:** Aircraft noise, Ray-tracing, Heuristic formulation, Outdoor propagation

## 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

<footnote>*Corresponding author: bill.kayser@ec-lyon.fr</footnote>
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 $\mathbf{x} = (x, y, z)$, the source trajectory is denoted $\mathbf{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}_0 \cdot \nabla p + \rho_0 c_0^2 \nabla \cdot \mathbf{v} = c_0^2 S_0 e^{-i\omega_0 t} \delta[\mathbf{x} - \mathbf{x}_s(t)], \quad (1) $$

$$ \frac{\partial \mathbf{v}}{\partial t} + \mathbf{V}_0 \cdot \nabla \mathbf{v} + \frac{1}{\rho_0} \nabla p + \mathbf{v} \cdot \nabla \mathbf{V}_0 = 0. \quad (2) $$

where $p$ and $\mathbf{v}$ are respectively the acoustic pressure and velocity, $\rho_0$ is the air density, $c_0$ is the sound speed, $\mathbf{V}_0$ is the base flow velocity, and $\omega_0 = 2\pi f_0$ with $f_0$ the source frequency. The source amplitude $S_0$ 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 + [\tilde{\beta} * p] = 0, \quad (3) $$

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

From the source trajectory, we define the source speed vector $\mathbf{V}_s = d\mathbf{x}_s/dt$, the Mach number vector $\mathbf{M}_s = \mathbf{V}_s/c_0$, the Mach number $M_s = |\mathbf{M}_s|$, and the time derivative of the Mach number vector $\dot{\mathbf{M}}_s = d\mathbf{M}_s/dt$. We introduce also the Mach number vector associated to the mean flow as $\mathbf{M}_0 = \mathbf{V}_0/c_0$. Finally, we define the source-receiver distance at the time $t$ as $R(t) = |\mathbf{R}(t)|$ with $\mathbf{R}(t) = \mathbf{x} - \mathbf{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 $\rho_0 = 1.2 \text{ kg m}^{-3}$ and $c_0 = 340 \text{ m s}^{-1}$. 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_0$ is constant in this section and $\mathbf{V}_0 = 0$. Without loss of generality, the source position is chosen as $\mathbf{x}_s(t) = (M_s c_0 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 $\varphi$ and $\theta$ 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)).

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{i\omega_0 e^{-i\omega_0 t}}{4\pi} \left[ \frac{e^{ik_0 R_{e1}}}{R_{e1}(1 - M_{r1})^2} + \frac{e^{ik_0 R_{e2}}}{R_{e2}(1 - M_{r2})^2} \right], \quad (4) $$

where $k_0 = \omega_0/c_0$ is the wave number and

$$ M_{r1} = \mathbf{M}_s \cdot \mathbf{R}_{e1} / R_{e1} = M_s \cos \varphi_1 \cos \theta_1, \quad (5) $$

$$ M_{r2} = \mathbf{M}_s \cdot \mathbf{R}_{e2} / R_{e2} = M_s \cos \varphi_2 \cos \theta_2, \quad (6) $$

are the projections of the Mach number vector in the source-receiver directions. Even if the source signal is
<figure><img src="image_1.png" /><figcaption>Figure 1. Geometry considered for a moving monopole in rectilinear motion above the ground.</figcaption></figure>

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_{e1} = \omega_0/(1 - M_{r1})$. Similarly, the frequency associated to the reflected wave is $\omega_{e2} = \omega_0/(1 - M_{r2})$. It differs from that associated to the direct wave $\omega_{e1}$, 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:

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

where

$$F(u) = 1 + iu\sqrt{\pi w}(u) \qquad (8)$$

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

$$\mathcal{R} = \frac{\sin \varphi_2 - \beta(\omega_{e2})}{\sin \varphi_2 + \beta(\omega_{e2})}, \qquad (9)$$

is the plane wave reflection coefficient, and

$$d = \frac{1}{2}(1 + i)\sqrt{\frac{k_0 R_{e2}}{1 - M_{r2}}}[\sin \varphi_2 + \beta(\omega_{e2})], \qquad (10)$$

is the numerical distance. Note that the ground surface admittance is evaluated at the Doppler frequency associated to the reflected wave $\omega_{e2}$ and not at the source frequency. In Reference [18], it was shown that evaluating the ground admittance at the source frequency $\omega_0$ instead of the Doppler frequency $\omega_{e2}$ 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)^{-2}$ 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, \qquad (11)$$

where $R_e = |\mathbf{R}_e| = R(t_e)$ is the source-receiver distance at the emission time and $\mathbf{R}_e = \mathbf{x} - \mathbf{x}_s(t_e)$. One has:

$$R_e = \sqrt{(x - M_s c_0 t_e)^2 + y^2 + (z \pm z_s)^2}. \qquad (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)}. \qquad (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
<figure><img src="image_2.png" /><figcaption>Figure 2. Geometry considered for a moving monopole in arbitrary motion above the ground.</figcaption></figure>

motion in free-field (see, e.g. [32]), but not above an absorbing ground.

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

$$p = -\frac{i\omega_0 e^{-i\omega_0 t}}{4\pi} \left[ \left(1 - M_{r1} + \frac{i\dot{M}_{r1}}{\omega_0}\right) \frac{e^{ik_0 R_{e1}}}{R_{e1}(1 - M_{r1})^3} + Q \left(1 - M_{r2} + \frac{i\dot{M}_{r2}}{\omega_0}\right) \frac{e^{ik_0 R_{e2}}}{R_{e2}(1 - M_{r2})^3} \right]. \quad (14)$$

In this formulation, $M_r = M_s \cdot 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{M}_s \cdot R_e / R_e$ is the projection of the vector $M_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 = 0$. From the emission time, the geometric quantities in equation (14) can be determined.

It is noteworthy that $\dot{M}_r$ becomes negligible when it is small compared to $\omega_0(1 - M_r)$. Denoting by $\tau_a$ a characteristic time of the source acceleration, this occurs if $M \ll \omega_0\tau_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.

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{i\omega_0 e^{-i\omega_0 t}}{4\pi} [A_1 e^{i\omega_0 \tau_1} + QA_2 e^{i\omega_0 \tau_2}]. \quad (15)$$

with $\tau$ the travel time from the source at the emission time to the receiver, i.e. $\tau = t - t_e$, and $A$ the amplitude, that
<figure><img src="image_3.png" /><figcaption>**Figure 3.** Geometry considered for a moving monopole in arbitrary motion in a inhomogeneous atmosphere above a ground.</figcaption></figure>

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 $Q$ are expressed as:

$$ \mathcal{R} = \frac{\sin \varphi_g - \beta(\omega_{e2})}{\sin \varphi_g + \beta(\omega_{e2})}, \quad (16) $$

$$ d = \frac{1}{2}(1 + i)\sqrt{\omega_{e2}\tau_2[\sin \varphi_g + \beta(\omega_{e2})]}, \quad (17) $$

with $\omega_{e2}$ the Doppler frequency associated to the reflected wave (see Eq. (20)), and $\varphi_g$ the angle of incidence of the ray on the ground. It should be remarked that $\varphi_g$ differs from $\varphi_2$ 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 $\mathbf{x}_r$ the ray position and by $\mathbf{n}$ the local unit vector normal to the wavefront, the six ray equations are written as (see, e.g. [11, 12, 33]):

$$ \frac{d\mathbf{x}_r}{dt_r} = c_0\mathbf{n} + \mathbf{V}_0 \quad (18) $$

$$ \frac{dn}{dt_r} = (\mathbf{q} \cdot \mathbf{n})\mathbf{n} - \mathbf{q}, \quad (19) $$

with $t_r$ the travel time along the ray and $\mathbf{q} = \nabla c_0 + (\nabla \mathbf{V}_0) \cdot \mathbf{n}$, where $(\nabla \mathbf{V}_0) \cdot \mathbf{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. $\mathbf{x}_r(t_r = 0) = \mathbf{x}_s(t_e)$ and at $t_r = \tau$, the eigenray reaches the receiver, i.e. $\mathbf{x}_r(t_r = \tau) = \mathbf{x}$. The ray equations are solved using the standard fourth order Runge-Kutta algorithm.

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

$$ \omega_e = \omega_0 \frac{1 + \mathbf{M}_0(0) \cdot \mathbf{n}(0)}{1 + [\mathbf{M}_0(0) - \mathbf{M}_s(0)] \cdot \mathbf{n}(0)}, \quad (20) $$

where $\mathbf{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, $\mathbf{M}_0(0)$ is the Mach number vector associated to the flow at the source position at the emission time and $\mathbf{M}_s(0)$ is the source Mach number vector at the emission time.

We have also $\sin\varphi_g = -\mathbf{n}_g \cdot \mathbf{e}_n$, with $\mathbf{n}_g$ is the unit vector normal to the incident wavefront on the ground and $\mathbf{e}_n$ the unit vector normal to the ground surface directed towards the air (for a flat ground, $\mathbf{e}_n = [0, 0, 1]$).

As we consider a monopole, the rays are parameterized using two angles $\varphi$ and $\theta$, 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 $\theta$ and $\varphi$; denoting $\mathbf{x}_\varphi = d\mathbf{x}_r/d\varphi$ and $\mathbf{x}_\theta = d\mathbf{x}_r/d\theta$, one has $S(t_r) = |\mathbf{x}_\varphi \times \mathbf{x}_\theta|$. The evolution of $\mathbf{x}_\varphi$ and $\mathbf{x}_\theta$ 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(\tau) = \sqrt{\frac{\rho_0(\tau)c_0(\tau)E(\varphi, \theta)}{S(\tau)}} \frac{1}{1 + \mathbf{M}_0(\tau) \cdot \mathbf{n}(\tau)} \quad (21) $$

with $E(\varphi, \theta)$ 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(\varphi, \theta) = \frac{|\cos \varphi|}{\rho_0(0)c_0(0)} \frac{[1 + \mathbf{M}_0(0) \cdot \mathbf{n}(0)]^2}{[1 + (\mathbf{M}_0(0) - \mathbf{M}_s(0)) \cdot \mathbf{n}(0)]^4} \quad (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 $\mathbf{V}_0$ by $-\mathbf{V}_0$ in equations (A.18) and (A.19). A coarse discretization of the ray parameters $\theta$ and $\varphi$ 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}(\theta, \varphi)$ for each source position. We then determine the local minima of $d_{\min}(\theta, \varphi)$. 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 $(\theta, \varphi)$ 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(\mathbf{x}) = \frac{1}{\pi^{3/2} B^3} \exp \left( -\frac{|\mathbf{x}|^2}{B^2} \right), \quad (23)$$

with $B$ the Gaussian width. Note that $S(\mathbf{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_0$, by:

$$p_B(\mathbf{x}, t) = \exp \left( -\frac{\omega_0^2 B^2}{4c_0^2} \right) p_0(\mathbf{x}, t). \quad (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{i\omega_0 e^{-i\omega_0 t}}{4\pi} \left[ \exp\left(-\frac{k_0^2 B^2}{4(1-M_{r_1})^2}\right) \left(1-M_{r_1} + \frac{i\dot{M}_{r_1}}{\omega_0}\right) \right. \\ \times \frac{e^{ik_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. \times \left(1-M_{r_2} + \frac{i\dot{M}_{r_2}}{\omega_0}\right) \frac{e^{ik_0 R_{e2}}}{R_{e2}(1-M_{r_2})^3} \right]. \quad (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 $\Delta x = \Delta y = \Delta z = 0.1$ m. The time step $\Delta t$ is determined from the Courant-Friedrichs-Lewy (CFL) number, so that $\text{CFL} = c_0 \Delta t / \Delta x = 0.5$. The width of the Gaussian source is set to $B = 2\Delta x / \sqrt{\log 2}$.
<figure><img src="image_4.png" /><figcaption>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 $f_0 = 500$ Hz, considering two receivers located at (40, 0, 0) m (left) and (40, 0, 5) m (right).</figcaption></figure>

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 \text{ m}; 210 \text{ m}] \times [-10 \text{ m}; 18 \text{ m}] \times [0; 15 \text{ m}]$. The number of spatial points in the 3D domain is $2201 \times 281 \times 151$, corresponding to about 93 million points. The simulation time is 4 s; with $\Delta t = 1.47 \times 10^{-4}$ 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:

$$ \text{SPL}(\mathbf{x}, t) = 20\log_{10}\left(\frac{|p(\mathbf{x}, t)|}{P_{\text{ref}}}\right) \quad (26) $$

with $p_{\text{ref}} = 2 \times 10^{-5}$ 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.

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}_s(t) = \mathbf{x}_s(t=0) + \mathbf{M}_s(t=0)c_0t + \frac{1}{2}\dot{\mathbf{M}}_s c_0 t^2, \quad (27) $$

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

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><img src="image_5.png" /><figcaption>**Figure 5.** Mach number (left), and trajectory (right) of the monopole. The receivers are placed at the positions represented by the colored dots.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 6.** Comparison of sound pressure level obtained with the FDTD solution and the heuristic solution at $m_1 = (40, 0, 0)$ m, $m_2 = (40, 25, 0)$ m, $m_3 = (40, 0, 2)$ m and $m_4 = (40, 25, 2)$ m, for both $f_0 = 100$ Hz (left) and $f_0 = 500$ Hz (right).</figcaption></figure>

### 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(t) \cos(2\pi t/T), \quad (28)$$

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

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

where $v_{z0} = 1$ m/s, $T = 1$ s, $r_s(t) = a_r t^2$ 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 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.

### 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_0 = (V_0, 0, 0)$. The flow profile is defined by:
<figure><img src="image_7.png" /><figcaption>**Figure 7.** Mach number (left), and trajectory (right) of the monopole for the 3D helical test case. The two receivers $m_1 = (0, 10, 0)$ m and $m_2 = (0, 10, 5)$ m are represented by the colored dots.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 8.** Comparison of sound pressure level obtained with the FDTD solution (plain lines) and the heuristic solution (dashed lines) at $m_1 = (0, 10, 0)$ m and $m_2 = (0, 10, 5)$ m, for both $f_0 = 100$ Hz (left) and $f_0 = 500$ Hz (right).</figcaption></figure>

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

with $A_u = 2$ m s⁻¹ and $z_0 = 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 $\mathbf{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_0 = 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 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_0(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

$\mathbf{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 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_0 = 200$ Hz and $f_0 = 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
<figure><img src="image_9.png" /><figcaption>**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 $f_0 = 500$ Hz is moving at constant height $z_s = 2$ m and Mach number $M_s = 0.1$.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**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 $V_0(z)$ is shown in the right figure.</figcaption></figure>

does not allow for estimating SPL in the shadow zone since no rays would reach the receiver under these conditions.

## 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 $\mathbf{x}_s(t = 0) = (0, 0, 250)$ m, $\mathbf{M}_s(t = 0) = (0.295, 0, 0.05)$ and $\mathbf{M}_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.

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_0 = 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 \text{ kPa s m}^{-2}$. A third simulation is
<figure><img src="image_11.png" /><figcaption>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 $f_0 = 200$ Hz (left) and $f_0 = 500$ Hz (right).</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>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.</figcaption></figure>

performed for a moving atmosphere with a log wind profile with $A_u = 2 \text{ m s}^{-1}$ (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[-\alpha(\omega_c)s]$, with $\alpha$ 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 $\alpha$ 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} \left| \frac{dx_r}{dt_r} \right| dt_r \qquad (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.

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.
<figure><img src="image_13.png" /><figcaption>**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).</figcaption></figure>

# 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 aircraft 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 Université 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 Modélisation et de Calcul en Sciences de l'Ingenieur et de l'Information) of Ecole Centrale de Lyon and PSMN (Pole Scientifique de Modélisation Numerique) of ENS de Lyon, members of FLMSN (Fédération Lyonnaise de Modélisation 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.

# References

1. M. Albert, P. Bousquet, D. Lizarazu; Ground effects for aircraft noise certification (2017-3845), in: 23rd AIAA/CEAS Aeroacoustics Conference, Denver, CO, USA, 5-9 June, American Institute of Aeronautics and Astronautics, 2017, AIAA paper 2017-3845.

2. P. Bousquet, V.P. Blandeau, Feasibility of determining aircraft certification noise levels using ground plane micro-
phone measurements, in: AIAA AVIATION 2021 Forum, August 2–6, Virtual event, American Institute of Aeronautics and Astronautics, 2021, AIAA paper 2021-2159.

3. E. Nesbitt, J. Lan, S. Hunkler: Microphone acoustic characteristics for aircraft flyover testing, in: AIAA Aviation 2020 Forum, June 15–19, Virtual event, American Institute of Aeronautics and Astronautics, 2020, AIAA paper 2020-2613.

4. J.M. Giannakis: Evaluation of a correction factor for flyover-noise ground plane microphones, in: AIAA Aviation 2020 Forum, June 15–19, Virtual event, American Institute of Aeronautics and Astronautics, 2020, AIAA paper 2020-2612.

5. V.P. Blandeau, P. Bousquet: A new plate design to improve the accuracy of aircraft exterior noise measurements on the ground, in: AIAA AVIATION 2021 Forum, August 2–6, Virtual event, American Institute of Aeronautics and Astronautics, 2021, AIAA paper 2021-2158.

6. T. Norum, C. Liu: Point source moving above a finite impedance reflecting plane – experiment and theory, Journal of the Acoustical Society of America 63, 4 (1978) 1069–1073.

7. S. Oie, R. Takeuchi: Sound radiation from a point source moving in parallel to a plane surface of porous material, Acta Acustica united with Acustica 48, 3 (1981) 123–129.

8. M. Ochmann: Exact solutions for sound radiation from a moving monopole above an impedance plane, Journal of the Acoustical Society of America 133, 4 (2013) 1911–1921.

9. K.M. Li, Y. Wang: On the three-dimensional sound fields from a moving monopole source above a non-locally reacting ground, Journal of the Acoustical Society of America 147, 4 (2020) 2581–2596.

10. K. Attenborough, T. Van Renterghem: Predicting outdoor sound, CRC Press, Boca Raton, 2021.

11. S.M. Candel: Numerical solution of conservation equations arising in linear wave theory: application to aeroacoustics, Journal of Fluid Mechanics 83, 3 (1977) 465–493.

12. A.D. Pierce: Acoustics: an introduction to its physical principles and applications, Springer International Publishing, 2019.

13. H. Bian, Q. Tan, S. Zhong, X. Zhang: Efficient computation of broadband noise propagation using Gaussian beam tracing method, Journal of the Acoustical Society of America 151, 5 (2022) 3387–3397.

14. L. Dallois, P. Blanc-Benon, D. Juvé: A wide-angle parabolic equation for acoustic waves in inhomogeneous moving media: applications to atmospheric sound propagation, Journal of Computational Acoustics 9, 2 (2001) 477–494.

15. V.E. Ostashev, J. Colas, D. Dragna, D.K. Wilson: Phase-preserving narrow- and wide-angle parabolic equations for sound propagation in moving media, Journal of the Acoustical Society of America 155, 2 (2024) 1086–1102.

16. M. Hornikx, R. Waxler, J. Forssén: The extended Fourier pseudospectral time-domain method for atmospheric sound propagation, Journal of the Acoustical Society of America 128 (2010) 1632–1646.

17. M. Hornikx, D. Dragna: Application of the Fourier pseudospectral time-domain method in orthogonal curvilinear coordinates for near-rigid moderately curved surfaces, Journal of the Acoustical Society of America 138, 1 (2015) 425–435.

18. D. Dragna, P. Blanc-Benon, F. Poisson: Modeling of broadband moving sources for time-domain simulations of outdoor sound propagation, AIAA Journal 52, 9 (2014) 1928–1939.

19. D. Dragna, P. Blanc-Benon: Towards realistic simulations of sound radiation by moving sources in outdoor environments, International Journal of Aeroacoustics 13, 5–6 (2014) 405–426.

20. H.H. Brouwer: A ray acoustics model for the propagation of aircraft noise through the atmosphere, International Journal of Aeroacoustics 13, 5–6 (2014) 363–383.

21. P. Schäfer, M. Vorländer: Atmospheric ray tracing: an efficient, open-source framework for finding eigenrays in a stratified, moving medium, Acta Acustica 5 (2021) 26.

22. F. Yunus, D. Casalino, F. Avallone, D. Ragni: Toward inclusion of atmospheric effects in the aircraft community noise predictions, Journal of the Acoustical Society of America 150, 2 (2021) 759–768.

23. C. Wu, S. Redonnet: Aircraft noise impact prediction with incorporation of meteorological effects, Transportation Research Part D: Transport and Environment 125 (2023) 103945.

24. M. Buret, K.M. Li, K. Attenborough: Optimisation of ground attenuation for moving sound sources, Applied Acoustics 67, 2 (2006) 135–156.

25. D. Dragna, P. Blanc-Benon: Sound radiation by a moving line source above an impedance plane with frequency-dependent properties, Journal of Sound and Vibration 349 (2015) 259–275.

26. Y. Wang, K.M. Li, D. Dragna, P. Blanc-Benon: On the sound field from a source moving above non-locally reacting grounds, Journal of sound and vibration 464 (2020) 114975.

27. V.E. Ostashev, D.K. Wilson, L. Liu, D.F. Aldridge, N.P. Symons, D. Marlin: Equations for finite-difference, time-domain simulation of sound propagation in moving inhomogeneous media and numerical implementation, Journal of the Acoustical Society of America 117, 2 (2005) 503–517.

28. M. Buret: New analytical Models for outdoor moving sources of sound, Doctoral Dissertation, Open University, Milton Keynes, 2002.

29. K. Attenborough, K.M. Li, K.V. Horoshenkov: Predicting outdoor sound, CRC Press, London, 2006.

30. V.E. Ostashev, D.K. Wilson: Acoustics in moving inhomogeneous media, CRC Press, London, 2015.

31. M. Roger, Sound radiation by moving surfaces and the Green's functions technique, in: R. Camussi (Ed.), Noise sources in turbulent shear flows: fundamentals and applications, Springer, Vienna, 2013, pp. 73–116.

32. P.M. Morse, K.U. Ingard: Theoretical acoustics, Princeton University Press, Princeton, 1968.

33. J. Scott, P. Blanc-Benon, O. Gainville: Weakly nonlinear propagation of small-wavelength, impulsive acoustic waves in a general atmosphere, Wave Motion 72 (2017) 41–61.

34. K. Attenborough, I. Bashir, S. Taherzadeh: Outdoor ground impedance models, Journal of the Acoustical Society of America 129, 5 (2011) 2806–2819.

35. C. Bogey, C. Bailly: A family of low dispersive and low dissipative explicit schemes for flow and noise computations, Journal of Computational physics 194, 1 (2004) 194–214.

36. J. Colas, A. Emmanuelli, D. Dragna, P. Blanc-Benon, B. Cotté, R.J.A.M. Stevens: Wind turbine sound propagation: comparison of a linearized Euler equations model with parabolic equation methods, Journal of the Acoustical Society of America 154, 3 (2023) 1413–1426.

37. R. Troian, D. Dragna, C. Bailly, M.-A. Galland: Broadband liner impedance eduction for multimodal acoustic propagation in the presence of a mean flow, Journal of Sound and Vibration 392 (2017) 200–216.

38. D. Komatitsch, R. Martin: An unsplit convolutional perfectly matched layer improved at grazing incidence for the seismic wave equation, Geophysics 72, 5 (2007) SM155–SM167.

39. ISO9613-1:1993: Acoustics – Sound attenuation in free field – Part 1: atmospheric absorption calculation, Technical report, International Standards Organization, Genève, 1993.
40. E. Salomons, D. van Maercke, J. Defrance, F. de Roo: The harmonoise sound propagation model, Acta Acustica united with Acustica 97, 1 (2011) 62–74.

41. O. Gainville: Modelisation de la propagation atmosphérique des ondes infrasonores par une méthode de tracé de rayons non-linéaire (“Numerical modelling of atmospheric infra-sound propagation using a nonlinear ray-tracing method”), PhD thesis No. 2008-07, Ecole Centrale de Lyon, Lyon, 2008.

42. D. Blokhintzev: The propagation of sound in an inhomogeneous and moving medium I, Journal of the Acoustical Society of America 18, 2 (1946) 322–328.

# 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 $\mathbf{x}_{\varphi}$, $\mathbf{x}_{\theta}$, $\mathbf{n}_{\varphi}$ and $\mathbf{n}_{\theta}$ in order to determine the amplitude along the rays. The equations for $\mathbf{x}_{\varphi}$ and $\mathbf{n}_{\varphi}$ write:

$$ \frac{d\mathbf{x}_{\varphi}}{dt_r} = (\mathbf{x}_{\varphi} \cdot \nabla c_0)\mathbf{n} + c_0\mathbf{n}_{\varphi} + (\mathbf{x}_{\varphi} \cdot \nabla)\mathbf{V}_{0}, \quad (A.1) $$

$$ \frac{d\mathbf{n}_{\varphi}}{dt_r} = (\mathbf{q}_{\varphi} \cdot \mathbf{n})\mathbf{n} + (\mathbf{q} \cdot \mathbf{n}_{\varphi})\mathbf{n} + (\mathbf{q} \cdot \mathbf{n})\mathbf{n}_{\varphi} - \mathbf{q}_{\varphi}, \quad (A.2) $$

with $\mathbf{q}_{\varphi} = (\mathbf{x}_{\varphi} \cdot \nabla)\nabla c_0 + [(\mathbf{x}_{\varphi} \cdot \nabla)\nabla\mathbf{V}_0] \cdot \mathbf{n} + (\nabla\mathbf{V}_0) \cdot \mathbf{n}_{\varphi}$. The equations for $\mathbf{x}_{\varphi}$ and $\mathbf{n}_{\varphi}$ are similar to equations (A.1) and (A.2) with the subscript $\varphi$ replaced by $\theta$.

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

$$ \frac{dx_{ri}}{dt_r} = c_0 n_i + V_{0i} \quad (A.3) $$

$$ \frac{dn_i}{dt_r} = q_j n_j n_i - q_i, \quad (A.4) $$

$$ \frac{dx_{\varphi i}}{dt_r} = x_{\varphi j} \frac{\partial c_0}{\partial x_j} n_i + c_0 n_{\varphi i} + x_{\varphi j} \frac{\partial V_{0i}}{\partial x_j}, \quad (A.5) $$

$$ \frac{dn_{\varphi i}}{dt_r} = q_{\varphi j} n_j n_i + q_j n_{\varphi j} n_i + q_i n_j n_{\varphi i} - q_{\varphi i}, \quad (A.6) $$

with the following expressions for the components of vectors $\mathbf{q}$ and $\mathbf{q}_{\varphi}$:

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

$$ q_{\varphi i} = x_{\varphi j} \frac{\partial^2 c_0}{\partial x_i \partial x_j} + n_j x_{\varphi k} \frac{\partial^2 V_{0j}}{\partial x_i \partial x_k} + n_{\varphi j} \frac{\partial V_{0j}}{\partial x_i}. \quad (A.8) $$

The ray variables are initialized with the relations:

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

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

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

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

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

$$ \mathbf{x}_{\theta}(t_r = 0) = 0. \quad (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 $\mathbf{e}_n$. For a flat horizontal ground, one has $\mathbf{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 = \mathbf{x}_r, \quad (A.15) $$

$$ \mathbf{n}' = \mathbf{n} - 2(\mathbf{n} \cdot \mathbf{e}_n)\mathbf{e}_n, \quad (A.16) $$

$$ \mathbf{x}'_{\varphi} = \mathbf{x}_{\varphi} - 2(\mathbf{x}_{\varphi} \cdot \mathbf{e}_n)\mathbf{e}_n, \quad (A.17) $$

$$ \mathbf{x}'_{\theta} = \mathbf{x}_{\theta} - 2(\mathbf{x}_{\theta} \cdot \mathbf{e}_n)\mathbf{e}_n, \quad (A.18) $$

$$ \mathbf{n}'_{\varphi} = \mathbf{n}_{\varphi} - 2(\mathbf{n}_{\varphi} \cdot \mathbf{e}_n)\mathbf{e}_n - \frac{\mathbf{x}_{\varphi} \cdot \mathbf{e}_n}{c_0 \mathbf{n} \cdot \mathbf{e}_n} \left[ \frac{d\mathbf{n}}{dt} - \frac{d\mathbf{n}'}{dt} - 2 \left( \frac{d\mathbf{n}}{dt} \cdot \mathbf{e}_n \right) \mathbf{e}_n \right], \quad (A.19) $$

$$ \mathbf{n}'_{\theta} = \mathbf{n}_{\theta} - 2(\mathbf{n}_{\theta} \cdot \mathbf{e}_n)\mathbf{e}_n - \frac{\mathbf{x}_{\theta} \cdot \mathbf{e}_n}{c_0 \mathbf{n} \cdot \mathbf{e}_n} \left[ \frac{d\mathbf{n}}{dt} - \frac{d\mathbf{n}'}{dt} - 2 \left( \frac{d\mathbf{n}}{dt} \cdot \mathbf{e}_n \right) \mathbf{e}_n \right]. \quad (A.20) $$

Note that $d\mathbf{n}/dt$ and $d\mathbf{n}'/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 $\mathbf{M}_s$ in a uniformly moving atmosphere with Mach number vector $\mathbf{M}_0$. First, we use the change of coordinates $\mathbf{x}' = \mathbf{x} - \mathbf{M}_0 c_0 t$ and $t' = t$. Equations (1) and (2) write:

$$ \frac{\partial p}{\partial t'} + \rho_0 c_0^2 \nabla' \cdot \mathbf{v} = c_0^2 S_0 e^{-i\omega_0 t'} \delta[\mathbf{x}' - (\mathbf{M}_s - \mathbf{M}_0)c_0 t'], \quad (B.1) $$

$$ \frac{\partial \mathbf{v}}{\partial t'} + \frac{1}{\rho_0} \nabla' p = 0. \quad (B.2) $$

where $\nabla' = \partial/\partial x'$. In the system of coordinates $(\mathbf{x}', t')$, the problem is equivalent to that of a monopole moving at Mach number $\mathbf{M}_s - \mathbf{M}_0$ in an homogeneous atmosphere at rest, for which an analytical solution is available. Back in the original system of coordinates $(\mathbf{x}, t)$, the acoustic pressure can be written as [31]:
$$p(\mathbf{x}, t) = - \frac{i\omega_0 e^{-i\omega_0 t}}{4\pi} \frac{e^{ik_0 R_e}}{R_e (1 - M_{s0} \cos \psi_e)^2} \quad (B.3)$$

where $R_e$ and $\cos\psi_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) \quad (B.4)$$

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

with $R = \mathbf{x} - \mathbf{M}_s c_0 t$, $R = |\mathbf{R}|$, $\cos\psi = |\mathbf{R} \cdot (\mathbf{M}_s - \mathbf{M}_0)|/(RM_{s0})$, $M_{s0} = |\mathbf{M}_s - \mathbf{M}_0|$, 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(t_r)c_0(t_r)} [1 + \mathbf{M}_0(t_r) \cdot \mathbf{n}(t_r)]^2 S(t_r) = \text{uste} = E(\varphi, \theta). \quad (B.6)$$

We can then deduce the pressure amplitude with:

$$A(t_r) = \sqrt{\frac{E(\varphi, \theta) \rho_0(t_r) c_0(t_r)}{S(t_r) [1 + \mathbf{M}_0(t_r) \cdot \mathbf{n}(t_r)]}} \quad (B.7)$$

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

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

$$\mathbf{x}_r(\tau) = \mathbf{x}_S(t_r) + (\mathbf{M}_0 + \mathbf{n})c_0\tau = \mathbf{x}_S(\tau) + (\mathbf{M}_0 - \mathbf{M}_s + \mathbf{n})c_0\tau, \quad (B.8)$$

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

$$\mathbf{n}_{\varphi}(\tau) = \mathbf{n}_{\varphi}(0), \quad (B.10)$$

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

$$\mathbf{x}_{\varphi}(\tau) = \mathbf{n}_{\varphi}(0)c_0\tau, \quad (B.12)$$

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

One can then deduce the raytube cross-section: $S(\tau) = |\mathbf{x}_{\varphi} \times \mathbf{x}_{\theta}| = (c_0\tau)^2|\cos\varphi|$.

The ray solution writes:

$$p(\mathbf{x}, t) = - \frac{i\omega_0 e^{-i\omega_0 t}}{4\pi} A(\tau) e^{i\omega_0 \tau}. \quad (B.14)$$

By identification with equation (B.3), one deduces that $R_e = c_0\tau$. The ray trajectory in equation (B.8) can be rewritten as $\mathbf{R} = (\mathbf{M}_0 - \mathbf{M}_s + \mathbf{n})R_e$. By projecting $\mathbf{R}$ on $\mathbf{M}_0 - \mathbf{M}_s$, one finds with equation (B.5) the equivalence $M_{s0}\cos\psi_e = \mathbf{n} \cdot (\mathbf{M}_s - \mathbf{M}_0)$.

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

$$A(\tau) = \frac{1}{c_0 \tau [1 + \mathbf{n} \cdot (\mathbf{M}_0 - \mathbf{M}_s)]^2} \quad (B.15)$$

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

$$E(\varphi, \theta) = \frac{|\cos \varphi|}{\rho_0 c_0} \frac{[1 + \mathbf{M}_0 \cdot \mathbf{n}]^2}{[1 + \mathbf{n} \cdot (\mathbf{M}_0 - \mathbf{M}_s)]^4} \quad (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).

Cite this article as: Kayser B. Dragna D. & Blanc-Benon P. 2024. Heuristic solution for the acoustic radiation of a moving monopole in an inhomogeneous and moving atmosphere. Application to aircraft noise. Acta Acustica, 8, 62. https://doi.org/10.1051/aacus/2024048.