Mathematical Analysis of Modified BEM-FEM Coupling Approach for 3D Electromagnetic Levitation Problem

In electromagnetic analysis, the finite element and boundary element methods jointly known as 'FEM-BEM coupling' is applied for numerically solving levitation problem based on eddy current. The main focus behind this coupled analysis method is to determine the dynamic characteristic of the levitating body in the presence of a magnetic field. An innovative 3D structure is developed that couples Lagrangian description and BEM-FEM coupling method for this purpose. The coupling methodology is based on the boundary conditions on the common boundaries between FEM and BEM sub-domains. Subsequent coding has been developed to simulate the problem in the MATLAB environment. An example similar to TEAM (Testing Electromagnetic Analysis Methods) workshop problem 28 has been used to study the efficiency of code for computationally inexpensive analysis.



There are no comments yet.


page 1

page 2

page 3

page 4


Caccioppoli-type estimates and ℋ-Matrix approximations to inverses for FEM-BEM couplings

We consider three different methods for the coupling of the finite eleme...

Hybrid coupling of finite element and boundary element methods using Nitsche's method and the Calderon projection

In this paper we discuss a hybridised method for FEM-BEM coupling. The c...

An optimization-based strategy for peridynamic-FEM coupling and for the prescription of nonlocal boundary conditions

We develop and analyze an optimization-based method for the coupling of ...

A multi-grid sampling multi-scale method for crack initiation and propagation

In this study, a multi-grid sampling multi-scale (MGSMS) method is propo...

Deep Convolutional Neural Networks to Predict Mutual Coupling Effects in Metasurfaces

Metasurfaces have provided a novel and promising platform for the realiz...

On the Lagrangian-Eulerian Coupling in the Immersed Finite Element/Difference Method

The immersed boundary (IB) method is a non-body conforming approach to f...
This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

I Introduction

The modelling of dynamics of a conducting object during external magnetic field’s presence can be calculated by solving the equation of motion considering the magnetic forces for induced eddy currents. Traditional difficulties in this area are the analysis of the motion of conducting body, re-meshing techniques, force calculation, weak versus strong electromechanical coupling, efficient time-stepping schemes etc. [1],[2]. TEAM (Testing Electromagnetics Analysis Method) workshop problems 9 and 17 manage moving bodies with external magnetic field[5]. In contrast, TEAM workshop problem 28 is a transient problem with electromechanical coupling. The proposed model for evaluating levitating setup is similar to TEAM workshop problem 28.

In general, the discretization of a domain by finite elements which contain movable parts is not straightforward to handle[2]. Even a small displacement can create huge mesh distortion; subsequent re-meshing of the entire domain is necessary for tackling the continuous displacement of the moving parts. This shortcoming could be removed [3] by the BEM-FEM coupling approach. The advantage of this approach is meshing pattern does not alter during the motion. It gets shifted according to the levitating object due to BEM application and there is no need for entire domain re-meshing. The behavior of the levitated body is described by Maxwell’s equations and other constitutive equations. Solving the Electromechanical problem by those equations in three dimensional space is a time-consuming CPU task, consequently barely resorted [4]. The ongoing increment of computation capability and constantly changing numerical techniques have dramatically diminished the need for time-consuming calculations. During evaluation, ’BEM-FEM’ coupling code has been developed in MATLAB environment. Laplace equation is used as governing equation in the air (BEM) domain for simulation. Results have been compared with measured values from Multiphysics software ”COMSOL.”An analysis without detailed mathematical treatment of the BEM approach is provided in [5]

to calculate the magnetic vector potential at the boundary of the conducting medium. In this paper, an extension of previous work

[5], based on mathematical ground is provided.

Ii Electromagnetic Formulation

Maxwell’s equations that describes all classical electromagnetic phenomena, are given by (1)-(4).


Displacement current is omitted for sake of simplifications.Above (3) and (4) displayed an crucial role for development of transient electromagnetic problem.Useful identity of vector calculus gives (5)


A’ and V’ are magnetic vector potential and electric scalar potential respectively. Object under study has no residual magnetism. Accounting constitutive relations in (6) and (7).


Here, , , and are impressed current density, magnetic field strength, conductivity and velocity respectively. Substituting (5) in (7) gives,


Simplification of (8) yields


