Nonlinear Vibration Study Based on Uncertainty Analysis in MEMS Resonant Accelerometer

This paper aims to develop a resonant accelerometer for high-sensitivity detection and to investigate the nonlinear vibration of the MEMS resonant accelerometer driven by electrostatic comb fingers. First, a nonlinear vibration model of the resonator with comb fingers in a MEMS resonant accelerometer is established. Then, the nonlinear and nonlinear stiffness coefficients are calculated and analyzed with the Galérkin principle. The linear natural frequency, tracking error, and nonlinear frequency offset are obtained by multi-scale method. Finally, to further analyze the nonlinear vibration, a sample-based stochastic model is established, and the uncertainty analysis method is applied. It is concluded from the results that nonlinear vibration can be reduced by reducing the resonant beam length and increasing the resonant beam width and thickness. In addition, the resonant beam length and thickness have more significant effects, while the resonant beam width and the single concentrated mass of comb fingers have little effect, which are verified by experiments. The results of this research have proved that uncertainty analysis is an effective approach in nonlinear vibration analysis and instructional in practical resonant accelerometer design.


Introduction
The MEMS (micro-electro-mechanical system) resonant accelerometer can directly convert acceleration into frequency output. It benefits from having small size, light weight, low power, low cost, high resolution, good stability, wide dynamic range and quasi-digital output [1] etc. It is widely used in inertial navigation, seismic detection, intelligent robots and other fields, which has become an important development direction of MEMS sensors. However, when it is vibrating with large amplitude, hard spring nonlinearity of the resonator arises due to the small size of the resonator, which reduces the precision and even causes the sensor to be out of order. Therefore, the research on the nonlinear dynamics of the MEMS resonant accelerometer is of great significance to improve the precision and ensure its normal running.
In recent years, there have been various reports on the nonlinear research of the MEMS resonant accelerometer. In terms of temperature characteristics, Zhang proposed a nonlinear dynamic model of the resonant beam of the differential resonant accelerometer based on Hamilton's principle under varying temperature conditions, which could completely describe and analyze the nonlinear behavior of the resonant cavity [2]. Defoort et al. used the nonlinear amplitude frequency coupling effect to compensate for the resonator's passive temperature, which reduced silicon's temperature coefficient of frequency (TCF) to a level comparable with that of an AT-quartz resonator [3]. Then, Shin et al. also used the nonlinear amplitude frequency effect to improve the bias stability of resonant accelerometer in a large temperature range. In terms of mechanical coupling [4], Gusso proposed and theoretically studied the nonlinear damping mechanism of the transverse vibration of the double-clamped beam resonator [5]. Zou and Seshia optimized the bias voltage of the resonator by using the phase feedback oscillator circuit, and then improved the noise performance of the MEMS resonant accelerometer in the linear and nonlinear state in the wide band [6]. Juillard et al. studied the properties of nonlinearly operated weakly coupled resonators (WCRs) for resonant sensing applications nonlinear operation of the weak coupling resonator [7]. Lu et al. designed the nonlinear digital gain adjustment for rapid establishment of resonance oscillation and linearity improvement of MEMS vibratory gyroscopes [8].
In terms of process materials, Agarwal et al. established and verified the amplitude frequency dependence (A-F) effect model of MEMS resonators and studied the influence and mechanism of nonlinearity on the frequency stability [9]. Mahmoodi et al. undertook a comprehensive analysis and detailed comparative study on the nonlinearity of two types of microcantilever sensors actuated via a piezoelectric ZnO layer, which showed that the nonlinear relation between the stress and strain in some piezoelectric materials had a considerable effect on the sensor [10]. In 2019, Hashemi Kachapi et al. used Gurtin-Murdoch surface/interface theory to analyze the nonlinear vibration and frequency response of double wall piezoelectric nano resonators based on cylindrical nano shells [11]. Behbahani et al. used methods stemming from the ring dynamics [12] to tune the nonlinearity due to manufacturing imperfection in resonators [13,14], and similar frequency tuning method [15][16][17] can be applied to MEMS gyroscopes and accelerometers. These reports mainly used the methods of nonlinear amplitude frequency effect, phase oscillation feedback circuit, material performance optimization and modeling the manufacturing imperfections to reduce the influence of nonlinear vibration of resonators.
As for the uncertainty analysis method, it was first applied to sheetpile cofferdam design in 1987 [18], and then Padmanabhan and Pitchumani used stochastic model to investigate the influence of the uncertainty in the process and the material on the nonisothermal filling process [19]. After model improvement [20], the sample-based stochastic model was used to study the influence of uncertainty on the variability of refractive index, residual stress, maximum tension and defect concentration in the optical fiber stretching process [21][22][23]. At that time, stochastic model had been applied to safety assessments of technological systems [24], thermosetting-matrix composites fabrication [25], proton exchange membrane (PEM) fuel cells [26]. In 2012, Peng et al. developed the sample-based stochastic model to analyze the influence of the uncertainty of parameters on the solid-liquid-vapor phase change of metal particles and identified the laser fluence had dominant effects [27]. Since then, more applications of the stochastic model have been found in the design optimization of resonators [28], thermal damage of living biological tissues by laser irradiation [29], and fluctuation parameters on flow stability [30]. Therefore, the uncertainty method will be a mighty tool to study the nonlinear dynamics of the MEMS resonant accelerometer.
Shi and Fan et al. applied uncertainty method to investigate the effects of different uncertain parameters in electro-thermal excited MEMS resonant sensor, which demonstrated convincingly that the DC excitation voltage had dominant effects [31]. However, it is more complicated to analyze the nonlinear dynamics of the MEMS resonant accelerometer driven by comb.
In this paper, the nonlinear vibration model of resonant beam which is driven by comb fingers in a MEMS resonant accelerometer is established. After being deduced by the Galérkin principle and multi-scale method, the relation between the equivalent stiffness of the resonator and the geometric parameters of the accelerometer is discussed. In the other hand, the influence of pairs of comb fingers on the nonlinear effect of resonator is analyzed. Furthermore, the connection between the natural frequency-tracking error of the accelerometer caused by the nonlinear vibration and the measured acceleration and the geometric parameters of the resonant beam is investigated. Given the uncertainty distribution of geometric parameters and single concentrated mass of the comb fingers due to fabricating errors, the sample-based stochastic model will be applied to investigate influence on vibrating nonlinearity including linear natural frequency, nonlinear frequency offset and their ratio. Finally, we will design a circuit experiment to verify the validity of the uncertainty method.

