Finite element formulation and analysis of a functionally graded Timoshenko beam subjected to an accelerating mass including inertial effects of the mass

This study describes a new finite element method that can be used to analyse transverse and axial vibrations of a Functionally Graded Material (FGM) beam under an accelerating / decelerating mass. The differential equations of the FGM beam are obtained using First-order Shear Deformation Theory (FSDT). In these equations, the interaction terms of mass inertia are derived from the second-order exact differentiation of displacement functions with respect to mass contact point. The FGM beam is made of two different materials (Steel and Alumina Al 2 O 3 ), which vary in thickness with a power law. Including the effects of neutral axis shift and mass inertia, the proposed method can be used when the dynamic behaviour of the FGM Timoshenko beams is to be analysed in transverse and axial directions, depending on the interaction with the acceleration of the moving loads. After validating this work with literature studies, new investigations and findings are presented for both moving load and mass assumptions. In addition, the obtained results of Timoshenko Beam (TBT) and Euler Bernoulli beam theory (EBT) are compared for FGM beams with various speeds and accelerations of moving mass.


INTRODUCTION
In recent years, there has been a new class of composite materials known as Functionally Graded Materials (FGMs), where different components of materials are graded continuously under a force law. This nonhomogeneous composite new material is an alternative to stress concentrations in conventional composite structures where the laminates have very different mechanical and thermal properties. In FGM, the proportions of constituent materials are slowly changed volumetrically to remove the interfacial effects of thermal and mechanical stresses. In some engineering applications, materials are expected to withstand high temperatures under severe loading conditions. Because of better thermal resistance, the properties of the ceramic component can withstand high temperature environments, while the metal component in the functionally graded materials (FGMs) provides more mechanical strength to reduce fracture performance probability. In defense, aviation and aerospace, automotive, transportation systems and similar engineering applications, the FGMs have the ability to combine many different properties such as thermal and mechanical strength, structural lightness and wear resistance. It is expected that the use of FGMs in engineering applications will increase as the production opportunities with the developing technology increase. A wide range of FGM application systems are support systems such as beams. In recent years, many scientists have studied the static and vibrational behaviour of FGM beams (Aydogdu and Taskin, 2007;Chakraborty et al., 2003;Ebrahimi et al., 2016;Eltaher et al., 2013;Jafari-Talookolaei et al., 2018;Kang and Li, 2010;Khan et al., 2016;Kim and Lee, 2014;Kosmatka, 1995;Lien et al., 2017;Mashat et al., 2014;Masoodi and Moghaddam, 2015;Murín et al., 2010;Nguyen and Nguyen, 2015;Rezaiee-Pajand and Masoodi, 2016;Rezaiee-Pajand and Rajabzadeh-Safaei, 2018;Sina et al., 2009;Taati, 2018;Taeprasartsit, 2013;Yokoyama, 1996;Zhang and Zhou, 2008). Analysis of systems subjected to moving loads, an ongoing long-standing problem, and some studies have investigated the transverse responses of the limited number of FGM beams to moving forces, in the transverse and axial directions, neglecting the effects of inertia (Kadivar and Mohebpour, 1998;Nguyen et al., 2013;Rajabi et al., 2013;Şimşek and Al-shujairi, 2017;Şimşek and Kocatürk, 2009). Other studies of moving loads at constant and variable speeds have also investigated the transverse responses of beams and plates with uniform cross-section and homogeneous materials (Kadivar and Mohebpour, 1998;Michaltsos, 2002;Wang, 2009;Yoshida and Weaver, 1971 inertia effect of mass, some analytical and finite element studies of moving mass for beams and plates are given in (Cifuentes, 1989;Esen, 2013;Meirovitch, 1967;Wu, 2005).
The calculation methods given in the literature for FGM beams in general are for low speed load movements and simplified conditions where the convective acceleration effect of moving load mass is neglected. In this study, a more realistic modelling and analysis method of such systems is provided without neglecting the Coriolis, inertia and centrifugal effects of the moving load accelerating along the deflected shape of the FGM beams. The natural vibration frequencies are important in the dynamic behavior of FGM and other beams. For this reason, any factor that affects or alters the frequency of the beam in the dynamic interaction is also very important. This is also the reason why the inertia effect of mass should be involved in the dynamic interaction of moving mass and beam, and in this study, this behavior is examined in this context in order to better understand the dynamic behavior. In addition to transverse vibration, axial vibration is also induced by inertia effects of the mass of the load in the conditions of rapid acceleration and panic braking, which have been often neglected in the literature. Moreover, this work also presents a detailed comparison of moving load and moving mass assumptions for FGM beams, including differences between Timoshenko (TBT) and Euler-Bernoulli beam theories (EBT).
2 MATHEMATICAL MODELLING Figure 1 shows a FGM beam of length L exposed to a mass, moving from left to right with a variable speed. The modelling of the finite elements over the length is considered to have n beam elements with n + 1 nodes. The two nodes of each beam element have three displacements and forces at each node. Where xp is the time-dependent global position of the moving mass, and xm is the local coordinate of the mass measured from the left end of the element. A rectangular section with width b and height h and a material for the beam is considered to be a FGM consisting of two basic materials and it is assumed that effective material properties are distributed in the direction of thickness (z-direction) according to a power law as: In Eq. (1) P(z) signifies the effective material properties, including mass density, Young's modulus, shear modulus and Poisson ratios; Pt and Pb are the material properties of the material at the top and bottom surfaces, respectively; Vt and Vb respectively be the volume fractions of the materials at the top and bottom surfaces; n is the non-negative power-law index, expresses the distribution of the constituent materials. Figure 2 depicts the change in volume fraction along the beam thickness for various index values of the FGM beam. Figure 1: A thick FGM beam under the influence of a moving point mass.
Due to variation of the effective Young's modulus, the neutral axis is no longer at the mid-plane thus, the shift of the neutral axis can be decided by solving the following equation (Eltaher et al., 2013;Zhang and Zhou, 2008).
Where z0 is the geometric neutral axis shift and h is the thickness of the beam, the position of the physical neutral axis depends on the ratio between the Et / Eb of the Young modules and the power law index n of the material components. Figure 2 shows the effect of the power law index n on the neutral axis shift of the FGM beam for various values of the Et / Eb ratio. In the figure, the shift increases with increasing Et / Eb ratio, but decreases as n index increases. Here indices t and b, respectively, stand for top surface (ceramic) and bottom surface (metal). In the formulation, the following assumptions will be assumed ( Figure 1): inertia effects of mass will take part in the calculations, and the moving mass will always be in contact with the beam, separation not allowed. When a FGM beam is considered in a coordinate system (x, z) where the x-axis coincides with the physical neutral axis of the beam as shown in Fig. 1, by adopting the FSDT and by including the axial force at the contact point with the beam, the displacements at any point of the beam become: Where ( , ) u x t , ( , ) w x t and ( , ) x t  are the axial, transverse displacement and rotation of the cross section at point x on the physical neutral axis, Considering FSDT the strains and stresses related to the displacement field in Eq. (3) are given by Where  is the surface correction factor. From Eq. (4), the strain energy of the beam is given by Where D11 is the axial, D12 is the axial-rotational, D22 is the bending, and D33 is the shear rigidities, respectively and computed as The kinetic energy resulting from the equation (3) is Where The potential energy of the moving mass including convective acceleration is given by Where g is gravitational acceleration, δ(.) is the Dirac delta function, with x0 and v0 are initial position and velocity while a denotes the acceleration of the mass. Using Hamilton's principle, Eqs. (5), (7) and (9) lead to the resulting coupled governing differential equations and the geometric and natural boundary conditions for the hinged-hinged beam

Finite element formulation
When the static case of Eq. (11) solved, then using the derived shape functions u  , w  and   (given in Appendix A-A.1-A.4:), the nodal displacements of a two node Timoshenko beam element with six Dof are (Kosmatka, 1995;Yokoyama, 1996): are the displacements at node 2. With shape functions, the axial and transverse displacement and the rotation in the element are found from their nodal displacements using  are vectors of shape functions for u, w and θ, respectively Including geometric stiffness effect from axial force, and placing the shape functions in Eq. (14) into Eq. (6), the following total strain energy for all elements in equal length is obtained in the form where n symbolizes the total number of elements, and k is the element stiffness matrix, which consisting of the axial stiffness matrix k uu , coupling stiffness matrix k u , bending stiffness matrix k  , shear stiffness matrix k  , and geometric stiffness matrix G k , defined as follows The sing of the k G is negative for acceleration and positive for deceleration of the mass. The kinetic energy can be specified with the help of equation (14) as To include the effects of mass inertia, considering the convective acceleration the following contact force terms from Eq. (11) in axial and transverse directions should be evaluated. They are: These force terms effect on the beam element s on which the moving mass is located (Fig.1), at time t. Except for the axial force p m a and the gravitational transverse force p m g , when the velocity and acceleration of the mass is quite high, the inertia, Coriolis and centripetal force components can be high enough to effect of the interaction of the mass with the beam in transverse and axial directions. In such a case, for including these effects of the moving mass, Eq. (19) is rewritten using the same interpolating functions at the contact point of the mass and the nodal displacements of the Timoshenko beam element given in Eq. (14). Where, the position of the moving mass on the element s is  14), and then substituting into Eq. (19), finally rearranging the terms, the following time dependent second order differential matrix equation is obtained.
Where m , c and k are, respectively, the additional mass, damping and stiffness matrices, and f is the loading force vector of the beam element s (as given in Appendix A -A.6).

