Amplitude Stochastic Response of Rayleigh Beams to Randomly Moving Loads

TL;DR

This paper analyzes the nonlinear stochastic response of Rayleigh beams subjected to randomly moving loads, providing analytical estimates and numerical validation of amplitude behavior under stochastic load velocities.

Contribution

It introduces an analytical framework combining Fourier transform, residues, and stochastic averaging to study nonlinear beam responses to stochastic moving loads.

Findings

Mean amplitude exhibits nonmonotonic behavior with load velocity variability.

Analytical results align well with numerical simulations.

Load number influences the system's response significantly.

Abstract

We consider the problem of the nonlinear response of a Rayleigh beam to the passage of a train of forces moving with stochastic velocity. The Fourier transform and the theory of residues is used to estimate the mean-square amplitude of the beam, while the stochastic averaging method gives the stationary probability density function of the oscillations amplitude. The analysis shows that the effect of the load random velocities is highly nonlinear, leading to a nonmonotonic behavior of the mean amplitude versus the intensity of the stochastic term and of the load weight. The analytic approach is also checked with numerical simulations. The effect of loads number on the system response is numerically investigated.

Tables

Table 1. Table 1: Coefficients of stochastic Fourth- order Runge Kutta.

Coefficient	Value	Coefficient	Value
${\alpha }_1$ ${\alpha }_2$ ${\alpha }_3$ ${\alpha }_4$ $a_{21}$ $a_{31}$ $a_{32}$	$0.25001352164789$ $0.67428574806272$ $-0.00831795169360$ $0.08401868181222$ $0.66667754298442$ $0.63493935027993$ $0.00342761715422$	$a_{41}$ $a_{42}$ $a_{43}$ $q_1$ $q_2$ $q_3$ $q_4$	$-2.32428921184321$ $2.69723745129487$ $0.29093673271592$ $3.99956364361748$ $1.64524970733585$ $1.59330355118722$ $0.26330006501868$

Equations

$$\begin{array}[]{l}{M_{b}}\frac{{{\partial^{2}}w\left({x,t}\right)}}{{\partial{t^{2}}}}-{R_{a}}\frac{{{\partial^{4}}w\left({x,t}\right)}}{{\partial{t^{2}}\partial{x^{2}}}}+C\frac{{\partial w\left({x,t}\right)}}{{\partial t}}+\\ \\ \,EI\frac{{{\partial^{2}}}}{{\partial{x^{2}}}}\left[{\frac{{{\partial^{2}}w\left({x,t}\right)}}{{\partial{x^{2}}}}\,\left({1-\frac{3}{2}\,{{\left({\frac{{\partial w\left({x,t}\right)}}{{\partial x}}}\right)}^{2}}}\right)}\right]=\\ \\ \,=P\sum\limits_{i=1}^{{N_{v}}}{{\varepsilon_{i}}\delta\left[{x-{x_{i}}(t-{t_{i}})}\right]}\end{array}$$

$$\begin{array}[]{l}\left.w\left({x,t}\right)\right|_{x=0,L}=0,\,\,\,\left.{\frac{{\partial^{2}w\left({x,t}\right)}}{{\partial x^{2}}}}\right|_{x=0,L}=0.\\ \\ \left.w\left({x,t}\right)\right|_{t=0}=0\,\,,\,\left.\frac{{\partial w\left({x,t}\right)}}{{\partial t}}\right|_{t=0}=0.\end{array}$$

$$\begin{array}[]{l}\frac{{d{x_{i}}(t-{t_{i}})}}{{dt}}={v_{i}}(t-{t_{i}})={v_{0}}+{\sigma_{v}}{\xi_{i}}(t-{t_{i}}),\\ \\ 0\leq{x_{i}}(t-{t_{i}})\leq L.\end{array}$$

$$\begin{array}[]{l}\langle{v_{i}}(t-{t_{i}})\rangle={v_{0}}\,,\,\,\,\,\langle{\xi_{i}}(t-{t_{i}})\rangle=0,\\ \\ \langle{\xi_{i}}(t-{t_{i}}){\xi_{j}}(t-{t_{i}})\rangle=0\,\,\,\,\,\,{\rm{for}}\quad i\neq j,\\ \\ \langle{\xi_{i}}(t-{t_{i}}){\xi_{i}}(t-{t_{i}}+\zeta)\rangle=\,\,\gamma_{v}^{2}\delta(\zeta).\end{array}$$

$$w\left(x,t\right)=\sum _{n=1}^{\infty }{q}_n\left(t\right)\sin \left(\frac{n\pi x}{L}\right)$$

$${\chi }_n=\frac{q_n}{r},\tau =\frac{\tilde{v}t}{L},v=\frac{v_0}{\tilde{v}}$$

$$\begin{array}[]{l}\ddot{\chi}(\tau)+\lambda\dot{\chi}(\tau)+\chi(\tau)+\beta{\chi^{3}}(\tau)=\\ =\Gamma\sum\limits_{i=1}^{{N_{v}}}{{\varepsilon_{i}}\sin\left[{\Omega\tau+\gamma W_{i}\left(\tau-\tau_{i}\right)}\right]}\end{array}$$

