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

https://doi.org/10.1007/s00466-021-01987-6 ORIGINAL PAPER

Formulation and experimental validation of space-fractional

Timoshenko beam model with functionally graded materials effects

Paulina Stempin¹·Wojciech Sumelka¹

Received: 20 October 2020 / Accepted: 30 January 2021 / Published online: 13 March 2021

© The Author(s) 2021

Abstract

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.

Keywords Timoshenko beam·Non-local model·Fractional calculus·Microbeam bending Abbreviations

FGM Functionally graded materials

s-FEBB Space-fractional Euler–Bernoulli beam model

s-FTB Space-fractional Timoshenko beam

SE Scale effect

RVE Representative volume element

BD Body dimensions

LS Length scale

FDO Fractional differential operator CTB Classical Timoshenko beam

1 Introduction

The first works on scale effect (SE), which manifest depen- dence 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@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–7] and theoretical concepts [8–13] for SE analysis were considerably developed and moreover SE phenomena were utilised successfully in the modern industry [14–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 represen- tative volume element (RVE) of a specific material to body dimensions (BD).

From the theoretical side, when constructing a mathemat- ical model for the description of mechanical phenomena, one should choose certain mathematical objects proper to the experimental observation scale [19–22]. Herein, we operate on meso/macro level, therefore for SE modelling a phe- nomenological 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–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 intro- duced 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–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 appropri- ate governing equations and enriched numerical algorithm are provided. The discussion also includes a study of the influence of non-locality parameters for different bound- ary 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. Section2 deals with the s-FTB definition. Section3is devoted to the numerical scheme and parametric study. Section4provides experimen- tal validation and finally Sect.5concludes the paper.

2 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 fol- lowing definition of small strain

εi j =1 2^α−_f ¹

_α Dxj

ui(x)+D^α

xi

uj(x)

, (1)

whereui are the components of the displacement vector,x is a spatial variable, whereas the term D( . )^α denotes the Riesz–Caputo fractional derivative [35],

Dα xj

u_i(x)=1 2

(2−α) (2)

C

xj−fD^α_xju_i(x)+(−1)^{n C}xjD^α_x

j+fu_i(x)

, (2)

with the left-side and right-side Caputo derivatives

C

x_j−fD^α_xjui(x)= 1 (n−α)

_xj

x_j−f

u⁽_iⁿ⁾(τ)

(xj−τ)^α−n+1dτ, (3)

C

x_jD^α_xj+fui(x)= −1 (n−α)

x_j+f u⁽_iⁿ⁾(τ)

(xj−τ)^α−n+1dτ,(4)

whereis 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 bound- aries, has been kept. These two parameters α andf are regarded as associated with microstructure [37] and respon- sible for SE mapping.

In the presented study one considers the static bending behavior in the x1x3-plane, therefore, keeping the assump- tion that the cross-section is infinitely rigid in its own plane and remains plane after deformation. In consequence, the displacement field takes the form

u₁(x₁,x₂,x₃)=x₃2, u₂(x₁,x₂,x₃)=0, u₃(x₁,x₃,x₃)= ¯u₃. (5) whereu¯3= ¯u3(x1)is the rigid body translation of the cross- section in 3-rd axis direction of the coordinate system and 2=2(x1)is the rigid body rotation (positive keeping the right-hand rule). Rotation2depends on the Riesz–Caputo fractional derivative with the proportionality factor^α−_f ¹and, by comparison to the already developed s-FEBB [28], is extended by an additional rotation due to the fractional shear deformationγ13,

2= −^α−_f ¹D^α

x₁u¯3+γ13. (6)

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 con- sistency 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

ε11=x₃^α−1f

−D^α

x1

^α−1f

Dα x1u¯₃

+D^α

x1

γ13

, ε13=ε31=1

2

−^α−1f

Dα x1u¯₃+γ13

α

Dx3

x₃+D^α

x1u¯₃ ^α−1f =1

2 γ13.

(7) It is so because in Eq. (72) ^α−_f ¹D^α

x3

x3=1, thenε13 =¹₂γ13. The corresponding stresses are

σ11 =x3^α−_f ¹

−D^α

x1

^α−_f ¹D^α

x1u¯3

+D^α

x1

γ13

E(x2,x3),

σ13 =Gε13, (8)

where E(x2,x3) is Young’s modulus, G = G(x2,x3) =

E(x2,x3)

( +ν) is Kirchhoff modulus andνis Poisson’s ratio. Based

on the above results, the bending momentM2and the shear forceV3can be expressed as

M₂=M₂(x₁)=

A

x₃σ11dA and V₃=V₃(x₁)=

Aσ13dA. (9) and, by introducing Eq. (8), as

M₂=

A

x₃σ11dA=^α−_f ¹

−D^α

x1

^α−_f ¹D^α

x1u¯₃ +D^α

x1

γ13

(E I)^∗, (10)

V₃=

Aσ13dA=k(G A)^∗ γ13, (11)

wherekis the shear correction factor,(E I)^∗= _AE(x2,x3) x₃² dA and (G A)^∗ = _AG(x2,x3) dA. The relationship amongV3, M2 and distributed load p3 = p3(x1)(for the derivation of following formulas from the principle of vir- tual work see Appendix and [28]) is

V3=D^α

x₁

M2^α−_f ¹

, (12)

p3= −D^α

x1

V3^α−_f ¹

. (13)

Moreover, the shear forceV3, by using Eqs. (10) and (12), can be also expressed as

V3=D^α

x₁

²_f^α−2

−D^α

x₁

^α−_f ¹D^α

x₁u¯3

+D^α

x₁

γ13

(E I)^∗. (14)

Finally, substituting Eq. (14) in Eq. (13), and comparing Eq. (14) and Eq. (11) we obtain the fractional equations gov- erning the bending of s-FTB

⎧⎪

⎨

⎪⎩ Dα x₁

^α−_f ¹D^α

x₁

²_f^α−2D^α

x₁

^α−_f ¹D^α

x₁u¯3

(E I)^∗=p3+D^α

x₁

^α−_f ¹D^α

x₁

²_f^α−2D^α

x₁

γ13

(E I)^∗, Dα

x1

²_f^α−2D^α

x1

^α−_f ¹D^α

x1u¯3

E I−D^α

x1

²_f^α−2D^α

x1

γ13

(E I)^∗+k(G A)^∗ γ13 =0. (15)

It is worth noting that when the shear deformation is neglected, i.e.γ13≈0, the rotation Eq. (6) takes the form 2= −^α−_f ¹D^α

x₁u¯3, (16)

and consequently, Eq. (15) reduces to equation governing the bending of a s-FEBB model [28]

Dα x₁

^α−_f ¹D^α

x₁

²_f^α−2D^α

x₁

^α−_f ¹D^α

x₁u¯3

(E I)^∗= p3. (17)

Moreover, forα=1, the rotations Eqs. (6) and (16) take the classical form for Timoshenko and Euler–Bernoulli beam theories, respectively

2= −du¯3

dx1 +γ13 and 2= −du¯3

dx1, (18)

and the governing equations Eqs. (15)1and (17) also reduce to the classical local descriptions

d⁴u¯3

dx₁⁴(E I)^∗=p3+d³γ13

dx₁³ (E I)^∗ and d⁴u¯3

dx₁⁴(E I)^∗=p3, (19) Timoshenko and Euler–Bernoulli beam theories, respec- tively.

3 Numerical study

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–40] was applied to approximate Caputo fractional derivativesD^α

