Numerical analysis of nonlinear degenerate parabolic problems with application to eddy current models

This paper deals with the numerical analysis for a family of nonlinear degenerate parabolic problems. The model is spatially discretized using a finite element method; an implicit Euler scheme is employed for time discretization. We deduce sufficient conditions to ensure that the fully discrete problem has a unique solution and to prove quasi-optimal error estimates for the approximation. Finally, we propose a nonlinear degenerate parabolic problem that arises from electromagnetic applications in conductive nonlinear magnetic media and deduce its solubility and convergence by using the developed abstract theory, including some numerical results to confirm the obtained theoretical results.


Introduction
This paper is devoted to the fully discrete finite element approximation analysis for a class of nonlinear degenerate parabolic problems.In this work, we consider that a degenerate parabolic equation (see, for instance [1,2]) is an abstract evolution equation of the form d dt (Ru(t)) + Au(t) = f (t), (1.1) where R is a linear, bounded and monotone operator, and A is an operator which is not necessarily linear.The study for these problems is a relevant topic for the numerical analysis of partial differential equations, in particular for computing the eddy currents in electromagnetic field theory (see, for instance, [3][4][5][6][7][8]).For the case, which R is timedependent, we can refer to reader to [9].
The main purpose of this paper is to extend the results obtained in [9] where the operator A is assumed to be linear.To this purpose, we assume some reasonable conditions inspired in physical properties that allow to deal with ferromagnetic conducting materials.In this kind of materials, if the hysteresis effects and anisotropies are neglected, we can suppose that the reluctivity is a scalar function that has a nonlinear dependence on the absolute value of the magnetic induction; see, for instance, [10][11][12].
It is important to mention that the number of works referring to eddy current problems involving ferromagnetic conducting materials is very low.Among the papers dedicated to the study of this problem, we can cite the papers [10] and [13], which were based on the so-called multiharmonic-approach.In both cases, the authors apply the truncated Fourier series expansions to approximate numerical solutions of the eddy current problem.On the other hand, in [14] was presented a T − ψ formulation for a nonlinear eddy current model, and more recently, in [15], a nonlinear 2D transient magnetic fields with drop excitations was proposed.The numerical analysis for these two works was studied by using a finite element approximation and an implicit time discretization scheme.Besides, in [12], a FEM/BEM coupling was analyzed for a 3D nonlinear eddy current formulation based on a time-primitive of the electric field.Other important studies about numerical nonlinear degenerate parabolic problems have been proposed in [16,17].
This paper allows us to obtain convergence of conforming approximations in space (typically, the family of finite dimensional subspaces is defined by finite elements) and a backward Euler method for the time-approximation, under weaker assumptions on the solution.Following the approach proposed in [1], we obtain error estimates for the fully discrete scheme by assuming only a time-regularity of a suitable projection of the solution, which is a more natural regularity assumption for the solution of degenerate parabolic problems.In fact, this work is intended as a first step toward the analysis for nonlinear 3D transient eddy current models that happen in the presence of ferromagnetic materials, by considering variational formulations in terms of a timeprimitive of the electric field, which makes that the obtained problems are degenerate (see [18][19][20]).
The outline of the paper is as follows: In Sect.2, we summarize some results from [2] concerning nonlinear degenerate parabolic equations.In Sect.3, we propose a fully discrete approximation of the last problem by using finite dimensional subspaces to the space approximation and a backward Euler scheme in time.The results ensuring the quasi-optimal convergence of the approximation for numerical solution are shown in Sect. 4. Furthermore, the application of the theory to an eddy current model is studied in Sect.5, where we deduce its solvability and theoretical convergence by using the developed abstract theory.Finally, we report some numerical results that confirm the expected convergence of the method according to the theory.

A nonlinear degenerate parabolic problem
Let X and Y be two real separable Hilbert spaces such that X ⊂ Y with a continuous and dense embedding.Let us denote by X and Y the topological dual spaces of X and Y respectively, where X is obtained by using Y as a pivot space 1 .Besides, we consider a linear and bounded operator 2 R : Y → Y and an operator A : X → X which is not necessarily a linear operator.Let T > 0.Then, the nonlinear degenerate parabolic problem consists in given f ∈ L 2 (0, T ;X ) and u 0 ∈ X : where we denote by •, • X and •, • Y the duality pairings between X and Y in theirs corresponding dual spaces.In similar way, we will denote (•, •) X and (•, •) Y the inner products on X and Y respectively and • X , • Y the corresponding norms.
The following result shows sufficient conditions to obtain the existence and uniqueness of the solution for Problem 1, but first, we need to recall some definitions.

