## Abstract

This paper presents a method to simultaneously predict the elastic modulus, axial load, and boundary conditions of a nanoelectromechanical system (NEMS) beam from a minimum of two measured natural frequencies. The proposed method addresses the challenges of the inverse problem at the nano scale, which include high natural frequencies, small geometric beam dimensions, and measurements limited to natural frequencies. The method utilizes a finite element model of an Euler–Bernoulli beam under axial loading to predict the response of the beam with axial loading and flexible boundary conditions. By expressing the finite element model in terms of dimensionless beam parameters, the proposed method may be applied to nano scale beams while maintaining numerical stability of the finite element equation of motion. With the stabilized finite element model, the NEMS beam properties are predicted by iterating through values of dimensionless beam parameters until the normalized error between predicted and measured natural frequencies is minimized. A key feature of the proposed method is the simultaneous prediction of the elastic modulus during the iterative search, resulting in a reduction of the search space and significant computational savings. Additionally, the proposed method readily accommodates an arbitrary number of measured natural frequencies without the reformulation of procedures and analyses. Numerical examples are presented to illustrate the proposed method’s ability to predict the elastic modulus, axial load, and boundary conditions. The proposed method is applied to experimental measurements of a NEMS beam, where the normalized error between predicted and measured natural frequencies is reduced below 10−3.

## 1 Introduction

Nanoelectromechanical system (NEMS) devices have recently gained significant attention due to their ability to actuate and sense on the nanoscale [1]. Common NEMS devices include nanomechanical resonators such as nanotubes, cantilevers, doubly clamped beams, membranes, and more. With relevant dimensions in the deep submicron, NEMS devices offer access to an unprecedented set of dynamic properties including very high resonance frequencies in the megahertz to gigahertz range, mechanical quality factors in the thousands, and active masses in the picograms (10−12 g). Given this extraordinary set of parameters, a NEMS resonator can respond to even the slightest changes in the environment with extreme sensitivity [2]. As a result, NEMS devices are currently being developed as ultrasensitive sensors of physical quantities, such as mass [39], charge [10], and viscosity [11]. Complimentary to this effort, NEMS devices that are optimized for operation in water, referred to as bio-NEMS, are being envisioned for sensing biological quantities and phenomena [1214]. Such unprecedented sensing capabilities enable NEMS devices to serve as lab-on-a-chip biosensors, with applications that include mass spectrometry, label-free detection of biological molecules, and more [1].

To accurately sense the aforementioned perturbations, properties and parameters of the NEMS resonator must be precisely known. These properties include the elastic modulus, prestress, and boundary conditions of the resonator. Accurate properties of the NEMS resonator are paramount for an accurate nominal model. Such models are heavily relied upon when predicting the changes in response due to physical perturbations of a resonator. For example, when performing mass sensing, an accurate model of the nominal NEMS resonator is needed to predict the mass and position of a molecule adsorbed to the physically perturbed resonator [7]. In this example, the measured change in response is attributed to a change in the nominal model, which then characterizes the added mass. Consequently, to sense and measure perturbations of a NEMS resonator, it is critical to precisely know the properties of the resonator for an accurate nominal model.

It is typically very challenging, if not impossible, to predict the exact parameters of a NEMS resonator. There are several reasons for this. First, variations in the material properties of the semiconducting wafer, unnoticeable at the macroscopic scale, will lead to large changes in the parameters of a single NEMS device. For instance, due to randomness of ion implantation, wafers will have spatial stress variations at the length scales of NEMS [15,16]. Likewise, spatial fluctuations in the processing parameters during nanofabrication will result in significantly different parameters for two NEMS devices on the very same chip. Even after fabrication, the device properties will be highly susceptible to environmental factors such as temperature, given the fact that small perturbations lead to big changes in the parameters of NEMS. Moreover, boundary conditions of the NEMS device are typically hard to model: undercuts due to fabrication, mismatch between epilayers, and deposited films, which all complicate the boundary conditions. In summary, more modeling work is needed in the NEMS domain, and this paper addresses this need by presenting a new method for predicting initially uncertain properties of NEMS resonators from measured natural frequencies. Specifically, the present work focuses on beams with an axial load and flexible boundary conditions modeled by translational and rotational stiffnesses.

Axial loading is often present in NEMS beams and the thin films they are fabricated from due to the difference in thermal expansion between the substrate and the deposited material [17,18]. Axial loading has also been observed in macroscale beams, with examples such as cables on a bridge, diagonal braces of a truss, and struts and ties of a space truss structure. Methods have been proposed that predict the axial load of beams from measured displacements using analyses based on the solution to the equation of motion of axially loaded beams [1921]. However, these methods require displacement measurements from at least five sensors, which is infeasible for NEMS beams where measurements are often limited to natural frequencies. Alternatively, sensitivity-based methods have been proposed that predict the dependence of the beam’s modal properties on axial loading [2225]. With the computed sensitivities, the methods iteratively predict the uncertain beam properties until the error between the measured and the predicted natural frequencies is minimized. However, these methods assume the elastic modulus or boundary conditions to be known.

Stachiv et al. [26] presented a more general method, which predicted the elastic modulus, axial loading, density, and thickness of doubly clamped NEMS beams. By analyzing the dependence of the in-plane and out-of-plane modes on the axial load, Stachiv et al. developed a method to simultaneously predict the initially uncertain beam properties. Their method, however, required measurements from beams with different applied prestress forces and in-plane and out-of-plane fundamental frequencies. Recently, Pratab and Behera [27] presented a method to simultaneously predict the elastic modulus and axial load of buckled microelectromechanical system beams from measured natural frequencies. Although their method only requires measurements of a single even and odd bending mode, it is specific to buckled beams as the method utilizes the invariance of even modes to compressive loading.

To overcome the limitations of previous works, a new method is proposed that simultaneously predicts the elastic modulus, axial load, and boundary conditions of a beam from measured natural frequencies. The proposed method has the significant advantage of requiring only two measured natural frequencies. In addition, the method is general as it may be applied to beams that have either a tensile or a compressive load and any boundary condition. The proposed method utilizes a finite element model (FEM) of an Euler–Bernoulli beam to formulate a generalized eigenvalue problem that relates the axial load and boundary conditions to the beam’s natural frequencies. After algebraic manipulation of the equations of motion, the generalized matrices are expressed in terms of dimensionless beam parameters that are related to the axial load and boundary stiffnesses, while the eigenvalues are related to a dimensionless frequency parameter. The method iterates through values of the dimensionless axial load and stiffness parameters, while simultaneously predicting an elastic modulus by comparing the evaluated dimensionless frequency to measured natural frequencies. The elastic modulus, axial load, and boundary conditions are then predicted by iterating over the dimensionless parameters until the normalized error between measured and predicted natural frequencies is minimized. The proposed method is illustrated with numerical simulations, where the uncertain beam properties are accurately predicted from simulated measurements. Properties of a NEMS beam are then predicted from experimentally measured natural frequencies.

The remainder of this paper is organized as follows. In the following section, the FEM of an Euler–Bernoulli beam is presented. With the developed FEM, an analysis that evaluates the modal properties of the beam is presented in Sec. 3. The analysis introduces a dimensionless frequency parameter and dimensionless beam parameters that are related to the axial load and boundary stiffnesses. In Sec. 3.1, the dimensionless frequency’s dependence on the dimensionless beam parameters is examined to investigate the effect axial loading and boundary conditions have on natural frequencies and to illustrate potential limitations of the proposed method. In Sec. 4, the proposed method is presented, followed by an assortment of numerical examples in Sec. 5 intended to illustrate the accuracy of the method. The proposed method is then applied to experimental measurements from a NEMS resonator in Sec. 6, where ideal and flexible boundary conditions are considered.

## 2 Finite Element Model

To begin, we first consider the continuous model with the equation of motion of a beam with axial loading
$EI∂4w∂x4+ρA∂2w∂t2−P∂2w∂x2=f(x,t)$
(1)
where x is the location on the beam, t is time, w is the lateral beam displacement about the y-axis such that the beam’s thickness is along the direction of vibration, E is the elastic modulus, I is the moment of inertia of the beam cross section about the y-axis, $ρ$ is the mass density, A is the cross-sectional area, P is the axial force, and f(x, t) is the external force per unit length. The axial force is tensile when P is positive and compressive when P is negative. In the present work, the beam is assumed to be undamped. Although it is out of the scope of this paper, damping may be readily added by including a dissipative force in the equation of motion.
To predict the solution of Eq. (1), FEM of an Euler–Bernoulli beam with axial loading is considered. The FEM in consideration is depicted in Fig. 1, where the beam of length L is discretized into Q beam elements, such that each beam element has a length
$l=LQ$
(2)
The equation of motion for a beam element in the frequency domain is as follows:
$(k+kax+(iω)2m)x=f$
(3)
where x is the elemental displacement vector, f is the elemental force vector, k is the elemental bending stiffness matrix, kax is the elemental axial loading stiffness matrix, and m is the elemental mass matrix. The displacement vector is as follows:
$x={wnθnlwn+1θn+1l}$
(4)
where wn and $θn$ are the lateral displacement and slope, respectively, of the nth node. The elemental force vector is as follows:
$f={FnMn/lFn+1Mn+1/l}$
(5)
where Fn and Mn are the applied force and moment, respectively, of the nth node.
Fig. 1
Fig. 1
Close modal
A cubic function in x is used to model the lateral displacement of the beam element as it is an appropriate choice for a four degree-of-freedom (DOF) element, such as the beam element described in Eqs. (3)(5). Furthermore, the cubic displacement function satisfies the beam differential equation and the conditions of displacement and slope continuity at nodes shared by two elements [28]. As a result, the elemental beam matrices are found to be [29,30]
$k=EIl3[126−12664−62−12−612−662−64]$
(6a)
$kax=P30l[363−36334−3−1−36−336−33−1−34]$
(6b)
$m=ρAl420[1562254−1322413−35413156−22−13−3−224]$
(6c)
When deriving the elemental beam matrices in Eq. (6), the rotational DOFs and applied moments were intentionally expressed in terms of the elemental length l in Eqs. (4) and (5), respectively. This was done so that each elemental beam matrix in Eq. (6) is expressed as a product of a scalar and a matrix, where all dependencies on material properties, dimensions, and axial loading reside in the scalar. This expression will be useful for the analysis in Sec. 3.
With the elemental beam matrices in Eq. (6), the global finite element matrices are assembled to model the transverse vibration of the entire beam [30]. The resulting equation of motion is as follows:
$(K+Kax+(iω)2M)X=F$
(7)
where K is the global bending stiffness matrix, Kax is the global axial bending stiffness matrix, and M is the global mass matrix. Similarly, X is the global displacement vector of the beam and F is the global force vector. In Sec. 2.1, boundary conditions of the beam are discussed.

### 2.1 Boundary Conditions.

As depicted in Fig. 1, boundary conditions on the translational and rotational degrees-of-freedom are modeled in the FEM at the left and right ends of the beam to simulate the physical boundary conditions of the structure. Commonly, boundary conditions are modeled by either fixing degrees-of-freedom or allowing degrees-of-freedom to be free such that no force or moment is applied. In this paper, flexible boundary conditions are considered to more accurately model the physical boundary conditions of the NEMS beam. With flexible boundary conditions, when the ends of the beam translate and rotate, a force and moment are applied, respectively. Flexible boundary conditions are modeled with translational and rotational springs at the ends of the beam, as depicted in Fig. 2. For the present work, the boundary conditions are assumed to be symmetric such that the translational and rotational stiffnesses at the left end are equal to those at the right end.

