Document Type : Research Paper
Authors
^{1} Department of Mechanical Engineering, Heritage Institute of Technology, Kolkata, 700107, India
^{2} College of Food Technology, IGKV, Raipur, 492012, India
^{3} Department of Mechanical Engineering, Jadavpur University, Kolkata, 700032, India
Abstract
Keywords
Main Subjects
Higher Order Approximations for Bending of FG Beams
Using BSpline Collocation Technique
^{a} Department of Mechanical Engineering, Heritage Institute of Technology, Kolkata, 700107, India
^{b} College of Food Technology, IGKV, Raipur, 492012, India
^{c} Department of Mechanical Engineering, Jadavpur University, Kolkata, 700032, India
KEYWORDS 

ABSTRACT 
Functionally graded material; Hamilton principle; Bspline collocation; Power law; Axial stress; Shear stress 
In the present study, a functionally graded cantilever beam has been analyzed to observe its deformation behavior and stress variations. While the material properties (modulus of elasticity, modulus of rigidity, and density) have been varied along the height of the beam, Poisson’s ratio has been considered a constant. The governing equations have been derived using Hamilton’s Principle in the framework of higherorder beam theory. The derived equations are then simplified to a single equation, which is similar to that of isotropic beams. However, the work is extended to include a few higherorder terms and to study the effect of the incorporation of these terms on the resulting FG beam behavior. The development of a single governing equation for studying the statics and dynamics of an FG beam with the incorporation of higherorder terms is a unique part of the report. The solution of the governing equation is approached using approximate methods; in this work, the Bspline collocation technique is used to arrive at the results. A sixthorder Bspline basis function is used as an approximating polynomial, and the Greville abscissa has been used to generate collocation points. The obtained results have been verified against standard literature to find a satisfactory match. The results include comparative plots for normalized bending and transverse shear stresses, with and without the inclusion of higherorder terms. The results are then extended to study the effect of material index on the deformation and stress behavior of FG beams. The effect of aspect ratio on results is also studied for comparison of various beam theories. 
Since the mid1980s, with the advent of Functionally Graded Materials (FGMs), we have been witnessing a new era in the field of material technology. FGMs belong to a class of advanced materials that have continuous variation in properties along a desired direction and in a desired fashion. Compositions (volume of constituents) and hence the properties gradually change over the volume of such materials, resulting in a corresponding change in the properties of the material that is different from either of the parent materials. Functionally graded materials eliminate the sharp interfaces existing in composite materials and structures, where failure is usually initiated. Due to their customized behavior, FGMs may have very wide applications. If their manufacturing cost is reduced by improving the processes, then it may revolutionize the design process as a whole. A thorough overview of FGMs, their manufacturing techniques, modeling and design, and applications can be found in [1–5]. In [6], an exhaustive review of the modeling and analysis of functionally graded materials and their applications has been reported.
Although considerable research on functionally graded materials has been reported since their conceptualization, most of the work in the area of functionally graded structures (beams and plates) has been done only in the last two decades.
The literature search has been categorized into EulerBernoulli (Classical beam theory, CBT), Timoshenko (First order shear deformation theory, FSDT), and ReddyBickford (Higher order shear deformation theory, HSDT) theories so as to have a clear picture of the developments. Conventionally, CBT is suitable for slender beams with large values of aspect ratio, and as such, the effect of transverse shear deformation can be ignored. In FSDT, a constant state of transverse shear strain with respect to thickness is assumed; however, this theory requires a shear correction factor to accommodate for vanishing shear stresses at the top and bottom fibers of the beam. It has been established that this theory gives satisfactory results for small aspect ratios as well. Secondorder shear deformation theory is rarely adopted, and theories higher than thirdorder shear deformation theory are rarely used because, compared to the efforts required, the accuracy gained is difficult to justify. The thirdorder beam theory, also known as higherorder beam theory, adopts a displacement field such that it develops a parabolic shear strain across any crosssection, requiring no shear correction factor. Later, new theories like Carrera Unified Formulation (CUF) have also emerged to generalize the theories to a wider spectrum. The same concept has been adopted and reported by many researchers in the case of FG beams as well, and a categorical description is put forth in successive paragraphs.
Shankar [7] has studied the behavior of the FG EulerBernoulli beam using the elasticity approach. This work is considered the most fundamental in the realm of FG beams. In this work, an exponential variation of the modulus of elasticity across the thickness is considered, keeping Poisson’s ratio constant, monitoring the variation in stress and displacements, and comparing the behavior of slender and stubby beams under different loading conditions. In [8], Shankar extended his work to include a temperature gradient across the width of the beam while keeping material property variation only dependent on position and independent of temperature, and for the same problem, the author in [9] used the method of Fourier analysis combined with Galerkin’s method for solution. While Shankar has reported the study of the Euler Bernoulli beam, the study has also been extended to include shear deformation for short, stubby beams whose aspect ratio is less than 5. Shanker has reported that the CBT approximates long, slender FG beams within satisfactory limits.
Reddy and Chin [10] and Reddy [11] have studied and presented for the first time the modeling of FG plates on the basis of the shear deformable theory with consideration of the nonlinearities raised by large deflection in addition to the variation in material properties across the thickness. A new beam element has been developed by Chakraborty et al. [12], which can be applied to a bimaterial shear deformable beam having intermediate layers with various functional gradations and is hence useful to determine the mechanical response of functionally graded beams subjected to static and dynamic loads in a thermal environment.
Li [13] has modified the differential equations representing linear and angular deformations of an FG Timoshenko beam by a unique unified method using a common parameter, thus reducing the three displacement variables into a single fourthorder equation. The effects of rotary inertia and shear deformation have been included, and the entire formulation may be analytically reduced to Rayleigh and EulerBernoulli beams after neglecting a few terms.
Static and vibration analysis using higherorder deformation theory for FG beams has been reported in references [14–20]. The approaches to solutions in [14–15, 17–20] are various analytical techniques, while computational techniques were approached in [16]. Li et al. [20] have extended the formulation of [13] to higherorder beam theory for FG beams. In their work, authors have used the elasticity approach to arrive at discrete governing equations for the beam, and then by converting the domain variables into a function of an independent variable, the three equations of motion are simplified to a single equation. The approach of the authors in the paper is unique, and it is necessary to find only one polynomial function for approximating the beam behavior. However, a few terms have been ignored for brevity in the reported paper [20].
The theory of CUF is relatively new in the framework of beam modeling [21]. Carrera et al. have formulated a unified formulation (CUF) using a generic Norder approximation for displacement variables of the crosssection; CBT and FSDT can be considered as particular cases of the CUF. In [22], Giunta et al. applied the theory of CUF to FG beams. In a recent report, Nouri et al. [55] applied the CUF method to the static and buckling analysis of thick FG plates using the Finite Strip method.
In addition to the above, a number of cases pertinent to geometric nonlinearity, elastic foundation, buckling, and dynamic analysis have also been reported in recent literature. However, only a few cases have been considered here for brevity. The effect of geometric nonlinearity in the form of von Karman relations in displacement terms has been accounted for in [23, 24]. Due to the variation in material properties across the crosssection, the neutral axis does not coincide with the center of gravity of the crosssection. The variation of the position of the neutral axis with material gradation has been dealt with [25, 26]. An analysis of the FG beam resting on Winkler’s foundation and the solution using the finite element technique has been reported in [27, 28].
In recent years, the Bspline collocation technique has been increasingly used to find solutions to structural problems as well. Bsplines are widely accepted standards owing to their flexibility and computational efficiency. The collocation technique is quite simple as compared to the finite element method and easy to apply to many problems involving differential equations. A comparative study between the Bspline collocation technique, Finite element, and Finite difference techniques for a twoparameter singularly perturbed boundary value problem has been reported in [29]. The method of Bspline collocation gives better approximation and convergence than the other two techniques. As a thumb rule, we may presume that Bspline collocation may be preferred for comparatively simpler geometries like beams and plates. The technique has been used successfully in fluid flow, onedimensional heat and mass transfer [30], and conduction radiation problems [31]. A novel method of spline collocation for structures was devised and implemented by Bert et al. [32] in 1995 to find approximate solutions for beams and plates to study structural behavior under various loadings. Sun et al. [33] have studied the application of Bspline collocation to problems of linear elasticity using orthogonal cubic Bspline functions. Hsu et al. [34, 35] carried out work on vibration analysis of elastic foundationsupported nonuniform beams and pretwisted beams using spline collocation with uniform knot span in 2009. The effect of nonuniform knot span in the spline collocation procedure applied to a structural problem has been studied by Wu et al. [36]. The author has utilized the concept of Radial Spline functions, which include the knot spacing as a variable, applied the same to beam problems, and showed that nonuniform knot spacing gives a better approximation. The multiplicity of the knots may be an additional parameter in the analysis as it reduces the continuity of the spline at the knots. The effect of multiple knots in a structural problem is considered by Provatidis [37] through finite element analysis using cubic Bspline shape functions. Recently, NURBSbased Isogeometric (IG) analysis are being used by many researchers [38] as a means for refined formulation of Bsplines to formulate shape functions for the determination of approximate solutions to complex problems. Auricchio et al. [39] published a work on the theoretical analysis of a few mathematical problems using the NURBSbased Isogeometric Collocation technique. Reali [40] applied the IG collocation technique for the first time on slender beams and plates and also illustrated the inherent potential of the IG collocation technique in finding the solution of various differential equations. Lieu et al. [41] have applied the IG analysis to studying the bending and free vibration of inplane, bidirectional, variablethickness FG plates. In another report, Lieu and Lee [42] used the IG analysis approach to optimize the material distribution and shape variation for the dynamic analysis of multidirectional functionally graded (MFG) plates. In [43], authors have used the IG approach for reliabilitybased optimization of MFG plates.
In this paper, the Bspline collocation technique using sixthorder Bspline shape functions is used for finding approximate solutions to the governing differential equations for the static analysis of the functionally graded shear deformable beam problem. The authors have extended the work of Li et al. [20] to include a few higherorder terms to obtain a single fourthorder governing equation for FG beams in the framework of higherorder beam theory. While Li et al. [20] have adopted the theory of elasticity approach, in this paper the derivations are made using Hamilton’s principle. A comparative study of the effect of considering higherorder terms in formulation on the behavior of beams is studied and reported in this paper. Subsequently, comparative studies for the deflection and stress behavior of the FG cantilever beams for different cases of aspect ratio and material gradation subjected to uniform loading conditions have been reported.
In this section, a detailed discussion has been presented on the process of forming the governing equation and its solution methodology.
In the present work, a beam has been considered which is made of a functionally graded material. Let us consider, the length of the beam is ‘L’ and it has a rectangular crosssection of width ‘b’ and height ‘h’. It is shown in Figure 1 below. The xaxis is oriented in the axial direction along the midplane of the unbent beam, and the positive zaxis is directed upwards and perpendicular to the xaxis.
Let axial and transverse deflections of the beam be represented by ‘u’ and ‘w, respectively, whereas angular or rotational deflection of the crosssection at the midplane is represented by ‘𝜙’. To account for the variable values at midplane a subscript of value 0 has been used in this work. Further, in the present work, ‘w’ and ‘γ’ have been used to denote or represent the transverse deflection and the shear deformation respectively which are assumed to be functions of ‘x’ and are uniform for a crosssection. In this work, higher order beam theory is considered according to which the top and bottom surfaces of the beam are free from shear traction. To arrive at such boundary conditions, we assume the following conditions for shear deformation ( ) [20]:

(1) 
Fig. 1. Diagram of a Functionally Graded Beam
with coordinate system
For small deformations, we can write [20]:

(2) 
Integrating both sides with respect to ‘z’, the expression for deformation ‘u’ is obtained as

(3) 

(4) 

(5) 
The above Eq. (5) is obtained by substituting, and in Eq. (2). Knowing that we can write:

(6) 

(7) 
where,
, , ,
,
To derive the governing equations of motion for an FG beam, Hamilton’s Principle is applied as shown below.

(8) 

(9) 
Substituting the expressions for U, K, and V in Eq. (8) and assuming,
We get the following governing equation

(10) 
Integrating by parts and categorically separating the variables, the following governing equations and boundary conditions are derived.

(11) 

(12) 

(13) 
Boundary Conditions :

(14) 

(15) 

(16) 
To simplify the above equations, the following substitutions have also been incorporated:
Substituting in expressions for N_{x} and and simplification thereafter, we get

(17) 

(18) 

(19) 

(20) 
The following assumptions have been made for the simplification of the above equations.
and
(where, s =1,2,4)
The expression for N_{x} (Eq. 17) is used to eliminate u_{o} to arrive at the following expressions. Further, it is assumed that the resultant axial force will be zero (for statically determinate cases), so by substituting N_{x}=0 in the above equations obtained after eliminating u_{0}, we arrive at final governing equations for an FG beam in the framework of higherorder beam theory.

