Higher Order Approximations for Bending of FG Beams
Using B-Spline Collocation Technique

1. Biswas ^a , D. Mahapatra ^b* , S. C. Mondal ^c, S. Sarkar ^c

^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

1.     Introduction

Since the mid-1980s, 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 Euler-Bernoulli (Classical beam theory, CBT), Timoshenko (First order shear deformation theory, FSDT), and Reddy-Bickford (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. Second-order shear deformation theory is rarely adopted, and theories higher than third-order shear deformation theory are rarely used because, compared to the efforts required, the accuracy gained is difficult to justify. The third-order beam theory, also known as higher-order beam theory, adopts a displacement field such that it develops a parabolic shear strain across any cross-section, 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 Euler-Bernoulli 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 bi-material 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 fourth-order equation. The effects of rotary inertia and shear deformation have been included, and the entire formulation may be analytically reduced to Rayleigh and Euler-Bernoulli beams after neglecting a few terms.

Static and vibration analysis using higher-order 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 higher-order 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 N-order approximation for displacement variables of the cross-section; 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 non-linearity, 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 cross-section, the neutral axis does not coincide with the center of gravity of the cross-section. 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 B-spline collocation technique has been increasingly used to find solutions to structural problems as well. B-splines 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 B-spline collocation technique, Finite element, and Finite difference techniques for a two-parameter singularly perturbed boundary value problem has been reported in [29]. The method of B-spline collocation gives better approximation and convergence than the other two techniques. As a thumb rule, we may presume that B-spline collocation may be preferred for comparatively simpler geometries like beams and plates. The technique has been used successfully in fluid flow, one-dimensional 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 B-spline collocation to problems of linear elasticity using orthogonal cubic B-spline functions. Hsu et al. [34, 35] carried out work on vibration analysis of elastic foundation-supported non-uniform beams and pre-twisted beams using spline collocation with uniform knot span in 2009. The effect of non-uniform 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 non-uniform 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 B-spline shape functions. Recently, NURBS-based Isogeometric (IG) analysis are being used by many researchers [38] as a means for refined formulation of B-splines 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 NURBS-based 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 in-plane, bi-directional, variable-thickness 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 reliability-based optimization of MFG plates.

In this paper, the B-spline collocation technique using sixth-order B-spline 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 higher-order terms to obtain a single fourth-order governing equation for FG beams in the framework of higher-order 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 higher-order 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.

2.     Formulation

In this section, a detailed discussion has been presented on the process of forming the governing equation and its solution methodology.

2.1.  Derivation of Governing Equations

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 cross-section of width ‘b’ and height ‘h’. It is shown in Figure 1 below. The x-axis is oriented in the axial direction along the mid-plane of the unbent beam, and the positive z-axis is directed upwards and perpendicular to the x-axis.

Let axial and transverse deflections of the beam be represented by ‘u’ and ‘w, respectively, whereas angular or rotational deflection of the cross-section at the mid-plane is represented by ‘𝜙’. To account for the variable values at mid-plane 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 cross-section. 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]:

Fig. 1. Diagram of a Functionally Graded Beam
 with co-ordinate system

For small deformations, we can write [20]:

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

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

where,

, , ,

,

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

Substituting the expressions for U, K, and V in Eq. (8) and assuming,

We get the following governing equation

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

Boundary Conditions :

To simplify the above equations, the following substitutions have also been incorporated:

Substituting  in expressions for N_x and  and simplification thereafter, we get

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₀, we arrive at final governing equations for an FG beam in the framework of higher-order beam theory.

where,

Li et al. [20] have also derived similar equations, but they have neglected the higher-order 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 higher-order beam theory, taking into account the higher-order terms.

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

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 higher-order 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.

2.2.  Stress Analysis of Functionality Graded Beam using Higher Order Beam Theory

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.

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 cross-section. The governing differential equation for an FG beam is represented by Eq. (25) while the boundary conditions for a cantilever beam are:

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 B-spline Collocation technique. The variation of the modulus of elasticity is made according to the power law distribution given by:

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:

Fig. 2. Modulus of Elasticity variation as per power law index, β along the height.

The solution of Eq. (25) is approached using the B-spline 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:

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

2.3.  B-Spline Basis Functions

B-splines are piece-wise polynomials that can be used to approximate a solution to a mathematical problem. They are made up of linear combinations of b-spline basis functions. It is illustrated in the parametric space 't'. They depend on a set of non-decreasing 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 open-type knot vector referred from [47]:

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
B-spline basis of the ‘k^th’ order as below:

Equation 35 requires an initiation in the form of a first-order basis function, which is defined as follows [47]:

Using the above basis functions and a series of (n+1) control points (B₀, B₁... B_n,), we can define a  B-spline curve as:

The above B-spline 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₀, t₁, t₂……, t_n+k+1] the points are determined using the relation:

where ‘n’ is the degree of the B-spline function and t_i are the knots. The number of points obtained may be judiciously selected as per the requirement of the problem.

2.4.  B-Spline Collocation Technique

The B-spline 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 B-spline curve and the number of boundary conditions. In the present case, a sixth-order B-spline 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 B-spline collocation technique; a detailed discussion of the process is reported by Mahapatra et al. [46-48] for reference. A MATLAB code developed for the above formulation is used to obtain computational results.

3.     Results and Discussions

In the present work, new higher-order 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.

3.1.  Verification of Results

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(a-b).

As can be observed from Figure 3, the results of this study match exactly those of Li et al. [20] while ignoring the higher-order terms that contain E₆. The significance of the terms containing E₆ on the behavior of FG beams is explored in successive paragraphs.

A comparison of normalized axial and shear stresses for the higher-order beam theory used in [20] with the higher-order beam theory of the present work (considering E₆) reveals that the effect of E₆ 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)

3.2.  Other Numerical Experiments

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 three-dimensional plots that give better visualization and understanding. Figure 6(a-c) 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.

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₁ is assumed to be zero. To distinguish the current study and Li [20] for higher order beam theory, the terms containing E₆ 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 Figure-7(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 higher-order 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 Figure-7(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).

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 in-depth observation, the study has been categorized into two cases;

Case-A: when ceramic is on the top surface, Case-B: 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 in-depth observation, the study has been categorized into two cases; Case-A: when ceramic is on the top surface, Case-B: 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 Figure-8a that for thick stubby beams, there is a clear difference in the maximum transverse displacement corresponding to Euler, Timoshenko, and higher-order 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 Figure-8b. 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. Case-B, 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 Figure-7c&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

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

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 Figure-9 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 Set-B 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 Set-B for moduli ratio 7 while the compressive stresses are on the lower side for the same set of materials. On the contrary, for Set-B 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 Set-B 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

4.     Conclusions

Using Hamilton’s principle, a higher-order beam formulation is derived to obtain a single governing equation for functionally graded beams by incorporating higher-order 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 B-spline collocation. A cantilever functionally graded beam with grading along the cross-section 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 higher-order 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 High-Order 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

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

Table A2. Normalized Shear Stresses at fixed end beam

Table B1. Normalized Shear Stresses at mid span of simply supported beam

Table B2. Normalized Shear Stresses at mid span of simply supported beam

References

1. Miyamoto Y., Kaysser W. Rabin B. H., Kawasaki A., 1999. Functionally Graded Materials- Design Processing and Applications. Springer.
2. Suresh, S. and Mortensen, A., 1998. Fundamentals of Functionally Graded Materials. The Institute of Materials, IOM Communications Ltd., London.
3. Mahamood, R.M., Akinlabi, E.T., Shukla, M. and Pityana, S.L., 2012. Functionally graded material: an overview.
4. Udupa, G., Rao, S.S. and Gangadharan, K.V., 2014. Functionally graded composite materials: an overview. Procedia Materials Science, 5, 1291-1299.
5. Kieback, B., Neubrand, A. and Riedel, H., 2003. Processing techniques for functionally graded materials. Materials Science and Engineering: A, 362(1-2), 81-106.
6. Katunin, A., Wachla, D., Santos, P. and Reis, P.N., 2023. Fatigue life assessment of hybrid bio-composites based on self-heating temperature. Composite Structures, 304, p.116456.
7. Sankar, B.V., 2001. An elasticity solution for functionally graded beams. Composites Science and Technology, 61(5), 689-696.
8. Sankar, B.V. and Tzeng, J.T., 2002. Thermal stresses in functionally graded beams. AIAA journal, 40(6), 1228-1232.
9. Zhu, H. and Sankar, B.V., 2004. A combined Fourier series–Galerkin method for the analysis of functionally graded beams.  Appl. Mech., 71(3), pp.421-424.
10. Reddy, J.N. and Chin, C.D., 1998. Thermomechanical analysis of functionally graded cylinders and plates. Journal of thermal Stresses, 21(6), 593-626.
11. Reddy, J.N. and Chin, C.D., 1998. Thermomechanical analysis of functionally graded cylinders and plates. Journal of thermal Stresses, 21(6), 593-626.
12. Chakraborty, A., Gopalakrishnan, S. and Reddy, J., 2003. A new beam finite element for the analysis of functionally graded materials. International journal of mechanical sciences, 45(3), 519-539.
13. Li, X.F., 2008. A unified approach for analyzing static and dynamic behaviors of functionally graded Timoshenko and Euler–Bernoulli beams. Journal of Sound and vibration, 318(4-5), 1210-1229.
14. Aydogdu, M. and Taskin, V., 2007. Free vibration analysis of functionally graded beams with simply supported edges. Materials & design, 28(5), pp.1651-1656.
15. Zhong, Z. and Yu, T., 2007. Analytical solution of a cantilever functionally graded beam. Composites Science and Technology, 67(3-4), pp.481-488.
16. Kadoli, R., Akhtar, K. and Ganesan, N., 2008. Static analysis of functionally graded beams using higher order shear deformation theory. Applied mathematical modelling, 32(12), pp.2509-2525.
17. Benatta, M.A., Mechab, I., Tounsi, A. and Bedia, E.A., 2008. Static analysis of functionally graded short beams including warping and shear deformation effects. Computational Materials Science, 44(2), 765-773.
18. Sina, S.A., Navazi, H.M. and Haddadpour, H., 2009. An analytical method for free vibration analysis of functionally graded beams. Materials & Design, 30(3), pp.741-747.
19. Şimşek, M., 2010. Fundamental frequency analysis of functionally graded beams by using different higher-order beam theories. Nuclear Engineering and Design, 240(4), pp.697-705.
20. Li, X.F., Wang, B.L. and Han, J.C., 2010. A higher-order theory for static and dynamic analyses of functionally graded beams. Archive of Applied Mechanics, 80, pp.1197-1212.
21. Carrera, E., Giunta, G. and Petrolo, M., 2011. Beam structures: classical and advanced theories. John Wiley & Sons.
22. Giunta, G., Belouettar, S. and Carrera, E., 2010. Analysis of FGM beams by means of classical and advanced theories. Mechanics of Advanced Materials and Structures, 17(8), pp.622-635.
23. Tahani, M., Torabizadeh, M.A. and Fereidoon, A., 2006. Nonlinear analysis of functionally graded beams. Journal of Achievements in Materials and Manufacturing Engineering, 18(1-2),315-318.
24. Reddy, J., 2011. Microstructure-dependent couple stress theories of functionally graded beams. Journal of the Mechanics and Physics of Solids, 59(11),2382-2399.
25. Zhang, D.G., 2013. Nonlinear bending analysis of FGM beams based on physical neutral surface and high order shear deformation theory. Composite Structures, 100, pp.121-126.
26. Yaghoobi, H. and Fereidoon, A., 2010. Influence of neutral surface position on deflection of functionally graded beam under uniformly distributed load. World Applied Sciences Journal, 10(3), pp.337-341.
27. Mohanty, S.C., Dash, R.R. and Rout, T., 2011. Parametric instability of a functionally graded Timoshenko beam on Winkler's elastic foundation. Nuclear Engineering and Design, 241(8), pp.2698-2715.
28. Mohanty, S.C., Dash, R.R. and Rout, T., 2012. Static and dynamic stability analysis of a functionally graded Timoshenko beam. International Journal of Structural Stability and Dynamics, 12(04), p.1250025.
29. Kadalbajoo, M.K. and Yadaw, A.S., 2011. Finite difference, finite element and b-spline collocation methods applied to two parameter singularly perturbed boundary value problems1. Jnaiam, 5, pp.163-180.
30. Chawla, T.C., Leaf, G. and Chen, W., 1975. A collocation method using B-splines for one-dimensional heat or mass-transfer-controlled moving boundary problems. Nuclear Engineering and Design, 35(2), pp.163-180.
31. Chawla, T.C. and Chan, S.H., 1979. Solution of radiation-conduction problems with collocation method using B-splines as approximating functions. International Journal of Heat and Mass Transfer, 22(12),1657-1667.
32. Bert, C.W. and Sheu, Y., 1996. Static analyses of beams and plates by spline collocation method. Journal of engineering mechanics, 122(4), pp.375-378.
33. Sun, W., 2001. B-spline collocation methods for elasticity problems. Scientific Computing and Applications, 2001, pp.133-141.
34. Hsu, M.H., 2009. Vibration analysis of non-uniform beams resting on elastic foundations using the spline collocation method. Journal of Applied Science and Engineering, 12(2), pp.113-122.
35. Hsu, M.H., 2009. Vibration analysis of pre-twisted beams using the spline collocation method. Journal of Marine Science and Technology, 17(2),4.
36. Wu, L.Y., Chung, L.L. and Huang, H.H., 2008. Radial spline collocation method for static analysis of beams. Applied mathematics and computation, 201(1-2), pp.184-199.
37. Provatidis, C., 2014. Finite element analysis of structures using-continuous cubic B-splines or equivalent Hermite elements. Journal of Structures, 2014.
38. Cottrell, J.A., Hughes, T.J. and Bazilevs, Y., 2009. Isogeometric analysis: toward integration of CAD and FEA. John Wiley & Sons.
39. Auricchio, F., Da Veiga, L.B., Hughes, T.J.R., Reali, A. and Sangalli, G., 2010. Isogeometric collocation methods. Mathematical Models and Methods in Applied Sciences, 20(11),2075-2107.
40. Reali, A. and Gomez, H., 2015. An isogeometric collocation approach for Bernoulli–Euler beams and Kirchhoff plates. Computer Methods in Applied Mechanics and Engineering, 284, pp.623-636.
41. Lieu, Q.X., Lee, S., Kang, J. and Lee, J., 2018. Bending and free vibration analyses of in-plane bi-directional functionally graded plates with variable thickness using isogeometric analysis. Composite Structures, 192, pp.434-451.
42. Lieu, Q.X. and Lee, J., 2019. An isogeometric multimesh design approach for size and shape optimization of multidirectional functionally graded plates. Computer Methods in Applied Mechanics and Engineering, 343,407-437.
43. Lieu, Q.X. and Lee, J., 2019. A reliability-based optimization approach for material and thickness composition of multidirectional functionally graded plates. Composites Part B: Engineering, 164, pp.599-611.
44. Johnson, R.W., 2005. Higher order B-spline collocation at the Greville abscissae. Applied Numerical Mathematics, 52(1), pp.63-75.
45. Hayun, S., Paris, V., Mitrani, R., Kalabukhov, S., Dariel, M.P., Zaretsky, E. and Frage, N., 2012. Microstructure and mechanical properties of silicon carbide processed by Spark Plasma Sintering (SPS). Ceramics International, 38(8), pp.6335-6340.
46. Mahapatra, D., Sanyal, S. and Bhowmick, S., 2019. An approximate solution of functionally Graded Timoshenko beam using B-spline collocation method. Journal of Solid Mechanics, 11(2), pp.297-310.
47. Mahapatra D., Sanyal S., Bhowmick, S., 2019. B-spline collocation-based approximations for structural problems of fourth order. UPB Scientific Bulletin, Series-D; 81(4), pp.213-228
48. Mahapatra, D., Sanyal, S. and Bhowmick, S., 2020. An Investigation of Stress and Deformation Behavior of Functionally Graded Timoshenko Beams subjected to Thermo-Mechanical Load. Mechanics Of Advanced Composite Structures, 7(1),157-176.
49. Chaabani, H., Mesmoudi, S., Boutahar, L. and Bikri, K.E., 2022. Buckling of porous FG sandwich plates subjected to various non-uniform compressions and resting on Winkler–Pasternak elastic foundation using a finite element model based on the high-order shear deformation theory. Acta Mechanica, 233(12), pp.5359-5376.
50. Chaabani, H., Mesmoudi, S., Boutahar, L. and El Bikri, K., 2022. A high-order continuation for bifurcation analysis of functionally graded material sandwich plates. Acta Mechanica, 233(6),2125-2147.
51. Chaabani, H., Mesmoudi, S., Boutahar, L. and El Bikri, K., 2023. A high-order finite element continuation for buckling analysis of porous FGM plates. Engineering Structures, 279, p.115597.
52. Mesmoudi, S., Askour, O., Rammane, M., Bourihane, O., Tri, A. and Braikat, B., 2022. Spectral Chebyshev method coupled with a high order continuation for nonlinear bending and buckling analysis of functionally graded sandwich beams. International Journal for Numerical Methods in Engineering, 123(24), pp.6111-6126.
53. Mesmoudi, S., Hilali, Y., Rammane, M., Askour, O. and Bourihane, O., 2023. Highly efficient mesh-free approach to simulate the non-linear bending analysis of fg porous beams and sandwich beams with fg face sheets. Thin-Walled Structures, 185, p.110614.
54. Zahari, K., Hilali, Y., Mesmoudi, S. and Bourihane, O., 2022, December. Review and comparison of thin and thick FGM plate theories using a unified buckling formulation. In Structures (Vol. 46), 1545-1560.
55. Nouri, Z., Sarrami-Foroushani, S., Azhari, F. and Azhari, M., 2022. Application of Carrera unified formulation in conjunction with finite strip method in static and stability analysis of functionally graded plates. Mechanics of Advanced Materials and Structures, 29(2), pp.250-266.