Naturwissenschaften (2001) 88:183–192 DOI 10.1007/s001140100215

R E V I E W A RT I C L E

Stefan Kurz · Ulrich Becker · Harald Maisch

Dynamic simulation of electromechanical systems: from Maxwell’s theory to common-rail diesel injection

Published online: 10 May 2001 © Springer-Verlag 2001

Abstract This paper describes the state-of-the-art of dynamic simulation of electromechanical systems. Electromechanical systems can be split into electromagnetic and mechanical subsystems, which are described by Maxwell’s equations and by Newton’s law, respectively. Since such systems contain moving parts, the concepts of Lorentz and Galilean relativity are briefly addressed. The laws of physics are formulated in terms of (partial) differential equations. Numerical methods ultimately aim at linear systems of equations, which can be solved efficiently on digital computers. The various discretization methods for performing this task are discussed. Special emphasis is placed on domain decomposition as a framework for the coupling of different numerical methods such as the finite element method and the boundary element method. The paper concludes with descriptions of some applications of industrial relevance: a high performance injection valve and an electromechanical relay.

Introduction The nineteenth century produced many prominent scientists. Possibly the greatest of these was the physicist, James Clerk Maxwell. His work formed the basis of modern physics when he united electricity and magnetism into the concept of the electromagnetic field. This paved the way for the theoretical description of many electromagnetic and electromechanical products which are still manufactured in industry nowadays. In 1927, the first mechanical injection pump for diesel engines was manufactured by the Robert Bosch company. Since then, a tremendous amount of innovation has taken place. Today’s increasing demands for lower fuel consumption, cleaner exhaust gases and quieter running S. Kurz (✉) · U. Becker · H. Maisch Robert Bosch GmbH Stuttgart, Postfach 10 60 50, 70049 Stuttgart, Germany e-mail: [email protected] Tel.: +49-711-8117222, Fax: +49-711-8116301

of the diesel engine can no longer be met by using mechanically controlled injection systems. Very high injection pressures, as well as accurate control of the injection process, are required. These requirements are met very well by the common-rail fuel injection system, as shown in Fig. 1. The injection pressure is generated by a highpressure pump. The highly compressed fuel is stored in a tubular rail awaiting injection. The appropriate timing of injection and the amount of fuel required are calculated by the control unit and the fuel is injected into each cylinder by a magnetic actuated valve. Extremely precise fuel metering is achieved by minimal valve switching times, which imposes extraordinarily high performance requirements on the valves. Modeling and simulation are key elements in assuring the fast and successful design of electromechanical systems. On the geometric level, electromagnetic field computation is an indispensable tool in the design process. Since many practical problems can only be solved incompletely or not at all by analytical methods, numerical methods have gained increasing acceptance during recent decades. The successful solution of electromagnetic field problems by numerical methods has partly been enabled by powerful computer hardware. Fast processors

Fig. 1 The common-rail injection system: from right to left the control unit, the high pressure pump, the tubular rail and the four injectors are shown

184

and large memory capacities allow the efficient solution of large linear systems of equations, while the improved graphics capabilities are essential for creating userfriendly software. The computer implementation and comprehensive application of the theory arising from the classical physics of the nineteenth century have been pushed ahead by improvements in the numerical algorithms themselves. The starting point for the detailed discussion in the following section is the physical theory centered around Maxwell’s equations. The next section is devoted to the numerical methods, where special emphasis is placed on domain decomposition as a framework for the coupling of different methods such as the finite element method and the boundary element method. The numerical methods for the dynamic simulation of electromechanical systems presented here have been successfully applied to a variety of problems, the high performance magnetic valve of the common-rail injection system being one of them. The final section briefly describes some results for a magnetic valve and an electromechanical relay.

where the operators E and M describe the dynamic characteristics of the electromagnetic and mechanical subsystems. qE, qM are the electromagnetic and mechanical state variables, respectively, and the dot denotes the time derivative. sE is the electromagnetic source term (e.g., prescribed currents and voltages, permanent magnetizations) and sM the mechanical source term (e.g., spring pre-stress, gravitational forces). The mutual coupling of both subsystems stems from the fact that the operator E depends on qM and the operator M on qE: each change in the mechanical variables qM causes a change in the electromagnetic variables qE and vice versa. The coupled solution of Eqs. 1 and 2 (whether analytical or numerical) is a difficult task in general, because both equations exhibit nonlinear behavior, for example, due to iron saturation or the quadratic dependence of the electromagnetic forces on the fields. A numerical solution of Eqs. 1 and 2 inevitably requires a time discretization, where the time is subdivided into individual, sufficiently small time-steps ∆t: tn+1=tn+∆t, n=0, 1, 2, ...

Physical theory Electromagnetic and mechanical subsystems An electromechanical system can always be separated into electromagnetic and mechanical subsystems which interact with each other (see Fig. 2). The electromagnetic subsystem is governed by Maxwell’s equations at the low frequency limit. The mechanical subsystem obeys Newton’s law. The interaction between both subsystems is twofold: firstly, the electromagnetic forces and torques exerted on the mechanical subsystem cause deformation and acceleration; secondly, the motion affects the electromagnetic subsystem by motion-induced electromotive forces (EMF) and by a time-variant geometry. The geometry being time-variant means that the spatial distribution of material parameters as well as the location of material interfaces are changeable. Only moving rigid bodies are treated in this paper. A formal description of the electromechanical system in Fig. 2 is given by the equation: E(qE, q˙E, qM)=sE

