A teeth-discretized electromechanical model of a traveling-wave ultrasonic motor

This paper develops an electromechanical TWUSM (traveling-wave ultrasonic motor) model combining the driving circuit with the motor itself. An equivalent circuit model substituted for the piezoelectric ceramics is designed in the driving circuit model to obtain accurate input currents and powers. Then teeth discretization is implemented in the stator–rotor contact model, which can calculate the interaction forces more accurately. After building the complete model of TWUSM, a typical startup–stopping process is divided into five stages by evaluating the changes in contact status and driving forces. Finally, the fitness of transient responses of the rotor speed increasing from 55 % to 89 % shows that the proposed model fits better than the one without teeth discretization, and the experimental tests under various driving parameters verify the effectiveness of the model.

1 Introduction

The TWUSM (traveling-wave ultrasonic motor) is the piezoelectric actuator which is excited near the resonant frequencies. The polarized piezoelectric ceramics of TWUSM are usually actuated by two alternating singular voltages. The stator particles change their displacements in space elliptical motion (Liu et al., 2019; Chen et al., 2018; Shi et al., 2018) with the deformations of piezoelectric elements. Several teeth are distributed along the circumferential direction to amplify the driving effect and deformed in accordance with the traveling wave (Renteria-Marquez et al., 2018). Therefore these motors present many significant merits compared to electromagnetic motors, for instance, rapid response, simple structure, and the capability of miniaturization (Zhang et al., 2016).

Because of the combination of piezoelectric actuation and friction drive, the TWUSM model attracts many researchers' attention. In general, their models can be divided into three types. First and most typical is the model stemming from Hagood and McFarland (1995), who assume the vibrating stator to be a two-freedom lumped spring-mass-damping system. However, the teeth discretization is ignored by only assuming the contact model to be continuous springs covering the ring area. The second type is proposed by Giraud et al. (2004), who imitate the d–q decomposition from the three-phase alternating current motor. Similarly, the ultrasonic motor is decomposed in the d–q coordinate, where d means the modal value and q represents the torque value. Jing (2015) develops Giraud's model by evaluating the modal vibration trajectory of the stator. However, the modal cannot be controlled accurately, like the external input such as the current or torque within the electromagnet motor. The third type is the semi-analytical model developed by Hagedorn et al. (1993) and Chen and Zhao (2005): they divide the stator into several parts according to the respective shape functions. Similarly, Bolborici et al. (2014) and Arturo (2016) model the stator and the rotor through the finite-volume-method model. Though the detailed model can analyze more microscopic details, it consumes vast computation sources and cannot achieve convergence all the time.

Besides these deficiencies, the above researchers simplify the model by idealizing the input signals as the ideal sinusoidal signals instead of deriving them from the driving circuit. The simplification not only deviates from reality, but also fails to detect the real-time currents. In this paper, a hybrid mechatronic model is proposed by combining the electrical system with the mechanical system. Besides, a more straightforward discretization strategy is adopted in the hybrid model to solve the problems of insufficient or excessive discretization, which can not only guarantee the accuracy, but also reduce the simulation nodes and calculation sources.

This paper is organized as follows. In Sect. 2, the principle of the ultrasonic motor is introduced. Then a hybrid model including the electrical circuit and the ultrasonic motor is built in the third section. In Sect. 4, the transient response and the contact status are studied under different torques to explore the microcosmic law. In Sect. 5, the integrated test system is built, and the simulation results under different parameters are verified from the comparisons of the experimental results.

2 Structure of the TWUSM system

Figure 1 illustrates the working mechanism of the TWUSM. Piezoelectric ceramics are actuated by the two-phase sinusoidal voltages (U_A and U_B) determined by the amplitude (U_m), frequency (f) and phase difference (α), as shown in Fig. 1a. It is evident that the two-channel input currents I_A and I_B can be informed by linking to the electrical network of the piezoelectric ceramics. Then the stator vibrates with the amplitude ξ during the transformation from the electrical energy to the mechanical energy, just like Fig. 1b. The modal responses of phase A and phase B are characterized as q_A and q_B, respectively. The stator's circumferential rotation propelled by the traveling wave drives the rotor via the friction interaction with the rotor. In Fig. 1c, F_pre represents the preload force, and F_z means the vertical force acting on the rotor. As shown in Fig. 1d, the output torque T_out is producing to overcome the applied torque T_load. Moreover, the whole rotor's mass and inertial are M_r and J_r, respectively.

Figure 1 The block diagram of the TWUSM: (a) the driver; (b) the vibration and traveling wave of the stator; (c) the stator–rotor contact; (d) the rotor with the externally applied load.

3 The analytical model of the TWUSM system

The TWUSM prototype investigated in this paper is USR60-S3 (Shinsei Corp. Ltd, Japan), which works in the ninth vibration mode (N=9). So the natural vibration bending mode of the ring plate is denoted by B₀₉. In USR60-S3, 90 teeth are distributed along the circumferential direction of the stator to improve overall driving capacities (Lin et al., 2002).

