On the dynamic response of bi-directional functionally graded nanobeams under moving harmonic load accounting for surface effect

This paper presents an investigation of the dynamic behavior of bi-directionally functionally graded (BDFG) micro/nanobeams excited by a moving harmonic load. The formulation is established in the context of the surface elasticity theory and the modified couple stress theory to incorporate the effects of surface energy and microstructure, respectively. Based on the generalized elasticity theory and the parabolic shear deformation beam theory, the nonclassical governing equations of the problem are obtained using Lagrange’s equation accounting for the physical neutral plane concept. The material properties of the beam smoothly change along both the axial and thickness directions according to power-law distribution, accounting for the gradation of the material length scale parameter and the surface parameters, i.e., residual surface stress, two surface elastic constants, and surface mass density. Using trigonometric Ritz method (TRM), the trial functions denoting transverse, axial deflections, and rotation of the cross sections of the beam are expressed in sinusoidal form. Then, with the aid of Lagrange’s equation, the system of equations of motion are derived. Finally, Newmark method is employed to find the dynamic responses of BDFG subjected to a moving harmonic load. To validate the present formulation and solution method, some comparisons of the obtained fundamental frequency and dynamic response with those available in the literature are performed. A parametric study is performed to extensively explore the impact of the key parameters such as the gradient indices in both directions, moving speed, and excitation frequency of the acting load on the dynamic response of BDFG nanobeams. The obtained results can serve as a guideline for assessing the multi-functional and optimal design of micro/nanobeams acted upon by a moving load.

