Method for optimizing the magnetic circuit of a linear generator using FEM simulations

Within the area of permanent magnet electrical machines, there is an ongoing focus on replacing the rare earth permanent magnets with alternatives. An option is hard ferrites, commonly used in other applications. The relatively low coercive field strength of the ferrite magnets makes irreversible demagnetization an area that should not be neglected. In this paper, a methodology is proposed for the optimization of a slow-moving linear generator simulated in a finite element environment. The no-load phase voltage is maximized while accounting for iron saturation and permanent magnet irreversible demagnetization. This demagnetization is considered when the translator is alongside either the stator or air. The inclination angle between magnetization and magnetic field strength is accounted for by adjusting the intrinsic coercivity for each element of the permanent magnets. Characteristics for the magnet grades Y30 and Y40 are used in the optimization process. The velocity of the translator is set to resemble a speed common to wave power applications. Commercial finite element software is used together with two optimization algorithms: the genetic algorithm and the particle swarm optimization. The results of these optimization algorithms reach similar optimal solutions for the considered objective function, assuring a result close to a global maximum. The results also show a great difference in the optimal geometry for the two magnet grades and highlight the need to account for irreversible demagnetization when designing generators with ferrite magnets.


I. INTRODUCTION
There is public interest in renewable energy, and energy extraction from ocean waves is no exception. For a directly connected point absorber, a linear generator is a suitable solution for extracting energy from waves. Optimization methods for generators are more common among rotating machines because of their broad range in application and the historically more prominent establishment. When it comes to optimization of linear generators, there are, in comparison, few studies made. Among these are optimizations for minimized active power loss for tubular generators and total active volume for the following: rectangular generators; 1 maximized rectified output voltage with a reduced cogging force; 2,3 maximized power-to-weight ratio; 4,5 minimized cogging force; 6 minimized detent force; 7,8 reaching natural frequency of the moving mechanism by adjusted generator damping force; 9 maximized efficiency and minimized total loss; 10 efficiency in regard to cogging force and output voltage; 11 maximized electrical power output with a minimized steel core volume; 12 magnetic flux density distribution in regard to saturation; 13 maximized induced voltage per unit cost; 14 reduced thrust ripple with increased thrust force; 15 and back electromotive force and total harmonic distortion. 16 The choice of optimization algorithm varies between the reviewed studies. Among these studies are the use of the following: the genetic algorithm; 7,12 scatter search; 1 particle swarm optimization; 1 a hybridization between the genetic algorithm and particle swarm optimization; 2,3 and gradient based local optimization algorithm r-algorithm. 4,5 In some studies, the optimal parameters are derived based on a parametric sweep over certain value intervals. 6,[8][9][10][11]13,16 Neither one of the linear generator optimizations above has explicitly regarded the optimization for ferrite permanent magnets (PMs) [with the exception of a study where ferrites (C11) were compared to an earlier optimized geometry for Nd-Fe-B 14 ]. The ARTICLE scitation.org/journal/adv majority consider only rare earth PMs. [1][2][3][4][5][6][7][9][10][11] Because of the unstable prices of rare earth PMs lately, the possibility of substituting the rare earth PMs with hard ferrites is given increasing interest. As a result of the transition to magnets with a much lower remanence, the generator design must be altered to keep the flux density in the air gap at a similar magnitude. This is done by means of flux concentration. In addition to the lower remanence, there is a transition to magnets with a lower coercive field strength. This causes the PMs to be more susceptible to irreversible demagnetization. Most conventional optimization methods focus on output power or torque/force for machines with rare earth magnets. When using weaker magnets, the increased risk of irreversible demagnetization should not be overlooked. Studies show that it is not only the antiparallel component of the field strength in regard to the magnetization that can give rise to an irreversible demagnetization: [17][18][19] there can still be irreversible demagnetization when this inclination is at right angles. The behavior of the inclination angle is primarily investigated for rare earth PMs 17,18 but similar deviations are also shown for, e.g., SrFe 11.6 0 19 ferrites, which have a remanence close to the saturation magnetization. 19 Design optimization for PMs prone to irreversible demagnetization should account for demagnetization risks including, if possible, this inclined angle dependency. An earlier study of rotating electrical machines with a span of PM magnetic properties uses a methodology slightly resembling the one presented here. 20 The rotor geometry is then optimized in regard to torque using the Matlab ® function fminbnd. This function finds a minimum solution to a continuous function for a single parameter. In that study, the irreversible demagnetization is only determined based on the magnetic flux density component parallel to the magnetization.
To the extent of the authors' knowledge, no study of generator design optimization algorithm has included the inclination angle dependency of the irreversible PM demagnetization, in the least not when using rare earth free PMs, such as ferrites and certainly not for a linear generator.
In this paper, a proposed methodology for optimization of a slow-moving linear generator is described and demonstrated; genetic algorithm and particle swarm optimization are used to find an optimal magnetic circuit to maximize the no-load voltage of the generator in finite element (FE) simulation software while accounting for demagnetization of the PMs and saturation of the translator and stator iron. Although the optimized solution in itself is desirable, the focus of this paper is the methodology of finding such optimum.