Definition 2.1 Given Z be a Hilbert space, let us consider an operator G
1 Each element in Y must be identified as an element of X by the Riesz representation theorem in Y , i.e. y, x X := (y, x) Y for any y ∈ Y and x ∈ X .Here, •, • X is the duality pairing between X and X , and (•, •) Y is the inner product in Y . 2 In this context, an operator is understood as a function between normed vector spaces.
• G is hemicontinuous if for every u, v ∈ Z the function t → G(u + tv), v Z is continuous.
The following result shows sufficient conditions to obtain the existence and uniqueness of solution for Problem 1: Theorem 2.1 Assume that R is a monotone and symmetric operator.Futhermore, we suppose the operator A is monotone and hemicontinuous, and there is a constant Then, there exists a solution of Problem 1.In addition, if A is a strictly monotone operator, the solution of Problem 1 is unique.

Fully discrete approximation for the nonlinear degenerate parabolic problem
In this section, we present the fully discrete approximation for the nonlinear degenerate parabolic problem which was introduced in the previous section.To this end, we assume that R and A are operators that satisfy the sufficient conditions given in Theorem 2.1 to guarantee the existence and uniqueness of solution of Problem 1.
The fully discrete approximation will be obtained by using the finite-element method in space and a backward-Euler scheme in time.Let {X h } h>0 be a sequence of finite-dimensional subspaces of X and let t n := n t, n = 0, . . ., N , be a uniform partition of [0, T ] with a time-step t := T /N .
For any finite sequence {θ n : n = 0, . . ., N } we denote Let u 0,h ∈ X h a given approximation of u 0 .The fully discrete approximation of Problem 1 reads as follows.
We can note that in each step n = 1, . . ., N , u n h is computed as the solution of the following problem: find u n h ∈ X h such that where A and F n are defined by We will use the theorem of Browder and Minty ([21, Theorem 26.A]) to deduce the existence and uniqueness of solution of Problem 2 for each n = 1, . . ., N .Since F n is linear and bounded, we need to prove that A : X h → X h is monotone, hemicontinuous and coercive.We will only prove the coerciveness, because the other properties are fulfilled due to that X h is a subset of X .In fact, by recalling that R is monotone, for any v ∈ X h , we have Thus, the coerciveness of A follows from (2.1).Consequently, we have proved the following result.

Theorem 3.1 If R and A are operators that satisfy the sufficient conditions given in Theorem 2.1. Then, the fully discrete Problem 2 has a unique solution u n
h ∈ X h with n = 1, . . ., N .

Error estimates for the fully discrete approximation
In this section, we will deduce some error estimates for the fully discrete approximation.To this aim, we start by recalling that R : Y → Y is a linear and bounded operator that is monotone and symmetric.Hence, if we denote the Riesz isomorphism by Y : Y → Y and define R := −1 Y R, we deduce that R : Y → Y is a linear and bounded operator.Moreover, R is monotone (i.e., ( Rv, v) Y ≥ 0 for any v ∈ Y ) and symmetric (i.e., ( Rv, w) Y = (v, Rw) Y for any v, w ∈ Y ).Thus, in virtue of [1] we can define Y + as the orthogonal space of Y 0 := ker R and let us consider Y 1/2 + be the completion of Y + under the norm given by v + := R1/2 v Y .
Next, we consider the orthogonal projection operator P + : Y → Y + defined by Therefore, Rv = R P + v for any v ∈ Y and furthermore On the other hand, it is clear that Furthermore, recalling that R is monotone and symmetric, it satisfies the Cauchy-Schwarz type inequality given by Thus, we obtain The previous inequality allows us to extend R : + → Y by using standard arguments of continuity and density.Furthermore, we have Next, we proceed to introduce a convenient splitting of the approximation error to obtain the estimates.To do that, we assume the assumptions of Theorem 2.1.and suppose that the solution of Problem 1 satisfies where V is a subspace of X .Besides, we assume that there exists a linear and bounded operator h : V → X h .