x1( . )iat the node

x₁ⁱby the sum of first-order derivatives at nodesx₁^i−m÷x₁^i+m (see the Fig.1) with appropriate weight coefficients marked asB,C^aandC^b, according to the Eq. (20)

Dα

x₁( . )i =h^1−αA

B( . )i−m+

i−1

j^a=i−m+1

C^a( . )j^a +2( . )i +

i+m−1 j^b=i+1

C^b( . )_jb +B( . )i+m

, (20)

where

m=mi =(f)i/h ≥2, A= (2−α)

2(2)(3−α), B=(m−1)^2−α−(m+α−2)m^1−α, C^a=(i−j^a+1)^2−α−2(i−j^a)^2−α+(i−j^a−1)^2−α, C^b=(j^b−i+1)^2−α−2(j^b−i)^2−α+(j^b−i−1)^2−α, (21)

and(f)i =f(x₁ⁱ).

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 nodex₁ⁱ are given by

Di(E I)^∗= p3(x₁ⁱ)+Ci(E I)^∗, Ci(E I)^∗−Bi(E I)^∗+k(G A)^∗ γ13(xⁱ₁)=0, V3(xⁱ₁)=(−Ci +Bi)(E I)^∗, or V3(x₁ⁱ)=k(G A)^∗ γ13(x₁ⁱ),

M2(x₁ⁱ)=(^α−_f ¹)i(−Bi+Ai)(E I)^∗, 2(x₁ⁱ)= −(^α−_f ¹)iAi+γ13(x₁ⁱ),

(22)

where Di =D^α

x1

(^α−_f ¹)i

Dα x1

(²_f^α−2)i

Dα x1

(^α−_f ¹)i

Dα

x1u¯3(x₁ⁱ)

=D^α

x1

(^α−_f ¹C)i

=h^1−αA

B(^α−_f ¹C)_i−m +

i−1

j^a=i−m+1

C^a(^α−_f ¹C)j^a +2(^α−_f ¹C)i+

i+m−1 j^b=i+1

