# Nanobeams and AFM Subject to Piezoelectric and Surface Scale Effects

Authors: Julio Claeyssen; Rosemaira Copetti; Leticia Tonetto; Joao Carvalho
Journal: Advances in Mathematical Physics
 (2018)

---

## Abstract

Vibration dynamics of elastic beams that are used in nanotechnology, such as atomic force microscope modeling and carbon nanotubes, are considered in terms of a fundamental response within a matrix framework. The modeling equations with piezoelectric and surface scale effects are written as a matrix differential equation subject to tip-sample general boundary conditions and to compatibility conditions for the case of multispan beams. We considered a quadratic and a cubic eigenvalue problem related to the inclusion of smart materials and surface effects. Simulations were performed for a two stepped beam with a piezoelectric patch subject to pulse forcing terms. Results with Timoshenko models that include surface effects are presented for micro- and nanoscale. It was observed that the effects are significant just in nanoscale. We also simulate the frequency effects of a double-span beam in which one segment includes rotatory inertia and shear deformation and the other one neglects both phenomena. The proposed analytical methodology can be useful in the design of micro- and nanoresonator structures that involve deformable flexural models for detecting and imaging of physical and biochemical quantities.


---

## Body

## 1. Introduction

Atomic force microscopy (AFM) is a scanning probe microscopy (SPM) technique to obtain images of surface topography at the atomic scale, in a noninvasive manner, from a wide variety of samples on a scale from angstroms to 100 microns [1]. A typical AFM consists of a microcantilever with a sharp tip, a sample positioning system, a detection system, and a control system. The associated length scales are sufficiently small to call the applicability of classical continuum models into question [2]. In this work, we seek to develop a vibration dynamics framework for beams that include smart materials and subject to surface effects. This framework is also considered in the case of a two-span beam in which the first segment is governed by the Timoshenko model and the second segment is an Euler-Bernoulli beam model [3].Recently, new generations of active microcantilevers have included piezoelectric materials locally attached at the microbeam with the role of sensors and/or actuators linear and nonlinear [4–9], among others. This has led to the study of multispan beams for AFM. The inclusion of smart materials layers will modify material properties between neighboring segments, producing discontinuities and fulfilling compatibility conditions for the continuity of the displacement and rotation of the beam and for the equilibrium of bending moment and shear at the discontinuity points.Surface effects often play a significant role in the physical properties of micro- and nanosized materials and structures. Since the atoms within a very thin layer near surfaces experience a different local environment form that is experienced by atoms in bulk, the physical properties and mechanical response of surfaces will be distinct from those of bulk materials. Contrary to macroscopic structures, surface effects can strongly influence the stress and deformation properties of nanodevices. This latter is due to the increasing ratio between surface/interface area and volume. For instance, for a bulk deformation energyUb≈Ed3 (E elastic stiffness) and Us≈γd2 (γ surface energy per area), we shall have Us>Ub with distance d<γ/E. In [10], it has been shown that with a modified Timoshenko beam model the dynamics of nanoscaled tubes with surface effects deviate considerably from results obtained with classic theories.In some models, considering effective properties (D) of nanosized beams when compared with standard continuum mechanics (Dc), the relative error (D-Dc)/Dc scale is found as αS/Ed, where α is a geometrical constant depending upon the nanobeam considered, S is a surface elastic constant, E is a bulk elastic modulus, and d is a length defining the size of the structural element [11].  Here we shall consider free vibration of a nanoscale continuum beam models which incorporates surface energy [12–14]. The search for time exponential solutions leads to a quadratic and cubic eigenvalue problems associated with a second-order modal differential system and a singular third-order differential system, respectively.Simulations have been performed for triple span beams with and without a patch in the first segment and with loads in the first and last segment with tip-sample interaction. These loads included a uniform pulse load and its second derivative associated with the moment due to a piezoelectric patch. It is observed how the inclusion of piezoelectric materials absorbs vibrations when compared with classical multispan beams. Also, the robustness of the impulse response method with varying parameters and its influence in the behavior of the responses are observed.The size dependence effects in Timoshenko microbeams with surface effects (TMB) have been simulated for the nondimensional natural frequency and compared with those of the classical Timoshenko beam model (TB). It has been observed that for beam length on the order of nanometer to microns, the difference between natural frequencies is apparent, and by increasing the length of the microbeam, the results tend to Timoshenko classical theory; that is, the surface effects are significant only in nanoscale. This same behavior was observed in [14] for a microbeam simply supported.This work involves a methodology that may be adequately adapted to study recent advances in micro-/nanosized structures that are intensively used to design advanced micro-/nanosensors and molecular transportation systems devices for various engineering and medical applications. These sensors due to their ultrahigh resonant frequencies are important in sensitive sensing, molecular transportation, molecular separation, high-frequency signal processing, and biological imaging. Some deposited processes in the formation of thin-film composite membranes can be similar to types of multispan vibrating beams. However, the complex configuration of materials used in real world devices requires the study of material properties stemming from their atomically thin layered structures. Chemical deposited processes in the manufacturing of thin films on semiconductors can involve low temperatures with silicon rich nitride. Conventional structural analysis methods assume ideal structures (free of irregularities) but material and/or geometrical variations in a structure may result in drastic changes in its dynamical behaviors. Diffusion- and reaction-controlled interfacial polymerization is an important and practical topic that is beyond our scope [15–19].This paper is organized as follows. In Section2, we consider the Euler-Bernoulli multispan beam model, their compatibility conditions, and the boundary conditions in a matrix form. In Section 3, the dynamic response of the matrix Euler-Bernoulli multispan beam model subject to tip-sample interactions and external forcing is given in terms of the distributed matrix impulse response. The case of a microcantilever with a piezoelectric layer is discussed in Section 4. In Section 5 are presented the results of a double-span beam, obtained by the expansion of the first segment of the triple piezoelectric beam. In Sections 6 and 7 is considered the eigenanalysis of Timoshenko beam models with surface effects and their comparison with the Timoshenko classical model, for micro- and nanoscale. In Section 8, we consider the multispan Timoshenko beam model. AFM-based nanoscale processing with continuum surrounding media such as that found in biology and nanomachining applications [20–22] suggests observing frequency effects that arise with an academic two span beam model in which one segment includes rotatory inertia and shear deformation and the other one neglects both effects.


## 2. Flexural Vibrations Using the Euler-Bernoulli Multispan Beam Model

We consider a multispan microcantilever of lengthL composed of N segments in which the displacement vj(t,x) in the segment xj-1≤x≤xj is governed by the Euler-Bernoulli model(1)mj∂2vjt,x∂t2+kj∂4vjt,x∂x4=fjt,x,where mj=ρjAj is the mass per unit length, Aj is the cross section area, ρj is the mass density, kj=EjIj is the flexural stiffness, and fj(t,x) is the transverse dynamic load at the j-th station, j=1:N. Here x0=0 and XN=L.


### 2.1. Matrix Formulation

The stepped Euler-Bernoulli model (1) can be written as a second-order differential equation(2)Mv¨+Kv=F,with diagonal matrix coefficients (3)M=M10⋱0MN,K=EI∂4∂x4=EI10⋱0EIN∂4∂x4,where(4)Mj=ρjAj,EIj=EjIj.Here the displacement and the dynamic load v(t,x) and F(t,x), at the multispan beam, are given in vector form as (5)vt,x=v1t,xv2t,x⋮vNt,x,F=f1t,xf2t,x⋮fNt,x.


### 2.2. Boundary and Internal Conditions for a Multispan Cantilever Beam

For a cantilever beam, the boundary conditions can be written in a compact matrix form(6)B0v0=0,BLvL=0,where(7)B0=10000100,BL=00100001,(8)v0=v1t,0v1′t,0v1′′t,0v1′′′t,0,vL=vNt,LvN′t,LvN′′t,LvN′′′t,L.At the pointsxi, i=1,2,…,N-1, located between two consecutive segments, the compatibility conditions for the displacement, slope, bending moment, and shear lead to the compatibility matrix relationship(9)C1,ivit,xi=C2,ivi+1t,xi,with(10)C1,i=1000010000100001,C2,i=1000010000αi0000αi,vi=vit,xivi′t,xivi′′t,xivi′′′t,xi,being αi=EIi+1/EIi for i=1,…,N-1.


## 2.1. Matrix Formulation

The stepped Euler-Bernoulli model (1) can be written as a second-order differential equation(2)Mv¨+Kv=F,with diagonal matrix coefficients (3)M=M10⋱0MN,K=EI∂4∂x4=EI10⋱0EIN∂4∂x4,where(4)Mj=ρjAj,EIj=EjIj.Here the displacement and the dynamic load v(t,x) and F(t,x), at the multispan beam, are given in vector form as (5)vt,x=v1t,xv2t,x⋮vNt,x,F=f1t,xf2t,x⋮fNt,x.


## 2.2. Boundary and Internal Conditions for a Multispan Cantilever Beam

For a cantilever beam, the boundary conditions can be written in a compact matrix form(6)B0v0=0,BLvL=0,where(7)B0=10000100,BL=00100001,(8)v0=v1t,0v1′t,0v1′′t,0v1′′′t,0,vL=vNt,LvN′t,LvN′′t,LvN′′′t,L.At the pointsxi, i=1,2,…,N-1, located between two consecutive segments, the compatibility conditions for the displacement, slope, bending moment, and shear lead to the compatibility matrix relationship(9)C1,ivit,xi=C2,ivi+1t,xi,with(10)C1,i=1000010000100001,C2,i=1000010000αi0000αi,vi=vit,xivi′t,xivi′′t,xivi′′′t,xi,being αi=EIi+1/EIi for i=1,…,N-1.


## 3. Forced Vibrations

The dynamic response of a forced multispan Euler-Bernoulli modelMv¨+Kv=F given in (2), subject to homogeneous boundary conditions, could be described in terms of the N×N diagonal matrix impulse response or initial-value Green function h(t,x,ξ) satisfying the initial-value problem(11)Mh¨t,x,ξ+Kht,x,ξ=0,0<x<L,0<ξ<L,t>0,h0,x,ξ=0,Mht0,x,ξ=δx-ξI,B0h0=0,BLhL=0,C1,ihit,xi-,ξ=C2,ihit,xi+,ξ,where I is the N×N identity matrix and C1,i,C2,i, hi are as given in (10), by using the i-th diagonal component of h.To work on the frequency domain, we need to introduce transfer functionH(s,x,ξ) as being the Laplace transform of h(t,x,ξ) with respect to time (For inputs eλtF(x) we seek outputs eλtV(x) where V(x)=HF(x) with H defined in (17)). Thus(12)s2M+KHs,x,ξ=δx-ξI,0<x,ξ<L,B0H0=0,BLHL=0,C1,iHit,xi-,ξ=C2,iHit,xi+,ξ.It turns out thath(t-τ,x,ξ) acts as an integrating factor in Lagrange’s adjoint method for the forced equation (2). Multiplying (2) by h(t-τ,x,ξ) and integrating by parts, we obtain the dynamic response(13)vt,x=∫0t∫0Lhtτ,x,ξMvoξ+hτ,x,ξMv1ξdξdτ+∫0t∫0Lht-τ,x,ξFτ,ξdξdτ,wherevo(x)=v(0,x) and v1(x)=vt(0,x) are given initial conditions.The procedure mentioned above is also related to the Riemann function method for integrating partial differential equations. Dynamic responses have been considered in the field of control of distributed systems and in elastodynamics in connection with vibrations and cracking problems [23–26].In practice, when computing the convolution integral, which corresponds to the forced responsev(t,x), with null initial values v(0,x)=0,vt(0,x)=0, we actually have v(t,x)=vh(t,x)+vp(t,x), where vh(t,x) is a free vibration introduced by the system, whose initial values are a priori unknown, and vp(t,x) is a particular response. When vp(t,x) can be determined by other means, then those initial values can be supplied as vh(0,x)=-vp(0,x) and (vh)t(0,x)=-(vp)t(0,x) [27]. Thus, the induced system free response, due to a permanent response vp(t,x), is given by(14)vht,x=-∫0Lht,x,ξMv˙p0,ξdξ-∫0L∂ht,x,ξ∂tMvp0,ξdξ.


