To focus our attention on the FG beam, a straight beam structure with circular cross-section is taken into consideration, as shown in Figure 1. For convenience, the origin o is chosen as the central point of the circular cross-section on the left end, the x-axis is the centroid line of cross-sections, and the Cartesian coordinate system (x, y, z) is established. Meanwhile, the cylindrical coordinate system (x, r, θ) (see Figure 1) is introduced as well, which have the same origin o and the same x-axis described before. R is the radius of the circular cross-section. So, the coordinate transformation relationships are z = r ⁡ cos ⁡ θ , y = r ⁡ sin ⁡ θ , r = z 2 + y 2 (1)

To exclude the coupling of bending and torsion it will be assumed that the applied loads are symmetric about the x-z plane.

It is assumed that the material properties of FG beam vary along the radial direction. For simplicity, the power-law variation of Young’s modulus can be expressed as E ( r ) = ( E 2 − E 1 ) ( r R ) k + E 1 (2) where E ₁ and E ₂ denote the Young’s modulus of the metal (material-1) and ceramic (material-2) constituents of the FG beam respectively, k is the non-negative variable parameter (power-law exponent). Poisson’s ratio is assumed to be a constant, and this can be found in the relevant references (Pradhan and Chakraverty, 2014; Pei et al., 2019). Thus, the Young’s modulus of metal with respect to r/R is depicted for various k, as shown in Figure 2.

It can be seen from Figure 2 that k = 0 represents the homogeneous material, k = 1 represents the linearly distributed material and k = 5 represents the exponential distribution material which changes rapidly in r-direction.

Firstly, in order to facilitate theoretical derivation, the basic assumptions and definitions of beam problem are introduced.

Assumption 1: stress assumption.

Based on the elasticity theory for bending of a column, it is reasonable to assume that the following stress components vanish (Huang et al., 2013) σ z = σ y = τ z y = 0 (3)

The Assumption 1 adopted for columns is totally different from the plane stress assumption (Aydogdu, 2009; Şimşek, 2010; Pei et al., 2018) for beams with rectangular cross-section.

Assumption 2: displacement field assumption.

