# Simulations in a Digital Twin of an Electrical Machine

Karim Cherifi  Volker Mehrmann  Philipp Schulze

###### Abstract

Digital twins have become popular for their ability to monitor and optimize a
process or a machine during its lifetime using simulations and sensor data. In
this paper, we focus on the challenge of the implementation of accurate and
real time simulations for digital twins in the context of electrical machines.

In general, this involves not only computational models for the
electromagnetic aspects, but also mechanical and thermal effects need to be
taken into account. We address mathematical tools that can be employed to
carry out the required simulations based on physical laws as well as surrogate
or data-driven models. One of those tools is a model hierarchy of very fine to
very course models as well model reduction which is required for obtaining
real-time simulations.

We discuss in detail the coupling of electromagnetic, mechanical, and thermal
models of an electrical machine to obtain a simulation model which is able to
describe the interaction of those different physical components. In this
context, a very promising setting is provided by energy-based formulations
within the port-Hamiltonian framework, which has received much attention in
the past years, especially in the context of multiphysics modeling. We present
such port-Hamiltonian formulations for the considered electromagnetic,
mechanical, and thermal models as well as for the coupled overall system.

Keywords— Digital twins, electrical machines, systems modeling, port-
Hamiltonian systems, partial differential equations.

##  1 Introduction: Numerical Simulation and Requirements for a Digital Twin

When an electrical machine is equipped with a large number of sensors, then it
produces a lot of data that can be collected, stored, and used [IMR21, Mel18].
A digital twin is a tool that uses this data combined with numerical
simulation models to accompany its physical twin in an optimal way throughout
its life cycle [RSK20, BVC20, FK21]. By means of the phenomenal advancements
in both hardware and software in recent years in terms of data processing,
computation power, and data-driven modeling methods, digital twins are
becoming a reality that can be useful in many ways. This includes the
following functions: storage of data history and feedback to design,
documentation and design optimization, condition monitoring, and predictive
maintenance.

A digital twin can be used to continuously monitor the behavior of the machine
and detect any deviations in the operating conditions of individual components
or of the overall system. In addition, a digital twin can improve predictive
maintenance by continuously predicting failures of system components using
detailed life-time simulations of the critical components in the machine. This
simulation should be adapted to the operating and usage conditions that may
affect the original life-time model. This can be done using data-driven
methods and data assimilation. A successful detection and warning in the event
of abnormal behavior allows the replacement of faulty components before their
breakdown and can prevent damage to other components. Some of these functions
need accurate simulation models that can replicate and predict he operation of
the physical twin in real time. The simulation models can also be used to
optimize the machine parameters for the respective operating environment based
on the knowledge gained from operation. For example, motor-specific monitoring
of the efficiency under parameter variations allows to identify optimized
operating points online while the machine is running or offline before turning
on the machine based on the system parameters and stored data.

In this context, high fidelity simulations that can be obtained in real time
are required by some functions of the digital twin. However, in most
conventional simulation software, it is not possible to obtain a simulation in
real time, thus the models have to simplified. For security reasons, these
simplifications and resulting decisions have to be understood, justified, and
documented. This is one of the main challenges of digital twins as discussed
in [RSK20]. One example of a mathematically sound simplification method that
can be used is model order reduction (MOR), see for instance [HHW18] for a
corresponding discussion in the context of digital twins and [BGTQ+21a,
BGTQ+21b, BGTQ+21c] for a general overview. One may also rely on data-driven
surrogate models obtained from measurement data [FFLMM22].

Altogether, having models of different scales and accuracy allows us to build
a catalogue of models including a hierarchy from complex to simple (and
faster) models as illustrated in Figure 1.

The physical system is typically described by a partial differential equation
(PDE) together with boundary conditions and may include some algebraic
constraints, arising from interconnection and interface conditions or balance
laws.

The PDE is then discretized in space using a grid hierarchy of coarse to very
fine grids that allow for accuracy adjustment. Space discretization results in
a system of ordinary differential equations (ODEs) or differential-algebraic
equations (DAEs) if we add the constraints.

The resulting dynamical model can be further simplified using reduced order
modeling. In this way, the accuracy of the resulting simulation can be
adjusted based on a user-defined tolerance used in the model reduction
algorithm. Then, depending on the application, one can use a suitable model
based on given accuracy and computation time requirements.

On the other hand, one can describe the dynamics of the physical system by
means of a data-driven realization using measurements from the physical
system. A data-driven realization can also be used to model a predefined
lumped parameter model. Alternatively, the lumped parameters can be estimated
from the discretized PDE. However, these methods typically do not provide any
guarantees about the estimation of the error.

Figure 1: Catalogue of simulation models including a system hierarchy.

In addition, we also require the simulation to be adaptive to the electrical
machine specifications and parameters. We should be able to have a simulation
setup where a new simulation can be obtained if some electrical machine part
is changed or the electrical machine itself is changed in a factory for
example. This is particularly true when dealing with high power electric
motors that can have power up to 100 MW [RE20]. These motors are generally
custom-made and thus models have to be adapted to each new motor produced. In
addition, one has to consider effects that may affect its lifetime. For
example, the vibrations that result from the operation of these huge machines
can end up damaging the machine and shorten its lifetime [TD99, SG04].

The simulation models that we consider in the construction of a digital twin
for an electrical machine can be divided into three main model families:
electrical models, mechanical models and thermal models. In this paper we
describe the basic equations describing the three model families and how these
can be incorporated in a digital twin practically, cf. sections 2–4. In
addition, we identify the main challenges facing the implementation of a
digital twin for electrical machines. Since all of these simulations belong to
the same physical twin, they are naturally coupled. This alone is a real
challenge for most simulation technologies, since one needs a multi-physics
simulation framework to obtain an accurate model of the electrical machine.

A recently very successful systems theoretic framework to deal with
multiphysics system that couple different pyhsical phenomena is that of energy
based modeling via port-Hamiltonian systems. We will employ this approach and
present in section 5 a fully coupled mathematical model accounting for
electromagnetism, mechanics, and thermodynamics in a corresponding port-
Hamiltonian formulation.

##### Mathematical Notation

Throughout this manuscript, we use bold font for vectors and vector-valued
functions. We denote the set of real numbers by $\mathbb{R}$, the standard
Euclidean space of dimension $n$ by $\mathbb{R}^{n}$, and the Euclidean inner
product of two vectors $\boldsymbol{v},\boldsymbol{w}\in\mathbb{R}^{n}$ by
$\boldsymbol{v}\cdot\boldsymbol{w}$. Furthermore, we use $\mathbb{R}^{m,n}$
for the space of $m\times n$ matrices and
$\mathbb{R}^{n_{1},n_{2},n_{3},n_{4}}$ for fourth-order tensors of size
$n_{1}\times n_{2}\times n_{3}\times n_{4}$. The tensor product between two
matrices $A\in\mathbb{R}^{m,n}$ and $B\in\mathbb{R}^{p,q}$ is denoted by
$A\otimes B\in\mathbb{R}^{m,n,p,q}$ and defined via $[A\otimes
B]_{ijkl}\vcentcolon=[A]_{ij}[B]_{k\ell}$ for $i=1,\ldots,m$, $j=1,\ldots,n$,
$k=1,\ldots,p$, and $\ell=1,\ldots,q$. Finally, for the Frobenius inner
product of two matrices $A,B\in\mathbb{R}^{m,n}$ we use the notation
$\left\langle
A,B\right\rangle_{\mathrm{F}}\vcentcolon=\mathop{\mathrm{tr}}\left(A^{\top}B\right)$,
where $\mathop{\mathrm{tr}}\left(\cdot\right)$ denotes the trace.

##  2 Electrical Simulation

In digital twins, we require a hierarchy of models where the appropriate model
is chosen based on a trade-off between simulation accuracy and computation
time. In the electromagnetic simulation, we start with the most fundamental
equations based on Maxwell equations up to the simplest surrogate models based
on equivalent circuits of the electrical machine. Although most of the
following equations hold for all rotating machines, there are still some
significant differences between synchronous and asynchronous induction
machines. In the following, as a use case we focus on induction motors when
appropriate.

###  2.1 Fundamental Equations and Finite Element Simulation

In order to study the transient time domain behavior of rotating machines, we
consider the Maxwell equations for homogeneous isotropic materials described
by

$\displaystyle\nabla\cdot\boldsymbol{E}$
$\displaystyle=\frac{\rho_{\mathrm{c}}}{{\epsilon}},$ (1a)
$\displaystyle\nabla\cdot\boldsymbol{B}$ $\displaystyle=0,$ (1b)
$\displaystyle\nabla\times\boldsymbol{E}$
$\displaystyle=-\frac{\partial\boldsymbol{B}}{\partial t},$ (1c)
$\displaystyle\nabla\times\boldsymbol{B}$
$\displaystyle=\mu(\boldsymbol{J}+{\epsilon}\frac{\partial\boldsymbol{E}}{\partial
t}),$ (1d)

where $\boldsymbol{E}$ is the electric field, $\boldsymbol{B}$ the magnetic
flux density, $\boldsymbol{J}$ the total electric current density, and
$\rho_{\mathrm{c}}$ the total electric charge density. Furthermore, $\epsilon$
and $\mu$ are the permittivity and the permeability, respectively. These
depend on the material and in free space they just reduce to the vacuum
permittivity ${\epsilon}_{0}$ and the vacuum permeability ${\mu}_{0}$,
respectively. We emphasize that (1) does not take into account any effects
resulting from the motion of the material, i.e., the coupling to the mechanics
via velocity-dependent effects is neglected in this section. An overall
coupled model which also accounts for such effects is presented in section
5.2.

Although the Maxwell equations can be applied in various ways in order to
study dynamics of rotating machines, we present below how they are most
commonly used and specifically how they are employed in commercial software
such as Ansys Maxwell. One example is the configuration shown in figure 2.

Figure 2: Example 2D configuration of an induction machine in Ansys software

The following equations focus on the 2D case and are derived mainly from the
references [Sal95, IMR21, Bia05]. The 3D case is not discussed but most of the
concepts discussed in what follows are valid for the 3D case and can be
generalized from the 2D case.

First, we focus on the electrical side. The solution for the electric
potential $\Phi(x_{1},x_{2})$ on a plane represented by $x_{1}$ and $x_{2}$
(2D case) is derived from the Maxwell equations resulting in the equation

$\nabla\cdot({\epsilon}\nabla\Phi)=-\rho_{\mathrm{c}}.$ (2)

This equation is solved using finite element methods (FEM) for a scalar
potential function $\Phi$. For more details refer to [Sal95].

Given the electric potential, one can compute the static field
$\boldsymbol{E}$, the current density $\boldsymbol{J}$, the electric flux
density $\boldsymbol{D}$ and electric flux linkage $\psi_{\mathrm{e}}$ in a
surface $S$ as

$\displaystyle\boldsymbol{E}$ $\displaystyle=-\nabla\Phi,$ (3a)
$\displaystyle\boldsymbol{J}$ $\displaystyle={\sigma}\boldsymbol{E},$ (3b)
$\displaystyle\boldsymbol{D}$ $\displaystyle={\epsilon}\boldsymbol{E},$ (3c)
$\displaystyle\psi_{\mathrm{e}}$
$\displaystyle=\int_{S}{\boldsymbol{E}\cdot\mathrm{d}\boldsymbol{S}},$ (3d)

where $\sigma$ is the conductivity specific to the material. The current is
then simply the surface integral of the electric current density. From the
electric field, the generated mechanical power and resulting torque are
derived, see for instance [CM84].

Another quantity of interest is the magnetic vector potential $\boldsymbol{A}$
that is used to compute the magnetic flux density $\boldsymbol{B}$, the
magnetic field $\boldsymbol{H}$ (for linear materials), and the magnetic flux
linkage $\psi_{i}$ in each coil as presented in

$\displaystyle\boldsymbol{B}$ $\displaystyle=\nabla\times\boldsymbol{A},$ (4a)
$\displaystyle\boldsymbol{H}$ $\displaystyle=\frac{\boldsymbol{B}}{\mu},$ (4b)
$\displaystyle\psi_{i}$
$\displaystyle=L_{ii}I_{i}+\sum\limits_{j=1}^{n}{L_{ij}I_{j}}+\psi_{\mathrm{m}},$
(4c)

where $L_{ii}$ is the self inductance in a single coil, $L_{ij}$ is the mutual
inductance between the coils and $\psi_{\mathrm{m}}$ is the permanent magnetic
flux linkage. The derivation of the magnetic flux linkage is discussed in
detail in [MHCG16] where the infinite dimensional equation is discretized and
solved efficiently using model order reduction. The reduction of the size of
the model resulting from the discretization of equations described above is
mandatory in our context of a digital twin to achieve real time computations
and real time adaptability as required by the digital twin. The work done in
[MHCG16] represents only one example how to increase the speed of the
simulation.

In the static or quasi-stationary case, the magnetic vector potential
$\boldsymbol{A}$ which is the quantity of interest here, can be derived in the
2D case from the Maxwell equations resulting in the diffusion equation

$\boldsymbol{J}={\nabla}\times(\frac{1}{\mu}{\nabla}\times\boldsymbol{A}),$
(5)

where $\boldsymbol{A}$ is the magnetic vector potential and $\boldsymbol{J}$
is the current density field. Here $\boldsymbol{A}$ and $\boldsymbol{J}$ are
assumed to be $x_{3}$-directed and independent of $x_{3}$ as they are a result
of the $x_{1}$\- and $x_{2}$-directions [Sal95]. Given an excitation
$\boldsymbol{J}$, the FEM solver finds the solution of (5) for
$\boldsymbol{A}$. Typically $\boldsymbol{J}$ is considered to be having three
parts: one due to the applied source, another due to the induced electric
field produced by time varying magnetic flux, and the third due to motion-
induced voltage. This introduces a time varying term for $\boldsymbol{A}$ and
describes the problem of the eddy current effects. Eddy currents are currents
induced within conductors by a changing magnetic field in the conductor. The
eddy current will itself create a magnetic field that cause energy loss in the
rotating machine. In the case of induction machines, eddy currents are
responsible for producing the main torque. Eddy currents and corresponding
equations describe the coupling between the electrical and the magnetic field.

In order to compute the eddy currents, we start with rewriting the Maxwell
equation (1c) using equation (4a) as

$\nabla\times\boldsymbol{E}=-\frac{\partial\nabla\times\boldsymbol{A}}{\partial
t},\\\ $ (6)

resulting in

$\boldsymbol{E}=-\frac{\partial\boldsymbol{A}}{\partial t}-\nabla\Phi.\\\ $
(7)

This reformulation results in an equation that describes the coupling between
electric and magnetic fields and their mutual effects based on (3b), (5), and
(7) as

$\nabla\times(\frac{1}{\mu}\nabla\times\boldsymbol{A})=-\sigma\frac{\partial\boldsymbol{A}}{\partial
t}-\sigma\nabla\Phi.\\\ $ (8)

This is discussed in more details in [Szü01]. In the 2D case, only the
$x_{1}x_{2}$-plane is considered taking into account the symmetry of the
machine. In the electrical engineering literature, the $x_{1}x_{2}$-plane is
commonly referred to as the $xy$-plane while $x_{3}$ is referred to as $z$ as
shown in Figure 3. As shown in [Bia05], due to the symmetry assumption, the
flux density and magnetic field strength vectors have only non-zero components
in $x_{1}$\- and $x_{2}$-direction, i.e., $\boldsymbol{B}=(B_{1},B_{2},0)$ and
$\boldsymbol{H}=(H_{1},H_{2},0)$. The resulting magnetic vector potential and
the current density vectors have an $x_{3}$-axis component only, which are
$\boldsymbol{A}=(0,0,A_{3})$ and $\boldsymbol{J}=(0,0,J_{3})$ respectively.

When applying equations (5) and (8) in the context of a rotating electrical
machine, one has to consider boundary conditions. This is discussed
extensively in [Bia05]. In the special case of asynchronous induction motors,
the geometrical and magnetic symmetry allows to limit the analysis domain to a
portion of the motor section as shown in figure 3.

Figure 3: Example of an Ansys Maxwell simulation exploiting symmetry of an
induction machine

This portion is primarily determined by the stator and rotor slot numbers, as
well as the stator winding. Another important parameter is the number of
magnetic pole pairs considered. These are are periodically arranged, so
analysis can be carried out only on one pole pair.

When a stator symmetry exists in its geommetry for example in the poles of the
machine, i.e. the number of slots is a multiple of the number of poles, the
study can be carried out only on one single pole. There are usually no
problems with geometrical symmetry, especially if the slot number is a
multiple of the pole number. On the other hand, due to the stator winding,
special attention must be paid to magnetic symmetry. The analysis domain can
be condensed by imposing appropriate periodic boundary conditions. In
addition, special operating conditions allow Dirichlet or Neumann boundary
conditions to be assigned [Bia05]. Nevertheless, we consider here the general
case in normal operating condition.

A Dirichlet boundary condition $A_{3}=0$ is assigned to the stator’s external
circumference, causing the flux lines to remain confined within the stator
yoke. However, in case of high saturation, this may be dangerous. In that
case, it is better to have some air outside the yoke in the model.
Furthermore, because the shaft is generally excluded from the analysis domain,
the Dirichlet boundary condition is also applied between the shaft and the
rotor lamination.

Considering the boundary conditions in the context of an asynchronous machine
and both geometric and magnetic symmetry conditions, Dirichlet boundary
conditions are assigned on the external boundary of the stator and Neumann
boundary conditions $\frac{\partial A_{3}}{\partial x_{3}}=0$, in the normal
direction, with respect to the $x_{1}x_{2}$-plane in this context, are
assigned on the remaining boundary of the sector constituted of either one
pole, a pair of poles or a set of pole pairs. This condition is equivalent to
applying a given value to the derivative of $\boldsymbol{A}$ normal to the
boundary, resulting in flux lines (which are the integral of the normal
component of the magnetic field $\boldsymbol{B}$) with a given incidence angle
with the boundary.

Finally, for other types of electrical machines and also in a magnetically
asymmetric context, one can still use periodic boundary conditions assigned
along the lateral boundary of the sector. When these are taken into account,
the remaining external circumference is subjected to Dirichlet boundary
conditions.

In terms of initial conditions, the machine is typically simulated starting
from rest with an assigned current or current density assigned to specific
coils.

More details on the space and time discretization are discussed in detail in
[Bia05, Sal95] along with the solution of the resulting discretized finite
dimensional system. If, in addition, we assume that all time-varying
electromagnetic quantities in the problem have the form

$Z(t)=Z_{\mathrm{m}}\sin(\omega t+\theta)\\\ $ (9)

where $Z_{\mathrm{m}}$ is the magnitude, $\omega$ is the angular frequency and
$\theta$ is the phase shift, then equation (8) can be rewritten as

$\nabla\times(\frac{1}{\mu}\nabla\times\boldsymbol{A})=(\sigma+j\omega\epsilon)(-j\omega\boldsymbol{A}-\nabla\Phi).\\\
$ (10)

