Spectral Element-Based Multi-Physical Modeling Framework for Axisymmetric Wireless Power Transfer Systems

This paper concerns a multi-physical modeling framework based on the spectral element method (SEM) for axisymmetric wireless power transfer systems. The modeling framework consists of an electromagnetic and a thermal model. The electromagnetic model allows for eddy currents in sourceand non-source regions to be included in the analysis. The SEM is a numerical method, which is particularly advantageous in 2D problems for which the skin-depth is several orders of magnitude smaller compared to the object dimensions and complex geometrical shapes are absent. The SEM applies high-order trial functions to obtain the approximate solution to a boundary-value problem. To that end, the approximation is expressed as an interpolation at a set of nodal points, i.e., the nodal representation. The trial functions are Legendre polynomials, which reduces the complexity of the formulation. Furthermore, numerical integration is performed through Gaussian quadratures. In order to verify the SEM, a benchmark system is modeled using both the SEM and a finite element-based commercial software. The differences in the SEM solutions, i.e., magnetic vector potential and temperature distribution, and the discrepancies in essential post-processing quantities are assessed with respect to the finite element solutions. Additionally, the computational efforts of both methods are evaluated in terms of the sparsity, number of degrees of freedom, and non-zero elements.


Introduction
Wireless power transfer (WPT) by means of an inductive coupling is a popular and mature technology. Correspondingly, inductive WPT systems are widely applied in applications that require reliable and durable power transfer to rotary parts, for example, robotics and battery charging applications [1,2], as well as electrical machines and actuators [3,4]. Such a system generally employs a cylindrical transformer, which consists of a stationaryand rotary side separated by a (small) air gap. For the purpose of achieving a high magnetic coupling, a magnetic core, which provides a low reluctance path to the magnetic flux, is often added on either or both sides of the system. Furthermore, a high electrical frequency is generally applied, which benefits the efficiency and reduces the volume [5]. Consequently, gallium-nitride (GaN) transistors are emerging in WPT applications, since frequencies in the range of several MHz are realized compactly and efficiently [6]. The electrical frequency of a WPT system based on an inductive coupling typically ranges from tens of kHz up to several MHz [7]. In order to accommodate for high electrical frequencies, the magnetic core material is generally ferrite, which, as a result of the low electrical conductivity, has low losses in a wide frequency range, i.e., up to 3 MHz for energy conversion applications. Furthermore, the coils typically consist of litz wire, such that the losses caused by the skinand proximity-effect are minimized. Alternatively, foil windings are often applied in case the required effective copper cross-section is high [8]. In small transformers, foil windings are popular as a result of the simplicity of the winding operations [9]. is desired. For this purpose, the SEM is applicable as well. In this case, the heat sources and temperature-dependent material properties can be enforced on the same grid, i.e., without the need for interpolation, within the same software environment. Alternatively, an interpolation step can be included, such that the models can apply a different degree of the approximation, which potentially reduces the computational effort. The formulation of the SEM for 2D Cartesian steady-state thermal problems is introduced in [27]. However, in order to apply axisymmetry, modifications are necessary. Succinctly, a multi-physical modeling framework for axisymmetric WPT systems based on the SEM is an imminent research topic. In particular, in case the analysis involves (high-frequency) eddy currents in both source and non-source regions.
Multi-physical modeling frameworks are emerging in a wide variety of research areas, e.g., earth and nuclear science, fluid-structure interaction, selective laser melting, and fuel cells [34][35][36][37][38][39]. In earth science, non-linear phenomena occur that result from multi-physics coupling across multiple scales from the quantum level to the scale of the earth. Moreover, the time scale varies in the range of femtoseconds to billions of years. In [34], a first overview is presented of the current approaches to coupling the scales using methods borrowed from other disciplines. Two trends in the approaches to this type of problem have been reported. The first trend is a statistical mechanics-based upscaling approach that considers discrete energy interactions in a hierarchical nesting of calculations that can be seen as the equivalent of quantum mechanical simulations extended to larger scales. In the second trend, the problem is approached from a thermodynamically-based continuum framework. With regard to nuclear science, in [35], the Picard and Jacobian-free Newton-Krylor (JFNK) methods are developed based on the high-temperature gas-cooled reactor simulator TINTE, which is widely used to analyze the transient behavior of such a reactor. The numerical results have indicated that both algorithms, as a result of the accuracy and stability advantages, achieve a higher computational performance with respect to the original semi-implicit coupling algorithm in TINTE. Furthermore, in the field of nuclear science, ref. [36] employs a multi-physics coupling of SCALE/TRACE for neutronic analysis and spent fuel validation of boiling-water reactors spent fuel isotopics. The methodology has been applied to two boiling-water reactor assemblies, discharged from the Fukushima Daini-2 unit. The neutronic analysis on the assembly level has indicated that the coupling process is verified. Moreover, validation of the coupling process has been carried out through spent fuel isotopics data, which demonstrated good results for uranium isotopes and satisfactory results for other isotopes. On the subject of fluid-structure interaction, ref. [37] introduces a partitioned Newton-based method for solving non-linear coupled systems arising in the numerical approximation of such problems. The work proposes a new strategy to implement linearized solvers for the evaluation the various Jacobians involved in Newton's algorithm. The main contribution lies in the establishment of the expressions for these Jacobians. With regard to selective laser melting, ref. [38] introduces a mesoscale discrete element method and computational fluid dynamics combined simulation framework for the simulation of such processes. The application of the framework has demonstrated a successful layer-by-layer simulation. Furthermore, multi-physical modeling frameworks are applied in the research area concerned with fuel cells. In [39], a multi-physical and electrochemical coupling model has been established for a protonic ceramic fuel cell. Experimental validation of the model has been carried out and a good agreement between the simulation and experimental results has been obtained.
In this paper, a multi-physical modeling framework based on the SEM for axisymmetric WPT systems is presented. The framework includes an electromagnetic and a thermal model, of which the former incorporates eddy currents in both source and non-source regions. The formulation and implementation of both models are discussed and verified. Furthermore, the computational effort is evaluated. The modeling framework is well-suited for the design analysis of (high-frequency) axisymmetric WPT systems, in which complex geometrical shapes are absent and various electrically conductive materials are present.

