A Fully-Discrete Finite Element Approximation for the Eddy Currents Problem

The eddy current model is obtained from Maxwell’s equations by neglecting the displacement currents in the Ampère-Maxwell’s law. The so-called “A, V − A potential formulation” is nowadays one of the most accepted formulations to solve the eddy current equations numerically, and Bíró & Valli have recently provided its well-posedness and convergence analysis for the time-harmonic eddy current problem. The aim of this paper is to extend the analysis performed by Bíró & Valli to the general transient eddy current model. We provide a backward-Euler fully-discrete approximation based on nodal finite elements and we show that the resulting discrete variational problem is well posed. Furthermore, error estimates that prove optimal convergence are settled.


Introduction
In applications related to electrical power engineering (see for instance [1]) the displacement currents existing in a metallic conductor are negligible compared to the conduction current.In such situations the displacement currents can be dropped from Maxwell equations to obtain a magneto-quasistatic submodel usually called eddy current problem; see for instance [2,Chapter 8].From the mathematical point of view, this

|112
Ingeniería y Ciencia submodel provides a reasonable approximation to the solution of the full Maxwell system in the low frequency range [3].
Among the numerical methods used to approximate eddy current equations, the finite element method (FEM) and methods combining FEM and boundary element method (FEM-BEM) are the most extended, see, for instance, the recent book by Alonso & Valli [4] for a survey on this subject including a large list of references.In the applied mathematical literature, we can find several recent papers devoted to the numerical analysis of the 3D time dependent eddy current model, in bounded domains as well as in unbounded domains by using FEM and FEM-BEM methods: Meddahi & Selgas [5], Ma [6], Acevedo et al. [7], Kang & Kim [8], Prato et al. [9], Acevedo & Meddahi [10], Bermudez et al. [11], Camaño & Rodríguez [12].The main differences among these works are the selected unknowns to compute the electromagnetic field, the topological assumptions on the conductor domain and the boundary conditions when the problem is solved in a bounded domain.
The aim of this paper is to analyze a finite element fully-discrete approximation for the time-dependent eddy current problem in a bounded domain, based in a formulation in terms of a vector magnetic potential and an electric scalar potential.Numerical experiments showing the efficiency of this approach, were reported by Bíró & Preis (1989) [13], and nowadays, it is the basis of several commercial codes to compute the solution of the eddy current model.Bíró & Valli (2007) [14] have studied that potential formulation for the time-harmonic eddy current model and have proved its well-posedness and theoretical convergence in a general geometric situation.However, a similar analysis for the transient case has not been realized yet.Although, Kang & Kim (2009) in [8] have recently provided quasi-optimal error estimates, showing the theoretical convergence of the method in the case of a simply-connected conductor and by assuming homogeneous conditions for the electromagnetic fields on the boundary of the conductor, these assumptions are very restrictive for many real applications (e.g., metallurgical electrodes [15] and power transformers [16]).Moreover, the decay conditions on the fields at infinity (see [3]) allow to assume that the electromagnetic field is weak far away from the conductor and not on its boundary, as it was assumed by Kang & Kim when the artificial homogeneous boundary conditions are supposed.
In order to improve the results obtained by Kang & King for more realistic applications, we consider a multiply-connected conductor and we opt for a typical approach: to restrict the eddy current equations to a sufficiently large computational domain containing the region of interest and impose a convenient artificial homogeneous boundary conditions for the electromagnetic fields on its border.We provide a backward-Euler fully-discrete approximation based on nodal finite elements to approximate the solution of the resultant model.This fully-discrete scheme solves in each time an elliptic problem and we use the techniques used in [14] and [17] to prove the ellipticity of its bilinear form.Furthermore, by using this ellipticity we define projection operators to the discrete finite element subspaces and obtain quasi-optimal error estimates, which allows us to approximate the typical physical variables of interest of the eddy current problem: the eddy currents in the conductor and the magnetic induction in the computational domain.
The outline of the paper is as follows: In section 2, we summarize some results concerning tangential traces in H(curl; Ω) and recall some basic results for the study of time-dependent problems.In Section 3, we introduce the eddy current model in a bounded domain and deduce a potential-based formulation (the so-called "A, V − A potential formulation") for the time-dependent eddy current problem.In Section 4, we obtain a variational formulation for the problem and its fully-discrete approximation scheme is analyzed in Section 5. Finally, the results presented in Section 6 prove that the resulting fully discrete scheme is convergent in an optimal way.We end this paper by summarizing its main results in Section 7.
We also recall that Hilbert space and that C ∞ (Ω) 3 is dense in H(curl; Ω) (see, for instance, [18,Theorem 3.26]).Tangential traces of functions in H(curl; Ω) are also understood even in the case of polyhedral domains thanks to the recent results given by Buffa & Ciarlet [19,20] and Buffa, Costabel & Sheen [21].We give here a brief summary of these fundamental tools.We begin by considering the space endowed with the standard norm in L 2 (Γ) 3 .We define the tangential trace γ τ : The previous traces can be extended by continuity to H 1 (Ω) 3 .
Since in the eddy current problem the physical domain requires to be split into two non-overlapping domains: the conductor and the insulator domain (see Section 3.1 below), we end this subsection by recalling the following classical result that will be used often in the rest of the paper.

