Development of a generic framework for lumped parameter modeling

2.1 Mathematical model

The geometry is discretized in a number of bidirectional elementary reluctance blocks as shown in Figure 1. R _x1, R _x2, and R _y1, R _y2 are the reluctances on the x- and y-axis branches of the block element, respectively. Compared to previous studies, elements are allowed to be attached with two kinds of materials, which enables establishing a more precise model for sloping sides in machines.

where T _x and T _y are shape functions in the x- and y-direction, determined by the shape and size of elements. μ _i and μ _j are different permeabilities assigned to the element.

Figure 1

Reluctance block elements.

The circuit solver is developed based on Kirchhoff’s laws [11]. It can be established for a network containing N _n nodes and N _b branches as follows:

where Φ _B is the vector of the fluxes on branches, Φ B s is the vector of the fluxes generated by current sources, such as a PM and current. C is a connectivity matrix of N _n × N _b size, and the item C _ij is defined as follows:

According to Ohm’s law, Φ _B and Φ B s can be written as:

where U, P, and F are the voltage vector, the diagonal matrix of permeance, and the magneto-motive force (MMF) source vector, respectively. And the conversion between the voltage U and the node potential V can be expressed as:

Consequently, equation (2) is transformed into

MMF sources due to the PM need to be placed on the branches of the magnetization direction in the elements through all layers of PM areas, where each source is calculated as shown in equation (7) [12]. H _c is the PM magnetic coercivity, and l _mag is the length along the magnetization direction for each element (Figure 2).

To obey Ampere’s circuital law, MMF sources due to current need to model the distribution of MMF in the air gap first, which is a function of spatial position according to the winding distribution, denoted as F _g(x,t). φ _A(x), φ _B(x), and φ _C(x) are defined as the functions of MMF distribution in the air gap, when the unit current is provided for each phase. x is the x-axis coordinate of each element. Combining the spatial position function with the three-phase input current, we have

where F _A, F _B, and F _C are the functions of MMF distribution changing with time. By summing them up, the expression of F _g(x,t) is

Figure 2

Modeling of the PM unit.

Figure 3 shows the MMF distribution in the air gap, when setting t = 0 s, θ = 0°. Then, the MMF value contributed by each element should be calculated using the proportion to the whole current area at the corresponding point.

where h _w is the height of the winding area and l _l,e is the length along the magnetization direction for each element.

Figure 3

MMF distribution in the air gap.

The main principle of motion handling is to avoid regenerating the mesh in the studied domain with no or little compromise in precision. An interpolation method is employed here, based on the continuity of the scalar magnetic potential and the continuity of the flux density normal component at the sliding interface (Figure 4).

The subscripts s and r are used to denote the values from the stator side and the rotor side. Equation (11) shows a linear interpolation on the interface, and for the interpolation of k orders, the formulation is as follows:

By using the projection against a specific direction, we can obtain the normal component of the flux density of interface nodes from the connecting branch components.

Figure 4

Interface between the rotor and the stator.

Apart from the sliding interface, the RN is also defined by the Dirichlet boundary condition and the periodic or anti-periodic condition (Figure 5). For node n _Dk , we have

For node i and node j, no matter what kind of boundary condition is assigned, they are viewed as adjacent elements, so they are connected by a branch, denoted as b. When using the periodic condition, we create the connectivity matrix as per the definition in equation (3). When using the anti-periodic condition, element C _jb in the connectivity matrix is set to 1 rather than −1.

Figure 5

Boundary condition.

The Newton–Raphson method is a technique used for solving nonlinear equations numerically [13]. This method uses an iterative process to approach one root of a function, and this root depends on the initial value. For equation (2), we can write

knowing that

where operator ∘ is the element-wise product.

Branch b connects node i and node j, and by considering each part in two adjacent elements of a unique branch separately, the derivative of permeance with respect to node potentials is obtained from the corresponding components based on the expression of permeance calculation.

where l _b is the length of branch b. Note that the Dirichlet boundary condition and interpolation formulation at the sliding interface modify some values in the coefficient matrix of linear equations, but they have no effect on the calculation of the Jacobian matrix. Accordingly, these items are set to zero after getting the Jacobian matrix using equation (13).

2.2 Tool structure

The framework is developed in four layers to provide genericity, flexibility, and extensibility (Figure 7). The frontal layer is a designing part serviced by widely used mesh tools. Apart from commercial software like ANSYS^® and FLUXTM, GMSH is an open-source toolkit with a powerful and fast mesh function. In this study, an Excel worksheet constructed in a specific format that describes the structure of the machine is used as an input file (Figure 6).

