Next Article in Journal
The Performance Analysis of an Integrated CO2 Refrigeration System with Multi-Ejectors Installed in a Supermarket
Next Article in Special Issue
Performance Estimation: CCL WPT Topologies with Helical Coils
Previous Article in Journal
Assessment of Environmental Loads in the Life Cycle of a Retail and Service Building
Previous Article in Special Issue
Generalized Circuit Model of Shielded Capacitive Power Transfer
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

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

Department of Electrical Engineering, Electromechanics and Power Electronics, Eindhoven University of Technology, 5600 MB Eindhoven, The Netherlands
*
Author to whom correspondence should be addressed.
Energies 2022, 15(9), 3145; https://doi.org/10.3390/en15093145
Submission received: 30 March 2022 / Revised: 16 April 2022 / Accepted: 22 April 2022 / Published: 25 April 2022
(This article belongs to the Special Issue Modelling of Wireless Power Transfer II)

Abstract

:
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 source- and 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.

1. 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 stationary- and 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 skin- and 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].
A wide variety of WPT systems for rotary applications have been designed and analyzed in the literature [3,10,11,12,13,14,15,16,17]. For the purpose of designing and analyzing such systems, a modeling framework, which consists of an electromagnetic and a thermal model, is imperative to predict the characteristics of the system. In order to significantly reduce the computational effort, the cylindrical geometry is often reduced to an axisymmetric problem, e.g., as done in (part of) the analysis presented in [3,10,13,14,16,17]. Generally, the electromagnetic model employs the finite element method (FEM), e.g., the works presented in [2,3,5,17]. Similarly, in order to analyze the temperature distribution of the system, the FEM is often applied, for instance [13,16,18].
A WPT system is often a subsystem of a unified whole. Consequently, the application of shielding techniques may be necessary, such that undesirable effects, caused by the penetration of the magnetic field into surrounding objects, are mitigated. For example, in [17], a WPT system is designed for a satellite application. In order to minimize interference caused by the air gap fringing-fluxes, an open-ended cylinder acting as a shield surrounds the cylindrical transformer. The incorporation of surrounding electrically conductive materials, e.g., as a result of shielding techniques, into the modeling framework is crucial, since the performance might degrade significantly [19]. However, for high-frequency applications in which skin-depths are orders of magnitude smaller compared to the overall model dimensions, the locally different requirements on the mesh density may cause meshing problems in the electromagnetic model [20,21]. Therefore, electrically conductive materials present in the direct vicinity of a WPT system are often replaced by surface impedance boundary conditions (SIBCs). This type of model reduction introduces a (small) discrepancy, which is dependent on the relation of the radius to the skin-depth [22,23]. For this reason, SIBCs are not applicable to any arbitrary problem. Furthermore, SIBCs are unable to obtain the magnetic field distribution inside the electrically conductive region. In coil regions, analytical models, e.g., such as presented in [8,24], are often applied for the approximation of the coil ac resistance. However, such approximations might lead to significant underestimations of the ac resistance [25]. Consequently, the accuracy of the thermal model and efficiency calculation is compromised. Alternatively, for these types of problems, higher-order methods can be applied, such as the spectral element method (SEM), which, in comparison to the FEM, provides superior accuracy at the expense of geometrical domain flexibility [26,27].
The SEM is a numerical method, which can be applied in order to obtain approximate solutions to boundary-value problems. The SEM was originally applied to problems occurring in fluid dynamics and meteorology [28]. Over the past decades, the application of the method has spread to a wide variety of linear and non-linear engineering and mathematical problems [29], e.g., wave propagation, structural analysis, astrophysics, financial engineering, and more recently electromagnetic field modeling [27,30], including eddy current problems [16,20,31]. The formulation of the SEM for 2D Cartesian magneto-static and thermal problems is presented in [27], of which the former has been applied to a linear synchronous actuator in [32]. In [20], the SEM has been applied for the modeling of induced losses in a high-speed rotating cylinder. Alternatively, approaches that combine the SEM and FEM are observed in the literature, e.g., in [21] a coupled SEM-FEM approach has been applied to model an induction machine, where the former is used to model rotor eddy current losses and the latter is applied to the remainder of the domain. Another example is observed in [33], where a coupled SEM-FEM approach is applied to investigate magneto-convection for transformer cooling. In the context of high-frequency axisymmetric WPT systems, in [31], the SEM has been applied for the modeling of the losses that are induced in an electrically conductive and permeable cylinder, which surrounds the domain. The results have demonstrated that, compared to the FEM, the SEM obtains a higher accuracy per degree of freedom (DoF) and lower computation time. Similarly, in [16], the SEM has been applied in order to design such a system, whereas the FEM has been applied to analyze the temperature distribution. However, in none of the aforementioned works have eddy currents in coil regions, i.e., consisting of foil windings, been incorporated into the modeling technique. Moreover, a thermal model to complement the modeling framework 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.

