Both surface elasticity and surface stress can result in changes of resonant frequencies of a micro/nanostructure. There are infinite combinations of surface elasticity and surface stress that can cause the same variation for one resonant frequency. However, as shown in this study, there is only one combination resulting in the same variations for two resonant frequencies, which thus provides an efficient and practical method of determining the effects of both surface elasticity and surface stress other than an atomistic simulation. The errors caused by the different models of surface stress and mode shape change due to axial loading are also discussed.
The application of the ansatz that nanostructure=bulk+surface  in continuum mechanics leads to the so-called core–shell model [2–11], in which the core is the bulk and the shell is the surface layer. In a surface layer, the total surface stress (τ) is given as follows [2,4,5,8–13]: 1.1where ϵ is the dimensionless strain and Cs is the surface modulus. Here, τ is the result of charge redistribution as the electrons respond to the effects of terminating a solid at a surface . By the thermodynamics definition, τ is a tensor associated with the reversible work to elastically stretch a pre-existing surface . We see that τ consists of two parts: σ and Csϵ; σ, which is strain independent, is often referred to as surface stress [16–18]; Csϵ, which is strain dependent, is often referred to as surface elasticity [17–19]. Surface elasticity is due to the formation of a surface layer that has a different elastic property from that of a bulk [3–5,20]. Surface relaxation  and sometimes surface reconstruction  are the two profound mechanisms responsible for the formation of a surface layer. Both σ and Cs have units of force per unit length (N m−1) compared with those of force per unit area (N m−2) for the bulk stress and modulus. Here, σ and Cs can be either positive or negative , which can thus either stiffen or soften a micro/nanostructure. In dynamics, the effect of the structure stiffness change (either stiffening or softening) is embodied in the shifts of the structure resonant frequencies.
Lagowski et al.  found that the resonant frequency of a microcantilever deviate significantly from those predicated by the axial load-free beam theory; they proposed that surface stress is the mechanism causing the resonant frequency shifts. Subsequently, Gurtin et al.  disputed Lagowski's explanation; they argued that the resonant frequency is independent of surface stress and that surface elasticity is the only mechanism responsible for the resonant frequency shifts. A rigorous three-dimensional analysis [17,18] shows that surface stress does change the structure stiffness. The effect of surface stress can be viewed as residual stress [9,10], which can change the structure bending stiffness , the morphology of the nanomaterial , the resonant frequencies , etc. Recently, many theoretical investigations [7,8,18,23] show that both surface elasticity and surface stress have an effect on the structure stiffness; however, the effect of surface elasticity is too small to explain the experimental observations. The effects of surface elasticity and surface stress are both size dependent, which stand out when the structure dimensions are small. A dimensional analysis is presented in this study to show how surface elasticity and surface stress, together with the structure dimensions, impact on the resonant frequency variations. Our dimensional analysis shows that, for a cylinder-like nanowire with length L and diameter D, the effect of surface stress is amplified by a factor of (L/D)2 compared with that of surface elasticity. When the magnitudes of σ and Cs are comparable with each other, the effect of surface stress on a slender structure stands out because (L/D)2 is a large number. The previous studies formulate the forward problem: Cs and σ are given or derived by atomistic simulation and then used in a model to predicate the nanowire buckling, bending and resonant frequency [2–11,24–31]. Resonant frequencies can be experimentally obtained with relatively high accuracies [13,16,32,33]. In practice, the following inverse problem is thus encountered: how to determine the effects of surface elasticity and of surface stress by measuring the resonant frequency shifts?
Lagowski et al.  and McFarland et al.  found surface stress by measuring the shifts of one resonant frequency and by assuming zero surface elasticity. Similarly, Gurtin et al.  found surface elasticity by proving that surface stress has no impact on the resonant frequency. Gavan et al.  were the first to use the two resonant frequencies of a microcantilever to determine the effects of surface elasticity and surface stress. Again, Gavan et al.  extracted the fitting values of surface stress/surface elasticity by assuming those of surface elasticity/surface stress to be zero, which excludes the general case that both surface elasticity and surface stress are non-zero . It is also noted that there is about 20% difference between the fitting values of surface elasticity for the first and second modes given by Gavan et al. . Surface stress can also be determined by the routine process of a static bending test via the Stoney formula or its modified forms [13,23,32,34–36]. However, the implicit assumption used in the static bending test is that the effect of surface elasticity is negligible. The static bending test alone cannot determine the effects of both surface elasticity and surface stress. Furthermore, static bending is induced by a differential surface stress [32,34–36]. For a micro/nanostructure with a symmetric surface layer [3–5,19], there is no differential surface stress and thus no bending, in which the Stoney formula cannot be applied . Similarly, other static tests such as nanoindentation and the three-point bending test can measure the variation of the nanowire mechanical properties owing to surface effects. However, these tests cannot differentiate the contributions to the variation by surface elasticity and surface stress. According to Song et al. , so far, there is no experimental method that can determine both non-zero surface elasticity and surface stress at the same time. This study presents an efficient method to uniquely determine the effects of surface elasticity and surface stress by using the fact that surface elasticity and surface stress have different impacts on different resonant frequencies. By solving this inverse problem, we actually provide a viable experimental method to determine the effects of both surface elasticity and surface stress. Because the surface stress effect is modelled as an axial load on the structure, two different models arise: the concentrated load model [3–8,13,16,24,35–37] and the distributed load model [35–37]. The difference between these two models is also discussed.
2. Model development
Figure 1a is a schematic diagram to illustrate the concentrated load model, in which the effect of surface stress is modelled as a concentrated load F applied at the beam free end. The concentrated load model leads to the following governing equation [3–8,16,24,35–37]: 2.1where m is the beam mass per unit length and m=ρA (ρ and A are the mass density and cross-sectional area of the beam, respectively); w=w(x, t) is the beam deflection. The concentrated load F is due to the strain-independent surface stress of σ, and F=πσD for a circular beam (D is the diameter of a circular beam). For a rectangular beam with a width b and a thickness h, F=2σb=2σA/h (A=bh) [13,16,19], which only includes the load induced by the surface stress on the upper and lower surfaces. We see that ∂2w/∂x2 is the curvature of the upper/lower surface and F is in the horizontal direction; −F∂2w/∂x2 is a distributed transverse load given by the Young–Laplace (YL) formula . Therefore, the concentrated load model of equation (2.1) is also referred to as the YL model [28,31]. Here, positive F is tensile and negative is compressive. J* is the beam bending rigidity with the presence of a surface layer, which is calculated as the following: 2.2where E* is the beam's effective Young's modulus. The in-plane surface stress in the upper and lower surfaces is usually in two directions. For example, the upper and lower surfaces experience the surface stress in both x and y directions; E*=E/(1−ν) is to account for this biaxial loading scenario [18,35]. Here, E and ν are the beam's Young modulus and Poisson ratio, respectively. If a beam is (assumed) to bend in a cylindrical shape , E*=E/(1−ν2) [14,16,23]; or simply, E*=E [4,5,13,19]. We have I=bh3/12 for a rectangular beam and I=πD4/64 for a circular beam. We denote t as the surface layer thickness, as seen in figure 1c. When t≪b and t≪h for a rectangular beam, or t≪D for a circular beam, equation (2.2) recovers the following results given by He & Lilley [4,5]: 2.3For the concentrated load model, or say, the YL model to apply, the following three conditions need to be satisfied: (i) , (ii) (for a rectangular beam) b≫h, and (iii) the surface stress in the transverse direction is neglected . Because t is very small, the first condition can only be violated when the nanostucture is so small, and the surface stress effect is huge enough to cause a deflection with very large curvature. The second condition, in other words, is to say that the in-plane distributed load induced by the side surfaces can be ignored. Song et al.  developed a model incorporating the effects of side surfaces to overcome the second condition. Olsson & Park  developed a three-dimensional model incorporating the effects of the transverse surface stress to overcome the third condition. The modified YL models developed by Song et al.  and Olsson & Park  in essence change the expressions of J* and F in equation (2.1). As shown later, our method is to detect the variations of J* and F by measuring the shifts of the beam resonant frequencies. Therefore, the method can be applied to the concentrated load model and modified YL models alike.
The boundary conditions of a cantilever beam under a concentrated load at its free end are the following [24,35,36,40]: 2.4It is worth pointing out that the concentrated load F appears in the fourth boundary condition, which determines the resonant frequencies and mode shapes of a beam [24,36,40]. We see that ∂3w/∂x3(L, t)=0 has been (incorrectly) taken as a boundary condition for this concentrated load model of a nanowire [4,5], which should be responsible for the abnormal resonant frequency deviation behaviour of a cantilever as presented in fig. 2 of He & Lilley .
On the other hand, Finot et al.  argued that surface stress should be viewed as the sum of two contributions: one is an axial force per unit length and the other is a moment per unit cross section. Using the Gurtin–Murdoch surface elasticity theory, Wang et al.  also showed that surface stress effectively exerts a distributed load. The distributed load model shown in figure 1b gives the following governing equation [35,36]: 2.5where s=F/L is a uniformly distributed load. The boundary conditions are as follows: 2.6According to Lachut & Sader , as ‘the cantilever free end is unrestrained with zero net force’, the above concentrated load model as shown in figure 1a, equations (2.1) and (2.4) are ‘therefore in direct violation of Newton's third law’. In comparison, the distributed load model indeed guarantees the zero net force of the cantilever free end, which has also been shown to better model the beam static deflection under differential surface stress loading . For the concentrated load model, Gurtin et al.  argued that owing to membrane stretching, the two surface layers (of a very thin beam, i.e. the h≪b case) store an energy of , which, in terms of the force balance, exerts a distributed transverse load of 2bσ∂2w/∂x2. The concentrated load owing to surface stress is F=2bσ, which also generates a membrane-stretching energy of and a distributed transverse load of −F∂2w/∂x2. We have −F∂2w/∂x2+2bσ∂2w/∂x2=0, and Gurtin et al.  thus concluded that surface stress has no influence on the beam resonant frequency. Here, the energy contribution of the surface stress in Gurtin's model has been accounted for twice: as and as . Gurtin's membrane-stretching energy of treats the surface stress effect as a concentrated load. In comparison, for the distributed load model, the membrane-stretching energy is calculated as .
By introducing ξ=x/L, and W=w/L, the governing equation of the concentrated load model, equation (2.1), is non-dimensionalized as follows: 2.7where the dimensionless quantities Δ and N are defined as follows: 2.8Here, Δ indicates the bending rigidity ratio of the surface layer to the bulk; N is the ratio of the concentrated load to the beam (transverse) stiffness. The size-dependent properties of Δ and N are also noted. Compared with Δ, N has an amplification factor of (L/h)2 for a rectangular beam and (L/D)2 for a circular beam; L/h or L/D have been identified as important geometric parameters in the buckling and resonant frequencies of a nanowire [25,26,30]. Surface elasticity (Δ) in essence changes the effective Young modulus of a micro/nanostructure . For a nanocomposite with a spherical shell structure, the change of the effective Young modulus can be significant . However, because of the amplification factor of (L/h)2, the effect of surface elasticity can be remarkably reduced when compared with that of surface stress for a slender beam structure, as discussed later.
The boundary conditions of equation (2.4) are now non-dimensionalized as follows: 2.9By assuming W(ξ, τ)=V (ξ)eiωτ (ω is the dimensionless circular frequency), V (ξ) can be solved from equation (2.7) as the following: 2.10where ϕT, ϕ0 and ϕC are the mode shapes when the concentrated axial load is tensile, zero and compressive, respectively. Here, C1, C2, C3 and C4 are four unknown constant; f0, f1 and f2 are the quantities defined as follows: 2.11By substituting equation (2.10) into equation (2.9), an eigenvalue problem is formed and the eigenfrequencies (resonant frequencies), ω, can thus be computed. Clearly, as seen in equation (2.10), N/(1+Δ) has a direct impact on the beam mode shapes and eigenfrequencies.
For the distributed load model, the governing equation of equation (2.5) and boundary conditions of equation (2.6) are non-dimensionalized as follows: 2.12 2.13To compute the eigenfrequencies of equation (2.12), the Galerkin method [24,43] is applied, which assumes the following form for W(ξ, τ): 2.14where H is the mode number and is the jth mode shape of a uniform cantilever beam with the zero axial concentrated load (N=0) as given in equation (2.10) . Clearly, the mode shape of also satisfies the boundary conditions of equation (2.13). Substituting equation (2.14) into equation (2.12), time and integrating from 0 to 1, the following eigenvalue problem is formed: 2.15K and M are the H×H matrices of stiffness and mass, which are given as follows: 2.16There is no damping in equations (2.7) and (2.12), which is also the case in many studies [3,4,6,7,13,18,19,33]. The presence of damping reduces the resonance, and damping can be determined by the so-called half-power method from the frequency response curve obtained by experiment . However, the high quality factor (i.e. small damping) is a much sought-after property in many applications of micro/nanoresonators, which can significantly enhance the sensitivity . Therefore, the shifts of resonant frequencies owing to small damping can be ignored in many dynamic models of microstructure vibration in air or vacuum. However, when a microstructure vibrates in liquid  or in contact with a viscous material , damping plays an important role and cannot be ignored. The eigenfrequency computation of a damped system using the above Galerkin method is presented by Zhang & Murphy . It also needs to be kept in mind that damping, or say, the quality factor, is different for different modes . It is also worth mentioning that equations (2.7) and (2.12) are confined to the linear theory of surface elasticity. The finite deformation effects can be important in some cases . The higher-order terms of strain owing to finite deformation will result in a nonlinear governing equation such as the Duffing equation .
3. Results and discussions
Figure 2 examines how the first mode shape given in equation (2.10) changes as the axial concentrated load N varies when Δ=0. Three values of N=−2.46, N=0 and N=10 are taken. Timoshenko's  mode shape of 4(ξ2/2−ξ3/3+ξ4/12) is also presented for comparison (the factor of 4 is to normalize the shape). Here, N=−π2/4≈−2.46 is the dimensionless concentrated buckling load of a cantilever beam . If N is smaller than this critical value, the beam enters the post-buckling regime and equation (2.7) cannot apply. The mode shape of N=0 is very close to Timoshenko's, which is actually the static deflection shape of a cantilever under a uniformly distributed load . Clearly, as seen in figure 2, the change of mode shape is significant as N varies. This shape variation can cause significant errors in the eigenfrequency computation of the cantilever. Therefore, when the Rayleigh–Ritz or Galerkin method is used for eigenfrequency computation, the mode shape should be recalculated whenever N varies [24,36,40]. Figure 3 shows the second and third mode shape changes with variation of N. Clearly, the changes of higher-mode shapes with variation of N are much less than those of the first mode shape.
Figure 4 plots the first eigenfrequency (ω1) calculated by the concentrated load model, the distributed load model and Timoshenko's load model. Timoshenko  gave the following approximate first eigenfrequency expression of the concentrated load model: 3.1where is the first eigenfrequency of a cantilever with zero Δ and N. The first three eigenfrequencies with zero Δ and N are given by Chang & Craig  as 3.2Because of the softening effect of compressive load, the effective stiffness of a structure becomes zero at the buckling load , which causes the first eigenfrequency to be zero. As mentioned above, Gurtin et al.  argued that surface stress has no impact on the beam resonant frequencies. However, Gurtin & Murdoch  concluded that the compressive surface stress can cause buckling. It is noted that, at N=−2.46 of the buckling load, Timoshenko's ω1 is not zero. Timoshenko's eigenfrequency is close to our concentrated load model when the magnitude of N is relatively small, and the difference enlarges as the magnitude of N increases. Because the Rayleigh–Ritz method is used in the derivation by Timoshenko , the reason for the enlarging error between equation (3.1) and our concentrated load model is that Timoshenko's mode shape deviates more and more from ours as the magnitude of N increases, as shown in figure 2. At N=−2.46, ω1 of the distributed load model is not zero, either. The reason has already been given by Timoshenko  that for a cantilever under a distributed load of s=F/L, its effect on ω1 is as if a concentrated load of 7F/20 is applied at the free end. Compared with the concentrated load, the distributed load has a much less softening effect when N is negative and a much less stiffening effect when N is positive. Therefore, as shown in figure 4, ω1 of the distributed load model is smaller/larger than that of the concentrated load when N is tensile/compressive. When surface stress (N) is tensile, Olsson & Park  also found that both the concentrated load (YL) model and the modified YL model of Song et al.  have significantly higher ω1 values than those obtained by an atomistic simulation. The difference between the concentrated load model, distributed load model and Timoshenko's model of equation (3.1) vanishes when N=0. We see that ω1 values predicated by the three models are the same at N=0, and (N, ω1)=(0,3.516) is the intersection point of the three curves, which is marked as a circle in figure 4. Figure 5 plots the second and third eigenfrequencies as N varies. Again, both ω2 and ω3 of the distributed load model are larger/smaller than those of the concentrated load model when N is compressive/tensile; the intersection points marked as circles are (N, ω2)=(0,22.034) for the second eigenfrequency curve and (N, ω3)=(0,61.697) for the third eigenfrequency curve.
Now let us discuss how to use the shifts of eigenfrequencies to determine the surface elasticity of Δ and the surface stress of N. As shown above, the concentrated load model underestimates/overestimates the eigenfrequencies when N is compressive/tensile. Here, only the distributed model is used for eigenfrequency computation. The following quantities are taken from Gavan et al.'s  experimental data of a SiNx rectangular microcantilever: L=100 μm, b=8 μm, h=100 nm, Cs≈1170 Nm−1, σ≈1 Nm−1 and E=300 GPa. We see that Δ and N are thus calculated as Δ=0.235 and N=0.81 using equation (2.8). Substituting these two Δ and N values into equation (2.15), the first three eigenfrequencies are obtained as follows: 3.3Clearly because of the positive Δ and N, the beam stiffness is enhanced and the eigenfrequencies thus all increase when compared with those in equation (3.2). In the beam resonance test, Δ and N are unknown; the eigenfrequencies are extracted from the beam frequency response curves [13,16,32,33]. Therefore, using the eigenfrequencies to determine Δ and N forms an inverse problem. A similar inverse problem is also encountered in the case of the resonant frequency shifts induced by surface stress and mass .
Figure 6 plots the variation of the first eigenfrequency, ω1, as a function of Δ and N, which is a tilted plane. ω1 increases monotonically with the increase of both Δ and N. The level plane is the one with the fixed first eigenfrequency value of ω1=4.065. The intersection of these two planes is all the combinations of Δ and N, which results in the same first eigenfrequency of ω1=4.065. The intersection is a line, which also indicates that these combinations are infinite.
Figure 7 plots the variation of the second eigenfrequency, ω2, as a function of Δ and N. Compared with figure 6, the plane tilts much more along the Δ-axis and much less along the N-axis. Physically, this means that Δ and N have different impacts on different eigenfrequencies, which is also the mechanism that leads us to use the eigenfrequency shifts to determine the effects of surface elasticity (Δ) and surface stress (N). From the viewpoint of the forward problem, Park & Klein  found that the given values of surface elasticity and surface stress have different impacts on the first and second bending resonant frequencies; Dorignac et al.  also noted that the variation of the first eigenfrequency is the most sensitive to N and the variations of the higher eigenfrequencies are much less sensitive to N. In figure 7, the level plane is the one with the fixed value of ω2=24.628. The intersection line of the two planes is the combinations of Δ and N, which result in ω2=24.628. For any given Δ and N, each eigenfrequency is uniquely determined by equation (2.15). As an inverse problem, for a given eigenfrequency, there are infinite combinations of Δ and N. However, for two given eigenfrequencies, the combinations of Δ and N are two lines intersecting each other, which is used here to uniquely determine the combination of Δ and N. When the two lines of figures 6 and 7, which are obtained by intersecting two planes, are projected on to the Δ−N plane, the two lines intersect, as shown in figure 8. The intersection point is marked as a circle, which happens to be exactly (Δ, N)=(0.235,0.81).
Atomistic simulation is often used to obtain Cs and σ, which are determined by the underlying lattice structure and interatomic potential [1,26,52]. For example, CAgs=1.22 Nm−1 and σAg=0.89 Nm−1 [1,5], with EAg=76 GPa for Ag(001); CAus=−3.6 Nm−1 and σAu=1.4 Nm−1 [1,4], with EAu=37 GPa for Au(100); CSis=−11.5 Nm−1 and σSi=−0.5 Nm−1 [2,52], with ESi=130 GPa for Si(100). For a nanowire with dimensions of L=0.5 μm and D=40 nm, the corresponding Δ and N can be calculated using equation (2.8) as ΔAg=3.2×10−3 and NAg=2.9276, ΔAu=−1.95×10−2 and NAu=9.4595, ΔSi=−1.77×10−2 and NSi=−0.9615. For SiNx calculated above, ΔSiNx/NSiNx≈0.29, and as seen in figure 6, Δ and N contributions to the first eigenfrequency variation are comparable. For the nanowires of the above three materials, ΔAg/NAg≈1.1×10−3, ΔAu/NAu≈−2.1×10−3 and ΔSi/NSi≈1.84×10−2. As seen in equation (2.8), Δ and N are both size-dependent quantities; at the same time, L=0.5 μm and D=40 nm are relatively small dimensions. Therefore, for those materials that have very small Δ/N ratios such as silver, gold and silicon, the effect of surface stress is the dominant factor influencing the micro/nanobeam eigenfrequencies [4,16]. As given in equation (2.8), Δ∝CsD−1 and N∝σD−1(L/D)2. For a beam structure, L/D is around 10 or larger. Because of this additional amplification factor of (L/D)2, the magnitude of Cs/σ needs to be around 102 or larger for the surface elasticity effect to be comparable with that of surface stress. For example, Cs/σ≈103 for SiNx, as computed above. Shenoy  computed Cs and σ for several face centred cubic metals, whose |Cs/σ| ratios range around 1–102. By checking the two ratios of Cs/σ and L/D, the effects of surface elasticity and surface stress can be estimated. Similarly, Chiu & Chen  derived an analytical expression for the buckling load of a nanowire, which indicates that the surface stress contribution to the buckling load is also amplified by a factor of (L/D)2 compared with that of surface elasticity. It is also necessary to bear in mind that in the modified YL models [28,31], the transverse surface stress and surface stress from side surfaces can significantly enhance the effective surface elasticity effect by increasing the nanowire bending stiffness.
It is noted that Cs and σ calculated by atomistic simulation can be significantly different for the same material. For example, the surface stress of Au(100) is calculated as σAu=1.4 Nm−1  and σAu=4.56 Nm−1 using the first-principle approach and σAu=1.79 Nm−1 using the semi-empirical approach . For a nanowire of L=0.5 μm and D=40 nm, the corresponding N values are NAu=9.4595, NAu=30.8109 and NAu=12.0946, which are the values large enough to have a significant impact on the beam eigenfrequencies. Using the above method, it is not difficult to tell which surface stress value is the correct one.
Two models of the resonant frequencies of a micro/nanocantilever beam with the presence of surface elasticity and surface stress are presented and compared. The concentrated load model violates the boundary conditions at the cantilever free end, which also underestimates/overestimates the resonant frequencies when the load is compressive/tensile. The distributed load model reserves the net zero force boundary condition at the free end. Surface elasticity and surface stress play different roles in the variation of the beam resonant frequencies. Mathematically, the effect of surface elasticity is embodied in the fourth-order differential term of the governing equation; the effect of surface stress is in the second-order differential term (both the concentrated load and distributed load models) and the first-order differential term (the distributed load model only). That surface elasticity and surface stress have different impacts on different resonant frequencies is used as a mechanism to uniquely determine their effects. The method of using the shifts of different resonant frequencies to determine the effects of surface elasticity and surface stress is only applied to the distributed load model and a cantilever beam in this study. However, the method can also be applied to the concentrated load model, modified YL models and a beam with different boundary conditions. The method is shown to be accurate and provides a new experimental approach for determining surface elasticity and surface stress.
Y.Z. acknowledges support from the National Natural Science Foundation of China (nos. 110212622, 11023001 and 11372321) and Chinese Academy of Sciences (no. KJCX2-EW-L03).
- Received July 10, 2013.
- Accepted August 23, 2013.
- © 2013 The Author(s) Published by the Royal Society. All rights reserved.