Spectral element method modeling of eddy current losses in high-frequency transformers

This paper concerns the modeling of eddy current losses in conductive materials in the vicinity of a high-frequency transformer; more specifically, in two-dimensional problems where a high ratio between the object dimensions and the skin-depth exists. The analysis is performed using the Spectral Element Method (SEM), where high order Legendre–Gauss–Lobatto polynomials are applied to increase the accuracy of the results with respect to the Finite Element Method (FEM). A convergence analysis is performed on a two-dimensional benchmark system, for both the SEM and FEM. The benchmark system consists of a high-frequency transformer confined by a conductive cylinder and is free of complex geometrical shapes. Two different objectives are investigated. First, the discretizations at which the relative error with respect to a reference solution is minimized are compared. Second, the discretizations at which the trade-off between computational effort and accuracy is optimized are compared. The results indicated that by applying the SEM to the two-dimensional benchmark system, a higher accuracy per degree of freedom and significantly lower computation time are obtained with respect to the FEM. Therefore, the SEM is proven to be particularly useful for this type of problem.


Introduction
Wireless Power Transfer (WPT) by means of an inductive coupling is widely applied in applications where physical electrical contact is problematic or undesirable, for example in aerospace, biomedical, and robotics applications [1]. To improve the performance, a high-frequency power-supply is often applied. Therefore, gallium-nitride transistors have gained popularity in WPT applications, because these transistors allow for switching frequencies in the range of several MHz, low switching losses, and a high power density [2,3].
Adding a high-frequency power-supply to an application imposes several challenges in the electromagnetic modeling, especially in situations where restrictions on the field penetration into surrounding objects are of importance; for example, in solid rotor, as well as permanent magnet high-speed machines [4][5][6] and in order to comply with electromagnetic compatibility safety standards [7]. The Finite Element Method (FEM) is often applied for the estimation of the resulting losses [8,9]. For high-frequency applications in which skin-depths are orders of magnitude smaller in comparison to the overall model dimensions, meshing problems arise, which lead to locally different requirements of the mesh density within the model [5]. In order to minimize the mesh size, elongated mesh elements are preferred, which compromise the accuracy of the FEM. Alternatively, for this type of problem, higher order methods can be applied, such as the Spectral Element Method (SEM), in which elongated elements do not compromise the accuracy of the results [5,10].
The SEM has been widely applied in various research areas, such as fluid dynamics or waveguide analysis [11,12]. In [5], a single spectral element was used to model the eddy current losses in a high-speed rotating cylinder. In [10], the formulation of the SEM for the magnetic field and temperature distribution in electric machines was discussed, and in [13], the application of the method to a linear synchronous machine was presented. In [14], a coupled FEM and SEM simulation was introduced, where the non-linear SEM was used for modeling of the rotor and the FEM was applied for the rest of the model. Another example of a coupled approach was provided in [15], where transformer cooling was investigated using coupled FEM and SEM. A detailed discussion on the implementation of the SEM was presented in [16], and the mathematical formulation of eddy current problems was provided in [17][18][19]. The modeling approach presented in this paper is valid for two-dimensional problems. However, the method can be extended to three-dimensional problems. In various research areas, such as computational seismology, nanodevice applications, and waveguide analysis, the SEM has been successfully implemented and applied to three-dimensional problems [12,20,21].
In this paper, it is demonstrated that by applying the SEM to a benchmark system, a higher accuracy per degree of freedom (d.o.f.) is obtained, which results in a significant reduction of the computational effort as compared to the FEM. These features prove to be useful for efficient and accurate solving of eddy current loss problems with high ratios between the object dimensions and skin-depth. The formulation of the method, as well as the implementation of the specific transformer model are discussed and analyzed.

