Analytical solution of squeeze film and thermoelastic damping in bulk resonators

Analytical solution of squeeze film and thermoelastic damping in bulk resonators

Sensors and Actuators A 297 (2019) 111530 Contents lists available at ScienceDirect Sensors and Actuators A: Physical journal homepage: www.elsevier...

5MB Sizes 0 Downloads 16 Views

Sensors and Actuators A 297 (2019) 111530

Contents lists available at ScienceDirect

Sensors and Actuators A: Physical journal homepage: www.elsevier.com/locate/sna

Analytical solution of squeeze film and thermoelastic damping in bulk resonators Yuanlong Cao, Ke Yu, Kai Feng, Hanqing Guan, Yihua Wu, Kai Zhang ∗ State Key Laboratory of Advanced Design and Manufacturing for Vehicle Body, Hunan University, Changsha, 410082, China

a r t i c l e

i n f o

Article history: Received 17 May 2019 Received in revised form 26 July 2019 Accepted 29 July 2019 Available online 15 August 2019 Keywords: MEMS Squeeze film damping Thermoelastic damping Quality factor Improved Fourier series method

a b s t r a c t Bulk micro-resonators exhibit potential applications for timing and frequency control based on their highquality factor and high frequency. However, performance of bulk resonators is difficult to predict due to the multi-physics nature of system including electrical actuation and vibration behavior. An analytical model of in-plane vibration of bulk resonators that consider effects of squeeze film damping (SFD) and thermoelastic damping (TED) is presented by a modified Fourier series method. The squeeze gas pressure, which generates the equivalent elastic and damping coefficients of gas, is determined by the Reynolds equation using a Green’s function method. The anchor, which links the fixed position and resonator, and the linear electrostatic force are simplified as equivalent but separate springs. Vibration displacement and temperature distributions are calculated by Fourier series coefficients. Numerical results of the effects of SFD and TED on the resonator are solved with a parametric analysis and verified by published data and simulation values of COMSOL® . The coupling effect of SFD and TED and their relative contribution on the bulk resonator are also evaluated. The presented model can guide to explore characteristics of bulk resonators in the presence of SFD and TED. © 2019 Elsevier B.V. All rights reserved.

1. Introduction Silicon micro-resonators have attracted wide attention in recent decades due to their high-quality factor, high frequency, and compatibility with modern integrated circuits, which is applied to accelerometers [1], sensors and timing devices [2]. Among the key parameters, the damping dissipation of resonators is a significant factor that influences system performance and development. For instance, in resonant sensors, high damping will suppress resonant response and reduce the sensitivity and resolution of system [3]. In most resonators, the major mechanisms of energy dissipation include squeeze film damping (SFD), thermoelastic damping (TED), and anchor and intrinsic losses [4,5]. Anchor loss is the vibration energy of resonators dissipated through their supports [6], which can be minimized by an excellent design of structure and fixing anchors in reasonable positions, such as nodes of mode shape [7,8]. Intrinsic loss, on the other hand, cannot be eliminated because they depend on the structural material and geometric property. Moreover, the system damping is dominated by intrinsic damping when operating pressure is near vacuum [9]. Squeeze film damping is caused by the energy dissipation of gas viscous damping when the contained gas is sucked into (or squeezed out of) the control domain. SFD can be modeled by the Reynolds equation [10] in which rarefied effects of the operating pressure are accounted for the variable viscosity [11,12]. Darling, Hivick and Xu [13] utilized Green’s function method to determine the elastic and damping components of the squeeze gas in cases of arbitrary venting conditions. Pandey and Pratap [14] provided solutions for the elastic and damping components of cantilever resonators around flexural modes. The predicted values showed that the quality factor is improved obviously because the SFD fells by 84% from the first mode to the second mode. Nayfeh and Younis [3] utilized perturbation method to determine pressure distribution in cases of 2D plate deflection and simulated the dynamic characteristics of beam resonators to predict the quality factor of SFD. The results showed that air damping is dominant when the operating pressure is high. Li, Fang and Wu [15] utilized molecular dynamics (MD) method to calculate the squeeze-film damping in free molecular regime. The energy dissipation is caused by the collisions between the microplate and the molecules of squeeze air.

∗ Corresponding author. E-mail addresses: [email protected] (Y. Cao), [email protected] (K. Zhang). https://doi.org/10.1016/j.sna.2019.111530 0924-4247/© 2019 Elsevier B.V. All rights reserved.

2

Y. Cao et al. / Sensors and Actuators A 297 (2019) 111530

Nomenclature Lateral area of the bulk resonator A cseff Effective damping coefficient of pressure per unit area Fs Effective squeezed-film force per unit area Fe Effective eletromechanical force per unit area Nominal gap spacing go H Thickness of the bulk resonator Ix0 ,In Number of elastic point support on edge at x = 0 and number of element Normal support stiffness of T-shape tether per unit area kT kt Tangential support stiffness of T-shape tether per unit area ke Effective stiffness of eletromechanical force per unit area kseff Effective elastic coefficient of pressure per unit area Pa Ambient pressure Entropy generation per unit volume S T Temperature function U, V Displacement function in the x direction and y direction DC voltage Vdc Ener Stored energy Dissipated energy of squeezed-film damping Es Ether Dissipated energy of thermoelastic damping ε0 Permittivity of free space εr Relative permittivity of the gap medium respect to the free space Angular frequency ω ˇl Hypothetical vibrational amplitude eff Equivalent viscosity of the fluid  Poisson’s ratio Coefficient of thermal expansion ˛  Thermal conductivity Cv Heat capacity A, B, C, E, F, H, K , M, P and Q Coefficient matrix of governing equation and boundary conditions

Thermoelastic damping (TED), which is first realized by Zener [16], is a type of irreversible intrinsic damping that is caused by the flow dissipation of temperature when the compression or expansion of the structure generate thermal gradient to flow inside structure [17]. Many analytical models have recently proposed the Fourier heat conduction theory for analyzing TED in resonators [18,19]. Prabhakar and Vengallatore [20] presented the expression of TED in infinite series by a Green’s function where the equation of heat conduction is solved by flexural mode shapes. Duwel, Candler, Kenny and Varghese [21] presented two approaches to analyze TED in micromechanical resonators, namely, the fully coupled thermomechanical equations and the uncoupled thermal and mechanical dynamic equations. Jiang and Li [22] calculated the thermal conduction equation with the deflection function of the clamped microplates and presented an analytical expression of TED. Guo, Yi and Pourkamali [23] applied an eigenvalue formulation and finite element method to study the effect of vent microplates on thermoelastic damping in which the quality factor of attached vent sections is 2000 more times higher than that of primary structure. Micromechanical resonators are mainly divided into flexural and bulk resonators [24]. Typical flexural resonators are slender and cantilever beams, and their vibrations are only in the thickness direction [3,25]. In fact, flexural resonator development is limited by relatively low natural frequency and high energy consumption [9]. However, bulk resonators employ the whole microplate as a massspring–damper system in which structural vibrations are only in the in-plane direction. They can yield high-resonant frequency and high Q to improve phase noise performance and measurement resolution of resonant frequencies [26]. Hao, Pourkamali and Ayazi [27] explored the dynamic characteristics of disk resonators and offered predictions by constructing a mechanical model. From the presented model, the electrical circuit parameters were to discuss the temperature coefficients of frequency and performance. Grogg and Tekin et al. [28] studied the in-plane vibrations of bulk resonators that consisted of multiple parallel beams. Then, the experimental behavior of bulk resonators is compared with those of disk and plate resonators. Lee and Seshia [29] presented a technique to excite the different resonant modes but with similar electrode configurations of a bulk resonator. The function was implemented by controlling the electrode position of the DC voltage, and then the equivalent circuit model was proposed to extract the motion resistance, resonant frequency, and quality factor of the bulk resonator [30]. However, most of the studies on resonators were based on experiments or equivalent circuit models. Numerical methods, which are widely used to calculate the performance of flexural resonators, have not been applied to in-plane vibrations of bulk resonators. Moreover, the performance predictions for SFD and TED of bulk resonators have not been fully explored because of the mechanical complexity of in-plane vibrations, squeeze film, and temperature. Nonetheless, for the in-plane vibrations of bulk resonators, several technical papers are available in the literature. The energy-based methods, such as the Rayleigh–Ritz technique, were used to investigate the in-plane vibrations of thin plates with classical boundary conditions. Gorman [31] introduced a superposition method to analyze the in-plane vibrations of thin plates with general and elastic boundary conditions. Du and Li et al. [32] developed a modified Fourier series method for in-plane vibrations with elastically restrained edges. Then, the approach was applied to the vibration characteristics of thin plates with multiple elastically point-supported edges [33]. Liu and Xing [34] proposed a separation of variables method for in-plane vibrations of orthotropic plates in cases wherein at least two of its edges are the simple support.

