Rectangular Closed Double Magnetic Circuit O ﬀ ering Ultra-Long Stroke for Ultra-Low-Frequency Vibration Exciter

: High-performance magnetic circuit o ﬀ ering uniform magnetic ﬂux density (MFD) along ultra-long stroke is the key to develop a vibration exciter for ultra-low-frequency (ULF) vibration calibration. In this paper, a rectangular closed double magnetic circuit (RCDMC) o ﬀ ering ultra-long stroke up to 1.2 m is modeled and optimized. In order to overcome the modeling di ﬃ culty arising from the long stroke, a high-accuracy theoretical model is established taking advantage of the structural symmetry of the RCDMC through lumped parameter magnetic equivalent circuit method. Matrix equations are derived based on Kirchho ﬀ ’s law and solved by iteration calculation to deal with the strong nonlinear characteristics of the yoke material. The deviations between the model and ﬁnite element method (FEM) analysis results are less than 1% for non-saturated yokes and ~10% for saturated yokes. Theoretically, an MFD up to 122 mT and an acceleration waveform harmonic distortion (AWHD) as low as 0.45% are achieved through model-based optimization. Experiments are carried out using an RCDMC prototype assembled in a horizontal vibration exciter. The experimental results show that an MFD of 102 mT and an AWHD of 0.27% along 1.2 m stroke are achieved, making the proposed RCDMC a solution for ULF vibration exciter.


Introduction
Calibration of ultra-low-frequency (ULF) (below 0.1 Hz) vibration is urgently needed in various fields including aerospace engineering, precision machining and manufacturing, structural dynamic analysis, earthquake forecast, sea-wave monitoring and warning, mineral exploration, etc. [1][2][3]. However, it is quite difficult due to a variety of reasons, especially the lack of the high-precision vibration exciter that can generate ULF standard vibration with very low acceleration waveform harmonic distortion (AWHD) [4]. For the development of the high-precision low-frequency electromagnetic vibration exciter, high-performance magnetic circuit offering uniform magnetic flux density (MFD) along the long stroke (above 1 m) is essential [5], which is also the technological know-how of a couple of vibration exciter manufacturers in the world. Although the performance advantages of closed double magnetic circuits (CDMCs) have been demonstrated [6][7][8], on the one hand, it places extremely high demands on permanent magnet (PM) material and the precision of machining and assembly, and on the other hand, theoretical modeling for the structure design is demanding.
Compared to short stroke magnetic circuits, great difficulty in theoretical modeling is presented in developing an ultra-long stroke CDMC for an ULF metrological vibration exciter. The reluctance caused by the yokes, which is ignorable in short-stroke magnetic circuits [9][10][11][12], becomes a key factor limiting the uniformity of MFD distribution. The MFD curve which is approximately a straight line in short-stroke magnetic circuits is of an obvious basin shape. On the other hand, to achieve high and uniform MFD along long stroke and low acceleration waveform harmonic distortion, the requirement for modeling precision is much higher than normal electrodynamic actuators. Currently, there is very little literature on high-accuracy theoretical models for long-stroke CDMCs. He et al. proposed a lumped parameter magnetic equivalent circuit (LPMEC) model [7] to analyze a cylindrical CDMC based on Kirchhoff's law and the superposition theorem. As the reluctance in the yoke is difficult to express, and the MFD varies along the radial direction, the established model is limited. Meanwhile, the strong nonlinear characteristics of yokes and magnetic flux leakage in the CDMC bring extra difficulty in solving model equations [10,[13][14][15][16][17]. Therefore, only qualitative conclusions rather than precise quantitative parameters or characteristic curves were obtained.
In this paper, a patented rectangular closed double magnetic circuit (RCDMC) [18,19] offering ultra-long stroke for ULF vibration exciter is proposed. A high-accuracy theoretical modeling method is proposed through the LPMEC method taking advantage of structural symmetry. The model is solved by iteration calculation and verified through comparison with the finite-element method (FEM). Model-based optimization and experimental verification are carried out. Finally, a vibration exciter with a stroke up to 1.2 m is realized.

