Free vibration characteristics of sectioned unidirectional/bidirectional functionally graded material cantilever beams based on finite element analysis

Advancements in manufacturing technology, including the rapid development of additive manufacturing (AM), allow the fabrication of complex functionally graded material (FGM) sectioned beams. Portions of these beams may be made from different materials with possibly different gradients of material properties. The present work proposes models to investigate the free vibration of FGM sectioned beams based on one-dimensional (1D) finite element analysis. For this purpose, a sample beam is divided into discrete elements, and the total energy stored in each element during vibration is computed by considering either Timoshenko or Euler-Bernoulli beam theories. Then, Hamilton’s principle is used to derive the equations of motion for the beam. The effects of material properties and dimensions of FGM sections on the beam’s natural frequencies and their corresponding mode shapes are then investigated based on a dynamic Timoshenko model (TM). The presented model is validated by comparison with three-dimensional (3D) finite element simulations of the first three mode shapes of the beam.


Introduction
Functionally graded materials (FGMs) are advanced material systems characterized by progressively varying properties, which allow them to be used in a wide range of engineering * Citation: VIET, N. V., ZAKI, W., and WANG, Q. Free vibration characteristics of sectioned unidirectional/bidirectional functionally graded material cantilever beams based on finite element analysis. applications [1] . In particular, FGMs can be designed to meet specific requirements in terms of corrosion resistance, weight reduction, toughness, strength, wear resistance, heat isolation, etc. These requirements may be needed in extreme environments such as nuclear power generation [2] , aerospace applications [3] , engines [4][5] , and biomedical engineering [6] .
FGM-like properties are found in the fabrics of living beings such as bamboo, human teeth, and bone, but do not typically occur in conventional engineering materials. The fabrication of FGMs has recently become more accessible with the development of the additive manufacturing (AM) technique [7] . Indeed, AM allows the fabrication of sectioned FGM structures of relatively large sizes [8] with the ability to change the grading direction and the material type used within the geometry. Such structures are of critical importance in applications where different parts of the structure may be subject to significantly variable environmental conditions [8] . Moreover, FGMs can be utilized to allow gradual transitions among joined materials of significantly different properties. Such transition is desirable to prevent strong gradients of mechanical properties, for example, which can result in dangerous stress-concentrations.
FGM-sectioned beams have been increasingly investigated recently, mainly with experimental methods. Indeed, Bobbio et al. [8] presented an approach to fabricate a beam graded from Ti-6Al-4V ("V" stands for the chemical element "vanadium") to invar 36 (w(Fe) = 64% and w(Ni) = 36%) based on AM. The phase, composition, microstructure, and micro hardness were indicated as functions of positions in the FGM. The component was divided into two monolithic elastic regions separated by an FGM. Lima et al. [6] proposed a method to successfully fabricate graded Ti-Nb-Zr alloys using laser engineered net shaping. The sample was divided into five graded material sections with a view of applying the technique to current and future biomedical applications, such as the fabrication of orthopedic implants. Carroll et al. [9] used directed energy deposition to produce an FGM consisting of 304L stainless steel incrementally graded to inconel 625. The position-dependent properties of the FGM, including micro hardness, microstructure, and chemical and phase composition were determined by energy dispersive spectroscopy, microscopy, and X-ray diffraction. The work utilizes thermodynamic computation and experimental procedures for the fabrication and evaluation of the FGM. Fabrications of FGM and FGM-sectioned beams have been carried out by several research groups in recent years. Liu et al. [10] fabricated a novel Ti/Al lightweight graded beam using direct laser deposition, with the composition continuously changing from ϕ(Ti6Al4V) = 100% to ϕ(AlSi10Mg) = 100% along the length. Bobbio et al. [11] used AM to fabricate a sectioned beam consisting of FGM stainless steel 304L to V to Ti-6Al-4V. Bobbio et al. [12] later presented thermodynamic computations and experimental observations of an FGM of V to invar 36 beam based on AM. Zhang et al. [13] fabricated an SS316LIN625 FGM sectioned beam by powder-fed directed energy deposition. Hu et al. [14] fabricated FGM beams made from aluminum alloy using directional solidification. Rezapoor et al. [15] fabricated an Fe-TiC FGM wear resistant coating on the CK45 steel foundation by plasma spray and characterized the properties of the coating.
The development of dynamic models for FGM beams has been carried out by numerous scholars from theoretical and finite element modeling backgrounds [16][17][18][19][20][21] . In this regard, Wang and Li [16] investigated the free vibration characteristics of FGM beams using the Levinson beam theory. The governing equations of vibration were derived by directly integrating the stress form of the governing equations, including higher-order shear strain. The effects of the length-to-depth ratio, gradient index, and boundary conditions on the vibrational behavior of the beam were reported. Li et al. [17] studied the free vibration of FGM beams using classical and first-order shear deformation beam models. The derivation of the motion equations for the beams considered the shear deformation and the axial, transversal, rotational, and axialrotational inertial forces. This work presented a simple and useful model to compute the natural frequencies of FGM beams without dealing with the tension-bending coupling problem. Li et al. [18] proposed a model to study the free vibration of post-buckled FGM beams coupled with piezoelectric layers under both voltage and heat load, based on the Euler-Bernoulli beam theory.
The authors reported that the three lowest frequencies of the pre-buckled beam decreased with the increase in temperature, whereas those of a buckled beam increased monotonically. In addition, the tension force induced by the piezoelectric layers was found to efficiently increase the critical buckling temperature and the natural frequency. Hadji et al. [19] employed a fourvariable refined plate model for studying the free vibration of rectangular FGM sandwich plates. The governing equation was derived in this case based on Hamilton's principle, and the Navier theory was used to solve the system of governing equations. Cao and Gao [20] introduced a model based on the asymptotic development method to study the free vibration of non-uniform axial FGM beams. The developed model is simple and efficient to analytically study the free vibration of non-uniform axial FGM beams. Moreover, the model is applicable to tapered beams with randomly varying cross width. Lal et al. [21] studied the second-order statistics by standard deviation and mean of normalized nonlinear transverse dynamic central deflection responses of an un-damped elastically supported FGM beam with surface-bonded piezoelectric layers under the action of a moving load. The effects of temperature increments, volume fraction exponents, moving loads and velocity, nonlinearity, slenderness ratio, and foundation parameters on the FGM beam were reported. Simsek [22] introduced a dynamic model for bidirectional FGM beams using Euler and Timoshenko theories, in which the derivation of motion equations was carried out by using Lagrange's method. Simsek and Kocaturk [23] further presented an analytical model to study the forced and free vibrations of an FGM beam in presence of a concentrated moving harmonic load. Kapuria et al. [24] studied the dynamic properties of a layered FGM beam using theoretical and experimental approaches.
Despite the importance of the references cited here, none of them considered the development of dynamic models for additive manufactured FGM sectioned beams, wherein multiple sections with different materials or grading of properties may exist along the beam length. The present work focuses on the development and validation of such models. In addition to application to beams with explicit sectioning, the model is applicable to beams in which the properties of the material can evolve under the influence of an external factor. This includes, for example, superelastic shape memory beams in which the local stiffness depends on the applied load, which may result in phase transformation leading to the formation of multiple solid phase regions with varying material properties [25][26][27] . The model proposed in this paper is based on one-dimensional (1D) finite element modeling of sectioned FGM beams including the Euler-Bernoulli model (EBM) and the Timoshenko model (TM). Each FGM section is divided into elements of finite length, and equations are derived for the kinetic and potential energies stored in each element. The derived equations are then plugged into Hamilton's equation to obtain the equations of motion of the beam. The accuracy of the two models (EBM and TM) is evaluated by comparison with a three-dimensional (3D) finite element eigenvalue analysis of the first three mode shapes of the beam. The solution for the eigenvalues considering damping in the beam is also presented. Moreover, the influence of material properties and length of the FGM sections on the natural frequencies of the beam is investigated for the first four modes.
In this manuscript, Section 2 is dedicated to the development of the governing equations of motion using Timoshenko and Euler-Bernoulli beam theories, followed by 3D finite element model (FEM) of a sectioned FGM beam in Section 3. Section 4 then presents results obtained using the proposed models with comparison with 3D FEM simulations, followed by conclusions in Section 5.

