Nội dung

Vietnam Journal of Mechanics, VAST, Vol. 41, No. 1 (2019), pp. 63 – 78 DOI: https://doi.org/10.15625/0866-7136/12781 VIBRATION UNDER VARIABLE MAGNITUDE MOVING DISTRIBUTED MASSES OF NON-UNIFORM BERNOULLI–EULER BEAM RESTING ON PASTERNAK ELASTIC FOUNDATION T. O. Awodola∗ , S. A. Jimoh, B. B. Awe Federal University of Technology, Akure, Nigeria ∗ E-mail: toawodola@futa.edu.ng Received: 15 July 2018 / Published online: 29 October 2018 Abstract. The dynamic response to variable magnitude moving distributed masses of simply supported non-uniform Bernoulli–Euler beam resting on Pasternak elastic foundation is investigated in this paper. The problem is governed by fourth order partial differential equation with variable and singular coefficients. The main objective of this work is to obtain closed form solution to this class of dynamical problem. In order to obtain the solution, a technique based on the method of Galerkin with the series representation of Heaviside function is first used to reduce the equation to second order ordinary differential equations with variable coefficients. Thereafter the transformed equations are simplified using (i) The Laplace transformation technique in conjunction with convolution theory to obtain the solution for moving force problem and (ii) finite element analysis in conjunction with Newmark method to solve the analytically unsolvable moving mass problem because of the harmonic nature of the moving load. The finite element method is first used to solve the moving force problem and the solution is compared with the analytical solution of the moving force problem in order to validate the accuracy of the finite element method in solving the analytically unsolvable moving mass problem. The numerical solution using the finite element method is shown to compare favorably with the analytical solution of the moving force problem. The displacement response for moving distributed force and moving distributed mass models for the dynamical problem are calculated for various time t and presented in plotted curves. Keywords: moving mass; finite element; Newmark method; Pasternak elastic foundation; Galerkin’s method; resonance. 1. INTRODUCTION Force vibration of elastic bodies (stretched string, spring mass system, rods, etc.) have been extensively studied by several authors [1–11]. The vibrations may be due to (i) a force (load) which is a function of the space coordinates only or (ii) a force which c 2019 Vietnam Academy of Science and Technology 64 T. O. Awodola, S. A. Jimoh, B. B. Awe varies in both space and time. Such forces can either be of constant magnitude or variable magnitude. The present work concerns the effects of a force of variable magnitude moving at constant speed on an elastic body. In particular, the elastic body under consideration is the beam. It should, however, be mentioned from the onset, that such an elastic body (long and thin or stubby) is normally considered as a one-dimensional body [9–11] whose physical properties (stiffness, mass, length) are described with reference to a single dimension, the position along the elastic axis. Consequently, the partial differential equation describing the motion of such an elastic body is made up of only two independent variables, distance along the axis and time. Limiting the discussions to that of a moving force on a finite beam, Timoshenko [2], Inglis [12] and Muscolino [13] considered the problem of transverse oscillations of a beam subjected to a harmonic force moving with a uniform velocity. They assumed that the beam is simply supported. An analysis of the effect of such a moving force on the beam is given. Recently, Steele [14] investigated the effect of this moving force on beams to a unit force moving at a uniform velocity. The effect of this moving force on beams with and without an elastic foundation are analyzed. Some of the previous works involving non-uniform beams include that of Wu [15] studied the dynamic responses of multi-span non-uniform beams under moving loads using the transfer matrix method. Dogush [16] also studied dynamic behavior of multispan non-uniform beams traversed by a moving load at constant and variable velocities using both modal analysis and direct methods. Ahmadian et al. [17] investigated the analysis of a variable cross-section beam subjected to a moving concentrated force and mass using finite element method. However, the above research works on both uniform and non-uniform beams are impactful, but moving loads have been idealized as moving concentrated loads which act at certain points on the structure and along a single line in space. That is, the moving load is modelled as a lumped load. In practice, it is known that loads are actually distributed over a small segment or over the entire length of the structural member as they traverse the structure. Such moving loads are termed uniform distributed loads. Concentrated forces are mere mathematical idealization, which cannot be found in the real world, where surface forces act over an area and body forces act within volume. We also remark at this juncture, that only long thin uniform beam (called Euler’s beam) resting on one parametric foundation or bi-parametric foundation that is not harmonic in nature were considered. Thus, the present investigation is concerned with the vibration under variable magnitude moving distributed masses of non-uniform Bernoulli–Euler beam resting on Pasternak elastic foundation. The vital aspects of inertia terms are considered. Specifically, the elastic properties of the beam such as the flexural rigidity and the mass per unit length of the beam are assumed not constant. That is, the beam is of non-uniform cross-section and mass contains negligible damping. 2. THEORY AND FORMULATION OF THE PROBLEM In this study, the problem of a non-uniform Bernoulli–Euler beam and carrying a mass M is investigated. The beam’s properties such as moment of inertia I and the mass per unit length µ of the beam remained changing along the span length L. The transverse displacement V ( x, t) of the beam travelling at a uniform velocity as shown in Fig. 1 is Vibration under variable magnitude moving distributed masses of non-uniform Bernoulli–Euler beam resting on. . . 65 Fig. 1. Geometry of a variable moving masses of non-uniform Bernoulli–Euler beam resting on Pasternak foundation given as   ∂2 ∂2 V ( x, t) ∂2 V ( x, t) ∂2 V ( x, t) ∂2 V ( x, t) EI ( x ) + µ( x ) − N0 + K0 V ( x, t)− G0 = P( x, t), 2 2 2 2 ∂x ∂x ∂t ∂x ∂x2 (1) where t is the time coordinate, µ( x ) is the variable mass per unit length of the beam, EI ( x ) is the variable flexural stiffness, x is the spatial coordinate, K0 is the foundation stiffness, G0 is the shear modulus, N0 is the axial force and P( x, t) is the uniform distributed load acting on the beam. In this problem, the distributed load moving on the beam under consideration has mass commensurable with the mass of the beam. Consequently, the load inertia is not negligible but significantly affects the behaviour of the dynamical system. Thus, the distributed load P( x, t) takes the form   j 1 d2 V ( x, t) , P( x, t) = cos(ωt) ∑ Mi gH [ x − f (t)] 1 − g dt2 i =1 (2)  2 2 d2 ∂2 d f ( t ) ∂2 d2 f ( t ) ∂ d f (t) ∂ = 2 +2 + + , dt2 ∂t dt ∂x∂t dt ∂x2 dt2 ∂x where g denotes the acceleration due to gravity, d2 is a convective acceleration operadt2 ∂2 is the support beam’s acceleration at the point of contact with the moving mass, ∂t2   d f ( t ) ∂2 d f ( t ) 2 ∂2 is the well-known Coriolis acceleration, is the centripetal acceldt ∂x∂t dt ∂x2 d2 f ( t ) ∂ eration of the moving mass and is the acceleration component in the vertical dt2 ∂x direction when the moving load is not constant. In the same vein, for constant velocity c the direction/distance travelled by the load on the beam at any given instance of the time t is given as tor, f (t) = ci t, where x0 represent the point of application of force P( x, t) at any instant time t = 0. (3) 66 T. O. Awodola, S. A. Jimoh, B. B. Awe Moreover, the moving load is assumed to be of mass, M and time t is assumed to be limited to that interval of time within which the mass M is on the beam. i.e. 0 ≤ f (t) ≤ L, (4) and H [ x − f (t)] is the Heaviside function, which is a typical engineering function made to measure engineering application involving function that are either “on” or “off” and it is defined as (  1, x ≥ f (t). 1, x > ct. H (x) = H [ x − f (t)] = (5) 0, x ≤ ct. 0, x < f (t). As an example, let the variable moment of inertia I and the mass per unit length of the beam be defined, respectively, as [18] 3    πx πx , µ( x ) = µ0 1 + sin , (6) I ( x ) = I0 1 + sin L L where I0 and µ0 are constant moment of inertia and constant mass per unit length of the corresponding uniform beam respectively. To this end, substituting Eqs. (2), (3) and (6) into (1), after some simplification and rearrangement yields    EI0 2πx πx 3πx ∂4 V ( x, t) 6πEI0 πx 2πx 10 − 6 cos + 15 sin − sin + 5 cos + 4 sin 4 4 L L L 4L L L ∂x  3   2 2 3πx ∂ V ( x, t) 3π EI0 3πx 2πx πx ∂ V ( x, t) − cos + 3 sin + 8 cos − 5 sin L L L L ∂x3 4L2 ∂x2   2 πx ∂ V ( x, t) ∂2 V ( x, t) ∂2 V ( x, t) + µ0 1 + sin − N + K V ( x, t ) − G 0 0 0 L ∂t2 ∂x2 ∂x2   j j 2 ∂2 V ( x, t) ∂2 V ( x, t) 2 ∂ V ( x, t ) + 2c + cos(ωt) ∑ Mi H ( x − ci t) + c = Mi g cos ωtH ( x − ci t). i ∑ i ∂x∂t ∂t2 ∂x2 i =1 i =1 (7) The boundary conditions of the above problem are assumed to be arbitrary, that is, it can take any form of the classical boundary conditions. The initial conditions however without any loss of generality is given by ∂V ( x, 0) = 0. ∂t Eq. (7) forms the fundamental equation of the dynamic problem. V ( x, 0) = (8) 2.1. Solution procedure Eq. (7) is a non-homogeneous partial differential equation with variable coefficients. Evidently, the method of separation of variables is inapplicable as a difficulty arises in getting separate equations whose functions are function of a single variable. Thus, we resort to a modification of the approximate method best suited for solving diverse problem in dynamics of structures generally referred to as Galerkin’s Method. Therefore, we use the Galerkin’s method described in Oni and Awodola [19, 20] to reduce the fourth order Vibration under variable magnitude moving distributed masses of non-uniform Bernoulli–Euler beam resting on. . . 67 partial differential equation to a sequence of second order ordinary differential equation. Thus, yield a solution of the form n V ( x, t) = ∑ Ym (t)Um (x), (9) m =1 where Um ( x ) is chosen as a suitable kernel of the Galerkin’s method in (9) such that the boundary conditions given are satisfied. It is remarked at this juncture that the beam under consideration is assumed to have general boundary conditions at the edges x = 0 and x = L, therefore, we choose an appropriate selection of function for beam problem i.e beam shapes. Thus, the mth normal mode of vibration of the beam Um ( x ) = sin λm x λm x λm x λm x + Am cos + Bm sinh + Cm cosh , L L L L (10) is chosen such that the boundary conditions are satisfied. The kernel is chosen as Uk ( x ) = sin λk x λ x λ x λ x + Ak cos k + Bk sinh k + Ck cosh k , L L L L (11) where in (10) and (11), λm and λk are the mode frequency. Am , Bm , Cm , Ak , Bk and Ck are constants which are obtained by substituting (6) and (7) into appropriate boundary condition. Therefore, substituting Eqs. (9) into (7), yields n EI0  2πx πx 3πx  00 πx  Um ( x )Ÿm (t) + Um ( x )Ym (t) 10 − 6 cos + 15 sin − sin L 4µ0 L L L m =1 6πEI0  3πx πx 2πx 3πx  000 3π 2 EI0  2πx + 3 sin 5 cos + 4 sin − cos Um ( x )Ym (t) + 8 cos 4µ0 L L L L L L 4µL2  00 N0 00 πx K0 G0 00 Um ( x )Ym (t) − − 5 sin U ( x )Ym (t) + Um ( x )Ym (t) − U ( x )Ym (t) L µ m µ µ m ∑ h 1 + sin j  0 + ∑ Mi cos ωt H ( x − ci t)Um ( x )Ÿm (t) + 2ci H ( x − ci t)Um ( x )Ẏm (t) i =1 j i  00 + c2i H ( x − ci t)Um ( x ) − ∑ Mi g cos ωtH ( x − ci t) = 0. i =1 (12) In order to determine an expression for Ym (t), we shall consider a mass M travelling with uniform velocity c along the x-coordinate. The solution for any arbitrary umber of moving masses can be obtained by superposition of the individual solution since the governing equation is linear. Therefore, for the single mass M1 , it is required that the expression on the left hand side of Eq. (12) is orthogonal to the function Uk ( x ). Thus, using Eqs. (10) and (11) in (12), yields   cos ωt g cos ωt ∗ ∗ ∗ ∗ 2 ∗ I0 Ÿm (t) + I1 Ym (t) + M I2 Ÿm (t) + 2cI3 Ẏm (t) + c I4 Ym (t) = MI50 , (13) µ0 µ0 68 T. O. Awodola, S. A. Jimoh, B. B. Awe where I0∗ n ∑ = m =1 0 I1C n EI0 4µ0 1 + sin Z L πx  Um ( x )Uk ( x )dx, L I1∗ = I1A + I1B + I1C − I1D + I1E − IF , (14) πx 3πx  iv 2πx + 15 sin − sin Um ( x )Uk ( x )dx, ∑ L L L m =1 0 Z L πx 2πx 3πx  000 6πEI0 n + 4 sin − cos = Um ( x )Uk ( x )dx, 5 cos 4µ0 L m∑ L L L =1 0 Z L 3πx 3π 2 EI0 n 2πx πx  00 = 3 sin + 8 cos − 5 sin Um Uk ( x )dx, 4µ0 L2 m∑ L L L =1 0 I1A = I1B Z L I1D = n Z L N µ0 m =1 0 n Z L ∑ 10 − 6 cos 00 Um Uk ( x )dx, I1E = K0 µ0 n ∑ Z L m =1 0 Um Uk ( x )dx, I1F = G0 µ0 n ∑ Z L m =1 0 (15) (16) (17) 00 Um Uk ( x )dx, (18) I2∗ = I4∗ = ∑ m =1 0 n Z L ∑ m =1 0 H ( x − ct)Um Uk ( x )dx, 00 H ( x − ct)Um Uk ( x )dx, I3∗ = I50 = n ∑ Z L m =1 0 Z L 0 0 H ( x − ct)Um Uk ( x )dx, H ( x − ct)Uk ( x )dx. (19) (20) Using the property of Heaviside function, it can be expressed in series form given by [13] i.e. ∞ sin(2n + 1)πx cos(2n + 1)πct 1 − 2n + 1 π n =0 ∞ cos(2n + 1)πx sin(2n + 1)πct . 2n + 1 n =0 (21) Thus, in view of (14)–(20) and (21), it can be shown that 1 1 H ( x − ct) = + 4 π ∑ I ∗ (m, k ) ε 0 cos ωt Ÿm (t) + 1∗ Ym (t) + ∗ I0 (m, k ) I0 (m, k ) Lψ1A (m, k ) + L π ∞ cos(2n + 1)πct ∗ I5 (m, k ) 2n + 1 n =0 ∑   sin(2n + 1)πct ∗ L ∞ cos(2n + 1)πct ∗ × I6 (m, k) Ÿm (t)+ 2c Lψ2A (m, k)+ ∑ I7 (m, k ) ∑ 2n + 1 π n =0 2n + 1 n =0 "  L ∞ sin(2n + 1)πct ∗ L ∞ cos(2n + 1)πct ∗ I8 (m, k ) Ẏm (t) + c2 Lψ3A (m, k ) + ∑ I9 (m, k ) − ∑ π n =0 2n + 1 π n =0 2n + 1 #)  L ∞ sin(2n + 1)πct ∗ MgL cos ωt h − ∑ I10 (m, k ) Ym (t) = × − cos λk x + Ak sin λk x π n =0 2n + 1 µλk I0∗ (m, k ) λ ct λ ct λ ct λ ct i + Bk cosh λk x + Ck sinh λk x + cos k − Ak sin k − Bk cosh k − Ck sinh k , L L L L (22) L − π ∞ ( ∑ Vibration under variable magnitude moving distributed masses of non-uniform Bernoulli–Euler beam resting on. . . 69 where M . (23) µL Eqs. (22) is now the fundamental equation governing the dynamic problem. This coupled non-homogeneous second order ordinary differential equation holds for all variants of classical boundary conditions. It follows that two special cases of Eq. (22) arise, namely, the moving force and moving mass problems. ε0 = 2.2. Non-uniform Bernoulli–Euler beam traversed by moving distributed force for simply supported end condition In this section, an approximate model of the differential equation describing the response of the elastic structure is obtained by neglecting inertia terms. i.e. Γ0 = 0 and we shall limit our example on simply supported end condition. In this case, the displacement and the bending moment vanish. Thus ∂2 Vm ( L, t) ∂2 Vm (0, t) = 0 = , ∂x2 ∂x2 (24) Um (0) = 0 = Um ( L), ∂2 Um (0) ∂2 Um ( L) = 0 = , ∂x2 ∂x2 (25) Uk (0) = 0 = Uk ( L), ∂2 Uk (0) ∂2 Uk ( L) = 0 = . ∂x2 ∂x2 (26) Vm (0, t) = 0 = Vm ( L, t), and hence for normal modes which implies that Applying (25) and (26) on (10), we have Am = Ak = Bm = Bk = Cm = Ck = 0, λm = mπ, λk = kπ, (27) as the mode frequencies, and mπx kπx , Uk ( x ) = sin , L L as the mode functions. Using (24)–(28) in (22), yields h i Ÿm (t) + ω 2f Ym (t) = Pk cos ωt cos θk t + Rk , Um ( x ) = sin (28) (29) where ω 2f = I1∗ (m, k ) , I0∗ (m, k ) Pk = MgL , µkπ I0∗ (m, k ) θk = kπc , L Rk = −(−1)k , −4mkL , m 6= k π [1 − (m + k )2 ][1 − (m − k )2 ] ∗ I0 (m, k ) = L 4m2 L   , m≡k  − 2 π (1 − 4m2 ) I1∗ (m, k ) = I1A − I1B − I1C + I1D + I1E + I1F , (30)     (31) 70 T. O. Awodola, S. A. Jimoh, B. B. Awe I1A I1B I1C  4 4 h i −60mkL 12mkL   m π EI0 , +  4µ0 L4 π [1 − (m + k)2 ][1 − (m − k)2 ] π [9 − (m + k)2 ][9 − (m − k)2 ] = h i 4 4 4m2 L 4m2 L   m π EI0 5L −  , + 2 2 π (1 − 4m ) 3π (9 − 4m ) 4µ0 L4  3 3 h
i m+k m+k m−k m−k   3m π EI0 − − 5 + ,  8µ0 L3 9 − ( m + k )2 9 − ( m − k )2 1 − ( m + k )2 1 − ( m − k )2 = h 2m 3 3 10m i   3m π EI0  − , 8µ0 L3 9 − 4m2 1 − 4m2  2 4 h i −36mkL 20mkL   3m π EI0 , +  2 2 2 2 4 π [9 − (m + k) ][9 − (m − k) ] π [1 − (m + k) ][1 − (m − k) ] 4µ0 L = h −4mkL 2 4 20m2 L i   3m π EI0  + , π (9 − 4m2 ) π (1 − 4m2 ) 4µ0 L4 m 6= k m ≡ k. (32) m 6= k m ≡ k. (33) m 6= k m ≡ k. (34) m2 π 2 N K0 L m2 π 2 G0 , I1E = , I1F = . (35) 2µ0 L 2µ0 2µ0 L Therefore, solving (29) using the method of Laplace transforms and convolution with the initial conditions (8), we have I1D = n Pk V ( x, t) = ∑ 2ω f m =1 ( ωf (ω 2f − ω 2 )(ω 2f − Ω21 )(ω 2f − Ω22 ) h 2Rk (ω 2f − Ω21 )(ω 2f − Ω22 )(cos ωt − cos ω f t) i mπx . +(ω 2f − ω 2 )(ω 2f − Ω22 )(cos Ω1 t − cos ω f t) + (ω 2f − ω 2 )(ω 2f − Ω21 )(cos Ω2 t − cos ω f t) sin L (36) (36) above represents the transverse displacement response to a distributed force moving at constant velocity of a non-uniform simply supported Bernoulli–Euler beam resting on Pasternak elastic foundation. 2.3. Non-uniform Bernoulli–Euler beam traversed by moving distributed mass In this section, the solution to the entire equation (22) is sought when no terms of the coupled differential equation is neglected. Evidently, an exact solution to this equation is not possible. All conventional methods break down, even the popular Struble’s technique [21] could not handle it because of the variability of the magnitude of the moving load. Hence we resort to using finite element method (FEM) to model the structure. Then Newmark numerical method of integration is used to solve the resulting semi-discrete time-dependent equation to obtain the desired responses. 2.3.1. Finite Element Method (FEM) The finite element technique assumes that the unknown transverse deflection of the non-uniform beam, V ( x, t), can be represented approximately by a set of piecewise continuous functions which are defined over a finite number of sub-regions called elements Vibration under variable magnitude moving distributed masses of non-uniform Bernoulli–Euler beam resting on. . . 71 and composed of the numerical values of the unknown deflection within the region. Thus, the first step involved in the technique, consist of dividing the spatial solution domain of the non-uniform beam, which happened to be the length of the beam in this case, into a number of sub-domains known as finite elements. These elements are joined to each other at selected points called nodes. Next, the weak or variational form of the governing equation (1) is constructed as follows: Consider a typical element of length L so that its domain λe = (0, L). Substituting (2) and (3) into (1), we have   ∂2 V ( x, t) ∂2 V ( x, t) ∂2 ∂2 V ( x, t) ∂2 V ( x, t) − N + K V ( x, t ) − G EI ( x ) + µ ( x ) 0 0 0 ∂x2 ∂x2 ∂t2 ∂x2 ∂x2    2 j 2 ∂2 V ( x, t) ∂ V ( x, t) 2 ∂ V ( x, t ) + 2ci . = cos(ωt) ∑ Mi H ( x − ci t) g − + ci ∂t2 ∂x∂t ∂x2 i =1 (37) In order to solve Eq. (37), we shall consider a mass M travelling with uniform velocity c along the x-coordinate. The solution for any arbitrary number of moving masses can be obtained by superposition of the individual solution since the governing equation is linear. Therefore, for the single mass M1 , let W ( x ) be Galerkin’s weight function. Multiplying Eq. (37) by the weight function and integrate over the domain λe , after some simplification and rearrangement yields Z Le ∂2 V ( x, t) ∂2 W ( x ) dx + ∂x2 ∂x2 Z Le µ( x ) ∂2 V ( x, t) W ( x )dx − ∂t2 Z Le ∂2 V ( x, t) W ( x )dx ∂x2 0 0 0 Z Le Z Le Z Le ∂2 V ( x, t) + K0 V ( x, t)W ( x )dx − Mg cos(ωt) H ( x − ct)W ( x )dx + M cos(ωt) H ( x − ct) W ( x )dx ∂t2 0 0 0 Z Le Z Le ∂2 V ( x, t) ∂2 V ( x, t) + 2Mc cos(ωt) W ( x )dx + Mc2 cos(ωt) W ( x )dx H ( x − ct) H ( x − ct) ∂x∂t ∂x2 0 0 ∂W ∂W + W ( Le ) B3e − W (0) B1e − B4e + Be = 0, e ∂x x= L ∂x x=0 2 (38) EI ( x ) ( N0 + G0 ) where    2  h i ∂ ∂2 V ( x, t) ∂ V ( x, t) e λ= EI ( x ) , φ = EI ( x ) , B = λW ( x ) k ∂x ∂x2 ∂x2 Le 0 h ∂W ( x ) i − φ ∂x Le 0 , (39) λ is the shear force, φ is the bending moment and Bke , (k = 1, 2, . . . , 4) are the extremely important and necessary four boundary terms, two at each of the end nodes of the element. Furthermore, it can be readily shown that Z Le 0 H ( x − ct) f ( x )dx = Z Le ct f ( x ). (40) 72 T. O. Awodola, S. A. Jimoh, B. B. Awe Thus, Eq. (38) becomes Z Le 0 + EI ( x ) Z Le 0 ∂2 V ( x, t) ∂2 W ( x ) dx + ∂x2 ∂x2 Z Le 0 µ( x ) ∂2 V ( x, t) W ( x )dx − ∂t2 K0 V ( x, t)W ( x )dx − Mg cos(ωt) + 2Mc cos(ωt) Z Le 2 ∂ V ( x, t) ∂x∂t ∂W e e e + W ( L ) B3 − W (0) B1 − ∂x ct Z Le ct B4e + ∂W ∂x 0 x =0 Z Le Z Le 2 ∂ V ( x, t) ct ∂2 V ( x, t) W ( x )dx ∂x2 ∂2 V ( x, t) W ( x )dx ∂t2 (41) ( N0 + G0 ) W ( x )dx + M cos(ωt) W ( x )dx + Mc2 cos(ωt) x = Le Z Le ∂x2 ct W ( x )dx B2e = 0. Eq. (41) is the desired weak form of the variable magnitude moving distributed masses of non-uniform Bernoulli–Euler beam resting on elastic foundation. Therefore, we seek an approximate solution over the element under consideration and thereby construct the corresponding shape function. To this end, it is assumed that the unknown deflection V ( x, t) could be expressed approximately as V ( x, t) ≈ Vn ( x, t) = H1 ( x )V1 (t) + H2 ( x )V2 (t) + H3 ( x )V3 (t) + H4 ( x )V4 (t) 4 = ∑ Hk (x)Vk (t) = { H }{V (t)}, j = 1, 2, 3, 4 (42) k =1 where Hj ( x ) are called Hermite cubic shape functions and Vk (t) are the modal deflection functions and H is a row vector defined as   [ H ] = H1 ( x ), H2 ( x ), H3 ( x ), H4 ( x ) . (43) Using the procedures involved in constructing the Hermit-cubic interpolation functions in [22], yields 3x2 2x3 x2 x3 3x2 2x3 x2 x3 + , H = x − + , H = − , H = − + , (44) 2 3 4 h2 h3 h h2 h2 h3 h h2 where x is the spatial coordinate. Now substituting Eqs. (42)–(44) into the weak form (41), after some simplification and rearrangement gives    e      K V (t) + C e V̇ (t) + Me V̈ (t) + f e + Qe = 0. (45) H1 = 1 − The matrix equation (45) is the governing equation describing the behavior of a typical  finite element of the non-uniform beam traversed by a harmonic moving load. K e is the     element stiffness matrix, Me is the element mass matrix, C e is the element centripetal   matrix, f e is the force vector and Qe is the element boundary term vector. The next step is assembling of the element equations. The procedure for assembling various matrices and vectors for several beam elements which constitute a mesh is well discussed in [23, 24]. Hence the governing assembled equation of motion describing the dynamic behavior of the moving load problem with Pasternak foundation is        K V (t) + C V̇ (t) + M V̈ (t) = F , (46)