*Article* **Memory Efficient Method for Electromagnetic Multi-Region Models Using Scattering Matrices**

#### **C. H. H. M. Custers \*, J. W. Jansen and E. A. Lomonova**

Department of Electrical Engineering, Electromechanics and Power Electronics, Eindhoven University of Technology, 5612 AZ Eindhoven, The Netherlands; j.w.jansen@tue.nl (J.W.J.); e.lomonova@tue.nl (E.A.L.)

**\*** Correspondence: c.h.h.m.custers@tue.nl

Received: 25 October 2018; Accepted: 7 November 2018; Published: 9 November 2018

**Abstract:** This paper describes the scattering matrix approach to obtain the solution to electromagnetic field quantities in harmonic multi-layer models. Using this approach, the boundary conditions are solved in such way that the maximum size of any matrix used during the computations is independent of the number of regions defined in the problem. As a result, the method is more memory efficient than classical methods used to solve the boundary conditions. Because electromagnetic sources can be located inside the regions of a configuration, the scattering matrix formulation is developed to incorporate these sources into the solving process. The method is applied to a 3D electromagnetic configuration for verification.

**Keywords:** scattering matrix; Fourier analysis; permanent magnet machines; analytical modeling

#### **1. Introduction**

The design and optimization of electrical machines and other electromagnetic applications necessitate accurate models of electromagnetic fields to precisely predict their performance. As these systems and their structures get more complex, it is becoming more difficult to accurately model, e.g., eddy current effects [1,2] or hysteresis effects [3] in both the spatial and time domain. To model electromagnetic field distributions and the related forces and power losses, the finite element method (FEM) is often used because of its ability to produce accurate results when it is correctly utilized. However, the method can be demanding in terms of memory and relatively slow in terms of computation time. Therefore, semi-analytical models have been proposed over the years for increasingly complex structures in both 2D and 3D. One of the semi-analytical models is the harmonic modeling technique [4–7], which uses a Fourier basis to describe the solutions to the electromagnetic field quantities. By including position dependent material properties, such as the electrical conductivity [8] and the permeability [9–11], into the solutions, the spatial distribution of electromagnetic fields can also be accurately modeled for complex devices.

In many electromagnetic configurations, accurate results are obtained using a relatively low number of harmonics. However, for more complex structures, with detailed features compared to the device size, the number of harmonics has to be increased to retain accuracy. Computationally, this leads to an increase in the required memory. In addition, as an electromagnetic configuration has to be divided into regions (or layers), the number of regions influences the required memory. As a result, the number of harmonics that can be considered is limited when a relatively large number of regions has to be incorporated. Consequently, the advantage in terms of memory of the harmonic model in comparison to FEM is reducing. The scattering matrix approach [12–16] reformulates a multiple region electromagnetic problem in such a way that it can be expressed as a single region with inand out-going fields. Consequently, the dependency of the computational memory on the number of regions is removed. All matrices that are manipulated with this method have a maximum size of *L* × *L*,

where *L* is the total number of harmonics. In [12], where the scattering matrices are used to solve a high-frequent electromagnetic problem, the source is an incident wave on one side of the problem. However, to apply the scattering formulation in low-frequent problems in which sources (currents or magnets) are present in regions, the scattering formulation has to be extended to incorporate the related source terms.

In this paper, the scattering matrix approach for harmonic models of electromagnetic configurations is described. Using this approach, the memory required to obtain the solutions of the model is significantly reduced. The method presented in [12] is extended to take electromagnetic sources inside regions into account. The method is verified by application to a 3D electromagnetic configuration presented in [10], where the permeability varies in both periodic directions as a function of position.

#### **2. Model Assumptions**

The mathematical models of electromagnetic configurations that can be solved with the method described in this paper are subject to a number of assumptions. The considered models are set up in the Cartesian coordinate system, where, for 2D configurations in one direction, periodicity is assumed. For 3D configurations, two directions are assumed periodic. Because of the periodicity, the solution in the periodic directions is expressed by Fourier series. Based upon the electromagnetic material properties, such as the permeability, permittivity and conductivity, and the electromagnetic field sources, the model is divided into regions in the non-periodic direction. In the non-periodic direction, the properties of a region may not vary, whereas, in the periodic direction, they are allowed to change as a function of position. Lastly, any configuration needs to have an air region at both sides, extending to plus and minus infinity in the non-periodic direction, for the applied method.

#### **3. Scattering Matrix Formulation**

