1 Introduction

The micro-scaled functionally graded (MFG) material is a new class of composite structures. Not the same as conventional composite materials, such as laminate and sandwich materials, the properties of MFG materials, such as the elasticity modulus, shear modulus, material density, and electrical conductivity, change continuously and smoothly, obeying a specific distribution function. Due to the smooth continuity of material properties, the stress concentrations as well as severe material failure may be exempted^[1-3]. Actually, the characteristics of MFG materials have been found in some structures in nature, such as pomelo peel^[4], bones^[5], and arapaima gigas scales^[6], and a better understanding of the statical and dynamical characteristics of MFG materials may be helpful to synthesize new materials applied for micro-nano-electro-mechanical systems (MEMS/NEMS), such as bio-inspired coatings for implants^[7], the micro-components in the form of shape memory alloy thin films^[8-9], nano-/micro-shaft^[10-12], atomic force micro-scopes^[13], and micro-electrically actuated devices including bent beam actuators^[14-16], and thermal buckling type actuators^[17-22].

There are many experimental phenomena showing that the statical and dynamical characteristics of micro-scaled continuous elements shall be highly size-dependent^[23-27]. These scaling phenomena can be predicted well by using a series of refined continuum mechanics models^[28-34]. The stiffening and softening effects of the micro-/nano-structures were studied^[35-40]. The effects of boundary conditions were also studied based on the refined continuum mechanics models^[38]. Since the modified couple stress theory was proposed by Yang et al.^[31], there have been many papers discussing the size effect of the MFG beams by utilizing this theory^[41-48]. A large number of papers discussing the scaling statical and dynamical characteristics of functionally graded (FG) materials have been published recently in terms of different non-classical continuum mechanics^{[12, 33, 34, 41, 49-69]}. Besides, Li et al.^[70] investigated the static and dynamic behaviors of axially FG micro-beams based on the nonlocal strain gradient theory. Li and Hu^[71-73] studied torsional behaviors of micro-bars based on nonlocal elasticities as well as nonlinear characteristics of FG beams based on the nonlocal strain gradient theory. These works focuses on thickness-directional or axially FG beams/plates.

It can be observed that these studies mentioned above are related to FG continuous components whose material properties vary in only one direction. However, in advanced machines, these properties may be expected to change in more than one direction, which means that conventional FG continuous components may not be useful in the design of such structures. For example, in modern aerospace shuttles and crafts, the temperature varies in two or three directions, i.e., on the outer surface along the top of the fuselage and through the thickness inside the plane^[74-75]. For this requirement, it is of great significance to develop bi-directional FG materials whose material properties vary multi-directionally. Based on the conventional continuum mechanics, several researches^[76-79] studied the thermoelastic behaviors of bi-directional FG materials. Nemat-Alla^[77] found that bi-directional FG materials enable more reduction of thermal stresses by choosing appropriate values of gradient parameters. Cho and Ha^[80] optimized the bi-directional distributions of heat resisting FG materials. Qian and Batra^[81] designed the bi-directional FG plate based on a higher-order shear and normal deformable plate theory for optimal natural frequencies. As the distribution of material properties is very essential for designing structures made of bi-directional FG materials, further studies of mechanical properties and behaviors of these components are of fundamental demands. Nie and Zhong^[82-83] studied the bending and vibration problems of bi-directional FG plates with material properties having an exponential-law variation in two directions. Lyu et al.^{[75, 84]} proposed a semi-analytical method for multi-directional FG beams and plates. Deng and Cheng^[85] studied the dynamical characteristics of a bi-directional FG Timoshenko beam. Static behaviors of bi-directional FG micro-beams were studied based on the modified couple stress theory^[86-87]. Nguyen et al.^[88] carried out a vibration analysis of bi-dimensional FG Timoshenko beams excited by a moving load. Șimșek^[89-90] studied the buckling, free and forced vibration behaviors of a bi-directional FG Timoshenko beam under various boundary conditions. The free vibration and buckling behaviors of bi-directional FG Euler-Bernoulli nano-beams based on the nonlocal theory were analyzed^[91-92]. Shafiei and Kazemi^[93] studied the buckling problems of bi-directional FG beams based on the Euler-Bernoulli beam theory. Among these engineering problems, the buckling problem, as one of the most significant structural problems, is an essential problem in the design of bi-directional FG micro-structures. As we know, buckling is a kind of instability leading to a failure mode. If a bi-directional MFG element is subject to a high compressive stress in its strong plane, the structure has a tendency to perform a sideway failure by buckling in its weaker plane. Besides, as important as its disadvantages, buckling behaviors have been used to design actuators^{[17-18, 94]}. Thus, we can also benefit from buckling behaviors of micro-structures.