(1)

M(qM, q˙M, qE)=sM

(2)

Fig. 2 Decomposition of an electromechanical system into a electromagnetic and a mechanical subsystem which interact with each other

(3)

A time integration scheme must be able to obtain the valt+∆t t t ues qt+∆t E , qM at t+∆t from the values qE, qM at t. Different strategies for this task can be found in the literature. The “strong” electromagnetic coupling solves Eqs. 1 and 2 simultaneously (Ren and Razek 1994; Nysveen and Nilssen 1997), where difficulties arise from the abovementioned nonlinearities. A simplification can be achieved by means of predictor–corrector type methods. These methods try to predict the state variables at t+∆t with a simple (low order explicit) method. Starting from the predicted values, the equations are solved with a more sophisticated (higher order implicit) integration scheme. Both steps are repeated iteratively until some accuracy criterion is met. If the predicted values q¯M of the electromagnetic variables for the solution of the mechanical equation, an effective decoupling can be achieved. The corrector step can be omitted if the timestep is chosen to be so small that neither the geometry nor the electromagnetic forces change significantly during ∆t (Nicolet et al. 1992): the mechanical variables q¯M will be regarded as constant for the solution of the electromagnetic equation and vice versa. This algorithm is called “weak” electromagnetic coupling (Brauer and Ruehl 1994; Henrotte et al. 1994; Hsieh 1995) and is depicted in Fig. 3. First, Eq. 1 will be solved with the moving bodies bet . The electromaging fixed in the positions given by qM netic forces and torques are computed from the resulting state qEt+∆t and regarded as constant for the solution of Eq. 2. Finally, the position of the moving bodies is upt+∆t . Completely independent algodated according to qM rithms can be employed for both subsystems. The communication between them relies on the exchange of force and position information only. Accordingly, the rest of this section is organized as follows: first the electromagnetic subsystem described

185

n · [B ]=0

(11)  where ρs and j s are possible surface charge and current densities. The brackets denote the jump of their argument when a surface of discontinuity is passed in the direction of the normal n. The conditions expressed in Eqs. 9 and 10 are valid for interfaces at rest and have to be supplemented by additional terms in case of moving interfaces (Costen and Adamson 1965). The system of Maxwell’s equations (Eqs. 4, 5, 6, 7) cannot be solved until the fields D and H are known in terms of E and B . These connections are known as constitutive relations: D =D (E ,B ) (12)

H =H (E ,B )

(13)

In addition, for conducting media there is the generalized Ohm’s law j=j(E ,B ) (14)

Fig. 3 “Weak” electromechanical coupling by means of alternating solution of the electromagnetic and mechanical equations

by Maxwell’s equations is presented, then the mechanical subsystem is examined. The electromagnetic subsystem In 1873, J.C. Maxwell published his famous “Treatise” (Maxwell 1954) where he summarized all classical electromagnetic phenomena in a set of four equations, now known as Maxwell’s equations: r (4) div D = ρ r r r ∂D (5) curl H – =j ∂t r r (6) curl E + ∂B = 0 ∂t r (7) div B = 0 where D is the electric displacement, H the magnetic field, E the electric field and B the magnetic induction. ρ and j are the free charge and current densities in the considered medium (Jackson 1975). From Maxwell’s Eqs. 4–7 the following interface conditions can be derived: (8) n · [D ]=ρs n × [H ]=js n × [E ]=0

(9) (10)

The general notation used in Eqs. 1–14 indicates that these relations are not necessarily simple. Anisotropic, nonlinear and even hysteretic behavior may occur. However, for the description of electromechanical systems the following simple form of the constitutive relations shown in Eqs. 12–14 is sufficient in many cases: (15) D =ε0E H =ν (x,|B |)B j=σ(x)E

(16) (17)

where ε0 is the permittivity of free space, ν the reluctivity, and σ the conductivity of the considered material point. The reluctivity ν is the reciprocal of the permeability µ and the notation of Eq. 16 means that ν might depend on the position x and the magnetic induction. This allows us to take into account inhomogeneous materials and magnetic saturation. ν(|B |) can be computed easily from the usual B(H) curves. Sometimes it is convenient to split Eq. 16 according to: (18) H =ν 0B–M , M =(ν 0–ν (x,|B |))B where ν0=1/µ0 is the reluctivity of free space and M is the magnetization. Equation 17 is the differential form of Ohm’s law for bodies at rest. The homogeneous equations (Eqs. 6 and 7) can be solved by introducing the electric scalar potential ϕ and the magnetic vector potential A according to: r r ∂ A E = – gradϕ – (19) ∂t B =curlA (20) The potentials A and ϕ are not uniquely defined by Eqs. 19 and 20. A pair of potentials (A,ϕ) can always be chosen such that r ∂ϕ =0 divA + 12 (21) c ∂t

186

