Bifurcation Analysis of a Micro-Machined Gyroscope with Nonlinear Stiffness and Electrostatic Forces

The bifurcation of the periodic response of a micro-machined gyroscope with cubic supporting stiffness and fractional electrostatic forces is investigated. The pull-in phenomenon is analyzed to show that the system can have a stable periodic response when the detecting voltage is kept within a certain range. The method of averaging and the residue theorem are employed to give the averaging equations for the case of primary resonance and 1:1 internal resonance. Transition sets on the driving/detecting voltage plane that divide the parameter plane into 12 persistent regions and the corresponding bifurcation diagrams are obtained via the singularity theory. The results show that multiple solutions of the resonance curves appear with a large driving voltage and a small detecting voltage, which may lead to an uncertain output of the gyroscope. The effects of driving and detecting voltages on mechanical sensitivity and nonlinearity are analyzed for three persistent regions considering the operation requirements of the micro-machined gyroscope. The results indicate that in the region with a small driving voltage, the mechanical sensitivity is much smaller. In the other two regions, the variations in the mechanical sensitivity and nonlinearity are analogous. It is possible that the system has a maximum mechanical sensitivity and minimum nonlinearity for an appropriate range of detecting voltages.


Introduction
Micro-machined gyroscopes are widely used because of their small size, light weight, simple structure, low manufacturing cost, and high reliability [1]. The performance of micro-machined gyroscopes is closely related to their dynamical characteristics. The existence of nonlinear factors makes their dynamical response more complex, so the nonlinear dynamics of micro-machined gyroscopes have attracted the attention of researchers over the last decades.
Braghin et al. [2] investigated a lumped parameter model of a MEMS gyroscope by using semi-analytical and numerical methods. The existence of hardening nonlinear stiffness was confirmed through experiments. Martynenko et al. [3] studied a MEMS gyroscope with a ring resonator supported by a flexible torsion system by considering the material nonlinearity. The results showed that the nonlinear elastic properties led to additional errors, unstable branches of resonance curves, and quenching. Ref. [4] also analyzed the nonlinear vibrations of an MEMS tuning fork gyroscope with an elastic rod. The influence of the slowly changing angular velocity was taken into consideration. Li et al. [5,6] considered the nonlinear dynamics of a machined gyroscope with nonlinear stiffness and delayed feedback control and analyzed the effects of system parameters and feedback gains. Mojahedi et al. [7] studied the static and dynamical characteristics of micro-/nano-mechanical gyroscopes with nonlinear electrostatic forces and intermolecular forces. The effects of parameters on static/dynamic deflection, natural frequency, and instability were discussed numerically. Nitzan et al. [8] discussed the influence of cubic nonlinearity on the error of rate-integrated gyroscopes. Hamed et al. [9] considered an MEMS gyroscope with linear and nonlinear parametric excitations. Considering the parametric resonance and 1:1 internal resonance, the steady-state response and its stability were discussed. With the existence of parameter uncertainties and noise measurement, comparisons between optimal linear feedback control and active control via negative nonlinear cubic velocity feedback were given. Shang et al. [10,11] considered a machined gyroscope with parametric excitation and cubic nonlinearity. The Hopf bifurcation condition was obtained, and a numerical simulation was used to study the driving and detecting amplitudes. Awrejcewicz et al. [12] studied a micromechanical gyroscope with torsion joints, considering the effects of stiffness and joint nonlinearity, and performed a multi-scale analysis of the steady-state response of the system in the case of primary and internal resonance.
Electrostatic force is also an important nonlinearity in these systems. Micro-machined gyroscopes are usually driven and detected by capacitors. The nonlinear effect of the electrostatic force of the driving and detecting capacitors cannot be negligible. Considering the nonlinearity of electrostatic forces, Lestev et al. [13] studied the steady-state response of a micro-machined gyroscope. Tsai et al. [14] investigated a micro-machined gyroscope with nonlinear stiffness and electrostatic forces, and showed the unstable region on the driving/detecting frequency plane. Then, the periodic solution and chaos were analyzed. Kacem et al. [15] established a continuum model of a resonant gyroscope. The steady-state response of the system was studied based on the Galerkin method and the perturbation method. Lajimi et al. [16] studied a micro-machined gyroscope consisting of a beam and a rigid body, and solved the equations with a multi-scale method. The response was analyzed for both amplitude modulation and frequency modulation. Ref. [17] also investigated the mechanical-thermal noise of a micro-machined gyroscope. Sharma et al. [18] numerically and experimentally analyzed the pull-in phenomenon of a micro-machined gyroscope and discussed the effects of the dynamic pull-in voltage and the measured angular velocity. Tehrani et al. [19] experimentally described the performance of an MEMS gyroscope with different forms of mechanical and electrostatic nonlinearities, and found that a low-angle random walk can be achieved even though the gyroscope's drive mode exhibits high amplitude-frequency dependence, and that the bias instability is largely independent of the operating regime. Ouakad [20] investigated the influences of the size effects on the pull-in phenomena and the primary resonance response of an MEMS gyroscope model given by the modified coupled-stress non-classical mechanics theory.
A micro-machined gyroscope system is a typical nonlinear vibration system with 1:1 internal resonance [8,10,12]. This 1:1 internal resonance is one of the important research directions of nonlinear dynamics [21,22]. It occurs in numerous mechanical systems, such as cables [23], beams [24,25], circular or square plates [26,27], etc. Many investigations, including those of periodic response [23], backbone curves [21,22], bifurcation characteristics [28], and chaos [29], have been done on this subject. The singularity theory is an effective method for bifurcation investigation. Based on this theory, the transition sets on the unfolding parameter plane and the corresponding bifurcation diagrams with different topological structures can be given [30]. Chen and Langford [31] studied a generally non-linear Mathieu equation and obtained six typical bifurcating response curves by setting the system parameters as the unfolding parameters. By using the singularity theory, Qin et al. [23] studied the bifurcation of an elastic cable with 1:1 internal resonance. Hou et al. [32] investigated the bifurcation modes of a periodic solution in a duffing system under constant force and harmonic excitation. Li et al. [33] considered the bifurcation of a coupled pitch-roll ship model with quadratic and cubic nonlinearity.
Generally speaking, the electrostatic force of a micro-machined gyroscope in the detecting direction is a fraction of the displacement; it has conventionally been developed into a Taylor series of the displacement in recent studies. However, Taylor expansion is accurate only in the neighborhood of the expansion point. For micro-machined gyroscopes, the amplitude of the detecting direction tends to be close to the gap of the capacitance to improve its sensitivity. In this case, the conclusions must be inaccurate with Taylor expansion. Moreover, most of the existing studies have discussed the influence of parameters on the response. Few studies have focused on the singularity analysis of micro-machined gyroscope systems.
In this paper, we propose an investigation of the bifurcation of a micro-machined gyroscope with nonlinear stiffness and electrostatic forces. In Section 2, the Lagrange equation is used to establish the dynamical equations of the micro-machined gyroscope. The pull-in phenomenon in the detecting direction is analyzed in Section 3. In Section 4, the method of averaging is employed to obtain the analytical solutions without Taylor expansion. Then, the influences of the driving and detecting voltages on the bifurcation characteristics of the amplitude-frequency response are discussed with the singularity theory in Section 5. Finally, the conclusion summarizes the work. Our innovation is the singularity analysis of a micro-machined gyroscope system with consideration of the fractional electrostatic force without Taylor expansion. Figure 1 gives a model of a micro-machined gyroscope. Figure 1a shows the schematic representation of the micro-gyroscope, and Figure 1b shows the equivalent spring-mass model. In Figure 1b, m is the proof mass of the sensitive structure. x a −o−y a is the inertial coordinate system, and x−o−y is the driving/detecting coordinate system of the micro-machined gyroscope. The origin o is located at the center of the proof mass, as the system is in static balance, as shown in Figure 1b. x is the driving direction, y is the detecting direction, and z is the direction of vertical paper facing outwards. θ is the angular displacement of the gyroscope around the z-axis, and ω z is the angular velocity being measured. As mentioned in Ref. [34], the angular velocity ω z is assumed to be constant. Therefore, Euler's second law of motion can be omitted in the considerations. Based on the relations between the displacements of the driving/detecting coordinate system and the inertial coordinate system shown in Figure 1c, the displacements of the proof mass in the inertial coordinate system can be given as x a = x cos ω z t − y sin ω z t, y a = x sin ω z t + y cos ω z t.