Fig. 2
Fig. 2
Close modal
With flexible boundary conditions, the beam is modeled with the equation of motion
$(K+Kax+Kt+Kr+(iω)2M)X=F$
(8)
where Kt and Kr are the global translational and rotational boundary condition stiffness matrices, respectively. With the displacement and force vector defined in Eqs. (4) and (5), the global boundary condition stiffness matrices are as follows:
$Kt=kt[10⋯0000⋯00⋮⋮⋱⋮⋮00⋯1000⋯00]$
(9a)
$Kr=krl2[00⋯0001⋯00⋮⋮⋱⋮⋮00⋯0000⋯01]$
(9b)
where kt is the translation spring constant and kr is the rotational spring constant. It should be noted that nonsymmetric boundary conditions may be modeled by letting the springs at the left and right ends of the beam be different. This would result in four global boundary conditions stiffness matrices, two for each end to model the translational and rotational stiffnesses, where each matrix would have only one nonzero element.

## 3 Evaluating Modal Properties

In this section, the beam FEM introduced in Sec. 2 is used to evaluate modal properties of the beam. Specifically, the dimensionless frequencies and mode shapes of the beam will be found by evaluating a generalized eigenvalue problem. The analysis provides a straightforward alternative to numerically solving the characteristic equation when computing the natural frequencies of the beam, which will be essential for the proposed method.

To begin, free vibration is considered by setting F = 0, which simplifies Eq. (8) to
$(K+Kax+Kt+Kr−ω2M)X=0$
(10)
Recalling Eqs. (6) and (9) and assuming that the beam is homogeneous with constant cross-sectional area and that all elements are the same size with equal length l, the equation of motion in Eq. (10) may be rewritten as follows:
$(EIl3K~+PlK~ax+ktK~t+krl2K~r−ω2(ρAl)M~)X=0$
(11)
where the global FEM matrices are expressed explicitly in terms of their equivalent scalar-matrix product. Here, the matrices $K~$, $K~ax$, $K~t$, $K~r$, and $M~$ are the global FEM beam matrices with all material properties, dimensions, and axially loading factored out. These matrices may be considered to be dimensionless FEM matrices that model the dynamics of a beam but are independent of properties specific to the structure under consideration. Dividing both sides of Eq. (11) by (EI/l3) results in
$(K~+Pl2EIK~ax+ktl3EIK~t+krlEIK~r−ω2ρAl4EIM~)X=0$
(12)
which may be further simplified by substituting in Eq. (2) and introducing dimensionless parameters, such that
$(K~+αQ2K~ax+βQ3K~t+χQK~r−Ω2Q4M~)X=0$
(13)
where α, $β$, $χ$, and Ω are dimensionless beam parameters defined as follows:
$α=PL2EI$
(14a)
$β=ktL3EI$
(14b)
$χ=krLEI$
(14c)
$Ω=ωρAL4EI$
(14d)
Here, α is considered to be the dimensionless axial load parameter that determines the effect that the axial load has on the response of the beam. Similarly, $β$ and $χ$ are the dimensionless translational and rotational stiffness parameters, respectively. The dimensionless parameter Ω is considered to be the dimensionless frequency parameter [31,32], as it is related to the natural frequencies of the beam in Eq. (14d). The formulation in Eq. (13) provides a novel method for modeling an FEM Euler–Bernoulli beam with dimensionless beam parameters and FEM matrices. This formulation is significant, as it offers numerical stability when modeling beams on the nanoscale. In addition, the formulation in Eq. (13) will allow for uncertain beam properties to be readily predicted through the dimensionless beam properties and the evaluation of the dimensionless frequency parameter in the proposed method.
To solve for the dimensionless frequency parameters and natural frequencies of the beam, Eq. (13) is rearranged in the form of a generalized eigenvalue problem:
$Aϕn=λnBϕn$
(15)
where A and B are the generalized matrices that are related to the factored FEM matrices by
$A=K~+αQ2K~ax+βQ3K~t+χQK~r$
(16a)
$B=M~$
(16b)
and $λn$ and $ϕn$ are the eigenpairs of the system. Here, $ϕn$ are the dimensionless mode shapes of the beam and $λn$ are related to the dimensionless frequencies by
$λn=Ωn2Q4$
(17)
When the generalized eigenvalue problem in Eq. (15) is evaluated, the dimensionless frequencies are found by
$Ωn=λn1/2Q2$
(18)
which are related to the dimensional natural frequencies of the beam by
$ωn=ΩnEIρAL4$
(19)
With the dimensionless FEM in Eq. (16), the generalized eigenvalue problem in Eq. (15) may be readily evaluated for any length scale while maintaining numerical stability. This model will be essential for the proposed method when predicting properties of NEMS beams.

From Eqs. (15)(17), it is clear that the dimensionless frequency parameters, and hence natural frequencies, depend on the axial load and boundary conditions of the beam that are modeled with the dimensionless parameters α, $β$, and $χ$. In Sec. 3.1, the effect of axial loading and boundary conditions on the dimensionless frequencies is examined to determine reasonable ranges for such parameters with respect to the proposed method.

### 3.1 Dimensionless Frequency.

Consider the case where the beam is unloaded, such that α = 0. The dimensionless frequencies for the first four modes, excluding the rigid-body mode, of beams with ideal boundary conditions, which consist of fixed and free degrees-of-freedom, are listed in Table 1. For this example and the purpose of examining trends, the boundary conditions were modeled by setting $β=χ=0$ in Eq. (16a), thus reducing Eqs. (15) and (16) to the generalized eigenvalue problem of an ideal free-free beam. The various ideal boundary conditions were then modeled by fixing the appropriate degrees-of-freedom with the method presented by Wu et al. [33]. With the imposed boundary conditions, the generalized eigenvalue problem was evaluated and the dimensionless frequencies were found with Eq. (18). The results in Table 1 are in agreement with expected values that are often found by numerically solving the characteristic equation [30]. In Sec. 3.1.1, the effect of including an axial tensile load is examined.

Table 1

Dimensionless frequency of first four modes, excluding the rigid-body mode, for beams with different ideal boundary conditions and axial loading