The total current flowing in the conductor $\boldsymbol{I_{\mathrm{T}}}$
equals the integral of the current density over the cross-section of the
conductor $\Omega$ resulting in the equation

$\boldsymbol{I_{\mathrm{T}}}=\int\limits_{\Omega}{(\sigma+j\omega\epsilon)(-j\omega\boldsymbol{A}-\nabla\Phi)d\Omega}.$
(11)

######  Assumption 2.1.

For (10) to hold, all quantities must have the same frequency $\omega$, but
can have different phase angles $\theta$. If a quantity is not purely
sinusoidal, it can be decomposed into its sinusoidal harmonics and solved
separately for each harmonic (in the linear case) using the superposition
principle.

The solution of equations (8) or (10) requires a discretization scheme [MMM08,
Isl10].

The computation time for a complex problem may take a few hours especially in
the 3D case. Such computing times cannot be tolerated in a digital twin in
real operation. For this, it is necessary to use model reduction to be closer
to the real time requirements. Some work in this direction has been proposed
[Cod15, SI13, MSK+21] using for example proper orthogonal decomposition (POD)
[SI13] or Gappy POD and the discrete empirical interpolation method (DEIM)
[MSK+21]. In addition, model reduction was used in unsymmetrical problems for
rotating machines [BLRS18] when small perturbations in the geometry or the
material parameters are introduced during the mass production process and have
to be taken into account in the model treated by the digital twin for the
specific electrical machine used in the factory.

Model uncertainty is another challenge that has to be taken into account in
digital twins. This is typically done by incorporating the data gathered from
the physical twin to update the parameters in the model of the digital twin.
An application in this direction can be found in [KKH+20].

In terms of simulation software, Ansys Maxwell is one of the industrial
standards in solving the equations presented in this subsection. It makes the
design and analysis of these electrical machines easy with preset models of
different motors classes and provides also a readily available visualization
interface. However, the software lacks flexibility and transparency in terms
of model reduction that is an essential part of digital twins. More recently,
Ansys presented Ansys Twin Builder to answer this question. However, an effort
can still be made to use the wide range of model reduction techniques existing
in the literature, especially for coupled systems. In addition, a model
hierarchy is also required to have the best model for a specific application,
see for instance [BDD+04, DHL+21, LFW21] for some examples of model
hierarchies in different application fields from engineering and physics. In a
digital twin where model uncertainty and fast computations have to be taken
into account, Ansys has to be more open and flexible concerning model
extraction and manipulation.

Open-source software, such as e.g. FEniCS, can also be used to solve the
equations mentioned above. The mathematical model would then be available and
could be adapted, reduced, or coupled with other systems. Although this gives
more freedom, it requires the implementation of the models for different types
of machines and their coupled models. In addition, Ansys, still has more
attractive visualization tools than most open-source software implementations
for electrical machines.

In the following, we present an energy based system theoretic port-Hamiltonian
(pH) formulation for the full set of Maxwell equations including displacement
current. We first rewrite the Maxwell equations (1c,1d) based on (3b) in the
form

$\displaystyle\frac{1}{\mu}\frac{\partial\boldsymbol{B}}{\partial t}$
$\displaystyle=-\frac{1}{\mu}\nabla\times\boldsymbol{E},$ (12a)
$\displaystyle{\epsilon}\frac{\partial\boldsymbol{E}}{\partial t}$
$\displaystyle=\frac{1}{\mu}\nabla\times\boldsymbol{B}-\sigma\boldsymbol{E},$
(12b)

which may also be written as

$\begin{bmatrix}\varepsilon I_{3}&0\\\
0&\frac{1}{\mu}I_{3}\end{bmatrix}\begin{bmatrix}\dot{\boldsymbol{E}}\\\
\dot{\boldsymbol{B}}\end{bmatrix}=\begin{pmatrix}\begin{bmatrix}0&\frac{1}{\mu}\nabla\times\\\
-\frac{1}{\mu}\nabla\times&0\end{bmatrix}-\begin{bmatrix}\sigma I_{3}&0\\\
0&0\end{bmatrix}\end{pmatrix}\begin{bmatrix}\boldsymbol{E}\\\
\boldsymbol{B}\end{bmatrix}.$ (13)

If we also consider the algebraic constraints (1a) and (1b), then the complete
system may be written as

$\begin{bmatrix}\epsilon I_{3}&0\\\ 0&\frac{1}{\mu}I_{3}\\\ 0&0\\\ 0&0\\\
\end{bmatrix}\begin{bmatrix}\dot{\boldsymbol{E}}\\\
\dot{\boldsymbol{B}}\end{bmatrix}=\begin{pmatrix}\begin{bmatrix}0&\frac{1}{\mu}\nabla\times&0&\epsilon\nabla\\\
-\frac{1}{\mu}\nabla\times&0&\nabla&0\\\ 0&\nabla\cdot&0&0\\\
\epsilon\nabla\cdot&0&0&0\\\ \end{bmatrix}-\begin{bmatrix}\sigma
I_{3}&0&0&0\\\ 0&0&0&0\\\ 0&0&0&0\\\ 0&0&0&0\\\
\end{bmatrix}\end{pmatrix}\begin{bmatrix}\boldsymbol{E}\\\ \boldsymbol{B}\\\
0\\\ 0\end{bmatrix}+\begin{bmatrix}0\\\ 0\\\ 0\\\ -1\\\
\end{bmatrix}\rho_{\mathrm{c}},$ (14)

or in a more compact form as

$\mathcal{E}_{\mathrm{em}}\dot{\boldsymbol{z}}_{\mathrm{em}}=(\mathcal{J}_{\mathrm{em}}-\mathcal{R}_{\mathrm{em}})\boldsymbol{e}_{\mathrm{em}}(\boldsymbol{z}_{\mathrm{em}})+\mathcal{B}_{\mathrm{em}}u_{\mathrm{em}}$

with

$\displaystyle\mathcal{E}_{\mathrm{em}}$
$\displaystyle\vcentcolon=\begin{bmatrix}\epsilon I_{3}&0\\\
0&\frac{1}{\mu}I_{3}\\\ 0&0\\\ 0&0\end{bmatrix}\in\mathbb{R}^{8,6},\quad$
$\displaystyle\boldsymbol{z}_{\mathrm{em}}$
$\displaystyle\vcentcolon=\begin{bmatrix}\boldsymbol{E}\\\
\boldsymbol{B}\end{bmatrix},\quad\mathcal{R}_{\mathrm{em}}\vcentcolon=\begin{bmatrix}\sigma
I_{3}&0&0&0\\\ 0&0&0&0\\\ 0&0&0&0\\\
0&0&0&0\end{bmatrix}\in\mathbb{R}^{8,8},\quad\boldsymbol{e}_{\mathrm{em}}(\boldsymbol{z}_{\mathrm{em}})\vcentcolon=\begin{bmatrix}\boldsymbol{E}\\\
\boldsymbol{B}\\\ 0\\\ 0\end{bmatrix}\in\mathbb{R}^{8},$
$\displaystyle\mathcal{B}_{\mathrm{em}}$
$\displaystyle\vcentcolon=\begin{bmatrix}0\\\ 0\\\ 0\\\ -1\\\
\end{bmatrix}\in\mathbb{R}^{8},\quad$ $\displaystyle u_{\mathrm{em}}$
$\displaystyle\vcentcolon=\rho_{\mathrm{c}}.$

Furthermore, $\mathcal{J}_{\mathrm{em}}$ is a linear operator on a space of
sufficiently smooth functions mapping from the spatial domain $\mathcal{X}$ to
$\mathcal{V}_{\mathrm{em}}\vcentcolon=\mathbb{R}^{3}\times\mathbb{R}^{3}\times\mathbb{R}\times\mathbb{R}$
and is defined via

$\displaystyle\mathcal{J}_{\mathrm{em}}$
$\displaystyle\vcentcolon=\begin{bmatrix}0&\frac{1}{\mu}\nabla\times&0&\epsilon\nabla\\\
-\frac{1}{\mu}\nabla\times&0&\nabla&0\\\ 0&\nabla\cdot&0&0\\\
\epsilon\nabla\cdot&0&0&0\end{bmatrix}.$

When also considering prescribed Dirichlet boundary values
$\boldsymbol{u}_{\mathrm{em,b}}$ for the state $\boldsymbol{z}_{\mathrm{em}}$,
we may write the governing equations together with the boundary conditions in
the form

$\begin{bmatrix}\mathcal{E}_{\mathrm{em}}\\\
0\end{bmatrix}\dot{\boldsymbol{z}}_{\mathrm{em}}=\begin{bmatrix}\mathcal{J}_{\mathrm{em}}-\mathcal{R}_{\mathrm{em}}\\\
-\mathcal{U}_{\mathrm{em}}\end{bmatrix}\boldsymbol{e}_{\mathrm{em}}(\boldsymbol{z}_{\mathrm{em}})+\begin{bmatrix}\mathcal{B}_{\mathrm{em}}&0\\\
0&I_{6}\end{bmatrix}\begin{bmatrix}u_{\mathrm{em}}\\\
\boldsymbol{u}_{\mathrm{em,b}}\end{bmatrix},$ (15)

where the linear operator $\mathcal{U}_{\mathrm{em}}$ takes sufficiently
smooth functions mapping from $\mathcal{X}$ to $\mathcal{V}_{\mathrm{em}}$ and
returns functions mapping from $\partial\mathcal{X}$ to $\mathbb{R}^{6}$ via
the rule

$\mathcal{U}_{\mathrm{em}}\left(\begin{bmatrix}\boldsymbol{\eta}_{1}\\\
\boldsymbol{\eta}_{2}\\\ \eta_{3}\\\
\eta_{4}\end{bmatrix}\right)\vcentcolon=\begin{bmatrix}\left.\boldsymbol{\eta}_{1}\right|_{\partial\mathcal{X}}\\\
\left.\boldsymbol{\eta}_{2}\right|_{\partial\mathcal{X}}\end{bmatrix}.$

One can show that the special form (15) of the governing equations indeed
corresponds to a port-Hamiltonian formulation as introduced in [DHL+21, sec.
9.2]. To this end, we consider the Hamiltonian

$\displaystyle\mathcal{H}_{\mathrm{em}}(\boldsymbol{E},\boldsymbol{B})$
$\displaystyle\vcentcolon=\int\limits_{\mathcal{X}}\frac{1}{2}\left(\epsilon\boldsymbol{E}^{\top}\boldsymbol{E}+\frac{1}{\mu}\boldsymbol{B}^{\top}\boldsymbol{B}\right)\,\mathrm{d}\mathcal{X},$
(16)

which represents the energy of the electromagnetic field, see for instance
[Jac99, eq. (6.106)]. The essential properties which characterize the port-
Hamiltonian structure of (15) are the identity

$\mathcal{E}_{\mathrm{em}}^{\top}\boldsymbol{e}_{\mathrm{em}}(\boldsymbol{z}_{\mathrm{em}})=\frac{\delta\mathcal{H}_{\mathrm{em}}}{\delta\boldsymbol{z}_{\mathrm{em}}}(\boldsymbol{z}_{\mathrm{em}})$
(17)

and the fact that

$\beta_{\mathrm{em}}\left(\begin{bmatrix}\boldsymbol{\eta}_{1}\\\
\boldsymbol{\eta}_{2}\\\ \eta_{3}\\\
\eta_{4}\end{bmatrix},\begin{bmatrix}\mathcal{J}_{\mathrm{em}}-\mathcal{R}_{\mathrm{em}}\\\
-\mathcal{U}_{\mathrm{em}}\end{bmatrix}\begin{bmatrix}\boldsymbol{\eta}_{1}\\\
\boldsymbol{\eta}_{2}\\\ \eta_{3}\\\
\eta_{4}\end{bmatrix}\right)=-\int\limits_{\mathcal{X}}\sigma\boldsymbol{\eta}_{1}^{\top}\boldsymbol{\eta}_{1}\,\mathrm{d}\mathcal{X}\leq
0$ (18)

holds for all sufficiently smooth
$(\boldsymbol{\eta}_{1},\boldsymbol{\eta}_{2},\eta_{3},\eta_{4})\colon\mathcal{X}\to\mathcal{V}_{\mathrm{em}}$,
where the bilinear form $\beta_{\mathrm{em}}$ is defined via

$\displaystyle\beta_{\mathrm{em}}\left(\begin{bmatrix}\boldsymbol{\eta}_{1}\\\
\boldsymbol{\eta}_{2}\\\ \eta_{3}\\\
\eta_{4}\end{bmatrix},\begin{bmatrix}\boldsymbol{\xi}_{1}\\\
\boldsymbol{\xi}_{2}\\\ \xi_{3}\\\ \xi_{4}\\\ \boldsymbol{\xi}_{5}\\\
\boldsymbol{\xi}_{6}\end{bmatrix}\right)$
$\displaystyle\vcentcolon=\int\limits_{\mathcal{X}}\left(\sum_{i=1}^{2}\boldsymbol{\eta}_{i}^{\top}\boldsymbol{\xi}_{i}+\sum_{j=3}^{4}\eta_{j}\xi_{j}\right)\,\mathrm{d}\mathcal{X}+\int\limits_{\partial\mathcal{X}}\mathcal{Y}_{\mathrm{em}}(\boldsymbol{\eta}_{1},\boldsymbol{\eta}_{2},\eta_{3},\eta_{4})^{\top}\begin{bmatrix}\boldsymbol{\xi}_{5}\\\
\boldsymbol{\xi}_{6}\end{bmatrix}\,\mathrm{d}S$
$\displaystyle=\int\limits_{\mathcal{X}}\left(\sum_{i=1}^{2}\boldsymbol{\eta}_{i}^{\top}\boldsymbol{\xi}_{i}+\sum_{j=3}^{4}\eta_{j}\xi_{j}\right)\,\mathrm{d}\mathcal{X}+\int\limits_{\partial\mathcal{X}}\left(\frac{1}{\mu}\boldsymbol{\xi}_{6}\times\boldsymbol{\eta}_{1}+\epsilon\eta_{4}\boldsymbol{\xi}_{5}+\eta_{3}\boldsymbol{\xi}_{6}\right)^{\top}\boldsymbol{\nu}\,\mathrm{d}S.$

Here, $\boldsymbol{\nu}$ denotes the outer unit-normal vector field on the
boundary $\partial\mathcal{X}$ and the linear operator
$\mathcal{Y}_{\mathrm{em}}$ takes sufficiently smooth functions mapping from
$\mathcal{X}$ to $\mathcal{V}_{\mathrm{em}}$ and returns functions mapping
from $\partial\mathcal{X}$ to $\mathbb{R}^{6}$ via the rule

$\mathcal{Y}_{\mathrm{em}}\left(\begin{bmatrix}\boldsymbol{\eta}_{1}\\\
\boldsymbol{\eta}_{2}\\\ \eta_{3}\\\
\eta_{4}\end{bmatrix}\right)\vcentcolon=\begin{bmatrix}\left.\epsilon\eta_{4}\right|_{\partial\mathcal{X}}\boldsymbol{\nu}\\\
\left.\left(\frac{1}{\mu}\varepsilon\boldsymbol{\eta}_{1}+\eta_{3}I_{3}\right)^{\top}\right|_{\partial\mathcal{X}}\boldsymbol{\nu}\end{bmatrix},$

where $\varepsilon\in\mathbb{R}^{3,3,3}$ denotes the third-order tensor
corresponding to the Levi–Civita symbol, see for instance [RHB06, sec. 26.8].
An important property of a port-Hamiltonian system is that its structure
encodes an energy balance. More precisely, if (17) and (18) are satisfied, we
obtain the dissipation inequality

$\displaystyle\frac{\mathrm{d}\mathcal{H}_{\mathrm{em}}}{\mathrm{d}t}$
$\displaystyle=\int\limits_{\mathcal{X}}\left(\dot{\boldsymbol{E}}^{\top}\frac{\delta\mathcal{H}_{\mathrm{em}}}{\delta\boldsymbol{E}}(\boldsymbol{z}_{\mathrm{em}})+\dot{\boldsymbol{B}}^{\top}\frac{\delta\mathcal{H}_{\mathrm{em}}}{\delta\boldsymbol{B}}(\boldsymbol{z}_{\mathrm{em}})\right)\,\mathrm{d}\mathcal{X}=\beta_{\mathrm{em}}\left(\boldsymbol{e}_{\mathrm{em}}(\boldsymbol{z}_{\mathrm{em}}),\begin{bmatrix}\mathcal{E}_{\mathrm{em}}\dot{\boldsymbol{z}}_{\mathrm{em}}\\\
0\end{bmatrix}\right)$ (19)
$\displaystyle=\beta_{\mathrm{em}}\left(\boldsymbol{e}_{\mathrm{em}}(\boldsymbol{z}_{\mathrm{em}}),\begin{bmatrix}\mathcal{J}_{\mathrm{em}}-\mathcal{R}_{\mathrm{em}}\\\
-\mathcal{U}_{\mathrm{em}}\end{bmatrix}\boldsymbol{e}_{\mathrm{em}}(\boldsymbol{z}_{\mathrm{em}})+\begin{bmatrix}\mathcal{B}_{\mathrm{em}}&0\\\
0&I_{6}\end{bmatrix}\begin{bmatrix}u_{\mathrm{em}}\\\
\boldsymbol{u}_{\mathrm{em,b}}\end{bmatrix}\right)$
$\displaystyle=-\int\limits_{\mathcal{X}}\sigma\boldsymbol{E}^{\top}\boldsymbol{E}\,\mathrm{d}\mathcal{X}+\beta_{\mathrm{em}}\left(\boldsymbol{e}_{\mathrm{em}}(\boldsymbol{z}_{\mathrm{em}}),\begin{bmatrix}\mathcal{B}_{\mathrm{em}}u_{\mathrm{em}}\\\
\boldsymbol{u}_{\mathrm{em,b}}\end{bmatrix}\right)$
$\displaystyle=-\int\limits_{\mathcal{X}}\sigma\boldsymbol{E}^{\top}\boldsymbol{E}\,\mathrm{d}\mathcal{X}+\int\limits_{\mathcal{X}}\underbrace{0}_{=\vcentcolon
y_{\mathrm{em}}}\cdot
u_{\mathrm{em}}\,\mathrm{d}\mathcal{X}+\int\limits_{\partial\mathcal{X}}\underbrace{\mathcal{Y}_{\mathrm{em}}(\boldsymbol{e}_{\mathrm{em}}(\boldsymbol{z}_{\mathrm{em}}))^{\top}}_{=\vcentcolon\boldsymbol{y}_{\mathrm{em,b}}^{\top}}\boldsymbol{u}_{\mathrm{em,b}}\,\mathrm{d}S$
$\displaystyle=-\int\limits_{\mathcal{X}}\sigma\boldsymbol{E}^{\top}\boldsymbol{E}\,\mathrm{d}\mathcal{X}+\int\limits_{\partial\mathcal{X}}\frac{1}{\mu}\left(\boldsymbol{B}\times\boldsymbol{E}\right)^{\top}\boldsymbol{\nu}\,\mathrm{d}S\leq\int\limits_{\partial\mathcal{X}}\frac{1}{\mu}\left(\boldsymbol{B}\times\boldsymbol{E}\right)^{\top}\boldsymbol{\nu}\,\mathrm{d}S.$