Principle of RCDMC
As shown in Figure 1, the proposed structure of RCDMC consists of two PMs, one center yoke, two outer yokes and two end yokes. The two PMs are magnetized horizontally in opposite directions, thus the magnetic poles of the two PMs in contact with the center yoke are of the same polarity, which is the south pole in Figure 1. In this way, in the proposed RCDMC there are only one pole pair (two PMs) and one stator tooth pair, which are protrusive in the middle of outer yokes. Two long air gaps (LAGs) are formed between the stator tooth pair and the center yoke. In this structure, magnetic flux lines can form a set of closed loops through the center yoke, air gaps, outer yokes and one end yoke. Uniform MFD distribution is therefore formed in the two LAGs, along which the Al2O3 ceramic coil skeleton with coil winding wrapped can be moved with driving current flowing in coil winding.

Theoretical Modeling by LPMEC Method
In the 2D planar model as shown in Figure 1b, only the upper half loop is theoretically modeled taking advantage of its symmetry. The RCDMC is modeled assuming that there is no saturation occurring in the yokes, and the magnetic flux leakage outside the RCDMC is ignored. This assumption is proved reasonable by later optimization design results.

LPMEC Model
According to existing research, the magnetic flux leakage inside a CDMC is inevitable, cannot be ignored and is difficult to be quantitatively analyzed [9,13]. In this research, an LPMEC model without considering the magnetic flux leakage is first established, as shown in Figure 2. The long air gap is divided into N equal segments along the axial direction, the center yoke and outer yoke are respectively divided into N + 1 segments along the axial direction: N − 1 equal segments in the middle and two segments left at the two ends.

Theoretical Modeling by LPMEC Method
In the 2D planar model as shown in Figure 1b, only the upper half loop is theoretically modeled taking advantage of its symmetry. The RCDMC is modeled assuming that there is no saturation occurring in the yokes, and the magnetic flux leakage outside the RCDMC is ignored. This assumption is proved reasonable by later optimization design results.

LPMEC Model
According to existing research, the magnetic flux leakage inside a CDMC is inevitable, cannot be ignored and is difficult to be quantitatively analyzed [9,13]. In this research, an LPMEC model without considering the magnetic flux leakage is first established, as shown in Figure 2. The long air gap is divided into N equal segments along the axial direction, the center yoke and outer yoke are respectively divided into N + 1 segments along the axial direction: N − 1 equal segments in the middle and two segments left at the two ends. The magnetic flux leakage is divided into three parts, as shown in Figure 3. The reluctance of part I, Rmf, corresponds to PM leakage flux; the reluctance of part II, Roc, corresponds to outer-tocenter-yoke leakage flux; the reluctance of part III, Rtf, corresponds to tooth-to-center-yoke leakage flux. The Rmf is obtained by modeling magnetic flux flowing in the air gap, as depicted in Figure 4. For the case of hm < Ts + δ, as shown in Figure 4a, the magnetic flux leakage Rmf is modeled with a circular-arc model, the whole magnetic flux leakage by PM, ml φ , is an integral of differential flux flowing through differential width dx, and can be expressed as where F(x) is the magneto-motive force and is expressed as is the reluctance corresponding to width dx, and is expressed as The magnetic flux leakage is divided into three parts, as shown in Figure 3. The reluctance of part I, R mf , corresponds to PM leakage flux; the reluctance of part II, R oc , corresponds to outer-to-center-yoke leakage flux; the reluctance of part III, R tf, corresponds to tooth-to-center-yoke leakage flux.

Theoretical Modeling by LPMEC Method
In the 2D planar model as shown in Figure 1b, only the upper half loop is theoretically modeled taking advantage of its symmetry. The RCDMC is modeled assuming that there is no saturation occurring in the yokes, and the magnetic flux leakage outside the RCDMC is ignored. This assumption is proved reasonable by later optimization design results.

