A New Beam Model for Simulation of the Mechanical Behaviour of Variable Thickness Functionally Graded Material Beams Based on Modified First Order Shear Deformation Theory

There are many beam models to simulate the variable thickness functionally graded material (FGM) beam, each model has advantages and disadvantages in computer aided engineering of the mechanical behavior of this beam. In this work, a new model of beam is presented to study the mechanical static bending, free vibration, and buckling behavior of the variable thickness functionally graded material beams. The formulations are based on modified first order shear deformation theory and interpolating polynomials. This new beam model is free of shear-locking for both thick and thin beams, is easy to apply in computation, and has efficiency in simulating the variable thickness beams. The effects of some parameters, such as the power-law material index, degree of non-uniformity index, and the length-to-height ratio, on the mechanical behavior of the variable thickness FGM beam are considered.


Introduction
Beams are used widely in engineering fields, such as aerospace, mechanical, and civil engineering, and so on. Researchers have proposed various analytical and numerical methods to analyze static bending, free vibration, and buckling of beams. Chakrabortyet al. [1] developed a new beam element model based on first order shear deformation theory (FSDT) for static bending, free vibration of thermo-elastic FGM beams. Khan et al. [2] simulated the static bending and free vibration analysis of FGM beams using the finite element method (FEM) and zig-zag theory. Wang et al. [3] and Lee et al. [4] applied Euler-Bernoulli beam theory to analyze the free vibration of FGM beams. The nonlinear bending of a two-directionally FGM beam was considered by Li et al. [5] using the generalized differential quadrature method (GDQM). Jafari et al. [6] used analytical approximation for nonlinear vibration analysis of Euler-Bernoulli beams. Alshorbagy et al. [7] used FEM to analyze free vibration of a Euler-Bernoulli beam made of FGM. Oreh et al. [8] employed FEM for stability and free vibration analysis of a Timoshenko beam. Li et al. [9] applied an analytical solution to investigate relations between buckling loads of functionally graded Timoshenko and homogeneous Euler-Bernoulli beams. Farhatnia et al. [10] used first order shear deformation theory and GDQM for buckling analysis of FGM provide difficulties in large strain analysis; for such an analysis, the displacement/pressure (or u/p formulation) elements are more reliable and effective [37,39].
FGM is microscopically inhomogeneous composites fabricated from a mixture of ceramics and metals. FGM plates, beams, and shells can be considered as a multi-coated thin-walled composite mechanical components structures, making it possible to obtain specific mechanical properties, such as superior stiffness-to-weight ratio, high strength, and high damping. The damping efficiency of the coating systems was discussed by some researchers [40][41][42][43][44][45]. In these works, the theoretical and experimental results show that the coating structure obtained the best balance between the strength and the damping capacity, which is important for reducing the mechanical vibrations. Due to its outstanding features, the FGM has been extensively used in engineering applications, including nuclear power plants, spacecraft, turbine engines, and automotive and biomechanical applications [46][47][48][49][50][51]. Because of an increasing use of the FGM for many engineering applications nowadays, further studies for mechanical behaviours of structures made of FGM are necessary.
The main purpose of this contribution is to develop a new beam model based on modified first order shear deformation theory. The proposed beam model is simple in its formulation and free of shear locking without reduced integration. The proposed beam model is used for static bending, free vibration, and buckling analysis of a variable thickness beam made of FGM to demonstrate its accuracy and effectiveness.
The organization of this paper is as follows. In Section 2, a brief review of first order shear deformation theory is presented. In Section 3, a modified first order shear deformation for the beam model is given. Next, in Section 4, based on modified first order shear deformation and interpolating polynomials, a new beam model is developed. In addition, the element stiffness matrix, element mass matrix, and element force vector are constructed. Section 5 focuses on verification and numerical analysis for static bending, free vibration, and buckling responses of variable thickness FGM beams. In addition, the effects of some parameters on the mechanical behavior of the FGM beam are investigated. In the conclusion section, some discussions and highlights of the proposed beam model are given.

