Nonlinear Dynamic Response of Nanocomposite Microbeams Array for Multiple Mass Sensing

Abstract

A nonlinear MEMS multimass sensor is numerically investigated, designed as a single input-single output (SISO) system consisting of an array of nonlinear microcantilevers clamped to a shuttle mass which, in turn, is constrained by a linear spring and a dashpot. The microcantilevers are made of a nanostructured material, a polymeric hosting matrix reinforced by aligned carbon nanotubes (CNT). The linear as well as the nonlinear detection capabilities of the device are explored by computing the shifts of the frequency response peaks caused by the mass deposition onto one or more microcantilever tips. The frequency response curves of the device are obtained by a pathfollowing algorithm applied to the reduced-order model of the system. The microcantilevers are described by a nonlinear Euler-Bernoulli inextensible beam theory, which is enriched by a meso-scale constitutive law of the nanocomposite. In particular, the microcantilever constitutive law depends on the CNT volume fraction suitably used for each cantilever to tune the frequency bandwidth of the whole device. Through an extensive numerical campaign, the mass sensor sensitivity estimated in the linear and nonlinear dynamic range shows that, for relatively large displacements, the accuracy of the added mass detectability can be improved due to the larger nonlinear frequency shifts at resonance (up to 12%).

1. Introduction

In the past decades MEMS have been subject to disruptive developments in several fields of engineering both in research and applications (mainly, signal processing, actuation, and sensing). The lower power consumption and higher sensitivity are distinct advantages with respect to comparable macro devices. A general introduction about MEMS devices can be found in [1] with an overview of the fundamental concepts pertaining to linear and nonlinear dynamics of such systems. In recent years, a great deal of studies have sought to advance our understanding of the nonlinear dynamic behavior of nano/microresonators, which represent the main components of actual NEMS/MEMS devices [2,3,4,5]. Nonlinear dynamic aspects of electrostatically actuated microcantilevers are discussed in [6] while dynamics of the same systems subject to piezoelectric actuation are studied in [7].

Accurate models still need to achieve better agreement between analytical/numerical predictions and experiment results, and explore the use of new materials and technical capabilities. Initial imperfections of the structure, together with electrostatic actuation and nonlinear damping effects, are discussed in [8,9]. Moreover, systems with feedback loop control are investigated to enhance the sensitivity of the microdevices [10,11]. By restricting the attention to micromass sensing, the literature is becoming increasingly richer due to the present interest in gas detection, chemical agent detection, biological molecules and pollutant particles identification. In mass/gas detection, different physical principles can be used to measure the different target substances. In [12,13,14] some electrical resistive detectors are investigated. In these studies the functioning of a semiconductor placed in an integrated circuit aimed to measure target molecules such as SO${}_2$ or NO${}_2$ in order to perform environmental control or to measure the pollution associated with combustion chamber exhaust gases. Moreover, an example of a MEMS gas sensor using Knudsen force for hydrogen detection can be found in [15,16]. Therefore, MEMS gas sensors using the mechanical features of micro oscillators are of current interest. An overview of the working principles of these devices can be found in [17]. A micro mass sensor based on a microcantilever is discussed in [18]. This micromechanical resonator is used as gas sensor exploiting a working principle based on the detection of the beam natural frequency shift due to mass deposition onto the tip of the microbeam. In addition, similar works addressing microcantilevers used as sensor resonators to detect the chemical concentration of matter are found in [19,20]. On the other hand, a detection procedure is discussed in [21,22] where monitoring of the amplitudes of the localized vibration modes that vary when the system mass is perturbed was performed. A device made of two microcantilevers weakly coupled through an overhang working with localized modes is discussed in [23] while the use of these tools in highly viscous fluids is explored in [24,25]. Different kinds of electrostatic actuation are addressed in [26,27] to achieve large motions as well as higher sensitivity. In the field of NEMS, the performance of standalone carbon nanotube oscillators for mass detection is investigated in [28,29].

Numerous studies have been carried out to explore new mass detection strategies exploiting nonlinear effects. A parametrically excited clamped-clamped beam with cubic stiffness nonlinearity is explored as micromass sensor in [30]. The principle is based on the passage from the stable region to the unstable region across the transition curves of the Strutt diagrams. Also the jumps between the nonresonant and resonant branches occurring at the fold bifurcations along the frequency response curves are used as detection principle in [31,32]. Moreover, an electrostatically actuated microcantilever dynamic mass sensor was analytically investigated and experimentally tested in [33] to indirectly measure the added mass and its position along the span exploiting the frequency shift. The nonlinear dynamic behavior of an electrostatically actuated microbeam was studied in [34] to compare the different kinds of micromass detections according to the jumps at the bifurcation points, the frequency sweeps, and the frequency shifts of the resonant peaks. A specific study on the static and dynamic effects of the pull-in phenomena for a MEMS gas sensor can be found in [35]. For an extensive and up-to-date review, the reader can refer to [36], where a comprehensive categorization of micro/nanomechanical resonators and their engineering applications is provided, including a classification of NEMS/MEMS devices based on linear/nonlinear, single/array, symmetric/asymmetric, together with frequency shift-based and amplitude shift-based resonators.

In the present work, a single-input single-output (SISO) system is investigated making use of a path following algorithm optimized to tackle several degrees of freedom. The baseline idea of the SISO system can be found in [37] where a dynamic lumped system is fabricated to work in the linear range, and then successfully experimented for detection purposes, with experimental results in good agreement with a previous theoretical contribution of the same author [38]. The working principle of the device consists in the use of the shuttle mass motion both to excite the sensors and to read off the overall dynamic response of the system which includes the vibrations of all microcantilevers.

