## Abstract

In this article, we present a new contact resonance atomic force microscopy-based method utilizing a square, plate-like microsensor to accurately estimate viscoelastic sample properties. A theoretical derivation, based on Rayleigh–Ritz method and on an “unconventional” generalized eigenvalue problem, is presented and a numerical experiment is devised to verify the method. We present an updated sensitivity criterion that allows users, given a set of measured in-contact eigenfrequencies and modal damping ratios, to select the best eigenfrequency for accurate data estimation. The verification results are then presented and discussed. Results show that the proposed method performs extremely well in the identification of viscoelastic properties over broad ranges of nondimensional sample stiffness and damping values.

## 1 Introduction

Contact resonance atomic force microscopy (CR-AFM), originally pioneered by Rabe et al. [1], uses the measured freely vibrating and in-contact resonance frequencies of a coupled microsensor-sample system to estimate material properties of the sample. Contact resonance (CR) spectroscopy AFM techniques are an ideal candidate for nanoscale quantitative mechanical property extraction due to their avoidance of nonlinear interactions present in other AFM measurement techniques. CR-AFM operates primarily in the linear net-repulsive regime of the tip-sample interaction force curve. CR-AFM has been used to measure both elastic and viscoelastic sample properties [26]. Recently, Aureli et al. [7] showed that by using a square, plate-like cantilevered microsensor, as opposed to the traditional rectangular, narrow beam AFM cantilevered sensor, several advantages arise. These advantages include an increase of sensing modes available in a given frequency bandwidth and the ability to optimally place the sensing tip within the plate domain to maximize measurement sensitivity.

In this study, we expand the method presented in Ref. [7], which was formulated to test purely elastic samples, to include viscoelastic materials. With this addition, the aforementioned measurement benefits of the proposed sensor can be extended to include the measurement of many viscoelastic materials of interest, including, for example, polymer composites and biological materials. In addition, with the inclusion of damping in the system, it becomes possible to characterize and account for hydrodynamic forces present in liquid measurement environments, such as those encountered when conducting AFM measurements of live biological samples ex vivo.

To incorporate viscoelastic effects, Hamilton’s least action principle is used to create a variational formulation of the problem. The strain energy of the plate is calculated following the assumptions of the Kirchhoff–Love plate equations. The sample contact, originally modeled in Ref. [7] as a linear spring element, is replaced with a Kelvin–Voigt element. We include an additional term in the variational formulation of the method that accounts for the nonconservative work done by the new viscous damper added to the model. A Rayleigh–Ritz method is used to transform the problem into a generalized eigenvalue problem (EVP) with complex eigenvalues. A nontraditional interpretation of the eigenvalue problem leads to an efficient method for the estimation of the sample viscoelasticity. Specifically, given the measured in-contact natural frequencies and modal damping ratios of the coupled sensor-sample system, the solution of the generalized EVP yields estimates of the viscoelastic material properties, including stiffness and damping.

Following the derivation of the theoretical model, a numerical experiment, using the finite element method, is created to verify the method. Given a set of measured in-contact eigenmodes, a criterion is proposed to select which of these frequencies should be used to produce the most accurate material property estimates. In addition, we discuss methods to extract the model input data from experimentally measured frequency response functions of the system. Finally, the verification results are presented over a broad range of parameters of interest, and future research directions of the work are discussed.

## 2 Theoretical Framework

### 2.1 Governing Equations.

In this study, we extend the theoretical formalism developed in our previous paper [7], which is concerned with purely elastic response of the sample, to incorporate the effect of viscoelasticity in the sample mechanical behavior. As in Ref. [7], we will consider the free vibrations of a thin rectangular cantilevered Kirchhoff–Love plate as illustrated in Fig. 1. The plate is clamped at its edge at x = 0 and free at the edges x = Lx and y = ±Ly/2, see Fig. 1. In-plane dimensions of the plate are its length Lx and width Ly, and the plate thickness is indicated with h. We assume that the mechanical constitutive behavior of the plate is linear elastic, isotropic and homogeneous, with Young’s modulus E and Poisson’s ratio $ν$, and we denote its mass density per unit volume with $ρ$.

Fig. 1
Fig. 1
Close modal

The plate is in contact through an attachment point of coordinates (xs, ys) with a spring-damper element, which models the Kelvin–Voigt viscoelastic behavior of a sample being analyzed via contact resonance AFM. For simplicity, both the elastic and the viscous behaviors are assumed linear and fully described by a single constant each, that is, the sample stiffness ks and the sample damping coefficient cs. We restrict our study to small oscillations, so that the sample is assumed to react to displacements of the attachment point in the z-direction only.

The governing equations of the free vibrations of the plate supported by the viscoelastic sample are as follows [710]
$D∇4w(x,y;t)+ρhw¨(x,y;t)=−[ksw(x,y;t)+csw˙(x,y;t)]δ(x−xs)δ(y−ys)$
(1)
where w(x, y; t) is the plate deflection in the z-direction, t is the time variable, $D=Eh3/(12(1−ν2))$ is the plate stiffness modulus, $∇4(∙)$ is the bi-Laplacian operator, $δ(∙)$ is the Dirac delta distribution, and a superimposed dot denotes time derivative. Boundary conditions are given as follows [710]:
$w=0,∂w∂x=0atx=0$
(2a)
$∂3w∂x3+(2−ν)∂3w∂x∂y2=0,∂2w∂x2+ν∂2w∂y2=0atx=Lx$
(2b)
$∂3w∂y3+(2−ν)∂3w∂y∂x2=0,∂2w∂y2+ν∂2w∂x2=0aty=±Ly2$
(2c)
We scale the aforementioned equations by dividing all lengths by Lx, so that now the domain P of the plate is described as $P={(x~,y~)∈[0,1]×[−ℓ/2,ℓ/2]}$, where = Ly/Lx indicates the plate aspect ratio and a superimposed tilde denotes dimensionless variables. The plate deflection in the z-direction is scaled as $w/A0=w~$, where A0 is some characteristic displacement amplitude. Following Ref. [7], we introduce the parameter $Ω=D/(ρhLx4)$ that incorporates material and geometric properties of the plate. We further nondimensionalize the sample stiffness and damping by defining the usual cantilever stiffness $kc=3ℓ(1−ν2)D/Lx2$ and introducing the parameters
$α=kskc,β=csΩ/kc$
(3)
Note that kc is in effect the stiffness of an equivalent Euler–Bernoulli cantilever beam with length Lx, width Ly, thickness h, and Young’s modulus E, see also Ref. [7]. While keeping the time variable dimensional, the scaled form of Eq. (1) now reads
$∇~4w~+Ω−2w~¨=−(αw~+βΩ−1w~˙)[3ℓ(1−ν2)]δ~(x~−x~s)δ~(y~−y~s)$
(4)
where we have omitted the arguments of $w~$ for ease of notation.

