General 3D Analytical Method for Eddy-Current Coupling with Halbach Magnet Arrays Based on Magnetic Scalar Potential and H-Functions

: Rapid and accurate eddy-current calculation is necessary to analyze eddy-current couplings (ECCs). This paper presents a general 3D analytical method for calculating the magnetic ﬁeld distributions, eddy currents, and torques of ECCs with different Halbach magnet arrays. By using Fourier decomposition, the magnetization components of Halbach magnet arrays are determined. Then, with a group of H-formulations in the conductor region and Laplacian equations with magnetic scalar potential in the others, analytical magnetic ﬁeld distributions are predicted and veriﬁed by 3D ﬁnite element models. Based on Ohm’s law for moving conductors, eddy-current distributions and torques are obtained at different speeds. Finally, the Halbach magnet arrays with different segments are optimized to enhance the fundamental amplitude and reduce the harmonic contents of air-gap ﬂux densities. The proposed method shows its correctness and validation in analyzing and optimizing ECCs with Halbach magnet arrays. The analytical models with 3D MSP and H-formulations are The of four Halbach are The magnetic ﬁeld distributions, eddy currents, torques, fundamental amplitudes, and harmonics of air-gap ﬂux in different PM methodology, X.L.; software, X.H.; validation, Q.H., S.Y. and M.C.; formal analysis, X.L.; investigation, X.H., Q.H., S.Y. and M.C.; resources, X.L.; data curation, X.L.; writing—original draft preparation, X.L. and X.H.; writing—review and editing, P.J., Q.H., S.Y. and M.C.; visualization, X.H.; supervision, P.J.; project administration,


Introduction
As a kind of speed-regulating device based on the principle of asynchronous transmission, cylindrical and radial eddy-current couplings (ECCs) with permanent magnets (PMs) have the advantages of high efficiency, energy saving, low maintenance cost, simple installation, and adaptation to harsh environments [1]. Therefore, they can realize power transmission between motor, fan, and pump loads with satisfactory operational characteristics without a mechanical connection [2,3]. A reliable theoretical model for evaluating magnetic field distributions and torque performances is desirable in a PM ECC design. There are three possible approaches to establishing such a model, which is a solution to the moving conductor eddy-current problem. The first approach is to apply numerical methods, such as the finite element method (FEM), which can consider geometric details, the nonlinearity of magnetic materials, and complicated boundary conditions [4]. The 2D FEM is preferred in the early design stage because of time savings [5]. In [6], a radial 2D FEM solver combined with an axial one enables the prediction of the eddy-current paths in three dimensions with the associated losses. The 3D FEM with a Dirichlet boundary condition iteration is introduced to shorten the solving times [7]. In [8], the authors consider the hysteresis behavior of ferromagnetic material by the 3D FEM. However, the calculation of the FEM is time consuming and cannot provide significant relationships between structural parameters and performances [9].
The second approach is the lumped parameter magnetic circuit method (LPMCM), which requires less programming and computing time. An MEC-based analytical method is outlined to calculate a PM ECC's eddy current and reaction flux [10]. In [11], the magnetic field and eddy-current modulations by iron teeth in a conductor disk are considered by dividing the magnetic flux paths in air gaps into several branches. In [12], the authors calculate an adjustable-speed PM ECC under an asymmetric magnetic field. However, the mentioned LPMCMs can only calculate one PM magnetization direction, namely either the tangential direction or parallel direction, which makes optimization of the magnets impossible [13,14].
As the third approach, the analytical methods derived from simplified boundary conditions and Maxwell's equations are more insightful than the FEM and more accurate than the LPMCM [15][16][17][18][19][20][21][22][23][24][25]. The two most general methods in analytical calculations are magnetic vector potential (MVP) and magnetic scalar potential (MSP). The subdomain technique is used to calculate the 2D magnetic flux density [18]. In [19], the authors calculate the magnetic field distribution of an ECC equipped with disk magnets and high-temperature superconducting pancake coils by assuming an infinite radius of the ferromagnetic yokes. Two-dimensional MSP is also used in an axial field magnetic coupler by assuming the yokes have an infinite radius [20]. In [21,22], 3D MSP is used in the magnet region as well as air gaps for a cylindrical ECC with two-segment Halbach magnets and axial-flux ECCs. In [23], 3D MSP and modified BESSEL functions are used to analyze an air-cored linear and a rotary permanent magnet actuator. In [24], an ECC model directly considers the radial edge effects and the curvature effects on the torque prediction without the need for any correction factor. However, a general analytical analysis with 3D MSP in a cylindrical ECC with Halbach PMs considering the harmonic reduction and optimization has not been presented yet.
This paper presents a 3D analytical analysis of cylindrical ECCs with different Halbach magnet arrays. The coupling analytical models with 3D MSP and H-formulations are established. The magnetization components of four typical Halbach magnet arrays are determined. The magnetic field distributions, eddy currents, torques, fundamental amplitudes, and harmonics of air-gap flux density in different PM arrays are deduced. In addition, the Halbach magnet arrays are optimized to reduce harmonics by selecting the appropriate polar arc coefficient and inclination angle. The 3D simulation and verification are carried out in Cartesian coordinates. Figure 1 shows the structure of an ECC with a Halbach magnet array. It consists of two rotors: an outer rotor and an inner one. The conductor is on the surface of the former, and the Halbach magnet array is on the latter's surface. The second approach is the lumped parameter magnetic circuit method (LPMCM), which requires less programming and computing time. An MEC-based analytical method is outlined to calculate a PM ECC's eddy current and reaction flux [10]. In [11], the magnetic field and eddy-current modulations by iron teeth in a conductor disk are considered by dividing the magnetic flux paths in air gaps into several branches. In [12], the authors calculate an adjustable-speed PM ECC under an asymmetric magnetic field. However, the mentioned LPMCMs can only calculate one PM magnetization direction, namely either the tangential direction or parallel direction, which makes optimization of the magnets impossible [13,14].