and numerical techniques to investigate the static and dynamic problems of FGM structures, whose material properties are graded in only one direction, i.e., transversely functionally graded in the thickness direction (TFG), or axially functionally graded in the length direction (AFG), (see Refs. [4][5][6]).
In the last decades, the design and analysis of nanomaterials and nanostructures have been extensively increased because of their superior mechanical, thermal, and electrical properties. Nanostructures such as nanorods, nanobeams, nanoplates, and nanoshells have recently been used in many applications in micro/nanoscale devices such as sensors and actuators. For modelling of micro/nanoscale structures, the size effect on their mechanical and physical properties should be considered. Since the classical continuum mechanics is size-independent, it cannot capture the effect of small size [7][8][9][10][11][12]. In this regard, the molecular dynamic approaches and size-dependent continuum mechanics are employed to include the small-scale effect in microand nanostructures [13,14]. Numerical simulations based on the molecular dynamic approaches are, however, generally complex and computationally expensive, especially for complicated structures with a huge number of atoms, and thus the nonclassical size-dependent continuum mechanics is extensively used. Several different nonclassical continuum mechanics theories including additional material length scale parameter(s) have been developed and employed to model the size-dependent phenomenon in miniaturized systems, such as the strain gradient theory "SGT" ( [15,16]), couple stress theory "CST" ( [17][18][19]), nonlocal continuum elasticity theory "NET" ( [20]), modified strain gradient theory "MSGT" ( [9]), nonlocal strain gradient theory ( [21]), modified couple stress theory "MCST" ( [22]), and surface elasticity theory "SET" ( [23,24]). Critical reviews on modelling of micro/nanostructures based on the different nonclassical continuum mechanics theories have been presented by [25,26]. The main difficulty encountered in the nonclassical continuum theories is the determination of the microstructural material-dependent length scale parameters. The MCST proposed by [22] has the advantage of involving only one additional microstructure-dependent length scale parameter for isotropic materials. Here, the MCST will be employed to capture the microstructure effect on the predicted response.
One of the most important characteristics of nanoscale structures is their extremely increasing ratio of surface area-to-bulk volume. As a result, the energy associated with the surface atoms is different from that of the bulk. Unlike in classical continuum mechanics, the surface energy cannot be neglected as it contributes significantly to the total energy. This is because the effects of the surface residual stress (surface tension) and surface elasticity are of the most important size-dependent effects in the analysis of the mechanical behavior of the nanoscale structures [46]. Gurtin and Murdoch ([23,24]) proposed a surface elasticity theory to approximate the contribution of surface energy by assuming an elastic body with elastic surface layers of zero thickness that is perfectly bonded to the surface of the bulk continuum. The Gurtin and Murdoch surface elasticity theory (GM-SET) has been successfully employed in conjunction with the other nonclassical continuum theories in the analysis of homogeneous and one-dimensional FG structures, i.e., GM-SET combined with the DNET ( [47][48][49]), GM-SET combined with the SGT ( [50]), GM-SET combined with the DNSGT ( [51]), GM-SET combined with the MCST ( [52][53][54][55][56][57][58][59][60][61][62][63][64]), and GM-SET combined with the MCST in the framework of DNET ( [65,66]). Based on these studies, it has been shown that incorporating the influence of surface energy may show a stiffness-hardening or stiffness-softening of the studied nanostructures depending on whether the signs of the surface elastic constants are positive or negative.
Considering the surface effect on BDFG nanobeams, only few models have been recently developed. Lal and Dangi [67] adopted the differential nonlocal elasticity theory together with GM-SET to investigate the vibration response of BDFG Timoshenko nanobeam using the differential quadrature method (DQM). Power law was used to express the gradation of the material properties of the bulk continuum in the thickness and length directions, whereas the surface elastic constants were assumed to vary in thickness direction only, and the nonlocal parameter was assumed to be spatially-independent. Shanab and Attia [68][69][70] studied the bending, buckling, and vibration characteristics of uniform and tapered BDFG Euler-Bernoulli and Timoshenko micro/nanobeams in the framework of the MCST and GM-SET. Unlike [67], all the material properties of the bulk and surface continuums were assumed to vary in both the thickness and length directions via power law.
Furthermore, analysis of structures' response under moving mass and moving load is very important for many practice engineering applications, such as bridges, tunnels, and rail. The vibration problems of a homogenous beam subjected to moving loads were extensively studied by [71][72][73][74][75]. The dynamic performance of one-directional FGM beams undergoing a moving load was studied based on the material gradation in the thickness direction by [76][77][78][79][80][81][82] or the material gradation in the length direction by [83][84][85]. Taking the features of BDFG into account,Şimşek [86] investigated the free and forced vibration responses of BDFG Timoshenko beams acted upon by moving loads using Ritz and the implicit time integration methods. Exponential functions were adopted to express the material variation in both thickness and length directions. Utilizing finite element method (FEM) and Newmark method, the vibration of a BDFG Timoshenko beam under a moving load was investigated by [87] assuming that the material properties were varied in both thickness and length directions via power law functions. Yang et al. [88] studied the vibration characteristics of a tapered BDFG Timoshenko beam subjected to a moving harmonic force employing the meshfree boundary-domain integral equation in conjunction with 2D elasticity theory. The exponential distribution was assumed for the variation of Young's modulus and mass density in transverse and axial directions. Employing Ritz method and Gram-Schmidt orthogonalization procedure, [89] studied the free and forced vibrations of functionally graded graphene nanoplatelet-reinforced beams due to multiple moving loads based on third-order shear deformation beam theory. Recently, [90] utilized FEM with the aid of Newmark method to explore the dynamic response of a sandwich Timoshenko beam with BDFG face layers under a moving load. Power law functions were adopted to express the material properties of the skin layers. This work was extended by [91] to explore the effect of a partial Pasternak support on the dynamic response of BDFG sandwich beams under a moving mass using a quasi-3D theory. Based on the third-order shear deformation theory, the influences of Coriolis and centrifugal forces on the dynamic response of an inclined BDFG sandwich beam were considered by [92]. For the double-BDFG porous Timoshenko beam system subjected to moving load, its vibration performance was studied by [93] adopting Ritz and Newmark methods. Based on Timoshenko beam theory, [94] utilized Ritz method with the aid of Newmark method to explore the free and dynamic responses of a sandwich beam with TFG face layers under harmonic moving loads with elastic foundation for different boundary conditions. Using Chebyshev collocation method, and based on the third-order shear deformation theory, the free vibration of TFG beam and sandwich plates resting on an elastic foundation was investigated by [95] and [96], respectively. Recently, [97] utilized the Gram-Schmidt procedure to generate the shape functions for the Ritz method, and then the Newmark method was used to explore the dynamic response of TFG sandwich plates under multiple moving loads.
To capture the scale phenomena of micro/nanoscale beams acted upon by moving loads, [98,99] adopted the DNET to study the dynamic response of simply supported TFG Euler-Bernoulli nanobeams subjected to a constant moving load. Based on DNET, the influences of surface energy and viscoelastic foundation on the steady-state response of Euler-Bernoulli nanobeam in a thermal environment and subjected to a moving concentrated load were examined by [100] using the multiple scales method. In [101], the forced vibration response of TFG nanobeams resting on Winkler-Pasternak foundation and under a uniform harmonic dynamic load was investigated employing a higher-order shear deformation beam theory in the context of DNET. In the context of MCST, [102] used the FEM and Wilson-theta method to study the forced vibration of nonuniform BDFG microbeams resting on a linear elastic foundation and acted upon by a moving harmonic load/mass. The influence of various models of the material gradation was presented. In the framework of parabolic shear deformation theory (PSDPT), [103] employed the MCST to study the vibration response of BDFG porous microbeams excited by a moving harmonic load using FEM. For even and uneven porosity distributions, power law functions were adopted to model the material variation in both thickness and length directions. The authors [104] adopted the MCST to explore the effects of thermal rise and moving load on the vibration response of a BDFG microbeam using FEM. Temperature-dependent material properties were assumed with a power law distribution in both the thickness and length directions. In the framework of NSGT, the dynamic performance of TFG nanobeams under a moving load was studied by [105,106].
Based on the above-mentioned review, it is observed that the dynamic performance of BDFG nanobeams is still limited. The main objective of the present study is to develop an integrated microstructure-surface energy-based model for the dynamic response of BDFG nanobeams under a moving harmonic load based on a higher-order shear deformation beam theory, for the first time. Effects of microstructure and surface energy are captured via the MCST and GM-SET, respectively. The material properties of the bulk and surface continuums of the nanobeam are varying gradually along both the thickness and length directions according to power law functions. Based on the generalized elasticity theory, the formulation accounts for the physical neutral axis concept. Employing Lagrange's equation, the nonclassical equations of motion and boundary conditions are derived. Ritz method in conjunction with Newmark method is used to obtain the dynamic deflection of a BDFG nanobeam under moving load. Both the present formulation and solution procedure are verified by comparing the predicted results with previous studies. By a parametric study, the effects of gradient indices Velocity Fig. 1 Illustration of a BDFG nanobeam with surface layers exposed to a moving load in both thickness and length directions, velocity of moving load, excitation frequency, microstructure, and surface energy on the dynamic response of simply supported BDFG nanobeams are explored in detail.

