Abstract
The constant-frequency double-rotor generator has the potential to be used in the marine energy conversions, such as tides and marine current power generation system due to the advantages. However, there are two different rotational speed movement boundaries in the constant-frequency double-rotor generator, especially the inner moving part boundary between the inner wound rotor and outer permanent magnet rotor. This complicates the finite element analysis and greatly increases the analysis time. In this article, a new finite element analysis method based on double boundaries interpolation multiple reference frame is proposed. The principle of the interpolation multiple reference frame method is introduced, and the analysis and comparison with the traditional method are carried out. The analysis results show that the interpolation multiple reference frame method simplifies the preprocessing procedure, reduces computation time, ensures high calculation accuracy, and has the potential to be applied in all kinds of multiple mechanical port machines.
Keywords
Introduction
The constant-frequency double-rotor (CFDR) generator is a kind of variable-speed constant-frequency (VSCF) generator, which may be used in the application of ocean energy harvesting. 1 The CFDR generator outputs a fixed frequency electrical power, so it can be directly integrated with the grid.2,3 The CFDR generator has the advantages of no gearbox, high efficiency, and flexible regulation of active and reactive power. Then the CFDR generator is very suitable for large-scale conversion of marine energy.
The finite element method (FEM) is one of the most commonly used methods in electrical machine design and analysis. The static as well as the transient FEM are commonly used to analyze the performance of the electric machine. The time-stepping FEM based on field-circuit-mechanical equation coupling is one of the most effective transient FEMs. In electrical machines, there is a relative movement between the stator and rotor. The meshing must be executed in every step of the finite element analysis, which will take a lot of computing time. And the shape of the new mesh will be distorted, which may lead to reduced calculation accuracy. Hence, how to deal with the motion problem of the electrical machine is very important.
In previous literatures, there are several effective methods to treat the movement of the electrical machines. The conception of an air-gap element for the dynamic analysis of the electromagnetic field in electric machines is proposed in Abdel-Razek et al. 4 In this air-gap element method, the entire air-gap region is one macro element and its field is calculated analytically using a Fourier series expansion of the vector potential. The movement is accomplished by a simple recalculation of the Fourier series without any modification of the mesh structure, resulting in simple and efficient time-stepping.5–7 The results can be very accurate because of the high order of the field representation in the air-gap region. 8 But for the insertion of dense matrix blocks into the sparse finite element matrices, the computational time is substantially increased.9,10 The moving band method is proposed in Davat et al. 11 There is no need of the special elements or coupling techniques; no dense blocks are generated in the system matrix and the system dimension is not increased. 12 For this reason, the moving band method should be superior to the air-gap element method in terms of computational speed. However, re-meshing the air-gap region is inevitable, and thus the numbering as well as the amount of nodes in the mesh are not constant. The conditioning of the system matrix is not maintained and preprocessing must be rerun when the mesh changes. Also, because the elements in the air-gap are geometrically distorted or re-meshed to accommodate arbitrary movement, the results often have an oscillating error component. 5 It has been shown that this problem can be drastically reduced by using higher order elements in the moving band. 13 However, the use of higher order elements increases the density of the element, and the computational time will substantial increase as the case of the air-gap element method.
A hybrid finite element–boundary integral method is proposed in Salon and Schneider. 14 Using this method, users can represent any region by either finite elements or boundary elements. There are many advantages, such as the number of equations can be reduced by using boundary, open boundary problems are easily formulated, and non-linear regions can be represented by finite elements. 15 However, this method increases the bandwidth of the system matrix and increases the computational time. An interpolation slip surface method is proposed in Yan et al. 16 It effectively treats the single movement boundary. But for the multiple mechanical port machines, there are different speed movement boundaries, which complicate the preprocessing procedure and increase the analysis time. In this article, a new finite element analysis method based on double boundaries interpolation multiple reference frame (DBIMRF) method suitable for the CFDR generator is proposed. The proposed method simplifies the preprocessing procedure, reduces computation time, ensures high calculation accuracy, and has the potential to be applied in all kinds of multiple mechanical port machines. The DBIMRF method could be implemented in ANSYS by ANSYS Parametric Design Language. Through the combination of the IMRF method and ANSYS, the simulation time may be greatly decreased.
DBIMRF method
During the design of the electrical machines, the simplest way to analyze the static electromagnetic field is re-meshing the machine domain in every step of the movement, and it can obtain the most accurate calculation result. However, this may greatly increase the simulation time. The moving band method is commonly used in the transient electromagnetic field simulation, which only needs re-meshing the air-gap between the rotor and the stator in every step, so that it can reduce the simulation time, but the re-meshing will produce distorted elements and the distorted element will lead to the calculation accuracy decreasing.
The slip surface method
The slip surface method is proposed in Yan et al. 16 It does not need to re-mesh and would not change the shape of the element. So it ensures the calculation accuracy, greatly reduces the simulation time, and it is easy to realize. 17
The theory of the slip surface method is as follows. First, two independent coordinate systems, namely, a static coordinate system and a rotating coordinate system, are used for the stator and the rotor, respectively. Second, the electromagnetic field equations of the stator and the rotor are set up, respectively. Then the equations of stationary parts and rotating parts are solved in parallel through the movement boundary condition. Finally, the whole field solution can be obtained. The specific implementation tips are presented as shown in Figure 1.