$$\begin{array}[]{l}\Omega=\frac{{\pi{v_{0}}}}{{\tilde{v}}}\,,\,\,\,\Gamma=\frac{{2P{L^{3}}}}{{rEI{\pi^{4}}}}\,\,\,\,,\,\,\,\beta=-\frac{3}{8}{\left({\frac{{\pi r}}{L}\,}\right)^{2}}\,\,,\\ \\ \lambda=\frac{{C{L^{3}}}}{{{\pi^{2}}\sqrt{EI\left({{L^{2}}{M_{b}}+{R_{a}}{\pi^{2}}}\right)}}}\,\,\,\,{R_{a}}={M_{b}}{r^{2}}\,\,,\\ \\ \tilde{v}={\pi^{2}}\sqrt{\frac{{EI}}{{{L^{2}}{M_{b}}+{R_{a}}{\pi^{2}}}}}\,\,,\,\,\,r=\sqrt{\frac{I}{S}}\,\,,\,\,\,\,\,\,\gamma=\pi\sigma_{v}\,\,,\\ \\ \tau_{i}=\frac{{\tilde{v}t_{i}}}{L}.\end{array}$$

$$\ddot{\chi }(\tau )+\lambda \dot{\chi }(\tau )+\chi (\tau )+\beta {\chi }^3(\tau )=\Gamma \sin \left[\Omega \tau +\gamma W\left(\tau \right)\right]$$

$$\begin{array}[]{l}{C_{\eta}}(\tau,\tau^{\prime})=\langle\eta(\tau)\eta(\tau^{\prime})\rangle=\\ \\ =\frac{{{\Gamma^{2}}}}{2}\exp\left({-\frac{{{\gamma^{2}}{|\tau-\tau^{\prime}|}}}{2}}\right)\cos\Omega(\tau-\tau^{\prime}).\end{array}$$

$$\ddot{\chi }(\tau )+\lambda \dot{\chi }(\tau )+{\omega }^2\chi (\tau )=\eta (\tau )$$

$${\omega }^2$$

$$\eta (\tau )$$

$$\begin{array}[]{l}TF\left({\chi(\tau)}\right)\,=\chi(\omega^{\prime})\,,\\ \\ TF\left({\eta(\tau)}\right)\,=\eta(\omega^{\prime})\,,\\ \\ TF\left({{\chi^{(n)}}(\tau)}\right)\,={(i\omega^{\prime})^{n}}TF\left({\chi(\tau)}\right)\end{array}$$

$$\chi ({\omega }^{\prime })=\frac{\eta ({\omega }^{\prime })}{{\omega }^2-{{\omega }^{\prime }}^2+i\lambda {\omega }^{\prime }}$$

$$⟨\eta ({\omega }^{\prime })\eta ({\omega }^{\prime \prime })⟩=\delta \left({\omega }^{\prime }+{\omega }^{\prime \prime }\right)S_{\eta }({\omega }^{\prime })$$

$$S_{\eta }({\omega }^{\prime })=\frac{{\left(\Gamma \gamma \right)}^2}{2\pi }\left[\frac{1}{4{\left({\omega }^{\prime }-\Omega \right)}^2+{\gamma }^4}+\frac{1}{4{\left({\omega }^{\prime }+\Omega \right)}^2+{\gamma }^4}\right]$$

$$\begin{array}[]{l}\langle\chi(\omega^{\prime})\chi(\omega^{\prime\prime})\rangle\,\,\,\,=\,\,\,\frac{\begin{array}[]{l}\langle\eta(\omega^{\prime})\eta(\omega^{\prime\prime})\rangle\\ \end{array}}{{\left({{\omega^{2}}-{{\omega^{\prime}}^{2}}+i\lambda\omega^{\prime}}\right)\left({{\omega^{2}}-{{\omega^{\prime\prime}}^{2}}+i\lambda\omega^{\prime\prime}}\right)}}\\ \\ \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,=\frac{\begin{array}[]{l}\delta\left({\omega^{\prime}+\omega^{\prime\prime}}\right)\,\,{S_{\eta}}(\omega^{\prime})\\ \end{array}}{{\left({{\omega^{2}}-{{\omega^{\prime}}^{2}}+i\lambda\omega^{\prime}}\right)\left({{\omega^{2}}-{{\omega^{\prime\prime}}^{2}}+i\lambda\omega^{\prime\prime}}\right)}}\end{array}$$

$$A^2=⟨{\chi }^2(\tau )⟩$$