Boundary
conditions
α = 0α = 1α = 10α = 50α = 100α = 500α = 1000α = 10,000
Ω1 = 9.8696
Ω2 = 39.4784
Ω3 = 88.8264
Ω4 = 157.9137
Ω1 = 10.3575
Ω2 = 39.9753
Ω3 = 89.3250
Ω4 = 158.4129
Ω1 = 14.0038
Ω2 = 44.1965
Ω3 = 93.6931
Ω4 = 162.8369
Ω1 = 24.3082
Ω2 = 59.4346
Ω3 = 111.0471
Ω4 = 181.1972
Ω1 = 32.9298
Ω2 = 74.2050
Ω3 = 129.5098
Ω4 = 201.8120
Ω1 = 70.9381
Ω2 = 145.9375
Ω3 = 228.6993
Ω4 = 322.3252
Ω1 = 99.8349
Ω2 = 202.5758
Ω3 = 310.9929
Ω4 = 427.6101
Ω1 = 314.3143
Ω2 = 629.5576
Ω3 = 946.6544
Ω4 = 1266.520
Ω1 = 22.3733
Ω2 = 61.6728
Ω3 = 120.9034
Ω4 = 199.8595
Ω1 = 23.4507
Ω2 = 62.5490
Ω3 = 121.6736
Ω4 = 200.5703
Ω1 = 31.3892
Ω2 = 69.8892
Ω3 = 128.3823
Ω4 = 206.8526
Ω1 = 52.8630
Ω2 = 95.1260
Ω3 = 154.3553
Ω4 = 232.5673
Ω1 = 70.1390
Ω2 = 118.3247
Ω3 = 180.9463
Ω4 = 260.7555
Ω1 = 144.9682
Ω2 = 225.5369
Ω3 = 315.2119
Ω4 = 416.3045
Ω1 = 202.0838
Ω2 = 309.3531
Ω3 = 423.8021
Ω4 = 547.4178
Ω1 = 629.5080
Ω2 = 946.4871
Ω3 = 1266.124
Ω4 = 1589.284
Ω1 = 22.3733
Ω2 = 61.6728
Ω3 = 120.9034
Ω4 = 199.8595
Ω1 = 22.6464
Ω2 = 62.0450
Ω3 = 121.3117
Ω4 = 200.2883
Ω1 = 24.9574
Ω2 = 65.2921
Ω3 = 124.9250
Ω4 = 204.1063
Ω1 = 33.1763
Ω2 = 77.9997
Ω3 = 139.8260
Ω4 = 220.2625
Ω1 = 40.9932
Ω2 = 91.2643
Ω3 = 156.3955
Ω4 = 238.8921
Ω1 = 77.9113
Ω2 = 160.2803
Ω3 = 251.0408
Ω4 = 353.3270
Ω1 = 106.5809
Ω2 = 216.2878
Ω3 = 332.0612
Ω4 = 456.5178
Ω1 = 320.7300
Ω2 = 642.4135
Ω3 = 965.9985
Ω4 = 1292.423
Ω1 = 3.5160
Ω2 = 22.0345
Ω3 = 61.6972
Ω4 = 120.9019
Ω1 = 4.1102
Ω2 = 22.7566
Ω3 = 62.3205
Ω4 = 121.4914
Ω1 = 7.1675
Ω2 = 28.2944
Ω3 = 67.6582
Ω4 = 126.6700
Ω1 = 13.1150
Ω2 = 43.8214
Ω3 = 87.0217
Ω4 = 147.3205
Ω1 = 17.6105
Ω2 = 56.4993
Ω3 = 105.5011
Ω4 = 169.1492
Ω1 = 36.8506
Ω2 = 112.5007
Ω3 = 193.7980
Ω4 = 283.9516
Ω1 = 51.3544
Ω2 = 155.4788
Ω3 = 263.7737
Ω4 = 378.7869
Ω1 = 158.6856
Ω2 = 476.5181
Ω3 = 795.7319
Ω4 = 1117.240
Ω1 = 15.4182
Ω2 = 49.9649
Ω3 = 104.2477
Ω4 = 178.2697
Ω1 = 15.7868
Ω2 = 50.3923
Ω3 = 104.6977
Ω4 = 178.7317
Ω1 = 18.7601
Ω2 = 54.0847
Ω3 = 108.6638
Ω4 = 182.8365
Ω1 = 28.2268
Ω2 = 68.0614
Ω3 = 124.7580
Ω4 = 200.0577
Ω1 = 36.5823
Ω2 = 82.1562
Ω3 = 142.2974
Ω4 = 219.6758
Ω1 = 74.2657
Ω2 = 152.8053
Ω3 = 239.4464
Ω4 = 337.3108
Ω1 = 103.0991
Ω2 = 209.2184
Ω3 = 321.2171
Ω4 = 441.6689
Ω1 = 317.4898
Ω2 = 635.9209
Ω3 = 956.2297
Ω4 = 1279.343
Ω1 = 15.4182
Ω2 = 49.9649
Ω3 = 104.2477
Ω4 = 178.2697
Ω1 = 16.2749
Ω2 = 50.6718
Ω3 = 104.8926
Ω4 = 178.8810
Ω1 = 22.4012
Ω2 = 56.6107
Ω3 = 110.5203
Ω4 = 184.2890
Ω1 = 38.3118
Ω2 = 77.1052
Ω3 = 132.4469
Ω4 = 206.5280
Ω1 = 51.1261
Ω2 = 96.0450
Ω3 = 155.0711
Ω4 = 231.0744
Ω1 = 107.4790
Ω2 = 185.1940
Ω3 = 271.4863
Ω4 = 368.9759
Ω1 = 150.5600
Ω2 = 255.4300
Ω3 = 366.8290
Ω4 = 486.9893
Ω1 = 471.7514
Ω2 = 787.7684
Ω3 = 1106.052
Ω4 = 1427.492
Ω1 = 5.5933
Ω2 = 30.2258
Ω3 = 74.6389
Ω4 = 138.7913
Ω1 = 5.8612
Ω2 = 30.6321
Ω3 = 75.0797
Ω4 = 139.2481
Ω1 = 7.8368
Ω2 = 34.0648
Ω3 = 78.9352
Ω4 = 143.2935
Ω1 = 13.2454
Ω2 = 46.2359
Ω3 = 94.1508
Ω4 = 160.0307
Ω1 = 17.6680
Ω2 = 57.8092
Ω3 = 110.1907
Ω4 = 178.7345
Ω1 = 36.8599
Ω2 = 112.7474
Ω3 = 194.9002
Ω4 = 286.8104
Ω1 = 51.3588
Ω2 = 155.5970
Ω3 = 264.3142
Ω4 = 380.2390
Ω1 = 158.6860
Ω2 = 476.5290
Ω3 = 795.7823
Ω4 = 1117.378
Ω1 = 2.4674
Ω2 = 22.2066
Ω3 = 61.6850
Ω4 = 120.9027
Ω1 = 2.9250
Ω2 = 22.7011
Ω3 = 62.1830
Ω4 = 121.4016
Ω1 = 5.5463
Ω2 = 26.7432
Ω3 = 66.4973
Ω4 = 125.8033
Ω1 = 11.3780
Ω2 = 40.0433
Ω3 = 83.0018
Ω4 = 143.7449
Ω1 = 15.9006
Ω2 = 52.0941
Ω3 = 99.8676
Ω4 = 163.4250
Ω1 = 35.2106
Ω2 = 107.6868
Ω3 = 186.1385
Ω4 = 273.9868
Ω1 = 49.7342
Ω2 = 150.6643
Ω3 = 255.9103
Ω4 = 368.1306
Ω1 = 157.0990
Ω2 = 471.7618
Ω3 = 787.8168
Ω4 = 1106.184
Ω1 = 9.8696
Ω2 = 39.4784
Ω3 = 88.8264
Ω4 = 157.9137
Ω1 = 10.3575
Ω2 = 39.9753
Ω3 = 89.3250
Ω4 = 158.4129
Ω1 = 14.0038
Ω2 = 44.1965
Ω3 = 93.6931
Ω4 = 162.8369
Ω1 = 24.3082
Ω2 = 59.4346
Ω3 = 111.0471
Ω4 = 181.1972
Ω1 = 32.9298
Ω2 = 74.2050
Ω3 = 129.5098
Ω4 = 201.8120
Ω1 = 70.9381
Ω2 = 145.9375
Ω3 = 228.6993
Ω4 = 322.3252
Ω1 = 99.8349
Ω2 = 202.5758
Ω3 = 310.9929
Ω4 = 427.6101
Ω1 = 314.3143
Ω2 = 629.5576
Ω3 = 946.6544
Ω4 = 1266.520
Ω1 = 5.5933
Ω2 = 30.2258
Ω3 = 74.6389
Ω4 = 138.7913
Ω1 = 6.5991
Ω2 = 30.9884
Ω3 = 75.3094
Ω4 = 139.4172
Ω1 = 12.1336
Ω2 = 37.0879
Ω3 = 81.0816
Ω4 = 144.9247
Ω1 = 23.7382
Ω2 = 55.8711
Ω3 = 102.5205
Ω4 = 167.0796
Ω1 = 32.6296
Ω2 = 72.0629
Ω3 = 123.5789
Ω4 = 190.8295
Ω1 = 70.8762
Ω2 = 145.4485
Ω3 = 227.0920
Ω4 = 318.6875
Ω1 = 99.8040
Ω2 = 202.3290
Ω3 = 310.1673
Ω4 = 425.6855
Ω1 = 314.3112
Ω2 = 629.5328
Ω3 = 946.5707
Ω4 = 1266.322
Boundary
conditions
α = 0α = 1α = 10α = 50α = 100α = 500α = 1000α = 10,000
Ω1 = 9.8696
Ω2 = 39.4784
Ω3 = 88.8264
Ω4 = 157.9137
Ω1 = 10.3575
Ω2 = 39.9753
Ω3 = 89.3250
Ω4 = 158.4129
Ω1 = 14.0038
Ω2 = 44.1965
Ω3 = 93.6931
Ω4 = 162.8369
Ω1 = 24.3082
Ω2 = 59.4346
Ω3 = 111.0471
Ω4 = 181.1972
Ω1 = 32.9298
Ω2 = 74.2050
Ω3 = 129.5098
Ω4 = 201.8120
Ω1 = 70.9381
Ω2 = 145.9375
Ω3 = 228.6993
Ω4 = 322.3252
Ω1 = 99.8349
Ω2 = 202.5758
Ω3 = 310.9929
Ω4 = 427.6101
Ω1 = 314.3143
Ω2 = 629.5576
Ω3 = 946.6544
Ω4 = 1266.520
Ω1 = 22.3733
Ω2 = 61.6728
Ω3 = 120.9034
Ω4 = 199.8595
Ω1 = 23.4507
Ω2 = 62.5490
Ω3 = 121.6736
Ω4 = 200.5703
Ω1 = 31.3892
Ω2 = 69.8892
Ω3 = 128.3823
Ω4 = 206.8526
Ω1 = 52.8630
Ω2 = 95.1260
Ω3 = 154.3553
Ω4 = 232.5673
Ω1 = 70.1390
Ω2 = 118.3247
Ω3 = 180.9463
Ω4 = 260.7555
Ω1 = 144.9682
Ω2 = 225.5369
Ω3 = 315.2119
Ω4 = 416.3045
Ω1 = 202.0838
Ω2 = 309.3531
Ω3 = 423.8021
Ω4 = 547.4178
Ω1 = 629.5080
Ω2 = 946.4871
Ω3 = 1266.124
Ω4 = 1589.284
Ω1 = 22.3733
Ω2 = 61.6728
Ω3 = 120.9034
Ω4 = 199.8595
Ω1 = 22.6464
Ω2 = 62.0450
Ω3 = 121.3117
Ω4 = 200.2883
Ω1 = 24.9574
Ω2 = 65.2921
Ω3 = 124.9250
Ω4 = 204.1063
Ω1 = 33.1763
Ω2 = 77.9997
Ω3 = 139.8260
Ω4 = 220.2625
Ω1 = 40.9932
Ω2 = 91.2643
Ω3 = 156.3955
Ω4 = 238.8921
Ω1 = 77.9113
Ω2 = 160.2803
Ω3 = 251.0408
Ω4 = 353.3270
Ω1 = 106.5809
Ω2 = 216.2878
Ω3 = 332.0612
Ω4 = 456.5178
Ω1 = 320.7300
Ω2 = 642.4135
Ω3 = 965.9985
Ω4 = 1292.423
Ω1 = 3.5160
Ω2 = 22.0345
Ω3 = 61.6972
Ω4 = 120.9019
Ω1 = 4.1102
Ω2 = 22.7566
Ω3 = 62.3205
Ω4 = 121.4914
Ω1 = 7.1675
Ω2 = 28.2944
Ω3 = 67.6582
Ω4 = 126.6700
Ω1 = 13.1150
Ω2 = 43.8214
Ω3 = 87.0217
Ω4 = 147.3205
Ω1 = 17.6105
Ω2 = 56.4993
Ω3 = 105.5011
Ω4 = 169.1492
Ω1 = 36.8506
Ω2 = 112.5007
Ω3 = 193.7980
Ω4 = 283.9516
Ω1 = 51.3544
Ω2 = 155.4788
Ω3 = 263.7737
Ω4 = 378.7869
Ω1 = 158.6856
Ω2 = 476.5181
Ω3 = 795.7319
Ω4 = 1117.240
Ω1 = 15.4182
Ω2 = 49.9649
Ω3 = 104.2477
Ω4 = 178.2697
Ω1 = 15.7868
Ω2 = 50.3923
Ω3 = 104.6977
Ω4 = 178.7317
Ω1 = 18.7601
Ω2 = 54.0847
Ω3 = 108.6638
Ω4 = 182.8365
Ω1 = 28.2268
Ω2 = 68.0614
Ω3 = 124.7580
Ω4 = 200.0577
Ω1 = 36.5823
Ω2 = 82.1562
Ω3 = 142.2974
Ω4 = 219.6758
Ω1 = 74.2657
Ω2 = 152.8053
Ω3 = 239.4464
Ω4 = 337.3108
Ω1 = 103.0991
Ω2 = 209.2184
Ω3 = 321.2171
Ω4 = 441.6689
Ω1 = 317.4898
Ω2 = 635.9209
Ω3 = 956.2297
Ω4 = 1279.343
Ω1 = 15.4182
Ω2 = 49.9649
Ω3 = 104.2477
Ω4 = 178.2697
Ω1 = 16.2749
Ω2 = 50.6718
Ω3 = 104.8926
Ω4 = 178.8810
Ω1 = 22.4012
Ω2 = 56.6107
Ω3 = 110.5203
Ω4 = 184.2890
Ω1 = 38.3118
Ω2 = 77.1052
Ω3 = 132.4469
Ω4 = 206.5280
Ω1 = 51.1261
Ω2 = 96.0450
Ω3 = 155.0711
Ω4 = 231.0744
Ω1 = 107.4790
Ω2 = 185.1940
Ω3 = 271.4863
Ω4 = 368.9759
Ω1 = 150.5600
Ω2 = 255.4300
Ω3 = 366.8290
Ω4 = 486.9893
Ω1 = 471.7514
Ω2 = 787.7684
Ω3 = 1106.052
Ω4 = 1427.492
Ω1 = 5.5933
Ω2 = 30.2258
Ω3 = 74.6389
Ω4 = 138.7913
Ω1 = 5.8612
Ω2 = 30.6321
Ω3 = 75.0797
Ω4 = 139.2481
Ω1 = 7.8368
Ω2 = 34.0648
Ω3 = 78.9352
Ω4 = 143.2935
Ω1 = 13.2454
Ω2 = 46.2359
Ω3 = 94.1508
Ω4 = 160.0307
Ω1 = 17.6680
Ω2 = 57.8092
Ω3 = 110.1907
Ω4 = 178.7345
Ω1 = 36.8599
Ω2 = 112.7474
Ω3 = 194.9002
Ω4 = 286.8104
Ω1 = 51.3588
Ω2 = 155.5970
Ω3 = 264.3142
Ω4 = 380.2390
Ω1 = 158.6860
Ω2 = 476.5290
Ω3 = 795.7823
Ω4 = 1117.378
Ω1 = 2.4674
Ω2 = 22.2066
Ω3 = 61.6850
Ω4 = 120.9027
Ω1 = 2.9250
Ω2 = 22.7011
Ω3 = 62.1830
Ω4 = 121.4016
Ω1 = 5.5463
Ω2 = 26.7432
Ω3 = 66.4973
Ω4 = 125.8033
Ω1 = 11.3780
Ω2 = 40.0433
Ω3 = 83.0018
Ω4 = 143.7449
Ω1 = 15.9006
Ω2 = 52.0941
Ω3 = 99.8676
Ω4 = 163.4250
Ω1 = 35.2106
Ω2 = 107.6868
Ω3 = 186.1385
Ω4 = 273.9868
Ω1 = 49.7342
Ω2 = 150.6643
Ω3 = 255.9103
Ω4 = 368.1306
Ω1 = 157.0990
Ω2 = 471.7618
Ω3 = 787.8168
Ω4 = 1106.184
Ω1 = 9.8696
Ω2 = 39.4784
Ω3 = 88.8264
Ω4 = 157.9137
Ω1 = 10.3575
Ω2 = 39.9753
Ω3 = 89.3250
Ω4 = 158.4129
Ω1 = 14.0038
Ω2 = 44.1965
Ω3 = 93.6931
Ω4 = 162.8369
Ω1 = 24.3082
Ω2 = 59.4346
Ω3 = 111.0471
Ω4 = 181.1972
Ω1 = 32.9298
Ω2 = 74.2050
Ω3 = 129.5098
Ω4 = 201.8120
Ω1 = 70.9381
Ω2 = 145.9375
Ω3 = 228.6993
Ω4 = 322.3252
Ω1 = 99.8349
Ω2 = 202.5758
Ω3 = 310.9929
Ω4 = 427.6101
Ω1 = 314.3143
Ω2 = 629.5576
Ω3 = 946.6544
Ω4 = 1266.520
Ω1 = 5.5933
Ω2 = 30.2258
Ω3 = 74.6389
Ω4 = 138.7913
Ω1 = 6.5991
Ω2 = 30.9884
Ω3 = 75.3094
Ω4 = 139.4172
Ω1 = 12.1336
Ω2 = 37.0879
Ω3 = 81.0816
Ω4 = 144.9247
Ω1 = 23.7382
Ω2 = 55.8711
Ω3 = 102.5205
Ω4 = 167.0796
Ω1 = 32.6296
Ω2 = 72.0629
Ω3 = 123.5789
Ω4 = 190.8295
Ω1 = 70.8762
Ω2 = 145.4485
Ω3 = 227.0920
Ω4 = 318.6875
Ω1 = 99.8040
Ω2 = 202.3290
Ω3 = 310.1673
Ω4 = 425.6855
Ω1 = 314.3112
Ω2 = 629.5328
Ω3 = 946.5707
Ω4 = 1266.322