Magnetization of Halbach Magnet Arrays
As the third approach, the analytical methods derived from simplified boundary conditions and Maxwell's equations are more insightful than the FEM and more accurate than the LPMCM [15][16][17][18][19][20][21][22][23][24][25]. The two most general methods in analytical calculations are magnetic vector potential (MVP) and magnetic scalar potential (MSP). The subdomain technique is used to calculate the 2D magnetic flux density [18]. In [19], the authors calculate the magnetic field distribution of an ECC equipped with disk magnets and high-temperature superconducting pancake coils by assuming an infinite radius of the ferromagnetic yokes. Two-dimensional MSP is also used in an axial field magnetic coupler by assuming the yokes have an infinite radius [20]. In [21,22], 3D MSP is used in the magnet region as well as air gaps for a cylindrical ECC with two-segment Halbach magnets and axial-flux ECCs. In [23], 3D MSP and modified BESSEL functions are used to analyze an air-cored linear and a rotary permanent magnet actuator. In [24], an ECC model directly considers the radial edge effects and the curvature effects on the torque prediction without the need for any correction factor. However, a general analytical analysis with 3D MSP in a cylindrical ECC with Halbach PMs considering the harmonic reduction and optimization has not been presented yet.
This paper presents a 3D analytical analysis of cylindrical ECCs with different Halbach magnet arrays. The coupling analytical models with 3D MSP and H-formulations are established. The magnetization components of four typical Halbach magnet arrays are determined. The magnetic field distributions, eddy currents, torques, fundamental amplitudes, and harmonics of air-gap flux density in different PM arrays are deduced. In addition, the Halbach magnet arrays are optimized to reduce harmonics by selecting the appropriate polar arc coefficient and inclination angle. The 3D simulation and verification are carried out in Cartesian coordinates. Figure 1 shows the structure of an ECC with a Halbach magnet array. It consists of two rotors: an outer rotor and an inner one. The conductor is on the surface of the former, and the Halbach magnet array is on the latter's surface.     a cylindrical ECC's circumferential and axial cutaways with a threesegment Halbach permanent magnet. The Halbach arrays are distributed along the circumferential direction. The copper is a little longer than the PMs for higher PM utilization rates in the axial direction.

Magnetization of Halbach Magnet Arrays
the inner radius of the inner yoke The divided regions are listed in Table 1.