As we know, the Euler-Bernoulli beam theory is the simplest one, which is applicable only for slender beams, since it neglects the effect of transverse shear strain. In order to take the effect of transverse shear deformation into consideration, the Timoshenko beam theory, or the first-order shear deformation beam theory, is applicable for the thick beam. The Timoshenko theory admits a non-zero transverse shear strain, but it needs to introduce an additional parameter (a shear correction factor) to correct the calculated results. In order to overcome this disadvantage, higher-order beam shear theories, such as the Reddy beam theory^[95], the Touratier beam theory^[96], the Karama beam theory^[97], the Soldatos beam theory^[98], the Aydogdu beam theory^[99], were proposed for the analysis of beams and plates. Higher-order shear beam theories are widely used to investigate the statical and dynamical behaviors of FG small-scaled beam structures (see Refs. [42]-[43] and [100]-[104] for further details on this topic). Besides, Pydah and Batra^[86] studied the static bending behaviors of circular bi-directional FG micro-beams by employing a logarithmic shear assumption.

Additionally, in our previous studies, we have discussed the buckling problems of uniform micro-beams^[105], conventional thickness-directional FG micro-beams^[72], axially FG micro- beams^[70], as well as torsional vibration problems of bi-directional FG micro-beams^[73]. But it still lacks of articles discussing the issue on the bi-directional FG beams considering higher-order shear deformations. These facts motivate this study, and the novel contributions are summarized as follows:

(ⅰ) The modified couple stress theory as well as various higher-order shear beam theories is used to formulate the mechanical model of through-thickness and through-length MFG beams.

(ⅱ) A generalized differential quadrature method (GDQM) is implemented for calculating the critical buckling loads of bi-directional MFG beams.

(ⅲ) The effects of through-thickness and through-length FG indexes, the length scale parameter, as well as different boundary conditions on the critical buckling loads are studied.

The works listed above will be organized as follows. In Section 2, we assume the Young's modulus of the bi-directional FG material to obey an exponential distribution function in both length and thickness directions. Section 3 simply reviews the modified couple stress theory and considers a displacement filed of a micro-structure-dependent beam theory that can include different kinds of higher-order shear assumptions as well as the two familiar beam theories (namely, Euler-Bernoulli and Timoshenko beam theories) as special cases. Besides, by considering the modified couple stress theory and the principle of minimum potential energy, the governing equations and boundary conditions for a micro-structure-dependent beam theory are derived to study the buckling behaviors of bi-directional FG beams. A numerical solution procedure, on the basis of the GDQM, is used to calculate the numerical results of the bi-directional FG beams in Section 4. Then in Section 5, the effects of two exponential FG indexes, the length scale parameter and different boundary conditions on the critical buckling loads are studied in detail. Finally, important conclusions are summarized.

2 Bi-directional FG material

A bi-directional FG material is applied for a micro-beam with the length L, as shown in Fig. 1. In previous works^{[51, 73, 106]}, the authors and co-workers discussed an FG beam with a circular cross-section. In this paper, we consider a constant rectangular cross-section, whose width and thickness are represented by b and h, respectively. The material properties E_LB, E_LT, E_RB, and E_RT represent Young's moduli at the left-bottom, left-top, right-bottom, and right-top points, respectively. The variation of Young's modulus at the arbitrary point in the beam body is assumed to be expressed as a product of two functions in length and thickness directions,

(1)

where f(x) and g(z) can be arbitrary functions. The effect of Poisson's ratio on the deformation of non-homogeneous FG beams has been proved to be much less than that of Young's modulus^[107], and then Poisson's ratio is assumed as a constant in the material body.

For a special case, Young's modulus of the bi-directional FG material is assumed to obey an exponential distribution function in both length and thickness directions^{[75, 83, 85, 108-109]}. This distribution law was also verified experimentally^[110-111] and was applied in previous studies^[112-113]. The FG index along the length direction is represented by α, and the FG index along the thickness is represented by β. Thus, the distribution functions f(x) and g(z) can be expressed as

(2)

With (2), Young's modulus of the material can be expressed as

(3)

Therefore, the relations of material properties at four corners can be obtained by substituting x=0 or L and z=-h/2 or h/2 into the above equation,

