Development of Mathematical Models in Explicit Form for Design and Analysis of Axial Flux Permanent Magnet Synchronous Machines

: This article proposes a methodology for the design of double-sided coreless axial flux permanent magnet synchronous machines, which is based on a developed model for calculating the axial component of the magnetic flux density in the middle of the distance between opposite permanent magnets, which also represents the middle of the stator. Values for different geometric parameters represent the input data for the mathematical model in explicit form. The input data are calculated by using a simplified finite element method (FEM), which means that calculations of simplified 3D models are performed. The simplified model consists of two rotor disks with surface mounted permanent magnets and air between them, instead of stator windings. Such a simplification is possible due to similar values of permeability of the air and copper. For each simplified model of the machine the axial component of the magnetic flux density is analyzed along a line passing through the center of opposite permanent magnets and both rotor disks. Values at the middle of the distance between the opposite permanent magnets are the lowest and are therefore selected for the input data at different stator, rotor disks and permanent magnets (PM) thicknesses. Such input data enable the model to consider the nonlinearity of materials.


Introduction
Stators of axial flux permanent magnet synchronous machines (AFPMSMs) can be constructed with or without ferromagnetic cores [1].
Coreless topologies are used for low-and medium-power generators [2]. This article deals with a coreless double-sided AFPMSM with surface-mounted permanent magnets (PMs) on the massive ferromagnetic material of two external rotor disks.
The advantages of coreless topologies are the absence of a cogging torque [3] and the use of a massive instead of a laminated iron core in rotor discs. Due to the absence of losses in the ferromagnetic core of the stator, this topology can operate at a higher efficiency [4] compared to the conventional machines [3,5,6].
Torque and induced voltage sizes are mainly limited by: 1. outer dimensions of the machine and its mass and 2. current density in the windings.
The maximum permissible external dimensions of the machine limit the space for the installation of coils and permanent magnets, and the maximum allowed temperature limits the current density in the windings [5,7], whereby the torque and induced voltage depend on the active volume of the copper of the winding, the volume of permanent magnets and of the current density in the windings [8].
The volume of PMs determines the magnetic flux density in the air gap, and is determined by their thickness, by the number of poles of the machine, and by the external and internal radii of the active winding volume of the machine [9].
Magnetic flux density in the air gap also depends on the thickness of the stator, as it determines the distance between the opposite PMs. This distance together with the volume (dimensions) of PMs also influence the size of the attraction forces between the rotor disks [10], with which the thickness of the rotor disks is determined [2].
Induced voltage also depends on the magnetic flux density in the air gap, more precisely from the axial component of the magnetic flux density which passes between the opposite permanent magnets [11]. This dependence is written with the help of Faraday's induction law based on a variable magnetic flux which perpendicularly passes through the coil [3].
From the above we can conclude that the magnetic flux density in the air gap is one of the key factors [12] to produce torque and the induced voltage [13], since it is directly related to all the dimensions and materials in the machine. Therefore, the calculation of the magnetic flux density is a priority task in the design of the machine.
The distribution of the magnetic field in the AFPMSM can be determined using analytical and numerical methods among which the finite element method (FEM) is most commonly used. The advantage of calculations using the FEM is the higher accuracy of the results, but the weakness is in the longer required time for individual calculations [9].
Analytical methods are not time-consuming, but they have a disadvantage in that it is difficult to accurately capture all three dimensions of the distribution of the magnetic field in the machine [14].
There are many different approaches for the analysis of magnetic fields in the scientific literature. Here, only methods that are based on the middle radius of the machine, (i.e., PMs) are presented, since the methodology presented in this article is also based on the average radius of the machine.
The authors of [15] present the analytical model of an ironless one-sided AFPMSM, based on the distribution of the magnetic field in the medium radius. To determine the magnetic scalar potential, Laplace equations are solved in a Cartesian coordinate system using the Fourier-type method. To confirm the method and the accuracy of the results 2D analysis with FEM is also performed, which indicates that the deviation of the analytically obtained results is within 3% compared to the results obtained with FEM.
In [16], the analytical model of a coreless AFPMSM is presented for the calculation of the magnetic field in the air gap, which is based on the distribution of the magnetic field in the average radius of the machine (PMs). The authors note that the effects of curvature and the marginal (edge) effects of the machine on the magnetic flux must be considered. Azzouzi and others in [17] represent a unique solution for the calculation of a magnetic field with an average radius of the slotted AFPMSM under no load conditions. The effects of curvature and marginal (edge) effects are considered by a function, which is derived from the 3D FEM analysis.
Due to the shortcomings of the models based on the distribution of the magnetic field at the middle radius, the authors use "multilayer" models. This means that 2D models are calculated on different radii, thus considering the effects of curvature. To calculate quantities such as machine torque, the accumulated results of all layers [18][19][20] are considered.
In order to determine the distribution of the magnetic flux density in a coreless AFPMSM, a "multilayer" analytical model is proposed in [17,18], where the authors also considered the radial dependence of the excitation (the constant width of the winding side). In [21], an analytical model for the calculation of the magnetic field in the air gap of a slotted AFPMSM is presented. The particularity of the method is that the model is assembled for a slotless machine and these are later considered through the magnetic conductivity function. Iron losses were separately calculated by means of a nonlinear reluctance circuit. This article is the continuation of research presented in [22,23] and proposes a new methodology for design of double-sided coreless AFPMSM, which is based on a developed model for calculating the axial component of the magnetic flux density in the middle of the stator with two air gaps (in the middle of the distance between opposite permanent magnets). The input data for the mathematical model are summarized by an approximation method and include the effects of stator thicknesses, rotor thicknesses, PM thicknesses as well as the nonlinearity of the rotor disks.