The divergence of current density is also zero.


Substituting (10) into (8) gives (11).

Fig. 1: structure of the considered domain[5]

The considered domain is decomposed as shown in Fig. 1 into different conducting region like , a region free of eddy current (air domain). In multiply connected region , there is no conducting body but it may contain a current source, The BEM treats the surrounding air space. The following continuity conditions are valid on coupling boundaries for magnetic vector potential, magnetic field intensity. Current density is continuous across the coupling boundary.


Iii Electromechanical System Modeling

Fig. 2: Node moving from x to x+ within t[5]

In the general case of movable boundaries, at least a portion of the meshing structure must move as a whole [7] with the boundary to account for the interface conditions. Fig. 2 shows a node which is part of a moving mesh. The temporal variation of the vector potential at this node is mentioned by [3].


Assuming nodal elements, the temporal variation of the magnetic vector potential values gives an approximation for the total time derivative rather than the partial time derivative [3]. Hence


substituting (17) in (9), we have


similarly substituting (17) in (11)


Where =

From (18)-(19), it gives that moving the mesh [8] and replacing V by automatically takes into account the subsequent motional E.M.F till [3]. This is a scenario for 2-D problems and 3-D problems with small perturbation. In those cases, it is completely sufficient to move the mesh connected to the respective moving bodies for all electromagnetic effects due to motion.

Iv Fem-Bem Coupling Method

The application of the Galerkin method for the weak integral form of (18)-(19) and using nodal elements and boundary conditions yields the coupled matrix of FEM-BEM method.FEM-BEM method across coupling boundaries is given by (20).


In order to apply Galerkin in (18) and (19), we multiply (18) with a suitable weight function.


where, a = 1, 2, 3.



Application of boundary condition,vector identity, weight function and further simplification yields (23).


Simplification gives (24).


It can be further simplified again using vector analysis, and taken into account and and it can be written as follows-


The fundamental solution of scalar Laplace equation in 3-D is given by following equation -


It gives therefore


By substituting (27) in (25) gives (28).


Where, ; with being called as edge factor and is the solid angle at the point rim area inside. Quantity c(r) is equal to 1 if the field point belongs to the interior point. If it lies on the smooth part of the boundary, then it is equal to , and it takes 0 value if the point lies outside. Till now, the boundary element method-based integral equation has been derived whose solution can be obtained by clustering of the boundary and selection of appropriate basis function A and [9]. The edge is in several non-core elements decomposed.The basic condition for coupling along boundary [3]- [7].


The vector element A and interpolated to (boundary element) to evaluate the integrals (30) numerically [5]. The functions used here are the same as in the finite element method. The standard Lagrangian method and subsequent transformation reduced the dimension of the problem by one unit. The matrix format of (30) follows (31).


Similarly the FEM domain written in matrix form gives (32).


Where in domain


[1ex]   =

[1ex]   =

Now combination of (31) and (32) with above condition yields the local matrix of FEM-BEM coupling.


Where k, c, t are respectively element stiffness matrix, damping matrix and boundary matrix.The element equation for are obtained for obtaining the global matrix for FEM-BEM coupling [9]


Now, the quantity under first and second bracket in (34) is denoted by and respectively where refers to matrix equivalent of finite element and (34) is a linear differential equation of first order[9] .

V Modified Approach of Fem-Bem Coupling

The traditional approach of BEM-FEM coupling in 3D proceeds as stated in the previous section. The simultaneous implementation of the FEM and BEM approach raises the computational load and execution time. According to the modified approach proposed here, FEM has to be implemented once during the whole process in the area encircled by the internal circle as shown in Fig. 5-6. At external circle in Fig. 5-6, boundary condition are nearly zero as stated earlier. FEM provides as the input of BEM and due to lack of magnetization vector in Aluminium body. The output of FEM remains almost constant during the levitation process. FEM has been adopted once in the internal circle domain, and output is obtained on the boundary of the internal circle as in Fig.5-6. BEM is used as the propagation tool to obtain the required result at the boundary of the metallic object. In order to make the calculations simpler, the entire 3D object has been projected to 2D, and analysis, computations have been done in 2D. It is analogous to the realization that the planar view of the levitation process is observed and the depth of the body has been taken as unity. In planar view, the governing equation for FEM and BEM sub-domains is converted to 2D Laplace and 2D poison’s equation, respectively.

