Non-linear modes of vibration of Timoshenko nanobeams under electrostatic actuation

Accepted Manuscript

Non-Linear Modes of Vibration of Timoshenko Nanobeams Under Electrostatic Actuation Marco Alves, Pedro Ribeiro PII: DOI: Reference:

S0020-7403(17)30038-3 10.1016/j.ijmecsci.2017.06.003 MS 3710

To appear in:

International Journal of Mechanical Sciences

Received date: Accepted date:

5 January 2017 4 June 2017

Please cite this article as: Marco Alves, Pedro Ribeiro, Non-Linear Modes of Vibration of Timoshenko Nanobeams Under Electrostatic Actuation, International Journal of Mechanical Sciences (2017), doi: 10.1016/j.ijmecsci.2017.06.003

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

ACCEPTED MANUSCRIPT

Non-Linear Modes of Vibration of Timoshenko Nanobeams Under Electrostatic Actuation

CR IP T

Marco Alves Department of Mechanical Engineering, Imperial College London South Kensington Campus, SW7 2AZ London, UK [email protected]

Pedro Ribeiro DEMec/INEGI, Faculdade de Engenharia, Universidade do Porto, 4200-465 Porto, Portugal [email protected]

AN US

Abstract

PT

ED

M

The non-linear dynamical behaviour of electrostatically actuated nanobeams of rectangular cross section is investigated using a p-version finite element derived by Galerkin’s method, where the partial differential equations of motion are reduced into a finite dimensional system of non-linear ordinary differential equations in the time domain. Timoshenko’s beam theory is applied, as well as the beam analogue of von K´ arm´ an’s plate theory in order to take into account the geometrical non-linearity. The formulation considers non-local effects which affect the inertia of the system, the non-linear stiffness terms and the electrostatic force. Considering several harmonics for the periodic solution, the harmonic balance method is used to transform the ordinary differential equations into algebraic equations of motion in the frequency domain, which are then solved by an arc-length continuation method. Convergence studies show that accurate results are achieved with a reasonably low number of degrees of freedom, and the different terms related to the effects considered in the proposed model are validated with results published in the literature. The importance of the shear deformation and the rotary inertia in this problem is investigated, as well as the influence in the dynamic response of the electrostatic force, fringing fields and non-local effects, combined with the geometrical non-linearity. It is found that different combinations of these effects lead to different outcomes in the system dynamics, changing the natural frequencies and, to a smaller extent, the mode shapes, and leading to hardening, softening or even the combination of both. Furthermore, non-linear phenomena such as internal resonances and bifurcations are studied.

AC

CE

Keywords: Timoshenko nanobeams; Electrostatic; Non-local; Non-linear modes of vibration; Fringing field; Internal resonance.

ACCEPTED MANUSCRIPT 1 Introduction

1

1. Introduction With the evolution of high precision manufacturing techniques and with the pressure to produce even smaller devices, MEMS are being miniaturized day by day, and are now deep in the sub-micrometer range giving birth to nanoelectromechanical systems (NEMS). These electromechanical devices are envisioned for many applications as well as to access interesting and new regimes in physics and engineering, opening a vast number of exploratory research areas. Possible applications for NEMS are mass and force sensing, nano-switches, nano-actuators, and molecules or cell detection, which can be important in specific disease diagnosis. Knowing the resonant frequencies and other dynamic characteristics is key to the success of NEMS in the aforementioned applications [1, 2].

CR IP T

Understanding the motion of MEMS and NEMS is a challenging task, since they can experience deflections that are large when compared to their dimensions, introducing geometric non-linearities. Furthermore, many devices are actuated by an electrostatic force [2, 3], which offers a well-controlled force over the displacement and requires very low current, but is a highly non-linear function of the distance between the two electrodes involved, further enriching the dynamic behaviour of MEMS and NEMS.

AC

CE

PT

ED

M

AN US

Atomistic simulations, such as molecular dynamics, have been extensively used in nano mechanics [4, 5]. In this type of approach, Newton’s equations of motion for all atoms of a system are numerically solved for a certain time interval, in order to obtain an overall view of the dynamical evolution of the system in study. Despite the potentiality of the atomistic simulations, they cannot be used to model a complicated system due to the high computational resources that would be necessary to carry out the simulations [6, 7]. Continuum modelling techniques do not suffer from this handicap and classical elasticity theories (beam, plate, shell) have been applied in attempts to replicate the mechanical behaviour of components used in nanotechnology [8]. However, a few analyses indicate that these classical elasticity theories are not always able to correctly predict the behaviour of nanostructures, because they miss small scale effects [9, 10]. In order to solve this problem, non-local continuum theories were developed, taking into account the interatomic forces, and the internal length scales were introduced into the constitutive equations. In particular, the non-local elasticity theory of Eringen [11] has been used as the basis of many studies, of which only some will be mentioned in this text. Reddy [12] reformulated the classical, first and higher order beam theories using the non-local differential constitutive relations of Eringen, and studied the deflection and the fundamental natural frequency of vibration of nanobeams. Wang et al. [13] applied the non-local elastic constitutive equations and Timoshenko’s beam theory to study the influence of the non-local parameter in the vibration of single-walled carbon nanotubes (SWCNT). Roque et al. [14] applied the non-local elasticity theory of Erigen and a meshless method to study bending, buckling and free vibrations of Timoshenko nanobeams. Ansari and Sahmani [15] computed the natural frequencies of CNTs with different types of chirality and boundary conditions. The natural frequencies computed using non-local continuum beam theories agreed with the ones computed by molecular dynamics simulations. Yang an Lin [16] presented analytical solutions for free vibration of Timoshenko, non-local, nanobeams. By comparisons with molecular dynamics simulations, it was verified that the non-local model provides more accurate values for the natural frequencies of vibration than the local model. Non-local beam theories were applied to compute the linear natural frequencies of vibration of Timoshenko’s and Rayleigh’s beams in, respectively, [17] and [18]. In [19], non-local rod, beam and plate non-local models were employed to obtain linear modes of vibration and frequency response functions. Small-scale effects on the vibration behaviour and stability regions of nanoscale beams with axial velocity were studied in [20]. Generally, non-local beam theories are based on constitutive equations written in differential form; Eptaimeros et al. [21] defined the energy functional of non-local beams in an integral form resorting to Hamilton’s principle, and analysed the differences between this formulation, the differential non-local beam theory and the classical beam theory. Research has also been carried out in the statics and dynamics of electrostatically actuated MEMS and NEMS, with most studies in the sub-micron range addressing CNTs. This is, for example, the case of reference [22], where non-local elasticity theory was used to investigate the deflection and instability of electrostatically actuated CNTs. Due to their different geometry, the expression for the electrostatic force on beams of rectangular cross section differs from the one for CNTs; furthermore, in the former, fringing field effects may be important [23, 24, 25]. Pull-in and dynamic stability of electrostatically actuated micro-beams were analysed by Krylov et al. in [24] and [25]. Kacem et al. [26, 27] investigated pull-in and the frequency response functions in the non-linear regime of NEMS where the moving electrode was a beam of rectangular cross section. The models employed took into account fringing field effects. Ouakad and Younis [28] investigated the dynamic behaviour of double clamped micro-machined arches, also of rectangular cross sections, when actuated by a small DC electrostatic load superimposed to an AC harmonic load. The natural frequencies and mode shapes of the arches were calculated for various values of DC voltages

ACCEPTED MANUSCRIPT 2 Theoretical Model and Solution

2

and initial rises. None of the works on vibrations of electrostatically actuated beams of rectangular cross section analysed the importance of non-local effects.

CR IP T

The modes of vibration provide key information for the correct understanding of the dynamic behaviour of any system. One of the reasons why it is important to know the modes of vibration is that the main resonances that take place in forced responses occur in the neighbourhood of the natural frequencies of vibration, often with the oscillations dominated by a specific mode of vibration. In linear conservative systems, a mode of vibration is defined by a natural frequency and a mode shape that remain unchanged during the vibration cycle. In a non-linear system, like the ones of interest in this paper, this is no longer true. In this work, a non-linear mode is interpreted as a periodic oscillation, not necessarily harmonic, where the correspondent shape and frequency are amplitude dependent. When the amplitudes of vibration tend to low values, the non-linear mode of vibration tends to the linear mode. Another important aspect to consider in the study of non-linear modes is the occurrence of internal resonances [29, 30, 31]: the natural frequencies change with the vibration amplitude and two or more natural frequencies may become commensurable, creating conditions for the interaction of different modes of vibration. It is important to study internal resonances not only in order to understand the dynamics of a specific system, but also because internal resonances may be used in biological and chemical sensors, as well as in RF filters [32].