Modeling Approach
To model the eddy current problem, the formulation of the SEM for a generic partial differential equation is considered: where c is a material property, ϕ is the physical quantity of interest, and s is the source descriptions. In order to model eddy currents, s takes the form of: where k is a different material property. In order to solve the differential equation using the SEM, the investigated geometry is divided into rectangular elements. In each element, the nodal Galerkin method is applied to (1), meaning that the residuals of the approximated solution are projected onto a set of basis functions, therefore obtaining the weak form. After applying Green's first identity, the following expression is obtained: where β is the basis functions and Ω ξ,η is the computational domain of the element with local coordinate axes ξ and η [10]. The basis functions applied in the SEM are based on the Legendre polynomials, which are given by: in which, L 0 = 1 and L 1 = ξ (4) where L N is the Legendre polynomial of degree N. The degree of the Legendre polynomial directly relates to the number of grid points in an element. The grid points correspond to the roots of the Legendre-Gauss-Lobatto (LGL) polynomials, which for degree N are given by: where P LGL are the roots, calculated through an iterative process. The integrals shown in (3) are evaluated numerically by Gaussian quadratures ( where f j is the value of the function f at the corresponding root (ξ j ). The Lagrangian interpolation polynomials are used to reconstruct the polynomial basis and interpolate the approximated function: where l j is the Lagrangian basis function given by: The derivative of the basis functions from (3) is also computed on the roots using the Lagrangian basis functions from (7), which results in a derivative matrix. The entries of the derivative matrix are computed as: where l j is the derivative of the Lagrangian basis function. The derivative matrix replaces the derivative operators in (3). Finally, the system of linear equations for each element is obtained and assembled in a global matrix [10,16]. For the modeling of eddy current losses in conductive materials in the frequency domain, the partial differential equation shown in (1) takes the form of: where ν is the magnetic reluctivity, A θ and J θ are the vector potential and current density in the circumferential direction, respectively, ω is the electrical frequency, and σ is the electric conductivity. The eddy current losses in a conductive region are calculated as: where J e is the induced current density and V is the volume. The integral in (11) is evaluated numerically by the Gaussian quadratures of the corresponding region.

Benchmark System
The benchmark system considered in this paper consists of a pot core transformer, which is operated at an electrical frequency of 1 MHz. The transformer is confined by a conductive cylinder, which is directly positioned in the fringing-flux path parallel to the air gap. As a result, eddy current losses are induced in the conductive cylinder. An overview of the investigated domain in the three-dimensional space is shown in Figure 1, and the various components of the geometry are indicated in the figure. In each side of the core, a small indentation is present, such that the winding is able to enter and exit the magnetic core. The effect of the indentation on the transferred power is negligible [22]. Therefore, the benchmark system is modeled as an axisymmetric problem. An overview of the investigated domain in the two-dimensional space is shown in Figure 2, and the corresponding geometrical parameters and physical quantities are shown in Table 1. For this type of problem, (3) is typically solved for rϕ, as opposed to directly solving for ϕ [23,24]. Consequently, a negligible offset (r 0 ) is present in the geometry, such that the zero-axis is excluded from the domain. The investigated domain is bounded by a zero Dirichlet boundary condition. The benchmark system is modeled using both the SEM and FEM. A comparison between the two methods is made by performing a convergence analysis and refinement optimization.
Mesh density points

SEM Model
The division of the domain into rectangular elements and the polynomial degree discretization for the SEM model are shown in Figure 2a. The polynomial degree in the domain is discretized into three different parameters. One parameter (N z ) is assigned to the axial (z) direction, whereas two parameters (N r1 and N r2 ) are assigned to the radial (r) direction. A distinction between non-conductive and conductive regions is made because in the conductive region, the magnetic vector potential is expected to have a higher gradient compared to the rest of the geometry. Therefore, a high polynomial degree is required for a locally more accurate approximation.

