# Air Gap Shape Optimization for Minimizing Proximity Losses in an Inductor

Stephane Gaydier, Itaï Zehavi, Théodore Cherrière, Peter Gangl

► To cite this version:

Stephane Gaydier, Itaï Zehavi, Théodore Cherrière, Peter Gangl. Air Gap Shape Optimization for Minimizing Proximity Losses in an Inductor. IEEE Transactions on Magnetics, 2026, pp.1-1.
10.1109/TMAG.2026.3651849. hal-05467608

HAL Id: hal-05467608

https://cnrs.hal.science/hal-05467608v1

Submitted on 20 Jan 2026

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 CC BY 4.0 - Attribution - International License
# Air Gap Shape Optimization
## for Minimizing Proximity Losses in an Inductor

S. Gaydier<fnref n="1" />,<fnref n="2" />,<fnref n="3" /> I. Zehavi <fnref n="4" />, T. Cherrière<fnref n="5" /> and P. Gangl<fnref n="1" />,<fnref n="6" />

In air-gapped inductors, the magnetic flux bulges outward from the gap - a phenomenon called the fringing effect. When these flux lines cut across the conductor windings, they generate undesired eddy currents that lead to substantial AC losses and reduced efficiency. To mitigate these losses, we develop a constrained shape optimization, aiming at finding an optimal air gap profile that minimizes conductor losses and keeps the inductance at a desired value. Our framework leverages established techniques such as complex permeability homogenization for loss modeling, shape sensitivity analysis coupling the adjoint method with the so-called Hilbertian regularization-extension method. The optimized air gap profile is curved to concentrate the flux lines away from the (fixed) conductors. Our implementation uses the open-source finite-element toolbox NGSolve and the Mmg remeshing library; all the codes are freely available.

*Index Terms*—Inductor, AC losses, Parametric optimization, Shape optimization

## I. INTRODUCTION

OPTIMAL design is, by now, a common practice for the conception of electromagnetic devices to gain in efficiency, mass, and cost. Among them, inductors are part of the basic components in electrical engineering: they are widely used in filters and power electronics for DC-DC conversion to mitigate current ripple and obtain better signal quality. It is often preferable to reduce the inductor size to keep a low mass by increasing the switching frequency, which comes with additional losses in the switches, magnetic material, and conductors [1]. This work focuses on inductor designs that contain an air gap to prevent saturation of the soft magnetic core [2]. With the frequency rise, the magnetic flux lines coming out of the air gap due to the fringing effect cross the conductors, which creates eddy currents and additional losses known as proximity losses [3]. Therefore, such air-gapped inductors have to be designed carefully, using, for instance, optimization algorithms.

The first possibility is to optimize the conductors directly. It is possible to attenuate the fringing flux with a shield [4], but other losses appear in the shield. Some work has been done on their positioning within the winding window [5]. Clearly, the wires should be located as far as possible from the air gap. The design can be further refined by optimizing the shape of the conductors. It is well known that flat foil conductors should be oriented in parallel to magnetic flux lines to minimize the proximity effect. From this principle, [6], [7] propose to build bent foils that follow the natural curvature of the flux using additive manufacturing, resulting in minimal losses for a given magnetic core. Since the shape of the conductor has a very small influence on the flux lines, the conceptual guidelines of optimal conductor design seem to be well understood.

The second component to optimize is the magnetic core. Many approaches are available in the literature. On the one hand, a distributed air gap obtained by using a segmented core or soft magnetic composite is the best solution to avoid the fringing effect; however, it is the most expensive. On the other hand, discrete air-gapped cores are easy to manufacture and cheap [8]. The additional losses due to the fringing effect can be mitigated by a proper design of the air gap profile, which becomes possible due to the development of additive manufacturing [9]. Previous work [10] has optimized the air gap profile, showing that the best air gap shape seems to be a cone, which keeps the flux lines away from the conductors. However, the gradient-free optimization algorithm used by [10] limits the number of control points and therefore the resolution of the air gap profile. Indeed, the CPU cost required to compute the derivative of the objective functions using finite differences increases linearly with the number of parameters. To avoid computing finite difference approximations of the sensitivities, one can use the adjoint method. Using shape derivatives eliminates the need for intrinsic parametrization of the air gap profile [11], [12].

This paper provides a concrete application of the adjoint method and shape derivatives to minimize losses in an inductor due to the fringing effect. We compare three approaches: a gradient-free and a gradient-based method that employ control points and a gradient-based adaptive mesh deformation strategy. Section II exposes the optimization problem and the sensitivity computation. The algorithms and their implementation details are provided in Section III. Our code is written in Python using the NGSolve package and is freely available on [13]. The performances and results of our approaches are compared in Section IV, and Section V concludes this work.