Specific implementation tips of slip surface method.
As shown in Figure 1, a movement boundary that consists of two overlapping arcs belonging to the stator and the rotor, respectively, is set up inside the air-gap between the stator and the rotor. There are n equidistant nodes on each arc, and each arc is divided into n − 1 segments. In the initial position, the nodes on the two arcs coincide. The radian of each curve segment or its integral multiple is equal to the rotor rotation step size. So the nodes belonging to the stator and the rotor on the movement boundary are overlapped in every step and the meshes of the stator and the rotor do not need to be changed.
The slip surface method can reduce computation time and ensure high calculation accuracy. In the design of the CFDR generator, this method can improve the efficiency and accuracy very well, but for the transient performance simulation, it is not suitable.
Figure 2 shows the structure of the CFDR generator. It consists of three parts: wound stator, inner wound rotor, and outer permanent magnet (PM) rotor. The freely rotating outer rotor is located between the stator and the inner rotor. There are two air-gaps and two movement boundaries in CFDR generator where the inner movement boundary is between two rotating parts. For the different rotational speed of two movement boundaries, the time step size is not easy to determine. The nodes on the boundaries of slip surface method might be too dense in analyzing the CFDR generator. It will lead to longer calculation time, lower calculation accuracy, and even not convergence. So the slip surface method must be improved for the analysis of CFDR generator.

Structure of CFDR generator.
The DBIMRF method
The first-order linear triangle-shaped element is selected in the finite element model. The triangle-shaped element has six nodes, as shown in Figure 3. The vector magnetic potential (AZ) of the node m is equal to (AZi + AZk)/2. According to the principle of finite element, the AZ of the point on the side of element satisfies linear equation, and the AZ of node on the movement boundary can be similar to the AZ of the point on the side of element. So the Lagrangian interpolation method can be used to deal with the movement boundary and do not need to request the space step length and node distance matching.

Triangle-shaped element.
Figure 4 shows that there are two movement boundaries, L1 and L2. For the existence of L2 which is between two rotating parts, the difficulty of dealing with the boundary problem is increased. In order to reduce the complexity and difficulty of programming, improve the portability of the program; double boundaries method is used.

