## Abstract

The large amplitude vibrations of functionally graded (FG) beams consisting of metal rich layers at the bottom, ceramic rich layers at the top, and a concentrated mass at the mid-span have been studied using coupled displacement field method. Unlike traditional methods, the coupled displacement field method reduces the 2n undetermined coefficients problem, one each for total rotation and transverse displacement distribution of the beam at n modes, to n undetermined coefficients using a coupling equation obtained from the minimization of potential energy principle. A suitable admissible function having single undetermined coefficient has been assumed for total rotation distribution and the corresponding transverse displacement distribution of the beam has been obtained at each mode for hinged-hinged and clamped-clamped FG beams. The equations of motion for large amplitude vibrations of FG beams at each mode in terms of the undetermined coefficients are derived from the conservation of total energy principle. The free vibration problem is solved using harmonic balance method whereas the forced vibration problem is solved using the Newmark-β method to obtain the time response of the undetermined coefficients and the dynamic response of the beam has been computed from the modal superposition method. The proposed coupled displacement field approach has been successfully applied for the first time to study the large amplitude vibrations of FG beams with suitable validations, and the influence of power law index, slenderness ratio, harmonic load, and concentrated mass has been investigated.

## 1 Introduction

Functionally graded materials (FGMs) are advanced type of composite materials which are made up of two or more inorganic materials such as metals, ceramics, and polymers. The material properties of FGMs are varied in two or more desired directions by varying the volume fractions of their material constituents. They are generally preferred over traditional composites due to their ability to withstand high temperatures and thermal stresses, unlike traditional composites which undergo de-lamination and de-bonding under high thermal stresses. The concept of FGMs was first introduced in mid 1980s by Japanese scientists and since then, their applications have been considered in wide range of engineering fields such as automobile, aerospace, biomedical, defense, and aircraft to eliminate the disadvantages of traditional composites. Due to the wide range of applications of FGMs, it is important to analyze the FGM structures and their dynamic behavior in different mechanical loading and boundary conditions.

The works related to the dynamic behavior of FGM structures are extensively reported in the literature since the introduction of FGMs. The initial works on the dynamic behavior of FGM structures were carried out by Loy et al. [1] and Pradhan et al. [2] on the vibration characteristics of functionally graded (FG) cylindrical shells for different material distributions and boundary conditions. Kim [3] carried out the temperature dependent vibration analysis of FG rectangular plate. Zhao et al. [4] used element free *kp*-Ritz method to carry out the free vibration analysis of FG plates. Baferani et al. [5] developed an exact analytical method for the study of free vibrations of thin rectangular FG plates based on classical plate theory. Taheri et al. [6] proposed an improved isogeometric analysis approach for the study of free vibration characteristics of FG structures. These are few of the many works related to the dynamic behavior of FG cylindrical shells and plates.

However, the dynamic behavior of FG beams is of particular interest in the present study and many related works are available in the literature. Aydogdu and Taskin [7] performed the free vibration analysis of simply supported FG beams using different higher order shear deformation theories and classical beam theories. Alshorbagy et al. [8] used finite element method to study the free vibration behavior of FG Euler–Bernoulli beam under different boundary conditions. Şimşek and Reddy [9] used a new higher order beam theory and the modified couple stress theory to study the bending and free vibration of functionally graded micro-beams. Celebi and Tutuncu [10] employed an exact plane elasticity approach for free vibration analysis of FG beams to obtain exact natural frequencies.

Simsek and Kocaturk [11] studied the free and forced vibrations FG beam having simply supported edges subjected to moving harmonic loads. Azadi [12] carried out the forced vibration analysis of FG beam using finite element method considering the temperature dependency of the material properties and damping effects. Nguyen and Bui [13] developed a finite element method for the dynamic analysis of FG beam subjected to harmonic loads in thermal environment using higher order hierarchical Timoshenko beam element. Zibdeh and Hannah [14] developed an analytical method to obtain the dynamic response of simply supported FG beam under the assumptions of Euler–Bernoulli beam theory. However, in all these works, the FG beam has been assumed to undergo linear vibration behavior. Although it is true for small amplitude vibration, at large amplitudes, the FG beam exhibit non-linear vibration behavior.

