Analytical Analysis of Magnetic Levitation Systems with Harmonic Voltage Input

Abstract

In this paper, a new analytical method using Lagrange equations for the analysis of magnetic levitation (MagLev) systems is proposed, using Thomson’s jumping ring experiment. The method establishes the dependence of the primary and induced currents, and also the equilibrium height of the levitating object on the input voltage through the mutual inductance of the system. The mutual inductance is calculated in two ways: (i) by employing analytical formula; (ii) through an improved semi-empirical formula based on both measurements and analytical results. The obtained MagLev model was analyzed both analytically and numerically. Analytical solutions to the resulting equations were found for the case of a dynamic equilibrium. The numerical results obtained for the dynamical model under transient operation show a close correspondence with the experimental results. The good precision of the analytical and numerical results demonstrates that the developed method can be effectively implemented.

1. Introduction

The principles of magnetic levitation (MagLev) are used in many important applications, such as actuators, power interrupters for eliminating electric arcs, electromagnetic shock absorbers in electrical vehicles, hybrid suspension systems with active control, levitation of superconductor materials, sensors and electromagnetic mass drivers [1]. One of the perspective usages of the MagLev systems in the near future is related to electromagnetic suspension systems and damping shock absorbers. Automobiles and trucks have shock absorbers with similar geometry to that of Thomson’s ring system. These systems are used to damper vibrations generated by imperfections of road surfaces [2,3,4,5].

Thus, it is an important challenge to be able mathematically and numerically model, analyze, control and design MagLev devices similar to Thomson’s ring system. There are also highly accurate density measurements and density-based detection methods applied to MagLev systems with similar geometries. In [6,7,8], an experimental study of a single-ring MagLev system was presented, opening up a wide operational space and enabling object manipulation and density-based measurements.

In the literature, there are works focused on the numerical analysis of MagLev systems similar to that of Thomson’s ring. Vilchis et al. [9] presented factors that affect 2D finite element method (FEM) transient simulation accuracy for an axisymmetric model of a Thomson-coil actuator. Yannan Zhou et al. [10] numerically studied a Thomson-coil actuator applied in a circuit breaker in a high-voltage direct current system. Tezuka et al. [11] proposed a novel magnetically-levitating motor, able to control five degrees of freedom. In all these systems, analytical or hybrid solutions can be implemented, resulting in accurate and low-demanding computational models.

Dynamic models (DM) of the MagLev systems have been previously developed. TJ-Teo et al. [12] presented a new electromagnetic actuator which has two degrees of freedom (linear and rotary motions). Consequently, the developed MagLev system presents two mutual inductances as functions of linear and rotatory positions. Additionally, Nagai et al. [13] realized a compact actuator system by a sensor-free approach. They presented the circuit equations for the solenoid, where the mutual inductance was approximated as a second-order polynomial function.

Here, we present a new and rigorous analytical solution to the dynamic equations of the MagLev systems, wherein a mutual inductance, which is the key parameter of electromagnetic force, is calculated analytically. This way, we obtain a method that retains the simplicity of lumped-parameter models plus the accuracy of analytical methods within a field modeling context. Our approach was validated using experimental test results obtained from a Thomson’s jumping ring experiment.

2. Model

In this work, a system known as the Thomson’s jumping ring is studied, where an aluminum ring is used as a levitating object. The experimental setup is shown in Figure 1. The system is composed of a primary copper coil of 1140 turns, and an aluminum ring that moves freely up and down along a laminated ferromagnetic core.

The equivalent circuit depicted in Figure 2 models this MagLev system.

The governing differential equations that govern the dynamics of the circuit represented in Figure 2 can be obtained from the Lagrangian

$$\begin{array}{c}\hfill \mathcal{L}=\underset{\mathrm{kinetic}\mathrm{energy}}{\underbrace{\frac{l_1{i}_1^2}{2}+\frac{l_2{i}_2^2}{2}+M{i}_1{i}_2+\frac{m{\dot{z}}^2}{2}}}-\underset{\mathrm{potential}\mathrm{energy}}{\underbrace{mgz}},\hfill \end{array}$$

(1)

where z is the ring height, and $i_1$ and $i_2$ are the electric currents in the primary coil and aluminum ring, respectively. The parameters $l_1$ and $l_2$ are the primary and secondary circuits’ inductances; M is the mutual inductance between the coil and ring; and m and g are the aluminum ring mass and gravitational constant, respectively. Notice that the mutual inductance $M\left(z\right)$ is a function of the distance z between the coil and ring [14]. Thus, the Lagrange equations are as follows:

$$\begin{array}{ccc}& & \frac{d}{dt}\frac{\partial \mathcal{L}}{\partial i_1}-\frac{\partial \mathcal{L}}{\partial q_1}=V\left(t\right)-R_1{i}_1,\hfill \end{array}$$