3.1 Modeling of the driving circuit

The driving circuit is supposed to generate pure sinusoidal waves for each fragment of the piezoelectric ceramics. Figure 2 illustrates the framework of the driving circuit actuated by three parameters (U_m, F_s, α). The original sinusoid waves (U_A0, U_B0) are achieved by the two-channel waveform generators. Subsequently, the amplifying circuits and the transformers are assigned to produce the driving voltages. In terms of the capacitive characteristics of the piezoelectric segments, two inductors are placed in parallel to reduce the reactive power dissipation. Furthermore, in order to obtain the input voltages and currents, the circuit equivalent model (Fig. 3) is adopted to substitute for the piezoelectric ceramics (Mojallali et al., 2007).

Figure 2 The framework of the driving circuit.

Figure 3 Equivalent circuit of the piezoelectric ceramics (one phase).

As shown in Fig. 3, C_p is the clamping capacitor of piezoelectric ceramics, whereby R_m, C_m and L_m represent equivalent resistance, equivalent capacitance and equivalent inductance, respectively. R_d means the resistance loss during the energy dissipation (Lu et al., 2011, 2020). The right element is the equivalent voltage ϕ_v, which describes the modal force derived from the stator–rotor contact model discussed in the next section. Finally, the current source I_T represents the variation with the external load. It can be calculated from the proportion of the speed reduction caused by the applied torque. Its equation can be depicted as

$$\begin{array}{}\text{(1)}& I_T={\displaystyle \frac{\sqrt{\mathrm{2}}\mathit{\sigma }{k}_{\mathrm{c}}{R}_{\mathrm{o}}^{\mathrm{2}}{T}_{\mathrm{load}}}{N(h_{\mathrm{s}}+h_{\mathrm{p}})}}.\end{array}$$

Here, R_o is the middle radius of the stator; h_s and h_p are the thickness of the stator and the piezoelectric layer. k_c is the mechatronic coupling coefficient which describes the speed drop versus the applied torque (T_load). When the maximum voltages and currents are U_am, U_bm, I_am and I_bm, the active input power of the two-phase piezoelectric ceramics can be calculated as

$$\begin{array}{}\text{(2)}& {\displaystyle }{P}_{\mathrm{A}}=\mathrm{0.5}{U}_{\mathrm{am}}{I}_{\mathrm{am}}\cos \left({\mathit{\varphi }}_{\mathrm{a}}\right),\text{(3)}& {\displaystyle }{P}_{\mathrm{B}}=\mathrm{0.5}{U}_{\mathrm{bm}}{I}_{\mathrm{bm}}\cos \left({\mathit{\varphi }}_{\mathrm{b}}\right).\end{array}$$

3.2 Modeling of the ultrasonic motor

3.2.1 Vibration model of the stator with a piezoelectric ring

In order to simplify the calculation, the vibration system (the piezoelectric elements and the stator) can be characterized as a 2-degree-of-freedom spring-mass-damping system. The modal coordinates (q_A, q_B) and the dynamic functions can be depicted as

$$\begin{array}{}\text{(4)}& {\displaystyle }q=\left[\begin{array}{l}{q}_{\mathrm{A}}\left(t\right)\\ q_{\mathrm{B}}\left(t\right)\end{array}\right]=\left[\begin{array}{l}\mathit{\xi }\cos \left(wt\right)\\ \mathit{\xi }\sin \left(wt\right)\end{array}\right],\text{(5)}& {\displaystyle }\begin{array}{rl}& m_{\mathrm{o}}{\ddot{q}}_{\mathrm{A}}+d_{\mathrm{o}}{\dot{q}}_{\mathrm{A}}+k_{\mathrm{o}}{q}_{\mathrm{A}}=k_{\mathrm{c}}{U}_{\mathrm{A}}+\mathit{\epsilon }{k}_{\mathrm{c}}{U}_{\mathrm{B}}+F_{\mathrm{d}\mathrm{1}},\\ & m_{\mathrm{o}}{\ddot{q}}_{\mathrm{B}}+d_{\mathrm{o}}{\dot{q}}_{\mathrm{B}}+k_{\mathrm{o}}{q}_{\mathrm{B}}=k_{\mathrm{c}}{U}_{\mathrm{B}}+\mathit{\epsilon }{k}_{\mathrm{c}}{U}_{\mathrm{A}}+F_{\mathrm{d}\mathrm{2}},\end{array}\end{array}$$

where w=2πF_s denotes the angular driving frequency, m_o, d_o and k_o represent modal mass, modal damping and modal stiffness, respectively, and ε is the imbalance coefficient between two-phase signals. F_d1 and F_d2 are the respective modal forces in both the tangential direction and the vertical direction. Besides, if λ is the wavelength of the traveling wave and $k=\mathrm{2}\mathit{\pi }/\mathit{\lambda }$, $\cos k\stackrel{⌢}{x}$ and $\sin k\stackrel{⌢}{x}$ are the eigenmode functions. When Φ(r) is the radial variance of this mode shape, the orthogonal modal shape function can be represented as