In mathematics for differential equations, a boundary value problem is a differential equation together with a set of additional constraints, called the boundary conditions. A solution to a boundary value problem is a solution to the differential equation, satisfying the boundary conditions. Dirichlet (or first-type) boundary condition is a type of boundary condition when imposed on an ordinary or a partial differential equation, it gives those values which need to take on along the boundary of the domain. Neumann (or second-type) boundary condition is a type of boundary condition imposed on an ordinary or partial differential equation. It specifies derivative based values of a solution on the domain’s boundary. Robin boundary conditions are a weighted combination of Dirichlet boundary conditions and Neumann boundary conditions (mixed boundary condition). In contrast to mixed boundary conditions, only one type of boundary conditions are applied on any portion of the boundary traditionally: either the function value is specified or the normal derivative but not both.


For some non-zero constants, x and y and a predefined function defined on . Here, is the unknown solution defined on , and denotes the normal derivative at the boundary. More generally, x and y are functions rather than constants.In the BEM sub-domain the governing equation is given by [10].


Fundamental solution of ’’ plays an pivotal role in formation of boundary integral equation and it is given by (37).Where (x, y) is the source point and is the field point.


From Gauss Divergence theorem and (36) it can be stated that


and are two arbitrary solutions of Laplace’s equation in region R bounded by simple closed curved C. and are normal derivatives of the functions. Fundamental solution and the general solution satisfy (36).Hence,


In (39), the fundamental solution possesses a discontinuity at , . Hence,

Fig. 3: Modified approach when a field lie in the interior[5]

If R then C has to be replaced by C [10]. is a circle of center radius has center as shown in Fig. 3. This is because along with its partial derivatives are well defined in the annular region between and .

Fig. 4: Modified approach when a field lie on the boundary [5]

If C lies on the boundary, then C has to be changed by as in Fig. 4. If has the center and radius , then D is the part of C that lies outside . is that of that is inside R. Actually to solve the equations numerically, we need to discretize the boundary integral equation. That is


Where is any arbitrary point in and

In this approximation procedure of and on the boundary condition of C is required. Here we have taken constant approximation, i.e., it is assumed that and functions are constants over where . More precisely it can be said that and where and are values of and at midpoints of . Here, is one arc of the boundary or simply a boundary element. Now, the boundary integral equation boils down to (41).


Here, the interior region ’R’ is the annulus created by two circles. is the boundary of a doubly connected annular domain. Hence, discretizing into 80 boundary segments where precisely 40 of those lie on the boundary of the smaller circle and the other 40 elements lie on the other one. Constant shape functions have been taken to simplify the formulation process in (41) and it leads to (42).





As it can be seen from (42) both and are needed on every individual boundary element . Here is the Dirichlet boundary condition and is the Neumann boundary condition. Generally, both are not given for any node; hence, to apply (42), we have to find missing Dirichlet conditions where only Neumann conditions are given and vice versa. Basically, they are complementary boundary conditions. In order to get , complementary boundary conditions have to be obtained by projecting (42) on to the boundary . Hence from definition of , we get (43)


In a matrix form (43) can be written as (44).


In (44), lies on the smooth part of the boundary should be taken as . As it has been mentioned above, both Dirichlet and Neumann conditions are needed at every node. But generally, one of the conditions is given on each node. In (42), all and can be calculated by known formula. Hence it can be solved for N unknowns by N equations. As some node (that is for some k) is unknown, and at some other node is unknown. Now coming to this case, it is easy to obtain values at every node. So, the unknown quantities are values (the Neumann conditions).From (44), it can be observed that, 80 equations and 80 unknown values. The equation in (44) can be re-written as (45).


For M=1, 2, 3,……, N

Once the equation in (45) is solved for the unknown , , ,……, , the values of and over the element are given by and respectively, where K=1, 2, 3,……, N. Then, =1 is used to compute the value of the function at any interior point. Subsequently, and values need to be computed. The points on the element are defined by (46)


Where is the length of the and [, ] = is the normal vector to pointing away from R. For (x,y) , it can be found that

= =

Using (46), it can be stated that



Distance ’S’ can be expressed as