(21) 

(22) 
where,
Li et al. [20] have also derived similar equations, but they have neglected the higherorder terms considering E6 and I6. The derivation in their work has been approached using the theory of elasticity approach, and hence the scope of incorporation of these (E6 and I6) terms could not be made. Moreover, the boundary conditions are also affected. In the present work, the derivation of the above governing equation is followed from the first principles using Hamilton’s principle, which also enables us to derive the relevant expressions for boundary conditions. The above formulation can be compared with that of [20], and the above equations are reduced to the same equations by ignoring E6 and I6.
A significant finding of [13, 20] is to reduce the above two equations into a single governing equation by considering an independent parameter ‘F’ such that domain variables i.e., w and are functions of ‘F’ such that it satisfies the governing equations. In the present work, a similar substitution is made after appropriate modification as below and substituted in the governing equation to obtain a single expression for higherorder beam theory, taking into account the higherorder terms.

(23) 
After substituting Eq. (23) to Eq. (21 and 22) the governing equation for FG beams is simplified to:

(24) 
where,
It can be observed that the above equation is exactly similar to the equation derived in reference [41] when the terms E6 and I6 are neglected. The above equation represents the single governing equation for FG beams in the framework of higherorder beam theory. Once the parameter ‘F’ is calculated, the other dependent variables can be easily determined, which can be used to calculate the stresses and strains in the beam for various loading and boundary conditions. We may also calculate the shear forces and bending moments caused by the external loads using the equations discussed in previous paragraphs.
Additionally, by substituting c1 = 0 and c2 = 0 in the above equation, we arrive at Timoshenko beam theory for FG beams [13], which can be further reduced to Rayleigh and Euler beam theory by considering infinite shear modulus. To distinguish the current study and work of Li [20] for higher order beam theory, the terms containing E6 are multiplied with an arbitrary constant whose value is taken as zero or one; when it is zero, it suffices to Li’s work, while its value is one it represents the present study.
Ignoring the time derivatives in the governing equation, we remain with a highly simplified equation for analyzing the behavior of beams under different loading and boundary conditions.