(2)

$$\begin{array}{ccc}& & \frac{d}{dt}\frac{\partial \mathcal{L}}{\partial i_2}-\frac{\partial \mathcal{L}}{\partial q_2}=-R_2{i}_2,\hfill \end{array}$$

(3)

$$\begin{array}{ccc}& & \frac{d}{dt}\frac{\partial \mathcal{L}}{\partial \dot{z}}-\frac{\partial \mathcal{L}}{\partial z}=-\beta \dot{z},\hfill \end{array}$$

(4)

where $R_1$ and $R_2$ are the coil and ring resistances; $q_1$ and $q_2$ are the coil and ring electric charges, respectively; and $\beta$ is the friction coefficient. The right-hand side terms of the Lagrange equations are generalized forces, which are considered as independent terms, associated with the generalized coordinate. Substitution of the Lagrangian (1) into Lagrange equations results in the following system of differential equations (see also [15,16]):

$$\begin{array}{ccc}& & l_1\frac{d{i}_1}{dt}=V\left(t\right)-R_1{i}_1-\frac{d\left(M\left(z\right)i_2\right)}{dt}\hfill \end{array}$$

(5)

$$\begin{array}{ccc}& & l_2\frac{d{i}_2}{dt}=-R_2{i}_2-\frac{d\left(M\left(z\right)i_1\right)}{dt},\hfill \end{array}$$

(6)

$$\begin{array}{ccc}& & m\ddot{z}=i_1{i}_2\frac{dM\left(z\right)}{dz}-\beta \dot{z}-mg.\hfill \end{array}$$

(7)

3. Formula for the Mutual Inductance

Mutual inductance between a coil and aluminum ring is a function of the system geometry and the magnetic medium properties between two coils. It can be calculated in two ways: numerically (using FEM) or analytically. In this section, the linear mutual inductance is analytically calculated based on the results of [17], where an infinite isotropic core is considered. According to [17], the mutual inductance can be divided into two parts:

$$\begin{array}{c}\hfill M=M_{\mathrm{air}}+M_{\mathrm{core}},\end{array}$$

(8)

where $M_{\mathrm{air}}$ is the mutual inductance in air, i.e., in the absence of the core, and $M_{\mathrm{core}}$ is the contribution of the ferromagnetic core to the whole inductance. The first part of the mutual inductance can be represented in the form (see [17]):

$$\begin{array}{c}\hfill M_{\mathrm{air}}={\mu }_0\frac{N_a}{w_a{h}_a{w}_b{h}_b}\underset{a}{\overset{a+h_a}{\int }}dA\underset{b}{\overset{b+h_b}{\int }}dB\underset{-w_a}{\overset{0}{\int }}d{z}_a\underset{z}{\overset{z+w_b}{\int }}d{z}_b\times \frac{2\sqrt{AB}}{f}\left[\left(1-\frac{f^2}{2}\right)K\left(f\right)-E\left(f\right)\right],\hfill \end{array}$$

(9)

where

$$\begin{array}{c}\hfill f=\sqrt{\frac{4r_a{r}_b}{{\left(z_b-z_a\right)}^2+{\left(r_a+r_b\right)}^2}}\end{array}$$

(10)

is a dimensionless parameter; $K\left(f\right)$ and $E\left(f\right)$ are complete elliptic integrals of the first and second kind, respectively [18]. Here, $(r_a,z_a)$ and $(r_b,z_b)$ are pairs of cylindrical coordinates inside the coil and ring (subindices a and b correspond to the coil and aluminum ring, respectively). Namely, $r_a$ and $r_b$ are radii-vectors and $z_a$ and $z_b$ positions on the core. Integration with respect to these variables signifies that the whole volume of the coil and ring is involved in the calculation of the mutual inductance between them. Parameters $w_a$ and $w_b$ are the coil and ring widths; $h_a$ and $h_b$ are their heights; and a and b are their internal radii (see Figure 3). $N_a$ is the number of turns in the coil. The mutual inductance caused by the core is as follows:

$$\begin{array}{ccc}\hfill M_{\mathrm{core}}& =& 2{\mu }_0\frac{N_a}{w_a{h}_a{w}_b{h}_b}\underset{0}{\overset{\infty }{\int }}\frac{dk}{k^2}\hfill \\ & \times & \left\{\cos \left(k(z+w_b)\right)+\cos \left(k(z+w_a)\right)-\cos \left(kz\right)-\cos \left(k(z+w_a+w_b)\right)\right\}\hfill \\ & \times & \frac{\pi }{2k}\left\{(a+h_a)[K_1\left(k(a+h_a)\right)L_0\left(k(a+h_a)\right)\right.+K_0\left(k(a+h_a)\right)L_1\left(k(a+h_a)\right)]\hfill \\ & -& \left.a\left[K_1\left(ka\right)L_0\left(ka\right)+K_0\left(ka\right)L_1\left(ka\right)\right]\right\}\frac{\pi }{2k}\left\{(b+h_b)[K_1\left(k(b+h_b)\right)L_0\left(k(b+h_b)\right)\right.\hfill \\ & +& K_0\left(k(b+h_b)\right)L_1\left(k(b+h_b)\right)]\left.-b\left[K_1\left(kb\right)L_0\left(kb\right)+K_0\left(kb\right)L_1\left(kb\right)\right]\right\}\hfill \\ & \times & \frac{{\mu }_{r}k{I}_1\left(\beta R\right)I_0\left(kR\right)-\beta I_0\left(\beta R\right)I_1\left(kR\right)}{{\mu }_{r}k{I}_1\left(\beta R\right)K_0\left(kR\right)-\beta I_0\left(\beta R\right)K_1\left(kR\right)},\hfill \end{array}$$

(11)

where $\beta =\sqrt{k^2+j\omega \mu \sigma }$, $\mu ={\mu }_0{\mu }_r$ is the absolute permeability of the core material, ${\mu }_r$ is its relative permeability and $\sigma$ is core conductance.

The total mutual inductance is a function of the ring position z; i.e., $M=M\left(z\right)$. Equations (9) and (11) are too complicated to be implemented in the dynamical analysis of the MagLev system. Equations (9) and (11), and FEM calculations, are too complicated and time-consuming to be implemented in the dynamical analysis of the MagLev system, especially in real-time calculations for an electromagnetic suspension system. Therefore, equations (9) and (11) must be simplified. In Figure 4, analytical behavior of the total inductance of the ring height is depicted. It can be seen in Figure 4 that the total mutual inductance of the system can be approximated by exponential function as follows:

$$\begin{array}{c}\hfill M\left(z\right)={\alpha }_1{e}^{-{\beta }_{1}z}+{\alpha }_2{e}^{-{\beta }_{2}z},\end{array}$$

(12)

where ${\alpha }_1$, ${\alpha }_2$, ${\beta }_1$ and ${\beta }_2$ are some empirical parameters to be obtained by fitting function (12) to the analytical result. The reason for the use of (12), apart from reproducing qualitatively theoretical results (9) and (11), is that it provides a correct asymptotic behavior of the mutual inductance at large distances. A polynomial approximation would diverge at high distances. The numerical values of the parameters involved in (12) can be obtained from (9) and (11) through the widely-used least squares method [19]. As a result, the following values can be obtained: ${\alpha }_1=0.000464$ H, ${\alpha }_2=0.000088$ H, ${\beta }_1=4.77562$ m⁻¹ and ${\beta }_2=22.74612$ m⁻¹. The behavior of the simplified formula for mutual inductance is presented in Figure 4. The precision of the approximation by (12) is such that both curves (the analytical and simplified) are indistinguishable. Another advantage of the approximation (12) is that it can be easily implemented in subsequent analytical calculations.

4. Analytical Solution

In this section, the asymptotic solution of the system of Equations (5)–(7) is found in steady state conditions. The smallness of the mutual inductance (see parameters ${\alpha }_1$ and ${\alpha }_2$) makes possible the implementation of the perturbation theory (PT). This fact is emphasized using a small dimensionless parameter $\epsilon \ll 1$, which is included in the set of Equations (5)–(7) as follows:

$$\begin{array}{ccc}& & l_1\frac{d{i}_1}{dt}+R_1{i}_1=V\left(t\right)-\epsilon M\frac{d{i}_2}{dt}-\epsilon i_2\frac{dM}{dt},\hfill \end{array}$$

(13)

$$\begin{array}{ccc}& & l_2\frac{d{i}_2}{dt}+R_2{i}_2=-\epsilon M\frac{d{i}_1}{dt}-\epsilon i_1\frac{dM}{dt},\hfill \end{array}$$

(14)

$$\begin{array}{ccc}& & \ddot{z}=\epsilon \frac{i_1{i}_2}{m}\frac{dM}{dz}-\frac{\beta }{m}\dot{z}-{\epsilon }^{2}g,\hfill \end{array}$$

(15)

where the multiplier ${\epsilon }^2$ next to the last term on the right-hand side of Equation (15) has been placed in order to provide the existence of a bounded steady state solution to the equation system (13)–(15).

A harmonic input voltage in the primary circuit is assumed:

$$\begin{array}{c}\hfill V\left(t\right)=\frac{1}{2}\left(V{e}^{j\omega t}+V^{*}{e}^{-j\omega t}\right)\end{array}$$

(16)

where V and $V^{*}$ are complex amplitudes ($V^{*}$ is the complex conjugate of V). Notice that Equation (16) gives a cosine-form input voltage. Indeed, complex amplitude V can be represented in the form: $V=\left|V\right|e^{j\varphi }$, where $\varphi$ is the phase. Then, $V^{*}=\left|V\right|e^{-j\varphi }$. Substitution of these expressions into Equation (16) yields: $V\left(t\right)=\left|V\right|\cos (\omega t+\varphi )$. Although the input voltage is of the cosine form, higher harmonics in the currents $i_1\left(t\right)$ and $i_2\left(t\right)$ and the height $z\left(t\right)$ appear due to the fact that Equations (13)–(15) are nonlinear. According to PT, they can be expanded into a power series with respect to the small parameter $\epsilon$ as follows:

$$\begin{array}{ccc}& & i_k\left(t\right)=i_{k,0}\left(t\right)+\epsilon i_{k,1}\left(t\right)+{\epsilon }^2{i}_{k,2}\left(t\right)+\mathcal{O}\left({\epsilon }^3\right),\hfill \end{array}$$

(17)

$$\begin{array}{ccc}& & z\left(t\right)=z_0\left(t\right)+\epsilon z_1\left(t\right)+{\epsilon }^2{z}_2\left(t\right)+\mathcal{O}\left({\epsilon }^3\right),\hfill \end{array}$$

(18)

where $k=1,2$, and $\mathcal{O}\left({\epsilon }^3\right)$ refers to the well-known Landau big-$\mathcal{O}$ notation. $i_{k,n}\left(t\right)$ (where $k=1,2$) and $z_n\left(t\right)$ are the nth order PT contributions to the primary and secondary currents and the ring coordinate (height), respectively.

The substitution of the expansions (17) and (18) into the system (13)–(15) leads to an infinite chain of equations. According to PT, it can be truncated at a certain term. As a result, the zero order of PT yields:

$$\begin{array}{ccc}& & l_1\frac{d{i}_{1,0}}{dt}+R_1{i}_{1,0}=V\left(t\right),\hfill \end{array}$$

(19)

$$\begin{array}{ccc}& & l_2\frac{d{i}_{2,0}}{dt}+R_2{i}_{2,0}=0,\hfill \end{array}$$

(20)

$$\begin{array}{ccc}& & {\ddot{z}}_0=-\frac{\beta }{m}{\dot{z}}_0.\hfill \end{array}$$

(21)

The general solution to Equation (21) is

$$\begin{array}{c}\hfill z_0=C_1+C_2{e}^{-\beta t/m},\end{array}$$

(22)

where $C_1$ and $C_2$ are constants. In can be seen in (22) that the unique steady state solution to Equation (21) is

$$\begin{array}{c}\hfill z_0=\mathrm{const}.\end{array}$$

(23)

This solution represents the equilibrium position of the ring. The explicit form of $z_0$ will be found below for higher orders of PT.

After taking into account the solution (23), the system of Equations (13)–(15) in the first order of PT becomes:

$$\begin{array}{c}\hfill l_1\frac{d{i}_{1,1}}{dt}+R_1{i}_{1,1}=-M\left(z_0\right)\frac{d{i}_{2,0}}{dt},\hfill \end{array}$$

(24)

$$\begin{array}{c}\hfill l_2\frac{d{i}_{2,1}}{dt}+R_2{i}_{2,1}=-M\left(z_0\right)\frac{d{i}_{1,0}}{dt},\hfill \end{array}$$

(25)

$$\begin{array}{c}\hfill {\ddot{z}}_1=\frac{i_{1,0}{i}_{2,0}}{m}{M}^{\prime }\left(z_0\right)-\frac{\beta }{m}{\dot{z}}_1.\hfill \end{array}$$

(26)

Finally, second order of PT leads to:

$$\begin{array}{c}\hfill l_1\frac{d{i}_{1,2}}{dt}+R_1{i}_{1,2}=-M\left(z_0\right)\frac{d{i}_{2,1}}{dt}-M^{\prime }\left(z_0\right)z_1\frac{d{i}_{2,0}}{dt}-M^{\prime }\left(z_0\right)\dot{z_1}{i}_{2,0},\end{array}$$

(27)

$$\begin{array}{c}\hfill l_2\frac{d{i}_{2,2}}{dt}+R_2{i}_{2,2}=-M\left(z_0\right)\frac{d{i}_{1,1}}{dt}-M^{\prime }\left(z_0\right)z_1\frac{d{i}_{1,0}}{dt}-M^{\prime }\left(z_0\right){\dot{z}}_1{i}_{1,0},\end{array}$$