Physical Model of the Micro-Machined Gyroscope and Its Mathematical Description
(1) The velocity of the proof mass can be obtained from the derivative of Equation (1) with respect to time asẋ Then, the kinetic energy of the proof mass is The supporting beams of the micro-machined gyroscope have cubic nonlinear stiffness, and the elastic potential energy can be obtained as where k x , k y , α x , α y are the linear and nonlinear stiffness coefficients in the driving and detecting directions, respectively. The Rayleigh dissipation function of the system is where c x , c y are the damping coefficients in the driving and detecting directions, respectively. According to the Lagrange equation, the motion equations of the system are given as where −2mω zẏ and 2mω zẋ are the Coriolis forces, which are also known as the gyro forces of the gyroscope. F x and F y denote the electrostatic forces of the driving and detecting directions, respectively. The driving and detecting structures are shown in Figure 1d,e, respectively. When the proof mass has a displacement x, the two capacitances in the driving direction are where n is the number of combs, l is the overlap length of the driving combs, a is the thickness of the driving comb, and s d is the clearance between the combs. η is the dielectric constant. The driving voltages are Then, the electrostatic force of the driving direction can be given as Similarly, when the proof mass has a displacement y, the two capacitances in the detecting direction are where l l and l w are the length and thickness of the detecting electrode. s s is the capacitance gap. Then, the electrostatic force of the detecting direction can be given as In the stiffness term of Equation (6), mω z 2 can be ignored, as k x mω z 2 , k y mω z 2 . By setting X = x/s s , Y = y/s s , τ = ω 0 t, the motion equations can be rewritten in the non-dimensional form as where 1 is a small parameter, and e X , e Y are the non-dimensional parameters relevant to the driving and detecting voltages, respectively. Moreover, since s s is a certain value, Equation (12) has a constraint Y ≤ 1.