$$\begin{array}{}\text{(6)}& {\mathrm{\Phi }}_{\mathrm{w}}=\left[\begin{array}{c}{\mathrm{\Phi }}_{\mathrm{wA}}(r,\stackrel{⌢}{x})\\ {\mathrm{\Phi }}_{\mathrm{wB}}(r,\stackrel{⌢}{x})\end{array}\right]=\left[\begin{array}{c}\mathrm{\Phi }\left(r\right)\cos k\stackrel{⌢}{x}\\ \mathrm{\Phi }\left(r\right)\sin k\stackrel{⌢}{x}\end{array}\right].\end{array}$$

Thus, the vibration response of the stator stuck with the piezoceramics can be obtained by multiplying Eqs. (4) and (6); the result is depicted as

$$\begin{array}{}\text{(7)}& w_{\mathrm{0}}(r,\stackrel{⌢}{x},t)={\mathrm{\Phi }}_{\mathrm{w}}^{T}q=\mathrm{\Phi }\left(r\right)\mathit{\xi }\cos (k\stackrel{⌢}{x}-\mathit{\omega }t).\end{array}$$

Finally, the tangential speed of the surface point can be derived. If h is the half-thickness of the stator and the middle radius of the stator is r_p, whose radial variance is defined as ψ_av=Φ(r_p), then the stator speed can be expressed as

$$\begin{array}{}\text{(8)}& v_{\mathrm{s}}(r_{\mathrm{p}},\stackrel{⌢}{x},t)=-h{\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}{\displaystyle \frac{\partial w_{\mathrm{0}}(\stackrel{⌢}{x},t)}{\partial \stackrel{⌢}{x}=-h\mathit{\xi }k{\mathit{\psi }}_{\mathrm{av}}\cos (k\stackrel{⌢}{x}-wt)}}.\end{array}$$

3.2.2 Contact model with teeth discretization

Figure 4 illustrates the contact schematics of the motor with teeth discretization, and the stator–rotor contact model is built on the assumption that the rotor is rigid over all the contact areas, whereas the friction layer can be modeled as a series of linear springs. The yellow points mean the contact borders with the zone [$-x_{\mathrm{0}},x_{\mathrm{0}}$], and the demarcation points between the traction zone and the braking zone are −x₁ and x₁.

Figure 4 The contact schematic of the TWUSM.

According to Fig. 4 and Eq. (9), the coordinates of contact border x₀ and stick point x₁ yield

$$\begin{array}{}\text{(9)}& {\displaystyle }{x}_{\mathrm{0}}={\displaystyle \frac{\mathrm{1}}{k}}\arccos \left({\displaystyle \frac{h-z\left(t\right)}{{\mathit{\psi }}_{\mathrm{av}}\mathit{\xi }}}\right),\text{(10)}& {\displaystyle }{x}_{\mathrm{1}}={\displaystyle \frac{\mathrm{1}}{k}}\arccos \left({\displaystyle \frac{{\dot{\mathit{\theta }}}_{\mathrm{r}}{R}_{\mathrm{0}}^{\mathrm{2}}}{kh{\mathit{\psi }}_{\mathrm{av}}\mathit{\xi }\mathit{\omega }}}\right),\end{array}$$

where z(t) is the real-time vertical displacement of the rotor and ${\dot{\mathit{\theta }}}_{\mathrm{r}}$ means the rotor speed. When the joint stiffness of the dispersed springs is K_f, the unit pressure distributing from the contact interface to the stator can be given by

$$\begin{array}{}\text{(11)}& f\left(x\right)=K_{\mathrm{f}}\mathrm{\Delta }x=K_{\mathrm{f}}{\mathit{\psi }}_{\mathrm{av}}\mathit{\xi }\left[\cos \right(kx)-\cos (k{x}_{\mathrm{0}}\left)\right].\end{array}$$

As shown in Fig. 5, every tooth can be simplified as a cuboid whose length and width are R_c and f, respectively. Since USR60 operates in ninth mode with 90 teeth, 10 teeth scatter in every wavelength. The coordinate plane OYZ is built on the central axis of the middle tooth, and the x axis is assigned along the traveling wave direction.

Figure 5 The teeth discretization implemented in the stator.

In succession, the left teeth edges are recorded as A_j, while the right ones are assigned as B_j. As shown in Fig. 5, the code numbers are named for every tooth, and only half of No. 1 and No. 11 are included in one wavelength. We define the gap width as d. Therefore every tooth position can be determined by its edge, whose horizontal coordinate and vertical coordinate are depicted as