(28)

$$\begin{array}{c}\hfill {\ddot{z}}_2=\frac{1}{m}\left[i_{1,0}{i}_{2,0}{M}^{″}\left(z_0\right)z_1+\left(i_{1,0}{i}_{2,1}+i_{1,1}{i}_{2,0}\right)M^{\prime }\left(z_0\right)\right]-\frac{\beta }{m}{\dot{z}}_2-g.\end{array}$$

(29)

The input voltage (16), in the dynamic equilibrium state, produces periodic currents $i_k\left(t\right)$ and a periodic coordinate $z\left(t\right)$, so that each term of the expansions (17) and (18) can be in turn expanded respectively in Fourier series as follows:

$$\begin{array}{c}\hfill i_{k,n}\left(t\right)=\sum _{p=-\infty }^{\infty }{i}_{k,n}^{\left[p\right]}{e}^{jp\omega t},k=1,2\end{array}$$

(30)

$$\begin{array}{c}\hfill z_n\left(t\right)=\sum _{p=-\infty }^{\infty }{z}_n^{\left[p\right]}{e}^{jp\omega t},\end{array}$$

(31)

where $i_{k,n}^{\left[p\right]}$ and $z_n^{\left[p\right]}$ are complex amplitudes. Equations (17) and (18) should be substituted into Equations (19)–(29), and the amplitudes of their respective harmonics should be equated in both parts of each equation. As a result, we will obtain a system of algebraic Equations that connect Fourier coefficients $i_{k,n}^{\left[p\right]}$ and $z_n^{\left[p\right]}$ with the input voltage amplitude $\left|V\right|$. Substitution of (30) into Equations (19) and (20) yields:

$$\begin{array}{ccc}& & \sum _{p=-\infty }^{\infty }\left(jp\omega l_1+R_1\right)i_{1,0}^{\left[p\right]}{e}^{jp\omega t}=\frac{1}{2}\left(V{e}^{j\omega t}+V^{*}{e}^{-j\omega t}\right),\hfill \end{array}$$

(32)

$$\begin{array}{ccc}& & \sum _{p=-\infty }^{\infty }\left(jp\omega l_2+R_2\right)i_{2,0}^{\left[p\right]}{e}^{jp\omega t}=0.\hfill \end{array}$$

(33)

The solutions to these equations are:

$$\begin{array}{ccc}& & i_{1,0}^{\left[p\right]}=\frac{1}{2}\frac{V{\delta }_{p,1}+V^{*}{\delta }_{p,-1}}{R_1+jp\omega l_1},\hfill \end{array}$$

(34)

$$\begin{array}{ccc}& & i_{2,0}^{\left[p\right]}=0,\hfill \end{array}$$

(35)

where ${\delta }_{p,l}$ is the symbol of Kronecker. After substituting these results into (30), the form of the currents in the time domain in the zero order of PT can be obtained:

$$\begin{array}{ccc}& & i_{1,0}\left(t\right)=\mathbf{Re}\left\{\frac{V{e}^{j\omega t}}{R_1+j\omega l_1}\right\},\hfill \end{array}$$

(36)

$$\begin{array}{ccc}& & i_{2,0}\left(t\right)=0.\hfill \end{array}$$

(37)

Similarly, first and second orders of PT are calculated as:

$$\begin{array}{ccc}& & i_{1,1}\left(t\right)=0,\hfill \end{array}$$

(38)

$$\begin{array}{ccc}& & i_{2,1}\left(t\right)=\mathbf{Im}\left\{\frac{\omega M\left(z_0\right)V{e}^{j\omega t}}{(R_1+j\omega l_1)(R_2+j\omega l_2)}\right\},\hfill \end{array}$$

(39)

$$\begin{array}{ccc}& & z_1=\mathrm{const},\hfill \end{array}$$

(40)

$$\begin{array}{ccc}& & i_{1,2}\left(t\right)=-\mathbf{Re}\left\{\frac{{\omega }^2{M}^2\left(z_0\right)V{e}^{j\omega t}}{{(R_1+j\omega l_1)}^2(R_2+j\omega l_2)}\right\},\hfill \end{array}$$

(41)

$$\begin{array}{ccc}& & i_{2,2}\left(t\right)=z_1\mathbf{Im}\left\{\frac{\omega M^{\prime }\left(z_0\right)V{e}^{j\omega t}}{(R_1+j\omega l_1)(R_2+j\omega l_2)}\right\}.\hfill \end{array}$$

(42)

Substitution of (30) into Equation (29) with solutions (36)–(42) taken into account results in the following equation after some simplifications:

$$\begin{array}{c}\hfill {\ddot{z}}_2=-\frac{j\omega M\left(z_0\right)M^{\prime }\left(z_0\right)}{4m}\{\frac{V^2{e}^{j2\omega t}}{{(R_1+j\omega l_1)}^2(R_2+j\omega l_2)}\hfill \\ \hfill -\frac{{V^{*}}^2{e}^{-j2\omega t}}{{(R_1-j\omega l_1)}^2(R_2-j\omega l_2)}-\frac{2j\omega l_2{\left|V\right|}^2}{(R_1^2+{\omega }^2{l}_1^2)(R_2^2+{\omega }^2{l}_2^2)}\}-\frac{\beta }{m}{\dot{z}}_2-g\hfill \end{array}$$

(43)

Equation (43) has a solution that linearly increases with time. To avoid this instability, it is necessary to equate the constant part on the right-hand side of (43) to zero. As a result, the following equation can be obtained:

$$\begin{array}{c}\hfill -\frac{{\omega }^2{l}_{2}M\left(z_0\right)M^{\prime }\left(z_0\right){\left|V\right|}^2}{2m(R_1^2+{\omega }^2{l}_1^2)(R_2^2+{\omega }^2{l}_2^2)}-g=0.\end{array}$$

(44)

This equation associates the equilibrium height $z_0$ with the input voltage amplitude $\left|V\right|$. In order to establish the dependence of the ring position on the input voltage amplitude, the simplified formula for mutual inductance (12) should be substituted into Equation (44). Therefore, we come across with the following transcendental equation:

$$\frac{{\omega }^2{l}_2{\left|V\right|}^2}{2m(R_1^2+{\omega }^2{l}_1^2)(R_2^2+{\omega }^2{l}_2^2)}\left({\alpha }_1{e}^{-{\beta }_1{z}_0}+{\alpha }_2{e}^{-{\beta }_2{z}_0}\right)\times \left({\alpha }_1{\beta }_1{e}^{-{\beta }_1{z}_0}+{\alpha }_2{\beta }_2{e}^{-{\beta }_2{z}_0}\right)=g,$$

(45)

which can be numerically solved with respect to the ring position.

After taking into account (44), the solution to Equation (43) can be found in the form:

$$\begin{array}{ccc}\hfill z_2\left(t\right)& =& -\frac{M\left(z_0\right)M^{\prime }\left(z_0\right)}{4}\mathbf{Re}\left\{\frac{V^2{e}^{2j\omega t}}{{(R_1+j\omega l_1)}^2(R_2+j\omega l_2)(\beta +2jm\omega )}\right\}.\hfill \end{array}$$

(46)

The approximate time-behaviors of the currents $i_1\left(t\right)$, $i_2\left(t\right)$ and the ring height $z\left(t\right)$ are finally as follows:

$$\begin{array}{ccc}\hfill i_1\left(t\right)& =& \mathbf{Re}\left\{\frac{V{e}^{j\omega t}}{R_1+j\omega l_1}\left(1-\frac{{\omega }^2{M}^2\left(z_0\right)}{(R_1+j\omega l_1)(R_2+j\omega l_2)}\right)\right\},\hfill \end{array}$$

(47)

$$\begin{array}{ccc}\hfill i_2\left(t\right)& =& \mathbf{Im}\left\{\frac{\omega M\left(z_0\right)V{e}^{j\omega t}}{(R_1+j\omega l_1)(R_2+j\omega l_2)}\right\},\hfill \end{array}$$

(48)

$$\begin{array}{ccc}\hfill z\left(t\right)& =& z_0-\frac{M\left(z_0\right)M^{\prime }\left(z_0\right)}{4}\mathbf{Re}\left\{\frac{V^2{e}^{2j\omega t}}{{(R_1+j\omega l_1)}^2(R_2+j\omega l_2)(\beta +2jm\omega )}\right\}.\hfill \end{array}$$

(49)

5. Simulations and Comparison with Experimental Results

5.1. Equilibrium Position of the Ring

The experiment was designed to validate the dynamic model (5)–(7) and the obtained equations (12) and (52) for mutual inductance. The dynamic simulations were performed by solving the system of Equations (5)–(7) employing the standard Runge–Kutta method. The values of the parameters used in the simulations are given in

Table 1. Parameters of the model.

Parameter	Symbol	Value
Coil resistance	$R_1$	$10.200$$\Omega$
Coil inductance	$l_1$	$36.303$$\mathrm{m}$$\mathrm{H}$
Ring resistance	$R_2$	$61.277$ μΩ
Ring inductance	$l_2$	$58.200$ μH
Ring mass	m	$0.053$$\mathrm{k}$$\mathrm{g}$
Friction coefficient	$\beta$	$0.404$$\mathrm{N}$$\mathrm{s}$m−1
Voltage amplitude	V	120–180 V
Frequency	f	60 Hz
Parameters of the	${\alpha }_1$	$464.000$ μH
mutual inductance	${\alpha }_2$	$88.000$ μH
(analytical formula)	${\beta }_1$	$4.776$ m−1
	${\beta }_2$	$22.746$ m−1