2. Materials and Methods

2.1. 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 stainless-steel 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 boundary-value 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.

2.2. 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 Section 2.3 and Section 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 Section 2.3.4 and Section 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 Section 2.3.5 and Section 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 Section 2.3.6 and Section 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.

2.3. 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:
A ˜ θ ( r , z ) e θ = r A θ ( r , z ) e θ ,
where A θ and A ˜ θ 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:
c · ν r c A ˜ θ = J θ , 0 + j ω σ 1 r A ˜ θ χ A ˜ θ ,
J θ , 0 = σ I r Ω r 1 σ d r d z ,
χ A ˜ θ = Ω r 1 σ A ˜ θ d r d z r Ω r 1 σ d r d z ,
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, χ A ˜ θ 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:
c · ν r c A ˜ θ = 1 r j ω σ A ˜ θ .
Finally, in the third class of region, neither a source nor an electrical conductivity are assigned, thus (5) reduces to:
c · ν r c A ˜ θ = 0 .
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.

2.3.1. 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:
ψ x = k = Φ ^ k φ k x ,
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 w 2 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:
L N + 1 ( ξ ) = 2 N + 1 N + 1 ξ L N ( ξ ) N N + 1 L N 1 ( ξ ) , in which , L 0 = 1 and L 1 = ξ ,
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:
ψ ξ = k = 0 Φ ^ k φ k ξ , in which , Φ ^ k = k + 1 2 1 1 ψ ξ φ k ξ d ξ , φ k ξ = L k ξ .
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 } k = 0 N , can be calculated using the values of the function at a finite number of selected points, ψ ξ j j = 0 N , typically these points, { ξ j } j = 0 N , 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 j = 0 N and { Φ ^ k } k = 0 N . 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:
Q ψ = j = 0 N ψ ξ j q j = a b ψ ξ d ξ + ε ,
where ψ is the function to integrate, { ξ j } j = 0 N and { q j } j = 0 N 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:
ξ j = L N ξ = 0 , j = 0 , 1 , , N ,
L N ξ = L N + 1 ξ L N 1 ξ ,
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]:
q j = 2 N ( N + 1 ) ( L N ( ξ j ) ) 2 .

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:
ψ ξ = j = 0 N ψ j j ξ ,
where ψ j is the value of the function to approximate at node ξ j and j ξ are the Lagrange interpolating polynomials, which are given by:
j ξ = w j ξ ξ j j = 0 N w j ξ ξ j ,
where w j are the weights, which are given by [40]:
w j = 1 n = 0 n j N ξ j ξ n .
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:
ψ ξ i = j = 0 N ψ j j ξ i = j = 0 N D i , j ψ j , i = 0 , 1 , , N ,
where D i , j is the derivative matrix, which is given by:
D i , j = w j w i 1 ξ i ξ j , i j ,
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:
D i , i = j = 0 j i N 1 D i , j , i = 0 , 1 , N .
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.

2.3.2. 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:
Ω ν r c A ˜ θ · n ϕ ξ , η d l + Ω c ϕ ξ , η · ν r c A ˜ θ d ξ d η + j ω Ω σ 1 r A ˜ θ χ A ˜ θ ϕ ξ , η d ξ d η = Ω J θ , 0 ϕ ξ , η d ξ d η ,
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:
ϕ ξ , η = i , j = 0 N ϕ i , j i ξ j η ,
A ˜ θ ξ , η = n , m = 0 N φ n , m n ξ m η ,
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:
n = 0 N φ n , j q j ν r D n i ξ + m = 0 N φ i , m q i ν r D j m η + j ω σ r φ i , j q i q j j ω σ χ q i q j = q i q j J θ , 0 , i , j = 1 , 2 , , N 1 , in which , χ = r 1 σ φ i , j q i q j r n , m = 0 N r 1 σ q n q m , J θ , 0 = σ I r n , m = 0 N r 1 σ q n q m .
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.