$$\begin{array}[]{l}{A^{2}}\,\,\,=\,\int{\frac{{d{\omega_{1}}}}{{2\pi}}}\frac{{d{\omega_{2}}}}{{2\pi}}\prec\chi({\omega_{1}})\chi({\omega_{2}})\succ{e^{i\left({{\omega_{1}}+{\omega_{2}}}\right)\tau}}\\ \\ =\int{\frac{{d{\omega_{1}}}}{{2\pi}}}\frac{{d{\omega_{2}}}}{{2\pi}}\frac{{\delta\left({{\omega_{1}}+{\omega_{2}}}\right){S_{\eta}}({\omega_{1}}){e^{i\left({{\omega_{1}}+{\omega_{2}}}\right)\tau}}}}{{\left({{\omega^{2}}-\omega_{1}^{2}+i\lambda{\omega_{1}}}\right)\left({{\omega^{2}}-\omega_{2}^{2}+i\lambda{\omega_{2}}}\right)}}\\ \\ =\,\frac{{{{\left({\Gamma\gamma}\right)}^{2}}}}{{2\pi}}\int{\frac{{d{\omega_{1}}}}{{2\pi}}\frac{{2{\gamma^{4}}+8\left({\omega_{1}^{2}+{\Omega^{2}}}\right)}}{{\left[{4{{\left({{\omega_{1}}-\Omega}\right)}^{2}}+{\gamma^{4}}}\right]\left[{4{{\left({{\omega_{1}}+\Omega}\right)}^{2}}+{\gamma^{4}}}\right]\left[{{{\left({{\omega^{2}}-\omega_{1}^{2}}\right)}^{2}}+{\lambda^{2}}\omega_{1}^{2}}\right]}}}\\ \\ \,\,\,\,\,=\,\frac{{{{\left({\Gamma\gamma}\right)}^{2}}}}{{4\pi\lambda}}.\frac{{\left[{{\gamma^{4}}+4\left({{\omega^{2}}+{{\left({\lambda-\Omega}\right)}^{2}}}\right)}\right]}}{{{\omega^{2}}\left[{{\gamma^{8}}+16{{\left({{\omega^{2}}+\left({\lambda-\Omega}\right)\Omega}\right)}^{2}}+4{\gamma^{4}}\left({{\lambda^{2}}+2\left({{\omega^{2}}+{\Omega^{2}}-\lambda\Omega}\right)}\right)}\right]}}\end{array}$$

$${\Theta }_4{A}^8+{\Theta }_3{A}^6+{\Theta }_2{A}^4+{\Theta }_1{A}^2={\Theta }_0$$

$${\Theta }_0={\Gamma }^2\left[{\gamma }^6+4{\gamma }^2\left(1+{\left(\lambda -\Omega \right)}^2\right)\right]$$

$${\Theta }_1=4\pi \lambda \times$$

$$\times \left[{\gamma }^8+4\left({\lambda }^2+2\left(1+{\Omega }^2-\lambda \Omega \right)\right){\gamma }^4-\frac{3\beta {\Gamma }^2}{4\pi \lambda }{\gamma }^2\right]+$$

$$64\pi \lambda \Omega \left[2\lambda +\Omega \left({\Omega }^2+{\lambda }^2-2-2\lambda \Omega \right)\right]$$

$${\Theta }_2=3\lambda \pi \beta \left[{\gamma }^8+4\left(4+{\lambda }^2+2\Omega \left(\Omega -\lambda \right)\right){\gamma }^4\right]+$$

$$48\lambda \pi \beta \left[3+{\Omega }^4+\Omega \left(\Omega {\lambda }^2-2\lambda {\Omega }^2+4\left(\lambda -\Omega \right)\right)\right]$$

$${\Theta }_3=18\lambda \pi {\beta }^2\left[{\gamma }^4+6+4\Omega \left(\lambda -\Omega \right)\right]$$

$${\Theta }_4=27\lambda \pi {\beta }^3$$

