Vibration analysis of Euler–Bernoulli nanobeams embedded in an elastic medium by a sixth-order compact finite difference method

Accepted Manuscript

Vibration analysis of Euler–Bernoulli nanobeams embedded in an elastic medium by a sixth-order compact finite difference method S.A. Mohamed , R.A. Shanab , L.F. Seddek PII: DOI: Reference:

S0307-904X(15)00543-0 10.1016/j.apm.2015.08.019 APM 10708

To appear in:

Applied Mathematical Modelling

Received date: Revised date: Accepted date:

6 September 2013 4 June 2015 26 August 2015

Please cite this article as: S.A. Mohamed , R.A. Shanab , L.F. Seddek , Vibration analysis of Euler– Bernoulli nanobeams embedded in an elastic medium by a sixth-order compact finite difference method, Applied Mathematical Modelling (2015), doi: 10.1016/j.apm.2015.08.019

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 Highlights

Vibration analysis of a nonlocal Euler–Bernoulli beam embedded in an elastic medium. Pasternak elastic foundation model - Nonlocal differential elasticity of Eringen. Sixth-order accuracy schemes for governing equation and boundary conditions. The proposed 6th order scheme is simple and outperforms similar existing methods. Parametric study for nonlocal parameter, slenderness ratios, and boundary conditions.

AC

CE

PT

ED

M

AN US

CR IP T

    

1

ACCEPTED MANUSCRIPT

Vibration analysis of Euler–Bernoulli nanobeams embedded in an elastic medium by a sixthorder compact finite difference method S. A. Mohamed, R. A. Shanab (b) and L.F. Seddek (c) Department of Engineering Mathematics, Faculty of Engineering, Zagazig University, Egypt

CR IP T

Department of Engineering Mathematics, Faculty of Engineering, Zagazig University, Zagazig, Egypt. e-mail: [email protected] Department of Engineering Mathematics, Faculty of Engineering, Zagazig University, Zagazig, Egypt. email: [email protected] Abstract

ED

M

AN US

This paper presents an efficient sixth-order finite difference discretization for vibration analysis of a nonlocal Euler–Bernoulli beam embedded in an elastic medium. The Pasternak elastic foundation model is utilized to represent the surrounding elastic medium. Nonlocal differential elasticity of Eringen is exploited to reveal the nun-locality effect of nanobeams. Sixth-order accuracy schemes are developed for discretization of both the governing equation and boundary conditions. Sixth-order accuracy schemes are derived for simply supported, clamped and free boundary conditions. Numerical results include comparison with exact solutions and with previously published works are presented for the fundamental frequencies. In addition, numerical results are presented to figure out the effects of nonlocal parameter, slenderness ratios, and boundary conditions on the dynamic characteristics of the beam.

1. Introduction

PT

Key Words: nanobeam, 6th order, Fundamental frequencies, nonlocal, Pasternak foundation.

CE

In recent years, nanobeams and carbon nanotubes (CNTs) have attracted worldwide attention due to their outstanding mechanical, chemical and electrical properties [1–5]. Unfortunately, the classical continuum theories are deemed to fail for these nanostructures, because the

AC

length dimensions at nano-scale are often sufficiently small [6]. Nonlocal elasticity theory, which was first developed by Eringen [7–10], has been proposed to predict the accurate behavior of these nanostructures by considering size-dependent effects. According to this theory, the stress state at a reference point is considered as a function of strain states of all points in the body. This theory involves the information about the long-range forces between atoms and introduces the internal length scale in the constitutive equation as a material parameter.

2

ACCEPTED MANUSCRIPT Reddy [11] reformulated various available beam theories, including the Euler–Bernoulli, Timoshenko, Reddy, and Levinson beam theories by using the nonlocal differential constitutive relations of Eringen to study analytically bending, vibration, and buckling behaviors of nanobeams. Aydogu [12] proposed a generalized nonlocal beam theory to study bending, buckling, and free vibration of nanobeams using the nonlocal constitutive equations of Eringen in the formulations. Both Reddy [11] and Aydogu [12] obtain analytic solution for a simply supported nanobeam in the form of infinite power series and reported results of the

CR IP T

first 100 terms. Ansari et al. [13] utilized a high order compact finite difference method for discretization of a nonlocal beam embedded in an elastic medium to obtain the fundamental frequencies corresponding to various sets of boundary conditions. Eltaher et al. [14] presented an efficient finite element model for vibration analysis of a nonlocal Euler– Bernoulli beam.

AN US

In the present paper, the vibration analysis of nonlocal nanobeams is investigated. Euler– Bernoulli beam theory incorporated with nonlocal differential equation of Eringen [7–10] is used to derive the nonlocal differential equations of motion. 6th order compact finite difference schemes are derived for the governing differential equations and different

M

boundary conditions.

The paper is organized in six sections. Section 2 presents the mathematical formulation for

ED

the nonlocal Euler–Bernoulli beam model to get the governing equation. The basic formulations and derivation of the proposed sixth-order finite difference discretization for the beam is presented in Section3. Section 4 shows the derivation of the six-order difference

PT

equations for the three commonly used boundary conditions, namely, simply supported, clamped and free. Section 5 is devoted to the numerical results to demonstrate the accuracy of

CE

the current numerical method. Finally, Section 6 contains the conclusions.

AC

2. Mathematical formulation Consider a straight uniform beam with the length L and a rectangular cross-section of thickness h. A coordinate system (x, y, z) is introduced on the central axis of the beam, where the x axis is taken along the length of the beam, the y axis in the width direction and the z axis along the height direction. Also, the origin of the coordinate system is selected at the left end of the beam (Fig. 1). 𝑍

0

𝑥

ℎ 𝑏

𝐿 3

ACCEPTED MANUSCRIPT

Figure (1): A straight uniform nanobeams with rectangular cross-section. It is assumed that the deformations of the beam take place in the 𝑥– 𝑧 plane, so the displacement components (𝑢1 , 𝑢2 , 𝑢3 ) along the axis (𝑥, 𝑦, 𝑧) are only dependent on the x and z coordinates and time t. In a general form, the following displacement field can be written: 𝑢1 (𝑥, 𝑧, 𝑡) = −𝑧

𝜕𝑤(𝑥,𝑡) 𝜕𝑥

,

𝑢2 (𝑥, 𝑧, 𝑡) = 0 , 𝑢3 (𝑥, 𝑧, 𝑡) = 𝑤(𝑥, 𝑡)

(1)

CR IP T

where w is the transverse displacement of the beam. The governing equation relevant to Euler–Bernoulli beam theory can be expressed in terms of the displacement as 𝜕4 𝑤

𝜕2 𝑤

𝜕4 𝑤

−𝐸𝐼 𝜕𝑥 4 = 𝑄 + 𝜌𝐴 𝜕𝑡 2 − 𝜌𝐼 𝜕𝑥 2 𝜕𝑡 2

(2)

AN US

where 𝐸 , 𝐼 , 𝜌, 𝐴 and 𝑄 are the elastic modulus, moment of inertia, mass density, cross sectional area and transverse load of the nanobeam, respectively. The elastic surrounding is simulated using the Pasternak foundation model which considers both normal pressure and transverse shear stress. So the loading corresponding to this type of foundation model yields 𝜕2 𝑤

𝑄 = 𝐾𝑤 𝑤 − 𝐾𝑠 𝜕𝑥 2 ,

(3)

M

where 𝐾𝑤 is the Winkler modulus parameter corresponding to normal pressure and 𝐾𝑠 is the Pasternak modulus parameter relevant to transverse shear stress.

ED

According to Eringen’s nonlocal elasticity theory, the stress at a point x in a body depends not only on the strain at point x but also on those at all other points of the body. Thus, the nonlocal stress tensor 𝜎 at point x is expressed as

PT

𝜎 = ∫𝑉 𝛼(|𝑥̅ − 𝑥|, 𝜏) 𝑇(𝑥̅ )𝑑𝑥̅

(4)

AC

CE

where T(x̅) is the classic stress tensor at point x̅ given by T(x̅) = 𝐶(𝑥): 𝜖(𝑥). Here, ϵ(x) is the strain tensor, C(x) is the fourth-order elasticity tensor and ': ' denotes the ‘double-dot product’, 𝛼(|𝑥̅ − 𝑥|, 𝜏) is the nonlocal modulus or attenuation function incorporating into the constitutive equations the nonlocal effects at the reference point x produced by the local strain 𝑒 𝑎 at the source 𝑥̅ . Also, |𝑥̅ − 𝑥| is the Euclidean distance, and 𝜏 = 0𝑙 is defined as small scale factor where 𝑒0 is a constant to adjust the model to match the reliable results by experiments or other models, a is internal characteristic length, and 𝑙 is the external length. For a beam structure, the sizes of thickness and width are much smaller than the size of length. Therefore, the integral constitutive relations Eq.(2) can be represented in an equivalent differential form as, (1 − 𝜏 2 𝑙 2 ∇2 )𝜎 = 𝑇

(5)

4

ACCEPTED MANUSCRIPT By adding the elastic medium terms Eq.(3) and using Eq. (5) to incorporate the nonlocal effect, the governing equation given by Eq. (2) for the nonlocal Euler–Bernoulli beam model is obtained as 𝜕4 𝑤

𝜕2 𝑤

−(𝜇𝐾𝑠 + 𝐸𝐼) 𝜕𝑥 4 + (𝜇𝐾𝑤 + 𝐾𝑠 ) 𝜕𝑡 2 − 𝐾𝑤 𝑤 = 𝜌𝐴

𝜕2 𝑤 𝜕𝑡 2

𝜕4 𝑤

𝜕6 𝑤

− (𝜇𝜌𝐴 + 𝜌𝐼) 𝜕𝑥 2 𝑡 2 + 𝜇𝜌𝐼 𝜕𝑥 4 𝑡 2 (6)

where μ = e20 a2 is the nonlocal parameter used in the nonlocal continuum model. 3. Basic Formulations and Numerical Procedure

CR IP T

The transverse displacement of the nanobeams can be assumed in the following generalized form: 𝑤 = 𝑊(𝑥)𝑒 𝑖𝜔𝑡

(7)

By substituting eq. (7) into eq. (6) we get rid of time and that yields

𝜌𝐼)𝜔2

𝑑2 𝑊(𝑥) 𝑑𝑥 2

𝑑4 𝑊(𝑥) 𝑑𝑥 4

+ (𝜇𝐾𝑤 + 𝐾𝑠 )

− 𝜇𝜌𝐼𝜔2

𝑑2 𝑊(𝑥) 𝑑𝑥 2

𝑑4 𝑊(𝑥) 𝑑𝑥 4

3.1 The sixth order compact scheme

− 𝐾𝑤 𝑊(𝑥) = −𝜌𝐴𝜔2 𝑊(𝑥) + (𝜇𝜌𝐴 +

AN US

−(𝜇𝐾𝑠 + 𝐸𝐼)

(8)

The general approach for developing high-order compact differencing schemes is to utilize

M

the governing differential equation to help approximate truncation error terms. Compact schemes have been developed for a variety of differential equations, see e.g., [15-19] where

ED

supporting numerical studies showing the higher-order rates of convergence were presented. In general, a direct discretization of fourth-order partial derivative uses at least five grid

points [20-25].

PT

points. However, the fourth-order partial derivative can be also discretized using three grid

CE

In this Section, we present the details of deriving a sixth order finite difference scheme for Eq.(8). Let the beam shown in Fig.(1) be divided uniformly into 𝑁 + 1 intervals by the nodes 𝐿

𝑥𝑖 = 𝑖ℎ, 𝑖 = 0, 1, 2, . . . , 𝑁 + 1, where ℎ = 𝑁+1 . These nodes include the interior ones

AC

𝑥𝑖 , 𝑖 = 1, 2, . . . , 𝑁 and the end nodes 𝑥0 = 0 and 𝑥𝑁+1 = 𝐿. It is easy to use Taylor's expansion to derive the following symmetric finite difference approximations to 6-order accuracy for 𝑑2 𝑊(𝑥) 𝑑𝑥 2 𝑑4 𝑊(𝑥) 𝑑𝑥 4

1

= 𝛿𝑥2 𝑊𝑖 + 90 ℎ4 1

= 𝛿𝑥4 𝑊𝑖 − 6 ℎ2

𝑑6 𝑊(𝑥) 𝑑𝑥 6

𝑑6 𝑊(𝑥) 𝑑𝑥 6

𝑑2 𝑤 𝑑4 𝑤 𝑑𝑥 2

, 𝑑𝑥 4 at an interior node i , (2 ≤ 𝑖 ≤ 𝑁 − 1),

+ 𝑂(ℎ6 ), 1

− 80 ℎ4

𝑑8 𝑊(𝑥) 𝑑𝑥 8

(9) + 𝑂(ℎ6 )

5

(10)

ACCEPTED MANUSCRIPT where 𝛿𝑥2 𝑊𝑖 =

−𝑊𝑖−2 +16𝑊𝑖−1 −30𝑊𝑖 +16𝑊𝑖+1 −𝑊𝑖+2

𝛿𝑥4 𝑊𝑖 =

𝑊𝑖−2 −4𝑊𝑖−1 +6𝑊𝑖 −4𝑊𝑖+1 +𝑊𝑖+2

(11)

12ℎ2

(12)

ℎ4

In view of Eqs.(9-12), eq. (8) may be approximated at point 𝑖 as

CR IP T

(𝜇𝜌𝐼𝜔2 − (𝜇𝐾𝑠 + 𝐸𝐼))𝛿𝑥4 𝑊𝑖 + ((𝜇𝐾𝑤 + 𝐾𝑠 ) − (𝜇𝜌𝐴 + 𝜌𝐼)𝜔2 )𝛿𝑥2 𝑊𝑖 + (𝜌𝐴𝜔2 − 𝐾𝑤 )𝑊𝑖 + 𝜏𝑖 = 0 (13) where 𝜏𝑖 is the truncation error and is as follows 1

𝜏𝑖 = − 80 ℎ4 (𝜇𝜌𝐼𝜔2 − (𝜇𝐾𝑠 + 𝐸𝐼)) ( 1 6

2

(𝜇𝜌𝐼𝜔 − (𝜇𝐾𝑠 + 𝐸𝐼))} (

𝑑6 𝑊(𝑥) 𝑑𝑥 6

𝑑8 𝑊(𝑥) 𝑑𝑥 8

1

) + ℎ2 {90 ((𝜇𝐾𝑤 + 𝐾𝑠 ) − (𝜇𝜌𝐴 + 𝜌𝐼)𝜔2 ) − 𝑖

6

) + 𝑂(ℎ ) 𝑖

(14)

(

𝑑6 𝑊(𝑥) 𝑑𝑥 6

𝑑8 𝑊(𝑥) 𝑑𝑥 8

) =

((𝜇𝜌𝐴+𝜌𝐼)𝜔 2 −(𝜇𝐾𝑤 +𝐾𝑠 ))𝛿𝑥4 𝑊𝑖 −(𝜌𝐴𝜔2 −𝐾𝑤 )𝛿𝑥2 𝑊𝑖 𝜇𝜌𝐼𝜔 2 −(𝜇𝐾𝑠 +𝐸𝐼)

𝑖

) =

(15)

𝑑6 𝑊(𝑥) ) −(𝜌𝐴𝜔 2 −𝐾𝑤 )𝛿𝑥4 𝑊𝑖 𝑑𝑥6 𝑖

((𝜇𝜌𝐴+𝜌𝐼)𝜔 2 −(𝜇𝐾𝑤 +𝐾𝑠 ))(

(16)

M

(

AN US

By differentiating Eq.(8) w.r.t 𝑥, twice to get the 6th derivative of W and four times to get the 8th derivative, hence

𝜇𝜌𝐼𝜔 2 −(𝜇𝐾𝑠 +𝐸𝐼)

𝑖

ED

Substituting with Eqs. (15, 16) into eq. (14) then inserting the result into eq. (13) yields (𝜇𝜌𝐼𝜔2 − (𝜇𝐾𝑠 + 𝐸𝐼))𝛿𝑥4 𝑊𝑖 + ((𝜇𝐾𝑤 + 𝐾𝑠 ) − (𝜇𝜌𝐴 + 𝜌𝐼)𝜔2 )𝛿𝑥2 𝑊𝑖 + (𝜌𝐴𝜔2 − 𝐾𝑤 )𝑊𝑖 + − 𝐾𝑤 )𝛿𝑥4 𝑊𝑖 + ℎ2

((𝜇𝜌𝐴+𝜌𝐼)𝜔2 −(𝜇𝐾𝑤 +𝐾𝑠 ))𝛿𝑥4 𝑊𝑖 −(𝜌𝐴𝜔2 −𝐾𝑤 )𝛿𝑥2 𝑊𝑖

PT

1 4 ℎ (𝜌𝐴𝜔2 80

1

(𝜇𝜌𝐴 + 𝜌𝐼)𝜔2 ) − (𝜇𝜌𝐼𝜔2 − (𝜇𝐾𝑠 + 6

𝜇𝜌𝐼𝜔2 −(𝜇𝐾𝑠 +𝐸𝐼) 1 𝐸𝐼)) − ℎ2 ((𝜇𝜌𝐴 + 80

1

{90 ((𝜇𝐾𝑤 + 𝐾𝑠 ) −

𝜌𝐼)𝜔2 − (𝜇𝐾𝑤 + 𝐾𝑠 ))} = 0

(17)

AC

CE

Now, substituting for 𝛿𝑥2 𝑊𝑖 and 𝛿𝑥4 𝑊𝑖 from Eqs.(11,12) into Eq.(17) and due to the symmetry of these central difference schemes, one can conclude that the result would also be symmetric in the sense that it has identical coefficients for the pairs Wi+1 , Wi−1and Wi+2 , Wi−2 . Rearranging the result as a polynomial in 𝜔2 yields (𝑎2 𝑊𝑖−2 + 𝑎1 𝑊𝑖−1 + 𝑎0 𝑊𝑖 + 𝑎1 𝑊𝑖+1 + 𝑎2 𝑊𝑖+2 ) 𝜔4 + (𝑏2 𝑊𝑖−2 + 𝑏1 𝑊𝑖−1 + 𝑏0 𝑊𝑖 + 𝑏1 𝑊𝑖+1 + 𝑏2 𝑊𝑖+2 ) 𝜔2 + (𝑐2 𝑊𝑖−2 + 𝑐1 𝑊𝑖−1 + 𝑐0 𝑊𝑖 + 𝑐1 𝑊𝑖+1 + 𝑐2 𝑊𝑖+2 ) = 0 (18) The values of coefficients 𝑎0 , 𝑎1 , 𝑎2 , 𝑏0 , 𝑏1 , 𝑏2 , 𝑐0 , 𝑐1 , and 𝑐2 of Eq.(18) are given in appendix A. In fact Eq.(18) can be used as a 6th order difference equation only at regular nodes (2 ≤ 𝑖 ≤ 𝑁 − 1). However, for nodes 𝑖 = 1 and 𝑖 = 𝑁, fictitious nodes have to be added and boundary conditions are applied at end nodes 𝑖 = 0 and 𝑖 = 𝑁 + 1. Details of deriving 6th order 6

ACCEPTED MANUSCRIPT schemes at nodes 𝑖 = 1 and 𝑖 = 𝑁 for different boundary conditions are given in the next section. 3.2 The algebraic eigenvalue problem After applying the boundary conditions, one obtains a homogeneous system of equation, which corresponds to a generalized eigenvalue problem: (𝐴𝜔4 + 𝐵𝜔2 + 𝐶){𝑊} = 0

(19)

𝐴(𝑖, 𝑖 − 2) = 𝐴(𝑖, 𝑖 + 2) = 𝑎2 𝐵(𝑖, 𝑖 − 2) = 𝐵(𝑖, 𝑖 + 2) = 𝑏2 𝐶(𝑖, 𝑖 − 2) = 𝐶(𝑖, 𝑖 + 2) = 𝑐2

AN US

𝐴(𝑖, 𝑖) = 𝑎0 , 𝐴(𝑖, 𝑖 − 1) = 𝐴(𝑖, 𝑖 + 1) = 𝑎1 , 𝐵(𝑖, 𝑖) = 𝑏0 , 𝐵(𝑖, 𝑖 − 1) = 𝐵(𝑖, 𝑖 + 1) = 𝑏1 , 𝐶(𝑖, 𝑖) = 𝑐0 , 𝐶(𝑖, 𝑖 − 1) = 𝐶(𝑖, 𝑖 + 1) = 𝑐1 ,

CR IP T

where {𝑊} = [𝑊1 𝑊2 ⋯ 𝑊𝑀−1 𝑊𝑀 ]𝑇 is the vector of unknowns and 𝐴, 𝐵, and 𝐶 are 𝑀 × 𝑀 matrices. In fact, 𝑀 = 𝑁 for simply supported-simply simply-supported and clamped-clamped beams since 𝑊 = 0 at both ends. However, 𝑀 = 𝑁 + 1 for a clampedfree beam. According to Eq.(18), the nonzero elements in row 𝑖 , 2 ≤ 𝑖 ≤ 𝑁 − 1 of matrix A, B and C are given by:

4. Boundary conditions

(20)

M

In this section, 6th order difference equations are derived for the three commonly used boundary conditions, namely simply supported, clamped and free. In Fig.(2), the left end of the beam is shown with a fictitious node at 𝑥 = −ℎ. W0

W1

W2

W3

W4

x=-h

0

h

2h

3h

4h

ED

W-1

PT

Figure (2): real and fictitious nodes at left end of a beam.

CE

Now, coefficient of 𝜔4 in Eq.(18) can be written for node 𝑖 = 1 (𝑥 = ℎ) as 𝑎2 𝑊−1 + 𝑎1 𝑊0 + 𝑎0 𝑊1 + 𝑎1 𝑊2 + 𝑎2 𝑊3 ,

(21)

AC

with similar coefficients for 𝜔2 and 𝜔0 . In the next subsections, the relation between 𝑊−1and other unknowns are derived and appropriate coefficients in matrices 𝐴, 𝐵, and 𝐶 corresponding to irregular nodes are computed for different boundary conditions.

4.1

simply supported boundary condition

At the boundary point (𝑖 = 0), 𝑊0 = 0 ,

𝜕2 𝑊

| = 0 then using Taylor's expansion, it is easy

𝜕𝑥 2 0

to prove that

7

ACCEPTED MANUSCRIPT 𝑊1 + 𝑊−1 = 2𝑊0 + ℎ2

𝜕2 𝑊 | 𝜕𝑥 2 0

ℎ 4 𝜕4 𝑊 | 𝜕𝑥 4 0

+ 12

ℎ 6 𝜕6 𝑊

ℎ 4 𝜕4 𝑊

ℎ 6 𝜕6 𝑊

+ 360 𝜕𝑥 6 | + ⋯ = 12 𝜕𝑥 4 | + 360 𝜕𝑥 6 | + ⋯ 0

0

0

Our approach here is to get 6th order difference scheme by making use of the governing 𝜕2 𝑊

equation at this end point. Substituting 𝑊0 = 0 , 𝜕𝑥 2 | = 0 in Eq. (8) yields 0

important to notice that if Eq. (8) is differentiated twice, one can prove that similarly all higher even derivatives vanish for Accordingly,

𝜕4 𝑊

| = 0. It is

𝜕𝑥 4 0 𝜕6 𝑊

| = 0 and

𝜕𝑥 6 0

simply supported boundary condition.

CR IP T

𝑊−1 = −𝑊1 ,

(22)

produces not only 6th order scheme but an exact one. Substituting this value for 𝑊−1 and setting 𝑊0 = 0 in Eq.(21), then (𝑎0 − 𝑎2 )𝑊1 + 𝑎1 𝑊2 + 𝑎2 𝑊3 = 0

AN US

And hence, the nonzero coefficients in first row of matrix A can be set as 𝐴(1,1) = 𝑎0 − 𝑎2 , 𝐴(1,2) = 𝑎1 and 𝐴(1,3) = 𝑎2 .

(23)

The coefficients of the first row in matrices B, C are defined similarly. 4.2

Clamped boundary condition

scheme can be derived as

| = 0 then using Taylor's expansion a 6th order

𝜕𝑥 0

1

ED

5

𝜕𝑊

M

At the boundary point (𝑖 = 0), 𝑊0 = 0 ,

𝑤−1 = 10𝑤1 − 5𝑤2 + 3 𝑤3 − 4 𝑤4 + 𝒪(ℎ6 )

(24)

PT

Substituting 𝑤0 = 0, 𝑤−1from Eq.(24) in Eq.(21), the coefficient of 𝜔4 at node 1 is given by 8

1

CE

(10 𝑎2 + 𝑎0 )𝑊1 + (−5𝑎2 + 𝑎1 ) 𝑊2 + 𝑎2 𝑊3 − 𝑎2 𝑊4 + 𝒪(ℎ6 ) = 0, 3 4 which produces the following entries in matrix A 8

1

AC

𝐴(1,1) = 𝑎0 + 10𝑎2 , 𝐴(1,2) = (𝑎1 − 5𝑎2 ) , 𝐴(1,3) = 3 𝑎2 and 𝐴(1,4) = − 4 𝑎2

4.3

Free boundary condition

Consider a beam with free right end (𝑥 = 𝐿). At this end, 𝑊(𝐿) = 𝑊𝑁+1 is unknown but 𝜕2 𝑊 𝜕𝑥 2

𝜕3 𝑊

| = 0 , 𝜕𝑥 3 | = 0 and two fictitious nodes are required as shown in Fig.(3). 𝐿

𝐿

wN-2

wN-1 wN wN+3

𝑥 = (𝑁 − 2)ℎ (𝑁 − 1)ℎ 𝑁ℎ

wN+1

wN+2

𝐿

𝐿+ℎ

wN+3 𝐿 + 2ℎ

Figure (3): Two fictitious nodes at right free end 8

(25)

ACCEPTED MANUSCRIPT Using Taylor's expansion, the following 6th order finite difference schemes representing the free boundary conditions can be derived 𝑑2 𝑤

| = 𝑑𝑥 2

−𝑤𝑁+3 +16𝑤𝑁+2 −30𝑤𝑁+1 +16𝑤𝑁 −𝑤𝑁−1 12ℎ2

𝐿

𝑑3 𝑤

| = 𝑑𝑥 3

ℎ4 𝑑6 𝑤

+ 90 𝑑𝑥 6 | + 𝒪(ℎ6 ) = 0

−𝑤𝑁+3 −𝑤𝑁+2 +10𝑤𝑁+1 −14𝑤𝑁 +7𝑤𝑁−1 −𝑤𝑁−2 4ℎ3

𝐿

(26)

𝐿

+

ℎ3

𝑑6 𝑤

8

𝑑𝑥 6 𝐿

| + 𝒪(ℎ4 ) = 0

(27)

1

𝑤𝑁+2 = 17 (40 𝑤𝑁+1 − 30𝑤𝑁 + 8𝑤𝑁−1 − 𝑤𝑁−2 ) + 𝒪(ℎ6 ) 1

CR IP T

Solving Eqs.(26, 27) for 𝑤𝑁+2, 𝑤𝑁+3, yields

𝑤𝑁+3 = 17 (130 𝑤𝑁+1 − 208𝑤𝑁 + 111𝑤𝑁−1 − 16 𝑤𝑁−2 ) + 𝒪(ℎ6 ) At end Node N+1 (free)

AN US

𝑎2 𝑤𝑁−1 + 𝑎1 𝑤𝑁 + 𝑎0 𝑤𝑁+1 + 𝑎1 𝑤𝑁+2 + 𝑎2 𝑤𝑁+3 = 0

(28) (29)

(30)

Substituting from Eqs.(28,29), the nonzero elements in row (𝑁 + 1) of matrix A can be written as 𝐴(𝑁 + 1, 𝑁 − 2) =

1, 𝑁 − 1) =

128𝑎2 +8𝑎1 , 17

𝐴(𝑁 + 1, 𝑁) =

−208𝑎2 −13𝑎1 , 17

(31)

M

𝐴(𝑁 + 1, 𝑁 + 1) =

−16𝑎2−𝑎1 , 𝐴(𝑁 + 17 130𝑎2 +40𝑎1 +17𝑎0 17

At the node next to free end (Node N)

ED

𝑎2 𝑤𝑁−2 + 𝑎1 𝑤𝑁−1 + 𝑎0 𝑤𝑁 + 𝑎1 𝑤𝑁+1 + 𝑎2 𝑤𝑁+2 = 0

Substituting from Eqs.( 28,29), the nonzero elements in row (𝑁 + 1) of matrix A is given by 17

, 𝐴(𝑁, 𝑁 − 1) =

PT

16𝑎2

8𝑎2 +17𝑎1 17

, 𝐴(𝑁, 𝑁) =

CE

𝐴(𝑁, 𝑁 − 2) =

−30𝑎2 +17𝑎0 17

, 𝐴(𝑁, 𝑁 + 1) =

40𝑎2 +17𝑎1 17

(32)

5. Numerical results

AC

To study the dynamical analysis of nanobeam, the polynomial eigen value problem Eq.(19) must be solved for 𝜔2 which has 2𝑀 solutions. As a first step, matrices A, B and C are formed by introducing appropriate coefficients using Eq.(20) corresponding to regular interior points and appropriate coefficients corresponding to specific boundary conditions (Eq.(23) for a simply supported end, Eq.(25) for a clamped end, and Eqs.(31,32) at a free end). However, solving the eigenvalue problem (Eq.(19)), the resulting eigen frequencies were inaccurate. To analyze Eq.(19), its coefficients (norms of matrices A, B and C) were computed showing large variation in norms values with ratios of order of 104 and 108 . Practically, in such situations, some scaling procedure is needed. We suggest a natural scale

9

ACCEPTED MANUSCRIPT

ρA

factor such that Eq.(19) is solved for the non-dimensional frequency ω ̅ = L2 √ EI × ω rather than and ω , hence Eq.(19) is replaced by (𝐴𝜔 ̅ 4 + 𝐵̅ 𝜔 ̅ 2 + 𝐶̅ ){𝑊} = 0

(33)

2

4

𝜌𝐴 𝜌𝐴 where 𝐵̅ = (𝐿2 √ 𝐸𝐼 ) 𝐵 and 𝐶̅ = (𝐿2 √ 𝐸𝐼 ) 𝐶.

CR IP T

To demonstrate the accuracy of the proposed sixth-order discretization scheme, the polynomial eigen value problem Eq.(33) is solved for 𝜔 ̅ which has 4𝑀 eigenvalues. In general the most important are the fundamental ones corresponding to the few smallest eigen values. The three cases of boundary conditions are considered and numerical results have been presented. The following properties for the nanobeams considered are: 𝐸 = 30 × ℎ4

106 , 𝜌 = 1 and 𝑏 = ℎ, 𝐼 = 12 , 𝐴 = ℎ2 .

AN US

For the purpose of comparison with Ansari et al.[13], we have to choose ℎ = 1 and 𝐿 varies in the range 20–50. However, for the purpose of comparison with Eltaher et al.[14], 𝐿 = 10 and ℎ varies in the range 0.1–0.5

4.1.1 Exact solution

ED

4.1 simply supported nanobeam

M

In case of simply supported boundary condition, the results obtained from the proposed 6th order method are also compared with those from the exact solution.

PT

In case of simply supported nanobeam, the boundary conditions are

𝑑2 𝑊 𝑑𝑥 2

= 0, 𝑊 = 0 at both

CE

ends and hence an exact solution for transverse displacement 𝑊of the nanobeam can be assumed in the form 𝑚𝜋𝑥

𝑊(𝑥, 𝑡) = ∑∞ 𝑚=1 𝑊𝑚 sin (

𝐿

) 𝑒 𝑖𝜔𝑚𝑡

(34)

AC

where 𝜔𝑚 is the value of the mth frequency mode. Substituting Eq. (34) into Eq. (8) and solving for 𝜔𝑚 , an explicit formula for the fundamental frequencies of a simply supported nanobeam is obtained as 2 𝜔𝑚 =

(𝑚𝜋)4 𝜇𝐾𝑠 +(𝑚𝜋)4 𝐸𝐼+(𝑚𝜋)2 𝜇𝐾𝑤 𝐿2 +(𝑚𝜋)2 𝜇𝐾𝑠 𝐿2 +𝐾𝑤 𝐿4 𝜌((𝑚𝜋)2 𝐿2 𝐼+(𝑚𝜋)2 𝐴𝜇𝐿2 +𝐴𝐿4 +(𝑚𝜋)4 𝐼𝜇)

(35)

To demonstrate accuracy and stability of the present compact finite difference scheme, different numbers of interior nodes (𝑁 = 5, 8, 16) are used in the discretization of the nanobeam. The results are compared with Ansari et.al. [13] and the exact values given by Eq.(35) in Table (1). It is clear that the present method converges more rapidly to the exact values than [13]. This is expected because 6th order formulas for both second and fourth order 10

ACCEPTED MANUSCRIPT derivatives (Eqs.(9,10)) are introduced in the discretization of the governing equation. Precisely, the third term in the right hand side of Eq.(10) is responsible to raise the accuracy compared with [13]. Additionally, the derived discretization formulas Eqs.(22-23) represent the simply supported boundary conditions exactly. The present results coincide (to 4 decimal places at least) with the exact solution when 𝑁 = 16. In addition, accurate results (to 2 decimal places) are obtained even with small number of nodes (𝑁 = 5). Table 1: Comparison of non-dimensional fundamental frequencies 𝜔 ̅1 of present method, Ansari et.al. [13] and exact solution for simply supported nanobeam (𝐿 = 20, ℎ = 1). 𝐾𝑤 = 𝐾𝑠 = 0

20 0 1 2 3 4

present [13] present [13] present [13] present [13] present [13]

𝐾𝑤 = 20, 𝐾𝑠 = 2

N:5

8

16

Exact

9.8594 9.9333 9.7397 9.8130 9.6243 9.6969 9.5129 9.5849 9.4053 9.4767

9.8595 9.8818 9.7400 9.7623 9.6247 9.6470 9.5135 9.5357 9.4060 9.4282

9.8595 9.8628 9.7400 9.7436 9.6248 9.6286 9.5136 9.5175 9.4062 9.4103

N:5

9.8595 9.7400 9.6248

9.9241 9.9975 9.8053 9.8780 9.6906 9.7628 9.5800 9.6515 9.4731 9.5440

CR IP T

𝜇

8

16

9.9242 9.9464 9.8055 9.8277 9.6910 9.7132 9.5806 9.6027 9.4739 9.4959

9.9242 9.9275 9.8056 9.8091 9.6911 9.6948 9.5807 9.5846 9.4740 9.4781

AN US

𝐿/ℎ

9.5136 9.4062

Exact

9.9242 9.8056 9.6911 9.5807 9.4740

ED

M

In Table 2, the first five fundamental eigenvalues computed with different number of interior grid points 'N' are reported and compared with the exact ones. Results show that the present method converges rapidly to exact solution for different frequency modes even if 𝐾𝑤 and 𝐾𝑠 have non-zero values. It is noted also that values of foundation modulus 𝐾𝑤 and 𝐾𝑠 have little effect on the fundamental frequencies. In addition, smaller mesh sizes (greater “N”) is needed to compute higher frequency modes accurately.

PT

Table 2: Comparison of the fundamental five modes of non-dimensional frequencies of present method and exact solution for simply supported nanobeam at 𝜇 = 2, 𝐿 = 10, ℎ = 1. 𝜔 ̅

CE

𝐿/ℎ

AC

10

𝐾𝑤 = 0,𝐾𝑠 = 0

N:5

8

𝜔 ̅1 𝜔 ̅2 𝜔 ̅3 𝜔 ̅4 𝜔 ̅5

8.9809 28.8035 48.7856 62.6479 69.8399

8.9823 28.9929 50.8870 70.2060 85.0730

𝜔 ̅1 𝜔 ̅2 𝜔 ̅3 𝜔 ̅4 𝜔 ̅5

8.9854 28.8049 48.7865 62.6486 69.8405

8.9868 28.9942 50.8878 70.2065 85.0735

16

20

8.9826 8.9826 29.0339 29.0359 51.3955 51.4199 72.5838 72.7061 91.5706 91.9526 𝐾𝑤 = 20, 𝐾𝑠 = 2 8.9871 29.0353 51.3963 72.5843 91.5711

11

8.9871 29.0373 51.4207 72.7066 91.9531

30

Exact

8.9826 29.0370 51.4342 72.7785 92.1817

8.9826 29.0373 51.4381 72.7977 92.2427

8.9871 29.0384 51.4350 72.7791 92.1821

8.9871 29.0387 51.4388 72.7983 92.2431

ACCEPTED MANUSCRIPT

CR IP T

Next, the effect of nonlocal parameter 𝜇 and slenderness ratio 𝐿/ℎ on the free (𝐾𝑠 = 𝐾𝑤 = 0) vibration characteristics for simply supported beam is considered. In Table 3, the first five fundamental modes are summarized and compared with the exact ones and those given by Eltaher et al. [14] where an efficient finite element procedure was used. Although the computations in both the present 6th order method and [14] use the same number of unknowns (𝑁 = 50), the present method is more efficient. It is observed that as the nonlocal parameter increased the frequencies decreased. It is also observed from Table 3 and other results not reported here that as 𝐿/ℎ increased the frequency increased slightly. The present results are in good agreement with the exact solution. Table 3: Effects of (𝐿/ℎ) and nonlocal parameter on the fundamental frequencies for simply supported beam (𝐿 = 10, ℎ = 0.5, 0.1).

1

2

3

4

100

PT

5

CE

0

AC

1

2

3

4

5

𝜔 ̅1 9.8595 9.8595 9.8798 9.4062 9.4062 9.4238 9.0102 9.0102 9.0257 8.6604 8.6604 8.6741 8.3483 8.3483 8.3606 8.0678 8.0678 8.0789 9.8692 9.8692 9.87 9.4155 9.4155 9.4162 9.0191 9.0191 9.0197 8.6689 8.6689 8.6695 8.3566 8.3566 8.3571 8.0757 8.0757 8.0762

𝜔 ̅2 39.3171 39.3171 39.646 33.2910 33.2911 33.4887 29.3905 29.3905 29.5252 26.6023 26.6023 26.7012 24.4818 24.4818 24.5582 22.7990 22.7990 22.86 39.4719 39.4719 39.4849 33.4222 33.4222 33.4301 29.5062 29.5063 29.5117 26.7071 26.7071 26.7111 24.5782 24.5783 24.5814 22.8888 22.8888 22.8914

𝜔 ̅3 88.0158 88.0158 89.7046 64.0510 64.0515 64.6783 52.8208 52.8213 53.1629 45.9760 45.9765 46.1957 41.2481 41.2486 41.4027 37.7314 37.7319 37.8467 88.7936 88.7936 88.8595 64.6170 64.6175 64.6429 53.2876 53.2881 53.3024 46.3822 46.3828 46.3922 41.6126 41.6131 41.6199 38.0648 38.0653 38.0705

𝜔 ̅4 155.3782 155.3785 160.904 96.7476 96.7506 97.9538 76.1935 76.1964 76.748 64.8657 64.8684 65.1889 57.4406 57.4431 57.652 52.0933 52.0956 52.2422 157.8099 157.8099 158.0185 98.2615 98.2646 98.3149 77.3857 77.3887 77.4133 65.8806 65.8834 65.8986 58.3394 58.3420 58.3525 52.9084 52.9108 52.9187

AN US

0

present exact [14] present exact [14] present exact [14] present exact [14] present exact [14] present exact [14] present exact [14] present exact [14] present exact [14] present exact [14] present exact [14] present exact [14]

M

𝜇

ED

𝐿/ℎ 20

12

𝜔 ̅5 240.6313 240.6328 217.0355 129.2157 129.2269 131.0051 98.7662 98.7761 99.4839 83.0066 83.0153 83.3891 72.9795 72.9874 73.2108 65.8837 65.8909 66.0331 246.4867 246.4868 246.9973 132.3588 132.3706 132.4498 101.1685 101.1790 101.2145 85.0255 85.0349 85.056 74.7545 74.7630 74.7775 67.4861 67.4939 67.5046

ACCEPTED MANUSCRIPT 4.2

Clamped-clamped nanobeam

Since there is no exact solution for fundamental frequencies of a clamped-clamped nanobeam, the results of the present 6th order method will be compared with existing works. Table 4 shows the convergence of the present method in comparison with Ansari et al. [13] where the fundamental frequency ω ̅ 1 is computed on different numbers of grid points "N". Good agreement is noticed and the rapid convergence of the present method is also observed.

𝜇

𝐾𝑤 = 0, 𝐾𝑠 = 0

20 0 1 2 3 4

present [13] present [13] present [13] present [13] present [13]

𝐾𝑤 = 20, 𝐾𝑠 = 2

N:5

8

16

20

30

N:5

22.4816 22.5368 22.1498 22.2074 21.8309 21.8909 21.5242 21.5865 21.2290 21.2935

22.3671 22.3778 22.0306 22.0440 21.7082 21.7240 21.3989 21.4171 21.1020 21.1223

22.3458 22.3466 22.0078 22.0114 21.6843 21.6904 21.3743 21.3827 21.0765 21.0873

22.3451 22.3454 22.0068 22.0102 21.6839 21.6891 21.3730 21.3813 21.0760 21.0858

22.3447 22.3448 22.0066 22.0095 21.6830 21.6884 21.3729 21.3805 21.0752 21.0851

22.5101 22.5652 22.1787 22.2362 21.8602 21.9201 21.5539 21.6161 21.2591 21.3236

8

16

20

30

22.3957 22.4064 22.0597 22.0730 21.7377 21.7535 21.4288 21.4470 21.1323 21.1526

22.3744 22.3753 22.0369 22.0405 21.7138 21.7200 21.4042 21.4126 21.1069 21.1177

22.3737 22.3741 22.0358 22.0392 21.7134 21.7186 21.4029 21.4112 21.1064 21.1162

22.3734 22.3734 22.0357 22.0386 21.7125 21.7179 21.4028 21.4105 21.1056 21.1154

AN US

𝐿/ℎ

CR IP T

Table 4: Comparison of non-dimensional fundamental frequencies 𝜔 ̅1 of present method and [13] for clamped-clamped nanobeam (𝐿 = 20, ℎ = 1).

ED

M

Next, the free vibration of the clamped-clamped nanobeam is considered and the five fundamental frequencies is reported in Table (5) with corresponding results of Altaher et al. [14] for different nonlocal parameter 𝜇 and 𝐿/ℎ ratios. The present results are in good agreement with the previous work [14]. As shown in table 5, as the nonlocal parameter increased the frequencies decreased.

PT

Table 5: Effects of (𝐿/ℎ) and nonlocal parameter on the fundamental frequencies for clamped-clamped beam (𝐿 = 10, ℎ = 0.5, 0.1). 𝜇 0

CE

𝐿/ℎ 20

AC

1 2 3 4 5

100

0 1 2

present [14] present [14] present [14] present [14] present [14] present [14] present [14] present [14] present

𝜔 ̅1 22.3447 22.4022 21.0751 21.1236 19.9954 20.0368 19.0632 19.099 18.2482 18.2795 17.5280 17.5556 22.3721 22.3744 21.1077 21.1096 20.0313

𝜔 ̅2 61.3790 61.9872 50.6788 51.0178 44.1031 44.3229 39.5507 39.7063 36.1630 36.2795 33.5166 33.6074 61.6609 61.6847 50.9707 50.9844 44.3827

13

𝜔 ̅3 119.6760 122.2778 84.6363 85.5011 69.1806 69.6261 60.0021 60.2766 53.7517 53.938 49.1426 49.2768 120.8531 120.9536 85.6707 85.7081 70.0819

𝜔 ̅4 196.3774 204.8657 118.7677 120.2577 93.0371 93.6814 79.0692 79.4257 69.9717 70.1931 63.4408 63.587 199.7147 200.0042 121.2332 121.3058 95.0531

𝜔 ̅5 290.6602 216.3147 151.7723 153.7997 115.7132 116.4608 97.2369 97.5988 85.5270 85.7194 77.2551 77.3578 298.2217 298.89 156.5060 156.6286 119.4346

ACCEPTED MANUSCRIPT [14] present [14] present [14] present [14]

3 4 5

20.033 19.1013 19.1028 18.2877 18.289 17.5685 17.5696

44.392 39.8152 39.822 36.4131 36.4184 33.7538 33.7581

70.1033 60.8098 60.8244 54.4905 54.5015 49.8281 49.8369

95.0923 80.8175 80.8443 71.5376 71.5581 64.8720 64.8888

119.5018 100.4097 100.4574 88.3429 88.3807 79.8148 79.8467

CR IP T

Since the analytical fundamental frequencies of the clamped-clamped beam is not available, the convergence rate of the developed numerical method can be computed by the following formula [23] ℎ

Convergence rate =

log(∆𝜔 ℎ ⁄∆𝜔 2 ) log(2)

(36)

AN US

where ∆𝜔 ℎ/2 = 𝜔 ℎ/2 − 𝜔ℎ denotes the difference in two computed eigenvalues, 𝜔ℎ and 𝜔 ℎ/2 using grid sizes ℎ and ℎ/2, respectively. The convergence rates for the fundamental frequencies of a clamped-clamped nanobeam are shown in Table 6. It is important to mention that the rates reported in Table 6 are computed for the solution of the generalized eigenvalue problem Eq. (19) but not for the differential equation Eq. (8) which is 6th order. The reported rates are affected by the presence of the numerical error in the solution for the eigenvalues.

4.2

𝜔 ̅2 43.82926 44.07373 44.10153 44.10323 44.10333

Rate

ED

Rate

3.78 4.35 4.26

PT

𝜔 ̅1 19.97708 19.99408 19.99532 19.99538 19.99538

3.14 4.03 4.10

𝜔 ̅3 67.41981 68.97837 69.16917 69.1818 69.18255

Rate

3.03 3.92 4.07

𝜔 ̅4 86.45322 92.29031 92.99234 93.04194 93.04495

Rate

3.06 3.82 4.04

𝜔 ̅5 98.69522 113.7438 115.5877 115.7269 115.7355

Rate

3.03 3.73 4.02

CE

Mesh size L/8 L/16 L/32 L/64 L/128

M

Table 6: Convergence rates using formula (36) of the fundamental frequencies for clampedclamped beam (𝐿 = 10, ℎ = 0.5, 𝜇 = 2, 𝐾𝑤 = 0, 𝐾𝑠 = 0).

Clamped-free nanobeam

AC

Again, as no exact fundamental frequencies is available for clamped-free nanobeams, we compare our results with Ansari et al. [13] for the fundamental frequency in Table (7) and Altaher et al. [14] for the first five free fundamental frequencies in Table (8). Table 7: Comparison of non-dimensional fundamental frequencies 𝜔 ̅1of present method and [13] for clamped-free nanobeam (𝐿 = 20, ℎ = 1). 𝜇 0 1

𝐾𝑤 = 0, 𝐾𝑠 = 0 present [13] present

𝐾𝑤 = 20, 𝐾𝑠 = 2

N:5

8

16

20

30

N:5

8

16

20

30

3.5164 3.4241 3.5202

3.5163 3.4792 3.5201

3.5163 3.5069 3.5201

3.5163 3.5103 3.5201

3.5163 3.5136 3.5201

3.6939 3.6061 3.6976

3.6939 3.6585 3.6975

3.6939 3.6849 3.6974

3.6939 3.6881 3.6974

3.6939 3.6913 3.6974

14

ACCEPTED MANUSCRIPT

2 3 4

[13] present [13] present [13] present [13]

3.4276 3.5241 3.4312 3.5280 3.4348 3.5319 3.4385

3.4829 3.5240 3.4866 3.5278 3.4903 3.5317 3.4941

3.5107 3.5239 3.5145 3.5277 3.5183 3.5316 3.5221

3.5141 3.5239 3.5179 3.5277 3.5217 3.5316 3.5255

3.5174 3.5239 3.5212 3.5277 3.5250 3.5316 3.5289

3.6095 3.7013 3.6129 3.7050 3.6164 3.7087 3.6198

3.6620 3.7011 3.6656 3.7048 3.6691 3.7085 3.6727

3.6885 3.7011 3.6921 3.7047 3.6957 3.7084 3.6993

3.6917 3.7011 3.6953 3.7047 3.6989 3.7084 3.7026

3.6949 3.7011 3.6985 3.7047 3.7022 3.7084 3.7058

2 3 4 5 100

0 1 2 3

PT

4

𝜔 ̅2 22.0040 22.1106 20.6445 20.7329 19.4718 19.5464 18.4453 18.5092 17.5359 17.5913 16.7223 16.7707 22.0333 22.0375 20.6783 20.6817 19.5084 19.5111 18.4834 18.4857 17.5747 17.5767 16.7611 16.7629

CE

5

𝜔 ̅1 3.5163 3.5177 3.5316 3.533 3.5472 3.5486 3.5632 3.5646 3.5796 3.581 3.5964 3.5978 3.5160 3.5161 3.5313 3.5314 3.5469 3.547 3.5630 3.563 3.5794 3.5795 3.5963 3.5963

𝜔 ̅3 61.4044 62.2043 50.7623 51.1979 44.2710 44.5479 39.8220 40.0146 36.5507 36.6923 34.0293 34.1377 61.6854 61.7171 51.0534 51.0695 44.5507 44.5603 40.0881 40.0945 36.8045 36.8089 34.2728 34.2759

𝜔 ̅4 119.6755 108.9708 84.6179 85.6564 69.0957 69.61.39 59.8153 60.1276 53.4381 53.645 48.6855 48.8302 120.8522 120.9744 85.6528 85.6894 69.9963 70.0139 60.6192 60.6294 54.1686 54.1751 49.3572 49.3616

AN US

1

present [14] present [14] present [14] present [14] present [14] present [14] present [14] present [14] present [14] present [14] present [14] present [14]

M

𝜇 0

ED

𝐿/ℎ 20

CR IP T

Table 8: Effects of (𝐿/ℎ) and nonlocal parameter 𝜇 on the fundamental frequencies for clamped-free beam (𝐿 = 10, ℎ = 0.5, 0.1). 𝜔 ̅5 196.3825 122.6045 118.8038 107.3971 93.1273 93.7548 79.2597 79.5971 70.2987 70.4918 63.9303 64.0432 199.7172 200.0512 121.2702 121.3306 95.1465 95.1713 81.0167 81.0295 71.8816 71.8887 65.3892 65.3932

AC

In Table (7), the fundamental frequency for different values of foundation constants are computed on grids with different number of interior nodes N. Our results are in good agreement with [13]. Moreover, it is observed that they converge faster than those reported in [13] with more accurate results on grids with smaller numbers of points. With respect to comparison with Altaher et al. [14], Table (8) summarizes the results of the first five frequency modes for different nonlocal parameter 𝜇 and 𝐿/ℎ ratios. Most of our results are close to those of [14] but deviate for higher frequency modes at low values of 𝜇 (𝜇 = 0, 1) especially for 𝐿/ℎ = 20. This feature is also noticed for other boundary conditions as reported in Tables (3 and 5). 6. Conclusions

15

ACCEPTED MANUSCRIPT Vibration analysis of a nonlocal Euler–Bernoulli beam embedded in an elastic medium is considered. The nonlocal effect of nanobeams is introduced using the nonlocal differential elasticity of Eringen. 6th order accurate finite difference discretization schemes are derived for each of the governing differential equation and different types of boundary conditions. Fundamental frequencies of the nanobeams computed by the proposed scheme are compared with those available in literature. Moreover, in case of simply supported nanobeams, comparisons with the exact ones are presented to demonstrate its accuracy and convergence. The following conclusions can be derived from the presented results.

AN US

CR IP T

1. The proposed 6th order scheme is simple and outperforms similar existing methods. 2. Type of boundary conditions has great influence on the values of fundamental frequencies. 3. As the nonlocal parameter increased, the values of fundamental frequencies decreased for all boundary conditions except for clamped-free beam where the first fundamental frequency increased slightly. 4. The Pasternak elastic foundation coefficients have little influence on the fundamental frequencies which increase slightly as these coefficients increase. Similar 6th order finite difference discretization schemes can be developed to analyze bending and buckling of nanobeams under various loading conditions.

M

Appendix A

ED

The components of the coefficients matrices A, B and C in Eq.(19) corresponding to the interior points can be expressed as follows. For matrix A 𝑎1 = −

16ℎ4

𝑎3 = −

PT

𝜌2 (17(𝜇𝐴+𝐼)𝐴ℎ6 +(51𝐼 2 +195𝜇𝐴𝐼+51𝜇 2 𝐴2 )ℎ4 −360(𝜇𝐴+𝐼)𝜇𝐼ℎ2 −2160𝜇 2 𝐼2 ) 6ℎ4

𝜌2 (17(𝜇𝐴+𝐼)𝐴ℎ6 +(204𝐼 2 +420𝜇𝐴𝐼+204𝜇 2 𝐴2 )ℎ4 +720(𝜇𝐴+𝐼)𝜇𝐼ℎ2 −8640𝜇2 𝐼 2 )

CE

𝑎2 =

𝜌2 (85(𝜇𝐴+𝐼)𝐴ℎ6 +(204𝐼 2 −540𝜇𝐴𝐼+204𝜇 2 𝐴2 )ℎ4 −2160(𝜇𝐴+𝐼)𝜇𝐼ℎ2 −8640𝜇 2 𝐼2 )

96ℎ4

(A1) (A2) (A3)

AC

For matrix B

135((𝜇𝐾𝑠 +𝐸𝐼)(𝜇𝐴+𝐼)+𝜇𝐼(𝜇𝐾𝑤 +𝐾𝑠 ))ℎ2 +1080(𝐸𝐼+𝜇𝐾𝑠 )𝜇𝐼

𝑏1 = −𝜌 { 𝜌 4

2

8

60((𝜇𝐾𝑠 +𝐸𝐼)(𝜇𝐴+𝐼)+𝜇𝐼(𝜇𝐾𝑤 +𝐾𝑠 ))ℎ2 +720(𝐸𝐼+𝜇𝐾𝑠 )𝜇𝐼 ℎ4

}−

(A4) 17 6

{65𝜇(𝐴𝐾𝑠 + 𝐼𝐾𝑤 ) + 34(𝜇 2 𝐴𝐾𝑤 + 𝐼𝐾𝑠 ) + 31𝐸𝐼𝐴}

𝑏3 = 𝜌 { 𝜌

85

} + 16 𝜌{2𝜇𝐴𝐾𝑤 + 𝐴𝐾𝑠 + 𝐼𝐾𝑤 }ℎ2 −

{135𝜇(𝐴𝐾𝑠 + 𝐼𝐾𝑤 ) − 102(𝜇 2 𝐴𝐾𝑤 + 𝐼𝐾𝑠 ) + 237𝐸𝐼𝐴}

𝑏2 = 𝜌 { 𝜌

ℎ4

15((𝜇𝐾𝑠 +𝐸𝐼)(𝜇𝐴+𝐼)+𝜇𝐼(𝜇𝐾𝑤 +𝐾𝑠 ))ℎ2 −360(𝐸𝐼+𝜇𝐾𝑠 )𝜇𝐼 2ℎ4

{35𝜇(𝐴𝐾𝑠 + 𝐼𝐾𝑤 ) + 34(𝜇 2 𝐴𝐾𝑤 + 𝐼𝐾𝑠 ) + 𝐸𝐼𝐴} 16

𝜌{2𝜇𝐴𝐾𝑤 + 𝐴𝐾𝑠 + 𝐼𝐾𝑤 }ℎ2 − (A5)

17

} + 96 𝜌{2𝜇𝐴𝐾𝑤 + 𝐴𝐾𝑠 + 𝐼𝐾𝑤 }ℎ2 + (A6)

ACCEPTED MANUSCRIPT For matrix C 𝑐1 = 16

ℎ4

(𝐾𝑠 𝐾𝑤 + 𝜇𝐾𝑤2 )ℎ2

𝑐2 = − 17 6

96

(A7)

60((𝐸𝐼+𝜇𝐾𝑠 )(𝐾𝑠 +𝜇𝐾𝑤 )ℎ2 +6(𝐸𝐼+𝜇𝐾𝑠 )2 ) ℎ4

1

+ 6 {51(𝐾𝑠2 + 𝜇 2 𝐾𝑤2 ) + 195𝜇𝐾𝑤 𝐾𝑠 + 93𝐾𝑤 𝐸𝐼} +

(𝐾𝑠 𝐾𝑤 + 𝜇𝐾𝑤2 )ℎ2

𝑐3 = − 17

1

+ 4 {−51(𝐾𝑠2 + 𝜇 2 𝐾𝑤2 ) + 135𝜇𝐾𝑤 𝐾𝑠 + 237𝐾𝑤 𝐸𝐼} −

(A8)

15((𝐸𝐼+𝜇𝐾𝑠 )(𝐾𝑠 +𝜇𝐾𝑤 )ℎ2 −12(𝐸𝐼+𝜇𝐾𝑠 )2 ) 2ℎ4

1

− 8 {17(𝐾𝑠2 + 𝜇 2 𝐾𝑤2 ) + 35𝜇𝐾𝑤 𝐾𝑠 + 𝐾𝑤 𝐸𝐼} −

(𝐾𝑠 𝐾𝑤 + 𝜇𝐾𝑤2 )ℎ2

References

(A9)

W. D. Zhang, Y. Wen, S. M. Liu, W. C. Tjiu, G. Q. Xu, L. M. Gan, Synthesis of

AN US

[1]

CR IP T

85

135((𝐸𝐼+𝜇𝐾𝑠 )(𝐾𝑠 +𝜇𝐾𝑤 )ℎ2 +4(𝐸𝐼+𝜇𝐾𝑠 )2 )

vertically aligned carbon nanotubes on metal deposited quartz plates, Carbon 40 (2002) 1981–1989. [2]

L. Liu, Y. Zhang, Multi-wall carbon nanotube as a new infrared detected material, Sensors and Actuators A, 116 (2004) 394–397.

X. B. Yan, X. J. Chen, B. K. Tay, K. A. Khor, Transparent and flexible glucose via

M

[3]

layer-by-layer assembly of multi-wall carbon nanotubes and glucose oxidase,

[4]

ED

Electrochemistry Communications 9 (2007) 1269–1275. C. Zhao, Y. Song, J. Ren, X. Qu, A DNA nanomachine induced by single-walled

[5]

PT

carbon nanotubes on gold surface, Biomaterials 30 (2009) 1739–1745. C. Qin, J. Shen, Y. Hu, M. Ye, Facile, attachment of magnetic nanoparticles to carbon

CE

nanotubes via robust linkages and its fabrication of magnetic nanocomposites, Composites Science and Technology 69 (2009) 427–431. R. Ansari, S. Sahmani, Bending behavior and buckling of nanobeams including surface

AC

[6]

stress effects corresponding to different beam theories, International Journal of Engineering Science 49 (2011) 1244–1255.

[7]

A. C. Eringen, Nonlocal polar elastic continua, Int. J. Eng. Sci. 10 (1972) 1–16.

[8]

A. C. Eringen, D. Edelen, On nonlocal elasticity, Int. J. Eng. Sci. 10 (3) (1972) 233– 248.

[9]

A. C. Eringen, On differential equations of nonlocal elasticity and solutions of screw dislocation and surface waves, J. Appl. Phys. 54 (1983) 4703–4710.

17

ACCEPTED MANUSCRIPT [10]

A. C. Eringen, Nonlocal Continuum Field Theories, Springer-Verlag, New York, 2002.

[11] J. N. Reddy, Nonlocal theories for bending, buckling and vibration of beams, International Journal of Engineering Science 45 (2007) 288–307. [12] M. Aydogu, A general nonlocal beam theory: its application to nanobeam bending, buckling and vibration, Phys. E 41 (2009) 1651–1655. [13] R. Ansari, R. Gholami, K. Hosseini, S. Sahmani, A sixth-order compact finite

CR IP T

difference method for vibrational analysis of nanobeams embedded in an elastic medium based on nonlocal beam theory, Mathematical and Computer Modelling 54 (2011) 2577–2586. [14]

M. A. Eltaher, A. E. Alshorbagy, F.F. Mahmoud, Applied Mathematical Modelling 37 (2013) 4787–4797.

AN US

[15] M. M. Gupta, A fourth-order Poisson solver, J. Comput. Phys., 55 (1984) 166-172. [16] W. F. Spotz and G. F. Carey, Higher order compact scheme for the steady stream function vorticity equations. Int J Numer Meth Eng, 38 (1995) 3479-3512. [17] M. gupta, J. zhang, High accuracy multigrid solution of the 3D Convection Diffusion

M

equation, Applied Mathematics and Computation 113 (2000) 249–274. [18] Z.F. Tian , S.Q. Dai, High-order compact exponential finite difference methods for

ED

convection-diffusion type problems, Journal of Computational Physics, 220 (2007) 952–974.

[19] Y. Fu, Compact fourth-order finite difference schemes for Helmholtz equation with

PT

high wave numbers, J. Comput. Math., 26 (2008) 98-111. [20] J. W. Stephenson, Single cell discretizations of order two and four for biharmonic

CE

problems, Journal of Computational Physics, 55 (1984), 65-80. [21] N. Mohamed, S.A. Mohamed, L.F. Seddek, Exponential higher-order compact scheme

AC

for 3D steady convection–diffusion problem, Applied Mathematics and Computation, 232 (2014), 1046–1061.

[22] M. Ben-Artzi, J. Croisille, D. Fishelov, S. Trachtenberg, A pure-compact scheme for the streamfunction formulation of Navier-Stokes equations, Journal of Computational Physics, 205 (2005), 640-664. [23] Z. F. Tian, P. X. Yu, An efficient compact difference scheme for solving the streamfunction formulation of the incompressible Navier-Stokes equations Journal of Computational Physics, 230 (2011), 6404-6419.

18

ACCEPTED MANUSCRIPT [24] P. X. Yu, Z. F. Tian, Compact computations based on a stream function-velocity formulation of two-dimensional steady laminar natural convection in a square, Physical review E, 85 (2012), 036703-12. [25] P. X. Yu, Z. F. Tian, A compact stream function-velocity scheme on nonuniform grids for the 2D steady incompressible Navier-Stokes equations, Computers and

AC

CE

PT

ED

M

AN US

CR IP T

Mathematics with Applications, 66 (2013), 1192-1212.

19