A mixed finite element method for computing three-dimensional time-domain electromagnetic fields in strongly inhomogeneous media

(1)

674 IEEE TRANSACTIONS ON MAGNETICS, VOL. 26, NO. 2, MARCH 1990 A MIXED FINITE ELEMENT METHOD FOR COMPUTING THREE-DIMENSIONAL TIME-DOMAIN ELECTROMAGNETIC FIELDS IN STRONGLY

INHOMOGENEOUS MEDIA Gerrit Mur

Delft University of Technology, Faculty of Electrical Engineering, Laboratory-of Electromagnetic Research,

P.O. Box 5031, 2600CA. Delft, The Netherlands, Telephone 31-15-786294. Abstract

A mixed finite-element method is presented for the computer modeling of three-dimensional time- dependent (transient) electromagnetic fields in elec- trically and/or magnetically strongly inhomogeneous

time-invariant media. In the spatial domain

of

compu-

tation both Edge and Cartesian elements are used for the electric field strength as well as the magnetic field strength and, for both expansions, the corresponding program decides locally what type of element to be used for obtaining the best results both as regards computational efficiency and final accuracy. Both implicit and explicit methods can be used for carrying out the integration along the time axis.

Introduction

In the present paper an accurate finite-element method for solving the time-dependent electromagnetic field equations in strongly inhomogeneous media is presented. By formulating the problem in terms of the electric and/or the magnetic field strength in stead

of (vector)potentials we can avoid the loss of accura-

cy that results from the differentiations required for

computing the electric or the magnetic field strength from these potentials. This approach has the addition- al advantage of avoiding problems due to multiple connectedness and the choice of guages for ensuring uniqueness, that are met when using potentials C1.2, 31. In our method the electric and the magnetic field

strengths are computed si~ultaneously, methods of this

type are referred to as mixed methods C41. This method was chosen to be able to treat the electric and the magnetic field strength in a symmetric way for obtaining the best overall accuracy. When using the finite-element method for the computer modeling of three-dimensional electromagnetic fields in strongly inhomogeneous media, it is necessary to use elements that automatically account for the continuity of the tangential field components across interfaces where the relevant constitutive coefficients jump, and allow for a jump in the normal components C5,61. Elements satisfying these conditions are referred to as Edge elements, the classical elements that make all field components continuous across interfaces, and that are accurate only in homogeneous media, are referred to as Cartesian elements. In C61 a finite-element method for computing time-harmonic electromagnetic fields in strongly inhomogeneous media was presented that uses both Edge and Cartesian elements for obtaining optimum results both as regards computational efficiency and as regards accuracy.

As regards the discretization in space, the method presented in this paper was developed along similar lines as the one described in C61. It uses computationally optimum combinations of Cartesian and Edge elements for the expansion of both the electric and the magnetic field strengths. The optimum expansions are obtained by letting the program decide, for each combination of two adjacent domains and for both the expansion of the electric field strength and the expansion of the magnetic field strength, which type of element will, locally, be used. For small differences in electromagnetic medium properties it will use Cartesian elements. For large differences in the electrical properties of the media, Edge elements will be

used for the expansion of the electric field strength;

similarly, for large differences in the magnetic prop-

erties of the media, Edge elements will be used for the expansion of the magnetic field strength. The degree of inhomogeneity above which Edge elements are used, i.e. a strong inhomogeneity, is user-defined and depends on the final accuracy the user aims at. of coupled first-order ordinary differential equations are derived that can be used for computing the evolution in time of the expansion coefficients by either using an implicit or an explicit method.

Using the method of weighted residuals, systems

Cartesian expansion functions and Edge expansion functions

---

The geometrical domain in which the finite- element method is applied, is subdivided into adjoin-

ing tetrahedra. To specify a position in space, we employ the coordinates [x,y,z} with respect to the

Cartesian reference frame with origin 0 and three