2.3.3. 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:
M ξ , η = 1 2 1 ξ Γ 4 η + 1 + ξ Γ 2 η + 1 η Γ 1 ξ + 1 + η Γ 3 ξ 1 4 1 ξ 1 η Γ 1 1 + 1 + η Γ 3 1 + 1 + ξ 1 η Γ 1 1 + 1 + η Γ 3 1 ,
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:
M ξ , η = R ξ , η e r + Z ξ , η e z ,
where R and Z are the mapping components corresponding to the r- and z-direction, respectively. The covariant basis vectors are defined as:
a 1 = M ξ = R ξ e r + Z ξ e z ,
a 2 = M η = R η e r + Z η e z ,
a 3 = e θ ,
where R ξ , Z ξ , R η , and Z η are the derivatives of the mapping components with respect to the ξ and η coordinates. The derivatives are given by:
M ξ = 1 2 Γ 2 η Γ 4 η + 1 η Γ 1 ξ ξ + 1 + η Γ 3 ξ ξ 1 4 1 η Γ 1 1 Γ 1 1 + 1 + η Γ 3 1 Γ 3 1 , M η = 1 2 1 ξ Γ 4 η η + 1 + ξ Γ 2 η η + Γ 3 ξ Γ 1 ξ 1 4 1 ξ Γ 3 1 Γ 1 1 + 1 + ξ Γ 3 1 Γ 1 1 .
The contravariant basis vectors are given by:
J d a 1 = a 2 × a 3 = Z η e r R η e z ,
J d a 2 = a 3 × a 1 = Z ξ e r + R ξ e z ,
where J d is the determinant of the Jacobian, which is given by:
J d = a i · a j × a k = R ξ Z η Z ξ R η .
Furthermore, the gradient of a function under mapping is given by:
c ψ = 1 J d Z η ψ ξ Z ξ ψ η e r + R η ψ ξ + R ξ ψ η e z ,
where ψ is an arbitrary function. The covariant flux components along the ξ - and η -axis are distinguishable (33), which are given by the following expressions, respectively,
F 1 = 1 J d Z η ψ ξ Z ξ ψ η ,
F 2 = 1 J d R η ψ ξ + R ξ ψ η .
The contravariant flux components are given by:
F 1 = Z η 2 + R η 2 J d ψ ξ Z ξ Z η + R ξ R η J d ψ η ,
F 2 = Z ξ Z η + R ξ R η J d ψ ξ + Z ξ 2 + R ξ 2 J d ψ η ,
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:
c A ˜ θ · c ϕ = 1 J d F 1 ϕ ξ + F 2 ϕ η .
Furthermore, the governing PDE given in (20) contains the boundary normals, which, under mapping, are given by:
n 1 = | J d | J d Z η e r R η e z Z η 2 + R η 2 ,
n 2 = | J d | J d Z ξ e r + R ξ e z Z ξ 2 + R ξ 2 .
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:
c A ˜ θ · n i = F n = | J d | J d Z η 2 + R η 2 F 1 , for Γ 2 and Γ 4 , | J d | J d Z ξ 2 + R ξ 2 F 2 , for Γ 1 and Γ 3 ,
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:
Ω ^ ν r F n ϕ ξ , η d l + Ω ^ ν r F 1 ϕ ξ , η ξ + F 2 ϕ ξ , η η d ξ d η + j ω Ω ^ σ 1 r A ˜ θ χ A ˜ θ ϕ ξ , η J d d ξ d η = Ω ^ J θ , 0 ϕ ξ , η J d d ξ d η , in which , J θ , 0 = σ I r Ω ^ r 1 σ J d d ξ d η , χ A ˜ θ = Ω ^ r 1 σ A ˜ θ J d d ξ d η r Ω ^ r 1 σ J d d ξ d η ,
where Ω ^ is the computational domain.

2.3.4. 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:
E u = b , in which , E = A + C K W T I k , u = φ χ k , b = y 0 k ,
where A and C are the stiffness and mass matrix, respectively, K and W are the matrices required to include the χ A ˜ θ 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 χ A ˜ θ 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:
ξ = ξ 0 ξ 1 ξ N T ,
η = η 0 η 1 η M T .
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:
Q = q ξ q η T ,
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],
L D ξ = I   D ,
L D η = D   I ,
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:
F 1 = diag vec ν Z η 2 + R η 2 J d 1 L D ξ diag vec ν Z ξ Z η + R ξ R η J d 1 L D η ,
F 2 = diag vec ν Z ξ Z η + R ξ R η J d 1 L D ξ + diag vec ν Z ξ 2 + R ξ 2 J d 1 L D η ,
J d = R ξ Z η Z ξ R η ,
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:
A * = L D ξ T diag vec R 1 diag vec Q F 1 + L D η T diag vec R 1 diag vec Q F 2 ,
C * = diag vec σ R 1 J d Q j ω ,
k k = vec σ R 1 J d Q sum vec σ R 1 J d Q ,
w k = vec σ R 1 J d Q j ω
y * = vec J θ , 0 J d Q , in which ,
J θ , 0 = σ I R sum vec σ R 1 J d Q ,
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:
u = E 1 b .
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 quasi-minimal residual method, can be applied in order to solve the linear system of equations.

2.3.5. 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:
y i = A ˜ θ , 0 ,
where A ˜ θ , 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:
f n F n , 0 = vec sign J d Z η 2 + R η 2 F 1 , 0 , for Γ 2 and Γ 4 , vec sign J d Z ξ 2 + R ξ 2 F 2 , 0 , for Γ 1 and Γ 3 ,
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:
y i = r i 1 f n i F n , 0 q i + vec J θ , 0 J d Q i ,
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.