Development of governing equations of motion
The FGM beam is considered to consist of n individual FGM sections in the direction of the length. The beam is represented in Fig. 1 with respect to the Cartesian system of coordinates Oxyz.
The beam of length L has a uniform rectangular cross section of height h and width b. The sections denoted by FGM-m, m ∈ (1, 2, · · · , n), are shown in Fig. 1. Each section consists of a different bidirectional FGM.

EBM
Based on the EBM, the transverse and axial displacements in the beam, denoted by w and u, respectively, at any location x and time t, can be computed as where w 0 (x, t) and u 0 (x, t) represent the transverse and axial displacements, respectively, of points located on the middle plane of the beam. The axial strain, ε xx , can be obtained from the displacement field using the following relation: Equations (1) and (2) can be transformed into the following matrix forms: where the subscripts ", x" and ", xx" indicate the first and second partial derivatives with respect to the coordinate x, respectively. Without loss of generality, the analysis is carried out by considering one arbitrary beam section FGM-m (see Fig. 2) of length l m . In Fig. 2, x m is the local longitudinal coordinate along FGM-m, such that x m = 0 at the interface between beam sections FGM-(m − 1) and FGM-m, and Oxyz is the local coordinate system for the element m. The section FGM-m is discretized into finite elements of uniform length l e . One such element, between nodes im and jm, is shown in Fig. 2(b). Each node i is described by using three degrees of freedom (DOFs), i.e., the axial displacement u i , the vertical displacement w i , and the slope θ i .
During vibration, the energy content of a single beam element is the sum of its potential and kinetic energies. The potential energy stored in one element of the beam is described in terms of the normal strain and stress as With the strain-stress relationship and the homogeneity of the beam in the direction z, the potential energy is further expressed as where E em represents Young's modulus of elements in section m.
In this work, the material properties within an FGM section are considered to vary in both vertical and axial directions according to the following exponential rule [23] : where R em is a material property representing either Young's modulus or the density of the considered element im-jm, R 0m is the reference value of the material property at the lowest corner on the left side (x m = 0, y m = −h/2) of the FGM section FGM-m, g my and g mx are gradient indices expressing variations along the vertical and longitudinal directions, respectively, and x im is the coordinate of the node im corresponding to the local coordinate x m (see Fig. 2(a)) of FGM-m. It can be noted from Eq. (7) that the FGM section becomes homogenous when the gradient indices are zero, i.e., for g my = g mx = 0. For the purpose of further developing the 1D FEM of the beam, the axial and transverse displacement fields between two element nodes are expressed as linear and cubic functions, respectively, in terms of the local coordinate x. The displacement field can, therefore, be written as where the following shape functions are introduced: Taking the first and second derivatives of u 0 and w 0 in Eq. (8) with respect to x, one has where are the first and second derivatives of the shape functions with respect to x, and δ em (t) = u i w i θ i u j w j θ T j represents the nodal degrees of freedom, which are functions of time t. From Eq. (9), the normal strain in a particular element can be expressed as By plugging Eq. (10) into Eq. (6), the potential energy of the element FGM-m is expressed as In addition, the stiffness of an element belonging to the section FGM-m is now expressed as The kinetic energy stored in an element of FGM-m is given by From Eq. (3), the transverse and axial displacement fields and slope of any point on the midplane of the element are further described in terms of shape functions as By correlating Eqs. (7), (13), and (14), the kinematic energy of the element is written as where the element's mass matrix M em is further expressed as The equations of motion of the beam, based on the Euler-Bernoulli theory, can now be derived by using Hamilton's method. In the interest of brevity, the details of further derivations are not included here. Instead, the reader is refered to the next section, in which more general derivations by utilizing a TM are presented. These derivations can be easily specialized to the case of an Euler-Bernoulli beam.

