Vibrations of the Euler–Bernoulli Beam Under a Moving Force based on Various Versions of Gradient Nonlocal Elasticity Theory: Application in Nanomechanics

In the classical (local) theory of elasticity, the stress at a given point depends only on the strain at the same point. Many theories have been developed based on this assumption for various types of structures such as rods, beams, plates, and shells. In turn, the experiments associated with nanotechnology demonstrate that the local continuous theory cannot predict the behavior of nanoscale structures. In such structures, the size effect takes place, and for this reason, the nanomaterials are better described by nonlocal continuous theory. The theory of nonlocal continuous mechanics assumes that the stress at a particular point is a function of strains (stresses) at all points in the continuum. Nonlocal elasticity was initiated in articles [1, 2, 3, 4]. In the past few years several problems have been solved using nonlocal continuous elasticity theory, in particular problems connected with the buckling and vibration of beams [5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19]. A dynamic response of the nanotube subjected to a moving nanoparticle is an interesting and important problem. The vibration of different types of nanostructures, such as nanotubes, double-walled carbon nanotubes, and nanoplates, under a moving load are considered in articles [20, 21, 22, 23, 24, 25, 26, 27, 28, 29]. The problem of molecular modeling in nanostructured materials have been considered, among others, in the articles [30, 31]. To describe the impact of small-sized nonlocal properties of materials in the dynamics of the abovementioned nanostructures, the stress gradient model is most often used, although the strain gradient model is also used. It is worth remembering that in experimental studies, strain but not stress are measured. Hence, it follows that comparative theoretical and numerical studies should be conducted simultaneously for both gradient models.

In this article, we study and compare two models of vibration of the Euler–Bernoulli beam under a moving force based on two different gradient versions of the nonlocal theory of elasticity, namely, the nonlocal Eringen’s model, in which the strain is a function of stress gradient, and the nonlocal Aifantis’s model, in which the stress is a function of strains gradient. A comparative analysis of the integral stress-driven model versus strain-driven model is presented in the article [33]. We study the dynamic response of a finite, simply supported beam under a moving force. The force is moving along the beam with a constant velocity. We will present the particular solutions in the form of an infinite series and some solutions in a closed form. We also present the numerical results. The presented theoretical solutions and numerical results can be used in the dynamic analysis of nanotube structures under the excitation of a moving nanoparticle with a weight of mg modeled by a moving constant force P= mg Moreover, the movement of a nanovehicle on the nanobeam can be modeled using a moving force. In the recent years, intensive theoretical and experimental research in the field of nanotechnology has been carried out. The progress in the technology of nanodevices such as nanovehicles is presented, among others, in articles [32,33]. This progress in nanotechnology generates similar problems as observed in civil and mechanical engineering structures loaded with a moving load [34].

The outline of the article is given as follows: in Section 2, the dynamic equation of the Euler–Bernoulli beam subordinate, according to the stress gradient model, is derived. In Section 3, the vibration of the Euler–Bernoulli beam under a force moving with a constant velocity is considered. The classical solution has the form of an infinite series. In this section, it has been shown that the aperiodic part of the solution can also be presented in a closed form instead of in an infinite series. In Section 4, the dynamic response because of the force moving on the Euler–Bernoulli beam, based on the strain gradient nonlocal elasticity model, is considered. In Section 5, the comparative numerical analysis of the dynamic response of the beam because of the moving force for two nonlocal models of the Euler–Bernoulli beam is presented. Finally, conclusions and comments are drawn in Section 6.

In the well-known nonlocal elasticity theory proposed by Eringen [4], the nonlocal constitutive relationship is defined as (1)[1−(e0a)2∇2]σij(x)=λεrr(x)δij+2μεij(x),\left[ {1 - \left( {e_0 a} \right)^2 \nabla ^2 } \right]\sigma _{ij} (x) = \lambda \varepsilon _{rr} \left( x \right)\delta _{ij} + 2\mu \varepsilon _{ij} \left( x \right), where ∇² is the Laplacian, λ and μ are the Lame constants, the symbol δ_ij denotes the Kronecker e₀ is a material constant, and a is an internal characteristic length.

For a one-dimensional element, Eq. (1) takes the following form: (2)[1−(e0a)2∂2∂x2]σ(x)=Eε(x),[ {1 - \left( {e_0 a} \right)^2 \frac{{\partial ^2 }}{{\partial x^2 }}} ]\sigma \left( x \right) = E\varepsilon (x), where E is Young’s modulus.

Using Eq. (2) and the well-known relationship between the strain and the transverse deflection of the Euler–Bernoulli beam, (3)ε(x,t)=zr(x,t)=−z∂2w(x,t)∂x2,\varepsilon \left( {x,t} \right) = \frac{z}{{r\left( {x,t} \right)}} = - z\frac{{\partial ^2 w\left( {x,t} \right)}}{{\partial x^2 }}, one obtains the relationship between the beam deflection and the bending moment (4)EI∂2w(x,t)∂x2=−M(x,t)+(e0a)2∂2M(x,t)∂x2,EI\frac{{\partial ^2 w\left( {x,t} \right)}}{{\partial x^2 }} = - M\left( {x,t} \right) + \left( {e_0 a} \right)^2 \frac{{\partial ^2 M\left( {x,t} \right)}}{{\partial x^2 }}, where I is the moment of inertia and M is the resultant bending moment, which is defined by (5)M(x,t)=∫Aσ(x,t)zdA.M\left( {x,t} \right) = \int\limits_A {\sigma (x,t)zdA.}

