Abstract
This paper provides an analytical method for calculating the electromagnetic field of the tubular linear permanent magnet eddy current brake. The method of separating variables based on Fourier series is used to solve the boundary value problem, and the analytical solution of electromagnetic field is obtained. The braking force is obtained by calculating the Lorentz force. Then the accuracy of the analytical method is verified by comparing the analytical solutions with the solutions obtained by finite element method (FEM). In order to improve braking performance in a limited space, the parameter analysis is carried out to investigate the influence of various structural parameters on the braking characteristics. According to the results of parameter analysis, the range of design parameters is determined. Finally, a multi-objective optimization is applied to find the optimal design parameters of the eddy current brake. The braking force increases by 27% and the volume of permanent magnet decreases by 7.7% after optimization.
Introduction
Eddy current brakes are widely used in many industrial applications, such as aerospace applications, energy recovery system, high speed train braking systems, automobile suspension, and auxiliary braking devices, due to their advantages of simple structure, no mechanical contact, high reliability and easy maintenance [1–5]. Eddy current brakes can be divided into three types: electric excitation [6–8], permanent magnet [9–11], and hybrid excitation [12,13] eddy current brakes according to different excitation sources. The permanent magnet eddy current brake does not require additional power, is small in size and more reliable. Compared with planar linear eddy current brakes, tubular linear eddy current brakes have lower magnetic leakage, higher utilization of permanent magnets and no transverse edge effect [14].
The research methods of eddy current brake include analytical method and finite element method (FEM). The FEM can be applied to various complex structures, but it is time consuming. For some eddy current brakes with relatively regular structures, analytical methods can be used to analyze them. Bae et al. applied the Biot-Savart law and the Lorentz force formula to solve the magnetic flux density, then investigated the characteristics of eddy current damping when a permanent magnet was placed in a conductive tube. However, the predictability of this model is not accurate at high excitation frequency because the reaction of eddy current is not considered [15]. Ebrahimi et al. also studied the electromagnetic force of the tubular eddy current brake through Biot-Savart law and the Lorentz force formula, and considered the skin effect on the electromagnetic force. However, the eddy current reaction field was ignored and the effect of soft magnetic materials was not considered in the derivation [16]. In order to accurately account the eddy current reaction field, the electromagnetic field was calculated by solving Maxwell’s equations in some papers. Wang et al. deduced the magnetic field distribution of the tubular linear permanent magnet machines according to the magnetic vector potential and the cylindrical coordinate formula. Then he compared it with the finite element analysis to verify the results extensively [17]. However, only static results are given. Bianchi expressed the electromagnetic field problem by means of the magnetic scalar potential, then solved the Poisson’s equation and Laplace’s equation using separation variable method. In order to reduce the difficulty of solving the Poisson equation, they introduced a convenient approximation in the permanent magnet magnetization [18]. Magnetic scalar potentials can only be used in current-free regions, therefore cannot account the eddy current reaction field. Park et al. used the same approximation to express magnetization, then calculated the static air gap magnetic flux density of a tubular linear synchronous machine [19]. Amara et al. predicted the eddy current loss in the moving armature of a tubular permanent magnet machine using equivalent current sheets [20]. Although simplified models with equivalent current sheets may accurately predict the performance of some specific machines, Atallah et al. indicated that using equivalent current sheets can imply a loss of accuracy [21].
Most researches focus on the modeling process of eddy current brakes, then analyze the braking force characteristics and the influence of various structural parameters on the braking force characteristics [22–24]. According to different working conditions, there are different requirements for the performance of ECB. In order to play a better role in specific application scenarios, it is necessary to coordinate and match ECB comprehensively. However, the structural parameters of the ECB have not been further optimized to achieve the desired braking force characteristics in those researches. In addition, there are more than one indicator for evaluating the braking force characteristics. Multi-objective optimization techniques can be used to deal with multi-objective optimization problems.
In this paper, the reaction of the eddy current field is fully considered, and the analytical model of a tubular linear permanent magnet eddy current brake is established based on the magnetic vector potential and the cylindrical coordinate system. The residual magnetism model is used to represent the permanent magnet instead of the equivalent current sheets model because of the equivalent current sheets model is not accurate for all mechanical topologies. Then the boundary value problem of partial differential equations is solved by the method of separating variables and Fourier series. Analytical expressions of electromagnetic fields under quasi-static conditions are obtained. The accuracy of the analytical model is verified by comparison with the results calculated by the FEM. Through harmonic analysis, it is found that the braking force with sufficient accuracy can be obtained by calculating the 1, 5 and 9 orders. When the harmonic order increases, it is necessary to increase the number of digits of numerical calculation to obtain accurate results. Considering that the optimization process requires a lot of iterative calculations and the performance of ECB is difficult to evaluate with a single index, multi-objective optimization of the ECB is carried out by establishing a surrogate model and combining with GA.
Analytical model
Analytical modeling