holds, where c=1/√ε0µ0≈3·108m/s is the velocity of light in vacuum. The relation expressed by Eq. 21 is called the Lorentz condition. Using Maxwell’s equations (Eqs. 4 and 5), the constitutive relations (Eqs. 15 and 18), the potential ansatz (Eqs. 19 and 20), and the Lorentz condition (Eq. 21), the following equations can be derived: r r r 2 = – µ0 jtot ∆ A – 12 ∂ A (22) 2 c ∂t

ρ ∂2ϕ (23) ∆ϕ – 12 2 = – ε0 c ∂t In Eq. 22, jtot denotes the total current density which is the sum of the free current density j and the magnetization current density jmag. r r r r r jtot = j + jmag = j + curl M (24) r r r r = j + curl(ν 0 – ν ( x, B ))curl A Equations 22 and 23 are the inhomogeneous wave equations for the potentials A and ϕ and completely equivalent to Maxwell’s equations (Eqs. 4, 5, 6, 7). The currents and charges on the right sides of Eqs. 22 and 23 must not be prescribed independently. They have to obey the principle of charge conservation: r ∂ρ (25) div j + =0 ∂t which follows from Maxwell's equations by taking the divergence of Eq. 5 and using Eq. 4. Equation 25 implies the interface condition r r (26) n⋅ j = 0

[]

Interface conditions for the potentials follow from Eqs. 19, 20 and 21 and read r A =0 (27)

[]

[ϕ ] = 0

(28) i.e. the potentials have to be continuous across interfaces. For technical electromechanical applications, the electromagnetic phenomena are usually low frequency so that wave propagation can be neglected. This approximation of electromagnetism goes under the name quasistatic (Penfield and Haus 1967). A change of the current and charge distribution propagates instantaneously throughout the whole space, which corresponds to the limiting case c→∞. This approximation is justified if the dimensions of the considered apparatus are much smaller than a wavelength of a uniform plane wave at the frequency in question. The most general quasi-static system can be decomposed into a quasi-static electric field system and a quasi-static magnetic field system, which are complementary to each other, the latter also sometimes being called quasistationary. The splitting of a general quasi-static system separates the capacitive effects (which are ascribed to the electric field system) from the inductive effects (which are ascribed to the magnetic field system). From the application’s point of

Fig. 4 Sample problem which consists of an exciting coil driven by the current i0, a conducting magnetic body Ωc with boundary Γ=∂Ωc and the surrounding air domain Ωa

view, actuators which rely on magnetic forces (e.g. magnetic valves) can be described by a quasi-static magnetic field system. Conversely, actuators based on electric forces (e.g. piezo-electric devices) should be described by a quasi-static electric field system. Quasistatic electric field systems will not be discussed further in this paper. The quasi-static magnetic field system can be introduced formally by ε0→0, which implies D →0, ρ→0, c→∞. The displacement current caused by a time-variant electric field is disregarded. This kind of problem is also called magnetic diffusion or an eddy current problem. What remains from the general theory outlined above can be explained best by means of a simple model problem as shown in Fig. 4. It consists of an exciting coil driven by the current i0, a conducting magnetic body Ωc with boundary Γ=∂Ωc and the surrounding air domain Ωa. In the domain Ωc quasi-static Maxwell’s equations are: r r curl H = j (29) r ∂Br (30) curl E + =0 ∂t r (31) div B = 0 with the constitutive relations r r r r r r H = ν ( x, B ) B = ν 0 B – M r r j = σE

(32) (33)

where a constant conductivity σ throughout Ωc has been assumed for simplicity. In the external region Ωa the field equations are r r (34) curl H = j0 r (35) div B = 0 where j0 denotes the impressed exciting current density in the coil and r r H = ν0 B (36) The problem is compounded by regularity conditions at infinity and the interface conditions described in Eqs. 9, 10, 11 and 26 on Γ. In the usual case of finite conductivity σ and continuous time dependency there are no surface currents, hence js=0 in Eq. 9. The calculation of the electric field E

187

in the external region Ωa is not required to determine the eddy currents in Ωc (Albanese and Rubinacci 1990). The eddy current problem can be formulated in many different ways, for instance in terms of the magnetic or the electric field variables or by introducing various potentials (Kriezis et al. 1992). If the potentials A and ϕ are used in connection with the Coulomb condition div A =0 [which is the limiting case of the Lorentz condition (Eq. 21) for c→∞] one receives from Eqs. 29–36 a set of vector and scalar Poisson equations: rr r ∆ A = – µ0 jtot in Ωc (37)

∆ϕ = 0 in Ωc rr r ∆ A = – µ0 j0 in Ωa

(38) (39)

