A Time-Domain Finite-Difference Method for Bending Waves on Infinite Beams on an Elastic Foundation

Abstract

To model the vibration and structure-borne sound excitation and propagation of a railway rail, it can be modeled as an infinite beam on an elastic foundation. Existing analytical or numerical models are either formulated in the frequency domain or consider only finite beams in the time domain. Therefore, a time-domain approach for bending wave propagation on an effectively infinite beam on an elastic foundation is proposed. The approach makes use of an implicit finite-difference method that allows for varying properties of the beam and the foundation along the length of the beam. Strategies for an efficient discretization are discussed. The method is validated against existing analytical models for a single layer and two layers, as well as continuous and discrete support. The results show very good agreement, and it can be concluded that the proposed method can be seen as a versatile method for simulating the behavior of a beam on different kinds of elastic foundations.

1. Introduction

The vibration of a rail makes an essential contribution to railway rolling noise and has to be considered in any rolling noise model. One possible approach to modeling the dynamic behavior of a rail is a beam on an elastic foundation. This is appropriate as long as only frequencies below approximately 1.5 kHz [1] and the excitation and propagation of bending waves are considered. The source of railway rolling noise is the interaction of wheels and rails at their point of contact, depending on the time and location. The contact between a rolling wheel and a rail can be accounted for as dynamic moving load.

While a moving load can be seen as an inherently time-dependent process, most available models for moving loads on beams are formulated in the frequency domain. The models proposed in [2,3,4,5] may serve as examples in which the response of infinite beams on a continuous or discrete elastic foundation was considered in the frequency domain. Frequency-domain models for wave propagation are essentially linear models. However, aside from purely linear excitation and propagation of structure-borne noise in a rail, nonlinear effects may be also important for both rolling noise generation [6,7,8] and propagation [9].

Nonlinear effects may still be considered in frequency-domain models in some special cases in which appropriate transformations are available. For example, the influence of nonlinearity in the viscoelastic foundations on the dynamic response of a beam can be analyzed [9] by using such an approach. However, if a model takes into account multiple nonlinear effects, including those in the wheel–rail contact, it profits from an analysis in the time domain. A time-domain model for the dynamic response of a finite beam resting on a nonlinear viscoelastic foundation was reported in [10], where a Runge–Kutta method based on a Galerkin approximation was used. However, for a rail, it is more appropriate if it is considered as an (effectively) infinite beam.

The elastic foundation of a beam represents the pads, sleepers, and ballasts or different constructions supporting the rail. Depending on the railway track type considered and the fidelity of the model, a number of different approaches to modeling the elastic foundation are available. Continuous single- and double-layer foundations [8] were considered, as well as equidistant discrete supports (e.g., [7,11,12,13]). Some models also allow for variation of the properties, such as randomized spacing of the supports [3,14] or varying support stiffness [14].

Higher-fidelity models for bending waves in beams on elastic foundations can be obtained with the finite element method (FEM). Frequency-domain FEM approaches for the analysis of beams on elastic foundations can be found, e.g., in [1,15] (FEM use of periodic structure theory), [16] (waveguide FEM), and [17] (semi-analytical FEM). To use finite element models in the time domain, frequency-domain models are implemented with a time-stepping technique, e.g., [18,19], or with the moving element method, e.g., [20]. Thus, it is generally possible to include nonlinear effects, as shown for a finite beam [21,22] and for a continuously supported beam [23]. However, one problem of time-domain FEM analysis is the considerable computational cost connected with it. This limits its applicability, especially for parameter studies or comprehensive models where multiple wheels and their interactions with a realistic track are considered.

One alternative numerical method is the finite difference method (FDM). One advantage of this method is the inherent formulation in the time domain, which allows a simple consideration of the time-dependent behavior of the different components of the system. The governing wave equations for bending waves are of the fourth order and, thus, very different from the second-order differential equations that are considered in most applications of the FDM. While the basis for the application of the FDM to fourth-order differential equations is available [24,25], only a few applications for finite beams have been reported, e.g., in [26,27]. No application to infinite beams on elastic foundations has yet been documented. Therefore, the feasibility of the FDM approach in this case and the validity of the computed results remain open questions.

The present research addresses this question by introducing the finite difference method for bending waves on infinite beams on different types of elastic foundations in the time domain and analyzing the results. This can serve as a basis for considering arbitrary foundation properties that vary along the beam and for the inclusion of nonlinear effects of both excitation and foundation in the analysis. This enables the treatment of time-dependent interactions of the excitation force from the wheel–rail contact with rail vibrations, and even between multiple contacts. The present paper, however, restricts itself to non-moving loads acting on a beam on an elastic foundation.