One of the key aspects of the proposed design is the use of CNT nanocomposite material (i.e., polymeric matrix hosting aligned CNTs forests) for the microcantilever detectors. Due to their excellent mechanical, thermal, electrical, and tribological properties, compared with their lightweight nature, CNT-based nanocomposites reached in the last decades a widespread use in various engineering applications (see, e.g., [39,40,41,42,43] specifically for reinforced polymers). For sensor applications, their ability of specific and stable detections of mechanical and electrochemical properties was extensively explored in the literature: for an up-to-date survey on this, see for instance [44,45]. An experimental and numerical study about the response of such nanostructured materials is given in [46,47,48], where the accuracy of the constitutive modeling here employed was also successfully validated.

On the other hand, the analytical investigation via a perturbation technique of the design guidelines for SISO systems can be found in [49]. In the present work a different approach based on the eigenvalue analysis for a five-dof system is presented. Moreover, contrary to the conventional way of tuning the detector frequencies by changing their length, the proposed design consists of an array of equal-length beams having different CNT volume fractions which leads to easily tunable frequencies. The shuttle mass is a mass-spring-damper system while the microcantilevers clamped to the shuttle mass are modeled as nonlinear Euler-Bernoulli beams, enhanced by the constitutive relationship derived at meso-scale level for CNT-nanocomposite beams as proposed in [50]. The equations of motion are discretized via the Galerkin method and the system sensitivity is investigated for different cases such as the case of the system driven by a harmonic excitation or by an electrostatic actuation as well as the case of the system with one or more working sensor beams. Finally, the sensitivity of the device is computed in the linear and nonlinear dynamic ranges proving that the nonlinearities can enhance the detection capabilities.

2. Equations of Motion for the Nanocomposite Microcantilevers

The SISO system investigated in the present study is depicted in  Figure 1.

The flexural motion of the jth microcantilever is described by the nonlinear Euler-Bernoulli beam theory: for full details, we refer the interested reader to [51]. Let s and t be the arclength along the beam baseline and time, respectively, and let ${\mathit{\partial }}_s$ and ${\mathit{\partial }}_t$ denote differentiation with respect to s and t, respectively. The nonlinear equation of motion for the jth undamped microcantilever is given by the following integral-partial differential equation:

$$\begin{array}{cc}\hfill & \rho A({\partial }_{tt}{v}_{j}cos{\theta }_j-{\partial }_{tt}{u}_{j}sin{\theta }_j)-{\kappa }_j{N}_j^{\ast }+{\partial }_{ss}{M}_j\hfill \\ \hfill & +{\kappa }_j{\int }_s^l\left[\rho A({\partial }_{tt}{u}_{j}cos{\theta }_j+{\partial }_{tt}{v}_{j}sin{\theta }_j)-{\kappa }_j{\partial }_s{M}_j\right]d\xi =0,\hfill \end{array}$$

(1)

where $N_j^{*}$ is the normal force at the tip of the cantilever, which is due to the presence of the added mass at the tip denoted by $m_j$; $\rho A$ is the beam mass per unit length, with l being its span, while $u_j$ and $v_j$ are the axial and transverse (i.e., with respect to the fixed directions ${\mathit{e}}_3$ and ${\mathit{e}}_1$) displacements of the cantilevers, respectively. Furthermore, ${\theta }_j$ denotes the counterclockwise rotation of the beam cross section, ${\kappa }_j$ is the flexural curvature, and $M_j=E{J}_j^{\mathrm{eq}}{\kappa }_j$ is the elastic bending moment, proportional to the equivalent flexural stiffness $E{J}_j^{\mathrm{eq}}$ of the jth nanocomposite cantilever cross section. In particular, according to [47,50,52,53], such equivalent flexural stiffness is derived by a 3D continuum homogeneization model which provides the following expression:

$$E{J}_j^{\mathrm{eq}}={\widehat{EJ}}^{\mathrm{eq}}({\varphi }_{\mathrm{C}},E_{\mathrm{M}},E_{\mathrm{C}},{\nu }_{\mathrm{M}},{\nu }_{\mathrm{C}}),$$

(2)

where the equivalent stiffness depends on the volume fraction ${\varphi }_{\mathrm{C}}$ of the CNTs embedded in the polymeric hosting matrix, and of the isotropic elastic moduli of the two materials, Young’s and Poisson’s moduli $E_{\mathrm{C}}$ and ${\nu }_{\mathrm{C}}$, respectively for the CNTS, Young’s and Poisson’s moduli $E_{\mathrm{M}}$ and ${\nu }_{\mathrm{M}}$, respectively for the hosting polymeric matrix.

The boundary conditions for the jth microcantiliver considering the presence of the added mass $m_j$ at its tip are:

$$\begin{array}{cc}\hfill & v_j(0,t)=x\left(t\right),{\partial }_s{v}_j(0,t)=0,{\partial }_{ss}{v}_j(l,t)=0,\hfill \\ & E{J}_j^{\mathrm{eq}}{\partial }_{sss}{v}_j(l,t)=m_j{\partial }_{tt}{v}_j(l,t).\hfill \end{array}$$

(3)

Note that the dynamic boundary condition involving the tip mass $m_j$ is obtained via linearization of the fully nonlinear equation of motion for the lumped tip mass. Such equation enforces the balance of linear momentum at the cantilever tip expressed as