Equations 37 and 38 are the limiting case of the wave equations (Eqs. 22 and 23) for ε0→0. It can be shown that ρ/ε0→0 in Eq. 23 for constant conductivity σ, but the formulation can be easily extended to the nonconstant case σ(x) (Kurz et al. 1999). What are the interface conditions for the potentials A and ϕ? Due to the potential ansatz (Eqs. 19 and 20), the homogeneous Maxwell’s equations are fulfilled automatically and so are the interface conditions (Eqs. 10 and 11). The remaining interface conditions are as Eqs. 9 and 26. Condition 9 gives a constraint on curl A and has to be supplemented by a constraint on div A to ensure a unique solution (Biro and Preis 1989; Albertz et al. 1996) resulting in r r r r (40) n × ν curl A + n ν 0 div A = 0 on Γ  Condition 26 in connection with j =0 in Ωa ensures that no eddy current can flow out of the conducting domain. When the A–ϕ formulation (Eqs. 37, 38, 39, 40) is applied, the potentials play the role of the electromagnetic state variables qE. The equations for the eddy current problem have been derived so far without taking account of motion. The theory of electrodynamics for moving media was formulated by Minkowski (1908) 3 years after Einstein’s enunciation of the special theory of relativity (Einstein 1905). It is well known that Maxwell’s equations (Eqs. 4, 5, 6, 7) are covariant under Lorentz transformation. This means that the laws of electrodynamics have the same form in all inertial frames. When dealing with eddy current problems the propagation velocity c is infinitely high and therefore the Lorentz transformation reduces to the classical Galilean transformation. Quasi-static Maxwell’s equations (Eqs. 29, 30, 31) together with the constitutive relations (Eqs. 32 and 33) are subject to Galilean relativity (Bossavit 1981). For the analysis of eddy current problems with moving bodies, there is freedom in the choice of the reference frame (Landau and Lifschitz 1984). If there is only one conducting body, it is thus advisable to use its rest frame for the solution of the electromagnetic problem (Muramatsu et al. 1992). Because the body is at rest in this frame by definition, all terms that depend explicitly on the velocity can be dropped from the equations. For

[

] [

]

this reason, it was not necessary to take motion into account. In a general situation with more than one conducting body, it is possible to use as many different frames as necessary to describe each body in its rest frame. The individual rest frame descriptions can be coupled together by “frame hopping” using the Galilean transformations and the interface conditions (Kurz et al. 1999). The mechanical subsystem Only rigid body motions are considered here. The most general motion of a rigid body can be broken up into a translation of the center of gravity C associated with three degrees of freedom, and a rotation with respect to C associated with another three degrees of freedom. The translational part is governed by Newton’s law: r r r r r (41) + F , r˙ = v mv˙ = F mech

mag

where m is the mass and r and v¯ are location and velocity of C, respectively. The total force in Eq. 41 has been split into a mechanical force Fmech (e.g., spring or gravitational force), which is internal to the mechanical subsystem and a magnetic force Fmag which acts across the interface of the subsystems (see Fig. 2). The rotational part is governed by Euler’s equation: rr r˙ r rr r r r (42) J ⋅ ω + ω × J ⋅ ω = Tmech + Tmag where J is the tensor of inertia, ω the angular velocity, Tmech the mechanical and Tmag the magnetic torque with respect to C. All quantities in Eq. 42 have to be expressed in components of the body’s rest frame. From Newton’s law (Eq. 41) we can solve directly for the location vector r, but to conclude the instantaneous position of a body from Euler’s equation (Eq. 42) an appropriate set of parameters has to be introduced. A conventional approach is to use Euler angles (Ψ,Θ,Φ) (Greenwood 1988). The union of vectors r,v,ω, and Euler angles (Ψ,Θ,Φ) of the moving bodies can be identified as mechanical state variables qM. The equations of motion (Eqs. 41–42) mentioned so far do not yet constitute the complete description of the mechanical subsystem. They have to be complemented by force and torque laws of the form: r r r r (43) Fmech = Fmech (qM ), Tmech = Tmech (qM ) and by magnetic coupling forces and torques: r r r r Fmag = Fmag (qE ), Tmag = Tmag (qE )

(44)

From the practical point of view, Eq. 43 represents spring characteristics, damping forces, and so on. The magnetic forces and torques can be obtained by the Maxwell stress tensor, which is related to the law of momentum conservation for the electromagnetic field (Jackson 1975). With this method, a normal force density fn and a tangential force density ft are calculated on the boundary Γ of the considered part. The net force and torque are then computed by integration over that boundary.

188

Numerical methods Overview The starting point for the application of numerical methods is a physical description of the electromagnetic subsystem, which is usually given by partial differential equations (e.g., Eqs. 37, 38, 39). All numerical methods ultimately aim at the solution of linear systems which can be done efficiently on a computer. Various steps are necessary to accomplish this task: 1. A time discretization is introduced and the time derivative of the magnetic vector potential in Eq. 19 (which is the only time derivative in the A –ϕ formulation) is replaced by finite differences in time. The simplest suitable method is the backward Euler method, where the time derivative is discretized as: r r r ∂A t ≈ A(tn +1 ) – A(tn ) (45) ∆t ∂t n +1 2. The remaining spatial derivatives have to be tackled by some space discretization scheme. Various methods are known and a vast amount of literature exists. Only the most important methods will be discussed in this paper. The finite element method (FEM) appears to be the most powerful and flexible numerical method (Davies 1980; Dhatt and Touzot 1984). ● The boundary element method (BEM) is the most recent of the major numerical methods (Brebbia 1984; Brebbia et al. 1984). It is very well suited for many electromagnetic problems. ●

3. At this point we are left with a nonlinear system of equations to be solved at each time-step. The nonlinearity is a consequence of the magnetic saturation and demands for iterative solution. To this end, the Newton-Raphson method (NR) (Shyamkumar and Cendes 1988) is very common. With the help of the NR method, the nonlinear problem is replaced by a series of linear problems which are obtained by linearization in the neighborhood of the respective approximate solution. 4. Finally a linear system has to be solved within each nonlinear iteration step. This can be done directly (e.g., by Gaussian elimination or matrix factorization) or once again iteratively (Axelsson 1996).