Equation of Motion of the entire System
The motion equation for the multi-degree of freedom damped structural system shown in Fig.1 is as follows Where M, C, and K are total mass, damping and stiffness matrices, respectively, q  , q  and q are acceleration, velocity, and displacement vectors, respectively. Also, F is the general external force vector at time t of the system. If there is a mass moving with a constant acceleration on the structure, mass, stiffness matrices of the all elements are obtained from Eqs.
(16 and 18))) and then the element matrices of the beam element (s) are modified by adding the property matrices given in Eq. (20). The instantaneous force vector F is also time dependent and the coefficients of the total force vector are zero except for the nodal forces of the beam element s. Thus, the overall force vector of the complete system becomes: effects of the mass

RESULTS AND DISCUSSION WITH NUMERICAL SOLUTIONS
The following non-dimensional frequency parameter is defined as In the literature, when dealing with the moving mass problem, the moving mass is generally accepted as the moving load and the inertia effects of the mass are neglected, and the fundamentals frequency parameter is calculated according to the above formula. According to this assumption, the frequency parameters are given according to the slenderness of the beams in Table 1 in comparison for the same materials of FGM beam made of Al and Al2O3 used in previous studies (NGUYEN et al., 2013;Sina et al., 2009) . However, it may be an acceptable approach if the velocity of the load is too small and the mass of the moving load is too small relative to the mass of the beam. Because, as shown in Figure 3, as the mass ratio increases, the frequency parameter is greatly influenced by the inertia effects of mass. Figure 3 shows the variation of the fundamental frequency parameter according to different power law exponent n and Et / Eb ratio and mass ratio (mass of the load / mass of the beam), and the position of the mass on the beam. As the mass ratio approaches 1, the fundamental frequency decreases by an average of 50% and, of course, this rate varies depending on the position of the mass on the beam. Figure 4 shows how the first four natural frequencies change relative to the position of the mass and mass ratio ε. The most important point to note here is the nodal points of the corresponding vibration mode at the 2nd, 3rd and 4th natural frequencies and there is no frequency change when the mass passes through these points. The shifts of physical neutral axis are given in Table 2.   The well-known velocity parameter for the fundamental mode is the ratio of the loading frequency (ω = πV / L) of the moving load to the fundamental natural frequency ω1 of the beam. If the velocity is constant, this term does not change, but in the case of the load motion with an acceleration, the velocity is not constant, so the loading frequency changes. When the load is considered as a moving force, the fundamental frequency is constant at ω1 for constant speeds. However, when the effect of mass inertia is included, the ω1 varies depending on the position and the mass of the load. (see Figure  3). In this work, the reason this issue is referred to and taken care of is that the frequency parameter is given only according to the material composition and the beam dimensions in most literature studies. Design and calculation should be done considering this fact in the FGM applications where the mass of moving load is high, and the travelling velocity is not constant.
Example 1. To compare the present method with the same assumptions as the others, we consider a simple supported isotropic beam-plate under m = 0.4485 kg, F = 4.4 N as a moving load investigated in the literature in both cases, adding the inertial effects of mass to the account and accepting the mass as a moving load. The size and material properties of the plate are the same as those preferred in Meirovitch 1967, i.e. lx = 10.36 cm; ly = 0.635 cm, h = 0.635 cm; E = 206.8 GPa, ρ = 10686.9 kg / m 3 . In Table 3, dynamic amplification factors (DAF), defined as the ratio of maximum dynamic displacement to maximum static displacement, are compared with numerical, analytical and experimental results using different theories available in the literature. For the Timoshenko beam (column 3) the results obtained with the present method appear to be close to the results of the analytical solution (Meirovitch, 1967) and the first order shear deformation theory (FSDT) (Michaltsos, 2002) . In addition, the current study and moving mass studies in the literature have been compared to show differences between moving mass and moving load assumptions, and the results of (Sharbati and Szyszkowski, 2011;Wu, 2005) and the present study are very close for moving mass assumption. However, unlike the moving load (motion force) approach which neglects the effects of mass inertia, the DAF has increased significantly at some high speeds, but at very low speeds, all results are very close. (1) Present moving load (TBT), (2) Present moving load (EBT), (3) Moving load, analytical from (Meirovitch, 1967), (4) Moving load, third order shear deformation theory (TSDT) from (Mohebpour et al., 2011), (5) Moving load (FSDT) from (Michaltsos, 2002), (6) Present moving mass (TBT), (7) Present moving mass (EBT), (8) Moving mass (EBT) from (Sharbati and Szyszkowski, 2011), (9) moving mass (EBT) from (Wu, 2005).
Example 2: let us consider the following FGM beam investigated by (Nguyen et al., 2013) and (Şimşek and Kocatürk, 2009), whose properties are shown in Table 4. Here, beam length L = 20 m, height h = 0.9 m, width 0.4 m and power lawexponent n = 0.2, 0.3, 0.5, 1, 2, 0, 1000. The calculated physical properties of the beam based on the power law-exponent are given in Table 5. It is assumed here that the beam is full steel on the base surface and full Al2O3 (99.9%) on the top surface. In columns 2 and 3, the effective Young's modulus's and densities are given according to power-law exponent. In column 4 the term Φ gives the relative significance of the shear deformations to the bending deformations, by setting Φ = 0 for neglecting shear deformation, displacement functions Eq. (A.3) are exact for beam element of EBT. The loading is 100 kN, and the corresponding moving mass, mp in column 7 while mass of the beam mb is given in 6.  Dynamic amplification factors (DAFs) of the maximum midpoint displacements versus velocity of the mass for different power-law exponents considering moving load and moving mass assumptions, are given in Figure 5. Where DAF Max (w(L/2, t))/ w0 is the ratio of the maximum displacement occurring at midpoint of the beam to the static displacement occurring when the mass is statically at the midpoint of the beam, where w0 =FL 3 /48EzIz. In Table 6, the maximum DAFs and corresponding mass velocities are given in the assumption of moving mass and moving load. There are small differences in DAFs compared to the given in (Şimşek and Kocatürk, 2009) even if the corresponding velocities are the same or very close. These small differences are due to the same ratio of Poisson's of alumina and steel, which is not clarified in (Şimşek and Kocatürk, 2009) and but they are accepted as the same in (Şimşek and Al-shujairi, 2017). The Poisson ratio of the alumina is v = 0.23 and it is 0.30 for the steel. Furthermore, the static deformation must be calculated considering the FGM beam, not the normal beam. Because, as seen in Table 5, the effective-Young's modulus, effective density, and beam mass vary according to the power-law exponent. This change affects both static and dynamic displacements. The other differences are the inclusions of axial coupling and rotary inertia. In this study, in addition to the shear effect, in the TBT (Timoshenko beam theory), the effects of the rotary inertia and axial coupling are fully applied. No effect has been neglected in this study, with the aim of setting up a basic reference to the future literature studies and researchers without neglecting the full effects. As can be seen from Fig. 5, the acceptance of moving mass leads to displacements which are twice as high as the moving load assumption. At small velocities below 25 m / s, both admissions yield close results. At high speeds, calculation of actual displacements with moving load assumption may give misleading results. As can be seen in Eqs. (A.6), the components of k vary linearly with the mass of the load while vary in proportion to the square of the velocity. At the same time, if the Coriolis effect is involved in the calculation, the components of the matrix c change linearly with both mass and velocity; and the effect of inertia is accounted for, the components matrix m vary depending on the mass. This negligence can be made at low speeds but must not be neglected at high and very high speeds.
Nowadays, especially in the case of high-speed train transport, the speed of transport vehicles has reached 500 km/h, as is the case of the Japanese bullet train. For moving load and moving mass assumptions; to compare the effects of low and very high speed, Fig. 6 shows the results obtained at velocities of v = 25, 50, 75, 100, 150 and 200 m / s. In tables 7 to 10, the normalized maximum displacement and corresponding mass positions in these graphs are presented.