Y. Cao et al. / Sensors and Actuators A 297 (2019) 111530

3

Fig. 1. Schematic of an electrically actuated bulk resonator.

Fig. 2. Schematic of the change in SFD and TED caused by electrical loading.

In this paper, an analysis model is presented to study the effect of squeeze film and thermoelastic damping on bulk resonators. In-plane vibration of the bulk resonator, which is described as a thin plate with three elastic constraints, is analyzed by a modified Fourier series method. The expression of gas pressure determined by the linearized Reynolds equation is derived by the Green’s function method. Then, the equivalent elastic and damping coefficients of gas are calculated by pressure distribution. The anchor, which links the fixed position to the resonator, and the linear electrostatic force are also simplified as equivalent but separate springs. The influence of TED on bulk resonator is investigated by that temperature distribution, which is expressed by expansion Fourier coefficients of in-plane displacements, participates in equilibrium equations and boundary conditions. Numerical results of the effects of SFD and TED are solved with a parametric analysis and verified by experimental data and simulation values of COMSOL® . The coupling effect of SFD and TED and their relative contributions on the bulk resonator are also evaluated. 2. Theoretical formulation A bulk resonator is described as a combination of a square plate and eight fixed electrodes in Fig. 1. The square plate is restrained by T-shape anchors at four corner points [24]. Eight electrodes are fixed around the plate and divided into driving and sensing electrodes. The isotropic thermoelastic solid is excited by an electric load along x direction, which is composed by the voltage difference between the resonator with DC voltage Vdc and the drive electrode with AC voltage Vac [35,36]. For in-plane vibration of the bulk resonator, the z-axis is identified on the basis of all uniform strains and no loads [21]. When the electrostatic driving force generated by the capacitive transducer is applied, structural expansion or compression of the bulk resonator change gas gap between bulk resonator and electrodes. Any gap fluctuation will generate squeezing film forces in Fig. 2 that includes two components: viscous damping force when the gas is sucked into or squeezed out of the plate domain and the elastic force when the gas film is expanded or compressed. Moreover, the TED from the heat flow is caused by the expansion and compression of the structure. Similar to the SFD, the TED also has two components: the elastic part and the damping part. The fundamental equation of TED based on the Fourier’s law can be presented as a part of the governing equations [37].

4

Y. Cao et al. / Sensors and Actuators A 297 (2019) 111530

2.1. Governing differential equations For the in-plane vibration of a bulk resonator, governing equations of the isotropic thermoelastic solid are given by a modified Fourier’s law of heat transfer [38] and force balance constraints [19]. 2

2

2

∂ U (1 − ) ∂ U (1 + ) ∂ V ∂T (1 − 2 ) + + U¨ − ˛(1 + ) = 2 2 E ∂x2 ∂ y2 ∂x∂y ∂x 2

2

(1)

2

∂ V (1 − ) ∂ V (1 + ) ∂ U ∂T (1 − 2 ) + + V¨ = − ˛(1 + ) 2 2 2 2 E ∂y ∂x ∂x∂y ∂y 2

(

(2)

2

∂ ∂ E ∂U˙ ∂V˙ ˛2 (1 + )E T+ T ) − (C + ( To )T˙ − ˛To + )=0 2 2 1 −  ∂x (1 − 2)(1 − ) ∂x ∂y ∂y

(3)

where U and V are in-plane displacements in the x- and y-directions, T denotes the structural temperature, and , C , , and To are thermal conductivity, heat capacity, density, and initial temperature, respectively. Fig. 2 shows boundary conditions of the microplate. There is no flow of heat across the boundary surfaces of the resonator [39]. The displacements and temperature are expressed by modified Fourier series and written as: U(x, y, t) = u(x, y)ejωt =

 ∞ ∞ 

∞ 

Amn cos am x cos bn y + 1b (y)

m=0 n=0 ∞ 

+ 1a (x)

∞ 

cn cos bn y + 2a (x)

n=0

+ 1a (x)

m=0



bm cos am x

m=0

ejωt

dn cos bn y

(4)

n=0

V (x, y, t) = v(x, y)ejωt = ∞ 

∞ 

am cos am x + 2b (y)

 ∞ ∞ 

∞ 

Bmn cos am x cos bn y + 1b (y)



m=0 n=0 ∞



gn cos bn y + 2a (x)

n=0

hn cos bn y

∞ 

em cos am x + 2b (y)

m=0

fm cos am x

m=0

(5)

ejωt

n=0

T (x, y, t) = (x, y)ejωt =

∞ ∞  

Tmn cos am x cos bn yejωt

(6)

m=0 n=0

where am = m/a and bn = n/b, in which a and b are the length and width of the bulk resonator, respectively. Note that four single Fourier series on the right-side of Eq. (4) is equal to the derivatives of the displacements at corresponding edges, which shifts the possible discontinuity or jump from the x derivative of the displacement function to certain single Fourier series. They guarantee a continuous derivative of the displacement at the edge and improve the convergence of the standard Fourier series [32,40]. Similar conditions are also applied in the Eq. (5).



1a (x) = aς x (ς x − 1)2 , 2a (x) = aς 2x (ς x − 1), ς x = x/a

(7)

1b (x) = bς y (ς y − 1)2 , 2b (x) = bς 2y (ς y − 1), ς y = y/a

To obtain unknown expansion coefficients, all the sine terms and polynomials in Eqs. (1)–(3), which are calculated by substituting Eqs. (4)–(7) into the governing equations, are transformed into cosine series [32], as shown in Appendix. Finally, Eqs. (1)–(3) are expressed as: ( 2am + +( +

1− 2 1− 1− 1− bn )Amn + (ˇ1n 2am − 1n )am + (ˇ2n 2am − 2n )bm + ( ˛1m 2bn − ε1m )cn 2 2 2 2

1− 1 +   1 +  1 +  p q p p Bpq ap bq am bn + ep 1n ap am + fp 2n ap am ˛2m 2bn − ε2m )dn − 2 2 2 2

1 +  2

( 2bn +

p

q

gq 1m bq bn +

q

1 +  2

2

q

q



hq 2m bq bn − ˛(1 + )

q

p

(8)

p

p

Tpn ap am −

p

ω2 cL2

(Amn + am ˇ1n + bm ˇ2n + cn ˛1m + dn ˛2m ) = 0

1− 2 1− 1− 1− am )Bmn + ( ˇ1n 2am − 1n )em + ( ˇ2n 2am − 2n )fm + (˛1m 2bn − ε1m )gn 2 2 2 2

+(˛2m 2bn − +

1 + 

q

1− 1 +   1 +  1 +  p q p p ε2m )hn − Apq ap bq am bn + ap 1n ap am + bp 2n ap am 2 2 2 2 p

q

cq 1m bq bn +

1 +  2

q

q

q

dq 2m bq bn − ˛(1 + )

 q

p

p

q

Tmq bq kbn −

ω2 cL2

(Bmn + em ˇ1n + fm ˇ2n + gn ˛1m + hn ˛2m ) = 0

(9)

Y. Cao et al. / Sensors and Actuators A 297 (2019) 111530



+



jω˛T0 E

Tmn =

˛2 (1 + )ET0