For all problems at hand, the Fourier coefficients of describing the electromagnetic fields have the same form. For the remainder of the paper, the non-periodic direction is the *z*-direction. The *z*-dependent solution describing the coefficients **u** and **v** of particular field components of a region *i* is given by

$$\mathbf{u}\_{i} = \mathbf{Q}\_{\mathbf{u}i} \left( e^{\lambda\_{i}z} \mathbf{c}\_{i}^{+} + e^{-\lambda\_{i}z} \mathbf{c}\_{i}^{-} \right) + \mathbf{p}\_{u,i'} \tag{1}$$

$$\mathbf{v}\_{i} = \mathbf{Q}\_{\mathbf{v}\_{i}} \left( e^{\lambda\_{i} z} \mathbf{c}\_{i}^{+} - e^{-\lambda\_{i} z} \mathbf{c}\_{i}^{-} \right) + \mathbf{p}\_{v, i\prime} \tag{2}$$

where λ is a vector with eigenvalues and **Qu** and **Qv** are eigenvector matrices. Combined, λ and **Q** describe the propagation of field components inside a region. The particular solutions described by **p***<sup>u</sup>* and **p***<sup>v</sup>* relate source terms, such as coil currents and magnetization, which are known in advance, to the solutions **u** and **v**. All these variables are determined by the properties of a region that is under consideration. The way to obtain these variables for a variety of problems has been described in e.g., [8–10]. The unknown coefficients **c**<sup>+</sup> *<sup>i</sup>* and **c**<sup>−</sup> *<sup>i</sup>* in (1) and (2) are determined by applying continuity boundary conditions between regions.

The continuity boundary conditions force both **u** and **v** to be continuous on an interface between adjacent regions. Often in harmonic models used for electrical machines, all boundary conditions are gathered in a large matrix which forms a system of linear equations,

where *h* specifies the height of a boundary in the *z*-direction. By solving the system of equations, all unknowns can be obtained. A drawback of obtaining the unknowns in this way, is that the size of the matrix **E** depends on the number of harmonics that is considered for the problem and, on the number of regions that build up the model. Hence, with a certain amount of available computational memory, the number of harmonics has to be compromised when regions are added to the problem.

To develop a more memory efficient manner of solving, the scattering matrix approach presented in [12] is applied. With this solving method, the unknowns of all regions between the top and bottom air region of the problem are eliminated. For the considered model formulations, the solution in the non-periodical direction is described by two waves, where one wave is traveling in the positive and the other in the negative *z*-direction, as presented in (1) and (2). The scattering matrix, visually represented in Figure 1, couples the incoming waves of a region to the outgoing waves of that region. The sub-matrices **S**<sup>11</sup> and **S**<sup>22</sup> represent the transfer, in negative and positive *z*-direction respectively, from incoming to outgoing waves and **S**<sup>12</sup> and **S**<sup>21</sup> represent the reflection at each side of the region. The method described in [12] is adapted, since sources can be located inside a region. The source terms in scattering vector form are represented by a vector **t** which describes the contribution of the sources to the outgoing waves of the region. By combining all scattering matrices of the regions in the model, a scattering matrix for electromagnetic configuration under consideration can be obtained.

**Figure 1.** Visual representation of the scattering formulation for a single region.

To derive the sub-matrices of the scattering matrix, the continuity boundary conditions between a region with the index *i* and the index *i* + 1 are written in matrix form

$$\begin{aligned} \begin{bmatrix} \mathbf{c}\_{i}^{+} \\ \mathbf{c}\_{i}^{-} \end{bmatrix} &= \begin{bmatrix} \mathbf{E}\_{i}\mathbf{A}\_{i} & \mathbf{E}\_{i}\mathbf{B}\_{i}\mathbf{E}\_{i+1}^{-1} \\ \mathbf{B}\_{i} & \mathbf{A}\_{i}\mathbf{E}\_{i+1}^{-1} \end{bmatrix} \cdot \begin{bmatrix} \mathbf{c}\_{i+1}^{+} \\ \mathbf{c}\_{i+1}^{-} \end{bmatrix} \\ &+ \begin{bmatrix} \mathbf{E}\_{i} & \mathbf{0} \\ \mathbf{0} & \mathbf{I} \end{bmatrix} \cdot \begin{bmatrix} \mathbf{d}\_{1,i} \\ \mathbf{d}\_{2,i} \end{bmatrix} .\end{aligned} \tag{4}$$

where **I** is the identity matrix and

$$\mathbf{A}\_{i} = \frac{1}{2} \left( \mathbf{Q} \mathbf{u}\_{i}^{-1} \; \mathbf{Q} \mathbf{u}\_{i} + \mathbf{Q} \mathbf{v}\_{i}^{-1} \; \mathbf{Q} \mathbf{v}\_{i} + 1 \right),\tag{5}$$

$$\mathbf{B}\_{i} = \frac{1}{2} \left( \mathbf{Q}\_{\mathbf{u}\_{i}}{}^{-1} \mathbf{Q}\_{\mathbf{u}i+1} - \mathbf{Q}\_{\mathbf{v}\_{i}}{}^{-1} \mathbf{Q}\_{\mathbf{v}i+1} \right),\tag{6}$$

$$\mathbf{E}\_i = \mathbf{c}^{\lambda\_i h\_i},\tag{7}$$

$$\mathbf{d}\_{1,i} = \frac{1}{2} \begin{bmatrix} \mathbf{Q}\_{\mathbf{u}\_{i}}{}^{-1} & \mathbf{Q}\_{\mathbf{v}\_{i}}{}^{-1} \end{bmatrix} \cdot \begin{bmatrix} \mathbf{P}\_{u,i+1} - \mathbf{P}\_{u,i} \\ \mathbf{P}\_{v,i+1} - \mathbf{P}\_{v,i} \end{bmatrix} \tag{8}$$

$$\mathbf{d}\_{2,i} = \frac{1}{2} \begin{bmatrix} \mathbf{Q} \mathbf{u}\_i^{-1} & -\mathbf{Q} \mathbf{v}\_i^{-1} \end{bmatrix} \cdot \begin{bmatrix} \mathbf{p}\_{u,i+1} - \mathbf{p}\_{u,i} \\ \mathbf{p}\_{v,i+1} - \mathbf{p}\_{v,i} \end{bmatrix} \tag{9}$$

where *hi* is the height of region *i* in the *z*-direction. The scattering relation, denoted with index *i* − 1, between the top region with index 0 and an arbitrary region *i* is given by

$$
\begin{bmatrix} \mathbf{c}\_{i}^{+} \\ \mathbf{c}\_{0}^{-} \end{bmatrix} = \begin{bmatrix} \mathbf{S}\_{11,i-1} & \mathbf{S}\_{12,i-1} \\ \mathbf{S}\_{21,i-1} & \mathbf{S}\_{22,i-1} \end{bmatrix} \cdot \begin{bmatrix} \mathbf{c}\_{0}^{+} \\ \mathbf{c}\_{i}^{-} \end{bmatrix} + \begin{bmatrix} \mathbf{t}\_{1,i-1} \\ \mathbf{t}\_{2,i-1} \end{bmatrix} . \tag{10}
$$

The scattering matrices, denoted by **S***<sup>i</sup>* and **t***i*, for the scattering relation between region 0 and *i* + 1 can now be defined, thereby removing the unknowns of region *i* from the problem as shown in Figure 1.

This relation is given by

$$
\begin{bmatrix} \mathbf{c}\_{i+1}^{+} \\ \mathbf{c}\_{0}^{-} \end{bmatrix} = \begin{bmatrix} \mathbf{S}\_{11,i} & \mathbf{S}\_{12,i} \\ \mathbf{S}\_{21,i} & \mathbf{S}\_{22,i} \end{bmatrix} \cdot \begin{bmatrix} \mathbf{c}\_{0}^{+} \\ \mathbf{c}\_{i+1}^{-} \end{bmatrix} + \begin{bmatrix} \mathbf{t}\_{1,i} \\ \mathbf{t}\_{2,i} \end{bmatrix} \tag{11}
$$

where the sub-matrices of **S***<sup>i</sup>* and vectors **t**1,*<sup>i</sup>* and **t**2,*<sup>i</sup>* can be computed by

$$\begin{aligned} \mathbf{S}\_{11,i} &= \left(\mathbf{A}\_{i}\mathbf{E}\_{i}^{-1}\mathbf{S}\_{12,i-1}\mathbf{B}\_{i}\right)^{-1}\left(\mathbf{E}\_{i}^{-1}\mathbf{S}\_{11,i-1}\right), \\ \mathbf{S}\_{12,i} &= \left(\mathbf{A}\_{i}\mathbf{E}\_{i}^{-1}\mathbf{S}\_{12,i-1}\mathbf{B}\_{i}\right)^{-1} \\ &\quad \left(\mathbf{E}\_{i}^{-1}\mathbf{S}\_{12,i-1}\mathbf{A}\_{i}\mathbf{E}\_{i+1}^{-1} - \mathbf{B}\_{i}\mathbf{E}\_{i+1}^{-1}\right), \\ \mathbf{S}\_{21,i} &= \mathbf{S}\_{22,i-1}\mathbf{B}\_{i}\mathbf{S}\_{11,i} + \mathbf{S}\_{21,i-1}, \\ \mathbf{S}\_{22,i} &= \mathbf{S}\_{22,i-1}\mathbf{B}\_{i}\mathbf{S}\_{12,i} + \mathbf{S}\_{22,i-1}\mathbf{A}\_{i}\mathbf{E}\_{i+1}^{-1}, \end{aligned} \tag{12}$$

and

$$\begin{split} \mathbf{t}\_{1,i} &= \left(\mathbf{A}\_{i}\mathbf{E}\_{i}^{-1}\mathbf{S}\_{12,i-1}\mathbf{B}\_{i}\right)^{-1} \\ & \quad \left(\mathbf{E}\_{i}^{-1}\mathbf{S}\_{12,i-1}\mathbf{d}\_{2,i} + \mathbf{E}\_{i}^{-1}\mathbf{t}\_{1,i-1} - \mathbf{d}\_{1,i}\right), \\ \mathbf{t}\_{2,i} &= \mathbf{S}\_{22,i-1}\mathbf{B}\_{i}\mathbf{t}\_{1,i} + \mathbf{S}\_{22,i-1}\mathbf{d}\_{2,i} + \mathbf{t}\_{2,i-1}. \end{split} \tag{13}$$

In these equations, stability is ensured by the formulation where only decaying exponential functions are used. Starting from an initial scattering matrix **S**0, which is an identity matrix, and vector **t**0, which contains zeros, the scattering formulation of all regions can be found by repeating (12) and (13). Now, let the air region on the bottom of the problem be denoted with index *I*, so a total of *I* + 1 regions is assumed. The unknowns of the top and bottom region are now related through

$$
\begin{bmatrix} \mathbf{c}\_I^+ \\ \mathbf{c}\_0^- \end{bmatrix} = \mathbf{S}\_{I-1} \cdot \begin{bmatrix} \mathbf{c}\_0^+ \\ \mathbf{c}\_I^- \end{bmatrix} + \mathbf{t}\_{I-1}. \tag{14}
$$

Because all field components of the bottom region have to disappear when *z* goes to minus infinity, the vector of unknowns **c**− *<sup>I</sup>* has to equal zero. When applying the same analysis on the top region, it is obtained that also **c**<sup>+</sup> <sup>0</sup> has to equal zero. The non-zero unknowns of the top and bottom air region are then, using (14), calculated through

$$
\begin{bmatrix} \mathbf{c}\_I^+ \\ \mathbf{c}\_0^- \end{bmatrix} = \mathbf{t}\_{I-1}. \tag{15}
$$

The unknowns of all other regions can be computed from the unknowns **c**− <sup>0</sup> and **<sup>c</sup>**<sup>+</sup> *<sup>I</sup>* , by consequently applying (4).

#### **4. Application to an Electromagnetic Configuration**

The scattering approach is applied to the electromagnetic configuration shown in Figure 2, where a rectangular slab of magnetic material is placed above three permanent magnets. The configuration is divided into five regions in the *z*-direction, as depicted in Figure 2. The model is periodic in both the *x*- and *y*- direction and the permeability inside Regions 2 and 4 varies as a function of position. The configuration was analyzed in [10], where also the expressions for the magnetic field components have been obtained and all geometric properties have been given. The magnetic properties are given in Table 1. In [10], however, the unknowns where obtained by solving a linear system of equations as presented in (3).

**Figure 2.** Electromagnetic configuration consisting of three magnets with a slab of magnetic material above [10].


**Table 1.** Model parameters.

Because only the magneto-static fields have to be considered for the configuration, the magnetic scalar potential, *ψ*, is introduced to obtain the expressions for the magnetic field quantities. In the boundary conditions, the continuity of this scalar potential and the *z*-component of the magnetic flux density are forced. Therefore, **u** = **bz** and **v** = *ψ* in (1) and (2), respectively. The expressions for **bz** and *ψ* are described in [10]. Note that, for a problem where also an electric field exists, **u** and **v** would consist of the tangential components of the *H*- - and - *E*-field, respectively, and their size would be twice that compared to the current magneto-static problem. Hence, the expressions for **u** and **v** and their size depend on the problem at hand.

The number of harmonics that is used to calculate the solution to the configuration of Figure 2 is denoted by *L*. The number of harmonics used in the *x*- and *y*-directions is denoted by *N* and *M*, respectively, and thus the total number of harmonics *L* = *NM*. With the solving method of (3), the matrix that has to be inverted has a size of 8*L* ×8*L*, while with the scattering approach the maximum size of any matrix that is inverted is equal to *L* × *L*. This means that the matrix size is reduced by 87.5%. Moreover, with the classical solving method, adding a region to the problem would mean that the matrix **E** becomes larger by 2*L* in both directions. In the scattering approach, the maximum size of any matrix would be *L* × *L*, independent of the number of regions that makes the method especially beneficiary when the number of regions in a model is relatively large. The solving method using (3) is performed in Matlab R2017a using the command 'A\b', which means that the method of inversion is automatically chosen. It should be noted that properties reducing the computational effort of the inversion of a matrix, such as zeros in the off-diagonal entries and matrix symmetry are limited in the boundary condition matrix of (3), in contrast to the matrices inverted in (5), (6), (8), (9), (12) and (13).

In Figure 3a, the magnetic field strength, computed with the scattering matrix approach with a total of 6561 harmonics, in the center of Region 2 is depicted. The difference with the result computed by the finite element method (FEM) is shown in Figure 3b. In the FEM analysis, a total of 2,384,628 second-order mesh elements is used. It has been verified that, with the applied mesh, the error of the FEM (compared to a FEM result obtained with a denser mesh) has converged to an error less than 1%. In [10], a total of 1089 harmonics was used to obtain the same result. The maximum relative error with respect to FEM with 6561 harmonics is equal to 9.2%, where, computed with 1089 harmonics, the error was equal to 17.4%. For problems where the electromagnetic properties of a region are varying with a high spatial frequencies, the addition of harmonics can be necessary to obtain accurate field results.

Using the hardware available (quadcore Intel Core i7-4790 with 32 GB RAM), the computation time for the 'classical' solving method for a total of 1089 harmonics is equal to 20.4 s. With the developed scattering matrix approach, the computation time, for the same number of harmonics, is equal to 12.6 s. Hence, the number of calculations is increased in the scattering matrix approach, however, because all manipulations are performed on matrices of size *L* × *L*, the computation time is not significantly increased or even decreased compared to the classical solving method. On the same hardware, the FEM analysis takes around 125 s to mesh and solve the electromagnetic problem. The semi-analytical approach obtains the results for the electromagnetic configuration of Figure 2 faster than the FEM; however, the solving time of the semi-analytical approach depends on the number of harmonics and number of regions in the model. This means that, for models where the geometry contains high spatial frequencies (relatively small details with respect to the periodic width), and the number of harmonics needs to be large to accurately compute the fields, and the advantage over FEM regarding computation time will decrease.

**Figure 3.** Magnetic field strength in the *z*-direction on an *xy*-plane in the center of the region with magnetic material (Region 2). (**a**) calculated *Hz* with the scattering matrix approach; (**b**) difference in calculated *Hz* between the scattering matrix approach and the finite element method.

#### **5. Conclusions**

The paper presents the scattering matrix approach to solve the boundary conditions of a multiple region electromagnetic problem, with Fourier based solutions. Compared to the classical solving method, where all boundary conditions are collected in a single matrix, the scattering approach is more memory efficient. The boundary conditions for each region are reformulated, so that all unknowns, except the ones from the top and bottom region, can be removed from the problem. As a result, the maximum size of any matrix that is used during calculations is determined by the number of harmonics and, more importantly, independent of the number of regions inside the problem. Therefore, the scattering approach is especially beneficial for problems with a relatively large number of regions. Sources located inside a certain region of the electromagnetic configuration are taken into account in the described scattering formulation. Furthermore, due to the general formulation used, the approach can be used for a variety of electromagnetic problems with continuity boundary conditions between regions.

**Author Contributions:** The theory presented in this paper has been developed by C.H.H.M.C. The analysis of the results has been performed in cooperation with J.W.J. and E.A.L. The paper was written by C.H.H.M.C. and contributions and improvements to the content have been made by J.W.J. and E.A.L.

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

#### **References**


© 2018 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

*Article*