Movement boundaries and boundary conditions of CFDR generator.
The implementation of DBIMRF method is as follows. Three independent coordinate systems, namely, the Euler coordinate system (static coordinate system) and two Lagrangian coordinate systems (servo) are set up for the stator, inner rotor, and outer rotor, respectively. The magnetic fields of the three parts are solved, respectively. And the magnetic fields are coupled through the vector magnetic potential (AZ) of the nodes on the movement boundaries. The nodes, boundary conditions, and movement boundaries of the CFDR generator in the IMRF method are shown in Figure 4. The moving boundary L1 is composed of two overlapping circular arcs, which belong to the stator and the rotor, respectively. Lagrange interpolation coupling method is used to couple AZ of the nodes by periodic boundary condition, as
where K1 and K2 are the Lagrangian interpolation coefficient. So it is easy in preprocessing where the space step length has no relationship with the node distance.
Not same as L1, the boundary L2 between the outer and the inner rotor is not treated as the active boundary. Using double moving boundary method, two movement boundaries, L3 and L4, are defined. Then the air-gap between the outer and the inner rotor is divided into three parts. The part between L3 and L4 is treated as a static part where the boundary of the moving part and the static part is realized. In this case, L3 and L4 may be treated in the way as L1. Using this method, the finite element preprocessing procedure is greatly simplified.
The DBIMRF method is the improvement of the slip surface method in the application of the transient electromagnetic field analysis. The distance between the nodes is determined by the space step length as the CFDR generator running in the steady state. So the high accuracy of the CFDR static performance analysis and the effectiveness of transient performance analysis all may be achieved.
Element stiffness matrix
The shape and size of the element are not changed in every step by using the DBIMRF method. Because the shape and size of the element are not changed, the element stiffness matrix will not be changed. So the calculation accuracy will not be reduced and it can be proved as follows.
Figure 3 shows that the numbers of three vertices of the triangle-shaped element are i, j, and k, which are arranged in counterclockwise direction. And their coordinates are (xi, yi), (xj, yj), and (xk, yk), respectively. The side lengths of the three sides are li, lj, and lk, respectively
where Se is the area of the element and µ is the permeability of the element material. Then the element stiffness matrix can be obtained as equation (5)
Similarly, the other elements of the element stiffness matrix can be obtained
From equations (6)–(9), it can be seen that the element stiffness matrix is only dependent on the shape and size of the element, no matter with the position and rotation angle of the element. So, using the slip surface and DBIMRF method, the element stiffness matrix will not be changed in every step. So the global stiff matrix will not be changed and the high calculation accuracy can be achieved.
Electromagnetic field analysis
The static electromagnetic field analysis
In order to verify the effectiveness of DBIMRF method, the IMRF method in static electromagnetic field simulation is used for analyzing a specific CFDR generator model. The main dimensions of the CFDR generator are shown in Table 1. Considering the symmetry of CFDR generator structure, the geometrical model of half CFDR generator is established as shown in Figure 5, where the material properties are displayed by different colors.
The main dimensions of the CFDR generator.
OD: outer diameter; ID: inner diameter.

Geometrical model of CFDR generator.
The accuracy of FEM largely depends on the level of meshing. For achieving high accuracy calculation, high level of meshing for the CFDR generator is adopted, as shown in Figure 6. Then, according to Figure 4, the boundary conditions are implemented on the finite element model, as shown in Figure 7.

Mesh model of CFDR generator.

Finite element model of CFDR generator with boundary conditions.
The flux distributions and magnetic flux density of the CFDR generator only under PM excitation are shown in Figure 8. It is shown that most of the PM flux goes through the stator and inner rotor; meanwhile, the magnetic field of the corresponding stator and rotor tooth is saturated, and the excitation of the generator is provided by the PM. The result is consistent with the actual situation.

Magnetic field of CFDR generator: (a) magnetic flux distribution in CFDR generator and (b) magnetic flux density in CFDR generator.
The calculation accuracy and time of the IMRF method are studied and compared with the re-meshing method. The same number and shape of element are selected. Here, the number of the element is 48,363. For comparing fairly, two methods are carried out on the same computer, which has a configuration as shown in Table 2.
Computer configuration.
The un-alignment of boundary nodes may influence the calculation accuracy of the IMRF method. In order to analyze calculation accuracy of the IMRF method in the disadvantage conditions, the rotor rotation angle is set not equal to the integer times of node spacing. In the static electromagnetic field analysis, only the no-load electromotive force (EMF) of the inner rotor is studied. So the inner rotor is seen as a static part, and only the outer rotor rotates. Using the IMRF method and the re-meshing method, the no-load EMF at rated speed can be obtained, respectively. Figure 9 shows the no-load EMF waveforms of phase A winding of the inner rotor. It can be seen that the result agrees very well. Figure 10 shows the relative error caused by the DBIMRF method. The relative error is calculated based on the result of the re-meshing method as truth value. According to Figure 10, it can be seen that the impact caused by nodes un-alignment is very small and can be ignored in the DBIMRF method. So the IMRF method can keep the high calculation accuracy.

No-load EMF of phase A winding of the inner rotor.