### 3.1. Frequency Response

Harmonic and piecewise linear forcing are of interest in frequency analysis. When seeking a response of the same type, the transfer function is introduced. Given the harmonic input(15)ft,x=eiωtvx,we have the harmonic output response(16)vpt,x=eiωtHiωvx,where(17)Hiωvx=∫0LHiω,x,ξvξdξ.In particular, for a concentrated force at a pointx=a, of spatial amplitude v(x)=v(x)δ(x-a), we have the permanent response(18)vpt,x=eiωtHiω,x,ava.With the initial values vp(0,ξ)=H(iω,ξ,a)v(a), v˙p(0,ξ)=iωvp(0,ξ), the induced free response is given by(19)vht,x=-∫0Lrt,x,ξ,ωMHiω,a,ξvadξ,where(20)r=htt,x,ξ+iωht,x,ξ.For a pulse amplitude(21)vx=voux-a-ux-b,where u(x) is the unit step function, the permanent response turns out(22)vpt,x=eiωt∫abHiω,x,ξvodξ.As before, by substituting the initial values in (14), the induced free response will now be(23)vht,x=-∫0Lrt,x,ξ,ωMHiω,x,ξvξdξ,with r given as in (20).


### 3.2. Closed Form of the Transfer Function for a Multispan Cantilever Beam

The Green function of a multispan Euler-Bernoulli can be determined by using an appropriate solution basis in each segment. From (12), H(s,x,ξ) is a diagonal matrix, with entries Hj that satisfy the differential equation(24)Xivx-εj4Xx=δx-ξ,where(25)εj4=-ρjAjs2EjIj,j=1,2,⋯,N.The functions hx-xj-1,h′x-xj-1,h′′x-xj-1,h′′′x-xj-1 constitute a solution basis of the equation X(iv)(x)-εj4X(x)=0. Here(26)hx=sinh⁡εjx-sin⁡εjx2εj3is the solution of the initial-value problem(27)hivx-εj4hx=0,h0=0,h′0=0,h′′0=0,h′′′0=1.By denoting(28)h0x=h′′′x=12cosh⁡εjx+cos⁡εjx,h1x=h′′x=12εsinh⁡εjx+sin⁡εjx,h2x=h′x=12εj2cosh⁡εjx-cos⁡εjx,h3x=hx=12εj3sinh⁡εjx-sin⁡εjx,in each segment we have(29)Hs,x,ξ=H00⋯00H1⋯0⋮⋮⋮⋮00⋯HN-1,for j=0,..,N-1(30)Hjs,x,ξ=h0x-xjHs,xj,ξ+h1x-xjHxs,xj,ξ+h2x-xjHxxs,xj,ξ+h3x-xjHxxxs,xj,ξ+h3x-xjux-ξ,xj≤x≤xj+1.In order to determine H(s,x,ξ), we need to find the initial-value vectors(31)H=Hs,xj,ξHxs,xj,ξHxxs,xj,ξHxxxs,xj,ξ.By using the compatibility conditions(32)C1,iHs,xj-,ξ=C2,iHs,xj+,ξ,j=1,2,⋯,N-1,given similarly to (9), and by deriving (30) with respect to x.The above process will be carried out for the case of a two stepped cantilever beam with the intermediate discontinuity pointx1=l1.The boundary conditions and the initial values ofh(x) lead to a simplification of H in the first and last segments. We have two null initial values at x=0, and we can use the boundary conditions at x=L as initial values, instead of those at xN-1. This is possible by changing the solution basis at the last segment in (30) to the one generated by h(L-x).Thus, for0≤ξ≤L(33)Hs,x,ξ=h2xHxxs,0,ξ+h3xHxxxs,0,ξ+h3x-ξux-ξ,0≤x≤x1,Hs,x,ξ=h0L-xHs,L,ξ+h1L-xHxs,L,ξ+h3ξ-xuξ-x,x1≤x≤L.By applying the compatibility conditions(34)Hs,l1-,ξ=Hs,l1+,ξ,Hxs,l1-,ξ=Hxs,l1+,ξ,E1I1Hxxs,l1-,ξ=E2I2Hxxs,l1+,ξ,E1I1Hxxxs,l1-,ξ=E2I2Hxxxs,l1+,ξ,we get the algebraic system(35)h2l1h3l1-h0L-xx=l1-h1L-xx=l1h2′l1h3′l1-h0′L-xx=l1-h1′L-xx=l1h2′′l1h3′′l1-E2I2E1I1h0′′L-xx=l1-E2I2E1I1h1′′L-xx=l1h2′′′l1h3′′′l1-E2I2E1I1h0′′′L-xx=l1-E2I2E1I1h1′′′L-xx=l1H′′s,0,ξH′′′s,0,ξHs,L,ξH′s,L,ξ=hl1-ξul1-ξ-hl1-ξul1-ξh′l1-ξul1-ξ-h′l1-ξul1-ξE2I2E1I1h′′ξ-xx=l1uξ-l1-h′′l1-ξul1-ξE2I2E1I1h′′′ξ-xx=l1uξ-l1-h′′′l1-ξul1-ξ,for determining Hj(s,x,ξ) from (30).Remark. The case of axially loaded stepped beams has been investigated in [28] with nonclassical conditions. There, the Green function was obtained by working with a solution basis that considers hyperbolic and trigonometric functions in the general solution. These functions depend upon simple roots of a complete quadratic polynomial due to the inclusion of a second derivative spatial term in (24). When the axial load is removed, the solution basis reduces to the one considered for particular cases of vibrating beams. We observe that through a limit process the use of h(x), which involves a division by root parameters, can also handle the case of repeated roots, which is of interest in the static case or critical frequencies.


## 3.1. Frequency Response

Harmonic and piecewise linear forcing are of interest in frequency analysis. When seeking a response of the same type, the transfer function is introduced. Given the harmonic input(15)ft,x=eiωtvx,we have the harmonic output response(16)vpt,x=eiωtHiωvx,where(17)Hiωvx=∫0LHiω,x,ξvξdξ.In particular, for a concentrated force at a pointx=a, of spatial amplitude v(x)=v(x)δ(x-a), we have the permanent response(18)vpt,x=eiωtHiω,x,ava.With the initial values vp(0,ξ)=H(iω,ξ,a)v(a), v˙p(0,ξ)=iωvp(0,ξ), the induced free response is given by(19)vht,x=-∫0Lrt,x,ξ,ωMHiω,a,ξvadξ,where(20)r=htt,x,ξ+iωht,x,ξ.For a pulse amplitude(21)vx=voux-a-ux-b,where u(x) is the unit step function, the permanent response turns out(22)vpt,x=eiωt∫abHiω,x,ξvodξ.As before, by substituting the initial values in (14), the induced free response will now be(23)vht,x=-∫0Lrt,x,ξ,ωMHiω,x,ξvξdξ,with r given as in (20).


## 3.2. Closed Form of the Transfer Function for a Multispan Cantilever Beam

The Green function of a multispan Euler-Bernoulli can be determined by using an appropriate solution basis in each segment. From (12), H(s,x,ξ) is a diagonal matrix, with entries Hj that satisfy the differential equation(24)Xivx-εj4Xx=δx-ξ,where(25)εj4=-ρjAjs2EjIj,j=1,2,⋯,N.The functions hx-xj-1,h′x-xj-1,h′′x-xj-1,h′′′x-xj-1 constitute a solution basis of the equation X(iv)(x)-εj4X(x)=0. Here(26)hx=sinh⁡εjx-sin⁡εjx2εj3is the solution of the initial-value problem(27)hivx-εj4hx=0,h0=0,h′0=0,h′′0=0,h′′′0=1.By denoting(28)h0x=h′′′x=12cosh⁡εjx+cos⁡εjx,h1x=h′′x=12εsinh⁡εjx+sin⁡εjx,h2x=h′x=12εj2cosh⁡εjx-cos⁡εjx,h3x=hx=12εj3sinh⁡εjx-sin⁡εjx,in each segment we have(29)Hs,x,ξ=H00⋯00H1⋯0⋮⋮⋮⋮00⋯HN-1,for j=0,..,N-1(30)Hjs,x,ξ=h0x-xjHs,xj,ξ+h1x-xjHxs,xj,ξ+h2x-xjHxxs,xj,ξ+h3x-xjHxxxs,xj,ξ+h3x-xjux-ξ,xj≤x≤xj+1.In order to determine H(s,x,ξ), we need to find the initial-value vectors(31)H=Hs,xj,ξHxs,xj,ξHxxs,xj,ξHxxxs,xj,ξ.By using the compatibility conditions(32)C1,iHs,xj-,ξ=C2,iHs,xj+,ξ,j=1,2,⋯,N-1,given similarly to (9), and by deriving (30) with respect to x.The above process will be carried out for the case of a two stepped cantilever beam with the intermediate discontinuity pointx1=l1.The boundary conditions and the initial values ofh(x) lead to a simplification of H in the first and last segments. We have two null initial values at x=0, and we can use the boundary conditions at x=L as initial values, instead of those at xN-1. This is possible by changing the solution basis at the last segment in (30) to the one generated by h(L-x).Thus, for0≤ξ≤L(33)Hs,x,ξ=h2xHxxs,0,ξ+h3xHxxxs,0,ξ+h3x-ξux-ξ,0≤x≤x1,Hs,x,ξ=h0L-xHs,L,ξ+h1L-xHxs,L,ξ+h3ξ-xuξ-x,x1≤x≤L.By applying the compatibility conditions(34)Hs,l1-,ξ=Hs,l1+,ξ,Hxs,l1-,ξ=Hxs,l1+,ξ,E1I1Hxxs,l1-,ξ=E2I2Hxxs,l1+,ξ,E1I1Hxxxs,l1-,ξ=E2I2Hxxxs,l1+,ξ,we get the algebraic system(35)h2l1h3l1-h0L-xx=l1-h1L-xx=l1h2′l1h3′l1-h0′L-xx=l1-h1′L-xx=l1h2′′l1h3′′l1-E2I2E1I1h0′′L-xx=l1-E2I2E1I1h1′′L-xx=l1h2′′′l1h3′′′l1-E2I2E1I1h0′′′L-xx=l1-E2I2E1I1h1′′′L-xx=l1H′′s,0,ξH′′′s,0,ξHs,L,ξH′s,L,ξ=hl1-ξul1-ξ-hl1-ξul1-ξh′l1-ξul1-ξ-h′l1-ξul1-ξE2I2E1I1h′′ξ-xx=l1uξ-l1-h′′l1-ξul1-ξE2I2E1I1h′′′ξ-xx=l1uξ-l1-h′′′l1-ξul1-ξ,for determining Hj(s,x,ξ) from (30).Remark. The case of axially loaded stepped beams has been investigated in [28] with nonclassical conditions. There, the Green function was obtained by working with a solution basis that considers hyperbolic and trigonometric functions in the general solution. These functions depend upon simple roots of a complete quadratic polynomial due to the inclusion of a second derivative spatial term in (24). When the axial load is removed, the solution basis reduces to the one considered for particular cases of vibrating beams. We observe that through a limit process the use of h(x), which involves a division by root parameters, can also handle the case of repeated roots, which is of interest in the static case or critical frequencies.


## 4. A Cantilever Beam with a Piezoelectric Layer in AFM

