Boundary Element Procedure for 3D Electromagnetic Transmission Problems with Large Conductivity

We consider the scattering of time periodic electro-magnetic fields by metallic obstacles, the eddy current problem. In this interface problem different sets of Maxwell equations must be solved in the obstacle and outside, while the tangential components of both electric and magnetic fields are continuous across the interface. We describe an asymptotic procedure, which applies for large conductivity and reflects the skin effect in metals. The key to our method is to introduce a special integral equation procedure for the exterior boundary value problem corresponding to perfect conductors. The asymptotic procedure leads to a great reduction in complexity for the numerical solution since it involves solving only the exterior boundary value problem. Furthermore we introduce a new fem/bem coupling procedure for the transmission problem and consider the implementation of the Galerkin elements for the perfect conductor problem and present numerical experiments.


Introduction
We present asymptotic expansions with respect to inverse powers of conductivity for the electrical and magnetical fields and report the algorithm of MacCamy and Stephan [12] which allows to compute the expansion terms of the electrical field in the exterior domain by solving sucessively only exterior problems (so-called perfect conductor problems) with different data on the interface between conductor (metal) and isolator (air). We solve these exterior problems numerically by applying the Galerkin boundary element method to first kind boundary integral equations which were originally introduced by MacCamy and Stephan in [11]. This system of integral equations on the interface Σ results from a single layer potential ansatz for the electrical field and has unknown densities namely a vector field and a scalar function on Σ which we approximate with lower order Raviart Thomas elements and continous piecewise linear functions on a regular, triangular mesh on Σ As in the two dimensional case, investigated by Hariharan [4,5] and MacCamy and Stephan [13], the asymptotic procedure gives for the computation of the solution of the transmission problem a great reduction in complexity since it involves solving only the exterior problem and furthermore only a few expansion terms must be computed. We describe in detail how to implement the boundary element method for the perfect conductor problem. As an alternative to the asymptotic expansions for the solution of the transmission problem we introduce a new finite element/boundary element Galerkin coupling procedure which converges quasi-optimally in the energy norm.