II. METHOD
The linear generator is simulated using Comsol Multiphysics ® .
Since the optimization algorithms are running from Matlab Optimization Toolbox TM , the Comsol model is built within Matlab using the LiveLink TM feature between the two programs. The method for this paper is, therefore, divided into three subsections: one containing general information for the linear generator, one that handles the environment of the FE model in Comsol, and one that handles the optimization methods implemented in Matlab. Since the Comsol model is built and modified through Matlab, there is some overlap between the later two subsections.

A. Linear generator
The linear generator used in this paper is of a rectangular shape. This is to coincide with the research done by the authors' research group at Uppsala University, 21,22 where a set of such rectangular blocks forms the generator as a whole, as illustrated in Fig. 1. The authors' research group has, among other things, built and evaluated ferrite linear generators. 22,23 Buried PMs are used with alternating polarity and pole shoes placed in between, as can be seen in Fig. 2. Comparing the linear generator with buried magnets to the more commonly used surface mounted magnets, the length in the direction of magnetization of the PM is restricted by the pole pitch for the buried magnet, while for the surface mounted magnet, the length perpendicular to the magnetization is restricted. For buried magnets, the length perpendicular to the magnetization is free to increase. The same freedom also applies to the pole shoes. Throughout the paper, height is used in the direction of PM magnetization and width in the direction perpendicular to the PM magnetization. In the model, the translator consists of PMs with iron in between. In reality, a translator interior is needed to fixate the pole shoes and PMs. This translator interior should be non-magnetic with similar permeability to air. The constant flux density of the interior in the moving reference frame results in negligible induced electric field and eddy currents. There is a small area of air moving with the translator and yet another larger area of stationary air to allow magnetic flux on both sides of the translator. This is visualized in Fig. 3.
The number of slots per pole and phase, q, that links the translator pole pitch to the stator slot pitch is assumed constant. In this paper, q is set to 6/5 and a winding pattern similar to earlier studies for wave energy converters is used. 24 Between the cables wound in the stator, there is a layer of air of dimensions similar to the insulation used earlier. 24 Fixed parameter values for the generator can be seen in Table I.  A change in the characteristics of an electrical machine may have a large impact on its performance. Certain changes, such as winding patterns, can appear both stochastic and non-gradient to the objective function. For every increase in the parametric degree of freedom, there is an even further need to increase the number of iterations to find an optimal solution. In consideration of this, the study only handles three parameters that are allowed to change independently. These independent parameters are the width of the PMs, PM width ; the ratio of PM height to pole pitch, PM % = PM height /τp; and PS % , implying the partial removal of the pole shoe at the far left side of the translator that is used to force the flux to the right. 25,26 The remaining iron of the pole shoe is given as a percentage PS % of PM width , where PS width = PS % ⋅PM width . A geometric example of the parameters being optimized can be seen in Fig. 2.
The PM grades used in this paper are accessed from online and can be seen in Table II. 27 The average values are used for every magnitude.