Theory and formulation
Consider a straight uniform BDFG microbeam with length L, width b, and thickness h, as demonstrated in Fig. 1, in a Cartesian coordinate system (x, y, z) that denotes the geometrical neutral plane (midplane). The beam is excited by a moving load with a constant velocity v. The BDFG beam is composed of a mixture of ceramic and metallic constituents, where the lowermost (x 0, z −h/2) and uppermost (x L, z h/2) surfaces of the nanobeam are pure metal "m" and pure ceramic "c", respectively. The effective material properties describing the bulk and surface continuums of BDFG can be expressed via power law in both axial and transverse directions, such that [40,41,102]: where the superscripts "s" and "B" refer, respectively, to the surface and bulk continuums of the beam. E B , ν, and ρ B are, respectively, Young's modulus, Poisson's ratio, and the mass density of the bulk continuum, l represents the material length-scale parameter describing the microstructure effect on the bulk continuum, and τ s , λ s , μ s , and ρ s are the surface residual stress, elastic parameters, and surface mass density of the surface continuum. k x and k z are the gradient indices in the axial and transverse directions, respectively. The classical Lamé's parameters of the bulk continuum are given by Ignoring the Poisson's effect yields λ B (x, z) + 2μ B (x, z) ≡ E B (x, z) as adopted by [43,44,103]. In Eqs. (1)(2)(3), the axially (AFG) and transversely (TFG) functionally graded distributions can be recovered by setting k z 0 and k x 0, respectively. An homogeneous beam made of pure metal constituent is obtained when k x k z 0.
Due to the nonsymmetric distribution of the material properties of the BDFG beam about its midplane, the geometrical neutral plane (GNP) does not coincide with the physical neutral plane (PNP), as demonstrated in Fig. 1. The deviation between the positions of GNP and PNP is given by [40,41,68,69] e n (x) Unlike the TFG beam, the position of the physical neutral axis "e n " in a BDFG beam depends on the x-coordinate. For simplicity, Cartesian coordinates are used to approximate the z n -coordinate instead of the curvilinear coordinates.

