Technical Article

Homogenization and Model Order Reduction
for Fast FE analysis
Abstract — When sing the standard FEM, it is still challenging toanalyze the eddy currents in each strand of a Litz wire and alsoeach magnetic particle in soft magnetic composite. Because theradii of the strand and magnetic particle are much smaller thanthe device size, fine FE discretization necessary for the eddy
current analysis results in quit a large equation system. This is a
so-called multi-scale problem. In this article, it is shown that thehomogenization method, which models material with finestructure as homogenous one with effective material constants, isfairly effective for the analysis of such multi-scale problems.
Moreover, the standard FE analysis of electric apparatus suchas motors and inductors need too long computational time toperform dynamic simulation coupled with power circuits. This
article shows that the model order reduction can not only reduce
the FE equations to accelerate the analysis but also allows us tosynthesize the equivalent circuits directly from the FE equations
which can easily be coupled with external circuits.

I. INTRODUCTION
The FE analysis of electric and electronic apparatus such as
motors, transformers and inductors is now indispensable for
their research and development. We can perform the design
optimization and also three-dimensional field computations
using FEM within acceptable duration of time. The standard
FEM would need, however, unacceptably long computationaltime to analyze materials composed of multiple elements which
are much smaller than the overall material size. For example,
when we analyze the eddy currents in a Litz wire which iscomposed of a number of stranded wires using the standard
FEM, the resultant equation systems would be far too large tobe solved in reasonable duration of time under usualcomputational environments. We encounter similar difficultywhen computing the eddy current losses in the soft magnetic
composite (SMC) which consists of many conductive magnetic
particles coated by insulation layers. This article will show thatthe homogenization method is fairly effective for the eddy
current analysis of these problems considering the skin andproximity effects. In this method, material such as multi-turncoil, Litz wire and SMC is modeled as a homogenous material
with macroscopic complex permeability. This method allows us,
therefore, to perform the FE analysis with rather coarse FE
meshes, and to reduce the computational time.
The standard FE analysis is also too time consuming toperform dynamic simulations of electric apparatus such as
motors and inductors coupled with power circuits. For thisreason, they are often analyzed with behavior modeling or
equivalent circuits to reduce the computational time. Thesemodels can be realized by hardware to perform HILS
(hardware-in-the-loop-simulation). However, it would be
difficult to evaluate eddy current and hysteresis losses using
these approaches. Moreover, the computational accuracy issometimes unsatisfactory. In this article, it will be shown thatthe model order reduction (MOR) can be the alternative to theseconventional methods. The MOR approach converts theoriginal FE equation to reduced one with small number of
unknowns. Moreover, MOR enables to synthesize theequivalent circuit of electric apparatus directly from its FE
model.

In the next section, the homogenization method based on themacroscopic complex permeability will be formulated and theresults of FE analysis of multi-turn coils and SMC will bereported. Moreover, the integral equation approach based on thecomplex permeability to perform three-dimensional analysis ofLitz wire will be described. In the third section, the MOR
techniques will be surveyed and numerical results will bereported. Moreover, the synthesis of equivalent circuits basedon MOR will be discussed. The last section provides theconclusions and outlooks.
II. HOMOGENIZATION METHODS
A. Analysis of multi-turn and Litz wire coils [1]When analyzing electric apparatus which has rich higher
harmonics in the coil windings, the eddy current effects not onlyin the magnetic cores but also in wires have to be considered.
There are two effects relevant to eddy currents: skin andproximity effects. The skin effect means that the eddy currents
parallel to the wire axis are induced by the self-current so thatthe current concentrates near the wire surface. On the other hand,the proximity effect means that the anti-parallel eddy currents
are induced by the magnetic fluxes generated by the ambient
currents, and consequently the total current distribution
becomes non-symmetric in the wire cross section. To reducethese effects, the Litz wire composed of fine twisted strands is
widely used. Figure 1 depicts the Litz wire coil used for awireless power transfer device.
In the conventional FE analysis of eddy currents in themulti-turn coils and in Litz wire coils, the wires have to be
subdivided into so fine elements that they are smaller than theskin depth. This results in large number of unknowns in the FE
equation. Especially in the three-dimensional analysis, it would
be difficult to solve the resultant large equation system evenusing high-performance computers. To circumvent thisproblem, the homogenization method which will be described
below is fairly effective.

Fig.1 Litz wire coil
To reduce the eddy current losses, the Litz wire coil is formed frommany stranded fine wires.

Fig.2 Stranded wire and isolated wire in time-harmonic uniform field
A wire included in a bundle of stranded wires carrying alternating currents
is modeled as an isolated wire immersed in a uniformed field generated by
the surrounding wire currents.

�−𝑗𝑗𝑗𝑗𝑗𝑗𝑗𝑗, 𝜁𝜁 = 𝑘𝑘𝑘𝑘. The first and second terms in (1) express thefields relevant to the skin and proximity effects. The proximity
effect means that the alternating fields generated by thesurrounding wires induce the anti-parallel eddy currents in thewire as shown in Fig.3. These eddy currents, which can be
represented by a magnetic dipole moment m, generate magnetic
fields. Hence the eddy currents and magnetic dipole moment
have to be determined self-consistently. Moreover, because theinduced eddy currents are anti-parallel to each other, thedistribution of the total current becomes non-symmetric asshown in Fig.3 unlike the result of the skin effect. This nonsymmetry is represented by the second terms in (1).
The fields outside the wire can be expressed as
𝐸𝐸𝑧𝑧 = �𝑞𝑞1+ 𝑗𝑗𝑗𝑗𝑗𝑗𝐵𝐵0 � cos 𝜃𝜃 ,𝑟𝑟

𝐻𝐻𝜃𝜃 = −Fig.3 Skin (left) and proximity (right) effects
Left: the circular magnetic flux (red) generated by the self-current (blue)
induces the eddy currents (orange). As a result, the total current (yellow)
concentrates near the wire surface. Right: The magnetic flux generated bysurrounding currents induce the anti-parallel eddy currents. As a result, thetotal current (yellow) has the crescent-shaped distribution.

𝑟𝑟 ≥ 𝑎𝑎(2a)𝑟𝑟 ≥ 𝑎𝑎(2b)− 1⁄𝑎𝑎 𝑝𝑝11
� � � = 𝑗𝑗𝑗𝑗𝑗𝑗𝐵𝐵0 � �
1⁄𝑎𝑎 𝑞𝑞1
1

(3)

𝑗𝑗
𝑞𝑞1
�− 2 + 𝑗𝑗𝜔𝜔𝐵𝐵0 � cos 𝜃𝜃 ,
𝑟𝑟
𝜔𝜔𝜇𝜇0

The first and second terms in (2) express the dipole field due tothe proximity effect and the uniform external field, respectively.
Note here that the components independent of 𝜃𝜃, which will not
be involved in the following discussions, are omitted forsimplicity. The coefficient 𝑝𝑝0 can be determined from theAmpere law for a given total current 𝐼𝐼 : 𝑝𝑝0 =𝑘𝑘𝑘𝑘 ⁄(2𝜋𝜋𝜋𝜋𝜋𝜋𝐽𝐽1 (𝛼𝛼)) , 𝛼𝛼 = 𝑘𝑘𝑘𝑘 , while 𝑝𝑝1 , 𝑞𝑞1 can be determinedfrom continuity in the tangential component of E and H, that is�

𝐽𝐽1 (𝛼𝛼)′𝛼𝛼𝐽𝐽1 (𝛼𝛼)⁄𝜇𝜇𝑟𝑟

The dipole fields, which represent the diamagnetic property
due to the anti-parallel eddy currents along the wire, can also be
expressed in terms of the magnetic dipole 𝒎𝒎 in the form
Fig.4 Homogenization of stranded wires
The stranded wire is modeled as a homogenous material with the complex
macroscopic permeability given by (8) defined on the wire cross section.
For the three-dimensional analysis, the permeability tensor [𝜇𝜇] =diag[〈𝜇𝜇̇ 〉, 〈𝜇𝜇̇ 〉, 𝜇𝜇] is introduced, where 𝜇𝜇 is the original wire permeability.

Let us consider a bundle of wires composed of many parallel
round conductors, as shown in Fig.2, which carry alternating
currents of angular frequency 𝜔𝜔. The effect of twisting will beconsidered later. Let us pay attention to one of the conductors
shown in the right. The magnetic field 𝑩𝑩0 in which thisconductor is immersed would approximately be uniform, atleast locally. This assumption is justified from the well-known
fact that the magnetic field in a cylindrical hole in a cylindricalconductor carrying a uniform current is uniform [2].
The quasi-static fields inside and outside the wire can beanalytically determined by solving the Maxwell equation. The
fields inside the conductor can be written as
𝐸𝐸𝑧𝑧 = 𝑝𝑝0 𝐽𝐽0 (𝜁𝜁) + 𝑝𝑝1 𝐽𝐽1 (𝜁𝜁) cos 𝜃𝜃 ,𝐻𝐻𝜃𝜃 =

