Solving the Coriolis Vibratory Gyroscope Motion Equations by Means of the Angular Rate B-Spline Approximation

Abstract

The exact solution of the movement equation of the Coriolis vibratory gyroscope (CVG) with a linear law of variation of the angular rate of rotation of the base is given. The solution is expressed in terms of the Weber functions (the parabolic cylinder functions) and their asymptotic representations. On the basis of the obtained solution, an analytical solution to the equation of the ring dynamics in the case of piecewise linear approximation of an arbitrary angular velocity profile on a time grid is derived. The piecewise linear solution is compared with the more rough piecewise constant solution and the dependence of the error of such approximations on the sampling step in time is estimated numerically. The results obtained make it possible to significantly reduce the number of operations when it is necessary to study long-range dynamics of oscillations of the system, as well as quantitatively and qualitatively control the convergence of finite-difference schemes for solving the movement equations of the Coriolis vibratory gyroscope.

1. Introduction

The functioning of Coriolis vibratory gyroscopes (CVG) is based on the effect of the precession of elastic waves excited in axisymmetric shells [1,2,3]. Taking into account the fact that the main part of the vibration energy of the elastic shell corresponds to the region adjacent to its edge, the ring-shaped model is the basic one for studying the dynamics of the CVG of any configuration.

It is especially simple to analyze the dynamics of a ring-shaped CVG resonator in the case of a constant angular rate. With certain dependencies, by reducing to the known second-order differential equations with variable coefficients, the solutions of the equations of dynamics can also be found analytically [4]. However, firstly, the form of these solutions can be quite cumbersome (expressed in terms of special functions). Secondly, not every law of rotation of the base allows one to obtain well-studied analytically differential equations. In particular, this applies to such practically important cases when the analytical dependence changes smoothly or abruptly with time, as well as if it is specified not in analytical but in tabular form.

Therefore, the numerical methods of solution are the most universal [5]. One of the most significant drawbacks of numerical methods is the need to find a compromise between the speed and accuracy of difference schemes. It is not possible to obtain a solution at some time moment without evaluation of the solutions at all previous time grid knots. To construct a continuous solution by the discrete set of values, additional approximation is required. Finally, it is almost impossible to get physical interpretation of the numerical solutions obtained by means of the finite difference method.

To preserve the advantages of analytical methods, which make it possible to quickly calculate the parameters of motion at any predetermined moment of time according to given initial conditions, and at the same time use the universality of numerical methods, the following numerical–analytical approach is proposed in this work. Its essence lies in the piecewise approximation of an arbitrary dependence with the subsequent exact solution of problems on each elementary segment. In this case, it is enough just to calculate the values of the linear displacement and acceleration of the points of the ring at the end of a subinterval, which are at the same time the initial conditions for the subsequent segment. In this paper, we first consider the simplest piecewise constant approximation of the angular rate. Then we move on to consider the piecewise linear approximation of an arbitrary angular rate, based on the analytical solution for a linear angular velocity on each subinterval, taking into account situations requiring the use of asymptotic expansions. The latter are important in the case of high-intensity dynamics of CVG [6,7]. In the experimental section, we compare both types of solution and show the benefits of the piecewise linear approximation.

2. Statement of the Problem

Unlike the variety of microelectromechanical systems (MEMS) designs based on lumped masses for sensing Coriolis acceleration [3], the ring-shaped CVG [8,9,10], as well as the hemispherical resonator gyroscope (HRG) [2] or the cylindrical resonator gyroscope (CRG) [11], are examples of devices using sensitive elements with distributed masses. An in-plane standing wave is excited in an elastic ring-shaped resonator oscillating along axes X and Y with primary displacement $x_1$ (Figure 1). In the absence of external rotation applied to the resonator, the primary wave has four nodes with zero ring displacement, which are located at 45 deg between axes X and Y. After applying a rotation $\overrightarrow{\Omega }$ around the sensitive axis, perpendicular to the plane of the ring, the primary wave begins to move around the ring due to the effect of the Coriolis inertia forces. As a result of this movement, secondary oscillations appear at the node points, thus making it possible to measure external angular rate projection to the sensitivity axis $\Omega$. This phenomenon of standing wave inertia in rotating axisymmetrical solid bodies was first observed by G. Bryan [12] and sometimes is called the Bryan effect.

