# Volume integral eddy current problem as a Friedrichs' system

Anthony Torosyan, Bertrand Bannwarth, Olivier Chadebec, Jean-Michel

Guichon, Gérard Meunier, Jean-René Poirier

► To cite this version:

Anthony Torosyan, Bertrand Bannwarth, Olivier Chadebec, Jean-Michel Guichon, Gérard Meunier, et al.. Volume integral eddy current problem as a Friedrichs' system. 25th International Conference on the Computation of Electromagnetic Fields (COMPUMAG'2025), Jun 2025, Napoli, Italy. hal-05135458

HAL Id: hal-05135458

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

Submitted on 30 Jun 2025

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
# Volume integral eddy current problem as a Friedrichs' system

Anthony Torosyan¹, Bertrand Bannwarth¹, Olivier Chadebec¹, Jean-Michel Guichon¹,
Gérard Meunier¹ and Jean-René Poirier²

¹Univ. Grenoble Alpes, CNRS, Grenoble INP, G2Elab, 38000 Grenoble, France

²LAPLACE, Université de Toulouse, CNRS, INPT, UPS, Toulouse, France

The paper introduces Discontinuous Galerkin discretization of a volume integral eddy current problem by recasting it into a Friedrichs' system. A numerical example is given at the end.

**Index Terms**—Volume integral equations, Discontinuous Galerkin, Eddy currents problem.

## I. INTRODUCTION

THE Volume Integral Method is a general framework for solving electromagnetic problems. It permits to discretize only the active parts of the geometry but leads to full matrix systems (i.e. with almost only non zeros coefficients). For eddy current problems, the formulation introduced in [1] contains a volume integral part and is usually discretized with the Partial Element Equivalent Circuit (PEEC) method. It is a conforming approximation with a circuit approach to ensure the condition on the divergence of the current density.

The theory of Friedrichs systems' is a general framework for the discretization of systems of order one partial differential equations. It was introduced by Friedrich [2] in 1958. In the 2000s, Ern and Guermond in a series of three papers [3], [4] and [5] give a general theory for the Discontinuous Galerkin discretization of these systems in the Finite Element context. In the paper, we show that the theory of Friedrichs' systems can be used to describe the volume integral eddy current problem and leads almost naturally to a Discontinuous Galerkin discretization. To our knowledge, this new framework is original and allows non-conforming meshes with basis functions of any order on each element. Moreover, this discretization could be used efficiently to develop compression methods for the volume integral part by tensors representation (e.g. tensor train [6]) in order to deal with bigger problems and to obtain more accuracy.

## II. THEORETICAL FRAMEWORK

### A. The volume integral eddy current problem

Let $\Omega \subset \mathbb{R}^3$ be a bounded Lipschitz domain. Assume $a_0 \in L^2(\Omega; \mathbb{C}^3)$ (vector potential) is given. Following the work of [1], we want to find $(\mathbf{j}, u) \in H(\mathrm{div}, \Omega) \times H^1(\Omega)$ (current density, scalar potential) such that