Benchmark System
For the purpose of establishing, applying, and verifying the multi-physical modeling framework based on the SEM, a benchmark system is considered, which is representative of a WPT system for a rotary application. An overview of the benchmark system in 3D space is shown in Figure 1. The primary and secondary sides of the system both consist of a ferrite magnetic core, in which a copper foil winding is situated. The transformer is fully enclosed by a stainless-steel shaft and an aluminum housing, which acts as a shield. The electrical frequency is 1 MHz, thus (high-frequency) eddy current losses are induced in all electrically conductive materials, i.e., the coils, the aluminum housing, and the stainlesssteel shaft. Furthermore, hysteresis losses are present in both ferrite magnetic cores. The eddy current losses in the magnetic cores are neglected, since the electrical conductivity of ferrite is negligible compared to aluminum and stainless-steel. As a result of the various loss components, the temperature of the domain rises. The aluminum housing transfers heat by means of convection and radiation to an infinitely large surrounding. The cylindrical geometry of the benchmark system, shown in Figure 1, is reduced to an axisymmetric problem, such that the computational effort is significantly reduced. In order to obtain the electromagnetic field solution and temperature distribution with the SEM, the investigated domain is divided into multiple elements of matching material properties. Furthermore, the boundary conditions are applied, such that a boundaryvalue problem is obtained. In the case of the electromagnetic model, the investigated domain is bounded by a zero Dirichlet boundary condition. In the thermal model, both convection and radiation to an infinitely large surrounding are present on the domain boundary, thus a Robin boundary condition is applied. In both cases, continuous boundary conditions are applied to the interior boundaries between the elements. An overview of the benchmark system with axisymmetry is shown in Figure 2. The various boundary conditions, geometrical parameters, and physical quantities for both the electromagnetic and thermal domain are indicated in the figure. The corresponding numerical values of the geometrical parameters and physical quantities are given in Table 1. Unless otherwise specified, the relative permeability is equal to one. A negligible offset is present in the geometry, such that the zero-axis is excluded from the domain and singularities in the solution are avoided. The design of the ferrite cores present in the benchmark system has been derived in [16]. However, in this paper, the primary and secondary coil consist of foil windings, which, for simplicity, employ two turns on both sides of the system.

Boundary conditions
Material properties Figure 2. Overview of the benchmark system with axisymmetry, including indication of the boundary conditions, element discretization, geometrical parameters, and material properties in both the electromagnetic and thermal domain.

Global Overview of the SEM
In this paper, an electromagnetic and a thermal model are established based on the SEM for axisymmetric WPT systems, such as the benchmark system, which constitutes the multi-physical modeling framework. To that end, the equations governing the electromagnetic and thermal behavior of the benchmark system are introduced in Sections 2.3 and 2.4, respectively. Furthermore, in these sections, the numerical approximations of the governing equations using the SEM are discussed. Essentially, the SEM consists of six main processes, which are summarized as follows: 1.
An initialization step is performed, in which the nodes where the field variable is approximated are derived. Furthermore, the matrices required to compute integrals and derivatives are derived, which are essential components for the approximation of the governing equations. These aspects are further detailed in Section 2.3.1.

2.
The higher-order trial functions employed by the SEM to approximate the solution, achieve the optimal convergence rate in a [−1, 1] 2 square domain. The approximation of the governing equations for such a domain is discussed in Section 2.3.2. In reality, such a domain is scarcely ever observed. Therefore, a mapping is applied, such that the convergence rate is maintained, while the geometrical flexibility of the method is greatly enhanced. In Section 2.3.3, the mapping functions are presented and the weak form of the governing equations under mapping is given for the electromagnetic model. In the case of the thermal model, the latter is given in Section 2.4.1, whereas the mapping functions are equivalent.

3.
In order to solve the linear system of equations, the problem is expressed in its equivalent matrix form, which is further detailed in Sections 2.3.4 and 2.4.1 for the electromagnetic and thermal model, respectively. 4.
The boundary conditions are applied, which is the last step before the system of equations is solved. The various boundary conditions occurring in an axisymmetric WPT system are discussed in Sections 2.3.5 and 2.4.2 for the electromagnetic and thermal model, respectively.

5.
The linear system of equations is solved, yielding the solution of the field variable at the nodes in the investigated domain. 6.
Finally, various post-processing calculations are performed, e.g., to obtain the various loss components or the average temperature. The post-processing calculations, relevant to axisymmetric WPT systems, are discussed in Sections 2.3.6 and 2.4.3 for the electromagnetic and thermal model, respectively.
A flowchart of the six main processes is shown in Figure 3. The flowchart and the processes are further detailed in Section 2.5. The multi-physical modeling framework, consisting of the electromagnetic and thermal model, is verified by comparing the solutions from the benchmark system to a reference solution, which has been obtained with a commercial FEM software package. The method of verification is further discussed in Section 2.6, whereas the results are presented in Section 3.

6.
Start Figure 3. Flowchart of the six main processes in the SEM.

Electromagnetic Formulation
The benchmark system, shown in Figure 2, consists of three distinctive regions, which are characterized according to either the electrical conductivity or imposed current source:

1.
The electrically conductive source regions, i.e., the copper foil windings situated in the primary and secondary core.

2.
The electrically conductive non-source regions, i.e., the aluminum housing, which surrounds the design space, and the stainless-steel shaft.