As a whole from mathematical point of view


and can be written as


It can be mathematically proved that if lies on the boundary element on which we are integrating, then it creates a singularity in both the integrals. So two different cases have to be considered for each integral.
- case 1:
In this case, doesn’t lie on the boundary element on which integration is performed. Analytically this integration can be stated as in (51)


Case 2: In this case lies on the boundary element but doesn’t lie on the boundary on which we are integrating. Hence,


- case 1:
In this case doesn’t lie on the boundary element on which we are performing integration. So integrating analytically, it can be stated that


Case 2:
In this case lies on the boundary element but doesn’t lie on the boundary on which we are integrating. Hence


Based on multiple parameters of the system, which has been detailed in section II, a one-time FEM on the entire domain via ”COMSOL MULTIPHYSICS” software has been performed. Using the FEM, magnetic vector potential values on both circles, i.e. at , k=1, 2, ……, 80 have been obtained. The output can be represented as

The values obtained through COMSOL is exported to the MATLAB environment. Using the FEM-BEM coupling condition, the output of FEM has been used as the input of BEM.

Then, substituting (44) in (43) in for k=1, 2, 3, ………, 80 have been obtained. After getting this complementary solution set, at any point in the interior part of the annulus can be calculated by (56)


In a more abridged form, (55) can be written as


Where R. Hence, in further calculating the magnetic vector potential at any interior point of the domain, the computation process will solely depend upon BEM and will be free from the adaptive domain discretization process. In the modified approach, subsequent coding has been developed in a MATLAB environment to calculate magnetic vector potential at different points of the aluminum plate.

Vi Levitating System Description

The simulation setup of the problem follows Fig. 5. A cuboidal plate of aluminum of conductivity (=3.7m), having a depth of 1 mm, is placed above one cylindrical coil (incoming and outgoing) having 960 turns. Coils are shown by domains that carry a homogeneous azimuthal current density for numerical purposes. The coils are made of copper wire with 1.2 mm diameter and contain insulating layers. The levitation height ’Y’ denotes to the distance between the midpoint of the lower edge of the plate and the upper edge of the current-carrying area (Y = 0) for t 0 the plate lies above the coils at a distance of 18 mm due to the thickness of the winding form. Centre of gravity of the system at (x = 0 mm, y = 18 mm). For calculation purposes, two fictitious circles are defined; the inner one has a radius of 15 mm, and the outer ring has a radius of 100 mm. After some instant due to oscillations for the applied force, the plate shifts 2 mm upward and 2 mm left-wise and rotated at an angle of 10 (refer to final position) as shown in Fig. 6 (not to scale).For t 0, sinusoidal current i(t) flows in the coil segment in opposite directions.

Fig. 5: Magnetic levitation set up at t = 0[5]
Fig. 6: Disturbed levitation set up after some instant[5]

The interaction of the induced eddy current on current transferring conductor is omitted for easy calculation purpose. For this reason, the current in (57) can be regarded as impressed.

Vii Results and Discussions

Fig. 7: Magnetic vector potential distribution at initial condition [5]
Fig. 8: Magnetic vector potential distribution at some instant after applying disturbance [5]
Measured Calculated Error Average
value (T m) value (T m) error
-0.0065 -0.0066 1.53% 7.47%
-0.0025 -0.0028 12%
0.0024 0.0025 4.16%
0.0063 0.0062 1.58%
-0.0090 -0.0085 5.55%
-0.0046 -0.0040 13.04%
0.0042 0.0037 11.90%
0.0090 0.0081 10%
Measured Calculated Error Average
value (T m) value (T m) error
-0.0045 -0.0045 0% 5.87%
-0.0013 -0.0013 0%
0.0026 0.0030 15.38%
0.0067 0.0069 2.98%
-0.0067 -0.0063 5.97%
-0.0031 -0.0025 19.35%
0.0044 0.0044 0%
0.0092 0.0089 3.26%
Fig. 9: Local error analysis at initial time[5]
Fig. 10: Local error analysis after displacement [5]
Fig. 11: Relative analysis of average error with number of nodes [5]

