An Inverse Method to Predict NEMS Beam Properties From Natural Frequencies

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⁻¹² 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 [3–9], 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 [12–14]. 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 [19–21]. 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 [22–25]. 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.

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.

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.

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.

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.

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 Ω₁ 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 α > 10⁴. 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 10⁰ ≤ α ≤ 10⁵.

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 Ω₁ 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 2k_t, 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, Ω₁ = 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$⁠.

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 Ω₁ 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 Ω₁ = 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 Ω₁ = 9.8696 to Ω₁ = 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$⁠.

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 Ω₁ 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.

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.

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.1–5.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.

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.

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.

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 k_t = 0.0171 N/m with $log10(β)=2$⁠. To better interpret these results, the logarithm of the normalized error in predicted frequency log₁₀(ɛ) 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.

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.

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.

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.

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⁻⁴ ≤ α ≤ 10⁶ 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.

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⁻³. 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 [43–46]. 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 [47–49].

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 k_t. 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⁻⁴ ≤ α ≤ 10⁶ 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⁻⁴ 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.

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⁻⁴. By modeling the homogeneous silicon beam with translational springs, the original error of ɛ = 3.1678 × 10⁻³ 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⁻³ in predicted natural frequencies. This error may be further reduced to ɛ = 5.2461 × 10⁻⁴ by modeling the beam with flexible boundary conditions and the corresponding predicted beam properties.

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).

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.

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⁻⁵. 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⁻². 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⁻³ 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.

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⁻³ when the beam was modeled with ideal boundary conditions. The normalized error was further reduced to ɛ = 5.2461 × 10⁻⁴ 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.

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.