ED

M

AN US

The non-local theory of Eringen has been applied to investigate non-linear modes of vibration of CNTs [33, 34, 35, 36]. However, of those works only in [36] were non-harmonic vibrations and electrostatic actuation considered. The study performed in the latter reference is based on a Bernoulli-Euler type p-version finite element, which cannot be applied to beams with small thickness to length relations. Furthermore, only CNTs are considered in [36] and, as already pointed out, the electrostatic force in CNTs and in beams with rectangular cross section differ, with the latter beams affected by fringing fields that can be very important [23]. In this paper, a p-version finite element based on Timoshenko’s beam theory and on Eringen’s non-local theory is derived. In order to investigate the nonlinear modes of vibration, free, periodic vibrations of beams subjected to electrostatic forces due to DC voltages are examined. The equations of motion are solved by the harmonic balance method (HBM) and by an arc-length continuation method [30, 37]. It is known, mostly due to research at the macro level, that the p-version finite element has to its advantage the fact that it requires a small number of degrees of freedom for accuracy [31, 35, 36, 38]. It will be shown here that the Timoshenko, non-local, p-element for beams under electrostatic actuation still benefits from this feature. Focus is given to beams of rectangular cross section and, for the first time, the effect of fringing fields on the non-linear modes of vibration of non-local beams is analysed. Softening and hardening spring effects, internal resonances and bifurcations are found; their variation with the DC voltage and with the non-local effects is discussed. The influence of the gap size and of the length of the beams on their dynamics is investigated; Timoshenko and Bernoulli-Euler based non-local, non-linear, models are compared.

PT

2. Theoretical Model and Solution 2.1. Erigen’s non-local elasticity theory

AC

CE

Classical continuum models can describe the evolution of a macroscopic scaled system. However, the length scales - such as lattice spacing between atoms - associated with nanostructures are so small that their effect should be considered in a continuum model, in order that it can be applied, one needs to consider the effect of small lengths, such as lattice spacing between atoms, which makes a consistent model formulation very challenging. To overcome this challenge, the non-local elasticity theory of Erigen has been proposed [11]. This theory takes into account the remote forces between atoms, which causes the stress field at a specific point to depend not only on the strain at that point but also on the strains at all the surrounding points of the domain. The non-local stresses σij are given by σij =

Z

V

K(|x0 − x|, τ )tij (x0 )dx0 .

(1)

where tij (x0 ) is the macroscopic stress tensor at point x’, the attenuation function K(|x0 − x|, τ ) expresses the non-local effects at point x caused by a strain at point x’, with |x0 − x| representing the distance in Euclidean norm and τ is a material constant that depends on the internal and external characteristic lengths. It is possible to represent the constitutive relation presented in Eq. (1) in a differential form [11],

ACCEPTED MANUSCRIPT 2.2

Partial differential equations of motion

3

(1 − µ∇2 )σ = t

µ = (e0 a)2

(2)

where e0 is a material constant, a is the internal characteristic length and ∇2 is the Laplace operator. In the case of Timoshenko beams, the non-local constitutive relations of Eq. (2) can be approximated to a one dimensional form as [12]:

∂ 2 σxx = Exx , ∂x2 ∂ 2 σxz −µ = Gγxz . ∂x2

(3a)

σxz

(3b)

AN US

2.2. Partial differential equations of motion

CR IP T

σxx − µ

(a)

(b)

M

Figure 1: a) Axis and displacement components of the beam element; b) Schematic layout of a double clamped nanobeam under electrostatic actuation.

PT

ED

Based on Timoshenko’s beam theory, the displacements of an arbitrary point in the beam are given as:   u1 (x, z, t) = u(x, t) + zφ(x, t) u2 (x, z, t) = 0   u3 (x, z, t) = w(x, t)

(4)

CE

where φ(x, t), u(x, t) and w(x, t) are, respectively, the rotation of the cross section, the axial and transverse displacement of the point x on the mid-plane of the beam.

AC

Since the transverse deflection is large compared with the beam thickness (leading to geometrical non-linearity) and yet very small when compared to the length L (w << L), the beam analogue of von K´ arm´ an’s plate theory can be applied. The strain-displacement relationship can be expressed as follows ∂u(x, t) 1  ∂w(x, t) 2 ∂φ(x, t) + +z , ∂x 2 ∂x ∂x ∂w(x, t) γxy (x, t) = + φ(x, t). ∂x xx (x, t) =

(5a) (5b)

The equations of motion can be derived by considering Hamilton’s principle, which can be expressed as Z

t2

t1

(δT − δU )dt +

Z

t2

t1

δWnc dt = 0,

(6)

ACCEPTED MANUSCRIPT 2.2

Partial differential equations of motion

4

where T is the kinetic energy, U is the strain energy and Wnc represents the work done by non-conservative forces applied on the system. The potential energy U is given by Z Z

U=

L

(σxx xx + σxz γxz )dAdx,

(7)

A

where A is the cross-section area of the beam, σxx and σxz are the normal and shear stresses, respectively. By replacing the strain-displacement relations of the Eq. (5) in Eq. (7), the potential energy can be written as follows Z Z "

h ∂u

CR IP T

#  ∂w  ∂φ 1  ∂w 2 i + σxx z σxx + + σxz + φ dAdx ∂x 2 ∂x ∂x ∂x L A # Z " h  ∂w  ∂φ ∂u 1  ∂w 2 i 1 + Mx Nx + + Qx + φ dx. = 2 L ∂x 2 ∂x ∂x ∂x

1 U= 2

(8)

where the normal force Nx , the transverse shear force Qx and the bending moment Mx are calculated from Z

σxx dA,

Qx =

A

σxz dA

Mx =

A

The kinetic energy T , can be calculated from 1 T = 2

Z

Z " L

ρA

 ∂u 2 ∂t

Z

σxx zdA.

(9)

A

AN US

Nx =

+ ρA

 ∂w 2 ∂t

+ ρI

 ∂φ 2 ∂t

#

(10)

M

where ρ is the mass density of the beam material and I is the second moment of area of the cross section.

ED

Since the nanobeam is going to be submitted to an electrostatic actuation, an external force fe must be considered. The virtual work of the non-conservative forces, δWnc , is defined by δWnc =

Z

fe δwdx.

(11)

L

AC

CE

PT

Substituting Eq. (8), (10), (11) into Eq. (6), integrating by parts and setting the variational coefficients δu, δw, δφ to zero, one arrives at the equations of motion in terms of the generalized stresses,

∂Nx = 0, ∂x ∂Qx ∂  ∂w  ∂2w + N + fe = ρA 2 , ∂x ∂x ∂x ∂t ∂2φ ∂Mx − Qx = ρI 2 . ∂x ∂t

(12a) (12b) (12c)

In order to simplify the following equations, and since the longitudinal accelerations are much smaller than the ones that occur along the other directions, the longitudinal inertia term was neglected in Eq. (12a). The set of equations of motion presented in Eq. (12) are the same as those for the classical Timoshenko beam theory; however, for a non-local theory, the generalized stresses are different due to the non-local constitutive relation (Eq. (3)), and can be written as,

ACCEPTED MANUSCRIPT 2.3

p-version finite element equations of motion

5

" # ∂u 1  ∂w 2 ∂ 2 Nx = EA , Nx − µ + ∂x2 ∂x 2 ∂x  ∂w  ∂ 2 Qx Qx − µ = k GA + φ , s ∂x2 ∂x 2 ∂φ ∂ Mx = EI , Mx − µ 2 ∂x ∂x

(13a) (13b) (13c)

AN US

CR IP T

where ks is the shear correction factor used in Timoshenko’s beam theory, in order to compensate for the error in assuming a constant shear strain through the thickness of the beam. By substituting Eq. (12) into Eq. (13) the explicit expressions of the non-local generalized stress resultants can be derived: # " ∂u 1  ∂w 2 , (14a) + Nx = EA ∂x 2 ∂x " #  ∂w  ∂ 2  ∂w  ∂fe ∂3w Qx = ks GA − Nx − + φ + µ ρA , (14b) ∂x ∂x∂t2 ∂x2 ∂x dx " # ∂3φ ∂2w ∂  ∂w  ∂φ + µ ρI + ρA 2 − Mx = EI Nx − fe . (14c) ∂x ∂x∂t2 ∂t ∂x ∂x Replacing the stress resultants by the their relation with the generalized displacements, the non-linear and non-local equations of motion are finally derived,