In dimensionless form, the system of motion equations for free oscillations of such perfect ring-shaped CVG without dissipation can be written as [3]

$$\begin{array}{l}{\ddot{x}}_1+\left({\omega }^2-{\Omega }^2\right)x_1+2K\Omega {\dot{x}}_2+K\dot{\Omega }{x}_2=0,\\ {\ddot{x}}_2+\left({\omega }^2-{\Omega }^2\right)x_2-2K\Omega {\dot{x}}_1-K\dot{\Omega }{x}_1=0,\end{array}$$

(1)

where $x_1$ and $x_2$ are the primary and secondary displacements; $\omega$ is the natural frequency; $\Omega$ is the angular rate projection; and $K$ is the CVG’s coupling factor (the Bryan coefficient) corresponding to the excited mode of oscillations.

To obtain the unique solution to Equation (1), the initial conditions at time moment $t_0$ must be given:

$$\begin{array}{l}{x}_1(t_0)={\alpha }_1,\hspace{1em}{\dot{x}}_1(t_0)={\beta }_1;\\ x_2(t_0)={\alpha }_2,\hspace{1em}{\dot{x}}_2(t_0)={\beta }_2.\end{array}$$

(2)

Unlike a number of other works, in Equation (1) we do not assume that ${\Omega }^2\approx 0$ and/or $\dot{\Omega }\approx 0$. Instead, we will solve Equation (1) taking into account both the centrifugal forces corresponding to terms ${\Omega }^2{x}_i$ and the variable angular rate $\Omega =\Omega (t)$.

In the matrix form, Equation (1) can be written as

$$\ddot{x}+\left({\omega }^2-{\Omega }^2\right)x+2K\Omega J\dot{x}+K\dot{\Omega }Jx=0,$$

(3)

where

$$x=\left(\begin{array}{l}{x}_1\\ x_2\end{array}\right),\hspace{1em}J=\left(\begin{array}{cc}0& -1\\ 1& 0\end{array}\right).$$

(4)

According to [13], we introduce a variable phase $\gamma (t)$ and change the variables as

$$\begin{array}{l}{x}_1=y_1\cos 2\gamma -y_2\sin 2\gamma ,\\ x_2=y_1\sin 2\gamma +y_2\cos 2\gamma ,\end{array}$$

or

$$x=Ay,\hspace{1em}A=\left(\begin{array}{cc}\cos 2\gamma & -\sin 2\gamma \\ \sin 2\gamma & \cos 2\gamma \end{array}\right).$$

(5)

Obviously,

$$\frac{\mathrm{d}A}{\mathrm{d}\gamma }=-2\left(\begin{array}{cc}\sin 2\gamma & \cos 2\gamma \\ -\cos 2\gamma & \sin 2\gamma \end{array}\right)$$

(6)

and

$$\frac{{\mathrm{d}}^{2}A}{\mathrm{d}{\gamma }^2}=-4A,\hspace{1em}{A}^{\mathrm{T}}\frac{\mathrm{d}A}{\mathrm{d}\gamma }=2J,\hspace{1em}{A}^{\mathrm{T}}J\frac{\mathrm{d}A}{\mathrm{d}\gamma }=-2I,$$

(7)

where I is the identity matrix.

Differentiation gives us

$$\dot{x}=A\dot{y}+\dot{\gamma }\frac{\mathrm{d}A}{\mathrm{d}\gamma }y,\hspace{1em}\ddot{x}=A\ddot{y}+2\dot{\gamma }\frac{\mathrm{d}A}{\mathrm{d}\gamma }\dot{y}+\ddot{\gamma }\frac{\mathrm{d}A}{\mathrm{d}\gamma }y-4{\dot{\gamma }}^{2}Ay.$$