<footnote n="1">¹Schneider Electric, Global Technology 38320 Eybens, France, stephane.gaydier@g2elab.grenoble-inp.fr</footnote>

<footnote n="2">²Université Grenoble Alpes, CNRS, Grenoble INP, G2Elab, 38000 Grenoble, France</footnote>

<footnote n="3">³Université Grenoble Alpes, INRIA, CNRS, Grenoble INP, LJK, 38000 Grenoble, France</footnote>

<footnote n="4">⁴Electrical Engineering Department, ENS Paris-Saclay, Université Paris-Saclay, 91190 Gif-sur-Yvette, France</footnote>

<footnote n="5">⁵Université Paris-Saclay, CentraleSupélec, Sorbonne Université, CNRS, GeePs, 91192 Gif-sur-Yvette, France</footnote>

<footnote n="6">⁶Johann Radon Institute for Computational and Applied Mathematics (RICAM), ÖAW, 4040 Linz, Austria</footnote>
## II. MATHEMATICAL MODEL

### A. Physical equations

We consider a simple model for the magnetic properties within a 2D region $\Omega = \Omega_a \cup \Omega_f \cup \Omega_c$, in the time-harmonic regime, at frequency $f = 50$ kHz, as illustrated in Figure 1. Here,

* $\Omega_a$ is the air region where the magnetic permeability is $\mu_0 = 4\pi \times 10^{-7}$ H/m;

* $\Omega_f$ is filled with a soft magnetic material with permeability $\mu_f = 1000\mu_0$. The iron losses are neglected in this region;

* $\Omega_c$ represents electrical conductors. Its complex-valued permeability comes from homogenized models depending on the frequency and configuration of the wires. For the sake of simplicity, we take a simple model $\mu_c = \mu_0 \exp(-i\delta)$ with $\delta = 0.1$ rad for the frequency $f = 50$ kHz, $i$ being the unit imaginary number i.e $i^2 = -1$. We underline each complex-valued variable here and throughout the article. More elaborate models are possible [14], [15].

The state equation then reads

$$ \begin{cases} -\operatorname{div}(\underline{\nu}\nabla\underline{a}) = j, \\ \underline{a} = 0 \text{ on } \Gamma_D \text{ and } \frac{\partial\underline{a}}{\partial n} = 0 \text{ on } \Gamma_N, \end{cases} \quad (1) $$

where the reluctivity $\underline{\nu}$ equals $\frac{1}{\mu_c}$ in $\Omega_c$, $\frac{1}{\mu_f}$ in $\Omega_f$ and $\frac{1}{\mu_0}$ in $\Omega_a$, $\underline{a}$ is the out-of-the-plane component of the magnetic vector potential such that $\underline{b} = \begin{bmatrix} 0 & 1 \\ -1 & 0 \end{bmatrix} \nabla\underline{a}$, $j$ the real-valued current density that is the phase reference and assumed to be sinusoidal, $\Gamma_D$ and $\Gamma_N$ the portions of the boundary depicted by Figure 1. In this device, the AC losses come from the proximity losses occurring in the conductor. It is accounted by the imaginary part of the permeability [16]; more specifically, it equals

$$ P(\underline{a}) = sl_z \pi f \int_{\Omega_c} \operatorname{Im} \left( \frac{1}{\mu_c} \right) |\nabla \underline{a}|^2 dx, \quad (2) $$

while the inductance expressed from the reactive power reads

$$ L(\underline{a}) = \frac{sl_z}{I^2} \int_{\Omega} \operatorname{Re}(\underline{\nu}) |\nabla \underline{a}|^2 dx, \quad (3) $$

where $s = 4$ is the symmetry coefficient, $l_z = 1$ cm is the z-thickness (fringing is neglected in the z direction), and $I$ is the electric current amplitude set to 2 A with 200 turns. This model is discretized and solved using NGSolve.

### B. Parametrization and minimization problem

Let $h$ be a parametrization of the air gap profile. For example, $h$ may stand for a set of control points (as in our first two optimization approaches) or a NURBS representation of the interfaces [17]. Alternatively, $h$ can be the shapes $(\Omega_a, \Omega_f, \Omega_c)$ of the regions at play, or it can be their numerical representation, which is the collection of the mesh nodes (as in our last optimization approach). Thus, the following objects, some of which are depicted in Fig. 1, depend on $h$:

