Approximate Solutions for Bending of Beams and Buckling of Columns Made of Functionally Graded Materials

Abstract

In this research, approximate solutions are developed to solve functionally graded material (FGM) used in beams/columns. FGM is considered in the transverse direction of beams/columns using power/exponential/sigmoidal functions considering different material combinations. Approximate solutions are first developed for FGM beams using iterative techniques based on averaging the elasticity modulus in beam’s tension/compression zones to determine different structural outputs like deflection, slope, shear, and moment diagrams. Then outputs are compared to classical Euler bending theory to validate the adopted approximate methods. In addition to beam bending problems, approximate solutions for FGM column buckling problems are developed. The methods adopted in solution are Rayleigh’s quotient, Timoshenko’s quotient, and Rayleigh-Ritz method. These methods are compared to the analytical classical Euler buckling solution. MATLAB program is adopted in the solution. The results of this study shed light on the importance of approximate solutions in solving FGM bending/buckling problems. For bending problems, the approximate method resulted in trivial error (~ 0%) when compared to analytical solution. As for buckling problems, Rayleigh-Ritz method was the most accurate in calculating critical buckling load with error less than 0.75%, while Rayleigh method led to significant error (~ 22%).

Introduction

The development of new materials is considered an engineering concern over time. The characteristics of the new engineering products are important for their use (Karamanli [1]). The stability of the new product is considered in engineering (Jarachi et al. [2]). The problems associated with new products depend on the properties of different original material components and distribution of constituents, or gradation, in different geometric directions. Functionally graded materials (FGMs) are characterized by the variation of composition and structure gradually over volume, resulting in corresponding changes in the properties of the materials. Therefore, FGMs can be designed for specific function or applications such as FG beam bending and FG column buckling. FGM was first introduced in 1984 in Japan in aerospace project (Niino et al. [3]). The project required production of a material with thickness 10 mm and tolerate outside/inside temperatures of 2,000 K/1,000 K. The efforts on the development of FGM continued thereafter. The activities of FGM development included civil engineering, mining, chemical, and biomedical products (Zhang et al. [4]). FGMs provided the opportunity to produce smooth materials with different chemical, mechanical, and physical properties such as Poisson’s ratio, Young’s modulus, coefficient of thermal expansion, density, and shear modulus. Moreover, the gradual and its variability provided new properties of FGM including stiffness ratio, corrosion resistivity, thermal conductivity, hardness, and specific heat (Allahyarzadeh et al. [5]).

Among the other works of functionally classified materials, various types of gradient composites were analyzed [6]. In addition, the authors reviewed reported potential applications. Spatial gradients in polymeric material structures due to varying chemical, mechanical, biomedical, and transport properties were studied by [7]. They indicated that structural gradients in the polymeric system could lead to desired gradient in single or combination of properties, as possible applications could be considered.

FGM activities in Japan were reviewed through different research projects, 1984–1997 by [8]. Many research interests were present at that time to develop new FGM to be used for different purposes, such as: thermal barriers, improve mechanical strength and thermal conductivity, relaxation of thermal stresses and energy conversion materials. Referring to (1987–1991) a five-year project in Japan by [9], which presented an overview that aimed to develop new FGMs that could relax thermal stresses effectively in high thermal environments.

Functionally classified materials are innovative composite materials that allow the deformation problems of classic composite materials to be overcome, especially if subjected to high thermal and mechanical stresses. For this reason, metal and ceramic are often combined to have a ductile, very resistant material with high thermal capacity. Development of a new type of ceramic metal assemblage which is designed for hot structural components, as a functionally gradient material is conducted by [10]. An overview of thirty-year research works available in the literature on FGM manufacturing methods, applications, and future challenges is conducted by [11].

It can be noticed from the study by [12] on boundary conditions for thermally loaded FG beam, the actual neutral surface is adopted as a base for FG Timoshenko beam analysis. Two distinct thermal loading cases are considered in the analysis based on original or simplified boundary conditions. It was found that the original boundary condition–based linear analysis for evaluation of equilibrium path of system yielded inaccurate trend.

Bending behavior of FG nanobeams via local/nonlocal stress gradient theory of elasticity is investigated by [13] for general boundary conditions using numerical analysis. Considering two loading cases, a parametric study is carried out.

Vibration analysis of FG nanobeams using surface stress-driven model is carried out by [14, 15]. Application of the surface stress-driven nonlocal theory of elasticity for the study of bending response of FG cracked nanobeams is conducted by [16]. FGM is now widely used to create nanoscale devices where other very important aspects of functionally classified materials come into play [17, 18].

Other works have also studied the effects of porosity in FGM as well as the hygrothermal effects coupled with porosity [19, 20, 21–22]. Nonlinear analysis of buckling, free vibration, and dynamic stability for the piezoelectric FG beams in thermal environment is carried out by [23]. Another research study on dynamic stability of nanobeams, which is based on nonlocal Timoshenko theory considering surface effects is done by [24].

In civil engineering, the material of beam/column affects bending of beams/buckling of columns [25]. In addition, the type of loading has a great effect on the response of the member [26, 27, 28, 29–30]. The results of the buckling load factor depend on the structure of the column and the type of FGMs used [2]. The wide variety of civil engineering products of FGM have increased the problems associated with the calculations. One problem is associated with stability of FGM products [31]. The stability of FGM structural members is related to several outputs like buckling stress, natural frequencies, and distribution of displacements [32]. This requires determining the mechanical/physical characteristics of any new FGM member and finding solutions for problems associated with tolerability and force distribution in complex structures [33]. Usually, the classical methods to solve FGM structural elements stability problems may be complicated since they were derived for isotropic materials with assumptions and simplifications.

Bending of beams has been investigated using different analytical methods such as Euler-Bernoulli method, Timoshenko method [1, 34, 35], and Reddy-Bickford method [1]. For buckling problem, different analytical methods were used to find solutions for FGM columns. The most familiar buckling result was calculated using Timoshenko beam theory (TBT). For example, the TBT is used to solve buckling of FG circular columns and shear deformation [36]. Another example where the TBT is used to find buckling that resulted from exposure to different temperatures is presented in [37]. Many other researchers followed these methods to solve FGM column problems [1, 2, 38, 39].

To simplify the analysis of FGM elements (beams/columns), approximate solutions should be developed. This can be done by following any of the approximate methods in structural analysis. For beam bending problems, the approximate solution can be developed by averaging the elasticity modulus above and below neutral axis. This can be done by trying several iterations with preliminary assumptions. As for column buckling problems, approximate solutions can be developed following three main methods: Rayleigh’s quotient, Timoshenko’s quotient, and Rayleigh-Ritz method. These methods can be efficient in solving FGM problems and reducing calculation time and complexity.

The aim of this research is to derive approximate solutions to FGM beam/column problems. To achieve this aim, analytical solutions to bending/buckling problems will be first introduced. Then, several approximate solutions will be derived and validated by comparing its outputs with the analytical solutions. This research will answer if approximate solutions be simple/accurate when analyzing bending of FG beams/buckling of FG columns? Are approximate solutions for bending/buckling problems valid under different boundary conditions/material properties? Several factors will be taken into consideration such as beams/columns dimensions, type of materials used, and boundary conditions. A MATLAB program is adopted to solve the equations and to compare the findings of this study with analytical methods to validate the results.

Analytical Solutions for FGM Members

Assumptions, Approximations, and Methodology

In analyzing beams and slender columns, the following assumptions are made:

1. The material of the element is homogeneous and isotropic.

2. Plan sections perpendicular to the neutral axis before bending remain plan and perpendicular to the neutral axis after bending

3. The stress–strain curve is identical in tension and compression.

4. No local type of instability is present.

5. The effect of transverse shear is negligible that is: shear deformation is negligible.

6. No applicable initial curvature exists.

7. The loads and bending moments act in a plane passing through a principal axis of inertia of the cross-section.

8. Hook’s law holds (stress is linearly proportional to strain).

9. The deflections are small as compared to the cross-sectional dimensions.

10. The loads are assumed coplanar, and applied quasi-statically, axial or transverse.

Even though there is normal stress concentration at regions of application the loads, normal stress in the transverse direction is neglected. As in this work, most researchers use this assumption in the analysis of FG beams since inclusion of normal stress and its variation in transverse direction has a minor effect on the solution. Poisson’s effect can be easily included in analysis through the definition of flexural rigidity. The variation of Poisson’s ratio in the transverse direction for FG beams/columns can be included in the analysis. In the current study, Poisson’s effect is not considered. Also, shear deformation is not included in the present approximate analysis of FG beams/columns. Euler-Bernoulli beam theory is used as another limitation of the current study. Unlike other research where Timoshenko beam theory (TBT), Reddy-Bickford beam theory (RBT), higher order shear deformation beam theories are used in the analysis. The effect of shear modulus of elasticity of materials used and their distribution function in the direction of material grading should be more pronounced for short beams/columns. Long beams and slender columns where length to thickness ratio is ten or more are used in the current study.

Even though, it is a well-known fact that elastic and geometric centers do not coincide for FG beams/columns in the transverse direction, many researchers provide solutions for differential equations of equilibrium considering bending about mid-axis rather than neutral axis. That is stresses and strains at the mid-axis are zero. This approximation is not part of the present study, where the actual location of neutral axis is found, and equilibrium of forces and moments are satisfied. In the present study, the main approximation is due to averaging the modulus of elasticity in the compression and tension zones to evaluate compressive and tensile forces and bending moment that satisfy equilibrium and compatibility of strains.