(25) 
The above equation is very similar to the conventional beam equation, and we can use direct methods for obtaining solutions so as to study the beam behavior. However, for future studies pertinent to complex geometric and loading conditions, approximate solutions have also been recorded. Before proceeding with the analysis of the study, we must verify the results for the above formulation.
Here static analysis has been performed on a cantilever beam built of functionally graded material that has material property gradation across the crosssection. The governing differential equation for an FG beam is represented by Eq. (25) while the boundary conditions for a cantilever beam are:

(26) 
The modulus of elasticity follows a power law distribution across the height of the beam, and Poisson’s ratio is assumed to be a constant. The solution to the above problem is approached using the Bspline Collocation technique. The variation of the modulus of elasticity is made according to the power law distribution given by:

(27) 
The subscript ‘β’ refers to the power law index whereas the subscripts ‘b’ and ‘t’ represent the bottom and top fiber respectively. The variation of the modulus of elasticity for various values of the power law index is plotted in Figure 2. The variation of the modulus of rigidity will be similar to the elasticity modulus due to the assumption of constant Poisson’s ratio. The modulus of rigidity is calculated from the relation:

(28) 
Fig. 2. Modulus of Elasticity variation as per power law index, β along the height.
The solution of Eq. (25) is approached using the Bspline collocation technique to calculate parameter ‘F’ which is then used to calculate the transverse deflection ‘w’ and slope ‘𝜙’. Axial and shear stresses are calculated as discussed in the following steps.
To generalize the results, expressions for can be obtained in terms of , as shown below. From Eq. (17) we can write:

(29) 
Integrating both sides, assuming that the resultant axial force vanishes, and using the boundary conditions, we can obtain the equation for u_{0}, which can be further substituted in Eq. (5) to obtain axial deformation, strains, and stresses as below.

(30) 

(31) 

(32) 