(1 − )[k( 2am + 2bn ) + jω(Cv + q

Bmq bq kbn − 1n em − 2n fm +

×( )]

) (1 − 2)(1 − q gq bq ˛1m kbn +

q

q

p

Apn ap am +

p



5 p

ap ap ˇ1n am +



p

p

bp ap ˇ2n am − 1m cn − 2m dn

p

(10)

q

hq bq ˛2m kbn )

q



where ω denotes the angular frequency and cL = E/ (1 − 2 ) is the longitudinal wave speed. Eq. (10) shows that the expansion coefficient of temperature can be expressed by expansion coefficients of in-plane displacements. And Eq. (10) is rewritten in matrix form. T = T P P+T C C

(11)

Where P = {a0 , ..., aM , b0 , ..., bm , c0 , ..., cN , d0 , ..., dN , e0 , ..., eM , f0 , ..., fM , g0 , ..., gN , h0 , ..., hN }

T

(12)

C = [A00 , ..., Am0 , ..., Amn , ..., AMN , B00 , ..., Bm0 , ..., Bmn , ..., BMN ]T

(13)

By extracting the coefficients of the formulas and substituting Eqs. (11)–(13) into Eqs. (8) and (9), governing equations can be given in matrix form.



A11

A12

A21

A22





C+

B11

B12

B21

B22





P+

T A1 T A2



(T P P+T C C) −

ω2 cL2



E 11

E 12

E 21

E 22





C+

F 11

F 12

F 21

F 22

 P

=0

(14)

2.2. Boundary conditions for the bulk resonator Boundary conditions for the in-plane vibration are restrained by T-shape anchors, squeezing gas pressure, and electric load in Fig. 2, which are depicted by a set of linear springs at each edge. When the stiffness coefficient of spring is zero along certain edge, this edge is free. It is assumed that each edge of the bulk resonator is specified in terms of two independent linear springs (i.e., tangential and normal). For the T-shape anchor, each point support is described in terms of an independent spring and divided into two equal parts in the x- and y-directions. Similar to T-shape anchor, the stiffness and damping coefficients of the SFD and electric load are specified in terms of two independent directions. When different coefficients are selected, different boundary conditions can be determined. In accordance with thermoelastic theory, stress boundary conditions for elastically restrained edges are given as:

∂U ∂V ∂U ∂V + , knx0 U + Fsx0 − Fex0 = + − ˛(1 + )T ∂y ∂x ∂x ∂y

At

x = 0, ktx0 V =

At

x = a, ktx1 V = −

At

y = 0, kty0 U =

At

y = b, kty1 U = −

(15)

∂U ∂V ∂U ∂V − , knx1 U + Fsx1 − Fex1 = − − + ˛(1 + )T ∂y ∂x ∂x ∂y

(16)

∂U ∂V ∂U ∂V + , kny0 V + Fsy0 − Fey0 =  + − ˛(1 + )T ∂y ∂x ∂x ∂y

(17)

∂U ∂V ∂U ∂V − , kny1 V + Fsy1 − Fey1 = − − + ˛(1 + )T ∂y ∂x ∂x ∂y

(18)

Where E, , and ˛ are Young’s modulus, Poisson’s ratio, and the coefficient of thermal expansion, respectively; kn and kt represent the support stiffness coefficients of the T-shape anchor for the normal and tangential directions, respectively; Fs is the equivalent force due to the pressure field of the squeeze gas between electrodes and resonator; and Fe is the electromechanical coupling equivalent force. 2.2.1. Squeeze film damping For small displacement of the bulk resonator, the nonlinear Reynolds equation [41] due to inertial and viscous effect of the fluid can be linearized [42] and applied to solve the pressure of squeeze film between plate side and electrodes. For small variations in local pressure and gap spacing at x = 0, and p = Pa + ıp, the dimensionless Reynolds equation is expressed as [11,13]: 2

2

∂ P ∂ P + = ∂ y2 ∂z 2



∂P ∂G + ∂t ∂t



(19)

where P = ıp/Pa is the normalized local pressure variation. G = U(0, y, t)/g is the normalized local gap variation and.  = 12eff /Pa g 2 . refers to the squeeze number. eff signifies the equivalent viscosity of the fluid [43]. The boundary conditions of squeeze film pressure in Fig. 1 is equal to the ambient pressure at the bottom and top lateral surfaces of the gap region, such that P(y, 0, t) = P(y, H, t) = 0, and the pressure gradient is zero at the side edges: ∂P(0, z, t)/∂y = ∂P(b, z, t)/∂y = 0. The pressure function is obtained by the Green’s function method [13,44]. ıp(y, z, t) = + ˇ2m



 m=0 

 m =0

2

n =odd

bm ) × cos 