Kinematics and constitutive relations
In this study, a general shear deformation theory is used to express the kinematics of the BDFG nanobeam, in which the displacement field is given by where u and w are the displacement components of the midplane along, respectively, x and z directions, and φ is the transverse shear strain, and t denotes time. Accounting for the physical neutral plane concept, the shear-strain function R(z n ) R(z) − r n (x), where Various beam theories can be satisfied by appropriate section of the functions f (z n ) and R(z n ). Adopting the third-order parabolic shear deformable beam theory (PSDBT), [107], In the framework of generalized elasticity theory in combination with the modified couple stress theory, [22], the infinitesimal Green strain tensor ε, classical Cauchy stress tensor σ B , symmetric curvature tensor χ, and the deviatoric part of the couple stress tensor m are given by [55,65]: where u and θ represent the displacement field and rotation field vectors, respectively. l denotes the material length scale parameter (MLSP) of the material, which captures the size-dependent effect of the material microstructure. Here, the gradation of the MLSP in both thickness and length directions of the BDFG beam is considered, as depicted in Eq. (2). Based on the Gurtin-Murdoch surface elasticity theory, the surface layer of a bulk material is assumed to be of zero thickness and fulfills different constitutive equations involving the surface parameters, i.e., residual surface stress τ s and the surface elastic constants λ s and μ s , [23,24]. According to this theory, the in-plane (σ s ) and out-of-plane (σ s nα ) components of the surface stress tensor are given by Gurtin and Murdoch [23,24] as follows, respectively: In Eqs. (11), signs (+) and (−) stand for the upper and lower surface layers of the BDFG beam, respectively. n i denotes the components of the outward unit normal vector n to the beam lateral surface. The nonzero components of the Green strain and the symmetric curvature tensor in accordance with Eqs. (6), (8), (9.1), and (10.1) can be obtained as where subscript x represents the derivative w.r.to x. It should be mentioned that the derivative ∂ 2 R(z n ) ∂z∂ x will be zero.
In light of Eqs. (9.2) and (10.2), the nonzero components of the classical Cauchy stress tensors and the deviatoric part of the couple stress tensor are [40,41] Besides, the nonzero components of the surface stresses in the BDFG nanobeam based on PSDBT can be obtained by substituting Eqs. (6) and (12) into Eq. (11), such that where n y and n z denote the components of the unit outward normal vector n to the beam lateral surface in the y and z directions, respectively, i.e., with θ the angle between the y-axis and the normal vector n, then n y and n y are given by, respectively, cosθ and sinθ . The subscript "s" represents the direction of the unit tangent vector s on the boundary of the beam. σ s xz represents the in-plane shear stress component on the two lateral surfaces (parallel to the xz plane) of a beam with a uniform rectangular cross section. When θ 0, the values of σ s zx and σ s xz , equal those of σ s sx and σ s xs , respectively. It is important to emphasize the fact that the in-plane shear stress tensor defined by Eq. (11.1) in not symmetric, and thus σ s xz and σ s zx are not equal [53]. Many authors employed the Gurtin-Murdoch surface elasticity theory claiming a symmetric in-plane shear stress tensor in their variational models.

Variational formulation
In the context of the generalized elasticity theory in combination with the surface elasticity theory and modified couple stress theory, the total strain energy (U) of an isotropic elastic BDFG deformed nanobeam can be written as follows [55,65]: where A and ∂ A are, respectively, the cross-sectional area and the boundary of the beam. Substitution of Eqs. (14)(15)(16) into Eq. (17) yields Substituting Eqs. (14)(15)(16) and (19) into Eq. (18) and with some mathematical manipulations, the total strain energy of the BDFG nanobeam including the microstructure and surface energy effects can be obtained in terms of the displacement field as follows: It is worth noting that the stiffnesses B B x x (x) disappear in case of neutral axis analysis. The kinetic energy of the BDFG nanobeams accounting for the surface effects can be expressed as in which the mass moments of inertia, incorporating the mass density of the bulk and surface continuums, are A general form for the virtual work done by the forces applied on the current beam in the context of the modified couple stress theory and surface elasticity theory can be expressed as [53,65]: where f and f c are, respectively, the body force resultant and body couple resultant, both per unit volume, t and s are, respectively, the traction resultant and surface couple resultant, both per unit area. In the absence of an applied compressive force, assuming the externally applied harmonic load moves with a constant speed ignoring its inertial effect, and employing zero initial conditions, the virtual work can be obtained as where v is the velocity of the moving harmonic load, δ represents the Dirac-delta function, f u and q are, respectively, the distributed loads in x-and z-directions per unit length along the x-axis, f c is the y-component of the body couple per unit length along the x-axis. N , V , M c , and M nc are, respectively, the applied axial force and lateral force, classical external bending moment, and the nonclassical external bending moment due to the couple stress exerted at the two ends of the beam, and a c is defined as The applied external moving harmonic load is given by in which P 0 and are, respectively, the amplitude and frequency of the applied moving load. The location of the moving load at any instant, measured from the left end of the beam, is defined as Finally, the total energy of the BDFG nanobeam is given as