(33) 
Bsplines are piecewise polynomials that can be used to approximate a solution to a mathematical problem. They are made up of linear combinations of bspline basis functions. It is illustrated in the parametric space 't'. They depend on a set of nondecreasing coordinates in parametric space known as a knot vector and are defined in terms of the order or degree of the curve. For the present work, we have assumed the following opentype knot vector referred from [47]:
(34) 
in Eq. 34, the knot vector is determined using (n+1) control points and the order of spline ‘k’. A recursive relation has been defined for the
Bspline basis of the ‘k^{th}’ order as below:
(35) 
Equation 35 requires an initiation in the form of a firstorder basis function, which is defined as follows [47]:
(36) 
Using the above basis functions and a series of (n+1) control points (B_{0}, B_{1}... B_{n},), we can define a Bspline curve as:
(37) 
The above Bspline curve (Eq. 37) may be forced to satisfy the governing equation at discrete points called collocation points. There are various methods to determine the collocation points, however, in the current work, the points are calculated using the Greville abscissa [44] approach. In this approach, for a knot vector, T = [t_{0}, t_{1}, t_{2}……, t_{n+k+1}] the points are determined using the relation:
(38) 
where ‘n’ is the degree of the Bspline function and t_{i }are the knots. The number of points obtained may be judiciously selected as per the requirement of the problem.
The Bspline basis functions of sixth order
(k = 6, resulting in a polynomial of degree 5) are selected as per Eq. 37. It is used as a trail solution for function ‘F’ and is substituted in Eq. 25. The resulting equation is then forced to satisfy the governing equation at collocation points given by Greville abscissa.
It can be assessed that the number of collocation points required will be equal to the difference between the number of control points of the Bspline curve and the number of boundary conditions. In the present case, a sixthorder Bspline function with no intermediate knots is used. Hence, we have six coefficients to be determined. As there are 4 boundary conditions available, we require only 2 collocation points to develop six linear equations that can be solved easily to find the coefficients.
In the next section, the solution is approached using the Bspline collocation technique; a detailed discussion of the process is reported by Mahapatra et al. [4648] for reference. A MATLAB code developed for the above formulation is used to obtain computational results.
In the present work, new higherorder unified formulations have been derived for the determination of deflection, axial stress, and shear stress, as mentioned in the previous section. Now the results generated using these formulae have been validated, as discussed below.
To verify the results with standard literature references, we consider the data as discussed in a numerical example of a cantilever beam in Li et al. [20]. The considered beam is of length L = 0.5m, depth h = L/3, and width of unit length. Aluminum has been considered at the bottom and silicon carbide has been considered at the top as the constituent materials for the functionally graded material for this beam, whose properties have been mentioned below.
A constant value of Poisson’s ratio is assumed here as 0.3. The value of the power law index is assumed to be equal to 1. It has been assumed that at the top of the beam, a pressure load of 1kN/m is being imposed uniformly. The axial and shear stresses are normalized as follows:
Fig. 3. Verification of present work, a) normalized axial stress
b) normalized shear stresses (ignoring E6)
The normalized axial and shear stresses are calculated at different locations in the span of the beam and are presented in Figures 3(ab).
As can be observed from Figure 3, the results of this study match exactly those of Li et al. [20] while ignoring the higherorder terms that contain E_{6}. The significance of the terms containing E_{6} on the behavior of FG beams is explored in successive paragraphs.
A comparison of normalized axial and shear stresses for the higherorder beam theory used in [20] with the higherorder beam theory of the present work (considering E_{6}) reveals that the effect of E_{6} is negligible for axial stresses; however, from Figure 4, it can be clearly observed that the effect of consideration of the higher order terms has a moderate effect (approximately 4.083%) on shear stresses and a very small difference of 0.836% in the normalized axial stresses.
Fig. 4. Comparison plots for the effect of higher order terms on shear stress (b=1)
In this section, the behavior of FG beams for different cases of material distribution (using different index values) and the effect of shear modulus on the stresses are discussed. In Figures 5a & b, the normalized axial and shear stresses at the fixed end of the beam are respectively plotted for different values of the power law index ranging from (0, 0.2, 1, 2, and 5). Hence, a comparative study of axial and shear stresses in FG beams with an equivalent isotropic beam can be observed. The asymmetry in the material composition results in behavioral variations, as reported in Figure 5. With the increase in index value, the maximum axial stress (at the top and bottom fibers) varies considerably as compared to the isotropic case. The maximum shear stress also increases with increasing index up to a certain maximum, then decreases.
Fig. 5. Effect of material gradation on stresses (fixed end),
a) Normalized axial stresses b) Normalized shear stresses
The study of the effect of material gradient on the deflection and stresses of an FG beam can be further extended using threedimensional plots that give better visualization and understanding. Figure 6(ac) reports the 3D plots corresponding to the beam and can be referred to for the study of the effect of material index on the beam behavior.
(a) 
(b) 
(c) 
Fig. 6. 3D plots to visualize the effect of material gradation on beam behavior a) Normalized axial stress,
b) Normalized shear stress, c) Deflection
In Figure 7, comparative plots for different beam theories are reported for axial and shear stresses for different values of aspect ratio. For Euler’s theory, the shear modulus is considered to be infinite while for Timoshenko's theory, the value of c_{1} is assumed to be zero. To distinguish the current study and Li [20] for higher order beam theory, the terms containing E_{6} are multiplied with an arbitrary constant whose value is taken as zero/ one; when it is zero, it suffices to Li’s work while its value is one for the present study. The normalized transverse deflection is calculated using the relation:
In Figure7(a&b), the deflection of the cantilever corresponding to various beam theories is plotted for two aspect ratios of 3 and 10. It can be observed that for L/h=3, there is a significant difference in deflection corresponding to Euler’s theory, while, Li’s work and present deflections are almost equal, differences being negligible. Timoshenko's deflection only marginally varies from the higherorder cases. However, for L/h=10, all types of deflections are practically indistinguishable. The results are in line with the previous similar reported studies.
It can be observed from Figure7(c&d) that the axial stresses calculated by all the beam theories are nearly equal for both aspect ratios. The normalized axial stress values are proportionally increasing for increasing aspect ratios while the variation is of the same profile. However, shear stresses are significantly different for the various theories considered in Figure 6 (e&f). The shear stresses are zero for Euler’s theory due to the infinite value of shear modulus while for other theories, the observed values of shear stresses are parabolically obtained with zero shear at top and bottom surfaces. The values of shear stresses obtained using higher order beam theory [20] are slightly higher than that of Timoshenko theory and the results obtained using the present work are in between the stresses obtained by the Timoshenko beam [13] and higher order beam theory reported in [20].
The average percentage difference between the results, with and without the inclusion of higher order terms (i.e., present work and cited literature) for normalized axial stress is 0.8355% and for normalized shear stress, it is 4.083% (For details, please refer to appendix).



(a) 



(A) 



(b) 


(B) 



(c) 



(C) 