𝑗𝑗𝑗𝑗[𝑝𝑝 𝐽𝐽 (𝜁𝜁) − 𝑝𝑝1 𝐽𝐽1′ (𝜁𝜁) cos 𝜃𝜃],𝜔𝜔𝜔𝜔 0 1

𝑟𝑟 ≤ 𝑎𝑎𝑟𝑟 ≤ 𝑎𝑎

(1a)

(1b)

where (𝑟𝑟, 𝜃𝜃, 𝑧𝑧) is the cylindrical coordinate in which z-axis islocated along the wire axis, 𝐽𝐽0 and 𝐽𝐽1 denote the zeroth and firstorder Bessel functions, prime denotes derivative and 𝑘𝑘 =𝐴𝐴𝑧𝑧 = −𝜇𝜇0 𝑚𝑚cos 𝜃𝜃
2𝜋𝜋 𝑟𝑟(4)

By comparing the first term of (2a) and (4) where 𝐸𝐸𝑧𝑧 = −𝑗𝑗𝑗𝑗𝐴𝐴𝑧𝑧 ,we have the magnetization
𝑀𝑀 =𝑚𝑚
𝐵𝐵0 𝜇𝜇𝑟𝑟 𝐽𝐽1 (𝛼𝛼) − 𝛼𝛼𝐽𝐽1 ′(𝛼𝛼)= 2𝜋𝜋 𝑎𝑎22𝜇𝜇0 𝜇𝜇𝑟𝑟 𝐽𝐽1 (𝛼𝛼) + 𝛼𝛼𝐽𝐽1 ′(𝛼𝛼)
𝜋𝜋𝑎𝑎

(5a)

By taking the static limit in (5), we obtain the well-known result
for the magnetostatic field as
𝑀𝑀 = 2𝐵𝐵0 𝜇𝜇𝑟𝑟 − 1
𝜇𝜇0 𝜇𝜇𝑟𝑟 + 1(5b)Here let us define the complex permeability 𝜇𝜇̇ so that (5b)holds at any frequency. To do so, we replace 𝜇𝜇𝑟𝑟 with 𝜇𝜇̇ 𝑟𝑟 in (5b)and compare (5a) with (5b) to obtain𝜇𝜇̇ 𝑟𝑟 (𝜔𝜔) = 𝜇𝜇𝑟𝑟

𝐽𝐽1 (𝛼𝛼)=𝛼𝛼𝐽𝐽1′ (𝛼𝛼)𝜇𝜇𝑟𝑟
𝐽𝐽0 (𝛼𝛼)𝛼𝛼−1𝐽𝐽1 (𝛼𝛼)(6)

Now the proximity effect in a wire can be computed using (6)without fine FE discretization. However, FE discretization of
the stranded wires and air region among them would still giverise to long computational time. To circumvent this problem,
we use the homogenization method which models material
composed of multiple elements as homogenous one.

It is known that the magnetostatic property of SMC can beevaluated at good accuracy [3] using the homogenizedpermeability given by〈𝜇𝜇𝑟𝑟 〉 = 1 +𝜂𝜂(𝜇𝜇𝑟𝑟 − 1)
1 + 𝑁𝑁(1 − 𝜂𝜂)(𝜇𝜇𝑟𝑟 − 1)(7)

where 𝜂𝜂 and 𝑁𝑁 denote the volume fraction (or filling rate) anddiamagnetic constant which equals 1/2 for the round wire.
Equation (7) is referred to as the Ollendorff formula [4] whose
derivation is given in [1]. The Clausius-Mossotti formula basedon the molecule field, often used for computation ofhomogenized permittivity of composite dielectric materials,
and also the Maxwell-Garnett formula, used in microwave
engineering, are shown equivalent to (7). The discussions onthese equivalence are included in Ref.[1]. What we have to do
to obtain the macroscopic complex permeability 〈𝜇𝜇̇ 〉 is just to
substitute 𝜇𝜇̇ 𝑟𝑟 in (6) into 𝜇𝜇𝑟𝑟 in (7). The result is〈𝜇𝜇̇ 𝑟𝑟 (𝜔𝜔)〉 = 1 +(8)

Re〈𝜇𝜇̇ 𝑟𝑟 〉1.20

real part of complex permeability2𝜂𝜂(𝜇𝜇̇ 𝑟𝑟 (𝜔𝜔) − 1)
2 + (1 − 𝜂𝜂)(𝜇𝜇̇ 𝑟𝑟 (𝜔𝜔) − 1)

1.00
0.80
0.60
0.40
0.20
0.00

0.1

1
nomalized wire radius

0.1

1

0.00

10

𝑎𝑎⁄𝛿𝛿

10

Fig.6 Frequency dependence of 50-turned coil impedance
An enameled copper wire, radius 0.15 mm, cover thickness 0.03 mm iswound around a bobbin, radius 30 mm, manufactured by a 3D printer.

Now the stranded wire composed of fine wires can be modeled
as a homogenous material as shown in Fig.4. Note there that〈𝜇𝜇̇ (𝜔𝜔)〉 is defined on the wire cross-section, whereas thepermeability in the wire-axis direction can be approximated by
the original permeability 𝜇𝜇 because the external magnetic flux
along the wire axis would have little effects.
The frequency characteristics for different volume fractions
are plotted in Fig.5, where the abscissa denotes the wire radiusnormalized by the skin depth 𝛿𝛿 = �2⁄𝜇𝜇𝜇𝜇𝜇𝜇 . The real andimaginary parts represent the permeability possessing thediamagnetic property due to the induced eddy currents and theeddy current loss, respectively. The former monotonously
decreases with frequency, while the latter has a peak. This
property is also found in the Debye relaxation of dielectricmaterial. The reader would wonder that the curves for rather
high values of 𝜂𝜂, which cannot be realized when using the roundwires, are included in Fig.5. Such a high filling rate can berealized when using the wire with square-shaped cross section
whose diamagnetic constant N is also approximately 1/2.
For the test of this homogenization method, the impedanceof 50-turned coil, shown in Fig.6, has been experimentallymeasured and analyzed using the proposed method varying
frequency for comparison. In the analysis, we solve themagnetostatic equation for the vector potential given byimaginary part

-0.10
-0.20

(9)

𝑅𝑅0 𝛼𝛼𝛼𝛼0 (𝛼𝛼) 2
𝑗𝑗𝑗𝑗
|𝐼𝐼|
� 𝜇𝜇|𝑯𝑯|2 𝑑𝑑𝑑𝑑 +4 𝐽𝐽1 (𝛼𝛼)2(10)on the plane containing the cross section of the coil. We use
FEM where the permeability in the coil region is set to 𝜈𝜈 =1⁄〈𝜇𝜇̇ (𝜔𝜔)〉 whereas 𝜈𝜈 = 1⁄𝜇𝜇0 in the other regions. From thecomputed complex power,

-0.30
-0.40
-0.50

∇ ⋅ (𝜈𝜈∇𝐴𝐴z ) = −𝐽𝐽z

Im〈𝜇𝜇̇ 𝑟𝑟 〉normalized wire radius𝑎𝑎⁄𝛿𝛿Fig.5 Profiles of 〈𝜇𝜇̇ 𝑟𝑟 〉 with different filling rates
The abscissa is the ratio of the wire radius a to the skin depth δ. The
original wire permeability 𝜇𝜇𝑟𝑟 is set to 1. The red, blue, orange, purple
curves are plotted for 𝜂𝜂 = 1, 0.75,0.5, 0.25.

𝑃𝑃̇ =

Ω

we compute the impedance 𝑍𝑍̇ = 2Re�𝑃𝑃̇ �⁄𝐼𝐼 2 , where 𝑅𝑅0 is the
DC resistance of the wire. The real parts of the first and second
terms in (10) represent the eddy current losses due to theproximity and skin effects, respectively. Note here that the eddy
current losses and also its diamagnetic effects can be evaluatedfrom magnetostatic analysis where the permeability in the coilregion is simply replaced with (8). Of course, if there are otherconductors in which eddy currents flow, we have to solve thequasi-static electromagnetic equation instead of (9). However,also in this case, the coil region still can be modeled by (9).
From the results plotted in Fig.6, we conclude that the proposed
method is sufficiently accurate.
Until now, we have considered parallel conductors. If thetwisting pitch of stranded wires is sufficiently longer than thewire radius a, the above-mentioned method is valid. Unless thisis the case, we have to consider the twisting effect. One of thepossible ways to do so is to introduce the tensorial macroscopic
complex permeability whose entities are determined from thewire direction. Namely, assuming that the stranded wires are
locally parallel to each other, the permeability tensor is written
as [𝜇𝜇] = diag[〈𝜇𝜇̇ (𝜔𝜔)〉, 〈𝜇𝜇̇ (𝜔𝜔)〉, 𝜇𝜇] in the local coordinates
where the third axis is parallel to the wire. In the FE equation,
this locally defined tensor is rewritten in the global coordinates.
There is another way to treat the twisted wires; a conducting
curved wire is modeled as a curved line on which onedimensional integral equation is formulated. The latter method
will be mentioned below.
B. One-dimensional Integral Equation [5]As we have seen above, the proximity effect in a wire ismodeled by the complex permeability (6), by which the FE
discretization of the wire cross section is no more necessary.
This means that the each wire can be modeled as a curved lineon the plane perpendicular to which 𝜇𝜇̇ and magnetization M
generated by the anti-parallel eddy currents due to the proximity
effect are defined (see Fig. 7). The proximity effect can beanalyzed by determining M distributed along the curved line
self-consistently. If we formulate the integral equation for M
along the curved line, we do not need to consider the air region
among the twisted wires. We know the magnetostatic integral
equation for M in the magnetic material which is described inRef. [6]. By extending this integral equation, we obtain the onedimensional integral equation given by𝑴𝑴(𝒙𝒙)= 𝝉𝝉 × ��� ℊ𝑴𝑴(𝒙𝒙′ ) 𝑑𝑑𝑣𝑣 ′ + 𝑯𝑯(𝒙𝒙′ )� × 𝝉𝝉�𝜇𝜇̇ 𝑟𝑟 (𝜔𝜔) − 1
Ω𝑤𝑤ℊ𝑴𝑴(𝒙𝒙) = −(𝑴𝑴(𝒙𝒙) ⋅ 𝑹𝑹)𝑹𝑹
1 𝑴𝑴(𝒙𝒙)� 3 −3�
4𝜋𝜋 𝑅𝑅𝑅𝑅5(11a)

𝑅𝑅0 𝛼𝛼𝛼𝛼0 (𝛼𝛼) 2
𝑗𝑗𝑗𝑗
|𝐼𝐼|
� 𝑱𝑱∗ ⋅ 𝑨𝑨 𝑑𝑑𝑑𝑑 +4 𝐽𝐽1 (𝛼𝛼)2(12)𝑴𝑴(𝒙𝒙′ ) × 𝑹𝑹 𝑱𝑱(𝒙𝒙′ )
𝜇𝜇0
� �
+
� 𝑑𝑑𝑣𝑣 ′4𝜋𝜋 𝛺𝛺𝑤𝑤
𝑅𝑅3
𝑅𝑅

(13)

Ω𝑤𝑤

where 𝐼𝐼, 𝑱𝑱 denote the total current and corresponding current
density, and A is computed from𝑨𝑨(𝒙𝒙) =(a) parallel wires(b) multiple strand
Fig.7 Modeling of stranded wires as curved lines.
The magnetization M distributed along the curved lines can be
determined by solving (11).

(11b)

which describes the proximity effect in the stranded wires,
where 𝛕𝛕, Ω𝑤𝑤 , 𝑯𝑯 denote the unit vector parallel to the wire, wire
domain and external magnetic field, respectively, and 𝑹𝑹 = 𝒙𝒙 −𝒙𝒙′. Because M exists on the plane perpendicular to 𝛕𝛕, we takethe vector triple product in (11a). To numerically solve (11), we
discretize the curved line into line segments to each of which M
is assigned. The discretization of (11) leads to a system of
equations containing a full matrix. This equation is can beeffectively solved by iterative methods such as Jacobi andGauss-Seidel methods.
After solving (11), we can evaluate the complex power from𝑃𝑃̇ =Figure 8 shows the eddy current distribution in a multiple strand,
shown in Fig.7 (b), computed from 𝑴𝑴 obtained by solving (11).
Moreover, it is shown that the eddy current losses due to theproximity effect computed by this method and homogenizationbased FEM mentioned in II.A are in good agreement [5]. It isstressed here that the method based on the one-dimensional
integral equation can analyze Litz wire three-dimensional
structure. Moreover, this method is also valid even when there
exit magnetic and conducting materials (see [5] for details). If
there exit a number of twisted wires, say, more than one
hundred, the simple discretization of (11) would impose large
computational burden with respect to time and memory. We
would introduce the fast multipole expansion or H-matrix
method to solve this system.

Fig.8 Eddy current distribution in multiple strand composed of 49 wires
After solving (11), the magnetic field 𝑩𝑩0 outside the wire is determinedfrom 𝑩𝑩0 = 𝜇𝜇0 (𝜇𝜇̇ 𝑟𝑟 + 1)𝑴𝑴⁄2(𝜇𝜇̇ 𝑟𝑟 − 1). Then, the unknown coefficient 𝑝𝑝1
in (1a) can be determined from (3). The eddy currents in this figure are
finally computed from 𝑱𝑱 = σ𝑬𝑬.

C. Time-domain Analysis [7]We consider here how we apply the above-mentioned
homogenization method to time-domain problems. This isrequired when we analyze, for example, the windings used inwireless power transfer devices which include nonlinear power
circuits. Moreover, the present homogenization method canalso be applied to SMC, as will be mentioned later. When we
consider the magnetic nonlinearity in the magnetic particles in
SMC, we also have to perform time-domain analysis.

The simplest way to do so would be to employ the
convolution integral including the inverse Laplace transform of
the reciprocal of 〈𝜇𝜇̇ 〉. The direct inverse Laplace transformation
seems, however, uneasy. To circumvent this problem, let us
consider the expansion of α 𝐽𝐽0 (𝛼𝛼)⁄𝐽𝐽1 (𝛼𝛼) into contined fractions.
The result is
𝜇𝜇̇ 𝑟𝑟 (𝜔𝜔) = 𝜇𝜇𝑟𝑟

2−

1

𝛼𝛼 2
𝛼𝛼 2
4−
𝛼𝛼 2
6−
8−⋯

(14)

When we truncate the continued fraction in (14), we obtain therational fraction which is referred to as the Padé approximation
as follows:𝜇𝜇̇ 𝑟𝑟 = 𝜇𝜇𝑟𝑟

𝑐𝑐0 + 𝑐𝑐1 𝑠𝑠 + 𝑐𝑐2 𝑠𝑠 2 + ⋯ + 𝑐𝑐𝑞𝑞−1 𝑠𝑠 𝑞𝑞−1
𝑑𝑑0 + 𝑑𝑑1 𝑠𝑠 + 𝑑𝑑2 𝑠𝑠 2 + ⋯ + 𝑑𝑑𝑞𝑞 𝑠𝑠 𝑞𝑞(15)

where we replaced 𝑗𝑗𝑗𝑗 with 𝑠𝑠 . Substitution of (15) into the
Offendorff formula (7) results in the partial fraction expansionof the form〈𝜈𝜈̇𝑟𝑟 〉 =

𝑞𝑞

1
𝑘𝑘𝑖𝑖
= 𝜈𝜈∞ + �〈𝜇𝜇̇ 𝑟𝑟 〉𝑠𝑠 + 𝑝𝑝𝑖𝑖
𝑖𝑖=1

(16)

Now we can easily perform the inverse Laplace transform of(16) to obtain𝑞𝑞ℒ −1 〈𝜈𝜈̇𝑟𝑟 〉 = 𝜈𝜈∞ 𝛿𝛿(𝑡𝑡) + � 𝑘𝑘𝑖𝑖 𝑒𝑒 −𝑝𝑝𝑖𝑖𝑡𝑡
𝑖𝑖=1

D. Eddy Current Analysis of Soft Magnetic Composites [9, 10]
The soft magnetic composite, SMC, is composed of small
iron particles embedded in insulation medium. Due to thisstructure, global eddy current cannot flow in SMC, and thus theeddy current loss can be suppressed. Moreover, SMC has
another merit over the conventional silicon steel plates thatthree-dimensional shapes can easily be formed. For this reason,
SMC has widely been used for inductors and reactors. However,as the driving frequency increases, the local eddy currents
flowing in the iron particles come to give rise to non-negligible
losses, which is also the case for the silicon steel sheets. Figure
9 shows the magnetic field distributions in an isolated ironparticle for different frequencies.
The homogenization method mentioned above can also be
applied to analysis of SMC considering the skin and proximity
effects. Let us consider an isolated iron particle of radius a
immersed in the uniform time-harmonic field 𝑩𝑩0 𝑒𝑒 −𝑗𝑗𝑗𝑗𝑗𝑗 .
Introducing the polar coordinates (𝑟𝑟, 𝜃𝜃, 𝜙𝜙) around the center ofthe particle, we can obtain the fields by analytically solving thequasi-static Maxwell equations as follows:
𝐴𝐴𝜙𝜙 = 𝑝𝑝𝐵𝐵0 𝑖𝑖1 (𝜁𝜁)sin𝜃𝜃, 𝑟𝑟 ≤ 𝑎𝑎

𝑞𝑞 𝑟𝑟
𝐴𝐴𝜙𝜙 = 𝐵𝐵0 � 2 + � sin𝜃𝜃, 𝑟𝑟 ≥ 𝑎𝑎𝑟𝑟
2𝑖𝑖1 (𝜁𝜁) =1(sin 𝜁𝜁 − 𝜁𝜁 cos 𝜁𝜁)𝜁𝜁 2(d) imaginary (𝛿𝛿⁄𝑎𝑎 = 5)

(19)

As in the same way for the stranded wires, the coefficients canbe determined from the continuity condition in the electric andmagnetic field on the particle surface, that is𝑖𝑖1 (𝛼𝛼)�−𝑖𝑖1 (𝛼𝛼) + sin(𝛼𝛼)𝑎𝑎 1−𝑎𝑎−2 𝑝𝑝−2 � �𝑞𝑞 � = �2�
2
𝑎𝑎

(20)

where 𝛼𝛼 = 𝑘𝑘𝑘𝑘.
The magnetic field generated by a magnetic moment m inthree-dimensional space can be written as
𝜇𝜇0 𝑚𝑚sin 𝜃𝜃
4𝜋𝜋 𝑟𝑟 2(21)

3𝑚𝑚
3𝐵𝐵0 (2𝜇𝜇 + 𝜇𝜇0 )𝑖𝑖1 (𝛼𝛼) − 𝜇𝜇0 sin 𝛼𝛼=
32𝜇𝜇0 (𝜇𝜇 − 𝜇𝜇0 )𝑖𝑖1 (𝛼𝛼) + 𝜇𝜇0 sin 𝛼𝛼
4𝜋𝜋𝑎𝑎

(22)

3 𝜇𝜇 − 𝜇𝜇0
𝜇𝜇0 𝜇𝜇 + 2𝜇𝜇0

(23)

𝐴𝐴𝜙𝜙 =(c) real (𝛿𝛿⁄𝑎𝑎 = 5)(18b)where 𝑝𝑝, 𝑞𝑞 are coefficients, 𝑘𝑘 = �−𝑗𝑗𝑗𝑗𝑗𝑗𝑗𝑗, 𝜁𝜁 = 𝑘𝑘𝑘𝑘 , and 𝑖𝑖1
denotes the first-order modified Bessel function defined by(17)

Applying (17) to the convolution integral, we can perform thetime-domain analysis. In particular, we employ the RC method[8] by which the convolution integral can recursively beevaluated with reference to the field at just one previous time
step. The authors have found that this approach gives numericalresults which are in good agreement with those computed by
the conventional FE analysis which requires very finediscretization.

(18a)

The first term in (18b) also expresses the dipole field, generated
by the eddy currents in the particle. Hence, by comparing (18a)and (21), we find 𝑚𝑚 = 4𝜋𝜋𝐵𝐵0 𝑞𝑞 ⁄𝜇𝜇0 . The magnetization is thengiven by𝑀𝑀 =When taking the static limit 𝜔𝜔 → 0 in (22), we obtain the wellknown result(c) real (𝛿𝛿⁄𝑎𝑎 = 0.35)(d) imaginary (𝛿𝛿⁄𝑎𝑎 = 0.35)

Fig.9 Magnetic flux lines around a magnetic conductive particle
immersed in uniform time-harmonic magnetic field, where 𝑎𝑎, 𝛿𝛿denote the particle radius and skin depth [11].

𝑀𝑀 =We again define the complex permeability 𝜇𝜇̇ (𝜔𝜔) so that (23)holds at any frequencies. From (22), (23), we finally obtain

𝜇𝜇̇ 𝑟𝑟 (𝜔𝜔) =2𝜇𝜇𝑟𝑟 (𝛼𝛼 − tan 𝛼𝛼)(1 − 𝛼𝛼 2 )tan 𝛼𝛼 − 𝛼𝛼This result has already been reported in Ref. [12]. By
substituting (24) into the Ollendorff formula (7) where we set
𝑁𝑁 = 1⁄3 , we obtain the macroscopic complex permeability〈𝜇𝜇̇ (𝜔𝜔)〉.
When we consider the magnetic saturation in the ironparticles in SMC, we have to perform time-domain analysis. In
this case, the inverse Laplace transform based on the Padé
approximation given by (17) would be invalid because of themagnetic nonlinearity in 𝜇𝜇 . Instead, we compute the timeresponse of SMC using the Cauer-type ladder circuit, shown in
Fig.10, which is directly synthesized from 〈𝜇𝜇̇ (𝜔𝜔)〉 . In this
synthesis, we expand (24) into a continued fraction using
tan 𝑧𝑧
=
𝑧𝑧

1−

1

𝑧𝑧 2
𝑧𝑧 2
3−
5−⋯

R1

(24)

R2
L2

L1

L3

…

Fig.10 Cauer circuit
In the frequency domain, this circuit is equivalent to the macroscopic
complex permeability of SMC given by (26). This circuit can be used
for time domain analysis of SMC cores as well as laminated steel plates.
Moreover, electric apparatus such as inductors and induction heating
device can be modeled by the Cauer circuit as mentioned in III.D. In
this case, the DC resistance of the winding is inserted as R0 in left of L1.
B [T]

(25)

By substituting the resultant continued fraction into theOllendorff formula, we have the rational polynomial of the form(15) for 〈𝜇𝜇̇ (𝜔𝜔)〉. Using the Euclidian algorithm, the rational
polynomial can be expressed in the form of continued fraction〈𝜇𝜇̇ 〉 =𝐿𝐿1 − 𝑅𝑅

1

𝑠𝑠

1

−

1

1
𝐿𝐿2 + ⋯

(26)(a) with considering eddy current in SMC coreThis continued fraction corresponds to the Cauer circuit shownin Fig. 10. It is known [13, 14] that 𝐿𝐿1 in the Cauer circuit isrelevant to the primal magnetic field generated by an externalcurrent which usually overwhelms the fields generated by eddy
currents. For this reason, by introducing proper nonlinear
relationship between the flux and current in 𝐿𝐿1 , we can take themagnetic saturation in SCM particles into consideration. The
circuit equation for the Cauer circuit in Fig.10, which is
assigned in each finite element in the SMC region, is coupledwith the FE equation. In Ref. [10], the time response of SMC
computed by this approach has been shown to agree well withthat computed by the unit-cell method. Figure 11 shows themagnetic distribution in the SMC core with and without
considering the eddy currents in the iron particles. We can see
the clear difference between these distributions.

(b) without considering eddy current in SMC coreFig.11 Magnetic fields in SMC core [see also 10]

E. Laminated Steel Sheets
We can also obtain the analytical expression of the macroscopic
complex permeability 〈𝜇𝜇̇ 〉 for an infinitely wide iron sheet,shown in Fig.12, as follows:𝜇𝜇̇ 𝑟𝑟 (𝜔𝜔) = 𝜇𝜇𝑟𝑟

tan 𝛼𝛼𝛼𝛼(27)where 𝛼𝛼 = 𝑘𝑘𝑘𝑘, 𝑘𝑘 = �−𝑗𝑗𝑗𝑗𝑗𝑗𝑗𝑗, which is obtained by solving the
one-dimensional quasi-static Maxwell equation
𝑑𝑑 2 𝐻𝐻𝑧𝑧+ 𝑘𝑘 2 𝐻𝐻𝑧𝑧 = 0
𝑑𝑑𝑥𝑥 2(28)

Because the diamagnetic constant 𝑁𝑁 nearly equals zero for theiron sheet, the Ollendorff formula (7) reduces to 〈𝜇𝜇̇ 𝑟𝑟 〉 ≈ 𝜂𝜂𝜇𝜇̇ 𝑟𝑟 .
The Cauer circuit shown in Fig.11 can also be synthesized from(27) [13]. It has been shown that the same Cauer circuit is

y
x

2d
z

Fig.12 A steel plate included in laminated steel sheets
The eddy current field can be determined by solving (23).

B [T]derived by expanding the magnetic field in (28) using the
Legendre polynomials 𝑃𝑃𝑘𝑘 as [15]𝐵𝐵𝑧𝑧 (𝑡𝑡, 𝑦𝑦) =�

𝑘𝑘=0,2,4…

𝑏𝑏𝑘𝑘 (𝑡𝑡)𝑃𝑃𝑘𝑘 (𝑦𝑦⁄𝑑𝑑 )(29)

Moreover, the time-domain analysis [16] of the laminated steel
sheets based on the expansion (29) has been shown equivalent
to the circuit analysis using the Cauer circuit [15]. It isconjectured from this fact that in general the Cauer circuits
synthesized from the complex permeability are always
equivalent to those obtained from the field expansion using
appropriate orthogonal polynomials. Furthermore, the synthesisof the Cauer circuits from quasi-static electromagnetic fields inarbitrary shaped domain is discussed in [17]. The Cauer circuit
for the arbitrary quasi-static fields can also be synthesized using
MOR, which will be described in III.D.
F. Unit-cell Approach
There is yet another homogenization method to analyzematerials composed of multiple elements such as multi-turn coiland SMC, which is the unit-cell approach [18-21]. In thisapproach, material composed of multiple elements is assumed
to consist of the same-sized unit cells, as shown in Fig.13. The
field in the unit cell is analyzed using, e.g., FEM under theboundary conditions considering the spatial periodicity. Thismethod is advantageous over the method based on the complex
permeability when the conductor included in the unit cell, suchas a wire and particle, has arbitrary shapes. The drawback in theunit-cell method is that we have to perform FE analyses fordifferent filling rates, frequencies and material constants unlike
the method based on the complex permeability. In thefrequency-domain analysis, the latter is, thus, advantageouswhen we consider simple shapes, such as circle and sphere, andalso square and cuboid which can be approximated by theformer two shapes. We can also obtain the macroscopic
complex permeability 〈𝜇𝜇̇ (𝜔𝜔)〉 for arbitrary shaped materials
based on the unit-cell approach with aid of MOR, which will bementioned in III.E.
In the time-domain analysis using the unit-cell approach, wehave to solve the one or two-dimensional quasi-static Maxwell
equation in each finite element at all the time steps. This seems
time consuming. If the system is linear, the RC method with
Padé approximation is fairly effective because no field analysis
is required. If the system is nonlinear, the Cauer-circuit
approach, which is equivalent to the expansion with the
Legendre polynomials, would be effective.

III. MODEL ORDER REDUCTION
A. Overview
In the design of electric apparatus, we often encounter thesituations in which we have to the FE equation for many times
changing material constants, shapes, configurations,frequencies and so on. In the design optimization based on FEM,
the FE equation has to be solved for different parameters. When
we solve the time-domain problems, we also repeatedly solve
the FE equation at different time steps. MOR is quite useful to
reduce the computational time in these situations.
One of the aims of MOR is to reduce the degree of freedoms(DoFs) in a given system [22]. For example, let us consider ahumanoid which includes a lot of actuators. It would be uneasy
to adequately control all the DoFs included in the humanoid ina coordinated manner. However, all the DoFs in various joints
in a human body seem not to be linearly independent. Thismeans that DoFs can be reduced without losing abilities torepresent natural motion patterns. In the proper orthogonal
decomposition (POD), which will be described in the following
subsection, important DoFs are systematically extracted toreduce the system equation.
Another aim of MOR is to construct a simple expression for
the input-output relation for a given system. When the relation
is written in a rational polynomial in the s-domain, thiscorresponds to the Padé approximation. The rational
polynomials can be obtained from POD and also via tridiagonalization of the system matrix using, e.g., the Padé via
Lanczos (PVL) method whose details will be described in III.C.
The rational polynomials can also be obtained directly from themeasured data or results of frequency sweep using, e.g., thevector fitting [23], and also by curve fitting based on
optimization algorithms such as genetic algorithm (GA) andsimulated annealing [24]. These MOR approaches, referred toas macro-modeling in microwave engineering [25], areschematically shown in Fig.14.
Analytical solution, Homogenization
PVL, PRIMA etc.
Equation

Reduced Eq.
POD

Field analysis
Snapshots

Data X =[𝒙𝒙1 , … , 𝒙𝒙𝑚𝑚 ]Freq. sampling

Eigenvalue
problem

Vector fitting
Curve fitting using GA etc.

Rational
polynomial

Equivalent
Circuits

Fig.14 Methods in model order reduction

material composed ofmultiple elements

unit cells

Fig.13 Unit-cell approach
It is assumed that the material composed of multiple elements
consists of the same-sized unit cells. The field in each unit cell isnumerically analyzed under adequate boundary conditions.
Fig.15 Reduction of DoFs by POD

B [T]

B. Proper Orthogonal Decomposition (POD)
Let us consider the linearly independent data 𝒙𝒙𝑖𝑖 ∈ ℝ𝑛𝑛 , 𝑖𝑖 =1,2 … , 𝑚𝑚 which are represented by dots in Fig. 15. Assumingthat they are obtained from 𝒙𝒙𝑖𝑖 = 𝒙𝒙′𝑖𝑖 − 𝝁𝝁, where 𝒙𝒙′𝑖𝑖 , 𝝁𝝁 are theoriginal data and their mean, the mean of 𝒙𝒙𝑖𝑖 is set to zero.
Moreover, we assume that 𝑛𝑛 ≫ 𝑚𝑚 , which would hold, forexample, when FE solutions with 𝑛𝑛 DoFs are snapshotted at 𝑚𝑚
points. If 𝒙𝒙𝑖𝑖 distribute along an axis, they can be identified byreferring their one-dimensional coordinate along the axis shownin Fig.15. The direction vector w of this coordinate axis can befound in such a way that the variation 𝜎𝜎 2 of 𝒙𝒙𝑖𝑖 projected on w
becomes maximum. This can be mathematically expressed as𝜎𝜎 2 = 𝒘𝒘t S𝒘𝒘 → max.,(a) conventional FEM

sub. to ‖𝒘𝒘‖2 = 1(30)

X = [𝒙𝒙1 , 𝒙𝒙2 , … , 𝒙𝒙𝑚𝑚 ] ∈ ℝ𝑛𝑛×𝑚𝑚(31)

𝒘𝒘𝑡𝑡 S𝒘𝒘→ max.
𝒘𝒘𝑡𝑡 𝒘𝒘(32)

where S denotes the variance-covariance matrix defined byS=

1
XX t ,
𝑚𝑚

The matrix X is called the data matrix, which is included in
Fig.14. The eigenvalue and corresponding eigenvectors of S are
here denoted by λ1 ≥ 𝜆𝜆2 ≥ ⋯ ≥ 𝜆𝜆𝑚𝑚 , 𝒘𝒘1 , 𝒘𝒘2 , … , 𝒘𝒘𝑚𝑚 . Note that
because the rank of S is 𝑚𝑚 , it has 𝑛𝑛 − 𝑚𝑚 zero eigenvalues,
which are omitted. The problem (30) is equivalent to
𝑅𝑅 =(b) POD
Fig.16 Magnetic fields in iron core for weak current
B [T]It can be shown that 𝑅𝑅 ≤ λ1 , and, 𝑅𝑅 = 𝜆𝜆1 when 𝒘𝒘 = 𝒘𝒘1 . It is,therefore, concluded that 𝒘𝒘 which satisfies (30) is just 𝒘𝒘1 .
If 𝒙𝒙𝑖𝑖 distribute within a three-dimensional ellipsoid, they
can be identified by the coordinates along the elliptic axes. The
direction vectors, orthogonal to each other, which span thethree-dimensional space, are given by 𝒘𝒘1 , 𝒘𝒘2 , 𝒘𝒘3 . Moreover,the elliptic radii along the axes are given by �𝜆𝜆𝑖𝑖 . If λ1 ≥ 𝜆𝜆2 ≫𝜆𝜆3 , then it can be said that 𝒙𝒙𝑖𝑖 almost distribute within a twodimensional ellipsoid (ellipse). In this case, 𝒙𝒙𝑖𝑖 can be
expressed by the linear combination of the two major direction
vectors. In general, if there are dominant k eigenvalues, we havethe approximation
𝒙𝒙 ≈ W𝒚𝒚,(a) conventional FEM

W ∈ ℝ𝑛𝑛×𝑘𝑘

where W = [𝒘𝒘1 , 𝒘𝒘2 , … , 𝒘𝒘𝑘𝑘 ]. The eigenvalues and eigenvectors
can be obtained by the spectral decomposition of X t X or
performing the singular value decomposition of X. Since thevector 𝒙𝒙 is expressed by the linear combination of theorthogonal bases, this method is called POD. The basis vectors
are constructed from the snapshotted data in POD, while they
are constructed by tri-diagonalization of the system matrix in
PVL as will be shown in the next sub-section.
In POD applied to FE analysis, we perform frequency ortime-domain analysis to obtain the snapshots 𝒙𝒙𝑖𝑖 of theelectromagnetic field [26, 27]. Then we compute W from X to
reduce the FE equation, K𝒙𝒙 = 𝒃𝒃, K ∈ ℂ𝑛𝑛×𝑛𝑛 , toW t KW𝒚𝒚 = W t 𝒃𝒃(b) POD
Fig.17 Magnetic fields in iron core for strong current

(33)

(34)

Now the number of unknowns is reduced from n to k. If 𝑘𝑘 ≪ 𝑛𝑛,
the computational time for solution of the FE equation could
drastically be reduced. Figures 16 and 17 show the magnetic
field distributions computed by solving the original quasi-static
FE equation and its reduced counterpart. They are in good
agreement when the magnetic induction is weak. There exist,

however, small deteriorations between them when the magnetic
induction is saturated.
When there are strong magnetic saturations in the magnetic
core, it would be difficult to express the field as the linear
combination of the form (33). Several remedies have been
proposed for this problem. In Ref.[28, 29], POD is applied only
to the linear domain which does not include saturating magnetic
cores. In Ref.[30], the strongly saturated portions are extracted
so that the field is interpolated by them in the discrete empiricalinterpolation method.
When analyzing electromagnetic fields around moving
objects, we have to perform a lot of FE computations fordifferent positions. Although the computational time could be
reduced by employing POD, we would need many basis vectors
to obtain accurate results because the fields have significant
variations for different configurations. The block-MOR basedon POD has been proposed to circumvent this difficulty [31]. In
this method, all the possible configurations are classified into
several gropes, and the different set of basis vectors isconstructed for the configurations belonging to each grope.
That is, the snapshots are subdivided into several blocks: X =[X1 , X2 , … ]. The basis matrix W𝑖𝑖 is then constructed for X𝑖𝑖 . The
reduced FE equation W𝑖𝑖t KW𝑖𝑖 𝒚𝒚 = W𝑖𝑖t 𝒃𝒃 is solved depending on
the configuration. This method has extended to apply it to theanalysis of rotating machines [32], which is shown in Fig. 18.
In the design optimization based on population-based
methods such as GA and immune algorism, the field equation
has to be solved for many times for fitness evaluations. To
reduce the computational cost, POD has been applied to those
optimizations [34]. In this method, the snapshots are taken byrandomly changing the value of design parameters to construct
the basis matrix W. The design optimization is carried out inwhich the reduced FE equation is solved. It has been shown thatthe computational time is reduced to about 10% when applying
this method to antenna-shape optimization. Moreover, POD canalso be used to synthesize the equivalent circuit from an FE
model of electric apparatus. This will be mentioned in III.D.
C. Padé via Lanczos
In PVL [35], the Padé approximation of the transfer function is
obtained on the basis of the Lanczos process which is widelyused for eigenvalue computations. This method can avoid ill
conditioning in the asymptotic waveform evaluation which

directly determines the expansion coefficients of the transfer
function. This section is devoted to formulation of PVL which
is applied to the analysis of quasi-static electromagnetic fields.
Let us consider the quasi-static electromagnetic field
coupled with electric circuits which are governed by𝜅𝜅𝜕𝜕𝑨𝑨+ rot(𝜈𝜈rot 𝑨𝑨) = 𝐼𝐼𝒋𝒋𝜕𝜕𝜕𝜕
𝑅𝑅𝑅𝑅 + 𝐿𝐿

𝑑𝑑𝑑𝑑 𝑑𝑑Φ
+= 𝑉𝑉
𝑑𝑑𝑑𝑑 𝑑𝑑𝑑𝑑

(35b)

𝑑𝑑𝒙𝒙
= 𝒍𝒍𝑉𝑉𝑑𝑑𝑑𝑑(36a)where 𝑨𝑨, 𝜈𝜈, 𝜅𝜅, 𝑅𝑅, 𝐿𝐿, 𝐼𝐼, 𝑉𝑉 are the vector potential, magnetic
reluctivity, electric conductivity, resistance, inductance, current,
voltage and Φ = ∫ 𝑨𝑨 ⋅ 𝒋𝒋𝑑𝑑𝑑𝑑 is the magnetic flux. Moreover, 𝒋𝒋
denotes the current direction vector defined by 𝒋𝒋 = J⁄I, where
𝑱𝑱 is the current density. Equation (35a) is discretized by FEM.
The resultant FE equation as well as circuit equation (35b) canbe written as a state and output equations as follows:
K𝒙𝒙 + N𝐼𝐼 = 𝒍𝒍𝑡𝑡 𝒙𝒙(36b)

𝐻𝐻(𝑠𝑠) = 𝒍𝒍𝑡𝑡 (K + 𝑠𝑠N)−𝟏𝟏 𝒍𝒍(37)

𝐻𝐻(𝑠𝑠0 + 𝜎𝜎) = 𝒍𝒍𝑡𝑡 (I − 𝜎𝜎A)−𝟏𝟏 𝒓𝒓(38)

where K, N denote the FE matrices, 𝒙𝒙 ∈ ℝ𝑛𝑛 is the unknown
vector composed of the magneto-motive forces (integration of
𝑨𝑨 along element edges) and I, and 𝒍𝒍 = [0,0, . . ,0,1]𝑡𝑡 .
The Laplace transform of (36a) leads to the transfer function
given bywhich corresponds to the admittance function. We expand (37)around an expansion point 𝑠𝑠0 to obtainwhere

A = −(K + 𝑠𝑠0 N)−𝟏𝟏 N
𝒓𝒓 = (K + 𝑠𝑠0 N)−𝟏𝟏 𝒍𝒍

Conventional POD

θ=θ1

θ=θ4

Block-MOR
θ=θ2

θ=θ1

block
ブロック11(θ1-θ2)

θ=θ3

block 2ブロック2
(θ2-θ3)

(35a)

θ=θ4

block 3ブロック3
(θ3-θ4)

Fig.18 POD approaches for motor analysis [33]

(39a)

(39b)

It can be shown that A has real eigenvalues λ1 ≥ 𝜆𝜆2 ≥ ⋯ ≥ 𝜆𝜆𝑛𝑛
and diagonalizable with its orthonormal eigenvectors U =[𝒖𝒖1 , 𝒖𝒖2 , … , 𝒖𝒖𝑛𝑛 ] : A = UΛU t , where Λ = diag[ 𝜆𝜆1 , … , 𝜆𝜆𝑛𝑛 ] . By
diagonalization of A, (38) becomes
𝐻𝐻(𝑠𝑠0 + 𝜎𝜎) = 𝒇𝒇𝑡𝑡 (I − 𝜎𝜎Λ)−𝟏𝟏 𝒈𝒈
𝑛𝑛

=�
𝑖𝑖=1

𝑓𝑓𝑖𝑖 𝑔𝑔𝑖𝑖
1 − 𝜎𝜎𝜆𝜆𝑖𝑖

(40)

⎤
⎥
⎥
⎥
⎥
𝛼𝛼𝑞𝑞 ⎦(41a)

where 𝑓𝑓𝑖𝑖 = U t 𝒍𝒍, 𝑔𝑔𝑖𝑖 = U t 𝒓𝒓. Now it is possible to easily computethe admittance from the rational function (40) at any frequency
around 𝑠𝑠0 . It needs, however, large computational cost of order
𝑂𝑂(𝑛𝑛3 ) to perform the diagonalization when the number of
unknowns n is large as usual in FE analysis.
For this reason, we employ the Lanczos process whichconstructs the tri-diagonal matrix of order 𝑞𝑞 , 𝑞𝑞 ≪ 𝑛𝑛 , whichapproximates A . Since A is unsymmetric, we use the biLanczos process, which generates the tri-diagonal matrices
𝛼𝛼
⎡ 1
⎢𝜌𝜌2
T𝑞𝑞 = ⎢ 0
⎢
⎢
⎣
𝛼𝛼
⎡ 1
⎢𝜂𝜂2
�𝑞𝑞 = ⎢ 0
T
⎢
⎢
⎣

𝛽𝛽2
𝛼𝛼2
𝜌𝜌3

0
𝛽𝛽2
𝛼𝛼3

𝛾𝛾2
𝛼𝛼2
𝜂𝜂3

0
𝛾𝛾2
𝛼𝛼3

⋱

⋱

and the basis vectors which satisfy

⎤
⎥
⎥
⎥
⎥
𝛼𝛼𝑞𝑞 ⎦

AV𝑞𝑞 = V𝑞𝑞 T𝑞𝑞 + �0,0, … , 𝒗𝒗𝑞𝑞+1 �𝜌𝜌𝑞𝑞+1

(41a)

(42a)

�𝑞𝑞 + �0,0, … , 𝒘𝒘𝑞𝑞+1 �𝜂𝜂𝑞𝑞+1
At W𝑞𝑞 = W𝑞𝑞 T(42b)

W𝑞𝑞𝑡𝑡 �AV𝑞𝑞 𝒚𝒚 − 𝜆𝜆V𝑞𝑞 𝒚𝒚� = D𝑞𝑞 �T𝑞𝑞 𝒚𝒚 − 𝜆𝜆𝒚𝒚� = 0(43)

where Vq = �𝒗𝒗1 , 𝒗𝒗2 , … , 𝒗𝒗𝑞𝑞 � , Wq = �𝒘𝒘1 , 𝒘𝒘2 , … , 𝒘𝒘𝑞𝑞 � . The
generated vectors are bi-orthogonal, that is, 𝒗𝒗𝑖𝑖 ⋅ 𝒘𝒘𝑖𝑖 = 𝑑𝑑𝑗𝑗 , 𝒗𝒗𝑖𝑖 ⋅�𝑞𝑞𝑡𝑡 D𝑞𝑞 ,𝒘𝒘𝑗𝑗 = 0, 𝑖𝑖 ≠ 𝑗𝑗. From (42), it can be proved that D𝑞𝑞 T𝑞𝑞 = T
D𝑞𝑞 = diag[𝑑𝑑1 , 𝑑𝑑2 , … , 𝑑𝑑𝑞𝑞 ].
The generated tri-diagonal matrix is an approximation of in
the following sense: we consider the eigenvalue problem A𝒙𝒙 =𝜆𝜆𝒙𝒙 . We assume that 𝒙𝒙 can be approximated by the linear
combination of 𝒗𝒗𝑖𝑖 , 𝑖𝑖 = 1,2, … , 𝑞𝑞, that is, 𝒙𝒙 = V𝑞𝑞 𝒚𝒚. Moreover,
we apply the Galerkin method to the eigenvalue problem toobtainHence by solving the small eigenvalue problem T𝑞𝑞 𝒚𝒚 = 𝜆𝜆𝒚𝒚, we
obtain q approximated eigenvalues of A. Moreover, it can beshown that the transfer function (38) is expressed in terms of
T𝑞𝑞 as follows:−𝟏𝟏
𝐻𝐻(𝑠𝑠0 + 𝜎𝜎) = 𝒍𝒍 𝒓𝒓𝒆𝒆1𝑡𝑡 �I − 𝜎𝜎T𝑞𝑞 � 𝒆𝒆1
𝑡𝑡

(44)

Diagonalization of T𝑞𝑞 in (44) leads to the partial fraction
expansion, which is equivalent to the Padé approximation. Thisis the result of PVL.

D. Synthesis of equivalent circuits using MOR
It is shown in II.D that the Cauer circuit equivalent to themacroscopic complex permeability of SMC can be derived torealize time-domain analysis. The Cauer circuit equivalent to
the FE model of electric apparatus can also be derived on thebasis of MOR. Such an equivalent circuit can be analyzed much
faster the FE equation. Moreover, it is possible to couple theequivalent circuits with the power and control circuits, and toanalyze the overall performance using circuit simulators. In
Ref.[36], the Cauer equivalent circuits for an air core coil andalso an inductor are synthesized using PVL taking the eddy
currents in the coil windings into account. Moreover, in
Ref.[37], the Cauer equivalent circuit for an induction heating(IH) device is synthesized using POD.
Here let us consider the synthesis of the Cauer circuit fromthe FE model of an IH device [38], shown in Fig.19, using PVL.
In this IH device, the alternating current, 300 Hz, carried by thecoil generates the magnetic flux which penetrates into the
conducting plate. The eddy current induced by this magnetic
flux generates Joule losses which increase the temperature of
the conducting plate. By applying PVL to the FE equation (36a)which governs the IH device, we obtain the admittance function
expressed in the rational polynomial of the form (44). As we dofor 〈𝜇𝜇̇ (𝜔𝜔)〉 in II.C, we apply the Euclidean algorithm to rewrite
the rational polynomial as a continued fraction. From this result,
we can directly synthesize the Cauer circuit shown in Fig.10,
where 𝑅𝑅0 is inserted in left of 𝐿𝐿1 , which is equivalent to the FEmodel of the IH device shown in Fig.19. Because the iron corehas magnetic saturation, we introduce the nonlinear relationshipbetween current and flux in 𝐿𝐿1 of the Cauer circuit. In thiscircuit, 𝑅𝑅0 represents the DC resistance of the winding coil,while 𝑅𝑅𝑘𝑘 , 𝐿𝐿𝑘𝑘 , 𝑘𝑘 ≥ 2 represent the resistance and inductancerelevant to the eddy current loss and fluxes generated by theeddy currents in the conducting plate. Figure 20 shows the time
responses of the coil current in the IH device obtained from thesynthetized Cauer circuit and corresponding FE equation. The
amplitude of the alternating input voltage is set to 1V and 2V
for Fig.20 (a) and (b), respectively. When the amplitude is 1V,
the current is almost sinusoidal, while it includes higher
harmonics when the amplitude is increased to 2V due to themagnetic saturation in the iron core. In spite of existence in the
waveform distortion in Fig.20 (b), the Cauer circuit in which
the magnetic nonlinearity is introduced in 𝐿𝐿1 , agrees well withthat computed by FEM. Once the Cauer equivalent circuit issynthesized, simulation and design of the IH device coupledwith the external control and power circuits can effectively be
performed.

conducting plate
Fig.19 Simple model of induction heating device

FEM

1000

Linear_Cauer circuit

Present method

800

Current [A]

600
400
200
0
-200
-400
-600

0

0.01

0.02

0.03

0.04

Time [sec.]

0.05(a) input voltage: 1V

5000

Current [A]

4000
3000
2000
1000

VI. CONCLUSIONS AND OUTLOOKIn this article, the homogenization method and MOR have beenshown effective for reduction of computational cost in FE
analysis of electromagnetic fields. The following problems
remain for future studies.
(1) There is no experimental verification for the numericalmethod for modeling of Litz wires mentioned in II.B. In
particular, the dependence of eddy current loss and coil
impedance on the twisting has to be studied by experimentand numerical computations.
(2) It is known that there exist electric and magnetic contacts
among the magnetic particles used in SMC. It is necessary
to take this effect into account for homogenization of SMC.
(3) It has been shown that the Cauer equivalent circuit can beconstructed for laminated steel sheets, SMC cores well as
stationary electric apparatus such as inductors and IH
devices. It remains unclear if this method can be extended
to modeling of rotating machines.
(4) The homogenization-based analysis of multi-turn coils and
SMC cores can be performed based on the Padé
approximation and Cauer equivalent circuit. The resulted
algorithm is, however, rather complicated. Simplification ofthese methods would be required.
V. REFERENCES

0
-1000
-2000

FEM
0

0.01

linear_Cauer circuit
0.02

0.03

Time [sec.]

Present method
0.04

0.05

(b) input voltage: 2V
Fig.20 Transient currents in windings in IH device“linear_Cauer circuit” and “Present method” represent the resultsfrom the Cauer circuit, shown in Fig.11, without and withintroducing the magnetic nonlinearity in L1.

E. Homogenization using MOR
It is shown in II.A and II.D that the macroscopic complex
permeability 〈𝜇𝜇̇ (𝜔𝜔)〉 of stranded wires and SMC can beobtained by substituting the complex permeability 𝜇𝜇̇ (𝜔𝜔) ,expressed in a closed form, into the Ollendorff formula (7). Thismethod is valid when the inclusion shape is simple. It is,however, still possible to obtain 〈𝜇𝜇̇ (𝜔𝜔)〉 for arbitrary shaped
inclusions such as ellipsoidal magnetic particles and wire withrectangular-shaped cross section using MOR. There are two
possibilities to do so. In the first method, the magnetization 𝑀𝑀
induced by eddy currents in an isolated inclusion immersed intime-harmonic magnetic field is analyzed using, e.g., FEM. The
magnetization 𝑀𝑀 , which can be formulated as a transfer
function, is expressed as a rational function of s using either
POD or PVL, which then is substituted into the Ollendorff
formula (7) to obtain 〈𝜇𝜇̇ (𝜔𝜔)〉 . In this computation, thedemagnetization constant 𝑁𝑁 of the inclusion has to be
numerically determined. This method has been proposed anddiscussed in [39]. In the second method, the magnetization of
the inclusion confined in a unit cell is analyzed using FEM andis then expressed as the rational function of s using either POD
or PVL. In this case, 〈𝜇𝜇̇ (𝜔𝜔)〉 is obtained without introducing 𝑁𝑁
although field computations have to be carried out to know its
dependence on the filling factor 𝜂𝜂 . This will be discussed
elsewhere.