By considering the equilibrium of an infinitesimal element of length dx and by applying D’Alembert’s principle, the equilibrium equation takes the following form: (6)−∂2M(x,t)∂x2+ρA∂2w(x,t)∂t2=p(x,t), - \frac{{\partial ^2 M(x,t)}}{{\partial x^2 }} + \rho A\frac{{\partial ^2 w(x,t)}}{{\partial t^2 }} = p(x,t), where p (x, t) is the load process, ρ is the material density, and A is the cross-sectional area.

Let us consider a beam of finite length L simply supported on both ends.

Considering the set of differential equations (4) and (6), the boundary conditions for a simply supported beam have the following form: (7)w(0,t)=w(L,t)=0,w(0,t) = w(L,t) = 0,(8)M(0,t)=M(L,t)=0.M(0,t) = M(L,t) = 0.

After combining Eqs. (4) and (6), one obtains a differential equation beam in the following form: (9)EI∂4w(x,t)∂x4+ρA∂2w(x,t)∂t2−−(e0a)2ρA∂4w(x,t)∂x2∂t2=p(x,t)−(e0a)2∂2p(x,t)∂x2.\matrix{{EI\frac{{\partial ^4 w(x,t)}}{{\partial x^4 }} + \rho A\frac{{\partial ^2 w(x,t)}}{{\partial t^2 }}} \\ { - \left( {e_0 a} \right)^2 \rho A\frac{{\partial ^4 w(x,t)}}{{\partial x^2 \partial t^2 }} = p(x,t) - \left( {e_0 a} \right)^2 \frac{{\partial ^2 p(x,t)}}{{\partial x^2 }}.} }

In this case, the boundary conditions have the following forms: (10)w(0,t)=w(L,t)=0,w(0,t) = w(L,t) = 0,(11)∂2w(x,t)∂x2|x=0=∂2w(x,t)∂x2|x=L=0.\left. {{{\partial ^2 w(x,t)} \over {\partial x^2 }}} \right|_{x = 0} = \left. {{{\partial ^2 w(x,t)} \over {\partial x^2 }}} \right|_{x = L} = 0.

The boundary condition (11) requires some comment. This is correct in a simply supported beam according to the classical theory of elasticity. In the Eringen model, it is not fulfilled for various loads [see (40)]. It is widely accepted in dynamic issues because of the possibility of presenting a solution in the form of a Fourier series. Some additional comment on this problem is given in the final conclusions.

Let us consider the vibrations of a beam excited by a point force P moving with a constant velocity v as in Figure 1

Figure 1

In this case, the load function in Eqs. (6) and (9) has the following form: (12)p(x,t)=Pδ(x−vt),p(x,t) = P\delta (x - vt), where δ(.) is the Dirac delta.

After introducing the dimensionless variables, (13)ξ=xL, T=vtL, ξ∈[0,1], T∈[0,1],\xi = {x \over L},\;\;\;T = {{vt} \over L},\;\;\;\xi \in [0,1],\;\;\;T \in [0,1],

Eqs. (4), (6), and (9) take the following forms: (14)EI∂2w(ξ,T)∂ξ2+L2M(ξ,T)−μe2L2∂2M(ξ,T)∂ξ2=0,EI\frac{{\partial ^2 w(\xi ,T)}}{{\partial \xi ^2 }} + L^2 M(\xi ,T) - \mu _e^2 L^2 \frac{{\partial ^2 M(\xi ,T)}}{{\partial \xi ^2 }} = 0,(15)−∂2M(ξ,T)∂ξ2+Aρv2∂2w(ξ,T)∂T2=PLδ(ξ−T), - \frac{{\partial ^2 M(\xi ,T)}}{{\partial \xi ^2 }} + A\rho v^2 \frac{{\partial ^2 w(\xi ,T)}}{{\partial T^2 }} = PL\delta (\xi - T), and (16)=P0δ(ξ−T)−μe2P0∂2δ(ξ−T)∂ξ2, = P_0 \delta (\xi - T) - \mu _e^2 P_0 \frac{{\partial ^2 \delta (\xi - T)}}{{\partial \xi ^2 }}, where σ2=ρAv2L2EI=ηg2λ2\sigma ^2 = \frac{{\rho Av^2 L^2 }}{{EI}} = \eta _g^2 \lambda ^2, ηg=vvg\eta _g = \frac{v}{{v_g }}, vg=Eρv_g = \sqrt {\frac{E}{\rho }}, λ=Lr\lambda = \frac{L}{r}, r=IAr = \sqrt {\frac{I}{A}} , P0=PL3EIP_0 = \frac{{{PL}^3 }}{{{EI}}}, μe=e0aL\mu _e = \frac{{e_0 a}}{L}.

Instead of solving Eq. (16), it is convenient to solve the set of the two equations (14) and (15). In the latter case, we obtain both the solution for the vertical displacement and the bending moment.

The boundary conditions have the following form: (17)w(0,T)=w(1,T)=0,w(0,T) = w(1,T) = 0,(18)M(0,T)=M(1,T)=0.M(0,T) = M(1,T) = 0.

Let the initial conditions take the following form: (19)w(ξ,0)=0, ∂w(ξ,T)∂T|T=0=0.w(\xi ,0) = 0,\,\,\,\,\,\,\,\,\left. {\frac{{\partial w(\xi ,T)}}{{\partial T}}} \right|_{T = 0} = 0.

The solution of Eqs. (14) and (15) for the boundary conditions (17) and (18) is assumed to be in the form of the sine series (20)w(ξ,T)=∑n=1∞yn(T)sin nπξ,w(\xi ,T) = \sum\limits_{n = 1}^\infty {y_n (T)\sin \,n\pi \xi ,} (21)M(ξ,T)=∑n=1∞sn(T)sin nπξ.M(\xi ,T) = \sum\limits_{n = 1}^\infty {s_n (T)\sin \,n\pi \xi .}