#### 3.1.1 Effect of Axial Force.

The effect of axial loading on the dimensionless frequency may be observed in Table 1, which lists the dimensionless frequencies for the first four modes of the previously described beams with different values of the dimensionless axial loading parameter α. The dimensionless frequencies, and thus, the dimensional natural frequencies tend to increase as the dimensionless axial load parameter α increases. This trend is ubiquitous for all observed boundary conditions and modes. It is well known that a beam begins to behave like a string when either the tensile force P becomes very large or the flexural rigidity EI becomes very small [31,32]. This trend is illustrated in Fig. 3 with a plot adapted from Wei et al. [34]. The figure plots the dimensionless frequency for the fundamental mode Ω1 versus the dimensionless axial load parameter for a string, which behaves as Ωn = nπα1/2, and three beams with different boundary conditions. From Fig. 3, at small values of α, the effect of axial loading is negligible and the dimensionless frequencies closely resemble those of a beam with no loading. This may be confirmed by comparing Fig. 3 and Table 1. However, as α becomes larger, the axial load begins to dominate and the response of the beam becomes independent of flexural rigidity. This trend is depicted in Fig. 3 where the three beams closely follow the same dependence of a string for α > 104. The proposed method is developed and intended for beams with observable dependence on axial loading and flexural rigidity. Examples and results presented later will be for beams with a dimensionless axial loading parameter 100 ≤ α ≤ 105.

Fig. 3
Fig. 3
Close modal

#### 3.1.2 Effect of Flexible Boundary Conditions.

In this section, the effect of flexible boundary conditions on the dimensionless frequency and response of the beam will be considered. To better examine these effects, an unloaded beam is considered such that α = 0. First, a beam with translational springs at both ends, modeled with the dimensionless parameter $β$ and $K~t$, and fixed rotations at the ends is considered. The dimensionless frequency of the fundamental mode Ω1 is plotted versus the dimensionless translational stiffness parameter $β$ in Fig. 4(a). At small $β$, the low stiffness of the boundary springs dominates and the response follows that of a mass-spring system, where the mass is the mass of the entire beam and the effective stiffness is 2kt, which results in the relation $Ω1=2β$. This trend is also observed in Fig. 4(b), which plots the fundamental mode shapes for different values of $β$. From the figure, for $β=101$, the mode is largely dominated by the extension of the boundary springs. As $β$ increases, the flexural rigidity begins to dominate and the beam begins to experience more bending and less extension at the boundary springs. This trend continues until $β$ becomes very large and the response follows that of a beam with fixed–fixed boundary conditions. This is observed in Fig. 4(a), where the dimensionless frequency approaches the dimensionless frequency of a fixed–fixed beam, Ω1 = 22.3733, as listed in Table 1. This trend and dependence on $β$ are followed for all modes, as depicted in Fig. 4(c), which plots the dimensionless frequency for the first four modes. Considering the beams of interest, the present work will be limited to translational stiffnesses captured in the range $101≤β≤105$.

Fig. 4
Fig. 4
Close modal

Finally, an unloaded beam with rotational springs at both ends, modeled with $χ$ and $K~r$, and fixed translation at the ends is considered. The dimensionless frequency of the fundamental mode Ω1 is plotted versus the dimensionless rotational stiffness parameter $χ$ in Fig. 5(a). As would be expected, at low values of $χ$, the presence of the rotational springs is negligible and the beam behaves like a beam with pinned–pinned boundary conditions such that Ω1 = 9.8696, as listed in Table 1. However, as $χ$ increases, the dimensionless frequency is observed to depend more strongly on the rotational spring until the beam begins to behave like a fixed–fixed beam. This is depicted in Fig. 5(a) with a transition region from Ω1 = 9.8696 to Ω1 = 22.3733. This trend is also exhibited in all modes as observed in Fig. 5(b), which plots the dimensionless frequency for the first four modes versus the dimensionless rotational stiffness. In comparison to the translational stiffness, the rotational stiffness is observed to have a weaker effect on the response of the beam and dimensionless frequency. The present work will be limited to rotational stiffnesses captured in the range $10−2≤χ≤104$.

Fig. 5
Fig. 5
Close modal

## 4 Proposed Method

In this section, the proposed method for simultaneously predicting the elastic modulus, axial load, and boundary conditions from measured natural frequencies is presented. The method assumes that the dimensions and density of the beam are known and at least two natural frequencies have been measured. The method predicts the initially uncertain beam properties by finding values of the dimensionless parameters, α, $β$, and $χ$, that minimize the error between the measured and predicted natural frequencies. The predicted natural frequencies are found with the dimensionless frequencies Ωn, which are evaluated with the generalized eigenvalue problem. A scaling factor is then computed by comparing the fundamental dimensionless frequency to the measured fundamental frequency. It should be noted that during this step, the elastic modulus of the beam is simultaneously predicted. With the scaling factor, the natural frequencies of the beam are computed and a normalized error between the predicted and measured natural frequencies is evaluated. The uncertain beam properties are predicted by finding the dimensionless parameters that minimize the normalized error. The proposed method is outlined in detail with the following steps.

For each choice of α, $β$, and $χ$:

1. Compute eigenvalues $λn$ by evaluating the generalized eigenvalue problem in Eq. (15) with generalized matrices defined in Eq. (16). With the eigenvalues, compute the dimensionless frequencies Ωn with Eq. (18).

2. Compute the scaling factor using the measured fundamental frequency $ω1,act$ and evaluated dimensionless fundamental frequency Ω1 with the relationship
$s=ω1,actΩ1$
(20)
It should be noted that by computing the scaling factor s, the elastic modulus of the beam is predicted by virtue of Eq. (14d) and the relationship between the dimensionless frequency and dimensional natural frequency. Comparing Eqs. (14d) and (20), the scaling factor is related by $s=EI/ρAL4$.
3. Compute predicted natural frequencies with scaling factor and dimensionless frequency with
$ωn,pred=sΩn$
(21)
4. Compute normalized error between predicted and measured natural frequencies with
$ε=rms(ωact−ωpred)rms(ωact)$
(22)
where $ωact$ and $ωpred$ are vectors of the same length containing measured and predicted natural frequencies, respectively, and rms(·) is the root-mean-square of the vector argument.

The initially uncertain beam properties are predicted by iterating through values of α, $β$, and $χ$ in search for values that minimize the error in Eq. (22). When the error between predicted and measured natural frequencies is minimized, the beam properties are predicted with the following steps.

For values $α^$, $β^$, and $χ^$ that minimize the normalized error in Eq. (22):

1. Compute eigenvalues $λ^n$ by evaluating the generalized eigenvalue problem in Eq. (15) with generalized matrices defined in Eq. (16). With the eigenvalues, compute the dimensionless frequencies $Ω^n$ with Eq. (18).

2. Compute the scaling factor $s^$ using the measured fundamental frequency $ω1,act$ and dimensionless fundamental frequency $Ω^1$ with Eq. (20).

3. Evaluate predicted elastic modulus with computed scaling factor $s^$ and known dimensions and density with
$E^=s^2(ρAL4I)$
(23)
4. With predicted elastic modulus and dimensionless parameters $α^$, $β^$, and $χ^$, evaluate axial load and boundary conditions from
$P^=α^(E^IL2)$
(24a)
$k^t=β^(E^IL3)$
(24b)
$k^r=χ^(E^IL)$
(24c)

The proposed method is advantageous to previous methods as it simultaneously predicts the elastic modulus from the scaling factor while iterating through the dimensionless parameters. This will significantly reduce the computational cost of any search algorithm and make brute force scans more feasible, as there will be one less parameter that must be optimized. Additionally, although the method was derived for a beam with boundary conditions depicted in Fig. 2, the FEM and the proposed method may be easily adapted to predict properties for a wide range of beams. This may be done by imposing known ideal boundary conditions with the method presented by Wu et al. [33]. Furthermore, although it is not within the scope of this paper, flexible boundary conditions that are not symmetric may be modeled and predicted by introducing additional stiffness terms in Eq. (8) that will allow for the variation in stiffness between left and right ends.

### 4.1 Required Measurements.

In the proposed method, a minimum of two measured natural frequencies is required. Since the scaling factor is computed with a measured natural frequency in Eq. (20), an additional measured frequency is needed to compute the normalized error with Eq. (22). If only one measured natural frequency is available, then the normalized error will be independent of the iterated values α, $β$, and $χ$, and equal to ɛ = 0. By measuring a second natural frequency, an error between predicted and measured natural frequencies is characterized and the uncertain beam parameters may be accurately predicted. Although only two natural frequencies are required, the proposed method generally performs better when more measurements are available. The dependence of the method’s accuracy on the number of measured natural frequencies will be examined with numerical and experimental examples. It should also be noted that any natural frequency and dimensionless frequency may be used to compute the scaling factor s, since the relationship in Eq. (14d) is valid for all modes. In Sec. 5, the proposed method is illustrated with numerical examples that consider NEMS beams and a range of uncertain parameters.

## 5 Numerical Example