Materials and Methods
AFPMSMs can be designed using analytical or numerical methods (FEM) [19]. FEM is used for solving complex problems, where the method usually uses magnetic potentials and divides two-dimensional (2D) or three-dimensional (3D) models of the machine to a limited number of elements.
The initial model of the machine can be determined, with a known magnetic circle and current in the windings (dimensions of PMs, rotor disks, air gap and stator), using a set of general equations, also known as "Sizing equations" [9,[24][25][26][27].
This model is then used as a starting point for more accurate calculations of individual parameters using FEM [28,29].
Using FEM for designing AFPMSMs is oftentimes consuming both in the construction of the model itself and in the computation time [3,[30][31][32][33]. So far, published scientific papers deal with the solution of this problem through the development of analytical methods [18,27,[29][30][31][32][33][34][35][36], whose derivations are often mathematically demanding and require certain simplifications and assumptions [37][38][39], which on one side influences the accuracy of the analytical model and on the other side enables a faster calculation by means of derived equations for the magnetic flux density in an explicit form [40].
The presented methodology uses simplified FEM calculations for the initial dataset. A simplified 3D model was constructed in Ansys Maxwell software, where instead of the stator only air was included. This simplification was chosen due to the similarity of the permeabilities of air and copper.
Simplified calculations with FEM were performed at various stator thicknesses with two air gaps, with particular interest in the values of the axial component of the magnetic flux density in the middle of the stator, since at this point the magnetic flux density is the lowest, which gives us reliable starting points for further and more-accurate calculations.
In the first step of the design, the types of PMs and their shapes and dimensions are usually chosen, as they have the greatest influence on the magnetic flux density in the air gap. The density of the magnetic flow in the air gap is mainly influenced by the thickness of the stator, the width of the two air slots, and thus the distance between the TM between the two disks and the thickness of the rotor discs.

Simplified FEM Calculations
The initial FEM analysis of the machine is represented by simplified FEM calculations aimed at determining the size of the axial component of the magnetic flux density along a line that extends between both rotor discs through the center of gravity of PMs. The simplification is that the calculations were performed without considering the windings; only the air between the opposite rotor disks with surface-mounted PM was considered, since the permeability of air and copper is almost the same. Figure 1a shows the initial AFPMSM model consisting of two rotor disks, 20 permanent magnets and a 100-mm virtual air gap (distance between PMs on opposite rotor disks), which consists of 20 layers. For each layer of the virtual air gap, a simplified calculation with FEM was performed, analyzing the axial component of the magnetic flux density across the virtual line passing through both the rotor disks, the virtual air gap, and the PMs on opposite rotor disks (through the center of gravity). The position of the line is shown in Figure 1b.    Simplified FEM calculations were performed with limitations, such as the condition that stator (virtual air gap) thickness is not larger the 4 times the PM thickness (ds ≤ 4hTM). Initial simplified FEM calculations were performed for the 3-phase AFPMSM with the data shown in Table 1 [41] and again repeated for the AFPMSM with the data shown in Table 2.

Least Square Approximation Method
Values for Bz_min were used to produce a polynomial using a least square approximation (LSA) method. This is a mathematical procedure that can find a curve that fits best to a known set of given points by minimizing the sum of the squares of the offsets ("the residuals") of the points from the curve [42]. The sum of the squares of the offsets was used instead of the offset absolute values because this allows the residuals to be treated as a continuous differentiable quantity [43]. Vertical least squares fitting proceeds by finding the sum of the squares of the vertical deviations R 2 of a set of n data points [43], which is presented in Equation (1).

( )
where R is the residual, yi is the data point calculated using simplified FEM, f is the fitting function, xi is an independent variable of the fitting function, and a1, a2, and an are coefficients of the fitting function. In our case the polynomial was chosen as a fitting function. Different order functions were used in order to determine the best fitting polynomial. Since three parameters were considered (stator thickness, rotor thickness and PM thickness) a function with 3 unknowns must be included (Figure 3).
Using the procedure described above, a mathematical model in explicit form was obtained for the calculation of the axial component of magnetic flux density in the middle of the stator, where PM thickness (from 5 to 10 mm), stator (up to 4 times PM thickness) and rotor disk thickness (from 4 to 15 mm) was considered (Equation (2) where Bz is the axial component of magnetic flux density in the middle of the stator, d is the thickness of the stator with two air gaps, hPM is the thickness of the PMs and dFe is the thickness of the rotor disks.
As is was mentioned above, the mathematical model is derived from the input data obtained by the simplified FEM, which was conducted with limitations, such as the condition that stator (virtual air gap) thickness is not larger than 4 times the PM thickness (ds ≤ 4 hTM) and with the PM span selected in such way that it covers 120° electrical. These limitations are selected on the basis of FEM analysis results, which are shown in Section 3.

Design and Analysis of AFPMSM
Designing AFPMSMs, as well as other typologies of drives, has steps. The first step is usually determined by requirements or constraints such as required torque size, speed, maximum permissible dimensions of the machine, cooling method and maximum induced voltage. Therefore, at the beginning of this process, the required dimensions of the machine, especially the internal diameter of the installation of permanent magnets and windings, the external diameter of the machine (rotor and stator) and the axial length of the machine must be specified as accurately as possible.
To determine these parameters, the standard approach is to use the so-called sizing equations [7,26]. Two types of sizing equations can be found in the scientific literature [44].
The first type is a classical equation (Equation (3) The emphasis of this article is on the second type of sizing equation, which refers to the connection between the electromagnetic torque and the basic geometric, electrical and magnetic parameters (Equation (4)) [45].

( )
In Equation (4), Tem is electromagnetic torque, Bz the axial component of magnetic flux density, Ai the average electric load on the radius that matches the inner radius of PMs, Ro the outer PM radius, Ri the inner PM radius and λ the ratio between the inner and outer PM radius.
This equation was first introduced in [45], and since then it has been widely used [9,21,46], since it is generally useful. From this equation we can detect the magnetic shear stress when divided by the active surface of the rotor. In the presented form, it achieves greater accuracy in the case of a trapezoidal distribution of the magnetic flux density and at higher values of electrical current through areas where the magnetic flux density is maximal. It is also possible to determine the torque density by dividing the equation with the mass of active machine parts (Nm/kg) or the total machine volume (Nm/m 3 ).
After completing the evaluation of the main geometric parameters, the design process can be continued according to the standard steps [9,44], which are presented below.
The number of the machine's synchronous revolutions according to Equation (5), can be determined: where ns is the synchronous speed of the machine, f the network frequency and p the number of pole pairs. It can be seen that the higher number of poles of the machine is required for lower synchronous rotations.
The winding factor (kw) is determined by the geometric method, from the ratio between the geometric (egeometric) and arithmetic (earithmetic) value of the induced voltage (Equation (6)) [23]. where lFe is the magnetic flux path in the iron, dag the air gap thickness (between rotor and stator), dFe the rotor disk thickness and μr the relative permeability of steel.
Axial component of magnetic flux density (Bz) is determined by Equation (8 (8) where Br is the remanent magnetic flux density, dag the air gap thickness (between rotor and stator), ds the stator thickness, μrec the recoil permeability, hPM the PM thickness and ksat the saturation factor.
Ratio between PM pitch and pole pitch (αi) is independent of the radius and can be written as Equation (9) where τp is the pole pitch and τPM the PM pitch. Magnetic flux (Φ), induced voltage and electromagnetic torque can be determined using Equations (10)-(12), respectively [9]: where Bz is the axial component of magnetic flux density, p the number of poles, Do and Di the outer and inner PM diameter, E the induced voltage, f the network frequency, N1 the number of turns per coil, kw the winding factor, Φ the magnetic flux, Tem the electromagnetic torque, αi the PM pitch-topole pitch ratio, m the number of phases, and Ia the effective value of stator current.

Influence of Different Geometrical Parameters on Characteristics of the Machine
In this section, firstly the influence of PM thickness and rotor disk thickness on the axial component of magnetic flux density is calculated by using simplified FEM. Additionally, the influence of PM span and stator thickness on torque and induced voltage of the machine was analyzed by using the 3D FEM. Figure 4 shows the influence of PM thickness on the Bz size in the middle of the stator for the machine described in Table 1. This analysis was carried out using simplified FEM. It can be seen in Figure 4 that the PM thickness has a significant influence on the size of the axial component of magnetic flux density and on its degree of changing with the increase in the distance between PMs on opposite rotor disks. Figure 5 shows the influence of rotor disk thickness on the Bz size in the middle of the stator for the machine described in Table 1. This analysis was also carried out using the simplified FEM.  Figure 5 shows that rotor disks thickness also influences the size of the axial component of magnetic flux density up to a specific thickness. This is the reason that the rotor disk thickness has to be carefully selected. When designing the machine, namely rotor disks, the thickness has to be carefully selected in order to withstand the pull forces from the PMs (to keep the air gap constant and prevent the deflection of rotor disks) on one hand and on the other hand it influences the total weight of the machine and its efficiency.
Since some characteristics of the machine, namely waveforms, cannot be analyzed by the simplified FEM, a regular 3D FEM analysis was performed where PM span and stator thickness influence were analyzed. Figures 6 and 7 show the influence of PM span on the torque and induced voltage for the machine described in Table 1 and Figure 8 shows the influence of stator thickness on the torque size of the machine.     Figure 6 shows that PM span has a positive influence of the torque size of the machine, since the torque value increases with the increase in PM span. As we can see in Figure 7, the PM span increase has a different influence on the induced voltage of the machine. Its value does not change significantly, but its waveform does. We can see that for higher PM spans the induced voltage waveform has a flat top and for lower PM spans it has a more triangle-shaped waveform. It can be noticed that for 24° PM span the induced voltage has a sinusoidal waveform. That is the reason for all the input data for the derivation of the mathematical model being obtained with the simplified FEM calculations of a model with a 24° PM span (it covers 120° electrical).
When we derived the mathematical model, we also stated that the ratio between stator and PM thickness can be up to 1:4. Figure 8 shows the reason for this limitation, since the torque value increase rate reduces with the increase in the stator thickness.

Design and Analysis of AFPMSM
Two coreless AFPMSM were designed (described in Tables 1 and 2) by three methods, namely an analytical method (Equations (6)-(13)), and a combination of an analytical method and mathematical models in explicit form and FEM (Ansys Maxwell 3D).
Results of the design by analytical method are shown in Table 3, where different stator, rotor and PM thicknesses were considered. Table 4 shows the results of combination of the analytical method and mathematical models in explicit form and Table 5 shows the results gained by FEM as well as the comparison of the results of all three methods.  Table 5 shows that all three used methods produce similar results with minor deviations at certain points. For the selected parameters of the AFPMSM, the matching of the results of the torque size of the analytical method and FEM is 102.98% and between the mathematical model in explicit form and FEM is 98.54%.
Matching of the induced voltage between the analytical method and FEM is 93.62% and between the mathematical model in explicit form and FEM is 97.91%.
The procedure of designing was repeated for the 20-pole AFPMSM, described in Table 2 with the same rotor disk thicknesses used for the AFPMSM in Table 1. A new mathematical model with new input data for different pole number was produced for this purpose by the same procedure as presented in the previous section. Table 6 shows the comparison of the torque size, calculated by using three different methods for designing AFPMSMs. For the selected parameters, the matching of the results between the analytical method and FEM is 132.14% and between the mathematical model in explicit form and FEM is 95.43%. It can be seen that the mathematical model in explicit form produces more accurate results for PM thicknesses hPM = 5 mm and hPM = 7 mm in comparison with the analytical method. The reason for that can be found in the fact that the analytical method does not consider the PM span, but it works under the assumption that the PM surface and air gap surface are the same.
For hPM = 10 mm, neither the analytical method nor mathematical model in explicit form produce accurate results compared to the gained FEM results. The reason for this is in the magnetic flux leakage, since the PM thickness is larger than the distance between adjacent PMs on the same rotor disk (hPM > lPM). These dimensions are also shown in Figure 9.  In the case of hPM = 10 mm, the magnetic flux leakage between adjacent PMs on the same rotor disk increases (Path 3). This leakage is also increased by increasing the stator thickness, due to the magnetic flux path enclosing in Path 3 and the decreasing share of magnetic flux in Path 1.
In Table 6, it can also be seen that the developed mathematical model in explicit form produces accurate results for the stator thickness ds = 20 mm, which represents twice the PM thickness. For larger stator thicknesses magnetic flux leakage increases and the developed mathematical model in

Laboratory Measurements
Laboratory measurements were performed on the prototype AFPMSM with the data presented in Table 1. Measurements of the torque size for different loadings were performed as well as the induced voltage measurements at different rotation speeds. Table 7 shows the measured values of induced voltages in individual phases of the machine and average value for all three phases at different rotation speeds. Table 7. Induced voltages as different rotation speeds.  Table 8 shows the comparison of average measured values of induced voltages and calculated values using the analytical method and developed mathematical model in explicit form.  Table 9 shows the comparison of measured values of torque at different electrical current values and comparison with the calculated values, using the analytical method and developed mathematical model in explicit form.  Tables 8 and 9 show that the developed mathematical model is accurate, since the matching of the measured values of torque and induced voltage are in good agreement with the calculated results. Some deviation is present in the matching of the results which may be caused by the measuring equipment. For this purpose, we measured the axial component of magnetic flux density in the middle of the distance between PMs on opposite rotor disks.
To measure the axial component of magnetic flux density (Bz) in the middle of the distance between PMs on opposite rotor disks we developed a special tool that enables the PMs to stay in their position and measure the values at exact positions. Two tools were modeled and 3D printed (one 17mm thick and one 22-mm thick). Figure 10a shows the model of the tool with the location of the opening where the measuring probe (Magnetometer KOSHAWA 5 with the magnetic probe) will be inserted and Figure 10b shows the printed tool. The measuring probe is inserted into the opening, which is printed to the middle radius of the PMs.   Results of the measurements of the axial component of magnetic flux density (Bz) for different rotor disk thicknesses and two stator thicknesses are shown and compared to the calculated values (using FEM and the developed mathematical model) in Table 10.  Table 10 also shows the accuracy of the presented methodology, since the average matching between the measured values and calculated values is 96.92%.

Conclusions
The article presents a complete procedure for the development of the mathematical models in explicit form for the calculation of Bz in the middle of the stator of a double-sided coreless AFPMSM. The developed model simultaneously considers the influence of the thicknesses of rotor disks, PMs and stators along with the nonlinearity of the materials.
The FEM, of course, also considers the nonlinearity of the materials, but it is often a timeconsuming method, since in the changes of geometric parameters it is necessary to compile and calculate a new model of the machine. The developed mathematical model consists of input data, which are determined using simplified FEM calculations and already consider the different geometrical parameters of the machine (ds, hPM and dFe), so they do not require new long-lasting calculations.
In addition to the results of calculations with the FEM, we can see that the measurements also confirm the accuracy of the methodology.