$$\begin{array}{}\text{(12)}& {\displaystyle }\begin{array}{rl}& {\stackrel{⌢}{X}}_{\mathrm{A}j}={\displaystyle \frac{(\mathrm{2}n-\mathrm{13})f}{\mathrm{2}N}}+{\displaystyle \frac{(n-\mathrm{6})d}{N}},\\ & {\stackrel{⌢}{X}}_{\mathrm{B}j}={\displaystyle \frac{(\mathrm{2}n-\mathrm{11})f}{\mathrm{2}N}}+{\displaystyle \frac{(n-\mathrm{6})d}{N}},\end{array}\text{(13)}& {\displaystyle }\begin{array}{rl}& Z_{\mathrm{A}j}=h+{\mathit{\psi }}_{\mathrm{av}}\mathit{\xi }hk\cos (wt-k{\stackrel{⌢}{X}}_{\mathrm{A}j}),\\ & {\mathrm{Z}}_{\mathrm{B}j}=h+{\mathit{\psi }}_{\mathrm{av}}\mathit{\xi }hk\cos (wt-k{\stackrel{⌢}{X}}_{\mathrm{B}j}).\end{array}\end{array}$$

Therefore, there are four cases when comparing the vertical coordinates with the stator half-thickness h. For convenience, the real contact boundaries are symbolized as s₁, s₂ with the following four cases.

• a.

  If Z_Aj<h and Z_Bj<h, both edges of the tooth are apart from the contact area, and this tooth makes no contributions to the force generation.

• b.

  If Z_Aj>h and Z_Bj>h, both the edges of the tooth are located in the contact area; therefore, $s_{\mathrm{1}j}={\stackrel{⌢}{X}}_{\mathrm{A}j}$, $s_{\mathrm{2}j}={\stackrel{⌢}{X}}_{\mathrm{B}j}$.

• c.

  If Z_Aj>h and Z_Bj<h, the left edge of the tooth lies in the contact zone, and the intersection between the contact area and the tooth is just the right contact border; then, $s_{\mathrm{1}j}={\stackrel{⌢}{X}}_{\mathrm{A}j}$, s_2j=x₀.

• d.

  If Z_Aj<h and Z_Bj>h, the right edge of the tooth is located in the contact region; similarly, $s_{\mathrm{1}j}=-x_{\mathrm{0}}$, $s_{\mathrm{2}j}={\stackrel{⌢}{X}}_{\mathrm{B}j}$.

However, not all the particles of the contact region serve as the valid driving points. Whether the particles play the driving role needs the following detailed further discussions of stick–slip regional distribution, which also can be divided into four cases via the comparison between the contact borders and stick points.

• a.

  If $s_{\mathrm{1}j}<-x_{\mathrm{1}}<x_{\mathrm{1}}<s_{\mathrm{2}j}$, both the stick points are in the contact area, the driving zone equals $[-x_{\mathrm{1}},x_{\mathrm{1}}]$, and the blocking zones are $[s_{\mathrm{1}j},-x_{\mathrm{1}}]$ and [x₁, s_2j].

• b.

  If $-x_{\mathrm{1}}<s_{\mathrm{1}j}<x_{\mathrm{1}}<s_{\mathrm{2}j}$, the driving zone is [s_1j, x₁], while the blocking zone is [x₁, s_2j].

• c.

  If $s_{\mathrm{1}j}<-x_{\mathrm{1}}<s_{\mathrm{2}j}<x_{\mathrm{1}}$, the driving zone is $[-x_{\mathrm{1}},s_{\mathrm{2}j}]$, while the blocking zone is $[s_{\mathrm{1}j},-x_{\mathrm{1}}]$.

• d.

  If $-x_{\mathrm{1}}<s_{\mathrm{1}j}<s_{\mathrm{2}j}<x_{\mathrm{1}}$, all the particles share the driving effects.

Based on the above classifications, the interaction forces between the stator and the rotor become more accurate. First of all, the friction force F_T can be derived from the integral operation of the whole teeth. Therefore F_T(j) for every tooth can be expressed as Eqs. (15) and (16):

