*Article* **Assembly of Copolymer and Metal***−***Organic Framework HKUST-1 to Form Cu2–xS/CNFs Intertwining Network for Efficient Electrocatalytic Hydrogen Evolution**

**Yuanjuan Bai 1,2,\* , Yanran Li 2, Gonggang Liu <sup>1</sup> and Jinbo Hu 1,\***


**Abstract:** The construction of complex intertwined networks that provide fast transport pathways for ions/electrons is very important for electrochemical systems such as water splitting, but a challenge. Herein, a three dimensional (3-D) intertwined network of Cu2–xS/CNFs (x = 0 or 0.04) has been synthesized through the morphology-preserved thermal transformation of the intertwined PEG-*b*-P4VP/ HKUST-1 hybrid networks. The strong interaction between PEG chains and Cu2+ is the key to the successful assembly of PEG-*b*-P4VP nanofibers and HKUST-1, which inhibits the HKUST-1 to form individual crystalline particles. The obtained Cu2–xS/CNFs composites possess several merits, such as highly exposed active sites, high-speed electronic transmission pathways, open pore structure, etc. Therefore, the 3-D intertwined hierarchical network of Cu2–xS/CNFs displays an excellent electrocatalytic activity for HER, with a low overpotential (η) of 276 mV to reach current densities of 10 mA cm<sup>−</sup>2, and a smaller Tafel slope of 59 mV dec−<sup>1</sup> in alkaline solution.

**Keywords:** assembly; metal-organic frameworks; hydrogen evolution reaction; Cu2–xS

#### **1. Introduction**

Electrochemical water splitting is a critical energy conversion process for producing clean and sustainable hydrogen, which is composed of two half-cell reactions: oxygen evolution reaction (OER) and hydrogen evolution reaction (HER) [1–3] Electrolysis is a process that consumes electricity, therefore a catalyst is needed to reduce the potential. Pt-based electrocatalysts exhibit the best performance for H2 evolution in strongly acidic electrolytes, however their HER activities are substantially diminished under alkaline conditions [4,5] Consequently, considerable attempts have been devoted to developing sustainable, highly efficient, and non-precious electrocatalysts to meet a target of Pt-based catalysts replacement. In recent years, transition metal sulfides (TMSs) have been widely investigated and have demonstrated their potential as HER catalysts due to their high catalytic activity and chemical stability [6] Metal-organic frameworks (MOFs) are an intriguing class of porous crystalline materials constructed by the coordination of metal ions or clusters with organic linkers [7–11] For example, HKUST-1({[Cu3(C9H3O6)2(H2O)3]}n, aka Cu-BTC) is one of the very first permanently porous MOFs, which has been widely studied for multiple applications [12–14] The microporosity and tunable functionality of MOFs make them ideal template precursors to fabricate various transition metal-based carbon composites, including TMSs, by means of the pyrolysis process under a specific atmosphere [6,15,16] However, MOFs may undergo structural collapse during pyrolysis, resulting in the dramatic decrease in surface area, the wreck of well-defined MOF pore/channel structures, and the uneven distribution of active components, which significantly reduces the electrochemical performances of the MOFs derived materials. Therefore,

**Citation:** Bai, Y.; Li, Y.; Liu, G.; Hu, J. Assembly of Copolymer and Metal−Organic Framework HKUST-1 to Form Cu2–xS/CNFs Intertwining Network for Efficient Electrocatalytic Hydrogen Evolution. *Nanomaterials* **2021**, *11*, 1505. https://doi.org/ 10.3390/nano11061505

Academic Editors: Mohammad Malikan, Shahriar Dastjerdi and Danil N. Dybtsev

Received: 13 May 2021 Accepted: 3 June 2021 Published: 7 June 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

the performance of the existing MOF derivatives' electrocatalysts is still struggles to match precious-metal-based materials.

One of the best ways to solve the above issues is the assembly of crystalline MOF nanoparticles into well-aligned one, two, or three dimensional (1-, 2- or 3-D) superstructures, used as an appropriate template as both scaffold and directing agent for the epitaxial growth of MOFs [17–19] After a pyrolysis process, the MOFs superstructures would be converted into transition-metal-based materials with ordered stacking and porous nanostructure. The choice of the template is vital to the success of the assembly. Commonly used templates include metal-based materials (e.g., Te nanowire [20]), carbon materials (carbon nanotubes (CNTs), carbon nanofiber (CNFs) and so on) and polymers [21–23] Among of them, 1-D polymer nanowires or nanofiber templates with unique high surface-to-volume ratio, adjustable size, and surface modifiability attract wide research interests [24,25] In particular, 1-D polyacrylonitrile (PAN) nanowires substrate prepared by electrospinning is the most widely used in many research studies [23,26] For instance, Centrone et al. reported a microwave irradiation approach to grow MIL-47 on 1-D PAN nanowires [27] Han et al. proposed that the uniform and stable growth of MOFs on PAN nanowires can be enhanced through two types of chemical modification methods [23] These studies have shown that the surface-exposed functional groups play a key role in the MOFs' growth on the PAN substrate, which act as binding sites for metal species of the target MOFs. Nevertheless, due to the difficulties in controlling the number of functional groups, and the compatibility between functional groups and MOFs, MOFs may fall off and be unevenly distributed on the surface of PAN.

In 2012, Chen and coworkers reported 1-D poly (ethylene glycol)-*b*-poly(4-vinylpyridine) (PEG-*b*-P4VP) nanowires composed of a PEG shell and slightly crosslinked P4VP core [28] When the aspect ratio of nanowires is high enough, the 3-D intertwining network can be formed. The follow-up research confirmed that the PEG-*b*-P4VP intertwined nanowires are excellent templates for the nucleation and growth of small MOF crystals due to their high aspect ratio (an average diameter of 30 nm and length up to several microns) and original abundant nucleation sites for metal ions [29] In the electrochemical application, the construction of complex intertwined conductive networks is desired, which would provide fast transport pathways for mass and charge, and then enhance the performance of catalysts [30,31] It is likely that a TMSs/carbon composite intertwining network could be prepared by the vulcanization of the as-prepared PEG-*b*-P4VP@MOFs network; but this hypothesis remains unexplored.

Herein, we have successfully constructed a novel intertwined Cu2–xS/CNFs (x=0 or 0.04; CNFs = carbon nanofibers) network. The self-assembling strategy is used to control the growth of HKUST-1 nanocrystals on the 1-D super long PEG-*b*-P4VP nanofibers to form PEG-*b*-P4VP@HKUST-1 composites under room temperature. The PEG shell of the PEG-*b*-P4VP nanofibers provide numerous nucleation sites for HKUST-1, while the P4VP core with a positive charge repels the metal ions (precursor of the HKUST-1). As a result, hybridization of the PEG-*b*-P4VP nanofibers by the HKUST-1 occurs selectively in the PEG shell. Afterwards, the Cu2–xS/CNFs composite materials with an intertwined network structure can be prepared with un-changed morphology by using in-situ thermal calcination PEG-*b*-P4VP@ HKUST-1 precursors. Compared with Cu2–xS/C composites derived from individual HKUST-1 crystal, Cu2–xS/CNFs significantly enhance the HER performance, which is strongly related to the novel network superstructure.

#### **2. Results and Discussion**

Figure 1 shows the forming process of the 3-D network of Cu2–xS/CNFs nanocomposites. First, the core-crosslinked polymeric linear-like micelles with a PEG shell and a slightly crosslinked P4VP core, designated as PEG-*b*-P4VP nanofibers (NFs), is synthesized in a water/methanol mixed solvent according to a reported method [28]. Second, as the assembly substrate, PEG-*b*-P4VP NFs are dispersed in mixed solvent of water and N, N-Dimethylformamide (DMF) containing Cu(CH3COO)2·H2O (Cu(OAC)2). During this

process, Cu2+ ions were favorably absorbed onto the PEG chain surface by the electrostatic interaction. Along to the previous step, trimesic acid (BTC) ligands solutions are added to assemble with Cu2+ ions via the coordination interactions [12–14]. Then, a uniform PEG/HKUST-1 hybrid shell is generated in-situ on the surface of the P4VP core, resulting in a PEG-*b*-P4VP@HKUST-1 composites with core-shell structure. Ultimately, the PEG-*b*-P4VP@HKUST-1 composites are converted into a 3-D hierarchical network of Cu2–xS/CNFs after the sulfurization reaction between PEG-*b*-P4VP@HKUST-1 and thiourea in an argon flow. The Cu2–xS/CNFs hierarchical network have many advantages for HER in alkaline solution including: (1) the continuous and conductive network of the CNFs with a hierarchically porous structure can enable fast charge and mass transfer; and (2) the high active specific surface areas of the Cu2–xS/CNFs offer abundant electrocatalytic active sites that are easily accessible in HER.

**Figure 1.** Schematic diagram of the synthesis of 3-D Cu2–xS/CNFs hierarchical network electrode.

Field emission scanning electron microscope (FESEM) and transmission electron microscopy (TEM) in Figure 2a–c reveal that the PEG-*b*-P4VP nanofibers have a 3-D feature assembled by a series of stacked 1-D entangled ultrafine nanofibers with uniform diameters of about 30 nm and extra-long lengths above 2 μm. From the typical FESEM and TEM images in Figure 2d–f, we can see the structure of PEG-*b*-P4VP@HKUST-1 composites comprise a tightly crowded nanosized HKUST-1 crystal layer encapsulating the PEG*b*-P4VP nanofibers with uniform size and well-defined shape. Notably, the generated composites inherit the intertwining network structure of their precursors. No unassembled irregular HKUST-1 crystals are observed, suggesting that the nucleation and growth of MOF crystals are localized on the surface of the PEG-*b*-P4VP templates. In addition, it is found that the polymer chains penetrate the crystalline structure of HKUST-1, and the particle size was less than 100 nm. Consequently, the interaction is PEG-*b*-P4VP nanofibers and HKUST-1 crystals, as it is extended from the core surface to the deep crystalline structure, which is important for the stability of HKUST-1 nanoparticles. The calculation result shows that the HKUST-1 content in PEG-b-P4VP/HKUST-1 composites is around 91 wt % (experimental section, Supporting Information). Because the individual composite nanofibers were dispersible in the suspension, the resultant architecture of the 3-D network shows a loose morphology (Figure S1). Figure 3 shows the X-ray diffraction (XRD) patterns of the PEG-b-P4VP/HKUST-1 composites. Based on the HKUST-1 structure data [12–14], the observed pattern and the simulated pattern show high similarity, confirming the formation of pure crystalline HKUST-1. It is found that the packing density of the HKUST-1 crystals be tuned by simply adjusting precursor concentration (Figure S2). In other words, the number of nucleation sites are controllable. As the addition of Cu2+ and organic ligands increased, the number of HKUST-1 nanoparticles on the surface of PEG-*b*-P4VP nanofibers

gradually increased. In the end, the 0.24 g of Cu(OAC)2 and 0.21 g of BTC ligands were selected to apply for defined PEG-*b*-P4VP@HKUST-1.

**Figure 2.** Morphology of the as-prepared (**a**–**c**) PEG-*b*-P4VP nanofibers and (**d**–**f**) PEG-*b*-P4VP@HKUST-1 hybrids as observed by SEM and TEM with different magnifications. The inset TEM image in Figure 2c shows the diameter of PEG-*b*-P4VP nanowire is approximately 30 nm.

**Figure 3.** (**a**) The structure diagram of HKUST-1. (**b**) XRD pattern of the PEG-*b*-P4VP@HKUST-1 composites.

Characterization using XRD confirmed the structural identity and phase purity of the as-synthesized samples. As shown in Figure 4a, the XRD peaks of the final product after thermal treatment are well matched to those of the typical crystalline structures of Cu2S (JCPDS card no 09-0328) and Cu1.96S phase (JCPDS card no 29-0578). Therefore, the value of x can be determined to be 0 or 0.04. SEM (Figure 4b) and TEM (Figure 4c) images show that the 3-D intertwining network structure can be well-maintained for the Cu2–xS /CNFs even after calcination with the presence of thiourea under 400 ◦C. Closer observation of the Cu2–xS /CNFs (Figure 4d) reveals that Cu2–xS nanoparticle with coarse surfaces are uniformly packed along the ultra-long CNFs. No scattered Cu2–xS nanoparticles are observed, suggesting that the nanoparticles are sturdily attached to the CNFs. High-resolution TEM (HRTEM) (Figure 4e) observes clearly resolved and welldefined lattice fringes, revealing the high crystallinity in agreement with XRD results. The distance between the adjacent lattice planes is 0.339 nm, corresponding to standard spacing of (302) plane of Cu2S. By contrast, for individual HKUST-1 crystals prepared from the assembly of Cu(OAc)2 and BTC ligands without the PEG-*b*-P4VP template, the collapse of the porous structure and random aggregation of the nanoparticles occurred after sulfurization (Figure S3). The product of the individual HKUST-1 crystals after

sulfurization treatment is denoted as Cu2–xS/C. The X-ray spectrometry (EDS) spectrum of the Cu2–xS/CNFs (Figure S4a) shows Cu and S signals attributable to Cu2–xS and C signals to CNFs. The corresponding EDS element mapping images (Figure S4b–e) demonstrate the homogeneous distribution of Cu, S, and C elements in the hierarchical nanostructure.

