Vibrations of an Elastic Beam Subjected by Two Kinds of Moving Loads and Positioned on a Foundation having Fractional Order Viscoelastic Physical Properties

Abstract

The present chapter investigates both the effects of moving loads and of stochastic wind on the steady-state vibration of a first mode Rayleigh elastic beam. The beam is assumed to lay on foundations (bearings) that are characterized by fractional-order viscoelastic material. The viscoelastic property of the foundation is modeled using the constitutive equation of Kelvin-Voigt type, which contain fractional derivatives of real order. Based to the stochastic averaging method, an analytical explanation on the effects of the viscoelastic physical properties and number of the bearings, additive and parametric wind turbulence on the beam oscillations is provided. In particular, it is found that as the number of bearings increase, the resonant amplitude of the beam decreases and shifts towards larger frequency values. The results also indicate that as the order of the fractional derivative increases, the amplitude response decreases. We are also demonstrated that a moderate increase of the additive and parametric wind turbulence contributes to decrease the chance for the beam to reach the resonance. The remarkable agreement between the analytical and numerical results is also presented in this chapter.

1. Introduction

There is a large amount of vehicles passing through in-service bridges every day, while sizable wind blows on the bridge decks. Vibrations caused by the service loads is of great theoretical and practical significance in civil engineering. In this chapter, it follows a list, by no means exhaustive, of research related to this kind engineer problem. To start with, Xu et al. [1] have explored the basic dynamics interaction between suspension bridges and the combined effects of intense wind and a single moving train; however, the interaction between the wind and the train dynamics has been altogether neglected. In this limit the suspension bridge response is dominated by wind force. The coupled dynamic analysis of vehicle and cable-stayed bridge system under turbulent wind has also been recently conducted by Xu and Guo [2] under the other limit of low wind speed. In the same view, the both effects of turbulent wind and moving loads on the brigde response are investigated numerically by Chen and Wu [3]. Another interesting results related to the problem of the dynamics of bridges subjected to the combined dynamic loads of vehicles and wind are presented in Refs. [4, 5, 6]. To summarize: from the standpoint of bridge engineering the disturbances, either due to wind (in low or high speed limits), passage of heavy loads (single massive trains or disordered aggregates of smaller freight carriers, result in a complex interaction with the bridge vibrations. However, the elastic properties of the bridges are enhanced by the insertion of bearings (the part ranging between the bridge deck and the piers) as a possible protection against severe earthquakes. For if one wants the bearings to protect the bridge, they should isolate the structure from ground vibrations and/or transfer the load to the foundation [7]. Noticed that the bearings can be constituted by some elastic or viscoelastic material. In the literature, the dynamics analysis of bridges with elastic bearings to moving loads has received limited attention. nevertheless, some authors like Yang et al. [8], Zhu and Law [9], Naguleswaran [10] and Abu Hilal and Zibdeh [11] have adressed a very interesting resultats about this subjet. There are investigated the pros, and the cons, of the elastic bearings.

The bearings can also be constituted by some viscoelastic materials (such as elastomer)[12]. Therefore, The viscoelastic property of the materials may be modelled by using the constitutive equation of Kelvin-Voigt type, which contain fractional derivatives of real order. In this Chapter we aim to investigate first the pros, and the cons, of the viscoelastic bearings and second the turbulence effect of the wind actions on the response of beam. To accomplish our goal some methods (analytical [13, 14, 15] and numerical [13, 14, 15, 16]) are used.

2. Structural system model

2.1 Mathematical modelling

In this chapter, a simply supported Rayleigh beam [17, 18] of finite length L with geometric nonlinearities [19, 20], subjected by two kinds of moving loads (wind and train actions) and positioned on a foundation having fractional order viscoelastic physical properties is considered as structural system model and presented in Figure 1 .

Figure 1.

Sketch of a beam-foundation system subjected to wind actions and series of moving forces.

As demonstrated in (Appendix A) and in Refs. [16, 19, 21], the governing equation for small deformation of the beam-foundation system is given by:

$$\begin{array}{l}\mathit{\rho S}\frac{{\partial }^{2}w\left(x,t\right)}{\partial t^2}+\mathit{EI}\frac{{\partial }^{4}w\left(x,t\right)}{\partial x^4}-\frac{3}{2}\mathit{EI}\frac{{\partial }^2}{\partial x^2}\left[\frac{{\partial }^{2}w\left(x,t\right)}{\partial x^2}{\left(\frac{\partial w\left(x,t\right)}{\partial x}\right)}^2\right]-\mathit{\rho I}\frac{{\partial }^{4}w\left(x,t\right)}{\partial x^2\partial t^2}+\mu \frac{\partial w\left(x,t\right)}{\partial t}\\ -\frac{\mathit{ES}}{2L}\frac{{\partial }^{2}w\left(x,t\right)}{\partial x^2}\underset{0}{\overset{L}{\int }}{\left(\frac{\partial w\left(x,t\right)}{\partial x}\right)}^2\mathrm{d}x+\sum _{j=1}^{N_P}\left(k_j+c_j{D}_t^{{\alpha }_j}\right)w\left(x,t\right)\delta \left[x-\frac{\mathit{jL}}{N_P+1}\right]=F_{\mathit{ad}}\left(x,t\right)+\\ P\sum _{i=0}^{N-1}{\epsilon }_i\delta \left[x-x_i\left(t-t_i\right)\right].\end{array}$$ E1

In which $\mathit{\rho S}$ , $\mathit{EI}$ , $\mathit{\rho I}$ , $\mu$ , $w\left(x,t\right)$ , are the beam mass per unit length, the flexural rigidity of the beam, the transverse Rayleigh beam coefficient, the damping coefficient and the transverse displacement of the beam at point $x$ and time $t$ , respectively. In Eq. (1), $\mathit{\rho S}\frac{{\partial }^{2}w\left(x,t\right)}{\partial t^2}$ represents the inertia term of the beam per unit length, $\mathit{\rho I}\frac{{\partial }^{4}w\left(x,t\right)}{\partial x^2\partial t^2}$ is the rotary inertia force of the beam per unit length, $\mu \frac{\partial w\left(x,t\right)}{\partial t}$ is the damping force of the beam per unit length, ${\sum }_{j=1}^{N_P}\left(k_j+c_j{D}_t^{{\alpha }_j}\right)w\left(x,t\right)\delta \left[x-\frac{\mathit{jL}}{N_P+1}\right]$ is the foundation-beam interaction force (per unit length of the beam’s axis), $\mathit{EI}\frac{{\partial }^{4}w\left(x,t\right)}{\partial x^4}$ and $\frac{3}{2}\mathit{EI}\frac{{\partial }^2}{\partial x^2}\left[\frac{{\partial }^{2}w\left(x,t\right)}{\partial x^2}{\left(\frac{\partial w\left(x,t\right)}{\partial x}\right)}^2\right]$ are the linear and the nonlinear term of the rigidity of the beam essentially due to Euler Law. The nonlinear term is obtained by using the Taylor expansion of the exact formulation of the curvature up to the second order [19, 22, 23]. $\frac{\mathit{ES}}{2L}\frac{{\partial }^{2}w\left(x,t\right)}{\partial x^2}\underset{0}{\overset{L}{\int }}{\left(\frac{\partial w\left(x,t\right)}{\partial x}\right)}^2\mathrm{d}x$ is the inplane tension of the beam [19, 20]. The terms on the right-hand side of Eq. (1) are used to describe the wind an train actions over the beam. In particular, the first term $F_{\mathit{ad}}\left(x,t\right)$ is the aerodynamic force given after some derivations by [24, 25, 26]:

$$F_{\mathit{ad}}\left(x,t\right)=\frac{1}{2}{\rho }_a{\mathit{bU}}^2\left[A_0+\frac{A_1}{U}\frac{\partial w\left(x,t\right)}{\partial t}+\frac{A_2}{U^2}{\left(\frac{\partial w\left(x,t\right)}{\partial t}\right)}^2\right],$$ E2

where ${\rho }_a$ is the air mass density, $b$ is the beam width, $A_j$ (j = 0, 1, 2) are the aerodynamic coefficients ( $A_0$ = 0.0297, $A_1$ = 0.9298, $A_2$ = -0.2400) [24]. $U$ is the wind velocity which can be decomposed as $U=\overline{u}+u\left(t\right)$ , where $\overline{u}$ is a constant (average) part representing the steady component and $u\left(t\right)$ is a time varying part representing the turbulence. It is assumed in this work that ( $\overline{u}\gg u\left(t\right)$ ).

According to Figure 1 , the boundary conditions of the beam are considered as [27]

$$w\left(0,t\right)=w\left(L,t\right)=0;\frac{{\partial }^{2}w\left(0,t\right)}{\partial x^2}=\frac{{\partial }^{2}w\left(L,t\right)}{\partial x^2}=0.$$ E3

The next section deals with the reduction of the main equation Eq. (1).

2.2 Reduced model equation

According to the Galerkin’s method [27, 28] and by taking into account the boundary conditions of the beam, the solution of the partial differential Eq. (1) is given by:

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

where $q_n\left(t\right)$ are the amplitudes of vibration and $sin\left(\mathit{n\pi x}/L\right)$ are modal functions solutions of the beam linear natural equation with the associated boundary conditions. It is convenient to adopt the following dimensionless variables:

$${\chi }_n=\frac{q_n}{l_r},\tau ={\omega }_{0}t,\xi =\frac{u}{u_c},$$ E5

the single one-dimensional modal equation with ${\chi }_n=\chi \left(\tau \right)$ is given as:

$$\begin{array}{l}\ddot{\chi }\left(\tau \right)+\left(2\lambda +{\vartheta }_1\right)\dot{\chi }\left(\tau \right)+\chi \left(\tau \right)+{\mathit{\beta \chi }}^3\left(\tau \right)+\eta \sum _{j=1}^{N_p}\left(k_j+c_j{\omega }_0^{{\alpha }_j}{D}_{\tau }^{{\alpha }_j}\right)\chi \left(\tau \right){sin}^2\left(\frac{\mathit{j\pi }}{N_p+1}\right)\\ ={\vartheta }_2{\dot{\chi }}^2\left(\tau \right)+{\vartheta }_0+\left({\theta }_0+{\theta }_1\dot{\chi }\left(\tau \right)\right)\xi \left(\tau \right)+f_0\sum _{i=0}^{N-1}{\epsilon }_{i}sin\mathrm{\Omega }\left[\tau -i\frac{d{\omega }_0}{v}\right].\end{array}$$ E6

With

$$\begin{array}{l}\mathrm{\Omega }=\frac{\mathit{\pi v}}{L{\omega }_0},f_0=\frac{2{\mathit{PL}}^3}{l_r\mathit{EI}{\pi }^4},\eta =\frac{2L^3}{\mathit{EI}{\pi }^4},\beta =\frac{l_r^2}{4}\left[\frac{S}{I}-\frac{3}{2}{\left(\frac{\pi }{L}\right)}^2\right],\\ {\vartheta }_1=\frac{{\rho }_a{\mathit{bL}}^3{A}_1\overline{U}}{2{\pi }^2\sqrt{\mathit{EI\rho }\left[L^{2}S+I{\pi }^2\right]}},{\vartheta }_0=\frac{2{\rho }_a{\mathit{bA}}_0{\overline{U}}^2{L}^4}{{\mathit{EIl}}_r{\pi }^5},{\vartheta }_2=\frac{4{\rho }_a{\mathit{bL}}^2{A}_2{l}_r}{3\rho \left[L^{2}S+I{\pi }^2\right]},\\ {\theta }_0=\frac{2\overline{U}{\rho }_a{\mathit{bL}}^4{U}_c{A}_0}{{\mathit{EIl}}_r{\pi }^6},{\theta }_1=\frac{{\rho }_a{\mathit{bL}}^3{A}_1{U}_c}{2{\pi }^2\sqrt{\mathit{EI\rho }\left[L^{2}S+I{\pi }^2\right]}},\lambda =\frac{\mu L^3}{2{\pi }^2\sqrt{\mathit{EI\rho }\left[L^{2}S+I{\pi }^2\right]}}.\end{array}$$ E7

And

$${\omega }_0=\frac{{\pi }^2}{L}\sqrt{\frac{\mathit{EI}}{\rho \left(L^{2}S+I{\pi }^2\right)}},l_r=\frac{L}{2}.$$ E8

According to Refs.[29, 30], Eq. (6) becomes:

$$\begin{array}{l}\ddot{\chi }\left(\tau \right)+\left(2\lambda +{\vartheta }_1\right)\dot{\chi }\left(\tau \right)+\chi \left(\tau \right)+{\mathit{\beta \chi }}^3\left(\tau \right)+\eta \sum _{j=1}^{N_p}\left(k_j+c_j{\omega }_0^{{\alpha }_j}{D}_{\tau }^{{\alpha }_j}\right)\chi \left(\tau \right){sin}^2\left(\frac{\mathit{j\pi }}{N_p+1}\right)\\ ={\vartheta }_2{\dot{\chi }}^2\left(\tau \right)+{\vartheta }_0+\left({\theta }_0+{\theta }_1\dot{\chi }\left(\tau \right)\right)\xi \left(\tau \right)+F_{0N}sin\left(\mathrm{\Omega }\tau \right)-G_{0N}cos\left(\mathrm{\Omega }\tau \right).\end{array}$$ E9

Where

$$\begin{array}{l}{F}_{0N}=P_0\left[1+\frac{2sin{\tilde{\tau }}_{0}sin\left(\left(N-1\right){\tilde{\tau }}_0\right)}{1-cos\left(2{\tilde{\tau }}_0\right)}cos\left(N{\tilde{\tau }}_0\right)\right],{\tilde{\tau }}_0=\frac{\mathit{d\pi }}{2L},\\ \\ G_{0N}=\frac{2P_{0}sin{\tilde{\tau }}_{0}sin\left(\left(N-1\right){\tilde{\tau }}_0\right)}{1-cos\left(2{\tilde{\tau }}_0\right)}sin\left(N{\tilde{\tau }}_0\right).\end{array}$$ E10

3. Analytical explanation of the model

3.1 Effective analytical solution of the problem

In order to directly evaluate the response of the beam, the stochastic averaging method [13, 14, 15] is first applied to Eq. (6), then the following change in variables is introduced:

$$\chi \left(\tau \right)=a_0+a\left(\tau \right)cos\psi ,\dot{\chi }\left(\tau \right)=-\mathrm{\Omega }a\left(\tau \right)sin\psi ,\psi =\mathrm{\Omega }\tau +\varphi \left(\tau \right),$$ E11

Substituting Eq. (11) into Eq. (9) we obtain:

$$\left\{\begin{array}{l}\dot{a}cos\psi -a\dot{\phi }sin\psi =0\\ \dot{a}sin\psi -a\dot{\phi }cos\psi =-\frac{1}{\mathrm{\Omega }}\left[M_1\left(a,\psi \right)+M_2\left(a,\psi \right)\right].\end{array}\right.$$ E12

where

$$\begin{array}{l}{M}_1\left(a,\psi \right)=F_{0N}sin\left(\psi -\phi \right)-G_{0N}cos\left(\psi -\phi \right)+\left(2\lambda +{\vartheta }_1\right)a\mathrm{\Omega }sin\mathit{\psi a}-\frac{1}{4}\beta a^{3}cos3\psi \\ -\left[1+\eta \sum _{j=1}^{N_p}{k}_j{sin}^2\left(\frac{\mathit{j\pi }}{N_p+1}\right)+3\beta a_0^2+\frac{3}{4}\beta a^2-{\mathrm{\Omega }}^2\right]cos\psi +\frac{a^2}{2}\left[{\vartheta }_2{\mathrm{\Omega }}^2-3\beta a_0\right]cos2\psi ,\\ \\ M_2\left(a,\psi \right)=-\eta \sum _{j=1}^{N_p}{c}_j{\omega }_0^{{\alpha }_j}{sin}^2\left(\frac{\mathit{j\pi }}{N_p+1}\right)D_{\tau }^{{\alpha }_j}\left(acos\psi \right)+\left({\theta }_0-{\theta }_1\mathrm{\Omega }asin\psi \right)\xi \left(\tau \right).\end{array}$$ E13

According to Eq. (13) The derivatives of the generalized amplitude $a$ and phase $\varphi$ could be solved as:

$$\left\{\begin{array}{l}\dot{a}=-\frac{1}{\mathrm{\Omega }}\left[M_1\left(a,\psi \right)+M_2\left(a,\psi \right)\right]sin\psi \\ \\ a\dot{\varphi }=-\frac{1}{\mathrm{\Omega }}\left[M_1\left(a,\psi \right)+M_2\left(a,\psi \right)\right]cos\psi \end{array}\right.$$ E14

$a_0$ satisfies the following non-linear equation:

$$\beta a_0^3+\left[1+\frac{3}{2}\beta a^2+\eta \sum _{j=1}^{N_p}{k}_j{sin}^2\left(\frac{\mathit{j\pi }}{N_p+1}\right)\right]a_0={\vartheta }_0-\frac{1}{2}{\vartheta }_2{\mathrm{\Omega }}^2{a}^2.$$ E15

Then, one could apply the stochastic averaging method [13, 14, 15] to Eq. (15) in time interval [0, T].

$$\left\{\begin{array}{l}\dot{a}=-\underset{T\to \infty }{lim}\frac{1}{T}\underset{0}{\overset{T}{\int }}\frac{1}{\mathrm{\Omega }}\left[M_1\left(a,\psi \right)+M_2\left(a,\psi \right)\right]sin\mathit{\psi d\psi }\\ a\dot{\phi }=-\underset{T\to \infty }{lim}\frac{1}{T}\underset{0}{\overset{T}{\int }}\frac{1}{\mathrm{\Omega }}\left[M_1\left(a,\psi \right)+M_2\left(a,\psi \right)\right]cos\mathit{\psi d\psi }\end{array}\right.$$ E16

According to this method, one could select the time terminal T as $T=2\pi /\mathrm{\Omega }$ in the case of periodic function ( $M_1\left(a,\psi \right)$ ), or $T=\infty$ in the case of aperiodic one ( $M_2\left(a,\psi \right)$ ). Accordingly, one could obtain the following pair of first order differential equations for the amplitude a( $\tau$ ) and the phase $\varphi \left(\tau \right)$ :

$$\begin{array}{l}\dot{a}=-\left(2\lambda -{\vartheta }_1\right)\frac{a}{2}-\frac{1}{2}\mathit{\eta a}\sum _{j=1}^{N_p}{c}_j{\omega }_0^{{\alpha }_j}{\mathrm{\Omega }}^{{\alpha }_j-1}{sin}^2\left(\frac{\mathit{j\pi }}{N_p+1}\right)sin\left(\frac{{\alpha }_j\pi }{2}\right)+\frac{1}{2\mathrm{\Omega }}\left[G_{0N}sin\phi -F_{0N}cos\phi \right]\\ +\frac{{\mathit{\pi \theta }}_1^{2}a}{8}\left[3S_{\xi }\left(2\mathrm{\Omega }\right)+2S_{\xi }\left(0\right)\right]+\frac{{\mathit{\pi \theta }}_0^2}{2{\mathrm{\Omega }}^{2}a}{S}_{\xi }\left(\mathrm{\Omega }\right)+\sqrt{\frac{{\mathit{\pi \theta }}_0^2}{{\mathrm{\Omega }}^2}{S}_{\xi }\left(\mathrm{\Omega }\right)+\frac{{\mathit{\pi \theta }}_1^2{a}^2}{4}\left[S_{\xi }\left(2\mathrm{\Omega }\right)+2S_{\xi }\left(0\right)\right]}{\xi }_1\left(\tau \right),\end{array}$$ E17

and

$$\begin{array}{l}a\dot{\phi }=\frac{a}{2\mathrm{\Omega }}\left[1+3\beta a_0^2-{\mathrm{\Omega }}^2+\frac{3}{4}\beta a^2+\eta \sum _{j=1}^{N_p}\left(k_j+c_j{\omega }_0^{{\alpha }_j}{\mathrm{\Omega }}^{{\alpha }_j}cos\left(\frac{{\alpha }_j\pi }{2}\right)\right){sin}^2\left(\frac{\mathit{j\pi }}{N_p+1}\right)\right]\\ +\frac{1}{2\mathrm{\Omega }}\left[G_{0N}cos\phi +F_{0N}sin\phi \right]-\frac{{\mathit{\pi \theta }}_1^{2}a}{4}{\mathrm{\Psi }}_{\xi }\left(2\mathrm{\Omega }\right)+\sqrt{\frac{{\mathit{\pi \theta }}_0^2}{{\mathrm{\Omega }}^2}{S}_{\xi }\left(\mathrm{\Omega }\right)+\frac{{\mathit{\pi \theta }}_1^2{a}^2}{4}{S}_{\xi }\left(2\mathrm{\Omega }\right)}{\xi }_2\left(\tau \right).\end{array}$$ E18

Here $S_{\xi }\left(\mathrm{\Omega }\right)$ and ${\mathrm{\Psi }}_{\xi }\left(\mathrm{\Omega }\right)$ are the cosine and sine power spectral density function, respectively [31]:

$$\begin{array}{l}{S}_{\xi }\left(\Omega \right)=\underset{-\infty }{\overset{+\infty }{\int }}R\left(\zeta \right)cos\left(\mathit{\Omega \tau }\right)\mathrm{d}\zeta =2\underset{0}{\overset{+\infty }{\int }}R\left(\zeta \right)cos\left(\mathit{\Omega \tau }\right)\mathrm{d}\zeta =2\underset{-\infty }{\overset{0}{\int }}R\left(\zeta \right)cos\left(\mathit{\Omega \tau }\right)\mathrm{d}\zeta ,\\ {\mathrm{\Psi }}_{\xi }\left(\Omega \right)=2\underset{0}{\overset{+\infty }{\int }}R\left(\zeta \right)sin\left(\mathit{\Omega \tau }\right)\mathrm{d}\zeta =-2\underset{-\infty }{\overset{0}{\int }}R\left(\zeta \right)sin\left(\mathit{\Omega \tau }\right)\mathrm{d}\zeta ,\\ \underset{-\infty }{\overset{+\infty }{\int }}R\left(\zeta \right)sin\left(\mathit{\Omega \tau }\right)\mathrm{d}\zeta =0;R\left(\zeta \right)=\mathrm{E}\left[\xi \left(\tau \right)\xi \left(\tau +\zeta \right)\right].\end{array}$$ E19

In this work, $\xi \left(\tau \right)$ is assumed to be an harmonic function with constant amplitude ${\sigma }_i$ and random phases ${\gamma }_i{B}_i\left(\tau \right)+{\theta }_i$ . So, according to Refs. [31, 32, 33] the following model of $\xi \left(\tau \right)$ has been chosen:

$$\xi \left(\tau \right)=\sum _{i=1}^m{\sigma }_{i}cos\left[{\omega }_i\tau +{\gamma }_i{B}_i\left(\tau \right)+{\theta }_i\right],$$ E20

this model of the turbulent component of the wind $\xi \left(\tau \right)$ amounts to a bounded or cosine-Wiener noise, whose spectral density is given by:

$${\Phi }_{\xi }\left(\omega \right)=\sum _{i=1}^m\frac{{\sigma }_i^2{\gamma }_i^2\left({\omega }^2+{\omega }_i^2+{\gamma }_i^4/4\right)}{4\pi \left[{\left({\omega }^2-{\omega }_i^2-{\gamma }_i^4/4\right)}^2+{\gamma }_i^2{\omega }^2\right]}.$$ E21

The next sections of this chapter will presents the analytical developments that we have made in order to express the beam response as a function of the system parameters. Then, let’s start with the case where the beam is subjected to the moving loads only.

3.2 Analytical estimate of the beam response under moving loads only

We first consider system (1) with only deterministic moving loads ( $F_{\mathit{ad}}\left(x,t\right)=0$ ) neglecting wind effects on the beam. If ${\vartheta }_1={\theta }_0={\theta }_1=0$ , Eqs. (15) and (16) become:

$$\dot{a}=-\mathit{\lambda a}-\frac{1}{2}\mathit{\eta a}\sum _{j=1}^{N_p}{c}_j{\omega }_0^{{\alpha }_j}{\Omega }^{{\alpha }_j-1}{sin}^2\left(\frac{\mathit{j\pi }}{N_p+1}\right)sin\left(\frac{{\mathit{\pi \alpha }}_j}{2}\right)+\frac{1}{2\mathrm{\Omega }}\left[G_{0N}sin\varphi -F_{0N}cos\varphi \right],$$ E22

and

$$\begin{array}{l}a\dot{\phi }=\frac{a}{2\mathrm{\Omega }}\left[1-{\mathrm{\Omega }}^2+3\beta a_0^2+\frac{3}{4}\beta a^2+\eta \sum _{j=1}^{N_p}\left(k_j+c_j{\omega }_0^{{\alpha }_j}{\mathrm{\Omega }}^{{\alpha }_j}cos\left(\frac{{\mathit{\pi \alpha }}_j}{2}\right)\right){sin}^2\left(\frac{\mathit{j\pi }}{N_p+1}\right)\right]\\ +\frac{1}{2\mathrm{\Omega }}\left[G_{0N}cos\phi +F_{0N}sin\phi \right].\end{array}$$ E23

By substituting $a=A,\varphi =\mathrm{\Phi }$ and $\dot{a}=0$ , $\dot{\varphi }=0$ in Eqs. (22) and (23), algebraic manipulations give for the steady-state vibrations of the system response $A$ the following non-linear equation:

$$\frac{9}{16}{\beta }^2{A}^6-\frac{3}{2}\beta {\mathrm{\Theta }}_1\left({\alpha }_j\right)A^4+\left[{\mathrm{\Theta }}_1^2\left({\alpha }_j\right)+{\mathrm{\Theta }}_2^2\left({\alpha }_j\right)\right]A^2=F_{0N}^2+G_{0N}^2,$$ E24

with

$$\begin{array}{l}{\mathrm{\Theta }}_1\left({\alpha }_j\right)={\mathrm{\Omega }}^2-1-3\beta a_0^2-\eta \sum _{j=1}^{N_p}\left(k_j+c_j{\omega }_0^{{\alpha }_j}{\mathrm{\Omega }}^{{\alpha }_j}cos\left(\frac{{\mathit{\pi \alpha }}_j}{2}\right)\right){sin}^2\left(\frac{\mathit{j\pi }}{N_p+1}\right),\\ {\mathrm{\Theta }}_2\left({\alpha }_j\right)=2\mathrm{\Omega }\lambda +\eta \sum _{j=1}^{N_p}{c}_j{\omega }_0^{{\alpha }_j}{\mathrm{\Omega }}^{{\alpha }_j}sin\left(\frac{{\mathit{\pi \alpha }}_j}{2}\right){sin}^2\left(\frac{\mathit{j\pi }}{N_p+1}\right).\end{array}$$ E25

The stability of the steady-state vibration of the system response is investigated by using the method of Andronov and Witt [34] associated to the Routh-Hurwitz criterion [35]. Thus, the steady-state response is asymptotically stable if Eq. (26) is satisfied and unstable if Eq. (27) is satisfied:

$${\left(\frac{{\Theta }_2\left({\alpha }_j\right)}{2\mathrm{\Omega }}\right)}^2+\frac{1}{4\mathrm{\Omega }}\left[\frac{3\beta }{4}{A}^2-{\Theta }_1\left({\alpha }_j\right)\right]\times \left[\frac{9\beta }{4}{A}^2-{\Theta }_1\left({\alpha }_j\right)\right]>0,$$ E26

$${\left(\frac{{\Theta }_2\left({\alpha }_j\right)}{2\mathrm{\Omega }}\right)}^2+\frac{1}{4\mathrm{\Omega }}\left[\frac{3\beta }{4}{A}^2-{\Theta }_1\left({\alpha }_j\right)\right]\times \left[\frac{9\beta }{4}{A}^2-{\Theta }_1\left({\alpha }_j\right)\right]<0.$$ E27

The trivial solution of Eq. (15) is $a_0=0$ .

What about the case where the beam is subjected to the stochastic wind loads?

3.3 Approximate solution of the beam response subjected to wind loads only

In this case ( $F_{\mathit{ad}}\left(x,t\right)\ne 0$ ) and $F_{0N}=G_{0N}=0$ , Eqs. (22) and (23) become:

$$\begin{array}{l}\mathit{da}=\left\{-\left(2\lambda -{\vartheta }_1\right)\frac{a}{2}-\frac{1}{2}\mathit{\eta a}\sum _{j=1}^{N_p}{c}_j{\omega }_0^{{\alpha }_j}{\mathrm{\Omega }}^{{\alpha }_j-1}{sin}^2\left(\frac{\mathit{j\pi }}{N_p+1}\right)sin\left(\frac{{\alpha }_j\pi }{2}\right)+\frac{{\mathit{\pi \theta }}_1^{2}a}{8}\left[3S_{\xi }\left(2\mathrm{\Omega }\right)+2S_{\xi }\left(0\right)\right]\right\}\mathit{d\tau }\\ +\frac{{\mathit{\pi \theta }}_0^2}{2{\mathrm{\Omega }}^{2}a}{S}_{\xi }\left(\mathrm{\Omega }\right)\mathit{d\tau }+\sqrt{\frac{{\mathit{\pi \theta }}_0^2}{{\mathrm{\Omega }}^2}{S}_{\xi }\left(\mathrm{\Omega }\right)+\frac{{\mathit{\pi \theta }}_1^2{a}^2}{4}\left[S_{\xi }\left(2\mathrm{\Omega }\right)+2S_{\xi }\left(0\right)\right]}{\mathit{dW}}_1\left(\tau \right),\end{array}$$ E28

and

$$\begin{array}{l}\mathit{d\phi }=\frac{1}{2\mathrm{\Omega }}\left\{1+3\beta a_0^2-{\mathrm{\Omega }}^2+\frac{3}{4}\beta a^2+\eta \sum _{j=1}^{N_p}\left(k_j+c_j{\omega }_0^{{\alpha }_j}{\mathrm{\Omega }}^{{\alpha }_j}cos\left(\frac{{\alpha }_j\pi }{2}\right)\right){sin}^2\left(\frac{\mathit{j\pi }}{N_p+1}\right)\right\}\mathit{d\tau }\\ -\frac{{\mathit{\pi \theta }}_1^2}{4}{\mathrm{\Psi }}_{\xi }\left(2\mathrm{\Omega }\right)\mathit{d\tau }+\sqrt{\frac{{\mathit{\pi \theta }}_0^2}{{\mathrm{\Omega }}^2{a}^2}{S}_{\xi }\left(\mathrm{\Omega }\right)+\frac{{\mathit{\pi \theta }}_1^2{a}^2}{4}{S}_{\xi }\left(2\mathrm{\Omega }\right)}{\mathit{dW}}_2\left(\tau \right).\end{array}$$ E29

Here $W_1\left(\tau \right)$ and $W_2\left(\tau \right)$ are independent normalized Weiner processes. In order to evaluate the effects of wind parameters on the system response, we derive an evolution equation for the Probability Density Function (PDF) of the variable amplitude $a\left(\tau \right)$ . The Fokker-Planck equation corresponding to the Langevin (Eq. (28)) reads:

$$\begin{array}{c}\frac{\partial P\left(a,\tau \right)}{\partial \tau }=-\frac{\partial }{\partial a}\\ \left[\left(-\left(2\lambda -{\vartheta }_1\right)\frac{a}{2}-\frac{1}{2}\mathit{\eta a}\sum _{j=1}^{N_p}{c}_j{\omega }_0^{{\alpha }_j}{\Omega }^{{\alpha }_j-1}{sin}^2\left(\frac{\mathit{j\pi }}{N_p+1}\right)sin\left(\frac{{\alpha }_j\pi }{2}\right)+\frac{{\mathit{\pi \theta }}_0^2}{2{\Omega }^{2}a}{S}_{\xi }\left(\Omega \right)\right)P\left(a,\tau \right)\right]\\ -\frac{\partial }{\partial a}\left[\frac{{\mathit{\pi \theta }}_1^{2}a}{8}\left[3S_{\xi }\left(2\Omega \right)+2S_{\xi }\left(0\right)\right]P\left(a,\tau \right)\right]+\frac{1}{2}\left(\frac{{\mathit{\pi \theta }}_0^2}{{\Omega }^2}{S}_{\xi }\left(\Omega \right)+\frac{{\mathit{\pi \theta }}_1^2{a}^2}{4}\left[S_{\xi }\left(2\Omega \right)+2S_{\xi }\left(0\right)\right]\right)\frac{{\partial }^{2}P\left(a,\tau \right)}{\partial a^2}.\end{array}$$ E30

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

$$P_s\left(a\right)=\mathit{Na}{\left({\Gamma }_0+a^2{\Gamma }_1\right)}^{-\left(Q+1\right)},$$ E31

where

$$\begin{array}{ll}{\Gamma }_0& =\frac{{\mathit{\pi \theta }}_0^2}{{\Omega }^2}{S}_{\xi }\left(\Omega \right),{\Gamma }_1=\frac{{\mathit{\pi \theta }}_1^2}{4}\left[S_{\xi }\left(2\Omega \right)+2S_{\xi }\left(0\right)\right],Q=\frac{{\Gamma }_1-2{\Gamma }_2}{2{\Gamma }_1},\\ {\Gamma }_2& =-\frac{1}{2}\left(2\lambda -{\vartheta }_1\right)-\frac{1}{2}\eta \sum _{j=1}^{N_p}{c}_j{\omega }_0^{{\alpha }_j}{\Omega }^{{\alpha }_j-1}{sin}^2\left(\frac{\mathit{j\pi }}{N_p+1}\right)sin\left(\frac{{\mathit{\pi \alpha }}_j}{2}\right)+\frac{{\mathit{\pi \theta }}_1^2}{8}\left[3S_{\xi }\left(2\Omega \right)+2S_{\xi }\left(0\right)\right].\end{array}$$ E32

Above $N$ is a normalization constant that guarantees $\underset{0}{\overset{\infty }{\int }}{P}_s\left(a\right)\mathrm{d}a=1$ .

What about the case where the beam is subjected to the both moving vehicles and stochastic wind loads?

3.4 Approximate solution of the beam responses subjected to the both moving loads

Finally, the case where the beam is subjected to the series of lumped loads and the wind actions is investigated. For the analytical purposes, we assume that the beam is linear and it is submitted to only the additive effects of the wind loads. Thus, Eqs. (22) and (23) become:

$$\begin{array}{l}a=\left[-\left(2\lambda -{\vartheta }_1\right)\frac{a}{2}-\frac{1}{2}\mathit{\eta a}\sum _{j=1}^{N_p}{c}_j{\omega }_0^{{\alpha }_j}{\mathrm{\Omega }}^{{\alpha }_j-1}{sin}^2\left(\frac{\mathit{j\pi }}{N_p+1}\right)sin\left(\frac{{\alpha }_j\pi }{2}\right)+\frac{{\mathrm{\Gamma }}_0}{a}\right]\mathit{d\tau }\\ +\frac{1}{2\mathrm{\Omega }}\left[G_{0N}sin\phi -F_{0N}cos\phi \right]\mathit{d\tau }+\sqrt{{\mathrm{\Gamma }}_0}{\mathit{dW}}_1\left(\tau \right),\end{array}$$ E33

and

$$\begin{array}{l}\mathit{d\phi }=\frac{1}{2\mathrm{\Omega }}\left[1-{\mathrm{\Omega }}^2+\eta \sum _{j=1}^{N_p}\left(k_j+c_j{\omega }_0^{{\alpha }_j}{\mathrm{\Omega }}^{{\alpha }_j}cos\left(\frac{{\alpha }_j\pi }{2}\right)\right){sin}^2\left(\frac{\mathit{j\pi }}{N_p+1}\right)\right]\mathit{d\tau }\\ +\frac{1}{2\mathrm{\Omega }a}\left[G_{0N}cos\phi +F_{0N}sin\phi \right]\mathit{d\tau }+\frac{1}{a}\sqrt{{\mathrm{\Gamma }}_0}{\mathit{dW}}_2\left(\tau \right).\end{array}$$ E34

The averaged Fokker-Planck-Kolmogorov equation associated with the previous Itô Eqs. (33) and (34) is

$$\frac{\partial P\left(a,\varphi ,\tau \right)}{\partial \tau }=-\frac{\partial }{\partial a}\left({\overline{a}}_{1}P\left(a,\varphi ,\tau \right)\right)-\frac{\partial }{\partial \varphi }\left({\overline{a}}_{2}P\left(a,\varphi ,\tau \right)\right)+\frac{1}{2}\frac{{\partial }^2}{\partial a^2}\left({\overline{b}}_{11}P\left(a,\varphi ,\tau \right)\right)+\frac{1}{2}\frac{{\partial }^2}{\partial {\varphi }^2}\left({\overline{b}}_{22}P\left(a,\varphi ,\tau \right)\right),$$ E35

where

$$\begin{array}{l}{\overline{a}}_1=-\left(2\lambda -{\vartheta }_1\right)\frac{a}{2}-\frac{1}{2}\mathit{\eta a}\sum _{j=1}^{N_p}{c}_j{\omega }_0^{{\alpha }_j}{\Omega }^{{\alpha }_j-1}{sin}^2\left(\frac{\mathit{j\pi }}{N_p+1}\right)sin\left(\frac{{\alpha }_j\pi }{2}\right)+\frac{1}{2\Omega }\left[G_{0N}sin\varphi -F_{0N}cos\varphi \right]+\frac{{\mathrm{\Gamma }}_0}{a}\\ {\overline{a}}_2=\frac{1}{2\mathrm{\Omega }}\left[1-{\Omega }^2+\eta \sum _{j=1}^{N_p}\left(k_j+c_j{\omega }_0^{{\alpha }_j}{\Omega }^{{\alpha }_j}cos\left(\frac{{\alpha }_j\pi }{2}\right)\right){sin}^2\left(\frac{\mathit{j\pi }}{N_p+1}\right)+\frac{1}{a}\left[G_{0N}cos\varphi +F_{0N}sin\varphi \right]\right]\\ {\overline{b}}_{11}={\Gamma }_0\\ {\overline{b}}_{22}=\frac{{\Gamma }_0}{a^2}.\end{array}$$ E36

Applying the solution procedure proposed by Huang et al. [36], one obtains the following exact stationary solution

$$P_s\left(a,\varphi \right)=N^{\prime }aexp\left\{\frac{{\Gamma }_2^{\prime }}{{\Gamma }_0}{a}^2-\frac{a}{\mathrm{\Omega }\left({\mathrm{\Gamma }}_0^2+d_0^2\right)}\left[\left(d_0{G}_{0N}+F_{0N}{\Gamma }_0\right)cos\varphi +\left(d_0{F}_{0N}-G_{0N}{\Gamma }_0\right)sin\varphi \right]\right\}$$ E37

where $N^{\prime }$ is a normalization constant and

$$\begin{array}{l}{\Gamma }_2^{\prime }=-\frac{1}{2}\left(2\lambda -{\vartheta }_1\right)-\frac{\eta }{2}\sum _{j=1}^{N_p}{c}_j{\omega }_0^{{\alpha }_j}{\Omega }^{{\alpha }_j-1}{sin}^2\left(\frac{\mathit{j\pi }}{N_p+1}\right)sin\left(\frac{{\alpha }_j\pi }{2}\right),d_0=-\frac{{\Gamma }_0{\Gamma }_3}{{\Gamma }_2^{\prime }},\\ {\Gamma }_3=\frac{1}{2\Omega }\left[1-{\Omega }^2+\eta \sum _{j=1}^{N_p}\left(k_j+c_j{\omega }_0^{{\alpha }_j}{\Omega }^{{\alpha }_j}cos\left(\frac{{\alpha }_j\pi }{2}\right)\right){sin}^2\left(\frac{\mathit{j\pi }}{N_p+1}\right)\right].\\ \end{array}$$ E38

4. Numerical analysis of the model

All the parameters concerning the chosen models of beam, of foundation and of the aerodynamic force are presented in Ref. [21]. Theses parameters clearly help to calculated the dimensionless coefficients defined in Eq. (7). It is well know that the validation of the results obtained through analytical investigation is guarantee by the perfect match with the results obtained through numerical simulations. thus, The numerical scheme used in this chapter is based on the Grunwald-Letnikov definition of the fractional order derivative Eq. (39). [37, 38, 39, 40] and the Newton-Leipnik algorithm [37, 38]

$$D_{\tau }^{\alpha }\left[\chi \left({\tau }_{n_f}\right)\right]\approx h^{-\alpha }\sum _{l=0}^{n_f}{C}_l^{\alpha }\chi \left({\tau }_{n_f-l}\right),$$ E39

where $h$ is the integration step and the coefficients $C_l^{\alpha }$ satisfy the following recursive relations:

$$C_0^{\alpha }=1,C_l^{\alpha }=\left(1-\frac{1+\alpha }{l}\right)C_{l-1}^{\alpha }.$$ E40

Now we display in some figures the effects of the main parameters of the proposed model. For example, Figure 2(a) shows the effect of the number of the bearings on the amplitude of vibration of the beam. This graph also shows a comparison between the results from the mathematical analysis (curve with dotted line) and the results obtained from numerical simulation of Eq. (9) using Eq. (39). The match between the results shows a good level of precision of the approximation made in obtaining Eq. (24). This figure also reveals that the vibration amplitude of the beam decreases and the resonance frequency of the system increases as the number of bearings increases. Figure 2(b) , shows the effect of loads number on the beam response. It is observed that as the value of $N_v$ increases, the amplitude of vibration at the resonant state merely increases.

Figure 2.

Effects of the number of the bearings $N_p$ (a) and number of moving loads $N_v$ (b) on the amplitude response of the beam when driving frequency $\Omega$ : Analytical curves (-), numerical curves (o). All The parameters are given in [21].

Looking at the effects of the order of the fractional derivative $\alpha$ on the amplitude of the beam, we obtain the graph of Figure 3 . Small value of $\alpha$ leads to large value of the maximum vibration amplitude. It is also clearly shown that the system is more stable for the highest order of the derivative. The multivalued solution appears for small order and disappears progressively as the order increases. The resonance (a peak of the amplitude) appears as the parameter $k_0^{\prime }$ increases, see Figure 3(a)-(d) . The good match between the analytical and the numerical results gives a validation of the approximations made.

Figure 3.

Effect of the fractional order on the amplitude of the beam A when the dimensionless stiffness $k_o^{\prime }$ varying: Analytical curves (-), numerical curve (o). All The parameters are given in [21].

Also, the stochastic analysis has allowed to estimate the probabilistic distribution as a consequence of the wind random effects. The beam response, and more specifically the stationary probability density function $P_s\left(a\right)$ of its amplitude $a$ , can also be retrieved ( Figures 4 and 5 ). This type of analysis indicates that as the additive wind turbulence parameter increases, the peak value of the probability density function decreases and progressively shifts toward larger amplitude values, while the average center position stays in the same position. Thus, the additive ( ${\theta }_0$ ) and parametric ( ${\theta }_1$ ) wind turbulence decreases the chance for the beam to quickly reach the amplitude resonance. It is also demonstrated that the PDF has only one maximum situated in the vicinity of $a_m=0.2$ .

Figure 4.

Stationary probability distribution function as function of the amplitude for different values of the additive wind turbulence parameter. Analytical curves (-), numerical curve (o). All The parameters are given in [21].

Figure 5.

Stationary probability distribution function as function of the amplitude for different values of the parametric wind turbulence parameter. Analytical curves (-), numerical curve (o). All The parameters are given in [21].

We have plotted curves Figure 6(a) and (b) that presents the stationary probability distribution function $P_s\left(a,\varphi \right)$ versus the amplitude $a$ and the phase $\varphi$ . This graph just confirm the results obtained in Figures 4 and 5 and the highest value of the PDF is more visible.

Figure 6.

Stationary probability distribution function versus the amplitude $a$ and the phase $\varphi$ for (a) ${\theta }_0=0.09$ and (b)  ${\theta }_0=0.2$ . All The parameters are given in [21].

Figure 7 presents the times histories of the maximum vibration of the beam. The case where the beam is subjected to moving loads ( $a$ ), to wind actions ( $b$ ) and to the both wind and loads ( $c$ ).

Figure 7.

Vibration amplitude of the beam χ(τ) as function of the time τ. All The parameters are given in [21].

5. Conclusions

In this Chapter we have revised some aspects of the response of viscoelastic foundation of bridges to a simplified model of moving loads and wind random perturbations. The results have been compared to the numerical solution for the modal equation, obtained with the deterministic and stochastic version of Newton-Leipnik algorithm. The analysis has begun modeling the steady-state vibration of the beam suspensions made of a fractional-order viscoelastic material. The resulting mathematical model consists of a component for the beam, and the Kelvin-Voigt foundation type containing fractional derivative of real order, as well as a stochastic term to account for wind pressure. We have highlighted the simplifications. Perhaps the most significant, in the very model formulation, has been the assumption that the load passage consists of concentrated masses, spatially periodic and moving at constant speed. This simplification is crucial to reduce the full partial differential equation to a single smode model — Wind load is modeled as the aerodynamic force related to the wind that blows orthogonally to the beam axis with random velocity. The whole system has then been modeled with a partial differential equation that can be reduced to a one-dimensional modal equation. The beam response under moving and/or stochastic wind loads has been estimated analytically assuming that the first mode contains the essential information, and using the stochastic averaging method. The analysis has therefore some limitations, namely the limited values of parameters that have been explored, and to have retained the first mode only in the Galerkin’s method. Also, the vehicles’ train has been (over?) simplified as a simple periodic drive. With this limitations in mind, let us summarize the main findings. To start with, in this framework it is possible to investigate how the main parameters of the moving loads and of the bearings affect the beam response, and especially how the driving frequency, the loads number, the stiffness coefficient, fractional-order of the viscosity term and the number of bearings affect the dynamic behaviour of the beam. The resonance phenomenon and the stability in the beam system strongly depends on the stiffness and fractional-order of the derivative term of the viscous properties of the bearings. There are a number of quantitative results that are worth mentioning. Firstly, as the number of moving loads increases, the resonant amplitude of the beam increases as well. Secondly, it has been established that as the number of bearings increases, the resonant amplitude decreases and, more importantly, shifts toward larger frequency values. Thirdly, the system response becomes more stable as the order of the derivative increases, for the multivalued solution only appears for the smallest order and quickly disappears as the order increases. All the above results are a consequence of the analysis of the oscillations. However, the stochastic analysis has allowed to estimate the probabilistic distribution as a consequence of the wind random effects. The beam response, and more specifically the stationary probability density function of its amplitude, can also be retrieved. This type of analysis indicates that as the additive wind turbulence parameter increases, the peak value of the probability density function decreases and progressively shifts toward larger amplitude values, while the average center position stays in the same position. Thus, the additive and parametric wind turbulence decreases the chance for the beam to quickly reach the amplitude resonance.

Numerical simulations have confirmed these predictions. This behavior is depicted in Figures 2 – 5 , which have practical implications, on which we would like to comment. To make an example, the beam system frequency $\Omega$ displays the bridge response as the vehicles speed changes, see Eq. (7). In principle it could be possible to avoid large oscillations by controlling the speed of the freight vehicles, albeit in practice it is more realistic to set a maximum speed at which the bridge should be crossed. In other words, to keep resonance at bay, it is necessary to set a speed limit below the resonance insurgence. Analogously, one could think to limit the vehicles number, not to have a minimum number of vehicles across the bridge. We conclude the parameters analysis, namely the stiffness and the viscoelastic properties of the foundations, noticing that such parameters can be optimized with an appropriated tuning, see e.g., Figure 3 , or the analogous indications that stems from the results of Figures 4 and 5 for the wind features ${\theta }_0$ and ${\theta }_1$ .

By way of conclusion, let us summarize that the special properties of the viscoelastic foundations and of the time dependent perturbations, vehicles and wind, interact. As a result also the construction and management parameters are not to be considered independent procedures, for they are deeply interwoven if safe transportation is to be guaranteed.

Acknowledgments

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.

To deal with the modelling, let us consider the dynamic equilibrium of a beam element of length $\mathit{dx}$ ; $w=w\left(x,t\right)$ and $\theta =\theta \left(x,t\right)$ be the transversal displacement and the angle of rotation of the beam element respectively. We denote the internal bending moment by $M$ , the internal shear force by $V$ , the inplane tension due to the inplane strain, issue of the assumed negligible longitudinal displacement of the beam by $T$ , the foundation-beam interaction force (per unit length of the beam’s axis) by $Q_F\left(x,t\right)$ and the external distributed loading by $F_{\mathit{ad}}\left(x,t\right)$ and $f\left(x,t\right)$ .

Setting the vertical forces on the element equal to the mass times acceleration gives:

$$\frac{\partial V}{\partial x}=Q_F\left(x,t\right)-f\left(x,t\right)-F_{\mathit{ad}}\left(x,t\right)+\mathit{\rho S}\frac{{\partial }^{2}w\left(x,t\right)}{\partial t^2}$$ E41

While summing moments produces:

$$\frac{\partial M}{\partial x}=V-\mathit{\rho I}\frac{{\partial }^2\theta \left(x,t\right)}{\partial t^2}-T\frac{\partial w\left(x,t\right)}{\partial x}$$ E42

For small rotation $\theta \left(x,t\right)\approx \frac{\partial w\left(x,t\right)}{\partial x}$ , Eq. (42) becomes:

$$\frac{\partial M}{\partial x}=V-\mathit{\rho I}\frac{{\partial }^{3}w\left(x,t\right)}{\partial t^2\partial x}-T\frac{\partial w\left(x,t\right)}{\partial x}$$ E43

Combining Eq. (41) and Eq. (43) then yields:

$$\frac{{\partial }^{2}M}{\partial x^2}=Q_F\left(x,t\right)-f\left(x,t\right)-F_{\mathit{ad}}\left(x,t\right)+\mathit{\rho S}\frac{{\partial }^{2}w\left(x,t\right)}{\partial t^2}-\mathit{\rho I}\frac{{\partial }^{4}w\left(x,t\right)}{\partial t^2\partial x^2}-T\frac{{\partial }^{2}w\left(x,t\right)}{\partial x^2}$$ E44

From the geometry of the deformation, and using Hooke’s law ${\sigma }_x=E{\epsilon }_x$ , one can show that (see reference [19]):

$$\begin{array}{l}M=\frac{\mathit{EI}}{R}=-\mathit{EI}\frac{\frac{{\partial }^{2}w\left(x,t\right)}{\partial x^2}}{{\left[1+{\left(\frac{\partial w\left(x,t\right)}{\partial x}\right)}^2\right]}^{\frac{3}{2}}}\approx -\mathit{EI}\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]+O\left({\left(\frac{\partial w\left(x,t\right)}{\partial x}\right)}^2\right)\\ \\ \approx -\mathit{EI}\frac{{\partial }^{4}w\left(x,t\right)}{\partial x^4}+\frac{3}{2}\mathit{EI}\frac{{\partial }^2}{\partial x^2}\left[\frac{{\partial }^{2}w\left(x,t\right)}{\partial x^2}{\left(\frac{\partial w\left(x,t\right)}{\partial x}\right)}^2\right]+O\left({\left(\frac{\partial w}{\partial x}\right)}^2\right)\end{array}.$$ E45

where the Taylor expansion of the inverse of the radius of curvature $\left(\frac{1}{R}\right)$ up to the second order is carried out. According to the assumed negligible longitudinal displacement of the beam, the tension in the beam $T$ can be determined as (see the details of their derivation in Ref.[19]).

$$T=\frac{\mathit{ES}}{2L}\underset{0}{\overset{L}{\int }}{\left(\frac{\partial w\left(x,t\right)}{\partial x}\right)}^2\mathrm{d}x$$ E46

Finally taking into account the dissipation ( $\mu \frac{\partial w\left(x,t\right)}{\partial t}$ ), putting Eq. (44), Eq. (45) and Eq. (46) together gives the new desired result (Eq. (1) of the manuscript)

$$\begin{array}{l}\mathit{\rho S}\frac{{\partial }^{2}w\left(x,t\right)}{\partial t^2}-\mathit{\rho I}\frac{{\partial }^{4}w\left(x,t\right)}{\partial x^2\partial t^2}+\mathit{EI}\frac{{\partial }^{4}w\left(x,t\right)}{\partial x^4}+\mu \frac{\partial w\left(x,t\right)}{\partial t}-\frac{3}{2}\mathit{EI}\frac{{\partial }^2}{\partial x^2}\left[\frac{{\partial }^{2}w\left(x,t\right)}{\partial x^2}{\left(\frac{\partial w\left(x,t\right)}{\partial x}\right)}^2\right]\\ \\ -\frac{\mathit{ES}}{2L}\frac{{\partial }^{2}w\left(x,t\right)}{\partial x^2}\underset{0}{\overset{L}{\int }}{\left(\frac{\partial w\left(x,t\right)}{\partial x}\right)}^2\mathrm{d}x+Q_F\left(x,t\right)=F_{\mathit{ad}}\left(x,t\right)+f\left(x,t\right)\end{array},$$ E47

where:

$$\begin{array}{l}f\left(x,t\right)=P\sum _{i=0}^{N-1}{\epsilon }_i\delta \left[x-x_i\left(t-t_i\right)\right]\\ \\ Q_F\left(x,t\right)=\sum _{j=1}^{N_P}\left(k_j+c_j{D}_t^{{\alpha }_j}\right)w\left(x,t\right)\delta \left[x-\frac{\mathit{jL}}{N_P+1}\right]\end{array}$$