Active microcantilever beams due to their structural flexibility and sensitivity to atomic and molecular forces have received increased attention in a variety of nanoscale sensing and measuring applications, including atomic force microscopy, thermal scanning microscopy, and biomass sensing.In [9], a microcantilever model was proposed for studying atomic force microscope. It was formulated as a three stepped beam with a piezoelectric layer patch in the first segment and the other two segments were simple beams with different cross-sectional areas according to Figure 1(a).Figure 1
Schematic representation of the (a) triple span and (b) double-span piezoelectric beams.
(a) (b)The governing equations can be written as a system of three second-order partial differential equations with constant coefficients(36)m1∂2u1t,x∂t2+EI1∂4u1t,x∂x4=∂2Mpt,x∂x2,0≤x≤l1,m2∂2u2t,x∂t2+EI2∂4u2t,x∂x4=0,l1<x≤l2,m3∂2u3t,x∂t2+EI3∂4u3t,x∂x4=0l2<x≤L,subject to cantilever boundary conditions (7) and compatibility conditions (9) at the discontinuity points x1,x2. The parameters mj and (EI)j for each segment of the beam are given by(37)m1=ρA1=ρpwptp+ρbwb1tb,m2=ρA2=ρbwb1tb,m3=ρA3=ρbwb2tb,(38)EI1=EpIpp+EbIpb+EbIb1,EI2=EbIb1,EI3=EbIb2.The parametersρp, ρb are the mass density of piezoelectric material and mass density of beam material; wp, wb1, and wb2 are the piezoelectric layer width, the width of the second, and third beam segments, respectively; tp and tb are the piezoelectric layer thickness and beam thickness. Ep and Eb are the Young modulus of piezoelectric layer material and Young modulus of second and third segment material. Ib1, Ib2, Ipp, Ipb are the beam moment of inertia, layer moment of inertia, and layer-beam moment of inertia and are calculated as Ib1=wb1tb3/12, Ib2=wb2tb3/12, Ipp=zn2tpwp+((1/3)tp3+(1/2)tbtp2+(1/4)tb2tp)wp-zn(tp2+tbtp)wp; Ipb=zn2tbwb1. And zn is the neutral axis of the beam on the composite portion expressed by(39)zn=12Eptpwptp+tbEptpwp+Ebtbwb1.The values of these parameters [9] are presented in Table 1.Table 1
Parameters of piezoelectric multispan beam.
Parameter Numeric Value Unity LengthL 486×10-6 m Lengthl1 325×10-6 m Lengthl2 360×10-6 m Widthwb1 250×10-6 m Widthwb2 55×10-6 m Widthwp 130×10-6 m Thicknesstp 4×10-6 m Thicknesstb 4×10-6 m Densityρb 2330 Kg/m3 Densityρp 6390 Kg/m3 Young’s modulusEb 105 GPa Young’s modulusEp 104 GPaThe distributed cross-sectional momentMp(x,t)=Mpo(t)(1-u(x-l1)) in the first segment induces a forcing term. Its time amplitude is given by(40)Mpot=12wpEpd31tb+tpvt,where d31 and v(t) are the coefficients of the converse piezoelectric effect and the applied voltage, respectively.


### 4.1. Forced Responses

For homogeneous boundary conditions and null initial data, we have from (13) the forced response(41)vt,x=∫0t∫0Lht-τ,x,ξFτ,ξdξdτ,Ft,x=∂2Mpt,x∂x200,where ∂2Mp(τ,ξ)/∂ξ2=-Mpoτδ′(ξ-l1). Thus(42)vt,x=∫0t∫0Lht-τ,x,ξ∂2Mpτ,ξ∂ξ2dξdτ=∫0tMpoτ∂h∂ξt-τ,x,l1dτ.By assuming a discrete complete spectrum of natural frequencies ωj and the normal mode property, the Galerkin method can be used for determining approximate forced vibrations. We have the modal series(43)ht,x,ξ=∑n=1∞Xnx,ξsin⁡ωntωn,0≤x,ξ≤L,and the spatial Green function(44)Ht,x,ξ=∑n=1∞Xnx,ξ1s2+ωn2,where(45)Xnx,ξ=X1nxX1nξ0⋯00X2nxX2nξ⋯0⋮⋮⋮⋮00⋯XN-1nxXN-1nξ.Here X(x)(46)Xx=X1xX2x⋮XNxis the spatial amplitude of exponential solutions v(t,x)=eλtX(x) of the unforced Euler-Bernoulli model Mv¨+Kv=0 that satisfy the boundary conditions (6), and the compatibility conditions (9) at the discontinuity points x1 and x2.The modal approximation(47)∂h∂ξt-τ,x,ξ=∑n=1m∂Xn∂ξx,ξsin⁡ωnt-τωnhas been used in simulations with and without the inclusion of a piezoelectric layer. Expansions were truncated with a small numberm of terms, usually between 5 and 10 [30].In Table2 are showed the first six natural frequencies obtained in this work for a triple span beam with piezoelectric patch and without it. The first three frequencies are compared with those obtained experimentally in [9]. There is an agreement with the methodology of this work and it can be observed that the elimination of piezoelectric layers diminishes the frequency values. In Figure 2 are shown the corresponding vibration modes; in the first segment the modes of piezoelectric beam have smaller amplitude.Table 2
Comparative natural frequencies between cases of double and triple span beams with and without piezoelectric patch on the first segment.
Case Piezoelectric path Natural frequenciesωi (Hz) ω 1 ω 2 ω 3 ω 4 ω 5 ω 6 Triple span [9]∗ with piezoelectric path 52.3 203.0 382.5 - - - Triple span (Figure1(a)) with piezoelectric path 50.6 209.2 393.7 910.7 1412.5 2016.7 Double span (Figure1(b)) with piezoelectric path 46.2 222.9 376.2 875.3 1537.5 1948.6 Double/Triple span without piezoelectric path 27.3 142.9 320.1 618.9 1073.0 1579.4  ∗Experimentally obtainedFigure 2
First six beam vibration modes: beam with piezoelectric layer (black line); beam (gray line) 1st segment (dash)[0,l1], l1=325×10-6m; 2nd segment (dot) [l1,l2], l2=360×10-6m; 3rd segment (solid) [l2,L], L=486×10-6m.The effects of the inclusion of a piezoelectric patch are appreciated by applying a harmonic rectangular pulse load in the first segment and in the third segment for a three stepped beam with and without piezoelectric layer. It is observed that the spatial amplitude diminishes with the inclusion of a piezoelectric patch when the pulse is positioned in the piezoelectric layer segment; this effect does not appear when the pulse is positioned in the third segment (Figure3).Figure 3
Left: forced responsev(t,x) with fixed t for f1(x)=sin(ωt)(u(x-a)-u(x-b)), a=100×10-6m, and b=150×10-6m. Right: forced response v(t,x) with fixed t for f3(t,x)=sin⁡(ωt)(u(x-c)-u(x-d)), c=430×10-6m, and d=440×10-6m. Solid black line: beam with piezoelectric layer. Dashed gray line: beam without piezoelectric layer.In Figure4 are compared the forced responses due to forcing terms only in the first segment. The first output tracks the forcing input f1(t,x)=Mp0(t)(1-u(x-l1)). For the original forcing term proposed in [9] f1(t,x)=∂2Mpt,x/∂x2=-Mp0tδ′(x-l1), we can observe that the effect of doublet due to derivation of Heaviside function is to concentrate the output at the point x=l1. Also, the amplitudes of forced response are smaller for three stepped beam with piezoelectric layer.Figure 4
Left: forced responsev(t,x) for f1(t,x)=v(t)(1-u(x-l1)). Right: forced response v(t,x) for f1(t,x)=v(t)(∂2Mp(x)/∂x2). Beam without piezoelectric layer (dashed line) and beam with piezoelectric layer (solid line), v(t)=sin(ωt).


## 4.1. Forced Responses

For homogeneous boundary conditions and null initial data, we have from (13) the forced response(41)vt,x=∫0t∫0Lht-τ,x,ξFτ,ξdξdτ,Ft,x=∂2Mpt,x∂x200,where ∂2Mp(τ,ξ)/∂ξ2=-Mpoτδ′(ξ-l1). Thus(42)vt,x=∫0t∫0Lht-τ,x,ξ∂2Mpτ,ξ∂ξ2dξdτ=∫0tMpoτ∂h∂ξt-τ,x,l1dτ.By assuming a discrete complete spectrum of natural frequencies ωj and the normal mode property, the Galerkin method can be used for determining approximate forced vibrations. We have the modal series(43)ht,x,ξ=∑n=1∞Xnx,ξsin⁡ωntωn,0≤x,ξ≤L,and the spatial Green function(44)Ht,x,ξ=∑n=1∞Xnx,ξ1s2+ωn2,where(45)Xnx,ξ=X1nxX1nξ0⋯00X2nxX2nξ⋯0⋮⋮⋮⋮00⋯XN-1nxXN-1nξ.Here X(x)(46)Xx=X1xX2x⋮XNxis the spatial amplitude of exponential solutions v(t,x)=eλtX(x) of the unforced Euler-Bernoulli model Mv¨+Kv=0 that satisfy the boundary conditions (6), and the compatibility conditions (9) at the discontinuity points x1 and x2.The modal approximation(47)∂h∂ξt-τ,x,ξ=∑n=1m∂Xn∂ξx,ξsin⁡ωnt-τωnhas been used in simulations with and without the inclusion of a piezoelectric layer. Expansions were truncated with a small numberm of terms, usually between 5 and 10 [30].In Table2 are showed the first six natural frequencies obtained in this work for a triple span beam with piezoelectric patch and without it. The first three frequencies are compared with those obtained experimentally in [9]. There is an agreement with the methodology of this work and it can be observed that the elimination of piezoelectric layers diminishes the frequency values. In Figure 2 are shown the corresponding vibration modes; in the first segment the modes of piezoelectric beam have smaller amplitude.Table 2
Comparative natural frequencies between cases of double and triple span beams with and without piezoelectric patch on the first segment.
Case Piezoelectric path Natural frequenciesωi (Hz) ω 1 ω 2 ω 3 ω 4 ω 5 ω 6 Triple span [9]∗ with piezoelectric path 52.3 203.0 382.5 - - - Triple span (Figure1(a)) with piezoelectric path 50.6 209.2 393.7 910.7 1412.5 2016.7 Double span (Figure1(b)) with piezoelectric path 46.2 222.9 376.2 875.3 1537.5 1948.6 Double/Triple span without piezoelectric path 27.3 142.9 320.1 618.9 1073.0 1579.4  ∗Experimentally obtainedFigure 2
First six beam vibration modes: beam with piezoelectric layer (black line); beam (gray line) 1st segment (dash)[0,l1], l1=325×10-6m; 2nd segment (dot) [l1,l2], l2=360×10-6m; 3rd segment (solid) [l2,L], L=486×10-6m.The effects of the inclusion of a piezoelectric patch are appreciated by applying a harmonic rectangular pulse load in the first segment and in the third segment for a three stepped beam with and without piezoelectric layer. It is observed that the spatial amplitude diminishes with the inclusion of a piezoelectric patch when the pulse is positioned in the piezoelectric layer segment; this effect does not appear when the pulse is positioned in the third segment (Figure3).Figure 3
Left: forced responsev(t,x) with fixed t for f1(x)=sin(ωt)(u(x-a)-u(x-b)), a=100×10-6m, and b=150×10-6m. Right: forced response v(t,x) with fixed t for f3(t,x)=sin⁡(ωt)(u(x-c)-u(x-d)), c=430×10-6m, and d=440×10-6m. Solid black line: beam with piezoelectric layer. Dashed gray line: beam without piezoelectric layer.In Figure4 are compared the forced responses due to forcing terms only in the first segment. The first output tracks the forcing input f1(t,x)=Mp0(t)(1-u(x-l1)). For the original forcing term proposed in [9] f1(t,x)=∂2Mpt,x/∂x2=-Mp0tδ′(x-l1), we can observe that the effect of doublet due to derivation of Heaviside function is to concentrate the output at the point x=l1. Also, the amplitudes of forced response are smaller for three stepped beam with piezoelectric layer.Figure 4
Left: forced responsev(t,x) for f1(t,x)=v(t)(1-u(x-l1)). Right: forced response v(t,x) for f1(t,x)=v(t)(∂2Mp(x)/∂x2). Beam without piezoelectric layer (dashed line) and beam with piezoelectric layer (solid line), v(t)=sin(ωt).


## 5. Double-Span Microbeam with Material and Geometric Discontinuities