LPMEC Model
According to existing research, the magnetic flux leakage inside a CDMC is inevitable, cannot be ignored and is difficult to be quantitatively analyzed [9,13]. In this research, an LPMEC model without considering the magnetic flux leakage is first established, as shown in Figure 2. The long air gap is divided into N equal segments along the axial direction, the center yoke and outer yoke are respectively divided into N + 1 segments along the axial direction: N − 1 equal segments in the middle and two segments left at the two ends. The magnetic flux leakage is divided into three parts, as shown in Figure 3. The reluctance of part I, Rmf, corresponds to PM leakage flux; the reluctance of part II, Roc, corresponds to outer-tocenter-yoke leakage flux; the reluctance of part III, Rtf, corresponds to tooth-to-center-yoke leakage flux. The Rmf is obtained by modeling magnetic flux flowing in the air gap, as depicted in Figure 4. For the case of hm < Ts + δ, as shown in Figure 4a, the magnetic flux leakage Rmf is modeled with a circular-arc model, the whole magnetic flux leakage by PM, ml φ , is an integral of differential flux flowing through differential width dx, and can be expressed as where F(x) is the magneto-motive force and is expressed as ′ is the reluctance corresponding to width dx, and is expressed as The R mf is obtained by modeling magnetic flux flowing in the air gap, as depicted in Figure 4. For the case of h m < T s + δ, as shown in Figure 4a, the magnetic flux leakage R mf is modeled with a circular-arc model, the whole magnetic flux leakage by PM, φ ml , is an integral of differential flux flowing through differential width dx, and can be expressed as where F(x) is the magneto-motive force and is expressed as F(x) = B r x / µ m ; R (x) is the reluctance corresponding to width dx, and is expressed as R (x)= πx / (2µ 0 Ldx). Appl. Sci. 2020, 10, x FOR PEER REVIEW 4 of 13 For the case of hm > Ts + δ, as shown in Figure 4b, the whole magnetic flux leakage Rmf can be divided into a circular-arc model and a straight-line model.
Thus, ml φ can be calculated by the sum of two integrals The Rtf can be obtained by calculating its corresponding permeance Ptf. The circular-arc straightline model is one of the most satisfactory models for modeling the flux flowing in an air gap as depicted in Figure 5. Ptf is an integral of permeances of differential width dx, the length of each differential element is δ + πx/2. Thus, Ptf can be expressed as Noting the reciprocal relationship between reluctance and permeance, Rtf can be given as  For the case of h m > T s + δ, as shown in Figure 4b, the whole magnetic flux leakage R mf can be divided into a circular-arc model and a straight-line model. R (x) can be expressed as Thus, φ ml can be calculated by the sum of two integrals Finally, the equivalent reluctance R mf can be calculated by dividing the magneto-motive force of PM by the whole flux leakage, i.e., R mf = F m /φ ml .
The R oc can be easily calculated as The R tf can be obtained by calculating its corresponding permeance P tf . The circular-arc straight-line model is one of the most satisfactory models for modeling the flux flowing in an air gap as depicted in Figure 5. P tf is an integral of permeances of differential width dx, the length of each differential element is δ + πx/2. Thus, P tf can be expressed as Appl. Sci. 2020, 10, x FOR PEER REVIEW 5 of 13 Next, an improved LPMEC model taking the magnetic flux leakage in the RCDMC into consideration is obtained by adding Rmf, Roc and Rtf to the LPMEC model in Figure 2, as shown in Figure 6. This improved model can yield accurate analytical results of the MFD in long air gaps. Noting the reciprocal relationship between reluctance and permeance, R tf can be given as Next, an improved LPMEC model taking the magnetic flux leakage in the RCDMC into consideration is obtained by adding R mf , R oc and R tf to the LPMEC model in Figure 2, as shown in Figure 6. This improved model can yield accurate analytical results of the MFD in long air gaps. Next, an improved LPMEC model taking the magnetic flux leakage in the RCDMC into consideration is obtained by adding Rmf, Roc and Rtf to the LPMEC model in Figure 2, as shown in Figure 6. This improved model can yield accurate analytical results of the MFD in long air gaps. The N + 7 closed loops Li (i = 1, 2, … N + 7) are defined with the direction marked in Figure 6. According to Kirchhoff's law, the following equations can be obtained by where, Fm is the magneto-motive force of the PMs, and Fm = Brhm /μm; i φ (i = 1, 2, …, N + 7) is the magnetic flux passing through the ith closed loop; Rm is the reluctance of the PMs; R1i, R3i (i = 0, 1, …, N+1) are the reluctances of the ith axial segment of the center yoke and outer yoke, respectively; R2i (i = 1, 2, …, N) is the reluctance of ith segment of the air gap; R41, R42 are the reluctances of the end yokes. The various reluctances in the improved LPMEC model can be calculated with equations listed in Table 1, in which μy is the permeability of yokes which varies at different positions of yokes. The N + 7 closed loops L i (i = 1, 2, . . . N + 7) are defined with the direction marked in Figure 6. According to Kirchhoff's law, the following equations can be obtained by where, F m is the magneto-motive force of the PMs, and F m = B r h m /µ m ; φ i (i = 1, 2, . . . , N + 7) is the magnetic flux passing through the ith closed loop; R L i (i = 1, 2, ..., N + 7) is the sum of reluctances in the ith closed loop, i.e., is the intersected reluctances of the ith and (i + 1)th closed loops, i.e., Herein R m is the reluctance of the PMs; R 1i , R 3i (i = 0, 1, . . . , N + 1) are the reluctances of the ith axial segment of the center yoke and outer yoke, respectively; R 2i (i = 1, 2, . . . , N) is the reluctance of ith segment of the air gap; R 41 , R 42 are the reluctances of the end yokes.
The various reluctances in the improved LPMEC model can be calculated with equations listed in Table 1, in which µ y is the permeability of yokes which varies at different positions of yokes.

