Finite Element Modelling for Static and Free Vibration Response of Functionally Graded Beam

A 1D Finite Element model for static response and free vibration analysis of functionally graded material (FGM) beam is presented in this work. The FE model is based on efficient zig-zag theory (ZIGT) with two noded beam element having four degrees of freedom at each node. Linear interpolation is used for the axial displacement and cubic hermite interpolation is used for the deflection. Out of a large variety of FGM systems available, Al/SiC and Ni/Al2O3 metal/ceramic FGM system has been chosen. Modified rule of mixture (MROM) is used to calculate the young’s modulus and rule of mixture (ROM) is used to calculate density and poisson’s ratio of FGM beam at any point. The MATLAB code based on 1D FE zigzag theory for FGM elastic beams is developed. A 2D FE model for the same elastic FGM beam has been developed using ABAQUS software. An 8-node biquadratic plane stress quadrilateral type element is used for modeling in ABAQUS. Three different end conditions namely simply-supported, cantilever and clampedclamped are considered. The deflection, normal stress and shear stress has been reported for various models used. Eigen Value problem using subspace iteration method is solved to obtain un-damped natural frequencies and the corresponding mode shapes. The results predicted by the 1D FE model have been compared with the 2D FE results and the results present in open literature. This proves the correctness of the model. Finally, mode shapes have also been plotted for various FGM systems.


Latin American Journal of Solids and
high-transverse shear stresses occur at the interface during fabrication and operation.This leads to delamination at the interface and poor load-bearing performance.Developments in materials engineering lead to a new type of composites with smoothly and continuously varying thermos-mechanical properties that are called functionally graded materials (FGMs).Functionally graded materials (FGMs) provide an elegant solution to this problem, wherein the abrupt change in composition of material from one layer to another is replaced by layers of gradually varying microstructure and composition.FGM structures have various advantages over the composite laminates, such as smaller stress concentrations, smaller thermal residual stresses and the possibility to achieve specific for different applications.This concept has found many potential applications like thermal and corrosion barriers, medical implants, lightweight armour material with high-ballistic efficiency, etc. FGMs are made by combining different materials through complex processes such as powder metallurgy methods, physical and chemical vapor deposition, solid freeform fabrication etc. Detailed process for manufacturing FGM's and various areas of application are presented by Rasheedat M. et al.
Among the various models which have been proposed to predict the effective elastic properties of two-phase materials, the self-consistent model of Hill (1965), the Mori and Tanaka (1973), the linear rule of mixtures and in Jones (1999) the intermediate rule of mixtures provide simple and convenient ways for predicting the overall response.For kinematic modelling, several one-dimensional (1D) beam and 2D plate theories such as Cho and Oden (2000) applied CLT (Classical Laminate Theory) for thermal stress analysis of infinite functionally graded beams and plates with layer-wise compositional change, Nguyen et al. (2008) found the effect of transverse shear stresses by energy equivalence using FGM models based on FSDT theory.The transverse shear stresses are obtained using the membrane stresses and various equilibrium equations and Sina et al. (2009) developed a new beam theory different from the traditional first-order shear deformation beam theory is used to analyze free vibration of functionally graded beams.The beam properties is varied through the thickness according to simple power law distribution in terms of volume fraction of material constituents.It is assumed that the lateral normal stress is zero and the governing equations of motion are derived using Hamilton's principle.Reddy (1983) developed a higher order shear deformation theory of composites.This theory accounts for the parabolic distribution of the transverse strains through the thickness of the structure.This theory has been found to be predicting deflection and stresses more accurately than FSDT.Cheng and Batra (2000) exploited Reddy's third-order plate theory to study buckling and steady state vibrations of a simply supported functionally gradient isotropic polygonal plate resting on a Winkler Pasternak elastic foundation and subjected to uniform in-plane hydrostatic loads.Sankar (2001) presented an elasticity solution for simply supported FG beams subjected to sinusoidal transverse loading.Kapuria et al. (2008) have presented an efficient zigzag model and its experimental validation for thermo-elastic static and free vibration response of Al/SiC and Ni/Al2O3 FGM beams.Amal et al. (2011) gave free vibration characteristics of functionally graded beam with material graduation in axially or transversally through the thickness based on the power law and the obtained results are compared with previously published work.Nuttawit and Variddhi (2012) applied the differential transformation method (DTM) to investigate free vibration of functionally graded beams supported by arbitrary boundary conditionsand the material properties of beams are assumed to obey the power law distribution.H. Yaghoobi et al. (2010) investigated Bending analysis of a functionally graded simply supported beam for varied neutral surface subjected to a uniformly distributed load.Structures 13 (2016) 690-714 Results indicate that position of neutral surface is very important in functionally graded materials.Kadoli et al. (2008) implemented displacement field based on higher order shear deformation theory to study the static behavior of functionally graded metal-ceramic (FGM) beams under ambient temperature.FGM beams with variation of volume fraction of metal or ceramic based on power law exponent are considered.Najeeb and Alam (2012) presented a one dimensional finite element model using an efficient layerwise (zigzag) theory for the dynamic analysis of laminated beams integrated with piezoelectric sensors and actuators.Mohanty et al. (2012) presented the evaluation of static and dynamic behavior of functionally graded ordinary (FGO) beam and functionally graded sandwich (FGSW) beam for pined-pined end condition.The variation of material properties along the thickness is assumed to follow exponential and power law.A finite element method is used assuming first order shear deformation theory for the analysis.Furqan and Naushad (2013) assessed higher order theory of laminated beams under static mechanical loads.The Third order theory and First order shear deformation theory are assessed by comparison with the exact two-dimensional elasticity solution of the simply-supported beam.Mehta et al. (2013) used finite element method in modelling the dynamic behavior of FGM to determine its natural frequency.The properties in the functionally graded material are assumed to vary according to power law.The natural frequencies were obtained for FG beams under various boundary conditions including Clamped-Fixed, Simply supported-Fixed, Clamped-Clamped, Simply supported-Simply-supported, and Clamped-Simply supported.Shi-rong Li et al. (2014) studied the free vibration of functionally graded material (FGM) beams based on both the classical and the first-order shear deformation beam theories.Nguyen et al. (2014) presented the analytical solutions for the static analysis of the transversely or axially functionally graded beams with tapered cross-section.The elastic modulus of the beam varies according to the power form.