The remainder of the paper is organized as follows. First, the beam model is presented. Then, the finite difference method is introduced and a technique for representing an infinite beam is explained. Finally, computations are compared with analytical results for four different cases of elastic foundations.

2. Model of a Beam on an Elastic Foundation

In the present paper, an infinite Euler–Bernoulli beam on an elastic foundation is considered, as shown in Figure 1. While more sophisticated beam models, such as the Timoshenko beam, may generally be more appropriate in some cases, this simple model already serves the purpose of demonstrating the applicability of the FDM. The following differential equation describes the bending vibrations of the Euler–Bernoulli beam on an elastic foundation [5,26]:

$$\begin{array}{c}B\frac{{\partial }^{4}u(x,t)}{\partial x^4}+m_r^{\prime }\frac{{\partial }^{2}u(x,t)}{\partial t^2}+d_r\left(x\right)\frac{\partial u(x,t)}{\partial t}=q(x,t)-F_s(x,t).\end{array}$$

(1)

The transverse deflection u depends on both x—the coordinate along the beam axis—and the time t. B is the bending stiffness, and $m_r^{\prime }$ is the beam mass per unit length. q represents the excitation force per unit length, and $d_r$ is the viscous damping coefficient per unit length of the beam. The force $F_s$ represents the effect of the foundation. Equation (1) is a fourth-order partial differential equation, and its solutions correspond to bending waves on the beam. Because (1) is just an approximation, its solutions may exhibit nonphysical behaviors. This is especially true for high frequencies, where the phase speed of the waves grows above every limit. This is a major challenge for the numerical treatment in the time domain. However, the same is also true for more sophisticated beam theories. As long as the bending wavelength is greater than six times the height of the beam, good agreement with more sophisticated beam theories can be obtained with the Euler–Bernoulli beam theory [28].

The effect of an elastic foundation (Figure 2) can be modeled as

$$F_s(x,t)=s_p\left(x\right)(u(x,t)-u_s(x,t))+d_p\left(x\right)\left(\frac{\partial u(x,t)}{\partial t}-\frac{\partial u_s(x,t)}{\partial t}\right)$$

(2)

where $s_p$ and $d_p$ are the stiffness and the viscous damping coefficient per unit length, respectively, and $u_s=0$ in the case of a single-layer elastic foundation. For a two-layer elastic foundation, $u_s$ is the transverse deflection of the intermediate layer and is given by

$$\begin{array}{c}{m}_s^{\prime }\frac{{\partial }^2{u}_s(x,t)}{\partial t^2}=s_p\left(x\right)u(x,t)-(s_b\left(x\right)+s_p\left(x\right))u_s(x,t)+d_p\left(x\right)\frac{\partial u(x,t)}{\partial t}\hfill \\ \hfill -(d_p\left(x\right)+d_b\left(x\right))\frac{\partial u_s(x,t)}{\partial t}\end{array}$$

(3)

with the mass per unit length $m_s^{\prime }$ of the intermediate layer, the stiffness per unit length $s_b$, and the viscous damping coefficient per unit length $d_b$ of the lower elastic layer. The viscous damping coefficients and stiffnesses per unit length may vary with the location. Whereas for continuously supported beams, they are constant along the length, discrete supports can be modeled by setting them to zero at all locations, except at the supports. This is illustrated in Figure 3 for the case of equally spaced supports that are a distance $L_s$ apart.

With the system of differential equations resulting from (1)–(3), the beam deflections can be calculated for both the single-layer and two-layer elastic foundation and continuous or discrete support. This analysis leads to some characteristic frequencies that determine the frequency response of these systems. For the single-layer supported beam, the frequency

$${\omega }_0=\sqrt{\frac{s_p}{m_r^{\prime }}}$$

(4)

corresponds to the resonance frequency of a spring–mass system when the mass of the beam per unit length is supported by the stiffness per unit length. This frequency is a cut-on frequency, as free wave propagation occurs only for frequencies above ${\omega }_0$. For frequencies well above this value, the wave propagation tends to be the same as that for an unsupported beam [8].

For the two-layer support, two resonance frequencies are found. The frequency ${\omega }_1=\sqrt{\frac{s_b}{m_s^{\prime }}}$ characterizes the resonance of the mass of the intermediate layer on the stiffness $s_b$. An anti-resonance occurs at ${\omega }_2=\sqrt{\frac{s_p+s_b}{m_s^{\prime }}}$. Cut-on frequencies also exist for the two-layer supports:

$${\omega }_{c1/2}=\sqrt{\frac{{\omega }_0^2+{\omega }_2^2}{2}\pm \sqrt{\frac{{({\omega }_0^2+{\omega }_2^2)}^2}{4}-{\omega }_0^2{\omega }_1^2}}.$$

(5)

