Analytical and experimental study of the dynamics of a micro-electromechanical resonator based digital-to-analog converter

Micro-beams are the most commonly used structures in micro-electromechanical systems (MEMS) resonators [1, 2]. In particular, clamped-clamped beams, a key resonator component, have been extensively explored for many practical applications such as in band-pass filters [3–5], oscillators, digital logic gates [6–11], and memory devices [12, 13]. The static and dynamic behaviors of resonant beams are complicated, especially under the nonlinear electrostatic-force actuation. Thus, it is necessary to understand the behavior of resonant beams in such cases.

Considerable effort has been made to understand the complex static and dynamic behaviors of MEMS structures, especially beams. A review of recent advances in this field is presented in [14]. Mestrom et al [15] studied the softening and hardening behaviors of a clamped-clamped beam resonator through a combined analytical-numerical and experimental approach. The obtained results show softening or hardening nonlinear dynamics depending on excitations. Mestrom et al [16] also utilized the clamped-clamped beam resonator to predict the measured resonator response and the nonlinear dynamic steady-state behaviors under different parameters. A good quantitative match was achieved between simulation and experimental results. Venstra et al [17] studied the clamped-clamped mechanical resonators in water through a magnetomotive drive. The results demonstrate that the magnetomotive technique can be used to drive and detect micromechanical resonator vibrations. Azizi et al [18] investigated the nonlinear chaotic dynamics of an electrostatically and piezo-electrically actuated clamped-clamped micro beam. The results indicate that an appropriate DC source to piezoelectric layer can passively control the chaotic response of the micro system. Hu et al [19] investigated the snap-through, pull-in behaviors and dynamics of an initially curved micro clamped-clamped beam under an electro-dynamical actuation. Kumar et al [20] studied the nonlinear dynamics and the internal resonances of a clamped-clamped beam MEMS resonator by using Galerkin based reduced-order model and a finite element method. Their results can serve as guidelines for developing controllers that regulate the displacement of driven mode by its interaction with the coupled mode. Ghayesh et al [21] analytically presented an investigation on an electrically actuated MEMS clamped-clamped resonator and studied the pull-in instability, force-response, and dynamic behavior of the system through the time histories and phase-plane portraits. Li et al [22] studied the nonlinear modal interactions of electrically actuated clamped-clamped micro beams and the effects of antisymmetric modes on the nonlinear dynamics. The results show that the nonlinear modal interactions can transfer the energy from one mode to a nearby mode.

As mentioned above, the Galerkin method has been widely used to investigate the dynamics of MEMS structures. However, under low damping conditions and for complex electrostatic actuation, the Galerkin method suffers since it can only catch part of the dynamic state (solution). Furthermore, simulations using this method are time-consuming and do not converge easily. To alleviate this, a combined Galerkin-shooting technique can be used [1]. Most resonance curves, especially near the bifurcation point [1, 23], can be captured by this powerful method. Nayfeh et al [24] and Younis et al [25] adopted a shooting method along with the Galerkin method to study the nonlinear dynamics of a novel MEMS switch excited around the primary resonance and undergoing large motion. Yu et al [26] adopted the combined Galerkin-shooting technique to analyze the electromechanical post-buckling response of a clamped-clamped beam actuator subjected to axial residual stresses and a symmetric electrostatic field. Using same technique, Emam et al [27] experimentally and analytically investigated the nonlinear behaviors of a clamped-clamped buckled beam near subharmonic resonance and investigate the stability and bifurcations based on a combined method of Galerkin-shooting. Nayfeh and Younis [28] utilized the combination of the shooting technique and a reduced-order model to analyze and simulate the dynamics of an electrically actuated micro-beam, near primary, sub-harmonic, and super-harmonic resonances.

Digital-to-analog converters (DACs) are basic elements for electronics and microcontrollers in MEMS. Recent, efforts have been made to realize DACs using MEMS technology motivated by the need to reduce power consumption, which has a more promising digital signal processing application in MEMS devices [29]. Yeh et al [30] demonstrated a technique to create a 3-bit mechanical DAC by using lever arms and thermal actuator arrays. It can take digital electrical signals and produce mechanical displacements at the output which can be used in applications such as micro-optics, mechanical-computing, and micro-robotics. Toshiyoshi et al [31] discussed the design and implementation of a 4 bit micromachined electrostatic actuator-based MEMS DAC in which a set of local digital displacements are converted into an analog output displacement. Pandiyan et al [32] used a bending beam electro thermal compliant actuator to design and simulate an MEMS-based DAC. The proposed DAC generates a mechanical displacement at the output from the digital input based on the principle of weighted stiffness which is analogous to binary-weighted resistor DAC.

The investigated DACs above based on different micro-structures are complicated consisting of arrays of many asymmetric and symmetric structures with multi-bent beam actuators. As a result, these kinds of reported DACs structures occupy large areas. Also, some of the reported DACs need DC currents for their operation, which results in high energy consumption. The proposed DAC here consists of a single straight micro clamped-clamed beam with multi-partial electrodes, which is simple and easy to fabricate. Also, this DAC does not need a DC current, hence it consumes less power. The proposed DAC can be used to internet of things (IoT) applications that demand moderate speed and lower power consumption.

In this paper, we analytically and experimentally investigate the statics and dynamics of the MEMS resonator-based 4-bit DAC [33, 34] in depth. The device consists of a clamped-clamped beam, separated from various partial electrodes of different air-gap sizes. The digital inputs are DC voltages, applied to four partial electrodes, which tune the beam's resonance frequency using the electrostatic softening effect. Our objective is to develop a theoretical model, predict all the relevant aspects of the experimental response, and optimize the device performance. First, we derive the dynamic equations of the device, and then we investigate the static deflection, eigenvalue problem (EVP), optimization, and dynamic response of the micro-resonator based on the Galerkin-shooting procedure. Lastly, we examine two measurement samples of the DAC under the nonlinear operation. The forward and backward frequency sweep on each digital case is conducted, and excellent agreements are found between the experimental and theoretical results.