Solution procedure
In this Section, the Lagrange's equations are used to obtain the system of equations of motion. For this purpose, trigonometric Ritz method (TRM) is applied first. The displacement functions w(x, t) u(x, t), and φ(x, t) are approximated by series of trigonometric functions that satisfy the geometric simply supported boundary conditions as follows: where A r (t), B r (t), and C r (t) are the unknown time-dependent coefficients to be determined. The admissible trigonometric functions are defined as follows: Substituting Eqs. (31,32) into Eq. (30), and then using Lagrange's equations given by Eq.
in which Applying the Lagrange equations yields the following system of equations of motion: where the elements of the stiffness matrix [K] and the mass matrix [M] are, respectively, defined as Due to the presence of the three terms C s ∂φ ∂ x in the total strain energy, Eq. (20), an internal right-hand side is obtained when applying the Lagrange equations. These terms are nonzero when including the effect of the surface residual stress of the BDFG nanobeam, i.e., when τ s 0 and k z 0 as seen from Eq. (22.5). So, these terms act as self-excitation loading and cause deformation of the nanostructure at no external load [63].
The force vector F SE represents a self-excitation force resulting from the surface energy contribution, i.e., and the force vector's F t components due to the external moving harmonic load are given by The system of equations of motion, Eq. (35), is solved by using the implicit time integration Newmark method ( [108,109]). After obtaining the time-dependent coefficients {A r (t), B r (t), C r (t)}, the displacements, velocities, and accelerations of the beam at the considered point and time are determined for any time t between For better numerical analysis of the forced vibration, the following nondimensional quantities are introduced: where w(x, t) represents the normalized dynamic deflections of the BDFG beam, D 0 is the deflection of a pure metal beam under a point load P 0 at mid-span (D 0 P 0 L 3 /48E m I ), I bh 3 /12. ω 1 is the dimensionless fundamental frequency of free vibration analysis. v indicates the dimensionless velocity of the external load, where V c is the critical velocity (the velocity at which the maximum magnitudes of the maximum deflections occur). τ denotes the dimensionless time, and 0 < τ < 1. It is worth mentioning that, when τ is larger than one, the load leaves away from the BDFG microbeam. In addition, for the Newmark procedure, 2000-time increments are adopted for each dynamic response.   Table 1 compares the predicted maximum dynamic magnification factors, i.e., peak values of the maximum normalized deflections w p and the corresponding absolute velocities v p of the simply supported TFG SUS304/Al 2 O 3 beam with the results reported by Simsek and Kocatürk [76] and Nguyen et al. [87], for Euler-Bernoulli beam theory (EBT) and Timoshenko beam theory (TBT), respectively, at f 0. In the same example, the variation of the maximum normalized deflection at the center of a simply supported TFG beam with the power-law exponent k z for f 40 and 80 rad/s is also predicted and compared with by [76], as shown in Fig. 2. It is depicted that the present results well agree with the published one for the dynamic response of a TFG beam based on the classical formulation.
In the second example, the dynamic response of the simply supported BDFG microbeam predicted by the present model and that obtained by Zhang and Liu [103] is compared. The material properties of SUS304/Al 2 O 3 mentioned above are used with a constant material length scale parameter of l m l c 0.5h. The beam dimensions are b 1 μm, h 10 μm, and L 20h. Based on the MCST, the dimensionless central deflections are plotted in Fig. 3 at a dimensionless moving velocity of v 0.1. Under these conditions, the maximum dynamic magnification factors, i.e., peak values of the maximum normalized deflections w p , and the corresponding dimensionless velocities v p are obtained and compared as provided in Table 2. A very good agreement between the present results with those of [103] is observed from Fig. 3 and Table 2.
Finally, the third validation example compares the predicted dimensionless fundamental frequency (ω 1 ) of the present PSDBT simply supported BDFG with that reported by Shafiei et al. [110], as illustrated in Fig. 4. In this example, Young's modulus, Poisson ratio, and mass density are, respectively, E m 201.04 GPa, ν m   3800 kg/m 3 for the ceramic material. The geometrical parameters of the microbeam are h 15 μm, b h, and L 100h, and the material length scale parameter is taken as l m l c 0.25h. It is observed from Fig. 4 that there is a good agreement between the present model for the free vibration response of BDFG microbeams. Considering the comparison presented in Figs. 2, 3, and 4 and Tables 1 and 2, it can be concluded that the developed model and solution procedure can accurately predict the dynamic response of BDFG beams subjected to a moving load.