Latin American Journal of Solids and
The present study deals in developing a 1D FE model based on Zigzag theory for FGM elastic beams using MATLAB as a mathematical tool.The model is used to compute the deflection, stress and shear stress at various values of side to thickness ratio and upto first five fundamental frequency.A 2D FE model for the same elastic FGM beam has been developed using ABAQUS software.The results predicted by the 1D FE model have been compared with the results obtained using 2D FE model.These predicted results have also been compared with the results presented in the open literature.This is done to prove the correctness of the FE model developed.A three layer Al/SiC beam, five layer Ni/Al2O3 and ten layer Ni/Al2O3 beam is modelled for the analysis purpose.Mode shapes has also been plotted for these FGM system under various boundary conditions viz.simply supported, clamped-clamped and cantilever.

MATERIAL MODELLING
For the bi-material FGM systems, the volume fractions Vc and Vm of the ceramic and the metal as taken by Kapuria et. al (2006) are assumed to vary along z-direction according to following power law: Latin American Journal of Solids and Structures 13 (2016) 690-714 Where M is a non-negative real number called the inhomogeneity parameter.The variation of Vm across z is plotted for M = 4.The volume fraction for a layer is computed at the center of the respective layer.
The effective material properties for a layer are computed using different averaging methods based on volume fraction at the centre of layer.The effective material elastic modulus E of a layer is computed using the modified rule of mixtures (MROM), which was originally proposed for cemented carbides by Tomota et al (1976) and later adopted for FGM by many researchers Cho and Oden (2000).
According to this approach, the two phase material is treated as an isotropic composite for which uniaxial stress σ, strain ɛ, are related to constituent stresses , andɛ , ɛ as, Where q is Stress to Strain transfer ratio between two phases.The value of q for Al/SiC has been experimentally determined by Kapuria et. al (2008) to be 91.6 GPa, which is used in this work for computing the elastic modulus of Al/SiC FGM samples.For Ni/Al2O3 system, the value of q is taken as 4.5 GPa as used by Kapuria et. al (2008).Now the Young's modulus E of the FGM material can be calculated explicitly using MROM as, For this work, other material property i.e.Poisson's ratio ν and density ⍴ have been calculated using linear rule of mixture (ROM) as: Latin American Journal of Solids and Structures 13 (2016) 690-714 Three Functionally graded beams are considered in this paper, namely beam (a), beam (b) and beam (c).All three beams are combination of metal and ceramic with bottom of the beam being metal rich and top of the beam being ceramic rich.Beam (a) is a three layered fgm beam of Al/SiC and the volume fraction each layer is pre-defined as shown in figure 2. Beam (b) is a 5 layered fgm beam of Ni/Al2O3 and in this beam also the pre-defined volume fraction is shown in the figure 2. Beam (c) is a 10 layered fgm beam of Ni/Al2O3 and is bit different from the other two beams.In beam (c) each layer is of equal thickness i.e. 0.1h and the bottom most layer is of pure metal.Moreover, the volume fraction of the remaining nine layers is computed using equation 1. Young's modulus and Density of each layer is calculated using MROM and ROM respectively.

