The modal analysis of disc resonator gyroscope using isogeometric analysis

Qingsi Cheng1 , Guolai Yang2 , Chen Lin3 , Jianli Ge4 , Quanzhao Sun5

1, 2, 3, 4, 5School of Mechanics, Nanjing University of Science and Technology, Nanjing, 210094, China

2Corresponding author

Journal of Vibroengineering, Vol. 19, Issue 1, 2017, p. 468-476. https://doi.org/10.21595/jve.2016.17431
Received 19 July 2016; received in revised form 8 November 2016; accepted 21 November 2016; published 15 February 2017

Copyright © 2017 JVE International Ltd. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
Creative Commons License
Table of Contents Download PDF References
Cite this article
Views 188
Reads 79
Downloads 937
Abstract.

This paper presents a free vibration analysis of the disc resonator gyroscope (DRG) based on the Isogeometric Analysis (IGA). The non-uniform rational B-splines (NURBS) basis functions are employed to represent both the exact geometry of DRG and displacement fields. It provides a unified process between the geometric design and analysis which significantly reduces the design cycle by eliminating the re-meshing process. Moreover, the modeling process of IGA and the numerical example of three-circle cascaded model with four different refinement schemes are presented and the computational accuracy and efficiency of the present numerical method are verified by comparing with the standard finite method (ABAQUS).

Keywords: NURBS, IGA, modal analysis, disc resonator gyroscope (DRG).

1. Introduction

With the widespread application of microelectromechanical system (MEMS) gyros especially in automotive safety and military tactical, researchers keep trying to find the high performance MEMES gyros [1-5]. The Hemispherical Resonator Gyroscope (HRG) is known as one of the highest performance gyroscopes thanks to inherent advantages of a degenerate wineglass mode operation such as thermal stability, electrostatic tuning capability, and high quality factor [6, 7]. However, the early resonator designs set stringent requirements on complex control strategy such as close-loop force feedback and mode-tuning, ultra-low noise electronics design, high cost 3D assembly technique in improving the capacitive sensitivity. In addition, the gyro bias was not stable due to the high dynamic load and imprecision of the bonded joints when running long. Given this background, the conception of the Disc Resonator Gyroscope (DRG) was proposed which yielded a compact planar micro-machined design with central support and carrying no critical loads [8-10].

The optimization of structural parameters and material selection has made contributions to achieve much stable performance, and this work is still in progress. The current approach of modal analysis is usually based on the FEM software including three steps: building the geometric model, transforming the geometry into mesh and making the modal analysis. It means that in every optimization cycle, the model will be modified once and re-meshed to analyze which brings a lot of repetitive work.

The concept of Isogeometric Analysis was firstly proposed by Hughes et al. [11] in 2005. It seeks to unify the fields of the computer-aided design (CAD) and the finite element analysis (FEA). The root idea is that the basis used to exactly model the geometry will also serve as the basis for the solution space of the numerical method. On this base, the model used for analysis is automatically created from CAD geometry and the procedure from geometry to mesh is skipped which reduces the repetitive work. On the other hand, the mesh refinement in IGA does not change the geometry which ensures the unification of geometry and analysis [12-14]. The application of IGA to structural vibration was firstly introduced by Cottrell et al. [15]. Since then, a various number of researches employing IGA approach have been presented to solve the dynamic problem in a wide range of research areas.

In this context, the advantages of IGA approach are employed to determine the intrinsic frequencies and the mode shapes for the DRG. When designing a DRG, the overall structure will not change except the size and material parameters. It provides convenience for the optimization because the procedure can be parameterized by using IGA. Moreover, the NURBS can exactly describe the circular geometry even at coarse mesh level. The main body of the DRG is composed of rings which means the IGA will offer an effective and accurate calculation.

In this work, we don’t focus on how to design a DRG with good performance but the modal simulation of DRG using IGA which is the basis of optimization. This paper is organized as follows:

Section 2 gives a brief introduction to the NURBS discretization and the free vibration analysis. Section 3 presents the structure of DRG and the analysis model. Section 4 shows the simulation results and comparison and Section 5 concludes this paper.

2. NURBS-based discretization and Free vibration analysis

Here we briefly review the concept of NURBS. For further details and extensive references, the reader is referred to Piegl and Tiller and Cottrell et al. [16-20].

A NURBS curve Cu can be described as:

(1)
C u = i = 1 n R i , p u P i ,
(2)
R i , p u = N i , p u w i i = 1 n N i , p u w i ,

where Ri,pu are the NURBS basis functions, Pi are the control points, wi are the weights of the control points, U=u1=0,,ui,,un+p+1=1Tis the knot vector in the parameter space, p is the polynomial order [21].

