Numerical verification of the hydrodynamic reconstruction method for contact resonance atomic force microscopy

In the field of atomic force microscopy (AFM), accurately predicting the hydrodynamic forces acting on an AFM cantilever is of paramount importance. Of particular interest is predicting the surface-coupled hydrodynamic loading of the system. There are a myriad of advantages to performing AFM imaging while immersed in a liquid. These advantages include: the reduction of van der Waals forces,^1,2 the elimination of capillary forces,^3,4 and the ability to image biological samples in domains that more closely resemble their native environments.⁵ Additionally, the ability to accurately predict hydrodynamic forces acting on microstructures will benefit a wide range of applications including: microelectromechanical system (MEMS) design and operation,⁶ energy harvesting,^7–9 and biomimetic propulsion.^9–11 

The present work of estimating the hydrodynamic forces acting on cantilevered structures begins with Stokes, who developed the theory to estimate the forces acting on an infinite cylinder oscillating with small amplitude in an incompressible fluid.¹² From this work, one can show that the natural frequency of oscillating structures drops dramatically when moving from air to liquid environments. This is due to the increased added mass effect of the much denser liquid compared to the density of air. Tuck¹³ extended the theory of Stokes and provided a boundary integral technique to calculate the hydrodynamic force components of an infinitely thin ribbon. Using a similar boundary integral formulation as Tuck, Green and Sader¹⁴ calculated the hydrodynamic loading on an infinite cylinder immersed in an inviscid fluid and considered the effect of a nearby rigid wall. Tung et al.¹⁵ provided semi-analytical compact solutions of Green and Sader’s theory and verified them using three-dimensional fluid-structure interaction computations. The aforementioned theories consider the hydrodynamic effects of cantilevered structures oscillating near to, but not physically coupled to, a nearby substrate. In contact resonance (CR) AFM, and other contact modes of AFM, the cantilever is in physical contact with the substrate. Additionally, the mechanical properties of the substrate are unknown and are the primary parameters to be measured and characterized. CR AFM utilizes the surface-coupled cantilever resonances to infer the material properties of the sample under test. Estimation of the hydrodynamic forces can be done accurately in proximity to a surface, however, additional complexities and uncertainties are introduced while in contact with a surface of unknown properties.

When a cantilever that is freely vibrating near a surface, in one or more of its eigenmodes, is brought into repulsive contact with that surface the eigenfrequencies can shift a substantial amount. This shift in eigenfrequencies can greatly change the hydrodynamic forces (which are frequency dependent) acting on the system, see Fig. 1. Several previous works have tried to address this issue. Mirman and Kalinin¹⁶ and Payam¹⁷ considered the near surface effects but used semi-analytical hydrodynamic formulae derived for the unbounded fluid case. Ploscariu and Szoskiewicz¹⁸ proposed an experimental reconstruction method for inertial fluid loading which required extensive information about the cantilever and the experimental conditions. Furthermore, this method did not consider the near surface effects and did not account for physical cantilever-sample coupling effects. To overcome this knowledge gap, Tung et al.¹⁹ presented a Hydrodynamic Reconstruction Method (HRM) which can accurately characterize both the intertial and viscous hydrodynamic forces in the CR AFM system. This method does not require a priori knowledge of the geometry and material properties of the cantilever or the liquid it is immersed in. The Hydrodynamic Reconstruction Method has been validated experimentally in both air and water for one particular sample stiffness and cantilever geometry. In this paper, we present a numerical verification of the Hydrodynamic Reconstruction Method for various cantilever geometries and sample stiffnesses. A fully three-dimensional fluid-structure interaction (FSI) simulation has been developed with ADINA²⁰ software to accomplish this goal. We show that the HRM can accurately predict the hydrodynamic forces in CR AFM systems for a wide range of cantilever geometries and sample stiffnesses.

For a microcantilever submerged in a fluid oscillating in one or more of its eigenmodes near to, but not mechanically coupled to, a nearby rigid surface, Tung et al.¹⁹ showed that the hydrodynamic forces acting on the system can be discretely calculated using three experimental quantities. These quantities are the in-vacuo natural frequency ω_n, the wet natural frequency ω_wet, and the wet quality factor of vibration Q_wet. More specifically, the complex-valued hydrodynamic function Γ(Re) can be calculated discretely at each eigenmode of vibration using the aforementioned experimental quantities, where Re is the unsteady Reynolds number given by $Re=ρfωb24μf$⁠, ρ_f is the mass density of the fluid, ω is the frequency of oscillation of the microcantilever, b is the microcantilever width, and μ_f is the shear viscosity of the fluid. The hydrodynamic function Γ(Re) can be decomposed into real and imaginary parts Γ(Re) = Γ_r(Re) − iΓ_i(Re). The real part of the hydrodynamic function Γ_r is related to inertial forces and the imaginary part of the hydrodynamic function Γ_i is related to viscous forces in the fluid. For a given eigenmode of vibration, the discrete values of the hydrodynamic function are given by:¹⁹ 

and

The term χ is a non-dimensional parameter given by $χ=π4ρfb2ρA$⁠, where A is the cross-sectional area of the beam and ρ is the mass density of the beam. χ is typically a constant in AFM experiments. Additionally, so long as χ is a constant, its value has no effect on the hydrodynamic reconstruction method (HRM).¹⁹ 

To experimentally reconstruct the hydrodynamic function, the AFM cantilever is brought very close to the substrate, to capture surface-coupled hydrodynamic effects, and is then excited at several of its eigenmodes. For each eigenmode, the wet natural frequency ω_wet and the wet quality factor Q_wet are measured. The in-vacuo natural frequency ω_n can be measured in vacuum or approximated by using the natural frequency measured in air far from the sample surface. Once a suitable number of discrete points of the hydrodynamic function are calculated, a regression analysis is performed that accounts for the asymptotic nature of the hydrodynamic function. The general form of the regression functions for the real and imaginary parts of the hydrodynamic function are:¹⁹ 

and

Figure 2 depicts the procedure for using the HRM. First, the hydrodynamic function is measured at each of the discrete uncoupled eigenmodes of the system while vibrating near to, but uncoupled from, the sample surface using the measurement quantities described above. These data points are shown as green squares in Fig. 2. Next, Eqs. (3) and (4) are used in a regression analysis to determine the best fit to the data. Finally, with the best fit determined, the hydrodynamic function can be interpolated at any point within the data set for a given contact resonance frequency (Reynold’s number) of the original AFM system. The blue circles depict the hydrodynamic function estimates for a given set of in-contact natural frequencies.

With the hydrodynamic function successfully reconstructed for a given cantilever, the hydrodynamic forces acting on the cantilever can be accurately predicted at any arbitrary frequency. This predictive capability is especially useful in situations such as contact resonance (CR) AFM, where experimental measurement of the hydrodynamic function becomes impossible when coupled to a surface of unknown mechanical properties. Tung et al.¹⁹ experimentally validated the HRM for a single cantilever geometry in water for a single sample stiffness. In the following sections, we present a numerical verification of the HRM for several cantilever geometries and sample stiffnesses in a liquid water environment.

To conduct the numerical verification study, two distinct tip-less cantilever geometries were chosen. The first cantilever geometry chosen, which will henceforth be named “Cantilever 1”, was composed of silicon and had a length of 197 μm, a width of 29 μm, and a thickness of 2 μm. The mass density ρ and Young’s modulus E of the cantilever was assigned as ρ = 2320 kg ⋅ m⁻³ and E = 169 GPa. The second cantilever geometry, “Cantilever 2”, was also composed of silicon and had a length of 49 μm, a width of 29 μm, and a thickness of 2 μm. The mass density and Young’s modulus was identical to Cantilever 1. In both cantilevers, the Poisson’s ratio was assigned as ν = 0.25. The geometry of Cantilever 1 was chosen to mimic common AFM cantilevers that are regularly used in the field. Additionally, the geometry of Cantilever 1 was chosen such that we could easily benchmark our simulation data with previously published results by Basak et al.²¹ The geometry of Cantilever 2 was chosen to probe the limits of the HRM method. To simplify the hydrodynamic theory used in analysis it is often assumed that cantilever length greatly exceeds its width. Cantilever 2 has an aspect ratio closer to unity, therefore, the system hydrodynamics can no longer be considered two-dimensional. Cantilever 2 will test the ability of the HRM to capture three-dimensional hydrodynamic forces.

With each cantilever geometry, five distinct sample stiffnesses were selected to emulate various CR AFM sample contact conditions. In the simulations performed, sample contact is approximated by attaching linear springs at the free (distal) edge of the cantilever. The first stiffness was chosen to represent a soft contact. The second, third, and fourth stiffnesses were chosen such that the in-vacuo frequency sensitivity to contact stiffness changes of the first, second, and third eigenmodes were maximized, respectively. The frequency sensitivity to contact stiffness changes is defined as $∂fi∂k$⁠, where f_i is the i^th eigenfrequency of the cantilever-spring system and k is the sample stiffness. The fifth stiffness was chosen to represent an extremely stiff contact.

The numerical verification was conducted in three distinct steps. These steps are repeated for each cantilever geometry. The first step is to determine the first four in-vacuo natural frequencies and mode shapes of the cantilever. This was accomplished by performing a modal analysis of each geometry using ADINA software. The next step is to simulate the cantilever vibrating near to, but uncoupled from, a nearby rigid surface while immersed in fluid. The final step is to simulate the cantilever vibrating while coupled to the sample surface while immersed in fluid.

In the simulation, the cantilever is submerged in liquid water with a dynamic viscosity η = 8.59 × 10⁻⁴ kg ⋅ m⁻¹ ⋅ s⁻¹, a mass density ρ = 997 kg ⋅ m⁻³ and a bulk modulus K = 2.2 GPa. The cantilever is inclined at an angle of 11° relative to the rigid surface. More specifically, we define a unit vector ê_x parallel to the width direction of the cantilever originating from the centroid of the fixed end surface of the cantilever, a unit vector ê_y parallel to the length direction of the cantilever originating from the fixed end, and a rigid sample surface normal $n^$⁠. The cantilever is first oriented in relation to the rigid surface such that $(êx×êy)∥n^$⁠. Next, the cantilever is inclined such that angle between ê_y and the surface is 11°. The distance from the cantilever to the surface is set such that the minimum distance of a vector originating from the sample surface, parallel to the surface normal, and terminating at the free (distal) end of the cantilever is 7.5 μm and represents the minimum distance between the surface and any point on the cantilever.

Next, the transient ringdown of the cantilever is simulated for each of the four calculated in-vacuo mode shapes. For each simulation, the material points of the cantilever are given initial velocity conditions proportional to that mode shape. The maximum tip velocity is scaled to 10 nm ⋅ s⁻¹. From the ringdown data, the wet natural frequency ω_$w$et and wet quality factor Q_$w$et can be calculated using a single degree of freedom harmonic oscillator model and regression analysis, see Basak et al.²¹ Using these data, the discrete values of the hydrodynamic function can be calculated using Eq. (1) and Eq. (2). Using the discrete values of the hydrodynamic function, an estimate of the hydrodynamic function for arbitrary Reynolds number (frequency) values can be generated using Eqs. (3) and (4), see Fig. 2.

The final step of the verification process is to calculate the wet natural frequency and wet quality factor of the in-contact system. The cantilever is submerged in liquid water and oriented to the rigid surface as described previously. Linear springs are added to the free end of the cantilever that provide restoring forces proportional to the transverse displacement of the cantilever. The spring-coupled cantilever is simulated in-vacuo to determine the natural frequencies and mode shapes in the absence of fluid loading. Unique mode shapes are determined for each distinct spring value assigned. These in-vacuo mode shapes are then used as initial velocity conditions in the fully coupled FSI simulation, as described above.

From the spring-coupled liquid ringdown, the wet natural frequency ω_$w$et and wet quality factor Q_$w$et are determined, as described previously. For each spring value and mode shape, ω_$w$et, Eq. (1), and Eq. (3) are used to predict the system’s natural frequency in the absence of hydrodynamic effects ω_n. This estimate is directly proportional to the sample stiffness. A direct comparison can be made between the ω_n predicted using ω_wet and Eq. (3) and the ω_n calculated for the coupled in-vacuo case. If these two values of ω_n are identical, then the HRM has perfectly predicted the hydrodynamic loading present in the system. Additionally, the hydrodynamic damping in the system Q_$w$et can be predicted using ω_$w$et from the simulation/experiment, the previously predicted ω_n, Eq. (2), and Eq. (4). In the case of our numerical simulation, the only damping present in the system is due to hydrodynamics. Thus, a direct comparison between the Q_$w$et measured from the in-contact ringdown data and the Q_$w$et predicted using the HRM can be made. This is in contrast to an experimental case – where additional unknown damping would be introduced when the tip is in-contact with the sample surface.

Figure 3 shows the percentage difference between the numerically simulated and measured in-vacuo sample-coupled natural frequency of Cantilever 1 vs. the sample-coupled natural frequency predicted using the HRM for a range of nondimensional sample stiffnesses α. Here, we have defined the nondimensional stiffness as $α=kkc$ where k is the assigned sample stiffness and $kc=3EIL3$ is the static stiffness of the cantilever. It can be seen that the prediction error of the HRM is less than 8% for the entire range of α for the first three eigenmodes of vibration. Figure 4 shows the percentage difference between the in-liquid measured sample-coupled quality factor vs. the quality factor predicted by the HRM for Cantilever 1. The prediction error is less than 11% for all modes and nondimensional stiffness values.

Tabulated results can be seen in Table I. For each nondimensional stiffness value simulated, the quality factor measured from the simulation output Q_measured is reported along with the quality factor predicted by the HRM Q_predicted. Additionally, the wet natural frequency f_wet, the predicted in-vacuo natural frequency (with hydrodynamic forces removed) $fnpredicted$⁠, and the numerically simulated and measured in-vacuo natural frequency $fnmeasured$ are reported. Recall that f_wet represents the natural frequency extracted from the in-liquid, sample-coupled ringdown of the system and is affected by the hydrodynamic forces present in the system. In CR AFM, a measurement of the sample-coupled natural frequency, in the absence of fluid loading effects, must be obtained to generate accurate estimates of the sample’s material properties. We can conclude from the data that the HRM can be successfully used in situations where the cantilever geometry is similar to Cantilever 1. The nondimensional stiffness values probed reveal that the HRM prediction accuracy remains high across the entire range of experimentally achievable sample stiffnesses.

Figure 5 shows the percentage difference between the numerically simulated and measured in-vacuo sample-coupled natural frequency of Cantilever 2 vs. the sample-coupled natural frequency predicted using the HRM for a range of nondimensional sample stiffnesses α. The prediction accuracy of the HRM stays approximately below 10% for modes 2 and 3, however the prediction error is much greater for mode 1 with a maximum error of approximately 18%. It can be seen that the prediction error decreases with each subsequently higher eigenmode.

Figure 6 shows the percentage difference between the in-liquid measured sample-coupled quality factor vs. the quality factor predicted by the HRM for Cantilever 2. The prediction error is very large, with a maximum of approximately 50%. The same trend of increasing prediction accuracy with increasing mode that was observed in the natural frequency data is not seen in the quality factor data. We believe that the discrepancy in the prediction error, for both natural frequency and quality factor, is due to reduced dimensions of Cantilever 2. The length to width ratio of Cantilever 1 is approximately 6.8, while the length to width ratio of Cantilever 2 is approximately 1.7. Cantilever 2 behaves much more like a plate than an Euler-Bernoulli beam. Figure 7 shows the third transverse eigenmode calculated for Cantilever 1 and 2 – anticlastic curvature is apparent in Cantilever 2. The anticlastic motion of Cantilever 2 increases the complexity of the system’s overall hydrodynamic loading. Furthermore, the primary assumption in two-dimensional hydrodynamic theory that the cantilever be very long relative to its width is also violated by the geometry of Cantilever 2. These additional complications are not well captured by the HRM. It is clear from the numerical analysis conducted that the HRM should be used in cases where the cantilever is long relative to its width. Tabulated results for Cantilever 2 can be seen in Table II.