The methodology of the current study can be summarized by the flowchart of Fig. 1.

[See PDF for image]

Fig. 1

Flowchart of the work presented in this study

In this study, approximate solutions of bending of beams and buckling of columns made of FGMs are developed. Several examples are solved with different dimensions, boundary conditions, materials, and material gradations. Derivations of approximate methods, according to Euler-Bernoulli theory, are presented. Numerical routines, computations, and evaluation of results are performed using MATLAB. The approximate solutions assure the fact that neutral axis and mid-axis do not coincide. Therefore, the main step in the analysis of the cross-section of the member is to locate the neutral axis and the strain function accordingly. Equilibrium of axial tensile and compressive forces in addition to equilibrium of moments are satisfied. It should be noted here that in the solution adopted by many researchers that the neutral axis is assumed at mid-section of FG beams/columns. Even though the differential equation of equilibrium is solved exactly, the flexural rigidity has an error due to the shift of neutral axis from the mid-axis. This error can be relatively high for certain gradation of FG materials. It should be noted that assuming neutral axis at mid-section yields to loss of equilibrium of internal axial forces.

Neutral axis location should be determined by satisfying equilibrium of axial forces, then equilibrium of moment is satisfied and therefore flexural rigidity is evaluated. In this study, neutral axis is located using iterative technique to satisfy equilibrium of axial forces in an approximate manner. That is evaluating forces in compression and in tension using average values of modulus of elasticity above and below the neutral axis. The strain distribution is found by satisfying moment equilibrium at the section. Stress distribution is approximated to be linear in tension and in compression. Therefore, the approximate solution has an advantage of simplicity in evaluating the different structural parameters, unlike the exact solution that involves mathematical complexity due to the actual nonlinear distribution of the modulus of elasticity and therefore the nonlinear stress distribution. Due to the main similarities in the structural analysis of reinforced concrete beams/columns, those made of FG materials, structural engineers would find it easy and efficient to adopt the present solution or develop new solutions. New reinforced concrete materials can be developed in a way that can be analyzed as FG material using proper gradation functions.

Approximate solutions can be used to analyze FG beams/columns for their relatively good accuracy that is good for most practical purposes. Empirical formulas, design charts, interaction diagrams, etc. can be generated based on approximate solutions to facilitate the design of beams/columns made of FGMs.

In summary, the present study aims to highlight the importance of using approximate solutions. It assures the fact the neutral axis and mid-axis do not coincide and satisfy equilibrium of forces and moments, and compatibility of axial strains. It gives a better insight about the effect of different modulus of elasticity distribution functions and their impact on strain and stress functions in addition to internal forces and moments. Using MATLAB, several involving various parameters are solved and the results are solved and the results are compared to those in literature.

Beams Subjected to Bending

The major equations for FG beams are derived and used in the MATLAB. Consider a beam made of two materials as shown in Figs. 2 and 3. The beam is simply supported and subjected to a distributed load of intensity q₀ Fig. 2.

[See PDF for image]

Fig. 2

FG beam under uniformly distributed loading of length L, width b, and depth h

[See PDF for image]

Fig. 3

Strain, modulus of elasticity and stress distribution for FG beam in the transverse direction

For static equilibrium, the summation of axial forces acting on the beam must be zero. Hence, the following equation can be written:

1

$$\underset{\frac{-h}{2}}{\overset{\frac{h}{2}}{\int }}\sigma \left(z\right).dA=\underset{\frac{-h}{2}}{\overset{\frac{h}{2}}{\int }}\epsilon \left(z\right).E\left(z\right).b.dz=0$$where “$\sigma \left(z\right)$” is the stress at depth z, “dA” is the area of a small element in the beam as shown in Fig. 3, “E(x)” is the elasticity modulus at depth z, “$\epsilon \left(x\right)$” is the strain at depth z, and “b” is the beam width.

Another equation is related to the summation of internal moments in the section which must be equal to the external applied moment. In this case the moment is considered at the midspan and equals to $\frac{q_0{L}^2}{8}$. . Therefore, the following equation can be written:

2

$$\underset{\frac{-h}{2}}{\overset{\frac{h}{2}}{\int }}z.\epsilon \left(z\right).E\left(z\right).b.dz=\frac{q_0.L^2}{8}$$

Figure 3 shows strain, modulus of elasticity, and stress distributions in the transverse direction of the beam. Neutral axis of zero strain and zero stress is located Z₀ distance from top extreme fiber. Subscripts 1,2 denote top and bottom fibers respectively for material (MAT 1, MAT 2), strain, modulus of elastity, and stress values. The strain distribution is according to Euler-Bernoulli Theory.

According to Fig. 3, the compatibility of strain can be used as follows:

3

$$\frac{\epsilon \left(z\right)}{\widehat{z}}=\frac{{\epsilon }_1}{z_0}\to \epsilon \left(z\right)=\frac{\widehat{z}.{\epsilon }_1}{z_0}$$where “${\epsilon }_1$” is the maximum strain at the upper metal surface, “$z_0$” is the depth of neutral axis measured from top. “$\widehat{z}$” can be written as $z+z_0-\frac{h}{2}$. , substituting this in Eq. 3 leads to the following relation:

4

$$\epsilon \left(z\right)=\frac{{\epsilon }_1}{z_0}\left(z+z_0-\frac{h}{2}\right)$$

Another relation can be used which is the power law. This law relates the elasticity modulus at depth z “E(z)” to the elasticity modulus of both materials (upper and bottom surface materials) and reflects the uniformity omaterial distribution through the depth. The power law can be written as follows:

5

$$E\left(z\right)=E_2+\left(E_1-E_2\right){\left(\frac{z}{h}+\frac{1}{2}\right)}^k$$where “$E_1$” is the elasticity modulus of the upper surface material, and “$E_2$” is the elasticity modulus of the lower surface material. By substituting Eqs. 4 and 5 in Eqs. 1 and 2, and for k = 1, they become:

6

$$\underset{\frac{-h}{2}}{\overset{\frac{h}{2}}{\int }}\left[E_2+\left(E_1-E_2\right)\left(\frac{z}{h}+\frac{1}{2}\right)\right]\left[\frac{{\epsilon }_1}{z_0}\left(z+z_0-\frac{h}{2}\right)\right]b.dz=0$$

7

$$\underset{\frac{-h}{2}}{\overset{\frac{h}{2}}{\int }}z.\left[E_2+\left(E_1-E_2\right)\left(\frac{z}{h}+\frac{1}{2}\right)\right]\left[\frac{{\epsilon }_m}{z_0}\left(z+z_0-\frac{h}{2}\right)\right]b.dz=\frac{q_0.L^2}{8}$$

Solving Eqs. 6 and 7 leads to the magnitude of Z₀ and ${\epsilon }_1$. Therefore, the stress at level z can be determined as follows:

8

$$\sigma \left(z\right)=E\left(z\right).\epsilon \left(z\right)$$

As for deflection equations for FGB, they will be considered at the mid depth of the beam (z = 0). Therefore, the strain at this level will be:

9

$${\epsilon }_{\mathit{xx}}=\frac{d\theta }{\mathit{dx}}=-z\frac{d^{2}y}{d{x}^2}$$

According to Hook’s law:

10

$${\sigma }_{\mathit{xx}}=E\left(z\right).{\epsilon }_{\mathit{xx}}=-\left[E_2+\left(E_1-E_2\right){\left(\frac{z}{h}+\frac{1}{2}\right)}^k\right].z\frac{d^{2}y}{d{x}^2}$$

Substituting Eq. 5 in Eq. 7 leads to:

11

$$M_{\mathit{xx}}=-\underset{-\frac{h}{2}}{\overset{\frac{h}{2}}{\int }}\left[E_2+\left(E_1-E_2\right){\left(\frac{z}{h}+\frac{1}{2}\right)}^k\right]·z^2\frac{d^{2}y}{d{x}^2}dz=-D_{\mathit{xx}}\frac{d^{2}y}{d{x}^2}$$where D_xx is:

12

$$D_{\mathit{xx}}=\underset{-\frac{h}{2}}{\overset{\frac{h}{2}}{\int }}\left[E_2+\left(E_1-E_2\right){\left(\frac{z}{h}+\frac{1}{2}\right)}^k\right]·z^{2}dz$$

This will lead to the following differential equation:

13

$$\frac{d^2}{d{x}^2}\left(D_{\mathit{xx}}\frac{d^{2}y}{d{x}^2}\right)=w\left(x\right)$$where w(x) is the load function, “σ_xx” is the normal stress along the longitudinal axis of the beam (at z = 0), and “${\epsilon }_{\text{xx}}$” is the normal strain along the longitudinal axis of the beam (at z = 0).

Other methods can be used to determine the elasticity modulus of FGM beams such as exponential and sigmoidal methods. According to exponential method E(z) will be as follows:

14

$$E\left(z\right)=E_2{e}^{-\partial \left(1-\frac{2z}{h}\right)}$$

As for the sigmoidal method, E(z) will be as follows:

15