Basic spaces for time dependent problems
Since we will deal with a time-domain problem, besides the Sobolev spaces defined above, we need to introduce spaces of functions defined on a bounded time interval (0, T ) and with values in a separable Hilbert space V , whose norm is denoted here by ∥ • ∥ V .We use the notation C 0 ([0, T ]; V ) for the Banach space consisting of all continuous functions f : We also consider the space L 2 (0, T ; V ) of classes of functions f : (0, T ) → V that are Böchner-measurable and such that

|118
Ingeniería y Ciencia Furthermore, we will use the space where d dt f is the (generalized) time derivative of f ; see, for instance, [22,Section 23.5].

Eddy current problem
We consider a standard eddy current problem: to determine the electromagnetic fields induced in a three-dimensional conducting domain by a given source time-dependent compactly-supported current density.The eddy current problem is in principle posed in the whole space.However, we restrict it to a bounded computational domain containing both, the conductor and the support of the source current, such that adequate boundary conditions can be imposed on its boundary.To this aim, we choose the geometry of the computational domain as simple as possible (e.g., simply connected with a connected boundary).
Let Ω c ⊂ R 3 be the conducting domain and let us assume that it is an open and bounded set with boundary Γ c .Let Ω ⊂ R 3 be a simply connected bounded domain with a connected boundary Γ, such that Ω c ⊂ Ω.We suppose that both, Ω and Ω c are Lipschitz domains and we denote by n and n c the outward unit normal vectors to Ω and Ω c , respectively.Let Ω d := Ω \ Ω c be the subdomain of Ω occupied by dielectric material, which includes the support of the source current J d (see Figure 1).

The electric and magnetic fields
solutions of a submodel of Maxwell's equations obtained by neglecting the displacement currents (see, for instance, [13]): Let us remark that because of H(•, t) must belong to H(curl; Ω) a.e. in [0, T ], Lemma 1 implies the magnetic field has to satisfy the coupling conditions, Equation (8).
The magnetic permeability µ and the conductivity σ are bounded functions satisfying: The data of the problem are the source current density J d ∈ L 2 (0, T ; L 2 (Ω) 3 ), for which we assume for a.e.t ∈ (0, T ) and the initial magnetic field H 0 ∈ H(curl; Ω).

The A, V − A potential formulation
We now recall a classical formulation of the eddy current problem in terms of two potentials: a magnetic vector potential A and an electric scalar potential V .We refer to Bíró & Preis [13] for a detailed discussion, which also includes numerical tests showing the efficiency of this approach.

If we define
we obtain Hence, from (3) there follows that div Then, from (3), (5) and recalling that suppJ Therefore, the Equation (17) implies that On the other hand, using the Equation (3) together with the Equations (10) and (17), there follows that σ Moreover, from the Equations ( 5) and ( 10), we have curl The Equations ( 11)-( 15) together with the Equations ( 18)-( 21) will be collected in the potential formulation.Hence, we are led to the following formulation of problem (3)-( 9) in terms of the potentials A : div curl ) and satisfying Remark 3.1.Equation ( 33) is obtained from the Equations ( 7), ( 9) and (10), and the Equation ( 32) is a consequence of (16).Furthermore, the initial condition A 0 in (31) must satisfy and since Ω is a simply connected set, it can be characterized as the ing.cienc., vol.9, no.17, pp.111-145, enero-junio.2013.