$$(N_j^{*}{\mathbf{b}}_3+Q_j^{*}{\mathbf{b}}_1)=-m_j\left[{\partial }_{tt}{u}_j(l,t){\mathbf{e}}_3+{\partial }_{tt}{v}_j(l,t){\mathbf{e}}_1\right].$$

(4)

where $({\mathbf{b}}_1,{\mathbf{b}}_3)$ are the beam cross section-fixed unit vectors providing the current cross section orientation and the direction normal to it while $({\mathbf{e}}_1,{\mathbf{e}}_3)$ indicate the fixed horizontal and vertical directions, respectively. Therefore, the jth tension $N_j^{*}$ and shear force $Q_j^{*}$ can be obtained as

$$\begin{array}{c}{N}_j^{*}=-m_j\left[cos{\theta }_j(l,t){\partial }_{tt}{u}_j(l,t)+sin{\theta }_j(l,t){\partial }_{tt}{v}_j(l,t)\right],\end{array}$$

(5)

$$\begin{array}{c}{Q}_j^{*}=-m_j\left[-sin{\theta }_j(l,t){\partial }_{tt}{u}_j(l,t)+cos{\theta }_j(l,t){\partial }_{tt}{v}_j(l,t)\right].\end{array}$$

(6)

The subsequent linearization of Equation (6) yields

$$Q_j^{*}=-m_j{\partial }_{tt}{v}_j(l,t).$$

(7)

By further expanding Equation (1) in Taylor series up to the third order adopting the moment-curvature law proposed in [52,53], and further developed in [50], expressing the longitudinal motion $u_j$ as a function of the transverse motion $v_j$, the Euler-Bernoulli equations of motion governing the transverse dynamics of the jth microbeam with the tip mass $m_j$ are cast in the form

$$\begin{array}{cc}\hfill & \rho A{\mathit{\partial }}_{tt}{v}_j\left(1-\frac{1}{2}{\left({\mathit{\partial }}_s{v}_j\right)}^2\right)+\rho A{\mathit{\partial }}_s{v}_j{\int }_0^s\left[{\mathit{\partial }}_s{v}_j{\mathit{\partial }}_{stt}{v}_j+{\left({\mathit{\partial }}_{st}{v}_j\right)}^2\right]d\xi \hfill \\ \hfill & +{\mathit{\partial }}_{ss}{v}_j{\int }_s^l\rho A\left[{\mathit{\partial }}_s{v}_j{\mathit{\partial }}_{tt}{v}_j-{\int }_0^{\xi }\left({\mathit{\partial }}_s{v}_j{\mathit{\partial }}_{stt}{v}_j+{\left({\mathit{\partial }}_{st}{v}_j\right)}^2\right)dz\right]d\xi \hfill \\ \hfill & +{\mathit{\partial }}_{ss}{M}_j-{\mathit{\partial }}_{ss}{v}_j{\int }_s^l{\mathit{\partial }}_{ss}{v}_j{\mathit{\partial }}_s{M}_{j}d\xi \hfill \\ \hfill & +m_j{\mathit{\partial }}_{ss}{v}_j\left[\left({\mathit{\partial }}_s{v}_j\right){|}_l\left({\mathit{\partial }}_{tt}{v}_j\right){|}_l+{\int }_0^l\left({\mathit{\partial }}_s{v}_j{\mathit{\partial }}_{stt}{v}_j+{\left({\mathit{\partial }}_{st}{v}_j\right)}^2\right)d\xi \right]=0,\hfill \end{array}$$

(8)

where $v_j$ is the transverse deflection of the beam, the constitutive law for the bending moment is

$$M_j=E{J}_j^{\mathrm{eq}}{\kappa }_j=E{J}_j^{\mathrm{eq}}\left[{\mathit{\partial }}_{ss}{v}_j\left(1+\frac{1}{2}{\left({\mathit{\partial }}_s{v}_j\right)}^2\right)\right]\mathrm{with}{\kappa }_j={\mathit{\partial }}_{ss}{v}_j\left(1+\frac{1}{2}{\left({\mathit{\partial }}_s{v}_j\right)}^2\right).$$

Finally, the equation of motion for the shuttle mass is given by

$$M{\mathit{\partial }}_{tt}x+c{\mathit{\partial }}_{t}x+Kx+\sum _{j=1}^{N_b}{\mathit{\partial }}_s{M}_j{|}_{s=0}=fcos\mathsf{\Omega }t,$$

(9)

where x is the displacement of the shuttle mass M, c is the damping coefficient of the dashpot, k the spring constant, $f\left(t\right)=Fcos\mathsf{\Omega }t$ with F and $\mathsf{\Omega }$ being the excitation amplitude and frequency, respectively. The last term is the sum the shear forces of the $N_b$ beams at the roots attached to the shuttle mass, on accounting for the relationship between shear force and moment gradient $Q=-{\mathit{\partial }}_{s}M$. In fact, the associated boundary conditions for the beams are the kinematic relations for the clamp and the mechanical boundary conditions at the tip

$$v_j(0,t)={\mathit{\partial }}_s{v}_j(0,t)=x\left(t\right)=0$$

and

$$M_j(l,t)=0,Q_j(l,t)=-m_j({\mathit{\partial }}_{tt}x+{\mathit{\partial }}_{tt}{v}_j){|}_{s=l}.$$

Expressed in terms of displacements, the mechanical the boundary conditions at the tip become

$$E{J}_j^{\mathrm{eq}}{\mathit{\partial }}_{ss}{v}_j(l,t)=0,E{J}_j^{\mathrm{eq}}{\mathit{\partial }}_{sss}{v}_j(l,t)=m_j({\mathit{\partial }}_{tt}x+{\mathit{\partial }}_{tt}{v}_j){|}_{s=l}.$$