TM
Based on the TM, the axial strain ε xx and shear strain γ xy can be obtained as [28] ⎧ In this work, the shear correction factor for the rectangular cross section using the Timoshenko beam theory is selected as 5/6. Using the Timoshenko theory, each element has ten DOFs as shown in Fig. 3. The displacement and slope of each element are described in terms of ten DOFs by means of the following equations: where δ em = (u i , w i , θ i , w i,x , θ i,x , u j , w j , θ j , w j,x , θ j,x ) T , and the shape functions for one element are The equations of the shape functions are obtained by using the following assumed forms of the slope and displacements [29] : where a q , b q , and c q are constant coefficients of the cubic polynomial function expressing the mode shape of deflection and slope in terms of x [30] . The displacement field of the Timoshenko beam can be written as The displacements and slopes of the Timoshenko beam can be described in terms of the DOFs as ⎛ The potential energy for one finite element FGM-m is written as where G em is the shear modulus, given by G em = E em /(2 + 2v). The notation v stands for Poisson's ratio, chosen here to be equal to 0.3. Using Eqs. (17), (18), and (19) for normal and shear strains and then substituting them into the expression of the potential energy in Eq. (24), the following form is obtained for the potential energy: where the mass matrix is The kinetic energy of one Timoshenko finite element is a combination of transverse, axial, and rotational motions, which can be expressed as By plugging the expressions of the displacements and slopes in Eqs. (21), (22), and (23) into Eq. (27), the kinetic energy is further written as In addition, the stiffness of an element belonging to the section FGM-m is expressed as Hamilton's principle applied to δ em (t) between two specified times t 1 and t 2 in an element is given by Inserting Eqs. (25) and (28) into Eq. (30) and then assembling all equations of motion of the considered elements throughout the beam, the following equation is obtained: where M and K are the global mass and stiffness matrices found respectively by assembling the local mass and stiffness matrices of all the beam elements. The notation δ Σ is a combination of all nodal DOFs of the entire beam. In a general case that considers damping, Eq. (31) is rewritten as follows: where C is the global damping matrix. Following Ref. [30], Eq. (32) is transformed to the general matrix form, where δ 0Σ is the amplitude of the displacement δ Σ . By solving Eq. (33) numerically for eigenvalues and eigenvectors, the natural frequencies and mode shapes of the FGM-sectioned beam are, respectively, determined.