$$\left\{\begin{array}[]{l}\dot{\chi}(\tau)=y(\tau)\\ \\ \dot{y}(\tau)=-\lambda y(\tau)-\chi(\tau)-\beta{\chi^{3}}(\tau)+\Gamma\sin\left({z(\tau)}\right)\\ \\ \dot{z}(\tau)=\Omega+\gamma\xi(\tau)\end{array}\right.$$

$$\begin{array}[]{l}{\chi_{k+1}}={\chi_{k}}+{\alpha_{1}}{K_{1}}+{\alpha_{2}}{K_{2}}+{\alpha_{3}}{K_{3}}+{\alpha_{4}}{K_{4}}\\ {y_{k+1}}={y_{k}}+{\alpha_{1}}{L_{1}}+{\alpha_{2}}{L_{2}}+{\alpha_{3}}{L_{3}}+{\alpha_{4}}{L_{4}}\\ {z_{k+1}}={z_{k}}+{\alpha_{1}}{M_{1}}+{\alpha_{2}}{M_{2}}+{\alpha_{3}}{M_{3}}+{\alpha_{4}}{M_{4}}\\ \\ {k_{1}}={y_{k}}\,\Delta\tau\\ {L_{1}}=\left[{-\lambda{y_{k}}-{\chi_{k}}-\beta\chi_{k}^{3}+\Gamma\sin\left({{z_{k}}}\right)}\right]\Delta\tau\\ {M_{1}}=\Omega\Delta\tau+{\zeta_{1}}\\ \\ {K_{j}}=\left[{{y_{k}}+{Y_{j}}}\right]\Delta\tau\\ {L_{j}}=[-\lambda\left({{y_{k}}+{Y_{j}}}\right)-\left({{\chi_{k}}+{X_{j}}}\right)-\beta{{\left({{\chi_{k}}+{X_{j}}}\right)}^{3}}+\\ +\Gamma\sin\left(z_{k}+Z_{j}\right)]\Delta\tau\\ {M_{j}}=\Omega\,\Delta\tau+{\zeta_{j}}\end{array}$$

$$\begin{array}[]{l}{X_{j}}=\sum\limits_{i=1}^{j-1}{{a_{j-1}}{K_{i}}\,\,,\,\,\,}{Y_{j}}=\sum\limits_{i=1}^{j-1}{{a_{j-1}}{L_{i}}\,\,,\,\,{Z_{j}}=\sum\limits_{i=1}^{j-1}{{a_{j-1}}{M_{i}}\,,}}\\ \\ {\zeta_{1}}={r_{1}}\sqrt{2{q_{1}}\gamma\Delta\tau}\,\,,\,\,\,{\zeta_{j}}={r_{j}}\sqrt{2{q_{j}}\gamma\Delta\tau}\,\,;\,\,j=2,\,3,\,4\end{array}$$

$$\begin{array}[]{l}\chi(\tau_{k})=\chi_{k}+O\left(\Delta{\tau^{5}}\right),\,\,y(\tau_{k})=y_{k}+O\left(\Delta{\tau^{5}}\right),\\ z(\tau_{k})=z_{k}+O\left(\Delta{\tau^{5}}\right)\end{array}$$

$$\begin{array}[]{l}\chi(\tau)=a(\tau)\cos\theta,\,\,\,\,\,\,\,\,\,\dot{\chi}(\tau)=-\Omega a(\tau)\sin\theta,\\ \theta=\Omega\tau+\varphi(\tau).\end{array}$$

Full text

††journal: To appear on Nonlinear Dynamics

1 Introduction

The problem of transverse vibrations of beams subject to the passage of randomly moving loads is of considerable practical interest in the dynamics of structures, for it has a wide range of applications in civil, mechanical and aircraft industries. This problem has been especially studied in the context of the bridges behavior; for example, Fryba Fryba1 ; Fryba2 has investigated the response of a simply supported beam using an Euler-Bernoulli model subject to a single moving force, as well as the response to continuous random forces, in particular the effects of the constant speed and damping on the beam dynamics. Zibdeh Zibdeh has investigated the vibrations of a simply supported elastic beam under the action of a point load moving with random non constant velocity, while the beam is also subject to axial deterministic forces. Closed form solutions for the mean and variance of the response have also benn obtained. To take into account the rotary and high frequency motion of beam elements, Chang Chang proposed to treat the deterministic and random vibration analysis of a Rayleigh-Timoshenko beam on an elastic foundation with modal analysis to compute the dynamic responses of the structure (such as the displacement and bending moment) and some statistical responses (such as the mean square values of the dynamic displacement and the mean bending moment). With the same method Argento et al. Argento have studied the response of a rotating Rayleigh beam with different boundary conditions subject to an axially accelerating distributed surface line load.

The above mentioned works deal with bridge beams vibrations caused by a random stream of moving forces, assuming that the loads move with an exactly determined velocity. A step further is to consider the velocities as stochastic variables Sniady1 ; Sniady2 ; Sniady3 . If the beam is loaded by stochastically moving loads, the problem is more complicated and, generally, only numerical methods allow to retrieve the resulting vibrations Bryja

The aim of this work is to give an analytical approach to characterize the probabilistic features of the nonlinear beam response, namely the mean-square amplitude and the probability density function due to loads moving with stochastic velocities. The effects of the number of moving loads, the mean velocity and the intensity of the velocity disturbances on the dynamics of the beam are also analysed.

We believe that we have made progresses in three directions: (a) Modelling the dynamic response of a Rayleigh two way beam subjected to a train of forces moving with stochastic velocity in such a manner that is, at least in the lowest order, tractable; (b) Estimating the mean square amplitude of the beams by using the Fourier transform and the theory of residues method; (c) Exploiting the stochastic averaging method to retrieve the stationary probability density function and its dependence on the intensity of the stochastic velocity. The latter result is, in our view, the most striking effect we have found.

The paper is organized as follows. After this Introduction, the physical model and the algorithm for numerical simulations are presented in Section 2. In Section 3, the results of the Fourier analysis of the modal equation is described. Then, the results of the stochastic averaging method, the amplitude of the oscillations and the stationary probability density function of the amplitude for the noisy modal equations are retrieved in Section 4. Finally, Section 5 is devoted to concluding remarks.

2 The Beam Model

This Section is devoted to the presentation of the system (Subsect. 2.1), and to the corresponding reduced modal equations (Subsect. 2.2).

2.1 Mathematical modelling

Let us consider a beam of finite length $L$, that is a nonlinear elastic structure. In particular, our attention is paid to the geometric nonlinearities described due to the Euler-Bernoulli law, that states that the bending moment of the beam is proportional to the change in the curvature produced by the load fertis ; oumbe ; Duffing . To take into account the high frequency motion of the beam, a Rayleigh beam correction (up to the second order of the bending angle) seon is used to refine the theory of Euler-Bernoulli beam for high frequency motion. When the governing equation for the vertical displacement of the beam incorporates the Rayleigh term into the analysis, the correction affects the summing moments produced in the simple Euler-Bernoulli theory.

Vibrations of the beam are caused by a set of point forces of constant amplitudes P, the inter-arrival times are different, deterministic variables $t_{i}$ and the forces are moving along the beam with stochastic velocities ${\it v_{i}}$ (see Fig. 1).

Considering the classical viscous damping for the viscosity materials and Newton second law of motion, for an infinitesimal element of the beam, the equation of motion for the small deformations: $\theta\left(x,t\right)\approx\partial w\left(x,t\right)/\partial x$ (here $w(x,t)$ is the transverse defection of the beam at point $x$ and time $t$) is obtained as

[TABLE]

with the boundary and initial conditions

[TABLE]

Here $EI$ denotes the flexural rigidity of the beam, $M_{b}$ the beam mass per unit length, $C$ the damping coefficient, $R_{a}$ the transverse Rayleigh beam coefficient and $\delta[...]$ the Dirac delta function. The position of the $i^{th}$ force at the time $t$ reads ${x_{i}}\left({t-{t_{i}}}\right)$, where $t_{i}=(i-1)d/v_{0}$ is the deterministic arriving time of the $i^{th}$ load at the beam. The average spacing loads is $d$ , $v_{0}$ the average speed of the moving loads, and $N_{v}$ the number of the applied loads. Eq. (1) contains a nonlinear term due to a possible high deflection motion of the beam fertis ; oumbe ; Duffing . To facilitate a compact representation of the equations, a window function $\varepsilon_{i}$ is employed Nikkhoo : ${\varepsilon_{i}}=0$ when the load has left the beam and ${\varepsilon_{i}}=1$ while the load is crossing the beam.

A more realistic and useful model of highway traffic loads takes into account random arrival times, in the form, for examples, of Poisson process Tung1 or renewal counting process Tung2 , to represent the vehicular traffic fluctuations. Load fluctuations are thus assumed to be a set of point forces of constant amplitudes whose inter-arrival times are different, for the forces are moving with stochastic velocities.

Thus, let us consider that the velocities $v_{i}(t-{t_{i}})$ are Gaussian distributed around the average speed $v_{0}$ Sniady1 . These velocities are a function of the difference between the time $t$ and and the instant at which the vehicles enter the bridge, $t_{i}=(i-1)d/v_{0}$

[TABLE]

Here $v_{i}(t-{t_{i}})$ is the stochastic velocity of the $i^{th}$ force, $\sigma_{v}$ its standard deviation and $\xi_{i}(t-{t_{i}})$ the velocity disturbances which it is assumed to be independent and stationary white noise random processes; i.e.

[TABLE]

The brackets $\langle...\rangle$ denote the time average and $\gamma_{v}=v_{0}\sigma_{v}$.

In the presence of the stochastic term there is a finite probability that the fluctuations produce a negative velocity, that become sizeable when fluctuations $\sigma_{v}$ are comparable with the average speed $v_{0}$. We do not exclude such negative velocities, that we interpret as a vehicle moving in the opposite direction, thus assuming that the beam represents in fact a two-way bridge.

2.2 Modal equations

If one takes into account the boundary conditions, the transversal deflection $w(x,t)$ for the simply-supported beam can be represented in a series form as

[TABLE]

Here, $q_{n}\left(t\right)$ is the amplitude of the $n^{th}$ mode, and $\sin\left(n\pi x/L\right)$ is the solution of the eigen value problem, that depends on the boundary conditions of the free oscillations of the beam. It is convenient to adopt the following dimensionless variables:

[TABLE]

The equivalent stochastic dimensionless modal equation is obtained substituting Eq. (3) and Eq. (5) into Eq. (1) and considering the first mode only. We have limited the analysis to the first mode, inasmuch the first mode of vibrations is expected to carry most of the energy, and therefore one hopes that it could suffice to obtain a first estimate of the system behavior. If we indicate with $W_{i}(\tau-\tau_{i})$ a unit Wiener stochastic process:

[TABLE]

with

[TABLE]

The ansatz (5) of the modal analysis is but the simplest method to deal with this intrinsically nonlinear problem, for the nonlinear terms of Eq. (1) only results in the cubic nonlinearity of Eq. (7). However, the method is capable both to retain (to a limited and approximated extent) the essential nonlinearity of the full partial differential equation, and at the same time allows for an analytical treatment, as we shall see in the following. Also, the roughness of the method calls for detailed numerical simulations to validate the analysis. As we shall see in the following, numerical simulations support the usefulness of the method.

Eq. (7) amounts to a stochastic Duffing oscillator which describes the unbounded or catastrophic motion of the beam. The catastrophic behavior of the beam is related to the configuration of the potential of the system, as described in details in Ref. nana

3 Modal Analysis

It has been shown by Iwankiewicz and Śniady Iwankiewicz that the influence of the free vibration (${\varepsilon_{i}}=0$) on the probabilistic characteristics of the structure response is negligibly small when the speed of the moving forces is below $130$ km/h (about $36$ m/s). Therefore, for simplicity, it is assumed in the first case (${\varepsilon_{i}}=0$) that the dynamic response function is equal to zero ($\chi(\tau)=0$) and for the second case (${\varepsilon_{i}}=1$) the response function is calculated from the equation (considering the case of a single moving load):

[TABLE]

It is observed that the right hand side of Eq.(9) is a harmonic function with constant amplitude and random phases (mathematically equivalent to frequency fluctuations of a nonmonochromatic drive Filatrella02 ), and therefore it amounts to a bounded or sine-Wiener noise $\eta(\tau)$ Bobryk05 , whose covariance is given by

[TABLE]

The response of the stochastic Eq. (9) is obtained using the Fourier transform associated with the residues theorem Datta , as will be shown in the following.

3.1 Fourier Analysis

The equivalent linearization of Eq.(9) is

[TABLE]

where

[TABLE]

and $A$ is the root mean-square displacement. Using the Fourier Transformations one obtains

[TABLE]

(here $i^{2}=-1\,\,\,$ and $TF(.)$ is the Fourier transform operator). The solution of Eq. (11) in Fourier space reads

[TABLE]

From Eq. (10), one have

[TABLE]

where ${S_{\eta}}(\omega^{\prime})$ is the spectral density of the noise $\eta(\tau)$ defined by

[TABLE]

Hence,

[TABLE]

The mean-square amplitude can be calculated as

[TABLE]

It follows from the above definitions that:

[TABLE]

This is the main result of this part: the dependence of the oscillations amplitude upon the beam parameters (damping $\lambda$, natural frequency $\omega_{0}=1$ and nonlinear component $\beta$) and of the loads traffic (loads weights intensity $\Gamma$, average velocity $\Omega$, and velocity fluctuation intensity $\gamma$). In fact substituting Eq. (12a) into Eq. (19), the amplitude can be found from the roots of

[TABLE]

with

[TABLE]

3.2 Algorithm of stochastic numerical simulations

To numerically treat the random process of the nonlinear Eq. (9), the stochastic Fourth-order Runge Kutta (RK) algorithm developed by Kasdin Kasdin is used. Introducing the new variables y($\tau$) and z($\tau$), the system (9) can be rewritten as

[TABLE]

The RK solution for $\chi_{k}$, $y_{k}$ and $z_{k}$ are given by the following set of equations:

[TABLE]

where

[TABLE]

This numerical method can be extended to simulate the case of several moving loads. Here, $r_{n}$ (n=1,…,4) are the Gaussian white noise from random numbers $x_{m}$ (m=1,…,8) by using the Box -Mueller algorithm Knuth . The step $\Delta\tau$ here used is $\Delta\tau=0.001$ and the results are averaged over $100$ realizations. The coefficients (see Table 1) $\alpha_{i}$, $q_{i}$ and $a_{ji}$ are chosen in the deterministic case to ensure that $\chi_{k}$, $y_{k}$ and $z_{k}$ simulate the solution $\chi(\tau)$, $y(\tau)$ and $z(\tau)$ with error of order $\Delta{\tau^{5}}$. That is,

[TABLE]

These coefficients are given in Table 1.

3.3 Numerical Analysis

To check the validity of the analytical estimates, we have compared some analytical results with numerical simulations. The physical parameters of the beam used throughout the paper are: $L=34.0{m}$, $E=3.0\times{10^{10}}N/{m^{2}}$, $I=3.07{m^{4}}$, ${M_{b}}=11400.0{K_{g}}/m$, $C=350.5N.S/m$, $S=0.02{m^{2}}$, and the dimensionless parameters are: $\beta={\rm{-0}}{\rm{.52}}$, $\lambda={\rm{0}}{\rm{.00078}}$, $\Gamma={\rm{0}}{\rm{.0022}}$.

The mean response of the amplitude of the oscillators, as estimated with the modal equations method, Eq. (20), compared with numerical simulations of Eq. (9) obtained through the stochastic RK algorithm is displayed in Fig. 2(a). The Figure demonstrates that there is a resonant velocity $v\simeq 0.32$ where the amplitude of the oscillations increases, in both the theoretical prediction and the numerical simulations.

The analysis of a single load cannot be extended to the multiple vehicles case described by Eq.(7), inasmuch the superposition principle of a linear system is not established for the nonlinear dynamics system (7). As the present analysis cannot be rigorously extended to multiple load cases, the limits of validity of the extension can only be numerically verified, as shown in Fig. 2(b). In the Figure it is observed that the resonance is made more pronounced by the increase of the number of loads $N_{v}$. Another interesting feature of Fig. 2(b) is the reentrant behavior. For some values below the resonant velocity $v=0.32$ there are as many as three different amplitudes that correspond to the same speed. It is interesting to notice that the peculiar behavior is observed in both the analytical treatment for the noise value $\gamma=0.1$, see Fig.  2(a). At this level fluctuations are comparable with the resonant speed $v_{0}=0.32$ in Fig. 2(a). Physically, it corresponds to rare but sizeable negative velocities. Instead, in Fig. 2(b) fluctuations are so small that it practically represents the case of a one way bridge. We conclude that the analytical treatment captures the main effect, as the agreement is fairly good and the resonance is correctly captured by the analysis in a wide range of parameters.

The reentrant behavior is further analyzed in Fig. 3, where the mean-square amplitude as a function of the noise intensity for different values of the loads speed $\it v$ is displayed. In Fig. 3 the peak of the oscillations, marked by an oval, occurs at relatively high noise values ($\gamma\simeq 1.4$) for low speed ($v=0.1$) speeds, decreases at lower noise levels ($\gamma\simeq 0.4$) for higher speed ($v=0.35$), and moves again to higher fluctuations ($\gamma\simeq 0.6$) with a further increase of the speed ($v=0.4$). Thus, not only the amplitude has a maximum response for a finite value of the noise intensity, but also the peak position exhibits a nontrivial behavior. In general these features can be ascribed to the peculiar character of bounded noise in model Eq. (9) Li12 . We note that the amplitude decrease occurs for fluctuation levels that are comparable or higher than the average speed, and therefore they basically describe a two way bridge.

There are the values of $\gamma$ for which the amplitude is almost constant at the maximum; a further increase of the noise amplitude $\gamma$ leads to a decrease of the amplitude (see Fig. 3).

Beside the reentrant effect, it is noticed that the noise intensity tends to broad the resonance, and the number of loads increases both the mean square amplitude and the width of the peak. These features are naturally expected on physical grounds; more surprising is the effect of the single load weight on the mean square amplitude, as shown in Fig. 4. Increasing the load of each moving weight, one observes that the mean amplitude also shows a reentrant behavior with multiple solutions below the resonance speed (Fig. 4). It is significant to note that when the speed is below the resonance, one observes critical values of the weight load where we have a bifurcation of the system. On a general ground this could be again attributed to the presence of bounded noise, that can produce bimodal distributions Bobryk05 . More specifically, the effect can be related to the special features of the bounded noise spectral density (16), combined with the nonlinear amplitude-frequency relation of Eq.(20). These features give rise to the intricate effects above described, inasmuch the spectral content determined by Eq. (16) is, at variance with the uniform case of white noise, peaked around some frequencies. Thus, the frequency of the peaks depends upon the noise intensity: increasing the noise intensity the harmonic content is altered. This is the physical origin of the reentrant behavior observed in Fig. 3: when the noise intensity is increased it first moves the harmonic content towards the resonance of the beam, increasing the amplitude of the oscillations. The effect, however, saturates and eventually the term at the denominator $\gamma$ dominates.

In Figs. 2,3,4 the mean amplitude of the oscillations of the beam has been displayed. A more complete description in terms of the full distribution is presented in the following Section, where the probability distribution of the representative coordinate $\chi$ is analysed.

4 Stationary probability distribution of a catastrophic monostable system

In this Section, a Duffing oscillator with bounded Gaussian noise described by Eq. (9) is considered. In particular, the stationary probability distribution analytically with the method of stochastic averaging method Stratonovich ; Zhu ; Roberts ; Chamgoue is sought. To do so, it is assumed that the noise intensity is small and the following change of variables is introduced,

[TABLE]

Substituting Eqs. (26) into Eq. (9), we obtain

[TABLE]

The steady-state response for the case of perfect periodicity is considered first, $\gamma=0$, and Eqs. (27) become

[TABLE]

By applying the standard averaging method Nayfeh , Eqs. (28) reduce to

[TABLE]

The steady-state solutions of Eqs. (29) can be found by putting ${\it a=a_{0}},\varphi=\varphi_{0}$ and $\dot{a}=0$, $\dot{\varphi}=0$, this leads to the following result:

[TABLE]

This equation has more than one steady-state solution for some parameters. The variation of steady-state response ${\it a_{0}}$ as a function of the speed $v$ is compared with the numerical simulation of Eq. (9) and shown in Fig. 5. It can be seen from this figure that the deterministic response predicted by the standard averaging method is in good agreement with that obtained by the numerical simulations. In particular the resonance is correctly captured by the analysis. The time response of the system (9) and phase plot are shown in Fig. 6 for the noiseless case $\gamma=0$, and for $v=0.32$, $N_{v}=1$, $\lambda=0.0078$. Clearly, the response is periodic and the phase trajectory is a limit cycle.

Next, the stationary response of system (9) in the noisy case is determined, $\gamma\neq 0$. To do so, Eqs. (27) are rewritten as:

[TABLE]

To apply the method Stratonovich ; Zhu ; Roberts ; Chamgoue , one average at the frequency $\Omega$ and the following pair of stochastic equations for amplitude a($\tau$) and the phase $\varphi(\tau)$ is obtained:

[TABLE]

Here $W_{1}(\tau)$ and $W_{2}(\tau)$ are independent normalized Weiner processes. Clearly, Eq. (32) governing the amplitude $a(\tau)$ does not depend on $\varphi(\tau)$; thus, one can consider the probability density for $a(\tau)$, rather than a joint density for $a(\tau)$ and $\varphi(\tau)$. The probability density function $P(a(\tau),\tau|a(\tau_{0}),\tau_{0})$ for the amplitude is governed by the Fokker-Planck-Kolmogorov equation:

[TABLE]

In the stationary case, $P:{{\partial P\left({a,\tau}\right)}}/{{\partial\tau}}=0$, the solution of Eq. (33) is:

[TABLE]

where

[TABLE]

Combining Eq. (34) and Eqs. (35), we get

[TABLE]

where $N$ has been determined by the normalization condition:

[TABLE]

Moreover, one can find the peaks of the distribution Eq. (36) in the points where ${\partial P\left({a}\right)}/{\partial a}=0$:

[TABLE]

As the amplitude of the oscillations in the ansatz Eq. (26) is defined positive, it is concluded that, the distribution should exhibit just a peak. In Fig. 7 the stationary probability density as a function of the amplitude $a$ for different noisy velocity intensity $\gamma$ is displayed, compared with numerical solutions of Eq. (9). The numerical results confirm that the distribution has only one maximum. It is also evident that the analytic solution predicts a peak of the probability distribution at higher values than the numerical solution. The discrepancy, however, decreases as noise increases, and the agreement is fairly good when fluctuations are very large (compared to the average speed $v=0.25$). We therefore conclude that there is a better agreement for a two way bridge.

The results obtained in this section confirm that an increase of the noise intensity $\gamma$ leads to a decrease of the amplitude of the oscillations, Eq.(38). This is in qualitative agreement with the findings of the Fourier analysis of Section 3, as it confirms the counter-intuitive role of noise, that tends to decrease the amplitude of the oscillations. However, it is underlined that the stochastic averaging does not capture the amplitude of the oscillations at low noise values.

5 Concluding remarks

An analytic approach to the dynamics of a Rayleigh beam under the effect of loads moving with stochastic speed has been considered. To do so, one has first reduced the full partial differential equation to an effective nonlinear one-dimensional equation. For this reduced equation, the random terms have been modelled as equal weights moving with a stochastic velocity, thus neglecting other effects as the spread of the weights. Within these limitations, it has been investigated how the main parameters of the moving loads affect the beam response, and especially how the intensity of the random component of the loads velocity affects the dynamic behaviour of the beam. Applying the Fourier transform and the theory of residues a mean amplitude equation has been obtained. It has been found that velocity noise first causes an increase of the mean-square amplitude of the beam oscillations, and after a maximum a further increase of the noise causes a decrease of the mean-square displacement, a behavior that closely follows stochastic resonance phenomena Gammaitoni98 ; Wellens04 . In fact, high levels of noise appear to be paradoxically beneficial for the beam safety, inasmuch the amplitude only increases up to a special noise value where it exhibits a maximum damage to the beam wearing. Such counterintuitive effect of the random disturbance has been extensively debated in signal detection Ando01 and biological systems Mcdonnell11 , and are generally connected with suboptimal detection observable as the observed amplitude of the oscillations Addesso12 . However, at variance with electronic signal detection where stochastic resonance is extensively studied and can be possibly exploited, the beneficial effect of noise is of course of difficult application in engineering, for the addition of further noise sources is unpractical, if not just impossible. We notice, however, that the results indicate that negative velocities, corresponding to a two way bridge, might prove safer than a one way beam. From the physical point of view one ascribe the above behavior to the change in the harmonic content of the bounded noise that is altered, in a nontrivial manner, by an increase of the stochastic intensity.

The analysis has also allowed to estimate the response of the beam to the loads and the other system parameters. It has been thus demonstrated that the number of moving loads considerably influences the mean-square vibration amplitude of the beam and the appearance of multiple amplitude solutions for some range of the loads weight. Some numerical simulations of the modal equation, obtained with a stochastic version of the fourth-order Runge-Kutta algorithm, have confirmed the analytic results.

Finally, the stationary probability density function of the amplitude for the noisy modal equation has been analytically derived using a stochastic averaging method and compared to the numerical solution. It has been found that the agreement is fairly good, and that also this approach confirms that the most probable position to find the oscillator decreases when the intensity of the stochastic fluctuations of the velocity increase. However, stochastic averaging seems to be not able to capture the initial increase of the oscillations amplitude.

By way of short summary, we have tried to retrieve what is the effect of the load random velocities on the behavior of the mean amplitude versus the parameters. Some nonlinear dependences are new but not fully unexpected: the number of vehicles changes the resonant velocity (Figs. 2,5). More interesting is the unexpected influence of the load, that is not just nonlinear, but reentrant (Fig. 4). To single out what we regard as the most relevant finding, we focus on the non-monotonic effect of the noise intensity, Fig. 3, that discloses an unexpected beneficial role of noise intensity on the bridge stability.

We would like to conclude mentioning some open issues. First, it could prove interesting to understand the discrepancy between the stochastic averaging approach and the Fourier Transform approach. Second, the noise effect here considered, only in terms of velocity fluctuations, could be complemented with an analysis of other random sources, as the individual weights of the loads. Third, the deterministic underlying equation could possibly be treated with another approach to improve the deterministic approximation.

Acknowledgements.

G.F. acknowledges partial financial support from MIUR through PON R&C (Programma Operativo Nazionale ”Ricerca e Competitività”) 2007-2013 under Grant No. PON NAFASSY, PONa3_00007. Part of this work was completed during a research visit of Prof Nana Nbendjo at the University of Kassel in Germany. He is grateful to the Alexander von Humboldt Foundation for financial support within the Georg Forster Fellowship.