A study of three-dimensional edge and corner problems using the neBEM solver

Abstract

The previously reported neBEM solver has been used to solve electrostatic problems having three-dimensional edges and corners in the physical domain. Both rectangular and triangular elements have been used to discretize the geometries under study. In order to maintain very high level of precision, a library of C functions yielding exact values of potential and flux influences due to uniform surface distribution of singularities on flat triangular and rectangular elements has been developed and used. Here we present the exact expressions proposed for computing the influence of uniform singularity distributions on triangular elements and illustrate their accuracy. We then consider several problems of electrostatics containing edges and singularities of various orders including plates and cubes, and L-shaped conductors. We have tried to show that using the approach proposed in the earlier paper on neBEM and its present enhanced (through the inclusion of triangular elements) form, it is possible to obtain accurate estimates of integral features such as the capacitance of a given conductor and detailed ones such as the charge density distribution at the edges/corners without taking resort to any new or special formulation. Results obtained using neBEM have been compared extensively with both existing analytical and numerical results. The comparisons illustrate the accuracy, flexibility and robustness of the new approach quite comprehensively.

Introduction

One of the elegant methods for solving the Laplace/Poisson equations (normally an integral expression of the inverse square law) is to set up the boundary integral equations (BIE) which lead to the moderately popular boundary element method (BEM). In the forward collocation version of the BEM, surfaces of a given geometry are replaced by a distribution of point singularities such as source/dipole of unknown strengths. The strengths of these singularities are obtained through the satisfaction of a given set of boundary conditions that can be Dirichlet, Neumann or of the Robin type. The numerical implementation requires considerable care [1] because it involves evaluation of singular (weak, strong and hyper) integrals. Some of the notable two-dimensional (while all the devices are three-dimensional by definition, useful insight is often obtained by performing a two-dimensional analysis) and three-dimensional approaches used to evaluate the singular integrals are discussed in [1], [2], [3], [4], [5], [6], [7] and the references in these papers. It is well-understood that many of the difficulties in the available BEM solvers stem from the assumption of nodal concentration of singularities which leads to various mathematical difficulties and to the infamous numerical boundary layers [8], [9], [33] when the source is placed very close to the field point ([2] and Refs. [4], [5], [6] therein). While mathematical singularities (that occur when the source and field points coincide) have been shown to be artifacts, several techniques have been used to remove difficulties related to physical or geometrical singularities (that occur when boundaries are degenerate, i.e., geometrically singular, or due to a jump in boundary conditions) such as Gaussian quadrature integration, mapping techniques for regularization, bicubic transformation, nonlinear transformation and dual BEM techniques [8]. The last technique seems to be a popular one and capable of dealing with a relatively wide range of similar problems.

Departing from the approaches mentioned in the above references and many more to be mentioned below, we had shown in an earlier paper [10] that many of these problems can be eliminated or reduced if we adopt a new paradigm in which the elements are endowed with singularities distributed on them, rather than assuming the singularities to be concentrated at specific nodal points. Despite a large body of literature, closed form analytic expressions for computing the effects of distributed singularities are rare [11], [12], complicated to implement and, often, valid only for special cases [13], [14], [15]. For example, in [11], the integration of the Green function to compute the influence of a constant source distribution is modified to an “n-plane” integration. The evaluation of this integration involves co-ordinate transformations and the resulting expressions are rather complicated. In [12], the Gauss–Bonnet concept is used in which the panel is projected onto a unit sphere and the solid angle is determined from the sum of the induced angles. The procedure and the resulting expressions are neither simple, nor easy to implement in a computer code. In fact, possibly due to these difficulties, these approaches have remained relatively unpopular and even in very recent papers it is maintained that for evaluating the influence due to source distributed on triangular elements in a general case, one must apply non-analytic procedures [14]. Thus, for solving realistic but difficult problems involving, for example, sharp edges and corners or thin or closely spaced elements, introduction of special formulations (usually involving fairly complicated mathematics, once again) becomes a necessity [8], [16], [33]. These drawbacks are some of the major reasons behind the relative unpopularity of the BEM despite its significant advantages over domain approaches such as the finite-difference and finite-element methods (FDM and FEM) while solving non-dissipative problems [17], [18].