Below ${\omega }_{c1}$, there is no wave propagation. For frequencies ${\omega }_{c1}<\omega <{\omega }_2$, free wave propagation can be observed. A blocked region for ${\omega }_2<\omega <{\omega }_{c2}$ follows. For frequencies above ${\omega }_{c2}$, the wave propagation also tends to that on an unsupported beam [8].

In the frequency domain, damping effects are often considered by using a complex-valued stiffness $s(1+j\eta )$, where $\eta$ is the damping loss factor, which is usually taken with no weak frequency dependence. In the time domain, however, it is not feasible to use complex-valued quantities. Therefore, the loss factor is converted into the viscous damping coefficient by using

$$d=\frac{\eta s}{\omega }.$$

(6)

Only for a monofrequent response at frequency $\omega$, the loss factor and the viscous damping coefficient result in fully identical damping effects [8].

3. Finite Difference Method for Infinite Beams on Foundations

To determine the beam deflections $u(x,t)$ from (1), (2), and (3) by using finite differences, the beam needs to be discretized into N equally long segments. The number of nodes is $N+1$ with a local step size $\Delta x$, as shown in Figure 4. The time is also discretized into discrete intervals $\Delta t$. Hence, the solution domain $x,t$ is covered with a rectangular grid with uniform spacing $\Delta x$ and uniform spacing $\Delta t$.

Equations (1) and (2) can be rewritten as

$$\begin{array}{c}\frac{{\partial }^{2}u(x,t)}{\partial t^2}+\frac{(d_r\left(x\right)+d_p\left(x\right))}{m_r^{\prime }}\frac{\partial u(x,t)}{\partial t}-\frac{d_p\left(x\right)}{m_r^{\prime }}\frac{\partial u_s(x,t)}{\partial t}=\\ -\frac{B}{m_r^{\prime }}\frac{{\partial }^{4}u(x,t)}{\partial x^4}-\frac{s_p\left(x\right)}{m_r^{\prime }}u(x,t)+\frac{s_p\left(x\right)}{m_r^{\prime }}{u}_s(x,t)+\frac{q(x,t)}{m_r^{\prime }}\end{array}$$

(7)

This equation contains a fourth derivative in space, as well as first and second derivatives with respect to time. In the FDM, these derivatives are approximated using finite differences (an overview of the FDM is available in [29]). For the fourth derivative in space in (7), this yields [24]

$$\begin{array}{c}\hfill {\left.\frac{{\partial }^4}{\partial x^4}u(x,t)\right|}_{x=i\Delta x}\approx \frac{1}{\Delta x^4}(u_{i-2}\left(t\right)-4u_{i-1}\left(t\right)+6u_i\left(t\right)-4u_{i+1}\left(t\right)+u_{i+2}\left(t\right)).\end{array}$$

(8)

Here and in the following, the lower indices denote the numbers of the grid points such that $u_i\left(t\right)=u(i\Delta x,t)$. In the matrix representation, (8) can be rewritten as

$$\begin{array}{c}\hfill \frac{{\partial }^4}{\partial x^4}u(x,t)\approx \frac{1}{\Delta x^4}\underline{\underline{D}}\underline{u}\end{array}$$

(9)

with the $(N+1)\times (N+1)$ pentadiagonal matrix and the $N+1$ vector

$$\underline{\underline{D}}=\left(\begin{array}{cccccc}\hfill 6& \hfill -4& \hfill 1& \hfill & \hfill \cdots & \hfill 0\\ \hfill -4& \hfill 6& \hfill -4& \hfill 1& \hfill & \hfill ⋮\\ \hfill 1& \hfill -4& \hfill 6& \hfill -4& \hfill 1& \\ \hfill & \hfill \ddots & \hfill \ddots & \hfill \ddots & \hfill \ddots & \hfill \ddots \\ \hfill ⋮& \hfill & \hfill 1& \hfill -4& \hfill 6& \hfill -4\\ \hfill 0& \hfill \cdots & \hfill & \hfill 1& \hfill -4& \hfill 6\end{array}\right)\mathrm{and}\underline{u}=\left(\begin{array}{c}{u}_0\hfill \\ u_1\hfill \\ ⋮\hfill \\ u_{N-1}\hfill \\ u_N\hfill \end{array}\right),$$

(10)

respectively. The time derivatives can be approximated in a similar fashion:

$$\begin{array}{cc}\hfill \frac{\partial }{\partial t}u(x,t)& \approx \frac{1}{\Delta t}(u^{n+1}-u^n)\hfill \end{array}$$

(11)

$$\begin{array}{cc}\hfill \frac{{\partial }^2}{\partial t^2}u(x,t)& \approx \frac{1}{\Delta t^2}(u^{n+1}-2u^n+u^{n-1})\hfill \end{array}$$