Reluctances Equations Reluctances Equations
Thus, the matrix equation that gives the relation between reluctance, magnetic flux and magneto-motive force is expressed as where R is the reluctance matrix, the dimension of which is (N + 7) × (N + 7), and can be expressed as I is the magnetic flux vector, the dimension of which is N + 7, and can be expressed as ; F is the magneto-motive force vector, the dimension of which is N + 7, and can be expressed as Finally, the MFD in the long air gap can be expressed as

Iteration Solving Method
The solving of the LPMEC model in the form of matrix equation is quite difficult because of the strong nonlinearity of yoke characteristics [13][14][15][16], which is the main limit during the theoretical modeling of magnetic circuits. In order to solve R, I and F and thus to calculate the MFD in the air gap, the key is to determine the relative permeability of pure iron µ r at different positions of yokes. As a result that the yokes are made of pure iron material with strong nonlinear H-µ r curve, direct solving of the model equations is not possible.
In this research, an iterative solving method is proposed to realize the rapid and accurate solution of the matrix equation through iteration calculation. The relative permeability of pure iron µ r as a function of the magnetic field strength H is tested by the National Institute of Metrology (NIM) of China and depicted in Figure 7. The specimen has an outer diameter of 50 mm, an inner diameter of 40 mm and a thickness of 10 mm. The relationship between µ r and H is expressed as follows through polynomial fitting.
where H 0 (µ max r ) is the magnetic field strength when the relative permeability µ r reaches maximum. It can be seen that the polynomial fitting result of the H-µ r curve fits the actual testing results very well.  The procedure of the iterative solving calculation is summarized in Figure 8. First, the dimensions and material property parameters of the RCDMC are initialized. Second, an initial value for the relative permeability of yokes is used, and the reluctance matrix R is calculated. Magnetic fluxes, MFDs and magnetic field strengths in each element of the model are then computed through solving the matrix equation. The values of relative permeability r μ in each element are updated for the next iteration step according to the Hr μ fitting curve expressed as Equation (11). Then the matrix of reluctances and flux densities is recomputed. The procedures above are repeated continually until the difference of MFD between two iteration steps is less than a required precision index (such as 0.1 mT), and the iteration is regarded as converged.  The procedure of the iterative solving calculation is summarized in Figure 8. First, the dimensions and material property parameters of the RCDMC are initialized. Second, an initial value for the relative permeability of yokes is used, and the reluctance matrix R is calculated. Magnetic fluxes, MFDs and magnetic field strengths in each element of the model are then computed through solving the matrix equation. The values of relative permeability µ r in each element are updated for the next iteration step according to the H-µ r fitting curve expressed as Equation (11). Then the matrix of reluctances and flux densities is recomputed. The procedures above are repeated continually until the difference of MFD between two iteration steps is less than a required precision index (such as 0.1 mT), and the iteration is regarded as converged.  The procedure of the iterative solving calculation is summarized in Figure 8. First, the dimensions and material property parameters of the RCDMC are initialized. Second, an initial value for the relative permeability of yokes is used, and the reluctance matrix R is calculated. Magnetic fluxes, MFDs and magnetic field strengths in each element of the model are then computed through solving the matrix equation. The values of relative permeability r μ in each element are updated for the next iteration step according to the Hr μ fitting curve expressed as Equation (11). Then the matrix of reluctances and flux densities is recomputed. The procedures above are repeated continually until the difference of MFD between two iteration steps is less than a required precision index (such as 0.1 mT), and the iteration is regarded as converged.