This shows that the energy of the electromagnetic field within the spatial
domain $\mathcal{X}$ may only change due to dissipation and due to the energy
flow through the boundary. The variables $y_{\mathrm{em}}$ and
$\boldsymbol{y}_{\mathrm{em,b}}$ introduced in (19) represent the dual port
variables to the inputs $u_{\mathrm{em}}$ and
$\boldsymbol{u}_{\mathrm{em,b}}$, respectively. These are defined such that

$\int\limits_{\mathcal{X}}y_{\mathrm{em}}\cdot
u_{\mathrm{em}}\,\mathrm{d}\mathcal{X}\quad\text{and}\quad\int\limits_{\partial\mathcal{X}}\boldsymbol{y}_{\mathrm{em,b}}^{\top}\boldsymbol{u}_{\mathrm{em,b}}\,\mathrm{d}S$

correspond to the energy change caused by the distributed input
$u_{\mathrm{em}}$ and the boundary input $\boldsymbol{u}_{\mathrm{em,b}}$,
respectively. In particular, we have $y_{\mathrm{em}}=0$, which reflects that
the input $u_{\mathrm{em}}=\rho_{\mathrm{c}}$ has no effect on the energy of
the electromagnetic field. In the pH literature it is common to explicitly add
the definitions of these dual port variables to the pH formulation, see for
instance [DHL+21].

It remains to show that (17) and (18) indeed hold. To this end, we first
compute the variational derivatives of the Hamiltonian which are given by

$\frac{\delta\mathcal{H}_{\mathrm{em}}}{\delta\boldsymbol{E}}(\boldsymbol{z}_{\mathrm{em}})=\epsilon\boldsymbol{E},\quad\frac{\delta\mathcal{H}_{\mathrm{em}}}{\delta\boldsymbol{B}}(\boldsymbol{z}_{\mathrm{em}})=\frac{1}{\mu}\boldsymbol{B}.$

From this we immediately obtain that (17) is satisfied. To verify (18), we
compute

$\displaystyle\beta_{\mathrm{em}}\left(\begin{bmatrix}\boldsymbol{\eta}_{1}\\\
\boldsymbol{\eta}_{2}\\\ \eta_{3}\\\
\eta_{4}\end{bmatrix},\begin{bmatrix}\mathcal{J}_{\mathrm{em}}-\mathcal{R}_{\mathrm{em}}\\\
-\mathcal{U}_{\mathrm{em}}\end{bmatrix}\begin{bmatrix}\boldsymbol{\eta}_{1}\\\
\boldsymbol{\eta}_{2}\\\ \eta_{3}\\\ \eta_{4}\end{bmatrix}\right)$
$\displaystyle=\int\limits_{\mathcal{X}}\left(\frac{1}{\mu}\left(\boldsymbol{\eta}_{1}^{\top}\nabla\times\boldsymbol{\eta}_{2}-\boldsymbol{\eta}_{2}^{\top}\nabla\times\boldsymbol{\eta}_{1}\right)+\epsilon(\boldsymbol{\eta}_{1}^{\top}\nabla\eta_{4}+\eta_{4}\nabla\cdot\boldsymbol{\eta}_{1})+\boldsymbol{\eta}_{2}^{\top}\nabla\eta_{3}+\eta_{3}\nabla\cdot\boldsymbol{\eta}_{2}-\sigma\boldsymbol{\eta}_{1}^{\top}\boldsymbol{\eta}_{1}\right)\,\mathrm{d}\mathcal{X}$
$\displaystyle\quad-\int\limits_{\partial\mathcal{X}}\left(\frac{1}{\mu}\boldsymbol{\eta}_{2}\times\boldsymbol{\eta}_{1}+\epsilon\eta_{4}\boldsymbol{\eta}_{1}+\eta_{3}\boldsymbol{\eta}_{2}\right)^{\top}\boldsymbol{\nu}\,\mathrm{d}S$
$\displaystyle=\int\limits_{\mathcal{X}}\nabla\cdot\left(\frac{1}{\mu}\boldsymbol{\eta}_{2}\times\boldsymbol{\eta}_{1}+\epsilon\eta_{4}\boldsymbol{\eta}_{1}+\eta_{3}\boldsymbol{\eta}_{2}\right)\,\mathrm{d}\mathcal{X}-\int\limits_{\mathcal{X}}\sigma\boldsymbol{\eta}_{1}^{\top}\boldsymbol{\eta}_{1}\,\mathrm{d}\mathcal{X}-\int\limits_{\partial\mathcal{X}}\left(\frac{1}{\mu}\boldsymbol{\eta}_{2}\times\boldsymbol{\eta}_{1}+\epsilon\eta_{4}\boldsymbol{\eta}_{1}+\eta_{3}\boldsymbol{\eta}_{2}\right)^{\top}\boldsymbol{\nu}\,\mathrm{d}S$
$\displaystyle=-\int\limits_{\mathcal{X}}\sigma\boldsymbol{\eta}_{1}^{\top}\boldsymbol{\eta}_{1}\,\mathrm{d}\mathcal{X},$

which holds for all sufficiently smooth
$(\boldsymbol{\eta}_{1},\boldsymbol{\eta}_{2},\eta_{3},\eta_{4})\colon\mathcal{X}\to\mathcal{V}_{\mathrm{em}}$.

######  Remark 2.2.

We emphasize that port-Hamiltonian formulations are in general not unique and
this is in particular true for (15). This non-uniqueness also includes the
choice of state variables. An alternative port-Hamiltonian formulation for the
Maxwell’s equations could for instance be derived using the state variables
$\boldsymbol{E}$ and $\boldsymbol{H}$.

######  Remark 2.3.

In general, we may distinguish two different kinds of boundary conditions in
the framework of port-Hamiltonian partial differential-algebraic equations
(PDAEs): those which encode the change of the Hamiltonian due to boundary
flows and those which are not directly linked with the Hamiltonian, but which
may be important for ensuring that the PDAE has a unique solution. The first
kind of boundary conditions is especially important for a system-theoretic
description, since the corresponding boundary port variables may be used for a
power-preserving coupling between several subsystems which exchange energy
through the subsystem boundaries. For instance, in (15) the boundary port
variable
$\boldsymbol{u}_{\mathrm{em,b}}=\left.\boldsymbol{z}_{\mathrm{em}}\right|_{\partial\mathcal{X}}$
corresponds to prescribed Dirichlet values of the electric and the magnetic
field at the boundary. It may be used together with the dual port variable
$\boldsymbol{y}_{\mathrm{em,b}}$ to couple the system in a power-conserving
way to other subsystems, for instance in the case, where $\partial\mathcal{X}$
corresponds to an interface between different electromagnetic fields. In
general, we may also add other boundary conditions like Neumann boundary
conditions to the pH formulation (15), without destroying the pH structure.

To summarize, in the simulation of electrical rotating machines, the following
assumptions are generally used:

  * •

Linear materials are assumed to obtain the linear relations in equations (3)
and (4) that include polarization and magnetization microscopic effects.

  * •

When considering only electromagnetic simulations, the mechanical effects such
as velocity-dependent effects and deformation are neglected.

  * •

If all considered inputs are sinusoidal with the same frequency then one can
simplify the computation of the current density and total current in (10) and
(11).

###  2.2 Steady State Equivalent Circuit Model

It is popular in electrical engineering to use lumped-parameter and surrogate
models. These parameters are represented by impedance values, voltage source
and current sources that constitute an equivalent circuit of the machine.
These calculation models are extremely fast but do not have any guarantees on
error bounds. This is the reason why these models are placed at the bottom of
our hierarchy of models. The surrogate models mainly used in steady state
electromagnetic simulation are based on equivalent circuit models of the
electrical machine. These models are specific to the type of machines being
studied. For asynchronous AC machines, the equivalent circuit can be derived
through a series of simplifications and approximations [Mel18, Ch. 4].

Assuming that the stator current results in a magnetomotive force (MMF)
$F_{1}$ and the rotor current yields an MMF $F_{2}$, both MMFs rotate
synchronously. The resulting magnetic field is just the vectorial sum of the
two MMFs. If the magnetic circuit is assumed to be linear, then the magnetic
fields, the fluxes and the electromotive forces (EMFs) from the stator and
rotor currents can be summed up. If in addition, we add the resistive and
inductive leakage voltage drops in the windings, the resulting equivalent
circuit for the steady state operating behavior of an asynchronous machine on
a symmetrical multiphase network is represented in figure 4.

$U_{1}$$I_{1}$$R_{1}$$X_{1}$$X_{3}$$I^{\prime}_{2}$$X^{\prime}_{2}$$\frac{R^{\prime}_{2}}{s}$$U_{2}$
Figure 4: Equivalent circuit model of an induction machine

The corresponding equations can be derived as

$\displaystyle{{U}_{1}}={{R}_{1}}\cdot{{I}_{1}}+j{{X}_{1}}\cdot{{I}_{1}}+j{{X}_{3}}\cdot\left({{I}_{1}}-I_{2}^{{}^{\prime}}\right),$
(20) $\displaystyle U_{2}^{{}^{\prime}}=-\frac{R_{2}^{{}^{\prime}}}{s}\cdot
I_{2}^{{}^{\prime}}-jX_{2}^{{}^{\prime}}\cdot
I_{2}^{{}^{\prime}}+j{{X}_{3}}\cdot\left({{I}_{1}}-I_{2}^{{}^{\prime}}\right),$
(21)

where ${U}_{1}$ is the input voltage, ${I}_{1}$ and $I_{2}^{{}^{\prime}}$ are
the stator and the rotor currents respectively with the prime indicating
quantities referred to the stator winding. It is important to note that
sometimes, depending on the machine type, $U_{2}^{{}^{\prime}}$ is set to zero
(short circuited) which simplifies the computations. $s$ is the slip that
describes the relation between the rotor speed $n$ and the synchronous speed
$n_{\mathrm{s}}$ of the stator flux as

$s=\frac{n_{\mathrm{s}}-n}{n_{\mathrm{s}}}.$ (22)

All the impedance values are chosen to mimic the interactions between the
stator and the rotor. These values are generally estimated using frequency
measurements from the machine [Sal95]. In view of this, the equivalent circuit
model can be considered to be a surrogate model.

Alternatively, the impedance values can be generated based on some simple FEM
calculations that fit into the time requirements of the digital twin [MMS+20].

The equations (20) are generally solved in phasor form resulting in a set of
linear equations [Nür52, Mel18]. The solution gives approximate values for the
currents in the stator and the rotor which are then used to compute power and
torque. This provides a simple and fast method to compute the currents and can
be directly implemented in the digital twin. This also allows for a fast
computation of the losses as discussed in the following section.

###  2.3 Computation of Losses

One important factor in the study of electrical machines and their efficiency
is the total power loss of the machine. This represents all the wasted power
during the conversion from electrical energy to mechanical energy in the case
of motors and the conversion from mechanical to electrical energy in the case
of generators. Reducing this loss results in an increased efficiency of the
machine. For this reason, monitoring the losses during the lifetime of the
machine within the digital twin can help to optimize the design of new
machines for a reduced loss and an increased efficiency. In some cases, it is
even possible to increase the efficiency online by changing the flux and
current to reduce the total loss $P_{\mathrm{V}}$ which is computed as the sum

$P_{\mathrm{V}}=P_{\mathrm{Fe}}+P_{\mathrm{R}}+P_{\mathrm{w}1}+P_{\mathrm{w}2}+P_{\mathrm{Z}}$
(23)

where $P_{\mathrm{Fe}}$ is the iron core loss (sometimes also called magnetic
loss), i.e., the loss that occurs in the stator and rotor cores due to
hysteresis and eddy currents. The friction losses resulting from the
mechanical rotation are summarized in $P_{\mathrm{R}}$. More details on these
losses can be found in [YBK+17]. $P_{\mathrm{w}1}$ is the stator winding loss
and $P_{\mathrm{w}2}$ is the rotor winding loss. These two losses are the
result of the inherent resistive properties of windings and conductors.
$P_{\mathrm{Z}}$ are all the additional losses that are difficult to quantify
and model. Typically these are taken as a global percentage of the total
power. More details can be found in [IMR21].

For the computation of the losses in the stator and the rotor, the
corresponding currents are required as shown in the equations

$P_{\mathrm{w}1}=N_{\mathrm{p}}R_{1}I_{1}^{2},\quad
P_{\mathrm{w}2}=N_{\mathrm{p}}R_{2}^{{}^{\prime}}I_{2}^{{}^{\prime}2},$ (24)

where $N_{\mathrm{p}}$ indicates the number of phases of the system. The
currents are typically computed directly based on the equivalent circuit
described in the previous subsection. It is important to note that the
material properties change with increased temperatures, hence the resistance
changes with the temperature of the conductor. Typically, these changes are
neglected when using this type of simulation models. Instead, a typical
operating temperature is used to compute for example $R_{1}$ and
$R_{2}^{{}^{\prime}}$.

On the other hand, using computational models based on discretizing the PDEs,
as discussed in section 2.1, may be computationally prohibitively expensive.
If in addition one wants to do parameter optimization, then repeated
simulations for different parameter configurations are required, rendering the
determination of a solution in real-time within the digital twin impossible.
One approach for avoiding such difficulties is to construct a parametric model
that could be reduced using model order reduction and used in the optimization
step. In addition, depending on the model reduction technique used, one may
also obtain error bounds. Alternatively, one could use the equivalent circuit
model, however, the results may not be as accurate and are typically without
error estimates.

###  2.4 Condition Monitoring and Predictive Maintenance

While the electrical machine is operating, the simulation models in the
digital twin are used to validate the measurements from the machine and detect
any deviations from the normal operation of the machine. This is why the
simulation models have to be accurate and fast in order to adapt to any
changes in the parameters. In addition to model validation and fault
detection, the digital twin has to be able to predict faults and when a
maintenance should be planned.

As a use case, let us consider one of the most common faults that lead to a
break down of rotating machines, i.e., faults in the insulation system [Sch90,
GALH08, Tle08]. In order to have a good prediction of faults in the insulation
system, we need on the one hand a good lifetime model of the insulation
material and on the other hand we should be able to detect the causes of
failure in the insulation system.

The study of the insulation system and its lifetime in an electrical machine
is a challenging task. The only reasonable models derived from physics are on
a very small scale level, resulting in a very complex full model if that is to
be used in all parts of the machine.

Thus this constitutes a very good use case that showcases how a digital twin
can improve predictive maintenance by using the data collected on the machine
to construct a surrogate model.

The main challenge in such an approach is to obtain enough data that include
information about the failure of the insulation system. Since the acquisition
of this data is challenging as it may require provoking failures in the
insulation system which may be very expensive or risky, instead one could use
material specimen to generate the necessary data, usually a stator bar or
coil. Some results in this direction are available in [Sch18]. However, these
theoretical results are subject to deviations in practice. In addition, it is
hard to validate these kind of models.

The data gathered by a digital twin may allow for more accurate surrogate
models that furthermore can be improved over time as the digital twin
assimilates the new data coming into its data flow system.

The second aspect to consider is to detect the causes of complete failure in
parts of the insulation system. The main cause for this is partial discharge,
[Mey62, Küc09, Sch16], which are local electric discharges in an insulation
system when the field strength exceeds what the insulating material can
withstand. A repetition of such discharges damages the insulating material
permanently. The digital twin can be used to detect and determine
approximately where partial discharges are happening. A warning is then sent
to the operator in order to plan in advance a maintenance before a fault
happens that leads to a stop in the operation of the machine.

It is also important to note that keeping track of the temporal trend of
partial discharges is crucial for interpreting the data. This is done by
monitoring online the voltage at each phase and measuring the partial
discharge at each phase.

This approach can be used to determine in which phase the partial discharge
did happen and one may even determine approximately where the partial
discharge happened by triangulating the distance to the failure point based on
the distances to the other measurement points.

The detection of the evolution of partial discharges can also be done by
routinely connecting coupling capacitors to each phase of the input voltage. A
data acquisition unit is then connected to the coupling capacitors which
filters the signal to keep only the partial discharge. The measured partial
discharge data is then interpreted using surrogate modeling and machine
learning to classify the partial discharges and to know when and where a
partial discharge occurs [AVR+11]. In its simplest form, one would just
compare the measured data with user-defined thresholds. A relevant
international standard, that specifies how to monitor partial discharges on
rotating machines, is presented in IEC 60034-27-2.

##  3 Mechanical Simulation

An important step of the design process of electrical machines is the
determination of maximal mechanical stresses, since they limit (among other
effects) the power of an electrical machine, cf. [GZP21]. Furthermore,
vibration and noise analysis are typical mechanical calculations performed for
an electrical machine in order to reduce acoustic noise emissions, cf.
[BSKDD10, GWL06, JLCLM95, PLF12, RHKH96, XHC18]. Both aspects are discussed in
section 3.1, whereas we address the relationship between the electromagnetic
torque and the resulting rotational speed in section 3.2 via Newton’s second
law for rotation.

###  3.1 Analysis of Displacements and Vibrations

In general, in kinematics one is interested in deformations and motions of a
physical body with respect to a reference configuration
$\mathcal{X}\subset\mathbb{R}^{3}$. A configuration or deformation is
described by a mapping from the reference configuration $\mathcal{X}$ to
$\mathbb{R}^{3}$ and a motion is a time-dependent family of deformations, see
for instance [MH83] for more details. The simplest examples for motions are
rigid body motions such as translations or rotations. In general, the time
evolution of the motion is governed by the balance of momentum which, for a
solid continuous body, can be written in Lagrangian coordinates as

$\rho(\boldsymbol{x})\frac{\partial^{2}\boldsymbol{\phi}}{\partial
t^{2}}(t,\boldsymbol{x})=\nabla\cdot
P(t,\boldsymbol{x})+\rho(\boldsymbol{x})\boldsymbol{F}(t,\boldsymbol{x}),$
(25)

where $\boldsymbol{x}\in\mathcal{X}$ is an arbitrary point in the reference
configuration, $\boldsymbol{\phi}\colon\mathbb{R}_{\geq
0}\times\mathcal{X}\to\mathbb{R}^{3}$ is the motion of the considered body,
$P\colon\mathbb{R}_{\geq 0}\times\mathcal{X}\to\mathbb{R}^{3,3}$ denotes the
first Piola–Kirchhoff stress tensor, $\rho\colon\mathcal{X}\to\mathbb{R}_{>0}$
is the density in the reference configuration, and
$\boldsymbol{F}\colon\mathbb{R}_{\geq 0}\times\mathcal{X}\to\mathbb{R}^{3}$ is
an external body force, cf. [MH83, ch. 1–2]. Here, the divergence of $P$ is
defined via

$[\nabla\cdot P]_{i}\vcentcolon=\sum_{j=1}^{3}\frac{\partial P_{ij}}{\partial
x_{j}}\quad\text{for }i=1,\ldots,3.$