We consider a double-span beam of lengthL=486×10-6m, formed for two segments with different rectangular transversal section area, according to Figure 1(b). For simulation purposes, we consider the same parameters utilized in [9] for the three stepped beam with a piezoelectric patch in the first segment, and now for the double-span beam, obtained by removing the second segment of the first three stepped beam, which means extending the piezoelectric layer from l1 to l2. In our methodology, this amounts to consider the limit case when l1 tends to l2, so that the flexural stiffness (EI)1 and mass density (ρA)1 will be valid on the first segment [0,l2], since l1 was extended to l2, and on the new second segment [l2,L] we will have (EI)3 and (ρA)3, described in (37) and (38).The first six natural frequencies of this beam, with and without piezoelectric patch, were given in the third and fourth rows of Table2 and were obtained by the modal method. We observed that the reduction of three segments for two segments, maintaining the others material and geometric properties, influences the magnitude of the natural frequencies. When the piezoelectric patch is removed, then the frequencies also decrease. The results were the same as those of the case of three stepped beam when removing the piezoelectric patch.In Figure5, by using the methodology described from (24) to (35), is presented a 3D graphic of the Green function H(s,x,ξ) of the two stepped piezoelectric beam, for the sixth natural frequency. The amplitude of this Green function has peaks at the same natural frequencies obtained by the modal method (Table 2, third row).Figure 5
Left: frequency functionH(iω6,x,ξ). Right: modal frequency response H(iω,xo,ξo) with the first three natural frequencies: 46.2KHz, 222.9KHz, and 376.2Khz at a given point.


## 6. Surface Elasticity and Residual Surface Tension in the Timoshenko Beam Model

The Timoshenko beam model has been modified with the inclusion of the residual tension and surface elastic modulus [14]. By using the energy method with the elastic strain energy of surface induced in the potential energy and residual surface tensions, expressed by Laplace-Young equations, it was obtained the modified Timoskenko model(48)ρAwtt-κGAwxx+κGA¯ψx=ft,x,ρIψtt-EI¯ψxx-κGAwx-ψ=qt,x,for the flexural deflection w(t,x) of the beam and the rotation angle ψ(t,x) of the cross section of the beam. Here f(t,x) denotes a transverse dynamic load and q(t,x) a moment load. The values(49)κGA¯=κGA-τu+τbb,EI¯=EI+2bh2Esare the effective curvature effect and flexural rigidity, respectively. Here τu and τb denote the upper and lower surfaces residual tensions and Es is a surface elastic modulus. Without these later parameters, the model (48) that includes surface effects becomes the classical Timoshenko model. As usual, Young’s modulus, moment of inertia of the cross section area, shear modulus, mass density area, and shear deformation factor of the beam are denoted as E, I, G, ρ, and κ, respectively. It is assumed that the beam has length L and area A=2bh, where b, 2h are the width and thickness of the beam and I=2bh3/3. The boundary conditions are those of a cantilever beam or subject to balance of the moment and shear at the free end x=L [31].


### 6.1. Matrix Formulation

The coupled Timoshenko model (48) can be written as a second-order differential equation with matrix coefficients(50)Mv¨+Kv=F,where(51)v=wt,xψt,x,F=ft,xqt,x,M=ρA00ρI,K=K2∂2∂x2+K1∂∂x+K0,with(52)K2=-κGA00-EI¯,K1=0κGA¯-κGA0,K0=000κGA.For a microcantilever beam of length L, we have the boundary conditions in a matrix formulation(53)1001vt,0+0000vxt,0=0,000-1vt,L+0110vxt,L=0,or in a more compact form(54)Avt,0+Bvxt,0=0,Pvt,L+Qvxt,L=0.


### 6.2. Eigenanalysis

The search of exponential solutions(55)vt,x=eλtvx,vx=wxψx,of the unforced Timoshenko model(56)M∂2v∂t2+Kv=0,subject to general separated homogeneous boundary conditions (54), results in that v is the general solution of a second-order matrix differential equation and is given by(57)vx=hxc1+h′xc2,for constant 2×1 vectors c1 and c2. Here h(x) is the 2×2 matrix solution of the initial-value problem(58)Mh′′x+Ch′x+Kλhx=0,h0=0,Mh′0=I,where 0 denotes the 2×2 null matrix and I the 2×2 identity matrix. The matrix coefficients are(59)M=-κGA00-EI¯,C=0κGA¯-κGA0,Kλ=ρAλ200λ2ρI+κGA.By using the initial values ofh(x) in (57) and the clamped boundary condition v(0)=0, b, it turns out that c2=0. Thus we have to determine λ so that(60)vx=hx,λc1satisfies the boundary condition at the free end x=L. By assuming homogeneous boundary conditions, we have the nonlinear eigenmatrix problem(61)Uλc=PhL,λ+Qh′L,λc1=0.From this, it turns out the characteristic equation(62)Δλ=det⁡U=0.We thus have the reduced system(63)UDc1=0,c1=c11c12T,with(64)UD=ad′′L-ad′′′L+λ2cd′L-bmd′′′L+λ2ed′Lad′′L-λ2cdL-amd′′L.where(65)dx=δsenh⁡εx-εsen⁡δxabmεδδ2+ε2is the solution of the initial-value problem(66)abmdivx+-aeλ2-cλ2bm-a2+amad′′x+cλ2a+cλ4edx=0,d0=d′0=d′′0=0,abmd′′′0=1,with(67)ε=12-2g2+2g4+4r4¯,δ=122g2+2g4+4r4¯,(68)g2=gm2¯+s¯2,gm2¯=-ebm+caλ2,s¯2=1bmam-a,r4¯=-cλ2a+eλ2abm.Here(69)a=κGA,c=ρA,e=ρI,am=κGA¯=κGA-τu+τbb,bm=EI¯=EI+2bh2Es.From (62), the natural frequencies λ=iω will satisfy(70)Δiω=-ω4ced′L2+ω2acdLd′′L-ae+cbmd′Ld′′′L+a2-aamd′′L2-abmd′′′L2=0.The fundamental matrix responseh(x) is obtained in a similar way to the classical case [32] (71)hx=a+eλ2dx-bmd′′x-amd′xad′x-ad′′x+cλ2dx.We observe from Figure6 that for frequencies not above the critical frequency, the fundamental response d(x) has lower amplitude when considering surface effects, while for higher frequencies the oscillatory behavior of d(x) follows the classical one.Figure 6
Behavior of the fundamental solutiond(x), according to critical natural frequency ωc=a/e, using L=5×10-9 and Al[111] beam configuration, including surface effects (dashed line) and classical case (solid line).For a microcantilever beam described by the Timoshenko model with surface effects, the size dependence in the natural frequency of Timoshenko classical model and Timoshenko model including surface effects is illustrated in Figure7. The solutions based on classical Timoshenko beam theory and Timoshenko beam theory including surface effects are denoted by TB and TMB, respectively. The fundamental natural frequencies are normalized to fundamental frequency of cantilever Euler-Bernoulli beam. In this figure are considered the parameters utilized in [14] for the same purposes. The parameters of surface elasticity and residual surface tension can be determined by molecular dynamics simulations or experiments. Residual surface stresses can be either positive or negative, depending on the crystallographic structure [33]. For an anodic alumina Al (Young’s modulus E=70GPa, Poisson’s ratio ν=0.3 and ρ=2700kg/m3) are considered two types of crystallographic direction(72)Al100:Es=-7.9253N/mandτ=0.5689N/m,Al111:Es=5.1882N/mandτ=0.9108N/m.Figure 7
Influence of surface effects and size dependence on the normalized fundamental natural frequency of the microcantilever for 2h=0.2L, b=0.4L, andκ=5/6.In Figure7, the size dependence effects in the nondimensional natural frequency of TMB microbeams in comparison to the classical TB are illustrated. We can observe that for beam length on the order of nanometer to microns, the difference between natural frequencies is apparent and, by increasing the length of the microbeam, the results tend to Timoshenko classical theory; that is, the surface effects are significant only in nanoscale. This same behavior was observed in [14] for a microbeam simply supported.


## 6.1. Matrix Formulation

The coupled Timoshenko model (48) can be written as a second-order differential equation with matrix coefficients(50)Mv¨+Kv=F,where(51)v=wt,xψt,x,F=ft,xqt,x,M=ρA00ρI,K=K2∂2∂x2+K1∂∂x+K0,with(52)K2=-κGA00-EI¯,K1=0κGA¯-κGA0,K0=000κGA.For a microcantilever beam of length L, we have the boundary conditions in a matrix formulation(53)1001vt,0+0000vxt,0=0,000-1vt,L+0110vxt,L=0,or in a more compact form(54)Avt,0+Bvxt,0=0,Pvt,L+Qvxt,L=0.


## 6.2. Eigenanalysis

The search of exponential solutions(55)vt,x=eλtvx,vx=wxψx,of the unforced Timoshenko model(56)M∂2v∂t2+Kv=0,subject to general separated homogeneous boundary conditions (54), results in that v is the general solution of a second-order matrix differential equation and is given by(57)vx=hxc1+h′xc2,for constant 2×1 vectors c1 and c2. Here h(x) is the 2×2 matrix solution of the initial-value problem(58)Mh′′x+Ch′x+Kλhx=0,h0=0,Mh′0=I,where 0 denotes the 2×2 null matrix and I the 2×2 identity matrix. The matrix coefficients are(59)M=-κGA00-EI¯,C=0κGA¯-κGA0,Kλ=ρAλ200λ2ρI+κGA.By using the initial values ofh(x) in (57) and the clamped boundary condition v(0)=0, b, it turns out that c2=0. Thus we have to determine λ so that(60)vx=hx,λc1satisfies the boundary condition at the free end x=L. By assuming homogeneous boundary conditions, we have the nonlinear eigenmatrix problem(61)Uλc=PhL,λ+Qh′L,λc1=0.From this, it turns out the characteristic equation(62)Δλ=det⁡U=0.We thus have the reduced system(63)UDc1=0,c1=c11c12T,with(64)UD=ad′′L-ad′′′L+λ2cd′L-bmd′′′L+λ2ed′Lad′′L-λ2cdL-amd′′L.where(65)dx=δsenh⁡εx-εsen⁡δxabmεδδ2+ε2is the solution of the initial-value problem(66)abmdivx+-aeλ2-cλ2bm-a2+amad′′x+cλ2a+cλ4edx=0,d0=d′0=d′′0=0,abmd′′′0=1,with(67)ε=12-2g2+2g4+4r4¯,δ=122g2+2g4+4r4¯,(68)g2=gm2¯+s¯2,gm2¯=-ebm+caλ2,s¯2=1bmam-a,r4¯=-cλ2a+eλ2abm.Here(69)a=κGA,c=ρA,e=ρI,am=κGA¯=κGA-τu+τbb,bm=EI¯=EI+2bh2Es.From (62), the natural frequencies λ=iω will satisfy(70)Δiω=-ω4ced′L2+ω2acdLd′′L-ae+cbmd′Ld′′′L+a2-aamd′′L2-abmd′′′L2=0.The fundamental matrix responseh(x) is obtained in a similar way to the classical case [32] (71)hx=a+eλ2dx-bmd′′x-amd′xad′x-ad′′x+cλ2dx.We observe from Figure6 that for frequencies not above the critical frequency, the fundamental response d(x) has lower amplitude when considering surface effects, while for higher frequencies the oscillatory behavior of d(x) follows the classical one.Figure 6
Behavior of the fundamental solutiond(x), according to critical natural frequency ωc=a/e, using L=5×10-9 and Al[111] beam configuration, including surface effects (dashed line) and classical case (solid line).For a microcantilever beam described by the Timoshenko model with surface effects, the size dependence in the natural frequency of Timoshenko classical model and Timoshenko model including surface effects is illustrated in Figure7. The solutions based on classical Timoshenko beam theory and Timoshenko beam theory including surface effects are denoted by TB and TMB, respectively. The fundamental natural frequencies are normalized to fundamental frequency of cantilever Euler-Bernoulli beam. In this figure are considered the parameters utilized in [14] for the same purposes. The parameters of surface elasticity and residual surface tension can be determined by molecular dynamics simulations or experiments. Residual surface stresses can be either positive or negative, depending on the crystallographic structure [33]. For an anodic alumina Al (Young’s modulus E=70GPa, Poisson’s ratio ν=0.3 and ρ=2700kg/m3) are considered two types of crystallographic direction(72)Al100:Es=-7.9253N/mandτ=0.5689N/m,Al111:Es=5.1882N/mandτ=0.9108N/m.Figure 7
Influence of surface effects and size dependence on the normalized fundamental natural frequency of the microcantilever for 2h=0.2L, b=0.4L, andκ=5/6.In Figure7, the size dependence effects in the nondimensional natural frequency of TMB microbeams in comparison to the classical TB are illustrated. We can observe that for beam length on the order of nanometer to microns, the difference between natural frequencies is apparent and, by increasing the length of the microbeam, the results tend to Timoshenko classical theory; that is, the surface effects are significant only in nanoscale. This same behavior was observed in [14] for a microbeam simply supported.


