Sound Radiation from Railway Wheels including Ground Reflections: A half-space formulation for the Fourier Boundary Element Method

Current models for the acoustic radiation from railway wheels assume free field radiation. However, slab tracks are increasingly used for new railway lines. The acoustically hard surface of those tracks makes a re-evaluation of the free field assumption relevant, as such a surface can affect the radiation efficiency of an acoustic radiator. The wheel as the acoustic radiator is most conveniently described in a cylindrical coordinate system, thus making use of its axisymmetry. While this is a viable solution for the structural vibrations, for instance by using the curved Waveguide Finite Element formulation, the axisymmetry breaks when including a reflective plane in the calculation of the acoustic radiation. A convenient method to include an infinitely large, reflective plane is by using half-space Green’s functions in combination with the Boundary Element method. This method can be formulated in cylindrical coordinates using the Fourier series BEM (FBEM). However, the FBEM has not yet been combined with half-space Green’s functions. This paper provides a half-space formulation for the FBEM, which enables e.g. the evaluation of sound radiation of railway wheels over reflective surfaces. Finally, it is shown that the assumption of free field radiation for railway wheels is valid, as there is no major contribution of the reflective plane to the radiation efficiency of the wheel. The developed method is validated against laboratory measurements as well as analytical models. © 2020 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY license ( http://creativecommons.org/licenses/by/4.0/ )


Introduction
A solid has axial symmetry when it can be created by rotating a planar geometry around an axis. The axisymmetry of objects is often used to downscale numerical models, as simplifying the geometry to its planar representation considerably decreases the degrees of freedom in the system. The comparatively more elaborate element formulation of axisymmetric elements pays off in decreased calculation times. This downscaling is utilised in the following for calculating both structural vibrations as well as sound radiation. Standard axisymmetric finite elements are well established, see e.g. [1] , and part of most standard finite element packages [2,3] .
However, there is a large body of research using a curved waveguide finite element (WFE) formulation for axisymmetric bodies, especially in the field of predicting tyre vibrations and noise [4][5][6][7] . A summary of the curved WFE method can be

Method
In the following, the waveguide finite element method as well as the boundary element method for axisymmetric bodies are summarised. A half-space formulation for solving the acoustic radiation problem in Fourier-expanded cylindrical coordinates is introduced. Considerations regarding its efficiency are presented and a validation of the method is shown.

WFE formulation
A summary of the method is stated here for convenience. A comprehensive overview can be found in [6,8] . Describing the three-dimensional geometry by its constant cross-section enables the use of a two dimensional element formulation and drastically reduces the degrees of freedom in the system. The cross-section is defined in the (x, r) -plane with the radial direction r. The behaviour in the third dimension, which in case of an axisymmetric body refers to the circumferential direction θ , is described by propagating, decaying waves.
The displacement u at any location on the discretised body is described by u (x, r, θ ) = N (ξ , η) u (θ ) (1) with N (ξ , η) a matrix whose entries are shape functions N i and u containing all nodal displacements. For elastic materials without viscosity, Hamilton's principle is written δL = δ(U − K + A ) = 0 (2) with the Lagrangian L, the total strain energy U, the kinetic energy K and the potential energy of the loading A, where δ denotes "the first variation of". For a harmonic motion, Eq. (2) can be rewritten as with the volume of the structure V, the material stiffness matrix D and density ρ. The variable ω represents the circular frequency and f is an external force density. The superscript H denotes the complex conjugate transpose.
The material strain vector can be expressed on element e in terms of the nodal displacements u e as shown in Appendix A . Using this definition, Hamilton's principle can be reformulated to with the stiffness and mass matrix A mn and M and the weighted force vector f , respectively, which are obtained by solving a e,mn , m e and f e for each element and assembling to the global matrices and vector and with S e the surface of the element. The integral and the 'first variation of' vanishes after integration by parts of δL with respect to the θ -coordinate. The equation of motion is expressed in spatial and wavenumber domain as follows, (10) with K 2 = −A 11 , K 1 = A 01 − A 10 and K 0 = A 00 . Since only integer wavenumbers can propagate in circumferential direction, it is possible to prescribe the wavenumber κ and solve the resulting linear eigenvalue problem for ω. Physically, this corresponds to solving for the cross-sectional modes of the body for each circumferential order. Assuming the equation system is solved for N circumferential orders and S cross-sectional modes are obtained for each order, then each mode at eigenfrequency ω n,s is described by the corresponding eigenvector U n,s . The nodal displacement u for each circumferential order can be obtained by modal superposition [6,Ch. 4] , with the loss factor η and the modal mass m n,s δ s,t m n,s = U H n,s MU n,t (12) in which δ s,t is the Kronecker delta.
The nodal displacements of the FE-nodes at the boundary of the structure can be projected on the boundary normal direction at each node. Multiplication with j ω produces the surface normal velocity which serves as the input to the boundary element calculation.

FBEM Formulation
The Kirchhoff-Helmholtz integral equation describes the relation of the pressure amplitude on a boundary to the pressure field in a surrounding fluid. It is the basis for the Boundary Element Method (BEM) and is expressed as Fig. 1. Axisymmetric body over a reflecting plane. R is the path length between the source position P s and the observer position P obs . Notice the total reflected path length is R = R 1 + R 2 . The surface S is created by rotating the generator L around the symmetry axis. z plane is the z-coordinate of the intersection of the plane and with z-axis. The z-axis is oriented normal both to the plane and the symmetry axis.
with the Green's function of the considered problem G (r ; r 0 ) , the pressure p( r 0 ) and surface normal velocity v n ( r 0 ) at a source point P s on the boundary S, the density of the fluid ρ. The symbol ∂ n represents the derivative with respect to the outgoing normal direction on the fluid boundary n . The coefficient c(r ) is 1 / 2 for an observer point P obs on a smooth boundary, 1 in the fluid and 0 otherwise.
The Fourier Boundary Element Method (FBEM) utilises the axisymmetry of a body similar to WFEM, assuming that the sound field in circumferential direction can be described by a superposition of circumferential orders and therefore expanded in Fourier series. However, the reflection from a non axisymmetric plane is not easily described in this cylindrical coordinate system. The following sections present the necessary steps leading to the numerical implementation of half-space Green's functions in cylindrical coordinates described by the basis vectors (x b , r b , θ b ) .

FBEM formulation of the physical quantities
Consider the arbitrary axisymmetric body shown in Fig. 1 radiating over an infinite plane with a given reflection coefficient R p .
In cylindrical coordinates the half space Green's function G hs is with the first term corresponding to the free field Green's function and the second term describing the reflected field. Notice that this simplified formulation of the half-space Green's function is only a valid approximation for three special cases: (i) an infinite normal impedance at the boundary, i.e. R p = 1 , and (ii) for a perfectly soft boundary, i.e. R p = −1 [9] . The trivial case (iii) R p = 0 describes the situation in which no boundary is present. The distances R and R are written and R = ( x − x 0 ) 2 + ( r − r 0 ) 2 + 2 rr 0 ( 1 + cos ( θ r + θ r 0 ) ) + 4 z plane r cos ( θ r ) + r 0 cos ( θ r 0 ) + z plane (17) with the wavenumber k in circumferential direction.
The surface integral in Eq. (14) becomes d S 0 = r 0 d l 0 d θ r 0 , where d l 0 represents a line segment on the generator L . The necessary circumferential periodicity of all variables leads to the Fourier expansion in which the subscript * denotes that P obs can be located both on the boundary and in the domain. The reflected field G r (r , r 0 ) can not be directly expanded on θ r − θ r 0 . A 2D Fourier Transform has to be applied on G r (r , r 0 ) such that G hs (r , r 0 ) with the coefficients H m and H r nm , which are dependent on x, r, x 0 and r 0 , Analogously, At this point, solving p(x, r, θ r ) is equivalent to solving P * m (x, r) for every order m . This means that the initial 3D problem can be replaced by an infinite sum of 2D boundary problems if the the Fourier coefficient amplitudes are known.

Rewriting of the Kirchhoff-Helmholtz integral equation
where the orthogonality property between complex exponentials can be used. The integral over θ r 0 yields m = n, therefore 5 With that, the Fourier Boundary Integral equation can be expressed in a 3D half-space as which can be written as if the observer point P obs is located on a smooth boundary. The efficiency of the calculation of H m , H m , H r mp and H r mp is fundamental to this method and is discussed in the following subsection.

Fast fourier transform computations
H m and H m can be evaluated numerically by means of Gaussian quadrature. However, since they need to be computed for each order, Kuijpers [14] proposed a method using the Fast Fourier Transform (FFT) algorithm, enabling the computation of all necessary orders at once. An expression for the minimum required number of Fourier samples n FFT for the circumferential coordinate θ r − θ r 0 is proposed to evaluate the free field contribution G (r , r 0 ) in the Green's function. A relative error not exceeding 10 −3 is obtained in [14] with the expression , an oscillation criteria describing the impact of e −j kR for high frequencies and/or large difference (R max − R min ) .
R max and R min are, respectively, the maximum and the minimum distance between P obs (θ r ) and P s (θ r 0 ) when considering every possible difference (θ r − θ r 0 ) . Once G (r , r 0 ) is known for a given couple source/observer points for n FFT values, one can obtain H m by selecting the coefficient at the sought position in the FFT. The same process is applied to ∂ n G (r , r 0 ) to get This approach has been adapted for the calculation of the 2D Fourier Transform of G r (r , r 0 ) . Since G r (r , r 0 ) is of the same form as G (r , r 0 ) , the same criteria can be used to evaluate the minimum number of Fourier Transform samples for each angle ( n FFT ,r and n FFT , r 0 ) for a given accuracy. The steepness and oscillation criteria for the Fourier Transform over θ r are noted c s,r and c o,r while the one over θ r 0 are noted c s, r 0 and c o, r 0 . Fig. 2 shows the function G r (r , r 0 ) over both angular coordinates. It is apparent that e.g. for a transformation over θ r , the slope of G r (r , r 0 ) varies depending on the value of θ r 0 . Thus, the criteria are functions of θ r 0 and write c s,r (θ r 0 ) = R r, max (θ r 0 ) /R r, min (θ r 0 ) and c o,r (θ r 0 ) = k (R r, max (θ r 0 ) − R r, min (θ r 0 )) . R r, max (θ r 0 ) and R r, min (θ r 0 ) are, respectively, the relative maximum and minimum of R when varying θ r as shown in Fig. 3 . In order to obtain an expression for n FFT ,r (c s,r , c o,r ) = n FFT ,r (θ r 0 ) , G r (r , r 0 ) is computed for a large number of pairs of source and observer points so that c s,r ∈ [1 , 9] and orders from 0 to 500 using 5000 integration points. This computation is then used as a reference to determine the specific minimum number of Fourier Transform points n FFT ,r for each pair (c s,r , c o,r ) to achieve convergence. The condition for the evaluation of the convergence is a relative error not exceeding 10 −2 for all Fourier coefficients whose value is at least 1% of the maximum Fourier coefficient. These parameter triplets (n FFT ,r , c s,r , c o,r ) are finally used in a curve fitting procedure. The found expression is similar to the one found in [14] .
Since the dependency of G r (r , r 0 ) is the same on both angles, the expression of n FFT , r 0 (c s, r 0 , c o, r 0 ) has to be of the same form.
To enable the use of the fast radix-2 FFT algorithm, the next higher number that is a power of two should be used. Wider ranges for c s,r ( [1 , 1100] ) and c o,r ([0,2000]) have been studied afterwards to evaluate the robustness of n FFT ,r and n FFT , r 0 . The expressions of n FFT ,r and n FFT , r 0 never underestimate the necessary number of Fourier Transform samples. However, when the steepness criterion is large, (e.g. larger than 50, which would be the case when a large body is located very close to the ground or when the observer point is far from the ground), n FFT ,r and n FFT , r 0 start to overestimate the real value. A more complex expression has been developed to limit this overestimation, Eq. (29) is thus optimised for solving on the boundary of an object that is under 2 m in diameter and not much closer than one centimetre to the ground. When solving for field points, Eq. (30) is likely more computationally efficient. Fig. 4 represents n FFT ,r (θ r 0 ) and n FFT , r 0 (θ r ) for an arbitrarily chosen pair of source and observer point. This case is a representative example to show how to compute the 2D Fourier transform by computing consecutive 1D Fourier transforms. This proved to be computationally advantageous in this case due to reduced memory usage. The steps leading to the Fourier coefficients for this case are: • find min ( max (n FFT ,r (θ r 0 )) , max (n FFT , r 0 (θ r ))) to select first transformation direction. Here (see Fig. 4 ) n FFT ,r = 52 for θ r 0 = 0 , so select θ r • create θ r from max (n FFT ,r (θ r 0 )) equally spaced values in [0 , 2 π [ , Fig. 4. n FFT for different angles for the chosen example of source and receiver points.
• compute n FFT , r 0 (θ r ) for each value of θ r and create the corresponding θ r 0 , • evaluate G r (r , r 0 ) for each value of θ r and its corresponding values of θ r 0 , • calculate the Fourier transform of G r (r , r 0 ) over θ r 0 and transform it over θ r .
Singularities can occur in the computation of H m and H m . The integration points on the line integral 27 need to be deliberately chosen to avoid this. In addition to this, a convergence study about the combination of the Fourier transform and the Gaussian quadrature for singular kernels has been made via a regularisation method, as shown in Subsection 2.2.5 .

Discretisation process
The generator is discretised using quadratic line elements. Any point on the generator can be described by the following coordinates where N is a matrix containing the second order iso-parametric shape functions, x e and r e contain the nodal coordinates of the element to which r 0 belongs and ξ is the local element coordinate. P m (x 0 , r 0 ) and V m (x 0 , r 0 ) can now be expressed by the same formulation with p e m and v e m containing the nodal values of the corresponding element e . The line segment d l 0 is expressed as a function of d ξ with the Jacobian J e . With that, the element formulation in Eq. (27) can be written in the discretised form The vectors p m , p p , v m and v p collect the nodal values on the generator. Assembling the entries of the element matrices

G e m (r ) , H e m (r ) , G r e mp (r ) and H r e mp (r ) produces the global matrices G m (r ) , H m (r ) , G r mp (r ) and H r mp (r )
. P m at any point of the generator is now linked to its amplitude at every node of the generator. A collocation scheme is used, with P obs being placed successively on each node on the generator. Eq. (35) is written successively for each location, which leads to the system of equations where each row of G m is equal to G m (r ) evaluated at the corresponding node, and likewise for H * m , G r mp and H r mp . C is a diagonal matrix collecting the c coefficients. Since Eq. (36) yields a coupling between the Fourier coefficients P m , it needs to be written for every circumferential order m = [1 · · · M] which results in a solvable system of equations with H mm = C + H * m + H * r mm and G mm = G m + G r mm . The pressure on the boundary can now be solved from Eq. (37) inserting v m for m = [1 · · · M] as computed by the WFEM. Once the pressure is solved on the boundary, field points evaluation can be proceeded from 8 Since acoustic modes are coupled due to the reflective plane, the degree of accuracy in the computation of each one of them depends on the total number of acoustic modes solved for. This means that even if only M circumferential structural modes are excited, it might be necessary to solve for more than M circumferential acoustic modes. Thus a convergence study should be done to determine the number of circumferential acoustic modes necessary for a desired precision. This can be done simply by adding null vectors in the right hand side of Eq. 37 , which corresponds to adding Fourier coefficients of the normal velocity equal to zero.

Regularisation of singular integrals
Singular integral kernels occur in the BE formulation during the collocation process on each element where the observer point is placed. Integrating over such kernels leads to R → 0 in Eq. 15 . By selecting an even number of Gauss points in the quadrature, the function is not evaluated at a singularity, instead the numerical approximation of the integral approaches the correct value with an increasing number of Gauss points. To test this, every singular kernel in the collocation scheme has been regularised by introducing a constant as follows: The radiated sound power has been evaluated for a monopole source using the regularised formulation and compared to the analytical model. The difference between the analytical and the BE model has been evaluated for 4, 10 and 20 points in the Gaussian quadrature, for different values of . This is shown in Fig. 5 . For all three cases, the solution converges for decreasing . The difference is, expectedly, higher for a smaller number of Gauss points in the integration scheme.
Precautions have to be taken also in the computation of G r mp (r ) and H r mp (r ) since singular integrals arise when solving for a field point of coordinates x = x 0 and r ∈ [ −2 z plane − r 0 , −2 z plane + r 0 ] . Fig. 6 illustrates such a situation, where two   singularities occur (marked by "o"). In theses cases, when combining the 2D Fourier Transform with the line integral over the generator, R → 0 in Eq. 15 . A regularisation procedure similar to the one shown for the free field singularities should be used.

Validation of the Method
The validation is conducted in two steps, a validation of the combined WFE and free-field FBEM method, and a validation of the half-space axisymmetric FBEM method. All the following results are obtained using, respectively, 4 and 20 Gauss points to compute regular and regularized singular integrals over the generator. For the first part, a directivity measurement on an axisymmetric metal disk was conducted. The disk with the outside diameter of 22 cm, a central hole with 6 cm diameter and 8 mm thickness was suspended in an anechoic environment as shown in Fig. 7 . A shaker was used to apply an axial excitation at its outer diameter. The force input and acceleration at the input position were measured. A microphone was then manually positioned in 5 degree increments around the structure, where the 0 to 180 degree axis corresponds to the rotational symmetry axis of the structure. Fig. 8 shows the measured and calculated input mobility for the structure. The agreement is considered sufficient. The energy input into the structure is high at the resonance frequencies of the body. Due to the contribution of the shaker noise to the sound pressure at all frequencies, which is not included in the model, deviations of the measured and simulated sound pressure level are found, especially in proximity to the shaker. Thus, the comparison focuses on the resonance frequencies, at which the sound radiation from the structure dominates. Fig. 9 shows the comparison of the sound pressure level for some resonance frequencies. In general, these match closely. However, some shifts in angle can be observed, which might be due to imprecisions in the manual positioning. The increasingly complex radiation patterns require a higher accuracy and resolution of the positions with increasing frequencies. The second part of the validation considers the half-space formulation. Here, a numerical solution is compared to analytical models. The researched geometry is a breathing sphere. A semi-circle with a = 1 m diameter functions as the generator. A reflective plane is introduced with a distance of 0.51 m from the axis of rotation. To test the effect of the reflection, a reflection factor of -1 is used, effectively phase-shifting the mirror-image by π rad. The effect is shown in Fig. 10 as function of the wavenumber in air relative to the diameter of the sphere ka . For low ka, the structure follows the predicted radiation efficiency of a dipole. This is expected considering the phase relation between the sphere and its mirror image as well as the large wavelength compared to the size of the radiator. For large ka, the influence of the reflection decreases. Thus, towards high ka, the radiation efficiency approaches that of the breathing sphere.

Application on a railway wheel
The presented methods are applied on a common railway wheel with a straight web. The wheel geometry is of type BA093 as e.g. used in the noise measurement car (SMW) of DB Systemtechnik, described in [25] . A medium worn profile is assumed with a rolling radius of 0.47 m.

Analysis of the structural dynamics
First, the structural response of the wheel is evaluated. Fig. 11 shows the geometry and discretisation using 9-node, isoparametric, quadratic elements. As in [20] , the extended geometry of the axle is neglected and a rigid connection on the wheel hub is assumed, (dashed rectangle in Fig. 11 ). The dispersion relation for the presented wheel is shown in Fig. 12 .
In the considered frequency range up to 10 kHz, 79 modes are found with up to 13 nodal diameters. Axial, radial and circumferential modes with an increasing number of nodal circles exist for each number of nodal diameters. The mobility at the contact location is also presented in Fig. 12 . Expectedly, the axial mobility is mainly determined by the axial modes with a low number of nodal circles. Analogously, radial mode shapes have a dominant influence on the radial mobility. However both directions are coupled due to the non-symmetry of the wheel. The vibration behaviour below 2 kHz is dominated by the axial modes with up to four nodal diameters and up to one nodal circles (see Fig. 13 a and 13 b, with the first radial mode shape ( Fig. 13 d) occuring at 1.3 kHz.

Radiation ratios of the free and half-space model
The radiated sound power is evaluated for a harmonic unit force at the contact node, in axial and radial direction, respectively. An analysis similar to [20] is performed by evaluating the radiation ratio . 13. The mode shapes of the modes labelled in Fig. 12 . The notation is ( n, m, a/r) with the number of nodal diameters n, the number of nodal circles m and the main direction of motion, axial or radial, respectively.   with the spatially ( ) and temporally ( ) averaged squared velocity v on the surface S, the impedance of a plane wave in air ρ 0 c 0 and the radiated sound power P rad . The radiation ratio has been evaluated for a wheel in a free field as well as for a half space setup with the lowest point of the wheel 20 cm above an acoustically rigid ground, representing for example the wheel on a rail over a slab track. Fig. 14 shows this radiation ratio for six orders for radiation in free field and excitation in axial direction. The low-frequency radiation ratio follows the function described in [20] , f 2 n +4 for order n .
The behaviour of these orders can therefore be approximated by the radiation of multipoles, with σ tending towards unity once the wavelength in air is in the same magnitude or smaller than the radiating object. The radiation ratio of six orders for a radial force input is shown in Fig. 15 . The first orders follow the multipole approximation described in [20] for axial modes, f 2 n +4 . This means that the radial force input produces a motion in the axial direction large enough to dominate the  radiation. Higher orders follow the predicted pattern of f 2 n +2 for low frequencies, meaning that for those, the radial force input leads to the dominant radiation by radial modes. Next, the half-space solution is analysed. A reflecting surface can be close to the outer diameter of the wheel, as it is the case for embedded rails for trams. For a slab track system, this distance between the top of the rail and the slab can be around 20 cm to 30 cm. An acoustically hard reflective plane ( R p = 1 ) is introduced at different heights below the wheel to quantify this effect. Due to the inherent coupling of wavenumbers in the half-space formulation, an analysis of the radiation ratio per order is not meaningful. However, the combined radiation ratio of all included nodal diameters can be compared. Fig. 16 shows this radiation ratio for a radial excitation. The influence of the plane on the radiation ratio is visible mainly at frequencies below about 170 Hz. Considering that here, the radiation ratio is below 0.01 and decreasing towards lower frequencies, the overall effect of the reflective plane is comparatively small.

Approximation of the radiation ratio by independent mirror sources
The half-space solution with a rigid reflective surface can alternatively be modelled by including a mirror source. For both the original and the mirror source, the sound pressure at a defined location can be calculated individually, without the other source present. Adding these independently calculated sound pressures can be an approximation of the solution of the coupled system. This approximation is compared to the coupled half-space and the free-space solution of the BE model. Fig. 17 shows this comparison. In the frequency range of up to about 150 Hz, this simplification provides a better model of the radiation ratio than only using the wheel in free space. The approximation gets closer to the half-space model for lower frequencies (note that the frequency axis is shown down to 20 Hz). Above about 150 Hz, the approximation matches the radiation ratio of the wheel both in half-space and free space.

Sound pressure level at field points
As seen in the previous sections, the radiation ratio is not affected significantly by the reflective plane. However, the sound pressure level at distinct points in the sound field can see a larger impact of the reflective plane. In the following, the sound pressure level created by a unit harmonic force at the wheel-rail contact point is evaluated at three positions next to the wheel. The points are chosen according to ISO 3095 [26] , assuming that the wheel is located on the rail close to the defined measurement positions. Fig. 18 shows the evaluation of the sound pressure level at the defined field points for the wheel in free and half-space as well as for the approximation by two independent sources as described in Section 3.3 . In general, the interference pattern of the sound radiated by the wheel and the reflective plane creates a sound field which is (as expected) different from that of a free wheel. For frequencies up to about 100 Hz, the expected difference of about 6 dB between the radiation in free and in half-space is visible. Here, using two independent sources for the model is an acceptable approximation. For positions (a) and (c) in the figure, the approximation is fairly close up to 300 Hz.

Conclusions
The purpose of this paper is to evaluate the influence of sound radiation from railway wheels over reflective planes. This is of practical interest when researching the influence of the acoustically hard surface of slab tracks on the source strength of wheel radiation. To examine this influence, the existing Fourier Boundary Element method (FBEM) has been expanded to comprise half space Green's functions. An axisymmetric formulation of the Waveguide Finite Element method (WFEM) has been used for the structural vibrations.
The proposed method of including half space Green's functions in the FBEM has been implemented and validated against numerical models. The combined model of WFEM for structural vibrations and FBEM for sound radiation showed accurate predictions when comparing to laboratory measurements of a real structure. The validated model then successfully reproduced the predictions on the radiation characteristics of a railway wheel in free field from [20] . Furthermore, it was shown that the overall radiation efficiency of a typical railway wheel is low below 150 Hz.
The effect of a reflecting plane on the source strength of a railway wheel was shown to be rather small. Consequently, the assumption used in literature (e.g. [20] ) is valid above about 150 Hz. These results confirm that the models currently used to estimate the source strength from railway wheels on ballasted track are also valid approximations for wheels on slab track and even tram wheels on embedded rails.
It is further found that, even though the radiation characteristic of the wheel remains largely unchanged when introducing the reflective plane, the interference pattern of the sound radiated by the wheel and the reflective plane creates a significantly different sound field which can be observed when measuring the pressure at distinct positions. A hybrid model introduced to approximate the half-space solution was able to produce acceptable results for the source strength but could not reliably predict the sound pressure level in the sound field.
A limitation of the proposed method is the range of researched geometries. The computational cost increases for large bodies that are close to the reflective plane. The presented algorithm has been optimized for bodies in the range of about 1 m diameter and frequencies up to 10 kHz.

Declaration of Competing Interest
None.