- Email: [email protected]
An E-based splitting method with a non-matching grid for eddy current equations with discontinuous coefficients
Applied Mathematics and Computation 170 (2005) 462–484 www.elsevier.com/locate/amc
An E-based splitting method with a non-matching grid for eddy current equations with discontinuous coefficients Tong Kang
a,*
, Kwang Ik Kim
b
a
b
Department of Applied Mathematics, School of Sciences, Communication University of China, Beijing 100024, PR China Department of Mathematics, Pohang University of Science and Technology, Pohang, 790-784, Republic of Korea
Abstract A new decoupled finite element method based on solving a vector and a scalar from the splitting of the electric field is presented to approximate time-dependent eddy current equations with discontinuous coefficients in a three-dimensional polyhedral domain. By introducing a Lagrangian multiplier, this method realizes an optimal control for the interface problem. A non-matching finite element grid on the interface is considered and an optimal energy-norm error estimate in finite time is obtained. 2005 Elsevier Inc. All rights reserved. Keywords: Eddy current equations; Interface problem; E-based splitting method; Edge and nodal elements; Non-matching grid; Error estimate
*
Corresponding author. E-mail addresses: [email protected] (T. Kang), [email protected] (K.I. Kim).
0096-3003/$ - see front matter 2005 Elsevier Inc. All rights reserved. doi:10.1016/j.amc.2004.12.020
T. Kang, K.I. Kim / Appl. Math. Comput. 170 (2005) 462–484
463
1. Introduction The fundamental task in eddy current computation is the following: Given a time-dependent solenoidal excited current and a conductor of prescribed shape, determine the induced eddy currents in the conductor. The eddy current model emerges from Maxwells equations by formally dropping the displacement current, which amounts to neglecting capacitive effects (space charges) and provides a reasonable approximation in the low frequency range and in the presence of high conductivity (see [1,2,6–10]). In [4], Ciarlet and Zou first studied both the nodal finite element method and the edge finite element method for Maxwell equations. Later, this idea was applied to time-dependent eddy current equations in [9] and a full discrete finite element method based on solving a vector and a scalar from splitting of E was suggested. Not only appropriate boundary conditions and the divergencefree properties of both the current density and the magnetic flux are considered in the splitting formulation, but this scheme is also a decoupled one between vector and scalar unknowns. The error estimate shows that it has a good precision. In [3], Chen et al. discussed finite element methods for Maxwell equations with discontinuous coefficients by introducing Lagrangian multipliers. In this paper, we are going to extend this method to eddy current equations with discontinuous coefficients in the conductor and present an E-based splitting finite element method with a non-matching grid. This paper is organized as follows. In Section 2, the eddy current interface problem and its E-based splitting formulations are described. In Section 3, its variational form based on an optimal control formulation is given and some characteristics of its weak solution are studied. In Section 4, a fully discrete decoupled scheme with non-matching grids is proposed and finally in Section 5, its optimal energy-norm error estimate is derived.
2. E-based splitting formulations This paper is concerned with the following eddy current equations in the space occupied by the conductor: 8 curl H ¼ rE þ Js in X ð0; T Þ; > > < ð2:1Þ curl E ¼ oðlHÞ in X ð0; T Þ; ot > > : divðlHÞ ¼ 0 in X ð0; T Þ: Here X R3 is a simply-connected Lipschitz polyhedral domain. We suppose that the magnetic permeability l and the conductivity r are discontinuous across an interface C X, where C is the boundary of a simply-connected Lipschitz polyhedral domain X1 with X1 X and X2 ¼ X n X1 . Without loss of
464
T. Kang, K.I. Kim / Appl. Math. Comput. 170 (2005) 462–484
generality, we consider only the case with l and r being two piecewise constant functions in the domain X, i.e., l1 in X1 ; r1 in X1 ; l¼ r¼ l2 in X2 ; r2 in X2 ; where li, and ri(i = 1, 2) are positive constants (see Fig. 1). Excitation is provided by a divergence-free source current Js. Here we assume that suppðJs Þ X and there is no flux of Js through oX. Let n be the unit outward normal vector on oX or on the interface C to oX1. Then, it is known that magnetic field must satisfy the jump conditions across C given by ½lH n ¼ 0;
½H n ¼ 0:
ð2:2Þ
Throughout this paper, the jump of any function A across C is defined as ½A :¼ A2 jC A1 jC with Ai = AjXi,i = 1, 2. For the sake of simplicity, let the boundary conditions for (2.1) be supplemented with H n ¼ 0 on oX ð0; T Þ;
ð2:3Þ
and due to no flux of Js through oX, we have rE n ¼ 0
on oX ð0; T Þ:
Fig. 1. Simply-connected polyhedral domains X1 and X2.
ð2:4Þ
T. Kang, K.I. Kim / Appl. Math. Comput. 170 (2005) 462–484
465
We also assume that the initial condition for (2.1) is Hðx; 0Þ ¼ H0 ðxÞ in X;
ð2:5Þ
which satisfies div (lH0) = 0. Further, in the third equation of (2.1), we introduce a vector A defined by lH ¼ curl A;
ð2:6Þ
so that we have oA r/; ð2:7Þ E¼ ot where / is an arbitrary scalar function. In this paper, we prescribe /1 = /2 = 0 on C with /i ¼ /jXi , i = 1, 2. In finite element analysis, we usually append to [A] · n = 0 on the interface to guarantee [curl A] Æ n = [lH] Æ n = 0 (see [8]). Now, using the above notations, the interface problem (2.1)–(2.5) is taken as the following E-based splitting formulations: ( r oA þ rr/ þ curl l1 curl A ¼ Js in X1 [ X2 ð0; T Þ; ot ð2:8Þ div r oA þ rr/ ¼ 0 in X [ X ð0; T Þ 1 2 ot with the interface and the boundary condition, respectively, given by 1 curl A n ¼ 0; ½A n ¼ 0; /1 ¼ /2 ¼ 0 on C ð0; T Þ; l
1 oA curl A n ¼ 0; r þ rr/ n ¼ 0 on oX ð0; T Þ: l ot ð2:9Þ Further, we have the initial conditions for (2.8) Aðx; 0Þ ¼ A0 ðxÞ;
and
/ðx; 0Þ ¼ /0 ðxÞ ¼ 0
in X;
ð2:10Þ
where A0 is the solution of the problem defined by 8 < curl 1 curl A0 ¼ curl H0 in X; l : 1 curl A n ¼ 0; A n ¼ 0 on oX: l
0
0
3. Variational formulation For solving the system (2.8)–(2.10), we are going to consider non-matching meshes on the interface C that allow the triangulations in X1 and X2 to be generated independently. This brings some difficulty to the convergence analysis since the resulting finite element space will be nonconforming for A. Thus we deal with the constraint [A] · n = 0 on C by the Lagrangian multiplier approach (see [3]).
466
T. Kang, K.I. Kim / Appl. Math. Comput. 170 (2005) 462–484
Firstly, we introduce some notations that will be used throughout the paper. Let Lp(0, T;X) denote the set of all strongly measurable functions u(t,Æ) from [0, T] into the Hilbert space X such that (RT p kuðtÞkX dt < 1 for 1 6 p < 1; 0 sup06t6T kuðtÞkX < 1
for
p ¼ 1:
s
We say u 2 H (0, T;X) (s is a positive integer number) if and only if u, are in L2(0, T;X). Let C(0, T;X) denote the space of continuously differentiable functions from (0, T) into X. We also define the spaces n o H ðcurl; XÞ ¼ v 2 L2 ðXÞ3 ; curl v 2 L2 ðXÞ3 ; s ou ; . . . ; oosut ot
H 0 ðcurl; XÞ ¼ fv 2 H ðcurl; XÞ; v n ¼ 0 on oXg; n o 3 3 H a ðcurl; XÞ ¼ v 2 H a ðXÞ ; curl v 2 H a ðXÞ ; a > 0 with the norms 1=2 2 2 kvk0;curl;X ¼ kvk0;X þ kcurl vk0;X and 1=2 2 2 : kvka;curl;X ¼ kvka;X þ kcurl vka;X Hereafter kÆk0,X denotes the L2(X)3-norm (or the L2(X)-norm for scalar functions) and kÆks,X denotes the norm of the Sobolev space Hs(X)3 (or Hs(X) for scalar functions) for s > 0. Additionally, (Æ,Æ)X stands for either the inner product in L2(X)3 (respectively, L2(X)) or the duality pairing. The above definitions for X are similarly defined for X1 and X2. For the remainder of this paper, C and Ci (i = 0, 1, 2) will always represent generic constants independent of the time step and the mesh size. We now denote the spaces X1 = H(curl;X1), X2 = H0(curl;X2), Y 1 ¼ H 1C;0 ðX1 Þ and Y 2 ¼ H 10 ðX2 Þ, where H 0 ðcurl; X2 Þ ¼ fv 2 H ðcurl; X2 Þ; v n ¼ 0 H 1C;0 ðX1 Þ ¼ fw 2 H 1 ðX1 Þ; w ¼ 0 H 10 ðX2 Þ ¼ fw 2 H 1 ðX2 Þ; w ¼ 0
on oX2 n Cg;
on Cg; on oX2 g;
and set X = X1 · X2 and Y = Y1 · Y2. Secondly, to establish an appropriate variational form for the system (2.8)– (2.10), we need to use a few important mathematical analysis tools adopted from [3]. Let 3
T ðCÞ ¼ fs 2 H 1=2 ðCÞ ; 9v 2 H 0 ðcurl; XÞ such that v n ¼ s
on Cg:
T. Kang, K.I. Kim / Appl. Math. Comput. 170 (2005) 462–484
467
Then, it is not difficult to see that T(C) is a Banach space under the norm defined by kskT ðCÞ ¼ inffkvk0;curl;X ; v 2 H 0 ðcurl; XÞ and v n ¼ s
on Cg:
Further, for any s 2 T(C), we define Z Z s; w1;C ¼ v curl w dx curl v w dx 8w 2 X 1 ; X1
s; w2;C ¼
Z
ð3:1Þ
ð3:2Þ
X1
curl v w dx
Z
X2
v curl w dx 8w 2 X 2 ;
ð3:3Þ
X2
where v 2 H0(curl;X) so that v · n = s on C. Here, (3.2) and (3.3) are independent of the choice of v 2 H(curl;X) where v · n = s on C. We also provide Lemmas 3.1 and 3.2 which are given in [3]. Lemma 3.1. For any s 2 T(C), we have the equality s; w1;C s; w2;C kskT ðCÞ ¼ sup : kwkX 1 X 2 w2X 1 X 2 In practice, Lemma 3.1 is rather inconvenient since it uses information from both X1 and X2 to define the norm on T(C). But noting that s; w1;C s; w1;C ksk1;C ¼ sup ; ksk2;C ¼ sup ð3:4Þ kwkX 1 kwkX 2 w2X 1 w2X 2 are also norms of T(C), we have the following lemma. Lemma 3.2. The norms kÆk1,C and kÆk2,C are equivalent to kÆkT(C). 2 Y in the second Now, taking any A 2 X in the first equation of (2.8) and / equation of (2.8) as test functions, and then integrating by parts leads to the following weak formulations: Problem (I). Find (A, /, p) 2 H1(0, T;X) · L2(0, T;Y) · L2(0, T;T(C)) such that (
) 2 X oA 1 þ rr/; AÞXi þ curl A; curl A ðr ot l Xi i¼1 þ p; A2;C p; A1;C ¼
2 X
ðJs ; AÞXi
8A 2 X ;
ð3:5Þ
i¼1
2 X oA 2 Y; þ rr/; r/ r ¼ 0 8/ ot Xi i¼1
ð3:6Þ
468
T. Kang, K.I. Kim / Appl. Math. Comput. 170 (2005) 462–484
A; p2;C A; p1;C ¼ 0
8p 2 T ðCÞ:
ð3:7Þ
The system (3.5)–(3.7) is consistent in the finite element scheme with nonmatching grids on C and can be derived from the interface problem (2.8)– (2.10) on the basis of an optimal control formulation. In order to analyze the solution of the Problem (I), we first study the existence and uniqueness of the solution for the following problem: Problem (II). Find (U, p) 2 H1(0, T;X) · L2(0, T;T(C)) such that (
) 2 X oU 1 ; UÞXi þ curl U; curl U ðr þ p; U2;C p; U1;C ot l Xi i¼1 2 X ¼ ðJs ; UÞXi 8U 2 X ; ð3:8Þ i¼1
U; p2;C U; p1;C ¼ 0
8p 2 T ðCÞ:
ð3:9Þ
Furthermore, we only need to discuss the following stationary variational problem, whose result can be extended to Problem (II) by the standard analytic method: Problem (III). For a given f 2 L2(X)3, find (Q, p) 2 X · T(C) such that 2 n o X ðai curl Q; curl QÞXi þ ðbi Q; QÞXi þ p; Q2;C p; Q1;C i¼1
¼
2 X ðf; QÞXi
8Q 2 X ;
ð3:10Þ
i¼1
Q; p2;C Q; p1;C ¼ 0 8p 2 T ðCÞ;
ð3:11Þ
where ai and bi are piecewise positive constants in Xi for i = 1, 2. Theorem 3.3. There exists a unique solution (Q, p) 2 X · T(C) for Problem (III). Proof. First, define a bilinear form a:X · X ! R such that 2 n o X aðu; vÞ ¼ ðai curl u; curl vÞXi þ ðbi u; vÞXi for u; v 2 X : i¼1
Then it is obvious that aðu; uÞ P a0 kuk2X for some constant a0 > 0.
ð3:12Þ
T. Kang, K.I. Kim / Appl. Math. Comput. 170 (2005) 462–484
469
Next, we verify the inf-sup condition, i.e., there exists a constant C > 0 such that sup B2X
s; B2 2;C s; B1 1;C P CkskT ðCÞ kBkX
8s 2 T ðCÞ;
ð3:13Þ
where B = Bi in Xi for i = 1, 2. Let B 2 H(curl;X1) be a solution of the problem given by ðcurl B; curl BÞX1 þ ðB; BÞX1 ¼ s; B1;C
8B 2 H ðcurl; X1 Þ:
Further, we define ( B in X1 ; e B¼ 0 in X2 :
ð3:14Þ
ð3:15Þ
e 2 X and Then, clearly B e ¼ kBk k Bk X 0;curl;X1 :
ð3:16Þ
Thus, it follows from (3.14) that e 2 2;C s; B e 1 1;C ¼ s; B1;C ¼ kBk s; B 0;curl;X1 2
which in turn yields, together with Lemma 3.2, e 2 2;C s; B e 1 1;C s; B ¼ kBk0;curl;X1 ¼ ksk1;C P CkskT ðCÞ : e k Bk X
Therefore the conclusion of this theorem follow from (3.12) and (3.13).
h
From the result of Theorem 3.3, we conclude that the solution (U, p) of Problem (II) exists and is unique. Further, by using an appropriate application of the Greens formula, we find that the Lagrange multiplier p satisfies the relation given by 1 p ¼ curl U n l
in T ðCÞ ð0; T Þ:
Now, for a solution U of the Problem (II) and a given / 2 L2(0, T;Y), let w = /dt P2 and A = U $w. Then, applying divJs = 0 with suppðJs Þ X, we have i¼1 ðJs ; r/ÞXi ¼ 0. Hence, replacing U ¼ A and U ¼ r/ for any 2 X Y , respectively, in (3.8) and (3.9), and then noting the definition ðA; /Þ n ¼ 0 on C, we see that A,/ and p satisfy Problem of T(C) with r/ n ¼ r/ (I), which in turn leads to the following theorem:
470
T. Kang, K.I. Kim / Appl. Math. Comput. 170 (2005) 462–484
Theorem 3.4. The solution (A, /, p) of Problem (I) exists, but only the Lagrange multiplier p is unique. Further, 1 p ¼ curl A n l
in T ðCÞ ð0; T Þ:
ð3:17Þ
Remark 3.1. The boundary condition l1 curl A n ¼ 0 on oX · (0, T) in the problem (2.8)–(2.10) is consistent with the definition of the trace space T(C). 4. Fully discrete decoupled finite element scheme with non-matching grids We now propose a finite element method to solve Problem (I), which allows a non-matching finite element grids on the interface C. Let Th1 and Th2 be shape regular triangulations of X1 and X2, respectively. Then, they induce naturally two finite triangulations Ch1 and Ch2 on C. Let Ch0 be an another shape regular triangulation over C. Note here that Chi ,i = 0, 1, 2, are allowed to be different for each other. However, we make the following reasonable assumption: H 1. Each triangle in Ch1 and Ch2 must be contained in some triangles of Ch0 . We introduce the Ne´de´lec H(curl, Xi)-conforming edge element space (see [11]) defined by X hi ¼ fvh 2 X i ; vh ¼ aK þ bK x on K; 8K 2 Thi g for i ¼ 1; 2; where aK and bK are two constant vectors. Then, it is known that any function vh 2 X hi is uniquely determined by the degrees of freedom in the set ME(v) of the moments on each element K 2 Thi , which is defined by Z M E ðvÞ ¼ v jds; eis an edge of K ; e
where j is the unit vector along the edge. Then, for any v 2 Hs(Xi)3 with curl v 2 Lp(Xi)3 (i = 1, 2), we can define an interpolation ph v 2 X hi , where s > 1/2 and p > 2. Here, phv has the same degrees of freedom (defined by ME(v)) as v on each K in Thi . Now, let T h0 ðCÞ be the Ne´de´lec T(C)-conforming edge element space defined by T h0 ðCÞ ¼ fsh 2 T ðCÞ; sh ¼ ðas þ bs nÞ n on any s 2 Ch0 ; as ; bs 2 R3 g: We also define the finite element spaces given by Y hi ¼ fuh 2 Y i ; uh jK 2 P1 ; 8K 2 Thi g for i ¼ 1; 2; where P1 is the space of linear polynomials. Further, let Ph be the standard interpolating operator and set X h ¼ X h1 X h2 and Y h ¼ Y h1 Y h2 . Then, we will assume the inf-sup condition:
T. Kang, K.I. Kim / Appl. Math. Comput. 170 (2005) 462–484
471
H 2. There exists a constant C* > 0 independent of h0, h1 and h2 such that sh ; whi i;C sup P C ksh kT ðCÞ 8sh 2 T h0 ðCÞ for i ¼ 1; 2: ð4:1Þ whi 2X hi kwhi k0;curl;Xi The assumption (H2) indicates that the mesh Ch0 should be coarse enough compared with the mesh Th1 or Th2 in order to stabilize the effect of the introduced Lagrangian multiplier. In subsection 4.4 of [3], by using a general compactness argument, (4.1) is verified to be valid at least when the mesh h1 or h2 is suitably small compared with h0. Further, for u 2 Ha(curl;Xi) and w 2 H1+a(Xi) (i = 1, 2, a > 1/2), we have the following interpolation error estimates (see [5]) ku ph uk0;Xi þ kcurl ðu ph uÞk0;Xi 6 Chai kuka;curl;Xi ;
ð4:2Þ
kw Ph wk0;Xi þ hi kw Ph wk1;Xi 6 Chi1þa kwk1þa;Xi :
ð4:3Þ
Now, let us divide the time interval (0, T) into M equally-spaced subintervals by using nodal points 0 = t0 < t1 < < tM = T with tn = ns and s = T/M. Denote the n-th subinterval by In = (tn1, tn]. For u 2 C(0, T;L2(X)3) and w 2 C(0, T;L2(X)), we define un = u(x, tn) and wn = w(x, tn) for 0 6 n 6 M. For Hilbert space X, we also denote the space ‘p(X) = {u = (u1, u2, . . . , uM):uk 2 X, 1 6 k 6 M} equipped with the norm 8
M1 1=p P kþ1 p > > < kuk‘p ðX Þ ¼ s ku kX for 1 6 p < 1; k¼0 ð4:4Þ > > : kuk‘1 ðX Þ ¼ max kukþ1 kpX for p ¼ 1: 06k6M1
Further, let kos u þ vk‘2 ðX Þ ¼
2 !1=2 M 1 X ukþ1 uk kþ1 þv s s k¼0
ð4:5Þ
X
for u,v 2 ‘2(X)3 if u0 is prescribed. Now, we give the following estimates (see [5]) for X = Ha(curl;Xi) (i = 1, 2, a P 0) which will be used in Section 5: nþ1 Z tnþ1 2 ou u un 2 dt 8u 2 H 1 ð0; T ; X Þ; 61 ð4:6Þ ot s s tn X X nþ1 Z tnþ1 2 2 u 2un þ un1 2 o u 61 dt 2 s tn1 ot2 X s X
8u 2 H 2 ð0; T ; X Þ:
ð4:7Þ
Finally, we are in the position to introduce a discrete version of Problem (I).
472
T. Kang, K.I. Kim / Appl. Math. Comput. 170 (2005) 462–484
Problem (IV). For n = 0, 1 ,. . . , M 1, let A0h ¼ ph A0 and /0h ¼ Ph /0 , then we have the following two steps for a discrete version of Problem (I): Step 1 Find ðAhnþ1 ; pnþ1 h Þ 2 X h T h0 ðCÞ such that (
) 2 X Anþ1 Anh 1 nþ1 h curl Ah ; curl Ah ; Ah r þ l s Xi Xi i¼1 þ phnþ1 ; Ah 2;C pnþ1 h ; Ah 1;C 2 n o X n ¼ ðJnþ1 ; A Þ ðrr/ ; A Þ 8Ah 2 X h ; h h s h Xi Xi
ð4:8Þ
i¼1
ph 2;C Ahnþ1 ; ph 1;C ¼ 0 Ahnþ1 ;
8ph 2 T h0 ðCÞ;
ð4:9Þ
Step 2 Find /hnþ1 2 Y h such that 2 X i¼1
rr/hnþ1 ; r/ h Xi
2 X Anþ1 Anh h ; r/ h ¼ r s Xi i¼1
2 Y h: 8/ h
ð4:10Þ
Next, we give the theorem which provides the existence and the uniqueness for the Problem (IV). nþ1 Theorem 4.1. Under the assumptions (H1) and (H2), the solution ðAnþ1 h ; ph Þ of (4.8) and (4.9) exists and is unique.
Proof. We only need to verify the inf-sup condition, i.e., for any sh 2 T h0 ðCÞ, there exists a constant C > 0 such that sup Bh 2X h
sh ; Bh2 2;C sh ; Bh1 1;C P Cksh kT ðCÞ ; kBh kX
where Bh ¼ Bhi in Xi for i = 1,2. Without loss of generality, suppose that i = 1 in (4.1). Then there exists wh1 2 X h1 such that sh ; wh1 1;C P C ksh kT ðCÞ : kwh1 k0;curl;X1
ð4:11Þ
Further, let Bh1 2 X h1 be a solution of the problem given by ðcurl Bh1 ; curl Bh1 ÞX1 þ ðBh1 ; Bh1 ÞX1 ¼ sh ; Bh1 1;C
8Bh1 2 X h1 :
ð4:12Þ
Then, replacing wh1 and Bh1 by Bh1 , respectively, in (4.11) and (4.12) as test functions, and then using (4.11) and (4.12), we get
T. Kang, K.I. Kim / Appl. Math. Comput. 170 (2005) 462–484
473
C ksh kT ðCÞ 6 kBh1 k0;curl;X1 6 ksh k1;C 6 Cksh kT ðCÞ : Now define Bh1 e Bh ¼ 0
in X1 ; in X2 :
Then we have e h k ¼ kBh k kB X 1 0;curl;X1
ð4:13Þ
and it follows from (4.13) that e h 2;C sh ; B e h 1;C ¼ sh ; Bh 1;C ¼ kBh k2 sh ; B 2 1 1 1 0;curl;;X1 which in turn leads to e h 2;C sh ; B e h 1;C sh ; Bh 1;C sh ; B 2 1 1 ¼ ¼ kBh1 k0;curl;X1 e kBh1 k0;curl;X1 k B h kX P C ksh kT ðCÞ : Therefore, (4.8) and (4.9) have a unique solution (Ahnþ1 ; pnþ1 h ), which is the conclusion of this theorem. h
5. Finite element error estimate In this section, an energy-norm error estimate is derived for the fully discrete finite element scheme (4.8)–(4.10). For the simplicity of the notation, let b:X · T(C) ! R be the bilinear form defined by bðQ; pÞ ¼ Q; p2;C Q; p1;C
8ðQ; pÞ 2 X T ðCÞ:
Further, define the elliptic projection operator Ph:X · T(C) ! Xh · Th0(C) such that, for any (Q, p) 2 X · T(C), we have (
) 2 X 1 curlðP h Q QÞ; curl Qh ðrP h Q rQ; Qh ÞXi þ l Xi i¼1 þ bðQh ; P h p pÞ ¼ 0; ph Þ ¼ 0 bðP h Q Q; h Þ 2 X h T h0 ðCÞ. Let Qh ¼ rwh 2 X h for any wh 2 Yh. Noting for all ðQh ; p $wh · n = 0 on C, we see from the definition of the spaces T(C) and Th0(C) that 2 X i¼1
ðrP h Q rQ; rwh ÞXi ¼ 0:
ð5:1Þ
474
T. Kang, K.I. Kim / Appl. Math. Comput. 170 (2005) 462–484
For the error analysis, we need a solution function A to be defined also in the time interval [s, T]. This can be done by extending A with some regularity from the time interval [0, T] to the time interval [s, T]. Thus, we shall always implicitly assume that A is well defined in terms of time variable t on the interval [s, T]. Now, we give the theorem for error estimates in terms of the projection operator Ph. Theorem 5.1. Under the assumptions (H1) and (H2), let (An+1, /n+1, pn+1) be the solution of Problem (I) at time t = tn+1 for n = 0, 1, . . . , M 1. Suppose that A 2 H 3 ð0; T ; H a ðcurl; X1 Þ H a ðcurl; X2 ÞÞ; / 2 H 2 ð0; T ; Y \ H 1þa ðX1 Þ H 1þa ðX2 ÞÞ; for some a > 1/2. Then we have 2 X
kAnþ1 P h Anþ1 k20;curl;Xi 6
i¼1
2 X
2a C i h2a i þ C 0 h0 ;
ð5:2Þ
i¼1
2 2 X Anþ1 An Anþ1 An Ph s s
0;curl;Xi
i¼1
6
2 X
2a C i h2a i þ C 0 h0 ;
ð5:3Þ
i¼1
2 2 X ðAnþ1 An Þ ðAn An1 Þ ðAnþ1 An Þ ðAn An1 Þ Ph s s
0;curl;Xi
i¼1
6
2 X
2a C i sh2a i þ C 0 sh0 :
i¼1
ð5:4Þ Proof. First, using the fact that p ¼ l1 curl A n :¼ T n, we can show by the first equation of (2.8) that, ! 2 oA 2 2 2 kcurl Tka;X1 6 C þ k/k1þa;X1 þ kJs ka;X1 ; ot a;X1 2 curl oT ot
a;X1
2 2 curl o T ot2
a;X1
2 2 ! 2 2 o/ oJs o A 6C 2 þ þ ; ot a;X1 ot 1þa;X1 ot a;X1 3 2 2 2 2 2 ! o A o / o Js : 6C þ þ ot3 ot2 ot2 a;X1 1þa;X1 a;X1
T. Kang, K.I. Kim / Appl. Math. Comput. 170 (2005) 462–484
475
Then, for T 2 H2(0, T;Ha(curl, X1)), we have ! 2 oA 2 2 2 2 kTka;curl;X1 6 C kAka;curl;X1 þ k/k1þa;X1 þ þ kJs ka;X1 ; ot a;X1 2 oT 6C ot a;curl;X1 2 2 o T 6C ot2 a;curl;X1
2 2 2 ! 2 2 oA o/ oJs o A þ þ 2 þ ; ot ot 1þa;X1 ot a;X1 ot a;X1 a;curl;X1 2 2 2 2 3 2 2 2 ! o A o / o A o Js þ þ þ : ot2 ot2 ot3 ot2 a;curl;X1 1þa;X1 a;X1 a;X1
Next, we introduce a triangulation Th0 in X1 whose restriction on C coincides with Ch0 . Let X h0 be the Ne´de´lec H(curl, X1)-conforming edge element over the mesh Th0 . Then, from the definition of Th0(C), we know that T h0 ðCÞ ¼ fvh n; vh 2 X h0 g. Further, for the projection operator Ph, we have 2 X kAnþ1 P h Anþ1 k20;curl;Xi þ kpnþ1 P h pnþ1 k2T ðCÞ i¼1
6
2 X i¼1
2
C i inf kAnþ1 Qh k0;curl;Xi þ C 0 Qh 2X h
inf
qh 2T h0 ðCÞ
2
kpnþ1 qh kT ðCÞ ;
ð5:5Þ
2 nþ1 2 2 X Anþ1 An p pn Anþ1 An pnþ1 pn P P þ h h s s s s T ðCÞ 0;curl;Xi i¼1 2 X Anþ1 An 2 Qh 6 C i inf Qh 2X h s i¼1 nþ1 2 0;curl;Xi n p p qh þ C 0 inf ð5:6Þ ; qh 2T h0 ðCÞ s T ðCÞ 2 2 X ðAnþ1 An Þ ðAn An1 Þ ðAnþ1 An Þ ðAn An1 Þ P h s s 0;curl;Xi i¼1 nþ1 nþ1 n n n1 2 ðp pn Þ ðpn pn1 Þ ðp p Þ ðp p Þ Ph þ s s T ðCÞ 2 2 X ðAnþ1 An Þ ðAn An1 Þ Qh 6 C i inf Qh 2X h s i¼1 nþ1 2 0;curl;Xi n n n1 ðp p Þ ðp p Þ qh þ C 0 inf ð5:7Þ : qh 2T h0 ðCÞ s T ðCÞ Now, making use of (4.2) and then (4.6) and (4.7) to the first terms on the right hand side of the inequalities in 5.5, 5.6 and 5.7 yields
476
T. Kang, K.I. Kim / Appl. Math. Comput. 170 (2005) 462–484 2
2
nþ1 inf kAnþ1 Qh k0;curl;Xi 6 Ch2a ka;curl;Xi ; i kA
ð5:8Þ
Qh 2X hi
2 nþ1 n inf A sA Qh
Qh 2X hi
6
Ch2a i
2 max oA
06t6T
0;curl;Xi
Anþ1 An 2 6 Ch2a i s
a;curl;Xi
ð5:9Þ
; a;curl;Xi
ot
nþ1 2 ðA An Þ ðAn An1 Þ Qh inf Qh 2X h s nþ1 2 An Þ ðAn An1 Þ 2a ðA 6 Chi s
0;curl;Xi
a;curl;Xi
6 Csh2a i
Z T tnþ1 2 2 2a o A dt 6 Csh i ot2 tn1 s a;curl;Xi
Z
2 2 o A dt: ot2 a;curl;Xi
ð5:10Þ
Further, by Lemma 3.2, the second terms on the right hand side of the inequalities in 5.5, 5.6 and 5.7 become inf
qh 2T h0 ðCÞ
2
2
kpnþ1 qh kT ðCÞ 6 C inf kTnþ1 vh k0;curl;X1 vh 2X h0
2
nþ1 6 Ch2a ka;curl;X1 ; 0 kT
nþ1 2 p pn q inf h qh 2T h ðCÞ s
T ðCÞ
ð5:11Þ
nþ1 2 T Tn v 6 C inf h vh 2X h s 0
nþ1 2 Tn 2a T 6 Ch0 s a;curl;X1 2 oT 6 Ch2a ; 0 max 06t6T ot a;curl;X1
0;curl;X1
ð5:12Þ
nþ1 2 ðp pn Þ ðpn pn1 Þ qh inf qh 2T h ðCÞ s
T ðCÞ
0
nþ1 2 ðT Tn Þ ðTn Tn1 Þ vh 6 C inf vh 2X h s
0;curl;X1
0
nþ1 ðT Tn Þ ðTn Tn1 Þ2 6 Ch2a 0 s
a;curl;X1
6 Csh2a 0
Z
T s
2 2 o T dt: ot2 a;curl;X1 ð5:13Þ
Hence, using (5.5)–(5.13), the results of this theorem follow.
h
T. Kang, K.I. Kim / Appl. Math. Comput. 170 (2005) 462–484
477
Now, we state the theorem on the error estimate which is relevant to the solutions of Problem (I) and Problem (IV). Theorem 5.2. Under the assumptions (H1) and (H2), let (A, /, p) and (Ah, /h, ph) be the solutions of Problem (I) and Problem (IV), respectively. Suppose that for some a > 1/2, A 2 H 3 ð0; T ; X \ H a ðcurl; X1 Þ H a ðcurl; X2 ÞÞ; / 2 H 2 ð0; T ; Y \ H 1þa ðX1 Þ H 1þa ðX2 ÞÞ: Then we have
2 2 2 X X 2 ðos Ah þ r/h Þ oA þ r/ þ kcurl ðAh AÞk‘1 ðL2 ðXi ÞÞ 2 2 ot ðL ðX ÞÞ ‘ i¼1 i¼1 i 6 Cs þ
2 X
2a C i h2a i þ C 0 h0 :
i¼1
Proof. Using the backward difference scheme at t = tk+1 for 3.5 ,3.6,3.7, we obtain ( )
2 X Akþ1 Ak 1 kþ1 kþ1 curl A ; curlA ;A r þ þ ðrr/ ; AÞXi l s Xi Xi i¼1 þ bðA; pkþ1 Þ ¼
2 2 X X ðJkþ1 ; AÞ ðrRkþ1 ; AÞXi ; s Xi i¼1
ð5:14Þ
i¼1
( ) 2 2 X X Akþ1 Ak kþ1 ; ; r/ r þ ðrr/ ; r/Þ ðrRkþ1 ; r/Þ ¼ Xi Xi s Xi i¼1 i¼1 ð5:15Þ pÞ ¼ 0 bðAkþ1 ;
ð5:16Þ
for all ðA; /; pÞ 2 X Y T ðCÞ, where Rkþ1 ¼
oA Akþ1 Ak 1 ðtkþ1 Þ ¼ ot s s
Z
tkþ1
ðt tk Þ tk
o2 A dt ot2
ð5:17Þ
with kRkþ1 k20;Xi 6 Cs
Z
tkþ1 tk
2 2 o A ot2 dt: 0;Xi
ð5:18Þ
478
T. Kang, K.I. Kim / Appl. Math. Comput. 170 (2005) 462–484
¼/ 2 Y h and Let A ¼ Ah 2 X h , / p¼ ph 2 T h0 ðCÞ and then subtracting h (5.14)–(5.16), respectively, from (4.8)–(4.10) at t = tk+1, we obtain ( 2 X ðAkþ1 Akþ1 Þ ðAkh Ak Þ ; Ah r h þ ðrrð/kh /kþ1 Þ; Ah ÞXi s Xi i¼1 )
1 kþ1 þ curl ðAkþ1 A Þ; curlA pkþ1 Þ þ bðAh ; pkþ1 h h h l Xi ¼
2 X ðrRkþ1 ; Ah ÞXi ;
ð5:19Þ
i¼1
( ) 2 X ðAkþ1 Akþ1 Þ ðAkh Ak Þ kþ1 kþ1 h ;r/h r þ ðrrð/h / Þ;r/h ÞXi s Xi i¼1 ¼
2 X Þ ; ðrRkþ1 ;r/ h Xi i¼1
ð5:20Þ
Akþ1 ; ph ¼ 0: b Akþ1 h
ð5:21Þ
¼ Akþ1 P h Akþ1 , gkþ1 ¼ /kþ1 Ph /kþ1 and g~kh ¼ /kh Further, let Hkþ1 h h h h kþ1 kþ1 k ¼ s2 gkþ1 in Ph / . Then, if we replace Ah ¼ sðHh Hh Þ in (5.19) and / h h (5.20), it follows from the definition of the projection operator Ph and the equality (5.1) that ( 2 X pffiffiffi kþ1 r H Hk 2 þ s rr~ gk ; Hkþ1 Hk h
h
h
h
0;Xi
h Xi
i¼1
kþ1 1 þs curl Hkþ1 Hkh h ; curl Hh l ¼
2 X
b kþ1 ; Hkþ1 Hk sr R h h
i¼1
) Xi
Xi
þ s rA
kþ1
kþ1
rP h A
; Hkþ1 h
Hkh Xi
;
ð5:22Þ 2 n X pffiffiffi kþ1 2 o kþ1 s rrg r Hh Hkh ; srgkþ1 þ h h Xi 0;Xi i¼1
¼
2 X kþ1 kþ1 kþ1 b kþ1 ; srgkþ1 sr R þ s rA rP A ; srg ; h h h Xi i¼1
Xi
ð5:23Þ
T. Kang, K.I. Kim / Appl. Math. Comput. 170 (2005) 462–484
479
where kþ1
P h Akþ1 Þ ðAk P h Ak Þ þ rð/kþ1 Ph /kþ1 Þ þ Rkþ1 : s ð5:24Þ
b kþ1 ¼ ðA R
Now, adding (5.22) and (5.23) and then multiplying both hand sides by 2 for the use of the fact that a2 + 2ab + b2 = (a + b)2 which is true for any real a and b, we have ( 2 X pffiffiffi kþ1 pffiffiffi 2 r H Hk þ srgkþ1 2 þ r Hkþ1 Hk þ sr~ gk : h
h
h
h
0;Xi
h
h
0;Xi
i¼1
) pffiffiffi k 2 pffiffiffi kþ1 2 kþ1 1 kþ1 k þ s rrgh 0;Xi s rr~gh Þ 0;Xi þ 2s curl Hh ; curl Hh Hh l Xi 2 X b kþ1 ; Hkþ1 Hk þ srgkþ1 2 sr R ¼ h h h i¼1
kþ1
þ2s rA
rP h A
kþ1
Xi
;Hkþ1 Hkh þ srgkþ1 h h Xi
:
Then, making some rearrangements and using Cauchy-Schwarzs inequality, we obtain ( 2 X pffiffiffi kþ1 pffiffiffi r H Hk þ srgkþ1 2 þ r Hkþ1 Hk þ sr~gk 2 h
i¼1
h
h
pffiffiffi kþ1 2 b r sR sr~ gkh
0;Xi
h
0;Xi
h
h
0;Xi
pffiffiffi kþ1 2 b þ r sR srgkþ1 h 0;Xi )
1 kþ1 curl Hkþ1 Hkh Þ h ; curlðHh l Xi 2 X b kþ1 ; Hkþ1 Hk þ sr~ ¼ 2 sr R gkh h h þ2s
Xi
i¼1
þ2s rAkþ1 rP h Akþ1 ; Hkþ1 Hkh þ srgkþ1 h h Xi 2 X pffiffiffi b kþ1 2 6 rs R i¼1
0;Xi
pffiffiffi 2 þ rðHkþ1 Hkh þ sr~gkh Þ0;Xi h
pffiffiffi 2 þ4s2 r Akþ1 P h Akþ1
0;Xi
1 pffiffiffi k kþ1 2 þ r Hkþ1 H þ srg : h h h 0;Xi 2
Further, 2dividing both hand sides by s and using the inequality 2 a(a b) P a /2b /2 which holds for any real a and b, we have
480
T. Kang, K.I. Kim / Appl. Math. Comput. 170 (2005) 462–484
( 2 2 pffiffiffi kþ1 kþ1 k 2 X pffiffiffi s kþ1 r Hh Hh þ rgkþ1 þ s b r R rg h h 2 s 0;Xi 0;Xi i¼1 2 2 ) 1 p1ffiffiffi curl Hk þ curl Hkþ1 h h pffiffiffi l l 0;Xi 0;Xi 2 X pffiffiffi kþ1 pffiffiffi b kþ1 2 2 b 6 s r R r~ gkh þ s r R 0;Xi
i¼1
0;Xi
pffiffiffi 2 þ4s r Akþ1 P h Akþ1
0;Xi
ð5:25Þ
:
Next, we will estimate the right hand side of the inequality (5.25) by using 4.2) and (4.3) and the Theorem 5.1. First we consider that 2 X b kþ1 b k 2 R R 0;Xi
i¼1
k kþ1 kþ1 2 X A P h Ak A P hA 6 s i¼1 k 2 A P h Ak Ak1 P h Ak1 s
0;curl;Xi
2 X r /kþ1 Ph /kþ1 r /k Ph /k 2 þ 0;Xi i¼1
þ
2 2 X X Rkþ1 2 þ Rk 2 0;Xi 0;Xi i¼1
6
2 X
i¼1
2a C i sh2a i þ C 0 sh0 þ
i¼1
þ
2 X
kþ1 / /k 2 C i h2a i 1þa;Xi
i¼1
2 X
Z
tkþ1
Cis
i¼1
6
2 X
tk1
2 2 o A ot2 dt 0;Xi
2a C i sh2a i þ C 0 sh0 þ
i¼1
2 X
C i sh2a i
Z
T 0
i¼1
2 o/ dt ot 1þa;Xi
tkþ1 2 2 o A dt þCs ot2 tk1 0;X Z
6
2 X i¼1
C i sh2a i
þ
C 0 sh2a 0
þ Cs
Z
tkþ1
tk1
2 2 o A ot2 dt: 0;X
T. Kang, K.I. Kim / Appl. Math. Comput. 170 (2005) 462–484
481
Then 2 X pffiffiffi b kþ1 b k 2 r R R þ r Ph /kþ1 Ph /k
0;Xi
i¼1
2 X b kþ1 b k 2 6 Ci R R
0;Xi
i¼1
2 þ Ph /kþ1 /k 1;Xi
2 2 2 X o A 2 C i /kþ1 /k 1;Xi ot2 dt þ tk1 0;X i¼1 i¼1 ! Z Z tkþ1 2 2 2 2 t kþ1 X o/ o A 2a 2a 6 C i shi þ C 0 sh0 þ Cs ot dt : ot2 dt þ tk1 tk 1;X 0;X i¼1 6
2 X
2a C i sh2a i þ C 0 sh0 þ Cs
Z
tkþ1
Consequently, using the inequality (a + b)2 6 (1 + s)a2 + (1/s + 1)b2 which holds for any real a, b and s > 0 with s small enough, we obtain 2 2 X pffiffiffi b kþ1 gkh r R r~ 0;Xi
i¼1
2 pffiffiffi kþ1 X 2 pffiffiffi b k b k þ r Ph /kþ1 Ph /k b ¼ R r R rgkh þ r R
6
i¼1 2 X i¼1
0;Xi
pffiffiffi k 2 b rgk ð1 þ sÞ r R h
0;Xi
2 CX pffiffiffi b kþ1 b k 2 þ r R R þ r Ph /kþ1 Ph /k s i¼1 0;Xi 2 2 X pffiffiffi k X 2 2a b rgk 6 ð1 þ sÞ þ C i h2a r R i þ C 0 h0 h
þC
Z
0;Xi
i¼1 tkþ1 2
tk1
2 Z o A dt þ C ot2 0;Xi
tkþ1 tk
i¼1
2 o/ dt: ot
ð5:26Þ
1;Xi
In addition, we have Akþ1 P Akþ1 Ak P Ak 2 h h 6 Ci s 0;Xi i¼1 i¼1 0;curl;Xi Z 2 2 2 2 t kþ1 X X 2 o A dt þ C i rð/kþ1 Ph /kþ1 Þ0;Xi þ Cis ot2 t 0;Xi k i¼1 i¼1 Z tkþ1 2 2 2 2 X X o A 2a 2a 2a kþ1 2 C i hi þ C 0 h0 þ C i hi / þ Cs 6 ot2 dt 1þa;Xi tk 0;X i¼1 i¼1 Z tkþ1 2 2 2 X o A 2a C i h2a ð5:27Þ 6 i þ C 0 h0 þ Cs ot2 dt;
2 X pffiffiffi b kþ1 2 rR
i¼1
2 X
tk
0;X
482
T. Kang, K.I. Kim / Appl. Math. Comput. 170 (2005) 462–484
and 2 2 X X pffiffiffi 2a rðAkþ1 P h Akþ1 Þ2 6 C i h2a i þ C 0 h0 : 0;Xi i¼1
ð5:28Þ
i¼1
Thus, making use of (5.26)–(5.28) and 5.25 becomes ( 2 2 pffiffiffi kþ1 kþ1 k 2 X pffiffiffi s kþ1 r Hh Hh þ rgkþ1 þ s b r R rg h h 0;Xi 2 s 0;Xi i¼1 pffiffiffi k 2 b rgk s r R h
0;Xi
6
2 X i¼1
2 2 ) 1 1 kþ1 k þ pffiffiffi curl Hh pffiffiffi curl Hh l l 0;Xi 0;Xi
pffiffiffi k 2 b rgk C i s2 r R h
0;Xi
þ
2 X
2a C i sh2a i þ C 0 sh0
i¼1
! Z tkþ1 2 2 Z tkþ1 2 o/ o A þ Cs ot dt : ot2 dt þ tk1 tk 1;X 0;X
ð5:29Þ
For the simplicity of the notation, we define R0 ¼
oA0 A0 AðsÞ ; s ot
b 0 ¼ ðA0 ph A0 Þ ðAðsÞ ph AðsÞÞ þ rð/0 Ph /0 Þ þ R0 : R s Summing both hand sides of (5.29) over k = 0,1, . . . , n and then applying the discrete Gronwalls inequality, we obtain 2 X
(
i¼1
2 n X pffiffiffi Hkþ1 Hkh kþ1 h þ rgh s r s
0;Xi
k¼0
2 1 nþ1 þpffiffiffi curl Hh l 0;Xi 6
2 X i¼1
þ
i¼1
0;Xi
C i h2a i
þ
C 0 h2a 0
0;Xi
)
pffiffiffi 0 2 b rg0 C i s r R h
2 X
pffiffiffi nþ1 2 b þ s r R rghnþ1
þ Cs
Z
þ
2 X i¼1
T s
2 1 0 Ci curl H p ffiffiffi h l 0;Xi
2 2 Z T o A dt þ ot2 0 0;X
! 2 o/ dt : ð5:30Þ ot 1;X
T. Kang, K.I. Kim / Appl. Math. Comput. 170 (2005) 462–484
483
Noting g0h ¼ 0 and then using (4.2) and Theorem 5.1, we get the inequalities 2 2 X X pffiffiffi 0 pffiffiffi 0 2 b rg0 Þk2 6 b k sk rð R sk r R 0;Xi h 0;Xi i¼1
i¼1
6
2 X
C i sh2a i
þ
C 0 sh2a 0
þ Cs
2
i¼1
6
2 X
2 2 o A ot2 dt s 0;X
Z
0
2a 2 C i sh2a i þ C 0 sh0 þ Cs :
i¼1
2 2 X 1 pffiffiffi curl H0 h l
6
0;Xi
i¼1
2 X
C i ðkph A0 A0 k20;curl;Xi þ kA0 P h A0 k20;curl;Xi Þ
i¼1
6
2 X
2a C i h2a i þ C 0 h0 ;
i¼1
Hence, by the definitions (4.4) and (4.5), (5.30) becomes 2 X
2
kos Hh þ rgh k‘2 ðL2 ðXi ÞÞ þ
i¼1
2 X
2
kcurl Hh k‘1 ðL2 ðXi ÞÞ
i¼1
6 Cs þ
2 X
2a C i h2a i þ C 0 h0 :
i¼1
Finally, applying the triangle inequality, we have
2 2 2 X X ðos Ah þ r/h Þ oA þ r/ þ kcurlðAh AÞk2‘1 ðL2 ðXi ÞÞ 2 2 ot ‘ ðL ðXi ÞÞ i¼1 i¼1 6 Cs þ
2 X
2a C i h2a i þ C 0 h0 :
i¼1
This completes the proof of this theorem.
h
Remark 5.1. From (2.6) and (2.7), we see that l1 curl Anþ1 and h n Anþ1 nþ1 nþ1 h Ah s r/nþ1 may be replaced by H and E , respectively. Hence, h h h the conclusion of Theorem 5.2 can be rewritten as 2 X
2
kHh Hk‘1 ðL2 ðXi ÞÞ 6 Cs þ
i¼1 2 X i¼1
2 X
2a C i h2a i þ C 0 h0 ;
i¼1
2
kEh Ek‘2 ðL2 ðXi ÞÞ 6 Cs þ
2 X i¼1
2a C i h2a i þ C 0 h0 :
484
T. Kang, K.I. Kim / Appl. Math. Comput. 170 (2005) 462–484
Acknowledgment Supported by Com2MaC-KOSEF, POSTECH BSRI Research fund 2004 and the scientific research fund of the State Administration of Radio, Film & TV of China.
References [1] R. Albanese, G. Rubinacci, Formulation of the eddy-current problem, IEE Proceedings 137 (1990) 16–22. [2] H. Ammari, A. Buffa, J.-C. Ne´de´lec, A justifiction of eddy currents model for the Maxwell equations, SIAM J. Appl. Math. 60 (2000) 1805–1823. [3] Z. Chen, Q. Du, J. Zou, Finite element methods with matching and nonmatching meshes for Maxwell equations with discontinuous coefficients, SIAM J. Numer. Anal. 37 (2000) 1542– 1570. [4] P. Ciarlet, J. Zou, Finite element convergence for the Darwin model to Maxwells equations, RAIRO Math. Modelling Numer. Anal. 31 (1997) 213–249. [5] P. Ciarlet, J. Zou, Fully discrete finite element approaches for time-dependent Maxwells equations, Numerische Mathematik 82 (1999) 193–219. [6] R. Hiptmair, Symmetric coupling for eddy current problems, SIAM J. Numer. Anal. 40 (2002) 41–65. [7] Q. Hu, J. Zou, a non-overlapping domain decomposition method for Maxwells equations in three dimensions, SIAM J. Numer. Anal. 42 (2003) 1682–1708. [8] J. Jin, The Finite Element Method in Eletromagnectics, second ed., John Wiley & Sons, INC., North-Holland, 2002. [9] T. Kang, K.I. Kim, An E-based splitting finite element method for time-dependent eddy current equations, in press. [10] C. Ma, Studies on A/ methods for time-dependent electromagnetic fields (in Chinese), Ph.D. thesis, Chinese Academy of Sciences, 2003. [11] J.-C. Ne´de´lec, Mixed finite elements in R3 , Numer. Math. 35 (1980) 315–341.
An E-based splitting method with a non-matching grid for eddy current equations with discontinuous coefficients Tong Kang
a,*
, Kwang Ik Kim
b
a
b
Department of Applied Mathematics, School of Sciences, Communication University of China, Beijing 100024, PR China Department of Mathematics, Pohang University of Science and Technology, Pohang, 790-784, Republic of Korea
Abstract A new decoupled finite element method based on solving a vector and a scalar from the splitting of the electric field is presented to approximate time-dependent eddy current equations with discontinuous coefficients in a three-dimensional polyhedral domain. By introducing a Lagrangian multiplier, this method realizes an optimal control for the interface problem. A non-matching finite element grid on the interface is considered and an optimal energy-norm error estimate in finite time is obtained. 2005 Elsevier Inc. All rights reserved. Keywords: Eddy current equations; Interface problem; E-based splitting method; Edge and nodal elements; Non-matching grid; Error estimate
*
Corresponding author. E-mail addresses: [email protected] (T. Kang), [email protected] (K.I. Kim).
0096-3003/$ - see front matter 2005 Elsevier Inc. All rights reserved. doi:10.1016/j.amc.2004.12.020
T. Kang, K.I. Kim / Appl. Math. Comput. 170 (2005) 462–484
463
1. Introduction The fundamental task in eddy current computation is the following: Given a time-dependent solenoidal excited current and a conductor of prescribed shape, determine the induced eddy currents in the conductor. The eddy current model emerges from Maxwells equations by formally dropping the displacement current, which amounts to neglecting capacitive effects (space charges) and provides a reasonable approximation in the low frequency range and in the presence of high conductivity (see [1,2,6–10]). In [4], Ciarlet and Zou first studied both the nodal finite element method and the edge finite element method for Maxwell equations. Later, this idea was applied to time-dependent eddy current equations in [9] and a full discrete finite element method based on solving a vector and a scalar from splitting of E was suggested. Not only appropriate boundary conditions and the divergencefree properties of both the current density and the magnetic flux are considered in the splitting formulation, but this scheme is also a decoupled one between vector and scalar unknowns. The error estimate shows that it has a good precision. In [3], Chen et al. discussed finite element methods for Maxwell equations with discontinuous coefficients by introducing Lagrangian multipliers. In this paper, we are going to extend this method to eddy current equations with discontinuous coefficients in the conductor and present an E-based splitting finite element method with a non-matching grid. This paper is organized as follows. In Section 2, the eddy current interface problem and its E-based splitting formulations are described. In Section 3, its variational form based on an optimal control formulation is given and some characteristics of its weak solution are studied. In Section 4, a fully discrete decoupled scheme with non-matching grids is proposed and finally in Section 5, its optimal energy-norm error estimate is derived.
2. E-based splitting formulations This paper is concerned with the following eddy current equations in the space occupied by the conductor: 8 curl H ¼ rE þ Js in X ð0; T Þ; > > < ð2:1Þ curl E ¼ oðlHÞ in X ð0; T Þ; ot > > : divðlHÞ ¼ 0 in X ð0; T Þ: Here X R3 is a simply-connected Lipschitz polyhedral domain. We suppose that the magnetic permeability l and the conductivity r are discontinuous across an interface C X, where C is the boundary of a simply-connected Lipschitz polyhedral domain X1 with X1 X and X2 ¼ X n X1 . Without loss of
464
T. Kang, K.I. Kim / Appl. Math. Comput. 170 (2005) 462–484
generality, we consider only the case with l and r being two piecewise constant functions in the domain X, i.e., l1 in X1 ; r1 in X1 ; l¼ r¼ l2 in X2 ; r2 in X2 ; where li, and ri(i = 1, 2) are positive constants (see Fig. 1). Excitation is provided by a divergence-free source current Js. Here we assume that suppðJs Þ X and there is no flux of Js through oX. Let n be the unit outward normal vector on oX or on the interface C to oX1. Then, it is known that magnetic field must satisfy the jump conditions across C given by ½lH n ¼ 0;
½H n ¼ 0:
ð2:2Þ
Throughout this paper, the jump of any function A across C is defined as ½A :¼ A2 jC A1 jC with Ai = AjXi,i = 1, 2. For the sake of simplicity, let the boundary conditions for (2.1) be supplemented with H n ¼ 0 on oX ð0; T Þ;
ð2:3Þ
and due to no flux of Js through oX, we have rE n ¼ 0
on oX ð0; T Þ:
Fig. 1. Simply-connected polyhedral domains X1 and X2.
ð2:4Þ
T. Kang, K.I. Kim / Appl. Math. Comput. 170 (2005) 462–484
465
We also assume that the initial condition for (2.1) is Hðx; 0Þ ¼ H0 ðxÞ in X;
ð2:5Þ
which satisfies div (lH0) = 0. Further, in the third equation of (2.1), we introduce a vector A defined by lH ¼ curl A;
ð2:6Þ
so that we have oA r/; ð2:7Þ E¼ ot where / is an arbitrary scalar function. In this paper, we prescribe /1 = /2 = 0 on C with /i ¼ /jXi , i = 1, 2. In finite element analysis, we usually append to [A] · n = 0 on the interface to guarantee [curl A] Æ n = [lH] Æ n = 0 (see [8]). Now, using the above notations, the interface problem (2.1)–(2.5) is taken as the following E-based splitting formulations: ( r oA þ rr/ þ curl l1 curl A ¼ Js in X1 [ X2 ð0; T Þ; ot ð2:8Þ div r oA þ rr/ ¼ 0 in X [ X ð0; T Þ 1 2 ot with the interface and the boundary condition, respectively, given by 1 curl A n ¼ 0; ½A n ¼ 0; /1 ¼ /2 ¼ 0 on C ð0; T Þ; l
1 oA curl A n ¼ 0; r þ rr/ n ¼ 0 on oX ð0; T Þ: l ot ð2:9Þ Further, we have the initial conditions for (2.8) Aðx; 0Þ ¼ A0 ðxÞ;
and
/ðx; 0Þ ¼ /0 ðxÞ ¼ 0
in X;
ð2:10Þ
where A0 is the solution of the problem defined by 8 < curl 1 curl A0 ¼ curl H0 in X; l : 1 curl A n ¼ 0; A n ¼ 0 on oX: l
0
0
3. Variational formulation For solving the system (2.8)–(2.10), we are going to consider non-matching meshes on the interface C that allow the triangulations in X1 and X2 to be generated independently. This brings some difficulty to the convergence analysis since the resulting finite element space will be nonconforming for A. Thus we deal with the constraint [A] · n = 0 on C by the Lagrangian multiplier approach (see [3]).
466
T. Kang, K.I. Kim / Appl. Math. Comput. 170 (2005) 462–484
Firstly, we introduce some notations that will be used throughout the paper. Let Lp(0, T;X) denote the set of all strongly measurable functions u(t,Æ) from [0, T] into the Hilbert space X such that (RT p kuðtÞkX dt < 1 for 1 6 p < 1; 0 sup06t6T kuðtÞkX < 1
for
p ¼ 1:
s
We say u 2 H (0, T;X) (s is a positive integer number) if and only if u, are in L2(0, T;X). Let C(0, T;X) denote the space of continuously differentiable functions from (0, T) into X. We also define the spaces n o H ðcurl; XÞ ¼ v 2 L2 ðXÞ3 ; curl v 2 L2 ðXÞ3 ; s ou ; . . . ; oosut ot
H 0 ðcurl; XÞ ¼ fv 2 H ðcurl; XÞ; v n ¼ 0 on oXg; n o 3 3 H a ðcurl; XÞ ¼ v 2 H a ðXÞ ; curl v 2 H a ðXÞ ; a > 0 with the norms 1=2 2 2 kvk0;curl;X ¼ kvk0;X þ kcurl vk0;X and 1=2 2 2 : kvka;curl;X ¼ kvka;X þ kcurl vka;X Hereafter kÆk0,X denotes the L2(X)3-norm (or the L2(X)-norm for scalar functions) and kÆks,X denotes the norm of the Sobolev space Hs(X)3 (or Hs(X) for scalar functions) for s > 0. Additionally, (Æ,Æ)X stands for either the inner product in L2(X)3 (respectively, L2(X)) or the duality pairing. The above definitions for X are similarly defined for X1 and X2. For the remainder of this paper, C and Ci (i = 0, 1, 2) will always represent generic constants independent of the time step and the mesh size. We now denote the spaces X1 = H(curl;X1), X2 = H0(curl;X2), Y 1 ¼ H 1C;0 ðX1 Þ and Y 2 ¼ H 10 ðX2 Þ, where H 0 ðcurl; X2 Þ ¼ fv 2 H ðcurl; X2 Þ; v n ¼ 0 H 1C;0 ðX1 Þ ¼ fw 2 H 1 ðX1 Þ; w ¼ 0 H 10 ðX2 Þ ¼ fw 2 H 1 ðX2 Þ; w ¼ 0
on oX2 n Cg;
on Cg; on oX2 g;
and set X = X1 · X2 and Y = Y1 · Y2. Secondly, to establish an appropriate variational form for the system (2.8)– (2.10), we need to use a few important mathematical analysis tools adopted from [3]. Let 3
T ðCÞ ¼ fs 2 H 1=2 ðCÞ ; 9v 2 H 0 ðcurl; XÞ such that v n ¼ s
on Cg:
T. Kang, K.I. Kim / Appl. Math. Comput. 170 (2005) 462–484
467
Then, it is not difficult to see that T(C) is a Banach space under the norm defined by kskT ðCÞ ¼ inffkvk0;curl;X ; v 2 H 0 ðcurl; XÞ and v n ¼ s
on Cg:
Further, for any s 2 T(C), we define Z Z s; w1;C ¼ v curl w dx curl v w dx 8w 2 X 1 ; X1
s; w2;C ¼
Z
ð3:1Þ
ð3:2Þ
X1
curl v w dx
Z
X2
v curl w dx 8w 2 X 2 ;
ð3:3Þ
X2
where v 2 H0(curl;X) so that v · n = s on C. Here, (3.2) and (3.3) are independent of the choice of v 2 H(curl;X) where v · n = s on C. We also provide Lemmas 3.1 and 3.2 which are given in [3]. Lemma 3.1. For any s 2 T(C), we have the equality s; w1;C s; w2;C kskT ðCÞ ¼ sup : kwkX 1 X 2 w2X 1 X 2 In practice, Lemma 3.1 is rather inconvenient since it uses information from both X1 and X2 to define the norm on T(C). But noting that s; w1;C s; w1;C ksk1;C ¼ sup ; ksk2;C ¼ sup ð3:4Þ kwkX 1 kwkX 2 w2X 1 w2X 2 are also norms of T(C), we have the following lemma. Lemma 3.2. The norms kÆk1,C and kÆk2,C are equivalent to kÆkT(C). 2 Y in the second Now, taking any A 2 X in the first equation of (2.8) and / equation of (2.8) as test functions, and then integrating by parts leads to the following weak formulations: Problem (I). Find (A, /, p) 2 H1(0, T;X) · L2(0, T;Y) · L2(0, T;T(C)) such that (
) 2 X oA 1 þ rr/; AÞXi þ curl A; curl A ðr ot l Xi i¼1 þ p; A2;C p; A1;C ¼
2 X
ðJs ; AÞXi
8A 2 X ;
ð3:5Þ
i¼1
2 X oA 2 Y; þ rr/; r/ r ¼ 0 8/ ot Xi i¼1
ð3:6Þ
468
T. Kang, K.I. Kim / Appl. Math. Comput. 170 (2005) 462–484
A; p2;C A; p1;C ¼ 0
8p 2 T ðCÞ:
ð3:7Þ
The system (3.5)–(3.7) is consistent in the finite element scheme with nonmatching grids on C and can be derived from the interface problem (2.8)– (2.10) on the basis of an optimal control formulation. In order to analyze the solution of the Problem (I), we first study the existence and uniqueness of the solution for the following problem: Problem (II). Find (U, p) 2 H1(0, T;X) · L2(0, T;T(C)) such that (
) 2 X oU 1 ; UÞXi þ curl U; curl U ðr þ p; U2;C p; U1;C ot l Xi i¼1 2 X ¼ ðJs ; UÞXi 8U 2 X ; ð3:8Þ i¼1
U; p2;C U; p1;C ¼ 0
8p 2 T ðCÞ:
ð3:9Þ
Furthermore, we only need to discuss the following stationary variational problem, whose result can be extended to Problem (II) by the standard analytic method: Problem (III). For a given f 2 L2(X)3, find (Q, p) 2 X · T(C) such that 2 n o X ðai curl Q; curl QÞXi þ ðbi Q; QÞXi þ p; Q2;C p; Q1;C i¼1
¼
2 X ðf; QÞXi
8Q 2 X ;
ð3:10Þ
i¼1
Q; p2;C Q; p1;C ¼ 0 8p 2 T ðCÞ;
ð3:11Þ
where ai and bi are piecewise positive constants in Xi for i = 1, 2. Theorem 3.3. There exists a unique solution (Q, p) 2 X · T(C) for Problem (III). Proof. First, define a bilinear form a:X · X ! R such that 2 n o X aðu; vÞ ¼ ðai curl u; curl vÞXi þ ðbi u; vÞXi for u; v 2 X : i¼1
Then it is obvious that aðu; uÞ P a0 kuk2X for some constant a0 > 0.
ð3:12Þ
T. Kang, K.I. Kim / Appl. Math. Comput. 170 (2005) 462–484
469
Next, we verify the inf-sup condition, i.e., there exists a constant C > 0 such that sup B2X
s; B2 2;C s; B1 1;C P CkskT ðCÞ kBkX
8s 2 T ðCÞ;
ð3:13Þ
where B = Bi in Xi for i = 1, 2. Let B 2 H(curl;X1) be a solution of the problem given by ðcurl B; curl BÞX1 þ ðB; BÞX1 ¼ s; B1;C
8B 2 H ðcurl; X1 Þ:
Further, we define ( B in X1 ; e B¼ 0 in X2 :
ð3:14Þ
ð3:15Þ
e 2 X and Then, clearly B e ¼ kBk k Bk X 0;curl;X1 :
ð3:16Þ
Thus, it follows from (3.14) that e 2 2;C s; B e 1 1;C ¼ s; B1;C ¼ kBk s; B 0;curl;X1 2
which in turn yields, together with Lemma 3.2, e 2 2;C s; B e 1 1;C s; B ¼ kBk0;curl;X1 ¼ ksk1;C P CkskT ðCÞ : e k Bk X
Therefore the conclusion of this theorem follow from (3.12) and (3.13).
h
From the result of Theorem 3.3, we conclude that the solution (U, p) of Problem (II) exists and is unique. Further, by using an appropriate application of the Greens formula, we find that the Lagrange multiplier p satisfies the relation given by 1 p ¼ curl U n l
in T ðCÞ ð0; T Þ:
Now, for a solution U of the Problem (II) and a given / 2 L2(0, T;Y), let w = /dt P2 and A = U $w. Then, applying divJs = 0 with suppðJs Þ X, we have i¼1 ðJs ; r/ÞXi ¼ 0. Hence, replacing U ¼ A and U ¼ r/ for any 2 X Y , respectively, in (3.8) and (3.9), and then noting the definition ðA; /Þ n ¼ 0 on C, we see that A,/ and p satisfy Problem of T(C) with r/ n ¼ r/ (I), which in turn leads to the following theorem:
470
T. Kang, K.I. Kim / Appl. Math. Comput. 170 (2005) 462–484
Theorem 3.4. The solution (A, /, p) of Problem (I) exists, but only the Lagrange multiplier p is unique. Further, 1 p ¼ curl A n l
in T ðCÞ ð0; T Þ:
ð3:17Þ
Remark 3.1. The boundary condition l1 curl A n ¼ 0 on oX · (0, T) in the problem (2.8)–(2.10) is consistent with the definition of the trace space T(C). 4. Fully discrete decoupled finite element scheme with non-matching grids We now propose a finite element method to solve Problem (I), which allows a non-matching finite element grids on the interface C. Let Th1 and Th2 be shape regular triangulations of X1 and X2, respectively. Then, they induce naturally two finite triangulations Ch1 and Ch2 on C. Let Ch0 be an another shape regular triangulation over C. Note here that Chi ,i = 0, 1, 2, are allowed to be different for each other. However, we make the following reasonable assumption: H 1. Each triangle in Ch1 and Ch2 must be contained in some triangles of Ch0 . We introduce the Ne´de´lec H(curl, Xi)-conforming edge element space (see [11]) defined by X hi ¼ fvh 2 X i ; vh ¼ aK þ bK x on K; 8K 2 Thi g for i ¼ 1; 2; where aK and bK are two constant vectors. Then, it is known that any function vh 2 X hi is uniquely determined by the degrees of freedom in the set ME(v) of the moments on each element K 2 Thi , which is defined by Z M E ðvÞ ¼ v jds; eis an edge of K ; e
where j is the unit vector along the edge. Then, for any v 2 Hs(Xi)3 with curl v 2 Lp(Xi)3 (i = 1, 2), we can define an interpolation ph v 2 X hi , where s > 1/2 and p > 2. Here, phv has the same degrees of freedom (defined by ME(v)) as v on each K in Thi . Now, let T h0 ðCÞ be the Ne´de´lec T(C)-conforming edge element space defined by T h0 ðCÞ ¼ fsh 2 T ðCÞ; sh ¼ ðas þ bs nÞ n on any s 2 Ch0 ; as ; bs 2 R3 g: We also define the finite element spaces given by Y hi ¼ fuh 2 Y i ; uh jK 2 P1 ; 8K 2 Thi g for i ¼ 1; 2; where P1 is the space of linear polynomials. Further, let Ph be the standard interpolating operator and set X h ¼ X h1 X h2 and Y h ¼ Y h1 Y h2 . Then, we will assume the inf-sup condition:
T. Kang, K.I. Kim / Appl. Math. Comput. 170 (2005) 462–484
471
H 2. There exists a constant C* > 0 independent of h0, h1 and h2 such that sh ; whi i;C sup P C ksh kT ðCÞ 8sh 2 T h0 ðCÞ for i ¼ 1; 2: ð4:1Þ whi 2X hi kwhi k0;curl;Xi The assumption (H2) indicates that the mesh Ch0 should be coarse enough compared with the mesh Th1 or Th2 in order to stabilize the effect of the introduced Lagrangian multiplier. In subsection 4.4 of [3], by using a general compactness argument, (4.1) is verified to be valid at least when the mesh h1 or h2 is suitably small compared with h0. Further, for u 2 Ha(curl;Xi) and w 2 H1+a(Xi) (i = 1, 2, a > 1/2), we have the following interpolation error estimates (see [5]) ku ph uk0;Xi þ kcurl ðu ph uÞk0;Xi 6 Chai kuka;curl;Xi ;
ð4:2Þ
kw Ph wk0;Xi þ hi kw Ph wk1;Xi 6 Chi1þa kwk1þa;Xi :
ð4:3Þ
Now, let us divide the time interval (0, T) into M equally-spaced subintervals by using nodal points 0 = t0 < t1 < < tM = T with tn = ns and s = T/M. Denote the n-th subinterval by In = (tn1, tn]. For u 2 C(0, T;L2(X)3) and w 2 C(0, T;L2(X)), we define un = u(x, tn) and wn = w(x, tn) for 0 6 n 6 M. For Hilbert space X, we also denote the space ‘p(X) = {u = (u1, u2, . . . , uM):uk 2 X, 1 6 k 6 M} equipped with the norm 8
M1 1=p P kþ1 p > > < kuk‘p ðX Þ ¼ s ku kX for 1 6 p < 1; k¼0 ð4:4Þ > > : kuk‘1 ðX Þ ¼ max kukþ1 kpX for p ¼ 1: 06k6M1
Further, let kos u þ vk‘2 ðX Þ ¼
2 !1=2 M 1 X ukþ1 uk kþ1 þv s s k¼0
ð4:5Þ
X
for u,v 2 ‘2(X)3 if u0 is prescribed. Now, we give the following estimates (see [5]) for X = Ha(curl;Xi) (i = 1, 2, a P 0) which will be used in Section 5: nþ1 Z tnþ1 2 ou u un 2 dt 8u 2 H 1 ð0; T ; X Þ; 61 ð4:6Þ ot s s tn X X nþ1 Z tnþ1 2 2 u 2un þ un1 2 o u 61 dt 2 s tn1 ot2 X s X
8u 2 H 2 ð0; T ; X Þ:
ð4:7Þ
Finally, we are in the position to introduce a discrete version of Problem (I).
472
T. Kang, K.I. Kim / Appl. Math. Comput. 170 (2005) 462–484
Problem (IV). For n = 0, 1 ,. . . , M 1, let A0h ¼ ph A0 and /0h ¼ Ph /0 , then we have the following two steps for a discrete version of Problem (I): Step 1 Find ðAhnþ1 ; pnþ1 h Þ 2 X h T h0 ðCÞ such that (
) 2 X Anþ1 Anh 1 nþ1 h curl Ah ; curl Ah ; Ah r þ l s Xi Xi i¼1 þ phnþ1 ; Ah 2;C pnþ1 h ; Ah 1;C 2 n o X n ¼ ðJnþ1 ; A Þ ðrr/ ; A Þ 8Ah 2 X h ; h h s h Xi Xi
ð4:8Þ
i¼1
ph 2;C Ahnþ1 ; ph 1;C ¼ 0 Ahnþ1 ;
8ph 2 T h0 ðCÞ;
ð4:9Þ
Step 2 Find /hnþ1 2 Y h such that 2 X i¼1
rr/hnþ1 ; r/ h Xi
2 X Anþ1 Anh h ; r/ h ¼ r s Xi i¼1
2 Y h: 8/ h
ð4:10Þ
Next, we give the theorem which provides the existence and the uniqueness for the Problem (IV). nþ1 Theorem 4.1. Under the assumptions (H1) and (H2), the solution ðAnþ1 h ; ph Þ of (4.8) and (4.9) exists and is unique.
Proof. We only need to verify the inf-sup condition, i.e., for any sh 2 T h0 ðCÞ, there exists a constant C > 0 such that sup Bh 2X h
sh ; Bh2 2;C sh ; Bh1 1;C P Cksh kT ðCÞ ; kBh kX
where Bh ¼ Bhi in Xi for i = 1,2. Without loss of generality, suppose that i = 1 in (4.1). Then there exists wh1 2 X h1 such that sh ; wh1 1;C P C ksh kT ðCÞ : kwh1 k0;curl;X1
ð4:11Þ
Further, let Bh1 2 X h1 be a solution of the problem given by ðcurl Bh1 ; curl Bh1 ÞX1 þ ðBh1 ; Bh1 ÞX1 ¼ sh ; Bh1 1;C
8Bh1 2 X h1 :
ð4:12Þ
Then, replacing wh1 and Bh1 by Bh1 , respectively, in (4.11) and (4.12) as test functions, and then using (4.11) and (4.12), we get
T. Kang, K.I. Kim / Appl. Math. Comput. 170 (2005) 462–484
473
C ksh kT ðCÞ 6 kBh1 k0;curl;X1 6 ksh k1;C 6 Cksh kT ðCÞ : Now define Bh1 e Bh ¼ 0
in X1 ; in X2 :
Then we have e h k ¼ kBh k kB X 1 0;curl;X1
ð4:13Þ
and it follows from (4.13) that e h 2;C sh ; B e h 1;C ¼ sh ; Bh 1;C ¼ kBh k2 sh ; B 2 1 1 1 0;curl;;X1 which in turn leads to e h 2;C sh ; B e h 1;C sh ; Bh 1;C sh ; B 2 1 1 ¼ ¼ kBh1 k0;curl;X1 e kBh1 k0;curl;X1 k B h kX P C ksh kT ðCÞ : Therefore, (4.8) and (4.9) have a unique solution (Ahnþ1 ; pnþ1 h ), which is the conclusion of this theorem. h
5. Finite element error estimate In this section, an energy-norm error estimate is derived for the fully discrete finite element scheme (4.8)–(4.10). For the simplicity of the notation, let b:X · T(C) ! R be the bilinear form defined by bðQ; pÞ ¼ Q; p2;C Q; p1;C
8ðQ; pÞ 2 X T ðCÞ:
Further, define the elliptic projection operator Ph:X · T(C) ! Xh · Th0(C) such that, for any (Q, p) 2 X · T(C), we have (
) 2 X 1 curlðP h Q QÞ; curl Qh ðrP h Q rQ; Qh ÞXi þ l Xi i¼1 þ bðQh ; P h p pÞ ¼ 0; ph Þ ¼ 0 bðP h Q Q; h Þ 2 X h T h0 ðCÞ. Let Qh ¼ rwh 2 X h for any wh 2 Yh. Noting for all ðQh ; p $wh · n = 0 on C, we see from the definition of the spaces T(C) and Th0(C) that 2 X i¼1
ðrP h Q rQ; rwh ÞXi ¼ 0:
ð5:1Þ
474
T. Kang, K.I. Kim / Appl. Math. Comput. 170 (2005) 462–484
For the error analysis, we need a solution function A to be defined also in the time interval [s, T]. This can be done by extending A with some regularity from the time interval [0, T] to the time interval [s, T]. Thus, we shall always implicitly assume that A is well defined in terms of time variable t on the interval [s, T]. Now, we give the theorem for error estimates in terms of the projection operator Ph. Theorem 5.1. Under the assumptions (H1) and (H2), let (An+1, /n+1, pn+1) be the solution of Problem (I) at time t = tn+1 for n = 0, 1, . . . , M 1. Suppose that A 2 H 3 ð0; T ; H a ðcurl; X1 Þ H a ðcurl; X2 ÞÞ; / 2 H 2 ð0; T ; Y \ H 1þa ðX1 Þ H 1þa ðX2 ÞÞ; for some a > 1/2. Then we have 2 X
kAnþ1 P h Anþ1 k20;curl;Xi 6
i¼1
2 X
2a C i h2a i þ C 0 h0 ;
ð5:2Þ
i¼1
2 2 X Anþ1 An Anþ1 An Ph s s
0;curl;Xi
i¼1
6
2 X
2a C i h2a i þ C 0 h0 ;
ð5:3Þ
i¼1
2 2 X ðAnþ1 An Þ ðAn An1 Þ ðAnþ1 An Þ ðAn An1 Þ Ph s s
0;curl;Xi
i¼1
6
2 X
2a C i sh2a i þ C 0 sh0 :
i¼1
ð5:4Þ Proof. First, using the fact that p ¼ l1 curl A n :¼ T n, we can show by the first equation of (2.8) that, ! 2 oA 2 2 2 kcurl Tka;X1 6 C þ k/k1þa;X1 þ kJs ka;X1 ; ot a;X1 2 curl oT ot
a;X1
2 2 curl o T ot2
a;X1
2 2 ! 2 2 o/ oJs o A 6C 2 þ þ ; ot a;X1 ot 1þa;X1 ot a;X1 3 2 2 2 2 2 ! o A o / o Js : 6C þ þ ot3 ot2 ot2 a;X1 1þa;X1 a;X1
T. Kang, K.I. Kim / Appl. Math. Comput. 170 (2005) 462–484
475
Then, for T 2 H2(0, T;Ha(curl, X1)), we have ! 2 oA 2 2 2 2 kTka;curl;X1 6 C kAka;curl;X1 þ k/k1þa;X1 þ þ kJs ka;X1 ; ot a;X1 2 oT 6C ot a;curl;X1 2 2 o T 6C ot2 a;curl;X1
2 2 2 ! 2 2 oA o/ oJs o A þ þ 2 þ ; ot ot 1þa;X1 ot a;X1 ot a;X1 a;curl;X1 2 2 2 2 3 2 2 2 ! o A o / o A o Js þ þ þ : ot2 ot2 ot3 ot2 a;curl;X1 1þa;X1 a;X1 a;X1
Next, we introduce a triangulation Th0 in X1 whose restriction on C coincides with Ch0 . Let X h0 be the Ne´de´lec H(curl, X1)-conforming edge element over the mesh Th0 . Then, from the definition of Th0(C), we know that T h0 ðCÞ ¼ fvh n; vh 2 X h0 g. Further, for the projection operator Ph, we have 2 X kAnþ1 P h Anþ1 k20;curl;Xi þ kpnþ1 P h pnþ1 k2T ðCÞ i¼1
6
2 X i¼1
2
C i inf kAnþ1 Qh k0;curl;Xi þ C 0 Qh 2X h
inf
qh 2T h0 ðCÞ
2
kpnþ1 qh kT ðCÞ ;
ð5:5Þ
2 nþ1 2 2 X Anþ1 An p pn Anþ1 An pnþ1 pn P P þ h h s s s s T ðCÞ 0;curl;Xi i¼1 2 X Anþ1 An 2 Qh 6 C i inf Qh 2X h s i¼1 nþ1 2 0;curl;Xi n p p qh þ C 0 inf ð5:6Þ ; qh 2T h0 ðCÞ s T ðCÞ 2 2 X ðAnþ1 An Þ ðAn An1 Þ ðAnþ1 An Þ ðAn An1 Þ P h s s 0;curl;Xi i¼1 nþ1 nþ1 n n n1 2 ðp pn Þ ðpn pn1 Þ ðp p Þ ðp p Þ Ph þ s s T ðCÞ 2 2 X ðAnþ1 An Þ ðAn An1 Þ Qh 6 C i inf Qh 2X h s i¼1 nþ1 2 0;curl;Xi n n n1 ðp p Þ ðp p Þ qh þ C 0 inf ð5:7Þ : qh 2T h0 ðCÞ s T ðCÞ Now, making use of (4.2) and then (4.6) and (4.7) to the first terms on the right hand side of the inequalities in 5.5, 5.6 and 5.7 yields
476
T. Kang, K.I. Kim / Appl. Math. Comput. 170 (2005) 462–484 2
2
nþ1 inf kAnþ1 Qh k0;curl;Xi 6 Ch2a ka;curl;Xi ; i kA
ð5:8Þ
Qh 2X hi
2 nþ1 n inf A sA Qh
Qh 2X hi
6
Ch2a i
2 max oA
06t6T
0;curl;Xi
Anþ1 An 2 6 Ch2a i s
a;curl;Xi
ð5:9Þ
; a;curl;Xi
ot
nþ1 2 ðA An Þ ðAn An1 Þ Qh inf Qh 2X h s nþ1 2 An Þ ðAn An1 Þ 2a ðA 6 Chi s
0;curl;Xi
a;curl;Xi
6 Csh2a i
Z T tnþ1 2 2 2a o A dt 6 Csh i ot2 tn1 s a;curl;Xi
Z
2 2 o A dt: ot2 a;curl;Xi
ð5:10Þ
Further, by Lemma 3.2, the second terms on the right hand side of the inequalities in 5.5, 5.6 and 5.7 become inf
qh 2T h0 ðCÞ
2
2
kpnþ1 qh kT ðCÞ 6 C inf kTnþ1 vh k0;curl;X1 vh 2X h0
2
nþ1 6 Ch2a ka;curl;X1 ; 0 kT
nþ1 2 p pn q inf h qh 2T h ðCÞ s
T ðCÞ
ð5:11Þ
nþ1 2 T Tn v 6 C inf h vh 2X h s 0
nþ1 2 Tn 2a T 6 Ch0 s a;curl;X1 2 oT 6 Ch2a ; 0 max 06t6T ot a;curl;X1
0;curl;X1
ð5:12Þ
nþ1 2 ðp pn Þ ðpn pn1 Þ qh inf qh 2T h ðCÞ s
T ðCÞ
0
nþ1 2 ðT Tn Þ ðTn Tn1 Þ vh 6 C inf vh 2X h s
0;curl;X1
0
nþ1 ðT Tn Þ ðTn Tn1 Þ2 6 Ch2a 0 s
a;curl;X1
6 Csh2a 0
Z
T s
2 2 o T dt: ot2 a;curl;X1 ð5:13Þ
Hence, using (5.5)–(5.13), the results of this theorem follow.
h
T. Kang, K.I. Kim / Appl. Math. Comput. 170 (2005) 462–484
477
Now, we state the theorem on the error estimate which is relevant to the solutions of Problem (I) and Problem (IV). Theorem 5.2. Under the assumptions (H1) and (H2), let (A, /, p) and (Ah, /h, ph) be the solutions of Problem (I) and Problem (IV), respectively. Suppose that for some a > 1/2, A 2 H 3 ð0; T ; X \ H a ðcurl; X1 Þ H a ðcurl; X2 ÞÞ; / 2 H 2 ð0; T ; Y \ H 1þa ðX1 Þ H 1þa ðX2 ÞÞ: Then we have
2 2 2 X X 2 ðos Ah þ r/h Þ oA þ r/ þ kcurl ðAh AÞk‘1 ðL2 ðXi ÞÞ 2 2 ot ðL ðX ÞÞ ‘ i¼1 i¼1 i 6 Cs þ
2 X
2a C i h2a i þ C 0 h0 :
i¼1
Proof. Using the backward difference scheme at t = tk+1 for 3.5 ,3.6,3.7, we obtain ( )
2 X Akþ1 Ak 1 kþ1 kþ1 curl A ; curlA ;A r þ þ ðrr/ ; AÞXi l s Xi Xi i¼1 þ bðA; pkþ1 Þ ¼
2 2 X X ðJkþ1 ; AÞ ðrRkþ1 ; AÞXi ; s Xi i¼1
ð5:14Þ
i¼1
( ) 2 2 X X Akþ1 Ak kþ1 ; ; r/ r þ ðrr/ ; r/Þ ðrRkþ1 ; r/Þ ¼ Xi Xi s Xi i¼1 i¼1 ð5:15Þ pÞ ¼ 0 bðAkþ1 ;
ð5:16Þ
for all ðA; /; pÞ 2 X Y T ðCÞ, where Rkþ1 ¼
oA Akþ1 Ak 1 ðtkþ1 Þ ¼ ot s s
Z
tkþ1
ðt tk Þ tk
o2 A dt ot2
ð5:17Þ
with kRkþ1 k20;Xi 6 Cs
Z
tkþ1 tk
2 2 o A ot2 dt: 0;Xi
ð5:18Þ
478
T. Kang, K.I. Kim / Appl. Math. Comput. 170 (2005) 462–484
¼/ 2 Y h and Let A ¼ Ah 2 X h , / p¼ ph 2 T h0 ðCÞ and then subtracting h (5.14)–(5.16), respectively, from (4.8)–(4.10) at t = tk+1, we obtain ( 2 X ðAkþ1 Akþ1 Þ ðAkh Ak Þ ; Ah r h þ ðrrð/kh /kþ1 Þ; Ah ÞXi s Xi i¼1 )
1 kþ1 þ curl ðAkþ1 A Þ; curlA pkþ1 Þ þ bðAh ; pkþ1 h h h l Xi ¼
2 X ðrRkþ1 ; Ah ÞXi ;
ð5:19Þ
i¼1
( ) 2 X ðAkþ1 Akþ1 Þ ðAkh Ak Þ kþ1 kþ1 h ;r/h r þ ðrrð/h / Þ;r/h ÞXi s Xi i¼1 ¼
2 X Þ ; ðrRkþ1 ;r/ h Xi i¼1
ð5:20Þ
Akþ1 ; ph ¼ 0: b Akþ1 h
ð5:21Þ
¼ Akþ1 P h Akþ1 , gkþ1 ¼ /kþ1 Ph /kþ1 and g~kh ¼ /kh Further, let Hkþ1 h h h h kþ1 kþ1 k ¼ s2 gkþ1 in Ph / . Then, if we replace Ah ¼ sðHh Hh Þ in (5.19) and / h h (5.20), it follows from the definition of the projection operator Ph and the equality (5.1) that ( 2 X pffiffiffi kþ1 r H Hk 2 þ s rr~ gk ; Hkþ1 Hk h
h
h
h
0;Xi
h Xi
i¼1
kþ1 1 þs curl Hkþ1 Hkh h ; curl Hh l ¼
2 X
b kþ1 ; Hkþ1 Hk sr R h h
i¼1
) Xi
Xi
þ s rA
kþ1
kþ1
rP h A
; Hkþ1 h
Hkh Xi
;
ð5:22Þ 2 n X pffiffiffi kþ1 2 o kþ1 s rrg r Hh Hkh ; srgkþ1 þ h h Xi 0;Xi i¼1
¼
2 X kþ1 kþ1 kþ1 b kþ1 ; srgkþ1 sr R þ s rA rP A ; srg ; h h h Xi i¼1
Xi
ð5:23Þ
T. Kang, K.I. Kim / Appl. Math. Comput. 170 (2005) 462–484
479
where kþ1
P h Akþ1 Þ ðAk P h Ak Þ þ rð/kþ1 Ph /kþ1 Þ þ Rkþ1 : s ð5:24Þ
b kþ1 ¼ ðA R
Now, adding (5.22) and (5.23) and then multiplying both hand sides by 2 for the use of the fact that a2 + 2ab + b2 = (a + b)2 which is true for any real a and b, we have ( 2 X pffiffiffi kþ1 pffiffiffi 2 r H Hk þ srgkþ1 2 þ r Hkþ1 Hk þ sr~ gk : h
h
h
h
0;Xi
h
h
0;Xi
i¼1
) pffiffiffi k 2 pffiffiffi kþ1 2 kþ1 1 kþ1 k þ s rrgh 0;Xi s rr~gh Þ 0;Xi þ 2s curl Hh ; curl Hh Hh l Xi 2 X b kþ1 ; Hkþ1 Hk þ srgkþ1 2 sr R ¼ h h h i¼1
kþ1
þ2s rA
rP h A
kþ1
Xi
;Hkþ1 Hkh þ srgkþ1 h h Xi
:
Then, making some rearrangements and using Cauchy-Schwarzs inequality, we obtain ( 2 X pffiffiffi kþ1 pffiffiffi r H Hk þ srgkþ1 2 þ r Hkþ1 Hk þ sr~gk 2 h
i¼1
h
h
pffiffiffi kþ1 2 b r sR sr~ gkh
0;Xi
h
0;Xi
h
h
0;Xi
pffiffiffi kþ1 2 b þ r sR srgkþ1 h 0;Xi )
1 kþ1 curl Hkþ1 Hkh Þ h ; curlðHh l Xi 2 X b kþ1 ; Hkþ1 Hk þ sr~ ¼ 2 sr R gkh h h þ2s
Xi
i¼1
þ2s rAkþ1 rP h Akþ1 ; Hkþ1 Hkh þ srgkþ1 h h Xi 2 X pffiffiffi b kþ1 2 6 rs R i¼1
0;Xi
pffiffiffi 2 þ rðHkþ1 Hkh þ sr~gkh Þ0;Xi h
pffiffiffi 2 þ4s2 r Akþ1 P h Akþ1
0;Xi
1 pffiffiffi k kþ1 2 þ r Hkþ1 H þ srg : h h h 0;Xi 2
Further, 2dividing both hand sides by s and using the inequality 2 a(a b) P a /2b /2 which holds for any real a and b, we have
480
T. Kang, K.I. Kim / Appl. Math. Comput. 170 (2005) 462–484
( 2 2 pffiffiffi kþ1 kþ1 k 2 X pffiffiffi s kþ1 r Hh Hh þ rgkþ1 þ s b r R rg h h 2 s 0;Xi 0;Xi i¼1 2 2 ) 1 p1ffiffiffi curl Hk þ curl Hkþ1 h h pffiffiffi l l 0;Xi 0;Xi 2 X pffiffiffi kþ1 pffiffiffi b kþ1 2 2 b 6 s r R r~ gkh þ s r R 0;Xi
i¼1
0;Xi
pffiffiffi 2 þ4s r Akþ1 P h Akþ1
0;Xi
ð5:25Þ
:
Next, we will estimate the right hand side of the inequality (5.25) by using 4.2) and (4.3) and the Theorem 5.1. First we consider that 2 X b kþ1 b k 2 R R 0;Xi
i¼1
k kþ1 kþ1 2 X A P h Ak A P hA 6 s i¼1 k 2 A P h Ak Ak1 P h Ak1 s
0;curl;Xi
2 X r /kþ1 Ph /kþ1 r /k Ph /k 2 þ 0;Xi i¼1
þ
2 2 X X Rkþ1 2 þ Rk 2 0;Xi 0;Xi i¼1
6
2 X
i¼1
2a C i sh2a i þ C 0 sh0 þ
i¼1
þ
2 X
kþ1 / /k 2 C i h2a i 1þa;Xi
i¼1
2 X
Z
tkþ1
Cis
i¼1
6
2 X
tk1
2 2 o A ot2 dt 0;Xi
2a C i sh2a i þ C 0 sh0 þ
i¼1
2 X
C i sh2a i
Z
T 0
i¼1
2 o/ dt ot 1þa;Xi
tkþ1 2 2 o A dt þCs ot2 tk1 0;X Z
6
2 X i¼1
C i sh2a i
þ
C 0 sh2a 0
þ Cs
Z
tkþ1
tk1
2 2 o A ot2 dt: 0;X
T. Kang, K.I. Kim / Appl. Math. Comput. 170 (2005) 462–484
481
Then 2 X pffiffiffi b kþ1 b k 2 r R R þ r Ph /kþ1 Ph /k
0;Xi
i¼1
2 X b kþ1 b k 2 6 Ci R R
0;Xi
i¼1
2 þ Ph /kþ1 /k 1;Xi
2 2 2 X o A 2 C i /kþ1 /k 1;Xi ot2 dt þ tk1 0;X i¼1 i¼1 ! Z Z tkþ1 2 2 2 2 t kþ1 X o/ o A 2a 2a 6 C i shi þ C 0 sh0 þ Cs ot dt : ot2 dt þ tk1 tk 1;X 0;X i¼1 6
2 X
2a C i sh2a i þ C 0 sh0 þ Cs
Z
tkþ1
Consequently, using the inequality (a + b)2 6 (1 + s)a2 + (1/s + 1)b2 which holds for any real a, b and s > 0 with s small enough, we obtain 2 2 X pffiffiffi b kþ1 gkh r R r~ 0;Xi
i¼1
2 pffiffiffi kþ1 X 2 pffiffiffi b k b k þ r Ph /kþ1 Ph /k b ¼ R r R rgkh þ r R
6
i¼1 2 X i¼1
0;Xi
pffiffiffi k 2 b rgk ð1 þ sÞ r R h
0;Xi
2 CX pffiffiffi b kþ1 b k 2 þ r R R þ r Ph /kþ1 Ph /k s i¼1 0;Xi 2 2 X pffiffiffi k X 2 2a b rgk 6 ð1 þ sÞ þ C i h2a r R i þ C 0 h0 h
þC
Z
0;Xi
i¼1 tkþ1 2
tk1
2 Z o A dt þ C ot2 0;Xi
tkþ1 tk
i¼1
2 o/ dt: ot
ð5:26Þ
1;Xi
In addition, we have Akþ1 P Akþ1 Ak P Ak 2 h h 6 Ci s 0;Xi i¼1 i¼1 0;curl;Xi Z 2 2 2 2 t kþ1 X X 2 o A dt þ C i rð/kþ1 Ph /kþ1 Þ0;Xi þ Cis ot2 t 0;Xi k i¼1 i¼1 Z tkþ1 2 2 2 2 X X o A 2a 2a 2a kþ1 2 C i hi þ C 0 h0 þ C i hi / þ Cs 6 ot2 dt 1þa;Xi tk 0;X i¼1 i¼1 Z tkþ1 2 2 2 X o A 2a C i h2a ð5:27Þ 6 i þ C 0 h0 þ Cs ot2 dt;
2 X pffiffiffi b kþ1 2 rR
i¼1
2 X
tk
0;X
482
T. Kang, K.I. Kim / Appl. Math. Comput. 170 (2005) 462–484
and 2 2 X X pffiffiffi 2a rðAkþ1 P h Akþ1 Þ2 6 C i h2a i þ C 0 h0 : 0;Xi i¼1
ð5:28Þ
i¼1
Thus, making use of (5.26)–(5.28) and 5.25 becomes ( 2 2 pffiffiffi kþ1 kþ1 k 2 X pffiffiffi s kþ1 r Hh Hh þ rgkþ1 þ s b r R rg h h 0;Xi 2 s 0;Xi i¼1 pffiffiffi k 2 b rgk s r R h
0;Xi
6
2 X i¼1
2 2 ) 1 1 kþ1 k þ pffiffiffi curl Hh pffiffiffi curl Hh l l 0;Xi 0;Xi
pffiffiffi k 2 b rgk C i s2 r R h
0;Xi
þ
2 X
2a C i sh2a i þ C 0 sh0
i¼1
! Z tkþ1 2 2 Z tkþ1 2 o/ o A þ Cs ot dt : ot2 dt þ tk1 tk 1;X 0;X
ð5:29Þ
For the simplicity of the notation, we define R0 ¼
oA0 A0 AðsÞ ; s ot
b 0 ¼ ðA0 ph A0 Þ ðAðsÞ ph AðsÞÞ þ rð/0 Ph /0 Þ þ R0 : R s Summing both hand sides of (5.29) over k = 0,1, . . . , n and then applying the discrete Gronwalls inequality, we obtain 2 X
(
i¼1
2 n X pffiffiffi Hkþ1 Hkh kþ1 h þ rgh s r s
0;Xi
k¼0
2 1 nþ1 þpffiffiffi curl Hh l 0;Xi 6
2 X i¼1
þ
i¼1
0;Xi
C i h2a i
þ
C 0 h2a 0
0;Xi
)
pffiffiffi 0 2 b rg0 C i s r R h
2 X
pffiffiffi nþ1 2 b þ s r R rghnþ1
þ Cs
Z
þ
2 X i¼1
T s
2 1 0 Ci curl H p ffiffiffi h l 0;Xi
2 2 Z T o A dt þ ot2 0 0;X
! 2 o/ dt : ð5:30Þ ot 1;X
T. Kang, K.I. Kim / Appl. Math. Comput. 170 (2005) 462–484
483
Noting g0h ¼ 0 and then using (4.2) and Theorem 5.1, we get the inequalities 2 2 X X pffiffiffi 0 pffiffiffi 0 2 b rg0 Þk2 6 b k sk rð R sk r R 0;Xi h 0;Xi i¼1
i¼1
6
2 X
C i sh2a i
þ
C 0 sh2a 0
þ Cs
2
i¼1
6
2 X
2 2 o A ot2 dt s 0;X
Z
0
2a 2 C i sh2a i þ C 0 sh0 þ Cs :
i¼1
2 2 X 1 pffiffiffi curl H0 h l
6
0;Xi
i¼1
2 X
C i ðkph A0 A0 k20;curl;Xi þ kA0 P h A0 k20;curl;Xi Þ
i¼1
6
2 X
2a C i h2a i þ C 0 h0 ;
i¼1
Hence, by the definitions (4.4) and (4.5), (5.30) becomes 2 X
2
kos Hh þ rgh k‘2 ðL2 ðXi ÞÞ þ
i¼1
2 X
2
kcurl Hh k‘1 ðL2 ðXi ÞÞ
i¼1
6 Cs þ
2 X
2a C i h2a i þ C 0 h0 :
i¼1
Finally, applying the triangle inequality, we have
2 2 2 X X ðos Ah þ r/h Þ oA þ r/ þ kcurlðAh AÞk2‘1 ðL2 ðXi ÞÞ 2 2 ot ‘ ðL ðXi ÞÞ i¼1 i¼1 6 Cs þ
2 X
2a C i h2a i þ C 0 h0 :
i¼1
This completes the proof of this theorem.
h
Remark 5.1. From (2.6) and (2.7), we see that l1 curl Anþ1 and h n Anþ1 nþ1 nþ1 h Ah s r/nþ1 may be replaced by H and E , respectively. Hence, h h h the conclusion of Theorem 5.2 can be rewritten as 2 X
2
kHh Hk‘1 ðL2 ðXi ÞÞ 6 Cs þ
i¼1 2 X i¼1
2 X
2a C i h2a i þ C 0 h0 ;
i¼1
2
kEh Ek‘2 ðL2 ðXi ÞÞ 6 Cs þ
2 X i¼1
2a C i h2a i þ C 0 h0 :
484
T. Kang, K.I. Kim / Appl. Math. Comput. 170 (2005) 462–484
Acknowledgment Supported by Com2MaC-KOSEF, POSTECH BSRI Research fund 2004 and the scientific research fund of the State Administration of Radio, Film & TV of China.
References [1] R. Albanese, G. Rubinacci, Formulation of the eddy-current problem, IEE Proceedings 137 (1990) 16–22. [2] H. Ammari, A. Buffa, J.-C. Ne´de´lec, A justifiction of eddy currents model for the Maxwell equations, SIAM J. Appl. Math. 60 (2000) 1805–1823. [3] Z. Chen, Q. Du, J. Zou, Finite element methods with matching and nonmatching meshes for Maxwell equations with discontinuous coefficients, SIAM J. Numer. Anal. 37 (2000) 1542– 1570. [4] P. Ciarlet, J. Zou, Finite element convergence for the Darwin model to Maxwells equations, RAIRO Math. Modelling Numer. Anal. 31 (1997) 213–249. [5] P. Ciarlet, J. Zou, Fully discrete finite element approaches for time-dependent Maxwells equations, Numerische Mathematik 82 (1999) 193–219. [6] R. Hiptmair, Symmetric coupling for eddy current problems, SIAM J. Numer. Anal. 40 (2002) 41–65. [7] Q. Hu, J. Zou, a non-overlapping domain decomposition method for Maxwells equations in three dimensions, SIAM J. Numer. Anal. 42 (2003) 1682–1708. [8] J. Jin, The Finite Element Method in Eletromagnectics, second ed., John Wiley & Sons, INC., North-Holland, 2002. [9] T. Kang, K.I. Kim, An E-based splitting finite element method for time-dependent eddy current equations, in press. [10] C. Ma, Studies on A/ methods for time-dependent electromagnetic fields (in Chinese), Ph.D. thesis, Chinese Academy of Sciences, 2003. [11] J.-C. Ne´de´lec, Mixed finite elements in R3 , Numer. Math. 35 (1980) 315–341.