Region Range of z-Direction
Range of x-Direction I Inner yoke region As the curvature effect can be ignored [23], a circumferential topology can be substituted by an equivalent linear structure, where z = r and x = Rmθ. Then, the relative speed vm is given by: For simplification of the analysis, we adopt the following assumptions: (1) The relative recoil permeability of the permanent magnet region is μr = 1.
(2) The conductivity of the outer yoke is zero, and the permeability is infinite. The material of the inner yoke is iron, and the relative permeability is μiron.
(3) In region IV, the conductivity of copper σ is constant.
By the separate variable method, the residual magnetization vector M  is decomposed as: The divided regions are listed in Table 1.
As the curvature effect can be ignored [23], a circumferential topology can be substituted by an equivalent linear structure, where z = r and x = R m θ. Then, the relative speed v m is given by: For simplification of the analysis, we adopt the following assumptions: (1) The relative recoil permeability of the permanent magnet region is µ r = 1.
(2) The conductivity of the outer yoke is zero, and the permeability is infinite. The material of the inner yoke is iron, and the relative permeability is µ iron . (3) In region IV, the conductivity of copper σ is constant.
The flux density → B and magnetic field intensity → H in regions I, II and III, are given by: By the separate variable method, the residual magnetization vector → M is decomposed as: The amplitude of → M can be given by: Generally, different Halbach arrays are employed in industrial applications, as in Figure 3. The authors have analyzed the magnetic fields of those typical topologies based on 2D magnetic scalar potential [25]. This paper focuses on the 3D analytical analysis of cylindrical ECCs, where the eddy currents are considered. Figure 3 shows the 3D topologies with different Halbach magnet arrays.
The amplitude of M  can be given by: Generally, different Halbach arrays are employed in industrial applications, as in Figure 3. The authors have analyzed the magnetic fields of those typical topologies based on 2D magnetic scalar potential [25]. This paper focuses on the 3D analytical analysis of cylindrical ECCs, where the eddy currents are considered. Figure 3 shows the 3D topologies with different Halbach magnet arrays.
Region Ⅳ Region Ⅲ Region Ⅱ Region І  Figure 4 shows the magnetization components of the arrays, where δx = |sinθ|, δy = |cosθ|. Mx is an odd function for all four typical Halbach magnet arrays, and Mz is an even function along the x-direction. In the z-direction, as the relative velocity is zero, all regions are extended periodically by the imaging method, and Mz is a homogeneous harmonic even function, as shown in Figure 4e.  Figure 4 shows the magnetization components of the arrays, where δ x = |sinθ|, δ y = |cosθ|. M x is an odd function for all four typical Halbach magnet arrays, and M z is an even function along the x-direction. In the z-direction, as the relative velocity is zero, all regions are extended periodically by the imaging method, and M z is a homogeneous harmonic even function, as shown in Figure 4e.
By utilizing the Fourier series method, M x and M z can be expressed by: and: A xk and A zk can be given in Equations (7)-(10) within different topologies.  By utilizing the Fourier series method, Mx and Mz can be expressed by: and: Axk and Azk can be given in Equations (7)-(10) within different topologies. In the two-segment Halbach array: In the two-segment Halbach array: In the three-segment Halbach array: In the four-segment Halbach array: In the ideal Halbach array:

Magnetic Field and Eddy-Current Distribution
Since the permeability of an iron core is assumed to be infinite, the Dirichlet boundary condition where z = R s is satisfied. On both sides of the boundary z = R i , the magnetic field lines are confined mainly in region I. On the boundary z = R i , the Neumann boundary condition can be used approximately. The periodic boundary condition is employed in all directions. Boundary conditions in Figure 2 are given as: Solving the partial differential of MSP φ i (i = I, II, III), the magnet field intensity components in different regions can be found as φ i is expressed by the Laplacian equation as: Then, in regions I and III, the MSP is: In region II, the general solution is given as: In region IV, the time-invariant steady-state is considered. The magnetic field satisfies the static approximation of Maxwell's equations as: Ohm's law for a moving conductor with reference to the stationary frame is expressed as: From Equations (16) and (17), the eddy-current problem in region IV is reduced to a magnetostatic problem, which can be solved using an H-formulation as: By using the method of separation of variables to solve Equations (13) and (19), subject to the distributions of magnetization components and the boundary conditions, general solutions in the four regions are expressed as Equations (20)-(34).
In iron region I: In PM region II: In air-gap region III: In copper region IV: By integrating the magnetic field, the electromagnetic torque exerted on the rotor is as follows:

Comparison with the Finite Element Analysis
In order to validate the proposed model, the analytical results are compared with those obtained by the finite element analysis based on three-dimensional (3D) models, using the commercial software package Ansoft Maxwell. The major parameters of the model used in the validation are given in Section 2. The governing equations and boundary conditions in the finite element simulations are: In the penalty function, the relative permeability of iron µ iron is defined in Equation (37) to make the coefficient matrix symmetrical.
Considering the structural periodicity of the proposed ECC, only a pair of magnetic poles of this volume is needed to be analyzed. Figure 5a shows the field and grid of the ECC. The mapping mesh with a hexahedron shape is adopted, where the numbers of nodes are 151,040, and elements are 142,506. The material's nonlinear properties were considered in the finite element analysis to demonstrate the analytical model's effectiveness and limitations. Since the relative velocity exists in the x-direction, a periodic boundary is used along the direction of motion. Figure 5b shows the eddy-current distribution of the ECC in region IV for a slip speed of 300 rpm. Figure 6 compares the analytically predicted flux densities and the FEM results of the model shown in Figure 3a in the region IV at z = Rc + 1/1000, region III at z = Rm + 3/1000, region II at z = Rm − 1/1000, and region I at z = Rr − 1/1000, respectively. The red line and symbol represent flux densities Br when the slip speed of the ECC is 300 rpm, and the blue line and symbol represent flux densities Bs when the ECC is static.  The material's nonlinear properties were considered in the finite element analysis to demonstrate the analytical model's effectiveness and limitations. Since the relative velocity exists in the x-direction, a periodic boundary is used along the direction of motion. Figure 5b shows the eddy-current distribution of the ECC in region IV for a slip speed of 300 rpm. Figure 6 compares the analytically predicted flux densities and the FEM results of the model shown in Figure 3a in the region IV at z = R c + 1/1000, region III at z = R m + 3/1000, region II at z = R m − 1/1000, and region I at z = R r − 1/1000, respectively. The red line and symbol represent flux densities Br when the slip speed of the ECC is 300 rpm, and the blue line and symbol represent flux densities Bs when the ECC is static. The material's nonlinear properties were considered in the finite element analysis to demonstrate the analytical model's effectiveness and limitations. Since the relative velocity exists in the x-direction, a periodic boundary is used along the direction of motion. Figure 5b shows the eddy-current distribution of the ECC in region IV for a slip speed of 300 rpm. Figure 6 compares the analytically predicted flux densities and the FEM results of the model shown in Figure 3a in the region IV at z = Rc + 1/1000, region III at z = Rm + 3/1000, region II at z = Rm − 1/1000, and region I at z = Rr − 1/1000, respectively. The red line and symbol represent flux densities Br when the slip speed of the ECC is 300 rpm, and the blue line and symbol represent flux densities Bs when the ECC is static.   Figure 7 presents the eddy-current density distributions along the x-direction in the middle of the copper at z = R c + 3/1000 and y = l/2, obtained by proposed analytical Equations (32)-(34) at different speeds. Figure 8 shows the 3D eddy-current density components at z = R c + 3/1000, in which the x-direction component of the induced current is null at y = ±h/2.  Figure 7 presents the eddy-current density distributions along the x-direction in the middle of the copper at z = Rc + 3/1000 and y = l/2, obtained by proposed analytical Equations (32)-(34) at different speeds. Figure 8 shows the 3D eddy-current density components at z = Rc + 3/1000, in which the x-direction component of the induced current is null at = ± / 2 y h .  Figure 9 shows the change in the electromagnetic torque in the analytical method when the slip is changed from 0 to 1 (the rated speed is 1000 rpm). The rated slip and torque of the ECC are 0.02 and 2850 Nm, respectively.   Figure 7 presents the eddy-current density distributions along the x-direction in the middle of the copper at z = Rc + 3/1000 and y = l/2, obtained by proposed analytical Equations (32)-(34) at different speeds. Figure 8 shows the 3D eddy-current density components at z = Rc + 3/1000, in which the x-direction component of the induced current is null at = ± / 2 y h .    Figure 9 shows the change in the electromagnetic torque in the analytical method when the slip is changed from 0 to 1 (the rated speed is 1000 rpm). The rated slip and torque of the ECC are 0.02 and 2850 Nm, respectively.  Figure 7 presents the eddy-current density distributions along the x-direction in the middle of the copper at z = Rc + 3/1000 and y = l/2, obtained by proposed analytical Equations (32)-(34) at different speeds. Figure 8 shows the 3D eddy-current density components at z = Rc + 3/1000, in which the x-direction component of the induced current is null at = ± / 2 y h .  Figure 9 shows the change in the electromagnetic torque in the analytical method when the slip is changed from 0 to 1 (the rated speed is 1000 rpm). The rated slip and torque of the ECC are 0.02 and 2850 Nm, respectively.

Harmonics and Magnetic Field Optimization
Since there are only odd harmonics of the flux density distribution of the Halbach magnet array in the x and z directions, the harmonics can be easily calculated by the analytical prediction method by taking k and n as odd integers. The amplitude of the air-gap flux density is given by Equation (38), which is obtained from Equation (28) by setting x = 0 and cos(kπy/H) = 1.
Generally, to reduce cogging torque and speed fluctuation, the air gap magnetic field distribution in ECCs is required to be sinusoidal. The THD of the magnetic field distributions is given by: The model shown in Figure 3a was analyzed by an analytical method; Figure 10 shows its variation in the THD and the FA of its flux density in region III when α r changed from 0 to 1. The model shown in Figure 3b was analyzed by the analytical method; Figure 11a shows the contour of the THD and the FA of its flux density in region III when α r changed from 0 to 1, and θ changed from 0 to 90 • . The model shown in Figure 3c was analyzed by the analytical method; Figure 11b shows its contour of the THD and the FA of its flux density in region III when α r changed from 0 to 1, and θ changed from 0 to 90 • .
Taking the model shown in Figure 3b as an example, Table 2 lists the preferred and worst results of the THD and FA. When α r = 0.35 and θ = 54 • , the preferred THD is 0.1270, and the preferred FA F 1 is about 1.1051 T. When α r is changed to zero, and θ is changed to 90, the FA decreases rapidly. the x-direction length ratio of the z-direction magnetized PM to τ.
The model shown in Figure 3b was analyzed by the analytical method; Figure 11a shows the contour of the THD and the FA of its flux density in region III when αr changed from 0 to 1, and θ changed from 0 to 90°. The model shown in Figure 3c was analyzed by the analytical method; Figure 11b shows its contour of the THD and the FA of its flux density in region III when αr changed from 0 to 1, and θ changed from 0 to 90°. Figure 11. The predicted contours of THD (black numbers) and FA (red numbers) of (a) three-segment and (b) four-segment Halbach magnet array. αr represents the x-direction length ratio of the z-direction magnetized PM to τ. θ represents the angle of the angular magnetic PM. Figure 11. The predicted contours of THD (black numbers) and FA (red numbers) of (a) threesegment and (b) four-segment Halbach magnet array. α r represents the x-direction length ratio of the z-direction magnetized PM to τ. θ represents the angle of the angular magnetic PM.  Figure 12 compares the flux densities and harmonics of the four Halbach magnet arrays. Figure 12a,b compare the flux densities and harmonics when θ = 45 • and the Halbach arrays are divided into equal parts. Figure 12c,d compare the flux densities and harmonics after optimization. The ideal Halbach magnet array has the smallest harmonic and the largest FA. Figure 13 shows the electromagnetic torques of ECCs with the four Halbach magnetic arrays. The analytical optimization takes less than 10 s to obtain the preferred and worst results, which is much less than the FEM method.
Energies 2021, 14, x FOR PEER REVIEW Taking the model shown in Figure 3b as an example, Table 2 lists the prefer worst results of the THD and FA. When αr = 0.35 and θ = 54°, the preferred THD i and the preferred FA F1 is about 1.1051 T. When αr is changed to zero, and θ is cha 90, the FA decreases rapidly.  Figure 12 compares the flux densities and harmonics of the four Halbach ma rays. Figure 12a,b compare the flux densities and harmonics when θ = 45° and the arrays are divided into equal parts. Figure 12c,d compare the flux densities and ha after optimization. The ideal Halbach magnet array has the smallest harmonic largest FA. Figure 13 shows the electromagnetic torques of ECCs with the four magnetic arrays. The analytical optimization takes less than 10 s to obtain the p and worst results, which is much less than the FEM method.

Conclusions
According to Maxwell's equation, this paper solved the expressions of magnetic field distribution of four typical Halbach arrays. The expressions that have different distribution factors are similar. Three-dimensional analytical magnetic distributions are given using H-formulations in the conductor region and Laplacian equations with magnetic scalar potential in other regions. Eddy currents in the conductor region are obtained from the Ampere law. The results of this method agree with those of the finite element analysis. The comparison verifies the low computational time and good accuracy of the proposed 3D analytical analysis model. The harmonic analysis of the flux densities in the air-gap region is given. Through optimization, the THD of the flux densities decreases, and the amplitude of fundamental waveforms in the air-gap region increases.
The following conclusions can be drawn: (1) Based on the 3D analytical analysis, we directly express the 3D flux densities and intensities in the iron region, copper region, PM region, and air-gap region of the typical Halbach topologies and the eddy current in the copper region and the torque of ECCs. All of the analytical results of the flux densities, eddy currents, and torques are verified by the FEM. The proposed analytical model has high precision and efficiency with less calculation time. (2) The optimization of the parameters of the typical Halbach topologies is very fast, which can reduce the harmonic of ECC, increase the amplitude of a fundamental wave, and improve the torque efficiently.