2.3.6. 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:
B = 1 r A ˜ θ z e r + 1 r A ˜ θ r e z .
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,
B r = R 1 J d 1 Z η mat M × N L D ξ φ Z ξ mat M × N L D η φ ,
B z = R 1 J d 1 R η mat M × N L D ξ φ + R ξ mat M × N L D η φ ,
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:
P c = V c h f α | | B | | β d V ,
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:
P c = sum vec c h f α B m β R Q J d 2 π ,
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:
P j = V | J θ | 2 σ d V ,
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:
J θ = J θ , 0 + j ω σ R 1 sum vec R 1 σ mat M × N φ J d Q sum vec σ R 1 J d Q j ω σ R 1 mat M × N φ ,
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:
J θ = j ω σ R 1 mat M × N φ .
Subsequently, in the equivalent matrix form, Joule effect losses are calculated as:
P j = sum vec σ 1 | J θ | 2 R J d Q 2 π .
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 post-processing 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:
Ψ = N t A V A θ d V ,
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:
Ψ = N t sum vec mat M × N φ J d Q 2 π sum vec J d Q .
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.

2.4. 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:
· k T = q v ,
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:
1 r c · r k c T = q v .

2.4.1. Nodal Galerkin Approximation

The weak form of (74) is obtained by applying the same approach as is discussed in Section 2.3.2 and Section 2.3.3. Consequently, the weak form of (74) is given by:
Ω ^ r k F n ϕ ξ , η d l + Ω ^ r k F 1 ϕ ξ , η ξ + F 2 ϕ ξ , η η d ξ d η = Ω ^ r q v ϕ ξ , η J d d ξ d η .
Similar to Section 2.3.4, the system of linear equations governed by (75) is expressed in the equivalent matrix form according to:
E u = b ,
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:
E * = L D ξ T diag vec R diag vec Q F 1 + L D η T diag vec R diag vec Q F 2 ,
b * = vec R S J d Q ,
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 .

2.4.2. 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:
b i = r i f n i F n , 0 q i + vec R S J d Q i ,
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:
b i = r i f n i h c v T f q i + vec R S J d Q i + r i f n i ε σ s b u i 4 T 0 4 q i ,
where h c v is the convection coefficient, T f and T 0 are the fluid and ambient temperature, respectively, ε is the emissivity coefficient, σ s b 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:
H cv = diag r i f n i h c v q i .

2.4.3. 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:
q = k T r e r + T z e z .
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,
Q q , r = K th J d 1 Z η mat M × N L D ξ u Z ξ mat M × N L D η u ,
Q q , z = K th J d 1 R η mat M × N L D ξ u + R ξ mat M × N L D η u .
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:
T a = 1 A A T d A ,
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:
T a = sum vec mat M × N u J d Q sum vec J d Q .

2.5. Detailed Overview of the SEM Implementation

The formulation of the SEM presented in Section 2.3 and Section 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.
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.

2.6. 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.

2.6.1. 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).

2.6.2. 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 skin-depth 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 post-processed 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.

2.6.3. 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:
Δ x = x s x f ,
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:
ϵ = | f s f f | f f ,
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.

3. 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 post-processed 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].

4. 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., flux-linkage, 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.

Author Contributions

Writing—original draft, K.B.; writing—review and editing, D.C.J.K. and E.A.L. All authors have read and agreed to the published version of the manuscript.

Funding

This research received no external funding.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Not applicable.

Conflicts of Interest

The authors declare no conflict of interest.

Abbreviations

The following abbreviations are used in this manuscript:
2DTwo-dimensional
3DThree-dimensional
acAlternating current
DoFDegree of freedom
FEMFinite element method
JFNKJacobian-free Newton–Krylor
LGLLegendre–Gauss–Lobatto
PDEsPartial differential equations
SEMSpectral element method
SIBCsSurface impedance boundary conditions
WPTWireless power transfer