FEM Verification
The FEM method is used to verify the accuracy of the proposed modeling method. A FEM simulation is conducted using commercial FEM software Maxwell. FEM simulation results and LPMEC model analysis results of MFD in the 1.4 m long air gap, which is depicted as the blue line P1P2 in Figure 1b, are compared.
The results given by FEM and the proposed LPMEC model are shown in Figure 9. The parameters used for comparison are l g = 1400 mm, h m = 50 mm, T m = 150 mm, T y = 75 mm, T s = 25 mm, l se = 150 mm and δ = 20 mm. For the mathematical modeling process of the RCDMC structure, N is set to 1400, which means the 1.4 m air gap is divided equally into 1400 segments with each spaced 1 mm. It can be seen that the results given by the LPMEC model agree with the FEM results very well both in curve shape and values. Only a little difference at the end regions of curves is observed owing to the fringing effect [14]. The calculation results based on the LPMEC model without considering the magnetic flux leakage, which are shown in Figure 2, are also given for comparison. It can be seen that the results agree with the FEM results in curve shape, while the values are~8% higher than the FEM results. Therefore, it can be concluded that the modeling of the magnetic flux leakage in the RCDMC is very successful and brings significant improvement in precision to the LPMEC model.
The FEM method is used to verify the accuracy of the proposed modeling method. A FEM simulation is conducted using commercial FEM software Maxwell. FEM simulation results and LPMEC model analysis results of MFD in the 1.4 m long air gap, which is depicted as the blue line P1P2 in Figure 1b, are compared.
The results given by FEM and the proposed LPMEC model are shown in Figure 9. The parameters used for comparison are lg = 1400 mm, hm = 50 mm, Tm = 150 mm, Ty = 75 mm, Ts = 25 mm, lse= 150 mm and δ = 20 mm. For the mathematical modeling process of the RCDMC structure, N is set to 1400, which means the 1.4 m air gap is divided equally into 1400 segments with each spaced 1 mm. It can be seen that the results given by the LPMEC model agree with the FEM results very well both in curve shape and values. Only a little difference at the end regions of curves is observed owing to the fringing effect [14]. The calculation results based on the LPMEC model without considering the magnetic flux leakage, which are shown in Figure 2, are also given for comparison. It can be seen that the results agree with the FEM results in curve shape, while the values are ~8% higher than the FEM results. Therefore, it can be concluded that the modeling of the magnetic flux leakage in the RCDMC is very successful and brings significant improvement in precision to the LPMEC model. where B(l) is the MFD in the air gap along the x-direction.
The cases with different geometric parameters {hm, Ts, lse, δ} are listed in Table 2. The corresponding results of the FEM and the proposed LPMEC model are listed in Table 3. It can be seen that the deviations of Bm, Bave and Bu between FEM and LPMEC models are all less than 1%. This further proves that the improved LPMEC model has high accuracy, and is valid and suitable for the design and optimization of ultra-long stroke RCDMC.
Further research shows that the deviations increase as Ty decreases, e.g., when Ty = 50 mm, and all other comparison parameters are the same, the deviations of Bm, Bave and Bu are 5%, 6.5%, 10.5%, respectively. The reason is that during the calculation of Rmf, Roc and Rtf in Section 3.1, it is assumed that there is no saturation occurring in the yokes, and the permeability of yokes is far bigger than that of air. When Ty decreases to a certain extent, this assumption is not valid anymore. During the optimization design in Section 5, saturation in yokes is by all means avoided. The B ave and B u are defined as where B(l) is the MFD in the air gap along the x-direction.
The cases with different geometric parameters {h m , T s , l se , δ} are listed in Table 2. The corresponding results of the FEM and the proposed LPMEC model are listed in Table 3. It can be seen that the deviations of B m , B ave and B u between FEM and LPMEC models are all less than 1%. This further proves that the improved LPMEC model has high accuracy, and is valid and suitable for the design and optimization of ultra-long stroke RCDMC.  Further research shows that the deviations increase as T y decreases, e.g., when T y = 50 mm, and all other comparison parameters are the same, the deviations of B m , B ave and B u are 5%, 6.5%, 10.5%, respectively. The reason is that during the calculation of R mf , R oc and R tf in Section 3.1, it is assumed that there is no saturation occurring in the yokes, and the permeability of yokes is far bigger than that of air. When T y decreases to a certain extent, this assumption is not valid anymore. During the optimization design in Section 5, saturation in yokes is by all means avoided.