In this section, numerical examples are presented to illustrate the accuracy of the method with cases in which actual beam parameters are known. In all examples presented, a NEMS beam with the dimensions and density listed in Table 2 is considered. Here, the density is that of silicon nitride [35]. The aspect ratios are (b/L) = 0.019 and (h/L) = 0.002 and are within the acceptable range for which the Euler–Bernoulli FEM beam model is valid. In the examples, that follow, the highest mode examined is the fifth mode. Consequently, the FEM beam model is discretized into Q = 200 elements to satisfy the general criterion of six elements per wavelength [36]. In Secs. 5.15.3, three different examples are examined that consider different initially uncertain parameters. The actual beam properties for all examples and cases are listed in Table 3.

Table 2

Dimensions and density of NEMS beam

 L Length 50 μm b Width 950 nm h Thickness 100 nm I Second moment of area 7.9167 × 10−29 m4 $ρ$ Density 3100 kg/m3
 L Length 50 μm b Width 950 nm h Thickness 100 nm I Second moment of area 7.9167 × 10−29 m4 $ρ$ Density 3100 kg/m3
Table 3

Actual beam properties for numerical examples

CaseE (GPa)αP (nN)$log10(β)$kt (N/m)$χ$kr (nN m)
Ax. 12701085.5
Ax. 227050427.5
Ax. 3270100855
Ax. 42705004275
Ax. Tr. 127010085520.0171
Ax. Tr. 227010085530.1710
Ax. Tr. 32702502137.520.0171
Ax. Tr. 42702502137.530.1710
Ax. Ro. 1270100855500.0214
Ax. Ro. 2270100855750.0321
Ax. Ro. 32702502137.5500.0214
Ax. Ro. 42702502137.5750.0321
CaseE (GPa)αP (nN)$log10(β)$kt (N/m)$χ$kr (nN m)
Ax. 12701085.5
Ax. 227050427.5
Ax. 3270100855
Ax. 42705004275
Ax. Tr. 127010085520.0171
Ax. Tr. 227010085530.1710
Ax. Tr. 32702502137.520.0171
Ax. Tr. 42702502137.530.1710
Ax. Ro. 1270100855500.0214
Ax. Ro. 2270100855750.0321
Ax. Ro. 32702502137.5500.0214
Ax. Ro. 42702502137.5750.0321

### 5.1 Uncertain Elastic Modulus and Axial Load.

In the first example, a beam with ideal fixed–fixed boundary conditions and an uncertain elastic modulus and axial load is considered. The boundary conditions are assumed to be known and modeled with $β=χ=0$ in Eq. (16a) and the method outlined in Ref. [33]. Four cases, Ax. 1–Ax. 4, with different actual values of the dimensionless axial load parameter, α, are examined. To predict the elastic modulus and axial load, the steps outlined in Sec. 4 were employed. A global optimization search was performed with matlab’s GlobalSearch [37] with the local solver fmincon and the interior point algorithm [38] to iterate through values of α. The objective function to be minimized was the normalized error ɛ described in Eq. (22), where the vectors $ωact$ and $ωpred$ were composed of the first two simulated and predicted natural frequencies, respectively. Again, it is important to note that the proposed method only requires the axial load parameter α to be optimized, as the elastic modulus is simultaneously predicted from the scaling factor s, which is evaluated at each iteration.

The dimensionless axial load parameter $α^$ that minimized the normalized error in Eq. (22), the predicted elastic modulus $E^$, and axial load $P^$ are listed in Table 5 for all cases examined. The normalized error between the actual and predicted properties is also listed in the table. From Table 5, it is observed that the predicted elastic modulus and axial load agree well with the actual properties, as all errors are well below $1%$ and an exact prediction was found in case Ax. 4. Also important to observe is the normalized error between the actual and predicted natural frequencies for the first two modes listed in Table 4, where the error is denoted as ɛa. To recall earlier in the example, this error was used as the objective function that was minimized to predict the initially uncertain beam properties. Consequently, a low value is expected for this error. For a better representation of the accuracy of the proposed method and predicted properties, the normalized error in predicted natural frequencies for higher modes, not included in the optimization search, is computed. This error is labeled as ɛb in Table 4 and is the normalized error in natural frequencies for the third to fifth mode. From the results, it is observed that the error over the higher modes ɛb is slightly larger than the error over the lower modes ɛa. This increase in error is due to the fact that the optimization search decreases the error between actual and predicted frequencies at only the first two modes. Thus, the predicted beam properties minimize the error at those two particular frequencies, as opposed to the error at higher frequencies. Although the error over the higher frequencies ɛb is slightly larger than the error over the fit frequencies ɛa, it is important to note that they are still reasonably small. These results indicate that the predicted beam properties accurately model the response of the beam at higher modes that were not included in the initial fit and highlight the proposed method’s ability to predict beam properties from as little as two measured natural frequencies.

Table 4

Normalized error between actual and predicted natural frequencies for modes one and two, ɛa, and modes three to five, ɛb

Caseɛaɛb
Ax. 13.6503 × 10−47.1217 × 10−4
Ax. 22.4340 × 10−76.2649 × 10−7
Ax. 37.0953 × 10−62.2150 × 10−5
Ax. 40.00000.0000
Caseɛaɛb
Ax. 13.6503 × 10−47.1217 × 10−4
Ax. 22.4340 × 10−76.2649 × 10−7
Ax. 37.0953 × 10−62.2150 × 10−5
Ax. 40.00000.0000
Table 5

Predicted beam properties and corresponding normalized error for numerical examples

Case$E^$ (GPa)Error$α^$$P^$ (nN)Error$log10(β^)$$k^t$ (N/m)Error$χ^$$k^r$ (nN · m)Error
Ax. 1269.53021.7400 × 10−310.089386.11327.1715 × 10−3
Ax. 2269.99951.7684 × 10−650.0002427.50071.5705 × 10−6
Ax. 3270.01967.2657 × 10−599.9893854.97103.3928 × 10−5
Ax. 4270.00000.0000500.000042750.0000
Ax.Tr. 1249.90537.4425 × 10−2108.8095861.08237.1138 × 10−32.03441.7134 × 10−21.9799 × 10−3
Ax.Tr. 2276.76892.5070 × 10−295.3850835.98722.2237 × 10−22.99281.7242 × 10−18.3162 × 10−3
Ax.Tr. 3249.74377.5023 × 10−2271.07412143.80262.9486 × 10−32.03411.7111 × 10−26.3219 × 10−4
Ax.Tr. 4272.74801.0178 × 10−2247.42312136.99902.3438 × 10−42.99521.7083 × 10−11.0044 × 10−3
Ax.Ro. 1275.74802.1289 × 10−299.5219869.02701.6406 × 10−239.32001.7167 × 10−21.9686 × 10−1
Ax.Ro. 2271.41285.2326 × 10−399.8245857.96543.4683 × 10−369.61332.9915 × 10−26.6966 × 10−2
Ax.Ro. 3276.56262.4306 × 10−2248.45372175.91161.7970 × 10−235.72281.5643 × 10−22.6818 × 10−1
Ax.Ro. 4271.19864.4391 × 10−3249.59762143.53302.8225 × 10−369.72282.9939 × 10−26.6236 × 10−2
Case$E^$ (GPa)Error$α^$$P^$ (nN)Error$log10(β^)$$k^t$ (N/m)Error$χ^$$k^r$ (nN · m)Error
Ax. 1269.53021.7400 × 10−310.089386.11327.1715 × 10−3
Ax. 2269.99951.7684 × 10−650.0002427.50071.5705 × 10−6
Ax. 3270.01967.2657 × 10−599.9893854.97103.3928 × 10−5
Ax. 4270.00000.0000500.000042750.0000
Ax.Tr. 1249.90537.4425 × 10−2108.8095861.08237.1138 × 10−32.03441.7134 × 10−21.9799 × 10−3
Ax.Tr. 2276.76892.5070 × 10−295.3850835.98722.2237 × 10−22.99281.7242 × 10−18.3162 × 10−3
Ax.Tr. 3249.74377.5023 × 10−2271.07412143.80262.9486 × 10−32.03411.7111 × 10−26.3219 × 10−4
Ax.Tr. 4272.74801.0178 × 10−2247.42312136.99902.3438 × 10−42.99521.7083 × 10−11.0044 × 10−3
Ax.Ro. 1275.74802.1289 × 10−299.5219869.02701.6406 × 10−239.32001.7167 × 10−21.9686 × 10−1
Ax.Ro. 2271.41285.2326 × 10−399.8245857.96543.4683 × 10−369.61332.9915 × 10−26.6966 × 10−2
Ax.Ro. 3276.56262.4306 × 10−2248.45372175.91161.7970 × 10−235.72281.5643 × 10−22.6818 × 10−1
Ax.Ro. 4271.19864.4391 × 10−3249.59762143.53302.8225 × 10−369.72282.9939 × 10−26.6236 × 10−2

As previously mentioned, the first two natural frequencies were used to evaluate the objective function. These frequencies were chosen with the consideration that lower modes are often easier to excite and measure. However, it is worth noting that the proposed method will perform best when using frequencies that are most sensitive to changes in the parameters of interest. Specifically, measured frequencies that will produce large gradients in normalized error will be advantageous for predicting beam properties with global optimization searches or brute force scans. For example, Fig. 6 plots the normalized error in predicted natural frequency for case Ax. 1 versus dimensionless axial load parameter α for different pairs of fit natural frequencies. For each case, the lowest frequency was used to compute the scaling factor s. From Fig. 6, it is clear that for each case analyzed, a minimum in error occurs at α = 10, corresponding to the actual beam property for case Ax. 1. These results imply that any set of measured natural frequencies may be compared to predicted frequencies to compute an error, which when minimized will accurately predict the properties of a beam. However, Fig. 6 also illustrates the different sensitivities for various pairs of fit frequencies. Relatively large gradients in error are observed when the first natural frequency is included in the fit. Furthermore, the gradient is observed to increase with an increase in the second fit frequency. These results suggest that using nonconsecutive natural frequencies is advantageous for generating a well-defined error profile. It is important to note that these results are specific to this example and the beam’s dependence on the axial load. To determine optimal sets of fit natural frequencies, detailed sensitivity studies [39] are required, which are out of the scope of this paper. In all examples, to follow, consecutive natural frequencies starting at the fundamental are used to compute the normalized error between predicted and measured natural frequencies.

Fig. 6
Fig. 6
Close modal

### 5.2 Uncertain Elastic Modulus, Axial Load, and Translational Stiffness.

In the examples presented in this section, the NEMS beam under consideration has an uncertain elastic modulus, axial load, and translational stiffness at the boundaries. The rotational degrees-of-freedom are taken to be fixed at both ends and assumed to be known. The actual values of the elastic modulus, axial load, and translational stiffness for four different cases, Ax. Tr. 1–Ax. Tr. 4, are listed in Table 3. Similar to the previous example, the uncertain parameters are predicted using a global optimization search that optimizes the dimensionless parameters α and $β$, while simultaneously predicting the elastic modulus. Specifically, the logarithm of the dimensionless translational stiffness $log10(β)$ was set as a search parameter to account for the relatively large numerical range. For this example, the first two natural frequencies were used to compute the objective function in Eq. (22).