<figure><img src="image_2.png" /><figcaption>Fig. 1: Simulation space reduced to one quarter by symmetry.</figcaption></figure>

* The computational domain partition $\Omega = \Omega_{a,h} \cup \Omega_{f,h} \cup \Omega_{c,h}$, which represents the air, magnetic, and conductor regions, respectively;

* The reluctivity $\underline{\nu}_h$, defined as $\frac{1}{\mu_c}$ in $\Omega_{c,h}$, $\frac{1}{\mu_f}$ in $\Omega_{f,h}$, and $\frac{1}{\mu_0}$ in $\Omega_{a,h}$;

* The magnetic vector potential $\underline{a}_h$, which is the solution to the variational problem

$$ \begin{gathered} \text{Find } \underline{a}_h \in H := \{\underline{a} \in H^1(\Omega) \mid \underline{a} = 0 \text{ on } \Gamma_D\} \text{ s.t.} \\ \forall \underline{v} \in H, \int_{\Omega} \underline{\nu}_h(\nabla \underline{a}_h) \cdot (\nabla \underline{v})^* dx = \int_{\Omega_{c,h}} j \underline{v}^* dx \quad (4) \end{gathered} $$

where $H^1(\Omega) = \{\underline{a} \in L_2(\Omega) \mid \nabla \underline{a} \in L_2(\Omega)\}$ is a Sobolev space and “*” is the notation for complex conjugates.

The optimal design problem at hand reads:

$$ \min_{h} P(h, \underline{a}_h) \quad \text{subject to} \quad \left\{ \begin{array}{l} L(h, \underline{a}_h) = 1 mH \\ \underline{a}_h \text{ solution of (4)} \end{array} \right. \quad (5) $$

In our gradient-based methods, we solve (5) using a gradient descent algorithm that is based on the derivatives of the PDE-constrained functionals $P(h, \underline{a}_h)$, $L(h, \underline{a}_h)$. Their calculation is realized thanks to the adjoint method.

### C. Adjoint method

In shape and topology optimization, the design variable $h$ often lives in high-dimensional space, making the use of finite difference methods to compute the derivatives of $P$ and $L$ unaffordable. Fortunately, the adjoint method [11], [12] offers an efficient alternative. It introduces two auxiliary variables, called adjoint states, denoted in this work as $\underline{p}_P$ and $\underline{p}_L$. (We omit their dependency on $h$ for notational simplicity.)

In these adjoint equations, the right-hand sides involve the derivatives of $P$ and $L$ with respect to the state variable, and the bilinear operator at the left-hand side is the adjoint (Hermitian conjugate) of the (linearized) operator of the state equation.

In the context of our work, the adjoint state $\underline{p}_P$ solves the variational equation

$$ \forall \underline{v} \in H, \int_{\Omega} \underline{\nu}_h^* \nabla \underline{v} \cdot \nabla \underline{p}_P dx = - \frac{\partial P}{\partial \underline{a}}(\underline{a})(\underline{v}), \quad (6) $$

and $\underline{p}_L$ is defined analogously with
$$ \forall \underline{v} \in H, \int_{\Omega} \underline{v}_h^* \nabla \underline{v} \cdot \nabla \underline{p}_L dx = - \frac{\partial L}{\partial \underline{a}}(\underline{a})(\underline{v}). \quad (7) $$

These two equations involve straightforward (complex-valued) Fréchet-derivatives of real-valued functions of a complex-valued argument. The next part provides explicit formulas for the sensitivities of P and L to a change in h involving $\underline{p}_P$ and $\underline{p}_L$.

### D. Shape derivatives

The sensitivities to a variation of the air gap profile can be evaluated using shape derivatives [11], [12]. In a nutshell, the shape derivative of a functional $g(\Omega)$, depending on a shape $\Omega \subset \mathbb{R}^2$, is a linear form defined on a deformation vector field $\phi \in W^{1,\infty}(\mathbb{R}^2) := \{\phi \in L^\infty(\mathbb{R}^2) | \nabla\phi \in L^\infty(\mathbb{R}^2)\}$, illustrated in 2, by the limit

$$ dg(\Omega)(\phi) = \lim_{t \to 0} \frac{g(\Omega + t\phi) - g(\Omega)}{t}. $$

The shape derivatives of the AC losses and the inductance read

$$ dP(\Omega)(\phi) = \operatorname{Re} \left( sl_z \pi f \int_{\Omega_c} \operatorname{Im} \left( \frac{1}{\mu_c} \right) (A_\phi \nabla \underline{a}) \cdot (\nabla \underline{a})^* dx \right. \\ \left. + \int_{\Omega} \underline{v}(A_\phi \nabla \underline{a}) \cdot \nabla \underline{p}_P^* dx - \int_{\Omega_c} j \operatorname{div} \phi \underline{p}_P^* dx \right) \quad (8) $$

and

$$ dL(\Omega)(\phi) = \operatorname{Re} \left( \frac{sl_z}{I^2} \int_{\Omega} \underline{v}(A_{\phi} \nabla \underline{a}) \cdot (\nabla \underline{a})^* dx \right. \\ \left. + \int_{\Omega} \underline{v}(A_{\phi} \nabla \underline{a}) \cdot \nabla \underline{p}_L^* dx - \int_{\Omega_c} j \operatorname{div} \phi \underline{p}_L^* dx \right) \quad (9) $$

where $\underline{p}_P$ and $\underline{p}_L$ are the adjoint states and are the solutions of the respective equations (6) and (7) and with the real $2 \times 2$ matrix $A_\phi = \operatorname{div}\phi I_2 - \nabla\phi - \nabla\phi^\top$ ($I_2$ is the $2 \times 2$ identity matrix).

Appendix A offers motivation for the introduction of the adjoint equations (6) and (7) and some proof elements for the derivation of the shape derivatives (8) and (9).

## III. IMPLEMENTATION DETAILS

Our shape optimization workflow is sketched in Figure 6. The detailed flowcharts for our particular optimization methods can be found in appendix B. The present section details a few additional important issues.

### A. Augmented Lagrangian

Our treatment of the constrained optimization problem (5) hinges on a simple augmented Lagrangian strategy [18].

The augmented Lagrangian method converts (5) into a series of problems indexed by $n \in \mathbb{N}$,

$$ \min_{h} J(h, \ell^{n}, b^{n}), \text{ where} $$

$$ J(h, \ell, b) := P(h) + \ell(L(h) - L_0) + \frac{b}{2}(L(h) - L_0)^2, \quad (10) $$

<figure><img src="image_3.png" /><figcaption>Fig. 2: Shape variation from boundary perturbation</figcaption></figure>

with $L_0 = 1 \, \text{mH}$. The equality constraint of (5) has been eliminated. To keep this presentation light, we omitted the dependence of P and L on $\underline{a}_h$.

After each iteration of the optimization loop, we update the $\ell$ and $b$ parameters according to the rule (with $\alpha > 1$ and $b_{\text{target}} > 0$):

$$ \ell^{n+1} = \ell^n + b^n (L(h) - L_0) \quad (11) $$

$$ \text{and } b^{n+1} = \begin{cases} \alpha b^n & \text{if } b^n < b_{\text{target}}, \\ b^n & \text{otherwise.} \end{cases} \quad (12) $$

### Algorithm 1 Augmented Lagrangian Method

**Initialization:** $h^0, \ell^0, b^0$.

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

Find a descent direction $\phi^n$ for $h \mapsto J(h, l^n, b^n)$ as in III-B.

Calculate $h^{n+1} = h^n - \tau^n \phi^n$ where $\tau^n$ is small enough such that $J(h^{n+1}, l^n, b^n) < J(h^n, l^n, b^n)$.

Calculate $\ell^{n+1}$ and $b^{n+1}$ following (11) and (12).

**end for**

**Return:** $h^n$.

### B. Hilbertian extension-regularization

A descent direction for J is defined as a deformation field $\phi$ that satisfies

$$ dJ(\Omega)(\phi) < 0, $$

where $dJ(\Omega)$ can be calculated as a combination of $dP(\Omega)$ and $dL(\Omega)$.

Such a $\phi$ is not immediately apparent from the linear forms (8) and (9). To obtain one, we can solve the following auxiliary problem:

$$ \text{Find } \phi \in V \text{ such that } \forall \hat{\phi} \in V, \langle \phi, \hat{\phi} \rangle_V = dJ(\Omega)(\hat{\phi}), \quad (13) $$

$$ \text{with } \langle \phi, \hat{\phi} \rangle_V := \alpha^2 \langle \nabla\phi, \nabla\hat{\phi}^* \rangle_{L^2(\Omega)^{2 \times 2}} + \langle \phi, \hat{\phi}^* \rangle_{L^2(\Omega)^2} $$

where $V$ is a subspace of $H^1(\Omega)^2$ with suitable Dirichlet boundary and $\alpha > 0$ is a regularization parameter with the dimension of a length [12]. Therefore, $-\phi$ is a direction of descent because

$$ dJ(\Omega)(-\phi) = -\langle \phi, \phi \rangle_V < 0. $$

This procedure has a regularizing effect that can be controlled by the choice of $V$ and its Hermitian product and can be interpreted as a convolution filter applied on $\phi$. We set $\alpha$ to twice the maximum mesh element size. This regularization
ensures a globally smooth deformation field $\phi$, preserving mesh quality and preventing distortions such as inverted or collapsed elements.

### C. Dirichlet conditions on the deformation field

Imposing Dirichlet boundary conditions in the space $V$ ensures that some components of the descent direction vanish on portions of the interfaces. The source term of the magneto-harmonic equation (4), in the domain $\Omega_c$, must stay the same during all the optimization iterations. We also fix interfaces that are not part of the air-gap profile, and we restrict the deformation of the latter to its vertical component. This lowers the chance of flipping triangles when we deform the mesh using our "adaptive mesh deformation approach," which is described in the next part.

## IV. COMPARISON OF OPTIMIZATION APPROACHES

In this section, the results of our optimizations are given and compared to a reference design shown in Figure 1. Three methodologies are assessed: a classical gradient-free algorithm, a gradient-based algorithm that employs control points, and a third method that employs adaptive mesh deformations and is also gradient-based. A summary is provided in Table I and a flowchart for the gradient-based methods is given in Appendix B.

### A. Gradient-free optimization with control points

The easiest option to optimize the air gap profile is to parametrize a set of control points on the interface, and then use a black-box solver from a standard library to optimize their position. For reference, we implemented this test case using the differential evolution algorithm from the `scipy.minimize` package [19], with a population size of 15 individuals and a maximum iteration number of 100. The constraint is considered to be satisfied if $L = 1 \text{ mH} \pm 5\%$. The parameters `atol` and `tol` are both set to 0.1. We observe that the computational cost seems to increase linearly with the number of control points (tested until 20), with almost no improvement on the losses. Figure 4a shows a design with 10 control points that optimized using differential evolution. To accelerate the computation, a gradient-based approach can be used.

### B. Gradient-based optimization with control points

Again, we start by defining four (or more) control points $p_i$ as shown in Figure 3. An iteration of the augmented Lagrangian loop starts by creating a mesh according to the control point positions. This allows us to solve the magneto-harmonic equation (1) using the finite element method and calculate the cost function that results from (10). Then we check the convergence of the optimization and eventually halt. If not, we evaluate the discrete gradient of the Lagrangian with respect to the positions of the control points. This is done by applying the displacement field induced by a variation of each control point position to the shape derivative (understood as a linear form, see (II-D))

<figure><img src="image_4.png" /><figcaption>Fig. 3: Control points optimization setup</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 4: Optimization results with 10 control points.</figcaption></figure>

This leads to a parametric optimum. Note that artifacts may occur without regularization. In particular, it is necessary to enforce the symmetry condition by imposing the boundary to be orthogonal to the symmetry axis; otherwise, numerical instabilities lead to undesired local optima. The orthogonality is ensured by choosing a descent direction with equal values associated with the control point on the symmetry axis and its immediate neighbor. The design that emerged from the gradient-based optimization with 10 control points is shown in 4b and its optimization results are given in Table I. The resulting air gap profile is shaped to carry high flux away from the conductors, halving AC losses, and leads to lower losses than the one obtained with differential evolution for a fraction of its computational cost.

On the one hand, the number of iterations is significantly reduced compared with the differential evolution approach, and does not depend on the number of optimization variables. On the other hand, the AC losses values are also lower, and decrease even more when the number of control points increases; an example with 30 control points is shown in Figure 5a. Therefore, the best design should be obtained considering a continuous air gap profile.

### C. Adaptive mesh deformation optimization

The optimization follows the same path as before, but instead of creating a new mesh at each loop iteration, which was expensive, we start by initializing the mesh that will be deformed. Moreover, with this method, $\phi$ is a deformation field that takes values at each mesh node, substantially increasing the number of effective control points and leading to greater design diversity.

Furthermore, the mesh quality may drastically deteriorate at the end of the loop after we deform it. To remedy this, we use the MMG software [20] to remesh every five iterations of
<figure><img src="image_6.png" /><figcaption>Fig. 5: Final results of the high-resolution optimizations.</figcaption></figure>

the optimization loop to keep the mesh quality acceptable. We set the smallest and largest mesh sizes with a maximum ratio of 10 and refine the elements when the interfaces are curved. For that to happen, we had to convert the Netgen mesh data structure into an MMG data structure and vice versa.

The optimized design in Figure 5b shows a clean curved air gap profile as suggested in [10], achieving even lower AC losses, converging with a similar computation time as the gradient-based control points approach.

The obtained air gap profile is similar, though not identical, to the profile obtained with 30 control points shown in Figure 5a. Regarding computational efficiency, the mesh adaptive method required approximately twice as many magnetostatic finite element solves to converge, as step sizes were constrained to ensure mesh integrity. However, this was offset by the fact that each individual iteration was roughly two times cheaper, owing to the elimination of the meshing step. Consequently, both gradient-based methods completed in under 3 minutes, which is significantly faster than the zero-order method that requires over an hour. Ultimately, this technique yielded slightly better performance in terms of AC losses. It remains an open question whether this improvement stems from superior optimizer convergence, the greater number of degrees of freedom afforded by the mesh deformation, or other numerical effects.

During optimization, the two components of the air gap tend to close together and form tips. To get around that, we halt the optimization when the inductance crosses the $L_0 = 1$ mH value after some arbitrary number of iterations before the air gap closes. The ferromagnetic material may eventually become saturated as a side-effect, but this is not the case for our design at the selected operating point, where, luckily, the induction field stays between 0 and 1 T. Of course, this would become an issue as soon as we increase the current intensity.

The adaptive mesh deformation approach can be improved by considering the ferromagnetic saturation using nonlinear modeling, which could enable a proper optimality stopping criterion. However, in this case, the complex permeability homogenization must be dropped since it relies on the linearity of the magneto-harmonic equation.

## V. CONCLUSION

This paper explains the methodologies of control points and adaptive mesh deformation for shape optimization, effectively

<table><caption>TABLE I: Optimization results (the number of iterations of the differential evolution increases with the number of points).</caption><thead><tr><th>Optimization Type</th><th>Magn. solves</th><th>P (W)</th><th>L (mH)</th></tr></thead><tbody><tr><td>Reference</td><td>-</td><td>13.16</td><td>1.00</td></tr><tr><td>Differential evolution (avg)</td><td>6385</td><td>7.85</td><td>1 mH ± 5%</td></tr><tr><td>Gradient-based control points</td><td>310</td><td>3.80</td><td>1.00</td></tr><tr><td>Adaptive mesh deformation</td><td>604</td><td>3.76</td><td>1.00</td></tr></tbody></table>

reducing proximity losses while maintaining the desired induction. The use of a gradient method through shape differentiation enables it to be significantly faster than differential evolution, which is a popular, ready-to-use optimization algorithm. Adaptive mesh deformation achieved slightly superior results with a similar computational cost but requires remeshing. Since gradient methods are more difficult to implement in general, we provide the Python codes [13] relying on NGSolve and the MMG software. In future work, we plan to consider the nonlinearity of iron as well, its losses, together with conductor design optimization as in [7] using the same framework.

## ACKNOWLEDGMENT

The work of P. Gangl was partly supported by the Collaborative Research Centre CREATOR of the German Research Foundation (DFG) and the Austrian Science Fund (FWF) at Technische Universität Darmstadt. Technische Universität Graz, Johannes Kepler Universität Linz and Johann Radon Institute for Computational and Applied Mathematics (RICAM) Linz, under Grant 492661287/TRR361 (DFG) and Grant 10.55776/F90 (FWF). The work of P. Gangl was additionally supported in part by the State of Upper Austria.

## REFERENCES

[1] C. W. T. McLyman, Transformer and Inductor Design Handbook. CRC Press, Mar. 2004. [Online]. Available: http://dx.doi.org/10.1201/9780203913598

[2] F. Battal, "A comparative analysis of core material and gap sizing effect on the high-power inductor design," Engineering Science and Technology, an International Journal, vol. 63, p. 102001, Mar. 2025. [Online]. Available: http://dx.doi.org/10.1016/j.jestch.2025.102001

[3] R. Ridley, "Proximity loss in magnetics windings," Switching Power Magazine, vol. Designers Series XIII, 2005

[4] T. Ewald and J. Biela, "Analytical model and concept study of an air gap fringing shield to reduce losses in gapped inductors," IEEE Transactions on Power Electronics, vol. 39, no. 10, pp. 13487-13496, 2024.

[5] R. Jensen and C. Sullivan, "Optimal core dimensional ratios for minimizing winding loss in high-frequency gapped-inductor windings," in 18th IEEE APEC, 2003, pp. 1164-1169 vol.2.

[6] N. Simpson and P. H. Mellor, "Additive manufacturing of shaped profile windings for minimal ac loss in electrical machines," in 2018 IEEE Energy Conversion Congress and Exposition (ECCE), 2018, pp. 5765-5772.

[7] N. Simpson, C. Tighe, and P. Mellor, "Design of high performance shaped profile windings for additive manufacture," in 2019 IEEE Energy Conversion Congress and Exposition (ECCE), 2019, pp. 761-768.

[8] M. Beraki, J. P. Trovão, and M. Perdigão, "Comprehensive comparison and selection of magnetic materials for powertrain dc-dc converters," IET Electrical Systems in Transportation, vol. 10, no. 2, p. 125-134, Jun. 2020. [Online]. Available: http://dx.doi.org/10.1049/iet-est.2019.0113

[9] T.-W. Chang, P.-W. Huang, C.-H. Shih, A. Hung, C.-C. Wang, H.-F. Ko, and M.-C. Tsai, "Additive manufacturing high current power inductor," IEEE Transactions on Magnetics, vol. 60, no. 9, pp. 1-5, 2024.

[10] D. I. Zaikin, S. Jonasen, and S. L. Mikkelsen, "An air-gap shape optimization for fringing field eddy current loss reductions in power magnetics," IEEE Trans. Power Electron., vol. 34, no. 5, pp. 4079-4086, 2019.

[11] P. Gangl, K. Sturm, M. Neunteufel, and J. Schöberl, "Fully and semi-automated shape differentiation in NGSolve," Struct. Multidiscip. Optim., vol. 63, no. 3, pp. 1579-1607, 2021.
[12] G. Allaire, C. Dapogny, and F. Jouve, “Shape and topology optimization,” in Handbook of Numerical Analysis. Elsevier, 2021, vol. 22, pp. 1–132. [Online]. Available: https://linkinghub.elsevier.com/retrieve/pii/S1570865920300181

[13] S. Gaydier, I. Zehavi, T. Cherrière, and P. Gangl, “Stephlemedef/inductance-compumag: v0.1.2,” 2025, github repository. [Online]. Available: https://zenodo.org/doi/10.5281/zenodo.14631426

[14] X. Nan and C. Sullivan, “A two-dimensional equivalent complex permeability model for round-wire windings,” in 2005 IEEE 36th Power Electronics Specialists Conference, 2005, pp. 613–618.

[15] K. Niyomsatian, J. Gyselinck, and R. V. Sabariego, “Closed-form complex permeability expression for proximity-effect homogenisation of litz-wire windings,” IET Sci. Meas. Technol., vol. 14, no. 3, pp. 287–291, 2020.

[16] O. Moreau, L. Popiel, and J. Pages, “Proximity losses computation with a 2d complex permeability modelling,” IEEE Transactions on Magnetics, vol. 34, no. 5, p. 3616–3619, 1998. [Online]. Available: http://dx.doi.org/10.1109/20.717854

[17] M. Merkel, P. Gangl, and S. Schöps, “Shape optimization of rotating electric machines using isogeometric analysis,” IEEE Transactions on Energy Conversion, vol. 36, no. 4, pp. 2683–2690, 2021.

[18] J. Nocedal and S. J. Wright, Numerical optimization, ser. Springer series in operations research and financial engineering. New York, NY: Springer, 2006.

[19] P. Virtanen, R. Gommers, and T. E. t. Oliphant, “SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python,” Nature Methods, vol. 17, pp. 261–272, 2020.

[20] C. Dapogny, C. Dobrzynski, and P. Frey, “Three-dimensional adaptive domain remeshing, implicit domain meshing, and applications to free and moving boundary problems,” Journal of Computational Physics, vol. 262, pp. 358–378, Apr. 2014, publisher: Elsevier BV. [Online]. Available: https://linkinghub.elsevier.com/retrieve/pii/S0021999114000266

# APPENDIX A

## DERIVATION OF SHAPE DERIVATIVE WITH CÉA’S METHOD

We follow Jean Céa’s method to obtain the volume form of the shape derivative; see [11] and [12] for a more comprehensive explanation. The volume form of the shape derivatives for the AC losses P and the inductance L can be derived by differentiating perturbed Lagrangians with respect to the shape parameter t at t = 0.

Let $\mathcal{L}_P$ and $\mathcal{L}_L$ denote the Lagrangians associated with P and L, respectively. They are constructed to satisfy the following:

1) The magnetostatic state equation (4) is enforced as a constraint using $\underline{p}$ as a Lagrange multiplier.

2) The Lagrangians must be defined on a fixed reference domain to avoid differentiating over moving domains. We use the substitutions $\underline{a} \circ T_t^{-1}$ and $\underline{p} \circ T_t^{-1}$ with the change of variables $y = T_t(x)$, where $\overline{T_t}(x) = x + t\phi(x)$ and $\phi \in W^{1,\infty}(\mathbb{R}^2)$.

3) For any $\underline{p}$, evaluating the Lagrangians at the magnetostatic state $\underline{a}_t$ for the deformed subdomains must recover the original objective and constraint functions:

$$ 
\begin{aligned}
P(t) &= \mathcal{L}_P(t, \underline{a}_t \circ T_t^{-1}, \underline{p} \circ T_t^{-1}), \\
L(t) &= \mathcal{L}_L(t, \underline{a}_t \circ T_t^{-1}, \underline{p} \circ T_t^{-1}).
\end{aligned} 
\quad (14) $$

We set $\xi(t) = \det(\nabla T_t)$ and define $A(t) = \xi(t) (\nabla T_t)^{-1} (\nabla T_t)^{-T}$. The Lagrangians are then given by

$$ 
\begin{aligned}
\mathcal{L}_P(t, \underline{a} \circ T_t^{-1}, \underline{p} \circ T_t^{-1}) \\
&=sl_z \pi f \int_{\Omega_c} \mathrm{Im} \left( \frac{1}{\mu_c} \right) (A(t) \nabla \underline{a}) \cdot (\nabla \underline{a})^* dx \\
&\quad + \int_{\Omega} \underline{\nu}(A(t) \nabla \underline{a}) \cdot (\nabla \underline{p})^* dx - \int_{\Omega_c} \xi(t) (j \circ T_t) \underline{p}^* dx
\end{aligned} 
$$

and

$$ 
\begin{aligned}
\mathcal{L}_L(t, \underline{a} \circ T_t^{-1}, \underline{p} \circ T_t^{-1}) \\
&= \frac{sl_z}{I^2} \int_{\Omega} \mathrm{Re}(\underline{\nu})(A(t)\nabla \underline{a}) \cdot (\nabla \underline{a})^* dx \\
&\quad + \int_{\Omega} \underline{\nu}(A(t)\nabla \underline{a}) \cdot (\nabla \underline{p})^* dx - \int_{\Omega_c} \xi(t)(j \circ T_t) \underline{p}^* dx.
\end{aligned} 
$$

The core principle of Céa’s method is to introduce the adjoint states $\underline{p}_P$ and $\underline{p}_L$ so that $(\underline{a}, \underline{p}_L)$ and $(\underline{a}, \underline{p}_L)$ are stationary points of $(\underline{a}, \underline{p}) \mapsto \mathcal{L}_P(0, \underline{a}, \underline{p})$ and $(\underline{a}, \underline{p}) \mapsto \mathcal{L}_L(0, \underline{a}, \underline{p})$ respectively. This reflects the KKT conditions:

$$ \frac{\partial \mathcal{L}}{\partial \underline{a}} = 0, \quad \frac{\partial \mathcal{L}}{\partial \underline{p}} = 0, $$