3.
The electrically non-conductive non-source regions, i.e., the ferrite cores and air gaps.
Consequently, the electromagnetic behavior of the benchmark system is governed by three different partial differential equations (PDEs). The imposed current sources are assumed to be sinusoidally time-varying, thus the PDEs are represented in the frequency domain. The formulation and implementation of the SEM introduced in this paper is based on [27,30,40], in which primarily 2D problems in the Cartesian coordinate-system of various physical domains are considered, e.g., magneto and electro quasi-static. Therefore, the transformation discussed in [41] is applied to all three PDEs, such that the Del operator acting on the field variable is evaluated as if the coordinate-system were Cartesian. Consequently, the formulation and implementation of the SEM discussed in [27,30,40], is conveniently modified to problems with axisymmetry. As a result of the transformation, a change of the field variable is introduced, which is given by: where A θ andÃ θ are the magnetic vector potential and the modified variant in the azimuthal direction, respectively, and r and z are the radial and axial coordinate, respectively.
In the first class of regions, an additional DoF must be incorporated, which forces the average value of the eddy current density to be zero. The derivation of this additional DoF for axisymmetric problems is discussed in [42,43]. In conjunction with the transformation from [41], the governing PDE for the first class of regions is given by: where ∇ c is the Del operator in the Cartesian coordinate-system, where (r, z, θ) corresponds to (x, y, z), ν is the magnetic reluctivity, i.e., the inverse of the permeability, J θ,0 is the imposed current density, ω and σ are the electrical frequency and conductivity, respectively, χ Ã θ is the additional DoF, I is the imposed current, and Ω is the physical domain of the respective source region. In the second class of regions, (3) and (4) are always equal to zero, hence, in such a region, (2) reduces to: Finally, in the third class of region, neither a source nor an electrical conductivity are assigned, thus (5) reduces to: Succinctly, the electromagnetic behavior of the benchmark system is governed by (2)-(6). As a result of applying the transformation from [41], the inverse of the radial coordinate is present in all three PDEs. Consequently, in order to avoid singularities at the symmetry axis, i.e., r = 0, a negligible offset has to be added to the physical domain.

Spectral Approximation
In all spectral methods, the approximation is based on the expansion of a function in terms of an infinite sequence of orthogonal trial functions, which is given by: where ψ(x) is the function to approximate, k is the number of expansions,Φ k are the expansion coefficients, and ϕ k (x) is an orthogonal trial function. The accuracy of spectral methods increases for an increasing number of trial functions applied to (7). The series converges, i.e., N → ∞, in case the infinite set of trial functions, {ϕ k (x)} ∞ k=0 , forms a basis for the space of functions that are square integrable in the interval [a, b] with respect to a weight w(x), which is typically represented by L 2 w (a, b). The most common example is the expansion of periodic functions in the Fourier series. The trial functions are periodic on the interval [0, 2π] and are orthogonal with respect to the weight function w(x) = 1. Therefore, a basis for L 2 (0, 2π) is obtained and any square integrable function can be represented by (7) [40]. For the purpose of providing an approximate solution, the infinite series is truncated. In case the function ψ is infinitely smooth and both the function and all its derivatives are periodic, the discrepancy between ψ and the kth order truncated the Fourier series decays faster than any inverse power of k. This behavior is referred to as spectral accuracy of the Fourier method [29].
The characteristic of spectral accuracy is also achievable for smooth non-periodic functions, on the condition that the trial functions are chosen properly. In non-periodic problems, the trial functions are orthogonal polynomial series. Similar to Fourier based methods, polynomial spectral methods start by the construction of an orthogonal basis for square integrable functions, in which the functions to approximate are expanded. Generally, these bases are generated using the Sturm-Liouville theorem [26,29,40]. The series coefficients demonstrate spectral accuracy if the trial functions are eigenfunctions of singular Sturm-Liouville problems. On the interval [−1, 1], polynomial solutions of singular Sturm-Liouville problems are Jacobi polynomials, e.g., Chebyshev and Legendre polynomials. A valuable property of the Legendre polynomials is that the weight function is unity, i.e., w(x) = 1; thus, simplifying the evaluation of integrals [40]. Furthermore, the unit weight function simplifies the integration by parts in Galerkin formulations of second order differential equations [44]. Therefore, Legendre polynomials are adopted, which are further detailed in the next paragraph.

Legendre Polynomials
Legendre polynomials satisfy the three-term recursion given by: where L N is the Legendre polynomial with N expansions and ξ is the local coordinate in the reference domain, i.e., the interval [−1, 1] [40]. The expansion of any function ψ ∈ L 2 [−1, 1] in terms of the Legendre polynomials, is given by: in which, The expansion coefficientsΦ k are dependent on the values of the function ψ in the physical domain. The expansion coefficients can seldom be calculated exactly. Alternatively, a finite number of approximate expansion coefficients, {Φ k } N k=0 , can be calculated using the values of the function at a finite number of selected points, ψ ξ j N j=0 , typically these points, {ξ j } N j=0 , are the nodes of a high-precision quadrature formula. This operation defines a discrete transformation between the set of values of the function at the quadrature nodes and the set of approximate expansion coefficients, i.e., a discrete transformation between ψ ξ j N j=0 and {Φ k } N k=0 . In case the corresponding quadrature formula is properly selected, the discrete transformation yields an interpolation of the values of the function at the quadrature nodes. Moreover, if the characteristics of spectral accuracy are maintained by the discrete transformation, the interpolant series can be used instead of the truncated series in order to approximate a function [29]. These aspects are further detailed in the paragraphs below and are essential to provide an approximation to the governing PDE.

Quadratures
A quadrature rule defines a formula for the approximation of the integral of a function, which is given by: where ψ is the function to integrate, {ξ j } N j=0 and {q j } N j=0 are the set of nodes and quadrature weights, respectively, and ε is the error. The Gauss quadrature rules are exact for the largest order polynomials [40]. Furthermore, the Gauss quadrature rules retain the characteristics of spectral accuracy [44]. As a result of the approximation by Jacobi polynomials, the Jacobi Gauss Quadrature has to be employed. In this case, the nodes are calculated as the zeros, i.e., roots, of the trial function [40]. However, the zeros of the Legendre polynomials, which are employed in this paper, are not located at the end-points. This feature is particularly disadvantageous for the implementation of boundary conditions, since the implementation of a boundary condition involves an interpolation step. Alternatively, for problems with Legendre weighted integrals, the Gauss-Lobatto rule can be applied, which includes the end-points. In this case, the nodes are calculated from the zeros of the Legendre-Gauss-Lobatto (LGL) polynomial, which results in: where L N is the LGL polynomial. The nodes are calculated through an iterative process. The zeros of the LGL polynomial are located at the end-points, regardless of the number of expansions. The corresponding quadrature weights are calculated according to [40]: Interpolation Spectral methods approximate a square integrable function as a finite expansion in orthogonal trial functions by the application of either one of two different techniques. The first technique is to truncate the infinite series. The second technique is to apply interpolation of a function at a set of nodes by Lagrange interpolating polynomials [40]. In this case, the characteristic of spectral accuracy is maintained in case the interpolation points are Gauss-type quadrature points corresponding to Jacobi polynomials [44]. Consequently, the approximation of a function is expressed as a Lagrange interpolant at the LGL nodes, which is given by: where ψ j is the value of the function to approximate at node ξ j and j (ξ) are the Lagrange interpolating polynomials, which are given by: where w j are the weights, which are given by [40]: The expansion coefficients correspond exactly to the values of the function at the nodes. Therefore, this approximation approach is referred to as a nodal basis, since each basis function, i.e., trial function, reproduces the value of the polynomial at a single specific node in the domain [29].