Taking the form of Eq. (1), the related reference coordinate and displacement located in the control point can be defined as:

(3)
x = A = 1 n c p R A x A ,
(4)
d = A = 1 n c p R A d A ,

where ncp is the number of control points, RA is the NURBS basis function corresponding to control point A.

For 2D problem, the Green strain is expressed as:

(5)
ε = B d ,
(6)
B = J - 1 R A u 1 0     R A u 2     0 R A u 2 R A u 1 T ,
(7)
J = ( x , y ) ( u 1 , u 2 ) = i = 1 n c p R i u 1 x i i = 1 n c p R i u 2 x i i = 1 n c p R i u 1 y i i = 1 n c p R i u 2 y i ,

where B is the strain-displacement matrix, J is the Jacobian mapping physical and parametric domains. (xi,yi) are the coordinates of control points.

Thus, for a body defined in the domain Ω, we can get the stiffness matrix and the mass matrix as:

(8)
K = Ω B T D B t d x d y ,
(9)
M = ρ Ω B T B t d x d y ,

where D is the elastic matrix, ρ is the material density and t is the thickness.

To compute the natural frequency, we ignore the damping matrix and the load vector. Thus, the global equation of motion for the DRG system can be written as:

(10)
M U ¨ + K U = 0 ,

where M is the global mass matrix, K is the global stiffness matrix and U={u11,u12,,un1,un2} is the global displacement vector (the superscript of ui1,2 is the direction of DOF and n is the number of control points).

By assuming U=eiωtU-, the motion equation is:

(11)
K - ω 2 M U - = 0 .

Then, it can be translated into the eigenvalue problem:

(12)
K - ω 2 M = 0 ,

where U- is the associated displacement vector and the solution of Eq. (12) is the natural frequency.

3. DRG model

The DRG is constructed by a series of rings and connected with spokes. Fig. 1(a) gives the picture of real product and Fig. 1(b) shows a design schematic diagram of DRG [7]. In the actual model, there is a fixed solid disk in the center as shown in Fig. 1(b). For the analysis model, we fixed the most inside spokes instead and set it as the displacement boundary condition as shown in Fig. 2(a). The details of the control points and the mesh are shown in Fig. 2(b) and Fig. 2(c).

Fig. 1. The structure of DRG: a) shows a real model, for b), the solid blue lines represent spokes that are fixed, while the dashed red lines represent spokes that are variable

The structure of DRG: a) shows a real model, for b), the solid blue lines represent spokes  that are fixed, while the dashed red lines represent spokes that are variable

a) SEM image of epi-sealed DRG

The structure of DRG: a) shows a real model, for b), the solid blue lines represent spokes  that are fixed, while the dashed red lines represent spokes that are variable

b) The design sketch

To obtain the analysis model in the traditional FEA, we should draw the model first and transform it into mesh then. If we want to change the number of the rings or other parameters, these two steps have to be repeated. Here we give a detailed introduction of the modeling process using IGA which is parameterized.

The tensor product structure of the parameter space of a patch makes it poorly suited for representing complex, multiply connected domains [12]. The DRG can be handled by using multiple patches as shown in Fig. 2(c).

Fig. 2. The sketchy analysis model: a) is the diagram, b) shows the control points (blue points), c) shows the multiple patches and the green lines constitute the mesh

The sketchy analysis model: a) is the diagram, b) shows the control points (blue points),  c) shows the multiple patches and the green lines constitute the mesh

a)

The sketchy analysis model: a) is the diagram, b) shows the control points (blue points),  c) shows the multiple patches and the green lines constitute the mesh

b)

The sketchy analysis model: a) is the diagram, b) shows the control points (blue points),  c) shows the multiple patches and the green lines constitute the mesh

c)

Table 1. The number of patches

Parameter
Value
Number of rings (Nr)
N
Number of patches inside the odd ring (N1)
8 * N + M o d N , 2 / 2
Number of patches inside the even ring (N2)
8 * N - M o d N , 2 / 2
Number of patches for all the rings (N3)
32 * N
Number of total patches (NT)
N 1 + N 2 + N 3

Table 2. Design parameters

Parameter
Value
Internal diameter of the inner ring (R0)
566 mm
Ring width (nominal, Wrn )
7 mm
Electrode gap size (Wg)
3 mm
Spoke width (Ws)
7 mm
Variable spoke angle (α)
22.5°
Young’s modulus (E)
169000 MPa
Poisson’s ratio (γ)
0.3
Density (ρ)
2.33E6 Kg/m3