**Figure 4.** (**a**) XRD pattern, (**b**) SEM, (**c**,**d**) TEM and (**e**) HRTEM images of Cu2–xS/CNFs composites (inset image is the intensity plot of d-spacing for the (302) plane of Cu2–xS in e).

Electrochemical water splitting is a surface chemical process. The surface of materials plays a major role in determining catalysts behaviors. Consequently, the elemental compositions and valence states of Cu2–xS/CNFs network are studied in detail by X-ray photoelectron spectroscopies (XPS). Figure 5a displays clear a survey spectrum, which confirms the successful synthesis of Cu2–xS/CNFs. Regarding the Cu(2p) XPS spectrum as exhibited in Figure 5b, two strong peaks for Cu 2p orbit can be assigned to the Cu 2p3/2 and Cu 2p1/2 [32,33]. The binding energy of Cu 2p3/2 can be fitted into two peaks simulated at 932.7 and 934.6 eV [32,33]. The shakeup satellite peaks appearing at 943.7 eV might result from the partial surface oxidation of the samples when they made contact with air. In the S 2p region of the high resolution XPS spectra in Figure 5c, the binding energies (BEs) at 162.8 and 161.8 eV are ascribed to S 2p3/2 and S 2p1/2, respectively [34]. The two doublet peaks simulated at 163.7 and 164.9 eV correspond to the BEs of S-S bond. This result provided further evidence of the surface oxidation of Cu2–xS nanoparticles caused by air exposure. Raman spectroscopy in Figure S5 is used to investigate the nature of the carbon of Cu2–xS/CNFs. The two peaks located at 1351 and 1582 cm−<sup>1</sup> match well with the D (disordered carbon) and G (graphitized carbon) bands of carbon, respectively. The D-band to G-band intensity ratio (ID/IG) of the samples is calculated to be around 1.07, implying the limited graphitization degree of Cu2–xS/CNFs.

**Figure 5.** XPS survey spectra (**a**) of Cu2–xS/CNFs. High-resolution (**b**) Cu 2p and (**c**) S 2p XPS spectra for Cu2–xS/CNFs.

The electrocatalytic HER performance of Cu2–xS/CNFs is tested in 1.0 M KOH in a typical three-electrode setup with a scan rate of 5 mV s−1. Pt/C, bare glassy carbon (GC) electrode, and Cu2–xS/C composites are also examined for comparison. Figure 6a shows the linear sweep voltammetry (LSV) curves on the reversible hydrogen electrode (RHE) scale after iR correction. The bare GC electrode shows nearly no catalytic activity, while Pt/C exhibits excellent activity for HER. Cu2–xS/C composites show the limited HER performance with an overpotential of 432 mV at the current density of 10 mA cm−2. However, the as-prepared Cu2–xS/CNFs exhibit much superior HER activity compared to Cu2–xS/C composites, and only demand overpotentials of 276 and 337 mV to achieve current densities of 10, and 50 mA cm<sup>−</sup>2, respectively, fully demonstrating the importance of hierarchical intertwining network structures in the catalyst. The Tafel slope is an important indicator of electron-transfer kinetics and the rate-determining step in the HER process. The linear portions of the Tafel plots are fitted to the Tafel equation: η = b log j + a, where j is the current density and b is the Tafel slope. As shown in Figure 6b, Pt/C, Cu2–xS/C composites, and Cu2–xS/CNFs has Tafel slopes of 30, 102, and 59 mV dec−1, respectively. Based on the kinetic mechanism for alkaline HER, the HER over Cu2–xS/CNFs follows a Volmer−Heyrovsky mechanism, indicating that the Heyrovsky process (electrochemical desorption of hydrogen atoms) is the rate determining step [4,5,35] Electrochemical stability is another vital criterion to evaluate the new catalysts. As observed in Figure 6c, the Cu2–xS/CNFs retained 94% of its initial HER activity after a 12 h test, while the overpotential increase of the Cu2–xS/C electrode under the same condition was as high as 24% after only 8 h electrocatalysis. Meanwhile, the cyclic voltammetry (CV) durability tests of the Cu2–xS/CNFs electrode for HER in alkaline media is carried out at a scan rate of 100 mV s−<sup>1</sup> (Figure 6d). After 2000 CV sweeps, the polarization curve shows a negligible difference compared with the initial curve, suggesting superior stability of Cu2–xS/CNFs in the long-term electrochemical process.

**Figure 6.** (**a**) LSV polarization curves and (**b**) corresponding Tafel plots of the Cu2–xS/CNFs, Cu2–xS/CNFs, 20% Pt/C, and a bare glassy carbon electrode derived from the polarization curves. (**c**) Chronopotentiometric curves at 10 mA cm−<sup>2</sup> for Cu2–xS/CNFs and Cu2–xS/C composites over 10 h. (**d**) Accelerated HER polarization curves of Cu2–xS/CNFs. (**e**) Nyquist plots of Cu2–xS/CNFs and Cu2–xS/C composites at 300 mV overpotential in 1 M KOH (inset image is the magnified Nyquist plot of Cu2–xS/CNFs). (**f**) Capacitive current densities at 0.2 V vs. RHE as a function of scan rate of Cu2–xS/CNFs and Cu2–xS/C composites.

To further study the influence factors of catalytic kinetics, electrochemical impedance spectroscopy (EIS) and cyclic voltammetry (CV) are conducted. EIS analysis is performed to investigate the electrode/electrolyte interface properties of two catalysts with different structures. Figure 6e shows the Nyquist plots of the Cu2–xS/CNFs and Cu2–xS/C electrode at an overpotential of 300 mV (vs. RHE) in the frequency range of 100 kHz to 0.1 Hz. It can be observed that the charge-transfer resistance (Rct) is only 6.4 Ω for Cu2–xS/CNFs,

but 58 Ω for the Cu2–xS/C electrode, in accordance with the HER results. This reflects the highly faradaic efficiency and fast electron transfer of the catalysts during the reaction. The smaller value of Rct of Cu2–xS/CNFs electrocatalysts may be ascribed to their 3-D hierarchical intertwining network structure, which increased the contact of the active sites with the electrolyte, leading to a significant acceleration of the interfacial electrocatalytic reactions. We further tested the electrochemical double layer capacitances (Cdl) (Figure 6f) of the catalysts by a simple cyclic voltammetry method, to relate the catalytic activity with the electrochemical surface area (ECSA) (Figure S6). The values are measured to be 29.3 and 6.6 mF cm−2, for Cu2–xS/CNFs and Cu2–xS/C, respectively, revealing the higher active surface area of Cu2–xS/CNFs. The large ECSA indicates that Cu2–xS/CNFs exposes higher accessible active sites, which is one of the possible reasons for the excellent HER performance.

#### **3. Conclusions**

In conclusion, Cu2–xS/CNFs with 3-D hierarchical intertwining network structures have successfully been designed and fabricated. They are derived from a 3-D network of PEG-*b*-P4VP@HKUST-1. The key to the successful assembly of the PEG-*b*-P4VP and HKUST-1 is the strong interaction between Cu2+ ions and PEG chains. Benefiting from the unique hierarchical structure and uniformly distributed active sites, the as-prepared Cu2–xS/CNFs Cu2–xS/CNFs exhibit high intrinsic HER activity in alkaline medium. The overpotential is 276 and 337 mV at the current of 10 and 50 mA cm<sup>−</sup>2, respectively. The Tafel slope is calculated to be 59 mV/decade, and the high activity can be maintained for more than 12 h. The assembly strategy of HKUST-1 and PEG-*b*-P4VP herein can be extended to the hybrid of other MOF-based materials and copolymer, which shows great potential, not only for catalysts, but also for gas sensors, energy storage, and environmental science.

**Supplementary Materials:** The following are available online at https://www.mdpi.com/article/10 .3390/nano11061505/s1, Figure S1: SEM and optical images of PEG-*b*-P4VP/HKUST-1; Figure S2: TEM images of PEG-*b*-P4VP@HKUST-1 composites at different weight of Cu(OAc)2 and BTC; Figure S3: SEM images of individual HKUST-1 crystals and Cu2–xS/C composites; Figure S4: EDS spectrum of the Cu2–xS/CNFs and the corresponding element mapping image; Figure S5: Raman spectrum of Cu2–xS/CNFs; Figure S6: CV curves of the Cu2–xS/CNFs and Cu2–xS/C composites with different scan rates.

**Author Contributions:** Conceptualization, methodology, writing—original draft preparation, data curation, writing—review and editing were done by Y.B. and J.H.; software, formal analysis, investigation, visualization were done by G.L. and Y.L.; validation, supervision, funding acquisition and project administration were done by Y.B. and Y.L. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the National Natural Science Foundations of China (no. 22005348, 21908251), China Postdoctoral Science Foundation (no. 2020M670970), Hunan Provincial Natural Science Foundation of China (no. 2020JJ5961), Scientific Research Project of Education Department of Hunan Province (no. 19C1915), and research startup foundation of Central South University of Forestry and Technology (no. 2017YJ003).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Data available in a publicly accessible repository.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


## *Article* **Nonlinear Free and Forced Vibrations of a Hyperelastic Micro/Nanobeam Considering Strain Stiffening Effect**

**Amin Alibakhshi 1, Shahriar Dastjerdi <sup>2</sup> , Mohammad Malikan 3,\* and Victor A. Eremeyev 3,4**


**Abstract:** In recent years, the static and dynamic response of micro/nanobeams made of hyperelasticity materials received great attention. In the majority of studies in this area, the strain-stiffing effect that plays a major role in many hyperelastic materials has not been investigated deeply. Moreover, the influence of the size effect and large rotation for such a beam that is important for the large deformation was not addressed. This paper attempts to explore the free and forced vibrations of a micro/nanobeam made of a hyperelastic material incorporating strain-stiffening, size effect, and moderate rotation. The beam is modelled based on the Euler–Bernoulli beam theory, and strains are obtained via an extended von Kármán theory. Boundary conditions and governing equations are derived by way of Hamilton's principle. The multiple scales method is applied to obtain the frequency response equation, and Hamilton's technique is utilized to obtain the free undamped nonlinear frequency. The influence of important system parameters such as the stiffening parameter, damping coefficient, length of the beam, length-scale parameter, and forcing amplitude on the frequency response, force response, and nonlinear frequency is analyzed. Results show that the hyperelastic microbeam shows a nonlinear hardening behavior, which this type of nonlinearity gets stronger by increasing the strain-stiffening effect. Conversely, as the strain-stiffening effect is decreased, the nonlinear frequency is decreased accordingly. The evidence from this study suggests that incorporating strain-stiffening in hyperelastic beams could improve their vibrational performance. The model proposed in this paper is mathematically simple and can be utilized for other kinds of micro/nanobeams with different boundary conditions.

**Keywords:** hyperelastic micro/nanobeam; extended modified couple stress theory; strain-stiffening effect; nonlinear frequency response

#### **1. Introduction**

For many decades, vibration analysis of mechanical structures was a major topic among scientists [1–10]. Over recent decades, a surge of interest in studying hyperelastic materials was shown. The main characteristic of hyperelastic materials is that their strainstress diagram is nonlinear and may undergo large deformations [11–13]. Hyperelastic materials play a vital role in soft systems and structures, e.g., soft robotics [14], human organs [15,16], soft actuators [17,18], soft sensors [19,20], and soft energy harvesters [19–22]. Data from previous studies show that various mechanical structures such as beams, plates, membranes, and shells were made of hyperelastic materials [23–35]. It was reported that hyperelastic beams are an appropriate candidate to fabricate systems with high performance. For this reason, this paper focuses on beams made of hyperelastic materials. In light

**Citation:** Alibakhshi, A.; Dastjerdi, S.; Malikan, M.; Eremeyev, V.A. Nonlinear Free and Forced Vibrations of a Hyperelastic Micro/Nanobeam Considering Strain Stiffening Effect. *Nanomaterials* **2021**, *11*, 3066. https:// doi.org/10.3390/nano11113066

Academic Editor: Kosei Ueno

Received: 24 October 2021 Accepted: 12 November 2021 Published: 14 November 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

of the applications and properties of hyperelastic beams, it is becoming extremely difficult to ignore their investigation in different situations.

There are numerous published works that investigate dependent and time-independent responses of hyperelastic beams. For example, the nonlinear postbuckling of a hyperelastic beam-like structure was investigated by Lubbers et al. [36]. They employed the neo-Hookean hyperelastic model in conjunction with empirical tests and the finite element technique in their study. Wang and coworkers [37] studied nonlinear vibration of hyperelastic beams utilizing time history diagrams and frequency responses, who employed a compressible neo-Hookean constitutive law. He and coworkers [38] developed the Euler–Bernoulli beam model in a new finite strain framework to model a neo-Hookean hyperelastic beam. Xu and Liu [39] proposed an improved method to dynamically explore the response of a beam-like hyperelastic structure, where a Yeoh model was utilized to capture the material nonlinearity. Nonlinear dynamic characteristics of a soft hyperelastic beam were investigated by Wang et al. [40], employing a compressible neo-Hookean model and variational approach. Wang and Zhu [41] studied the nonlinear oscillation of a harmonically excited hyperelastic beam. They utilized the frequency-amplitude response, time histories, and a compressible neo-Hookean model in their investigation. The finite bending of a beam made of hyperelastic materials was analyzed by Bacciocchi and Tarantino [42]. They utilized a compressible Mooney-Rivlin hyperelastic material model to physical nonlinearity of the beam. Dadgar–Rad and Sahraee [43], by considering the incompressibility condition, investigated the large deformation response of a beam made of hyperelastic materials, where a neo-Hookean model was employed as the hyperelastic constitutive model. Bacciocchi and Tarantino [44] conducted a finite anticlastic bending analysis of hyperelastic beams using two hyperelastic models, namely Mooney–Rivlin and Saint Venant–Kirchhoff. Lanzoni and coworker [45] studied the nonuniform bending of a beam made of the hyperelastic beam, taking the Mooney–Rivlin into account. The large deformation response of hyperelastic beams was explored by Dadgar–Rad and Firouzi [46]. They incorporated Fung's quasilinear viscoelasticity theory and Mooney-Rivlin model.

