Mathematical and numerical analysis of a transient magnetic model with voltage drop excitations

Mathematical and Numerical Analysis of a Transient Magnetic Model with Voltage Drop Excitations A. Bermúdeza,b,c, M. Piñeiroa, P. Salgadoa,c,∗ aDepartamento de Matemática Aplicada, Facultade de Matemáticas, Universidade de Santiago de Compostela, E-15782 Santiago de Compostela, Spain bInstituto de Matemáticas (IMAT), Facultade de Matemáticas, E-15782 Santiago de Compostela, Spain cInstituto Tecnológico de Matemática Industrial (ITMATI), Edif. Instituto Investigaciones Tecnológicas, planta -1, E-15782 Santiago de Compostela, Spain Abstract This paper deals with the mathematical and numerical analysis of a nonlinear 2D transient magnetic model when the source data are given in terms of the voltage drop excitations in conductors and the remanent magnetic flux for permanent magnets. The formulation consists of a distributed nonlinear magnetostatic model with time appearing as a parameter, and a circuit equation linking currents and voltage drops. This last equation is used to express the problem as an implicit ODE system whose operator involves the resolution of the distributed model. The model is spatially discretized using a finite element method and an implicit Euler scheme is employed for time discretization. We perform the mathematical analysis of the problem at both the continuous and discrete levels and obtain an error estimate that is illustrated with some numerical results. Keywords: Transient magnetic, Nonlinear partial differential equation, Finite element approximation, Voltage drops 1. Introduction The objective of this paper is the mathematical and numerical analysis of a nonlinear transient magnetic model defined in a two-dimensional domain, with sources given in terms of the potential drops in conductors and the remanent fluxes of permanent magnets. This model arises, for instance, in the simulation of electric machines and, in particular, of permanent magnet synchronous motors (PMSM). In this kind of devices, the magnetic core is usually laminated orthogonally to the direction of the currents traversing the coils. Moreover, eddy current losses are often neglected in permanent magnets, so that these regions are modelled as non-conducting; eventually, a posteriori formulas could be used to estimate such losses (see, for instance, [17, 21]). Both of the above simplifications allow us to build a 2D transient magnetic model in a cross section of the device, the stator coils being the only conducting part. These coils are generally composed by stranded wires carrying a uniformly distributed current density. The mathematical model used to simulate these conductors strongly depends on the kind of the imposed source; see, for instance, [5]. Indeed, if the source data are given in terms of the current traversing the wires, the problem reduces to solving a nonlinear magnetostatic problem at each time step, and thus time appears as a parameter. However, in the case where the potential drops are given, the distributed magnetostatic model has to be coupled with a circuit equation linking currents and voltage drops. In this paper, we focus on this last case because the model offers challenges from a mathematical and numerical point of view, as detailed below. ∗Corresponding author Email address: [email protected] (P. Salgado) Here, we give a first step towards the analysis of the genuine physical problem, as we do not consider the motion of the machine, what would lead to a much more difficult problem; see, for instance, [8] for a case incorporating also induction effects. Our mathematical model will be obtained from the low-frequency approximation of Maxwell equations, without taking any eddy current effects into consideration. Therefore, we will deal with an integro-differential problem coupling an elliptic partial differential equation, written in terms of the magnetic vector potential, with the circuit equations relating currents and voltage drops in stranded conductors. The partial differential equations are nonlinear due to the presence of ferromagnetic materials in the cores which usually have a strongly nonlinear magnetic behavior. In the literature, we can find several references dealing with the analysis of low-frequency electromagnetic models coupled with circuit equations. For example, in [15], the authors study the well-posedness of a threedimensional field/circuit nonlinear problem in the presence of eddy currents and provide error estimates for time discretization. In [11], the authors deal with a 3D field/circuit linear model, focusing only on the continuous formulation. Alternatively, field/circuit models also fit in the framework of differential algebraic systems of equations (DAE), usually when using finite integration techniques for the spatial discretization; see, for instance, [4, 3]. Finally, we also highlight the results presented in [10], where we can find a study of some classes of differential algebraic systems of equations in an abstract framework, in particular covering the case of systems of DAE coupled with partial differential equations (PDE). However, we deal with a system of elliptic partial differential equations coupled with a vector ordinary differential equation in terms of time that is not covered by the previous results. As discussed above, we focus on a model that does not consider eddy current effects. Firstly, we obtain an integro-differential problem arising from the coupling of the Maxwell system of equations with the circuit equations relating currents and voltage drops in stranded conductors. To perform its mathematical analysis, it is written as a nonlinear system of implicit ordinary differential equations in terms of the currents traversing the coils, which are functions of time. The operator defining this system expresses the so-called flux linkages per unit length in the coils in terms of the currents traversing them, via the resolution of some 2D magnetostatics problems. The properties of this operator are deduced directly from results already existing in the literature (specifically, those appearing in [14]). To perform the numerical approximation of the continuous problem we propose an Euler-implicit scheme for the ODE, combined with a finite element method for the approximation of the involved distributed operator. Some convergence results are obtained for this numerical scheme. However, for the numerical implementation, we use the alternative approach proposed in [5] which consists in eliminating the unknown currents from the system by means of the circuit equations. This idea is also exploited in the theoretical analysis of the eddy current model performed in [11]. As a consequence, we need to prove an equivalence result between the implemented scheme and the discrete problem theoretically analysed. The paper is organized as follows. In Section 2 we present the 2D nonlinear transient magnetic model in a cross-section transversal to the device, written in terms of the magnetic vector potential; moreover, we express the problem as a system of implicit ODE, and perform its mathematical analysis in the continuous case. In Section 3 we introduce the finite element discretization of the magnetostatics problem involved in the definition of the ODE operator. In Section 4 we propose an implicit Euler scheme for the discretization of the system of ODE and prove an error estimate for its solution. In Section 5, we show some numerical results for a test with analytical solution to illustrate the obtained convergence results. Finally, in appendix Appendix A we prove the equivalence between the analysed problem and the implemented one; appendix Appendix B contains the analytical expression of the magnetic vector potential corresponding to the analytical test used for the numerical results. 2. Mathematical Analysis of the Continuous Problem In this section we state a 2D transient magnetic problem that arises in the mathematical modeling of laminated magnetic media, with sources given in terms of the potential drops per unit length in conductors and the remanent fluxes in permanent magnets. A similar model, without permanent magnets, has been studied in [5] from a computational point of view. Here, we will also present and analyze the continuous formulation. 2 Figure 1: Illustration of the subdomains. 2.1. A two-dimensional transient magnetic model Let us assume that the current density sources Jhave non-null component only in the zspace direction and that this component does not depend on z, i.e., J=Jzez, with Jz=Jz(x, y, t). We also assume that the geometry and the magnetic field Hare invariant along the z−direction, and that all materials are magnetically isotropic. In this case, under an appropriate decay of fields at infinity (see [16]), the magnetic field H, and then the magnetic induction B, have only components on the xy−plane and both are independent of z. Since we are interested in using a finite element method for the numerical solution, we will restrict ourselves to a bounded domain. Thus, let us consider a 2D convex bounded domain Ω, with Lipschitz continuous boundary, containing a cross-section transversal to the device. For a given current density J∈L2(Ω)3, we seek H∈H(curl,Ω) and B∈H(div,Ω) such that: curl H=Jin Ω,(2.1) div B= 0 in Ω,(2.2) B·n= 0 on ∂Ω,(2.3) where the boundary condition means there is no magnetic flux through the boundary. This model is completed below with the constitutive law relating the magnetic field to the flux density. Let us assume that Ωis composed of the following open subsets (see fig. 1) • a magnetically linear subdomain Ω0, • non-magnetic connected conductors Ωn,n= 1, . . . , Nc, • a permanent magnet region Ωpm and • a nonlinear ferromagnetic core Ωnl. We further assume that the boundaries of the conductors, ∂Ωn,n= 1, . . . , Nc, are mutually disjoint and do not touch the boundary of Ω, and also that the same is true for the boundaries of the connected components of the permanent magnet region. We notice that all parts of the domain are non-conducting except for the non-magnetic conductors, which support the current density J. Moreover, we will use the notation Ωc:= ∪Nc n=1Ωnand U:= {Ω0,Ωc,Ωpm,Ωnl}. In this framework, vector fields Hand Bare linked by the constitutive relations: H=ν0Bin Ω0∪Ωc, H=νpmB−νpmBrin Ωpm, H=eν(|B|)Bin Ωnl, 3 where ν0is the vacuum magnetic reluctivity, Bris the remanent flux density in the permanent magnets and νpm : Ωpm −→ R+is the magnetic reluctivity in the permanent magnets. In principle, both the magnetically linear subdomain and the nonlinear core subdomain may have several parts with different magnetic reluctivities. However, for the sake of simplicity, we have assumed there is only one for each of these two subdomains and that the magnetic reluctivity of the linear one is that of the vacuum, ν0. We assume that Brhas only components on the xy−plane and both are independent of the z−coordinate and that νpm ∈L∞(Ωpm), also being uniformly bounded from below by a positive constant. Furthermore, we define the global magnetic reluctivity function ν: Ω ×R+ 0−→ R+as ν(x;s) :=    ν0if x∈Ω0∪Ωc, νpm(x)if x∈Ωpm, eν(s)if x∈Ωnl. Let us make the following assumptions on the nonlinear reluctivity eν:R+ 0−→ R+, ∃ν1, ν2>0 : ν1≤eν(s)≤ν2a.e. in R+ 0,(2.4) ∃Meν>0 : |eν(p)p−eν(q)q| ≤ Meν|p−q| ∀p, q ∈R+ 0,(2.5) ∃αeν>0:(eν(p)p−eν(q)q)(p−q)≥αeν|p−q|2∀p, q ∈R+ 0.(2.6) We notice that the above assumptions on the reluctivity can be derived from the natural properties of the physical BH-curves corresponding to nonlinear magnetic materials (see [12, Chapter 2]). Finally, we will suppose that all conductors are stranded, which makes it possible to assume that the current density is uniform and given by Jz,n(t) = in(t) meas(Ωn), n = 1, . . . , Nc, where in(t)denotes the total current across Ωnat time t. Actually, for each conductor, the source can be given in terms of either the current or the potential drop per unit length in the z-direction, and we will focus on the latter alternative. In order to solve the described two-dimensional model, it is convenient to introduce a magnetic vector potential because it leads to solving a scalar problem instead of a vector one. Since Bis divergence-free, there exists a so-called magnetic vector potential Asuch that B=curl A. Under the assumptions above, we can choose a magnetic vector potential that does not depend on zand does not have either xor y components, i.e., A=A(x, y, t)ez(see, for instance, [7]). Next, we will see how to include the potential drops per unit length as sources of our formulation. For the sake of simplicity, we will assume that the electric conductivity σis constant for all conductors, but otherwise the development below can be applied with no significant change (see [5]). Let us denote by E the electric field. From Faraday’s law in the conducting domain, ∂B ∂t +curl E=0in Ωc, the invariance under translation in the z-direction hypothesis and the axial direction of the currents, we deduce that there exist Ncscalar potentials vn,n= 1, . . . , Nc, unique up to a constant, such that ∂A ∂t +E=−grad vnin Ωn×R, n = 1, . . . , Nc. Taking into account the assumptions on Jand Ohm’s law, J=σE, we deduce that Ein conductors has non-null component only in the zspace direction which is spatially constant in each Ωn,n= 1, . . . , Nc. Moreover, since A=Aez, we have −grad vn=−∂vn ∂z ezin each Ωn,n= 1, . . . , Nc. As a consequence, the 4 above equation reduces to ∂A ∂t +Ez=−cn(t)in Ωn, n = 1, . . . , Nc,(2.7) where cn(t) := −∂vn ∂z (t)is the potential drop per unit length in direction zin conductor Ωn,n= 1, . . . , Nc. Multiplying eq. (2.7) by the electric conductivity, integrating on each Ωn,n= 1, . . . , Nc, and taking Ohm’s law into account we deduce d dt ZΩn σA(x, y, t)dxdy +in(t) = −cn(t)σmeas(Ωn), n = 1, . . . , Nc. Thus, if voltage drops are given in conductors, then the problem to be solved is the following: Problem 1. Given c(t)∈RNc,Br(x, y)and a vector of initial currents i0∈RNc, find A(x, y;t)and i(t)∈RNcfor every t∈[0, T]satisfying i(0) = i0, −div(ν0grad A) = 0 in Ω0,(2.8) −div(ν0grad A) = in(t) meas(Ωn)in Ωn, n = 1, . . . , Nc,(2.9) −div(νpmgrad A) = −div νpm (Br)⊥in Ωpm,(2.10) −div(eν(|grad A|)grad A) = 0 in Ωnl,(2.11) [ν(·;|grad A|)grad A·n]Γ=νpm (Br)⊥·npm, 0, if Γ⊂∂Ωpm, otherwise,(2.12) A= 0 on ∂Ω,(2.13) and, for every t∈(0, T], d dt ZΩn σA(x, y, t)dxdy +in(t) = −cn(t)σmeas(Ωn), n = 1, . . . , Nc.(2.14) In the above equations, [·]Γdenotes the jump across any interface Γ,(Br)⊥:= −Br yex+Br xey,nis a unit normal vector to interface Γand npm is a unit normal vector to ∂Ωpm pointing outside Ωpm.We observe that eqs. (2.8) to (2.11) follow from eq. (2.1) and that boundary condition (2.13) implies eq. (2.3). Remark 2. Notice that the jump discontinuity in eq. (2.12) follows from the transmission condition [H×n]Γ=0, which, at the same time, follows directly from the regularity of H, as long as there are no surface currents on Γ. The variational formulation associated to eddy currents problem eqs. (2.8) to (2.13) can be obtained using classical techniques, resulting in Problem 3. Given c(t)∈ C([0, T ])Nc,i0∈RNcand Br∈L2(Ωpm)3, find A(t)∈H1 0(Ω) for every t∈[0, T] and i(t)∈ C0,1([0, T])Ncsatisfying i(0) = i0, ZΩ ν(x;|grad A(x, t)|)grad A(x, t)·grad W(x) = Nc X n=1 ZΩn in(t) meas(Ωn)W(x) + ZΩpm νpm(x) (Br)⊥(x)·grad W(x), 5 for every W∈H1 0(Ω) and t∈[0, T], and d dt ZΩn σA(t) + in(t) = −cn(t)σmeas(Ωn), n = 1, . . . , Ncin (0, T ]. In the next section, we will write Problem 3 as a nonlinear implicit system of ordinary differential equations in order to prove that it is well-posed. 2.2. Transient magnetic problem as a system of ODE Let F:RNc−→ RNcbe the nonlinear operator defined as F(i) := ZΩ1 σA, . . . , ZΩNc σA!T ∈RNc, with Athe solution of the nonlinear magnetostatics problem: Problem 4. Given i∈RNcand Br∈L2(Ωpm)3, find A∈H1 0(Ω) such that ZΩ ν(x;|grad A(x)|)grad A(x)·grad W(x) = Nc X n=1 ZΩn in meas(Ωn)W(x) + ZΩpm νpm(x) (Br)⊥(x)·grad W(x), for every W∈H1 0(Ω). Let us notice that the integrals characterising the components of Fare related to the so-called flux linkages per unit length since the latter are defined, for conductor Ωn,n= 1, . . . , Nc, as 1 σmeas(Ωn)ZΩn σA. Moreover, we notice that we are using the following notation convention: we denote vector fields with uppercase bold letters, vectors in Rnwith lowercase bold letters and vector operators with calligraphic bold letters. Theorem 5. Problem 4 has a unique solution. Proof. The proof of this theorem follows directly from the results presented in [14]. Indeed, let B: H1 0(Ω) −→ H−1(Ω) be the operator defined by ⟨B(A), W ⟩=ZΩ ν(·;|grad A|)grad A·grad W for every W∈H1 0(Ω). Under conditions (2.4) to (2.6), operator Bis strongly monotone and Lipschitz continuous, with constants α=C−2 P F min ν0, ν1 pm, αeνand M= max ν0,||νpm||L∞(Ωpm),3Meν, respectively, CP F being the Poincaré-Friedrichs inequality constant (see [12]). Concerning the right-hand side, since functions in meas(Ωn)χΩnbelong to L2(Ω) for n= 1, . . . , Nc, (χKbeing the indicator function of set K), and Br∈L2(Ωpm)3, then the operator associated to the right-hand side is in H−1(Ω). 6 From this theorem, we deduce that operator Fis well defined and therefore we can rewrite Problem 3 as Problem 6. Given c(t)∈ C([0, T ])Ncand i0∈RNc, find i(t)∈ C0,1([0, T ])Ncsuch that i(0) = i0and d dtF(i(t)) + i(t) = −(c1(t)σmeas(Ω1), . . . , cNc(t)σmeas(ΩNc))Tin (0, T ]. Remark 7. We notice that, due to the definition of operator F, it is obvious that Problems 3 and 6 are equivalent. Theorem 8. Operator Fis strongly monotone and globally Lipschitz continuous in RNcwith respective constants CSM and CLto be defined below. Proof. Let i1,i2∈RNcbe given and A1, A2∈H1 0(Ω) be the associated solutions to Problem 4, respectively. Then, Fnij=RΩnσAjfor j= 1,2and n= 1, . . . , Nc. Let us consider the inner product in RNcdefined as follows k1∗k2:= Nc X n=1 k1 nk2 n σmeas(Ωn), with ∥·∥∗the associated norm. First, we will prove that Fis strongly monotone: ⟨B(A1)−B(A2), A1−A2⟩= Nc X n=1 ZΩn i1 n−i2 n meas(Ωn)(A1−A2) = i1−i2∗Fi1−Fi2. Since Bis strongly monotone with constant α, α||A1−A2||2 H1(Ω) ≤ ⟨B(A1)−B(A2), A1−A2⟩=i1−i2∗Fi1−Fi2.(2.15) Now, taking W∈H1 0(Ω) such that ZΩn σW =i1 n−i2 n, n = 1, . . . , Nc,and ||W||H1(Ω) ≤Ci1−i2∗(2.16) for some C > 0independent of i1,i2, we get ⟨B(A1)−B(A2), W ⟩= Nc X n=1 ZΩn i1 n−i2 n meas(Ωn)W=i1−i22 ∗. Thus, taking into account that Bis Lipschitz continuous, i1−i22 ∗≤M||A1−A2||H1(Ω)||W||H1(Ω) ≤MC||A1−A2||H1(Ω) i1−i2∗. Replacing in eq. (2.15) we get α M2C2i1−i22 ∗≤α||A1−A2||2 H1(Ω) ≤i1−i2∗Fi1−Fi2, and then Fis a strongly monotone operator globally in RNcin the ||·||∗-norm with constant α/M2C2. We define CSM the corresponding constant in the usual norm. We notice that we can take W∈H1 0(Ω) verifying eq. (2.16). Indeed, let e Ω := Ω \Ωc,Cn:= i1 n−i2 n meas(Ωn), 7 n= 1, . . . , Nc, and gn∈H1/2(∂Ωn)with gn(x) = Cnfor every x∈∂Ωn,n= 1, . . ., Nc. Then, ZΩn σCn=i1 n−i2 n. Moreover, let us consider the following Dirichlet problem:            Given gn∈H1/2(∂Ωn),n= 1, . . . , Nc, find f W∈H1 0(e Ω) such that −∆f W= 0 in e Ω, f W=gnon ∂Ωn, n = 1, . . . , Nc. This problem is well-defined and ||f W||H1( e Ω) ≤Ci1−i2∗, with Cindependent of i1,i2. We can define W∈H1 0(Ω) by W:= f Win e Ω, Cnin Ωn, n = 1, . . . , Nc. Now, we will show that Fis Lipschitz continuous. Indeed, since ZΩn σ(A1−A2)2 ≤σ2||A1−A2||2 L1(Ωn)≤σ2meas(Ωn)||A1−A2||2 L2(Ωn)2 ≤σ2meas(Ωn)||A1−A2||2 H1(Ωn), then, for every n= 1, . . . , Nc, Fi1−Fi2∗≤ Nc X n=1 σ||A1−A2||2 H1(Ωn)!1/2 ≤C1||A1−A2||H1(Ω). Finally, since ||A1−A2||2 H1(Ω) ≤1 αi1−i2∗Fi1−Fi2 ≤1 αi1−i2∗Fi1−Fi2∗≤C1 αi1−i2∗||A1−A2||H1(Ω), we conclude that Fi1−Fi2∗≤C2 1 α||i1−i2||∗, and therefore Fis Lipschitz continuous globally in RNcin the ||·||∗-norm with constant C2 1/α. We define CLthe corresponding constant in the usual norm. From the last theorem, applying a result by E. H. Zarantonello (see [19], Theorem 25.B), we deduce that Corollary 9. Fis invertible and its inverse F−1is Lipschitz continuous with Lipschitz constant equal to 1 CSM . This result allows us to rewrite Problem 6 in the following way: 8 Problem 10. Given c(t)∈ C([0, T ])Ncand i0∈RNc, find i(t)∈ C0,1([0, T ])Ncsuch that i(0) = i0and d dtℓ(t) + F−1(ℓ(t)) = −(c1(t)σmeas(Ω1), . . . , cNc(t)σmeas(ΩNc))Tin (0, T ], with ℓ(t) = F(i(t)) for every t∈[0, T]. Remark 11. Problems 6 and 10 can be defined with lower regularity assumptions on the source data c(t). For instance, if c(t)is Lebesgue-measurable in [0, T]and |c(t)|is bounded by a Lebesgue integrable function, both problems have a unique absolutely continuous solution, i(t)∈AC([0, T]), fulfilling the differential equation almost everywhere in [0, T]. Furthermore, most of the results presented in this paper can also be proved under these assumptions. Theorem 12. Problem 6 has a unique solution i(t)∈ C0,1([0, T])Ncsuch that ∥i(t)∥ ≤ ∥i0∥+T CSM ∥i0∥+σmax n=1,...,Nc{meas(Ωn)}∥c∥L2(0,T )eT/CSM for every t∈[0, T]. Proof. Since F−1is globally Lipschitz continuous in RNc, from Theorem 2.15 in [2] we conclude that Problem 10 has a unique solution i=F−1(ℓ), with ℓ∈ C1([0, T ])Nc. Therefore, Problem 6 has a unique solution i∈ C0,1([0, T])Nc. Moreover, integrating the equation appearing in Problem 10 in (0, t), ℓ(t)−ℓ(0) = −Zt 0F−1(ℓ(s)) + (c1(s)σmeas(Ω1), . . . , cNc(s)σmeas(ΩNc))Tds for every t∈[0, T ]. Thus, ∥ℓ(t)−ℓ(0)∥ ≤ Zt 0F−1(ℓ(s))ds +Zt 0(c1(s)σmeas(Ω1), . . . , cNc(s)σmeas(ΩNc))Tds ≤Zt 0F−1(ℓ(s)) −F−1(ℓ(0))ds +Tσmax n=1,...,Nc{meas(Ωn)}∥c∥L2(0,T )+∥i0∥ ≤Zt 0 1 CSM ∥ℓ(s)−ℓ(0)∥ds +Tσmax n=1,...,Nc{meas(Ωn)}∥c∥L2(0,T )+∥i0∥. Then, taking Gronwall’s inequality into account (see, for instance, [13], Lemma 1.4.1), ∥ℓ(t)−ℓ(0)∥ ≤ Tσmax n=1,...,Nc{meas(Ωn)}∥c∥L2(0,T )+∥i0∥eT/CSM . Now, since i(t) = F−1(ℓ(t)), ∥i(t)∥−∥i0∥ ≤ ∥i(t)−i0∥=F−1(ℓ(t)) −F−1(ℓ(0))≤1 CSM ∥ℓ(t)−ℓ(0)∥ ≤T CSM σmax n=1,...,Nc{meas(Ωn)}∥c∥L2(0,T )+∥i0∥eT/CSM . 9 Similarly, for the second term on the right-hand side of eq. (4.3) we have ∆t ℓ X k=1 DF(i0)−Fh(i0),i(tk)−ik hE ≤ε2 2∥F(i0)−Fh(i0)∥2+1 2ε2∆t ℓ X k=1 i(tk)−ik h 2 ≤Ch2εε2 2∥i0∥2+||Br||2 L2(Ωpm)3+1 2ε2∆t ℓ X k=1 i(tk)−ik h 2 .(4.7) Moreover, regarding the third term on the right-hand side of eq. (4.3), writing again eq. (4.4) and using summation by parts, we obtain ∆t ℓ X k=1 Ztk 0 i∆t(t)−i(t)dt, i(tk)−ik h =*Ztℓ 0 i∆t(t)−i(t)dt, ∆t ℓ X k=1 i(tk)−ik h+− ℓ−1 X k=1 *Ztk+1 tk i∆t(t)−i(t)dt, ∆t k X m=1 (i(tm)−im h)+. (4.8) Using Young’s inequality in the first term of eq. (4.8), we deduce *Ztℓ 0 i∆t(t)−i(t)dt, ∆t ℓ X k=1 i(tk)−ik h+ ≤ε3 2Ztℓ 0 i∆t(t)−i(t)dt 2 +1 2ε3∆t ℓ X k=1 i(tk)−ik h 2 ≤ε3T 2∥i∆t−i∥2 L2(0,T )Nc+1 2ε3∆t ℓ X k=1 i(tk)−ik h 2 , for every ε3>0. Concerning the second term in eq. (4.8), we have ℓ−1 X k=1 *Ztk+1 tk i∆t(t)−i(t)dt, ∆t k X m=1 (i(tm)−im h)+ ≤ ℓ−1 X k=1 ∥i∆t−i∥L2(tk,tk+1)Nc√∆t∆t k X m=1 (i(tm)−im h) ≤ ℓ−1 X k=1 β1 2∥i∆t−i∥2 L2(tk,tk+1)Nc+ ∆t ℓ−1 X k=1 1 2β1∆t k X m=1 (i(tm)−im h) 2 ≤β1 2∥i∆t−i∥2 L2(0,T )Nc+ ∆t ℓ−1 X k=1 1 2β1∆t k X m=1 (i(tm)−im h) 2 ,(4.9) for every β1>0. Finally, we can bound analogously the fourth term on the right-hand side of eq. (4.3), 16 obtaining ∆t ℓ X k=1 Ztk 0 g(t)−g∆t(t)dt, i(tk)−ik h≤ε4T+β2 2∥g−g∆t∥2 L2(0,T )Nc +1 2ε4∆t ℓ X k=1 i(tk)−ik h 2 + ∆t ℓ−1 X k=1 1 2β2∆t k X m=1 (i(tm)−im h) 2 ,(4.10) for every ε4, β2>0. Using eqs. (4.5) to (4.10) in eq. (4.3) and rearranging the terms appropriately, we conclude CSM −1 2ε1ℓ X k=1 ∆ti(tk)−ik h 2+1 2−1 2ε2−1 2ε3−1 2ε4∆t ℓ X k=1 i(tk)−ik h 2 ≤ε3T+β1 2∥i∆t−i∥2 L2(0,T )Nc+ε4T+β2 2∥g−g∆t∥2 L2(0,T )Nc +h2εC1ε2+C2Tε1 2∥i0∥2+||Br||2 L2(Ωpm)3 +h2εC2ε1T 2∥g∥2 L2(0,T )Nc+ ∆t1 2β1 +1 2β2ℓ−1 X k=1 ∆t k X m=1 (i(tm)−im h) 2 , for every ε1, ε2, ε3, ε4, β1, β2>0and ℓ= 1, . . . , M. Thus, taking ε1>1 2CSM >0,ε2=ε3=ε4= 3ε1 (1−CSM )ε1+1 >0and β1=β2=2ε1 2CSM ε1−1>0, ℓ X k=1 ∆ti(tk)−ik h 2+∆t ℓ X k=1 i(tk)−ik h 2 ≤C∥i∆t−i∥2 L2(0,T )Nc+∥g−g∆t∥2 L2(0,T )Nc+Ch2ε∥i0∥2+||Br||2 L2(Ωpm)3+∥g∥2 L2(0,T )Nc + ∆t ℓ−1 X k=1 ∆t k X m=1 i(tk)−ik h 2 . Now, using the discrete Gronwall inequality (see Lemma 1.4.2 in [13]), we conclude M X k=1 ∆ti(tk)−ik h 2≤ M X k=1 ∆ti(tk)−ik h 2+∆t M X k=1 i(tk)−ik h 2 ≤C∥i∆t−i∥2 L2(0,T )Nc+∥g−g∆t∥2 L2(0,T )Nc+h2ε∥i0∥2+||Br||2 L2(Ωpm)3+∥g∥2 L2(0,T )Nc. 5. Numerical Results In this section we report some numerical results obtained from a Fortran code that solves a problem equivalent to Problem 22, allowing us to confirm the convergence result stated in Theorem 25. This equivalence is proved in appendix Appendix A. At each time step, the nonlinearity is solved by means of the fixed-point algorithm proposed in [5]. To this end, we have solved an academic problem built from the analytical test presented in [6] for a linear case. In our setting, we consider sources given in terms of time17 Figure 2: Sketch of the domain Ω(left). Coarsest mesh (right). dependent voltage drops per unit length and replace the linear material with a permanent magnet and a nonlinear core. In fig. 2-left we show the problem domain Ωthat includes the cross sections of two coaxial copper wires, Ω1and Ω2, separated by a permanent magnet, Ωpm, and a ferromagnetic core, Ωnl. We assume that the core and the magnet are non-conducting, and that the copper domains carry a uniformly distributed current density, i.e., they are stranded conductors. Let us consider a cylindrical coordinate system (ρ, θ, z), with eρ,eθand ezthe corresponding local orthonormal basis. We assume that the zaxis is orthogonal to the domain at point O. To apply the 2D transient magnetic model analysed in this paper, we suppose that the current density of the sources is uniform in each of them and orthogonal to the computational domain. More precisely, J=Jz(ρ, t)ez=             i(t) πR2 1 ezin (0, R1)×[0, T], 0in (R1, R3)×[0, T], −i(t) π(R2 4−R2 3)ezin (R3, R4)×[0, T]. In this case, if the magnetic constitutive law in the permanent magnet is of the form H=νpmB−νpmBr with a remanent flux Br=Breθ,Br∈R, then all fields are independent of the azimuthal variable and the solution to the magnetostatics problem eqs. (2.1) to (2.3) is H=Hθ(ρ, t)eθ=                  ρ i(t) 2πR2 1 eθin (0, R1)×[0, T], i(t) 2πρeθin (R1, R3)×[0, T], i(t)1 2πρ +(R2 3−ρ2) 2π(R2 4−R2 3)ρeθin (R3, R4)×[0, T]. 18 In the ferromagnetic core, we will consider the nonlinear constitutive magnetic law given by Bθ=µ0Hθ+2Js πatan π(µr−1)µ0Hθ 2Js.(5.1) We notice that, from Corollary 2.2 in [12], it can be seen that the corresponding nonlinear reluctivity function satisfies eqs. (2.4) to (2.6). Following the same arguments as in [5], and using the notation γ:= (µr−1)µ0 4Js, the expression of the solution to the magnetostatics problem (2.8)–(2.13) can be obtained by integrating Bθin space (see appendix Appendix B). In particular, it can be seen that the magnetic vector potential vanishes at R4for every t∈[0, T ]. This property allows us to have a conductor, Ω2, that touches the boundary of the whole domain. Indeed, if we had considered a domain Ω0representing the air surrounding the device, the solution Awould be identically zero there. Moreover, the expression of the potential drops per unit length in Ω1and Ω2can also be analytically computed using eq. (2.7), obtaining c1(t) = i′(t) 8πν0−i′(t) 2πν0ln R2 R1+ν0 νpm ln R3 R2+R2 4 R2 4−R2 3 ln R4 R3 −Jsγi′(t) πln γ2i(t)2+R2 2 γ2i(t)2+R2 1−i(t) σπR2 1 , c2(t) = i′(t) π(R2 4−R2 3)2ν0R2 3R2 4 2ln R4 R3−R4 4−R4 3 8+i(t) σπ(R2 4−R2 3). For the numerical computations, we have used the geometrical data R1= 0.5m, R2= 0.75 m, R3= 1 m and R4= 1.25 m. Moreover, the copper coils electrical conductivity σis equal to 5.7×107(Ohm m)−1and the magnetic reluctivity of the vacuum ν0=1 4π×107H−1m; the material of the permanent magnet is characterised by νpm = 0.95ν0and the remanent flux density Br=Breθby Br= 1.3T. Moreover, we have considered µ0=1 ν0,µr= 5000 and Js= 1.75 T in the nonlinear material law of the ferromagnetic core. The considered source in the coils is the potential drop per unit length obtained for a current i(t) = 3000 cos(2πft)A with a frequency f= 50 Hz. Finally, the initial currents in the conductors are i0,1= 3000 A and i0,2=−3000 A, respectively. These initial conditions allow us to obtain the steady state current i(t)from the beginning of the simulation. We solve the problem in a source cycle (that is, in the time interval [0, T ] = [0,0.02] seconds) with several successively refined meshes and time steps, starting from the mesh shown in fig. 2-right and a step size ∆t=T 40 . We have computed the errors by comparing the numerical solutions with the analytical one given by i(t) = i(t),−i(t)T. Specifically, we have computed the relative error for currents {im h}M m=1 in the L2(0, T)-norm, that is, E∆t h:= PM m=1 ∆t|i(tm)−im h|21/2 PM m=1 ∆t|i(tm)|21/2. Table 1 shows these relative errors at different levels of discretization. We notice that, when we take a time step small enough, an O(h2)error decay can be observed (see last row in table 1). On the other hand, considering a mesh size small enough allows us to show the expected convergence order in time O(∆t)(see single-framed values in the last column in table 1). We notice that the continuous solution is such that A(t)|Ωn∈H2(Ωn),n= 1,2, for every t∈[0, T]. Thus, the corresponding part of the convergence order proved in Theorem 25 is O(h), which is less than the one numerically obtained. The improvement of this 19 hh 2 h 4 h 8 h 16 ∆t0.1138 0.0806 0.0774 0.0772 0.0771 ∆t 20.0895 0.0438 0.0390 0.0388 0.0388 ∆t 40.0831 0.0279 0.0199 0.0195 0.0195 ∆t 80.0815 0.0222 0.0105 0.0098 0.0098 ∆t 16 0.0812 0.0206 0.0062 0.0050 0.0049 ∆t 32 0.0811 0.0201 0.0045 0.0027 0.0025 ∆t 64 0.0811 0.0200 0.0040 0.0016 0.0012 ∆t 128 0.0811 0.0200 0.0039 0.0012 0.0006 ∆t 256 0.0811 0.0200 0.0038 0.0011 0.0004 ∆t 512 0.0811 0.0200 0.0038 0.0011 0.0003 Table 1: Relative errors E∆t h. order is due to the fact that the norm used in eq. (3.6) is || · ||H1(∪Nc n=1Ωn), while the L2∪Nc n=1Ωn-norm could have been used. In the L2∪Nc n=1Ωn-norm, the magnetostatics problem converges with order O(h2) for this particular example. Once the convergence order is checked, we illustrate in one single figure the simultaneous dependence on hand ∆tof the error for current iin the L2(0, T )-norm by choosing initial coarse values for both discretization step-sizes and, for each successively refined mesh, we take the value of ∆tproportional to h2 (see the double-framed values in table 1). fig. 3 shows a log-log plot of the corresponding relative errors E∆t h versus the number of degrees of freedom (d.o.f.). The slope of the curve shows again the convergence order O(h2+ ∆t). Appendix A. An Equivalence Result between Two Fully Discrete Schemes In Section 4, we analysed the numerical convergence for Problem 22. However, following [5], we have implemented a different numerical scheme for the same problem. In this section, we will prove the equivalence between the two discretizations. Let us consider the following discrete problem, which is the one used for the implementation: Problem 26. Given c(t)∈ C([0, T ])Nc,i0∈RNcand Br∈L2(Ωpm)3, find Am h∈ L0 h(Ω),m= 1, . . . , M, 20 102103104105 10−4 10−3 10−2 10−1 100 Relative error Number of d.o.f. O(h2+∆t) convergence Relative error Figure 3: E∆t hversus d.o.f. (log-log scale); ∆t=Ch2. such that ZΩ ν(·;|grad Am h|)grad Am h·grad Wh+1 ∆t Nc X n=1 ZΩnZΩn σAm h1 meas(Ωn)Wh =1 ∆t Nc X n=1 ZΩnZΩn σAm−1 h1 meas(Ωn)Wh− Nc X n=1 ZΩn σcn(tm)Wh+ZΩpm νpm (Br)⊥·grad Wh, for every Wh∈ L0 h(Ω), with A0 h∈ L0 h(Ω) the solution to the weak formulation ZΩ ν(·;|grad A0 h|)grad A0 h·grad Wh= Nc X n=1 ZΩn i0,n meas(Ωn)Wh+ZΩpm νpm (Br)⊥·grad Wh, for every Wh∈ L0 h(Ω). Theorem 27. Let im h∈RNc,m= 0, . . . , M, be the solution to Problem 22 and Am h∈ L0 h(Ω),m= 1, . . . , M, defined as the solution to Problem 13. Then, Am h∈ L0 h(Ω),m= 1, . . . , M, are a solution to Problem 26. Proof. Since Problem 13 has a unique solution, we deduce that Fh,n (im h) = ZΩn σAm h, n = 1, . . . , Nc, m = 1, . . . , M. Furthermore, since im h,m= 0, . . . , M, is the solution to Problem 22, im h,n =−1 ∆tFh,n (im h) + 1 ∆tFh,n im−1 h−cn(tm)σmeas(Ωn) =−1 ∆tZΩn σAm h+1 ∆tZΩn σAm−1 h−cn(tm)σmeas(Ωn) for n= 1, . . . , Nc,m= 1, . . . , M. Replacing these expressions in Problem 13 we conclude that Am h∈ L0 h(Ω), m= 1, . . . , M, are a solution to Problem 26. Theorem 28. Let Am h∈ L0 h(Ω),m= 1, . . . , M, be a solution to Problem 26. Let us define i0 h:= i0and 21 im h∈RNc,m= 1, . . . , M, such that im h,n := −1 ∆tZΩn σAm h+1 ∆tZΩn σAm−1 h−cn(tm)σmeas(Ωn),(A.1) n= 1, . . . , Nc. Then, im h,m= 0, . . . , M, are the solution to Problem 22. Proof. Let e Am h∈ L0 h(Ω),m= 0, . . . , M, be the solution to Problem 13 with the currents defined in eq. (A.1). In particular, e A0 h=A0 h. Furthermore, taking the definitions of im h,m= 0, . . . , M, into account, fields e Am h∈ L0 h(Ω),m= 1, . . . , M, are also the solutions to the following problems: ZΩ ν(·;|grad e Am h|)grad e Am h·grad Wh+1 ∆t Nc X n=1 ZΩnZΩn σAm h1 meas(Ωn)Wh =1 ∆t Nc X n=1 ZΩnZΩn σAm−1 h1 meas(Ωn)Wh− Nc X n=1 ZΩn σcn(tm)Wh+ZΩpm νpm (Br)⊥·grad Wh, for every Wh∈ L0 h(Ω). By subtracting to the above equalities those in Problem 26, we deduce that ZΩν(·;|grad e Am h|)grad e Am h−ν(·;|grad Am h|)grad Am h·grad Wh= 0 for every Wh∈ L0 h(Ω),m= 1, . . . , M. Since L0 h(Ω) ⊂H1 0(Ω), we can rewrite the last equality in the following way: DBe Am h−B(Am h), WhE= 0 for every Wh∈ L0 h(Ω),m= 1, . . . , M. In particular, taking Wh=e Am h−Am h, and since operator Bis strongly monotone, we obtain 0 = DBe Am h−B(Am h),e Am h−Am hE≥Me Am h−Am h 2 H1(Ω) , and therefore e Am h=Am hfor m= 1, . . . , M. Finally, taking the definition of Fhinto account, we have, Fh,n (im h) = ZΩn σe Am h=ZΩn σAm h, n = 1, . . . , Nc, m = 1, . . . , M, and then Fh(im h)+∆tim h=Fhim−1 h−(c1(tm)σmeas(Ω1), . . . , cNc(tm)σmeas(ΩNc))T, m= 1, . . . , M. Hence im hare the solution to Problem 22. Remark 29. In Theorems 27 and 28, we have seen that given a solution to Problem 26 we can compute the corresponding solution to Problem 22, and vice versa. Moreover, from the proof of the last theorem it can be deduced that Problem 26 has a unique solution. Indeed, given two solutions Am,1 h, Am,2 h∈ L0 h(Ω), m= 1, . . . , M, to Problem 26, let im,1 h,im,2 h∈RNc,m= 0, . . . , M, be the corresponding solutions to Problem 22, built as indicated in Theorem 28. Moreover, let e Am,1 h,e Am,2 h∈ L0 h(Ω),m= 1, . . . , M, be the solutions to Problem 13 corresponding to these currents. In the above proof, we have seen that e Am,1 h=Am,1 h and e Am,2 h=Am,2 h. Therefore, since Problem 22 is well-posed, im,1 h=im,2 h,m= 1, . . . , M. Consequently, e Am,1 h=e Am,2 h, and thus Am,1 h=Am,2 h,m= 1, . . . , M. 22 Appendix B. Analytical Expression of the Solution to the Numerical Example In this appendix we state the analytical expression of the solution to the magnetostatics problem appearing in Section 5, obtained from the magnetic flux density given in eq. (5.1): In (0, R1)×[0, T], A(ρ, t) = I(t) 2πν0−ρ2 2R2 1 + ln R2 R1+ν0 νpm ln R3 R2+R2 4 R2 4−R2 3 ln R4 R3 +2Js πR2atan γI(t) R2−R1atan γI(t) R1+γI(t) 2ln γ2I(t)2+R2 2 γ2I(t)2+R2 1+Br(R3−R2). In (R1, R2)×[0, T], A(ρ, t) = I(t) 2πν0−1 2+ ln R2 ρ+ν0 νpm ln R3 R2+R2 4 R2 4−R2 3 ln R4 R3 +2Js πR2atan γI(t) R2−ρatan γI(t) ρ+γI(t) 2ln γ2I(t)2+R2 2 γ2I(t)2+ρ2+Br(R3−R2). In (R2, R3)×[0, T], A(ρ, t) = I(t) 2πν0−1 2+ν0 νpm ln R3 ρ+R2 4 R2 4−R2 3 ln R4 R3+Br(R3−ρ). In (R3, R4)×[0, T], A(ρ, t) = I(t) 2πν0ρ2−R2 4 2(R2 4−R2 3)+R2 4 R2 4−R2 3 ln R4 ρ. Acknowledgements Work partially supported by FEDER and Xunta de Galicia (Spain) under grant GRC2013–014, by Ministerio de Economía y Competitividad (Spain) under the research project ENE2013–47867–C2–1–R, and by Ministerio de Educación, Cultura y Deporte (Spain) under grant FPU13/03409. References References [1] Abdulle, A., Huber, M. E., 2016. Error estimates for finite element approximations of nonlinear monotone elliptic problems with application to numerical homogenization. Numer. Methods Partial Differential Eq. 32 (3), 955–969. [2] Bagagiolo, F., 2016. Ordinary differential equations. www.science.unitn.it/∼bagagiol/noteODE.pdf. [3] Baumanns, S., 2012. Coupled electromagnetic field/circuit simulation: Modeling and numerical analysis. Ph.D. thesis, Mathematisch-Naturwissenschaftlichen Fakultät. [4] Benderskaya, G., 2007. Numerical methods for transient field-circuit coupled simulations based on the finite integration technique and a mixed circuit formulation. Ph.D. thesis, Fachbereich Elektrotechnik und Informationstechnik. [5] Bermúdez, A., Óscar Domínguez, Gómez, D., Salgado, P., 2013. Finite element approximation of nonlinear transient magnetic problems involving periodic potential drop excitations. Comput. Math. Appl. 65, 1200–1219. [6] Bermúdez, A., Rodríguez, R., Salgado, P., 2008. A finite element method for the magnetostatic problem in terms of scalar potentials. SIAM J Numer Anal 46 (3), 1338–1363. [7] Bossavit, A., 1999. Eddy currents in dimension 2: voltage drops. In: Int. Symp. on Theoret. Electrical Engineering (Proc. ISTET’99, W. Mathis, T. Schindler, eds). University Otto-von-Guericke (Magdeburg, Germany), pp. 103–107. 23 [8] Buffa, A., Maday, Y., Rapetti, F., 2001. A sliding mesh-mortar method for a two dimensional eddy current model of electric engines. ESAIM: Math Model Num Anal 35 (2), 191–228. [9] Heise, B., 1994. Analysis of a fully discrete finite element method for a nonlinear magnetic field problem. SIAM J Numer Anal 31 (3), 745–759. [10] Matthes, M., 2012. Numerical analysis of nonlinear partial differential-algebraic equations: A coupled and an abstract systems approach. Ph.D. thesis, Mathematisch-Naturwissenschaftlichen Fakultät. [11] Nicaise, S., Tröltzsch, F., 2013. A coupled maxwell integrodifferential model for magnetization processes. Math Nachr 287 (4), 432–452. [12] Pechstein, C., 2004. Multigrid-newton-methods for nonlinear magnetostatic problems. Ph.D. thesis, Johannes Kepler University Linz, Austria. [13] Quarteroni, A., Valli, A., 1994. Numerical Approximation of Partial Differential Equations. Springer Verlag. [14] Roubícek, T., 2013. Nonlinear Partial Differential Equations with Applications. Springer Basel. [15] Slodivka, M., Vrábel’, V., 2017. Existence and uniqueness of a solution for a field/circuit coupled problem. ESAIM: Mathematical Modelling and Numerical Analysis 51 (3), 1045–1061. [16] Touzani, R., Rappaz, J., 2013. Mathematical Models for Eddy Currents and Magnetostatics with Selected Applications. Springer. [17] Ugalde, G., Almandoz, G., Poza, J., González, A., 2009. Computation of iron losses in permanent magnet machines by multi-domain simulations. In: 2009 13th European Conference on Power Electronics and Applications. IEEE, pp. 1–10. [18] Xu, J., 1996. Two-grid discretization techniques for linear and nonlinear pdes. SIAM J Numer Anal 33 (5), 1759–1777. [19] Zeidler, E., 1990. Nonlinear Functional Analysis and its Applications II/B. Nonlinear Monotone Operators. SpringerVerlag, New York. [20] Zenísek, A., 1990. The finite element method for nonlinear elliptic equations with discontinuous coefficients. Numerische Mathematik 58 (1), 51–77. [21] Zhu, Z., Ng, K., Schofield, N., Howe, D., 2004. Improved analytical modelling of rotor eddy current loss in brushless machines equipped with surface-mounted permanent magnets. IEE Proceedings - Electric Power Applications 151 (6), 641–650. 24