Parameters of the	$\gamma$	$1377.91$$\mathrm{V}$m−1
mutual inductance	$\delta$	$41.04$ V
(semi empirical)

.

The parameters of the model were obtained by means of laboratory measurements and manufacturing data [16].

In Figure 5 a schematic diagram illustrating the experimental setup used to measure the aluminum ring height is shown. A camera was installed in front of the Thomson’s ring system at a distance $T>z$, registering the ring position in a slow-motion mode 120 fps (frames per second).

In Figure 6, the experimental and analytical behavior of the ring height vs. the RMS voltage is presented. It can be observed in Figure 6 that, in general, both curves are correlated. Nevertheless, a complete match was not attained. This may have been due to the approximation of an infinitely long core that was adopted to obtain the Equations (9) and (11). In the real experimental setup, as shown in Figure 1, we dealt with a core of a finite height. In this real case, a less number of magnetic lines were concentrated in the core due to dispersion through the upper cross-section of the core. As a result, the magnetic flux through the ring was less than in the case of the infinitely long core. Therefore, the inductance was less than that theoretically calculated in (9) and (11).

5.2. Semi-Empirical Formula for the Mutual Inductance

The analytical calculation of the realistic inductance corresponding to the experimental setup (see Figure 1) is a complicated mathematical task due to its nontrivial geometry. Nevertheless, a different strategy can be adopted in order to estimate the mutual inductance between the coil and ring. In Figure 6, it can be seen that the experimental dependence of the ring position on the RMS voltage is almost linear within the interval of variation of the ring position. This fact can be used to determine the semi-empirical form of mutual inductance. First, assuming the linear behavior of the RMS voltage shown in Figure 6, we can write:

$$\begin{array}{c}\hfill V_{\mathrm{rms}}=\gamma z+\delta ,\end{array}$$

(50)

where the values of $\gamma$ and $\delta$ can be estimated using the least squares method. The results of estimation are shown in Table 1. The relative $L_2$-error of this estimation calculated according to the formula:

$$\begin{array}{c}\hfill \mathrm{Error}=\sqrt{\frac{{\sum }_{k=1}^{10}{\left(V_k-\gamma z_k-\delta \right)}^2}{{\sum }_{k=1}^{10}{V}_k^2}}\times 100\%,\hfill \end{array}$$

is about $2.87\%$ which is considered good enough. The resulting empiric RMS voltage behavior is shown in Figure 7, where a good match between both the experimental and empiric curves can be observed.

Consequently, Equation (44) can be represented in the form of the following differential equation:

$$\begin{array}{c}\hfill \frac{d{M}^2\left(z\right)}{dz}=-\frac{2mg(R_1^2+{\omega }^2{l}_1^2)(R_2^2+{\omega }^2{l}_2^2)}{{\omega }^2{l}_2{V}_{\mathrm{rms}}^2\left(z\right)},\end{array}$$

(51)

where $V_{\mathrm{rms}}=\left|V\right|/\sqrt{2}$ is the RMS voltage. By substituting the empirical formula (50) into Equation (51), this equation can be easily solved with respect to the mutual inductance by integrating with the boundary condition $M\left(z\right)\to 0$ as $z\to \infty$. The result of integration of Equation (51) is the following:

$$\begin{array}{c}\hfill M\left(z\right)=\sqrt{\frac{2mg(R_1^2+{\omega }^2{l}_1^2)(R_2^2+{\omega }^2{l}_2^2)}{{\omega }^2{l}_2\gamma \left(\gamma \right|z|+\delta )}}.\end{array}$$

(52)

The obtained semi-empirical mutual inductance (solid line) and its analytical counterpart (dashed line) as functions of the distance z are shown in Figure 8 (notice that it is not possible to implement a sensor to directly measure the mutual inductance). It can be observed that, as expected, the semi-empirical inductance obtained based on the experiments is lower than that provided by analytical calculations.

Of course, in Equation (50), instead of a linear function, we could use a polynomial of order N:

$$\begin{array}{c}\hfill V\left(z\right)=\sum _{k=0}^N{\delta }_k{z}^k,\end{array}$$

(53)