### 2.2 The Rayleigh–Ritz Eigenvalue Problem.

By following our original derivation in Ref. [7], we now develop a Rayleigh–Ritz method [11,12] to estimate the complex eigenvalues of Eq. (4). We assume the following ansatz for the motion of the plate $w(x,y;t)=W(x,y)eλt$, where W(x, y) indicates the complex displacement amplitude and $λ∈C$ is a separation constant. This important difference from Ref. [7] is necessary to capture both frequency and damping ratio for the separable solutions of Eq. (4). Note, in fact, that the system studied in Ref. [7] is undamped so it was therein assumed that $λ=±iω$.

Corresponding to this ansatz, the kinetic energy of the plate is thus [7,13]
$T=12λ2ρhLx2(A0eλt)2∫PW~2dx~dy~$
(5)
where $W~=W/A0$ indicates the complex displacement amplitude W scaled with the characteristic displacement amplitude A0. Similarly, neglecting the contribution of strain energy due to membrane and shear effects [10], the potential elastic strain energy due to plate bending Ub is written as follows [7,8]:
$Ub=12D(A0eλt)2Lx2∫P[(W~,x~x~)2+(W~,y~y~)2+2νW~,x~x~W~,y~y~+2(1−ν)(W~,x~y~)2]dx~dy~$
(6)
where a comma subscript indicates differentiation with respect to the indicated variable.
The potential (elastic) energy due to the sample deformation Us is written as follows:
$Us=12ks(A0eλt)2W~(x~s,y~s)2$
(7)
and the work of the nonconservative forces due to the sample viscosity Wncs is given by
$Wncs=12csλ(A0eλt)2W~(x~s,y~s)2$
(8)
In the following, we will refer exclusively to dimensionless variables, so we will drop the superimposed tilde for ease of notation. We use the energy expressions in Eqs. (5)(8) to derive the characteristic equation for the eigenvalue problem of the damped plate vibration in a matrix form. We use Hamilton’s least action principle [11,12] with the inclusion of the work of the nonconservative forces to obtain
$∫t1t2[δT−δUb−δUs−δWncs]dt=0$
(9)
where the symbol $δ$ indicates the first variation. All the variations in Eq. (9) can be found explicitly worked out in Ref. [7], except $δWncs$ which reads $δWncs=csλ(A0eλt)2W~(x~s,y~s)δW~(x~s,y~s)$.
Now, we assume that the solution W(x, y) can be expressed in terms of a weighted series of admissible functions Φj(x, y) to be specified later, that is, $W(x,y)=∑j=1NΦj(x,y)qj$, with qj some weights. Then, $δW(x,y)=∑i=1NΦi(x,y)δqi$ since Φi(x, y) are independent and the qi are arbitrary. These representations for W(x, y) and $δW(x,y)$ are substituted into Hamilton’s principle in Eq. (9). By cancelling the nonzero factor $(A0eλt)2$ throughout and collecting the coefficients of each $δqi$ and setting them equal to zero separately, Eq. (9) thus reduces to the N × N matrix problem:
$[(DLx2)K+(ks+csλ)S′+(λ2ρhLx2)M]q=0$
(10)
The matrices in Eq. (10) have been derived in Ref. [7] and are reported below for completeness. Specifically, $K$ is the (nondimensional) symmetric bending stiffness matrix, whose (i, j)-entry Kij is obtained from Eq. (6) as follows:
$Kij=∫P[Φ,xxjΦ,xxi+Φ,yyjΦ,yyi+ν(Φ,xxjΦ,yyi+Φ,yyjΦ,xxi)+2(1−ν)(Φ,xyjΦ,xyi)]dxdy$
(11)
Similarly, $S′$ is the (nondimensional) symmetric sample viscoelasticity matrix, whose (i, j) entry Sij is given by Eq. (7) as follows:
$Sij′=Φj(xs,ys)Φi(xs,ys)$
(12)
Finally, $M$ is the (nondimensional) symmetric mass matrix, whose (i, j)-entry Mij is given by Eq. (5) as follows:
$Mij=∫PΦjΦidxdy$
(13)
In Eq. (10), the column vector $q$ contains the entries qj, that is, the frequency-dependent “modal” coefficients of w in the basis of the functions Φj. As in Ref. [7], we divide through Eq. (10) by the term $(D/Lx2)$ and redefine the sample viscoelasticity matrix as $S=[3ℓ(1−ν2)]S′$. Thus, Eq. (10) is rewritten in the final form of the eigenvalue problem that will be used in this work as follows:
$[K+(α+μβ)S+μ2M]q=0$
(14)
where $μ=λ/Ω$, and α, $β$, and Ω are defined earlier. It is remarkable that this equation is structurally very similar to the eigenvalue problem in Ref. [7], which could be rederived from Eq. (14) by letting $β=0$ (i.e., undamped system) and $μ=±iω/Ω$ (i.e., following an ansatz of purely harmonic plate vibrations). A detailed description of the properties of the system matrices can be found in Ref. [7]. Observe that, in the absence of plate-sample contact (“unsprung” case), $α=β=0$ and Eq. (14) reduces to the usual generalized eigenvalue problem for the isolated plate $[K+μ¯2M]q¯=0$, where overbars refer to “unsprung” quantities, and $μ¯=±iω¯/Ω$, with $ω¯$ being one of the natural frequencies of the unsprung plate.
If the matrices $K$, $S$, and $M$ are known, along with the scalars α, $β$, and Ω, Eq. (14) represents a quadratic eigenvalue problem [13,14], which has been studied in detail in a large body of works, for both direct and inverse problems. One of the popular techniques that proves effective in attacking this problem is, for example, the method by Frazer, Duncan, and Collar [13,15]. We observe here that the effective stiffness matrix of the sprung system is given by $[K+αS]$, the effective damping matrix is $[βS]$, and the effective mass matrix is simply $M$, with eigenvalue μ. The system in Eq. (14) is thus rewritten in terms of 2N × 2N matrices as follows:
$[−[K+αS]00M][qμq]=μ[[βS]MM0][qμq]$
(15)
with the important advantage that the so-constructed 2N × 2N system matrices are symmetric, and the resulting modes are orthogonal. The 2N eigenvalues μ extracted from this process appear in complex conjugate pairs, for which the imaginary parts $Im[μ]$ describe the natural frequencies of the associated modes and their real parts $Re[μ]$ are related to modal damping. Eigenvectors from Eq. (15) contain information on the mode shapes of the system, but the number of degrees-of-freedom is doubled in the process.

