Model Order Reduction for Multi-Strand Electrical Machine Windings

We develop an efficient model order reduction (MOR) method for the finite-element (FE) simulation of an electrical machine with a multi-strand stator winding. The reduced model is independent of the connection of the strands. Therefore, simulation results of only a single slot are required to compute a reduced basis, which can be used for different winding configurations. Numerical simulations of a high-speed induction motor in time-harmonic and time-dependent regimes show significant savings in terms of memory and solving time of the proposed method compared with a standard FE approach.


I. INTRODUCTION
T HE demand to electrify mobility requires efficient and suitable hardware to meet certain requirements, such as high-power output with respect to the size and weight, sustainability, and so on.One trend to meet those requirements is to supply electrical machines with high frequency.Despite the benefits of an increased frequency, this can also cause higher losses, in particular, in windings.To keep the losses at a low level and in order to construct power dense electrical machines, it is important to know already in the design stage where the losses occur, how high they are, and how they are distributed across the conductors the winding is made of.In the mid-frequency range, litz wires provide many benefits and are, therefore, a common choice nowadays to construct windings.For the designer, on the other hand, litz wires cause a prohibitive computational load in the simulations if a classical finite-element (FE) approach is used.Various analytical, numerical, and combined simulation methods have been proposed in the past to estimate losses in litz wires and machine windings made thereof.The developed simulation methods target to mitigate the computational burden due to the required fine FE mesh to resolve the strands in litz wires.If certain assumptions are made, analytical methods allow for fast estimation of skin-effect and proximity losses in litz wires (see [1], [2], [3], [4], [5], and the references therein).According to [6], circulating currents often contribute a significant amount to the winding losses in electrical machines, which add up on top of the skin-effect and proximity losses.The analytical model developed in [7] was capable of estimating losses due to circulating currents in litz wires for frequencies up to 5 kHz in high-power medium-frequency transformers.
Analytical models often require a lot of assumptions to hold or are even oversimplified.Numerical models, on the Manuscript received 12  other hand, are much more flexible and general and can be adopted to estimate the losses in litz wires.To mitigate the computational burden due to the fine mesh in the strands, different approaches have been proposed in the literature.One option is to replace the strands by a uniform current density and compute the losses in a postprocessing step as is done in [8], [9], [10], and [11].These methods are often referred to as hybrid methods, since the flux density obtained from the FE simulation is used in analytical formulas to estimate the losses.
Although such methods are computationally efficient and may catch the proximity and skin effects accurately, they cannot account for the circulating currents, since circuit equations are not solved for individual conductors.In addition, the effect of the eddy currents on the field solution is neglected.Homogenization methods, such as [12], [13], and [14], can include the skin and proximity effects in the field solution, but they also neglect circulating currents.Note that also the previously mentioned approach, where a uniform current density is used, can be interpreted as a homogenization method (see [14]).Another method that is often found in the literature to estimate losses in litz wires is the partial element equivalent circuit method [15], [16].This method can also capture circulating currents; however, this method is only practical for simple geometries, and the equivalent circuit equations introduce large dense blocks in the final system matrix.Approaches, which are capable of computing circulating currents, can be found in [17], [18], [19], and [20], where the authors also provided some applications to electrical machines.In [18] and [20], the slot domains are considered as linear time-invariant systems, and hence, the responses to unit impulses at the boundaries are computed.In [18], the solutions to the unit impulses are collected in snapshot matrices and used as the representations of the basis functions in the slots.In [20], on the other hand, the Lagrange multipliers enforce certain compatibility constraints at the interfaces, and the voltages were expressed by the responses to the input impulses involving convolutions.These expressions were utilized in assembling a reduced linear system that has to be solved for the magnetic vector potential and the currents.This method Fig. 1.Cross section of the considered induction machine with its mesh.has been successfully applied to time-stepping simulations of an induction motor with a stranded winding.In both methods, many local systems need to be solved as there are boundary nodes and strands in a slot, which can be quite large for fine meshes.However, in [18], coupling nodes were introduced to reduce at least the number of solutions related to the unit impulses at the boundary.
In this article, we present a new general, fully numerical method that is capable of estimating the winding losses and computing the exact distribution of the currents within the strands of litz wires in the windings of an electrical machine at a reduced computational cost.To reduce the number of degrees of freedom (dofs) in the slots, we use a local randomized model order reduction (MOR) technique (see [21], [22]).The reduced basis is computed from local solutions of the problem in one slot.The derived local problem is independent of the strand connections, and therefore, the computed reduced basis can be used for all the slots in the machine.The computation of the reduced basis vectors is done in a preprocessing stage, which is very fast.By numerical experiments, we demonstrate that the reduced system can be solved up to 16 times faster than a high-fidelity (HF) FE problem with almost the same accuracy.