Structure of the tubular linear permanent magnet eddy current brake with Halbach array.
Figure 1 shows a schematic of a tubular linear permanent magnet eddy current brake with Halbach array. The Halbach array can increase the magnetic flux density on one side and decrease it on the other. Therefore, the magnetic permeability of the support tube with Halbach array has a small impact on the braking force. The permanent magnets are installed on the support tube. The conductive tube is covered by back iron, which can significantly enhance the air gap magnetic flux density. When the support tube moves, eddy current is generated in the conductive tube. Then the eddy current interacts with the magnetic field to generate a braking force.

Simplified analytical model.
Figure 2 shows the simplified analytical model. The analysis will be performed in a cylindrical coordinate system (r, φ, z) fixed to the permanent magnet. r
0 and r
1 are the inner and outer radii of the support tube; r
2 is the outer radii of the permanent magnets; r
3 and r
4 are the inner and outer radii of the nonmagnetic conductive tube. b and c are the axial length of axial magnetization and radial magnetization of the permanent magnets. τ is the pole pitch. h
m
, δ and d are the radial thickness of permanent magnets, air gap and nonmagnetic conductive tube, respectively. In the particular machine under consideration, the permanent magnets are mounted on a high strength steel tube. In order to obtain an analytical solution for the electromagnetic field, the following assumptions are made: The z axial end effect is ignored so that the electromagnetic field is periodic distributed along the z axis. The relative permeability of ferromagnetic materials is constant and the radial thickness of back iron is infinite. For it has been observed that this approximation does not modify the performance prediction of the eddy current brake when the radial thickness of the back iron is large [25]. The current density distribution is supposed directed only along the z axis. Therefore, the electromagnetic field distribution can be described as a two-dimensional problem.
The analytical model is divided into five sub-regions. Assume that only region IV has a non-zero conductivity. The assumption that the conductivity of region V is zero will be explained in Section 2.2. So in regions I, II, III and V,
The constitutive relation of permanent magnet can be described by residual magnetism

Waveform of residual magnetism.
The rectangular waveform can be expressed by Fourier series
Therefore, the governing equations of region II is
Since all the electromagnetic quantities are periodic with the pole pitch τ, the magnetic vector potential can be written as [25]
Solve the above ordinary differential equation. Then the magnetic vector potential of each region described by the modified Bessel equations are obtained.
The boundary conditions ((16)) and ((17)) indicate that the radial component of magnetic flux density vanishes at r
0 and infinity. Then the equations ((18))–((21)) ensure the continuity of the normal component of magnetic flux density at the interface of different regions. Equation ((22))–((25)) ensure the continuity of the tangential component of the magnetic field strength. Substituting ((13))–((15)) into the boundary conditions results in a system that can be written as ((26)).
According to the solution of the magnetic vector potential, the magnetic flux density and eddy current distribution in the nonmagnetic conductive tube can be obtained
Finally, the Lorentz force can be obtained
The analytical model is verified by comparing the analytical results with numerical solution obtained by the FEM. The calculations are performed using the parameters given in Table 1. The axial length of an actual permanent magnet array is finite. In order to eliminate end effects, the braking force calculated by the finite element model comes from one period of the permanent magnet array. The FEM analysis is implemented using the COMSOL Multiphysics software. Choose the AC/DC module to solve the Maxwell equations formulated using vector magnetic potential to establish a two-dimensional axisymmetric finite element model. The local finite element model is shown in Fig. 4. The setting of constitutive relations is consistent with the analytical model. Give the conductor tube and the back iron a very small acceleration to calculate the braking force under quasi-static. Figure 5 shows the magnetic field distribution at 0 m/s and 20 m/s. It can be seen that in the static state, the magnetic induction lines pass through the conductor cylinder almost perpendicularly from the permanent magnet magnetized radially. When the eddy current brake is working, due to the reaction of the eddy current field, the magnetic field distribution is distorted.

Local finite element model.
Parameters of the eddy current brake used in the calculations

Magnetic flux density distribution at the speed of 0 m/s (a) and 20 m/s (b).
Figure 6 compares the analytical results and the finite element results of B r in the middle of the nonmagnetic conductive tube (r = r 3 + d∕2), whereas Fig. 7 compares the B z . The analytical results and the finite element results very good agreement. It can be observed that, consistent with the results shown in Fig. 5, the magnetic flux density of the static magnetic field forms a platform at the permanent magnet magnetized radially. The reaction of the eddy current field will cause the platform to disappear into an extreme point, showing the effect of field suppression.

Comparison of magnetic flux density B r as functions of z at r = r 3 + d∕2.

Comparison of magnetic flux density B z as functions of z at r = r 3 + d∕2.
The curves of current density as functions of v at different points in the middle of nonmagnetic conductive tube are shown in Fig. 8. The analytical and finite element curves also show excellent consistency. Due to the distortion of magnetic induction intensity distribution, the eddy current distribution has also changed, and the peak position of the eddy current has shifted.

Comparison of current density as functions of v at different points (r = r 3 + d∕2).

Comparison of braking force as functions of v.
Figure 9 is braking force as the function of speed. Since the analytical expression of the magnetic field distribution contains a lengthy combination of Bessel functions, the braking force is obtained by numerical integration. As the harmonic order increases, in order to obtain accurate numerical solutions of the boundary condition equations, the number of digits in the numerical calculation process must be increased. Otherwise, the result may be wrong. This is discussed in detail in the next section. At low speeds, the eddy current damping coefficient is roughly linear. Due to the reaction of eddy current, the magnetic field distribution is distorted, the damping coefficient changes from linear to non-linear, and the braking force has a peak value. This analytical method can accurately describe this physical phenomenon. Comparing the calculation results of the analytical method with the results of the FEM shows that the two methods are in good agreement. Therefore, the accuracy of the analytical method is verified. The quantitative comparison of the error between the braking force obtained by the analytical method and the braking force obtained by the FEM in Fig. 8 is discussed in Section 2.3.
When the conductivity of region V is 11.2 MS/m, the braking force as the function of speed is shown in Fig. 10. It can be seen that eddy currents mainly appear in region IV, so the assumption that the conductivity of region V is zero is reasonable.

Comparison of braking force of different tubes.
The analytical model based on Fourier series exhibits a similar problem. The accuracy of the analytical model is related to the amount of harmonic orders used. Since the computational time of the model mainly depends on the amount of harmonic orders, considering the computation cost, the amount of harmonic orders to be calculated cannot be infinite. When the analytical model is actually used, the amount of harmonic orders calculated should be determined according to the numerical accuracy and computational time. Moreover, the amount of harmonic orders required to calculate the magnetic flux density, eddy current density, and braking force may be different.
A method to reduce computational time is to perform a preliminary harmonic analysis of the machine under study to exclude a series of harmonic orders. This reduces the amount of harmonic orders that need to be calculated, thus reducing computational time. As shown in Fig. 3, the spatial distribution of the remanence is symmetrical about r axis.
Obviously, the harmonic orders can exclude even terms according to (33) and (34). As can been observed from Fig. 11, the main harmonic orders in the radial magnetic flux density are first, fifth and ninth orders. The harmonic amplitudes of the third, seventh, eleventh etc. are very small and can be ignored.
The relative difference between analytical (ANA) and FEM results are quantified using the normalized root mean square deviation (NRMSD) [26].

Harmonic distribution of radial magnetic flux density in a conductive tube.

NRMSD of B r and F z at different amount of harmonic order.
Figure 12 shows the NRMSD of braking force is only 0.083% with the highest order reaching nine, whereas the NRMSD of radial magnetic flux density is 0.15% with the highest order reaching twenty-five. It can also been seen that the contribution of third, seventh, eleventh etc. harmonics to the result can be neglected.
For an increased harmonic number, solving the equations for the boundary conditions results in ill-conditioned system of equations; hence, the solution becomes inaccurate [27]. In addition, there are inevitably rounding errors in the numerical calculation process. The expressions of the unknown coefficient
For (36), sufficient number of digits are needed to obtain correct results. The larger the harmonic order, the more digits are required, otherwise the result of (36) is zero. For the solution of system (26), an increase in speed can also lead to inaccurate solutions. As the harmonic order and speed increase, in order to obtain an accurate solution, the number of digits in the numerical calculation process must be increased. However the increase in the number of digits will undoubtedly increase the time consumed for computation.
In order to verify the accuracy of analytical method and FEM, a small prototype is manufactured. The dominating specifications of the prototype are summarized in Table 2. The experimental device is depicted in Fig. 13. The air-hammer hits the slider to give speed to the moving part of the ECB. Figure 14 shows the test results of the small prototype. It can be seen that the test results are in good agreement with the FEM results. The accuracy of the FEM is validated. The analytical method does not consider the end effect, so the calculated result is larger than that of the FEM.

Experimental device of the testing system.
Specifications of small prototype

Comparison between FEM, analytical, and experimental results.
The reasonable design of the eddy current brake can save materials and improve performance. Optimizing structural parameters requires a lot of calculations. Considering the computing cost, an approximate model is established instead of the original model.
Parametric analysis
The performance indicators for evaluating eddy current brakes include damping coefficient at low speed, peak braking force, critical speed (speed corresponding to peak braking force), damping density, and braking force density, where damping density is defined as the ratio of the damping coefficient at low speed over the permanent magnets volume and braking force density is defined as the ratio of the peak braking force over the permanent magnets volume.
The performance of the eddy current brake is related to many structural parameters, such as the axial length of axially magnetized and radially magnetized permanent magnets and the radial thickness of permanent magnets and nonmagnetic conductive tube. It is necessary to analyze the relationship between structural parameters and braking characteristics, before optimizing structural parameters. Then the design parameters and corresponding value ranges are determined according to the results of the parameter analysis.
As shown in Fig. 15, the increase of d leads to the increase of damping coefficient c. Nevertheless, the damping coefficient will not increase indefinitely, for an increase in d means an increase in magnetic resistance, and the existence of skin effect. On the other hand, the critical speed and the peak braking force are reduced with the increased eddy current reaction.
In Fig. 16, the braking force curves are shown for different permanent magnets thickness h m . It can be seen that the increase in h m results in an increase in the peak braking force and an increase in the damping coefficient. However, there are optimum values of h m for the braking force density and the damping density, because the magnitude of the increase in braking force density and damping density will gradually decrease.

Braking force as functions of v at different conductor plate thicknesses.

Braking force as functions of v at different permanent magnets thicknesses.
As can be seen from Fig. 17, b clearly has the optimal value. The increase in b increases both the magnetomotive force and the reluctance. Properly increasing b can effectively, but not infinitely increase the air gap magnetic flux density.
It can been observed in Fig. 18 that the relationship between c and braking force characteristics is similar to that between h m and braking force characteristics. Both are mainly changes in the braking force characteristics caused by changing the volume of the radially magnetized permanent magnet and the size of the eddy current region. Therefore, there is also an optimal value of c for the braking force density and damping density.

Braking force as functions of v at different axial lengths of axial magnetized permanent magnets.

Braking force as functions of v at different axial lengths of radial magnetized permanent magnets.
In the previous section, the influence of design parameters on braking force characteristics was analyzed. Different work scenarios require suitable braking force characteristics. In the particular eddy current brake considered, assuming the maximum working speed does not exceed 13 m/s. So 13 m/s is selected as the working speed, which is the end of the profile of the moving speed, related to braking force, during actual work. The entire profile of braking force needs to be considered when optimizing the design. Therefore, this paper chooses the damping coefficient at low speed and the braking force at working speed to characterize the overall profile. The damping coefficient represents the speed at which the braking force rises, and the braking force corresponding to the working speed represents the braking force at the end of the profile. The working speed of the original structure is lower than the critical speed, that is, the braking force at the working speed of the original structure is the maximum braking force during actual work. Therefore, the entire profile of the moving speed, related to braking force, can be estimated based on the damping coefficient and the braking force under the working speed. The braking force at working speed and the damping coefficient can be increased by reducing some evaluation indicators, such as the critical speed and the peak braking force.
The sensitivities of the braking force at the working speed and the damping coefficient with respect to b and c are shown in Figs 19 and 20. It can be seen that compared to b, the braking force and damping coefficient are more sensitive to c. In the specific combination of b and c, the braking force is inversely related to b, which must be avoided.