B. Setup of the FE model
The linear generator is set up in a 2D quasi-static magnetic field environment within Comsol Multiphysics, solving for the magnetic vector potential. Quadratic discretization is used for the magnetic vector potential. Because of the repetition of the winding pattern every 36 slots, a periodic condition is introduced with a repetition of the magnetic vector potential Aup(x) = A down (x). The number of PM poles required is, therefore, p = Q/mq = 10, where Q is the number of slots, m is the number of phases, and q is the earlier mentioned slots per pole and phase. The Dirichlet boundary condition with the magnetic vector potential A = 0 is set on the outer boundaries of the air on either side. The motion of the translator is allowed using a moving mesh. To allow for the flux to cross the boundary between the translator and the stator, one can introduce a perfect magnetic conductor-interface (n × H = 0). 28 Instead, a modification of the continuous boundary formulation is used (A dest. = Asource), where the magnetic vector potential is actively mapped within the continuity formulation to the appropriate position as in Fig. 3. In Fig. 4, the magnetic vector potential and the norm of the magnetic flux density are plotted. The same formulation is also used on the left moving boundary. An example of the mesh used in this paper can be seen in Fig. 3(b). A too coarse mesh will cause irregularities in the objective function. Such irregularities should be minimized to avoid unwanted exploitation in the optimization process. The air gap is split into several layers to further enhance the mesh.
The PMs are modeled in Comsol using the remanent flux density (Br) formulation as follows: A relative permeability greater than 1 allows for reversible changes in the magnetization of the PMs, given by the susceptibility χ = ∂M ∂H = μr − 1. The relative permeability is derived from the maximum energy product with the assumption that the B-H curve is linear above the point of maximum energy product in the second quadrant. This linear behavior of certain ferrites can be seen for Y30 from the manufacturer, 27 for Y33BH 29 and is expressed in the literature. 30 This yields a maximum of the energy product BHmax at B = B r 2 and H = B−B r μ 0 μ r = − B r 2μ 0 μ r . This, in turn, gives the relative permeability as a function of remanent flux density and BHmax stated as follows: This expression for the relative permeability in combination with ferrites 31 and rare earth 32 PMs has been used in previous research. Combining this expression with the material data in Table II does, however, give rather different magnetic flux density values for the knee-point due largely to the big difference in the intrinsic coercivity. This is visualized in Fig. 6. Initial simulations are computed where the permeability is defined as either a tensor (μr,x = μr,z ≠ μr,y) or a constant value. The difference in time consumption is roughly 12% more for simulations with tensor (μr,x = μr,z = 1) with little influence on the result. With no data on μr,x and μr,z and to save time, the permeability here is assumed constant. Silicon Steel NGO 50PN270 is used as the material for the stator and translator iron. This material is accessed from the Comsol material database. The B-H curve and the μ 0 M-H curve are shown in Fig. 5. From the figure, it can be seen that the iron reaches saturation at roughly 1.7 T. It can also be seen that there is a zero crossing at zero applied field in (b). In other words, there are no hysteresis losses of the iron in the model. Eddy currents are also neglected by setting a zero conductivity of the iron. With an electrical conductivity for ferrite magnets that is far less than that of iron, the conductivity for the PMs is also set to zero. A comparison is made between the simulated formulation in Comsol and an in-house software program earlier experimentally verified, 33 showing a difference in peak air gap magnetic flux density values of less than 2% for surface mounted Y40 magnets. Experimental verification of Comsol for a multi-pole ferrite generator has earlier been presented for a rotating machine. 34 The no-load voltage can either be computed directly from Comsol or by (3). The latter originates from Ampere's law ∇ × E = − ∂(∇×A) ∂t and Faraday's law for the electromotive force For a 2D simulation with the conductors in z, the average electric field over one conductorĒ cond is given bȳ Assuming homogeneous conductors with no skin effects, the induced no-load voltage for one phase is, thus, where l stack is the out-of-plane stack length of the generator and is here set to 1 m. The fraction S tot,phase /S cond corresponds to the number of conductors per phase. Every turn of the closed integral of

ARTICLE scitation.org/journal/adv
Faraday's law is comprised of two conductors of length l stack . A comparison is done between the induced phase voltage derived from integration over the conductors based on (3) and the built-in Comsol function, showing errors in the root mean square value of less than 0.625% between the two approaches. Because of the negligible difference between the two, the built-in Comsol function for the no-load voltage is used for the optimization algorithms. The temperature is kept at a constant temperature of 20 ○ C throughout the simulations.