The Inverse Square Law Exact Solutions (ISLES) library developed in conjunction with the nearly exact BEM (neBEM) solver [10], in contrast, is capable of truly modeling the effect of distributed singularities precisely and, thus, is not limited by the proximity of other singular surfaces or their curvature or their size and aspect ratio. The library consists of analytic solutions for both potential and flux due to uniform distribution of singularity on flat rectangular and triangular elements. These close-form exact solutions, termed as foundation expressions, are in the form of algebraic expressions that are long but without complications and are fairly straight-forward to implement in a computer program. In deriving these foundation expressions, while the rectangular elements were allowed to be of any arbitrary size [10], [19], the triangular element was restricted to be a right-angled triangle of arbitrary size [20], [21], [22]. Since any real geometry can be represented through elements of the above two types (or by the triangular type alone), this library has allowed us to develop the neBEM solver that is capable of solving three-dimensional potential problems involving arbitrary geometry. It may be noted here that any non-right-angled triangle can be easily decomposed into two right-angled triangles. Thus, the right-angled triangles considered here, in fact, can take care of any three-dimensional geometry.

A set of particularly difficult problems to be dealt with by BEM is one that contains corners and edges and, in this work, we will attempt to solve several problems belonging to this set. The perfectly conducting bodies studied here are unit square plate, L-shaped plate, cube, L-shaped volume and two rectangular plates meeting at various angles and creating an edge. Besides being interesting and difficult, these solutions can have significant applications in micro-electromechanical systems (MEMS), nano-devices, atomic force microscopy (AFM), electro-optical elements, micro-pattern gas detectors (MPGD) and many other disciplines in science and technology. For these problems, it is important to study integral features such as the capacitance of the conductors, as well as detailed features such as the charge density, potential and flux on various surfaces of these objects including regions close to the geometric singularities. While several approaches including FDM, FEM, BEM and its variants such as the surface charge method (SCM) and various implementations of the Monte-Carlo technique (often coupled with Kelvin transformation) have been used to study these problems, only the latter two approaches, namely, BEM and Monte-Carlo technique are found to possess the precision necessary to model the curiously difficult electrostatics with acceptable levels of accuracy [23]. The volume discretization methods are known to be unsuitable because of the open nature of the problem and the inadequate representation of edge and corner singularities. Methods using Kelvin inversion (or quadratic inversion), although accurate, have been found incapable of handling planar problems. It may be noted here that despite the usefulness of two-dimensional analysis, there are an overwhelming number of problems that need to be addressed in three dimensions. As a result, several interesting approaches have been developed to analyze edge and corner related problems in complete three dimensions, without even the assumption of axial symmetry. In order to maintain applicability in the most general scenarios, in this work we will deal with the problems of edge and corner as truly three-dimensional objects even when comparing the results with two-dimensional analytic ones.

The problem of estimation of capacitance of square plate and cube raised to unit volt has been studied by an especially large number of workers using entirely different approaches. In fact, these have been considered to be some of the major unsolved problems of electrostatics, of which a solution is said to have been given by Dirichlet and subsequently lost. One of the more popular numerical approaches used to explore these problems is the BEM/SCM [13], [24], [25], [26], [27], [28]. Some of the solution attempts are more than a century old and yielded quite acceptable results. The later studies [13], [27], [28] used the mesh refinement technique and extrapolation of N (the number of elements used to discretize a given body) to infinity in order to arrive at more precise estimates of the capacitance. In order to carry out this extrapolation, uniform charge density scenario has been maintained so that the form of charge distribution on individual segments becomes irrelevant. According to [28], it is justified to use uniform charge density on individual elements because increase in complexity through the use of non-uniform charge density ultimately does not lead to computational advantage. In [28], the author mentions that for the cube, the element sizes are made such that the charge on each element remains approximately a constant (independent of its distance from an edge) since this arrangement is found to give the most accurate results.

The problem of estimation of the order of singularities at edges and corners of different nature is strongly coupled with the problem of estimation of integral properties such as the capacitance of conductors of various shapes. Thus, this problem, which also has importance in relation to other areas of science and technology as discussed earlier, has attracted the attention of a large number of workers as well. Here, fortunately, some analysis has been possible using purely theoretical tools [29], [30], at least for two-dimensional cases. In [31], the authors used a singular perturbation technique to obtain the singularity index at inside and outside corners of a sectorial conducting plate. Similarly, corner singularity exponents were numerically obtained in [28]. According to [32], it was possible to achieve accuracy of one in million through the use of FEM approximations for both electrodes and surface charge density, in addition to proper handling of edge and corner singularities. In this investigation, the Fichera's theorem was used to correctly describe the peculiarities of surface charge density behavior in the vicinity of the electrode ribs and tips.