Relative error caused by the double boundaries IMRF method.
There are 90 steps which are set up in analysis, and the calculation times by using IMRF method and re-meshing method are 42.5 and 68 min, respectively. It is shown that the DBIMRF method may greatly reduce the calculation time. And with the step number increasing, the proportion of the saved time will be increased.
Time-stepping FEM simulation
In this section, the DBIMRF method is used in the time-stepping FEM to analysis the no-load start-up process of the CFDR generator. The electromagnetic field and the mechanical movement equation should be coupled during the calculation. There are two coupling methods: direct coupling method and indirect coupling method.18,19 The direct method directly plugs the angular velocity into the electromagnetic field finite element equation (10), which will result in asymmetrical matrix and lead to the difficulty of solving finite element equations. So the indirect coupling method is adopted in this article
where µ is the reluctivity, A is the magnetic vector potential, V is the velocity, and J is the current density in the coils.
The basic idea of indirect coupling method is as follows. The position of the rotor at kth step is replaced by the position of the rotor at (k − 1)th step. The coefficient matrix of the field-circuit coupling finite element equations will keep constant during the solution. After the kth step calculation of the electromagnetic field is finished, the torque can be obtained and the rotor position of kth step will be obtained by combining the rotor mechanical motion equation. Then the continued calculation of (k+1)th step might be carried out in the same way. By this method, the solution of the field-circuit-mechanical motion coupled equations will be obtained.
The mechanical motion equation of the outer PM rotor is
where the Ω is the rotating speed of the outer PM rotor, θ is the position of the PM rotor, Jm is the rotational inertia of the PM rotor, Tem is the electromagnetic torque, and TL is the load torque.
When the space step size Δθ is fixed, forward-difference (explicit Euler method) is used for the discretization of the mechanical motion equation (12). Assuming that the rotor rotates at a constant velocity in every step, the mechanical motion equation is
Fixed space step size is suitable for slip surface method, but it is not suitable for analyzing the start-up process of CFDR generator, since the generator runs in a wide range of speed. The transient process in a short time is concerned in transient analysis, and the fixed space step size limited the time step, that is inconvenient.
When the time step size is fixed, backward-difference (implicit Euler method) is used for the discretization of the mechanical motion equation (11). Assuming that the rotor rotates as uniformly accelerated motion in every step, the mechanical motion equation can be expressed as
After knowing the electromagnetic torque, the current speed Ω(n) can be calculated by equation (13), and then the Δθ(n), which is the space step size of the next step, can be obtained. The position of the rotor can also be calculated. Through this method, the dynamic change in electromagnetic torque and current of the CFDR generator may be calculated.
The DBIMRF method adopts the fixed time step size, which overcomes the difficulty in selection of the time step size and does not cause any increase in the calculation time and the complexity of the pre-treatment procedure. This effectively guarantees the calculation accuracy. The analysis result of the CFDR generator start-up process is shown in Figure 11. For verifying correctness of the result, the result is compared with high-order moving band method, which has been proved highly accurate. It can be seen from Figure 11 that the accuracies of the two methods are closing. So the IMRF method guarantees the calculation accuracy. And according to Figure 11, the transient stability of the CFDR generator is well where the speed, torque, and current only have slight perturbation around the ideal value. The CFDR generator has the potential to be applied in ocean energy harvesting.

Analysis result of start-up process: (a) out rotor speed, (b) inner rotor current, and (c) torque.
Conclusion
The CFDR generator has two rotors with different rotation speeds and brings difficulties to the electromagnetic field finite element analysis. To deal with this problem, a DBIMRF method is introduced and studied in this article. Due to cancellation of the constraint that the nodes on the movement boundaries must be constantly overlapped, a fixed time step size is adopted without considering the space step size. The space step size is obtained by calculation. It is verified that the method solves the difficulty in choosing the time step without increasing the density of the elements, greatly shortens the calculation time, and effectively guarantees the calculation accuracy.
Footnotes
Academic Editor: Fabrizio Marignetti
Declaration of conflicting interests
The authors declare that they have no conflict of interest.
Funding
This work was supported, in part, by the National Key Basic Research Program of China (973 Program of China) (2013CB035603) and the National Natural Science Foundation of China (51137001, 51320105002).