Working Principle
The MEMS resonant accelerometer mainly includes the mass block, support beams, resonators, drive units and detection units. Each resonator consists of two identical resonant beams. The working principle of the MEMS resonant accelerometer is shown in Figure 1a. The resonator works in the resonant state, and the acceleration acts on the mass block. When the natural frequency of the resonator changes under the inertial force along the axial direction, the drive unit excites the resonator to vibrate and maintains the resonant state. Meanwhile, the detection unit detects the vibration signal, and the drive unit tracks the natural frequency of the resonator controlled by the closed-loop feedback control circuit. It guarantees the frequency of the excitation force is always consistent with the natural frequency of the resonator. Meanwhile, the detection unit detects the vibration signal, and the drive unit tracks the natural frequency of the resonator controlled by the closed-loop feedback control circuit, which guarantees the frequency of the excitation force is always consistent with the natural frequency of the resonator. According to the vibration signal detected by the detection unit, the change of natural frequency can be obtained, and then the measured acceleration can be converted. Figure 1b shows a schematic of the drive/detection units. When the driving voltage (AC and DC bias voltage) is applied to the fixed comb fingers, the alternating electric field will generate a lateral driving force. Under the action of the lateral driving force, the active comb fingers will generate reciprocating vibration relative to the fixed comb fingers. When the frequency of the driving voltage is consistent with the natural frequency of the active comb fingers, the active comb fingers will resonate, and the resonant frequency of the resonator will be obtained through the detection.   Figure 2 shows the equivalent mechanical model of the resonator single beam in the MEMS resonant accelerometer proposed in Figure 1. Neglecting the moment of inertia of the comb fingers, the comb finger is simplified as a particle attached to the double-clamped resonant beam. Taking the nonlinear geometric effect into account, the equation of motion of the resonant beam is obtained in line with Euler-Bernoulli beam modeling when the electrostatic force is considered. The lateral deflection w(x i ,t) is