Effect of the Sudden Acceleration and Deceleration
The effects of sudden acceleration on systems such as missile launchers and bridges are worth examining today as they may be in the case of sudden accelerations of missiles and projectiles or sudden braking of high-speed vehicles. To examine this effect and analyse the interaction of the beam with the mass, the results of sudden acceleration and deceleration situations are given in the following section. The same constant acceleration and deceleration rate of 600 m/s 2 are considered where the velocity of the mass is increased from zero to 154.86 m/s and decreased from this velocity to zero. In the case of acceleration, the mass is at the left end of the beam and the velocity is zero. In the case of braking, the velocity of the mass is zero at the right end of the beam. Figure 6 shows that normalized time dependent mid-point displacements of the FGM beam for v=25, 50, 75 and 150 m/s and for different power-law exponent n. Figure 7 shows the effect of constant acceleration on normalized mid-point displacements of different material constituents; and for EBT and TBT. In both theories, the largest displacements are occurred in the beams made of full steel material. The reason for this is that in the case of full steel, the mass of the beam is at its greatest value. Large mass causes large mass matrix m and high displacement due to low natural frequency in dynamic load. With the same geometry, the FGM compound, which is large in effective modulus of elasticity and low in effective density, makes the smallest displacement. In fact, this is one of the reasons why FGM materials are preferred. As can be seen from the figures, EBT has a lower displacement than TBT. One of the reasons for this is the greater the stiffness matrix in EBT and the stiffer the material is. Another reason is the effect of axial coupling and rotary inertia which does not considered in EBT. This is better understood from the following axial displacement graphs in Figure 8. For this reason, TBT modelling gives more accurate results in high-speed moving load applications. At low speeds and low dynamic loads, EBT can be applied instead of TBT as the results will be close to each other. We have chosen the example of high acceleration with the aim of revealing the difference between these two theories. In the case of acceleration and deceleration, considering the relative position of the mass on the beam, the contact point xp of the mass varies with a quadratic function of time while the velocity of the mass varies linearly, therefore, according to the relative position, it is expected that the displacement graphs of both situations will be very different. at rapid acceleration, the mass starts to accelerate from the left end of the beam, and before it reaches a quantity that will magnify the beam vibrations, the mass travels a great deal over the beam. Only after the mass passes the midpoint of the beam, the vibrational amplitudes grow. In the case of deceleration, the amplitudes increase earlier as the initial velocity and acceleration are higher. Another point is that the early excited beam has already begun to vibrate. The vibration amplitude tends to decrease with a fluctuation although the velocity of the mass rapidly decays over time. The decrease in amplitude is not as fast as the decrease in mass velocity. These events can be observed in both EBT and TBT models.
. These events can be observed in both EBT and TBT models.