Differentiation
Similar to interpolation at a set of nodes, the derivative of a function at a set of nodes is approximated in the Lagrange form according to: where D i,j is the derivative matrix, which is given by: where the weights w j and w i are calculated according to (16). The diagonal elements are equal to zero, since the derivative of a constant is zero. In order to enforce this condition, the following condition is enforced on the diagonal: The diagonal elements from (19) are not exactly equal to zero; however, the effect of rounding errors is minimized [40]. The derivative matrix is an essential part for the approximation of a partial differential equation by the SEM, as is demonstrated in the next section.

Nodal Galerkin Approximation in the Reference Domain
Spectral methods are typically based on the Galerkin approximation, which shares a lot of similarities with finite element methods based on the Galerkin approximation. The main difference being the global nature of the approximation in the SEM, as opposed to the localized approximation in the FEM. The nodal representation of the Galerkin method is often preferred, since the expansion coefficient can rarely be computed exactly. In this section, the nodal Galerkin approximation for the governing PDE of the first class of regions, shown in (2), is presented in the reference domain, i.e., (ξ, η) ∈ [−1, 1] 2 . The nodal Galerkin approximation provides an approximate solution to the governing partial differential equation through the weak form. Furthermore, Green's first identity is applied to the weak form, such that the boundary and surface integrals are separated, which realizes the implementation of the boundary conditions through the manipulation of the boundary integral. Consequently, (2) is transformed into: where Ω is the physical domain, n is the unit vector normal to the boundary, and φ is the test function. The test function is a function from the same set of basis functions that are used to approximate the solution and satisfies the boundary conditions. In this case, the physical domain coincides with the computational domain, thus, the radial coordinate r is equivalent to the local coordinate ξ. In a nodal approximation, the test and trial function are expressed in the Lagrange form, which is given by, respectively: where φ i,j and ϕ n,m are the values of the test and trial function at the nodes, respectively, (ξ) and (η) are the Lagrange interpolating polynomials for the corresponding ξ and η axis, respectively, which are calculated according to (15). Equations (21) and (22) apply the same order, and thus, operate on the same set of grid points. The integrals in (20) are evaluated by the Gauss-Lobatto quadratures of (13), whereas the derivatives of the Lagrange interpolating polynomials are calculated according to (18) and (19). Under the assumption that Dirichlet boundary conditions are applied, the boundary integral vanishes, since the test function satisfies the boundary conditions, i.e., φ i,j = 0 [40]; thus, (20) is rewritten into: in which, Equation (23) represents the nodal Galerkin approximation on the reference element for the first class of regions in the benchmark system and describes a linear system for the interior nodal values. By the omission of the imposed current density and additional DoF terms, the description for the second class of regions is obtained. Similarly, by the omission of all electrically conductive terms, the description for the third class of regions is obtained.

Mapping Functions for Non-Squared Elements
The derivation of the nodal Galerkin approximation discussed in the previous section is valid if the physical domain coincides exactly with the computational domain, i.e., a [−1, 1] 2 square domain. However, in practice, electromagnetic and thermal problems consist of a wide variety of dimensions and shapes. In order to incorporate these geometrical variations, an algebraic map, i.e., r = (r, z) = M(ξ, η), between the coordinates in the physical domain (r, z) and the coordinates in the computational domain (ξ, η) is employed. The algebraic map transforms a point in the computational domain to a point in the physical domain. The most common method to perform a transformation between a quadrilateral in the physical domain and the reference square in the computational domain is called transfinite interpolation. The corresponding procedure is discussed in [40]. In this section, a brief overview of the essential equations is provided.
The mapping function, which performs transfinite interpolation, is given by: where Γ 1-4 are the four sides of an arbitrary quadrilateral. In an axisymmetric WPT system, such as the benchmark system, the physical domain typically consists of straight-lined rectangular elements. However, the mapping function is valid for any smooth curved quadrilateral. As a result of the mapping between the domains, the governing PDEs of the formulation are transformed as well. Typically, the transformation from the computational domain to the physical domain is known, whereas the inverse of the transformation often yields an intractable task. In order to circumvent this complication, the covariant and contravariant basis vectors, which are used to describe directions in the physical space, are introduced. The covariant basis vector a i is oriented tangentially to a coordinate line, whereas the contravariant basis vector a i is oriented normally to a coordinate line.
In an axisymmetric problem, the mapping between the domains is given by: where R and Z are the mapping components corresponding to the rand z-direction, respectively. The covariant basis vectors are defined as: where R ξ , Z ξ , R η , and Z η are the derivatives of the mapping components with respect to the ξ and η coordinates. The derivatives are given by: The contravariant basis vectors are given by: where J d is the determinant of the Jacobian, which is given by: Furthermore, the gradient of a function under mapping is given by: where ψ is an arbitrary function. The covariant flux components along the ξ-and η-axis are distinguishable (33), which are given by the following expressions, respectively, The contravariant flux components are given by: in which the determinant of the Jacobian is derived from (32). The weak form of the governing PDE in the reference domain, given by (20), consists of the dot product between the gradients of the test function and the solution. Under mapping, this term is derived from (33) and can be expressed in terms of the contravariant flux components, i.e., (36) and (37), which results in: Furthermore, the governing PDE given in (20) contains the boundary normals, which, under mapping, are given by: Consequently, from (33), (39) and (40), the dot product between the gradient and the boundary normals is derived, which in terms of the contravariant flux components results in: where F n are the fluxes normal to the boundary. Subsequently, the governing PDE of the nodal Galerkin approximation under mapping is derived by the substitution of (38) and (41) into (20). Furthermore, the integration on the physical domain is modified according to dΩ = J d dξdη, thus, the governing PDE yields: in which, whereΩ is the computational domain.

Nodal Galerkin Approximation in Matrix Form
In order to obtain the electromagnetic field solution with the SEM, the investigated domain is divided into multiple elements of matching material properties, as shown in Figure 2. Furthermore, the boundary conditions are applied to each element, such that a boundary-value problem is obtained, which is governed by (42). The boundary conditions are imposed through the boundary integral of (42), which is discussed in the next section. In order to solve the system of linear equations, the problem is expressed in its equivalent matrix form, which is given by: where A and C are the stiffness and mass matrix, respectively, K and W are the matrices required to include the χ Ã θ term, I k is a k × k identity matrix, k is the number of conductors, i.e., regions having both an imposed current and electrical conductivity, ϕ is the column vector containing the values of the electromagnetic field solution at the nodes, χ k is a column vector containing k values of the χ Ã θ term, y is the column vector containing the source terms, and 0 k is a column vector containing k zeros. The assembly of the various sub-matrices of (43) is discussed in this section.
In each element, a 2D grid of nodes is created in which the electromagnetic field solution is approximated. The grid consists of the roots of the LGL polynomial, which are calculated by forcing (12) to zero. The nodes along the ξ and η axes (in the computational domain) are stored in two separate vectors: Hence, a [N + 1] × [M + 1] grid is created in each element, in which the solution is approximated. Similarly, the quadratures, calculated from (13), are stored in vectors as well, from which a quadrature matrix for each element is created, given by: where q ξ and q η are the column vectors of the quadratures corresponding to the ξ-and η-axis, respectively. Furthermore, the derivative matrix is calculated according to (18) and (19). In order to find the partial derivative of a 2D function with respect to ξ, each row is multiplied with the derivative matrix. The partial derivative with respect to η is obtained by the multiplication of the rows with the derivative matrix. Therefore, the 2D derivative matrices are constructed through the Kronecker multiplication between the identity and derivative matrix, thus the partial derivatives with respect to ξ and η are given by, respectively [27,30], where I and D are the identity and derivative matrices, respectively. Subsequently, the contravariant flux components and the determinant of the Jacobian, i.e., (36), (37) and (32), respectively, are expressed in matrix form, which results in: where vec(. . .) is a vectorization operation, which returns a column vector of the input matrix, diag(. . .) returns a diagonal matrix where the elements of the input vector are on the diagonal of the output matrix, and •, •2 , and •−1 are the Hadamard product, power, and inverse, respectively, i.e., an element-wise matrix operation, ν is the reluctivity matrix, R ξ and R η are the derivative matrices of the mapping in the r-direction with respect to the ξ-and η-coordinates, respectively. Similarly, Z ξ and Z η correspond to the z-direction. By the substitution of (46)-(51) into (42), the sub-matrices of (43) are derived, which results in: in which, where R is a matrix containing the radial coordinates in the physical domain, σ is the conductivity matrix, sum(. . .) calculates the sum of the elements in the input vector, the subscript * indicates a sub-matrix or vector for each element, i.e., the matrices A, and C, as well as the vector y consist of the sub-matrices and sub-vector for each element. Similarly, the subscript k refers to the sub-vector for each conductor, from which the matrices K and W are constructed. By the substitution of (52)-(57) into (43), and the implementation of the boundary conditions, the linear system of equations can be solved according to: In case the problem is sufficiently small, a direct solver can be applied. Alternatively, iterative solving methods, e.g., the generalized minimum residual method or the quasiminimal residual method, can be applied in order to solve the linear system of equations.

Boundary Conditions
In order to obtain the electromagnetic field solution at the nodes, boundary conditions have to be considered. In an electromagnetic problem, three different types of boundary conditions are identified, i.e., continuous, Dirichlet, and Neumann boundary conditions. In this section, the application and implementation of the aforementioned boundary conditions are discussed.

Continuous Boundary Condition
The continuous boundary condition is an important asset for the purpose of establishing a multi-domain approach, since elements couple with each other in case the solution and test functions are continuous at the nodes shared between elements [40]. Consequently, the solution at the shared nodes has to be equal, as well as the values of the test functions. In the equivalent matrix representation of (43), this continuity between elements is implemented in the global matrix, i.e., E, by the summation of the matrix entries that have shared nodes. Therefore, the shared boundary nodes of two neighboring elements are represented by a single row or column, depending on the orientation of the boundary. Furthermore, the continuity of the first normal derivative of the solution is established by equating the boundary integral of (42) on the shared edges of neighboring elements.

Dirichlet Boundary Condition
The Dirichlet boundary condition defines the solution along a boundary. In order to include the Dirichlet boundary condition in the equivalent matrix representation of (43), the entries of the source vector y i , which correspond to the boundary nodes where the condition is valid, are modified according to: whereÃ θ,0 is the imposed solution. Subsequently, the entries in the global matrix E corresponding to the boundary nodes are replaced by a diagonal filled with ones.

Neumann Boundary Condition
The Neumann boundary condition defines the first derivative of the potential normal to the boundary. In order to implement the Neumann boundary condition, the boundary integral of (42) is used. Consequently, the contravariant flux components from (41) are replaced by the imposed value, which in the equivalent matrix representation yields: where F n,0 , F 1,0 , and F 2,0 are the imposed values. Therefore, the corresponding entries of the source vector y i are given by: where r is a column vector containing the radial coordinates, q is the quadrature column vector corresponding to either the ξ or η axis, depending on the orientation of the boundary condition, and the subscript i refers to the entries corresponding to the boundary nodes where the condition is imposed. As opposed to the Dirichlet boundary condition, the global matrix E remains unchanged. In electromagnetic problems, the Neumann boundary condition is typically used to mimic axes of symmetry, or to represent a boundary with a soft-magnetic material having an infinite relative permeability.

