Finite Element Analysis of Functionally Graded Beams using Different Beam Theories

The present study deals with buckling, free vibration, and bending analysis of Functionally Graded (FG) and porous FG beams based on various beam theories. Equation of motion and boundary conditions are derived from Hamilton’s principle, and the finite element method is adopted to solve problems numerically. The FG beams are graded through the thickness direction, and the material distribution is controlled by power-law volume fraction. The effects of the different values of the power-law index, porosity exponent, and different boundary conditions on bending, natural frequencies and buckling characteristics are also studied. A new function is introduced to approximate the transverse shear strain in higher-order shear deformation theory. Furthermore, shifting the position of the neutral axis is taken into account. The results obtained numerically are validated with results obtained from ANSYS and those available in the previous work. The results of this study specify the crucial role of slenderness ratio, material distribution, and porosity condition on the characteristic of FG beams. The deflection results obtained by the proposed function have a maximum of six percent difference when the results are compared with ANSYS. It also has better results in comparison with the Reddy formulae, especially when the beam becomes slender.


Introduction
A composite material is made of two or more constituent material with different mechanical properties. This new material has physical and chemical characteristics, unlike that of the individual components. In laminate composite structures, isotropic elastic layers are joined together to provide mechanical and advanced material properties. The typical problem with laminate composite is the concentration of stress at the site of separation of the different layers, which causes cracks and the delamination phenomenon. In functionally graded materials (FGMs), because the changes from one material to another are trivial, there is no delamination [1]. FGMs are categorized as composite materials that have contiguous conversion in the properties of materials from one plane to another, thus reducing the stress concentration existed in conventional composites [2]. FGMs have several potential advantages that made their use more common in comparison with laminated composites [3]. These advantages are including reducing in-plane and transverse stresses along with thickness, proportional distribution of residual stress, improved thermal properties, greater fracture and corrosion resistance, and reducing stress concentration factors [1].These features have led to their widespread use in various scientific and engineering applications, such as mechanical, structural, aerospace, nuclear, armory and, etc. FGMs are typically made of isotropic components (e.g., metals and ceramics). FGMs are also used as The present work deals with the study of static, buckling, and free vibration characteristic of unidirectional FG beams and porous FG beams. Hamilton's principle is used to obtain the equations of motion and the essential boundary conditions for different beam theories. The FEM is used for numerical solving of various FG and PFG beam problems. FGM and PFGM beams for the various parameters like length to thickness ratio, power-law and porosity index, and boundary conditions are studied. The accuracy and effectiveness of this paper are verified by a comparison between the results obtained by ANSYS software and those available in previous research. A SOLID-186 element having three degrees of freedom per node has been employed in the ANSYS software. The functionally graded material beam with a uniform variation of the material property through the thickness is estimated as a laminated section containing a number of isotropic layers. The power law is used to determine material properties in each layer. Ten by forty mesh and twenty number of layers are found to give good accuracy from convergence studies. As a part of this study, a new polynomial function is introduced to approximate the shear strain. In the previous study, the effect of the neutral axis position is not significant for various analyses. Moreover, the effect of the power-law index on the shear correction factor is not considered for different analyses, including bending, free vibration, and buckling.

Constitutive Relation
Various prevalent methods are existed methods for acquiring the effective mechanical properties of materials. These methods are the rule of mixture, the Murray-Tanaka method, and Hill's own adaptive approach [25]. In the present paper, the rule of the mixture is used to obtain the effective material properties. Where , Lh and b are the length, the height, and width of the beam, respectively.

Rule of Mixture
In Figure 1, it can be seen a general FGM beam, which is made of two different materials (e.g., steel and ceramic), in which the mechanical properties are changing smoothly through the thickness. In Figure 2-a, the effect of the various power-law index is shown. According to the power-law distribution, mechanical properties of the FG beam can be defined as: Where; , bt ff show the material properties of the beam at the bottom and the upper face of the beam. Also, the parameter n depicts the non-negative power-law index, which relates to the distribution of material properties along with the thickness of the beam. The component ,

Kinematic
Consider a beam in Figure 1, a global coordinate system is assumed, the x-coordinate coincide with the beam axis, z-coordinate is taken along the thickness, and y-coordinate is along the width of the beam. The displacement field of the beam can be expressed as follows based on different theories:

Euler-Bernoulli Beam Theory (EBBT)
The governing displacement equation of Euler-Bernoulli beam with 3 degrees of freedom per node is given by: Displacement field equation in matrix form is:  By assuming Hook's law, the strain-stress relationship describes as below: Using Hamilton's principle and applying the variational method to Equations 3 and 5, the general displacement field of the beam is:

First-order Shear Deformation Beam Theory (FSDBT)
The axial and transverse displacement equations of a FG beam with three degrees of freedom per node based on First-order shear deformation theory are expressed as: Displacement field equation in the matrix form: Using Equation 9.a, the only non-zeros strain-displacement field equations are given by: By considering Hook's law, and using Equation 9.a, the stress field equation describes as below: Where;  is the shear correction factor, which indicates the variation of shear stress thorough the beam thickness. In composite materials, the shear correction factor  is not constant and depends on both the cross-section shape and the distribution of material along with the beam thickness [7]. Herein, for the sake of brevity, the mathematical formulation of the shear correction factor is neglected; therefore, for more information, the [7] can be seen. In Table 1, some shear correction factor is calculated and shown. Where; , SK and V  are strain energy, kinetic energy, and work done by external forces, respectively. 1 t and 2 t are the initial and final times, respectively.

Higher-order Shear Deformation Beam Theory (HSDBT)
The general displacement equations of Higher-order shear deformation beam with 4 degrees of freedom per node are given by: 15.a can be rewritten in the matrix form as follows: Where; () fz is a function which approximates the shear strain through thickness. () fz  denotes the derivative of () fzwith respect to the z .
Using Equation 15.a, the strain field equation of HSDT is of the form of: By assuming Hook's law, and using Equation 16, the stress field equations describe as below: Using Hamilton's principle and applying the vibrational method to Equations 16 and 17, the general displacement field of the beam is: Where; , SK

The Position of the Neutral Axis
In FG materials, due to the variation of the mechanical properties, the position of the neutral axis is changing [11]. The effect of the power-law index and material properties is shown in Figure 3. The position of the neutral axis can be computed by using the following formulation:

Finite Element Formulation
There are various numerical methods used for solving engineering problems [26][27][28]. In the present paper, the Finite Element Method (FEM) is used to solve the governing equations. The equations of motion in the previous sections, are numerically solved by the FEM for bending, buckling, and free vibration problems. FEM formulation of each beam theories describe separately as follow:

Euler-Bernoulli Beam Theory (EBBT)
In EBBT, the Lagrange shape function approximates the in-plane displacement. Also, the Hermite cubic shape function is used for the estimation of the transverse displacement and rotation.

First-order Shear Deformation Beam Theory (FSDBT)
In FSDBT, the Lagrange shape function is used for both in-plane and transverse displacements and also rotation.

Third-order Shear Deformation Beam Theory (TSDBT)
In TSDT, the Lagrange shape function estimated the in-plane displacement, and the Hermit shape function estimated the transverse displacement and rotation, respectively.
The element Stiffness, mass, and geometric stiffness matrix, respectively can be calculated as below: The components , ee KMand e G K denote the element stiffness, mass, and geometric stiffness matrix, respectively. Here the Gaussian quadrature is used to solve the integration. Where Ng is indicate the number of gauss points. The construction of finite element analysis is demonstrated as follows: q at the center of the beam is considered. The deflection of the mid-point of the present function is compared with Reddy's third-order shear deformation theory and also the exact solution obtained in [29].

Structures of FEM code for FG beam analysis
The results in Table 2  nel indicates the number of elements used for beam mesh.

Numerical Results
In this section, FG beams are analyzed differently under various boundary conditions, including clamped-free (C-F), simply support-simply support (S-S), and clamped-clamped (C-C). For more detailed information on boundary conditions, Simsek (2010) [9] to be seen. Functionally graded material of beam is a composition of aluminum (AL) (as metal) and alumina (Al2O3) (as ceramic). Two different kinds of material are used for numerical analysis. The following material properties are used for numerical modeling. For bending and free vibration, material 1 from Ref. [16] is used. For buckling problems, material 2 from Kehya and Turan (2016) [8] is utilized.  The properties alter according to the power-law. Thus, the upper surface is pure alumina, while the bottom surface is pure aluminum. The results of buckling analysis are compared with those available in Kehya and Turan (2016) [8]. Bending and free vibration results of FG beams are compared with those obtained and normalized from ANSYS solid 186 (see Figure 4). Simply supported condition for brick element gives inaccurate results; therefore, analyzing with simply support boundary condition is neglected. The FG beams are designed by solid 186 (20 nodes 3D element) with 20 layers. The cross-section mesh is 40 10  , and the length of the beam is meshed with 30 elements (see Figure 4-a).

