Coupling of finite and boundary elements for transient eddy current problems

MATHEMATICAL ANALYSIS AND NUMERICAL MATHEMATICS VOLUME: 15 |NUMBER: 2 |2017 |JUNE Coupling of Finite and Boundary Elements for Transient Eddy Current Problems Dalibor LUKAS1,2, Petr KACOR3, Lubomir IVANEK4, Veleslav MACH3, Christian SCHEFFLER5 1Department of Applied Mathematics, Faculty of Electrical Engineering and Computer Science, VSB–Technical University of Ostrava, 17. listopadu 15/2172, 708 33 Ostrava, Czech Republic 2IT4Innovations National Supercomputing Center, VSB–Technical University of Ostrava, 17. listopadu 15/2172, 708 33 Ostrava, Czech Republic 3Department of Electrical Power Engineering, Faculty of Electrical Engineering and Computer Science, VSB–Technical University of Ostrava, 17. listopadu 15/2172, 708 33 Ostrava, Czech Republic 4Department of General Electrical Engineering, Faculty of Electrical Engineering and Computer Science, VSB–Technical University of Ostrava, 17. listopadu 15/2172, 708 33 Ostrava, Czech Republic 5Fraunhofer Institute for Machine Tools and Forming Technology IWU, Reichenhainer Strasse 88, 09126 Chemnitz, Germany dalibor.luk[email protected], petr.k[email protected], lubomir.iv[email protected], veleslav.mac[email protected], christian.sc[email protected] DOI: 10.15598/aeee.v15i2.2262 Abstract. A symmetric coupling of methods of finite and boundary elements for numerical solution of transient eddy current problems is described. This is an essential step in modelling of electromagnetic forming of metalic sheets. The finite element method is employed in the conducting region of the metalic sheet. The boundary element method relies on the StrattonChu representation formula and it models the electromagnetic field in the air including its decay at infinity. We impose external currents by the Biot-Savart law. Keywords Boundary elements, eddy current problem, finite elements. 1. Introduction Eletromagnetic forming of metalic sheets relies on generating pulses of eddy currents, which are imposed by a surrounding coil. This gives rise to the Lorentz forces that are pushing the metalic sheet against a form. In order to analyze and later optimize this metalurgical process we shall model the transient eddy current problem and propose a numerical method that gives accurate enough results. This is the aim of the present paper. Other parts of the model such as contact mechanics, plasticity, and eventually thermal distribution shall be treated elsewhere. We consider a domain Ωint ⊂R3occupied by the metalic sheet and the exterior Ωext := R3\Ωint. The transient eddy current problem reads as follows: for i∈ {int,ext}and (x, t)∈Ωi×R+compute the distributions of the magnetic strength density and the electric intensity H(x, t) := (Hint(x, t)x∈Ωint, Hext(x, t)x∈Ωext, (1) E(x, t) := (Eint(x, t)x∈Ωint, Eext(x, t)x∈Ωext, respectively, that satisfy the low frequency case of Maxwell’s equations ∂ ∂t Hi(x, t) + 1 µ0curl Ei(x, t) = 0, curl Hi(x, t)−σi(x)Ei(x, t) = Ji(x, t), div Hi(x, t) = 0, div Ei(x, t) = 0.(2) Here, µ0>0is the permeability of air, σint >0is the conductivity of the metallic sheet, σext := 0,Jext is the impressed current density, and Jint := 0. The equations are completed by the transmission conditions: for c 2017 ADVANCES IN ELECTRICAL AND ELECTRONIC ENGINEERING 280 MATHEMATICAL ANALYSIS AND NUMERICAL MATHEMATICS VOLUME: 15 |NUMBER: 2 |2017 |JUNE (x, t)∈∂Ωint ×R+ Hext(x, t)−Hint(x, t) = 0, Eext(x, t)−Eint(x, t) = 0,(3) the decay conditions: for |x|→∞and t∈R+ |Eext(x, t)|,|Hext(x, t)| → 0,(4) and the initial conditions: for x∈R3 H(x, 0) = E(x, 0) = 0.(5) There are several approaches to formulate the model in the sense of distributions, which is the best-known concept allowing for geometrical as well as material singularities or jumps. Basically, there are potentialbased formulations [11] and [12], magnetic-field-based (H-based) formulations [6] and [15] and electric-fieldbased (E-based) formulations [1], [2], [3], [4], [5], [10] and [17]. We prefer the latter approach, with which we are experienced [13] and [14]. The rest of the paper is organized as follows: In Sec. 2. we present the E-based variational formulation in Ωint and a finite element discretization. In Sec. 3. we recall Stratton-Chu representation and boundary element method (BEM) in the exterior. Section 4. is devoted to Hiptmair’s symmetric FEMBEM coupling. In Sec. 5. numerical results are presented. We give conclusions in Sec. 6. 2. E-Based FEM After applying curl to the first equation in Eq. (2), ∂/∂t to the second one, and adding both we arrive (up to µ0) at the E-based formulation of Eq. (2): for (x, t)∈Ωi×R+ σi∂ ∂t Ei(x, t) + 1 µ0curl curl Ei(x, t) = −∂ ∂t Ji(x, t), div Ei(x, t)=0.(6) A variational formulation of Eq. (6) was introduced and analyzed in [4]. It reads as follows: find Eint ∈ V := L2(0, T),H(curl; Ωint)∩ H1(0, T),H−1(curl; Ωint)such that µ0σZΩint ∂ ∂tEint(x, t)·v(x)dx | {z } =:hM(∂tEint ),vi +ZΩint curl Eint(x, t)·curl v(x)dx | {z } =:hA(Eint ),vi −ZΓ γNEint(x, t)·γDv(x)dS(x) | {z } =:hf(γNEint ),vi = 0,(7) for all v∈H(curl; Ωint). Here, Γ := ∂Ωint,nis the outer unit normal vector to Ωint, and we define the following Dirichlet and Neumann traces, respectively, γDv(x) := n(x)×(v(x)×n(x)) , γNu(x) := curl u(x)×n(x).(8) The formulation is completed by a boundary condition and the initial condition. We approximate Sobolev space H(curl; Ωint)using the lowest-order Nedelec-I finite elements [16]. We search for a piecewise polynomial approximation Eint(x, t)≈ n X i=1 ei(t)ϕi(x),(9) and arrive at the system of ordinary differential equations (ODEs): for t∈R+ M eint0(t) + A eint(t) + f(γNEint)(t) = 0, eint(0) = 0,(10) where we denote by Mand Athe so-called conductivity matrix and permittivity matrix, respectively. Yet f is to be specified. We can solve the ODEs analytically as far as we are able to find the eigenvalues λand the eigenvectors of the matrix pencil A−λM. This can be typically done for n≤103. Otherwise, we have to employ a time-integration scheme. Note that in a pure FEM the domain Ωint has to be actually extended by a large portion of Ωext so that the support of Jis included. On the boundary of this extended domain the electric field is assumed to vanish, thus, the boundary term fdisappears. On the righthand side there is an extra term related to −∂tJ. 3. E-Based BEM In the exterior domain we follow the approach of Hiptmair [10]. We employ Stratton-Chu representation formula: for (x, t)∈Ωext ×R+ Eext(x, t) = W(γDEext(y, t))(x) −e V(γNEext(y, t))(x) + N(−µ0J0 t(y, t))(x),(11) where e V(λ(y))(x) := ZΓ λ(y)1 4π|x−y|dS(y),(12) is the vectorial single-layer operator, W(u(y))(x) := e V(n(y)×u(y))(x),(13) is the Maxwell double-layer operator, and N(g(y))(x) := ZΩext g(y)1 4π|x−y|dS(y),(14) c 2017 ADVANCES IN ELECTRICAL AND ELECTRONIC ENGINEERING 281 MATHEMATICAL ANALYSIS AND NUMERICAL MATHEMATICS VOLUME: 15 |NUMBER: 2 |2017 |JUNE is the Newton potential (Biot-Savart law). Note that in general there are two additional terms in the formula, which vanish in our formulation. Applying γDto the Stratton-Chu formula leads to the first-kind boundary integral equation: for (x, t)∈ Γ×R+ γDEext −γDW(γDEext) | {z } =:−B(γDEext)(x,t) +γDe V(γNEext) | {z } =:V(γNEext)(x,t) =γDN(−µ0J0 t) | {z } =:c(x,t) .(15) Applying γNto Stratton-Chu formula gives rise to the second-kind boundary integral equation: for (x, t)∈ Γ×R+ γNEext =γNW(γDEext)(x, t) | {z } =:−D(γDEext)(x,t) −γNe V(γNEext) | {z } =BT(γNEext)(x,t) +γNN(−µ0J0 t) | {z } =:b(t) .(16) The boundary integral equations Eq. (15) and Eq. (16) are again understood in the Sobolev variational framework [7] and [8], which allows a stable boundary element discretization. The discrete space consists of tangential traces of Nedelec-I elements [16], the so-called stream functions. 4. FEM-BEM Coupling We need FEM to properly model the transient behaviour of eddy currents. However, FEM can approximate decay condition Eq. (4) only at the high cost of additional volume discretization of a large portion of Ωext. On the other hand, BEM models the decay condition by definition, but it suffers from modelling of transient fields in Ωint. Fortunately, there is a natural coupling of FEM and BEM, cf. [9]. It allows us to get rid of the unknown Neumann boundary data γNEint in Eq. (7). From Eq. (15) we can eliminate the Neumann data of the exterior field γNEext =V−1c+B(γDEext).(17) Plugging the latter to Eq. (16) we arrive at a boundary integral equation with the exterior Steklov-Poincare (Dirichlet-to-Neumann) operator S γNEext = −D+BTV−1B | {z } =:S γDEext −b−BTV−1c | {z } =:d .(18) Now the latter and transmission condition Eq. (3) replaces the boundary term in Eq. (7) hM(∂tEint),vi+hA(Eint),vi+hS(γDEint),γDviΓ =−hd,γDviΓ.(19) Finally, we employ the finite element discretization Eq. (9) and arrive at ODEs: for t∈R+ M eint0(t)) + A+ (IΓ)TS IΓ | {z } =:K eint(t) = −d(t), eint(0) = 0,(20) where IΓis the restriction to the boundary degrees of freedom (identity matrix completed by zeros). The resulting FEM-BEM system of ODEs can be analytically integrated in time eint(t) = − n X i=1 Zt 0 vi·d(τ) e−λi(t−τ)dτvi,(21) where K vi=λiM vi,kvikM= 1.(22) 5. Numerical Results We consider an axisymmetric setup of a coil and an aluminium plate disc, see Fig. 1. The radius of the disc is 8 cm and the disc is 2 mm thin. It is placed 2 mm above the coil. The coil is modeled by 3 line circular turns of radii rk∈ {2.1,3.7,5.3}cm. Hence, we replace the Newton potential by the following dimensionallyreduced Biot-Savart law: N(x) := −µ0I 4π 3 X k=1 2π Z0 (−sin t, cos t, 0) kx−rk(cos t, sin t, 0)krkdt. (23) The excited current pulse has the amplitude I:= 100 kA. The shape g(t)is half of the sine at frequency f:= 8.33 kHz, g(t) := (sin(2π f t), t ≤1 2f, 0, t ≥1 2f.(24) We employ the particular solution approach: find Eint(x, t) = Eint 0(x, t) + g0(t)N(x)so that Eint 0∈ V, Eint 0(x, 0) = 0, and h(∂t+K)Eint 0,vi=−g00(t)hN,vi ∀v∈H(curl; Ωint), which is the counterpart to Eq. (19). In Fig. 2 we depict a comparison of eddy current distribution and Lorentz forces computed by FEM and FEM-BEM methods. The numbers of uknowns were 9751 in case of FEM and 800 in case of FEM-BEM. The difference in the Lorentz force magnitudes is shown in Fig. 3. The FEM-BEM results are, by definition, more precise. c 2017 ADVANCES IN ELECTRICAL AND ELECTRONIC ENGINEERING 282 MATHEMATICAL ANALYSIS AND NUMERICAL MATHEMATICS VOLUME: 15 |NUMBER: 2 |2017 |JUNE 01 01 0011 0011 01 01 000000000000000000000000000000000000000 000000000000000000000000000000000000000 000000000000000000000000000000000000000 000000000000000000000000000000000000000 000000000000000000000000000000000000000 000000000000000000000000000000000000000 000000000000000000000000000000000000000 111111111111111111111111111111111111111 111111111111111111111111111111111111111 111111111111111111111111111111111111111 111111111111111111111111111111111111111 111111111111111111111111111111111111111 111111111111111111111111111111111111111 111111111111111111111111111111111111111 Fig. 1: Geometry of the example. On the left figure a sketch of the device is depicted. The metalic plate (solid rectangle) is pushed against the form (hatched object on the top). Outward and inward orientation of the currents in the three circular turns is depicted with circles and crosses, respectively. The right figure shows the situation in 3D. Fig. 2: Comparison of eddy current distributions and Lorentz forces calculated by FEM (top figure) and FEM-BEM (bottom figure) at the half-period of the pulse. c 2017 ADVANCES IN ELECTRICAL AND ELECTRONIC ENGINEERING 283 MATHEMATICAL ANALYSIS AND NUMERICAL MATHEMATICS VOLUME: 15 |NUMBER: 2 |2017 |JUNE Fig. 3: Evolution of the maximal Lorentz force density computed by FEM (blue line) and FEM-BEM (red line) methods. 6. Conclusion We presented a coupling of FEM and BEM for solution of transient eddy current problems that arise in the process of electromagnetic forming of metalic sheets. In our forthcoming work we shall complete the model with contact mechanics, plasticity, and thermal field distribution. Acknowledgment The work is supported by the European Union, the Ministry of Industry and Trade, Czech Republic, under the OPEIC project No. CZ.01.1.02/0.0/0.0/15_007/0002298, and by the Federal Ministry of Economics and Energy, Germany, under the IGF project 173 EBR. References [1] ALONSO, A. A mathematical justification of the low-frequency heterogeneous timeharmonic Maxwell equations. Mathematical Models and Methods in Applied Sciences. 1999, vol. 9, iss. 3, pp. 475–489. ISSN 1793-6314. DOI: 10.1142/S0218202599000245. [2] ACEVEDO, R., S. MEDDAHI and R. RODRIGUEZ. An E-based mixed formulation for a time-dependent eddy current problem. Mathematics of Computation. 2009, vol. 78, iss. 268, pp. 1929–1949. ISSN 1088-6842. [3] ACEVEDO, R. and S. MEDDAHI. An Ebased mixed FEM and BEM coupling for a time-dependent eddy current problem. IMA Journal of Numerical Analysis. 2011, vol. 31, iss. 2, pp. 667–697. ISSN 0272-4979. DOI: 10.1093/imanum/drp049. [4] BUFFA, A., H. AMMARI and J. C. NEDELEC. A justification of eddy currents model for the Maxwell equations SIAM Journal on Applied Mathematics. 2006, vol. 60, iss. 5, pp. 1805–1823. ISSN 0036-1399. DOI: 10.1137/S0036139998348979. [5] ARNOLD, L. and B. HARRACH. A unified variational formulation for the parabolic-elliptic eddy current equations. SIAM Journal on Applied Mathematics. 2012, vol. 72, iss. 2, pp. 558–576. ISSN 0036-1399. DOI: 10.1137/110831477. [6] BERMUDEZ, A., D. GOMEZ, R. RODRIGUEZ and P. VENEGAS. Numerical analysis of a transient non-linear axisymmetric eddy current model Computers and Mathematics with Applications. 2015, vol. 70, iss. 8, pp. 1984–2005. ISSN 08981221. DOI: 10.1016/j.camwa.2015.08.017. [7] BUFFA, A. and P. CIARLET. On traces for functional spaces related to Maxwell’s equations. Part I: An integration by parts formula in Lipschitz polyhedra. Mathematical Methods in the Applied Sciences. 2000, vol. 24, iss. 1, pp. 9–30. ISSN 1099-1476. DOI: 10.1002/1099-1476(20010110)24:1<9::AIDMMA191>3.0.CO;2-2. [8] BUFFA, A. and P. CIARLET. On traces for functional spaces related to Maxwell’s equations Part II: Hodge decompositions on the boundary of Lipschitz polyhedra and applications. Mathematical Methods in the Applied Sciences. 2000, vol. 24, iss. 1, pp. 31–48. ISSN 1099-1476. DOI: 10.1002/1099-1476(20010110)24:1<31::AIDMMA193>3.0.CO;2-X. [9] CONSTABEL, L. Symmetric Methods for the Coupling of Finite Elements and Boundary Elements (Invited contribution). In: Boundary Elements IX: Mathematical and Computational Aspects. Berlin: Springer, 1987, pp. 411–420. ISBN 978-3-662-21908-9. DOI: 10.1007/978-3-66221908-9_26. [10] HIPTMAIR, R. Symmetric Coupling for Eddy Current Problems. SIAM Journal on Numerical Analysis. 2002, vol. 40, iss. 1, pp. 41–65. ISSN 1095-7170. DOI: 10.1137/S0036142900380467. [11] KUHN, M. and O. STEINBACH. Symmetric coupling of finite and boundary elements for exterior magnetic field problems. Mathematical Methods in the Applied Sciences. 2002, vol. 25, iss. 5, pp. 357– 371. ISSN 1099-1476. DOI: 10.1002/mma.286. c 2017 ADVANCES IN ELECTRICAL AND ELECTRONIC ENGINEERING 284 MATHEMATICAL ANALYSIS AND NUMERICAL MATHEMATICS VOLUME: 15 |NUMBER: 2 |2017 |JUNE [12] L’EPLATTENIER, P., G. COOK, C. ASHCRAFT, M. BURGER, J. IMBERT and M. WORSWICK. Introduction of an Electromagnetism Module in LS-DYNA for Coupled Mechanical-Thermal-Electromagnetic Simulations. Steel Research. 2009, vol. 80, iss. 5, pp. 351– 358. ISSN 1869-344X. DOI: 10.2374/SRI08SP152. [13] LUKAS, D., K. POSTAVA, O. ZIVOTSKY and J. PISTORA. Optimization of electromagnet for high-field polar magneto-optical microscopy. Journal of Magnetism and Magnetic Materials. 2010, vol. 41, iss. 32, pp. 1471–1474. ISSN 0304-8853. DOI: 10.1016/j.jmmm.2009.07.040. [14] LUKAS, D., K. POSTAVA and O. ZIVOTSKY. A shape optimization method for nonlinear axisymmetric magnetostatics using a coupling of finite and boundary elements. Mathematics and Computers in Simulation. 2012, vol. 82, iss. 10, pp. 1721–1731. ISSN 0378-4754. DOI: 10.1016/j.matcom.2011.01.015. [15] MEDDAHI, S. and V. SELGAS. An H-based FEM-BEM formulation for a time dependent eddy current problem. Applied Numerical Mathematics. 2008, vol. 58, iss. 8, pp. 1061–1083. ISSN 01689274. DOI: 10.1016/j.apnum.2007.04.002. [16] NEDELEC, J. C. Mixed finite elements in R3.Numerische Mathematik. 1980, vol. 35, iss. 3, pp. 315–341. ISSN 0945-3245. DOI: 10.1007/BF01396415. [17] REN, Z. and A. RAZEK. New technique for solving three-dimensional multiply connected eddy-current problems. IEE Proceedings A - Physical Science, Measurement and Instrumentation, Management and Education. 1990, vol. 137, iss. 3, pp. 135–140. ISSN 0143702X. DOI: 10.1049/ip-a-2.1990.0021. About Authors Dalibor LUKAS was born in Ostrava, Czech Republic. He received his Ph.D. from VSB–Technical University of Ostrava in 2003. His research interests include numerical analysis and optimization for partial differential equations. Petr KACOR was born in Ostrava, Czech Republic. He received his Ph.D. from VSB–Technical University of Ostrava in 2003. His research interests include electromagnetic and thermal field analysis of electrical machines and devices. Lubomir IVANEK was born in Frydek Mistek, Czech Republic. He received his Ph.D. from Czech Technical University in Prague in 1993. His research interests include electromagnetism, electromagnetic wave propagation and antennas. Veleslav MACH was born in Ostrava, Czech Republic. He received his Ph.D. from VSB–Technical University of Ostrava in 1997. His research interests include electromagnetism and high-voltage techniques. Christian SCHEFFLER was born in KarlMarx-Stadt, Germany. He received his M.Sc. from Technische Universitaet Chemnitz in 2009. His research interests include electromagnetic forming. c 2017 ADVANCES IN ELECTRICAL AND ELECTRONIC ENGINEERING 285