Nonlinear Vibration Model of Resonator with Comb Fingers
The boundary conditions of the double-clamped resonant beam are q(t) is the electrostatic force generated by a single driving comb finger, which is defined as where ρ is the material density. L,B and H represent the length, width and thickness of the resonant beam, respectively. E is the young's modulus of the material. I is the moment of inertia. N a is the axial load caused by the measured acceleration. N r is the residual stress. m c is the single concentrated mass of the comb finger on the resonant beam. N is the pairs of the comb fingers on the resonant beam.x i is the distance from coordinate of the ith concentrated comb mass to the resonant beam end andt is time. c is the damping coefficient c is simplified to be frequency-independent in the neighborhood of resonance and the value is determined by the measured quality factor. δ is the Delta function, which is used to describe the position of the driving moments. ε 0 is the permittivity in vacuum. b e is the width of comb finger. g is the gap between two plates of comb capacitance. U p and U d are the dc voltage applied to the structure electrode.ω is the working frequency. Upon Selecting dimensionless parameters where r = √ I/BH is the radius of rotation of cross section of resonant beam, the moment of inertia is I = BH 3 /12, andω n is the linear natural frequency of the resonant beam. Equation (1) and the boundary conditions are changed as The Galérkin principle is a method to discretize the partial differential equation into a reduced-order model. Assuming the lateral deflection w(x, t) of beam is where φ i (x) is the ith mode of vibration of the resonant beam, and u i (x) is the generalized coordinate corresponding to the ith mode [32]. φ 1 (x) is first-order linear undamped mode function of the resonant beam. The solutions of Equations (5) and (6) can be expressed by mode of the resonant function. Through multiplying and integrating from 0 to 1, a reduced-order model in the form of ordinary differential equations can be obtained. Since the distribution of electrostatic force is symmetrical about the midpoint of the resonant beam and the frequency of electrostatic force is close to the first natural frequency, it can be considered that the resonant beam vibrates approximately according to the first mode. According to the undamped free vibration equation of the resonant beam, the first-order vibration mode satisfies the following relationship Replacing the first-order vibration mode function and Equation (8) into (5) can obtain Multiplying both sides of Equation (9) by φ 1 (x) and integrating both sides of Equation (9) with respect to x from 0 to 1, one obtains where the mechanical quality factor of the resonant beam first-order modal vibration Q, the amplitude of equivalent excitation force F eq , equivalent linear stiffness coefficient k 1 and nonlinear stiffness coefficient k 3 are expressed as Equation (12) is the ratio of nonlinear stiffness coefficient to linear stiffness coefficient, which can reflect the nonlinear degree of resonant beam vibration.
Supposing that when the accelerometer receives positive acceleration, the resonators receive tensile axial force, and when the accelerometers receive negative acceleration, the resonators receive compressive axial force. According to the determined structural parameters of the resonators with comb fingers, the ratio of nonlinear to linear stiffness coefficient k 3 /k 1 with the measured acceleration and the pairs of comb fingers N can be obtained from Equation (12), as shown in the Figure 3. Although increasing the pairs of comb fingers on the resonant beam increases the additional mass of the resonant beam, it can be seen from Figure 3 that increasing the pairs of comb fingers has small effect on the reduction of nonlinearity. The approximate analytical solution of Equation (10) can be obtained by using the multi-scale method. According to Figure 3 and Equation (12), the nonlinear stiffness coefficient k 3 1, which can be designated as the small parameter ε in the multi-scale method. Generally, the resonator in resonant accelerometer works in vacuum environment to obtain high mechanical quality factor Q, so 1/Q 1 can be expressed as The resonant beam vibrates approximately at their natural frequency when the accelerometer works, so let ω = 1 + εσ, where σ is the detuning parameter. When the resonant beam vibrates according to the natural frequency, a small excitation amplitude can cause a large vibration of the resonant beam, so let F eq = εK, where ω n can be obtained from Equation (12) and considering the relationship between k 3 and ε ω n = − The linear natural frequency f rn is So far, Equation (10) can be rewritten as Using the multi-scale method, the amplitude frequency response equation and the phase frequency characteristic equation of the resonant beam are obtained aŝ where γ is the phase-shift andα is the amplitude of resonant beam.