The unknown of (25) is the motion $\boldsymbol{\phi}$, whereas the density
$\rho$ and the external force $\boldsymbol{F}$ are assumed to be given. For
instance, in the context of electrical machines, the latter term usually
includes electromagnetic forces, see for instance [XHC18]. The stress tensor
$P$ is specified by means of constitutive laws. For instance, a homogeneous
and elastic material can be described by a constitutive law of the form

$P(t,\boldsymbol{x})=P_{\mathrm{c}}\left(\frac{\partial\boldsymbol{\phi}}{\partial\boldsymbol{x}}(t,\boldsymbol{x})\right),$

which is characterized by the constitutive function
$P_{\mathrm{c}}\colon\mathbb{R}^{3,3}\to\mathbb{R}^{3,3}$, cf. [MH83, p. 8].
The tensor $\frac{\partial\boldsymbol{\phi}}{\partial\boldsymbol{x}}$ is
called the deformation tensor and denoted with $F$ in the following. If one
considers viscous dissipation, then the constitutive law of a homogeneous
material takes the form

$P(t,\boldsymbol{x})=P_{\mathrm{c}}\left(F(t,\boldsymbol{x}),\frac{\partial
F}{\partial t}(t,\boldsymbol{x})\right),$ (26)

with constitutive function
$P_{\mathrm{c}}\colon\mathbb{R}^{3,3}\times\mathbb{R}^{3,3}\to\mathbb{R}^{3,3}$,
cf. [Ant05, sec. 13.9]. Finally, (25) is closed by initial conditions

$\boldsymbol{\phi}(0,\boldsymbol{x})=\boldsymbol{\phi}_{0}(\boldsymbol{x}),\quad\frac{\partial\boldsymbol{\phi}}{\partial
t}(0,\boldsymbol{x})=\boldsymbol{v}_{0}(\boldsymbol{x})$

with initial deformation $\boldsymbol{\phi}_{0}$ and velocity
$\boldsymbol{v}_{0}$, as well as by appropriate boundary conditions. Typical
boundary conditions are displacement boundary conditions of the form

$\boldsymbol{\phi}(t,\boldsymbol{x})=\boldsymbol{\phi}_{\mathrm{b}}(t,\boldsymbol{x})$

for $\boldsymbol{x}\in\Gamma_{1}$ with prescribed boundary deformation
$\boldsymbol{\phi}_{\mathrm{b}}$ and traction boundary conditions of the form

$P(t,\boldsymbol{x})\boldsymbol{\nu}(\boldsymbol{x})=\boldsymbol{\tau}_{\mathrm{b}}(t,\boldsymbol{x})$
(27)

for $\boldsymbol{x}\in\Gamma_{2}$ with
$\Gamma_{1}\cup\Gamma_{2}=\partial\mathcal{X}$,
$\Gamma_{1}\cap\Gamma_{2}=\emptyset$, outer unit-normal vector field
$\boldsymbol{\nu}$ on the boundary $\partial\mathcal{X}$, and a prescribed
boundary traction $\boldsymbol{\tau}_{\mathrm{b}}$. For some examples of such
boundary conditions, we refer to [MH83, p. 11 ff.].

######  Remark 3.1.

In this work, we only consider constitutive laws of the form (26). However, we
emphasize that in general $P_{\mathrm{c}}$ may not only depend on the values
of $F=\frac{\partial\boldsymbol{\phi}}{\partial\boldsymbol{x}}$ and
$\frac{\partial F}{\partial t}=\frac{\partial^{2}\boldsymbol{\phi}}{\partial
t\partial\boldsymbol{x}}$, but on $\boldsymbol{\phi}$ as a function. Thus,
also integrals or higher-order derivatives of $\boldsymbol{\phi}$ may occur in
the constitutive laws. Furthermore, the stress may not only depend on the
motion $\boldsymbol{\phi}$, but also on the temperature or the time, and for
non-homogeneous materials even on the position $\boldsymbol{x}$, cf. [MH83,
ch. 3]. In the context of electrical machines the temperature dependence of
the stress has for example been considered in [SBK+16, SS18]

In the following, we consider the special case

$P_{\mathrm{c}}(F_{1},F_{2})\vcentcolon=\frac{1}{2}F_{1}\left(\mathbb{H}(F_{1}^{\top}F_{1}-I_{3})+\mathbb{V}(F_{2}^{\top}F_{1}+F_{1}^{\top}F_{2})\right)$
(28)

of the constitutive law (26), which corresponds to a Kelvin–Voigt model for
the material, see for instance [GM21, LLT13]. Here,
$\mathbb{H},\mathbb{V}\in\mathbb{R}^{3,3,3,3}$ are the so-called instantaneous
or Hookean elasticity tensor and the viscosity tensor, respectively.
Furthermore, we restrict ourselves to the special case of a homogeneous and
isotropic material such that $\mathbb{H}$ and $\mathbb{V}$ can be specified as

$\mathbb{H}=\lambda_{1}I_{3}\otimes
I_{3}+2\lambda_{2}\mathbb{I}\quad\text{and}\quad\mathbb{V}=\left(\zeta_{1}-\frac{2}{3}\zeta_{2}\right)I_{3}\otimes
I_{3}+2\zeta_{2}\mathbb{I},$ (29)

cf. [LLT13]. Here $\lambda_{1},\lambda_{2}\in\mathbb{R}_{\geq 0}$ are the so-
called Lamé parameters, $\zeta_{1},\zeta_{2}\in\mathbb{R}_{\geq 0}$ are the
bulk and the shear viscosity, respectively, $I_{3}\in\mathbb{R}^{3,3}$ is the
$3\times 3$ identity matrix, $\otimes$ denotes the Kronecker product, and
$\mathbb{I}\in\mathbb{R}^{3,3,3,3}$ is the fourth-order identity tensor, i.e.,
$\mathbb{I}$ is the unique fourth-order tensor satisfying $\mathbb{I}A=A$ for
all $A\in\mathbb{R}^{3,3}$.

Using the special constitutive law (28) with $\mathbb{H}$ and $\mathbb{V}$ as
in (29), we may derive a port-Hamiltonian formulation of the governing
equations. To this end, we consider the velocity
$\boldsymbol{v}=\frac{\partial\boldsymbol{\phi}}{\partial t}$ and the
deformation tensor $F$ as variables which satisfy

$\displaystyle\rho\dot{\boldsymbol{v}}$
$\displaystyle=\frac{1}{2}\nabla\cdot\left(F\left(\mathbb{H}(F^{\top}F-I_{3})+\mathbb{V}\left(\left(\frac{\partial\boldsymbol{v}}{\partial\boldsymbol{x}}\right)^{\top}F+F^{\top}\frac{\partial\boldsymbol{v}}{\partial\boldsymbol{x}}\right)\right)\right)+\rho\boldsymbol{F},$
(30) $\displaystyle\dot{F}$
$\displaystyle=\frac{\partial\boldsymbol{v}}{\partial\boldsymbol{x}}.$

Furthermore, we introduce the associated Hamiltonian as

$\mathcal{H}_{\mathrm{m}}(F,\boldsymbol{v})=\int\limits_{\mathcal{X}}\left(\frac{\rho}{2}\boldsymbol{v}^{\top}\boldsymbol{v}+\frac{1}{8}\left\langle\mathbb{H}(F^{\top}F-I_{3}),F^{\top}F-I_{3}\right\rangle_{\mathrm{F}}\right)\,\mathrm{d}\mathcal{X},$

i.e., it is given by the sum of the kinetic and the elastic energy, cf.
[LLT13, sec. 3.3]. Based on this Hamiltonian, we may write the governing
equations (30) together with prescribed values for the boundary traction
$\boldsymbol{u}_{\mathrm{m,b}}=\boldsymbol{\tau}_{\mathrm{b}}$, cf. (27), as

$\begin{bmatrix}\mathcal{E}_{\mathrm{m}}\\\
0\end{bmatrix}\dot{\boldsymbol{z}}_{\mathrm{m}}=\begin{bmatrix}\mathcal{J}_{\mathrm{m}}-\mathcal{R}_{\mathrm{m}}(\boldsymbol{z}_{\mathrm{m}})\\\
-\mathcal{U}_{\mathrm{m}}(\boldsymbol{z}_{\mathrm{m}})\end{bmatrix}\boldsymbol{e}(\boldsymbol{z}_{\mathrm{m}})+\begin{bmatrix}\mathcal{B}_{\mathrm{m}}&0\\\
0&I_{3}\end{bmatrix}\begin{bmatrix}\boldsymbol{u}_{\mathrm{m}}\\\
\boldsymbol{u}_{\mathrm{m,b}}\end{bmatrix}.$ (31)

Here,
$\boldsymbol{z}_{\mathrm{m}}\vcentcolon=\left[\begin{smallmatrix}\boldsymbol{v}\\\
F\end{smallmatrix}\right]$ denotes the state and $\mathcal{E}_{\mathrm{m}}$ is
a linear mapping from
$\mathcal{V}_{\mathrm{m}}\vcentcolon=\mathbb{R}^{3}\times\mathbb{R}^{3,3}$ to
itself defined via

$\mathcal{E}_{\mathrm{m}}\vcentcolon=\begin{bmatrix}\rho I_{3}&0\\\
0&\mathbb{I}\end{bmatrix}.$ (32)

Furthermore, $\mathcal{J}_{\mathrm{m}}$ and
$\mathcal{R}_{\mathrm{m}}(\boldsymbol{z}_{\mathrm{m}})$ are linear operators
on a space of sufficiently smooth functions mapping from $\mathcal{X}$ to
$\mathcal{V}_{\mathrm{m}}$ and are defined via

$\displaystyle\mathcal{J}_{\mathrm{m}}$
$\displaystyle\vcentcolon=\begin{bmatrix}0&\nabla\cdot\\\
\frac{\partial}{\partial\boldsymbol{x}}&0\end{bmatrix},\quad\mathcal{R}_{\mathrm{m}}(\boldsymbol{z}_{\mathrm{m}})\vcentcolon=-\frac{1}{2}\begin{bmatrix}\nabla\cdot\left(F\left(\mathbb{V}\left(\left(\frac{\partial}{\partial\boldsymbol{x}}\cdot\right)^{\top}F+F^{\top}\frac{\partial}{\partial\boldsymbol{x}}\cdot\right)\right)\right)&0\\\
0&0\end{bmatrix}.$

The linear operator $\mathcal{U}_{\mathrm{m}}(\boldsymbol{z}_{\mathrm{m}})$
takes sufficiently smooth functions mapping from $\mathcal{X}$ to
$\mathcal{V}_{\mathrm{m}}$ and returns functions mapping from
$\partial\mathcal{X}$ to $\mathbb{R}^{3}$ via the rule

$\mathcal{U}_{\mathrm{m}}(\boldsymbol{z}_{\mathrm{m}})\left(\begin{bmatrix}\boldsymbol{\eta}_{1}\\\
\eta_{2}\end{bmatrix}\right)\vcentcolon=\left.\left(\eta_{2}+\frac{1}{2}F\left(\mathbb{V}\left(\left(\frac{\partial\boldsymbol{\eta}_{1}}{\partial\boldsymbol{x}}\right)^{\top}F+F^{\top}\frac{\partial\boldsymbol{\eta}_{1}}{\partial\boldsymbol{x}}\right)\right)\right)\boldsymbol{\nu}\right|_{\partial\mathcal{X}}.$

The linear operator $\mathcal{B}_{\mathrm{m}}$ maps from $\mathbb{R}^{3}$ to
$\mathcal{V}$ and is defined as

$\mathcal{B}_{\mathrm{m}}\vcentcolon=\rho\begin{bmatrix}I_{3}\\\
0\end{bmatrix}.$

Finally, $\boldsymbol{e}_{\mathrm{m}}(\boldsymbol{z}_{\mathrm{m}})$ and
$\boldsymbol{u}_{\mathrm{m}}$ are given by

$\boldsymbol{e}_{\mathrm{m}}(\boldsymbol{z}_{\mathrm{m}})\vcentcolon=\begin{bmatrix}\boldsymbol{v}\\\\[4.0pt]
\frac{\delta\mathcal{H}_{\mathrm{m}}}{\delta
F}(\boldsymbol{z}_{\mathrm{m}})\end{bmatrix}\quad\text{and}\quad\boldsymbol{u}_{\mathrm{m}}\vcentcolon=\boldsymbol{F}.$
(33)

To check that (30) coincides with (31), we first compute the variational
derivatives of the Hamiltonian which are given by

$\displaystyle\frac{\delta\mathcal{H}_{\mathrm{m}}}{\delta\boldsymbol{v}}(\boldsymbol{z}_{\mathrm{m}})=\rho\boldsymbol{v},\quad\frac{\delta\mathcal{H}_{\mathrm{m}}}{\delta
F}(\boldsymbol{z}_{\mathrm{m}})$
$\displaystyle=\frac{1}{2}F\left(\mathbb{H}\left(F^{\top}F-I_{3}\right)\right),$
(34)

where the latter expression follows from the calculation (86) in appendix A.
Consequently, the first block row in (31) reads

$\displaystyle\begin{bmatrix}\rho\dot{\boldsymbol{v}}\\\
\dot{F}\end{bmatrix}=\begin{bmatrix}\frac{1}{2}\nabla\cdot\left(F\left(\mathbb{H}\left(F^{\top}F-I_{3}\right)+\mathbb{V}\left(\left(\frac{\partial\boldsymbol{v}}{\partial\boldsymbol{x}}\right)^{\top}F+F^{\top}\frac{\partial\boldsymbol{v}}{\partial\boldsymbol{x}}\right)\right)\right)+\rho\boldsymbol{F}\\\
\frac{\partial\boldsymbol{v}}{\partial\boldsymbol{x}}\end{bmatrix}$

which coincides with (30). Furthermore, the second block row in (31) is
equivalent to

$\underbrace{\left(\frac{1}{2}F\left(\mathbb{H}\left(F^{\top}F-I_{3}\right)\right)+\frac{1}{2}F\left(\mathbb{V}\left(\left(\frac{\partial\boldsymbol{v}}{\partial\boldsymbol{x}}\right)^{\top}F+F^{\top}\frac{\partial\boldsymbol{v}}{\partial\boldsymbol{x}}\right)\right)\right)}_{=P}\boldsymbol{\nu}\Bigg{|}_{\partial\mathcal{X}}=\boldsymbol{u}_{\mathrm{m,b}}$

and, thus, coincides with the traction boundary condition.

Similarly as in section 2.1, an important property of the port-Hamiltonian
formulation (31) is that its structure encodes an energy balance for the
Hamiltonian $\mathcal{H}_{\mathrm{m}}$. The essential properties are the
identity

$\mathcal{E}_{\mathrm{m}}(\boldsymbol{z}_{\mathrm{m}})^{*}\boldsymbol{e}_{\mathrm{m}}(\boldsymbol{z}_{\mathrm{m}})=\frac{\delta\mathcal{H}_{\mathrm{m}}}{\delta\boldsymbol{z}_{\mathrm{m}}}(\boldsymbol{z}_{\mathrm{m}})$
(35)

and the fact that

$\beta_{\mathrm{m}}\left(\begin{bmatrix}\boldsymbol{\eta}_{1}\\\
\eta_{2}\end{bmatrix},\begin{bmatrix}\mathcal{J}_{\mathrm{m}}-\mathcal{R}_{\mathrm{m}}(\boldsymbol{z})\\\
-\mathcal{U}_{\mathrm{m}}\end{bmatrix}\begin{bmatrix}\boldsymbol{\eta}_{1}\\\
\eta_{2}\end{bmatrix}\right)=-\frac{1}{2}\int\limits_{\mathcal{X}}\left\langle\mathbb{V}\left(\left(\frac{\partial\boldsymbol{\eta}_{1}}{\partial\boldsymbol{x}}\right)^{\top}F+F^{\top}\frac{\partial\boldsymbol{\eta}_{1}}{\partial\boldsymbol{x}}\right),F^{\top}\frac{\partial\boldsymbol{\eta}_{1}}{\partial\boldsymbol{x}}\right\rangle_{\mathrm{F}}\,\mathrm{d}\mathcal{X}\leq
0$ (36)

holds for all sufficiently smooth
$(\boldsymbol{\eta}_{1},\eta_{2})\colon\mathcal{X}\to\mathcal{V}_{\mathrm{m}}$,
where the bilinear form $\beta_{\mathrm{m}}$ is defined via

$\displaystyle\beta_{\mathrm{m}}\left(\begin{bmatrix}\boldsymbol{\eta}_{1}\\\
\eta_{2}\end{bmatrix},\begin{bmatrix}\boldsymbol{\xi}_{1}\\\ \xi_{2}\\\
\boldsymbol{\xi}_{3}\end{bmatrix}\right)$
$\displaystyle\vcentcolon=\int\limits_{\mathcal{X}}\left(\boldsymbol{\eta}_{1}^{\top}\boldsymbol{\xi}_{1}+\left\langle\eta_{2},\xi_{2}\right\rangle_{\mathrm{F}}\right)\,\mathrm{d}\mathcal{X}+\int\limits_{\partial\mathcal{X}}\mathcal{Y}_{\mathrm{m}}(\boldsymbol{\eta}_{1},\eta_{2})^{\top}\boldsymbol{\xi}_{3}\,\mathrm{d}S\vcentcolon=\int\limits_{\mathcal{X}}\left(\boldsymbol{\eta}_{1}^{\top}\boldsymbol{\xi}_{1}+\left\langle\eta_{2},\xi_{2}\right\rangle_{\mathrm{F}}\right)\,\mathrm{d}\mathcal{X}+\int\limits_{\partial\mathcal{X}}\boldsymbol{\eta}_{1}^{\top}\boldsymbol{\xi}_{3}\,\mathrm{d}S.$

The non-positivity of the integral expression in (36) follows from the
calculation (88)–(89) in appendix A. If (35) and (36) are satisfied, we obtain
the dissipation inequality