The predicted beam properties are listed in Table 5 with their corresponding normalized error. Although the errors are slightly larger than those observed in the previous example, the results suggest that the proposed method is able to accurately predict the elastic modulus, axial load, and boundary conditions simultaneously. From the results in Table 5, all beam properties were predicted with an error less than $10%$. Furthermore, for cases Ax. Tr. 2 and 4, all beam properties were predicted with an error less than $5%$. From Table 5, relatively large errors in the elastic modulus were observed in cases Ax. Tr. 1 and Ax. Tr. 3., which corresponded to the cases where the actual translational stiffness was kt = 0.0171 N/m with $log10(β)=2$. To better interpret these results, the logarithm of the normalized error in predicted frequency log10(ɛ) is plotted versus the dimensionless axial load parameter α and transnational stiffness parameter $β$ in Figs. 7(a)7(d) for four different cases. Comparing Figs. 7(a) and 7(c) with Figs. 7(b) and 7(d), less variation in error and a shallower contour is observed for cases Ax. Tr. 1 and Ax. Tr. 3. These trends suggest that the predicted beam parameters in cases Ax. Tr. 1 and Ax. Tr. 3 result in small errors in predicted natural frequencies but do not represent the global minimum. As a result, small errors in predicted beam properties may arise when an optimization search is employed.

Fig. 7
Fig. 7
Close modal

The results in Table 5 illustrate that the proposed method is able to simultaneously predict three unknown parameters, elastic modulus, axial load, and translational stiffness, from as little as two natural frequencies. As discussed in Sec. 4.1, the two frequencies are required to characterize a normalized error between measured and predicted natural frequencies. From Figs. 7(a)7(d), a global minimum occurs when the uncertain beam properties are evaluated at values equivalent to the actual properties. This unique global minimum allows for such beam properties to be predicted. However, from Figs. 7(a)7(d), additional local minima exist resulting in the observed error in predicted properties. To obtain a smoother error profile and reduce the number of local minima, the number of measured natural frequencies included in the fit should be increased. This trend is illustrated in Figs. 7(e)7(h), which plot the logarithm of the normalized error computed using the first three natural frequencies for cases Ax. Tr. 1–Ax. Tr. 4. Comparing Figs. 7(a)7(d) with Figs. 7(e)7(h), the increase in fit natural frequencies results in a more well-defined global minimum. This can be clearly observed in the results for cases Ax. Tr. 1 and Ax. Tr. 3, where the two frequency fit, in Figs. 7(a) and 7(c), resulted in many local minima and very small gradients near the global minimum. In comparison, the three frequency fit, in Figs. 7(e) and 7(g), resulted in fewer local minima and a clearer global minimum. These results suggest that increasing the number of fit natural frequencies will give rise to error profiles that will better characterize the accurate beam properties. As a result, the proposed method generally performs better when more measured natural frequencies are available to compare to predicted frequencies. Additional results and analysis on the effect the number of fit natural frequencies has on the proposed method will be presented in Sec. 6.

### 5.3 Uncertain Elastic Modulus, Axial Load, and Rotational Stiffness.

In the final numerical example, the NEMS beam under consideration has an uncertain elastic modulus, axial load, and rotational stiffness at the boundaries. The translational degrees-of-freedom are taken to be fixed at both ends and assumed to be known. The actual beam properties for four different cases, Ax. Ro. 1–Ax. Ro. 4, are listed in Table 3. In this example, the first two natural frequencies of the beam are taken to be measured, and the normalized error in Eq. (22) is minimized with a global optimization search that optimizes the dimensionless axial load parameter α and rotational stiffness parameter $χ$. The predicted properties and their corresponding error are listed in Table 5.

From the results, it is observed that the elastic modulus and axial load are predicted accurately for all cases, with errors well below 5%. However, it is apparent that the proposed method with the global optimization search has difficulty predicting the rotational stiffness, with a relatively large error of $26.818%$ in case Ax. Ro. 3. The difference in errors is primarily due to the fact that the response of the beam is fairly insensitive to changes in the rotational stiffness. This trend may be observed in Fig. 5(b), where the dimensionless frequencies Ωn experience relatively small change with change in the dimensionless rotational stiffness $χ$. Additionally, similar trends may be observed in Fig. 8, which plots the logarithm of the normalized error in predicted natural frequencies versus the dimensionless axial load and rotational stiffness parameters. From the contour plots, the error appears to be relatively insensitive to changes in the dimensionless rotational stiffness with very shallow error regions and small variations. As a result, the proposed method has a difficult time finding the accurate dimensionless rotational stiffness parameter from the error in predicted natural frequencies. From the numerical examples presented in this section, it is concluded that the proposed method is limited to accurately predicting the elastic modulus, axial load, and translational boundary conditions. When predicting the rotational stiffness, the proposed method is capable of reducing the error of the predicted natural frequencies. However, due to the observed insensitivity, the predicted rotational stiffness will be prone to relatively large errors.

Fig. 8
Fig. 8
Close modal

## 6 Experimental Example

In this section, the proposed method is applied to an example with experimental measurements taken from a NEMS beam with uncertain beam properties. A scanning electron microscope (SEM) image of the beam is shown in Fig. 9(a) with a top-view schematic of the beam, not drawn to scale, depicted in Fig. 9(b). The dimensions and density of the NEMS beam are identical to those listed in Table 2. With the fabrication procedure presented by Ari et al. [40], the NEMS beam was fabricated with a standard top–down approach starting with a 100-nm low-pressure chemical vapor deposition (LPCVD) silicon nitride-coated silicon wafer. The exact elastic modulus and axial load of the film and beam are unknown. In the experiment, the beam was excited using electrothermal actuation by applying an AC voltage to the gold electrode [41]. The applied voltage resulted in localized heating and a nonuniform expansion, attributed to the difference in linear thermal coefficients of the gold electrode and silicon nitride, which resulted in the flexing of the beam. For the first 11 modes, the displacement of the beam was measured optically with a path-stabilized Michelson interferometer at 800 discrete frequencies near resonance. A modal fit was then performed to find the first 11 natural frequencies of the NEMS beam [42]. In the following sections, the proposed method is used to predict the elastic modulus, axial load, and boundary conditions of the beam with the measured natural frequencies.

Fig. 9
Fig. 9
Close modal

### 6.1 Ideal Boundary Conditions.

First, the NEMS beam is modeled with ideal fixed–fixed boundary conditions, and initially, uncertain elastic modulus and axial load are simultaneously predicted with the proposed method. The FEM beam model was discretized with Q = 200 elements to again satisfy the general criterion of six elements per wavelength. In this model, the entire 50 μm NEMS beam, which includes sections attached to the gold electrodes, as pictured in Fig. 9(b), is assumed to be homogeneous with a constant elastic modulus, density, axial load, and cross-sectional area. Similar to the numerical example in Sec. 5.1, the boundary conditions are modeled by setting $β=χ=0$ and fixing the translational and rotational degrees-of-freedoms at both ends. As outlined in Sec. 4, the proposed method then iterated through values of the dimensionless axial load parameter α in search for a value that minimized the normalized error between predicted and measured natural frequencies. The normalized error was evaluated using all 11 measured natural frequencies. For this example, the iterated values of the dimensionless axial load parameter were evaluated from a set of discrete values in the range of 10−4 ≤ α ≤ 106 to ensure the global minimum was captured within reasonable values. The evaluated error is plotted in Fig. 10(a) where a clear global minimum is observed at $α^=1.4940×103$. With $α^$ and the steps outlined in Sec. 4, the predicted elastic modulus and axial load were evaluated with Eqs. (23) and (24a), and their values are listed in Table 6.

Fig. 10
Fig. 10
Close modal
Table 6

Predicted beam properties of experimental NEMS beam for ideal and flexible boundary conditions

Boundary conditions$E^$ (GPa)$α^$$P^(μN)$$β^$$k^t$ (N/m)
Fixed–fixed194.63881493.99359.2083
Translational spring212.7958953.00567.93654.3545 × 1048.0899
Boundary conditions$E^$ (GPa)$α^$$P^(μN)$$β^$$k^t$ (N/m)
Fixed–fixed194.63881493.99359.2083
Translational spring212.7958953.00567.93654.3545 × 1048.0899

To evaluate the accuracy of the predicted properties, the predicted natural frequencies are compared with the measured natural frequencies in Fig. 10(b), where good agreement is observed. The resulting normalized error in the predicted natural frequencies is ɛ = 3.12678 × 10−3. These results suggest that the predicted axial load and elastic modulus accurately model the response of the NEMS beam. For further validation, the predicted elastic modulus may be compared to values reported in the literature for LPCVD silicon nitride, which vary between 146 and 290 GPa [4346]. From the results in Table 6, the predicted elastic modulus $E^=194.6388GPa$ is of reasonable magnitude with respect to previously reported values. Additionally, from the predicted axial load, the normal tensile stress of the beam is predicted to be $σ=96.9297MPa,$ which falls within the range of reported values of the tensile stress of low stress silicon nitride wafers, which typically have an upper bound of 1000 MPa [4749].

### 6.2 Flexible Boundary Conditions.

To improve the model of the NEMS beam and agreement between measured and predicted natural frequencies, the beam is modeled with translational springs at both ends. Specifically, the homogeneous section of silicon between the two gold electrodes with length L = 45 μm, pictured in Fig. 9(c), is modeled with two translational springs and fixed rotations at both ends. The resulting FEM depicted in Fig. 9(d) reduces previous assumptions by modeling a section of the beam that is approximately homogeneous in material properties and dimensions, while accounting for the stiffness presented by the cantilever ends with translational springs kt. It should be noted that rotational springs are not included because it was found that the response of the beam was highly insensitive to the added rotational stiffnesses, as previously observed in Sec. 5.3. Specifically, it was found that the error between predicted and measured natural frequencies did not vary with changes in the rotational spring constant. Furthermore, it was observed that including rotational springs in the model did not improve the accuracy of the model and predicted natural frequencies. As a result, only translational springs are considered.

The elastic modulus, axial load, and translational stiffness are simultaneously predicted by iterating through values of the dimensionless axial load parameter α and translational stiffness parameter $β$. To begin, an initial coarse scan was performed with discrete values in the range of 10−4 ≤ α ≤ 106 and $101≤β≤106$. These values were chosen based on the preliminary studies examined in Sec. 3.1 with an increase in the upper bound to fully capture the error profile. The logarithm of the normalized error in predicted natural frequency, defined in Eq. (22), is plotted versus α and $β$ in Fig. 11(a). For this example, the ranges of the iterated values for α and $β$ were refined three times to better capture the global minimum of the normalized error. This refinement process is pictured in Fig. 11(a), where a minimum error of ɛ = 5.2461 × 10−4 was found at values $α^=953.0056$ and $β^=4.3545×104$. With $α^$ and $β^$, the elastic modulus, axial load, and translational stiffness were predicted and are listed in Table 6.

Fig. 11
Fig. 11
Close modal

The accuracy of the proposed method and predicted beam properties may be examined by comparing the measured and predicted natural frequencies in Fig. 11(b). The predicted natural frequencies are in better agreement than those depicted in Fig. 10(b), with an error of ɛ = 5.2461 × 10−4. By modeling the homogeneous silicon beam with translational springs, the original error of ɛ = 3.1678 × 10−3 obtained from modeling the entire beam with ideal boundary conditions was reduced by an order of magnitude. Additionally, the predicted value of the elastic modulus, $E^=212.7958GPa$, still falls within the previously mentioned range of reported values. The normal tensile stress found from the predicted axial load may also be compared to the silicon wafer specifications, where $σ=83.5421MPa$ is in agreement with the upper bound of 1000 MPa. Although it is difficult to validate the predicted translational stiffness, $k^t=8.0899N/m$, the predicted mode shapes may be observed in Fig. 11(c). Here, the predicted mode shapes follow the general trend of mode shapes for a beam with fixed–fixed boundary conditions. However, it is clear that both ends do not go to zero, which is expected considering the physical structure and model in Figs. 9(c) and 9(d). To recall, the translational springs model the portion of the NEMS beam with the gold electrode, which vibrate at each mode. Consequently, the predicted mode shapes support the expectations that energy would be distributed between bending the beam and extending the boundary springs.