$$ \left\{ \begin{array}{l} \rho \mathbf{j} + i\omega\mu g * \mathbf{j} + \nabla u = -i\omega a_0, \quad \text{on } \Omega \\ \nabla \cdot \mathbf{j} = 0, \quad \text{on } \Omega \\ \mathbf{j} \cdot \mathbf{n} = 0, \quad \text{on } \partial\Omega \end{array} \right. \quad (1) $$

at a given frequency $\nu > 0$ (i.e. $\omega := 2\pi\nu$) and with $\rho, \mu \in L^\infty(\Omega; \mathbb{R})$ (resistivity, permeability) such that

$$ \mathrm{ess\,inf}_{\mathbf{x} \in \Omega} \rho(\mathbf{x}) := \rho_{\#} > 0, \quad \mathrm{ess\,inf}_{\mathbf{x} \in \Omega} \mu(\mathbf{x}) := \mu_{\#} > 0. \quad (2) $$

We denote by $*$ the convolution and we define the Green kernel of the Laplace equation in $\mathbb{R}^3$ by $g(\mathbf{x}, \mathbf{y}) := 1/(4\pi|\mathbf{x} - \mathbf{y}|^2)$ ($|\cdot|_2$ is the euclidean norm). The normal vector field of the boundary of $\Omega$ is denoted by $\mathbf{n}$.

**Remark II.1.** (Constraint for the well-posedness). In the problem (1), we need to ensure $\langle u \rangle_{\Omega} := |\Omega|^{-1} \int_{\Omega} u \, dx = 0$ where $|\Omega|$ is the measure associated with $\Omega$ (see [5] for system with partial coercivity). For the rest, we introduce the space $H_*(\Omega) := \{v \text{ such that } v \in H^1(\Omega), \langle v \rangle_{\Omega} = 0\}$.

### B. Friedrichs' formalism

In this part, we show that our problem (1) is naturally handled by the theory developed by Ern and Guermond (see [3], [4] and [5]).

We start by denoting $v := (\mathbf{j}; u) \in V$ with $V := H(\mathrm{div}, \Omega) \times H_*(\Omega)$ (we call it the graph space) and $L := L^2(\Omega; \mathbb{C}^4)$. We define also two orthogonal projectors $P_j, P_u \in \mathbb{C}^{4 \times 4}$ such that $P_j(v) := (\mathbf{j}, 0), P_u(v) := (0, u)$.

We can rewrite the system (1) as $A(v) = f$ with the linear operator $A: V \to L$ and $f \in L$. We have $A(v) := K(P_j(v)) + D(v)$ with $K(P_j(v)) := \rho P_j(v) + i\omega\mu g * P_j(v)$ and $D(v) := \sum_{k \in \{1, 2, 3\}} D_k \partial_k v$ such that

$$ D_k := \begin{bmatrix} \mathbb{O}_{3,3} & e_k \\ (e_k)^T & 0 \end{bmatrix}, \quad \forall k \in \{1, 2, 3\} \quad (3) $$

where $\mathbb{O}_{3,3}$ is the zero matrix of $\mathbb{R}^3$ and $e_k \in \mathbb{R}^3$ the $k$-th vector of the cartesian basis. The matrix fields $D_k$ are constants so $D$ is bounded in $L$ and such that $D_k^T = D_k$ for $k \in \{1, 2, 3\}$. The operator $K$ is also bounded in $L$ because $G(v) := g * v$ is continuous from $L^2(\Omega; \mathbb{C})$ into $L^2(\Omega; \mathbb{C})$ and by assumptions on $\rho, \mu$.

We introduce the formal adjoint of $A$ as $\tilde{A}(v) := K^*(P_j(v)) - D(v)$ where $K^*(v) := \rho v - i\omega\mu g * v$ is the adjoint of $K$.

From the problem (1) we define the bilinear form $a: V \times L \to \mathbb{C}$ and the linear form $l: L \to \mathbb{C}$ respectively such that $a(v, w) := (A(v), w)_L$ and $l(w) := (f, w)_L$.

In order to give a meaning to traces on the boundary for functions in $V$, we define the operator $N: V \to V'$ such that $\langle N(v), w \rangle_{V',V} := (A(v), w)_L - (v, \tilde{A}(w))_L$. Then the operator $N$ can be expressed in term of the traces on $H(\mathrm{div}, \Omega)$ and $H^1(\Omega)$ (see [4]). In the following we assume enough regularity on the solution so we can give a meaning
to $N$ in $L^2(\partial\Omega) := L^2(\partial\Omega; \mathbb{C}^{4\times4})$. Thus, we define the fields $\mathcal{N} \in L^\infty(\Omega; \mathbb{C}^{4\times4})$ such that $<N(v), w >_{V',V} = (Nv, w)_{L(\partial\Omega)}$. With this settings we have the coercivity up to $\mathcal{N}$. To fix this (see [3]), we define $\mathcal{M} \in L^\infty(\Omega; \mathbb{C}^{4\times4})$ where the new operator $\mathcal{M} - \mathcal{N}$ respects the boundary conditions (there are others assumptions, see [3]). In particular, the following matrix fields are suitable for our problem:

$$ \mathcal{N} := \begin{bmatrix} \mathbb{O}_{3,3} & \mathbf{n} \\ \mathbf{n}^T & 0 \end{bmatrix}, \quad \mathcal{M} := \begin{bmatrix} \mathbb{O}_{3,3} & \mathbf{n} \\ -\mathbf{n}^T & 0 \end{bmatrix} \quad (4) $$

This leads to the following weak problem (With $V_0 := \{v$ such that $v \in V, v \in \ker(\mathcal{M} - \mathcal{N})\}$):

$$ \left\{ \begin{array}{l} \text{Find } \xi \in V_0 \text{ such that} \\ a(\xi, v) = l(w), \quad \forall w \in L \end{array} \right. \quad (5) $$

**Theorem II.2.** (Well-posedness). The problem (5) is well-posed.

*Proof.* The framework for this proof is the Banach-Nečas-Babuška theorem. Note that in our case, we have only the partial coercivity on $\mathcal{P}_j$. So the Poincaré-Steklov inequality is needed to have a bound on $\mathcal{P}_u$ (see [5] for the proof). $\square$

It is also possible to show that if we take the bilinear form $a^{bc}(v, w) := (A(v), w)_L + \frac{1}{2}((\mathcal{M} - \mathcal{N})v, w)_{L(\partial\Omega)}$, we impose weakly the boundary conditions and we solve the problem (1).

### C. Discretization in Discontinuous Galerkin

We assume that $\Omega$ is discretized by a shape regular affine mesh denoted by $\mathcal{Q}_h$ with $h > 0$ a refinement parameter (more general mesh can be considered). We denote by $\mathcal{F}_h$ the set of all the faces of the mesh such that $\mathcal{F}_h = \mathcal{F}_h^\partial \cup \mathcal{F}_h^i$ with $\mathcal{F}_h^\partial$, $\mathcal{F}_h^i$ respectively the sets of all boundary and interface faces. For every $F \in \mathcal{F}_h$, we define the oriented normal vector as $\mathbf{n}_F$.

We introduce the broken finite element space $P_k^b := \{v_h \in L^\infty(\Omega; \mathbb{C}^4) \text{ such that } v_{h|Q} \in \mathbb{P}_Q^k, \forall Q \in \mathcal{Q}_h\}$ with $\mathbb{P}_Q^k$ a basis of polynomials of degree $k \ge 0$ on $Q \in \mathcal{Q}_h$. Our discrete space is defined as $V_h := P_p^b \times P_q^b \not\subset V$, with $p, q \ge 0$. The jump and average operators are defined respectively by $[[v]]_F := v|_{\mathcal{Q}_p} - v|_{\mathcal{Q}_q}$ and $\{v\}_F := \frac{1}{2}(v|_{\mathcal{Q}_p} + v|_{\mathcal{Q}_q})$ with $F = \partial\mathcal{Q}_p \cap \partial\mathcal{Q}_q \in \mathcal{F}_h^i$.

To design a suitable discrete bilinear form, we define a discrete flux by reproducing the behaviour of our problem on an element $Q \in \mathcal{Q}_h$. Then, we sum over each mesh cells and we obtain the discrete bilinear form

$$ a_h^{bc}(v_h, w_h) := (A_h(v_h), w_h)_L + \frac{1}{2}((\mathcal{M} - \mathcal{N})v_h, w_h)_{L(\partial\Omega)} - \sum_{F \in \mathcal{F}_h^i} (\mathcal{N}|_F [[v_h]]_F, \{w_h\}_F)_{L(F)} \quad (6) $$

In order to improve the error analysis (see [3]) and to better ensure the boundary conditions, we add the penalization on the interface $\gamma_h(v_h, w_h) := \sum_{F \in \mathcal{F}_h^i} (S^i|_F [[v_h]]_F, [[w_h]]_F)_{L(F)}$ and we define the new operator $\mathcal{M}^\partial := \mathcal{M} + S^\partial$ with the penalty fields $S^i, S^\partial \in L^\infty(\Omega; \mathbb{C}^{4\times4})$ (see [3] for the assumptions on

<figure><img src="image_2.png" /><figcaption>Fig. 1. Comparison between MIPSE and the Discontinuous Galerkin discretization for a regular mesh with 400 quadrangles. The parameters $\mu$ and $\rho$ are respectively the permeability of the air and the resistivity of the copper.</figcaption></figure>

$S^i$ and $S^\partial$.

In our case, these two penalty fields can be chosen as

$$ S^\partial = \begin{bmatrix} \eta_1 \mathbf{n} \otimes \mathbf{n} & \mathbb{O}_{d,1} \\ \mathbb{O}_{d,1}^T & 0 \end{bmatrix}, \quad S^i|_F = \begin{bmatrix} \eta_2 \mathbf{n}_F \otimes \mathbf{n}_F & \mathbb{O}_{d,1} \\ \mathbb{O}_{d,1}^T & \eta_3 \end{bmatrix} \quad (7) $$

where $\eta_1, \eta_2, \eta_3 > 0$ are user dependant constants which can depend to the size of mesh cells.

This leads to the discrete problem (with $a_h^{stb}(v_h, w_h) := a_h^{bc}(v_h, w_h) + \gamma_h(v_h, w_h)$):

$$ \left\{ \begin{array}{l} \text{Find } \xi_h \in V_h \text{ such that} \\ a_h^{stb}(\xi_h, w_h) = l(w_h), \quad \forall w_h \in V_h \end{array} \right. \quad (8) $$

**Remark II.3.** The problem (8) is consistent by construction.

**Theorem II.4.** (Well-posedness). The discrete problem (8) is well-posed.

*Proof.* It consists to prove the inf-sup conditions with a suitable norm $||\cdot||_{V_h}$ (see [3]). $\square$

### III. NUMERICAL EXPERIMENTS

Let $\Omega := [0, 3.10^{-2}]^2 \times [0, 35.10^{-6}]$ (in meters), $V_h = P_1^b \times P_1^b$ with Lagrange polynomials and $\eta_1 = \eta_3 := 10$, $\eta_2 := 100$. Then, we compare our results to the conforming approach [1] in (Fig.1). We obtain less than 1% of error on the computation of Joule losses for all frequencies considered.

### ACKNOWLEDGMENT

This work is supported by the ANR TensorVim project 22-CE50-0028.

### REFERENCES

[1] Nguyen, T. T., Meunier, G., Guichon, J. M., Chadebec, O., Nguyen, T. S. (2014). An integral formulation for the computation of 3-D eddy current using facet elements. IEEE transactions on magnetics, 50(2), 549-552.

[2] Friedrichs, K. O. (1958). Symmetric positive linear differential equations. Communications on Pure and Applied Mathematics, 11(3), 333-418.

[3] Ern, A., Guermond, J. L. (2006). Discontinuous Galerkin methods for Friedrichs' systems. I. General theory. SIAM journal on numerical analysis, 44(2), 753-778.

[4] Ern, A., Guermond, J. L. (2006). Discontinuous Galerkin methods for Friedrichs' systems. Part II. second-order elliptic PDEs. SIAM journal on numerical analysis, 44(6), 2363-2388.

[5] Ern, A., Guermond, J. L. (2008). Discontinuous Galerkin methods for Friedrichs' systems. Part III. Multifield theories with partial coercivity. SIAM journal on numerical analysis, 46(2), 776-804.

[6] Oseledets, I. V. (2011). Tensor-train decomposition. SIAM Journal on Scientific Computing, 33(5), 2295-2317.