(12)

In this notation, the upper indices define the time step such that $u{\left(x\right)}^n=u(x,n\Delta t)$.

There are various possible approaches to solving the finite difference version of (7), where all derivatives are replaced by (9), (11), and (12). In this paper, the implicit Crank–Nicolson scheme was chosen because it has been shown that it yields stable solutions for the Euler–Bernoulli beam wave equation and positive values of B and $m_r^{\prime }$ [24,26]. The general idea of the Crank–Nicolson scheme is to replace all parts except the time derivatives at instant n with a weighted sum of its values at different instants. Thus, an implicit formulation arises where a system of equations has to be solved to produce the values of u at the next instant [30]. In the present case, these terms on the right-hand side of (7) are evaluated at instants $n+1$ and $n-1$, which leads to

$$\begin{array}{c}\frac{1}{\Delta t^2}(u_i^{n+1}-2u_i^n+u_i^{n-1})+\frac{(d_{ri}+d_{pi})}{m_r^{\prime }}\frac{1}{\Delta t}(u_i^{n+1}-u_i^n)-\frac{d_{pi}}{m_r^{\prime }}\frac{1}{\Delta t}(u_{si}^{n+1}-u_{si}^n)=\hfill \\ -\frac{B}{m_r^{\prime }}\frac{1}{\Delta x^4}\frac{1}{2}(u_{i-2}^{n+1}-4u_{i-1}^{n+1}+6u_i^{n+1}-4u_{i+1}^{n+1}+u_{i+2}^{n+1}\\ +u_{i-2}^{n-1}-4u_{i-1}^{n-1}+6u_i^{n-1}-4u_{i+1}^{n-1}+u_{i+2}^{n-1})\\ \hfill -\frac{s_{pi}}{m_r^{\prime }}\frac{1}{2}(u_i^{n-1}+u_i^{n+1})+\frac{s_{pi}}{m_r^{\prime }}\frac{1}{2}(u_{si}^{n-1}+u_{si}^{n+1})+\frac{q_i^n}{m_r^{\prime }}\end{array}$$

(13)

$d_{ri},d_{pi}$, and $s_{pi}$ indicate the stiffness and damping per unit length at node i. $q_i^n$ describes the force per unit length at the node i at time step n. The separation of time steps $n-1$, n, and $n+1$ results in the implicit representation of (7):

$$\begin{array}{c}\frac{B\Delta t^2}{2m_r^{\prime }\Delta x^4}\left(u_{i-2}^{n+1}-4u_{i-1}^{n+1}+6u_i^{n+1}-4u_{i+1}^{n+1}+u_{i+2}^{n+1}\right)+\hfill \\ \left(1+\frac{\Delta t}{m_r^{\prime }}(d_{ri}+d_{pi})+\frac{\Delta t^2}{2m_r^{\prime }}{s}_{pi}\right)u_i^{n+1}+\left(-\frac{\Delta t}{m_r^{\prime }}{d}_{pi}-\frac{\Delta t^2}{2m_r^{\prime }}{s}_{pi}\right)u_{si}^{n+1}=\\ -\left(1+\frac{\Delta t^2}{2m_r^{\prime }}{s}_{pi}\right)u_i^{n-1}-\frac{B\Delta t^2}{2m_r^{\prime }\Delta x^4}\left(u_{i-2}^{n-1}-4u_{i-1}^{n-1}+6u_i^{n-1}-4u_{i+1}^{n-1}+u_{i+2}^{n-1}\right)\\ \hfill +\left(2+\frac{\Delta t}{m_r^{\prime }}(d_{ri}+d_{pi})\right)u_i^n+\left(-\frac{\Delta t}{m_r^{\prime }}{d}_{pi}\right)u_{si}^n+\frac{\Delta t^2}{2m_r^{\prime }}{s}_{pi}{u}_{si}^{n-1}+\frac{\Delta t^2{q}_i^n}{m_r^{\prime }}\end{array}$$

(14)

which allows the computation of the result for instant $n+1$ from the values in the past at instants n.

For a two-layer elastic foundation, the finite difference version of (3) is needed in an implicit numerical form as well:

$$\begin{array}{c}\left(1+\frac{\Delta t}{m_s^{\prime }}(d_{pi}+d_{bi})+\frac{\Delta t^2}{2m_s^{\prime }}(s_{pi}+s_{bi})\right)u_{si}^{n+1}+\left(-\frac{\Delta t}{m_s^{\prime }}{d}_{pi}-\frac{\Delta t^2}{2m_s^{\prime }}{s}_{pi}\right)u_i^{n+1}=\hfill \\ \hfill \left(2+\frac{\Delta t}{m_s^{\prime }}(d_{pi}+d_{bi})\right)u_{si}^n+\left(-1-\frac{\Delta t^2}{2m_s^{\prime }}(s_{pi}+s_{bi})\right)u_{si}^{n-1}-\frac{\Delta t}{m_s^{\prime }}{d}_{pi}{u}_i^n+\frac{\Delta t^2}{2m_s^{\prime }}{s}_{pi}{u}_i^{n-1}.\end{array}$$