C^b(^α−_f ¹C)_jb+B·(^α−_f ¹C)i+m

, (23)

Ci =D^α

x1

(²_f^α−2)i

Dα x1

(^α−_f ¹)i

Dα

x1u¯3(xⁱ1)

=D^α

x1

(²_f^α−2B)i

=h^1−αA

B(²_f^α−2B)i−m

+

i−1

j^a=i−m+1

C^a(²_f^α−2B)_j^a +2(²_f^α−2B)_i+

i+m−1 j^b=i+1

C^b(²_f^α−2B)_jb+B(²_f^α−2B)_i+m

, (24)

Bi =D^α

x₁

(^α−_f ¹)i

Dα x₁u¯3(x₁ⁱ)

=D^α

x₁

(^α−_f ¹A)i

=h^1−αA

B(^α−_f ¹A)_i−m +

i−1

j^a=i−m+1

C^a(^α−_f ¹A)j^a +2(^α−_f ¹A)i+

i+m−1 j^b=i+1

C^b(^α−_f ¹A)_jb +B(^α−_f ¹A)i+m

, (25)

Ai =D^α

x1u¯3(x1ⁱ)=h^1−αA

Bu¯_i−m+

i−1

j^a=i−m+1

C^au¯_ja+2u¯_i +

i+m−1 j^b=i+1

C^bu¯_jb+Bu¯_i+m

(26) Ci =D^α

x1

(^α−_f ¹)i

Dα x1

(²_f^α−2)i

Dα x1

γ13(x1ⁱ)

=D^α

x1

(^α−_f ¹B)i

=h^1−αA

B(^α−_f ¹B)i−m

+

i−1

j^a=i−m+1

C^a(^α−_f ¹B)_j^a+2(^α−_f ¹B)_i +

i+m−1 j^b=i+1

C^b(^α−_f ¹B)_jb +B(^α−_f ¹B)_i+m

, (27)

Fig. 1 Discretization of the analysed beam of lengthL—

homogeneous grid: real nodes x₁⁰÷x₁ⁿ; fictitious nodes x₁⁻⁸÷x⁻₁¹andxⁿ₁⁺¹÷xⁿ₁⁺⁸

Bi =D^α

x₁

(²_f^α−2)i

Dα x₁

γ13(x₁ⁱ)

=D^α

x₁

(²_f^α−2A)i

=h^1−αA

B(²_f^α−2A)i−m

+

i−1

j^a=i−m+1

C^a(²_f^α−2)j^aAj^a +2(²_f^α−2A)i+

i+m−1 j^b=i+1

C^b(²_f^α−2A)_jb +B(²_f^α−2A)i+m

, (28)

Ai =D^α

x1

γ13(x₁ⁱ)=h^1−αA

Bγi−m+

i−1

j^a=i−m+1

C^aγj^a+2γi+

i+m−1 j^b=i+1

C^bγ_jb +Bγi+m

. (29)

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 Table1,

( . )_i =

−( . )i+N1−N2 +( . )i+N1+N2

1

2N2h, (30) whereN1andN2determine which finite difference scheme has been applied (see Table1).

One should conclude that similarly to Eq. (19), forα=1 Eq. (22) returns for N1 = 0 and N2 = ¹₂ to the classical central difference scheme (A=¹₂,B=C^a=C^b=0)

i = −−ui−1/2+ui+1/2

h +γi, (31)

Mi =

−ui−1−2ui+ui+1

h² +−γi−1/2+γi+1/2

h

(E I)^∗, (32)

Vi =

−−ui−3/2+3ui−1/2−3ui+1/2+ui+3/2

h³ +γi−1−2γi+γi+1

h²

(E I)^∗, (33)

and

ui−2−4ui−1+6ui−4ui+1+ui+2

h⁴ (E I)^∗= pi +−γi−3/2+3γi−1/2−3γi+1/2+γi+3/2

h³ (E I)^∗, (34)

wherei =2(x₁ⁱ),Mi = M2(x₁ⁱ),Vi =V3(x₁ⁱ), pi = p3(x₁ⁱ),γi =γ13(x₁ⁱ)andui = ¯u3(x₁ⁱ).

Applied numerical methods have resulted in fictitious nodes (x₁⁻⁸÷x₁⁻¹ andx₁ⁿ⁺¹÷x₁ⁿ⁺⁸) outside the beam in addition to real nodes (x₁⁰÷x₁ⁿ) - see Fig.1. The application of the variable length scalef(x), which is decreasing at the boundaries [36], results in only 8 fictitious nodes (x₁⁻⁸÷x₁⁻¹, x₁ⁿ⁺¹÷xⁿ₁⁺⁸) on each side of the beam. These points are elim- inated in final set of equations, by the analogy to the approach presented in [36], by equating of the finite difference approx- imation