$$E(z)=\begin{array}{cc}\left(1-\frac{1}{2}{\left(\frac{\left[\frac{h}{2}\right]-z}{\frac{h}{2}}\right)}^p\right)E_2+\left(\frac{1}{2}{\left(\frac{\frac{h}{2}-z}{\frac{h}{2}}\right)}^p\right)E_1,& 0\le z\le \frac{h}{2}\\ \left(\frac{1}{2}{\left(\frac{\left[\frac{h}{2}\right]+z}{\frac{h}{2}}\right)}^p\right)E_2+\left(1-\frac{1}{2}{\left(\frac{\frac{h}{2}+z}{\frac{h}{2}}\right)}^p\right)E_1,& -\frac{h}{2}\le z\le 0\end{array}$$where “$\partial$” is the exponential gradation, “p” is the power coefficient.

Columns Subjected to Buckling

The solution for buckling problems in the case of FG columns is similar in form to isotropic columns. For the case of simply supported column (Fig. 4a), the lateral displacement and buckling load will be according to [40] as follows:

16

$$y(x)=d_z\left(cos\left(\sqrt{\frac{-P}{D_{\mathit{xx}}}x}\right)+\frac{1-cos\left(L\sqrt{-P/D_{\mathit{xx}}}\right)}{sin\left(L\sqrt{-P/D_{\mathit{xx}}}\right)}sin\left(\sqrt{\frac{-P}{D_{\mathit{xx}}}x}\right)-1\right)$$

17

$$P=\frac{{\pi }^2{D}_{\mathit{xx}}}{L^2}$$

[See PDF for image]

Fig. 4

a Simply supported and b cantilever beams subjected to axial compressive force P

As for the case of a cantilever column (Fig. 4b), the buckling load according to [40] will be as follows:

18

$$P=\frac{{\pi }^2{D}_{\mathit{xx}}}{4L^2}$$where

19

$$d_z=\frac{\left(E_1-E_2\right)hk}{2\left(k+2\right)\left(E_1+k{E}_2\right)}$$

Approximate Solutions for Bending and Buckling Problems

Beams Subjected to Bending

The derivation of approximate solution of beam bending is described in Fig. 5. This method is based on using average value for elasticity modulus of the upper surface and bottom surface materials. First, Using the compatibility of strains, the following equation can be written:

20

$${\epsilon }_2=\frac{h-Z_0}{Z_0}{\epsilon }_1$$

[See PDF for image]

Fig. 5

FG beam cross-section and strain, modulus of elasticity, approximate modulus of elasticity and approximate stress distributions in the transverse direction

By taking the equilibrium of forces (C and T), the following equation can be reached:

21

$${\left(h-Z_0\right)}^2=Z_0^2\frac{f_1{E}_1}{f_2{E}_2}$$

By taking the equilibrium of moment at the centerline, the following relation can be determined between Z₀ and ${\epsilon }_1$ as follows:

22

$${\epsilon }_1=\frac{3q_0{L}^2}{8f_1{E}_{1}bh{Z}_0}$$

The terms f₁E₁ and f₂E₂ can be determined as follows:

23

$$f_1{E}_1=\frac{{\int }_{\frac{h}{2}-Z_0}^{h/2}E\left(z\right)dz}{Z_0}$$

24

$$f_2{E}_2=\frac{{\int }_{-\frac{h}{2}}^{\frac{h}{2}-Z_0}E\left(z\right)dz}{h-Z_0}$$

E(z) can be specified using power, exponential, or sigmoidal functions as described earlier.

Figure 5 shows again the cross-section of FG beam and strain, modulus of elasticity, approximate modulus of elasticity and approximate stress distributions in the transverse direction.

Using a specialized software like the MATLAB, strain values and therefore stress and deformation functions can be determined as illustrated in the flowchart, Fig. 6. The neutral axis location is found using iterative technique. Approximate average values of moduli of elasticity distribution in compression and tension are adopted. Equilibrium of forces in the axial directions and moments in addition of compatibility of strains conditions are all satisfied accordingly.

[See PDF for image]

Fig. 6

Flowchart to determine the location of neutral axis approximately using iterative technique

The axial strain at top fiber as a function of the longitudinal direction is given by:

25

$${\epsilon }_1\left(x\right)=\frac{3q_0\left(Lx-x^2\right)}{2f_1{E}_1{Z}_{0}bh}$$

Therefore, for a simply supported FG beam, the curvature, slope and deflection are given by the following three equations respectively:

26

$$\frac{d^{2}W}{d{x}^2}=\frac{{\epsilon }_1\left(x\right)}{Z_0}=\frac{3q_0\left(Lx-x^2\right)}{2f_1{E}_1{Z}_0^{2}bh}$$

27

$$\frac{\mathit{dW}}{\mathit{dx}}=\frac{3q_0}{2f_1{E}_1{Z}_0^{2}h}\left(\frac{L{x}^2}{2}-\frac{x^3}{3}-L^3/12\right)$$

28

$$W\left(x\right)=\frac{q_0}{12b{f}_1{E}_1{h}^2{Z}_0^2}\left(-x^4+2L{x}^3-L^{3}x\right)$$

After locating the neutral axis at Z₀ distance from extreme top fiber using iterative technique, the approximate internal compressive (C) and tensile (T) forces, in addition to bending moment (M) at the middle section are evaluated according to the following:

29

$$C=\frac{1}{2}{f}_1{E}_1{\epsilon }_1{Z}_0$$

30

$$T=\frac{1}{2}{f}_2{E}_2{\epsilon }_2\left(h-Z_0\right)$$

31

$$M=\frac{2}{3}\left(C{Z}_0+T\left(h-Z_0\right)\right)$$

Columns Subjected to Buckling

There are many approximate solutions for buckling problems. Three main methods will be described: the Timoshenko approach, the Rayleigh approach, and the Rayleigh-Ritz approach.

First let’s derive Timoshenko’s approach equations [41]. Consider the column shown in Fig. 7. Assume “z” the deflection in the starting shape, “y” the additional deflection from the starting shape to the final shape, which means the total final deflection y + z, and “EI” is the flexural stiffness.

[See PDF for image]

Fig. 7

The parameters considered in deriving Timoshenko’s quotient method [41]

Assume that the bending moment has a value of $S=R.z$ at the starting shape. This will lead to $dS=dR.z$. Since an infinitesimal element with length dx is used in the derivation, the following relation can be written to determine the deflection:

32

$$-\frac{d^{2}y}{d{x}^2}dx=\left(\frac{R.z}{\mathit{EI}}+\frac{d^{2}z}{d{x}^2}\right)dx$$

Since the deflection at support levels is zero, this leads to the following relation:

33

$$dU=\underset{0}{\overset{l}{\int }}\left(\frac{R.z}{\mathit{EI}}+\frac{d^{2}z}{d{x}^2}\right)dx.z.dR=0$$where “dU” is the complimentary energy increment for the infinitesimal element. This leads to the following relation for the force “R”:

34

$$R=-\frac{{\int }_0^{l}z.\frac{d^{2}z}{d{x}^2}dx}{{\int }_0^l\frac{z^{2}dx}{\mathit{EI}}}$$

By choosing the parabola as a starting shape:

35

$$z=\frac{4c}{l^2}\left(lx-x^2\right)$$

The R will have the following value:

36

$$R=\frac{10EI}{l^2}$$

Timoshenko improved Eq. 34 to take the following form:

37

$$P=\frac{{\int }_0^l{\left(\frac{\mathit{dz}}{\mathit{dx}}\right)}^{2}dx}{{\int }_0^l\frac{z^{2}dx}{\mathit{EI}}}$$where “P” is Timoshenko’s variant of the energy method.

Now let’s derive the main equations for Ryleigh method [42]. Consider the column in Fig. 8 for the derivation. At initial stages and as the load is increased gradually before the buckling phenomena takes place, the strain energy due to axial deformation will be as follows:

38

$$E_s=\underset{0}{\overset{l}{\int }}\frac{1}{2}EA{\epsilon }^{2}dx$$

[See PDF for image]

Fig. 8

The parameters considered in deriving Rayleigh’s quotient method [42]

After buckling, bending strain energy is added as follows:

39

$$E_s=\underset{0}{\overset{l}{\int }}\frac{1}{2}EA{\epsilon }^{2}dx+\underset{0}{\overset{l}{\int }}\frac{1}{2}EI{k}^{2}dx$$

After buckling, the change in strain energy will result from the bending component only:

40

$$\mathrm{\Delta }{E}_s={\int }_0^l\frac{1}{2}EI{k}^{2}dx={\int }_0^l\frac{1}{2}EI{\left(\frac{d^{2}w}{{\mathit{dx}}^2}\right)}^{2}dx$$

This strain energy is equal to the work performed by the concentrated load (P) when the column is transferred from the straight position to the bending position. This work is $W=F\times u_F$ is horizontal displacement can be expressed in t vertical deflection using Pythagoras theorem as shown in Fig. 9.

[See PDF for image]

Fig. 9

The relation between horizontal and vertical deflection (Rayleigh’s quotient) [42]

Then the total horizontal displacement will be equal to:

41

$$u_f=\underset{0}{\overset{l}{\int }}\frac{1}{2}{\left(\frac{\mathit{dw}}{\mathit{dx}}\right)}^{2}dx$$

Equating strain energy due to bending to the potential energy done by axial compressive force results in the following equation:

42