(15)

where $s_{bi}$ and $d_{bi}$ indicate the stiffness and damping coefficient per unit length of the intermediate layer at the node i, respectively.

To solve both (14) and (15) together for any grid point i, it is useful to join them into one system of equations:

$$\begin{array}{c}\left(\begin{array}{cc}\underline{\underline{A_{11}}}& \underline{\underline{A_{12}}}\\ \underline{\underline{A_{21}}}& \underline{\underline{A_{22}}}\end{array}\right){\left(\begin{array}{c}\underline{u}\\ {\underline{u}}_s\end{array}\right)}^{n+1}=\left(\begin{array}{cc}\underline{\underline{B_{11}}}& \underline{\underline{B_{12}}}\\ \underline{\underline{B_{21}}}& \underline{\underline{B_{22}}}\end{array}\right){\left(\begin{array}{c}\underline{u}\\ {\underline{u}}_s\end{array}\right)}^n+\hfill \\ \hfill \left(\begin{array}{cc}\underline{\underline{C_{11}}}& \underline{\underline{C_{12}}}\\ \underline{\underline{C21}}& \underline{\underline{C_{22}}}\end{array}\right){\left(\begin{array}{c}\underline{u}\\ {\underline{u}}_s\end{array}\right)}^{n-1}+{\left(\begin{array}{c}\frac{\Delta t^2\underline{q_i}}{m_r^{\prime }}\\ 0\end{array}\right)}^n.\end{array}$$

(16)

The sub-matrices contained here can be found in Appendix A. Note that all sub-matrices are diagonal, except $\underline{\underline{A_{11}}}$ and $\underline{\underline{C_{11}}}$, which are pentadiagonal. This is an advantageous property for the quick solution of (16). The support mass deflection ${\underline{u}}_s$ and the force per unit length $\underline{q}$ are, respectively, in vector form, equivalent to $\underline{u}$ in (10). For a beam with a single-layer support, only (14) is necessary, omitting all terms containing $u_s$. Generally, all stiffness and damping values may be updated after each time step, depending on the deflection history up until that time step. This would enable the modeling of nonlinear behavior of the foundation. A similar technique could be implemented for a nonlinear force excitation. However, in the following analysis, it is assumed that both foundation and excitation exhibit strictly linear behavior.

For the solution procedure of the FDM, it can be assumed that at the outset, the deflection and its time derivative (the velocity) along the beam are given. With these initial conditions, the unknown deflections of time step $n+1$ can be calculated. In this way, it is possible to solve (16) for each forthcoming time step.

Given the grid spacing in space and time, a characteristic value $r=\frac{B\Delta t^2}{m_r^{\prime }\Delta x^4}$ can be defined for the present analysis. For a spacing $\Delta x$ that is too small, oscillations of the solution can occur even with the implicit FDM because the Crank–Nicolson scheme is not L-stable. The Crank–Nicolson scheme for parabolic difference equations has a good numerical accuracy only for small values of r ($r\le 6$) [31]. Therefore, the spatial spacing of the grid is related to the time step:

$$\begin{array}{c}\hfill \Delta x=b\sqrt[4]{\frac{B}{6m_r^{\prime }}}\sqrt{\Delta t}.\end{array}$$

(17)

The constant b must be valid ($b\ge 1$).

One problem with the application of the finite difference method to infinite systems is that an infinite number of grid points N would be required. Therefore, methods have to be developed to mimic the infinite properties of the beam with a finite number of grid points. One possible approach is to equip the boundaries of a finite system in a way in which all wave energy is absorbed once it travels toward the boundaries. Artificially increased damping near the boundaries is a possible method for accomplishing this. This damping absorbs the energy of the bending wave in such a way that virtually no part of the wave is reflected from the boundaries.

Since there should be no deflections at the boundaries ($u_0=0$/ $u_{N+1}=0$), no specific precaution for the formulation of a boundary condition needs to be taken in $\underline{\underline{D}}$. To apply the damping, the spatial domain is divided into boundary domains and a calculation domain; see Figure 5. In the boundary domains, the inner damping coefficient $d_r$ increases slowly towards the boundaries over the length $l_B$. To ensure that the damping of the boundary connects fluently to the damping of the computational domain without a leap, a function is needed whose derivative disappears when it equals zero. Therefore, functions of the form

