Formulation and experimental validation of space-fractional Timoshenko beam model with functionally graded materials effects

In this study, the static bending behaviour of a size-dependent thick beam is considered including FGM (Functionally Graded Materials) effects. The presented theory is a further development and extension of the space-fractional (non-local) Euler–Bernoulli beam model (s-FEBB) to space-fractional Timoshenko beam (s-FTB) one by proper taking into account shear deformation. Furthermore, a detailed parametric study on the influence of length scale and order of fractional continua for different boundary conditions demonstrates, how the non-locality affects the static bending response of the s-FTB model. The differences in results between s-FTB and s-FEBB models are shown as well to indicate when shear deformations need to be considered. Finally, material parameter identification and validation based on the bending of SU-8 polymer microbeams confirm the effectiveness of the presented model.


Introduction
The first works on scale effect (SE), which manifest dependence on the answer of the system (e.g. deformation of a material body -the crux of this paper) in relation to its dimensions, date to the investigations of Leonardo da Vinci In Honor of Professor Tomasz Łodygowski on the Occasion of His 70th Birthday.
B Wojciech Sumelka wojciech.sumelka@put.poznan.pl 1 Institute of Structural Analysis, Poznan University of Technology, Piotrowo 5 Street, 60-965 Poznan, Poland at the very beginning of XVI century [1,2]. Since that time both experimental techniques [3][4][5][6][7] and theoretical concepts [8][9][10][11][12][13] for SE analysis were considerably developed and moreover SE phenomena were utilised successfully in the modern industry [14][15][16][17][18]. The main message is that from the standpoint of mechanics, which constitutes the central point of presented considerations, one can say that each structure, due to the complexity of material over different scales of observation, reveals SE. On the other hand, it is clear that the strength of SE phenomena is different depending on the analysed case and proportional to the ratio of the representative volume element (RVE) of a specific material to body dimensions (BD).
From the theoretical side, when constructing a mathematical model for the description of mechanical phenomena, one should choose certain mathematical objects proper to the experimental observation scale [19][20][21][22]. Herein, we operate on meso/macro level, therefore for SE modelling a phenomenological approach is used, thus in consequence RVE to BD ratio is mapped utilising so-called length scale (LS) parameter (it is clear that LS meaning is different depending on certain theory [15,[23][24][25][26]). To be precise, the developed s-FTB theory is defined in the framework of space-Fractional Continuum Mechanics (s-FCM) [27,28], where LS is introduced through the fractional differential operator (FDO) and furthermore, an additional parameter which controls SE is introduced, namely order of FDO. Finally, because of the complex nature of mechanical properties through the beam thickness FGM concept is also used [29][30][31][32].
This paper is a continuation of previous studies focused on the development of non-local beam models [27,28]. The currently formulated space-Fractional Timoshenko beam model, compared to the space-Fractional Euler-Bernoulli beam model, takes into account the shear effect to extend modelling range for thick beams. The definition of rotation has been therefore extended by additional rotation resulting from non-local shear deformations. Hence, the appropriate governing equations and enriched numerical algorithm are provided. The discussion also includes a study of the influence of non-locality parameters for different boundary conditions on the static beam response to bending, a comparison study between s-FTB and s-FEBB models, and identification and validation of s-FTB model parameters for the experimental data of the SU-8 polymer microcantilever bending test (including FGM effects).
The paper is structured as follows. Section 2 deals with the s-FTB definition. Section 3 is devoted to the numerical scheme and parametric study. Section 4 provides experimental validation and finally Sect. 5 concludes the paper.

Theory
The space-Fractional Timoshenko beam (s-FTB) theory is an extension of the space-Fractional Euler-Bernoulli beam (s-FEBB) theory [28] to include shear deformation and make it suitable for thick beams as mentioned in the introduction. As in the previous study, the non-locality is taken into account, based on the fractional elasticity concept [33,34], by the following definition of small strain where u i are the components of the displacement vector, x is a spatial variable, whereas the term α D( . ) denotes the Riesz-Caputo fractional derivative [35], with the left-side and right-side Caputo derivatives where is an Euler gamma function, n = [α] + 1, and [.] denotes the integer part of a real number, α ∈ (0, 1] is the order of FDO, and f is LS i.e. the surrounding affecting the considered material point. The concept of variable length scale [36] f = f (x), as function decreasing at the boundaries, has been kept. These two parameters α and f are regarded as associated with microstructure [37] and responsible for SE mapping.
In the presented study one considers the static bending behavior in the x 1 x 3 -plane, therefore, keeping the assumption that the cross-section is infinitely rigid in its own plane and remains plane after deformation. In consequence, the displacement field takes the form is the rigid body translation of the crosssection in 3-rd axis direction of the coordinate system and 2 = 2 (x 1 ) is the rigid body rotation (positive keeping the right-hand rule). Rotation 2 depends on the Riesz-Caputo fractional derivative with the proportionality factor α−1 f and, by comparison to the already developed s-FEBB [28], is extended by an additional rotation due to the fractional shear deformation γ 13 , Herein, it should be emphasised that on one hand side the assumption Eq. (6) reflects the influence of microstructure on beam cross section rotation, but simultaneously acts as a consistency condition. Namely, it allows to obtain proper relation between the component of the fractional Cauchy strain ε 13 and the fractional shear deformation γ 13 . The last statement causes that in limit case when s-FTB reduces to s-FEBB ε 13 = 0 which is fundamental for s-FEBB [28]. Next, using Eqs (1) and (6) the nonzero fractional Cauchy strains are It is so because in Eq.
x 3 = 1, then ε 13 = 1 2 γ 13 . The corresponding stresses are where 2(1+ν) is Kirchhoff modulus and ν is Poisson's ratio. Based on the above results, the bending moment M 2 and the shear force V 3 can be expressed as 13 dA. (9) and, by introducing Eq. (8), as where k is the shear correction factor, The relationship among V 3 , M 2 and distributed load p 3 = p 3 (x 1 ) (for the derivation of following formulas from the principle of virtual work see Appendix and [28]) is Moreover, the shear force V 3 , by using Eqs. (10) and (12), can be also expressed as Finally, substituting Eq. (14) in Eq. (13), and comparing Eq. (14) and Eq. (11) we obtain the fractional equations governing the bending of s-FTB It is worth noting that when the shear deformation is neglected, i.e. γ 13 ≈ 0, the rotation Eq. (6) takes the form and consequently, Eq. (15) reduces to equation governing the bending of a s-FEBB model [28] α Moreover, for α = 1, the rotations Eqs. (6) and (16) Timoshenko and Euler-Bernoulli beam theories, respectively.

Discretization
Equation (15) has been solved utilising the numerical method [38]. The beam was discretized in n intervals of length h (see Fig. 1). The trapezoidal rule [38][39][40] was applied to approximate Caputo fractional derivatives where and . Based on the Eq. (20), the numerical representation of the distributed load Eq. (15), shear force Eq. (14), bending moment Eq. (10) and rotation Eq. (6) at node x i 1 are given by where Fig. 1 Discretization of the analysed beam of length Lhomogeneous grid: real nodes x 0 1 ÷ x n 1 ; fictitious nodes It should be pointed out that the first derivatives ( . ) in Eqs. (23 ÷ 29) are approximated by forward, backward or central difference formulae at the relevant nodes according to Eq. (30) and Table 1, where N 1 and N 2 determine which finite difference scheme has been applied (see Table 1). One should conclude that similarly to Eq. (19), for α = 1 Eq. (22) returns for N 1 = 0 and N 2 = 1 2 to the classical central difference scheme and 1 ÷ x n+8 1 ) outside the beam in addition to real nodes (x 0 1 ÷ x n 1 ) -see Fig. 1. The application of the variable length scale f (x), which is decreasing at the boundaries [36], results in only 8 fictitious nodes (x −8 1 ) on each side of the beam. These points are eliminated in final set of equations, by the analogy to the approach presented in [36], by equating of the finite difference approximation Table 1 The applied finite difference schemes evaluated at node x i 1 where N 1 and N 2 determine which finite difference scheme has been applied (see also Eq. (30) and Fig. 1) • with the central and forward schemes for the fourth order derivative of displacement at nodes x −6 1 ÷ x −1 1 and for the third order derivative of strain at nodess x −4 1 ÷ x 1 1 , • with the central and backward schemes for the fourth order derivative of displacement at nodes x n+1 1 ÷ x n+6 1 and for the third order derivative of strain at nodes x n−1

Parametric study
This section highlights the influence of the material parameters α and f on the bending behaviour of the beam and compares the fractional and classical approaches to show the ability of taking into account the SE. A comparison of s-FTB and s-FEBB is provided as well to emphasize the need of considering the shear effect when the beam is thick in relation to its length. The examples of beam with following data are thereby considered: beam length L = 100 μm, width a = 10 μm and height b = 50 μm of rectangular cross-section, homogenized Young's modulus E = 10 GPa, Poisson's ratio ν = 0.2 and h = 0.1 μm. The shear correction factor for rectangular cross-section is assumed k = 5/6. In each of the examples, the beam was loaded with a concentrated force: at the mid-span for simply supported, fixed, and propped cantilever schemes, and at the free end for a cantilever scheme (see Fig. 2). The point load P = 10 μN is introduced by equivalent continuous load [28,41] where k 1 = 100 and k 2 are dimensionless parameters,  Table 2. The effect of non-locality was investigated for the following parameters: α ∈ {0.8, 0.6} and max f ∈ {0.001L, 0.10L, 0.2L} with a symmetric distribution that is smoothly decreasing at the boundaries (see Fig. 3), described by the following function [42] f ( where The results of this study, presented in Fig. 4, show the effect of a small scale on the deflection. It can be observed that parameters α and f influence the stiffness of the beam. With a decrease in α or an increase in max f , the deflection of the simply supported and cantilever beam increases, while for propped cantilever and fixed beams the stiffening or softening effect is observed for certain values of these parameters. For length scale f , small in relation to L ( max f = 0.001L), the effect of non-locality disappears and the results are identical to the classic local formulation, as expected.
Moreover, in the above examples, the length to height ratio is L/b = 2, and it is visible that for this beam geometry, the deflections predicted by s-FTB and s-FEBB differ significantly. Beam deflections predicted by s-FTB and s-FEBB theories in the wider range of L/b = 2 ÷ 26 are compared in Fig. 5. The maximum deflections are higher for s-FTB than the s-FEBB model, and similarly to the classic approach, the differences are noticeable for thick beams and negligible for slender beams. Moreover, the inclusion of the non-locality makes the shear effect more significant than in the classic approach. It is especially visible for the fixed beam scheme, where, for example, for the ratio L/b = 2, the maximum deflection for the Timoshenko beam is about 3.85 times greater than for the Euler-Bernoulli beam in the classical approach (α = 1) and 4.85 in fractional approach with α = 0.6 and max f = 0.1L. Therefore, the shear effect for small scale thick beam definitely should not be ignored and should even be considered for a higher ratio L/b than in local theory. However, for a higher L/b ratio, the ratio of deflections u s−FT B u s−F E B B reaches a value of 1.0. For this reason, for long beams, the s-FTB model can be reduced to the s-FEBB model without losing the correctness of the results -which is analogical to local Timoshenko and Euler-Bernoulli beams.

Experimental validation
To show the effectiveness of the developed model, it was validated based on the microcantilever bending experiment presented in [43]. The tested samples were made of a SU-8 polymer with the following material parameters: the bulk value of Young's modulus E b = 6.9 GPa and Poisson's ratio ν = 0.22. These cantilevers had a length in the range of L = 80 ÷ 850 μm and a rectangular cross-section with a width of a = 82 μm and a = 122 μm, and a height of b = 8.4 μm and b = 14.4 μm. These data allow for the analysis of non-local effects depending not only on the crosssection dimensions but also on the length of the beam. As shown in Fig. 8, decreasing the length of the beam causes a decrease in modulus of elasticity compared to the bulk value. Conversely, the decrease of the cross-section dimensions is manifested by the increase of Young's modulus. Moreover, we took this increase into account by analogy to the results obtained for iPP+0.2% PACS polymer samples. Figure 6 shows the distribution of elastic modulus in cross-section measured for samples from the nanoindentation test, while Young's modulus from the tensile test is 1.75 ± 0.1 GPa (results derived from [44]). This map of Young's modulus indicates that it is not constant in the whole cross-section, but there is a core with other characteristics. Based on this information, we include the change of elastic modulus in x 3 direction by the core-shell (CS) model (see Fig. 7) with a stiffer core (E c = f E b ) of thickness c and a bulk shell (E b ), where f means the ratio of the core Young's modulus E s and bulk Young's modulus E b . The stiffness of the core-shell model can be expressed by the effective bending stiffness where E e is the effective Young's modulus and I e , I c and I b are the moments of inertia of the effective (homogeneous) cross-section, the core, and the bulk shell respectively, I e = ab 3 12 , I c = ac 3 12 , Then, the effective Young's modulus is It was found that for SU-8 polymer rectangular samples the height of core zone is c = 7.0 μm, which has a Young modulus f = 1.91 times higher than the bulk value, regardless of sample size. Subsequently, using the Eq.   Table 2 Boundary conditions applied for the static schemes presented in Fig. 2 (see also Eq. (30) and Fig. 1)

Beam type Boundary conditions
Simply supportedū 3 The validation results are shown in Fig. 8, where the elastic modulus was calculated according to Timoshenko beam theory for a cantilever with a point load at the free end It is observed, that the s-FTB model results agree very well with the measurement data. In addition, the results from the classical Timoshenko beam (CTB) model are plotted in Fig. 8, demonstrating that the classical approach is not suitable for describing the behavior of small scale beams.

Conclusions
In the presented work, the space-Fractional Timoshenko beam theory has been developed from space-Fractional Euler-Bernoulli beam theory by including the shear deformation. The effect of non-locality for bending of different beam types was studied and the results of the space-Fractional Timoshenko and space-Fractional Euler-Bernoulli beam theories were compared. The validation of the model was also provided. From these analyzes, the following conclusions are drawn: • the numerical algorithm has been enriched with the shear effect while keeping the possibility of including any boundary conditions, any transverse load, and variable length scale, • parameters α and f control the stiffness of the fractional beam, • the shear effect is more significant in non-local beams, therefore it should be considered even with more slender beams than in a local approach, • the current model is adequate for modeling small scale short beams, • the validation confirms the applicability of the presented model.    (43) and introducing Eqs. (10), (11) and α−1 (44) and δW is the variation in the external work δW = L p 3 δū 3 dL.
To obtain the equilibrium equations, we look for a minimum of the functional J which, using the fractional Euler-Lagrange equation [45,46], provides The boundary conditions for selected beam types are presented in Table 3.