Asymptotic expansion for large conductivity and skin effect
Let Ω − be a bounded region in R 3 representing a metallic conductor and Ω + := R 3 \Ω − . Ω + representing air. The parameters ε 0 , µ 0 , σ 0 denote permittivity, permeability and conductivity. is assumed to have zero conductivity in Ω + with ε, µ, σ in Ω − . Let the incident electric and magnetic fields, E 0 and H 0 , satisfy Maxwell's equations in air. The total fields E and H satisfy the same Maxwell's equations as E 0 and H 0 in Ω + but a different set of equations in Ω − . Across the interface Σ := ∂Ω − = ∂Ω + , which is assumed to be a regular analytic surface, the tangential components of both E and H are continuous. E − E 0 and H − H 0 represent the scattered fields. All fields are time-harmonic with frequency ω. As in [12] we neglect conduction (displacement) currents in air (metal). Then, with appropriate scaling, the eddy current problem is (see [22,25]).
Problem (P αβ ): Given α > 0 and β > 0, find E and H such that; (1) Here α 2 = ω 2 µ 0 ε 0 and β 2 = ωµσ −iω 2 µε are dimensionless parameters, and β 2 = ωµσ > 0 if displacement currents are neglected in metal (ε = 0). The subscript T denotes tangential component and the superscripts plus and minus denote limits from Ω + and Ω − . At higher frequencies the constant β is usually large leading to the perfect conductor approximation. Formally this means solving only the Ω + equation and requiring that E T = 0 on Σ. If we let E and H denote the scattered fields, we obtain Problem (P α∞ ): Given α > 0, find E and H such that; Remark 1. There exists at most one solution of problem (P αβ ) for any α > 0 and 0 < β ≤ ∞ (see [18]).
We are interesting in an asymptotic expansion of the solution of problem (P αβ ) with respect to inverse powers of conductivity. With τ denoting the distance from Σ measured into Ω − along the normal to Σ the expansions reads: Here E n and H n are independent of β which is proportional to √ σ. The exponential in (5) and (6) represents the skin effect. Next we present from [12] these expansions for the half-space case where the various coefficients can be computed recursively. Note E 0 and H 0 in (3) and (4) is simply the perfect conductor approximation, that is, the solution of (P α∞ ). E n and H n in (3) and (4) can be calculated successively by solving a sequence of problems of the same form as (P α∞ ) but with boundary values determined from earlier coefficients. The E n and H n in (5) and (6) are obtained by solving ordinary differential equations in the variable x 3 . For the ease of the reader we present here for the half-space case Ω + = R 3 + i.e. x 3 > 0 and Ω − = R 3 − i.e. x 3 < 0 a formal procedure to compute E n , H n which was given by MacCamy and Stephan [12]. They substitute in (3)-(6) into (P αβ ) for Σ = R 2 and equate coefficients of β −n . Here we give a short description of their approach. Let χ = e √ −iβx 3 and decompose fields F into tangential and normal components with orthogonal component F ⊥ = e 3 × F, and unit vectors e i (i = 1, 2, 3). Then one computes with the surface gradient grad T for the rotation and Now setting E n = E n + ℓ n e 3 one obtains for and (11) Hence, equating coefficients of β 2 and β, respectively yields ℓ 0 ≡ 0, Now the gauge condition div E 0 = 0 implies ℓ 1 ≡ 0 and div E 0,x 3 = 0, hence E 1,x 3 = 0 and MacCamy and Stephan obtain in [12] with ℓ 1 = 0, h 0 = 0 E 0 = 0: and and For x 3 > 0, we have with curl E = H yields Equating coefficients of β −n one finds in x 3 > 0 curl E 0 = H 0 , curl E n = H n , n ≥ 0, (and correspondlying due to curl H = α 2 E) With the above relations the recursion process goes as follows. First one use (6.10) for n = 0 and (6.13), in [12], to conclude that is just the solution of (P α∞ ) which we can solve by the boundary integral equation procedure introduce in MacCamy and Stephan and revisited belov. But from (1) 3 we obtain Now the right side of (17) is known and easily computed. Then (1) 3 and (17) yield Therefore by (6.10), in [12], we have a new again solvable problem for (E 1 , H 1 ) which is just like (P α∞ ), that is but with new boundary values for E T as given by (18).
For the complete algorithm see [12].
A comparison with Peron's results (see Chapter 5 in [20]) shows that W cd Furthermore we see that the first terms in the asymptotic expansion of the electrical field for a smooth surface Σ derived by Peron coincide with those for the half-space x 3 = 0 investigated by MacCamy and Stephan, namely ℓ 0 = w 0 = 0, ℓ 1 = w 1 = 0, E 0 = W cd 0 = 0. Remark 3. Since due Theorem 5 in Chapter 3 of [19] there exists only one solution of the electromagnetic transmission problem for a smooth interface this solution can be compute by the boundary integral equation procedure below, when we assume that (19) holds. Then for the electrical field E obtained via the boundary integral equation system we have that in the tubular region Ω ± (δ) = {x ∈ Ω ± , dist(x, Σ) < δ} there holds for the remainders E is(cd) m obtained by truncating (3) and (5) Next we describe the integral equation procedure for (P αβ ) and (P α∞ ) from [12,25]. Throughout the section we require that This methods, like others, are based on the Stratton-Chu formulas from [22]. To describe these we need some notation. We will let n denote the exterior normal to Σ. Given any vector field v defined on Σ we have where v T , which lies in the tangent plane, is the tangential component of v.
We define the simple layer potential V κ for density ψ (correspondingly for a vector field) for the surface Σ by For a vector field v on Σ we define V κ (v) by (21) with v replacing ψ. We collect in the following lemma some of the well-known results about the simple layer potential V κ .
For problem (1) 2 , in Ω − the Stratton-Chu formula gives Similarly, for problem (1) For given n × H, n × E and n · E (23) yield a solution of (P α∞ ). But we know only n × E. The standard treatment of (P α∞ ) starts from (23), sets n × H = 0 and n · E = 0 and replaces −n × E by an unknown tangential field L yielding Then the boundary condition yields an integral equation of the second kind for L in the tangent space to Σ. The method (24) is analogous to solving the Dirichlet problem for the scalar Helmholtz equation with a double layer potential. But having found L it is hard to determine H T , or equivalently n × H, on Σ. Note calculating n × H on Σ involves finding a second normal derivative of V α (L). The method in [12] for (P α∞ ) is analogous to solving the scalar problems with a simple layer potential (see [9]). MacCamy and Stephan use (23) but this time they set n × E = 0 and replace n × H and n · E by unknowns J and M . Thus they take If they can determine J then in this case they can use Remark 4 to determine n × H, hence H T on Σ.
With the surface gradient grad T ψ = (grad ψ) T on Σ, the boundary condition in (1) and (25) imply, by continuity of V α , We note that for any field v defined in a neighbourhood of Σ one can define the surface divergence div T by As shown in [12, Lemma 2.3], there holds for any differentiable tangential field Setting divE = 0 on Σ yields therefore with (25)