Fig. 7. Comparison of various results according to various beam theories for L/h ratio 3 &10 – (a, A) Beam Deflection for L/h=3 &10, (b, B) Normalized Axial Stress for L/h=3&10, and (c, C) Normalized Shear Stress for L/h=3&10
To continue with the numeric experiments, normalized transverse deflection is studied as a function of the material gradient index. The maximum deflection that is obtained is a function of index value as reported in the literature. Hence to observe such effects in the present case along with the effect of aspect ratio using various beam theories is reported in Figure 8. For an indepth observation, the study has been categorized into two cases;
CaseA: when ceramic is on the top surface, CaseB: the metal surface is on top with the external uniformly distributed load acting on the top surface in the downward direction. It can be expected that changing the loading surface will also have a substantial effect on the behavior of the beam for a given material composition.
To continue with the numeric experiments, normalized transverse deflection is studied as a function of the material gradient index. The maximum deflection that is obtained is a function of index value as reported in the literature. Hence to observe such effects in the present case along with the effect of aspect ratio using various beam theories is reported in Figure 8. For an indepth observation, the study has been categorized into two cases; CaseA: when ceramic is on the top surface, CaseB: the metal surface is on top with the external uniformly distributed load acting on the top surface in the downward direction. It can be expected that changing the loading surface will also have a substantial effect on the behavior of the beam for a given material composition.
It can be clearly observed from Figure8a that for thick stubby beams, there is a clear difference in the maximum transverse displacement corresponding to Euler, Timoshenko, and higherorder theories. This difference is small for lower gradients, and increases as the value of index β increases for L/h=3. This difference in normalized transverse displacement reduces as the aspect ratio increases; for L/h =10, the difference is practically insignificant as in Figure8b. Another noticeable observation in Figure 8 (a&b) is that the slope of the curve (rate of increase in maximum deflection) decreases as the material gradient increases. This can be explained by the fact that the percentage of metal (less stiff component) increases as the value of β increases, (when β=0 then the material at the top surface is 100%, for β=1, both components are 50% and further increase in β increases the composition of metal). It can be observed that the effect of the material gradient is more significant for values up to β=2; the effect is less significant for higher index values.
In the next study, i.e. CaseB, the beam is inverted to expose the metal part to the external load applied at the top surface. For β=0, i.e. for isotropic material made of metal, the deflection is maximum and as the stiffer component fraction increases with the value of β the deflection decreases considerably. This can be observed in Figure7c&d for aspect ratios 3 and 10. Similar to the previous case the difference in the Euler and higher order displacements is significant for lower aspect ratios and is practically insignificant for aspect ratio greater than 10.




Fig. 8. Dimensionless maximum transverse
deflection as a function of material gradient index for
Case A: Metal at base and Ceramic at top and
Case B: Metal at top and Ceramic at base.
Table 1. Material Properties
Sl. 
Material 
Modulus of 
Poisson’s 
Ref. 
1. 
Aluminum 
70GPa 
0.30 
[41] 
2. 
Silicon Nitride 
348GPa 
0.24 
[41] 
3. 
Silicon Carbide 
440GPa 
0.27 
[40] 
4. 
Steel 
210GPa 
0.30 
[11] 
In this section, three discrete combinations of materials have been selected to observe the variation in beam behavior due to an arbitrary ratio of material properties. The material properties of selected materials are mentioned in Table 1. The combination of materials with the ratio of material properties is mentioned in Table 2.
Table 2. Combination of materials
Set 
Materials 
Ratio of Modulus of Elasticity 
SetA 
Aluminum+ Silicon Nitride 
(348GPa /70GPa) ⁓5 
SetB 
Aluminum+ Silicon Carbide 
(440GPa /70GPa) ⁓7 
SetC 
Aluminum+ Steel 
(210GPa /70GPa) ⁓3 
The three combinations can be observed to have a discrete ratio of modulus of elasticity
(3, 5, and 7). The behavior of such beams corresponding to β=1 has been plotted in Figure9 to observe the variations and relate the same to the ratio of material properties. The behavior of such beams corresponding to β=1 under uniformly distributed load for a comparative study is mentioned in Figure 9.
The deflection of the cantilever beam is plotted in Figure 9a. From the figure, it can be observed that the composition corresponding to SetB is stiffer and that validates the elasticity moduli ratio. However, when we observe the stresses, the variation is as shown in Figure 9b for axial stress and 9c for shear stress. The stresses are found to be higher for higher moduli ratios. In the present case, the normalized tensile stresses at the top surface are higher for SetB for moduli ratio 7 while the compressive stresses are on the lower side for the same set of materials. On the contrary, for SetB with moduli ratio 3, the maximum tensile stress at the top fiber is lower, and compressive stress at the lower surface is higher than the other sets. Another observation from Figure 9b is that the steepness (curvature) of the axial stress curves increases with an increase in the moduli ratio. Figure 9c shows the normalized shear stress plots for the three combinations of materials. It is clear that shear stress is higher for SetB with a higher moduli ratio, which can also be inferred from the axial stress values. This discussion is pertinent to the design plan of FG beams.
Fig. 9. Comparative plots for various combinations of materials with aluminum a) Normalized Transverse Deflection b) Normalized Axial Stress
c) Normalized Shear Stress
The study is further extended to include simply supported boundary conditions (bending moment and deflection are zero at the supports) for the FG beam. In Figure 9, deflection and stress plots corresponding to simply supported boundary conditions are presented. It is observed that the differences in axial stress and shear stress are similar in trend as for cantilever support. The percentage difference in axial stresses for a cantilever beam is about 1.2% for an isotropic beam and 1.9% for an FG beam. On the other hand, the difference in shear stress is 7.7% for the isotropic beam and 8.0% for the FG beam (details can be referred to in the appendix)