C. Setup of the optimization in Matlab
In this study, there are two different optimization algorithms that are used in parallel: the genetic algorithm and the particle swarm optimization. These optimization algorithms are accessed through Matlab as functions. They are both population based algorithms and do not depend upon the derivative of the solution. The genetic algorithm is made to mimic the behavior of a population through generations and revolves around the properties of the children. The properties can either be mutated with random changes from one parent, a crossover between two parents or, if being the child with best score, directly passed onto the next generation without alteration. 35,36 The particle swarm optimization was initially developed to mimic the behavior of birds but transitioned more to resemble that of particles moving through space. The velocities of the particles are based on its own best location and the best location of the whole swarm. 37 One of the benefits of the genetic algorithm is that it gives the user the possibility of using integer values. This is advantageous when, e.g., varying the number of coils per slot or changing between several winding patterns. Using more than one optimization algorithm in parallel allows a comparison between optimal solutions. This study has the benefit of comparing the performance of the optimization algorithms as only non-integer parameters are used.
Comsol is operated through Matlab by the LiveLink feature. For every iteration, new parameter values are generated from the optimization algorithms and used in the Comsol model. The results are then extracted and used for the objective function. The objective function chosen for this paper is the median no-load peak voltage of the three phases. It is possible to couple the magnetic field interface to electric circuits. This will, however, increase the computational time severely. Considering the total time in Table IV makes this  impractical. A higher flux density in the silicon steel increases the no-load voltage. An increased field strength will still slightly increase the flux density even when the steel is saturated, with a slight increase in the no-load voltage. To account for iron saturation, a reduction factor kr,sat to the objective function is introduced. This reduction factor takes effect after a certain magnetic flux density is reached, and as the flux increases further, the reduction factor restricts the objective function. In this paper, kr,sat has a constant value of 1 until the flux density in the iron reaches 1.7 T. Between 1.7 T and 1.8 T, the reduction factor has a linear decline toward 0.
It is also necessary to account for the risk of irreversible PM demagnetization. The formulation of the magnetic flux density used in Comsol is stated in (1). This is illustrated by the red dotted line in Fig. 6 and assumes a constant susceptibility of the PM. At some point, however, there will be a non-linear relation between the magnetization and the magnetic field strength. There are several approaches to take this into account. One way is to approximate the B-H curve as two straight lines intercepted at a knee-point. 38 This is illustrated by the blue lines (dark gray in grayscale) in Fig. 6. The lower of the connected lines is defined as the line that crosses both B = 0 at the coercive field strength H = H c,b and B = μ 0 H c,j at the intrinsic coercivity H = H c,j when the magnetization is zero. Another approach is to model the B-H curve as an analytical exponential function. 39 This is given by the following equation and is illustrated by the yellow lines (light gray in grayscale) in Fig. 6: where K 1 defines the "sharpness" of the knee-point. The exponential function is used in this study to express the non-linear demagnetization of the PMs. No adequate demagnetization curves were found for Y30 and Y40 to validate the knee-point. Therefore, the sharpness is set where good agreement was shown for higher remanence magnets 40 with the value K 1 = −1.5 × 10 −4 m A −1 . The K 1 value can be modified to correspond to experimental data for every grade, but for this paper, it is assumed constant. K 2 is defined in (5) as a function of intrinsic coercivity H c,j , K 1 , and the conversion factor E = 1 T, Once in the nonlinear region of the magnetization curve, the B-H curve will recoil to a reduced Br. This requires the solver to actively update Br, which is something that is highly time-consuming and, in turn, unpractical for an optimization algorithm with several hundreds of iterations. Therefore, a second reduction factor is introduced in the post-processing and takes effect when the demagnetization of the PM exceeds the point where the exponential function ARTICLE scitation.org/journal/adv deviates from the linear function. This point is illustrated by the asterisk in Fig. 6 and is defined as the point when the exponential function differs by more than 1 mT to the straight line given by (1). Below this point, the PM will begin to recoil toward a lower Br. Comparing (1) and (4), it is easily seen that the difference between the two equations is given by This relationship between (1) and (4) is only valid when the relative permeability and Br are the same for both equations. Updating Br in the time dependent simulation will not change Br in (4) but will cause a recoil to a reduced Br in (1), making (6) valid along the major hysteresis loop where the parameters of (4) are first defined.
After that, the change in Br must be accounted for. The minor loop related to the reduced Br is often approximated by a linear recoil permeability of similar value to the relative permeability of the major loop.
The literature reveals an intrinsic coercivity dependency in regard to the angle between the magnetic field strength and magnetization, which is not compliant with a perfect demagnetization protection at right angles. 17,41,42 One such case with large angles between the magnetic field strength and magnetization is buried linear generators where the flux-concentrating topology forces the flux to change direction close to the border of the PM. The sample size of 4 mm diameter and 20 mm height was used in pulsed-field measurements in one study, 17 hence on a scale that is much smaller than that of the simulation of an electrical machine. To account for the angular dependency, one can estimate alignment distribution functions to regard the misalignment of the easy axis. 41 Another way is to use a polynomial function to directly compensate the angular dependence with an increased magnitude of the intrinsic coercivity. 17 The polynomial function in (7) is used to set the intrinsic coercivity for every node of the PMs, where γ is the angle between the magnetic field strength and the magnetization, H ang c,j = H c,j (1 + a 1 γ + a 2 γ 2 + a 3 γ 3 ).
Using this polynomial to account for the angular deviation of the intrinsic coercivity, the magnitude of the field strength is looked at rather than just the field strength in the antiparallel direction.
Comparing the inclination dependency on the intrinsic coercivity at 20 ○ C of the ferrites SrFe 11.6 0 19 19 to the rare earth magnets Nd 14.2 B 6.2 Co 1.0 Fe bal. and Nd 14.2 Dy 0.3 B 6.2 Co 1.0 Fe bal. 18 shows little difference when the remanence of each PM is close to the saturation magnetization. With no access to Y30 or Y40 data to account for the inclination, similar per unit behavior is assumed 17 and α 1 , α 2 , and α 3 are set accordingly (a 1 = 3.17 × 10 −4 deg −1 , a 2 = −3.38 × 10 −5 deg −2 , and a 3 = 1.37 × 10 −6 deg −3 ). The demagnetization curves for different H ang c,j can be seen in Fig. 7 together with the straight line given by (1) and the positions on the curve where the exponential function starts to deviate from the straight line.
The factor k r,demag accounting for demagnetization looks at how much of the PMs (given as a percentage) that exceeds the point given by the asterisk for every node in the FE simulation. The equation for each node is given as |Δeqs| ≤ 1 mT, 0 (above asterisk) |Δeqs| > 1 mT, 1 (beneath asterisk). The demagnetization quota, derived as 1 S PM ∫ S C node dS and given as a percentage, is then used in the reduction factor k r,demag . The reduction factor is expressed as the linear decline from 1 to 0 when the demagnetization quota is increased from 0% to 10%. The overestimation of the simulated no-load voltage caused by the lack of updating to a reduced Br in the corners is negligible compared to the reduced value of the objective function caused by the reduction factor. One thing to emphasize is that neither of the reduction factors influences the computation of the FE model. They only influence the objective function in the post-processing. The spatial distributions (e.g., H ang c,j /H c,j , Δeqs, and C node ) are defined in the Comsol environment as either functions (upper and lower limits) or variables (no limits).
Unlike a rotating machine, the translator is often assembled before it is inserted into the stator. 43 In addition, partial overlap of the translator might be allowed. 22 For the no-load case, the demagnetization is more severe when the translator is surrounded by air because of the forced increase in reluctance path. A previous study on demagnetization of ferrites in a spoke-type configuration during a short circuit concluded that the worst demagnetization case was when the magnet was placed in air during assembly. 44 In regard to the difference in reluctance path, there are two computations of the Comsol model for every iteration of the objective function: one that computes the no-load phase voltages, the saturation factor and demagnetization factor with a stator adjacent; and one that computes the demagnetization factor when the translator is surrounded only by air.
The objective function is given by (9), where τp is the pole pitch and N PM is the number of magnets, f obj = − V no-load,peak kr,satk r,demag τpN PM .