Pull-In Analysis of the Detecting Equation
It is known that the pull-in phenomenon, which means that stable equilibriums of the system cease to exist, takes place when the detecting voltage is greater than a critical value. The proof mass is contacted with the detecting electrode, and the system cannot work. Therefore, we focus on pull-in analysis of the motion equation of the detecting direction in this section. By ignoring items with ε and setting Y = Y 1 , Y = Y 2 , the motion equation in the detecting direction can be written in the form of state equations as Obviously, Equation (13) is a Hamiltonian system, and the equilibriums of the equation are determined by Equation (14). Figure 2 shows the bifurcation diagram of the equilibriums. The fork bifurcation occurs with the change in e Y . e Y =1 is the bifurcation point. As e Y > 1, the system has an unstable equilibrium (0,0), and as e Y < 1, the system has three equilibriums P 0 (0, 0), Here, P 0 is the center, and the other two equilibriums P 1 , P 2 are saddle points. In order to deeply understand the dynamics of Equation (13), the potential energy function and the Hamiltonian function can be obtained as  Figure 3 shows the phase trajectories and potential energy curves of the system. As e Y ≥ 1, the response of the undisturbed system will tend to infinity no matter what the initial value is. For the case of e Y < 1, the phase trajectory of the undisturbed system is a closed curve, which means that the system has a periodic response if the initial value is between two heteroclinic trajectories in the phase plane. When the initial value is given outside the two heteroclinic trajectories, the response of the system will also tend to infinity. The area between the two heteroclinic trajectories is the security region of the undisturbed system. As e Y ≥ 1, the potential energy curve has only a maximum point, and it is obviously impossible to have stable motion. As e Y < 1, the potential energy curve has a potential that is close to Y 1 = 0. If the amplitude of the system response is less than Y s , the motion is periodic.
Since the system contains a constraint Y 1 ≤ 1, if the response of the system is unstable, finally, Y 1 = 1, the capacitance gap in the detecting direction is zero, and the proof mass and the detecting electrode are closely attached, which is commonly known as pull-in. Figure 4 shows the influences of the detecting voltage on the security region of the undisturbed system. The area between the two heteroclinic trajectories is the security region. It can be seen that the security region of the undisturbed system decreases with the detecting voltage.