If we set every spoke as a patch, there are 8 patches inside every ring which can be divided into 32 patches. Table 1 lists the number of patches associated with the number of rings. Table 2 gives the parameters set in the simulation model of DRG.

The internal diameter of each ring can be derived based on the ring width and the electrode gap size, then, the control points’ coordinates of all the patches can also be calculated. Due to the model is completely based on the parameters given in Table 1 and Table 2, we only need to change the parameters at the beginning of the model process and the model modification will be automatically generated.

Fig. 3 shows the analysis models with different number of rings (10 rings and 35 rings) and the only thing to change is the parameter Nr in the whole modeling process. For the convenience of calculation and the modeling in ABAQUS, we choose the model of 3 rings in the numerical simulation as an example.

Fig. 3. a) is the model of 10 rings, b) is the model of 35 rings

a) is the model of 10 rings, b) is the model of 35 rings

a)

a) is the model of 10 rings, b) is the model of 35 rings

b)

4. Numerical results and comparison

In order to verifies the IGA model, ABAQUS is selected to compare the simulation results between IGA and FEM. In the comparison, four different meshes are employed in each method, the first ten natural frequencies and mode shapes are computed and we assume that the result of highly refined model (387104 nodes) in ABAQUS is the standard value. The results are separately shown in Table 3 and Table 4 (The unit of all the frequencies is Hz).

Table 3. The first ten natural frequencies using ABAQUS

Mode
DOF
12918
28478
55036
199652
Standard value
1
589.9
7449.8
7817.9
8132.7
8159.0
2
667.9
7751.4
8159.5
8498.7
8565.3
3
1014.0
7792.4
8167.2
8498.9
8565.9
4
1157.3
8449.7
8929.2
9256.1
9321.6
5
1263.6
8506.7
8958.3
9348.6
9428.2
6
1491.5
8515.3
8966.6
9350.1
9428.7
7
1752.4
9032.2
9550.8
9896.8
9966.7
8
1957.9
9076.6
9553.5
9899.2
9966.9
9
2164.1
9349.3
9892.4
10353.0
10450.0
10
2307.0
9365.1
9896.1
10354.0
10450.0

For further comparison, we choose the value of the frequency of fundamental, 2nd order, 4th order and 10th order harmonics, respectively. The correlation curves are shown in Fig. 4(a), (b), (c) and (d). Because of the large span of the abscissa, we take the format of log 2. The first ten mode shapes are shown in Fig. 5. The relative error of the second frequency between IGA and FEM is given in Table 5.

Table 4. The first ten natural frequencies using IGA

Mode
DOF
4176
6480
8640
15360
Standard value
1
8599.8
8327.9
8284.7
8104.0
8159.0
2
8925.7
8666.2
8633.1
8507.0
8565.3
3
8935.5
8674.3
8640.1
8512.1
8565.9
4
9695.7
9479.6
9462.9
9225.0
9321.6
5
9701.6
9483.6
9464.8
9436.5
9428.2
6
9931.6
9564.1
9528.4
9438.4
9428.7
7
10598.6
10243.0
10181.5
9646.9
9966.7
8
10602.3
10325.3
10280.0
10017.6
9966.9
9
10644.1
10485.4
10474.8
10453.0
10450.0
10
10731.6
10488.7
10477.5
10460.9
10450.0

Fig. 4. The correlation curves in four different modes

The correlation curves in four different modes

a) Mode 1

The correlation curves in four different modes

b) Mode 2

The correlation curves in four different modes

c) Mode 5

The correlation curves in four different modes

d) Mode 10

Obviously, we find that, the DOF in IGA is much less than FEM to approach the standard value. In the finite element method, the DOF will directly affect the solution time. Therefore, it proves that the computational efficiency of IGA is higher than FEM. Table 5 also shows the high accuracy of IGA when the DOF is relatively less.

Table 5. The relative error of the frequency of second harmonic between IGA and FEM

DOF
IGA
DOF
FEM
4176
4.21 %
12918
–92.2 %
6480
1.18 %
28478
-9.5%
8640
0.78 %
55036
–4.74 %
15360
–0.68 %
199652
–0.78 %

Fig. 5. The first ten mode shapes. a)-j) show the mode shapes from 1st to 10th, the left one of each part is the result of IGA (drawn by Tecplot), the right one is the result of ABAQUS

 The first ten mode shapes. a)-j) show the mode shapes from 1st to 10th, the left one of each part is the result of IGA (drawn by Tecplot), the right one is the result of ABAQUS