$\displaystyle\frac{\mathrm{d}\mathcal{H}_{\mathrm{m}}}{\mathrm{d}t}$
$\displaystyle=\int\limits_{\mathcal{X}}\left(\dot{\boldsymbol{v}}^{\top}\frac{\delta\mathcal{H}_{\mathrm{m}}}{\delta\boldsymbol{v}}(\boldsymbol{z}_{\mathrm{m}})+\left\langle\dot{F},\frac{\delta\mathcal{H}_{\mathrm{m}}}{\delta
F}(\boldsymbol{z}_{\mathrm{m}})\right\rangle_{\mathrm{F}}\right)\,\mathrm{d}\mathcal{X}=\beta_{\mathrm{m}}\left(\boldsymbol{e}_{\mathrm{m}}(\boldsymbol{z}_{\mathrm{m}}),\begin{bmatrix}\mathcal{E}_{\mathrm{m}}\dot{\boldsymbol{z}}_{\mathrm{m}}\\\
0\end{bmatrix}\right)$
$\displaystyle=\beta_{\mathrm{m}}\left(\boldsymbol{e}_{\mathrm{m}}(\boldsymbol{z}_{\mathrm{m}}),\begin{bmatrix}\mathcal{J}_{\mathrm{m}}-\mathcal{R}_{\mathrm{m}}(\boldsymbol{z}_{\mathrm{m}})\\\
-\mathcal{U}_{\mathrm{m}}(\boldsymbol{z}_{\mathrm{m}})\end{bmatrix}\boldsymbol{e}_{\mathrm{m}}(\boldsymbol{z}_{\mathrm{m}})+\begin{bmatrix}\mathcal{B}_{\mathrm{m}}&0\\\
0&I_{3}\end{bmatrix}\begin{bmatrix}\boldsymbol{u}_{\mathrm{m}}\\\
\boldsymbol{u}_{\mathrm{m,b}}\end{bmatrix}\right)$
$\displaystyle=\beta_{\mathrm{m}}\left(\boldsymbol{e}_{\mathrm{m}}(\boldsymbol{z}_{\mathrm{m}}),\begin{bmatrix}\mathcal{B}_{\mathrm{m}}\boldsymbol{u}_{\mathrm{m}}\\\
\boldsymbol{u}_{\mathrm{m,b}}\end{bmatrix}\right)-\frac{1}{2}\int\limits_{\mathcal{X}}\left\langle\mathbb{V}\left(\left(\frac{\partial\boldsymbol{v}}{\partial\boldsymbol{x}}\right)^{\top}F+F^{\top}\frac{\partial\boldsymbol{v}}{\partial\boldsymbol{x}}\right),F^{\top}\frac{\partial\boldsymbol{v}}{\partial\boldsymbol{x}}\right\rangle_{\mathrm{F}}\,\mathrm{d}\mathcal{X}$
$\displaystyle=\int\limits_{\mathcal{X}}\underbrace{\rho\boldsymbol{v}^{\top}}_{=\vcentcolon\boldsymbol{y}_{\mathrm{m}}^{\top}}\boldsymbol{u}_{\mathrm{m}}\,\mathrm{d}\mathcal{X}+\int\limits_{\partial\mathcal{X}}\underbrace{\mathcal{Y}_{\mathrm{m}}(\boldsymbol{e}_{\mathrm{m}}(\boldsymbol{z}_{\mathrm{m}}))^{\top}}_{=\vcentcolon\boldsymbol{y}_{\mathrm{m,b}}^{\top}}\cdot\,\boldsymbol{u}_{\mathrm{m,b}}\,\mathrm{d}S-\frac{1}{2}\int\limits_{\mathcal{X}}\left\langle\mathbb{V}\left(\left(\frac{\partial\boldsymbol{v}}{\partial\boldsymbol{x}}\right)^{\top}F+F^{\top}\frac{\partial\boldsymbol{v}}{\partial\boldsymbol{x}}\right),F^{\top}\frac{\partial\boldsymbol{v}}{\partial\boldsymbol{x}}\right\rangle_{\mathrm{F}}\,\mathrm{d}\mathcal{X}$
$\displaystyle=\int\limits_{\mathcal{X}}\left(\rho\boldsymbol{v}^{\top}\boldsymbol{F}-\frac{1}{2}\left\langle\mathbb{V}\left(\left(\frac{\partial\boldsymbol{v}}{\partial\boldsymbol{x}}\right)^{\top}F+F^{\top}\frac{\partial\boldsymbol{v}}{\partial\boldsymbol{x}}\right),F^{\top}\frac{\partial\boldsymbol{v}}{\partial\boldsymbol{x}}\right\rangle_{\mathrm{F}}\right)\,\mathrm{d}\mathcal{X}+\int\limits_{\partial\mathcal{X}}\boldsymbol{v}^{\top}P\boldsymbol{\nu}\,\mathrm{d}S\leq\int\limits_{\mathcal{X}}\rho\boldsymbol{v}^{\top}\boldsymbol{F}\,\mathrm{d}\mathcal{X}+\int\limits_{\partial\mathcal{X}}\boldsymbol{v}^{\top}P\boldsymbol{\nu}\,\mathrm{d}S.$

This shows that the Hamiltonian $\mathcal{H}_{\mathrm{m}}$ may only change due
to viscous friction, the mechanical power
$\rho\boldsymbol{v}^{\top}\boldsymbol{F}$ caused by the body force
$\rho\boldsymbol{F}$, and due to the traction-related energy flow through the
boundary.

It remains to show that (35) and (36) indeed hold. The former property follows
directly from (32)–(34). For the latter one, we compute

$\displaystyle\beta_{\mathrm{m}}\left(\begin{bmatrix}\boldsymbol{\eta}_{1}\\\
\eta_{2}\end{bmatrix},\begin{bmatrix}\mathcal{J}_{\mathrm{m}}-\mathcal{R}_{\mathrm{m}}(\boldsymbol{z})\\\
-\mathcal{U}_{\mathrm{m}}\end{bmatrix}\begin{bmatrix}\boldsymbol{\eta}_{1}\\\
\eta_{2}\end{bmatrix}\right)$
$\displaystyle=\int\limits_{\mathcal{X}}\left(\boldsymbol{\eta}_{1}^{\top}\left(\nabla\cdot\eta_{2}+\frac{1}{2}\nabla\cdot\left(F\left(\mathbb{V}\left(\left(\frac{\partial\boldsymbol{\eta}_{1}}{\partial\boldsymbol{x}}\right)^{\top}F+F^{\top}\frac{\partial\boldsymbol{\eta}_{1}}{\partial\boldsymbol{x}}\right)\right)\right)\right)+\left\langle\eta_{2},\frac{\partial\boldsymbol{\eta}_{1}}{\partial\boldsymbol{x}}\right\rangle_{\mathrm{F}}\right)\,\mathrm{d}\mathcal{X}$
$\displaystyle\quad-\int\limits_{\partial\mathcal{X}}\boldsymbol{\eta}_{1}^{\top}\left(\eta_{2}+\frac{1}{2}F\left(\mathbb{V}\left(\left(\frac{\partial\boldsymbol{\eta}_{1}}{\partial\boldsymbol{x}}\right)^{\top}F+F^{\top}\frac{\partial\boldsymbol{\eta}_{1}}{\partial\boldsymbol{x}}\right)\right)\right)\boldsymbol{\nu}\,\mathrm{d}S$
$\displaystyle=\int\limits_{\mathcal{X}}\nabla\cdot\left(\left(\eta_{2}+\frac{1}{2}F\left(\mathbb{V}\left(\left(\frac{\partial\boldsymbol{\eta}_{1}}{\partial\boldsymbol{x}}\right)^{\top}F+F^{\top}\frac{\partial\boldsymbol{\eta}_{1}}{\partial\boldsymbol{x}}\right)\right)\right)^{\top}\boldsymbol{\eta}_{1}\right)\,\mathrm{d}\mathcal{X}$
$\displaystyle\quad-\frac{1}{2}\int\limits_{\mathcal{X}}\left\langle
F\left(\mathbb{V}\left(\left(\frac{\partial\boldsymbol{\eta}_{1}}{\partial\boldsymbol{x}}\right)^{\top}F+F^{\top}\frac{\partial\boldsymbol{\eta}_{1}}{\partial\boldsymbol{x}}\right)\right),\frac{\partial\boldsymbol{\eta}_{1}}{\partial\boldsymbol{x}}\right\rangle_{\mathrm{F}}\,\mathrm{d}\mathcal{X}$
$\displaystyle\quad-\int\limits_{\partial\mathcal{X}}\boldsymbol{\eta}_{1}^{\top}\left(\eta_{2}+\frac{1}{2}F\left(\mathbb{V}\left(\left(\frac{\partial\boldsymbol{\eta}_{1}}{\partial\boldsymbol{x}}\right)^{\top}F+F^{\top}\frac{\partial\boldsymbol{\eta}_{1}}{\partial\boldsymbol{x}}\right)\right)\right)\boldsymbol{\nu}\,\mathrm{d}S$
$\displaystyle=-\frac{1}{2}\int\limits_{\mathcal{X}}\left\langle
F\left(\mathbb{V}\left(\left(\frac{\partial\boldsymbol{\eta}_{1}}{\partial\boldsymbol{x}}\right)^{\top}F+F^{\top}\frac{\partial\boldsymbol{\eta}_{1}}{\partial\boldsymbol{x}}\right)\right),\frac{\partial\boldsymbol{\eta}_{1}}{\partial\boldsymbol{x}}\right\rangle_{\mathrm{F}}\,\mathrm{d}\mathcal{X}=-\frac{1}{2}\int\limits_{\mathcal{X}}\left\langle\mathbb{V}\left(\left(\frac{\partial\boldsymbol{\eta}_{1}}{\partial\boldsymbol{x}}\right)^{\top}F+F^{\top}\frac{\partial\boldsymbol{\eta}_{1}}{\partial\boldsymbol{x}}\right),F^{\top}\frac{\partial\boldsymbol{\eta}_{1}}{\partial\boldsymbol{x}}\right\rangle_{\mathrm{F}}\,\mathrm{d}\mathcal{X},$

which holds for all sufficiently smooth
$(\boldsymbol{\eta}_{1},\eta_{2})\colon\mathcal{X}\to\mathcal{V}_{\mathrm{m}}$.
Here, we used the product rule (87) in appendix A for the second equality.

For practical calculations, often linearized equations of motion are
considered by employing an assumption of small displacements or strains. In
the context of static components like the stator in an electric machine, the
small displacements are considered with respect to the reference
configuration. On the other hand, when considering rotating components like
the rotor of an electric machine, the small discplacements are considered with
respect to a rigid body rotation. Both cases are addressed separately in
sections 3.1.1 and 3.1.2. Furthermore, in section 3.1.3 we address previous
studies applying model order reduction for mechanical simulations in the
context of electrical machines.

####  3.1.1 Dynamics of Non-Rotating Components

The constitutive law (28) especially satisfies $P_{\mathrm{c}}(I_{3},0)=0$,
i.e., the trivial motion $\boldsymbol{\phi}(t,\boldsymbol{x})=\boldsymbol{x}$
has zero stress. Consequently, a particular solution of (25) is given by
$\boldsymbol{\phi}(t,\boldsymbol{x})=\boldsymbol{x}$ which corresponds to the
external force $\boldsymbol{F}=0$. In many works concerned with a mechanical
analysis of the static components of an electrical machine, the displacements
are assumed to be small and, thus, linearized equations of motions are
considered, see e.g. [SBK+16, Wer17]. More precisely, the governing equations
(25) are linearized around the solution
$\boldsymbol{\phi}(t,\boldsymbol{x})=\boldsymbol{x}$. This is done by
substituting the ansatz
$\boldsymbol{\phi}(t,\boldsymbol{x})=\boldsymbol{x}+\epsilon\boldsymbol{u}(t,\boldsymbol{x})$,
$\boldsymbol{F}(t,\boldsymbol{x})=\epsilon\widecheck{\boldsymbol{F}}(t,\boldsymbol{x})$,
differentiating the resulting equation with respect to $\epsilon$, and setting
$\epsilon=0$. The resulting linearized equation reads

$\rho(\boldsymbol{x})\frac{\partial^{2}\boldsymbol{u}}{\partial
t^{2}}(t,\boldsymbol{x})=\nabla\cdot\left(\widetilde{\mathbb{H}}\frac{\partial\boldsymbol{u}}{\partial\boldsymbol{x}}+\widetilde{\mathbb{V}}\frac{\partial^{2}\boldsymbol{u}}{\partial
t\partial\boldsymbol{x}}\right)(t,\boldsymbol{x})+\rho(\boldsymbol{x})\widecheck{\boldsymbol{F}}(t,\boldsymbol{x}),$
(37)

where $\widetilde{\mathbb{H}}$ and $\widetilde{\mathbb{V}}$ are given by

$\widetilde{\mathbb{H}}=\lambda_{1}I_{3}\otimes
I_{3}+\lambda_{2}(\mathbb{I}+\mathbb{T})\quad\text{and}\quad\widetilde{\mathbb{V}}=\left(\zeta_{1}-\frac{2}{3}\zeta_{2}\right)I_{3}\otimes
I_{3}+\eta(\mathbb{I}+\mathbb{T})$ (38)

and $\mathbb{T}\in\mathbb{R}^{3,3,3,3}$ denotes the unique fourth-order tensor
satisfying $\mathbb{T}A=A^{\top}$ for all $A\in\mathbb{R}^{3,3}$. This
linearization procedure is for instance illustrated in [MH83, p.9 f.] for the
case without dissipation, i.e., $\zeta_{1}=\zeta_{2}=0$. The major model
assumptions leading to (37) may be summarized as follows:

  * •

The material may be described by a constitutive law of the form (28) with
$\mathbb{H}$ and $\mathbb{V}$ as in (29).

  * •

The displacement $\boldsymbol{\phi}-\boldsymbol{x}$ and the external force
$\boldsymbol{F}$ are assumed to be small.

After semi-discretization in space by a Galerkin projection based on finite
element basis functions, the semi-discretized governing equations are of the
form

$M\ddot{\boldsymbol{s}}(t)+D\dot{\boldsymbol{s}}(t)+K\boldsymbol{s}(t)=\boldsymbol{f}(t),$
(39)

where $\boldsymbol{s}(t)\in\mathbb{R}^{N}$ denotes the vector of displacements
at time $t$, $M\in\mathbb{R}^{N,N}$ is the mass matrix, $D\in\mathbb{R}^{N,N}$
the damping matrix, $K\in\mathbb{R}^{N,N}$ the stiffness matrix, and
$\boldsymbol{f}(t)\in\mathbb{R}^{N}$ the vector of external forces at time
$t$, cf. [MC02].

Instead of deriving (39) from a semi-discretization of a partial differential
equation, it is also possible to arrive at a system like (39) by modeling the
flexible body by a mass-spring-damper system, see for instance [CP03, Ch. 9].

Sometimes, researchers also follow a hybrid approach, where the mass and the
stiffness matrix are derived from a finite element discretization and the
damping matrix is for instance obtained by a surrogate data based ansatz, or
the assumption of Rayleigh damping, i.e., as a linear combination of the mass
and the stiffness matrix, cf. [Bat14, sec. 9.3.3]. The coefficients of this
linear combination may for example be tuned based on data from a digital twin
and updated during the lifetime of the machine, see for instance [KPW21]. It
is also common to completely neglect damping effects and consider the undamped
system

$M\ddot{\boldsymbol{s}}(t)+K\boldsymbol{s}(t)=\boldsymbol{f}(t),$

see e.g. [HSHS92] or [Bat14, sec. 9.3.2].

####  3.1.2 Dynamics of Rotating Components

In the case of rotating components like the rotor of an electrical machine,
the assumption of small displacements as discussed in the last subsection is
often not valid, since the rotation leads to large displacements.
Nevertheless, in a co-rotating frame, the displacements may still be
considered to be small. In the following, we assume the case of a rigid-body
rotation in the $x_{1}x_{2}$-plane with constant rotation velocity $\omega$.
This rotation may be described by the rotation matrix

$H(t)=\begin{bmatrix}\cos(\omega t)&-\sin(\omega t)&0\\\ \sin(\omega
t)&\cos(\omega t)&0\\\ 0&0&1\end{bmatrix}.$ (40)

Especially, we note that the constitutive law (28) satisfies
$P_{\mathrm{c}}(H(t),\dot{H}(t))=0$ for all $t\geq 0$, i.e., the rigid body
rotation described by $\boldsymbol{\phi}(t,\boldsymbol{x})=H(t)\boldsymbol{x}$
is stress-free. This follows from (28) and the facts that $H(t)$ is orthogonal
and

$\dot{H}(t)^{\top}H(t)=\omega\underbrace{\begin{bmatrix}0&1&0\\\ -1&0&0\\\
0&0&0\end{bmatrix}}_{=\vcentcolon J}$

is skew-symmetric for any $t\geq 0$. Consequently, a particular solution of
(25) is given by $\boldsymbol{\phi}(t,\boldsymbol{x})=H(t)\boldsymbol{x}$
which corresponds to the external force
$\boldsymbol{F}(t,\boldsymbol{x})=\ddot{H}(t)\boldsymbol{x}$. Substituting the
ansatz
$\boldsymbol{\phi}(t,\boldsymbol{x})=H(t)(\boldsymbol{x}+\epsilon\boldsymbol{u}(t,\boldsymbol{x}))$,
$\boldsymbol{F}(t,\boldsymbol{x})=\ddot{H}(t)\boldsymbol{x}+\epsilon\widecheck{\boldsymbol{F}}(t,\boldsymbol{x})$,
differentiating the resulting equation with respect to $\epsilon$, setting
$\epsilon=0$, and multiplying from the left by $H(t)^{\top}$ leads to the
linearized equation

$\rho(\boldsymbol{x})\left(\frac{\partial^{2}\boldsymbol{u}}{\partial
t^{2}}(t,\boldsymbol{x})+2\omega J\frac{\partial\boldsymbol{u}}{\partial
t}(t,\boldsymbol{x})-\omega^{2}\boldsymbol{u}(t,\boldsymbol{x})\right)=\nabla\cdot\left(\widetilde{\mathbb{H}}\frac{\partial\boldsymbol{u}}{\partial\boldsymbol{x}}+\widetilde{\mathbb{V}}\frac{\partial^{2}\boldsymbol{u}}{\partial
t\partial\boldsymbol{x}}\right)(t,\boldsymbol{x})+\rho(\boldsymbol{x})H(t)^{\top}\widecheck{\boldsymbol{F}}(t,\boldsymbol{x})$
(41)

with $\widetilde{\mathbb{H}}$ and $\widetilde{\mathbb{V}}$ as defined in (38).
A similar linearization approach has been considered in [GMM07] in the context
of linear elasticity. Similarly as in the previous subsection, the major model
assumptions leading to (41) may be summarized as follows:

  * •

The material may be described by a constitutive law of the form (28) with
$\mathbb{H}$ and $\mathbb{V}$ as in (29).

  * •

The relative displacement $\boldsymbol{\phi}-H\boldsymbol{x}$ and the relative
external force $\boldsymbol{F}-\ddot{H}\boldsymbol{x}$ are assumed to be
small.

Semi-discretization in space via a Galerkin finite element method leads to the
semi-discretized system

$M\ddot{\boldsymbol{s}}(t)+(D+2\omega
G)\dot{\boldsymbol{s}}(t)+(K-\omega^{2}Z)\boldsymbol{s}(t)=\boldsymbol{f}(t),$
(42)

where $\boldsymbol{s}(t)\in\mathbb{R}^{N}$ is the displacement vector in the
co-rotating coordinate system at time $t$, $G\in\mathbb{R}^{N,N}$ the
gyroscopic matrix, $Z\in\mathbb{R}^{N,N}$ the centrifugal stiffness matrix,
and $\boldsymbol{f}(t)\in\mathbb{R}^{N}$ the vector of external forces in the
co-rotating coordinate system at time $t$, cf. [Kir16] for an alternative
derivation via an Euler–Lagrange approach. Also, a lumped-parameter modeling
approach may lead to a similar equation system as the one in (42), cf. [Gen05,
ch. 4]. Alternatively, one can use a data-driven method to estimate the values
of the parameters in the lumped-parameter model as in [DGB19, ch. 3].