CONCLUSIONS
In this study, a realistic model of a FGM Timoshenko beam under the influence of accelerated mass is established. The inertia effect of moving masses which are not adequately addressed in the literature studies for FGM beams is widely presented for different speeds and accelerations. It is understood that the dynamic behaviour of the beam is not overly affected by inertia effects at low and constant load speeds, as indicated in the literature for low speed civil engineering applications. However, at high speeds and in the case of acceleration or deceleration, it has been observed that the dynamic behaviour of the FGM beam changes significantly due to effect of the mass of the load. In some cases, the DAF increased by a factor of 2 due to the natural frequency change in the beam under the influence of the moving mass. In the analyses made, it was observed that the fundamental frequency changed by 50% (figures 5 and 6) when the mass of the moving load was equal to the mass of the beam. Thus, the decrease in fundamental frequency will reduce the speed parameter of the load. The speed parameter is the ratio of the influence frequency (ω=πv/L) of the moving mass to the actual fundamental frequency of the beam. The speed parameter of the fundamental frequency is defined as the ratio of ω/ω1, and the maximum response occurs when this ratio is equal to 1. It was also emphasized that design and calculation should be done considering the frequency variation in FGM applications where the mass and velocity of the load is high, and the speed of motion is not constant.
Modelling of the FGM beam with a moving mass have been extensively studied in the assumptions of TBT, EBT, moving mass and moving load at constant and variable moving speeds of the mass, and results of the analysis have been given in detail. The effects of sudden acceleration / deceleration cases on transverse and axial vibration of the FGM beam were also presented.