# Space-time shape optimization of rotating electric machines

Alessio Cesarano, Charles Dapogny, Peter Gangl

► To cite this version:

Alessio Cesarano, Charles Dapogny, Peter Gangl. Space-time shape optimization of rotating electric machines. Mathematical Models and Methods in Applied Sciences, 2024, 34 (14), pp.2647-2708. 10.1142/S0218202524500568. hal-04452223

HAL Id: hal-04452223

https://hal.science/hal-04452223v1

Submitted on 12 Feb 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
SPACE-TIME SHAPE OPTIMIZATION OF ROTATING ELECTRIC MACHINES

A. CESARANO<fnref n="1" />, C. DAPOGNY², P. GANGL<fnref n="1" />

¹ Johann Radon Institute of Computational and Applied Mathematics, Altenberger Straße 69, 4040 Linz

² Univ. Grenoble Alpes, CNRS, Grenoble INP<fnref n="1" />, LJK, 38000 Grenoble, France.

**ABSTRACT.** This article is devoted to the shape optimization of the internal structure of an electric motor, and more precisely of the arrangement of air and ferromagnetic material inside the rotor part with the aim to increase the torque of the machine. The governing physical problem is the time-dependent, non linear magneto-quasi-static version of Maxwell's equations. This multiphase problem can be reformulated on a 2d section of the real cylindrical 3d configuration; however, due to the rotation of the machine, the geometry of the various material phases at play (the ferromagnetic material, the permanent magnets, air, etc.) undergoes a prescribed motion over the considered time period. This original setting raises a number of issues. From the theoretical viewpoint, we prove the well-posedness of this unusual non linear evolution problem featuring a moving geometry. We then calculate the shape derivative of a performance criterion depending on the shape of the ferromagnetic phase via the corresponding magneto-quasi-static potential. Our numerical framework to address this problem is based on a shape gradient algorithm. The non linear time periodic evolution problems for the magneto-quasi-static potential is solved in the time domain, with a Newton-Raphson method. The discretization features a space-time finite element method, applied on a precise, meshed representation of the space-time region of interest, which encloses a body-fitted representation of the various material phases of the motor at all the considered stages of the time period. After appraising the efficiency of our numerical framework on an academic problem, we present a quite realistic example of optimal design of the ferromagnetic phase of the rotor of an electric machine.

CONTENTS

<table><tbody><tr><td>1. Introduction</td><td>2</td></tr><tr><td>2. Modeling of the physical behavior of an electric motor</td><td>3</td></tr><tr><td>2.1. The eddy current problem for the 3d Maxwell's equations</td><td>3</td></tr><tr><td>2.2. The 2d magneto-quasi-static equation</td><td>5</td></tr><tr><td>2.3. Description of the geometric structure of the motor</td><td>6</td></tr><tr><td>2.4. Notations and mathematical formulation of the magneto-quasi-static problem</td><td>8</td></tr><tr><td>3. Generalities about shape optimization and treatment of an academic problem</td><td>9</td></tr><tr><td>3.1. Shape optimization in a nutshell</td><td>9</td></tr><tr><td>3.2. Calculation of the shape derivative of a "simple" functional of the domain</td><td>11</td></tr><tr><td>4. Calculation of the shape derivative in the magneto-quasi-static context</td><td>15</td></tr><tr><td>4.1. Mathematical preliminaries about the magneto quasi-static evolution problem</td><td>15</td></tr><tr><td>4.2. A few properties of the solution to the magneto-quasi-static problem</td><td>16</td></tr><tr><td>4.3. Calculation of the shape derivative of a physical objective function J(Ω) in the linear magneto-quasi-static setting</td><td>18</td></tr><tr><td>4.4. Extensions of these results to the non linear setting</td><td>25</td></tr><tr><td>5. Numerical methods</td><td>26</td></tr><tr><td>5.1. Sketch of the shape optimization workflow</td><td>26</td></tr><tr><td>5.2. Resolution of the state and adjoint equations</td><td>27</td></tr><tr><td>5.3. Evaluation of the shape derivative J'(Ω)(θ) and calculation of a descent direction</td><td>27</td></tr><tr><td>6. Numerical examples</td><td>28</td></tr><tr><td>6.1. An academic, spatially one-dimensional example</td><td>29</td></tr><tr><td>6.2. Optimization of the structure of a rotating electric machine</td><td>29</td></tr></tbody></table>

<footnote n="1">¹Institute of Engineering Univ. Grenoble Alpes</footnote>
7. **Conclusion and perspectives** 33

Appendix A. **A few useful formulas in shape calculus** 33

Appendix B. **Main results about space-time variational problems** 34

B.1. **A few words about time-space functional spaces** 34

B.2. **Preliminary results about the well-posedness of some non linear problems** 35

B.3. **Proof of Theorem 4.1** 37

References 39

# 1. INTRODUCTION

Electric motors are devices meant to convert electric energy into mechanical energy. For multiple reasons that are related to burning scientific and societal challenges, they have recently aroused a tremendous enthusiasm in the academic and industrial communities. Notably, as they can achieve much superior yield to thermal engines, their systematic deployment would allow for decisive energy savings, not depending on expensive and increasingly scarce fossil fuels. They are also regarded as a promising energy production means in the perspective of the current environmental crisis, as their carbon footprint can be inferior to that of thermal motors by up to 80%. We refer to classical reference books such as [24, 43] for a more exhaustive presentation of electric machines.

Anticipating a little on the more complete description of Section 2, an electric motor is composed of an external, static part called stator, and an internal, rotating part – the rotor. An alternating electric current is injected through the coils placed in the stator, thus generating a magnetic field which sets the rotor in motion. The resulting mechanical work is transmitted to a shaft attached to its core and can be used directly or stored, see Fig. 1 below for an illustration. The efficiency of this workflow strongly depends on the arrangement of the costly components of the rotor, in particular the ferromagnetic material and the permanent magnets (usually made of rare-earth elements), and the optimization of their layout presages tremendous performance improvements, see [55] about this issue.

Over the past decades, optimal design has been a thriving field of research at the crossroads of mathematics, physics and engineering; multiple frameworks are available to address such issues, with competing assets and drawbacks. Among them, shape optimization aims to improve an initial guess through smooth variations of the boundaries or material interfaces at play; this relies on the information contained in the “shape gradients” of the objective and constraint functionals of the problem – see e.g. [2, 3, 23, 36, 54, 67] and [20, 31, 11, 58] for numerical implementations. On the other hand, topology optimization methods affect the connectivity of a design, for instance by leveraging topological derivatives, indicating where holes can be fruitfully nucleated inside a shape [5, 6, 56, 57], or via a suitable reformulation of the physical situation and the optimal design problem under scrutiny in terms of density functions [10, 66]. Surprisingly enough, such techniques have only been applied recently in the physical context of electric motors. Among the contributions in this direction and without aiming exhaustivity, let us mention the articles [34, 35, 52] using density-based topology optimization strategies, and those [15, 49, 61] relying on the level set method; see also [44] about a coupling between the level set method and a body-fitted meshed representation of the machine enabling more accurate finite element simulations. We refer to [26, 51] for recent surveys of optimal design techniques in the field of electromagnetism, and to [12, 73] for a particular focus on their application to electric machines.

The specific setting of electric motors raises several challenges. At first, although the boundary value problem governing the behavior of the motor is a reduction of the full 3d Maxwell’s equations to its 2d cross-section, it features a non linear material law which is essential to capture the decisive saturation effect of the ferromagnetic constituent of the rotor [19]. The physical problem at hand is also time-dependent; moreover, as the rotor is in rotation, the optimal design criterion depends on all the stages of this periodic motion; this feature is fairly awkward in view of optimal design because of the tremendous computational cost incurred, see e.g. [40, 65]. To the best of our knowledge and with the exception of the recent contribution [29], all the existing studies about the shape or topology optimization of electric motors resort to an additional simplification referred to as the magneto-static approximation of Maxwell’s equations: the physical problem
is thereby reduced to a series of decoupled static (elliptic and non linear) problems posed on the various
rotated configurations of the motor, which can be solved independently from one another. Unfortunately,
this convenient reformulation does not allow to account for subtle and realistic temporal effects such as eddy
currents.

The present article is a natural continuation to the series [30, 53]. We wish to optimize the repartition of the constituent materials of an electric motor; yet, contrary to the prevailing practice in the literature, we rely on the genuine magneto-quasi-static description of the physical situation. The boundary value problem at stake can no longer be reduced to a collection of independent stationary problems: it is a non linear, time-dependent problem of mixed elliptic-parabolic type [8], equipped with time periodicity conditions, in which the optimization criterion brings into play all the rotated geometric configurations of the internal structure of the motor.

From the theoretical viewpoint, our first contribution is to prove the well-posedness of this physical evolution problem, which leverages techniques from the theory of non linear partial differential equations. We next analyze shape optimization problems where the performance functionals at play depend on the solution to this problem. From the numerical viewpoint, the specific nature of the considered time-dependent problem raises the need for an adapted strategy. Our framework relies on a space-time finite element discretization, see e.g. [46, 68] in the context of linear problems, [71] for a quasilinear problem, or more recently [28] for an application in the context of electric machines. In a nutshell, the time variable is treated as if it were an additional space variable, thus converting a time-dependent problem posed in a *d*-dimensional domain into a static problem on a (*d* + 1)-dimensional space time cylinder; the latter is solved by applying the finite element method on a mesh of the latter. While this approach suffers from an increase in the dimension of the problem, it enjoys a number of unique assets over classical time-stepping methods (see e.g. [69]):

* The adaptive refinement of the mesh and the parallelization of the solution of the evolution problem can be realized in space and time, jointly;

* The adjoint problem involved in the shape derivative of the optimization criterion can be solved concurrently with the state equation, thus allowing for a further level of parallelization;

* Shapes evolving over the considered time period are static in the space-time domain; hence, evolution problems involving such moving domains can be treated by “standard”, time independent numerical methods once a space-time mesh adapted to the moving geometry is available;

* Time periodicity conditions are straightforward to enforce by direct identification of the degrees of freedom on the bottom and top sections of the space-time cylinder.

The remainder of this article is organized as follows. In Section 2, we present the physical and mathemati-
cal settings of the boundary value problem governing the physics of motors, as well as the shape optimization
problem considered in this setting. The next Section 3 is devoted to an academic, albeit instructive pre-
liminary situation, that of the calculation of the shape derivative of a “simple” functional, depending on
a domain Ω via the integral of a given function over a collection of deformed versions of Ω. In Section 4,
we outline a few mathematical features of the considered magneto-quasi-static evolution problems, and we
detail the calculation of a shape derivative in this context. After discussing a few details of our numerical
implementation in Section 5, we show in Section 6 two numerical examples illustrating the previous devel-
opments; a few conclusions and perspectives of our work are then given in Section 7. The article ends with
two Appendices A and B, where a few “classical” technical results are recalled, and the main steps of the
proof of the well-posedness of the considered magneto-quasi-static problem are sketched.

## 2. MODELING OF THE PHYSICAL BEHAVIOR OF AN ELECTRIC MOTOR

This section introduces the physical and mathematical aspects of the optimal design of electric motors. After a few generalities about these devices and their modeling in terms of the Maxwell's equations in Section 2.1, we discuss the approximate 2d magneto-quasi-static setting in Section 2.2. The internal structure of the considered motors is presented in Section 2.3, and we eventually set the mathematical framework and the main notations used throughout the article in the next Section 2.4.
## 2.1. The eddy current problem for the 3d Maxwell's equations

This section describes, in an informal manner, the general operation of an electric motor and its mathematical formulation via the eddy current problem stemming from the 3d Maxwell's equations; our main references about this topic are Chap. 2 in [26], [62] and [72].

The motor under scrutiny has a cylindrical structure and its transverse section is sketched in Fig. 1. It is made of an inner rotating part called rotor, and a static outer part, the stator; these regions are separated by a thin air gap and the core of the rotor is connected to a transmission shaft. The stator is an arrangement of ferromagnetic material and coils made of thin copper wires; the rotor is composed of ferromagnetic material, air inclusions and several permanent magnets.

In the classical language of electromagnetism, for which we refer to [33, 39] and whose main concepts and notations are used in this presentation, the device is activated when a time-dependent electric current $J_i$ is powered into the coils of the stator. According to Ampère's law, this generates a time-dependent magnetic field $H$ in the whole region of interest, which in turn induces an electric field $E$ by Faraday's law. The resulting Lorentz force sets in motion the charge carriers contained in the rotor and thereby the attached shaft; this mechanical work is eventually used or stored.

In order to put this rough sketch into equations, let us first note that, as opposed to light waves, the present application falls in the regime of low-frequency electromagnetism, so that the displacement currents, i.e. the time derivative of the magnetic induction $D$, can be neglected, see [62] for an intuitive explanation of this fact and [17] for a mathematical justification. The magnetic field $H$ generated around the motor is then related to the total density of current $J$ via the following magneto-quasi-static version of Ampère's law:

$$ (2.1) \qquad \mathrm{curl}(\mathbf{H}) = \mathbf{J}. $$

On the other hand, this field is related to the magnetic induction (or magnetic flux density) $B$ via the following constitutive law:

$$ (2.2) \qquad \mathbf{H} = \nu(|\mathbf{B}|)\mathbf{B} - \mathbf{M}. $$

Here, the coefficient $\nu$ represents the magnetic reluctivity, i.e. the inverse of the magnetic permeability; it expresses how the material develops a magnetic field in the presence of a magnetization force. The vector field $M$ denotes the permanent magnetization, i.e. the density of permanent magnetic dipoles within the motor; this field vanishes outside the permanent magnets of the rotor. The reluctivity $\nu$ takes different values depending on the material phase – a dependence with respect to the spatial location which is omitted in (2.2) for simplicity. It assumes constant values $\nu_a, \nu_m, \nu_c$ in air, permanent magnets and copper, respectively. By contrast, in ferromagnetic materials, it is a function of the intensity $|\mathbf{B}|$ of the magnetic induction: for small values of $H$, the material amplifies the magnetic flux density $B$, whereas for high values of $H$, the relation between both quantities is similar to that occurring in void. We refer to [59] about the calibration of such a reluctivity function from physical measurements.

According to Gauss's law, the magnetic induction $B$ is solenoidal, i.e.

$$ (2.3) \qquad \mathrm{div}(\mathbf{B}) = 0. $$

In turn, the time variations of $B$ induce an electric field $E$ according to Faraday's law:

$$ (2.4) \qquad \mathrm{curl}(\mathbf{E}) = -\frac{\partial \mathbf{B}}{\partial t}. $$

This model is completed by Ohm's law, which expresses the density $J$ of current inside the region of interest as the superposition of the density $J_i$ of current impressed in the coils and that produced by the Lorentz force acting on the charge carriers:

$$ (2.5) \qquad \mathbf{J} = \mathbf{J}_i + \sigma(\mathbf{v} \times \mathbf{B} + \mathbf{E}), $$

where $\mathbf{v}$ is the velocity field of the rotating part and $\sigma$ is the electric conductivity, which takes a different, constant value inside each material phase.

The combination of these relations allows to express all the above quantities in terms of a single vector field – the vector potential $\mathbf{A}$ – and to characterize the latter by a partial differential equation. Indeed, at first, (2.3) implies the existence of a vector field $\mathbf{A}$ such that $\mathbf{B} = \mathrm{curl}(\mathbf{A})$, which is determined up to the addition of the gradient of a scalar field. Judging from (2.4), the vector potential $\mathbf{A}$ can be selected as the
unique such vector field satisfying $\mathbf{E} = -\frac{\partial \mathbf{A}}{\partial t}$. Now injecting the constitutive relation (2.2) and the expression (2.5) of the current density into Ampère's law (2.1), we obtain the three-dimensional magneto-quasi-static version of Maxwell's equation, also referred to as the eddy current problem:

$$ (2.6) \qquad \sigma \left( \frac{\partial \mathbf{A}}{\partial t} - \mathbf{v} \times \mathbf{curl}(\mathbf{A}) \right) + \mathbf{curl}(\nu(|\mathbf{curl}(\mathbf{A})|)\mathbf{curl}(\mathbf{A})) = \mathbf{J}_i + \mathbf{curl}(\mathbf{M}). $$

The partial differential equation (2.6) is complemented by suitable time and spatial boundary conditions. As far as the former are considered, our study concerns the permanent regime of the motor, where the material distribution retrieves its initial configuration and all physical quantities return to their initial values after one rotation. The equation (2.6) is then equipped with a time periodic condition over the considered time period $[0, T]$. As regards boundary conditions, at each time $t \in [0, T]$, the usual transmission conditions are assumed at the interfaces between the various material phases, and no magnetic flux leaves the computational domain:

$$ (2.7) \qquad \mathbf{B} \cdot \mathbf{n} = 0 \text{ on the outer boundary of the stator,} $$

where $\mathbf{n}$ is the unit normal vector field to this boundary, pointing outward the device.

## 2.2. The 2d magneto-quasi-static equation

As we have mentioned, the full 3d motor under scrutiny has a cylindrical structure. Taking advantage of its relatively large size along the shaft axis, that we suppose without loss of generality to be the line passing through the origin, oriented along the third coordinate vector $\mathbf{e}_3$, and of the invariance of its geometry in this direction, it is customary to approximate the three-dimensional problem (2.6) with a simpler two-dimensional version, posed in the cross-section of the device [42, 70].

This reduction is justified since the impressed current density $\mathbf{J}_i$ and the magnetization $\mathbf{M}$ are of the form

$$ \mathbf{J}_i(t, x_1, x_2, x_3) = (0, 0, f(t, x_1, x_2)), \quad \mathbf{M}(t, x_1, x_2, x_3) = (M_1(t, x_1, x_2), M_2(t, x_1, x_2), 0) $$

and the velocity $\mathbf{v}$ representing the rotation of the rotor reads

$$ \mathbf{v}(t, x_1, x_2, x_3) = (v_1(t, x_1, x_2), v_2(t, x_1, x_2), 0). $$

In turn, the vector potential $\mathbf{A}$ is oriented along the axis $\mathbf{e}_3$ and its values depend only on the spatial position in the cross-section $D$, i.e.

$$ \mathbf{A}(t, x_1, x_2, x_3) = (0, 0, u(t, x_1, x_2)), $$

for a suitable function $u(t, x_1, x_2)$.

Under these assumptions, the eddy current equation (2.6) reduces to a 2d partial differential equation for the scalar field $u$ which is posed in the cross-section $D$ of the machine; an elementary calculation indeed yields:

$$ (2.8) \qquad \sigma \left( \frac{\partial u}{\partial t} + \mathbf{v} \cdot \nabla u \right) - \operatorname{div}(\nu(|\nabla u|)\nabla u) = f - \operatorname{div} \mathbf{M}^{\perp}, \quad x \in D, t \in (0, T). $$

Here, we have dropped the bold font for vector fields when they lie inside the cross-section plane of the motor; we have set

$$ M(t, x_1, x_2) := (M_1(t, x_1, x_2), M_2(t, x_1, x_2)), \quad v(t, x_1, x_2) := (v_1(t, x_1, x_2), v_2(t, x_1, x_2)), $$

and we use the notation $y^{\perp} = (-y_2, y_1)$ for the 90° clockwise rotate of a two-dimensional vector $y = (y_1, y_2)$.

Eventually, the 2d equation (2.8) is complemented with the homogeneous Dirichlet boundary condition

$$ u = 0 \text{ on } \partial D, $$

inherited from the no-flux assumption (2.7) made in the complete 3d model of Section 2.1. As far as temporal boundary conditions are concerned, the time periodic setting implies that:

$$ \forall x = (x_1, x_2) \in D, \quad u(0, x) = u(T, x). $$

**Remark 2.1.** The mathematical analyses conducted in this article can be adapted straightforwardly to handle other time and boundary conditions about $u$, and notably the more classical initial condition $u(0, \cdot) = u_0$. Since the time periodic setting is the most relevant one for our purpose, we do not discuss those alternative instances for brevity.
<figure><img src="image_2.png" /><figcaption>FIGURE 1. Two-dimensional cross-section D of an electric machine of the form considered in Section 2: the ferromagnetic material is depicted in red, air is in dark blue, the coils are in yellow, and the permanent magnets are in light blue.</figcaption></figure>

In engineering applications, the performance of the motor is often measured in terms of the average torque induced by the rotation of the motor shaft over the time period (0, T). This quantity Tor(u) is evaluated in terms of the potential u by the so-called Arkkio's method [63]:

$$ (2.9) \qquad \mathrm{Tor}(u) := \frac{1}{T} \frac{L \nu_a}{(r_s - r_r)} \int_0^T \int_{\Sigma} Q(x_1, x_2) \nabla u \cdot \nabla u \, dx dt $$

where L is the length of the true, 3d electric machine in the e3 direction, Σ is an annulus with inner and outer radii $r_r < r_s$ lying in the air gap between the rotor and the stator, and the 2 × 2 matrix $Q(x)$ is defined by:

$$ (2.10) \qquad \forall x = (x_1, x_2) \in D, \quad Q(x_1, x_2) = \frac{1}{\sqrt{x_1^2 + x_2^2}} \begin{pmatrix} x_1 x_2 & \frac{x_2^2 - x_1^2}{2} \\ \frac{x_2^2 - x_1^2}{2} & -x_1 x_2 \end{pmatrix}. $$

## 2.3. Description of the geometric structure of the motor

In this section, we get into more specifics about the constituent material phases of the motor and we present the corresponding model governing the expressions of the conductivity and reluctivity coefficients $\sigma$ and $\nu$ in the formulation of the evolution problem (2.8). In the first Section 2.3.1, we describe the internal structure of the considered motor when it is at rest; we then discuss its time evolution due to the rotation of the rotor in Section 2.3.2.

### 2.3.1. Structure of the motor at rest

The structure of the 2d cross-section D reflects the three-dimensional arrangement of the motor introduced in the previous Section 2.1; it features three disjoint annulus-shaped regions:

$$ \overline{D} = \overline{D}_{\mathrm{rot}} \cup \overline{D}_{\mathrm{gap}} \cup \overline{D}_{\mathrm{stat}}, $$

where the open sets $D_{\mathrm{rot}}$, $D_{\mathrm{gap}}$ and $D_{\mathrm{stat}}$ stand for the rotor, the separating air gap and the stator, respectively. The stator $D_{\mathrm{stat}}$ consists of

* A region $D_{\mathrm{stat,f}}$ filled with ferromagnetic material;

* A region $D_{\mathrm{stat,a}}$ made of air;

* The coils $D_{\mathrm{stat,c}}$ featuring copper wires.

The part $D_{\mathrm{rot}}$ is composed of:
* A ferromagnetic core, occupying the region $\Omega$;

* Air, occupying the region $\Omega_a$;

* Permanent magnets represented by the domain $D_{\text{mag}}$.

In this study, the interface $\Gamma := \partial\Omega \cap \partial\Omega_a$ between the phases of the rotor made of ferromagnetic material and air is the center of attention when it comes to optimal design, while the region $D_{\text{mag}}$ occupied by the permanent magnets and the stator $D_{\text{stat}}$ are not subject to optimization, see Fig. 1 for an illustration of this structure.

The conductivity and reluctivity coefficients $\sigma$ and $\nu$ have different expressions in the various phases featured in these decompositions. In order to emphasize their dependence on the actual shape $\Omega \subset D_{\text{rot}}$ of the optimized ferromagnetic core, they are labeled with an $\Omega$ subscript.

$$ (2.11) \quad \sigma_{\Omega}(x) = \begin{cases} \sigma_f & x \in \Omega \cup D_{\text{stat,f}}, \\ \sigma_a & x \in \Omega_a \cup D_{\text{gap}} \cup D_{\text{stat,a}}, \\ \sigma_m & x \in D_{\text{mag}}, \\ \sigma_c & x \in D_{\text{stat,c}}, \end{cases} \quad \text{and} \quad \nu_{\Omega}(x, s) = \begin{cases} \hat{\nu}(s) & x \in \Omega \cup D_{\text{stat,f}}, \\ \nu_a & x \in \Omega_a \cup D_{\text{gap}} \cup D_{\text{stat,a}}, \\ \nu_m & x \in D_{\text{mag}}, \\ \nu_c & x \in D_{\text{stat,c}}. \end{cases} $$

In practice, the value of the electric conductivity $\sigma_a$ in the air region $\Omega_a \cup D_{\text{gap}} \cup D_{\text{stat,a}}$ is negligible; moreover, due to the laminated structure of the part filled with ferromagnetic material, which is made of iron sheets, and since the copper wires typically used in electric machines are insulated, the values of $\sigma$ are also negligible in $\Omega \cup D_{\text{stat,f}}$ and $D_{\text{stat,c}}$:

$$ \sigma_f = \sigma_a = \sigma_c = 0. $$

The reluctivity $\nu$ takes constant, positive values $\nu_a$, $\nu_m$ and $\nu_c$ in air, in the magnets, and in the wires, respectively. However, it is a function $\hat{\nu} : \mathbb{R}_+ \to \mathbb{R}_+$ of the amplitude of the gradient of the potential $u$ in the regions $\Omega$ and $D_{\text{stat,f}}$ occupied by ferromagnetic material. As evidenced in [59], the physical properties of the latter imply that the mapping $s \mapsto \hat{\nu}(s)s$ in (2.11) is strongly monotone and Lipschitz continuous, i.e.

$$ (2.12) \quad \text{There exists } 0 < \nu \le \bar{\nu} < \infty \text{ s.t. } \forall s_1, s_2 \in \mathbb{R}_+, \quad \begin{cases} (\hat{\nu}(s_1)s_1 - \hat{\nu}(s_2)s_2)(s_1 - s_2) \ge \nu(s_1 - s_2)^2, \\ |\hat{\nu}(s_1)s_1 - \hat{\nu}(s_2)s_2| \le \bar{\nu}|s_1 - s_2|. \end{cases} $$

In particular, (2.12) implies that $\hat{\nu}$ is uniformly bounded away from 0 and $\infty$:

$$ \forall s \in \mathbb{R}_+, \quad \nu \le \hat{\nu}(s) \le \bar{\nu}. $$

### 2.3.2. Description of the motion of the rotor

The operation of the motor is analyzed through a representative time period $(0, T)$; during the latter, $D_{\text{rot}}$ moves according to the smooth velocity field $v : (0, T) \times \overline{D_{\text{rot}}} \to \mathbb{R}^d$. The corresponding flow $\varphi : (0, T) \times D_{\text{rot}} \to \mathbb{R}^d$ is the solution to the following ordinary differential equation:

$$ (2.13) \quad \forall x \in D_{\text{rot}}, \quad \begin{cases} \frac{\partial \varphi}{\partial t}(t, x) = v(t, \varphi(t, x)) & \text{on } (0, T), \\ \varphi(0, x) = x; \end{cases} $$

we also let the shorthand $\varphi_t \equiv \varphi(t, \cdot)$. It follows from the standard theory of ordinary differential equations that:

$$ (2.14) \quad \bullet \quad \varphi \text{ is smooth over } [0, T] \times \overline{D_{\text{rot}}}; $$

* For each $t \in [0, T]$, the mapping $x \mapsto \varphi(t, x)$ is a diffeomorphism from $\overline{D_{\text{rot}}}$ onto itself;

We also assume that the considered motion of the rotor is T-periodic, and so it holds:

$$ (2.15) \quad \forall x \in D_{\text{rot}}, \quad \varphi(0, x) = \varphi(T, x) = x. $$

Eventually, for convenience in the mathematical analysis, we extend $\varphi_t$ to a smooth diffeomorphism of $D$ such that:

$$ (2.16) \quad \varphi_t(D_{\text{gap}}) = D_{\text{gap}}, \quad \text{and} \quad \varphi_t \equiv \text{Id} \text{ on } D_{\text{stat}}. $$

This motion of the rotor induces an evolution in time of the constituent material phases of $D_{\text{rot}}$. Let $\Omega(t) := \varphi_t(\Omega)$ and $\Gamma(t) := \varphi_t(\Gamma)$ (resp. $D_{\text{mag}}(t) = \varphi_t(D_{\text{mag}})$, etc.) denote the deformed versions of $\Omega$ and
$\Gamma$ (resp. $D_{\text{mag}}$, etc.) induced by $\varphi_t$. From the mathematical viewpoint, the material coefficients $\sigma$ and $\nu$ in (2.21) depend on time, as:

$$ (2.17) \qquad \forall x \in D, s \in \mathbb{R}_+, \quad \sigma_{\Omega(t)}(\varphi_t(x)) = \sigma_{\Omega}(x) \text{ and } \nu_{\Omega(t)}(\varphi_t(x), s) = \nu_{\Omega}(x, s). $$

In practice, the cross-section $D$ of the motor is a two-dimensional disk with center 0 and $\varphi_t$ corresponds to the rotation with origin 0 and angle $\alpha(t) := \frac{2\pi t}{T}$:

$$ (2.18) \qquad \varphi_t(x) = R_{\alpha(t)} \begin{pmatrix} x_1 \\ x_2 \end{pmatrix}, \quad \text{where } R_{\alpha(t)} := \begin{pmatrix} \cos \alpha(t) & -\sin \alpha(t) \\ \sin \alpha(t) & \cos \alpha(t) \end{pmatrix}. $$

With this definition, a simple calculation reveals that the velocity field $v$ is divergence-free and that the Jacobian determinant of $\varphi_t$ identically equals 1:

$$ (2.19) \qquad \forall t \in [0, T], x \in D, \text{ div } v(t, x) = 0, \text{ and } \det \nabla \varphi_t(x) = 1. $$

**Remark 2.2.** Even though this very specific setting where the space dimension $d$ equals 2 and the motion $\varphi(t, x)$ of the rotor accounts for a rotation is the relevant one in practice, the general case where $d$ is arbitrary and $\varphi_t$ is an arbitrary smooth diffeomorphism of $D$ satisfying (2.14) to (2.16), with associated velocity field $v(t, x) = \frac{\partial \varphi}{\partial t}(t, \varphi_t^{-1}(x))$, does not pose any additional difficulty in the analysis. For this reason, in the subsequent developments, we retain generic notations for $d$ and $\varphi$ insofar as possible.

## 2.4. Notations and mathematical formulation of the magneto-quasi-static problem

The following notations are adopted throughout the remainder of this article.

* We denote by $\text{Id} : \mathbb{R}^d \to \mathbb{R}^d$ the identity mapping of $\mathbb{R}^d$ and by $\text{I} \in \mathbb{R}^{d \times d}$ the identity matrix with size $d$.

* For any time $t \in [0, T]$, $n_{\Omega(t)}$ (resp. $n$) stands for the unit normal vector to the deformed interface $\Gamma(t)$ (resp. to the interface at rest $\Gamma$), pointing outward $\Omega(t)$ (resp. outward $\Omega$).

* We denote by

$$ Q = (0, T) \times D, \quad \text{and} \quad Q_{\Omega} := \{ (t, x) \in Q, x \in \Omega(t) \} $$

the (open) space-time cylinders induced by $D$ and $\Omega$ over the time period $(0, T)$.

* Let $\alpha : D \to \mathbb{R}$ be a smooth quantity on $\bar{\Omega}$ and $D \setminus \Omega$, which is possibly discontinuous across the interface $\Gamma$; we denote by

$$ \alpha^{\pm}(x) = \lim_{\substack{t \to 0 \\ t \to 0}} \alpha(x \pm tn(x)) $$

the one-sided limits of $\alpha$ at $x \in \Gamma$, from outside and inside $\Omega$, respectively.

* In particular, let $u : D \to \mathbb{R}$ be a smooth function on $\bar{\Omega}$ and $D \setminus \Omega$, which is possibly discontinuous at the interface $\Gamma$. For any $x \in \Gamma$, we denote by

$$ \frac{\partial u^{\pm}}{\partial n}(x) = \lim_{\substack{t \to 0 \\ t \to 0}} \nabla u(x \pm tn(x)) \cdot n(x) \text{ and } \nabla_{\Gamma} u^{\pm}(x) = \nabla u(x \pm tn(x)) - \frac{\partial u^{\pm}}{\partial n}(x)n(x), $$

the one-sided normal derivatives and tangential gradients of $u$ at $x$.

* Let $w : D \to \mathbb{R}^d$ be a smooth enough vector field. We denote by

$$ w_{\Gamma} : \Gamma \to \mathbb{R}^d, \quad w_{\Gamma}(x) = w(x) - (w(x) \cdot n(x))n(x) $$

the tangential part of $w$ on $\Gamma$.

* Let $u : Q \to \mathbb{R}$ be a smooth enough scalar function of time and space; then,

- $\nabla u(t, x)$ denotes the gradient of $u$ with respect to the space variable $x$ only;

- $\frac{\partial u}{\partial t}(t, x)$ is the partial derivative of $u$ with respect to the time variable, i.e. the derivative of the partial mapping $t \mapsto u(t, x)$ for fixed $x \in D$;
− du/dt (t, x) is the total (or material) derivative of u according to the velocity field v(t, x), that is:

$$
(2.20) \qquad \frac{du}{dt}(t, x) = \frac{\partial u}{\partial t}(t, x) + v(t, x) \cdot \nabla u(t, x);
$$

in other terms, $\frac{du}{dt}(t, x)$ is the derivative of the mapping $t \mapsto u(t, \varphi_t(y))$ evaluated at time $t$ and
point $y = \varphi_t^{-1}(x)$.

* Let $w = (w_1, \dots, w_d) : Q \to \mathbb{R}^d$ be a (smooth enough) vector-valued function of space and time. We denote by $\nabla w$, or often $[\nabla w]$, the derivative of $w$ with respect to the spatial variable, i.e. the $d \times d$ matrix field with entries:

$$
[\nabla w]_{ij} = \frac{\partial w_i}{\partial x_j}, \quad i, j = 1, \dots, d.
$$

Summarizing the previous discussions and taking these notations into account, the transverse component
of the vector potential inside the motor, which we denote by $u_{\Omega}: Q \to \mathbb{R}$ to emphasize its dependence on
the phase $\Omega$ to be optimized, satisfies the following evolution problem:

$$
(2.21) \qquad \left\{
\begin{array}{ll}
\sigma_{\Omega(t)} \frac{du_{\Omega}}{dt} - \text{div} \left( \nu_{\Omega(t)}(x, |\nabla u_{\Omega}|) \nabla u_{\Omega} \right) = f & \text{in } (0, T) \times D, \\
u_{\Omega}(t, x) = 0 & \text{for } t \in (0, T), x \in \partial D, \\
u_{\Omega}(0, x) = u_{\Omega}(T, x) & \text{for } x \in D_{\text{mag}},
\end{array}
\right.
$$

where the coefficients $\sigma_{\Omega(t)}$ and $\nu_{\Omega(t)}$ are given by (2.11). The precise mathematical setting of this problem and its well-posedness are discussed in Section 4.1 below.

**Remark 2.3.**

* In passing from the formulation (2.8) to (2.21), we have absorbed the term $\operatorname{div}(M^{\perp})$ at the right-hand side of the problem (2.8) into the source $f$. This convenient notational simplification is slightly abusive though, as both terms $f$ and $\operatorname{div}(M^{\perp})$ are not exactly of the same mathematical nature: while $f(t, \cdot)$ will typically belong to $L^2(D)$ in our mathematical analysis, the quantity $M^{\perp}(t, \cdot)$ will be smooth on $D_{\text{mag}}$ and vanish outside this set, so that $\operatorname{div}(M^{\perp})(t, \cdot) \in H^{-1}(D)$. Nevertheless, the subsequent developments are easily adapted to handle this term, see Remark 4.4 about this point.

* The subsequent mathematical discussions rely on the fact that the domain $D_{\text{mag}}$ is non empty. If it were to be the case, (2.21) would simply boil down to decoupled elliptic problems.

* As we have mentioned, the non linearity of the boundary value problem (2.21) (notably, the dependence of the reluctivity on |∇u|) is key in the correct modeling of the physical problem. However, the linear setting, where

$$
\hat{v}(s) \equiv \nu_f \text{ for some } \nu_f > 0,
$$

will sometimes be considered in the mathematical setting, as it allows for a less technical exposition
of the main ideas.

3. GENERALITIES ABOUT SHAPE OPTIMIZATION AND TREATMENT OF AN ACADEMIC PROBLEM

In this section, we enter into the optimal design framework of this article. After introducing the model shape
optimization problem under scrutiny in Section 3.1, we detail in Section 3.2 the calculation of a “simple”
shape derivative as a handful preliminary to our subsequent analyses.

3.1. Shape optimization in a nutshell

In this article, we aim to optimize the shape of the phase $\Omega$ made of ferromagnetic material with respect to
a physical criterion, that is

$$
(3.1) \qquad \min_{\Omega \subset D_{\text{rot}}} J(\Omega).
$$

The performance criterion under scrutiny $J(\Omega)$ depends on $\Omega$ via the transverse component of vector potential $u_{\Omega}(t, x)$ of the magnetic flux density, solution to the evolution problem (2.21). We consider a model objective function of the form

$$
(3.2) \qquad J(\Omega) = \int_Q j(u_{\Omega}(t, x)) \, dx dt,
$$
where $j: \mathbb{R} \to \mathbb{R}$ is a smooth function, satisfying the growth conditions:

$$ \exists C > 0, \forall u \in \mathbb{R}, \quad |j(u)| \le C(1+|u|)^2, \quad |j'(u)| \le C(1+|u|), \quad \text{and} \quad |j''(u)| \le C. $$

Note that different quantities of interest could be considered without much changes to our analysis. Notably, $J(\Omega)$ could also depend on the gradient of $u_\Omega$, as in the torque functional $\text{Tor}(u_\Omega)$ defined in (2.9), see Section 6. Besides, constraints (e.g. on the volume of $\Omega$) could be added to the problem (3.1) without much change to the subsequent theoretical developments.

The numerical solution of the problem (3.1) usually calls for the derivative of the objective function $J(\Omega)$ with respect to the domain – a notion which can be defined in various ways. Here, we rely on Hadamard’s boundary variation method, see e.g. [36, 54, 67] or [2] for a recent overview. In a nutshell, variations of a shape $\Omega \subset D_{\text{rot}}$ are considered under the form

$$ \Omega_{\theta} := (\text{Id} + \theta)(\Omega), \quad \theta \in W^{1,\infty}(\mathbb{R}^d, \mathbb{R}^d), \quad \|\theta\|_{W^{1,\infty}(\mathbb{R}^d, \mathbb{R}^d)} < 1, $$

that is, $\Omega_\theta$ is a version of $\Omega$ whose points have been displaced according to the “small” vector field $\theta$, see Fig. 2.

<figure><img src="image_3.png" /><figcaption>FIGURE 2. The variation $\Omega_\theta$ of $\Omega$ featured in the method of Hadamard is obtained by displacing its points according to the “small” vector field $\theta$.</figcaption></figure>

A function $F(\Omega)$ of the domain is then called shape differentiable at a particular domain $\Omega$ provided the underlying mapping $\theta \mapsto F(\Omega_\theta)$, from $W^{1,\infty}(\mathbb{R}^d, \mathbb{R}^d)$ into $\mathbb{R}$, is Fréchet differentiable at $\theta = 0$. The corresponding shape derivative $F'(\Omega)(\theta)$ gives rise to the following expansion:

$$ (3.3) \qquad F(\Omega_\theta) = F(\Omega) + F'(\Omega)(\theta) + o(\theta), \quad \text{where} \quad \frac{|o(\theta)|}{||\theta||_{W^{1,\infty}(\mathbb{R}^d, \mathbb{R}^d)}} \xrightarrow{\theta \to 0} 0. $$

In particular, a deformation $\theta \in W^{1,\infty}(\mathbb{R}^d, \mathbb{R}^d)$ such that $F'(\Omega)(\theta) < 0$ is a descent direction for $F$ from $\Omega$, since for a small enough descent step $\tau > 0$,

$$ F(\Omega_{\tau\theta}) = F(\Omega) + \tau F'(\Omega)(\theta) + o(\tau), $$

i.e. a “small” deformation of $\Omega$ according to $\theta$ produces a design $\Omega_{\tau\theta}$ which is “better” with respect to the criterion $F(\Omega)$.

In the present application, we only optimize the interface $\Gamma$ between the phases $\Omega$, $\Omega_a$ of the rotor, respectively made of ferromagnetic material and air. Accordingly, the deformations $\theta$ involved in Hadamard’s method ought to be restricted to the subspace

$$ \Theta_{\text{ad}} := \{ \theta \in W^{1,\infty}(\mathbb{R}^d, \mathbb{R}^d), \theta = 0 \text{ on } \overline{D}_{\text{stat}}, \theta \cdot n = 0 \text{ on } \partial D_{\text{rot}} \cup \partial D_{\text{mag}} \}. $$

Let us end this section with a few informal words about the general structure of the shape derivative of a “smooth enough” function of the domain $F(\Omega)$, such as those considered in the sequel; we refer to e.g. §5.9 in [36] for more details about this subject. Usually, several equivalent formulas are available for $F'(\Omega)(\theta)$. A
rigorous calculation usually yields the so-called volume (or distributed) expression of this derivative, which is of the form:

$$ (3.4) \qquad F'(\Omega)(\theta) = \int_D (M_\Omega : \nabla\theta + s_\Omega \cdot \theta) \, dx, $$

where $M_\Omega: D \to \mathbb{R}^{d \times d}$ and $s_\Omega: D \to \mathbb{R}^d$ are matrix and vector fields depending on $\Omega$, $u_\Omega$ and the function $F(\Omega)$, notably via a so-called adjoint state $p_\Omega$, solution to a problem similar to (2.21), with a different right-hand side. Often, (3.4) can be turned into an equivalent, surface expression:

$$ (3.5) \qquad F'(\Omega)(\theta) = \int_{\Gamma} v_{\Omega}(\theta \cdot n) \, ds, $$

featuring a scalar field $v_\Omega: \partial\Omega \to \mathbb{R}$. This alternative structure encodes the intuitive fact that $F'(\Omega)(\theta)$ solely depends on the normal component of the deformation $\theta$ – roughly speaking, vector fields taking only tangential directions on the boundary $\partial\Omega$ only account for a “reparametrization” of $\Omega$. The surface form (3.5) is also helpful for numerical purpose, as it immediately reveals a descent direction for $F(\Omega)$. Indeed, letting $\theta = -v_\Omega n$ on $\partial\Omega$ immediately yields $F'(\Omega)(\theta) < 0$.

Mathematically, the surface formula (3.5) for $F'(\Omega)(\theta)$ is often obtained from its volume counterpart (3.4) after integration by parts, assuming some regularity from the boundary $\partial\Omega$ and the functions $u_\Omega, p_\Omega$, see e.g. the proofs of Lemma 3.1 and Proposition 4.1 below.

## 3.2. Calculation of the shape derivative of a “simple” functional of the domain

We now present the calculation of the derivative of a “simple” functional $J(\Omega)$, depending on the shape $\Omega$ via the integration of a fixed, smooth function $f(t, x)$ on the transformed versions $\Omega(t) = \varphi_t(\Omega)$ of $\Omega$ for $t \in (0, T)$. We seize this opportunity to introduce a few technical preliminaries useful throughout this article in Section 3.2.1, before proceeding to the calculation, properly speaking, in Section 3.2.2.

### 3.2.1. Technical preliminaries

For any vector field $\theta \in \Theta_{ad}$, let us introduce the corresponding space-time deformation mapping $\Theta: Q \to (0, T) \times \mathbb{R}^d$

$$ (3.6) \qquad \Theta(t, x) := (\Theta_t(t, x), \Theta_x(t, x)) = (t, \varphi_t \circ (\mathrm{Id} + \theta) \circ \varphi_t^{-1}(x)), $$

whose notation omits its dependence on $\theta$ for simplicity. The space-time cylinder $Q_{\Omega_\theta}$ induced by $\Omega_\theta$ then reads:

$$ Q_{\Omega_\theta} := \Theta(Q_\Omega) = \{ (t, x) \in Q, x \in \varphi_t(\Omega_\theta) \}. $$

**Remark 3.1.** Loosely speaking, the space-time deformation mapping $\Theta$ transmits the effect of $\theta$ on the configuration at rest $\Omega$ (which is deformed into $\Omega_\theta$) to the rotated configurations $\Omega(t) = \varphi_t(\Omega)$ (which become $\varphi_t(\Omega_\theta)$). Contrary to the works [13, 14], we do not aim to find a time-dependent shape, but rather optimize the design of its configuration at rest with respect to a measure of performance depending on its known deformations $\varphi_t(\Omega)$. The action of $\Theta$ on the space-time cylinder $Q_\Omega$ is exemplified in Fig. 3.

The derivative of $\Theta(t, x)$ with respect to the variables $(t, x)$ is the following $(d+1) \times (d+1)$ matrix field $F(\theta): Q \to \mathbb{R}^{(d+1) \times (d+1)}$:

$$ F(\theta) = \begin{pmatrix} F_{tt}(\theta) & F_{tx}(\theta) \\ F_{xt}(\theta) & F_{xx}(\theta) \end{pmatrix}, $$

whose $1 \times 1$ and $1 \times d$ blocks $F_{tt}(\theta)$ and $F_{tx}(\theta)$ read:

$$ F_{tt}(\theta)(t, x) = \frac{\partial \Theta_t}{\partial t}(t, x) = 1, \quad F_{tx}(\theta)(t, x) = \nabla \Theta_t(t, x) = 0, $$
<figure><img src="image_4.png" /><figcaption>FIGURE 3. One deformation $\theta$ of the shape $\Omega$, with the associated deformation $\Theta$ of the space-time cylinder $Q_{\Omega}$.</figcaption></figure>

and whose $d \times 1$ and $d \times d$ blocks $F_{xt}(\theta)$ and $F_{xx}(\theta)$ are given by:

$$ 
\begin{aligned}
F_{xt}(\theta)(t, x) &= \frac{\partial \Theta_x}{\partial t}(t, x) \\
&= \frac{\partial \varphi_t}{\partial t}((\mathrm{Id} + \theta)(\varphi_t^{-1}(x))) + \left[\nabla \varphi_t((\mathrm{Id} + \theta)(\varphi_t^{-1}(x)))\right] \left[(\mathrm{Id} + \nabla \theta)(\varphi_t^{-1}(x))\right] \frac{\partial}{\partial t}(\varphi_t^{-1}(x)),
\end{aligned}
$$

$$ 
\begin{aligned}
F_{xx}(\theta)(t, x) &= \nabla \Theta_x(t, x) \\
&= \left[\nabla \varphi_t((\mathrm{Id} + \theta)(\varphi_t^{-1}(x)))\right] \left[(\mathrm{Id} + \nabla \theta)(\varphi_t^{-1}(x))\right] \left[\nabla \varphi_t^{-1}(x)\right].
\end{aligned}
$$

A simple calculation immediately yields the expression of the matrix field $F(\theta)^{-T}$:

$$ (3.7) \qquad F(\theta)^{-T} = \begin{pmatrix} 1 & b(\theta)^T \\ 0 & F_{xx}(\theta)^{-T} \end{pmatrix}, \quad \text{where } b(\theta) : Q \to \mathbb{R}^d \text{ is defined by } b(\theta) := -[F_{xx}(\theta)]^{-1} F_{xt}(\theta). $$

Let now $u : Q \to \mathbb{R}$ be a smooth enough function; the derivative of the composite function $v := u \circ \Theta : Q \to \mathbb{R}$ satisfies:

$$ \left( \frac{\partial u}{\partial t} \quad \nabla u \right) \circ \Theta = F(\theta)^{-T} \left( \frac{\partial v}{\partial t} \quad \nabla v \right), $$

and so

$$ (3.8) \qquad \left( \frac{\partial u}{\partial t} \right) \circ \Theta = \frac{\partial v}{\partial t} + b(\theta) \cdot \nabla v, \quad \text{and} \quad (\nabla u) \circ \Theta = [F_{xx}(\theta)]^{-T} \nabla v. $$

We now list a series of useful technical identities in the calculation of shape derivatives. These can be rigorously established along the lines of Chap. 5 in [36]; for brevity, we limit ourselves to formal statements and computations which in particular omit the functional setting.
• The vector field $b(0): Q \to \mathbb{R}^d$ equals:

$$
\begin{align*}
b(0)(t, x) &= -F_{xt}(0)(t, x) \\
&= -\frac{\partial \varphi_t}{\partial t}(\varphi_t^{-1}(x)) - [\nabla \varphi_t(\varphi_t^{-1}(x))] \frac{\partial}{\partial t}(\varphi_t^{-1}(x)) \\
&= -\frac{\partial}{\partial t}(\varphi_t(\varphi_t^{-1}(x))) \\
&= 0,
\end{align*}
$$

where we have used the fact that

$$
F_{xx}(0)(t, x) = [\nabla \varphi_t(\varphi_t^{-1}(x))] [\nabla \varphi_t^{-1}(x)] = \nabla(\varphi_t \circ \varphi_t^{-1}(x)) = \mathbf{I}.
$$

• It follows easily from the above two relations that the mapping $\theta \mapsto b(\theta)$ is Fréchet differentiable at $\theta = 0$, with derivative:

$$
(3.9) \qquad b'(0)(\theta) = -F'_{xt}(0)(\theta).
$$

• The derivatives of the mappings $\theta \mapsto F_{xt}(\theta)$ and $\theta \mapsto F_{xx}(\theta)$ at $\theta = 0$ read:

$$
\begin{equation}
\begin{split}
(3.10) \quad F'_{xt}(0)(\theta)(t,x) = \left[ \nabla \left( \frac{\partial \varphi_t}{\partial t} \right) (\varphi_t^{-1}(x)) \right] \theta(\varphi_t^{-1}(x)) &+ \left[ \nabla^2 \varphi_t(\varphi_t^{-1}(x)) \theta(\varphi_t^{-1}(x)) \right] \frac{\partial}{\partial t} (\varphi_t^{-1}(x)) \\
&+ \left[ \nabla \varphi_t(\varphi_t^{-1}(x)) \right] \left[ \nabla \theta(\varphi_t^{-1}(x)) \right] \frac{\partial}{\partial t} (\varphi_t^{-1}(x)),
\end{split}
\tag{and}
\end{equation}
$$

$$
(3.11) \quad F'_{xx}(0)(\theta)(t,x) = [\nabla^2 \varphi_t(\varphi_t^{-1}(x))\theta(\varphi_t^{-1}(x))] [\nabla \varphi_t^{-1}(x)] + [\nabla \varphi_t(\varphi_t^{-1}(x))] [\nabla \theta(\varphi_t^{-1}(x))] [\nabla \varphi_t^{-1}(x)].
$$

• The determinant $m(\theta) : Q \to \mathbb{R}$ of the space-time change of variables induced by $\Theta$ in (3.6),

$$
(3.12) \qquad m(\theta) := |\det F(\theta)|,
$$

is a Fréchet differentiable function of $\theta$ at 0 and:

$$
(3.13) \quad m'(0)(\theta)(t,x) := \text{tr} \left( \left[ \nabla^2 \varphi_t(\varphi_t^{-1}(x)) \theta(\varphi_t^{-1}(x)) \right] \left[ \nabla \varphi_t^{-1}(x) \right] \right) + (\text{div}\theta)(\varphi_t^{-1}(x)).
$$

• Let $f: Q \to \mathbb{R}$ be a smooth scalar function. The mapping $\theta \mapsto f \circ \Theta$ is Fréchet differentiable at $\theta = 0$, and its derivative $f_1(\theta)$ reads:

$$
(3.14) \qquad f_1(\theta)(t,x) := \left( \left[ \nabla \varphi_t(\varphi_t^{-1}(x)) \right]^T \nabla f(t,x) \right) \cdot \theta(\varphi_t^{-1}(x)).
$$

• Likewise, let $w: Q \to \mathbb{R}^d$ be a smooth vector field; then the mapping $\theta \mapsto w \circ \Theta$ is Fréchet differentiable at $\theta = 0$, and its derivative $w_1(\theta)$ reads

$$
(3.15) \qquad w_1(\theta)(t,x) := [\nabla w(t,x)] [\nabla \varphi_t(\varphi_t^{-1}(x))] \theta(\varphi_t^{-1}(x)).
$$

• The matrix-valued mapping $A(\theta) : Q \to \mathbb{R}^{d \times d}$ defined by

$$
(3.16) \qquad A(\theta) := m(\theta) F_{xx}^{-1}(\theta) F_{xx}^{-T}(\theta);
$$

is Fréchet differentiable at $\theta = 0$, with derivative:

$$
(3.17) \qquad A'(0)(\theta) = m'(0)(\theta)I - F'_{xx}(0)(\theta) - F'_{xx}(0)(\theta)^T.
$$

**Remark 3.2.** When the motion $\varphi_t$ is induced by a rotation as in (2.18), some of the above expressions can be conveniently simplified. Indeed, it holds in this case:

$$
\frac{dR_{\alpha(t)}}{dt} = \alpha'(t) \begin{pmatrix} -\sin(\alpha(t)) & -\cos(\alpha(t)) \\ \cos(\alpha(t)) & -\sin(\alpha(t)) \end{pmatrix} \quad \text{and} \quad \left[ \nabla v(t, x) \right] = \alpha'(t) \begin{pmatrix} 0 & -1 \\ 1 & 0 \end{pmatrix},
$$

and so:

$$
F'_{xx}(0)(\theta) = R_{\alpha(t)}[\nabla\theta(\varphi_t^{-1}(x))]R_{-\alpha(t)},
$$

$$
F'_{xt}(0)(\theta) = \left( \frac{dR_{\alpha(t)}}{dt} \right) \theta(\varphi_t^{-1}(x)) + R_{\alpha(t)} \left[ \nabla \theta(\varphi_t^{-1}(x)) \right] \left( \frac{dR_{-\alpha(t)}}{dt} \right),
$$
and for any smooth enough vector field $w: Q \to \mathbb{R}^d$,

$$w_1(\theta) = [\nabla w(t, x)] R_{\alpha(t)} \theta(\varphi_t^{-1}(x)).$$

### 3.2.2. An academic shape optimization problem

In this section, we calculate the derivative of the following shape functional:

$$(3.18) \qquad J(\Omega) = \int_{Q_{\Omega}} f(t, x) \, dxdt,$$

depending on the domain $\Omega$ via the integral of a smooth function $f \in C^{\infty}([0, T] \times \mathbb{R}^d)$ over the space-time cylinder $Q_{\Omega}$, which involves all the deformed version $\Omega(t)$ of $\Omega$, $t \in (0, T)$. Again, this quite academic question is a useful preliminary to the treatment of the more “physical” shape functionals described in Section 4.

**Lemma 3.1.** The functional $J(\Omega)$ in (3.18) is shape differentiable at any bounded Lipschitz shape $\Omega \subset \mathbb{R}^d$, and its shape derivative reads, in volume form:

$$(3.19) \quad J'(\Omega)(\theta) = \int_{Q_{\Omega}} \operatorname{tr} \left( \left[ \nabla^2 \varphi_t(\varphi_t^{-1}(x)) \theta(\varphi_t^{-1}(x)) \right] \left[ \nabla \varphi_t^{-1}(x) \right] \right) f(t, x) \, dxdt \\ + \int_{Q_{\Omega}} (\operatorname{div}\theta)(\varphi_t^{-1}(x)) f(t, x) \, dxdt + \int_{Q_{\Omega}} \left[ \nabla \varphi_t(\varphi_t^{-1}(x)) \right]^T \nabla f(t, x) \cdot \theta(\varphi_t^{-1}(x)) \, dxdt,$$

or equivalently, in surface form:

$$(3.20) \quad J'(\Omega)(\theta) = \int_{\partial\Omega} v_{\Omega} \theta \cdot n \, ds, \quad \text{where } v_{\Omega}(x) := \int_{0}^{T} |\det \nabla \varphi_t(x)| f(t, \varphi_t(x)) \, dt.$$

*Proof.* For any small enough deformation $\theta \in \Theta_{\text{ad}}$, the definition of $J(\Omega_{\theta})$ reads:

$$J(\Omega_{\theta}) = \int_{Q_{\Omega_{\theta}}} f(t, x) \, dxdt.$$

Using a change of variables based on the mapping $\Theta$ in (3.6) and the fact that, by definition $\Theta^{-1}(Q_{\Omega_{\theta}}) = Q_{\Omega}$, we obtain:

$$J(\Omega_{\theta}) = \int_{Q_{\Omega}} m(\theta)(f \circ \Theta)(t, x) \, dxdt,$$

where we recall the notation $m(\theta) = |\det F(\theta)|$. Combining the expansions (3.13) and (3.14) of the two terms in the above integrand, we see that the mapping $\theta \mapsto J(\Omega_{\theta})$ is differentiable at $\theta = 0$, with derivative

$$J'(\Omega)(\theta) = \int_{Q_{\Omega}} \operatorname{tr} \left( \left[ \nabla^2 \varphi_t(\varphi_t^{-1}(x)) \theta(\varphi_t^{-1}(x)) \right] \left[ \nabla \varphi_t^{-1}(x) \right] \right) f(t, x) \, dxdt \\ + \int_{Q_{\Omega}} (\operatorname{div}\theta)(\varphi_t^{-1}(x)) f(t, x) \, dxdt + \int_{Q_{\Omega}} \left[ \nabla \varphi_t(\varphi_t^{-1}(x)) \right]^T \nabla f(t, x) \cdot \theta(\varphi_t^{-1}(x)) \, dxdt,$$

which is the announced formula (3.19).

To infer the surface expression (3.20), we first express (3.19) as nested integrals over $(0, T)$ and $\varphi_t(\Omega)$ and we use a change of variables in terms of the mapping $\varphi_t$ in the spatial integral to obtain

$$(3.21) \quad J'(\Omega)(\theta) = \int_{0}^{T} \int_{\Omega} |\det \nabla \varphi_t(x)| \operatorname{tr} \left( \left[ \nabla^2 \varphi_t(x) \theta(x) \right] \left[ \nabla \varphi_t^{-1}(\varphi_t(x)) \right] \right) f(t, \varphi_t(x)) \, dxdt \\ + \int_{0}^{T} \int_{\Omega} |\det \nabla \varphi_t(x)| (\operatorname{div}\theta)(x) f(t, \varphi_t(x)) \, dxdt + \int_{0}^{T} \int_{\Omega} |\det \nabla \varphi_t(x)| \left[ \nabla \varphi_t(x) \right]^T \nabla f(t, \varphi_t(x)) \cdot \theta(x) \, dxdt.$$

We now proceed according to the strategy summarized in e.g. [2]: integrating by parts in all the terms of the above expression containing derivatives of $\theta$ yields an expression of the form:

$$J'(\Omega)(\theta) = \int_{\Omega} M_{\Omega} \cdot \theta \, dx + \int_{\Gamma} v_{\Omega}(\theta \cdot n) \, ds + \int_{\Gamma} t_{\Omega} \cdot \theta_{\Gamma} \, ds,$$
where $M_\Omega: D \to \mathbb{R}^d$, $v_\Omega: \Gamma \to \mathbb{R}$ and $t_\Omega: \Gamma \to \mathbb{R}^d$ depend on the function $f$ and on $\Omega$. If we expect that $J'(\Omega)(\theta)$ admits an expression with the desirable structure (3.5), the vector fields $M_\Omega$ and $t_\Omega$ must vanish. That this indeed holds is eventually verified by an elementary, albeit tedious calculation.

In the present case, an integration by parts in the second integral on the right-hand side of (3.21) (which is the only one involving derivatives of $\theta$) produces the following expression:

$$J'(\Omega)(\theta) = \int_0^T \int_\Gamma |\det \nabla \varphi_t(x)| f(t, \varphi_t(x)) \theta(x) \cdot n(x) \, ds(x)dt + r(\theta),$$

where $r(\theta)$ is a collection of integrals posed on the domain $\Omega$ involving only $\theta$ (and not its derivatives). A simple, albeit tedious calculation which is omitted for brevity reveals that $r(\theta)$ actually vanishes, which allows to conclude. □

## 4. CALCULATION OF THE SHAPE DERIVATIVE IN THE MAGNETO-QUASI-STATIC CONTEXT

In this section, we turn to the calculation of the shape derivative of a “physical” objective function $J(\Omega)$ which depends on the shape $\Omega$ via the potential $u_\Omega$, solution to the magneto-quasi-static problem (2.21) where $\Omega$ represents the phase occupied by ferromagnetic material. After sketching some mathematical preliminaries about this evolution problem in Section 4.1 and a few technical features about its solution in Section 4.2, we proceed to the calculation of the shape derivative $J'(\Omega)(\theta)$ in Section 4.3, under the simplifying assumption that the reluctivity $\nu$ of the ferromagnetic material is constant, in which case (2.21) becomes linear. We eventually provide the formal extension of this result in the general situation of a non linear material in Section 4.4.

### 4.1. Mathematical preliminaries about the magneto quasi-static evolution problem

The analysis of the evolution problem (2.21) brings into play suitable functional spaces for functions depending on the time and space variables. For the sake of convenience, a few basic facts about these are recalled in Appendix B.1; we refer to e.g. §1 in Chap. XVIII of [21], or Chap. 23 in [74] for more exhaustive presentations.

Let us consider the functional space $L^2(0, T; H_0^1(D))$, whose dual $L^2(0, T; H_0^1(D))^*$ can be identified with $L^2(0, T; H^{-1}(D))$, see (B.2). We also introduce

$$W := \left\{ u \in L^2(0, T; H_0^1(D)) \text{ s.t. } \sigma_{\Omega(t)} \frac{\partial u}{\partial t} \in L^2(0, T; H^{-1}(D)) \right\},$$

which is a Hilbert space when equipped with the norm

$$\|u\|_W^2 := \int_0^T \|u(t, \cdot)\|_{H_0^1(D)}^2 dt + \int_0^T \left\| \sigma_{\Omega(t)} \frac{\partial u}{\partial t}(t, \cdot) \right\|_{H^{-1}(D)}^2 dt.$$

Intuitively, $W$ is a variant of the more classical space $W(0, T; H_0^1(D), L^2(D))$ defined in (B.3), which is adapted to handle the fact that the coefficient $\sigma_{\Omega(t)}$ vanishes outside the region $D_{\text{mag}}(t) \subset D$, see (2.11). The properties of $W(0, T; H_0^1(D), L^2(D))$ summarized in Lemma B.1 are also verified by $W$, as the proof of this result can be straightforwardly adapted to the present case. In particular, any element $u \in W$ induces a continuous mapping $[0, T] \ni t \mapsto u(t, \cdot) \in L^2(D_{\text{mag}})$, which allows to introduce the closed subspace $W_{\text{per}}$ of $W$ defined by:

$$(4.1) \qquad W_{\text{per}} := \{u \in W \text{ s.t. } u(t=0, \cdot) = u(t=T, \cdot) \text{ on } D_{\text{mag}}\}.$$

**Remark 4.1.**

* As the conductivity $\sigma_{\Omega(t)}$ equals $\sigma_m > 0$ on $D_{\text{mag}}(t)$ and vanishes outside $D_{\text{mag}}(t)$, the time derivative $\frac{\partial}{\partial t}(u(t, \varphi_t(\cdot)))$ of a function $u \in W$ belongs to $L^2(0, T; H^1(D_{\text{mag}})^*)$.

* The time periodicity condition satisfied by functions $u \in W_{\text{per}}$ only holds in $D_{\text{mag}}$, as an application of Lemma B.1 (ii).
Using this language, the magneto-quasi-static problem (2.21) can be expressed in variational form:

$$ (4.2) \quad \text{Search for } u_{\Omega} \in W_{\text{per}}, \text{ s.t. } \forall w \in L^2(0, T; H_0^1(D)), $$

$$ \int_Q \sigma_{\Omega(t)} \left( \frac{\partial u_{\Omega}}{\partial t} + v \cdot \nabla u_{\Omega} \right) w \, dxdt + \int_Q \nu_{\Omega(t)}(x, |\nabla u_{\Omega}|) \nabla u_{\Omega} \cdot \nabla w \, dxdt = \int_Q fw \, dxdt, $$

or equivalently:

$$ (4.3) \quad \text{Search for } u_{\Omega} \in W_{\text{per}}, \text{ s.t. for a.e. } t \in (0, T), $$

$$ \forall w \in H_0^1(D), \quad \int_D \sigma_{\Omega(t)} \left( \frac{\partial u_{\Omega}}{\partial t} + v \cdot \nabla u_{\Omega} \right) w \, dx + \int_D \nu_{\Omega(t)}(x, |\nabla u_{\Omega}|) \nabla u_{\Omega} \cdot \nabla w \, dx = \int_D fw \, dx. $$

The following theorem deals with the well-posedness of the variational problem (4.2).

**Theorem 4.1.** Let $\varphi_t$ be a time-dependent diffeomorphism of $\bar{D}$ satisfying (2.14) to (2.16). For any source $f \in L^2(0, T; H^{-1}(D))$, the evolution problem (4.2) has a unique solution $u_{\Omega} \in W_{\text{per}}$, which has a Lipschitz dependence on $f$: there exists a constant $C > 0$ such that the solutions $u_1, u_2 \in W_{\text{per}}$ associated to the respective right-hand sides $f_1, f_2 \in L^2(0, T; H^{-1}(D))$ satisfy:

$$ \|u_1 - u_2\|_{W} \le C \|f_1 - f_2\|_{L^2(0,T;H^{-1}(D))}. $$

The proof of this result is sketched in Appendix B. Let us point out that the analysis of a similar statement about a linear version of (4.2) was presented recently in [28].

## 4.2. A few properties of the solution to the magneto-quasi-static problem

This section is devoted to a few useful technical facts about functions in the space $W$, and more precisely about the solution $u_{\Omega}$ to (4.2). Throughout, $\varphi_t$ is a smooth time-dependent diffeomorphism of $\bar{D}$ satisfying (2.14) to (2.16) and the coefficients $\sigma_{\Omega(t)}$ and $\nu_{\Omega(t)}$ are given by (2.11).

The first remark is an avatar of the so-called Reynolds (or transport) theorem:

**Lemma 4.1.** Let $u, p \in W$; the following identity holds in the sense of distributions on $(0, T)$:

$$ \frac{d}{dt} \left( \int_D \sigma_{\Omega(t)}(x) u(t, x) p(t, x) \, dx \right) = \int_D \sigma_{\Omega(t)}(x) \frac{du}{dt}(t, x) p(t, x) \, dx \\ + \int_D \sigma_{\Omega(t)}(x) u(t, x) \frac{dp}{dt}(t, x) \, dx + \int_Q \sigma_{\Omega(t)} \operatorname{div} v(t, x) u(t, x) p(t, x) \, dxdt, $$

where we recall the notation $\frac{du}{dt}$, $\frac{dp}{dt}$ in (2.20) for the total time derivatives of $u$ and $p$.

*Proof.* By Lemma B.1 (i), it is enough to prove that the desired identity holds at every $t \in (0, T)$ when $u$ and $p$ belong to the space $C^{\infty}([0, T], H_0^1(D))$. To achieve this, we first use a change of variables based on the mapping $\varphi_t$:

$$ \int_D \sigma_{\Omega(t)}(x) u(t, x) p(t, x) \, dx = \int_D |\det \nabla \varphi_t(x)| \sigma_{\Omega(t)}(\varphi_t(x)) u(t, \varphi_t(x)) p(t, \varphi_t(x)) \, dx \\ = \int_D |\det \nabla \varphi_t(x)| \sigma_{\Omega}(x) u(t, \varphi_t(x)) p(t, \varphi_t(x)) \, dx, $$

where the second line follows from (2.17). Taking derivatives with respect to time, we obtain:

$$ \frac{d}{dt} \left( \int_D \sigma_{\Omega(t)}(x) u(t, x) p(t, x) \, dx \right) = \int_D |\det \nabla \varphi_t(x)| \sigma_{\Omega}(x) \frac{du}{dt}(t, \varphi_t(x)) p(t, \varphi_t(x)) \, dx \\ + \int_D |\det \nabla \varphi_t(x)| \sigma_{\Omega}(x) u(t, \varphi_t(x)) \frac{dp}{dt}(t, \varphi_t(x)) \, dx \\ + \int_D |\det \nabla \varphi_t(x)| \operatorname{tr} \left( \nabla \varphi_t(x)^{-1} \nabla \left( \frac{\partial \varphi_t}{\partial t} \right)(x) \right) \sigma_{\Omega}(x) u(t, \varphi_t(x)) p(t, \varphi_t(x)) \, dx. $$
Here, we have used the usual formula for the derivative of the determinant mapping. Another change of variables now yields:

$$ \frac{d}{dt} \left( \int_D \sigma_{\Omega(t)}(x) u(t, x) p(t, x) \, dx \right) = \int_D \sigma_{\Omega(t)}(x) \frac{du}{dt}(t, x) p(t, x) \, dx + \int_D \sigma_{\Omega(t)}(x) u(t, x) \frac{dp}{dt}(t, x) \, dx \\ + \int_D \mathrm{tr} \left( \nabla \varphi_t (\varphi_t^{-1}(x))^{-1} \nabla \left( \frac{\partial \varphi_t}{\partial t} \right) (\varphi_t^{-1}(x)) \right) \sigma_{\Omega(t)}(x) u(t, x) p(t, x) \, dx. $$

The desired identity eventually results from the following elementary calculation and the definition (2.13) of the velocity field $v(t, x)$:

$$ \begin{aligned} \nabla \left( \frac{\partial \varphi_t}{\partial t} \right) (\varphi_t^{-1}(x)) \nabla \varphi_t (\varphi_t^{-1}(x))^{-1} &= \nabla \left( \frac{\partial \varphi_t}{\partial t} \right) (\varphi_t^{-1}(x)) \nabla (\varphi_t^{-1}(x)) \\ &= \nabla \left( \frac{\partial \varphi_t}{\partial t} (\varphi_t^{-1}(x)) \right) \\ &= \nabla v(t, x). \end{aligned} $$

□

**Corollary 4.1.** The following relation holds, for all $u, p \in W$:

$$ \int_Q \sigma_{\Omega(t)} \frac{du}{dt} p \, dx dt = \int_D \left( \sigma_{\Omega(T)} u(T, x) p(T, x) - \sigma_{\Omega(0)} u(0, x) p(0, x) \right) \, dx \\ - \int_Q \sigma_{\Omega(t)} \frac{dp}{dt} u \, dx dt - \int_Q \sigma_{\Omega(t)} (\mathrm{div} v) u p \, dx dt. $$

*Proof.* For a.e. $t \in (0, T)$, Lemma 4.1 states that:

$$ \int_D \sigma_{\Omega(t)}(x) \frac{du}{dt}(t, x) p(t, x) \, dx = \frac{d}{dt} \left( \int_D \sigma_{\Omega(t)}(x) u(t, x) p(t, x) \, dx \right) \\ - \int_D \sigma_{\Omega(t)}(x) u(t, x) \frac{dp}{dt}(t, x) \, dx - \int_D \sigma_{\Omega(t)} \mathrm{div} v(t, x) u(t, x) p(t, x) \, dx. $$

In particular, the mapping $t \mapsto \int_D \sigma_{\Omega(t)} u(t, x) p(t, x) \, dx$ is absolutely continuous on $(0, T)$. The result follows immediately from integration of this identity over the time period $(0, T)$. □

The variational problem (4.2) implicitly encompasses jump conditions for $u_{\Omega}(t, \cdot)$ at the interface $\Gamma(t)$ where the conductivity and the reluctivity are discontinuous; these are the subject of the next lemma:

**Lemma 4.2.** Let the right-hand side $f$ belong to the space $L^2(0, T, L^2(D))$; the corresponding solution $u_{\Omega}(t, \cdot)$ to (4.2) is continuous across $\Gamma(t)$:

$$ (4.4) \qquad \forall t \in (0, T), \quad x \in \Gamma, \quad u_{\Omega}^{-}(t, \varphi_{t}(x)) = u_{\Omega}^{+}(t, \varphi_{t}(x)), $$

and so are its normal flux and tangential gradient:

$$ (4.5) \qquad \forall t \in (0, T), \quad x \in \Gamma, \quad \left( \nu_{\Omega} \nabla u_{\Omega}(t, \varphi_t(x)) \right)^{-} \cdot n_{\Omega(t)}(\varphi_t(x)) = \left( \nu_{\Omega} \nabla u_{\Omega}(t, \varphi_t(x)) \right)^{+} \cdot n_{\Omega(t)}(\varphi_t(x)), $$

$$ (4.6) \qquad \forall t \in (0, T), \quad x \in \Gamma, \quad \nabla_{\Gamma(t)} u_{\Omega}^{-}(t, \varphi_{t}(x)) = \nabla_{\Gamma(t)} u_{\Omega}^{+}(t, \varphi_{t}(x)). $$

Finally, it holds:

$$ (4.7) \qquad \forall t \in (0, T), \quad x \in \Gamma, \quad \frac{du_{\Omega}^{-}}{dt}(t, \varphi_t(x)) = \frac{du_{\Omega}^{+}}{dt}(t, \varphi_t(x)). $$

*Proof.* Since $u_{\Omega} \in L^2(0, T; H_0^1(D))$, the function $u_{\Omega}(t, \cdot)$ belongs to $H_0^1(D)$ for a.e. $t \in (0, T)$, and so (4.4) holds true.

Let now $t \in (0, T)$ be given; by considering smooth test functions $w \in C_c^{\infty}(D)$ with compact support inside $\Omega(t_0)$ or $D \setminus \Omega(t_0)$ in (4.3), then integrating by parts, we see that:

$$ (4.8) \qquad \sigma_{\Omega(t)} \frac{du_{\Omega}}{dt}(t, x) - \mathrm{div} \left( \nu_{\Omega}(t) (x, |\nabla u_{\Omega}(t, x)|) \nabla u_{\Omega}(t, x) \right) = f(t, x) $$
in the sense of distributions in $\Omega(t)$ or $D \setminus \overline{\Omega}(t)$. Repeating the calculation with arbitrary smooth test functions $w \in C_c^\infty(\bar{D})$ and using (4.8), we now obtain:

$$ \int_{\Gamma(t)} \left( \nu_{\Omega(t)} \frac{\partial u_{\Omega}}{\partial n} \right)^{-} w \, ds - \int_{\Gamma(t)} \left( \nu_{\Omega(t)} \frac{\partial u_{\Omega}}{\partial n} \right)^{+} w \, ds = 0, $$

and (4.5) follows immediately.

Eventually, (4.6) is obtained by taking the tangential gradient in (4.4), and taking time derivatives in (4.4) yields the last jump relation (4.7). □

**Remark 4.2.** In the following, we denote by $\nabla_{\Gamma(t)}u_{\Omega}(t, \varphi_t(x))$ and $\nu_{\Omega(t)}\nabla u_{\Omega}(t, \varphi_t(x)) \cdot n_{\Omega(t)}(\varphi_t(x))$ the common values of the one-sided limits of the tangential gradient and normal flux of $u_{\Omega}(t, \cdot)$ through $\Gamma(t)$.

### 4.3. Calculation of the shape derivative of a physical objective function $J(\Omega)$ in the linear magneto-quasi-static setting

In this section, we detail the calculation of the shape derivative of the physical objective function $J(\Omega)$ in (3.2) involving the magneto-quasi-static potential $u_{\Omega}$, with a particular emphasis on the arguments which are quite specific to the present physical context. To keep technicality at a minimum, our analysis is conducted in the linear case where the reluctivity $\hat{\nu}$ of the ferromagnetic material is a constant function $\hat{\nu} \equiv \nu_f$, see Remark 2.3. The needed adaptations to the general case of a non constant function $\hat{\nu}$ are elementary, albeit tedious; they are discussed in the next Section 4.4. Throughout the remainder of this section, we rely on the following shorthand in notation: we denote by $\Omega_1 := \Omega$ the optimized phase of the rotor filled by ferromagnetic material, and by $\Omega_2 := D \setminus \bar{\Omega}$ its complement in $D$. Accordingly, we denote by $\nu_1, \sigma_1$ the reluctivity and conductivity inside $\Omega_1$, and by $\nu_2, \sigma_2$ the piecewise constant functions matching the repartition (2.11).

The evolution problem (2.21) then rewrites:

$$ (4.9) \qquad \left\{ \begin{array}{ll} \sigma_{\Omega(t)} \frac{du_{\Omega}}{dt} - \text{div}(\nu_{\Omega(t)} \nabla u_{\Omega}) = f & \text{in } Q, \\ u_{\Omega}(t, x) = 0 & \text{for } t \in (0, T), x \in \partial D, \\ u_{\Omega}(0, x) = u_{\Omega}(T, x) & \text{for } x \in D_{\text{mag}}, \end{array} \right. $$

where

$$ \sigma_{\Omega(t)}(x) = \begin{cases} \sigma_1 & \text{if } x \in \Omega(t), \\ \sigma_2 & \text{otherwise,} \end{cases} \quad \text{and} \quad \nu_{\Omega(t)}(x) = \begin{cases} \nu_1 & \text{if } x \in \Omega(t), \\ \nu_2 & \text{otherwise.} \end{cases} $$

**Proposition 4.1.** Let the source term $f: Q \to \mathbb{R}$ be smooth; the functional $J(\Omega)$ in (3.2) is shape differentiable at any bounded and Lipschitz shape $\Omega \subset D_{\text{rot}}$, and its shape derivative reads, under volume form:

$$ (4.10) \quad \forall \theta \in \Theta_{\text{ad}}, \quad J'(\Omega)(\theta) = \int_Q m'(0)(\theta) j(u_{\Omega}) \, dxdt \\ \qquad + \int_Q \sigma_{\Omega(t)} \left( m'(0)(\theta) \frac{du_{\Omega}}{dt} - [F'_{xx}(0)(\theta)] v \cdot \nabla u_{\Omega} + v_1(\theta) \cdot \nabla u_{\Omega} + b'(0)(\theta) \cdot \nabla u_{\Omega} \right) p_{\Omega} \, dxdt \\ \qquad + \int_Q \nu_{\Omega(t)} [A'(0)(\theta)] \nabla u_{\Omega} \cdot \nabla p_{\Omega} \, dxdt - \int_Q (m'(0)(\theta)f + f_1(\theta)) p_{\Omega} \, dxdt, $$

where the quantities $A'(0)(\theta)$, $m'(0)(\theta)$, $F'_{xx}(0)(\theta)$ and $b'(0)(\theta)$ are respectively defined by (3.9), (3.11), (3.13) and (3.17), and $f_1(\theta)$ and $v_1(\theta)$ are the derivatives of the mappings $\theta \mapsto f \circ \Theta$ and $\theta \mapsto v \circ \Theta$, see (3.14) and (3.15).

The adjoint state $p_{\Omega}$ is the unique solution in $W_{\text{per}}$ to the evolution problem:

$$ (4.11) \qquad \left\{ \begin{array}{ll} -\sigma_{\Omega(t)} \frac{dp_{\Omega}}{dt} - \text{div}(\nu_{\Omega(t)} \nabla p_{\Omega}) - \sigma_{\Omega(t)}(\text{div}v)p_{\Omega} = -j'(u_{\Omega}) & \text{in } Q, \\ p_{\Omega}(t, x) = 0 & \text{for } t \in (0, T), x \in \partial D, \\ p_{\Omega}(0, x) = p_{\Omega}(T, x) & \text{for } x \in D. \end{array} \right. $$
It is equivalently characterized as the solution to the following variational problem:

$$
(4.12) \quad \text{Search for } p_{\Omega} \in W_{\text{per}} \text{ s.t. } \forall w \in L^2(0, T; H_0^1(D)), \\
\qquad \qquad \int_Q \left( -\sigma_{\Omega(t)} \frac{dp_{\Omega}}{dt} w + \nu_{\Omega(t)} \nabla p_{\Omega} \cdot \nabla w - \sigma_{\Omega(t)} (\text{div} v) p_{\Omega} w \right) dxdt = - \int_Q j'(u_{\Omega}) w dxdt.
$$

Assuming $\Omega$ to be smooth, the shape derivative $J'(\Omega)(\theta)$ admits the surface form:

$$
\begin{equation} \tag{4.13}
\begin{split}
\forall \theta \in \Theta_{\text{ad}}, \quad J'(\Omega)(\theta) &= \int_{\Gamma} v_{\Omega}(\theta \cdot n) \, \text{d}s, \text{ where} \\
& \qquad v_{\Omega}(x) := -(\sigma_2 - \sigma_1) \int_0^T |\det \nabla \varphi_t(x)| \frac{du_{\Omega}}{dt}(t, \varphi_t(x)) p_{\Omega}(t, \varphi_t(x)) \, \text{d}t \\
& \qquad \quad - (\nu_2 - \nu_1) \int_0^T |\det \nabla \varphi_t(x)| |\nabla_{\Gamma(t)} u_{\Omega}(t, \varphi_t(x)) \cdot \nabla_{\Gamma(t)} p_{\Omega}(t, \varphi_t(x))| \, \text{d}t \\
& \qquad + \left( \frac{1}{\nu_2} - \frac{1}{\nu_1} \right) \int_0^T |\det \nabla \varphi_t(x)| \left( \nu_{\Omega(t)} \nabla u_{\Omega}(t, \varphi_t(x)) \cdot n_{\Omega(t)}(\varphi_t(x)) \right) \left( \nu_{\Omega(t)} \nabla p_{\Omega}(t, \varphi_t(x)) \cdot n_{\Omega(t)}(\varphi_t(x)) \right) \, \text{d}t.
\end{split}
\end{equation}
$$

Remark 4.3.

* The well-posedness of the (linear) adjoint problem (4.11) follows from similar, although simpler arguments to those involved in the proof of Theorem 4.1 about the (non linear) state problem (2.21).

* Like $u_{\Omega}(t, \cdot)$, the adjoint state $p_{\Omega}(t, \cdot)$ satisfies jump relations at the interface $\Gamma(t)$; a similar analysis to that in the proof of Lemma 4.2 indeed yields:

$$
\forall t \in (0, T), x \in \Gamma, \quad p_{\Omega}^{+}(t, \varphi_t(x)) = p_{\Omega}^{+}(t, \varphi_t(x)),
$$

and

$$
\forall t \in (0, T), x \in \Gamma, \quad (\nu_{\Omega} \nabla p_{\Omega}(t, \varphi_t(x)))^{-} \cdot n_{\Omega(t)}(\varphi_t(x)) = (\nu_{\Omega} \nabla p_{\Omega}(t, \varphi_t(x)))^{+} \cdot n_{\Omega(t)}(\varphi_t(x)).
$$

Remark 4.4.

* The smoothness assumption about the source $f$ in the statement of Proposition 4.1 can be weakened in various ways. For instance, the statement holds mutatis mutandis if $f$ only belongs to $L^2(0, T; H^{-1}(D))$, but for each time $t \in (0, T)$, $f(t, \cdot)$ has support in the region $D_{\text{stat}}$ where deformations $\theta \in \Theta_{\text{ad}}$ vanish.

* As we have evoked in Remark 2.3, in realistic applications, the right-hand side $f$ contains a contribution of the form $\text{div}(L1_{D_{\text{mag}}(t)})$ where the smooth vector field $L(t, x)$ stands for the 90° clockwise rotate $M^{\perp}$ of the magnetization $M$. In such case, the last integral in the right-hand side of (4.10) should be supplemented with:

$$
(4.14) \quad - \int_{Q_{D_{\text{mag}}}} \left( m'(0)(\theta)L + L_1(\theta) - [F'_{xx}(0)(\theta)]L \right) \cdot \nabla p_{\Omega} \,dxdt,
$$

where $L_1(\theta)$ is the first-order term in the asymptotic expansion of the vector field $L \circ \Theta$ at $\theta = 0$, see (3.15).

Proof of Proposition 4.1. However technically tedious, the proof follows a quite classical trail, as exposed for instance in [2, 54, 36]. We proceed in five steps.

Step 1: We prove that the mapping $\Omega \mapsto u_{\Omega}$ has a "Lagrangian derivative".

More precisely, we prove that the transported potential $\theta \mapsto \overline{u}_{\Omega}(\theta) := u_{\Omega_{\theta}} \circ \Theta$ is Fréchet differentiable from a neighborhood of 0 in $\Theta_{\text{ad}}$ into $W_{\text{per}}$. This starts from the variational problem satisfied by the magnetostatic quasi-static field $u_{\Omega_{\theta}} \in W_{\text{per}}$:

$$
(4.15) \quad \forall w \in L^2(0, T; H_0^1(D)), \quad \int_Q \sigma_{\varphi_t}(\Omega_\theta)(x) \frac{du_{\Omega_\theta}}{dt} w \,dxdt + \int_Q \nu_{\varphi_t}(\Omega_\theta)(x) \nabla u_{\Omega_\theta} \cdot \nabla w \,dxdt = \int_Q fw \,dxdt.
$$
A change of variables via $\Theta$ in the above integrals yields:

$$ (4.16) \quad \int_Q (\sigma_{\varphi_t(\Omega_\theta)} \circ \Theta_x) m(\theta) \left( \left( \frac{\partial u_{\Omega_\theta}}{\partial t} \right) \circ \Theta + (v \circ \Theta) \cdot (\nabla u_{\Omega_\theta} \circ \Theta) \right) (w \circ \Theta) \, dxdt \\ + \int_Q (\nu_{\varphi_t(\Omega_\theta)} \circ \Theta_x) m(\theta) \left( \nabla u_{\Omega_\theta} \circ \Theta \right) \cdot \left( \nabla w \circ \Theta \right) \, dxdt = \int_Q m(\theta)(f \circ \Theta)(w \circ \Theta) \, dxdt, $$

where we recall the expression (3.12) of $m(\theta)$. Now, from the definitions (2.17) of $\sigma_\Omega$ and $\nu_\Omega$, it holds:

$$ \sigma_{\varphi_t(\Omega_\theta)} \circ \Theta_x = \sigma_{\Omega(t)}, \quad \text{and} \quad \nu_{\varphi_t(\Omega_\theta)} \circ \Theta_x = \nu_{\Omega(t)}, $$

so that, by using the relations (3.8) and taking test functions of the form $w \circ \Theta^{-1}$, $w \in L^2(0, T; H_0^1(D))$, we see that the transported function $\overline{u_\Omega}(\theta)$ satisfies the variational problem:

$$ (4.17) \quad \forall w \in L^2(0, T; H_0^1(D)), \quad \int_Q \sigma_{\Omega(t)} m(\theta) \left( \frac{\partial \overline{u_\Omega}(\theta)}{\partial t} + b(\theta) \cdot \nabla \overline{u_\Omega}(\theta) + [F_{xx}^{-1}(\theta)](v \circ \Theta) \cdot \nabla \overline{u_\Omega}(\theta) \right) w \, dxdt \\ + \int_Q \nu_{\Omega(t)} [A(\theta)] \nabla \overline{u_\Omega}(\theta) \cdot \nabla w \, dxdt = \int_Q m(\theta)(f \circ \Theta)w \, dxdt, $$

where $A(\theta)$ and $b(\theta)$ are respectively defined by (3.16) and (3.7).

This problem can be reformulated as a more abstract equation for $\overline{u_\Omega}(\theta)$:

Search for $u \in W_{\text{per}} \text{ s.t. } \mathcal{F}(\theta, u) = 0,$

where we have defined the mapping $\mathcal{F}: \Theta_{\text{ad}} \times W_{\text{per}} \mapsto L^2(0, T; H^{-1}(D))$ by:

$$ (4.18) \quad \mathcal{F}(\theta, u) = \left[ \begin{array}{l} w \mapsto \int_Q \sigma_{\Omega(t)} m(\theta) \left( \frac{\partial u}{\partial t} + b(\theta) \cdot \nabla u + [F_{xx}^{-1}(\theta)](v \circ \Theta) \cdot \nabla u \right) w \, dxdt \\ \qquad + \int_Q \nu_{\Omega(t)} [A(\theta)] \nabla u \cdot \nabla w \, dxdt - \int_Q m(\theta)(f \circ \Theta)w \, dxdt \end{array} \right]. $$

An elementary verification shows that this equation fulfills the assumptions of the implicit function theorem at the point $(0, u_\Omega)$, see e.g. [45] Chap. I, Th. 5.9. Hence, the mapping $\theta \mapsto \overline{u_\Omega}(\theta)$ is Fréchet differentiable at $\theta = 0$. Its derivative, denoted by $\dot{u}_\Omega(\theta)$, is the so-called Lagrangian derivative of $u_\Omega$.

**Step 2:** We characterize the Lagrangian derivative $\dot{u}_\Omega(\theta) \in W_{\text{per}}$.

Relying on the information that the mapping $\theta \mapsto \overline{u_\Omega}(\theta)$ is differentiable at $\theta = 0$, taking derivatives in (4.17) shows that the function $\dot{u}_\Omega(\theta) \in W_{\text{per}}$ satisfies, for all $w \in L^2(0, T; H_0^1(D))$:

$$ (4.19) \quad \int_Q \sigma_{\Omega(t)} \frac{\mathrm{d}\dot{u}_\Omega(\theta)}{\mathrm{d}t} w \, dxdt + \int_Q \nu_{\Omega(t)} \nabla \dot{u}_\Omega(\theta) \cdot \nabla w \, dxdt = - \int_Q \sigma_{\Omega(t)} m'(0)(\theta) \frac{\mathrm{d}u_\Omega}{\mathrm{d}t} w \, dxdt \\ + \int_Q \sigma_{\Omega(t)} [F'_{xx}(0)(\theta)] v \cdot \nabla u_\Omega w \, dxdt - \int_Q \sigma_{\Omega(t)} v_1(\theta) \cdot \nabla u_\Omega w \, dxdt - \int_Q \sigma_{\Omega(t)} b'(0)(\theta) \cdot \nabla u_\Omega w \, dxdt \\ - \int_Q \nu_{\Omega(t)} [A'(0)(\theta)] \nabla u_\Omega \cdot \nabla w \, dxdt + \int_Q (m'(0)(\theta) f + f_1(\theta)) w \, dxdt. $$

This is a well-posed evolution problem for $\dot{u}_\Omega(\theta) \in W_{\text{per}}$.

**Step 3:** We calculate the derivative of the functional $J(\Omega)$ in (3.2).

A change of variables based on the mapping $\Theta$ in the definition of $J(\Omega_\theta)$ yields the following expression:

$$ J(\Omega_\theta) = \int_Q m(\theta) j(\overline{u_\Omega}(\theta)) \, dxdt, $$

which naturally brings into play the transported function $\overline{u_\Omega}(\theta)$. As we have proved the differentiability of the latter with respect to the deformation $\theta$, we now obtain, by taking derivatives in the above expression:

$$ (4.20) \quad J'(\Omega)(\theta) = \int_Q \left( m'(0)(\theta) j(u_\Omega) + j'(u_\Omega) \dot{u}_\Omega(\theta) \right) dxdt. $$
Step 4: We apply the adjoint method to derive the volume form of $J'(\Omega)(\theta)$.

The expression (4.20) is not suitable for practical calculations, as it depends in a complicated way on $\theta$, via the Lagrangian derivative $\dot{u}_{\Omega}(\theta)$, which is only known implicitly, as the solution to (4.19). In particular, (4.20) does not allow for the simple identification of a descent direction for $J(\Omega)$.

To overcome this classical issue, we rely on the adjoint method from optimal control theory, see e.g. [50] or [60] for an intuitive presentation. We introduce the adjoint state $p_{\Omega} \in W_{\text{per}}$ as the solution to the variational problem (4.12), and take $w = \dot{u}_{\Omega}(\theta)$ as test function in there to obtain:

$$ - \int_Q j'(u_{\Omega}) \dot{u}_{\Omega}(\theta) \, dxdt = \int_Q \left( -\sigma_{\Omega(t)} \frac{dp_{\Omega}}{dt} \dot{u}_{\Omega}(\theta) + \nu_{\Omega(t)} \nabla p_{\Omega} \cdot \nabla \dot{u}_{\Omega}(\theta) - \sigma_{\Omega(t)} (\text{div} v) p_{\Omega} \dot{u}_{\Omega}(\theta) \right) \, dxdt. $$

Using Corollary 4.1 and the fact that $\dot{u}_{\Omega}(\theta) \in W_{\text{per}}$ together with $\sigma_{\Omega(T)} = \sigma_{\Omega(0)}$ due to (2.15) and (2.16), it follows:

$$ \int_Q \left( \sigma_{\Omega(t)} \frac{du_{\Omega}(\theta)}{dt} p_{\Omega} + \nu_{\Omega(t)} \nabla p_{\Omega} \cdot \nabla \dot{u}_{\Omega}(\theta) \right) dxdt = - \int_Q j'(u_{\Omega}) \dot{u}_{\Omega}(\theta) dxdt. $$

Now introducing the variational characterization (4.19) of $\dot{u}_{\Omega}(\theta)$, we obtain: (4.21)

$$ \int_Q j'(u_{\Omega}) \dot{u}_{\Omega}(\theta) \, dxdt = \int_Q \sigma_{\Omega(t)} m'(0)(\theta) \frac{du_{\Omega}}{dt} p_{\Omega} \, dxdt - \int_Q \sigma_{\Omega(t)} [F'_{xx}(0)(\theta)] v \cdot \nabla u_{\Omega} p_{\Omega} \, dxdt \\ + \int_Q \sigma_{\Omega(t)} v_1(\theta) \cdot \nabla u_{\Omega} p_{\Omega} \, dxdt + \int_Q \sigma_{\Omega(t)} b'(0)(\theta) \cdot \nabla u_{\Omega} p_{\Omega} \, dxdt \\ + \int_Q \nu_{\Omega(t)} [A'(0)(\theta)] \nabla u_{\Omega} \cdot \nabla p_{\Omega} \, dxdt - \int_Q (m'(0)(\theta) f + f_1(\theta)) p_{\Omega} \, dxdt. $$

Injecting this expression into (4.20), we obtain the volume form (4.10) of the shape derivative of $J(\Omega)$.

Step 5: We infer the surface form (4.13) of $J'(\Omega)(\theta)$ from the volume form (4.10).

Let us write (4.10) under the form

$$ J'(\Omega)(\theta) = \sum_{i=1}^{4} I_i(\theta), $$

where $I_i(\theta)$, $i = 1, \dots, 4$ are the four integral quantities at the right-hand side of (4.10). Our derivation of the surface form of $J'(\Omega)(\theta)$ follows the methodology exposed in [2], which we have already exemplified in the proof of Lemma 3.1: we integrate by parts all the terms featured in the $I_i(\theta)$ which involve derivatives of $\theta$, and among the resulting terms, we only retain those featuring boundary integrals on $\Gamma$ of the normal component $\theta \cdot n$, since the other contributions are expected to cancel according to the prediction of the structure theorem for shape derivatives.

Let us treat the first integral $I_1(\theta)$; injecting the expression (3.13) of $m'(0)(\theta)$, we see that:

$$ I_1(\theta) = \int_Q (\text{div}\theta)(\varphi_t^{-1}(x)) j(u_{\Omega}(t, x)) \, dxdt + r(\theta). $$

Here and throughout the proof, $r(\theta)$ is a remainder that may change from one line to the next, of the form

$$ r(\theta) = \int_D s_{\Omega} \cdot \theta \, dx + \int_{\Gamma} t_{\Omega} \cdot \theta_{\Gamma} \, ds, $$

for some vector fields $s_{\Omega}: D \to \mathbb{R}^d$, $t_{\Omega}: \Gamma \to \mathbb{R}^d$ that depend on $u_{\Omega}$, $p_{\Omega}$ and the motion $\varphi_t$. Decomposing the integral in the above right-hand side as a nested integral over $(0, T)$ and $D$, and using a change of variables by $\varphi_t$ in the innermost integral, we obtain:

$$ \begin{aligned} I_1(\theta) &= \int_0^T \int_D |\det \nabla \varphi_t(x)| \text{div} \theta(x) j(u_{\Omega}(t, \varphi_t(x))) \, dxdt + r(\theta) \\ &= \int_0^T \int_{\partial D} |\det \nabla \varphi_t(x)| j(u_{\Omega}(t, \varphi_t(x))) (\theta \cdot n)(x) \, dxdt + r(\theta), \end{aligned} $$
where we have used an integration by parts from the first line to the second one, together with the continuity of $u_\Omega(t, \cdot)$ across $\Gamma(t)$, see Lemma 4.2. Since $\theta \in \Theta_{ad}$ vanishes on $\partial D$, it follows that:

$$I_1(\theta) = r(\theta).$$

Let us now turn to the second integral $I_2(\theta)$, which rewrites:

$$ \begin{aligned} I_2(\theta) &= \int_Q \sigma_{\Omega(t)} \left( m'(0)(\theta) \frac{du_{\Omega}}{dt} - [F'_{xx}(0)(\theta)] v \cdot \nabla u_{\Omega} + v_1(\theta) \cdot \nabla u_{\Omega} - F'_{xt}(0)(\theta) \cdot \nabla u_{\Omega} \right) p_{\Omega} \,dxdt \\ &= \int_0^T \int_D \sigma_{\Omega(t)}(x) \left( (\operatorname{div}\theta)(\varphi_t^{-1}(x)) \frac{du_{\Omega}}{dt}(t, x) - \left[\nabla\varphi_t(\varphi_t^{-1}(x))\right] \left[\nabla\theta(\varphi_t^{-1}(x))\right] \left[\nabla\varphi_t^{-1}(x)\right] v(t, x) \cdot \nabla u_{\Omega}(t, x) \right. \\ & \qquad \left. - \left[\nabla\varphi_t(\varphi_t^{-1}(x))\right] \left[\nabla\theta(\varphi_t^{-1}(x))\right] \frac{\partial}{\partial t}(\varphi_t^{-1}(x)) \cdot \nabla u_{\Omega} \right) p_{\Omega}(t, x) \,dxdt + r(\theta), \end{aligned} $$

where we have used the derivatives (3.10), (3.11) and (3.13). Again, the change of variables $x \mapsto \varphi_t(x)$ in the innermost integral yields:

$$ I_2(\theta) = \int_0^T \int_D \sigma_{\Omega}(x) \left| \det \nabla \varphi_t(x) \right| \left( (\operatorname{div}\theta)(x) \frac{du_{\Omega}}{dt}(t, \varphi_t(x)) - \left[\nabla \varphi_t(x)\right] \left[\nabla \theta(x)\right] \left[\nabla \varphi_t^{-1}(\varphi_t(x))\right] v(t, \varphi_t(x)) \cdot \nabla u_{\Omega}(t, \varphi_t(x)) \right. \\ \left. - \left[\nabla \varphi_t(x)\right] \left[\nabla \theta(x)\right] \frac{\partial \varphi_t^{-1}}{\partial t}(\varphi_t(x)) \cdot \nabla u_{\Omega}(t, \varphi_t(x)) \right) p_{\Omega}(t, \varphi_t(x)) \,dxdt + r(\theta), $$

and so

$$ I_2(\theta) = \int_0^T \int_D \sigma_{\Omega}(x) \left| \det \nabla \varphi_t(x) \right| \left( (\operatorname{div}\theta)(x) \frac{du_{\Omega}}{dt}(t, \varphi_t(x)) - \left[\nabla \theta(x)\right] \left[\nabla \varphi_t^{-1}(\varphi_t(x))\right] v(t, \varphi_t(x)) \cdot \left[\nabla \varphi_t(x)\right]^T \nabla u_{\Omega}(t, \varphi_t(x)) \right. \\ \left. - \left[\nabla \theta(x)\right] \frac{\partial \varphi_t^{-1}}{\partial t}(\varphi_t(x)) \cdot \left[\nabla \varphi_t(x)\right]^T \nabla u_{\Omega}(t, \varphi_t(x)) \right) p_{\Omega}(t, \varphi_t(x)) \,dxdt + r(\theta). $$

Decomposing the innermost integral as the sum of its contributions on $\Omega$ and $D \setminus \bar{\Omega}$, then integrating by parts with Proposition A.1, we obtain:

$$ \begin{aligned} I_2(\theta) = -(\sigma_2 - \sigma_1) \int_0^T \int_{\Gamma} \left| \det \nabla \varphi_t(x) \right| \frac{du_{\Omega}}{dt}(t, \varphi_t(x)) p_{\Omega}(t, \varphi_t(x)) (\theta \cdot n)(x) \,dsdt \\ + \int_0^T \int_{\Gamma} \left| \det \nabla \varphi_t(x) \right| \left( \left[ \nabla \varphi_t^{-1}(\varphi_t(x)) \right] v(t, \varphi_t(x)) \cdot n(x) + \frac{\partial \varphi_t^{-1}}{\partial t}(\varphi_t(x)) \cdot n(x) \right) \\ \left( \sigma_2 \nabla u_{\Omega}^{+}(t, \varphi_t(x)) - \sigma_1 \nabla u_{\Omega}^{-}(t, \varphi_t(x)) \right) \cdot \left( \left[ \nabla \varphi_t(x) \right] \theta(x) \right) p_{\Omega}(t, \varphi_t(x)) \,dsdt + r(\theta). \end{aligned} $$

Eventually, remarking that the definition (2.13) of the velocity $v(t, x)$ implies that

$$ \left[\nabla \varphi_t^{-1}(\varphi_t(x))\right] v(t, \varphi_t(x)) + \frac{\partial \varphi_t^{-1}}{\partial t}(\varphi_t(x)) = \frac{\partial}{\partial t}(\varphi_t^{-1}(\varphi_t(x))) = 0, $$

we obtain:

$$ I_2(\theta) = -(\sigma_2 - \sigma_1) \int_0^T \int_{\Gamma} \left| \det \nabla \varphi_t(x) \right| \frac{du_{\Omega}}{dt}(t, \varphi_t(x)) p_{\Omega}(t, \varphi_t(x)) (\theta \cdot n)(x) \,dsdt + r(\theta). $$
We next turn to the third integral $I_3(\theta)$ on the right-hand side of (4.10). At first, injecting the derivatives (3.11) and (3.17) into the defining formula, we have:

$$
\begin{align*}
I_3(\theta) &= \int_Q \nu_{\Omega(t)}[A'(0)(\theta)] \nabla u_{\Omega} \cdot \nabla p_{\Omega} \,dxdt \\
&= \int_0^T \int_D \nu_{\Omega(t)}(m'(0)(\theta)I - [F'_{xx}(0)(\theta)] - [F'_{xx}(0)(\theta)]^T) \nabla u_{\Omega} \cdot \nabla p_{\Omega} \,dxdt \\
&= \int_0^T \int_D \nu_{\Omega(t)}(x) \left( (\mathrm{div}\theta)(\varphi_t^{-1}(x))I - [\nabla\varphi_t(\varphi_t^{-1}(x))] [\nabla\theta(\varphi_t^{-1}(x))] [\nabla\varphi_t^{-1}(x)] \right. \\
&\qquad \left. - [\nabla\varphi_t^{-1}(x)]^T [\nabla\theta(\varphi_t^{-1}(x))]^T [\nabla\varphi_t(\varphi_t^{-1}(x))]^T \right) \nabla u_{\Omega}(t, x) \cdot \nabla p_{\Omega}(t, x) \,dxdt + r(\theta).
\end{align*}
$$

Hence, a change of variables based on the diffeomorphism $\varphi_t$ in the innermost integral produces:

$$
I_3(\theta) = \int_0^T \int_D \nu_\Omega(x) |\det \nabla \varphi_t(x)| \left( (\mathrm{div}\theta)(x) I - [\nabla \varphi_t(x)] [\nabla \theta(x)] [\nabla \varphi_t^{-1}(\varphi_t(x))] \right. \\
\left. - [\nabla \varphi_t^{-1}(\varphi_t(x))]^T [\nabla \theta(x)]^T [\nabla \varphi_t(x)]^T \right) \nabla u_\Omega(t, \varphi_t(x)) \cdot \nabla p_\Omega(t, \varphi_t(x)) \,dxdt + r(\theta),
$$

and so

$$
\begin{align*}
I_3(\theta) &= \int_0^T \int_D \nu_\Omega(x) |\det \nabla \varphi_t(x)| (\mathrm{div}\theta)(x) \nabla u_\Omega(t, \varphi_t(x)) \cdot \nabla p_\Omega(t, \varphi_t(x)) \,dxdt \\
&\quad - \int_0^T \int_D \nu_\Omega(x) |\det \nabla \varphi_t(x)| [\nabla \theta(x)] [\nabla \varphi_t^{-1}(\varphi_t(x))] \nabla u_\Omega(t, \varphi_t(x)) \cdot [\nabla \varphi_t(x)]^T \nabla p_\Omega(t, \varphi_t(x)) \,dxdt \\
&\quad - \int_0^T \int_D \nu_\Omega(x) |\det \nabla \varphi_t(x)| [\nabla \theta(x)] [\nabla \varphi_t^{-1}(\varphi_t(x))] \nabla p_\Omega(t, \varphi_t(x)) \cdot [\nabla \varphi_t(x)]^T \nabla u_\Omega(t, \varphi_t(x)) \,dxdt + r(\theta).
\end{align*}
$$

Decomposing each integral over $D$ as the sum of two integrals over $\Omega$ and $D \setminus \bar{\Omega}$, then using the integration by parts formulas in Proposition A.1, we obtain:

$$
\begin{align*}
I_3(\theta) = & - \int_0^T \int_{\Gamma} |\det \nabla \varphi_t(x)| (\nu_2 \nabla u_{\Omega}^+(t, \varphi_t(x)) \cdot \nabla p_{\Omega}^+(t, \varphi_t(x)) - \nu_1 \nabla u_{\Omega}^-(t, \varphi_t(x)) \cdot \nabla p_{\Omega}^-(t, \varphi_t(x))) (\theta \cdot n)(x) ds(x) dt \\
& + \int_0^T \int_{\Gamma} |\det \nabla \varphi_t(x)| \left( \nu_2 \left( [\nabla \varphi_t^{-1}(\varphi_t(x))] \nabla u_{\Omega}^+(t, \varphi_t(x)) \cdot n(x) \right) \left( [\nabla \varphi_t(x)]^T \nabla p_{\Omega}^+(t, \varphi_t(x)) \cdot \theta(x) \right) \right. \\
& \qquad \left. - \nu_1 \left( [\nabla \varphi_t^{-1}(\varphi_t(x))] \nabla u_{\Omega}^-(t, \varphi_t(x)) \cdot n(x) \right) \left( [\nabla \varphi_t(x)]^T \nabla p_{\Omega}^-(t, \varphi_t(x)) \cdot \theta(x) \right) \right) ds(x) dt \\
& + \int_0^T \int_{\Gamma} |\det \nabla \varphi_t(x)| \left( \nu_2 \left( [\nabla \varphi_t^{-1}(\varphi_t(x))] \nabla p_{\Omega}^+(t, \varphi_t(x)) \cdot n(x) \right) \left( [\nabla \varphi_t(x)]^T \nabla u_{\Omega}^+(t, \varphi_t(x)) \cdot \theta(x) \right) \right. \\
& \qquad \left. - \nu_1 \left( [\nabla \varphi_t^{-1}(\varphi_t(x))] \nabla p_{\Omega}^-(t, \varphi_t(x)) \cdot n(x) \right) \left( [\nabla \varphi_t(x)]^T \nabla u_{\Omega}^-(t, \varphi_t(x)) \cdot \theta(x) \right) \right) ds(x) dt + r(\theta),
\end{align*}
$$

and we now proceed to simplify each of the three integrals in the above right-hand side, which we denote by $J_1(\theta)$, $J_2(\theta)$ and $J_3(\theta)$, respectively.
At first, we introduce the tangential and normal components of the gradients of $u_\Omega(t, \cdot)$ and $p_\Omega(t, \cdot)$; using the jump relations from Lemma 4.2 and Remark 4.3, $J_1(\theta)$ rewrites:

$$J_1(\theta) = -(\nu_2 - \nu_1) \int_0^T \int_{\Gamma} |\det \nabla \varphi_t(x)| \nabla_{\Gamma(t)} u_\Omega(t, \varphi_t(x)) \cdot \nabla_{\Gamma(t)} p_\Omega(t, \varphi_t(x)) (\theta \cdot n)(x) \, ds(x)dt \\ - \left(\frac{1}{\nu_2} - \frac{1}{\nu_1}\right) \int_0^T \int_{\Gamma} |\det \nabla \varphi_t(x)| (\nu_{\Omega(t)} \nabla u_\Omega(t, \varphi_t(x)) \cdot n_{\Omega(t)}(\varphi_t(x))) (\nu_{\Omega(t)} \nabla p_\Omega(t, \varphi_t(x)) \cdot n_{\Omega(t)}(\varphi_t(x))) (\theta \cdot n)(x) \, ds(x)dt.$$

Let us next deal with $J_2(\theta)$; to this end, we observe that

$$|\det \nabla \varphi_t(x)| \nu_2 ( [\nabla \varphi_t^{-1}(\varphi_t(x))] \nabla u_{\Omega}^{+}(t, \varphi_t(x)) \cdot n(x) ) = \nu_2 ( \nabla u_{\Omega}^{+}(t, \varphi_t(x)) \cdot ( [\mathrm{com}(\nabla \varphi_t(x))] n(x) ) ),$$

where $\mathrm{com}(M)$ is the cofactor matrix of a $d \times d$ matrix. Hence, using Lemma A.1 and the jump relations from Lemma 4.2, it follows that

$$(4.22) \qquad \begin{aligned} |\det \nabla \varphi_t(x)| \nu_2 ( [\nabla \varphi_t^{-1}(\varphi_t(x))] \nabla u_{\Omega}^{+}(t, \varphi_t(x)) \cdot n(x) ) &= |\det \nabla \varphi_t(x)| \nu_1 ( [\nabla \varphi_t^{-1}(\varphi_t(x))] \nabla u_{\Omega}^{-}(t, \varphi_t(x)) \cdot n(x) ) \\ &= |\mathrm{com}(\nabla \varphi_t((x)))n(x)| (\nu_{\Omega(t)} \nabla u_{\Omega} \cdot n_{\Omega(t)}(t, \varphi_t(x))), \end{aligned}$$

and so:

$$(4.23) \quad J_2(\theta) = \int_0^T \int_{\Gamma} |\mathrm{com}(\nabla \varphi_t(x)) n(x)| (\nu_{\Omega} \nabla u_{\Omega} \cdot n_{\Omega(t)}(t, \varphi_t(x))) ( [\nabla \varphi_t(x)]^T \nabla p_{\Omega}^{+}(t, \varphi_t(x)) \cdot \theta(x) - [\nabla \varphi_t(x)]^T \nabla p_{\Omega}^{-}(t, \varphi_t(x)) \cdot \theta(x) ) ds(x)dt.$$

We now bring into play the jump relations satisfied by $p_\Omega$, summarized in Remark 4.3 to simplify the second parenthesis in the above integrand:

$$\begin{aligned} [\nabla \varphi_t(x)]^T \nabla p_{\Omega}^{+}(t, \varphi_t(x)) \cdot \theta(x) - [\nabla \varphi_t(x)]^T \nabla p_{\Omega}^{-}(t, \varphi_t(x)) \cdot \theta(x) & \\ &= (\nabla p_{\Omega}^{+}(t, \varphi_t(x)) - \nabla p_{\Omega}^{-}(t, \varphi_t(x))) \cdot ([\nabla \varphi_t(x)] \theta(x)) \\ &= \left(\frac{1}{\nu_2} - \frac{1}{\nu_1}\right) \left(\nu_{\Omega(t)} \nabla p_{\Omega}(t, \varphi_t(x)) \cdot n_{\Omega(t)}(\varphi_t(x))\right) ([\nabla \varphi_t(x)] \theta(x) \cdot n_{\Omega(t)}(\varphi_t(x))) \\ &= \left(\frac{1}{\nu_2} - \frac{1}{\nu_1}\right) \frac{|\det \nabla \varphi_t(x)|}{|\mathrm{com}(\nabla \varphi_t(x)) n(x)|} (\nu_{\Omega(t)} \nabla p_{\Omega}(t, \varphi_t(x)) \cdot n_{\Omega(t)}(\varphi_t(x))) (\theta \cdot n)(x). \end{aligned}$$

Hence, $J_2(\theta)$ equals:

$$J_2(\theta) = \left(\frac{1}{\nu_2} - \frac{1}{\nu_1}\right) \int_0^T \int_{\Gamma} |\det \nabla \varphi_t(x)| (\nu_{\Omega} \nabla u_{\Omega} \cdot n_{\Omega(t)}(t, \varphi_t(x))) (\nu_{\Omega} \nabla p_{\Omega} \cdot n_{\Omega(t)}(t, \varphi_t(x))) (\theta \cdot n)(x) \, ds(x)dt.$$

Applying a similar treatment to $J_3(\theta)$, it follows:

$$J_3(\theta) = -(\nu_2 - \nu_1) \int_0^T \int_{\Gamma} |\det \nabla \varphi_t(x)| \nabla_{\Gamma(t)} u_{\Omega}(t, \varphi_t(x)) \cdot \nabla_{\Gamma(t)} p_{\Omega}(t, \varphi_t(x)) (\theta \cdot n)(x) \, ds(x)dt \\ + \left(\frac{1}{\nu_2} - \frac{1}{\nu_1}\right) \int_0^T \int_{\Gamma} |\det \nabla \varphi_t(x)| (\nu_{\Omega} \nabla u_{\Omega}(t, \varphi_t(x)) \cdot n_{\Omega(t)}(\varphi_t(x))) (\nu_{\Omega} \nabla p_{\Omega}(t, \varphi_t(x)) \cdot n_{\Omega(t)}(\varphi_t(x))) (\theta \cdot n)(x) \, ds(x)dt \\ + r(\theta).$$

Eventually, the fourth and last integral $I_4(\theta)$ in the right-hand side of (4.21) equals:

$$I_4(\theta) = \int_Q (m'(0)(\theta)f + f_1(\theta)) p_{\Omega} \, dxdt = \int_0^T \int_D (\mathrm{div}\theta)(\varphi_t^{-1}(x)) f(t, x) p_{\Omega}(t, x) \, dxdt + r(\theta).$$

Since $p_\Omega(t, \cdot)$ is continuous across $\Gamma(t)$ (see again Remark 4.3), it follows like in the treatment of $I_1(\theta)$ that

$$I_4(\theta) = r(\theta).$$
Eventually, the desired surface expression (4.13) of the shape derivative $J'(\Omega)(\theta)$ is obtained by gathering the above expressions of the various integrals $I_i(\theta)$, $i = 1, \dots, 4$ and $J_j(\theta)$, $j = 1, \dots, 3$, and by checking that the integrals collected in the various remainders $r(\theta)$ do cancel, which is elementary, but tedious. This concludes the proof of the theorem. □

### 4.4. Extensions of these results to the non linear setting

In this section, we outline how the shape differentiability of the functional $J(\Omega)$ in (3.2), proved in Section 4.3, extends to the general situation where the function $\hat{v} : \mathbb{R}_+ \to \mathbb{R}_+$ characterizing the reluctivity of the ferromagnetic material is no longer a constant, in which case the state equation (2.21) for $u_\Omega$ becomes a non linear evolution problem.

We content ourselves with the statement of the volume form of the shape derivative of $J(\Omega)$, as it is the only information used in our numerical implementation, see Section 5.3. This formula results from a formal calculation similar to that involved in the proof of Proposition 4.1, see also [30] for a related analysis.

**Proposition 4.2.** The functional $J(\Omega)$ in (3.2), featuring the solution $u_\Omega$ to the non linear version of the evolution problem (2.21) with a reluctivity coefficient $\hat{v}$ as in (2.12), is shape differentiable at any bounded, Lipschitz shape $\Omega \subset D_{\text{rot}}$ and its shape derivative reads:

$$ (4.24) \quad \forall \theta \in \Theta_{\text{ad}}, \quad J'(\Omega)(\theta) = \int_Q m'(0)(\theta) j(u_\Omega) \,dx\,dt \\ \qquad + \int_Q \sigma_{\Omega(t)} \left( m'(0)(\theta) \frac{du_\Omega}{dt} - [F'_{xx}(0)(\theta)] v \cdot \nabla u_\Omega + v_1(\theta) \cdot \nabla u_\Omega + b'(0)(\theta) \cdot \nabla u_\Omega \right) p_\Omega \,dx\,dt \\ \qquad + \int_Q \left( \nu_{\Omega(t)}(x, |\nabla u_\Omega|) [A'(0)(\theta)] - \frac{\nu'_{\Omega(t)}(x, |\nabla u_\Omega|)}{|\nabla u_\Omega|} \left( [F'_{xx}(0)(\theta)]^T \nabla u_\Omega \cdot \nabla u_\Omega \right) I \right) \nabla u_\Omega \cdot \nabla p_\Omega \,dx\,dt \\ \qquad - \int_Q \left( m'(0)(\theta) f + f_1(\theta) \right) p_\Omega \,dx\,dt. $$

Here, $\nu'_{\Omega(t)}(x, s)$ stands for the derivative of the partial mapping $s \mapsto \nu_{\Omega(t)}(x, s)$, and $p_\Omega$ is the solution to the adjoint evolution problem:

$$ (4.25) \quad \text{Search for } p_{\Omega} \in W_{\text{per}} \text{ s.t. } \forall w \in L^2(0, T; H_0^1(D)), $$

$$ \int_Q \left( -\sigma_{\Omega(t)} \frac{\mathrm{d}p_{\Omega}}{\mathrm{d}t} w - \sigma_{\Omega(t)} (\mathrm{div} v) p_{\Omega} w + \left( \nu_{\Omega(t)}(x, |\nabla u_{\Omega}|) I + \frac{\nu'_{\Omega(t)}(x, |\nabla u_{\Omega}|)}{|\nabla u_{\Omega}|} \nabla u_{\Omega} \otimes \nabla u_{\Omega} \right) \nabla p_{\Omega} \cdot \nabla w \right) \mathrm{d}x\mathrm{d}t = - \int_Q j'(u_{\Omega}) w \mathrm{d}x\mathrm{d}t. $$

The quantities $b'(0)(\theta)$, $F'_{xx}(0)(\theta)$, $f_1(\theta)$, $v_1(\theta)$ and $A'(0)(\theta)$ are respectively given by (3.9), (3.11), (3.14), (3.15) and (3.17).

Again, note that the formula in (4.24) is supplemented by the terms (4.14) in the presence of permanent magnetization, see Remark 4.4.

**Remark 4.5.** As usual in the practice of the adjoint method, the defining problem for $p_\Omega$ is (the transpose of) the linearized version of the state problem (2.21) for $u_\Omega$. In particular, the well-posedness of (4.25) is not an immediate consequence of that of (2.21). In the present context, this property results from similar arguments to those in the proof of Theorem 4.1, once we have observed that the matrix field

$$ A(x) := \left( \hat{v}(|\nabla u_{\Omega}|) I + \frac{\hat{v}'(|\nabla u_{\Omega}|)}{|\nabla u_{\Omega}|} \nabla u_{\Omega} \otimes \nabla u_{\Omega} \right) $$

induces a linear and uniformly elliptic operator in $D_{\text{mag}}$. To see this, we first note that, by taking $s_2 = s$, $s_1 = s + h$ in (2.12) for $s > 0$ and sufficiently small $h \in \mathbb{R}$, then dividing both sides by $h^2$ and letting $h \to 0$, we obtain:

$$ (4.26) \qquad s\hat{v}'(s) + \hat{v}(s) \geq \underline{v}. $$
Now, for any unit vector $\xi \in \mathbb{R}^d$, the operator attached to $A(x)$ satisfies:

$$A(x)\xi \cdot \xi = \hat{\nu}(|\nabla u_{\Omega}|)|\xi|^2 + \frac{\hat{\nu}'(|\nabla u_{\Omega}|)}{|\nabla u_{\Omega}|} |\nabla u_{\Omega} \cdot \xi|^2.$$

Then, if $\hat{\nu}'(|\nabla u_{\Omega}|) \ge 0$, we obtain immediately that

$$A(x)\xi \cdot \xi \ge \nu|\xi|^2.$$

If on the contrary $\hat{\nu}'(|\nabla u_{\Omega}|) \le 0$ the Cauchy-Schwarz inequality together with (4.26) yield:

$$A(x)\xi \cdot \xi \ge \left( \hat{\nu}(|\nabla u_{\Omega}|) + \frac{\hat{\nu}'(|\nabla u_{\Omega}|)}{|\nabla u_{\Omega}|} |\nabla u_{\Omega}|^2 \right) |\xi|^2, \\ \ge \nu|\xi|^2,$$

which is the expected conclusion.

## 5. NUMERICAL METHODS

The shape optimization problems of interest in this work are of the form

$$(5.1) \qquad \min_{\Omega \subset D_{\text{rot}}} J(\Omega),$$

where $J$ depends on the 2d ferromagnetic phase $\Omega$ within (the cross-section of) the rotor $D_{\text{rot}}$ via the solution $u_{\Omega}$ to the (linear or non linear) magneto-quasi-static problem (2.21).

In this section, we describe the main numerical methods employed in our shape optimization framework. After sketching the overall strategy in Section 5.1, we outline in Section 5.2 the practical resolution of the state and adjoint evolution problems (2.21) and (4.11). We eventually discuss in Section 5.3 how the formulas provided by Propositions 4.1 and 4.2 are used in practice to identify a suitable descent direction for $J(\Omega)$.

### 5.1. Sketch of the shape optimization workflow

In a nutshell, we solve (5.1) thanks to a gradient descent algorithm, which improves the shape $\Omega$ by successive deformations along descent directions for the objective function $J(\Omega)$. The numerical discretization is based on a space-time finite element framework: the 3d space-time cylinder $Q = (0, T) \times D$ is consistently equipped with a tetrahedral mesh $\mathcal{K}$ which (approximately) resolves the moving phases $\Omega(t), \Omega_a(t), D_{\text{mag}}(t)$ etc. for all times $t \in [0, T]$. In our experiments, this space-time mesh is constructed by defining the moving space-time geometry, e.g., using the Open Cascade Technology (OCCT) geometry kernel$^2$, and creating an unstructured tetrahedral mesh for this geometry using the mesh generator Netgen/NGSolve$^3$.

The global strategy is outlined in Algorithm 1. In there, and throughout the following, we denote with an $n$ superscript the instances of the different objects at stake (the shape $\Omega$, the descent direction $\theta$, ...) at each iteration $n = 0, \dots$ of the process.

**Algorithm 1** Shape gradient algorithm for the solution of (5.1).

**Initialization:** 3d tetrahedral mesh $\mathcal{K}^0$ of $Q$ where the moving phases $\varphi_t(\Omega^0), \varphi_t(\Omega_a^0)$, etc. are discretized for $t \in (0, T)$.

**for** $n = 0, \dots$, until convergence **do**

(1) Solve the state (2.21) and adjoint (4.25) equations for $u_{\Omega^n}, p_{\Omega^n} : Q \to \mathbb{R}$.

(2) Evaluate the volume form of the shape derivative $J'(\Omega^n)(\theta)$ (4.24).

(3) Identify a descent direction $\theta^n : D \to \mathbb{R}^d$ for $J(\Omega)$ from $\Omega^n$ using the Hilbertian procedure (5.2).

(4) Select a suitably small descent step $\tau^n$ and calculate the space-time deformation $\Theta^n$ associated to $\tau^n \theta^n$ via (3.6).

(5) Modify the space-time mesh $\mathcal{K}^n$ according to $\Theta^n$ to obtain the new mesh $\mathcal{K}^{n+1}$.

**end for**

**return** Mesh $\mathcal{K}^n$ of $Q$, where the deformations $\Omega^n(t^k)$ of the optimized 2d shape $\Omega^n$ are explicitly discretized.

<footnote>$^2$https://dev.opencascade.org/doc/refman/html/index.html</footnote>

<footnote>$^3$https://ngsolve.org/</footnote>
Briefly, the algorithm starts with the datum of a tetrahedral mesh of the 3d space-time cylinder where the configurations of the various moving or fixed phases (the ferromagnetic phase $\phi_t(\Omega^0)$, the air phase $\phi_t(\Omega_a^0)$, etc.) are approximately discretized for $t \in [0, T]$.

At each iteration $n = 0, \dots$ of the process, the evolution problems for $u_{\Omega^n}$ and $p_{\Omega^n}$ are solved on the mesh $\mathcal{K}^n$ of the space-time cylinder via the space-time finite element method. This task is described in Section 5.2. This allows to evaluate the volume form of the shape derivative $J'(\Omega^n)$ and to extract a suitable descent direction $\theta^n$ from the latter via the so-called Hilbertian procedure. This operation is detailed in Section 5.3. Then, the space-time deformation mapping $\Theta^n$ induced by $\tau^n \theta^n$ is calculated, where $\tau^n$ is a “small” descent step, and it is used to deform the mesh $\mathcal{K}^n$ of $Q$ into a new mesh $\mathcal{K}^{n+1}$ where the updated shape $\Omega^{n+1} = (\mathrm{Id} + \tau^n \theta^n)(\Omega^n)$ (and the other constituent phases of the motor) is explicitly discretized.

## 5.2. Resolution of the state and adjoint equations

At any iteration $n$ of the execution of Algorithm 1, the space-time cylinder $Q = (0, T) \times D$ is equipped with a tetrahedral mesh $\mathcal{K}^n$ in which the deformed structures of the motor (in particular, the ferromagnetic phase $\Omega^n(t)$, the air phase $\Omega_a^n(t)$, etc.) for $t \in [0, T]$ are explicitly discretized. Hence, the motion of the structure of the motor induced by $\phi_t$ is directly encoded in the space-time mesh.

This mesh serves as the support for the numerical solution of the state and adjoint equations (2.21) and (4.25) for $u_{\Omega^n}$ and $p_{\Omega^n}$ via a space-time finite element method. These are discretized on the mesh $\mathcal{K}^n$ by means of continuous, piecewise linear finite elements as suggested in [68], see also [27] for an application to the parallel simulation of an electric machine. We mention that also higher order polynomial degrees could be used, see e.g. [47] for a related investigation in an optimal control setting. The non linearity in the magneto-quasi-static problem (2.21) for the state $u_{\Omega^n}$ is handled thanks to a Newton-Raphson method with damping in order to ensure a decrease of the residual. The time-periodic conditions can be treated in a straightforward way by identifying the degrees of freedom lying on the bottom and top sections of the mesh $\mathcal{K}^n$ of the space-time cylinder $Q$ provided that the bottom and top surfaces of the space-time cylinder feature an identical mesh.

Our numerical experiments rely on the finite element environment NGSolve [64]. Note that the precise setting at hand does not allow to take advantage of the capabilities of this package to compute automatically shape derivatives in the context of more “standard” shape optimization problems [31].

### Remark 5.1.

* As we have mentioned in the introduction, in principle, space-time finite element approaches allow for parallelization of the numerical resolution of an evolution problem in space-time by means of classical domain decomposition approaches, see e.g. [27]. In the same spirit, adaptive mesh refinement in space-time is possible [47]. Last but not least, the state and adjoint equations could actually be gathered into a single system and be solved in parallel, allowing for another level of parallelization when compared to more conventional time-stepping schemes. While we do solve the associated systems of linear equations in parallel, as discussed in Section 6.2.2, adaptivity and the combination of state and adjoint are not exploited in the present study, but they are interesting leads for future work.

* The present use of a space-time finite element method is possible since the physical problem under scrutiny arises in two space dimensions, thus resulting in a 3d space-time mesh. This strategy could be extended to the case of three space dimensions, provided numerical method for meshing (3 + 1)-dimensional space-time domains and solving related variational problems be available [47].

## 5.3. Evaluation of the shape derivative $J'(\Omega)(\theta)$ and calculation of a descent direction

At any iteration $n = 0, \dots$ of Algorithm 1, the reference to which is omitted in this section for notational simplicity, we rely on the volume form of the shape derivative $J'(\Omega)(\theta)$ to identify a descent direction for $J(\Omega)$. As discussed in e.g. [2, 32], the use of the volume form is a priori less convenient to achieve this purpose than that of the surface form. Nevertheless, it presents several advantages, such as the less stringent regularity demanded from $\Omega$, $u_\Omega$, $p_\Omega$, a better consistency between the theoretical and discrete shape optimization settings [37], and the possibility to conduct part of the evaluation in an automated way [31].

To exploit the volume form of $J'(\Omega)(\theta)$, we rely on the so-called “Hilbertian method”, for which we refer to e.g. [7, 18, 22]. Let $V$ be a Hilbert space of vector fields on $D$, equipped with an inner product $b(\cdot, \cdot)$,
which is continuously embedded in $\Theta_{ad}$. We solve the following auxiliary variational problem for the shape gradient $\theta$ of the derivative $J'(\Omega)$ associated to the Hilbert space $(V, b(\cdot, \cdot))$:

$$ (5.2) \qquad \text{Search for } \theta \in V \text{ s.t. for all } \eta \in V, \quad b(\theta, \eta) = J'(\Omega)(\eta), $$

If $\theta$ is the solution to (5.3), then $-\theta$ is a descent direction for $J(\Omega)$, since

$$ J'(\Omega)(-\theta) = -b(\theta, \theta) < 0. $$

Moreover, $\theta$ benefits from the features encoded in the space $V$, such as for instance a higher regularity than mere vector fields in $\Theta_{ad}$, etc.

Several choices are possible about $V$ and $b$, starting from Sobolev spaces $H^m(D)^d$, for $m$ large enough. One efficient and popular choice, albeit formal, is the following:

$$ (5.3) \qquad V = H_0^1(D)^d, \text{ with the inner product } b(\theta, \eta) := \int_D (\alpha(x)(\nabla\theta + \nabla\theta^\top) : \nabla\eta + \beta(x)\theta \cdot \eta) \,dx, $$

where $\alpha, \beta: D \to \mathbb{R}$ are smooth, positive and non negative weights, respectively. Intuitively, this choice, inspired by the linear elasticity system, produces shape gradient which “resemble” elastic displacements, and thereby cause as little compression as possible in the mesh. In the same spirit, one may also enrich $b(\cdot, \cdot)$ with a term attached to the Cauchy-Riemann equations, as a means to urge the solution $\theta$ to (5.2) to be (approximately) conformal, thus preserving (approximately) the angles of the mesh, see [38].

The numerical solution of the identification problem (5.2) raises the need to calculate the volume form of $J'(\Omega)(\theta)$ – by this, we mean, for a given shape $\Omega$, to calculate the quantity $J'(\Omega)(\theta)$ for a wide range of deformation fields $\theta$. To achieve this, we note that the shape derivatives (3.19), (4.10) and (4.24) of the functionals considered in this article can be written under the form

$$ J'(\Omega)(\theta) = \int_0^T \int_{\varphi_t(D)} \left( S_0(t, x) \cdot \theta(\varphi_t^{-1}(x)) + S_1(t, x) : \nabla\theta(\varphi_t^{-1}(x)) \right) dx dt, $$

for some vector and matrix fields $S_0: Q \to \mathbb{R}^d$ and $S_1: Q \to \mathbb{R}^{d \times d}$, see [48] for a general result in this spirit. Hence, a change of variables yields:

$$ \begin{aligned} J'(\Omega)(\theta) &= \int_0^T \int_D |\det\nabla\varphi_t(x)| (S_0(t, \varphi_t(x)) \cdot \theta(x) + S_1(t, \varphi_t(x)) : \nabla\theta(x)) \,dx \,dt \\ &= \int_D g_0(x) \cdot \theta(x) + g_1(x) : \nabla\theta(x) \,dx, \end{aligned} $$

with

$$ (5.4) \qquad g_0(x) = \int_0^T S_0(t, \varphi_t(x)) |\det \nabla \varphi_t(x)| \, dt, \quad \text{and} \quad g_1(x) = \int_0^T S_1(t, \varphi_t(x)) |\det \nabla \varphi_t(x)| \, dt. $$

Hence, in practice, when it comes to evaluate $J'(\Omega)(\theta)$, we first calculate the fields $g_0(x)$ and $g_1(x)$ on a triangular mesh $\mathcal{T}$ of the spatial domain $D$. For instance, $\mathcal{T}$ can be a mesh of the bottom surface of the space-time cylinder $Q$ enclosed in the total 3d mesh $\mathcal{K}$. We use a piecewise constant approximation of the functions $g_i$ ($i = 0, 1$) on $\mathcal{T}$: the constant value of $g_i$ on each triangle $T \in \mathcal{T}$ is obtained by evaluating the “vertical line integral” featured in (5.4) at the centroid of $T$; this is realized, in turn, by traveling the space-time mesh $\mathcal{K}$ of $Q$ and using a composite trapezoidal rule.

Eventually, once a descent direction $\theta$ is obtained for $J(\Omega)$ and a suitable descent step $\tau > 0$ is selected, the corresponding space-time deformation mapping

$$ \Theta(t, x) = (t, \varphi_t \circ (\mathrm{Id} + \tau\theta) \circ \varphi_t^{-1}(x)) $$

is used to update the space-time mesh $\mathcal{K}$, see again Fig. 3.

## 6. NUMERICAL EXAMPLES

In this section we discuss two numerical examples to illustrate the shape derivative formulas of Proposi- tions 4.1 and 4.2 and the implementation of our numerical strategy dedicated to shape optimization prob- lems involving moving domains. The first Section 6.1 deals with a rather academic situation where the
one-dimensional optimized shape is subjected to a motion which is slightly more complicated than a rota-
tion. We then turn in the Section 6.2 to the original motivation for this work, namely the optimization of
the 2d structure of a rotating electric machine.

6.1. An academic, spatially one-dimensional example

This first example arises in one space dimension ($d = 1$). The domain $D$ is the unit interval $(0, 1) \subset \mathbb{R}$ and we consider a simplified version of the structure described in Section 2.3: $D$ is solely made of two complementary phases $\Omega^1 := \Omega$ and $\Omega^2 := D \setminus \bar{\Omega}$, whose constituent materials have constant conductivity and reluctivity, set to $\sigma_1 = \nu_1 = 1$ and $\sigma_2 = \nu_2 = 5$, respectively.

This arrangement of $D$ is subjected to the deformation $\varphi_t(x) = x+tx^2$ over the time period $(0, T) = (0, 1)$.
The physical behavior of $\Omega$ is then accounted for by the solution $u_\Omega$ to the linear version of the magneto-
quasi-static equation (2.21), where the source term $f(t, x)$ reads:

$$
f(t, \varphi_t(x)) = (x - 0.4)(x - 0.6)\sqrt{\varphi_t(x)}(1 + t - \varphi_t(x))
$$

see Fig. 4 (a).

Admittedly, this setting does not exactly comply with the assumptions of Sections 2, 3.2 and 4: the mapping $t \mapsto \varphi_t(x)$ is not periodic, and so the domains $\Omega$, $D$ at initial time do not coincide with their respective final versions $\Omega(T)$, $D(T)$. However, a simple verification reveals that all the findings in there – and notably, the well-posedness Theorem 4.1 about (2.21) and the calculation of the shape derivative in Proposition 4.1 – remain valid in this context, provided the time-periodicity condition encompassed in $W_{\text{per}}$ is replaced with the more general version:

$$
u(0, x) = u(T, \varphi_T(x)), \quad x \in D.
$$

In this setting, we consider the shape optimization problem

$$
(6.1) \qquad \min_{\Omega \subset D} J(\Omega), \text{ where } J(\Omega) = \int_Q u_{\Omega}(t, x) \,dx\,dt
$$

where $u_\Omega$ solves (2.21) with $\sigma_1 = 10$, $\sigma_2 = 0$, $\nu_1 = 1$, $\nu_2 = 10$, recalling the notation of Section 4.3. This problem (6.1) is solved by means of the shape gradient Algorithm 1. The inner product $b(\cdot, \cdot)$ featured in the practice (5.3) of the Hilbertian method is induced by the Laplace operator: we solve (5.2) with $\alpha = 0.5$, $\beta = 0$.

The initial shape $\Omega^0$ is the interval $(0.4, 0.6)$, see Fig. 4 (b), and the solution $u_{\Omega^0}$ to (2.21) is depicted
in Fig. 4 (c). The deformation $\Theta^0$ applied to the space-time mesh at the first iteration is represented on
Fig. 5 (a). We apply 59 iterations of the gradient-based Algorithm 1 until the norm of the calculated
descent direction gets lower than a prescribed threshold of $10^{-9}$. The optimized design $\Omega^*$ as well as the
corresponding potential $u_{\Omega^*}$ are depicted in Fig. 5 (b) and (c), respectively. The value of $J(\Omega)$ decreases
from $5.091 \cdot 10^{-4}$ to $4.231 \cdot 10^{-4}$ in the process.

<figure><img src="image_5.png" /><figcaption>FIGURE 4. (a) Right-hand side $f(t,x)$; (b) Space-time cylinder $Q$ associated to the initial design $\Omega^0$ in the 1d academic example of Section 6.1; (c) Potential $u_{\Omega^0}$ associated to $\Omega^0$.</figcaption></figure>

6.2. Optimization of the structure of a rotating electric machine

In this section, we consider the optimization of the structure of an electric motor.
<figure><img src="image_6.png" /><figcaption>FIGURE 5. (a) Descent direction $\Theta^0$ from the initial design $\Omega^0$, after extension of $\theta^0$ to the space-time cylinder; (b) Space-time cylinder $Q_{\Omega^*}$ associated to the optimized design $\Omega^*$; (c) Potential $u_{\Omega^*}$ for the optimized design $\Omega^*$.</figcaption></figure>

### 6.2.1. Description of the geometric setting

The considered machine operates at steady state with a constant rotational speed of 600 rotations per minute. The geometry of its 2d cross-section $D$ is symmetric with respect to both coordinates axes, see again Fig. 1; hence, the time horizon $T$ is chosen as the duration of a 90 degree rotation, i.e. $T = 0.025$ s.

The conductivities of the various materials at play take the values:

$$ \sigma_f = \sigma_a = \sigma_c = 0, \quad \sigma_m = 10^6 \text{ A}^2 \cdot \text{s}^3 \cdot \text{kg}^{-1} \cdot \text{m}^{-3}, $$

and the reluctivity of air, copper and magnet equal:

$$ \nu_a = \nu_c = \frac{10^7}{4\pi} \text{ Am} \cdot \text{V}^{-1} \cdot \text{s}^{-1}, \quad \text{and} \quad \nu_m = \frac{\nu_a}{1.05} \text{ Am} \cdot \text{V}^{-1} \cdot \text{s}^{-1}. $$

The material behaviour of the ferromagnetic material is characterized by the nonlinear reluctivity function $\hat{\nu}$ defined by:

$$ \hat{\nu}(s) = \nu_a - (\nu_a - c_1)\exp(-c_2s^{c_3}), $$

with the parameters $c_1 = 200$, $c_2 = 0.001$, $c_3 = 6$, see Fig. 6 (a). We also assume that the magnet region $D_{\text{mag}}$ is the support of a permanent magnetization field $M$. The direction of $M$ is depicted in Fig. 6 (b) and its magnitude equals $|M| = \nu_m B_r A \cdot m^{-1}$, where $B_r = 1.216 \text{ V} \cdot \text{s} \cdot \text{m}^{-2}$ is the magnetic remanence.

The coil region $D_{\text{stat},c}$ of the stator is occupied by a large number of thin copper wires wound around the iron core of the machine. This phase is subdivided into six regions $\Omega_{U+}$, $\Omega_{U-}$, $\Omega_{V+}$, $\Omega_{V-}$, $\Omega_{W+}$, $\Omega_{W-}$ and the source current density $f$ is spatially homogeneous with smooth time variations in each of these:

$$ (6.2) \quad f(t, x) = (1_{\Omega_{U+}}(x) - 1_{\Omega_{U-}}(x)) I_U(t) + (1_{\Omega_{V+}}(x) - 1_{\Omega_{V-}}(x)) I_V(t) + (1_{\Omega_{W+}}(x) - 1_{\Omega_{W-}}(x)) I_W(t). $$

In this expression, the current densities read:

$$ I_U(t) = a \sin(4\alpha(t) + \psi), \quad I_V(t) = a \sin\left(4\alpha(t) + \psi + \frac{4}{3}\pi\right), \quad I_W(t) = a \sin\left(4\alpha(t) + \psi + \frac{2}{3}\pi\right), $$

where $\alpha(t) = \frac{\pi}{2} \frac{t}{T}$ and $\psi = \frac{\pi}{18}$ is called the current angle. Here, the amplitude $a$ is given by $a = I c / A$ where $I = 1555.64 \text{ A}$ is the product of the current intensity and the number of copper wires per coil, $c = 2.75$ is the stacking factor and $A \approx 1.8 \cdot 10^{-4} \text{ m}^2$ is the area of one coil.

In this application, the transformation $\varphi_t$ stands for a rotation of the inner part $D_{\text{rot}}$ at constant angular velocity $\alpha'(t) = \frac{\pi}{2T} \text{ s}^{-1}$, i.e.,

$$ \varphi_t(x) = \begin{cases} R_{\alpha(t)}x & \text{if } x \in D_{\text{rot}} \\ x & \text{if } x \in D_{\text{stat}}, \end{cases} \quad \text{where } R_{\alpha(t)} = \begin{pmatrix} \cos \alpha(t) & -\sin \alpha(t) \\ \sin \alpha(t) & \cos \alpha(t) \end{pmatrix}. $$

In this setting, we aim to optimize the average torque of the device by acting on the shape of the phase $\Omega \subset D_{\text{rot}}$ containing ferromagnetic material – and thereby on the shape of the air pockets near the magnets in the rotor domain – while keeping fixed the region $D_{\text{mag}}$ containing the permanent magnets:

$$ (6.3) \qquad \min_{\Omega \subset D} J(\Omega), \quad \text{where } J(\Omega) := -\text{Tor}(u_{\Omega}), $$

where the magneto-quasi-static potential $u_{\Omega}$ is the solution to the non linear evolution problem (2.21) and the torque $\text{Tor}(u)$ is defined in (2.9). Note that, strictly speaking, $J(\Omega)$ is not of the form considered
<figure><img src="image_7.png" /><figcaption>FIGURE 6. Magnetic reluctivity curve $s \mapsto \hat{v}(s)$ in the ferromagnetic material considered in the motor optimal design example of Section 6.2; (b) Zoom on the initial design $\Omega^0$.</figcaption></figure>

in Section 4, as it involves the gradient of the potential $u_\Omega$. Nevertheless, as noted in Section 3.1, our calculations performed can be adapted straightforwardly to handle this case.

### 6.2.2. Practical implementation details

At each iteration $n = 0, \dots, \text{Algorithm 1}$, the tetrahedral space-time mesh $\mathcal{K}^n$ of the cylinder $Q$ consists of 320,597 vertices and 1,758,666 elements. As discussed in Section 5.2, the non linear magneto-quasi-static problem (2.21) for $u_\Omega$ is solved by a Newton-Raphson method with damping, based on the space-time finite element framework. In the first optimization iteration $n = 0$, this method is initialized with the function 0, while at each subsequent iteration $n \ge 1$, the numerical solution $u_{\Omega^{n-1}}$ to (2.21) at the previous iteration is chosen instead, which reduces the number of necessary Newton iterations. The systems of linear equations involved in the procedure are solved in parallel with openMPI [25] and MUMPS [4], via the PETSc[9] library. We use NGSolve for both the meshing and the assembly of the finite element stiffness matrices and right-hand side vectors. Since, at the current state, NGSolve does not support MPI parallelization in combination with periodic boundary conditions, the meshing and finite element assembly operations can be carried out on one core only. We run the simulations on RICAM's RADON cluster<fnref n="4" /> by using a single node with 40 cores and 1TB of RAM. The solution of the state problem (2.21) at the first optimization iteration $n = 0$ requires 19 Newton iterations, for about one hour of computation. The needed CPU time to achieve this operation is reduced by half at the subsequent steps, as only 9 Newton iterations are needed. The (linear) adjoint problem (4.11) is solved in the same way within about two minutes.

The evaluation of the volume form of the shape derivative $J'(\Omega)(\theta)$ is carried out along the lines of Section 5.3. Since only the interface $\Gamma$ between the ferromagnetic material and air phases $\Omega$ and $\Omega_a$ is optimized, we restrict the resulting descent directions $\theta$ considered in the shape derivative (4.24) and the auxiliary problem (5.3) to the rotor domain $D_{\text{rot}}$ deprived of the magnet region $D_{\text{mag}}$, i.e. we actually solve the identification problem (5.2) with the Hilbert space $V = H_0^1(\Omega \cup \Omega_a)^2$ and the inner product in (5.3) with the weights $\alpha(x) \equiv 1$ and $\beta(x) \equiv 0$. In addition, following [38] we add a term of Cauchy-Riemann type to the bilinear form $b(\cdot, \cdot)$ which helps in preserving the mesh quality.

**Remark 6.1.** Note that the smoothness assumption about the source function $f$ in the derivation of Section 4 is not satisfied by (6.2). However, as noted in Remark 4.4, Propositions 4.1 and 4.2 still hold in the present case where $f(t, \cdot)$ has support only in $D_{\text{stat}}$.

<footnote n="4">⁴https://www oeaw ac at ricam hpc</footnote>
<figure><img src="image_8.png" /><figcaption>FIGURE 7. Solution $u_{\Omega^0}$ to the version of the magneto-quasi-static state problem (2.21) considered in Section 6.2, attached to the initial design $\Omega^0$ at times (a) $t = 0$; (b) $t = T/2$; (c) $t = T$.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>FIGURE 8. (a) Spatial descent direction $\theta^0$ at the initial iteration; (b) Corresponding space-time deformation $(0, \Theta^0(t, x) - x)$.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>FIGURE 9. Optimized design of the rotating machine considered in Section 6.2 at times (a) $t = 0$; (b) $t = T/2$, and (c) $t = T$.</figcaption></figure>
### 6.2.3. Description of the numerical results

The initial design $\Omega^0$ is sketched in Fig. 6 (b), and the values of the potential $u_{\Omega^0}(t, \cdot)$ at times $t = 0$, $t = T/2$ and $t = T$ are depicted in Fig. 7 in three cross sections of the space-time mesh. The first descent direction $\theta^0: D \to \mathbb{R}^d$ is shown in Fig. 8 (a) and the corresponding space-time vector field $(0, \Theta_x^0(t, x) - x)$ is represented at several times in Fig. 8 (b).

At each iteration of Algorithm 1, a suitable descent step $\tau^n$ is calculated via a line search procedure in which $\tau^n$ is decreased by half as long as the scaled deformation $\tau^n\theta^n$ does not allow for a decrease of the objective function $J(\Omega)$. This procedure is interrupted when $\tau^n$ becomes smaller than $10^{-10}$.

The optimized design resulting from this procedure after 14 iterations is depicted in Fig. 9 at the three different times $t = 0$, $t = T/2$ and $t = T$; the torque $\text{Tor}(u_\Omega)$ is increased from 522.94 N·m to 587.79 N·m in the process. We note that, at the final iteration, the norm of the vector field $\theta^n$ does not vanish, but it no longer allows to decrease the value of $J(\Omega)$. One possible explanation for this fact is that the components of $\Omega_a$ made of the air pockets close to the magnet regions $D_{\text{mag}}$ tend to get closer and eventually merge with the air gap region $D_{\text{gap}}$, which is not allowed by the homogeneous Dirichlet boundary condition on $\partial D_{\text{stat}}$ imposed in the identification problem (5.2) for $\theta$. In principle, this undesirable behavior could be prevented by adding a constraint about the mechanical stiffness of the device to the optimization problem (6.3), as in e.g. [16]. This issue will be addressed in a further work.

**Remark 6.2.** A careful look at Fig. 9 shows that the rotated versions of the optimized design at the considered times $t^k$ do not have exactly the same look. This numerical artifact is caused by the approximation of the extension of the deformation vector field $\theta$ into the space-time deformation mapping $\Theta$ in (3.6), which is realized on an unstructured mesh of $Q$.

## 7. CONCLUSION AND PERSPECTIVES

In this article, we have addressed a shape optimization problem related to the internal structure of an electric motor. This work departs from the related literature by the fact that the physical problem at stake is described by a version of the non linear magneto-quasi-static equation featuring a time evolving geometry – a rich setting, which captures fine, realistic effects. After proving the well-posedness of the physical problem and calculating the shape derivative of a generic optimization criterion in this context, we have proposed a shape gradient workflow based on the space-time finite element method. We have validated our numerical strategy on an academic problem and tackled a more realistic example in the physical context of interest.

The work in this article paves the way to multiple further investigations. On the one hand, the physical model, although already complicated enough, could be enriched to include multiphysics effects, for instance couplings between electromagnetic and thermal effects. This could require time horizons larger than one single rotation period, and could be addressed by combining the present space-time finite element method with time-stepping techniques. On a different note, constraints about the mechanical integrity and stability of the optimized design could be modeled with the help of structural mechanics, as suggested in [73]. From the numerical viewpoint, more robust methods could be implemented to deal with the update of the shape, such as the level set method [3] or one of its avatars [1, 20] retaining an exact, meshed description of the optimized phases. Eventually, optimal design strategies based on topological derivatives could also be considered in this context, see e.g. [29] for a related investigation.

**Acknowledgments.** The work of A. C. and P. G. is supported by the FWF funded project P32911 as well as the joint DFG/FWF Collaborative Research Centre CREATOR (CRC - TRR361/F90) at TU Darmstadt, TU Graz, RICAM and JKU Linz. This work was completed while C. D. was visiting the Laboratoire Jacques-Louis Lions from Université Paris-Sorbonne, whose hospitality is thankfully acknowledged.

## APPENDIX A. A FEW USEFUL FORMULAS IN SHAPE CALCULUS

The following integration by parts formula is a version of the usual Green's formula that is used repeatedly in the main part of the article.

**Proposition A.1.** Let $\Omega \subset \mathbb{R}^d$ be a Lipschitz bounded domain.
(i) Let $\theta \in W^{1,\infty}(\mathbb{R}^d, \mathbb{R}^d)$ and $a: \mathbb{R}^d \to \mathbb{R}$ be a smooth scalar field. Then:

$$ \int_{\Omega} \frac{\partial \theta_i}{\partial x_j} a \, dx = \int_{\partial \Omega} \theta_i a n_j \, ds - \int_{\Omega} \frac{\partial a}{\partial x_j} \theta_i \, dx. $$

(ii) Let $\theta \in W^{1,\infty}(\mathbb{R}^d, \mathbb{R}^d)$ and $a, b: \mathbb{R}^d \to \mathbb{R}^d$ be smooth vector fields. The following integration by parts formula holds:

$$ \int_{\Omega} \nabla \theta a \cdot b \, dx = \int_{\partial \Omega} (\theta \cdot b)(a \cdot n) \, ds - \int_{\Omega} (\mathrm{div} a)(\theta \cdot b) \, dx - \int_{\Omega} \nabla b a \cdot \theta \, dx $$

Let us also recall the following relation between the normal and tangent vectors to the boundary of a domain and those attached to a deformed configuration.

**Lemma A.1.** Let $\Omega \subset \mathbb{R}^d$ be a bounded Lipschitz domain, and let $\varphi: \mathbb{R}^d \to \mathbb{R}^d$ be a Lipschitz continuous homeomorphism with Lipschitz inverse. Then, for any tangential vector field $\tau: \partial\Omega \to \mathbb{R}^d$, the vector field $\nabla\varphi(\varphi^{-1}(y))\tau(\varphi^{-1}(y))$ is tangential to $\partial(\varphi(\Omega))$. Besides, the unit normal vector fields $n_\Omega$ and $n_{\varphi(\Omega)}$ are related by:

$$ \forall x \in \Omega, \quad n_{\varphi(\Omega)}(\varphi(x)) = \frac{1}{|\mathrm{com}(\nabla\varphi(x))n_{\Omega}(x)|} \mathrm{com}(\nabla\varphi(x))n_{\Omega}(x), $$

where $\mathrm{com}(M) \in \mathbb{R}^{d \times d}$ is the cofactor matrix of a $d \times d$ matrix $M$.

# APPENDIX B. MAIN RESULTS ABOUT SPACE-TIME VARIATIONAL PROBLEMS

This appendix is devoted to the proof of **Theorem 4.1**, which states the well-posedness of the version (2.21) of the non linear magneto-quasi-static equation featuring coefficients depending on the rotating distribution of materials within the domain $D$. Although the proof hinges on quite classical methods for the analysis of non linear evolution equations and mixed parabolic-elliptic equations, the exact setting considered in this article, featuring a time-dependent geometry, is new to the best of our knowledge (see however [8, 41] for related investigations).

After recalling a few facts about time-space functional spaces in **Appendix B.1** and providing the statement of some needed results about the well-posedness of non linear boundary value and evolution problem in **Appendix B.2**, we detail the proof of **Theorem 4.1**, properly speaking, in **Appendix B.3**.

## B.1. A few words about time-space functional spaces

This section gathers a few definitions and basic facts about the time-space functional spaces which are naturally involved in the mathematical formulation of evolution problems.

Let $V \subset H \subset V^*$ be an evolution triple, that is:

* $(V, \langle \cdot, \cdot \rangle_V)$ and $(H, \langle \cdot, \cdot \rangle_H)$ are real and separable Hilbert spaces, and $V^*$ is the topological dual of $V$; we denote by $\langle \cdot, \cdot \rangle_{V^*,V}$ the usual duality pairing between $V^*$ and $V$.

* The inclusion $V \subset H$ is continuous and dense.

Given a fixed time interval $[0, T]$, we define the space $L^2(0, T; V)$ by:

$$ (B.1) \qquad L^2(0, T; V) = \left\{ u : [0, T] \to V, \int_0^T \|u(t)\|_V^2 \, dt < \infty \right\}. $$

The latter is a Hilbert space when equipped with the natural inner product

$$ \langle u(t), v(t) \rangle_{L^2(0,T;V)} := \int_0^T \langle u(t), v(t) \rangle_V \, dt. $$

Its topological dual $L^2(0, T; V)^*$ is naturally identified to $L^2(0, T; V^*)$ via the pairing:

$$ (B.2) \qquad \langle f, u \rangle_{L^2(0,T;V^*),L^2(0,T;V)} = \int_0^T \langle f(t), u(t) \rangle_{V^*,V} \, dt. $$
**Definition B.1.** Let $V \subset H \subset V^*$ be an evolution triple, and let $Z$ be another Hilbert space. A function $u \in L^2(0, T; V)$ has a time derivative in the space $L^2(0, T; Z)$ if there exists $z \in L^2(0, T; Z)$ such that, for all real-valued test functions $\varphi \in C_c^\infty(0, T)$, it holds:

$$ \int_0^T u(t)\varphi'(t) \, dt = - \int_0^T z(t)\varphi(t) \, dt. $$

This derivative is then denoted by $\frac{\partial u}{\partial t} = z$.

Let $V \subset H \subset V^*$ be an evolution triple. The time-dependent problems considered in this article bring into play Hilbert spaces of the form

$$ (B.3) \qquad W(0, T; V, H) = \{u \in L^2(0, T; V) \text{ s.t. } \frac{\partial u}{\partial t} \in L^2(0, T; V^*)\}, $$

equipped with the natural norm

$$ \|u\|_{W(0,T;V,H)}^2 := \int_0^T \|u(t)\|_V^2 \, dt + \int_0^T \left\| \frac{\partial u}{\partial t}(t) \right\|_{V^*}^2 \, dt. $$

We recall the following facts about such functional spaces, see Lemmas 1 and 2 of Chap XVIII, §1.2 of [21].

### Lemma B.1.

(i) The subspace

$$ \left\{ u(t) = \sum_{k=0}^{K} t^k u_k, \ K \ge 0, \ u_k \in V \right\} \subset W(0, T; V, H) $$

of polynomial functions in the time variable with coefficients in $V$ is dense in $W(0, T; V, H)$. In particular, the space $C^\infty([0, T], V)$ of smooth $V$-valued functions on $[0, T]$ is dense in $W(0, T; V, H)$.

(ii) The space $W(0, T; V, H)$ is continuously embedded in the space $C([0, T]; H)$ of continuous $H$-valued functions on $[0, T]$.

The last point of this lemma is of crucial importance in practice, as it gives a rigorous meaning to the time (e.g. initial, periodic) conditions usually imposed to the solution of an evolution problem when it is sought in $W(0, T; V, H)$.

## B.2. Preliminary results about the well-posedness of some non linear problems

Let us first recall the following consequence of the Banach fixed point theorem, about the solution to a non linear stationary equation in a Hilbert space, whose operator satisfies suitable monotonicity and continuity conditions, see Th. 25 B. from [75].

**Theorem B.1.** Let $(V, \|\cdot\|_V)$ be a real Hilbert space and let $\mathcal{A}: V \to V^*$ be a (possibly non linear) mapping satisfying the following properties:

(i) $\mathcal{A}$ is strongly monotone, i.e. there exists $\alpha > 0$ such that:

$$ \forall u, v \in V, \quad \langle \mathcal{A}u - \mathcal{A}v, u - v \rangle_{V^*, V} \ge \alpha \|u - v\|_V^2. $$

(ii) $\mathcal{A}$ is Lipschitz continuous, i.e. there exists $L > 0$ such that:

$$ \forall u, v \in V, \quad \|\mathcal{A}u - \mathcal{A}v\|_{V^*} \le L \|u - v\|_V. $$

Then, for each $b \in V^*$, the equation

$$ (B.4) \qquad \mathcal{A}u = b $$

has a unique solution $u \in V$. This solution has a Lipschitz continuous dependence with respect to the data: there exists a constant $C > 0$ depending only on $\alpha$ and $L$ such that, if $u_1, u_2$ are the solutions to (B.4) with respective right-hand sides $b_1, b_2 \in V^*$, it holds:

$$ \|u_1 - u_2\|_V \le C \|b_1 - b_2\|_{V^*}. $$
In particular, the above result is the pivotal ingredient in the proof of the following lemma about the well-posedness of a stationary version of the considered evolution problem (2.21), see [59] for a similar argument:

**Lemma B.2.** Let $\mathcal{O}$ be a bounded Lipschitz domain in $\mathbb{R}^d$, let $\hat{\nu} : \mathbb{R}_+ \to \mathbb{R}_+$ be a function satisfying (2.12), and let $B \in C^\infty(\bar{\mathcal{O}}; \mathbb{R}^{d \times d})$ be a matrix-valued function such that $B(x)$ is invertible for all $x \in \bar{\mathcal{O}}$. Then for all data $f \in L^2(\mathcal{O})$ and $g \in H^{1/2}(\partial\mathcal{O})$, the boundary value problem

$$ (B.5) \qquad \left\{ \begin{array}{ll} -\mathrm{div}(\hat{\nu}(|B^T(x)\nabla u|)B(x)B^T(x)\nabla u) = f & \text{in } \mathcal{O}, \\ u = g & \text{on } \partial\mathcal{O}, \end{array} \right. $$

has a unique solution $u \in H^1(\mathcal{O})$ which is Lipschitz continuous with respect to the data: the solutions $u_1, u_2 \in H^1(\mathcal{O})$ associated to different data $(f_1, g_1), (f_2, g_2) \in L^2(\mathcal{O}) \times H^{1/2}(\partial\mathcal{O})$ satisfy the following inequality:

$$ \|u_1 - u_2\|_{H^1(\mathcal{O})} \le C (\|f_1 - f_2\|_{L^2(\mathcal{O})} + \|g_1 - g_2\|_{H^{1/2}(\partial\mathcal{O})}), $$

for a constant $C$ which depends only on the domain $\mathcal{O}$, the matrix $B$ and the coefficients $\underline{\nu}, \bar{\nu}$ in (2.12).

*Proof.* Let $\tilde{g} \in H^1(\mathcal{O})$ be such that $\tilde{g}|_{\partial\mathcal{O}} = g$. Letting the change of unknown functions $w = u - \tilde{g} \in H_0^1(\mathcal{O})$, the considered boundary value problem (B.5) rewrites:

$$ (B.6) \qquad \left\{ \begin{array}{ll} -\mathrm{div}(\hat{\nu}(|B^T(x)(\nabla w + \nabla \tilde{g})|)B(x)B^T(x)(\nabla w + \nabla \tilde{g})) = f & \text{in } \mathcal{O}, \\ w = 0 & \text{on } \partial\mathcal{O}. \end{array} \right. $$

Let us equip $H_0^1(\mathcal{O})$ with the inner product

$$ \langle u, v \rangle_{H_0^1(\mathcal{O})} := \int_{\mathcal{O}} \nabla u \cdot \nabla v \, dx, $$

and let us introduce the operator $\mathcal{A}: H_0^1(\mathcal{O}) \to H^{-1}(\mathcal{O})$ defined by:

$$ \mathcal{A}w = -\mathrm{div}(\hat{\nu}(|B^T(x)(\nabla w + \nabla \tilde{g})|)B(x)B^T(x)(\nabla w + \nabla \tilde{g})). $$

We aim to show that $\mathcal{A}$ is strongly monotone; to this end, let $u_1, u_2 \in H_0^1(\mathcal{O})$ be arbitrary functions; introducing the shorthand $U_i = B^T(x)(\nabla u_i + \nabla \tilde{g}) \in L^2(\mathcal{O})^d$, $i = 1, 2$, an elementary calculation yields:

$$ \begin{align*} \langle \mathcal{A}u_1 - \mathcal{A}u_2, u_1 - u_2 \rangle_{H^{-1}(\mathcal{O}), H_0^1(\mathcal{O})} & \\ &= \int_{\mathcal{O}} (\hat{\nu}(|U_1|)U_1 - \hat{\nu}(|U_2|)U_2) \cdot (U_1 - U_2) \, dx \\ &= \underline{\nu} \int_{\mathcal{O}} |U_1 - U_2|^2 \, dx + \int_{\mathcal{O}} ((\hat{\nu}(|U_1|) - \underline{\nu})U_1 - (\hat{\nu}(|U_2|) - \underline{\nu})U_2) \cdot (U_1 - U_2) \, dx \\ &= \underline{\nu} \int_{\mathcal{O}} |U_1 - U_2|^2 \, dx \\ &\quad + \int_{\mathcal{O}} ((\hat{\nu}(|U_1|) - \underline{\nu})|U_1|^2 + (\hat{\nu}(|U_2|) - \underline{\nu})|U_2|^2 - (\hat{\nu}(|U_1|) + \hat{\nu}(|U_2|) - 2\underline{\nu})U_1 \cdot U_2) \, dx \\ &\ge \underline{\nu} \int_{\mathcal{O}} |U_1 - U_2|^2 \, dx \\ &\quad + \int_{\mathcal{O}} ((\hat{\nu}(|U_1|) - \underline{\nu})|U_1|^2 + (\hat{\nu}(|U_2|) - \underline{\nu})|U_2|^2 - (\hat{\nu}(|U_1|) + \hat{\nu}(|U_2|) - 2\underline{\nu})|U_1||U_2|) \, dx \\ &= \underline{\nu} \int_{\mathcal{O}} |U_1 - U_2|^2 \, dx + \int_{\mathcal{O}} (\hat{\nu}(|U_1|)|U_1| - \hat{\nu}(|U_2|)|U_2|) (|U_1| - |U_2|) \, dx \\ &\quad - \underline{\nu} \int_{\mathcal{O}} ||U_1| - |U_2||^2 \, dx \\ &\ge \underline{\nu} \|U_1 - U_2\|_{L^2(\mathcal{O})}^2, \end{align*} $$
where we have used the Cauchy-Schwarz inequality to pass from the third to the fourth line and (2.12) to obtain the last line. Now, by definition of $U_1, U_2$ it holds:

$$||U_1 - U_2||^2_{L^2(\mathcal{O})^d} = \int_{\mathcal{O}} B(x)B^T(x)(\nabla u_1 - \nabla u_2) \cdot (\nabla u_1 - \nabla u_2) \, dx.$$

Besides, since $B$ is a smooth matrix-valued function on $\bar{\mathcal{O}}$ such that $B(x)$ is invertible for all $x \in \bar{\mathcal{O}}$, there exists a constant $\underline{\gamma} > 0$ such that

$$\forall x \in \bar{\mathcal{O}}, \forall \xi \in \mathbb{R}^d, \quad B(x)B^T(x)\xi \cdot \xi > \underline{\gamma}|\xi|^2.$$

We eventually infer from these facts that:

$$\langle Au_1 - Au_2, u_1 - u_2 \rangle_{H^{-1}(\mathcal{O}), H_0^1(\mathcal{O})} \geq \underline{\gamma} ||\nabla u_1 - \nabla u_2||^2_{L^2(\mathcal{O})^d},$$

which precisely expresses that strong monotonicity of $\mathcal{A}$.

A similar calculation shows that $\mathcal{A}: H_0^1(\mathcal{O}) \to H^{-1}(\mathcal{O})$ is a Lipschitz continuous operator, with a Lipschitz constant depending only on $\bar{\nu}$ and the matrix $B(x)$.

It follows from Theorem B.1 that the problem (B.6) has a unique solution $w$, with Lipschitz continuous dependence on the data function $f$, which readily implies the desired statement about (B.5). $\square$

Our study of (2.21) rests on results about the well-posedness of non linear evolution problems which stem from the general theory of maximal monotone operators. For simplicity, we limit ourselves with a relatively simple statement which is enough for our purpose, see e.g. Chap. 32 of [75] (and notably Th. 32.D).

**Theorem B.2.** Let $V \subset H \subset V^*$ be an evolution triple. For every $t \in (0, T)$, let $\mathcal{A}(t) : V \to V^*$ be a strongly monotone and Lipschitz continuous operator, whose monotonicity and Lipschitz constants do not depend on $t$. Then for each $f \in L^2(0, T; V^*)$, the evolution equation

$$\left\{ \begin{array}{ll} \frac{\partial u}{\partial t} + \mathcal{A}u = f & \text{in } (0, T), \\ u(0) = u(T), & \end{array} \right. $$

has a unique solution $u$ in $W(0, T; V, H)$.

**Remark B.1.** Strictly speaking, the statement of Th. 32.D in [75] features a non linear operator $\mathcal{A}$ which is independent of time. However, inspection of the proof of this result reveals that it carries over mutatis mutandis to the time-dependent case, provided the strong monotonicity and Lipschitz constants of $\mathcal{A}$ are independent of time.

### B.3. Proof of Theorem 4.1

We proceed in four steps. In the first three steps, we assume that there exists a solution $u \in W_{\text{per}}$ to the variational problem (4.2)-(4.3), and we construct an equivalent problem for $u$, see (B.9); the well-posedness of the latter results from the application of Theorem B.2. This procedure implies, in particular, that a solution $u$ to (4.2) is unique, if it exists. In the final step, we construct a solution to (4.2) from the solution to (B.9), thus concluding about the existence part of the statement.

**Step 1:** We use a change of variables to obtain an equivalent variational problem to (4.3) featuring a fixed arrangement of phases within $D$.

Let us assume that (4.2)-(4.3) has a solution $u$ in the space $W_{\text{per}}$ defined by (4.1). Then, for a.e. $t \in (0, T)$, the following identity is satisfied for all test functions $w \in H_0^1(D)$:

$$(B.7) \quad \int_D \sigma_{\Omega(t)}(x) \left( \frac{\partial u}{\partial t}(t, x) + v(t, x) \cdot \nabla u(t, x) \right) w(x) \, dx \\ \qquad + \int_D \nu_{\Omega(t)}(x, |\nabla u(t, x)|) \nabla u(t, x) \cdot \nabla w(x) \, dx = \int_D f(t, x) w(x) \, dx.$$

We introduce the transported function

$$\bar{u} \in L^2(0, T; H_0^1(D)), \quad \bar{u}(t, x) := u(t, \varphi_t(x)) \text{ for } t \in [0, T] \text{ and } x \in D.$$
By considering test functions of the form $\bar{w} = w \circ \varphi_t^{-1}$, $w \in H_0^1(D)$, a change of variables based on $\varphi_t$ in (B.7) yields, for a.e. $t \in (0, T)$:

$$ (B.8) \quad \forall w \in H_0^1(D), \quad \int_D \sigma_\Omega(x) \frac{\partial \bar{u}}{\partial t}(t,x)w(x) \,dx \\ \qquad + \int_D \nu_\Omega(x, |B_t^T(x)\nabla \bar{u}(t,x)|)B_t(x)B_t^T(x)\nabla \bar{u}(t,x) \cdot \nabla w(x) \,dx = \int_D f(t, \varphi_t(x))w(x) \,dx, $$

where, we have defined the smooth matrix-valued function:

$$ B_t(x) = \nabla \varphi_t^{-1}(x), $$

and we have used the calculus identities (2.17). The evolution problem (B.8) is of mixed parabolic-elliptic problem, since $\sigma_\Omega$ vanishes outside the subset $D_{\text{mag}}$ of $D$, where it takes the value $\sigma_m > 0$.

Step 2: We reformulate (B.8) as a parabolic evolution problem posed on $D_{\text{mag}}$.

By using test functions $w$ with compact support inside $D \setminus \overline{D_{\text{mag}}}$ in (B.8), we see that $\bar{u}$ satisfies:

$$ -\text{div}(\nu_{\Omega}(x, |B_t^T(x)\nabla \bar{u}(t,x)|)B_t(x)B_t^T(x)\nabla \bar{u}(t,x)) = f(t, \varphi_t(x)) \quad \text{in } D \setminus \overline{D_{\text{mag}}}. $$

On the other hand, an integration by parts in (B.8) reveals that the spatial restriction of $\bar{u}(t, \cdot)$ to $D_{\text{mag}}$, which we still denote by $\bar{u}$, satisfies the following problem:

$$ (B.9) \quad \forall w \in H^1(D_{\text{mag}}), \quad \sigma_m \int_{D_{\text{mag}}} \frac{\partial \bar{u}}{\partial t}(t,x)w(x) \,dx + \langle A\bar{u}, w \rangle_{H^1(D_{\text{mag}})^*, H^1(D_{\text{mag}})} \\ \qquad + \langle B\bar{u}, w \rangle_{H^1(D_{\text{mag}})^*, H^1(D_{\text{mag}})} = \int_{D_{\text{mag}}} f(t, \varphi_t(x))w(x) \,dx. $$

Here, we have introduced the linear operator $A: H^1(D_{\text{mag}}) \to H^1(D_{\text{mag}})^*$ and the non linear operator $B: H^1(D_{\text{mag}}) \to H^1(D_{\text{mag}})^*$ defined by, respectively:

$$ \forall u, w \in H^1(D_{\text{mag}}), \quad \langle Au, w \rangle_{H^1(D_{\text{mag}})^*, H^1(D_{\text{mag}})} = \nu_m \int_{D_{\text{mag}}} B_t(x)B_t^T(x)\nabla u(t,x) \cdot \nabla w(x) \,dx, $$

and

$$ \forall u, v \in H^1(D_{\text{mag}}), \quad \langle Bu, w \rangle_{H^1(D_{\text{mag}})^*, H^1(D_{\text{mag}})} = \\ - \int_{\partial D_{\text{mag}}} \left( \nu_{\Omega}(x, |B_t^T(x)\nabla \mathcal{L}_t u(t,x)|)B_t(x)B_t^T(x)\nabla \mathcal{L}_t u(t,x) \right) \cdot n(x) w(x) \,ds(x), $$

where $n$ is the unit normal vector to $\partial D_{\text{mag}}$ pointing outward $D_{\text{mag}}$ and for $z \in H^1(D_{\text{mag}})$, $\mathcal{L}_t z \in H^1(D \setminus \overline{D_{\text{mag}}})$ is the unique solution to the boundary value problem:

$$ (B.10) \quad \left\{ \begin{array}{ll} -\text{div}(\nu(x, |B_t^T(x)\nabla \mathcal{L}_t z(x)|)B_t(x)B_t^T(x)\nabla \mathcal{L}_t z(x)) = f(t,x) & \text{in } D \setminus \overline{D_{\text{mag}}}, \\ \mathcal{L}_t z = z & \text{on } \partial D_{\text{mag}}, \\ \mathcal{L}_t z = 0 & \text{on } \partial D. \end{array} \right. $$

Let us then analyze $A$ and $B$.

* The linear operator $A$ is continuous and monotone; more precisely, it holds:

$$ \forall u \in H^1(D_{\text{mag}}), \quad \langle Au, u \rangle_{H^1(D_{\text{mag}})^*, H^1(D_{\text{mag}})} \geq C \|\nabla u\|_{L^2(D_{\text{mag}})}^2. $$

Here and throughout the rest of the proof, $C$ stands for a positive constant which may change from one line to the other but is independent of time. The existence of such a constant in the above inequality is guaranteed by the smoothness and invertibility of the mapping $\varphi$.

* The operator $B$ is Lipschitz continuous, with a Lipschitz constant independent of $t$, as readily follows from the application of Lemma B.2 to the non linear operator $\mathcal{L}_t: H^1(D_{\text{mag}}) \to H^1(D \setminus \overline{D_{\text{mag}}})$.

Moreover, for all elements $u_1, u_2 \in H^1(D_{\text{mag}})$, several integrations by parts followed by a calculation identical to that conducted in the proof of Lemma B.2 yield:

$$ \langle Bu_1 - Bu_2, u_1 - u_2 \rangle_{H^1(D_{\text{mag}})^*, H^1(D_{\text{mag}})} \geq C \|\nabla \mathcal{L}_t u_1 - \nabla \mathcal{L}_t u_2\|_{L^2(D \setminus \overline{D_{\text{mag}}})}^2. $$
Since the mappings $u \mapsto ||\nabla u||_{L^2(D)^d}$ and $u \mapsto ||u||_{H^1(D \setminus \overline{D_{\text{mag}}})}$ are equivalent norms on the subspace of $H^1(D \setminus \overline{D_{\text{mag}}})$ functions with null trace on $\partial D$, we obtain:

$$ \langle \mathcal{B}u_1 - \mathcal{B}u_2, u_1 - u_2 \rangle_{H^1(D_{\text{mag}})^*, H^1(D_{\text{mag}})} \geq C ||u_1 - u_2||_{H^1(D \setminus \overline{D_{\text{mag}}})}^2 \\ \geq C ||u_1 - u_2||_{H^{1/2}(\partial D_{\text{mag}})}^2, $$

where the second line follows from the trace inequality.

The combination of both points reveals that the sum $A + B$ is Lipschitz continuous; it is also strongly monotone since, for all functions $u_1, u_2 \in H^1(D_{\text{mag}})$, it holds:

$$ \langle (A + B)u_1 - (A + B)u_2, u_1 - u_2 \rangle_{H^1(D_{\text{mag}})^*, H^1(D_{\text{mag}})} \geq C (||\nabla(u_1 - u_2)||_{L^2(D_{\text{mag}})^d}^2 + ||u_1 - u_2||_{H^{1/2}(\partial D_{\text{mag}})}^2) \\ \geq C ||u_1 - u_2||_{H^1(D_{\text{mag}})}^2, $$

as follows from (an avatar of) Poincaré's inequality. Moreover, the Lipschitz and strong monotonicity constants of $(A + B)$ are independent of time.

### Step 3: We use the abstract Theorem B.2.

The application of Theorem B.2 reveals that the variational problem (B.9) has a unique solution $\bar{u} \in W(0, T; H^1(D_{\text{mag}}), L^2(D_{\text{mag}}))$ such that $\bar{u}(t = 0, \cdot) = \bar{u}(t = T, \cdot)$. This almost readily implies that the solution $u \in W_{\text{per}}$ to the original evolution problem (4.2) is unique, when it exists. Indeed, if $u_i, i = 1, 2$ are two such solutions, let us define the functions $\bar{u}_i \in L^2(0, T; H_0^1(D))$ and $z_i \in L^2(0, T; H^1(D_{\text{mag}}))$ by:

$$ \bar{u}_i(t, x) = u_i(t, \varphi_t(x)), \text{ and } z_i(t, \cdot) = \bar{u}_i(t, \cdot)|_{D_{\text{mag}}}. $$

As a result of the previous steps, both functions $z_i$ belong to $W(0, T; H^1(D_{\text{mag}}), L^2(D_{\text{mag}}))$ and satisfy the evolution equation (B.9) with the time periodic condition $z_i(t = 0, \cdot) = z_i(t = T, \cdot)$, $i = 1, 2$. As a result, $z_1 = z_2$, i.e. for a.e. $t \in [0, T]$, $\bar{u}_1(t, \cdot)$ and $\bar{u}_2(t, \cdot)$ coincide on $D_{\text{mag}}$. Moreover, since $\bar{u}_i(t, \cdot)|_{D \setminus \overline{D_{\text{mag}}}} = \mathcal{L}_t z_i$ on $D \setminus \overline{D_{\text{mag}}}$, it immediately follows that $\bar{u}_1 = \bar{u}_2$ in $L^2(0, T; H_0^1(D))$, so that $u_1 = u_2$.

Summarizing, we have proved that there exists at most one solution $u \in W_{\text{per}}$ to (B.9).

### Step 4: We construct a solution $u \in W_{\text{per}}$ to (B.9).

This task essentially relies on the previous three steps. Let $z \in W(0, T; H^1(D_{\text{mag}}), L^2(D_{\text{mag}}))$ be the unique solution to the well-posed variational problem (B.9) equipped with time periodic boundary conditions:

$$ z(t = 0, \cdot) = z(t = T, \cdot) \text{ in } L^2(D_{\text{mag}}). $$

In particular, $z(t, \cdot)$ is an $H^1(D_{\text{mag}})$ function for a.e. $t \in [0, T]$, and so we may introduce the unique solution $z_{\text{ext}}(t, \cdot) = \mathcal{L}_t z(t, \cdot) \in H^1(D \setminus \overline{D_{\text{mag}}})$ to the version of (B.10) featuring Dirichlet data $z(t, \cdot)$ on $\partial D_{\text{mag}}$. Let now $\bar{u}(t, \cdot) \in H_0^1(D)$ be defined by:

$$ \bar{u}(t, x) = \begin{cases} z(t, x) & \text{if } x \in D_{\text{mag}}, \\ z_{\text{ext}}(t, x) & \text{if } x \in D \setminus \overline{D_{\text{mag}}}. \end{cases} $$

By construction, $\bar{u}$ belongs to $L^2(0, T; H_0^1(D))$, and it is one solution to (B.8).

It is also easily seen that $\sigma_{\Omega} \frac{\partial \bar{u}}{\partial t} \in L^2(0, T; H^{-1}(D))$, since $\frac{\partial z}{\partial t} \in L^2(0, T; H^1(D_{\text{mag}})^*)$ and for a.e. $t \in [0, T]$ and for any test function $\psi \in H_0^1(D)$,

$$ \left\langle \sigma_{\Omega} \frac{\partial \bar{u}}{\partial t}(t, \cdot), \psi \right\rangle_{H^{-1}(D), H_0^1(D)} = \sigma_m \left\langle \frac{\partial z}{\partial t}(t, \cdot), \psi \right\rangle_{H^1(D_{\text{mag}})^*, H^1(D_{\text{mag}})}. $$

Eventually, let $u \in W_{\text{per}}$ be defined by

$$ u(t, x) = \bar{u}(t, \varphi_t^{-1}(x)), \quad t \in [0, T], \quad x \in D. $$

By reusing the calculations conducted in the first step, we see that $u$ is one solution to (B.9). The Lipschitz dependence of this function on the data function $f$ results immediately from the Lipschitz dependence of the solutions to the problems (B.9) and (B.10) with respect to their right-hand sides, and we omit the details for brevity. This terminates the proof. $\square$
REFERENCES

[1] G. ALLAIRE, C. DAPOGNY, AND P. FREY, *Shape optimization with a level set based mesh evolution method*, Computer Methods in Applied Mechanics and Engineering, 282 (2014), pp. 22–53.

[2] G. ALLAIRE, C. DAPOGNY, AND F. JOUVE, *Shape and topology optimization*, in Geometric partial differential equations, part II, A. Bonito and R. Nochetto eds., Handbook of Numerical Analysis, vol. 22, (2021), pp. 1–132.

[3] G. ALLAIRE, F. JOUVE, AND A.-M. TOADER, *Structural optimization using sensitivity analysis and a level-set method*, Journal of computational physics, 194 (2004), pp. 363–393.

[4] P. AMESTOY, I. DUFF, J.-Y. L'EXCELLENT, AND J. KOSTER, *Mumps: A general purpose distributed memory sparse solver*, 01 2000, pp. 121–130.

[5] S. AMSTUTZ, *An introduction to the topological derivative*, Engineering Computations, 39 (2022), pp. 3–33.

[6] S. AMSTUTZ AND H. ANDRÄ, *A new algorithm for topology optimization using a level-set method*, Journal of Computational Physics, 216 (2006), pp. 573–588.

[7] H. AZEGAMI AND Z. C. WU, *Domain optimization analysis in linear elastic problems: approach using traction method*, JSME international journal. Ser. A, Mechanics and material engineering, 39 (1996), pp. 272–278.

[8] F. BACHINGER, U. LANGER, AND J. SCHÖBERL, *Numerical analysis of nonlinear multiharmonic eddy current problems*, Numerische Mathematik, 100 (2005), pp. 593–616.

[9] S. BALAY, S. ABHYANKAR, M. F. ADAMS, S. BENSON, J. BROWN, P. BRUNE, K. BUSCHELMAN, E. CONSTANTINESCU, L. DALCIN, A. DENER, V. EİJKHOUT, J. FAIBUSSOWITSCH, W. D. GROPP, V. HAPLA, T. ISAAC, P. JOLIVET, D. KARPEEV, D. KAUSHIK, M. G. KNEPLEY, F. KONG, S. KRUGER, D. A. MAY, L. C. MCINNES, R. T. MILLS, L. MITCHELL, T. MUNSON, J. E. ROMAN, K. RUPP, P. SANAN, J. SARICH, B. F. SMITH, S. ZAMPINI, H. ZHANG, H. ZHANG, AND J. ZHANG, *PETSc/TAO users manual*, Tech. Rep. ANL-21/39 - Revision 3.20, Argonne National Laboratory, 2023.

[10] M. P. BENDSOE AND O. SIGMUND, *Topology optimization: theory, methods, and applications*, Springer Science & Business Media, 2013.

[11] S. BLAUTH, *cashocs: A computational, adjoint-based shape optimization and optimal control software*, SoftwareX, 13 (2021), p. 100646.

[12] G. BRAMERDORFER, J. A. TAPIA, J. J. PYRHÖNEN, AND A. CAVAGNINO, *Modern electrical machine design optimization: Techniques, trends, and best practices*, IEEE Transactions on Industrial Electronics, 65 (2018), pp. 7672–7684.

[13] R. BRÜGGER AND H. HARBRECHT, *On the reformulation of the classical stefan problem as a shape optimization problem*, SIAM Journal on Control and Optimization, 60 (2022), pp. 310–329.

[14] R. BRÜGGER, H. HARBRECHT, AND J. TAUSCH, *On the numerical solution of a time-dependent shape optimization problem for the heat equation*, SIAM Journal on Control and Optimization, 59 (2021), pp. 931–953.

[15] O. BRUN, O. CHADEBEC, P. FERROUILLAT, I. NIYONZIMA, Z. LUO, Y. LE FLOCH, J. SIAU, F. VI, AND L. GERBAUD, *A level-set based topology optimization for maximizing the torque of switch reluctance machines*, in 24th International Conference on the Computation of Electromagnetic Fields (Compumag'2023)), 2023.

[16] O. BRUN, O. CHADEBEC, P. FERROUILLAT, I. NIYONZIMA, J. SIAU, L. GERBAUD, F. VI, AND Y. LE FLOCH, *Magneto-mechanical free shape optimization of synchronous reluctance machines*, in 17th International Workshop on Optimization and Inverse Problems in Electromagnetism (OIPE'2023), 2023.

[17] A. BUFFA, H. AMMARI, AND J.-C. NÉDÉLEC, *A justification of eddy currents model for the maxwell equations*, SIAM Journal on Applied Mathematics, 60 (2000), pp. 1805–1823.

[18] M. BURGER, *A framework for the construction of level set methods for shape optimization and reconstruction*, Interfaces and Free boundaries, 5 (2003), pp. 301–329.

[19] B. D. CULLITY AND C. D. GRAHAM, *Introduction to magnetic materials*, John Wiley & Sons, 2011.

[20] C. DAPOGNY AND F. FEPPON, *Shape optimization using a level set based mesh evolution method: an overview and tutorial*, Comptes Rendus. Mathématique, 361 (2023), pp. 1267–1332.

[21] R. DAUTRAY AND J.-L. LIONS, *Evolution problems i, volume 5 of mathematical analysis and numerical methods for science and technology*, 1992.

[22] F. DE GOURNAY, *Velocity extension for the level-set method and multiple eigenvalues in shape optimization*, SIAM journal on control and optimization, 45 (2006), pp. 343–367.

[23] M. C. DELFOUR AND J.-P. ZOLÉSIO, *Shapes and geometries: metrics, analysis, differential calculus, and optimization*, SIAM, 2011.

[24] A. E. FITZGERALD, C. KINGSLEY, AND S. D. UMANS, *Electric machinery*, McGraw-Hill, 2003.

[25] E. GABRIEL, G. E. FAGG, G. BOSILCA, T. ANGSKUN, J. J. DONGARRA, J. M. SQUYRES, V. SAHAY, P. KAMBADUR, B. BARRETT, A. LUMSDAINE, R. H. CASTAIN, D. J. DANIEL, R. L. GRAHAM, AND T. S. WOODALL, *Open MPI: Goals, concept, and design of a next generation MPI implementation*, in Proceedings, 11th European PVM/MPI Users' Group Meeting, Budapest, Hungary, September 2004, pp. 97–104.

[26] P. GANGL, *Sensitivity-based topology and shape optimization with application to electrical machines*, PhD thesis, Ph. D. thesis, Johannes Kepler University Linz, 2016.

[27] P. GANGL, M. GOBRIAL, AND O. STEINBACH, *A parallel space-time finite element method for the simulation of an electric motor*, Tech. Rep. 35, RICAM, 2022.

[28] P. GANGL, M. GOBRIAL, AND O. STEINBACH, *A space-time finite element method for the eddy current approximation of rotating electric machines*, 2023. https://arxiv.org/abs/2307.00278.

[29] P. GANGL AND N. KRENN, *Topology optimization of a rotating electric machine by the topological derivative*, PAMM, 23 (2023), p. e202200052.
[30] P. GANGL, U. LANGER, A. LAURAIN, H. MEFTAHI, AND K. STURM, *Shape optimization of an electric motor subject to nonlinear magnetostatics*, SIAM Journal on Scientific Computing, 37 (2015), pp. B1002-B1025.

[31] P. GANGL, K. STURM, M. NEUNTEUFEL, AND J. SCHÖBERL, *Fully and semi-automated shape differentiation in NGSolve*, Structural and Multidisciplinary Optimization, 63 (2020), pp. 1579-1607.

[32] M. GIACOMINI, O. PANTZ, AND K. TRABELSI, *Volumetric expressions of the shape gradient of the compliance in structural shape optimization*, arXiv preprint arXiv:1701.05762, (2017).

[33] D. J. GRIFFITHS, *Introduction to Electrodynamics*, Cambridge University Press, 2017.

[34] F. GUO AND I. P. BROWN, *Simultaneous magnetic and structural topology optimization of synchronous reluctance machine rotors*, IEEE Transactions on Magnetics, 56 (2020), pp. 1-12.

[35] F. GUO, M. SALAMEH, M. KRISHNAMURTHY, AND I. P. BROWN, *Multimaterial magneto-structural topology optimization of wound field synchronous machine rotors*, IEEE Transactions on Industry Applications, 56 (2020), pp. 3656-3667.

[36] A. HENROT AND M. PIERRE, *Shape Variation and Optimization*, EMS Tracts in Mathematics Vol. 28, 2018.

[37] R. HIPTMAIR, A. PAGANINI, AND S. SARGHEINI, *Comparison of approximate shape gradients*, BIT Numerical Mathematics, 55 (2015), pp. 459-485.

[38] J. A. IGLESIAS, K. STURM, AND F. WECHSUNG, *Two-dimensional shape optimization with nearly conformal transformations*, SIAM Journal on Scientific Computing, 40 (2018), pp. A3807-A3830.

[39] J. D. JACKSON, *Classical electrodynamics*, John Wiley & Sons, 2007.

[40] D. KALISE, K. KUNISCH, AND K. STURM, *Optimal actuator design based on shape calculus*, Mathematical Models and Methods in Applied Sciences, 28 (2018), p. 2667-2717.

[41] M. KOLMBAUER, *The multiharmonic finite element and boundary element method for simulation and control of eddy current problems*, PhD thesis, Ph. D. thesis, Johannes Kepler University Linz, 2012.

[42] J. KOLOTA AND S. STEPIEN, *Analysis of 2d and 3d finite element approach of a switched reluctance motor*, Przegla d Elektrotechniczny, 87 (2011), pp. 188-190.

[43] P. C. KRAUSE, O. WASYNCZUK, S. D. SUDHOFF, AND S. D. PEKAREK, *Analysis of electric machinery and drive systems*, vol. 75, John Wiley & Sons, 2013.

[44] E. KUCI, M. JANSEN, AND O. COULAUD, *Level set topology optimization of synchronous reluctance machines using a body-fitted mesh representation*, Structural and Multidisciplinary Optimization, 64 (2021), pp. 3729-3745.

[45] S. LANG, *Fundamentals of differential geometry*, vol. 191, Springer Science & Business Media, 2012.

[46] U. LANGER, S. E. MOORE, AND M. NEUMÜLLER, *Space-time isogeometric analysis of parabolic evolution problems*, Computer Methods in Applied Mechanics and Engineering, 306 (2016), pp. 342-363.

[47] U. LANGER AND A. SCHAFELNER, *Adaptive space-time finite element methods for parabolic optimal control problems*, Journal of Numerical Mathematics, 30 (2022), p. 247-266.

[48] LAURAIN, ANTOINE AND STURM, KEVIN, *Distributed shape derivative via averaged adjoint method and applications*, ESAIM: M2AN, 50 (2016), p. 1241-1267.

[49] J. LEE AND S. WANG, *Topological shape optimization of permanent magnet in voice coil motor using level set method*, IEEE transactions on magnetics, 48 (2012), pp. 931-934.

[50] J. LIONS, *Optimal control of systems governed by partial differential equations*, Grundlehren der mathematischen Wissenschaften, Springer-Verlag, 1971.

[51] F. LUCCHINI, R. TORCHIO, V. CIRIMELE, P. ALOTTO, AND P. BETTINI, *Topology optimization for electromagnetics: A survey*, IEEE Access, 10 (2022), pp. 98593-98611.

[52] B. MA, J. ZHENG, G. LEI, J. ZHU, P. JIN, AND Y. GUO, *Topology optimization of ferromagnetic components in electrical machines*, IEEE Transactions on Energy Conversion, 35 (2019), pp. 786-798.

[53] M. MERKEL, P. GANGL, AND S. SCHÖPS, *Shape optimization of rotating electric machines using isogeometric analysis*, IEEE Transactions on Energy Conversion, 36 (2021), pp. 2683-2690.

[54] F. MURAT AND J. SIMON, *Sur le contrôle par un domaine géométrique*, Pré-publication du Laboratoire d'Analyse Numérique, (76015), (1976).

[55] P. NEITTAANMÄKI, M. RUDNICKI, AND A. SAVINI, *Inverse problems and optimal design in electricity and magnetism*, vol. 35, Oxford University Press, 1996.

[56] A. A. NOVOTNY AND J. SOKOŁOWSKI, *Topological derivatives in shape optimization*, Springer Science & Business Media, 2012.

[57] ——, *An introduction to the topological derivative method*, Springer, 2020.

[58] A. PAGANINI AND F. WECHSUNG, *Fireshape: a shape optimization toolbox for firedrake*, Structural and Multidisciplinary Optimization, 63 (2021), pp. 2553-2569.

[59] C. PECHSTEIN AND B. JÜTTLER, *Monotonicity-preserving interproximation of B-H-curves*, Journal of Computational and Applied Mathematics, 196 (2006), pp. 45-57.

[60] R.-E. PLESSIX, *A review of the adjoint-state method for computing the gradient of a functional with geophysical applications*, Geophysical Journal International, 167 (2006), pp. 495-503.

[61] X. REN, A. THABUIS, A. BELAHCEN, AND Y. PERRIARD, *Topology optimization for coils of electric machine with level-set method*, in 2019 22nd International Conference on Electrical Machines and Systems (ICEMS), IEEE, 2019, pp. 1-4.

[62] A. A. RODRÍGUEZ AND A. VALLI, *Eddy Current Approximation of Maxwell Equations: Theory, Algorithms and Applications*, vol. 4, Springer Science & Business Media, 2010.

[63] N. SADOWSKI, Y. LEFEVRE, M. LAJOIE-MAZENC, AND J. CROS, *Finite element torque calculation in electrical machines while considering the mouvement*, Magnetics, IEEE Transactions on, 28 (1992), pp. 1410-1413.
[64] J. SCHÖBERL, C++11 implementation of finite elements in NGSolve, Tech. Rep. 30, Institute for Analysis and Scientific Computing, 2014.

[65] V. H. SCHULZ, M. SIEBENBORN, AND K. WELKER, *Structured inverse modeling in parabolic diffusion problems*, SIAM Journal on Control and Optimization, 53 (2015), pp. 3319–3338.

[66] O. SIGMUND AND K. MAUTE, *Topology optimization approaches*, Structural and Multidisciplinary Optimization, 48 (2013), pp. 1031–1055.

[67] J. SOKOLOWSKI AND J.-P. ZOLÉSIO, *Introduction to shape optimization*, Springer, 1992.

[68] O. STEINBACH, *Space-time finite element methods for parabolic problems*, Comput. Meth. Appl. Math., 15 (2015), pp. 551–566.

[69] V. THOMÉE, *Galerkin Finite Element Methods for Parabolic Problems*, Springer Berlin Heidelberg, 2006.

[70] H. TORKAMAN AND E. AFJEI, *Comprehensive study of 2-d and 3-d finite element analysis of*, Journal of Applied Sciences, 8 (2008), pp. 2758–2763.

[71] I. TOULOPOULOS, *Numerical solutions of quasilinear parabolic problems by a continuous space-time finite element scheme*, SIAM Journal on Scientific Computing, 44 (2022), pp. A2944–A2973.

[72] R. TOUZANI AND J. RAPPAZ, *Mathematical Models for Eddy Currents and Magnetostatics: With Selected Applications*, Springer Netherlands, 2014.

[73] B. WANG, *Topology optimization of electric machines: A review*, in 2022 IEEE Energy Conversion Congress and Exposition (ECCE), IEEE, 2022, pp. 1–8.

[74] E. ZEIDLER, *Nonlinear functional analysis and its applications: II/A: linear monotone operators*, Springer Science & Business Media, 2013.

[75] ——, *Nonlinear functional analysis and its applications: II/B: Nonlinear monotone operators*, Springer Science & Business Media, 2013.