After substituting expressions (20) and (21) into Eqs. (14) and (15) and by using the orthogonality method, we obtain the set of ordinary differential equations (22)EI(nπ)2yn(T)−L2[1+μe2(nπ)2]sn(T)=0,{EI}\left( {n\pi } \right)^2 y_n (T) - L^2 [ {1 + \mu _e^2 \left( {n\pi } \right)^2 } ]s_n (T) = 0,(23)d2yn(T)dT2+(nπ)2Aρv2zn(T)=2PLAρv2sin nπT.\frac{{d^2 y_n (T)}}{{dT^2 }} + \frac{{\left( {n\pi } \right)^2 }}{{A\rho v^2 }}z_n (T) = \frac{{2PL}}{{A\rho v^2 }}{\rm{sin}}\,n\pi T.

After eliminating the function z_n (T), one obtains (24)d2yn(T)dT2+ωn2yn(T)=2PLAρv2sin nπT,\frac{{d^2 y_n (T)}}{{dT^2 }} + \omega _n^2 y_n (T) = \frac{{2PL}}{{A\rho v^2 }}{\rm{sin}}\,n\pi T, where the natural frequency in the dimensionless variables (13) is equal to (25)ωn2=EI(nπ)4Aρv2L2[1+(nπ)2μe2]=(nπ)4ηg2λ2[1+(nπ)2μe2].\omega _n^2 = \frac{{EI\left( {n\pi } \right)^4 }}{{A\rho v^2 L^2 [ {1 + \left( {n\pi } \right)^2 \mu _e^2 } ]}} = \frac{{\left( {n\pi } \right)^4 }}{{\eta _g^2 \lambda ^2 [ {1 + \left( {n\pi } \right)^2 \mu _e^2 } ]}}.

The initial condition has the following form: (26)yn(0)=0, dyn(T)dT|T=0=0.y_n (0) = 0,\,\,\,\,\,\,\,\,\,\,\,\left. {\frac{{dy_n (T)}}{{dT}}} \right|_{T = 0} = 0.

Finally, the solutions of the system of Eqs. (14) and (15) are the sums of the particular integrals w_A(ξ, T), M_A(ξ, T), and the general integrals w_S(ξ, T) and M_S(ξ, T) have the following forms: (27)w(ξ,T)=wA(ξ,T)+wS(ξ,T),w(\xi ,T) = w_A (\xi ,T) + w_S (\xi ,T),(28)M(ξ,T)=MA(ξ,T)+MS(ξ,T),M(\xi ,T) = M_A (\xi ,T) + M_S (\xi ,T), where (29)wA(ξ,T)=2PL3EI∑n=1∞[1+(nπ)2μe2]sin nπT sin nπξ(nπ)2[(nπ)2(1−σ2μe2)−σ2],w_A (\xi ,T) = \frac{{2PL^3 }}{{EI}}\sum\limits_{n = 1}^\infty {\frac{{[ {1 + \left( {n\pi } \right)^2 \mu _e^2 } ]{\rm{sin}}\,n\pi T\,{\rm{sin}}\,n\pi \xi }}{{\left( {n\pi } \right)^2 [ {\left( {n\pi } \right)^2 \left( {1 - \sigma ^2 \mu _e^2 } \right) - \sigma ^2 } ]}},} (30)wS(ξ,T)=−2PL3EI∑n=1∞[1+(nπ)2μe2]sin ωnT sin nπξωn(nπ)[(nπ)2(1−σ2μe2)−σ2],w_S (\xi ,T) = \frac{{ - 2PL^3 }}{{EI}}\sum\limits_{n = 1}^\infty {\frac{{[ {1 + \left( {n\pi } \right)^2 \mu _e^2 } ]{\rm{sin}}\,\omega _n T\,{\rm{sin}}\,n\pi \xi }}{\omega _n{\left( {n\pi } \right)^2 [ {\left( {n\pi } \right)^2 \left( {1 - \sigma ^2 \mu _e^2 } \right) - \sigma ^2 } ]}},} (31)MA(ξ,T)=2PL∑n=1∞sin nπT sin nπξ(nπ)2(1−σ2μe2)−σ2,M_A (\xi ,T) = 2PL\sum\limits_{n = 1}^\infty {\frac{{{\rm{sin}}\,n\pi T\,{\rm{sin}}\,n\pi \xi }}{{\left( {n\pi } \right)^2 \left( {1 - \sigma ^2 \mu _e^2 } \right) - \sigma ^2 }},} (32)MS(ξ,T)=−2PL∑n=1∞nπ sin ωnT sin nπξωn[(nπ)2(1−σ2μe2)−σ2].M_S (\xi ,T) = - 2PL\sum\limits_{n = 1}^\infty {\frac{{n\pi \,{\rm{sin}}\,\omega _n T{\rm{sin}}\,n\pi \xi }}{{\omega _n [ {\left( {n\pi } \right)^2 \left( {1 - \sigma ^2 \mu _e^2 } \right) - \sigma ^2 } ]}}.}

Let us note that when the condition (33)π2(1−σ2μe2)−σ2=0\pi ^2 \left( {1 - \sigma ^2 \mu _e^2 } \right) - \sigma ^2 = 0 is fulfilled, solutions (29)–(32) tend to infinity.

Thus the resonance velocity v_ecr is equal to (34)vecr=πLEIAρ(1+π2μe2)=πλ1+π2μe2vg.v_{{ecr}} = \frac{\pi }{L}\sqrt {\frac{{{EI}}}{{A\rho \left( {1 + \pi ^2 \mu _e^2 } \right)}}} = \frac{\pi }{{\lambda \sqrt {1 + \pi ^2 \mu _e^2 } }}v_g .