Remark 1
The subspace V can be understood as a functional space with a major regularity where the solution belongs in.For example, if X = H 1 ( ) (i.e., u ∈ L 2 (0, T ; H 1 0 ( ))), we can take V := H 1 0 ( ) ∩ H 1+s ( ) with s > 0 (i.e., in this case, we would assume u ∈ H 1 (0, T ; H 1 0 ( ) ∩ H 1+s ( ))).On the other hand, there are exist different examples of operators h satisfying the previous conditions.
1.In general, we can take V := X and define two orthogonal projection operators It is easy to see that these operators satisfy the following estimates for any w ∈ X: 2. In several cases, the space V allows us to define an interpolation operator I h : V → X h such that there exists a function γ (h) tending to zero as h goes to zero, such that for every w ∈ V it holds that In particular, if X := H 1 0 ( ), V := H 1 0 ( ) ∩ H 1+s ( )) and X h is given by the usual family of Lagrange finite element subspaces, we can consider the Lagrange interpolant operator L h : H 1 0 ( ) ∩ H 1+s ( )) → X h , which satisfies (see, for instance, [22]) Now, we define the error and consider its splitting where Furthermore, assuming + ) we denote Moreover, we can easily check that .
Proof Let n ∈ {1, . . ., N }, k ∈ {1, . . ., n} and v ∈ X h .Then, from Problem 1 and Problem 2, we obtain the error equation Thus, by using (4.1) and Furthermore, due to that By testing this previous identity with v = σ k h ∈ X h , we obtain Now, since A is a Lipschitz and strongly monotone operator, we have and On the other hand, using the fact that R is symmetric and monotone, it follows that Hence, using (4.11)-(4.13) in (4.10) together with Young's inequality, we deduce From (4.3) and Young's inequality, we have On the other hand, using (4.2) and Young's inequality Therefore, by replacing the previous inequalities in (4.14) and using the fact that R is a bounded operator, we obtain Thus, summing over k, we obtain Then, we have Therefore, by using the discrete Gronwall's Lemma (see, for instance, [23, Lemma 1.4.2]),we obtain Finally, by using the last inequality to estimate the second term in the right-hand term of (4.15), we have established the result.
Below, we show the theorem that allows to obtain the error estimates for nonlinear degenerate parabolic problems.

Theorem 4.1 Assume the conditions of Lemma 4.1. Let u be the solution to Problem 1 and u n
h , n = 1, . . ., N , that to Problem 2. If u ∈ H 1 (0, T ;V ) with P + u ∈ H 2 (0, T ;Y + ) then, there exists a constant C > 0, independent of h and t, such that Proof We estimate each term on the right hand of the inequality given in Lemma 4.1.
To start, we write σ 0 h = e 0 h − ρ 0 h and using the fact that R is symmetric and monotone to obtain On the other side, a Taylor's expansion shows that (4.17) Now, we can see that On the other hand, we can easily check that Cauchy-Schwarz inequality implies Consequently, using Lemma 4.1 and (4.16)-(4.18), the result is established.
Remark 2 If we use the operator 2,h given by (4.4) to define ρ h , then by using (4.5), it follows that Hence, we can replace these two inequalities in the estimate given in Theorem 4.1 to obtain a Céa-like result for the approximation error u(t n ) − u n h .

Application to the eddy current problem
Let us recall that eddy currents are usually modeled by the low-frequency Maxwell's equations [11, chapter 8].The aim for the eddy current problem is to determine the eddy currents induced in a three-dimensional conducting domain ˆ c by a given time dependent compactly supported current density J.The eddy current problem consists in the following system of equations which gets the connection between the magnetic field H, the magnetic intensity field B, and the electric field E: where σ represents the electric conductivity in the conductor.Furthermore, we assume that the relationship between the magnetic induction and intensity of the magnetic field is nonlinear.More exactly, in similar way as in [10,12,14], we suppose that where ν : R + → R represents the magnetic reluctivity in ferromagnetic materials.Besides, we assume that the condutivity is a piecewise smooth real-valued function satisfying Zlamal [3] has proposed a solution of a particular case of the eddy current system (5.1) by solving the following two-dimensional nonlinear degenerate parabolic problem for a given data function where the physical parameters σ is independent of x 3 .
It is easy to deduce the relationship between the eddy current problem (5.1) and the nonlinear degenerate parabolic equation Problem 3. In fact, the following proposition allows to write a particular solution of the eddy current model in terms of the solution u of Problem 3.