h ∂u 1  ∂w 2 i ∂w  EA + ∂x ∂x 2 ∂x

(15a) !#

(15b)

(15c)

PT

ED

M

 ∂ 2 u ∂w ∂ 2 w  EA + = 0, ∂x2 ∂x ∂x2 "  ∂ 2 w ∂φ  ∂4w ∂3 ∂2w + − µ ρA − ρA 2 − ks GA ∂t ∂x2 ∂x ∂x2 ∂t2 ∂x3 " # h ∂u 1  ∂w 2 i ∂ ∂w  ∂ 2 fe − + = fe − µ 2 , EA ∂x ∂x ∂x 2 ∂x ∂x   2 2  2  ∂ φ ∂w ∂ ∂ φ EI 2 − ks GA + φ = ρI 2 φ − µ 2 . ∂x ∂x ∂t ∂x

CE

2.3. p-version finite element equations of motion

AC

The vector formed by the displacement components u, w and φ can be expressed as the product of the shape functions matrix and the vector of generalized nodal displacements q,    g(x)  u(x, t)   w(x, t) =  0   φ(x, t) 0

   qu (t)   f(x) 0  qw (t) .   qφ (t) 0 t(x) 0

0

(16)

The number of longitudinal, transverse and rotational shape functions are defined by pl , po and pt respectively. Although there are cases where it is necessary to use more than one p-element [39], a single element is sufficient to represent the beams investigated in this paper. Therefore, only the shape functions that satisfy the geometric boundary conditions are included in the model. In the examples that follow, hinged-hinged and, above all, double clamped beams are considered; nevertheless, the p-version finite element can be applied to analyse beams with other boundary conditions, as clamped-free [38]. The shape functions used are derived from the Rodrigues’ form of the Legendre polynomials [31, 40, 41, 42].

ACCEPTED MANUSCRIPT 2.3

p-version finite element equations of motion

6

CR IP T

Inserting expression (16) in the partial differential equations of motion (15), and applying Garlerkin’s method the mass and stiffness matrices can be derived, and the equations of motions can be written as follows,   0    Ku 0 0 0 0 0 ¨u (t)  q     0 ¨w (t) +  0 K0γ22 K0γ23  +  q 0 Mw + Mµw   q ¨φ (t) 0 0 Mφ + Mµφ 0 K0γ32 K0γ33 + K0w (17)       0 K1uw 0 !  qu (t)   0    1 qw (t) = fe 2Kuw + K1µuw K2w + K2µw 0       qφ (t) 0 0 0 0

Mw and Mφ are the local mass matrices associated with the transverse and rotational direction, whereas Mµw and Mµφ are the non-local mass matrices associated to the same directions. In the stiffness matrices, the numbers in superscript indicate the order of dependency of the matrices with the transverse generalised displacement qw (t), where the superscript 0 refers to a constant matrix, superscript 1 refers to a linear dependency, and superscript 2 refers to a quadratic dependency. The subscripts of the stiffness matrices indicate the displacement they refer to, and if non-local effects are considered. The only exception is the matrices K0γ which refer to linear shear stiffness components. The local matrices are given in [43], whereas the non-local matrices are defined in Appendix A

AN US

From the equation of motion corresponding to the longitudinal direction, and in order to reduce the number of equations of motion, one can derive an expression to relate qu with qw given by 1 qu (t) = −(K−1 u Kuw )qw (t)

Electrostatic Force

(18)

The expression that gives the distributed electrostatic force applied in the nanobeam of rectangular cross section, includes a first order fringing field effect correction [44, 45, 24, 23], and is given by, d + w ε0 bV 2  1 + 0.65ζ 2(d + w)2 b

M

fe = −

(19)

CE

PT

ED

where ε0 denotes the permittivity of vacuum, V is the voltage applied, d is the gap between the beam and the electrodes and b represents the width of the beam. The electrostatic force is highly non-linear, and in order to consider its effect in the p-version model, a Taylor expansion, considering the first five terms of the series, was applied: " 2  ε0 bV 2 1 0.65ζ 0.65ζ  0.65ζ  2 3 fe ≈ + − w + + w + 2 d2 bd d3 bd2 d4 bd3 # (20)   0.65ζ  0.65ζ  3 4 4 5 −w + +w + d5 bd4 d6 bd5

AC

The electrostatic force also originates both local and non-local stiffness matrices after Galerkin’s method is applied, which are defined as well in Appendix A. Taking advantage of the relation presented in Eq.(18), and considering all the terms due to the electrostatic force, one derives two systems of equations of motion that are coupled and depend only on qw and qφ .   1 [Mw + Mµw ]¨ qw (t) + K0γ23 qφ (t) + [K0γ22 + K0f e + K0µf e ]qw (t) − [2K1uw + K1µuw ][K−1 K ] u uw qw (t)   (21a) + K2w + K2muw + K1f e + K2f e + K3f e + K1µf e + K2µf e + K3µf e qw (t) + vf e = 0 [Mφ + Mµφ ]¨ qφ (t) + K0γ32 qw (t) + [K0γ33 + Kw ]qφ (t) = 0

(21b)

ACCEPTED MANUSCRIPT 2.4

Solution method

7

The equations of motion contain quadratic, cubic and quartic non-linear terms. Grouping all linear mass terms, all linear stiffness terms and all non-linear stiffness terms for both equations presented above, one can rewrite them as follows,   ¨ w (t) + KLφ1 qφ (t) + KLw1 + KNL(qw (t)) qw (t) + vf e = 0 (22a) Mlw q {z } | KNLw

¨ φ (t) + KLφ2 qφ (t) + KLw2 qw (t) = 0 Mlφ q

CR IP T

2.4. Solution method

(22b)

Since only periodic solutions of the equations of motion are of interest, the HBM finds its natural application, and the generalized solutions can be written in a Fourier series, 3

qw (t) =

u0 X + (ui cos(iωt)) 2 i=1 7

u4 X + (ui cos((i − 4)ωt)) 2 i=5

AN US

qφ (t) =

(23a) (23b)

M

where the vectors ui are the coefficients of the harmonics. One constant term and three harmonics are considered for the periodic solution so the model can accurately represent the periodic oscillations, and also detect typical non-linear phenomena such as internal resonances and bifurcation points. The equations of motion contain odd and even polynomial type of non-linear terms, and since no type of damping was considered, only the cosine terms of ¨w, q ¨ w in the Fourier expansion are necessary. Differentiating qw and qφ twice in order to t, and inserting qw , qφ , q Eq. (21) one arrives to non-differential equations of the form, F(u0 , u1 , u2 , u3 , u4 , u5 , u6 , u7 , ω, t) = 0

(24)

PT

ED

In order to eliminate the cosine terms and the time dependency of the equations, by taking advantage of the orthogonality of the trigonometric functions, the equations are multiplied by each cosine term and then integrated over one period of vibration T . Z 2 T F(u0 , u1 , u2 , u3 , u4 , u5 , u6 , u7 , ω, t) cos(iωt)dt, i = 0, 1, 2, 3 (25) T 0

CE

After applying the HBM, a system of non-linear algebraic equations is obtained, and is given by (−ω 2 MHBM + KHBM )uHBM + fHBM = fel

(26)

where MHM B , KHBM , fHBM , fel of Eq. (26) are presented in Appendix B.

AC

The vector fHBM contains non-linear terms and long expressions would be required in order to present them in detail. In short, each component Fi,HBM of the vector fHBM , is given by, Z 2 T Fi,HBM (ω, u0,1,2,3 ) = Knl qw (cos(iωt))dt, i = 1, 2, 3 (27) T 0 The vector F0,HBM is calculated by Eq. (27), but divided by 2. Matrix Knl contains all non-linear terms of Eq. (22a) and can be separated according to the order of the non-linear dependence of its terms with the solution uHB . The non-linear equations of motion in the frequency domain are given by,   − ω 2 M + Kl + Knl (uHB ) uHB = 0 which are then solved by an arc-length continuation method [30, 46].

(28)

ACCEPTED MANUSCRIPT 3 Numerical Results

8

3. Numerical Results 3.1. Convergence Analysis

po = pt = 2, pl = 4;

po = pt = 5, pl = 7;

po = pt = 8, pl = 10;

po = pt = 12,

ED

Figure 2: Convergence analysis; pl = 14.

M

AN US

CR IP T