technical electromechanical problems. There are always conducting and/or magnetic bodies and the surrounding space filled with air, which is not conducting or magnetic. An example of DD of this kind can be seen in Fig. 4, where the whole space has been divided into two subdomains; a domain Ωc which contains conducting and/or magnetic media and a domain Ωa which represents the surrounding space. Only impressed currents (like the exciting current i0) are admissible there. Both subdomains are connected by the common boundary Γ, where the interface conditions in Eqs. 27 and 40 apply. The subdomain Ωa stretches to infinity, where the magnetic vector potential is assumed to decay asymptotically ~1/r. The principle of DD gives the freedom to choose the most appropriate discretization method for each subdomain separately. A combination of FEM and BEM, the BEM–FEM coupling, has proven to be of great advantage for electromechanical systems (Fetzer 1995; Kurz 1998) and has been used for the applications presented later. Inhomogeneous and nonlinear media as well as magnetic saturation can be easily treated within the framework of the FEM. For this reason, the subdomain Ωc is discretized by the FEM, that is further subdivision into individual small elements. Unbounded domains, however, cannot be handled naturally by the FEM, because the discretization is restricted to a limited number of elements. In contrast, the BEM is preferable for homogeneous linear regions which may extend to infinity, the asymptotic condition being fulfilled exactly. Therefore the subdomain Ωa is treated by the BEM which requires a discretization of the boundary Γ only. No discretization for the interior of the subdomain Ωa is necessary. This simplifies the pre-processing considerably and, more importantly, avoids difficulties with the continuously changing positions of the moving parts. A discretization of electromechanical devices with the FEM only would not be easy. Even small displacements would cause a distortion of the mesh, while large displacements would require frequent remeshing or other tricky features, such as “sliding elements” (Tani et al. 1998). With the BEM–FEM approach, these drawbacks can be circumvented because the structure of the mesh does not change during the motion. The solution of the overall problem can be obtained efficiently by iteration. The basic idea is to match the solutions in both subdomains with respect to the interface conditions (Eqs. 27 and 40). Usually, one or both conditions is inherent to the solution strategy, while the other one has to be enforced by the iteration. If both conditions are met sufficiently well in the sense of a suitable error criterion, the overall problem is regarded as solved.

Domain decomposition The principle of domain decomposition (DD) means dividing the computational domain into subdomains, which are coupled to each other by interface conditions (Smith et al. 1996). Domain decomposition is preferable if subdomains with rather different physical properties can be identified. This is usually the case when dealing with

The finite element method (FEM) Let us assume that the governing equation of an electromagnetic system is of the form E(qE)=sE

(46)

189

in the domain Ωc. When Eq. 37 is considered, E is the Laplacian, sE the total current density, and qE the unknown vector potential. The weighted residual method can be applied to derive a linear system from Eq. 46. To do so, the unknown function qE is approximated by a linear combination of N basis functions φ1,φ2,...,φN with constant coefficients ci: N

qE ≈ qEN = ∑ ciφi .

(47)

i =1

To determine the coefficients ci so that qNE is a reasonable approximation of qE, Eq. 46 is enforced in a general sense, i.e.

ψ j , E(qEN )

Ωc

= ψ j , sE

Ωc

j = 1, 2, K N,

,

(48)

where ψj is a set of weighting functions and the notation <ψ,φ>Ωc indicates an inner product defined by

ψ ,φ

Ωc

= ∫ ψφdΩ .

(49)

Ωc

Assuming that the operator E is linear, we receive from Eqs. 47 and 48: N

∑ ψ j , E(φi )

i =1

Ωi

ci = ψ j , sE

Ωc

,

j = 1, 2, K, N.

(50)

The set of equations given by Eq. 50 can be written in matrix form as [K] {c}={f}

(51)

where [K] is a N×N matrix with elements kji=<ψj,E(φi)>Ωc, {c} is the column vector of the coefficients c1,c2,...,cN and {f} is a column vector with elements fj=<ψj,sE>Ωc. For second-order elliptic problems, the integral expression for the coefficients kji can be integrated by parts by applying Green’s formula. This case is of prime interest because it poses the same continuity requirements on ψ and φ. This means that identical function spaces can be used for both kinds of functions (Strang and Fix 1973). If the weighting functions are chosen to be identical to the basis functions (ψj=φj) the procedure is called Galerkin’s method. The FEM can be regarded as a systematic and flexible means of constructing the basis functions φi. To this end, the problem domain Ωc is subdivided into individual small elements Ωe of simple geometry (e.g., tetrahedrons or hexahedrons). If standard node-based elements are considered, each element type is equipped with a set of node points and corresponding low order polynomials, the shape functions. Within each element Ωe the unknown function qE is represented by a linear combination of the shape functions. The global representation (Eq. 47) can be constructed by assembling all local representations. The basis functions φi obtained that way differ from zero only inside a small number of elements, so the matrix [K] in Eq. 51 contains only a relatively small number of non-zero entries, it is a sparse matrix. [K] is also symmetric if E is a self-adjoint operator, e.g. the Laplacian.

There are various solution techniques which take advantage of these properties. While direct solvers based on elimination are less favorable, because they produce additional matrix fill-in, iterative solvers such as ICCG (= incomplete Cholesky conjugate gradient) are widely used (Jin 1993). Recently, multigrid techniques are ˇ emerging in computational electromagnetism (Cingoski et al. 1999; Hiptmair 1999; Kaltenbacher et al. 2000; Schinnerl et al. 2000). Multigrid solvers behave asymptotically optimally: the number of operations to solve Eq. 51 as well as the memory requirements scale ~N for N→∞, where N is the number of degrees of freedom. Apart from the conventional node-based finite elements where the unknown functions are approximated by interpolating their nodal values, there are additional classes of special finite elements which are well suited for the representation of vector fields. The degrees of freedom of edge elements are associated with the edges of the mesh rather than the nodes. These elements, known also as Whitney elements, guarantee the tangential continuity of a vector field across interfaces without imposing normal continuity. Edge elements have gained a lot of significance for the solution of electromagnetic problems (Barton and Cendes 1987; Kameari 1990; Mur 1994). There are also facet elements with the property that the normal component of a vector field is continuous across facets. From the theoretical point of view, nodal, edge and facet elements are closely related to differential forms (Bossavit 1998). The boundary element method (BEM) We look once more at the governing equation (Eq. 46) of the electromagnetic system and assume sE=0 for simplicity. If the electromagnetic field qE is considered on the boundary Γ of the domain Ωa, two types of boundary data can be identified: the Dirichlet data u=D(qE) and the Neumann data λ=N(qE). When Eq. 39 is considered, E is the Laplacian, sE=0 the impressed exciting current density, qE the vector potential, D the restriction of a function to the boundary and N its normal derivative. By theoretical considerations, a boundary integral equation can be derived which provides a relation between the Dirichlet data u and the Neumann data λ. We obtain for smooth boundaries Γ r r r r r r r u(ξ ) (52) = ∫Γ u ∗ ( x, ξ )λ ( x )dΓ – ∫Γ λ ∗( x, ξ ) u( x )dΓ , 2 where ξ is the observation point, x the integration point and u* (x,ξ), λ* (x,ξ) are kernel functions that result from the application of the operators D and N to the Green’s function G (for details see Jackson 1975). Contributions of impressed currents sE could easily be added to Eq. 52. Taking into account an exciting coil in the BEM domain (as in Fig. 4) would require the computation of its source field by means of Biot-Savart’s law (Jackson 1975). The discretization of the integral equation (Eq. 52) follows the same lines as in the FEM case. The boundary

190

Γ is subdivided into small elements. Due to the reduction of the dimensionality we have 2-D elements (e.g., triangles or quadrilaterals) for 3-D problems. The BEM discretization can easily be obtained by restricting the FEM discretization of Ωc to Γ provided that the BEM–FEM strategy is applied. Basis functions φi are constructed from the elements’ shape functions and the functions λ and u are represented by linear combination: M

M

i =1

i =1

λ = ∑ λiφi , u = ∑ uiφi ,

(53)

where M is the number of basis functions. If standard node-based elements are used, M equals the number of nodal points ξ i on the boundary. The nodal values of λ and u are denoted by λi and ui, respectively. If Eq. 53 is put into Eq. 52 and the integral equation is enforced at the nodal points ξ=ξ j,j=1,2,..., M we end up with a linear system of equations. This procedure is called pointwise collocation. The linear system can be written in matrix form as [G]{λ}=[H]{u}

Fig. 5 3-D finite element mesh of a magnetic valve, which consists of a moving armature, a slotted core and the exciting coil

(54)

where [G] and [H] are M×M matrices and {λ} and {u} are the column vectors of the nodal values λ1,λ2...,λM and u1,u2,...,uM, respectively. The disadvantage of the BEM is that, in contrast to the sparse (and usually symmetric and positive definite) matrices of the FEM, the resulting matrices in the BEM are full and usually nonsymmetric as in the case of pointwise collocation. The solution of the corresponding linear systems is frequently done by using direct solvers. However, direct solvers perform rather poorly in this case, since the memory requirement to store the full matrix scales ~M2 and the operation count even scales ~M3 for M→∞. To overcome these drawbacks, iterative solution algorithms in connection with matrix compression techniques have been proposed (Bebendorf and Rjasanow 2000; Buchau et al. 2000). With these techniques, asymptotically optimal BEM implementations can be achieved. There are many more BEM procedures in literature. It is worth mentioning that the BEM does not necessarily result in unsymmetrical matrices (Costabel 1987; Heise and Kuhn 1996) and that BEM discretizations based on edge and facet elements are gaining importance in low frequency electromagnetic computations (Hiptmair 2000).

Applications Magnetic valve The 3-D finite element mesh of a magnetic valve is displayed in Fig. 5. To minimize eddy current effects, slots have been inserted into the core. Therefore this problem cannot be regarded as axisymmetric. Due to the coupling between finite elements and boundary elements, only material parts have to be meshed with finite elements but

Fig. 6 Current in the coil and armature position versus time for the magnetic valve, both normalized

not the surrounding space. The surface of the structure will be covered automatically by boundary elements. The coupled approach leads to major benefits as explained earlier, namely simplified meshing, inclusion of the asymptotic condition at infinity, no additional difficulties due to the motion of the armature. Two different kinds of results are gained by the numerical computations, namely local and global quantities. Local quantities are the spatial and temporal distributions of different fields, best viewed with contour or vector plots. Among these local quantities are the magnetic field and the eddy current density. Global quantities do not depend on the spatial coordinates. The current in the coil or the magnetic force acting on the armature are typical examples. The evolution of the current in the coil and of the computed position of the armature are displayed in Fig. 6. The coil is fed by an impressed current which runs towards its maximum value exponentially with a time constant of 40 µs. The armature is loaded by hydraulic and spring forces, which can be modeled by a pre-stressed linear spring. After approximately 70 µs the magnetic force exceeds the pre-stress, the armature starts to move and reaches its final position after a response time of 225 µs. Practical experience showed that in spite of the relatively coarse mesh, computed global quantities agree with measured values within 10%, provided that quadratic shape functions have been used and the eddy

191

Fig. 8 Current in the coil and air gap width versus time during the switching process of the relay, both normalized Fig. 7 Computed spatial distribution of the magnetic induction inside the magnetic parts of the relay at t=0.41 ms Table 1 Mesh and computational data for the magnetic valve Number of finite/boundary elements Memory requirements Number of nonlinear iterations per time-step Number of time-steps CPU time (DEC ALPHA workstation)

2722/1218 Approx. 320 MB Average 2.7 44 2.5 h

Table 2 Mesh and computational data for the electromechanical relay Number of finite/boundary elements Memory requirements Number of nonlinear iterations per time-step Number of time-steps CPU time (DEC ALPHA workstation)

1216/1304 Approx. 540 MB Average 1.8 395 79 h

current pattern is not too complex, as in this case. This is sufficiently accurate for many design decisions at the price of a reasonable numerical expense. Some mesh and computational data are collected in Table 1. Electromechanical relay The numerical analysis of an electromechanical relay (Rischmüller et al. 2000) serves as a second example. The magnetic circuit of the relay consists of a magnetic yoke, a moving armature and the exciting coil (see Fig. 7). The switch stack and an additional spring acting on the armature do not influence the magnetic conditions and need therefore not to be discretized. However, in the mechanical part of the model they have been accounted for by an equivalent spring characteristic. Some mesh and computational data can be seen in Table 2. A 12 V DC voltage is applied to the coil at t=0. A magnetic field will build up in the magnetic circuit and therefore a magnetic force will be exerted on the armature. The computed distribution of the magnetic induction at t= 0.41 ms is shown in Fig. 7. The evolution of the current in the coil and of the air gap width have also been calculated and plotted in

Fig. 8. As expected, the current in the coil increases monotonously at first. After about 2.5 ms the magnetic force exceeds the value of the spring force and the relay begins to close. The decrease in the air gap width is accompanied by a strong increase in the inductance of the magnetic circuit. Because of the increasing inductance and also because of the build-up of the magnetic field, an electrical voltage will be induced in the coil, acting against further increase of the current according to Lenz’s law. This reaction is so pronounced that the current may even decrease. At t=4.5 ms the armature reaches its stopping point, the closing process is finished and the current finally starts to rise again. It converges towards its final value, limited by the Ohmic resistance of the coil. This effect shows that it is not possible to study the electromagnetic and mechanical subsystems separately. A correct treatment has to take into account their mutual interaction (see Fig. 2), as demonstrated successfully by this computation.

References Albanese R, Rubinacci G (1990) Formulation of the eddy-current problem. IEE Proc A 137:16–22 Albertz D, Dappen S, Henneberger G (1996) Calculation of the 3D non-linear eddy current field in moving conductors and its application to braking systems. IEEE Trans Magn 32:768–771 Axelsson O (1996) Iterative solution methods. Cambridge University Press, Cambridge Barton ML, Cendes ZJ (1987) New vector finite elements for three-dimensional magnetic field computation. J Appl Phys 61:3919–3921 Bebendorf M, Rjasanow S (2000) Matrix compression for the radiation heat transfer in exhaust pipes. In: Sändig A-M, Schiehlen W, Wendland W (eds) Multifield problems. Springer, Berlin Heidelberg New York, pp l83–191 Biro O, Preis K (1989) On the use of the magnetic vector potential in the finite element analysis of three-dimensional eddy currents. IEEE Trans Magn 25:3145–3159 Bossavit A (1981) On the numerical analysis of eddy-current problems. Comput Methods Appl Mech Eng 27:303–318 Bossavit A (1998) Computational electromagnetism. Academic Press, San Diego, Calif. Brauer J, Ruehl J (1994) 3D coupled electromagnetic and structural finite element analysis of motional eddy current problems. IEEE Trans Magn 30:3288–3291 Brebbia C (1984) The boundary element method for engineers, 2nd edn. Pentech, London