Figure 6

Format of describing the structure of the machine in EXCEL.

Figure 7

Implementation of the tool.

Then, the framework reads the generated file into the circuit generator layer using an adaptor, which decodes the mesh information according to the pre-setting rules. The Excel worksheet in Figure 6 contains dimensions and confined mesh instructs for ExcelAdaptor to decide the concrete mesh scheme. DomainManager helps to manage the geometry properties of the studied domain and generate the size data for each element. Coupled with object Winding and object MaterialManager, DomainManager is processed by object CircuitGenerator, which provides the circuit descriptor encapsulating detailed information of branches that are required in object Assembler during iterations. MaterialManager offers the function of maintaining properties defined by users and helps to manage different types of materials especially nonlinear materials. Also, an interface is opened to accept the customer defined B–H curve fitting function. By default, the program uses an analytic saturation curve with knee adjustment.

where

The implementation of object Winding depends on different configurations of winding, which means that distributed windings and concentrated windings ought to be modeled in two ways. CircuitGenerator calls the unified interface to get the function F _g(x,t) in the air gap when producing an MMF distribution model in the initial state to accelerate the solving process.

Figure 8 shows the main process of simulation at each time step in the circuit solver layer, when users launch a new analysis. Due to the nonlinear properties of materials, the circuit solver layer is developed in two types of iterative methods: fixed-point method and Newton–Raphson method. The Newton–Raphson method approaches a sufficiently precise value much faster than the fixed-point method because of using the tangent of the graph, but sometimes it is unable to reach a convergence as stable as the others. So in the program, both methods are provided. The construction of the Jacobian matrix has been described in the last section. The core solvers of fields are divided in single modules so that users are capable of building their own solver on these stable functions.

Figure 8

Process of the circuit solver layer.

The circuit solver layer inquires the upper layer to obtain the updated magnetic state through the specified interface at each step, which ensures the decoupling between two layers so that users are capable of designing their own modules by utilizing the well-encapsulated modules. Overloading of the default motion and boundary processing function is supported in the program, so a more effective motion handling method like the hybrid analytical method, which performs better in calculating the magnetic field in the air gap, is obtained.

All fields output by the circuit solver layer are registered in the DomainManager and exported to the CSV format, a type of file supported by many programs for post-processing.

3.1 Structure presentation

The structure in Figure 9 is a classical 12/10 FSPM machine with concentrated windings, which has been analyzed in many previous cases. Although the dimensions of the machine were chosen to simplify the modeling study, the values can be set as any number in actual application. Table 1 gives the main machine characteristics.

Figure 9

Studied structure.

3.2 Validation study

The comparisons shown in Figures 10–14 are between the FEM model meshed with 18,756 nodes and the RN model meshed with 4,800 nodes. Local quantities, such as the tangential and normal component of air gap flux density (Figures 10 and 12); global quantities, such as flux linkage and induced voltage versus time (Figures 11 and 13); torque; and iron loss (Figure 14) are evaluated by the classical method used in FEM. Flux linkage is the sum of fluxes (in the r-axis) passing through all block elements located at slot mid-height spanning one coil pitch. The electromotive force is obtained by the differentiation of the phase flux. The Maxwell stress tensor method is used for the computation of torque. For iron loss estimation, the experimentally verified iron loss model is used here [14]:

where P _h, P _c, and P _e are, respectively, static hysteresis loss, classical eddy current loss, and excess loss; k _h, k _c, and k _e are the coefficients of the corresponding loss component; f and B _m are the frequency and amplitude of the fundamental flux density; and β is an empirical parameter obtained from experimental measurement.

Figure 10

Air-gap magnetic flux density component comparison, I _m = 0 A. (a) Normal component and (b) tangential component.

Figure 11

Flux linkage and EMF comparison, I _m = 0 A. (a) Flux linkage and (b) induced voltage.

Figure 12

Air-gap magnetic flux density component comparison, I _m = 10 A. (a) Normal component and (b) tangential component.

Figure 13

Flux linkage and EMF comparison, I _m = 10 A. (a) Flux linkage and (b) Induced voltage.

Figure 14

Torque and iron loss versus input current. (a) Torque versus current amplitude. (b) Iron loss versus current amplitude.

As the FEA model is used as a reference, it is indicated that the RN model achieves a very good agreement with fewer nodes and faster computation speed. The errors of the flux linkages and the voltage waveforms could be obtained by comparing their RMS values obtained from both approaches, and the maximum error is lower than a few percent. For the torque and the core loss, the error is obtained by directly comparing the values of these quantities, and again the error is lower than a few percent (Table 2).