International Journal of Mechanical Sciences 79 (2014) 75–83
Contents lists available at ScienceDirect
International Journal of Mechanical Sciences journal homepage: www.elsevier.com/locate/ijmecsci
A Timoshenko beam element based on the modified couple stress theory M.H. Kahrobaiyan a, M. Asghari a, M.T. Ahmadian a,b,n a b
School of Mechanical Engineering, Sharif University of Technology, Tehran, Iran Center of Excellence in Design, Robotics and Automation (CEDRA), Sharif University of Technology, Tehran, Iran
art ic l e i nf o
a b s t r a c t
Article history: Received 2 June 2013 Received in revised form 7 November 2013 Accepted 14 November 2013 Available online 22 November 2013
Since the classical continuum theory is neither able to evaluate the accurate stiffness nor able to justify the size-dependency of micro-scale structures, the non-classical continuum theories such as the modified couple stress theory have been developed. In this paper, a new comprehensive Timoshenko beam element has been developed on the basis of the modified couple stress theory. The shape functions of the new element are derived by solving the governing equations of modified couple stress Timoshenko beams. Subsequently, the mass and stiffness matrices are developed using energy approach and Hamilton’s principle. The formulations of the modified couple stress Euler–Bernoulli beam element and also classical Timoshenko and Euler–Bernoulli beam elements can be recovered from the original formulations of the new Timoshenko beam element. By two examples, it is indicated that how the new beam element can be applied to deal with the real-case problems. The static deflection of a short microbeam and pull-in voltage of an electrostatically actuated microcantilever made of silicon are evaluated by employing the new beam element and the results are compared to the experimental data as well as the classical FEM results. It is observed that the results of the new beam element are in good agreement with the experimental findings while the gap between the classical FEM and experimental results is notable. & 2013 Elsevier Ltd. All rights reserved.
Keywords: Modified couple stress theory Finite element method Timoshenko beam Size-dependency Length scale parameter Microbeam
1. Introduction
Study of the mechanical behavior of conducting polymer electromechanical actuators (CPEA): Metz et al. [6].
Micro/nano-scale mechanical components, such as micro/ nano-beams, are the major building blocks of micro/nano electromechanical systems (MEMS/NEMS) [1,2] and atomic force microscopes (AFMs) [3,4]. Hence, investigating the mechanical behavior of such components has always been an important issue among the researchers. Due to complications existing in micro/nano-scale systems such as the presence of the complex forces like electrostatic, Casimir, Van Der Waals and capillary forces, complex geometry or some other issues like existence of the squeeze film damping, the exact analytical solutions may not be achieved for the behavior of the mechanical components in many cases; so, some approaches other than analytical one are required. One of the most popular approaches is the finite element method (FEM). The FEM is utilized by many researchers in order to investigate the mechanical behavior of micro/nano-scale systems. For example:
Analysis of piezoelectric cantilever type beam actuators: Wu et al. [5].
n
Corresponding author. Tel.: þ 98 21 6616 5503. E-mail address:
[email protected] (M.T. Ahmadian).
0020-7403/$ - see front matter & 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.ijmecsci.2013.11.014
Investigation of the static behavior and pull-in voltage of
electrostatically actuated cantilever microswitches: Coutu et al. [7]. Analysis of the mechanical response of microswitches with piezoelectric actuation: Cahpius et al. [8] Investigating the dynamic pull-in of an electrostatically actuated micro/nano-plate considering geometrical nonlinearities and fluid pressure by employing a nine-node plate element: Tajalli et al. [9]. Modeling the MEMS subjected to electrostatic forces by developing a non-conforming element: Rochus et al. [10].
It is noted that all the above-mentioned works are based on the elements developed by the classical continuum theory. The experimental observations have indicated that the classical continuum mechanics not only underestimates the stiffness of micro/nano-scale structures but also is incapable of justifying the size-dependency observed in these structures [11–13]; noted that the size-dependency is a peculiar phenomenon in which the normalized mechanical quantities of micro-scale structures that the classical continuum theory predicts to be independent of the structure size, significantly changes by the size. Hence, during past
76
M.H. Kahrobaiyan et al. / International Journal of Mechanical Sciences 79 (2014) 75–83
years, some non-classical continuum theories such as the couple stress theory have been emerged, developed, modified and employed to study the mechanical behavior of the micro-scale structures. The couple stress theory has been introduced and developed by some researchers such as Koiter, Ejike, Mindlin and Tiersten in early 1960s [14–16]. In this theory beside the two classical material constants (i.e. the elastic modulus and Poison’s ratio) additional material parameters are appearing which enable the theory to predict and model the size-dependency in micro-scale structures. Asghari et al. developed a Timoshenko beam model based on this theory to investigate the size-dependent static behavior of microbeams [17]. They concluded that the bending stiffness of the new model is greater than those evaluated based on the classical Timoshenko beam theory. Yang et al. [18] performed a modification on the couple stress theory and presented the modified couple stress theory. They utilized the equilibrium equation of moments of couples in addition to two classical equilibrium equations i.e. the equilibrium equation of forces and moment of forces. Soon after that, this new theory became a popular non-classical theory. Many researchers utilized the modified couple stress theory to develop beam and plate models as well as investigate the size-dependent phenomena in microsystems. Some of these works on developing beam and plate models can be listed as below:
Linear homogenous Euler–Bernoulli beam model by Park and Gao [19] and Kong et al. [20].
Linear homogenous Timoshenko model by Ma et al. [21]. Nonlinear homogenous Euler–Bernoulli beam model by Xia et al. [22] and Kahrobaiyan et al. [23].
Nonlinear homogenous Timoshenko beam model by Asghari et al. [24]
Linear functionally graded Euler–Bernoulli and Timoshenko
important structural element that has been widely investigated in the literature [35–40], a modified couple stress Timoshenko beam element is developed in this paper. The new beam element is a comprehensive beam element that the formulations of the modified couple stress Euler–Bernoulli beam element as well as the classical Timoshenko and Euler–Bernoulli beam elements can be achieved by letting some parameter to zero in the original formulations. The shape functions of the new beam element are derived by directly solving the static equilibrium equations of modified couple stress Timoshenko beams. The stiffness and mass matrices are developed on the basis of the Hamilton’s principle. Some examples are prepared to indicate that how the newly developed beam element can apply to the real-case problems and by comparing the results of the new beam element with the experimental data, it is indicated that the new beam element can successfully capture the size-dependency unlike the classical beam elements. In addition, it is observed that the results of the new beam element are in good agreement with the experimental results whereas the gap between the experimental and the classical FEM outcomes is considerable. The first example deals with the static deflection of a short microcantilever subjected to a concentrated force at its free end. In this example, the results of the modified couple stress and classical Timoshenko and Euler– Bernoulli beam elements are compared to the experimental data and it is observed that the new modified couple stress Timoshenko beam element has the best agreement with experimental findings while the gap between the classical FEM and experimental results is significant. In the second example, the static pull-in voltage of an electrostatically actuated microcantilever made of silicon is determined utilizing the new beam element and the present results are compared to the experimental and the classical FEM results. It is observed that the results of the new beam element are in good agreement with the experimental findings unlike the results of the classical beam element.
beam models by Asghari et al. [25,26].
Linear homogenous Kirchhoff plate model by Tsiatas [27]. In addition to developing beam and plate models, mechanical behavior of microsystems have also been investigated and analyzed based on the modified couple stress theory. Some of these works can be expressed as
Investigating the vibration of fluid-conveying microtubes by Wang [28].
Analyzing the buckling of micro-tubules by Fu and Zhang [29]. Studying the dynamic characteristics of atomic force micro
scopes (AFMs) by Kahrobaiyan et al. [30]. Investigating the size-dependent static behavior of electrostatically actuated microcantilevers and micro-bridges by Rahaeifard et al. [31,32]. Due to:
vast applications micro-scale components such as microbeams in MEMS/NEMS,
necessity of employing the non-classical continuum theories in
order to capture the size-dependency and evaluate reliable stiffness for micro-scale structures, necessity of utilizing the FEM in MEMS/NEMS because of the complexities in these systems,
developing the structural finite elements based on the nonclassical continuum theories seems to be essential. Recently, non-classical bar and Euler–Bernoulli beam elements are developed on the basis of the strain gradient elasticity, a non-classical continuum theory [33,34]. Since Timoshenko beam element is an
2. Preliminaries Strain energy U of an elastic continuum occupying the volume of V modeled by the modified couple stress theory can be mentioned as [18] Z 1 U¼ ðs ε þmij χ ij Þ dV ði; j ¼ 1; 2; 3Þ; ð1Þ 2 V ij ij where sij , εij , mij and χ ij refer to the components of the classical stress and strain tensors, the symmetric part of the couple stress tensor and the symmetric part of the curvature (or rotation gradient) tensor respectively. The strain energy of a modified couple stress continuum consists of two parts: (1) a classical part, i.e. ð1=2Þsij εij and (2) a non-classical part, i.e. ð1=2Þmij χ ij . For an isotropic material, the strain tensor can be related to the stress tensor via Hooke’s law as 1 1 ν εij ¼ ð1 þ νÞsij νskk δij ¼ sij skk δij ; ð2Þ E 2μ 1þν where E and μ represent the elastic (Young’s) modulus and the shear modulus respectively and ν stands for Poisson’s ratio (noted that E ¼ 2ð1 þ νÞμ ). In addition, the couple stress tensor is related to the curvature tensor as [18] 2
mij ¼ 2l μχ ij ;
ð3Þ
in which l denotes the material length scale parameter, an additional material parameter enabling the theory to capture the size-dependency. The strain tensor can be expressed as the
M.H. Kahrobaiyan et al. / International Journal of Mechanical Sciences 79 (2014) 75–83
symmetric part of the displacement gradient tensor ui;j as εij ¼
1 ðu þ uj;i Þ: 2 i;j
ð4Þ
Moreover, the symmetric part of the curvature (rotation gradient) tensor χ ij can be expressed as [18,19] χ ij ¼
1 ðθ þθj;i Þ; 2 i;j
ð5Þ
where θi , which refers to the components of the rotation vector, can be related to the components of the displacement vector field ui as θ¼
1 ∇ u 3 θi ¼ 12 Ε ijk uk;j ; 2
ð6Þ
in which θ and u, respectively, denote the rotation and displacement vectors and ∇ and Ε ijk , respectively, represent the Nabla and Permutation symbols.
3. Developing the stiffness and mass matrices In this section, the stiffness and mass matrices of the new sizedependent modified couple stress Timoshenko beam element are derived. To that end, the appropriate shape functions are developed by directly solving the governing equations of the modified couple stress Timoshenko beam element and the mass and stiffness matrices are obtained by establishing the potential and kinetic energies of the new beam element and then utilizing Hamilton’s principle. In addition, it is indicated that the new beam element is a comprehensive beam element so that the results of the classical Timoshenko beam element, the modified couple stress Euler–Bernoulli beam element and the classical Euler–Bernoulli beam element can be recovered from the formulations of the newly developed modified couple stress Timoshenko beam element. In other words, at the end of this section, it is shown that the mass and stiffness matrices of the new modified couple stress Timoshenko beam element will reduce to those of the modified couple stress Euler– Bernoulli beam element as the ratio of the beam element length to the gyration radius of beam element cross-section increases which consequently guarantees that the shear-locking phenomenon will not happen for the new beam element. Moreover, it is pointed out that as the dimensions of the beam increase, i.e. the ratio of the beam cross-section gyration radius to the material length scale parameter increases, the formulations of the present element approach the formulations of the classical beam elements. In order to develop the new beam element, the displacement field of a Timoshenko beam model (see Fig. 1) can be expressed as u1 ¼ zψ ðx; tÞ;
u2 ¼ 0;
u3 ¼ wðx; tÞ;
ð7Þ
77
where u1 , u2 and u3 represent the displacements along x-, y- and z-axes respectively, z denotes the displacement of an arbitrary point from the neutral axis, ψðx; tÞ stands for the rotation angle of the beam cross-section and wðx; tÞ refers to the lateral deflection of the beam. Substituting Eq. (7) into Eqs. (2)–(6), the non-zero components of the strain, stress, curvature and couple stress tensors as well as the nonzero component of the rotation vector are obtained as [17,21] ∂ψ 1 ∂w ε11 ¼ z ; ε13 ¼ ε31 ¼ ψ ; ∂x 2 ∂x 1 ∂w 1 ∂ψ ∂2 w ; χ 12 ¼ χ 21 ¼ þ 2 ; θ2 ¼ ψ þ 2 ∂x 4 ∂x ∂x ∂ψ ∂w ψ ; s11 ¼ Eε11 ¼ Ez ; s13 ¼ 2με13 ¼ s31 ¼ 2με13 ¼ μ ∂x ∂x 2 β ∂ψ ∂ w þ m12 ¼ 2βχ 12 ¼ m21 ¼ 2βχ 21 ¼ : ð8Þ 2 ∂x ∂x2 Now, substitution of Eq. (8) into Eq. (1) gives the strain (potential) energy U of the new Timoshenko beam element as Z Z 1 L U¼ ðε11 s11 þ 2ε13 s13 þ 2χ 12 m12 Þ dA dx 2 0 A ( 2 2 2 ) Z 2 1 L ∂ψ ∂w μl A ∂ψ ∂2 w ψ þ þ 2 ¼ EI þ μA dx; ð9Þ 2 0 ∂x ∂x 4 ∂x ∂x in which L, A and I respectively denote the length, cross-section area and cross-section area moment of inertia of new beam element and R noted that I ¼ A z2 dA. Furthermore, the kinetic energy T of the beam element is expressed as # Z Z " 1 L ∂u1 2 ∂u2 2 ∂u3 2 T¼ dA ρ þ þ 2 0 A ∂t ∂t ∂t " # 2 Z L 2 1 ∂ψ ∂w ¼ ρ I þA dx; ð10Þ 2 0 ∂t ∂t where ρ refers to the density of the beam element. It is noted that the kinetic energy of the Timoshenko beam element comprises two parts: (1) the rotary kinetic energy caused by rotation of the beam crossRL sections: ð1=2Þ 0 ρIð∂ψ=∂tÞ2 dx and (2) the transitional kinetic energy RL caused by the lateral deformations of the beam: ð1=2Þ 0 ρ Að∂w=∂tÞ2 dx. The work W of the external loads exerted to the beam element can be written as ) Z L Z L ( w )T ( Fðx; tÞ W¼ ðFðx; tÞw Mðx; tÞψÞ dx ¼ dx; ð11Þ Mðx; tÞ ψ 0 0 where Fðx; tÞ and Mðx; tÞ, respectively, represent the external distributed lateral force-per-unit-length and bending moment-per-unitlength. It is noted that the reason of appearing the minus sign “ ” in Eq. (11) is that the direction of the applied distributed couple is in opposition to the direction of the cross-section rotation angle as it is indicated in Fig. 1. In order to obtain the governing equations of the new Timoshenko beam element, Hamilton’s principle can be utilized as follows: Z t2 δðT U þ WÞ dt ¼ 0: ð12Þ t1
Substitution of Eqs. (9)–(11) into Eq. (12) gives the governing equations as 2 ∂2 ψ ∂w μAl ∂2 ψ ∂3 w ∂2 ψ ψ þ þ ð13Þ Mðx; tÞ ¼ ρI 2 ; EI 2 þ μA ∂x 4 ∂x ∂x2 ∂x3 ∂t 2 2 ∂ w ∂ψ μAl ∂3 ψ ∂4 w ∂2 w þ Fðx; tÞ ¼ ρA 2 : þ μA 4 ∂x2 ∂x ∂x3 ∂x4 ∂x Fig. 1. A Timoshenko beam model: kinematic parameters, loadings and coordinate system.
ð14Þ
Now, in order to develop the stiffness and mass matrices of the new beam element, consider a two-node beam element depicted in Fig. 2
78
M.H. Kahrobaiyan et al. / International Journal of Mechanical Sciences 79 (2014) 75–83
Z
L
f¼
Νw Ν
0
T (
Fðx; tÞ Mðx; tÞ
ψ
)
Z
L
dx ¼
fðΝw ÞT Fðx; tÞ ðΝψ ÞT Mðx; tÞg dx:
0
ð23Þ By employing Hamilton’s principle, i.e. substitution of Eqs. (18)(20) into Eq. (12), the governing equation of the new Timoshenko beam element is derived as € Kδ ¼ f; Μδþ
Fig. 2. A two-node Timoshenko beam element: nodal degrees of freedom, geometry and coordinate system.
and two degrees of freedom are assigned to each node: (1) lateral deflection and (2) rotation angle of cross-section. Hence, the nodal displacement vector δ of the beam element is expressed as 9 8 w > > > > > 1> > = < ψ1 > ; ð15Þ δ¼ w > > > > > 2> > > :ψ ;
ð24Þ
It is now clear that Μ and K respectively denote the mass and stiffness matrices of the new beam element and f represents the nodal force vector of the element. By obtaining appropriate shape functions, the mass and stiffness matrices of the new beam element can be derived. Hereafter, the procedure of developing the shape functions is presented. In order to obtain the shape functions, the static governing equations (equilibrium equations) are considered in which the external loads are neglected. Afterwards, the equilibrium equations are solved and appropriate boundary conditions are applied. As the applied boundary conditions, it is assumed that the element has the lateral deflection of w1 and rotation angle of ψ 1 at the first node (x ¼ 0) and lateral deflection of w2 and rotation angle of ψ 2 at the second node (x ¼ L). By applying the aforementioned boundary conditions, the shape functions are achieved. The equilibrium equations of the new beam element can be obtained from Eqs. (13) and (14) by letting ∂=∂t ¼ 0 as follows: EI
2 2 2 d ψ dw μAl d dw ψ þ ¼ 0; þμA ψþ 2 2 dx dx 4 dx dx
ð25Þ
2
where wi and ψ i (i ¼ 1; 2) respectively stand for the lateral deflection and cross-section rotation angle of node i in the beam element. The displacement wðx; tÞ and rotation angle fields ψðx; tÞ can be expressed as " w # ( ) N14 w ¼ δ41 ; ð16Þ Nψ14 ψ 21
w
Nw 2
Nw 3
Nw 4 ;
ψ
Νψ ¼ ½ N 1
Nψ2
Nψ3
N ψ4 :
ð17Þ
θ ¼ 2θ2 ¼ ψ þ
1 T δ Kδ; 2
ð18Þ
dw ; dx
γ ¼ 2ε13 ¼ 2ε31 ¼
dw ψ; dx
ð27Þ
the equilibrium equations mentioned in Eqs. (25) and (26) can be rewritten with respect to the new variables as follows: ! 2 d μAl μAγ θ″ ¼ 0; ð28Þ dx 4
It is noted that these shape functions will be derived by directly solving the static governing equations (equilibrium equations) of the new beam element that will be discussed later. Substituting Eq. (16) into Eqs. (9)–(11) leads to U¼
ð26Þ
noted that in these equations, the external loads are neglected. By introducing the following new variables:
24
in which Nw and Nψ represent the shape function matrices whose components are the appropriate shape functions of the new beam element: Νw ¼ ½ N 1
2 3 d dw μAl d dw μA ψ ¼ 0; ψ þ dx dx dx 4 dx3
2
EI μAl þ 2 4
! θ″
EI γ″ þ μAγ ¼ 0; 2
ð29Þ
where the prime symbol refers to differentiation with respect to x. Integrating Eq. (28) leads to 2
1 T _ T ¼ δ_ Μδ; 2
ð19Þ
W ¼ δT f;
ð20Þ
where dot symbol refers to differentiation with respect to time. In addition, w T w Z L ( ψ T ψ ∂N ∂N ∂N ∂N þ μA Nψ Nψ K¼ EI ∂x ∂x ∂x ∂x 0 ) T 2 μAl ∂Nψ ∂2 Nw ∂Nψ ∂2 Nw þ þ þ dx; ð21Þ 4 ∂x ∂x ∂x2 ∂x2 Z Μ¼
L 0
ρIðNψ ÞT Nψ dx þ
Z
L 0
ρAðNw ÞT Nw dx;
ð22Þ
μAγ
μAl θ″ ¼ c1 ; 4
ð30Þ
in which c1 is a constant. Substitution of γ from Eq. (30) into Eq. (29) results in 2
l ð4Þ 1 c1 θ ð1 þ αÞθ″ ¼ ; 2 8 EI
ð31Þ
where α¼
2 2 μAl μ l ¼ ; E r EI
ð32Þ
in which r denotes the gyration pffiffiffiffiffiffiffiradius of the new beam element cross-section defined as r ¼ I=A. Performing a non-dimensional analysis on Eq. (31) will be helpful. To that end, a dimensionless parameter is introduced
M.H. Kahrobaiyan et al. / International Journal of Mechanical Sciences 79 (2014) 75–83
here: x~ ¼ x=L. Substitution of x~ in Eq. (31) yields 2 4 2 1 l d θ d θ 2c1 L2 ; ð1 þ αÞ 2 ¼ 4 4 L dx~ EI dx~
ð33Þ
As it can be seen in Eq. (33), since the length scale parameter is in the order of a few microns for most of the materials and consequently the ratio of the length scale parameter to the beam 4 length is very small, the coefficient of d θ=dx~ 4 , i.e. ð1=4Þðl=LÞ2 , will 2 be negligible comparing with the coefficient of d θ=dx~ 2 , i.e. ð1 þ αÞ. 4 Hence, by neglecting ð1=4Þðl=LÞ2 ðd θ=dx~ 4 Þ, Eq. (33) reduces to θ″ ¼
2c1 ; EIð1 þαÞ
the shape functions of the new modified couple stress Timoshenko beam element expressed in Eqs. (16) and (17) are derived as
6 x x x 3 x Nψ1 ¼ 1 ; N ψ2 ¼ 1 1 ; Nψ3 Lð1 þ φÞ L L L ð1 þ φÞ L
x 6 x x 3 x 1 ; N ψ4 ¼ 1 1 ; ð42Þ ¼ Lð1 þ φÞ L L L ð1 þ φÞ L 3
x2
x 1 x 2 ; Nw 3 φ 2 1þφ L L L 3
x2
x L x 2 ; Nw ¼ ð4 þ φÞ þ ð2 þ φÞ 3 2ð1 þ φÞ L L L 2
3 1 x x x 3 ; Nw ¼ 2 þφ 4 1 þφ L L L 3
x2
x L x 2 : ¼ þ ðφ 2Þ φ 2ð1 þ φÞ L L L
Nw 1 ¼ 1þ
ð34Þ
then by double-integrating Eq. (34), one can express θ¼
c1 x2 þc2 x þ c3 ; EIð1 þαÞ
ð35Þ
in which c2 and c3 are the constants of integration. Substituting Eq. (34) into Eq. (30) yields
γ¼
c1 1 þ α=2 : μA 1 þ α
ð36Þ 2
MT:I: ¼
6 6 6 26 210ð1 þφÞ 6 4
ð70φ2 þ 147φ þ 78Þ
ρAL
6 6 ρI 6 MR:I: ¼ 30Lð1 þ φÞ2 6 4
L 2 4ð35φ þ77φ þ 44Þ
ð35φ2 þ 63φ þ 27Þ
4L ð35φ2 þ 63φ þ 26Þ
L2 2 4 ð7φ þ14φ þ 8Þ
L 2 4ð35φ þ 63φ þ 26Þ 2
ð70φ þ 147φ þ78Þ
Lð3 15φÞ
36
L2 ð10φ2 þ 5φ þ 4Þ
Lð3 15φÞ
36
36 Symm:
Lð3 15φÞ
1 2
w¼
Z ðγ þ θÞdx ¼
1 c1 1 þ α=2 c 1 x3 c2 x þ x2 þ c 3 x þ c 4 : 2 μA 1 þα 3EIð1 þ αÞ 2
By applying the appropriate boundary conditions, i.e. wðx ¼ 0Þ ¼ w1 ;
ψðx ¼ 0Þ ¼ ψ 1 ;
wðx ¼ LÞ ¼ w2
and
ψðx ¼ LÞ ¼ ψ 2 : ð39Þ
the four constant coefficients of Eqs. (37) and (38), i.e. ci ; i ¼ 1; 2; ::; 4, can be determined as c1 ¼
EIð1 þ αÞ 3
12w1 6Lψ 1 þ 12w2 6Lψ 2 ;
L ð1 þ φÞ ðψ ψ 1 Þ c1 L þ ; c2 ¼ 2 2 L EIð1 þ αÞ c1 1 þ α=2 ; c4 ¼ w1 ; c3 ¼ 2ψ 1 þ μA 1 þ α
ð40Þ
where φ¼
12EI μAL
2
1þ
α E r 2 α ¼ 12 1þ : 2 μ L 2
ð41Þ
By substituting the constants from Eq. (40) into Eqs. (37) and (38),
ð45Þ
3 ð46Þ
where “symm.” implies the symmetric nature of stiffness and mass matrices. In addition, MT:I: and MR:I: respectively denote the tensors of the transitional and rotary inertia that together make the total mass matrix of the Timoshenko beam element M as M ¼ MT:I: þ MR:I: :
ð38Þ
3
7 2 L4 ð7φ2 þ 14φ þ 6Þ 7 7 7; 4L ð35φ2 þ 77φ þ 44Þ 7 5 L2 2 4 ð7φ þ 14φ þ 8Þ
7 L2 ð5φ2 5φ 1Þ 7 7; Lð3 15φÞ 7 5 2 L ð10φ2 þ 5φ þ 4Þ
Now, having θ and γ, the rotation angle ψ and lateral deflection w can be obtained by substitution of Eqs. (35) and (36) into Eq. (27) as ! 1 1 c1 c1 1 þ α=2 2 x þ c2 x þ c3 ψ ¼ ðθ γÞ ¼ ; ð37Þ 2 2 EIð1 þ αÞ μA 1 þ α
ð43Þ
Substituting Eqs. (42) and (43) into Eq. (17), one can obtain the shape function matrices, i.e. Nw and Nψ . Subsequently, substitution of Nw and Nψ into Eqs. (21) and (22) gives the stiffness and mass matrices as 2 3 12 6L 12 6L 6 7 ð4 þ φÞL2 6L ð2 φÞL2 7 EIð1 þ αÞ 6 6 7; ð44Þ K¼ 3 6 12 6L 7 L ð1 þ φÞ4 5 Symm: ð4 þ φÞL2
Symm: 2
79
ð47Þ
The new modified couple stress Timoshenko beam element is a comprehensive beam element that recovers the formulations of modified couple stress Euler–Bernoulli beam element, classical Timoshenko beam element and classical Euler–Bernoulli beam element. By letting α ¼ 0 in Eq. (44), the formulations of the classical Timoshenko beam element can be achieved [40] noted that the condition of α ¼ 0 happens either when one utilizes the classical continuum theory so considers l ¼ 0 in the formulations or when the dimensions of the structure are large, e.g. in macro scales, and consequently the ratio of the length scale parameter to the gyration radius of the beam cross-section l=r becomes negligible. Moreover, the formulations of a modified couple stress Euler–Bernoulli beam element can be obtained by letting φ ¼ 0 in Eqs. (44) and (45). In addition, letting α ¼ 0 and simultaneously φ ¼ 0 in Eqs. (44) and (45), the stiffness and mass matrices of classical Euler–Bernoulli beam element can be derived [40]. It is noted that φ ¼ 0 occurs when the ratio of the beam length to the gyration radius of the beam cross-section is large. Hence, it is concluded that for beams whose length is relatively long, the formulations of the new beam element approach the formulations of Euler–Bernoulli beam elements which
80
M.H. Kahrobaiyan et al. / International Journal of Mechanical Sciences 79 (2014) 75–83
consequently guarantee that the shear-locking phenomenon will not happen for the new beam element.
300 250 200
4. Examples
150
In this section, it is indicated that how the new beam element can be employed in order to solve the real-case problems. The results obtained by utilizing the new beam element are compared to the experimental data as well as the classical FEM outcomes. It is observed that the gap between the classical FEM and experimental results are considerable whereas the results of the new beam element are in good agreements with the experimental observations.
100 50 0
0
200
400
600
800
1000
1200
1400
1600
1800
4.1. A microcantilever under a concentrated force Consider a microcantilever with length Lb having a uniform rectangular cross-section with height h and width b subjected to concentrated force P at its free end (see Fig. 3). The aforementioned microcantilever is modeled by the new beam element and the results are compared to the experimental results extracted from the work done by Lam et al. [12]. They conducted a bending test on a micro-cantilever made of epoxy with the following mechanical properties: the elastic modulus: E ¼ 1:44 GPa and Poisson ratio: ν ¼ 0:38. It is noted that the length scale parameter of epoxy is evaluated to be l ¼ 17:6 μm [19]. Moreover, in the aforementioned bending test, the geometrical properties of the micro-cantilever are reported as follows: the thickness: h ¼ 38 μm, the width: b ¼ 0:235 mm and the ratio of the length to the thickness: Lb =h ¼ 10 [12]. In order to bend the micro-cantilever, Lam et al. [12] used a nano-loading system (Hysitron Triboindenter) to produce a concentrated force at the free end of the cantilever. After that, they graphically reported the experimentally measured end forces versus the end deflections. In order to validate the new non-classical beam element, the aforementioned micro-cantilever is modeled by using 10 new beam elements. The stiffness matrices of the elements are assembled and the boundary conditions, i.e. wð0Þ ¼ ψð0Þ ¼ 0, are applied. The end force applied to the microcantilever is depicted in Fig. 4 versus the end deflection of the microcantilever, wLb . In this figure, the results obtained by the FEM based on the newly established beam elements have been compared to the experimental results and the results obtained by applying the classical beam elements. The results indicate that the outcomes based on the new beam elements are in good agreement with the experimental results while there is a significant difference between the experimental and classical FEM results. It is also observed that the results of the modified couple stress Timoshenko beam element are closer to the experimental data than the results of the modified couple stress Euler–Bernoulli beam element. Since for the aforementioned microcantilever, the ratio of the beam length to the beam height is 10, the above-mentioned observation is justifiable. It is noted that in Fig. 4, in order to evaluate the Euler–Bernoulli-beam-
Fig. 3. A microcantilever subjected to a concentrated force: geometry, loading and coordinate system.
Fig. 4. End force versus the end deflection of the microcantilever for different beam elements.
Fig. 5. An electrostatically actuated microcantilever: geometry and configuration.
element-based results, φ ¼ 0 is considered in formulations of the new beam element (Eqs. (44)–(46)). It is inferred that in order to model the micro-scale structures, employing the non-classical elements is essential. The good agreement between the current and the experimental results implies that the newly developed beam elements are valid, reliable and can be successfully employed to deal with the mechanical problems in micron and sub-microns scales. It is seen that the non-classical beam elements predict that the beams are stiffer than those of the classical beam theories.
4.2. Pull-in voltage of an electrostatically actuated microcantilever Here, the newly developed beam element is employed to evaluate the pull-in voltage of electrostatically microbeams. Consider an electrostatically actuated microcantilever with length Lb an initial distance d from a fixed substrate as shown in Fig. 5. The cross-section of the beam is assumed to be rectangular with width b and height h and the applied voltage is denoted by V. As the applied voltage increases, the attractive electrostatic force between the fixed substrate and the microcantilever increases and so does the deflection of the microcantilever, as the consequence. At a certain voltage, the instability begins and the microcantilever tends to collapse to the fixed substrate. The aforementioned instability is well-known in the literature as the pull-in phenomenon and the corresponding voltage is recognized as the pull-in voltage denoted by V P . Here, the static pull-in voltage of a microcantilever made of silicon is determined by employing the newly developed beam elements. For the microbeam depicted in Fig. 5, the distributed electrostatic force
M.H. Kahrobaiyan et al. / International Journal of Mechanical Sciences 79 (2014) 75–83
per unit length is expressed as ε0 bV 2 d w ; 1 þ0:65 FðxÞ ¼ b ðd wÞ2
81
ð48Þ
where ε0 ¼ 8:854 10 12 represents the vacuum permittivity and V stands for the applied voltage between the beam and the substrate. Substituting F from Eq. (48) into Eq. (23) and considering M to be zero, the nodal force vector can be calculated. Assembling the stiffness matrix and force vector of all elements and also satisfying the boundary conditions, i.e. wð0Þ ¼ 0 and ψð0Þ ¼ 0, the equation of the static deflection of the electrostatically actuated microcantilever is determined as follows: ~ 1 f~ ; δ~ ¼ K
ð49Þ
where “ ~ ” symbol refers to the assembled version of the respected quantity. Eq. (49) is a nonlinear equation in which the ~ force vector f~ is a nonlinear function of the displacement vector δ. This equation can be solved using an iterative method as follows. At first, the electrostatic load on the un-deformed microbeam (i.e. 0 δ~ ¼ 0 which yields w0 ðxÞ ¼0) is calculated as ε0 bV 2 d : ð50Þ F 0 ðxÞ ¼ 1 þ 0:65 2 b d
Fig.6. Comparing the present and the experimental results of static pull-in voltage for silicon microbeams fabricated in 110 direction.
0
Given F 0 ðxÞ, the force vector f and consequently the first estima1 tion of the nodal displacement vector (i.e. δ~ ) and subsequently the displacement field (i.e. w1 ðxÞ) can be calculated using Eq. (49). Now, the next estimation of the electrostatic load and force vector can be obtained by substituting w1 ðxÞ into Eq. (48) and again the next estimation of the nodal displacements vector and displacement field can be found. This procedure will be stopped when the convergence is observed or pull-in instability is happened. The convergence criterion is considered as error i o error desired ; in which ~ i ~ i 1 δ δ error i ¼ ; ~ i δ
ð51Þ Fig. 7. Comparing the present and the experimental results of static pull-in voltage for silicon microbeams fabricated in 010 direction.
error desired ¼ 10 8 ;
ð52Þ
and pull-in happens if wmax Z1. In order to indicate the advantages of the new beam element, the results obtained using the new beam element are compared with those reported in the experimental research performed by Osterberg [41]. The specifications of the microbeam tested by Osterberge [41] are presented in Table 1. It is noted that the microbeams were fabricated in two different directions of silicon crystal: 110 direction; i.e. the length of the beam is along the 110 direction and the side plane of the beam normal to 110 direction of silicon crystal and 010 direction; i.e. the length of the beam is along the 010 direction and the side plane of the beam normal to 010 direction of silicon crystal [41,42]. Figs. 6 and 7 compare the results of the pull-in voltage evaluated by the new beam element with those evaluated based on the classical FEM and also the experimental observations [41] Table 1 Specifications of the silicon microcantilevers tested by Osterberg [41]. Specification
Group 1
Group2
Crystal direction along beam length Elastic modulus along beam length [43] Poisson's ratio in side plane of the beam [43] The range for length ( Lb) Height (h) width (b) Distance from the base (d)
110 169.2 GPa 0.239 75–250 μm 2.94 μm 50 μm 1.05 μm
010 130.4 GPa 0.177 75–225 μm 2.94 μm 50 μm 1.05 μm
for silicon microbeams respectively in 110 and 010 directions. It is noted that in order to generate the graphs of Figs. 6 and 7, the length scale parameter of the silicon is considered to be l110 ¼0.58 mm and l010 ¼0.71 mm in crystal planes normal to 110 and 010 directions respectively. These length scale parameters are evaluated by matching the numerical results with the experimental data which is in agreement with the results of Rahaeifard et al. [31]. The length scale parameter is in fact a material characteristic which relates the couple stresses to the curvature of the continuum. This parameter is different for each material and can be determined by performing some standard experimental tests such as micro-bending tests and micro-torsion tests [11–13]. In microbending (micro-torsion) tests, the clamped-free samples, i.e. microbeams with rectangular cross-sections (circular micro-bars) with different thicknesses (diameters) are subjected to the bending (torsion) loads. Afterward, the graphs of normalized bending (torsional) rigidity versus the sample thickness (diameter) are delineated and by curve-fitting of these graphs with the results obtained on the bases of the modified couple stress theory, the length scale parameter can be determined [11–13]. In a recent work, the length scale parameter of nickel, aluminum and copper are determined by performing micro-bending and micro-torsion tests in which the loading process continued until the yield point of the samples. Afterwards, the length scale parameters are obtained by comparing the experimentally obtained yield loads with those derived on the basis of the modified couple stress theory [44]. In addition, by comparing the modified couple stress theory with dislocation theory, a relation is developed between
82
M.H. Kahrobaiyan et al. / International Journal of Mechanical Sciences 79 (2014) 75–83
the length scale parameter and the dislocation based physical quantities such as barrier strength of the boundary for slip transmission, Burger vector length and shear modulus in that work [44]. Moreover, the length scale parameter of high-strength concretes is related to the grain-size of these materials [45]. Figs. 6 and 7 indicate that the classical FEM underestimates the pull-in voltage of the microcantilever. On the other hand, the pullin voltages evaluated by the new beam element are indeed in good agreement with the experimental observations [41]. Hence, it comes to conclusion that the new modified couple stress beam elements can reduce the gap between the experimental observations and simulation results and therefore the necessity of using the non-classical-continuum-based FEM such as the present modified couple stress based beam element will be comprehensible.
5. Conclusion The classical continuum theory not only underestimates the stiffness of the micro-scale structures, but also is unable to justify the size-dependent mechanical behavior observed in these structures. So, the non-classical continuum theories such as the modified couple stress theory have recently become very popular in order to deal with the modeling of micro-scale structures. Accordingly, developing the new finite elements based on the non-classical continuum theories seems to be essential to achieve that goal. In this paper, a new Timoshenko beam element is developed based on a non-classical continuum theory, the modified couple stress theory. The shape functions are obtained by solving the governing equations and applying the proper boundary conditions. In addition, the stiffness and mass matrices are derived employing an energy-based approach. It is observed that the new modified couple stress Timoshenko beam element is a comprehensive beam element that the formulations of the modified couple stress Euler– Bernoulli beam element, classical Timoshenko beam element and classical Euler–Bernoulli beam element can be achieved from the original formulations. The mass and stiffness matrices of the modified couple stress Euler–Bernoulli beam element are obtained from the respected matrices of the new element when the ratio of the beam length to the gyration radius of the beam cross-section increases. In addition, the mass and stiffness matrices of the classical Timoshenko beam element are derived from those of the new Timoshenko beam element when the ratio of the gyration radius of the beam cross-section to the material length scale parameter increases. Moreover, the formulations of the classical Euler–Bernoulli beam element are achieved from the present formulations when both the aforementioned ratios increase. By some examples, it is indicated that how the new beam element can be utilized in practice. In first example, the static deflection of a short microcantilever is determined by employing the new Timoshenko beam element and the result is compared to the experimental data as well as the results obtained by utilizing the modified couple stress Euler–Bernoulli beam element and classical beam elements. It is observed that the gap between the experimental and classical FEM results is notable whereas the results of the modified couple stress beam elements are in good agreement with the experimental outcomes. It is noted that among the modified couple stress beam elements, the new Timoshenko beam element result is closer to the experimental findings than the new Euler–Bernoulli beam element result which confirms that the Timoshenko beam elements are more suitable for modeling the short beams. In second example, the pull-in voltage of an electrostatically actuated microcantilever made of silicon is evaluated using the newly developed beam elements and the results are compared to the classical FEM results and experimental data. It is observed that the new FEM results are in excellent agreement with
the experimental data while the gap between the classical FEM and experimental results are significant. It is also observed that as the size of the beam increases, the difference between the new and the classical FEM results decreases.
References [1] Attia P, Tremblay G, Laval R, Hesto P. Characterisation of a low-voltage actuated gold microswitch. Mater Sci Eng B 1998;51:263–6. [2] Moeenfard H, Ahmadian MT. Analytical modeling of bending effect on the torsional response of electrostatically actuated micromirrors. Optik – Int J Light Electron Opt 2012;124(12):1278–86, http://dx.doi.org/10.1016/j.ijleo.2012.06.025. [3] Kahrobaiyan MH, Rahaeifard M, Ahmadian MT. Nonlinear dynamic analysis of a V-shaped microcantilever of an atomic force microscope. Appl Math Modell 2011;35:5903–19. [4] Kahrobaiyan MH, Ahmadian MT, Haghighi P, Haghighi A. Sensitivity and resonant frequency of an AFM with sidewall and top-surface probes for both flexural and torsional modes. Int J Mech Sci 2010;52:1357–65. [5] Wu DH, Chien WT, Yang CJ, Yen YT. Coupled-field analysis of piezoelectric beam actuator using FEM. Sens Actuators A 2005;118:171–6. [6] Metz P, Alici G, Spinks GM. A finite element model for bending behaviour of conducting polymer electromechanical actuators. Sens Actuators A 2006;130:1–11. [7] Coutu RA, Kladitis PE, Starman LA, Reid JR. A comparison of micro-switch analytic, finite element, and experimental results. Sens Actuators A 2004;115:252–8. [8] Chapuis F, Bastien F, Manceau JF, Casset F, Charvet PL. FEM modelling of Piezoactuated microswitches. In: Seventh international conference on thermal, mechanical and multiphysics simulation and experiments in micro-electronics and micro-systems. EuroSime 2006, IEEE; 2006. p. 1–6. [9] Tajalli SA, Moghimi Zand M, Ahmadian MT. Effect of geometric nonlinearity on dynamic pull-in behavior of coupled-domain microstructures based on classical and shear deformation plate theories. Eur J Mech A Solids 2009;28:916–25. [10] Rochus V, Rixen D, Golinval JC. Non-conforming element for accurate modelling of MEMS. Finite Elem Anal Des 2007;43:749–56. [11] Fleck NA, Muller GM, Ashby MF, Hutchinson JW. Strain gradient plasticity: theory and experiment. Acta Metall Mater 1994;42:475–87. [12] Lam DCC, Yang F, Chong ACM, Wang J, Tong P. Experiments and theory in strain gradient elasticity. J Mech Phys Solids 2003;51:1477–508. [13] Stölken JS, Evans AG. A microbend test method for measuring the plasticity length scale. Acta Mater 1998;46:5109–15. [14] Koiter WT. Couple-stresses in the theory of elasticity: I and II. Proc Nederl Akad Wetensch Ser B 1964;67:17–29. [15] Mindlin RD, Tiersten HF. Effects of couple-stresses in linear elasticity. Arch Ration Mech Anal 1962;11:415–48. [16] Ejike UBCO. The plane circular crack problem in the linearized couple-stress theory. Int J Eng Sci 1969;7(9):947–61. [17] Asghari M, Kahrobaiyan MH, Rahaeifard M, Ahmadian MT. Investigation of the size effects in Timoshenko beams based on the couple stress theory. Arch Appl Mech 2011;81(7):863–74. [18] Yang F, Chong ACM, Lam DCC, Tong P. Couple stress based strain gradient theory for elasticity. Int J Solids Struct 2002;39(10):2731–43. [19] Park SK, Gao XL. Bernoulli–Euler beam model based on a modified couple stress theory. J Micromech Microeng 2006;16(11):2355–9. [20] Kong S, Zhou S, Nie Z, Wang K. The size-dependent natural frequency of Bernoulli–Euler micro-beams. Int J Eng Sci 2008;46:427–37. [21] Ma HM, Gao XL, Reddy JN. A microstructure-dependent Timoshenko beam model based on a modified couple stress theory. J Mech Phys Solids 2008;56:3379–91. [22] Xia W, Wang L, Yin L. Nonlinear non-classical microscale beams: static bending, postbuckling and free vibration. Int J Eng Sci 2010;48(12):2044–53. [23] Kahrobaiyan MH, Asghari M, Rahaeifard M, Ahmadian MT. A nonlinear strain gradient beam formulation. Int J Eng Sci 2011;49(11):1256–67. [24] Asghari M, Kahrobaiyan MH, Ahmadian MT. A nonlinear Timoshenko beam formulation based on the modified couple stress theory. Int J Eng Sci 2010;48 (12):1749–61. [25] Asghari M, Ahmadian MT, Kahrobaiyan MH, Rahaeifard M. On the sizedependent behavior of functionally graded micro-beams. Mater Des 2010;31 (5):2324–9. [26] Asghari M, Rahaeifard M, Kahrobaiyan MH, Ahmadian MT. The modified couple stress functionally graded Timoshenko beam formulation. Mater Des 2011;32(3):1435–43. [27] Tsiatas GC, Yiotis AJ. A microstructure-dependent orthotropic plate model based on a modified couple stress theory. Recent developments in boundary element methods. Southampton: WIT Press; 2010; 295–308 (A volume to honour Professor John T. Katsikadelis). [28] Wang L. Size-dependent vibration characteristics of fluid-conveying microtubes. J Fluids Struct 2010;26(4):675–84. [29] Fu Y, Zhang J. Modeling and analysis of microtubules based on a modified couple stress theory. Physica E 2010;42:1741–5. [30] Kahrobaiyan MH, Asghari M, Rahaeifard M, Ahmadian MT. Investigation of the size-dependent dynamic characteristics of atomic force microscope microcantilevers based on the modified couple stress theory. Int J Eng Sci 2010;48 (12):1985–94.
M.H. Kahrobaiyan et al. / International Journal of Mechanical Sciences 79 (2014) 75–83
[31] Rahaeifard M, Kahrobaiyan MH, Asghari M, Ahmadian MT. Static pull-in analysis of microcantilevers based on the modified couple stress theory. Sens Actuators A 2011;171(2):370–4. [32] Rahaeifard M, Kahrobaiyan MH, Ahmadian MT, Firoozbakhsh K. Sizedependent pull-in phenomena in nonlinear microbridges. Int J Mech Sci 2012;54(1):306–10. [33] Kahrobaiyan MH, Asghari M, Ahmadian MT. Longitudinal behavior of strain gradient bars. Int J Eng Sci 2013;66 and 67:44–59. [34] Kahrobaiyan MH, Asghari M, Ahmadian MT. Strain gradient beam element. Finite Elem Anal Des 2013;68:63–75. [35] Thomas DL, Wilson JM, Wilson RR. Timoshenko beamfinite elements. J Sound Vib 1973;31(3):315–30. [36] Lees AW, Thomas DL. Unified Timoshenko beam finite element. J Sound Vib 1982;80(3):355–66. [37] DAWE DJ. A finite element for the vibration analysis of Timoshenko beams. J Sound Vib 1978:11–2060(l) 1978:11–20. [38] Davis R, Henshell RD, Warburton GB. A Timoshenko beam element. J Sound Vib 1972;22(4):475–87.
83
[39] Rakowski J. The interpretation of the shear locking in beam elements. Comput Struct 1990;37(5):769–76. [40] Friedman Z, Kosmatka JB. An improved two-node Timoshenko beam finite element. Comput Struct 1993;47(3):473–81. [41] Osterberg PM., Electrostatically actuated micromechanical test structure for material property measurement [PhD dissertation]. Massachusetts Institute of Technology; 1995. [42] Osterberg PM, Senturia SD. M-TEST: a test chip for MEMS material property measurement using electrostatically actuated test structures. J Microelectromech Syst 1997;6:107–18. [43] Gupta RK., Electrostatic pull-in test structure design for in-situ mechanical property measurement of microelectromechanical systems (MEMS) [PhD dissertation]. Massachusetts Institute of Technology; 1997. [44] Kahrobaiyan MH, Rahaeifard M, Ahmadian MT. A size-dependent yield criterion. Int J Eng Sci 2014;74:151–61. [45] Sherafatnia K, Kahrobaiyan MH, Farrahi GH. Size-dependent energy release rate formulation of notched beams based on a modified couple stress theory. Eng Fract Mech, http://dx.doi.org/10.1016/j.engfracmech.2013.12.001, in press.