II. COMPUTATIONAL MODEL
In this section, we present the computational model of the 2-MW, two-pole, 12-kr/min high-speed induction machine from [23].The model is derived with the aim of accurately estimating the resistive losses in the stator winding, accounting for both skin and proximity effects and the circulating currents in the parallel strands.We use the 2-D AVI formulation, see [24], which couples the magnetoquasistatic magnetic vector potential formulation to circuit equations for solid conductors.For the discretization of the problem, we use Lagrange FEs and polynomial basis functions of first order on a triangular mesh.Denoting the cross section of the solid-rotor induction machine in Fig. 1 by , the AVI formulation after spatial discretization is given by  where the matrices S and M are the usual stiffness and mass matrices, respectively; and D S , C E , and L denote the current source, the back EMF, and the connection matrices, respectively.The diagonal matrix R consists of the strand resistances, and R ew and L ew are the end-winding resistance and end-winding inductance, respectively.The end-winding resistance is estimated as follows: where r slot is the average radius of the slot, τ is the coil pitch angle (τ = π for the full-pitch two-pole winding in this case), ℓ ov is the overhang of the end winding, and R co is the resistance of the winding inside the core region calculated based on the number and area of the conductors.The value for L ew = 7•10 −5 H is taken from [25], where it was estimated by an analytical method.The operator (∂/∂t) denotes the partial derivative operator with respect to time, and I denotes the identity matrix.The vectors a, u, and i ′ contain the coefficients of the discretized out-of-plane component of the magnetic vector potential a, and the voltage drops in the strands and the loop currents, respectively.The voltage source is denoted by the vector u s .The entries of the matrices are given by where ν is the reluctivity, σ is the electrical conductivity, ℓ z is the core length of the electrical machine, and the values of m , m = 1, . . ., n v , denote the cross sections of the strands of the wires.The numbers n, n v ∈ N denote the number of basis functions ϕ i of the discretization space and the number of conductor regions in the cross section, respectively.The number of current loops (number of columns of L) is denoted by n i .Given a domain , we denote by | | the area of .
For the purpose of building the reduced model, we decouple the slots from the rest of the machine and allow that the slots can be modeled independently of the stator.Using the mortar method [26], we ensure that the solution is continuous across the stator-slot interfaces.The mortar conditions are enforced via the matrix B. The matrix B is of size n λ × n, where n λ ∈ N denotes the number of continuity constraints.So, the system we are going to solve is where In the time-harmonic regime, the time derivative in (2) becomes a multiplication with jω, where j is the imaginary unit and ω denotes the fundamental frequency.In order to get a real-valued time harmonic AVI system, we split the real and the imaginary parts of the variables where the entities get the Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.indices R and I , which denote the real and imaginary parts of the entities, respectively.Problem (2) will denote the HF problem in all the remaining sections for both the (real-valued) time-harmonic case and the time-dependent case.