Sensitivity of F z with respect to b and c.

Sensitivity of c d with respect to b and c.
In this study, the optimization is performed using back-propagation (BP) neural network surrogate model, because the FEM is time consuming, the unknown coefficients in the analytical expressions are too verbose and complex, and (32) it is difficult to give an analytical expression. The genetic algorithm optimization solving strategy based on the back-propagation neural network surrogate model is shown in Fig. 21.

Genetic algorithm optimization solving strategy based on the back-propagation neural network surrogate model.
Table 3 shows the design variables and corresponding ranges. The inner diameter of the permanent magnet and the thickness of the back iron remain unchanged according to the actual work needs. The optimization goal is to increase the braking force at the working speed and damping coefficient, while reducing the volume of the permanent magnet. In order to ensure that the back iron is not magnetically saturated, the air gap magnetic flux density B
g
is limited to less than 1.46 T. So the multi-objective optimization problem can be described as
Design variables
Sample points are obtained using the optimal Latin square method. Then set up a parametric finite element model to obtain training samples. The accuracy of the numerical results is guaranteed by comparison with the analytical results. Given that there are many local optimal value of the objective function, GA is used for global optimization. GA has better global searching ability compared with derivative-based optimization methods, which can effectively search for the optimal solution in the design space without becoming trapped in the local optimal solution.
20 sets of samples are selected using the optimal Latin hypercube experiment design method and the finite element model and then they are used to check the prediction accuracy of the surrogate models. The prediction accuracy measured by the determination coefficient R
2 can be obtained by following equation:
The Pareto optimal front obtained by multi-objective optimization calculation is illustrated in Fig. 22. According to the preference of actual work requirements, the optimization parameters are selected from the Pareto optimal solution set. Comparison of results before and after optimization are shown in Table 4. Figure 23 compares the braking force characteristics before and after the eddy current brake optimization. The entire profile of the moving speed after optimization (related to the braking force) is located above the profile before optimization. The braking force has been improved throughout the working speed range. It can be found that the optimized eddy current brake increases the damping coefficient and the braking force at general working speed by reducing the peak braking force and the critical speed. The optimization results show that the volume of the permanent magnet is reduced by 7.7%, whereas the braking force at general working speed is increased by 27%. The damping density and braking force density increased by 90.8% and 37.7% respectively.

Pareto optimal front of the optimization result.

Braking force characteristics before and after optimization.
Optimization results
In such devices, the radial and axial components of the Lorentz force are very high. Thus these forces should be considered [28,29]. After the optimization of the ECB in this paper, the overall braking force characteristics have been improved. At the same time, it is necessary to ensure that the mechanical strength of the optimized ECB meets the requirements of use. Using the COMSOL Multiphysics software, the AC/DC module and the structural mechanics module are coupled to calculate the stress distribution of the ECB. The material properties of conductive tube and back iron used in ECB are presented in Table 5. The axial force density and radial force density distribution of the conductive tube are shown in the Fig. 24. The conductive tube and the back iron are interference fit, and the interference is 0.4 mm. The coefficient of friction between the two is 0.17. The stress distribution of the inner tube is shown in Fig. 25. The maximum stress is less than the yield stress, so the optimized ECB is safe.
Material properties of conductive tube and back iron

The axial force and radial force density distribution of the conductive tube.

The stress distribution at the end of the conductive tube.
Based on the Fourier series and the method of separating variables, an analytical model in a cylindrical coordinate system is established in this paper. The results of the analytical model are in good agreement with the results of the finite element method. Through harmonic analysis, it can be found that the error of braking force is only 0.083% by calculating only the first, fifth and ninth order harmonics. A multi-objective optimization design of eddy current brake is carried out by combining BP neural network surrogate model and GA. After optimization, the braking force of the eddy current brake increased by 27%, the damping coefficient increased by 76%, and the volume of the permanent magnet decreased by 7.7%.
Footnotes
Acknowledgements
This work was supported by the National Natural Science Foundation of China under Project 11572158 and 51705253, and the Postgraduate Research & Practice Innovation Program of Jiangsu Province under Grant KYCX19_0334.
Appendix
The expressions of constants