Approximate Analytical Solution of the Micro-Machined Gyroscope System
The performance of micro-machined gyroscope is closely related to its periodic response. In this section, the solutions are analytically solved with the method of averaging. (12) can be rewritten as The solutions of Equation (18) are assumed as [35] where A 1 and A 2 are the amplitudes, and θ 1 and θ 2 are phase angles in the x and y directions, respectively, which can be determined by the following averaging equations: π −π f 1 sin φ 1 dφ 1 , π −π f 1 cos φ 1 dφ 1 , where where σ 1 and σ 2 are the two detuning parameters that describe the nearness of primary resonance and 1:1 internal resonance: Note that f 1 and f 21 are polynomials, and f 22 is a fraction. The integration of Equation (20) is easy to solve, except for the fractional part. Thus, we set The integrand off A is It is obvious thatf A = 0, since Equation (22) is an odd function about φ 2 and the integral interval [−π, π] is symmetric about zero. On the calculation of f A , by setting z = exp(iφ 2 ), one has According to the residue theorem [36], substituting Equation (23) into f A gives where where z k are the isolated singularities of f (z) within the unit circle. Due to A 2 ≤ 1, Equation (25) has two secondary poles in the unit circle, which are Substituting Equation (26) into Equation (24) leads to Then, the averaging equations are obtained as By setting A 1 = A 2 = θ 1 = θ 2 = 0 and eliminating the trigonometric function terms in Equation (28), the bifurcation equations for the amplitude A 1 , A 2 can be given as To study the stability of the periodic response, we set A 1 = A 10 + p 1 , θ 1 = θ 10 + p 2 , A 2 = A 20 + p 3 , θ 2 = θ 20 + p 4 , where A 10 , A 20 , θ 10 , θ 20 are the steady-state solutions and p i (i = 1, 2, 3, 4) are small perturbations. Substituting them into Equation (28), and ignoring the high-order terms of p i , one has where p = [p 1 , p 2 , p 3 , p 4 ] T , and M is the Jacobian matrix of the periodic response. The expressions of its elements are given in the Appendix A. The stability of the steady-state solutions can be determined by the eigenvalues of M. The steady solution is stable if and only if the real part of all eigenvalues is negative. A positive real part of the eigenvalue leads to unstable solutions. The steady solution loses its stability as one real eigenvalue changes its sign, and the saddle-node bifurcation occurs. A pair of complex conjugate eigenvalues passing the imaginary axis of the complex plane at the bifurcation point leads to Hopf bifurcation [37]. If there are no special instructions, the values of the parameters in this paper are chosen as ξ X = ξ Y = 15, γ X = γ Y = 50, Ω z = 10, σ 2 = 50, ε = 0.001. Figure 5 shows the comparisons between the analytical solutions and the numerical solutions. The numerical solutions are given by the fourth-order Runge-Kutta method. It can be observed that the analytical solutions are in good agreement with the numerical solutions. Based on Equations (29) and (30) and the stability analysis, Figure 6 provides the results based on both the residue theorem and Taylor expansion. It is shown that the results obtained by the residue theorem are more accurate than the results based on Taylor expansion.
The resonance curves with different values of the driving and detecting voltages are shown in Figure 7. The figure shows that the topological structures of the curves are greatly different when E X , E Y are varied. Therefore, it is necessary to analyze the bifurcation characteristics to illustrate the effects of the driving and detecting voltages on the resonance curves.