The outline of the rest of the paper is as follows. Section 2 describes the modeling procedure. The static response, EVP, and optimization are presented in sections 3 and 4, respectively. Section 5 experimentally and analytically illustrates the frequency response for the forward and backward frequency sweeps and examines the effects of the biasing voltage of the beam on the DAC. Finally, conclusions remarks are summarized in the last section.

A picture of the proposed resonator DAC is shown in figure 1(a). The device is fabricated by MEMSCAP Inc foundry through SOIMUMPs process [35] using a surface micromachining process on a silicon on an insulator wafer. The process starts by defining electrodes with gold by lifting off a Cr/Au bi-layer (20 nm/500 nm), and then the silicon device layer (25 µm thick) is patterned by deep reactive ion etching to define the beam, the anchors, and corner electrodes. Finally, the structure is released by etching the SiO₂ in using vapor Hydrogen Fluoride (HF). The used wafer has the following layer thicknesses: (a) silicon thickness: 25 ± 1 µm; (b) oxide thickness: 1 ± 0.05 µm; (c) handle wafer (substrate) thickness: 400 ± 5 µm.

Zoom In Zoom Out Reset image size

The DAC device consists of an in-plane-clamped-clamped micro beam based resonator (length l= 500 µm, lateral width W= 2.66 µm, thickness h= 25 µm) with a DC voltage (V_Beam), and side partial electrodes with various air-gap widths as shown in figure 1(b). The digital inputs of the DAC are DC voltages, applied to four corner electrodes marked as D₃, D₂, D₁, and D₀ in figure 1(a). The actuation of the micro beam resonator is realized through the different digital input combinations (DC source) and AC voltage components. The resonator here is driven to vibrate on the first vibration mode.

The input vector D₃ D₂ D₁ D₀ is a binary number where D₃ (with the smallest gap) is the most significant bit and D₀ (with the largest gap) is the least significant bit. Each bit can be either '1' (on) or '0' (off), which are defined as '40' V or '0' V, respectively. Thus, there are 16 combinations (2⁴ = 16) for the inputs ranging from '0000' to '1111' representing the decimal numbers from 0 to 15. The resonance frequency of the beam should increase as the decimal equivalent value of the applied digital input increases.

The equations of motion governing the transverse deflection $\hat w$ of a clamped-clamped beam under the partial electrode actuation can be described as

$$\begin{align} EI{\text{ }}\hat w{^{{\prime} {\prime}}}{^{{\prime} {\prime}}} & + \rho A{\text{ }}\ddot{\hat{w}} + \hat c{\text{ }}\dot{\hat{w}}\; \nonumber \\ & = \;\left[\hat N + \frac{{EA}}{{2l}}\int_0^l {{{(\hat w(\hat x,\hat t)^{\prime} )}^2}} {\text{d}}\hat x \right]{\text{ }}\hat w{^{{\prime} {\prime}}}{\text{ }} \nonumber \\ & \quad + \left[\frac{1}{2}\varepsilon b\frac{{{{({V_{{\text{D}}1}} - {V_{{\text{Beam}}}})}^2}}}{{{{({d_1} - \hat w)}^2}}} - \frac{1}{2}\varepsilon b\frac{{{{({V_{{\text{D}}0}} - {V_{{\text{Beam}}}})}^2}}}{{{{({d_0} + \hat w)}^2}}}\right] \nonumber \\ &\quad\quad [U({{\hat a}_5} - \hat x) - U({{\hat a}_4} - \hat x)] \nonumber \\ & \quad + \left[\frac{1}{2}\varepsilon b\frac{{{{({V_{{\text{Beam}}}} - {V_{{\text{AC}}}}\cos (\Omega \hat t))}^2}}}{{{{(d - \hat w)}^2}}} - \frac{1}{2}\varepsilon b\frac{{{V_{{\text{Beam}}}}^2}}{{{{(d + \hat w)}^2}}}\right] \nonumber \\ &\quad\quad [U({{\hat a}_3} - \hat x) - U({{\hat a}_2} - \hat x)]\; \nonumber \\ & \quad + \left[\frac{1}{2}\varepsilon b\frac{{{{({V_{{\text{D3}}}} - {V_{{\text{Beam}}{} }})}^2}}}{{{{({d_3} - \hat w)}^2}}} - \frac{1}{2}\varepsilon b\frac{{{{({V_{{\text{D2}}}} - {V_{{\text{Beam}}}})}^2}}}{{{{({d_2} + \hat w)}^2}}}\right] \nonumber \\ &\quad\quad [U({{\hat a}_1} - \hat x) - U({{\hat a}_0} - \hat x)]] \end{align} \tag{ 1 }$$

with the following boundary conditions:

$$\begin{equation}\hat w(0,\hat t) = \hat w(l,\hat t) = 0{\text{ }}{\text{and}}{\text{ }}{\hat{w}}{^{^{\prime}}}(0,\hat t) = {\text{ }}\hat w{^{^{\prime}}}(l,\hat t) = 0\end{equation} \tag{ 2 }$$

where the prime indicates the derivative in space ${\hat w^{\prime} }(\hat x,\hat t) = \frac{{\partial \hat w}}{{\partial \hat x}}$ and the dot indicates the derivative in time $\dot{\hat{w}}(\hat x,\hat t) = \frac{{\partial \hat w}}{{\partial \hat t}}$.

Here, x is the position along the beam length, E and I are the elastic modulus and the moment of inertia, b is the width of the beam, h is the beam thickness, $\hat c$ is the viscous damping coefficient, $\hat N$ is the axial load, A is the cross-section area of the beam, l is the length of the beam, d is the gap between the beam and drive or sense electrode. The air gaps between the corner electrodes and beam are d₀, d₁, d₂, and d₃, respectively. Also, U(x) is a step function. The distance of each partial electrode position to the left starting point (o) of the resonator in figure 1(a) is ${\hat a_i}$ (i = 0, 1, 2, 3, 4, 5) respectively. The detailed values are illustrated in table 1.

Table 1. The specifications of the resonator.