192 Brebbia C, Telles J, Wrobel L (1984) Boundary element techniques. Springer, Berlin Heidelberg New York Buchau A, Huber C, Rieger W, Rucker W (2000) Fast BEM computations with the adaptive multilevel fast mutipole method. IEEE Trans Magn 36:680–684 ˇ Cingoski V, Tsubota K, Yamashita H (1999) Investigation of the efficiency of the multigrid method for finite element electromagnetic field computations using nested meshes. IEEE Trans Magn 35:3751–3753 Costabel M (1987) Symmetric methods for the coupling of finite elements and boundary elements. In: Brebbia C, Wendland W, Kuhn G (eds) Boundary elements IX. Springer, Berlin Heidelberg New York, pp 411–420 Costen R, Adamson D (1965) Three-dimensional derivation of the electrodynamic jump conditions and momentum-energy laws at a moving boundary. Proc IEEE 53:1181–1196 Davies A (1980) The finite element method – a first approach. Oxford University Press, Oxford Dhatt G, Touzot G (1984) The finite element method displayed. Wiley, New York Einstein A (1905) Zur Elektrodynamik bewegter Körper. Ann Phys 17:891–921 Fetzer J (1995) Die Lösung statischer und quasistationärer elektromagnetischer Feldprobleme mit Hilfe der Kopplung der Methode der finiten Elemente und der Randelementmethode. PhD thesis, University of Stuttgart Greenwood D (1988) Principles of dynamics, 2nd edn. PrenticeHall, Englewood Cliffs, N.J. Heise B, Kuhn M (1996) Parallel solvers for linear and nonlinear exterior magnetic field problems based upon coupled FE/BE formulations. Computing 56:237–258 Henrotte F, Nicolet A, Hédia H, Genon A, Legros W (1994) Modelling of electromechanical relays taking into account movement and electric circuits. IEEE Trans Magn 30:3236–3239 Hiptmair R (1999) Multigrid method for Maxwell’s equations. SIAM J Numer Anal 36:204–225 Hiptmair R (2000) Symmetric coupling for eddy current problems. Technical Report 148, SFB 382, University of Tübingen Hsieh K (1995) A Lagrangian formulation for mechanically, thermally coupled electromagnetic diffusive processes with moving conductors. IEEE Trans Magn 31:604–609 Jackson J (1975) Classical electrodynamics, 2nd edn. Wiley, New York Jin J (1993) The finite element method in electromagnetics. Wiley, New York Kaltenbacher M, Reitzinger S, Schöberl J (2000) Algebraic multigrid method for solving 3D nonlinear electrostatic and magnetostatic field problems. IEEE Trans Magn 36:1561–1564 Kameari A (1990) Calculation of transient 3-D eddy currents using edge elements. IEEE Trans Magn 26:466–469

Kriezis E, Tsiboukis T, Panas S, Tegopoulos J (1992) Eddy currents: theory and applications. Proc IEEE 80:1559–1589 Kurz S (1998) Die numerische Behandlung elektromechanischer Systeme mit Hilfe der Kopplung der Methode der finiten Elemente und der Randelementmethode. PhD thesis, University of Stuttgart Kurz S, Fetzer J, Lehner G, Rucker W (1999) Numerical analysis of 3D eddy current problems with moving bodies using BEMFEM coupling. Surv Math Ind 9:13l–l50 Landau L, Lifschitz E (1984) Electrodynamics of continuous media, 2nd edn. (Course in theoretical physics, vol VIII) Pergamon, Oxford Maxwell J (1954) A treatise on electricity and magnetism. Dover Publications, New York Minkowski H (1908) Die Grundgleichungen für die elektromagnetischen Vorgänge in bewegten Körpern. Göttinger Nachrichten, pp 53–116 Mur G (1994) Edge elements, their advantages and their disadvantages. IEEE Trans Magn 30:3552–3557 Muramatsu K, Nakata T, Takahashi N, Fujiwara K (1992) Comparison of coordinate systems for eddy current analysis in moving conductors. IEEE Trans Magn 28:1186–1189 Nicolet A, Delincé F, Genon A, Legros W (1992) Finite elementsboundary elements coupling for the movement modeling in two-dimensional structures. J Phys III France 2:2035–2044 Nysveen A, Nilssen R (1997) Time domain simulation of magnetic systems with a general moving geometry. IEEE Trans Magn 33:1394–1397 Penfield P Jr, Haus H (1967) Electrodynamics of moving media. MIT Press, Cambridge, Mass. Ren Z, Razek A (1994) A strong coupled model for analysing dynamic behaviours of non-linear electromechanical systems. IEEE Trans Magn 30:3252–3255 Rischmüller V, Haas M, Kurz S, Rucker W (2000) 3D transient analysis of electromechanical devices using parallel BEM coupled to FEM. IEEE Trans Magn 36:1360–1363 Schinnerl M, Schöberl J, Kaltenbacher M (2000) Nested multigrid methods for the fast numerical computation of 3D magnetic fields. IEEE Trans Magn 36:1557–1560 Shyamkumar BB, Cendes ZJ (1988) Convergence of iterative methods for nonlinear magnetic field problems. IEEE Trans Magn 24:2585–2587 Smith B, Bjørstad P, Gropp W (1996) Domain decomposition. Cambridge University Press, Cambridge Strang G, Fix G (1973) An analysis of the finite element method. Prentice-Hall, Englewood Cliffs, N.J. Tani K, Yamada T, Kawase Y (1998) A new technique for 3-D dynamic finite element analysis of electromagnetic problems with relative movement. IEEE Trans Magn 34:3371–3374