ARTICLE scitation.org/journal/adv
The product τpNPM is the total length for the simulated generator needed to enforce periodicity. With dimensionless reduction factors, the unit of the objective function f obj is given as V m −1 . The optimization algorithms will try to find the minimum of the objective function min(f obj ). This is the reason for the negative sign of (9). The two reduction factors will thus limit the maximum flux density (kr,sat) in the iron and the minimum flux density (k r,demag ) in the PMs.
The flowchart for the entire optimization process can be seen in Fig. 8. The genetic algorithm and particle swarm optimization are both metaheuristic and cannot guarantee a solution that is a global optimum. Once the search is complete (20 iterations through 20 generations), parametric sweeps are done in the proximity of the optimal solutions to see if the solution is at a global maximum.

III. RESULTS
The optimal designs for Y30 and Y40 for the given objective function are shown in Figs. 9 and 10 and presented together with their achieved objective functions in Table III. The related demagnetization quotas and maximum flux densities can be found in Table IV. There is a notable difference in optimal design for the two PM grades. The low magnitude of the intrinsic coercivity in the Y30 PM grade restricts the design of the PM to be high with short width. It also prohibits the PM to be directly adjacent to air in the direction of magnetization, dismissing any partial removal of the iron material in the pole shoes. Comparing the maximum flux density from Table IV to the B-H curve for iron in Fig. 5, it is evident that the optimal design for Y30 is very far from reaching any iron saturation. On the contrary, the iron is not well utilized as the maximum value for the flux density is less than half of what is allowed before the reduction factor for saturation starts affecting the objective function. The Y40 PM grade does, on the other hand, allow some partial removal of the iron to allow further flux concentration through the air gap. These PMs are also allowed to be wider, resulting in an increased total flux entering the stator. From the tables, it can be seen that the flux density for Y40 reaches the threshold of 1.7 T, where kr,sat begins to influence the objective function. There is also less difference in the severity of demagnetization quota when comparing the translator being surrounded by stator or air. In comparison to the Y30 PM, higher flux densities are reached in the translator iron rather than in the stator. Both Y30 and Y40 reach similar values for the maximum demagnetization quota for the optimal solutions. Parametric sweeps are done in the proximity of the optimal solution for each optimized parameter. These distributions are presented in Fig. 11. Comparing the two optimization algorithms, it can be seen that they reach similar values for the objective function, strongly suggesting a solution close to the global optimum.
Surface distributions for field strength H, inclination angle γ, normalized intrinsic coercivity H ang c,j /H c,j , and the demagnetization distribution C node for the optimal solutions are shown in Figs. 12-15. The scales in Figs. 14 and 15 are the same to illustrate the difference in magnitude. The inclination of H gives rise to an angle γ, which, in turn, increases the magnitude of the intrinsic coercivity H ang c,j when compared to the magnitude of the field strength H. From the figures, it can be seen that the inclination reaches values up to roughly 45 ○ . Equation (7) yields, at this angle, the normalized intrinsic coercivity H ang c,j /H c,j = 1.0707. This value, in combination with the magnitude of H, indicates an increased susceptibility to demagnetization compared to only considering the antiparallel value of H, with a value of 1 √ 2 |H| at 45 ○ . From Fig. 15(a) of Y40 surrounded by air, it is evident that the field strength H is more distributed in the left corners.
The low inclination angles close to the left boundary of the pole shoe result in a lower magnitude of the intrinsic coercivity H ang c,j . The relatively low magnitude of the intrinsic coercivity in combination with the more distributed field strength results in an area that is more susceptible to irreversible demagnetization.

ARTICLE scitation.org/journal/adv
The scores of the objective function for different generations can be seen in Fig. 16 for each optimization algorithm and PM grade, respectively. Particle swarm optimization reaches greater values in the early generations than the genetic algorithm. The improvement per generation seems, however, to stagnate in the later generations and converge to a value similar to the genetic algorithm. Other combinations of iterations and generations may show different behavior, and although the particle swarm optimization algorithm reaches adequate values earlier, the genetic algorithm has the clear advantage of working well with "integer only" parameters.

IV. CONCLUSION
An optimization method for the magnetic circuit of a linear generator has been presented and demonstrated. In Sec. III, it is evident that the two optimization algorithms reach similar results. It is also clear that the optimal geometry is highly dependent on the magnet grade, specifically on the intrinsic coercivity. In this paper, certain parameters are set based on earlier studies using higher graded permanent magnets. Although the "optimal solutions" for Y30 and Y40 by themselves are interesting, the key result of this paper is the need to not only look at the magnetic flux density and output power but also the risk for demagnetization when designing electrical machines using permanent magnets that are more prone to irreversible demagnetization. In regard to the low magnetic flux density reached for Y30, it can, however, be concluded that Y30 is a too weak magnet to work well for this kind of application.