Multi-material topology optimization of a ﬂ ux switching machine

. This paper investigates the topology optimization of the rotor of a 3-phase ﬂ ux-switching machine with 12 permanent magnets located within the stator. The objective is to ﬁ nd the steel distribution within the rotor that maximizes the average torque for a given stator, permanent magnets, and electrical currents. The optimization algorithm relies on a density method based on gradient descent. The adjoint variable method is used to compute the sensitivities ef ﬁ ciently. Since the rotor topology depends on the current feedings, this approach is tested on several electrical periods and returns alternative topologies. Then, the method is extended to the multi-material case and applied to optimize the non-magnet part of the stator. When dealing with 3 phases, the algorithm returns the reference topology as well as a theoretical machine with no return conductor according to the set current angle. To illustrate the creativity of the method, the optimization is ﬁ nally performed with a single-phase and returns a new topology.


Introduction
Topology optimization is a conception tool that makes no use of geometric parametrization.It was first developed in mechanical engineering by [1], then introduced in electrical engineering by [2], and has gained interest from engineers in the last decade with the development of additive manufacturing [3].Among various techniques, such as the level-set method [4], or the phase-field approach [5] (see [6] for a recent overview), density-based approaches [1,7] are the most popular.The principle of this technique is explained hereafter.
The geometry to be optimized is represented by a so-called density field q, which can be interpreted as pixelation on N e elements once discretized spatially, as shown in Figure 1.A density value of 0 represents air, and 1 represents iron in the corresponding mesh.The topology optimization problem then reads as find q opt ¼ arg min f q ð Þ subject to q 2 ½0; where f is the objective function to minimize (such as the cost, the mass, or in the present case the opposite of the average torque), and q is the vector of optimization variables, which contains the density of each mesh element.
To avoid solving a combinatorial problem that may be intractable, density methods introduce intermediate materials (also called gray materials) associated with density values strictly between 0 and 1.Therefore, the material properties can be continuously interpolated between the actual candidate materials to use a fast gradient-based optimization algorithm.
Although the presence of gray materials is not critical during optimization, they should not appear in the final geometry.Indeed, gray materials do not necessarily have a proper physical interpretation or may represent, in the best case, a microstructure that is complex to manufacture [8].A common solution for removing gray materials proposed by the literature involves penalizing the material interpolation, such as the Solid Isotropic Materials with Penalization (SIMP) scheme [8,9].Yet topology optimization problems are generally ill-posed, and other numerical artifacts can occur, such as checkerboards that need regularization (e.g., numerical filtering [10]).Density-based approaches have been applied to many magnetics problems, from the simple C-core electromagnets [11,12] to 3-phase synchronous reluctant machines with several materials and physics [13,14].However, apart from [15], no such work has studied unconventional structures such as Flux Switching Machines (FSM), where the unbiased creativity of topology optimization may be the most useful.
This article aims to apply a density-based topology optimization methodology to FSM.First proposed by [16], this type of machine contains a passive rotor with the field inductor located in the stator only, making its cooling easier [17] and suitable for high-speed applications.Several types of FSM exist in the literature; see [18].The chosen test case is a 3-phase Permanent Magnet Flux Switching (PMFS) machine with 12 permanent magnets evenly spaced within the stator, adapted from [19] and shown in Figure 2.
This paper is structured as follows.First, Section 2 recalls the magnetostatics equations.The optimization algorithm is then detailed in Section 3. Next, two numerical applications exploring unconventional topologies are presented and discussed.On the one hand, classical rotor optimizations are conducted under various electrical frequencies in Section 4. On the other hand, stator optimization, which requires a multi-material framework, explores the topologies obtained with 1 and 3-phase current feeding in Section 5.In conclusion, Section 6 summarizes the essential results and draws perspectives on this work.

Physical problem
Since density-based topology optimization will be used, we first set up the interpolations of material properties in Section 2.1 in the Kennelly convention [20].These quantities were then injected in the physical equation of magnetostatics recalled in Section 2.2, which is further discretized with the 2D Finite Element Method (FEM).