ensuring the state and adjoint equations hold.

Consequently, the shape derivatives are then obtained by differentiating (14) with respect to t at t = 0 and applying the chain rule:

$$ 
\begin{aligned}
dP(\Omega)(\phi) &= \frac{\partial P(t)}{\partial t}\bigg|_{t=0} = \frac{\partial \mathcal{L}_P}{\partial t}(0, \underline{a}, \underline{p}_P) \\
&\quad + \underbrace{\frac{\partial \mathcal{L}_P}{\partial \underline{a}}(0, \underline{a}, \underline{p}_P)}_{=0} \left( \frac{\partial(\underline{a}_t \circ T_t^{-1})}{\partial t} \bigg|_{t=0} \right) \\
&= \frac{\partial \mathcal{L}_P}{\partial t}(0, \underline{a}, \underline{p}_P),
\end{aligned} 
$$

and similarly for $dL(\Omega)(\phi)$.

The final expressions (8) and (9) are obtained by substituting the identities $\xi'(0) = \mathrm{div}(\phi)$ and $A'(0) = A_\phi = \mathrm{div}(\phi)I_2 - \nabla\phi - \nabla\phi^T$.

# APPENDIX B

## FLOWCHART

<figure><img src="image_7.png" /><figcaption>Fig. 6: Air gap profile shape optimization’s flowchart</figcaption></figure>