Modeling of field singularities at dielectric edges using grid based methods

Electric field singularities at sharp metallic edges or at a dielectric contact line can be described analytically by asymptotic expressions. The a priori known form of the field distribution in the vicinity of these edges can be used to construct numerical methods with improved accuracy. This contribution focuses on a modified Finite Integration Technique and on a Discontinuous Galerkin Method with singular approximation functions. Both methods are able to handle field singularities at perfectly electric conducting as well as at dielectric edges. The numerical accuracy of these methods is investigated in a number of simulation examples including static and dynamic field problems.


Introduction
One of the fundamental principles of electrodynamics states that the electromagnetic field energy within a finite domain is finite. This remains valid even when the electromagnetic field within the domain becomes singular as may be the case at perfectly conducting edges and at dielectric contact lines. Electromagnetic field singularities are restricted, according to the finite energy principle, to be no stronger than ρ −1+χ , where χ > 0 and ρ is the distance from the edge. Asymptotic expressions for the singular fields at a perfectly conducting edge have been early presented in the literature (Meixner, 1972;Hurd, 1976). More than a decade later (Olyslager, 1994), singular field solutions for the general case of a contact edge between a number of dielectric and/or magnetic materials were found.
The idea of using the a priori known asymptotic behavior of singular fields in numerical simulations to improve numerical accuracy appears natural. This approach has been proposed by several authors in the context of different dis-Correspondence to: C. Classen (classen@tet.upb.de) cretization methods. In Mur (1981) singular correction terms are incorporated into the Finite Difference Time Domain (FDTD) method for the time domain modeling of high frequency problems, see also Beilenhoff and Heinrich (1993), Przybyszewski and Mrozowski (1998). In the context of Finite Element Methods (FEM), several approaches based on the use of specialized scalar and vector basis functions incorporating the singular field behavior have been proposed (Webb, 1988;Graglia and Lombardi, 2004). Recent developments include the Extended Finite Element Method (XFEM) (Moës et al., 1999) and the Partition of Unity Finite Element Method (PUFEM) (Melenk and Babuška, 1996).
The main difficulty with most of the discretization methods using singularity correction techniques is the increased numerical and implementation complexity. This is, e.g., the case for the FEM where specialized singular functions, fullfilling global continuity conditions, must be used in the approximation. The continuity condition imposes an important constraint which limits the flexibility and, as indicated by numerical results (Chahine et al., 2006), the accuracy of the method. The application of this approach for geometrically complicated problems (e.g. for sharp edges interfacing at several dielectric and perfectly conducting domains) is difficult. Furthermore, from the numerical point of view, it may be advantageous to apply a standard discretization method which includes a simple to implement (probably less accurate but numerically more efficient) technique for singular field correction.
In this contribution, two new numerical approaches to handle field singularities are proposed. The first one consists in a modification of the Finite Integration Technique (FIT) (Weiland, 1977) for right angle corner edges in Cartesian grid problems. The second approach uses a flexible singular basis approximation on unstructured grids based on the Discontinuous Galerkin (DG) method. The numerical accuracy of these approaches is investigated for a number of singular field problems with different types of field singularities.

Singular field solutions
In the direct vicinity of edges of perfectly conducting -or dielectric wedges, the asymptotic form of the electromagnetic field solution is given by The singularity index t obeys the edge condition and has to be in the range t ∈ (0,1). It is possible to determine t exactly. A simple analytical procedure to compute t in the general case of a contact edge interfacing to a number of dielectric, magnetic and perfectly conducting materials can be found in (Olyslager, 1994).