Fig. 10. Comparison of the isotropic and FG beam behavior with/without higher order terms a) deflection,
b) normalized axial stress, c) normalized shear stress
Using Hamilton’s principle, a higherorder beam formulation is derived to obtain a single governing equation for functionally graded beams by incorporating higherorder terms that were not considered in previous similar work. The unification of three developed equations into a single equation is done by considering a parameter that is a function of domain variables.
The solution to the governing equation is approached using the approximation technique of Bspline collocation. A cantilever functionally graded beam with grading along the crosssection is studied for its behavioral properties under the action of a distributed mechanical load. The variation of modulus of elasticity and modulus of rigidity is considered by power law, and Poisson’s ratio is considered to be a constant. The collocation points are calculated using Greville abscissa. It is observed that the new higherorder terms have a moderate effect (approximately 4%) on the shear stress calculations compared to the similar report cited in the literature and only 0.8% on normalized axial stresses. The normal stresses are negligibly affected by the new terms considered for cantilever beams. The difference in values corresponding to simply supported end conditions is slightly on the higher side; 1.8% for axial stress and 8% for shear stresses. The study is then extended to explore the effect of material gradient, aspect ratio, and moduli ratio on transverse deformation, normalized axial stresses, and normalized shear stress variations.
As a future scope, the present model may be compared with relatively newer techniques like Carrera Unified Formulation and the HighOrder Continuation Model for the study of FG structural elements. Thermal stresses and dynamic analysis are other areas in which the present study can be extended.
Nomenclature
X 
Variables should appear in the first column with the definition in the second column 
t & Me 
One and twoletter abbreviations should appear in italics 
Sim 
Threeletter abbreviations should not appear in italics 
Fr 
Dimensionless numbers and parameters don't appear in italics 
Conflicts of Interest
The author declares that there is no conflict of interest regarding the publication of this manuscript. In addition, the authors have entirely observed the ethical issues, including plagiarism, informed consent, misconduct, data fabrication and/or falsification, double publication and/or submission, and redundancy.
Appendixes
Table A1. Normalized Axial Stresses at fixed end of beam
z/h 
b=0 
b =0.2 
b =0.5 
b =1 
b =2 
b =5 

Present 
[20] 
Present 
[20] 
Present 
[20] 
Present 
[20] 
Present 
[20] 
Present 
[20] 

0.5 
8.9257 
9 
2.3328 
2.3506 
3.1346 
3.156 
4.2939 
4.3195 
5.6848 
5.7161 
6.6325 
6.6777 
0.4 
7.1905 
7.2 
6.6903 
6.7002 
5.8678 
5.8761 
5.0416 
5.047 
5.0153 
5.018 
5.6359 
5.6382 
0.3 
5.422 
5.4 
5.7254 
5.708 
5.7149 
5.701 
5.2249 
5.2128 
4.5748 
4.5619 
4.6219 
4.6041 
0.2 
3.6285 
3.6 
4.2798 
4.2546 
4.7926 
4.7706 
4.8367 
4.8169 
4.1543 
4.1358 
3.614 
3.5928 
0.1 
1.8184 
1.8 
2.5698 
2.553 
3.3559 
3.3406 
3.8741 
3.8594 
3.5418 
3.5276 
2.6344 
2.6211 
0 
0 
0 
0.6827 
0.6837 
1.5219 
1.5236 
2.339 
2.3403 
2.5253 
2.5254 
1.6611 
1.6625 
0.1 
1.8184 
1.8 
1.3317 
1.3117 
0.6374 
0.616 
0.2373 
0.2594 
0.8958 
0.917 
0.5549 
0.5749 
0.2 
3.6285 
3.6 
3.4392 
3.4084 
3.0696 
3.0358 
2.4204 
2.3831 
1.5498 
1.5094 
1.0406 
1.0014 
0.3 
5.422 
5.4 
5.6133 
5.5895 
5.7328 
5.706 
5.6191 
5.5872 
5.0056 
4.966 
3.8218 
3.7749 
0.4 
7.1905 
7.2 
7.8316 
7.8432 
8.5906 
8.6039 
9.3394 
9.353 
9.6539 
9.6647 
8.9657 
8.9649 
0.5 
8.9257 
9 
10.0746 
10.1607 
11.6101 
11.7118 
13.5576 
13.6805 
15.662 
15.8176 
18.2803 
18.4926 
Table A2. Normalized Shear Stresses at fixed end beam
z/h 
b =0 
b =0.2 
b =0.5 
b =1 
b =2 
b =5 