## 7. Timoshenko Model with Surface Energy in Thick Nanobeams

In [13], the surface effects are included in the Timoshenko beam model following Gurtin-Murdoch continuum theory [34], in which is considered an elastic surface with zero thickness fully bonded to its bulk material. This elastic surface adds a set of specific constitutive equations relative to distinct material properties of the surface and surface energy effects. Let us denote an isotropic material beam, with rectangular transversal section, length L, width b, thickness H=2h, moment of inertia I=2bh3/3, and A=2bh which is the transversal sectional area. The equations for the bulk are the same constitutive equations of classical theory of elasticity. The government equations of transverse and angular displacement, w(t,x) and ϕ(t,x), respectively, are then given by(73)ρA+ρ0s∗wtt-κGA+τ0s∗wxx+κGAϕx=0,ρI+ρ0I∗ϕtt-EI+2μ0+λ0I∗ϕxx-2νIρ0Hwxtt+2νIτ0Hwxxx-κGAwx-ϕ=0.In a compact notation we have(74)Mv¨+Kv=0,where(75)v=wϕ,M=M0+M1∂∂x,M0=cm00em,M1=00-I00,andK is a sum of matrix differential operators(76)K=K3∂3∂x3+K2∂2∂x2+K1∂∂x+K0,(77)K3=00I10,K2=-am00-bm,K1=0a-a0,K0=000a,with cm=ρA+ρ0s∗,am=κGA+τ0s∗,a=κGA,em=ρI+ρ0I∗,I0=2νIρ0/H,I1=2νIτ0/H.The material parameters of bulk areE, ν, and ρ, which represent Young’s modulus, Poisson’s ratio, and mass density, respectively. Also, we have beam surface parameters λ0, μ0, τ0, ρ0, Lamé constants, surface residual stress, and mass density of surface. I∗ is the moment of inertia relative to the surface I∗=2bh2+4h3/3, and for rectangular transversal section s∗=2b. κ is the shear correction coefficient.Rewriting (74), we have(78)cm∂2∂t2-am∂2∂x2a∂∂x-I0∂3∂t2∂x+I1∂3∂x3-a∂∂xem∂2∂t2-bm∂2∂x2+av=0.Considering exponential solutions in the formv=eλt+βxw, then(79)Lv=0where(80)L=cmλ2-amβ2aβ-I0λ2β+I1β3-aβemλ2-bmβ2+a.In order to guaranteew nonzero, we must have det(L)=0. Thus we obtain the characteristic equation(81)Δβ,λ=cmemλ4+aI0-bmcm-amemβ2+acmλ2+a2-aamβ2+ambm-aI1β4=0.


### 7.1. Singular Eigenvalue Problem

By assuming exponential solutionsv(t,x)=eλtw(x),w(x)=W(x)Ψ(x) in (74), it follows the eigenvalue problem(82)λ2M+Kwx=0,which can be written as a singular third-order differential matrix equation(83)Rw′′′x+Mw′′x+Cλw′x+Kλwx=0,with(84)R=00I10,M=-am00-bm,Cλ=0a-λ2I0+a0,Kλ=λ2cm00λ2em+a.We observe that the nature of the matrix R characterizes the problem (83) as being singular. Following the spectral methodology, we seek exponential solutions w(x)=eβxv, v≠0, of (83). It follows that(85)Γβ,λv=β3R+β2M+βCλ+Kλv=0,where(86)Γβ,λ=-β2amcmλ2βaβ3I1-βI0λ2+a-β2bm+emλ2+a.For each λ, non-zero solutions v can be obtained whenever β is a root of the characteristic polynomial(87)Pβ,λ=det⁡β3R+β2M+βCλ+Kλ.The characteristic equation(88)Δβ,λ=Pβ,λ=β4+g2λβ2-r4λ=0,with(89)g2λ=b2b0=I0a-bmcm-amemλ2+a2-aamambm-aI1,r4λ=-b4b0=-λ2acm+cmemλ2ambm-aI1has the roots β1,2=±ε and β3,4=±iδ, where(90)ε=-g2λ+g2λ2+4r4λ2,δ=g2λ+g2λ2+4r4λ2.The case of simple roots can be characterized from the conditionP′=2β3+g2(λ),β≠0. This excludes the possibility of β=0 to be a simple root. Moreover, we can determine the eigenvectors satisfying (85). We summarize these results in the following.Lemma 1.
β is a simple root of P(β,λ) if and only if β≠0 and β2≠-g2(λ)/2. For simple roots, the eigenvectors satisfying (85) are given by(91)vj=v1jβj2am-λ2cmβjav1j,where v1j is an arbitrary constant.
For the case of distinct rootsβ, the general solution(92)wx=c1eβ1tv1+c2eβ2tv2+c3eβ2tv3+c4eβ4tv4=c1c2c3c4c1ε2am-λ2cmaεc2ε2am-λ2cm-aεc3i-δ2am-λ2cmaδc4i-δ2am-λ2cm-aδeεxe-εxeδixe-δix,and can be also presented in the form(93)wx=C1C2C3C4C2KaC1KaC4Kb-C3Kbcosh⁡εxsinh⁡εxcos⁡δxsin⁡δx,where C1=c1+c2, C2=c1-c2, C3=c3+c4, C4=(c3-c4)i and(94)Ka=amε2-λ2cmaε,Kb=amδ2-λ2cmaδ.The above representation ofw(x) in terms of hyperbolic and trigonometric functions follows by introducing in (92) the basic change(95)eεxe-εxeiδxe-iδx=11001-100001i001-icosh⁡εxsinh⁡εxcos⁡δxsin⁡δx.


#### 7.1.1. Repeated Roots Analysis

From the given characterization of simple roots, it is clear that the repeated roots can be onlyβ=0 and β2=-g2(λ)/2. The conditions P′′′=12β=0 and P(iv)=12 exclude the possibility of having nonzero third roots. In such a case, β=0 will be a quadruple root. More precisely, we have the following.Lemma 2.
( 1 ) If β=0 is a root of P(β,λ), then it is a double root if g2(λ)≠0, or quadruple if g2(λ)=0.( 2 ) If β=-g2(λ)/2 and β=--g2(λ)/2 are roots, with g2(λ)≠0, then they will be double roots.( 3 ) The characteristic polynomial P(β,λ) does not have triple roots.From the above lemmas, critical frequencies values will arise for natural frequenciesλ=iω.(i)β=0 is a root of P(β,iω), if and only if r4(ω)=0. That is, r4(iω)=0 for ω=0 or ω=ωc, with ωc=a/em.(ii)β=±-g2(ω)/2, g2(ω)≠0, is a root of P(β,iω) if and only if (g2(iω))2+4r4(ω)=0, with r4(ω)≠0. We have g2(ω)=0 when ω=ωg=(a2-aam)/(bmcm+amem-I0a).


### 7.2. Frequency Equations

Frequency equations will arise when applying boundary conditions as given below.


#### 7.2.1. Simply Supported Beam

In this case we have the following boundary conditions, atx=0(96)w0=0,ω2I0w0-bmϕ′0+I1w′′0=0,and at x=L(97)wL=0,ω2I0wL-bmϕ′L+I1w′′L=0.Using the first boundary condition at x=0, we obtain C2+C4=0⇒C4=-C2. From the second boundary condition we have(98)I1w′′0+bmϕ′0=0,-C2I1δ2+ε2+C2bmKaδ+Kbε=0,replacing Ka and Kb and using that ambm-I1a>0 [13], ε2+δ2≠0, results in C2=C4=0.From the boundary conditions atx=L, written in matrix form(99)100ω2I0I1-bmWLW′′Lψ′L=0the system Uc=0 is obtained, described by(100)100ω2I0I1-bmsin⁡δLsinh⁡εL-sin⁡δLδ2ε2sinh⁡εL-Kasin⁡δLδεKbsinh⁡εLC1C3=0.Sincebmam-I1a>0 and ε2+δ2≠0, the frequency equation (101)Δω=det⁡U=I1a-bmamε2+δ2asinh⁡εLsin⁡δL=0reduces to(102)sinh⁡εLsin⁡δL=0.


#### 7.2.2. Cantilever Beam

In this case the boundary conditions are given by(103)w0=0ψ0=0,amw′L-aψL=0-bmψ′L+I1w′′L+I0ω2wL=0.Similarly to the previous case, by using the boundary conditions at x=0, we obtain(104)Wx=C1sin⁡δx-KaKbsinh⁡εx+C2cos⁡δx-cosh⁡εx,Ψx=C1Kacos⁡δx-cos⁡εx-C2Kasin⁡δx+Kbsinh⁡εx,The boundary conditions in matrix form at x=L(105)0am0-a0ω2I00I10-bmWLW′LW′′LψLψ′LC1C2=0lead to a system Uc=0 from which we obtain the characteristic equation.In Table3 and Figure 8, we present numerical results about natural frequencies and mode shapes. Both cases of supported-supported (S-S) and cantilever (C) boundary conditions are considered. The simulation was done with a silicon beam having the following material and surface parameters [13]: (106)E=107GPa,ν=0.33,ρ=2.33×103kg/m3,ρ0=3.17×10-7kg/m2,μ0=-2.7779N/m,λ0=-4.4939N/m,τ0=0.6056N/m.Table 3
Comparison between natural frequencies of Timoshenko beam models including and not including surface effects, considering nano- and microscale, of a Silicon beam. Nanobeam geometric parameters:L=50×10-9m, H=2h=6×10-9m, and b=3×10-9m. Microbeam geometric parameters: L=50×10-6m, H=2h=6×10-6m, and b=3×10-6m.
Nano Micro BC Model Natural frequenciesωi(GHz) Natural frequenciesωi(MHz) ω 1 ω 2 ω 3 ω 4 ω 1 ω 2 ω 3 ω 4 S-S TMB 7.0808 25.063 51.327 82.922 7.1996 27.010 55.669 89.756 TB 7.1998 27.012 55.673 89.763 7.1998 27.012 55.673 89.763 TMB [13] 7.08 25.07 51.33 82.92 - - - - TB [13] 7.20 27.02 55.70 89.76 - - - -  C TMB 3.0168 14.674 36.484 64.514 2.5976 15.273 39.297 69.725 TB 2.597 15.273 39.300 69.730 2.5971 15.273 39.300 69.730 TMB [13] 3.02 14.67 36.38 64.46 - - - - TB [13] 2.60 15.28 39.31 69.74 - - - -Figure 8
Mode shapes of simply supported beam (first column) and cantilever beam (second column), classical case (solid line) and including surface effects (dashed line).Comments(i) From Table3 and Figure 8, it is observed that the results agree with those found in the literature [13, 14].(ii) When the beam length increases, from nano- to microscale, the results tend to the classical results.(iii) Natural frequencies can increase or decrease compared to classical results, depending on elastic surface constants signs and boundary conditions (cantilever beam is more influenceable than simply supported beam).(iv) According to Figure8, the mode shape component W(x) of a simply supported beam with the inclusion of surface effects is the same of a classical beam. Perceptible differences occur for the cantilever beam case, principally for the first one mode shape, and just in the nanoscale.


## 7.1. Singular Eigenvalue Problem

By assuming exponential solutionsv(t,x)=eλtw(x),w(x)=W(x)Ψ(x) in (74), it follows the eigenvalue problem(82)λ2M+Kwx=0,which can be written as a singular third-order differential matrix equation(83)Rw′′′x+Mw′′x+Cλw′x+Kλwx=0,with(84)R=00I10,M=-am00-bm,Cλ=0a-λ2I0+a0,Kλ=λ2cm00λ2em+a.We observe that the nature of the matrix R characterizes the problem (83) as being singular. Following the spectral methodology, we seek exponential solutions w(x)=eβxv, v≠0, of (83). It follows that(85)Γβ,λv=β3R+β2M+βCλ+Kλv=0,where(86)Γβ,λ=-β2amcmλ2βaβ3I1-βI0λ2+a-β2bm+emλ2+a.For each λ, non-zero solutions v can be obtained whenever β is a root of the characteristic polynomial(87)Pβ,λ=det⁡β3R+β2M+βCλ+Kλ.The characteristic equation(88)Δβ,λ=Pβ,λ=β4+g2λβ2-r4λ=0,with(89)g2λ=b2b0=I0a-bmcm-amemλ2+a2-aamambm-aI1,r4λ=-b4b0=-λ2acm+cmemλ2ambm-aI1has the roots β1,2=±ε and β3,4=±iδ, where(90)ε=-g2λ+g2λ2+4r4λ2,δ=g2λ+g2λ2+4r4λ2.The case of simple roots can be characterized from the conditionP′=2β3+g2(λ),β≠0. This excludes the possibility of β=0 to be a simple root. Moreover, we can determine the eigenvectors satisfying (85). We summarize these results in the following.Lemma 1.
β is a simple root of P(β,λ) if and only if β≠0 and β2≠-g2(λ)/2. For simple roots, the eigenvectors satisfying (85) are given by(91)vj=v1jβj2am-λ2cmβjav1j,where v1j is an arbitrary constant.
For the case of distinct rootsβ, the general solution(92)wx=c1eβ1tv1+c2eβ2tv2+c3eβ2tv3+c4eβ4tv4=c1c2c3c4c1ε2am-λ2cmaεc2ε2am-λ2cm-aεc3i-δ2am-λ2cmaδc4i-δ2am-λ2cm-aδeεxe-εxeδixe-δix,and can be also presented in the form(93)wx=C1C2C3C4C2KaC1KaC4Kb-C3Kbcosh⁡εxsinh⁡εxcos⁡δxsin⁡δx,where C1=c1+c2, C2=c1-c2, C3=c3+c4, C4=(c3-c4)i and(94)Ka=amε2-λ2cmaε,Kb=amδ2-λ2cmaδ.The above representation ofw(x) in terms of hyperbolic and trigonometric functions follows by introducing in (92) the basic change(95)eεxe-εxeiδxe-iδx=11001-100001i001-icosh⁡εxsinh⁡εxcos⁡δxsin⁡δx.


### 7.1.1. Repeated Roots Analysis

From the given characterization of simple roots, it is clear that the repeated roots can be onlyβ=0 and β2=-g2(λ)/2. The conditions P′′′=12β=0 and P(iv)=12 exclude the possibility of having nonzero third roots. In such a case, β=0 will be a quadruple root. More precisely, we have the following.Lemma 2.
( 1 ) If β=0 is a root of P(β,λ), then it is a double root if g2(λ)≠0, or quadruple if g2(λ)=0.( 2 ) If β=-g2(λ)/2 and β=--g2(λ)/2 are roots, with g2(λ)≠0, then they will be double roots.( 3 ) The characteristic polynomial P(β,λ) does not have triple roots.From the above lemmas, critical frequencies values will arise for natural frequenciesλ=iω.(i)β=0 is a root of P(β,iω), if and only if r4(ω)=0. That is, r4(iω)=0 for ω=0 or ω=ωc, with ωc=a/em.(ii)β=±-g2(ω)/2, g2(ω)≠0, is a root of P(β,iω) if and only if (g2(iω))2+4r4(ω)=0, with r4(ω)≠0. We have g2(ω)=0 when ω=ωg=(a2-aam)/(bmcm+amem-I0a).


## 7.1.1. Repeated Roots Analysis

From the given characterization of simple roots, it is clear that the repeated roots can be onlyβ=0 and β2=-g2(λ)/2. The conditions P′′′=12β=0 and P(iv)=12 exclude the possibility of having nonzero third roots. In such a case, β=0 will be a quadruple root. More precisely, we have the following.Lemma 2.
( 1 ) If β=0 is a root of P(β,λ), then it is a double root if g2(λ)≠0, or quadruple if g2(λ)=0.( 2 ) If β=-g2(λ)/2 and β=--g2(λ)/2 are roots, with g2(λ)≠0, then they will be double roots.( 3 ) The characteristic polynomial P(β,λ) does not have triple roots.From the above lemmas, critical frequencies values will arise for natural frequenciesλ=iω.(i)β=0 is a root of P(β,iω), if and only if r4(ω)=0. That is, r4(iω)=0 for ω=0 or ω=ωc, with ωc=a/em.(ii)β=±-g2(ω)/2, g2(ω)≠0, is a root of P(β,iω) if and only if (g2(iω))2+4r4(ω)=0, with r4(ω)≠0. We have g2(ω)=0 when ω=ωg=(a2-aam)/(bmcm+amem-I0a).


## 7.2. Frequency Equations

Frequency equations will arise when applying boundary conditions as given below.


### 7.2.1. Simply Supported Beam

In this case we have the following boundary conditions, atx=0(96)w0=0,ω2I0w0-bmϕ′0+I1w′′0=0,and at x=L(97)wL=0,ω2I0wL-bmϕ′L+I1w′′L=0.Using the first boundary condition at x=0, we obtain C2+C4=0⇒C4=-C2. From the second boundary condition we have(98)I1w′′0+bmϕ′0=0,-C2I1δ2+ε2+C2bmKaδ+Kbε=0,replacing Ka and Kb and using that ambm-I1a>0 [13], ε2+δ2≠0, results in C2=C4=0.From the boundary conditions atx=L, written in matrix form(99)100ω2I0I1-bmWLW′′Lψ′L=0the system Uc=0 is obtained, described by(100)100ω2I0I1-bmsin⁡δLsinh⁡εL-sin⁡δLδ2ε2sinh⁡εL-Kasin⁡δLδεKbsinh⁡εLC1C3=0.Sincebmam-I1a>0 and ε2+δ2≠0, the frequency equation (101)Δω=det⁡U=I1a-bmamε2+δ2asinh⁡εLsin⁡δL=0reduces to(102)sinh⁡εLsin⁡δL=0.


### 7.2.2. Cantilever Beam

In this case the boundary conditions are given by(103)w0=0ψ0=0,amw′L-aψL=0-bmψ′L+I1w′′L+I0ω2wL=0.Similarly to the previous case, by using the boundary conditions at x=0, we obtain(104)Wx=C1sin⁡δx-KaKbsinh⁡εx+C2cos⁡δx-cosh⁡εx,Ψx=C1Kacos⁡δx-cos⁡εx-C2Kasin⁡δx+Kbsinh⁡εx,The boundary conditions in matrix form at x=L(105)0am0-a0ω2I00I10-bmWLW′LW′′LψLψ′LC1C2=0lead to a system Uc=0 from which we obtain the characteristic equation.In Table3 and Figure 8, we present numerical results about natural frequencies and mode shapes. Both cases of supported-supported (S-S) and cantilever (C) boundary conditions are considered. The simulation was done with a silicon beam having the following material and surface parameters [13]: (106)E=107GPa,ν=0.33,ρ=2.33×103kg/m3,ρ0=3.17×10-7kg/m2,μ0=-2.7779N/m,λ0=-4.4939N/m,τ0=0.6056N/m.Table 3
Comparison between natural frequencies of Timoshenko beam models including and not including surface effects, considering nano- and microscale, of a Silicon beam. Nanobeam geometric parameters:L=50×10-9m, H=2h=6×10-9m, and b=3×10-9m. Microbeam geometric parameters: L=50×10-6m, H=2h=6×10-6m, and b=3×10-6m.
Nano Micro BC Model Natural frequenciesωi(GHz) Natural frequenciesωi(MHz) ω 1 ω 2 ω 3 ω 4 ω 1 ω 2 ω 3 ω 4 S-S TMB 7.0808 25.063 51.327 82.922 7.1996 27.010 55.669 89.756 TB 7.1998 27.012 55.673 89.763 7.1998 27.012 55.673 89.763 TMB [13] 7.08 25.07 51.33 82.92 - - - - TB [13] 7.20 27.02 55.70 89.76 - - - -  C TMB 3.0168 14.674 36.484 64.514 2.5976 15.273 39.297 69.725 TB 2.597 15.273 39.300 69.730 2.5971 15.273 39.300 69.730 TMB [13] 3.02 14.67 36.38 64.46 - - - - TB [13] 2.60 15.28 39.31 69.74 - - - -Figure 8
Mode shapes of simply supported beam (first column) and cantilever beam (second column), classical case (solid line) and including surface effects (dashed line).Comments(i) From Table3 and Figure 8, it is observed that the results agree with those found in the literature [13, 14].(ii) When the beam length increases, from nano- to microscale, the results tend to the classical results.(iii) Natural frequencies can increase or decrease compared to classical results, depending on elastic surface constants signs and boundary conditions (cantilever beam is more influenceable than simply supported beam).(iv) According to Figure8, the mode shape component W(x) of a simply supported beam with the inclusion of surface effects is the same of a classical beam. Perceptible differences occur for the cantilever beam case, principally for the first one mode shape, and just in the nanoscale.


## 7.2.1. Simply Supported Beam

In this case we have the following boundary conditions, atx=0(96)w0=0,ω2I0w0-bmϕ′0+I1w′′0=0,and at x=L(97)wL=0,ω2I0wL-bmϕ′L+I1w′′L=0.Using the first boundary condition at x=0, we obtain C2+C4=0⇒C4=-C2. From the second boundary condition we have(98)I1w′′0+bmϕ′0=0,-C2I1δ2+ε2+C2bmKaδ+Kbε=0,replacing Ka and Kb and using that ambm-I1a>0 [13], ε2+δ2≠0, results in C2=C4=0.From the boundary conditions atx=L, written in matrix form(99)100ω2I0I1-bmWLW′′Lψ′L=0the system Uc=0 is obtained, described by(100)100ω2I0I1-bmsin⁡δLsinh⁡εL-sin⁡δLδ2ε2sinh⁡εL-Kasin⁡δLδεKbsinh⁡εLC1C3=0.Sincebmam-I1a>0 and ε2+δ2≠0, the frequency equation (101)Δω=det⁡U=I1a-bmamε2+δ2asinh⁡εLsin⁡δL=0reduces to(102)sinh⁡εLsin⁡δL=0.


## 7.2.2. Cantilever Beam

In this case the boundary conditions are given by(103)w0=0ψ0=0,amw′L-aψL=0-bmψ′L+I1w′′L+I0ω2wL=0.Similarly to the previous case, by using the boundary conditions at x=0, we obtain(104)Wx=C1sin⁡δx-KaKbsinh⁡εx+C2cos⁡δx-cosh⁡εx,Ψx=C1Kacos⁡δx-cos⁡εx-C2Kasin⁡δx+Kbsinh⁡εx,The boundary conditions in matrix form at x=L(105)0am0-a0ω2I00I10-bmWLW′LW′′LψLψ′LC1C2=0lead to a system Uc=0 from which we obtain the characteristic equation.In Table3 and Figure 8, we present numerical results about natural frequencies and mode shapes. Both cases of supported-supported (S-S) and cantilever (C) boundary conditions are considered. The simulation was done with a silicon beam having the following material and surface parameters [13]: (106)E=107GPa,ν=0.33,ρ=2.33×103kg/m3,ρ0=3.17×10-7kg/m2,μ0=-2.7779N/m,λ0=-4.4939N/m,τ0=0.6056N/m.Table 3
Comparison between natural frequencies of Timoshenko beam models including and not including surface effects, considering nano- and microscale, of a Silicon beam. Nanobeam geometric parameters:L=50×10-9m, H=2h=6×10-9m, and b=3×10-9m. Microbeam geometric parameters: L=50×10-6m, H=2h=6×10-6m, and b=3×10-6m.
Nano Micro BC Model Natural frequenciesωi(GHz) Natural frequenciesωi(MHz) ω 1 ω 2 ω 3 ω 4 ω 1 ω 2 ω 3 ω 4 S-S TMB 7.0808 25.063 51.327 82.922 7.1996 27.010 55.669 89.756 TB 7.1998 27.012 55.673 89.763 7.1998 27.012 55.673 89.763 TMB [13] 7.08 25.07 51.33 82.92 - - - - TB [13] 7.20 27.02 55.70 89.76 - - - -  C TMB 3.0168 14.674 36.484 64.514 2.5976 15.273 39.297 69.725 TB 2.597 15.273 39.300 69.730 2.5971 15.273 39.300 69.730 TMB [13] 3.02 14.67 36.38 64.46 - - - - TB [13] 2.60 15.28 39.31 69.74 - - - -Figure 8
Mode shapes of simply supported beam (first column) and cantilever beam (second column), classical case (solid line) and including surface effects (dashed line).Comments(i) From Table3 and Figure 8, it is observed that the results agree with those found in the literature [13, 14].(ii) When the beam length increases, from nano- to microscale, the results tend to the classical results.(iii) Natural frequencies can increase or decrease compared to classical results, depending on elastic surface constants signs and boundary conditions (cantilever beam is more influenceable than simply supported beam).(iv) According to Figure8, the mode shape component W(x) of a simply supported beam with the inclusion of surface effects is the same of a classical beam. Perceptible differences occur for the cantilever beam case, principally for the first one mode shape, and just in the nanoscale.


## 8. Multispan Timoshenko Beams

The same formulation, which was presented in Section2 for Euler-Bernoulli multispan beams, can be used with the Timoshenko multispan beam model. We replace the scalar components that appear in the block matrices defined after (1), by the 2×2 matrix components Mj and Kj(107)Mj=ρjAj00ρjIj,Kj=-κjGjAj∂2∂x2κjGjAj∂∂x-κjGjAj∂∂x-EjIj∂2∂x2+κjGjAj.The components vj become 2×1 vectors, vj=wj(t,x)ϕj(t,x), where wj(t,x) is the flexural displacement and ϕj(t,x) is the angular displacement. Similarly, the components of Fj are external loads and moments, respectively, fj and qj. Then(108)v=w1ϕ1w2ϕ2⋮wNϕN,F=f1q1f2q2⋮fNqN.General boundary conditions can be written as(109)B0w1t,0=n0,BLwNt,L=nL,where(110)B0=a11b11c11d11a12b12c12d12,n0=n1n2,BL=a21b21c21d21a22b22c22d22,nL=n3n4,where nj corresponds to devices or interaction of extremities with the external medium, when they occur, and(111)wjt,x=wjt,xϕjt,xwj′t,xϕj′t,x,j=1:N.The continuity conditions for displacement, gyro, flexural moment, and shear force in the transversal section discontinuity pointsx=xj, j=1:N-1, including intermediate devices can be written as(112)C1,iwit,xi=C2,iwi+1t,xi,i=1:N-1,where(113)C1,i=e11if11ig11ih11ie21if21ig21ih21ie31if31ig31ih31ie41if41ig41ih41i,C2,i=e12if12ig12ih12ie22if22ig22ih22ie32if32ig32ih32ie42if42ig42ih42i,and the vector wi defined in (111).In what follows, we shall consider double-span beams in which the first segment is modeled by the Timoshenko beam model and the second segment can be modeled by a Timoshenko (TB) or an Euler-Bernoulli (EB) model.TB-TB. In this case, wj(t,x) and ϕj(t,x) satisfy for j=1,2 the equations(114)ρjAj∂2wjt,x∂t2-κjGjAj∂2wjt,x∂x2-∂ϕjt,x∂x=0,ρjIj∂2ϕjt,x∂t2-EjIj∂2ϕjt,x∂x2-κjGjAj∂wjt,x∂x-ϕjt,x=0.For a simply supported beam we have(115)w1t,0=0,ϕ1′t,0=0,w2t,L=0,ϕ2′t,L=0,and for a cantilever beam(116)w1t,0=0,ϕ1t,0=0,ϕ2′t,L=0,w2′t,L-ϕ2t,L=0.The compatibility conditions, at the pointx=x1, without including intermediate devices, are given by(117)w1t,x1=w2t,x1,ϕ1t,x1=ϕ2t,x1,E1I1ϕ1′t,x1=E2I2ϕ2′t,x1,κ1G1A1w1′t,x1-ϕ1t,x1=κ2G2A2w2′t,x1-ϕ2t,x1.The spatial amplitude of a free transversal vibration, at each segment,vj(t,x)=eλtXj(x), where Xj(x)=Wj(x)Ψj(x), for j=1,2, satisfies(118)λ2Mj+λCjXjx+KjXjx=0,or in a differential form(119)MjXj′′+CjXj′+KjλXj=0,with(120)Mj=-κjGjAj00-EjIj,Cj=0κjGjAj-κjGjAj0,Kj=λ2ρjAj00ρjIjλ2+κGA.The boundary conditions, obtained from (109), can be conveniently written as(121)B0X10=0,BLX2L=0,and the compatibility conditions, obtained from (112), at x=x1(122)C1,1X1x1=C2,1X2x1,Xjx=XjxXj′x.where(123)Xx=X1x,0≤x≤x1,X2x,x1<x≤L.We choose the fundamental basis at each segment, considering convenient translations,(124)Ψ1x=h1x,h1′x,Ψ2x=h2L-x,h2′L-x,with hj(x) being the 2×2 matrix solution of the initial-value problem(125)Mjhj′′x+Cjλhj′x+Kjλhjx=0,j=1,…,N.Mjhj′0=I,hj0=0.The matrix fundamental solutionhj is given by [25](126)hjx=ajdjx+bjdj′′x-sjdjxsjdjxrjdjx-sjdj′′x,where(127)aj=κjGjAj+λ2ρjIj,bj=-EjIj,rj=λ2ρjAj,sj=κjGjAj.Then(128)X1x=h1xc11+h1′xc12,X2x=h2x-x1c21+h2′x-x1c22.Replacing (128) in the boundary and compatibility conditions we obtained a system Uc=0, with U=BΦ,(129)B=B00000C1,1-C2,10000BL8×16.For a cantilever beam we have(130)B0=10000100,BL=00010-110.The matrixΦ is composed of the elements of the solution basis(131)Φ=Φ100Φ1l100Φ2l10Φ2L16×8,where(132)Φjx=hjx-ajhj′x-ajhj′x-ajhj′′x-aj,j=1,2,a1=0,a2=L.TB-EB. We now consider that the movement in the second segment is governed by the Euler-Bernoulli beam model (TB-EB). For simplicity, we shall consider Euler-Bernoulli compatibility conditions at the discontinuity point.In this case, for0≤x≤x1 the transversal displacement wT and angular displacement ϕT are governed by the Timoshenko equations(133)ρ1A1∂2wTt,x∂t2-κ1G1A1∂2wTt,x∂x2-∂ϕTt,x∂x=0,ρ1I1∂2ϕTt,x∂t2-E1I1∂2ϕTt,x∂x2-κ1G1A1∂wTt,x∂x-ϕTt,x=0,while for x1≤x≤L, the displacement wE follows the Euler-Bernoulli model(134)ρA2∂2wEt,x∂t2+EI2∂4wEx∂x4=0.The problem can be also described in matrix form asMv¨+Kv=0, where the block matrix coefficients are(135)M=M1T00M2E,K=K1T00K2E,v=v1Tv2E=wTt,xϕTt,xwEt,x,with M1T, C1T, and K1T being the 2×2 matrix coefficients corresponding to Timoshenko equations and M2E, K2E are scalar components corresponding to Euler-Bernoulli equations.The general boundary conditions are given inx=0 according to Timoshenko model, in terms of the state vector wT. And in x=L the general boundary conditions are given as Euler-Bernoulli model; that is,(136)B0wTt,0=0,wTt,x=wTt,xϕTt,xwT′t,xϕT′t,x,BLwEt,L=0,wEt,x=wEt,xwE′t,xwE′′t,xwE′′′t,x.In the discontinuity pointx=x1, point of section transversal area change, we have chosen compatibility conditions referring to Euler-Bernoulli beam model, having(137)wTt,x1=wEt,x1,wT′t,x1=wE′t,x1,E1I1wT′′t,x1=E2I2wE′′t,x1,E1I1wT′′′t,x1=E2I2wE′′′t,x1.The eigensolutions are of the type v=eλtX(x) with(138)Xx=X1TxX2Ex,X1T=WTxΨTx,X2E=WEx.The eigenvalue problem (λ2M+K)X(x)=0 is equivalent to the differential system(139)MXT′′x+CXT′x+KλXTx=0,0≤x≤x1,XEivx-ε24XEx=0,x1<x≤L,whereM, C, and K were defined in (120) and ε2 was described in (25).The boundary conditions can be written as(140)B~0X~T0=0,BLX~E0=0,where(141)X~Tx=WTxΨTxWT′xΨT′xWT′′xWT′′′x,X~E=WExWE′xWE′′xWE′′′x.The compatibility conditions (137) can be written as(142)100000-10000010000-10000001000-α0000001000-αX~x1=0,with respect to the state span vector X~(x)=X~T(x)X~E(x).The solutionsXT(x) and XE(x) can be written using the solutions basis, which were introduced in the previous sections:(143)XTx=WTxΨTx=hxc1c2+h′xc3c4,XEx=dL-xc5+d′L-xc6+d′′L-xc7+d′′′L-xc8.By replacing (143) in the boundary conditions and compatibility equations, we obtain a system Uc=0, c=[c1c2c3c4c5c6c7c8]T and the characteristic equation Δ(λ)=det(U)=0. As before, we decompU as U=Bϕ, where(144)ϕ=Φ100Φ1x100Φ2x10Φ2L20×8,with(145)ϕ1=hx2×2h′x2×2h′x2×2h′′x2×2h~x2×2h~′x2×26×4,h~x=h′′x1,1h′′x1,2h′′′x1,1h′′′x1,2,ϕ2=dL-xd′L-xd′′L-xd2′′′L-xd′L-xd′′L-xd′′′L-xd′′′′L-xd′′L-xd′′′L-xdivL-xdvL-xd′′′L-xdivL-xdvL-xd2viL-x4×4,(146)B=B~02×6000Bx14×10000BL2×4.


### 8.1. Simulations