FEM Model
The mesh density in the FEM model is parametrized into five different minimum element sizes (mesh points). The smallest mesh points (M s1 , M s2 , and M s3 ) are assigned to the conductive cylinder and air gap regions, whereas the largest mesh point (M l ) is assigned to the domain boundary. Furthermore, an intermediate mesh point (M m ) is assigned to part of the transformer core and winding area. The division of the mesh points is shown in Figure 2b. The FEM model is solved using commercial software (Altair Flux) [24]. Triangular mesh elements are created in the non-conductive regions, whereas rectangular mesh elements are employed in the conductive region. The FEM requires at least two elements per skin-depth in the radial direction of the conductive region, such that the second order elements provide an accurate approximation of the exponential variation occurring in this region [24]. Therefore, rectangular elements are preferred because the contrasting mesh requirements in the radial and axial direction can be realized by creating elongated elements. As a result, a lower number of d.o.f. is required for meshing the conductive region compared to a mesh consisting of triangular elements [24]. In the FEM model, two mesh elements per skin-depth are ensured regardless of the size of the mesh elements in the rest of the model.

Assumptions
Both the transformer core and conductive cylinder are assumed to have linear material properties. The current density values and corresponding phase angles were obtained from an FEM simulation, in which the primary side was supplied with a peak voltage of 48.0 V, the secondary side was connected to a load of 23.5 Ω, the primary side had four turns, and the secondary side had seven turns. The current densities and phase angles were equal for both methods and remained unchanged in the convergence analysis.
In the interest of significantly reducing the mesh requirements in the conductive cylinder, model reduction was applied in the cylinder; only a radial thickness equal to three-times the skin-depth was taken into account. The effect on the induced losses is shown in Figure 3, in which the absolute value of the relative error with respect to the reference, i.e., the loss for the full radial thickness of the conductive cylinder, is shown as a function of the modeled radial thickness. The absolute value of the relative error with respect to the reference ( ) is given by: where P e is the calculated eddy current loss and P * e is the reference. The results were generated by an FEM simulation, in which at every iteration, the modeled radial thickness of the conductive cylinder was increased by one skin-depth. In case the modeled radial thickness was equal to three-times the skin-depth, a relative error with respect to the reference of less than 0.1% was obtained, therefore justifying the assumption.

Convergence Analysis and Refinement Optimization
The two methods were compared to each other by performing a convergence analysis on the benchmark system. In the convergence analysis, p-refinement was applied in the SEM model, i.e., the accuracy of the results was improved by increasing the polynomial degree. The p-refinement was applied according to the discretization shown in Table 2, which resulted in a total of 9000 combinations being evaluated and a number of d.o.f. in the range of 56-17,091. In the FEM model, h-refinement was applied, i.e., the accuracy of the results was improved by reducing the characteristic length of the mesh elements, resulting in the number of elements being increased. The discretization of the h-refinement applied in the convergence analysis is shown in Table 2, which resulted in a total of 34,300 combinations being evaluated and a number of d.o.f. in the range of 937-27,528. The coarsest and densest FEM meshes occurring in the convergence analysis are shown in Figures 4 and 5, respectively. For visualization purposes, the mesh is only shown in the most critical part of the model, i.e., the region where the fringing-flux around the outer transformer teeth induces the eddy current loss. This region is also indicated in Figure 2b by the blue dotted rectangle. Furthermore, the mesh in the non-conductive and conductive regions are shown in separate figures. For every combination of por h-refinement evaluated in the convergence analysis, the number of d.o.f. and the losses in the conductive cylinder were obtained. Furthermore, the absolute value of the relative error with respect to the reference solution was calculated using (12), where the FEM solution having the highest number of d.o.f. was chosen as the reference. Additionally, the time required for meshing and solving of the models was stored. Both models were evaluated on a single core of the same machine using a direct solving technique. Two Pareto fronts were generated from the obtained solutions. The first front consisted of the absolute relative error with respect to the reference and the number of d.o.f. In the second front, the latter was changed to the computation time (consisting of the time required for meshing and solving of the models).
The optimum discretization of the polynomial degree and mesh points was determined by applying the weighted sum method, given by: minimize: where: {w 1 , w 2 } ∈ (0, 1) subject to: where x is the set of design variables (the polynomial degree and mesh points),f n,t is the normalized number of d.o.f. in the first case, or the computation time in the second case, a weighting factor of w 1 was applied, andf is the normalized absolute value of the relative error with respect to the reference having weighting factor w 2 . The quantities were normalized with respect to the minimum and maximum values occurring in the dataset generated by the corresponding method. In this analysis, two different situations were investigated. First, the weighting factors were set to zero and one, respectively, which resulted in the converged solution. Second, both weighting factors were set to 0.5, which resulted in an equally-weighted trade-off between the absolute value of the relative error with respect to the reference and the number of d.o.f. for the first front, or the computation time for second front, respectively. Table 2. Discretization of the polynomial degree (SEM) and mesh points (FEM) evaluated in the convergence analysis.