Bifurcation Analyses
It can be obtained from the Equation (30) that Substituting it into Equation (29) gives The above equation is an algebraic equation about A 2 . Its solutions may bifurcate when parameters such as σ 1 and σ 2 are varied, since it is nonlinear. In this section, we have conducted a bifurcation analysis of Equation (32). For simplicity, σ 1 in Equation (32) is regarded as the bifurcation parameter, and σ 2 is fixed in the following discussion.
Then, the singularity theory [30,38] is employed to analyze the effects of E X and E Y on the bifurcation characteristics of Equation (32). The transition sets of Equation (32) are defined as D = B ∪ H ∪ DL, where B, H, DL denote the bifurcation sets, the hysteresis sets, and the double limit point sets, respectively. Their expressions are given in Table 1. Table 1. Expressions of the transition sets of the singularity theory.

Transition Sets Expressions
Bifurcation sets B , σ 1 ) = 0, i = 1, 2 and A 2(1) = A 2(2) Figure 8 shows the transition sets on the driving-detecting voltage plane, which divide the plane into 12 persistent regions. Figure 9 presents the bifurcation diagrams corresponding to the persistent regions in Figure 8. The corresponding values of E X and E Y are given in Table 2. It can be seen that, corresponding to region A, there is only one solution of A 2 for different σ 1 . However, in region B, corresponding to one interval of σ 1 , the system may have three solutions. In the bifurcation diagram of region C, multiple solutions appear in two intervals of σ 1 . In region D, a different solution branch appears, and in region E, folding happens in the added solution branch.  Table 2. Values of E X and E Y corresponding to subplots A through L of Figure 9. In order to explain the differences in the bifurcation diagrams corresponding to regions F through K, we defined four key points as points a through d. They are all turning points of the bifurcation curves, as shown in Figure 9. Figure 10 shows the partial enlargements of the bifurcation curves. We can see in the bifurcation diagram of region F that points a and b do not exist. However, corresponding to regions G and I, point d does not exist. The differences in G and I are that point c in the curves of region G is on the left side of point a, and point c in the curves of region I is between a and b. Corresponding to regions H, J, and K, points a through d all exist, while in the curves of J, point c is between a and b, and the curves of regions H and K are to the left of point a. In region K, point d is between a and b, and point d of regions J and H is on the right of point b. The different positions of these turning points may cause different jumping phenomena as the driving frequency is increased or decreased. Figure 9 also gives the comparisons between the numerical solutions and the analytical solutions. It can be seen that the numerical solutions are in good agreement with the analytical solutions, except for the case in which the value of σ 1 is large. This is because in the calculation of Equation (32), we set Ω = 1, but actually, Ω = 1 + εσ 1 .