The two span beam models considered have been compared with results of Euler-Bernoulli and Timoshenko double-span beam models in dimensional and unidimensional formulations.Dimensional Model. We assume a double-span beam, with transversal section discontinuity due to variable thickness. This beam has constant width, but different thickness. The geometric and materials properties of this beam are presented in Table 4.Table 4
Parameters of the double-span beam [29].
Parameters Numeric value Unity Lengthl1 0,254 m Lengthl2 0,140 m Widthw 0,02545 m Thicknesst1 0,01905 m Thicknesst2 0,00549 m Mass densityρ 2830 K g / m 3 Young’s modulusE 71,7 G P a Moment of inertiaI1 1 , 4633 × 1 0 - 8 m 4 Moment of inertiaI2 3 , 502 × 1 0 - 10 m 4We present results about natural frequencies, of three cases: two Euler-Bernoulli segments (EB-EB), two Timoshenko segments (TB-TB), and coupled case (TB-EB), in Table5.Table 5
Comparison of natural frequencies from the double-span beam. EB-EB: two segments of Euler-Bernoulli beam, TB-TB: two segments of Timoshenko beam, and TB-EB: first segment using Euler-Bernoulli beam model and second segment using Timoshenko beam model.
Case Model Natural frequenciesωi (Hz) ω 1 ω 2 ω 3 ω 4 S-S EB-EB 101.77 782.09 1300.56 2874.69 TB-EB 101.74 777.85 1291.65 2816.98 TB-TB 101.64 775.85 1285.65 2796.89  C-F EB-EB 139.56 326.69 1170.27 1797.41 TB-EB 139.39 325.84 1154.90 1771.78 TB-TB 139.31 325.45 1150.82 1761.06  F-F EB-EB 292.42 1181.28 1804.01 3603.04 TB-EB 292.14 1775.94 1775.8 3512.95 TB-TB 291.77 1167.89 1786.63 3483.40 [Ref] 286-291 1159-1165 1759-1771 -We observe from Table5 that the values of natural frequencies of the coupled beam (TB-EB) are between the respective values obtained by TB-TB and EB-EB cases, wherein the EB-EB natural frequencies, as expected, are greater than the others, in all boundary conditions cases presented: S-S (simply supported), C-F (clamped free), and F-F (free-free). Moreover, in the case of free-free boundary condition, as expected, TB-TB produces the most similar results to those obtained experimentally in [29].Dimensionless Model. We now consider the two span models with classical dimensionless variables. For the Euler-Bernoulli beam model, the spatial problem in a dimensionless form is given by(147)d4X¯jXjdXj-Ωj4X¯jXj=0,where the dimensionless parameters areXj=xj/L,X¯j(Xj)=Xj/L,Ωj4=ρjAjω2/EjIj,Ωj4=μjΩ14.The Timoshenko model in a dimensionless matrix form is given by(148)ξj001Yj′′XΨj′′X+0-ξjξj1Yj′XΨj′X+Ωj200-ξj+ηjΩj2YjXΨjX=00,with the dimensionless parameters defined by(149)X=xL;YjX=YjxL;ΨjX=ψjx;Ωj2=ρAjω2l4EIj;ηj=IjAjL2;ξj=κGjAjL2EIj=κ2ηj1+νj.Then we have the dimensionless frequencies by each segment(150)Ωj2=ρAjω2L4EIjTimoshenko,Ωj4=ρAjω2L4EIjEuler-Bernoulli.These frequencies can be written as a relation ofΩ1. Ωj=(μj)1/2Ω1, for Timoshenko beam, and Ωj=μj1/4Ω1 for Euler-Bernoulli beam, being μj=(ρA)j/(ρA)1(EI)1/(EI)j.In Table6 we present results about dimensionless natural frequencies; we have considered variations in the position of the discontinuity point Rj=xj/L and variations of thickness of the second segment; constant width is considered, according to the following configurations.Table 6
The effects of the change of discontinuity pointR1=x1/L in the dimensionless natural frequencies for a cantilever double span beam with variations of thickness.
Configuration A Configuration B Case Model Dimensionless natural frequenciesΩ[i] Ω [ 1 ] Ω [ 2 ] Ω [ 3 ] Ω [ 4 ] Ω [ 1 ] Ω [ 2 ] Ω [ 3 ] Ω [ 4 ] R 1 = 0.2 TB-TB 1.4093 3.4619 5.6243 7.4545 1.7184 4.0974 6.5148 8.7633 TB-EB 1.4098 3.4721 5.6637 7.5438 1.7193 4.1198 6.5990 8.9637 EB-EB 1.4101 3.4769 5.6854 7.6159 1.7209 4.1381 6.6554 9.077  R 1 = 0.5 TB-TB 1.9972 3.5642 5.9015 8.2377 2.0150 4.0913 6.8545 9.1581 TB-EB 1.9985 3.5725 5.9372 8.3063 2.0162 4.1061 6.9121 9.3086 EB-EB 2.0004 3.5901 5.9831 8.5546 2.0198 4.1331 7.050 9.5862  R 2 = 0.8 TB-TB 2.0982 4.7587 6.6143 9.1855 1.9918 4.7325 7.3808 9.7778 TB-EB 2.1000 4.7673 6.6436 9.2371 1.9933 4.7395 7.4083 9.8619 EB-EB 2.1046 4.8267 6.7614 9.6670 1.9973 4.8113 7.6241 10.329Configuration A: double-span beam with variable thickness and constant width, t1=0.08, t2=0.03, w1=w2=0.03.Configuration B: double-span beam with variable thickness and constant width, t1=0.08, t2=0.05, w1=w2=0.03.Comparative results of dimensionless natural frequencies are presented in Table6.We observe from Table6 that the values of dimensionless natural frequencies of the coupled beam (TB-EB) are between the respective values obtained by TB-TB and EB-EB cases, wherein the EB-EB natural frequencies, as expected, are greater than the others. Also, we can observe that when the thickness of the second segment increases, the natural frequencies become greater, mainly the higher frequencies (third and fourth frequencies), except for the case R1=0.8 in which the first and second frequencies become greater. When R1 goes from 0.2 to 0.5, all frequencies increase, but when R1 goes from 0.5 to 0.8, this occurs only for Configuration A, and for Configuration B the first one frequency does not increase in this case.


## 8.1. Simulations

The two span beam models considered have been compared with results of Euler-Bernoulli and Timoshenko double-span beam models in dimensional and unidimensional formulations.Dimensional Model. We assume a double-span beam, with transversal section discontinuity due to variable thickness. This beam has constant width, but different thickness. The geometric and materials properties of this beam are presented in Table 4.Table 4
Parameters of the double-span beam [29].
Parameters Numeric value Unity Lengthl1 0,254 m Lengthl2 0,140 m Widthw 0,02545 m Thicknesst1 0,01905 m Thicknesst2 0,00549 m Mass densityρ 2830 K g / m 3 Young’s modulusE 71,7 G P a Moment of inertiaI1 1 , 4633 × 1 0 - 8 m 4 Moment of inertiaI2 3 , 502 × 1 0 - 10 m 4We present results about natural frequencies, of three cases: two Euler-Bernoulli segments (EB-EB), two Timoshenko segments (TB-TB), and coupled case (TB-EB), in Table5.Table 5
Comparison of natural frequencies from the double-span beam. EB-EB: two segments of Euler-Bernoulli beam, TB-TB: two segments of Timoshenko beam, and TB-EB: first segment using Euler-Bernoulli beam model and second segment using Timoshenko beam model.
Case Model Natural frequenciesωi (Hz) ω 1 ω 2 ω 3 ω 4 S-S EB-EB 101.77 782.09 1300.56 2874.69 TB-EB 101.74 777.85 1291.65 2816.98 TB-TB 101.64 775.85 1285.65 2796.89  C-F EB-EB 139.56 326.69 1170.27 1797.41 TB-EB 139.39 325.84 1154.90 1771.78 TB-TB 139.31 325.45 1150.82 1761.06  F-F EB-EB 292.42 1181.28 1804.01 3603.04 TB-EB 292.14 1775.94 1775.8 3512.95 TB-TB 291.77 1167.89 1786.63 3483.40 [Ref] 286-291 1159-1165 1759-1771 -We observe from Table5 that the values of natural frequencies of the coupled beam (TB-EB) are between the respective values obtained by TB-TB and EB-EB cases, wherein the EB-EB natural frequencies, as expected, are greater than the others, in all boundary conditions cases presented: S-S (simply supported), C-F (clamped free), and F-F (free-free). Moreover, in the case of free-free boundary condition, as expected, TB-TB produces the most similar results to those obtained experimentally in [29].Dimensionless Model. We now consider the two span models with classical dimensionless variables. For the Euler-Bernoulli beam model, the spatial problem in a dimensionless form is given by(147)d4X¯jXjdXj-Ωj4X¯jXj=0,where the dimensionless parameters areXj=xj/L,X¯j(Xj)=Xj/L,Ωj4=ρjAjω2/EjIj,Ωj4=μjΩ14.The Timoshenko model in a dimensionless matrix form is given by(148)ξj001Yj′′XΨj′′X+0-ξjξj1Yj′XΨj′X+Ωj200-ξj+ηjΩj2YjXΨjX=00,with the dimensionless parameters defined by(149)X=xL;YjX=YjxL;ΨjX=ψjx;Ωj2=ρAjω2l4EIj;ηj=IjAjL2;ξj=κGjAjL2EIj=κ2ηj1+νj.Then we have the dimensionless frequencies by each segment(150)Ωj2=ρAjω2L4EIjTimoshenko,Ωj4=ρAjω2L4EIjEuler-Bernoulli.These frequencies can be written as a relation ofΩ1. Ωj=(μj)1/2Ω1, for Timoshenko beam, and Ωj=μj1/4Ω1 for Euler-Bernoulli beam, being μj=(ρA)j/(ρA)1(EI)1/(EI)j.In Table6 we present results about dimensionless natural frequencies; we have considered variations in the position of the discontinuity point Rj=xj/L and variations of thickness of the second segment; constant width is considered, according to the following configurations.Table 6
The effects of the change of discontinuity pointR1=x1/L in the dimensionless natural frequencies for a cantilever double span beam with variations of thickness.
Configuration A Configuration B Case Model Dimensionless natural frequenciesΩ[i] Ω [ 1 ] Ω [ 2 ] Ω [ 3 ] Ω [ 4 ] Ω [ 1 ] Ω [ 2 ] Ω [ 3 ] Ω [ 4 ] R 1 = 0.2 TB-TB 1.4093 3.4619 5.6243 7.4545 1.7184 4.0974 6.5148 8.7633 TB-EB 1.4098 3.4721 5.6637 7.5438 1.7193 4.1198 6.5990 8.9637 EB-EB 1.4101 3.4769 5.6854 7.6159 1.7209 4.1381 6.6554 9.077  R 1 = 0.5 TB-TB 1.9972 3.5642 5.9015 8.2377 2.0150 4.0913 6.8545 9.1581 TB-EB 1.9985 3.5725 5.9372 8.3063 2.0162 4.1061 6.9121 9.3086 EB-EB 2.0004 3.5901 5.9831 8.5546 2.0198 4.1331 7.050 9.5862  R 2 = 0.8 TB-TB 2.0982 4.7587 6.6143 9.1855 1.9918 4.7325 7.3808 9.7778 TB-EB 2.1000 4.7673 6.6436 9.2371 1.9933 4.7395 7.4083 9.8619 EB-EB 2.1046 4.8267 6.7614 9.6670 1.9973 4.8113 7.6241 10.329Configuration A: double-span beam with variable thickness and constant width, t1=0.08, t2=0.03, w1=w2=0.03.Configuration B: double-span beam with variable thickness and constant width, t1=0.08, t2=0.05, w1=w2=0.03.Comparative results of dimensionless natural frequencies are presented in Table6.We observe from Table6 that the values of dimensionless natural frequencies of the coupled beam (TB-EB) are between the respective values obtained by TB-TB and EB-EB cases, wherein the EB-EB natural frequencies, as expected, are greater than the others. Also, we can observe that when the thickness of the second segment increases, the natural frequencies become greater, mainly the higher frequencies (third and fourth frequencies), except for the case R1=0.8 in which the first and second frequencies become greater. When R1 goes from 0.2 to 0.5, all frequencies increase, but when R1 goes from 0.5 to 0.8, this occurs only for Configuration A, and for Configuration B the first one frequency does not increase in this case.


## 9. Conclusions

This paper formulated in matrix form the time domain determination of forced responses and in the frequency domain the obtention of modes of elastic models that can be used in AFM when subjected to tip-sample interactions.Forced responses are determined by convolution of the input load with the time domain impulse matrix function. The corresponding matrix transfer function and modes of a multispan cantilever beam were determined in terms of solution basis which have the same shape and are generated by a fundamental solution.It was observed that the spatial amplitude diminishes with the inclusion of a piezoelectric patch when the pulse is positioned in the piezoelectric layer segment; this effect does not appear when the pulse is positioned in the third segment. The reduction of three segments for two segments, maintaining the others material and geometric properties, influences the magnitude of the natural frequencies. When the piezoelectric patch is removed, then the frequencies decrease.We worked with two Timoshenko beam models including surface effects through different approaches and changing the classical Timoshenko equations in different manners. For both models were presented the results which evidence that, in agreement with the literature [13, 14], for beam length on the order of nanometer to microns, the difference between natural frequencies is apparent and, by increasing the length of the microbeam, the results tend to Timoshenko classical theory; that is, the surface effects are significant only in nanoscale.Simulations were performed with double-span beams in which the first segment is modeled by the Timoshenko beam model and the second segment can be modeled by a Timoshenko (TB) or an Euler-Bernoulli (EB) model. Consistent results of TB-EB beams were obtained, in dimensional and dimensionless forms, when compared with EB-EB and TB-TB cases. The value of frequencies for this model (TB-EB) is located between the natural frequencies of uniform double-span models TB-TB and EB-EB.



---

*Source: 9268973-2018-10-18.xml*