Results from earlier studies demonstrate that few researchers addressed the modelling of hyperelastic beams with the strain-stiffening effect. Furthermore, previous studies have notably investigated a beam-like hyperelastic structure on a large scale and have not considered the hyperelastic beams in micro/nanoscales. However, fabrication of such beams in smaller scales was feasible, and hence analyzing hyperelastic micro/nanobeam and proposing more sophisticated theories should be developed for such structures. A challenging problem that arises in this domain is accurate modelling for hyperelasticity in micro/nanoscales. More specifically, in nanoscale, it is necessary to capture the size effect. Because hyperelastic materials may undergo large deformation and large rotation, these conditions should be considered on micro/nanoscale. One of the problems that it investigates in hyperelasticity is the strain-stiffening effect. This effect may improve or limit the performance of hyperelastic micro/nanobeams. Therefore, incorporating strain-stiffening with simple mathematical modelling in micro/nanoscale is essential. Specifically, to our knowledge, no study has considered large deformation, strain-stiffening, and moderate rotation for hyperelastic micro/nanobeams.

This paper aims to propose a sophisticated model for a micro/nanobeam made of hyperelastic materials that incorporate the small-scale and strain-stiffening effects of nonlinear elasticity. The nonlinear equations of motion are derived via Hamilton's principle and an extended von-Kármán theory. The frequency-amplitude plot and nonlinear resonance plot are presented by considering different system parameters. The results are discussed in detail, and influential parameters on free and forced vibrations of the hyperelastic micro/nanobeam are identified.

#### **2. Governing Equations**

The schematic view of the hyperelastic micro/nanobeam is illustrated in Figure 1, where the length, width, and height of the beam are denoted by *L*, *b*, and *d*, respectively. A clamped-clamped boundary condition is assumed to the beam, and a harmonic transverse mechanical load is applied to it. It is considered that the length of the beam is much greater than the depth. In addition, the shear deformation and rotary inertia are neglected. Thus, we use the Euler–Bernoulli (E-B) beam theory to define the displacement field.

The displacement field for the beam is established according to the Euler–Bernoulli beam equation, namely,

$$\begin{aligned} u\_x &= -z \frac{\partial w(x,t)}{\partial x} \\ u\_y &= 0 \\ u\_z &= w(x,t) \end{aligned} \tag{1}$$

where *w*(*x*, *t*) stands for the transverse displacement of any point on the neutral axis. The strain-displacement relations originated for the Euler–Bernoulli beam theory are modelled based on an extended von Kármán equation, in which large deformation, moderate rotation, and transverse strain are included, namely [47,48]

**Figure 1.** Schematic representation of a clamped-clamped hyperelastic beam.

$$\begin{aligned} \varepsilon\_1 &= \frac{1}{2} \left( \frac{\partial w}{\partial \mathbf{x}} \right)^2 - z \frac{\partial^2 w}{\partial \mathbf{x}^2} \\ \varepsilon\_3 &= \frac{1}{2} \left( \frac{\partial w}{\partial \mathbf{x}} \right)^2 \end{aligned} \tag{2}$$

Other components of the extended von Kármán equation are equal to zero.

The strain energy of the hyperelastic micro/nanobeam is decomposed into two parts, i.e., the potential due to the hyperelasticity and the potential due to small-scale effects.

For hyperelastic materials, a strain energy function is used to obtain the strain energy of the system. Numerous hyperelastic strain energy functions can capture the strain stiffening, for instance, the standard Gent, the Arruda–Boyce, and modified versions of the Standard Gent model [49–51]. In this work, for simplicity, a standard Gent model is considered, in which the strain-stiffening effect is incorporated, namely [52,53]

$$\Psi\_1 = \frac{\mu}{2} \left[ (I\_1 - \mathfrak{Z}) + \frac{1}{2} \frac{1}{f\_m} (I\_1 - \mathfrak{Z})^2 + \dots + \frac{1}{(n+1)f\_m^n} (I\_1 - \mathfrak{Z})^{n+1} \right] \tag{3}$$

where *μ* is the shear modulus; *I*<sup>1</sup> denotes the first invariant of the deformation tensor; *Jm* is a dimensionless parameter that is called the stiffening parameter.

For simplicity, the second-order expansion of the standard Gent model is utilized, such that

$$\Psi\_1 = \frac{\mu}{2} \left[ (I\_1 - 3) + \frac{1}{2} (I\_1 - 3)^2 \right] \tag{4}$$

The first principal invariant of the right Cauchy–Green deformation tensor in terms of the extended von Kármán strains is formulated as [54]

$$I\_1 = \mathfrak{Z}(\varepsilon\_1 + \varepsilon\_2 + \varepsilon\_3) + \mathfrak{Z} \tag{5}$$

Substituting Equation (2) into Equation (5), the first principal invariant is reformulated as

$$I\_1 = 2\left[\left(\frac{\partial w}{\partial \mathbf{x}}\right)^2 - z\frac{\partial^2 w}{\partial \mathbf{x}^2}\right] + 3\tag{6}$$

Substituting Equation (6) into Equation (4), the Gent strain energy function as a function of transverse displacement is obtained below

$$\Psi\_1 = \int\_0^L \left[ \mu A \left( \frac{\partial w}{\partial \mathbf{x}} \right)^2 + \frac{\mu A}{f\_m} \left( \frac{\partial w}{\partial \mathbf{x}} \right)^4 + \frac{\mu}{f\_m} I \left( \frac{\partial^2 w}{\partial \mathbf{x}^2} \right)^2 \right] \mathbf{d} \mathbf{x} \tag{7}$$

It is mentioned that Equation (7) was obtained by considering the following relations

$$\begin{aligned} I &= \int\_A z^2 \, \mathrm{d}y \mathrm{d}z = b \frac{d^3}{12} \\ A &= \int\_A \, \mathrm{d}y \mathrm{d}z = bd \\ 0 &= \int\_A z \, \mathrm{d}y \mathrm{d}z \end{aligned} \tag{8}$$

In Equation (8), *A* is the cross-section area, and *I* represents the second moment of the cross-section.

The potential of the small-scale effect is considered through the use of an extended modified couple stress theory, such that [47]

$$\Psi\_2 = \frac{1}{2} \left( 2\mu A \ell^2 \right) \int\_0^L \left( \frac{\partial^2 w}{\partial x^2} \right)^2 \mathrm{d}x \tag{9}$$

where is a length-scale parameter.

Comparing Equation (9) with previous studies, for the moderate rotation, a coefficient 2 appears in the equation in comparison to the small rotation [55].

Finally, the total strain energy of the hyperelastic micro/nanobeam is written as

$$
\mathcal{U}\_s = \Psi\_1 + \Psi\_2 \tag{10}
$$

The moving beam generates the kinetic energy in the system, which is formulated as

$$\mathrm{d}I\_k = \frac{1}{2}\rho A \int\_0^L \left(\frac{\partial w}{\partial t}\right)^2 \mathrm{d}t \tag{11}$$

where *ρ* stands for the mass-density of the hyperelastic beam.

The transverse applied periodic loading does the work of the following form

$$\mathcal{W}\_{\rm F} = \int\_0^L F \cos(\omega t) w \,\mathrm{d}x \tag{12}$$

in which *F* is the amplitude and *ω* indicates the excitation frequency. The work generated from the viscous damping is expressed as

$$\mathcal{W}\_{\rm D} = -c\_{\rm D} \int\_{0}^{L} \frac{\partial w}{\partial t} w \, \mathrm{d}x \tag{13}$$

where *cD* is the viscous damping coefficient.

To derive boundary conditions and governing equation, Hamilton's principle is utilized, namely

$$
\delta \int\_{t\_1}^{t\_2} [\mathcal{U}\_k - \mathcal{U}\_\mathcal{S}] \mathbf{d}t + \delta \int\_{t\_1}^{t\_2} [\delta \mathcal{W}\_F + \delta \mathcal{W}\_D] \mathbf{d}t = 0 \tag{14}
$$

Substituting Equations (10)–(13) into Equation (14), we obtain the following equations

$$\rho A \frac{\partial^2 w}{\partial t^2} + C\_D \frac{\partial w}{\partial t} + \frac{2\mu I}{\int\_{\mathcal{U}} \frac{\partial^4 w}{\partial x^4}} + 2\mu A \ell^2 \frac{\partial^4 w}{\partial x^4} - 2\mu A \frac{\partial^2 w}{\partial x^2} - \frac{12\mu A}{\int\_{\mathcal{U}} \left(\frac{\partial w}{\partial x}\right)^2} \frac{\partial^2 w}{\partial x^2} = F \cos(\omega t) \tag{15}$$

and boundary conditions for the double-clamped micro/nanobeam

$$w(0) = 0, \; w(L) = 0, \; \frac{\text{d}w(0)}{\text{d}x} = 0, \; \frac{\text{d}w(L)}{\text{d}x} = 0 \tag{16}$$

The above equations are made dimensionless to simplify and generalize the vibration analysis of the micro/nanobeam. The following nondimensional quantities are introduced, such that

$$\begin{aligned} \mathfrak{X} &= \frac{\mathfrak{x}}{L}, \mathfrak{d} = \frac{w}{L}, \mathfrak{f} = t\sqrt{\frac{\mu I}{\rho AL^4}}, \mathfrak{f} = \frac{cL^4}{\mu I}\sqrt{\frac{\mu I}{\rho AL^4}}, \,\,\Omega = \omega\sqrt{\frac{\rho AL^4}{\mu I}}\\ \eta\_1 &= \frac{2\mu A\ell^2}{\mu I}, \,\eta\_2 = -\frac{2\mu AL^2}{\mu I}, \,\,\mathfrak{f} = -\frac{12\mu AL^2}{\mu I \mu}, \,\,\widetilde{\mathfrak{f}} = \frac{FL^3}{\mu I} \end{aligned} \tag{17}$$

Utilizing the above equations, the dimensionless partial differential equation governing the transverse vibration of the beams is obtained as (the hat notation is omitted for convenience).

$$\frac{\partial^2 w}{\partial t^2} + c \frac{\partial w}{\partial t} + \frac{1}{f\_m} \frac{\partial^4 w}{\partial x^4} + \eta\_1 \frac{\partial^4 w}{\partial x^4} + \eta\_2 \frac{\partial^2 w}{\partial x^2} + \beta \left(\frac{\partial w}{\partial x}\right)^2 \frac{\partial^2 w}{\partial x^2} = F \cos(\Omega t) \tag{18}$$

Subsequently, the boundary conditions become

$$w(0) = 0, \; w(1) = 0, \; \frac{\text{d}w(0)}{\text{d}x} = 0, \; \frac{\text{d}w(1)}{\text{d}x} = 0 \tag{19}$$

The system is continuous, and therefore there are infinite modes of vibration. In this paper, the first mode is considered only, with the aid of the separation of variable technique and the Galerkin method. Based on the separation of variable technique, we assume the transverse response is approximated as

$$w(\mathbf{x}, t) = \mathcal{W}(\mathbf{x})q(t) \tag{20}$$

in which *q*(*t*) is the time-dependent coordinate of vibration; *W*(*x*) stands for the mode shape of a double-clamped beam that is given below [56]

$$\mathcal{W}(\mathbf{x}) = \sqrt{\frac{2}{3}} \left[ 1 - \cos(2|\pi \,\mathbf{x}) \right] \tag{21}$$

The function expressed in Equation (21) satisfies conditions in Equation (19).

According to the Galerkin method, Equation (20) is substituted in Equation (18), and the resulting equation is multiplied by Equation (21), and integration over [0 1] is taken, which results in ..

$$
\ddot{q} + c\dot{q} + \omega\_0^2 q + \alpha q^3 = f \cos(\Omega t) \tag{22}
$$

In which

$$
\omega\_0 = \left(\frac{\int\_0^1 \left\{\eta\_1 W^{\prime\prime\prime\prime} W + \frac{1}{\int\_0^1 W^{\prime\prime\prime} d\mathbf{x}} W^{\prime\prime\prime} W^{\prime} W\right\} dx}{\int\_0^1 W^2 dx}\right)^{\frac{1}{2}}
$$

$$
\mathfrak{a} = \frac{\int\_0^1 (\beta W^{\prime\prime} W^{\prime} W^{\prime\prime} W) dx}{\int\_0^1 W^2 dx} \tag{23}
$$

$$
f = \frac{\int\_0^1 (FW) dx}{\int\_0^1 W^2 dx}
$$