(8)

Substituting Equation (8) in Equation (3), after left multiplication by matrix $A^{\mathrm{T}}$, we get

$$\ddot{y}+\left(4\dot{\gamma }-2K\Omega \right)J\dot{y}+\left(2\ddot{\gamma }-K\dot{\Omega }\right)Jy+\left(4K\Omega \dot{\gamma }-4{\dot{\gamma }}^2+{\omega }^2-{\Omega }^2\right)y=0.$$

(9)

Assuming

$$\gamma (t)=\frac{K}{2}{\displaystyle \underset{t_0}{\overset{t}{\int }}\Omega (\tau )\mathrm{d}\tau },$$

(10)

we obtain the system of ordinary differential equations (ODEs)

$$\ddot{y}+{\omega }^2\left[1-\left(1-K^2\right)\frac{{\Omega }^2}{{\omega }^2}\right]y=0,$$

(11)

which can be solved independently for each scalar function $y_1,\hspace{0.17em}{y}_2$.

In previous works, researchers usually supposed that ${\Omega }^2/{\omega }^2<<1$ and neglect the second term in square brackets in Equation (11) even if the angular rate is not constant [3]. In this case, the solution of Equation (11) is simple and obtained in the form of trigonometric functions. However, if we need to study high-intensive dynamic processes with variable angular rates [5,6], we cannot omit the term $\left(1-K^2\right){\Omega }^2/{\omega }^2$ and must solve ODE (11) with a variable coefficient.

Considering that

$$\dot{\gamma }(t)=\frac{K}{2}\Omega (t),$$

(12)

initial conditions for the new variables will take the form

$$\begin{array}{l}{y}_1(t_0)={\alpha }_1,\hspace{1em}{\dot{y}}_1(t_0)={\beta }_1-K\Omega (t_0){\alpha }_1;\\ y_2(t_0)={\alpha }_2,\hspace{1em}{\dot{y}}_2(t_0)={\beta }_2-K\Omega (t_0){\alpha }_2.\end{array}$$

(13)

In particular, for the most important second form of in-plane oscillations with four nodes and four antinodes (Figure 1), the angle between the axes of secondary displacements, $x_1$ and $x_2$, is equal to $\pi /4$, and we have $K=-0.8$, $\omega =6/\sqrt{5}$.

Each ODE of system (11) takes the form

$$\ddot{y}+\frac{36}{5}\left(1-\frac{1}{20}{\Omega }^2\right)y=0.$$

(14)

It should be noted that for the HRG, CRG, or other CVG models it is possible to obtain movement equations analogous to Equation (11), but with corresponding coefficients before the term ${\Omega }^2/{\omega }^2$.

The ODE (14) can be easily solved exactly for the case of a constant angular rate. If the angular rate is variable, it is possible to make its rough piecewise constant (stepwise) approximation with recurrent redefinition of initial conditions (13). This idea is explained in the next section.

3. Stepwise Approximation

In the case of $\Omega (t)={\Omega }_0$, Equation (14) has two linearly independent solutions

$$y^{(1)}(t)=\cos \left(\frac{3}{5}t\sqrt{20-{\Omega }_0^2}\right),\hspace{1em}{y}^{(2)}(t)=\sin \left(\frac{3}{5}t\sqrt{20-{\Omega }_0^2}\right).$$

(15)

The variable phase (10) is equal to

$$\gamma (t)=-\frac{2}{5}{\Omega }_0\cdot (t-t_0).$$

(16)

The solution to Equation (14) has the form

$$y(t)=C_1{y}^{(1)}(t)+C_2{y}^{(2)}(t),$$

(17)

with undetermined coefficients $C_1,C_2$ found from initial conditions (13).

Let the angular rate $\Omega =\Omega (t)$ be defined on interval $t\in [t_0,T]$. Introduce the uniform mesh