Material interpolation
There are at least two material properties of interest when dealing with a magnetostatics problem: the magnetic polarization of materials m1 and the current density.In the case of soft magnetic material, m is collinear to the flux density b and 0 when b = 0, so that it can be expressed in terms of a scalar reluctivity m(b) Two possibilities are generally reported in the literature to interpolate the magnetic properties when dealing with soft magnetic materials.One can interpolate either the magnetic reluctivity m as in [12,21,22], or the magnetic permeability l = m À1 as in [2,23,24], among many others.However, an additional magnetic quantity to interpolate (the remanence) should be introduced when dealing with polarized materials, as in [24,25].We choose an alternative, proposed in [26], which interpolates the magnetic polarization m of a FeCo alloy (AFK1 from Aperam [27]) directly.The corresponding BH curve is plotted in Figure 3.Its interpolation reads in the general case of n m materials where P m is a penalization function associated with the magnetic polarization, x i a shape function defined in [28] associated with the i-vertex of a convex interpolation domain D, and m i the magnetic polarization of the i-material.To address the topology optimization of the stator, the current density should also be interpolated with a similar scheme: where P j is a penalization function associated with the current density, and j i , is the out-of-plane component of the current density of the material i.The chosen domain D as well as the penalization functions are given in Sections 4 and 5, and an example of interpolation is drawn in Figure 4.More details on the interpolation scheme and associated algorithms can be found in [29].

Physical equations
From the static Maxwell's equations, one can write the magnetostatics equation for a 2D problem [26]: with mx and my the x and y components of the interpolated magnetic polarization m, jz the current density, and a z the z-component of the magnetic vector potential a, related to the flux density by the formula b = r Â a.
To be solved numerically, the strong equation ( 5) should be written in a weak form and discretized using, for instance, the FEM.Then, the physical problem reads as follows where K is the finite element matrix, a the vector containing the discretized degrees of freedom a z , and s right-hand side related to source terms.Figure 5 draws the mesh used in this article.Homogeneous Dirichlet boundary conditions confined the flux inside the machine, and a locked-step sliding band technique [30] is applied to emulate the rotation of the rotor without remeshing.The nonlinear FEM system ( 6) is solved with the Newton-Raphson scheme: with r n = Ka n À s(q, a n ) is the residual of ( 6) at iteration n.After convergence, the torque can be computed from the flux density field within the airgap e by Arkkio's method [31]: where b r , b h are the cylindrical components of the flux density, and R s and R r are the stator's inner radius and the rotor's outer radius, respectively.The axial length of the machine is normalized to L = 1 m.

Optimization algorithm
The purpose of the optimization is to solve the problem (1), i.e. to find an optimal density vector q that minimizes a given objective function f.In this article, the objective function is set to the opposite of the average torque (i.e., the average torque should be maximized), computed with (8) over 720 angular positions along one mechanical turn: As the density vector q contains many components (at least one per mesh element), the optimization algorithm should rely on the sensitivity d q f.The optimization flowchart is given in Figure 6.The efficient computation of d q f components is nontrivial since f rarely depends on q explicitly.Actually, f rather depends on the system's physical state, such as the magnetic field b in the present case.The latter depends   implicitly on the density vector q through the physical equation ( 6), which involves the material interpolations (3) and (4).

Sensitivities computation
The Adjoint Variable Method (AVM) is the most efficient approach to compute first-order derivatives from many variables.Once the FEM system ( 6) is solved, the sensitivities are computed using the adjoint state k, which is the solution of the following adjoint system dr da then the sensitivities with respect to each design variable q i can be calculated with Note that ( 10) is linear and should be solved only once per iteration.Consequently, the average single-thread computing time per optimization iteration of the sensitivities (41s) is almost negligible compared to the one of the nonlinear FEM (352s).This high computing time is caused by the high number of nodes and angular positions for the calculation.

Update
To accelerate the optimization to the real materials, each component d T of the descent direction associated with the mesh element T is normalized as follows: with d q T f the sensitivity vector of f to the part of the global density vector associated with the element T .A simplified trust-region algorithm [32] adapts the step size heuristically according to a quality indicator, which compares the measured variation of the objective function Df between two consecutive iterations with its variation predicted by the linearized model.It is defined as where Dq is the variation of the density vector between two consecutive iterations.If q is too low, the iteration is rejected, and the step size is reduced.If q is big enough, the iteration is accepted, and the step size can also be increased.The flowchart in Figure 7 gives more details on this algorithm that is applied in the following sections.
Then, q n+1 is projected orthogonally to remain within the bounds of D. The algorithm stops if the step a becomes smaller than 5 Â 10 À4 , or after 100 iterations.