III. COMPUTATION OF THE REDUCED BASIS
In this section, we describe the randomized local MOR method we used to lower the number of field unknowns in the stator slots.We present the method for the real-valued timeharmonic problem derived from (2) and extend the method then to the time-dependent regime in Section V. We consider a subdomain that is the union of slot elements of one slot plus some extra layers of elements around, see Fig. 2, where the green elements extend the slots to a subdomain.The local system for the subdomain is extracted from the main system and is given by where The entries with ∼ on top are obtained by a restriction of the corresponding variables from the main system without ∼.
Note that (3) does not depend on the connection between the strands.Furthermore, we consider S as not being dependent on a, since the material laws in the slot are linear.The two sources that cause a magnetic field in the slot region are excitations gR and ũ I at the boundary of the subdomain and voltage differences ũ R and gI in the strands.The function space where the Dirichlet data gR and gI are taken from should represent all possible restrictions of the solution to the subdomain boundary, independent of the number and location of the slots in the whole machine as well as the connection of the strands.This motivates the use of random, but continuous boundary conditions, which are generated according to [ where the values of σ i are the singular values of X * and ε tol > 0 is a prescribed tolerance.Once the reduced basis vectors for the slots are computed, we define as W ∈ R 2n RED ×2n , where n RED := n + n v + n i − n slots + N R + N I and n slots ∈ N is the number of dofs in the slots, the matrix that is used to replace the FE basis in the slots with the reduced basis vectors.
The reduced matrix and right-hand side for the time-harmonic problem are then obtained from A and f of original system as follows:

NUMERICAL RESULTS FOR THE TIME-HARMONIC PROBLEM
In this section, we report numerical results for the nonlinear time-harmonic problems, which were carried out on a Lenovo ThinkPad P1 with 16 11th Gen Intel Core i7-11850H @ 2.50-GHz processors.The proposed method was implemented in MATLAB, and the meshes were generated using GMSH [27].The conductivity σ was set to 0 in air regions and the laminated stator, 4 MS/m in the solid rotor, and to 25.8 MS/m in the copper regions.The supply voltage and supply frequency were 660 V and 200 Hz, respectively.The rotor slip was set to 0.005%, and 36 parallel strands were modeled in the slot.The subdomain to compute the reduced basis for the slots was the union of the slot elements and two extra layers of elements around the slot.The parameters required to create the boundary conditions were N f = 10, α = 2, and u max = 0.002; see Algorithm 1.The tolerance ε tol for determining the number of reduced basis vectors was 10 −6 .For ũ * = (u * ,1 , . . ., u * ,n v ) ⊤ , we considered two different choices as follows. 1) 2) u * ,1 , . . ., u * ,n v ∈ U(−1, 1).We refer to them as V 1 and V 2, respectively.The distribution U denotes the uniform one over the interval [−1, 1].The HF problem had 443 154 dofs.The HF solution with nonlinear material laws can be seen in Fig. 3.
A Newton-Raphson implementation required around 126 s to solve the HF system.Using a linear local system, it took about 5.13 s to compute the snapshot matrices with 50 random boundary conditions.The reduced systems with 42 402 (V 1) and 45 102 (V 2) dofs were solved in about 22.1 and 70.2 s, respectively.The corresponding relative ℓ 2 errors in the solutions were about 0.01 and 6 • 10 −3 , respectively.Since the model with V 1 is solved three times faster with satisfying error, we continue with this approach.Moreover, numerical experiments indicated that we require about three times more random samples for V 2 to work properly.This is probably the case due to the fact that more variations in the data for the voltages are possible; however, we can see that the flexibility that V 2 provides is not necessary for our problem.
In Table I, we show the solving times (stime) as well as the number of dofs and the density of the matrices for both the HF and the reduced systems for an increasing number of strands per slot.The speedup increases from about 6 to 13 if we increase the number of strands per slot from 16 to 100.The reduced system with 100 strands was about 1/23 times as large as the related HF system, and the relative error in the solution was at about 1.74 • 10 −2 measured in the standard Euclidean norm.The number of required basis functions for the reduced systems for the different numbers of strands was more or less equivalent, as can be seen in Table II.For comparison, we show in Table III the number of dofs and the solving times for a Schur complement method [28], where we eliminated the dofs inside the slots.We can see that the proposed method is up to four times faster than the Schur complement method.
The next set of experiments shows that it is possible that we can build a reduced order model for varying parameters.Without loss of generality we restrict ourselves to variations of the copper conductivity σ c .We considered a machine with 36 strands in each slot.To create the reduced basis vectors, we solved system (3), where σ c ∈ [19.8, 33.8] MS/m was chosen randomly.The rest of the parameters remained the same.We solved the local problem 400 times.The number of   reduced basis vectors was eight for both the real and imaginary parts, and the number of reduced dofs was 42 546.We did five actual simulations with the obtained reduced basis, where we chose σ c ∈ {20.0, 23.0, 26.0, 29.0, 33.0} MS/m.The reduced systems were solved in about 18.0 s, and the errors between the solutions of the HF and the reduced models were between 0.01 and 0.06.