$$\begin{array}{c}\hfill d_r\left(x_{bc}\right)=\frac{d_{r,bc}}{{l_B}^{\alpha }}{x}_{bc}^{\alpha }\end{array}$$

(18)

were chosen with exponents $\alpha >1$. $x_{bc}$ is the coordinate of the boundary domains and $l_B=n_B\Delta x$, $n_b$ is the number of grid points used for the boundary domains (see Figure 5). $d_{r,bc}$ is the maximum value of the damping coefficient in the boundary domains.

The parameters $\alpha$, $l_B$, and $d_{r,bc}$ should be defined to accurately simulate the properties of an infinite beam while requiring less computing time. Therefore, $l_B$ should be as small as possible because, otherwise, the increases in the grid size N result in greater lengths of computing time. Therefore, a parametric study was conducted to determine these parameters. The parametric study revealed that the exponent $\alpha$ should be larger than 7 to obtain low values of $l_B$ (Figure A1). With the approach $d_{r,bc}=\frac{r}{2}\frac{m_r^{\prime }}{dt}$, good simulations of the infinite beam can be obtained over a wide range of different $m_r^{\prime }$ and, thus, different sound propagation velocities; see Figure A2.

If the analysis is carried out for a time interval $[0,t_{end}]$, in principle, the result allows the consideration of frequency components down to

$$f_{min,time}=\frac{1}{t_{end}}.$$

(19)

For other finite difference models of wave propagation, it was found that accurate solutions were obtained with $\Delta x⪅\lambda /10$, where $\lambda$ was the wavelength [32]. Considering ${\lambda }_{B_{min}}$ as the shortest bending wavelength of interest of a free infinite Euler–Bernoulli beam, the upper frequency limit determined by the spatial distance is given by:

$$\begin{array}{c}\hfill \Delta x⪅{\lambda }_{B_{min}}/g\end{array}$$

(20)

$$\begin{array}{c}\hfill f_{max,space}\le \frac{2\pi \sqrt{\frac{B}{m_r^{\prime }}}}{g^2\Delta x^2}.\end{array}$$

(21)

In the present analysis, sufficiently precise solutions were already found with the constant $g\ge 4$. Inserting Equation (17) into Equation (20) yields

$$\begin{array}{c}\hfill f_{max,space}\le \frac{2\pi \sqrt{6}}{b^2{g}^2\Delta t}.\end{array}$$

(22)

The upper frequency limit is also directly dependent on b. The constant b should preferably be set at the lowest possible value ($b=1$) for a broad frequency range, but this also increases the computational time.

The amplitude of the deflection is well approximated at high frequencies with 10 time steps per period. This results in the following upper frequency limit determined through time sampling:

$$\begin{array}{c}\hfill f_{max,time}\le \frac{1}{10\Delta t}.\end{array}$$

(23)

The lower frequency of Equation (21) or Equation (23) indicates the maximum frequency for the calculation procedure. Which of these two frequencies determines the upper limit depends on the choice of b and g.

4. Results and Discussion

In the following, the applicability of the proposed FDM procedure to rail-like applications is examined. Different beam models are set up by using the parameters given in

Table 1. Parameters of continuously supported beams.

Beam
Beam bending stiffness	B	6.42	MN/m2
Beam mass per unit length	$m_r^{\prime }$	60	kg/m
Viscous damping coefficient per unit length	$d_r$	0	Ns/m2
Continuous Support
Stiffness per unit length in the first layer	$s_p$	300	MN/m2
Stiffness per unit length in the second layer	$s_b$	100	MN/m2
Viscous damping coefficient per unit length in the first layer	$d_p$	30,000	Ns/m2
Viscous damping coefficient per unit length in the second layer	$d_b$	80,000	Ns/m2
Support mass per unit length	$m_s^{\prime }$	250	kg/m
Discrete Support
Stiffness of the first layer	$s_p$	180	MN/m
Stiffness of the second layer	$s_b$	60	MN/m
Viscous damping coefficients of the first layer	$d_p$	18,000	Ns/m
Viscous damping coefficients of the second layer	$d_b$	48,000	Ns/m
Support mass	$m_s^{\prime }$	150	kg
Support distance	$L_s$	0.6	m
Layer length	$L_l$	$\Delta x$

; these can be taken as representative for a railway rail; see [8] (Chap. 3). Unless otherwise specified, the FDM models all use the parameters from

Table 2. Numerical parameters.

