POD-based reduced-order model of an
eddy-current levitation problem

arXiv:1710.08180v1 [] 23 Oct 2017

MD Rokibul Hasan, Laurent Montier, Thomas Henneron, and Ruth V. Sabariego

Abstract The accurate and efficient treatment of eddy-current problems with movement is still a challenge. Very few works applying reduced-order models are available in the literature. In this paper, we propose a proper-orthogonal-decomposition
reduced-order model to handle these kind of motional problems. A classical magnetodynamic finite element formulation based on the magnetic vector potential is used
as reference and to build up the reduced models. Two approaches are proposed. The
TEAM workshop problem 28 is chosen as a test case for validation. Results are
compared in terms of accuracy and computational cost.

1 Introduction
The finite element (FE) method is widely used and versatile for accurately modelling
electromagnetic devices accounting for eddy current effects, non-linearities, movement,... However, the FE discretization may result in a large number of unknowns,
which maybe extremely expensive in terms of computational time and memory. Furthermore, the modelling of a movement requires either remeshing or ad-hoc techniques. Without being exhaustive, it is worth mentioning: the hybrid finite-element
boundary-element (FE-BE) approaches [1], the sliding mesh techniques (rotating
machines) [2] or the mortar FE approaches [3].
Physically-based reduced models are the most popular approaches for efficiently
handling these issues. They extract physical parameters (inductances, flux linkages,...) either from simulations or measurements and construct look-up tables covMD Rokibul Hasan · Ruth V. Sabariego
KU Leuven, Dept. Electrical Engineering (ESAT), Leuven & EnergyVille, Genk, Belgium, e-mail:
rokib.hasan@esat.kuleuven.be, ruth.sabariego@esat.kuleuven.be
Laurent Montier · Thomas Henneron
Laboratoire d’Electrotechnique et d’Electronique de Puissance, Arts et Metiers ParisTech, Lille,
France, e-mail: laurent.montier@ensam.eu, thomas.henneron@univ-lille1.fr

1

2

MD Rokibul Hasan, Laurent Montier, Thomas Henneron, and Ruth V. Sabariego

ering the operating range of the device at hand [4, 5]. Future simulations are performed by simple interpolation, drastically reducing thus the computational cost.
However, these methods depend highly on the expert’s knowledge to choose and
extract the most suitable parameters.
Mathematically-based reduced-order (RO) techniques are a feasible alternative,
which are gaining interest in electromagnetism [6]. RO modelling of static coupled
system has already been implemented in [7,8]. Few RO works have addressed problems with movement (actuators, electrical machines, etc.) [9–11].
In [9], authors consider a POD-based FE-BE model electromagnetic device comprising nonlinear materials and movement. Meshing issues are avoided but the system matrix is not sparse any more, increasing considerably the cost of generating
the RO model. In [10], a magnetostatic POD-RO model of a permanent magnet
synchronous machine is studied. A locked step approach is used, so the mesh and
associated number of unknowns remains constant. A POD-based block-RO model
is proposed in [11,12], where the domain is split in linear and nonlinear regions and
the ROM is applied only to the linear part.
In this paper, we consider a POD-based FE model of a levitation problem, namely
the Team Workshop problem 28 (TWP28) [5, 13] (a conducting plate above two
concentric coils, see Fig. 1). The movement is modelled with two RO models based
on: 1) FE with automatic remeshing of the complete domain; 2) FE with constraint
remeshing, i.e., localized deformation of the mesh around the moving plate, hereafter referred to as mesh deformation. Both models are validated in the time domain
and compared in terms of computational efficiency.

az (0/1)
-0.000304

-9.04e-05

az (0/1)
0.000123

Y
ZX

-0.000304

-9.04e-05

0.000123

Y
ZX

Fig. 1: 2D axisymmetric mesh of TWP28: aluminium plate above two concentric
coils (12.8 mm clearance). Real part of the magnetic flux density. Left: automatic
remeshing of the full domain, Right: mesh deformation of sub-domain around plate
with nodes fixed at it’s boundaries (except axes).

POD-based reduced-order model of an eddy-current levitation problem

3

2 Magnetodynamic levitation model
Let us consider a bounded domain Ω = Ωc ∪ ΩcC ∈ R3 with boundary Γ . The conducting and non-conducting parts of Ω are denoted by Ωc and ΩcC , respectively.
The (modified) magnetic-vector-potential (a−) magnetodynamic formulation (weak
form of Ampère’s law) reads: find a, such that
(ν curla, curl a′ )Ω + (σ ∂t a, a′ )Ωc + hn̂ × h, a′iΓ = ( js , a′ )Ωs , ∀a′

(1)

with a′ test functions in a suitable function space; b(t) = curla(t), the magnetic
flux density; js (t) a prescribed current density and n̂ the outward unit normal vector
on Γ . Volume integrals in Ω and surface integrals on Γ of the scalar product of
their arguments are denoted by (·, ·)Ω and h·, ·iΩ . The derivative with respect to
time is denoted by ∂t . We further assume linear isotropic and time independent
materials with magnetic constitutive law, so that the magnetic field is h(t) = ν b(t)
(reluctivity ν ) and electric constitutive law, given by induced eddy current density
j(t) = σ e(t), (conductivity σ ) where, electric field e(t) = −∂t a(t). Assuming a rigid
Ωc (no deformation) and a purely translational movement (no rotation, no tilting),
the electromagnetic force appearing due to the eddy currents in Ωc can be modelled
as a global quantity with only one component (vertical to the plate). If Ωc is nonmagnetic, Lorentz force can be used:
Fem (t) =

Z

Ωc

j(t) × b(t) dΩc =

Z

Ωc

−σ ∂t a(t) × curla(t) dΩc .

(2)

The 1D mechanical equation governing the above described levitation problem
reads:
m ∂t v(t) + ξ v(t) + ky(t) + mg = Fem (t)
(3)
where unknown y(t) is the center position of the moving body in the vertical direction, v(t) = ∂t y(t) is the velocity of the moving body, m is the mass of the moving
body, g is the acceleration of gravity, ξ is the scalar viscous friction coefficient,
k is the elastic constant. We apply the backward Euler method to solve (3). The
moving body displacement of system (3) results from the ensuing electromagnetic
force generated by system (1) and thus affects the geometry. Given that, the dynamics of the mechanical equation is much slower than the electromagnetic equation,
if the time-step is taken sufficiently small, one can decouple the equations. Under
this condition, the electromagnetic and mechanical equations can be solved alternatively rather than simultaneously by the weak electromechanical coupling algorithm
of [14]. We adopt this approach.

4

MD Rokibul Hasan, Laurent Montier, Thomas Henneron, and Ruth V. Sabariego

3 POD-based model order reduction
The proper orthogonal decomposition (POD) is applied to reduce the matrix system
resulting from the FE discretisation of (1):
A∂t x(t) + Bx(t) = C(t) .

(4)

where x(t) ∈ RN×1 is the time-dependent column vector of N unknowns, A, B ∈
RN×N are the matrices of coefficients and C(t) ∈ RN×1 is the source column vector.
Furthermore, the system (4) is discretized in time by means of the backward Euler
scheme. A system of algebraic equation is obtained for each time step from tk−1 to
tk = tk−1 + ∆ t, ∆ t the step size. The discretized system reads:
[A∆ t + B]xk = A∆ t xk−1 + Ck

(5)

with A∆ t = ∆At , xk = x(tk ) the solution at instant tk , xk−1 = x(tk−1 ) the solution at
instant tk−1 , Ck the right-hand side at instant tk .
In RO techniques, the solution vector x(t) is approximated by a vector xr (t) ∈
M×1
R
within a reduced subspace spanned by Ψ ∈ RN×M , M ≪ N,
x(t) ≈ Ψ xr (t) ,

(6)

with Ψ an orthonormal projection operator generated from the time-domain full
solution x(t) via snapshot techniques [15].
Let us consider the snapshot matrix, S = [x1 , x2 , . . . , xM ] ∈ RN×M from the set
of solution xk for the selected number of time steps. Applying the singular value
decomposition (SVD) to S as,
(7)
S = U ΣV T .
where Σ contains the singular values, ordered as σ1 > σ2 > . . . > 0. We consider
Ψ = U r ∈ RN×r , that corresponds to the truncation (r first columns, which has
larger singular values than a pre-defined error tolerance ε ) with orthogonal matrices
U ∈ RN×r and V ∈ RM×r . Therefore, the RO system of (5) reads
[A∆r t + Br ] xrk = Ar∆ t xrk−1 + Ckr ,

(8)

with Ar∆ t = Ψ T A∆ t Ψ , Br = Ψ T BΨ and Cr = Ψ T C [16].

3.1 Application to an electro-mechanical problem with movement
3.1.1 RO modelling with automatic remeshing technique
In case of automatic remeshing, we transfer results from the source meshk−1 to the
new target meshk by means of a Galerkin projection, which is optimal in the L2 -

POD-based reduced-order model of an eddy-current levitation problem

5

norm sense [17]. Note that, this projection is limited to the conducting domain, i.e.
the plate, as it is only there that we need to compute the time derivative. The number
of unknowns per time step tk varies and the construction of the snapshot matrix S is
not straightforward. As the solution at tk is supported on its own mesh, the snapshot
vectors xk have a different size. They have to be projected to a common basis using
a simple linear interpolation technique before being assembled in S and getting the
projection operator Ψ . The procedure becomes thus extremely inefficient.

3.1.2 RO modelling with mesh deformation technique
The automatic remeshing task is replaced by a mesh deformation technique, limited
to a region around the moving body (see, e.g., the box in Fig. 2). Therefore, in this
case, the remeshing is done by deforming the initial mesh, which is generated with
the conducting plate placed at, e.g., y0 (avoiding bad quality elements), see Fig. 2.
The mesh elements only inside the sub-domain can be deformed (shrink/expand),
see Fig. 3 and the nodes at the boundary of the sub-domain are fixed. The surrounding mesh does not vary. In our test case, we assume a vertical force (neglect the other
two components) in (2), therefore, the mesh elements only deform in the vertical direction and the nodes are fixed at the boundary of the sub-domain (not at the axes
due to the axisymmetry). The size of the sub-domain (a×b) is determined by the
extreme positions of the moving body. In our validation example, the minimum position (3.8 mm) is given by the upper borders of the coils and the maximum position
(22.3 mm) could be estimated by means of a circuital model, e.g. [5]. The number
of unknowns per time step remains now constant so the construction of matrix S is
direct.
r
ymax

b

y0
ymin
a

Fig. 2: Sub-domain for deformation: plate position at y0 = 12.8 mm (initial mesh).

r
ymax
y
b

ymin
a

Fig. 3: Sub-domain for deformation: plate position at y = 20 mm. Mesh elements
under the plate are expanded and above the plate are shrinked.

6

MD Rokibul Hasan, Laurent Montier, Thomas Henneron, and Ruth V. Sabariego

Algorithm 1: Automatic remeshing
Input : snapshot vectors
{Sc } ← {xk } ∈ Rn(k)×1
time steps
{tk }, k ∈ [1, ..., K]
A∆ t , B, C, tolerance ε
m ≤ n(k) snapshot vectors
Output: displacement yk
1

2

3
4
5
6

7
8

9

10

11
12

13
14
15
16

y0 = initial position, ∆ y0 = 0
//Time resolution
for k ← 1 to K do
//Magnetics
generate matrices A∆ tk , Bk , Ck
find length of Ck ∈ Rn(k)×1
S p = 0 ∈ Rn(k)×m
S p ← projection of {Sc } to n(k)
rank subspace
SVD of S p = U Σ V T
Ψk = U (:, 1 . . . r) with r such
that σ (i)/σ (1) > ε , ∀i ∈ [1 . . . r]
Ar∆ tk =ΨkT A∆ tk Ψk ,
Brk
Ckr

=ΨkT BkΨk ,
=ΨkT Ck

solve


A∆r tk + Brk xrk = Ckr + Ar∆ tk xrk−1
xk ≈ Ψk xrk

Algorithm 2: Mesh deformation
Input : snapshot matrix
S = [x1 , . . . , xm ] ∈ Rn×m ,
xk ∈ Rn×1
time steps
{tk }, k ∈ [1, ..., K]
A∆ t , B, C, tolerance ε
m ≤ n snapshot vectors
Output: displacement yk
1
2
3
4

5

6

7

8

9
10

compute force Fk
//Mechanics
compute displacement yk
update ∆ yk = yk − yk−1
remesh with yk
end

11
12
13
14

y0 = initial position, ∆ y0 = 0
get initial mesh
SVD of S = U Σ V T
Ψk = U (:, 1 . . . r) with r such that
σ (i)/σ (1) > ε , ∀i ∈ [1 . . . r]
//Time resolution
for k ← 1 to K do
//Magnetics
generate matrices A∆ tk , Bk , Ck
Ar∆ tk =Ψ T A∆ tk Ψ ,
Brk =Ψ T BkΨ ,
Ckr =Ψ T Ck
solve


A∆r tk + Brk xrk = Ckr + Ar∆ tk xrk−1

xk ≈ Ψ xrk
compute force Fk

//Mechanics
compute displacement yk
update ∆ yk = yk − yk−1
deform mesh with yk
end

4 Application example
We consider TWP28: an electrodynamic levitation device consisting of a conducting
cylindrical aluminium plate (σ = 3.47 · 107 S/m, m = 0.107 Kg, ξ = 1) above two
coaxial exciting coils. The inner and outer coils have 960 and 576 turns respectively.
Note that, if we neglect the elastic force, the equilibrium is reached when the Fem is

POD-based reduced-order model of an eddy-current levitation problem

7

1N. At t = 0, the plate rests above the coils at a distance of 3.8 mm. For t ≥ 0, a timevarying sinusoidal current (20 A, f = 50 Hz) is imposed, same amplitude, opposite
directions [13]. Assuming a translational movement (no rotation and tilting) we can
use an axisymmetric model. A FE model is generated as reference and origin of
the RO models. We have time-stepped 50 periods (100 time steps per period and
step size 0.2 ms), discretization that ensures accuracy and avoids degenerated mesh
elements during deformation.

4.1 RO modelling with automatic remeshing full domain
In case of full domain remeshing, the first 1500 time steps (300 ms) of the simulation, that correspond to the first two peaks (2P) in Fig. 4, are included in the snapshot
matrix.
Three POD-based RO models are constructed based on the r number of first
singular value modes greater than a prescribed error tolerance ε , that is set manually
observing the singular values decay curve of the snapshot matrix (see in Table 1).
The smaller the prescribed ε , the bigger the size of the RO model will be (size of
RO3>RO2>RO1).
Table 1: L2 -relative errors of RO models on levitation height for 2P (automatic
remeshing).
RO models
RO1
RO2
RO3

M
ε
1085 10−6
1403 10−11
1411 10−15

rel. error
1.25 · 10−1
1.03 · 10−2
2.45 · 10−6

The displacement and relative error of the full and RO models are shown in Fig. 4.
Accurate results have been achieved with the truncated basis models: RO2 and RO3,
with fix size per time step M = 1403 and 1411. This approach is completely inefficient, as the maximum number of unknowns we have in the full model is 1552.

4.2 RO modelling with mesh deformation of a sub-domain
The choice of the sub-domain to deform the mesh is a non-trivial task: it should
be as small as possible while ensuring a high accuracy. From our reference FE
solution [13], by observing the minimum and maximum levitation height of the
plate, we fixed the sub-domain size along the y−axis between ymin = 1.3 mm and
ymax = 29.3 mm, distances measured from the upper border of the coils. The size
along the x−axis has a minimum equal to the radius of the plate, i.e. r = 65 mm.
This value is however not enough due to fringing effects. We have taken different

MD Rokibul Hasan, Laurent Montier, Thomas Henneron, and Ruth V. Sabariego
levitation (mm)

8

20

Full
RO1
RO2
RO3

15
10
5
−1

relative error

10

−3

10

−6

10

−9

10

0

100

200

300

400

500

time (ms)

Fig. 4: Displacement (up) and relative error (down) between full and RO models.

size along the x−axis: 1.5r, 2r, 3r (97.5, 130 and 195 mm), measured from the axis
(Fig. 2). The meshed boxes yield 1921, 1836 and 1780 number of unknowns.
The relative errors in time shown in Table 2 decrease with the increasing subdomain lengths/box sizes considered. We have therefore chosen to further analyse
the RO results obtained with a box length along x of 195 mm (3r). The discretization
is kept constant for all RO models computation.
Table 2: L2 -relative errors of RO models on levitation height for 1P (mesh deform).
sub-domain lengths (mm) M = 7
M = 35
97.5
8.24 · 10−2 6.14 · 10−4
130
5.71 · 10−2 1.90 · 10−4
195
4.53 · 10−3 3.73 · 10−5

The first 800 time steps (160 ms) of the simulation, that correspond to the first
peak (1P), are taken in the snapshot matrix in order to generate the projection basis Ψ . In the snapshot matrix, the most important time step solutions are included,
which found as optimum selection for approximating the full solution. Then the basis is truncated as Ψ = U r (r first columns) by means of prescribed error tolerance
(ε = 10−5, 10−8 ). The basis are truncated for 1P to get RO models of size M = 7
and 35.
From Fig. 5, it can be observed that, RO model already shows very good argument with only M = 7 truncated basis, which is generated from the snapshot matrix that incorporates first peak (1P). The accuracy of RO models does not improve
significantly with the addition of following transient peaks (2P) into the snapshot
matrix, but the accuracy certainly improves with M. Hence, with M = 35 the full
and RO curves are indistinguishable. The accuracy of the RO models can also be
observed from the L2 -relative errors figure.

POD-based reduced-order model of an eddy-current levitation problem

levitation (mm)

20

9

Full
RO, M=7 1P
RO, M=7 2P
RO, M=35 1P
RO, M=35 2P

15
10
5
−2

relative error

10

−4

10

−6

10

−8

10

0

100

200

300

400

500

time (ms)

Fig. 5: Displacement (up) and relative error (down) between full and RO models for
195 mm sub-domain length.

With regard to the computation time (5000 time steps), the RO with M = 7, can
be solved less than an hour, which is 3.5 times faster than the full-domain automatic
remeshing approach, where the major time consuming part is to project the Ψ on
a same dimensional basis as the system coefficient matrices, to reduce the system
in each time step. Be aware that the computation is not optimized, performed on a
laptop, (Intel Core i7-4600U CPU at 2.10 GHz) without any parallelization.

5 Conclusion
In this paper, we have proposed two approaches for POD-based RO models to treat
a magnetodynamic levitation problem: automatic remeshing and mesh deformation
of a sub-domain around a moving body. The RO model is completely inefficient
with automatic remeshing technique, as the computational cost is nearly expensive
as the classical approach. The approach with sub-domain deformation to limit the
influence of the movement on the RO model construction has proved accurate and
efficient (low computational cost). We have shown results for three different subdomain sizes, the bigger the sub-domain the higher the accuracy. Further, computationally efficient RO modelling of such parametric model is ongoing research.

References
1. Sabariego, R., Gyselinck, J., Dular, P., Geuzaine, C., Legros W.: Fast multipole acceleration
of the hybrid finite-element/boundary-element analysis of 3-D eddy-current problems. IEEE

10

MD Rokibul Hasan, Laurent Montier, Thomas Henneron, and Ruth V. Sabariego

Trans. Magn. vol. 40, no. 2, pp. 1278–1281, (2004).
2. Boualem, B., Piriou, F.: Numerical models for rotor cage induction machines using finite
element method. IEEE Trans. on Magn., vol. 34, no. 5, pp. 3202–3205, (1998).
3. Rapetti, F.: An overlapping mortar element approach to coupled magneto-mechanical problems. Math. and Comput. in Sim., vol. 88, no. 8, pp. 1647–1656, (2010).
4. Liu, Z., Liu, S., Mohammed, O. A.: A Practical Method for Building the FE-Based Phase
Variable Model of Single Phase Transformers for Dynamic Simulations. IEEE Trans. Magn.,
vol. 43, no. 4, pp. 1761–1764, (2007).
5. Lee, S. M., Lee, S. H., Choi, H. S., Park, I. H.: Reduced Modeling of Eddy Current-Driven
Electromechanical System Using Conductor Segmentation and Circuit Parameters Extracted
by FEA. IEEE Trans. Mag., vol. 41, no. 5, pp. 1448–1451, (2005).
6. Schilders, W.H., Van der Vorst, H.A. and Rommes, J.: Model order reduction: theory, research
aspects and applications. Springer-Verlag, (2008).
7. Yue, Y., Feng, L., Meuris, P., Schoenmaker, W., Benner, P.; Application of Krylov-type Parametric Model Order Reduction in Efficient Uncertainty Quantification of Electro-thermal Circuit Models. In Proc. of the Prog. In Electromag. Res. Symp. (PIERS), pp. 379–384, (2015).
8. Banagaaya, N., Feng, L., Meuris, P., Schoenmaker, W., Benner, P.; Model order reduction of
an electro-thermal package model. IFAC-PapersOnLine, vol. 48, no. 1, pp. 934–935, (2015).
9. Albunni, M.N., Rischmuller, V., Fritzsche, T., Lohmann, B.: Model-order reduction of moving
nonlinear electromagnetic devices. IEEE Trans. Mag., vol. 44, no. 7, pp. 1822–1829, (2008).
10. Henneron, T., Clénet, S., Model order reduction applied to the numerical study of electrical
motor based on POD method taking into account rotation movement. Int. J. Numer. Model.,
vol. 27, no. 3, pp. 485–494, (2014).
11. Sato, T., Sato, Y., Igarashi, H.: Model order reduction for moving objects: fast simulation of
vibration energy harvesters. COMPEL, vol. 34, no. 5 pp. 1623–1636, (2015).
12. Schmidthäusler, D., Schöps and Clemens, M.: Linear subspace reduction for quasistatic field
simulations to accelerate repeated computations. IEEE Trans. Magn., vol. 50, no. 2, article
# 7010304, (2014).
13. Karl, H., Fetzer, J., Kurz, S., Lehner, G., Rucker, W. M.: Description of TEAM workshop
problem 28: An electrodynamic levitation device. In Proc. of the TEAM Workshop, Graz,
Austria, pp. 48–51, (1997).
14. Henrotte, F., Nicolet, A., Hedia, H., Genon, A., Legros, W.; Modelling of electromechanical
relays taking into account movement and electric circuits. IEEE Trans. Magn. vol. 30, no. 5,
pp. 3236–3239, (1994).
15. Sato, Y., Igarashi, H.: Model reduction of three-dimensional eddy current problems based on
the method of snapshots. IEEE Trans. Magn., vol. 49, no. 5, pp. 1697–1700, (2013).
16. Hasan, MD R., Sabariego, R. V., Geuzaine, C., Paquay, Y.; Proper orthogonal decomposition
versus Krylov subspace methods in reduced-order energy-converter models. In Proc. of IEEE
Inter. Ener. Conf., pp. 1–6, (2016).
17. Parent, G., Dular, P., Ducreux, J.P. and Piriou, F.; Using a Galerkin projection method for
coupled problems. IEEE Trans. Magn., vol. 44, no. 6, pp. 830–833, (2008).