References

  1. Zhang, C.; Lin, D.; Hui, S.Y.R. Ball-Joint Wireless Power Transfer Systems. IEEE Trans. Power Electron. 2018, 33, 65–72. [Google Scholar] [CrossRef]
  2. Han, H.; Mao, Z.; Zhu, Q.; Su, M.; Hu, A.P. A 3D Wireless Charging Cylinder With Stable Rotating Magnetic Field for Multi-Load Application. IEEE Access 2019, 7, 35981–35997. [Google Scholar] [CrossRef]
  3. Maier, M.; Parspour, N. Operation of an Electrical Excited Synchronous Machine by Contactless Energy Transfer to the Rotor. IEEE Trans. Ind. Appl. 2018, 54, 3217–3225. [Google Scholar] [CrossRef]
  4. Xia, Q.; Yan, L. Application of Wireless Power Transfer Technologies and Intermittent Energy Harvesting for Wireless Sensors in Rotating Machines. Wirel. Power Transf. 2016, 3, 93–104. [Google Scholar] [CrossRef]
  5. Chen, X.Y.; Jin, J.X.; Zheng, L.H.; Wu, Z.H. A Rotary-Type Contactless Power Transfer System Using HTS Primary. IEEE Trans. Appl. Supercond. 2016, 26, 5501406. [Google Scholar] [CrossRef]
  6. Qian, W.; Zhang, X.; Fu, Y.; Lu, J.; Bai, H. Applying Normally-Off GaN HEMTs for Coreless High-Frequency Wireless Chargers. CES Trans. Electr. Mach. Syst. 2017, 1, 418–427. [Google Scholar] [CrossRef]
  7. Covic, G.A.; Boys, J.T. Modern Trends in Inductive Power Transfer for Transportation Applications. IEEE J. Emerg. Sel. Top. Power Electron. 2013, 1, 28–41. [Google Scholar] [CrossRef]
  8. van den Bossche, A.; Valchev, V.C. Inductors and Transformers for Power Electronics; Taylor & Francis Group: New York, NY, USA, 2005. [Google Scholar]
  9. Kulkarni, S.V.; Khaparde, S.A. Transformer Engineering; CRC Press: Boca Raton, FL, USA, 2013. [Google Scholar]
  10. Trevisan, R.; Costanzo, A. A 1-kW Contactless Energy Transfer System Based on a Rotary Transformer for Sealing Rollers. IEEE Trans. Ind. Electron. 2014, 61, 6337–6345. [Google Scholar] [CrossRef]
  11. Ditze, S.; Endruschat, A.; Schriefer, T.; Rosskopf, A.; Heckel, T. Inductive Power Transfer System with a Rotary Transformer for Contactless Energy Transfer on Rotating Applications. In Proceedings of the IEEE International Symposium on Circuits and Systems (ISCAS), Montreal, QC, Canada, 22–25 May 2016; pp. 1622–1625. [Google Scholar] [CrossRef]
  12. Zhu, X.; Lin, B.; Liu, L.; Luan, Y. Power Transfer Performance and Cutting Force Effects of Contactless Energy Transfer System for Rotary Ultrasonic Grinding. IEEE Trans. Ind. Electron. 2016, 63, 2785–2795. [Google Scholar] [CrossRef]
  13. Bastiaens, K.; Krop, D.C.J.; Jumayev, S.; Lomonova, E.A. Optimal Design and Comparison of High-Frequency Resonant and Non-Resonant Rotary Transformers. Energies 2020, 13, 929. [Google Scholar] [CrossRef] [Green Version]
  14. Raminosoa, T.; Wiles, R.H.; Wilkins, J. Novel Rotary Transformer Topology with Improved Power Transfer Capability for High-Speed Applications. IEEE Trans. Ind. Appl. 2020, 56, 277–286. [Google Scholar] [CrossRef]
  15. Ruviaro, M.; Runcos, F.; Sadowski, N.; Borges, I.M. Analysis and Test Results of a Brushless Doubly Fed Induction Machine with Rotary Transformer. IEEE Trans. Ind. Electron. 2012, 59, 2670–2677. [Google Scholar] [CrossRef]
  16. Bastiaens, K.; Krop, D.C.J.; Curti, M.; Jumayev, S.; Lomonova, E.A. Coupled Design Approach for an Integrated GaN-Based Wireless Power Transfer Rotating System. Int. J. Appl. Electromagn. Mech. 2021, 64, S51–S64. [Google Scholar] [CrossRef]
  17. Wang, F.; Feng, T.; Chen, X. Design and Optimisation of a Wireless Power Transfer System for Satellite Application. IET Power Electron. 2019, 12, 2586–2598. [Google Scholar] [CrossRef]
  18. Godbehere, J.; Hopkins, A.; Yuan, X. Design and Thermal Analysis of a Rotating Transformer. In Proceedings of the IEEE International Electric Machines Drives Conference (IEMDC), San Diego, CA, USA, 12–15 May 2019; pp. 2144–2151. [Google Scholar] [CrossRef]
  19. Campi, T.; Cruciani, S.; Feliziani, M. Magnetic Shielding of Wireless Power Transfer Systems. In Proceedings of the International Symposium on Electromagnetic Compatibility, Tokyo, Japan, 12–16 May 2014; pp. 422–425. [Google Scholar]
  20. de Gersem, H. Spectral-Element Method for High-Speed Rotating Cylinders. COMPEL 2009, 28, 730–740. [Google Scholar] [CrossRef]
  21. De Gersem, H. Combined Spectral-Element, Finite-Element Discretization for Magnetic-Brake Simulation. IEEE Trans. Magn. 2010, 46, 3520–3523. [Google Scholar] [CrossRef]
  22. Jayasekera, K.A.S.N.; Ciric, I.R. Evaluation of Surface Impedance Models for Axisymmetric Eddy-Current Fields. IEEE Trans. Magn. 2007, 43, 1991–2000. [Google Scholar] [CrossRef]
  23. Wagner, B.; Renhart, W.; Magele, C. Error Evaluation of Surface Impedance Boundary Conditions with Magnetic Vector Potential Formulation on a Cylindrical Test Problem. IEEE Trans. Magn. 2008, 44, 734–737. [Google Scholar] [CrossRef]
  24. Stein, A.L.; Kyaw, P.A.; Sullivan, C.R. Wireless Power Transfer Utilizing a High-Q Self-Resonant Structure. IEEE Trans. Power Electron. 2019, 34, 6722–6735. [Google Scholar] [CrossRef]
  25. Pereira, A.; Sixdenier, F.; Raulet, M.A.; Lefebvre, B.; Burais, N. Comparison Between Numerical and Analytical Methods of AC Resistance Evaluation for Medium-Frequency Transformers: Validation on a Prototype and Thermal Impact Analysis. Can. J. Electr. Comput. Eng. 2017, 40, 101–109. [Google Scholar] [CrossRef]
  26. Shen, J.; Tang, T.; Want, L.L. Spectral Methods Algorithms, Analysis and Applications; Springer: Berlin/Heidelberg, Germany, 2011. [Google Scholar]
  27. Curti, M.; Van Beek, T.A.; Jansen, J.W.; Gysen, B.L.J.; Lomonova, E.A. General Formulation of the Magnetostatic Field and Temperature Distribution in Electrical Machines Using Spectral Element Analysis. IEEE Trans. Magn. 2018, 54, 1–9. [Google Scholar] [CrossRef]
  28. Trefethen, L.N. Spectral Methods in Matlab; Society for Industrial and Applied Mathematics: Philadelphia, PA, USA, 2000. [Google Scholar]
  29. Canuto, C.; Hussaini, M.Y.; Quarteroni, A.; Zang, T.A. Spectral Methods Fundamentals in Single Domains; Springer: Berlin/Heidelberg, Germany, 2006. [Google Scholar]
  30. Curti, M. A Spectral Element Based Framework for Efficient Computation of Electromechanical Problems. Ph.D. Thesis, Eindhoven University of Technology, Eindhoven, The Netherlands, 2019. [Google Scholar]
  31. Bastiaens, K.; Curti, M.; Krop, D.C.J.; Jumayev, S.; Lomonova, E.A. Spectral Element Method Modeling of Eddy Current Losses in High-Frequency Transformers. Math. Comput. Appl. 2019, 24, 28. [Google Scholar] [CrossRef] [Green Version]
  32. Curti, M.; Paulides, J.J.H.; Lomonova, E.A. Magnetic Modeling of a Linear Synchronous Machine With the Spectral Element Method. IEEE Trans. Magn. 2017, 53, 8112006. [Google Scholar] [CrossRef] [Green Version]
  33. Zanella, R.; Nore, C.; Bouillault, F.; Cappanera, L.; Tomas, I.; Mininger, X.; Guermond, J. Study of Magnetoconvection Impact on a Coil Cooling by Ferrofluid with a Spectral/Finite-Element Method. IEEE Trans. Magn. 2018, 54, 4600104. [Google Scholar] [CrossRef]
  34. Regenauer-Lieb, K.; Veveakis, M.; Poulet, T.; Wellmann, F.; Karrech, A.; Liu, J.; Hauser, J.; Schrank, C.; Gaede, O.; Fusseis, F.; et al. Multiscale Coupling and Multiphysics Approaches in Earth Sciences: Applications. J. Coupled Syst. Multiscale Dyn. 2013, 1, 281–323. [Google Scholar] [CrossRef] [Green Version]
  35. Zhang, H.; Guo, J.; Lu, J.; Li, F.; Xu, Y.; Downar, T.J. An Assessment of Coupling Algorithms in HTR Simulator TINTE. Nuclear Sci. Eng. 2018, 190, 287–309. [Google Scholar] [CrossRef]
  36. Price, D.; Radaideh, M.I.; Mui, T.; Katare, M.; Kozlowski, T. Multiphysics Modeling and Validation of Spent Fuel Isotopics Using Coupled Neutronics/Thermal-Hydraulics Simulations. J. Coupled Syst. Multiscale Dyn. 2020, 2020, 2764634. [Google Scholar] [CrossRef]
  37. Fernández, M.Á.; Moubachir, M. A Newton Method Using Exact Jacobians for Solving Fluid–Structure Coupling. Comput. Struct. 2005, 83, 127–142. [Google Scholar] [CrossRef] [Green Version]
  38. Ao, X.; Liu, J.; Xia, H.; Yang, Y. A Numerical Study on the Mesoscopic Characteristics of Ti-6Al-4V by Selective Laser Melting. Materials 2022, 15, 2850. [Google Scholar] [CrossRef]
  39. Yan, D.; Wang, W.; Li, R.; Jiang, S.; Lu, L.; Levtsev, A.; Chen, D. Multi-Physical and Electrochemical Coupling Model for the Protonic Ceramic Fuel Cells with H+/e/O2− Mixed Conducting Cathodes. Appl. Sci. 2022, 12, 3889. [Google Scholar] [CrossRef]
  40. Kopriva, D.A. Implementing Spectral Methods for Partial Differential Equations; Springer Science + Business Media B.V.: Dordrecht, The Netherlands, 2009. [Google Scholar]
  41. Winslow, A. Numerical Calculation of Static Magnetic Fields in an Irregular Triangular Mesh; UCR-7784; University of California: Livermore, CA, USA, 1964. [Google Scholar]
  42. Touzani, R.; Rappaz, J. Mathematical Models for Eddy Currents and Magnetostatics; Springer Science + Business Media B.V.: Dordrecht, The Netherlands, 2014. [Google Scholar]
  43. Bermúdez, A.; Gómez, D.; Salgado, P. Mathematical Models and Numerical Simulation in Electromagnetism; Springer International Publishing AG: Cham, Switzerland, 2014. [Google Scholar]
  44. van de Vosse, F.N.; Minev, P. Spectral Element Methods: Theory and Applications; EUT Report; Department of Mechanical Engineering, Eindhoven University of Technology: Eindhoven, The Netherlands, 1996. [Google Scholar]
  45. Bertotti, G. Hysteresis in Magnetism; Academic Press, Inc.: San Diego, CA, USA, 1998; ISBN 978-0-12-093270-2. [Google Scholar]
  46. Altair Engineering, Inc., Meylan, France. User Guide Flux 2019.1. 2018. Available online: https://www.altair.com/flux/ (accessed on 29 March 2022).
  47. Davis, T.A. UMFPACK User Guide. 2019. Available online: https://fossies.org/linux/SuiteSparse/UMFPACK/Doc/UMFPACK_UserGuide.pdf (accessed on 29 March 2022).