Rotor optimization
The presented topology optimization algorithm is first applied to FSM rotor.The optimization domain is shown in Figure 8.In order to explore new structures, different electrical frequencies are injected within the FSM.The optimization setup is presented in Section 4.1, and the results in Section 4.2.

Parametrization
There are only two candidate materials to distribute inside the rotor: soft magnetic steel (FeCo alloy AFK1 from Aperam [27]) and air.Therefore, the interpolation domain D is a 1D line segment that joins both materials and supports the optimization variables contained in q; air corresponding to q = 0 and steel to q = 1.The density interpolates the magnetic polarization linearly:  Because there is no current neither in air nor steel, the current density inside the rotor jz is set to 0. However, the current density inside a stator conductor i, following the enumeration in Table 1, reads: where J = 10 A/mm2 is the current density amplitude, N is the number of electrical periods during one complete mechanical rotation, h m is the mechanical angle of the rotor, and w the electric angle set to 0°for the rotor optimization.The N values tested between 1 and 20, corresponding to different electrical frequencies.

Results
The final average torques of optimized machines are plotted in Figure 10.The highest torque is obtained for N = 10, corresponding to the current feeding of the reference machine shown in Figure 2. Therefore, a very similar rotor design shown in Figure 11a is obtained and develops about the same torque (hT 10 i = 2244 Nm/m) as the reference one (hT ref i = 2236 Nm/m) under the same current feeding.The results indicated in Figure 10 show that N = 4 and N = 16 can produce a lower torque than the reference but more significant than the rest of the N values.The associated designs are drawn in Figures 11c and 11e, and their convergence curves are compared to the one obtained with an unsuitable value N = 13 in Figures 12a and 12b.Although they reach a lower torque than the N = 10 structure, their topologies fundamentally differ from the reference rotor and are unusual for a human designer.Indeed, they contain disconnected parts, illustrating the ability of topology optimization to find innovative designs.
The algorithm fails to find meaningful structures for the rest of the N values, which may indicate the impossibility of designing suitable rotors under these conditions.An example of such a "bad" optimized design for N = 13 is drawn in Figure 13, which leads to the conjecture that there exist suitable rotors only for N = 4 + 6n, n 2 N.

Stator optimization
Density-based methods can be extended to multi-material problems, such as optimizing an FSM stator.Following the results of the previous section, we impose N = 10 with the reference rotor geometry shown in Figure 2. The magnets' positions and orientations are also fixed, given the optimization zone drawn in Figure 14.

Three-phase machine
The optimization of a 3-phase stator is first addressed.Since the permanent magnets are imposed, there are 8 candidate materials (steel, air, and conductors3 ) listed in Table 1.A multi-material interpolation formalism is necessary to address this problem.To handle this large amount of materials, the interpolations (3) and ( 4) are supported by a multidimensional space D. We chose D as a convex polytope that supports shape functions set {x i } i2s0,7t [28], represented in Figure 15.The diamond-shaped domain allows the distribution of the conductors onto the same plane with respect to their electrical phase.The basis function values are computed with an optimized version of the code given in [33], available at [34].
The chosen penalization functions P m and P j are adapted from [35] and read: The optimization is then performed for several current angles w 2 [0°, 60°].The resulting designs are given in Figure 16 with their average torques.The current angle strongly influences the resulting topologies; when it is incorrectly set, it may return non-manufacturable designs with no return conductors.Intermediate w values promote single conductors in each slot, as there is a factor ffiffi ffi 3 p

2
' 0:87 between the total current amplitude in sliced conductors and plain conductor slots, as illustrated in Figure 17.
The structure obtained with w = 0°is very similar to the reference shown in Figure 2.However, the conductors are larger, carrying more current for the same current density, set to J = 10 A/mm 2 .Consequently, the average torque of the optimized design is higher.Note that the order of conductors is inverted within the same slot.As topology optimization considers only electromagnetism performance, it places the conductors according to their electrical phase, resulting in distributed windings.Figure 2 proposes concentrated windings, probably for manufacturability reasons that are still ongoing research topics in topology optimization.