$$\underset{0}{\overset{l}{\int }}\frac{1}{2}EI{\left(\frac{d^{2}w}{d{x}^2}\right)}^{2}dx=F\times \underset{0}{\overset{l}{\int }}\frac{1}{2}{\left(\frac{\mathit{dw}}{\mathit{dx}}\right)}^{2}dx$$where “F” is the buckling load. By assuming a displacement field w(x), we can find a buckling load which is usually higher than the real buckling load. The better the assumed displacement field the closer the buckling load will be to the real value.

As for Rayleigh-Ritz approa, the energy functional F is defined as the difference between the maximum kinetic energy and the maximum potential energy [43]. The deflection function W is chosen in terms of trial functions Øi as follows:

43

$$W=\sum _{i=1}^N{c}_i{\mathrm{\varnothing }}_i$$where “c_i” is an arbitrary coefficient and “Ø_i” is a function that can take any form (polynomial, exponential, trigonometric, etc.). The factor c_i can be determined by minimizing the energy functional F by takingartial derivatives with respect to c_i$\frac{\partial F}{\partial c_i}=0$. As for the trial function it can be expressed as follows:

44

$${\mathrm{\varnothing }}_i=f_i\prod _{J=1}^{\mathit{Ne}}{B}_J^{\mathit{SP}},i=1\to N$$where “f_i” is ith term of the complete set of polynomials, “Ne” is the edges number, “B_J” is the boundary expression of the jth edge, and “SP” is the suitable power.

For a beam, the maximum strain energy will be equal to:

45

$$E_s=\frac{\mathit{EI}}{2}\underset{0}{\overset{l}{\int }}{\left(\frac{d^{2}w}{d{x}^2}\right)}^{2}dx$$

The maximum kinetic energy will be equal to:

46

$$K_e=\frac{\rho A{\omega }^2}{2}\underset{0}{\overset{l}{\int }}{w}^{2}dx$$where “ρ” is the density, “A” is tharea, “s$\omega$” is the circular frequency. Assuming X = x/l and W = w/l, and equating maximum strain and kinetic energies we obtain:

47

$$\underset{0}{\overset{l}{\int }}{\left(\frac{d^{2}W}{d{X}^2}\right)}^{2}dX-{\mathrm{\Omega }}^2\underset{0}{\overset{l}{\int }}{W}^{2}dX=0$$where

48

$$\mathrm{\Omega }=\omega L^2\sqrt{\frac{\rho A}{\mathit{EI}}}$$

This will lead to the following eigen value problem:

49

$$\left[{\left[K\right]}_{N\times N}-{\mathrm{\Omega }}^2{\left[M\right]}_{N\times N}\right]{\left[C\right]}_{Nx1}=0$$where

50

$$k_{\mathit{ij}}=\underset{0}{\overset{l}{\int }}\frac{d^2{\mathrm{\varnothing }}_i}{d{X}^2}\times \frac{d^2{\mathrm{\varnothing }}_j}{d{X}^2}dX$$

51

$$m_{\mathit{ij}}=\underset{0}{\overset{l}{\int }}{\mathrm{\varnothing }}_i{\mathrm{\varnothing }}_{j}dX$$

MATLAB will be adopted to model each study case. Modeling inputs for every case can be provided through coding such as, span length, cross-sectional dimensions, applied loads, material gradation, etc. Then all required mathematical calculations for both analytical and approximate solutions can be entered into MATLAB using special functions such as integration functions, derivation functions, and differential equation functions. This will lead to the outputs to be extracted for each case. Different results can be gathered such as displacement, slope, bending moment, shear force, critical buckling load, etc. These results will come in different forms (tables/diagrams).

There are several variables/parameters that are considered in modeling and analysis of FG beams/columns. These parameters/variables defined in MATLAB include dimensions (length, depth and width), load pattern, boundary conditions, and material combinations. Defining these variables, MATLAB will analyze each problem using both analytical and approximate methods defined. Then, outputs can be compared to validate approximate solution methods and prove their efficiency in solving FGM problems.

Analysis and Discussion of Results

Dimensions of beams/columns are length L = 1, 2, 4 m, width b = 0.1 m, and depth h = 0.1 m. Several span—depth ratios will be considered in this study (L/h = 10, 20, and 40). A uniformly distributed load q₀ = 0.01 MN/m is considered for the different examples solved. Boundary condition of simply supported and cantilever FG beams and columns are considered in the analysis. Three material combinations Three material combinations are considered: Aluminum (Al), E₁ = 70 GPa—Silicon (Si), E₂ = 151 GPa, metal (SuS304), E₁ = 208 GPa—Zirconia (ZrO2_,), E₂ = 168 GPa, and Titanium alloy (Ti–6Al–4 V), E₁ = 110 GPa—Silicon nitride (Si3N4), E₂ = 322 GPa. The effect of different material gradation functions is considered: power, exponential, and sigmoidal, along with different power values: 0, 1, and 10. The findings from both analytical and approximate methods were compared to prove the efficiency of approximate solutions.

FG Beams

For different material combinations, the neutral axis position that is dimensionless values Z₀/h, are listed in Table 1, for different values of power (k = 0, 0.5, 1, 2, 5). The largest values for modulus of elasticity for all combinations is considered to be at the top.

Table 1. Dimensionless values for location of neutral axis (Z₀/h) for different power (k) and different material combinations

Material combinations	k
	0	0.5	1	2	5
Al–Si	0.5	0.4565	0.4389	0.4304	0.4423
SuS304–ZrO2	0.5	0.4863	0.4822	0.4816	0.4864
Ti-6Al–4 V–Si3N4	0.5	0.4438	0.4183	0.4022	0.4132

It can be noticed from Table 1 that, as the difference between the moduli of elasticity increases, the deviation of location of neutral axis from mid-axis increases. Material combination SuS304-ZrO2 has the least deviation of about 3.68% and Ti–6Al–4 V–Si3N4 has the most of about 19.56%. For all material combinations, as k increases the deviation increases up to a certain value of k < 5 then starts decreasing after k > 5 since the fast growth of elasticity function becomes bounded by a narrow region. Therefore, it would have a minor effect on axial force associated with the elasticity distribution which is smooth and almost linear for most of cross-Section (0.5 < k < 1).

Assuming neutral axis at the middle of the column and using the location of neutral axis at its actual position, Table 2 lists different estimations of flexural rigidity for the two cases (D_mid, D_act), in addition to the value of flexural rigidity according to the present approximate solution (D_app), for different material combinations. The power gradation function is used with k = 2. Therefore, D_mid, D_act, D_app are the flexural rigidities according to two the analytical solutions, in addition to the present approximate solution. The percentages difference between the various values are also listed in Table 2.

Table 2. Flexural rigidity values (D_mid, D_act, D_app) and percentage differences for different material combinations

Material combinations	Dmid (MN.m2)	Dact (MN.m2)	% Diff	Dapp (MN.m2)	% Diff, mid	% Diff, act
Al–Si	0.9883	1.0353	4.76	0.9901	0.18	− 4.37
SuS304–ZrO2	1.5167	1.5228	0.40	1.6031	5.39	5.27
Ti–6Al–4 V–Si3N4	1.9767	2.1495	8.78	1.9942	0.13	− 8.16

According to the present approximate solution, the top and bottom widths (b_t, b_b) for an equivalent isotropic T-section beam are calculated and presented in Table 3. They are calculated for the two cases of f₁E₁ and f₂E₂ values of moduli of elasticity. It is clear that when the modulus of elasticity for equivalent T-beam is f₁E₁, the width at the top remains at 100 mm, while the width at the bottom becomes less than 100 mm since the isotopic T-beam is stronger than the original FG beam. When using f₂E₂ modulus of elasticity for equivalent T-section beam, the bottom width remains at 100 mm, while the top width increases since the equivalent isotropic T-beam is weaker than the original FG beam. Using the same method, the top and bottom widths for equivalent T-section isotropic beam can be found for any value of modulus of elasticity. In this case, both widths would increase or decrease. For example, for E = 151 MPa, the values are reduced b_t = 89.67 mm and b_b = 70.06 mm, while for E = 70 MPa, the values are increased b_t = 193.42 mm and b_b = 151.14 mm.

Table 3. Top and bottom widths (b_t, b_b) of equivalent isotropic T-section beam, having different moduli of elasticity f₁E₁ and f₂E₂, for different material combinations according to the present approximate solution

Material combinations	Z0 (MN)	f1E1 (GPa)	bt (mm)	bb (mm)	f2E2 (GPa)	bt (mm)	bb (mm)
Al–Si	46.92	135.4	100.00	78.14	105.55	127.98	100.00
SuS304–ZrO2	48.97	200.5	100.00	92.09	184.90	108.50	100.00
Ti–6Al–4 V–Si3N4	46.05	279.3	100.00	72.86	203.90	137.25	100.00

Different Beam Theories

The values of maximum deflection, at mid-span of simply supported and at free end of cantilever FG beams, for different span to depth ratios, 10, 20, 40 are listed in Table 4. Present approximate solution values are compared with those of Euler Bernoulli Beam Theory EBT, Timoshenko and Goodier, TG [16], Higher Order Shear Deformation Theory, HSDT [11]. The two material constituents considered for the FG beam are: Aluminum (Al) with elasticity modulus of 70 GPa and Silicon (Si) with elasticity modulus of 151 GPa. The percentages difference of the approximate solution are: 4.37%, 1.99%, and 1.96% for the simply supported FG beam with L/h = 40, from EBT, HSDT, and TG respectively. For the cantilever FG beam the percentages difference are a little higher form HSDT and TG, 2.88% and 3.53, while the difference from EBT is almost the same.