The two sets of predicted beam properties listed in Table 6 may also be compared with one another. As the predicted beam properties were found using two different beam models with different boundary conditions, it is expected that the predicted beam properties will be different. This is observed in Table 6, where the inclusion of translational springs resulted in an increase in elastic modulus and decrease in axial load. From Eq. (8), these changes may be thought of as a change in the distribution of stiffness contributions. By modeling the boundary conditions with translational springs, the boundaries effectively became less stif,f which then resulted in the observed increase in elastic modulus and decrease in axial load. When analyzing such results, the main thing to note is that for each boundary condition, the predicted beam properties minimized the error between predicted and measured natural frequencies, as observed in Figs. 10(a) and 11(a). Hence, the proposed method predicts beam properties that most accurately predict natural frequencies for the model under consideration. With the predicted beam properties, the beam modeled with ideal fixed–fixed boundary conditions offers a good representation of the NEMS beam with a low error of ɛ = 3.1678 × 10−3 in predicted natural frequencies. This error may be further reduced to ɛ = 5.2461 × 10−4 by modeling the beam with flexible boundary conditions and the corresponding predicted beam properties.

### 6.3 Efficiency of Proposed Method.

The proposed method offers major computational savings as it simultaneously predicts the elastic modulus from the determined scaling factor, thus reducing the dimension of the search space. Such reduction has profound advantages when performing brute force scans, such as those performed in the previous experimental examples. For example, in Sec. 6.1, the elastic modulus and axial load were predicted by scanning over 500 discrete values of the dimensionless axial load parameter α. The total central processing unit (CPU) time to perform this search was measured to be 18.0400 s. This time would significantly increase if the proposed method is not employed. Considering this same example, a standard brute force scan that scans over discrete values of both the elastic modulus and axial load would increase the one-dimensional search space to a two-dimensional search space. The required CPU time would then increase proportionally with the number of values evaluated for the elastic modulus. For example, a standard search that evaluates 500 discrete values of the elastic modulus and 500 discrete values of the axial load would require approximately 500 times more computing time than the proposed method. This would approximately result in a search that requires 2.5 h of CPU time. The reduction in search space becomes even more beneficial when considering a beam with flexible boundary conditions, where the elastic modulus, axial load, and translational spring constant must be predicted. With the proposed method employed in Sec. 6.2, 151 values of the dimensionless axial load parameter, α, and 151 values of the dimensionless translational stiffness parameter $β$ were scanned. This initial coarse scan required 741.1000 s of CPU time. If the proposed method was not employed, a three-dimensional search would have been required. For this case, if 151 values of the elastic modulus were also scanned, the required CPU time would significantly increase to approximately 31 h. As a result, the proposed method offers major computational savings due to its reduction in the dimension of the search space. Results presented are from tests performed on a machine with a 2.9 GHz dual-core Intel Core i5 processor (Intel, Santa Clara, CA).

### 6.4 Number of Measured Frequencies.

Finally, it is valuable to examine the effect the number of measured natural frequencies has on the predicted beam properties. For this study, the ideal boundary conditions described in Sec. 6.1 are used to model the NEMS beam, and the elastic modulus and axial load are simultaneously predicted with the proposed method. In this example, the beam properties are predicted for an increasing number of measured natural frequencies. First, the proposed method is employed while only using the first two natural frequencies to compute the normalized error in Eq. (22). Additional measured natural frequencies were then continuously included in the analysis until the beam properties were predicted using all measured frequencies from the first 11 modes. The predicted elastic modulus $E^$ and dimensionless axial load parameter $α^$ from the various simulations are plotted in Fig. 12(a). The results indicate that the value of the predicted beam parameters vary as more measured frequencies are included in the fitting procedure until nine or more measured natural frequencies are fit. These results suggest that measuring more natural frequencies is advantageous for accurately predicting beam properties.

Fig. 12
Fig. 12
Close modal

Additional results may be observed in Fig. 12(b), which plots the normalized error in predicted natural frequencies. The error is computed for each simulation, where the error over the fit modes and all modes are plotted. From the plot, when a small number of measured natural frequencies is used to perform the fit, the normalized error over the fit modes is small while the error over all 11 modes is relatively large. From Fig. 12(b), the error over all modes first increases when the number of fit modes increases from two to three. This initial increase is due to the relatively low error obtained for the two frequency fit. For this case, the error over all modes is dominated by the good agreement of the first two fit modes and its extremely low error of 1.6813 × 10−5. However, after this initial increase, as the number of frequencies included in the fit increases, the error over all modes decreases and the predicted beam properties become more accurate. From Fig. 12(b), it is observed that the error over all 11 modes is largest when only three measured frequencies are used in the fit, with an error of ɛ = 2.7820 × 10−2. Although this error is reasonably small, it may be reduced by including more natural frequencies in the fit, where it is observed that a minimum error of ɛ = 3.1678 × 10−3 is obtained when nine or more measured natural frequencies are used. As a result, although the proposed method only requires two measured natural frequencies, a better model and prediction of the beam is made when more measured natural frequencies are provided.

## 7 Conclusion

A method that simultaneously predicts the elastic modulus, axial load, and boundary conditions of a beam from at least two measured natural frequencies is presented. The method utilizes a FEM model of an Euler–Bernoulli beam under axial loading to evaluate a generalized eigenvalue problem that is expressed in terms of dimensionless beam parameters and frequencies. The beam properties are predicted by iterating through the dimensionless beam parameters in search for values that minimize the normalized error between predicted and measured natural frequencies. The proposed method is advantageous as it only requires two measured natural frequencies and simultaneously predicts the elastic modulus, which reduces the number of search parameters. The effect of axial loading and flexible boundary conditions on the natural frequencies of the beam was investigated. It was found that the response of the beam has a relative weak dependence on the rotational boundary stiffness. Numerical examples were presented, which validated the method’s accuracy in predicting the elastic modulus, axial load, and translational stiffness of a beam. However, relatively large errors were found when predicting rotational stiffnesses due to the beam’s insensitivity to such property. Results from an experimental example with a NEMS beam highlighted the methods ability to predict beam properties, which reduced the normalized error in predicted natural frequencies to ɛ = 3.1678 × 10−3 when the beam was modeled with ideal boundary conditions. The normalized error was further reduced to ɛ = 5.2461 × 10−4 when the beam was modeled with translational springs. Since the analysis and FEM are formulated from dimensionless parameters, the proposed method may be used for beams with any length scale while maintaining numerical stability.

## Acknowledgment

This material is based on the work supported by the National Science Foundation Graduate Research Fellowship Program (Grant No. DGE-1840990). Partial support from the National Science Foundation (Grant No. CBET-1604075) is also acknowledged. Any opinions, findings, and conclusions or recommendations expressed in this material are those of the authors and do not necessarily reflect the views of the National Science Foundation. The authors thank Professor M. Selim Hanay of Bilkent University, Turkey, for help with NEMS device fabrication.

## References