Calculation end time	${\mathit{t}}_{\mathit{e}\mathit{n}\mathit{d}}$	0.5	s
Time increment	$\Delta t$	${10}^{-5}$	s
Grid parameter	b	1
Local step size	$\Delta x$	0.036	m
Boundary-domain exponent	$\alpha$	10
Number of boundary grid points	$n_B$	${10}^3$
Force parameter	$\sigma$	$0.7·{10}^{-4}$	s
Force parameter	a	$0.5·{10}^2$

. The overall length of the calculation domain is 80 m, and both boundary domains are chosen to have a length of 36.5 m, which results in an overall length of the rail of 153 m.

To investigate the propagation of bending waves in time and space, a broadband transient excitation was used. A pulse (first derivative of the Gaussian function) served as the force excitation $q(x,t)$, which excited the beam at a fixed location $x_F$ in the calculation domain (Figure 6):

$$\begin{array}{c}q(x,t)=a\frac{t_g}{{\sigma }^2}{e}^{\frac{-t_g^2}{{\sigma }^2}}\delta (x-x_F)\\ t_g=t-4\sigma \end{array}$$

(24)

The arbitrary parameters a and $\sigma$ govern the amplitude and duration of the pulse, respectively. $\delta$ is the Dirac delta.

In the case of a discretely supported beam with two layers, the computation required half a minute of CPU time on an Intel^® Core™i5-7200U CPU at 2.50 GHz for $5·{10}^4$ time steps.

Figure 7 shows the response of a beam on a single-layer elastic foundation at a distance of 10 m from the excitation point. The result is given for different time step sizes and associated spatial step sizes according to (17). The expected dispersive behavior of the bending wave can be observed. The high-frequency components and, thus, smaller wavelengths propagated faster than lower-frequency components. Therefore, bending waves with small wavelengths passed the receiver point $x_t$ first. Over time, the wavelengths increased. From the results, it turned out that with a smaller time step size (and smaller spatial step sizes), shorter-wavelength components could be computed more accurately. These short-wavelength components traveled faster due to the dispersion. Therefore, the first high-frequency wave components for a time step size of $5·{10}^{-6}$ s occurred earlier than those for the other, larger time step sizes. The energy originally introduced by the pulse was distributed over a larger frequency range for the smaller time step size than for the larger time step sizes. This explains why the maximum amplitudes were somewhat larger for the larger time steps. However, for lower-frequency components (and, thus, at a later time), the results are very similar to each other.

From a frequency-domain perspective, it would be interesting to know to the upper frequency at which a model produces reliable results. Analytical models are available that allow for a comparison. For a continuously supported beam on single-layer elastic foundation, the mobility $Y_{b,s}$ is given by [8]

$$\begin{array}{cc}\hfill Y_{b,s}\left(x_a\right)& =\frac{\omega }{4B{k}_{B,p}^3}(e^{-i{k}_{B,p}\left|x_a\right|}-i{e}^{-k_{B,p}\left|x_a\right|})\end{array}$$

(25)

$$\begin{array}{cc}\hfill k_{B,p}^2& =\sqrt{\frac{{\omega }^2{m}_r^{\prime }-s_p-j\omega d_p}{B}}\end{array}$$

(26)

where $x_a$ is the distance from the force point and $k_{B,p}$ is the bending wave number.

Figure 8 shows the results for the point mobility ($x_a=0$) at the excitation point as a function of frequency. To compute any mobility using the FDM method, the beam is excited at one point (24), and the response of the beam at one or more receiving points is computed. The response and the excitation signal are then transformed using a discrete Fourier transform. The mobility is estimated from the Fourier transforms $U\left(\omega \right)$ and $Q\left(\omega \right)$ using

$$Y=\omega \frac{U}{Q}.$$

(27)

For each time step size, the upper frequency limits could be estimated using (23). Figure 8 shows that for $\Delta t=2.5·{10}^{-5}$ s, the deviation from the analytical benchmark increased rapidly above the maximum frequency of 4 kHz. However, for shorter time step sizes and the associated smaller spatial step sizes, a much better agreement with the analytical solution could be observed. Since a time step size of $\Delta t=1·{10}^{-5}$ s gave satisfactory results in the frequency range of interest, this step size was used for all remaining computations.

Figure 9 shows the deflection of the beam after 0.05 s. The dispersive propagation that is visible in Figure 7 is superposed here by the resonance of the support. The decrease in the deflection amplitude in the boundary domain is obvious. As intended, the wave was damped and no reflection occurred.

To quantitatively analyze the effect of the proposed damped boundary domain method, the point mobility of an undamped, unsupported infinite beam was considered. It was given by [28]

$$\begin{array}{c}{Y}_b=\frac{\omega (1-j)}{2B{k}_B^3}\end{array}$$

(28)

$$\begin{array}{c}{k}_B^2=\sqrt{\frac{{\omega }^2{m}_r^{\prime }}{B}}.\end{array}$$