According to another recent work [33], low-order polynomials used to represent the corners and edges lead to errors in estimation of the derivatives of the potential, and that is the crux of the problem. Despite its advantages (no prior knowledge of singular elements and the order of singularity is necessary) and good convergence characteristics, the mesh refinement approach has been mentioned to be less accurate than two other methods, namely, (1) singular elements, and (2) singular functions. Among these, beforehand knowledge of the location and behavior of the singularity in terms of the order of singularity is required for (1). The singular function approach (2) performs the best since it uses both the above and also information on singularity profile (corresponding to the eigenvalues and eigenvectors of the given geometry). Unfortunately, this approach requires complex formulation and is generally restricted to two-dimensional. Since the singular element approach uses the location information and only the order of the singularity, it is more versatile and popular. Thus, [33] implements the method of singular elements in three dimensions while attempting to use proper shape function/interpolation function to model correct singular behavior at corners and edges. It is shown that singular elements are needed to accurately capture the behavior at singular regions, such as sharp corners and edges, where standard elements fail to give an accurate result. Unfortunately, singular elements must be defined before it is possible to apply either the singular function or the singular element approach. This is a non-trivial task since for a realistic device, there can easily be thousands, if not millions, of elements involved. According to [34], the manual classification of boundary elements based on their singularity conditions is an immensely laborious task if not outright impractical. In [34], the authors developed an algorithm to extract the regions where singularity arises by querying the geometric model for convex edges based on geometric information of the model. The associated nodes of the boundary elements on these edges were then retrieved and categorized according to different types of singularity configuration. The algorithm developed was implemented in the PATRAN command language (PCL) on the MSC/PATRAN platform [35] which is an industry standard finite-element pre- and post-processor allowing a high degree of customization. In order to determine the order of singularity, two-dimensional results have been used directly [29], [36] for edges under the assumption that the point lies sufficiently far from any corner. In fact, in [34], any point other than the vertices has been considered to be far enough. This approach is quite unlikely to be very accurate, especially when the dimensions of the devices are small enough so that no point may even be considered to be far enough. Later, while discussing the results obtained using our proposed approach, we will return to this issue again. A study on the effect of bias ratio (a ratio of the largest element length to the smallest element length, similar to r-mesh refinement) has led the authors [33] to state that while for normal elements this ratio should be around 4:1 for a given problem, the singular element approach works better with lower bias ratio.

Besides the above procedures to identify the singular elements and their orders, attempts have also been made to place nodes at optimal positions while solving problems involving edges and corners. For example, in fracture mechanics [37], for handling singularities of the order 0.5, mid-side nodes in quadratic elements are normally shifted to quarter point positions. Several other techniques such as least square, constrained displacement of side nodes adjacent to a crack-tip [38] have also been used to find the optimum position of the mid-side nodes for handling singularities of other orders. In [39], the authors have subsequently shown that this approach of shifting mid-side nodes is likely to produce singularities of order 0.5 only and it is not suitable to impose arbitrary singularity in isoparametric elements by simply shifting side nodes to assumed positions.

In addition to BEM approaches, work has also been carried out using radically different approaches, for example, the Monte-Carlo methods. An impressive array of work exists providing very accurate estimates of capacitance and variation of charge density near edges and corners [23], [40], [41], [42], [43], [45], [46]. According to [23], the use of BEM introduces unnatural estimate of the charge density distribution—not for shapes with smooth contours (disks or spheres), but for plates and cubes. It has been mentioned that the situation is worst for the corner singularities of plates in which case there are no other surfaces present (as in a cube) to weaken the order of singularity at the corner. We will discuss this issue when results for square plates are presented below. The other point is that the BEM cannot satisfy the boundary conditions at or, at least, close to the edge because, the collocation points in BEM do not match the boundary of the device. In fact, it is mentioned that the solutions become unstable when the collocation points are shifted away from the centroid of the elements. This notion, as mentioned in the earlier paragraph, is not without counter-examples. Moreover, this is a point that we plan to take up in a future communication where we hope to illustrate that for the proposed formulation, shifting the collocation points to non-centroid locations does not lead to numerical instabilities. In addition, the approach of extrapolating capacitance has been criticized on the ground that they do not match for different amounts of shift. It has been mentioned that no formal error analysis exists for methods other than FDM and the extrapolation is purely empirical in nature. It has been observed that the apparent high accuracy may be illusory in nature and citing [47], it has been emphasized that situation can become even worse by attempting extrapolation of results obtained using non-equivalent meshes.