An effect which is not captured by the linearized theory leading to (42) is
the so-called centrifugal stiffening effect, see for instance [SVQ87]. This
effect leads to a higher bending stiffness of a rotating beam in comparison to
a non-rotating beam, cf. [BS02]. To take this effect into account, usually a
term of the form $\omega^{2}K_{\mathrm{G}}\boldsymbol{s}(t)$ with
$K_{\mathrm{G}}\in\mathbb{R}^{N,N}$ is added to the left-hand side of (42),
see [Kir16, SVQ87].

In [CG17], the authors investigated electromagnetic noise within an electrical
machine by considering a model of the form (42) for the rotor, with an
additional geometric stiffness matrix and without viscous damping effects,
i.e., they neglected the $D$ term.

In [ANS10], the authors investigate the stability of a whirling rotor by a
Jeffcott rotor model with one complex-valued degree of freedom or,
equivalently, with two real-valued degrees of freedom of the form

$m\ddot{\boldsymbol{s}}(t)+d\dot{\boldsymbol{s}}(t)+k\boldsymbol{s}(t)=\boldsymbol{f}(t,\boldsymbol{s}),$

where $m$ is the rotor mass, $d$ the damping coefficient, $k$ the stiffness
constant, $\boldsymbol{f}$ the time- and state-dependent external forces, and
$\boldsymbol{s}(t)\in\mathbb{R}^{2}$ the coordinate vector of the displacement
at time $t$, see also [Gen05, sec. 2.4]. This model does neither take
structural damping into account nor gyroscopic forces. The force vector
$\boldsymbol{f}$ includes forces caused by the eccentricity of the rotor and
electromagnetic forces. The latter ones are dependent on the displacements and
in [ANS10] the authors calculate these forces by adding two further complex-
valued differential equations to the system. The resulting overall
differential equation system is linear and the stability is investigated by
computing the poles of the corresponding transfer function.

####  3.1.3 Model Order Reduction for the Mechanical Simulation

Especially, when using full 3D models for the mechanical behavior of the rotor
and the stator, the simulation times are very high and prevent any real-time
applications. As pointed out in the last two sections, then model order
reduction is a useful tool to reduce the computational burden. In [SBB14], the
authors use a Krylov subspace method to obtain a reduced-order model (ROM) for
simulating the stator vibrations within an electric motor. To this end, the
authors consider the mapping from an input vector which corresponds to
external forces to an output vector which coincides with the displacement
vector. The Arnoldi algorithm, see e.g. [LSY98], is then applied to the
transfer function of the undamped system to determine a suitable Krylov
subspace and a corresponding projection matrix, which is afterwards used to
obtain the reduced-order model via projection of the original damped system.
The article also compares this MOR approach with the classical mode
superposition method, where the projection matrix is based on the eigenmodes,
and observe a superior performance of the reduced-order model obtained by the
Krylov technique. Similarly, in [EGM+19] a Krylov subspace method has been
proposed for constructing a reduced-order model for investigating the noise
radiated from the stator. In contrast to [SBB14], they consider only a few
pre-defined displacements as output vector and they apply the Krylov technique
directly to the damped system by employing a second-order Arnoldi algorithm. A
model reduction scheme for the rotor dynamics has been proposed in [Kot19]
based on the POD method, where the POD basis captures the dynamics for a
certain range of the rotor speed $\omega$. For further model reduction
approaches for mechanical systems in general, we refer to the references in
[BBD+19, LZZ+21].

The simulation of the mechanical full-order model (FOM) is often based on
commercial software. In particular, Ansys Mechanical is frequently used, cf.
[BSKDD10, EGM+19, Kot19], but also others such as NASTRAN, e.g. in [PLF12], or
Abaqus FEA, cf. [CG17]. In [SS18], where also the coupling to the thermal
field is considered, the authors use the software COMSOL Multiphysics to
investigate the stresses caused by thermal effects. In some works, open source
software packages are employed, as for instance Elmer which is used in
[SBK+16]. Some of the mentioned commercial software packages have also built-
in model reduction tools. However, to allow for more flexibility in the choice
of the MOR scheme and to be able to set up a model hierarchy including error
bounds, it may be more effective to export the coefficient matrices from the
software and, afterwards, perform MOR based on these extracted matrices. For
Ansys, this extraction and subsequent MOR has been, for example, demonstrated
in [RK06] and for Abaqus in [Her08]. Furthermore, the software LiveLink allows
to export the system matrices from COMSOL Multiphysics to MATLAB, which has
been exploited in [HDH+21] for model order reduction.

###  3.2 Calculation of the Rotational Speed

In section 3.1.2 we have considered the case of a rotor with constant
rotational speed. Especially, in the start-up and in the shut-down phase as
well as when passing from one operation speed to another, the rotation speed
is not constant, but instead a result of an equilibrium of torques. A simple
model for simulating the rotational speed of the rotor shaft based on the
electromagnetic torque and the load torque is given by

$\Theta\dot{\omega}(t)=-\mu_{\mathrm{f}}\omega(t)+T_{\mathrm{e}}(t)-T_{\mathrm{L}}(t),$
(43)

where $\omega$ denotes the motor speed, $\Theta$ the torque of inertia,
$\mu_{\mathrm{f}}$ the friction coefficient, $T_{\mathrm{e}}$ the
electromechanical torque, and $T_{\mathrm{L}}$ the load torque, see for
instance [DHK21]. Neglecting the friction in (43) leads to the simplified
undamped model

$\Theta\dot{\omega}(t)=T_{\mathrm{e}}(t)-T_{\mathrm{L}}(t),$ (44)

cf. [HRM20, sec. 2.20.3].

##  4 Thermal Simulation

In this section the simulation of the temperature distribution within an
electrical machine is discussed. There are several reasons which lead to a
need for temperature simulations in electrical machines. For instance, high
temperatures usually lead to high resistances in copper which result in turn
in an increase of copper losses and, thus, it has also an effect on the
overall efficiency. Also ageing and damage of the windings can be facilitated
by an excessive thermal load, cf. [YBK+17]. Thus, a thermal analysis by means
of numerical simulations is helpful for a suitable design of the cooling
system which prevents excessively high temperatures in the machine. A higher
temperature than expected could also be a sign of a fault in the system and
thus tracking the temperature trend in the machine can help in predictive
maintenance.

The most relevant approaches for theremal simulation can roughly be subdivided
into two classes, see for instance [Gra92, YBK+17]. Approaches which are based
on the numerical discretization of the heat conduction equation are discussed
in section 4.1, whereas so-called _thermal equivalent circuits_ , which are
based on lumped-parameter representations, are addressed in section 4.2.

###  4.1 Energy Balance and Numerical Discretization Schemes

The time evolution of the temperature field is primarily governed by the
balance of energy, which reads in general

$\rho(\boldsymbol{x})\frac{\partial u}{\partial
t}(t,\boldsymbol{x})=-\nabla\cdot\boldsymbol{q}(t,\boldsymbol{x})+\frac{1}{2}\left\langle
S(t,\boldsymbol{x}),\frac{\partial}{\partial
t}(F^{\top}F)(t,\boldsymbol{x})\right\rangle_{\mathrm{F}}+\rho(\boldsymbol{x})r(t,\boldsymbol{x}).$
(45)

Here, $u\colon\mathbb{R}_{\geq 0}\times\mathcal{X}\to\mathbb{R}$ denotes the
specific internal energy, $\boldsymbol{q}\colon\mathbb{R}_{\geq
0}\times\mathcal{X}\to\mathbb{R}^{3}$ the heat flux vector,
$S\vcentcolon=F^{-1}P$ the second Piola–Kirchhoff stress tensor, $F$ the
deformation tensor introduced in the previous section,
$\left\langle\cdot,\cdot\right\rangle_{\mathrm{F}}$ the Frobenius inner
product, and $r\colon\mathbb{R}_{\geq 0}\times\mathcal{X}\to\mathbb{R}$ the
specific heat supply, cf. [MH83, p. 145]. Moreover, one can show the identity

$\frac{1}{2}\left\langle S,\frac{\partial}{\partial
t}(F^{\top}F)\right\rangle_{\mathrm{F}}=\left\langle P,\frac{\partial
F}{\partial t}\right\rangle_{\mathrm{F}},$ (46)

where $P$ denotes the first Piola–Kirchhoff stress tensor as introduced in the
previous section, cf. [MH83, p. 144]. Furthermore, the specific internal
energy $u$ may be expressed in terms of the Helmholtz free energy $\psi$, the
temperature $\vartheta$, and the specific entropy $s$ via

$u=\psi+\vartheta s.$ (47)

The system is then closed via constitutive equations for $\psi$, $s$, and
$\boldsymbol{q}$, whereas the temperature $\vartheta$ is the unknown of the
energy equation. To be consistent with the considerations from section 3.1, we
consider again a Kelvin–Voigt material which results in a constitutive law
$\psi=\psi_{\mathrm{c}}(\vartheta,F)$ of the form

$\psi_{\mathrm{c}}(\vartheta,F)=\psi_{0}(\vartheta)+\frac{1}{8\rho}\left\langle\mathbb{H}(F^{\top}F-I_{3}),F^{\top}F-I_{3}\right\rangle_{\mathrm{F}},$
(48)

where $\psi_{0}$ is a given function of the temperature representing the free
energy in a deformation-free state and $\mathbb{H}$ is the elasticity tensor
as defined in (29), cf. [GM21, sec. 6.3]. The specific entropy can be derived
from the free energy via the relation

$s=-\frac{\partial\psi_{\mathrm{c}}}{\partial\vartheta}=-\frac{\partial\psi_{0}}{\partial\vartheta},$
(49)

cf. [GM21, sec. 3]. Furthermore, for the heat flux $\boldsymbol{q}$ we use
Fourier’s law

$\boldsymbol{q}=-\kappa\nabla\vartheta$ (50)

where $\kappa\in\mathbb{R}^{3,3}$ is the thermal conductivity matrix, cf.
[Dim11, sec. 7.1.2]. Finally, by substituting (47),
$\psi=\psi_{\mathrm{c}}(\vartheta,F)$, (49), (50), (46), and
$P=P_{\mathrm{c}}(F,\frac{\partial F}{\partial t})$ into (45) and by using the
constitutive equations (48) and (28), we arrive at the heat conduction
equation

$\rho(\boldsymbol{x})c(\vartheta(t,\boldsymbol{x}))\frac{\partial\vartheta}{\partial
t}(t,\boldsymbol{x})=\nabla\cdot(\kappa\nabla\vartheta)(t,\boldsymbol{x})+\left\langle\widetilde{\mathbb{V}}\left(F(t,\boldsymbol{x})^{\top}\frac{\partial
F}{\partial
t}(t,\boldsymbol{x})\right),F(t,\boldsymbol{x})^{\top}\frac{\partial
F}{\partial
t}(t,\boldsymbol{x})\right\rangle_{\mathrm{F}}+\rho(\boldsymbol{x})r(t,\boldsymbol{x}),$
(51)

where

$c(\vartheta)\vcentcolon=-\vartheta\psi_{0}^{\prime\prime}(\vartheta)$ (52)

is the specific heat at constant strain, cf. [BH93, eq. (8)], and
$\widetilde{\mathbb{V}}$ is as defined in (38). The model assumptions leading
to (51) are essentially the constitutive laws (28) and (48)–(50) with
$\mathbb{H}$ and $\mathbb{V}$ as in (29).

The heat source due to mechanical friction enters (51) via the second term on
the right-hand side and other heat sources enter the equation via the $r$
term, which may in general also depend on the temperature. For instance, when
an electric current flows through the material, this is usually accompanied by
Joule heating, which corresponds to a heat source of the form

$\rho r_{\mathrm{J}}=\boldsymbol{J}\cdot\boldsymbol{E},$

where $\boldsymbol{J}$ denotes the current density and $\boldsymbol{E}$ the
electric field, cf. section 2. In many cases, a linear relationship between
$\boldsymbol{J}$ and $\boldsymbol{E}$ may be assumed leading to
$\widehat{\boldsymbol{J}}=\sigma\widehat{\boldsymbol{E}}$, where $\sigma$
denotes the electrical conductivity, see for instance [LL60, § 20]. Another
origin of heat sources is a time-varying magnetic field which leads e.g. to
hysteresis losses, cf. [Sud14, sec. 6.2], [EM20, sec. 10.3.1], [Fer16, sec.
10.3], and eddy current losses, cf. [Sud14, sec. 6.1]. For a more extensive
discussion of the different origins of losses within electrical machines, we
refer to [YBK+17, sec. 4].

The partial differential equation (51) is closed by an initial condition of
the form

$\vartheta(0,\boldsymbol{x})=\vartheta_{0}(\boldsymbol{x})\quad\text{for
}\boldsymbol{x}\in\mathcal{X}$

with initial temperature profile $\vartheta_{0}$ and by boundary conditions of
the form

$\displaystyle\vartheta(t,\boldsymbol{x})$
$\displaystyle=\vartheta_{\mathrm{b}}(t,\boldsymbol{x})\quad\text{for
}\boldsymbol{x}\in\Gamma_{1},$ (53)
$\displaystyle\boldsymbol{q}(t,\boldsymbol{x})^{\top}\boldsymbol{\nu}(\boldsymbol{x})$
$\displaystyle=q_{\mathrm{b}}(t,\boldsymbol{x})\quad\text{for
}\boldsymbol{x}\in\Gamma_{2},$ (54) $\displaystyle
w_{1}\boldsymbol{q}(t,\boldsymbol{x})^{\top}\boldsymbol{\nu}(\boldsymbol{x})+w_{2}\vartheta(t,\boldsymbol{x})$
$\displaystyle=\gamma_{\mathrm{b}}(t,\boldsymbol{x})\quad\text{for
}\boldsymbol{x}\in\Gamma_{3},$ (55)

with $\Gamma_{1}\cup\Gamma_{2}\cup\Gamma_{3}=\partial\mathcal{X}$,
$\Gamma_{i}\cap\Gamma_{j}=\emptyset$ for $i\neq j$, prescribed boundary
temperature $\vartheta_{\mathrm{b}}$, prescribed boundary heat flow
$q_{\mathrm{b}}$, and outer unit-normal vector field $\boldsymbol{\nu}$ on the
boundary $\partial\mathcal{X}$, cf. [MH83, p. 211]. Furthermore, the Robin
boundary condition (55) is determined by the prescribed function
$\gamma_{\mathrm{b}}$ and by the weighting parameters
$w_{1},w_{2}\in\mathbb{R}$. Typical in the context of electrical machines are
Neumann boundary conditions of the form (54) and Robin boundary conditions of
the form (55), see for instance [CZTF00, FZQ+19]. One approach for simulating
the temperature distribution within an electrical machine is based on
discretizing the heat conduction equation. The space discretization is usually
performed via the finite element method, cf. [FZQ+19, CZTF00]. In [ZCC18] the
authors use a coupled FEM-circuit model for describing the heat transfer
within an electrical machine, i.e., the finite element method is used for some
components of the electrical machine and a thermal equivalent circuit, cf.
section 4.2, for other components.

Since the heat transfer from the motor components to the surrounding coolant
depends also on the fluid flow of the coolant, numerical simulation of the
heat transfer is often also coupled to a simulation of the coolant’s fluid
dynamics. For instance, in [KHH+09] the authors use a computational fluid
dynamics (CFD) approach for simulating the coupled heat and fluid flow under
steady state conditions. The coupling of steady-state heat transfer, fluid
flow, and electromagnetics was addressed in [ZRH+12] using the finite element
method.

We close this subsection by presenting a port-Hamiltonian formulation of the
heat conduction equation (51) with a prescribed boundary value
$u_{\mathrm{th,b}}=q_{\mathrm{b}}$ for the heat flux as in (54). To this end,
we consider the Hamiltonian

$\mathcal{H}_{\mathrm{th}}(\vartheta)\vcentcolon=\int\limits_{\mathcal{X}}\rho\left(\psi_{0}(\vartheta)-\vartheta\psi_{0}^{\prime}(\vartheta)\right)\,\mathrm{d}\mathcal{X},$
(56)

which represents the internal energy of the system, cf. (47)–(49). Then, the
heat conduction equation (51) together with the heat flux boundary condition
may be written as

$\begin{bmatrix}\mathcal{E}_{\mathrm{th}}(\vartheta)\\\
0\end{bmatrix}\dot{\vartheta}=\begin{bmatrix}\mathcal{J}_{\mathrm{th}}(\vartheta)\\\
-\mathcal{U}_{\mathrm{th}}(\vartheta)\end{bmatrix}e_{\mathrm{th}}+\begin{bmatrix}\mathcal{B}_{\mathrm{th}}&0\\\
0&1\end{bmatrix}\begin{bmatrix}\boldsymbol{u}_{\mathrm{th}}\\\
u_{\mathrm{th,b}}\end{bmatrix}.$ (57)

with $\mathcal{E}_{\mathrm{th}}(\vartheta)\vcentcolon=\rho c(\vartheta)$,
$\mathcal{U}_{\mathrm{th}}(\vartheta)\vcentcolon=-\left.\boldsymbol{\nu}^{\top}\kappa\nabla\vartheta\right|_{\partial\mathcal{X}}$,
$e_{\mathrm{th}}=1$,

$\mathcal{B}_{\mathrm{th}}\vcentcolon=\begin{bmatrix}1&\rho\end{bmatrix},\quad\text{and}\quad\boldsymbol{u}_{\mathrm{th}}\vcentcolon=\begin{bmatrix}\left\langle\widetilde{\mathbb{V}}\left(F^{\top}\dot{F}\right),F^{\top}\dot{F}\right\rangle_{\mathrm{F}}\\\
r\end{bmatrix}.$

Furthermore, $\mathcal{J}_{\mathrm{th}}(\vartheta)$ is a linear operator on a
space of sufficiently smooth functions mapping from $\mathcal{X}$ to
$\mathbb{R}$ and is defined via

$\displaystyle\mathcal{J}_{\mathrm{th}}(\vartheta)\eta$
$\displaystyle\vcentcolon=\nabla\cdot\left(\eta\kappa\nabla\vartheta\right)+(\kappa\nabla\vartheta)^{\top}\nabla\eta.$

The definitions of $\mathcal{U}_{\mathrm{th}}(\vartheta)$ and
$e_{\mathrm{th}}$ as well as (50) immediately yield that the second block row
in (57) coincides with the boundary condition for the heat flux specified in
(54). Also the equivalence of the heat conduction equation (51) and the first
block row of (57) follows directly from the definitions of
$\mathcal{E}_{\mathrm{th}}(\vartheta)$,
$\mathcal{J}_{\mathrm{th}}(\vartheta)$, $\mathcal{B}_{\mathrm{th}}$, and
$\boldsymbol{u}_{\mathrm{th}}$.

Similarly as in the previous sections, the structure of the port-Hamiltonian
formulation (57) encodes an energy balance for the Hamiltonian
$\mathcal{H}_{\mathrm{th}}$. The essential properties are the identity

