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, remeshing techniques, force calculation, weak versus strong electromechanical coupling, efficient timestepping 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 remeshing of the entire domain is necessary for tackling the continuous displacement of the moving parts. This shortcoming could be removed [3] by the BEMFEM 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 remeshing. 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 timeconsuming CPU task, consequently barely resorted [4]. The ongoing increment of computation capability and constantly changing numerical techniques have dramatically diminished the need for timeconsuming calculations. During evaluation, ’BEMFEM’ 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).
(1) 
(2) 
(3) 
(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)
(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).
(6) 
(7) 
Here, , , and are impressed current density, magnetic field strength, conductivity and velocity respectively. Substituting (5) in (7) gives,
(8) 
Simplification of (8) yields
(9) 
The divergence of current density is also zero.
(10) 
Substituting (10) into (8) gives (11).
(11) 
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.
(12) 
(13) 
(14) 
(15) 
Iii Electromechanical System Modeling
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].
(16) 
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
(17) 
substituting (17) in (9), we have
(18) 
similarly substituting (17) in (11)
(19) 
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 2D problems and 3D 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 FemBem 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 FEMBEM method.FEMBEM method across coupling boundaries is given by (20).
(20) 
(21) 
where, a = 1, 2, 3.
(22) 
Where,
Application of boundary condition,vector identity, weight function and further simplification yields (23).
(23) 
Simplification gives (24).
(24) 
It can be further simplified again using vector analysis, and taken into account and and it can be written as follows
(25) 
The fundamental solution of scalar Laplace equation in 3D is given by following equation 
(26) 
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 methodbased 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 noncore elements decomposed.The basic condition for coupling along boundary [3] [7].
(29) 
(30) 
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).
(31) 
Similarly the FEM domain written in matrix form gives (32).
(32) 
Where in domain
=
[1ex] =
[1ex] =
Now combination of (31) and (32) with above condition yields the local matrix of FEMBEM coupling.
(33) 
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 FEMBEM coupling [9]
(34) 
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 FemBem Coupling
The traditional approach of BEMFEM 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. 56. At external circle in Fig. 56, 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.56. 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 subdomains 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 firsttype) 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 secondtype) 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.
(35) 
For some nonzero 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 subdomain the governing equation is given by [10].
(36) 
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.
(37) 
From Gauss Divergence theorem and (36) it can be stated that
(38) 
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,

(39) 
In (39), the fundamental solution possesses a discontinuity at , . Hence,
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 .
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
(40) 
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).
(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).
(42) 
Where,
[0.5ex]
[0.5ex]
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)
(43) 
(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 rewritten as (45).
(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)
(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
= =
Distance ’S’ can be expressed as
(47) 
Where
[1ex]
[1ex]
As a whole from mathematical point of view
(48) 
and can be written as
(49) 
(50) 
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)
(51) 
Case 2: In this case lies on the boundary element but doesn’t lie on the boundary on which we are integrating. Hence,
(52) 
 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
(53) 
Case 2:
In this case lies on the boundary element but doesn’t lie on the boundary on which we are integrating. Hence
(54) 
Based on multiple parameters of the system, which has been detailed in section II, a onetime 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 FEMBEM 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)
(55) 
In a more abridged form, (55) can be written as
(56) 
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 currentcarrying 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 leftwise 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.
(57) 
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
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% 
A comparative analysis of the magnetic vector potential obtained from Multiphysics software COMSOL using remeshing strategies as per Fig.78 and computed results from MATLAB are shown in Table III obtained with the help of developed BEMFEM code. Relative error analysis graphs as in Fig. 910 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 tensorbased method during comparison of hardware and softwarebased 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 currentcarrying 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 FEMBEM 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 BEMFEM 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.
References
 [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, “Threedimensional transient BEMFEM coupled analysis of electrodynamic levitation problems,” IEEE Trans. Magn., vol. 32, no. 3, pp. 10621065, 1996.
 [3] J. Fetter et al, “Transient BEMFEM coupled analysis of 3D electromechanical systems: a watch stepping motor driven by a thin wire coil,” IEEE Trans. Magn., vol. 34, no. 5, pp. 31543157, 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. 13411345, July 2000.
 [5] A. Jena, S. Sarkar and MashuqunNabi, ”A modified BEMFEM coupling approach for 3D electromagnetic levitation problem,” 2015 Annual IEEE India Conference (INDICON), New Delhi, 2015, pp. 16.
 [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. 20352044, 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 BEMFEM coupling,” IEEE Trans. Magn., vol. 34, no. 5, pp. 30683073, 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 BEMFEM coupling approach,” IEEE Trans. Magn., vol. 35, no. 3, pp. 17931796, 1999.
 [9] S. Kurz, J. Fetzer, G. Lehner and W. M. Rucker, “The application of the BEMFEM coupling method for the treatment of threedimensional nonlinear shielding problems of lowfrequency fields using the example of the TEAM problem 21,” Archive for Electrical Engineering, vol. 80, no. 2, pp. 91104, 1997.
 [10] W. T. Ang, A Beginner’s Course in Boundary Element Method, Universal Publisher, 2007.
Comments
There are no comments yet.