The functions w_A(ξ, T) and M_A(ξ, T) describe aperiodic vibrations and satisfy the nonhomogeneous differential equations. (14) and (15). These functions do not satisfy the initial conditions of motion (19). The functions w_S(ξ, T) and M_S(ξ, T) correspond to free vibrations of the beam and satisfy the homogeneous differential equations (14), (15), and (P = 0), and together with the aperiodic functions, the initial conditions of motion are satisfied. Now we will present the aperiodic solutions w_A(ξ, T) and M_A(ξ, T) given by expressions (29) and (31) in closed analytical form.

Let us notice an important fact that these functions are not only solutions to the system of partial differential equations (14) and (15) but also to the system of ordinary equations (see [41, 42, 43]) (35)EId2wA(ξ,T)dξ2+L2MA(ξ,T)−μe2L2d2MA(ξ,T)dξ2=0,{EI}\frac{{d^2 w_A (\xi ,T)}}{{d\xi ^2 }} + L^2 M_A (\xi ,T) - \mu _e^2 L^2 \frac{{d^2 M_A (\xi ,T)}}{{d\xi ^2 }} = 0,(36)−d2MA(ξ,T)dξ2+Aρv2d2wA(ξ,T)dξ2=PLδ(ξ−T), - \frac{{d^2 M_A (\xi ,T)}}{{d\xi ^2 }} + A\rho v^2 \frac{{d^2 w_A (\xi ,T)}}{{d\xi ^2 }} = PL\delta (\xi - T), for the boundary conditions (17) and (18).

The variable T in Eq. (36) is the only parameter that describes the location of the moving force on the beam. After solving the set of Eqs. (35) and (36) using, for example, the Laplace transformation, we can obtain the functions w_A(ξ, T) and M_A(ξ, T) in a closed form instead of in a series. The closed form of the solutions depends on the velocity of the moving force.

In the case of σμe<1 (v<vgμeλ)\sigma \mu _e < 1\,\,\,( {v < \frac{{v_g }}{{\mu _e \lambda }}}), the solutions have the following form: (37)wA(ξ,T)=−PL3EIσ2ξ(1−T)+PL3(1+β2μe2)EIσ2sin β(1−T) sin βEβsinβξ≤T,\matrix{ {w_A (\xi ,T) = - \frac{{PL^3 }}{{EI\sigma ^2 }}\xi (1 - T) + \frac{{PL^3 \left( {1 + \beta ^2 \mu _e^2 } \right)}}{{EI\sigma ^2 }}\frac{{{\rm{sin}}\,\beta (1 - T)\,{\rm{sin}}\,\beta E}}{{\beta \sin \beta }}} \\ {\xi \le T,} }(38)wA(ξ,T)=−PL3EIσ2(1−ξ)T+2PL3(1+β2μe2)EIσ2sin βT sin β(1−ξ)βsinβξ≥T,\matrix{{w_A (\xi ,T) = - \frac{{PL^3 }}{{EI\sigma ^2 }}(1 - \xi )T + \frac{{2PL^3 \left( {1 + \beta ^2 \mu _e^2 } \right)}}{{EI\sigma ^2 }}\frac{{{\rm{sin}}\,\beta T\,{\rm{sin}}\,\beta (1 - \xi )}}{{\beta \sin \beta }}} \\ {\xi \ge T,}} and (39)MA(ξ,T)=PLβσ2sin β(1−T) sin βξsin β ξ≤T,M_A (\xi ,T) = PL\frac{\beta }{{\sigma ^2 }}\frac{{{\rm{sin}}\,\beta (1 - T)\,{\rm{sin}}\,\beta \xi }}{{{\rm{sin}}\,\beta }}\,\,\,\,\xi \le T,(40)MA(ξ,T)=PLβσ2sin βT sin β(1−ξ)sin β ξ≥T,M_A (\xi ,T) = PL\frac{\beta }{{\sigma ^2 }}\frac{{{\rm{sin}}\,\beta T\,{\rm{sin}}\,\beta (1 - \xi )}}{{{\rm{sin}}\,\beta }}\,\,\,\,\xi \ge T, where β2=σ21−σ2μe2\beta ^2 = \frac{{\sigma ^2 }}{{1 - \sigma ^2 \mu _e^2 }}.