V. EXTENSION TO THE TIME-DEPENDENT CASE
In this section, we extend the procedure from Section III to the time-dependent case.The AVI system for the time-dependent regime with the continuity constraints enforced by the matrix B is given by (2).The backward Euler time-stepping scheme is applied.We refer to (2) as the HF problem in the time-dependent setting.The procedure to get the reduced basis vectors is similar as above; i.e., we consider one subdomain, which is a slot plus some extra layers of mesh elements.The local problem we have to solve is given by where we use analogously as above Algorithm 1 to obtain g and the variants V 1 for ũ.The procedure to compute the reduced basis vectors for the slot regions is summarized in Algorithm 2. Analogously, as in Section III, for the timeharmonic case, we define a matrix W ∈ R n RED ×n , where, in this section, n RED := n + n v + n i − n slots + N , where n slots is defined as in Section III and N is the number of reduced Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.basis functions in the time-dependent case computed by ( 4), that replaces the basis in the slots.The reduced time-dependent problem is obtained from (2) as follows: where Algorithm 2 Compute a Reduced Basis for the Time-Dependent Problem Require: time interval (t s , T ], where T > t s ≥ 0, t h = (t 0 , . . ., t n t ), where t 0 = t s and t n t = T and .) X = ∅ 3.) Compute g according to Algorithm 1 4.) Draw ũ ∈ U(−1, 1) n v for t = 1, . . ., n t do 5.) Solve (5) for ãt and set X ← X ∪ ãt end for 6.) Compute the SVD decomposition X = U V T 7.) Get the most significant singular values N via (4) 8.) Y ← N most significant left singular values of X end for 9.) Compute the SVD decomposition Y = U V T 10.) Get the most significant singular values N via (4) 11.) R ← N most significant left singular values of Y

VI. NUMERICAL RESULTS FOR THE TIME-DEPENDENT PROBLEM
In this section, we present numerical results for the timedependent problem.The hardware and the parameters were almost the same as already stated in the beginning of Section IV.We considered 120 time steps per fundamental period and initialized the simulation with the solution of the nonlinear time-harmonic problem.As seen in Fig. 4, the supply voltage was gradually changed from sinusoidal to a rectangular wave corresponding to converter supply.Three periods with n t = 360 time steps in total were simulated, and all the results were recorded during the last period.
We report the performance of the method only for V 1 but for a varying number of parallel strands.We used 16, 36, 64, and 100 parallel strands for our numerical experiments.
In Table IV, we show the number of dofs, the solving times, and the error of the solutions between the HF and reduced models.The error e is measured in the norm where ) ⊤ is the solution vector without the Lagrange multipliers of the HF problem (2) at time step i and y H,i is computed from W x RED , with x RED being the solution to the reduced problem (5), and removing the coefficients for the reduced λ.The norm ∥•∥ ℓ 2 is the standard Euclidean norm.
We observe from Table IV that the number of dofs can be reduced significantly, and the solving times for the reduced problems drop by factors between 2 and 4. The errors were about the same order for every number of strands.The number of basis vectors for each case can be seen in Table II.This table shows that the solutions in the slot over all the timesteps can also be represented by a low number of basis functions.
In the following, we compare the phase currents and the winding losses for the model with 100 strands per slot.The top graph of Fig. 5 shows a very good agreement between the phase currents of the HF and the reduced model.In the bottom graph, we see the variations of the strand currents, represented by the gray shaded areas, and the average strand currents of phase a, represented by a blue solid line.The eddy-current losses computed from the solutions of the HF and reduced problems were about 47.3 and 46.3 kW, respectively.So, the losses computed from the solution of the reduced problem underestimate the losses from the HF models only slightly.

VII. IMPLEMENTATION CONSIDERATIONS
In this section, we describe the implementation steps of the proposed method, which can be used in any commercial software package.The implementation of the MOR method can be divided into three steps: preprocessing, solving, and postprocessing.In the following, we describe the implementation steps for the time-stepping method.

A. Preprocessing
This step considers the setup of the reduced basis.1) Use Algorithm (2) to obtain the reduced basis for one slot.

B. Solving
In this step the reduced linear system is solved.1) Use the reduced basis to project the HF data in the slots into a reduced space spanned by the reduced basis functions.2) Solve (5) to get the reduced solution x RED .

C. Postprocessing
The solution of the reduced model must be prolongated into the HF space to compute the losses.
1) Prolongate the reduced solution into the HF space by computing x = W x RED .
2) Use x to compute the losses P ed by where t * and t * are the start and end times of the last period.