mutually perpendicular base vectors (ix=il, iy=i2, i =

i 1 of unit length each. The time coordinate is deno-

ted by t.

The position vectors of the vertices of a partic- ular tetrahedron T are {ro,rl,r2,r3]; the outwardly

directed vectorial areas of the faces of T are (Ao,Al,

A2,A3} C5l and the volume of T is denoted by V. Let rb

be the position vector of the barycenter of T, then the linear functions 9, that equal unity when r

-

ri I (1-0,1,2,3) and are zero at the remaining three vertices of T can be expressed as

3

$l(r) = 114

-

(r

-

rb).Ai/3V. (1)

Introducing the vectors .)'!V and V(E) through i. (j=1,2,3), and

local Cartesian expansion functions W)'(

-

J i. j V(.E) = -A /3V (j==0,1,2,3), the J J j are given by isj

and local the Edge expansion functions ?':U are given

1J by

where the factor a

-

Ir -r

1,

denoting the length

of the edge joining ri and rj, was introduced and for

making W(E) dimensionless and for normalizing it for improving the condition of the final system of equations. Note that both the Cartesian and the Edge ex-

pansion functions are consistently linear functions of

all spatial coordinates. Consistently linear expansion functions are preferred to expansion functions that

are not consistently linear [7,81 because they yield

more accurate results with less computational effort C5I.

magnetic field strength H are locally expanded as

i.j j

1.j

When r E T, the electric field strength E and the

(2)

615

( 4 )

( 5 )

where e .(t) and hi' (t) are time dependent expan-

sion coefficients. In ( 4 ) and (5) the summation over j and j1 runs over values that depend on the type of element. Since, in general, the materials in a given configuration will, as regards their electric and magnetic properties, not have the same degree, or kind, of inhomogeneity, the distribution of Cartesian and Edge elements used for the expansion of the electric field strength will differ from the distribution used for the expansion of the magnetic field strength.

In our notation we have stressed this difference by

using primed subscripts for the expansion of the magnetic field strength. Note, however, that both expansions are related to the same subdivision of the domain of computation into tetrahedra.

The expansion functions given above are all of the local type, i.e. they apply to one tetrahedron only. We now introduce the global expansion functions as the functions that consist of all local expansion functions that are related to the same expansion coeffi- cient. Global Cartesian expansion functions therefore consist of all local Cartesian expansion functions that are related to a Cartesian component in a common vertex of a number of adjoining tetrahedra (the sim-

plicial star [SI of this vertex). Global Edge expan-

sion functions consist of all local Edge expansion functions that are related to an Edge expansion coef- ficient related to a common edge of a number of adjoining tetrahedra (the simplicial star of this edge). From this it follows that the simplicial star (the support) of global Cartesian expansion functions is much larger than the simplicial star of global Edge expansion functions; this fact causes Edge expansion functions to yield more unknowns which makes them computationally less efficient. This larger number of unknowns may, however, be compensated by the increased sparsity of the system matrices, which is another con- sequence of difference in the support of the global expansion functions.

-

T h e 2 t e m of differential equations in the expansion coefficients

i,J ,jl

The time-domain electromagnetic field equations are written in the form

atEE + oE

-

V x H =

-

Je, _{( 6 )}

atpH + V x E =

-

Ke, ( 7 )

where Je and Ke denote external sources of electric and magnetic current, respectively. The system of differential equations in the expansion coefficients is obtained by applying to ( 6 ) and ( 7 ) the method of weighted residuals, where the weighting functions

W(C'E)(r), that are used for weighting (61, and

W(y*E)(r), that are used for weighting (7), are taken

equal to the corresponding expansion functions for E and H, respectively. In this way we obtain

P19 P .qf

where all double summations have, for brevity, been replaced by a single summation sign, and where the double subscripts, indicating the relevant local expansion, have been replaced by a single greek sub-

script related to the global expansions (a and a') and

the global weighting ( 6 and 5'). Assuming the medium properties to be constant in the interior of each tetrahedron and assuming the external source distribu- tions to vary linearly in the interior of each tetrahedron, all spatial integrations implied in ( 8 ) and

( 9 ) can be carried out analytically [ l o ] . For linear media, the resulting system of equations can, in matrix form, be written as

M

a

e + H o e

-

M# =

-

je,

?Jath + ME" =

-

ke. ( 1 1 )

(10) E t

Note that ME, MO and M are square matrices, whereas

I$, and ME are, in general, not square. Together with

the appropriate initial and boundary conditions, ( 1 0 ) and ( 1 1 ) are a coupled system of ordinary linear differential equations from which the evolution in time of the expansion coefficients can be obtained by, for instance, using a single step time marching scheme [ l l ] . In a first attempt to solve the electromagnetic- field problem, however, we have replaced the above fully implicit formulation by a simpler explicit one. Explicit method of solution

For obtaining an explicit formulation of the problem

it is general practice to use a so-called lumping procedure [121. Lumping implies choosing a lower-order integration rule for evaluating the spatial integra-

tions so as to obtain approximate matrices that are

diagonal because of which the formulation becomes ex-

plicit (in our case a diagonalization of the matrices

M E , M and M is required). Rather than choosing lower-order integration rules we will obtain the diagonalization by using modified weighting functions ~ ( ~ ' ~ ? r ) . These modified weighting functions are

i,j

r.)/ai where 6 ( r ) denotes a spatial Dirac delta

function and where ri is taken as the limit towards

the vertex i, taken from the inside of the relevant tetrahedron. The integrations in ( 8 ) and ( 9 ) can be evaluated exactly by using the definition of the spatial Dirac delta function and an explicit formulation

of the problem is obtained. Our approach in deriving

an explicit formulation can be interpreted as replac- ing the Galerkin-type weighting of the approximated form of (6) and (7) over the relevant simplicial stars by collocation, i.e. by requiring the approximated

form of these equations to be satisfied exactly at, or

in the immediate vicinity of, the relevant vertices.

The explicit integration in time is carried out by

evaluating e at t = to + nAt, n = 1,2,

...

and evaluating h at t = to + (n-ll2)At. The differentiations with respect to time can be carried out as central differences and we obtain a result that can be written as

(3)

676

I n a l l computations t h e mesh c o n s i s t s of 8*8*8 cubes

of e q u a l s i z e t h a t a r e each s u b d i v i d e d i n t o 6 tetrahedra. The comvutatlons are c a r r i e d o u t f o r t h e en+l = (HE + A t H0/2)-1((ME

-

A t H0/2)en

(12) time i n t e r v a l to = 0 I t 5 2.5*10-9s = tend, u s i n g 50

-

A t I4-l (MEen + ken). (13) s t e p s i n time (At = 0.5*10-10s). The i n i t i a l condl-

t i o n s f o r the magnetic f i e l d s t r e n g t h a t t = -At12 are t a k e n from ( 1 4 ) , t h e i n i t i a l c o n d i t i o n s f o r t h e e l e c - t r i c f i e l d s t r e n g t h a t t = 0 a r e computed u s i n g ana- hn+l 12 = hn-l /2

Note t h a t the i n i t i a l c o n d i t i o n s eo and h-1,2 are r e q u i r e d a t t=to and t a t 0 - A t , r e s p e c t i v e l y .

Numerical r e s u l t s

To d e m o n s t r a t e t h e a c c u r a c y of our method, we have a p p l i e d t h e FEMAXT code, developed by u s u s i n g t h e methods d e s c r i b e d i n the p r e s e n t paper, t o a problem f o r which t h e s o l u t i o n is known a n a l y t i c a l l y . T h i s problem was first s o l v e d by u s i n g (12) and (13) under t h e c o n d i t i o n t h a t Edge expansion f u n c t i o n s a r e used f o r the expansion of t h e e l e c t r i c f i e l d s t r e n g t h when t h e r e l a t i v e c o n t r a s t i n t h e numerical v a l u e s of E and/or a i n two a d j a c e n t tetrahedra exceeds 10%. S i m i l a r l y , Edge expansion f u n c t i o n s are used for t h e expansion of t h e magnetic f i e l d s t r e n g t h when t h e r e l a t i v e c o n t r a s t i n the numerical v a l u e s of 11 i n two a d j a c e n t t e t r a h e d r a exceeds 10%. C a r t e s i a n e l e m e n t s a r e used when c o n t r a s t s less t h a n 10% are found. For comparison, t h e same problem i s a l s o s o l v e d by u s i n g C a r t e s i a n elements only. The c o n f i g u r a t i o n used f o r t e s t i n g t h e method is t h e c u b i c , s o u r c e - f r e e , r e g i o n D (-0.5

<

x

<

0.5~1, 0

<

y

<

lm, 0

<

z

<

lm} c o n s i s t i n g

of two homogeneous p a r t s with d i f f e r e n t medium p r o p e r t i e s v i z . D, (-0.5

<

x

<

Om, 0

<

y

<

I m , 0

<

z

<

Im) w i t h the p r o p e r t i e s of a vacuum (El=EO,

q q 0 ,

ol-O} and D2=D\D1 with the medium p r o p e r t i e s [ ~ 2 = 1 0 ~ o , u2-u0, a2=0.001S/m} ( t h i s c h o i c e a p p l i e s t o t h e r e a l -

l y t i c a l e x p r e s s i o n s f o r t h e e l e c t r i c f i e l d s t r e n g t h . The l a t t e r e x p r e s s i o n s are a l s o used f o r p r e s c r i b i n g t h e boundary c o n d i t i o n s f o r t h e e l e c t r i c f i e l d s t r e n g t h a t t h e o u t e r boundary o f t h e domain of

computation, t h e magnetic f i e l d s t r e n g t h a t the o u t e r boundary remains u n s p e c i f i e d .

I n Fig. 1 t h e l o c a l r e l a t i v e e r r o r i n t h e s o l u t i o n a t t = tend is p l o t t e d f o r a given c r o s s - s e c t i o n ( z = 0.5001) of the c o n f i g u r a t i o n , and f o r t h e case Edge e l e m e n t s are used a l o n g t h e p l a n e of d i s c o n t i n u i t y . I n F i g . 2 r e s u l t s a r e p l o t t e d for e x a c t l y t h e same problem t h a t was now s o l v e d u s i n g C a r t e s i a n e l e m e n t s only. I n t h e r e s u l t s p r e s e n t e d , t h e local r e l a t i v e error i n t h e s o l u t i o n i s d e f i n e d as t h e a b s o l u t e v a l u e

of t h e e r r o r a t a g i v e n p o i n t i n the domain of

computation and a t a c e r t a i n moment of time, d i v i d e d by t h e maximum a b s o l u t e v a l u e of the e x a c t s o l u t i o n anywhere i n t h e domain of computation and a t t h e same moment of time. I n T a b l e 1 t h e r o o t mean s q u a r e ( R M S ) e r r o r of a l l f i e l d components are p r e s e n t e d f o r t h e same problems. The RMS e r r o r of, f o r i n s t a n c e , Ex is d e f i n e d as

( 1 5 ) i s t i c case of a p l a n e i n t e r f a c e between a l o s s y d i e -

lectric and a

t h e c o n f i g u r a t i o n is chosen as a time-harmonic p l a n e wave i n D 1 , t h a t has i t s magnetic f i e l d p o l a r i z e d

The magnetic field strength in where S d e n o t e s _{mentioned above.}t h e _{A s}c r o s s - s e c t i o n of _{r e g a r d s the computation}t h e c o n f i g u r a t i o n _times_for the two problems, we n o t e t h a t t h e y were about t h e same. t h e e f f e c t s of t h e larger number of unknowns

-

a l o n g t h e z-axis. Because o f t h e d i s c o n t i n u i t y a t x =

0 . w e a l s o have a r e f l e c t e d wave i n D. and a e l i m i n a t e d by t h e e f f e c t s of t h e i n c r e a s e d s p a r s i t y of

t h a t are r e q u i r e d when u s i n g Edge e l e m e n t s being

I t r a n s m i t t e d wave i n D2. A f i e l d of t h i s t y p e can, u s i n g complex a r i t h m e t i c , be r e p r e s e n t e d as Hx = O s HY = O * HZ = Re(exp(jwt-Yi.r) + R exp(jwt+Yr.r)) when x

<

0, and t HZ = Re(T exp(jwt-Y . r ) ) when x 5 0, ( 1 4 )

where Re( ) i n d i c a t e s t h a t t h e real p a r t has t o be taken. W r i t i n g Yi = Y i s i , where si d e n o t e s the u n i t v e c t o r i n t h e d i r e c t i o n of p r o p a g a t i o n of t h e i n c i d e n t wave, w e have chosen s: = si = 2-'12 and s: = 0. W i t h

Y

t h i s c h o i c e , ( 1 4 ) r e p r e s e n t s a p l a n e wave i n D1 tra-

t h e m a t r i c e s . From t h e p l o t s it is observed t h a t , c o n t r a r y t o Edge e l e m e n t s , C a r t e s i a n e l e m e n t s produce e x t r e m e l y l a r g e e r r o r s i n a r e g i o n n e a r t h e i n t e r f a c e

of d i s c o n t i n u i t y . From o t h e r r e s u l t s it was observed t h a t t h e w i d t h of t h e r e g i o n i n which t h e r e s u l t s show l a r g e e r r o r s i n c r e a s e s i n time. The f a c t t h a t t h e u s e

of Edge e l e m e n t s g r e a t l y improves t h e accuracy i s confirmed by the r e s u l t s f o r t h e RMS-errors g i v e n i n Table 1 . Regarding t h e s t a b i l i t y of t h e p r e s e n t e x p l i c i t method we n o t e t h a t I t is n o t u n c o n d i t i o n a l l y s t a b l e . Numerical experiments show t h a t Edge e l e m e n t s are more s e n s i t i v e t o i n s t a b i l i t i e s t h a n C a r t e s i a n elements. When problems have t o be s o l v e d f o r l a r g e i n t e r v a l s of time, i t w i l l be n e c e s s a r y t o u s e an i m p l i c i t method for s o l v i n g (IO) and (111, i n o r d e r t o cope w i t h t h e s t a b i l i t y problem.

A l l computations have been c a r r i e d o u t a t a VAX- S t a t i o n 2000, u s i n g t h e SEPRAN f i n i t e element package C131.

T a b l e 1 . R e l a t i v e e r r o r s i n $ when u s i n g C a r t e s i a n v e l l i n g i n the x-y p l a n e and having a n a n g l e of i n c i -

dence of 4 5 O on t h e i n t e r f a c e . Using ( 6 ) and ( 7 ) to- g e t h e r w i t h t h e c o n t i n u i t y c o n d i t i o n s a t x = 0, ex- p r e s s i o n s f o r t h e unknown q u a n t i t i e s i n ( 1 4 ) can be e l e m e n t s t o g e t h e r w i t h Edge e l e m e n t s o r u s i n g C a r t e s i a n e l e m e n t s only. --^-^---^---_^_____I___^ RMS-error i n % o b t a i n e d e a s i l y . I n t h e same way e x p r e s s i o n s f o r E

and E can be o b t a i n e d , EZ

=

0. Due t o t h e d i s c o n t i n - u i t y i n t h e medium p r o p e r t i e s Ex w i l l be d i s c o n t i n u o u s Elements

_----.---___._^___

HZ H E X E Y E Z H X Y Y used

---

---- ---.---

a t x = 0, the remaining f i e l d components b e i n g c o n t i n - Cart.+ Edge 2-20 3 - 3 6 0.38 0.45 0.59 3.30

uous throughout the domain of computation. The f i e l d C a r t . o n l y 40.7 4.07 0.04 0.12 0.03 3.98

(4)

671

Conclusion

---

We have shown that a mixed finite element method that uses Edge elements together with Cartesian elements for computing three-dimensional time-domain electromagnetic fields in strongly inhomogeneous media, yields results that are, as regards accuracy, superior to results that are obtained using Cartesian elements only. This conclusion does not only apply to the simple configuration for which we have tested the method, but also to configurations containing

arbitrary inhomogeneities. The reason for this is that

in the method a decision is made, for each combination of adjacent tetrahedra, what type of elements should be used for obtaining accurate results. Thus the ability of the method to model arbitrary inhomogene-

ities is limited only by the ability of the mesh-

generator used to accurately model the configuration geometrically. As regards stability an implicit method

is required for being able to solve problems of this

type for large intervals of time, for relatively short intervals of time the explicit method can be employed.

Finally we note that our method for computing electric

and/or magnetic field strengths is more accurate than

the often used (vector)potential methods because, for computing the electric and/or the magnetic field strength from these potentials, the latter require a numerical differentiation that causes a large loss of accuracy in the final results.

References

[I] C. R. I. Emson and C. W. Trowbridge, "Transient 3D eddy currents using modified magnetic vector potentials and magnatic scalar potential",

Transactions on Magnetics, Vol. MAG-24, No. 1,

January 1988, pp. 86-89.

C21 C. R. I. Emson and J. Simkin, "An optimal method for 3D eddy currents", IEEE Transactions on Magnetics, Vol. MAG19, No. 6, November 1983, pp. 2450-2452.

[3] P. J. Leonard and D. Rodger, "Finite element scheme for transient 3D eddy currents",

IEEE-

--

Transactions on Magnetics, Vol. MAG24, NO. 1 , January 1988, pp. 90-93.

[ll] 0. C. Zienkiewicz and R. L. Taylor, The finite

element method, Vol. 1, Basic formulation and

linear problems, London, McGraw-Hill Book Company. 1989, p. 319.

[5] G. Mur-and A. T. de Hoop, "A finite element method for computing three-dimensional

electromagnetic fields in inhomogeneous media," -IEEE Transactions on Magnetics, -_-I vol. MAG-21, No. 6, November 1985, pp. 2188-21 91.

[6] G. Mur, "Optimum choice of finite elements for comuuting three-dimensional electromagnetic fieids i n inhomogeneous media," IEEE Transactions on Magnetics, Vol. MAG-24, No. 1 , January 1988, DD. 330-333.

C71 M. Li Barton and 2. J. Cendes, "New vector finite elements for three-dimensional magnetic field computation," J. Appl. Phys,, vol. 61, No. 8, April 1987, pp. 3919-3921.

[8] A. Bossavit and J. C. VBritB, "The 'Trifou' code: Solving the 3-D Eddy-currents problem by using H as state variable," IEEE Transactions on Magnet- ics, Vol. MAG-19, NO. 6, 1983, pp. 2465-2470. C91 T L . Naber, Topological methods in Euclidean spaces, Cambridge, Cambridge University Press, 1980, D. 62.

1101 A. H; Stroud, Approximate integration of multiple integrals, Englewood Cliffs, Prentice-Hall Inc., 1971. ChaDter 7. _~.

[ll] 0. C. Zienkiewicz, W. L. Wood and N. W. Hine, "A unified set of single step algorithms, Part 1 :

General formulation and application," In_t. j.

numer. methods eng., Vol. 20, (19841, pp. 1529-

1552.

Computational aspects, Vo1.3, Englewood Cliffs, Prentice-Hall Inc., 1984, p.273.

cl21 G. F. Carey and J. T. Oden, Finite elements,

Cl31 A. Segal, SEPRAN, Sepra Analyes, User Manual, Leidschendam, The Netherlands, Sepra, 1984.

1 y1

'$

>. 0.8 0.6 0.4 0.2 0

ror i n Ex i n plane z= O.500100Et00 i n % of problem: femaxld-l

)

/

\

0

i

0 I CONTOUR KEY

k Z l

0 . 3 4 1 E t 0 1

I

il

::::::::I

-0.4 - 0 . 2 0 0 . 2 x-axis 0 . 4

Fig.1. Contour plot of the local relative error in Ex,

Edge elements, RMS(E ) = 2.28, 4137 unknowns.

1 lo .I 1 >r 0.8 0.6 0.4 0.2 0

:ror i n Ex i n plane z- IEtOI

I

?,% of problem: fenaxldl0 0 . 4 3 2 E t 0 2 0 . 6 1 9 E t 0 2 I 1 0 . 2 x - a x i s 0 . 4 Fig.2. Contour plot of the local relative error in Ex