1.
Eom
,
K.
,
Park
,
H. S.
,
Yoon
,
D. S.
, and
Kwon
,
T.
,
2011
, “
Nanomechanical Resonators and Their Applications in Biological/Chemical Detection: Nanomechanics Principles
,”
Phys. Rep.
,
503
(
4–5
), pp.
115
163
. 10.1016/j.physrep.2011.03.002
2.
Ekinci
,
K. L.
, and
Roukes
,
M. L.
,
2005
, “
Nanoelectromechanical Systems
,”
Rev. Sci. Instrum.
,
76
(
6
), p.
061101
. 10.1063/1.1927327
3.
Ekinci
,
K. L.
,
Huang
,
X. M. H.
, and
Roukes
,
M. L.
,
2004
, “
Ultrasensitive Nanoelectromechanical Mass Detection
,”
Appl. Phys. Lett.
,
84
(
22
), pp.
4469
4471
. 10.1063/1.1755417
4.
Yang
,
Y. T.
,
Callegari
,
C.
,
Feng
,
X. L.
,
Ekinci
,
K. L.
, and
Roukes
,
M. L.
,
2006
, “
Zeptogram-Scale Nanomechanical Mass Sensing
,”
Nano Lett.
,
6
(
4
), pp.
583
586
. 10.1021/nl052134m
5.
Jensen
,
K.
,
Kim
,
K.
, and
Zettl
,
A.
,
2008
, “
An Atomic-Resolution Nanomechanical Mass Sensor
,”
Nat. Nanotechnol.
,
3
(
9
), p.
533
. 10.1038/nnano.2008.200
6.
Naik
,
A. K.
,
Hanay
,
M. S.
,
Hiebert
,
W. K.
,
Feng
,
X. L.
, and
Roukes
,
M. L.
,
2009
, “
Towards Single-Molecule Nanomechanical Mass Spectrometry
,”
Nat. Nanotechnol.
,
4
(
7
), p.
445
. 10.1038/nnano.2009.152
7.
Hanay
,
M. S.
,
Kelber
,
S.
,
Naik
,
A. K.
,
Chi
,
D.
,
Hentz
,
S.
,
Bullard
,
E. C.
,
Colinet
,
E.
,
Duraffourg
,
L.
, and
Roukes
,
M. L.
,
2012
, “
Single-Protein Nanomechanical Mass Spectrometry in Real Time
,”
Nat. Nanotechnol.
,
7
(
9
), p.
602
. 10.1038/nnano.2012.119
8.
Hanay
,
M. S.
,
Kelber
,
S. I.
,
O’Connell
,
C. D.
,
Mulvaney
,
P.
,
,
J. E.
, and
Roukes
,
M. L.
,
2015
, “
Inertial Imaging With Nanomechanical Systems
,”
Nat. Nanotechnol.
,
10
(
4
), p.
339
. 10.1038/nnano.2015.32
9.
Yuksel
,
M.
,
Orhan
,
E.
,
Yanik
,
C.
,
Ari
,
A. B.
,
Demir
,
A.
, and
Hanay
,
M. S.
,
2019
, “
Nonlinear Nanomechanical Mass Spectrometry at the Single-Nanoparticle Level
,”
Nano Lett.
,
19
(
6
), pp.
3583
3589
. 10.1021/acs.nanolett.9b00546
10.
Cleland
,
A. N.
, and
Roukes
,
M. L.
,
1998
, “
A Nanometre-Scale Mechanical Electrometer
,”
Nature
,
392
(
6672
), pp.
160
162
. 10.1038/32373
11.
Kara
,
V.
,
Yakhot
,
V.
, and
Ekinci
,
K. L.
,
2017
, “
Generalized Knudsen Number for Unsteady Fluid Flow
,”
Phys. Rev. Lett.
,
118
(
7
), p.
074505
. 10.1103/PhysRevLett.118.074505
12.
Villa
,
M. M.
, and
Paul
,
M. R.
,
2009
, “
Stochastic Dynamics of Micron-Scale Doubly Clamped Beams in a Viscous Fluid
,”
Phys. Rev. E Stat. Nonlinear Soft Matter Phys.
,
79
(
5 Pt 2
), p.
056314
. 10.1103/PhysRevE.79.056314
13.
Kara
,
V.
,
Sohn
,
Y. I.
,
Atikian
,
H.
,
Yakhot
,
V.
,
Loncar
,
M.
, and
Ekinci
,
K. L.
,
2015
, “
Nanofluidics of Single-Crystal Diamond Nanomechanical Resonators
,”
Nano Lett.
,
15
(
12
), pp.
8070
8076
. 10.1021/acs.nanolett.5b03503
14.
Tamayo
,
J.
,
Kosaka
,
P. M.
,
Ruz
,
J. J.
,
San Paulo
,
A.
, and
Calleja
,
M.
,
2013
, “
Biosensors Based on Nanomechanical Systems
,”
Chem. Soc. Rev.
,
42
(
3
), pp.
1287
1311
. 10.1039/C2CS35293A
15.
Treml
,
R.
,
Kozic
,
D.
,
Zechner
,
J.
,
Maeder
,
X.
,
Sartory
,
B.
,
Gänser
,
H. P.
,
Schöngrundner
,
R.
,
Michler
,
J.
,
Brunner
,
R.
, and
Kiener
,
D.
,
2016
, “
High Resolution Determination of Local Residual Stress Gradients in Single- and Multilayer Thin Film Systems
,”
Acta Mater.
,
103
, pp.
616
623
. 10.1016/j.actamat.2015.10.044
16.
Wang
,
S.
,
Chen
,
J.
,
Li
,
D.
,
Huang
,
Y.
,
Li
,
Z.
, and
Zhang
,
W.
,
2006
, “
Evaluating Interface Effect on Stresses in Thin Films by a Local Curvature Metrology With High Accuracy and Resolution
,”
2006 1st IEEE International Conference on Nano/Micro Engineered and Molecular Systems
,
Zhuhai, China
,
Jan. 18–21
, pp.
1513
1516
.
17.
Fang
,
W.
, and
Lo
,
C.-Y.
,
2000
, “
On the Thermal Expansion Coefficients of Thin Films
,”
Sens. Actuators., A
,
84
(
3
), pp.
310
314
. 10.1016/S0924-4247(00)00311-3
18.
Kim
,
M.-G.
,
Kanatzidis
,
M. G.
,
Facchetti
,
A.
, and
Marks
,
T. J.
,
2011
, “
Low-Temperature Fabrication of High-Performance Metal Oxide Thin-Film Electronics Via Combustion Processing
,”
Nat. Mater.
,
10
(
5
), pp.
382
388
. 10.1038/nmat3011
19.
Li
,
S.
,
Reynders
,
E.
,
Maes
,
K.
, and
De Roeck
,
G.
,
2013
, “
Vibration-Based Estimation of Axial Force for a Beam Member With Uncertain Boundary Conditions
,”
J. Sound Vib.
,
332
(
4
), pp.
795
806
. 10.1016/j.jsv.2012.10.019
20.
Rebecchi
,
G.
,
Tullini
,
N.
, and
Laudiero
,
F.
,
2013
, “
Estimate of the Axial Force in Slender Beams With Unknown Boundary Conditions Using One Flexural Mode Shape
,”
J. Sound Vib.
,
332
(
18
), pp.
4122
4135
. 10.1016/j.jsv.2013.03.018
21.
Maes
,
K.
,
Peeters
,
J.
,
Reynders
,
E.
,
Lombaert
,
G.
, and
De Roeck
,
G.
,
2013
, “
Identification of Axial Forces in Beam Members by Local Vibration Measurements
,”
J. Sound Vib.
,
332
(
21
), pp.
5417
5432
. 10.1016/j.jsv.2013.05.017
22.
Greening
,
P. D.
, and
Lieven
,
N. A. J.
,
2003
, “
,”
J. Sound Vib.
,
260
(
1
), pp.
101
115
. 10.1016/S0022-460X(02)00902-1
23.
Park
,
S.
,
Choi
,
S.
,
Oh
,
S.-T.
,
Stubbs
,
N.
, and
Song
,
H.-C.
,
2006
, “
Identification of the Tensile Force in High-Tension Bars Using Modal Sensitivities
,”
Int. J. Solids Struct.
,
43
(
10
), pp.
3185
3196
. 10.1016/j.ijsolstr.2005.06.089
24.
Kim
,
B. H.
, and
Park
,
T.
,
2007
, “
Estimation of Cable Tension Force Using the Frequency-Based System Identification Method
,”
J. Sound Vib.
,
304
(
3–5
), pp.
660
676
. 10.1016/j.jsv.2007.03.012
25.
Bahra
,
A. S.
, and
Greening
,
P. D.
,
2011
, “
Identifying Multiple Axial Load Patterns Using Measured Vibration Data
,”
J. Sound Vib.
,
330
(
15
), pp.
3591
3605
. 10.1016/j.jsv.2011.02.024
26.
Stachiv
,
I.
,
Kuo
,
C.-Y.
,
Fang
,
T.-H.
, and
Mortet
,
V.
,
2016
, “
Simultaneous Determination of the Residual Stress, Elastic Modulus, Density and Thickness of Ultrathin Film Utilizing Vibrating Doubly Clamped Micro-/Nanobeams
,”
,
6
(
4
), p.
045005
. 10.1063/1.4947031
27.
Pratap
,
R.
, and
Behera
,
A. R.
,
2018
, “
Simultaneous Determination of Young’s Modulus and Residual Stress in PECVD a-SiC From Postbuckling Vibration of MEMS Beams
,”
ECS Trans.
,
86
(
16
), pp.
87
100
. 10.1149/08616.0087ecst
28.
Logan
,
D. L.
,
2011
,
A First Course in the Finite Element Method
, 5th ed.,
CL Engineering
,
Pacific Grove, CA
.
29.
Sa¸kar
,
G.
,
2013
, “
The Effect of Axial Force on the Free Vibration of an Euler-Bernoulli Beam Carrying a Number of Various Concentrated Elements
,”
Shock Vib.
,
20
(
3
), pp.
357
367
. 10.1155/2013/735061
30.
Rao
,
S. S.
,
2016
,
Mechanical Vibrations
, 6th ed.,
Pearson
,
.
31.
Bokaian
,
A.
,
1988
, “
Natural Frequencies of Beams Under Compressive Axial Loads
,”
J. Sound Vib.
,
126
(
1
), pp.
49
65
. 10.1016/0022-460X(88)90397-5
32.
Bokaian
,
A.
,
1990
, “
Natural Frequenices of Beams Under Tensile Axial Loads
,”
J. Sound Vib.
,
142
(
3
), pp.
481
498
. 10.1016/0022-460X(90)90663-K
33.
Wu
,
B.
,
Xu
,
Z.
, and
Li
,
Z.
,
2007
, “
A Note on Imposing Displacement Boundary Conditions in Finite Element Analysis
,”
Commun. Numer. Methods Eng.
,
24
(
9
), pp.
777
784
. 10.1002/cnm.989
34.
Wei
,
X.
,
Chen
,
Q.
,
Xu
,
S.
,
Peng
,
L.
, and
Zuo
,
J.
,
2009
, “
Beam to String Transition of Vibrating Carbon Nanotubes Under Axial Tension
,”
,
19
(
11
), pp.
1753
1758
35.
Petersen
,
K. E.
,
1982
, “
Silicon As a Mechanical Material
,”
Proc. IEEE
,
70
(
5
), pp.
420
457
. 10.1109/PROC.1982.12331
36.
Marburg
,
S.
,
2002
, “
Six Boundary Elements Per Wavelength: Is That Enough?
J. Comput. Acoust.
,
10
(
01
), pp.
25
51
. 10.1142/S0218396X02001401
37.
Ugray
,
Z.
,
Lasdon
,
L.
,
Plummer
,
J.
,
Glover
,
F.
,
Kelly
,
J.
, and
Marti
,
R.
,
2007
, “
Scatter Search and Local NLP Solvers: A Multistart Framework for Global Optimization
,”
INFORMS J. Comput.
,
19
(
3
), pp.
328
340
. 10.1287/ijoc.1060.0175
38.
Byrd
,
R. H.
,
Gilbert
,
J. C.
, and
Nocedal
,
J.
,
2000
, “
A Trust Region Method Based on Interior Point Techniques for Nonlinear Programming
,”
Math. Program.
,
89
(
1
), pp.
149
185
. 10.1007/PL00011391
39.
,
J. E.
,
,
M.
, and
Friswell
,
M. I.
,
2011
, “
The Sensitivity Method in Finite Element Model Updating: A Tutorial
,”
Mech. Syst. Sig. Process.
,
25
(
7
), pp.
2275
2296
. 10.1016/j.ymssp.2010.10.012
40.
Arı
,
A. B.
,
Çağatay Karakan
,
M.
,
Yanık
,
C.
,
Kaya
,
I. I.
, and
Selim Hanay
,
M.
,
2018
, “
Intermodal Coupling as a Probe for Detecting Nanomechanical Modes
,”
Phys. Rev. Appl.
,
9
(
3
), p.
034024
. 10.1103/PhysRevApplied.9.034024
41.
Bargatin
,
I.
,
Kozinsky
,
I.
, and
Roukes
,
M. L.
,
2007
, “
Efficient Electrothermal Actuation of Multiple Modes of High-Frequency Nanoelectromechanical Resonators
,”
Appl. Phys. Lett.
,
90
(
9
), p.
093116
. 10.1063/1.2709620
42.
Maia
,
N.
,
1997
,
Theoretical and Experimental Modal Analysis
,
Research Studies Press
,
Hertfordshire, UK
.
43.
Carlotti
,
G.
,
Colpani
,
P.
,
Piccolo
,
D.
,
Santucci
,
S.
,
Senez
,
V.
,
Socino
,
G.
, and
Verdini
,
L.
,
2002
, “
Measurement of the Elastic and Viscoelastic Properties of Dielectric Films Used in Microelectronics
,”
Thin Solid Films
,
414
(
1
), pp.
99
104
. 10.1016/S0040-6090(02)00430-3
44.
Khan
,
A.
,
Philip
,
J.
, and
Hess
,
P.
,
2004
, “
Young’s Modulus of Silicon Nitride Used in Scanning Force Microscope Cantilevers
,”
J. Appl. Phys.
,
95
(
4
), pp.
1667
1672
. 10.1063/1.1638886
45.
Chuang
,
W. H.
,
Luger
,
T.
,
Fettig
,
R. K.
, and
Ghodssi
,
R.
,
2004
, “
Mechanical Property Characterization of LPCVD Silicon Nitride Thin Films at Cryogenic Temperatures
,”
J. Microelectromech. Syst.
,
13
(
5
), pp.
870
879
. 10.1109/JMEMS.2004.836815
46.
Liu
,
X.
,
Metcalf
,
T. H.
,
Wang
,
Q.
, and
,
D. M.
,
2011
, “
Elastic Properties of Several Silicon Nitride Films
,”
MRS Proceedings
, Vol.
989
,
San Fransisco, CA
.
47.
Gardeniers
,
J. G. E.
,
Tilmans
,
H. A. C.
, and
Visser
,
C. C. G.
,
1996
, “
LPCVD Silicon–Rich Silicon Nitride Films for Applications in Micromechanics, Studied With Statistical Experimental Design
,”
J. Vac. Sci. Technol. A
,
14
(
5
), pp.
2879
2892
. 10.1116/1.580239
48.
Habermehl
,
S.
,
1998
, “
Stress Relaxation in Si-Rich Silicon Nitride Thin Films
,”
J. Appl. Phys.
,
83
(
9
), pp.
4672
4677
. 10.1063/1.367253
49.
Temple-Boyer
,
P.
,
Rossi
,
C.
,
Saint-Etienne
,
E.
, and
Scheid
,
E.
,
1998
, “
Residual Stress in Low Pressure Chemical Vapor Deposition SiNx Films Deposited From Silane and Ammonia
,”
J. Vac. Sci. Technol. A
,
16
(
4
), pp.
2003
2007
. 10.1116/1.581302