Present 
[20] 
Present 
[20] 
Present 
[20] 
Present 
[20] 
Present 
[20] 
Present 
[20] 

0.5 
0 
0 
0 
0 
0 
0 
0 
0 
0 
0 
0 
0 
0.4 
0.5143 
0.54 
0.4119 
0.4318 
0.3121 
0.3272 
0.2403 
0.2527 
0.2424 
0.2562 
0.3573 
0.3788 
0.3 
0.9143 
0.96 
0.8101 
0.8491 
0.6828 
0.7158 
0.5487 
0.5769 
0.4803 
0.5077 
0.636 
0.6743 
0.2 
1.2 
1.26 
1.1299 
1.1844 
1.025 
1.0745 
0.8795 
0.9248 
0.7384 
0.7805 
0.8417 
0.8924 
0.1 
1.3714 
1.44 
1.3492 
1.4143 
1.2956 
1.3582 
1.1874 
1.2485 
1.0167 
1.0746 
0.9915 
1.0512 
0 
1.4286 
1.5 
1.4547 
1.5248 
1.4635 
1.5342 
1.4266 
1.5 
1.2906 
1.3641 
1.1156 
1.1828 
0.1 
1.3714 
1.44 
1.4367 
1.506 
1.5038 
1.5765 
1.5517 
1.6315 
1.5105 
1.5966 
1.247 
1.3221 
0.2 
1.2 
1.26 
1.2879 
1.35 
1.3954 
1.4629 
1.5171 
1.5952 
1.6025 
1.6939 
1.3901 
1.4739 
0.3 
0.9143 
0.96 
1.0022 
1.0505 
1.1196 
1.1737 
1.2774 
1.3431 
1.4679 
1.5515 
1.4618 
1.5498 
0.4 
0.5143 
0.54 
0.5744 
0.6021 
0.6596 
0.6915 
0.7868 
0.8273 
0.9831 
1.0391 
1.1952 
1.2672 
0.5 
0 
0 
0 
0 
0 
0 
0 
0 
0 
0 
0 
0 
Table B1. Normalized Shear Stresses at mid span of simply supported beam
z/h 
b =0 
%diff 
b =1 
%diff 

Li et al. 
Present 
Li et al. 
Present 

0.5 
2.2933 
2.3243 
0.0135176 
1.0799 
1.1055 
0.023706 
0.4 
1.7847 
1.8095 
0.0138959 
1.2617 
1.2671 
0.00428 
0.3 
1.3094 
1.328 
0.014205 
1.3032 
1.2911 
0.009285 
0.2 
0.8591 
0.8715 
0.0144337 
1.2042 
1.1845 
0.016359 
0.1 
0.4254 
0.4316 
0.0145745 
0.9649 
0.9502 
0.015235 
0 
0 
0 
0 
0.5851 
0.5864 
0.002222 
0.1 
0.4254 
0.4316 
0.0145745 
0.0649 
0.087 
0.340524 
0.2 
0.8591 
0.8715 
0.0144337 
0.5958 
0.5585 
0.062605 
0.3 
1.3094 
1.328 
0.014205 
1.3968 
1.365 
0.022766 
0.4 
1.7847 
1.8095 
0.0138959 
2.3383 
2.3519 
0.005816 
0.5 
2.2933 
2.3243 
0.0135176 
3.4201 
3.5431 
0.035964 
Table B2. Normalized Shear Stresses at mid span of simply supported beam
z/h 
b=0 
%diff 
b =1 
%diff 

Li et al. 
Present 
Li et al. 
Present 

0.5 
0 
0 
0 
0 
0 
0 
0.4 
0.27 
0.2443 
0.0951852 
0.1263 
0.114 
0.097387 
0.3 
0.48 
0.4343 
0.0952083 
0.2885 
0.2602 
0.098094 
0.2 
0.63 
0.57 
0.0952381 
0.4624 
0.4171 
0.097967 
0.1 
0.72 
0.6514 
0.0952778 
0.6242 
0.5631 
0.097885 
0 
0.75 
0.6786 
0.0952 
0.75 
0.6766 
0.097867 
0.1 
0.72 
0.6514 
0.0952778 
0.8158 
0.7359 
0.097941 
0.2 
0.63 
0.57 
0.0952381 
0.7976 
0.7195 
0.097919 
0.3 
0.48 
0.4343 
0.0952083 
0.6715 
0.6058 
0.097841 
0.4 
0.27 
0.2443 
0.0951852 
0.4137 
0.3732 
0.097897 
0.5 
0 
0 
0 
0 
0 
0 
References