In Equation (23), *ω*<sup>0</sup> indicates dimensionless linear natural frequency.

#### **3. Solution Method**

This section is divided into two parts. In the first one, the forced vibration is solved using the Multiple Scales Method (MSM) [57], and in the second one, the free vibration is solved using Hamilton Approach (HA).

#### *3.1. Forced Vibration Solution*

To implement the MSM, the forced vibration equation, Equation (22), is converted to a perturbated form by introducing the following parameters

$$
\mathfrak{c} = 2\mathfrak{e} \ c\_{d\prime} \ \mathfrak{a} = \varepsilon \mathfrak{a}\_{1\prime} \ f = \mathfrak{e} \ f\_1 \tag{24}
$$

where *ε* is a dimensionless quantity that measures the strength of the nonlinearity of the beam and is called the gauge parameter.

Substituting Equation (24) into Equation (22), we obtain

$$
\ddot{q} + 2\varepsilon \, c\_d \dot{q} + \omega\_0^2 q + \varepsilon a\_1 q^3 = \varepsilon \, f\_1 \cos(\Omega t) \tag{25}
$$

In line with the MSM, the original time is replaced with new time scales as *Tn* = *ε<sup>n</sup> t*; *n* = 1, 2, . . . and therefore, the ODE is converted to a PDE.

New differential operators based on new time scales are *Dn* = *∂*/*∂Tn*, and original time first and second derivatives in terms of these operators are expressed as

$$\begin{aligned} \frac{\text{d}}{\text{d}t} &= D\_0 + \varepsilon \, D\_1 + \varepsilon^2 \, D\_2 + \dots \\ \frac{\text{d}^2}{\text{d}t^2} &= D\_0 \, ^2 + 2 \, \varepsilon \, D\_0 D\_1 + \varepsilon^2 \left( D\_1^2 + 2 \, D\_0 D\_2 \right) + \dots \end{aligned} \tag{26}$$

The governing equation includes a nonlinear cubic term. Therefore, a first-order perturbation approximation is accurate enough, such that

$$q = q\_0 + \varepsilon \, q\_1 \tag{27}$$

*qn*, *n* = 0, 1 are independent of the gauge parameter *ε*. For this reason, we can equal the coefficient of each power of *ε* to zero.

Combining Equations (25)–(27), and equating coefficients of *ε*<sup>0</sup> and *ε*<sup>1</sup> to zero, the following PDEs are attained

Coefficients of *ε*<sup>0</sup>

$$D\_0^2 q\_0 + \omega\_0^2 q\_0 = 0\tag{28}$$

Coefficients of *ε*<sup>1</sup>

$$D\_0^2 q\_1 + \omega\_0^2 q\_1 = -2D\_0 D\_1 q\_0 - 2D\_0 q\_0 - \alpha\_1 q\_0^{\;3} + f\_1 \cos(\Omega t) \tag{29}$$

The solution of Equation (28) takes the following form

$$q\_0 = A(T\_1)e^{i\
u\_0 T\_0} + \overline{A}(T\_1)e^{-i\
u\_0 T\_0} \tag{30}$$

in which *A*(*T*1) is a complex-valued function and *A*(*T*1) is its complex conjugate.

Substituting Equation (30) into Equation (29), the following equation is obtained as

$$D\_0^2 q\_1 + \omega\_0^2 q\_1 = \left[ -3\omega\_1 A^2 \overline{A} - 2\dot{\alpha}\_d A \omega\_0 - 2\dot{\omega}\omega\_0 \frac{d\mathbf{A}}{dT\_1} \right] e^{i\omega\_0 T\_0} + f\_1 \cos(\Omega t) + \text{CC} + \text{NST} \tag{31}$$

In the above equation, the terms inside the box bracket shows secular terms, CC stands for complex conjugates of previous terms, and NST is an abbreviation for terms with higher degrees of *e<sup>i</sup> <sup>ω</sup>*0*T*<sup>0</sup> (nonsecular terms).

By equating secular terms to zero, the frequency-amplitude relation can be obtained. However, the external loading can also give rise to secure terms. This fact is considered

in two states, i.e., the primary resonance and the secondary resonance. In this paper, the primary resonance is analyzed, which states that

$$
\Omega = \omega\_0 + \varepsilon \,\sigma \tag{32}
$$

Writing the trigonometric function in Equation (31) and using Equation (32), we obtain

$$3\alpha\_1 A^2 \overline{A} + 2ic\_d A \omega\_0 + 2i\omega\_0 \frac{d\mathbf{A}}{dT\_1} - \frac{1}{2} f\_1 e^{i\sigma T\_1} = 0\tag{33}$$

The complex-valued function *A* is written as

$$A = \frac{1}{2}ae^{i\theta}, \ \overline{A} = \frac{1}{2}ae^{-i\theta} \tag{34}$$

in which *a* and *θ* are the amplitude and phase, which are functions of *T*1.

Substituting Equation (34) into Equation (33) and then separating the resulting equation into real and imaginary parts yields

Imaginary parts:

$$\frac{d\mathbf{a}}{d\mathbf{T}\_1} = -ac\_d + \frac{1}{2\omega\_0} f\_1 \sin(\sigma T\_1 - \theta) \tag{35}$$

Real parts:

$$a\frac{d\theta}{dT\_1} = \frac{3}{8\omega\_0}a\_1a^3 - \frac{1}{2\omega\_0}f\_1\cos(\sigma T\_1 - \theta) \tag{36}$$

Equations (35) and (36) are converted to an autonomous equation by introducing *γ* = (*σT*<sup>1</sup> − *θ*), which results in

$$\frac{\text{da}}{\text{dT}\_1} = -ac\_d + \frac{1}{2\omega\_0} f\_1 \sin(\gamma) \tag{37}$$

$$a\frac{d\gamma}{dT\_1} = \sigma a - \frac{3}{8\omega\_0}\alpha\_1 a^3 + \frac{1}{2\omega\_0}f\_1 \cos(\gamma)\tag{38}$$

A bounded response is acquired while da dT1 <sup>=</sup> *<sup>a</sup>* <sup>d</sup>*<sup>γ</sup>* <sup>d</sup>*T*<sup>1</sup> = 0, whereby one can obtain

$$ac\_d = \frac{1}{2\omega\_0} f\_1 \sin(\gamma) \tag{39}$$

$$
\sigma a - \frac{3}{8\omega\_0} \alpha\_1 a^3 = -\frac{1}{2\omega\_0} f\_1 \cos(\gamma) \tag{40}
$$

After some mathematical manipulation and using the fact sin2(*γ*) + cos2(*γ*) = 1, we obtain the frequency-amplitude response as

$$\left[c\_d a\right]^2 + \left[\sigma a - \frac{3}{8\omega\_0} a\_1 a^3\right]^2 = \left[\frac{1}{2\omega\_0} f\_1\right]^2\tag{41}$$

#### *3.2. Free Vibration Solution*

In this subsection, the nonlinear frequency of the micro/nanobeam with neglecting the external force and damping is obtained via Hamilton's approach. The initial conditions for the vibration of the hyperelastic beam are expressed as

$$q(0) = a\_{0\prime} \,\,\dot{q}(0) = 0\tag{42}$$

where *a*<sup>0</sup> stands for the maximum amplitude of the vibration. Based on Hamilton's principle, the nonlinear frequency is derived as [58]

$$
\omega\_{nl} = \sqrt{\omega\_0^2 + \frac{49}{70} \alpha a\_0^2} \tag{43}
$$

#### **4. Result and Discussion**

The effects of several parameters such as the stiffening parameter, the length scale parameter, and forcing amplitude and damping on the frequency response and nonlinear frequency of the system are analyzed. The material and geometrical parameters of the hyperelastic microbeam are given in Table 1.



#### *4.1. Frequency Response*

Figure 2 depicts the influence of the gauge parameter *ε* on the frequency response under the following parameter *f*<sup>1</sup> = 0.5, - = 0, *cd* = 0.004, and *Jm* = 100. As the gauge parameter *ε* is decreased, the nonlinearity of the system increases. Mathematically speaking, with the decrease of *ε*, the value of nonlinear terms in the equation of motion becomes higher in comparison to the value of linear terms. Depending on the accuracy, an arbitrary value for *ε* can be adopted, which in this paper it is chosen as *ε* = 1.

**Figure 2.** Influence of gauge parameter (*ε*) on frequency response of system. Systems parameters: -= 0; *cd* = 0.004; *Jm* = 100; *f*<sup>1</sup> = 0.5.

Illustrated in Figure 3 is the influence of the damping coefficient *cd* on the frequency response of the system while considering the following parameters - = 0, and *Jm* = 100. From the figure, it is concluded that increasing the damping coefficient decreases the response amplitude of the hyperelastic micro/nanobeam. The damping in hyperelastic materials mainly originates from the viscosity of matter. In the remaining part of the numerical simulation, as a test case, the damping coefficient is taken as *cd* = 0.004.

**Figure 3.** Influence of damping coefficient (*cd*) on frequency response of system. Systems parameters: -= 0; *Jm* = 100; *f*<sup>1</sup> = 0.5.

Figure 4 represents the impact of the length scale parameter on the nonlinear resonant vibration of the hyperelastic beam. As seen, increasing the size effect, the response amplitude decreases, and the hardening nonlinearity becomes weaker. This result is in agreement with that shown in the literature for linear materials. Obtaining the accurate small-length scall parameter in the experimental test is a crucial task for engineers. Finding an exact value for the length scale parameter for the hyperelastic beam in the experimental test should be carried out to improve the performance of hyperelastic microbeams.

**Figure 4.** Influence of length scale parameter (-) on frequency response of system. Systems parameters: *Jm* = 100; *f*<sup>1</sup> = 0.5; *cd* = 0.004.

The influence of the stiffening parameter *Jm* on the frequency-amplitude plot is shown in Figure 5. It is concluded that as the stiffening parameter is decreased, the hardening nonlinearity gets stronger. When the stiffening parameter is equal to *Jm* = ∞, i.e., the conversion of the Gent model to the neo-Hookean model, the system's response is linear. It is noted that if the stiffening parameter is smaller, the strain-stiffening effect is stronger. As reported by Amabili, a stiffening parameter in a range *Jm* = 30 − 100 stands for rubber materials, and values less than them stand for biological tissues [12].

**Figure 5.** Influence of stiffening parameter (*Jm*) on frequency response of system. Systems parameters: -= 0; *f*<sup>1</sup> = 0.5; *cd* = 0.004.

The influence of the amplitude of the external loading *f*<sup>1</sup> on the resonant characteristics of the hyperelastic micro/nanobeam is analyzed in Figure 6. Increasing *f*<sup>1</sup> the response amplitude increases, and the frequency response becomes wider. Moreover, the forcing amplitude cannot alter the nonlinear nature of the system and only quantitatively alter the resonant behaviour.

**Figure 6.** Influence of forcing amplitude (*f*1) on frequency response of system. Systems parameters: -= 0; *Jm* = 100; *cd* = 0.004.

We analyze the influence of the strain-stiffening parameter on the force-response in Figure 7. The system parameters are - = 0; *σ* = 0.05; *cd* = 0.004. We can see that by increasing the value of the strain-stiffening parameter, a higher value of forcing amplitude is required to cause the jump phenomenon. Moreover, by increasing the strain-stiffening parameter, the system becomes stable and for the neo-Hookean model.

We show the influence of the length-scale parameter on the force-response in Figure 8. With the inclusion of the effect of size, the jump phenomenon arises for higher values of forcing amplitude.

**Figure 7.** Influence of stiffening parameter (*Jm*) on force response of system. Systems parameters: -= 0; σ = 0.05; *cd* = 0.004.

**Figure 8.** Influence of stiffening parameter (-) on force response of system. Systems parameters: *Jm* = 100; σ = 0.05; *cd* = 0.004.

#### *4.2. Nonlinear Frequency*

The previous section demonstrates the results for the forced vibration of the Gent hyperelastic beam. Herein, the nonlinear frequency of the system given in Equation (42) is evaluated.

Illustrated in Figure 9 is the nonlinear frequency versus the maximum amplitude when - = 0, and *Jm* = 100. It is found that by increasing the maximum amplitude *a*<sup>0</sup> the nonlinear frequency increases.

As depicted in Figure 10, the nonlinear frequency for variations of the length of the beam is presented. As the length is increased, the dimensionless nonlinear frequency increases accordingly.

The nonlinear frequency versus the stiffening parameter *Jm* is presented in Figure 11. Increasing *Jm*, the nonlinear frequency decreases.

As depicted in Figure 12, the nonlinear frequency versus the length scale parameter is presented. As the size effect is increased, the nonlinear frequency increases accordingly.

**Figure 9.** Influence of maximum amplitude (*a*0) on nonlinear frequency of system. Systems parameters: -= 0; *Jm* = 100, *L* = 30 μm.

**Figure 10.** Influence of length of micro/nanobeam (*L*) on nonlinear frequency of system. Systems parameters: -= 0; *Jm* = 100; *a*<sup>0</sup> = 0.5.

**Figure 11.** Influence of stiffening parameter (*Jm*) on nonlinear frequency of system. Systems parameters: *L* = 30 μm; -= 0; *a*<sup>0</sup> = 0.5.

**Figure 12.** Influence of length scale parameter on nonlinear frequency of system. Systems parameters: *L* = 30 μm; *Jm* = 100; *a*<sup>0</sup> = 0.5.

#### **5. Discussion on the Strain-Stiffening**

Rubber-like materials can be deformed by stretching. In the beginning, we can stretch rubbers easily, but if the stretch is large enough, the stretching process becomes difficult. This is due to the strain-stiffening effect in rubber-like materials. The strain-stiffening is a nonlinear behavior that is seen even in soft biological materials such as liver and brain tissue [59]. We can use this property in hyperelastic materials so as to evade damage. The strain-stiffening can also be connected to the molecular-statistical point of view in nonlinear elasticity. The stiffening parameters *Jm* in the Gent model relates to the number of rigid links in a single chain *N* using *Jm* = 3(*N* − 1). *N* is also called the classical number of Kuhn segments [60]. The results of Figures 5, 7 and 11, can also be interpreted based on molecular-statistical point of view. We see that altering *Jm*, the number of segments changes accordingly. Therefore, this change affects the frequency/force response of the hyperelastic microbeam. Taken together, the results of this paper can help researchers who would like to analyze the hyperelastic microbeam via molecular-statistical hyperelastic models such as generalized neo-Hookean model.

#### **6. Conclusions**

In this paper, nonlinear, free, and forced oscillations of a hyperelastic micro/nanobeam were investigated with the inclusion of the small-scale effect, strain-stiffening effect, and moderate rotation. A developed Euler–Bernoulli beam theory was utilized to model the beam, and the energies and works that appeared in the system were formulated. The equation of motion was derived using Hamilton's principle and the Galerkin decomposition method. Frequency-amplitude curves and the nonlinear natural frequency diagrams were illustrated by analytically solving the equation of motion. This paper concludes that:


**Author Contributions:** Conceptualization, A.A. and S.D.; methodology, A.A.; software, A.A.; validation, A.A. and S.D.; formal analysis, A.A.; investigation, M.M. and V.A.E.; resources, A.A. and M.M.; data curation, A.A.; writing—original draft preparation, A.A. and S.D.; writing—review and editing, A.A., S.D., M.M. and V.A.E.; visualization, A.A., S.D.; supervision, V.A.E.; project administra-

tion, M.M. and V.A.E.; funding acquisition, S.D. and M.M. All authors have read and agreed to the published version of the manuscript.

**Funding:** V.A. Eremeyev acknowledges the financial support of the Ministry of Education and Science of the Russian Federation as part of the program of the Mathematical Center for Fundamental and Applied Mathematics under the agreement №075-15-2019-1621.

**Data Availability Statement:** Data sharing is not applicable.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


## *Article* **Hyperelastic Microcantilever AFM: Efficient Detection Mechanism Based on Principal Parametric Resonance**

**Amin Alibakhshi <sup>1</sup> , Sasan Rahmanian 2, Shahriar Dastjerdi <sup>3</sup> , Mohammad Malikan 4,\* , Behrouz Karami <sup>5</sup> , Bekir Akgöz <sup>3</sup> and Ömer Civalek 3,6**


**Abstract:** The impetus of writing this paper is to propose an efficient detection mechanism to scan the surface profile of a micro-sample using cantilever-based atomic force microscopy (AFM), operating in non-contact mode. In order to implement this scheme, the principal parametric resonance characteristics of the resonator are employed, benefiting from the bifurcation-based sensing mechanism. It is assumed that the microcantilever is made from a hyperelastic material, providing large deformation under small excitation amplitude. A nonlinear strain energy function is proposed to capture the elastic energy stored in the flexible component of the device. The tip–sample interaction is modeled based on the van der Waals non-contact force. The nonlinear equation governing the AFM's dynamics is established using the extended Hamilton's principle, obeying the Euler–Bernoulli beam theory. As a result, the vibration behavior of the system is introduced by a nonlinear equation having a time-dependent boundary condition. To capture the steady-state numerical response of the system, a developed Galerkin method is utilized to discretize the partial differential equation to a set of nonlinear ordinary differential equations (ODE) that are solved by the combination of shooting and arc-length continuation method. The output reveals that while the resonator is set to be operating near twice the fundamental natural frequency, the response amplitude undergoes a significant drop to the trivial stable branch as the sample's profile experiences depression in the order of the picometer. According to the performed sensitivity analysis, the proposed working principle based on principal parametric resonance is recommended to design AFMs with ultra-high detection resolution for surface profile scanning.

**Keywords:** atomic force microscopy; hyperelastic microcantilever; softening resonance; non-contact cantilever; shooting and arc-length continuation method; developed Galerkin method

#### **1. Introduction**

Atomic force microscopy (AFM) is a device used to study materials' properties and surface structure in micro- and nanometer dimensions [1–3]. Flexibility, multiple potential signals, and the ability to operate the device in different modes have enabled researchers to study AFMs at different levels and under different environmental conditions. This device makes it possible to examine conductive or insulating surfaces, soft or hard, cohesive or

**Citation:** Alibakhshi, A.; Rahmanian, S.; Dastjerdi, S.; Malikan, M.; Karami, B.; Akgöz, B.; Civalek, Ö. Hyperelastic Microcantilever AFM: Efficient Detection Mechanism Based on Principal Parametric Resonance. *Nanomaterials* **2022**, *12*, 2598. https:// doi.org/10.3390/nano12152598

Academic Editor: James Murray Hill

Received: 9 June 2022 Accepted: 26 July 2022 Published: 28 July 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

powdery, biological and organic or inorganic. The general architecture of an AFM includes a microcantilever with a sharp tip at the free end contacting with the surface of samples [4]. The AFMs are applicable devices that have been widely used in so many scientific purposes, such as chemistry [5], materials science [6], medical and biophysics [7], and biological and colloidal applications [8].

Such as many structures that are influenced by the time [9–12], the operation of AFMs may be time-dependent and, consequently, under dynamic conditions. Thus, most previous studies in this field have focused on the dynamic performance of AFMs. The chaotic motion and bifurcation of a tapping mode AFM were examined by Yagasaki [13]. Bahrami and Nayfeh [14] explored the nonlinear oscillation of a tapping mode AFM by using a highdimensional Galerkin discretization technique and a four-step Adams–Moulton method. They modeled the microcantilever based on the Euler–Bernoulli beam theory. The flexural and torsional vibrations of an AFM were investigated by Kahrobaiyan et al. [15], who performed a sensitivity analysis. Arafat and co-investigators [16] analyzed the frequency–amplitude response of a contact mode AFM analytically using the multiple time-scale method. Dastjerdi and Abbasi [17] studied the free vibration of a cracked AFM, incorporating the influence of crack size and its location. They utilized the transfer matrix method to solve the cracked system. Also, they concluded that existing a crack (with a specific size and location) on the AFM cantilever can be beneficial for controlling some unwanted phenomena. Mahmoudi et al. [18] studied the resonance of a non-contact AFM using harmonic balance and multiple time-scale methods. Saeidi et al. [19] explored forced vibrations of an AFM, taking the temperature and contact effects into account. Ahmadi and co-workers [20] studied free and forced oscillations of AFM with rectangular and V-shaped and employed the finite element analysis. Kouchaksaraei and Bahrami [21] proposed a new multifrequency excitation for AFM in a non-contact regime. They analyzed the sensitivity to the Hamaker constant and the initial tip–sample distance to enhance the performance of the AFM. In most previous studies investigating dynamic characteristics of AFMs, the structural materials have been a linear type. This means that the linear elasticity governs their behavior with small strain and deformation.

There has been an increasing interest in soft and hyperelastic materials in recent decades, because they are lightweight, cheap, compatible with many flexible structures, and show ease of fabrication. The relation between the strain and stress for hyperelastic materials is nonlinear. Besides, the strain in hyperelastic materials is moderate or large, and they undergo large deformation. These are the main differences between hyperelastic materials and linear elastic materials. Hyperelastic materials are structural materials of many electromechanical systems such as actuators, sensors, and energy harvesters. More specifically, hyperelastic materials have been widely used for electro-active-based systems [22–25]. Up to now, hyperelastic materials have been designed in different geometries, e.g., beams [26–29], plates [30,31], and shells [32,33]. Reviewing the previous studies shows a growing body of literature on using hyperelastic materials in different systems. Due to these advantages of hyperelastic materials and their broad applications, this paper proposes an AFM made of such materials, unique in simulation and application. Most recently, the dielectric elastomer actuators have been proposed [34] for AFM, and because these kinds of actuators are hyperelastic, exploring hyperelastic-based AFM may be of significant importance.

As is reported in the literature review, it can be seen that there is a lack of nonlinear study on vibration analysis of AFM cantilevers made of hyperelastic material. Consequently, thorough research on nonlinear frequency analysis of AFM cantilevers has been conducted based on the hyperelasticity approach. This paper is structured as follows. Section 2 derives the equation of motion based on a mathematical model containing material and geometrical nonlinearities, size effects, and finite rotations. In Section 3, the governing equation is discretized with a developed Galerkin decomposition method and solved through the shooting method combined with the arc-length continuation method. Finally, in Section 4, numerical results have been provided to study different dynamical behaviors of the AFM cantilever under the assumed conditions.

#### **2. Mathematical Modeling**

Figure 1 demonstrates the schematic view of an AFM cantilever simulated in this paper. The system consists of a hyperelastic microcantilever with a sharp tip at the free end. A Cartesian coordinate system (*x*, *y*, *z*) is employed to define the configuration of the AFM. The base vectors in *x*, *y*, *z*−directions are *e*1, *e*2, and *e*3, respectively. The microcantilever's length, thickness, and width are denoted via *L*, *d*, and *b*, respectively.

**Figure 1.** The schematic view of an AFM's cantilever scanning a sample. (**a**) 3D view of the AFM. (**b**) 2D view of the cantilever.

Figure 1a is a 3D schematic and perhaps does not provide more details. On the other hand, Figure 1b gives more details. The AFM consists of a microcantilever mounting on a base with piezoelectric patches as the excitation source (see Figure 1b). As depicted in Figure 1b, *Z*ˆ is the distance between the tip and the sample when the microbeam is located at rest, and *w*(*x*, *t*) is the transverse displacement of the microcantilever with respect to the base frame. The transverse base excitation is *Zb*(*t*) = *d*<sup>0</sup> sin(*ωt*), in which *d*<sup>0</sup> and *ω* are the amplitude and frequency, respectively. The total deflection is *u*(*x*, *t*) = *w*(*x*, *t*) + *Zb*(*t*). Regarding these parameters, the instantaneous tip–sample separation is indicated by *<sup>z</sup>*ˆ(*t*) = *<sup>Z</sup>*<sup>ˆ</sup> <sup>−</sup> *<sup>w</sup>*(*L*, *<sup>t</sup>*) <sup>−</sup> *Zb*(*t*).

The governing equations are derived under the following assumptions: (1) the effect of rotary inertia and shear deformation is disregarded, i.e., the microbeam is modeled based on the Euler-Bernoulli beam theory. (2) both geometric and material nonlinearities are considered. (3) the rotation of the microcantilever is considered to be a moderate type, which means that *<sup>∂</sup>w*/*∂<sup>x</sup>* is moderate. In other words, the rotation is of order <sup>√</sup>, in which 1 stands for a small parameter (for more details on deriving moderate rotation see [35]). (4) the size effect is considered using the modified couple stress theory. (5) the axial displacement is neglected.

