Examination of Beam Theories for Buckling and Free Vibration of Functionally Graded Porous Beams

Abstract

This paper examines the accuracy and effectiveness of various beam theories in predicting the critical buckling loads and fundamental frequencies of functionally graded porous (FGP) beams whose material properties change continuously across the thickness. The beam theories considered are classical beam theory (CBT), first-order shear deformation beam theory (FSDBT), third-order shear deformation beam theory (TSDBT), and the broken-line hypothesis-based shear deformation beam theory (BSDBT). Governing equations for those beam theories are formulated by using the Hamilton’s principle and are then solved by means of the generalised differential quadrature method. Finite element simulation solutions are provided as reference results to assess the predictions of those beam theories. Comprehensive numerical results are presented to evaluate the influences of the porosity distribution and coefficient, slenderness ratio, and boundary condition on the difference between theoretical predictions and simulation results. It is found that the differences significantly increase as the porosity coefficient rises, and this effect becomes more noticeable for the rigid beam with a smaller slenderness ratio. Nonetheless, the results produced by the BSDBT are always the closest to simulation ones. The findings in this paper will contribute to the establishment of more refined theories for the mechanical analysis of FGP structures.

1. Introduction

Porous materials are a class of substances that contain numerous interconnected or enclosed pores or voids inside, with the key features of high porosity, low density, large specific surface area, and high specific strength. These characteristics contribute to their lightweight, good sound and heat insulation, high permeability, and energy adsorption capacity, making porous materials suitable for widespread applications in aerospace, automotive, construction, and other industries [1,2,3]. The development of porous materials has therefore become an active area of research, with scientists exploring new methods to prepare and modify porous materials with desired properties for specific applications.

One of recent advances in porous materials is the functionally graded porous materials (FGPMs), in which the porosity (pore morphology or density) shows continuous change along certain direction(s) so as to meet multifunctional requirements. FGPMs have a wide range of potential applications in aerospace, automotive, and biomedical engineering. For example, Deng et al. [4] designed a leaf vein-inspired heat exchanger filled with graded porous structures, exhibiting an improved heat transfer efficiency compared to the uniform counterpart. Zhang et al. [5] proposed a bi-graded foam-filled structure to improve the crashworthiness of vehicles. Perez-Boerema et al. [6] designed a femoral prosthesis characterised by functionally graded porous (FGP) structures as the bone tissue substitute. When such advanced materials are assembled into engineering systems as structural components (e.g., beams, plates, shells, etc.), their mechanical analysis is of great importance for further applications of FGPMs in practical engineering. So far, considerable research efforts have been dedicated to the buckling and free vibration analyses of FGP beams [7,8]. Chen et al. [9,10,11,12] systematically evaluated the buckling and vibration behaviours of FGP beams and indicated that their mechanical performances can be further improved via the symmetric nonuniform porosity distribution. Wattanasakulpong et al. [13,14] conducted the free vibration analysis of FGP beams based on the Chebyshev collocation method. Ebrahimi et al. [15,16] successively investigated effects of rotation and thermal loading on the vibration behaviour of FGP beams. Al Rjoub and Hamad [17] examined the free vibration of FGP beams using the transfer matrix method. Jamshidi and Arghavan [18] studied the optimal distributions of porosity for buckling and free vibration behaviour of FGP beams. Anirudh et al. [19] presented a comprehensive analysis on the buckling, bending, and free vibration of FGP beams by the finite element (FE) approach. Heshmati and Daneshmand [20] assessed the influence of different profile variations of vibrational properties of nonuniform FGP beams. Liu et al. [21] performed a thermal–mechanical coupling buckling analysis of FGP sandwich beams. Askari et al. [22] dealt with the free vibration of FGP beams with piezoelectric layers. Jena et al. [23] analysed the vibration behaviour of FGP beams resting on an elastic foundation. The buckling and free vibration problems of FGP beams were also revisited recently by other researchers based on different beam theories and solution methods [24,25,26,27,28], which are not elaborated here further to avoid repetition.

It is noted that various beam theories are used for buckling and vibration analyses of FGP beams in the abovementioned studies, categorised as the classical beam theory (Euler–Bernoulli beam theory, CBT), first-order shear deformation beam theory (Timoshenko beam theory, FSDBT), and higher-order shear deformation beam theories (HSDBTs), among which the HSDBTs include the widely used third-order shear deformation beam theory (TSDBT) of Reddy, sinusoidal shear deformation beam theory (SSDBT), and exponential shear deformation theory (ESDBT).

Table 1. Statistics on the beam theories used in the literature (the blank cell indicates a duplication of the content above).

Theories		Years	References	Mechanical Behaviours
CBT		2016	[16]	vibration
		2017	[17]
		2021	[22]
FSDBT		2015	[9]	buckling, bending
			[13]	vibration
			[15]
		2016	[10]
		2017	[11] [17]	buckling, vibration vibration
		2018	[18]	buckling, vibration
		2019	[20]	vibration
		2021	[22]
		2022 2023	[12] [27]	buckling, vibration
HSDBT	TSDBT	2018 2021	[14] [22]	vibration
		2022	[24]	buckling, bending, vibration
		2023	[25]
		2024	[26] [27]	buckling, vibration vibration
	SSDBT	2019	[19]	buckling, bending, vibration
			[21]	buckling
		2021	[22]	vibration
		2021	[23]
	ESDBT	2021	[22]	vibration

presents the statistics on the beam theories used in the literature. Those theories have proven their accuracy and effectiveness in analysing the mechanical behaviours of isotropic, homogenous, laminated composite and FGM beams. Distinct deformation hypotheses of those theories contribute to the differences of solutions in predicting mechanical behaviours, and the differences become negligible when the beam is slender enough (normally the slenderness ratio is over 10). However, Magnucki and Stasiewicz [29] pointed out earlier that the FGP beam with the maximum/minimum porosity in the middle/outer surfaces will behave like a sandwich beam, and the CBT and FSDBT cannot correctly determine the displacements of the cross-section of the beam. Instead, they proposed a new shear deformation beam theory by making use of the broken-line hypothesis (BSDBT) applied to three-layered structures. Unfortunately, there is an absence of a comprehensive analysis of the accuracy and effectiveness of the BSDBT, and comparative studies between conventional beam theories (CBT, FSDBT, and HSDBT) and the BSDBT in buckling and free vibration analyses of FGP beams have not been reported yet.