$$\begin{array}{}\text{(14)}& {\displaystyle }{F}_T\left(j\right)=-\mathit{\mu }\underset{R_e}{\overset{R_{\mathrm{2}}}{\int }}\underset{s_{\mathrm{1}j}}{\overset{s_{\mathrm{2}j}}{\int }}\mathrm{sgn}(\left|v_{\mathrm{s}}\right|-\left|v_{\mathrm{r}}\right|)f\left(x\right)\mathrm{d}x\mathrm{d}y,\text{(15)}& {\displaystyle }{F}_T\left(j\right)=\left\{\begin{array}{l}\mathrm{2}\mathrm{\Gamma }\left(x_{\mathrm{1}}\right)-\mathrm{2}\mathrm{\Gamma }(-x_{\mathrm{1}})+\mathrm{\Gamma }\left(s_{\mathrm{1}j}\right)-\mathrm{\Gamma }\left(s_{\mathrm{2}j}\right)\\ s_{\mathrm{1}j}<-x_{\mathrm{1}}<x_{\mathrm{1}}<s_{\mathrm{2}j}\\ \mathrm{2}\mathrm{\Gamma }\left(x_{\mathrm{1}}\right)-\mathrm{\Gamma }\left(s_{\mathrm{1}j}\right)-\mathrm{\Gamma }\left(s_{\mathrm{2}j}\right)\\ -x_{\mathrm{1}}<s_{\mathrm{1}j}<x_{\mathrm{1}}<s_{\mathrm{2}j}\\ \mathrm{2}\mathrm{\Gamma }\left(x_{\mathrm{1}}\right)-\mathrm{\Gamma }\left(s_{\mathrm{1}j}\right)-\mathrm{\Gamma }\left(s_{\mathrm{2}j}\right)\\ s_{\mathrm{1}j}<-x_{\mathrm{1}}<s_{\mathrm{2}j}<x_{\mathrm{1}}\\ \mathrm{\Gamma }\left(s_{\mathrm{2}j}\right)-\mathrm{\Gamma }\left(s_{\mathrm{1}j}\right)\\ -x_{\mathrm{1}}<s_{\mathrm{1}j}<s_{\mathrm{2}j}<x_{\mathrm{1}}\end{array}\right.,\end{array}$$

with

$$\begin{array}{}\text{(16)}& \mathrm{\Gamma }\left(x\right)={\displaystyle \frac{\mathit{\mu }{K}_{\mathrm{f}}{\mathit{\psi }}_{\mathrm{av}}\mathit{\xi }(R_{\mathrm{2}}-R_e)\left[\sin \right(kx)-kx\cos (k{x}_{\mathrm{0}}\left)\right]}{k}}.\end{array}$$

As shown in Eq. (18), the forces acting on the rotor along the axial direction can be calculated as F_z(j) without any classifications.

$$\begin{array}{}\text{(17)}& F_z\left(j\right)=\underset{R_e}{\overset{R_{\mathrm{2}}}{\int }}\underset{s_{\mathrm{1}j}}{\overset{s_{\mathrm{2}j}}{\int }}f\left(x\right)\mathrm{d}x\mathrm{d}y={\displaystyle \frac{\mathrm{\Gamma }\left(s_{\mathrm{2}j}\right)-\mathrm{\Gamma }\left(s_{\mathrm{1}j}\right)}{\mathit{\mu }}}\end{array}$$

Moreover, the modal forces F_d1 and F_d2 consist of the vertical forcing vectors (F_dn1, F_dn2) and the tangential forcing vectors (F_dt1, F_dt2). The vertical parts can be calculated with the combination of modal matric ${\mathrm{\Phi }}_{\mathrm{w}}^t$ and f(x). It can be expressed as Eqs. (19) and (20).

$$\begin{array}{}\text{(18)}& {\displaystyle }\begin{array}{rl}{F}_{\mathrm{dn}}\left(j\right)& =\left[\begin{array}{c}{F}_{\mathrm{dn}\mathrm{1}}\left(j\right)\\ F_{\mathrm{dn}\mathrm{2}}\left(j\right)\end{array}\right]=-\underset{R_e}{\overset{R_{\mathrm{2}}}{\int }}\underset{s_{\mathrm{1}j}}{\overset{s_{\mathrm{2}j}}{\int }}{\mathrm{\Phi }}_{\mathrm{w}}^{\mathrm{t}}f\left(x\right)\mathrm{d}x\mathrm{d}y\\ & =\mathrm{\Psi }\left(s_{\mathrm{2}j}\right)-\mathrm{\Psi }\left(s_{\mathrm{1}j}\right)\end{array}\text{(19)}& {\displaystyle }\mathrm{\Psi }\left(x\right)={\displaystyle \frac{\mathit{\mu }{K}_{\mathrm{f}}{\mathit{\psi }}_{\mathrm{av}}\mathit{\xi }(R_{\mathrm{2}}-R_e)\left[\sin \right(kx)-kx\cos (k{x}_{\mathrm{0}}\left)\right]}{k}}\end{array}$$

Similarly, the tangential forcing vectors from different areas can be read as Eqs. (21)–(25):

$$\begin{array}{}\text{(20)}& {\displaystyle }{F}_{\mathrm{dt}}\left(j\right)=-\mathit{\mu }h\underset{R_e}{\overset{R_{\mathrm{2}}}{\int }}\underset{s_{\mathrm{1}j}}{\overset{s_{\mathrm{2}j}}{\int }}{\displaystyle \frac{\partial {\mathrm{\Phi }}_w^t}{\partial \stackrel{⌢}{x}f\left(x\right)\mathrm{sgn}(\left|v_{\mathrm{s}}\right|-\left|v_{\mathrm{r}}\right|)}}\mathrm{d}x\mathrm{d}y\text{(21)}& {\displaystyle }\begin{array}{rl}{F}_{\mathrm{dt}}\left(j\right)& =\left[\begin{array}{c}{F}_{\mathrm{dt}\mathrm{1}}\left(j\right)\\ F_{\mathrm{dt}\mathrm{2}}\left(j\right)\end{array}\right]\\ & =\left\{\begin{array}{l}\mathrm{2}M\left(x_{\mathrm{1}}\right)-\mathrm{2}M(-x_{\mathrm{1}})+M\left(s_{\mathrm{1}j}\right)-M\left(s_{\mathrm{2}j}\right)\\ s_{\mathrm{1}j}<-x_{\mathrm{1}}<x_{\mathrm{1}}<s_{\mathrm{2}j}\\ \mathrm{2}M\left(x_{\mathrm{1}}\right)-M\left(s_{\mathrm{1}j}\right)-M\left(s_{\mathrm{2}j}\right)\\ -x_{\mathrm{1}}<s_{\mathrm{1}j}<x_{\mathrm{1}}<s_{\mathrm{2}j}\\ M\left(s_{\mathrm{1}j}\right)+M\left(s_{\mathrm{2}j}\right)-\mathrm{2}M(-x_{\mathrm{1}})\\ s_{\mathrm{1}j}<-x_{\mathrm{1}}<s_{\mathrm{2}j}<x_{\mathrm{1}}\\ M\left(s_{\mathrm{2}j}\right)-M\left(s_{\mathrm{1}j}\right)\\ -x_{\mathrm{1}}<s_{\mathrm{1}j}<s_{\mathrm{2}j}<x_{\mathrm{1}}\end{array}\right.\end{array},\end{array}$$

with

$$\begin{array}{}\text{(22)}& {\displaystyle }M\left(x\right)=\mathit{\mu }n{K}_{\mathrm{f}}h{\mathit{\psi }}_{\mathrm{av}}\mathit{\xi }(R_{\mathrm{2}}-R_e)\left[\begin{array}{c}{q}_{\mathrm{2}}{G}_{\mathrm{1}}-q_{\mathrm{1}}{G}_{\mathrm{2}}\\ -q_{\mathrm{1}}{G}_{\mathrm{1}}-q_{\mathrm{2}}{G}_{\mathrm{2}}\end{array}\right],\text{(23)}& {\displaystyle }{G}_{\mathrm{1}}=\left({\displaystyle \frac{\sin \mathrm{2}kx}{\mathrm{4}}}+{\displaystyle \frac{kx}{\mathrm{2}}}-\sin kx\cos k{x}_{\mathrm{0}}\right),\text{(24)}& {\displaystyle }{G}_{\mathrm{2}}=\left(-{\displaystyle \frac{\cos \mathrm{2}kx}{\mathrm{4}}}+\cos kx\cos k{x}_{\mathrm{0}}\right).\end{array}$$

Finally, since the number of traveling waves is 9, the force results of the whole motor can be derived by summarizing the above calculated results. In the end, the whole force vector can be depicted as

$$\begin{array}{}\text{(25)}& \left[\begin{array}{c}{F}_n\\ F_z\\ F_{\mathrm{dt}\mathrm{1}}\\ F_{\mathrm{dn}\mathrm{1}}\\ F_{\mathrm{dt}\mathrm{2}}\\ F_{\mathrm{dn}\mathrm{2}}\end{array}\right]=\mathrm{9}\sum _{j=\mathrm{1}}^{\mathrm{10}}\left[\begin{array}{c}{F}_n\left(\mathrm{j}\right)\\ F_z\left(\mathrm{j}\right)\\ F_{\mathrm{dt}\mathrm{1}}\left(\mathrm{j}\right)\\ F_{\mathrm{dn}\mathrm{1}}\left(\mathrm{j}\right)\\ F_{\mathrm{dt}\mathrm{2}}\left(\mathrm{j}\right)\\ F_{\mathrm{dn}\mathrm{2}}\left(\mathrm{j}\right)\end{array}\right].\end{array}$$

3.2.3 The whole motor model

In terms of the whole TWUSM, the vertical displacement and the rotational velocity of the rotor can be described in the third and fourth functions in Eq. (27), where d_z represents the vertical damping and d_r means the rotational damping. Finally, the final analytical model can be depicted as

$$\begin{array}{}\text{(26)}& \begin{array}{rl}& m_{\mathrm{o}}{\ddot{q}}_{\mathrm{A}}+d_{\mathrm{o}}{\dot{q}}_{\mathrm{A}}+k_{\mathrm{o}}{q}_{\mathrm{A}}=k_{\mathrm{c}}{U}_{\mathrm{A}}+\mathit{\epsilon }{k}_{\mathrm{c}}{U}_{\mathrm{B}}+F_{\mathrm{dt}\mathrm{1}}+F_{\mathrm{dn}\mathrm{1}},\\ & m_{\mathrm{o}}{\ddot{q}}_{\mathrm{B}}+d_{\mathrm{o}}{\dot{q}}_{\mathrm{B}}+k_{\mathrm{o}}{q}_{\mathrm{B}}=k_{\mathrm{c}}{U}_{\mathrm{B}}+\mathit{\epsilon }{k}_{\mathrm{c}}{U}_{\mathrm{A}}+F_{\mathrm{dt}\mathrm{2}}+F_{\mathrm{dn}\mathrm{2}},\\ & {\mathrm{m}}_{\mathrm{r}}{\ddot{W}}_{\mathrm{r}}+d_z{\dot{W}}_{\mathrm{r}}=F_Z-F_N,\\ & J_{\mathrm{r}}{\ddot{\mathit{\theta }}}_{\mathrm{r}}+d_{\mathrm{r}}{\dot{\mathit{\theta }}}_{\mathrm{r}}=T_{\mathrm{r}}-T_{\mathrm{load}}.\end{array}\end{array}$$

It should be emphasized that the bridges between the circuit model and the motor model are the equivalent electrical voltages (φ_va,φ_vb), which can be expressed as Eq. (28). The efficiency η_e from the input power to the output power can be depicted as Eq. (29), which is used to evaluate the energy utilization factors with variable loads.

$$\begin{array}{}\text{(27)}& {\displaystyle }\begin{array}{rl}& {\mathit{\phi }}_{\mathrm{va}}=F_{\mathrm{d}\mathrm{1}}=F_{\mathrm{dn}\mathrm{1}}+F_{\mathrm{dt}\mathrm{1}}\\ & {\mathit{\phi }}_{\mathrm{vb}}=F_{\mathrm{d}\mathrm{2}}=F_{\mathrm{dn}\mathrm{2}}+F_{\mathrm{dt}\mathrm{2}}\end{array}\text{(28)}& {\displaystyle }{\mathit{\eta }}_{\mathrm{e}}={\displaystyle \frac{P_{\mathrm{out}}}{P_{\mathrm{in}}}}={\displaystyle \frac{\mathrm{2}{T}_{\mathrm{r}}{\dot{\mathit{\theta }}}_{\mathrm{r}}}{P_{\mathrm{A}}+P_{\mathrm{B}}}}\end{array}$$

4 Simulation results

The total model covering the elements in Fig. 1 is built in the Simulink platform, with the parameters listed in Appendix A. In order to obtain more microscopic properties, simulations are implemented on the transient response without load and the contact status under different torques.

4.1 Transient response without load

Figure 6 displays the startup–stopping response of the TWUSM when the input signals are the sinusoid waves containing 800 periods and the amplitude U_m and the frequency F_s are 1 V and 43 kHz, respectively. The simulating signals cover the main electromechanical parameters discussed in Eqs. (10), (11) and (27).

Figure 6 Transient response of the unloaded motor with the driving signals containing 800 periods of activation: (a) the two-phase input voltages, (b) the vertical force applied in the rotor, (c) two-phase input currents, (d) the output torque, (e) the two-phase vibration response, (f) the tangential modal forces, (g) the revolving speed, (h) the vertical modal forces, (i) the rotor angle, and (j) the contact parameters.

There are five stages, including the pre-static stage [T₀, T₁], dynamic-friction and fluctuation stage [T₁, T₂], stabilized stage [T₂, T₃], vibration decay stage [T₃, T₄] and self-locking stage [T₄, T₅]. At the pre-static stage, the force generation between the stator and the rotor cannot overcome the static friction; therefore, the rotor stays motionless. Once the motor comes into the dynamic-friction stage, the vibration amplitude ξ (Fig. 6d) increases, accompanied by the ascending modal forces. The contact status changes from full contact to partial contact. However, the effective driving zone (Fig. 6j) shrinks, and the number of teeth involved in contact or driving also changes periodically, which causes the continuous fluctuation of the input currents (Fig. 6b). When the motor steps into the steady stage, the vibration amplitude and the driving zone become steady. Until the moment that the driving signals are withdrawn, the contact zone extends again and returns to the full-contact status before locking the rotor. At this time, the sustaining driving zone supports the final decaying process. The results also indicate that the proposed comprehensive integration model expands our understanding of the internal microcosmic law of the TWUSM, which is useful for the motor's stepwise position control.

4.2 Contact status with load

Figure 7 proposes the contact status of the teeth in a 40^∘ circumferential range when the torque is 0.5 N m. It can be observed that the stick points are all located inside the contact zones, and only No. 5, No. 6 and No. 7 of the teeth are in full-contact status. The other teeth are out of contact and make no contribution to force generation. Besides, Fig. 7b shows the contact parameters when the torque is increased. The driving zone becomes wider, while the length of the contact zone remains constant, which means more teeth are added to the driving range to overcome the increased torque.

Figure 7 The contact status: (a) the external torque is 0.5 N m; (b) the contact parameters under different torques.

5 Experimental setup and results

5.1 Experimental setup

Figure 8 displays the integrated measurement system that consists of the driving circuit, the mechanical platform, the FPGA board (National Instruments Corp, USA), and the dSPACE1103 control board (dSPACE Corp, Germany). There is an incremental encoder AFS60A (SICK Corp, Germany) which has 65 536 lines. The external load is generated by a torque motor with the maximum torque 1 N m. The dSPACE1103 also generates the input parameters (U_m, F_s, α) which output to the FPGA (field programmable gate array) board to generate the sinusoidal signals (U_A0  U_B0). When the motor is actuated, the currents (I_A1,I_B1) of the piezoelectric ceramics are measured by hall sensors (Zhonghuo Sensing Corp, PRC) and the voltages (U_A1, U_B1) are measured by voltage transformers (Zhonghuo Sensing Corp, PRC). Then, the signals are all processing in the FPGA board and transformed into the amplitude (I_am, I_bm, U_am, U_bm) and the phase differences (φ_a, φ_b). Therefore the input powers (P_A, P_B) can be obtained from Eqs. (2) and (3). Finally, we can get the motor efficiency according to Eq. (29). The experimental platform integrates the function of flexible adjustment and online calculation. Therefore parameter identification and performance evaluation can be realized in high efficiency.

Figure 8 The integrated measurement and control system of USR60-S3.

5.2 Experiment results

5.2.1 The startup response

As shown in Fig. 9, when the amplitude U_m is 1.3 V and the frequency F_s is 41.5 kHz, the startup velocity curves with or without teeth discretization are compared. As can be seen, the results with the teeth discretization method achieve better fitness. The fitness is 55 % if the teeth structure is ignored, and the fitness increases to 89 % when the teeth discretization is implemented. Moreover, the 2 ms delay which occurs in the results without the teeth discretization model disappears with the proposed model. This can be explained by eliminating the driving effect of the tooth gap by refining the interaction forces and the contact status. Furthermore, the contact area which has been limited to the teeth space rather than the whole ring is closer to reality. Furthermore, the startup speed responses actuated by various frequencies, and a certain amplitude, are displayed in Fig. 10a. To gain more detailed observations, the speed curves before 5 ms are displayed in Fig. 10b. All results under different frequencies show good agreements between the numerical results and experimental ones and the differences of the response process. Figure 10b indicates that the motor starts more quickly and obtains less overshoot with the increasing frequency.

Figure 9 The transient characteristics comparing the experiment results and simulation results with or without the teeth discretization.

Figure 10 Comparison between simulation results and experimental ones when the amplitude is 1 V: (a) the step response within 20 ms; (b) the step response within 5 ms.

5.2.2 Speed performances under different driving parameters

Due to the complexity of piezoelectric ceramics, the vibration conditions are different in diverse regions before or after the resonant points. In order to obtain the ideal working range, especially the frequency interval of the proposed integrated model, the rotor velocities of different amplitudes and frequencies are displayed in Fig. 11a and b.

Figure 11 Velocity response under different driving parameters: (a) speed versus amplitude; (b) speed versus frequency.

It can be seen from Fig. 11a that the linearity between the driving amplitudes and the rotor speed is evident on the whole. And Fig. 11b demonstrates that the speed gradually increases and then drops to a lower value when the frequency is changed from 40 to 44 kHz. Moreover, the frequencies located on the peak velocity decrease as the amplitude increases. This is due to the variation of the natural frequency of the stator caused by the softening nonlinearity. We can conclude that when the frequency comes close to the resonant frequency, the simulation model cannot describe the speed adequately, which may result from the simplification of the stator modal model, as the frequency higher than the resonant peak is usually chosen as the working range, where the fitness becomes better especially from 41 and 44 kHz. The effective fitting results further prove the feasibility of the proposed model.

In all, the above test of speed performances with two driving parameters verifies the validity of the model located in a specific working area, where the amplitude is higher than 0.7 V and the frequency is higher than 41 kHz.

5.2.3 The mechanical characteristics

In order to verify the model in the functions of external load, Fig. 12 illustrates the comparison of mechanical characteristics from the simulations and the experiments. As can be seen, the velocity–torque curve exists in good agreement within the frequencies from 42 to 43 kHz, while the results occur with a little error with the lower frequency. Also, we can observe that the efficiency becomes higher when the frequency is near the resonant value, and the optimal torque we should impose on the motor will be less than 0.5 N m.

Figure 12 Comparison of the mechanical characteristics between simulation and experimentation: (a) 0.9 V, 41.5 kHz; (b) 1 V, 41.5 kHz; (c) 1.3 V, 42 kHz; (d) 1.3 V, 42.5 kHz.

6 Conclusion

This paper presents an electromechanical hybrid model that combines the driving circuit and the TWUSM itself. The main work and contributions can be listed as follows.

1. A driving circuit model combining the circuit components with the load-dependent equivalent circuit model is proposed to simulate the real electric network, which not only supports the TWUSM model, but also helps to gain the input voltages and currents for the calculation of input power.

2. The teeth discretized method is employed to refine the contact status and interaction forces limited to the teeth space, which improves the accuracy of the rotor step response not only in the rising time, but also the steady value.

3. Model agreements are tested on the rotor speed under different parameters (amplitude, frequency, and torque) based on a multi-parameter test system, which demonstrates the feasibility and effectiveness within the ideal frequency working range higher than the resonant frequency.

Moreover, the proposed model achieves the observation of the microscopic characteristics like input currents and vibration response of the transient startup–stopping operation, which are of great significance for precise control of micro-stepping of the TWUSM in the future.

Appendix A

Table A1 The key simulation parameters of USR60.