Results and discussion
In this Section, numerical results are presented to investigate the effect of the gradient indices in length and thickness directions, velocity of the moving load, and excitation frequency on the dynamic response of simply supported BDFG nanobeams exposed to a moving load. The influences of microstructure and surface energy are considered. Unlike most of the existing models, the present model accounts for the exact location of the physical neutral axis, and the gradation of Poisson's ratio, material length scale parameter, residual surface stress, two surface elastic constants, and surface mass density. The BDFG beam is made of a mixture of aluminum (metal) and silicon (Si) materials. Young's modulus, Poisson's ratio, and the bulk mass density of the metallic constituent are, respectively, .1688 kg/m 2 , respectively, [70]. Since the value of the material length scale parameter (MLSP) differs from a material to another, it is not reasonable to assume a constant MLSP in FGMs. Unfortunately, there is  Figure 5 shows the dependency of the distances indicating the location of the physical neutral axis (e n and r n ) on the thickness and length gradient indices of the nanobeam. The influence of the thickness and length gradient indices on the variation of the maximum normalized dynamic deflection (dynamic magnification factor) at the center of the simply supported BDFG Al/Si nanobeam versus the dimensionless velocity of the moving load is illustrated in Figs. 6 and 7, respectively. To clearly explore the nonclassical effects due to the small-scale size, the present results are predicted based on the classical "CL" formulation when l 0 and τ s μ s λ s ρ s 0, couple stress "CS" formulation in the absence of the surface effect (τ s μ s λ s ρ s 0), surface energy "SE" formulation when the couple stress (microstructure) effect is ignored (l 0), and full-nonclassical "CSSE" formulation which incorporates the simultaneous effects of microstructure and the surface energy.
From Figs. 6 and 7, it is depicted that the maximum normalized dynamic deflection in the entire time history both increases and decreases, it then increases to reach the peak value when the moving load velocity reaches a certain value, and then after the maximum normalized dynamic deflection gradually decreases as the moving load velocity increases [113]. The velocity at which the maximum normalized dynamic deflection attains its peak value is referred to as the critical velocity of the BDFG beam [114]. At lower moving load velocities, the repeated increase and decrease of the maximum normalized dynamic deflection is due to the beam oscillations. Regardless of the velocity and the gradient indices, incorporating the nonclassical effects significantly reduces the predicted maximum normalized dynamic deflection. The highest maximum normalized dynamic deflection is associated with the classical elasticity-based formulation, followed with SE, CS, and CSSE formulations, which is due to the stiffness-hardening effect because of the couple stress and surface energy. A stiffer nanobeam resists a moving load better than a soft one. Generally, the trends of the maximum normalized dynamic deflection-dimensionless moving velocity curves are independent of the gradient indices and the nonclassical formulations. Peak values of the maximum normalized dynamic deflections (w p ) and their corresponding absolute and dimensionless critical velocities (v p and v p , respectively) are provided in Table 3 at various values of the gradient indices of the BDFG nanobeam and different formulations. It can also be discerned from Table 3 and Figs. 6 and 7 that increasing the gradation indices in thickness or length direction, k z or k x , respectively, shows a considerable increase in the maximum normalized dynamic deflections, which is observed for CL, CS, SE, and CSSE formulations. This is due to the fact that increasing k z and/or k x increases the volume fraction of metal constituent, and thus the effective rigidity of the beam decreases. Additionally, the influence of the gradient index in the length direction (k x ) on the dynamic response is much larger than that in the thickness direction (k z ). It is also depicted that the influence of the gradient indices is more pronounced when their values are low. Keeping k x 0.5, increasing k z from 0.0 to 1, shows an increase in w p by about 30.66, 17.78, 23.31, and 15.72%, and a decrease in v p by 14.32, 10.06, 12.60, and 9.98% for CL, CS, SE, and CSSE formulations, respectively, while increasing k z from 2 to 10 under the same conditions results in increasing w p by about 16.13, 11.74, 11.37, and 10.15% and reducing v p by 7.47, 6.56, 5.71, and 5.62% for CL, CS, SE, and CSSE formulations, respectively. On the other hand, increasing k x from 0.0 to 1, while keeping k z 0.5, increases w p by 30.49, 16.311, 23.39, and 14.59% and reduces v p by 16.32, 10.17, 9.43, and 10.88%, while increasing k x from 2 to 10 under the same conditions increases w p by 16.36, 12.63, 10.97, and 10.73% and reduces v p by 7.12, 5.39, 5.08, and 5.52% for CL, CS, SE, and CSSE formulations, respectively.
Furthermore, the results presented in Table 3 and Figs. 6 and 7 reveal that there is a considerable influence of the nonclassical microstructure and surface energy on the dynamic response of BDFG nanobeams. For the BDFG nanobeam under investigation, the highest and lowest nonclassical effects are associated with the CSSE and SE formulations, respectively. The maximum normalized dynamic deflections predicted using CS, SE, and CSSE formulations are, respectively, 0.395, 0.751, and 0.35 times that predicted using CL formulation when k x 0.5 and k z 1, and are, respectively, 0.382, 0.742, 0.338 when k x 1 and k z 0.5. Employing classical Table 3 Comparisons of the peaks of maximum normalized deflection at mid-span w p of the simply supported BDFG nanobeam and corresponding velocities at different gradation indices using various formulations k x k z Classical elasticity-based formulation (CL) Couple stress-based formulation (CS) Surface energy-based formulation (SE) Couple stress-surface energy-based formulation (CSSE)  elasticity-based formulation yields a considerable underestimation of the critical velocity v p . However, the impact of surface energy and couple stress on the dynamic response increases by increasing the bi-directional gradient indices. Also, the effect of the thickness gradient index on the contribution of surface energy and/or couple stress is slightly larger than that of the length gradient index.
For examining the mutual effect of the thickness and length gradient indices on the dynamic response of a BDFG nanobeam under moving load in more detail, the maximum normalized dynamic deflection w max as a function of the gradation indices k x and k z is shown in Figs. 8, 9, and10 for three different dimensionless moving load velocities, v 0.1, 0.6, and 0.8, respectively. It is seen that as the gradient indices k x and/or k z increase, w max increases, regardless of the velocity of moving load and the formulation type. When k x k z 0, the beam is made of purely ceramic "silicon", and thus the beam attains its lowest dynamic deflection, while the highest dynamic deflection is reached when k x k z 10 as the beam is almost pure metal. Based on the results in Figs. 6, 7, 8, 9, and 10, it can be extracted that the dynamic response of a BDFG nanobeam can be controlled by the selection of the gradient indices.  Table 4 records the numerical values of the dimensionless central deflection of a simply supported BDFG nanobeam for different values of the gradation indices and dimensionless velocities of the acting moving load. It is found that the gradient indices, moving load velocity, and the employed formulation have a significant role in the dynamic behavior of BDFG nanobeams. The gradation indices have a considerable influence on the amplitude of dynamic deflection, but they slightly affect the curve shapes of the time histories, regardless of the adopted formulation. It is indicated that the moving load velocity has a significant impact on the dynamic deflection response, and low velocities of the travelling load increase the number of vibration cycles of a BDFG nanobeam, which is attributed to the low ratio of the moving load velocity to the critical velocity. Figure 13 demonstrates the variations of the normalized dynamic deflection at the center of a simply supported BDFG nanobeam versus the dimensionless velocity of acting moving load and dimensionless time for k x k z 1. Various formulations are presented to show the influence of microstructure and surface energy on the predicted dynamic deflection. It is observed that the inclusion of couple stress and/or surface energy considerably decreases the predicted dynamic deflection during the time interval regardless of the moving load velocity. Employing CSSE formulation gives the smallest normalized dynamic deflection, whereas the classical formulation gives the largest one. For the different formulations, it is depicted that a faster moving load, i.e., higher moving velocity, takes less traveling time, and thus the BDFG nanobeam exhibits low fluctuations, and these fluctuations are almost not obvious the dimensionless velocity approaches 0.6.
The relation between the maximum normalized dynamic deflection at the mid-span of a BDFG nanobeam and the dimensionless velocity of moving load is shown in Fig. 14 at different values of the excitation frequency ratio, i.e., r f /ω 1 , of 0.1, 0.2, and 0.4. The gradation indices are fixed at k x =k z =1. For the classical and nonclassical formulations, it is obvious that while the dimensionless velocity increases, the trends of the dynamic deflection versus the dimensionless time are not influenced for different values of the excitation frequency ratio. Also, as the excitation frequency ratio increases, the peak value of the maximum normalized dynamic deflection and its corresponding dimensionless velocity are remarkably increased. Figure 15 presents the variation of the maximum normalized dynamic deflection with the excitation frequency for different gradient indices of a BDFG nanobeam employing various formulations. The dimensionless moving velocity is kept constant at v 0.2. It is noticed that the maximum dynamic deflection reaches its peak value at some value of the excitation frequency, which is the fundamental frequency of the BDFG nanobeam. It is noted that the fundamental frequency of the BDFG nanobeam reduces by increasing the gradation indices k x and k z . On the other hand, increasing the gradation indices yields a remarkable increase in the predicted peak values of the deflection curves, for both classical and nonclassical formulations. Inclusion of the microstructure and/or surface energy effects significantly reduces the maximum dynamic deflection and increases the fundamental frequency.
The required time for the moving load to reach the right end of the nanobeam, i.e., x p L, is t * π/(vω 1 ), and the corresponding load becomes P(t * ) P 0 sin(πr /v), where r /v represents the dimensionless time taken by the moving load to move across the entire nanobeam, Zhang and Liu [103]. In this regard, consider a BDFG nanobeam with k x 1 and k z 1 and subjected to a moving load with a dimensionless moving velocity of v 0.2. The time-history of the normalized dynamic deflection at the mid-span of the nanobeam is provided in Figs. 16 and 17 for, respectively, r ≤ 0.2 and r ≥ 0.2. When r ≤ 0.2, it is depicted from Fig. 16 that lower values of the excitation frequency ratios lead to lower values of the amplitudes of the dynamic deflections at the mid-span of the BDFG nanobeam. On the other hand, when r > 0.2, the resonance phenomenon is  Table 4 Comparisons of the maximum dimensionless deflection at mid-span Couple stress-surface energy-based formulation (CSSE)  observed when the fundamental frequency of the considered BDFG nanobeam equals the excitation frequency of the applied moving load, i.e., r f /ω 1 1. Therefore, the normalized mid-span dynamic deflection reaches it largest value in the resonance case. As the excitation frequency ratio increases more than 1, i.e., the fundamental frequency of the beam is lower than the excitation frequency of moving load, the periodic load decreases, and consequently, the amplitude of the dynamic deflection significantly decreases. When the excitation frequency ratio approaches 3, the periodic load is about zero, and thus, the dynamic deflection is almost vanishing.