$\mathcal{E}_{\mathrm{th}}(\vartheta)e_{\mathrm{th}}=\frac{\delta\mathcal{H}_{\mathrm{th}}}{\delta\vartheta}(\vartheta)$
(58)

and the fact that

$\beta_{\mathrm{th}}\left(\eta,\begin{bmatrix}\mathcal{J}_{\mathrm{th}}(\vartheta)\\\
-\mathcal{U}_{\mathrm{th}}(\vartheta)\end{bmatrix}\eta\right)=0$ (59)

holds for all sufficiently smooth $\eta\colon\mathcal{X}\to\mathbb{R}$, where
the bilinear form $\beta_{\mathrm{th}}$ is defined via

$\displaystyle\beta_{\mathrm{th}}\left(\eta,\begin{bmatrix}\xi_{1}\\\
\xi_{2}\end{bmatrix}\right)$
$\displaystyle\vcentcolon=\int\limits_{\mathcal{X}}\eta\xi_{1}\,\mathrm{d}\mathcal{X}+\int\limits_{\partial\mathcal{X}}\mathcal{Y}_{\mathrm{th}}(\eta)\xi_{2}\,\mathrm{d}S\vcentcolon=\int\limits_{\mathcal{X}}\eta\xi_{1}\,\mathrm{d}\mathcal{X}-\int\limits_{\partial\mathcal{X}}\eta\xi_{2}\,\mathrm{d}S.$

If (58) and (59) are satisfied, we obtain the dissipation inequality

$\displaystyle\frac{\mathrm{d}\mathcal{H}_{\mathrm{th}}}{\mathrm{d}t}$
$\displaystyle=\int\limits_{\mathcal{X}}\dot{\vartheta}\frac{\delta\mathcal{H}_{\mathrm{th}}}{\delta\vartheta}(\vartheta)\,\mathrm{d}\mathcal{X}=\beta_{\mathrm{th}}\left(e_{\mathrm{th}},\begin{bmatrix}\mathcal{E}_{\mathrm{th}}(\vartheta)\dot{\vartheta}\\\
0\end{bmatrix}\right)=\beta\left(e_{\mathrm{th}},\begin{bmatrix}\mathcal{J}_{\mathrm{th}}(\vartheta)\\\
-\mathcal{U}_{\mathrm{th}}(\vartheta)\end{bmatrix}e_{\mathrm{th}}+\begin{bmatrix}\mathcal{B}_{\mathrm{th}}&0\\\
0&1\end{bmatrix}\begin{bmatrix}\boldsymbol{u}_{\mathrm{th}}\\\
u_{\mathrm{th,b}}\end{bmatrix}\right)=\beta_{\mathrm{th}}\left(e_{\mathrm{th}},\begin{bmatrix}\mathcal{B}_{\mathrm{th}}\boldsymbol{u}_{\mathrm{th}}\\\
u_{\mathrm{th,b}}\end{bmatrix}\right)$
$\displaystyle=\int\limits_{\mathcal{X}}\underbrace{\begin{bmatrix}1&\rho\end{bmatrix}}_{=\vcentcolon\boldsymbol{y}_{\mathrm{th}}^{\top}}\boldsymbol{u}_{\mathrm{th}}\,\mathrm{d}\mathcal{X}+\int\limits_{\partial\mathcal{X}}\underbrace{\mathcal{Y}_{\mathrm{th}}(e_{\mathrm{th}})}_{=\vcentcolon
y_{\mathrm{th,b}}}\cdot\,u_{\mathrm{th,b}}\,\mathrm{d}S=\int\limits_{\mathcal{X}}\left(\left\langle\widetilde{\mathbb{V}}\left(F^{\top}\dot{F}\right),F^{\top}\dot{F}\right\rangle_{\mathrm{F}}+\rho
r\right)\,\mathrm{d}\mathcal{X}+\int\limits_{\partial\mathcal{X}}\boldsymbol{\nu}^{\top}\kappa\nabla\vartheta\,\mathrm{d}S.$

This shows that the internal energy may only change due to viscous friction,
the heat source $\rho r$, and due to the heat flux through the boundary.

It remains to show that (58) and (59) indeed hold. To this end, we first
compute the variational derivative of the Hamiltonian which is given by

$\frac{\delta\mathcal{H}_{\mathrm{th}}}{\delta\vartheta}(\vartheta)=\rho\left(\psi_{0}^{\prime}(\vartheta)-\psi_{0}^{\prime}(\vartheta)-\vartheta\psi_{0}^{\prime\prime}(\vartheta)\right)=\rho
c(\vartheta).$

From this we immediately obtain that (58) is satisfied. To verify (59), we
compute

$\displaystyle\beta_{\mathrm{th}}\left(\eta,\begin{bmatrix}\mathcal{J}_{\mathrm{th}}(\vartheta)\\\
-\mathcal{U}_{\mathrm{th}}(\vartheta)\end{bmatrix}\eta\right)$
$\displaystyle=\int\limits_{\mathcal{X}}\eta\left(\nabla\cdot\left(\eta\kappa\nabla\vartheta\right)+(\kappa\nabla\vartheta)^{\top}\nabla\eta\right)\,\mathrm{d}\mathcal{X}-\int\limits_{\partial\mathcal{X}}\eta^{2}\boldsymbol{\nu}^{\top}\kappa\nabla\vartheta\,\mathrm{d}S$
$\displaystyle=\int\limits_{\mathcal{X}}\nabla\cdot\left(\eta^{2}\kappa\nabla\vartheta\right)\,\mathrm{d}\mathcal{X}-\int\limits_{\partial\mathcal{X}}\eta^{2}\boldsymbol{\nu}^{\top}\kappa\nabla\vartheta\,\mathrm{d}S=0,$

which holds for all sufficiently smooth $\eta\colon\mathcal{X}\to\mathbb{R}$.

###  4.2 Thermal Equivalent Circuits

For practical temperature calculations for electrical machines, lumped-
parameter thermal networks are often used instead of an explicit spatial
discretization of the heat conduction equation, see for instance [BCS+09,
CHM16, DF04, GHGZ16, LS69, MRT91, PTP69] and the references therein. The main
idea of the thermal equivalent circuit modeling approach is to replace the
spatially distributed heat conduction equation by a discrete network of
temperature nodes. This is done by decomposing the spatial domain into several
subdomains and assuming a homogeneous temperature for each of these
subdomains. The resulting thermal equivalent circuit model shares analogies
with electrical resistor-capacitor circuits, where heat transfer between the
subdomains is modeled via a conductor including a thermal resistor and the
heat storage within each subdomain is described by a thermal capacitance.
Furthermore, there is a heat source term in each subdomain which corresponds
to a current source in the electrical circuit analogy, see for instance
[DT17]. A detailed derivation of thermal equivalent circuits based on the heat
conduction equation is for instance presented in [Sud14, Ch. 10].

$P_{1}$$C_{1}$$\Lambda_{10}$$P_{2}$$C_{2}$$\Lambda_{20}$$\Lambda_{12}$$\vartheta_{1}$$\vartheta_{2}$$\vartheta_{0}$$\vartheta_{0}$$\vartheta_{0}$$\vartheta_{0}$$\vartheta_{0}$$\vartheta_{0}$
Figure 5: Example of a thermal equivalent circuit with two subdomains

An example for a thermal equivalent circuit with two subdomains is depicted in
Figure 5. Each subdomain consists of a heat source term $P_{i}$ (current
source in circuit analogy), a thermal capacitance $C_{i}$ (a capacitor in
circuit analogy), and a thermal resistance $\Lambda_{i0}$ (resistor in circuit
analogy) for the heat transfer to the environment or to a cooling medium with
temperature $\vartheta_{0}$, which is assumed to be constant in the following.
The temperature $\vartheta_{i}$ of a subdomain corresponds in the circuit
analogy to the electric potential at the node where the three elements of the
subdomain are connected. Furthermore, the heat transfer between different
subdomains is also described by a thermal resistance, which is $\Lambda_{12}$
in the example from Figure 5.

The governing equations of a thermal equivalent circuit are derived by
applying Kirchhoff’s current law to the node of every subdomain. For a general
system of $N$ subdomains this yields an implicit ordinary differential
equation system of the form

$E\dot{\boldsymbol{z}}(t)=A\boldsymbol{z}(t)+\boldsymbol{P}(t)$ (60a) with
$\displaystyle E$ $\displaystyle\vcentcolon=\begin{bmatrix}C_{1}&&&\\\
&C_{2}&&\\\ &&\ddots&\\\
&&&C_{N}\end{bmatrix},\qquad\boldsymbol{z}\vcentcolon=\begin{bmatrix}\vartheta_{1}-\vartheta_{0}\\\
\vartheta_{2}-\vartheta_{0}\\\ \vdots\\\
\vartheta_{N}-\vartheta_{0}\end{bmatrix},\qquad\boldsymbol{P}=\begin{bmatrix}P_{1}\\\
P_{2}\\\ \vdots\\\ P_{N}\end{bmatrix},$ (60b) $\displaystyle A$
$\displaystyle=\Lambda_{0}+\Lambda,\qquad[\Lambda]_{ij}=\Lambda_{ij}=\Lambda_{ji}\quad\text{for
}i,j=1,\ldots,N,$ $\displaystyle\Lambda_{0}$
$\displaystyle=\begin{bmatrix}-\Lambda_{10}-\sum_{j=1}^{N}\Lambda_{1j}&&&\\\
&-\Lambda_{20}-\sum_{j=1}^{N}\Lambda_{2j}&&\\\ &&\ddots&\\\
&&&-\Lambda_{N0}-\sum_{j=1}^{N}\Lambda_{Nj}\end{bmatrix},$

cf. [Bac33]. We emphasize that this may also be formulated as a differential-
algebraic equation system if Kirchhoff’s current law is explicitly
incorporated as algebraic constraint. Furthermore, these thermal equivalent
circuit equations can be generalized by considering for instance temperature-
dependent coefficients or source terms, cf. [BK10, HRF17, Ric67].

######  Remark 4.1.

Usually, each of the $N$ subdomains is only connected by heat transfer to a
few neighboring subdomains which results in a sparse matrix $\Lambda$ and,
thus, also $A$ is sparse. In this case, the network can be characterized in a
concise manner via the incidence matrix of the corresponding graph, see for
instance [vdSM13].

The parameters determining the thermal equivalent circuit equations (60) can,
for example, be determined based on the heat transfer theory, which relies on
the geometry and material properties of the motor and possibly on empirical
correlations for the convective heat transfer, cf. [BCS05, QSDD14, SWM14].
Another approach is given by data-driven methods which estimate the parameters
based on experimental data, see e.g. [BP14, HPB14, KHL14, ZXL+19]. Also hybrid
approaches combining analytical techniques and experimental data are common,
see for instance [NGH21, SRW+10, WB15]. Alternatively, some authors propose to
use a numerical simulation based on the finite element method [AB08] or based
on the finite volume method [JTL08] for determining the coefficients of the
thermal equivalent circuit.

For a thermal steady-state analysis, the capacitances and the time-dependency
of the heat sources are neglected and the resulting steady-state equations
take the form

$0=A\boldsymbol{z}+\boldsymbol{P},$

where $A$, $\boldsymbol{z}$, and $\boldsymbol{P}$ are as defined in (60b), but
without dependency on time. Such a steady-state analysis based on thermal
equivalent circuits has for instance been considered in [BCLP03, CZC17, NRP08,
Rob69].

###  4.3 Model Order Reduction for the Thermal Simulation

When applying a numerical discretization scheme as discussed in section 4.1,
the resulting finite-dimensional model may be of very high dimension,
especially when a 3D model is considered. Depending on the complexity of the
thermal equivalent circuit, also the models arising from lumped-parameter
modeling as addressed in section 4.2 may be too expensive to be solved in real
time. To overcome these difficulties, several model order reduction techniques
have been applied in the literature. In [DdlPAO10], the authors first
construct a thermal equivalent circuit model and formulate it as an input-
output system, where the inputs correspond to the losses and the outputs are
temperatures at certain nodes of interest. Then, the system is discretized in
time and the balanced truncation MOR scheme is applied to achieve a real-time
simulation. Also in [Qi18] balanced truncation is used for reducing a thermal
equivalent circuit model of an electrical machine. A validation step
demonstrates that dimensions of the reduced order model of $10-20$ yield
sufficiently accurate results. In [BVTP00], the authors compare two different
MOR approaches.

One is a variant of the method presented in [EBH87], where the state of the
reduced-order model corresponds to an approximation of the output of the full-
order model (FOM). The reduced system matrices are determined based on a
minimization of the residual obtained when substituting the FOM output into
the state equation of the ROM. The other approach is based on a system
identification technique proposed in [PHV97], where a reduced model is
determined in modal coordinates, i.e., the reduced system of ordinary
differential equations is uncoupled. The coefficient matrices are determined
via an output error minimization based on simulated step response data. The
numerical results from [BVTP00] exhibit that the ROMs obtained by the
identification method are faster and more accurate than the ones obtained by
the method in [EBH87].

In [GCHH08], the authors consider balanced truncation and another MOR
technique which is based on pole-zero cancellations of the transfer function
of the FOM. They illustrate that both approaches yield low-dimensional and
accurate surrogate models based on numerical experiments. However, the
considered FOM has only few state variables, although the authors argue that
the framework may be also applied to systems of higher dimension. In contrast,
the authors of [ZPH13] consider as FOM a three-dimensional finite element
model, which is reduced via a modal truncation approach. To this end, they
split the set of eigenmodes into dynamic and static eigenmodes based on an
excitation criterion. Then, via neglecting the dynamics of the static
eigenmodes, they only need to compute the dynamic eigenmodes and project the
FOM onto their span. The numerical experiments indicate that a ROM of
dimension $30-40$ is sufficiently accurate and, moreover, allows for real-time
simulations.

In contrast to the mechanical simulation addressed in section 3.1.3, the
simulation of the thermal full-order models is more often done based on self-
implemented code, e.g., using MATLAB/Simulink, cf. [CZTF00, DdlPAO10, Qi18,
ZPH13]. Nevertheless, there are also dedicated software packages for creating
thermal models, see for instance [BCS+09] for an overview.

##  5 Multi-Physics Coupling and Discussion

In the previous sections, we discussed mathematical models for
electromagnetism, mechanics, and thermodynamics, separately. Subsequently, we
will address the interaction between these submodels by first summarizing and
discussing some existing works in this direction in section 5.1. Afterwards,
we present a mathematical model for the coupled system, which is based on
[Kov00], and derive a corresponding port-Hamiltonian formulation in section
5.2.

###  5.1 State of the Art and Discussion

In order to obtain an overall model, a coupling between the electromagnetic,
mechanical and thermal models is necessary [HDDGB99]. Typically only weak
coupling (serial coupling) is considered where for example the thermal model
is computed based on the results of the electromagnetic simulation [MFM20] or
the forces acting on the structural model are computed based on the
electromagnetic force densities computed from an FEM simulation [RBWS10].
These models generally do not take into account the real time coupling and
interaction between the multi-physics phenomena. At best, there is an
iteration between two domains where the results of one domain are fed to the
model of the other domain. Then the results of the latter one are fed to the
model for the first domain. This is repeated until the error level is
satisfactory. This cosimulation approach is carried out for example in [MFM20]
where the electromagnetic simulation resulting in the resistive and iron
losses are directly fed to a steady-state thermal simulation which determines
the actual operation temperature. It is then fed back to the electromagnetic
model for an update of the material properties. For copper losses, the
temperature dependence can also be directly included in the model.

In [VR03], the authors consider a weak coupling of thermal and electromagnetic
simulations using time-harmonic methods to solve the electromagnetic problem
coupled to the transient thermal simulation and in [PLF12] the authors
consider the coupling between electromagnetic and structural models. It is
suggested that the overall accuracy of these models could be improved by
considering a simultaneous coupling of the three domains together.

However, in most publications the research focus so far is mainly on a weak
coupling as pointed out in [XHC18, sec. 5.3]. Works considering strong
coupling especially in the 3D case are rare. In this direction, for several
reasons, port-Hamiltonian systems provide a good framework to solve this kind
of coupled multi-physics systems, see [MU22, OSMM01, Sch04, vSJ14] and the
references therein. Since energy (or power) is a common quantity in all
physical domains, one can model all the subsystems using energy as the common
quantity. and couple the systems via energy exchange. Furthermore, the power-
conserving interconnection of port-Hamiltonian systems results in another
port-Hamiltonian system. This means that one can model the subsystems
independently and then interconnect them and solve for the overall system.
Since pH systems are stable and passive by construction, the interconnection
in the pH way will preserve these properties.

As discussed, a compete coupling of the models for an electrical machine is
still an open problem. This is even more true in the context of the digital
twin where in addition to the requirement of having the overall coupled model,
the model needs to satisfy the restrictions as discussed in sections 1 and 2.

First of all, the computation time for simulation of the model has to be close
to real-time. This leads to the requirement of using model reduction. In this
way one ends up with a hierarchy of models that can be used based on the
application depending on the time and accuracy requirements. The coupled model
also needs to be adaptable to motor specifications and size and be feasible
for updating to changing operating conditions. An important tool in this
context is, for example, data assimilation based on the sensor data coming
from the machine. Finally, we need a complete software solution where the
coupled system can be used, reduced and visualized.

###  5.2 Port-Hamiltonian Formulation of the Coupled System

The mathematical models considered in the previous sections are either based
on purely electromagnetic or on purely mechanical and thermodynamic
considerations. In the following, we consider the completely coupled system
based on the general modeling framework presented in [Kov00, Ch. 15], see also
appendix A for more details on the derivation of the overall model. Also, the
model assumptions are addressed in appendix A and in particular they include
the omission of polarization and magnetization effects. As state variables, we
consider in the following the electric field $\boldsymbol{E}$, the magnetic
field $\boldsymbol{B}$, the velocity $\boldsymbol{v}$, the deformation tensor
$F$, and the temperature $\vartheta$, all in Lagrangian coordinates. However,
we emphasize that the choice of state variables is in general not unique and
one might also consider different sets of unknowns, for instance involving the
charge potential or the specific momentum. To summarize all state variables,
we use the notation

$\boldsymbol{z}\vcentcolon=\begin{bmatrix}\boldsymbol{E}\\\ \boldsymbol{B}\\\
\boldsymbol{v}\\\ F\\\ \vartheta\end{bmatrix}.$

The governing equations are given by the system