Table 4. Maximum deflection using different theories and present approximate solution for FG simply supported and cantilever beams

Theory and Present	Simply supported beam maximum deflection (mm)			Cantilever beam maximum deflection (mm)
	L/h = 10	L/h = 20	L/h = 40	L/h = 10	L/h = 20	L/h = 40
EBT	0.12577	2.01230	32.19679	1.20738	19.3181	309.089
HSDT [11]	0.12890	2.06238	32.99813	1.22622	19.6194	313.911
TG [16]	0.12893	2.06288	33.00608	1.21805	19.4889	311.827
Present	0.13151	2.10402	33.66656	1.26250	20.2000	323.200

Effect of Span—Depth (L/h) Ratio

Three span—depth (L/h) ratios were considered (10, 20 and 40). Figure 10a depicts the deflection diagram at each distance ratio (x/L) measured for the simply supported beam. As span—depth ratio increases, the deflection of simply supported beam at any point increases. For example, the maximum deflection for the beam with (L/h = 40) was 33.7 mm, while for (L/h = 10) the maximum deflection was 0.13 mm.

[See PDF for image]

Fig. 10

a Deflection, b slope, c bending moment, and d shear for simply supported beam with different L/h ratios

Similar trend can be realized for slope results as shown in Fig. 10b. As the L/h ratio increases, the maximum slope at supports increases. For instance, the maximum slope at supports for the beam with L/h = 40 was 0.03 rad while for the beam with L/h = 20 it was 0.0036 rad. As we get closer to the beam mid-span, the slope decreases until reaching zero at mid-span.

As for the interaction force (bending moment and shear force), results showed that as the L/h ratio increases, both bending moment and shear force increased. The beam with L/h = 40 had a maximum moment of 0.021 MN.m compared to 0.0014 MN.m for the beam with L/h = 10 as depicted in Fig. 10c, and a maximum shear of 0.022 MN compared to 0.0054 MN for the beam with L/h = 10 as depicted in Fig. 10d.

Figures 11a, b show the strain and stress distribution for simply supported beams with different span—depth ratio. As the span—depth ratio increased, the maximum stress values at top and bottom surfaces increased. This is due to higher moment effects for larger spans which led to higher normal stress magnitudes.

[See PDF for image]

Fig. 11

a Strain and b stress distributions for simply supported beam with different L/h ratios

The same effect was studied for cantilevered beams. Although the behavior of cantilevered beams is different due to different boundary conditions, the same trend was realized for different L/h ratio. According to Fig. 12a, as the L/h ratio increases, the maximum deflection at the tip increases. For example, the maximum deflection for the beam with (L/h = 40) was 323 mm compared to 1.26 mm for the beam with (L/h = 10). As for slope results (Fig. 12b), the maximum slope at the tip showed a similar trend as the maximum deflection. For instance, the maximum slope at the tip for the beam with (L/h = 40) was -0.115 rad compared to − 0.0146 rad for the beam with (L/h = 20).

[See PDF for image]

Fig. 12

a Deflection, b slope, c bending moment, and d shear for fix-supported beam with different L/h ratios

Regarding interaction forces, results showed that as L/h ratio increases, both bending moment and shear force increase. According to Fig. 12c, the bending moment at support had a value of − 0.087 MN.m for the beam with (L/h = 40) and − 0.0054 MN.m for the beam with (L/h = 10). This trend was similar to that for shear force at support, where the beam with (L/h = 40) had maximum shear of 0.043 MN, while the beam with (L/h = 10) had a maximum shear of 0.011 MN as shown in Fig. 12d.

As for stress and strain results, a similar trend was realized for the cantilevered beams as the simply supported beams. It was shown that as the span—depth ratio increased, the maximum stress and strain values at top and bottom surfaces increased as shown in Fig. 13a, b.

[See PDF for image]

Fig. 13

a Strain and b stress distributions for cantilevered beam with different L/h ratios

The approximate solution suggested for FG beams mentioned in the study conducted by [1] proved to be close to the analytical solution. Figure 14a through d depict the stress distribution along the height of FG beam in [1] compared to the suggested approximate solution. It can be noted that the approximate solution results in closer stress values to the analytical solution. This was proved for both simply supported and cantilevered beams with different span—depth ratio as shown in Fig. 14a through d.

[See PDF for image]

Fig. 14

Stress distribution for simply supported/cantilever FG beam in [1] with (L/h = 10, 20) using analytical and approximate solutions

As for deflection results, the stiffness coefficient or flexural rigidity, D_xx, was determined for both analytical and approximate solution cases. D_xx was close in both cases; $210\times {10}^6$ by analytical method and $9.207\times {10}^6$ by approximate method. Having these values to be too close means that the deflection will be similar in both cases since deflection equation depends on the coefficient D_xx as shown below:

52

$$\mathrm{\delta }\left(\text{x}\right)=\frac{{\text{q}}_0{\text{L}}^4}{24D_{\mathit{xx}}}\left(\frac{\text{x}}{\text{L}}-\frac{2{\text{x}}^3}{{\text{L}}^3}+\frac{{\text{x}}^4}{{\text{L}}^4}\right)\to \text{for simply supported}$$

53

$$\mathrm{\delta }\left(\text{x}\right)=\frac{{\text{q}}_0{\text{L}}^4}{24D_{\mathit{xx}}}\left(\frac{6{\text{x}}^2}{{\text{L}}^2}-\frac{4{\text{x}}^3}{{\text{L}}^3}+\frac{{\text{x}}^4}{{\text{L}}^4}\right)\to \text{for cantilever}$$

Effect of Material Function Used

Three functions were adopted to determine the respoe of FGB: power function, exponential function, and sigmoidal function. Table 5 summarizes the maximum displacement results determined by analytical and approxime methods for both beams; simply supported at mid-span and cantilever at free end with (L/h = 20). The values of the power function are calculated for k = 1. It can be clearly seen that both methods showed similar results with slight differences (less than 2%) regardless of the function used in describing material gradation and for different boundary conditions. This demonstrates the validity of the suggested approximate method to be used in determining different beam results.

Table 5. Analytical and approximate maximum deflection for FG simply supported and cantilever beam, according to different gradation functions

Function	Cantilever beam		Simply supported beam
	${\delta }_{max-Analytical}$. (mm)	${\mathrm{\delta }}_{\text{max}-\text{Approximate}}$(mm)	${\mathrm{\delta }}_{\text{max}-\text{Analytical}}$(mm)	${\mathrm{\delta }}_{\text{max}-\text{Approximate}}$(mm)
Power	21.716	.933	2.262	2.283
Exponential	21.692	22.126	2.260	2.284
Sigmoidal	21.707	21.967	2.260	2.283

Effect of Power Value (k)

Figures 15a through d show the effect of different power values on the stress and strain distribution in both simply supported beams and cantilevered beams with (L/h = 20). For a power value of k = 0, the material will be isotropic and has an equivalent elasticity modulus of (E₁) as defined earlier. As the power value increases, the concentration of E₁ material decreases and the concentration of E₂ material increases. This will lead to a stiff behavior for the beam. For example, for k = 0, the beam had higher stress compared to the case where k = 1 (equal distribution of materials). This was shown for both simply supported and cantilever beam. As the power value increases beyond 1, the maximum stress value decreased since the concentration of E₂ material is higher than E₁ material (E₂ > E₁). This was reflected also on the deflection results shown in Fig. 15e, where the displacement for the beam with power value = 10 was less than the beam with power values 0 and 1.

[See PDF for image]

Fig. 15

Strain, stress distributions, and deflection for simply supported/cantilever beam with different power values

Regarding approximate method results, Table 6 summarizes the maximum displacement for FG beams with different power values using both analytical and approximate methods. It can be noted that the suggested approximate method leads to maximum displacement values that are too close to those determined by the analytical method. This again proves the validity of the suggested approximate methods to be applied on FG bending problems.

Table 6. Analytical and approximate maximum deflection for different values of the power gradation function (k = 0, 1, and 10)

Function	Simply supported beam
	${\mathrm{\delta }}_{\text{max}-\text{Analytical}}$(mm)	${\mathrm{\delta }}_{\text{max}-\text{Approximate}}$(mm)
k = 1	2.262	2.283
k = 0	3.571	3.642
k = 10	1.850	1.876

Effect of Different Material Combination

Other material combinations were considered to have a better insight about the behavior of FG beams. The material parameters were taken from research conducted by [44]. The materials were Zirconia (ZrO₂), Metal (SuS304), Titanium alloy (Ti–6Al–4 V), and Silicon Nitride (Si3N4). These materials have different properties than the materials used before.

Figures 16a through f summarize the main results for this new material combination compared to the previous one for simply supported beams. In terms of deflection (Fig. 16a), results showed that new combinations showed less deflection than the one adopted before, (Al & Si). This is because the new materials have higher elasticity modulus than the previous ones (E_(ZrO2) = 168 GPa, E_(SuS304) = 208 GPa, E_(Si3N4) = 322 GPa, and E_{(Ti-6Al-4 V)} = 110 GPa). Similar trend was realized for other findings like slope (Fig. 16b), bending moment (Fig. 16c), shear force (Fig. 16d), strain (Fig. 16e) and stress (Fig. 16f). The reduction percentage reached 48.8% for maximum deflection, maximum slope, maximum bending moment, and maximum shear. Similarly, for stress and strain values, they were 30–35% less than those for the previous material combination used.