123|
unique solution of the boundary value problem (see [23, Theorem I.3.5.]): Moreover, A 0 can be computed by using a finite element approximation of the following mixed problem [24, Section 4]: find where In order to prove the well-posedness of this last problem, we can use the well-known Babuska-Brezzi theory for mixed problems (see for instance [25] and the references given there).Furthermore, it follows easily that the Lagrange multiplier λ vanishes identically.
Remark 3.2.The differential constraint ( 25) is called the Coulomb gauge condition and it is necessary to assure the uniqueness of potential A [13].
The Coulomb gauge condition is not easy to treat at the discrete level, because it is not simple to construct a suitable space of finite elements which are divergence-free.Here we will follow the ideas from [14], by using the addition of a div − div term to our variational formulation (see Section 4 bellow), which is equivalent to add a penalization term in the Ampère law (see, for instance, [13,14]).

Variational formulation
The aim of this section is to give a variational formulation of Problem, Equations( 22)- (33).To do this, we introduce the following functional spaces: where Ω 1 c , . . ., Ω m c are the connected components of Ω c .The spaces X and H 1 ♯ (Ω c ) are endowed respectively with the norms and Let Z ∈ X .By multiplying the Equations ( 22) and ( 24) by Z, integrating in Ω c and Ω d respectively and summing the resulting equations, we deduce that We notice now that from ( 22), ( 24), ( 28) and Lemma 1 it follows that Then, by integrating by parts (see (2)) and using (33), we obtain Hence, from (37) we have ing.cienc., vol.9, no.17, pp.111-145, enero-junio.2013.
On the other hand, by integrating by parts (1) and using ( 23) and ( 27) Summing the last two equations, we obtain the following variational formulation of Problem ( 22)- (33): ) and satisfying the initial condition where (u, w) σ := Remark 4.1.To our knowledge, the well-posedness of Problem (38)-(39) has not been proved yet.We have tried to apply the theoretical framework for parabolic problems arising from electromagnetism (see Zlámal [27,28]) and for classical degenerate parabolic problems (see Showalter [29, Chapter III]), without satisfactory results in both cases.However, by adapting the techniques from [27,28], it is easily seen that Problem (38)-( 39) has at most one solution, but this is not the case for the existence result.The main difficulty is that the bilinear form It is straightforward to verify that a solution of Problem (38)-( 39) is actually a solution of the strong form of the problem given by Equations ( 22)- (33).In particular, we can obtain the gauge condition div A = 0 as follows: for any t ∈ [0, T ], let ξ ∈ H 1 (Ω) be a weak solution of the compatible Neumann problem ∆ξ = div A(t) in Ω, ∂ξ/∂n = 0 on Γ.Hence, by testing (38) with Z := ∇ξ and u : 3. The idea of using an average µ * to impose the Coulomb gauge condition in the weak formulation is taken from [14].This is necessary because a finite element approximation based on a weak form, in which the term ∫ Ω µ −1 div Z div W is present, can be inefficient if the coefficient µ has jumps (see [30,Section 5.7.4]).

A fully discrete scheme
In what follows we assume that Ω and Ω c are Lipschitz polyhedra (we recall that Ω is simply-connected).Let {T h } h be a regular family of tetrahedral meshes of Ω such that each element K ∈ T h is contained either in Ω c or in Ω d .As usual, h stands for the largest diameter of tetrahedra K in T h .
We consider a uniform partition {t n := n∆t : n = 0, . . ., N } of [0, T ] with a step size ∆t := T N .For any finite sequence The fully-discrete version of Problem ( 38)-(39) reads as follows: for any (Z, u) ∈ X h × M h and satisfying where A 0,h ∈ X h is a suitable approximation of A 0 to obtain optimal error estimates (see Corollary 6.1 below).
In order to prove that Problem (40)-( 41) has a unique solution, we first notice that at each iteration step we need to find where Hence, the existence and uniqueness of solution of Problem (40)-(41) follows by combining the Lax-Milgram's Lemma and the following ellipticity result.
Lemma 2. There exists a constant C > 0 such that for any Z ∈ X and u ∈ H 1 ♯ (Ω c ).

Error estimates
In this section we will prove error estimates for our fully-discrete scheme.
To this end, we have to assume that Problem (38)-( 39) has a unique solution and consider the elliptic projection operators P h : X → X h and It is a simple matter to see that Lax-Milgram's Lemma and Lemma 2 imply that P h and Q h are well defined.Moreover, the well-known Cea's Lemma yields that there exist positive constants C 1 and C 2 , independent of h, such that From now on, let us introduce the following notations:

Ingeniería y Ciencia
Lemma 3.There exists a positive constant C, independent of h and ∆t, such that in the last identity and using the estimates we obtain ing.cienc., vol.9, no.17, pp.111-145, enero-junio.2013.
Then, by using the Cauchy-Schwartz inequality in the right hand side, it follows that 1 ∆t for k = 1, . . ., n.In particular, for k = 1, . . ., n.Then, summing over k, we get Hence, ] .
Let us now take ]