Therefore, this paper is devoted to examining the validity and accuracy of four beam theories, namely the CBT, FSDBT, TSDBT, and BSDBT, in analysing the buckling and free vibration behaviours of FGP beams with uniform and nonuniform distributions of porosity. Within the framework of Hamilton’s principle, governing equations of buckling and free vibration are deduced based on each individual beam theory and are all solved by using the generalised differential quadrature (GDQ) method. Comparative studies are conducted through comprehensive numerical examples, with a particular focus on the effects of porosity distribution and coefficient, beam geometry, and boundary conditions on the differences between results of those beam theories. In addition, the accuracy and effectiveness of those beam theories are evaluated by comparing the obtained buckling loads and natural frequencies with FE simulation results. Lastly, the underlying reasons for the discrepancies observed between the results of those theories and simulation are discussed. The present work will serve to better understand the differences and limitations of existing theories for buckling and free vibration analyses of FGP beams and contribute to the establishment of more refined theories for the mechanical analysis of FGP structures.

2. Functionally Graded Porous Beam

Figure 1 shows a porous beam of length L and thickness h defined in a rectangular coordinate system. Three porosity gradients that are most commonly investigated in the literature are considered in this study. Among those, patterns FGX/FGO imply that the porosity continuously decreases/increases from the middle to the upper and lower surfaces, while pattern UD means that the porosity is uniformly distributed in the beam. As a result, the elastic moduli either change or remain constant along the thickness direction, which are defined as follows.

$$E(z)=E_{\mathrm{s}}[1-e_0\phi (z)],G(z)=G_{\mathrm{s}}[1-e_0\phi (z)]$$

(1)

in which

$$\phi (z)=\{\begin{array}{l}\cos \left(\pi z/h\right)\hspace{1em}\hspace{1em}\hspace{1em} (\mathrm{FGX})\\ \cos \left(\left|\pi z/h\right|-\pi /2\right) \hspace{1em}(\mathrm{FGO})\\ {\phi }_0\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}(\mathrm{UD})\end{array}$$

(2)

E_s and G_s are the Young’s modulus and shear modulus of the solid material without pores, respectively; e₀ is the porosity coefficient defined by

$$e_0=1-E_{\mathrm{p}}/E_{\mathrm{s}}=1-G_{\mathrm{p}}/G_{\mathrm{s}}$$

(3)

where E_P and G_P are the elastic moduli of the layer with the highest porosity, namely the minimum values of E(z) and G(z), respectively. According to the relationship between Young’s modulus and shear modulus, the Poisson’s ratio is given by

$$\nu (z)=E(z)/2G(z)-1$$

(4)

which is apparently a constant for the entire beam.

Several micromechanical models [7] have been developed for the evaluation of mechanical properties of open- and closed-cell porous materials. Among those, the Gibson–Ashby model [30] is the most extensively applied due to its simple form and wide applicability. According to this model, the Young’s modulus of porous materials is related to the mass density as

$$\frac{E(z)}{E_{\mathrm{s}}}={\left(\frac{\rho (z)}{{\rho }_{\mathrm{s}}}\right)}^2$$

(5)

where ρ_s is the mass density of the solid material. By combining Equations (1) and (5), the mass density of the porous beam can be expressed as

$$\rho (z)={\rho }_{\mathrm{s}}\sqrt{1-e_0\phi (z)}$$

(6)

For direct comparisons, the total mass of porous beams with the same dimension but different porosity gradients should be equal, which requires

$${\displaystyle {\int }_{-h/2}^{h/2}{\rho }_{\mathrm{s}}\sqrt{1-e_0\phi (z)}}\mathrm{d}z={\displaystyle {\int }_{-h/2}^{h/2}{\rho }_{\mathrm{s}}\sqrt{1-e_0{\phi }_0}}\mathrm{d}z$$

(7)

from which the constant coefficient ${\phi }_0$ for the uniform porosity distribution can be derived as

$${\phi }_0=\frac{1}{e_0}\left\{1-{\left[\frac{1}{h}{\displaystyle {\int }_{-h/2}^{h/2}\sqrt{1-e_0\phi (z)}}\mathrm{d}z\right]}^2\right\}$$

(8)

For the sake of simplicity, the mass density of the porous beam is rewritten in a similar form with the elastic moduli in Equation (1) as

$$\rho (z)={\rho }_{\mathrm{s}}[1-e_{\mathrm{m}}\phi (z)]$$

(9)

in which e_m can be determined by comparing Equations (6) and (9) as

$$e_{\mathrm{m}}=\frac{1}{\phi (z)}\left[1-\sqrt{1-e_0\phi (z)}\right]$$

(10)

3. Theoretical Modelling

3.1. Beam Theories

In the present study, four beam theories, namely the CBT, FSDBT, TSDBT, and BSDBT, which are most often used in existing studies, are examined for buckling and free vibration analysis of FGP beams. Since different deformation assumptions are applied, the displacement fields of those beam theories differ from each other, as shown in Figure 2.

In the CBT, it is assumed that the transverse normal remains perpendicular to the midplane after deformation, which means that the transverse shear deformation is omitted. In this case, the displacement filed is given by

$$U(x,z,t)=-z\frac{\partial w(x,t)}{\partial x},W(x,z,t)=w(x,t)$$

(11)

In contrast, the FSDT takes into account the transverse shear deformation by adopting that the transverse normal needs not be perpendicular to the midplane but the cross-section remains a planar surface after deformation. Under this assumption, the displacement field of FSDBT is of the form

$$U(x,z,t)=z\psi (x,t),W(x,z,t)=w(x,t)$$

(12)

However, the assumption on the straightness and normality of a transverse normal after deformation is relaxed in the TSDBT and BSDBT by expanding the displacements as cubic and trigonometric functions, respectively, of the thickness coordinate. Accordingly, the displacement fields are expressed as

$$U(x,z,t)=z\psi (x,t)-\frac{4z^3}{3h^2}\left(\psi (x,t)+\frac{\partial w(x,t)}{\partial x}\right),W(x,z,t)=w(x,t)$$

(13)

for the TSDBT and