[See PDF for image]

Fig. 16

Deflection, slope, bending moment, shear, strain, and stress for simply supported beam with different material combinations

As for the approximate method suggested, Figs. 17a through d summarize strain and stress distribution results for FG beams with SuS304–ZrO2 and Ti–6Al–4 V-Si3N4. Both simply supported and cantilever beams were considered. It can be noted that the approximate solution led to a close strain and stress distribution to the analytical solution. The difference between the maximum stresses was trivial and reached 1% in simply supported cases.

[See PDF for image]

Fig. 17

Approximate vs analytical stress distribution for SuS304–ZrO2/ Ti–6Al–4 V–Si3N4 FG simply supported beam/cantilever beam

As for deflection curves, both methods led to similar curves since the stiffness constant D_xx (flexural rigidity) obtained from approximate solution, D_app = $1.60MN.m^2$. was close to that obtained from analytical solutions. For instance, D_xx for FG beam with SuS304–ZrO2 was $1.56MN.m^2$ obtained from the analytical solution [1]. As for FG beam with Ti–6Al–4 V–Si3N4 combination, D_app = $1.97MN.m$, while the flexural rigidity was $1.87MN.m$ as estimated by [1]. Again, this provethe validity of the suggested approximate solution for FG beams.

FG Columns

As mentioned earlier, two column cases were considered in the buckling analysis; the simply supported column and the cantilever column. Both cases were studied at three span—depth ratios: 10, 20, and 40. The compression load considered in this study was 1 MN. In addition, a comparison was done on the FG material models (power, exponential, and sigmoidal) to check the difference in the results among them. Furthermore, approximate solutions were compared to analytical solutions to confirm its validity for FG column buckling problems.

Effect of Span—Depth (L/h) Ratio

Three span—depth (L/h) ratios were considered in this study: 10, 20 and 40. Figure 18a depicts the displacement diagram at each distance ratio (x/L) measured for the simply supported column. As span—depth ratio increases, the deflection of simply supported column at any point increases. For example, the maximum displacement for the column with (L/h = 40) was 1.62 mm, while for (L/h = 10) the maximum deflection was 0.09 mm.

[See PDF for image]

Fig. 18

lateral displacement and buckling load for simply supported columns with different L/h ratios

As for the critical buckling load, it was shown that as the span—depth ratio increased, the buckling load decreased as depicted in Fig. 18b. The load decreased from 90 MN for L/h = 10 to 5.7 MN for L/h = 40.

Figure 19a depicts the displacement diagram at each distance ratio (x/L) measured for the cantilevered column. As span—depth ratio increases, the deflection of cantilevered column at any point increases. For example, the maximum displacement for the column with (L/h = 40) was 18.3 mm, while for (L/h = 10) the maximum deflection was 0.35 mm.

[See PDF for image]

Fig. 19

lateral displacement and buckling load for cantilevered columns with different L/h ratios

As for the critical buckling load, it was shown that as the span – depth ratio increased, the buckling load decreased as depicted in Fig. 19b. The load decreased from 22.7 MN for L/h = 10 to 1.4 MN for L/h = 40.

As for approximate methods adopted to determine the buckling load for FG columns, Figs. 20a through d show the critical buckling load determined by analytical method, Rayleigh method, Timoshenko method, and Rayleigh-Ritz method. It can be depicted that for both simply supported and cantilever columns, the Rayleigh-Ritz method has resulted in the most accurate critical buckling load compared to the analytical method. In average, Rayleigh quotient method led to an error percentage of 22% for simply supported column and 1.5% for cantilever column compared to the analytical solution. As for Timoshenko quotient method, the error percentage was 1.4% for simply supported columns and 0.8% for cantilevered columns. Finally, for Rayleigh-Ritz method, the error percentage was the least among all approximate methods and was around 0.7% for both boundary conditions.

[See PDF for image]

Fig. 20

Buckling load for simply supported/cantilever columns with (L/h = 10, 20) using different approximate method solutions

Effect of Material Function Used

Three functions were adopted to determine the response of FGB: power function, exponential function, and sigmoidal function. Table 7. Summarizes critical buckling load and maximum deflection results for both simply supported and cantilever columns with (L/h = 20). It can be clearly seen that all methods showed similar results with slight differences (less than 1%) for both simply supported and cantilevered columns. This demonstrates the validity for these methods to be used in determining different column results.

Table 7. Critical buckling load and maximum deflection of simply supported and cantilever columns using different material gradation functions

Function	Cantilever column		Simply supported column
	Pcr (MN)	${\delta }_{\mathit{max}}\left(\text{mm}\right)$	Pcr (MN)	${\delta }_{\mathit{max}}\left(\text{mm}\right)$
Power	5.6754	1.618022	22.7017	0.347325
Exponential	5.6816	1.615883	22.7264	0.3469301
Sigmoidal	5.6785	1.616952	22.7140	0.347127

The same columns were analyzed by approximate solution methods with different materials functions. Table 8 shows that all material functions have resulted in similar buckling loads for FG columns with error less than 1%. This proved the validity of different materials functions in approximate solution methods for buckling problems.

Table 8. Critical buckling load using approximate methods for different material gradation functions

Function	Pcr—Cantilever column (MN)			Pcr—Simply supported column (MN)
	RQ	TQ	RR	RQ	TQ	RR
Power	5.76	5.72	5.72	27.6	23.0	22.9
Exponential	5.78	5.73	5.73	27.8	22.2	23.0
Sigmoidal	5.77	5.73	5.72	27.7	23.0	23.1

Effect of Power Value (k)

The power value demonstrates the distribution of materials among the column length. For power constant = 0, the material will be isotropic with E = E₁. As for power = 1, both materials will have homogeneous distribution (material 1 will have 100% intensity at the bottom of the column and reduces gradually until reaching 0% at the top of the column). Figures 21a and b summarize lateral displacement and critical buckling load values for simply supported columns with power values = 0, 1, and 10.

[See PDF for image]

Fig. 21

lateral displacement and buckling load for simply supported columns with different power values

According to Fig. 21a, as the power (k) increases, the lateral displacement decreases. This is because the material with E₂ will have higher concentration in the column, and E₂ > E₁, which means the column will have higher stiffness since the (EI) will be greater as the value of k increases. This will lead to a less lateral displacement for the same applied compression force. This, on the other hand, means that the critical buckling load will increase as the power value increases since the column now had higher stiffness.

Approximate solution methods were applied on FG columns with different power values. Table 9 summarizes critical buckling load values for FG columns with power values = 0, 1, and 10. It can be noted that approximate method results were close to those obtained from the analytical method. In addition, the Rayleigh-Ritz method was the most accurate method compared to other approximate solution methods with error percentage of approximately 0.75%.

Table 9. Critical buckling load using approximate methods for different power values

Function	Pcr—Simply supported column (MN)
	Analytical	Rayleigh	Timoshenko	Rayleigh–Ritz
k = 0	14.4	17.6	14.6	14.5
k = 1	22.7	27.6	23.0	22.9
k = 10	27.8	33.4	28.0	27.9

Effect of Different Material Combination

Two other material combinations were investigated to validate approximate methods. The material combinations are metal (SuS304)—ceramic (ZrO2), and Titanium (Ti–6Al–4 V)—silicon nitride (Si3N4).

As the applied compressive force at mid-axis of FG column increases, the transverse deformation increases in the pre-buckling stage, because of the imperfectness due to eccentricity (e = h/2–Z₀), the distance between the mid-axis and the neutral axis. The maximum mid-span deformations of simply supported FG column made of different combinations of materials are listed in Table 10 Even though the flexural rigidity of Al–Si is less than that of SuS304–ZrO2, δ_max is more because it has higher eccentricity due to higher difference of moduli of elasticity (E₁–E₂). Ti–6Al–4 V–Si3N4 combination a higher eccentricity and a higher (E₁–E₂), therefore it has a much higher δ_max. Hence bending behavior is more pronounced before buckling. Values of mid-span deflection in Table 10 have been generated by super imposing the deflection due to moment (M = Pe) and moment (M = P ${\mathrm{\delta }}_{\left(\text{X}\right)\text{l}}$) using the Second Moment Area Theorem. For the three material combinations. This rigorous estimation yields a higher value.

Table 10. Maximum mid-span deflection due to increasing axial compressive force for different material combinations in the pre-buckling stage

P(MN)	δmax (mm)
	Al-Si	SuS304-ZrO2	Ti-6Al-4 V-Si3N4
2	2.019	0.343	1.247
4	4.714	0.769	2.716
6	8.086	1.278	4.404
8	12.134	1.869	6.324
10	17.158	2.543	8.443
(Pcr, Al-Si) = 10.218	17.415	–	–
12		3.548	10.773
14		4.140	13.364
(Pcr, SuS304-ZrO2) = 15.029	–	4.604	–
16			15.996
18			19.243
20			22.400
(Pcr, Ti-6Al-4 V-Si3N4) = 21.215	–	–	24.471

It can be noticed in Table 11, that the results of the present approximate method are in good agreement with the analytical solutions, however, they deviate more for the Ti–6Al–4 V–Si3N4 since the difference in moduli of material components is higher and therefore the gradation function has higher rate of change. Since all approximate methods yield values of P_cr, which are in marginal deviation from the analytical values, any of them can be considered in the stability analysis and design of FG slender columns subjected to axial compression. For the Raleigh quotient method, it is recommended to use functions that satisfy natural boundary conditions in addition the essential ones.