There are few works reported in the literature on the non-linear free and forced vibration analysis of the FG beam to study the large amplitude vibrations. Simsek [15] presented the non-linear vibration behavior of simply supported FG Timoshenko beam due to moving harmonic loads using von-Kármán’s non-linear strain-displacement relationship. The Lagrange’s equations were used to obtain the non-linear equations of motion and were solved using Newmark-β method combined with direct iteration method. Hemmatnezhad et al. [16] developed a Timoshenko beam theory based finite element formulation to study the large amplitude free vibration behavior of FG beams using von-Kármán type non-linear strain-displacement relationships. The non-linear free vibration analysis was carried out by Hosseini et al. [17] for size dependent FG micro-beams using higher order non-linear governing equations derived from Hamilton’s principle. The non-linear free vibration analysis was carried out for FG beams using different shear deformation theories by Xie et al. [18].

Rao et al. [19,20] studied the large amplitude free vibration behavior of both clamped-clamped and hinged-hinged isotropic beams using the coupled displacement field method. The geometric nonlinearity of the beams due to the effect of tension developed in the beam at large amplitudes given by Woinowsky-Krieger [21] was considered in these studies. The free vibration equation of motion was derived from the principle of conservation of total energy and was solved using harmonic balance method. The forced vibration analysis of non-linear systems has been performed using the principle of modal superposition by Kuran and Özgüven [22] and Ferhatoğlu et al. [23,24]. The aforementioned works are based on the non-linear vibrations of the homogenous material systems. However, modeling of FG beam using coupled displacement method in conjunction with harmonic balance method and model superposition to investigate large amplitude free and forced vibration behavior of the beam respectively is rarely found in the literature.

Herein, the study of non-linear vibrations of the Timoshenko FG beam at large amplitudes having concentrated mass at the center has been carried out for the first time using coupled displacement field. The hinged-hinged and clamped-clamped FG beams have been modeled using coupled displacement field method and the equation of motion for each mode has been derived by applying the principle of conservation of total energy. The harmonic balance method has been used to solve the free vibration equations of motion. The forced vibration equations of motion of each mode have been solved using the Newmark method and the modal superposition method has been applied to obtain the dynamic response of the system.

## 2 Functionally Graded Beam Material Gradation

An FG beam based on power law consists of pure metal at the bottom, pure ceramic at the top, and material gradation occurring between these layers in the direction of the beam thickness as shown in Fig. 1 has been considered.

*v*

_{c}is given by power law as

*h*is the beam height, $z\xaf$ is the distance of any point on the beam from the mid-plane, and

*k*is the power law index. The sum of the volume fractions of metal and ceramic is equal to one. Therefore, the volume fraction of the metal

*v*

_{m}can be expressed as

*P*of any given layer of the FG beam can be expressed as

*P*

_{m}and

*P*

_{c}represent the material properties of the metal and ceramic, respectively. Substituting the metal and ceramic volume fraction in the aforementioned equation, the equation can be rewritten as

## 3 Methodology

The large amplitude free and forced vibration analysis of FG beam having concentrated mass at the center has been considered in the present study. The methodology followed is summarized in flowchart as shown in Fig. 2. The kinematic relations of FG beam based on Timoshenko beam theory have been considered and the coupling equation has been obtained by applying the principle of minimization of potential energy on the FG beam. The admissible function for beam rotation distribution has been assumed and the corresponding transverse displacement distribution of the FG beam has been derived for each mode using coupled displacement field method for both hinged-hinged and clamped-clamped boundary conditions. Using the expressions obtained by coupled displacement field method, the strain and kinetic energy of the FG beam, work done by the tension developed in the beam, and the external forces acting on the beam have been computed for each mode. The equations of motion for large amplitude vibration of FG beam have been derived for each mode from the principle of conservation of total energy of the FG beam. The free vibration analysis of FG beam using harmonic balance method and forced vibration analysis of FG beam using Newmark method in conjunction with modal superposition method have been performed.

## 4 Coupling Equations

The kinematics of the FG beam based on Timoshenko beam theory has been considered and expressions for axial strain, shear strain, and strain energy have been derived. The principle of minimization of potential energy has been applied to the FG beam under static lateral load distribution to obtain the static equilibrium equation which couples the total rotation and transverse displacement distribution of FG beam.