where N is the number of experimental points in Figure 6—thereby obtaining voltage that better fits with the measured points. However, our aim was not only to reach a better fit in Figure 4, but also to solve Equation (51) analytically and obtain a suitable formula for calculating the mutual inductance. On the one hand, substitution of different values of N into Equation (53) shows that the higher the number N, the more complicated the solution to Equation (51). On the other hand, the use of polynomials to approximate functions is valid only within certain intervals. Out of these intervals, the divergence is higher by as much as the order N of the polynomial (53). Therefore, the implementation of polynomials in finding the mutual impedance (52) could lead to divergence for higher distances z. At the same time, Equation (6) provides an acceptable approximation; it results in Equation (51) being easy to integrate; the solution to Equation (51) (see Equation (52)) is an elementary function, which provides correct asymptotic behavior of the mutual inductance at long distances.

5.3. Dynamic Behavior of the Ring

The dynamic behavior of the ring in the transient case was obtained from numerical simulations of Equations (5)–(7) for two different values of the input voltage amplitude, (a) 180 and (b) 120 V (Figure 9), and two different expressions (analytical and semi-empirical) for the mutual inductance $M\left(z\right)$. The dashed and solid lines correspond to simulation based on Equations (5)–(7), whereas the line with points shows the experimental behavior of the ring height. The numerical results corresponding to the dashed-line (“Analytical”) were obtained using approximation (12) for the mutual inductance, based on the idealized model of an infinitely long core and analytical expressions (9) and (11). The numerical results represented by the solid-line (“Semi-empirical”) were computed based on the improved semi-empirical formula (52). The semi-empirical formula (52) approximates the mutual impedance much better than the idealized expression (12), which, in turn, has reflected in the fact that the “Semi-empirical” continuous lines better fit the experimental waveforms. Moreover, simulated results using (52) led to equilibrium state heights 0.0794 and 0.147 m for 120 and 180 V, respectively, which is an excellent match with the analytical results obtained from (44). The latter means that the obtained analytical solutions represent the steady-state behavior of the aluminum ring with a good precision.

According to the solutions (47)–(49), the interaction force between the coil and aluminum ring, averaged over an oscillation period, has the form:

$$\begin{array}{c}\hfill f_{\mathrm{AS}}=\overline{i_1{i}_2{M}^{\prime }\left(z\right)}=-\frac{{\omega }^2{l}_{2}M\left(z_0\right)M^{\prime }\left(z_0\right){\left|V\right|}^2}{2m(R_1^2+{\omega }^2{l}_1^2)(R_2^2+{\omega }^2{l}_2^2)},\end{array}$$

which according to (44) equals the gravitational force:

$$\begin{array}{c}\hfill f_{\mathrm{AS}}=mg.\end{array}$$

(54)

Similarly, the RMS currents ${i_1}_{\mathrm{rms}}$ and ${i_2}_{\mathrm{rms}}$ can be obtained from (47) and (48) up to the terms of the order ∼$M^2\left(z_0\right)$:

$$\begin{array}{ccc}& & {i_1}_{\mathrm{rms}}=\frac{V_{\mathrm{rms}}}{\sqrt{R_1^2+{\omega }^2{l}_1^2}}\left(1-\frac{{\omega }^2{M}^2\left(z_0\right)(R_1{R}_2-{\omega }^2{l}_1{l}_2)}{(R_1^2+{\omega }^2{l}_1^2)(R_2^2+{\omega }^2{l}_2^2)}\right),\hfill \end{array}$$

(55)

$$\begin{array}{ccc}& & {i_2}_{\mathrm{rms}}=\frac{M\left(z_0\right)V_{\mathrm{rms}}}{\sqrt{(R_1^2+{\omega }^2{l}_1^2)(R_2^2+{\omega }^2{l}_2^2)}}.\hfill \end{array}$$

(56)

6. Conclusions

In this work, an analytical method of solving the dynamical model of a MagLev system was proposed. The method permits finding the dynamical equilibrium of the Thomson’s jumping ring system. The precision of the calculation of the equilibrium state is closely related to the accuracy of the formula for estimating the mutual inductance between the coil and ring. A widely-known formula for the mutual inductance between two coils placed on an infinitely long cylindrical core was corrected and adapted for analytical calculation of the inductance of the system. However, the geometry considered in calculating the mutual inductance does not faithfully reproduce the considered MagLev system (in the Thomson’s jumping ring system the core is of a finite length). The analytical method of analysis of the considered MagLev system was essentially improved by implementing a new semi-empirical formula for the mutual inductance between the coil and aluminum ring, which was obtained based on measurements and the analytical solution for the equilibrium state. The obtained dynamical model presents with good precision the transient-state dynamics of the system. The resulting analytical solution—and particularly Equation (4)—can be effectively implemented in the dynamical control of the ring height. Consequently, the developed model can be implemented for electromagnetic shock absorbers and suspension systems.