(4)

Interestingly, when α=0, the bi-directional FG beam can be reduced to a conventional through-thickness FG beam. By simply controlling β=0, a bi-directional FG beam can be reduced to an axially FG beam. When both α and β vanish, we obtain a homogeneous beam structure.

3 Size-dependent governing equation for bi-directional FG beams

In this section, a micro-beam model will be formed for the bi-directional FG beam based on the modified couple stress theory.

3.1 Modified couple stress theory

The modified couple stress theory, first presented by Yang et al.^[31], states that the strain energy function should be dependent no longer only on the strain tensor but also on the curvature tensor, that is,

(5)

where the operator ":" denotes the double dot production of tensors, i. e., σ : ε = σ_ijε_ij, ε represents the conventional strain tensor, σ represents the Cauchy stress tensor, χ represents the curvature tensor, and m represents the couple stress tensor. These tensors are defined as follows:

(6)

Here, Lamé's constants λ and G are functions of x and z, respectively, and they can be expressed as

where

and l represents the length scale parameter in the modified couple stress theory to capture size-dependent behaviors. θ=0.5curl (u) represents the rotation vector. The parameter l is a material constant, and its value can be determined by matching the result of molecular dynamics, lattice models or experiment tests. Askes and Aifantis^[114] studied the relations between refined continuum theories and lattice models and suggested that the value of length scale parameter can be fitted to the size of material micro-cells.

3.2 Beam models based on general shear theory

The general shear theory can account for all the longitudinal and transverse displacements and their inertia as well as the rotation and its inertia to investigate the buckling behaviors of the bi-directional FG beam. Let u₁, u₂, and u₃ be the displacements along the x-, y-, and z- directions, respectively. The displacements u₁, u₂, and u₃ for a general shear beam theory can take the form as follows^[43]:

(7)

in which u and w are the axial and transverse displacements on the geometric middle plane, respectively. ϕ is a function of x which represents the effect of shear strain on the total rotation of the cross-section during bending motions. F(z) is a function of z and represents the distribution of shear stress along the thickness direction of the beam. Several classical forms of the function F(z) based on the well-known shear deformation theories are shown in Table 1. Thus, the formulas of the non-zero strains of the beam can be given by

(8)

According to Hook's law by (6), we obtain the non-zero components of the stress tensor as

(9)

With the help of (6) and (7), the non-zero components of rotation vector and curvatures can be expressed as

(10)

(11)

Therefore, the non-zero couple stresses can be obtained as

(12)

Substituting the strain and curvature tensors into the strain energy of the modified couple stress theory (5), we can obtain the variation of the strain energy δU stored in the bi-directional FG extensible shear-deformable micro-beam as

(13)

where the stress resultants can be expressed as

(14)

with the integral parameter A_ij and B_ij defined as follows:

(15)

Let P, Q, M_c, and M_nc be the external axial load, the shear force, the classical moment, and the non-classical moment, respectively. The variation of potential energy δW can be expressed by the following procedure^[115]:

(16)

The modified couple stress theory has been used to develop formulations for the elastic potential energy as well as the kinetic energy of the bi-directional FG beam. Introduce the principle of minimum potential energy, i.e., δΠ=δ (δU-δW)=0, and use the integration by parts. The Lagrange equation can be obtained and expressed as

(17)

Because the coefficients of δu, δϕ, and δw are not always equal to zero, we have

(18)

For the boundary conditions at the two ends, we can obtain the following boundary conditions from the principle of minimum potential energy:

(19)

Substituting the stress resultants (14) into (18) yields the size-dependent governing equations in terms of the displacement field,

(20)

3.3 Timoshenko beam theory

For most cases, the higher-order shear theory beam models are too complicated to solve, though applying these theories can avoid identifying the shear correction factor. The first-order shear beam theory (the Timoshenko beam theory) is a most widely used theory for an analysis of the beam structure. The governing equations and the corresponding boundary conditions for the Timoshenko beam model are deduced in the same way as the previous section. When considering the Timoshenko beam theory, we obtain the function F(z)=z. Thus, the integral parameters A_ij and B_ij have the following relations:

With the above equations, the equilibrium equations (18) can be simplified to

(21)

and the boundary conditions can be reduced to

(22)

Here, these stress resultants are redefined as

Unlike higher-order shear beam theories, a shear correction factor κ_s needs to be introduced for the resultant Q using the Timoshenko beam theory. Substituting these stress resultants into (21), the size-dependent governing equations based on the Timoshenko beam theory can be obtained as