*x*along the beam axis and distance

*z*from the neutral axis of the beam,

*w*is the transverse displacement,

*ψ*is the bending rotation, $\beta \xaf$ is the shear rotation,

*θ*is the total rotation, and

*x*,

*z*are the independent spatial variables which define any given generic point on the beam. The distance of the neutral axis from the mid-plane of the beam represented in Fig. 1 is obtained as

*e*is the distance between the neutral axis and mid-plane of the FG beam. If

*e*is positive, the neutral axis is above the mid-plane and is below when

*e*is negative. The expressions for axial strain

*ɛ*

_{x}and shear strain

*γ*

_{xz}of the beam are

The expression for strain energy *U* of the FG beam and expression for work done *W* by static lateral load distribution *p*(*x*) acting on the FG beam are given as follows:

*a*)

*b*)

*E*is Young’s modulus,

*G*is shear modulus,

*κ*is the shear correction factor that accounts for parabolic shear distribution,

*ν*is the Poisson’s ratio,

*b*is the beam width,

*h*is the beam height, and

*A*is the cross-sectional area of the beam.

## 5 Coupled Displacement Field Method

An FG beam consisting of concentrated point mass at the center has been considered for two different cases of boundary conditions: hinged-hinged (H-H) and clamped-clamped (C-C) at both the ends of the beam as shown in Fig. 3. The suitable admissible functions are assumed for beam rotation distribution of H-H and C-C functionally graded beams at each mode of vibration and corresponding coupled displacement fields for transverse displacement at each mode are obtained using the coupling equation (Eq. (15)).

A single admissible function with one undetermined coefficient for total rotation distribution *θ* of the FG beam at *i*th mode which satisfies the symmetric and essential boundary conditions at both the ends of the beam has been assumed for both H-H and C-C conditions as follows:

*a*

_{i}represents the undetermined coefficient proportional to maximum transverse displacement of the FG beam of length

*L*at

*i*th mode.

*θ*in the aforementioned equation and applying the boundary conditions for transverse displacement

*w*as

The coupled displacement field for transverse displacement distribution *w* and the transverse amplitude *A*_{i} of the FG beam at *i*th mode can be obtained for both hinged and clamped boundary conditions as follows:

The aforementioned method is known as coupled displacement field method which reduces two undetermined coefficients per mode to single undetermined coefficient problem *a*_{i} by coupling the transverse displacement *w*_{i} distribution with beam rotation *θ*_{i} distribution.

## 6 Equation of Motion

The admissible functions assumed for beam rotation *θ*_{i} and the coupled displacement fields obtained for transverse displacement *w*_{i} of the FG beam for each mode are substituted in the expressions for the strain energy, kinetic energy, and work done by tension developed in the FG beam. The principle of conservation of total energy has been applied and the equations of motion for the large amplitude vibrations of FG beam at each mode have been obtained for both C-C and H-H functionally graded beams.

The expressions for strain energy (*U*), kinetic energy (*T*), and work done by the tension (*W*^{T}) developed in the FG beam are listed as follows:

*m*is mass per unit length,

*J*

_{m}is the diametric moment of inertia per unit length,

*ρ*is the density of the FG beam, and

*M*is the concentrated mass located at the center of the FG beam.

*J*

_{m}=

*mr*

^{2}in Eq. (25), where

*r*is the radius of gyration of the FG beam, the kinetic energy of the FG beam can be rewritten as

*T*

_{a}is given by

*i*th mode under the influence of varying lateral force distribution

*f*(

*x*,

*t*), the below expression is obtained

*f*(

*x*,

*t*). Substituting the expression of strain energy, kinetic energy, and work done by tension developed in the FG beam in the aforementioned equation and taking

*β*=

*L*/

*r*as slenderness ratio, the equation has been obtained as

*a*)

*b*)

*c*)

*a*)

*b*)

*c*)

*i*th mode as

## 7 Large Amplitude Free Vibrations

*f*(

*x*,

*t*) = 0 as

*a*

_{i}=

*a*

_{i,m}sin(

*ω*

_{i,NL}

*t*), where

*ω*

_{i,NL}is the non-linear natural radian frequency at the given mode, and hence we get