[1] H. Igarashi, “Semi-analytical approach for finite element
analysis of multi-turn coil considering skin and proximity
effects,” IEEE Trans. Magn., vol.53, no.1, 7400107, 2017.
[2] L.D. Randau, L.P. Pitaevskii, E.M. Lifshitz,
Electrodynamics of Continuous Media, Second Edition,
Butterworth-Heinemann, 1979.
[3] Y. Ito, H. Igarashi, “Computation of macroscopic
electromagnetic properties of soft magnetic composite”,
IEEE Trans. Magn., 49(5), pp. 1953-1956, 2013.
[4] F. Ollendorff, “Magnetostatik der Massekerne,” Arch. f.
Eledtrotechnik., 25, pp. 436-447, 1931.
[5] S. Hiruma, H. Igarashi, “Fast three-dimensional analysisof eddy current in Litz wire using integral equation”,
presented at CEFC2016, to be published in IEEE Trans.
Magn.
[6] P.P Silvester, R.L. Ferrari, Finite elements for electricalengineers, Cambridge University Press, 1996,[7] Y. Sato, H. Igarashi, presented at Workshop “Advances inMultiscale Methods and Homogenization for Laminatesand Windings in Magnetic Fields,” 2016, to be submitted.
[8] R. Luebbers, et al., “A frequency-dependent finitedifference time-domain formulation for dispersivematerials,” IEEE Trans. Electromagnetic Compatibility,
vol.32, no. 3, pp.222-227, 1990.
[9] T. Sato, S. Aya, H. Igarashi, M. Suzuki, Y. Iwasaki, K.
Kawano, "Loss computation of soft magnetic composite
inductors based on interpolated scalar magnetic property,"
IEEE Trans. Magn., vol. 51, Issue 3, March, 2015, Art.
8001704[10] Y. Sato, H. Igarashi, “Time-domain analysis of soft
magnetic composite using equivalent circuit obtained via
homogenization”, presented at CEFC2016, submitted toIEEE Trans. Magn.
[11] H. Igarashi, et al., "Macroscopic magnetic property ofsoft magnetic composite considering eddy currents", (inJapanese) IEEJ Technical Paper, SA-14-075, RM-14-091,
2014

[12] A. Berthault, D. Rousselle, G. Zerah, “Magnetic
properties of Permalloy microparticles,”Journal ofMagnetism and Magnetic Materials, vol. 112, pp. 477-480,
1992.

[28]

D. Schmidthausler, S. Schops, M. Clemens, "Linear
subspace reduction for quasistatic field simulations toaccelerate repeated computations," IEEE Trans. Magn.
vol.50, no.2, Art. 7010304, 2014.

[13] E. J. Tarasiewicz, A. S. Morched, A. Narang, E. P. Dick,“Frequency dependent eddy current models for nonlinear
iron cores,” IEEE Trans. on Power Systems, vol.8 no.2 pp.
588-597, May, 1993.

[29]

Y. Sato, M. Clemens, H. Igarashi, "Adaptive subdomain
model order reduction with discrete empiricalinterpolation method for nonlinear magneto-quasi-static
problems," IEEE Trans. Magn., vol. 52, vol. 3, Art.
1100204, 2016.

[30]

T. Henneron, S. Clénet, Model order reduction of nonlinear magnetostatic problems based on POD and DEI
methods, IEEE Trans. Magn. vol.50, no.2, Art. 7000604,
2014

[31]

T. Sato, Y., Sato, H. Igarashi, "Model order reduction formoving objects: fast simulation of vibration energy
harvesters," COMPEL, vol. 34, no.5, 2015.

[32]

T. Shimotani, Y. Sato, T. Sato, Hajime Igarashi, "Fastfinite-element analysis of motors using block modelorder reduction," IEEE Trans. Magn., vol. 52, vol. 3, Art.
7207004, 2016.

[33]

T. Shimotani, Y. Sato, Hajime Igarashi, “Acceleration ofMotor Analysis Using Block Model Order Reduction”,(in Japanese), IEEJ Technical Report, SA-15-099, RM15-137, 2015.

[34]

Y. Sato, F. Campelo, H. Igarashi, "Fast shape
optimization of antennas using model order reduction,"
IEEE Trans. Magn. vol. 51, no. 3, 2015, Art. 7204304.

[35]

P. Feldmann and R. A. Freund, "Efficient linear circuit
analysis by Pade approximation via the Lanczos
process," IEEE Trans. Computer-Aided Design, vol. 14,
no. 5 pp. 639-649, 1995.

[36]

Y. Sato, H. Igarashi, "Generation of equivalent circuitfrom finite element model using model order reduction,"
IEEE Trans. Magn., vol. 52, no. 3, Art. 1100304, 2016.

[37]

T. Shimotani, Y. Sato, H. Igarashi, “ Direct synthesis ofequivalent circuits from reduced FE models using proper
orthogonal decomposition, ” COMPEL, vol. 35, no. 6,
pp. 2035-2044, 2016.
Y. Sato, T. Shimotani, H. Igarashi, "Synthesis of Cauerequivalent circuit based on model order reduction
considering nonlinear magnetic property", presented at
CEFC2016, submitted to IEEE Trans. Magn.

[14] T. Miyazaki, T. Mifune, T. Matsuo, Y. Shindo, Y.
Takahashi and K. Fujiwara, “Equivalent circuit modelingof dynamic hysteretic property of silicon steel under pulse
width modulation excitation,” Journal of Applied Physics,
117, 17D110, 2015.
[15] Y. Shindo, T. Miyazaki, T. Matsuo, "Cauer circuit
representation of the homogenized eddy-current field
based on the legendre expansion for a magnetic sheet,"IEEE Trans. Magn., vol.52, no.3, Art. 6300504, 2015.
[16] J. Gyselinck, R.V. Sabariego, P. Dular, “A nonlinear timedomain homogenization technique for laminated ironcores in three-dimensional finite-element models,” IEEE
Trans. Magn., vol. 42, pp. 763–766, 2006.
[17] Y. Shindo, A. Kameari, and T. Matsuo “High frequency
nonlinear modeling of magnetic sheets using polynomial
expansions for eddy-current field”, IEEJ Trans. Powerand Energy, Vol.137, No.3, to be published, 2017[18] A. Podoltsev and B. Lebedev, “Analysis of effective
resistance and eddy current losses in multiturn winding of
high-frequency magnetic components,” IEEE Trans.
Magn., vol. 39, no. 1, pp. 539–548, 2003.
[19] J. Gyselinck and P. Dular, “Frequency-domain
homogenization of bundles of wires in 2-D
magnetodynamic FE calculations,” IEEE Trans. Magn.,
vol. 41, no. 5, pp. 1416–1419, 2005.
[20] H. Waki and H. Igarashi, “Estimation of effective
permeability of magnetic composite materials”, IEEE
Trans. Magn., 41, no. 5, pp. 1520-1523, 2005.
[21] Ito and H. Igarashi, “Computation of macroscopic
electromagnetic properties of soft magnetic composite”,
IEEE Trans. Magn., 49, no.5, pp. 1953-1956, 2013.
[22]

W. H. Schilders, H. A. van der Vorst, J. Rommes (Eds.),
Model Order Reduction: Theory, Research Aspects andApplications, Springer, 2008

[23]

B. Gustavsen, A. Semlyen, "Rational approximation of
frequency domain responses by vector fitting," IEEE
Trans. Power Delivery, Vol. 14, No. 3, pp.1052-1061,
1999

[24]

T. Shimotani, Y. Sato, H. Igarashi, "Equivalent-circuit
generation from finite-element solution using proper
orthogonal decomposition," IEEE Trans. Magn., vol. 52,
vol. 3, Art. 7206804, 2016.

[25]

R. Achar, M.S. Nakhla, Simulation of high-speed
interconnects, Proceedings of the IEEE, Vol.89, no.5,
pp.693 - 728, 2002.

[26]

[27]

D. Schmidthäusler, M. Clemens, "Low-orderelectroquasistatic field simulations based on proper
orthogonal decomposition," IEEE Trans. Magn., vol.48,
no.2, pp.567-570, 2012
Y. Sato, H. Igarashi, "Model reduction of threedimensional eddy current problems based on the methodof snapshots," IEEE Trans. Magn., vol. 49, no.5,
pp.1697-1700, 2013

[38]

[39]

Y. Sato, H. Igarashi, "Homogenization method based onmodel order reduction of FE analysis of multi-turn coils",
presented at CEFC2016, submitted to IEEE Trans.
Magn.
AUTHORS NAME AND AFFILIATION

Hajime Igarashi, Yuki Sato, Shogo Fujita, Shingo Hiruma,
Yoshitsugu Otomo
Hokkaido University, Japan
Kita 14, Nishi 9, Kita-ku, Sapporo, 060-0814
igarashi@ssi.ist.hokudai.ac.jp
This review article has been published in International
COMPUMAG Society Newsletter, vol.24, no.1, pp.3-14, 2017.