(23)

3.4 Euler-Bernoulli beam theory

When analyzing a slender beam structure, the effect of transverse shear deformation can be neglected. In this case, the Euler-Bernoulli beam theory as the simplest beam theory becomes applicable to the research of the static and dynamical characteristics for bi-directional FG beams. For an Euler-Bernoulli beam, the function F(z) vanishes. Thus, the integral parameters A₀₁, A₁₁, A₀₂, B₁₁, B₁₂, and B₂₂ equal zero. The equilibrium equations are reduced to

(24)

with the reduced boundary conditions

(25)

Here, the stress resultants are

For the Euler-Bernoulli beam theory, the governing equations (20) can be expressed as

(26)

4 Numerical solution method

A numerical solution procedure, on the basis of a GDQM^[116-117], will be used to calculate the numerical results of the bi-directional FG beams. The axial displacement u and rotation ϕ can be interpolated using the Lagrangian shape functions, and the deflection w is interpolated using the Hermite shape function (because the derivation of the deflection exists in the essential classical boundary condition). In this regard, the displacement field can be approximated using the series expansions as follows:

(27)

where ℓ_j and ħ_j represent Lagrangian interpolation functions and Hermite interpolation functions, respectively. The Lagrangian interpolation functions can be expressed as

The expressions of Hermite interpolation functions are listed in Appendix A. Here, N represents the number of sampling points, which in this study can use the well-known Gauss-Lobatto-Chebyshev points,

With these sampling points, the rth-order derivatives of displacements at the arbitrary sampling point x_i can be expressed as

(28)

For the sake of clarity, the expressions of two kinds of interpolation functions and their derivations are listed in Appendix A. For the sake of simplification, the displacement fields and their derivations at sampling points are expressed in a vector form as

(29)

(30)

where

4.1 GDQM for higher-order shear beam model

(20) can be further expressed as

(31)

Substituting (30) into (31), the matrix formulations (at domain points, i. e., i=2, 3, …, N-1) can be obtained as

(32)

where is a diagonal matrix.

(32) can be rewritten as

(33)

where

The statement of the boundary points (i. e., i=1 or N) should obey the boundary condition and can be expressed as

By inserting the constraints of these boundary conditions into (33), we can solve the critical bucking force.

4.2 GDQM for Timoshenko beam model

For the Timoshenko beam model, the components of matrices K and P of (33) are expressed as follows:

In the case of the Timoshenko beam model, the boundary points (i. e., i=1 or N) can be selected as the following form to obey the boundary condition:

4.3 GDQM for Euler-Bernoulli beam model

In the case of the Euler-Bernoulli beam model, the matrix formulations (at domain points, i. e., i=2, 3, …, N-1) can be obtained as

where

and the boundary points (i. e., i=1 or N) can be expressed as

5 Numerical results and discussion 5.1 Verification of present results

In this subsection, the dimensionless critical buckling loads PL²/(EI) for homogeneous micro-beams with two simply supported ends based on various beam theories are calculated and compared with the results in Ref. [43] in order to verify the accuracy of the present GDQM solutions.

In Table 2, the effect of Poisson's ratio is neglected (i. e., ). In this subsection, Young's modulus E=1.44 GPa, Poison's ratio ν=0.38, the thickness of the beam h=17.6 μm, the length L=20 h, and the length scale parameter l=17.6 μm for the modified couple stress theory and l=0 μm for the classical theory. Particularly, for Timoshenko beams, the shear factor κ_s=5/6 is assumed as discussed by Deng and Cheng^[85]. As expected, good agreement can be reached between the present solutions and the results in the reference. This clearly indicates that the present solution method (GDQM) is reliable for analyzing the bucking problems of micro-beams.

5.2 Critical buckling load of bi-directional FG micro-beams for different beam theories