^{3}(

*ω*

_{i,NL}

*t*) in the aforementioned equation can be written as

*ω*

_{i,NL}in the aforementioned equation and substituting the expression in Eq. (37), non-linear natural radian frequency at amplitude ratio

*q*

_{i,m}=

*A*

_{i,m}/Θ ($\Theta =I/A$, where

*I*is area moment of inertia and

*A*is cross-sectional area) can be obtained as

## 8 Large Amplitude Forced Vibrations

*F*(

*t*) =

*F*

_{o}sinΩ

*t*has been considered to be acting at location

*x*=

*ξ*(constant) on the FG beam as shown in Fig. 4. For the case of moving load, the location (

*ξ*) of the load acting on the beam is a function dependent of time and can be defined accordingly. The concentrated load

*F*(

*t*) can be expressed as the distributed load

*f*(

*x*,

*t*) with the help of the popular Dirac Delta function

*δ*(

*x*) as

*a*

_{i}for $\u2200i\u2208[1,n]$ for each time-step from 0 s to 1 s (taking time-step, Δ

*t*= 0.001 s). From the principle of modal superposition, the lateral displacement of the beam

*w*and total beam rotation

*θ*at any given time-step

*t*=

*t*

_{o}and any given position

*x*=

*x*

_{o}is given by [26]

The total number of modes *n* is considered such that the results thus obtained converge to required accuracy. Thus, the forced vibration response of the FG beam subjected to harmonically varying concentrated lateral load can be obtained.

## 9 Validations

A code has been developed in the python programming language to carry out the large amplitude free and forced vibration analysis of an FG beam for hinged-hinged and clamped-clamped boundary conditions. The three-step validation of the developed method and code has been performed in this section. In the first step, the non-linear to linear frequency ratios of an isotropic beam have been compared with the published results to validate the coupled displacement field method and harmonic balance method used in the present study and the corresponding code developed in python. In the second step, non-linear to linear frequency ratios of an FG beam have been compared with the published results to verify the FG material modeling and validate the accuracy of the corresponding code. In the last step, the linear and non-linear maximum mid-span displacements of FG beam subjected to a moving point load have been obtained for different magnitudes of the load and have been compared with the published results to verify the application of Newmark method and modal superposition method used in the present study to obtain the large amplitude forced vibration response of an FG beam and validate the corresponding code.

The ratio of non-linear to linear radian frequency parameters of hinged-hinged and clamped-clamped isotropic steel beam at *β* = 50 for different amplitude ratios (ratio of mode amplitude to the radius of gyration of the beam cross section) has been computed. The results are in good agreement with the published results by Rao et al. [20], as tabulated in Table 1. This confirms the correctness of coupled displacement field and harmonic balance method formulation developed for hinged-hinged and clamped-clamped boundary conditions of the beam and the corresponding code.

The ratio of non-linear to linear radian frequency parameter of hinged-hinged and clamped-clamped FG beam for different amplitude ratios (ratio of mode amplitude to the radius of gyration of the beam cross section) has been computed at power law index *k* = 1/3 and *L*/*h* = 100. The computed results along with the results available in Taeprasartsit [27] are tabulated in Table 2. Taeprasartsit [27] used finite element formulation to study the geometric nonlinearity of thin Euler–Bernoulli functionally graded beam based on von-Kármán’s non-linear strain-displacement relationship. Given the present study is based on Timoshenko beam theory which studies the geometric nonlinearity of functionally graded beam by considering the effects of tension developed in the beam given by Woinowsky-Krieger [21], the results obtained are in good agreement with the literature.

The variation of linear and non-linear maximum mid-span displacements of hinged FG beam subjected to the moving point load at velocity 25 m/s (*ξ* = 25*t*, Ω = 0) has been obtained and validated with results available in Simsek [15] as shown in Fig. 5. It can be observed from the figure that the mid-span displacements obtained from present study are slightly higher compared to those from the literature. The slight difference is because the present study considers only the bending and transverse shear effects in the FG beams and the effects of tension developed due to geometric nonlinearity of the FG beams at large displacements given by Woinowsky-Krieger[21], whereas Simsek [15] considers the extension and extension-bending coupling effects in addition to the bending and transverse shear effects in the FG beams and assumes von-Kármán’s non-linear strain-displacement relationship to consider the geometric nonlinearity of FG beams at large displacements.

## 10 Numerical Results and Discussions

A functionally graded beam having concentrated mass at the center that is made up of metal (stainless steel) and ceramic (aluminum oxide) detailed in Sec. 2 has been considered in the present analysis. The material properties of metal and ceramic used in the analysis are tabulated in Table 3. The large amplitude free and forced vibration analysis of FG beam has been performed for hinged-hinged and clamped-clamped boundary conditions varying various FG beam parameters and amplitude ratios.

Table 4 shows the first mode of non-dimensional non-linear radian frequencies of H-H and C-C FG beams at *L*/*h* = 20 for different power law indices and amplitude ratios. As the maximum transverse amplitude ratio increases, the non-dimensional non-linear radian frequency of the FG beams increases due to increase in the effect of geometric nonlinearity of FG beams at higher amplitudes of vibration. As the power law index of FG beams increases, the non-dimensional non-linear radian frequencies decrease due to decrease in stiffness and increase in mass of the FG beams. The non-linear frequencies of C-C FG beam are higher compared to H-H FG beams as expected.

