Integration of the modified double layer potential of the vector boundary element method for eddy current problems

Abstract

The boundary element method for the eddy current problem (BEM-ECP) was proposed in a number of papers and is applicable to important tasks such as the problem of inductive heating and transmission of electromagnetic energy. BEM-ECP requires the construction of a system of linear algebraic equations in which the matrix is inherently dense and is constructed out of element matrices. For the process of the element matrix computation, two cases are normally considered: far-field interaction and near-field interaction, because the construction of element matrices requires integration of a singular function. In this article, we suggest a transform that allows computing the matrix components of the near-singular interaction part while implementing only the single and double layer potentials. The previously suggested modified double layer potential (MDLP) can be integrated by means of this transform, which simplifies the program implementation of BEM-ECP significantly. Solving model problems, we analyse the drawbacks of the previously suggested approach. This analysis includes the proof of the MDLP singularity that makes the integration of this potential a rather difficult task without the help of our transform. The previously suggested approach does not work well with surfaces that are not smooth. Our approach does consider such cases, which is its main advantage. We demonstrate this on the model problems with known analytical solutions.

1 Introduction

Boundary element method for the eddy current problem (BEM-ECP) is highly useful for the applications where the unbounded regions are involved and the electromagnetic energy transmission plays an important role. This method allows one to solve the Maxwell equations when eddy currents are taken into account. BEM-ECP was initially introduced in [17] written by Jörg Ostrowski and was also described in [7] by Jens Breuer for the problem of inductive heating and cooling.

BEM-ECP incorporates the Stratton–Chu formula first introduced in [21] as a valuable part of the electromagnetic wave diffraction theory. The most recent paper considering the same foundation for the boundary element method approach is [2] by Yang et al. However, the optimisation method mentioned in [2] deals with the case of the far-field interactions and the paper does not cover the computations of integrals with singular and near-singular issues. However, the only work that highlighted this approach, as we know, is Ostrowski’s [17] where the modified double layer potential was suggested to be applied to construct the components of the SLAE matrix of BEM-ECP for the near-singular cases. The convergence theorems were proved in [7] and [17] both for BEM-ECP and for the coupled approach used together with the finite element method correspondingly. Hence, this topic is not covered here.

In this article, we introduce a conceptually different method from what was initially proposed in the aforementioned works and highlight the drawbacks of the previously suggested method in [17]. We prove that modified double layer potential (MDLP) can be singular and hence it cannot be integrated by the means of the quadrature rules that do not take into account this singularity. We also show that MDLP is redundant and the matrix of BEM-ECP can be constructed for the singular and near-singular cases via the single and double layer potentials only.

The scalar single and double layer potentials can be extended to the vector-valued case in an obvious component-wise fashion. This leads to the significant simplification of the BEM-ECP implementation. A wide variety of methods consider the integration of these potentials, such as: Quadrature by Expansion (QBX) highlighted in [24, 25] and [1], the singularity extraction method described in [12] and the singularity cancellation methods [6, 23]. All of these methods can be applied to BEM-ECP with the suggested modification. In this article, we present the results for model problems to demonstrate the benefits of our approach.

We must mention that here we only consider the case of $\sigma \neq 0$ everywhere (conductivity is not zero in every domain). If this assumption is not true for some domains, we refer to [19] where we use the coupled approach with vector and scalar potentials. To consider the vector potential in [19], we use BEM-ECP with the enhancements we suggest in the present work.

2 Integral operators for eddy current problem

In this section, we give a brief description of BEM-ECP and the main integral operators involved in the computations. Also, we describe a special case of one integral operator defined in [17] and called MDLP.

2.1 Eddy current model and the boundary integral problem

Let us consider the case of domains with constant scalar material parameters of conductivity and magnetic permeability. The equations of the eddy current model for the time-harmonic electromagnetic field in every such domain have a form [7, 17]:

(2.1) \begin{equation}\nabla\times\nabla\times\vec{E}(\vec{x})+k_j^2 \vec{E}(\vec{x})=i \omega \mu_j \vec{I} (\vec{x}), \vec{x}\in\Omega_j, \Omega_j\subset\mathbb{R}^3, j = \overline{1, N},\end{equation}

where $\vec{E}$ is the electric field, $\left\{\Omega_j\right\}_{j=1}^{N}$ is a set of domains in $\mathbb{R}^3$ , $\Gamma_{m,j}$ is a piece of boundary surface between $\Omega_m$ and $\Omega_j$ , $k_j$ – complex wave number defined in $\Omega_j$ , $\omega$ – angular frequency, $\vec{n}_m$ – an exterior normal vector defined on $\partial\Omega_m$ , $\mu_j$ is the magnetic permeability, $\vec{I}$ – exciting current and i – imaginary one. Two additional requirements imposing continuity of tangential components of $\vec{E}$ on $\Gamma_{m,j}$ as well as continuity of tangential components of $\vec{H}$ are also applied to the model [7]. Note that in $\Omega_j$ the relation $\vec{H} = {i \over {\omega \mu_j}}\nabla\times\vec{E}$ is valid.

The wave number $k_j$ can be expressed in terms of magnetic permeability $\mu_j$ and conductivity $\sigma_j$ of $\Omega_j$ as follows:

(2.2) \begin{equation}k_j = \sqrt{i\omega \mu_j \sigma_j},\end{equation}

where the square root is taken in the sense that guarantees the positive real values for a purely imaginary argument. Note that $k_j \neq 0$ for every value of j.

Let a continuously differentiable vector function $\vec{U}$ be defined in a three-dimensional domain $\Omega$ with at least partially smooth boundary $\Gamma$ . The normal vector $\vec{n}$ is defined on $\Gamma$ and directed to the exterior of $\Omega$ . The Neumann trace operator $\gamma_{N}^{\Omega}$ , the Dirichlet trace operator $\gamma_{D}^{\Omega}$ and the normal trace $\gamma_{n}^{\Omega}$ are defined as follows [17, 7, 10]:

(2.3) \begin{equation}\gamma_{N}^{\Omega}\vec{U}(\vec{x}) =\lim_{\vec{r}\in\Omega,\vec{r}\to\vec{x}\in\Gamma}\left(\nabla_{\vec{r}}\times\vec{U}\left(\vec{r}\right)\right)\times\vec{n}(\vec{x}),\end{equation}

(2.4) \begin{equation}\gamma_{D}^{\Omega}\vec{U}(\vec{x}) =\lim_{\vec{r}\in\Omega,\vec{r}\to\vec{x}\in\Gamma}\left(\vec{n}(\vec{x})\times\vec{U}\left(\vec{r}\right)\right)\times\vec{n}(\vec{x}),\end{equation}

(2.5) \begin{equation}\gamma_{n}^{\Omega}\vec{U}(\vec{x}) =\lim_{\vec{r}\in\Omega,\vec{r}\to\vec{x}\in\Gamma}\left(\vec{n}(\vec{x})\cdot\vec{U}\left(\vec{r}\right)\right),\end{equation}

where the subscript of $\nabla_{\vec{r}}$ means that the function-argument is differentiated by $\vec{r}$ .

The foundation of the boundary element method for eddy current problem is the Stratton–Chu representation formula [7, p. 46] that gives the solution to (2.1) in terms of the mentioned trace operators:

(2.6) \begin{align}\vec{E}(\vec{x})&=\nabla_{\vec{x}}\times\int_{\Gamma}G_k(\vec{x},\vec{y})\left(\gamma_{D}^{\Omega}\vec{E}(\vec{y})\times\vec{n}(\vec{y})\right)\!ds_{\vec{y}} + \nonumber\\[3pt]& \quad +\int_{\Gamma}G_k(\vec{x},\vec{y})\left(\gamma_{N}^{\Omega}\vec{E}(\vec{y})\right)\!ds_{\vec{y}} + \nonumber\\[3pt]& \quad +\nabla_{\vec{x}}\int_{\Gamma}G_k(\vec{x},\vec{y})\left(\gamma_{n}^{\Omega}\vec{E}(\vec{y})\right)\!ds_{\vec{y}}.\end{align}

where $ds_{\vec{y}}$ means that integration variable is $\vec{y}$ and the integral is taken along a surface to which $\vec{y}$ belongs, similarly, $\nabla_{\vec{x}}$ designates the differentiation by the parameter $\vec{x}$ ; $G_k$ is the singular function of the Helmholtz equation [20, 21]:

(2.7) \begin{equation}G_k(\vec{x},\vec{y}) = \frac{e^{-k\|\vec{x} - \vec{y}\|}}{4\pi\|\vec{x} - \vec{y}\|}.\end{equation}

It must be noted that when $k\neq 0$ , the following relation between the normal trace and the Neumann trace takes place [7, p. 34, formula (4.10)]:

(2.8) \begin{equation}\gamma_{n}^{\Omega}\vec{E}(\vec{x}) =-{1 \over k^2} \nabla_{\vec{x}} \cdot \gamma_{N}^{\Omega}\vec{E}(\vec{x}),\end{equation}

so the two trace operators are connected. The relation (2.8) is valid for every $\vec{E}$ satisfying (2.1).

Let a vector function $\vec{u}(\vec{x})$ be defined for $\vec{x}\in\Gamma$ and $\vec{u}(\vec{x})\perp\vec{n}(\vec{x})$ for almost all $\vec{x}\in\Gamma$ . Using the trace operators one can write the integral operators $A_k$ and $B_k$ in accordance with [7] (see [7, p. 47, formula (5.10)], [7, p. 48, formula (5.11)] for the operator $B_k$ definitions; we define operator $A_k$ , that is denoted in [7] as $\tilde{A}_{\kappa}^{l}$ , in accordance with [7, p. 51, formula (5.26)]):

(2.9) \begin{align}\left(B_{k} \vec{u}\right)(\vec{x}):=\int\limits_{\Gamma}\gamma_{N,\vec{x}}^{\Omega}\left(G_k(\vec{x},\vec{y})\vec{u}(\vec{y})\right)\!ds_{\vec{y}},\vec{x}\in\Gamma,\end{align}

(2.10) \begin{align}\left(A_{k} \vec{u}\right)(\vec{x})& :=\int\limits_{\Gamma}\gamma_{D,\vec{x}}^{\Omega}\left(G_k\left(\vec{x},\vec{y}\right)\vec{u}(\vec{y})\right)\!ds_{\vec{y}}\nonumber\\[3pt] &\quad -\frac{1}{k^2}\nabla_{\vec{x}}\int\limits_{\Gamma}G_k(\vec{x},\vec{y})\nabla_{\vec{y}}\cdot\gamma_{N,\vec{y}}^{\Omega}\vec{u}(\vec{y})ds_{\vec{y}},\vec{x}\in\Gamma,\end{align}

For any two vector functions, $\vec{u}$ and $\vec{v}$ defined on and tangential to $\Gamma$ almost everywhere on $\Gamma$ , we introduce the scalar product [7]:

(2.11) \begin{equation}\langle \vec{u}, \vec{v} \rangle =\int\limits_{\Gamma}\vec{u}(\vec{x})\cdot\overline{\vec{v}}(\vec{x}) ds,\end{equation}

where $\overline{\cdot}$ means complex conjugation.

Suppose there are two function spaces [7]: $H_{\parallel}^{-\frac{1}{2}}\left(\mathrm{div}_{\Gamma},\Gamma\right)$ – the space of the Neumann trace image. It contains vector-valued functions tangential to $\Gamma$ and having continuous components normal to the edges of $\Gamma$ ; $H_{\perp}^{-\frac{1}{2}}\left(\mathrm{curl}_{\Gamma},\Gamma\right)$ – the space of the Dirichlet trace image. It also contains vector-valued functions tangential to $\Gamma$ , but in this case these functions have continuous components tangential to the edges of $\Gamma$ ; The detailed description of these spaces is beyond the scope of this work. For more information, see [7, 17] and [10].

Applying the Dirichlet trace operator to (2.6) and reformulating the resulting expression in terms of the integral operators (2.10) and (2.9), one can derive the following variational problem from [7, p. 51] using special properties of the integral operators listed in [7, p. 48]:

(2.12) \begin{equation}\langle A_k \gamma_N^{\Omega} \vec{E}, \vec{u}\rangle =\langle \gamma_D^{\Omega}\vec{E}, \left(\frac{1}{2}\mathrm{Id} + B_{\overline{k}}\right)\vec{u}\rangle, \forall \vec{u} \in H_{\parallel}^{-\frac{1}{2}}\end{equation}

(2.13) \begin{equation}\gamma_N^{\Omega}\vec{E}\in H_{\parallel}^{-\frac{1}{2}}\left(\mathrm{div}_{\Gamma},\Gamma\right),\gamma_D^{\Omega}\vec{E}\in H_{\perp}^{-\frac{1}{2}}\left(\mathrm{curl}_{\Gamma},\Gamma\right),\end{equation}

where Id is the identity operator.

The integrals involved in the computation of $A_k$ can be easily represented through the scalar single layer potential as it was demonstrated in [7] and [10]. The subsequent is concerned with the operator $B_k$ .

2.2 Modified double layer potential

To obtain a discrete version of (2.12), one has to define two finite subspaces: $W_{\parallel}^{K}\subset H_{\parallel}^{-\frac{1}{2}}\left(\mathrm{div}_{\Gamma},\Gamma\right)$ and $W_{\perp}^{M}\subset H_{\perp}^{-\frac{1}{2}}\left(\mathrm{curl}_{\Gamma},\Gamma\right)$ where K and M are the dimensions of the corresponding finite subspaces. Suppose we have to find a value for the following expression arising from (2.12):

\begin{equation*}\langle \vec{v}, B_{\overline{k}} \vec{u} \rangle,\end{equation*}

\begin{equation*}\vec{v}\in W_{\perp}^{K} \subset H_{\perp}^{-\frac{1}{2}}\left(\mathrm{curl}_{\Gamma},\Gamma\right), \vec{u}\in W_{\parallel}^{M} \subset H_{\parallel}^{-\frac{1}{2}}\left(\mathrm{div}_{\Gamma},\Gamma\right)\end{equation*}

Applying well known identities of vector calculus to the Neumann trace operator, we get similarly to [17, p. 42]:

\begin{align*}\langle \vec{v}, B_{\overline{k}} \vec{u} \rangle&=\left \langle \vec{v},\int\limits_{\Gamma}\gamma_{N,x}^{\Omega}\left(G_{\overline{k}}\left(\vec{x},\vec{y}\right)\vec{u}(\vec{y})\right)\!ds_{\vec{y}}\right\rangle \\[2pt]& =\left \langle \vec{v},\int\limits_{\Gamma}\left(\nabla_{\vec{x}}\times\left(G_{\overline{k}}\left(\vec{x},\vec{y}\right)\vec{u}(\vec{y})\right)\right)\times\vec{n}(\vec{x})ds_{\vec{y}}\right\rangle \\[2pt]& =\left \langle \vec{v},\int\limits_{\Gamma}\left(\nabla_{\vec{x}} G_{\overline{k}}\left(\vec{x},\vec{y}\right)\times\vec{u}(\vec{y})\right)\times\vec{n}(\vec{x})ds_{\vec{y}}\right\rangle \\[2pt]&=\left \langle \vec{v},\int\limits_{\Gamma}\left(-\nabla_{\vec{x}} G_{\overline{k}}\left(\vec{x},\vec{y}\right)\left(\vec{n}(\vec{x})\cdot\vec{u}(\vec{y})\right)+\vec{u}(\vec{y})\frac{\partial_{\vec{x}}G_{\overline{k}}(\vec{x},\vec{y})}{\partial \vec{n}(\vec{x})}\right)\!ds_{\vec{y}}\right\rangle. \end{align*}

For what follows next, the angular brackets are replaced with an integral in accordance with (2.11) and the order of integration is changed as it was done in [17]. This, as it will be shown in the subsequent, does not change the result and is a valid operation of mathematical analysis. However, there might be some numerical issues related to the change of integration order that will be highlighted afterwards.

(2.14) \begin{align} \langle \vec{v}, B_{\overline{k}} \vec{u} \rangle &=-\int\limits_{\Gamma}\left[\int\limits_{\Gamma}\vec{n}(\vec{x})\left(\vec{v}(\vec{x})\cdot\nabla_{\vec{x}}G_k(\vec{x},\vec{y})\right)\!ds_{\vec{x}}\right]\cdot\overline{\vec{u}}(\vec{y})ds_{\vec{y}} \nonumber\\[2pt] &\quad +\int\limits_{\Gamma}\left[\int\limits_{\Gamma}\vec{v}(\vec{x})\frac{\partial_{\vec{x}} G_k\left(\vec{x}, \vec{y}\right)}{\partial \vec{n}(\vec{x})}ds_{\vec{x}}\right]\cdot\overline{\vec{u}}(\vec{y})ds_{\vec{y}}.\end{align}

The integral in square brackets of the second term in (2.14) can be expressed via the double layer potential $K_{\Gamma}^{k}$ [20]:

(2.15) \begin{equation}\left(K_{\Gamma}^k f\right)(\vec{y}) :=\int\limits_{\Gamma}\frac{\partial_{\vec{x}} G_k\left(\vec{x}, \vec{y}\right)}{\partial \vec{n}(\vec{x})}f(\vec{x}) ds_{\vec{x}},\end{equation}

in a component-wise fashion for vector-valued functions.

For more information about the basic properties of $K_{\Gamma}^k$ that allow the mentioned above interchange of integrals in (2.14), see [20].

The other integral in square brackets can also be expressed in terms of MDLP $M_{\Gamma}^{k}$ introduced in [17] for each component of $\vec{n}$ :

(2.16) \begin{equation}\left(M_{\Gamma}^k \vec{p}\right)(\vec{y}):=\int\limits_{\Gamma} \nabla_{\vec{x}} G_k(\vec{x},\vec{y}) \cdot \vec{p}(\vec{x})ds_{\vec{x}}.\end{equation}

The component-wise approach mentioned here will be demonstrated in Section 4.

3 Analysis of MDLP

Here, we demonstrate the singularity of MDLP. To do so, we present the formula of integration by parts on a surface in the first subsection. In the second one, we use this formula to prove the mentioned property of MDLP.

3.1 Example

The singularity of $M_S^k$ can be demonstrated with help of a simple example. Let S be a square in OXY plane. The square belongs to the first quadrant on the OXY plane and the edge length is equal to 1. Also, let $k = 0$ and $\vec{u} = \left(1, 0, 0\right)$ in Cartesian coordinates. When $k = 0$ , to check the singularity of $M_S^k$ is a trivial task of calculus – this is the only place in the article where we allow this case to be true. It’s relatively easy to see that the same singularity occurs when $k \neq 0$ . Consider the following relation for arbitrary k:

(3.1) \begin{align} \left(M_{S}^k \vec{v}\right)(\vec{y}) & = \int\limits_{S} \nabla_{\vec{x}} G_k(\vec{x},\vec{y}) \cdot \vec{v}(\vec{x})ds_{\vec{x}} \nonumber\\[3pt] & =\int\limits_{S}\left(-\frac{e^{-k\|\vec{x} - \vec{y}\|}\left(\vec{x} - \vec{y}\right)}{4\pi\|\vec{x} - \vec{y}\|^3}- k \frac{e^{-k\|\vec{x} - \vec{y}\|}\left(\vec{x} - \vec{y}\right)}{4\pi\|\vec{x} - \vec{y}\|^2}\right)\cdot \vec{v}(\vec{x})ds_{\vec{x}} \nonumber\\[3pt]&= \int\limits_{S} e^{-k\|\vec{x} - \vec{y}\|}\left(-\frac{1 + k \|\vec{x} - \vec{y}\|}{4\pi\|\vec{x} - \vec{y}\|^3}\right)\left(\vec{x} - \vec{y}\right)\cdot \vec{v}(\vec{x})ds_{\vec{x}} \nonumber\\[3pt] &=\int\limits_{S} e^{-k\|\vec{x} - \vec{y}\|}\left(\frac{ e^{k\|\vec{x} - \vec{y}\|} - 1 - k \|\vec{x} - \vec{y}\|}{4\pi\|\vec{x} - \vec{y}\|^3}\right)\left(\vec{x} - \vec{y}\right)\cdot \vec{v}(\vec{x})ds_{\vec{x}} + \left(M_{\Gamma}^0 \vec{v}\right)(\vec{y}) \nonumber\\[3pt] & = \int\limits_{S} Q\left(\vec{x}, \vec{y}\right)\!ds_{\vec{x}} + \left(M_{\Gamma}^0 \vec{v}\right)(\vec{y}).\end{align}

One can see, that the function Q of $\vec{x}$ and $\vec{y}$ integrated in the first term of (3.1) is bounded when $\vec{x}\to\vec{y}$ . Hence, the first-term integral is also bounded as a function of $\vec{y}$ .

In Figure 2, we plot the function $\left(M_{S}^0 \vec{u}\right)(\vec{y})$ at the points taken along the line parallel to the direction of $\vec{u}$ (along OX) and splitting the square in half to demonstrate the growth at the edges. The line of the output and the square both can be seen in Figure 1.

Figure 1. Illustration of the example on the square of unit edge.

Figure 2. Demonstration of singularity for a particular case.

Here, we must mention that the integrand of the MDLP in (2.14) vanishes on a plane surface because $\vec{u}$ and $\vec{n}$ become orthogonal. The analysis below is valuable for the proper consideration of the edges and corners of surface geometry.

3.2 Integration by parts on a surface

Let S be a smooth bounded surface in $\mathbb{R}^3$ . Here and in the following we suppose that $S = \overline{S} = S \cup \partial S$ . Suppose that there exists a one-to-one mapping from $\mathbb{R}^2$ to S with this map’s first and second derivatives being continuous on $\mathbb{R}^2$ . Hence, every point on S is a function of two parameters $q^1$ and $q^2$ , and the two vectors:

\begin{equation*}\vec{r}_p = \frac{\partial \vec{r}}{\partial q^p}, \vec{r}\in S, p=1,2,\end{equation*}

are linearly independent.

The normal direction $\vec{n}$ to the surface S is naturally defined as a normalised cross product of these vectors:

\begin{equation*}\vec{n} = \frac{\vec{r}_1 \times \vec{r}_2}{\left|\vec{r}_1 \times \vec{r}_2\right|}.\end{equation*}

The triplet of vectors $\vec{r}_1$ , $\vec{r}_2$ and $\vec{n}$ constitutes a reference frame which is always right.

It is customary to introduce the vector $\vec{\tau}$ defined on $\partial S$ and tangential to it [9]:

\begin{equation*}\vec{\tau} = \sum\limits_{s = 1}^{2} \frac{d q^s}{d t} \vec{r}_s,\end{equation*}

where t is the parameter of length measured along the boundary curve. Finally, let $\vec{p}$ be defined in accordance with the formula:

(3.2) \begin{equation}\vec{p} = \vec{\tau} \times \vec{n},\end{equation}

notice that $\vec{p}$ is directed to the exterior of S.

Theorem 3.1. Let $\vec{a}$ be a continuously differentiable vector function defined on a smooth surface S and be tangential to it everywhere on S. If the surface S can be mapped to a subdomain Q of $\mathbb{R}^2$ such that the map is continuously differentiable then the equality below is true:

(3.3) \begin{equation}\int\limits_{S} \nabla\cdot\vec{a} ds = \int\limits_{\partial S} \vec{p} \cdot \vec{a} dt,\end{equation}

where t is the parameter of length of the boundary curve $\partial S$ .

Proof. See [3] or [4].

Theorem 3.2. Let u be the scalar function differentiable on S and $\vec{a}$ be just as before. Then

(3.4) \begin{equation}\int\limits_{S}\nabla u \cdot \vec{a} ds =-\int\limits_{S}u \nabla \cdot \vec{a} ds + \int\limits_{\partial S} u\vec{a} \cdot \vec{p} dt.\end{equation}

The relation (3.4) is known as the integration by parts formula.

Proof. This result directly follows from Theorem 3.1 and from the fact that

\begin{align*}\nabla\cdot\left(u \vec{a}\right) =\nabla u \cdot \vec{a} + u \nabla \cdot \vec{a}. \\[-30pt] \end{align*}

3.3 Singularity of the MDLP

Here we prove the fact, which is illustrated in Subsection 3.1, for a more general case. Let the surface S be as in the previous subsections. The integral in (2.16) is regular if the point-argument $\vec{y}$ is not on S. Otherwise, the integral in (2.16) is computed in the following sense: one must take a ball centered at $\vec{y}$ with radius r and subtract it from S. For the resulting integral, the limit should be computed when $r\to 0$ . Let’s denote this ball as $B_r(\vec{y})$ , then using formula (3.4) we get:

(3.5) \begin{align} \left(M_{S\backslash B_r(\vec{y})}^k\vec{u}\right)(\vec{y}) & = \int\limits_{S \backslash B_r(\vec{y})} \nabla_{\vec{x}} G_k(\vec{x},\vec{y}) \cdot \vec{u}(\vec{x})ds_{\vec{x}} \nonumber\\[3pt]&= - \int\limits_{S \backslash B_r(\vec{y})} G_k(\vec{x},\vec{y}) \big(\nabla_{\vec{x}} \cdot \vec{u}(\vec{x}) \big) ds_{\vec{x}} \nonumber\\ &\quad + \oint\limits_{ \partial \big( S \backslash B_r(\vec{y}) \big)}G_k(\vec{x}\left(t\right),\vec{y})\big(\vec{u}(\vec{x}\left(t\right)) \cdot \vec{p}(\vec{x}\left(t\right)) \big)dt, \end{align}

where $\vec{p}$ is the unit vector orthogonal to the normal vector $\vec{n}$ at the points of $\partial(S \backslash B_r(\vec{y}))$ and it is directed to the exterior of $S \backslash B_r(\vec{y})$ . The mentioned means that

(3.6) \begin{equation}\left(M_{S}^k\vec{u}\right)(\vec{y}) := \lim_{r\to 0} \left(M_{S\backslash B_r(\vec{y})}^k\vec{u}\right)(\vec{y}).\end{equation}

The use of formula (3.4) is justified because the function $\vec{u}$ can be expressed in terms of $\vec{r}_s, s = 1,2$ only, so the normal component of $\nabla_{\vec{x}}G_k(\vec{x},\vec{y})$ is multiplied by zero.

Theorem 3.3. Let S be a simple analytic smooth open surface. Let $\vec{y}\in\partial S$ and the boundary curve $\partial S$ be such that there exists a positive number $\varepsilon$ for which all the spheres centered at $\vec{y}$ with their radiuses being smaller than $\varepsilon$ have only two intersection points with $\partial S$ . Also, for the part of $\partial S$ inside the ball $B_{\varepsilon}(\vec{y})$ the diffeomorphism exists between the length of the curve measured from $\vec{y}$ to any point on $\partial S \cap B_{\varepsilon}(\vec{y})$ and the Euclidean distance measured between the same points. The diffeomorphism must be a monotonous function, that is: when the Euclidean distance decreases, the distance along the curve decreases as well. Under these conditions, if $\vec{u}$ is tangential to and differentiable on S and if the function $\vec{u}(\vec{x}) \cdot \vec{p}(\vec{x})$ is continuous when $\vec{x}\in \partial S \cap B_{\varepsilon}(\vec{y})$ and $\vec{u}(\vec{x}) \cdot \vec{p}(\vec{x})$ does not change its sign in $\partial S \cap B_{\varepsilon}(\vec{y})$ then the absolute value of the modified double layer potential computed at the point $\vec{y}$ is infinitely large.

Proof. The first term on the right hand side of formula (3.5) has a form of the single layer potential operator applied to the function $\nabla_{\vec{x}}\cdot \vec{u}(\vec{x})$ . The single layer potential is expressed as follows [20]:

(3.7) \begin{equation}\left(\tilde{V}_{S}^{k} \omega\right)\left(\vec{y}\right) =\int\limits_{S} G_k(\vec{x},\vec{y}) \omega (\vec{x}) ds_{\vec{x}},\end{equation}

where $\omega$ is a scalar function. The single layer operator is known to be weakly singular and bounded [20] when $\vec{y}\in\partial S$ and hence the limit:

(3.8) \begin{align} & \lim_{r\to 0}\int\limits_{\vec{x}\in S \backslash B_r(\vec{y})} G_k(\vec{x},\vec{y}) \big(\nabla_{\vec{x}} \cdot \vec{u}(\vec{x}) \big) ds_{\vec{x}} \nonumber\\[3pt] &\quad =\lim_{r\to 0}\left(\tilde{V}_{S\backslash B_r(\vec{y})}^k\nabla\cdot\vec{u}\right)(\vec{y}) = \left(\tilde{V}_{S}^k\nabla\cdot\vec{u}\right)(\vec{y})\end{align}

exists and is finite. For the second term on the right hand side of formula (3.5), let’s choose an arbitrary positive small number $\varepsilon$ . Using the addition property we have

(3.9) \begin{align} & \oint\limits_{ \vec{x} \in \partial \big(S \backslash B_r(\vec{y}) \big)}G_k(\vec{x},\vec{y})\big(\vec{u}(\vec{x}) \cdot \vec{p}(\vec{x}) \big)dt \nonumber \\[3pt] &= \int\limits_{ \vec{x} \in \partial\left( S \backslash B_r(\vec{y}) \right):\left|\vec{x} - \vec{y} \right| \geq \varepsilon} \bullet + \int\limits_{ \vec{x} \in \partial \left( S \backslash B_r(\vec{y}) \right):\left|\vec{x} - \vec{y} \right| < \varepsilon} \bullet.\end{align}

The second term on the right hand side of (3.9) is of a particular interest as the first term is just a regular integral that definitely has a finite value.

Let $\tau$ be the parameter of length from the point $\vec{x}$ to the point $\vec{y}$ and t be the distance between these points measured along the curve. By the conditions of the theorem, there exists a small positive real number $\varepsilon$ such that the function $\partial t \over \partial \tau$ does not change its sign as well as the function $\vec{u}\cdot\vec{p}$ . Also, the surface of $B_\varepsilon(\vec{y})$ has a radius $\varepsilon$ small enough so there are only two points of intersection between the curve $\partial S$ and the surface of $B_\varepsilon(\vec{y})$ . The same is true for the spherical surface of $B_r(\vec{y})$ because $r \to 0$ and for the fixed value of $\varepsilon$ we can consider the case of $r < \varepsilon$ . We denote the points of integration along the two different parts of the curve as $\vec{x}_1$ and $\vec{x}_2$ as a reminder that there are two points of intersection between the curve and $\partial B_\varepsilon(\vec{y})$ .

Let’s denote the points of intersection between the spherical surface of $B_r(\vec{y})$ and the curve $\partial S$ as $\vec{\rho}_1$ and $\vec{\rho}_2$ . The curve lying on S and connecting the points $\vec{\rho}_1$ and $\vec{\rho}_2$ along a curve-linear arch of the boundary $\partial \left( S \backslash B_r(\vec{y}) \right)$ we denote as $C_{arch}$ . The distance from the points of $C_{arch}$ to the point $\vec{y}$ remains equal to r.

$C_{arch}$ belongs to the surface of $B_r(\vec{y})$ . Hence, we can project this curve onto a ball of unit radius. We denote the projection curve as $C_{proj}$ . Let $\alpha$ be the parameter of length along the projection curve $C_{proj}$ . Then ${{d t} \over { d \alpha} }= r$ . Hence, we get

(3.10) \begin{align} & \lim_{r\to 0}\int\limits_{ \vec{x}\left(t\right) \in \partial \left( S \backslash B_r(\vec{y}) \right):\left|\vec{x} - \vec{y} \right| < \varepsilon}G_k(\vec{x}\left(t\right),\vec{y})\big(\vec{u}(\vec{x}\left(t\right)) \cdot \vec{p}(\vec{x}\left(t\right)) \big)dt \nonumber \\&= \lim_{r \to 0}\left(\frac{e^{-k r}}{4 \pi}\int\limits_{C_{proj}}\left(\vec{u}(\vec{x}(t(\alpha))) \cdot \vec{p}(\vec{x}(t(\alpha)))\right) d\alpha\right. \nonumber\\ &\quad + \left. \sum\limits_{j = 1}^{2}\int\limits_{r}^{\varepsilon}\frac{e^{-k\tau}}{4 \pi\tau}\left(\vec{u}(\vec{x}_j(\tau)) \cdot \vec{p}(\vec{x}_j(\tau))\right)\frac{\partial t}{\partial \tau}d\tau\right).\end{align}

The functions $\partial t \over \partial \tau$ , $\vec{u}$ and $\vec{p}$ are continuous at the point $\vec{y}$ and are not equal to zero. Also, $\partial t \over \partial \tau$ is bounded because the principal curvatures of S are bounded inside $B_\varepsilon(\vec{y})$ . As one can see, the limit exists for the first term in (3.10) because this integral is regular at $r = 0$ . The two remaining terms have the same sign by the conditions of the theorem. We arrive to the conclusion now:

(3.11) \begin{align} C_j & = \min_{\tau \in (0, \varepsilon)}\left|\frac{\partial t}{\partial \tau}e^{-k\tau}\left(\vec{u}(\vec{x}_j(\tau)) \cdot \vec{p}(\vec{x}_j(\tau))\right)\right| > 0, j=1,2, \nonumber\\C &=\min\!\left\{C_1,C_2\right\}\neq 0, \nonumber\\ & \left|\lim_{r\to 0}\sum\limits_{j=1}^{2}\int\limits_{r}^{\varepsilon}\frac{\partial t}{\partial \tau}\frac{e^{-k\tau}}{4 \pi\tau}\left(\vec{u}(\vec{x}_j(\tau)) \cdot \vec{p}(\vec{x}_j(\tau))\right)d\tau\right| > C\left|\int\limits_{0}^{\varepsilon}\frac{1}{\tau}d\tau\right| = \infty\\[-25pt] \nonumber\end{align}

Corollary 1. By choosing values of constants $\tilde{C}_1$ and $\tilde{C}_2$ one can find the upper bound of the form:

\begin{equation*}\tilde{C}_j = \max_{\tau \in (0, \varepsilon)}\left|\frac{\partial t}{\partial \tau}e^{-k\tau}\left(\vec{u}(\vec{x}_j(\tau)) \cdot \vec{p}(\vec{x}_j(\tau))\right)\right| > 0, j=1,2,\end{equation*}

\begin{equation*}\tilde{C} =2\max\!\left\{\tilde{C}_1,\tilde{C}_2\right\},\end{equation*}

(3.12) \begin{equation}\left|\sum\limits_{j=1}^{2}\int\limits_{r}^{\varepsilon}\frac{\partial t}{\partial \tau}\frac{e^{-k\tau}}{4 \pi\tau}\left(\vec{u}(\vec{x}_j(\tau)) \cdot \vec{p}(\vec{x}_j(\tau))\right)d\tau\right| < \tilde{C}\left|\int\limits_{r}^{\varepsilon}\frac{1}{\tau}d\tau\right|.\end{equation}

From this, we conclude that the MDLP has a logarithmic growth and therefore it’s integrable in a singular sense.

3.4 Weakly singular integral of MDLP

Let $\Omega$ be a subdomain in $\mathbb{R}^3$ . Also, let its boundary $\Gamma = \partial \Omega$ be composed out of the finite number of surface pieces that are parametrisable in the sense of Theorem 3.1:

(3.13) \begin{equation}\Gamma = \bigcup_{i = 1}^N S_i,\end{equation}

where N is the number of surface elements constituting $\Gamma$ . We note specifically that the intersection of any two sets $S_i$ and $S_j$ is either $\varnothing$ , a curve or a single point but never a surface.

We are interested in the integral:

(3.14) \begin{equation}\int\limits_{S_i} w(\vec{y})\left(M_{\Gamma}^{k} \vec{u}\right)(\vec{y})ds_{\vec{y}},\end{equation}

where w is a scalar function continuous on $S_i$ . Substituting (3.5) into (3.14) we get:

(3.15) \begin{align} & \quad \int\limits_{S_i} w(\vec{y})\left(M_{\Gamma}^{k} \vec{u}\right)(\vec{y})ds_{\vec{y}} \nonumber \\ & = \int\limits_{S_i} w(\vec{y})\left(-\lim_{r\to 0}\int\limits_{S \backslash B_r(\vec{y})} G_k(\vec{x},\vec{y}) \big(\nabla_{\vec{x}} \cdot \vec{u}(\vec{x}) \big) ds_{\vec{x}}\right)\!ds_{\vec{y}} + J,\end{align}

where

(3.16) \begin{equation}J = \int\limits_{S_i} w(\vec{y})\left(\lim_{r\to 0}\oint\limits_{ \partial \big( S_j \backslash B_r(\vec{y}) \big)}G_k(\vec{x}\left(t\right),\vec{y})\big(\vec{u}(\vec{x})\left(t\right) \cdot \vec{p}(\vec{x}\left(t\right)) \big)dt\right)\!ds_{\vec{y}}.\end{equation}

By Theorem 3.3, the integration on $S_i$ (3.16) is singular. Hence, there are issues with numerical integration of J because the most popular approaches to the numerical integration consider the integrated function to be at least continuous. A better way of computing J would be to interchange the integrals in (3.16), so that the integration problem can be reduced to the known ones. To do so, one can take a vicinity of $\partial S_j$ that we denote as $B_{\varepsilon}\left(\partial S_j\right)$ and $\partial B_{\varepsilon}\left(\partial S_j\right)$ is a curved pipe-like surface situated along the curve $\partial S_j$ with the diameter equal to $\varepsilon$ . Now, we can use the additivity of the integral in (3.16) taken over $S_i$ :

(3.17) \begin{align} J &= \int\limits_{S_i} w(\vec{y})\left(\lim_{r\to 0}\oint\limits_{ \partial \big( S_j \backslash B_r(\vec{y}) \big)}G_k(\vec{x},\vec{y})\big(\vec{u}(\vec{x}) \cdot \vec{p}(\vec{x}) \big)dt_{\vec{x}}\right)\!ds_{\vec{y}} \nonumber\\& =\int\limits_{S_i \backslash B_{\varepsilon}\left(\partial S_j\right)} w(\vec{y})\lim_{r\to 0}\oint\limits_{ \partial \big( S_j \backslash B_r(\vec{y}) \big)}F\left(\vec{x}\left(\vec{t}\right),\vec{y}\right)\! dtds_{\vec{y}} \nonumber\\&\quad +\int\limits_{S_i \bigcap B_{\varepsilon}\left(\partial S_j\right)} w(\vec{y})\lim_{r\to 0}\oint\limits_{ \partial \big( S_j \backslash B_r(\vec{y}) \big)}F\left(\vec{x}\left(\vec{t}\right),\vec{y}\right)\! dtds_{\vec{y}}, \nonumber \\ F(\vec{x},\vec{y}) &=G_k(\vec{x},\vec{y})\big(\vec{u}(\vec{x}) \cdot \vec{p}(\vec{x}) \big).\end{align}

The first term of summation in (3.17) has only regular integrals and hence the two integrals are interchangeable. The second term in (3.17) has a weakly singular external integral (computed for $S_i \bigcap B_{\varepsilon}\!\left(\partial S_j\right)$ ) and by the corollary (1) it converges to zero with $\varepsilon\to 0$ . Using (3.7), we get

(3.18) \begin{equation}J =\oint\limits_{ \partial S_j}\left(\vec{u}(\vec{x}\left(t\right)) \cdot \vec{p}(\vec{x}\left(t\right))\right)\left(\tilde{V}_{S_i}^{k} w\right)\left(\vec{x}\left(t\right)\right)\!dt.\end{equation}

4 Method discussion and implementation notes

The integration of MDLP, as we demonstrated in this paper, can be reduced to the integration of the single layer potential along a curve and on a surface. The techniques of the singular integration developed for numerical implementation are discussed in [14, 15, 8, 12] and [20]. The numerical integration of continuous functions is also needed. Such methods can be found in [16].

We can substitute (3.18) into (2.14), we get

(4.1) \begin{align} \langle \vec{v}, B_{\overline{k}} \vec{u} \rangle & =\int\limits_{\Gamma}\left[\int\limits_{\Gamma}\vec{v}(\vec{x})\left(\vec{n}(\vec{x})\cdot\nabla_{\vec{x}}G_k(\vec{x},\vec{y})\right)\!ds_{\vec{x}}\right]\cdot\overline{\vec{u}}(\vec{y})ds_{\vec{y}} \nonumber \\&\quad -\int\limits_{\Gamma}\left[\int\limits_{\Gamma}\vec{n}(\vec{x})\left(\vec{v}(\vec{x})\cdot\nabla_{\vec{x}}G_k(\vec{x},\vec{y})\right)\!ds_{\vec{x}}\right]\cdot\overline{\vec{u}}(\vec{y})ds_{\vec{y}} \nonumber \\ & =\int\limits_{\Gamma}\left[\int\limits_{\Gamma}\vec{v}(\vec{x})\left(\vec{n}(\vec{x})\cdot\nabla_{\vec{x}}G_k(\vec{x},\vec{y})\right)\!ds_{\vec{x}}\right]\cdot\overline{\vec{u}}(\vec{y})ds_{\vec{y}} \nonumber \\&\quad -\int\limits_{\Gamma}\sum\limits_{q=1}^{3}\left(M_{\Gamma}^{k}\left(\vec{n}\cdot\vec{c}_q\right)\vec{v}\right)(\vec{y})\left(\vec{c}_q\cdot\overline{\vec{u}}(\vec{y})\right)\!ds_{\vec{y}},\end{align}

where $\vec{c}_q$ are the Cartesian basis unit vectors. Applying the results of the previous chapters, we find

\begin{align*}& -\int\limits_{\Gamma}\sum\limits_{q=1}^{3}\left(M_{\Gamma}^{k}\left(\vec{n}\cdot\vec{c}_q\right)\vec{v}\right)(\vec{y})\left(\vec{c}_q\cdot\overline{\vec{u}}(\vec{y})\right)\!ds_{\vec{y}} \end{align*}

(4.2) \begin{align}\nonumber \\[-35pt]& = \int\limits_{\Gamma}\sum\limits_{q=1}^{3}\left(\tilde{V}_{\Gamma}^{k} \nabla \cdot\left(\left(\vec{n}\cdot\vec{c}_q\right)\vec{v}\right)\right)(\vec{y})\left(\vec{c}_q\cdot\overline{\vec{u}}(\vec{y})\right)\!ds_{\vec{y}} \nonumber\\[3pt] &\quad - \sum\limits_{i = 1}^N \oint\limits_{ \partial S_i}\sum\limits_{q = 0}^{3}\left(\left(\vec{n}(\vec{x})\cdot\vec{c}_q\right)\left(\vec{v}(\vec{x}) \cdot \vec{p}(\vec{x})\right)\right)\left(\tilde{V}_{\Gamma}^{k} \left(\vec{c}_q\cdot\overline{\vec{u}}\right)\right)(\vec{x})dt_{\vec{x}}.\end{align}

Formula (4.2) deals with integration of continuous functions on the integration domain. Hence, these integrals can be computed numerically with help of integration scheme applicable to the case of continuous functions.

5 About coupling with finite element method

To use the coupling with the finite element method (FEM), we exploited the concept of the Steklov–Poincare operator variational formulation presented in [7].

In order to define the Steklov–Poincare integral operator, the operators $N_k$ and $C_k$ were defined in [7, p. 47, formula (5.10)]:

(5.1) \begin{equation}\left(N_k \vec{u}\right)(\vec{x}) :=\gamma_{N,\vec{x}} \nabla_{\vec{x}} \times\int\limits_{\Gamma}G_k(\vec{x},\vec{y})\left(\vec{u}(\vec{y})\times\vec{n}(\vec{y})\right)\!ds_{\vec{y}},\vec{x}\in\Gamma,\end{equation}

(5.2) \begin{equation}\left(C_k \vec{u}\right)(\vec{x}) :=\gamma_{D,\vec{x}} \nabla_{\vec{x}} \times\int\limits_{\Gamma}G_k(\vec{x},\vec{y})\left(\vec{u}(\vec{y})\times\vec{n}(\vec{y})\right)\!ds_{\vec{y}} -\frac{1}{2}\vec{u}(\vec{x}),\vec{x}\in\Gamma,\end{equation}

If the function $\vec{u}$ is the solution of (2.1) then

(5.3) \begin{equation}\left(S_k \gamma_D^{\Omega} \vec{u}\right)(\vec{x}) =\gamma_N \vec{u},\end{equation}

where

(5.4) \begin{equation}S_k = N_k + \left(\frac{1}{2}\mathrm{Id} + B_k\right)A_k^{-1}\left(\frac{1}{2}\mathrm{Id} - C_k\right).\end{equation}

The existence of the inverse operator $A_k^{-1}$ is proved in [7]. The corresponding variational problem looks like this:

(5.5) \begin{equation}\left \langle S_k \gamma_D^{\Omega} \vec{u},\vec{v}\right\rangle=\left \langle \gamma_N \vec{u},\vec{v}\right\rangle,\forall \vec{v}\in H_{\perp}^{-\frac{1}{2}}\!\left(\mathrm{curl}_{\Gamma},\Gamma\right),\end{equation}

where $\vec{v}$ is the test function of the variational problem.

The coupling of the variational problem (5.5) with the finite element method for eddy current problems is not a difficult task and it was presented in a number of works such as [18] and [22].

The operators $N_k$ and $C_k$ do not involve computation of integrals that cannot be seen as special cases of the double layer potential and the single layer potential operators both computed either for the scalar functions or vector functions in a component-wise fashion. To find out more about these operators, see [7]. For more information about BEM-ECP, see [11].

6 Model problem

A model problem was presented in [17] to test convergence for the numerical approach implemented in the context of BEM-ECP. Here, we solve the same problem to test our approach.

Let $\Omega$ be a conducting ball of radius a with an infinitely thin coil of radius $b > a$ both centered at the same point as it is illustrated by Figure 3.

Figure 3. A ball with a coil centered at the same point.

The integral value of exciting current I is not zero at the points of the coil whereas the direction of current is continuously tangential to the coil. The spherical coordinate system is defined by three parameters: r, $\varphi$ , $\theta$ where r is the distance from the center of the ball, $\varphi$ is the angle of rotation around the axis perpendicular to the plane of the coil and $\theta$ is the angle between the mentioned axis and the radius vector $\vec{x}$ of the point where the value of $\vec{E}$ is desired. The equality $\theta = {\pi\over2}$ corresponds to the position of the coil. The analytic solution $\vec{E}$ has only one non-zero component $E_{\varphi}\left(r, \theta\right)$ in the spherical coordinate system. Let $E_{\varphi}^{I}$ be the desired component of the electric field $\vec{E}$ inside the ball and $E_{\varphi}^{O}$ – outside the ball. The solution can be expressed via the formulas [17, p. 87, formulas (6.28)–(6.34)]:

(6.1) \begin{equation}E_{\varphi}^{I} \! \left(r, \theta\right) =i \omega\sum\limits_{n = 0}^{\infty}E_n^{I}\cdot r^{-\frac{1}{2}} \cdot J_{2n+{3 \over 2}}\!\left(i r k\right) \cdot P_{2n + 1}^{1} (u),\end{equation}

(6.2) \begin{equation}E_{\varphi}^{O} \! \left(r, \theta\right) =i \omega\sum\limits_{n = 0}^{\infty}\left(E_n^{O} \cdot r^{-2n - 2} P_{2n + 1}^{1} (u) - E_{n}^{C} \cdot r^{2n+1} P_{2n + 1}^{1} (u)\right),\end{equation}

with

(6.3) \begin{equation}E_{n}^{O} = E_{n}^{C} a^{4n + 3}\frac{i a k \frac{J_{2n + {1 \over 2}}\left(i a k\right)}{J_{2n + {3 \over 2}}\left(i a k\right)} - 2n - 1 - \mu_{r} \left(2n + 2\right) }{ i a k \frac{J_{2n + {1 \over 2}}\left(i a k\right)}{J_{2n + {3 \over 2}}\left(i a k\right)} - 2n - 1 + \mu_r \left(2n + 1\right) },\end{equation}

(6.4) \begin{equation}E_{n}^{I} = \frac{E_{n}^{O} \cdot a^{-2n-{3 \over 2}} - E_{n}^{C} \cdot a^{2n + {3 \over 2}}}{J_{2n+{3\over2}}\left(i a k\right)}\end{equation}

(6.5) \begin{equation}E_{n}^{C} = \frac{I b \mu_0 \left(-1\right)^{n}\left(2n - 1\right)!!}{2^{n+2} \left(n + 1\right)! b^{2n + 2}},\end{equation}

(6.6) \begin{equation}u = \cos\!\left(\theta\right), \mu = \mu_r \mu_0,\end{equation}

where $\mu_0$ is the magnetic permeability of free space, $\mu$ is the magnetic permeability of the conducting ball, $P_{n}^{m}$ are the associated Legendre functions of the first kind, $J_\nu$ are the Bessel functions of the first kind.

The problem we solved contains two subdomains: the inside of the conducting ball and the outside for which we specified a conductivity value $\sigma_0 = 0.005 \left(\Omega {\rm m}\right)^{-1}$ made small enough to approximate the case of zero conductivity so that formula (2.8) would still make sense. The parameters of the model problem under consideration are: $I = 1000$ A, $\omega = 2\pi 10$ kHz, $a = 0.05$ m and $b = 0.065$ m. The parameters of the conducting ball are: $\sigma = 0.8 \cdot 10^6 \left(\Omega {\rm m}\right)^{-1}$ , $\mu_r = 10.0$ .

The coarse spherical mesh we used for the computation contains 234 nodes and 464 elements. Its finer counterpart we used to test the approximation contains 946 nodes and 1888 elements. The finest mesh we tested contains 3810 nodes and 7616 elements. These meshes are not subdivisions of one another – they were constructed via making two times more steps in a parameter space of the sphere with each refinement. The coarse, finer and finest meshes are illustrated in Figure 4.

Figure 4. coarse, finer and finest meshes.

We used two types of numerical experiments. In the first one, the matrix of the conducting ball was computed via the boundary element method. In the second one, the coupling with FEM was exploited.

For the approach involving FEM, we used volume meshes with radial condensation of volume element size directed to the spherical surface. The number of radial subdivisions $N_r$ applied to the volume mesh is equal to 4 for the coarse case, 8 for the finer case and 16 for the finest with the parameter of condensation q equal to 3, $\sqrt{3}$ , and $\sqrt[4]{3}$ respectively. The radius of the ball (denoted above as a) is equal to the following expansion:

(6.7) \begin{equation}a = C_r\sum\limits_{p = 0}^{p = N_r - 1} q^{p},\end{equation}

where the smallest (the first) term in the sum corresponds to the first step from the spherical surface made towards the ball centre (for $p = 0$ ); $C_r$ – the smallest step made from the spherical surface. The step can be deduced in accordance with (6.8).

(6.8) \begin{equation}C_r = \frac{a(1 - q)}{1 - q^{N_r}},\end{equation}

The illustration of the aforementioned steps is presented in Figure 5.

Figure 5. Illustration of volume subdivision steps.

We compare the smallest step $C_r$ from the spherical surface with the value of penetration depth $\delta$ computed in accordance with (6.9):

(6.9) \begin{equation}\delta = \sqrt{ {2} \over {\sigma \mu \omega} }.\end{equation}

In our case, this value is approximately equal to $1.78$ mm. $C_r$ for the coarse mesh is equal to $1.25$ mm, for the finer mesh – $0.458$ mm and for the finest – $0.198$ mm.

In these experiments, we implemented the incomplete edge-associated basis functions used both for FEM and BEM-ECP parts. Their edge-tangential components remain constant along the points of the edges they are associated with. For more information about the edge-associated vector basis applied to the FEM part, see [13, 5]. As for the case of the BEM-ECP part, see [17] and [7].

We used a coil mesh with a square section that has 2 mm edge length and 40 sections to approximate the round coil geometry. We used the boundary condition (6.10) applied to the thin coil in order to simulate the source of current. It has the following form [7]:

(6.10) \begin{equation}\left.\left(-{1 \over \mu_0 }\nabla \times \vec{E}\right) \times \vec{n}\right|_{\Gamma_c} =\left.i\omega\vec{H} \times \vec{n}\right|_{\Gamma_c},\end{equation}

where $\Gamma_c$ is the coil surface, $\vec{n}$ is the external normal vector defined on $\Gamma_c$ and the function $\vec{H}$ is chosen to be such that the equality:

(6.11) \begin{equation}\oint\limits_{C} \vec{H} \cdot \vec{dl} = I\end{equation}

is true for all the loops C circumventing the coil. On the spherical surface along a meridian, we computed the surface current $I^{s}$ from one pole to the other:

(6.12) \begin{equation}I^{s} = \sigma \left|\vec{E}\right| {\sqrt{2} \over 2}.\end{equation}

The numerical solutions and the analytical solutions are illustrated in Figures 6 and 8. The corresponding values of numerical discrepancies divided by the maximum absolute value of the analytical solution are presented in Figures 7 and 9.

Figure 6. surface current, BEM-ECP comparison.

Figure 7. values of relative discrepancy, BEM-ECP comparison.

Figure 8. surface current, BEM-ECP coupled with FEM.

Figure 9. values of relative discrepancy, BEM-ECP coupled with FEM.

We estimated the discrepancy numerically by taking $N_{\varphi}$ values of $\theta$ such that $\theta_n = n{\pi \over N_{\varphi}}$ where $n = 0..N_{\varphi}$ . We used these values to compute the relative error $\epsilon$ by formula (6.13).

(6.13) \begin{equation}\epsilon :=\sqrt{\frac{\sum\limits_{n = 1}^{N_{\varphi}}\left|I^{s}\left(a,\theta_n\right) - I_{num}^{s}\left(a,\theta_n\right)\right|^{2}}{\max\limits_{n = \overline{1, N_{\varphi}}}\left|I^{s}\!\left(a,\theta_n\right)\right|^2}},\end{equation}

where $I^{s}$ is the approximation computed as a truncated series in accordance with (6.1)–(6.5) and (6.12), $I_{num}^{s}$ is the numerical solution obtained via one of the mentioned approaches. These values are depicted in Table 1 for the mentioned numerical experiments.

Table 1. Values of $\epsilon$ computed for computational methods applied to three different meshes

The number of variables used in each mentioned computational approach is listed in Table 2. The number of variables used for the computation on coil is equal to 480 for all the cases listed in Table 2. It is constant because further subdivisions (with adjustments of the points to fit the coil geometry) of the coil mesh did not influence the results.

Table 2. Number of variables used for each computation

We performed the direct integration of the MDLP implemented in accordance with the formulas mentioned in [17] and we found no significant difference from our results. Hence, we introduce a primitive test with a box to clearly demonstrate that the problem does occur – the geometry is not smooth now.

The mentioned box is expressed as $\Omega = \left[\left[-0.5, 0.5\right], \left[-0.5, 0.5\right], \left[-0.5, 0.5\right]\right]$ (see Figure 10).

Figure 10. Illustration of the box example.

The test solution to (2.1) with $k = i$ is expressed as a function:

(6.14) \begin{align} \vec{u}(\vec{x}) &=0\cdot\vec{c}_1 + 0\cdot\vec{c}_2 + \operatorname{cos}\!\left(\vec{x}\cdot\vec{c}_1\right)\vec{c}_3, \nonumber\\[3pt] \nabla\times\vec{u}(\vec{x}) & =0\cdot\vec{c}_1 + \operatorname{sin}\!\left(\vec{x}\cdot\vec{c}_1\right)\cdot\vec{c}_2 + 0\cdot\vec{c}_3,\end{align}

where $\vec{c}_i$ – Cartesian basis functions.

Suppose the Neumann trace is known. Then, we solve the problem (5.5) to find an approximate value of the Dirichlet trace in a finite subspace $W_{\perp}^{M}$ using BEM-ECP.

The Dirichlet trace found by BEM-ECP is expressed in following terms:

(6.15) \begin{equation}\gamma_D^{\Omega} \vec{u} (\vec{x}) \approx\sum\limits_{i = 1}^{M} d_i \vec{f_i} (\vec{x}),\vec{f}_i \in W_{\perp}^{M}\end{equation}

where M – is the number of the basis functions (and their associated edges). As before, the levels of subdivision are denoted as x1, … , x8. For the cube, we used uniform meshes obtained by the subdivision of one another. The x1 mesh of the cube is made out of 12 triangles.

The numerical discrepancy for the Dirichlet trace is listed in Table 3 and computed similarly to (6.13) (we involve all the weights):

(6.16) \begin{equation}\operatorname{diff} =\left(\frac{\sum_{i = 1}^{M} \left(d_i - D_i\right)^2}{\sum_{i = 1}^{M} D_i^2}\right)^{\frac{1}{2}},\end{equation}

where $\operatorname{diff}$ – is the relative numerical discrepancy, $d_i$ – the decomposition weights (6.15) found via BEM-ECP with or without the integration by parts approach, $D_i$ – the weights obtained for the known solution via the least square method. In particular, $D_i$ are deduced as a result of the following minimisation problem:

(6.17) \begin{equation}\left \langle \sum\limits_{i = 1}^{M} D_i \vec{f_i} (\vec{x}) - \gamma_D^{\Omega} \vec{u} (\vec{x}),\sum\limits_{i = 1}^{M} D_i \vec{f_i} (\vec{x}) - \gamma_D^{\Omega} \vec{u} (\vec{x})\right\rangle \to \operatorname{min},\vec{f}_i \in W_{\perp}^{M}.\end{equation}

In (6.17), the scalar product (2.11) produces only a non-negative real number, so the minimisation problem makes sense.

Table 3. Numerical comparison for a box where $k = i$

As one can see in Table 3, the loss of convergence does indeed occur for fine meshes and is significantly larger than the one produced by the integration-by-parts approach.

We use the Stratton–Chu representation formula (2.6) to evaluate the solution on a segment along the OX line inside the box. This segment is coloured green in Figure 11.

Figure 11. Illustration of the box example.

Let us denote the solution computed through the Stratton–Chu representation formula as $\vec{U}$ – for that, we use the value of the Dirichlet trace computed via BEM-ECP ( $D_i$ weights for both the direct and by-parts integrations of the MDLP). The discrepancy of the solution given along the highlighted line in Figure 11 is illustrated in Figure 12.

Figure 12. Discrepancy for the finest mesh; $\vec{U}$ evaluated for both the direct and by-parts integrations of the modified double layer potential.

As one can see, the discrepancy is significantly larger for the case of direct integration of MDLP. The largest discrepancy is seen near the faces of the box.

We also plot the Dirichlet trace discrepancy, expressed as $\gamma_D\!\left(\vec{u}(\vec{x}) - \vec{U}(\vec{x})\right)\cdot\vec{c}_3$ , along the green segment belonging to the face with the normal vector directed opposite to the positive direction of OX and aligned with the direction of OZ; the illustration of the position of the output segment is given in Figure 13. The discrepancy itself is given in Figure 14.

Figure 13. Discrepancy for the finest mesh; $\vec{U}$ evaluated for both the direct and by-parts integrations.

Figure 14. Discrepancy for the finest mesh – values are given along the segment illustrated in Figure 13; $\vec{U}$ evaluated for both the direct and by-parts integrations of the modified double layer potential.

7 Conclusions

The numerical experiments presented in this work showed that our approach to MDLP integration based on the integration-by-parts formula provides numerical convergence for the supplied meshes.

Numerical experiments showed that our approach helps produce much better results at the edges and corners of the model geometry than one can produce with the direct integration of the MDLP.