Basic Equations of First Order Shear Deformation Theory of a Timoshenko Beam
The axial displacement, u, and the transverse displacement, w, at any point based on FSDT of a Timoshenko beam are given by: u(x, y, z) = −zψ(x) w(x, y) = w 0 (x) (1) where w 0 is the transverse displacement at the middle axial beam, and ψ is the angle of cross-section rotation.
In small deformation, we have: The bending moment resultant and shear force are taken by: in which: is flexural rigidity and shear rigidity, respectively, E is the Young's modulus and G = E/[2(1 + ν)] is the shear modulus, ν is Poison's ratio, k is the shear correction factor, and I and A are, respectively, the second moment of area and cross-section area of the beam. Equilibrium of moment and transverse forces leads to: Substituting Equations (3) and (4) into (5) results in:

Modification First Order Shear Deformation for the New Beam Model
Assuming that the total deflection consists of bending deflection and a contribution of transverse shear, the angles of the cross-section slope are the result of pure bending angles and shear angles as shown in Figure 1: the subscripts, b and s, denote the bending and shear deflection, respectively.

Modification First Order Shear Deformation for the New Beam Model
Assuming that the total deflection consists of bending deflection and a contribution of transverse shear, the angles of the cross-section slope are the result of pure bending angles and shear angles as shown in Figure 1: the subscripts, b and s, denote the bending and shear deflection, respectively. By substituting Equation (7) into Equation (6) it is possible to separate the variables of two different displacement fields: Equation ( bs w w w = + = , that solution is trivial since  is an arbitrary function. Therefore, Equation (8) and (9) have a realistic solution when: By substituting Equation (7) into Equation (6) it is possible to separate the variables of two different displacement fields: Equation (8) and Equation (9) have three unknown components, w b , w s , and θ, are satisfied if w b = −θ and w s = θ. That implies w b = −w s and w = w b + w s = 0, that solution is trivial since θ is an arbitrary function. Therefore, Equation (8) and (9) have a realistic solution when: Then: Substituting Equation (11) into Equation (9), the total deflection is:

The Variable Thickness FGM Beam Element Model
The variable thickness FGM beam with length, L, as depicted in Figure 2 is considered. Then: Substituting Equation (11) into Equation (9), the total deflection is:

The Variable Thickness FGM Beam Element Model
The variable thickness FGM beam with length, L , as depicted in Figure 2 is considered.
Modeling of variable thickness FGM beam.
The height of the beam varies along the axial of beam: The material of the beam is made of two partial materials, those are metal and ceramic. The ratio of values of materials assumes that it varies through the z -direction with the power-law [46][47][48][49][50]: (14) where , cm VV are, respectively, the volume fraction of the ceramic and metal, p is the gradient index of the volume fraction, and () h h x = is the thickness of the beam. The subscripts, , cm , denote the ceramic and metal constituents, respectively.
The effective material properties, , P such as Young's modulus, , E and mass density,  are expressed using the rule of mixture [46][47][48][49][50]: The Poisson's ratio,  , is assumed constant in this study. A two-node beam element with two degrees of freedom per node is expressed here. The bending deflection may be expressed as follows: bb w = Pa (16) where a is the coefficient vector with unknown terms i a of the approximation polynomial, and b P is bending polynomial, which may be demonstrated as: The height of the beam varies along the axial of beam: The material of the beam is made of two partial materials, those are metal and ceramic. The ratio of values of materials assumes that it varies through the z-direction with the power-law [46][47][48][49][50]: where V c , V m are, respectively, the volume fraction of the ceramic and metal, p is the gradient index of the volume fraction, and h = h(x) is the thickness of the beam. The subscripts, c, m, denote the ceramic and metal constituents, respectively. The effective material properties, P, such as Young's modulus, E, and mass density, ρ are expressed using the rule of mixture [46][47][48][49][50]: The Poisson's ratio, ν, is assumed constant in this study. A two-node beam element with two degrees of freedom per node is expressed here. The bending deflection may be expressed as follows: where a is the coefficient vector with unknown terms a i of the approximation polynomial, and P b is bending polynomial, which may be demonstrated as: where ξ = 2x−l e l e is a non-dimensional coordinate. The shear deflection is taken by: According to (11), the terms of shear polynomial are: The total deflection of the beam may be presented in the following form: The angles of rotation of the face sheet beam can be expressed as: Substituting the values of the node coordinate, ξ i into Equation (21) and Equation (22), the nodal displacement of beam element is obtained as: where δ T = w 1 ψ 1 w 2 ψ 2 and: The coefficient vector, a, can be determined from nodal displacement vector: Substituting Equation (25) into Equation (16) and (18), the bending and shear deflection can be presented as follows: The total deflection of the beam may be expressed as: where:

The Element Stiffness Matrix
The strain energy of bending deflection may be expressed as: where [κ] b κ b is the bending curvature, which can be presented in the form: Given: Then: By taking Equation (34) into account in Equation (30), the strain energy of bending is obtained as: The bending stiffness matrix of the beam element may be determined by: Substituting Equation (33) into Equation (36), yields: where: The strain energy of shear deflection may be expressed as: The shear strain vector is obtained from Equation (3) and Equation (7): Given: The strain energy of shear deflection may be expressed as: The shear stiffness matrix of the beam element is obtained as: By substituting Equation (33) into Equation (36), yields: In which: According to the account, the element stiffness matrix is: Because K e s depends on α = 4D/Sl 2 e , when the thickness of the beam is very small, then α ≈ 0 and K e s ≈ 0, and the transverse shear effects vanish, and, as a consequence, the proposed beam model is free of shear locking. For the functionally graded materials, the Poisson's ratio is smaller than 0.4, so this model has no volume locking.

The Element Mass Matrix
The kinetic energy of beam element is obtained as: According to Equation (48), the element mass matrix is obtained as: or: Given: The element mass matrix of the beam element is taken by:

The Element Geometric Stiffness Matrix
At the critical load, the beam takes the buckled form as shown in Figure 3. According to Equation (48), the element mass matrix is obtained as: (49) or: Given: The element mass matrix of the beam element is taken by:

The Element Geometric Stiffness Matrix
At the critical load, the beam takes the buckled form as shown in Figure 3.  According to Figure 3, the axial shortening of the beam can be taken as follows: The strain energy of an axial compressed load, Q , can be obtained as follows: According to Figure 3, the axial shortening of the beam can be taken as follows: The strain energy of an axial compressed load, Q, can be obtained as follows: As a result, the element geometric stiffness matrix may be expressed as follows: Given: Now, the element geometric stiffness matrix is hence given by:

The Element Load Vector
The work done by the transverse distribution load, q, is expressed as follows: The element load vector may be obtained as:

Static Bending Solution
For static bending analysis, the nodal displacements can be obtained by solving the following equation: where K, F are, respectively, the global stiffness matrix and global force vector of the beam.

Free Vibration Solution
The equation of motion for free vibration analysis of the beam is obtained as follows: where M is the global mass matrix of the beam. For free vibration analysis, assuming that: By substituting Equation (63) into Equation (62), the natural frequency is obtained by solving the following eigenvalue equation: with √ λ = ω denotes the natural frequency of the beam.

Buckling Solution
For buckling analysis and determination of the magnitude of a static compressive load that will produce beam buckling, the following eigenvalue equation will be achieved: where K g is the global geometric stiffness matrix of the beam. The lowest positive eigenvalue of this equation is the magnitude of the critical buckling load, Q cr , of the beam and the corresponding eigenvector is the deformed shape of the buckled beam.

Verification
To verify the proposed beam model, in this section, the comparison of static deflection of an isotropic beam with variable thickness and an FGM beam with constant thickness is investigated.
Firstly, a simple-simple supported (S-S) isotropic beam with length, L, and variable thickness is considered here, the beam is under uniform distribution load, q, and the material properties of the beam are E = 20.83 Msi, G = 3.71 Msi, ν = 0.44. The thickness of the beam varies along the x-direction by the following formula [36]: where h 0 is the height at the mid-span of the beam, λ is a small parameter, and n is the degree of non-uniformity. The non-dimensional mid-span deflection of the beam is calculated by the following formula [35]: The comparison of non-dimensional mid-span deflection of an S-S supported isotropic beam using the proposed beam model with the results of Zenkour [35] and numerical results using Abaqus FE software packages (SIMULIA, Zaltbommel, Netherlands) are given in Table 1. Secondly, an Al/Al 2 O 3 beam with a constant thickness under a uniform distribution load, q, is investigated. The material properties of aluminum (as metal) and alumina (as ceramic) are [14]: The non-dimensional deflection at the central point of the beam for different values of the length-to-height ratio, L/h, is given by [14]: where E m is the Young's modulus of metal. The comparison of non-dimensional deflection at the central point of the beam using the proposed beam model with the results of Thai et al. [14] (using the sinusoidal beam theory (SBT) and the hyperbolic beam theory (HBT)) and the results of Li et al. [20] (using analytical solutions) are listed in Table 2. According to the two above comparisons, it can be observed that the values obtained using the proposed beam model are in good agreement with the published data.

Numerical Results for Static Bending Analysis
In this section, we apply this model to explore the static bending response of a variable thickness FGM beam with length, L, the height at the left-hand end, h 0 , subject to uniform distribution load, q. The boundary conditions of the beam are simple-simple supported (S-S) and clamped-clamped supported (C-C).
where n is the degree of non-uniformity. The maximum non-dimensional transverse deflection of the beam is expressed by: The maximum non-dimensional transverse deflection of the beam with different values of index, p = 0, 0.5, 1, 2, 5, 10, index n = 0, 0.5, 1, 2 (n = 0 corresponding to the constant thickness beam), and the length-to-height ratio, L/h 0 = 10, 20, 50, 100, is listed in Table 3, Figures 4-7.
To illustrate the effects of the power-law index, p, and the degree of non-uniformity, n, on the bending response of a variable thickness FGM beam subjected to uniform load, the non-dimensional transverse deflections of the beam are given in Table 3 and plotted in Figures 4-7.
According to Table 3 and Figures 4-7, it shows that the maximum non-dimensional transverse deflection of the beam increases as a function of the power-law index, p. It means that the richer metal FGM beam is more flexible than the richer ceramic FGM beam. The maximum non-dimensional transverse deflection of the beam increases rapidly when the power-law index, p, increases in the range of 0 ÷ 1. In addition, when the length-to-height ratio, 0 / Lh , increases, the transverse deflection of the beam increases. The deflection affected by the boundary condition on the static bending response of the beam is shown in Figures 4 and 5. In general, the deformation shape of the variable thickness beam is not symmetric. The transverse deflection of the C-C supported beam is smaller than the one of the S-S supported beam. Figure 8 shows the distribution of the axial normal stress, x  , and shear stress, xy  , across the z -direction at the mid-span. It shows that when    The influence of index n on the maximum non-dimensional transverse deflection of the beam is shown in Table 3, Figures 4-7. By increasing index n, the maximum non-dimensional transverse deflection of the beam increases. When index n increases in the range of 0 ÷ 1, the transverse deflection of the beam increases rapidly. When n = 0 and n = ∞, the heights of the beam are constant, h(x) = h 0 and h(x) = h 0 /2, so that the maximum deflection position appears at the mid-span of the beam.
In addition, when the length-to-height ratio, L/h 0 , increases, the transverse deflection of the beam increases. The deflection affected by the boundary condition on the static bending response of the beam is shown in Figures 4 and 5. In general, the deformation shape of the variable thickness beam is not symmetric. The transverse deflection of the C-C supported beam is smaller than the one of the S-S supported beam.        Figure 8 shows the distribution of the axial normal stress, σ x , and shear stress, τ xy , across the z-direction at the mid-span. It shows that when p = 0, the axial normal stress, σ x , has linear distribution while the shear stress, τ xz , is constant along the height of the beam.

Verification
To confirm the accuracy of the proposed beam model, a comparison of frequencies of a cantilever isotropic tapered beam with length, L and different values of the taper ratio, c = 1 − h 1 /h 0 , is considered herein. The non-dimensional natural frequency of the beam is defined as [34]: where m 0 , I 0 are the mass per unit length and the flexural rigidity at the left-hand end of the beam, respectively. The comparison of non-dimensional natural frequencies using the proposed beam model with the results of Banerjee et al. [34] and the numerical results using Abaqus FE software packages is shown in Table 4. A (S-S) Al/Al 2 O 3 beam with constant thickness and different values of the ratio of L/h is considered. The material properties of Al (metal) and Al 2 O 3 (ceramic) are [14]: The non-dimensional natural frequencies of the beam are defined as [14]: Table 5 shows the first three non-dimensional natural frequencies of the beam using the proposed beam model and the results of Thai et al. [14] (using SBT and HBT). The present frequencies are in good agreement with the other published results.

Numerical Results of Free Vibration Analysis
In this section, the variable thickness beam, which is given in the static bending problem, is employed again herein. The non-dimensional frequencies are obtained as: The results of the free vibration analysis of the variable thickness FGM beam with different values of the length-to-height ratio, L/h 0 , power-law index, p, index n, and boundary conditions are reported in Tables 6 and 7 the stiffness of the FGM beam, which leads to a decrease in the non-dimensional frequencies of the beam. This is due to the fact that higher values of index p correspond to a high portion of the metal in comparison with the ceramic portion, and, consequently, the FGM beam becomes more flexible. According to Tables 6 and 7, and Figure 10, when increasing the index, n , the non-dimensional frequencies of the beam decrease. In fact, the fundamental non-dimensional frequency of the beam decreases rapidly when the index, n , increases in the range of 0 2.  Figure 11 plots the first four mode shapes of the S-S support variable thickness FGM beam with 0 / 10 Lh= , 0.5 p = , and 2. n = Figure 12 plots the first mode shapes of the variable thickness FGM beam with 0 / 10, 0.5 L h p == , and 2 n = for some boundary conditions. It can be seen that the mode shapes of the variable thickness FGM beam show greater differences from those of the constant thickness FGM beam. The amplitude of vibration at the slender end is higher than that at the other end.   Table 6. First three non-dimensional natural frequencies, ω i , i = 1, 2, 3, of S-S and C-C supported FGM beam with L/h 0 = 10 depending on index p and index n.

Verification
To verify the proposed beam model for buckling analysis, S-S and C-C isotropic beams with a constant thickness are considered. The length of the beam is L = 1, and the modulus of elasticity and Poisson's ratio of material are, respectively, 3

10
, The present buckling load results are compared with the analytical solutions (which are given in [8]), the results of Ferreira [52] and Oreh et al. [8], and numerical results using Abaqus FE software packages and are finally reported in Table 8.

Verification
To verify the proposed beam model for buckling analysis, S-S and C-C isotropic beams with a constant thickness are considered. The length of the beam is L = 1, and the modulus of elasticity and Poisson's ratio of material are, respectively, 3

10
, The present buckling load results are compared with the analytical solutions (which are given in [8]), the results of Ferreira [52] and Oreh et al. [8], and numerical results using Abaqus FE software packages and are finally reported in Table 8.  By considering Tables 6 and 7, the non-dimensional frequencies of the C-C beam are higher than the ones of the S-S beam. Furthermore, the non-dimensional frequencies decrease as a function of the length-to-height ratio, L/h 0 .
The influences of index p and the degree of non-uniformly, n, on the free vibration of the beam are considered here. Table 6 and Figure 9 show that the index, p, has a strong effect on the non-dimensional frequencies of the beam. It can be seen that increasing the index, p, will reduce the stiffness of the FGM beam, which leads to a decrease in the non-dimensional frequencies of the beam. This is due to the fact that higher values of index p correspond to a high portion of the metal in comparison with the ceramic portion, and, consequently, the FGM beam becomes more flexible. According to Tables 6 and 7, and Figure 10, when increasing the index, n, the non-dimensional frequencies of the beam decrease. In fact, the fundamental non-dimensional frequency of the beam decreases rapidly when the index, n, increases in the range of 0 ÷ 2. Figure 11 plots the first four mode shapes of the S-S support variable thickness FGM beam with L/h 0 = 10, p = 0.5, and n = 2. Figure 12 plots the first mode shapes of the variable thickness FGM beam with L/h 0 = 10, p = 0.5, and n = 2 for some boundary conditions. It can be seen that the mode shapes of the variable thickness FGM beam show greater differences from those of the constant thickness FGM beam. The amplitude of vibration at the slender end is higher than that at the other end.

Verification
To verify the proposed beam model for buckling analysis, S-S and C-C isotropic beams with a constant thickness are considered. The length of the beam is L = 1, and the modulus of elasticity and Poisson's ratio of material are, respectively, E = 10 3 GPa,ν = 0.333. Q cr = Q cr The present buckling load results are compared with the analytical solutions (which are given in [8]), the results of Ferreira [52] and Oreh et al. [8], and numerical results using Abaqus FE software packages and are finally reported in Table 8. Furthermore, an Al/Al 2 O 3 beam with a constant thickness subjected to a uniaxial load is investigated here. The material properties of the two components of the FGM beam are [9]: The length-to-height ratio is L/h = 5 and L/h = 10, and the boundary condition of the beam is S-S supported. The non-dimensional buckling load is defined as [9]: Table 9 shows the comparison between the present non-dimensional buckling load and the results of Li et al. [9] using the analytical solution. Table 9. Comparison of the nondimensional buckling load of the constant thickness FGM beam with the ratio of L/h = 5 and L/h = 10 for the S-S and C-C support condition.  The two comparisons indicate a good agreement between the present results by using this model and published results.

Numerical Results for Buckling Analysis
The variable thickness FGM beam, which is given in the static bending problem, is again considered herein. The beam is subject to a uniaxial compressed load, Q.
The non-dimensional critical load is obtained by: The effect of the power-law index, p, on the critical load of the variable thickness FGM beam is given in Table 10 and Figure 13. It implies that the power-law index, p, strongly affects the critical buckling load of the FGM beam in some range. When the power-law index, p, is increasing in the range of 0 ÷ 2, the critical load decreases rapidly, while p > 2, the index, p, has a slight effect. Table 10. Non-dimensional critical load, Q cr , of an S-S and C-C variable thickness FGM beam depending on the ratio of L/h 0 and index n.  In Table 10 and Figure 14, the critical load is a function of the degree of non-uniformity, n. In general, the critical load decreases when increasing index n. In addition, we can see again in Table 10 that the critical load of the C-C supported variable thickness beam is much higher than the S-S supported beam.
In Table 10 and Figure 14, the critical load is a function of the degree of non-uniformity, .
n In general, the critical load decreases when increasing index .
n In addition, we can see again in Table   10 that the critical load of the C-C supported variable thickness beam is much higher than the S-S supported beam. Figures 15 and 16 show the first buckling mode shapes of variable thickness FGM beam with

Conclusions
The paper established the modified first order shear deformation beam theory, and applied this theory to model the variable thickness FGM beam element. Using the proposed beam model, the static bending, free vibration, and buckling of the variable thickness FGM beam were investigated and simulated. This beam model can effectively simulate the mechanical behavior of FGM beams by comparing it with results of other methods and the numerical simulations using Abaqus FE software packages. Parameters that study the influence of geometry, materials, and boundary conditions on the mechanical behavior of the beam were considered.
The advantage of the proposed beam model is its simplicity in formulation, and high convergence and free shear-locking without reduced integration or selective reduced integration. This proposed beam model can be applied widely to simulate the mechanical behavior of FGM beams.