Post-Processing
Once the electromagnetic field solution at the nodes has been obtained, post-processing calculations are performed, such that quantities like the magnetic flux density are accessible. In this section, various post-processing calculations relevant to axisymmetric WPT systems, including the equivalent matrix form, are discussed.

Magnetic Flux Density
The magnetic flux density for an axisymmetric problem in terms of the modified magnetic vector potential is given by: By the inclusion of the appropriate sign convention and the radial coordinate of (62) into (33), the radial and axial components of the magnetic flux density are derived, which in the equivalent matrix form yields, respectively, where mat M×N (. . .) reshapes the input vector into a matrix of dimension M × N, which corresponds to (45) and (44), respectively. The post-processing calculation of the magnetic flux density is performed for each element.

Hysteresis Effect Losses
In soft-magnetic materials, the magnetic flux density can be used to approximate the iron losses. In a ferrite core, which is present in the benchmark system, these losses are primarily caused by the hysteresis effect. The corresponding losses are empirically determined according to Steinmetz's equation, which is given by: where c h , α, and β are empirical coefficients, f is the electrical frequency, and V is the volume of the corresponding region [45]. In the equivalent matrix form, (65) yields: where B m is a matrix containing the values of the magnitude of the magnetic flux density at the nodes. Similar to the calculation of the magnetic flux density, the hysteresis losses are evaluated per element. Therefore, in case multiple soft-magnetic regions are present, the individual results have to be summed in order to obtain the total losses within the domain.

Joule Effect Losses
In addition to the hysteresis effect, losses occur in an axisymmetric WPT system due to the Joule effect. In this case, the losses are calculated according to: where J θ is the total current density the azimuthal direction. In the equivalent matrix form, the current density in a source region is written as: where J θ,0 is given by (57). Consequently, in the first class of regions, the losses due to the Joule effect consist of the imposed and induced currents. In the second class of regions, the first two terms of (68) vanish and only the induced current causes losses, hence, the current density is given by: Subsequently, in the equivalent matrix form, Joule effect losses are calculated as: For multiple elements, the total losses due to the Joule effect are calculated as the sum of the individual contributions of each element.

Flux-Linkage
Besides the magnetic flux density and the loss components, another essential postprocessing quantity is the flux-linkage. The various inductances defining the equivalent circuit of a WPT system, are obtained from the flux-linkages. In a magnetically linear system, the flux-linkage is given by: where N t is the number of turns and A is the surface area of the corresponding region. Consequently, in the equivalent matrix form, the flux-linkage of a coil is obtained from: In case a coil consists of foil windings, each turn is modeled as a separate element, as shown in Figure 2; hence, the total flux-linkage is calculated as the sum of the individual contributions of each element.

Thermal Formulation
Apart from electromagnetic modeling, the SEM is also applicable for the thermal modeling of axisymmetric WPT system, such as the benchmark system. In steady-state conditions, the temperature distribution is governed by the heat conduction equation, which is given by where k is the thermal conductivity, T is the temperature, and q v is the volumetric heat source. Similar to Section 2.3, the governing PDE is modified according to the approach discussed in [41], such that the governing PDE is valid for an axisymmetric problem, whereas the Del operator is evaluated as if the coordinate-system were Cartesian. Therefore, (73) transforms into:

Nodal Galerkin Approximation
The weak form of (74) is obtained by applying the same approach as is discussed in Sections 2.3.2 and 2.3.3. Consequently, the weak form of (74) is given by: Similar to Section 2.3.4, the system of linear equations governed by (75) is expressed in the equivalent matrix form according to: where E is the global matrix, u is the column vector containing the values of the temperature at the nodes, and b is the column vector containing the source terms. The boundary integral of (75) is again used to apply the boundary conditions, the implementation of which for thermal problems is discussed in the next section. The global matrix and source vector of (76) have a sub-matrix and sub-vector for each element, respectively, which are given by: where S is a matrix containing the volumetric heat source values in the corresponding element. The matrices of the contravariant flux components, i.e., F 1 and F 2 , are given by (49) and (50), respectively, with the exception that the reluctivity matrix ν is replaced by the thermal conductivity matrix K th .

Boundary Conditions
The temperature at the nodes is obtained by solving the boundary-value problem governed by (75) in conjunction with the boundary conditions. In thermal problems, similar to Section 2.3.5, continuous, Dirichlet, and Neumann boundary conditions occur. Furthermore, the effects of convection and radiation can be included on a boundary through the Robin boundary condition, which is also present in the benchmark system. In this section, the application and implementation of boundary conditions occurring in thermal problems is discussed. The implementation of the continuous and Dirichlet boundary conditions is analogous to Section 2.3.5, and thus, are not repeated in this section.

Neumann Boundary Condition
Equivalent to electromagnetic problems, the Neumann boundary condition is implemented through the boundary integral of (75). Consequently, the corresponding entries of the source vector b i are given by: where f n is given by (60). Additionally, the global matrix E remains unchanged. An adiabatic boundary condition is imposed by forcing the first derivative of the temperature normal to the boundary to zero.

Robin Boundary Condition
In case both convection and radiation are present on a boundary, the Robin boundary condition is applied. In order to implement the condition, the corresponding entries of the source vector b i are modified according to: where h cv is the convection coefficient, T f and T 0 are the fluid and ambient temperature, respectively, ε is the emissivity coefficient, σ sb is the Stefan-Boltzmann constant, and u i is a column vector containing the temperature values at the boundary nodes where the condition is imposed. Consequently, in case radiation is present, an iterative or non-linear solver is required to solve the system of equations. In the context of this paper, an iterative solver is included, which, from an initial temperature rise, iterates until the temperature difference with respect to the previous iteration has converged to a value below a specified tolerance. Furthermore, in order to fully include convection on the boundary, the global matrix is summed with a diagonal matrix, which is given by:

Post-Processing
Once the temperature distribution in the domain has been obtained, a post-processing calculation is performed, such that the heat flux density is accessible, which, for an axisymmetric problem, is given by: Similar to the calculation of the magnetic flux density, the thermal flux density is obtained by including the sign convention and material properties of (82) into (33). Consequently, the radial and axial components of the thermal flux density are given by, respectively, Furthermore, the average temperature within an element is an important post-processing quantity, e.g., for the estimation of a uniform value for a temperature-dependent material property, which is given by: where A is the surface area of the corresponding element. In the equivalent matrix form, the average temperature in an element is calculated according to:

Detailed Overview of the SEM Implementation
The formulation of the SEM presented in Sections 2.3 and 2.4 has been implemented in MATLAB. A flowchart of the main processes executed in the software for an electromagnetic problem is shown in Figure 4, which is a detailed version of the flowchart shown in Figure 3. The processes shown in the figure are described as follows: 1.
The nodes and quadratures are generated for each element in the domain according to (11) and (13), respectively. Furthermore, the derivative matrices are constructed for each element according to (18) and (19).

2.
The mapping function and the corresponding derivatives are evaluated for each element according to (25) and (29), respectively. 3.
The sub-matrices and sub-vectors are constructed for each element in the domain according to (52)-(57). To that end, the quadrature matrix, partial derivative matrices, contravariant flux components, and the determinant of the Jacobian are obtained according to (46)-(51).

4.
The sub-matrices and sub-vector are assembled into the global matrix, E, and source vector, b, of (43). Furthermore, the boundary conditions are imposed according to the approach discussed in Section 2.3.5. 5.
The system of equations is solved to obtain the approximate solution at the nodes.
In the context of this paper, a direct solver is applied. Alternatively, for problems having a high number of non-zero elements, an iterative solver can be applied. In case the problem is non-linear, a non-linear solver, e.g., the Newton-Raphson or JFNK methods, has to be included to solve the non-linear state. 6.
After the solution has been obtained, the post-processing computations of Section 2.3.6 are executed. However, depending on the problem, additional quantities, e.g., the electromagnetic torque, can be added.
Generate nodes, quadratures, and derivative matrices Evaluate mapping function and obtain its derivative Generate sub-matrices and sub-vectors Assemble system of equations, including boundary conditions Solve system of equations to obtain approximate soltuion Obtain requested post-processing quantities 1.

5.
6. For the purpose of solving a thermal problem, the main processes are similar, unless radiation is present. In that case, an iterative solving approach is included, which iterates the fourth and fifth processes until the difference between the current and previous iteration is less than a specified tolerance.

Verification
In this section, the method for verification of the multi-physical modeling framework based on the SEM is introduced. The verification is carried out on the benchmark system from Section 2.1. The magnetic vector potential and temperature distribution in the domain are compared to a reference solution, which has been obtained with the FEM from a commercial software package. Furthermore, the discrepancies in the approximation of various key post-processing quantities, e.g., electromagnetic losses and the average coil temperature, are calculated with respect to the FEM. Additionally, the computational effort is assessed for both the SEM and FEM in terms of the number of DoF, number of non-zero elements, and the sparsity of the global matrices.

SEM Model
In order to obtain the solution at the nodes using the SEM, i.e., magnetic vector potential and temperature distribution, the domain is discretized into rectangular elements of matching material properties, as shown in Figure 2. The polynomial degree of the approximation is discretized into two different parameters. The first parameter is assigned along the axial direction of the foil windings and the adjacent elements, since the magnetic vector potential is expected to have a high gradient in these elements. The second parameter is assigned to the remainder of the model. In the magnetic model, the values are set to fifteen and ten, respectively, such that an accurate approximation of the magnetic vector potential is obtained. In the thermal model, the value is reduced to seven for both parameters, as a result of skin-depth effects being absent. In order to assess the computational effort, the number of DoF, number of non-zero elements, and the sparsity of the global matrices are stored for both the electromagnetic and thermal model.
Once the electromagnetic field solution at the LGL nodes has been obtained with the SEM, the post-processing computations are executed. The hysteresis losses in the magnetic cores, Joule losses in the electrically conductive materials, and the flux-linkage in the coils are calculated according to (66), (70) and (72), respectively. The individual loss component in each element is divided by the corresponding volume, such that the volumetric heat sources are obtained, which serve as the input for the thermal model. In the thermal model, the average coil temperature is post-processed from the temperature distribution at the nodes according to (86).

FEM Model
The FEM model is divided into rectangular elements equivalently to the SEM model, where in each element, a rectangular mesh is applied. At least two mesh layers per skindepth layer are ensured in electrically conductive materials, such that the second order elements are able to provide a sufficiently accurate approximation of the field variable. The mesh elements in the remainder of the model are tailored such that well-proportioned mesh elements are achieved in the entire domain, thus yielding a high mesh quality [46]. The same mesh is used for both the magnetic and thermal FEM models. Equivalent to the SEM magnetic model, the various losses and flux-linkages are obtained from the solution by performing post-processing operations. In the thermal model, the volumetric heat sources, which were obtained with SEM, are used as the input, such that the temperature distributions are compared using the same sources. The average coil temperature is postprocessed from the temperature distribution. The FEM model is solved using commercial software, i.e., Altair Flux [46]. Furthermore, the number of DoF, number of non-zero elements, and sparsity are stored for both models.

Method of Comparison
The rectangular FEM mesh in each element serves as the grid of nodes in which the solutions are compared. Firstly, the SEM solution at the LGL nodes is interpolated to the FEM grid according to (14). At each node in the FEM grid, the difference between the solutions is calculated as: where x is either the magnetic vector potential or the temperature, and the subscripts s and f refer to the SEM and FEM solution, respectively. Secondly, the absolute value of the relative discrepancy in the post-processing quantities, i.e., flux-linkages, Joule, and hysteresis effect losses, and average coil temperature, is calculated according to: where f s and f f are the post-processing quantities obtained with either the SEM or the FEM, respectively. The discrepancy in the losses is calculated per region, i.e., the primary and secondary coil, the stainless-steel shaft, the aluminum housing, and the primary and secondary ferrite core. Thirdly, the computational efforts of the SEM and FEM benchmark models are assessed in terms of the number of DoF, number of non-zero elements, and the sparsity of the global matrices. Additionally, the sparsity patterns of the SEM are shown. In the case of the FEM, these details are not accessible, since commercial software is applied.