Table 11. Critical buckling load for simply supported column with different material combinations using various approximate and analytical solutions

Material combinations	(Pcr)act (MN)	(Pcr)mid (MN)	(Pcr)app (MN)	(Pcr)TQ (MN)	(Pcr)QQ (MN)	(Pcr)RR (MN)
Al–Si	10.218	9.754	9.772	10.353	12.424	10.295
SuS304–ZrO2	15.029	14.969	15.822	15.228	18.273	15.142
Ti–6Al–4 V-Si3N4	21.215	19.509	19.485	21.495	25.794	21.374

It can be noted that the location of neutral axis and therefore the eccentricity is about 10% of the thickness for Ti–6Al–4 V–Si3N4 combination, which is more compared to that of SuS304–ZrO2, which is just about 2%. This relatively high difference in eccentricity is due to the high difference of moduli of elasticity of constituents, 212 GPa and 40 GPa for the above two combinations, respectively. The approximate solution deviates from the analytical solution due to the difference in eccentricity even though the flexural rigidities difference is much less.

The critical buckling load, P_cr, is calculated for different materials combinations of FG slender columns using various analytical and approximate solutions, Table 11. (P_cr)_act is the critical buckling load using neutral axis at its actual locations, while (P_cr)_mid is the value when adopting the analytical solution using the assumption that mid-axis is the neutral axis. (P_cr)_app is the critical buckling load using the present approximate solution. The Timoshenko quotient, Raleigh quotient, and the Raleigh-Ritz approximations of critical buckling load are (P_cr)_TQ, (P_cr)_RQ, and (P_cr)_RR respectively. The values of P_cr for the various solutions are according to a power function gradation in the transverse direction with power (k = 2). The columns have the same dimensions: L = 1 m, b = h = 0.1 m. Table 11 lists the values of P_cr for columns with simply supported boundary conditions.

Table 12 lists the values of Z₀, e, and the maximum deflection at mid-span just before buckling, δ_cr, for different material combinations. The values tabulated are for both analytical and approximate solutions.

Table 12. Values of Z₀, e, δ_cr (mm) using analytical and approximate solutions for different material combinations

Material combinations	Analytical			Approximate
	Z0 (mm)	e (mm)	δcr (mm)	Z0 (mm)	e (mm)	δcr (mm)
Al–Si	43.04	6.96	68.692	46.92	3.08	30.39
SuS304–ZrO2	48.16	1.84	18.166	48.97	1.03	16.46
Ti–6Al–4 V–Si3N4	40.22	9.78	96.525	46.05	3.95	38.98

Conclusions and Recommendations

Modeling and analysis of FG beams/columns are conducted to compare approximate methods and analytical methods for bending and buckling problems. It is concluded:

1. The suggested approximate solution method for bending FG beams proved to work efficiently. The method was tried for different boundary conditions (simply supported and fixed) and with different combination of materials. Different outputs were included in the comparison between the two methods like: bending moment, shear force, deflection, slope, stress distribution, and strain distribution. The error in all outputs was trivial (~ 0%) and proved the validity of the approximate method suggested for FG beams.

2. It was proven by both approximate and analytical solution methods that using high stiffness materials in FG beams has enhanced its structural behavior. This was shown by comparing maximum deflection and maximum moment for FG beams with different material combination, where beams with higher elasticity modulus materials had less maximum deflection and maximum moment compared to beams with less elasticity modulus materials.

3. The approximate solutions derived for buckling analysis of FG columns proved to work efficiently. This was checked by comparing the critical buckling load of different FG columns with different boundary/materials conditions. Among all approximate methods tried in this work, the Rayleigh–Ritz method proved to be the most accurate with an error percentage of less than 0.75%. On the other side, Rayleigh methods seemed to be the least accurate with an error percentage of 22%.

4. It was proved by both approximate and analytical solution methods that using high stiffness materials in FG columns has increased its critical buckling load, where FG columns with higher elasticity modulus materials had critical buckling load compared to columns with less elasticity modulus materials.

5. Approximate solutions proved to be effective in saving time/effort in analysis, as well as being good enough for practical purpose.

Based on the findings of this research, it is recommended to develop approximate solutions for dynamic problems on FG beams/columns and for other structural members such as FG walls/slabs. Also, Modeling FG beams/columns using finite element software to extract different complex outputs.

Many practical applications can be considered by the approximate method presented in this study. The method can be used to solve other material combinations and other distribution functions. It can be used to solve FG beams/columns subject to other types of loads; mechanical, thermal, static, dynamic, with various boundary conditions. For its simplicity and efficiency, structural and mechanical engineers can use the method in structural analysis and design of FG beams/columns. For civil engineers, the method can be implemented easily due to the basic similarities in the analysis methodology of reinforced concrete beam/column sections and that of FG adopted in this study. Design charts, empirical formulas, interaction diagrams, etc. can be generated using this method to be used in design of structural members made of FGMs. For certain steel distribution within fiber reinforced concrete, reinforced concrete members can be modeled as FG one. The approximate method presented provides the structural engineer with a fast technique to determine equivalent isotropic beams or columns to those made of FGMs by varying the cross-section dimensions or/and the moduli of elasticity.

The accuracy of the approximate method as results indicated, is sufficient for most practical purposes. Thus, it can be considered as a reliable method to be used to analyze and design FG beams/columns. The effect of the method provides a better insight of different parameters involves in the analysis from physical and mathematical points of view. The iterative technique used seems to be efficient and logically sound since convergence occurs fast for all different cases solved. The location of neutral axis that satisfies equilibrium and compatibility conditions is located for all cases attempted, even for higher order power, exponential, and sigmoidal functions.

Acknowledgements

The authors would like to express their gratitude to the Civil Engineering departments of both Isra University and the UAE University for making this research possible.

Author Contributions

All authors contributed to the preparation of the text and presentation of the results. They also read and approved the final manuscript.

Funding

The authors have not disclosed any funding.

Data Availability

Data sharing is not applicable to this article as no datasets were generated or analyzed during the current study.

Declaration

Conflict of interest

The authors declare no conflict of interest and the work is not submitted elsewhere.

Publisher's Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

References