Symbol	Value	Quantity	Value
d0	8.74 µm	${\hat a_0}$	8 µm
d1	6.24 µm	${\hat a_1}$	163.9 µm
d2	5.41 µm	${\hat a_2}$	171.9 µm
d3	3.33 µm	${\hat a_3}$	327.8 µm
d	7.075 µm	${\hat a_4}$	335.8 µm
l	500 µm	${\hat a_5}$	491.7 µm
h	25 µm	W	2.66 µm

Equation (1) is governed by 16 key parameters (d₁, d₂, d₃, d₄, d, and ${\hat a_0},{\hat a_1},{\hat a_2},{\hat a_3},{\hat a_4},{\hat a_5}$, V_D0, V_D1, V_D2, V_D3, V_Beam). To facilitate analysis of the equation, the following dimensionless variables are introduced

$$\begin{equation}w = \frac{{\hat w}}{{{d_3}}};\quad x = \frac{{\hat x}}{l};\quad t = \hat t\sqrt {\frac{{\rho bh{l^4}}}{{EI}}} ;\quad {a_i} = \frac{{{{\hat a}_i}}}{l}.\end{equation} \tag{ 3 }$$

We substitute equation (3) into equations (1) and (2), yielding the non-dimensional equation

$$\begin{align} w{^{^{\prime}}}{^{^{\prime}}}{^{^{\prime}}}{^{^{\prime}}} & + \ddot w + c{\text{ }}\dot w \nonumber \\ = & \left[{N_{{\text{non}}}} + {\alpha _1}\int_0^1 {{{(w(x,t)^{\prime} )}^2}} {\text{d}}x\right]{\text{ }}w{^{^{\prime}}}{^{^{\prime}}} \nonumber \\ & + {\alpha _2} \left[\frac{{{{({V_{{\text{D}}1}} - {V_{{\text{Beam}}}})}^2}}}{{{{(\frac{{{d_1}}}{{{d_3}}} - w)}^2}}} - \frac{{{{({V_{{\text{D}}0}} - {V_{{\text{Beam}}}})}^2}}}{{{{(\frac{{{d_0}}}{{{d_3}}} + w)}^2}}}\right] \nonumber \\ &\quad\quad [U({a_5} - x) - U({a_4} - x)] \nonumber \\ & + {\alpha _2} \left[\frac{{{{({V_{{\text{Beam}}}} - {V_{{\text{AC}}}}\cos (\Omega t))}^2}}}{{{{(\frac{d}{{{d_3}}} - w)}^2}}} - \frac{{{V_{{\text{Beam}}}}^2}}{{{{(\frac{d}{{{d_3}}} + w)}^2}}}\right] \nonumber \\ &\quad\quad [U({a_3} - x) - U({a_2} - x)] \nonumber \\ & + {\alpha _2} \left[\frac{{{{({V_{{\text{D3}}}} - {V_{{\text{Beam}}\,}})}^2}}}{{{{(1 - w)}^2}}} - \frac{{{{({V_{{\text{D}}2}} - {V_{{\text{Beam}}}})}^2}}}{{{{(\frac{{{d_2}}}{{{d_3}}} + w)}^2}}}\right] \nonumber \\ &\quad\quad [U({a_1} - x) - U({a_0} - x)]]. \nonumber \\ \end{align} \tag{ 4 }$$

The above is subjected to the non-dimensional boundary conditions:

$$\begin{equation}w(0,t) = w(1,t) = 0;\quad w{^{^{\prime}}}(0,t) = w{^{^{\prime}}}(1,t) = 0.\end{equation} \tag{ 5 }$$

The non-dimensional parameters appearing in equation (4) are

$$\begin{equation}{\alpha _1} = 6{\left( {\frac{{{d_3}}}{b}} \right)^2};{\alpha _{\text{2}}} = \frac{{6\varepsilon {l^4}}}{{E{b^3}{d_3}^3}};{N_{{\text{non}}}} = \frac{{{l^2}}}{{EI}}\hat N.\end{equation} \tag{ 6 }$$

The deflection of the clamped-clamped micro beam under the electrostatic force is expressed as a static normalized deflection w_s(x) due to the DC voltage component and a normalized dynamic component due to the AC source, denoted by w_d(x, t); that is

$$\begin{equation}w\left( {x,t} \right){ } = {w_{\text{s}}}\left( x \right) + {w_{\text{d}}}\left( {x,t} \right).\end{equation} \tag{ 7 }$$

To determine w_s(x), we set the time derivative and the AC forcing terms in equation (4) equal zero and obtain

$$\begin{align} {w_{\text{s}}}(x){^{^{\prime}}}{^{^{\prime}}}{^{^{\prime}}}{^{^{\prime}}} & - \left[{N_{{\text{non}}}} + {\alpha _1}\int_0^1 {{{({w_{\text{s}}}(x){^{^{\prime}}})}^2}} {\text{d}}x]{w_{\text{s}}}(x){^{^{\prime}}}{^{^{\prime}}}\right] \nonumber \\ & = + {\alpha _2}\left[\frac{{{{({V_{{\text{D1}}}} - {V_{{\text{Beam}}}})}^2}}}{{{{(\frac{{{d_1}}}{{{d_3}}} - {w_{\text{s}}})}^2}}} - \frac{{{{({V_{{\text{D0}}}} - {V_{{\text{Beam}}}})}^2}}}{{{{(\frac{{{d_0}}}{{{d_3}}} + {w_{\text{s}}})}^2}}}\right] \nonumber \\ & \quad [U({a_5} - x) - U({a_4} - x)] \nonumber \\ & + {\alpha _2}\left[\frac{{{V_{{\text{Beam}}}}^2}}{{{{(\frac{d}{{{d_3}}} - {w_{\text{s}}})}^2}}} - \frac{{{V_{{\text{Beam}}}}^2}}{{{{(\frac{d}{{{d_3}}} + {w_{\text{s}}})}^2}}}\right] \nonumber \\ & \quad [U({a_3} - x) - U({a_2} - x)] \nonumber \\ & + {\alpha _2}\left[\frac{{{{({V_{{\text{D}}3}} - {V_{{\text{Beam}}\,}})}^2}}}{{{{(1 - {w_{\text{s}}})}^2}}} - \frac{{{{({V_{{\text{D2}}}} - {V_{{\text{Beam}}}})}^2}}}{{{{(\frac{{{d_2}}}{{{d_3}}} + {w_{\text{s}}})}^2}}}\right] \nonumber \\ & \quad [U({a_1} - x) - U({a_0} - x)]] \end{align} \tag{ 8 }$$

with the associated boundary conditions

$$\begin{equation}{w_{\text{s}}}(0,t) = {w_{\text{s}}}(1,t) = 0;{w_{\text{s}}}^{\prime} (0,t) = {w_{\text{s}}}^{\prime} (1,t) = 0.\end{equation} \tag{ 9 }$$

By plugging equation (7) into equation (4), and only keeping linear terms of w_d after the electrostatic force is expanded in w_d by Taylor series expansions, we obtain

$$\begin{align} & {w_{\text{d}}}(x,t){^{^{\prime}}}{^{^{\prime}}}{^{^{\prime}}}{^{^{\prime}}} + {w_{\text{s}}}(x,t){^{^{\prime}}}{^{^{\prime}}}{^{^{\prime}}}{^{^{\prime}}} + {{\ddot w}_{\text{d}}}(x,t) + c{{\dot w}_{\text{d}}}(x,t) \nonumber \\ & = \left[{N_{{\text{non}}}} + {\alpha _1}\int_0^l {{{({w_{\text{d}}}(x,t){^{^{\prime}}} + {w_{\text{s}}}(x,t){^{^{\prime}}})}^2}} {\text{d}}x\right] \nonumber \\ &\quad\quad ({w_{\text{d}}}(x,t){^{^{\prime}}}{^{^{\prime}}} + {w_{\text{s}}}(x,t){^{^{\prime}}}{^{^{\prime}}}) \nonumber \\ & \quad + {\alpha _2}\left[{({V_{{\text{D1}}}} - {V_{{\text{Beam}}}})^2}{\text{ }}\left[\frac{1}{{{{(\frac{{{d_1}}}{{{d_3}}} - {w_{\text{s}}})}^2}}} + \frac{{2{w_{\text{d}}}}}{{{{(\frac{{{d_1}}}{{{d_3}}} - {w_{\text{s}}})}^3}}}\right]\right. \nonumber \\ &\quad\quad \left. - {({V_{{\text{D0}}}} - {V_{{\text{Beam}}}})^2}\left[\frac{1}{{{{(\frac{{{d_0}}}{{{d_3}}} + {w_{\text{s}}})}^2}}} - \frac{{2{w_{\text{d}}}}}{{{{(\frac{{{d_0}}}{{{d_3}}} + {w_{\text{s}}})}^3}}}\right]\right] \nonumber \\ &\quad\quad\quad [U({a_5} - x) - U({a_4} - x)] \nonumber \\ & \quad + {\alpha _2}\left[{({V_{{\text{Beam}}}} + {V_{{\text{AC}}}}\cos (\Omega t))^2}\left[\frac{1}{{{{(\frac{d}{{{d_3}}} - {w_{\text{s}}})}^2}}} + \frac{{2{w_{\text{d}}}}}{{{{(\frac{d}{{{d_3}}} - {w_{\text{s}}})}^3}}}\right]\right. \nonumber \\ &\quad\quad \left.- {V_{{\text{Beam}}}}^2{\text{ }}\left[\frac{1}{{{{(\frac{d}{{{d_3}}} + {w_{\text{s}}})}^2}}} - \frac{{2{w_{\text{d}}}}}{{{{(\frac{d}{{{d_3}}} + {w_{\text{s}}})}^3}}}\right]\right] \nonumber \\ &\quad\quad\quad [U({a_3} - x) - U({a_2} - x)] \nonumber \\ & \quad + {\alpha _2}\left[{({V_{{\text{D3}}}} - {V_{{\text{Beam}}\,}})^2}\left[\frac{1}{{{{(1 - {w_{\text{s}}})}^2}}} + \frac{{2{w_{\text{d}}}}}{{{{(1 - {w_{\text{s}}})}^3}}}\right]\right. \nonumber \\ &\quad\quad \left. - {({V_{{\text{D2}}}} - {V_{{\text{Beam}}}})^2}\left[\frac{1}{{{{(\frac{{{d_2}}}{{{d_3}}} + {w_{\text{s}}})}^2}}} - \frac{{2{w_{\text{d}}}}}{{{{(\frac{{{d_2}}}{{{d_3}}} + {w_{\text{s}}})}^3}}}\right]\right] \nonumber \\ &\quad\quad\quad [U({a_1} - x) - U({a_0} - x)]]. \nonumber \\ \end{align} \tag{ 10 }$$

By eliminating the terms representing the equilibrium position of equation (8) and the AC term from equation (10), this yields:

$$\begin{align} {w_{\text{d}}}(x,t){^{^{\prime}}}{^{^{\prime}}}{^{^{\prime}}}{^{^{\prime}}} & + {{\ddot w}_{\text{d}}}(x,t) + c{{\dot w}_{\text{d}}}(x,t) \nonumber \\ = & [{N_{{\text{non}}}} + {\alpha _1}\int_0^l {{{({w_{\text{s}}}(x){^{^{\prime}}})}^2}} {\text{d}}x]{w_{\text{d}}}(x,t){^{^{\prime}}}{^{^{\prime}}} \nonumber \\ & + {\alpha _1}\left[\int_0^l {({{({w_{\text{d}}}(x,t){^{^{\prime}}})}^2} + 2{w_{\text{d}}}(x,t){^{^{\prime}}}{w_{\text{s}}}(x){^{^{\prime}}}} ){\text{d}}x\right]\nonumber \\ & [{w_{\text{d}}}(x,t){^{^{\prime}}}{^{^{\prime}}} + {w_{\text{s}}}(x){^{^{\prime}}}{^{^{\prime}}}] \nonumber \\ & + {\alpha _2}\left[\frac{{2{{({V_{{\text{D1}}}} - {V_{{\text{Beam}}}})}^2}}}{{{{(\frac{{{d_1}}}{{{d_3}}} - {w_{\text{s}}})}^3}}}\right] + \frac{{2{{({V_{{\text{D0}}}} - {V_{{\text{Beam}}}})}^2}}}{{{{(\frac{{{d_0}}}{{{d_3}}} + {w_{\text{s}}})}^3}}}]{w_{\text{d}}}\nonumber \\ & [U({a_5} - x) - U({a_4} - x)] \nonumber \\ & + {\alpha _2}\left[\frac{{2V_{{\text{Beam}}}^{\text{2}}}}{{{{(\frac{d}{{{d_3}}} - {w_{\text{s}}})}^3}}} + \frac{{2{V_{{\text{Beam}}}}^2{\text{ }}}}{{{{(\frac{d}{{{d_3}}} + {w_{\text{s}}})}^3}}}\right]{w_{\text{d}}} \nonumber \\ & [U({a_3} - x) - U({a_2} - x)] \nonumber \\ & + {\alpha _2}\left[\frac{{2{{({V_{{\text{D}}3}} - {V_{{\text{Beam}}\,}})}^2}}}{{{{(1 - {w_{\text{s}}})}^3}}} + \frac{{2{{({V_{{\text{D2}}}} - {V_{{\text{Beam}}}})}^2}}}{{{{(\frac{{{d_2}}}{{{d_3}}} + {w_{\text{s}}})}^3}}}\right]{w_{\text{d}}} \nonumber \\ & [U({a_1} - x) - U({a_0} - x)]] \nonumber \\ \end{align} \tag{ 11 }$$

with the associated boundary condition

$$\begin{equation}{w_{\text{d}}}(0,t) = {w_{\text{d}}}(1,t) = 0;{w_{\text{d}}}^{\prime} (0,t) = {w_{\text{d}}}^{\prime} (1,t) = 0.\end{equation} \tag{ 12 }$$

In the following analysis, we only consider the 3-bit DAC while keeping the D₀ un-actuated to compare the analytical result with the experimental data reported in [33, 34]. The air gaps between the beam and the corner electrodes are adjusted based on the weight of each bit such that the more significant bits have stronger effects on the beam. This is realized by making the gap of D₀ the largest while making the gap of D₃ the smallest.

Next, we resort to the Galerkin procedure [1] to solve for the EVP and dynamic responses of the system by assuming the micro-beam deflection as

$$\begin{equation}w(x,t) = \sum\limits_{i = 1}^N {{u_i}(t){\phi _i}} (x)\quad \end{equation} \tag{ 13 }$$

where ϕ_i (x) is the ith mode shape of the micro-beam at zero voltage load and u_j is its modal amplitude. Equation (13) is then substituted into equation (11) for the EVP, and equation (4) for the forced vibration response, then we multiply the outcome by the mode shape, apply the orthogonality condition of the mode shapes, and integrate over the beam domain from 0 to l. A convergence study needs to be conducted to determine the minimum number of modes that are needed to capture the dynamic behavior accurately [1]. In this work, for simplicity, we use a single-mode to approximate the dynamic response. The obtained analytical results reach good agreement with the experimental ones. Thus, a single-mode yields sufficient accuracy.

Here, we examine the static response to better understand how static deflection can impact the system behavior. We explore and simulate the static response for each digital combination cases. For the static analysis, due to the complicated denominator of forcing terms in equation (8), the solution does not easily converge and the adoption of the Galerkin method is time-consuming. To improve the efficiency of simulations, and obtain more accurate results, we resort to a shooting method to solve such equations [36, 37].

Next, we solve the static response based on the shooting method with the boundary value problem. We shoot for the value of the complicated integral mid-plane stretching term of equation (8) iteratively by the Newton Raphson method until achieving convergence to within a very small predefined tolerance. Figure 2 shows the static deflections along the micro beam for eight different digital input combination cases from case-000 to case-111.

Zoom In Zoom Out Reset image size

From figure 2, the static deflection curve and the maximum deflection point of the micro beam for different digital input combinations are illustrated. We conclude that the static deflection curves vary for different digital input combinations. Also, we observe that the multi-actuated electrode configuration for the digital input can shift the maximum deflection along the micro beam. This is due to the fact that the generated electrostatic force by different digital combinations varies along the beam, and also due to the digital electrodes with different gap widths. Therefore, the maximum position along the micro beam varies with the digital input combinations.

For example, the digital case-010 (the bits D3 and D1 are off while the D2 is on) exhibits the maximum static deflection. The electrode D3 with the smallest air-gap generates the largest electrostatic force, and then the partial electrodes D3 and D1 on the same side of the beam work together to yield the higher amplitude response. However, for the digital case-111, the electrostatic force generated by the connected corner partial electrodes (D3 D2 D1) neutralize the ones from the micro beam, thus, almost zero amplitude is generated in this case. To examine case-100, one can see that D3 is digitally connected while the D2 and D1 are not connected, which means the electrostatic force generated by electrodes D2 and D3 works on the beam. But these forces mainly concentrate on the area around the anchor of either side of the micro beam, which bends the micro beam. Thus, the static deflection for case-100 is similar to a sinusoidal function.

We notice that, compared to the smallest air gaps (d₃ = 3.33 µm), the maximum deflection for every single digital case is very weak, particularly for the digital case-111. Thus, in other words, the static deflection of the micro beam at 40 V is almost zero which is negligible. As a result, we will neglect the static deflection w_s (x) when solving the EVP in the following section.

As pointed out previously, the static deflection is very low when digital partial electrodes 40 V are connected. Hence, to ease the computational cost, we solve the EVP in equations (11) and (12) without the damping term, while also neglecting the static deflection w_s(x). We apply a single-mode in the Galerkin model, and convert the partial differential equation into an ordinary differential equation. Then, we collect the coefficients of the linear term u(t); and take its square root, which yields the following expression for the first natural (resonance) frequency of the beam under electrostatic actuation:

$$\begin{align} {\small\begin{array}{l} {\text{Freq }}({\text{kHz}}) \\[8pt] \quad \displaystyle = \frac{{\sqrt { - 0.00684568{ }[\frac{{{{({V_{{\text{D}}0}} - {V_{{\text{Beam}}}})}^2}{d_3}^3}}{{{d_0}^3}} + \frac{{{{({V_{{\text{D}}1}} - {V_{{\text{Beam}}}})}^2}{d_3}^3}}{{{d_1}^3}}] - 0.00679926(\frac{{{{({V_{{\text{D}}2}} - {V_{{\text{Beam}}}})}^2}{d_3}^3}}{{{d_2}^3}} + {{({V_{{\text{D}}3}} - {V_{{\text{Beam}}}})}^2}) - \frac{{0.072397{V_{{\text{Beam}}}}^2{d_3}^3}}{{{d^3}}} + {{20.5004}^2}} }}{{2\pi T \times 1000}}.\end{array}}\end{align} \tag{ 14 }$$

Note that the resonance frequency of the resonator is a complicated function of eight key parameters. Based on equation (14), the resonance frequency of a resonator can be further analyzed and optimized.

To study the operation of the DAC in the linear domain, we experimentally test the DAC device under a large vacuum level. Figure 3(a) shows the measured frequency curves for different eight input combinations. For the resonator to work as a DAC, the frequency of the drive signal should be fixed at a selected operating frequency, for example, f₁ or f₂, and the output level should decrease or increase when a digital input is applied. Figure 3(b) shows the time response when the drive frequency is fixed around f₁. The resonance frequencies extracted from the experiment, figure 3(a), are plotted in figure 3(c), in which the frequency shifts range from 81.807 kHz for case-000 to 82.670 kHz for case-111. Based on the derived equation (14), the analytical resonance frequencies of the resonator for different input combinations are shown in figure 3(d). One can note good agreement between the analytical and experimental results.

Zoom In Zoom Out Reset image size

However, it is noticed that the experimentally obtained frequency response curves in figure 3(a) indicate considerable overlaps for the output under different combinations, which makes the amplitude levels created by different inputs at f₁ or f₂ very close and undistinguishable. Also, the range of the frequency-shifting for different electrode combinations is not large. It is difficult to separate each frequency curve, which is a basic requirement for the proper operation of the DAC [34]. Here, there are two explanations. One is the voltage generated by partial electrodes can partly neutralize the one from the micro beam. The effective power source to trigger the micro resonator, as a result, will be reduced if more partial electrodes are joined. Another one is that the four corner electrodes close to the anchor of the micro beam have less influence on its actuation/deformation.

To overcome this, one can optimize the sizes of air gaps to generate various electrostatic forces along the micro beam under different digital combination cases. By changing air gaps, the distribution of electrostatic force along the micro beam can be reallocated. The expected frequency shift for each case can be improved by optimizing the air gap sizes.

To optimize the air gaps, we introduce d_s as a basic scaling parameter, valued at d_s = 5 μm, for normalization ($w = \hat w/{d_{\text{s}}}$). This process will enable us to optimize all the gaps including d₁, d₂, d_3, and d. It is worthy to note that the new normalization-based d_s here is mainly responsible for the optimization part. Applying the same normalization procedure to the obtained equation (14), the resonance frequency can be rewritten as below:

$$\begin{align}& {\text{Freq}}({\text{kHz}}) \nonumber \\ & = \frac{{\sqrt { - 0.25107 \times {{10}^{ - 18}}(\frac{{{{({V_{{\text{D}}2}} - {V_{{\text{Beam}}}})}^2}}}{{{d_2}^3}} + \frac{{{{({V_{{\text{D}}3}} - {V_{{\text{Beam}}}})}^2}}}{{{d_3}^3}}) - 0.252784 \times {{10}^{ - 18}}{ }\frac{{{{({V_{{\text{D}}1}} - {V_{{\text{Beam}}}})}^2}}}{{{d_1}^3}} - 0.26823 \times {{10}^{ - 18}}\frac{{{V_{{\text{Beam}}}}^2}}{{{d^3}}} + {{20.5004}^2}} }}{{2\pi T \times 1000}}.\end{align} \tag{ 15 }$$

Here, we resort to multi-objective goal attainment problems [38] of all the gaps including d₁, d₂, d₃, and d. We optimize the gaps over the size range from 2 µm to 8 µm. Next, we define a target function and restricted conditions. The restrictions need to satisfy the following conditions: (a) the difference between each neighboring digital frequency is larger than 0.3 kHz in a specific domain; (b) the frequency difference between the last case (111) and first case (000) is not too far away; (c) and the frequency linearly and monotonically increase between each case. We optimize the frequency based on the equation (15). The outcome of the optimized gaps is d₁ = 5.47 µm, d₂ = 4.53 µm, d₃ = 3.58 µm, d= 4.84 µm. The simulation results based on these values are shown in figure 4. As noted in figure 4, the resonance frequencies for each digital case are linearly ordered.

Zoom In Zoom Out Reset image size

To further study the forced vibration response of the micro-beam, we apply the procedure of the Galerkin discretization with a single-mode. First, we substitute equation (13) into equation (4), multiply both sides of the outcome by ϕ_j, and then we integrate the results from the beam domain zero to one. The outcome yields:

$$\begin{align} {{\ddot u}_1}&(t) + c{{\dot u}_1}(t) + \omega _{{\text{non}}}^2{u_1}(t) \nonumber \\ & = {\alpha _1}\int_0^1 {{{({{\phi {^{^{\prime}}}}_1}(x){u_1}(t))}^2}} {\text{d}}x{\text{ }}\int_0^1 {{\phi _1}} (x){{\phi {^{{\prime} {\prime}}}}_1}(x){\text{d}}x{\text{ }}{u_1}(t) \nonumber \\ & \quad + {\alpha _2}[{({V_{{\text{D}}1}} - {V_{{\text{Beam}}}})^2}\int_0^1 {\frac{{{\phi _1}(x)}}{{{{(\frac{{{d_1}}}{{{d_3}}} - {\phi _1}(x){u_1}(t))}^2}}}} \nonumber \\ & \quad\quad [U({a_5} - x) - U({a_4} - x)]{\text{d}}x - {({V_{{\text{D}}0}} - {V_{{\text{Beam}}}})^2} \nonumber \\ & \quad\quad \int_0^1 {\frac{{{\phi _1}(x)}}{{{{(\frac{{{d_0}}}{{{d_3}}} + {\phi _1}(x){u_1}(t))}^2}}}} [U({a_5} - x) - U({a_4} - x)]{\text{d}}x \nonumber \\ & \quad + {\alpha _2}[{({V_{{\text{Beam}}}} - {V_{{\text{ac}}}}\cos (\Omega {\text{ }}t))^2}\int_0^1 {\frac{{{\phi _1}(x)}}{{{{(\frac{d}{{{d_3}}} - {\phi _1}(x){u_1}(t))}^2}}}} \nonumber \\ & \quad\quad [U({a_3} - x) - U({a_2} - x)]{\text{d}}x - {V_{{\text{Beam}}}}^2 \nonumber \\ & \quad\quad \int_0^1 {\frac{{{\phi _1}(x)}}{{{{(\frac{d}{{{d_3}}} + {\phi _1}(x){u_1}(t))}^2}}}} [U({a_3} - x) - U({a_2} - x)]{\text{d}}x \nonumber \\ & \quad + {\alpha _2}[{({V_{{\text{D3}}}} - {V_{{\text{Beam}}\,}})^2}\int_0^1 {\frac{{{\phi _1}(x)}}{{{{(1 - {\phi _1}(x){u_1}(t))}^2}}}} \nonumber \\ & \quad\quad [U({a_1} - x) - U({a_0} - x)]{\text{d}}x - {({V_{{\text{D2}}}} - {V_{{\text{Beam}}}})^2} \nonumber \\ & \quad \,\,\,\,\int_0^1 \!\!{\frac{{{\phi _1}(x)}}{{{{(\frac{{{d_2}}}{{{d_3}}} \!+ \!{\phi _1}(x){u_1}(t))}^2}}}} [U({a_1} - x) \!- \!U({a_0} - x)]{\text{d}}x. \end{align} \tag{ 16 }$$

Next, we rewrite equation (16) by assuming u₁ (t) = Y as below

$$\begin{align} \ddot{\text{Y}} & + c\dot {\text{Y}} + w_{{\text{non}}}^2Y \nonumber \\ & = {a_1}\int_0^1 {{{({\phi _1}(x)^{\prime} )}^2}} {\text{d}}x \int_0^1 {{\phi _1}} (x){\phi _1}(x)^{{\prime} {\prime}}\,\, {\text{d}}x\, {Y^3} \nonumber \\ & \quad + {a_2}[{({V_{{\text{D}}1}} - {V_{{\text{Beam}}}})^2}{F_1}(t) - {({V_{{\text{D}}0}} - {V_{{\text{Beam}}}})^2}{F_2}(t)] \nonumber \\ & \quad + {a_2}[{({V_{{\text{Beam}}}} - {V_{{\text{AC}}}}\cos (\Omega t))^2}{F_3}(t) - {V_{{\text{Beam}}}}^2{F_4}(t)] \nonumber \\ & \quad + {a_2}[{({V_{{\text{D}}3}} - {V_{{\text{Beam}}\,}})^2}{F_5}(t) - {({V_{{\text{D}}2}} - {V_{{\text{Beam}}}})^2}{F_6}(t)] \\[-12pt] \nonumber \end{align} \tag{ 17 }$$

in which Y is a modal coordinate amplitude. Above, the various function F_i (Y) are defined as

$$\begin{align} {F_1}(Y)& = \int_0^1 {\frac{{{\phi _1}(x)}}{{{{[\frac{{{d_1}}}{{{d_3}}} - {\phi _1}(x)Y]}^2}}}} [U({a_5} - x) - U({a_4} - x)]{\text{d}}x;{\text{ }}{F_2}(Y) \nonumber \\ & = \int_0^1 {\frac{{{\phi _1}(x)}}{{[\frac{{{d_0}}}{{{d_3}}} + {\phi _1}(x)Y){]^2}}}} [U({a_5} - x) - U({a_4} - x)]{\text{d}}x \nonumber \\ {F_3}(Y) & = \int_0^1 {\frac{{{\phi _1}(x)}}{{{{[\frac{d}{{{d_3}}} - {\phi _1}(x)Y]}^2}}}} [U({a_3} - x) - U({a_2} - x)]{\text{d}}x;{\text{ }}{F_4}(Y) \nonumber \\ & = \int_0^1 {\frac{{{\phi _1}(x)}}{{{{[\frac{d}{{{d_3}}} + {\phi _1}(x)Y]}^2}}}} ][U({a_3} - x) - U({a_2} - x)]{\text{d}}x \nonumber \\ {F_5}(Y) & = \int_0^1 {\frac{{{\phi _1}(x)}}{{{{[1 - {\phi _1}(x)Y]}^2}}}} [U({a_1} - x) - U({a_0} - x)]{\text{d}}x;{\text{ }}{F_6}(Y) \nonumber \\ & = \!\!\int_0^1 \!\!{\frac{{{\phi _1}(x)}}{{{{[\frac{{{d_2}}}{{{d_3}}}\! + \!{\phi _1}(x)Y]}^2}}}} ][U({a_1} - x) \!-\! U({a_0} - x)]{\text{d}}x. \end{align} \tag{ 18 }$$

Due to the term with Y in the denominator of the electrostatic force terms, the integral terms are not easily evaluated, which will cause computational problems [1, 39]. One approach to resolve this issue is by pre-calculating numerically these complicated terms F₁(Y), F₂(Y), F₃(Y), F₄(Y), F₅(Y), F₆(Y) of equation (18) for each step of ΔY = 0.05, and then tabulate and store them, as shown in the appendix. Once the table is built, it can be recalled for each time step of the needed numerical integration. Next, we can adopt the shooting technique to solve the forced vibration problem based on equation (17) combined with equation (18). It is noted that the shooting method [1] is a powerful numerical technique that is widely used to capture the periodic motions of the system, including the stable and unstable branches of the solutions. The stability of the periodic solution also can be assessed using the Floquet theory.

Here, to clarify the DAC function well in the linear domain, we solve the forced vibration response by the method of shooting and conduct the frequency sweeps around the primary resonance by using the optimized air gap sizes. One can notice that the frequency for each digital case in figure 5 shows linear response and shift to the higher frequency due to the different digital input ranging from the case-000 to case-111. When the DAC resonator is operated at any single driven frequency signal such as f₁ or f₂, the corresponding output response can be clearly distinguished. For example, the output level, at different driving signals f₁ or f₂, can significantly decrease or increase when subjected to the different digital input bits. Therefore, the optimized resonator-based DAC can function well in the linear domain compared to the previously proposed one of figure 3(a).

Zoom In Zoom Out Reset image size

It is worthy to mention that the proposed resonator-based DAC also can work in the nonlinear regime. Next, we can examine the operation of the proposed DAC in the nonlinear regime if we properly adjust the input power source and vacuum level (damping). Here, two experimental samples with very low damping are considered. It is noted that the frequency response curves are evaluated by calculating the average amplitude response value between maximum and minimum of the steady-state response for each case when conducting time integration.

Here, we compare the experimental data with the analytical results [34]. Figure 6(a) experimentally demonstrates the forward and backward frequency sweep of the resonator in the nonlinear regime for all digital combinations under −2 dBm AC input power level. Due to the different air gap sizes between the beam and digital input electrodes, distributions of the electrostatic load will be non-uniform along the micro beam and it will generate differences for each of the input combinations. In the nonlinear operation, one can notice that the amplitude response in the forward frequency sweep bends to higher values, nonlinear hardening behavior, due to the dominant mid-stretching nonlinearity over the softening electrostatic one. Figure 6(b) shows the S₂₁ parameter subjected to the digital inputs ranging from case-111 to case-000, in descending order. We can see that the digital case-000 input has the largest amplitude at operating point f₁ while case-111 has the smallest one. That is because the digital bits for case-111 generates the least influence on the actuation of the micro beam.

Zoom In Zoom Out Reset image size

The frequencies are swept under different digital input combinations and the inset for forward and backward frequency is shown in figure 7. The analytical results of figure 7 have good agreement with the experiment data of figure 6(a). To better understand the behavior at reached maximum excitation frequency point, we examine their eigenvalues (Floquet Multipliers) at that point. We found that the Multipliers of the system gradually approach the boundary of the unit circle when the excitation increases, and then the microsystem generates cyclic fold bifurcation phenomenon [1] at the maximum reached frequency.

Zoom In Zoom Out Reset image size

As mentioned, the shooting technique can capture all the periodic motions. Figure 8 depicts the appearance of three periodic solutions (for digital case-000) with different amplitude levels consisting of the upper branch (resonance stable branch, blue), the lower branch (non-resonance stable branch, magenta), and the intermediate branch (unstable branch, orange), which mainly depend on the initial conditions. As the frequency of excitation is increased or decreased, there are jumps in the response amplitude. For the forward excitation frequency sweep, we note that one Floquet multiplier approaches unity at the maximum amplitude of 3.59 µm, and also the upper stable branch and unstable branch are intersected at the point of the reached maximum frequency. This indicates a cyclic fold bifurcation [1].

Zoom In Zoom Out Reset image size

Figure 9 compares the dynamic response between theory and experiment at −12 dBm AC input signal power [30]. We notice that decreasing the AC source can lower the amplitude response. The frequency response curves experimentally and analytically show the same nonlinear hardening behavior. Thus, we observe a satisfactory agreement for the response among experimental and analytical results.

Zoom In Zoom Out Reset image size

Next, we notice that electrode D3 (with the smallest air gap), and V_Beam both play a significant role in the operation of the DAC. To clarify the influence of the electrode D3, we select the representative samples digital case-011 (represent the first four cases, D3 'off') and case-111 (represent the last four cases, D3 'on') in figures 10(a) and (b), respectively. From the simulated results, one can notice that lowering V_Beam can make the frequency-swept curve right-shifted and clearly distinguished due to the nonlinear softening effects. Also we note that V_Beam also can generate a big difference for case-011, compared to the case-111.

Zoom In Zoom Out Reset image size

In this work, we analytically and experimentally presented modeling, investigation, validation, and optimization of the MEMS resonator-based 3-bit DAC consisting of an in-plane clamped-clamped beam actuated by partial electrodes with different air-gap. First, we derived a dynamic equation of the proposed 3-bit DAC device, and the static response, linear EVPs, and dynamic response were solved then based on the method of the Galerkin-shooting technique. Second, to ensure the DAC function well in the linear domain, we optimized the size of air gaps. The optimized results indicate that the frequency responses ranging from the digital case '000' to '111' are monotonically increasing, linearly ordered, and significantly separated for the output at a single frequency. Third, we experimentally and analytically explored two samples of the DAC device under the nonlinear operation at −2 dBm and −12 dBm AC power levels, respectively. Forward and backward frequency sweeps were also investigated for all digital cases and were validated with measured results. The frequency response exhibit the characteristic bending toward higher frequencies and show the nonlinear hardening behaviors due to the high mid-plane stretching term. Moreover, we also notice that decreasing the AC power level from −2 dBm to −12 dBm can reduce the amplitude response and swept-range of the frequency curves. Last, we explored and discussed the effects of V_Beam on the DAC operation. Thus, the results suggest that the proposed modeling, simulations, and optimization analysis could be successfully implemented in the design of the DAC under various digital combinations. The rich nonlinear behavior with lower energy consumption could provide some high potential applications in IoT, such as logic, computation, sensing, and actuation.

We acknowledge the financial support from King Abdullah University of Science and Technology (KAUST).

The data that support the findings of this study are available upon reasonable request from the authors.

The direct integration of equation (16) can make the simulation very slow when we analyze the dynamic system [40]. To overcome this and improve the efficiency, first, we numerically integrate the complicated terms F_i (Y) and then store the values in a table choosing the step ΔY = 0.05. Here, the numerical evaluation F_i (Y) is reported.

Zoom In Zoom Out Reset image size