ZIG-ZAG THEORY FOR FGM BEAM
Consider a functionally graded rectangular solid beam.The width, thickness, and length of FGM beam are denoted by b, h, and a, respectively.The FGM substrate is modelled as a laminate of a number of perfectly bonded isotropic layers with different material properties, which vary as per given Latin American Journal of Solids and Structures 13 (2016) 690-714 gradient of volume fractions of the constituents.The axes along the length, width, and thickness are x, y, z, respectively.The loads given to the beam do not vary along thickness direction.The z-co-ordinates of the bottom surface of the k-th layer (numbered from bottom) is denoted as .The reference plane z=0 either passes through or is the bottom surface of the layer.For beam of small width, assumptions for mathematical simplifications of 1-D model are- Transverse normal stress neglected ( =0) (iii) Axial and transverse displacement considered independent of y.
For ZIGT, w is approximated by integrating the constitutive equations for ɛ by neglecting the contribution of .Thus, constitutive equation for ɛ is According to the basic assumption of the zigzag theory, For ZIGT the axial displacement u is assumed as a combination of a third order variation across the laminate thickness with a layerwise linear variation.So, for layer, Where and denote the translation and rotation variables of layer Now, Solids and Structures 13 (2016) 690-714

Equation of Motion
The Hamilton's principle is represented as below, Let , be the normal forces per unit area on the bottom and top surfaces of the beam in direction z.The extended Hamilton's principle for the elastic beam can be expressed as, using notation And integrating the term by parts for now a beam of length a under above mentioned loading condition.The principle equation given by Eq. ( 12) reduces to the form  This variational equation is expressed in terms of δ 0 u , δ 0 w , δ 0  to yield equations of motions and boundary conditions.
In the above equation, u and w can be expressed as with The inertia terms can be expressed as Latin American Journal of Solids and Structures 13 (2016) 690-714 with elements explicitly given as, [ 11 , The strain increments δ x  and zx  are related to virtual displacements by, with The strain energy term in variational equation Eq. ( 13) expressed as Where beam stress resultants are defined as The terms due to load on bottom and top surfaces are expressed as Latin American Journal of Solids and Structures 13 (2016) 690-714 At last, the terms due to load at the end of the beam are expressed as, Where the superscript * represents the mean value at the ends and x V = < zx  > Substituting all above calculated sub-parts in the original equation Eq. ( 13), we have

FINITE ELEMENT MODEL FOR FGM BEAM
A two noded element, having four degrees of freedom at each node, finite element model is developed for the analysis of elastic layered FGM beams under various kinds of mechanical loading employing 1-D zigzag theory presented earlier.
The stress resultants x N , x M , x P , x Q in eq. ( 24) can be related to 0 u , 0 w , 0  as: where, if A represents the beam stiffness matrix, and is given by, A = Now, the beam stress resultants from Eq. ( 23) can be expressed as below, The elements of above matrices are: Latin American Journal of Solids and Structures 13 (2016) 690-714 We can define stresses , strains , and stress resultant in a condensed form and call them as generalised stress ,generalised strain and generalised stress resultant respectively as, (32) also, and

Interpolation Function for FG Beam Element
As described above, a two noded element, having four degrees of freedom at each node, has been used for interpolating the mechanical variables for zigzag theory.Since the highest derivative of 0 u , 0 w , 0  in the variational equation are 0,x u , 0, xx w , 0 , x  , therefore to meet the convergence requirement of finite element method , 0 u and 0  must be C 0 continuous and 0 w should be C 1 -continuous at the element boundaries.Thus, 0 u and 0  have been interpolated using Lagrange interpolation function and 0 w has been interpolated using Hermite cubic interpolation function.Also, it is evident that as zx  is dependent on 0  only and not on 0 , x w , this interpolation scheme will not lead to shear locking.
Denoting the values of an entity (.) at node i by (.)i, 0 u , 0 w , 0  are interpolated in an element of length a as Where the elemental variable matrices 0 e u , 0 e w , 0 e  are given by The interpolation function N and N are given by Latin American Journal of Solids and Structures 13 (2016) 690-714 where, Now if we define the generalized displacement vector at element level as The generalized displacement vector û and generalized strain vector  can be expressed in terms of e U as: where,

System Equations and Boundary Conditions
Summation of all the elements of e T to the variational integral, the system equation is obtained as In which the M, K, P matrices are assembled from their elemental matrices respectively.In case of point loads applied at nodes are added to P, at locations corresponding to their degrees of freedom numbers.
The variationally consistent boundary conditions at the beam ends are obtained as, The geometric boundary conditions for various end conditions are as follows, Simply supported: x N = 0, (movable) 0 u = 0 (immovable) 0 w = 0, x M = 0, x P = 0 Clamped: Free: Latin American Journal of Solids and Structures 13 (2016) 690-714

ABAQUS MODELLING
The structure of various beams that are used are shown in Figure 2. The thickness of the beam is taken as unity and the length of the beam is taken corresponding to side-to thickness ratio.The global x-coordinate is taken along the length of the beam and the global y-coordinate is taken along the thickness.The model is developed using 20 elements along the axial direction and 1 elements in each layer along the thickness direction.Three different boundary conditions are taken into account viz.simply-supported, cantilever and clamped-clamped.The beam is analyzed for deflections and stresses under various boundary condition when the beam is subjected to load working along the Y-direction for various side-to-thickness ratios (a/h).The figure below shows the meshed model of beam (c).The center deflection and stresses are presented here in non-dimensional form.

RESULTS AND DISCUSSION
The validation of zigzag theory FE model developed is presented for static and free vibration response of FGM beams.The material properties and stress to strain transfer ratio for Al /SiC and Ni /Al2O3 systems is taken up from Kapuria et al. (2008).This ratio is used to predict effective elastic modulus of both the FGM systems using MROM.The static analysis and free vibration response of layered FGM beams, with ceramic content varying from 0 to 40%, has been presented here.The material properties of the basic constituents Al, SiC, Ni and Al2O3 of the FGM beams, considered for computing the effective properties of the layers, as used by Kapuria et. al (2008) are listed in Table 1.The young's modulus of the different layers of the beam are computed using MROM while density of the different layers of the beam are computed using ROM.An assessment of 1D Zigzag model for the static response of elastic functionally graded beams made of multiple layers of isotropic materials of layer-wise constant composition is presented in this chapter.Structures 13 (2016) 690-714 The following mechanical load cases are considered:

Latin American Journal of Solids and
1. Uniformly distributed load at the top surface of the beam (load case 1) Sinusoidally varying load over the top surface of the beam (load case 2) where x gives the location of any point on the beam and 'a' is the total span of the beam.The 0 p is a constant with magnitude 1.
The results that are presented have been non-dimensionalised using following expressions with, S = a/h, 0 Y = 127.8295GPa (For Ni/Al2O3) = 88.509GPa (For Al/SiC) The dimensionless entities are chosen such that their values are almost independent of S for different beams.

Static Response
The 2D FE results for w , x  , and zx  at points across the thickness are given in Table 2 and Table 3 for elastic beam (c) and for inhomogeneity parameter M = 4. Table 2 shows the result for load case 2 for span to thickness ratio S = 5, 10 and 40.The load is applied at the top surface of the beam.
The non-dimensionalized deflection, w is reported at the top surface.The observed error for simply supported beam (c) with load case 2 on comparison with Kapuria et al. (2006) is within 3% for thick beam with S = 5.For moderately thick beam, S = 10 error is within 7%.The 1D FE based on zig-zag theory and 2D FE Abaqus results are presented in Table 3 for beam (c) with cantilever and clamped-clamped end conditions for aspect ratios S = 5, 10 and 40 and inhomogeneity parameter M = 4. Table shows the results in bending for both the load cases i.e. sinusoidally varying and uniformly distributed load.Further, it can be seen from the table that the results are in close agreement with each other, which eventually shows the correctness of the FE model presented.
The Tables 4 and Table 5 present the results for 3 layered beam (a) with Al/SiC FGM system under sinusoidally varying load and uniformly distributed load for different boundary conditions.The inhomogeneity constant for each case is M = 4.A different FGM system is considered in Tables 6 and Table 7.The results are presented for 5 layered Ni / Al2O3 FGM beam (b) under sinusoidally varying load and uniformly distributed load on the top surface for different end conditions.The length to thickness ratio for both the load cases is taken as S = 10 with inhomogeneity parameter M = 4. Structures 13 (2016) 690-714 Finally, it can be observed that the results obtained using 1D FE model are in good agreement with the results obtained using 2-D FE abaqus, which proves the correctness of the 1-D zig-zag theory based FE model.

Free Vibration Response
The convergent FE results for natural frequencies are obtained by modelling the FGM beams for both of the considered FGM systems viz.Al/SiC and Ni/Al2O3 and for different boundary conditions.The predicted natural frequencies, n  of first few modes (using MROM and linear ROM based property estimates) obtained from 1D FE zigzag model are compared with 2D FE results.All the frequency values are given with unit Hz.An 8-node biquadratic plane stress quadrilateral type element, CPS8R is used for modeling in ABAQUS.Along the beam length 80 elements and 15 elements were seeded along thickness direction.The results for natural frequencies are validated with results presented in Kapuria et al. (2008).These results are presented in Tables 8 to Table 9.It can be seen that present beam model in conjunction with MROM with the value of q as 91.6 GPa (Kapuria et al. (2008)) predicts fundamental natural frequencies of 3 layer Al/SiC cantilever beams very accurately with maximum error of 1.35% and for clamped-clamped boundary condition with maximum error of 5.60%.
In case of 5 layer Ni/Al2O3 FGM system, beam (b) has been modelled using 1D FE as well as 2D FE.Same type of discretisation element (CPS8R) has been used with 80 elements along the length and 25 elements along thickness direction.The results have been found to be in close agreement with each other.The fundamental frequencies are compared with results of Kapuria et al. (2008) and found to have maximum error of 2.86% in case of cantilever and 12.9% maximum error in case of beam clamped at both the ends.Here the q value is taken to be 4.5 GPa (Kapuria et al. (2008)).Close agreement of the theoretical model results with the results presented in Kapuria et al. (2008) proves consistency and robustness of 1D FE model and thus using it, the fundamental frequencies of the two foresaid FGM systems for simply supported boundary conditions are presented in Table 11.From Table 11 it can be seen that the natural frequencies obtained using 1D FE model are in good agreement with 2D FE results for all three different beams.Also, a 10 layer Ni/Al2O3 beam model has been considered for frequency analysis.The results for 10 layer Ni/Al2O3 beam are presented in Table 10.All the 1D FE results are in good agreement with 2D FE results obtained using ABAQUS.Table 12 shows sensitivity analysis i.e. the effect of dividing the FG beams into different number of layers on natural frequencies using the developed 1-D FE model.Non-dimensional natural frequencies up to third mode of vibration are computed for span-to-thickness ratios 5 and 20.The beam used by Thai and Vo (2012) is divided into 8-layers, 10-layers and 12-layers of equal thickness and the results are compared with Thai and Vo (2012).The material properties and boundary conditions are taken from Thai and Vo (2012).Young's modulus and density for each layer are calculated using the linear rule of mixture at the center of each layer.Poisson's ratio is assumed to be constant.The inhomogeneity parameter (M) is taken as 5.The same results are also plotted in Figures 14 and15.It is evident from the graph that with the increase in number of layers the results approach to the exact solution presented by Thai and Vo (2012).Moreover, it is interesting to see that the variation in frequency for different layers increases slightly with increase in the mode of vibration

CONCLUSIONS
The developed 1D FE model yield fairly accurate results for static as well as free vibration response of functionally graded (FG) beam.In the assessment for static response, generally the results of 1D FE zig-zag theory based FE model and 2D FE ABAQUS are very close for the transverse deflection but for the stresses some fluctuation in the results is observed.In case of free vibration, natural frequencies and mode shapes have been obtained from both 1D FE as well as 2D FE sources, for different end constraints.The natural frequencies of all the FGM beam systems have been found to be in good agreement with the results from literature, which proves the correctness of the developed models.Mode shapes have been determined using ABAQUS and MATLAB.These are very accurate and almost coincide with each other.The natural frequencies in general increase with the increase of modes of vibration.Also that the material constitution has negligible effect on the general deflected shape of the beam under a particular end condition.The Sensitivity analysis in terms of number of Latin American Journal of Solids and Structures 13 (2016) 690-714 layers of FG beam shows that the computed results approach to the exact result as the number of layers are increased.FGM beam can also be modelled as layer wise using the developed 1D FE model.

Figure 1 :
Figure 1: Variation of volume fraction of metal across the non-dimensional thickness.

Figure 2 :
Figure 2: Geometry and Volume fraction variation through thickness for the FGM beams.
expressions for N and u f in the integral for P e , it yields, convenience, these elemental degrees of freedom have been arranged like below,
The mode shapes for mid-surface deflection for the first three flexural modes (n = 1, 2, 3) have been plotted using 2D FE model.The 2D FE mode shapes have been plotted using the values that are normalized with the maximum deflection values.They have also been obtained by using zigzag 1D FE analysis and are compared with the 2D FE results.They have been presented in Fig 5 to Fig 13 for moderately thick beam (S = 10), respectively for different boundary conditions.This comparison shows their close agreement with each other.Therefore they can be trustfully used for future references.

Figure 5 :
Figure 5: First three flexural mode shapes of 3 layer Al/SiC FGM beam for cantilever end condition.

Figure 6 :Figure 7 :
Figure 6: First Three Flexural mode shapes of 3 layered Al/SiC FGM beam for simply supported end condition.

Figure 8 :
Figure 8: First three flexural mode shapes of 5 layer Ni/Al2O3 FGM beam for cantilever end condition.

Figure 11 :
Figure 11: First three flexural mode shapes of 10 layer Ni/Al2O3 FGM beam for cantilever end condition.

Figure 12 :Figure 13 :
Figure 12: First three flexural mode shapes of 10 layer Ni/Al2O3 FGM beam for simply-supported end condition.

Figure 14 :
Figure 14: Sensitivity analysis of frequency in terms of number of layers for span to thickness ratio, S=5 of FG beam.

Figure 15 :
Figure 15: Sensitivity analysis of frequency in terms of number of layers for span to thickness ratio, S=20 of FG beam.

Table 1 :
Constituent Properties of materials used.

Table 2 :
Static response of 10 layer FGM beam (c) for simply supported end condition with M=4.

Table 4 :
Static response of 3 layer Al/SiC FGM beam (a) for simply supported end condition with M=4.

Table 5 :
Static response of 3 layer Al/SiC FGM beam (a) with M=4.

Table 6 :
Static response of FGM beam (b) for simply supported end condition with M=4.

Table 8 :
Natural Frequencies of 3 layer Al / SiC FGM beam under various end conditions for M =4.

Table 9 :
Natural Frequencies of 5 layer Ni / Al2O3 FGM beam under various end condition for M = 4.

Table 10 :
Natural Frequencies of 10 layer Ni / Al2O3 FGM beam under various end condition for M = 4.

Table 11 :
Natural Frequencies of simply supported FGM beams for M = 4.

Table 12 :
Non-dimensional natural frequencies of simply supported FGM beams for M = 5.