A mechanical–thermal noise analysis of a nonlinear microgyroscope

Mechanical Systems and Signal Processing 83 (2017) 163–175

Contents lists available at ScienceDirect

Mechanical Systems and Signal Processing journal homepage: www.elsevier.com/locate/ymssp

A mechanical–thermal noise analysis of a nonlinear microgyroscope S.A.M. Lajimi n, G.R. Heppler, E.M. Abdel-Rahman Department of Systems Design Engineering, University of Waterloo, Canada N2L 3G1

a r t i c l e i n f o

abstract

Article history: Received 14 April 2015 Received in revised form 31 May 2016 Accepted 6 June 2016 Available online 22 June 2016

The mechanical–thermal noise (MTN) equivalent rotation rate (Ωn) is computed by using the linear approximation of the system response and the nonlinear “slow” system. The slow system, which is obtained using the method of multiple scales, is used to identify the linear single-valued response of the system. The linear estimate of the noise equivalent rate fails as the drive direction stroke increases. It becomes imperative in these conditions to use a more complex nonlinear estimate of the noise equivalent rate developed here for the first time in literature. The proposed design achieves a high performance regarding noise equivalent rotation rate. & 2016 Elsevier Ltd. All rights reserved.

Keywords: Noise-equivalent rotation rate MEMS gyroscope Method of multiple scales

1. Introduction To hold a thermal equilibrium between a mechanical resonator and its surrounding environment, a fluctuating force should cause thermodynamic fluctuations in the oscillator's motion [1]. The mechanical–thermal (thermomechanical) noise is concerned with the molecular agitation where the (micro) structure and the surrounding fluid interact [2]. At the micro-scale, the mechanical–thermal noise becomes the dominant noise mechanism and requires special attention and careful characterization [2,3]. For a micro-/macro-mechanical resonator, the thermal noise is modeled by applying a root-mean-square (rms) force, Fn = 4 kB T c BW , to the structure [2–4]. Symbols are introduced in nomenclature. Considering a vibratory single-mass gyroscope, Annovazzi-Lodi and Merlo [3] and later Bao [5] computed the noiseequivalent rotation rate Ωn. Leland [6] studied the mechanical–thermal noise employing the method of averaging based on a the work of Lynch [7]. Leland considered a system of linear equations similar to that of Annovazzi-Lodi and Merlo [3], used the Laplace transform, and solved the resulting algebraic equations for the rotation rate in terms of noise force components and arrived at a relation to compute Ωn similar to the previous works. Currently, the microgyroscope industry requires the estimation of the thermomechanical noise effect and consequent noise equivalent rotation rate. For a beam-rigid-body MEMS gyroscope, we use a comprehensive mathematical model accurately modeling any beammass gyroscopes regardless of the size of the end-mass. Using the method of multiple-scales, we obtain the slow system, i.e. the modulation equations, which is then used for mechanical–thermal noise analysis. An approximate linear approach is used to obtain a closed-form relation for computing the Ωn of beam-mass-based gyroscopes to corroborate the results of the n

Corresponding author. E-mail addresses: [email protected], [email protected] (S.A.M. Lajimi).

http://dx.doi.org/10.1016/j.ymssp.2016.06.005 0888-3270/& 2016 Elsevier Ltd. All rights reserved.

164

S.A.M. Lajimi et al. / Mechanical Systems and Signal Processing 83 (2017) 163–175

Nomenclature kB c w (x , t ) L aB bM e

Vw gv g hw Iζζ E Aw M J M ηη

J B ξξ

J B ζζ

Qd

ωd

q^(t^ )

Δ

Boltzman constant (kB = 1.38 × 10−23 J K−1) damping coefficient in newton-seconds/meter flexural displacement in the z direction relative to the base frame length of beam width of beam thickness of the end rigid body eccentricity, the distance between the end of the cantilever beam and the center of mass of the rigid body DC voltage for the drive electrode initial distance between the rigid body and the sense electrode initial distance between the rigid body and each electrode ( g w = gv = g ) width of the drive electrode beam cross sectional second moment of area about ζ Young's modulus (E¼ 169 GPa) electrostatic forcing area for the drive electrode mass of the rigid body moment of inertia of the rigid body about the ηaxis beam's moment inertia of the beam about the ξ-axis beam's moment inertia of the beam about the ζ-axis quality factor in the drive direction natural frequency in the drive direction modal coordinate of the system response in the drive direction a beam cross sectional ratio, Δ = bB

electrostatic constant ν =

T BW

absolute temperature in Kelvin noise bandwidth defined by the mechanical resonator or the electronic filtering flexural displacement in the y direction relative to the base frame angular rotation rate about the x-axis thickness of beam length of the end rigid body width of the end rigid body DC voltage for the sense electrode permittivity coefficient (ϵ = 8.854 × 10−12 Fm−1) initial distance between the rigid body and the drive electrode width of the sense electrode width of each electrode (h w = hv = h) beam cross sectional second moment of area about η density (ρ ¼2330 kg/m3) electrostatic forcing area for the sense electrode moment of inertia of the rigid body about the ξaxis moment of inertia of the rigid body about the ζaxis beam's moment inertia of the beam about the η-axis 4 time constant κ = 12 ρ 2L

v (x , t )

Ω bB LM aM Vv

ϵ

gw hv h Iηη

ρ Av J M ξξ

J M ζζ

J B ηη

κ

E bB4 g 3

E bB

Qs

quality factor in the sense direction natural frequency in the sense direction

p^ (t^ )

modal coordinate of the system response in the sense direction b2 an aspect ratio constant α = B 2

ωs

α

B

6 ϵ h L4

ν

12L

nonlinear analysis for an appropriate set of parameters where the response is linear and single-valued. The proposed design of the micro-structure acquires a high performance by minimizing the Ωn.

2. Mathematical methods The schematic of the sensor is provided in Fig. 1a. A picture of one of our fabricated gyroscopes and a picture of the experimental setup are provided in Figs. 1b and c. Employing the extended Hamilton's principle, the governing equations of motion are obtained as [8–10]:

¨ + 2 m Ω v ̇ + m Ω̇ v − m Ω2w + c w ẇ − J EIηη w ⁗ + m w

Ω2w ″ − J

B ηη