a1)

 The first ten mode shapes. a)-j) show the mode shapes from 1st to 10th, the left one of each part is the result of IGA (drawn by Tecplot), the right one is the result of ABAQUS

a2)

 The first ten mode shapes. a)-j) show the mode shapes from 1st to 10th, the left one of each part is the result of IGA (drawn by Tecplot), the right one is the result of ABAQUS

b1)

 The first ten mode shapes. a)-j) show the mode shapes from 1st to 10th, the left one of each part is the result of IGA (drawn by Tecplot), the right one is the result of ABAQUS

b2)

 The first ten mode shapes. a)-j) show the mode shapes from 1st to 10th, the left one of each part is the result of IGA (drawn by Tecplot), the right one is the result of ABAQUS

c1)

 The first ten mode shapes. a)-j) show the mode shapes from 1st to 10th, the left one of each part is the result of IGA (drawn by Tecplot), the right one is the result of ABAQUS

c2)

 The first ten mode shapes. a)-j) show the mode shapes from 1st to 10th, the left one of each part is the result of IGA (drawn by Tecplot), the right one is the result of ABAQUS

d1)

 The first ten mode shapes. a)-j) show the mode shapes from 1st to 10th, the left one of each part is the result of IGA (drawn by Tecplot), the right one is the result of ABAQUS

d2)

 The first ten mode shapes. a)-j) show the mode shapes from 1st to 10th, the left one of each part is the result of IGA (drawn by Tecplot), the right one is the result of ABAQUS

e1)

 The first ten mode shapes. a)-j) show the mode shapes from 1st to 10th, the left one of each part is the result of IGA (drawn by Tecplot), the right one is the result of ABAQUS

e2)

 The first ten mode shapes. a)-j) show the mode shapes from 1st to 10th, the left one of each part is the result of IGA (drawn by Tecplot), the right one is the result of ABAQUS

f1)

 The first ten mode shapes. a)-j) show the mode shapes from 1st to 10th, the left one of each part is the result of IGA (drawn by Tecplot), the right one is the result of ABAQUS

f2)

 The first ten mode shapes. a)-j) show the mode shapes from 1st to 10th, the left one of each part is the result of IGA (drawn by Tecplot), the right one is the result of ABAQUS

g1)

 The first ten mode shapes. a)-j) show the mode shapes from 1st to 10th, the left one of each part is the result of IGA (drawn by Tecplot), the right one is the result of ABAQUS

g2)

 The first ten mode shapes. a)-j) show the mode shapes from 1st to 10th, the left one of each part is the result of IGA (drawn by Tecplot), the right one is the result of ABAQUS

h1)

 The first ten mode shapes. a)-j) show the mode shapes from 1st to 10th, the left one of each part is the result of IGA (drawn by Tecplot), the right one is the result of ABAQUS

h2)

 The first ten mode shapes. a)-j) show the mode shapes from 1st to 10th, the left one of each part is the result of IGA (drawn by Tecplot), the right one is the result of ABAQUS

i1)

 The first ten mode shapes. a)-j) show the mode shapes from 1st to 10th, the left one of each part is the result of IGA (drawn by Tecplot), the right one is the result of ABAQUS

i2)

 The first ten mode shapes. a)-j) show the mode shapes from 1st to 10th, the left one of each part is the result of IGA (drawn by Tecplot), the right one is the result of ABAQUS

j1)

 The first ten mode shapes. a)-j) show the mode shapes from 1st to 10th, the left one of each part is the result of IGA (drawn by Tecplot), the right one is the result of ABAQUS

j2)

5. Conclusions

In this work, the Isogeometric analysis based on NURBS is applied in the modal simulation for the disc resonator gyroscope (DRG). It provides two advantages:

• NURBS geometry is more efficient in describing circular structure and it provides higher accuracy and efficiency.

• It’s uniform between design and analysis which means there is no need to re-mesh the model in the optimization period. It can save a lot of repeated work in the optimization design.

Four different level of mesh refinements are computed by using IGA and ABAQUS in the numerical example. The results of first ten frequencies and mode shapes are analyzed to compare the efficiency and accuracy. The method of Isogeometric analysis shows great advantages.

The extension to three-dimension analysis is the ongoing research [22]. Also, we are trying to design a better performance DRG by using IGA and compare with physical experiment.