This subsection is devoted to analysis of the critical buckling loads of a bi-directional FG micro-beam with two simply supported ends. The geometrical parameters of the cross-section are the same as the last subsection while the distribution of Young's modulus is assumed to obey the exponential law described in (3). Suppose that the reference Young's modulus at the left-bottom edge of the beam structure is considered as E_LB=1.44 GPa, exponential law indexes α=1 and β=1. The solutions with different values of the length L and the length scale parameter l obtained with the GDQM are presented in Table 3. The dimensionless critical buckling load can be increased distinctly with the increase in the length scale parameter l/h. Thus, the influence of size effect (the coupled stress effect) on the stiffness of micro-beams is significant. With the Euler-Bernoulli beam theory, the critical buckling loads of micro-beams are independent of the total length of the beam. For different higher-order shear theories, the relative errors of the numerical solutions are sufficiently small. For Timoshenko beams, the solutions increase sharply with the increasing of the length of micro-beams. The relative errors between the Timoshenko and higher-order shear beams can be increased with the increase in the length scale parameter l and with the decrease in the beam length L. With consideration of a larger range of l/h, Fig. 2 plots the results for an FG beam with the thickness h=17.6 μm, the aspect ratio L/h=5, and FG indexes α=β=1. The variations of errors between different shear theories do not change with the value of the dimensionless parameter l/h. When the length to height ratio L/h < 10, the effect of shear deformation on the critical buckling loads becomes evident. Thus, when the length of micro-beam is short enough, the higher-order shear theory must be considered and used to calculate the critical buckling loads. With the decrease in the volume of micro-beams (i. e., L as well as h is decreasing), the effect of couple stress becomes significant, because l, as a material length scale parameter, is a constant and becomes comparable with h or L.

5.3 Effect of exponential law indexes

As discussed above, to obtain an accurate solution, higher-order shear assumptions should be adopted, and the solutions for different shear distribution functions are close enough. Thus, for convenience, only a Reddy beam theory is considered in this section.

Figure 3 plots the distribution of Young's modulus, as expressed in (3), for different values of axially FG index α and thickness-directional FG index β. As expected, for the cases α=0 and β=0, the micro-beam structure is reduced to the classical thickness-directional FG beam and the axially FG beam, respectively. Figure 4 plots Young's modulus at the point (x₀, z₀) varying with the modulus E_LB and FG indexes. The modulus changes differently with varying indexes α and β at different points. Then, we take a look at the distribution of Young's modulus in bi-directional FG micro-beams as shown in Fig. 5. Figure 5(a) plots the dimensionless Young's modulus distribution at z = -h/2 along the x-direction with different values of exponential law index α. At the arbitrary point of the bottom surface of beam structure, Young's modulus can be increased with the increasing of α. Similar disciplinarian can be found in Fig. 5(b) that the exponential law index β positively increases Young's modulus in the thickness direction. In other words, the increasing exponential law indexes α and β can make the bi-directional FG material get stiffer. Note that this phenomenon only occurs when the varying function (3) is adopted. It is reported that, for other distributions, the beam will get soften with the increasing FG index^[118].

Table 4 shows the numerical solutions with various values of the exponential FG indexes α and β for l=17.6 μm. Particularly, the bi-directional FG material can be reduced to a homogeneous material if we consider α=0 and β=0. The increasing value of β stiffens the micro-beam structure a little more than the increasing value of α, and the effects of two FG indexes on the critical buckling loads are both obvious.

5.4 Effect of couple stresses

In this section, we investigate the size effect on the buckling behaviors. Figure 6(a) shows that the critical load varies with different sizes and aspect ratios for bi-directional FG Reddy beams. Here, two exponential law indexes are identical and equal to 1, the couple stress effect scale parameter l = 17.6 μm, and the width of the beam structure is considered as b = 2l = 35.2 μm. As expected, the critical load decreases with the increasing of the beam length for a specified value of thickness h. Figure 6(b) plots the critical load versus the thickness, with different values of aspect ratio. In this figure, we can obviously compare the solutions for the modified couple stress theory with the solutions neglecting the effect of couple stresses. The effects of couple stresses can significantly stiffen the bi-directional FG beam structures in micro-scale (especially when the thickness of the beam becomes competitive with the material couple stress parameter l). Li^[35] and Li et al.^[36-37] have also reported similar results for micro-/nano-structures. It is interesting that, for a specific aspect, there is a minimum value of the critical load ratio (at the point with h/l being about 1.6). For the case h/l > 1.6, the critical load increases with increasing h (or the volume of beam), because the bending stiffness is mainly affected by the thickness. While h/l < 1.6, the scale effect on the stiffness becomes significant. Thus, the critical load increases with the decreasing volume (h and L).

In Figs. 7(b)-7(d), the solid curves plot the results of thickness-directional FG beams (i. e., α=0). Specially, the blue solid line in Fig. 7(a) represents the critical buckling loads for homogeneous beam structures. As shown in Figs. 7(a)-7(e), the critical buckling loads become large with the increasing ratio of the length scale parameter to the thickness l/h for various indexes α and β. Besides, since the two indexes α and β can stiffen the beam structure, the increasing indexes can increase the critical buckling loads of bi-directional FG micro-beams.