$$j=0,1,\dots ,N,\hspace{1em}h=\frac{T-t_0}{N},\hspace{1em}{t}_j=t_0+jh,\hspace{1em}{\Omega }_j=\Omega (t_j).$$

(18)

The stepwise approximation of $\Omega (t)$ is expressed in the following way:

$$\tilde{\Omega }(t)={\displaystyle \sum _{j=0}^{N-1}{\Omega }_j{B}_0\left(\frac{t}{h}-j\right)}.$$

(19)

where

$$B_0(t)=\{\begin{array}{cc}1,& t\in [0,1),\\ 0,& t\notin [0,1)\end{array}$$

(20)

are the zero-order B-splines [14].

On each interval $t\in [t_j,t_{j+1}]$, we have

$$\tilde{\Omega }(t)={\Omega }_j\hspace{1em}(j=0,1,\dots ,N-1).$$

(21)

The general algorithm is to recursively obtain a solution of ODE (14) with constant angular rate (21) on each of subintervals $[t_j,t_{j+1}]$. On the start subinterval $[t_0,t_1]$, the initial conditions (13) are used, and to determine the initial conditions on each subinterval $[t_j,t_{j+1}]$, $j=1,2,\dots ,N-2$, it is necessary to use the values calculated on the previous subinterval $[t_{j-1},t_j]$.

Instead of coefficients ${\Omega }_j$ in Equation (19), one can take other values providing slightly better accuracy in different metrics, for example, $({\Omega }_j+{\Omega }_{j+1})/2$. However, the order of approximation accuracy here remains the same, i.e., $O(h)$.

4. Piecewise Linear Approximation

Unfortunately, the stepwise approximation (19) having the 1st order of accuracy with respect to the time step h is not sufficient for precise and fast measurements. Therefore, the high-order representations of $\Omega (t)$ should be used, for example, the piecewise linear approximation having the 2nd order of accuracy with respect to h.

The main difficulty in using the piecewise linear approximation is in the fact that in previous works researchers did not study analytically the case of linear angular rate. In addition to the trivial case of the constant angular rate considered in Section 3, one can also find the analytical solutions for the cases of angular rate changing according to the parabolic law (in the form of the Bessel functions) and the harmonic and poly-harmonic laws (in the form of Mathieu functions) [4].

In this regard, we will further consider in more detail the solution for the case of a linear angular rate

$$\Omega (t)=ct+b.$$

(22)

Substituting Equation (22) in Equation (14), we obtain

$$\ddot{y}+\frac{36}{5}\left[1-\frac{1}{20}{(ct+b)}^2\right]y=0.$$

(23)

Using a linear change of variable, we reduce Equation (23) to the standard form of the Weber equation for parabolic cylinder functions with respect to the new variable $z$ [15]:

$$\frac{{\mathrm{d}}^{2}y}{\mathrm{d}{z}^2}-\left(\frac{1}{4}{z}^2+\mu \right)y=0,$$

(24)

where

$$z=\sqrt{\frac{6c}{5}}\left(t+\frac{b}{c}\right),\hspace{1em}\mu =-\frac{6}{c}.$$

Even and odd linearly independent solutions to Equation (24) have the form

$$y^{(1)}(\mu ,z)={\displaystyle \sum _{k=0}^{\infty }a{(\mu )}_{2k}\frac{z^{2k}}{(2k)!}},\hspace{1em}{y}^{(2)}(\mu ,z)={\displaystyle \sum _{k=0}^{\infty }a{(\mu )}_{2k+1}\frac{z^{2k+1}}{(2k+1)!}},$$

(25)

where

$$a{(\mu )}_0=a{(\mu )}_1=1,\hspace{1em}a{(\mu )}_{n+2}=\mu a{(\mu )}_n+\frac{1}{4}n(n-1)a{(\mu )}_{n-2}\hspace{1em}(n=2,3,\dots )$$

At large times ($t\to \infty$) or at small angular rate values ($\left|c\right|\to 0$), it is necessary to keep a sufficiently large number of terms of strongly oscillating infinite power series (25). In this regard, a variety of asymptotic representations should be used. In most cases of practical importance, the Darwin expansions [15] are suitable fast enough representations. For example, if $\mu$ is large and negative and $z$ is moderate, we can take

$$\begin{array}{c}{y}^{(1)}(\mu ,z)=\frac{2\sqrt{\Gamma \left(\frac{1}{2}-\mu \right)}}{{(2\pi )}^{1/4}}{\mathrm{e}}^{v_r(\mu ,z)}\cos \left[\frac{1}{4}\pi +\frac{1}{2}\pi \mu +\theta (\mu ,z)+v_i(\mu ,z)\right],\\ y^{(2)}(\mu ,z)=\frac{2}{{(2\pi )}^{1/4}\sqrt{\Gamma \left(\frac{1}{2}-\mu \right)}}{\mathrm{e}}^{v_r(\mu ,z)}\sin \left[\frac{1}{4}\pi +\frac{1}{2}\pi \mu +\theta (\mu ,z)+v_i(\mu ,z)\right],\end{array}$$

(26)

where

$$\begin{array}{c}Z=\sqrt{4\left|\mu \right|-z^2},\hspace{1em}\theta (\mu ,z)=\frac{1}{4}zZ+|\mu |\arcsin \frac{z}{2\sqrt{\left|\mu \right|}},\\ v_r(\mu ,z)=-\frac{1}{2}\ln Z-\frac{{\delta }_6}{Z^6}+\frac{{\delta }_{12}}{Z^{12}}-\dots ,\hspace{1em}{v}_i(\mu ,z)=\frac{{\delta }_3}{Z^3}-\frac{{\delta }_9}{Z^9}+\dots \end{array}$$

Some of the first coefficients ${\delta }_{3k}$ are given by

$$\begin{array}{l}{\delta }_3=\frac{1}{\mu }\left(\frac{z^2}{48}+\frac{1}{2}\mu z\right),\hspace{1em}{\delta }_6=\frac{3}{4}{z}^2-2\mu ,\\ {\delta }_9=\frac{1}{{\mu }^3}\left(-\frac{7}{5760}{z}^9-\frac{7}{320}\mu z^7-\frac{49}{320}{\mu }^2{z}^5+\frac{31}{12}{\mu }^3{z}^3-19{\mu }^{4}z\right),\\ {\delta }_{12}=\frac{153}{8}{z}^4-186\mu z^2+80{\mu }^2.\end{array}$$

The variable phase (10) is equal to

$$\gamma (t)=-\frac{2}{5}(t-t_0)\left[\frac{c}{2}(t+t_0)+b\right].$$

(27)

Now we apply the piecewise linear approximation of arbitrary $\Omega (t)$ on the regular mesh (18):

$$\tilde{\Omega }(t)={\displaystyle \sum _{j=0}^N{\Omega }_j{B}_1\left(\frac{t}{h}-j\right)},$$

(28)

where the linear B-splines are [13]

$$B_1(t)=\{\begin{array}{cc}1-\left|t\right|,& \left|t\right|\le 1,\\ 0,& \left|t\right|\ge 1.\end{array}$$

(29)

It can be seen that for $t\in [t_j,t_{j+1}]$

$$\tilde{\Omega }(t)=c_{j}t+b_j\hspace{1em}(j=0,1,\dots ,N-1),$$

(30)

where

$$c_j=\frac{1}{h}({\Omega }_{j+1}-{\Omega }_j),\hspace{1em}{b}_j=\frac{1}{h}({\Omega }_j{t}_{j+1}-{\Omega }_{j+1}{t}_j).$$

(31)

The common piecewise linear approximation algorithm is analogous to that for the stepwise approximation considered in Section 3.

Obviously, if the dependence of the angular rate is initially piecewise linear, then it is possible to obtain an exact solution through the Weber functions by appropriately choosing a non-uniform approximation grid. In this case, on subintervals $[t_j,t_{j+1}]$ with ${\Omega }_j={\Omega }_{j+1}$, we use the simplest analytical solutions (17) expressed in terms of trigonometric functions (15) and on subintervals with nonzero angular acceleration we apply the representation (17) expressed, depending on the sign of the acceleration and its absolute value as well as the length of the time interval, in terms of power series (25), Darwin expansions (26), or other asymptotic expansions [15].

5. Numerical Experiment

As an example we will find the solution to Equation (1) on interval $t\in [0,5]$, with an angular rate increasing according to the law [7]:

$$\Omega (t)=1-{\mathrm{e}}^{-0.5t}.$$

We take the following boundary conditions (2):

$$x_1(0)=1,\hspace{1em}{x}_2(0)={\dot{x}}_1(0)={\dot{x}}_2(0)=0.$$

Two types of approximations on the grid with $N=5$ (grid width $h=1$) are shown on Figure 2.

On Figure 3 the dynamics of the primary displacement $x_1$ is shown, evaluated by different approaches. The results were compared with solution ${\tilde{x}}_1$ (bold line) obtained by the high-order finite-difference procedure on the dense grid with width $\Delta =1\cdot {10}^{-3}$ [7]. Relative errors were computed according to the formula

$$\epsilon =\sqrt{\frac{{\displaystyle {\int }_0^5{\left[x_1(\tau )-{\tilde{x}}_1(\tau )\right]}^2\mathrm{d}\tau }}{{\displaystyle {\int }_0^5{\tilde{x}}_1^2(\tau )\mathrm{d}\tau }}}$$

It is seen that the piecewise linear approximation with Darwin representation provides much better accuracy ($\epsilon =0.013$) in comparison with the stepwise representation ($\epsilon =0.316$). In turn, it should be noted that the piecewise representation in the form of power series (25) is accurate only on relatively short time intervals and demonstrates instability at $t>>1$. Moreover, the latter approach is also much more computationally expensive (the number of terms in time series was 40) in comparison with the Darwin asymptotic expressions (26).

This example demonstrates the benefits of the novel approach based on linear piecewise approximation of an angular rate in comparison with rough method of stepwise approximation. We must mention that although analytical solutions obtained earlier for angular rates changing according to parabolic and harmonic laws [4] can also be applied for constructing the solution for arbitrary angular rate profiles, they require the use of more complicated approximation tools: parabolic B-splines and Fourier series, respectively. In the latter case it is difficult to approximate the non-smooth angular rate dynamics due to the necessity to take a large number of terms and the Gibbs effect. Taking these considerations into account we can conclude that the best way to arbitrary CVG’s dynamics simulation is the combination of both approaches proposed in our work (stepwise and linear piecewise), depending on the behavior of the angular rate.

6. Conclusions

For the first time, an analytical solution for the perfect CVG motion equations with linear angular rate of rotation is obtained. It complements the class of already known analytical solutions obtained for constant, parabolic law and harmonic and poly-harmonic angular rates [4]. On the base of this solution, an original approach using a piecewise approximation of an arbitrary angular rate is proposed, allowing one to obtain a simple recurrent sequence of analytical solutions on each linear segment by recurrent recalculation of initial conditions. Thus, due to such semi-analytical solutions, we can significantly expand the class of analytical solutions of CVG motion equations. In comparison to numerical solutions, the new approach has a high speed of computations, which makes it possible to accurately simulate the dynamics of oscillations over long time intervals without accumulating an error. Using the method of piecewise-linear approximation, it is also possible to study the accuracy and convergence of various finite-difference schemes. Applying suitable asymptotic expansions, one can study dynamic processes of both low and high intensity. The algorithm will be especially effective when approximating on an irregular grid with small width h on intervals with strong change in angular velocity (large angular accelerations) and great h on relatively “smooth” intervals. Our approach can easily be applied to solve analogous motion equations for other CVG models, for example, HRG or CRG.