Natural Frequency and Mode Shapes of Exponential Tapered AFG Beams on Elastic Foundation

A displacement based semi-analytical method is utilized to study non-linear free vibration and mode shapes of an exponential tapered axially functionally graded (AFG) beam resting on an elastic foundation. In the present study geometric nonlinearity induced through large displacement is taken care of by non-linear strain-displacement relations. The beam is considered to be slender to neglect the rotary inertia and shear deformation effects. In the present paper at first the static problem is solved through an iterative scheme using a relaxation parameter and later on the subsequent dynamic analysis is carried out as a standard eigen value problem. Energy principles are used for the formulation of both the problems. The static problem is solved by using minimum potential energy principle whereas in case of dynamic problem Hamilton’s principle is employed. The free vibrational frequencies are tabulated for exponential taper profile subject to various boundary conditions and foundation stiffness. The dynamic behaviour of the system is presented in the form of backbone curves in dimensionless frequency-amplitude plane and in some particular case the mode shape results are furnished.


Introduction
The ceramic-metal gradation technique in structural materials initiated by the Japanese material scientists in Sendai indicated the beginning of an era of advanced class of materials known as Functionally Graded Materials (FGMs). Since then various possibilites of using FGMs for different structural applications have been explored [1]. These FGMs, first invented in 1984, are increasingly being utilized in aerospace, aircraft, automobile and defense industrial applications. Their popularity is mainly due to their ability of mitigating the problem caused by the sudden change of thermomechanical properties as in the case of laminated composites [2,3]. Since the behavior of structural members made up with FGMs are of critical importance, they have received great attention from researchers in the last decades.
For functionally graded beams, gradient variation may be orientated in the cross-section or/and in the axial direction. For the former, there have been a large number of researches devoted to understand the static, buckling and dynamic behavior of the functionally graded structural components. Various methods are employed, including analytical method, semi-analytical method and numerical methods such as finite element method (FEM), differential quadrature method (DQM) and harmonic differential quadrature method (HDQM), various meshless methods, and quadrature element method (QEM) etc. Su et al. [4] formulated the dynamic stiffness method to investigate the free vibration behaviour of functionally graded beams. In their analysis derivation of the governing differential equations of motion was done utilizing Hamilton's principle. Kodali et al. [5] used displacement field based on higher order shear deformation theory to represent the static behavior of functionally graded metal-ceramic (FGM) beams under ambient temperature. Jin and Wang [6] established an N-node novel weak form quadrature functionally graded (FG) beam element based on the classical beam theory (CBT) and differential quadrature (DQ) rule. Yaghoobi and Fereidoon [7] carried out the bending analysis of simply supported FG beam under uniformly distributed load and commented on the influence of the neutral surface position on deflection. Jia et al. [8] studied the size effect on the free vibration of functionally graded micro-beams under the combined electrostatic force, temperature change and Casimir force based on Euler-Bernoulli beam theory and von Kármán geometric nonlinearity. Ebrahimi and Salari [9] represented the thermal effect on buckling and free vibrational behaviours of functionally graded (FG) size-dependent Timoshenko nanobeams loaded to an in-plane thermal loading by implementing a Navier type solution. Pradhan and Chakraverty [10] implemented classical and first order shear deformation beam theories to investigate the free vibration behavior of functionally graded material (FGM) beams which is subjected to different sets of boundary conditions. Static bending and elastic buckling analysis of shear deformable functionally graded porous beams on the assumption of Timoshenko beam theory were performed by Chen et al. [11]. Ebrahimi and Zia [12] investigated the large-amplitude nonlinear vibration behaviours of functionally graded (FG) Timoshenko beams made of porous material. Hemmatnezhad et al. [13] investigated the large-amplitude free vibration analysis of functionally graded beams by means of a finite element formulation, where Von-Karman type nonlinear strain-displacement relationships were emloyed.
As already stated, vast body of literature exists on FG beam which is graded in thickness direction. However, analysis of the axial functionally graded (AFG) structural element is the relatively newer domain, where there exists a few papers on different aspects of such components. Li et al. [14] analyzed the free vibration of axially inhomogeneous exponentially graded beams. A simple approach to calculate the critical buckling loads of beams with arbitrary axial inhomogeneity was developed by Huang and Luo [15]. Akgöz and Civalek [16] investigated vibration response of non-homogenous and non-uniform microbeams in conjunction with Bernoulli-Euler beam and modified couple stress theory. Shahba et al. [17] carried out free vibration and stability analysis of axially functionally graded tapered Timoshenko beams. The authors employed finite element technique to analyze the effects of taper ratio, elastic constraint, attached mass and the material nonhomogeneity on the natural frequencies and critical buckling load. Huang et al. [18] presented a new approach for investigating the vibration behaviors of axially functionally graded Timoshenko beams with non-uniform cross-section. Shahba and Rajasekaran [19] presented the free vibration and stability of axially functionally graded tapered Euler-Bernoulli beams through solving the governing differential equations of motion. Zeighampour and Beni [20] carried out the vibration of axially functionally graded material (AFGM) nanobeam by using strain gradient theory. Kien [21] performed the large displacement response of tapered cantilever beams made of axially functionally graded material by the finite element method. Free vibration analysis on axially functionally graded (AFG) tapered slender beams under different boundary conditions were represented through energy principle by Kumar et al. [22] and Sarkar and Ganguli [23].
Boundary conditions of a real physical system are seldom classical, hence simulation of beams with elastically restrained ends or beams resting on elastic foundations are important domains of research. Beam on elastic foundation is used to model a lot of engineering problems and has wide application in bio-mechanics, road, rail road, geo-technics and marine engineering. Extensive research work is performed on FG beam with classical boundary. So far, a few literatures are available on beams supported or resting on elastic foundation. Yas and Samadi [24] dealt with free vibrations and buckling analysis of nano-composite Timoshenko beams reinforced by single-walled carbon nanotubes (SWCNTs) resting on an elastic foundation. Simsek and Cansız [25] represented the dynamic responses of an elastically connected double-functionally graded beam system (DFGBS) carrying a moving harmonic load at a constant speed by using Euler-Bernoulli beam theory. Murin et al. [26] formulated a differential equation of the homogenized functionally graded material (FGM) beam deflection and its solution which used in free vibration analysis of the beams with polynomial continuous longitudinal and transversal variation of material properties. The FGM beams are considered to be resting on longitudinal variable (Winkler) elastic foundation. Simsek and Reddy [27] based on the modified couple stress theory (MCST), proposed a unified higher 10 IFSL Volume 9 order beam theory which contains various beam theories as special cases for buckling of a functionally graded (FG) micro-beam embedded in elastic Pasternak medium. Akgöz and Civalek [28] performed the bending response of non-homogenous micro-beams embedded in an elastic medium based on modified strain gradient elasticity theory in conjunctions with various beam theories. Kanani et al. [29] studied the large amplitude free and forced vibration of FG beam resting on nonlinear elastic foundation containing shearing layer and cubic nonlinearity. Niknam and Aghdam [30] made an attempt to obtain a closed form solution for both natural frequency and buckling load of nonlocal FG beams resting on nonlinear elastic foundation. Literature review reveals that the field of free vibration study of depth-wise functionally graded beams is explored comprehensively, while the AFG is the newer domain. AFG beam on elastic foundation incorporate a higher level of complexity. Hence, the present study is taken up with the objective of analyzing the large amplitude free vibration problem of axially functionally graded (AFG) tapered beams on elastic foundation with exponential taper profiles. Variation of material properties (elastic modulus and density) along the length of the beam is considered according to different functions. Effect of variation of foundation stiffness on the dynamic behaviour is also studied. The large amplitude free vibration behavior is presented as backbone curves in nondimensional amplitude-frequency plane and in some particular case the mode shape results are furnished.