3D FEM
The finite element software ABAQUS is used for the simulation of the 3D FGM cantilever beam structure. The dimensions and material properties are analogous to those used in the finite element simulations. To create the FGM layers in ABAQUS, the bottom-up mesh function built into the software is used to mesh the FGM layer into 20 layers of uniform thickness, and subsequently, uniform stiffness and density are assigned to each layer using the method outlined in Ref. [27], governing the evolution of the material properties. 3D linear brick elements with reduced integration (ABAQUS designation C3D8R) are utilized with a constant Poisson's ratio v = 0.3.

Results and discussion
The proposed 1D model is used to analyze a beam with three FGM sections, and the accuracy is validated by comparison with 3D FEM simulation results using ABAQUS. Subsequently, the effect of the lengths and material properties of different FGM sections on the natural frequencies is studied. The properties of the reference material are those of steel so that E 01 = E 02 = E 03 = 200 GPa and ρ 01 = ρ 02 = ρ 03 = 8 050 kg/m 3 . The considered beam length is L = 0.6 m, with individual FGM section lengths l 1 = l 2 = l 3 = 0.2 m. The gradient indices for the material of FGM-1 are assumed to be zero (g 1x = g 1y = 0), unless otherwise noted. Moreover, damping is assumed to be negligible in the present work.

Validation by 3D FEM simulation
For the purpose of validating the model, the gradient indices of the section 2 are assumed to be g 2x = 5, g 2y = 0 and those of the section 3 are such that g 3x = 0, g 3y = 5. Firstly, the first three mode shapes of the transverse displacements obtained from the EBM, TM, and 3D FEM are shown in Fig. 4. It is seen from the figure that the EBM, TM, and 3D FEM are in good agreement for the first two mode shapes that correspond to the lowest two frequencies. In addition, among the simulated modes, the third mode shape shows the biggest deviation between 3D FEM and the other two models. The difference can be explained by the influence of Poisson's ratio in 3D FEM simulation, which is not considered in 1D analysis. In addition, some minor error may result from the way that gradient material properties are specified in 3D FEM simulations using discrete layers. Table 1 shows the natural frequencies of the first four modes for two different sets of gradient indices g 2x = g 3y = 5 and g 2x = g 3y = 10 obtained from EBM, TM, and 3D FEM simulation. Generally, the agreement is good for the first two modes with a better fit obtained when lower gradient indices are considered. The deviation reaches 27% for the third mode and 25% for the fourth mode for the case g 2x = g 3y = 10. The error is calculated as the absolute percent difference between the actual value obtained using the TM and the one obtained using the finite element analysis, which is used here as reference. The reason of this difference is similar to that discussed in Fig. 4. Moreover, it is observed that the use of the Timoshenko theory gives better results, which may be attributed to the ability of this theory to better approximate actual mode shapes because of overall higher compliance compared with the Euler-Bernoulli method. To further validate the model, the non-dimensional natural frequencies of the first three modes obtained using the present model based on the EBM are compared with the results in Ref. [31]. The comparison, reported in Fig. 5, shows a very good fit. The material properties and gradient power equation are used in a manner similar to that in Ref. [31]. In particular, the beam is made of an axial FGM with one end fixed and the other free. The FGM itself consists of a graded mixture of steel and alumina. The simulation results used for this validation are taken from Tables 10, 11, and 12 in Ref. [31]. It can be seen from the figure that both prediction methods exhibit good agreement, with the maximal percent difference close to 9.5% for the natural frequency of the third mode.

Dynamic behavior
Before investigating the vibration behavior of the beam, a convergence analysis of the frequency in terms of the number of elements is carried out for the case of the TM and reported in Fig. 6. In the figure, the gradient indices are selected as g 2x = g 3y = 5. Convergence to within 1% relative error is achieved using 30 or more elements along the length. In this work, 36 elements in total are systematically utilized for accuracy. The dynamical behavior of the beam is reported in Fig. 6 for the first four modes and for the following two sets of gradient indices: g 2x = g 3y = 5 and g 2x = g 3y = 10. Since the section FGM-1 has a uniform and relatively low elastic stiffness compared with other sections, a significantly stiffer response is observed in all cases as one moves to the right of x = 0.2 m, which defines the boundary between the first two FGM beam sections. Compared with the response of a homogeneous and elastic beam with uniform stiffness and density equal to those in FGM-1, the slope varies less significantly through sections FGM-2 and FGM-3. Moreover, the variation of the slope across the boundary between FGM-2 and FGM-3 is clearly non-smooth because of the abrupt change in material properties at this location. Figure 7 shows the normalized deflection and slope for the first four mode shapes of the beam with sections made of different reference materials along the beam length. Two types of beams are considered, respectively, denoted by Type 1 and Type 2 in Fig. 7. In Type 1, the section FGM-2 is made of aluminum with E 02 = 69 GPa and the density ρ 02 = 2 700 kg/m 3 , and the section FGM-3 is made of steel with E 03 = 200 GPa and the density ρ 03 = 8 050 kg/m 3 as the reference material. In Type 2, the reference materials used are steel for the section , the , the , the , the , the , the , the , the , the , the , the , the , the , the , the , the Fig. 7 Mechanical behavior of the beam for the first four modes with sections made of different materials, (a) normalized transverse deflection for the first and second modes, (b) normalized transverse deflection for the third and fourth modes, (c) normalized slope for the first and second modes, and (d) normalized slope for the third and fourth modes FGM-2 and aluminum for the FGM-3. In both cases, gradient indices are chosen as g 2x = g 2y = g 3x = g 3y = 5, g 1x = g 1y = 0, and the reference material of the section FGM-1 is made of a homogeneous steel element. The figure shows different profiles of displacements and slopes for the two types of beams for the modes higher than one. In addition, Type 2 shows larger displacements and slopes, compared with Type 1, in the sections FGM-2 and FGM-3. It can be concluded from the figure that redistributing the material along the beam leads to change in mechanical behavior. For the materials considered in this section, a smaller deflection of the beam tip is observed when heavier and stiffer materials are located in the vicinity of the fixed end of the beam. The shape of the slope and deflection curves in Figs. 7 and 8 is abrupt at the interfaces between two neighboring beam sections. This is due to the non-smooth gradation of material properties across these interfaces, despite the smooth gradation within each section.