References

  1. Ahn C. H., et al. Geometric compensation of (100) single crystal silicon disk resonating gyroscope for mode-matching. Proceedings of 17th International Conference on Solid-State Sensors, Actuators, Microsystems, 2013, p. 1847-1850. [Publisher]
  2. Yazdi N., Ayazi F., Najafi K. Micromachined inertial sensor. Proceedings of IEEE, Vol. 86, Issue 8, 1998, p. 1640-1659. [Publisher]
  3. Shkel A. M. Type I and Type II micromachined vibratory gyroscopes. Position, Location, And Navigation Symposium, IEEE/ION IEEE, 2006, p. 586-593. [Search CrossRef]
  4. Shkel A. M. Microtechnology comes of age. GPS World, Vol. 22, Issue 9, 2011, p. 43-50. [Search CrossRef]
  5. Schleglov K. V., Challoner A. D. Isolated Planar Gyroscope with Internal Radial Sensing and Actuation. U.S. Patent 7,040,163, 2006. [Search CrossRef]
  6. Rozelle D. M. The hemispherical resonator gyro: from wineglass to the planets. Advances in the Astronautical Sciences, Vol. 134, 2009, p. 1157-1178. [Search CrossRef]
  7. Ahn C. H., Ng E. J., Hong V., et al. Mode-matching of wineglass mode disk resonator gyroscope in (100) single crystal silicon. Journal of Microelectromechanical Systems, Vol. 24, Issue 2, 2015, p. 343-350. [Publisher]
  8. Challoner D. A., Shcheglov K. V. Isolated Planar Mesogyroscope. US 7168318 B2, 2007. [Search CrossRef]
  9. Kubena R. L., Chang D. T. Disc Resonator Gyroscopes. US 7581443 B2, 2009. [Search CrossRef]
  10. Challoner A. D., Ge H. H., Liu J. Y. Boeing disc resonator gyroscope. Position, Location and Navigation Symposium, IEEE/ION IEEE, 2014, p. 504-514. [Publisher]
  11. Hughes T. J. R., Cottrell J. A., Bazilevs Y. Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement. Computer Methods in Applied Mechanics and Engineering, Vol. 194, 2005, p. 4135-4195. [Publisher]
  12. Cottrell J. A., Hughes T. J. R., Bazilevs Y. Isogeometric Analysis: Toward Integration of CAD and FEA. Wiley and Sons, Vol. 5, 2009, p. 299-325. [Publisher]
  13. Cottrell J. A., Hughes T. J. R., Bazilevs Y. Isogeometric Analysis. Wiley, 2009. [Publisher]
  14. Hughes T. J. R., Reali A., Sangalli G. Efficient quadrature for NURBS-based isogeometric analysis. Computer Methods in Applied Mechanics and Engineering, Vol. 199, 2010, p. 301-313. [Publisher]
  15. Cottrell J. A., Reali A., Bazilevs Y., et al. Isogeometric analysis of structural vibrations. Computer Methods in Applied Mechanics and Engineering, Vol. 195, Issue 41, 2006, p. 5257-5296. [Publisher]
  16. Piegl L. A., Tiller W. The NURBS Book. Springer, Berlin, 1995. [Publisher]
  17. Elguedj T., Bazilevs Y., Calo V. M., Hughes T. J. R. B and F projection methods for nearly incompressible linear and non-linear elasticity and plasticity using higher-order NURBS elements. Computer Methods in Applied Mechanics and Engineering, Vol. 197, 2008, p. 2732-2762. [Search CrossRef]
  18. Bazilevs Y., Calo V. M., Cottrell J. A., et al. Isogeometric analysis using T-splines. Computer Methods in Applied Mechanics and Engineering, Vol. 199, Issue 199, 2010, p. 229-263. [Publisher]
  19. Cottrell J. A., Hughes T. J. R., Reali A. Studies of refinement and continuity in isogeometric structural analysis. Computer Methods in Applied Mechanics and Engineering, Vol. 196, Issue 41, 2007, p. 4160-4183. [Publisher]
  20. Dörfel M. R., Jüttler B., Simeon B. Adaptive isogeometric analysis by local h -refinement with T-splines. Computer Methods in Applied Mechanics and Engineering, Vol. 199, Issues 5-8, 2010, p. 264-275. [Publisher]
  21. Cheng Q., Yang G., Lu J. An analysis of gear based on isogeometric analysis. Vibroengineering Procedia, Vol. 2, 2013, p. 17-22. [Search CrossRef]
  22. Xu G., Mourrain B., Duvigneau R., et al. Optimal analysis-aware parameterization of computational domain in 3D isogeometric analysis. Computer-Aided Design, Vol. 6130, Issue 4, 2010, p. 236-254. [Search CrossRef]

Cited By

Journal of Pressure Vessel Technology
Qingbo Yu, Guolai Yang
2019
Journal of Vibroengineering
Qingsi Cheng, Guolai Yang, Chencheng Yu, Quanzhao Sun
2017