In the case of σμe>1 (v>vgμeλ)\sigma \mu _e > 1\,\,\,( {v > \frac{{v_g }}{{\mu _e \lambda }}} ), the solutions have the following form: (41)wA(ξ,T)=−PL3EIσ2ξ(1−T)+PL3(1−β¯2μe2)EIσ2sinh β¯(1−T) sinh β¯ξβ¯ sinh β¯ξ≤T,\matrix{{w_A (\xi ,T) = - \frac{{PL^3 }}{{EI\sigma ^2 }}\xi (1 - T) + \frac{{PL^3 \left( {1 - \bar \beta ^2 \mu _e^2 } \right)}}{{EI\sigma ^2 }}\frac{{{\rm{sinh}}\,\bar \beta (1 - T)\,{\rm{sinh}}\,\bar \beta \xi }}{{\bar \beta\, {\rm{sinh}}\,\bar \beta }}} \\ {\xi \le T,}}(42)wA(ξ,T)=−PL3EIσ2(1−ξ)T+2PL3(1−β¯2μe2)EIσ2sinh β¯T sinh β¯(1−ξ)β¯ sinh β¯ξ≥T,\matrix{{w_A (\xi ,T) = - \frac{{PL^3 }}{{EI\sigma ^2 }}(1 - \xi )T + \frac{{2PL^3 \left( {1 - \bar \beta ^2 \mu _e^2 } \right)}}{{EI\sigma ^2 }}\frac{{{\rm{sinh}}\,\bar \beta T\,{\rm{sinh}}\,\bar \beta (1 - \xi )}}{{\bar \beta \,{\rm{sinh}}\,\bar \beta }}} \\ {\xi \ge T,}} and (43)MA(ξ,T)=−PLβ¯σ2sinh β¯(1−T) sinh β¯ξsinh β¯ ξ≤T,M_A (\xi ,T) = - PL\frac{{\bar \beta }}{{\sigma ^2 }}\frac{{{\rm{sinh}}\,\bar \beta (1 - T){\rm{sinh}}\,\bar \beta \xi }}{{{\rm{sinh}}\,\bar \beta }}\,\,\,\,\,\,\,\,\,\,\xi \le T,(44)MA(ξ,T)=−PLβ¯σ2sinh β¯T sinh β¯(1−ξ)sinh β¯ ξ≥T,M_A (\xi ,T) = - PL\frac{{\bar \beta }}{{\sigma ^2 }}\frac{{{\rm{sinh}}\,\bar \beta T{\rm{sinh}}\,\bar \beta (1 - \xi )}}{{{\rm{sinh}}\,\bar \beta }}\,\,\,\,\,\,\,\,\,\,\xi \ge T, where β¯2=σ2σ2μe2−1\bar \beta ^2 = \frac{{\sigma ^2 }}{{\sigma ^2 \mu _e^2 - 1}}.

Let us assume that the beam is loaded in a static way by a force in the point ξ₀. In this case, the displacement of the beam and the bending moment have the following form: (45)we(ξ,ξ0)=2PL3EI∑n=1∞[1+(nπ)2μe2]sin nπξ0 sin nπξ(nπ)4,w_e (\xi ,\xi _0 ) = \frac{{2PL^3 }}{{EI}}\sum\limits_{n = 1}^\infty {\frac{{[ {1 + \left( {n\pi } \right)^2 \mu _e^2 } ]{\rm{sin}}\,n\pi \xi _0 \,{\rm{sin}}\,n\pi \xi }}{{\left( {n\pi } \right)^4 }},} (46)Me(ξ0,ξ)=2PL∑n=1∞sin nπξ0 sin nπξ(nπ)2,M_e (\xi _0 ,\xi ) = 2PL\sum\limits_{n = 1}^\infty {\frac{{{\rm{sin}}\,n\pi \xi _0 \,{\rm{sin}}\,n\pi \xi }}{{\left( {n\pi } \right)^2 }},} or in a closed form (47)we(ξ,ξ0)=PL3{−(1−ξ0)ξ3+16(1−ξ0)ξ[1−(1−ξ0)2]+μe2(1−ξ0)ξ+[16(ξ−ξ0)3−μe2(ξ−ξ0)]H(ξ−ξ0)},w_e (\xi ,\xi _0 ) = PL^3 \{ { - \left( {1 - \xi _0 } \right)\xi ^3 + \frac{1}{6}\left( {1 - \xi _0 } \right)\xi [ {1 - \left( {1 - \xi _0 } \right)^2 } ] + \mu _e^2 (1 - \xi _0 )\xi + } {[ {\frac{1}{6}\left( {\xi - \xi _0 } \right)^3 - \mu _e^2 \left( {\xi - \xi _0 } \right)} ]H\left( {\xi - \xi _0 } \right)} \}, (48)Me(ξ0,ξ)=PL[ξ(1−ξ0)−(ξ−ξ0)H(ξ−ξ0)],M_e (\xi _0 ,\xi ) = {PL}[ {\xi (1 - \xi _0 ) - (\xi - \xi _0 )H(\xi - \xi _0 )} ],

The above solutions enable the dynamic effects caused by the movement of force when loaded statically and dynamically to be analyzed.

Let us return again to the relationship between stresses and strains. In this case, we assume for a one-dimensional system the strain gradient model in the form of (see [36]–[39]) (49)σ=E(1−δ2∂2∂x2)ε,\sigma = E( {1 - \delta ^2 \frac{{\partial ^2 }}{{\partial x^2 }}} )\varepsilon , where the constant δ represents nonlocal effects.

Taking into account Eq. (49), the relation between the bending moment M(x, t) and the beam displacement w(x, t) has the form

Taking into account equilibrium equation (6), the equation of motion for the strain gradient nonlocal Euler–Bernoulli beam model has the following form: (51)−EIδ2∂6w(x,t)∂x6+EI∂4w(x,t)∂x4+ρA∂2w(x,t)∂t2=p(x,t). - EI\delta ^2 \frac{{\partial ^6 w(x,t)}}{{\partial x^6 }} + EI\frac{{\partial ^4 w(x,t)}}{{\partial x^4 }} + \rho A\frac{{\partial ^2 w(x,t)}}{{\partial t^2 }} = p(x,t).

After introducing the dimensionless variable (13) and assuming that the load process is a point force moving with the constant velocity (12), the above equation can be presented in the following form: (52)−∂6w(ξ,T)∂ξ6+1μs2∂4w(ξ,T)∂ξ4+ηg2λ2μs2∂2w(ξ,T)∂T2=PL3μs2EIδ(ξ−T), - \frac{{\partial ^6 w(\xi ,T)}}{{\partial \xi ^6 }} + \frac{1}{{\mu _s^2 }}\frac{{\partial ^4 w(\xi ,T)}}{{\partial \xi ^4 }} + \frac{{\eta _g^2 \lambda ^2 }}{{\mu _s^2 }}\frac{{\partial ^2 w(\xi ,T)}}{{\partial T^2 }} = \frac{{PL^3 }}{{\mu _s^2 EI}}\delta (\xi - T), where μs=δL\mu _s = \frac{\delta }{L}.