(29)

where $k_B$ is the bending wave number. Figure 10 compares this result to that from the FDM if no boundary domain is used and to that of the FDM with a damped boundary domain. Without boundary domains, the beam was effectively finite, and thus, the eigenfrequencies and corresponding bending mode shapes governed the results. With the proposed boundary domains, the results were identical to the analytical results, except for some small deviations for very low and very high frequencies.

Figure 11 shows further results for the effects of the boundary domains, but for a continuously supported beam. Two cases were compared with different damping in the elastic foundation. If the damping was high, the boundary domains had less of an effect. With low damping and no boundary domains, the influence from modal behavior was essential, but only above the cut-on frequency ${\omega }_0$. Below the cut-on frequency, there was no modal behavior, since there was no free wave propagation.

The use of the damped boundary domains led to high agreement with the analytical results. Especially for weakly damped systems, they are essential.

The last step in the validation of the FDM results was a comparison of the mobilities with analytical solutions for all four foundations types. While (25) gives the mobility in the case of a single-layer continuous foundation, the analytical results for a continuous two-layer support beam are given by [8]

$$\begin{array}{cc}\hfill Y_{b,ms}\left(x_a\right)& =\frac{\omega }{4B{k}_{B,pb}^3}(e^{-i{k}_{B,pb}\left|x_a\right|}-i{e}^{-k_{B,pb}\left|x_a\right|})\end{array}$$

(30)

$$\begin{array}{cc}\hfill k_{B,pb}^2& =\sqrt{\frac{{\omega }^2{m}_r^{\prime }-\overline{s}}{B}}\end{array}$$

(31)

$$\begin{array}{cc}\hfill \overline{s}& =\frac{{\overline{s}}_p({\overline{s}}_b-m_s^{\prime }{\omega }^2)}{{\overline{s}}_p+{\overline{s}}_b-m_s^{\prime }{\omega }^2}\end{array}$$

(32)

where ${\overline{s}}_p=s_p+j\omega d_p$ and ${\overline{s}}_b=s_b+j\omega d_b$. The length $x_a$ is again the distance to the excitation point. An analytical method for the mobilities in the case of discrete supports is given in [3] for the Euler–Bernoulli beam. The results from this model were used as benchmark for the FDM here.

The point mobility gives the response at the excitation point and is important for the characterization of the structure-borne sound excitation. Numerical and analytical results for the point mobility are shown in Figure 12. The propagation of structure-borne sound is more appropriately described by the transfer mobility, which relates the response at a location other than the excitation point to the excitation force. The results in Figure 13 show the transfer mobility for a point at a distance of 10 m from the force point. In the case of discrete support, the excitation point is located in the middle between the two supports.

The FDM results for both the point and the transfer mobility agreed very well with the analytical results in all cases. The resonance frequency for the single-layer support given by (4) at 356 Hz showed itself clearly in the results for the point mobility in Figure 12. The same was true for the anti-resonance ${\omega }_2$ at 201 Hz for the two-layer cases. The resonance ${\omega }_1$ at 100 Hz was not noticeable due to the stronger influence of the cut-on frequency ${\omega }_{c1}$ at 90 Hz. For the two-layer support, the stop bands below ${\omega }_{c1}$ and between ${\omega }_2$ and ${\omega }_{c1}$ at 400 Hz showed themselves in the significantly lower transfer mobility (Figure 13).

The discretely supported beams exhibited a sharp peak for the pinned–pinned frequency at 1414 Hz, which was approximately the frequency corresponding to half of the bending wavelength of the distance between the supports. This frequency was also accurately predicted with the FDM. Only for higher frequencies that were well above the frequency range of interest in rolling noise models could minor differences from the analytical solution from [3] be found.

5. Conclusions

The finite difference approach for an Euler–Bernoulli beam on an elastic foundation that is proposed in this paper is able to handle wave propagation not only for a finite length of the beam, but also for infinite beams. The method was demonstrated to yield very good results for different kinds of elastic foundations. It can, thus, be seen as a versatile method for simulating the behavior of a beam on an elastic foundation—for instance, a rail in the context of rolling noise modeling. While the present paper is only concerned with examples that can also be handled using analytical methods in the frequency domain, the application of the method extends beyond such a use. Because it is a numerical method with spatial discretization, the parameters of the beam and the support can be easily varied along the length of the rail. Because it works in the time domain, the moving excitation in the form of a nonlinear rail–wheel contact and the nonlinear properties of the elastic foundation can be considered without further simplifying assumptions. This is a clear advancement past the present state of the art, as detailed in the introduction. The limitation that results from the use of the Euler–Bernoulli beam theory can be addressed by extending the approach and using sophisticated beam theories as the basis for the numerical approach.