Mathematical formulation
The present mathematical formulation is subdivided into two parts. A static solution is obtained first under external transverse loading and it is followed by a linear eigenvalue problem based on the previously obtained static solution. Formulation for both the static and dynamic part is based on energy methods, where, the basic principle is to extremize the total energy of the system in its equilibrium condition.
In the present study, an axially functionally graded tapered beam of length L is considered as shown in fig 1(a). The beam is further considered to be rested on an elastic foundation with different end conditions such as CC, CF, SS and CS respectively as shown in fig. 1(b). For schematic representation the figure shows the taper profile of the beam as linear. However, for the present work, exponential variation of thickness has been considered along axial direction and it is expressed through the equation given below.
Here t 0 is the thickness of the beam at the root of the beam and α is the taper parameter. ξ is the normalized axial coordinate, where normalization has been carried out by length of the beam, i.e., ξ = x/L. It is worth mentioning that the present formulation and solution methodology is capable of handling any taper pattern (both along thickness as well as width) mathematically expressible as function of the axial coordinate. The thickness of the beam is considered to be small compared to lengths, so the effect of shear deformation and rotary inertia may be neglected. So the Euler-Bernoulli hypothesis is considered. The beam material properties are considered to be graded in the axial direction (modulus of elasticity, E(x), and mass density, ρ(x)) according to different gradation functions. The elastic foundation is idealized as a series of parallel massless linear spring of same stiffness value (K).  Large displacement is taken into consideration for mathematical formulation and as a result geometric nonlinearity is incorporated into the system. In deriving the strain energy expressions nonlinear strain-displacement relations are used. It is well known that for large displacement analysis of beam both bending and stretching effects are taken into account. So, the strain energy expression for the system includes strain energy stored due to bending (U b ) and stretching effect (U m ) along with strain energy stored in the foundation (

IFSL Volume 9
Incorporating the strain-displacement relations the nonlinear strain energy of the system under transverse loading can be expressed as, Here, w and u are transverse and in-plane displacements of mid-plane of the beam. E(x), I(x) and A(x) are elastic modulus, moment of inertia about the mid-plane of the beam cross-section and cross-sectional area of the beam respectively. These parameters are functions of the axial coordinate and hence, introduce added complexity to the present formulation. A separate list of symbols used for the present formulation is provided as Appendix 1 at the end of the paper. The potential energy (V) due to externally applied transverse load P(x) is given by, Although the above energy functionals (Eq. (3) and (4)) are expressed with respect to dimensional axial coordinate (x), all computations are performed in non-dimensional domain (ξ).
The static problem is formulated from the minimization total potential energy of the system, which can be written as, δ(U + V) = 0.
(5) The energy functionals (U and V) can be evaluated from the assumed transverse and in-plane displacement fields. The expressions assumed for these displacements are in terms of unknown coefficients (c i ) and orthogonal admissible functions ( and  ).
The bending displacements of the beam are described by functions , while describe stretching of the mid-plane of the beam. Start functions for these displacement fields are assumed in such a way that they are continuous and differentiable within the domain and also satisfy the boundary conditions of the beam. The higher order functions are generated by using Gram-Schmidt orthogonalisation scheme from the selected start functions. Hence, an orthogonal set of kinematically admissible functions are obtained.
Substituting the appropriate expressions in Eq. (5), the governing differential equations of the system is written in matrix form as, Here, [K s ] is the stiffness matrix and {f} is the load vector. The form and elements of the matrix and vector are provided in the Appendix.
Once the solution of the static analysis is obtained, the stiffness matrix is completely known for the corresponding deflected configuration. Now, a small amplitude vibration is assumed to be taking place about the deflected position of the system. This problem is described by a standard eigenvalue analysis and the formulation is based on Hamilton principle, which is presented by, Here, δ is the variational operator and τ denotes time coordinate. T and U indicate the kinetic energy and strain energy of the system respectively. Potential energy (V) due applied load is International Frontier Science Letters Vol. 9 13 considered to be zero, as the system executes free vibration. The kinetic energy of the system is expressed as, where, and are differentiations with respect to time, ρ(x) is density of the beam material. The displacements w and u are to be assumed as approximate displacement fields, separable in time and space. They can be represented by linear combination of unknown coefficient (d i ) and orthogonal admissible function, and  , which are identical to those used for the static analysis.
Here, ω is the natural frequency of the system. nw and nu are number of functions for w and u respectively. d i is a new set of unknown coefficients that represent the eigenvector in the matrix form. Substituting the appropriate relations in Eq. (8) gives the set of governing equations in the matrix form. -

Solution Procedure
From the mathematical formulation for the static analysis it is clear that the elements of the stiffness matrix [K s ] are functions of the unknown parameter, c i . Hence the governing equation (Eq. (7)) is nonlinear in nature and cannot be solved directly. To solve the set of equations an iterative numerical technique is introduced and a direct substitution technique with successive relaxation scheme is utilized. The solution steps are given below, Step 1: The input parameters i.e. appropriate load, allowable error limit and the relaxation parameter are defined. It is assumed that the unknown coefficients have zero value at the initial stage.
Step 2: Utilizing the input parameter the stiffness matrix for bending and stretching (and thus the total stiffness matrix) are generated, along with the load vector. As the initial values for c i s are zero the initial stiffness matrix corresponds to a linear situation.
Step 3: A new set of unknown coefficients are evaluated by inversion of the stiffness matrix and subsequent pre-multiplication with the load vector. {c} (n + 1) = [K s ({c} (n) )] -1 {f}, where n represents the iteration counter.
Step 4: The calculated set of unknown coefficients is compared with the ones from previous iteration and if the error is above the predefined allowable error limit, the next iteration is performed with modified values of unknown coefficients. The modification takes into account a relaxation parameter to predict the guess for the next iteration.
Step 5: The modified coefficient values are used to generate the modified stiffness matrix and the process is repeated till convergence is achieved, i.e., the error value falls below the predefined limit. The corresponding flowchart for the solution procedure is given in Fig. 2.
The converged stiffness matrix from the static analysis is carried into the dynamic analysis, which is a standard eigenvalue problem. The solution to Eq. (11) is obtained by developing a Matlab code that utilizes the subroutine 'eig'. The square root of the eigenvalues gives the natural 14 IFSL Volume 9 frequency of the system at the deflected configuration, whereas the eigenvectors associated with these eigenvalues can be processed to plot the modeshapes of the vibrating system.

Results and discussions
In the present analysis, it is considered that the beam is supported by four different boundary condition, namely, CC, CS, SS, CF. These are the four possible combinations arising out of three different end conditions of the beam, namely, Clamped (C), Simply supported (S) and Free (F) end. For selecting the start functions for the transverse displacement (w) the above mentioned flexural boundary conditions are used. Similarly, the start function for axial displacement (u) is selected on the basis of the in-plane boundary conditions, which is taken as zero at the ends of the beam. These start functions for both transverse and in-plane displacements are shown in Table 1. Gram-Schmidt orthogonalization scheme is used to generate the higher order functions and the number of functions is taken as 8 for each displacement. The number of Gauss points to be used for generation of results is taken as 24.    In the present analysis the thickness of the beam is considered to be tapered with an exponential profile, as shown in Eq. (1). The value of the taper parameter (α) is kept constant throughout the analysis and the numerical value for the parameter is 0.510826. The material of the beam is taken to be axially functionally graded, which implies material properties such as elastic modulus and density vary along the axis of the beam. Three different material models are chosen for the paper as shown in Table 2. It is to be pointed out that Material 1 is in fact homogeneous in nature (constant elastic modulus and density) and it is included in the results for comparison purpose.
The AFG tapered beam is supported on elastic foundation, which is idealized as series of parallel linear spring of same stiffness. Five different values of the spring stiffness are considered here and are given in Table 3. The geometric dimensions and the material properties of the beam are given in Table 4. These numerical values are utilized to generate the results for the analysis. In the present work, emphasis is on investigating the dynamic behaviour of the system corresponding to changes in flexural boundary condition and spring stiffness.
The present formulation and solution technique is validated by comparison with established results already available in literature. A clamped-clamped homogeneous uniform beam is considered and the resulted backbone curve for fundamental mode is validated with the results published by Gupta et al. [31] as shown in Fig. 3. It should be mentioned here that to comply with the system analyzed by Gupta et al. [31] the spring stiffness is considered as 0 N/m. The two sets of results show similar trends and the values show satisfactory matching, hence, establishing the validity of the present method.

16
IFSL Volume 9  Natural frequencies for the first two modes (ω 1 and ω 2 ) for the beams are provided in Table 5 corresponding to various possible combinations of boundary conditions, spring stiffness and material models. It is observed that with the increase of stiffness of the foundation the natural frequencies in all cases increase. This is due to the fact that natural frequency is the function of mass and stiffness. While stiffness of the foundation is increased, overall stiffness of the system increases but the mass of the overall system remain same. It is also observed that for particular taper parameter, boundary condition and foundation stiffness value, the natural frequency is the highest for Material 1 and lowest for Material 2. It is also observed that for particular taper parameter, foundation stiffness value and material model, the natural frequency is the lowest for CF boundary condition as the system is less stiff under CF boundary condition and highest for CC boundary condition.  In order to study the effect of boundary conditions and material model on the large amplitude behaviour of AFG beams on elastic foundation, backbone curves for first, second, third and fourth mode are furnished for three different material models in Fig. 4-6. In this particular case taper parameter (α) is considered to be fixed at 0.510826 while, stiffness (K) is 10,000 N/m. The figures are plotted in non-dimensional frequency amplitude plane, where the ordinate is dimensionless amplitude (w max /t 0 ) and abscissa is normalized frequency (ω nl /ω l ). In the present study (w max /t 0 ) is taken as 2.0 for all cases.

18
IFSL Volume 9 International Frontier Science Letters Vol. 9 The fundamental frequency (ω l ) used to normalize the nonlinear frequencies are taken from Tables 5. For all the cases, stiffness of the beam increases with increasing load due to geometric nonlinearity present in the system. This increased stiffness causes increase in free vibration frequencies with increase in deflection of the beam, as can be observed from any of the figures. So, hardening type nonlinear behaviour is exhibited by the system for all combinations boundary conditions. First three mode shape of the system corresponding to three material models are provided in Fig.  7-9. Each of these figures consists of four different sub-plots for the boundary conditions. It is also worth pointing out that amplitude of vibration has an effect on the modeshape of the system. To study this aspect in more detail, two modeshape plots corresponding to linear (w max /t 0 = 0) and nonlinear (w max /t 0 = 2) frequencies are given for each of the vibration modes. It should also be noted that the amplitude of vibration for all the plots is normalized by the corresponding maximum deflection. It was observed that difference in linear and nonlinear mode shapes increase when the boundary condition changes from CC to CF, CS and SS due to the decreasing rigidity at the boundary. However no considerable change in the mode shapes could be identified for the different stiffness of the foundation perhaps due to normalization of the maximum displacement as shown in Fig. 10.

Conclusions
In the present analysis, large amplitude free vibration behaviour of axially functionally graded thin exponential taper beam with various end condition and material gradation is investigated. The beam is further assumed to be on elastic foundation and subjected to uniformly distributed load. Energy principle is applied for the mathematical formulation for both static and dynamic analysis. The methodology is general in nature as it can be applied for other type of material gradation and taper pattern. The obtained results are validated from previously published results and were found to be in good agreement. The results for the natural frequencies are tabutated for different combinations of boundary condition, foundation stiffness and material model. Results pertaining to various material types and flexural boundary conditions for a fixed taper parameter and foundation stiffness are furnished as backbone curve for the first four modes. For all combinations of the system parameters hardening type of nonlinearity is observed. The results for the mode shape plots are also furnished.