Static Analysis
In this section, the bending behaviour of the FG beam under point load 0 q is investigated. For this analysis, the material one is used. The boundary condition is clamped-free, simply supported, and clamped-clamped. Tables 3 to 5 contain the dimensionless maximum deflection of FG beams. It can be observed that the results obtained by CBT are underestimated. The FSDBT results are good where the beam is deep, and the shear is dominant. HSDBT results are accurate enough in a different situation. The CBT has better results where the shear stress has not prevailed (slender beam). The maximum error values for the higher-order shear deformation results are less than seven percent when compared with ANSYS results. The new proposed function has better results in a fixed-fixed situation. The deflection of the FG beam is increased, by increasing in power-law index and porosities since the stiffness of FG beams is decreased. Where; w and w are dimensionless and normal deflection of the FG beams, respectively. , and m E I L denote Young's modulus of the considered metal, second moment of inertia, and length of the FG beams, respectively. Also, 0 q indicates the point load.

Buckling analysis
In the following section (  [8].

2097
An obvious outcome here is that the results of the proposed model are approximately close to the results of Kehya and Turan (2016) [8]. Maximum differences of higher-order shear deformation results are almost five percent. By increasing the power-law index, the differences are increased. Increasing in Both power-law and porosity exponent has a significant effect on buckling loads. The buckling load decreases when the power-law and porosities are increased.

Free Vibration
The first three dimensionless frequencies  of C-F, S-S, and C-C of FG and PFG beams are illustrated in Tables 7 to 9, respectively. Various boundary conditions, length-to-depth ratio ( / ) Lhare considered. The effect of different values of the power-law index n and porosity index  on the vibration characteristic of the FG beams is investigated.
For free vibration analysis, the following eigen-value equation [30][31][32]. The results of different theories are compared with those obtained from ANSYS solid 186. The CBT overestimates the frequencies due to the neglect of shear deformation. The FSDBT loses accuracy when the FG beams become slender.
In most cases, results obtained from higher-order shear deformation theories have less than one percent differences in comparison with those obtained from ANSYS. The maximum difference in results is less than five percent. The results of the new proposed shear deformation have less difference with ANSYS when the FG beams become slender. By increasing in porosities, the FG beam vibrations decrease. An increase in the power-law index , n leads to decreasing the values of frequencies. While the power-law index tends to zero (full-ceramic), the frequency values are increased.
The first three mode shapes of FG beams are plotted in Figure 4-b for clamped-clamped FG beams of ( / 5) Lh  , 1 n  and 0.0 Where;  and  are dimensionless and normal frequency of the FG beams. , and eq eq Eh  denote equal density, equal Young's modulus of FG beams, and height of the FG beams, respectively.

Conclusions
In the present paper, static, buckling and free vibration of various FG beams by various beam theories are studied. Hamilton's principle is used for acquiring the equation of motion. Different boundary conditions, power-law index, and porosity conditions are assumed for estimating the behaviour of FG beams. The shift of the neutral axis position is taken into account in the analysis. A new polynomial function is introduced for approximating the shear strain along with the thickness. For validation of the results of the free vibration and bending results, FG beams are modeled in ANSYS (solid 186). The results specify the influences of the slenderness ratio (L/h), material distribution, porosity index, on the characteristics of the beam as follow:  The proposed shear strain function satisfies the stress-free boundary condition on the bottom and top surface of the beam. The numerical results show that significant accuracies can be reached using the new proposed shearstrain function. Therefore, the shear correction factors do not require, which is common in classic beam theories;  By increasing in slenderness ratio (L/h), the deflection, natural frequencies, and critical buckling load are increased, respectively;  Non-dimensional frequencies increase with a decrease in power-law and porosity index;  The stiffness of the FG beam is following a decreasing pattern by increasing the porosity and power-law index. Therefore, the deflection increases, while the critical buckling load is decreased; 2101  The CBT bending results is underestimated, while the buckling and free vibration results are overestimated, especially in the thick beam;  The FSDBT results are less inaccurate when the results are compared with HSDBT and ANSYS results. Due to the dependency of a shear correction factor to the power-law index, obtaining various shear correction factor is unendurable;  Both HSDBTs have accurate results when the results are compared with ANSYS results. The proposed model has less error, especially when the beam becomes slender. Therefore, the proposed function can be used in practical applications to obtain more accurate results.