Well-posedness for the nonlinear eddy current formulation
Let ˆ ⊂ R 3 be a simply connected and bounded set containing the closure of ˆ c and Supp J, with J as in Proposition 5.1.In order to obtain a weak formulation for Problem 3, we have to consider the projection of both sets ˆ and the conducting domain ˆ c onto the plane x 1 x 2 , that will be denoted and c , respectively.Then, given u c 0 ∈ L 2 ( c ) and J d ∈ L 2 (0, T ; L 2 ( )), by multiplying (5.4) with v ∈ H 1 0 ( ) and integrating by parts over , we obtain the following weak formulation for the Problem 3.
Remark 3 Since the electric permittivity is zero outside of the conductor, left hand term of (5.4) suggests that the companion initial condition u(0) should be only known in c .Consequently, a natural choice for the initial condition is u(0)| c = u c 0 , where u c 0 is a given data in L 2 ( c ).This fact is consistent with the structure of the abstract degenerate Problem 1 and the natural definition of operator R given by (5.7) below.
Our next goal is to prove the existence and uniqueness of solution of last problem.Thus, we proceed to fit Problem 4 in the abstract structure of Problem 1, so we have to define Y := L 2 ( ) and X := H 1 0 ( ) with their usual inner products.Then, we can easily deduce that these spaces satisfy the corresponding properties of Sect.Besides, we can notice that where χ c is the characteristic function of c .Furthermore, and In addition, we define the function f ∈ L 2 (0, T ;X ) given by Finally, we should notice that the initial condition to Problem 4 is equivalent to Ru(0) = Ru 0 in Y where u 0 = ũc 0 is an extension of u c 0 to the whole .Now, let us consider the following assumptions on the nonlinear reluctivity ν : ) (5.11) The last assumptions on the reluctivity can be derived from natural properties of the physical BH-curves (see [24,Section 2 ]).

Theorem 5.1 There exists a unique solution u of Problem 4.
Proof The conditions (5.9)-(5.11)imply that the nonlinear operator A defined in (5.6) is strongly monotone and Lipschitz continuous with constants α = C P F α ν and κ = 3M ν , where C P F is the Poincaré-Friedrichs constant, α ν and M ν are indicated in (5.10)- (5.11), respectively (see [24]).Besides, the operator R is monotone and symmetric.Finally, by using (5.3), noticing that and applying Theorem 2.1, we conclude the proof.

Remark 4 It is easy to see that
consequently u| c belongs to the space

Error estimates for the fully discrete nonlinear degenerate formulation
The fully discrete approximation for Problem 4 is obtained by using finite element subspaces to define the corresponding family of finite dimensional subspaces of X .To this aim, in what follows we assume that and c are polygonal domains.Let {T h } h be a regular family of triangles meshes of such that each element K ∈ T h is contained either in c or in d := \ c .As usual, h stands for the largest diameter of the triangles K in T h .We will use standard Lagrange finite element to define X h subspace of H 1 0 ( ), i.e., where C 0 ( ) is the space of scalar continuous functions defined on and P 1 is the set of polynomials of degree not greater than one.Then, the fully discrete approximation for the nonlinear degenerate parabolic formulation is given by Problem 2 by using the notation (5.6)-(5.8).More precisely, given u 0,h ∈ X h an approximation of u 0 , the fully discrete approximation of Problem 4 can be read as follows.
Thus, in similar way as in Theorem 5.1 by using the properties (5.9)-(5.11),we can guarantee the existence and uniqueness of solution u n h ∈ X h , n = 1, . . ., N , the fully dicrete solution of Problem 5. Furthermore, the following result is a direct consequence of Theorem 4.1 (see Remark 2).Theorem 5.2 If u 0 ∈ H 1 0 ( ) and u ∈ H 1 (0, T ; H 1 0 ( )) with u| c ∈ H 2 (0, T ; L 2 ( c )) then there exists a constant C > 0, independent of h and t, such that where Remark 5 order to satisfy the assumption of the previous theorem, it is necessary that u 0 ∈ H 1 0 ( ).If u c 0 ∈ H 1 0 ( c ), then u 0 can be defined as the extension by zero of u c 0 to the whole .If u c 0 ∈ H 1 ( c ) with u c 0 / ∈ H 1 0 ( c ), we can define Finally, to obtain the asymptotic error estimates, we need to consider the Sobolev space H 1+s (Q) for 0 < s ≤ 1, where Q is either c or d , and define the space Let L h denote the Lagrange interpolant.Then, if v ∈ X ∩ H 1 0 ( ) then L h v ∈ X h and (see, for instance, [22]) (5.12)where | • | H 1 ( ) denotes the usual semi-norm in H 1 ( ).Consequently, we have the following result.
Theorem 5.3 If u 0 ∈ H 1 0 ( ) and there exists 0 < s ≤ 1 such that u ∈ H 1 (0, T ; X ∩ H 1 0 ( )) with u| c ∈ H 2 (0, T ; L 2 ( c )).Then, there exists a constant C > 0 independent of h and t, such that Proof It is a direct consequence of Theorem 4.1 (see also Remark 1) and the interpolation error estimate (5.12).