Then, a simply supported FG Reddy beam with geometric parameters h=17.6 μm and L=20 h is considered. Young's modulus at the left-bottom edge of the beam is considered as E_LB=1.44 GPa. Figure 7 plots the numerical solutions of the critical buckling loads versus the dimensionless length scale parameter l/h for various values of exponential law indexes α and β.

5.5 Effect of boundary conditions

This section investigates the critical buckling load of bi-directional FG beams under different kinds of boundary conditions, i. e., simply supported-simply supported (S-S) and clamped-clamped (C-C) beams. Only the Reddy shear assumption is considered to plot the critical buckling loads. The geometric parameters and material properties are the same as the last subsection. Figure 8 shows the dimensionless critical buckling load PL²/(EI) varying with the dimensionless length scale parameter l/h for bi-directional FG Reddy beams (α = 1 and β = 1) and homogeneous beams (α = 0 and β = 0) under C-C and S-S boundary conditions. As shown in Fig. 8, the solutions of C-C and S-S beams increase with the dimensionless length scale parameter l/h and exponential FG indexes α and β. Because the stiffness of a C-C beam is higher than that of an S-S beam, the critical buckling loads of bi-directional FG beams under two clamped boundary conditions are larger than those of S-S bi-directional FG beams. What is more, the effect of the boundary conditions is much evident as well as the effect of l, α, and β.

6 Conclusions

In this paper, by applying the principle of minimum potential energy, the governing equations and boundary conditions for a micro-structure-dependent beam theory are derived to study the buckling behaviors of bi-directional FG beams. The proposed micro-beam theory incorporates the size effects based on the modified couple stress theory. The beam theory can be reduced to different kinds of higher-order shear assumptions as well as the two familiar beam theories (namely, Euler-Bernoulli and Timoshenko beam theories) as special cases. A numerical solution procedure, on the basis of a GDQM, is used to calculate the numerical results of the bi-directional FG beams. The numerical results show some conclusions as follows:

(ⅰ) The two exponential FG indexes α and β have a great effect on the critical buckling loads. It indicates that employing the bi-directional FG materials can significantly improve the stiffness of micro-scaled beam structures.

(ⅱ) The critical buckling loads, as well as the stiffness of bi-directional FG micro-beams, increase with the increase in the length scale parameter l and the decrease in the geometrical sizes (L and h). In other words, when the thickness h and length L decrease and fall into micro-scale, the stiffness of micro-beams increases obviously due to the couple stress effects. Thus, the appropriate geometric shape of micro-beams can be designed to reach the desired critical buckling load and to avoid the unstable failure.

(ⅲ) The errors among Euler-Bernoulli, Timoshenko, and higher-order shear beam models increase with the increase in the length scale parameter l and the decrease in the beam length L. When the length to height ratio L/h < 10, the effect of higher-order shear deformation becomes evident, and therefore the higher-order shear distribution functions must be adopted for short beam structures.

(ⅳ) The different boundary conditions have a great influence on the critical buckling loads. The boundary conditions affect the stiffness much evidently especially when the dimensionless length scale parameter l/h is large.

(ⅴ) The buckling analysis of different beam theories incorporating the size effects can be treated in a unified way by considering the proposed micro-structure-dependent beam theory.

Appendix A

This section presents the explicit expressions of the Lagrangian interpolation function and its derivatives, and simply reviews the Hermite interpolation functions ħ_j(x_i) for beam models^[119], where i=1, 2, …, N and j=1, 2, …, N+2. The Lagrangian interpolation functions ℓ_j(x_i) have the following properties:

Shu and Richards^[119] and Quan and Chang^[120] derived the arbitrary order derivatives of the Lagrangian interpolation functions ℓ_j(x_i) as

with

The Hermite interpolation functions ħ_j(x_i) and their derivations can be expressed as following cases:

(ⅰ) For the first sampling point (the end point), we have

where

(ⅱ) For the derivative of the deflection at the first sampling point (the end point), we have

where

(ⅲ) In the case of domain points (j=2, 3, …, N-1), we have

where

(ⅳ) For the last sampling point (the end point), we have

where

(ⅴ) For the derivative of the deflection at the last sampling point (the end point), we have

where