We rescale the absolute displacements $v_j$ of the jth beam and the shuttle mass displacement x by the length of the beams l, time by $1/{\omega }_c$ with ${\omega }_c^2=E{J}^{\mathrm{eq}}/\left(\rho A{l}^4\right)$, where $\rho A$ is the mass per unit length of the reference beam whose bending stiffness is denoted by $E{J}^{\mathrm{eq}}$. By retaining the symbols s and t as the beam nondimensional arclength and time, respectively, and using the prime for space derivatives and the over dot for time derivatives, the ensuing nondimensional Equation (8) become:

$$\begin{array}{cc}\hfill & \mu \ddot{x}\left(t\right)+\xi \dot{x}\left(t\right)+K^{*}x\left(t\right)+\sum _{j=1}^{N_b}{\alpha }_j{v}_j^{\prime \prime \prime }(0,t)=f^{*}cos{\mathsf{\Omega }}^{*}t,\hfill \\ \hfill & {\mu }_j{\ddot{v}}_j(s,t)+m_j^{*}\delta (s-1){\ddot{v}}_j(s,t)+{\xi }_j{\dot{v}}_j(s,t)+{\alpha }_j{v}_j^{\prime \prime \prime \prime }(s,t)\hfill \\ \hfill & +{\mu }_j{f}_j^I[v(s,t),\dot{v}(s,t),\ddot{v}(s,t)]+{\alpha }_j{f}_j^R\left[v(s,t)\right]\hfill \\ \hfill & +m_j^{*}{f}_j^{I*}[v(s,t),\dot{v}(s,t),\ddot{v}(s,t)]=0,\mathrm{for}j=1,\dots ,N_b,\hfill \end{array}$$

(10)

where

$$\begin{array}{cc}\hfill f_j^I[v,\dot{v},\ddot{v}]& =-\frac{1}{2}{\ddot{v}}_j{\left(v_j^{\prime }\right)}^2+v_j^{\prime \prime }{\int }_{s^{★}}^1\left[{\ddot{v}}_j{v}_j^{\prime }\right]d{\xi }^{★}+v_j^{\prime }{\int }_0^{s^{★}}\left[v_j^{\prime }{\ddot{v}}_j^{\prime }\right]d{\xi }^{★}\hfill \\ & -v_j^{\prime \prime }{\int }_{s^{★}}^1\left[{\int }_0^{s^{★}}\left[v_j^{\prime }{\ddot{v}}_j^{\prime }\right]d{\xi }^{★}\right]d{s}^{★}+v_j^{\prime }{\int }_0^{s^{★}}\left[{\left({\dot{v}}_j^{\prime }\right)}^2\right]d{\xi }^{★}\hfill \\ & -v_j^{\prime \prime }{\int }_{s^{★}}^1\left[{\int }_0^{s^{★}}\left[{\left({\dot{v}}_j^{\prime }\right)}^2\right]d{\xi }^{★}\right]d{s}^{★},\hfill \\ \hfill f_j^R\left[v\right]& ={\left(v_j^{\prime \prime }\right)}^3+\frac{1}{2}{v}_j^{\prime \prime \prime \prime }{\left(v_j^{\prime }\right)}^2+3v_j^{\prime \prime \prime }{v}_j^{\prime \prime }{v}_j^{\prime }-v_j^{\prime \prime }{\int }_{s^{★}}^1(v_j^{\prime \prime }{v}_j^{\prime \prime \prime }+\dots )d{\xi }^{★},\hfill \\ \hfill f_j^{I*}[v,\dot{v},\ddot{v}]& =v_j^{\prime \prime }\left[\left(v_j^{\prime }\right){|}_1\left({\ddot{v}}_j\right){|}_1-{\int }_0^1\left(v_j^{\prime }{\ddot{v}}_j^{\prime }+{\left({\dot{v}}_j^{\prime }\right)}^2\right)d\xi \right].\hfill \end{array}$$

In the governing equations, we introduced the following nondimensional parameters:

$$\begin{array}{c}\mu :=M/\left(\rho Al\right),{\mu }_j:=\rho A_j/\rho A,m_j^{*}:=m_j/\left(\rho Al\right),\\ \xi :=c/\left(\rho Al{\omega }_c\right),{\xi }_j:=c_j/\left(\rho A{\omega }_c\right),\\ K^{*}:=K/\left(\rho Al{\omega }_c^2\right),{\alpha }_j:=E{J}_j^{\mathrm{eq}}/\left(\rho A{\omega }_c^2{l}^4\right),\\ f^{*}:=f/\left(\rho A{l}^2{\omega }_c^2\right),{\mathsf{\Omega }}^{*}:=\mathsf{\Omega }/{\omega }_c,\end{array}$$

where $m_j$ is the tip mass of the jth beam, ${\xi }_j$ its damping coefficient, ($\mu$, $K^{*}$, $\xi$) are the nondimensional mass, spring constant and damping coefficient of the shuttle mass, $f^{*}$ and ${\mathsf{\Omega }}^{*}$ are the nondimensional excitation amplitude and frequency, respectively.

We further introduce the beam deflections $w_j$ relative to the shuttle mass according to $v_j(s,t)=w_j(s,t)+x\left(t\right)$. The equations of motion in terms of the relative beam deflections $w_j$ take the following form:

$$\begin{array}{cc}\hfill & \mu \ddot{x}\left(t\right)+\xi \dot{x}\left(t\right)+K^{*}x\left(t\right)+\sum _{j=1}^{N_b}{\alpha }_j{w}_j^{\prime \prime \prime }(0,t)=f^{*}cos{\mathsf{\Omega }}^{*}t,\hfill \\ \hfill & {\mu }_j({\ddot{w}}_j(s,t)+\ddot{x}\left(t\right))+m_j^{*}\delta (s-1)({\ddot{w}}_j(s,t)+\ddot{x}\left(t\right))+{\xi }_j{\dot{w}}_j(s,t)+{\alpha }_j{w}_j^{\prime \prime \prime \prime }(s,t)\hfill \\ & +{\alpha }_j{f}_j^R\left[w(s,t)\right]+{\mu }_j{f}_j^I[w(s,t),\dot{w}(s,t),\ddot{w}(s,t)]+m_j^{*}{f}_j^{I*}[w(s,t),\dot{w}(s,t),\ddot{w}(s,t)]\hfill \\ & +{\mu }_j{g}_j^I\left[w(s,t)\right]\ddot{x}\left(t\right)+m_j^{*}{g}_j^{I*}\left[w(s,t)\right]\ddot{x}\left(t\right)=0,\mathrm{for}j=1,\dots ,N_b,\hfill \end{array}$$

(11)

where

$$\begin{array}{c}\hfill g_j^I\left[w\right]=-\frac{1}{2}{\left(w_j^{\prime }\right)}^2+w_j^{\prime \prime }{\int }_{s^{★}}^1\left[w_j^{\prime }\right]d{\xi }^{★},g_j^{I*}\left[w\right]=w_j^{\prime \prime }\left(w_j^{\prime }\right){|}_1.\end{array}$$

3. Modal Analysis

The modal analysis of the SISO system here discussed is the first step towards the nonlinear analyses carried out on the Galerkin obtained reduced order models. Moreover, the variability of the eigenfrequencies of the SISO system reveals interesting preliminary results, thanks to the nanostructured nature of the materials constituting the microbeams.

Let the jth overall mode shape of the system be cast in the vector-valued form:

$${\mathbf{u}}_j^{\top }=[X_j,{\varphi }_{1,j}\left(s\right),{\varphi }_{2,j}\left(s\right),\dots ,{\varphi }_{N_b,j}\left(s\right)]=[X_j,{\varphi }_{k,j}\left(s\right)],k=1,\dots N_b.$$

(12)

For further details see Appendix A. The substitution of the jth mode in the nondimensional kth beam equation and shuttle mass equation yields

$$\begin{array}{c}-{\mu }_k{\omega }_j^2{\varphi }_{k,j}+{\alpha }_k{\varphi }_{k,j}^{\prime \prime \prime \prime }=0,\end{array}$$

(13)

$$\begin{array}{c}-\mu {\omega }_j^2{X}_j+K^{*}{X}_j+\sum _k{\alpha }_k{\varphi }_{k,j}^{\prime \prime \prime }\left(0\right)=0;\end{array}$$

(14)

the associated boundary conditions for the beams are

$$\begin{array}{c}{\varphi }_{k,j}\left(0\right)-X_j=0,{\varphi }_{k,j}^{\prime }\left(0\right)=0,\end{array}$$

(15)

$$\begin{array}{c}{\alpha }_k{\varphi }_{k,j}^{\prime \prime }\left(1\right)=0,{\alpha }_k{\varphi }_{k,j}^{\prime \prime \prime }\left(1\right)+m_k^{*}{\omega }_j^2{\varphi }_{k,j}\left(1\right)=0.\end{array}$$

(16)

We multiply Equation (13) by ${\varphi }_{k,i}$, integrate over $[0,1]$ and sum the resulting projected equations for all beams; we multiply Equation (14) by $X_i$ and sum the previously projected beam equations and the shuttle mass equation to obtain:

$$\begin{array}{cc}\hfill -{\omega }_j^2[\sum _{k=1}^{N_b}& {\mu }_k{\int }_0^1{\varphi }_{k,i}{\varphi }_{k,j}ds+\mu X_i{X}_j]+\sum _{k=1}^{N_b}{\alpha }_k{\int }_0^1{\varphi }_{k,i}{\varphi }_{k,j}^{\prime \prime \prime \prime }ds\hfill \\ \hfill & +K^{*}{X}_i{X}_j+X_i\sum _{k=1}^{N_b}{\alpha }_k{\varphi }_{k,j}^{\prime \prime \prime }\left(0\right).=0\hfill \end{array}$$

(17)

We further integrate by parts the bending terms in the beam equations to obtain

$${\int }_0^1{\varphi }_{k,i}{\varphi }_{k,j}^{\prime \prime \prime \prime }ds={\varphi }_{k,i}{\varphi }_{k,j}^{\prime \prime \prime }{{|}_{s=0}^1-{\varphi }_{k,i}^{\prime }{\varphi }_{k,j}^{\prime \prime }|}_{s=0}^1+{\int }_0^1{\varphi }_{k,i}^{\prime \prime }{\varphi }_{k,j}^{\prime \prime }ds.$$

(18)

Employing the boundary conditions (15) yields

$${\int }_0^1{\varphi }_{k,i}{\varphi }_{k,j}^{\prime \prime \prime \prime }ds=-{\varphi }_{k,i}\left(1\right){\omega }_j^2\frac{m_k^{*}}{{\alpha }_k}{\varphi }_{k,j}\left(1\right)-X_i{\varphi }_{k,j}^{\prime \prime \prime }\left(0\right)+{\int }_0^1{\varphi }_{k,i}^{\prime \prime }{\varphi }_{k,j}^{\prime \prime }ds.$$

(19)

Substituting (19) into (20) yields

$$\begin{array}{cc}\hfill \sum _{k=1}^{N_b}[& -{\mu }_k{\omega }_j^2{\int }_0^1{\varphi }_{k,i}{\varphi }_{k,j}ds-{\omega }_j^2{m}_k^{*}{\varphi }_{k,i}\left(1\right){\varphi }_{k,j}\left(1\right)-{\alpha }_k{X}_i{\varphi }_{k,j}^{\prime \prime \prime }\left(0\right)+{\alpha }_k{\int }_0^1{\varphi }_{k,i}^{\prime \prime }{\varphi }_{k,j}^{\prime \prime }ds]\hfill \\ \hfill & -\mu {\omega }_j^2{X}_i{X}_j+K^{*}{X}_i{X}_j+X_i\sum _{k=1}^{N_b}{\alpha }_k{\varphi }_{k,j}^{\prime \prime \prime }\left(0\right)=0,\hfill \end{array}$$

(20)

which gets simplified into

$$\begin{array}{cc}\hfill -{\omega }_j^2[& \sum _{k=1}^{N_b}{\mu }_k{\int }_0^1{\varphi }_{k,i}{\varphi }_{k,j}ds+\mu X_i{X}_j+\sum _{k=1}^{N_b}{m}_k^{*}{\varphi }_{k,i}\left(1\right){\varphi }_{k,j}\left(1\right)]\hfill \end{array}$$

(21)

$$\begin{array}{cc}\hfill & +\left[\sum _{k=1}^{N_b}{\alpha }_k{\int }_0^1{\varphi }_{k,i}^{\prime \prime }{\varphi }_{k,j}^{\prime \prime }ds+K^{*}{X}_i{X}_j\right]=0.\hfill \end{array}$$

(22)

so that the generalized mass and stiffness matrices can be defined as follows:

$$\begin{array}{cc}\hfill M_{ij}=[& \sum _{k=1}^{N_b}{\mu }_k{\int }_0^1{\varphi }_{k,i}{\varphi }_{k,j}ds+\mu X_i{X}_j+\sum _{k=1}^{N_b}{m}_k^{*}{\varphi }_{k,i}\left(1\right){\varphi }_{k,j}\left(1\right)]\hfill \end{array}$$

(23)

$$\begin{array}{cc}\hfill & K_{ij}=\left[\sum _{k=1}^{N_b}{\alpha }_k{\int }_0^1{\varphi }_{k,i}^{\prime \prime }{\varphi }_{k,j}^{\prime \prime }ds+K^{*}{X}_i{X}_j\right].\hfill \end{array}$$

(24)

Note that the normalization condition is set to be $M_{jj}=1$, which yields

$$\sum _{k=1}^{N_b}{\mu }_k{\int }_0^1{\varphi }_{k,j}{\varphi }_{k,j}ds+\mu X_j{X}_j+\sum _{k=1}^{N_b}{m}_k^{*}{\varphi }_{k,j}\left(1\right){\varphi }_{k,j}\left(1\right)=1.$$

(25)

This allows to compute the coefficients $\{a_1,\dots ,a_{N_b}\}$ of the eigenvectors. In particular, the $(N_b-1)$ coefficients $\{a_2,\dots ,a_{N_b}\}$ are first obtained as function of $a_1$ through the system $A_{ij}{a}_j=0$ (for $i=2,\dots ,N_b$), and finally $a_1$ is computed by imposing the previous normalization condition. According to this mass normalization, the frequencies are obtained as

$$\begin{array}{c}\hfill {\omega }_j^2=K_{jj}=\sum _{k=1}^{N_b}{\alpha }_k{\int }_0^1{\varphi }_{k,j}^{\prime \prime }{\varphi }_{k,j}^{\prime \prime }ds+K^{*}{X}_j{X}_j.\end{array}$$

(26)

Numerical results. The mode shapes of the system, together with the corresponding natural frequencies, are computed using the linear equations of motion and solving the associated eigenvalue problem discussed in the previous section. Here and henceforth, we will consider $N_b=4$ microcantilevers attached to the shuttle mass.

The shuttle stiffness is set according to the formula $k={\omega }^{2}M$, $\omega$ being a percentage of the lowest natural frequency of the sensors. The mechanical parameters are set to the following parameters: $k={10}^5\mathrm{N}\mathrm{/}\mathrm{m}$, $M=40\rho \mathrm{Al}$, $\rho =1236\mathrm{kg}/{\mathrm{m}}^3$, $h=b=50\mathsf{\mu }\mathrm{m}$, $E_{\mathrm{M}}=2.8\mathrm{GPa}$, $E_{\mathrm{C}}=970\mathrm{GPa}$, ${\nu }_{\mathrm{M}}=0.3$, ${\nu }_{\mathrm{C}}=0.1$, $l=270\mathsf{\mu }\mathrm{m}$. Here, $\rho$ is the mass density of the beam hosting matrix (i.e., epoxy), and h and b are the thickness and width of the microcantilevers cross section, respectively, so that $A=hb$ is the area and $J_{22}=1/12b{h}^3$ is the second moment of area of the microbeam cross section about the ${\mathbf{e}}_2$-axis (see  Figure 1).

The four microbeams are characterized by increasing volume fractions (from left- to right-beam), respectively ${\varphi }_{\mathrm{C}}=\{2,3,5,7\}\%$. The lowest 13 eigenvalues and the corresponding eigenmodes are depicted in  Figure 2; the 1st mode with ${\omega }_1=4.7423$ since it just involves the translational motion of the shuttle mass with no bending of the beams. The higher modes can be classified as local modes encompassing only the flexural motion of each cantilever separately. This indeed makes it possible to use the device as a multimass detector. In particular the 2nd mode, the first flexural of the first microcantilever, is associated with a frequency which is almost twice that of the 1st mode.

A parametric analysis is initially performed in order to investigate the sensitivity of the eigenfrequencies with respect to the CNT content variation and beam length. In  Figure 3 the first five frequencies are computed by varying the CNT volume fraction separately for each beam. Each curve shows that starting from certain values of the volume fraction, the five eigenfrequencies gain a plateaux, thus indicating a level of CNT content that does not affect the frequency properties of the system. This occurs for the shuttle mass eigenfrequency (i.e., the first one of the SISO system), for small values of CNT volume fraction (about 0.5%), for any beam considered. The same behaviour can be observed for each beam eigenfrequency, although this occurs for values of volume fraction ${\varphi }_{\mathrm{C}}$ larger and larger. In the same  Figure 3 we also indicate through dashed grid lines the values of ${\varphi }_{\mathrm{C}}$ selected for the next computations, thus proving to have an adequate frequency gap of the system.

The influence of the beam length is also considered in  Figure 4: the natural frequency of a beam made by a nanocomposite is compared with the natural frequency of a beam made of silicon in terms of the ratio ${\omega }_{\mathrm{cnt}}/{\omega }_{\mathrm{sil}}$ between the 1st natural frequency of the nanocomposite beam and that of the silicon beam, respectively. We explore in  Figure 4 how this ratio varies either by increasing the beam nominal length of $270\mathsf{\mu }\mathrm{m}$ by an amplification factor $r\in [0,4]$, and by increasing the CNT volume fraction, starting from a nominal value of 2% then amplified again by r. While the natural frequency ratio monotically increases as the CNT volume fraction increases, by considering longer and longer lengths of CNT beam made of the same mixture, the natural frequency ratio tends to a plateaux. In particular, a nanocomposite beam with a 2.8% CNT volume fraction has the same value of the 1st natural frequency of a silicon beam of equal length.

If one chooses to set the first frequency of a microbeam to a certain target frequency, the microbeam could be designed in terms of a suitable length and CNT volume fraction.  Figure 5 shows indeed how to reduce the beam length as the CNT volume fractions increases. In particular, to have the same first frequency a silicon-made microbeam would be twice as long as a microbeam of nanocomposite with very low CNT volume fraction; a 8%-CNT nanocomposite beam has the same first frequency of a silicon microbeam of equal length.

4. Frequency Response Analysis

The Galerkin method is adopted to reduce the nonlinear integral-partial-differential equations of motion into a reduced-order model (ROM). By setting the unknown vector

$$\mathbf{v}{(s,t)}^{\top }=[x\left(t\right),v_1(s,t),v_2(s,t),\dots ,v_{N_b}(s,t)],$$

(27)

the full system of equation can be cast as

$$\mathbf{I}\ddot{\mathbf{v}}+\mathbf{L}\mathbf{v}=-\mathbf{n}(\mathbf{v},\dot{\mathbf{v}},\ddot{\mathbf{v}})+\mathbf{f},$$

(28)

where $\mathbf{I}$ and $\mathbf{L}$ are the generalized inertia and stiffness operators, $\mathbf{n}$ is the vector collecting the nonlinear terms, and $\mathbf{f}$ is the vector of external forces. By expressing the solution as

$$\mathbf{v}{(s,t)}^{\top }={[q_1\left(t\right)X_1,\sum _{j=2}^{N_m}{q}_j\left(t\right){\mathbf{u}}_j\left(s\right)]}^{\top },$$

(29)

where we considered the jth mode ${\mathbf{u}}_j^{\top }=[X_1,{\varphi }_{k,j}\left(s\right)]$, for $k=1,\dots N_b$, and substituting it into (28) yields

$${\ddot{q}}_j\mathbf{I}{\mathbf{u}}_j+q_j\mathbf{L}{\mathbf{u}}_j=-\mathbf{n}(q_j{\mathbf{u}}_j,{\dot{q}}_j{\mathbf{u}}_j,{\ddot{q}}_j{\mathbf{u}}_j)+\mathbf{f},$$

(30)

where the summation convention has been adopted. Multiplying by ${\mathbf{u}}_i^{\top }$ and integrating over $[0,1]$ yields the Galerkin discretized equations

$${\ddot{q}}_j{\int }_0^1{\mathbf{u}}_i^{\top }\mathbf{I}{\mathbf{u}}_{j}ds+q_j{\int }_0^1{\mathbf{u}}_i^{\top }\mathbf{L}{\mathbf{u}}_{j}ds=-{\int }_0^1{\mathbf{u}}_i^{\top }\mathbf{n}(q_j{\mathbf{u}}_j,{\dot{q}}_j{\mathbf{u}}_j,{\ddot{q}}_j{\mathbf{u}}_j)ds+{\int }_0^1{\mathbf{u}}_i^{\top }\mathbf{f}ds,$$

(31)

which can be rewritten as

$$M_{ij}{\ddot{q}}_j+K_{ij}{q}_j+n_i=f_{i}cos\mathsf{\Omega }t$$

(32)

with

$$n_i={\int }_0^1{\mathbf{u}}_i^{\top }\mathbf{n}\left(q_j{\mathbf{u}}_j\right)ds,f_i={\int }_0^1{\mathbf{u}}_i^{\top }\mathbf{f}ds.$$