1. Karamanli, A. Analysis of bending deflections of functionally graded beams by using different beam theories and symmetric smoothed particle hydrodynamics. IJET; 2016; 2, 3 pp. 105-117. [DOI: https://dx.doi.org/10.19072/ijet.259394]

2. Jarachi, O; Abidi, M; Cherradi, T. Model identification for the evaluation of critical buckling load in reinforced concrete rectangular columns. Int. J. GEOMATE; 2019; 17, 63 pp. 103-110. [DOI: https://dx.doi.org/10.21660/2019.63.17223]

3. Niino, M; Hirai, T; Watanabe, R. The functionally gradient materials. J. Jap. Soc. Compos. Mat.; 1987; 13, pp. 257-264. [DOI: https://dx.doi.org/10.6089/jscm.13.257]

4. Zhang, N; Khan, T; Guo, H; Shi, S; Zhong, W; Zhang, W. Functionally graded materials: an overview of stability, buckling, and free vibration analysis. Adv. Mater. Sci. Eng. Hindawi; 2019; 2019, pp. 1-18. [DOI: https://dx.doi.org/10.1155/2019/1354150]

5. Allahyarzadeh, M; Aliofkhazraei, M; Rouhaghdam, A; Torabinejad, V. Gradient electrodeposition of Ni–Cu–W (alumina) nanocomposite coating. Mater. Des.; 2016; 107, pp. 74-81. [DOI: https://dx.doi.org/10.1016/j.matdes.2016.06.019]

6. M. B. Bever, P. E. Duwez, On gradient composites, preliminary reports, memoranda and technical notes of the ARPA materials summer conference, 117–140 (1970)

7. Shen, M; Bever, MB. Gradients in polymeric materials. J. Mater. Sci.; 1972; 7, 7 pp. 741-746. [DOI: https://dx.doi.org/10.1007/BF00549902]

8. Koizumi, MFGM. FGM activities in Japan. Compos. B Eng.; 1997; 28, 1–2 pp. 1-4. [DOI: https://dx.doi.org/10.1016/S1359-8368(96)00016-9]

9. Koizumi, M; Niino, M. Overview of FGM research in Japan. MRS Bull.; 1995; 20, 1995 pp. 19-24. [DOI: https://dx.doi.org/10.1557/S0883769400048867]

10. K. Wakashima, T. Hirano, M. Niino, Functionally gradient materials (FGM) architecture a new type of ceramic/metal assem- blage designed for hot structural components, Proc. Intl. Symp. Space Appl. Adv. Struct. Mater. (1990)

11. Saleh, B; Jiang, J; Fathi, R; Al-hababi, T; Xu, Q; Wang, L; Song, D; Ma, A. 30 years of functionally graded materials: an overview of manufacturing methods. Appl. Future Chall. Compos. Part B Eng.; 2020; 201, 108376. [DOI: https://dx.doi.org/10.1016/j.compositesb.2020.108376]

12. Dehrouyeh-Semnani, AM. On boundary conditions for thermally loaded FG beams. Int. J. Eng. Sci.; 2017; 119, pp. 109-127.3679759 [DOI: https://dx.doi.org/10.1016/j.ijengsci.2017.06.017]

13. Penna, R; Lambiase, A; Lovisi, G; Feo, L. Investigating hygrothermal bending behavior of FG nanobeams via local/nonlocal stress gradient theory of elasticity with general boundary conditions. Mech. Adv. Mater. Struct.; 2023; [DOI: https://dx.doi.org/10.1080/15376494.2023.2269938]

14. Lovisi, G; Feo, L; Lambiase, A; Penna, R. Application of surface stress-driven model for higher vibration modes of functionally graded nanobeams. Nanomaterials; 2024; 14, 350. [DOI: https://dx.doi.org/10.3390/nano14040350]

15. Feo, L; Lovisi, G; Penna, R. Free vibration analysis of functionally graded nanobeams based on surface stress-driven nonlocal model. Mech. Adv. Mater. Struct.; 2023; [DOI: https://dx.doi.org/10.1080/15376494.2023.2289079]

16. Lovisi, G. Application of the surface stress-driven nonlocal theory of elasticity for the study of the bending response of FG cracked nanobeams. Compos. Struct.; 2023; 324, 117549 pp. 0263-8223. [DOI: https://dx.doi.org/10.1016/j.compstruct.2023.117549]

17. Qiu, L; Zhu, N; Feng, Y; Michaelides, EE; Yła, GZ; Jing, D; Zhang, X; Norris, PM; Markides, CN; Mahian, OA. Review of recent advances in thermophysical properties at the nanoscale: FROM solid state to colloids. Phys. Rep.; 2020; 843, pp. 1-81. [DOI: https://dx.doi.org/10.1016/j.physrep.2019.12.001]

18. Łydz ba, D; Róz an ski, A; Sevostianov, I; Stefaniuk, D. Principle of equivalent microstructure in micromechanics and its connection with the replacement relations, thermal conductivity problem. Int. J. Eng. Sci.; 2019; 144, 103126. [DOI: https://dx.doi.org/10.1016/j.ijengsci.2019.103126]

19. Imani, A; Rabczuk, A; Friswell, T; Ian, MA. Finite element model for the thermo-elastic analysis of functionally graded porous nanobeams. Eur. J. Mech. A. Solids; 2019; 77, 103767.3973972 [DOI: https://dx.doi.org/10.1016/j.euromechsol.2019.04.002]

20. Ebrahimi, F; Jafari, A. A higher-order thermomechanical vibration analysis of temperature-dependent FGM beams with porosities. J. Eng.; 2016; 2016, 9561504. [DOI: https://dx.doi.org/10.1155/2016/9561504]

21. Penna, R; Feo, L; Lovisi, G; Fabbrocino, F. Hygro-thermal vibrations of porous FG nano-beams based on local/nonlocal stress gradient theory of elasticity. Nanomaterials; 2021; 11, 910. [DOI: https://dx.doi.org/10.3390/nano11040910]

22. Penna, R; Feo, L; Lovisi, G; Fabbrocino, F. Application of the higher-order Hamilton approach to the nonlinear free vibrations analysis of porous FG nano-beams in a hygrothermal environment based on a local/nonlocal stress gradient model of elasticity. Nanomaterials; 2022; 12, 2098. [DOI: https://dx.doi.org/10.3390/nano12122098]

23. Fu, Y; Wang, J; Mao, Y. Nonlinear analysis of buckling, free vibration and dynamic stability for the piezoelectric functionally graded beams in thermal environment. Appl. Math. Model.; 2012; 36, pp. 4324-4340.2929833 [DOI: https://dx.doi.org/10.1016/j.apm.2011.11.059]

24. Saffari, Sh; Hashemian, M; Toghraie, D. Dynamic stability of functionally graded nanobeam based on nonlocal Timoshenko theory considering surface effects. Physica. B; 2017; 520, pp. 97-105. [DOI: https://dx.doi.org/10.1016/j.physb.2017.06.029]

25. Alshafei, A. FE modeling and analysis of isotropic and orthotropic beams using first order shear deformation theory. Mater. Sci. Appl.; 2013; 4, 1 pp. 77-102. [DOI: https://dx.doi.org/10.4236/msa.2013.41010]

26. Temsah, Y; Jahami, A; Khatib, J; Sonebi, M. Numerical analysis of a reinforced concrete beam under blast loading. MATEC Web Conf.; 2018; 149, 02063. [DOI: https://dx.doi.org/10.1051/matecconf/201814902063]

27. Temsah, Y; Jahami, A; Khatib, J; Sonebi, M. Numerical derivation of ISO-damaged curve for a reinforced concrete beam subjected to blast loading. MATEC Web Conf.; 2018; 149, 02016. [DOI: https://dx.doi.org/10.1051/matecconf/201714902016]

28. Jahami, A; Temsah, Y; Khatib, J. The efficiency of using CFRP as a strengthening technique for reinforced concrete beams subjected to blast loading. Int. J. Adv. Struct. Eng.; 2019; 11, 4 pp. 411-420. [DOI: https://dx.doi.org/10.1007/s40091-019-00242-w]

29. Jahami, A; Temsah, Y; Khatib, J; Baalbaki, O; Darwiche, M; Chaaban, S. Impact behavior of rehabilitated post-tensioned slabs previously damaged by impact loading. Mag. Civ. Eng.; 2020; 93, 1 pp. 134-146. [DOI: https://dx.doi.org/10.18720/MCE.93.11]

30. Jahami, A; Temsah, Y; Khatib, J; Baalbaki, O; Kenai, S. The behavior of CFRP strengthened RC beams subjected to blast loading. Mag. Civ. Eng.; 2021; 103, 3 10309. [DOI: https://dx.doi.org/10.34910/MCE.103.9]

31. Matsunaga, H. Free vibration and stability of functionally graded circular cylindrical shells according to a 2D higher order deformation theory. Compos. Struct.; 2009; 88, 4 pp. 519-531. [DOI: https://dx.doi.org/10.1016/j.compstruct.2008.05.019]

32. Sofiyev, H. The vibration and stability behavior of freely supported FGM conical shells subjected to external pressure. Compos. Struct.; 2009; 89, 3 pp. 356-366. [DOI: https://dx.doi.org/10.1016/j.compstruct.2008.08.010]

33. Swaminathan, K; Naveenkumar, D. Higher order refined computational models for the stability analysis of FGM plates-analytical solutions. Eur. J. Mech.-A/Solids; 2014; 47, pp. 349-361.3233509 [DOI: https://dx.doi.org/10.1016/j.euromechsol.2014.06.003]

34. Huang, Y; Li, X. Buckling of functionally graded circular columns including shear deformation. Mater. Des.; 2010; 31, pp. 3159-3166. [DOI: https://dx.doi.org/10.1016/j.matdes.2010.02.032]

35. Huang, Y; Li, X. Buckling analysis of nonuniform and axially graded columns with varying flexural rigidity. J. Eng. Mech.; 2011; 137, 1 pp. 73-81. [DOI: https://dx.doi.org/10.1061/(ASCE)EM.1943-7889.000020]

36. Huang, Y; Wu, J; Li, X; Yang, L. Higher-order theory for bending and vibration of beams with circular cross section. J. Eng. Math.; 2013; 80, pp. 91-104.3055630 [DOI: https://dx.doi.org/10.1007/s10665-013-9620-2]

37. Kiani, Y; Eslami, M. Thermomechanical buckling of temperature dependent FGM beams. Lat. Am. J. Solids Struct.; 2013; 10, 2 pp. 223-246. [DOI: https://dx.doi.org/10.1590/S1679-78252013000200001]

38. Storch, J; Elishakoff, I. Buckling of axially graded columns: a fifth-order polynomial mode shape. AIAA J.; 2018; 56, 6 pp. 1-5. [DOI: https://dx.doi.org/10.2514/1.J056488]

39. Alshabatat, N. Optimal design of functionally graded material columns for buckling problems. J. Mech. Eng. Sci.; 2018; 12, 3 pp. 3914-3926. [DOI: https://dx.doi.org/10.15282/jmes.12.3.2018.11.0]

40. S. Taeprasartsit, A buckling analysis of perfect and imperfect functionally graded columns. In: Proceedings of the institution of mechanical engineers, part L: J. Mater. Design Appl. 226(1), 16–33 (2011), https://doi.org/10.1177/14644207114266

41. Westergaard, H. On the method of complementary energy. Trans. Am. Soc. Civ. Eng.; 1942; 107, 1 pp. 765-793. [DOI: https://dx.doi.org/10.1061/TACEAT.0005550]

42. I. Welleman, 2014. [Online] Icozct.tudelft.nl. Available at: <https://icozct.tudelft.nl/TUD_CT/CT3109/collegestof/arbeid_en_energie/files/Exercise8.pdf> [5 Accessed March 2021].

43. Kumar, Y. The Rayleigh-Ritz method for linear dynamic, static and buckling behavior of beams, shells and plates: a literature review. J. Vib. Control; 2017; 24, 7 pp. 1205-1227.3776284 [DOI: https://dx.doi.org/10.1177/10775463176947]

44. Kadoli, R; Akhtar, K; Ganesan, N. Static analysis of functionally graded beams using higher order shear deformation theory. Appl. Math. Model.; 2008; 32, 12 pp. 2509-2525. [DOI: https://dx.doi.org/10.1016/j.apm.2007.09.015]