2

   −4Pa jωˇl Amm + ˇ1m am ×(  gn  km n / + jω m=0 m =0 m=0

m y b

sin

n z H

ejωt

(20)

where km n = ( mb ) + ( nH ) , H is thickness of the resonator. ˇl is the peak of displacement at certain mode to normalize the displacement function of boundary edge. Eq. (20) shows that pressure ıp can be divided into real and imaginary parts to represent the elastic and

6

Y. Cao et al. / Sensors and Actuators A 297 (2019) 111530

damping effects of SFD, respectively. The equivalent elastic and damping coefficients are expressed by the mean values of the local elastic and damping of the side boundary [11].



kseff =



n =odd gn 2  2

m =0



(

 m=0

n=0



cseff =



(

m=0

n=0



m=0

×(



2 2 2 km  n / + ω



m=0

m =0

m=0



 

Amn cos bn y + 1b (y)



am + 2b (y)

km n /

8Pa n =odd gn 2  2



2 2 2 km  n / + ω

Amn cos bn y + 1b (y)

m =0





ω2

−8Pa

×(

m=0

m =0



am + 2b (y)

m=0

 m=0

am + ˇ2m

 m=0

bm ) cos

m y b



/ (21)

bm )

 m=0

Amm + ˇ1m

Amm + ˇ1m

 m=0

am + ˇ2m



 m=0

bm ) cos

m y / b

bm )

(22)

2.2.2. Electrical behavior The bulk resonator is surrounded by capacitive transducers used as driving and sensing [27,45]. For the in-plane vibration, electrical behavior is considered one-dimensional in the y- or x-direction from the view of a variable gap when the thickness of bulk resonator is small in Fig. 2. The electrostatic driving force generated by a capacitive transducer is given by: Fe =

ε0 A

2 ε0 A 2U(0, y, t) (1 + + ...)(Vdc − Vac ejωt ) g 2g 2

2

2(g − U(0, y, t))2

(Vdc − Vac ejωt ) =

(23)

where ε0 is the permittivity of free space [46]. An electric load is calculated by Eq. (23) in some boundary condition when AC component Vac is applied. The electric load is equal to zero when Vac is not used for any electrodes. 2.2.3. T-shape anchor description In the bulk resonator restrained by four T-shape anchors, the T-shape anchor is considered as an elastic point-support spring that links the fixed position to the resonator. Then, the support stiffness along normal and tangential direction knx0 and ktx0 [33] are expressed as: knx0 (y) =

Ix0 

i knx0 ı(y

i − yx0 ), ktx0 (y)

=

i=1

Ix0 

i i ktx0 ı(y − yx0 )

(24)

i=1

i i where Ix0 is the number of elastic point-support along edge, ı(x) is the Dirac delta function, while knx0 and ktx0 represent the stiffness at i . That variation yi can be any value means any fixed position of the constraint. Moreover, the stiffness in Eq. (24) support located at yx0 x0 is infinitely large, the anchor is clamped. The free anchor exists for case that the stiffness is zero. The remaining boundary conditions Eqs. (15)–(18) are expressed as similar formulas. Subsequently, substituting Eqs. (4)–(7) and (21)–(24) into Eqs. (15)–(18) and ignoring the nonlinear forcing terms [9] will lead to linearized boundary equations. By gathering coefficients of cosine terms, the relationship of the Fourier coefficients can be established as:

ı¯ am am + ı¯ am 1m

N 

gn + ı¯ am 2m

n Iy0 N  



dn

n

N 

hn −

i i i kty0 2a (xy0 ) cos am xy0

n

=

M N   s

ı¯ am bm + ı¯ am 1m

N 

gn (−1)n + ı¯ am 2m

dn (−1)n

n

+ı¯ am

M  N 

+

m

N 

cn (−1)n

Iy1 

n

 M

+ ı¯ am

M N  

s

(25) p Bpn ap am

(m = 0, ..., M)

n

i i i kty1 1a (xy1 ) cos am xy1

i=1

 Iy1

N

i i i kty1 2a (xy1 ) cos am xy1 =−

p

i cos am xy0

p

n

Bpn (−1)n ap am

i cos as xy0

Asn (−1)n

n

i i i kty1 cos as xy1 cos am xy1

(26)

i=1

(m = 0, ..., M)

n

Amt

Ix0 

t

M 

i kty0

i=1

i=1 N M   p

m



Asn



hn (−1)n +

Iy1

N

i i i kty0 1a (xy0 ) cos am xy0

i=1 Iy0

n

N 

n



cn

n

i=1

+

Iy0 N  

bm

i knx0

i cos bt yx0

i cos bn yx0

m

i=1 Ix0 

+

M 

am

Ix0 

i i i knx0 1b (yx0 ) cos bn yx0 + ı¯ bn

m

i=1 M 

i i i knx0 2b (yx0 ) cos bn yx0 + (kseffx0 + iwc seffx0 − kex0 )ı¯ bn (

Amn + ˇ1n

M 

m

i=1

= ı¯ bn cn + ı¯ bn 1n

M  m

em + ı¯ bn 2n

M  m

M N  

fm − ˛(1 + )ı¯ bn

M  m

Tmn

(n = 0, ..., N)

m

q

Bmq bq kbn

q

am + ˇ2n

M  m

bm )

(27)

Y. Cao et al. / Sensors and Actuators A 297 (2019) 111530 M N   m

Ix1 

Amt (−1)m

t

M 

+

i i i knx1 cos bt yx1 cos bn yx1 +

bm (−1)m

m

am (−1)m

Ix1 

m

i=1 Ix1 

M 

i i i knx1 1b (yx1 ) cos bn yx1

i=1 M 

i i i knx1 2b (yx1 ) cos bn yx1 + (kseffx1 + iwc seffx1 − kex1 )ı¯ bn (

+ˇ2n

m

bm (−1) ) = −ı¯ bn dn + ı¯ bn

M N  

m

m

+˛(1 + )ı¯ bn

M 

Amn (−1)m + ˇ1n

m

i=1 M 

7

Tmn (−1)m

Bmq (−1)

m

q bq kbn

− ı¯ bn 1n

M 

q

M 

am (−1)m

m m

em (−1) − ı¯ bn 2n

M 

m

(28) fm (−1)

m

m

(n = 0, ..., N)

m M  N  s

Bsn

Iy0 

n

i i i kny0 cos as xy0 cos am xy0 +

i i i kny0 1a (xy0 ) cos am xy0 + ı¯ am

+

hn

n

p

 N

i i i kny0 2a (xy0 ) cos am xy0 + (kseffy0 + iwc seffy0 − key0 )ı¯ am ( N 

cn + ı¯ am 2m

N 

n M N  

Bsn (−1)n

Iy1 

n

N 

+

dn −˛(1 + )ı¯ am

n

hn (−1)n

n

Bmn + ˛1m

+˛2m

i i i kny1 cos as xy1 cos am xy1 +

N 

gn (−1)n

n

Iy1 

p N 

(29)

i=1

Bmn (−1)n + ˛1m

n

n

hn )

n

i i i kny1 1a (xy1 ) cos am xy1 N 

N M  

N 

(m = 0, ..., M)

Tmn

i i i kny1 2a (xy1 ) cos am xy1 + (kseffy1 + iwc seffy1 − key1 )ı¯ am (

hn (−1)n ) = ı¯ am

+˛(1 + )ı¯ am

gn + ˛2m

n

i=1 N 

n

n

i=1 Iy1 

N 



p

Apn ap am

N

n

i=1

= ı¯ am em + ı¯ am 1m

M N  

i=1

Iy0

N

s

gn

n

i=1

 

Iy0 N  

p Apn (−1)n ap am − ı¯ am 1m

N 

n

N 

gn (−1)n

n

cn (−1)n − ı¯ am 2m

n

N 

(30)

dn (−1)n − ı¯ am fm

n

(m = 0, ..., M)

Tmn

n

ı¯ bn 1n

M 

am + ı¯ bn 2n

M 

m M N   m M 

+

m

q Amq bq kbn

+

q

am (−1)m + ı¯ bn 2n

M 

m

Ix1 

Bmt

t



m

i ktx0

bm (−1)m + ı¯ bn hn +

M 

N M   m

Bmt (−1)

t

i cos bt yx0

i cos bn yx0

em (−1)m

m

i=1 Ix1 

fm

i i i ktx0 2b (yx0 ) cos bn yx0

i=1

(31)

(n = 0, ..., N)

i=1

i i i ktx1 2b (yx1 ) cos bn yx1 = ı¯ bn

m



Ix0 M   m

Ix0 

m

fm (−1)m

M N  

i i i ktx0 1b (yx0 ) cos bn yx0

i=1

M  N  m

m M 

em

m

+ı¯ bn gn = ı¯ bn

ı¯ bn 1n

bm −

Ix0 M  

i i i ktx1 cos bt yx1 cos bn yx1

Ix1 

i i i ktx1 1b (yx1 ) cos bn yx1

i=1 q

Amq (−1)m bq kbn

(32)

q

(n = 0, ..., N)

i=1









Where ı¯ am = a 1 + ım0 /2, ı¯ bn = b 1 + ın0 /2, and ım0 is the Kronecker delta. M and N refer to the truncated number of all of the series coefficients. Substituting Eq. (11) into Eqs. (15)–(18), equations of boundary conditions are rewritten into matrices form. HP = QC+T Q (T P P+T C C)

(33)

Where

 Q =

qT11

qT21

qT31

qT41

qT51

qT61

qT71

qT81

qT12

qT22

qT32

qT42

qT52

qT62

qT72

qT82

T (34)

8

Y. Cao et al. / Sensors and Actuators A 297 (2019) 111530



a11

⎢ ⎢ a21 ⎢ ⎢a ⎢ 31 ⎢ ⎢ a41 ⎢ H=⎢ ⎢ a51 ⎢ ⎢ ⎢ a61 ⎢ ⎢a ⎣ 71 a81



b12

c13

d14

e15

f16

g17

h18

b22

c23

d24

e25

f26

g27

h28 ⎥

b32

c33

d34

e35

f36

g37

b42

c43

d44

e45

f46

g47

b52

c53

d54

e55

f56

g57

b62

c63

d64

e65

f66

g67

b72

c73

d74

e75

f76

g77

⎥ ⎥ ⎥ h38 ⎥ ⎥ h48 ⎥ ⎥ ⎥ h58 ⎥ ⎥ ⎥ h68 ⎥ ⎥ ⎥ h78 ⎦

b82

c83

d84

e85

f86

g87

h88

(35)

The coefficient matrix P can be eliminated by substituting Eq. (33) into Eq. (14). The final matrix form of governing equations can be rewritten as:



K−

ω2 cL2



M

C=0

(36)

Where K = A + B(H − T A T P )−1 (Q +T Q T C ) + T A (T P (H − T A T P )−1 (Q +T Q T C ) + T C )

(37)

M = E + F(H − T A T P )−1 (Q +T Q T C )

(38)

Eq. (36) represents a linearly distributed parameter model for the resonator under the coupled effect of SFD, TED, and linearized electrostatic forces. The complex eigenvalue can be determined from Eq. (36), while the physical mode shapes, pressure distribution, and temperature distribution can be determined by substituting the expansion Fourier coefficients of the eigenvector into Eqs. (4)–(6) and Eq. (20). 2.3. Quality factor The quality factor is defined as the ratio of stored energy and energy dissipation per cycle, which is an index for the resolution and phase noise of resonators. The stored energy can be calculated by: i Ener = mi v2i /2

(39)

where m is the equivalent mass and v is the velocity distribution. The maximum stored energy is acquired from the kinetic energy summation of all the mass elements, which can be obtained by dividing the structure into finite elements. The energy dissipation of SFD is caused by the damping-like behavior of pressure, whereas the spring-like behavior of pressure can only transmit energy. Energy dissipation per period, which is caused by the damping effect of SFD for an element, can be determined by multiplying the total force by its velocity.

 Esi

T

˙ Im(ıpi )dAi U(0, yi , t)dt

= 0

(40)

= Im(ıpi )dAi u(0, yi ) Finally, total dissipated energy is solved by overlaying the energy dissipation of all elements. The thermal strain is not in phase with stress and generates energy dissipation at one cycle [21,37]. The energy dissipation of TED per cycle is given:

 Eti = 0

T

2kTi ∇ 2 Ti dAi Hdt To

2k i ∇ 2 i = dAi H ωTo

(41)

Having the stored and dissipation energies per cycle, the quality factor of the SFD and TED and total quality factor can be defined as following formulas [24]. Q −1 =

Es Et + 2Ener 2Ener

(42)

Fig. 3 intuitively shows the coupling process of the boundary and governing equations, including SFD and TED. Firstly, two matrices of governing equations and boundary conditions can be derived by collecting coefficients of cosine terms obtained from the improved Fourier series method. Substituting boundary conditions into governing equations will determine eigenvalues and the expansion Fourier coefficients. Then, the total equivalent elastic and damping coefficients of squeeze film and temperature are substituted into two matrices to participate in the calculation process, which include the following four factors: equivalent squeeze film elastic and damping coefficients (their initial value is zero), equivalent electrostatic stiffness, equivalent support stiffness, and TED expressions. Finally, the iteration process cycles until the complex frequency reaches the desired precision and the Q-factor is calculated by Eq. (42).

Y. Cao et al. / Sensors and Actuators A 297 (2019) 111530

9

Fig. 3. Calculation flowchart for bulk resonator under the effects of SFD and TED. Table 1 Main parameters of bulk resonator [24]. Parameters

Values

Square side length and width Resonator thickness Initial gap Young’s modulus Poisson’s ratio Air density Ambient temperature Thermal expansion coefficient Thermal conductivity Special heat Permittivity of free space

500m 25m 2m 169 Gpa 0.28 2330 kg/m3 300 K 2.6·10−6 K-1 130 W/(m K) 1.631·106 J/(m3 K) 8.854·10−12 F/m

3. Results and discussions After several iteration processes in Fig. 3, the pressure profile inside the gap is stabilized and determined by known material and geometrical parameters in Table 1. As shown in Fig. 4, the fluctuation of gas pressure corresponds to the structural deformation of five modes. Local pressure decreases when the structural compression increases the gap, but local pressure increases due to structural expansion. For the extension mode at some point, all of the pressure distributions at different boundaries of bulk resonator are negative because the whole structure is compressed at the edge x = 0. For lame modes, the sinusoidal variation of the pressure is caused by the sinusoidal variation of the bulk structure. Moreover, profiles in structural deformation and pressure distributions become complex for high-order frequencies. Note that anchors in fixed positions are almost not compressed or stretched because they are in the node positions of the mode shape which generates almost no fluctuation of pressure for the lame modes at boundaries along z direction. However, the opposite deformation trend of adjacent boundaries may result in a stress concentration at a fixed location between anchors and bulk resonator. 3.1. Model verification To verify the numerical model, predicted results are compared with the published literature. Strictly speaking, none literature on bulk resonator is based on the fully coupling model. However, references [9,24,47] gave a conclusion that the damping of the resonator is dominated by mechanical damping of material when the ambient pressure near vacuum where the influence of pressure can be negligible. Thus, the measured quality factor in reference [24] in the case of high vacuum (i.e., 0.01 Torr) is regarded as the quality factors of mechanical damping. Although this quality factor includes the effect of TED, the influence of TED is very small which is verified by COMSOL in section

10

Y. Cao et al. / Sensors and Actuators A 297 (2019) 111530

Fig. 4. Predicted pressure profiles at x = 0 for (a) extension mode shape and (b–e) first four lame mode shapes.

Y. Cao et al. / Sensors and Actuators A 297 (2019) 111530

Fig. 5. Comparison of calculated and experimental results [24] in quality factor with pressure.

11

(a) Imaginary parts of structural frequencies(b) Quality factor of SFD.

Table 2 Comparative frequencies (MHz) of different modes.

fmea [24] fcal [24] a

fcal a

1 st lame

2nd lame

Extension

3rd lame

4th lame

8.27 7.24 7.53 7.53

16.54 14.47 15.05 15.05

8.73 9.42 9.70 9.36

22.58 22.58

30.11 30.11

represents simulated results from COMSOL and the fcal represents calculated results of the reference [24].

Table 3 Total quality factor for the resonator with different M = N under 760 Torr. M=N

Extension mode

1 st lame mode

2nd lame mode

14 18 22 26 30

5.83·104 5.81·104 5.79·104 5.79·104 5.78·104

3.24·104 3.23·104 3.22·104 3.21·104 3.21·104

17.22·104 17.23·104 17.33·104 17.33·104 17.32·104

3.3 of this paper. Therefore, the total quality factor can be legitimately depicted by the relation between the above measured quality factor at high vacuum and quality factor computed by the presented model at different pressures. Fig. 5 shows comparative results for three modes of the bulk resonator where the total quality factor decreases slightly with low pressure (pressure<1 Torr) and then decreases rapidly with high pressure. The reason is that the damping of pressure is proportional to ambient pressure in Eq. (22). The higher ambient pressure contributes to more energy dissipation caused by the effective damping of squeeze film pressure. Furthermore, the quality factor of extension mode is smaller than that of lame mode. This can be explained by mechanical damping of structural mode shapes in Fig. 4. The predicted results agree well with experimental results for three modes. The trivial deviation may be ascribed to a gap 15 ␮m between the electrodes which is not presented in the theoretical model. And the effect of the slit on pressure profile is also not reflected. For the resonator, the fact that the stored energy for the extensional frequency is larger than the dissipation energy of squeezed film damping generates a slight discrepancy between the predicted and measured quality factors. Table 2 summarizes the measured and calculated results for the resonant frequencies. Predicted lame frequencies are slightly less than the measured frequencies and agree well with structural natural frequencies by FE method. The above comparison results show that the proposed model is correct and can be used for parameter analysis. In the current calculations, displacement and temperature expressions retain the same number of interception and contain the first (M + 1)×(N + 1) terms. The convergence of the quality factor is examined. Table 3 shows total quality factor of system which is determined in the same way as that in Fig. 5 for three modes. The expected convergence behavior is observed from a sufficiently accurate solution when only a slight fluctuation is observed between adjacent results. According to predicted solutions, the truncated number M = N = 26 is applied in all the remaining calculations. 3.2. Squeeze film damping analysis 3.2.1. Effect of pressure on QSFD Fig. 6(a) shows imaginary parts of structural frequencies according to pressure while thermal damping is ignored by setting the temperature to zero. These curves shift to high value with the increase of ambient pressure, which indicate that the higher ambient pressure improves the damping force of squeeze gas [9]. The reason is that high pressure has more molecules colliding with each other than low pressure to dissipate energy of system. Fig. 6(b) shows the decrease of QSFD when ambient pressure increases, and this decrease is rapid when the ambient pressure is Pa<1 Torr. Then, the decline rate of QSFD is somewhat moderate as the pressure increases at 1 < Pa<760 Torr. Hence, high vacuum seal can significantly reduce energy consumption and improve the quality factor. In addition, QSFD increases with the

12

Y. Cao et al. / Sensors and Actuators A 297 (2019) 111530

Fig. 6. Parametric plot of (a) imaginary parts of structural frequencies and (b) quality factor of SFD for the extension and lame frequencies with respect to pressure. (a) Imaginary parts of structural frequencies(b) Quality factor of SFD. Table 4 Gas viscosities and mean free paths [46] at 300 K, 1 atm, and equivalent viscosities. Gas property

H2

He

N2

Air

O2

Ar

CO2

Viscosity, Pa·s Mean free path, nm Equivalent viscosity, Pa·s

8.89·10−6 126 1.74·10−6

20.33·10−6 198 2.56·10−6

17.78·10−6 67.7 5.93·10−6

18.71·10−6 68.7 6.17·10−6

20.61·10−6 72.7 6.50·10−6

22.72·10−6 72.7 7.17·10−6

15.13·10−6 45.4 6.71·10−6

Table 5 Dependence of QSFD on different gas viscosities for the extension and lame modes. Gas Extension 1 st lame 2nd lame 3rd lame 4th lame

H2

He 5

4.19·10 2.36·105 1.22·106 3.45·106 7.16·106

N2 5

4.94·10 2.63·105 1.54·106 4.33·106 8.86·106

Air 5

8.1·10 4.21·105 2.46·106 6.67·106 1.351·107

O2 5

8.29·10 4.31·105 2.51·106 6.8·106 1.378·107

Ar 5

8.53·10 4.45·105 2.57·106 6.98·106 1.415·107

CO2 5

8.99·10 4.71·105 2.7·106 7.33·106 1.486·107

8.67·105 4.53·105 2.61·106 7.09·106 1.437·107

increase in resonant frequency. The reason is that a part of fluid, which is squeezed out of the control region at fundamental frequency, is sucked back into the system to store energy when the structure vibrates at high frequencies. This is also the purpose of the development and advantages of bulk resonator compared with plate resonator in transverse vibration. 3.2.2. Effect of different gas media on QSFD Generally, the resonator operates in a certain gas medium. Gas viscosity is an important parameter of squeeze gas, which affects the performance of resonators. Here, hydrogen, helium, nitrogen, air, oxygen, argon, and carbon dioxide are considered to study the effect of gas media on the bulk resonator. In the analytical calculation, the ambient temperature and pressure are 300 K and 100 Torr, respectively. Table 4 shows the different gas properties and equivalent viscosities. Table 5 summarizes the QSFD of different gas media for the extension and first four lame modes. The equivalent viscosity of hydrogen leads to the minimum QSFD value, and argon obtains the best QSFD value which is twice as big as that of hydrogen. In conjunction with Table 4, a conclusion emerges that the gas with high equivalent viscosity will result in small energy dissipation and high-quality factor. This case directs us that argon with high equivalent viscosity is a good choice as operating medium of resonators. 3.2.3. Effect of the gap between the resonator and fixed electrodes on QSFD Several studies on the gap between electrodes and devices regarded that the coupling characteristics of resonators were related to the gap, and many methods have been developed to reduce the gap [48,49]. However, gap reduction can increase the equivalent damping of pressure. Fig. 7(a) shows that imaginary parts of structural frequencies decrease as the gap increases at one atmospheric pressure. Moreover, imaginary parts of different frequencies tend to be the same when the gap is greater than 6 ␮m. It means that the damping of SFD is affected by mode shapes only in small clearance which is less than 6 ␮m. In Fig. 7(b), the QSFD increases with the increase in gap dimensions. The QSFD also increases with the increase of resonant frequency. The reason is that only a few gas as a damping effect is squeezed out of the plate region. Thus, a higher gap size should be selected to increase the quality factor and reduce loss under the condition of satisfying the operating electrostatic force. Moreover, the large clearance is convenient for the fabricating process. 3.2.4. Effect of structural thickness on QSFD Fig. 8 explores the relationship between structural thickness and QSFD at Pa = 10 Torr. The quality factor first decreases to a minimum and then slightly rises as the plate thickness increases. For a thin bulk resonator, the fact that there is less gas in the space between resonator and electrodes reduces the damping of the squeezed gas. This condition leads to a corresponding reduction in the dissipated energy that is

Y. Cao et al. / Sensors and Actuators A 297 (2019) 111530

13

Fig. 7. (a) Imaginary parts of structural frequencies and (b) quality factor of SFD for the extension and lame frequencies as a function of the gap.

Fig. 8. Quality factor for lame and extension frequencies with structural thickness.

removed by the gas. As the thickness increases, a critical value will be reached so that the gas inside the space is only squeezed and plays an elastic role, and the gas at the boundary is squeezed out or sucked back and damps the system. The slight rise of QSFD may be due to the compressibility of the gas and the increase of structural stored energy. Additionally, the critical thickness decreases with the increase in resonant frequencies. 3.2.5. Effect of the structural length on QSFD Fig. 9 shows the effect of the various structural length on QSFD for the first lame frequency at Pa = 10 Torr. The result indicates that the bulk resonator with a small structural length can achieve a high QSFD . Furthermore, the change of QSFD decreases with the increase in resonator length. As the length of the structure increased, the contact area of squeeze gas also increased, such that additional gas is squeezed out of the plate domain. Additional squeeze gas flowed out to remove the energy when the frequency of the large-sized structure is reduced. 3.3. Thermoelastic damping analysis To verify the predicted model, in particular, numerical results of TED of bulk resonator for different resonant frequencies compare with that solved by COMSOL with main parameters in Table 1. The boundary conditions are shown in Fig. 2 and fluid forces in this analysis are not considered. The analysis module of the model in COMSOL is “thermal perturbation, eigenfrequency” where predefined mesh element size is “Extra Fine” including multiple preset parameters to build mesh element quality and size automatically: maximum element size is 19.4 ␮m, minimum element size 0.832 ␮m, maximum element growth rate 1.35, curvature factor 0.3, and resolution of narrow regions 0.85. The detailed process for calculating the quality factor of TED can be found in tutorial and demo app files of COMSOL official website. Note that the minimum size of element is smaller than the minimum size of the resonator and the rest is meshed automatically by the COMSOL, which ensure that the overall model has good mesh quality. Table 6 indicates that numerical results agree well with simulated results. The slight distinction is due to the existence of the energy dissipation of four anchors which is not taken into account in the theoretical model. Interestingly, quality factor of TED at extension frequency has minimum value although its frequency is not minimal. The fact that the whole structure of bulk resonator is stretched or compressed, which is different from the structure at lame frequencies, generates the same temperature trend in Fig. 10(a). This prevents the resulting temperature from being reabsorbed into the system, increasing the

14

Y. Cao et al. / Sensors and Actuators A 297 (2019) 111530

Fig. 9. Quality factor for the first lame frequency as a function of structural length.

(a) Extension frequency. (b) 1st lame frequency. (c) 2nd lame frequency.

Table 6 Comparative results of QSFD for different structural frequencies of bulk resonator. Quality factor

1 st lame

2nd lame

Extension

3rd lame

4th lame

COMSOL Predicted

6.08·1010 6.37·1010

2.70·1010 2.57·1010

2.85·107 2.95·107

1.98·1010 2.04·1010

13.2·109 9.76·109

Fig. 10. Temperature distribution of bulk resonator for (a) the extension and (b–c) first two lame frequencies.

Y. Cao et al. / Sensors and Actuators A 297 (2019) 111530

15

Fig. 11. QTED of different resonant modes as a function of ambient temperature.

Fig. 12. Parametric plot of QTED at first lame frequency with respect to structural length.

(a) 1st lame frequency. (b) 2nd lame frequency. (c) Extension frequency.

loss of energy. Furthermore, quality factor of TED decreases as the increase of lame frequencies. This fact can be explained by that the high-frequency vibration of the structure generates higher heat than that at low frequency in Fig. 10(b–c), which elevates the temperature of the structure close to each other and dissipates more energy into the environment. The slight fluctuation of temperature distribution at boundary may be related to the vibration wavelength of the structure. 3.3.1. Effect of ambient temperature on QTED Fig. 11 shows the variation in the QTED of the resonator with respect to ambient temperature. Findings show that the QTED decreases with the increase in ambient temperature, and that the dissipation energy from the TED increases with the increase in temperature. This change can be attributed to the thermal energy of the molecule of the structure and the number of intermolecular collisions that increase with the increase in ambient temperature to cause additional energy loss. For the extension vibration, the stretching of the resonator does not occur simultaneously with the compression of the resonator, but occurs sequentially. Therefore, all the energy produced by the vibration of the structure is dissipated. For the lame vibration, the stretching of the resonator exists simultaneously with the compression. The temperature of the whole structure is not consistent so that the energy from the local high temperature will flow to the position of local low temperature. 3.3.2. Effect of the structural length on QTED Generally, structural length as an important factor could change the vibration frequency of thin plate structure [9]. To study the effect of this factor on performance of resonator, comparison of different structural lengths at first lame frequency in Fig. 12 indicates that the higher QTED favors a large-sized resonator. The reason is that the in-plane vibration frequency of a large-sized resonator is smaller than that of a small-sized structure, and the large-sized structure produces lesser dissipation energy than the small-sized structure. Subsequently, a relatively large-sized structure should be selected to achieve a high-quality factor. 3.4. Fully coupled system analysis In a final part of study, the combined effect of two mechanisms of SFD and TED is performed to explore the way they interact each other. The mechanical damping is ignored. Fig. 13 shows the real and imaginary parts of the first two lame and extension frequencies with respect to pressure values which are selected by log axis of Fig. 14. The x axis is the change in real part of frequency. These curves show that the pressure in low pressure regime Pa<1 Torr give the small imaginary part of frequency, while the high pressure 1 < Pa<760 Torr

16

Y. Cao et al. / Sensors and Actuators A 297 (2019) 111530

Fig. 13. Real and imaginary parts of (a–b) the first two lame and (c) extension frequencies with respect to pressure.

Fig. 14. Relative contribution of SFD as a function of pressure (300 K).

influences the vibrational behavior of system: real and imaginary parts of frequencies. Fig. 14 shows that the relative contribution (in percentage relative fraction) of SFD with respect to the pressure, which is refined as the ratio of the total quality factor of SFD and TED which is estimated by Eq. (42) and the quality factor of SFD. Whereas, the remaining part that the distance from the curves in Fig. 14 to 100% at the top of the curve corresponds to the relative importance of TED over the total quality. Findings suggest that the TED is significant in the low pressure regime Pa<1 Torr. When ambient pressure is larger than 1 Torr, the contribution of TED is almost equal to zero and the proportion of SFD is equal to 1. These indicate that the effect of SFD on the quality factor of the resonator is absolute dominant in the high pressure regime Pa>1 Torr, which is ascribed to the proportional relation between the pressure and the equivalent damping of gas film as shown in Eq. (22). However, this conclusion does not apply to TED for the extension mode, which still affects the quality factor of

Y. Cao et al. / Sensors and Actuators A 297 (2019) 111530

17

resonator inside large pressure regime. Moreover, this effect on system diminishes because the SFD becomes dominate as the pressure increases. Furthermore, the relative contribution of SFD on the total quality factor gradually reduces with the increase of the frequencies of lame modes of the resonator. 4. Conclusions A general method was developed for the performance of bulk resonators in the presence of SFD and TED. The bulk resonator was described as an in-plane vibration plate where the governing relations were derived by the modified Fourier series method. The pressure distribution of squeeze gas in the gap between resonator and electrodes was calculated with the Reynolds equation through the Green’s function method, and the effective damping and spring coefficients were deduced. The anchor, which linked the fixed position and resonator, and the linearized electrostatic force were simplified as equivalent but separate springs. The predicted results were verified by the published results and simulation results of finite element. The study of SFD included the effects of various structural dimensions, gas media and ambient temperature on bulk resonator. The results showed that ambient pressure is a significant factor to influence the quality factor of SFD. Good vacuum sealing can better promote QSFD by almost five orders of magnitude than the atmospheric environment. Argon gas is a suitable choice when the device is not vacuum sealed. The quality factor of SFD for bulk resonator increases as the resonant frequencies rises. However, the quality factor of TED gives an opposite trend with the increase of structural frequency. Interestingly, the bulk resonator vibrated at the extension mode generates a minimum quality factor of TED. Acknowledgments The authors acknowledge the financial support of the National Natural Science Foundation of China (Nos. 51575170 and 51875185), the Foundation of Hunan Province (2018JJ1006), and the research program supported by the Science and Technology Commission of Shanghai Municipality (17DZ1201000). Appendix A. Supplementary data Supplementary material related to this article can be found, in the online version, at doi:https://doi.org/10.1016/j.sna.2019.111530. References [1] H. Luo, G. Zhang, L.R. Carley, G.K. Fedder, A post-CMOS micromachined lateral accelerometer, J. Microelectromech. Syst. 11 (2002) 188–195. [2] J.T.M. Van Beek, R. Puers, A review of MEMS oscillators for frequency reference and timing applications, J. Micromech. Microeng. 22 (2011), 013001. [3] A.H. Nayfeh, M.I. Younis, A new approach to the modeling and simulation of flexible microstructures under the effect of squeeze-film damping, J. Micromech. Microeng. 14 (2004) 170. [4] J. Yang, T. Ono, M. Esashi, Energy dissipation in submicrometer thick single-crystal silicon cantilevers, J. Microelectromech. Syst. 11 (2002) 775–783. [5] Y. Chen, B. Zhang, N. Zhang, M. Zheng, A condensation method for the dynamic analysis of vertical vehicle–track interaction considering vehicle flexibility, J. Vib. Acoust. 137 (2015), 041010. [6] S. Pourkamali, Z. Hao, F. Ayazi, VHF single crystal silicon capacitive elliptic bulk-mode disk resonators-part II: implementation and characterization, J. Microelectromech. Syst. 13 (2004) 1054–1062. [7] Y. Chen, B. Zhang, S. Chen, Model reduction technique tailored to the dynamic analysis of a beam structure under a moving load, Shock. Vib. 2014 (2014). [8] Y. Chen, J. Dagny, A. Peter, Underwater dynamic response at limited points expanded to full-field strain response, J. Vib. Acoust. 140 (2018), 051016. [9] P. Belardinelli, M. Brocchini, L. Demeio, S. Lenci, Dynamical characteristics of an electrically actuated microbeam under the effects of squeeze-film and thermoelastic damping, Int. J. Eng. Sci. 69 (2013) 16–32. [10] K. Feng, H.Q. Guan, Z.L. Zhao, T.Y. Liu, Active bump-type foil bearing with controllable mechanical preloads, Tribol. Int. (2017). [11] M. Bao, H. Yang, Squeeze film air damping in MEMS, Sens. Actuators A Phys. 136 (2007) 3–27. [12] Z. Guo, K. Feng, T. Liu, P. Lyu, T. Zhang, Nonlinear dynamic analysis of rigid rotor supported by gas foil bearings: effects of gas film and foil structure on subsynchronous vibrations, Mech. Syst. Signal. Process. 107 (2018) 549–566. [13] R.B. Darling, C. Hivick, J. Xu, Compact analytical modeling of squeeze film damping with arbitrary venting conditions using a Green’s function approach, Sens. Actuators A Phys. 70 (1998) 32–41. [14] A.K. Pandey, R. Pratap, Effect of flexural modes on squeeze film damping in MEMS cantilever resonators, J. Micromech. Microeng. 17 (2007) 2475–2484 (2410). [15] P. Li, Y. Fang, H. Wu, A numerical molecular dynamics approach for squeeze-film damping of perforated MEMS structures in the free molecular regime, Microfluid. Nanofluid. 17 (2014) 759–772. [16] C. Zener, Internal friction in solids. I. Theory of internal friction in reeds, Phys. Rev. 52 (1937) 230. [17] A.H. Nayfeh, M.I. Younis, Modeling and simulations of thermoelastic damping in microplates, J. Micromech. Microeng. 14 (2004) 1711. [18] C. Li, S. Gao, S. Niu, H. Liu, Thermoelastic coupling effect analysis for gyroscope resonator from longitudinal and flexural vibrations, Microsyst. Technol. 22 (2016) 1029–1042. [19] F.X. Zhou, S.R. Li, Y.M. Lai, Three-dimensional analysis for transient coupled thermoelastic response of a functionally graded rectangular plate, J. Sound Vib. 330 (2011) 3990–4001. [20] S. Prabhakar, S. Vengallatore, Theory of thermoelastic damping in micromechanical resonators with two-dimensional heat conduction, J. Microelectromech. Syst. 17 (2008) 494–502. [21] A. Duwel, R.N. Candler, T.W. Kenny, M. Varghese, Engineering MEMS resonators with low thermoelastic damping, J. Microelectromech. Syst. 15 (2006) 1437–1445. [22] P. Jiang, P. Li, A model for thermoelastic dissipation of the fully clamped rectangular microplates, in: International Conference on Mechanics and Mechatronics, 2015, pp. 237–244. [23] X. Guo, Y.B. Yi, S. Pourkamali, A finite element analysis of thermoelastic damping in vented MEMS beam resonators, Int. J. Mech. Sci. 74 (2013) 73–82. [24] M.S. Hajhashemi, A. Rasouli, B. Bahreyni, Performance optimization of high order RF microresonators in the presence of squeezed film damping, Sens. Actuators A Phys. 216 (2014) 266–276. [25] P. Belardinelli, S. Lenci, L. Demeio, Vibration frequency analysis of an electrically-actuated microbeam resonator accounting for thermoelastic coupling effects, Int. J. Dyn. Control 3 (2015) 157–172.

18

Y. Cao et al. / Sensors and Actuators A 297 (2019) 111530

[26] T. Mattila, J. Kiihamäki, T. Lamminmäki, O. Jaakkola, P. Rantakari, A. Oja, H. Seppä, H. Kattelus, I. Tittonen, A 12 MHz micromechanical bulk acoustic mode oscillator, Sens. Actuators A Phys. 101 (2002) 1–9. [27] Z. Hao, S. Pourkamali, F. Ayazi, VHF single-crystal silicon elliptic bulk-mode capacitive disk resonators-part I: design and modeling, J. Microelectromech. Syst. 13 (2004) 1043–1053. [28] D. Grogg, H.C. Tekin, N.D. Ciressan-Badila, D. Tsamados, M. Mazza, A.M. Ionescu, Bulk lateral MEM resonator on thin SOI with high Q -factor, J. Microelectromech. Syst. 18 (2009) 466–479. [29] J.E.-Y. Lee, A.A. Seshia, Binary excitation of a high-Q bulk acoustic microresonator by actuation polarity inversion, in: ASME 2008 International Mechanical Engineering Congress and Exposition, American Society of Mechanical Engineers, 2008, pp. 443–448. [30] J. Ey Lee, Y. Zhu, A.A. Seshia, A bulk acoustic mode single-crystal silicon microresonator with a high-quality factor, J. Micromech. Microeng. 18 (2008), 064001. [31] D.J. Gorman, Free in-plane vibration analysis of rectangular plates with elastic support normal to the boundaries, J. Sound Vib. 285 (2005) 941–966. [32] J. Du, W.L. Li, G. Jin, T. Yang, Z. Liu, An analytical method for the in-plane vibration analysis of rectangular plates with elastically restrained edges, J. Sound Vib. 306 (2007) 908–927. [33] J. Du, Z. Liu, W.L. Li, X. Zhang, W. Li, Free in-plane vibration analysis of rectangular plates with elastically point-supported edges, J. Vib. Acoust. 132 (2010) 1015–1016. [34] B. Liu, Y. Xing, Comprehensive exact solutions for free in-plane vibrations of orthotropic rectangular plates, Eur. J. Mech. – A/Solids 30 (2011) 383–395. [35] M. Palaniapan, L. Khine, Nonlinear behavior of SOI free-free micromechanical beam resonator, Sens. Actuators A Phys. 142 (2008) 203–210. [36] L. Shao, Nonlinear Vibration of Micromechanical Resonators, Dissertation, National University of Singapore, 2008. [37] R.B. Hetnarski, M.R. Eslami, Thermal Stresses – Advanced Theory and Applications, Springer, Netherlands, 2009. [38] A. Duwel, R.N. Candler, T.W. Kenny, M. Varghese, Engineering MEMS resonators with low thermoelastic damping, J. Microelectromech. Syst. 15 (2006) 1437–1445. [39] J.N. Sharma, D. Grover, Thermoelastic vibration analysis of Mems/Nems plate resonators with voids, Acta Mech. 223 (2012) 167–187. [40] W. Li, X. Zhang, J. Du, Z. Liu, An exact series solution for the transverse vibration of rectangular plates with general elastic boundary supports, J. Sound Vib. 321 (2009) 254–269. [41] K. Feng, Y. Cao, K. Yu, H. Guan, Y. Wu, Z. Guo, Characterization of a controllable stiffness foil bearing with shape memory alloy springs, Tribol. Int. 136 (2019) 360–371. [42] J.J. Blech, On isothermal squeeze films, J. Lubrication Technol. 105 (1983) 615. [43] T. Veijola, H. Kuisma, J. Lahdenperä, T. Ryhänen, Equivalent-circuit model of the squeezed gas film in a silicon accelerometer, Sensor Actuators A: Phys. 48 (1995) 239–248. [44] P. Morse, Herman Feshbach, Methods of Theoretical Physics, McGraw-Hill, 1999. [45] P. Belardinelli, S. Lenci, M. Brocchini, Modeling and analysis of an electrically actuated microbeam based on nonclassical beam theory, J. Comput. Nonlinear Dyn. 9 (2014), 031016. [46] C. Kavitha, M.G. Madhan, Study of squeeze film damping characteristics under different gas mediums in a capacitive MEMS accelerometer, J. Braz. Soc. Mech. Sci. Eng. 38 (2016) 241–252. [47] W.E. Newell, Miniaturization of tuning forks, Science 161 (1968) 1320–1326. [48] J.R. Clark, W.T. Hsu, M.A. Abdelmoneum, T.C. Nguyen, High-Q UHF micromechanical radial-contour mode disk resonators, J. Microelectromech. Syst. 14 (2005) 1298–1310. [49] B. Bahreyni, Fabrication and Design of Resonant Microdevices, William Andrew, 2008.

Biographies

Yuanlong Cao is currently pursuing the Ph.D. degree in College of Mechanical and Vehicle Engineering, Hunan University, China. His research interests include application and nonlinear dynamics of microresonators, electromagnetic filters, and antennas.

Ke Yu received the Bachelor of University of Electronic Science and Technology of China, 2016. He is currently pursuing the M.S. degree in College of Mechanical and Vehicle Engineering, Hunan University, China. His researches include design and manufacture of active gas bearing, signal analysis and processing, and experimental study on oil-free ACM.

Kai Feng is a professor and doctoral supervisor in Hunan University, China. He received the Ph.D. degree in The University of Tokyo, in 2009. He was a researcher at the University of Tokyo, Japan, in 2009–2010, and he is an assistant professor in Tokyo Institute of Technology, where he conducted research on high speed and precision gas lubrication technology for air foil bearings. He is an associate editor of Tribology Transactions. His research activities have focused on the high speed and precision gas lubrication technology; active and intelligent gas bearing technology; super high speed motor and control of low gravity air flotation system.

Y. Cao et al. / Sensors and Actuators A 297 (2019) 111530

19

Hanqing Guan is pursuing the Ph.D. degree in College of Mechanical and Vehicle Engineering, Hunan University, China. His researches focus on design, fabrication and characterization of active gas foil bearings, signal analysis and processing and rotor dynamics.

Yihua Wu received the Ph.D. degree from College of Mechanical and Vehicle Engineering, Hunan University, China, in 2019. His research includes design and control of porous tilting pad bearings, simulation of microgravity effects.

Kai Zhang received the Ph.D. degree in Huazhong University of Science and Technology, 2017, China. He is currently a lecturer in College of Mechanical and Vehicle Engineering, Hunan University, China. His current research includes design and characterization of air bearing with flexure pivot-tilting pads and rotor dynamics.