Based on Assumption 1, the displacement field of the column can be assumed in a general form as { u x ( x , y , z ) = f c ( x ) + ( z − z c ) ⋅ f 1 ( x ) + g ( y , z ) ⋅ f R ( x ) u y ( x , y , z ) = u y ( y , z ) u z ( x , y , z ) = u z ( x , 0,0 ) (4) where f _c (x), f ₁ (x) and f _R (x) are three unknown functions to be determined, g (y,z) is called the bidirectional warping function in (Geng et al., 2017) with two properties { g ( y , z ) | y = 0 , z = 0 = 0 ∂ g ( y , z ) / ∂ z | y = 0 , z = 0 = 0 (5)

It should be pointed out that, the displacement field assumption Eq. 4 in the present work is different from the displacement field in the theory for beams with rectangular cross-section (Geng et al., 2017; Pei et al., 2019) due to the bidirectional warping function g (y,z).

In the first of Eq. 4, z _c is the neutral point defined in Definition 1.

Definition 1: the definition of neutral point. (Pei et al., 2018) z c = ( 1 / B 0 ) ∫ A E ( r ) ⋅ z d A = 0 (6) where A is the area of cross-section and B 0 = ∫ A E ( r ) d A ≠ 0 (7) is the tensile rigidity. Unlike homogeneous material beams, the neutral axis of FG beam is often inconsistent with the centroid axis, except for some special distribution of material, such as the radial FG distribution.

In addition, based on the displacement field Eq. 4, the definition of basic variables in the beam theory is also introduced.

Definition 2: the definition of generalized displacement.

Based on Eq. 4, the variables in beam theory are defined as { u ¯ ( x ) = ( 1 / B 0 ) ∫ A E ( r ) ⋅ u x ( x , y , z ) d A ϕ ( x ) = ( 1 / B 2 ) ∫ A E ( r ) ⋅ ( z − z c ) ⋅ u x ( x , y , z ) d A w ( x ) = ( 1 / B 0 ) ∫ A E ( r ) ⋅ u z ( x , y , z ) d A ≡ u z ( x , 0,0 ) (8) where u ¯ ( x ) is the average stretch, ϕ(x) is the rotation of cross-section and w(x) is the deflection of the FG beam. Without loss of generality, the neutral point z _c is involved in the definition of ϕ(x), and B 2 = ∫ A E ( r ) ⋅ ( z − z c ) 2 d A ≠ 0 (9) is the flexural rigidity.

In order to derive the explicit expression of axial displacement, the shear stress free condition on the outer lateral surface is introduced as well which requires γ x r = 0 o n S (10) where S is the outer boundary S : S ( y , z ) = z 2 + y 2 ⇒ S ( y , z ) | S = R 2 (11)

According to Eq. 4, the radial shear strain γ _xr which can be expressed by γ x r = n y γ x y + n z γ x z = n z ( f 1 + w ′ ) + [ n y ∂ g ( y , z ) ∂ y + n z ∂ g ( y , z ) ∂ z ] ⋅ f R (12) where n _y and n _z are the cosines of outward normal on boundary S which are expressed { n y = ∂ S ∂ y n z = ∂ S ∂ z (13)

Substituting Eq. 13 into Eq. 12, the shear stress free condition is re-expressed as { n z ( f 1 + w ′ ) + [ n y ∂ g ( y , z ) ∂ y + n z ∂ g ( y , z ) ∂ z ] ⋅ f R } | S = 0 (14)

In order to make Eq. 14 be satisfied, considering a simple case, the partial differential equation is given by 1 n z [ n y ∂ g ( y , z ) ∂ y + n z ∂ g ( y , z ) ∂ z ] = S ( y , z ) (15)

It is worth mentioning that the right term of Eq. 15 can also be selected as sinusoidal function, exponential function, logarithmic function and other forms with respect to S (y, z), which will derive a variety of HSD models.

On the one hand, by solving partial differential equation Eq. 15, the explicit expression of g (y,z) can be derived as g ( y , z ) = ( z 3 + z y 2 ) 3 (16) and the third-order HSD model is obtained.

On the other hand, Eq. 14 is further recast as ( f 1 + w ′ ) + R 2 ⋅ f R ≡ 0 (17)

Based on Definition 1 Eq. 8 and making use of Eq. 17, the axial displacement of the FG beam can be represented by generalized displacements and arranged in the following orthogonal form u x ( x , y , z ) = u ¯ ( x ) + z ⋅ ϕ + g ∼ ( y , z ) ⋅ ( ϕ + w ′ ) (18) with { g ∼ ( y , z ) = [ λ ⋅ z − g ( y , z ) ] / ( R 2 − λ ) λ = ( 1 / B 2 ) ∫ A E ( r ) ⋅ z ⋅ g ( y , z ) d A (19)

In addition, we have following three properties to guarantee decoupling of the three terms in Eq. 18 { ∫ A E ( r ) ⋅ 1 ⋅ z d A = 0 ∫ A E ( r ) ⋅ 1 ⋅ g ∼ ( y , z ) d A = 0 ∫ A E ( r ) ⋅ z ⋅ g ∼ ( y , z ) d A = 0 (20)

It will be seen later that the orthogonal decomposition form of the axial displacement Eq. 18 is conducive to the complete decoupling of stretching, bending and warping of the FG beam, which lays a foundation for the simplification of the theoretical expression in this paper.

According to Eq. 18 and the last two of Eq. 4, the strains of the FG beam are further cast as { ε x = u ¯ ′ + z ⋅ ϕ ′ + g ∼ ( y , z ) ⋅ ( ϕ ′ + w ″ ) γ x z = ( 1 + ∂ g ∼ ( y , z ) / ∂ z ) ⋅ ( ϕ + w ′ ) γ x y = ( ∂ g ∼ ( y , z ) / ∂ y ) ⋅ ( ϕ + w ′ ) (21)

Based on the first of Eq. 21, the normal strain of present LBT can be decomposed into 3 parts, i.e., the uniform stretching, the linear normal strain and the warping strain. The linear normal strain (zϕ’) is the same as that obtained by the first-order model or TBT, while the warping strain a high-order deformation determined by bidirectional warping function g (y,z).

And the stresses are expressed as { σ x = E ( r ) ⋅ [ u ¯ ′ + z ⋅ ϕ ′ + g ∼ ( y , z ) ⋅ ( ϕ ′ + w ″ ) ] τ x z = G ( r ) ⋅ ( 1 + ∂ g ∼ ( y , z ) / ∂ z ) ⋅ ( ϕ + w ′ ) τ x y = G ( r ) ⋅ ( ∂ g ∼ ( y , z ) / ∂ y ) ⋅ ( ϕ + w ′ ) (22)

Similarly, the normal stress of present LBT can be decomposed into 3 parts corresponding to the three strains, which are the uniform stretching, the first-order normal stress and the warping stress, respectively.

According to the expression of Eq. 21, the generalized strains of FG beam can be defined as { ε ¯ = u ¯ ′ ( x ) κ = ϕ ′ ( x ) Γ = ϕ + w ′ Γ ′ = ϕ ′ + w ″ (23) And the generalized stresses (also called stress resultants) for the FG beam problem are defined as { N = ∫ A σ x ⋅ 1 d A M = ∫ A σ x ⋅ ( z − z c ) d A P = ∫ A σ x ⋅ ( g ∼ ( y , z ) − g ∼ c ) d A Q ∼ = ∫ A τ x z ⋅ 1 d A R z = ∫ A τ x z ⋅ [ ∂ g ∼ ( y , z ) ∂ z ] d A R y = ∫ A τ x y ⋅ [ ∂ g ∼ ( y , z ) ∂ y ] d A (24) where N, M and Q ∼ are, respectively, the membrane force, the moment and the shear force well known in Engineering, P is called the higher-order moment, R _z and R _y are two higher-order shear forces in the present work which are different from the R _z in the theory for beams with rectangular cross-section (Geng et al., 2017; Pei et al., 2019). It should be noted that the three factors, i.e., 1, z and g ∼ ( y ,   z ) , are extracted from the three terms in Eq. 18, respectively.

Substituting Eq. 22 into Eq. 24 and making use of Eq. 20, the uncoupled constitutive relations of the FG beam expressed by generalized stresses and strains are { N = B 0 ⋅ ε ¯ M = B 2 ⋅ κ Q = K s B s ⋅ Γ P = η B 2 ⋅ Γ ′ (25) with { Q = Q ∼ + R z + R y K s = K z + ξ z + ξ y B s = ∫ A G ( r ) d A (26) and { K z = ( 1 B s ) ∫ A G ( r ) ⋅ [ 1 + ∂ g ∼ ( y , z ) ∂ z ] d A ξ z = ( 1 B s ) ∫ A G ( r ) ⋅ [ 1 + ∂ g ∼ ( y , z ) ∂ z ] ⋅ [ ∂ g ∼ ( y , z ) / ∂ z ] d A ξ y = ( 1 B s ) ∫ A G ( r ) ⋅ [ ∂ g ∼ ( y , z ) ∂ y ] 2 d A η = ( 1 B 2 ) ∫ A E ( r ) ⋅ [ g ∼ ( y , z ) ] 2 d A (27) where B _s is the shear rigidity, K _z , η, ξ _z and ξ _y are the four non-dimensional rigidity coefficients, and K _s is like the shear coefficient K _T in TBT. It is interesting to see that, the four of Eq. 25 are respectively related to stretching, bending, shearing and higher-order bending which are in the decoupling form (Pei et al., 2018). Besides, these four coefficients will vary with the shape of the cross-section and the power-law exponent k.

By the way, the higher-order non-dimensional rigidity coefficients η, ξ _z and ξ _y are relatively small which have been proved in the theory for FG beams with rectangular cross-section (Pei et al., 2018; Pei et al., 2019), and this property still need to be verified in the present work.

Since the explicit expression of g (y,z) have already been obtained in Eq. 16, the non-dimensional rigidity coefficients in Eq. 27 are derived in the following form { K z = R 2 R 2 − λ − 4 3 ( R 2 − λ ) B 2 B 0 ξ z = λ R 2 ( R 2 − λ ) 2 − 4 R 2 − 5 λ 3 ( R 2 − λ ) 2 B 2 B 0 ξ y = λ 3 ( R 2 − λ ) 2 B 2 B 0 η = λ 2 ( R 2 − λ ) 2 + − 6 λ ∫ 0 R r 5 H ( r ) d r + ∫ 0 R r 7 H ( r ) d r 9 ( R 2 − λ ) 2 ∫ 0 R r 3 H ( r ) d r (28) with the following integral ∫ 0 R r n H ( r ) d r = ( H 2 − H 1 ) R n + 1 k + n + 1 + H 1 R n + 1 n + 1 w i t h n = 1,2,3 … (29)

For comparative purpose, aluminum is chosen as the material-1 and zirconia is the material-2, of which the material properties are listed in Table 1.

The non-dimensional rigidity coefficients in Eq. 28 are calculated and listed in Table 2. It is very interesting to notice the fact that the three higher-order rigidity coefficients (i.e. ξ _x , ξ _y and η) in Table 3 are indeed small quantities for the FG beam with the third-order HSD mode, which are usually ignored in traditional LBT (Huang et al., 2013; Huang, 2020). As for the shear coefficients K _z and K _s , we know that K _z is the shear coefficients in traditional LBT (Huang et al., 2013; Huang, 2020) and K _s is the shear coefficient corrected by the algebraic sum of ξ _x and ξ _y . For the homogeneous beam (when k = 0) with circular cross-section, there is no difference between the shear coefficient in the present LBT and K _T = 6/7 ˜ 0.857 in TBT. It is worth mentioning that, according to the change of material distribution, the shear coefficients K _z and K _s in LBT can be adjusted automatically, but the shear coefficient K _T in TBT cannot.

In the theory of elasticity, the principle of virtual work for elastic body can be written as δ W i n t E l a s t i c = δ W e x t E l a s t i c o r ∫ V ( σ x δ ε x + τ x z δ γ x z + τ x y δ γ x y ) d V = ∫ L q ( x ) ⋅ δ w ( x ) d x (30) where δ W i n t E is the virtual deformation energy and δ W e x t E is the external virtual work, q(x) is the distributed load with respect to x-coordinate.

On the other hand, according to the generalized strains and stresses of the FG beam, the principle of virtual work for FG beams can be written as 0 = ∫ L [ N δ u ¯ ′ + M δ ϕ ′ + Q δ ( ϕ + w ′ ) + P δ ( ϕ ′ + w ″ ) ] d x − ∫ L q ( x ) ⋅ δ w ( x ) d x (31)

Noticing Eqs. 21–24 and making use of the orthogonal properties Eq. 20, it is easy to prove the equivalence of Eq. 30 and Eq.(31). That is to say that the principle of virtual work of FG beam is variationally consistent with the one of elastic body which will lead to the HBT.

However, considering the characteristic of the small parameter η in Table 2 and ignoring the virtual deformation energy caused by high-order bending moment P, Eq. 31 will degenerate into an approximate form which is 0 = δ W i n t B e a m − δ W e x t B e a m = ∫ L [ N δ u ¯ ′ + M δ ϕ ′ + Q δ ( ϕ + w ′ ) ] d x − ∫ L q ( x ) ⋅ δ w ( x ) d x (32)

Via integration by parts, Eq. 32 yields 0 = [ N δ u ¯ ] | 0 L + [ M δ ϕ ] | 0 L + [ Q δ w ] | 0 L − ∫ L N ′ δ u ¯ d x − ∫ L [ M ′ − Q ] ⋅ δ ϕ d x − ∫ L [ Q ′ + q ( x ) ] ⋅ δ w d x (33)

From Eq. 33, the equilibrium equations are obtained as { δ u ¯ : N ′ = 0 δ ϕ : M ′ − Q = 0 δ w : Q ′ + q ( x ) = 0 (34) And the corresponding BCs (e.g. at the end of x = x ₀) are { g i v e n N | x = x 0 o r u ¯ | x = x 0 g i v e n M | x = x 0 o r ϕ | x = x 0 g i v e n Q | x = x 0 o r w | x = x 0 (35)

It can be seen from Eq. 34 that, equilibrium equations are frame independent, and the stretching, bending, higher-order bending are reciprocally uncoupled in the present framework with the help of the definitions of neutral point and generalized displacements, i.e., Eq. 6 and Eq. 8. In contrast, due to the coupled constitutive relations, the equilibrium equations will inevitably be coupled between bending and tension for the higher-order theories in (Aydogdu and Taskin, 2007; Thai and Vo, 2012).

Making use of Eq. 23 and Eq. 25, in terms of the generalized displacements, the governing equations of the variationally approximated lower-order theory are { B 0 u ¯ ″ ( x ) = 0 B 2 ϕ ‴ + q ( x ) = 0 B 2 ϕ ″ − K s B s ( ϕ + w ′ ) = 0 (36)

As we can see, the first of Eq. 36 is the separately governing equation for stretch. Using the second of Eq. 36 to eliminate ϕ, in terms of w, the third of Eq. 36 yields the governing equation for deflection as B 2 w ( 4 ) = q ( x ) − B 2 K s B s q ″ ( x ) (37)

After manipulations, it is not difficult to obtain the rotation as ϕ = − w ′ − B 2 K s B s ( w ‴ + 1 K s B s q ′ ( x ) ) (38)

Eq. 35 provides three pairs of different boundary conditions for Eq. 34. Among them, the first pair (2 in total) is also individually responsible for stretch governed by the first of Eq. 36 -a second-order ODE, while the last two pairs (4 in total) are individually for bending governed by Eq. 37 -a fourth-order ODE. Different from the coupled governing equations in (Huang, 2020), the present theory has concise mathematical form and clear physical meaning. Noticing that Eq. 34 and Eq. 35 are derived from the approximate principle of virtual work Eq. 32, compared with the variationally consistent HBT (Pei et al., 2018), the present LBT will bring great convenience to the calculation due to the reduction of the differential order of the governing equation.

In this section, the static bending problem of the FG beam will be studied to demonstrated the validity and accuracy of the present LBT. Considering the Clamped-Clamped (C-C) supported FG beam subjected to uniformly distributed load (q(x) = q ₀), based on Eq. 35, the C-C BCs can be expressed by x = − L / 2 : { u ¯ = 0 ϕ = 0 w = 0 a n d x = L / 2 : { u ¯ = 0 ϕ = 0 w = 0 (39)

Since Eqs. 37, 38 are in the same form as the governing equations of lower-order theory for homogeneous beams (Duan and Li, 2016b), it is convenient to obtain the theoretical solutions expressed as { w = q 0 x 4 24 B 2 + d 3 x 3 + d 2 x 2 + d 1 x + d 0 ϕ = − q 0 x 3 6 B 2 − q 0 x K s B s − 3 d 3 x 2 − 2 d 2 x − d 1 − 6 d 3 B 2 / K s B s (40) with { d 0 = q 0 L 2 ( K s B s L 2 + 48 B 2 ) ( 384 K s B s B 2 ) d 1 = 0 d 2 = − q 0 ( K s B s L 2 + 24 B 2 ) ( 48 K s B s B 2 ) d 3 = 0 (41)

In order to illustrate the correctness of the calculation results in this paper, the comparison of dimensionless maximum transverse deflections w (0)*100E ₁ R ³/q ₀ L ³ by present LBT, by Huang (Huang, 2020) and by TBT with K _T = 6/7 are listed in Table 3.

In addition, for the case when L/R = 5 and k = 0 or 5, the dimensionless transverse deflections w(x)* 100E ₁ R ³/q ₀ L ³ are compared, as shown in Figure 3.

It can be seen from Figure 3 that the results are the same for a homogeneous beam (i.e., k = 0). For the FG beam when k = 5, the deflection obtained by the present LBT is slightly larger than that obtained by Huang, and is also larger than that obtained by TBT.

The fact is that the bidirectional warping effect on the cross-section can be described based on the present LBT which is not taken into consideration in TBT because of the rigid cross-section hypothesis and can only be described in z-direction for beams with rectangular cross-section because of the one-variate warping function g(z). In view of this, it is very interesting to study the warping effect on the circular cross-section through the bidirectional warping function g (y, z).

Considering the cross-section of the FG beam with L/R = 10 at x = 0, the normal strains are described by the present LBT for k = 0 and k = 5 are shown in Figure 4 and Figure 5, respectively.

It can be seen from Figure 4 and Figure 5 that the normal strains are both symmetric about the z-axis and anti-symmetric about the y-axis, and this characteristic can also be found in radial shear strain and stresses (or in Figure 8 and Figure 9). It can be found from Eq. 21 that, in the pure bending problem, the normal strains are composed of two parts, i.e., the linear part and the warping part. Because the warping part is much smaller than the linear part, the normal strains are mainly determined by the linear part. The warping part obeys the variation law of cubic function in z-direction, and is significantly affected by the power-law exponent k. By comparing Figure 4 and Figure 5, it can be seen that there is an obvious transition area of material properties on the cross-section when k = 5, which leads to a more obvious difference between tensile and compressive deformation in the warping part of the normal strain.

Meanwhile, the normal stresses are described by the present LBT for k = 0 and k = 5 are shown in Figure 6 and Figure 7, respectively.

It can be seen from Figure 6 and Figure 7 that the stress distribution is obviously affected by the power-law exponent k. Compared with the results for k = 0, when k = 5 the normal stress is obviously eased, and the transition regions of tensile and compressive stresses move to the upper and lower boundaries respectively. Meanwhile, the warping parts of normal stresses have little effect on the normal stress distribution, and the warping part is also eased for k = 5.

Since the shear stress free condition on the outer lateral surface is adopted in the HSD model, the radial shear strain and stress state on the cross-section need to be discussed.

Based on Eq. 12, the non-dimensional radial shear strain and stress can be expressed as { γ x r ⋅ R 2 − λ 6 R 3 ( ϕ + w ′ ) = ( z R ) [ 1 − ( r R ) 2 ] τ x r ⋅ R 2 − λ 6 G 1 R 3 ( ϕ + w ′ ) = ( z R ) [ 1 − ( r R ) 2 ] ⋅ [ ( G 2 G 1 − 1 ) ( r R ) k + 1 ] (42)

Obviously, according to the first of Eq. 42, the non-dimensional radial shear strain on cross-section of FG beam is independent of k, which is shown in Figure 8. The peaks of the radial shear strain appear near the centroid of the upper and lower semi-circles. In addition, the non-dimensional radial shear stresses on cross-section of FG beam for k = 0 and 5 are shown in Figure 9.

It can be seen from Figure 8 and Figure 9 that the radial shear strain and stresses are free on the boundary of the cross-section. This further verifies the rationality of the HSD model. In Figure 9, compared with radial shear stress for k = 0, the peaks of the radial shear stress move to the upper and lower boundaries for k = 5, respectively.