For a simply supported beam, the boundary conditions have the following form: (53)w(0,T)=w(1,T)=0,w(0,T) = w(1,T) = 0,(54)∂2w(ξ,T)∂ξ2|ξ=0=∂2w(ξ,t)∂ξ2|ξ=1=0,\left. {\frac{{\partial ^2 w(\xi ,T)}}{{\partial \xi ^2 }}} \right|_{\xi = 0} = \left. {\frac{{\partial ^2 w(\xi ,t)}}{{\partial \xi ^2 }}} \right|_{\xi = 1} = 0,(55)∂4w(ξ,T)∂ξ4|ξ=0=∂4w(ξ,t)∂ξ4|ξ=1=0.\left. {\frac{{\partial ^4 w(\xi ,T)}}{{\partial \xi ^4 }}} \right|_{\xi = 0} = \left. {\frac{{\partial ^4 w(\xi ,t)}}{{\partial \xi ^4 }}} \right|_{\xi = 1} = 0.

The above partial differential equation (51) is of the sixth order and requires six boundary conditions. Compliance with dependencies (54) and (55) allows meeting (satisfy) the boundary condition (8) by taking into account that the bending moment determined by formula (50) contains the second- and fourth-order derivative of beam displacement.

The solution of Eq. (52) for the boundary conditions (53)–(55) is assumed to be in the form of a sine series (12), and after substituting into Eq. (30) and using the orthogonality method, we obtain a set of uncoupled equations (56)d2yn(T)dt2+ωsn2yn(T)=2P0ηg2λ2sin nπT,\frac{{d^2 y_n (T)}}{{dt^2 }} + \omega _{sn}^2 y_n (T) = 2\frac{{P_0 }}{{\eta _g^2 \lambda ^2 }}{\rm{sin}}\,n\pi T, where ωsn=(nπ)2ηgλ1+μs2(nπ)2\omega _{sn} = \frac{{\left( {n\pi } \right)^2 }}{{\eta _g \lambda }}\sqrt {1 + \mu _s^2 \left( {n\pi } \right)^2 }, P0=PL3EIP_0 = \frac{{PL^3 }}{{EI}}.

Finally, the solution of Eq. (52) is the sum of the particular integrals w_A(ξ, T) and the general integral w_S(ξ, T): (57)w(ξ,T)=wA(ξ,T)+wS(ξ,T),w(\xi ,T) = w_A (\xi ,T) + w_S (\xi ,T), where (58)wA(ξ,T)=2P0∑n=1∞sin nπT sin nπξ(nπ4)[1+μs2(nπ)2]−ηg2λ2(nπ)2,w_A (\xi ,T) = 2P_0 \sum\limits_{n = 1}^\infty {\frac{{{\rm{sin}}\,n\pi T\,{\rm{sin}}\,n\pi \xi }}{{\left( {n\pi ^4 } \right)[ {1 + \mu _s^2 \left( {n\pi } \right)^2 } ] - \eta _g^2 \lambda ^2 \left( {n\pi } \right)^2 }},} (59)wS(ξ,T)=−2P0∑n=1∞sin ωsnT sin nπξωsn{(nπ)3[1+μs2(nπ)2]−ηg2λ2(nπ)2}.w_S (\xi ,T) = - 2P_0 \sum\limits_{n = 1}^\infty {\frac{{{\rm{sin}}\,\omega _{sn} T\,{\rm{sin}}\,n\pi \xi }}{{\omega _{sn} \{ {\left( {n\pi } \right)^3[ {1 + \mu _s^2 \left( {n\pi } \right)^2 } ] - \eta _g^2 \lambda ^2 \left( {n\pi } \right)^2 } \}}}.}

Let us notice that when the condition (60)π2(1+π2μs2)−v2vg2λ2=0\pi ^2 ( {1 + \pi ^2 \mu _s^2 } ) - \frac{{v^2 }}{{v_g^2 }}\lambda ^2 = 0 is fulfilled, the solutions (58) and (59) tend to infinity. Thus the resonance velocity v_scr is equal to (61)vscr=πλ1+μs2π2vg.v_{{scr}} = \frac{\pi }{\lambda }\sqrt {1 + \mu _s^2 \pi ^2 } v_g .