### 2.3 Estimation of the Sample Viscoelasticity.

While Eq. (14) is typically used to calculate natural frequencies and modal damping ratios as described earlier, in this work, the model developed will be used for the identification of the sample stiffness parameter α and viscosity parameter $β$ in Eq. (3). Specifically, following in spirit the practice of CR-AFM measurements in which sample properties are estimated from a set of known frequency measurements, here we are concerned with the estimation of α and $β$ when $K$, $S$, and $M$ are known, along with Ω and a set of measured $λi$, with i = 1, …MN. As in Ref. [7], Ω can, in principle, be identified from measurements of the free unsprung vibrations of the plate.

Assuming that the parameter Ω is known and a complex eigenvalue $λi$ of Eq. (14) has been determined from an experiment or a physical measurement (see also Sec. 3), the eigenvalue problem is rearranged in the alternative form
$[(K+μi2M)−γ(−S)]qi=0$
(16)
where $γ=α+μiβ$ and $qi$ is the eigenvector associated to $λi$. Since Eq. (16) is satisfied identically, $γ$ can be interpreted as the generalized eigenvalue of the problem $(K+μi2M)qi=γ(−S)qi$, where the known matrix $(K+μi2M)$ is diagonalized with respect to $−S$. This new generalized eigenvalue problem will have N solutions for $γ$ for any μi. Following the argument in Ref. [7], we claim that there exists one and only one distinct finite solution and, therefore, only one physically admissible value for $γ$. Once this unique value is determined, the sample viscoelasticity can be unambiguously identified by noticing that $Re[γ]=α+Re[μi]β$ and $Im[γ]=Im[μi]β$, from which it follows that
$β=Im[γ]Im[μi]$
(17a)
$α=Re[γ]−Re[μi]Im[μi]Im[γ]$
(17b)
Finally, from α and $β$, the values of ks and cs are determined following the definitions in Eq. (3). The unique finite (and, therefore, smallest magnitude) eigenvalue of Eq. (16) can be easily computed for any available μi with widely available eigenvalue routines, such as for example those in lapack and mathworks matlab. Thus, each measured value of μi is uniquely associated with an estimate of α and $β$, that we indicate with $α^i$ and $β^i$, respectively. As in Ref. [7], we note that the proposed identification procedure can only succeed if $μi≠μ¯i$, that is, if there is indeed a measurable change in the eigenvalues of the problem between the unsprung and the in-contact conditions.
Given a set of measured in-contact and unsprung modal frequencies, we obtain a set of estimates $α^i$ and $β^i$. We must then determine which of these estimates are “best,” in a certain sense. To accomplish this, we define the modal sensitivity $σi$ of mode i as the quantity
$σi(α,β)=d|λi|dα$
(18)
where $λi$ should be considered functions of α and $β$. This quantity describes the change of the eigenvalue (magnitude) of a particular mode i in response to a unit increase of α. We posit that the estimates $α^i$ and $β^i$ are most reliable when the sensitivity $σi(α^i,β^i)$ is highest in the neighborhood of the particular mode i. Many different sensitivity metrics could be used to determine the suitability of the estimates, and the performance of a few of them will be analyzed in Sec. 4. As utilized in the extensive numerical campaign described later, we have found the sensitivity given in Eq. (18) to perform well, in general, over a broad range of viscoelasticity parameters and to outperform several other choices of sensitivity.

### 2.4 Some Implementation Remarks.

The calculation of the $K$, $S$, and $M$ matrices in Eqs. (11)(13) is based on the selection of suitable basis functions Φi(x, y) in the representations of W(x, y). Following our previous work [7], here we will again use basis functions separable in x and y and, for simplicity, the functions Φi(x, y) will be constructed as the Cartesian product of bases in the x-direction with bases in the y-direction that satisfy the essential boundary conditions. For the former set, we select the fixed-free eigenfunctions of an Euler–Bernoulli beam in the interval x ∈ [0, 1], that is [13],
$Xix(x)=[sin(βixx)−sinh(βixx)]−sin(βix)+sinh(βix)cos(βix)+cosh(βix)×[cos(βixx)−cosh(βixx)]$
(19)
with ix = 1, …, Nx, where Nx is the maximum number of modes retained in the x-direction. Here, $βix$ are the solutions of the characteristic equation $cosβcoshβ=−1$, whose first few values are {1.875, 4.694, 7.855, 10.996, 14.137}. For the latter set, we use a “re-scaled” version of the traditional Legendre polynomials [16] Pi, normalized as follows:
$Yiy(y)=2(iy−1)+1ℓP(iy−1)(2y/ℓ)$
(20)
for iy = 1, …, Ny, where Ny is the maximum number of modes retained in the y-direction. These polynomials represent an orthonormal and complete basis set in y ∈ [−/2, /2] with weight +1. Note that for iy = 1, 2, the terms $Y1(y)=1/ℓ$ and $Y2(y)=23/ℓ3y$ describe the even and odd rigid body modes of the cross section, respectively.

Construction of the Φi(x, y) functions as the Cartesian product of “modes” $Xix(x)$ and $Yiy(y)$, their numbering scheme, and further properties are discussed in Ref. [7]. The integrals in Eqs. (11) and (13) are approximated by using Gauss–Legendre quadrature points [17]. The complete details of the implementation are extensively discussed in Ref. [7]. We only desire to mention here that the proposed method results in a reduced-order model for the plate system, including only a relatively small number of modes (with much smaller matrix sizes than a comparably accurate finite element model) capable of accurate prediction of the dynamic characteristics of the system. As an example, when 20 modes are extracted from an N = (15 × 15) size model, calculation run-time is less than 2 s on a typical desktop personal computer, with approximately 75% of the time required to calculate the system matrices.

## 3 Numerical Experiments

### 3.1 Finite Element Analysis Simulations.

In lieu of CR-AFM experimental measurements, we thoroughly test our method for sample viscoelasticity identification via numerical experiments. Our numerical campaign is centered around a square plate-like cantilevered microscale sensor of dimensions 225 × 225 × 3 μm3, similar to one of the prototype plate sensors considered in Ref. [7]. The plate material is silicon, with Young’s modulus 169 GPa, Poisson’s ratio $ν=0.25$, and density $ρ=2320kg/m3$. The plate is connected with a linear viscoelastic element located at (xs, ys) = (0.9Lx, 0.4Ly). This location is close to, but not coincident with, the optimal location determined in Ref. [7] for maximum sensitivity for a square plate. While this location is chosen for its favorable sensitivity properties, the choice has to be regarded as purely illustrative, as the application of the method would be unchanged if a different attachment point were to be selected. The viscoelastic element can be assigned directly a stiffness ks and a viscous damping coefficient cs. A model for finite element analysis (FEA) is created for the system using the commercial software package ansys mechanical apdl v.17.0. Shell elements “SHELL63” are utilized to create the structure with a mapped quadrilateral mesh. The element size is selected as Lx/100 to ensure accuracy in the numerical results. Membrane effects are not considered. One combination element “COMBIN14” is used to model the viscoelastic sample. A comprehensive convergence study for the model was already presented in our previous work [7].

For each analyzed case, a known stiffness and damping are assigned to the linear viscoelastic element modeling, and the resulting “in-contact” natural frequencies and modal damping ratios of the system are then calculated via the software “MODAL” analysis solution routine. Because of the explicit presence of the damping in the viscoelastic element, the “DAMP” (damped system) solution method is selected. The output of this analysis consists of a set of complex eigenvalues, which essentially play the role of the $λi$ in Eq. (14). Thus, the output from the FEA is used in our algorithm described in Sec. 2 for estimating α and $β$. After the modal analysis is completed, we use the FEA model to conduct the identification of the system frequency response function (FRF) via the software “HARMIC" analysis solution routine. Specifically, we study the driving point FRF by applying a harmonic unit force F at the attachment point (xs, ys) and determining the harmonic displacement response p at the same point (xs, ys) over the range of frequencies f ∈ (0, fmax] of interest. This choice is representative and, certainly, not the only possible one. Further study is necessary to understand the potential of nondriving point FRFs in the context of viscoelasticity identification.

For the numerical study, we choose a measurement bandwidth of fmax = 2 MHz to represent frequency limitations present in real experimental systems, similar to Ref. [7]. We use 1000 equally spaced sampling points in the band 0–2 MHz to obtain a fine representation of the system FRF with complex displacement data every 2 kHz. Then, in the neighborhood of each of the resonance peaks identified from the FRF (and coincident with the output of the previous modal analysis), we conduct a second scan with another 100 frequency points to obtain a very fine FRF of the system in such a way as to obtain at least ten sample points in the 3dB bandwidth of each resonance peak, especially for lightly damped modes.

The entire analysis is automated, by calling the ansys solver on a prepared input file from within a matlab script. The script updates the analysis parameter with the current α and $β$ coefficients to be studied in the input file with which the ansys solver is called. The FEA output, consisting of modal data and FRF data, is then redirected to text files that are later re-imported within matlab for further postprocessing. Specifically, we consider 24 different combinations of parameters α (varying over five orders of magnitude to study very soft to very hard samples, when compared to the stiffness kc) and $β$ (varying over 4 orders of magnitude per each α to study very lightly to moderately damped systems). In particular, the selection of α and $β$ values is performed as follows. The six chosen values of α are 10−2, 10−1, 100, 101, 102, and 103 to explore a wide range of sample stiffness. For each value of α, physically meaningful values of $β$ are selected according to the following criterion:
$β=αΩ2πfmaxtanδ$
(21)
so that at the largest frequency of interest fmax, the ratio between the loss and storage modulus $tanδ$, see, for example, Ref. [18], of the sample is within physical limits practically observed in samples for CR-AFM. We select the following values for $tanδ$: 0.01, 0.1, 1.0, and 10 to explore very lightly to moderately damped systems. Equation (21) can also be rearranged as $cs(2πfmax)/ks=tanδ$, which further clarifies its physical meaning in terms of the viscoelastic properties of the sample.

Note that, in our previous analysis in Ref. [7], we have shown how the parameter Ω can be directly estimated from the eigenvalues of the unsprung, freely vibrating system from the FEA simulations. This procedure is analogous to the method used in traditional CR studies and, rather than trying to identify the individual geometric and material parameters separately, we choose to robustly identify the important governing parameter Ω. Therefore, in the remainder of the article, the value of Ω ≈ 1.5079 × 105 s−1, calculated from nominal geometric and material parameters of the plate system, will be presumed known.

### 3.2 Eigenvalue Extraction From the Frequency Response Function.

In a real experimental measurement, obtained, for example, by contact resonance atomic force microscopy, the frequency response function of the system is the primary measurand. In the proposed method, we require the complex eigenvalues of the in-contact system to generate material property estimates. Here, we show that the eigenvalues of the system can be effectively extracted from measured FRFs using a single-mode approximation, provided that the modes are not too closely spaced. Although, in the following, the procedure will be illustrated on a numeric dataset, it could be implemented seamlessly on actual experimental data.

As stated earlier, the free response of the plate is given by $w(x,y;t)=W(x,y)eλt$. Assuming that the modes of the system response are well spaced, we approximate each mode as a single-degree-of-freedom (SDOF) oscillator, whose free response is given by $p(t)=C+e(−ζωn+iωn1−ζ2)t$ + $C−e(−ζωn−iωn1−ζ2)t$, where C± are constants determined by the initial conditions, $ωn$ is the SDOF natural frequency, and $ζ$ is the damping ratio [13].

At any given location on the plate, the displacement frequency response to a harmonic point load F(t) at the same location (driving point FRF) is proportional to the SDOF response given in the frequency domain by
$|p^F^|=G(1−u2)2+(2ζu)2+N$
(22)
where $u=2πf/ωn$, G is an unknown gain, N is the noise floor of the measurement, and superimposed hats denote Fourier transformed variables. Given a measured FRF, the quantities $ωn$ and $ζ$ can be determined for each mode by fitting the measured data to Eq. (22). Then, for each mode, we equate the complex eigenvalue $λ$ to the SDOF assumed form. That is, $Re[λ]=−ζωn$ and $Im[λ]=±ωn1−ζ2$. Figure 2 shows extraction results for the representative case, where α = 100 and $β=12$. The eigenvalues obtained directly from the FEA modal analysis are identical to the eigenvalues extracted from the FRF generated in the FEA.
Fig. 2
Fig. 2
Close modal

Table 1 presents the extraction results for several combinations of assigned α and $β$. In general, our results demonstrate that, as long as the mode is adequately sampled (>10 sample points in the 3 dB bandwidth), the eigenvalues can be reliably calculated. Having shown the applicability of the extraction of eigenvalues from the system FRF, in the rest of this study, we will directly use the eigenvalues calculated through the FEA modal analysis routine as an input for our estimation procedure.

Table 1

Eigenvalues extracted from FRFs for several combinations of α and $β$ for various eigenmodes and comparison with results from FEA modal analysis

FRF fitFEA value% Difference
α$β$Mode$Im[λ]$$−Re[λ]$$Im[λ]$$−Re[λ]$$Im[λ]$$Re[λ]$
0.011.20 × 10−0340.657571.36 × 10−040.657571.36 × 10−040.00−0.06
0.11.20 × 10−0371.47432.10 × 10−051.47432.11 × 10−050.000.42
0.11.20 × 10−0230.516821.84 × 10−040.516821.83 × 10−040.00−0.21
11.20 × 10−0320.228452.21 × 10−040.228452.21 × 10−040.000.00
11.20 × 10−0210.103897.32 × 10−040.103867.36 × 10−040.030.47
101.20 × 10−0360.804412.68 × 10−040.804412.68 × 10−040.000.00
101.20 × 10−0271.33438.03 × 10−041.33438.04 × 10−040.000.12
1001.20 × 10−0251.17236.01 × 10−041.17236.01 × 10−040.00−0.02
1001.20 × 10+0110.142194.04 × 10−040.142194.04 × 10−040.000.00
10001.20 × 10+0051.24133.16 × 10−041.24133.16 × 10−040.000.01
FRF fitFEA value% Difference
α$β$Mode$Im[λ]$$−Re[λ]$$Im[λ]$$−Re[λ]$$Im[λ]$$Re[λ]$
0.011.20 × 10−0340.657571.36 × 10−040.657571.36 × 10−040.00−0.06
0.11.20 × 10−0371.47432.10 × 10−051.47432.11 × 10−050.000.42
0.11.20 × 10−0230.516821.84 × 10−040.516821.83 × 10−040.00−0.21
11.20 × 10−0320.228452.21 × 10−040.228452.21 × 10−040.000.00
11.20 × 10−0210.103897.32 × 10−040.103867.36 × 10−040.030.47
101.20 × 10−0360.804412.68 × 10−040.804412.68 × 10−040.000.00
101.20 × 10−0271.33438.03 × 10−041.33438.04 × 10−040.000.12
1001.20 × 10−0251.17236.01 × 10−041.17236.01 × 10−040.00−0.02
1001.20 × 10+0110.142194.04 × 10−040.142194.04 × 10−040.000.00
10001.20 × 10+0051.24133.16 × 10−041.24133.16 × 10−040.000.01

## 4 Results and Discussion

The results of the estimations, using the FEA calculated eigenvalues and Eqs. (16) and (17), are shown in Table 2, which reports the percent difference between the estimated values and assigned values in the two-dimensional (α, $β$) space explored.

Table 2

Predicted $α^$ and $β^$ versus assigned values α and $β$, with percent differences (%D)

α$β$Mode$α^$$β^$%D α%D $β$
0.011.20 × 10−0619.71 × 10−031.20 × 10−06−2.90.0
0.011.20 × 10−0519.71 × 10−031.20 × 10−05−2.90.0
0.011.20 × 10−0419.71 × 10−031.20 × 10−04−2.90.0
0.011.20 × 10−0319.71 × 10−031.20 × 10−03−2.90.0
0.11.20 × 10−0519.97 × 10−021.20 × 10−05−0.30.0
0.11.20 × 10−0419.97 × 10−021.20 × 10−04−0.30.0
0.11.20 × 10−0319.97 × 10−021.20 × 10−03−0.30.0
0.11.20 × 10−0219.97 × 10−021.20 × 10−02−0.30.0
11.20 × 10−0429.98 × 10−011.20 × 10−04−0.20.0
11.20 × 10−0329.98 × 10−011.20 × 10−03−0.20.0
11.20 × 10−0229.98 × 10−011.20 × 10−02−0.20.0
11.20 × 10−0129.99 × 10−011.20 × 10−01−0.10.0
101.20 × 10−0329.98 × 10+001.20 × 10−03−0.2−0.4
101.20 × 10−0229.98 × 10+001.20 × 10−02−0.2−0.4
101.20 × 10−0129.98 × 10+001.20 × 10−01−0.2−0.4
101.20 × 10+0041.01 × 10+011.21 × 10+000.90.5
1001.20 × 10−0259.87 × 10+011.17 × 10−02−1.3−2.6
1001.20 × 10−0159.87 × 10+011.17 × 10−01−1.3−2.6
1001.20 × 10+0059.93 × 10+011.17 × 10+00−0.7−2.5
1001.20 × 10+0119.86 × 10+011.13 × 10+01−1.4−5.5
10001.20 × 10−0128.47 × 10+028.62 × 10−02−15.3−28.1
10001.20 × 10+0028.47 × 10+028.61 × 10−01−15.3−28.2
10001.20 × 10+0128.55 × 10+028.64 × 10+00−14.5−28.0
10001.20 × 10+0221.40 × 10+037.61 × 10+0140.0−36.6
α$β$Mode$α^$$β^$%D α%D $β$
0.011.20 × 10−0619.71 × 10−031.20 × 10−06−2.90.0
0.011.20 × 10−0519.71 × 10−031.20 × 10−05−2.90.0
0.011.20 × 10−0419.71 × 10−031.20 × 10−04−2.90.0
0.011.20 × 10−0319.71 × 10−031.20 × 10−03−2.90.0
0.11.20 × 10−0519.97 × 10−021.20 × 10−05−0.30.0
0.11.20 × 10−0419.97 × 10−021.20 × 10−04−0.30.0
0.11.20 × 10−0319.97 × 10−021.20 × 10−03−0.30.0
0.11.20 × 10−0219.97 × 10−021.20 × 10−02−0.30.0
11.20 × 10−0429.98 × 10−011.20 × 10−04−0.20.0
11.20 × 10−0329.98 × 10−011.20 × 10−03−0.20.0
11.20 × 10−0229.98 × 10−011.20 × 10−02−0.20.0
11.20 × 10−0129.99 × 10−011.20 × 10−01−0.10.0
101.20 × 10−0329.98 × 10+001.20 × 10−03−0.2−0.4
101.20 × 10−0229.98 × 10+001.20 × 10−02−0.2−0.4
101.20 × 10−0129.98 × 10+001.20 × 10−01−0.2−0.4
101.20 × 10+0041.01 × 10+011.21 × 10+000.90.5
1001.20 × 10−0259.87 × 10+011.17 × 10−02−1.3−2.6
1001.20 × 10−0159.87 × 10+011.17 × 10−01−1.3−2.6
1001.20 × 10+0059.93 × 10+011.17 × 10+00−0.7−2.5
1001.20 × 10+0119.86 × 10+011.13 × 10+01−1.4−5.5
10001.20 × 10−0128.47 × 10+028.62 × 10−02−15.3−28.1
10001.20 × 10+0028.47 × 10+028.61 × 10−01−15.3−28.2
10001.20 × 10+0128.55 × 10+028.64 × 10+00−14.5−28.0
10001.20 × 10+0221.40 × 10+037.61 × 10+0140.0−36.6

Note: The eigenmode used for the prediction, selected via the maximum sensitivity criterion $σi(α,β)=d|λi|/dα$ in Eq. (18), is indicated.

Table 2 also indicates the specific mode used for the estimation, selected by the maximum sensitivity criterion given in Eq. (18). We see that the error of the estimation is below $3%$ for values of α ≤ 100 for all combinations of $β$. For α ≥ 1000, the error in the estimation begins to increase, becoming successively larger with increasing $β$. Interestingly, the mode number used does not increase monotonically when using the sensitivity criterion in Eq. (18). For example, for α = 100, $β=12$, the most sensitive mode is mode 1, while for α = 1000, the criterion in Eq. (18) selects mode 2.

The selection of the form in Eq. (18) for the sensitivity is inspired by the results in Ref. [7], in which, however, the eigenvalues are purely imaginary numbers since the system is undamped. Thus, it is of interest to evaluate whether other sensitivity parameters, leveraging the complex nature of the present eigenvalues, may offer a better estimation performance. To this aim, Table 3 presents the prediction results for α and $β$ using the mode with maximum sensitivity using $σi(α,β)=d|λi|/dβ$. We see that by using this sensitivity criterion, the predictions are less accurate compared with those using Eq. (18) in nearly all cases. In particular, the α estimates for assigned values of α ≤ 0.1 are extremely poor. On the other hand, the α estimates for assigned values of 1 ≤ α ≤ 100 are comparable to estimates using Eq. (18), except for the α = 100, $β=12$ case. In general, $β$ estimates perform well for assigned values of α ≤ 100. For large assigned values of α, the predictions for both α and $β$ diverge. In addition, and very interestingly, the mode selected by this criterion is much different than those selected by Eq. (18), possibly because of how the natural frequency of higher modes enters in the proposed definition of sensitivity.

Table 3

Predicted $α^$ and $β^$ versus assigned values $α$ and $β$, with percent differences (%D), using the eigenmode (indicated) with maximum sensitivity $σi(α,β)=d|λi|/dβ$

α$β$Mode$α^$$β^$%D α%D $β$
0.011.20 × 10−064−2.08 × 10−011.20 × 10−06$−2182.8$0.2
0.011.20 × 10−054−2.08 × 10−011.20 × 10−05$−2182.8$0.2
0.011.20 × 10−044−2.08 × 10−011.20 × 10−04$−2182.8$0.2
0.011.20 × 10−034−2.08 × 10−011.20 × 10−03$−2182.8$0.2
0.11.20 × 10−054−1.18 × 10−011.20 × 10−05$−218.4$0.2
0.11.20 × 10−044−1.18 × 10−011.20 × 10−04$−218.4$0.2
0.11.20 × 10−034−1.18 × 10−011.20 × 10−03$−218.4$0.2
0.11.20 × 10−024−1.18 × 10−011.20 × 10−02$−218.0$0.2
11.20 × 10−0429.98 × 10−011.20 × 10−04$−0.2$0.0
11.20 × 10−0329.98 × 10−011.20 × 10−03$−0.2$0.0
11.20 × 10−0229.98 × 10−011.20 × 10−02$−0.2$0.0
11.20 × 10−0129.99 × 10−011.20 × 10−01$−0.1$0.0
101.20 × 10−0329.98 × 10+001.20 × 10−03$−0.2$$−0.4$
101.20 × 10−0229.98 × 10+001.20 × 10−02$−0.2$$−0.4$
101.20 × 10−0129.98 × 10+001.20 × 10−01$−0.2$$−0.4$
101.20 × 10+0041.01 × 10+011.21 × 10+000.90.5
1001.20 × 10−0259.87 × 10+011.17 × 10−02$−1.3$$−2.6$
1001.20 × 10−0159.87 × 10+011.17 × 10−01$−1.3$$−2.6$
1001.20 × 10+0071.06 × 10+028.75 × 10−016.3$−27.1$
1001.20 × 10+0164.51 × 10+011.17 × 10+01$−54.9$$−2.2$
10001.20 × 10−0145.01 × 10+023.02 × 10−02$−49.9$$−74.8$
10001.20 × 10+0045.01 × 10+023.02 × 10−01$−49.9$$−74.8$
10001.20 × 10+0145.22 × 10+023.00 × 10+00$−47.8$$−75.0$
10001.20 × 10+0282.79 × 10+031.05 × 10+02178.8$−12.3$
α$β$Mode$α^$$β^$%D α%D $β$
0.011.20 × 10−064−2.08 × 10−011.20 × 10−06$−2182.8$0.2
0.011.20 × 10−054−2.08 × 10−011.20 × 10−05$−2182.8$0.2
0.011.20 × 10−044−2.08 × 10−011.20 × 10−04$−2182.8$0.2
0.011.20 × 10−034−2.08 × 10−011.20 × 10−03$−2182.8$0.2
0.11.20 × 10−054−1.18 × 10−011.20 × 10−05$−218.4$0.2
0.11.20 × 10−044−1.18 × 10−011.20 × 10−04$−218.4$0.2
0.11.20 × 10−034−1.18 × 10−011.20 × 10−03$−218.4$0.2
0.11.20 × 10−024−1.18 × 10−011.20 × 10−02$−218.0$0.2
11.20 × 10−0429.98 × 10−011.20 × 10−04$−0.2$0.0
11.20 × 10−0329.98 × 10−011.20 × 10−03$−0.2$0.0
11.20 × 10−0229.98 × 10−011.20 × 10−02$−0.2$0.0
11.20 × 10−0129.99 × 10−011.20 × 10−01$−0.1$0.0
101.20 × 10−0329.98 × 10+001.20 × 10−03$−0.2$$−0.4$
101.20 × 10−0229.98 × 10+001.20 × 10−02$−0.2$$−0.4$
101.20 × 10−0129.98 × 10+001.20 × 10−01$−0.2$$−0.4$
101.20 × 10+0041.01 × 10+011.21 × 10+000.90.5
1001.20 × 10−0259.87 × 10+011.17 × 10−02$−1.3$$−2.6$
1001.20 × 10−0159.87 × 10+011.17 × 10−01$−1.3$$−2.6$
1001.20 × 10+0071.06 × 10+028.75 × 10−016.3$−27.1$
1001.20 × 10+0164.51 × 10+011.17 × 10+01$−54.9$$−2.2$
10001.20 × 10−0145.01 × 10+023.02 × 10−02$−49.9$$−74.8$
10001.20 × 10+0045.01 × 10+023.02 × 10−01$−49.9$$−74.8$
10001.20 × 10+0145.22 × 10+023.00 × 10+00$−47.8$$−75.0$
10001.20 × 10+0282.79 × 10+031.05 × 10+02178.8$−12.3$

Table 4 presents the prediction results for α and $β$ using the mode with maximum sensitivity using $σi(α,β)=d|λi|/d|α+μβ|$. We see that, by using this criterion, the predictions are very similar to those using Eq. (18), except for cases of assigned values of α ≤ 0.01. As before, the modes selected by this criterion differ from the two previous criteria used. On the basis of these comparisons, thus, we recommend that the sensitivity in Eq. (18) be used to isolate the “best” prediction of the α and $β$ values from an array of available modal data.

Table 4

Predicted $α^$ and $β^$ versus assigned values $α$ and $β$, with percent differences (%D), using the eigenmode (indicated) with maximum sensitivity $σi(α,β),d|λi|/d|α+μβ|$

α$β$Mode$α^$$β^$%D α%D $β$
0.011.20 × 10−0628.72 × 10−031.20 × 10−06−12.80.0
0.011.20 × 10−0528.72 × 10−031.20 × 10−05−12.80.0
0.011.20 × 10−0428.72 × 10−031.20 × 10−04−12.80.0
0.011.20 × 10−0328.73 × 10−031.20 × 10−03−12.70.0
0.11.20 × 10−0519.97 × 10−021.20 × 10−05−0.30.0
0.11.20 × 10−0419.97 × 10−021.20 × 10−04−0.30.0
0.11.20 × 10−0319.97 × 10−021.20 × 10−03−0.30.0
0.11.20 × 10−0229.87 × 10−021.20 × 10−02−1.30.0
11.20 × 10−0429.98 × 10−011.20 × 10−04−0.20.0
11.20 × 10−0329.98 × 10−011.20 × 10−03−0.20.0
11.20 × 10−0229.98 × 10−011.20 × 10−02−0.20.0
11.20 × 10−0129.99 × 10−011.20 × 10−01−0.10.0
101.20 × 10−0329.98 × 10+001.20 × 10−03−0.2−0.4
101.20 × 10−0229.98 × 10+001.20 × 10−02−0.2−0.4
101.20 × 10−0129.98 × 10+001.20 × 10−01−0.2−0.4
101.20 × 10+0041.01 × 10+011.21 × 10+000.90.5
1001.20 × 10−0259.87 × 10+011.17 × 10−02−1.3−2.6
1001.20 × 10−0159.87 × 10+011.17 × 10−01−1.3−2.6
1001.20 × 10+0059.93 × 10+011.17 × 10+00−0.7−2.5
1001.20 × 10+0121.08 × 10+021.16 × 10+017.5−3.6
10001.20 × 10−0117.48 × 10+026.71 × 10−02−25.2−44.1
10001.20 × 10+0017.48 × 10+026.71 × 10−01−25.2−44.1
10001.20 × 10+0117.46 × 10+026.66 × 10+00−25.4−44.5
10001.20 × 10+0221.40 × 10+037.61 × 10+0140.0−36.6
α$β$Mode$α^$$β^$%D α%D $β$
0.011.20 × 10−0628.72 × 10−031.20 × 10−06−12.80.0
0.011.20 × 10−0528.72 × 10−031.20 × 10−05−12.80.0
0.011.20 × 10−0428.72 × 10−031.20 × 10−04−12.80.0
0.011.20 × 10−0328.73 × 10−031.20 × 10−03−12.70.0
0.11.20 × 10−0519.97 × 10−021.20 × 10−05−0.30.0
0.11.20 × 10−0419.97 × 10−021.20 × 10−04−0.30.0
0.11.20 × 10−0319.97 × 10−021.20 × 10−03−0.30.0
0.11.20 × 10−0229.87 × 10−021.20 × 10−02−1.30.0
11.20 × 10−0429.98 × 10−011.20 × 10−04−0.20.0
11.20 × 10−0329.98 × 10−011.20 × 10−03−0.20.0
11.20 × 10−0229.98 × 10−011.20 × 10−02−0.20.0
11.20 × 10−0129.99 × 10−011.20 × 10−01−0.10.0
101.20 × 10−0329.98 × 10+001.20 × 10−03−0.2−0.4
101.20 × 10−0229.98 × 10+001.20 × 10−02−0.2−0.4
101.20 × 10−0129.98 × 10+001.20 × 10−01−0.2−0.4
101.20 × 10+0041.01 × 10+011.21 × 10+000.90.5
1001.20 × 10−0259.87 × 10+011.17 × 10−02−1.3−2.6
1001.20 × 10−0159.87 × 10+011.17 × 10−01−1.3−2.6
1001.20 × 10+0059.93 × 10+011.17 × 10+00−0.7−2.5
1001.20 × 10+0121.08 × 10+021.16 × 10+017.5−3.6
10001.20 × 10−0117.48 × 10+026.71 × 10−02−25.2−44.1
10001.20 × 10+0017.48 × 10+026.71 × 10−01−25.2−44.1
10001.20 × 10+0117.46 × 10+026.66 × 10+00−25.4−44.5
10001.20 × 10+0221.40 × 10+037.61 × 10+0140.0−36.6

## 5 Conclusions

In this study, we have presented a new method, to be used in contact resonance atomic force microscopy, to accurately determine the viscoelastic properties of a sample using a plate-like microsensor. The theoretical formulation, based on the Rayleigh–Ritz method and on an “unconventional” generalized eigenvalue problem, was presented and a numerical experiment was devised to validate the method capability to identify viscoelastic parameters from modal data. The proposed method performed extremely well for nondimensional sample stiffnesses in the range 0.01 ≤ α ≤ 100 and nondimensional damping values in the range $1.2×10−6≤β≤120$. For values of α ≥ 1000, the prediction error became increasingly larger with successively larger values of $β$.

This study opens up several avenues for future research. First, future work should include validation of the method against experimental results. While we anticipate that the method will perform well in a practical implementation, additional study must be performed to determine the validity of the sensitivity metric used in Eq. (18) and whether better metrics exist to accurately select the modal data used in the predictions, given experimental test data rather than numerical surrogates. In addition, in the future work, we would like to address: the optimal placement of the sensing tip given damping, the effect of large material damping, and the implications of fluid–structure interactions in liquid environments and their effects on the sensor estimates.

## Acknowledgment

This material is based in part upon work supported by the National Science Foundation under grant CMMI-1660448 (to R. T.) and grant CMMI-1847513 (to M. A.).

## Conflict of Interest

There are no conflicts of interest.

## References

1.
Rabe
,
U.
, and
Arnold
,
W.
,
1994
, “
Acoustic Microscopy by Atomic-Force Microscopy
,”
Appl. Phys. Lett.
,
64
(
12
), pp.
1493
1495
. 10.1063/1.111869
2.
Hurley
,
D. C.
,
Kopycinska-Müller
,
M.
, and
Kos
,
A. B.
,
2007
, “
Mapping Mechanical Properties on the Nanoscale Using Atomic-Force Acoustic Microscopy
,”
J. Miner. Metals Mater. Soc.
,
59
(
1
), pp.
23
29
. 10.1007/s11837-007-0005-8
3.
Yuya
,
P.
,
Hurley
,
D.
, and
Turner
,
J. A.
,
2008
, “
Contact-Resonance Atomic Force Microscopy for Viscoelasticity
,”
J. Appl. Phys.
,
104
(
7
), p.
074916
. 10.1063/1.2996259
4.
Killgore
,
J. P.
,
Yablon
,
D. G.
,
Tsou
,
A. H.
,
Gannepalli
,
A.
,
Yuya
,
P. A.
,
Turner
,
J. A.
,
Proksch
,
R.
, and
Hurley
,
D. C.
,
2011
, “
Viscoelastic Property Mapping With Contact Resonance Force Microscopy
,”
Langmuir
,
27
(
23
), pp.
13983
13987
. 10.1021/la203434w
5.
Qian
,
W.
,
Li
,
W.
,
Nguyen
,
C.
,
Johnson
,
T. J.
, and
Turner
,
J. A.
,
2020
, “
Quantitative Nanoscale Measurements of the Thermomechanical Properties of Poly-Ether-Ether-Ketone (PEEK)
,”
J. Polym. Sci.
,
58
(
11
), pp.
1544
1552
. 10.1002/pol.20190274
6.
Reggente
,
M.
,
Angeloni
,
L.
,
Passeri
,
D.
,
Chevallier
,
P.
,
Turgeon
,
S.
,
Mantovani
,
D.
, and
Rossi
,
M.
,
2020
, “
Mechanical Characterization of Methanol Plasma Treated Fluorocarbon Ultrathin Films Through Atomic Force Microscopy
,”
Front. Mater.
,
6
, p.
338
. 10.3389/fmats.2019.00338
7.
Aureli
,
M.
,
Ahsan
,
S. N.
,
Shihab
,
R. H.
, and
Tung
,
R. C.
,
2018
, “
Plate Geometries for Contact Resonance Atomic Force Microscopy: Modeling, Optimization, and Verification
,”
J. Appl. Phys.
,
124
(
1
), p.
014503
. 10.1063/1.5038727
8.
Timoshenko
,
S.
, and
Woinowsky-Krieger
,
S.
,
1987
,
Theory of Plates and Shells
, 2nd ed.,
McGraw-Hill
,
New York
.
9.
Leissa
,
A. W.
,
1969
, “
Vibration of Plates
,”
Report No. NASA SP-160, Technical Report
10.
Boresi
,
A. P.
, and
Schmidt
,
R. J.
,
2003
,
, 6th ed.,
Wiley
,
New York
.
11.
Reddy
,
J. N.
,
1984
,
Energy and Variational Methods in Applied Mechanics: With an Introduction to the Finite Element Method
,
Wiley
,
New York
.
12.
Mura
,
T.
, and
Koya
,
T.
,
1992
,
Variational Methods in Mechanics
,
Oxford University Press
,
New York
.
13.
Meirovitch
,
L.
,
1967
,
Analytical Methods in Vibrations
,
Macmillan
,
New York
.
14.
Horn
,
R. A.
, and
Johnson
,
C. R.
,
1990
,
Matrix Analysis
,
Cambridge University Press
,
Cambridge, UK
.
15.
Frazer
,
R. A.
,
Duncan
,
W. J.
, and
Collar
,
A. R.
,
1957
,
Elementary Matrices
,
Cambridge University Press
,
New York
.
16.
Abramowitz
,
M.
, and
Stegun
,
I. A.
,
1967
,
Handbook of Mathematical Functions: With Formulas, Graphs, and Mathematical Tables
,
Dover Publications Inc.
,
New York
.
17.
Press
,
W. H.
,
Teukolsky
,
S. A.
,
Vetterling
,
W. T.
, and
Flannery
,
B. P.
,
1992
,
Numerical Recipes in Fortran
, 2nd ed.,
Cambridge University Press
,
Cambridge, UK
.
18.
Dowling
,
N. E.
,
2012
,
Mechanical Behavior of Materials: Engineering Methods for Deformation, Fracture, and Fatigue
,
Pearson
,
Boston, MA
.