FIT with singularity correction
The numerical field unknowns in the FIT formulation for electrostatics problems are given by the electric grid voltages e and dielectric grid fluxes d defined by respectively. At the discrete grid level, the relation between these quantities needs to be approximated by the matrix equation where M ;FIT is the permittivity matrix of the FIT formulation. In the case of staggered Cartesian grids (see Fig. 1a), a typical approximation for M ;FIT consists in a diagonal matrix with entries corresponding to the one-to-one relation between the k-th voltage e and the k-th flux d on the grid. In Eq. (5), y and z are the dual grid lengths in y-and z-directions (face A), respectively, and x is the primary grid length in x-direction. In the case of singular field problems, it has been early realized that approximation (5) leads to an extremely slow numerical convergence. The overall numerical accuracy, even at grid points far away from sharp metallic or dielectric edges is drastically reduced compared to regular field problems. A new variant of a singular field correction approach has recently been proposed in the context of FIT (Classen et al., 2010). Hereby, the idea is to derive a modified permittivity matrix by using the singular field behavior (1) in both integral expressions (3). Thus, instead of using the purely geometrical approximation (5), the modified coefficients in the material matrix (for each dual pair of edges and faces in the grid) become where the correction factor K k involves a numerical integration of the asymptotic expressions (1). Note that in this approximation, the analytically known singularity index, t, is used depending on geometry and on the electrical properties of the materials adjacent to the edge. Furthermore, this correction can be applied also at grid points not lying on the singular edge but which are sufficiently close to be influenced by the field singularity. The approach taken in this work is to apply the correction (6) within a small sphere surrounding the singular point (see Fig. 1b). The radius r of this sphere represents a free parameter of the formulation. The permittivity matrix M with singularity correction can be inserted as usual in the discrete FIT equations, e.g., for electrostatics as where S is the source operator, q are the volume charges and v the nodal potentials. Note, however, that the same approach can be used to derive a corrected permeability matrix for magnetic field and high frequency problems.

DG with singular basis functions
The local discontinuous Galerkin method is chosen in this work because of its flexibility in the choice of approximation spaces. Following the derivation by Arnold et al. (2002), in the electrostatic case, the potential V and the electric field E are approximated within every element K of the mesh by local approximation spaces P (K) and (K), respectively. The standard choice is P (K) = P p ; the space of polynomials of degree at most p, and (K) = (P p ) 2 . Now, let V s (ρ,ϕ) = ρ t (ϕ) denote a function corresponding to the asymptotic electric potential solution at the singular point.
Although being finite at every point, V s is referred to as singular function due to its low regularity. The approximation spaces for all elements lying completely or partially inside a region of radius r around the singular point can be easily enriched with additional basis functions of the type V s . Let < V s > denote the space of scalar multiples of V s , then for these elements the modified local approximation spaces are given bỹ Note that since the DG approximation is inherently discontinuous, the singular basis functions can be simply defined according to the asymptotic field expressions (1). There is no need to adapt these functions to a particular mesh in order to fulfill the continuity condition as is the case for the FEM. This flexibility represents also the main advantage of the method for the solution of singular field problems. Omitting details on the involved finite element spaces V h and h , the numerical fluxes are defined aŝ In Eqs. (10), (11), v h ∈ V h and τ h ∈ h approximate scalar and vector quantities, respectively. Furthermore, {·} denotes the mean value and · the jump whereas C 11 is a suitably chosen stabilization parameter, see also Castillo et al. (2000). For the numerical integration of terms involving singularities, higher order (or even adaptive) quadrature rules should be used. Using this approximation procedure, the discrete DG equations can be obtained as a matrix equation which is formally similar to Eq. (7): where S DG , M ;DG and v are the DG counterparts of the source, the permittivity matrix and the potential, respectively, and C DG is a stabilization matrix.

Numerical results
The FIT edge correction approach and the DG with singular basis functions are validated for a number of simulation setups.