Taking into account the relationships (50) and (57)–(59), the bending moment has the following form: (62)M(ξ,T)=MA(ξ,T)+MS(ξ,T),M(\xi ,T) = M_A (\xi ,T) + M_S (\xi ,T), where (63)MA(ξ,T)=2PL{∑n=1∞sin⁡nπTsin⁡nπξ(nπ)2[1+μS2(nπ)2]−ηg2λ2++ μS2∑n=1∞(nπ)2sin⁡nπTsin⁡nπξ(nπ)2[1+μS2(nπ)2]−ηg2λ2},\begin{array}{l} M_A (\xi ,T) = 2PL\left\{ {\sum\limits_{n = 1}^\infty {\frac{{{\rm{sin}}\,n\pi T\,{\rm{sin}}\,n\pi \xi }}{{\left( {n\pi } \right)^2[ {1 + \mu _s^2 \left( {n\pi } \right)^2 } ] - \eta _g^2 \lambda ^2 }} + } } \right. \\ \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\left. { + \mu _S^2 \sum\limits_{n = 1}^\infty {\frac{{\left( {n\pi } \right)^2 {\rm{sin}}\,n\pi T\,{\rm{sin}}\,n\pi \xi }}{{\left( {n\pi } \right)^2 [ {1 + \mu _S^2 \left( {n\pi } \right)^2 } ] - \eta _g^2 \lambda ^2 }}} } \right\}, \\ \end{array}(64)MS(ξ,T)=−2PL{∑n=1∞sin⁡ωsnTsin⁡nπξωsn{(nπ)[1+μS2(nπ)2]−ηg2λ2}++ μS2∑n=1∞(nπ)2sin⁡ωsnTsin⁡nπξωsn{(nπ)[1+μs2(nπ)2]−ηg2λ2}}.\begin{array}{l} M_S (\xi ,T) = 2PL\left\{ {\sum\limits_{n = 1}^\infty {\frac{{{\rm{sin}}\,\omega _{sn} T\,{\rm{sin}}\,n\pi \xi }}{{\omega _{sn} \{ {\left( {n\pi } \right)[ {1 + \mu _s^2 \left( {n\pi } \right)^2 } ] - \eta _g^2 \lambda ^2 } \}}} + } } \right. \\ \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\left. { + \mu _s^2 \sum\limits_{n = 1}^\infty {\frac{{\left( {n\pi } \right)^2 {\rm{sin}}\,\omega _{sn} T\,{\rm{sin}}\,n\pi \xi }}{{\omega _{sn} \{ {\left( {n\pi } \right)[ {1 + \mu _s^2 \left( {n\pi } \right)^2 } ] - \eta _g^2 \lambda ^2 } \}}}} } \right\}. \\ \end{array}

Let us assume as was the case in Section 3 that the beam is loaded in a static way with point ξ₀. In this case, the displacement of the beam and the bending moment have the following form: (65)w(ξ,ξ0)=2PL3EI∑n=1∞sin nπξ0 sin nπξ(nπ)4[1+μs2(nπ)2],w(\xi ,\xi _0 ) = 2\frac{{PL^3 }}{{EI}}\sum\limits_{n = 1}^\infty {\frac{{{\rm{sin}}\,n\pi \xi _0 {\rm{sin}}\,n\pi \xi }}{{\left( {n\pi } \right)^4 [ {1 + \mu _s^2 \left( {n\pi } \right)^2 } ]}},} (66)M(ξ,ξ0)=2PL{∑n=1∞sin⁡nπξ0sin⁡nπξ(nπ)2[1+μS2(nπ)2]++ μS2∑n=1∞(nπ)2sin⁡nπξ0sin⁡nπξ(nπ)2[1+μs2(nπ)2]}.\begin{array}{l} M (\xi ,\xi _0 ) = 2PL\left\{ {\sum\limits_{n = 1}^\infty {\frac{{{\rm{sin}}\,n\pi \xi _0 {\rm{sin}}\,n\pi \xi }}{{\left( {n\pi } \right)^2 [ {1 + \mu _s^2 \left( {n\pi } \right)^2 } ]}} + } } \right. \\ \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\left. { + \mu _S^2 \sum\limits_{n = 1}^\infty {\frac{{\left( {n\pi } \right)^2 {\rm{sin}}\,n\pi \xi _0 {\rm{sin}}\,n\pi \xi }}{{\left( {n\pi } \right)^2 [ {1 + \mu _S^2 \left( {n\pi } \right)^2 } ]}}} } \right\}. \\ \end{array}

By comparing the solutions presenting the form of the series for two models, it is easy to notice that, in the strain gradient model, the series converge faster than that in case of the stress gradient model.

The deflection and the bending moment for the two models of beams considered in the article, namely, the Eringen nonlocal model and the strain gradient model, were compared. The analysis was carried in the dimensionless variables for different parameters η_g, which describe the ratio of the force velocity to wave velocity, and also for different parameters μ_e, μ_s, which present nonlocal material properties. The slenderness of the beam is assumed to be equal to λ = 30. The displacement and the bending moment of the beam are shown in Figures 2–7 when the force is in the center of the span of the beam (T = 0.5). Figures 8–13 depict vibrations of the beam in the middle cross section, in the case of when the beam is excited by a moving force. In the figures, you can see various forms of the beam’s vibrations. The form of these vibrations depends on whether the force moves at a velocity smaller or greater than the critical velocity. Using Table 1, we can determine whether the force’s velocity is less (η < η_e/s,cr) than or greater (η > η_e/s,cr) than the critical velocity. Notice that when the force velocity is greater than the critical velocity, the displacement of the beam occurs only because of the moving force (Figs. 4, 6, 10, and 12). As can be seen from Table 1, the critical velocities for the stress gradient model are clearly smaller than those for the strain gradient model.

Figure 2

Figure 3

Figure 4

Figure 5

Figure 6

Figure 7

Figure 8

Figure 9

Figure 10

Figure 11

Figure 12

Figure 13

Figure 2 shows the displacement and the bending moment for the stress gradient model. As shown in Table 1, the force is moving at a velocity smaller than the critical velocity.

Figure 3 shows the displacement and the bending moment for the strain gradient model. As shown in Table 1, the force is moving at a velocity smaller than the critical velocity.

Figure 4 shows the displacement and the bending moment for the stress gradient model. In the case of μ_e = 0.2, the force is moving at a velocity smaller than the critical velocity; if μ_e = 0.4, the velocity of the force is slightly smaller than the critical speed; and when μ_e = 0.6, it is slightly higher than the critical velocity. Let us notice that the bending moment changes the sign if μ_e = 0.4 or μ_e = 0.6.

Figure 5 shows the displacement and the bending moment for the strain gradient model. As shown in Table 1, the force is moving at a velocity smaller than the critical velocity.