Optimization Objectives
Since acceleration is proportional to thrust force and MFD in the air gap, the value of average MFD along the long air gap, B ave , is the first optimization objective.
The non-uniformity of the MFD is the main contribution to the AWHD in a long stroke CDMC. As a result that the maximum acceleration generated by a CDMC is limited by its stroke, during the operation of a CDMC, a standard vibration waveform is always generated taking the middle point of the air gap along the x-direction as the original point. Therefore, the maximum vibration amplitude of the standard vibration waveform generated is (l g -l c )/2, wherein l g is the length of the air gap, l c is the length of the driving coil. According to Ampere law and the relationship between acceleration and displacement, the acceleration waveform with maximum amplitude can be expressed as where k(f ) equals 0.5(2π f ) 2 l g − l c /B ave when the vibration frequency is f, B(l) is the MFD in the air gap along the x-direction and θ is the phase of the standard vibration. The harmonic component of a can be calculated through spectrum analysis. Thus, the AWHD λ can be calculated by where a i is the magnitude of the ith harmonic component of the acceleration waveform.

Optimization of h m and l se
Dimensional optimization is also carried out to obtain large and uniformly distributed MFD as well as low AWHD. As shown in Figure 10, it can be seen that the AWHD has little variation as h m and l se change, while the average MFD along the air gap increases with a decreasing slope and decreases with a constant slope as h m and l se get larger, respectively. h m should be set as big as possible to increase the average MFD. However, taking into account the manufacturing capability of large block NdFeB, h m = 50 mm is set.
decreases with a constant slope as hm and lse get larger, respectively. hm should be set as big as possible to increase the average MFD. However, taking into account the manufacturing capability of large block NdFeB, hm = 50 mm is set.
The increasing of lse may decrease Roc according to Equation (4), and thus increase the magnetic flux leakage and decrease Bave. During the design of the RCDMC, four supporting blocks with a length of (lse -hm) in the x-direction and height of (Ts + δ) in the y-direction are assembled beside the ends of the stator teeth. The length of supporting blocks is set to 100 mm to support the mass of center yoke and outer yokes. Therefore, lse = 150 mm is set.

Optimization of δ and Ts
A performance-topology-based approach has been employed to conduct optimization of δ and Ts. Six hundred values of Bave and λ are calculated based on the LPMEC model to map the performance parameters of Bave and λ to the values of δ and Ts. The performance topology is charted as shown in Figure 11. The final selected design point is {δ = 20 mm, Ts = 25 mm} to balance between high Bave and low λ. After optimization, an MFD of 122 mT and an AWHD of 0.45% along 1.2 m stroke are achieved.

Optimization Design Results
Among the design parameters of the RCDMC, some are preset values due to the following design limitations: (1) Air gap length lg is determined by the desired length of stroke; (2) dimension in the z-direction, L, is determined by the desired thrust force; (3) dimension Tm is limited by the The increasing of l se may decrease R oc according to Equation (4), and thus increase the magnetic flux leakage and decrease B ave . During the design of the RCDMC, four supporting blocks with a length of (l se -h m ) in the x-direction and height of (T s + δ) in the y-direction are assembled beside the ends of the stator teeth. The length of supporting blocks is set to 100 mm to support the mass of center yoke and outer yokes. Therefore, l se = 150 mm is set.

Optimization of δ and T s
A performance-topology-based approach has been employed to conduct optimization of δ and T s . Six hundred values of B ave and λ are calculated based on the LPMEC model to map the performance parameters of B ave and λ to the values of δ and T s . The performance topology is charted as shown in Figure 11. The final selected design point is {δ = 20 mm, T s = 25 mm} to balance between high B ave and low λ. After optimization, an MFD of 122 mT and an AWHD of 0.45% along 1.2 m stroke are achieved.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 10 of 13 decreases with a constant slope as hm and lse get larger, respectively. hm should be set as big as possible to increase the average MFD. However, taking into account the manufacturing capability of large block NdFeB, hm = 50 mm is set. The increasing of lse may decrease Roc according to Equation (4), and thus increase the magnetic flux leakage and decrease Bave. During the design of the RCDMC, four supporting blocks with a length of (lse -hm) in the x-direction and height of (Ts + δ) in the y-direction are assembled beside the ends of the stator teeth. The length of supporting blocks is set to 100 mm to support the mass of center yoke and outer yokes. Therefore, lse = 150 mm is set.