|134
Ingeniería y Ciencia and consequently, using the Cauchy-Schwarz inequality on the right hand side, we deduce that ] .
Hence, using the inequality ] on the left hand side, leads to for k = 1, . . ., n. Summing over k, we have ] .

|136
Ingeniería y Ciencia there exists a constant C > 0, independent of h and ∆t, such that ] dt ing.cienc., vol.9, no.17, pp.111-145, enero-junio.2013. Hence, The regularity assumptions on A and v allow us commute the time derivative with P h and Q h , i.e,

∂ t (P h A(t))
Moreover, if we define ρ h 1 (t) := A(t) − P h A(t) and then from (45) it follows that On the other hand, the Cauchy-Schwarz inequality shows that

Proof. Let Π
) → M h be the standard scalar finite element Lagrange interpolant.The result follows by using previous theorem, the regularity assumptions on A and v, and the well-known approximation properties of Π h and Π h (see, for instance, [32]).Remark 6.2.Concerning this convergence result, it has to be noted that the spatial regularity of A (and in particular the regularity of A 0 ) is not ensured if Ω has reentrant corners or edges, namely, if it is a non-convex polyhedron (see [33], [34]).More important, in that case the space H 1 (Ω) 3 ∩ H(div; Ω) turns out to be a proper closed subspace of X and hence the nodal finite element approximate solution cannot approach an exact solution A ∈ L 2 (0, T ; X ) with A / ∈ L 2 (0, T ; H 1 (Ω) 3 ∩ H(div; Ω)), and convergence in X is lost.
However, the result we have proved here above ensures that the nodal finite element approximation is convergent either if the solution is regular (and this information could be available even for a non-convex polyhedron Ω) or if Ω is a convex polyhedron.Actually, in [17,Section 5] is underlined that if Ω is a convex polyhedron and µH ∈ H p (Ω) 3 with 0 < p ≤ 1, the vector potential A (which is characterized by ( 10)-( 12)) belongs to H 1+s (Ω) 3 for some 0 < s ≤ 1.Let us also note the assumption that Ω is convex is not a severe restriction, as in most real-life applications ∂Ω arises from a somehow arbitrary truncation of the whole space and hence, reentrant corners and edges of Ω can be easily avoided.
A cure for the lack of convergence of nodal finite element approximation for Maxwell equations in the presence of reentrant corners and edges has been proposed by Costabel and Dauge in [35], where they introduce a special weight in the penalization term, thus making it possible to use standard nodal finite elements in a numerically efficient way.Another alternative for the eddy current problem has been reported by Bíró [36]: edge elements are employed for the approximation of the potential A, without requiring that the Coulomb gauge condition be satisfied.However, a complete theory assuring the effectiveness of this last approach, even for the harmonic case, is not available.
To end this section, let us notice that we can approximate at each time t n the physical quantities of interest: the eddy currents σE(t n )

|140
Ingeniería y Ciencia in the conductor Ω c and the magnetic induction field µH(t n ) in the computational domain Ω.In fact, the identities ( 17) and (10)

Conclusions
We have proposed a fully-discrete approximation for the so-called "A, V − A potential formulation" of a transient eddy current problem in a bounded ing.cienc., vol.9, no.17, pp.111-145, enero-junio.2013.
domain, without topological restrictions on the conductor.The fullydiscrete scheme is obtained by a finite element space discretization of standard nodal finite elements and a time discretization based on a backward-Euler implicit scheme.
We have shown that at each step the fully-discrete scheme solves an elliptic problem and consequently the well-posedness of the discrete problem follows by using the well-known Lax-Milgram's Lemma.Furthermore, under usual regularity assumptions on the solution of the continuous problem, we obtain quasi-optimal error estimates, which allows us to approximate the typical physical variables of interest of the eddy current problem: the eddy currents in the conductor and the magnetic induction in the computational domain.

Figure 1 :
Figure 1: The geometrical setting of problem.In this example the conductor Ω c has two connected components Ω 1 c and Ω 2 c , Ω is the computational domain (the box) which contains the conductor and support of J d , and Ω d = Ω \ Ω c is the dielectric.
Finally, to get the error estimate for the magnetic induction approximation, we notice thatµH(t n ) − curl A n h = curl A(t n ) − curl A n h = curl e n 1for n