VIII. COMPARISON WITH PREVIOUS APPROACHES
The MOR methods proposed in [18] and [20] share some similarities to the suggested one in this article.In [18], the snapshot matrix containing the reduced basis vectors can be quite large due to a large number of boundary nodes and a large number of strands.The method in [20] uses the snapshot matrices to compute convolutions involving the inputs and responses to assemble and solve a reduced AVI problem.In both approaches, the problems related to the slot need to be solved consecutively, first with a unit impulse at each boundary node and then with a unit current for each strand.In our method, both field sources in the slot are set at once, and we demonstrated that we need to solve the local problem fewer times as there are boundary nodes.Moreover, the proposed method is capable of constructing vectors representing a reduced basis where material parameters may vary, as is demonstrated for the time-harmonic case, where we varied the electric conductivity.This means we have to compute only one reduced basis, which can be used for systems with different material parameters.This property is not obvious from the methods in [18] and [20] and would probably require some modifications of the algorithms.

IX. CONCLUSION AND OUTLOOK
In this article, an efficient MOR method is proposed to reduce the dofs in multi-strand electrical machine windings.Given that suitable parameters are used, the proposed method requires only little offline computation time and provides an accurate reduced model with good speedups in the online stage.The proposed method is independent of the strand connections.Hence, using only one subdomain to compute a reduced basis is sufficient for the whole machine.In particular, once a reduced basis is computed, it can be used to solve problems with different connection matrices L, provided that the number of strands in the slots remains the same.Furthermore, we compared the phase currents from the HF and the reduced models as well as the losses in the strands.The results show that the reduced model is indeed capable of providing a good estimate for the losses in the machine.The mortar method allows to create the slot geometry independent of the stator geometry; i.e., hanging nodes at the stator-slot interface do not cause any issues, and a reduced basis for the slot can be computed independent of the final rotor and stator geometries of the machine.One possible approach for modeling machines with twisted or transposed windings is to combine the proposed method with a multi-slice model where the 2-D problem is solved simultaneously in several slices along the axial direction [29], [30].

Fig. 2 .
Fig. 2. Discrete subdomain, conductors in orange, insulation and air regions in white, extra layers for the subdomain in green, and the surrounding stator material in gray.
22, Algorithm 2].To keep this article self-contained, we briefly recap the algorithm in Algorithm 1 for generating the boundary data.Since we do not make any assumptions about the voltage differences in the strands, we use random vectors for ũ R and ũ I , where the entries are sampled from a uniform distribution over the interval [−1, 1].The local system (3) is solved n S ∈ N times with n S different Dirichlet data gR , gI and voltage differences ũ R , ũ I .The local solutions ã R,i and ãI,i , i = 1, . . ., n S are collected in snapshot matrices X R = ( ã R,1 , . . ., ã R,n S ) and X I = ( ãI,1 , . . ., ãI,n S ), respectively, from which singular value decompositions are computed.The vectors for the basis of the reduced space are obtained by using the first N R and N I left singular vectors of the matrices X R and X I , respectively, where N R , N I ∈ N are determined by

Fig. 3 .
Fig. 3. Real part of the flux density of the nonlinear time-harmonic problem.

Fig. 5 .
Fig. 5. Above are the phase currents of the HF model (x marks) and the reduced model (solid lines) for 100 parallel loops, estimated losses: 47.3 and 46.3 kW, respectively.Below is the variation of the strand currents computed with the MOR model of phase a (gray shaded area) and the average current (blue solid line).

TABLE I NUMBER
OF DOFS, SOLVING TIMES (STIME), AND THE ERROR IN THE SOLUTION BETWEEN THE HF AND THE REDUCED MODELS IN THE TIME-HARMONIC CASE

TABLE II NUMBER
OF REDUCED BASIS FUNCTIONS FOR THE TIME-HARMONIC AND TIME-DEPENDENT PROBLEMS

TABLE III NUMBER
OF DOFS AND SOLVING TIMES (STIME) OF THE

TABLE IV NUMBER
OF DOFS, SOLVING TIMES (STIME), DENSITIES OF THE JACOBIANS, AND THE ERROR IN THE SOLUTION BETWEEN THE HF AND THE REDUCED MODELS IN THE TIME-STEPPING CASE