Precise computation of magnetic fields and static torque is fundamental for the performance analysis and electromagnetic design of harmonic drive gears. The primary computational challenge lies in accurately and efficiently modeling the complex, eccentric motion of the low-speed rotor. Traditional finite element analysis (FEA) is hampered by mesh distortion in the air gap, requiring computationally expensive re-meshing for each rotor position. This paper presents a novel hybrid methodology that seamlessly integrates a two-dimensional finite element model with an analytical double air-gap macro-element model. This approach effectively decouples the stator and rotor meshes, enabling unrestricted eccentric rotor motion while inherently accounting for the non-linear magnetic properties of the ferromagnetic materials.
The proposed hybrid model strategically partitions the eccentric air-gap region into three distinct layers. The two layers adjacent to the stator and rotor, characterized by uniform width, are treated as macro-elements. The middle layer, which contains the non-uniform permeance wave, along with the stator and rotor cores, are discretized using finite elements. This division is conceptually illustrated below, showing the physical structure of an eccentric magnetic harmonic drive gear with its key components and the proposed modeling layers.
The fundamental principle of the harmonic drive gear relies on the modulation of the magnetic field. The high-speed rotor, or “wave generator,” creates a rotating permeance wave in the air gap. This wave interacts with the magnetic field of the eccentric, low-speed permanent magnet rotor, which undergoes a planetary motion—rotating slowly about its own axis while its center revolves at high speed around the stator center. This interaction enables high torque transmission at a substantial gear ratio.
Mathematical Formulation of the Hybrid Model
The magnetic field within the FE regions (stator, rotor, and middle air-gap layer) is governed by Maxwell’s equations. For a 2D problem, using the magnetic vector potential A (with only a z-component, Az), the governing equation for regions containing no currents is derived from $$ \nabla \times ( u \nabla \times \mathbf{A} ) = 0 $$, where \( u \) is the magnetic reluctivity. In regions with permanent magnets, the equation becomes $$ \nabla \times ( u \nabla \times \mathbf{A} ) = \nabla \times ( u \mathbf{B}_r ) $$, where \( \mathbf{B}_r \) is the remanent flux density of the magnet. Applying the standard Galerkin weighted residual method and discretizing the domain with first-order triangular elements leads to the global matrix system:
$$ \mathbf{K}\mathbf{A} = \mathbf{f} $$
where \( \mathbf{K} \) is the stiffness matrix, \( \mathbf{A} \) is the vector of nodal magnetic vector potentials, and \( \mathbf{f} \) is the source vector from the permanent magnets. The elemental stiffness matrix coefficient \( k_{ij}^e \) and source vector component \( f_i^e \) for an element \( \Omega_e \) are given by:
$$ k_{ij}^e = \iint_{\Omega_e} u \left( \frac{\partial N_i}{\partial x} \frac{\partial N_j}{\partial x} + \frac{\partial N_i}{\partial y} \frac{\partial N_j}{\partial y} \right) \,dx\,dy $$
$$ f_i^e = \iint_{\Omega_e} u \left( B_{ry} \frac{\partial N_i}{\partial x} – B_{rx} \frac{\partial N_i}{\partial y} \right) \,dx\,dy $$
Here, \( N_i \) and \( N_j \) are the shape functions for nodes i and j, respectively.
The Double Air-Gap Macro-Element Model
The core innovation is the treatment of the two uniform air-gap layers as macro-elements. In these annular regions with constant inner radius \( R_q \) and outer radius \( R_p \), and no current sources, the magnetic vector potential satisfies Laplace’s equation, \( \nabla^2 A = 0 \). In polar coordinates, the general analytical solution can be expressed as a series expansion. To couple this analytical solution with the finite element mesh on its boundaries, a macro-element formulation is used. The potential at any point \( (r, \theta) \) within the macro-element is interpolated from the nodal potentials on its inner and outer boundaries:
$$ A(r, \theta) = \sum_{i=1}^{M_p} \alpha_i(r, \theta) A_{pi} + \sum_{j=1}^{M_q} \beta_j(r, \theta) A_{qj} $$
where \( A_{pi} \) and \( A_{qj} \) are the vector potentials at the \( M_p \) nodes on the outer boundary (connected to the stator or middle layer) and the \( M_q \) nodes on the inner boundary (connected to the rotor or middle layer), respectively. The shape functions \( \alpha_i \) and \( \beta_j \) are derived from the analytical solution and Fourier expansion to ensure continuity of the potential field. For the outer boundary shape function \( \alpha_i \), the expression is:
$$ \alpha_i(r, \theta) = g_{a,i} \frac{\ln(r/R_q)}{\ln(R_p/R_q)} + \sum_{n=1}^{\infty} \left[ g_{an,i} \cos(n\theta) + g_{bn,i} \sin(n\theta) \right] \frac{E_n(r, R_q)}{E_n(R_p, R_q)} $$
where \( E_n(u,v) = (u/v)^n – (v/u)^n \). The coefficients \( g_{a,i}, g_{an,i}, g_{bn,i} \) are determined by the angular positions of the boundary nodes to satisfy the nodal interpolation condition \( \alpha_i(r, \theta_k) = \delta_{ik} \). Similar expressions exist for the inner boundary shape functions \( \beta_j \).
The magnetic energy stored in a macro-element is:
$$ W_{age} = \frac{1}{2} \iint_{\Omega_{age}} u_0 |\nabla A|^2 \, d\Omega = \frac{L_{ef} u_0}{2} \int_{R_q}^{R_p} \int_{0}^{2\pi} \left[ \left( \frac{\partial A}{\partial r} \right)^2 + \frac{1}{r^2} \left( \frac{\partial A}{\partial \theta} \right)^2 \right] r \, d\theta \, dr $$
Minimizing this energy functional with respect to the nodal potentials \( A_{pi} \) and \( A_{qj} \) yields a set of linear equations of the form \( \mathbf{K}_{age} \mathbf{A}_{age} = 0 \), where \( \mathbf{K}_{age} \) is the macro-element stiffness matrix. This matrix, though full (non-sparse), is small and is directly assembled into the global finite element stiffness matrix \( \mathbf{K} \). This coupling creates the complete system equation for the entire harmonic drive gear domain, seamlessly blending the analytical and numerical domains.
Electromagnetic Torque Calculation
The electromagnetic torque transmitted between the rotors is calculated using the Maxwell Stress Tensor method. The advantage of using the macro-element formulation is that the torque can be derived as an explicit, path-independent function of the boundary nodal potentials, offering higher accuracy than integrating the discretized flux density from the FE mesh. The general torque formula for a cylindrical surface of radius \( r \) within a macro-element is:
$$ T_{em} = \frac{L_{ef} r^2}{u_0} \int_{0}^{2\pi} B_r(r, \theta) B_{\theta}(r, \theta) \, d\theta $$
where \( B_r = \frac{1}{r} \frac{\partial A}{\partial \theta} \) and \( B_{\theta} = -\frac{\partial A}{\partial r} \) are the radial and tangential flux density components. Substituting the analytical expressions for \( A \) from the macro-element model leads to a closed-form expression for the torque. For the outer macro-element (between stator and permeance wave), the torque on the stator \( T_1 \) is:
$$ T_1 = \frac{\pi L_{ef} u_0}{2} \sum_{n=1}^{\infty} \frac{n}{E_n(R_p, R_q)} \left[ \left( \sum_{i} g_{bn,i} A_{pi} \right) \left( \sum_{i} g_{an,i} A_{pi} \right) + \left( \sum_{j} h_{dn,j} A_{qj} \right) \left( \sum_{j} h_{cn,j} A_{qj} \right) \right] $$
A similar expression yields the torque \( T_2 \) on the rotor from the inner macro-element. The torque ratio \( T_2 / T_1 \) should theoretically equal the inverse ratio of the pole pairs, validating the force balance in the harmonic drive gear.
Handling Eccentric Motion and Improving Computational Speed
The decoupling of meshes via macro-elements naturally allows for free rotor motion. The planetary motion of the low-speed rotor, where it rotates by an angle \( \gamma \) about its own axis while its center revolves by \( p_r \gamma \) about the stator axis (with \( p_r \) being the rotor pole pairs), is implemented in two steps within the model.
However, a naive implementation would require re-calculating the dense macro-element matrix \( \mathbf{K}_{age} \) for every rotor position because the boundary node coordinates change. To solve this, an equivalent motion mode combined with the sliding surface technique is proposed. The reference frame is shifted to the center of the low-speed rotor. In this frame, the stator and the rotor appear to rotate in the same direction at different speeds. More critically, the interfaces between the macro-elements and the moving parts (stator and rotor) can now be treated as “sliding surfaces.”
In this technique, the nodes on the macro-element side of the interface are fixed in space. The nodes on the moving part’s side (e.g., the rotor surface) slide along this fixed interface. The magnetic vector potential at a moving node is assigned the value from the fixed node it currently coincides with. This relationship is managed by a cyclic permutation of node indices. Consequently, the geometry of the macro-element boundaries remains unchanged during rotation, and the matrix \( \mathbf{K}_{age} \) needs to be computed only once at the beginning of the simulation, leading to dramatic computational savings.
Numerical Validation and Prototype Testing
The proposed method was applied to analyze a dual-stage eccentric magnetic harmonic drive gear prototype. The key parameters for the second stage, which is the focus of the magnetic analysis, are summarized in the table below.
| Parameter | Value |
|---|---|
| Stator Outer Diameter | 103 mm |
| Rotor Inner Diameter | 55 mm |
| Stator/Rotor Pole Pairs (ps/pr) | 16 / 15 |
| Permanent Magnet Thickness | 3.5 mm |
| Remanent Flux Density (Br) | 1.25 T |
| Nominal Air-gap Length | 5 mm |
| Eccentricity | 4 mm |
| Axial Length (Lef) | 20 mm |
The magnetic field distribution computed by the hybrid model clearly shows the modulation effect. The radial flux density waveforms in the two air-gap macro-elements, when decomposed, show dominant space harmonics corresponding to the 16-pole stator field and the 15-pole rotor field, confirming the fundamental operating principle of the harmonic drive gear.
The static torque-angle characteristic was calculated. The peak holding torque for the second stage was found to be 19.16 N·m. The calculated torque ratio \( T_2/T_1 \) was approximately 0.94, closely matching the theoretical pole pair ratio of 15/16 ≈ 0.9375.
A prototype was built and tested. The experimental setup measured the steady-state output torque versus speed. At an input speed of 1820 rpm and an output speed of 100 rpm (corresponding to the gear ratio), the maximum sustained output torque was measured to be 17.2 N·m. This is reasonably aligned with the calculated static peak torque of 19.16 N·m, considering that the steady-state operating torque is necessarily lower than the theoretical static pull-out torque due to factors like minor misalignments, temperature, and dynamic effects.
A critical comparison was made between the proposed macro-element torque calculation and the traditional finite element method using a stress tensor integration path through the mesh. The results are compared in the table below, highlighting the stability and accuracy advantage of the macro-element approach, as it avoids errors from discretized field quantities.
| Method | Calculated Peak Torque | Waveform Stability | Key Characteristic |
|---|---|---|---|
| Proposed Hybrid Model (Macro-Element) | 19.16 N·m | Excellent | Analytical torque expression; High accuracy from nodal potentials. |
| Standard FEM (Mesh Integration) | 20.64 N·m | Good | Sensitive to mesh quality and integration path; Lower accuracy from discretized B-field. |
| Prototype Measurement (Steady-State) | 17.2 N·m | N/A | Practical operating value, lower than static peak. |
The computational efficiency of the equivalent motion mode was also validated. For a simulation involving 46 rotor positions, recalculating the macro-element matrices for each step took 7870 seconds. Using the proposed method with fixed macro-element boundaries and the sliding surface technique reduced the computation time to 2586 seconds, a time saving of over 65%.
Conclusion
This paper successfully developed and validated a hybrid finite element and double air-gap macro-element model for the comprehensive magnetic analysis of eccentric magnetic harmonic drive gears. The model effectively solves the critical challenge of modeling large eccentric motion by decoupling the stator and rotor meshes, allowing for unrestricted rotor movement without remeshing. The incorporation of analytical macro-elements not only facilitates this free motion but also provides a more accurate and stable method for calculating electromagnetic torque compared to standard finite element post-processing. The introduction of an equivalent motion mode coupled with a sliding surface technique dramatically improves computational efficiency by eliminating the need to recalculate macro-element coefficients during rotation. The close agreement between the calculated results and experimental measurements on a prototype harmonic drive gear confirms the validity, accuracy, and practical utility of the proposed methodology for the design and analysis of this promising class of high-torque, non-contact magnetic gears.