Natural Frequency-Tracking Error Caused by Nonlinear Vibration
When the phase-locked closed-loop circuit is used to track the natural frequency of the accelerometer, the phase-shift γ of the resonant beam is locked at π/2. Substituting γ = π/2 into Equation (20), the vibration frequency of the resonant beam is expressed aŝ The frequency-tracking error of the accelerometer can be obtained by Equation (21) The nonlinear frequency offset f o f f is Figure 4 presents the change curve of the natural frequency-tracking error of the accelerometer caused by the nonlinear vibration with the measured acceleration and the geometric parameters of the resonant beam. As seen from Figure 4, when the measured acceleration increases, the frequency-tracking error increases with the increase of resonant beam length L (Figure 4a) and decreases with the increase of resonant beam width B and thickness H (Figure 4b,c). Therefore, the nonlinear effect can be reduced by decreasing the length and increasing the width and thickness of the resonant beam.

Uncertainty Analysis Method
Based on the nonlinear vibration model of resonator, the nonlinear dynamics of MEMS resonant accelerometers is studied by using the sample-based stochastic model. Figure 5 illustrates the detailed procedure on how the sample-based stochastic modeling is realized based on the specific deterministic physical modeling for nonlinear vibration of resonator described with random selected stochastic instances. First, the input parameters are selected, and the degree of change is quantified. Then, through stochastic convergence analysis, a proper number of sample combinations of input parameters are determined, and the uncertainties of input parameters are calculated by the previously established deterministic physical model. Finally, the variability of output parameters is quantified according to the uncertainty of input parameters.
In this study, we focus on the uncertainty in four parameters: the length L, width B, thickness H and single concentrated mass of the comb finger m c of the resonant beam. The uncertainty of all the input parameters is assumed to obey the Gaussian distribution, which is commonly used to represent uncertain parameters of a physical model [20,26]. After determining the distribution of the input parameters, the Monte Carlo sampling method (MCS) is used to randomly select the input parameters from the given Gaussian distribution to obtain a sample combination. In the process of stochastic convergence analysis, with the increase of the number of samples, the mean value and standard deviation of input parameters will converge to the nominal mean value and standard deviation of Gaussian distribution, then the minimum sample value of input parameters will be obtained. To measure the uncertainty of input parameters, the coefficient of variance (COV) defined as σ/µ represents the uncertainty degree of the input parameters, where the average value (µ) is represented by the nominal value of the uncertainty parameters, and the standard deviation (σ) represents the variability of the input parameters.