Subplots
Micro-machined gyroscopes are not suitable for operation in multi-solution areas because the response may jump to other solutions, which may result in the gyroscope producing an incorrect output, as the system is disturbed. Figure 9 shows the multi-solution regions of the bifurcation curves, which are marked by the shaded areas. It can be seen in subplots D through L that most of the peak parts of the resonance curves are near multi-solution areas; thus, they do not meet the working requirements of the gyroscope. Only the bifurcation curves of regions A through C have large single-solution areas near the peak parts. Thus, in the following, we focus on the analysis of the response characteristics in regions A through C.
The effects of the driving and detecting voltages on the response are analyzed in regions A through C. The results are presented in Figure 11. In region A, the peak amplitude of the curves increases with the driving and the detecting voltages. In the bifurcation diagrams of region B, two saddle-node bifurcation points appear on the right side of the curve, which indicates that jump phenomena occur with the variation in σ 1 . The peak amplitude increases with the driving voltage and decreases with the detecting voltage. The effects of the driving and detecting voltages are opposite in regions A and B. Hence, the bifurcation analysis in this section is necessary. The peak frequency decreases with the driving and detecting voltages in region B. In region C, the saddle-node bifurcation points appear on both the left and right sides of the resonance curves, which means that jumping phenomena can take place two times-when the driving frequency is increased or decreased. The peak amplitude increases with the detecting voltage and decreases with the driving voltage.  The variations in the detecting amplitude with the measured angular velocity are the most important for micro-machined gyroscopes. Figure 12 shows the influences of the non-dimensional angular velocity Ω z on the detecting amplitude corresponding to different driving and detecting voltages. It can be observed in regions A through C, corresponding to any value of Ω z , that A 2 increases with the driving voltage. The effects of the detecting voltage are more complicated. In region A, and with a small value of Ω z in regions B and C, a medium detecting voltage leads to a large amplitude, while the curve of E Y = 20 changes sharply with a large value of Ω z in regions B and C. This change affects the linearity of the influence curve and is unfavorable for the operation of the gyroscope. Therefore, a small detecting voltage should be avoided in the B and C regions. In order to explain the influences of driving and detecting voltages on the response characteristics of the gyroscope clearly, mechanical sensitivity is considered. Mechanical sensitivity is the ratio of the detecting amplitude to the angular velocity [39,40]. However, for the nonlinear system, the ratio of the detecting amplitude to the angular velocity varies with the angular velocity. To calculate mechanical sensitivity, the least-square method is employed to fit the influence curves of the angular velocity linearly (shown in Figure 13). The fitting function is given as where S is defined as mechanical sensitivity. Meanwhile, in order to illustrate the difference between the fitting line and the original curve, we defined the expression of nonlinearity of the influence curves, which is given as follows.
13 Figure 13. Least-square fitting of the influence curve of the angular velocity; E X = 25, E Y = 30. Figures 14 and 15 show the variations in mechanical sensitivity and nonlinearity with driving and detecting voltages in regions A, B, and C. It can be seen that mechanical sensitivity and nonlinearity increase with the driving voltage on the whole. However, mechanical sensitivity is much smaller for the operation of the micro-gyroscope in region A compared with the results in regions B and C. Nonlinearity decreases monotonously with the increase in detecting voltage in region A. However, the variations of nonlinearity are non-monotonous in regions B and C. With the increase in detecting voltage, mechanical sensitivity increases until it peaks, and thenceforth, it decreases. The nonlinearity decreases initially, followed by an increase, and then again decreases as the detecting voltage increases.
Corresponding to E Y = 20 through 25 in regions B and C, the system may have the maximum mechanical sensitivity and minimum nonlinearity.

Conclusions
In this paper, the nonlinear dynamics of a micro-machined gyroscope system are presented with an approximate and numerical method that focuses on the effects of the driving and detecting voltages on the periodic motions and their bifurcations. Considering the nonlinear stiffness and the detecting electrostatic force, the dynamical equations of the micro-machined gyroscope are established with the Lagrange equation. The bifurcation of the equilibriums of the undisturbed system with the detecting voltage, as well as the pull-in effect, was studied, which showed that the system can have a stable periodic response only if the detecting voltage is kept within a certain range.
Through the method of averaging and the residue theorem, the dynamic equations with fractional terms were solved without Taylor expansion to find approximate analytical solutions. Numerical calculations verified the accuracy of the analytical solutions.
The influences of the driving and detecting voltages on the bifurcation characteristics of the system were studied through the singularity theory. Transition sets were given on the driving/detecting voltage plane, which was divided into 12 persistent regions. Analysis of the corresponding bifurcation diagrams showed that only three parameter regions met the operation requirements of the micro-machined gyroscope. In other regions that corresponded to a large driving voltage and small detecting voltage, most of the peak parts of the amplitude-frequency curves belonged to multi-solution areas.
The effects of the driving and detecting voltages on amplitude-frequency curves and the influence of these curves on angular velocity, as well as mechanical sensitivity and nonlinearity, were discussed for the three parameter regions. The results showed that in region A, which has a small driving voltage, mechanical sensitivity is much lower. In the other two regions, the variations in mechanical sensitivity and nonlinearity are analogous. It is possible that the system has a maximum mechanical sensitivity and minimum nonlinearity for an appropriate range of detecting voltages.

Data Availability Statement:
The data that support the findings of this study are available from the corresponding author, [H. Zhang], upon reasonable request.