$\displaystyle\epsilon_{0}\frac{\partial\boldsymbol{E}}{\partial t}$
$\displaystyle=\nabla\times\left(\boldsymbol{\mathcal{H}}_{\mathrm{c}}(\boldsymbol{z})\right)-\boldsymbol{J}_{\mathrm{c}}\left(\boldsymbol{\mathcal{E}}_{\mathrm{c}}(\boldsymbol{z}),\boldsymbol{B},\frac{\partial\boldsymbol{v}}{\partial\boldsymbol{x}},F,\vartheta,\nabla\vartheta\right)=\vcentcolon\boldsymbol{f}_{1}\left(\boldsymbol{z}\right),$
(61) $\displaystyle\frac{\partial\boldsymbol{B}}{\partial t}$
$\displaystyle=-\nabla\times\left(\boldsymbol{\mathcal{E}}_{\mathrm{c}}(\boldsymbol{z})\right)=\vcentcolon\boldsymbol{f}_{2}(\boldsymbol{z}),$
$\displaystyle\rho\frac{\partial\boldsymbol{v}}{\partial t}$
$\displaystyle=\nabla\cdot
P_{\mathrm{c}}\left(\boldsymbol{E},\boldsymbol{B},\boldsymbol{v},\frac{\partial\boldsymbol{v}}{\partial\boldsymbol{x}},F\right)+\rho\boldsymbol{F}+F^{-\top}\left(\epsilon_{0}\left(\left(\frac{\partial\boldsymbol{v}}{\partial\boldsymbol{x}}\right)^{\top}F^{-\top}(\boldsymbol{E}\times\boldsymbol{B})-\boldsymbol{E}\times\boldsymbol{f}_{2}(\boldsymbol{z})\right)-\boldsymbol{f}_{1}\left(\boldsymbol{z}\right)\times\boldsymbol{B}\right)$
$\displaystyle=\vcentcolon\boldsymbol{f}_{3}\left(\boldsymbol{z}\right),$
$\displaystyle\frac{\partial F}{\partial t}$
$\displaystyle=\frac{\partial\boldsymbol{v}}{\partial\boldsymbol{x}},$
$\displaystyle\rho c(\vartheta)\frac{\partial\vartheta}{\partial t}$
$\displaystyle=\nabla\cdot\left(P_{\mathrm{c}}\left(\boldsymbol{E},\boldsymbol{B},\boldsymbol{v},\frac{\partial\boldsymbol{v}}{\partial\boldsymbol{x}},F\right)^{\top}\boldsymbol{v}-\boldsymbol{q}_{\mathrm{c}}\left(\boldsymbol{\mathcal{E}}_{\mathrm{c}}(\boldsymbol{z}),\boldsymbol{B},\frac{\partial\boldsymbol{v}}{\partial\boldsymbol{x}},F,\vartheta,\nabla\vartheta\right)-\boldsymbol{\mathcal{E}}_{\mathrm{c}}(\boldsymbol{z})\times\boldsymbol{\mathcal{H}}_{\mathrm{c}}(\boldsymbol{z})\right)+\rho\left(\boldsymbol{v}^{\top}\boldsymbol{F}+r\right)$
$\displaystyle\quad-\rho\left\langle\frac{\partial\psi_{\mathrm{c}}}{\partial
F}(F,\vartheta),\frac{\partial\boldsymbol{v}}{\partial\boldsymbol{x}}\right\rangle_{\mathrm{F}}-\boldsymbol{v}^{\top}\boldsymbol{f}_{3}\left(\boldsymbol{z}\right)-\frac{1}{\det(F)}\left(\boldsymbol{E}^{\top}F^{\top}F\boldsymbol{f}_{1}\left(\boldsymbol{z}\right)+\frac{1}{\mu_{0}}\boldsymbol{B}^{\top}F^{\top}F\boldsymbol{f}_{2}(\boldsymbol{z})\right)$
$\displaystyle-\frac{1}{2}\sum_{i,j=1}^{3}\frac{\partial v_{i}}{\partial
x_{j}}\left(\epsilon_{0}\boldsymbol{E}^{\top}\frac{\partial}{\partial
F_{i,j}}\left(\frac{F^{\top}F}{\det(F)}\right)\boldsymbol{E}+\frac{1}{\mu_{0}}\boldsymbol{B}^{\top}\frac{\partial}{\partial
F_{i,j}}\left(\frac{F^{\top}F}{\det(F)}\right)\boldsymbol{B}\right)$

with

$\displaystyle\boldsymbol{\mathcal{H}}_{\mathrm{c}}(\boldsymbol{z})$
$\displaystyle\vcentcolon=\frac{1}{\det(F)}F^{\top}\left(\frac{1}{\mu_{0}}F\boldsymbol{B}-\epsilon_{0}\boldsymbol{v}\times(F\boldsymbol{E})\right),\quad\boldsymbol{\mathcal{E}}_{\mathrm{c}}(\boldsymbol{z})\vcentcolon=\frac{1}{\det(F)}F^{\top}\left(F\boldsymbol{E}+\boldsymbol{v}\times(F\boldsymbol{B})\right),$
(62) $\displaystyle
P_{\mathrm{c}}\left(\boldsymbol{E},\boldsymbol{B},\boldsymbol{v},\frac{\partial\boldsymbol{v}}{\partial\boldsymbol{x}},F\right)$
$\displaystyle\vcentcolon=\frac{1}{2}F\left(\mathbb{H}(F^{\top}F-I_{3})+\mathbb{V}\left(\left(\frac{\partial\boldsymbol{v}}{\partial\boldsymbol{x}}\right)^{\top}F+F^{\top}\frac{\partial\boldsymbol{v}}{\partial\boldsymbol{x}}\right)\right)+\frac{1}{\det(F)}F\left(\frac{1}{\mu_{0}}\boldsymbol{B}\boldsymbol{B}^{\top}+\epsilon_{0}\boldsymbol{E}\boldsymbol{E}^{\top}\right)$
$\displaystyle\quad+\frac{1}{\det(F)}\left(\epsilon_{0}\left((F\boldsymbol{E})\times(F\boldsymbol{B})\right)\boldsymbol{v}^{\top}-\frac{1}{2}\left(\epsilon_{0}\boldsymbol{E}^{\top}F^{\top}F\boldsymbol{E}+\frac{1}{\mu_{0}}\boldsymbol{B}^{\top}F^{\top}F\boldsymbol{B}\right)I_{3}\right)F^{-\top},$

$\psi_{\mathrm{c}}$ as in (48), and the functions
$\boldsymbol{J}_{\mathrm{c}},\boldsymbol{q}_{\mathrm{c}}\colon\mathbb{R}^{3}\times\mathbb{R}^{3}\times\mathbb{R}^{3,3}\times\mathbb{R}^{3,3}\times\mathbb{R}\times\mathbb{R}^{3}\to\mathbb{R}^{3}$
are assumed to satisfy

$\boldsymbol{J}_{\mathrm{c}}\left(\boldsymbol{\eta}_{1},\boldsymbol{\eta}_{2},\eta_{3},\eta_{4},\eta_{5},\boldsymbol{\eta}_{6}\right)^{\top}\boldsymbol{\eta}_{1}-\frac{1}{\eta_{5}}\boldsymbol{q}_{\mathrm{c}}\left(\boldsymbol{\eta}_{1},\boldsymbol{\eta}_{2},\eta_{3},\eta_{4},\eta_{5},\boldsymbol{\eta}_{6}\right)^{\top}\boldsymbol{\eta}_{6}\geq
0$ (63)

for all
$(\boldsymbol{\eta}_{1},\boldsymbol{\eta}_{2},\eta_{3},\eta_{4},\eta_{5},\boldsymbol{\eta}_{6})\in\mathbb{R}^{3}\times\mathbb{R}^{3}\times\mathbb{R}^{3,3}\times\mathbb{R}^{3,3}\times\mathbb{R}_{>0}\times\mathbb{R}^{3}$
with $\det(\eta_{4})>0$. This is for instance satisfied when
$\boldsymbol{J}_{\mathrm{c}}$ is chosen proportional to its first argument
with non-negative proportionality constant and $\boldsymbol{q}_{\mathrm{c}}$
proportional to its last argument with non-positive proportionality constant,
respectively. These linear relations are often referred to as Ohm’s law and
Fourier’s law, respectively.

Associated to the governing equations we consider the Hamiltonian

$\displaystyle\mathcal{H}(\boldsymbol{E},\boldsymbol{B},\boldsymbol{v},F,\vartheta)$
$\displaystyle\vcentcolon=\int\limits_{\mathcal{X}}\left(\frac{1}{2\det(F)}\left(\epsilon_{0}\boldsymbol{E}^{\top}F^{\top}F\boldsymbol{E}+\frac{1}{\mu_{0}}\boldsymbol{B}^{\top}F^{\top}F\boldsymbol{B}\right)\right.$
(64)
$\displaystyle\quad+\left.\rho\left(\psi_{\mathrm{c}}(F,\vartheta)-\vartheta\frac{\partial\psi_{\mathrm{c}}}{\partial\vartheta}(F,\vartheta)+\frac{1}{2}\boldsymbol{v}^{\top}\boldsymbol{v}\right)\right)\,\mathrm{d}\mathcal{X},$

which represents the total energy of the system, cf. appendix A. Based on this
Hamiltonian, we may write the governing equations together with Dirichlet
boundary conditions $u_{\mathrm{b}}$ for the total energy boundary flow

$\boldsymbol{\nu}^{\top}\left(\boldsymbol{q}_{\mathrm{c}}\left(\boldsymbol{\mathcal{E}}_{\mathrm{c}}(\boldsymbol{z}),\boldsymbol{B},\frac{\partial\boldsymbol{v}}{\partial\boldsymbol{x}},F,\vartheta,\nabla\vartheta\right)+\boldsymbol{\mathcal{E}}_{\mathrm{c}}(\boldsymbol{z})\times\boldsymbol{\mathcal{H}}_{\mathrm{c}}(\boldsymbol{z})-P_{\mathrm{c}}\left(\boldsymbol{E},\boldsymbol{B},\boldsymbol{v},\frac{\partial\boldsymbol{v}}{\partial\boldsymbol{x}},F\right)^{\top}\boldsymbol{v}\right)$
(65)

as

$\begin{bmatrix}\mathcal{E}(\boldsymbol{z})\\\
0\end{bmatrix}\dot{\boldsymbol{z}}=\begin{bmatrix}\mathcal{J}(\boldsymbol{z})\\\
-\mathcal{U}(\boldsymbol{z})\end{bmatrix}\boldsymbol{e}(\boldsymbol{z})+\begin{bmatrix}\mathcal{B}(\boldsymbol{z})&0\\\
0&1\end{bmatrix}\begin{bmatrix}\boldsymbol{u}\\\ u_{\mathrm{b}}\end{bmatrix}.$
(66)

Here, $\mathcal{E}(\boldsymbol{z})$ is a linear mapping from
$\mathcal{V}\vcentcolon=\mathbb{R}^{3}\times\mathbb{R}^{3}\times\mathbb{R}^{3}\times\mathbb{R}^{3,3}\times\mathbb{R}$
to itself defined via

$\mathcal{E}(\boldsymbol{z})\vcentcolon=\begin{bmatrix}\epsilon_{0}I_{3}&0&0&0&0\\\
0&I_{3}&0&0&0\\\ 0&0&\rho I_{3}&0&0\\\ 0&0&0&\mathbb{I}&0\\\ 0&0&0&0&\rho
c(\vartheta)\end{bmatrix}.$ (67)

Furthermore, $\mathcal{J}(\boldsymbol{z})$ is a linear operator on a space of
sufficiently smooth functions mapping from $\mathcal{X}$ to $\mathcal{V}$ and
is defined via

$\displaystyle\mathcal{J}(\boldsymbol{z})$
$\displaystyle\vcentcolon=\begin{bmatrix}0&0&0&0&\boldsymbol{f}_{1}(\boldsymbol{z})\\\
0&0&0&0&\boldsymbol{f}_{2}(\boldsymbol{z})\\\
0&0&0&0&\boldsymbol{f}_{3}(\boldsymbol{z})\\\
0&0&0&0&\frac{\partial\boldsymbol{v}}{\partial\boldsymbol{x}}\\\
-\boldsymbol{f}_{1}(\boldsymbol{z})^{\top}&-\boldsymbol{f}_{2}(\boldsymbol{z})^{\top}&-\boldsymbol{f}_{3}(\boldsymbol{z})^{\top}&-\left\langle\frac{\partial\boldsymbol{v}}{\partial\boldsymbol{x}},\cdot\right\rangle_{\mathrm{F}}&\mathcal{J}_{6,6}(\boldsymbol{z})\end{bmatrix},$
$\displaystyle\mathcal{J}_{6,6}(\boldsymbol{z})\eta$
$\displaystyle\vcentcolon=\nabla\cdot\left(\eta\left(P_{\mathrm{c}}\left(\boldsymbol{E},\boldsymbol{B},\boldsymbol{v},\frac{\partial\boldsymbol{v}}{\partial\boldsymbol{x}},F\right)^{\top}\boldsymbol{v}-\boldsymbol{q}_{\mathrm{c}}\left(\boldsymbol{\mathcal{E}}_{\mathrm{c}}(\boldsymbol{z}),\boldsymbol{B},\frac{\partial\boldsymbol{v}}{\partial\boldsymbol{x}},F,\vartheta,\nabla\vartheta\right)-\boldsymbol{\mathcal{E}}_{\mathrm{c}}(\boldsymbol{z})\times\boldsymbol{\mathcal{H}}_{\mathrm{c}}(\boldsymbol{z})\right)\right)$
$\displaystyle\quad+\left(P_{\mathrm{c}}\left(\boldsymbol{E},\boldsymbol{B},\boldsymbol{v},\frac{\partial\boldsymbol{v}}{\partial\boldsymbol{x}},F\right)^{\top}\boldsymbol{v}-\boldsymbol{q}_{\mathrm{c}}\left(\boldsymbol{\mathcal{E}}_{\mathrm{c}}(\boldsymbol{z}),\boldsymbol{B},\frac{\partial\boldsymbol{v}}{\partial\boldsymbol{x}},F,\vartheta,\nabla\vartheta\right)-\boldsymbol{\mathcal{E}}_{\mathrm{c}}(\boldsymbol{z})\times\boldsymbol{\mathcal{H}}_{\mathrm{c}}(\boldsymbol{z})\right)^{\top}\nabla\eta.$

The linear operator $\mathcal{U}(\boldsymbol{z})$ takes sufficiently smooth
functions mapping from $\mathcal{X}$ to $\mathcal{V}$ and returns functions
mapping from $\partial\mathcal{X}$ to $\mathbb{R}$ via the rule

$\mathcal{U}(\boldsymbol{z})\left(\begin{bmatrix}\boldsymbol{\eta}_{1}\\\
\boldsymbol{\eta}_{2}\\\ \boldsymbol{\eta}_{3}\\\ \eta_{4}\\\
\eta_{5}\end{bmatrix}\right)\vcentcolon=\left.\eta_{5}\boldsymbol{\nu}^{\top}\left(\boldsymbol{q}_{\mathrm{c}}\left(\boldsymbol{\mathcal{E}}_{\mathrm{c}}(\boldsymbol{z}),\boldsymbol{B},\frac{\partial\boldsymbol{v}}{\partial\boldsymbol{x}},F,\vartheta,\nabla\vartheta\right)+\boldsymbol{\mathcal{E}}_{\mathrm{c}}(\boldsymbol{z})\times\boldsymbol{\mathcal{H}}_{\mathrm{c}}(\boldsymbol{z})-P_{\mathrm{c}}\left(\boldsymbol{E},\boldsymbol{B},\boldsymbol{v},\frac{\partial\boldsymbol{v}}{\partial\boldsymbol{x}},F\right)^{\top}\boldsymbol{v}\right)\right|_{\partial\mathcal{X}}.$

Besides, the linear operator $\mathcal{B}(\boldsymbol{z})$ maps from
$\mathbb{R}^{3}\times\mathbb{R}$ to $\mathcal{V}$ and is defined as

$\mathcal{B}(\boldsymbol{z})\vcentcolon=\rho\begin{bmatrix}0&0\\\ 0&0\\\
0&0\\\ 0&0\\\ \boldsymbol{v}^{\top}&1\end{bmatrix}.$

Finally, $\boldsymbol{e}(\boldsymbol{z})$ and $\boldsymbol{u}$ are given by

$\boldsymbol{e}(\boldsymbol{z})\vcentcolon=\begin{bmatrix}\frac{1}{\det(F)}F^{\top}F\boldsymbol{E}\\\\[4.0pt]
\frac{\delta\mathcal{H}}{\delta\boldsymbol{B}}(\boldsymbol{z})\\\\[4.0pt]
\boldsymbol{v}\\\\[4.0pt] \frac{\delta\mathcal{H}}{\delta
F}(\boldsymbol{z})\\\\[4.0pt]
1\end{bmatrix}\quad\text{and}\quad\boldsymbol{u}\vcentcolon=\begin{bmatrix}\boldsymbol{F}\\\
r\end{bmatrix},$ (68)

where the $1$ in the last entry of $\boldsymbol{e}(\boldsymbol{z})$ denotes
the constant function from $\mathcal{X}$ to $\mathbb{R}$ whose value is one.

The definitions of $\mathcal{U}(\boldsymbol{z})$ and
$\boldsymbol{e}(\boldsymbol{z})$ immediately yield that the second block row
in (66) coincides with the boundary condition for the total energy boundary
flow specified in (65). Furthermore, since the last entry of
$\boldsymbol{e}(\boldsymbol{z})$ is just $1$, we obtain without further
computation that the first four block rows in

$\mathcal{E}(\boldsymbol{z})=\mathcal{J}(\boldsymbol{z})+\mathcal{B}(\boldsymbol{z})\boldsymbol{u}$
(69)

coincide with the first four equations in (61). To also check that the last
row of (69) coincides with the last equation in (61) we first compute the
variational derivatives of the Hamiltonian which are given by

$\displaystyle\frac{\delta\mathcal{H}}{\delta\boldsymbol{E}}(\boldsymbol{z})$
$\displaystyle=\frac{\epsilon_{0}}{\det(F)}F^{\top}F\boldsymbol{E},\quad\frac{\delta\mathcal{H}}{\delta\boldsymbol{B}}(\boldsymbol{z})=\frac{1}{\mu_{0}\det(F)}F^{\top}F\boldsymbol{B},\quad\frac{\delta\mathcal{H}}{\delta\boldsymbol{v}}(\boldsymbol{z})=\rho\boldsymbol{v},\quad\frac{\delta\mathcal{H}}{\delta\vartheta}(\boldsymbol{z})=\rho
c(\vartheta),$ (70) $\displaystyle\frac{\delta\mathcal{H}}{\delta
F_{ij}}(\boldsymbol{z})$
$\displaystyle=\frac{\epsilon_{0}}{2}\boldsymbol{E}^{\top}\frac{\partial}{\partial
F_{ij}}\left(\frac{F^{\top}F}{\det(F)}\right)\boldsymbol{E}+\frac{1}{2\mu_{0}}\boldsymbol{B}^{\top}\frac{\partial}{\partial
F_{ij}}\left(\frac{F^{\top}F}{\det(F)}\right)\boldsymbol{B}+\rho\frac{\partial\psi_{\mathrm{c}}}{\partial
F_{ij}}(F,\vartheta)$

for all $(i,j)\in\\{1,2,3\\}^{2}$. Consequently, we obtain for the last row of
(69)

$\displaystyle\rho c(\vartheta)\frac{\partial\vartheta}{\partial t}$