Table 1 The applied finite difference schemes evaluated at nodexⁱ₁whereN₁andN₂ determine which finite difference scheme has been applied (see also Eq. (30) and Fig.1)

Forward Backward Central Central

N₁=¹₂, N₁= −¹₂, N₁=0, N₁=0, N₂=¹₂ N₂=¹₂ N₂=¹₂ N₂=1

u_i xⁱ₁=x⁻⁸₁ x₁ⁱ=x₁ⁿ⁺⁸ x₁ⁱ =x₁^−7.5÷x^n+7.5₁ x₁ⁱ=x₁⁻⁷÷x²₁;x₁ⁿ⁻²÷x₁ⁿ⁺⁷ (^α−1f A)i xⁱ₁=x⁻⁶₁ x₁ⁱ=x₁ⁿ⁺⁶ x₁ⁱ =x₁⁻⁵÷x₁ⁿ⁺⁵

(²f^α−2B)i i

1=x₁⁻⁴ x₁ⁱ=x₁ⁿ⁺⁴ x₁ⁱ =x₁^−3.5÷xⁿ₁^+3.5 x₁ⁱ=x₁⁻³÷xⁿ₁⁺³ (^α−1f C)i xⁱ₁=x⁻²₁ x₁ⁱ=x₁ⁿ⁺² x₁ⁱ =x₁⁻¹÷x₁ⁿ⁺¹

γi xⁱ₁=x⁻₁⁶ x₁ⁱ=x₁ⁿ⁺⁶ x₁ⁱ =x₁^−5.5÷xⁿ₁^+5.5 x₁ⁱ=x₁⁻⁵÷x²₁;x₁ⁿ⁻²÷x₁ⁿ⁺⁵ (²_f^α−2A)_i xⁱ₁=x⁻₁⁴ x₁ⁱ=x₁ⁿ⁺⁴ x₁ⁱ =x₁⁻³÷x₁ⁿ⁺³

(^α−1f B)i xⁱ₁=x⁻²₁ x₁ⁱ=x₁ⁿ⁺² x₁ⁱ=x₁⁻¹÷xⁿ⁺¹₁

• with the central and forward schemes for the fourth order derivative of displacement at nodesx₁⁻⁶÷x₁⁻¹and for the third order derivative of strain at nodessx₁⁻⁴÷x₁¹,

ui−2−4ui−1+6ui−4ui+1+ui+2

h⁴ = ui−4ui+1+6ui+2−4ui+3+ui+4

h⁴ , forx₁⁻⁶÷x₁⁻¹,

−γi−2+2γi−1−2γi+1+γi+2

2h³ = −γi+3γi+1−3γi+2+γi+3

h³ , forx₁ⁱ =x₁⁻⁴÷x₁¹;

(35)

• with the central and backward schemes for the fourth order derivative of displacement at nodesx₁ⁿ⁺¹÷xⁿ₁⁺⁶and for the third order derivative of strain at nodesx₁ⁿ⁻¹÷x₁ⁿ⁺⁴

ui−2−4ui−1+6ui−4ui+1+ui+2

h⁴ = ui−4−4ui−3+6ui−2−4ui−1+ui

h⁴ , forx₁ⁱ =x₁ⁿ⁺¹÷x₁ⁿ⁺⁶,

−γi−2+2γi−1−2γi+1+γi+2

2h³ = −γi−3+3γi−2−3γi−1+γi

h³ , forx₁ⁱ =x₁ⁿ⁻¹÷x₁ⁿ⁺⁴.

(36)

3.2 Parametric study

This section highlights the influence of the material param- etersαandf 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 follow- ing data are thereby considered: beam lengthL =100μm, widtha = 10μm and heightb = 50μm of rectangular cross-section, homogenized Young’s modulusE =10 GPa, Poisson’s ratioν=0.2 andh =0.1μm. The shear correc- tion factor for rectangular cross-section is assumedk=5/6.

In each of the examples, the beam was loaded with a con- centrated 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 loadP=10μN is

introduced by equivalent continuous load [28,41]

p3(ξ)= k1k2

2 tanh(k1/2)

1

cosh²[k1(ξ−L1)]

P

L, (37)

where k1 = 100 and k2 are dimensionless parameters, k2 = (L1 +1)⁻²⁰⁰ +1 for L1 ∈ 0.0;0.5] and k2 = (−L1+2)⁻²⁰⁰+1 forL1∈(0.5;1.0],ξ =x1/Lis a dimen- sionless coordinate andL1is a point load position in relation to the beam length (L1 =0.5 for load at the midpoint and L1=1.0 for load at the end of beam). The boundary condi- tions are summarised in Table2. 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 dis- tribution that is smoothly decreasing at the boundaries (see Fig.3), described by the following function [42]