A comparative analysis of the magnetic vector potential obtained from Multi-physics software COMSOL using re-meshing strategies as per Fig.7-8 and computed results from MATLAB are shown in Table I-II obtained with the help of developed BEM-FEM code. Relative error analysis graphs as in Fig. 9-10 are also attached to know the behavior of error with a variation of nodal points. With an increasing number of node points on the boundary of the impressed current domain error, it can be seen that error is decreased at a sharp rate. This work is evaluated in Dell Inspiron 15 series computer operating with 2.1 GHz processor with 4 GB RAM. The interface between the aluminum plate and air coincides with coupling boundary =

and it has 8 nodes. The evaluated magnetic vector potential on the boundary of the aluminum body serves as input for tangential and normal force density and subsequently force and torque calculation for Maxwell stress tensor-based method during comparison of hardware and software-based result. The drawback of the work is that we assumed a ’constant approximation’ in gradient calculation. If we assume ’linear’ or ’quadratic approximation’, the error between the measured values and calculated values will be very less. A fictitious circle is drawn around a metal plate and current-carrying conductor domain. However, it may not be possible to draw a circle in between them due to constrain in intermediate distance.

Viii Conclusion

Magnetic vector potential on the boundary of the conducting surface has been calculated using the modified FEM-BEM coupling approach. A good agreement between the calculated results from two different software has been observed. The efficiency of the method has been illustrated by comparing the results obtained from nodal performance analysis. Previously, coding of BEM-FEM coupling approach has been done in FORTRAN environment, MATLAB based coding is newly introduced in this paper. Significant savings in computation time and storage requirements are achieved. Magnetic forces, normal force density, and tangential force density can be calculated from magnetic vector potential as per the standard formula from the value obtained on the boundary. This concept is extended to compute net force and torque by simply integrating force and torque over the boundary, assuming that the boundary is a coupled one. Then only the analysis of the dynamic behavior of the levitating body is possible.


  • [1] H. Karl, J. Fetzer, S. Kurz, G. Lehner and W. M. Rucker, “Description of TEAM Workshop Problem 28: An Electrodynamic Levitation device”.
  • [2] S. Kurz, J. Fetzer and G. Lehner, “Three-dimensional transient BEM-FEM coupled analysis of electrodynamic levitation problems,” IEEE Trans. Magn., vol. 32, no. 3, pp. 1062-1065, 1996.
  • [3] J. Fetter et al, “Transient BEM-FEM coupled analysis of 3-D electromechanical systems: a watch stepping motor driven by a thin wire coil,” IEEE Trans. Magn., vol. 34, no. 5, pp. 3154-3157, 1998.
  • [4] S. Suuriniemi, K. Forsman, L. Kettunen and J. Makinen, ”Computation of eddy currents coupled with motion,” IEEE Trans. Magn., vol. 36, no. 4, pp. 1341-1345, July 2000.
  • [5] A. Jena, S. Sarkar and Mashuq-un-Nabi, ”A modified BEM-FEM coupling approach for 3D electromagnetic levitation problem,” 2015 Annual IEEE India Conference (INDICON), New Delhi, 2015, pp. 1-6.
  • [6] A. Nicolet, F. Delincé, A. Genon and W. Legros, “Finite elements- boundary elements coupling for the movement modelling in two- dimensional structures,” J. Phys. III, vol. 2, pp. 2035-2044, 1992.
  • [7] S. Kurz, J. Fetzer, G. Lehenr, W. M. Rucker, “A novel formulation for 3D eddy current problems with moving bodies using a Lagrangian description and BEM-FEM coupling,” IEEE Trans. Magn., vol. 34, no. 5, pp. 3068-3073, 1998.
  • [8] J. Fetzer, S. Kurz, G. Lehner and W. M. Rucker, “Analysis of an actuator with eddy currents and iron saturation: comparison between a FEM and a BEM-FEM coupling approach,” IEEE Trans. Magn., vol. 35, no. 3, pp. 1793-1796, 1999.
  • [9] S. Kurz, J. Fetzer, G. Lehner and W. M. Rucker, “The application of the BEM-FEM coupling method for the treatment of three-dimensional nonlinear shielding problems of low-frequency fields using the example of the TEAM problem 21,” Archive for Electrical Engineering, vol. 80, no. 2, pp. 91-104, 1997.
  • [10] W. T. Ang, A Beginner’s Course in Boundary Element Method, Universal Publisher, 2007.