Single-phase machine
To illustrate the creativity of topology optimization, we now address an academic problem with no reference solution.More precisely, we want to optimize the stator of a 10/12 machine that has the same 10-tooth rotor and 12-magnet configuration, using only a single electrical phase, i.e., a positive and a negative conductor denoted as and .Single-phase FSMs are of interest for domestic and low-cost applications, and some classic topologies were compared in [36].
Since there are now only 4 candidate materials, the interpolation domain D is a tetrahedron drawn in Figure 18a, and the basis functions are the barycentric coordinates as illustrated in Figure 18b.The penalization functions are the same as in Section 5.1.
A parametric sweep on w was performed.The results are shown in Figure 19, and the corresponding torques are plotted in Figure 20.The case of w = 0°[60°] leads to manufacturable designs that reach 2200 Nm/m with a mixture of C-and E-core patterns [37,38].

Torque comparison
The dependence of the average torque with the current angle w is given in Figure 20 hereafter.
The rotor is almost insensitive to w because its optimized topology compensates for the electric phase shift by a mechanical rotation.This mechanical rotation is impossible for the stator optimization since the permanent magnet positions are imposed.For the threephase machine, the factor between the maximum and minimum average torque is 0:86 ' ffiffi 3 p 2 because of the sliced conductors, as predicted in Section 5.1.Since there are no sliced conductors within the stator slots for the single-phase machine structures, its average torque remains almost constant.
Even if it was not considered in the optimization, the torque ripple r ¼ maxðT ÞÀminðT Þ hT i is an important criterion when designing an electrical machine, and is given in Table 2.The corresponding instant torques of manufacturable machines (i.e., the structures obtained with w = 0°) are plotted in Figure 21.
The single-phase machine produces a pulsating torque with a high ripple but can still produce a significant average   torque.The two other optimized designs outperform the reference, and the best performance is obtained by optimizing the 3-phase stator.
In the end, the best-found topology is similar to the reference, but our topology optimization approach has explored other unusual structures.

Conclusion
This work presents and details a density-based topology optimization method applied to the torque maximization of a permanent magnet flux switching machine.Unconventional problems are addressed, such as the rotor optimization for several electrical frequencies and the stator optimization for different numbers of phases.Unusual structures were obtained, such as disconnected rotors or a hybrid single-phase stator made of C-and E-cores.This variety of results from an empty optimization domain highlights the creativity and unbiasedness of a topology optimization approach.The best-obtained topology is similar to the reference design but with better torque performances.In other applications, the optimal topology might differ from the known references.
In conclusion, topology optimization is a valuable tool to help design electrical machines that can be improved further.Future research will add the permanent magnet to the set of candidate materials so that the stator can be optimized in its entirety.Such work will be the first step to the topology optimization of a complete flux switching machine with both rotor and stator simultaneously.Additional objective functions can be added to tackle industrial problems such as ripple and back-electromotive force.Finally, other physics such as thermics and mechanics may be considered in order to obtain more ready-to-use designs.

Fig. 4 .
Fig. 4. Interpolation example of the current density amplitude between 7 materials.Note that in the real optimization problem in Section 5, the actual value of the current density at each electric angle is interpolated instead of its amplitude.

Fig. 3 .
Fig. 3. Anhysteretic BH curve of FeCo in logarithmic scale used for the computations.

Figure 12 .
Figure 12.Convergence curves during the optimization.a) Torque evolution, b) Step size evolution.

Figure 11 .
Figure 11.Resulting designs and flux maps with the highest torques for rotor optimization [15].

Figure 14 .
Figure 14.Design domain (gray) in case of the stator optimization.

Figure 15 .
Figure 15.Interpolation domain D that supports a set of shape functions in the 3-phase stator case.a) Material placement, b) Values of a shape function associated with B + .

Figure 17 .
Figure 17.Fresnel diagram of the current in sliced (black) and plain conductors (red).

Figure 18 .
Figure 18.Interpolation domain D supports a set of shape functions in the single-phase stator case.a) Materials placements, b) Values of a shape function associated with the FeCo steel.

19 .
Examples of single-phase optimized designs.

Figure 20 .
Figure 20.Average torque obtained for each w.

Table 2 .Figure 21 .
Figure 21.Instant torque obtained for the optimized designs for each w.

Table 1 .
Materials names and indices.