Results and Siscussions
According to the finite element analysis method [28], we assume the nominal mean values of L, B, H and m c are 500 µm, 40 µm, 4 µm, 4.66 × 10 −12 kg and the other relevant properties of the resonant beam in this paper are shown in Table 1.
To obtain the number of input samples Ns, the COV of each input parameters is set to be 0.04. Figure 6 shows the results of the stochastic convergence analysis of the mean values of input parameters L, B, H and m c . It can be seen that the mean values fluctuate significantly with the decrease of sample numbers but converge as the numbers of samples increases. When N s = 300, the mean values of the input parameters converge to within 1% for the distributions. In other words, the number of samples of N s = 300 may be sufficient so far, but it also depends on the standard deviations of input parameters. (d) The stochastic convergence analysis of the standard deviations of four input parameters are presented in Figure 7, which indicates that there is still a large fluctuation at 300 samples. The standard deviation values converge to within 1.3% for all input parameters while the number of samples exceeding 400.  The results plotted in Figure 9 show the stochastic convergence analysis of the standard deviations of output parameters. The standard deviation converges to within 0.63% for f rn , 0.03% for f o f f and 1.40% for f nol that all output parameters are in 2% when 400 samples are used. In consequence, the minimum number of samples of N s = 400 is obtained and will be used in the following study.   It should be noted that when the COV of one of the input parameters increases from 0.01 to 0.1, the COVs of other input parameters remain at 0.01. It can be seen from Figure 10a that the IQR of the linear natural frequency f rn significantly increases from about 12.64 kHz to 123.046 kHz with the COV of the length L increases from 0.01 to 0.1. Moreover, when the COV of thickness H increases from 0.01 to 0.1, the IQR of f rn also shows an observable increase from 14.436 kHz to 64.068 kHz. However, the increase of COV of other parameters has very little effect on IQR. The results are consistent with the f rn being a strong function of the length L and the thickness H.
The IQR of the nonlinear natural frequency f o f f as a function of the COV of various input parameters is presented in Figure 10b, which shows that the IQR of f o f f increases significantly from 1.562 × 10 −9 Hz to 1.326 × 10 −8 Hz with the COV of the thickness H increasing from 0.01 to 0.1. Meanwhile, when the COV of length L increases from 0.01 to 0.1, the IQR of f o f f also shows an observable increase from 1.464 × 10 −9 Hz to 8.022 × 10 −9 Hz. It is relatively unaffected by the increase in the COV of other uncertain parameters. Furthermore, the IQR of the scaling factor f nol is a function of the COVs of the input parameters, as depicted in Figure 10c, where it is seen that the COV of the length L and thickness H increasing from 0.01 to 0.1 result in a remarkable increase of the IQR of f nol . The results presented the thickness H and the length L of the resonant beam are the main factors that cause the nonlinear vibration of the MEMS resonant accelerometer, and the influence degree of the thickness H is similar to the length L. This is consistent with the result that the natural frequency-tracking error of the accelerometer caused by the nonlinear vibration changes with the measured acceleration and the geometric parameters of the resonant beam. In addition, the width B of the resonant beam and single concentrated mass of the comb finger m c have little influence on the nonlinearity of the MEMS resonant accelerometer, which is consistent with the pairs of comb fingers making little difference to reduce nonlinearity (Figure 3).