Results
The magnetic vector potential obtained with the SEM and the difference with respect to the FEM are shown in Figure 5. The modulus, real, and imaginary part of the magnetic vector potential in the investigated domain are shown in the figure. In all three cases, an accurate match between the methods is obtained. Negligible differences between the methods are situated on the edges of the elements, which are caused by Runge's phenomenon, i.e., oscillations that occur at the edges of the interval when using higher order polynomial interpolation over a set of equally spaced interpolation points [28].
The magnitude of the magnetic flux density obtained with the SEM, which is postprocessed from the magnetic vector potential, including the difference with respect to the FEM is shown in Figure 6. In this case, an accurate match is obtained within the elements; however, relatively large differences are observed on the edges, in particular at the corner points inside the secondary transformer core. These differences are, apart from Runge's phenomenon, caused by singularities that appear due to the abruptly changing line that describes the corner geometry. As a result, the derivative of the magnetic vector potential is discontinuous at the corner points. The singularities at the corner points can be averted by applying a rounded corner [27].
The resulting values of the post-processing quantities for both methods, as well as the absolute value of the relative discrepancy, are shown in Table 2. As a result of the accurate approximation of the magnetic vector potential, the flux-linkages and Joule losses are accurately approximated as well, a maximum discrepancy of 0.078% with respect to the FEM is observed. The impact of the corner singularities on the losses due to the hysteresis effect is limited, a discrepancy of less than 0.15% with respect to the FEM is observed.
The temperature distribution obtained with the SEM and the difference with respect to the FEM is shown in Figure 7. Similar to Figure 5, an accurate match between the methods is obtained and negligible differences due to Runge's phenomenon are observed on the edges. Consequently, the average primary and secondary coil temperatures, the results of which are included in Table 2, are accurately calculated, a maximum discrepancy of 0.081% with respect to the FEM is observed. The computational effort indices for both methods and models are given in Table 3. In the case of the electromagnetic model, the SEM reduces the number of DoF and non-zero elements by 88.3% and 83.0%, respectively. The sparsity is slightly reduced by the SEM model, i.e., 0.121%. In the case of the thermal model, the application of the SEM, with respect to the FEM, reduces the number of DoF and non-zero elements by 94.6% and 94.2%, respectively. In the case of the thermal model, the reductions are more significant as a result of the SEM applying a decreased order for the approximation, whereas the FEM mesh density is unchanged. The sparsity is again slightly reduced by the SEM model, i.e., 0.197%. Additionally, the sparsity patterns of the SEM global matrices are shown in Figure 8. For this type of matrix structure, i.e., sparse, square, unsymmetric, and with non-zero entries that are not confined in a narrow band, the MATLAB-based SEM framework applies the UMFPACK direct solver [47].

Conclusions
In this paper, a multi-physical modeling framework based on the SEM for axisymmetric WPT systems has been discussed. The framework consists of an electromagnetic and a thermal model, of which the former incorporates eddy currents in both source and non-source regions. The thermal formulation is based on the steady-state heat conduction equation, where heat transfer by means of convection and radiation can be applied on the boundaries of the domain.
The SEM introduced in this paper is a multi-domain method, which is based on the nodal Galerkin approximation with numerical integration. The nodal representation of the Galerkin method has been selected, since the expansion coefficients of a trial function can rarely be computed exactly. In order to attain the property of spectral accuracy, the trial functions are the Legendre polynomials in the [−1, 1] 2 domain; thus, a mapping between the physical and computational domain is employed. Furthermore, as a result of the Legendre polynomials being selected, the weight function is unity throughout the domain. The integration is performed by the application of Gauss quadrature rules, which retain the property of spectral accuracy and are exact for the largest order polynomials. The Gauss quadrature rules approximate an integral at a set of nodes, which are the zeros of the trial function. In order to include the end-points to the set of nodes, the Gauss-Lobatto quadrature rule is applied, which calculates the nodes from the zeros of the LGL polynomial. In the nodal form, the trial function is expressed in terms of Lagrange interpolating polynomials, where the interpolation points are the Gauss-Lobatto quadrature nodes, such that spectral accuracy is maintained. In order to solve the linear system of equations, the matrix form of the nodal Galerkin approximation has been introduced.
The implementation of the multi-physical modeling framework has been verified by the FEM. A benchmark system, which is representative of an axisymmetric WPT system, has been modeled using both the SEM and FEM. The resulting solutions, i.e., magnetic vector potential and temperature distribution, have been compared to each other. An accurate match between the methods, with the exception of negligible differences caused by the interpolation of the SEM solution to an equally spaced grid, i.e., Runge's phenomenon, has been observed for both solutions. Consequently, the post-processing quantities, i.e., fluxlinkage, the losses due the Joule effect, and the average coil temperature, are calculated with a maximum discrepancy of 0.081% with respect to the FEM. A slightly larger discrepancy with respect to the FEM of 0.15% is observed in the losses due to the hysteresis effect in the ferrite transformer core, which is caused by corner singularities. Therefore, the results of the comparison have verified the modeling framework. Additionally, the computational effort in terms of the number of DoF, number of non-zero elements, and the sparsity has been assessed for both methods and models of the benchmark system. The results have demonstrated that, in the case of the electromagnetic model, the SEM reduces the number of DoF and non-zero elements by 88.3% and 83.0% with respect to the FEM, respectively. In the case of the thermal model, reductions by 94.6% and 94.2% are observed. In both cases, the sparsity is slightly reduced, i.e., 0.121% and 0.197% for the electromagnetic-and thermal model, respectively. However, as a result of the significant reductions in the number of DoF and non-zero elements, the multi-physical modeling framework based on the SEM is expected to significantly reduce the computation time for this type of problem.

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

Abbreviations
The following abbreviations are used in this manuscript:

2D
Two-dimensional 3D Three-dimensional ac Alternating current DoF Degree of freedom FEM Finite element method JFNK Jacobian-free Newton-Krylor LGL Legendre-Gauss-Lobatto PDEs Partial differential equations SEM Spectral element method SIBCs Surface impedance boundary conditions WPT Wireless power transfer