PEC corner singularity
First tests are performed using a 2-D electrostatics problem containing a sharp perfectly conducting edge. The parameters in Fig. 2 are chosen as: a = 1 m, b = 0.5 m, V 0 = 1 V, α = π/2 and the boundary values at the surrounding box with edge length b are imprinted using the potential: Equation (13) represents the exact solution of the problem with β = 2π − α, and E n given by .
The singularity index t and the functions given in Eqs. (1) and (2), can be derived using the general formalism by Olyslager or -in this special case -directly by the derivative of Eq. (13). The problem is discretized with rectangles for the FIT and triangles for the DG method, respectively. The FIT results for different values of the parameter r are plotted in Fig. 3. The maximum norm, comparing numerical results to analytic results at the grid points, is plotted over the mesh step size h max . The green curve (+) shows the results of standard FIT without edge correction while the purple curve ( ) represents the results using a standard edge correction scheme as given in Beilenhoff and Heinrich (1993). Obviously the total error is decreased, although the order of convergence is not improved. The blue curve ( * ) represents the new edge correction approach which is only applied to the grid edges directly located at the singular point (as it is also done in the standard correction scheme). In this case, the results can be further improved. It can be clearly observed in these cases that no second order convergence can be achieved as it would be the case for a regular field problem. On the other hand, correcting all grid edges lying at least partially inside a fixed area with radius r (see Fig. 1b), second order convergence can be restored. As observed in Fig. 3, in order for this to occur, the edge length has to become smaller than the parameter r.
The results of the DG method with singular basis functions and first order polynomials are depicted in Fig. 4. Also here, a fixed area of enrichment with singular basis functions is needed to obtain a higher order of convergence. Using singular basis functions only in elements having a common point with the singularity (blue curve) ( * ), does not improve the order of convergence, see also Laborde et al. (2005).

Dielectric edge singularity
The setup consists in a dielectric loaded capacitor (see Fig. 1b). The total extension of the domain is 0.5 m in yand x-direction. At the top and bottom of the domain Neumann boundary conditions are imposed; on the left and right Dirichlet boundary conditions with the fixed potentials ±1V are applied. The dielectric inset has the permittivity r = 10. In this case the singularity is given by V s (ρ,ϕ) = ρ t (ϕ), with t ≈ 0.73.  It should be mentioned that both methods necessitate a slight modification when different materials surrounding the edge are present. The angular part of the singular function (ϕ) can be written in the form where the constant C, depending on the permittivities, the angular material distribution and the singularity index is different in each material. Therefore the singularity correction and the singular approximation functions have to be adapted to each material, refer to Olyslager (1994) for details. We compare the numerically computed electrical capacity to the one calculated using a very fine mesh and a higher order FEM simulation. The FIT results are presented in Fig. 5; the DG results for first order polynomials in Fig. 6. In both cases, the convergence behavior concerning mesh refinement is identical to the previous example. These results indicate that the described corrections are well suited for treating dielectric type singularities just as in the PEC case.

Dielectric contact singularity
A more complicated setup concerns singularities arising at the contact line between three or more dielectric materials. To avoid additional effects from a non-conformal geometry discretization in a rectangular FIT grid, this setup is investigated only for the DG method on a triangular grid. The case of three dielectric wedges, with relative permittivities 1 = 1, 2 = 80 and 3 = 4 is depicted in Fig. 7. For illustration purposes, also an exemplary triangular grid as well as the region used for DG singular basis enrichment are shown. The singularity is given by V s (ρ,ϕ) = ρ t (ϕ), with t ≈ 0.61, corresponding to a contact angle α = 120 • . The boundary conditions are identical to the previous example and a = 0.5 m. The reference result is obtained by a 2-D FEM calculation using a heavily refined mesh in the vicinity of the contact line. Also in this case, numerical results indicate that the error in the electric potential converges at second order when the DG approach with singular basis functions (within a fixed region) and first order polynomials is used (see Fig. 8).

3-D eigenvalue problem with reentrant corner
The last example is the solution of the curl-curl eigenvalue problem for a cavity with reentrant corner. The structure is a hollow, PEC bounded rectangular box (see where C FIT is the discrete curl matrix of FIT. The convergence of the first resonance frequency with respect to the grid size is shown in Fig. 10. As observed in the previous examples the standard correction technique ( ) and the new edge correction ( * ) improve the accuracy, meanwhile the order of convergence remains unchanged.
x y z Fig. 9. Arrowplot of the electric field distribution for the eigenvalue problem.

Conclusions
Two numerical approaches for the simulation of singular field problems have been presented. The first consists in a special modification of the material matrices of FIT, whereas the second represents a DG method with singular basis functions. Both methods lead to an immense improvement in numerical accuracy in the case of PEC as well as for dielectric edge singularities. Numerical simulations for a number of setups with typical electromagnetic field singularities show, in particular, that the optimal convergence order for regular field problems can be fully recovered when the respective singular corrections are applied within a small but finite region of influence surrounding the singularity.