Remark 6
The previous result shows that the solution of Problem 5 provides a suitable approximation for the physical variable magnetic field intensity B(t n ) in the threedimensional computational domain ˆ .More precisely, we can use the relationship given in (5.5), to define for any n = 1, . . ., N , and propose the following approximation Consequently, by using Corollary 5.3, we deduce the following quasi-optimal error estimates

Numerical results
In this subsection we present some numerical results obtained with a MATLAB code which implements the numerical method described in Problem 5, to illustrate the convergence with respect to the discretization parameters.To this aim, we write the results obtained for a test problem with a known analytical solution.

Implementation issues
To compute the solution of Problem 5 at each time step, it is necessary to solve the nonlinear system: find To Newton's method, we calculate the Gâteaux-Fréchet derivative of the operator A : X → X defined by Av, w X := ν(|∇v|)∇v • ∇w.Now, following the lines of [24, Section 2.4], the derivative of A with respect to v in a direction z is given by Now, even though the norm application is not differentiable at the null vector, the application w On the other hand, in the case ∇v = 0, we have Then, for a fixed time step t n , Newton's method generates a sequence defined by u k+1 = u k + δ k , where the Newton's direction δ k is computed by solving the linear system: until the solution satisfies a suitable stopping criterion.

A test with known analytical solution
We consider domains ˆ , ˆ c such that their respective projection and c onto the plane x 1 x 2 , are given by (see Fig. 1) and assume T := 1.The right hand side J d , is chosen so that u(x 1 , x 2 , t) = e −5π t sin(π x 1 ) sin(π x 2 ), is the solution to Problem 4. We have taken σ = σ 0 = 10 6 ( m) −1 in c , the electric conductivity.The numerical method has been applied with several successively refined meshes and time-steps.The computed approximate solution has been compared with the analytical one, by calculating the relative percentage error in time-discrete norm from Corollary 5.3 and Remark 6 given by 100 .
On the other side, we will consider the nonlinear constitutive magnetic law as in [25].This magnetization curve is called the Fröhlich-Kennelly model, that is given by  with parameters a, b > 0, μ 0 = 4π × 10 −7 Hm −1 the magnetic permeability, and where T c and T 0 denote Curie and the room temperature in Celsius.Now, from (5.2) and (5.

Percentage error
B H-curve is shown in Fig. 2 for the values of the parameters used in the numerical simulation.
The Table 1 shows the relative errors for B in the L 2 (0, ; L 2 ( ˆ ))-norm.We can notice that by taking a small enough time-step t, we can observe the behavior of the error with respect to the space discretization (see the row corresponding to t/64).On the other side, by considering a small enough mesh-size h, we can check the convergence order with respect t (see column corresponding to h/16).Hence, we conclude an order the convergence O(h + t) for B, which confirm the theoretical results proved in Corollary 5.3.Finally, Fig. 3 shows log-log plot of the error of B, versus number of degrees of freedom (d.o.f).To report this, we have used values of t proportional to h (see the values within boxes in Table 1).The slope of curve is clearly an order of convergence O(h + t).Unfortunately, our study does not allow to estimate the error for the variable ∂ t u = E.However, we compute the following percentage error: 100 by using the same data defined above.Table 2 shows that the time discretization error dominates the error for the space discretization, even for a small enough time-step t (see the row corresponding to t/64).In fact, an order O( t) can be observed for the finest mesh (see column corresponding to h/16).Finally, Fig. 4 shows log-log plot of the error of E, versus number of degrees of freedom (d.o.f).To report this, we have used values of t proportional to h (see the values within boxes in Table 2).The slope of curve confirms an order of convergence O( t) = O(h + t).

Fig. 3 Fig. 4
Fig. 3 Percentage discretization error curve for B versus number of d.o.f.(log-log scale)