(33)

By taking one trial function for each cantilever, a 5-dof overall ROM is obtained (i.e., 1 dof for the shuttle mass and 1 dof for each of the four microcantilevers). This approach was shown to be acceptable in [53,54] where a convergence analysis was carried out. We will also consider the parameters ${\zeta }_{\mathrm{M}}=0.01$ and $\zeta =0.011$ for the shuttle mass damping ratio and the microcantilevers damping ratio, respectively; all the other mechanical parameters are those reported in the previous section.

In order to obtain the periodic responses and the frequency response curves, we employed the adaptive pseudo-arclength pathfollowing procedure proposed in [54], developed by some of the authors and freely available on the Zenodo repository website: https://doi.org/10.5281/zenodo.6616482 (accessed on 8 May 2023). In particular, the search of periodic solutions is based on computing the fixed points of the Poincaré map whose Jacobian turns out to be the monodromy matrix, whose eigenvalues are the Floquet multipliers which dictate the stability of the periodic solutions [55,56]. The Jacobian matrix of the Poincaré map is calculated according to a finite difference approach in state space, and it works as iteration matrix. The interested reader can find further details in [54].

The dynamic behavior of the device is investigated by computing the frequency response curves for different device configurations, in particular, accounting for the case of additional masses $m_j$ at the tip of the jth microcantilever; unless otherwise stated, we set for $m_j^{*}=0.01$, i.e., 1% of the microcantilever mass ($\rho Al$ is about 835 pg). Finally, we took as characteristic frequency, the first bending frequency of a microbeam made of neat polymer, so that ${\omega }_c=345,755$ rad/s.

In the following, we present and discuss some numerical results aimed at demonstrating the sensitivity of the SISO system to the presence of additional tip masses. Note that in all figures reporting the frequency response curves of the SISO system, thin curves are referred to systems without tip masses, thick curves to systems with tip masses.

The frequency response curves reported in Figure 6, Figure 7 and Figure 8 are recovered in a frequency bandwidth close to the frequencies of the lowest bending modes of the microbeams. We tested four SISO configurations: a system without additional masses, a system where a tip mass is added to the first (from left) microbeam, another system with tip masses at both the first and the third microbeam, and finally a system with tip masses deposited on all microbeams. Such figures show that the frequency shift is bounded within a narrow frequency bandwidth close to the frequency peak of the microbeams with mass absorbed at the tip. The comparison also proves that the sensitivity increases with the frequency of the considered local vibration mode, and it is in line with similar results reported in the literature obtained for simple microstructures (e.g., microcantilevers or clamped-clamped microbeams), see [57,58].

Similar qualitative results can be found in Figure 9, Figure 10 and Figure 11 which report frequency response curves for the same four SISO configurations when the excitation frequency is increased until nearing to the frequency of the second bending mode of the microbeams. The nonlinearities of the system entail a higher sensitivity than that predicted by the linear model, thus proving that nonlinear effects are advantageous to ensure higher performance of the MEMS mass sensor system. In fact, 1% deposited mass produces a shift of the peak about 0.125 in nondimensional frequency around the first mode, while such a value is amplified about by ten times around the second mode.

Such a mass sensitivity is further investigated for increasing values of additional masses. We assume to add the tip mass just to the first (from left) microbeam. Figure 12 and Figure 13 show the frequency response curves limited to tip deflection of such a beam, with excitation frequency ranging in values around the first and the second mode, respectively.

Such figures not only reveal the expected trend of the resonance peak shift upon increasing the tip mass, but they also highlight a change in the length of the unstable branch of the response curve. In fact, around the first bending mode frequency, the shift in frequency occurs with a hardening nonlinear response which becomes more and more evident (see Figure 12); on the other hand, around the second bending mode frequency, the softening nonlinear response does not change much qualitatively (see Figure 13).

5. Conclusions

In this work the operating principles of a nonlinear MEMS multi-mass sensor are investigated. The device consists of an array of microcantilevers, made of carbon nanotube nanocomposite material, possessing a functionalized surface at their tips that can absorb molecules of target materials. The microcantilever sensors are clamped onto a shuttle mass which, in turn, is subject to electrostatic excitation. The system is treated as a single input-single output (SISO) system. Each microcantilever is modeled as a nonlinear Euler-Bernoulli beam undergoing a moderately large resonance motions.

The modal characteristics of the mechanical system are investigated carrying out a parametric study varying the CNT volume fraction in each microbeam. The frequency response and stability are investigated by path following the periodic solutions of the device near the lowest three bending modes of the microcantilevers.

The main finding is that the dynamic mass sensitivity measured as frequency shifts induced by the target mass absorbed at the cantilevers tip is enhanced by the nonlinear resonances of the microcantilevers modes, especially for higher modes. In particular, around the first bending mode frequency, we gain a frequency shift of the 12% in frequency by adding a 10% tip-mass, while for the second mode, the frequency shift is of 9% with 5% additional mass.

Although the experimental validation of the proposed model has not been considered, and future works will focus on this point, the accuracy of the nonlinear beam model for CNT-nanocomposites, adapted to the specific SISO system here theoretically designed, was widely validated in previous research contributions. In the present work, we showed how the fine tuning of the frequencies of the microcantilevers is positively effected by the carbon nanotube volume fraction in the nanocomposite material. Such innovative potentiality paves the way towards an effective design of multimass sensors since the microcantilevers may self-detect their motion exploiting the electrical conductivity of the CNT nanocomposites.