Figure 1. Overview of the benchmark system in 3D space.
Figure 1. Overview of the benchmark system in 3D space.
Energies 15 03145 g001
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.
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.
Energies 15 03145 g002
Figure 3. Flowchart of the six main processes in the SEM.
Figure 3. Flowchart of the six main processes in the SEM.
Energies 15 03145 g003
Figure 4. Flowchart of the main processes in the SEM implementation for an electromagnetic problem.
Figure 4. Flowchart of the main processes in the SEM implementation for an electromagnetic problem.
Energies 15 03145 g004
Figure 5. Magnetic vector potential in the investigated domain obtained with the SEM: (a) modulus, (b) real, and (c) imaginary part. The difference with respect to the FEM solution is shown in (df), respectively.
Figure 5. Magnetic vector potential in the investigated domain obtained with the SEM: (a) modulus, (b) real, and (c) imaginary part. The difference with respect to the FEM solution is shown in (df), respectively.
Energies 15 03145 g005
Figure 6. Magnitude of the magnetic flux density (a) obtained with the SEM. The difference with respect to the FEM solution is shown in (b).
Figure 6. Magnitude of the magnetic flux density (a) obtained with the SEM. The difference with respect to the FEM solution is shown in (b).
Energies 15 03145 g006
Figure 7. Temperature distribution (a) obtained with the SEM. The difference with respect to the FEM solution is shown in (b).
Figure 7. Temperature distribution (a) obtained with the SEM. The difference with respect to the FEM solution is shown in (b).
Energies 15 03145 g007
Figure 8. Sparsity patterns of the global matrices in the (a) electromagnetic and (b) thermal SEM model of the benchmark system.
Figure 8. Sparsity patterns of the global matrices in the (a) electromagnetic and (b) thermal SEM model of the benchmark system.
Energies 15 03145 g008
Table 1. Geometrical parameters and physical quantities of the benchmark.
Table 1. Geometrical parameters and physical quantities of the benchmark.
Geometrical Parameters
ParameterSymbolValueUnit
Offset radius of the geometry R x 1.00pm
Outer radius of the shaft R s 1.50mm
Inner radius of the core R 1 r 1.61mm
Inner radius of the coil area R 2 r 2.62mm
Inner radius of the air gap R 3 r 4.13mm
Outer radius of the coil area R 4 r 6.21mm
Outer radius of the core R 5 r 7.30mm
Inner radius of the housing R i 7.50mm
Radial thickness of the housing r h 1.00mm
Inner radius of the 1st turn r 1 s 2.92mm
Inner radius of the 2nd turn r 2 s 3.68mm
Inner radius of the 3rd turn r 3 p 4.95mm
Inner radius of the 4th turn r 4 p 5.74mm
Radial thickness of the air gap r a g 5.00 × 10 1 mm
Base height of the housing h e 5.00 × 10 1 mm
Height of the surrounding air h b 5.79mm
Base height of the core h c 1.02mm
Total height of the coil area h w , t 4.38mm
Height of the coil h w 4.18mm
Thickness of the foil winding d w 1.00 × 10 1 mm
Physical Quantities
ParameterSymbolValueUnit
Electrical frequencyf1.00MHz
Electrical conductivity of the housing σ h 3.77 × 10 7 S · m 1
Thermal conductivity of the housing k h 237W · m 1 · K 1
Electrical conductivity of the shaft σ s 1.45 × 10 6 S · m 1
Thermal conductivity of the shaft k s 16.3W · m 1 · K 1
Relative permeability of the core μ r , c 2.30 × 10 3 -
Thermal conductivity of the core k c 4.25W · m 1 · K 1
Electrical conductivity of the foil σ f 5.81 × 10 7 S · m 1
Thermal conductivity of the foil k f 385W · m 1 · K 1
Thermal conductivity of the insulator k w 1.00W · m 1 · K 1
Thermal conductivity of air k a 2.57 × 10 2 W · m 1 · K 1
Ambient temperature T 0 293K
Convection coefficient h c v 7.00W · m 2 · K 1
Emissivity coefficient ε 4.00 × 10 1 -
Imposed current (primary side) I p 4.950.00 A
Imposed current (secondary side) I s 4.9590.0 A
Hysteresis-loss coefficient c h 10.6W · s α · T β · m 3
Hysteresis-loss coefficient α 1.30-
Hysteresis-loss coefficient β 2.70-
Table 2. Post-processing quantities obtained with both the FEM and the SEM, including the absolute value of the relative discrepancy.
Table 2. Post-processing quantities obtained with both the FEM and the SEM, including the absolute value of the relative discrepancy.
QuantitySymbolFEMSEMUnitDiscrepancy [%]
Flux-linkage primary coil Ψ p 1.2711.270 μ Wb-t5.540 × 10 2
Flux-linkage secondary coil Ψ s 1.2481.248 μ Wb-t5.616 × 10 2
Joule losses primary coil P j , p 189.0189.0mW2.159 × 10 4
Joule losses secondary coil P j , s 156.3156.3mW1.160 × 10 2
Hysteresis losses primary core P c , p 7.0727.063mW1.387 × 10 1
Hysteresis losses secondary core P c , s 37.8737.82mW1.494 × 10 1
Eddy current losses shaft P j , s 55.4155.39mW3.216 × 10 2
Eddy current losses housing P j , h 14.8714.86mW7.786 × 10 2
Temperature primary coil T a , p 61.6461.63 C2.240 × 10 2
Temperature secondary coil T a , s 68.2268.17 C8.050 × 10 2
Table 3. Computational effort of both the electromagnetic and thermal model of the benchmark system for both the SEM and FEM.
Table 3. Computational effort of both the electromagnetic and thermal model of the benchmark system for both the SEM and FEM.
Electromagnetic
QuantityUnitSEMFEM
Number of DoF-1.738 × 10 4 1.488 × 10 5
Number of non-zero elements-3.999 × 10 5 2.347 × 10 6
Sparsity%99.8799.99
Thermal
QuantityUnitSEMFEM
Number of DoF-8.128 × 10 3 1.507 × 10 5
Number of non-zero elements-1.367 × 10 5 2.350 × 10 6
Sparsity%99.7999.99
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Bastiaens, K.; Krop, D.C.J.; Lomonova, E.A. Spectral Element-Based Multi-Physical Modeling Framework for Axisymmetric Wireless Power Transfer Systems. Energies 2022, 15, 3145. https://doi.org/10.3390/en15093145

AMA Style

Bastiaens K, Krop DCJ, Lomonova EA. Spectral Element-Based Multi-Physical Modeling Framework for Axisymmetric Wireless Power Transfer Systems. Energies. 2022; 15(9):3145. https://doi.org/10.3390/en15093145

Chicago/Turabian Style

Bastiaens, Koen, Dave C. J. Krop, and Elena A. Lomonova. 2022. "Spectral Element-Based Multi-Physical Modeling Framework for Axisymmetric Wireless Power Transfer Systems" Energies 15, no. 9: 3145. https://doi.org/10.3390/en15093145

APA Style

Bastiaens, K., Krop, D. C. J., & Lomonova, E. A. (2022). Spectral Element-Based Multi-Physical Modeling Framework for Axisymmetric Wireless Power Transfer Systems. Energies, 15(9), 3145. https://doi.org/10.3390/en15093145

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop