- Email: [email protected]
Boundary element method for multiply connected domains
Journal of Materials Processing Technology 161 (2005) 315–319
Boundary element method for multiply connected domains F. Hantila ∗,1 , M. Vasiliu 2 , M. Maricaru, A. Della Giacomo Universitatea POLITEHNICA din Bucuresti, Bucharest 060042, Romania
Abstract Two procedures in the boundary element method (BEM) for multiply connected domains are reported. The first procedure introduces cuts such that the multiply connected domain is transformed into a simply connected one, so that it allows the use of the scalar potential. Each cut introduces an additional term containing the “hall” current multiplied by the solid angle under which the cut is seen. The second procedure is based on a novel magnetic vector potential formulation. It uses edge elements and tree–cotree spanning. A zero normal component of this vector potential A and the condition for its line integral along any closed path on the boundary are imposed, such that the continuity of the normal component of the magnetic flux density be rigorously satisfied. The unknowns are the tangential components of × A and the tree edge element values. © 2004 Elsevier B.V. All rights reserved. Keywords: Boundary element method; Scalar and vector magnetic potential; Edge elements; Tree–cotree techniques
1. Introduction Boundary element method is often used for solving the magnetic field problems in free space and for obtaining the stiffness matrix for boundary conditions in both magnetostatic and eddy current problems. Media in motion may also be introduced – the stiffness matrix is time dependent, the electromagnetic field equations are expressed in terms of Lagrangean coordinates attached to the bodies. It is well known that the smaller number of unknowns used in scalar potential formulations constitutes an important advantage in comparison with the formulations based on vector potential ones. Unfortunately, use of scalar potential is accompanied by a lower accuracy in the more general case of inhomogeneous media, when large differences in the magnetic permeability values occur. Some difficulties also appear in the case of multiply connected regions. One way to overcome these difficulties is to separate the field component having Γ H dl = i and to use the reduced scalar potential in BEM. Another way, presented in this paper, is to use directly the scalar potential, together with cut contributions. ∗ 1 2
Corresponding author. E-mail address: [email protected] (F. Hantila). Member, IEEE. Student member, IEEE.
0924-0136/$ – see front matter © 2004 Elsevier B.V. All rights reserved. doi:10.1016/j.jmatprotec.2004.07.043
Inherently, a number of additional problems must be solved in order to obtain the stiffness matrix. Several vector potential formulations for the boundary integral method have been reported [1–3]. In most problems the vector potential A and the tangential component of × A are the unknowns and the discretization is performed by the nodal elements. Involving a greater number of unknowns and dealing with singularities of the surface integrals are major disadvantages. The aim of this paper is to introduce a novel BEM algorithm based on magnetic vector potential particularly suitable for multiply connected domains.
2. Scalar potential formulation The integral equation used in scalar BEM formulation for a simply connected domain Ω is given by: 1 ∂V (r ) (R · n ) αV (r) = dS − V (r ) dS + V0 3 ∂Ω R ∂n ∂Ω R (1) where V is the scalar potential; ∂Ω is the boundary of the domain Ω; α is the solid angle under which a small neighbourhood of Ω is seen from the observation point; r and r
316
F. Hantila et al. / Journal of Materials Processing Technology 161 (2005) 315–319
placed outside the facet. Supposing in addition that ∂V(r )/∂n is constant on each facet, Eq. (1) becomes: F ∂V dS αV (r) = ∂n k ωk R k=1
+
Fig. 1. Cut for a multiply connected domain.
F k=1
are the position vectors of the observation and source points, respectively; R = r − r ; n is the outward normal unit vector; V0 is the scalar potential of the given distribution of sources. For multiply connected domains, the right hand side of Eq. (1) includes V0 = k βk ik , where βk is the solid angle under which the cut σ k is seen from the observation point and ik is the “hall” current (Fig. 1).
3. Numerical solution of the scalar potential integral equation The boundary ∂Ω is approximated by a polyhedral surface with triangular facets. Assuming an arbitrary origin in the local coordinates of the facet the scalar potential can be interpolated linearly in the form (Fig. 2): V (ρ) = V- + T · ρ
(2)
where T =−
V- =
V1 R23 + V2 R31 + V3 R12 × n 2S
V1 S + V2 S + V3 S S
(3)
(4)
S is the facet area, V1 , V2 , V3 are the potentials of the facet nodes and V- is the potential in the origin, S is the area of the triangle defined by the origin and the edge R23 opposite to the node 1; it may have negative value, if the origin is
Fig. 2. A facet.
γk V- k − nk
R (Tk · ρ) 3 dS + V0 (5) R ωk
where γ k is the solid angle under which the facet is seen from the observation point. All integrals in Eq. (5) can be analytically evaluated, by dividing the facet into three triangles defined by the local origin (the projection of the observation point onto the facet) and the sides of the facet (Fig. 2). If the normal derivatives (∂V/∂n )k are known, we place the observation point in all nodes of the surface. We thus obtain a linear system with potentials Vk as unknowns. Since Eq. (1) is a Fredholm integral equation of the second kind, the linear system is well conditioned. The nodes belonging to the cut are treated differently: when they are taken into account by the facets placed above the cut, their potentials are increased by i/2 in comparison with the potentials of the same nodes as an observation point. For facets under the cut, their potentials are decreased with i/2.
4. Vector potential formulation The vector potential A satisfies the equation n R αA(r) = − × (∇ × A(r )) − 3 R R ∂Ω R × (n × A(r )) + 3 (n · A(r )) dS + 4πA0 (r) R (6) where A0 is the vector potential produced by the given distribution of currents. A verifies the gauge condition ·A. Since the tangential component At on the boundary uniquely defines the vector potential in the region, the normal component An and the tangential component ( × A)t are determined by At and thus Eq. (5) can be regarded as an equation in the unknowns An and ( × A)t , when At is given. The BEM based on Eq. (6) has two drawbacks [2]: it requires a large number of unknowns and presents a singularity in the integral of the third term. In what follows, we overcome the latter difficulty by imposing An = 0. Accordingly Eq. (6) transforms to: n αA(r) = − × (∇ × A(r )) dS ∂Ω R R + × (n × A(r )) dS + A0 (7) 3 ∂Ω R Obviously, the condition satisfied by the tangential component At cannot be enforced any longer. Instead, we impose a weaker condition, namely just that for the line integral of
F. Hantila et al. / Journal of Materials Processing Technology 161 (2005) 315–319
A along any closed path on ∂Ω. This represents in fact the actual boundary condition imposed on Bn , which determines uniquely the magnetic field in Ω. It should be remarked that this magnetic vector potential A is uniquely defined in Ω since ·A and × A at any point in Ω, as well as An on its boundary, are given.
5. Numerical solution of the vector potential integral equation
2π (−τ1i − τ2i + 2τ3i ) 3 F
(∇ × A)tk · (nk × R12i )
k=1
where Bnk is the normal component of the magnetic flux density, assumed to be constant over each facet ωk ; ρk is the barycentric position vector of an arbitrary point on ωk ; Sk is the area of ωk ; and τik , i = 1, 2, 3 is the line integral of A (i.e. the edge value of A) along the edge i of ωk opposed to the node of barycentric position vector ρik . Assuming a linear variation of A within a tetrahedral volume element, its value at any point can be expressed in terms of the edge values associated with the corresponding tetrahedron. Thus, the potential in Eq. (8), on the boundary, can be directly coupled with the vector potential outside Ω by employing edge elements. Using a tree–cotree spanning [5] of the graph of the edges on the boundary, the line integrals of A on the closed loops associated with the cotree edges are determined by the values of Bn from (9)
where τ C and τ B are vectors containing the line integrals of A on the cotree edges and on the tree edges, respectively; ΛCB is the tree–cotree connection matrix, and ϕC is the vector containing the magnetic fluxes on the surfaces bounded by the closed loops associated with the cotree edges. Eq. (9) allows to express the cotree edge values in terms of the tree edge values, which means that the number of edge unknowns is reduced to only the number N − 1 of tree edges, where N is the total number of nodes. Assuming that ( × A)t is also constant on each facet, we have in the local coordinates of the respective triangle ωk (∇ × A)tk = uk R12k + vk R23k
Bringing the observation point in Eq. (7) to the centres Gi of the triangles ωi on the boundary and projecting Eq. (7) on the edges R12i and R23i yields 2F equations. Taking into account Eq. (8), with ρk = 0, the projection on the edge R12i , for example, gives
=
Evaluation of the line integrals of A along closed paths on the boundary, required in the numerical solution of Eq. (7), can be performed efficiently by using edge elements [4]. The boundary ∂Ω is approximated by a polyhedral surface with triangular facets. On each such facet the tangential component of the vector potential can be interpolated linearly in the form 1 Bnk Atk = (τ1k ρ1k + τ2k ρ2k + τ3k ρ3k ) − ρk × nk 2Sk 2 (8)
τ C + ΛCB τ B = ϕC
317
(10)
where R12k = ρ2k − ρ1k , R23k = ρ3k − ρ2k , and uk and vk are unknowns associated with each of the F facets of the boundary.
+
ωk
dS R
F k=1k =i
Bnk (nk · R12i ) W k − RGik 2
dl Bnk (nk · RGik )R12i + 2 ∂ωk R Bnk − (nk × R12i ) W k − RGik Ωk + 4πA012i 2 (11) where Ωk is the solid angle under which the facet ωk is seen from the centre of the facet ωi and R = |rGi − r |,
RGik ≡ rGi − r Gk ,
A012i = A0 (r Gi ) · R12i
Wk ≡
1 (τ1 ρ1 + τ2k ρ2k + τ3k ρ3k ) 2Sk k k
It should be noticed that all the integrals in Eq. (11) could be evaluated analytically. In order to obtain the rest of N − 1 equations, there are two possibilities. One consists in placing the observation point at the centres of the facets connected to a node and imposing the condition that the sum of the fluxes of A through these facets be equal to zero, assuming that the normal components of A are constant over the facets. We enforce this condition for N − 1 nodes. The procedure amounts to solving a Fredholm equation of the first kind. The other possibility is to place the observation point on each edge i of the N − 1 edges of the tree and to impose the value of the line integral of A along the edge i to be τ i . The latter procedure is more convenient since it yields a Fredholm equation of the second kind, which can be written in the form αi τi =
F
(∇ × A)tk · (nk × ei )
k=1
+
Ti
F Ti k=1k∈{i} /
ωk
dS dli R
Bnk (nk · ei ) W k − RGk 2
Bnk + (nk · RGk )ei 2
∂ωk
dl R
318
F. Hantila et al. / Journal of Materials Processing Technology 161 (2005) 315–319
Bnk − (nk × ei ) W k − RGk Ωki dli 2 + 4π A0 · ei dli Ti
(12)
where ei is the unit vector specifying the orientation of the edge i, {i} is the set of facets containing the edge i, Ωki is the solid angle under which the facet ωk is seen from the integration point on the edge i. The integral over Ti along an edge i is computed numerically only when the triangle ωk is close to this edge. Otherwise, the observation point is placed in the middle of the edge and the product of the respective integrand and the length of that edge approximate the integral. The number of edge unknowns in Eqs. (11) and (12) is reduced to the number N − 1 of tree edges since the cotree edge values are eliminated by using Eq. (9).
Fig. 4. Normalized inductance L = L/(µ0 /4π) of a toroid of mean radius of 1 m and cross-section radius: () 1 cm; () 2 cm.
ductance has been computed for different radii of the toroid cross-section and for various refinements of the mesh. The normalized values of the inductance, L = L/(µ0 /4π), are plotted in Fig. 4 When the number of nodes is 169, for instance, we have F = 338 facets and only 846 unknowns, including the unknown values for the edges c and c , instead of a total of 1183 unknowns.
6. Results
6.2. Coupled conducting toroids vector potential
6.1. Perfectly conducting toroid vector potential
A more complex example consists of two perfectly conducting toroids with mutual arbitrary orientation (Fig. 5). The toroids have 1 m mean radius and circular cross-sections. The centre of the second toroid has the Cartesian coordinates (0; 2.5 m; 0.5 m). The origin is at the centre of the first toroid; z-axis being the axis of revolution. The inclination of the second toroid is obtained by a rotation of 15◦ around x-axis and by a second rotation by 30◦ around y-axis. Computed values for the normalized self-inductance of the first toroid and mutual inductance are plotted in Figs. 6 and 7, respectively. If the number of nodes per toroid is 400, for instance, then the total number of nodes is 800, and F = 1600 leading to 4002 unknowns, the total number of 2400 edge unknowns is reduced to 802 unknowns only, including the unknown values for the edges c and c . From Figs. 4 and 7 we observe that self-inductances are very close for the same discretization, as expected since the mutual inductance is comparatively small.
To illustrate the proposed method for multiply connected regions, we consider the region on the outside of a perfectly conducting toroid of mean radius of 1 m and circular crosssection. The aim is to compute the self-inductance of the toroid by vector potential formulation. The corresponding current carried by the toroid was determined assuming a magnetic flux through the toroid cut of 1 Wb. In this example, the procedure based on the Fredholm equation of the second kind was used. The set of the tree edges is completed with two cotree edges: the first one (c in Fig. 3) closes a loop along the toroid through which the magnetic flux, i.e. the line integral of A, is 1 Wb, and the second one (c in Fig. 3) closes a loop around the toroid, along which the line integral of A is zero. After computing the values of ( × A)t on the facets, the line integral of ( × A)t along a path around the toroid was determined, which gives the current carried by it. The in-
Fig. 3. Boundary mesh for a toroid: solid lines represent tree edges and thin lines, cotree edges.
Fig. 5. Two identical toroids of mean radius of 1 m and distance between centers of 2.55 m.
F. Hantila et al. / Journal of Materials Processing Technology 161 (2005) 315–319
319
with the middle cut of the previous second toroid (Fig. 5) gives Φ/(µ0 /4π) = 0.0214151, for 1 A. This value is comparable with normalized mutual inductance obtained by vector potential formulation for two toroids configuration.
7. Conclusions Fig. 6. Normalized self-inductance L = L/(µ0 /4π) of the first toroid; toroid cross-section radius: () 1 cm; () 2 cm.
Fig. 7. Normalized mutual inductance M = M/(µ0 /4π) for the two toroids system (Fig. 6); toroid cross-section radius: () 1 cm; () 2 cm.
6.3. Scalar potential For one toroid configuration, using scalar potential formulation, we obtain the equipotential lines shown in Fig. 8. The number of nodes and facets were 420 and 840, respectively. The magnetic flux through a surface that coincides
Both scalar and vector potentials can equally be used in BEM formulation for the magnetic field computation in multiply connected domains. Scalar potential approach leads to a smaller number of unknowns, but requires the derivative of Eq. (1) in order to calculate the magnetic fluxes. Unfortunately, singularities in integrals that occur when we have to compute the cut magnetic flux for inductance evaluation must be resolved. A future goal is to overcome this problem. Vector potential formulation avoids singularities. The large number of unknowns can be reduced by using edge elements and tree–cotree spanning so that both the computing time and the number of unknowns decrease substantially. Instead of imposing the boundary condition for the tangential component of A, the continuity of the normal component of B is imposed, which constitutes another advantage of the proposed procedure. Since the numerical results do not depend on the choice of the tree, the mesh topology can be adapted conveniently. In addition, the integrals defining the matrix entries can be evaluated analytically, leading to an increase of the accuracy of the results and the computational efficiency.
References
Fig. 8. Equipotential lines.
[1] S. Nath, W. Lord, T.J. Rudolphi, Three dimensional hybrid finiteboundary element model for eddy current NDE, IEEE Trans. Magn. 29 (1993) 1853–1856. [2] T. Onuki, S. Wakao, Novel boundary element analysis for 3-D eddy current problems, IEEE Trans. Magn. 29 (1993) 1520–1523. [3] J. Fetzer, St. Kurz, G. Lehner, Comparison of analytical and numerical integration technique for the boundary integrals in the BEM-FEM coupling considering TEAM workshop problem no.13, IEEE Trans. Magn. 33 (1997) 1227–1230. [4] R. Albanese, G. Rubinacci, Integral formulation for 3D eddy current computation using edge elements, IEE Proc. 135A (1988) 457– 462. [5] F. Hantila, I.R. Ciric, Magnetic Vector Potential Tree Edge Values for Boundary Elements, CEFC2002, June 16–19, 2002, Perugia, Italy.
Boundary element method for multiply connected domains F. Hantila ∗,1 , M. Vasiliu 2 , M. Maricaru, A. Della Giacomo Universitatea POLITEHNICA din Bucuresti, Bucharest 060042, Romania
Abstract Two procedures in the boundary element method (BEM) for multiply connected domains are reported. The first procedure introduces cuts such that the multiply connected domain is transformed into a simply connected one, so that it allows the use of the scalar potential. Each cut introduces an additional term containing the “hall” current multiplied by the solid angle under which the cut is seen. The second procedure is based on a novel magnetic vector potential formulation. It uses edge elements and tree–cotree spanning. A zero normal component of this vector potential A and the condition for its line integral along any closed path on the boundary are imposed, such that the continuity of the normal component of the magnetic flux density be rigorously satisfied. The unknowns are the tangential components of × A and the tree edge element values. © 2004 Elsevier B.V. All rights reserved. Keywords: Boundary element method; Scalar and vector magnetic potential; Edge elements; Tree–cotree techniques
1. Introduction Boundary element method is often used for solving the magnetic field problems in free space and for obtaining the stiffness matrix for boundary conditions in both magnetostatic and eddy current problems. Media in motion may also be introduced – the stiffness matrix is time dependent, the electromagnetic field equations are expressed in terms of Lagrangean coordinates attached to the bodies. It is well known that the smaller number of unknowns used in scalar potential formulations constitutes an important advantage in comparison with the formulations based on vector potential ones. Unfortunately, use of scalar potential is accompanied by a lower accuracy in the more general case of inhomogeneous media, when large differences in the magnetic permeability values occur. Some difficulties also appear in the case of multiply connected regions. One way to overcome these difficulties is to separate the field component having Γ H dl = i and to use the reduced scalar potential in BEM. Another way, presented in this paper, is to use directly the scalar potential, together with cut contributions. ∗ 1 2
Corresponding author. E-mail address: [email protected] (F. Hantila). Member, IEEE. Student member, IEEE.
0924-0136/$ – see front matter © 2004 Elsevier B.V. All rights reserved. doi:10.1016/j.jmatprotec.2004.07.043
Inherently, a number of additional problems must be solved in order to obtain the stiffness matrix. Several vector potential formulations for the boundary integral method have been reported [1–3]. In most problems the vector potential A and the tangential component of × A are the unknowns and the discretization is performed by the nodal elements. Involving a greater number of unknowns and dealing with singularities of the surface integrals are major disadvantages. The aim of this paper is to introduce a novel BEM algorithm based on magnetic vector potential particularly suitable for multiply connected domains.
2. Scalar potential formulation The integral equation used in scalar BEM formulation for a simply connected domain Ω is given by: 1 ∂V (r ) (R · n ) αV (r) = dS − V (r ) dS + V0 3 ∂Ω R ∂n ∂Ω R (1) where V is the scalar potential; ∂Ω is the boundary of the domain Ω; α is the solid angle under which a small neighbourhood of Ω is seen from the observation point; r and r
316
F. Hantila et al. / Journal of Materials Processing Technology 161 (2005) 315–319
placed outside the facet. Supposing in addition that ∂V(r )/∂n is constant on each facet, Eq. (1) becomes: F ∂V dS αV (r) = ∂n k ωk R k=1
+
Fig. 1. Cut for a multiply connected domain.
F k=1
are the position vectors of the observation and source points, respectively; R = r − r ; n is the outward normal unit vector; V0 is the scalar potential of the given distribution of sources. For multiply connected domains, the right hand side of Eq. (1) includes V0 = k βk ik , where βk is the solid angle under which the cut σ k is seen from the observation point and ik is the “hall” current (Fig. 1).
3. Numerical solution of the scalar potential integral equation The boundary ∂Ω is approximated by a polyhedral surface with triangular facets. Assuming an arbitrary origin in the local coordinates of the facet the scalar potential can be interpolated linearly in the form (Fig. 2): V (ρ) = V- + T · ρ
(2)
where T =−
V- =
V1 R23 + V2 R31 + V3 R12 × n 2S
V1 S + V2 S + V3 S S
(3)
(4)
S is the facet area, V1 , V2 , V3 are the potentials of the facet nodes and V- is the potential in the origin, S is the area of the triangle defined by the origin and the edge R23 opposite to the node 1; it may have negative value, if the origin is
Fig. 2. A facet.
γk V- k − nk
R (Tk · ρ) 3 dS + V0 (5) R ωk
where γ k is the solid angle under which the facet is seen from the observation point. All integrals in Eq. (5) can be analytically evaluated, by dividing the facet into three triangles defined by the local origin (the projection of the observation point onto the facet) and the sides of the facet (Fig. 2). If the normal derivatives (∂V/∂n )k are known, we place the observation point in all nodes of the surface. We thus obtain a linear system with potentials Vk as unknowns. Since Eq. (1) is a Fredholm integral equation of the second kind, the linear system is well conditioned. The nodes belonging to the cut are treated differently: when they are taken into account by the facets placed above the cut, their potentials are increased by i/2 in comparison with the potentials of the same nodes as an observation point. For facets under the cut, their potentials are decreased with i/2.
4. Vector potential formulation The vector potential A satisfies the equation n R αA(r) = − × (∇ × A(r )) − 3 R R ∂Ω R × (n × A(r )) + 3 (n · A(r )) dS + 4πA0 (r) R (6) where A0 is the vector potential produced by the given distribution of currents. A verifies the gauge condition ·A. Since the tangential component At on the boundary uniquely defines the vector potential in the region, the normal component An and the tangential component ( × A)t are determined by At and thus Eq. (5) can be regarded as an equation in the unknowns An and ( × A)t , when At is given. The BEM based on Eq. (6) has two drawbacks [2]: it requires a large number of unknowns and presents a singularity in the integral of the third term. In what follows, we overcome the latter difficulty by imposing An = 0. Accordingly Eq. (6) transforms to: n αA(r) = − × (∇ × A(r )) dS ∂Ω R R + × (n × A(r )) dS + A0 (7) 3 ∂Ω R Obviously, the condition satisfied by the tangential component At cannot be enforced any longer. Instead, we impose a weaker condition, namely just that for the line integral of
F. Hantila et al. / Journal of Materials Processing Technology 161 (2005) 315–319
A along any closed path on ∂Ω. This represents in fact the actual boundary condition imposed on Bn , which determines uniquely the magnetic field in Ω. It should be remarked that this magnetic vector potential A is uniquely defined in Ω since ·A and × A at any point in Ω, as well as An on its boundary, are given.
5. Numerical solution of the vector potential integral equation
2π (−τ1i − τ2i + 2τ3i ) 3 F
(∇ × A)tk · (nk × R12i )
k=1
where Bnk is the normal component of the magnetic flux density, assumed to be constant over each facet ωk ; ρk is the barycentric position vector of an arbitrary point on ωk ; Sk is the area of ωk ; and τik , i = 1, 2, 3 is the line integral of A (i.e. the edge value of A) along the edge i of ωk opposed to the node of barycentric position vector ρik . Assuming a linear variation of A within a tetrahedral volume element, its value at any point can be expressed in terms of the edge values associated with the corresponding tetrahedron. Thus, the potential in Eq. (8), on the boundary, can be directly coupled with the vector potential outside Ω by employing edge elements. Using a tree–cotree spanning [5] of the graph of the edges on the boundary, the line integrals of A on the closed loops associated with the cotree edges are determined by the values of Bn from (9)
where τ C and τ B are vectors containing the line integrals of A on the cotree edges and on the tree edges, respectively; ΛCB is the tree–cotree connection matrix, and ϕC is the vector containing the magnetic fluxes on the surfaces bounded by the closed loops associated with the cotree edges. Eq. (9) allows to express the cotree edge values in terms of the tree edge values, which means that the number of edge unknowns is reduced to only the number N − 1 of tree edges, where N is the total number of nodes. Assuming that ( × A)t is also constant on each facet, we have in the local coordinates of the respective triangle ωk (∇ × A)tk = uk R12k + vk R23k
Bringing the observation point in Eq. (7) to the centres Gi of the triangles ωi on the boundary and projecting Eq. (7) on the edges R12i and R23i yields 2F equations. Taking into account Eq. (8), with ρk = 0, the projection on the edge R12i , for example, gives
=
Evaluation of the line integrals of A along closed paths on the boundary, required in the numerical solution of Eq. (7), can be performed efficiently by using edge elements [4]. The boundary ∂Ω is approximated by a polyhedral surface with triangular facets. On each such facet the tangential component of the vector potential can be interpolated linearly in the form 1 Bnk Atk = (τ1k ρ1k + τ2k ρ2k + τ3k ρ3k ) − ρk × nk 2Sk 2 (8)
τ C + ΛCB τ B = ϕC
317
(10)
where R12k = ρ2k − ρ1k , R23k = ρ3k − ρ2k , and uk and vk are unknowns associated with each of the F facets of the boundary.
+
ωk
dS R
F k=1k =i
Bnk (nk · R12i ) W k − RGik 2
dl Bnk (nk · RGik )R12i + 2 ∂ωk R Bnk − (nk × R12i ) W k − RGik Ωk + 4πA012i 2 (11) where Ωk is the solid angle under which the facet ωk is seen from the centre of the facet ωi and R = |rGi − r |,
RGik ≡ rGi − r Gk ,
A012i = A0 (r Gi ) · R12i
Wk ≡
1 (τ1 ρ1 + τ2k ρ2k + τ3k ρ3k ) 2Sk k k
It should be noticed that all the integrals in Eq. (11) could be evaluated analytically. In order to obtain the rest of N − 1 equations, there are two possibilities. One consists in placing the observation point at the centres of the facets connected to a node and imposing the condition that the sum of the fluxes of A through these facets be equal to zero, assuming that the normal components of A are constant over the facets. We enforce this condition for N − 1 nodes. The procedure amounts to solving a Fredholm equation of the first kind. The other possibility is to place the observation point on each edge i of the N − 1 edges of the tree and to impose the value of the line integral of A along the edge i to be τ i . The latter procedure is more convenient since it yields a Fredholm equation of the second kind, which can be written in the form αi τi =
F
(∇ × A)tk · (nk × ei )
k=1
+
Ti
F Ti k=1k∈{i} /
ωk
dS dli R
Bnk (nk · ei ) W k − RGk 2
Bnk + (nk · RGk )ei 2
∂ωk
dl R
318
F. Hantila et al. / Journal of Materials Processing Technology 161 (2005) 315–319
Bnk − (nk × ei ) W k − RGk Ωki dli 2 + 4π A0 · ei dli Ti
(12)
where ei is the unit vector specifying the orientation of the edge i, {i} is the set of facets containing the edge i, Ωki is the solid angle under which the facet ωk is seen from the integration point on the edge i. The integral over Ti along an edge i is computed numerically only when the triangle ωk is close to this edge. Otherwise, the observation point is placed in the middle of the edge and the product of the respective integrand and the length of that edge approximate the integral. The number of edge unknowns in Eqs. (11) and (12) is reduced to the number N − 1 of tree edges since the cotree edge values are eliminated by using Eq. (9).
Fig. 4. Normalized inductance L = L/(µ0 /4π) of a toroid of mean radius of 1 m and cross-section radius: () 1 cm; () 2 cm.
ductance has been computed for different radii of the toroid cross-section and for various refinements of the mesh. The normalized values of the inductance, L = L/(µ0 /4π), are plotted in Fig. 4 When the number of nodes is 169, for instance, we have F = 338 facets and only 846 unknowns, including the unknown values for the edges c and c , instead of a total of 1183 unknowns.
6. Results
6.2. Coupled conducting toroids vector potential
6.1. Perfectly conducting toroid vector potential
A more complex example consists of two perfectly conducting toroids with mutual arbitrary orientation (Fig. 5). The toroids have 1 m mean radius and circular cross-sections. The centre of the second toroid has the Cartesian coordinates (0; 2.5 m; 0.5 m). The origin is at the centre of the first toroid; z-axis being the axis of revolution. The inclination of the second toroid is obtained by a rotation of 15◦ around x-axis and by a second rotation by 30◦ around y-axis. Computed values for the normalized self-inductance of the first toroid and mutual inductance are plotted in Figs. 6 and 7, respectively. If the number of nodes per toroid is 400, for instance, then the total number of nodes is 800, and F = 1600 leading to 4002 unknowns, the total number of 2400 edge unknowns is reduced to 802 unknowns only, including the unknown values for the edges c and c . From Figs. 4 and 7 we observe that self-inductances are very close for the same discretization, as expected since the mutual inductance is comparatively small.
To illustrate the proposed method for multiply connected regions, we consider the region on the outside of a perfectly conducting toroid of mean radius of 1 m and circular crosssection. The aim is to compute the self-inductance of the toroid by vector potential formulation. The corresponding current carried by the toroid was determined assuming a magnetic flux through the toroid cut of 1 Wb. In this example, the procedure based on the Fredholm equation of the second kind was used. The set of the tree edges is completed with two cotree edges: the first one (c in Fig. 3) closes a loop along the toroid through which the magnetic flux, i.e. the line integral of A, is 1 Wb, and the second one (c in Fig. 3) closes a loop around the toroid, along which the line integral of A is zero. After computing the values of ( × A)t on the facets, the line integral of ( × A)t along a path around the toroid was determined, which gives the current carried by it. The in-
Fig. 3. Boundary mesh for a toroid: solid lines represent tree edges and thin lines, cotree edges.
Fig. 5. Two identical toroids of mean radius of 1 m and distance between centers of 2.55 m.
F. Hantila et al. / Journal of Materials Processing Technology 161 (2005) 315–319
319
with the middle cut of the previous second toroid (Fig. 5) gives Φ/(µ0 /4π) = 0.0214151, for 1 A. This value is comparable with normalized mutual inductance obtained by vector potential formulation for two toroids configuration.
7. Conclusions Fig. 6. Normalized self-inductance L = L/(µ0 /4π) of the first toroid; toroid cross-section radius: () 1 cm; () 2 cm.
Fig. 7. Normalized mutual inductance M = M/(µ0 /4π) for the two toroids system (Fig. 6); toroid cross-section radius: () 1 cm; () 2 cm.
6.3. Scalar potential For one toroid configuration, using scalar potential formulation, we obtain the equipotential lines shown in Fig. 8. The number of nodes and facets were 420 and 840, respectively. The magnetic flux through a surface that coincides
Both scalar and vector potentials can equally be used in BEM formulation for the magnetic field computation in multiply connected domains. Scalar potential approach leads to a smaller number of unknowns, but requires the derivative of Eq. (1) in order to calculate the magnetic fluxes. Unfortunately, singularities in integrals that occur when we have to compute the cut magnetic flux for inductance evaluation must be resolved. A future goal is to overcome this problem. Vector potential formulation avoids singularities. The large number of unknowns can be reduced by using edge elements and tree–cotree spanning so that both the computing time and the number of unknowns decrease substantially. Instead of imposing the boundary condition for the tangential component of A, the continuity of the normal component of B is imposed, which constitutes another advantage of the proposed procedure. Since the numerical results do not depend on the choice of the tree, the mesh topology can be adapted conveniently. In addition, the integrals defining the matrix entries can be evaluated analytically, leading to an increase of the accuracy of the results and the computational efficiency.
References
Fig. 8. Equipotential lines.
[1] S. Nath, W. Lord, T.J. Rudolphi, Three dimensional hybrid finiteboundary element model for eddy current NDE, IEEE Trans. Magn. 29 (1993) 1853–1856. [2] T. Onuki, S. Wakao, Novel boundary element analysis for 3-D eddy current problems, IEEE Trans. Magn. 29 (1993) 1520–1523. [3] J. Fetzer, St. Kurz, G. Lehner, Comparison of analytical and numerical integration technique for the boundary integrals in the BEM-FEM coupling considering TEAM workshop problem no.13, IEEE Trans. Magn. 33 (1997) 1227–1230. [4] R. Albanese, G. Rubinacci, Integral formulation for 3D eddy current computation using edge elements, IEE Proc. 135A (1988) 457– 462. [5] F. Hantila, I.R. Ciric, Magnetic Vector Potential Tree Edge Values for Boundary Elements, CEFC2002, June 16–19, 2002, Perugia, Italy.