When using the p-version of the finite element method, it is important to know the minimum number of degrees of freedom to assure converged results. For this reason, in this section a convergence analysis is made with the model in its full potential, i.e considering all the effects contemplated in this study (non-local, electrostatic and fringing fields). Only convergence is analysed now, the peculiarities of the dynamic behaviour of the system are addressed later.

Figure 2 shows the evolution of the non-dimensional transverse displacement, w ¯ = w/h, computed at the beginning q

CE

PT

ρA of oscillation cycles, with the non-dimensional frequency of vibration of the non-linear system, Ω = ωL2 EI . It can bee seen that convergence has already been achieved when one considers po = pt = 8, pl = 10, and for this reason for all the following analysis presented in the paper these are the numbers of shape functions considered; the corresponding number of degrees of freedom is 16, due to the static condensation performed.

3.2. Model Validation

3.2.1. Validation of the non-local terms

AC

In order to verify the part of the model related with the non-local effects on the stiffness terms, two comparisons were carried out. Although the expression of the electrostatic force presented in section 2.2 is only valid for rectangular cross shaped nanobeams, in the absence of this force the p-version finite element developed can be applied to any symmetric cross shaped beam. This way, in the first example, a hinged-hinged carbon nanotube analysed in [47] was considered. The properties of the CNT are presented in Table 1, where L is the length of the CNT, t is the effective tube thickness, Rext is the external radius and η = eL0 a is the non-dimensional non-local parameter. The values considered for the shear correction factor ks , for the Poisson coefficient ν and the Young modulus E are not presented in [47], instead they can be found in [33] where the authors, who also used the Timoshenko beam theory to study the non-linear free vibration of carbon nanotubes, obtained them using MD simulations. In Figure 3, the frequency response curve obtained by Ghayesh [47] is presented, where the frequency is presented in a non-dimensional form by dividing each value by the first linear natural frequency (ω/ω1 ). The author of [47] considered a CNT with the properties given in Table 1, submitted to a harmonic excitation. The blue line represents

ACCEPTED MANUSCRIPT 3.2

Model Validation

9

Table 1: Properties of CNT for the validation of the non-linear/non-local model.

L [nm]

t [nm]

Rext [nm]

E [TPa]

ρ [Kg/ m3 ]

ν

ks

η

6

0.64

1.50

1.10

1300

0.19

0.563

0.033

Ghayesh [47]; backbone curve calculated at the midpoint of the CNT,

ED

Figure 3: Frequency response curve of the system, p-version model.

M

AN US

CR IP T

the backbone curve (curve which relates the natural frequencies of vibration with the vibration amplitude) obtained by the present model, which is expected to remain in the middle of the frequency response curve of the forced regime. Although the results are very similar, the agreement is not perfect at larger vibration amplitudes, possibly because a Bernoulli-Euler model was employed in [47], whilst a Timoshenko type approach is adopted here.

PT

In the second verification test of this section, the backbone curve computed by the presented model is compared with the one given in [34], for different values of the non-local parameter. The beam has a rectangular cross section, hinged at both ends and with the geometric properties presented in Table 2.

CE

Table 2: Geometric properties of the hinged-hinged beam considered for the validation of the non-linear/non-local terms.

L [µm]

h [µm]

b [µm]

510

1.5

100

AC

The values derived with the p-version model were compared to ones calculated by Simsek [34] and a good agreement was found. Figure 4 shows the variation of the non-linear frequency ratio with the amplitude of vibration for different values of the non-dimensional non-local parameter. It is evident that the non-local effects are gradually more influential on the non-linear frequency values as the vibration amplitude increases. 3.2.2. Validation of the electrostatic force terms When considering the effect of the electrostatic actuation, it is important to understand the static initial conditions of the problem. When tension is applied, the beam deforms while being attracted to the electrode. This deformation occurs until the internal elastic force of the beam balances the electrostatic one, and an equilibrium state is momently achieved. The deflection of this equilibrium state can be called the static deflection. It is important to know the static deflection because it is from this point that the natural modes of vibration are going to be studied. In order to validate the electrostatic force in the model developed, the static deflection of a nanobeam was calculated and compared to published results. The beam analysed has the properties presented in Table 3.

ACCEPTED MANUSCRIPT Model Validation

10

CR IP T

3.2

AN US

Figure 4: Variation of the non-linear frequency ratio with the non-dimensional amplitude (w/r) for different non-local parameter values Simsek [34]; η = 0; η = 0.1; - hinged-hinged nanobeam; η = 0.2; η = 0.3; η = 0.4. Table 3: Properties of the beam for the validation of the electrostatic force terms.

h [µm]

300

2

b [µm]

E [GPa]

ν

20

169

0.28

M

L [µm]

AC

CE

PT

ED

In Figure 5 the solution was calculated using four and five terms of the Taylor series for the electrostatic force, (Eq. 20). The results were compared to the solution obtained by Krylov [24], and accurate results were obtained with the developed model. The pull in voltage calculated by Krylov was 45.1 V , and the value calculated with

Figure 5: Deflection of the mid point of a CC beam with the applied tension; Taylor = 4; Taylor = 5.

Krylov [24];

the presented model was 44.3 V . The instability of equilibria of electrostatically driven MEMS/NEMS devices, commonly referred to as pull-in instability is encountered as a basic static instability mechanism limiting the operational range of these kind of devices in terms of deflection and voltage. The determination of the pull-in

ACCEPTED MANUSCRIPT 3.3

Analysis of the non-linear modes of vibration

11

voltage is obtained under the assumption that the increase of the voltage is slow so that inertia and damping related forces have no type of influence on the structure’s behaviour. 3.3. Analysis of the non-linear modes of vibration Diverse aspects that influence the non-linear modes of vibration of nanobeams with rectangular cross section under the action of electrostatic forces, are analysed in the ensuing sub-sections. 3.3.1. Non-local and voltage effects on the non-linear modes of vibration

CR IP T

In the following analysis, an extensive study of the natural frequencies and mode shapes of a nanobeam under a DC electrostatic force is presented. The nanobeam considered is double clamped, and the importance of the real values the properties employed is not very significant, since the results will always be shown in non-dimensional form. Nevertheless, the elastic properties considered for the following analysis are presented in Table 4, and the geometric properties are given in Table 5. These properties were chosen taking into consideration that one possible way to manufacture NEMS is to scale down silicon-based MEMS [1, 48, 49]. Table 4: Single crystal silicon properties.

ρ [kg/m3 ]

150

ν

AN US

E [GPa]

2300

0.17

Table 5: Nanobeam properties for the natural frequencies and mode shapes analysis.

h [nm]

b [nm]

d [nm]

400

10

10

40

M

L [nm]

ED

To study the non-linear frequencies, the system’s backbone curves are presented for two values of applied voltage V , and three values for the non-dimensional non-local parameter η. The adequate value of the non-local parameter varies from case study to case study [15, 50]; by varying it and employing non-dimensionalised values we provide an overall picture, which is not strictly case-dependent, on its influence. The non-linear natural frequencies q of

AC

CE

PT

ρA vibration and the transverse displacements are also presented in non-dimensional form, respectively: Ω = ωL2 EI and w ¯ = w/h. The transverse displacements are calculated at the centre of the beam and for t = 0, unless stated otherwise.

ACCEPTED MANUSCRIPT Analysis of the non-linear modes of vibration

12

AN US

CR IP T

3.3

AC

CE

PT

ED

M

Figure 6: Total amplitude of vibration of the transverse displacement versus the non-linear natural frequency, when V = 0 V; η = 0; η = 0.075; η = 0.1; static deflection.

Figure 7: Total amplitude of vibration of the transverse displacement versus the non-linear natural frequency, when V = 10 V; η = 0; η = 0.075; η = 0.1; static deflection.

Attending to Figures 6 and 7, it is possible to see the loss of symmetry of the backbone curves in relation to the horizontal axis, when voltage is applied to the system. This loss of symmetry is due to the initial deformed configuration of the beam about which the oscillations occur. In what the natural frequencies are concerned, when no voltage is considered, only hardening spring effect was found, i.e. as the amplitude of the oscillation reaches higher values, the midplane stretching effect stiffens the beam, causing the frequencies to increase. When the system was submitted to a voltage of 10 V, Figure 7, significant changes in the dynamic behaviour of the nanobeam were registered, with the appearance of softening spring effect and turning points, increasing the complexity of the

ACCEPTED MANUSCRIPT 3.3

Analysis of the non-linear modes of vibration

13

backbone curves.

PT

ED

M

AN US

CR IP T

The non-local effects have a significant influence in the backbone curves of an electrostatically actuated nanobeam. We start by analysing the natural frequencies which correspond to small oscillations about the static deflection in Figures 6 and 7. With the exception of the case where no voltage is applied in the system, we will not designate these frequencies as ”linear natural frequencies”, since the equilibrium condition corresponds to an already deflected beam, introducing geometric non-linear terms. One can see that the non-local parameter shifts these frequencies in the figures) to lower values, but this shift is much more pronounced in the absence (solutions represented by of electric tension. To understand the influence of non-local effects on the small amplitude oscillation natural frequencies, it is important to recall which terms of the equations of motion are altered by the non-locality. On one hand, the non-local parameter increases the inertia terms, causing a decrease in the natural frequencies, but, on the other hand, it affects the non-linear stiffness terms, by enhancing them and leading to a stiffness increase in the equilibrium configuration. In order to fully demonstrate that non-local effects tend to increase the non-linear stiffness terms of the equilibrium configuration, the static analysis shown in Figure 8 was performed, where one can see that gradually lower static deflections were obtained as the non-local parameter was increased. To obtain Figure 8, the non-local effects on the electrostatic forces were neglected. Otherwise, the electrostatic force would increase with the non-local parameter, not allowing us to investigate the effect of stiffness alone.

Figure 8: Static deflection evolution of a CC beam for different values of the non-dimensional non-local parameter, when V = 10 V; η = 0; η = 0.075; η = 0.1. Exceptionally, the non-local effects in the electrostatic force are not considered .

AC

CE

Considering now oscillations with increasing amplitudes, one sees that, because the non-local effects enhance the non-linear stiffness terms, they lead to more pronounced hardening or softening spring effects. To fully understand softening and hardening effect, one needs to consider the different stages of one vibration cycle. As we saw, the static deflection of the beam causes geometric non-linear effects, which introduce an initial longitudinal tension to the beam, increasing its stiffness. When the oscillation is in its upward motion, the beam is gradually being submitted to a lower attraction force, since it is farther to the static electrode, and at the same time the longitudinal tension decreases leading to lower stiffness values. As the cycle goes on, the beam passes through its horizontal position, and starts to deflect once again, which reintroduces the longitudinal tension, leading to a new increase in stiffness. When most part of oscillation cycle is dominated by a stiffness decrease, softening effect is detected. However, when most part of the cycle is dominated by an stiffness increase, hardening occurs. Since large vibration amplitudes lead to higher geometric non-linearities (higher axial tensions), softening is usually detected for low amplitude values, where the phenomena that lead to lower natural frequencies prevail. Consequently, when no voltage is applied, only spring hardening is detected as the effects that lead to lower frequency values are absent. In order to better understand the dynamic behaviour related to the backbone curves presented above, the evolution of the amplitude of each harmonic considered in the solution of the periodic response (Eq. 23) is presented for both analysed cases in Figures 9 and 10

ACCEPTED MANUSCRIPT Analysis of the non-linear modes of vibration

14

CR IP T

3.3

(b) 3rd harmonic

AN US

(a) 1st harmonic

Figure 9: Amplitude of the harmonics of the transverse displacement when V = 0 V, for two values of the non-local parameter; η = 0; η = 0.1.

M

When no voltage is considered, the only excited harmonics are the odd ones, and the constant term and the second harmonic have no influence in the system’s response. This occurs since when no electrostatic forces are considered, the system is represented by a non-linear equation of the Duffing type, i.e there are only odd non-linearities. When the electrostatic forces are considered, even non-linearities are originated, and the system is no longer represented by an equation of the Duffing type, thus leading to the appearance of the even harmonics on the system’s dynamic response. In Figure 9 it is clear that the first harmonic is the most dominant one, by comparing the difference of the amplitude range of both harmonics represented.

PT

ED

For the 10 V situation, turning points and softening behaviour were detected in the backbone curves of Figure 7. In Figure 10, the constant term has a strong influence for frequencies close to the linear one, justifying the asymmetric representation of Figure 7. Adding to this, the second harmonic, despite presenting small values of magnitude, also contributes to the asymmetric behaviour of the system. The first harmonic is still the most dominant one, however, the constant term has a relevant influence when low vibration amplitudes are considered. The non-local parameter, once again shifts the frequencies to lower values, and enhances both the softening and hardening effect.

AC

CE

In order to see the evolution of the shapes of the nanobeam with time, the transverse deflection given by Eq. (23 a) is considered. In Figure 11, the transverse oscillations of the nanobeam during half a period of vibration are shown, for a local and a non-local case, and for a point in each backbone curve that presents an amplitude of w ¯ = 0.25. Analysing the shape of the transverse displacement, both local and non-local representations are similar, although in the local case, the maximum displacements are slightly bigger than in the non-local case. Another difference between the local and non-local situation lies in their natural frequency of vibration: 3.3594 × 109 rad/s for the local case, and 3.2195 × 109 rad/s for the non-local one. With this analysis, it becomes clear that the shape assumed by the nanobeam during one vibration period is not the same as the first linear mode shape, due to the influence of the non-linear terms (geometric and electrostatic). 3.3.2. Influence of the beam’s length in the non-linear frequencies The evolution of the backbone curves is dependent on the relation between the attraction of the electrostatic force and the variations in the beam’s stiffness. One important geometric property that significantly affects the dynamic ¯ = L/h. of the system is the length of the nano beam, which will be here represented in a non dimensional form L ¯ (10, 40, 70), two The following analysis considers the properties presented in Tables 4 and 5, three values for L values for the non-local parameter (e0 a = 0 and 10 nm) for an applied voltage of 5 volt. Furthermore, a comparison between the results obtained with the Timoshenko (T) (present model) and the Bernoulli-Euler (B.E) beam theories

ACCEPTED MANUSCRIPT Analysis of the non-linear modes of vibration

15

CR IP T

3.3

(b) 1st harmonic

PT

ED

M

AN US

(a) constant term

(d) 3rd harmonic

CE

(c) 2nd harmonic

AC

Figure 10: Amplitude of the harmonics of the transverse displacement when V = 10 V, for two values of the non-local parameter; η = 0; η = 0.1.

is carried out.

In Figure 12, the effect of changing the nanobeam’s length is evident. The shortest beam is the one that presents a lower static deflection, and consequently only hardening is detected. When comparing the results from both beam theories, one can see significant differences for this beam configuration, since the Bernoulli-Euler theory neglects the transverse shear and the rotary inertia, and therefore is not accurate nor applicable for short beams. When the beam length is increased, one can see that both theories predict the same results, for both the local and non-local situation. Analysing the dynamics of the longest beam, one can see a behaviour that completely differs from the one of the shorter beams, as a significant larger static deflection is detected, and consequently a softening behaviour is detected. In what the non-local effects are concerned, and for the case of the shortest beam, one can see a different dynamic behaviour than the one of its local counterpart, since for Ω ≈ 26 an increase of the hardening rate is detected. This occurs due to a 1:3 internal resonance, where the third harmonic rapidly grows and presents a value

ACCEPTED MANUSCRIPT Analysis of the non-linear modes of vibration

16

CR IP T

3.3

AN US

(a) η = 0

(b) η = 0.1

AC

CE

PT

ED

M

Figure 11: Transverse component of the first transverse mode at different instants along half a period of vibration, when V = 10 V and w ¯ = 0.25.

(a) e0 a = 0 nm

(b) e0 a = 10 nm

Figure 12: Total amplitude of vibration of the transverse displacement versus the non-linear natural frequency, when V = 5 V for ¯ ¯ ¯ ¯ ¯ ¯ different beam lengths; (T) - L=10; (B.E) - L=10; (T) - L=40; (B.E) - L=40; (T) - L=70; (B.E) - L=70.

of opposite sign to the dominant harmonic, which is the first one. The Bernoulli-Euler based model predicts this internal resonance at totally different frequencies and amplitudes. 3.3.3. Influence of the gap size on the non-linear frequencies Analysing Figure 13, the influence of the gap size, which is here presented in a non dimensional form d¯ = d/h, in system dynamics is evident. The situation which considers a lower gap size presents a higher offset, as well as a

ACCEPTED MANUSCRIPT Analysis of the non-linear modes of vibration

17

AN US

CR IP T

3.3

Figure 13: Total amplitude of vibration of the transverse displacement versus the non-linear natural frequency, when V = 5 V; d¯ = 2; d¯ = 4; d¯ = 6; static deflection.

M

small region where no hardening or softening seem to prevail, as a balance between the effects that increase the natural frequencies and the ones that lead to their decrease is achieved. One can also see that a decrease in the gap size does not lead to a linear increase in the offset, as expected since the electrostatic force does not depend linearly on this parameter. 3.3.4. Effect of Fringing fields

ED

In order to see the influence of neglecting the fringe correction factor in Eq. (19), the backbone curve of Figure 7 is presented once again (ζ = 1), this time comparing it to the case where the fringing fields are neglected (ζ = 0).

AC

CE

PT

By analysing Figure 14, the error committed by neglecting the fringing fields is evident. When ζ = 0, the effective attraction force is lower, consequently leading to a lower amplitude offset. When the voltage applied is 5 V, one can already see some significant differences in the systems dynamical behaviour. However, increasing the applied voltage to 10 V, neglecting the fringing field effect leads to a totally different dynamic behaviour of the system. Instead of having softening, only spring hardening was detected and the natural frequencies for the same amplitude of displacement are extremely different. The effect of neglecting the fringing fields is evident in this beam configuration; it is related to the low width to thickness ratio, which causes the fringing fields at the beam’s ends to have a significant contribution to the electrostatic attraction force, when compared to the one of the uniform field that exists between the beam and the static electrode. In Figure 15 a much wider beam was considered. For the 5 V no difference is detected when the fringe effects are neglected, however when the voltage is raised to 10 V, the enhancement of the electrostatic force due to these fields leads to some noticeable differences in system’s dynamic behaviour. After these analysis, one can conclude that the fringing field effects can be neglected for wide non-local beams and when low voltages are applied in the system. For beams with small width to thickness ratio or for larger voltages, these fringed fields cannot be neglected as they significantly enhance the electrostatic force. 3.3.5. Bifurcation Analysis In this section, bifurcation points and secondary branches are studied. The beam considered for this analysis has the elastic properties presented in Table 4 and the geometric properties of Table 5. The influence of the applied voltage and non-local effects on the bifurcation points detected are also analysed.

ACCEPTED MANUSCRIPT Analysis of the non-linear modes of vibration

18

CR IP T

3.3

(b) V = 10 V

AN US

(a) V = 5 V

ζ = 0;

ζ = 1.

CE

PT

ED

M

Figure 14: Effect of neglecting the fringe effect in the backbone curve for ¯b = 1;

(a) V = 5 V

AC

Figure 15: Effect of neglecting the fringe effect in the backbone curve for ¯b = 10;

(b) V = 10 V ζ = 0;

ζ = 1.

In Figure 16, the backbone curves, calculated for a point at the centre of the beam (x = 0), are presented for different combinations of non-local and voltage effects. It is found that the non-local effects tend to shift the bifurcation point to lower amplitude values, for both cases of the applied voltage. This is explained by the facts that - as will be shown in the following paragraphs - the bifurcations are due to interactions between the first and higher order modes, and the latter are more altered by non-local effects than the former, because the distance between nodes decreases as the order of the mode increases. On the other hand, when voltage is changed, the bifurcation points seem to occur at the same frequency ratio, but for different values of amplitude, because the applied voltage introduces an asymmetric behaviour into the system. For all the cases presented in Figure 16, the bifurcation point detected is related to a 1:2 internal resonance, where energy exchanges between the first and second harmonic occur in the secondary branch. But the situations

ACCEPTED MANUSCRIPT Analysis of the non-linear modes of vibration

19

CR IP T

3.3

(b) V = 10 V, η = 0

PT

ED

M

AN US

(a) V = 0 V, η = 0

(d) V = 10 V, η = 0.1

CE

(c) V = 0 V, η = 0.1

AC

Figure 16: Main and secondary branches of the backbone curve of a CC beam for a point at the centre of the beam (x = 0), considering different values of V and η; main branch; secondary branch; bifurcation point.

where no voltage is applied qualitatively differ from the ones where it is. In the absence of electric tension, the system is symmetric and defined by an equation only with odd non-linearities: in the main branch, only the odd harmonics are excited. When a bifurcation occurs and if one follows the secondary branch, the second harmonic appears; therefore, there is a symmetry breaking that is best seen in a phase plane plot (not shown, for the sake of conciseness, a similar example can be seen in [51], for local beams). In addition to the loss of symmetry related to the time evolution, there is a loss of symmetry related with the shapes of vibration. The shape of the odd harmonics, which form the solutions in the main branch, is always symmetric with relation to the z axis. On the other hand, when the secondary branch emerges from the bifurcation point, the second harmonic that appears in the oscillations is connected with an anti-symmetric shape, as presented in Figure 17 (a). When voltage is applied, the dynamics of the system are governed by an equation with a constant term, odd and even non-linearities, thus the even harmonics are already present in the solutions of the main branch, as can be

ACCEPTED MANUSCRIPT 4 Conclusions

20

AN US

CR IP T

seen in Figure 10. However, the bifurcations found in Figures 16 (b) and (d) are still of the symmetry breaking type. To demonstrate this, the shape of the second harmonic is presented in Figure 17 (b), for a solution in the main branch just before the bifurcation occurs and for a solution in the secondary branch. For the point in the main branch, the second harmonic presents a symmetric shape (similar to the shape of first harmonic), however, due to the bifurcation, the shape of the second harmonic in the secondary branch is anti-symmetric.

(a) V = 0 V, η = 0

ω/ωbif = 0.9 (solution in the main branch; the second harmonic does not exist in the ω/ωbif = 1.1 (solution in the secondary branch).

M

Figure 17: Shape of the second harmonic; main branch when V = 0 V);

(b) V = 10 V, η = 0

ED

With or without electric tension, it is evident from Figure 17, that the symmetry breaking bifurcations found are due to a modal interaction involving the first and second modes of vibration.

4. Conclusions

CE

PT

In this work, a p-version finite element based on Timoshenko’s beam theory and on Erigen’s non-local elasticity theory was developed to study the non-linear static and dynamic behaviour of electrostatically actuated nanobeams with a rectangular cross section. The non-linearities considered in the model are related to the large vibration displacements and to the inherent non-linearity of the external electrostatic force. The non-linear equations of motion were solved by the harmonic balance method combined with an arc-length based continuation method, and several harmonics were considered in the Fourier expansion of the periodic solution of the equations of motion.

AC

The diverse parts of the model and the method of solution were validated by comparisons with the literature available. The non-local effects were considered in non-linear vibration analyses and accurate results were obtained when compared to published ones. The new approach was also positively tested in the computation of deflections due to electrostatic forces. Convergence tests demonstrated that accuracy is achieved with a reasonably small number of degrees of freedom. After the validation tests, static analysis were carried out. It was found that the non-local effects increase the non-linear stiffness terms, contributing to a reduction of the magnitude of deflections. However, when the nonlocal effects are considered in the electrostatic force, the effect of the enhanced stiffness is contradicted by the enhancement of the attraction force, which tends to increase static deflections. These static deflections also influence the dynamics, since they are related with initial stresses on the beam and with a deformed configuration about which the oscillations occur. The influence of the applied voltage and non-local effects in the dynamical behaviour of the system was carefully analysed. One found that electric tensions introduce an asymmetric behaviour, with the backbone curves of the

ACCEPTED MANUSCRIPT 4 Conclusions

21

system presenting an amplitude offset, but for low voltages only hardening spring effect occurred. As the voltage was increased to higher values, the system dynamics becomes more complex, and both softening and hardening were found. Non-local effects turn the asymmetric behaviour of the system more noticeable, and enhance both softening and hardening spring effects.

CR IP T

As the Timoshenko beam theory was used, the influence of the beam length in its dynamic behaviour was studied and it was found that, for the same applied voltage, longer beams have a more complex dynamical response, as they present higher deflections and therefore are submitted to higher variations of the electrostatic force during the oscillation cycle. For beams with smaller length to thickness ratio, local or non-local, Bernoulli-Euler theory erroneously predicted the backbone curve. When that ratio was equal to 10, a modal interaction was found solely in the non-local case, both by the Timoshenko and the Bernoulli-Euler model, but the error was large with the latter model. The influence of the distance between the nanobeam and the static electrode was also analysed, since this parameter, along with the voltage, affects significantly the electrostatic force. Softening behaviour was detected for lower values of the gap size and voltage, and as the gap size increased only hardening effect prevailed, even when high voltages were considered.

AN US

The influence of fringed fields in the dynamical behaviour of a nanobeam was also explored. The electrostatic force is related to the existing capacitance between the beam and the static electrode, and when the fringing field effect is considered, the actual existing capacitance is higher, thus leading to higher attraction forces. The influence of fringing fields is more noticeable for narrow beams, where fringing fields can even lead to the appearance of softening, and it becomes negligible as the beam width is increased.

ED

M

Symmetry breaking bifurcation points due to 1:2 internal resonances were detected; they gave birth to secondary branches where energy exchanges between the first and second modes of vibration took place. In all cases, the harmonics of the main branch were related to symmetric shapes but even harmonics related with antisymmetric shapes appeared in the secondary branch. It was found that non-local effects change the frequencies at which the bifurcation points occur, because the non-local effects are stronger in higher order modes, due to the fact that the higher the mode, the smaller the distance between nodes of vibration. It was also found that the applied voltage causes the bifurcations to appear for different deflection amplitudes, mainly because of the offset due to the attraction between the beam and the electrode. Furthermore, the applied voltage changes slightly the characteristics of the bifurcation. In fact, without electric tension, local and non-local systems presented only odd harmonics in the main branch, but, in the presence of electric tension, there are even harmonics already in the main branch and the loss of symmetry is solely related with the shapes the beam assumes during an oscillation cycle.

PT

Acknowledgement

AC

CE

The second author gratefully acknowledges the funding of Project NORTE-01-0145-FEDER-000022 - SciTech Science and Technology for Competitive and Sustainable Industries, cofinanced by Programa Operacional Regional do Norte (NORTE2020), through Fundo Europeu de Desenvolvimento Regional (FEDER).

ACCEPTED MANUSCRIPT 4 Conclusions

22

Appendix A

Mµφ =

2 ρIµ L

Z

1

−1 Z 1

T

∂g(ξ) ∂g(ξ) dξ ∂ξ ∂ξ

(29a)

T

∂f(ξ) ∂f(ξ) dξ ∂ξ −1 ∂ξ Z 1 2 T  ∂ 2 f(ξ) ∂t(ξ) T ∂ f(ξ) 16 q (t) K1µuw = 4 EAµ dξ w L ∂ξ 2 ∂ξ 2 ∂ξ −1 Z 1 2 T 2 ∂ 2 f(ξ) ∂ 2 f(ξ) T ∂ f(ξ) 16 q (t) dξ K2µw = 5 EAµ w L ∂ξ 2 ∂ξ 2 ∂ξ 2 −1 Z ε0 bLV 2  1 0.65ζ  1 vf e = − + f(ξ) dξ 4 d2 bd −1 Z 0.65ζ  1 ε0 bLV 2  2 + f(ξ) f(ξ)T dξ K0f e = − 4 d3 bd2 −1 Z   ε0 bLV 2  3 0.65ζ  1 K1f e = − + f(ξ) f(ξ)T f(ξ) qw (t) dξ 4 3 4 d bd −1 Z 1  2 2 0.65ζ ε bLV 4 0 T + K2f e = − f(ξ) f(ξ) f(ξ) q (t) dξ w 4 d5 bd4 −1 Z  3 0.65ζ  1 ε0 bLV 2  5 T + f(ξ) f(ξ) f(ξ) q (t) dξ K3f e = − w 4 d6 bd5 −1 Z T 0.65ζ  1 ∂f(ξ) ∂f(ξ) ε0 bµV 2  2 + dξ K0µf e = L d3 bd2 ∂ξ −1 ∂ξ Z T  ε0 bµV 2  3 0.65ζ  1 ∂f(ξ) ∂f(ξ)  1 Kµf e = + f(ξ) qw (t) dξ 4 3 L d bd ∂ξ −1 ∂ξ Z 1 T  2 2 ε0 bµV 4 0.65ζ ∂f(ξ) ∂f(ξ)  + f(ξ) q (t) dξ K2µf e = w L d5 bd4 ∂ξ −1 ∂ξ Z T 3 ε0 bµV 2  5 0.65ζ  1 ∂f(ξ) ∂f(ξ)  K3µf e = + f(ξ) qw (t) dξ 6 5 L d bd ∂ξ −1 ∂ξ 2 ρAµ L

AC

CE

PT

ED

M

AN US

CR IP T

Mµw =

(29b) (29c) (29d) (29e) (29f) (29g) (29h) (29i) (29j) (29k) (29l) (29m)

ACCEPTED MANUSCRIPT 4 Conclusions

23

Appendix B

0

0

0

0

0

0

Mlw

0

0

0

0

0

0

4Mlw

0

0

0

0

0

0

9Mlw

0

0

0

0

0

0

0

0

0

0

0

0

0

Mlφ

0

0

0

0

0

0

4Mlφ

0 0  KNLw  0    0   0  =  KLw2   0    0

0

0

0

0

0

0    0   0    0   0    0  0

0

(30a)

CR IP T

KHBM



0

9Mlφ

0

0

0

KLφ1

KNLw

0

0

0

KLφ1

0

0

KNLw

0

0

0

KLφ1

0

0

KNLw

0

0

0

0

0

0

KLφ2

0

0

KLw2

0

0 0

0

AN US

MHBM

 0 0   0  0  = 0  0   0

0

0

KLφ2

0

KLw2

0

0

0

KLφ2

0

KLw2

0

0

0

and the vectors fHBM , and fel are given by,



0    0   KLφ1    0   0    0 

(30b)

KLφ2

(31a)

AC

CE

PT

ED

M

  F0,HBM (ω, u0,1,2,3 )       F1,HBM (ω, u0,1,2,3 )        F2,HBM (ω, u0,1,2,3 )       F3,HBM (ω, u0,1,2,3 ) fHBM = 0         0         0       0     −vf e         0     0        0 fel = 0        0         0        0

0

References

(31b)

ACCEPTED MANUSCRIPT REFERENCES

24

References [1] Darrin J Young, Christian A Zorman, and Mehran Mehregany. MEMS/NEMS devices and applications. In Springer Handbook of Nanotechnology, pages 359–387. Springer, 2010. [2] K. L. Ekinci. Electromechanical transducers at the nanoscale: Actuation and sensing of motion in nanoelectromechanical systems (NEMS). Small, 1(8-9):786–797, 2005. [3] Mohammad I Younis. MEMS Linear and Nonlinear Statics and Dynamics, volume 20. Springer Science & Business Media, 2011.

CR IP T

[4] I. R. Fernandez, H. Fangohr, and A. Bhaskar. Normal modes of carbon nanotubes: Similarities and differences with their continuum counterpart. Journal of Physics: Conference Series, 26(1):131–134, 2006. [5] Kilho Eom, Harold S. Park, Dae Sung Yoon, and Taeyun Kwon. Nanomechanical resonators and their applications in biological/chemical detection: Nanomechanics principles. Physics Reports, 503(45):115 – 163, 2011. [6] D S´anchez-Portal, Emilio Artacho, J M Jos´e M Soler, Angel Rubio, and P Ordej´ on. Ab initio structural, elastic, and vibrational properties of carbon nanotubes. Physical Review B, 59(19):12678, 1999.

AN US

[7] S. Narendar, D. Roy Mahapatra, and S. Gopalakrishnan. Prediction of nonlocal scaling parameter for armchair and zigzag single-walled carbon nanotubes based on molecular structural mechanics, nonlocal elasticity and wave propagation. International Journal of Engineering Science, 49(6):509 – 522, 2011. [8] C.M. Wang, V.B.C. Tan, and Y.Y. Zhang. Timoshenko beam model for vibration analysis of multi-walled carbon nanotubes. Journal of Sound and Vibration, 294(45):1060 – 1072, 2006. [9] B. Arash and Q. Wang. A review on the application of nonlocal elastic models in modeling of carbon nanotubes and graphenes. Computational Materials Science, 51(1):303 – 313, 2012.

M

[10] Q. Wang. Wave propagation in carbon nanotubes via nonlocal continuum mechanics. Journal of Applied Physics, 98(12), 2005.

ED

[11] A. Cemal Eringen. On differential equations of nonlocal elasticity and solutions of screw dislocation and surface waves. Journal of Applied Physics, 54(9):4703–4710, 1983. [12] J. N. Reddy. Nonlocal theories for bending, buckling and vibration of beams. International Journal of Engineering Science, 45(2-8):288–307, 2007.

PT

[13] C M Wang, Y Y Zhang, and X Q He. Vibration of nonlocal Timoshenko beams. Nanotechnology, 18(10):105401, 2007.

CE

[14] C. M C Roque, A. J M Ferreira, and J. N. Reddy. Analysis of Timoshenko nanobeams with a nonlocal formulation and meshless method. International Journal of Engineering Science, 49(9):976–984, 2011.

AC

[15] R. Ansari and S. Sahmani. Small scale effect on vibrational response of single-walled carbon nanotubes with different boundary conditions based on nonlocal beam models. Communications in Nonlinear Science and Numerical Simulation, 17(4):1965–1979, 2012. [16] Yang Yang and C.W. Lim. Non-classical stiffness strengthening size effects for free vibration of a nonlocal nanostructure. International Journal of Mechanical Sciences, 54(1):57 – 68, 2012. [17] J.V.A. dos Santos and J.N. Reddy. Vibration of Timoshenko beams using non-classical elasticity theories. Shock and Vibration, 19:251 – 256, 2012. [18] MA Eltaher, Amal E Alshorbagy, and FF Mahmoud. Vibration analysis of Euler–Bernoulli nanobeams by using finite element method. Applied Mathematical Modelling, 37(7):4787–4797, 2013. [19] S. Adhikari, D. Gilchrist, T. Murmu, and M.A. McCarthy. Nonlocal normal modes in nanoscale dynamical systems. Mechanical Systems and Signal Processing, 6061:583 – 603, 2015.

ACCEPTED MANUSCRIPT REFERENCES

25

[20] Mousa Rezaee and Saeed Lotfan. Non-linear nonlocal vibration and stability analysis of axially moving nanoscale beams with time-dependent velocity. International Journal of Mechanical Sciences, 9697:36 – 46, 2015. [21] K.G. Eptaimeros, C. Chr. Koutsoumaris, and G.J. Tsamasphyros. Nonlocal integral approach to the dynamical response of nanobeams. International Journal of Mechanical Sciences, 115116:68 – 80, 2016. [22] Mir Masoud Seyyed Fakhrabadi, Abbas Rastgoo, and Mohammad Taghi Ahmadian. Size-dependent instability of carbon nanotubes under electrostatic actuation using nonlocal elasticity. International Journal of Mechanical Sciences, 80:144 – 152, 2014.

CR IP T

[23] P. M. Osterberg and S. D. Senturia. M-TEST: A test chip for MEMS material property measurement using electrostatically actuated test structures. Journal of Microelectromechanical Systems, 6(2):107–118, Jun 1997. [24] Slava Krylov. Lyapunov exponents as a criterion for the dynamic pull-in instability of electrostatically actuated microstructures. International Journal of Non-Linear Mechanics, 42(4):626–642, 2007. [25] Slava Krylov, Yuval Gerson, Tali Nachmias, and Uri Keren. Excitation of large-amplitude parametric resonance by the mechanical stiffness modulation of a microstructure. Journal of Micromechanics and Microengineering, 20(1):015041, 2010.

AN US

[26] N. Kacem, S. Hentz, D. Pinto, B. Reig, and V. Nguyen. Nonlinear dynamics of nanomechanical beam resonators: improving the performance of NEMS-based sensors. Nanotechnology, 20(27):275501, 2009. [27] N. Kacem, S. Baguet, S. Hentz, and R. Dufour. Computational and quasi-analytical models for non-linear vibrations of resonant MEMS and NEMS sensors. International Journal of Non-Linear Mechanics, 46(3):532 – 542, 2011. [28] Hassen M. Ouakad and Mohammad I. Younis. The dynamic behavior of MEMS arch resonators actuated electrically. International Journal of Non-Linear Mechanics, 45(7):704–713, 2010.

M

[29] W. Szemplinska-Stupnicka. The Behaviour of Non-linear Vibrating Systems. Kluwer Academic Publishers, 1990.

ED

[30] R. Lewandowski. Nonlinear free-vibrations of beams by the finite-element and continuation methods. Journal of Sound and Vibration, 170:577–593, 1994. [31] P. Ribeiro and M Petyt. Non-Linear Vibration of Beams With Internal Resonance By the Hierarchical FiniteElement Method. Journal of Sound and Vibration, 224(4):591–624, 1999.

PT

[32] Ashwin Vyas, Student Member, Dimitrios Peroulis, and Anil K Bajaj. A Microresonator Design Based on Nonlinear 1 : 2 Internal Resonance in Flexural Structural Modes. Journal of Microelectromechanical systems, 18(3):744–762, 2009.

CE

[33] J. Yang, L. L. Ke, and S. Kitipornchai. Nonlinear free vibration of single-walled carbon nanotubes using nonlocal Timoshenko beam theory. Physica E: Low-Dimensional Systems and Nanostructures, 42(5):1727– 1735, 2010.

AC

[34] M. Simsek. Large amplitude free vibration of nanobeams with various boundary conditions based on the nonlocal elasticity theory. Composites Part B: Engineering, 56:621–628, 2014. [35] A. Houmat. Nonlinear free vibration of non-prismatic single-walled carbon nanotubes by a non-local shear deformable beam p-element. Acta Mechanica, 227(4):1051–1065, 2016. [36] P. Ribeiro. Non-local effects on the non-linear modes of vibration of carbon nanotubes under electrostatic actuation. International Journal of Non-Linear Mechanics, 87:1 – 20, 2016. [37] R. Lewandowski. Computational formulation for periodic vibration of geometrically nonlinear structurespart 1: Theoretical background. International Journal of Solids and Structures, 34(15):1925 – 1947, 1997. [38] S. Stoykov and P. Ribeiro. Nonlinear forced vibrations and static deformations of 3D beams with rectangular cross section: The influence of warping, shear deformation and longitudinal displacements. International Journal of Mechanical Sciences, 52(11):1505–1521, 2010.

ACCEPTED MANUSCRIPT REFERENCES

26

[39] P Ribeiro. Hierarchical finite element analyses of geometrically non-linear vibration of beams and plane frames. Journal of Sound and Vibration, 246(2):225–244, 2001. [40] N. S. Bardell. Free vibration analysis of a flat plate using the hierarchical finite element method. Journal of Sound and Vibration, 151(2):263–289, 1991. [41] W.Han. The analysis of isotropic and laminated rectangular plates including geometrical non/linearity using the p-version finite element method. Ph.d. thesis, University of Southampton, 1993.

CR IP T

[42] Wanmin Han, Maurice Petyt, and Kuo-Mo Hsiao. An investigation into geometrically nonlinear analysis of rectangular laminated plates using the hierarchical finite element method. Finite Elements in Analysis and Design, 18(1-3):273–288, 1994. [43] P. Ribeiro. A p-version, first order shear deformation, finite element for geometrically non-linear vibration of curved beams. International Journal for Numerical Methods in Engineering, 61(15):2696–2715, 2004. [44] R. M C Mestrom, R. H B Fey, J. T M van Beek, K. L. Phan, and H. Nijmeijer. Modelling the dynamics of a MEMS resonator: Simulations and experiments. Sensors and Actuators, A: Physical, 142(1):306–315, 2008.

AN US

[45] R. M. C. Mestrom, R. H. B. Fey, K. L. Phan, and H. Nijmeijer. Simulations and experiments of hardening and softening resonances in a clamped-clamped beam MEMS resonator. Sensors and Actuators, A: Physical, 162(2):225–234, 2010. [46] R Lewandowski. Computational formulation for periodic vibration of geometrically nonlinear structurespart 2: numerical strategy and examples. International Journal of Solids and Structures, 34(15):1949–1964, 1997. [47] Mergen H. Ghayesh. Nonlinear size-dependent behaviour of single-walled carbon nanotubes. Applied Physics A: Materials Science and Processing, 117(3):1393–1399, 2014.

M

[48] Ronald F Gibson, Emmanuel O Ayorinde, and Yuan-Feng Wen. Vibrations of carbon nanotubes and their composites: a review. Composites Science and Technology, 67(1):1–28, 2007.

ED

[49] Samuel Guillon, Daisuke Saya, Laurent Mazenq, Sorin Perisanu, Pascal Vincent, Arnaud Lazarus, Olivier Thomas, and Liviu Nicu. Effect of non-ideal clamping shape on the resonance frequencies of silicon nanocantilevers. Nanotechnology, 22(24):245501, 2011. [50] W. H. Duan, C. M. Wang, and Y. Y. Zhang. Calibration of nonlocal scaling effect parameter for free vibration of carbon nanotubes by molecular dynamics. Journal of Applied Physics, 101(2):1–7, 2007.

AC

CE

PT

[51] S. Stoykov and P. Ribeiro. Nonlinear free vibrations of beams in space due to internal resonance. Journal of Sound and Vibration, 330(18-19):4574–4595, 2011.