The non-linear frequencies in the first row of Table 4 at amplitude ratio equal to zero are the linear frequencies which vary with power law index *k*. When the non-linear frequencies in each column are divided by corresponding linear frequency in the first row, it gives us the ratio of non-linear to linear frequencies varying with amplitude ratio. The ratio of non-linear to linear frequencies with the increase in the maximum transverse amplitude ratio at different power law indices of the H-H and C-C FG beams is plotted in Figs. 6(a) and 6(b) respectively. The ratio of non-linear to linear radian frequencies of FG beams increases as the nonlinearity of FG beams increases with increase in maximum transverse amplitude ratio. The rate of change of non-linear to linear frequency ratio with respect to the amplitude ratio increases with the increase in the amplitude ratio. It can be observed from the figures that the non-linear to linear frequency ratio increases with the increase in power law index from 0 to 0.5 and decreases with the increase in power law index from 1 to 5. This implies that the non-linear to linear frequency ratios initially increase and then decrease with the increase in the power law index of the FG beams. The ratios of non-linear to linear frequency values for the H-H FG beam are higher compared to the C-C FG beam.

The linear and non-linear vibration response at the mid-span displacements of FG beam at *k* = 1, *L*/*h* = 80, *L* = 2 m, and *b* = 0.04 m having H-H boundary conditions and subjected to harmonic load has been shown in Fig. 7. The harmonic load of magnitude 500 N and excitation frequency 50 rad/s have been applied at the mid-span of the beam. The maximum linear and non-linear mid-span displacements of FG beam have been shown in Fig. 7. It can be observed that the maximum non-linear mid-span displacement is lower compared to linear mid-span displacement of the FG beam as the stiffness increases due to the effect of geometric nonlinearity of the FG beam at large displacements.

The variation of maximum linear and non-linear mid-span transverse displacements with the magnitude of harmonic force applied has been plotted at Ω = 50, *k* = 1, *L* = 2 m, and *b* = 0.04 m for H-H and C-C FG beams as shown in Figs. 8(a) and 8(b), respectively. The linear mid-span transverse displacement of FG beams increases linearly with increase in the magnitude of harmonic load applied. The non-linear mid-span transverse displacement of FG beams deviates from the linear relationship with increase in the magnitude of harmonic load applied. As the magnitude of harmonic load applied increases, the amplitude of vibration of FG beams increases. The resistance offered to the deformation by FG beams is higher at larger displacement due to the geometric nonlinearity. Therefore, the non-linear displacements are smaller compared to linear displacements and deviation increases with increase in the magnitude of harmonic load applied.

The non-linear mid-span displacements are normalized by linear mid-span displacement obtained at the maximum harmonic load magnitude (1000 N). The variation of normalized mid-span displacements with harmonic load is plotted in Figs. 9–11. If the normalized mid-span displacement increases linearly from zero to one with increase in the magnitude of harmonic load, it implies that the effect of geometric nonlinearity in the beam due to tension is negligible. The increase in the deviation of the plot from the straight line implies that the effect of geometric nonlinearity in the beam is increased. The variation of normalized non-linear transverse displacement amplitudes at the mid-span with the magnitude of harmonic force applied has been plotted at Ω = 50, *k* = 1, *L* = 2 m, and *b* = 0.04 m for H-H and C-C FG beams for different slenderness ratios as shown in Figs. 9(a) and 9(b), respectively. The normalized non-linear mid-span displacement of FG beams decreases with increasing rate as the slenderness ratio of the FG beams increases. This implies that the effect of geometric nonlinearity of FG beams increases with the increase in the slenderness ratio. As the slenderness ratio of FG beams increases, the flexibility (less resistance to deformation) of the FG beams increases which result in higher amplitudes of vibrations, and thus the effect of geometric nonlinearity increases.

The change in nonlinearity of C-C FG beams to the variation of FG beam parameters can be best represented at *L*/*h* = 100 compared to *L*/*h* = 80 at which the deviation of non-linear from linear displacements is lower. Therefore, the H-H FG beam at *L*/*h* = 80 and C-C FG beam at *L*/*h* = 100 have been considered for further analysis in the present study. The normalized non-linear transverse displacement amplitudes at the mid-span with the increase in the magnitude of harmonic force applied have been plotted at Ω = 50, *L* = 2 m, and b = 0.04 m for H-H and C-C FG beams for different power law indices in Figs. 10(a) and 10(b) respectively. The normalized non-linear mid-span displacement of FG beams decreases with the increase in the power law index of the FG beams. This implies that the effect of geometric nonlinearity of FG beams increases with the increase in the power law index. As the power law index of FG beams increases, the stiffness of the FG beams decreases which results in higher amplitudes of vibrations at which the effect of geometric nonlinearity of FG beams is higher.

The normalized non-linear transverse displacement amplitudes at the mid-span with the increase in the magnitude of harmonic force applied have been plotted at Ω = 50, *k* = 1, *L* = 2 m, and b = 0.04 m for H-H and C-C FG beams for different concentrated masses at the mid-span as shown in Figs. 11(a) and 11(b) respectively. The normalized non-linear mid-span displacement of FG beam decreases as the concentrated mass at the mid-span of FG beams increases. This implies that the effect of geometric nonlinearity of FG beams increases with the increase in the concentrated mass. As the concentrated mass of FG beams increases, the natural frequency of the FG beams decreases which results in higher amplitudes of vibrations at which the effect of geometric nonlinearity of FG beams is greater.

## 11 Conclusions

In the present study, the large amplitude free and forced vibrations of FG beams having concentrated mass at the center have been studied using a coupled displacement field method. The coupled displacement field method reduces the number of unknowns from 2*n* to *n* undetermined coefficients and thus the overall computational time in obtaining the non-linear vibration responses of FG beams would be reduced. The FG beam rotation distribution is assumed, and the corresponding transverse displacement distributions are obtained for each mode. The equations of motion have been derived for each mode using the conservation of total energy principle and the corresponding linear and non-linear natural frequencies are obtained using harmonic balance method. The forced vibration equations of motion are solved using Newmark-*β* method and the dynamic response of the FG beam has been obtained from the principle of modal superposition.

The results of the present study indicate that the non-linear frequency of the FG beams increases as the power law index increases. The deviation of the non-linear vibration response from the linear vibration response increases with the increase in the magnitude of harmonic load applied, slenderness of the beam, concentrated mass, and power law index. However, the deviation of non-linear responses from linear responses is lower in case of C-C FG beam compared to H-H FG beam.

The future scope of the present work includes the incorporation of local effects such as extension and extension-bending coupling in the FG beam in the current formulation, the development of the current formulation for the large amplitude vibration analysis of FG rotor-bearing system, and also formulation of defects such as cracks, corrosion, and porosities in the FG beams and rotors to study their effects on the large amplitude vibration.

## Funding Data

This research received no specific grant from any funding agency in the public, commercial, or not-for-profit sectors.

## Conflict of Interest

There are no conflicts of interest. This article does not include research in which human participants were involved. Informed consent is not applicable. This article does not include any research in which animal participants were involved.

## Data Availability Statement

The authors attest that all data for this study are included in the paper.