Results
The absolute value of the relative error with respect to the reference as a function of the number of d.o.f. and the computation time for all evaluated data points is shown in Figure 6a,b, respectively. Additionally, the Pareto fronts and optima are indicated in the figures. The figures show that by applying the SEM, a higher accuracy per d.o.f. and a lower computation time were obtained compared to the FEM.
The converged and Pareto solutions for both methods are shown in Table 3. As a result of the FEM solution having the highest number of d.o.f. being the reference, the FEM was able to converge to zero; whereas, the SEM converged to a relative error with respect to the reference, approaching zero, while reducing the required number of d.o.f. and computation time by 76.5% and 75.1%, respectively. Additionally, an FEM solution is shown in Table 3, for which the difference between the relative errors was the smallest. In this case, the SEM reduced the number of d.o.f. and computation time by 62.2% and 70.0%, respectively, therefore proving the advantage of the SEM for this type of problem.
The solutions at the Pareto optima are also shown in Table 3. For the first Pareto optimum (the trade-off between the number of d.o.f. and the absolute value of the relative error), the SEM reduced the number of d.o.f. and computation time by 92.3% and 99.2%, respectively, whereas the relative error with respect to the reference was 0.604% higher. Compared to the converged solution, the number of d.o.f. and computation time were significantly reduced for both methods at the expense of an increased relative error with respect to the reference. The second Pareto optimum (the trade-off between the computation time and the absolute value of the relative error) presented a more attractive solution. For both the SEM and FEM optima, the error with respect to the reference solution was very small (less than 0.1%), while the number of d.o.f. and computation time were significantly reduced. Comparing the two solutions, both methods had an approximately equal relative error with respect to the reference, while the SEM reduced the number of d.o.f. and computation time by 79.5% and 98.2%, respectively. The discretization of the SEM and FEM model at the optima is shown in Table 4.   Table 4. Discretization of the polynomial degree (p-refinement) and mesh points (h-refinement) in the SEM and FEM model, respectively, at the two Pareto optima.

Conclusions
The SEM has been applied for the modeling of eddy currents in conductive materials in the vicinity of a high-frequency transformer. The performance of the method has been investigated by performing a convergence analysis on a two-dimensional benchmark system, consisting of a high-frequency transformer confined by a conductive cylinder. The results have indicated that the SEM is able to provide a fast and highly accurate estimation of the eddy current losses (less than 8.0 × 10 −4 % relative error with respect to the reference), while requiring 70.0% less d.o.f. and 62.2% less computation time with respect to the FEM. Therefore, the SEM provides a better alternative for this type of problem, i.e., two-dimensional problems where the skin-depth is several orders of magnitude smaller compared to the object dimensions and complex geometrical shapes are absent.

Author Contributions:
The results presented in this paper have been developed by K.B. and M.C. The analysis has been performed in cooperation with D.C.J.K, S.J. and E.A.L. The paper was written by K.B., and contributions and improvements to the content have been made by M.C., D.C.J.K, S.J. and E.A.L.

Conflicts of Interest:
The authors declare no conflict of interest.