FE/BE coupling
Next we present a coupling method for the interface problem (P αβ ) (see [1,2,6,7,19]). Integration by parts gives in Ω − for the first equation in (P αβ ) with γ N E = (curl E) × n, where K α is a smothing operator.
As shown in [12,Lemma 4.5] there exists a continuous map As shown in [11] the system of boundary operators on Σ (which is equivalent to (26) and (27)) is strongly elliptic as a mapping from where grad T (div T ) denote the surface gradient (surface divergence) and ∆ T the Laplace-Beltrami operator on Σ. Now, our fem/bem coupling method is based on the variational formulation: For given incident field E 0 on Σ find E ∈ H(curl, Ω − ), J ∈ H − 1 2 (Σ) and M ∈ H ∀v ∈ H(curl, Ω − ), j ∈ H − 1 2 (Σ), m ∈ H 1 2 (Σ). In order to formulate a conforming Galerkin scheme for (31) we take subspaces where A is the operator given by the left hand side in (31),

The Galerkin system (32) is uniquely solvable in
h and there exists C > 0, independent of h, where (E, J, M ) and (E h , J h , M h ) solve (31)-(32) respectively.
Proof. First we note that system (31) is strongly elliptic in X which follows by considering A as a system of pseudodifferential operators (cf. [11]). The only difference to [11] is that here we have additionally the first equation in (31). If we note ∆E = curlcurlE − graddivE and take divE = 0 we have that the principal symbol of A has the form (with |ξ| 2 = ξ 2 where (E 1 , E 2 ) = E T and E 3 is perpendicular to x 3 = 0.
Obviously the two subblocks are strongly elliptic (see [11] for the lower subblock). Assuming that (α, √ iβ) is not an eigenvalue of P αβ we have existence and uniqueness of the exact solution. Due to the strong ellipticity of A there exists a unique Galerkin solution and the a priori error estimate holds due to the abstract results by Stephan and Wendland [21]. 5 Galerkin procedure for the perfect conductor problem (P α∞ ) Next we consider the implementation of the Galerkin methods (see [3,19,23,25]) and present corresponding numerical experiments for the integral equations (26) and (27). These experiments are performed with the program package Maiprogs, cf. Maischak [15,17], which is a Fortran-based program package used for finite element and boundary element simulations [16]. Initially developed by M. Maischak, Maiprogs has been extended for electromagnetics problem by Teltscher [24] and Leydecker [10]. We will investigate the exterior problem (P α∞ ) by performing the integral equations procedure (26) and (27): Testing against arbitrary functions j ∈ H − 1 2 (Σ) and m ∈ H 1 2 (Σ) in (26) and (27), we get Partial integration in the second term of (35) 1 shows that the formulation (35) is symmetric: By definition of symmetric bilinear forms a, c, of the bilinear form b and linear form ℓ through the variational formulation has the form: We now proceed to finite dimensional subspaces R h ⊂ H − 1 2 (Σ) of dimension n and M h ⊂ H for all ψ k and ϕ l , 1 ≤ k ≤ n, 1 ≤ l ≤ m.

With matrices and vectors
A := (a(ψ i , ψ k )) i,k ∈ C n×n , (39) has also the form We have considered with {ψ i } n i=1 a basis of R h and {ϕ j } m j=1 a basis of M h . These functions, are chosen as piecewise polynomials. To win these bases, we consider suitable basis functions locally on the element of a grid, i.e. on each component grid.
with N elements, and let { ψ i } n i=1 and { ϕ j } m j=1 respectively bases on a square reference element Σ. The local basis functions on an element Σ k are each It should therefore be calculated first A := (a(ψ js , ψ iz )) iz ,js ∈ C n×n , where ψ js or ψ iz are the basics function of R h and We are dealing in this implementation with Raviart-Thomas basis functions. The transformation of these functions requires a Peano transformation ψ k,i = 1 a 2 ), detA k is calculated by detA k = (a 1 × a 2 ) · a 1 × a 2 a 1 × a 2 . The Peano-transformation of the local basis functions to the basic functions on the reference element then gives with x = a k + A k x and y = a l + A l y, and referent element Σ.
The calculation of the integrals with Helmholtz kernel G α is not exact. We consider the expansion of the Helmholtz kernel in a Taylor series. There holds The first term are singular for x = y and the correspondly integral are treated by analytic evaluation in Maiprogs, cf. Maischak [14,15,17] , but the integrals of all other summands can be calculated sufficiently well by Gaussian quadrature. We compute with ζ −1 ψ = ζ described above, and ζ −1 ϕ , the analogously defined map for the basic functions of M h . While a transformation of the scalar basis functions is not required, the transformation of the surface divergence of Raviart-Thomas elements is carried out by ∇ T · ψ k,i = with y = a k + A k y and x = a l + A l x.
The calculation of c(ϕ i , ϕ j ) is similar to the above-mentioned case.
The calculation of the right hand side appears simple at first glance, since there are no single layer potential terms. Howewere we must compute the right hand side with quadrature.
The quadrature of an integral over f on the reference element is determined by the quadrature points x x,y , and the associated weights w x,y = w x · w y , which are processed in x and y direction. We perform the two-dimensional quadrature as a combination of one-dimensional quadratures in each x and y direction, and we use here the weights from the already implemented one-dimensional quadrature formula. With n x quadrature points in x-direction, and n y quadrature points in y-direction, then the quadrature formula reads: The quadrature points on the square reference element and the corresponding weights for Gaussian quadrature are implemented in Maiprogs already. For triangular elements, we use Duffy transformation.
We will now calculate the right hand side in the Galerkin formulation, i.e. the linear form ℓ, applied to the bases functions ψ i , i = 1, . . . , n. The quadrature takes place on the reference element. We decompose the global into local basis functions and then use the Peano-transformation for the Raviart-Thomas functions. It is therefore Applying (45) leads with n x = n y := n to with x i,j = a k + A k x i,j . As before, the task is carried out by looping through all grid components, and the values are added to the entries for each of its base function. The electrical field can be calculated by We have for the first term in (47) with (38) 1 Then using Peano-transformation we have For the second term in (47) we have with The calculation of H ± T is done as follows (compare Remark 4 (v)) Example 1. Here, we consider one example to test the implementation. As domain we take the cube Ω − = [−2, 2] 3 , and we now want to test the Galerkin method in (37). We choose the wave number α = 0.1 (or α = 0.5, 1.5), and the exact solution and where n = (n 1 , n 2 , n 3 ) denotes the outer normal vector at a point on the surface Σ = ∪ 6 k=1 Σ k . We can write each term of equation (26) as: and Then, from (26), (54) and (55) holds We use different values of α for our investigation. In Table 1 we present the results of the errors in energy norm and L2-norm for α = 0.1, 0.5, 1.5 for the uniform h-version with polynomial degree p = 1. In Figures 1 and 2 we compare the h-version with different α. The exact norm is known by extrapolation for α = 0.1 is |C| = 8.580798, for α = 0.5 is |C| = 1.6171534, and for α = 1.5 is |C| = 1.8042380. Here C = Re E 0 T , J and C h = Re E 0 T , J h (see [8]). The exact L2-norm is known by extrapolation for α = 0. Let as compare our numerical convergence rates above for the boundary element methods obtained in the above example with the theoretical convergence rates predicted by Theorem 1. Note that we have implemented the boundary integral equation system (26), (27) and note the strongly elliptic system (30), where convergence is garanteed due to Theorem 1. Nevertheless our experiments show convergence for the boundary element solution, but with suboptimal convergence rates. Theorem 1 predicts (when Raviart-Thomas elements are used to approximate J and piecewise linear elements to approximate M ) a convergence rate of order η = 3 2 in the energy norm for smooth solutions J and M . Our computations depend on the parameter α which is a well-known effect with boundary integral equations where it may come to spurious eigenvalues diminishing the orders of the Galerkin approximations. Due to the cube Ω − = [−2, 2] 3 the numerical solution might become singular near the edges and corners of Ω − ; hence the Galerkin scheme converge suboptimally.  Table 1: Errors in L 2 -norm and energy norm with respect to the degrees of freedom for α = 0.1, 0.5, 1.5.

N DOF
Next, we apply the boundary element method above to compute the first terms in the asymptotic expansion of the electrical field considered in subsection 2 (Remark 1). In this way we obtain good results for the electrical field at some point away from the transmission surface Σ by only computing a few terms in the expansion.
Algorithm for the asymptotics of the eddy current problem: 1. First solve the exterior Problem (P α∞ ) by integral equations (26)  5. E = E 0 + β −1 E 1 + β −2 E 2 + R m , where E 0 is the solution of the step 1 and E 1 and E 2 are solutions of step 3.