Effects of gradient indices and section length
The effects of the gradient indices g 3y and g 3x of FGM-3 on the natural frequencies of the first four modes based on the TM are shown in Fig. 8, in which the gradient indices of FGM-2 are chosen as g 2x = g 2y = 5. It is seen from the figure that the frequency increases with a drop of both g 3x and g 3y in the first and second modes. For the third and fourth modes, the frequency increases with increasing g 3y and decreasing g 3x values. It can be seen that the axial gradient g 3x affects the frequency of the first four modes consistently. The phenomenon is due to the fact that increasing the gradient index g 3x leads to a gradually increasing density as one approaches the section tip, which results in a drop in the natural frequencies. Figure 9 demonstrates the effect of the gradient ratios g 2y /g 2x and g 3x /g 3y on the natural frequencies of the first four modes based on the TM. In this case, the gradient indices g 2x and g 3y are both chosen equal to 5.
It is obvious from the figure that deceasing both gradient index ratios significantly affects the frequencies of these first four modes. Figure 10 shows the influence of equal gradient indices g 2x = g 2y = g 3x = g 3y and length l 1 of the beam section FGM-1 on the natural frequency of the first four modes based on the TM. The mode frequency is found to decrease with increasing gradient index g 2x = g 2y = g 3x = g 3y where it is stable with the variation of length l 1 for the first and second modes. However, the frequency irregularly increases with the rise of l 1 for the third and fourth modes.

Conclusions
A dynamic model is developed to study the free vibration of a bidirectional FGM-sectioned beam by means of finite element analysis. The total energy of an element is computed by using both Euler-Bernoulli and Timoshenko beam theories. The obtained expressions are then plugged into Hamilton's principle to derive the equations of motion of the considered beam. Key observations regarding the proposed model are as follows.
(i) The EBM and TM are validated against 3D FEM simulations using ABAQUS. The results of the lowest two natural frequencies and their corresponding mode shapes obtained using all three prediction techniques are found to be in good agreement. However, the significant deviation is observed for the third mode shape and its corresponding natural frequency.
(ii) The gradient indices of the FGM sections significantly affect the modal behavior of the beam. It is found that increasing the gradient index, especially the one controlling the axial variation in material properties, leads to a decrease in the natural frequencies of the FGM sections and their corresponding mode shapes.
(iii) Exchanging different inclusion phases among FGM sections results in substantial change in the mode shapes of the beam. In particular, a smaller tip deflection is observed in the second, third, and fourth modes when the section made of the heavier material is the one toward the tip.
(iv) The variation of sectional length shows only a small effect on the natural frequencies and mode shapes of the beam.

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. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.