In reality, a sample subjected to probing may have an inhomogeneous surface, thus leading to surface depression or an increase in surface height. During probing, these increases in surface height may affect the performance of the hyperelastic AFM. As shown in Figure 2, the location of surface depression and increase in surface height assumed in the system are shown. In Figure 2, *δ* refers to any variation in the height of the sample profile with respect to the baseline, i.e., *Z*ˆ + *δ* = *Z*ˆ + 0 = *Z*ˆ.

**Figure 2.** Schematic representation of the increase and decrease in surface profile.

#### *2.1. Beam Theory with Finite Rotation and Deformation*

The components of the displacement vector *u* = *u*1*e*ˆ1 + *u*2*e*ˆ2 + *u*3*e*ˆ3 can be given as [36].

$$\begin{cases} \boldsymbol{u}\_1 = z \boldsymbol{\theta}\_{\mathcal{X}} & \boldsymbol{u}\_2 = \mathbf{0} \\ \boldsymbol{u}\_3 = \boldsymbol{w}(\mathbf{x}, t) \\ \boldsymbol{\theta}\_{\mathcal{X}} \equiv -\frac{\partial \boldsymbol{w}}{\partial \mathbf{x}} \end{cases} \tag{1}$$

In Equation (1), *θ<sup>x</sup>* is a defined parameter, namely, the minus of the slope of the transverse displacement.

The strain–displacement relations for the system are obtained according to the finite deformation. To this end, the von Kármán theory originated from Green's strain theory is used as follows [37]

$$\begin{cases} E\_{xx} = \varepsilon\_{xx}{}^{(0)} + z \varepsilon\_{xx}{}^{(1)} \\\ E\_{zz} = \varepsilon\_{zz}{}^{(0)} \end{cases} \tag{2}$$

in which