Optimization of δ and Ts
A performance-topology-based approach has been employed to conduct optimization of δ and Ts. Six hundred values of Bave and λ are calculated based on the LPMEC model to map the performance parameters of Bave and λ to the values of δ and Ts. The performance topology is charted as shown in Figure 11. The final selected design point is {δ = 20 mm, Ts = 25 mm} to balance between high Bave and low λ. After optimization, an MFD of 122 mT and an AWHD of 0.45% along 1.2 m stroke are achieved.

Optimization Design Results
Among the design parameters of the RCDMC, some are preset values due to the following design limitations: (1) Air gap length lg is determined by the desired length of stroke; (2) dimension in the z-direction, L, is determined by the desired thrust force; (3) dimension Tm is limited by the

Optimization Design Results
Among the design parameters of the RCDMC, some are preset values due to the following design limitations: (1) Air gap length l g is determined by the desired length of stroke; (2) dimension in the z-direction, L, is determined by the desired thrust force; (3) dimension T m is limited by the manufacturing capability of large block NdFeB magnets. Dimension T y is set to equal half the width of magnet T m , so that the magnetic circuits in the center yoke and "necks" of two outer yokes reach the same saturation level. The rest set of variables, {h m , T s , l se , δ}, is determined through optimization.
After model-based optimization, theoretically, an MFD of 122 mT and an AWHD of 0.45% along 1.2 m stroke can be achieved.

Experimental Setup
As shown in Figure 12, a well developed prototype of the optimized RCDMC has been assembled in a horizontal ULF vibration exciter. It should be noted that all the compenents, including the center yoke, outer yokes and end yokes have been vacuum annealed to guarantee the great permeability of yokes. PMs are assembled into the prototype with specially-designed assembly tools after ther are finish grinded and magnetized. The moving subassembly guided by a self-designed air bearing can slide freely along the air gap with coil winding, coil skeleton and working platform assembled. During the experiments, a GM55 gaussmeter (Shanghai Tindun Industry Co., Ltd, Shanghai, China) with a range of 200 mT, a resolution of 0.01 mT and an accuracy of ± 2.0% (DC) is used.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 11 of 13 manufacturing capability of large block NdFeB magnets. Dimension Ty is set to equal half the width of magnet Tm, so that the magnetic circuits in the center yoke and "necks" of two outer yokes reach the same saturation level. The rest set of variables, {hm, Ts, lse, δ}, is determined through optimization. After model-based optimization, theoretically, an MFD of 122 mT and an AWHD of 0.45% along 1.2 m stroke can be achieved.

Experimental Setup
As shown in Figure 12, a well developed prototype of the optimized RCDMC has been assembled in a horizontal ULF vibration exciter. It should be noted that all the compenents, including the center yoke, outer yokes and end yokes have been vacuum annealed to guarantee the great permeability of yokes. PMs are assembled into the prototype with specially-designed assembly tools after ther are finish grinded and magnetized. The moving subassembly guided by a self-designed air bearing can slide freely along the air gap with coil winding, coil skeleton and working platform assembled. During the experiments, a GM55 gaussmeter (Shanghai Tindun Industry Co., Ltd, Shanghai, China) with a range of 200 mT, a resolution of 0.01 mT and an accuracy of ± 2.0% (DC) is used. Figure 12. Prototype of RCDMC assembled in a horizontal ULF vibration exciter. Reproduced with permission from [19], IEEE, 2019.