By developing a model that incorporates the truly distributed nature of sources/doublets/vortices on surfaces of three-dimensional geometries, we have recently shown [10], [19], [22] that it is possible to use the same formulation for studying a very wide range of problems (multi-scale, involving multiple layers of dielectric materials) governed by the Poisson's equation. Recently, we have extended the new formulation with the capability of including triangular elements as an option for discretizing arbitrary three-dimensional bodies [20], [21]. Here, we present the expressions to evaluate the exact values of potential and fluxes at any arbitrary point due to uniform singularity distributed on right-angled triangular element. These expressions have been included as additional functions of the ISLES library and subsequently used in the neBEM solver in the manner usual to, probably, most of the BEM solvers available. Besides presenting the expressions, we have also shown results to illustrate the accuracy of the expressions under various circumstances.

In this report, we have presented studies on the electrostatic configuration of several three-dimensional bodies, all of which contain corners and/or edges. The classic benchmark problems of estimating the capacitance of a unit square plate and unit cube raised to unit volt have been addressed. Electrostatics of generic shapes such as L-shaped plate and L-shaped three-dimensional conductors has also been analyzed. In all the above cases, the order of the singularity distribution near the edges and corners has been estimated in addition to estimating the capacitance of the conductors themselves. The singularity distributions obtained have been compared with theoretical and numerical studies carried out by earlier workers. The variation of the singularity along an edge between two corners at its ends has been studied, probably for the first time. Finally, the well-known problem of electrostatic configuration of two planes intersecting at different angles has been addressed. The two-dimensional counter-part of this problem is known to have analytic solution and has even been discussed in several textbooks on electromagnetics [29], [30]. Although accessible analytically, this problem seems to have been rarely solved using numerical techniques [48]. This benchmark problem is known to be a difficult one and, in order to test the proposed approach under difficult circumstances, we have computed electrostatic properties of a three-dimensional analogue of this problem close to the point of intersection for a very wide range of angle of intersection. Following the above studies, we have come to the conclusion that the proposed approach is capable of solving critical multi-scale problems governed by the Poisson's equation in a rather straight-forward manner. While higher bit accuracy, improved evaluation of transcendental functions, adaptive mesh generation and parallelization is expected to be of significant help, no special mathematical treatment or new formulation has been found to be necessary to deal with problems involving corners/edges and extremely closely spaced surfaces.

According to [23], using BEM, it is difficult to obtain physically consistent results close to these geometric singularities. Wild variations in the magnitude of the charge density have been observed with the change in the degree of discretization, the reason once again being associated with the nodal model of singularities. In contrast, using neBEM, we have obtained very smooth variation close to corner. Presence of oscillations seemingly acceptable to [23] has not occurred. In fact, oscillations close to edges and corners considered in this work seem to indicate numerical inaccuracy and have been treated accordingly. In addition to the shape, the magnitudes of the charge density have been found to be consistently converging to physically realistic values. These results clearly indicate that since the foundation expressions of the solver are exact, it is possible to find the potential and flux accurately in the complete physical domain, including the critical near-field domain using neBEM. In addition, since the singularities are no longer assumed to be nodal and we have the exact expressions for potential and flux throughout the physical domain, the boundary conditions no longer need to be satisfied at special collocation points such as the centroid of an element. Although consequences of this considerable advantage are still under study, it is expected that this feature will allow neBEM to yield even more accurate estimates for problems involving corners and edges since it should be possible to generate an over-determined system of equations by placing extra collocation points near the edges/corners. This will also allow the method to satisfy the boundary conditions of a given geometry at its true boundaries.

It should be noted here that the exact expressions for triangular elements consist of a significantly larger number of mathematical operations than those for rectangular elements presented in [10]. Thus, for any solver based on the ISLES library, it is more economical to use a mixed mesh of rectangular and triangular elements using rectangular elements as much as possible.