Conclusions
In the present work, the dynamic response of bi-directional FG nanobeams acted upon by a moving concentrated harmonic load is investigated. A new nonclassical beam model is developed in the framework of higher-order shear deformation beam theory in conjunction with the modified couple stress theory and the Gurtin-Murdoch surface elasticity theory. The proposed model accounts for the bi-directional gradation of the material length scale parameter to capture the microstructure effect as well as the bi-directional gradation of the surface residual stress, two surface lamé constants, and the surface mass density to include the surface energy effect. Introducing the physical neutral plane concept, the equations of motion are derived using Lagrange's equations.
In the solution procedure, the trigonometric Ritz method is employed, and the dynamic response in the time domain is obtained by Newmark method. An extensive numerical parametric study is performed to examine the influence of some key parameters such as the gradation indices in thickness and length directions, moving load velocity, and the excitation frequency on the dynamic response of a simply supported BDFG nanobeam using the four different formulations. Some important conclusions can be summarized as follows: • The dynamic deflection employing the nonclassical formulations is considerably smaller than that predicted by the classical formulation, regardless of the gradient indices, due to the stiffening induced by surface energy. • The normalized dynamic deflection and the maximum deflection are significantly increased as the gradation indices increase, whereas the corresponding velocity is considerably decreased. • Incorporating the nonclassical surface energy and/or microstructure effects reduces the influence of the gradient indices. • The influence of the gradient index in thickness direction on the dynamic response is slightly lower than that in the length direction. • As the velocity of moving load increases, the predicted normalized maximum dynamic deflection increases first and then decreases as the velocity becomes larger than its critical value. • Increasing the velocity of the applied moving load reduces the number of the dynamic vibration cycles, whereas the number of the vibration cycles is almost unaffected by the gradient indices and the nonclassical contribution. • The predicted dimensionless mid-span deflection of a BDFG nanobeam value is increased as the excitation frequency increases and reaches its peak value as the fundamental frequency of a BDFG nanobeam equals the excitation frequency. Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.