$$\{\begin{array}{l}U(x,z,t)=-z\frac{\partial w(x,t)}{\partial x}+\frac{h}{\pi }\left[{\varphi }_1(x,t)\sin \left(\pi z/h\right)+{\varphi }_2(x,t)\sin \left(2\pi z/h\right){\cos }^2\left(\pi z/h\right)\right]\\ W(x,z,t)=w(x,t)\end{array}$$

(14)

for the BSDBT which is established by making use of the broken-line hypothesis applied to three-layered structures [29]. In the above equations, U and W denote displacements of the beam along the x and z directions, respectively; w represents the transverse displacement component in the midplane (z = 0); and ψ is the transverse normal rotation about the x-axis. ϕ₁ and ϕ₂ are the dimensionless functions of displacements, and t is the time. By setting ϕ₁ = ϕ₂ = 0, the BSDBT is reduced to the CBT.

According to the linear strain–displacement relations, the longitudinal tensile and transverse shear strains (ε_x, γ_xz) of different theories are given as

$$\mathrm{CBT}:{\epsilon }_x=-z\frac{{\partial }^{2}w}{\partial x^2},{\gamma }_{xz}=0$$

(15)

$$\mathrm{FSDBT}:{\epsilon }_x=z\frac{\partial \psi }{\partial x},{\gamma }_{xz}=\psi +\frac{\partial w}{\partial x}$$

(16)

$$\mathrm{TSDBT}:{\epsilon }_x=z\frac{\partial \psi }{\partial x}-c_1{z}^3\left(\frac{\partial \psi }{\partial x}+\frac{{\partial }^{2}w}{\partial x^2}\right),{\gamma }_{xz}=\psi +\frac{\partial w}{\partial x}-c_2{z}^2\left(\psi +\frac{\partial w}{\partial x}\right)$$

(17)

$$\mathrm{BSDBT}:\{\begin{array}{l}{\epsilon }_x=-z\frac{{\partial }^{2}w}{\partial x^2}+\frac{h}{\pi }\left[\frac{\partial {\varphi }_1}{\partial x}\sin \left(\pi z/h\right)+\frac{\partial {\varphi }_2}{\partial x}\sin \left(2\pi z/h\right){\cos }^2\left(\pi z/h\right)\right]\\ {\gamma }_{xz}={\varphi }_1\cos \left(\pi z/h\right)+{\varphi }_2\left[\cos \left(2\pi z/h\right)+\cos \left(4\pi z/h\right)\right]\end{array}$$

(18)

in which c₁ = 4/3h² and c₂ = 4/h².

3.2. Governing Equations

Consider an FGP beam with the rectangular coordinate system (x, z) as shown in Figure 1. The governing equations of the beam under a uniform compressive axial load F can be formulated by using Hamilton’s principle

$$\delta {\displaystyle {\int }_{t_1}^{t_2}\left[{\mathsf{\prod }}_K-\left({\mathsf{\prod }}_S+{\mathsf{\prod }}_W\right)\right]}\hspace{0.17em}\mathrm{d}t=0$$

(19)

in which the kinetic energy П_K, strain energy П_S, and virtual work П_W performed by the compressive load are given by

$${\mathsf{\prod }}_K=\frac{1}{2}{\displaystyle {\int }_0^L{\displaystyle {\int }_{-h/2}^{h/2}\rho }}\left[{\left(\frac{\partial U}{\partial t}\right)}^2+{\left(\frac{\partial W}{\partial t}\right)}^2\right]\mathrm{d}z\mathrm{d}x$$

(20)

$${\mathsf{\prod }}_S=\frac{1}{2}{\displaystyle {\int }_0^L{\displaystyle {\int }_{-h/2}^{h/2}(E{\epsilon }_x^2+G{\gamma }_{xz}^2)}}\mathrm{d}z\mathrm{d}x$$

(21)

$${\mathsf{\prod }}_W=-\frac{1}{2}{{\displaystyle {\int }_0^{L}F\left(\frac{\partial W}{\partial x}\right)}}^2\mathrm{d}x$$

(22)

Substituting Equations (20)–(22) into Equation (19) and integrating the parts over the beam thickness, the governing equations for each beam theory are then derived and written in terms of displacement components as

$$D_{11}\frac{{\partial }^{4}w}{\partial x^4}+F\frac{{\partial }^{2}w}{\partial x^2}=-I_0\ddot{w}+I_2\frac{{\partial }^2\ddot{w}}{\partial x^2}$$

(23)

for the CBT,

$$\kappa A_{55}\left(\frac{{\partial }^{2}w}{\partial x^2}+\frac{\partial \psi }{\partial x}\right)-F\frac{{\partial }^{2}w}{\partial x^2}=I_0\ddot{w}$$

(24)

$$D_{11}\frac{{\partial }^2\psi }{\partial x^2}-\kappa A_{55}\left(\frac{\partial w}{\partial x}+\psi \right)=I_2\ddot{\psi }$$

(25)

for the FSDBT, and

$$\begin{array}{l}\left(c_1{H}_{11}-{c_1}^2{J}_{11}\right)\frac{{\partial }^3\psi }{\partial x^3}-{c_1}^2{J}_{11}\frac{{\partial }^{4}w}{\partial x^4}+\left(A_{55}-2c_2{B}_{55}+{c_2}^2{D}_{55}\right)\left(\frac{\partial \psi }{\partial x}+\frac{{\partial }^{2}w}{\partial x^2}\right)-F\frac{{\partial }^{2}w}{\partial x^2}\\ =I_0\ddot{w}+\left(c_1{I}_4-{c_1}^2{I}_6\right)\frac{\partial \ddot{\psi }}{\partial x}-{c_1}^2{I}_6\frac{{\partial }^2\ddot{w}}{\partial x^2}\end{array}$$

(26)

$$\begin{array}{l}\left(D_{11}-2c_1{H}_{11}+{c_1}^2{J}_{11}\right)\frac{{\partial }^2\psi }{\partial x^2}+\left({c_1}^2{J}_{11}-c_1{H}_{11}\right)\frac{{\partial }^{3}w}{\partial x^3}+\left(2c_2{B}_{55}-A_{55}-{c_2}^2{D}_{55}\right)\left(\psi +\frac{\partial w}{\partial x}\right)\\ =\left(I_2-2c_1{I}_4+{c_1}^2{I}_6\right)\ddot{\psi }+\left({c_1}^2{I}_6-c_1{I}_4\right)\frac{\partial \ddot{w}}{\partial x}\end{array}$$

(27)

for the TSDBT. In the above equations, κ = 5/6 is the shear correction factor for the FSDBT; the overdot represents the partial derivative with respect to time; and the stiffness elements A_ij, B_ij, D_ij, F_ij, H_ij, and J_ij and inertia-related term I_i are defined as

$$\left(A_{11},B_{11},D_{11},F_{11},H_{11},J_{11}\right)={\displaystyle {\int }_{-h/2}^{h/2}E(z)(1,z,z^2,z^3,z^4,z^6)\hspace{0.17em}}\mathrm{d}z$$

(28)

$$(A_{55},B_{55},D_{55})={\displaystyle {\int }_{-h/2}^{h/2}G(z)(1,z^2,z^4)}\hspace{0.17em}\mathrm{d}z$$

(29)

$$I_i={\displaystyle {\int }_{-h/2}^{h/2}\rho (z){(z)}^i}\mathrm{d}z,\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}(i=0,1,2,\cdots ,6)$$

(30)

For the BSDBT, the governing equations are expressed as

$$E_{\mathrm{s}}{h}^3\left(C_1\frac{{\partial }^{4}w}{\partial x^4}-C_2\frac{{\partial }^3{\varphi }_1}{\partial x^3}-C_3\frac{{\partial }^3{\varphi }_2}{\partial x^3}\right)+F\frac{{\partial }^{2}w}{\partial x^2}={\rho }_{\mathrm{s}}{h}^3\left(B_1\frac{{\partial }^2\ddot{w}}{\partial x^2}-\frac{I_0}{h^3}\ddot{w}-B_2\frac{\partial {\ddot{\varphi }}_1}{\partial x}-B_3\frac{\partial {\ddot{\varphi }}_2}{\partial x}\right)$$

(31)

$$E_{\mathrm{s}}\left(C_2\frac{{\partial }^{3}w}{\partial x^3}-C_4\frac{{\partial }^2{\varphi }_1}{\partial x^2}-C_5\frac{{\partial }^2{\varphi }_2}{\partial x^2}\right)+\frac{E_{\mathrm{s}}}{2h^2\left(1+v\right)}\left(C_7{\varphi }_1+C_8{\varphi }_2\right)={\rho }_{\mathrm{s}}\left(B_2\frac{\partial \ddot{w}}{\partial x}-B_4{\ddot{\varphi }}_1-B_5{\ddot{\varphi }}_2\right)$$

(32)

$$E_{\mathrm{s}}\left(C_3\frac{{\partial }^{3}w}{\partial x^3}-C_5\frac{{\partial }^2{\varphi }_1}{\partial x^2}-C_6\frac{{\partial }^2{\varphi }_2}{\partial x^2}\right)+\frac{E_{\mathrm{s}}}{2h^2\left(1+v\right)}\left(C_8{\varphi }_1+C_9{\varphi }_2\right)={\rho }_{\mathrm{s}}\left(B_3\frac{\partial \ddot{w}}{\partial x}-B_5{\ddot{\varphi }}_1-B_6{\ddot{\varphi }}_2\right)$$

(33)

in which

$$B_i=\frac{1}{h}{\displaystyle {\int }_{-h/2}^{h/2}\left[1-e_m\phi (z)\right]{\xi }_i(z)\mathrm{d}z},C_i=\frac{1}{h}{\displaystyle {\int }_{-h/2}^{h/2}\left[1-e_0\phi (z)\right]{\xi }_i(z)\mathrm{d}z\hspace{0.17em}\hspace{0.17em}(i=0,1,2,\cdots ,9)}$$

(34)

By setting $\eta =z/h$, ξ_i in Equation (34) are given by

$$\begin{gathered} {\xi }_1={\eta }^2,{\xi }_2=\eta \sin \left(\pi \eta \right),{\xi }_3=\frac{\eta }{\pi }\sin \left(2\pi \eta \right){\cos }^2\left(\pi \eta \right),{\xi }_4=\frac{1}{{\pi }^2}{\sin }^2\left(\pi \eta \right) \\ {\xi }_5=\frac{1}{{\pi }^2}\sin \left(\pi \eta \right){\cos }^2\left(\pi \eta \right)\sin \left(2\pi \eta \right),{\xi }_6=\frac{1}{{\pi }^2}{\sin }^2\left(2\pi \eta \right){\cos }^4\left(\pi \eta \right),{\xi }_7={\cos }^2\left(\pi \eta \right) \\ {\xi }_8=\cos \left(\pi \eta \right)\left[\cos \left(2\pi \eta \right)+\cos \left(4\pi \eta \right)\right],{\xi }_9={\left[\cos \left(2\pi \eta \right)+\cos \left(4\pi \eta \right)\right]}^2 \end{gathered}$$

(35)

Two different end supports, namely hinged (H) and clamped (C), are considered in this study. The corresponding boundary conditions in terms of displacements are given in

Table 2. Boundary conditions for different beam theories.

Beam Theory	End Support	Boundary Conditions
CBT	Hinged	$w=0$$,\frac{{\partial }^{2}w}{\partial x^2}=0$
	Clamped	$w=0$$,\frac{\partial w}{\partial x}=0$
FSDBT	Hinged	$w=0$$,\frac{\partial \psi }{\partial x}=0$
	Clamped	$w=0$$,\psi =0$
TSDBT	Hinged	$w=0$$,\frac{{\partial }^{2}w}{\partial x^2}=0$$,\frac{\partial \psi }{\partial x}=0$
	Clamped	$w=0$$,\frac{\partial w}{\partial x}=0$$,\psi =0$
BSDBT	Hinged	$w=0$$,\frac{{\partial }^{2}w}{\partial x^2}=0$$,\frac{\partial {\varphi }_k}{\partial x}=0$ (k = 1, 2)
	Clamped	$w=0$$,\frac{\partial w}{\partial x}=0$$,{\varphi }_k=0$ (k = 1, 2)

.

3.3. Solution Procedure

Since it is difficult to obtain closed-form solutions, a numerical method such as the GDQ method [31,32] is adopted to solve the abovementioned differential governing equations. According to this method, the displacement components w, ψ, ϕ₁, and ϕ₂ and their m^th partial derivatives with respect to x are approximated by

$$\{\begin{array}{l}{\left\{w,\psi ,{\varphi }_1,{\varphi }_2\right\}|}_{x=x_i}={\displaystyle \sum _{j=1}^N{l}_j(x_i)}\left\{w_j,{\psi }_j,\hspace{0.17em}{\varphi }_{1j},\hspace{0.17em}{\varphi }_{2j}\right\}\\ {\frac{{\partial }^m}{\partial x^m}\left\{w,\hspace{0.17em}\psi ,\varphi ,{\varphi }_1,\hspace{0.17em}{\varphi }_2\right\}|}_{x=x_i}={\displaystyle \sum _{j=1}^N{G}_{ij}^{(m)}}\left\{w_j,\hspace{0.17em}{\psi }_j,{\varphi }_{1j},\hspace{0.17em}{\varphi }_{2j}\right\}\end{array}$$

(36)

where {w_j, ψ_j, ϕ_1j, ϕ_2j} are the values of {w, ψ, ϕ₁, ϕ₂} at x = x_j; l_j(x) are the Lagrange interpolation polynomials; and $G_{ij}^{\left(m\right)}$ are the weighting coefficients of the m^th derivatives and are determined by the recursive formula [33]. N is the total number of grid points which are spaced along the x-axis according to the following pattern:

$$x_i=\frac{1}{2}\left[1-\cos \frac{\pi (i-1)}{N-1}\right]\hspace{1em}(i=1,2,\dots ,N)$$

(37)

By applying the relationship in (36), the differential governing equations of each beam theory can be discretised as

$$D_{11}{\displaystyle \sum _{j=1}^N{G}_{ij}^{(4)}}{w}_j+F{\displaystyle \sum _{j=1}^N{G}_{ij}^{(2)}}{w}_j=-I_0{\ddot{w}}_i+I_2{\displaystyle \sum _{j=1}^N{G}_{ij}^{(2)}}{\ddot{w}}_j$$

(38)

for the CBT,

$$A_{55}({\displaystyle \sum _{j=1}^N{G}_{ij}^{(2)}}{w}_j+{\displaystyle \sum _{j=1}^N{G}_{ij}^{(1)}}{\psi }_j)-F{\displaystyle \sum _{j=1}^N{G}_{ij}^{(2)}}{w}_j=I_0{\ddot{w}}_i$$

(39)

$$D_{11}{\displaystyle \sum _{j=1}^N{G}_{ij}^{(2)}}{\psi }_j-A_{55}({\displaystyle \sum _{j=1}^N{G}_{ij}^{(1)}}{w}_j+{\psi }_i)=I_2{\ddot{\psi }}_i$$

(40)

for the FSDBT,

$$\begin{array}{l}\left(c_1{H}_{11}-{c_1}^2{J}_{11}\right){\displaystyle \sum _{j=1}^N{G}_{ij}^{(3)}{\psi }_j}-{c_1}^2{J}_{11}{\displaystyle \sum _{j=1}^N{G}_{ij}^{(4)}{w}_j}+\left(A_{55}-2c_2{B}_{55}+{c_2}^2{D}_{55}\right)\left({\displaystyle \sum _{j=1}^N{G}_{ij}^{(1)}{\psi }_j}+{\displaystyle \sum _{j=1}^N{G}_{ij}^{(2)}{w}_j}\right)\\ -F{\displaystyle \sum _{j=1}^N{G}_{ij}^{(2)}{w}_j}=I_0\ddot{w}+\left(c_1{I}_4-{c_1}^2{I}_6\right){\displaystyle \sum _{j=1}^N{G}_{ij}^{(1)}{\ddot{\psi }}_j}-{c_1}^2{I}_6{\displaystyle \sum _{j=1}^N{G}_{ij}^{(2)}{\ddot{w}}_j}\end{array}$$

(41)

$$\begin{array}{l}\left(D_{11}-2c_1{H}_{11}+{c_1}^2{J}_{11}\right){\displaystyle \sum _{j=1}^N{G}_{ij}^{(2)}{\psi }_j}+\left({c_1}^2{J}_{11}-c_1{H}_{11}\right){\displaystyle \sum _{j=1}^N{G}_{ij}^{(3)}{w}_j}\\ +\left(2c_2{B}_{55}-A_{55}-{c_2}^2{D}_{55}\right)\left({\psi }_i+{\displaystyle \sum _{j=1}^N{G}_{ij}^{(1)}{w}_j}\right)=\left(I_2-2c_1{I}_4+{c_1}^2{I}_6\right){\ddot{\psi }}_i+\left({c_1}^2{I}_6-c_1{I}_4\right){\displaystyle \sum _{j=1}^N{G}_{ij}^{(1)}{\ddot{w}}_j}\end{array}$$

(42)

for the TSDBT, and

$$\begin{array}{l}{E}_{\mathrm{s}}{h}^3\left(C_1{\displaystyle \sum _{j=1}^N{G}_{ij}^{(4)}}{w}_j-C_2{\displaystyle \sum _{j=1}^N{G}_{ij}^{(3)}}{\varphi }_{1j}-C_3{\displaystyle \sum _{j=1}^N{G}_{ij}^{(3)}}{\varphi }_{2j}\right)+F{\displaystyle \sum _{j=1}^N{G}_{ij}^{(2)}}{w}_j\\ ={\rho }_{\mathrm{s}}{h}^3\left(B_1{\displaystyle \sum _{j=1}^N{G}_{ij}^{(2)}}{\ddot{w}}_j-\frac{I_0}{h^3}{\ddot{w}}_i-B_2{\displaystyle \sum _{j=1}^N{G}_{ij}^{(1)}}{\ddot{\varphi }}_{1j}-B_3{\displaystyle \sum _{j=1}^N{G}_{ij}^{(1)}}{\ddot{\varphi }}_{2j}\right)\end{array}$$

(43)

$$\begin{array}{l}{E}_{\mathrm{s}}\left(C_2{\displaystyle \sum _{j=1}^N{G}_{ij}^{(3)}}{w}_j-C_4{\displaystyle \sum _{j=1}^N{{\rm K}}_{ij}^{(2)}}{\varphi }_{1j}-C_5{\displaystyle \sum _{j=1}^N{G}_{ij}^{(2)}}{\varphi }_{2j}\right)+\frac{E_{\mathrm{s}}}{2h^2\left(1+v\right)}\left(C_7{\varphi }_{1i}+C_8{\varphi }_{2i}\right)\\ ={\rho }_{\mathrm{s}}\left(B_2{\displaystyle \sum _{j=1}^N{G}_{ij}^{(1)}}{\ddot{w}}_j-B_4{\ddot{\varphi }}_{1i}-B_5{\ddot{\varphi }}_{2i}\right)\end{array}$$

(44)

$$\begin{array}{l}{E}_{\mathrm{s}}\left(C_3{\displaystyle \sum _{j=1}^N{G}_{ij}^{(3)}}{w}_j-C_5{\displaystyle \sum _{j=1}^N{G}_{ij}^{(2)}}{\varphi }_{1j}-C_6{\displaystyle \sum _{j=1}^N{G}_{ij}^{(2)}}{\varphi }_{2j}\right)+\frac{E_{\mathrm{s}}}{2h^2\left(1+v\right)}\left(C_8{\varphi }_{1i}+C_9{\varphi }_{2i}\right)\\ ={\rho }_{\mathrm{s}}\left(B_3{\displaystyle \sum _{j=1}^N{G}_{ij}^{(1)}}{\ddot{w}}_j-B_5{\ddot{\varphi }}_{1i}-B_6{\ddot{\varphi }}_{2i}\right)\end{array}$$

(45)

for the BSDBT. The boundary conditions in Table 2 can be handled in the same way, and their discretised forms are given in

Table 3. Discretised boundary conditions for different beam theories.

Beam Theory	End Support	Boundary Conditions
CBT	Hinged	$w_1=0$$,{\displaystyle \sum _{j=1}^N{G}_{2j}^{(2)}}{w}_j=0$ at x = 0 $w_N=0$$,{\displaystyle \sum _{j=1}^N{G}_{(N-1)j}^{(2)}}{w}_j=0$ at x = L
	Clamped	$w_1=0$$,{\displaystyle \sum _{j=1}^N{G}_{2j}^{(1)}}{w}_j=0$ at x = 0 $w_N=0$$,{\displaystyle \sum _{j=1}^N{G}_{(N-1)j}^{(1)}}{w}_j=0$ at x = L
FSDBT	Hinged	$w_1=0$$,{\displaystyle \sum _{j=1}^N{G}_{1j}^{(1)}}{\psi }_j=0$ at x = 0 $w_N=0$$,{\displaystyle \sum _{j=1}^N{G}_{Nj}^{(1)}}{\psi }_j=0$ at x = L
	Clamped	$w_1=0$$,{\psi }_1=0$ at x = 0 $w_N=0$$,{\psi }_N=0$ at x = L
TSDBT	Hinged	$w_1=0$$,{\displaystyle \sum _{j=1}^N{G}_{2j}^{(2)}}{w}_j=0$$,{\displaystyle \sum _{j=1}^N{G}_{1j}^{(1)}}{\psi }_j=0$ at x = 0 $w_N=0$$,{\displaystyle \sum _{j=1}^N{G}_{(N-1)j}^{(2)}}{w}_j=0$$,{\displaystyle \sum _{j=1}^N{G}_{Nj}^{(1)}}{\psi }_j=0$ at x = L
	Clamped	$w_1=0$$,{\psi }_1=0$$,{\displaystyle \sum _{j=1}^N{G}_{1j}^{(1)}}{w}_j=0$ at x = 0 $w_N=0$$,{\psi }_N=0$$,{\displaystyle \sum _{j=1}^N{G}_{Nj}^{(1)}}{w}_j=0$ at x = L
BSDBT	Hinged	$w_1=0$$,{\displaystyle \sum _{j=1}^N{G}_{2j}^{(2)}}{w}_j=0$$,{\displaystyle \sum _{j=1}^N{G}_{1j}^{(1)}}{\varphi }_{kj}=0$ at x = 0 $w_N=0$$,{\displaystyle \sum _{j=1}^N{G}_{(N-1)j}^{(2)}}{w}_j=0$$,{\displaystyle \sum _{j=1}^N{G}_{Nj}^{(1)}}{\varphi }_{kj}=0$ at x = L (k = 1, 2)
	Clamped	$w_1=0$$,{\displaystyle \sum _{j=1}^N{G}_{2j}^{(1)}}{w}_j=0$$,{\varphi }_{k1}=0$ at x = 0 $w_N=0$$,{\displaystyle \sum _{j=1}^N{G}_{(N-1)j}^{(1)}}{w}_j=0$$,{\varphi }_{kN}=0$ at x = L (k = 1, 2)

.

Substitution of the boundary conditions in Table 3 into the discretised governing equations of each theory leads to an algebraic system that governs the buckling and free vibration of the beam and can be written in a matrix form as

$$\mathbf{M}\ddot{\mathbf{d}}+\left(\mathbf{K}-F{\mathbf{K}}_{\mathrm{F}}\right)\mathbf{d}=\mathbf{0}$$

(46)

in which M is the mass matrix; K and K_F are the stiffness matrix and geometric stiffness matrix, respectively; and the unknown displacement vector d is composed of w_j, ψ_j, or ϕ_kj (k = 1, 2; j = 1, 2, …, N).

For the buckling analysis, the inertia terms are omitted and Equation (46) reduces to an eigenvalue equation as

$$\left(\mathbf{K}-F{\mathbf{K}}_{\mathrm{F}}\right)\mathbf{d}=\mathbf{0}$$

(47)

from which the critical buckling load can be obtained as the lowest positive eigenvalue.

For the free vibration analysis, the axial load F is absent and the displacement vector d is expressed as d = d*sin(Ωt). Equation (46) is then rewritten as

$$\left(\mathbf{K}-{\Omega }^2\mathbf{M}\right)\mathbf{d}=\mathbf{0}$$

(48)

where Ω is the natural frequency and can be determined through a standard eigenvalue algorithm.

4. Finite Element Simulation

FE simulation is a well-accepted tool for structural analysis, and the simulation results are used in this paper for direct comparison with those of different beam theories. To begin with, the FE model of the FG porous beam is created in the commercial software ABAQUS 2016 and meshed with eight-node reduced integration hexahedral elements (C3D8R) of controlled size, as shown in Figure 3. Due to the graded distribution of porosity, the material properties of the FG porous beam change continuously along the thickness direction, which is difficult to implement in the software. As an alternative, the continuous variation in material properties is approximated in a stepwise manner. To this end, the FG porous beam model is divided into a number of layers of equal thickness via the “offset faces” option, and the material properties are constant in each individual layer but change from layer to layer across the thickness. For the buckling analysis, an axial compressive load is applied at the reference point (RP) that couples with the right end of the beam. The step procedure named “Linear perturbation/Buckle” is selected to obtain the critical buckling load. For the free vibration analysis, the “Linear perturbation/Frequency” procedure is used to determine the natural frequency. Various numbers of layers were examined and it is established that a total number of 40 layers will suffice in simulating the FGP beam with continuous material distributions and furnishing converged simulation results.

5. Results and Discussion

In this section, we examine the validity and accuracy of the four beam theories in predicting the critical buckling load and fundamental frequency of FG porous beams by comparing them with their counterparts obtained from FE simulations. It is assumed that the porous beams of thickness h = 0.1 m are made from steel with material properties as follows: E_s = 200 GPa, vs. = 0.3, and ρ_s = 7850 kg/m³. Unless otherwise stated, the critical buckling load and fundamental frequency results are given in the dimensionless form as

$$f_{\mathrm{cr}}=\frac{100F_{\mathrm{cr}}}{h{E}_{\mathrm{s}}},\omega =10\Omega h\sqrt{\frac{{\rho }_{\mathrm{s}}}{E_{\mathrm{s}}}}$$

(49)

5.1. Convergence and Validation

Convergence analysis is first conducted, and the dimensionless buckling and free vibration results of the FG porous beams with different numbers of grid points are tabulated and compared in

Table 4. Dimensionless critical buckling loads and fundamental frequencies of the FGX porous beam varying with the total number of grid points (C-C, L/h = 10, e₀ = 0.99).

N	fcr				ω
	CBT	FSDBT	TSDBT	BSDBT	CBT	FSDBT	TSDBT	BSDBT
7	2.2697	1.7505	1.5349	1.0178	0.7041	0.6298	0.6242	0.1370
11	2.1117	1.7922	1.5666	0.0791	0.6982	0.6325	0.5932	0.4348
15	2.1116	1.7922	1.5666	0.2115	0.6982	0.6325	0.5886	0.5042
19	2.1116	1.7922	1.5666	0.4387	0.6982	0.6325	0.5873	0.5115
23	2.1116	1.7922	1.5666	0.7764	0.6982	0.6325	0.5870	0.5127
27	2.1116	1.7922	1.5666	1.2093	0.6982	0.6325	0.5869	0.5130
31	2.1116	1.7922	1.5666	1.2093	0.6982	0.6325	0.5869	0.5130

. It is seen that the present solutions are convergent for all beam theories when the total number of grid points is increased to N = 27. Hence, N = 27 is used in the subsequent calculations. Following that, validation analysis is performed, and the dimensionless critical buckling loads and fundamental frequencies of FG porous beams with different porosity gradients are calculated and compared in

Table 5. Comparison of dimensionless critical buckling loads and fundamental frequencies for the FG porous beams with different porosity distributions (C-C, L/h = 20, e₀ = 0.8).

Source	${\mathit{f}}_{\mathbf{cr}}=2(1-{\mathit{\nu }}^2){\mathit{F}}_{\mathbf{cr}}/\mathit{h}{\mathit{E}}_{\mathbf{s}}$			$\mathit{\omega }=\mathit{\Omega }\mathit{L}\sqrt{{\mathit{\rho }}_{\mathbf{s}}(1-{\mathit{\nu }}^2)/{\mathit{E}}_{\mathbf{s}}}$
	UD	FGX	FGO	UD	FGX	FGO
CBT	69.091	116.887	48.213	0.260	0.330	0.212
FSDBT	67.318	112.593	47.466	0.256	0.322	0.210
TSDBT	67.318	111.014	47.596	0.255	0.325	0.215
BSDBT	67.320	110.500	47.570	0.256	0.319	0.210
FE	67.557	110.808	47.527	0.256	0.327	0.216
Ref. [24]	67.103	110.322	47.520	0.255	0.326	0.215

with the FE simulation results, as well as those in the literature [24]. As can be observed, the present results agree well with the simulation and existing ones.

Nevertheless, Table 4 and Table 5 show that there exist differences between the results of the four beam theories, and the differences become much more significant for FGX porous beams with a high porosity (e₀ = 0.99). To understand this, the difference and accuracy of those beam theories in predicting the critical buckling load and fundamental frequency of FG porous beams are comprehensively examined through various numerical examples in the following sections.

5.2. Buckling Analysis

We first compare the buckling results predicted by the theories with those of FE simulation in

Table 6. Comparison of dimensionless critical buckling loads obtained from different beam theories and FE simulations (L/h = 10).

BCs	e0	FGX					UD					FGO
		CBT	FSDBT	TSDBT	BSDBT	FE	CBT	FSDBT	TSDBT	BSDBT	FE	CBT	FSDBT	TSDBT	BSDBT	FE
H-H	0.00	0.8225	0.8019	0.8019	0.8019	0.8020	0.8225	0.8019	0.8019	0.8019	0.8020	0.8225	0.8019	0.8019	0.8019	0.8020
	0.90	0.5547	0.5331	0.5210	0.5129	0.5119	0.3128	0.3049	0.3049	0.3049	0.3050	0.1684	0.1663	0.1667	0.1667	0.1655
	0.99	0.5279	0.5054	0.4856	0.4430	0.4097	0.2410	0.2349	0.2349	0.2349	0.2350	0.1030	0.1021	0.1022	0.1022	0.1024
H-C	0.00	1.6826	1.5908	1.5909	1.5909	1.5954	1.6826	1.5908	1.5909	1.5909	1.5954	1.6826	1.5908	1.5909	1.5909	1.5954
	0.90	1.1348	1.0400	0.9910	0.9600	0.9585	0.6398	0.6049	0.6049	0.6049	0.6067	0.3445	0.3352	0.3370	0.3366	0.3354
	0.99	1.0800	0.9816	0.9033	0.7584	0.6572	0.4929	0.4660	0.4661	0.4661	0.4674	0.2107	0.2066	0.2073	0.2072	0.2084
C-C	0.00	3.2900	2.9836	2.9840	2.9842	3.0004	3.2900	2.9836	2.9840	2.9842	3.0004	3.2900	2.9836	2.9840	2.9842	3.0004
	0.90	2.2188	1.9092	1.7637	1.6772	1.6738	1.2510	1.1345	1.1347	1.1347	1.1409	0.6736	0.6419	0.6478	0.6467	0.6464
	0.99	2.1116	1.7922	1.5666	1.2093	0.9847	0.9638	0.8741	0.8742	0.8742	0.8790	0.4119	0.3981	0.4005	0.4002	0.4040

for the FG porous beams with different boundary conditions and porosity coefficients. e₀ = 0 represents that the beam is homogeneous without pores. As expected, the porosity gradient FGX with fewer pores near the surfaces provides the highest stiffness and possesses the greatest critical buckling loads, followed by UD and FGO. Compared to other theories, the CBT always gives the highest predictions due to the fact that the transverse shear deformation is neglected in this theory. In contrast, the FSDBT, TSDBT, and BSDBT furnish very close critical buckling loads to the simulation results for the FGO and UD porous beams, regardless of the porosity value and boundary conditions. Thus, the following discussion centres on the FGX porous beams. From the results of FGX porous beams in Table 6, it is clear that the BSDBT predictions are closer to the simulation ones than other theories, and the difference between those theories and simulation results changes with respect to the porosity coefficient and boundary conditions.

For a better visualisation, the variations in the critical buckling load with respect to the porosity coefficient are plotted in Figure 4 for the FGX porous beams with different boundary conditions. The figure shows that the critical buckling load decreases but the difference between the results of theories and simulations rises as the porosity coefficient increases, and these effects become much more pronounced when the porosity coefficient exceeds 0.9. Nonetheless, the BSDBT predictions are always the closest to the simulation results, followed by those of the TSDBT, FSDBT, and CBT in order.

To evaluate the accuracy of those theories at higher porosity values (e₀ = 0.9, 0.98, and 0.99), the difference between the results obtained from theories and simulation are calculated by

$$difference=\frac{\left|f_{cr}(\mathrm{FE})-f_{cr}(\mathrm{Theory})\right|}{f_{cr}(\mathrm{FE})}\times 100\%$$

(50)

and compared in Figure 5, where results for e₀ = 0 are also provided as a reference. It is concluded from the figure that for the homogeneous beam (e₀ = 0), all the shear deformation beam theories provide very close predictions to the simulation results. However, for the high porosity values (e₀ ≥ 0.9), the difference between the results of theories and simulations significantly increases, and this effect turns out to be less noticeable for the soft (hinged) end supports. It should be noted that when the porosity coefficient e₀ is increased to 0.98, only the BSDBT gives close predictions to the simulation results, with a difference less than 10%. Unfortunately, when e₀ = 0.99, the differences go up to 22.8% and 15.4% for the C-C and H-C end supports, respectively.

Figure 6 displays the variation in critical buckling load with the slenderness ratio for the FGX porous beams with different boundary conditions. The critical buckling load and the difference between those beam theories declines as the slenderness ratio increases because a higher slenderness ratio leads to a lower stiffness and a weaker transverse shear effect. The differences between the results of various theories and the simulation at L/h = 6, 10, and 20 are further plotted and compared in Figure 7. It is seen that only the BSDBT provides very close predictions to the simulation ones with differences less than 10% when the slenderness ratio is increased to 10. The differences are further reduced to 3.02%, 1.84%, and 0.98% for the C-C, H-C, and H-H end supports in order when L/h = 20.

5.3. Free Vibration Analysis

We next investigate the accuracy and difference of those beam theories in predicting the fundamental frequencies of FG porous beams.

Table 7. Comparison of dimensionless fundamental frequencies obtained from different beam theories and FE simulations (L/h = 10).

BCs	e0	FGX					UD					FGO
		CBT	FSDBT	TSDBT	BSDBT	FE	CBT	FSDBT	TSDBT	BSDBT	FE	CBT	FSDBT	TSDBT	BSDBT	FE
H-H	0.00	0.2838	0.2802	0.2802	0.2801	0.2803	0.2837	0.2802	0.2802	0.2801	0.2803	0.2837	0.2802	0.2802	0.2801	0.2803
	0.90	0.2964	0.2907	0.2874	0.2850	0.2849	0.2228	0.2201	0.2201	0.2199	0.2201	0.1637	0.1627	0.1629	0.1628	0.1626
	0.99	0.3084	0.3019	0.2961	0.2829	0.2724	0.2088	0.2062	0.2062	0.2061	0.2062	0.1367	0.1361	0.1362	0.1362	0.1364
H-C	0.00	0.4429	0.4283	0.4285	0.4279	0.4294	0.4429	0.4283	0.4285	0.4279	0.4294	0.4429	0.4283	0.4285	0.4279	0.4294
	0.90	0.4626	0.4391	0.4273	0.4194	0.4188	0.3479	0.3363	0.3365	0.3361	0.3372	0.2556	0.2514	0.2522	0.2518	0.2523
	0.99	0.4814	0.4548	0.4338	0.3944	0.3608	0.3259	0.3151	0.3152	0.3149	0.3159	0.2135	0.2110	0.2115	0.2111	0.2123
C-C	0.00	0.6426	0.6054	0.6061	0.6061	0.6083	0.6426	0.6054	0.6061	0.6061	0.6083	0.6426	0.6054	0.6061	0.6061	0.6083
	0.90	0.6710	0.6127	0.5864	0.5710	0.5674	0.5046	0.4754	0.4760	0.4759	0.4777	0.3708	0.3600	0.3621	0.3618	0.3629
	0.99	0.6982	0.6325	0.5869	0.5131	0.4487	0.4728	0.4454	0.4459	0.4459	0.4475	0.3098	0.3033	0.3044	0.3044	0.3064

compares the fundamental frequency results of different beam theories and FE simulations for FG porous beams with different boundary conditions and porosity coefficients. Similar to the observations in buckling analysis, the FGX porous beam has the highest fundamental frequency and all shear deformation theories give very close predictions to simulation results for UD and FGO porous beams. Again, there exist certain differences between the results of theories and simulations for FGX porous beams, and the differences turn out to be extremely large at e₀ = 0.99.

To further illustrate the influence of the porosity coefficient, Figure 8 depicts the fundamental frequency versus porosity coefficient curves of FGX porous beams with different boundary conditions. According to the CBT, FSDBT, and TSDBT curves, the fundamental frequency slightly increases as the porosity coefficient rises. However, both the BSDBT and simulation results show that the frequency decreases with the porosity coefficient increasing, and this effect becomes much more noticeable when the porosity value e₀ increases from 0.9 to 0.99. The figure also shows that all the shear deformation theories give very close predictions to simulation results when e₀ < 0.9, but there are distinct differences between the results obtained from the theories and simulation when e₀ > 0.9. The differences at high porosity values (e₀ = 0.9, 0.98, and 0.99) are calculated and compared in Figure 9. Again, the difference is observably elevated when the porosity coefficient is increased to 0.99 and the end supports become more rigid. Nonetheless, the BSDBT always provides the closest predictions to simulation results, with the maximum differences of 6.90% and 14.35% at e₀ = 0.98 and 0.99, respectively.

Figure 10 shows the variation in the fundamental frequency with the slenderness ratio for the FGX porous beams with different boundary conditions. As observed in the buckling analysis, a higher slenderness ratio leads to a lower stiffness and thus results in smaller fundamental frequency and difference between the theories and simulation. For a better illustration, the differences between the results of various theories and simulation at L/h = 6, 10, and 20 are computed and compared in Figure 11. It is seen that for a thick FGX porous beam (L/h = 6), only the BSDBT furnishes reasonable predictions with the differences of 13.9%, 8.8%, and 3.32% for the C-C, H-C, and H-H end supports, respectively. As the slenderness ratio is increased to 20, all the shear deformation theories give close predictions with the maximum difference of 8.60%. It is also found that for the same beam, the difference of the frequency results is smaller than that of buckling results.

5.4. Analysis of Reasons for the Differences

The discrepancies between various theoretical predictions and simulation results can be attributed to different deformation hypotheses adopted in the beam theories, as mentioned in Section 3.1. The distinct displacement functions result in different normal and shear strain distributions in the cross-section of the beam, and they consequently give rise to the deviation of predictions. Figure 12 and Figure 13 display the normal and shear strain curves, respectively, of FGX porous beams with different porosity coefficients. It is seen that the differences between the strain curves of various theories and simulations become greater as the porosity increases. Particularly, when the porosity is increased to 0.99, the moduli at the top and bottom surfaces of the FGX beam are 100 times those at the midplane. In this case, the FGX beam will behave like a three-layered beam, and the beam theories cannot correctly determine the shear strain distribution in the cross-section. This is why all the theories lose their accuracy of prediction at e₀ = 0.99. Nonetheless, the strain curves of BSDBT are always closest to the simulation ones, which explains why the BSDBT solutions are invariably the closest to simulation ones.

6. Conclusions

In this paper, we have examined the validity and accuracy of various beam theories (CBT, FSDBT, BSDBT, and BSDBT) in predicting the buckling loads and fundamental frequencies of FG porous beams on the basis of the FE simulation results. Tabular and graphical results are presented to evaluate the influences of the porosity distribution and coefficient, slenderness ratio, and boundary conditions on differences between various theoretical predictions and simulation results. It is found that all beam theories give very close predictions for FGO and UD porous beams but show distinct deviations for FGX porous beams with high porosity. The differences between various theoretical and simulation results significantly increase as the porosity increases beyond 0.9 and the end supports become more rigid, but they decline as the slenderness ratio increases. The underlying reasons for the discrepancies observed between the theoretical predictions and simulation results have been analysed. Among those beam theories, the BSDBT always furnishes the closest predictions to simulation results because it is able to more adequately capture the deformations of the cross-section of FGX porous beams. Nonetheless, the BSDBT will also lose its accuracy when the FGX beam is fully clamped and the porosity is increased to 0.99. Those research findings are instrumental in developing enhanced theoretical models for the mechanical analysis of FGP structures.