MFD along Air Gaps
To measure the MFD along the air gaps, the sensitive part of the gaussmeter probe is placed in the air gap with the meter fixed on the moving subassembly. In this way, the MFD along the air gap can be sampled at equal intervals as the moving subassembly is controlled to travel slowly from one end to the other.
Measurement results of MFD along the upper and lower air gaps are shown in Figure 13. Compared with the results shown in Figure 9, the measured MFD is about 20% smaller than the theoretical model and FEM results. The main reason is that the remanence value of PM material is not exactly the same as assumed in the models, and magnetic flux leakage at the flank sides of the RCDMC is inevitable. Benefiting from the strict vacuum anneal process [20], the actual permeability of yokes is higher than the tested result of the specimen part in Figure 7, and the measurement result of uniformity of MFD is better than the theoretical model and FEM results. Due to the fringe effect, the MFD increases rapidly in the regions of 20 mm from both ends of the air gaps, which is consistent with the FEM results. Differing from the theoretical and FEM analysis results, the two MFD measurement curves along the upper and lower air gap respectively are not perfectly symmetrical, and there are differences between the two curves. This is due to the unequal thickness of the two air gaps in the prototype caused by machining and assembly errors, while material distortions also have

MFD along Air Gaps
To measure the MFD along the air gaps, the sensitive part of the gaussmeter probe is placed in the air gap with the meter fixed on the moving subassembly. In this way, the MFD along the air gap can be sampled at equal intervals as the moving subassembly is controlled to travel slowly from one end to the other.
Measurement results of MFD along the upper and lower air gaps are shown in Figure 13. Compared with the results shown in Figure 9, the measured MFD is about 20% smaller than the theoretical model and FEM results. The main reason is that the remanence value of PM material is not exactly the same as assumed in the models, and magnetic flux leakage at the flank sides of the RCDMC is inevitable. Benefiting from the strict vacuum anneal process [20], the actual permeability of yokes is higher than the tested result of the specimen part in Figure 7, and the measurement result of uniformity of MFD is better than the theoretical model and FEM results. Due to the fringe effect, the MFD increases rapidly in the regions of 20 mm from both ends of the air gaps, which is consistent with the FEM results. Differing from the theoretical and FEM analysis results, the two MFD measurement curves along the upper and lower air gap respectively are not perfectly symmetrical, and there are differences between the two curves. This is due to the unequal thickness of the two air gaps in the prototype caused by machining and assembly errors, while material distortions also have a major impact on that. The average MFD along the upper and lower air gaps is 102.44 mT and 102.79 mT, the non-uniformity of MFD is 0.21% and 0.26%, and AWHD is 0.24% and 0.27%, respectively.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 12 of 13 a major impact on that. The average MFD along the upper and lower air gaps is 102.44 mT and 102.79 mT, the non-uniformity of MFD is 0.21% and 0.26%, and AWHD is 0.24% and 0.27%, respectively.

Conclusions
To theoretically analyze the RCDMC offering ultra-long stroke up to 1.2 m for ULF vibration exciters, a high-accuracy theoretical modeling method is put forward through the LPMEC method, taking advantage of structural symmetry as well as easy reluctance lumping and expression of the RCDMC. The matrix model equation is solved by iteration calculation to deal with the strong nonlinearity of the yoke material and verified using FEM as a comparison reference. The deviations between the proposed LPMEC model and FEM results are less than 1%, proving the high accuracy of the model. After model-based optimization is conducted, an average MFD of 122 mT and an AWHD of 0.45% along 1.2 m stroke are achieved theoretically, while an average MFD of 102 mT and an AWHD of 0.27% along 1.2 m stroke are realized in the accordingly developed prototype. We can come to a conclusion that the proposed RCDMC and modeling method provide a good solution for high-performance ULF vibration exciters. Our future work will focus on the evaluation and integration of the vibration exciter in a 0.01 Hz ULF vibration calibration system.

Conclusions
To theoretically analyze the RCDMC offering ultra-long stroke up to 1.2 m for ULF vibration exciters, a high-accuracy theoretical modeling method is put forward through the LPMEC method, taking advantage of structural symmetry as well as easy reluctance lumping and expression of the RCDMC. The matrix model equation is solved by iteration calculation to deal with the strong nonlinearity of the yoke material and verified using FEM as a comparison reference. The deviations between the proposed LPMEC model and FEM results are less than 1%, proving the high accuracy of the model. After model-based optimization is conducted, an average MFD of 122 mT and an AWHD of 0.45% along 1.2 m stroke are achieved theoretically, while an average MFD of 102 mT and an AWHD of 0.27% along 1.2 m stroke are realized in the accordingly developed prototype. We can come to a conclusion that the proposed RCDMC and modeling method provide a good solution for high-performance ULF vibration exciters. Our future work will focus on the evaluation and integration of the vibration exciter in a 0.01 Hz ULF vibration calibration system.