Figure 6 shows the displacement and the bending moment for the stress gradient model. In all cases, the force is moving at a velocity higher than the critical velocity. Let us notice that the bending moment changes the sign. In the case μ_e = 0.4 or μ_e = 0.6, the beam displacement occurs only after the moving force.

Figure 7 shows the displacement and the bending moment for the strain gradient model. As shown in Table 1, the force is moving at a velocity smaller than the critical velocity.

As shown in Table 1, the force is moving at a velocity smaller than the critical velocity.

As shown in Table 1, the force is moving at a velocity smaller than the critical velocity.

In the case μ_e = 0.2, the force is moving at a velocity smaller than the critical velocity; if μ_e = 0.4, the velocity of the force is slightly higher than the critical speed, and when μ_e = 0.6, it is higher than the critical velocity. In the last case (for μ_e = 0.6), the beam displacement in the middle of the span occurs only after the force passes this point.

As shown in Table 1, the force is moving at a velocity smaller than the critical velocity.

As shown in Table 1, the force is moving at a velocity higher than the critical velocity.

As shown in Table 1, the force is moving at a velocity smaller than the critical velocity.

We considered the dynamic response of a finite, simply supported Euler–Bernoulli beam loaded by a force moving with a constant velocity based on two different gradient nonlocal elasticity models. Both in the Eringen model, in which the strain is a function of gradient stresses, and in the nonlocal model, in which the stress is a function of gradient strain, the theoretical solution for the beam displacements and the bending moments were obtained in the form of the sine series. For Eringen’s model, the particular integral was also obtained in a closed form. The convergence of the series in the gradient stress model (Eringen’s model) is much weaker than that in the gradient strain model. This is due to the fact that Eringen’s natural frequencies are lower than those of the gradient strain model. It means that the stiffness of the beam in the stress gradient model is lower than that in the strain gradient model. Moreover, the critical velocity in the gradient stress model is smaller than that in the gradient strain model. In the case of a stress gradient model, it was a bit difficult to determine the bending moment. This is due to the fact that it cannot be determined directly from formula (4) and also from the fact that the solution determining the bending moments is poorly convergent. These difficulties were overcome by solving the system of two differential equations (14) and (15) instead of Eq. (16) and using the closed solution for aperiodic vibrations. There is also a second reason for which the solution was determined from the system of Eqs. (14) and (15) and not the fourth-order differential equation (16). This approach to the problem allowed to present a solution for both displacement and bending moment in the form of the sine Fourier series. The boundary condition (11) requires some comment (see also [38]). Taking into account the dependencies (4) and (8), it follows that for a simply supported beam at both its ends, the displacement gradient is equal to the bending moment gradient and does not have to be zero in the general loading case. In a special case, when the beam is subjected to static action, concentrated force condition (11) is met. Therefore, assuming boundary condition (11) in the problem under consideration has some justification. The presented closed form solutions have an important meaning in the case when we consider the bending moment or shear force in the beam, particularly in the vicinity of the load point. The differential equation in the gradient strain model, which describes the vibration of the Euler–Bernoulli beam, is of the sixth order. For this reason, we need six boundary conditions instead of four, and it is more difficult to obtain the closed solution for a particular integral than for the stress gradient model and local elasticity model. It is worth emphasizing the significant difference between the two considered beam models. In the nonlocal stress model, the given strain at a given point generates the stresses with nonlocal distribution. In the second model, the stress at a given point generates the nonlocal deformations. In the numerical analysis, dimensionless parameters were used, such as the ratio of the force velocity to the wave propagation velocity, or the ratio of the parameter determining the non-local material properties to the beam span.

It is worth mentioning that in the stress gradient model in the beam, displacement curve discontinuity appears in the point of the concentrated load. This is due to the fact that in Eq. (16) apart from the Dirac delta function, the gradient of this delta also occurs. The second derivative of the Dirac delta function in the Euler–Bernoulli beam models is the angular dislocation of the displacement. It is the explanation for the discontinuity of the curve. For this model, the form of the closed solution depends on whether the speed of movement of the force is less than a certain limit value (v<vgμeλv < \frac{{v_g }}{{\mu _e \lambda }}) or greater (v>vgμeλv > \frac{{v_g }}{{\mu _e \lambda }}). A similar situation occurs in the case of string vibrations [42], the Timoshenko beam [41], and the sandwich beam [43], where the solution depends on whether the speed of force is lower or higher than the speed of propagation of the wave. Such a wave effect does not occur for the local elasticity Euler–Bernoulli beam model [41].

The obtained results for the analyzed gradient models differ in both quantitative and qualitative terms. As shown in [44], a similar situation occurs in the case of stationary stochastic vibrations of the beam. For stochastic excitation in the case of the stress gradient model, the nonlocal effect also depends on the load distribution along the beam. It is important whether the load is distributed over the entire length of the beam or is in the form of a concentrated force. As can be seen, the solutions in both models depend not only on the mechanical and structural properties of the materials but also on the dimensions of the beam (slenderness of the beam).

To sum up, the most important original elements of the presented work are: –

In the Eringen model, beam displacements (vibrations) and bending moments were determined from the system of equations (14) and (15) in the form of the sum of series rather than fourth-order differential equation, avoiding ambiguity associated with boundary conditions (11), and closed solutions were determined allowing to take into account the peculiarities of bending moments diagrams (“peak” in the place of concentrated force) and shear force (jump in the same place),

Critical velocities of the moving force were determined for both models and numerical analysis was carried out, emphasizing whether the force velocity is lower or higher than the critical velocity.

The presented problem of vibrations of the Euler–Bernoulli beam caused by a moving force can be generalized by adopting more complex gradient models, for example, by combining the two models adopted in the work or by adopting a higher-order gradient model (see also [37, 38, 39]).