EIζζ v⁗ + m v¨ − 2 m Ω ẇ − m Ω̇ w − m Ω2 v + cv v ̇ − J

B ζζ

Ω̇ v″ − J

B ηη

Ω2 v″ − J B ηη

¨ ″ = 0, w (1)

B ηη

Ω̇ w ″ − J B ζζ

v¨″ = 0, (2)

and the boundary conditions are given by

w = 0,

w ′ = 0,

v = 0,

v′ = 0 at x = 0, and

(3)

S.A.M. Lajimi et al. / Mechanical Systems and Signal Processing 83 (2017) 163–175

165

VDC Base L a Ω

Sense Electrode

e

2 ev

b

y

hv

aM

bM

z

x

VDC Sense Electrode VDC

LM Drive Electrode

VAC

hw

2 ew

Fig. 1. Schematic of the beam-rigid-body microgyroscope (MEMS Gyroscope). The figure is not drawn to scale. (a) An image of one of our fabricated gyroscopes (b) and a picture of the test setup (c).

¨ + ew ¨ ′)) + J (Ω̇ v′ + Ωv′̇ + w ¨ ′) + J (Ω2 w ′ − Ω v′̇ ) EIηη w ″ − M e(Ω2(w + e w ′) − Ω̇ (v + e v′) − 2 Ω(v ̇ + e v′̇ ) − (w M ηη

− J (Ω2 w ′ − Ω v′̇ ) = M ζζ

ϵ

ew A w Vw2

2(g w − w − e w ′)2

M ξξ

at x = L, (4)

EIζζ v″ − M e(Ω2(v + e v′) + Ω̇ (w + e w ′) + 2 Ω(ẇ + e ẇ ′) − (v¨ + ev¨′)) + J (v¨′ − Ω̇ w ′ − Ω ẇ ′) − J (Ω2 v′ + Ω ẇ ′) M ζζ

M ηη

ϵ e A (v(L, t ) + ev v′(L, t ))Vv2 at x = L, + J (Ω2 v′ + Ω ẇ ′ + Ω̇ w ′) = 2 v (gv − (v(L, t ) + e v′(L, t ))2)2 M ξξ

(5)

¨ + ew ¨ ′)) − J (Ω̇ v′ + Ω2w ′ + w ¨ ′) = − EIηη w ‴ + M (Ω2(w + ew ′) − Ω̇ (v + ev′) − 2Ω(v ̇ + e v′̇ )−(w B ηη

2(g w − w − e w ′)2

at x (6)

= L, EIζζ v‴ + M (Ω2(v + ev′) + Ω̇ (w + ew ′) + 2Ω(ẇ + e ẇ ′)−(v¨ + ev¨′)) − J (Ω2v′ + v¨′) − J B ζζ

=

ϵ

A w Vw2

Ω̇ w ′

B ηη

ϵ A v (v(L, t ) + ev v′(L, t ))Vv2 at x = L. (gv2 − (v(L, t ) + e v′(L, t ))2)2

(7)

The method of multiple scales is applied on a discretized model which is obtained by replacing the displacements with their single-mode approximations in the energy expressions and applying the Lagrange's equations to the resulting discretized Lagrangian [8–10]. The discretized model is given by:

¨ ̇ ^ (t^)p^ (̇ t^) + γ Ω ^ (̇ t^)p^ (t^) + (ΔΓ ″ − γ Ω ^ 2)q^(t^) = γ1,1 q^(t ) + c^q q^(t^) + γ1,2Ω 1,3 1,4

ν e^w γ1,5 Vw2 (1 − γ1,5 q^(t^))2

¨ ̇ ^ (t )q^(̇ t^) − γ Ω ^ (̇ t^)q^(t^) + (Δ3 Λ″ − γ Ω ^ 2)p^ (t^) = γ2,1 p^ (t^) + c^p p^ (t^) − γ2,2 Ω 2,3 2,4

2 2ν e^v γ2,5 Vv2 2 ^ ^ 2 2 (1 − γ2,5 p (t ) )

(8)

(9)

The explicit forms of the coefficients are given in Appendix A. We investigate the case of single sense electrode sensor. Separating the response into its static and dynamic components, q(t ) = qs + qd(t ) and p(t ) = ps + pd (t ), in Eqs. (8) and (9), and ^ for notational clarity). Substituting the static dropping the dynamic terms gives the static equations (we remove the () equations into the decomposed equations results in the dynamic equations in the form:

166

S.A.M. Lajimi et al. / Mechanical Systems and Signal Processing 83 (2017) 163–175

ωq2qd(t ) + q¨d(t ) = 2Q 14VAC Vwcos(tΩe ) 1 2 Q 14VAC cos(2tΩe ) − qd(t )2(Q 8Ωpḋ (t ) + Q 5qḋ (t ) + Q 2q¨d(t ) + Ω2(3Q 11qs + Q 10) + 3Q 13qs + Q 12) 2 − qd(t )(Ω(2Q 8qs + Q 7)pḋ (t ) + (2Q 5qs + Q 4 )qḋ (t )−(2Q 2qs + Q 1)q¨d(t ) − Ω2(3Q 11qs2 + 2Q 10qs + Q 9)) +

− (Q 11Ω2 + Q 13)qd(t )3 − Ω(Q 8qs2 + Q 7qs + Q 6)pḋ (t ) − (Q 5qs2 + Q 4qs + Q 3)qḋ (t ) ωp2pd (t )

(10)

+ p¨d (t ) = − (Ω2(3P11ps + P10) + P5pḋ (t ) + P2p¨d (t ) + P8Ωqḋ (t ) + 3P13ps − P12)pd (t )2 − pd (t ) (Ω(2P8ps + P7)qḋ (t ) + (2P5ps + P4 )pḋ (t )+(2P2ps + P1)p¨d (t ) + Ω2(3P11ps2 + 2P10ps + P9)) − (P11Ω2 + P13) pd (t )3 − Ω(P8ps2 + P7ps + P6)qḋ (t ) − (P5ps2 + P4ps + P3)pḋ (t )

(11)

The coefficients ωp, ωq, Qi and Pi for i = 1, 2, … , 14 are given in Appendix B. The detailed derivation of the modulation equations is provided in Appendix C using the method of multiple scales.

3. Mechanical–thermal noise analysis 3.1. Nonlinear analysis The mechanical–thermal force is associated with damping and appears in the same order as damping and the excitation forces, i.e in the third-order in the multiple-scales analysis. Given that the micro/nano-systems present very light damping and thus act as narrowband filters, we assume the noise force appears at the frequency of oscillation (which is equal to the effective natural frequency). For the highest sensitivity, the excitation frequency, the effective natural frequency in the drive direction, and the frequency of the sense direction are matched [11]. Therefore, the modified form of the third-order problem including the noise force in the sense direction is given by

∂ 20pd + pd ωp2 = Tncos(Ωe t ) − P13pd3 − P2∂ 20pd pd2 − 2P12pd pd − 6ps P13pd pd − 2P1∂1∂0pd pd − 4ps P2∂1∂0pd pd 3

3

1

1

1

2

1

2

1

1

1

1

1

− P1∂ 20pd pd − 2ps P2∂ 20pd pd − ∂12pd − P3∂0pd − ps P4 ∂0pd − ps2 P5∂0pd − ΩP6∂0qd − Ωps P7∂0qd 2

−

1

Ωps2 P8∂0qd 1

2

1

1

− 2∂2∂0pd − 2∂1∂0pd − 1

2

1

P1pd ∂ 20pd 2 1

1

−

1

1

1

2ps P2pd ∂ 20pd 2 1

(12)

^ where Tn = Fn(ϕ(1) + e ϕ′(1)) indicates the noise force. Following the same procedure as outlined in Appendix C, the third-order problem develops to 1 ∂ 20pd + ωq2 pd = ( 2 Tn eiσt2 − i(P8ps2 + P7ps + P6)ΩωqAq + ei δ t2((5H1P1 + 2H2P1+(10H1ps + 4H2ps + 3)P2)ωq2A¯ p 3

3

A p2 − (2H1P12+4H2P12 + (6H1ps + 12H2ps + 3)P13)A¯ p A p2 −iωq((P5ps2 + P4ps + P3)A p + 2Ȧ p )))eit0ωq + C . C .

(13)

And, the modulation equations in polar coordinates become

⎛1 ⎞ 3 1 ⎜ (G1 + 2G2)(3Q 13qs + Q 12⎟ + Q 13 − ωq2((5G1 + 2G2)(2Q 2qs + Q 1) + 3Q 2)) ⎠ 8 ⎝4 8 1 Ωωqa p(qs(Q 8qs + Q 7) + Q 6)sin(η2) + ωqaq(η1̇ − σ ) = Q 14VACVwcos(η1) 2 ⎛1 ⎞ 1 3 ⎜ (H1 + 2H2) + P13cos(η2⎟ − ωq2cos(η2)((5H1 + 2H2)(2P2ps + P1) + 3P2)) ⎠ 8 ⎝4 8 aq3 −

a p3 + ωq((δ − σ + η1̇ − η2̇ )cos(η2)a p − sin(η2)ȧp − =

1 a p(ps (P5ps + P4 ) + P3)sin(η2)) + (3P13ps + P12)cos(η2) 2

1 Tncos(η1) 2

(15)

⎛1 ⎞ 1 Ωωqa p(qs(Q 8qs + Q 7) + Q 6)cos(η2) + ωq⎜ aq(qs(Q 5qs + Q 4 ) + Q 3) + ȧq⎟ = Q 14VACVwsin(η1) ⎝2 ⎠ 2 ⎛1 ⎞ 1 3 ⎜ (H1 + 2H2) + P13sin(η2⎟ − ωq2sin(η2)((5H1 + 2H2)(2P2ps + P1) + 3P2)) ⎠ 8 ⎝4 8 ⎛1 ⎞ 1 a p3 + ⎜ (ps (P5ps + P4 ) + P3⎟cos(η2)a p + (δ − σ + η1̇ − η2̇ )sin(η2)a p + cos(η2)ȧp )ωq + (3P13ps + P12)sin(η2) + ⎠ ⎝2 2 Ω ωqaq(ps (P8ps + P7) + P6) =

(14)

1 Tnsin(η1) 2

(16)

(17)

S.A.M. Lajimi et al. / Mechanical Systems and Signal Processing 83 (2017) 163–175

167

In Cartesian coordinates system, while the first three modulation equations remain the same as Eq. (C.15), the last equation becomes

yṗ = f4

(18)

where

f4 =

1 1 1 Tn + yp (−ps (P5ps + P4 ) − P3) − Ωyq (ps (P8ps + P7) + P6) + (σ − δ )x p 2 ωq 2 2 ⎛ y2 (ω 2((5H + 2H )(2P p + P ) + 3P ) − 2(H + 2H )(3P p + P ) − 3P ) ⎞ q 1 2 2 s 1 2 1 2 13 s 12 13 p ⎟ + x p⎜⎜ ⎟ 8 ω q ⎝ ⎠ +

x p3(ωq2((5H1 + 2H2)(2P2ps + P1) + 3P2) − 2(H1 + 2H2)(3P13ps + P12) − 3P13) 8ωq

For practical purposes, it is desired to avoid a multi-valued response and the consequent instability. Therefore, the AC voltage is limited to the range resulting in a single-valued response in the frequency–response results.

3.2. Linear analysis The objective is to characterize the effect of MTN on the displacement in the sense direction for beam-based gyroscopes. The displacement in the drive axis is sufficiently large that the MTN effect is negligible in this direction. The virtual work by the thermal noise is computed at the center of end-rigid-body as

δWn = Fnδvc

(19)

where Fn represents the generalized noise force and vc the displacement at the center of end-mass. Including the thermal noise work, discretizing displacements using single-mode approximations, and computing Lagrange's equations, gives a similar equation to (8) with an additional generalized (modal) noise force term. The equations are linearized around the static position using a Taylor's expansion. The MTN is associated with the damping and appears in the dynamic equation. The system of equations are expressed as

q¨d(t ) + Dqqḋ (t ) + Kqqd(t ) = − Fqcos(tΩe )

(20)

^ p¨d (t ) + ΩCpqḋ (t ) + Dppḋ (t ) + Kppd (t ) = Fn(eϕ′(1) + ϕ(1))

(21)

where

Kq = Dq =

ΔΓ ″(qs(eψ ′(1) + ψ (1)) − 1)3 + 2eνVw2 (eψ ′(1) + ψ (1))2 Δ(qs(eψ ′(1) + ψ (1)) − 1)3(αΓ ′ + Γ + M (eψ ′(1) + ψ (1))2 + J22 ψ ′(1)2) Δ(αΓ ′ + Γ + M (eψ ′(1) + ψ (1))2 + J22 ψ ′(1)2)

2eνVACVw(eψ ′(1) + ψ (1)) Fq = − Δ(qs(eψ ′(1) + ψ (1)) − 1)2(αΓ ′ + Γ + M (eψ ′(1) + ψ (1))2 + J22 ψ ′(1)2) Kp = Dp = ^ Fn =

(22)

cq (23) (24)

Δ3 Λ″(ps (eϕ′(1) + ϕ(1)) − 1)3 + 2eνVv2(eϕ′(1) + ϕ(1))2 Δ(ps (eϕ′(1) + ϕ(1)) − 1)3(αΔ2 Λ′ + M (eϕ′(1) + ϕ(1))2 + J33 ϕ′(1)2 + Λ)

(25)

cp Δ(αΔ2 Λ′ + M (eϕ′(1) + ϕ(1))2 + J33 ϕ′(1)2 + Λ)

(26)

Fn 6eLBgv α 2

(27)

In Eq. (20), the effect of the Coriolis force on the drive direction is neglected which is acceptable by virtue of its small magnitude compared with the driving force. Inserting a solution of the form qd(t ) = aqcos(Ωe t − θ ) into Eq. (20), separating the coefficients of cos(Ωe t ) and sin(Ωe t ), and equating the coefficients to zero gives the following relations for the phase and amplitude:

θ = − arctan(

DqΩe Ωe2 − Kq

) (28)

168

S.A.M. Lajimi et al. / Mechanical Systems and Signal Processing 83 (2017) 163–175

aq =

Fq Dq2Ωe2

+ (Ωe2 − Kq)2

(29)

or

aq =

−2eνVACVw(eψ ′(1) + ψ (1)) Δ Dq2Ωe2 + (Ωe2 − Kq)2 (qs(eψ ′(1) + ψ (1)) − 1)2(αΓ ′ + Γ + M (eψ ′(1) + ψ (1))2 + J22 ψ ′(1)2)

(30)

Eq. (30) provides the necessary tool to measure the maximum amplitude for most of practical purposes where the response is limited to the linear range. According to the sense mode equation, (21), the Coriolis force is given by Fc = − ΩCpqḋ (t ). The maximum displacement rate in the drive direction, qḋ (t ) is Ωe aq . The amplitude is given in Eq. (30). The Ωn is determined by equating the thermal noise force and the Coriolis force, thus

^ Cp Ωn Ωe aq = Fn(eϕ′(1) + ϕ(1)) Solving Eq. (31) for

Ωn = −

(31)

Ωn gives ^ Fn(eϕ′(1) + ϕ(1))(αΔ2 Λ′ + M (eϕ′(1) + ϕ(1))2 + J33 ϕ′(1)2 + Λ)

aqΩe(2M (eψ ′(1) + ψ (1))(eϕ′(1) + ϕ(1)) + (−J11 + J22 + J33 )ψ ′(1)ϕ′(1) + 2Π )

(32)

Inserting Eq. (30) into Eq. (32) gives the closed-form relation for computing the noise-equivalent angular rotation rate, Ωn, for a beam-rigid body gyroscope (rotation rate sensor). Operating the gyroscope near the effective natural frequency in2 2 creases the sensitivity of the device, that is the sense amplitude increases. Therefore, for Ωe near Kq , the oscillation amplitude, Eq. (29), aq becomes aq =

Fq

Dq Ωe

.

The MTN equivalent rotation rate for beam-rigid-body gyro is thus given by

Ωn = −

^ ΓωdFn(eϕ′(1) + ϕ(1))(qs(eψ ′(1) + ψ (1)) − 1)2( 2eνVACQ dVw(eψ ′(1) + ψ (1)))

(33)

where

^ ^ ^ ^ ( = (αΓ ′ + Γ + M (e^ψ ′(1) + ψ (1))2 + J22 ψ ′(1)2)(αΛ′ + M (e^ϕ′(1) + ϕ(1))2 + J33 ϕ′(1)2 + Λ) ^ ^ ^ ^ ) = (2M (e^ψ ′(1) + ψ (1))(e^ϕ′(1) + ϕ(1)) − (J11 − J22 − J33 )ψ ′(1)ϕ′(1) + 2Π ) where the damping has been replaced with mq ωd /Q d with mq being the coefficient of the modal acceleration terms, that is γ1,1 in Appendix A. Eq. (33) derived here for the first time includes the effect of mode shapes and the end-mass parameters. The equation is of prime importance in the design procedure of beam-based rotation rate sensors. The noise equivalent rotation rate strongly depends on the eccentricity (e) and the rotary inertia of the end-mass indicating the importance of taking the dimensions of the end-mass into account. It is noted that the Ωn is inversely related to the quality factor in the drive direction, Qd and the AC voltage VAC . The static deflection qs is associated with the DC load Vw.

4. Numerical analysis To compute the noise equivalent (angular) rotation rate (Ωn) using the slow system, the pure noise induced displacement signal in the sense direction is computed using Eqs. (14)–(17). To this end, the angular rotation rate, Ω, and the excitation force, VAC , are set to zero and the steady-state response in the sense direction under noise excitation is computed. Having found the noise response, the equations of pure signal, no noise input, that is equations (C.15), are solved for Ωn, aq, η1, and η 2. The nonlinear analysis is employed to identify the onset of multi-valuedness where an unstable branch appears in the response, and the slow system including noise force is solved for pure noise response and the Ωn. The linear approximation of the noise equivalent rotation rate is computed from Eq. (32) or (33). For the same set of AC voltage, the linear estimation of Ωn is computed and compared with the result of nonlinear noise analysis. In practice, the linear analysis provides a reliable tool for MTN analysis and the optimization of the device for lower Ωn. The structural specification of the microbeam gyroscope are given in Table 1. The proposed design adheres to the standard SOIMUMPS process [12]. Initially, we perform a nonlinear analysis computing Eqs. (C.12) and (C.13), having found the solutions of equilibrium points of Eqs. (C.8)–(C.11). Having done such an analysis, see Figs. 2a and b, we find the appropriate region where the

S.A.M. Lajimi et al. / Mechanical Systems and Signal Processing 83 (2017) 163–175

169

Table 1 Dimensions. L

a

b

LM

aM

bM

LH

aH

bH

h

g

400 μm

25 μm

25 μm

80 μm

10 μm

25 μm

10 μm

10 μm

25 μm

25 μm

2 μm

response is linear and single-valued and therefore more desired for practical purposes. The red frequency–response plot indicating a higher center frequency is obtained by setting the DC voltages equal to 17.1 V and 3 V in the drive and sense directions. The middle (blue) curve is for the DC voltages 27.9 V and 22.4 V and the last curve operating around the lowest center frequency corresponds to the DC voltages 34.5 V and 30.3 V in the drive and sense directions, respectively. In Figs. 3a and b, the Ωn is plotted versus the AC load. The respective quality factors are 1000 in the drive direction and 1000 and 10 000 in the sense direction. It is observed that the linear analysis (the solid-line) and the nonlinear analysis (the dashed-line) are in excellent agreement while the nonlinear approximation is always larger than the linear one. It is noted that the Ωn is reduced with an increasing AC load. This presents the importance of operating at the highest possible AC voltage. On the other hand, increasing the quality factor has a positive effect on reducing Ωn. In Tables 2 and 3, the corresponding Ωn for the highest AC voltages for each set of results in Figs. 3a and b are presented. Defining the error (%) parameter as follows:

Percentage Error =

ΩnNonlinear − ΩnLinear ΩnNonlinear

× 100 (34)

we plot the error parameter in Figs. 4a and b. The percentage error is not significantly affected by the quality factor and varies mainly by the DC and AC voltages. Therefore, it is recommended to design a system with higher quality factor considering the desired bandwidth. An approach to reduce the damping and increase the quality factor is to package the sensor in vacuum and reduce the effect of squeeze film damping. The higher the bandwidth is the lower the quality factor. The static deflection is determined by the DC loads and affects the Ωn, see Eq. (32). It is seen that operating at higher DC voltages and therefore lower effective natural frequencies, increases Ωn. The conclusions are valid for any VAC , that is design a system with lower damping, higher quality factor, and operating at the lower VDC results in lower Ωn. The gyro achieves 0.86926°/s, 0.08694°/s, and 0.00881°/s noiseequivalent rotation rate for 105.232 Hz, 10.5232 Hz, and 1.05232 Hz bandwidth operating at voltages Vw = 17.1 V , Vv = 3 V , VAC = 2.8 V .

5. Conclusion A nonlinear approach using the method of multiple scales has been introduced to provide an accurate estimate of the noise-equivalent rotation rate. Gyroscopes based on a cantilever beam with a rigid-end-body and driven by a combined AC/ DC voltage have been our focus here. The linear approximation of the noise-equivalent rotation rate of the beam-based gyroscope designs considered here has been derived for the first time and provides a deeper understanding of the influencing parameters.

Fig. 2. The frequency–response plots for a set of DC voltages. The response remains linear and single-valued for the considered set of parameters. (For interpretation of the references to color in this figure caption, the reader is referred to the web version of this paper.)

170

S.A.M. Lajimi et al. / Mechanical Systems and Signal Processing 83 (2017) 163–175

Fig. 3. The noise equivalent (angular) rotation rate (Ωn) versus the AC voltage. The AC load is limited by the onset of multi-valuedness of the response. The solid-line indicates the results of nonlinear analysis and the dashed-line the linear analysis.

Table 2 Noise equivalent rotation rate (Ωn) for Q d = 1000 , Q s = 1000 .

Parameter

VAC (V) BW (Hz) Ωn (°/s) (Nonlinear analysis) Ωn (°/s) (Linear analysis)

Vw = 17.1 V

Vw = 27.9 V

Vw = 34.5 V

Vv = 3 V

Vv = 22.4 V

Vv = 30.3 V

2.8 105.232 0.869257 0.653696

1 104.929 1.45265 1.11747

0.6 104.672 1.88359 1.48059

As a one of the principal contributions in this paper we have provided, for the first time, a new method to compare the linear and the nonlinear perturbation methods of mechanical–thermal noise analysis. We have developed a new approach to compute the noise-equivalent rotation rate using the method of multiple scales. In addition it has been systematically demonstrated that the important parameters affecting the noise-equivalent rate for an actual gyroscope are, the bias voltage, the AC voltage, and the quality factor. In obtaining the noise-equivalent rotation

S.A.M. Lajimi et al. / Mechanical Systems and Signal Processing 83 (2017) 163–175

171

Table 3 Noise equivalent rotation rate (Ωn) for Q d = 1000 , Q s = 10 000 .

Parameter

VAC (V) BW (Hz) Ωn (°/s) (Nonlinear analysis) Ωn (°/s) (Linear analysis)

Vw = 17.1 V

Vw = 27.9 V

Vw = 34.5 V

Vv = 3 V

Vv = 22.4 V

Vv = 30.3 V

2.8

1

0.6

10.5232 0.0869372

10.4929 0.145678

10.4672 0.188383

0.0653696

0.110924

0.148059

Fig. 4. The percentage error versus the AC voltage. Changing the quality factor has an insignificant effect on the percentage error.

rate, we have shown that the higher the quality factor is, the lower the noise equivalent rotation rate will be. Similarly, the higher the bias voltage, and thus the larger the static deflection position are, the higher noise equivalent rotation rate will be. And finally as the AC voltage increases, the noise equivalent rotation rate reduces and the difference between the linear

172

S.A.M. Lajimi et al. / Mechanical Systems and Signal Processing 83 (2017) 163–175

and nonlinear approaches becomes larger. The nonlinear perturbation approach presented here provides a means to compute the noise-equivalent rotation rate regardless of the operating AC voltage. This allows the determination of the noise-equivalent rotation rate in the nonlinear regime that exists at high AC voltages. While the proposed design already reaches to a very low noise-equivalent rotation rate, given that the considered quality factor is not very high, if one chooses to operate at higher AC voltages the gyro achieves better performance.

Appendix A. The coefficients of the discretized model The coefficients of the discretized model are:

^ γ1,1 = Δ(Γ + α Γ ′ + J

^ ψ ′(1)2 + M (ψ (1) + e^ ψ ′(1))2),

M ηη

⎛ ⎞ ^ ^ ^ ^ γ1,2 = Δ(2 Π − ⎜⎜ J − J − J ⎟⎟ϕ′(1)ψ ′(1) + 2 M (ϕ(1) + e^ ϕ′(1))(ψ (1) + e^ ψ ′(1))), M M M ⎝ ξξ ηη ζζ ⎠ ^ ^ ′ ′( ′( γ1,3 = Δ(Π + αΠ + J ψ 1)ϕ 1) + M (ϕ(1) + e^ ϕ′(1))(ψ (1) + e^ ψ ′(1))), M ηη

^ γ1,4 = Δ(Γ − α Γ ′ + M (ψ (1) + e^ ψ ′(1))2 − J

^ ψ ′(1)2 + J

M ξξ

ψ ′(1)2),

M ζζ

γ1,5 = (ψ (1) + e^ψ ′(1)) ^ γ2,1 = Δ(Λ + α Δ2 Λ′ + J

ϕ′(1)2 + M (ϕ(1) + e^ ϕ′(1))2),

M ζζ

⎛ ⎞ ^ ^ ^ ^ γ2,2 = Δ(2 Π + ⎜⎜ J + J − J ⎟⎟ϕ′(1)ψ ′(1) + 2M (ψ (1) + e^ ψ ′(1))(ϕ(1) + e^ ϕ′(1))), ⎝M ηη M ζζ M ξξ ⎠ ⎛ ⎞ ^ ^ ^ γ2,3 = Δ(Π − αΠ′ + ⎜⎜ J − J ⎟⎟ψ ′(1)ϕ′(1) + M (ϕ(1) + e^ ϕ′(1))(ψ (1) + e^ ψ ′(1))), ⎝M ζζ M ξξ ⎠ ^ ^ ^ γ2,4 = Δ(Λ − α Δ2 Λ′ + M (ϕ + e^ ϕ′(1))2 − J ϕ′(1)2 + J ϕ′(1)2), M ξξ

M ηη

γ2,5 = ϕ(1) + e^ϕ′(1) and

Γ=

∫0

Λ=

∫0

1

1

ψ (x^)2dx^ ,

Γ′ =

∫0

ϕ(x^)2dx^ ,

Λ′ =

∫0

1

1

ψ ′(x^)2dx^ ,

Γ″ =

∫0

ϕ′(x^)2dx^ ,

Λ″ =

∫0

1

1

ψ ″(x^)2dx^ ,

ϕ″(x^)2dx^ ,

Π=

∫0

1

ψ (x^)ϕ(x^)dx^

The nondimensionalization has been achieved by defining

J ^ = w, w gw

v v^ = , gv

x x^ = , L

t t^ = , κ

M ^ M= , mL

^ = κ Ω, Ω

Appendix B. The coefficients of the dynamic equations The coefficients

ωp, ωq, Qi and Pi for i = 1, 2, …, 14 are given by

e e^ = , L

^ J M ξξ

=

M ξξ

, mL3

J ^ J M ηη

=

M ηη

, mL3

J ^ J M ζζ

=

M ζζ

mL3

S.A.M. Lajimi et al. / Mechanical Systems and Signal Processing 83 (2017) 163–175

Γ ″(3γ1,5qs − 1)

ωq2 =

γ1,1(γ1,5qs − 1)

Q4 = Q8 =

Q1 = −

,

−2Γγ1,5ωd

, Q d(γ1,5qs − 1)2 2 γ1,2γ1,5 2

γ1,1(γ1,5qs − 1)

Q 12 = −

Q5 =

,

2γ1,5Γ ″

Q9 =

, γ1,1(γ1,5qs − 1)2

2γ1,5

,

2

(γ1,5qs − 1) 2 Γγ1,5 ωd

2

Q d(γ1,5qs − 1) γ1,4

2

γ1,1(γ1,5qs − 1)

Q 13 =

2 γ1,5

Q2 =

(γ1,5qs − 1)2

,

Q6 =

,

Q 10 =

2 γ1,5 Γ″

γ1,1(γ1,5qs − 1)2

,

,

Q3 =

γ1,2 2

γ1,1(γ1,5qs − 1)

Q7 =

,

Q 11 =

2

γ1,1(γ1,5qs − 1)

,

−2γ1,2γ1,5

,

2γ1,3γ1,5

Q 14 =

Γωd Q d(γ1,5qs − 1)2

173

γ1,1(γ1,5qs − 1)2

,

2 −γ1,3γ1,5

γ1,1(γ1,5qs − 1)2

,

eνγ1,5 γ1,1(γ1,5qs − 1)2

for the drive equation, and

ωp2 = P4 = P8 = P12 =

Λ″(3γ2,5ps − 1) γ2,1(γ2,5ps − 1) −2Λγ2,5ωs

,

P1 = − P5 =

, γ2,1(γ2,5ps − 1)2

P9 =

Q s(γ2,5ps − 1) 2 −γ2,2γ2,5

−2γ2,5Λ″ 2

γ2,1(γ2,5ps − 1)

,

, (γ2,5ps − 1)2 2 Λγ2,5 ωs

,

2

2γ2,5

Q s(γ2,5ps − 1)2

P2 =

,

P6 =

γ2,4

, γ2,1(γ2,5ps − 1)2

P13 =

(γ2,5ps − 1)2

2

,

,

−γ2,2 γ2,1(γ2,5ps − 1)2

P10 =

2 γ2,5 Λ″

γ2,1(γ2,5ps − 1)

2 γ2,5

P3 = ,

2γ2,3γ2,5

, γ2,1(γ2,5ps − 1)2

P14 =

Λωs Q s(γ2,5ps − 1)2

P7 =

2γ2,2γ2,5 γ2,1(γ2,5ps − 1)2

P11 =

,

2 −γ2,3γ2,5

γ2,1(γ2,5ps − 1)2

,

eνγ2,5 γ2,1(γ2,5ps − 1)2

for the sense equation.

Appendix C. Method of multiple-scales To solve for the displacements using the method of multiple-scales, we replace the flexural displacements of the structure in the sense and drive directions with their third-order approximations [13]

qd(t ) = ϵ qd (t0, t1, t2) + ϵ2qd (t0, t1, t2) + ϵ3qd (t0, t1, t2)

(C.1)

pd (t ) = ϵ pd (t0, t1, t2) + ϵ2pd (t0, t1, t2) + ϵ3pd (t0, t1, t2)

(C.2)

1

1

2

2

3

3

where the time scales are defined as t0 = t , t1 = ϵ t , and t2 = ϵ2t. The parameter ϵ indicates the order of the terms. The forcing term, VAC , the rotation rate Ω, and the damping terms such that they appear in the third-order equations. Inserting (C.1) and (C.2) into (10) and (11) and collecting the coefficients of like powers of ϵ results in a sequence of three sets of equations which are solved consecutively to yield a third-order uniform expansion of the displacements [8]. An internal and an external detuning parameters are introduced to describe the primary and internal resonances as

Ωe = ωq + ϵ2σ

(C.3)

ωp = ωq + ϵ2δ

(C.4)

Eq. (C.3) characterizes the nearness of the excitation frequency Ωe to the drive natural frequency ωq, and Eq. (C.4) characterizes the nearness of two modal frequencies in the drive and sense directions, that is ωq and ωp. Substitution of solutions of the first-order and second-order problems into the third-order problem and dropping terms associated with higher frequencies results in [8]: ∂ 20qd + ωq2 qd = (((5G1 Q 1 + 2G 2 Q 1 + 10G1 qs Q 2 + 4G 2 qs Q 2 + 3Q 2)ωq2−2G1 Q 12 − 4G 2 Q 12 − (6G1 qs + 12G 2 qs + 3)Q 13)A¯ q 3

3

Aq2 − i(Q 5qs2 + Q 4qs + Q 3)ωqAq + e iσt 2Q 14VACVw −i(Q 8qs2 + Q 7qs + Q 6 )ωq Ω A p e i δ t 2 − 2iωqȦq )e it0ωq + C . C .

(C.5)

∂ 20pd + ωq2 pd = (−i(P8ps2 + P7ps + P6 )ΩωqAq + e i δ t 2((5H1P1 + 2H2P1+(10H1ps + 4H2ps + 3)P2)ωq2A¯ p A p2 3

3

− (2H1P12+4H2P12 + (6H1ps + 12H2ps + 3)P13)A¯ p A p2 −iωq((P5ps2 + P4ps + P3)A p + 2Ȧ p )))e it0ωq + C . C .

(C.6)

where Ȧq and Ȧp represent the time derivative of amplitudes Aq (t2) and Ap (t2) with respect to the slow time scale t2 and C . C .

174

S.A.M. Lajimi et al. / Mechanical Systems and Signal Processing 83 (2017) 163–175

stands for the complex conjugate of its preceding term. The coefficients of eit0ωq give rise to the secular terms and therefore are set to zero. To this end, the amplitudes, Aq (t2) and Ap (t2), are expressed in polar form as

Aq (t2) =

1 aq(t2)eiθq(t2) , 2

A p (t2) =

1 a p(t2)eiθp(t2) , 2

(C.7)

substituted in Eqs. (C.5)–(C.6), the real and imaginary parts are separately equated to zero resulting in the following four modulation equations:

⎛1 ⎞ 3 1 ⎜ (G1 + 2G2⎟(3Q 13qs + Q 12) + Q 13 − ωq2((5G1 + 2G2)(2Q 2qs + Q 1) + 3Q 2)) ⎠ ⎝4 8 8 aq3 −

1 Ωωqa p(qs(Q 8qs + Q 7) + Q 6)sin(η2) + ωqaq(η1̇ − σ ) = Q 14VACVwcos(η1) 2

(C.8)

⎛1 ⎞ 3 1 ⎜ (H1 + 2H2⎟ + P13cos(η2) − ωq2cos(η2)((5H1 + 2H2)(2P2ps + P1) + 3P2)) ⎠ 8 ⎝4 8 a p3 + ωq((δ − σ + η1̇ − η2̇ )cos(η2)a p − sin(η2)ȧp −

1 a p(ps (P5ps + P4 ) + P3)sin(η2)) + (3P13ps + P12)cos(η2) = 0 2

(C.9)

1 1 Ωωqa p(qs(Q 8qs + Q 7) + Q 6)cos(η2) + ωq( aq(qs(Q 5qs + Q 4 ) + Q 3) + aq̇ ) = Q 14VACVwsin(η1) 2 2

(C.10)

⎛1 ⎞ 3 1 ⎜ (H1 + 2H2⎟ + P13sin(η2) − ωq2sin(η2)((5H1 + 2H2)(2P2ps + P1) + 3P2)) ⎠ 8 ⎝4 8 ⎛1 1 a p3 + ⎜ (ps (P5ps + P4 ) + P3)cos(η2)a p + (δ − σ + η1̇ − η2̇ )sin(η2)a p + cos(η2)ȧp )ωq + (3P13ps + P12)sin(η2) + ⎝2 2 ⎞ Ω ωqaq(ps (P8ps + P7) + P6⎟ = 0 ⎠

(C.11)

where η2(t2) = θp(t2) − θq(t2) + δ t2, η1(t2) = σ t2 − θq(t2), and the time derives of the amplitudes and phases aq̇ , ȧp , η̇1, and η̇2 are computed with respect to the slow time scale t2. The general solutions, that is Eqs. (C.1) and (C.2), therefore become

⎛1 1 ⎞ qd(t ) = aq(t )cos(Ωet − η1(t )) − ⎜ Q 2qs − Q 1⎟aq(t )2cos(2Ωet − 2η1(t )) ⎝3 6 ⎠ ⎛ ⎞ ⎛1 ⎞ 1 3 1 ⎜ Q 13qs + Q 12⎟cos(2 Ωet − 2η1(t )) − Q 13qs − Q 12 ⎟ ⎜ ⎝ ⎠ 1 2 6 2 2 ⎟a (t )2 + ⎜ Q 2qs + Q 1 + ⎜ ⎟ q 2 ωq2 ⎜ ⎟ ⎝ ⎠ pd (t ) = a p(t )cos(Ωet − η1(t ) + η2(t )) + +

(2P2ps ωp2 − 3P13ps + P1ωp2 − P12) 2ωp2

(−2P2ps ωp2 + 3P13ps − P1ωp2 + P12)cos(2Ωet − 2η1(t ) + 2η2(t )) 6ωp2

(C.12)

a p(t )2 a p(t )2

(C.13)

To find the periodic solutions of the structure about the static equilibrium, the equilibrium points of the modulation equations are found by setting aq̇ (t ) = 0, ap(t ) = 0, η1(t ) = 0, and η2(t ) = 0 in (C.8)–(C.11). To re-express the modulation equation in the Cartesian coordinates, the amplitudes are described in the form:

Aq (t2) =

1 (xq(t2) − iyq(t2))eiλq (t2) , 2

A p (t2) =

1 (x p(t2) − iyp(t2))eiλp (t2) 2

(C.14)

resulting in

xq̇ = f1 ,

yq̇ = f2 ,

ẋp = f3 ,

yṗ = f4

(C.15)

S.A.M. Lajimi et al. / Mechanical Systems and Signal Processing 83 (2017) 163–175

175

where

f1 =

xq2yq (2(G1 + 2G2)(3Q 13qs + Q 12) + 3Q 13) − ωq2(((5G1 + 2G2)(2Q 2qs + Q 1) + 3Q 2)) 8ωq + −

yq3 (ωq2(−((5G1 + 2G2)(2Q 2qs + Q 1) + 3Q 2)) + 2(G1 + 2G2)(3Q 13qs + Q 12) + 3Q 13)

yq2 (ωq2((5G1 + 2G2)(2Q 2qs + Q 1) + 3Q 2) − 2(G1 + 2G2)(3Q 13qs + Q 12) − 3Q 13)

−

+

8ωq

xq3(ωq2((5G1 + 2G2)(2Q 2qs + Q 1) + 3Q 2) − 2(G1 + 2G2)(3Q 13qs + Q 12) − 3Q 13) 8ωq

+ (δ − σ )yp − +

1 xq(−qs(Q 5qs + Q 4 ) − Q 3) 2

Q V V 1 1 1 xp Ω(qs(Q 8qs + Q 7) + Q 6) − σyq f2 = 14 AC w + σxq − Ωyp (qs(Q 8qs + Q 7) + Q 6) + yq (−qs(Q 5qs + Q 4 ) − Q 3) 2 1 ωq 2 2

+ xq( +

+

8ωq

)

f3 =

1 x p(−ps (P5ps + P4 ) − P3) 2

x p2yp (ωq2(−((5H1 + 2H2)(2P2ps + P1) + 3P2)) + 2(H1 + 2H2)(3P13ps + P12) + 3P13) 1 Ωxq(ps (P8ps + P7) + P6) + 2 8ωq

yp3 (ωq2(−((5H1 + 2H2)(2P2ps + P1) + 3P2)) + 2(H1 + 2H2)(3P13ps + P12) + 3P13) 8ωq

f4 =

1 y (−ps (P5ps + P4 ) − P3) 2 p

yp2 (ωq2((5H1 + 2H2)(2P2ps + P1) + 3P2) − 2(H1 + 2H2)(3P13ps + P12) − 3P13) 1 ) Ωyq (ps (P8ps + P7) + P6) + (σ − δ )x p + x p( 2 8ωq x p3(ωq2((5H1 + 2H2)(2P2ps + P1) + 3P2) − 2(H1 + 2H2)(3P13ps + P12) − 3P13) 8ωq

The stability of stationary (equilibrium) points is determined by computing and evaluating the eigenvalues of the Ja∂F where F = f1 f2 f3 f4 and x = xq yq xp yp , resulting in a 4 ×4 matrix evaluated at the cobian matrix defined by J = ∂x equilibrium points. Equilibrium (stationary) solutions with corresponding eigenvalues with negative real parts are asymptotically stable while others with at least one positive eigenvalue represent an unstable solution.

{

}

{

}

References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13]

T.J. Hofler, S.L. Garrett, Thermal noise in a fiber optic sensor, J. Acoust. Soc. Am. 84 (1988) 471. T.B. Gabrielson, Mechanical–thermal noise in micromachined acoustic and vibration sensors, IEEE Trans. Electron Dev. 40 (5) (1993) 903–909. V. Annovazzi-Lodi, S. Merlo, Mechanical–thermal noise in micromachined gyros, Microelectron. J. 30 (12) (1999) 1227–1230. S.D. Senturia, Microsystem Design, Kluwer Academic Publishers, Boston, 2001. M.-H. Bao, Micro Mechanical Transducers: Pressure Sensors, Accelerometers and Gyroscopes, vol. 8, Elsevier B. V, Amsterdam, 2000. R.P. Leland, Mechanical–thermal noise in mems gyroscopes, Sens. J. IEEE 5 (3) (2005) 493–500. D. Lynch, Vibratory gyro analysis by the method of averaging, in: Proceedings of 2nd St. Petersburg Conference on Gyroscopic Technology and Navigation, St. Petersburg, pp. 26–34. S.A.M. Lajimi, Design, modeling, and nonlinear dynamics of a cantilever beam-rigid body microgyroscope (Ph.D. thesis), University of Waterloo, 2013. S.A.M. Lajimi, G.R. Heppler, E. Abdel-Rahman, A new cantilever beam rigid-body MEMS gyroscope: mathematical model and linear dynamics, in: 2nd International Conference on Mechanical Engineering and Mechatronics (ICMEM'13), August 8–9, Toronto, Canada, 2013, pp. 177–1–177–6. S.A.M. Lajimi, E. Abdel-Rahman, G. Heppler, Modeling and sensitivity analysis of a mems vibratory rotation rate sensor, in: Proceedings of SPIE 2014, Sensors and Smart Structures Technologies for Civil, Mechanical, and Aerospace Systems 2014, vol. 9061, 2014, pp. 90612I–90612I–11. V. Apostolyuk, Theory and design of micromechanical vibratory gyroscopes, in: C.T. Leondes (Ed.), MEMS/NEMS Handbook Techniques and Applications, Springer, New York, NY, 2006, pp. 173–195. A. Cowen, G. Hames, D. Monk, S. Wilcenski, B. Hardy, Soimumps Design Handbook. S. Lajimi, G. Heppler, E. Abdel-Rahman, Primary resonance of a beam-rigid body microgyroscope, Int. J. Non-Linear Mech. 77 (2015) 364–375.