$$\begin{cases} \begin{aligned} \varepsilon\_{xx}^{(0)} &= \frac{1}{2} \left( \frac{\partial w}{\partial x} \right)^{2} \\ \varepsilon\_{xx}^{(1)} &= \frac{\partial \theta\_{\mathbf{x}}}{\partial x} \\ \varepsilon\_{zz}^{(0)} &= \frac{1}{2} \theta\_{\mathbf{x}} \end{aligned} \tag{3}$$

In Equation (2), *Exx* is the strain in the *x*-direction, and *Ezz* stands for the strain in the *z*-direction.

The potential energy is expressed in terms of the strain energy function for hyperelastic materials. The strain energy can be formulated based on the displacement as mentioned above and strain components incorporating the finite deformation, finite rotation, and size effect as follows [38,39].

$$\Psi = \frac{1}{2} \left( a\_1 E\_{\text{xx}}^2 + a\_2 E\_{zz}^2 + a\_3 \left( \theta\_{\text{x}}^{'} \right)^2 + 2a\_4 E\_{\text{xx}} E\_{zz} \right) \tag{4}$$

in which

$$\begin{cases} a\_1 = a\_2 = 2\mu + \lambda\\ a\_4 = \lambda\\ a\_3 = 2\mu \ell^2 \end{cases} \tag{5}$$

In Equation (5), is the internal length-scale parameter capturing the size effect, which is defined as the square root of the ratio of the moduli of curvature to the shear [38,40] and is generally quantified as one-half or one-fourth of the thickness of the structure in theoretical analyses. It is noted that the value of this parameter from experimental evidence for hyperelastic microstructures seems to be rare. For more details on experimental methods of calculating -, you may refer to [41].

Moreover, in Equation (15), *μ* = *E*/2(1 + *ν*) stands for the shear modulus, and *λ* = [*Eν*]/[(1 + *ν*)(1 − 2*ν*)] is the Lamé's constant (*ν* is the Poisson's ratio). It is mentioned that Equation (5) has shown its applicability for hyperelastic structures in previous literature. The hyperelastic potential energy is therefore calculated as:

> *Uel* = Ψ*dV* (6)

where *V* is the volume of the microcantilever.

Substituting Equation (2) into Equation (4), then inserting it into Equation (6), the final form of the hyperelastic energy is presented as:

*V*

$$\begin{split} \mathcal{U}\_{el} = \int\_{0}^{L} \left[ \frac{1}{8} a\_1 A \left( \frac{\partial w}{\partial \mathbf{x}} \right)^4 + \frac{1}{2} a\_1 I \left( \frac{\partial^2 w}{\partial \mathbf{x}^2} \right)^2 + \frac{1}{8} a\_2 A \left( \frac{\partial w}{\partial \mathbf{x}} \right)^4 + \frac{1}{2} a\_3 A \left( \frac{\partial^2 w}{\partial \mathbf{x}^2} \right)^2 \\ \quad + \frac{1}{4} a\_4 A \left( \frac{\partial w}{\partial \mathbf{x}} \right)^4 \right] d\mathbf{x} \end{split} \tag{7}$$

With the cross-section area *A* = *bd* and the second moment of area *I* = *bd*3/12.

#### *2.2. Tip–Sample Interaction*

At the free end of the microcantilever, there is the interaction force between the tip and sample. In the present paper, as a test case, a non-contact tip–sample interaction so-called the van der Waals non-contact force is applied, i.e., [14],

$$F\_{\rm vdW} = -\frac{HR}{6\mathcal{E}^2} \tag{8}$$

in which *F*vdW stands for the van der Waals non-contact force, *H* is the Hamaker constant, and *R* is the radius of the spherical tip apex.

The system's potential energy due to the tip–sample interaction is *V*tp = − *F*vdW *dz*ˆ. By calculating this integration, the final form of the tip–sample interaction potential energy for the non-contact region will be formulated as [14]:

$$dL\_{\rm vdW} = -\frac{HR}{6\left[\hat{Z} - w(L\_\prime t) - Z\_b(t)\right]}\tag{9}$$

#### *2.3. Kinetic Energy*

The total kinetic energy of the AFM cantilever can be expressed below (*ρ* is the mass density of the microcantilever).

$$d\mathcal{U}\_k = \frac{1}{2} \int\_0^L \rho A \left[ \left( \frac{\partial w}{\partial t} \right) + \left( \frac{\partial Z\_b(t)}{\partial t} \right) \right]^2 dx \tag{10}$$

#### *2.4. Surrounding Damping Force*

The constitutive material of the AFM is hyperelastic, for example, it is rubbery or elastomeric. For such materials, viscoelasticity plays a major role that should be considered in the analysis. In this paper, a linear damping model is assumed, while other complicated types of damping can also be used.

The amount of work created by the viscous surrounding medium is formulated as (*cd* is the viscous damping coefficient).

$$W\_D = -c\_d \int\_0^L \left[ \left(\frac{\partial w}{\partial t}\right) + \left(\frac{\partial Z\_b(t)}{\partial t}\right) \right] w dx \tag{11}$$

#### *2.5. Hamilton's Principle and Equation of Motion*

The governing equation can be easily derived using variational approaches by obtaining the energies and works that appear in the system [42–44]. The partial differential equation governing the motion of the AFM is derived using the following Hamilton's principle as one of the variational approaches:

$$\int\_{t\_1}^{t\_2} [\mathcal{U}\_k - \mathcal{U}\_s] dt + \int\_{t\_1}^{t\_2} \delta \mathcal{W}\_D dt = 0 \tag{12}$$

where *Us* shows the total potential energy of the system, which for the non-contact region is: *Us* = *Uel* + *U*vdW.

Eventually, substituting the expression of the kinetic and potential energies and the work of damping force into Hamilton's principle gives the following formulation:

$$\begin{aligned} \rho A \frac{\partial^2 w}{\partial t^2} + \rho A & \quad \frac{\partial^2 Z\_b(t)}{\partial t^2} + c\_d \frac{\partial w}{\partial t} + c\_d \frac{\partial Z\_b(t)}{\partial t} - \frac{3}{2} a\_1 A \left(\frac{\partial^2 w}{\partial x^2}\right) \left(\frac{\partial w}{\partial x}\right)^2 \\ &+ a\_1 I \left(\frac{\partial^4 w}{\partial x^4}\right) - \frac{3}{2} a\_2 A \left(\frac{\partial^2 w}{\partial x^2}\right) \left(\frac{\partial w}{\partial x}\right)^2 + a\_3 A \left(\frac{\partial^4 w}{\partial x^4}\right) \\ &- 3 a\_4 A \left(\frac{\partial^2 w}{\partial x^2}\right) \left(\frac{\partial w}{\partial x}\right)^2 = 0 \end{aligned} \tag{13}$$

It is noteworthy that the mathematical definition of the boundary conditions is obtained from Hamilton's principle. Consequently, the boundary condition definition can be formulated as the following equations.

$$\begin{cases} w(0,t) = \frac{\partial w}{\partial \mathbf{x}}(0,t) = 0 \\ \frac{\partial^2 w}{\partial \mathbf{x}^2}(L,t) = 0 \\ a\_1 I \frac{\partial^3 w}{\partial \mathbf{x}^3}(L,t) + a\_3 A \frac{\partial^3 w}{\partial \mathbf{x}^3}(L,t) = -\frac{HR}{6 \left[2 - w(L,t) - g(t)\right]^2} \end{cases} \tag{14}$$

The boundary conditions for cantilevers state that the deflection and slope at the point *x* = 0 (at the fixed end) are equal to zero, and the bending moment and shear force at *x* = *L* (at the free end) are equal to zero. From Equation (14), it is observed that the difference between the boundary conditions of the cantilever and AFM appears in the shear force in such a way that for the AFM it is not equal to zero any longer and is nonhomogeneous and time-dependent due to tip–sample interaction.

#### *2.6. Nondimensionalization*

In this section, the equation of motion and boundary conditions are made dimensionless. Hence, the dimensionless parameters can be presented as follows [45]:

$$\begin{aligned} \mathbf{x}^\* &= \frac{\mathbf{x}}{\mathbf{L}}, w^\* = \frac{w}{\mathbf{Z}}, \overline{d}\_0 = \frac{d\_0}{\mathbf{Z}}, t^\* = t\sqrt{\frac{EI}{\rho AL^4}}, c = \frac{c\_d L^4}{EI}\sqrt{\frac{EI}{\rho AL^4}}, \Omega = \omega\sqrt{\frac{\rho AL^4}{EI}}\\ \mathbf{d}\_1 &= \frac{a\_1 I}{EI}, \mathbf{d}\_2 = \frac{a\_3 A}{EI} = \frac{2\mu A \hat{\mathbf{Z}}^2}{EI}, \mathbf{d}\_3 = \frac{3a\_1 A \hat{\mathbf{Z}}^2}{2EI}, \mathbf{d}\_4 = \frac{3a\_2 A \hat{\mathbf{Z}}^2}{2EI}, \mathbf{d}\_5 = \frac{3a\_4 A \hat{\mathbf{Z}}^2}{EI} \end{aligned} \tag{15}$$

The non-dimensional form of the equation of motion and boundary conditions are obtained as (the asterisk notation is disregarded for simplification):

$$\begin{aligned} \frac{\partial^2 w}{\partial t^2} + c \frac{\partial w}{\partial t} - \overline{d}\_0 \, \Omega^2 \sin(\Omega t) + c \overline{d}\_0 \, \Omega \cos(\Omega t) + \underbrace{(\mathbf{d}\_1 + \mathbf{d}\_2) \frac{\partial^4 w}{\partial x^4}}\_{=0} \\ - \underbrace{(\mathbf{d}\_3 + \mathbf{d}\_4 + \mathbf{d}\_5) \left(\frac{\partial^2 w}{\partial x^2}\right) \left(\frac{\partial w}{\partial x}\right)^2}\_{\beta} = 0 \end{aligned} \tag{16}$$

with the dimensionless boundary conditions for the non-contact region

$$\begin{aligned} w(0, t) &= \frac{\partial w}{\partial \mathbf{x}}(0, t) = \frac{\partial^2 w}{\partial \mathbf{x}^2}(1, t) = 0\\ (\underbrace{\mathbf{d}\_1 + \mathbf{d}\_2}\_{\alpha}) \frac{\partial^3 w}{\partial \mathbf{x}^3}(1, t) &= -\frac{\mathbf{d}\_6}{\left[ (\overline{\mathbf{z}}) - w(1, t) - \overline{d}\_0 \sin(\Omega t) \right]^2} \end{aligned} \tag{17}$$

#### **3. Discretization of the Governing Motion's Equation**

In this section, the governing partial differential equation is discretized by using the Galerkin method. Because the boundary condition is time-dependent, implementing the Galerkin method directly may be difficult. To solve such non-homogeneous boundary conditions, the following process is implemented.

It is assumed that the deflection has the following form [46]

$$w(\mathbf{x},t) = \mathcal{W}(\mathbf{x},t) + F(t) \; G(\mathbf{x}) \tag{18}$$

in which *F*(*t*) is obtained from non-homogeneous boundary conditions. *F*(*t*) can be expressed as the following formulation.

$$F(t) = -\frac{\mathbf{d}\_6}{a\left[\overline{Z} - w(1, t) - \overline{d}\_0 \sin(\Omega t)\right]^2} \tag{19}$$

In Equation (18), *G*(*x*) is an arbitrary function satisfying the following conditions:

$$G(0) = \frac{dG}{d\mathbf{x}}(0) = G(1) = \frac{dG}{d\mathbf{x}}(1) = \frac{d^2G}{d\mathbf{x}^2}(1) = 0,\\ \frac{d^3G}{d\mathbf{x}^3}(1) = 1\tag{20}$$

A suitable assumption for the function *G*(*x*) is chosen as

$$G(\mathbf{x}) = \frac{-1}{6}\mathbf{x}^2 + \frac{1}{2}\mathbf{x}^3 - \frac{1}{2}\mathbf{x}^4 + \frac{1}{6}\mathbf{x}^5 \tag{21}$$

Substituting Equation (18) into Equations (16) and (17), the governing equation is derived as:

$$\begin{array}{ll} \frac{\partial^2 \mathbf{W}}{\partial t^2} + c \frac{\partial \mathbf{W}}{\partial t} & + a \frac{\partial^4 \mathbf{W}}{\partial x^4} - \beta \left(\frac{\partial^2 \mathbf{W}}{\partial x^2}\right) \left(\frac{\partial \mathbf{W}}{\partial x}\right)^2 + G \frac{d^2 F}{dt^2} + cG \frac{dF}{dt} + aF \frac{d^4 G}{dx^4} \\ & - 2\beta F \frac{dG}{dx} \left(\frac{\partial^2 \mathbf{W}}{\partial x^2}\right) \left(\frac{\partial \mathbf{W}}{\partial x}\right) - \beta F^2 \left(\frac{dG}{dx}\right)^2 \left(\frac{\partial^2 \mathbf{W}}{\partial x^2}\right) \\ & - \beta F \left(\frac{d^2 G}{dx}\right) \left(\frac{\partial \mathbf{W}}{\partial x}\right)^2 - 2\beta F^2 \left(\frac{dG}{dx}\right) \left(\frac{d^2 G}{dx}\right) \left(\frac{\partial \mathbf{W}}{\partial x}\right) \\ & - \beta F^3 \left(\frac{dG}{dx}\right)^2 \left(\frac{d^2 G}{dx}\right) = \overline{d}\_0 \, \Omega^2 \sin(\Omega t) - c \overline{d}\_0 \, \Omega \cos(\Omega t) \end{array} \tag{22}$$

Also, the corresponding boundary condition for Equation (22) is as follows:

$$W(0,t) = W'(0,t) = W''(1,t) = W''''(1,t) = 0\tag{23}$$

The prime notation is derivative with respect to the dimensionless axial coordinate.

It is seen that the non-homogeneous boundary condition has been transformed into a homogeneous one. Now, the Galerkin method is applied to Equations (22) and (23). Since the contribution of the first mode is dominant in comparison to higher modes, the first mode is adopted in the present work.

By applying the variables separation technique, the response of the new variable *U*(*x*, *t*) is supposed to be as:

$$\mathcal{W}(\mathfrak{x},t) = q(t)\phi(\mathfrak{x})\tag{24}$$

in which *q*(*t*) denotes the generalized coordinate, and *φ*(*x*) stands for the eigenfunction for a clamped-free beam, given by

$$\begin{split} \phi(\mathbf{x}) &= \left\{ \cos(\beta\_1 \mathbf{x}) - \cos \mathbf{h}(\beta\_1 \mathbf{x}) \right\} \\ &- \left( \frac{\{\cos(\beta\_1) + \cos \mathbf{h}(\beta\_1)\}}{\{\sin(\beta\_1) + \sin \mathbf{h}(\beta\_1)\}} \left\{ \sin(\beta\_1 \mathbf{x}) - \sin \mathbf{h}(\beta\_1 \mathbf{x}) \right\} \right) \end{split} \tag{25}$$

where *β*<sup>1</sup> = 1.8751.

Substituting Equation (24) into Equation (22) and multiplying the resultant equation by *φ*(*x*), then integrating with respect to *x* from 0 to 1, the following governing equation in the time domain is derived as:

$$\begin{aligned} \ddot{q} + c\dot{q} + \omega\_0^2 q - r\_2 q^3 + r\_3 \ddot{F} + cr\_3 \dot{F} + r\_4 F - r\_5 F q^2 - r\_6 F^2 q - r\_7 F q^2 \\ -r\_8 F^2 q - r\_9 F^3 = r\_{10} \,\Omega^2 \sin(\Omega t) - cr\_{10} \,\Omega \cos(\Omega t) \end{aligned} \tag{26}$$

in which

$$\begin{array}{ll} r\_1 = \int\_0^1 \alpha \,\phi \phi^{\prime\prime\prime} dx = \omega\_0^2, & r\_2 = \int\_0^1 \beta \phi \phi^{\prime\prime} (\phi^{\prime})^2 dx \\ r\_3 = \int\_0^1 \phi \,G dx, & r\_4 = \int\_0^1 \alpha \,\phi \,G^{\prime\prime\prime} dx \\ r\_5 = \int\_0^1 2\beta \phi G^{\prime} \phi^{\prime} \phi^{\prime\prime} dx, & r\_6 = \int\_0^1 \beta \phi G^{\prime} G^{\prime} \phi^{\prime\prime} dx \\ r\_7 = \int\_0^1 \beta \phi \phi^{\prime} \phi^{\prime} G^{\prime} dx, & r\_8 = \int\_0^1 2\beta \phi \phi^{\prime} G^{\prime} G^{\prime} dx \\ r\_9 = \int\_0^1 \beta \phi G^{\prime} G^{\prime} G^{\prime} dx, & r\_{10} = \int\_0^1 \overline{d}\_0 \phi \, dx \end{array} \tag{27}$$

We now simplify the equation of motion by expanding the function *F*(*t*) around zero point, i.e.,

$$F(t) = k\_1 + k\_2 q(t) \tag{28}$$

in which

$$\begin{cases} k\_1 = -\frac{h\_1}{\left(\mathbb{Z} - \mathbf{g}(t)\right)^2}, \; k\_2 = -\frac{h\_2}{\left(\mathbb{Z} - \mathbf{g}(t)\right)^3}, \; \mathbf{g}(t) = \overline{d}\_0 \sin(\Omega t), \; h\_1 = \frac{d\_0}{a},\\ h\_2 = \frac{2d\_0\phi(1)}{a} \end{cases} \tag{29}$$

Substituting Equation (28) into Equation (26), we get

$$\dot{\Omega}M\ddot{q} + \mathcal{C}\dot{q} + \mathcal{K}\_l q + \mathcal{K}\_q q^2 + \mathcal{K}\_c q^3 = r\_{10}\,\Omega^2 \sin(\Omega t) - cr\_{10}\,\Omega\cos(\Omega t) \tag{30}$$

in which

*<sup>M</sup>* <sup>=</sup> <sup>1</sup> <sup>−</sup> *<sup>h</sup>*2*r*<sup>3</sup> (*Z*−*g*(*t*)) 3 *<sup>C</sup>* <sup>=</sup> *<sup>c</sup>* <sup>−</sup> *c h*2*r*<sup>3</sup> (*Z*−*g*(*t*)) 3 *KL* <sup>=</sup> *<sup>r</sup>*<sup>1</sup> <sup>+</sup> <sup>3</sup>*h*<sup>2</sup> <sup>1</sup>*h*2*r*<sup>9</sup> (*Z*−*g*(*t*)) <sup>7</sup> <sup>−</sup> *<sup>h</sup>*<sup>2</sup> 1*r*6 (*Z*−*g*(*t*)) <sup>4</sup> <sup>−</sup> *<sup>h</sup>*<sup>2</sup> 1*r*8 (*Z*−*g*(*t*)) <sup>4</sup> <sup>−</sup> *<sup>h</sup>*2*r*<sup>4</sup> (*Z*−*g*(*t*)) 3 *Kq* <sup>=</sup> <sup>3</sup>*h*1*h*<sup>2</sup> 2*r*9 (*Z*−*g*(*t*)) <sup>8</sup> <sup>−</sup> <sup>2</sup>*h*1*h*2*r*<sup>6</sup> (*Z*−*g*(*t*)) <sup>5</sup> <sup>−</sup> <sup>2</sup>*h*1*h*2*r*<sup>8</sup> (*Z*−*g*(*t*)) <sup>5</sup> + *<sup>h</sup>*1*r*<sup>5</sup> (*Z*−*g*(*t*)) <sup>2</sup> + *<sup>h</sup>*1*r*<sup>7</sup> (*Z*−*g*(*t*)) 2 *Kc* <sup>=</sup> <sup>−</sup>*r*<sup>2</sup> <sup>+</sup> *<sup>h</sup>*<sup>3</sup> 2*r*9 (*Z*−*g*(*t*)) <sup>9</sup> <sup>−</sup> *<sup>h</sup>*<sup>2</sup> 2*r*6 (*Z*−*g*(*t*)) <sup>6</sup> <sup>−</sup> *<sup>h</sup>*<sup>2</sup> 2*r*8 (*Z*−*g*[*t*]) <sup>6</sup> + *<sup>h</sup>*2*r*<sup>5</sup> (*Z*−*g*[*t*]) <sup>3</sup> + *<sup>h</sup>*2*r*<sup>7</sup> (*Z*−*g*(*t*)) 3 (31)

To distinguish between parametric and non-parametric resonances, the resultant equation is multiplied by *Z* − *g*(*t*) 9 , resulting in

$$\begin{aligned} \left(\mathcal{M}\left(\overline{\mathcal{Z}} - \operatorname{g}(t)\right)^{9}\overline{\boldsymbol{q}}\right) &+ \mathcal{C}\left(\overline{\mathcal{Z}} - \operatorname{g}(t)\right)^{9}\dot{\boldsymbol{q}} + \mathcal{K}\_{l}\left(\overline{\mathcal{Z}} - \operatorname{g}(t)\right)^{9}\boldsymbol{q} + \mathcal{K}\_{\boldsymbol{q}}\left(\overline{\mathcal{Z}} - \operatorname{g}(t)\right)^{9}\boldsymbol{q}^{2} \\ &+ \mathcal{K}\_{\boldsymbol{c}}\left(\overline{\mathcal{Z}} - \operatorname{g}(t)\right)^{9}\boldsymbol{q}^{3} \\ &= r\_{10}\,\Omega^{2}\sin(\Omega t)\left(\overline{\mathcal{Z}} - \operatorname{g}(t)\right)^{9} - cr\_{10}\,\Omega\cos(\Omega t)\left(\overline{\mathcal{Z}} - \operatorname{g}(t)\right)^{9} \end{aligned} \tag{32}$$

Simplifying the equations by assuming that *d*<sup>0</sup> is too small compared to *Z* = 1.

$$\mathbf{M}\_{\mathbf{1}}\ddot{\boldsymbol{q}} + \mathbf{C}\_{\mathbf{1}}\dot{\boldsymbol{q}} + \mathbf{K}\_{\mathbf{L}\mathbf{1}}\boldsymbol{q} + \mathbf{K}\_{\mathbf{q}\mathbf{1}}\boldsymbol{q}^{2} + \mathbf{K}\_{\mathbf{c}\mathbf{1}}\boldsymbol{q}^{3} = r\_{10}\,\Omega^{2}\sin(\Omega t) - cr\_{10}\,\Omega\cos(\Omega t) \tag{33}$$

in which

$$\begin{aligned} \mathbf{M\_1} &= 1 - h\_2 r\_3 + \left(6\overline{d}\_0 h\_2 r\_3 - 9\overline{d}\_0\right)\sin(\Omega t) \\ \mathbf{C\_1} &= c - c \, h\_2 r\_3 + \left(6\overline{d}\_0 c \, h\_2 r\_3 - 9\overline{d}\_0 c\right)\sin(\Omega t) \\ \mathbf{K\_{L1}} &= 3h\_1^2 h\_2 r\_9 - h\_1^2 r\_6 - h\_1^2 r\_8 - h\_2 r\_4 + r\_1 \\ &+ \left(-6\overline{d}\_0 h\_1^2 h\_2 r\_9 + 5\overline{d}\_0 h\_1^2 r\_6 + 5\overline{d}\_0 h\_1^2 r\_8 + 6\overline{d}\_0 h\_2 r\_4 - 9\overline{d}\_0 r\_1\right)\sin(\Omega t) \\ \mathbf{K\_{q1}} &= 3h\_1 h\_2^2 r\_9 - 2h\_1 h\_2 r\_6 - 2h\_1 h\_2 r\_8 + h\_1 r\_5 + h\_1 r\_7 \\ &+ \left(-3\overline{d}\_0 h\_1 h\_2^2 r\_9 + 8\overline{d}\_0 h\_1 h\_2 r\_6 + 6\overline{d}\_0 h\_1 h\_2 r\_8 - 7\overline{d}\_0 h\_1 r\_5 - 7\overline{d}\_0 h\_1 r\_7\right)\sin(\Omega t) \\ \mathbf{K\_{c1}} &= -r\_2 + h\_2 r\_5 + h\_2 r\_7 - h\_2^2 r\_6 - h\_2^2 r\_8 + h\_2^3 r\_9 \\ &+ \left(3\overline{d}\_0 h\_2^2 r\_6 + 3\overline{d}\_0 h\_2^2 r\_8 - 6\overline{d}\_0 h\_2 r\_5 - 6\overline{d}\_0 h\_2 r\_7 + 9\overline{d}\_0 r\_2\right)\sin(\Omega t) \end{aligned} \tag{34}$$

#### **4. Result and Discussion**

In this section, the numerical results for the system under consideration shown in Figures 1 and 2 are illustrated. Unless stated otherwise, the material and geometrical parameters of the AFM cantilever are listed in Table 1. Because the influence of the size effect has been well addressed in previous studies, the material length-scale parameter is taken as zero in the numerical simulation and only has been expressed in the mathematical formulation. The main aim of this section is to investigate the principal parametric resonance of the system [47].

**Table 1.** The material and geometrical parameters of the AFM cantilever.


As seen in Equations (33) and (34), the combination of van der Waals force and base-excitation results in time-dependent linear inertia, stiffness, damping terms, and harmonically time-varying nonlinear stiffness terms, creating parametric parameters excitations. Moreover, the AFM is affected by two external excitations having the same frequencies, which are equal to the frequency of the base excitation. It is observed that the frequency of all parametric excitations appearing in the linear, quadratic, and cubic nonlinearities is equal to that of the external excitation. In this work, the quality factor of the system is assumed to be 200, defining a low-damped vibrating system. Under this

condition, if the frequency of the parametric excitation varies twice the fundamental natural frequency of the AFM, principal parametric resonance is activated for small values of the base-excitation amplitude, *d*0. In other words, for the quality factor of 200, the activation level of parametric resonance is reached at small values of *d*<sup>0</sup> that can let us neglect the contribution of the time-varying term existing in the inertia term. Another point that is worth mentioning is that, while the frequency of the base excitation is close to two times the first natural frequency, the sub-harmonic resonance of the first mode can be activated by the external excitation; however, its activation level is much greater than that of the principal parametric resonance. Therefore, it can be concluded that this type of nonlinear resonance does not contribute to the system's response within the variation range of *d*0, which results in parametric resonance. With this in mind, the frequency-response behavior of the AFM is captured for different values of the base excitation amplitude and illustrated in Figure 3.

**Figure 3.** Frequency-response behavior of the AFM in the neighbor of its principal parametric resonance for different values of the base-excitation amplitude.

The figure shows that the frequency–response curves are composed of trivial stable, trivial unstable, nontrivial stable, and nontrivial unstable branches. The term *trivial* refers to the periodic response with zero amplitude; however, *nontrivial* phrase returns to the periodic orbit with non-zero amplitude. Moreover, the term *stable* clarifies that the system's states are absorbed by the periodic orbit after the system is disturbed; however, this is not the case for the *unstable* term, meaning that any small disturbance applied to the system's dynamics causes the system's states to get off the periodic orbit without the ability of returning back to it. The system's dynamics begin with stable zero-amplitude solutions and continue until primary Hopf bifurcation (sub-critical) occurs, and the stable branch loses its stability. Further increasing the parametric frequency, the unstable trivial solutions retrieve their stability at super-critical Hopf bifurcation points. In addition, for the smallest value of *d*<sup>0</sup> parameter close to the activation level, the stable and unstable nontrivial branches meet at the cyclic-fold bifurcation point at Ω ≈ 29.07; however, this is not the case for larger values of the parameter. As the parametric pump enhances, the resonance bandwidth becomes broader. The non-zero stable and unstable solutions meet at a displacement value

beyond the gap distance between the microbeam and the substrate. Furthermore, it can be inferred that the quadratic and cubic nonlinearities arising from the intermolecular force induce a softening effect on the steady-state dynamics of the AFM, making the nontrivial solution branches bend to the left. It should be mentioned that Wmax is the displacement of the cantilever tip. For better describing the frequency–response curves, it should be noted that there has been considered a set of two curves corresponding to each color. For instance, the blue plot drawn in Figure 3 is composed of two curves; the one which is marked by a star showing the stable solution, and the one that is mark-free shows the unstable solution branch corresponding to the same value of *d*0.

Figure 4 demonstrates the loci of the primary Hopf bifurcation points for different values of the excitation amplitude. The vertical axis shows the ratio between the excitation amplitude to its threshold, leading to parametric resonance activation. As seen in the figure, while this ratio is close to one from below, there is still no bifurcation in the system's dynamics. However, while this ratio is close to one from above, two bifurcation points are occurring at two frequencies so close to each other that the left ones return to the sub-critical, and the right ones introduce the position of the super-critical Hopf bifurcation points. As the amplitude of the parametric excitation grows, the distance between these two points increases. The grey area bounded by the fitted blue curve displays the parametric resonance region, known as the instability tongue. It is worth noting that if the parametric excitations are not accompanied by external excitation, the parametric noise squeezing effect is potent to happen while the excitation amplitude is below the activation level. However, this phenomenon cannot occur in the present dynamics, because the parametric excitation is collaborated by external excitation with the same excitation frequency. Inspecting the numerical values on the vertical axis, one can find that the parametric resonance bandwidth is extremely sensitive to the excitation amplitude.

**Figure 4.** Instability tongue. The loci of both sub- and super-critical primary Hopf bifurcation points define the boundary of the parametric resonance region, where the threshold value for the base excitation amplitude is *dth* = 0.222 nm.

Figure 5 illustrates the influence of the incompressibility (*ν*) on the frequency-response behavior of the AFM system by choosing two different values for Poisson's ratio, for *d*<sup>0</sup> = 0.23 nm. Increasing the value of Poisson's ratio causes the structure to become stiffer and more rigid, and therefore, its natural frequencies grow remarkably. With this in mind, the frequency of the parametric excitation is swept in the neighbour of twice the first natural frequency, which is evaluated for that specific Poisson's ratio. In order to have a precise evaluation of how the impact of Poisson's ratio influences the slope of the frequency–response curve, the displacement amplitude is plotted versus the difference between the excitation frequency and two times the natural frequency. It is found from the figure that, not only the resonance region becomes more expansive, but also the softening trend is intensified. At the same time, the resonator is made from a more incompressible material.

**Figure 5.** The impact of Poisson's ratio on the frequency-displacement behavior of the AFM, for *d*<sup>0</sup> = 0.23 nm. The dashed and solid lines are unstable and stable branches, respectively.

#### *Profile Height Detection Mechanism*

In this section, a sensing mechanism to detect the height of a sample's profile that is implemented based on the principal parametric resonance characteristics is proposed. First, it is shown that while the AFM is operating near its parametric resonance, the bifurcation points' position is highly affected by the gap distance between the microcantilever and its sample underneath. Figure 6 demonstrates the frequency-response behavior of the AFM for the excitation amplitude of *d*<sup>0</sup> = 0.23 nm. Here, the initial gap distance is set to be 60 nm. As seen in the figure, the sub-critical and super-critical Hopf bifurcation points shift to the left and to the right, respectively, while there is a 1 nm rise in the height of the sample profile (orange curve). An increment in the height of the sample profile causes the gap distance between the cantilever tip and the sample top surface to decrease. Therefore, the impact of van der Waals's force enhances, and consequently, the activation level decreases so that the resonance region becomes broader while the excitation amplitude is kept constant. This 1 nm rise in the height of the sample causes a 0.01 nondimensional frequency shift at primary Hopf bifurcation points. On the other hand, if the sample profile undergoes a 1 nm surface depression, extending the gap between the cantilever tip and the top surface of the sample, weakening the van der Waals force, the activation level of the parametric

resonance increases. Hence, the resonance region shrinks while the parametric pump is kept constant, meaning the sub-critical and super-critical Hopf bifurcation points shift to the right and the left, respectively (black curve). As observed in the figure, the frequency shift caused by the rise in the height of the sample profile is quite a bit smaller than that caused by surface depression. Here, the frequency shift at bifurcation points caused by a 1 nm profile depression is obtained at about 0.0125. In Figure 6, PH stands for primary Hopf bifurcation.

**Figure 6.** The influence of the surface profile variations on the primary Hopf bifurcation points near parametric resonance for the base excitation amplitude of 0.23 nm. The dashed and solid lines are unstable and stable branches, respectively.

The frequency-response behavior of the AFM device near the primary resonance of its first mode is illustrated in Figure 7. To monitor the steady-state dynamic response of the system near its primary resonance, the external excitation frequency is set to be changing in the vicinity of the first natural frequency, and the displacement amplitude of the periodic orbit is recorded for each specific value of the excitation frequency. As seen in the figure, the frequency-displacement behavior does not experience bifurcation for the excitation amplitude of *d*<sup>0</sup> = 0.15 nm; however, cyclic-fold bifurcation occurs beyond the gap distance for larger values of the base excitation amplitude, which is not meaningful, and we prevented presenting that result here. It is worth mentioning that for this case, the secondary parametric resonance of the first mode (Ω ≈ *ωn*) of oscillation is probable to happen because, as stated before, the frequency of the parametric terms is identical to that of the external stimulus. While the excitation frequency varies near the natural frequency of the first mode itself, the system has the potential to experience the combination of both primary resonance (due to external/direct excitation) and secondary parametric resonance (due to parametric excitation) of the first mode. However, the secondary parametric resonance (Ω ≈ *ωn*) requires a higher level of activation compared to the primary parametric resonance (Ω ≈ 2*ωn*). Therefore, it cannot be motivated for the proposed AFM, because the cantilever experiences tapping mode, which is caused by largeamplitude oscillation due to direct excitation. Inspecting the figure, it is seen that while there is a 1 nm increase in the height of the sample profile, the resonance amplitude grows

slightly (orange curve) compared to the case in which there is no variation in the surface height (blue curve). On the other hand, the displacement amplitude at resonance drops quite a bit while the surface of the sample undergoes a 1 nm depression. Conversely, these variations in the profile configuration do not result in a change in the primary resonance frequency. In this work, the authors propose an effective detection mechanism based on amplitude shift at bifurcation points that are highly sensitive to detecting the high variations of samples' surface profiles. The AFM needs to be run near its parametric resonance zone instead of the primary resonance region.

**Figure 7.** The influence of the surface profile variations on the resonance amplitude of the microcantilever near primary resonance, for *d*<sup>0</sup> = 0.15 nm.

Figure 8 depicts the frequency-displacement amplitude of the AFM near its parametric resonance for a parametric pump which is slightly above the activation level, *<sup>d</sup>*<sup>0</sup> *dth* <sup>=</sup> 1.001. For the case in which there is neither bulge nor surface depression on the surface profile (violet curve), a cyclic-fold bifurcation exists at Ω = 29.068. This dynamic behavior corresponds to the case where the cantilever tip is precisely on top of the sample area from which the initial gap distance, *Z*ˆ, is measured. This profile status is considered the surface baseline, provided that any variation in the height of the sample profile concerning this baseline is denoted by *δ*. As observed in the figure, the frequency–response curve shifts to the right once the distance between the microcantilever tip and the surface profile increases by 15 pm; namely, the AFM faces a surface depression of 15 pm on the sample profile. Because the AFM is already set to be operating at point A, this right shift in the steady-state behavior of the device causes a significant drop in the displacement amplitude, jumping from point A down to point B with zero amplitude. The AFM's sensitivity can be obtained by evaluating the dimensional displacement drop to the change in the profile height. Accordingly, the sensitivity of the proposed method is as follows:

$$S\_{D\delta}^{\text{w}} = \frac{\Delta \text{w}}{\Delta \delta} = \frac{0.71 \times 60 \text{ nm}}{15 \text{ pm}} = 2840 \frac{\text{amplitude} \,\text{(nm)}}{\delta \,\text{(nm)}} \tag{35}$$

**Figure 8.** The influence of the surface profile variations on the amplitude of the resonator's displacement at cyclic-fold bifurcation near parametric resonance, for *<sup>d</sup>*<sup>0</sup> *dth* = 1.001. The dashed and solid lines are unstable and stable branches, respectively.

On the other hand, the frequency–response curve shifts to the left while there is a 15 pm rise in the sample profile, so the resonator's response experiences a jump up to point C. Although this enhancement in the displacement amplitude is smaller than what is achieved for surface depression, the detection sensitivity is still considerable compared to the case where the cantilever is actuated near its primary resonance (comparing Figures 7 and 8). Similarly, for the case in which the distance between the AFM tip and the sample profile diminishes, the sensitivity of the device is obtained as follows,

$$S\_{R\delta}^{\text{uv}} = \frac{\Delta \text{w}}{\Delta \delta} = \frac{(0.8 - 0.71) \times 60 \text{ nm}}{15 \text{ pm}} = 360 \frac{\text{amplitude (nm)}}{\delta \text{ (nm)}} \tag{36}$$

Scientifically speaking, because of the softening behavior of the AFM, the response amplitude drops from the stable nontrivial branch down to the stable trivial branch. At the same time, there is a positive variation in the height of the sample (surface depression), leading to ultrahigh sensitivity. However, this is not the case while there are negative variations in the height of the sample profile, meaning that the system's response jumps from a stable nontrivial branch up to a new generated stable nontrivial solution branch, degrading the AFM sensitivity. This can be improved by designing a tuneable AFM that can show both softening and hardening behavior near its parametric resonance regime. It is worth noting that the most challenging part in achieving ultrahigh sensitivity in AFMs returns to the employment of the base excitation amplitude. The finer the resolution of the parametric pump (parametric excitation amplitude), the higher the possibility for the implementation of the proposed detection mechanism. This means that the suggested method requires the capability of changing the base excitation amplitude with fine resolution. The proposed method can be feasibly employed once this concern is responded to.

In this paper, a micropolar model Equation (4) suitable for hyperelastic materials was utilized. However, other hyperelastic models can be utilized for hyperelastic microcantilevers; for example, the Gent model [48], generalized neo-Hookean models [49], Gent-Gent model [50], etc.

The foremost aspect of the analysis in the present work is that the system operates in parametric resonance regime. We have identified such a regime for non-contact AFM that is common in many real setups. For other kinds of AFM, it is essential to seek the presence of parametric resonance in the system first. Then, the analysis of the current study can be developed for the system.

#### **5. Conclusions**

This work examined the nonlinear resonance of a hyperelastic-based cantilever AFM in the non-contact region so that the non-contact force is modeled using the van der Waals force. A vertical base excitation is used to vibrate the AFM. A hyperelastic model combined with the modified couple stress theory is proposed to incorporate the elastic energy stored in the microcantilever caused by deformation. The obtained equation of motion is discretized via a developed Galerkin method. Furthermore, an efficient detection mechanism based on principal parametric resonance, assessing surface depression and height increase of a sample profile, is proposed. The following concluding remarks can be deduced from the presented numerical results.


**Author Contributions:** Conceptualization, A.A. and S.D.; methodology, A.A., S.R. and S.D.; software, A.A. and S.R.; validation, S.D., M.M., B.K., B.A. and Ö.C.; formal analysis, A.A. and S.R.; investigation, M.M., B.K., B.A. and Ö.C.; resources, A.A. and S.D.; data curation, A.A. and S.R.; writing—original draft preparation, A.A. and S.D.; writing—review and editing, A.A., S.D., M.M., S.R., B.K., B.A. and Ö.C.; visualization, A.A., S.D. and S.R.; supervision, M.M., B.K., B.A. and Ö.C.; funding acquisition, S.D. and M.M.; project administration, Ö.C. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding. The APC was funded by S.D. and M.M.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** There is no data for sharing and all data are available within the paper.

**Conflicts of Interest:** There are no conflicts of interest.

#### **References**