The Equivalent Circuit Model
From the analysis of sample-based stochastic, it can be found that the length L and thickness H of the resonant beam have a great influence on the nonlinear vibration of the resonator. To verify this conclusion, the electromechanical equivalence principle forms a very effective base on which each term of the dynamic Equation (18) is converted into a one-to-one corresponding circuit unit and the equivalent circuit is established. Then, the equivalent circuit of the MEMS resonant accelerometer is established based on the study of its nonlinear vibration theory. In the circuit design, the excitation force only maintains the vibration of the resonant beam and has no influence on the vibration frequency, so the excitation term of Equation (18) is not considered.
The equivalent circuit model of the nonlinear vibration model of resonator with comb fingers is shown in Figure 11. It gives the circuit units corresponding to the terms in the differential Equation (18). u 1 (t) is the expression of first-order modal vibration of the resonant beam, which can be obtained by twice integrations ofü 1 (t). The −u 1 (t) term can be obtained by using proportional amplifier x as a feedback loop. The term of nonlinear stiffness coefficient k 3 u 3 1 (t) can be obtained by using proportional amplifier z and multiplier. The −u 1 (t) term corresponds to the proportional amplifier y. The three terms mentioned above are added by an adder to be theü 1 (t) term, then theu 1 (t) term and the u 1 (t) term can be obtained by integrator { and integrator | respectively. The nonlinear frequency offset is obtained by detecting the difference between the output frequency of u 1 (t) and the linear natural frequency calculated by theory.

Multiplier
Adder The circuit principle diagram of the nonlinear vibration model of resonator with comb fingers is demonstrated in Figure 12. The green box corresponds to the Proportional amplifier x, the pink box to the Proportional amplifier y, the yellow box to the Adder, the purple box to the Integrator {, the orange box to the Integrator |, the blue box to the Proportional amplifier z, and the red box to the Multiplier. The nonlinear coefficients corresponding to different geometric parameters are obtained by adjusting the feedback resistance R2 of the proportional amplifier in the blue box. Figure 13 is the diagram of experimental setup. The green box represents the experimental verification module. First, we assume the nominal mean values of L, B, H and m c are 500 µm, 40 µm, 4 µm and 4.66 × 10 −12 kg. Twelve points are sampled for each group of geometric parameters, and to ensure the comparability of the four group of data, one-fortieth of the nominal value of each change is taken. It should be noted that when studying the influence of one geometric parameter, keep the other three geometric parameters as nominal values. As mentioned above, the corresponding nonlinear coefficient k 3 will also change when the other three parameters are kept unchanged and only one parameter is changed. The change in k 3 can be achieved by adjusting the variable resistor shown in Figure 12. Secondly, the signal generator gives a signal to the experimental verification device to make it working in a resonant state. Then, adjusting R2 according to the change of geometric parameters, and the measuring instrument records the different frequencies of the output signal. Finally, through the data processing unit analyzes and processes the output frequency recorded by the measuring instrument, the variation of the detected frequency with geometric parameters is obtained. At the same time, the nonlinear frequency offset varying with geometric parameters can be obtained by comparing the detection frequency with the theoretical frequency.    Figure 14 illustrates the theoretical frequency and detected frequency varying with geometric parameters and the trend fitting curves of the nonlinear frequency offset f o f f . It can be seen from Figure 14a that as the length L increases from 437.5 µm to 575 µm, both the theoretical linear natural frequency and the nonlinear natural frequency of the circuit output show a decrease trend between 300 kHz and 700 kHz. Moreover, when the length L increases from 437.5 µm to 575 µm, the nonlinear frequency offset f o f f shows an observable increase from −150 kHz to 60 kHz in the fitting graph.   Figure 14b presents the theoretical frequency and detected frequency varying with the thickness H. For the convenience of observing the trend, the absolute value of the nonlinear frequency offset point after the nominal value of the resonator thickness H is taken. It can be observed that the theoretical frequency and the detected frequency of the circuit output show an increase trend from 150 kHz to 700 kHz with the thickness H increasing from 3.6 µm to 4.7 µm. Furthermore, the nonlinear frequency offset f o f f shows a remarkable increase from −85.156 kHz to 152.802 kHz in the fitting graph, which indicates that the nonlinear frequency offset is also strongly affected by the thickness H.

Experimental Results
The results plotted in Figure 14c show the theoretical frequency and detected frequency varying with the width B. It can be seen that as the width B increases from 36 µm to 47 µm, both the theoretical frequency and the detected frequency of the circuit output show a slightly upward trend from 380 kHz to 510 kHz. Meanwhile, the nonlinear frequency offset f o f f of the width B increases more tardily from −32 kHz to 5 kHz in the fitting graph than the length L and the thickness H.
Moreover, when the single concentrated mass of the comb finger m c increases from 4.0775 ng to 5.359 ng, as depicted in Figure 14d, it can be seen that neither the theoretical frequency nor the detected frequency of the circuit output changes significantly. The nonlinear frequency offset f o f f of the single comb finger mass m c varies from −0.416 kHz to 0.7912 kHz, which indicates that the single comb finger mass m c has little effect on the nonlinear frequency offset f o f f .
In conclusion, the detected frequency of the circuit output and the theoretical frequency change significantly with the increase of the length L and the thickness H, followed by B and m c . Meanwhile, the length L and the thickness H are the main factor affecting the nonlinear frequency offset f o f f , which is consistent with the results of the uncertainty analysis method.

Conclusions
• The nonlinear vibration model of resonator with comb fingers has been established. The nonlinear stiffness coefficient k 3 and the linear stiffness coefficient k 1 have been calculated and analyzed with the Galérkin principle. The linear natural frequency f rn , the nonlinear frequency offset f o f f and the scaling factor f nol are obtained by multi-scale method. • According to theoretical analysis, it is found that the pairs of comb fingers have little effect on the nonlinear vibration of the resonator. After further analyzing the relationship between the natural frequency-tracking error and the geometric parameters of the resonant beam, we can find that the frequency-tracking error increases with the increase of the length L of the resonant beam, and decreases with the increase of the width B and thickness H of the resonant beam. The nonlinear effect can be reduced by reducing the length L and increasing the width B and thickness H of the resonant beam. In the experimental verification, the detected frequency of the circuit output and the theoretical frequency change significantly with the increase of the length L and the thickness H, followed by B and m c . Meanwhile, the length L and the thickness H are the main factors affecting the nonlinear frequency offset f o f f , which is consistent with the results of the uncertainty analysis method. According to the above conclusions, to reduce the nonlinear characteristics of the sensor, four geometric parameters can be properly adjusted to increase the stiffness and obtain better structural performance.