For both Cantilever 1 and Cantilever 2, the prediction error for the in-vacuo natural frequency decreases as mode number increases. As the frequency of oscillation, and thus Reynolds number, is increased the rate of change of the real part of the hydrodynamic function Γ_r decreases. For successively higher eigenmodes, the total change in Γ_r due to the frequency shift from the un-coupled to coupled case decreases. This fact decreases the predictive error. Additionally, lower nondimensional stiffness values (α values) show less predictive error in both Cantilever 1 and Cantilever 2. This is due to the fact that for lower α values the extrapolation range of the HRM is reduced. In the limit of an alpha value of zero, the HRM would be 100% accurate.

To predict the in-vacuo natural frequency, only Γ_r is needed. To predict Q, both Γ_r and Γ_i are needed. In fact, the prediction of Q involves the quotient of these two quantities. For this reason, we believe that the prediction of Q has an increased uncertainty and although similar logic can be applied to the accuracy of Γ_i no trend can be seen for the Q prediction.

To further explore the effects of geometry on the HRM, we present results for a triangular cantilever. The triangular lever was created by merging two Cantilever 1 geometries rotated ninety degrees with respect to each other. The fixed end of the triangular cantilever was then cropped to provide a straight edge with fixed boundary conditions (see Fig. 8(a)). Additionally, the thickness of the structure was reduced to 0.8 micrometers to simulate a softer triangular lever. The triangular cantilever has the same material properties as Cantilever 1. The aforementioned HRM analysis was performed for one nondimensional stiffness value α = 2. Here, we have chosen a nondimensional stiffness value representative of an experimental case is which sensitivity of the first in-contact eigenmode of vibration is desired.²² We have foregone analysis of higher nondimensional stiffness values for the triangular cantilever.

Table III shows the computed results. Overall, the HRM provides excellent predictive capability for the first three eigenmodes of vibration for the triangular cantilever. The prediction error of the in-vacuo surface-coupled natural frequency is less than 5% for the entire range. However, the prediction error does not decrease monotonically as observed in both Cantilever 1 and Cantilever 2. One potential reason for this discrepancy is that the triangular cantilever geometry, when viewed in its component two-dimensional cross-sections, can be thought of as 2 two-dimensional cross-sections vibrating close to one another. Basak and Raman²³ have performed an in-depth analysis of this very topic. Basak and Raman showed that the hydrodynamic function, for the case of 2 cross-sections vibrating close to one another in an unbounded fluid, is dependent on the distance between the cantilever cross-sections. The triangular cantilever can be viewed as an ensemble of two-dimensional cross-sections of various separation distances. In our case, there is an additional complexity caused by the proximity of the rigid substrate. These two factors combined undoubtedly produce a hydrodynamic function that is much more complex than the form used in Eqs. (3) and (4). Despite these complexities, the HRM predictive capability for the triangular lever is still comparable to, if not better than, that of Cantilever 1.

A fully three-dimensional finite element analysis was conducted to verify the Hydrodynamic Reconstruction Method¹⁹ (HRM). The verification was performed for three cantilevers with distinct geometries for 3 separate transverse bending modes across a wide range of sample stiffnesses. This numerical verification supplements the limited experimental validation that was performed by Tung et al.,¹⁹ in which a single sample stiffness was tested. We have found that the HRM performs with a high degree of accuracy for cantilever geometries that are long relative to their width for a wide range of sample stiffnesses. Poor performance was discovered for cantilevers that are short relative to their width.

This material is based upon work supported by the National Science Foundation under Grant No. 1660448.

The software “ADINA” has been used throughout the work to perform the numerical fluid-structure interaction (FSI) computations. The numerical simulations are designed to recreate the experimental conditions present in a typical CR AFM experiment and to verify the HRM. The FSI computations are performed by defining two different computational domains: (i) the solid domain and (ii) the fluid domain.

The solid domain consists of the AFM cantilever, which is constructed with three-dimensional linearly elastic elements. Specifically, the element group “threedsolid” is used for all solid elements. Additionally, the small displacement assumption has been made for all simulations. Fixed boundary conditions are used on one end of the cantilever, while the other end is free to vibrate or constrained using linear springs, depending on the simulation. In the simulations, the geometry of the tip has been excluded.

The solid microcantilever is then surrounded by three-dimensional Navier-Stokes elements. This is achieved by creating a sphere of three-dimensional Navier-Stokes elements, from which a void matching the dimensions of the microcantilever is created. The solid microcantilever is then inserted into this void and is thereby surrounded by the fluid elements. The radius of the sphere has been chosen through a convergence study discussed below. The outer boundaries of the sphere are assigned no-slip conditions to approximate a rigid wall. A fluid-structure boundary condition is defined between the outer surfaces of the microcantilever and the adjacent void boundaries of Navier-Stokes elements. The resulting FSI simulation is conducted using a direct two-way coupling between the fluid and solid elements.

The simulations are run for three distinct dynamics cases: (i) the microcantilever oscillating in an unbounded fluid, (ii) the microcantilever oscillating near to, but not coupled to, a rigid surface, and (iii) the microcantilever oscillating near to and coupled to a rigid surface with linear springs anchored at the end of the microcantilever.

In order to increase our confidence in the numerical simulations conducted, we have performed benchmark simulations against Basak et al.²¹ Benchmark simulations of the transient ringdown of the first and second bending modes were conducted in an unbounded fluid domain. For each case, our predicted natural frequencies were within 1.5% and quality factors were within 0.5% of the benchmark data.

Additionally, benchmark simulations were conducted for 11° tilt orientations of the beam near to, but uncoupled from, the rigid sample surface for the first and second bending mode. Table IV shows the results of the benchmark simulations. Our prediction results are within 10% of the quality factor benchmark for gap heights above 10 μm and below 12% for gap heights of 7.5 μm. Basak et al.²¹ did not provide natural frequency vs. gap data for the tilted or un-tilted simulations. A thorough convergence study has been conducted for each gap height.

Following the work of Basak et al.²¹ we have performed convergence studies investigating the effect of the domain size, time step size, and mesh size on the final solutions. The convergence criterion we followed throughout our work is as follows:

where x_1,i is the displacement of the cantilever tip at i=1,…N instants of time during the first oscillation cycle and x_2,i are the displacements calculated at exactly the same instants of time for computations that are performed with a denser mesh, finer time step, or larger domain.

To begin the convergence study, we first study the effect of mesh size on the output eigenfrequencies of the solid model. The element size is systematically decreased until there is less than a 1% change in subsequent estimations of the transverse eigenfrequencies for the first 4 transverse eigenfrequencies of the system. This mesh is used as a starting point for the solid model in the FSI simulation.

The FSI convergence study has been conducted in an iterative fashion. The iterative convergence study was conducted for a cantilever oscillating in an unbounded fluid. First, the solid and fluid elements are systematically increased until the simulation output meets the convergence criterion defined above. Then, the time step is decreased until the convergence criterion is met. Finally, the domain size is increased until the convergence criterion is met. The process is then repeated until the convergence criterion is met simultaneously for the chosen mesh size, time step size, and domain size.

Domain independence is achieved when the radius of the simulation domain $R=1.232L$⁠, where L is the length of the microcantilever. Time step independence is achieved when the time step dt = T/50, where T is the period of cantilever oscillation. Mesh independence is achieved when using 272,200 fluid elements and 4050 solid elements. Mesh subdivisions were refined by adjusting the divisions of individual edges in the simulation domain. Free form meshing using Delaunay triangulation was used for the fluid elements.

For simulations in bounded domains, i.e. near to the substrate, the convergence criterion was checked for increasingly finer fluid mesh sizes. The initial mesh size calculated for the unbounded fluid case was insufficient to meet the convergence criterion for simulations conducted near to the surface. The fluid mesh density in the vicinity of the cantilever was increased accordingly to ensure simulation convergence.

The transient motion of the tip has been recorded for each time step of the simulation. From the tip displacement curve, we can calculate the quality factor and wet natural frequency of the system. As described in the body of the paper, we assign small initial velocities to the cantilever proportional to each eigenmode of the system and simulate the cantilever ringdown. Due to the small oscillations of the microcantilever, each eigenmode ringdown of the system can be analyzed using a damped single degree of freedom oscillator equation. Nonlinear regression is used to extract the wet natural frequency and quality factor of the system using the equation

where, x(t) is the tip displacement as a function of time, A is the amplitude of oscillation, ω_n is the wet natural frequency and Q_n is the quality factor. Specifically, the tip displacement curve has been fit in “MATLAB” using the Gauss-Newton algorithm with bisquare regression.